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