jplephinterp.pro

basic median filter simulation

  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, $
                  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 
  ;;1000 = TDB to TDT offset (s), returned in X component
  ;; This numbering scheme is 1-relative, to be consistent with the
  ;; Fortran version.  (units are seconds; derivative units are seconds/day)

  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)
          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:
  ;; -------------------------------------------------------
  ;; 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
          end
          'AU': begin
              au = info.au*info.c/1000d
              x = x / au
              y = y / au
              z = z / au
          end
          'LT-S': begin
              c = info.c / 1000d
              x = x / c
              y = y / c
              z = z / 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)
          end
          'KM/S': begin
              vx = vx * (1D/86400D)
              vy = vy * (1D/86400D)
              vz = vz * (1D/86400D)
          end
          'LT-S/S': begin
              c = info.c / 1000D
              vx = vx / (c*86400D)
              vy = vy / (c*86400D)
              vz = vz / (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
          end
          ELSE: message, 'ERROR: Unrecognized velocity units "'+vu+'"'
      endcase
  endif

CLEAN_RETURN:  
  return
end