;+
; NAME:
;   JPLEPHINTERP
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;   UPDATED VERSIONs can be found on my WEB PAGE: 
;      http://cow.physics.wisc.edu/~craigm/idl/idl.html
;
; PURPOSE:
;   Interpolate position and motion of planetary bodies (JPL Ephemeris)
;
; MAJOR TOPICS:
;   Planetary Orbits, Interpolation
;
; CALLING SEQUENCE:
;   JPLEPHINTERP, INFO, RAWDATA, T, X, Y, Z, [VX, VY, VZ, /EARTH, /SUN,
;         OBJECTNAME=, CENTER=, TBASE=, POSUNITS=, VELUNITS= ]
;
; DESCRIPTION:
;
;   JPLEPHINTERP interpolates the JPL DE200 or DE405 planetary
;   ephemeris to find the positions and motions of planetary bodies.
;
;   This routine is the second stage of a two-stage process to
;   interpolate the JPL ephemeris.  In this first stage, the file is
;   opened using JPLEPHREAD, and the relevant portions of the table
;   are read and stored into the two variables INFO and RAWDATA.  In
;   the second stage, the user actually interpolates the ephemeris for
;   the desired bodies and to the desired ephemeris time using
;   JPLEPHINTERP.
;
;   The only independent variable which must be specified is T, the
;   ephemeris time.  For low to moderate accuracy applications, T is
;   simply the conventional calendar date, expressed in Julian days.
;   See below for high precision applications.
;
;   Upon output, the position components of the desired body are
;   returned in parameters X, Y and Z, and if requested velocity
;   components are returned in parameters VX, VY and VZ.  Coordinates
;   are referred to the ephemeris's coordinate system: FK5 for
;   JPL-DE200 and ICRS for JPL-DE405.  By default, the origin of
;   coordinates is the solar system barycenter (SSB), unless another
;   origin is selected using the CENTER keyword.
;
;   Users must set the VELOCITY keyword to generate body velocities.
;   By default they are not generated.
;
;   Users can select the desired body by using either the EARTH or SUN
;   keywords, or the OBJECTNAME keyword.
;
;   By default, positions are returned in units of KM and velocities
;   in units of KM/DAY.  However, the output units are selectable by
;   setting the POSUNITS and VELUNITS keywords.
;
; High Precision Applications
;
;   If the required precision is finer than a few hundred meters, the
;   user must be aware that the formal definition of the ephemeris
;   time is the coordinate time of a clock placed at the solar system
;   barycenter (SSB).  If the user's time is measured by a clock
;   positioned elsewhere, then various corrections must be applied.
;   Usually, the most significant correction is that from the
;   geocenter to the SSB (see Fairhead & Bretagnon 1990; Fukushima
;   1995).  Not applying this correction creates an error with
;   amplitude ~170 nano-light-seconds ( = 50 m) on the earth's
;   position.  (see TDB2TDT)
;
;   For high precision, the user should also specify the TBASE
;   keyword.  TBASE should be considered a fixed epoch with respect to
;   which T is measured; T should be small compared to TBASE.
;   Internally, subtraction of large numbers occurs with TBASE first,
;   so truncation error is minimized by specifying TBASE.
;
; Nutations and Librations
;
;   This routine also provides information about earth nutations and
;   lunar librations, which are stored in the JPL ephemeris tables.
;   The POSUNITS and VELUNITS keywords do not affect these
;   computations.
;
;   Lunar librations in the form of three Euler angles are returned in
;   X, Y, Z, in units of radians, and their time derivatives are
;   returned in VX, VY, and VZ in units of radians per day.
;
;   The earth nutation angles psi (nutation in longitude) and epsilon
;   (nutation in obliquity) are returned in X and Y, in units of
;   radians.  Their time derivatives are returned in VX and VY
;   respectively.  The quantities returned in Z and VZ are undefined.
;
; Verification
;
;   The precision routine has been verified using JPLEPHTEST, which is
;   similar to the original JPL program EPHTEST.  For years 1950 to
;   2050, JPLEPHINTERP reproduces the original JPL ephemeris to within
;   1 centimeter.
;
; Custom Ephemerides
;
;   It is possible to make custom ephemerides using JPLEPHMAKE, or to
;   augmented an existing ephemeris with additional data.  In the
;   former case JPLEPHINTERP should automatically choose the correct
;   object from the table and interpolate it appropriately.
;
;   For augmented ephemerides, the object can be specified by name,
;   which works as expected, or by number, which has a special
;   behavior.  For augmented files only, the new objects begin at
;   number 100.
;
;
; PARAMETERS: 
;
;   INFO - structure returned by JPLEPHREAD.  Users should not modify
;          this structure.
;
;   RAWDATA - raw data array returned by JPLEPHREAD.  Users should not
;             modify this data array.
;
;   T - ephemeris time(s) of interest, relative to TBASE (i.e. the
;       actual interpolation time is (T+TBASE)).  May be a scalar or
;       vector.
;
;   X, Y, Z - upon return, the x-, y- and z-components of the body
;             position are returned in these parameters.  For
;             nutations and librations see above.
;
;   VX, VY, VZ - upon return, the x-, y- and z-components of the body
;                velocity are returned in these parameters, if the
;                VELOCITY keyword is set.  For nutations and
;                librations see above.
;
;
; KEYWORD PARAMETERS:
;
;   EARTH, SUN - set one of these keywords if the desired body is the
;                earth or the sun.  One of EARTH, SUN or OBJECTNAME
;                must be specified.
;
;   OBJECTNAME - a scalar string or integer, specifies the planetary
;                body of interest.  May take any one of the following
;                integer or string values.
;
;                   1 - 'MERCURY'     9 - 'PLUTO'
;                   2 - 'VENUS'      10 - 'MOON'  (earth's moon)
;                   3 - 'EARTH'      11 - 'SUN'
;                   4 - 'MARS'       12 - 'SOLARBARY' or 'SSB' (solar system barycenter)
;                   5 - 'JUPITER'    13 - 'EARTHBARY' or 'EMB' (earth-moon barycenter)
;                   6 - 'SATURN'     14 - 'NUTATIONS' (see above)
;                   7 - 'URANUS'     15 - 'LIBRATIONS' (see above)
;                   8 - 'NEPTUNE' 
;
;                For custom ephemerides, the user should specify the
;                object name or number.
;
;                For augmented ephemerides, the user should specify
;                the name.  If the number is specified, then numbers
;                1-15 have the above meanings, and new objects are
;                numbered starting at 100.
;
;   CENTER - a scalar string or integer, specifies the origin of
;            coordinates.  See OBJECTNAME for allowed values.
;            Default: 12 (Solar system barycenter)
;
;   VELOCITY - if set, body velocities are generated and returned in
;              VX, VY and VZ.
;              Default: unset (no velocities)
;
;   POSUNITS - a scalar string specifying the desired units for X, Y,
;              and Z.  Allowed values:
;                 'KM' - kilometers  (default)
;                 'CM' - centimeters
;                 'AU' - astronomical units
;                 'LT-S' - light seconds
;               If angles are requested, this keyword is ignored and
;               the units are always 'RADIANS'.
;
;   VELUNITS - a scalar string specifying the desired units for VX, VY
;              and VZ.  Allowed values:
;                 'KM/DAY' - kilometers per day  (default)
;                 'KM/S' - kilometers per second
;                 'CM/S' - centimeters per second
;                 'LT-S/S' or 'V/C' - light seconds per second or
;                     unitless ratio with speed of light, V/C
;                 'AU/DAY' - astronomical units per day
;
;   TBASE - a scalar or vector, specifies a fixed epoch against wich T
;           is measured.  The ephemeris time will be (T+TBASE).  Use
;           this keyword for maximum precision.
;
;
; EXAMPLE:
;
;   Find position of earth at ephemeris time 2451544.5 JD.  Units are
;   in Astronomical Units.
;
;   JPLEPHREAD, 'JPLEPH.200', pinfo, pdata, [2451544D, 2451545D]
;
;   JPLEPHINTERP, pinfo, pdata, 2451544.5D, xearth, yearth, zearth, $
;                 /EARTH, posunits='AU'
;     
;
; REFERENCES:
;
;   AXBARY, Arnold Rots.
;      ftp://heasarc.gsfc.nasa.gov/xte/calib_data/clock/bary/
;
;   HORIZONS, JPL Web-based ephermis calculator (Ephemeris DE406)
;      http://ssd.jpl.nasa.gov/horizons.html
;   
;   Fairhead, L. & Bretagnon, P. 1990, A&A, 229, 240
;
;   Fukushima, T. 1995, A&A, 294, 895
;
;   Standish, E.M. 1982, "Orientation of the JPL Ephemerides,
;      DE200/LE200, to the Dynamical Equinox of J2000", Astronomy &
;      Astrophysics, vol. 114, pp. 297-302.
;
;   Standish, E.M.: 1990, "The Observational Basis for JPL's DE200,
;      the planetary ephemeris of the Astronomical Almanac", Astronomy
;      & Astrophysics, vol. 233, pp. 252-271.    
;
; SEE ALSO
;   JPLEPHREAD, JPLEPHINTERP, JPLEPHTEST, TDB2TDT, JPLEPHMAKE
;   
; MODIFICATION HISTORY:
;   Written and Documented, CM, Jun 2001
;   Corrected bug in name conversion of NUTATIONS and LIBRATIONS, 18
;     Oct 2001, CM
;   Added code to handle custom-built ephemerides, 04 Mar 2002, CM
;   Fix bug in evaluation of velocity (only appears in highest order
;     polynomial term); JPLEPHTEST verification tests still pass;
;     change is of order < 0.5 cm in position, 22 Nov 2004, CM
;   Perform more validity checking on inputs; and more informative
;     outputs, 09 Oct 2008, CM
;   Allow SSB and EMB as shortcuts for solar system and earth-moon
;     bary center, 15 Oct 2008, CM
;   TBASE now allowed to be a vector or scalar, 01 Jan 2009, CM
;   VELFAC keyword gives scale factor between POSUNITS and VELUNITS, 
;     12 Jan 2009, CM
;   Add option VELUNITS='V/C' for unitless ratio with speed of light,
;     2012-10-02, CM
;
;  $Id: jplephinterp.pro,v 1.19 2012/10/02 11:32:59 cmarkwar Exp $
;
;-
; Copyright (C) 2001, 2002, 2004, 2008, 2009, 2012, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy and distribute unmodified copies for
; non-commercial purposes, and to modify and use for personal or
; internal use, is granted.  All other rights are reserved.
;-

pro jplephinterp_calc, info, raw, obj, t, x, y, z, vx, vy, vz, $
                       velocity=vel, tbase=tbase

  ; '$Id: jplephinterp.pro,v 1.19 2012/10/02 11:32:59 cmarkwar Exp $'

  if n_elements(tbase) EQ 0 then tbase = 0D
  ;; Number of coefficients (x3), number of subintervals, num of rows
  nc = info.ncoeff(obj)
  ns = info.nsub(obj)
  dt = info.timedel
  nr = info.jdrows
  jd0 = info.jdlimits(0) - tbase
  jd1 = info.jdlimits(1) - tbase

  ;; Extract coefficient data from RAW
  if obj EQ 11 then begin
      ;; Nutations have two components
      ii1 = info.ptr(obj)-1
      ii2 = ii1 + nc*ns*2L - 1
      coeffs = reform(dblarr(nc,3,ns,nr), nc,3,ns,nr, /overwrite)
      coeffs(0,0,0,0) = reform(raw(ii1:ii2,*),nc,2,ns,nr, /overwrite)
  endif else begin
      ;; All other bodies are done with three components
      ii1 = info.ptr(obj)-1
      ii2 = ii1 + nc*ns*3L - 1
      coeffs = reform(raw(ii1:ii2,*),nc,3,ns,nr, /overwrite)
  endelse

  ;; Decide which interval and subinterval we are in
  tint = (t-jd0)/dt        ;; Interval number (real)
  ieph = floor(tint)       ;; Interval number (index = int)
  tint = (tint-ieph)*ns    ;; Subinterval number (real)
  nseg = floor(tint)       ;; Subinterval number (index = int)
  ;; Chebyshev "x" (rescaled to range = [-1,1] over subinterval)
  tseg = 2D*(tint - nseg) - 1  

  ;; Below is an optimization.  If the time interval doesn't span an
  ;; ephemeris subinterval, then we can index the coefficient array by
  ;; a scalar, which is much faster.  Otherwise we maintain the full
  ;; vector-level indexing.
  mini = minmax(ieph) & minn = minmax(nseg)
  if mini(0) EQ mini(1) AND minn(0) EQ minn(1) then begin
      ieph = ieph(0)
      nseg = nseg(0)
  endif

  ;; Initialize the first two Chebyshev polynomials, which are P_0 = 1
  ;; and P_1(x) = x
  p0 = 1D
  p1 = tseg
  ;; Initial polynomials for Chebyshev derivatives, V_0 = 0, V_1(x) =
  ;; 1, V_2(x) = 4*x
  v0 = 0D
  v1 = 1D
  v2 = 4D*tseg
  tt = 2D*temporary(tseg)

  x  = 0D & y  = 0D & z  = 0D
  vx = 0D & vy = 0D & vz = 0D
  i0 = ieph*0 & i1 = i0 + 1 & i2 = i1 + 1

  ;; Compute Chebyshev functions two at a time for efficiency
  for i = 0, nc-1, 2 do begin
      if i EQ nc-1 then begin
          p1 = 0
          v1 = 0
      endif
      ii = i0 + i
      jj = i0 + ((i+1) < (nc-1))
      
      x = x + coeffs(ii,i0,nseg,ieph)*p0 + coeffs(jj,i0,nseg,ieph)*p1
      y = y + coeffs(ii,i1,nseg,ieph)*p0 + coeffs(jj,i1,nseg,ieph)*p1
      z = z + coeffs(ii,i2,nseg,ieph)*p0 + coeffs(jj,i2,nseg,ieph)*p1

      if keyword_set(vel) then begin
          vx = vx + coeffs(ii,i0,nseg,ieph)*v0 + coeffs(jj,i0,nseg,ieph)*v1
          vy = vy + coeffs(ii,i1,nseg,ieph)*v0 + coeffs(jj,i1,nseg,ieph)*v1
          vz = vz + coeffs(ii,i2,nseg,ieph)*v0 + coeffs(jj,i2,nseg,ieph)*v1

          ;; Advance to the next set of Chebyshev polynomials. For
          ;; velocity we need to keep the next orders around
          ;; momentarily.
          p2 = tt*p1 - p0
          p3 = tt*p2 - p1
          v2 = tt*v1 - v0 + 2*p1
          v3 = tt*v2 - v1 + 2*p2
          
          p0 = temporary(p2) & p1 = temporary(p3)
          v0 = temporary(v2) & v1 = temporary(v3)
      endif else begin
          ;; Advance to the next set of Chebyshev polynomials.  For no
          ;; velocity, we can re-use old variables.
          p0 = tt*p1 - temporary(p0)
          p1 = tt*p0 - temporary(p1)
      endelse
  endfor

  if keyword_set(vel) then begin
      vfac = 2D*ns/dt
      vx = vx * vfac
      vy = vy * vfac
      vz = vz * vfac
  endif

  return
end

pro jplephinterp_denew, info, raw, obj, t, x, y, z, vx, vy, vz, $
                        velocity=vel, tbase=tbase

  if n_elements(tbase) EQ 0 then tbase = 0D
  dt = info.timedel
  nr = info.jdrows
  jd0 = info.jdlimits(0)
  jd1 = info.jdlimits(1)
  c = info.c / 1000D
  cday = 86400D*info.c/1000D

  ;; Renormalize to fractional and whole days, so fractional
  ;; component is between -.5 and +.5, as needed by barycentering
  ;; approximation code.
  ti  = round(t)      ;; Delta Time: integer
  tbi = round(tbase)  ;; Base: integer
  
  tc = ti + tbi             ;; Total time: integer
  tt = (t-ti) + (tbase-tbi) ;; Total time: fractional

  tc = tc + round(tt)       ;; Re-round: integer
  tt = tt - round(tt)       ;; Re-round: fractional
  t2 = tt*tt                ;; Quadratic and cubic terms
  t3 = t2*tt

  ieph = tc - round(jd0)
  ;; Below is an optimization.  If the time interval doesn't span an
  ;; ephemeris subinterval, then we can index the coefficient array by
  ;; a scalar, which is much faster.  Otherwise we maintain the full
  ;; vector-level indexing.
  mini = minmax(ieph)
  if mini(0) EQ mini(1) then ieph = ieph(0)

  if obj EQ 3 then begin
      ;; Earth, stored as Taylor series coefficients per day
      x = (raw(0,ieph) + raw(3,ieph)*tt + 0.5D*raw(6,ieph)*t2 + $
           (raw(9,ieph)/6D)*t3)
      y = (raw(1,ieph) + raw(4,ieph)*tt + 0.5D*raw(7,ieph)*t2 + $
           (raw(10,ieph)/6D)*t3)
      z = (raw(2,ieph) + raw(5,ieph)*tt + 0.5D*raw(8,ieph)*t2 + $
           (raw(11,ieph)/6D)*t3)
      if keyword_set(vel) then begin
          vx = raw(3,ieph) + raw(6,ieph)*tt + 0.5D*raw(9 ,ieph)*t2
          vy = raw(4,ieph) + raw(7,ieph)*tt + 0.5D*raw(10,ieph)*t2
          vz = raw(5,ieph) + raw(8,ieph)*tt + 0.5D*raw(11,ieph)*t2
      endif
      x = reform(x, /overwrite)
      y = reform(y, /overwrite)
      z = reform(z, /overwrite)
 
  endif else if obj EQ 11 then begin
      ;; Sun, stored as daily components only
      
      x = reform(raw(12,ieph) + tt*0)
      y = reform(raw(13,ieph) + tt*0)
      z = reform(raw(14,ieph) + tt*0)
      if keyword_set(vel) then $
        message, 'ERROR: DENEW format does not provide solar velocity'

  endif else if obj EQ 1000 then begin

      tt = t - (jd0+jd1)/2D
      x = spl_interp(raw(15,*), raw(16,*), raw(17,*), tt)
      return
      
  endif else begin
      message, 'ERROR: DENEW format does not contain body '+strtrim(obj,2)
  endelse
end

pro jplephinterp, info, raw, t, x, y, z, vx, vy, vz, earth=earth, sun=sun, $
                  objectname=obj0, velocity=vel, center=cent, tbase=tbase, $
                  posunits=outunit0, velunits=velunit0, $
                  pos_vel_factor=velfac, $
                  xobjnum=objnum, decode_obj=decode

  if n_params() EQ 0 then begin
      message, 'USAGE: JPLEPHINTERP, info, rawdata, teph, x, y, z, '+$
        'vx, vy, vz, OBJECTNAME="body", /VELOCITY, CENTER="body", '+$
        'POSUNITS="units", VELUNITS="units", /EARTH, /SUN', /info
      return
  endif
  
  ;; The numbering convention for ntarg and ncent is:
  ;;   1 = Mercury            8 = Neptune
  ;;   2 = Venus              9 = Pluto
  ;;   3 = Earth             10 = Moon
  ;;   4 = Mars              11 = Sun
  ;;   5 = Jupiter           12 = Solar system barycenter
  ;;   6 = Saturn            13 = Earth-Moon barycenter
  ;;   7 = Uranus            14 = Nutations (longitude and obliquity; untested)
  ;;                         15 = Librations 
  ;; This numbering scheme is 1-relative, to be consistent with the
  ;; Fortran version.  (units are seconds; derivative units are seconds/day)
  ;;1000 = TDB to TDT offset (s), returned in X component

  sz = size(info)
  if sz(sz(0)+1) NE 8 then message, 'ERROR: INFO must be a structure'
  if ((info.format NE 'JPLEPHMAKE') AND $
      (info.format NE 'BINEPH2FITS') AND $
      (info.format NE 'DENEW')) then begin
      message, 'ERROR: ephemeris type "'+info.format+'" is not recognized'
  endif

  ;; Handle case of custom ephemerides
  if info.format EQ 'JPLEPHMAKE' then begin
      if n_elements(obj0) GT 0 then begin
          sz = size(obj0)
          if sz(sz(0)+1) EQ 7 then begin
              obj = strupcase(strtrim(obj0(0),2))
              wh = where(info.objname EQ obj, ct)
              if ct EQ 0 then $
                message, 'ERROR: '+obj+' is an unknown object'
              obj = wh(0) + 1
          endif else begin
              obj = floor(obj0(0))
              if obj LT 1 OR obj GT n_elements(info.objname) then $
                message, 'ERROR: Numerical OBJNAME is out of bounds'
          endelse

          ;; Interpolate the ephemeris here
          jplephinterp_calc, info, raw, obj-1, t, velocity=vel, $
            tbase=tbase, x, y, z, vx, vy, vz

          goto, COMPUTE_CENTER
      endif
      message, 'ERROR: Must specify OBJNAME for custom ephemerides'
  endif


  ;; ----------------------------------------------------------
  ;; Determine which body or system we will compute
  if n_elements(obj0) GT 0 then begin
      sz = size(obj0)
      if sz(sz(0)+1) EQ 7 then begin
          obj = strupcase(strtrim(obj0(0),2))
          case obj of 
              'EARTH':      obj = 3
              'SOLARBARY':  obj = 12
              'SSB':        obj = 12
              'EARTHBARY':  obj = 13
              'EMB':        obj = 13
              'NUTATIONS':  obj = 14
              'LIBRATIONS': obj = 15
              'TDB2TDT':    obj = 1000
              ELSE: begin
                  wh = where(info.objname EQ obj, ct)
                  if ct EQ 0 then $
                    message, 'ERROR: '+obj+' is an unknown object'
                  obj = wh(0) + 1
                  if obj GT 11 then obj = obj + 100 - 14
              end
          endcase
      endif else begin
          obj = floor(obj0(0))
      endelse
  endif else begin
      if NOT keyword_set(earth) AND NOT keyword_set(sun) then $
        message, 'ERROR: Must specify OBJNAME, EARTH or SUN'
  endelse
  if keyword_set(earth) then obj = 3
  if keyword_set(sun)   then obj = 11

  ;; If the caller is merely asking us to decode the objectnumber,
  ;; then return it now.
  objnum = obj
  if keyword_set(decode) then return
  
  jdlimits = info.jdlimits

  ;; -------------------------------------------------------
  ;; Handle case of de200_new.fits format
  if info.format EQ 'DENEW' then begin
      if objnum NE 3 AND objnum NE 11 AND objnum NE 1000 then $
        message, 'ERROR: DENEW ephemeris table does not support body #'+$
        strtrim(objnum,2)
      
      jplephinterp_denew, info, raw, objnum, t, x, y, z, vx, vy, vz, $
        velocity=vel, tbase=tbase

      if objnum GE 1000 then return
      goto, DO_UNIT
  endif

  ;; -------------------------------------------------------
  ;; Otherwise, construct the ephemeris using the Chebyshev expansion
  case obj of
      3: begin ;; EARTH (translate from earth-moon barycenter to earth)
          ;; Interpolate the earth-moon and moon ephemerides
          jplephinterp_calc, info, raw, 2, velocity=vel, tbase=tbase, $
            t, xem, yem, zem, vxem, vyem, vzem
          jplephinterp_calc, info, raw, 9, velocity=vel, tbase=tbase, $
            t, xmo, ymo, zmo, vxmo, vymo, vzmo
          emrat = info.emrat
          
          ;; Translate from the earth-moon barycenter to earth
          x = xem - emrat * xmo
          y = yem - emrat * ymo
          z = zem - emrat * zmo
          if keyword_set(vel) then begin
              vx = vxem - emrat * vxmo
              vy = vyem - emrat * vymo
              vz = vzem - emrat * vzmo
          endif
          
      end

      10: begin ;; MOON (translate from earth-moon barycenter to moon)
          jplephinterp_calc, info, raw, 9, t, velocity=vel, tbase=tbase, $
            x, y, z, vx, vy, vz
          ;; Moon ephemeris is geocentered.  If the center is
          ;; explicitly earth then return immediately.  Otherwise
          ;; follow the standard path via the solar barycenter.
          if n_elements(cent) GT 0 then begin
              jplephinterp, info, objectname=cent(0), tbase=tbase, $
                xobjnum=cent1, /decode_obj
              if cent1 EQ 3 then goto, DO_UNIT
          endif

          ;; Use solar barycenter via the earth-moon barycenter
          jplephinterp_calc, info, raw, 2, t, velocity=vel, tbase=tbase, $
            xem, yem, zem, vxem, vyem, vzem
          emrat = 1d - info.emrat
          x = xem + emrat * x
          y = yem + emrat * y
          z = zem + emrat * z
          if keyword_set(vel) then begin
              vx = vxem + emrat * vx
              vy = vyem + emrat * vy
              vz = vzem + emrat * vz
          endif              
      end

      12: begin ;; SOLARBARY
          x = t*0D & y = x & z = x
          vx = x   & vy = x & vz = x
      end

      13: begin ;; EARTHBARY
          jplephinterp_calc, info, raw, 2, velocity=vel, tbase=tbase, $
            t, x, y, z, vx, vy, vz
      end
      
      14: begin ;; NUTATIONS
          ;; X = PSI, Y = EPSILON, VX = PSI DOT, VY = EPSILON DOT 
          jplephinterp_calc, info, raw, 11, velocity=vel, tbase=tbase, $
            t, x, y, z, vx, vy, vz
          goto, CLEAN_RETURN
      end

      15: begin ;; LIBRATIONS
          jplephinterp_calc, info, raw, 12, velocity=vel, tbase=tbase, $
            t, x, y, z, vx, vy, vz
          goto, CLEAN_RETURN
      end

      1000: begin ;; TDT to TDB conversion
          x = tdb2tdt(t, deriv=vx, tbase=tbase)
          if n_elements(velunit0) GT 0 then begin
             ;; Special case of unit conversion when user asks for 
             ;; "per second"
             if strpos(strupcase(velunit0(0)),'/S') GE 0 then $
                vx = vx / 86400d
          endif

          goto, CLEAN_RETURN
      end

      else: begin
          ;; Default objects are derived from the index OBJNUM
          if obj GE 1 AND obj LE 11 then begin
              RESTART_OBJ:
              jplephinterp_calc, info, raw, obj-1, t, velocity=vel, $
                tbase=tbase, $
                x, y, z, vx, vy, vz
          endif else begin
              if info.edited AND obj GT 11 then begin
                  ;; Handle case of edited JPL ephemerides - they
                  ;; start at a value of 100, so shift them to the end
                  ;; of the JPL ephemeris columns
                  obj = obj - 100 + 14
                  if obj LE n_elements(info.objname) then $
                    goto, RESTART_OBJ
              endif
              message, 'ERROR: body '+strtrim(obj,2)+' is not supported'
          endelse
      end
  endcase

  ;; -------------------------------------------------------
  ;; Compute ephemeris of center, and compute displacement vector
  COMPUTE_CENTER:
  if n_elements(cent) GT 0 then begin
      jplephinterp, info, raw, t, x0, y0, z0, vx0, vy0, vz0, tbase=tbase, $
        objectname=cent, velocity=vel, posunits='KM', velunits='KM/DAY'
      x = temporary(x) - temporary(x0)
      y = temporary(y) - temporary(y0)
      z = temporary(z) - temporary(z0)
      if keyword_set(vel) then begin
          vx = temporary(vx) - temporary(vx0)
          vy = temporary(vy) - temporary(vy0)
          vz = temporary(vz) - temporary(vz0)
      endif
  endif

  DO_UNIT:
  
  velfac = 1d

  ;; -------------------------------------------------------
  ;; Convert positional units
  if n_elements(outunit0) GT 0 then begin
      pu = strupcase(strtrim(outunit0(0),2))
      case pu of
          'KM': km = 1 ;; Dummy statement
          'CM': begin
              x = x * 1D5
              y = y * 1D5
              z = z * 1D5
              velfac = velfac * 1D5
          end
          'AU': begin
              au = info.au*info.c/1000d
              x = x / au
              y = y / au
              z = z / au
              velfac = velfac / au
          end
          'LT-S': begin
              c = info.c / 1000d
              x = x / c
              y = y / c
              z = z / c
              velfac = velfac / c
          end
          ELSE: message, 'ERROR: Unrecognized position units "'+pu+'"'
      endcase
  endif

  ;; -------------------------------------------------------
  ;; Convert velocity units
  if n_elements(velunit0) GT 0 AND keyword_set(vel) then begin
      vu = strupcase(strtrim(velunit0(0),2))
      case vu of 
          'CM/S': begin
              vx = vx * (1D5/86400D)
              vy = vy * (1D5/86400D)
              vz = vz * (1D5/86400D)
              velfac = velfac / (1D5/86400D)
          end
          'KM/S': begin
              vx = vx * (1D/86400D)
              vy = vy * (1D/86400D)
              vz = vz * (1D/86400D)
              velfac = velfac / (1D/86400D)
          end
          'LT-S/S': begin ;; Light seconds per second
              c = info.c / 1000D
              vx = vx / (c*86400D)
              vy = vy / (c*86400D)
              vz = vz / (c*86400D)
              velfac = velfac / (c*86400D)
          end
          'V/C': begin ;; Unitless ratio V/C (same as LT-S/S
              c = info.c / 1000D
              vx = vx / (c*86400D)
              vy = vy / (c*86400D)
              vz = vz / (c*86400D)
              velfac = velfac / (c*86400D)
          end
          'KM/DAY': km = 1 ;; Dummy statement
          'AU/DAY': begin
              au = info.au*info.c/1000d
              vx = vx / au
              vy = vy / au
              vz = vz / au
              velfac = velfac * au
          end
          ELSE: message, 'ERROR: Unrecognized velocity units "'+vu+'"'
      endcase
  endif

CLEAN_RETURN:  
  return
end