C+
C NAME:
C       MessengerOrbit	
C PURPOSE:
C	Calculate Messenger orbit (position and velocity vectors)
C CATEGORY:
C	/gen/for/lib/ephem
C CALLING SEQUENCE:
	subroutine MessengerOrbit(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 Messenger 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,8)
     &          /1.001663657497907d+00,8.168951593474302d-02,2.563280663025649d+02,2453124.628101623151d0,
     &		 3.111249722587803d+02,6.845013980422952d+00,
     &           8.019447095476302d-01,2.658523304784657d-01,2.063396658772084d+00,2453456.712469597347d0,
     &           1.305447284592106d+02,2.664216629562539d+00,
     &           7.228282327538631d-01,2.463851525210938d-01,5.790967595762318d+01,2453985.644399113953d0,
     &           4.647221279089646d+01,8.212210576757776d+00,
     &           5.282648296522885d-01,4.182978788124355d-01,3.585284077838588d+02,2454202.959927953314d0,
     &           4.837960301496808d+01,7.403689052585330d+00,
     &           5.077205481037513d-01,3.828930538271369d-01,6.424844763538441d+00,2454488.927706247196d0,
     &           4.903726284801532d+01,6.920821620411199d+00,
     &           4.665958466483129d-01,3.520983387776521d-01,1.943142669334450d+01,2454754.833715905435d0,
     &           4.833131986457225d+01,7.003009660944397d+00,
     &           4.392148766241953d-01,3.109620288640121d-01,3.258748902171746d+01,2455109.302512918599d0,
     &           4.831869848958493d+01,7.008490679204085d+00,
     &           4.216103227199027d-01,2.700790583878330d-01,2.460386393404945d+01,2455637.410266774241d0,
     &           5.597165374661773d+01,7.889636530862079d+00/

	save			OrbitElement

	real			Element(6)

	double precision	JD
	double precision	JEpoch

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

	if	(jd .le. 2453220.5d0) then
	    call Say('MessengerOrbit','E','time','invalid (1)')
	else if (jd .le. 2453586.5d0) then
	    i = 1		! 2004/08/03 - 2005/08/04
        else if (jd .le. 2454033.5d0) then
            i = 2               ! 2005/08/04 - 2006/10/25
        else if (jd .le. 2454257.5d0) then
            i = 3               ! 2006/10/25 - 2007/06/06
        else if (jd .le. 2454480.5d0) then
            i = 4               ! 2007/06/06 - 2008/01/15
        else if (jd .le. 2454746.5d0) then
            i = 5               ! 2008/01/15 - 2008/10/07
        else if (jd .le. 2455104.5d0) then
            i = 6               ! 2008/10/07 - 2009/09/30
        else if (jd .le. 2455640.5d0) then
            i = 7               ! 2009/09/30 - 2011/03/20
        else if (jd .le. 2459658.5d0) then
            i = 8               ! 2011/03/20 - 2022/03/20
	else
	    call Say('MessengerOrbit','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