;+
; NAME:
;	big_eph
; PURPOSE:
;	Controls all available ephemeris data
; CATEGORY:
;	gen/idl/ephem
; CALLING SEQUENCE:
	FUNCTION big_eph, tt		, $
		jpls		= jpls		, $
		usnos		= usnos		, $
		mpcs		= mpcs		, $
		stars		= stars		, $
		m_min		= m_min 	, $
		m_max		= m_max 	, $
		body		= body		, $
		center		= center	, $
		to_sphere	= to_sphere	, $
		precess		= precess	, $
		degrees		= degrees	, $
		names		= names		, $
		nodist		= nodist	, $
		count		= count		, $
		clean		= clean 	, $
		onebody		= onebody	, $
		sc_loc		= sc_loc	, $
		state		= state 	, $
		_extra		= _extra
; INPUTS:
;	tt			array; type: time structure
;					list of times in UT
; OPTIONAL INPUT PARAMETERS:
;	Any combination of the following keywords determines which
;	bodies are used; if the selection of keywords selects a body
;	more than once, then the duplicate copy is removed from the list.
;
;	/jpls			if set, positions for Sun, Moon, and all 9 planets
;						are returned from the JPL ephemeris
;	/usnos			if set, positions of all asteroids from the USNO
;						ephemeris are returned
;	/mpcs			if set, positions for all minor planets and comets from
;						locally available ephemerides from the Minor Planet
;						Center are returned
;	/stars			if set, positions for all stars (brighter than magnitude
;						m_max) from the SMEI star catalogues are returned
;	m_max=m_max 	scalar: type: float; default: none
;						used only if /stars is set. Only stars with magnitudes
;						brighter than m_max are returned.
;	body=body		scalar; type: string; default: none
;						comma-separated list of bodies for which positions are
;						required. The names are compared against names in the
;						JPL, USNO, MPC and SMEI star databases. The comparison
;						is case-insensitive.
;	
;	center=center	scalar; type: string or integer; default: jpl_body(/sun)
;						any of the bodies in the JPL ephemeris, or the name of
;						of an Earth-orbiting satellite. This determines the origin
;						of the coordinates system used. See PROCEDURE.
;
;	/to_sphere		if set, spherical coordinates are returned. By default,
;						rectangular coordinates are used
;	_extra=_extra	passed to href=CvSky=
;						Should be used to set one of the to_* keywords allowed
;						for CvSky to determine the output coordinate frame.
;						Note that this means that CvSky determines what the
;						default coordinate frame for the output vectors are.
;						Currently the default in CvSky (if no to_* keywords are
;						specified) is ecliptic coordinates.
;	/precess		if set, coordinates are precessed to the equinox at time
;						'tt'. By default, J2000 coordinates are returned.
;	/degrees		if set, all angles will be in degrees. Default is radians
;
;	/clean			if set, the output arrays are reduced by removing all bodies
;						and all times for which no valid positions were obtained.
;	/onebody		can be used to remove the dummy second dimension if
;						only one body is specified, i.e. an array[3,ntime] is
;						returned, instead of an array[3,1,ntime]
; OUTPUTS:
;	Result			array[3,nbody,ntime]; type: double
;						positions of requested bodies. The distance scale is in AU.
;						By default the ephemeris returns J2000 ecliptic
;						rectangular coordinates (see keyword _extra).
;						This is modified with keywords /precess, /to_sphere and
;						_extra=_extra
;	tt				array; type: time structure
;						if /clean is set then this is the subset of the
;						input times for which valid positions were founds
;						if /clean is NOT set the input tt is not modified
; OPTIONAL OUTPUT PARAMETERS:
;	names=names 	array; type: string
;						names of bodies for which positions are returned
;	nodist=nodist	array; type: byte
;						0 indicates a planet, comet or asteroid
;						1 indicates a star
;						For stars the distance is arbitrarily set to one.
;	count=count 	scalar; type: integer
;						# bodies for which positions are returned.
; INCLUDE:
	@compile_opt.pro
; CALLS:
;	InitVar, IsType, destroyvar, where_common, BadValue, TimeSet
;	HOSOrbit, UlyssesOrbit, CvPrecess, CvSky, StereoAOrbit, StereoBOrbit
;	AngleRange
;	jpl_body, jpl_eph, usno_body, usno_eph, usno_close, mpc_body, mpc_eph
;	sgp4_eph, smei_star_info, SuperArray, smei_star_alias
;	TimeSystem, UlyssesOrbit, big_eph_boost, big_eph_clean
;	big_eph_short, sgp_body, sgp_alias
; PROCEDURE:
;	The 'center' keyword determines the origin of the coordinate system.
;	The default is heliocentric coordinates. This can be changed to
;	any of the bodies in JPL ephemeris, i.e. setting center=jpl_body(/earth)
;	forces geocentric coordinates.
;	Center can also be set to an Earth-orbiting spacecraft by specifying the
;	name of an available file with orbital elements. The spacecraft position
;	is taken into account using the sgp4_eph routine.
;	Currently the only satelite we have is Coriolis
;	NOTE: for stars the 'center' keyword is ignored
; MODIFICATION HISTORY:
;	JUL-2005, Paul Hick (UCSD/CASS)
;	OCT-2006, Paul Hick (UCSD/CASS)
;		If one of the ephemeris program would return a bad position vector
;		(e.g. when an MPC ephemeris is accessed outside the time range of the
;		MPC file) then the conversion to spherical coordinates would set the
;		angles to zero with only the radial distance staying bad. Fixed it
;		so the angles are now bad too.
;		The procedure now also keeps track of which bodies have already
;		been identified to avoid unnecessarily accessing an ephemeris.
;	MAR-2007, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Added /precess keyword to sgp4_eph calls to make sure J2000
;		coordinates are returned.
;-

CASE IsType(state,/defined) OF

0: BEGIN

	InitVar, tt   , TimeSystem(/silent)
	InitVar, jpls , /key
	InitVar, usnos, /key
	InitVar, mpcs , /key
	InitVar, stars, /key
	InitVar, body , ''

END

1: BEGIN

	tt     = state.time

	jpls   = state.jpls
	usnos  = state.usnos
	mpcs   = state.mpcs
	stars  = state.stars
	body   = state.body

	center = state.center
	m_min  = state.m_min
	m_max  = state.m_max

END

ENDCASE

InitVar, clean, /key

j2000 = TimeSet(jepoch=2000.0d0)
ntt = n_elements(tt)

; Spacecraft are Earth-orbiting. A file with TLEs should be available
; (is accessed by sgp4_eph).

IF IsType(center,/string) THEN BEGIN

	; If center is not a JPL body, then it better be an Earth-orbiting satellite

	pp = sgp_alias(strtrim(center,2))
	IF (where(pp EQ strlowcase(sgp_body())))[0] NE -1 THEN BEGIN
		sc_center = pp
		center = 'Earth'
	ENDIF

ENDIF

destroyvar, names, nodist

IF jpls THEN BEGIN

	list = (jpl_body(count=n))[0:n-1]

	pp = jpl_eph(tt,list,/location,center=center)

	big_eph_boost, ntt, rr, pp, names, list, nodist, 0B

ENDIF

IF usnos THEN BEGIN

	list = usno_body()

	pp = usno_eph(tt,list,/location,center=center,/silent)
	usno_close

	big_eph_boost, ntt, rr, pp, names, list, nodist, 0B

ENDIF

IF mpcs THEN BEGIN

	list = mpc_body()

	pp = mpc_eph(tt,list,/location,center=center,silent=2)

	big_eph_boost, ntt, rr, pp, names, list, nodist, 0B

ENDIF

IF stars THEN BEGIN

	list = smei_star_info(/get_struct,/silent,count=n,cat='*'); Get star structures

	tmp = smei_star_info(list,/get_mags)		; Extract smei magnitudes

	; Limit the range of magnitudes, if m_min and/or
	; m_max are set.

	CASE 1 OF
	IsType(m_min,/defined) AND IsType(m_max,/defined):	$
							tmp = where(m_min LE tmp AND tmp LE m_max, n)
	IsType(m_min,/defined): tmp = where(m_min LE tmp				 , n)
	IsType(m_max,/defined): tmp = where(				 tmp LE m_max, n)
	ELSE				  : tmp = lindgen(n)
	ENDCASE

	IF tmp[0] NE -1 THEN BEGIN

		list = list[tmp]						; Remaining structures

		pp = dblarr(3,n,/nozero)
		pp[0:1,*] = smei_star_info(list,degrees=degrees)
		pp[  2,*] = 1.0d0
		pp = SuperArray(pp,n_elements(tt),/trail)

		list = smei_star_info(list,/get_name)
		big_eph_boost, ntt, rr, pp, names, list, nodist, 1B

	ENDIF

ENDIF

IF body NE '' THEN BEGIN

	; Only bodies with recognizable (case-insensitive) names are picked up here,
	; so typos will make entries disappear.

	bodies = strtrim(strlowcase(strtok(body,',',/extract)),2)
	bodies = bodies[uniq(bodies,sort(bodies))]

	bodies_found = bytarr(n_elements(bodies))

	; The deep-space spacecraft routines return spherical
	; heliocentric ecliptic coordinates for the current equinox.
	; Needed are rectangular equatorial J2000 coordinates so the
	; JPL ephemeris can be used to shift the origin to 'center'.

	IF (where(bodies_found EQ 0))[0] NE -1 THEN BEGIN

		tmp = where_common(bodies,big_body('Helios B',/alias))

		IF tmp[0] NE -1 THEN BEGIN

			bodies_found[tmp] = 1

			pp = HOSOrbit(tt,hos=1,degrees=degrees)
			pp = CvPrecess(tt,j2000,pp,/from_ecliptic,degrees=degrees)
			pp = cv_coord(from_sphere=pp,/to_rect,degrees=degrees)

			IF IsType(center,/defined) THEN pp -= jpl_eph(tt,center,/location)

			big_eph_boost, ntt, rr, pp, names, 'Helios A', nodist, 0B

		ENDIF

	ENDIF

	IF (where(bodies_found EQ 0))[0] NE -1 THEN BEGIN

		tmp = where_common(bodies,big_body('Helios B',/alias))

		IF tmp[0] NE -1 THEN BEGIN

			bodies_found[tmp] = 1

			pp = HOSOrbit(tt,hos=2,degrees=degrees)
			pp = CvPrecess(tt,j2000,pp,/from_ecliptic,degrees=degrees)
			pp = cv_coord(from_sphere=pp,/to_rect,degrees=degrees)

			IF IsType(center,/defined) THEN pp -= jpl_eph(tt,center,/location)

			big_eph_boost, ntt, rr, pp, names, 'Helios B', nodist, 0B

		ENDIF

	ENDIF

	IF (where(bodies_found EQ 0))[0] NE -1 THEN BEGIN

		tmp = where_common(bodies,big_body('Stereo A',/alias))

		IF tmp[0] NE -1 THEN BEGIN

			bodies_found[tmp] = 1

			pp = StereoAOrbit(tt,degrees=degrees)
			pp = CvSky(j2000,from_ecliptic=pp, /to_equatorial, degrees=degrees, /silent)
			pp = cv_coord(from_sphere=pp, /to_rect, degrees=degrees)

			IF IsType(center,/defined) THEN pp -= jpl_eph(tt,center,/location)

			big_eph_boost, ntt, rr, pp, names, 'Stereo A', nodist, 0B

		ENDIF

	ENDIF

	IF (where(bodies_found EQ 0))[0] NE -1 THEN BEGIN

		tmp = where_common(bodies,big_body('Stereo B',/alias))

		IF tmp[0] NE -1 THEN BEGIN

			bodies_found[tmp] = 1

			pp = StereoBOrbit(tt,degrees=degrees)
			pp = CvSky(j2000,from_ecliptic=pp, /to_equatorial, degrees=degrees, /silent)
			pp = cv_coord(from_sphere=pp, /to_rect, degrees=degrees)

			IF IsType(center,/defined) THEN pp -= jpl_eph(tt,center,/location)

			big_eph_boost, ntt, rr, pp, names, 'Stereo B', nodist, 0B

		ENDIF

	ENDIF


	IF (where(bodies_found EQ 0))[0] NE -1 THEN BEGIN

		tmp = where_common(bodies,big_body('Ulysses',/alias))

		IF tmp[0] NE -1 THEN BEGIN

			bodies_found[tmp] = 1

			; UlyssesOrbit returns J2000 ecliptic coordinates.
			; We need J2000 equatorial.

			pp = UlyssesOrbit(tt,degrees=degrees)
			pp = CvSky(j2000,from_ecliptic=pp, /to_equatorial, degrees=degrees, /silent)
			pp = cv_coord(from_sphere=pp,/to_rect,degrees=degrees)

			IF IsType(center,/defined) THEN pp -= jpl_eph(tt,center,/location)

			big_eph_boost, ntt, rr, pp, names, 'Ulysses', nodist, 0B

		ENDIF

	ENDIF

	; The ephemerides return rect equat J2000 coord relative to center


	IF (where(bodies_found EQ 0))[0] NE -1 THEN BEGIN

		; Also check for Earth-orbiting spacecraft alongside JPL bodies.

		IF big_eph_short(sgp_alias(bodies),[jpl_body(),sgp_body()],names,add_list=list,found=bodies_found) NE 0 THEN BEGIN

			sc_index = where_common(list,sgp_body())
			IF sc_index[0] NE -1 THEN BEGIN 			; Earth-orbiting spacecraft present
				tmp = list[sc_index]					; Save names
				list[sc_index] = 'Earth'				; Set to Earth
			ENDIF

			pp = jpl_eph(tt,list,/location,center=center)

			IF sc_index[0] NE -1 THEN BEGIN
				pp = reform(pp,3,n_elements(list),ntt,/overwrite)
				list[sc_index] = tmp					; Restore names
				FOR i=0,n_elements(sc_index)-1 DO	$	; Add geocentric J2000 spacecraft location
					pp[*,sc_index[i],*] += sgp4_eph(tt,tmp[i],/location,/precess)
			ENDIF

			big_eph_boost, ntt, rr, pp, names, list, nodist, 0B

		ENDIF

	ENDIF

	IF (where(bodies_found EQ 0))[0] NE -1 THEN BEGIN

		IF big_eph_short(bodies,usno_body(),names,add_list=list,found=bodies_found) NE 0 THEN BEGIN

			pp = usno_eph(tt,list,/location,center=center,/silent)
			usno_close

			big_eph_boost, ntt, rr, pp, names, list, nodist, 0B

		ENDIF

	ENDIF

	IF (where(bodies_found EQ 0))[0] NE -1 THEN BEGIN

		IF big_eph_short(bodies,mpc_body(),names,add_list=list,found=bodies_found) NE 0 THEN BEGIN

			pp = mpc_eph(tt,list,/location,center=center,silent=2)

			big_eph_boost, ntt, rr, pp, names, list, nodist, 0B

		ENDIF

	ENDIF

	IF (where(bodies_found EQ 0))[0] NE -1 THEN BEGIN

		; Stars don't have a distance. Set it to 1 AU here, and set nodist=1

		bodies = strlowcase(smei_star_alias(bodies))
		n = big_eph_short(bodies,strtrim(smei_star_info(/get_name,/silent,cat='*')),names,add_list=list,found=bodies_found)

		IF n NE 0 THEN BEGIN

			pp = dblarr(3,n,/nozero)
			pp[0:1,*] = smei_star_info(name=list,degrees=degrees,silent=2,cat='*')
			pp[  2,*] = 1.0d0
			pp = SuperArray(pp,n_elements(tt),/trail)

			big_eph_boost, ntt, rr, pp, names, list, nodist, 1B

		ENDIF

	ENDIF

	IF (where(bodies_found EQ 0))[0] NE -1 THEN BEGIN

		IF big_eph_short(bodies,mpc_comets(),names,add_list=list,found=bodies_found) NE 0 THEN BEGIN

			pp = mpc_orbit_eph(tt,list,/location,center=center,silent=2,/comet)

			big_eph_boost, ntt, rr, pp, names, list, nodist, 0B

		ENDIF

	ENDIF

	IF (where(bodies_found EQ 0))[0] NE -1 THEN BEGIN

		IF big_eph_short(bodies,mpc_minor_planets(),names,add_list=list,found=bodies_found) NE 0 THEN BEGIN

			pp = mpc_orbit_eph(tt,list,/location,center=center,silent=2)

			big_eph_boost, ntt, rr, pp, names, list, nodist, 0B

		ENDIF

	ENDIF

ENDIF

; Removes bodies and times where no data are available.
; Note that this may modify the input tt array.

IF clean THEN big_eph_clean, tt, rr, names, nodist

; At this point all bodies with nodist=0 have rect equatorial J2000 coordinates
; For nodist=1 (stars) we now have spherical equatorial J2000 coordinates
; with the distance set to 1 AU.

count = n_elements(names)

CASE count EQ 0 OF
0: BEGIN

	InitVar, to_sphere	, /key
	InitVar, precess	, /key

	; Pick up planets, asteroids and minor planets (with nodist=0)

	pp = where(nodist EQ 0B, n)		; Planets, asteroids, minor planets

	IF n NE 0 THEN BEGIN
		IF IsType(sc_center,/defined) THEN BEGIN
			sc_loc = sgp4_eph(tt,sc_center,/location,/precess)
			rr[*,pp,*] -= SuperArray(sc_loc,n,after=1)
		ENDIF
	ENDIF

	to_rect = 1-to_sphere

	FOR i=0,count-1 DO BEGIN

		CASE nodist[i] OF
		0: BEGIN
			IF to_sphere THEN BEGIN
				rr[*,i,*] = cv_coord(from_rect=reform(rr[*,i,*]),/to_sphere,degrees=degrees)

				; Safety belt:
				; If rr values were bad, then cv_coord call sets the angles to zero, and only
				; then distance is bad. Set the angles bad here too.
				j = where(1-finite(rr[2,i,*]))
				IF j[0] NE -1 THEN rr[*,i,j] = BadValue(rr)
			ENDIF
		END
		1: IF to_rect THEN rr[*,i,*] = cv_coord(from_sphere=reform(rr[*,i,*]),/to_rect,degrees=degrees)
		ENDCASE

		; At this point we have equatorial J2000 coordinates (either rectangular or spherical)
		; Precess if necessary.

		IF precess THEN     $
			rr[0:2,i,*] = CvPrecess(j2000,tt,rr[0:2,i,*],degrees=degrees,rectangular=to_rect)

		; Now call CvSky to convert to the selected output coordinate frame
		; (the default is ecliptic coordinates)

		rr[0:2,i,*] = CvSky(tt,from_equator=rr[0:2,i,*],_extra=_extra,degrees=degrees,rectangular=to_rect,/silent)

		; Not sure anymore why this is necessary

		IF to_sphere THEN   $
			rr[0,i,*] = AngleRange(rr[0,i,*],degrees=degrees)

		IF count EQ 1 THEN BEGIN
			InitVar, onebody, /key
			IF onebody THEN rr = reform(rr,/overwrite)
		ENDIF

	ENDFOR

END

1: rr = replicate(BadValue(0.0d0),3)

ENDCASE

; Restore the input value of 'center'

IF IsType(sc_center,/defined) THEN center = sc_center

RETURN, rr  &  END