;+
; NAME:
;	jpl_eph
; PURPOSE:
;	Access JPL planetary ephemeris DE405
;	Reads the JPL ephemeris and provides position and velocity vectors of 'body' relative to
;	'center' at the time ET by interpolating on the ephemeris data
; CATEGORY:
;	smei/gen/idl/ephem, JPL Lunar and Planetary Ephemeris
; CALLING SEQUENCE:
	FUNCTION jpl_eph, UT, body, $
		center		= center	, $
		km			= km		, $
		loud		= loud		, $
		location	= location	, $
		speed		= speed		, $
		to_sphere	= to_sphere	, $
		to_ecliptic	= to_ecliptic,$
 		degrees		= degrees	, $
		precess		= precess	, $
		get			= get		, $
		noclose		= noclose	, $
		test		= test
; INPUTS:
;	UT		array[2,n]; type: double
;			scalar	  ; type: double
;			array[n]  ; type: double
;			array[n]  ; type: time structure
;				Ephemeris times specified as Julian date where interpolated positions are needed.
;
;				The input times are usually specified as a pair of numbers in days
;				Generally, ET[0] will be set to some convenient fixed epoch (specified as a
;				Julian day number) used as 'time origin' for the integrations, and ET[1] will
;				be the elapsed time since ET[0] (in days).
;
;				Any time ET[0]+ET[1] falling in the time covered by the ephemeris file is OK
;
;				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]
;
;				Special cases for the input argument are:
;				A double scalar is interpreted as ET = [ET,0]
;				A double array[n] (with n NOT equal 2) is interpreted as n pairs [ET[i],0], i=0,n-1
;				For an array of time structures the midnight, T0, preceding ET[0] is used as epoch;
;				then elapsed times relative to this time origin are calculate to set up the pairs
;				[T0,ET[i]] (i=0,n-1).
;
;	body	scalar		type: integer; default: 3 (Earth)
;				index of the 'target' body
;	center	scalar		type: integer: default: 0 (Sun)
;				index of the center body
;
;				The numbering convention for body and center is:
;				0 = Sun			 8 = Neptune
;				1 = Mercury		 9 = Pluto
;				2 = Venus		10 = Moon
;				3 = Earth		13 = Solar-system barycenter
;				4 = Mars		14 = Earth-Moon barycenter
;				5 = Jupiter		11 = Nutations (longitude and obliquity)
;				6 = Saturn		12 = Librations (if on file)
;				7 = Uranus
;
;				For nutations (body = 11) and librations (body = 12), center is
;				ignored.
;				jpl_state returns the Earth-Moon barycenter state as body 3
;				This is moved to state 14. Then states 14 and 10 (geocentric
;				Moon) are combined to provide the state for Earth.
; OPTIONAL INPUT PARAMETERS:
;	/km			if set Sun, planets and Moon states are returned in
;				km and km/s; default is AU and AU/day
;				For nutations and librations the angle unit is always radians
;	/loud		passed to jpl_init. Prints some info read from ephemeris files
;	/location	only return location of 'body'
;	/speed		only return speed of 'body'
;	/to_sphere	return location in spherical coordinates (longitude, latitude, distance)
;	/noclose	suppresses the call to jpl_close (closing the ephemeris file)
;	/precess
;	/to_ecliptic if set convert from equatorial to ecliptic coordinates
; OUTPUTS:
;	Result		array[6]	type: double
;				contains positions and velocity vector of 'body' relative to
;				'center'. Units depend on setting of /km for Sun, planets and
;				Moon. For librations units are radians and radians/day
;				For nutations only the first four words are set to nutation and
;				nutation rates with units radians and radians/day
;				If the ephemeris could not be accessed or ET was out of the
;				range covered by the ephemeris than R is filled with NaN values
;				(!values.d_nan).
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar, IsTime, TimeGet, TimeSet, TimeUnit, TimeOp, jpl_state
;	jpl_eph, jpl_body, CvPrecess, CvSky, AngleRange
; SEE ALSO:
;	jpl_close, jpl_interp, jpl_init
; COMMON BLOCKS:
	common	JPL_INFO, JPL_PNTR
; RESTRICTIONS:
;	DE405 covers JED 2305424.50  (1599 DEC 09)  to  JED 2525008.50  (2201 FEB 20)
; PROCEDURE:
; >	DE405 includes both nutations and libration, and covers
;
;	JED 2305424.50  (1599 DEC 09)  to  JED 2525008.50  (2201 FEB 20)
;
;	It is based upon the International Celestial Reference Frame (ICRF), which
;	apparently is within 0.01 arcseconds of the dynamical equator and
;	equinox of J2000.
; >	The ephemeris file is opened on the first call to jpl_eph.
; >	Unless /noclose is set, the ephemeris file is closed before returning.
;	The JPL_PNTR heap variable will still keep the last record (with 32 days
;	of ephemeris data) in memory. To destroy the heap variable an explicit
;	call jpl_close,/kill is required.
; MODIFICATION HISTORY:
;	AUG-1999, Paul Hick (UCSD/CASS)
;		Translation of JPL Fortran software provided with the ephemeris package
;	DEC-2002, Paul Hick  (UCSD/CASS)
;		Modified to accept time structures as input, in addition to pairs of
;		double precision numbers for a Julian date.
;	APR-2003, Paul Hick  (UCSD/CASS)
;		Added keywords to_sphere, degrees, speed and location
;	FEB-2005, Paul Hick  (UCSD/CASS; pphick@ucsd.edu)
;		Fairly substantial rewrite to minimize the time that the ephemeris file
;		is kept open.
;	APR-2009, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Bug fix. When body or center where set to 'SS-bary' or 'EM-bary' these
;		were mapped to integers 11 and 12, instead of 13 and 14.
;-

InitVar, get , /key
InitVar, body, jpl_body(/earth) 				; Default body is Earth

IsT = IsTime(UT)
NT  = n_elements(UT)
IF NOT IsT THEN $
	IF (size(UT,/dim))[0] EQ 2 THEN NT = NT/2	; If 1st dimension has 2 elements

IF NOT get THEN BEGIN

	InitVar, speed		, /key
	InitVar, location	, /key
	InitVar, to_sphere	, /key
	InitVar, to_ecliptic, /key
	InitVar, precess	, /key
	InitVar, noclose	, /key

	CASE IsT OF
	0: BEGIN
		ET = [UT]
		IF size(ET,/n_dim) EQ 1 THEN ET = [reform(ET,1,NT), dblarr(1,NT)]
		UT_ = TimeSet(jd=total(ET,1))
	END
	1: BEGIN
		unit = TimeUnit(/day)
		ET = TimeGet(UT[0],/botime,unit)
		ET = [replicate(TimeGet(ET,/jd,/scalar),1,NT),reform([TimeOp(/subtract,UT,ET,unit)],1,NT)]
		UT_ = UT
	END
	ENDCASE

	nbody = n_elements(body)
	PV = dblarr(6,nbody,NT)

	FOR N=0L,NT-1 DO PV[*,*,N] = jpl_eph(ET[*,N], body, center=center, km=km, loud=loud, /get,test=test)
	IF NOT noclose THEN jpl_close	; Closes ephemeris file (does NOT destroy JPL_PNTR)

	rectangular = 1-to_sphere

	J2000 = TimeSet(jepoch=2000.0d0)

	FOR ibody=0,nbody-1 DO BEGIN

		; The JPL ephemeris returns rectangular J2000 equatorial
		; (RA/dec) coordinates.

		; If /to_sphere is set, then convert to spherical coordinates.

		IF to_sphere THEN	$
			PV[0:2,ibody,*] = cv_coord(from_rect=reform(PV[0:2,ibody,*]), /to_sphere, degrees=degrees)

		CASE precess OF

		0: BEGIN

			; If /to_ecliptic is set then convert from equatorial
			; to ecliptic coordinates, giving ecliptic J2000 coordinates.

			IF to_ecliptic THEN	$
				PV[0:2,ibody,*] = CvSky(J2000,from_equator=PV[0:2,ibody,*],	$
					/to_ecliptic, degrees=degrees, rectangular=rectangular,	$
					/silent)
		END

		1: BEGIN

			; If /precess is set then precess coordinates from J2000 to
			; current equinox.

			PV[0:2,ibody,*] = CvPrecess(J2000,UT_, PV[0:2,ibody,*], degrees=degrees, rectangular=rectangular)

			; If /to_ecliptic is set then convert from equatorial
			; to ecliptic coordinates, giving ecliptic coordinates
			; for the equinox of date.

			IF to_ecliptic THEN	$
				PV[0:2,ibody,*] = CvSky(UT_,from_equator=PV[0:2,ibody,*],	$
					/to_ecliptic, degrees=degrees, rectangular=rectangular,	$
					/silent)

		END

		ENDCASE

		; If /to_sphere is set then make sure that RA/longitude are
		; in [0,360].

		IF to_sphere THEN	$
			PV[0,ibody,*] = AngleRange(PV[0,ibody,*], degrees=degrees)

	ENDFOR

	CASE 1 OF
	location: PV = PV[0:2,*,*]
	speed   : PV = PV[3:5,*,*]
	ELSE	:
	ENDCASE

	RETURN, reform(PV)		; Remove degenerate dimensions

ENDIF

; The remaining part is executed only with single times (NT=1)

IF n_elements(UT) EQ 1 THEN ET = [UT[0],0.0d0] ELSE ET = UT[0:1]
ET = double(ET)

S1 = ET[0]-0.5d0
I1 = floor(S1)
S2 = ET[1]
I2 = floor(S2)
ET = [I1,S1-I1,I2,S2-I2]

ET[0] += ET[2]+0.5d0
ET[1] += ET[3]

S1 = ET[1]
I1 = floor(S1)
ET[2:3] = [I1,S1-I1]
ET[0] += ET[2]

ET = [ET[0],ET[3]]

IF NOT jpl_init(ET, JPL_PNTR, loud=loud, test=test) THEN BEGIN
	message, /info, 'no ephemeris at '+TimeGet( TimeSet(jd=total(ET)), /ymd )
	RETURN, replicate(BadValue(0d0),6)
ENDIF

InitVar, km 	, /key							; Select output units
InitVar, center , jpl_body(/sun) 				; Default center is Sun

LST = intarr(15)

; jpl_body will map the Solar system barycenter (string 'SS-bary)
; to 11 and the Earth-Moon barycenter (string 'EM-bary' to 12).
; 2 is added to map these to 13 and 14, respectively.

CASE IsType(body,/string) OF
0: ibody = body
1: BEGIN
	ibody = jpl_body(body)
	n = where(ibody EQ 11 OR ibody EQ 12)
	IF n[0] NE -1 THEN ibody[n] += 2
END
ENDCASE

CASE IsType(center,/string) OF
0: icenter = center
1: BEGIN
	icenter = jpl_body(center)
	n = where(icenter EQ 11 OR icenter EQ 12)
	IF n[0] NE -1 THEN icenter[n] += 2
END
ENDCASE

LST[ibody  ] = 2	; LST[13] not used
LST[icenter] = 2

IF LST[ 3] EQ 2 THEN LST[10] = 2				; If Earth included then also include Moon
IF LST[10] EQ 2 THEN LST[ 3] = 2				; If Moon included then also include Earth
IF LST[14] EQ 2 THEN LST[ 3] = 2				; jpl_state returns the Earth-Moon barycenter in position 3

; Sun,9 planets,Moon,nutation,libration

PV = jpl_state(ET, JPL_PNTR, LST[0:12], km=km, /ssbary, loud=loud)

PV = [[PV],[0,0,0,0,0,0],[PV[*,3]]]				; Append solar system barycenter, Earth-Moon barycenter
EMRAT = (*JPL_PNTR).EMRAT

; The location of the Moon in 10 is relative to Earth.
; Convert to locations relative to SS barycenter for all bodies.
; Also fill 3 with location Earth relative to SS barycenter.

Moon2Earth = PV[*,10]							; Moon relative to Earth
PV[*, 3] = PV[*,14]-Moon2Earth/(1.0D0+EMRAT)		; Earth relative to SS barycenter
PV[*,10] = PV[*, 3]+Moon2Earth					; Moon relative to SS barycenter

; All bodies are now specified relative to the solar system barycenter
; Take difference with selected center (but don't touch nutation and libration)

LST = dindgen(15)
PV -= PV[*,icenter]#(LST NE 11 AND LST NE 12)

; If the Earth, Moon or EM barycenter are selected for the center then we
; recalculate these positions using only the Moon-Earth distance (using PV[*,14]
; as above ruins the precision in this case).

CASE icenter OF
 3: PV[*,[3,10,14]] =  Moon2Earth#[0.0D0,1.0D0,1.0D0/(1.0D0+EMRAT)]
10: PV[*,[3,10,14]] = -Moon2Earth#[1.0D0,0.0D0,EMRAT/(1.0D0+EMRAT)]
14: PV[*,[3,10,14]] = (Moon2Earth/(1.0D0+EMRAT))#[-1.0D0,EMRAT,0.0D0]
ELSE:
ENDCASE

RETURN, reform(PV[*,ibody])  &  END				; Extract the requested bodies.