;+ ; 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 ;===============================================================================