pro geoc2geod,hc,latc,hd,latd,flat=flat,req=req
; convert geocentric altitude and latitude (hc, latc) [km,deg]
; to geodetic altitude and latitude (hd,latd). Uses
; algorithm from Hedman, J. Spacecraft, vol 7, No. 8, aug 1970
; flat is the optional flattening term for the shape of the earth
;	by default, flat =1./298.25722356
; req is the mean equitorial radius
;	by default, req =6378.1370
; itmax is the maximum iterations allowed (default=10)
; tol is the convergence criteria for latitude in deg., default 1.e-5

common earth,re,e2,E,F

f=(1d)/298.25722356
re=6378.1370
if keyword_set(flat) then f=flat
if keyword_set(req) then re=req
one=1d
e2=one-(one-f)^2
b=re*sqrt(one-e2)
E=(hc+re)*cos(!dtor*latc)
F=(hc+re)*sin(!dtor*latc)


phi0=atan(re^2*F,b^2*E)
phi=newton(phi0,'minfun')
latd=!radeg*phi
nsub=re/sqrt(one-e2*sin(phi)^2)
xsub=nsub*cos(phi)
ysub=(nsub*(one-e2))*sin(phi)
biga=sqrt((E-xsub)^2+(F-ysub)^2)
bigb=xsub*(E-xsub)+ysub*(F-ysub)
hd=biga*(bigb)/abs(bigb)
return
end

function minfun,phi
common earth,re,e2,E,F
one=1d
return,re*e2/sqrt(one-e2*sin(phi)^2)-E/cos(phi)+F/sin(phi)
end