pro ephem, jdn, body, pv, status, posnid=pid, framid=fid, coorid=cid
;
; Given a double-precision Julian date and a solar system body
; identifier, this procedure returns the geocentric or heliocentric
; position and velocity of the body.
;
; B. Knapp, 1999-06-25, 2001-01-02, 2001-04-24 (32-bit & 64-bit)
;
; The input argument  JDN  may be a scalar or array, and the output
; argument  PV  will agree in kind.  The input argument  BODY
; must be an integer scalar identifying the solar system body
; (0=Sun, 1=Mercury, 2=Venus, 3=Earth, 4=Mars, 5=Jupiter, 6=Saturn,
; 7=Uranus, 8=Neptune, 9=Pluto, 10=Moon).
;
; The output position may be the true (posnid=1) or apparent (posnid=2,
; default) position; the reference frame may be Equatorial (framid=1, 
; default), Ecliptic (framid=2) or Heliocentric (framid=3), and the
; coordinates may be either spherical (coorid=1, default) or rectangular
; (coorid=2).  Angles are in degrees and distances in AU for planets or
; km for the Moon.  Each ouput position is a 6-vector with elements
; as follows:
;
;      I   SPHERICAL                   RECTANGULAR
;
;      0   RA or L (deg, 0 to 360)     X (AU)
;      1   Dec or B (deg, -90 to +90)  Y (AU)
;      2   R (AU, or km for Moon)      Z (AU)
;      3   dRA/dT or dL/dT (deg/day)   dX/dT (AU/day)
;      4   dDec/dT or dB/dT (deg/day)  dY/dT (AU/day)
;      5   dR/dT (AU/day or km/day)    dZ/dT (AU/day)
;
; The output argument  STATUS  is returned non-zero if an error
; occurred.
;
; The algorithms used are based on the VSOP87D Theory for the planets
; (P. Bretagnon & G. Francou, Astron. Astropyhs. 202 (1988), 309-315),
; and on the ELP2000-85 Theory for the Moon (M. Chapront-Touze &
; J. Chapront, Astron. Astrophys. 190 (1988), 342-352).  The output
; positions are double precision, and are quite accurate (on the order of
; one arc second for the sun and planets, ten arc seconds for the Moon).
;
; Print usage?
  if n_params() lt 2 then begin
     print, ' '
     print, ' Given a double-precision Julian date and a solar system'
     print, ' body identifier, ephem returns the geocentric or helio-'
     print, ' centric position and velocity of the body (0=Sun, 1='
     print, ' Mercury, 2=Venus, 3=Earth, 4=Mars, 5=Jupiter, 6=Saturn,'
     print, ' 7=Uranus, 8=Neptune, 9=Pluto, 10=Moon)'
     print, ' '
     print, '     ephem, jdn, body, pv, status, $'
     print, '        posnid=pid, $ ; 1=true, 2=apparent (default)'
     print, '        framid=fid, $ ; 1=Equatorial (default), 2=Ecliptic'
     print, '                      ; 3=Heliocentric'
     print, '        coorid=cid    ; 1=spherical (default), 2=rectangular'
     print, ' '
     return
  endif
;
; Keyword arguments?
  if n_elements( pid ) eq 0 then $
     pid_loc = 2L $
  else $
     pid_loc = long( pid )
;
  if n_elements( fid ) eq 0 then $
     fid_loc = 1L $
  else $
     fid_loc = long( fid )
;
  if n_elements( cid ) eq 0 then $
     cid_loc = 1L $
  else $
     cid_loc = long( cid )

; Cast inputs to types required by ephem.for
  id = long( body )
  djd = double( jdn )
;
; Define outputs
  status = 0L
;
; 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')")
;
; Call the Fortran service
  n = n_elements( jdn )
  if n eq 1 then begin
     pv = dblarr(6)
     result = call_external( dll, 'ephem_arg_', $
        djd[0], id, pid_loc, fid_loc, cid_loc, pv , status )
  endif else begin
     pv = dblarr(6,n)
     pvj = dblarr(6)
     for j=0l,n-1 do begin
        result = call_external( dll, 'ephem_arg_', $
           djd[j], id, pid_loc, fid_loc, cid_loc, pvj, status )
        pv[0,j] = pvj
        if status ne 0 then return
     endfor
  endelse
;
  return
  end