C C C THREE DIMENSIONAL LEAST SQUARES FIT, USING SPHERICAL HARMONICS. C C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(20,20),B(20,20),AO(20) DIMENSION ALG(10,10,36),ALO(10,36) DIMENSION VCM(20,20),VCO(20) DIMENSION PY(36) 1 FORMAT(I3) 50 FORMAT(5X,'A(',I2,',0)',2X,F13.8,36X,F15.10,/) 51 FORMAT(/,5X,'A(',I2,',',I2,')',2X,F13.8,5X,'B(',I2,',',I2,')',2X, *F13.8,10X,F15.10) 52 FORMAT(10X,'THREE DIMENSIONAL LEAST SQUARES FIT, USING SPHERICAL *HARMONICS',////) 53 FORMAT('1') 54 FORMAT(//,14X,'COSINE',21X,'SINE',19X,'VARIANCE',/) 55 FORMAT(12X,'COEFFICIENT',15X,'COEFFICIENT',////) C C FUNCTION TO BE APPROXIMATED C FA(X,Y)=(DSIN(Y)*DCOS(Y))**2*DSIN(X)*DCOS(X) C C NOTE. 5 DEGREE INTERVAL USED FOR NUMERICAL INTEGRATION C RHO=206264.80625D0/3600.D0 PI=3.14159265358979D0 READ(5,1) L AK=L**2/2.D0+L/2.D0+1.D0 WRITE(6,53) WRITE(6,52) WRITE(6,54) WRITE(6,55) C C DETERMINATION OF COEFFICIENTS A(N,M), B(N,M) C VR=0.D0 DO 10 NN=1,L DO 20 MM=1,NN A(NN,MM)=0.D0 B(NN,MM)=0.D0 T=0.D0 DO 80 IA=1,19 CALL ALEG(NN,MM,T,POLY) PY(IA)=POLY 80 T=T+5.D0/RHO MB=MM+NN MC=MB-MB/2*2+1 IF(MC.EQ.2) GOTO82 DO 81 K=1,17 81 PY(19+K)=PY(19-K) GOTO84 82 DO 83 K=1,17 83 PY(19+K)=-PY(19-K) 84 CONTINUE T=0.D0 DO 21 IA=1,36 AL=0.D0 ALG(NN,MM,IA)=PY(IA) DO 22 IB=1,72 A(NN,MM)=A(NN,MM)+FA(AL,T)*PY(IA)*DCOS(MM*AL)*DSIN(T) B(NN,MM)=B(NN,MM)+FA(AL,T)*PY(IA)*DSIN(MM*AL)*DSIN(T) AL=AL+5.D0/RHO 22 CONTINUE T=T+5.D0/RHO 21 CONTINUE VCM(NN,MM)=(2.D0*NN+1.D0)/(2.D0*PI)*FAC(NN-MM)/FAC(NN+MM) A(NN,MM)=VCM(NN,MM)*A(NN,MM)*(5.D0/RHO)**2 B(NN,MM)=VCM(NN,MM)*B(NN,MM)*(5.D0/RHO)**2 20 CONTINUE 10 CONTINUE C C DETERMINATION OF COEFFICIENTS A(N,0) C LL=L+1 DO 23 N=1,LL NL=N-1 AO(N)=0.D0 T=0.D0 DO 62 IA=1,19 CALL ALEG(NL,0,T,POLY) PY(IA)=POLY T=T+5.D0 62 CONTINUE MD=N-N/2*2+1 IF(MD.EQ.2) GOTO64 DO 63 K=1,17 63 PY(19+K)=PY(19-K) GOTO65 64 DO 66 K=1,17 66 PY(19+K)=-PY(19-K) 65 CONTINUE T=0.D0 DO 24 IA=1,36 AL=0.D0 ALO(N,IA)=PY(IA) DO 25 IB=1,72 AO(N)=AO(N)+FA(AL,T)*PY(IA)*DSIN(T) AL=AL+5.D0/RHO 25 CONTINUE T=T+5.D0/RHO 24 CONTINUE VCO(N)=(2.D0*NL+1.D0)/(4.D0*PI) AO(N)=AO(N)*VCO(N)*(5.D0/RHO)**2 23 CONTINUE C C CALCULATION OF POLYNOMIAL AND RESIDUALS C RES=0.D0 T=0.D0 DO 28 IA=1,36 AL=0.D0 DO 29 IB=1,72 SH=0.D0 DO 11 NN=1,L DO 26 MM=1,NN POLY=ALG(NN,MM,IA) SH=SH+A(NN,MM)*POLY*DCOS(MM*AL)+B(NN,MM)*POLY*DSIN(MM*AL) 26 CONTINUE 11 CONTINUE LL=L+1 DO 27 N=1,LL NL=N-1 POLY=ALO(N,IA) SH=SH+AO(N)*POLY 27 CONTINUE RES=RES+(FA(AL,T)-SH)**2 AL=AL+5.D0/RHO 29 CONTINUE T=T+5.D0/RHO 28 CONTINUE C C CALCULATION OF VARIANCE AND OUTPUT OF RESULTS C VAR=RES/(2592-AK) LL=L+1 DO 31 N=1,LL NL=N-1 VCO(N)=VCO(N)*VAR WRITE(6,50) NL,AO(N),VCO(N) 31 CONTINUE DO 36 NN=1,L DO 36 MM=1,NN VCM(NN,MM)=VCM(NN,MM)*VAR 36 WRITE(6,51) NN,MM,A(NN,MM),NN,MM,B(NN,MM),VCM(NN,MM) WRITE(6,53) STOP END C C ROUTINE TO CALCULATE LEGENDRE POLYNOMIALS C SUBROUTINE ALEG(NN,MM,T,POLY) IMPLICIT REAL*8 (A-H,O-Z) POLY=0.D0 IR=(NN-MM)/2 IRR=IR+1 PY=0.D0 DO 13 KK=1,IRR K=KK-1 PY=PY+(0.D0-1.D0)**K*FAC(2*NN-2*K)/(FAC(K)*FAC(NN-K)*FAC(NN-MM-2*K *))*(DCOS(T))**(NN-MM-2*K) K=K+1 13 CONTINUE IF(MM.EQ.0) GOTO14 POLY=(1.D0-(DCOS(T))**2)**(MM/2.D0)/(2.D0**NN)*PY RETURN 14 POLY=PY/(2.D0**NN) RETURN END C C ROUTINE TO CALCULATE FACTORIALS C FUNCTION FAC(N) IMPLICIT REAL*8 (A-H,O-Z) FAC=1.D0 M=N NA=N IF(M.EQ.0) RETURN DO 10 I=1,M FAC=FAC*NA NA=NA-1 10 CONTINUE RETURN END