;+
; NAME:
;	jpl_test
; PURPOSE:
;	Test program for JPL lunar and planetary ephemeris
; CATEGORY:
;	smei/gen/idl/ephem JPL Lunar and Planetary Ephemeris
; CALLING SEQUENCE:
	PRO jpl_test, newcomb=newcomb, aa_equ=aa_equ, aa_ecl=aa_ecl, convert=convert
; OPTIONAL INPUT PARAMETERS:
;	By default a test is run using the testp.405 file provided with the JPL software.
;
;	/newcomb		the JPL heliocentric position of Earth is compared with the
;					result from the function NewcombSun (which is known to reproduce
;					values from the Astronomical Almanac with a precision on the order
;					of arcseconds. The only correction to the JPL values is a precession
;					from J2000 to the equinox of date.
;	 /almanac		compares JPL and NewcombSun results with numbers from
;					the Astronomical Almanac
; OUTPUTS:
;	(none)
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	jpl_init, jpl_eph, jpl_close, InitVar, who_am_i, flt_read, TimeSet
;	AngleUnits, AngleRange, CvSky, TimeXAxis, PlotCurve
;	TimeOp, TimeUnit, ToDegrees, NewcombSun
; PROCEDURE:
;	Checks ephemeris calculations agains numbers in the test file
;		testpo.405 stored in the subdirectory 'jpl'.
; MODIFICATION HISTORY:
;	FEB-2000, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

InitVar, newcomb, /key
InitVar, aa_equ , /key
InitVar, aa_ecl , /key
InitVar, convert, /key

CASE 1 OF

convert: BEGIN

	status = flt_read( filepath(root=who_am_i(/dir),subdir='jpl','almanac.txt'), sun, /double )

	TT = TimeSet(jd=reform(sun[0,*]))

	lng = AngleUnits(from_sexa=sun[1: 3,*], /to_degree)
	lat = reform(sun[4,*])

	ra  = AngleUnits(from_time=sun[5: 7,*], /to_degree)
	dec = AngleUnits(from_sexa=sun[8:10,*], /to_degree, /singlesign)

	PosEqu = transpose([[ra],[dec]])
	PosEcl = CvSky(TT, from_equator=PosEqu, /to_ecliptic, /degrees, /silent)

	dlng = (reform(PosEcl[0,*])       -lng)*3600d0
	dlat = (reform(PosEcl[1,*])*3600d0-lat)

	twin,/del
	twin,xsize=4*200, ysize=3*200

	!p.multi=[0,2,2]

	unit = TimeUnit(/year)

	chars = 1.3

	TimeXAxis, TT, unit, upto=unit
	PlotCurve, TT, dlng, /oplot, /newyaxis, /silent, ytitle='d(ecl long) (arcsec)', chars=chars


	TimeXAxis, TT, unit, upto=unit
	PlotCurve, TT, dlat, /oplot, /newyaxis, /silent, ytitle='d(ecl lat) (arcsec)', chars=chars

	!p.multi=0

END

aa_equ: BEGIN

	status = flt_read( filepath(root=who_am_i(/dir),subdir='jpl','almanac.txt'), sun, /double )

	TT = TimeSet(jd=reform(sun[0,*]))

	;lng = AngleUnits(from_sexa=sun[1: 3,*], /to_degree)
	;lat = reform(sun[4,*]/3600d0)

	ra  = AngleUnits(from_time=sun[5: 7,*], /to_degree)
	dec = AngleUnits(from_sexa=sun[8:10,*], /to_degree, /singlesign)
	dis = reform(sun[11,*])

	AlmEqu = transpose([[ra ],[dec],[dis]])
	SunEqu = CvSky(TT,from_ecliptic=NewcombSun(TT, /degrees, /prec), /to_equator, /degrees, /silent)

	; The JPL ephemeris presumably returns RA and DEC for Epoch J2000.0
	; Precess the values to Epoch TT before

	JplEqu = jpl_eph(TT, jpl_body(/sun), center=jpl_body(/earth),/location,/to_sphere,/degrees)
	jpl_close, /kill
	JplEqu = CvPrecess(TimeSet(jepoch=2000.0d0),TT, JplEqu, /degrees)
	JplEqu[0,*] = AngleRange(JplEqu[0,*], /degrees)


	drasun  = (SunEqu[0,*]-AlmEqu[0,*])*3600d0
	drajpl  = (JplEqu[0,*]-AlmEqu[0,*])*3600d0

	ddecsun = (SunEqu[1,*]-AlmEqu[1,*])*3600d0
	ddecjpl = (JplEqu[1,*]-AlmEqu[1,*])*3600d0

	drsun   = (SunEqu[2,*]-AlmEqu[2,*])*ToDegrees()*3600d0
	drjpl   = (JplEqu[2,*]-AlmEqu[2,*])*ToDegrees()*3600d0

	twin,/del
	twin,xsize=4*200, ysize=3*200

	!p.multi=[0,2,2]

	;unit = TimeUnit(/day)
	unit = TimeUnit(/year)

	chars = 1.3

	TimeXAxis, TT, unit, upto=unit
	PlotCurve, TT, drasun, /oplot, /newyaxis, /silent, ytitle='d(RA) (arcsec)', chars=chars
	PlotCurve, TT, drajpl, /oplot, linestyle=2, /silent


	TimeXAxis, TT, unit, upto=unit
	PlotCurve, TT, ddecsun, /oplot, /newyaxis, /silent, ytitle='d(Dec) (arcsec)', chars=chars
	PlotCurve, TT, ddecjpl, /oplot, linestyle=2, /silent

	TimeXAxis, TT, unit, upto=unit
	PlotCurve, TT, drsun, /oplot, /newyaxis, /silent, ytitle='d(Dis) (arcsec)', chars=chars

	TimeXAxis, TT, unit, upto=unit
	PlotCurve, TT, drjpl, /oplot, /newyaxis, /silent, ytitle='d(Dis) (arcsec)', chars=chars

	!p.multi=0

END

aa_ecl: BEGIN


	status = flt_read( filepath(root=who_am_i(/dir),subdir='jpl','almanac.txt'), sun, /double )

	TT = TimeSet(jd=reform(sun[0,*]))

	lng = AngleUnits(from_sexa=sun[1: 3,*], /to_degree)
	lat = reform(sun[4,*]/3600.0d0)
	dis = reform(sun[11,*])

	AlmEcl = transpose([[lng],[lat],[dis]]) 	; Ecliptic coordinates from almanac
	SunEcl = NewcombSun(TT, /degrees, /prec)

	; The JPL ephemeris presumably returns RA and DEC for Epoch J2000.0
	; Precess the values to Epoch TT before

	JplEcl = jpl_eph(TT, jpl_body(/sun), center=jpl_body(/earth),/location,/to_sphere,/degrees,/precess)
	jpl_close, /kill
	JplEcl = CvSky(TT, from_equator=JplEcl, /to_ecliptic, /degrees, /silent)
	JplEcl[0,*] = AngleRange(JplEcl[0,*], /degrees)

	dlngsun = (SunEcl[0,*]-AlmEcl[0,*])*3600d0
	dlngjpl = (JplEcl[0,*]-AlmEcl[0,*])*3600d0

	dlatsun = (SunEcl[1,*]-AlmEcl[1,*])*3600d0
	dlatjpl = (JplEcl[1,*]-AlmEcl[1,*])*3600d0

	drsun   = (SunEcl[2,*]-AlmEcl[2,*])*ToDegrees()*3600d0
	drjpl   = (JplEcl[2,*]-AlmEcl[2,*])*ToDegrees()*3600d0

	twin,/del
	twin,xsize=4*200, ysize=3*200

	!p.multi=[0,2,2]

	;unit = TimeUnit(/day)
	unit = TimeUnit(/year)

	chars = 1.3

	;TimeXAxis, TT, unit
	;PlotCurve, TT, dlngsun, /oplot, /newyaxis, /silent, ytitle='d(Lng) (arcsec)', chars=chars
	;PlotCurve, TT, dlngjpl, /oplot, linestyle=2, /silent


	;TimeXAxis, TT, unit
	;PlotCurve, TT, dlatsun, /oplot, /newyaxis, /silent, ytitle='d(Lat) (arcsec)', chars=chars
	;PlotCurve, TT, dlatjpl, /oplot, linestyle=2, /silent

	;TimeXAxis, TT, unit
	;PlotCurve, TT, drsun, /oplot, /newyaxis, /silent, ytitle='d(Dis) (arcsec)', chars=chars

	;TimeXAxis, TT, unit
	;PlotCurve, TT, drjpl, /oplot, /newyaxis, /silent, ytitle='d(Dis) (arcsec)', chars=chars

	TimeXAxis, TT, unit, upto=unit
	PlotCurve, TT, dlngjpl, /oplot, /newyaxis, /silent, ytitle='d(Lng) (arcsec)', chars=chars


	TimeXAxis, TT, unit, upto=unit
	PlotCurve, TT, dlatjpl, /oplot, /newyaxis, /silent, ytitle='d(Lat) (arcsec)', chars=chars

	TimeXAxis, TT, unit, upto=unit
	PlotCurve, TT, drjpl, /oplot, /newyaxis, /silent, ytitle='d(Dis) (arcsec)', chars=chars

	!p.multi=0

END

newcomb: BEGIN

	;T0 = TimeSet(yr=2000, doy=266, hour=17, minute=14, sec=2, msec=970.91396d0)
	;T0 = TimeSet(yr=2000, doy=266)
	;dT = gridgen(361,range=180*[-1,1])

	T0 = TimeSet(yr=1960, doy=1)
	dT = gridgen(365L*80)

	;T0 = [TimeSet(yr=2000, doy=86),TimeSet(yr=2000, doy=86, hour=17, minute=14, sec=2, msec=970.91396d0)]
	;dT = 0 ;gridgen(361)

	dT = TimeSet(day=dT)

	TT = TimeOp(/add,T0,dT)
	;JD = TimeGet(TT, /jd)

	PP = NewcombSun(TT, /helio, /degrees, /prec)
	;PP = NewcombSun(TT, /degrees, /prec)

	;for i=0,n_elements(TT)-1 do	$
	;	print, [JD[i],AngleUnits(from_deg=PP[0,i],/to_sexa),PP[1,i]*3600d0]
	;return

	PP = CvSky(TT, from_ecliptic=PP, /to_equator, /degrees, /silent)

	;for i=0,n_elements(TT)-1 do	$
	;	print, JD[i],[AngleUnits(from_deg=PP[0,i],/to_time),AngleUnits(from_deg=PP[1,i],/to_sexa)]
	;return

	; The JPL ephemeris presumably returns RA and DEC for Epoch J2000.0
	; Precess the values to Epoch TT before
	cPP = cv_coord(from_sphere=PP, /to_rect, /degrees)

	QQ = jpl_eph(TT, /location, /to_sphere, /degrees, /precess)
	jpl_close, /kill
	cQQ = cv_coord(from_sphere=QQ, /to_rect, /degrees)
	QQ[0,*] = AngleRange(QQ[0,*], /degrees)


	dra  = (AngleRange(PP[0,*]-QQ[0,*], /degrees, /pi))*3600
	ddec = (PP[1,*]-QQ[1,*])*3600d0
	dr   = (PP[2,*]-QQ[2,*])*ToDegrees()*3600d0

	dp = cQQ-cPP
	dp = sqrt( total(dp*dp,1) )*ToDegrees()*3600d0

	print, 'residuals (mean and standard dev) in arcsec:'
	print, 'ra  residual:', mean(dra ), stddev(dra )
	print, 'dec residual:', mean(ddec), stddev(ddec)
	print, 'r   residual:', mean(dr  ), stddev(dr  )
	print, 'dis residual:', mean(dp  ), stddev(dp  )

	twin,/del
	twin,xsize=4*200, ysize=3*200

	!p.multi=[0,2,2]

	;unit = TimeUnit(/day)
	unit = TimeUnit(/year)

	chars = 1.3

	TimeXAxis, TT, unit, upto=unit
	PlotCurve, TT,  dra, /oplot, /newyaxis, /silent, ytitle='d(RA) (arcsec)', chars=chars

	TimeXAxis, TT, unit
	PlotCurve, TT, ddec, /oplot, /newyaxis, /silent, ytitle='d(Dec) (arcsec)', chars=chars

	TimeXAxis, TT, unit, upto=unit
	PlotCurve, TT, dr, /oplot, /newyaxis, /silent, ytitle='d(R) (arcsec)', chars=chars

	TimeXAxis, TT, unit, upto=unit
	PlotCurve, TT, dp, /ynozero, /oplot, /newyaxis, /silent, ytitle='distance (arcsec)', chars=chars

	!p.multi=0

END

ELSE: BEGIN

	NPT	= 50
	LINE = 0
	IMAP = [0,1,2,3,4,5,6,7,8,9,10,0,13,14,11,12]

	message, /info, ' JPL TEST-EPHEMERIS program. Last modified Feb 2005.'

	EPHFMT = '(15X,D10.1,3I3,D22.13)'
	ET     = 0d0
	NTARG  = 0
	NCTR   = 0
	NCOORD = 1
	XI     = 0.D0

	openr, /get_lun, iu, filepath(root=who_am_i(/dir),subdir='jpl','testpo.405')

	A = ''
	REPEAT readf, iu, A	UNTIL A EQ 'EOT'

	readf, iU, format=EPHFMT, ET,NTARG,NCTR,NCOORD,XI

	WHILE NOT eof(iU) DO BEGIN

		;405  1600.01.01 2305447.5  8  3  1     -26.3227808794400

		NCOORD = NCOORD-1
		NTARG = IMAP[NTARG]
		NCTR  = IMAP[NCTR ]
		R = jpl_eph(ET,NTARG,center=NCTR,/loud,/test,/noclose)

		DEL = abs(R[NCOORD]-XI)
		IF NTARG EQ IMAP[14] AND NCOORD EQ 3 THEN DEL = DEL/(0.23d0*(ET-2451545.d0))

		LINE = LINE+1
		IF LINE mod NPT EQ 0 THEN print, format='(I6,F10.1,3I5,3F20.13)',	$
		     LINE,ET,NTARG,NCTR,NCOORD,XI,R[NCOORD],DEL
						; Print WARNING if difference greater than tolerance.
		IF DEL GE 1.D-13 THEN print, format='(A,/,I6,F10.1,3I5,3F20.13)',	$
			'*****  WARNING : next difference >= 1.D-13  *****',LINE,ET,NTARG,NCTR,NCOORD,XI,R[NCOORD],DEL

		readf, iU, format=EPHFMT, ET,NTARG,NCTR,NCOORD,XI

	ENDWHILE

	jpl_close, /kill

	free_lun, iU

END

ENDCASE

RETURN  &  END