;+
; NAME:
;	StereoAOrbit
; PURPOSE:
;	Calculate position of Stereo A (Ahead)
; CATEGORY:
;	smei/gen/idl/ephem
; CALLING SEQUENCE:
	FUNCTION StereoAOrbit, 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= 9.509555173766652E-03 QR= 9.620747268989945E-01 IN= 1.306235954146905E-01
;	 OM= 2.380953113698972E+02 W = 3.355229542907300E+02 Tp=  2454210.608155566733
;	 N = 1.029594704279364E+00 MA= 2.620771466946642E+02 TA= 2.609996629128498E+02
;	 A = 9.713114668842425E-01 AD= 9.805482068694904E-01 PR= 3.496521480770163E+02
;
;	2007/02/01 -> 2007/02/28
;
;	2454147.500000000 = A.D. 2007-Feb-16 00:00:00.0000 (CT)
;	 EC= 6.279052167650651E-03 QR= 9.596197332195503E-01 IN= 1.213139724936570E-01
;	 OM= 2.169765459780858E+02 W = 5.374894575771112E+01 Tp=  2454265.510775187053
;	 N = 1.038608771871410E+00 MA= 2.374329737152352E+02 TA= 2.368291467561132E+02
;	 A = 9.656833090947861E-01 AD= 9.717468849700218E-01 PR= 3.466175231231067E+02
;
;	2007/03/01 -> 2007/03/31
;
;	2454175.500000000 = A.D. 2007-Mar-16 00:00:00.0000 (CT)
;	 EC= 6.226246274394636E-03 QR= 9.577018334129234E-01 IN= 1.234975641202794E-01
;	 OM= 2.149185709343068E+02 W = 7.828787777683904E+01 Tp=  2454286.911843170412
;	 N = 1.041813256894421E+00 MA= 2.439296648101152E+02 TA= 2.432909800026552E+02
;	 A = 9.637020798975118E-01 AD= 9.697023263821002E-01 PR= 3.455513717238894E+02
;
;	2007/04/01 -> 2007/10/01
;
;	2454253.500000000 = A.D. 2007-Jun-02 00:00:00.0000 (CT)
;	 EC= 5.917966070544746E-03 QR= 9.565996862733112E-01 IN= 1.251469150146348E-01
;	 OM= 2.143572837656983E+02 W = 9.239562034525594E+01 Tp=  2454299.754766134545
;	 N = 1.044099914070154E+00 MA= 3.117054026536591E+02 TA= 3.111966164343318E+02
;	 A = 9.622945125484443E-01 AD= 9.679893388235775E-01 PR= 3.447945882847869E+02
;
;	2007/10/01 -> now
;
;	2454466.500000000 = A.D. 2008-Jan-01 00:00:00.0000 (CT)
;	 EC= 5.655394382946006E-03 QR= 9.564978161375505E-01 IN= 1.253836830038676E-01
;	 OM= 2.142571411361414E+02 W = 9.225614784090398E+01 Tp=  2454299.536585085094
;	 N = 1.044680486714094E+00 MA= 1.744234215566499E+02 TA= 1.744859573944121E+02
;	 A = 9.619379546429810E-01 AD= 9.673780931484115E-01 PR= 3.446029715098184E+02
;
;	2011/01/21 -> 2012/01/21
;
;	2455583.500000000 = A.D. 2011-Jan-22 00:00:00.0000 (CT)
;	 EC= 5.649610774525156E-03 QR= 9.563712745374616E-01 IN= 1.256261457055498E-01
;	 OM= 2.140504362161804E+02 W = 9.348075068567380E+01 Tp=  2455678.680902264081
;	 N = 1.044896949233696E+00 MA= 2.605457655988884E+02 TA= 2.599079145507906E+02
;	 A = 9.618050989876956E-01 AD= 9.672389234379295E-01 PR= 3.445315830082727E+02
;
;	2012/01/21 -> now
;
;	2455948.500000000 = A.D. 2012-Jan-22 00:00:00.0000 (CT)
;	 EC= 5.705938568515530E-03 QR= 9.563266675804051E-01 IN= 1.257640058659871E-01
;	 OM= 2.139833070619798E+02 W = 9.406913424182217E+01 Tp=  2456023.716625969391
;	 N = 1.044881265762684E+00 MA= 2.814075566505216E+02 TA= 2.807657290779928E+02
;	 A = 9.618147233058820E-01 AD= 9.673027790313589E-01 PR= 3.445367543624463E+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=[01, 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.713114668842425d-01,9.509555173766652d-03,3.355229542907300d+02/dpm,2454210.608155566733d0-2451545.0d0],	$
		plane	 = [2.380953113698972d+02,1.306235954146905d-01]/dpm)

tt = TimeGet(TimeSet(yr=2007, month='feb',day=[01, 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	= [9.656833090947861d-01,6.279052167650651d-03,5.374894575771112d+01/dpm,2454265.510775187053d0-2451545.0d0],	$
		plane		= [2.169765459780858d+02,1.213139724936570d-01]/dpm)

tt = TimeGet(TimeSet(yr=2007, month='mar',day=[01, 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.637020798975118d-01,6.226246274394636d-03,7.828787777683904d+01/dpm,2454286.911843170412d0-2451545.0d0],	$
		plane		= [2.149185709343068d+02,1.234975641202794d-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	= [9.622945125484443d-01,5.917966070544746d-03,9.239562034525594d+01/dpm,2454299.754766134545d0-2451545.0d0], $
		plane		= [2.143572837656983d+02,1.251469150146348d-01]/dpm)

tt = TimeGet(TimeSet(yr=[2007,2012], 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	= [9.619379546429810d-01,5.655394382946006d-03,9.225614784090398d+01/dpm,2454299.536585085094d0-2451545.0d0], $
		plane		= [2.142571411361414d+02,1.253836830038676d-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	= [9.618050989876956d-01,5.649610774525156d-03,9.348075068567380d+01/dpm,2455678.680902264081d0-2451545.0d0], $
		plane		= [2.140504362161804d+02,1.256261457055498d-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	= [9.618147233058820d-01,5.705938568515530d-03,9.406913424182217d+01/dpm,2456023.716625969391d0-2451545.0d0], $
		plane		= [2.139833070619798d+02,1.257640058659871d-01]/dpm)
RETURN, p  &  END