convrt.f
80 lines
| 2.6 KiB
| text/x-fortran
|
FortranFixedLexer
| r0 | C $Id: convrt.f 3304 2011-01-17 15:25:59Z brideout $ | |||
| C | ||||
| SUBROUTINE CONVRT(I,GDLAT,GDALT,GCLAT,RKM) | ||||
| C | ||||
| C jmh - 11/79 ans fortran 66 | ||||
| C | ||||
| C Converts between geodetic and geocentric coordinates. the | ||||
| C reference geoid is that adopted by the iau in 1964. a=6378.16, | ||||
| C b=6356.7746, f=1/298.25. the equations for conversion from | ||||
| C geocentric to geodetic are from astron. j., vol 66, 1961, p. 15. | ||||
| C | ||||
| C Input: | ||||
| C I - 1 geodetic to geocentric, 2 geocentric to geodetic | ||||
| C Input, Output: | ||||
| C GDLAT - geodetic latitude (degrees) | ||||
| C GDALT - altitude above geoid (km) | ||||
| C GCLAT - geocentric latitude (degrees) | ||||
| C RKM - geocentric radial distance (km) | ||||
| C | ||||
| C .. Scalar Arguments .. | ||||
| DOUBLE PRECISION GCLAT,GDALT,GDLAT,RKM | ||||
| INTEGER I | ||||
| C .. | ||||
| C .. Local Scalars .. | ||||
| DOUBLE PRECISION A,A2,A4,A6,A8,AB2,C2CL,C4CL,CCL,CL2,COSBET, | ||||
| * COSLAT,DLTCL,DTR,EP2,GCL,GDL,RER,RGEOID,S2CL, | ||||
| * S4CL,S6CL,S8CL,SB2,SCL,SINBET,SINLAT,SL2,X,Y | ||||
| C .. | ||||
| C .. Intrinsic Functions .. | ||||
| INTRINSIC DATAN2,DCOS,DMIN1,DSIN,DSQRT | ||||
| C .. | ||||
| C .. Data statements .. | ||||
| C | ||||
| DATA A/6378.16D0/,AB2/1.0067397D0/,EP2/.0067397D0/ | ||||
| DATA DTR/.0174532925199D0/ | ||||
| C .. | ||||
| C | ||||
| IF (I.EQ.2) THEN | ||||
| C | ||||
| C .....geocentic to geodetic..... | ||||
| RER = RKM/A | ||||
| A2 = ((-1.4127348D-8/RER+.94339131D-8)/RER+.33523288D-2)/RER | ||||
| A4 = (((-1.2545063D-10/RER+.11760996D-9)/RER+.11238084D-4)/RER- | ||||
| * .2814244D-5)/RER | ||||
| A6 = ((54.939685D-9/RER-28.301730D-9)/RER+3.5435979D-9)/RER | ||||
| A8 = (((320.D0/RER-252.D0)/RER+64.D0)/RER-5.D0)/RER* | ||||
| * .98008304D-12 | ||||
| GCL = DTR*GCLAT | ||||
| CCL = DCOS(GCL) | ||||
| SCL = DSIN(GCL) | ||||
| S2CL = 2.D0*SCL*CCL | ||||
| C2CL = 2.D0*CCL*CCL - 1.D0 | ||||
| S4CL = 2.D0*S2CL*C2CL | ||||
| C4CL = 2.D0*C2CL*C2CL - 1.D0 | ||||
| S8CL = 2.D0*S4CL*C4CL | ||||
| S6CL = S2CL*C4CL + C2CL*S4CL | ||||
| DLTCL = S2CL*A2 + S4CL*A4 + S6CL*A6 + S8CL*A8 | ||||
| GDLAT = GCLAT + DLTCL/DTR | ||||
| GDALT = RKM - A/DSQRT(1.D0+EP2*SCL*SCL) | ||||
| ELSE | ||||
| C | ||||
| C .....geodetic to geocentric..... | ||||
| GDL = DTR*GDLAT | ||||
| SINLAT = DSIN(GDL) | ||||
| COSLAT = DCOS(GDL) | ||||
| SL2 = SINLAT*SINLAT | ||||
| CL2 = AB2*COSLAT | ||||
| CL2 = CL2*CL2 | ||||
| SINBET = SINLAT/DSQRT(SL2+CL2) | ||||
| SB2 = DMIN1(SINBET*SINBET,1.D0) | ||||
| COSBET = DSQRT(1.D0-SB2) | ||||
| RGEOID = A/DSQRT(1.D0+EP2*SB2) | ||||
| X = RGEOID*COSBET + GDALT*COSLAT | ||||
| Y = RGEOID*SINBET + GDALT*SINLAT | ||||
| RKM = DSQRT(X*X+Y*Y) | ||||
| GCLAT = DATAN2(Y,X)/DTR | ||||
| END IF | ||||
| RETURN | ||||
| C | ||||
| END | ||||
