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