/usr/share/gnudatalanguage/astrolib/geo2geodetic.pro is in gdl-astrolib 2018.02.16+dfsg-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 | ;+
; NAME:
; GEO2GEODETIC
;
; PURPOSE:
; Convert from geographic/planetographic to geodetic coordinates
; EXPLANATION:
; Converts from geographic (latitude, longitude, altitude) to geodetic
; (latitude, longitude, altitude). In geographic coordinates, the
; Earth is assumed a perfect sphere with a radius equal to its equatorial
; radius. The geodetic (or ellipsoidal) coordinate system takes into
; account the Earth's oblateness.
;
; Geographic and geodetic longitudes are identical.
; Geodetic latitude is the angle between local zenith and the equatorial plane.
; Geographic and geodetic altitudes are both the closest distance between
; the satellite and the ground.
;
; The PLANET keyword allows a similar transformation for the other
; planets (planetographic to planetodetic coordinates).
;
; The EQUATORIAL_RADIUS and POLAR_RADIUS keywords allow the
; transformation for any ellipsoid.
;
; Latitudes and longitudes are expressed in degrees, altitudes in km.
;
; REF: Stephen P. Keeler and Yves Nievergelt, "Computing geodetic
; coordinates", SIAM Rev. Vol. 40, No. 2, pp. 300-309, June 1998
;
; Planetary constants from "Allen's Astrophysical Quantities",
; Fourth Ed., (2000)
;
; CALLING SEQUENCE:
; ecoord=geo2geodetic(gcoord,[ PLANET=,EQUATORIAL_RADIUS=, POLAR_RADIUS=])
;
; INPUT:
; gcoord = a 3-element array of geographic [latitude,longitude,altitude],
; or an array [3,n] of n such coordinates.
;
;
; OPTIONAL KEYWORD INPUT:
; PLANET = keyword specifying planet (default is Earth). The planet
; may be specified either as an integer (1-9) or as one of the
; (case-independent) strings 'mercury','venus','earth','mars',
; 'jupiter','saturn','uranus','neptune', or 'pluto'
;
; EQUATORIAL_RADIUS : Self-explanatory. In km. If not set, PLANET's
; value is used.
; POLAR_RADIUS : Self-explanatory. In km. If not set, PLANET's value is
; used.
;
; OUTPUT:
; a 3-element array of geodetic/planetodetic [latitude,longitude,altitude],
; or an array [3,n] of n such coordinates, double precision.
;
; COMMON BLOCKS:
; None
;
; RESTRICTIONS:
;
; Whereas the conversion from geodetic to geographic coordinates is given
; by an exact, analytical formula, the conversion from geographic to
; geodetic isn't. Approximative iterations (as used here) exist, but tend
; to become less good with increasing eccentricity and altitude.
; The formula used in this routine should give correct results within
; six digits for all spatial locations, for an ellipsoid (planet) with
; an eccentricity similar to or less than Earth's.
; More accurate results can be obtained via calculus, needing a
; non-determined amount of iterations.
; In any case,
; IDL> PRINT,geodetic2geo(geo2geodetic(gcoord)) - gcoord
; is a pretty good way to evaluate the accuracy of geo2geodetic.pro.
;
; EXAMPLES:
;
; Locate the geographic North pole, altitude 0., in geodetic coordinates
; IDL> geo=[90.d0,0.d0,0.d0]
; IDL> geod=geo2geodetic(geo); convert to equivalent geodetic coordinates
; IDL> PRINT,geod
; 90.000000 0.0000000 21.385000
;
; As above, but for the case of Mars
; IDL> geod=geo2geodetic(geo,PLANET='Mars')
; IDL> PRINT,geod
; 90.000000 0.0000000 18.235500
;
; MODIFICATION HISTORY:
; Written by Pascal Saint-Hilaire (shilaire@astro.phys.ethz.ch), May 2002
; Generalized for all solar system planets by Robert L. Marcialis
; (umpire@lpl.arizona.edu), May 2002
; Modified 2002/05/18, PSH: added keywords EQUATORIAL_RADIUS and
; POLAR_RADIUS
;-
;================================================================================
FUNCTION geo2geodetic,gcoord,PLANET=planet, $
EQUATORIAL_RADIUS=equatorial_radius, POLAR_RADIUS=polar_radius
sz_gcoord = size(gcoord,/DIMEN)
if sz_gcoord[0] LT 3 then message, $
'ERROR - 3 coordinates (latitude,longitude,altitude) must be specified'
if N_elements(PLANET) GT 0 then begin
if size(planet,/tname) EQ 'STRING' then begin
choose_planet=['mercury','venus','earth','mars','jupiter','saturn', $
'uranus','neptune','pluto']
index=where(choose_planet eq strlowcase(planet))
index=index[0] ; make it a scalar
if index eq -1 then index = 2 ; default is Earth
endif else index = planet-1
endif else index=2
Requator = [2439.7d0,6051.8d0,6378.137D, 3397.62d0, 71492d0, $
60268.d0, 25559.d0, 24764.d0, 1195.d0]
Rpole = [2439.7d0, 6051.8d0, 6356.752d0, 3379.3845d0, 67136.5562d0, $
54890.7686d0, 24986.1354d0, 24347.6551d0, 1195.d0]
Re = Requator[index] ; equatorial radius
Rp = Rpole[index] ; polar radius
IF KEYWORD_SET(EQUATORIAL_RADIUS) THEN Re=DOUBLE(equatorial_radius[0])
IF KEYWORD_SET(POLAR_RADIUS) THEN Rp=DOUBLE(polar_radius[0])
e = sqrt(Re^2 - Rp^2)/Re
;f=1/298.257D ; flattening = (Re-Rp)/Re [not needed, here]
glat=DOUBLE(gcoord[0,*])*!DPI/180.
glon=DOUBLE(gcoord[1,*])
galt=DOUBLE(gcoord[2,*])
x= (Re+galt) * cos(glat) * cos(glon)
y= (Re+galt) * cos(glat) * sin(glon)
z= (Re+galt) * sin(glat)
r=sqrt(x^2+y^2)
s=(r^2 + z ^2)^0.5 * (1 - Re*((1-e^2)/((1-e^2)*r^2 + z^2))^0.5)
t0=1+s*(1- (e*z)^2/(r^2 + z^2) )^0.5 /Re
dzeta1=z * t0
xi1=r*(t0 - e^2)
rho1= (xi1^2 + dzeta1^2)^0.5
c1=xi1/rho1
s1=dzeta1/rho1
b1=Re/(1- (e*s1)^2)^0.5
u1= b1*c1
w1= b1*s1*(1- e^2)
ealt= ((r - u1)^2 + (z - w1)^2)^0.5
elat= atan(s1,c1)
elat=elat*180./!DPI
elon=glon
RETURN,[elat,elon,ealt]
END
;===============================================================================
|