C+
C NAME:
C	StereoAOrbit
C PURPOSE:
C	Calculate Stereo orbit (position and velocity vectors)
C CATEGORY:
C	/gen/for/lib/ephem
C CALLING SEQUENCE:
	subroutine StereoAOrbit(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, Say
C PROCEDURE:
C	Times have to be later than 2007/01/01 00 UT
C MODIFICATION HISTORY:
C	JUL-1998, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C-
	    integer	iYr
	    real	Doy
	    real	RR(3)
	    real	VV(5)

	! Copied from IDL routine stereoaorbit.pro:
	! semi-major axis (AU), eccentricity, argument of periheion (deg),
	! time of perihelion passage (JD), longitude of ascending node (deg),
	! and orbital inclination (deg).
	! All are J2000 elements.

	double precision	OrbitElement(6,5)
     &		/9.713114668842425d-01,9.509555173766652d-03,3.355229542907300d+02,2454210.608155566733d0,
     &		 2.380953113698972d+02,1.306235954146905d-01,
     &		 9.656833090947861d-01,6.279052167650651d-03,5.374894575771112d+01,2454265.510775187053d0,
     &		 2.169765459780858d+02,1.213139724936570d-01,
     &		 9.637020798975118d-01,6.226246274394636d-03,7.828787777683904d+01,2454286.911843170412d0,
     &		 2.149185709343068d+02,1.234975641202794d-01,
     &		 9.622945125484443d-01,5.917966070544746d-03,9.239562034525594d+01,2454299.754766134545d0,
     &		 2.143572837656983d+02,1.251469150146348d-01,
     &		 9.619379546429810d-01,5.655394382946006d-03,9.225614784090398d+01,2454299.536585085094d0,
     &		 2.142571411361414d+02,1.253836830038676d-01/

	save			OrbitElement

	real			Element(6)

	double precision	JD
	double precision	JEpoch

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

	if	(jd .le. 2454101.5d0) then
	    call Say('StereoAOrbit','E','time','invalid (1)')
	else if (jd .le. 2454132.5d0) then
	    i = 1		! 2007/01/01 - 2007/02/01
	else if (jd .le. 2454160.5d0) then
	    i = 2		! 2007/02/01 - 2007/03/01
	else if (jd .le. 2454191.5d0) then
	    i = 3		! 2007/03/01 - 2007/04/01
	else if (jd .le. 2454374.5d0) then
	    i = 4		! 2007/04/01 - 2007/10/01
	else if (jd .le. 2455958.5d0) then
	    i = 5		! 2007/10/01 - 2012/01/01
	else
	    call Say('StereoAOrbit','E','time','invalid (2)')
	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