C+
C NAME:
C	StereoOrbit
C PURPOSE:
C	Calculate the Stereo A and B orbit (position locations)
C CATEGORY:
C	Celestial mechanics
C CALLING SEQUENCE:
	subroutine StereoOrbit(nS,iYr,Doy,RR)
C INPUTS:
C	nS		1 - Stereo A
C			2 - Stereo B
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 CALLS:
C	Julian, KeplerOrbit
C INCLUDE:
C
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), JUL-2008 (UCSD/CASS; bjackson@ucsd.edu)
C-
	    integer	iYr
	    real	Doy
	    real	RR(3)
	    real	VV(5)

	double precision	JD
	double precision	JEpoch

	save			iPos

	real			OrbitElement1(6) /	! STEREO A
     &					0.962296479,	! Semi-major axis (AU)
     &					0.005911967,	! Eccentricity
     &					0.125191913,	! Inclination (deg)
     &					214.3493931,	! Lng ascending node (deg)
     &					306.88329,	! Lng of perihelion (deg)
     &                                  306.88329/      ! Mean lng on JD0
                                                        ! 92.53389979=arg perihelion

	real			OrbitElement2(6) /	! STEREO B
     &					1.042293101,	! Semi-major axis (AU)
     &					0.041831403,	! Eccentricity
     &					0.294262911,	! Inclination (deg)
     &					336.5674914,	! Lng ascending node (deg)
     &					483.29657,	! Lng of perihelion (deg)
     &                                  483.29657/      ! Mean lng on JD0
                                                        ! 146.7290772=arg perihelion

	double precision	JD01    /2454300.0d0/	! Following J2000
	double precision	JD02    /2454124.0d0/	! Following J2000	

	save			JD01, JD02, OrbitElement1, OrbitElement2

	if(iYr.lt.2007) then
	
	  print *, 'STEREO A or B orbit data not valid'
	  return
	  
	end if
	
	if(iYr.eq.2007.and.Doy.lt.90.0) print *, 'STEREO A or B orbit data inaccurate.'

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

	if(nS.lt.1.or.nS.gt.2) then
	
	  print *, 'Enter either STEREO A or STEREO B (1 or 2).'
	  
	end if

	if(nS.eq.1) then
	  call KeplerOrbit(JD-JD01,1.,0.,OrbitElement1,RR,VV)	! STEREO A
	else
	  call KeplerOrbit(JD-JD02,1.,0.,OrbitElement2,RR,VV)	! STEREO B	  
	end if

	return
	end