C+
C NAME:
C	UlyssesOrbit
C PURPOSE:
C	Calculate Ulysses orbit (position and velocity vectors)
C CATEGORY:
C	Celestial mechanics
C CALLING SEQUENCE:
	subroutine UlyssesOrbit(iYr,Doy,RR,VV)
C INPUTS:
C	iYr		integer		year
C	Doy		real		day of year
C OUTPUTS:
C	RR(3)		real		position vector: ecliptic longitude and
C					latitude (deg), radial distance (AU)
C	VV(5)		real		velocity vector: ecliptic longitude and
C					latitude (deg), magnitude, radial and
C					tangential velocity (AU/day)
C					(will be zero if the Ulysses orbital
C					data base is used; see PROCEDURE).
C CALLS:
C	Julian, KeplerOrbit, iFilePath, iFltArr
C INCLUDE:
	include		'sun.h'
	include		'filparts.h'
	include		'dirspec.h'
C PROCEDURE:
C >	If the data base with Ulysses orbital coordinates exists (file
C	$dat:UlyssesPos.txt) the position vector is obtained by
C	linear interpolation on the data base.
C >	If the data base does not exist, approximate orbital parameters
C	are used. These parameters were determined by least-square fitting
C	Ulysses orbital data from the years 1993 to 1996 (first polar passage).
C MODIFICATION HISTORY:
C	JUN-1997, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C-
	    integer	iYr
	    real	Doy
	    real	RR(3)
	    real	VV(5)

	character	cSay*12	/'UlyssesOrbit'/
	! Copied from IDL routine ulyssesorbit.pro: semi-major axis (AU),
	! eccentricity, argument of periheion (deg), time of perihelion
	! passage (JD), longitude of ascending node (deg), and orbital
	! inclination (deg).
	! The four sets are for the in-ecliptic phase, and for the three
	! perihelion passages
	! All are J2000 elements.

	double precision	OrbitElement(6,4)
     &				/9.030365512305314d+00,8.897464400966297d-01,7.973512530451710d+00,2448177.155379012693d0,
     &				 1.336294632570382d+01,1.997907686222203d+00,

     &				 3.373607542718132d+00,6.031503055745338d-01,3.588812465846812d+02,2449788.991885926574d0,
     &				 3.381812752049081d+02,7.913419307720035d+01,

     &				 3.373734110021353d+00,6.031132414101130d-01,3.589804732058623d+02,2452053.178729594219d0,
     &				 3.381356025082639d+02,7.912197709000114d+01,

     &				 3.402423448713773d+00,5.905943668771040d-01,3.592550173721452d+02,2454330.786102185957d0,
     &				 3.380855550555008d+02,7.864681204100637d+01/

	save			OrbitElement

	real			Element(6)

	double precision	JD
	double precision	JEpoch

	call Julian(0,iYr,Doy,JD,JEpoch)

	if	(JD .le. 2448180.5d0) then
	    call Say(cSay,'E','time','invalid (1)')
	else if (JD .le. 2448620.5d0) then
	    i = 1		! 1990/10/16 - 1991/12/30: ecliptic phase
	else if (JD .le. 2448700.5d0) then
	    call Say(cSay,'E','time','invalid (2)')
	else if (JD .le. 2450921.5d0) then
	    i = 2		! 1992/03/19 - 1998/04/18: first orbit
	else if (JD .le. 2453187.5d0) then
	    i = 3		! 1998/04/18 - 2004/07/01: second orbit
	else if (JD .le. 2454831.5d0) then
	    i = 4		! 2004/07/01 - 2008/12/31: third orbit
	else
	    call Say(cSay,'E','time','invalid (3)')
	end if

	Element(1) = sngl(OrbitElement(1,i))		! Semi-major axis
	Element(2) = sngl(OrbitElement(2,i))		! Eccentricity
	Element(3) = sngl(OrbitElement(6,i))		! Orbital inclination
	Element(4) = sngl(OrbitElement(5,i))		! Longitude of ascending node

	! Longitude of perihelion
	! = Longitude ascending node + argument of perihelion

	Element(5) = OrbitElement(5,i)+OrbitElement(3,i)

	! Mean longitude at JD0 (perihelion passage)
	! = longitude of perihelion

	Element(6) = Element(5)

	call KeplerOrbit(JD-OrbitElement(4,i),1.0,0.0,Element,RR,VV)

	return
	end