C+
C NAME:
C	comet67PCGorbit
C PURPOSE:
C	Calculate comet 67PCG (Rosetta) orbit (position and velocity vectors)
C CATEGORY:
C	Celestial mechanics
C CALLING SEQUENCE:
	subroutine comet67PCGorbit(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, iFilePath, iFltArr
C INCLUDE:
	include		'sun.h'
	include		'filparts.h'
	include		'dirspec.h'
C PROCEDURE:
C >	If the data base with comet 67PCG orbital coordinates exists (file
C	$dat:67PCGPos.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)
C	JULY-2014, Bernard V. Jackson (UCSD/CASS; bvjackson@ucsd.edu)
C-
	    integer	iYr
	    real	Doy
	    real	RR(3)
	    real	VV(5)

	character	cSay*12	/'UlyssesOrbit'/
	! Copied from IDL routine ulyssesorbit.f: semi-major axis (AU),
	! eccentricity, argument of periheion (deg), time of perihelion
	! passage (JD), longitude of ascending node (deg), and orbital
	! inclination (deg).
	! The four sets are for the in-ecliptic phase, and for the three
	! perihelion passages
	! All are J2000 elements.

	double precision	OrbitElement(6,1)
     &				/3.463242237470248d+00,6.410532328350245d-01,1.2773219225972d+01,2457247.531305577176d0,
     &				 5.015180617652601d+01,7.040848762748019d+00/

	save			OrbitElement

	real			Element(6)

	double precision	JD
	double precision	JEpoch

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

	if	(JD .le. 2448180.5d0) then
	    call Say(cSay,'E','time','invalid (1)')
	else 
	    i = 1		! 2014e
	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