SUBROUTINE FIELDG (DLAT,DLONG,ALT,TM,NMX,L,X,Y,Z,F) EQUIVALENCE (SHMIT(1,1),TG(1,1)) COMMON/FLDCOM/ST,CT,SPH,CPH,R,BT,BP,BR,B,NMAX COMMON /COEFFS/TG(18,18) DIMENSION G(18,18),GT(18,18),SHMIT(18,18),AID(15) IF (A.EQ.6378.165) IF(L)19,1,2 DATA MAXN,TZERO,J,K/10,1960.,0,0/,G/ X0.0,-30424.9,-1536.1,1300.9,957.6,-227.7,49.8,70.9,4.8,9.9,8*0.,57 X74.8,-2161.6,3000.2,-1987.0,802.8,359.5,60.7,-57.2,6.7,2.9,8*0.,-1 X949.8,204.3,1585.3,1290.4,502.6,231.3,4.5,5.6,-8.8,7.4,8*0.,-431.0 X,230.8,-130.0,871.2,-394.0,-31.2,-241.7,7.5,-13.8,-15.6,8*0.,152.0 X,-268.4,2.9,-250.5,271.4,-157.3,-1.2,-24.4,-3.3,11.4,8*0.,8.6,121. X2,-116.0,-110.4,79.9,-65.2,0.5,-1.5,7.1,11.1,8*0.,-11.9,102.8,60.9 X,-27.2,-12.4,-11.6,-109.1,14.1,-5.6,1.,8*0.,-54.,-24.4,-9.1,2.2,27 X.6,-21.1,-20.1,5.8,11.7, 9*0.,6.9,-12.2,5.8,-17.0,2.6,23.6,-2.5,-1 X6.0,6.4,1.6,8*0.,-22.0,15.6,5.1,-3.5,-1.8,9.6,12.1,0.2,-2.5,1.5,15 X2*0./ DATA GT/ X0.0,20.59,-29.07,2.66,-0.86,2.55,-0.70,11*0.,-3.94,6.02,1.21,-10.0 X3,1.94,-0.08,0.99,11*0.,-13.69,-15.78,-0.70,1.63,-1.17,1.53,0.85,1 X1*0.,6.49,2.93,-9.24,-1.30,-0.54,-0.43,2.11,11*0.,1.77,-1.54,3.18, X-5.48,-4.17,-0.72,1.57,11*0.,3.04,2.88,-1.86,1.25,0.80,1.64,-0.09, X11*0.,-1.39,0.12,1.53,-0.73,-0.06,0.45,0.06,209*0. / A=6378.165 FLAT=1.-1./298.3 A2=A**2 A4=A**4 B2=(A*FLAT)**2 A2B2=A2*(1.-FLAT**2) A4B4=A4*(1.-FLAT**4) IF (L)14,14,2 1 IF (TM-TLAST) 17,19,17 2 READ (3,3) J,K,TZERO,AID ! WAS UNIT 5 3 FORMAT (2I1,1X,F6.1,15A4,A3) L=0 WRITE (6,4) J,K,TZERO,AID 4 FORMAT (2I3,5X6HEPOCH=,F7.1,5X15A4,A3) MAXN=0 TEMP=0. 5 READ (3,6) N,M,GNM,HNM,GTNM,HTNM ! WAS UNIT 5 6 FORMAT(2I3,6F11.4) IF (N.LE.0) GOTO7 MAXN=(MAX0(N,MAXN)) G(N,M)=GNM GT(N,M)=GTNM TEMP= MAX1(TEMP, ABS(GTNM)) IF (M.EQ.1) GOTO5 G(M-1,N)=HNM GT(M-1,N)=HTNM GO TO 5 7 WRITE (6,8) 8 FORMAT (6H0 N M,6X1HG,10X1HH,9X2HGT,9X2HHT//) DO 12 N=2,MAXN DO 12 M=1,N MI=M-1 IF (M.EQ.1) GOTO10 WRITE (6,9) N,M,G(N,M),G(MI,N),GT(N,M),GT(MI,N) 9 FORMAT (2I3,4F11.4) GO TO 12 10 WRITE (6,11) N,M,G(N,M),GT(N,M) 11 FORMAT (2I3,F11.4,11X,F11.4) 12 CONTINUE WRITE (6,13) 13 FORMAT (1H1) IF (TEMP.EQ.0.) L=-1 14 IF (K.NE.0) GOTO17 SHMIT(1,1)=-1. DO 15 N=2,MAXN SHMIT(N,1)=SHMIT(N-1,1)*FLOAT(2*N-3)/FLOAT(N-1) SHMIT(1,N)=0. JJ=2 DO 15 M=2,N SHMIT(N,M)=SHMIT(N,M-1)*SQRT(FLOAT((N-M+1)*JJ)/FLOAT(N+M-2)) SHMIT(M-1,N)=SHMIT(N,M) 15 JJ=1 DO 16 N=2,MAXN DO 16 M=1,N G(N,M)=G(N,M)*SHMIT(N,M) GT(N,M)=GT(N,M)*SHMIT(N,M) IF (M.EQ.1) GOTO16 G(M-1,N)=G(M-1,N)*SHMIT(M-1,N) GT(M-1,N)=GT(M-1,N)*SHMIT(M-1,N) 16 CONTINUE 17 T=TM-TZERO DO 18 N=1,MAXN DO 18 M=1,N TG(N,M)=G(N,M)+T*GT(N,M) IF (M.EQ.1) GOTO18 TG(M-1,N)=G(M-1,N)+T*GT(M-1,N) 18 CONTINUE TLAST=TM 19 SINLA= SIN(DLAT/57.29578) RLONG=DLONG/57.29578 CPH= COS(RLONG) SPH= SIN(RLONG) IF (J.EQ.0) GOTO20 R=ALT+6371.2 CT=SINLA GO TO 21 20 SINLA2=SINLA**2 COSLA2=1.-SINLA2 DEN2=A2-A2B2*SINLA2 DEN= SQRT(DEN2) FAC=(((ALT*DEN)+A2)/((ALT*DEN)+B2))**2 CT=SINLA/ SQRT(FAC*COSLA2+SINLA2) R= SQRT(ALT*(ALT+2.*DEN)+(A4-A4B4*SINLA2)/DEN2) 21 ST= SQRT(1.-CT**2) NMAX=MIN0(NMX,MAXN) CALL FIELD Y=BP F=B IF (J) 22,23,22 22 X=-BT Z=-BR RETURN C TRANSFORMS FIELD TO GEODETIC DIRECTIONS 23 SIND=SINLA*ST- SQRT(COSLA2)*CT COSD= SQRT(1.0-SIND**2) X=-BT*COSD-BR*SIND Z=BT*SIND-BR*COSD RETURN END