;+
; NAME:
;	usno_eph
; PURPOSE:
;	Get position and velocity of asteroid from USNO AE98 ephemerides
; CATEGORY:
;	smei/gen/idl/ephem; USNO Asteroid Ephemeris
; CALLING SEQUENCE:
	FUNCTION usno_eph, UT, body, center=center, $
		location	= location	, $
		speed		= speed 	, $
		to_sphere	= to_sphere , $
		degrees 	= degrees	, $
		precess 	= precess	, $
		to_ecliptic = to_ecliptic,$
		get 		= get		, $
		silent		= silent
; INPUTS:
;	UT		scalar; type: double
;					Julian date when ephemerides are needed
;	body	scalar; type: string or integer; default: 'ceres'
;					asteroid name or asteroid number (see: usno_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)
; OUTPUTS:
;	Result	array[6] or array[6,*]; type: double
;					array[0:2,*] = position in AU
;					array[3:5,*] = velocity in AU/day
;					If the asteroid ephemeris could not be determined
;					(because the ephemeris file could not be located or because
;					the Julian day JD is outside the range of the ephemeris)
;					then all six 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, usno_body, usno_init, usno_read, usno_maket, jpl_eph
; COMMON BLOCKS:
;	common USNO_INFO, USNO_PNTR
; RESTRICTIONS:
;	If the center keyword is used then the jpl ephemeris is
;	initialized. Close the ephemeris using jpl_close.
; PROCEDURE:
; >	If the asteroid keyword is not specified than usno_eph calls itself
;	recursively to loop over all available asteroids.
; MODIFICATION HISTORY:
;	SEP-1999, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Modified from USNO AE98 C-software
;-

InitVar, get , /key
InitVar, body, 'ceres'
InitVar, silent, 0

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

	CASE IsT OF
	0: BEGIN
		jd = UT
		UT = TimeSet(jd=jd)
	END
	1: jd = TimeGet(UT, /jd, /scalar)
	ENDCASE

	; Will remove invalid asteroid names/numbers from 'body'

	body = usno_body(body)

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

	rectangular = 1-to_sphere

	J2000 = TimeSet(jepoch=2000.0d0)

	FOR ibody=0,nbody-1 DO BEGIN

		FOR N=0L,NT-1 DO	$
			 PV[*,ibody,N] = usno_eph(jd[N], body[ibody], center=center, /get, silent=silent)

		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 THEN $
				PV[0:2,ibody,*] = CvSky(J2000, from_equator=PV[0:2,ibody,*],$
					/to_ecliptic				, $
					degrees		= degrees		, $
					rectangular	= rectangular	, $
					/silent)

		END

		1: BEGIN

			PV[0:2,ibody,*] = CvPrecess(J2000, 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)

		END

		ENDCASE

		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

jd = UT

common USNO_INFO, USNO_PNTR

failed_vector = replicate(BadValue(0.0d0),6)

ast = usno_init(body,silent=silent)
IF ast EQ -1 THEN RETURN, failed_vector		; Failed to initialize asteroid

tmp = (*(USNO_PNTR[ast])).jdi
IF jd LT tmp THEN BEGIN
	message, /info, 'Julian day'+string(jd)+' is before start of ephemeris at'+string(tmp)
	RETURN, failed_vector
ENDIF

tmp = (*(USNO_PNTR[ast])).jdf
IF jd GT tmp THEN BEGIN
	message, /info, 'Julian day'+string(jd)+' is after end of ephemeris at'+string(tmp)
	RETURN, failed_vector
ENDIF

curjd	= (*(USNO_PNTR[ast])).curjd
curspan	= (*(USNO_PNTR[ast])).curspan

IF jd LT curjd THEN BEGIN					; Move pointer forward

	jdtmp = *((*(USNO_PNTR[ast])).jd) + *((*(USNO_PNTR[ast])).span)

	i = (where(jdtmp GE jd))[0]

	(*(USNO_PNTR[ast])).currec = i

	jdtmp = *((*(USNO_PNTR[ast])).jd   )
	order = *((*(USNO_PNTR[ast])).order) + 1

	j = (where(jdtmp EQ curjd))[0]
	counter = -( (8+2*2)*(j-i+1)+3*round(total(order[i:j]))*8 )

	tmpcounter = -round( total(8+2*2+3*order[i:j]*8) )
	IF tmpcounter NE counter THEN BEGIN
		print, counter, tmpcounter
		message,  'oops 1'
	ENDIF

	move_pointer = 'backward (rec='+strcompress(i,/rem)+')'

ENDIF ELSE IF jd GT curjd+curspan THEN BEGIN; Move pointer back

	jdtmp = *((*(USNO_PNTR[ast])).jd) + *((*(USNO_PNTR[ast])).span)
	order = *((*(USNO_PNTR[ast])).order) + 1

	i = (*(USNO_PNTR[ast])).currec+1
	j = (where(jdtmp GE jd))[0]

	CASE j GT i OF
	0: counter = 0L
	1: counter = (8+2*2)*(j-i)+3*round(total(order[i:j-1]))*8
	ENDCASE

	CASE j GT i OF
	0: tmpcounter = 0L
	1: tmpcounter = round( total(8+2*2+3*order[i:j-1]*8) )
	ENDCASE

	IF tmpcounter NE counter THEN BEGIN
		print, counter, tmpcounter
		message,  'oops 2'
	ENDIF

	(*(USNO_PNTR[ast])).currec = j

	move_pointer = 'forward (rec='+strcompress(i,/rem)+')'

ENDIF ELSE	$								; Don't move pointer

	counter = 0L


IF counter NE 0 THEN BEGIN
	iU = (*(USNO_PNTR[ast])).iU
	point_lun, -iU, pos
	point_lun,  iU, pos+counter
	IF silent LE 0 THEN message, /info, $
		'jd='+strcompress(jd,/rem)	  + $
		'; pntr '+move_pointer+' '+strcompress(pos,/rem)+'+'+strcompress(counter,/rem)+'='+strcompress(pos+counter,/rem)
	usno_read, ast
ENDIF

time = (jd-(*(USNO_PNTR[ast])).curjd)*2/(*(USNO_PNTR[ast])).curspan - 1

usno_maket, time, poscheb, velcheb

order = ( *(USNO_PNTR[ast])).curorder
coef  = ((*(USNO_PNTR[ast])).coef)[0:order,*]

tmp	  = replicate(1,3)
vector = [ total( (poscheb[0:order]#tmp)*coef, 1), total( (velcheb[0:order]#tmp)*coef, 1) ]

vector[3:5] *= 2.0d0/(*(USNO_PNTR[ast])).curspan

IF IsType(center, /defined) THEN vector -= jpl_eph(jd,center)

RETURN, vector  &  END