jplephinterp.pro

basic median filter simulation

;+
; NAME:
;   JPLEPHINTERP
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   [e-mail]
;   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
;
;   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' - light seconds per second
;                 '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
;
;  $Id: jplephinterp.pro,v 1.17 2009/01/02 17:49:59 craigm Exp $
;
;-
; Copyright (C) 2001, 2002, 2004, 2008, 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.17 2009/01/02 17:49:59 craigm 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