;+
; NAME:
;   LITMSOL
;
; 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:
;   Solve the light-time equation between two moving bodies
;
; MAJOR TOPICS:
;   Geometry, Physics
;
; CALLING SEQUENCE:
;   LITMSOL, T1, X1, Y1, Z1, T2, INFO2, RAW2, OBJ2, INFOSUN, RAWSUN, $
;            /RECEIVER, TBASE=, TOLERANCE=, POSUNITS=, MAXITER=, $
;            /NO_SHAPIRO
;
; DESCRIPTION:
;
;  The procedure LITMSOL solves the light time equation between two
;  moving bodies in the solar system.  Given the time and position of
;  reception or transmission of a photon, this equation determines the
;  time of transmission or reception at the other solar system body.
;  Since both bodies may be moving, the equation must be solved
;  iteratively.
;
;  The trajectories of solar system bodies must be described by either
;  a JPL ephemeris, or by a JPL-like ephemeris generated by
;  JPLEPHMAKE.  This routine calls JPLEPHINTERP.
;
;  The user specifies the known time and position of interaction as
;  T1, X1, Y1 and Z1, in units of POSUNITS.  The time of interaction
;  at the other body -- the solution to the light time equation -- is
;  returned as T2.  If the photon was *received* at time T1, then the
;  RECEIVER keyword should be set (in which case the transmission must
;  have occurred in the past).
;
;  Since the solution is iterative, the user may specify a solution
;  tolerance, and a maximum number of iterations.
;
;  If users wish to include the Shapiro time delay, which has a
;  maximum amplitude of approximately 250 usec, they must specify the
;  ephemeris of the Sun (INFOSUN, RAWSUN).  The Shapiro delay is the
;  extra general relativistic delay caused by the Sun's potential.
;
;
; INPUTS:
;
;   T1 - epoch of interaction, in Julian days, in the TDB timescale.
;        (scalar or vector)
;
;   X1, Y1, Z1 - coordinates of interaction, referred to the solar
;                system barycenter, in J2000 coordinates.  Units are
;                described by POSUNITS. (scalar or vector)
;
;   INFO2, RAW2 - ephemeris of other solar system body, returned by
;                 JPLEPHREAD or JPLEPHMAKE.
;
;   INFOSUN, RAWSUN - ephemeris of at least the Sun, as returned by
;                     JPLEPHREAD.  Only used of NO_SHAPIRO is not set.
;
;
; OUTPUTS:
;
;   T2 - upon output, epoch of interaction at the second solar system
;        body, in Julian days, in the TDB timescale.
;
;
; KEYWORD PARAMETERS:
;
;   RECEIVER - if set, then the epoch T1 is a reception of a photon.
;              Otherwise T1 is the epoch of transmission of a photon.
;
;   VX1, VY1, VZ1 - upon input, the body velocity at time T1, in
;                   VELUNITS units.  This information is required only
;                   if the SHAPIRO_DERIV is required.
;
;   X2, Y2, Z2 - upon return, the body position at time T2, in
;                POSUNITS units.
;   VX2, VY2, VZ2 - upon return, the body velocity at time T2, in
;                VELUNITS units.
;
;   TBASE - a fixed epoch time (Julian days) to be added to each value
;           of T1.  Since subtraction of large numbers occurs with
;           TBASE first, the greatest precision is achieved when TBASE
;           is expressed as a nearby julian epoch, T1 is expressed
;           as a small offset from the fixed epoch.  
;           Default: 0
;
;   POSUNITS - the units for positions, one of 'CM', 'KM', 'LT-S' or
;              'AU'.
;              Default: 'CM'
;   VELUNITS - the units for velocities (and Shapiro derivative).
;              Default: POSUNITS+'/S'
;
;   TOLERANCE - the solution tolerance, expressed in POSUNITS.
;               Default: 1000 CM
;
;   ERROR - upon return, a vector giving the estimated error in the
;           solution for each point, expressed in POSUNITS.  This
;           quantity should be less than TOLERANCE unless the number
;           of iterations exceeded MAXITER.
;
;   MAXITER - maximum number of solution iterations to be taken.
;             Default: 5
;   NITER - upon return, contains the actual number of iterations used.
;
;   SHAPIRO_CALC - method of calculating Shapiro delay, a string with
;                  one value of 'NONE', 'DELAY' or 'BOTH'.  NONE means
;                  do not calculate any Shapiro delay values.  DELAY
;                  means calculate Shapiro delay only.  BOTH means
;                  calculate the delay *and* its derivative with
;                  respect to time.  If SHAPIRO_CALC is set to
;                  DELAY or BOTH, then INFOSUN and RAWSUN must be
;                  specified.  If BOTH, then VX1, VY1 and VZ1 must
;                  also be specified. This keyword overrides
;                  NO_SHAPIRO.
;   NO_SHAPIRO - if set, then the Shapiro delay will not be accounted
;                for.  Use SHAPIRO_CALC instead.
;   SHAPIRO_DELAY - upon return, contains the Shapiro delay in
;                   seconds, if SHAPIRO_CALC is set to 'DELAY' or
;                   'BOTH'.
;   SHAPIRO_DERIV - upon return, contains the derivative of the
;                   Shapiro delay, in light seconds per time unit of
;                   velocity (SHAPIRO_CALC must be set to 'BOTH' to
;                   enable this calculation).  Note that you must
;                   supply VX1, VY1 and VZ1 to get the derivative
;                   value.
;
;
; EXAMPLE:
;
;
;
; SEE ALSO:
;
;   JPLEPHREAD, JPLEPHINTERP
;
;
; MODIFICATION HISTORY:
;   Written, 6 May 2002, CM
;   Documented, 12 May 2002, CM
;   Added TGUESS keyword, 29 May 2002, CM
;   Added ERROR and X/Y/ZOFF keywords, 25 Sep 2002, CM
;   Extensive revisions: addition of SHAPIRO_{CALC,DELAY,DERIV}
;     values; input VX1, VY1 and VZ1; output X2, Y2, Z2 and VX2 VY2
;     and VZ2; and VELUNITS keyword, 07 Mar 2007, CM
;   Allow user specified function to interpolate INFO2/RAW2 via
;     INTERP_FUNC keyword, 09 Oct 2008, CM
;
;  $Id: litmsol.pro,v 1.7 2008/10/10 00:50:19 craigm Exp $
;
;-
; Copyright (C) 2002, 2007, 2008, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-

pro litmsol, t1, x1, y1, z1, t2, info2, raw2, obj2, info, raw, $
             tbase=tbase, tolerance=tol0, posunits=posunits0, $
             receiver=receiver, maxiter=maxiter0, no_shapiro=noshap, $
             tguess=tguess, error=diff, niter=i, $
             shapiro_delay=shdelay, shapiro_deriv=shderiv, shapiro_calc=shcalc0, $
             vx1=vx1, vy1=vy1, vz1=vz1, velunits=velunits0, $
             x2=x2, y2=y2, z2=z2, vx2=vx2, vy2=vy2, vz2=vz2, $
             xoffset=xoff, yoffset=yoff, zoffset=zoff, $
             interp_func=intfunc0

  if n_elements(shcalc0) GT 0 then begin
      shcalc = strupcase(strtrim(shcalc0(0),2))
      if shcalc NE 'NONE' AND shcalc NE 'BOTH' AND $
         shcalc NE 'DELAY' then $
        message, 'ERROR: SHAPIRO_CALC value of "'+shcalc+'" was invalid.'
      if keyword_set(noshap) AND shcalc NE 'NONE' then $
        message, 'ERROR: SHAPIRO_CALC and NO_SHAPIRO did not agree.'
  endif else begin
      shcalc = 'DELAY'
      if keyword_set(noshap) then shcalc = 'NONE'
  endelse
      
  ;; Default position and velocity units
  if n_elements(posunits0) EQ 0 then begin
      posunits = 'CM'
  endif else begin
      posunits = strtrim(posunits0(0),1)
  endelse
  if n_elements(velunits0) EQ 0 then begin
      velunits=posunits+'/S'
  endif else begin
      velunits = strtrim(velunits0(0),1)
  endelse

  ;; Default tolerances
  if n_elements(tol0) EQ 0 then begin
      tol = 1000d     ;; 10 m tolerance
      posunits = 'CM'
  endif else begin
      tol = tol0(0)
  endelse

  ;; Default interpolation function: JPLEPHINTERP
  if n_elements(intfunc0) EQ 0 then begin
      intfunc = 'JPLEPHINTERP'
  endif else begin
      intfunc = strtrim(strupcase(intfunc0(0)),2)
  endelse

  if n_elements(xoff) EQ 0 then xoff = 0d
  if n_elements(yoff) EQ 0 then yoff = 0d
  if n_elements(zoff) EQ 0 then zoff = 0d

  if n_elements(maxiter0) EQ 0 then maxiter = 5L $
  else                              maxiter = floor(maxiter0(0))>2
      
  case posunits of 
      'CM':   clight = info2.c*1d2  ;; CM/S
      'KM':   clight = info2.c*1d-3 ;; KM/S
      'LT-S': clight = 1d           ;; LT-S/S
      'AU':   clight = 1d/info2.au  ;; AU/S
  endcase

  ;; Decide whether to compute derivative of shapiro delay
  if shcalc EQ 'BOTH' then begin
      calcvel = 1

      if n_elements(vx1) EQ 0 OR n_elements(vy1) EQ 0 OR $
        n_elements(vz1) EQ 0 then begin
          message, 'ERROR: You must specify VX1,VY1,VZ1 when computing '+$
            'the Shapiro derivative'
      endif
  endif
  if arg_present(vx2) OR arg_present(vy2) OR arg_present(vz2) then calcvel = 1

  ;; If any Shapiro calculations are required, pre-compute r1dot
  if shcalc NE 'NONE' then begin
      jplephinterp, info, raw, t1, xs, ys, zs, /sun, $
        velocity=calcvel, vxs, vys, vzs, $
        posunits=posunits, velunits=velunits, tbase=tbase
      r1 = sqrt((x1-xs)^2 + (y1-ys)^2 + (z1-zs)^2)
      if shcalc EQ 'BOTH' then $
        r1dot = ((x1-xs)*(vx1-vxs)+(y1-ys)*(vy1-vys)+(z1-zs)*(vz1-vzs))/r1
  endif

  ;; Use TGUESS if provided, otherwise estimate T2 as T1 to begin with
  dt0 = t1*0
  if n_elements(tguess) EQ n_elements(t1) then begin
      t2 = tguess 
      dtold = abs(tguess-t1)*86400d
      if keyword_set(receiver) then dtold = -dtold
  endif else begin
      t2 = t1
      dtold = 0d
  endelse


  ;; ==================== BEGIN ITERATION
  ct = 1L
  i = 0L
  while (ct GT 0) AND (i LT maxiter) do begin

      ;; Position and vel of body at estimated time T2
      call_procedure, intfunc, info2, raw2, t2, x2, y2, z2, objectname=obj2, $
        velocity=calcvel, vx2, vy2, vz2, $
        tbase=tbase, posunits=posunits, velunits=velunits

      ;; Compute distance from T1 to T2 in physical and light-second units
      dr = sqrt((x1-(x2+xoff))^2 + (y1-(y2+yoff))^2 + (z1-(z2+zoff))^2)
      dt = dr / clight

      ;; Add Shapiro delay
      if shcalc NE 'NONE' then begin

          ;; First iteration: recompute time of sun position, since
          ;; the actual sun position at time of closest approach will
          ;; be slightly different
          if i EQ 1 then begin
              ts = t1

              ;; Angle between S/C and sun, as seen from earth
              cos_th = ((x2-x1)*(xs-x1) + (y2-y1)*(ys-y1) + (z2-z1)*(zs-z1))/dr/r1
              whc = where(cos_th GT 0, ctc)
              ;; Correct for the light travel time to the time of
              ;; closest approach
              if ctc GT 0 then ts(whc) = (t1 + r1*cos_th/clight/86400d)(whc)

              ;; Recompute the sun position
              jplephinterp, info, raw, ts, xs, ys, zs, /sun, $
                velocity=calcvel, vxs, vys, vzs, $
                posunits=posunits, velunits=velunits, tbase=tbase
              r1 = sqrt((x1-xs)^2 + (y1-ys)^2 + (z1-zs)^2)
              if shcalc EQ 'BOTH' then $
                r1dot = ((x1-xs)*(vx1-vxs)+(y1-ys)*(vy1-vys)+(z1-zs)*(vz1-vzs))/r1
          endif

          
          ;; Compute Shapiro delay
          dt0 = 2*info.msol
          ;; Body-sun distance
          r2 = sqrt((x2-xs)^2 + (y2-ys)^2 + (z2-zs)^2)
          rsump = (r1+r2+dr+dt0*clight)
          rsumm = (r1+r2-dr+dt0*clight)
          shdelay = dt0*alog(rsump/rsumm)
          dt = dt + shdelay

          ;; Compute derivative of Shapiro delay
          if shcalc EQ 'BOTH' then begin
              r2dot = ((x2-xs)*(vx2-vxs)+(y2-ys)*(vy2-vys)+(z2-zs)*(vz2-vzs))/r2
              drdot = ((x2-x1)*(vx2-vx1)+(y2-y1)*(vy2-vy1)+(z2-z1)*(vz2-vz1))/dr
              shderiv = 2*dt0*(((r1+r2+dt0*clight)*drdot - (r1dot+r2dot)*dr) / $
                               (rsump*rsumm))
          endif
      endif

      ;; If 1 <-> 2 are switched, then reverse sign of DT
      if keyword_set(receiver) then dt = -dt

      ;; Correct T2 and compute how much change from previous iteration
      t2 = t1 + dt/86400d
      diff = abs(dtold-dt)*clight
      wh = where(diff GT tol, ct)

      ;; Prepare for next iteration
      dtold = dt
      i = i + 1
  endwhile

  return
end