C+
C NAME:
C	Time2PlanetOrbit
C PURPOSE:
C	Calculate positions and velocities for planets as a simple Kepler orbit
C CATEGORY:
C	Celestial mechanics
C CALLING SEQUENCE:
	subroutine Time2PlanetOrbit(iPlanet,tt,rr,vv)
C INPUTS:
C	iPlanet		integer		Planet identifier
C					1: Mercury
C					2: Venus
C					3: Earth/Moon
C					4: Mars
C					5: Jupiter
C					6: Saturn
C					7: Uranus
C					8: Neptune
C					9: Pluto
C	tt(2)		integer		time
C OUTPUTS:
C	rr(3)		real		position vector relative to m1:
C					1: ecliptic longitude (deg),
C					2: ecliptic latitude (deg),
C					3: radial distance (AU)
C	vv(5)		real		velocity vector relative to m1:
C					1: ecliptic longitude (deg),
C					2: ecliptic latitude (deg),
C					3: velocity (AU/day),
C					4: tangential velocity (AU/day),
C					5: radial velocity (AU/day)
C CALLS:
C	KeplerOrbit8, PRECESSION_APPROX
C RESTRICTIONS:
C	The orbit is strictly a two-body solution. No perturbations of
C	any kind are taken into account.
C PROCEDURE:
C	The orbital parameters are hardwired. The calculation is done
C	by a call to KeplerOrbit.
C MODIFICATION HISTORY:
C	May 1997, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C-
	    integer		iPlanet
	    integer		tt(2)
	    double precision	rr(3)
	    double precision	vv(5)

	parameter	(nPlanets=9)

	!-------
	! Average orbital elements and centennial rates of change (from Explanatory
	! Suppl. to the Astronomical Almanac, Univ. Science Books, 1992, p. 316).
	!	semi-major axis (AU),
	!	eccentricity,
	!	inclination of orbital plane to ecliptic plane (deg).
	!	ecliptic longitude of ascending node (deg),
	!	argument of perihelion (deg),
	!	mean longitude at epoch J2000 (deg),

	double precision	EE(6,nPlanets)	/
     &		 0.38709893d0, 0.20563069d0, 7.00487d0,  48.33167d0,  77.45645d0, 252.25084d0,
     &		 0.72333199d0, 0.00677323d0, 3.39471d0,  76.68069d0, 131.53298d0, 181.97973d0,
     &		 1.00000011d0, 0.01671022d0, 0.00005d0, -11.26064d0, 102.94719d0, 100.46435d0,
     &		 1.52366231d0, 0.09341233d0, 1.85061d0,  49.57854d0, 336.04084d0, 355.45332d0,
     &		 5.20336301d0, 0.04839266d0, 1.30530d0, 100.55615d0,  14.75385d0,  34.40438d0,
     &		 9.53707032d0, 0.05415060d0, 2.48446d0, 113.71504d0,  92.43194d0,  49.94432d0,
     &		19.19126393d0, 0.04716771d0, 0.76986d0,  74.22988d0, 170.96424d0, 313.23218d0,
     &		30.06896348d0, 0.00858587d0, 1.76917d0, 131.72169d0,  44.97135d0, 304.88003d0,
     &		39.48168677d0, 0.24880766d0,17.14175d0, 110.30347d0, 224.06676d0, 238.92881d0/

	!-------
	! Centennial rates for orbit elements (angles in arcsec/century)
        ! Presumably these convert from J2000 elements to elements at the current epoch.

	double precision	nRevs(nPlanets)	/415.0d0,162.0d0,99.0d0,53.0d0,8.0d0,3.0d0,1.0d0,0.0d0,0.0d0/

	double precision	dEE(6,nPlanets)	/
     &		 0.00000066d0, 0.00002527d0,-23.52d0,  -446.30d0,  573.57d0, 261628.29d0,
     &		 0.00000092d0,-0.00004938d0, -2.86d0,  -996.89d0, -108.80d0, 712136.06d0,
     &		-0.00000005d0,-0.00003804d0,-46.94d0,-18228.25d0, 1198.28d0,1293740.63d0,
     &		-0.00007221d0, 0.00011902d0,-25.47d0, -1020.19d0, 1560.78d0, 217103.78d0,
     &		 0.00060737d0,-0.00012880d0, -4.15d0,  1217.17d0,  839.93d0, 557078.35d0,
     &		-0.00301530d0,-0.00036762d0,  6.11d0, -1591.05d0,-1948.89d0, 513052.95d0,
     &		 0.00152025d0,-0.00019150d0, -2.09d0,  1681.40d0, 1312.56d0, 246547.79d0,
     &		-0.00125196d0, 0.00002514d0, -3.64d0,  -151.25d0, -844.43d0, 786449.21d0,
     &		-0.00076912d0, 0.00006465d0, 11.07d0,   -37.33d0, -132.25d0, 522747.90d0/

	!-------
	! Inverse planetary masses in units of the solar mass, Msun/Mplanet.
	! (from Astrophysical Data: Planets and Stars, K.R. Lang)

	double precision	InvMass(nPlanets)	/
     &		  60236.0000d2,
     &		 408525.1000d0,
     &		 328900.5550d0,
     &		3098710.0000d0,
     &		   1047.3492d0,
     &		   3497.9100d0,
     &		  22902.9400d0,
     &		  19412.2500d0,
     &		     13.0000d7/

	double precision	jd
	integer			t0(2)

	double precision	OrbitElement(6)

	call Time2JEpoch(2,tt,jd)	! Julian epoch - 2000
	jd = jd/100.0d0

        ! Set up elements at current epoch jd

	do i=1,2			! Semi-major axis, eccentricity
	    OrbitElement(i) = EE(i,iPlanet)+dEE(i,iPlanet)*jd
	end do

	do i=3,6			! Inclination, ascending node, arg. perihelion
	    OrbitElement(i) = EE(i,iPlanet)+dEE(i,iPlanet)/3600.0d0*jd
	end do
					! Mean longitude at (iYr,Doy)
	OrbitElement(6) = OrbitElement(6)+360.0d0*nRevs(iPlanet)*jd

        ! Time to current epoch is zero

        t0(1) = 0
        t0(2) = 0
	call Time2KeplerOrbit(t0,1.0d0,1.0d0/InvMass(iPlanet),OrbitElement,RR,VV)

	return
	end