C+
C NAME:
C	jpl_state
C PURPOSE:
C	Reads and interpolates ephemeris
C	(internal use by jpl_eph2 only)
C CATEGORY:
C	gen/for/lib
C CALLING SEQUENCE:
	subroutine jpl_state(et,lst,pv,bKm,bBary,EMR,bOK)
C INPUTS:
C	et(2)	double precision
C		Julian ephemeris epoch at which interpolation is wanted.
C		Any combination et(1)+et(2) insided the time span covered by the file
C		is allowed:
C		A. put entire epoch in et(1) and set et(2)=0
C		B. for maximum interpolatin accuracy, set et(1) to the most recent midnight
C		   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	lst(15)	integer
C		array specifying what interpolation is wanted for each of the bodies
C		=0, no interpolation
C		=1, position only
C		=2, position and velocity
C	The designation of the astronomical bodies is the same as used for arguments
C	NTARG and NCENT in jpl_eph
C OUTPUTS:
C	pv(6,15)	double precision
C		array with the interpolated quantities
C		THE BODY SPECIFIED BY lst(I) WILL HAVE ITS STATE IN THE ARRAY STARTING AT
C		pv(1,I).  (ON ANY GIVEN CALL, ONLY THOSE WORDS IN 'pv' WHICH ARE AFFECTED BY
C		THE FIRST 10 'lst' ENTRIES (AND BY lst(14) IF LIBRATIONS ARE ON THE FILE)
C		ARE SET.  THE REST OF THE 'pv' ARRAY IS UNTOUCHED.)  THE ORDER OF COMPONENTS
C		STARTING IN pv(1,I) IS: X,Y,Z,DX,DY,DZ.
C
C		ALL OUTPUT VECTORS ARE REFERENCED TO THE EARTH MEAN EQUATOR AND EQUINOX OF J2000
C		IF THE DE NUMBER IS 200 OR GREATER; OF B1950 IF THE DE NUMBER IS LESS THAN 200.
C
C		The  Moon state is always geocentric. THE OTHER NINE STATES ARE EITHER HELIOCENTRIC
C		OR SOLAR-SYSTEM BARYCENTRIC, DEPENDING ON THE SETTING OF bBary
C
C	pv(6,11)ARRAY CONTAINING THE BARYCENTRIC POSITION AND VELOCITY OF THE SUN.
C	pv(6,12)NUTATIONS, IF ON FILE, ARE PUT INTO pv(K,12) IF lst(12) IS 1 OR 2.
C		THE ORDER OF QUANTITIES IN pv IS:
C 			D PSI  (NUTATION IN LONGITUDE)
C			D EPSILON (NUTATION IN OBLIQUITY)
C			D PSI DOT
C			D EPSILON DOT
C	pv(6,13)LUNAR LIBRATIONS, IF ON FILE, ARE PUT INTO pv(K,13) IF lst(13) IS 1 OR 2.
C CALLS:
C	jpl_interp, CvR8, bOpenFile
C SEE ALSO:
C 	jpl_eph, jpl_init
C INCLUDE:
	implicit double precision (A-H,O-Z)

	include		'dirspec.h'
	include		'filparts.h'
	include		'openfile.h'
C PROCEDURE:
C MODIFICATION HISTORY:
C-
	    double precision et(2)
	    integer	lst(15)
	    double precision pv(6,15)
	    logical	bKm
	    logical	bBary
	    double precision EMR
	    logical	bOK

	character	cSay*9		/'jpl_state'/

	integer		IPT(3,13)
	integer		NCM(13)		/11*3,2,3/
	character	TTL(14,3)*6
	character	CNAM(400)*6
	double precision T(2)
	double precision BUF(1500)
	double precision SS(3)
	double precision CVAL(400)

	parameter	(NCOEFFS = 1018)
	parameter	(IOPN = OPN__TRYINPUT+OPN__NOPARSE+OPN__ONEPASS+OPN__STOP+OPN__READONLY+OPN__NOMESSAGE)

	logical		bOpenFile
	logical		bOSFindClose
	character	cSearch*(FIL__LENGTH)
	character	cFound *(FIL__LENGTH)

	integer		iU		/FIL__NOUNIT/
	integer		NRL		/0/
	save		iU, NRL, BUF, IPT, AU, SS, EMRAT, NR

	EMR = EMRAT				! Output argument

	if (NR .ne. NRL) then			! Read correct record if not in core
	    NRL = NR
	    read (iU,rec=NR,iostat=I) (BUF(K),K=1,NCOEFFS)
	    bOK = I .eq. 0
	    iU  = iFreeLun(iU)
	    if (.not. bOK) return		! Read error
	    call Say(cSay,'I','JPL','new ephemeris record read into core')
	    call CvR8(OS__UNIX,NCOEFFS,BUF)
	else
	    bOK = .TRUE.
	end if

	T(1) = ((et(1)-(dble(NR-3)*SS(3)+SS(1)))+et(2))/SS(3)

	if (bKm) then
	    T(2) = SS(3)*86400.D0
	    AUFAC= 1.D0
	else
	    T(2) = SS(3)
	    AUFAC= 1.D0/AU
	end if

	lst(11) = 2				! Always do the Sun

	do I=13,1,-1
	    if (lst(I)*IPT(2,I) .gt. 0) then
		call jpl_interp(BUF(IPT(1,I)),T,IPT(2,I),NCM(I),IPT(3,I),lst(I),pv(1,I))
		if (I .le. 11) then		! <=11: Planets, Moon and Sun
		    do J=1,2*NCM(I)
			pv(J,I) = pv(J,I)*AUFAC	! <= 9: Planets
			if (I .le. 9 .and. .not. bBary) pv(J,I) = pv(J,I)-pv(J,11)
		    end do
		end if
	    end if
	end do

	return
C+
C NAME:
C	jpl_init
C PURPOSE:
C	Initializes the JPL planetary ephemeris and opens an ephemeris file
C	(internal use by jpl_eph2 only)
C CALLING SEQUENCE:
	entry jpl_init(et,bOK)
C INPUTS:
C	et(2)	double precision	ephemeris time
C OUTPUTS:
C	bOK	logical			.TRUE.: ephemeris data available
C						(it is OK to call jpl_state)
C					.FALSE.: no ephemeris available.
C						(do not bother to call jpl_state)
C SIDE EFFECTS:
C	The variable NR (record for time et) is set here and will be used by jpl_state.
C	The ephemeris file DOES NOT need to be open. As long as the record
C	still in memory matches NR, jpl_state has all the information it needs.
C CALLS:
C	iFilePath, bOpenFile, iSearch, CvR8, CvI4, jpl_inside, iFreeLun
C SEE ALSO:
C 	jpl_state, jpl_interp, jpl_split
C PROCEDURE:
C	Entry point in href=jpl_state=
C	The ephemeris files are located in $EPHEM/jpl
C
C	iRecl = 1*2*NCOEFFS	! Multiplier for recordlength in open statement
C	if (cOpSys .ne. OS__VMS) iRecl = 4*2*NCOEFFS
C	open (iU, file=cFile, access='DIRECT', form='UNFORMATTED', recl=iRecl, status='OLD')
C MODIFICATION HISTORY:
C	???-????, Paul Hick (UCSD/CASS)
C	FEB-2005, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C		Substantial rewrite.
C-
	if (NRL .ne. 0) then

	    ! There still is a record in memory (even if the ephemeris
	    ! file has been close already by a call to jpl_close. Since jpl_state
	    ! will not have to read a new record this is OK.

	    NR = jpl_inside(et,SS)
	    bOK = NR .eq. NRL
	    if (bOK) return
	end if

	if (iU .ne. FIL__NOUNIT) then

	    ! The ephemeris file is still open. If the record required for et is in
	    ! this file than jpl_state will read it.

	    NR = jpl_inside(et,SS)
	    bOK = NR .ne. 0
	    if (bOK) return				! Correct ephemeris file already open
	    iU = iFreeLun(iU)				! Close wrong ephemeris file
	end if

	! Need to find and open the ephemeris file for time et

	bOK = .FALSE.
	I = iFilePath(cEnvi//'EPHEM',1,'jpl','unxp*.405',cSearch)

	I = 1
	do while (iSearch(I,cSearch,cFound) .eq. 1)	! Loop over all ephemeris files

	    ! Stop if the open fails (we know the file exists, so it is probably corrupt)

	    iRecl = 2*NCOEFFS				! Recordlength in 4-byte longwords
	    if (bOpenFile(IOPN,iU,cFound,iRecl)) then
		read (iU,rec=1) TTL,CNAM,SS,NCON,AU,EMRAT,((IPT(I,J),I=1,3),J=1,12),NUMDE,(IPT(I,13),I=1,3)
		read (iU,rec=2) CVAL

		NRL = 0

		call CvR8(OS__UNIX,3,SS    )
		call CvI4(OS__UNIX,1,NCON  )
		call CvR8(OS__UNIX,1,AU    )
		call CvR8(OS__UNIX,1,EMRAT )
		call CvI4(OS__UNIX,1,NUMDE )
		call CvI4(OS__UNIX,3*13,IPT)
		call CvR8(OS__UNIX,400,CVAL)
	    end if

	    ! Now that the ephemeris file is open, make sure it is the right one.
	    ! Note that since NRL has been set to zero, jpl_state is forced to read record NR.

	    NR = jpl_inside(et,SS)
	    bOK = NR .ne. 0

	    if (bOK) then
		if (bOSFindClose()) continue
		!write (*,'(/3F14.2/200(/2(A8,D24.16)))') SS,(CNAM(I),CVAL(I),I=1,NCON)
		!write (*,*) '  line -- jed --  t#   c#   x#   --- jpl value ---'//
     		!&			'   --- user value --    -- difference --'
     		return					! Correct ephemeris file is now open
	    end if					! (but NRL=0 so BUF is not loaded yet.

	    iU = iFreeLun(iU)				! Close wrong ephemeris file
	    I = 0
	end do

	return
C+
C NAME:
C	jpl_close
C PURPOSE:
C	Closes ephemeris file
C CALLING SEQUENCE:
	entry jpl_close()
C CALLS:
C	iFreeLun
C PROCEDURE:
C	Read explanation in jpl_eph2. Currently this function is useless.
C MODIFICATION HISTORY:
C	FEB-2005, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C-

	if (iU .ne. FIL__NOUNIT) iU = iFreeLun(iU)

	return
	end