;+
; NAME:
;	getsmeisources
; PURPOSE:
;	Extract SMEI lines of sight from orbital data files
; CATEGORY:
;	I/O
; CALLING SEQUENCE:
	FUNCTION getsmeisources, trange, $
		subdir				= subdir			, $
		use_good_sources	= use_good_sources	, $
		degrees 			= degrees			, $
		silent				= silent			, $
		point_sources 		= point_sources
; INPUTS:
;	trange			scalar or array[2]; type: time structure
;						time range (begin and end time)
;						a scalar is interpreted as [trange, trange+1 day]
;						(i.e. a one day time range starting at trange
; OPTIONAL INPUT PARAMETERS:
;	subdir=subdir	scalar; type: string; default='daily'
;						subdirectory in $DAT/smei where to look for data files
;	/degrees		if set pp will be in degrees (default: radians)
; OUTPUTS:
;	Result			scalar; type: long;
;						# ips sources in time range
; OPTIONAL OUTPUT PARAMETERS:
;	point_sources=point_sources
;			array[n]; type: structure
;				structure with lines of sight inside time range
;				(see href=vu_point_source=)
;		point_sources.name = string with source name
;		point_sources.pp   = 2-element array with RA and decl of source
;		point_sources.tt   = time structure with time of observation
;		point_sources.vv   = IPS velocity
;		point_sources.gg   = IPS g-level
;		point_sources
;		point_sources.xy   = used by PlotEarthSkymap to set pixel locations of sources
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	TimeFixYear, txt_read, TimeOp, TimeSet, TimeGet, TimeUnit
;	CvSky, NewcombSun, elongation, BadValue
;	CvPrecess, IsType, InitVar, vu_point_source
; PROCEDURE:
; MODIFICATION HISTORY:
;	AUG-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Derived from getnagoyasources
;-

InitVar, use_good_sources, ''
InitVar, silent , 0

dpm = ToDegrees(degrees=degrees)

CASE n_elements(trange) OF
0: trange = TimeSet(yr=[2003,2008], doy=1)
1: trange = [trange, TimeOp(/add,trange,TimeSet(/diff, day=1))]
2: IF IsType(trange,/string) THEN trange = TimeSet(trange)
ELSE: trange = [trange[0], trange[n_elements(trange)-1]]
ENDCASE

jd_range  = TimeGet(trange, /njd)
orb_range = smei_orbit_get( TimeGet(trange, /smei), /floor )
IF use_good_sources NE '' THEN orb_range[1] = orb_range[0]

CASE IsType(subdir,/defined) OF
0: BEGIN
	InitVar, subdir, 'points'
	smeidir = filepath(root=getenv('SSW_SMEI_DAT'),subdir='smei',subdir)
END
1: BEGIN
	CASE CheckDir(subdir) OF
	0: BEGIN
		InitVar, subdir, 'points'
		smeidir = filepath(root=getenv('SSW_SMEI_DAT'),subdir='smei',subdir)
	END
	1: smeidir = subdir
	ENDCASE
END
ENDCASE

nsources = 0L
destroyvar, point_sources

FOR orb=orb_range[0],orb_range[1] DO BEGIN			; Loop over yearly data files
													; Construct file name
	file = use_good_sources NE '' ? 'points.'+use_good_sources+'.txt' : 'points_*_'+string(orb,format='(I5.5)')+'.txt'
	file = filepath(root=smeidir,file)

	status = flt_read(file, data, silent=silent, /double); Read SMEI orbit file

	IF status THEN BEGIN

		tts	= TimeSet(smei=reform(data[3,*]))
		jds = TimeGet(tts, /njd)					; Time of obs in Julian days

		pps =  	     data[1:2,*]
		bbs = reform(data[  4,*])

		; Check for 'bad data'. For velocities check for value -999 or
		; quality indicator -999. For g-level check for value 0.0

		tmp = where(bbs LE -9999.0d0)
		IF tmp[0] NE -1 THEN bbs[tmp] = BadValue(bbs)


		; Pick up all sources inside time range with either V or g good

		isrc = where(jd_range[0] LE jds AND jds LE jd_range[1] AND finite(bbs), ni)

		IF ni NE 0 THEN BEGIN								; If obs inside time range found

			i = vu_point_source(i, lindgen(ni)	, /name, create=ni)
			i = vu_point_source(i, pps[*,isrc]	, /location, /degrees)
			i = vu_point_source(i, tts[isrc]	, /time)
			i = vu_point_source(i, bbs[isrc]/0.552, /brightness)

			CASE IsType(point_sources, /undefined) OF
			0: point_sources = [point_sources, i]
			1: point_sources = i
			ENDCASE

			nsources += ni

		ENDIF

	ENDIF

ENDFOR

; Precess from J2000 to current epoch

IF IsType(point_sources, /struct) THEN BEGIN
	tmp = CvPrecess(TimeSet(jepoch=2000.0d0), 		$
		vu_point_source(point_sources,/time),		$
		vu_point_source(point_sources,/location,degrees=degrees),	$
		degrees=degrees)
	point_sources =	vu_point_source(point_sources, tmp, /location, degrees=degrees)
ENDIF

dT = TimeOp(/subtract,trange[1],trange[0],TimeUnit(/hour))
IF dT/2*2 EQ dT THEN dT /= 2 ELSE dT /= 2.0
TT = TimeOp(/meant,trange[0],trange[1])

IF silent LE 1 THEN $
	message, /info, strcompress(nsources) +' SMEI lines of sight in '+	$
		TimeGet(TT,/ymd,upto=TimeUnit(/minute))+' [+/-'+strcompress(dT)+' hr]'

RETURN, nsources  &  END