;+ ; NAME: ; StereoBOrbit ; PURPOSE: ; Calculate position of Stereo B (Behind) ; CATEGORY: ; smei/gen/idl/ephem ; CALLING SEQUENCE: FUNCTION StereoBOrbit, T, degrees=degrees ; INPUTS: ; T array[n]; type: time structure ; times at which s/c positions are needed ; OPTIONAL INPUT PARAMETERS: ; /degrees if set then angles are output in degrees (default: radians) ; OUTPUTS: ; Result array[3,*]; type: float ; J2000 ecliptic longitude (deg/rad), latitude ; (deg/rad) and distance (AU). ; The angular units depend on the setting of /degrees ; INCLUDE: @compile_opt.pro ; On error, return to caller ; CALLS: ; TimeGet, TimeSet, ToDegrees, KeplerOrbit ; RESTRICTIONS: ; The precision is probably about 0.1 degrees in ecliptic longitude, ; provided orbital elements are added as needed (once or twice year?) ; PROCECURE: ; From http://ssd.jpl.nasa.gov/horizons.cgi ; (J2000.0 Orbital Elements) ; Elements cover four periods between 2007/01/01 and 2007/09/30: ; ; 2007/01/01 -> 2007/01/31 ; ; 2454115.500000000 = A.D. 2007-Jan-15 00:00:00.0000 (CT) ; EC= 2.717512429419864E-02 QR= 9.315766343405397E-01 IN= 2.062856692554812E-01 ; OM= 1.265504225590081E+02 W = 2.009817559818071E+02 Tp= 2453977.546381154098 ; N = 1.051787957703235E+00 MA= 1.450979550234674E+02 TA= 1.468311244708070E+02 ; A = 9.575995203295605E-01 AD= 9.836224063185813E-01 PR= 3.422743123871886E+02 ; ; 2007/02/01 -> 2007/02/28 ; ; 2454147.500000000 = A.D. 2007-Feb-16 00:00:00.0000 (CT) ; EC= 4.699285093170608E-02 QR= 9.871128055468245E-01 IN= 2.992287472926523E-01 ; OM= 3.392409078493881E+02 W = 1.283202740400959E+02 Tp= 2454109.145052516367 ; N = 9.349709833378660E-01 MA= 3.586076296454760E+01 TA= 3.917087846701419E+01 ; A = 1.035787408847744E+00 AD= 1.084462012148664E+00 PR= 3.850386872058772E+02 ; ; 2007/03/01 -> 2007/03/31 ; ; 2454175.500000000 = A.D. 2007-Mar-16 00:00:00.0000 (CT) ; EC= 4.372120282913373E-02 QR= 9.936734937715849E-01 IN= 3.079140638493386E-01 ; OM= 3.385303364716559E+02 W = 1.396338867376834E+02 Tp= 2454119.296489763539 ; N = 9.304977044270033E-01 MA= 5.229723725577440E+01 TA= 5.639454403113316E+01 ; A = 1.039104387456201E+00 AD= 1.084535281140816E+00 PR= 3.868897239479882E+02 ; ; 2007/04/01 -> 2007/10/01 ; ; 2454253.500000000 = A.D. 2007-Jun-02 00:00:00.0000 (CT) ; EC= 4.182627806223717E-02 QR= 9.986984397920566E-01 IN= 2.942094890131228E-01 ; OM= 3.365510702179411E+02 W = 1.467482264878249E+02 Tp= 2454124.112501636613 ; N = 9.262301226788787E-01 MA= 1.198425984820785E+02 TA= 1.238910903027097E+02 ; A = 1.042293706168792E+00 AD= 1.085888972545528E+00 PR= 3.886723085174492E+02 ; ; 2007/10/01 -> now ; ; 2454466.500000000 = A.D. 2008-Jan-01 00:00:00.0000 (CT) ; EC= 4.143096626093998E-02 QR= 9.996061655975762E-01 IN= 2.943195765731059E-01 ; OM= 3.364601578200525E+02 W = 1.467863556683245E+02 Tp= 2454512.857852065004 ; N = 9.255412455738544E-01 MA= 3.170938958574695E+02 TA= 3.137364080607738E+02 ; A = 1.042810825735152E+00 AD= 1.086015485872728E+00 PR= 3.889615959543680E+02 ; ; JDCT Epoch Julian Date, Coordinate Time ; EC Eccentricity, e ; QR Periapsis distance, q (AU) ; IN Inclination w.r.t xy-plane, i (degrees) ; OM Longitude of Ascending Node, OMEGA, (degrees) ; W Argument of Perifocus, w (degrees) ; Tp Time of periapsis (Julian day number) ; N Mean motion, n (degrees/day) ; MA Mean anomaly, M (degrees) ; TA True anomaly, nu (degrees) ; A Semi-major axis, a (AU) ; AD Apoapsis distance (AU) ; PR Orbital period (day) ; MODIFICATION HISTORY: ; AUG-2007, George Megally (UCSD/CASS; gmegally@ucsd.edu) ; AUG-2007, John Clover (UCSD/CASS; jclover@ucsd.edu) ; JUL-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu) ; Split up into seperate routines for A and B. ; Updated orbital elements from Horizon (almost, but not quite the ; same as the original ones). ; JUN-2009, John Clover (UCSD/CASS; jclover@ucsd.edu) ; Fixed time structure to cover the whole of the end of the day between ; elements. ;- dpm = ToDegrees(degrees=degrees) jd = TimeGet(T, /njd) p = replicate(BadValue(0.0),3,n_elements(jd)) tt = TimeGet(TimeSet(yr=2007, month='jan',day=[1,31.9999]), /njd) i = where( tt[0] LE jd AND jd LE tt[1] ) IF i[0] NE -1 THEN $ p[*,i] = KeplerOrbit(T[i], degrees=degrees , $ elements = [9.575995203295605d-01,2.717512429419864d-02,2.009817559818071d+02/dpm,2453977.546381154098d0-2451545.0d0], $ plane = [1.265504225590081d+02,2.062856692554812d-01]/dpm) tt = TimeGet(TimeSet(yr=2007, month='feb',day=[1,28.9999]), /njd) i = where( tt[0] LE jd AND jd LE tt[1] ) IF i[0] NE -1 THEN $ p[*,i] = KeplerOrbit(T[i], degrees=degrees , $ elements = [1.035787408847744d+00,4.699285093170608d-02,1.283202740400959d+02/dpm,2454109.145052516367d0-2451545.0d0], $ plane = [3.392409078493881d+02,2.992287472926523d-01]/dpm) tt = TimeGet(TimeSet(yr=2007, month='mar',day=[1,31.9999]), /njd) i = where( tt[0] LE jd AND jd LE tt[1] ) IF i[0] NE -1 THEN $ p[*,i] = KeplerOrbit(T[i], degrees=degrees , $ elements = [1.039104387456201d+00,4.372120282913373d-02,1.396338867376834d+02/dpm,2454119.296489763539d0-2451545.0d0], $ plane = [3.385303364716559d+02,3.079140638493386d-01]/dpm) tt = TimeGet(TimeSet(yr=2007, month=['apr','oct'], day=1), /njd) i = where( tt[0] LE jd AND jd LE tt[1] ) IF i[0] NE -1 THEN $ p[*,i] = KeplerOrbit(T[i], degrees=degrees , $ elements = [1.042293706168792d+00,4.182627806223717d-02,1.467482264878249d+02/dpm,2454124.112501636613d0-2451545.0d0], $ plane = [3.365510702179411d+02,2.942094890131228d-01]/dpm) tt = TimeGet(TimeSet(yr=[2007,2014], month=['oct','jan'], day=1), /njd) i = where( tt[0] LE jd AND jd LE tt[1] ) IF i[0] NE -1 THEN $ p[*,i] = KeplerOrbit(T[i], degrees=degrees , $ elements = [1.042810825735152d+00,4.143096626093998d-02,1.467863556683245d+02/dpm,2454512.857852065004d0-2451545.0d0], $ plane = [3.364601578200525d+02,2.943195765731059d-01]/dpm) tt = TimeGet(TimeSet(yr=[2011,2012], month='jan', day=1), /njd) i = where( tt[0] LE jd AND jd LE tt[1] ) IF i[0] NE -1 THEN $ p[*,i] = KeplerOrbit(T[i], degrees=degrees , $ elements = [1.042950389141085d+00,4.141689777858453d-02,1.469301650396323d+02/dpm,2455680.117963262834d0-2451545.0d0], $ plane = [3.364380048102094d+02,2.939527214019730d-01]/dpm) tt = TimeGet(TimeSet(yr=[2012,2022], month='jan', day=1), /njd) i = where( tt[0] LE jd AND jd LE tt[1] ) IF i[0] NE -1 THEN $ p[*,i] = KeplerOrbit(T[i], degrees=degrees , $ elements = [1.042929841928075d+00,4.146134713825472d-02,1.469129271961573d+02/dpm,2456069.114342078567d0-2451545.0d0], $ plane = [3.364246242628766d+02,2.937938987140259d-01]/dpm) RETURN, p & END