/usr/share/gnudatalanguage/astrolib/geodetic2geo.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 | ;+
; NAME:
; GEODETIC2GEO
;
; PURPOSE:
; Convert from geodetic (or planetodetic) to geographic coordinates
; EXPLANATION:
; Converts from geodetic (latitude, longitude, altitude) to geographic
; (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 (planetodetic to planetographic 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:
; gcoord = geodetic2geo(ecoord, [ PLANET= ] )
;
; INPUT:
; ecoord = a 3-element array of geodetic [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. Numeric scalar
; POLAR_RADIUS : Self-explanatory. In km. If not set, PLANET's value is
; used. Numeric scalar
;
; OUTPUT:
; a 3-element array of geographic [latitude,longitude,altitude], or an
; array [3,n] of n such coordinates, double precision
;
; The geographic and geodetic longitudes will be identical.
; COMMON BLOCKS:
; None
;
; EXAMPLES:
;
; IDL> geod=[90,0,0] ; North pole, altitude 0., in geodetic coordinates
; IDL> geo=geodetic2geo(geod)
; IDL> PRINT,geo
; 90.000000 0.0000000 -21.385000
;
; As above, but the equivalent planetographic coordinates for Mars
; IDL> geod=geodetic2geo(geod,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 geodetic2geo,ecoord,PLANET=planet, $
EQUATORIAL_RADIUS=equatorial_radius, POLAR_RADIUS=polar_radius
sz_ecoord = size(ecoord,/DIMEN)
if sz_ecoord[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]
;f=1/298.257D ; flattening = (Re-Rp)/Re
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
elat = DOUBLE(ecoord[0,*])*!DPI/180.
elon = DOUBLE(ecoord[1,*])
ealt = DOUBLE(ecoord[2,*])
beta=sqrt(1-(e*sin(elat))^2)
r=(Re/beta + ealt)*cos(elat)
z=(Re*(1-e^2)/beta + ealt)*sin(elat)
glat=atan(z,r)*180./!DPI
glon=elon
galt=sqrt(r^2+z^2) - Re
RETURN,[glat,glon,galt]
END
;===================================================================================
|