;+
; NAME:
; LITMSOL2
;
; 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, Dynamics
;
; CALLING SEQUENCE:
; LITMSOL2, T1, X1, Y1, Z1, T2, $
; FUNC2, INFO2, RAW2, FUNCTARGS=, FUNCTSAVE=, $
; /RECEIVER, TBASE=, TOLERANCE=, POSUNITS=, MAXITER=, $
; LIGHT_TIME=, TGUESS=, ERROR=, NITER=, $
; VX1=, VY1=, VZ1=, $
; X2=, Y2=, Z2=, VX2=, VY2=, VZ2=, $
; METHOD=, $
; DELAY_FUNCTION=, DELAY_ARG1=, DELAY_ARG2=, $
; DELAY_FUNCTARGS=
;
; DESCRIPTION:
;
; The procedure LITMSOL2 solves the light time equation between two
; moving bodies, A and B, in the solar system. Given the time and
; position of reception or transmission of a photon at A, this
; equation determines the time of transmission or reception at the
; other solar system body B. Since both bodies may be moving, the
; equation must be solved iteratively.
;
; The user must know the "A" endpoint of the ray, with time T1 and
; position X1,Y1,Z1. LITMSOL2 solves for the "B" endpoint time and
; position T2 and X2,Y2,Z2 by propagating a light ray from one to the
; other.
;
; The position of the "B" body must be described as an interpolatable
; function. The user function FUNC2 must calculate the position (and
; velocity) of the body at any applicable time T2, in the requested
; units.
;
; By default the body "A" is considered the transmitter and LITMSOL2
; calculates the time at which body "B" receives the ray. However,
; if /RECEIVER is set, then body "A" is considered the receiver, and
; LITMSOL2 calculates the time T2 in the past at which the ray must
; have been transmitted by body "B" in order to be received by "A" at
; time T1.
;
; LITMSOL2 is able to estimate the T2 knowing only the time and
; position at body "A". However, convergence may be faster if the
; TGUESS, METHOD and/or VX1,VY1,VZ1 keywords are used. By default,
; the initial guess for T2 is simply the same as T1. A better
; estimate can be passed in the TGUESS keyword.
;
; If velocity information is available, then LITMSOL2 can use a
; simple linear corrector method in order to speed convergence
; (i.e. Newton's method). The user should pass the velocity
; at time T1 in the VX1,VY1,VZ1 keywords, and METHOD='CORRECTOR'.
;
; The user may also specify a "delay" function which estimates any
; additional light propagation delays along the path based on the
; current estimates of the two ray endpoints. One such delay might
; be the "Shapiro" delay due to general relativity.
;
; Since the solution is iterative, the user may specify a solution
; tolerance, and a maximum number of iterations. An estimate of the
; solution error is returned in the ERROR keyword.
;
; USER FUNCTIONS
;
; The user must supply a function to interpolate the position of the
; body at time T, which is passed in parameter FUNC2. FUNC2, a
; scalar string, is the name of subroutine to call which must compute
; position of body at time T2. The calling convention is the same as
; JPLEPHINTERP, namely,
;
; PRO FUNC2, INFO2, RAW2, T2, X2, Y2, Z2, VX2, VY2, VZ2, $
; VELOCITY=, POSUNITS=, VELUNITS=, SAVE=, ...
;
; The variables INFO2 and RAW2 are described below. The variable T2
; is the requested time (TDB), and the position and velocity must be
; returned in X2,Y2,Z2, VX2,VY2,VZ2, with the requested units. The
; SAVE keyword can designate one keyword whose value will be returned
; to the calling routine. Any other keywords can be passed using the
; _EXTRA calling convention using the FUNCTARGS keyword.
;
; The user may also supply an optional function to compute an
; additional delay. The delay may be a function of the time and
; position of both points "A" and "B". For example, the "Shapiro
; delay" of photons in the solar potential is one such kind of delay.
; The calling convention is,
;
; DELAY = DELAY_FUNCTION(DELAY_ARG1, DELAY_ARG2, $
; T1, X1, Y1, Z1, T2, X2, Y2, Z2, $
; POSUNITS=, TBASE=, ...)
;
; The returned delay must be in seconds, with the sense that a
; positive value of DELAY indicates that the actual light travel time
; is *longer* than the classical geometric travel time.
;
; DELAY_ARG1, DELAY_ARG2 - can be any user-desired variables
; T1 - same as T1 passed to LITMSOL2
; X1,Y1,Z1 - same as passed to LITMSOL2
; T2 - trial T2 interaction time in TDB Julian days
; X2,Y2,Z2 - trial T2 interaction position, in POSUNITS
; POSUNITS, TBASE - same as passed to LITMSOL2
; ... additional keywords - passed via DELAY_FUNCTARGS
;
; 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)
;
; FUNC2 - a scalar string, is the name of subroutine to call which
; must compute position of body at time T2.
;
; INFO2, RAW2 - arguments to the FUNC2 interpolation function. At
; the very minimum, the INFO2 variable must be a
; structure of the form,
; INFO2 = {C: (speed of light in m/s), $
; AU: (1 AU in light-seconds), $
; ... other fields ... }
; The AU field is only required if POSUNITS EQ 'AU'.
;
; OUTPUTS:
;
; T2 - upon output, epoch of interaction at the second solar system
; body, in Julian days, in the TDB timescale.
;
; KEYWORD PARAMETERS:
;
; DELAY_FUNCTION - user function to compute extra delay factors
; based on the photon trajectory.
;
; DELAY_ARG1,DELAY_ARG2 - arguments to the DELAY_FUNCTION. These
; variables are not touched by LITMSOL2, but merely passed
; directly to DELAY_FUNCTION.
;
; DELAY_FUNCTARGS - a single structure containing additional keyword
; arguments passed to DELAY_FUNCTION using the _EXTRA method.
;
; 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.
;
; FUNCTARGS - a single structure containing additional keyword
; arguments passed to FUNC2 using the _EXTRA method.
;
; FUNCTSAVE - a named variable which will contain the results of
; the SAVE keyword when calling FUNC2 upon return.
;
; LIGHT_TIME - upon return, LIGHT_TIME is an array containing the
; estimated light time for each requested time.
;
; MAXITER - maximum number of solution iterations to be taken.
; Default: 5
;
; METHOD - solution method used, one of 'CONSTANT' or 'CORRECTOR'
; The 'CONSTANT' method uses simple iteration. The
; 'CORRECTOR' method uses a linear corrector to accelerate
; convergence by accounting for the line of sight velocity,
; but requires VX1, VY1, VZ1 to be passed.
; Default: 'CONSTANT'
;
; NITER - upon return, contains the actual number of iterations used.
;
; POSUNITS - the units for positions, one of 'CM', 'KM', 'LT-S' or
; 'AU'.
; Default: 'CM'
;
; RECEIVER - if set, then the epoch T1 is a reception of a photon.
; Otherwise T1 is the epoch of transmission of a photon.
;
; TGUESS - a vector of the same size as T1, containing an initial
; estimate of T2.
; Default: LITMSOL2 uses its own estimate based on T1.
;
; TOLERANCE - the solution tolerance, expressed in POSUNITS.
; Default: 1000 CM
;
; VX1, VY1, VZ1 - upon input, the body velocity at time T1, in
; VELUNITS units. This information is required only
; if the CORRECTOR method is used.
;
; VELUNITS - the units for velocities (and Shapiro derivative).
; Default: POSUNITS+'/S'
;
; X2, Y2, Z2, VX2, VY2, VZ2 - upon return, the body position and
; velocity at time T2, in units of POSUNITS and VELUNITS.
;
; EXAMPLE:
;
; SEE ALSO:
;
; JPLEPHREAD, JPLEPHINTERP, SHAPDEL
;
;
; MODIFICATION HISTORY:
; Major modifications, based on LITMSOL, 2009-01-05, CM
; Documented, 2009-05-12, CM
;
; $Id: litmsol2.pro,v 1.5 2010/05/04 21:01:52 craigm Exp $
;
;-
; Copyright (C) 2002, 2007, 2009, 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 litmsol2, t1, x1, y1, z1, t2, func2, info2, raw2, $
functargs=args2, functsave=save2, $
tolerance=tol0, posunits=posunits0, $
receiver=receiver, maxiter=maxiter0, light_time=ltm, $
tguess=tguess, error=diff, niter=i, $
vx1=vx1, vy1=vy1, vz1=vz1, $
x2=x2, y2=y2, z2=z2, vx2=vx2, vy2=vy2, vz2=vz2, $
method=method0, $
delay_function=delfunc, delay_arg1=delinfo, delay_arg2=delraw, $
delay_functargs=delargs
;; Default position and velocity units
if n_elements(posunits0) EQ 0 then begin
posunits = 'CM'
endif else begin
posunits = strtrim(posunits0(0),1)
endelse
velunits=posunits+'/S'
npts = n_elements(x1)
if (n_elements(y1) NE npts) OR (n_elements(z1) NE npts) then begin
message, 'ERROR: number of points in X1, Y1, Z1 do not match'
endif
usevel = 0
if n_elements(method0) EQ 0 then begin
method = 'CONSTANT'
endif else begin
method = strupcase(strtrim(method0(0),2))
endelse
if method EQ 'CORRECTOR' then begin
if (n_elements(vy1) NE npts) OR (n_elements(vz1) NE npts) then begin
message, 'ERROR: when METHOD="CORRECTOR" the number of velocity and position points must match'
endif
usevel = 1
endif
;; Default tolerances
if n_elements(tol0) EQ 0 then begin
tol = 1000d ;; 10 m tolerance
posunits = 'CM'
endif else begin
tol = tol0(0)
endelse
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
;; 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
endif else begin
t2 = t1
endelse
ltm = t2-t1
;; ==================== BEGIN ITERATION
ct = 1L
i = 0L
while (ct GT 0) AND (i LT maxiter) do begin
if arg_present(save2) OR n_elements(save2) GT 0 then begin
call_procedure, func2, info2, raw2, t1+ltm, $
x2, y2, z2, vx2, vy2, vz2, $
velocity=usevel, posunits=posunits, velunits=velunits, $
save=save2, _EXTRA=args2
endif else begin
call_procedure, func2, info2, raw2, t1+ltm, $
x2, y2, z2, vx2, vy2, vz2, $
velocity=usevel, posunits=posunits, velunits=velunits, $
_EXTRA=args2
endelse
;; Any systematic delays can be factored in as well (for example
;; the Shapiro delay). Must be computed in seconds
delay = 0
if n_elements(delfunc) GT 0 then begin
delay = call_function(delfunc, delinfo, delraw, $
t1, x1, y1, z1, t1+ltm, x2, y2, z2, $
posunits=posunits, tbase=tbase, $
_EXTRA=delargs)
delay = delay/86400d ;; Convert to days
endif
;; Compute distance from T1 to T2 in physical and light-second units
x21 = (x2-x1) & y21 = (y2-y1) & z21 = (z2-z1)
r21 = sqrt(x21^2 + y21^2 + z21^2)
p21 = 0
if method EQ 'CORRECTOR' then begin
;; Linear corrector factor
p21 = ((vx2-vx1)*x21 + (vy2-vy1)*y21 + (vz2-vz1)*z21)/r21
endif
if keyword_set(receiver) then begin
cfac = -clight
delay = -delay
endif else begin
cfac = clight
endelse
;; Linear corrector including line-of-sight velocity term
dtcor = (-ltm + r21/cfac/86400d + delay) / (1 - p21/cfac)
ltm = ltm + dtcor
diff = abs(dtcor)*clight*86400d
wh = where(diff GT tol, ct)
;; Prepare for next iteration
i = i + 1
endwhile
t2 = t1 + ltm
return
end