;+
; NAME:
;	mpc_orbit_eph
; PURPOSE:
;	Get position for MPC comets and minor planets from orbital elements
; CATEGORY:
;	smei/gen/idl/ephem
; CALLING SEQUENCE:
	FUNCTION mpc_orbit_eph, UT, body, $
		comet		= comet		, $
		center		= center	, $
		location	= location	, $
		speed		= speed 	, $
		to_sphere	= to_sphere , $
		degrees 	= degrees	, $
		precess 	= precess	, $
		to_ecliptic = to_ecliptic,$
		get 		= get		, $
		_extra		= _extra	, $
		silent		= silent
; INPUTS:
;	UT		scalar; type: double
;					Julian date when ephemerides are needed
;	body	scalar; type: string or integer; default: 'C/2004 Q2 (Machholz)'
;					body name or body number
;					(see: mpc_comets and mpc_minor_planets)
; 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)
; 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
;					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, AngleRange
;	mpc_comets, mpc_minor_planets, jpl_eph, jpl_body, CvPrecess
;	CvSky, who_am_i
; PROCEDURE:
;	The main input file is the list of orbital elements from
;	the Minor Planet Center
;		http://www.cfa.harvard.edu/iau/Ephemerides/Comets/Soft00Cmt.txt
; MODIFICATION HISTORY:
;	MAY-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

InitVar, silent , 0
InitVar, get 	, /key
InitVar, comet	, /key

CASE comet OF
0: InitVar, body, 'Ceres'
1: InitVar, body, 'C/2004 Q2 (Machholz)'
ENDCASE

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)

	; Remove invalid MPC names/number

	CASE comet OF
	0: body = mpc_minor_planets(body,_extra=_extra)
	1: body = mpc_comets	   (body,_extra=_extra)
	ENDCASE

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

	rectangular = 1-to_sphere

	FOR ibody=0,nbody-1 DO BEGIN

		PV[*,ibody,*] = mpc_orbit_eph(UT, body[ibody], center=center, $
			/get, degrees=degrees, comet=comet, _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=2000.0d0),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)

CASE comet OF
0: tmp = mpc_minor_planets(body,orbit=orbit,degrees=degrees,_extra=_extra)
1: tmp = mpc_comets	  (body,orbit=orbit,degrees=degrees,_extra=_extra)
ENDCASE

IF tmp NE '' THEN BEGIN
	; The orbital elements are for epoch J2000
	; The position vector is returned in AU; the velocity in AU/day

	vector = KeplerOrbit(UT,elements=orbit[0:3],plane=orbit[4:5],	$
				degrees=degrees, angularvelocity=velocity)

	; Convert from J2000 ecliptic to J2000 equtorial

	T0 = TimeSet(jepoch=2000.0d0)
	vector	 = CvSky(T0, from_ecliptic=vector  , degrees=degrees, /to_equatorial, /silent)
	velocity = CvSky(T0, from_ecliptic=velocity, degrees=degrees, /to_equatorial, /silent)

	; Convert to rectangular coordinates

	vector	 = cv_coord(from_sphere=vector  , degrees=degrees, /to_rect)
	velocity = cv_coord(from_sphere=velocity, degrees=degrees, /to_rect)

	; Combine position and velocity (same as for JPL ephemeris)

	vector = [vector,velocity]

	; Presumably we now have heliocentric equatorial J2000 coordinates
	; It 'center' is specified and 

	IF IsType(center, /defined) THEN BEGIN
		IF IsType(center,/string) THEN tmp = jpl_body(center) ELSE tmp = center
		IF tmp NE jpl_body(/sun) THEN	$
			vector -= jpl_eph(TimeGet(UT,/jd,/scalar),center)
	ENDIF

ENDIF

RETURN, vector  &  END