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