pro geod2geoc,hd,latd,hc,latc,flat=flat,req=req
; convert geodetic altitude and latitude (hd, latd) [km,deg]
; to geocentric altitude and latitude (hc,latc). 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

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

phi=!dtor*latd
N= re/sqrt((one-e2*sin(phi)^2))
E=(N+hd)*cos(phi)
F=(N*(1-e2)+hd)*sin(phi)

latc=!radeg*atan(E,F)
latc=!radeg*atan(F,E)
hc=sqrt((E^2+F^2))-re
return
end