;+
; NAME:
;	getootyasources
; PURPOSE:
;	Extract Ooty IPS observations from yearly data files
; CATEGORY:
;	I/O
; CALLING SEQUENCE:
	FUNCTION getootysources, 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/ooty 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 IPS sources inside time range
;		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.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=[1999,2002], 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)
yr_range = TimeGet(trange, TimeUnit(/year))
IF use_good_sources NE '' THEN yr_range[1] = yr_range[0]

CASE IsType(subdir,/defined) OF
0: BEGIN
	InitVar, subdir, 'daily'
	ipsdir = filepath(root=getenv('SSW_SMEI_DAT'),subdir='ooty',subdir)
END
1: BEGIN
	CASE CheckDir(subdir) OF
	0: BEGIN
		InitVar, subdir, 'daily'
		ipsdir = filepath(root=getenv('SSW_SMEI_DAT'),subdir='ooty',subdir)
	END
	1: ipsdir = subdir
	ENDCASE
END
ENDCASE

nsources = 0L
destroyvar, point_sources

FOR yr=yr_range[0],yr_range[1] DO BEGIN				; Loop over yearly data files

	file = use_good_sources NE '' ? 'ooty.'+use_good_sources+'.good' : 'ooty.'+strcompress(yr,/rem)
	file = filepath(root=ipsdir,file)		; Construct file name

	status = txt_read(file, data, silent=silent)	; Read yearly IPS file

	IF status THEN BEGIN

		srcs = strmid(data,18,8)

		yrs  = TimeFixYear( fix( strmid(data,2,2) ) )
		mons = fix  ( strmid(data,  4,2) )
		days = fix  ( strmid(data,  6,2) )
		uts  = float( strmid(data, 10,6) )

		vvs  = float( strmid(data, 79, 6) )
		ggs  = float( strmid(data, 91, 7) )
		pps  =  	  strmid(data, 47,18)

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

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

		tmp = where(ggs LE -99.0)
		IF tmp[0] NE -1 THEN ggs[tmp] = BadValue(ggs)

		tts = TimeSet(yr=yrs, mon=mons, day=days, hour=uts)	; Time of obs
		jds = TimeGet(tts, /njd)							; Time of obs in Julian days

		; 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(vvs) OR finite(ggs)), ni)

		IF ni NE 0 THEN BEGIN								; If obs inside time range found
			pps = reform(flt_string(pps[isrc]),3,2,ni)
			pps = AngleUnits(from_almanac=pps, /to_degrees)

			i = vu_point_source(i, srcs[isrc], /name, create=ni)
			i = vu_point_source(i, pps		 , /location, /degrees)
			i = vu_point_source(i, tts[isrc] , /time)

			i = vu_point_source(i, vvs[isrc] , /ips_speed		)
			i = vu_point_source(i, ggs[isrc] , /ips_glevel		)

			CASE yr EQ yr_range[0] OF
			0: point_sources = [point_sources, i]
			1: point_sources = i
			ENDCASE

			nsources += ni

		ENDIF

	ENDIF

ENDFOR

; Precess from B1950 to current epoch

IF IsType(point_sources, /struct) THEN BEGIN
	tmp = CvPrecess(TimeSet(bepoch=1950D0), 		$
		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])

say, strcompress(nsources) +' IPS sources in '+	$
	 TimeGet(TT,/ymd,upto=TimeUnit(/minute))+' [+/-'+strcompress(dT)+' hr]', threshold=1, silent=silent

RETURN, nsources  &  END