;+ ; 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. ; ; 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' ; ; 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 ; ; NO_SHAPIRO - if set, then the Shapiro delay will not be accounted ; for. If NO_SHAPIRO is set, then INFOSUN and RAWSUN ; do not need to be supplied. ; ; ; 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 ; ; $Id: litmsol.pro,v 1.4 2002/09/25 21:19:45 craigm Exp $ ; ;- ; Copyright (C) 2002, 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, $ xoffset=xoff, yoffset=yoff, zoffset=zoff if n_elements(posunits0) EQ 0 then begin posunits = 'CM' endif else begin posunits = strtrim(posunits0(0),1) endelse 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(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 if NOT keyword_set(noshap) then begin jplephinterp, info, raw, t1, xs, ys, zs, /sun, $ posunits=posunits, tbase=tbase r1 = sqrt((x1-xs)^2 + (y1-ys)^2 + (z1-zs)^2) endif 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 ct = 1L i = 0L while ct GT 0 AND i LT maxiter do begin jplephinterp, info2, raw2, t2, x2, y2, z2, objectname=obj2, $ tbase=tbase, posunits=posunits dr = sqrt((x1-(x2+xoff))^2 + (y1-(y2+yoff))^2 + (z1-(z2+zoff))^2) dt = dr / clight ;; Add Shapiro delay if NOT keyword_set(noshap) then begin dt0 = 2*info.msol r2 = sqrt((x2-xs)^2 + (y2-ys)^2 + (z2-zs)^2) dt = dt + dt0*alog((r1+r2+dr+dt0*clight)/(r1+r2-dr+dt0*clight)) endif if keyword_set(receiver) then dt = -dt t2 = t1 + dt/86400d diff = abs(dtold-dt)*clight wh = where(diff GT tol, ct) dtold = dt i = i + 1 endwhile return end