C+
C NAME:
C	PlanetOrbit
C PURPOSE:
C	Calculate positions and velocities for planets
C	as a simple Kepler orbit
C CATEGORY:
C	Celestial mechanics
C CALLING SEQUENCE:
	subroutine PlanetOrbit(iPlanet,iYr,Doy,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	iYr	integer		time (year and day of year) where ..
C	Doy	real		.. position should be calculated
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	Julian, KeplerOrbit, PRECESSION_APPROX
C RESTRICTIONS:
C	The orbit is strictly a two-body solution.
C	No perturbations of any kind are taken into account.
C PROCEDURE:
C	The orbital parameters are hardwired.
C	The calculation is done by KeplerOrbit.
C
C	Average orbital elements and centennial rates of change
C	(from Explanatory Suppl. to the Astronomical Almanac,
C	Univ. Science Books, 1992, p. 316). Orbital elements
C	 semi-major axis (AU),
C	 eccentricity,
C	 inclination of orbital plane to ecliptic plane (deg).
C	 ecliptic longitude of ascending node (deg),
C	 longitude of perihelion (deg),
C	 mean longitude at epoch J2000 (deg)
C MODIFICATION HISTORY:
C	May 1997, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C-
	    integer	iPlanet
	    integer	iYr
	    real	Doy
	    real	RR(3)
	    real	VV(5)

	parameter	(nPlanets=9)

	! Average orbital elements at Epoch J2000

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

	! Centennial rates for orbit elements (angles in arcsec/century)
	! Presumably this converts the J200 elements to elements at the
	! current epoch iYr, Doy

	integer	nRevs(nPlanets)	/415,162,99,53,8,3,1,0,0/

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

	! 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
	double precision	dJEpoch
	double precision	JEpoch
	double precision	JEpoch0		/2000.0d0/

	real			OrbitElement(6)

	call Julian(20,iYr,Doy,JD,JEpoch)
	!JD = JD-0.5d0

	dJEpoch = (JEpoch-JEpoch0)/100.0d0

	! Orbital elements relative to current epoch iYr,Doy

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

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

	! Time to current epoch is 0.0d0

	JD = 0.0d0
	call KeplerOrbit(JD,1.0,sngl(1.0d0/InvMass(iPlanet)),OrbitElement,RR,VV)

	return
	end