pro spacecraft_pv, spacecraft_id, jdn, pv, coorID=coorID, tle_path=path, $
    pv_uncertainty=dpv, interpolation_fraction=f, status=status
;
; Given a Norad (US SpaceCom) spacecraft id, and double-precision
; Julian Day Number and fraction, this procedure returns the ECI
; position and velocity of the spacecraft.
;
; B. Knapp, 2003-02-14
;
; The input argument  SPACECRAFT_ID  must be an integer scalar.
; The input argument  JDN  may be a scalar or array, and the output
; arguments  PV  and DPV will agree in kind.  The optional input
; argument tle_path specifies the full path of the directory in
; which the Norad TLE files are located; if it is omitted, the
; default path is /tools/knapp/dat/tle/.
;
; The coordinates used for the outputs may be either spherical (coorid=1)
; or rectangular (coorid=2, default).  Angles are in degrees and
; distances in km.  Each output position is a 6-vector with elements
; as follows:
;
;      I   SPHERICAL                            RECTANGULAR
;
;      0   Right Ascension (deg, 0 to 360)      X (km)
;      1   Declination (deg, -90 to +90)        Y (km)
;      2   R (km)                               Z (km)
;      3   dRA/dT (deg/sec)                     dX/dT (km/sec)
;      4   dDec/dT (deg/sec)                    dY/dT (km/sec)
;      5   dR/dT (km/sec)                       dZ/dT (km/sec)
;
; The output argument  STATUS  is returned non-zero if an error occurred.
;
; Print usage?
  if n_params() lt 2 then begin
    print,' '
    print,' Given a Norad (US Spacecom) spacecraft ID, and double-'
    print,' precision Julian Day Number and fraction,  this procedure'
    print,' returns the position and velocity of the spacecraft. Example'
    print,' spacecraft IDs are 21701 (UARS) or 27651 (SORCE).'
    print,' '
    print,' The coordinates used for the outputs may be either spherical'
    print,' (coorid=1) or rectangular (coorid=2, default).  Angles are in'
    print,' degrees and distances in km. Each output position is a 6-vector'
    print,' with elements as follows:'
    print,' '
    print,'      I   SPHERICAL                   RECTANGULAR'
    print,' '
    print,'      0   RA (deg, 0 to 360)          X (km)'
    print,'      1   Dec (deg, -90 to +90)       Y (km)'
    print,'      2   R (km)                      Z (km)'
    print,'      3   dRA/dT (deg/sec)            dX/dT (km/sec)'
    print,'      4   dDec/dT (deg/sec)           dY/dT (km/sec)'
    print,'      5   dR/dT (km/sec)              dZ/dT (km/sec)'
    print,' '
    print,' The output argument  STATUS  is returned non-zero if an error'
    print,' occurred.'
    print,' '
    print,'     spacecraft_pv, spacecraft_id, jdn, pv, $'
    print,'        coorID=coordID, $ ;1=spherical, 2=rectangular (default)'
    print,'        tle_path=path,  $ ;default is /tools/knapp/dat/tle/'
    print,'        pv_uncertainty = dpv, $'
    print,'        interpolation_fraction = f, $'
    print,'        status=status'
    print,' '
    return
  endif
;
; Keyword arguments?
  if n_elements(coorid) eq 0 then coorid = 2
  if n_elements(path) eq 0 then path = '/tools/knapp/dat/tle/'
  path=path_sep()+(STRJOIN(STRSPLIT(path, '\/',/EXTRACT), path_sep()))+path_sep()
;
; Cast inputs to types required by Fortran native routine
  id = long( spacecraft_id )
  djd = double( jdn )
  bpath = bytarr(1024)
  bpath[0] = byte(path)
  pathlen = strlen(path)
;
; Define outputs
  status = 0L
;

; Call the Fortran service
  n = n_elements( jdn )
  rad2deg = 180.d0/!dpi
  earth_radius = 6378.135d0
  if n eq 1 then begin
    xyz = dblarr(6)
    dxyz = dblarr(6)
    f = 0.d0
    result = call_external( !astron_lib_path, 'spacecraft_pv_arg_', $
      id, pathlen, bpath, djd[0], xyz, dxyz, f , status )
    if coorid eq  2 then begin
      pv = xyz*earth_radius
      pv[3:5] = pv[3:5]/60.d0
      dpv = dxyz*earth_radius
      dpv[3:5] = dpv[3:5]/60.d0
    endif else begin
      sphxyz,-1,pv,xyz,dpv,dxyz
      pv[0:1] = pv[0:1]*rad2deg
      pv[2] = pv[2]*earth_radius
      pv[3:4] = pv[3:4]*rad2deg/60.d0
      pv[5] = pv[5]*earth_radius/60.d0
      dpv[0:1] = dpv[0:1]*rad2deg
      dpv[2] = dpv[2]*earth_radius
      dpv[3:4] = dpv[3:4]*rad2deg/60.d0
      dpv[5] = dpv[5]*earth_radius/60.d0
    endelse
;
  endif else begin
    pv = dblarr(6,n)
    xyz = dblarr(6)
    dpv = dblarr(6,n)
    dxyz = dblarr(6)
    f = dblarr(n)
    fj = 0.d0
    for j=0l,n-1 do begin
      result = call_external( !astron_lib_path, 'spacecraft_pv_arg_', $
        id, pathlen, bpath, djd[j], xyz, dxyz, fj, status )
      if status ne 0 then return
      f[j] = fj
      if coorid eq 2 then begin
        pv[0,j] = xyz*earth_radius
        pv[3:5,j] = pv[3:5,j]/60.d0
        dpv[0,j] = dxyz*earth_radius
        dpv[3:5,j] = dpv[3:5,j]/60.d0
      endif else begin
        sphxyz,-1,sph,xyz,dsph,dxyz
        pv[0:1,j] = sph[0:1]*rad2deg
        pv[2,j] = sph[2]*earth_radius
        pv[3:4,j] = sph[3:4]*rad2deg/60.d0
        pv[5,j] = sph[5]*earth_radius/60.d0
        dpv[0:1,j] = dsph[0:1]*rad2deg
        dpv[2,j] = dsph[2]*earth_radius
        dpv[3:4,j] = dsph[3:4]*rad2deg/60.d0
        dpv[5,j] = dsph[5]*earth_radius/60.d0
      endelse
    endfor
  endelse
;
  return
  end