;+
; NAME:
;	mpc_eph
; PURPOSE:
;	Get position from MPC ephemerides
; CATEGORY:
;	smei/gen/idl/ephem
; CALLING SEQUENCE:
	FUNCTION mpc_eph, UT, body, $
		center		= center	, $
		location	= location	, $
		speed		= speed 	, $
		to_sphere	= to_sphere , $
		degrees 	= degrees	, $
		precess 	= precess	, $
		to_ecliptic = to_ecliptic,$
		source		= source	, $
		get 		= get		, $
		_extra		= _extra	, $
		silent		= silent
; INPUTS:
;	UT		scalar; type: double
;					Julian date when ephemerides are needed
;	body	scalar; type: string or integer; default: 'tempel_1'
;					body name or body number (see: mpc_body)
; OPTIONAL INPUT PARAMETERS:
;	center=center
;			scalar; type: integer; default: jpl_body(/sun)
;					by default, heliocentric coordinates are returned.
;					if center is set to a non-zero value between 1 and 9
;					then the coordinates are centered on the corresponding
;					planet (e.g. center=3 returns geocentric coordinates)
;	source=source
;			scalar; type: string; default: who_am_i(/dir)/mpc
;					By default MPC ephemeris files are looked for in subdirectory
;					mpc of the directory where this source code is located
;					Use this keyword to point to another directory.
; OUTPUTS:
;	Result	array[6] or array[6,*]; type: double
;					array[0:2,*] = position in AU
;					array[3:5,*] = velocity in AU/day
;					If the MPC ephemeris could not be determined
;					(because the ephemeris file could not be located or because
;					the time is outside the range of the ephemeris)
;					then all components are set to the double precision NaN
;					value !values.d_nan.
; OPTIONAL OUTPUT PARAMETERS:
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	IsType, InitVar, IsTime, TimeSet, TimeGet, TimeOp, TimeLimits, TimeUnit
;	TimeInterpol, AngleUnits, AngleRange, mpc_body, jpl_eph, jpl_body
;	CvPrecess, CvSky, who_am_i
; PROCEDURE:
;	Ephemeris files can be retrieved from
;		http://cfa-www.harvard.edu/iau/MPEph/MPEph.html
;	Fill out the form asking for geocentric coordinates:
;       - make sure "Return ephemerides" is set
;		- insert name of minor planet or comet
;		- insert ephemeris start date
;		- set number of dates to output
;       - set the ephemeris interval to 1 (= 1 day; units default to days)
;	Other fields should not matter. Click on "Get ephemeris'; then save the
;	resulting html file and put in the subdirectory 'mpc' of the directory
;	where this source code file is located (or, altenatively, into the directory
;	you intend to specify in the 'source' keyword to this routine.
;
;	The ephemeris record should look like this:
;	Date       UT      R.A. (J2000) Decl.    Delta     r     El.    Ph.   m1     Sky Motion
;	            h m s                                                            "/min    P.A.
;	2005 04 01 000000 13 21 22.5 +12 55 12   0.800   1.772  160.0  11.1  11.2    0.45    286.7
;
;	The ephemeris for a specified time is obtained by linear interpolation
;	on this list.
; MODIFICATION HISTORY:
;	JUL-2005, Paul Hick (UCSD/CASS)
;	JUL-2014, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Updated documentation
;-

InitVar, silent , 0
InitVar, get 	, /key
InitVar, body	, '9p_tempel'
InitVar, source , filepath(root=who_am_i(/dir),'mpc')

IsT = IsTime(UT)
NT = n_elements(UT)

IF NOT get THEN BEGIN

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

	IF not IsT THEN UT = TimeSet(jd=UT)

	body = mpc_body(body,count=nbody)		; Remove invalid MPC names/number, and match file names
	nbody = n_elements(body)

	PV = dblarr(6,nbody,NT)

	rectangular = 1-to_sphere

	FOR ibody=0,nbody-1 DO BEGIN

		PV[*,ibody,*] = mpc_eph(UT, body[ibody], center=center, /get, degrees=degrees, _extra=_extra, silent=silent)

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

		IF precess THEN		$
			PV[0:2,ibody,*] = CvPrecess(TimeSet(jepoch=2000d0),UT, PV[0:2,ibody,*], degrees=degrees, rectangular=rectangular)

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

		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

bad_vector = replicate(BadValue(0.0d0),6)
vector = SuperArray(bad_vector,NT,/trail)

IF flt_read(filepath(root=source,body+'.html'),eph,/double,silent=silent,atleast=10) THEN BEGIN

	; The ephemeris file is assumed to provide geocentric coordinates

	c_year =  0
	c_mon  =  1
	c_day  =  2
	c_hms  =  3

	c_ra   =  4		; Hours, minutes, seconds
	c_dec  =  7		; Degrees, arcminutes, arcseconds
	c_dis  = 10		; Geocentric distance in AU

	uday = TimeUnit(/day)
	hms  = round(reform(eph[c_hms,*]))

	tt_eph = TimeSet(	$
		yr  = round(reform(eph[c_year,*]))	, $	
		mon = round(reform(eph[c_mon ,*]))	, $
		day = round(reform(eph[c_day ,*]))	, $
		hr  = hms/10000 					, $
		minu= (hms mod 10000)/100			, $
		sec = hms mod 100					)

	nt_eph = n_elements(tt_eph)

	valid = where(	TimeOp(/subtract,UT,TimeLimits(tt_eph,/min),uday) GE 0 AND	$
					TimeOp(/subtract,UT,TimeLimits(tt_eph,/max),uday) LE 0 )

	IF valid[0] NE -1 THEN BEGIN

		vec_eph = dblarr(3,nt_eph,/nozero)

		vec_eph[0,*] = AngleUnits(from_time=eph[c_ra :c_ra +2,*],/to_angle,degrees=degrees)
		vec_eph[1,*] = AngleUnits(from_sexa=eph[c_dec:c_dec+2,*],/to_angle,degrees=degrees,/single)
		vec_eph[2,*] = eph[c_dis,*]

		vec_eph = cv_coord(from_sphere=vec_eph, /to_rect, degrees=degrees)

		tt = UT[valid]
		FOR i=0,2 DO vector[i,valid] = TimeInterpol(vec_eph[i,*], tt_eph, tt, uday, _extra=_extra)

		jd = TimeGet(tt, /jd, /scalar)

		CASE IsType(center, /defined) OF
		0: vector[*,valid] = jpl_eph(jd)+vector[*,valid]	; Change from geocentric to heliocentric
		1: BEGIN
			IF IsType(center,/string) THEN tmp = jpl_body(center) ELSE tmp = center
			IF tmp NE jpl_body(/earth) THEN vector[*,valid] = jpl_eph(jd)+vector[*,valid]-jpl_eph(jd,center)
		END
		ENDCASE 

	ENDIF

ENDIF

RETURN, vector  &  END