FUNCTION jpl_state, ET, JPL_PNTR, LST, km=km, ssbary=ssbary, loud=loud

;+
; NAME:
;	jpl_state
; PURPOSE:
; 	Reads and interpolates the JPL planetary ephemeris file
;	(internal use for jpl_eph only)
; CATEGORY:
;	smei/gen/idl/ephem; JPL Lunar and Planetary Ephemeris
; CALLING SEQUENCE:
;	R = jpl_state(ET,LST, km=km,ssbary=ssbary)
; INPUTS:
;	ET			array[2]; type: double (input is converted to double)
;					Julian ephemeris date where the interpolated positions are needed.
;					Any time ET[0]+ET[1] falling in the time covered by the ephemeris file is OK
;
;				Generally, ET[0] will be set to some convenient fixed epoch used as 'time
;				origin' for the integrations, and ET[1] will be the elapsed time since ET[0]
;
;				For maximum precision:
;				set ET[0] to the most recent midnight at or before the interpolation epoch and
;				set ET[1] to the fractional part of the day elapsed since ET[0]
;
;	LST		array[13]	type: integer
;						specifies what type of interpolation is needed for each body
;						= 0: no interpolation
;						= 1: position only
;						= 2: position and velocity
;						LST[0] (Sun) is internally set to 2
;
;				The designation of astronomical bodies in LST
;					0 = Sun			 7 = Uranus
;					1 = Mercury		 8 = Neptune
;					2 = Venus		 9 = Pluto
;					3 = Earth		10 = geocentric Moon
;					4 = Mars		11 = nutations (longitude and obliquity)
;					5 = Jupiter		12 = librations (if present on file)
;					6 = Saturn
; OPTIONAL INPUT PARAMETERS:
;	/km			if set the units for PV[0:10] are in km; default is AU
;	/ssbary		if set then PV[1:9] are returned in solar-system barycentric
;					coordinates; default is heliocentric
;	/loud		passed to jpl_init. Prints some info read from ephemeris files
; OUTPUTS:
;	PV		array[6,13]	type: double precision
;						requested interpolated quantities
;						The order of the components i: X,Y,Z,DX,DY,DZ.
;						Only those components are calculated as indicated by the LST array.
;						Everthing else is set to zero.
;				If the ephemeris file could not be accessed of JD is out of the range
;				covered by the ephemeris then PV will contain NaN values (!values.d_nan)
;
;				All output vectors are referenced to the Earth mea equator and equinox
;				of J2000 if the DE number is 200 or greater; of B1950 if the DE number
;				is less than 200.
;
;				PV[*,0]		the Sun state is always solar-system barycentric
;							[PLH: I think, see how its used in jpl_eph].
;				PV[*,1:9]	the planet states are heliocentric or solar system barycentric,
;							depending on the setting of ssbary
;							[PLH: actually I think PV[*,3] is the Earth-Moon barycenter;
;							see how its used in jpl_eph]
;				PV[*,10]	the Moon state is always geocentric
;
;				PV[0:3,11]	contains nutations and rates. The components are
;							dPSI		nutation in longitude (radians)
;							dEPSILON	nutation in obliquity (radians)
;							dPSI/dT		rates (radians/day)
;							dEPSILON/dT
;				PV[*,12]	lunar librations (if on file) if LST[12] is 1 or 2.
;							in radians and radians/day
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar, jpl_interp
; PROCEDURE:
; MODIFICATION HISTORY:
;	AUG-1999, Paul Hick (UCSD/CASS)
;		Based on JPL Fortran software
;	DEC-2002, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Bug fix. Fraction of day in ET[1] was not always processed properly.
;-

InitVar, km 	, /key
InitVar, ssbary , /key

NR = (*JPL_PNTR).NRL
SS = (*JPL_PNTR).SS

T = ((ET[0]-(double(NR)*SS[2]+SS[0]))+ET[1])/SS[2]

IF km THEN BEGIN			; Units: km/s
	T = [T,SS[2]*86400d0]
	AUFAC= 1.D0
ENDIF ELSE BEGIN			; Units: AU/day
	T = [T,SS[2]]
	AUFAC= 1.D0/(*JPL_PNTR).AU
ENDELSE

PV		= dblarr(6,13)
NCM 	= replicate(3,13)
NCM[11] = 2					; Nutations
LST[ 0] = 2					; Sun
IPT		= (*JPL_PNTR).IPT

FOR I=0,12 DO BEGIN
	IF LST[I]*IPT[1,I] GT 0 THEN BEGIN
		J = 2*NCM[I]-1
		PV[0:J,I] = jpl_interp((*JPL_PNTR).BUF[IPT[0,I]-1:*],T,IPT[1,I],NCM[I],IPT[2,I],LST[I])
		IF I LE 10 THEN BEGIN
			PV[0:J,I] = PV[0:J,I]*AUFAC
			IF I GT 1 AND NOT ssbary THEN PV[0:J,I] = PV[0:J,I]-PV[0:J,0]
		ENDIF
	ENDIF
ENDFOR

RETURN, PV  &  END