SUBROUTINE POT1(PHI,DLON,HT,UN,XI,ETA,DIST) 00052 IMPLICIT REAL*8 (A-H,O-Z) 00053 LOGICAL FIRST 00054 REAL*4 C,C0 00055 DIMENSION P(6) 00056 COMMON/CM/G1(3),G(3,3),CM3,CM2,CM1,C0,C(32760) 00057 DATA PI/3.141592653589793D0/ 00058 COMMON /ENTRY/ FIRST , NMAX 00059**2 C*************************************** 00060 C 00061 C INPUT 00062 C 00063 C PHI LATITUDE (GEODETIC) IN DEGREES 00064 C DLON LONGITUDE (WEST) IN DEGREES 00065 C HT HEIGHT IN METERS 00066 C 00067 C OUTPUT 00068 C 00069 C UN HEIGHT ANOMALY IN METERS 00070 C XI N-S DEFLECTION IN SEC OF ARC 00071 C ETA E-W DEFLECTION IN SECONDS OF ARC (WEST POSITIVE) 00072 C DIST GRAVITY DISTURBANCE IN MGALS 00073 C 00074 C*************************************** 00075 IF(FIRST) GO TO 300 00076 FIRST=.TRUE. 00077 DTR=PI/180.D0 00078 NEGN=-NMAX 00079 CM3=3.986005D+14 00080 CM2=6378137.D0 00081 CM1=0.D0 00082 C0=0.D0 00083 DO 10 I=1,32760 00084 C(I)=0.D0 00085 10 CONTINUE 00086 CALL LOADCS(NMAX) 00087 C0=0.D0 00088 CALL STORC(0,0,0.D0,0.D0) 00089 CALL STORC(1,0,0.D0,0.D0) 00090 CALL STORC(1,1,0.D0,0.D0) 00091 CALL STORC(2,0,0.D0,0.D0) 00092 CALL SETCM(NMAX) 00093 300 CONTINUE 00094 C CONVERT FROM WEST TO EAST LONGITUDE 00095 DLAM=360.D0+DLON 00096**2 C CONVERT FROM GEODETIC TO GEOCENTRIC LATITUDE 00097 DPHI=GCENLT( PHI,298.257D0) 00098 D=CM2+HT 00099 Z=D *DSIN(DPHI*DTR) 00100 R=D *DCOS(DPHI*DTR) 00101 X=R*DCOS(DLAM*DTR) 00102 Y=R*DSIN(DLAM*DTR) 00103 P(1)=DSQRT(X**2+Y**2) 00104 P(2)=DSQRT(X**2+Y**2+Z**2) 00105 COLAT=90.D0-DPHI 00106 P(3)=DCOS(COLAT*DTR) 00107 P(4)=DSIN(COLAT*DTR) 00108 P(5)=DSIN(DLAM*DTR) 00109 P(6)=DCOS(DLAM*DTR) 00110 C PASS NEGATIVE N SINCE POTENTIAL COEFFICIENTS HAVE BEEN 00111 C QUASI-NORMALIZED 00112 POT=POTCC(P,NEGN,1) 00113 C HEIGHT ANOMALY 00114 UN=POT/9.8D0 00115 C N-S 00116 XI = -G1(2)*206264.8D0 / 9.8D0 00117 C OUTPUT E-W DEFLECTIONS IN WEST SYSTEM 00118 C THUS LEAVE OFF NEGATIVE SIGN 00119 C TO CONVERT ETA TO WEST SYSTEM 00120 ETA = G1(1)*206264.8D0 / 9.8D0 00121 C GRAVITY DISTURBANCE 00122 DIST =-G1(3)*1.D+05 00123 RETURN 00124 END 00125 FUNCTION GCENLT(GDETLT,FLAT) 00126 IMPLICIT REAL*8 (A-H,O-Z) 00127 C FUNCTION TO CONVERT FROM GEODETIC LATITUDE TO GEOCENTRIC LATITUDE 00128 C INPUT 00129 C 00130 C GDETLT GEODETIC LATITUDE IN DEGREES 00131 C FLAT 1./FLATTENING EQUAL TO APPROXIMATELY 298.25 00132 C 00133 DATA DTR /.1745329251994330D-1/ 00134 GCENLT=DATAN(DTAN(GDETLT*DTR)*(1.D0-1.D0/FLAT)**2) / DTR 00135 RETURN 00136 END 00137 FUNCTION POTCC(PO,NMAX,ORDER) 00138 C 00139 C 00140 C GI REG.NO. 80039 AUTHOR -C.C.TSCHERNING, NOV 1980 IN ALGOL 00141 C -C.C.GOAD, MAR 1981 TRANSLATED TO FORTRAN 00142 C 00143 C REFERENCES: 00144 C (1) TSCHERNING, C.C.:ON THE CHAIN-RULE METHOD FOR COMPUTING 00145 C POTENTIAL DERIVATIVES. MANUSCRIPTA GEODAETICA, VOL.1, 00146 C PP. 125-141, 1976 00147 C 00148 C (2) GERSTL,M.:VERGLEICH VON ALGORITMEN ZUR SUMMATION VON 00149 C KUGELFLAECHENFUNKTIONEN. VEROEFFENTL. DER BAYER. KOMM. 00150 C F.D. INT. ERDMESSUNG DER BAYER. AKADEMIE DER WISSEN., 00151 C HEFT NR. 38, PP. 81-88, 1978. 00152 C THE PROCEDURE COMPUTES THE VALUE AND UP TO THE SECOND ORDER 00153 C DERIVATIVES OF THE POTENTIAL OF THE EARTH (W) OR OF ITS 00154 C CORRESPONDING AMONALOUS POTENTIAL(T). (THE COMPUTATION OF 00155 C THE SECOND ORDER DERIVATIVES HAS NOT YET BEEN IMPLEMENTED). 00156 C 00157 C THE POTENTIAL IS REPRESENTED BY A SERIES IN SOLID SPHERICAL 00158 C HARMONICS, WITH UN-NORMALIZED OR QUASI-NORMALIZED COEFFICIENTS. 00159 C THE CHAIN-RULE IS USED COMBINED WITH THE CLENSHAW ALGORITHM. 00160 C THE ARRAY C MUST HOLD THE COEFFICIENTS, C(0,0)=1.D0 FOR W AND 00161 C 0.0 FOR T, C(1)=C(1,0),C(2)=C(1,1),C(3)=S(1,1) ETC. UP TO 00162 C C((N+1)**2-1) = S(N,N). 00163 C 00164 C 00165 C PARAMETERS: 00166 C 00167 C (A) INPUT VALUES: 00168 C 00169 C NMAX 00170 C THE ABSOLUTE VALUES OF NMAX IS EQUAL TO THE MAXIMAL DEGREE AND 00171 C ORDER OF THE SERIES. NEGATIVE NMAX INDICATES THAT THE COEFFICIENTS 00172 C ARE QUASI-NORMALIZED. 00173 C 00174 C IORDER 00175 C THE MAXIMAL ORDER OF THE DERIVATIVES (< 2 P.T.). 00176 C 00177 C PO 00178 C ARRAY HOLDING POSITION INFORMATION. PO(6) 00179 C PO(1)=P, THE DISTANCE FROM THE Z (ROTATION) AXIS, 00180 C PO(2)=R, THE DISTANCE FROM THE ORIGIN, 00181 C PO(3),PO(4) COS AND SIN OF THE GEOCENTRIC POLAR ANGLE(COLATITUDE), 00182 C PO(5),PO(6) SIN AND COS OF THE LONGITUDE. 00183 C 00184 C C 00185 C C MUST BE DECLARED WITH BOUNDS (-3: (N+1)**2-1) WHEN THE 00186 C COEFFICIENTS ARE UN-NORMALIZED AND WITH BOUNDS (-3: (N+3)**2-2) 00187 C WHEN THE COEFFICIENTS ARE QUASI-NORMALIZED. C(1) TO C((N+1)**2-1) 00188 C CONTAIN THE COEFFICIENTS AND WE MUST HAVE 00189 C C(-3)=GM 00190 C C(-2)=A THE SEMI-MAJOR AXIS OF THE REF ELLIPSOID 00191 C C(-1)=THE ANGULAR VELOCITY (=0, WHEN DEALING WITH T). 00192 C 00193 C 00194 C SQUARE ROOT ARRAY 00195 C 00196 C C((N+1)**2+K) = SQRT(K), 0.LE.K.LE.2(ABS(N)+1)-1 WHEN N <0 00197 C 00198 C 00199 C MOD--APRIL,1981 00200 C WITH THE USE OF THE TABLE OF SQUARE ROOTS STORED IN ARRAY ROOT, 00201 C THE DIMENSION OF THE C ARRAY ONLY NEEDS TO BE (N+1)**2-1 FOR BOTH 00202 C UN-NORMALIZED AND NORMALIZED COEFFICIENTS. ARRAY ROOT MUST CONTAIN00203 C SQUARE ROOT OF 0 TO N+2 00204 C 00205 C 00206 C (B) RETURN VALUES 00207 C 00208 C G 00209 C THE RESULT IS STORED IN G AS FOLLOWS: 00210 C 00211 C G1(1)=DW/DX, G1(2)=DW/DY, G1(3)=DW/DZ 00212 C G(1,1)=DDW/DDX, G(1,2)=G(2,1)=DDW/DXDY, 00213 C G(1,3)=G(3,1)=DDW/DXDZ, G(2,2)=DDW/DDY, 00214 C G(2,3)=G(3,2)=DDW/DYDZ AND G(3,3)=DDW/DDZ 00215 C WHERE W MAY BE INTERCHANGED WITH T AND 00216 C VARIABLES X, Y, Z ARE THE CARTESIAN COORDINATES 00217 C IN A LOCAL (FIXED) FRAME WITH ORIGIN IN THE POINT 00218 C OF EVALUATION, X POSITIVE NORTH, Y POSITIVE EAST, 00219 C AND Z POSITIVE IN THE DIRECTION OF THE RADIUS 00220 C VECTOR, (CF. REF.(1),EQ (4) AND (5)). 00221 C THE VALUES OF W OR T WILL BE RETURNED IN POTCC. 00222 C 00223 IMPLICIT REAL*8 (A-H,O-Z) 00224 INTEGER CAPN,ORDER,CAPN21 00225 LOGICAL QUASI,DERIV1,DERIV2 00226 REAL*4 C,C0 00227 REAL*8 M21 00228 DIMENSION SML(181),CML(181),SMLP1(182),CMLP1(182),PO(6) 00229 COMMON/SQROOT/DZERO,ROOT(362) 00230 COMMON/CM/G1(3),G(3,3),CM3,CM2,CM1,C0,C(32760) 00231 EQUIVALENCE(SML(1),SMLP1(2)),(CML(1),CMLP1(2)) 00232 DATA IZ/0/ 00233 NAMELIST/DBUG/ CM2,CM1,C0, QUASI,CAPN,S,S2,MAX,CM3 00234 CAPN=NMAX 00235 P=PO(1) 00236 R=PO(2) 00237 T=PO(3) 00238 U=PO(4) 00239 SL=PO(5) 00240 CL=PO(6) 00241 QUASI=.FALSE. 00242 IF(CAPN.LT.0)QUASI=.TRUE. 00243 IF(QUASI)CAPN=-CAPN 00244 S=CM2/R 00245 S2=S**2 00246 CMLP1(1)=1.D0 00247 C CML(0)=1.D0 00248 SMLP1(1)=0.D0 00249 C SML(0)=0.D0 00250 M1=0 00251 DERIV1=.FALSE. 00252 IF(ORDER.GT.0)DERIV1=.TRUE. 00253 DERIV2=.FALSE. 00254 IF( ORDER.GT.1)DERIV2=.TRUE. 00255 C 00256 C SML(M) AND CML(M) ARE THE SINE AND COSINE OF M*LONGITUDE 00257 C 00258 C WRITE(6,DBUG) 00259 DO 1 M=1,CAPN 00260 SML(M)=SML(M1)*CL+CML(M1)*SL 00261 CML(M)=CML(M1)*CL-SML(M1)*SL 00262 1 M1=M 00263 C 00264 CAPN21=CAPN+CAPN+1 00265 VM=0.D0 00266 VXM=0.D0 00267 VYM=0.D0 00268 VZM=0.D0 00269 SQNM1=1.D0 00270 SQNPM1=1.D0 00271 KM=(CAPN+1)**2 00272 C 00273 C WE NOW USE THE CLENSHAW ALGORITHM, CF. REF.(2), EQ(9), 00274 C MODIFIED IN AN OBVIOUS WAY FOLLOWING REF.(1). 00275 C 00276 ITWO=2 00277 DO 7 IM=IZ,CAPN 00278 M=CAPN-IM 00279 IF(M.EQ.0)ITWO=1 00280 KM=KM-ITWO 00281 K=KM 00282 N21=CAPN21 00283 VY=0.D0 00284 VZ1=0.D0 00285 VZ=0.D0 00286 VY1=0.D0 00287 VX=0.D0 00288 VX1=0.D0 00289 V=0.D0 00290 V1=0.D0 00291 CM=CML(M) 00292 SM=SML(M) 00293 NM1=CAPN-M+2 00294 N1=CAPN+1 00295 NPM1=CAPN+M+2 00296 N=CAPN+1 00297 DO 5 IN=M,CAPN 00298 N=N-1 00299 NM2=NM1 00300 NM1=NM1-1 00301 NPM1=NPM1-1 00302 IF(.NOT.QUASI) GO TO 2 00303 SQNM2=SQNM1 00304 C SQNM1=C(NM1) 00305 SQNM1=ROOT(NM1) 00306 SQNPM2=SQNPM1 00307 SQNPM1=ROOT(NPM1) 00308 C SQNPM1=C(NPM1) 00309 SQ1=SQNM1*SQNPM1 00310 A1=S*N21/SQ1 00311 B2=-S2*SQ1/(SQNM2*SQNPM2) 00312 GO TO 3 00313 2 A1=S*N21/NM1 00314 B2=-S2*NPM1/NM2 00315 3 A1T=A1*T 00316 N21=N21-2 00317 CK=C(K) 00318 CK1=C(K+1) 00319 SMALLC=CM*CK+SM*CK1 00320 K=K-N21 00321 V2=V1 00322 V1=V 00323 V=V1*A1T+V2*B2+SMALLC 00324 IF(.NOT.DERIV1) GO TO 4 00325 D=-SM*CK+CM*CK1 00326 VX2=VX1 00327 VX1=VX 00328 VX=VX1*A1T+V1*A1*U+VX2*B2 00329 VY2=VY1 00330 VY1=VY 00331 VY=VY1*A1T+VY2*B2+D 00332 VZ2=VZ1 00333 VZ1=VZ 00334 VZ=VZ1*A1T+VZ2*B2-N1*SMALLC 00335 N1=N 00336 4 CONTINUE 00337 5 CONTINUE 00338 AUX=NPM1 00339 IF(QUASI)AUX=SQNPM1/SQNPM2 00340 M21=S*AUX 00341 IF(.NOT.DERIV1) GO TO 6 00342 VXM=VX+M21*(-T*VM+U*VXM) 00343 AUX=U 00344 IF(M.EQ.0)AUX=1.D0 00345 VYM=M*VY+M21*AUX*VYM 00346 VZM=VZ+M21*U*VZM 00347 C 00348 6 VM=V+M21*U*VM 00349 7 CONTINUE 00350 C 00351 C NOW ADD THE CONTRIBUTIONS FROM THE ROTATIONAL POTENTIAL 00352 C 00353 OM2=CM1**2 00354 S=CM3/R 00355 POTCC=S*VM+OM2*P**2/2 00356 S=S/R 00357 G1(1)=S*VXM-T*P*OM2 00358 G1(2)=S*VYM 00359 G1(3)=VZM*S+U**2*OM2*R 00360 RETURN 00361 END 00362 SUBROUTINE SETCM(CAPN) 00363 IMPLICIT REAL*8(A-H,O-Z) 00364 INTEGER CAPN 00365 REAL*4 C,C0 00366 COMMON/SQROOT/DZERO,ROOT(362) 00367 COMMON/CM/G1(3),G(3,3),CM3,CM2,CM1,C0,C(32760) 00368 C 00369 C THIS ROUTINE SETS THE SQUARE ROOT TABLE IN COMMON 00370 C CM AND CREATES QUASI-NORMALIZED COEFFICIENTS 00371 C FROM NORMALIZED COEFFICIENTS 00372 C 00373 C 00374 C APRIL 1981 CCG 00375 C THIS VERSION STORES SQUARE ROOT TABLE IN COMMON SQROOT 00376 C 00377 DZERO=0.D0 00378 DO 22 I=1,362 00379 22 ROOT(I)=DSQRT(DFLOAT(I)) 00380 G1(1)=0.D0 00381 G1(2)=0.D0 00382 G1(3)=0.D0 00383 SMALLC=1.D0 00384 IF(C0.NE.0.D0)SMALLC=1.D0/C0 00385 SQ2=DSQRT(2.D0) 00386 DO 200 N=1,CAPN 00387 N2=N+N 00388 S21=DSQRT(N2+1.D0) 00389 K=N**2 00390 C D IS THE QUASI-NORMALIZATION FACTOR FOR ZONAL TERMS 00391 D=SMALLC*S21 00392 C(K)=C(K)*D 00393 C GG IS THE QUASI-NORMALIZATION FACTOR FOR NON-ZONAL TERMS 00394 GG=D*SQ2 00395 DO 100 J=1,N 00396 KJ2=J+J+K 00397 C(KJ2-1)=C(KJ2-1)*GG 00398 C(KJ2)=C(KJ2)*GG 00399 100 CONTINUE 00400 200 CONTINUE 00401 RETURN 00402 END 00403 SUBROUTINE STORC(N,M,CNM,SNM) 00404 C 00405 C STORE INDIVIDUAL C AND S TERMS 00406 C 00407 IMPLICIT REAL*8 (A-H,O-Z) 00408 REAL*4 C,C0 00409 COMMON/CM/G1(3),G(3,3),CM3,CM2,CM1,C0,C(32760) 00410 C 00411 C SUM OF THE PREVIOUS NUMBER OF TERMS 00412 C 00413 J=(N-1)*(N+1) 00414 IF(M.EQ.0) GO TO 10 00415 K=M+M 00416 C(J+K)=CNM 00417 C(J+K+1)=SNM 00418 RETURN 00419 10 C(J+1)=CNM 00420 RETURN 00421 END 00422 SUBROUTINE LOADCS(NMAX) 00423 IMPLICIT REAL*8 (A-H,O-Z) 00424 1 READ(12, END=2)N,M,C,S IF(N.GT.NMAX) GO TO 1 00426 IF(N.LT.5)WRITE(6,55)N,M,C,S 00427 55 FORMAT(1X,6G20.12) 00428 1000 FORMAT(2I3,2G15.8) 00429**2 CALL STORC(N,M,C,S) 00430 GO TO 1 00431 2 CONTINUE 00432 RETURN 00433 END 00434