function eci_to_geodetic, jd, eci, spherical=spherical, debug=debug
;
; Given Julian Day Number and an ECI position, returns the
; geodetic position.
;
; B. Knapp, 2004-02-10
;
  rad2deg = 180.d0/!dpi
  if not keyword_set(spherical) then begin
;
;   ECI coordinates are rectangular
    lon = sidtim(jd)-atan(eci[1,*],eci[0,*])*rad2deg
    r = sqrt(eci[0,*]^2+eci[1,*]^2)
    z = eci[2,*]
  endif else begin
;
;   ECI coordinates are spherical (degrees)
    lon = sidtim(jd)-eci[0,*]
    phi = eci[1,*]/rad2deg
    r = eci[2,*]*cos(phi)
    z = eci[2,*]*sin(phi)
  endelse
;
  p = where(lon lt -180.d0, np)
  if np gt 0 then lon[p] = lon[p]+360.d0
  q = where(lon gt +180.d0, nq)
  if nq gt 0 then lon[q] = lon[q]-360.d0
  geodetic, r, z, lat, h
  if keyword_set(debug) then print, r, z, lat*rad2deg, h, $
    format="('r, z, lat, h:',4f20.14)"
;
; Return the correct type
  info = size(eci)
  if info[0] eq 0 then $
    return,[lon,lat*rad2deg,h] $
  else $
    return,transpose([[lon],[lat*rad2deg],[h]])
;
end