pro geodetic, r, z, phi, h
;
; Given Cartesian coordinates r (= sqrt(x^2+y^2)), z in kilometers, in a
; meridional plane of the Earth, this procedure returns the corresponding
; geodetic coordinates phi (latitude, radians) and h (geodetic height, km).
;
; B. Knapp, 2001-01-11
;
; Print usage?
  if n_params() lt 2 or n_elements(r) ne n_elements(z) then begin
     print, ' '
     print, ' Given Cartesian coordinates r (= sqrt(x^2+y^2)), and z'
     print, ' (both in km) in a meridional plane of the Earth, this'
     print, ' procedure returns the corresponding geodetic coordinates'
     print, ' phi (latitude, radians) and h (geodetic height, km).'
     print, ' '
     print, '     geodetic, r, z, phi, h'
     print, ' '
     return
  endif
;
; Use 32-bit or 64-bit shareable object library?
  mbits = ((!version.release ge 5.4) ? !version.memory_bits : 32)
  dll = getenv('bgkroot')+string(mbits,"('/lib/astron',i2,'.so')")
;
; Define outputs and then call the Fortran subroutine
  n = n_elements(r)
  r_loc = double(r)
  z_loc = double(z)
  if n eq 1 then begin
     phi = 0.d0
     h = 0.d0
     result = call_external(dll, 'geodetic_arg_', $
        r_loc[0], z_loc[0], phi, h)
  endif else begin
     phi = dblarr(n)
     h = dblarr(n)
     phi_j = 0.d0
     h_j = 0.d0
     for j=0l,n-1 do begin
        result = call_external(dll, 'geodetic_arg_', $
           r_loc[j], z_loc[j], phi_j, h_j)
        phi[j] = phi_j
        h[j] = h_j
     endfor
  endelse
;
  return
  end