;+
; NAME:
;	UlyssesOrbit
; PURPOSE:
;	Calculate position of Ulysses
; CATEGORY:
;	smei/gen/idl/ephem
; CALLING SEQUENCE:
	FUNCTION UlyssesOrbit, 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 results of this routine were compared with J2000 ecliptic
;	coordinates from http://cohoweb.gsfc.nasa.gov/helios/heli.html.
;	For the most part the agreement is smaller than 0.1 degrees,
;	but during the last two perihelion passages the discrepancy
;	in longitude and latitude rises to 1-2 degrees.
;	Adding more orbital elements where the discrepancies are
;	largest does not fix the problem
; PROCEDURE:
;	Orbital elements obtained herein are from JPL's HORIZONS System
;	(http://ssd.jpl.nasa.gov/?horizons)
;	(these are J2000 orbital elements)
;	
;	Earth-Jupiter segment (Oct 1990 - Dec 1991)
;
;	2448400.500000000 = A.D. 1991-May-24 00:00:00.0000 (CT)
;	 EC= 8.897464400966297E-01 QR= 9.956299449602790E-01 IN= 1.997907686222203E+00
;	 OM= 1.336294632570382E+01 W = 7.973512530451710E+00 Tp=  2448177.155379012693
;	 N = 3.632002010001063E-02 MA= 8.111881123489832E+00 TA= 1.165108914377567E+02
;	 A = 9.030365512305314E+00 AD= 1.706510107965035E+01 PR= 9.911888787745871E+03
;
;	First perihelion (1995/03/12):
;
;	2449788.500000000 = A.D. 1995-Mar-12 00:00:00.0000 (CT)
;	 EC= 6.031503055745338E-01 QR= 1.338815122439140E+00 IN= 7.913419307720035E+01
;	 OM= 3.381812752049081E+02 W = 3.588812465846812E+02 Tp=  2449788.991885926574
;	 N = 1.590604604763074E-01 MA= 3.599217603979999E+02 TA= 3.596037474470892E+02
;	 A = 3.373607542718132E+00 AD= 5.408399962997124E+00 PR= 2.263290316914575E+03
;
;	Second perihelion (2001/05/23):
;
;	2452053.500000000 = A.D. 2001-May-24 00:00:00.0000 (CT)
;	 EC= 6.031132414101130E-01 QR= 1.338990395270509E+00 IN= 7.912197709000114E+01
;	 OM= 3.381356025082639E+02 W = 3.589804732058623E+02 Tp=  2452053.178729594219
;	 N = 1.590515097125277E-01 MA= 5.109854306171003E-02 TA= 2.587555018718575E-01
;	 A = 3.373734110021353E+00 AD= 5.408477824772197E+00 PR= 2.263417685570354E+03
;
;	Third perihelion (2007/08/18):
;
;	2454330.500000000 = A.D. 2007-Aug-18 00:00:00.0000 (CT)
;	 EC= 5.905943668771040E-01 QR= 1.392971326172849E+00 IN= 7.864681204100637E+01
;	 OM= 3.380855550555008E+02 W = 3.592550173721452E+02 Tp=  2454330.786102185957
;	 N = 1.570440655511920E-01 MA= 3.599550693495362E+02 TA= 3.597836828567464E+02
;	 A = 3.402423448713773E+00 AD= 5.411875571254696E+00 PR= 2.292350231360063E+03
;
;    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:
;	SEP-1999, Paul Hick (UCSD/CASS)
;	JUL-2007, George Megally (UCSD/CASS; gmegally@ucsd.edu)
;	JUL-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Replaced all orbit elements by elements at perihelion
;		from JPL Horizon.
;-

dpm = ToDegrees(degrees=degrees)

jd = TimeGet(T, /njd)

p = replicate(BadValue(0.0),3,n_elements(jd))

; In-ecliptic phase

tt = TimeGet(TimeSet(yr=[1990,1991],month=['oct','dec'],day=[16,30]), /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.030365512305314d+00,8.897464400966297d-01,7.973512530451710d+00/dpm,2448177.155379012693d0-2451545.0d0], $
		plane		= [1.336294632570382d+01,1.997907686222203d+00]/dpm)

; First orbit

tt = TimeGet(TimeSet(yr=[1992,1998],month=['mar','apr'],day=[19,18]), /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	= [3.373607542718132d+00,6.031503055745338d-01,3.588812465846812d+02/dpm,2449788.991885926574d0-2451545.0d0], $
		plane		= [3.381812752049081d+02,7.913419307720035d+01]/dpm)

; Second orbit

tt = TimeGet(TimeSet(yr=[1998,2004],month=['apr','jul'],day=[18,01]), /njd)

i = where( tt[0] LT jd AND jd LE tt[1] )

IF i[0] NE -1 THEN	$
	p[*,i] = KeplerOrbit(T[i], degrees=degrees	 , $
		elements	= [3.373734110021353d+00,6.031132414101130d-01,3.589804732058623d+02/dpm,2452053.178729594219d0-2451545.0d0], $
		plane		= [3.381356025082639d+02,7.912197709000114d+01]/dpm)

tt = TimeGet(TimeSet(yr=[2004,2008],month=['jul','dec'],day=[01,31]), /njd)

; Third orbit

i = where( tt[0] LT jd AND jd LE tt[1] )

IF i[0] NE -1 THEN	$
	p[*,i] = KeplerOrbit(T[i], degrees=degrees	 , $
		elements	= [3.402423448713773d+00,5.905943668771040d-01,3.592550173721452d+02/dpm,2454330.786102185957d0-2451545.0d0], $
		plane		= [3.380855550555008d+02,7.864681204100637d+01]/dpm)

RETURN, p  &  END