C+
C NAME:
C	jpl_eph2
C PURPOSE:
C	Read the JPL planetary ephemeris and gives position and velocity
C	of body 'ntarg' with respect to body 'ncent'.
C CATEGORY:
C	gen/for/lib
C CALLING SEQUENCE:
	logical function jpl_eph2(et,ntarg_,ncent_,rrd,bKm,bBary)
C INPUTS:
C	et(2)	double precision	Ephemeris time specified as Julian day.
C			A. put entire epoch in et(1) and set et(2)=0 (equivalent to jpl_eph)
C			B. for maximum interpolatin accuracy, set et(1) to the most recent midnight
C			   at or before the interpolation time, and put the remaining fraction in et(2)
C			C. pick a fixed epoch for et(1) and put the elapsed time since that
C			   epoch in et(2)
C	ntarg	integer	number of the body whose position and velocity are needed
C	ncent	integer	number of the body used as the origin
C			The numbering convention for NTARG and NCENT is:
C			0 = SUN (same as 11)
C			1 = MERCURY		 9 = PLUTO
C			2 = VENUS		10 = MOON
C			3 = EARTH		11 = SUN
C			4 = MARS		12 = NUTATIONS (LONGITUDE AND OBLIQ)
C			5 = JUPITER		13 = LIBRATIONS, IF ON EPH FILE
C			6 = SATURN		14 = SOLAR-SYSTEM BARYCENTER
C			7 = URANUS		15 = EARTH-MOON BARYCENTER
C			8 = NEPTUNE		
C	(For nutations and librations used NTARG=12 and 13 respectively. NCENT is ignored)
C
C	bKm	logical	flag setting the units of the output (default: .FALSE.)
C			bKm	= .TRUE. : KM and KM/SEC
C		    		= .FALSE.: AU and AU/DAY
C			bKm determines the time unit for nutations and librations
C				The angles are always in radians.
C	bBary	logical	flag defining the output center (default: .FALSE.)
C			(affects only the 9 planets)
C			bBary	= .TRUE. : center is solar system barycenter
C				= .FALSE.: center is Sun
C OUTPUTS:
C	rrd(6)	double precision
C			Position (AU or km) and velocity (AU/day or km/s) of NTARG with
C			respect to NCENT
C			For librations the units are radians and radians/day
C			For nutations the first four will be set (units of radians and
C			radians/day).
C	jpl_eph2 logical	return value of call to jpl_init.
C				will be .FALSE. if input time was outside range of ephemeris
C				file or if there was a read error. Otherwise .TRUE.
C CALLS:
C	jpl_init, jpl_state, ArrR8Zero, ArrI4Zero, ArrR8Copy, Say, jpl_message
C RESTRICTIONS:
C	The ephemeris files are organized in records covering periods of 32 days.
C	When the ephemeris is accessed for a large number of times inside a time
C	range less than 32 days, the ephemeris file is accessed only once or twice,
C	reading in 32-days worth of data.
C
C	In this case it may be advantageous to close the ephemeris file
C	explicitly using jpl_close after each jpl_eph2 call. jpl_state will work
C	with the data in memory until a time outside the 32-day range is requested.
C	If this happens the appropriate ephemeris file is opened, and the required
C	record is read. This way the ephemeris file is open only for those brief
C	moments when a new 32-day period is needed. The disadvantage is a speed
C	penalty incurred because of repeated opening of files.
C
C	NOTE: currently the ephemeris file is closed after a record is read,
C	so a call to jpl_close is not needed. To keep the ephemeris file open
C	the call to iFreeLun in jpl_state needs to be commented out.
C PROCEDURE:
C	Currently two ephemeris files are accessed: 1950-2000 and 2000-2050.
C MODIFICATION HISTORY:
C	????, Paul Hick (UCSD/CASS)
C		Based on original JPL package. See:
C		ftp://ssd.jpl.nasa.gov/pub/eph/planets/fortran/
C	FEB-2005, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C		Allow 0 to represent Sun (in addition to 11)
C-
	implicit double precision (A-H,O-Z)

	    double precision et(2)
	    integer	ntarg_
	    integer	ncent_	    
	    double precision rrd(6)
	    logical	bKm
	    logical	bBary

	double precision jd(4)
	double precision tt(2)
	double precision pv(6,15)
	integer		lst(15)

	character	cSay *9			/'jpl_state'/

	call jpl_split(et(1)-0.5d0,jd(1))
	call jpl_split(et(2),jd(3))
	jd(1) = jd(1)+jd(3)+0.5d0
	jd(2) = jd(2)+jd(4)
	call jpl_split(jd(2),jd(3))
	jd(1) = jd(1)+jd(3)

	tt(1) = jd(1)
	tt(2) = jd(4)

	! Make sure that the correct 32-day period of data is accessible
	! (either in memory already or to be read from file).
	! Return with jpl_eph2 = .FALSE. if there is a problem.

	call jpl_init(tt,jpl_eph2)

	if (.not. jpl_eph2) then
	    call jpl_message(cSay,tt,'probably outside range of ephemeris')
	    return
	end if

	ntarg = ntarg_
	ncent = ncent_
	if (ntarg .eq. 0) ntarg = 11	! Sun is in 11
	if (ncent .eq. 0) ncent = 11

	call ArrR8Zero( 6,rrd)
	call ArrI4Zero(15,lst)

	if (ntarg .eq. 12 .or. ntarg .eq. 13) then! Nutations and librations
	    lst(ntarg) = 2
	    call jpl_state(tt,lst,pv,bKm,bBary,EMRAT,jpl_eph2)
	    if (.not. jpl_eph2) then
		call jpl_message(cSay,tt,'read error in ephemeris')
		return
	    end if
	    call ArrR8Copy(6,pv(1,ntarg),rrd)

	else if (ntarg .ne. ncent) then

	    do i=ntarg,ncent,ncent-ntarg	! Set up lst
		if (i .le. 10) lst( i) = 2	! Planets and Moon
		if (i .eq. 10) lst( 3) = 2	! Moon : Earth needed
		if (i .eq.  3) lst(10) = 2	! Earth: Moon needed
		if (i .eq. 15) lst( 3) = 2 	! Earth-Moon barycenter: Earth needed
	    end do

	    ! Make a barycentric call to jpl_state

	    call jpl_state(tt,lst,pv,bKm,.TRUE.,EMRAT,jpl_eph2)
	    if (.not. jpl_eph2) then
		call jpl_message(cSay,tt,'read error in ephemeris')
		return
	    end if
						! Solar system barycenter
	    if (ntarg .eq. 14 .or. ncent .eq. 14) call ArrR8Zero(6,pv(1,14))
						! Earth-Moon barycenter
	    if (ntarg .eq. 15 .or. ncent .eq. 15) call ArrR8Copy(6,pv(1,3),pv(1,15))
						! Earth (3) and Moon (10)
	    if (ntarg*ncent .eq. 30 .and. ntarg+ncent .eq. 13) then
		call ArrR8Zero(6,pv(1,3))

	    else
		if (lst(3) .eq. 2) then
		    do i=1,6
			pv(i,3) = pv(i,3)-pv(i,10)/(1d0+EMRAT)
		    end do
		end if

		if (lst(10) .eq. 2) then
		    do i=1,6
			pv(i,10) = pv(i,10)+pv(i,3)
		    end do
		end if

	    end if

	    do i=1,6
		rrd(i) = pv(i,ntarg)-pv(i,ncent)
	    end do

	end if

	return
	end