C+
C NAME:
C	StereoBOrbit
C PURPOSE:
C	Calculate Stereo B (position and velocity vectors)
C CATEGORY:
C	/gen/for/lib/ephem
C CALLING SEQUENCE:
	subroutine StereoBOrbit(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 stereoborbit.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.575995203295605d-01,2.717512429419864d-02,2.009817559818071d+02,2453977.546381154098d0,
     &		 1.265504225590081d+02,2.062856692554812d-01,
     &		 1.035787408847744d+00,4.699285093170608d-02,1.283202740400959d+02,2454109.145052516367d0,
     &		 3.392409078493881d+02,2.992287472926523d-01,
     &		 1.039104387456201d+00,4.372120282913373d-02,1.396338867376834d+02,2454119.296489763539d0,
     &		 3.385303364716559d+02,3.079140638493386d-01,
     &		 1.042293706168792d+00,4.182627806223717d-02,1.467482264878249d+02,2454124.112501636613d0,
     &		 3.365510702179411d+02,2.942094890131228d-01,
     &		 1.042810825735152d+00,4.143096626093998d-02,1.467863556683245d+02,2454512.857852065004d0,
     &		 3.364601578200525d+02,2.943195765731059d-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