;+
; NAME:
;	smei_sgp4_orbits
; PURPOSE:
;	Build list of start times for all orbits using TLE data base.
; CALLING SEQUENCE:
	PRO smei_sgp4_orbits, input_sgp4_file, new=new, smeitree=smeitree, silent=silent
; INCLUDE:
	@compile_opt.pro	; On error, return to caller
; INPUTS:
;	input_sgp4_file
;				existing file with start times.
;				By default the file in the SMEI tree is used
;				(if /new is NOT set)
; OPTIONAL INPUTS:
;	/new		creates a new file from scratch. In this case
;				input_sgp4_file is ignored
;	/smeitree 	The output file is created as $TUB/smei_sgp4_orbits.tmp
;				If /smeitree is set then the file in $TUB is copied to
;				the SMEI master tree and to the SMEI tree on smei.ucsd.edu
;				(requires that env $SMEIMASTER is defined)
; CALLS:
;	InitVar, txt_read, TimeUnit, TimeSet, TimeGet, TimeOp
;	sgp4_eph, who_am_i, sgp4_orbit_axis, big_eph, sphere_distance
; PROCEDURE:
;	sgp4_eph returns rectangular ECI coordinates relative to the
;	current equinox. Essentially this means equatorial coordinates
;	relative to the current equinox. So the declination is
;	effectively the geographic latitude.
;
;	Currently (March 2011) this procedure is called as part of
;	smei_pipeline.sh (executed hourly as cronjob on smeidb).
;
;	The file is used in IDL only by the procedure smei_coriolis.pro.
;	In Fortran the only function calling it is smei_orbit_info2.f.
;
;	Content of file:
;
;	1. SMEI orbit number
;	2. Orbit start time in format YYYY_DOY_hhmmssmss
;	   The orbital start time in format hh:mm:ss.mss
;	   Defined as the time Coriolis crosses a fixed geographic
;	   latitude (~ -31.1 degrees; close to the latitude of the
;	   SAA) going north.
;	3. Orbital period, defined as time difference between two
;	   successive orbital start times
;	4. Orbital period used in sgp4 software
;	   I think this is some sort of Keplerian orbital period;
;	   not particulary useful.
;	5. Precision of geographic latitude calculation
;
;	The remaining columns are related to the normal of the SMEI orbital
;	plane (pointing towards the Sun).
;
;	6. Solar elongation (in degrees) of the normal
;	7. RA (in degrees) of the fictitious Sun
;	   This progresses stricly linear with UT going through
;	   360 degrees in 365.25 days (one Julian year).
;	   This is approximate, but good enough for government work)
;	8. RA difference (J2000, in degrees) of the normal with the
;	   RA of the fictitious Sun.
;	9. Declination of the normal (J2000, in degrees)
; MODIFICATION HISTORY:
;	JAN-2006, Paul Hick (UCSD/CASS)
;	OCT-2006, Paul Hick (UCSD/CASS)
;		Replaced NewcombSun by big_eph calls
;	MAR-2007, Paul Hick (UCSD/CASS)
;		Eliminated big_eph. Now only sgp4_eph is used.
;	APR-2007, Paul Hick (UCSD/CASS)
;		Added check to avoid calculating start times past
;		the time of the last TLE available to sgp4_eph.
;	APR-2009, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Added column for elongation of normal to Coriolis
;		of orbital plane.
;	MAR-2011, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Added columns for RA,dec of orbit normal
;-

InitVar, new   , /key
InitVar, silent, 0

smei_sgp4_name = who_am_i(/procedure)	; Name of this procedure: smei_sgp4_orbits

usec  = TimeUnit(/sec )
umsec = TimeUnit(/msec)
uday  = TimeUnit(/day )

; smei_sgp4_file is the file with orbit numbers in the SMEI tree
; (this is the copy used by all IDL SMEI routine)

smei_sgp4_file = filepath(root=who_am_i(/dir),subdir='smei',smei_sgp4_name+'.txt')

; Temporary local copy of smei_sgp4_orbits. Once complete this file
; gets copied to the proper locations.

tmp_sgp4_file = filepath(root=getenv('TUB'),smei_sgp4_name+'.tmp')

IF NOT new THEN BEGIN
	InitVar, input_sgp4_file, smei_sgp4_file
	new = 1-txt_read(input_sgp4_file,content)
ENDIF

; Check the last TLE available. Extract the time for the last
; TLE. We do not want to go past this time in the calculation
; orbit start times. We only want start times that are bracketed
; by TLE; we DO NOT want to extrapolate.

smei_tle_file = filepath(root=who_am_i(/dir),subdir='sgp','sat27640.txt')
IF NOT txt_read(smei_tle_file, tle) THEN $
	message, 'no TLEs found: '+hide_env(smei_tle_file)

last_tle = (where( strpos(tle,'<End of file>') EQ 0 ))[0]
IF last_tle EQ -1 THEN	$
	message, 'TLE file should have <End of file> marker: '+hide_env(smei_tle_file)
last_tle -= 2

tle = tle[last_tle]

; Time of last TLE

tle = TimeSet(yr=2000+long(strmid(tle,18,2)),doy=double(strmid(tle,20,12)))

CASE new OF
0: BEGIN
	i = n_elements(content)

	tt_old = content[i-2]  				; One-but-last record
	tt_old = TimeSet(strmid(tt_old,8,18)); One-but-last orbit start time

	tmp = content[i-1]  				; Last record

	norbit = long	 (strmid(tmp, 0, 6)); Last orbit number
	tt		= TimeSet(strmid(tmp, 8,18)); Last orbit start time
	period1 = TimeSet(strmid(tmp,28,12)); Last orbital period
	dgeolat = double (strmid(tmp,54,12)); Last geolat uncertainty

	elo_orbnorm = double (strmid(tmp,66,10)); Last elongation of orbit normal
	ra_sun		= double (strmid(tmp,76,10) )
	dra_orbnorm = double (strmid(tmp,86,10) )
	dec_orbnorm = double (strmid(tmp,96,10) )

	content = content[0:i-2]			; Drop last record
END
1: BEGIN
	norbit = 286L						; Orbit near max orbital period
	tt = TimeSet('2003_021_023458127')

	period1 = !TheTime
	dgeolat = 0.0d0

	; Calculate the angle from the orbit normal to the Sun

	rsun = big_eph(tt,		  $				; Direction to sun (J2000)
		body	= 'sun'		, $
		center	= 'earth'	, $
		/to_sphere			, $
		/onebody			, $
		/to_equatorial		, $
		/degrees			, $
		/location			)
											; Direction orbit normal (J2000)
	rnorm = sgp4_orbit_axis(tt,/degrees,/to_sphere,/precess)
											; Elongation orbit normal
	elo_orbnorm = sphere_distance(rsun,rnorm,/degrees)

	ra_sun		= ra_fictitious_sun(tt,/degrees)
	dra_orbnorm	= AngleRange(rnorm[0]-ra_sun,/degrees,/pi)
	dec_orbnorm	= rnorm[1]

	destroyvar, content
END
ENDCASE

norbit_max = 50000L

geolat_ref = sgp4_eph(tt,'smei',/location)
geolat_ref = cv_coord(from_rect=geolat_ref,/to_sphere,/degrees)
geolat_ref = geolat_ref[1]				; Geographic latitude = declination
print, 'Reference geographic latitude is ', geolat_ref

dsec  = 10*[-1,1]						; +/- 10 second time window
dmsec = -20+indgen(41)					; +/- 20 msecond time window

IF silent LE 0 THEN message, /info, 'write '+hide_env(tmp_sgp4_file)
openw, /get_lun, iu, tmp_sgp4_file

; Rewrite everything but the last record.
; For a new file this means nothing is written.

FOR i=0L,n_elements(content)-1 DO printf, iu, content[i]

REPEAT BEGIN

	; Break if orbit start time is past time of last
	; SMEI TLE

	IF TimeOp(/subtract,tt,tle,usec) GE 0 THEN BEGIN
		IF silent LE 0 THEN message, /info, 'passed time of last TLE @ '+TimeGet(/_ydoy,tle,upto=usec)
		break
	ENDIF

	; Not sure exactly what period2 represents, but it is
	; very close to 7 seconds shorter than the orbital period
	; we want (time between crossings of given geographic
	; latitude).

	period2 = sgp4_orbit_period(tt)
	period2 = TimeSet(/diff,period2,uday)

	num_format  = '(I6,6A,E12.3,4F10.4)'
	time_format = 'hh:mm:ss.mss'

	printf, iu, format=num_format		, $
		norbit							, $
		'  ',TimeGet(tt,/_ydoy)			, $
		'  ',TimeGet(period1,format=time_format),	$
		'  ',TimeGet(period2,format=time_format),	$
		dgeolat							, $
		elo_orbnorm						, $
		ra_sun							, $
		dra_orbnorm						, $
		dec_orbnorm

	IF silent LE 0 THEN	$
		print, format=num_format		,$
			norbit						, $
			'  ',TimeGet(tt,/_ydoy)		, $
			'  ',TimeGet(period1,format=time_format),	$
			'  ',TimeGet(period2,format=time_format),	$
			dgeolat						, $
			elo_orbnorm					, $
			ra_sun						, $
			dra_orbnorm					, $
			dec_orbnorm

	tt_old = tt 							; Save last start time

	; First estimate of next start time is one period later
	; (but remember period2 underestimates by about 7 seconds.

	tt = TimeOp(/add,tt,period2)

	; Get s/c locations at times bracketing estimated start time
	; by +/-10 seconds. This should bracket the start time we want.

	ts = TimeOp(/add,tt,TimeSet(/diff,dsec,usec))

	geolat = sgp4_eph(ts,'smei',/location)
	geolat = cv_coord(from_rect=geolat, /to_sphere,/degrees)
	geolat = reform(geolat[1,*])

	; Interpolate on bracketing times to find time where lat is
	; equal to the ref latitude. tt is within milliseconds of
	; the time we want.

	dt = interpol(dsec,geolat,geolat_ref)
	tt = TimeOp(/add,tt,TimeSet(/diff,dt,usec))

	; Set up times around best estimate at msec spacing
	; for time window of +/-20 msec.

	tt = TimeOp(/add,tt,TimeSet(/diff,dmsec,umsec))

	geolat = sgp4_eph(tt,'smei',/location)
	geolat = cv_coord(from_rect=geolat, /to_sphere,/degrees)
	geolat = reform(geolat[1,*])

	; Pick the time where the spacecraft latitude is closest to
	; the reference latitude.
	; This defines the orbit start time with a precision of 0.5 msec

	dgeolat = min(abs(geolat-geolat_ref),ip)
	dgeolat = geolat[ip]-geolat_ref

	tt = tt[ip] 							; Orbit start time
	period1 = TimeOp(/subtract,tt,tt_old)	; Orbital period

	; Calculate the angle from the orbit normal to the Sun

	rsun = big_eph(tt,		  $				; Direction to sun (J2000)
		body	= 'sun'		, $
		center	= 'earth'	, $
		/to_sphere			, $
		/onebody			, $
		/to_equatorial		, $
		/degrees			, $
		/location			)
											; Direction orbit normal (J2000)
	rnorm = sgp4_orbit_axis(tt,/degrees,/to_sphere,/precess)
											; Elongation orbit normal
	elo_orbnorm = sphere_distance(rsun,rnorm,/degrees)

	ra_sun		= ra_fictitious_sun(tt,/degrees)
	dra_orbnorm	= AngleRange(rnorm[0]-ra_sun,/degrees,/pi)
	dec_orbnorm	= rnorm[1]

	norbit += 1 							; Increment orbit counter

ENDREP UNTIL norbit GT norbit_max

free_lun, iu

; Temp copy of sgp4 file is tmp_sgp4_file (in $TUB).
; If /smeitree is set, copies of this file are put in $SMEI,
; both on ips.ucsd.edu and smei.ucsd.edu 

InitVar, smeitree, /key
smei	    = getenv('SMEI')
smei_master = getenv('SMEI_MASTER')

IF smeitree AND smei NE '' AND smei_master NE '' THEN BEGIN

	hostname = getenv('HOSTNAME')

	; Copy to SMEI master tree. Use cp if we are running on the SMEI master host.
	; If not use scp (needs public key authentication)

	CASE strpos(smei_master,hostname) NE -1 OF
	0: tmp = 'scp -q '+tmp_sgp4_file+' '+smei_master+strmid(smei_sgp4_file,strlen(smei))
	1: tmp = 'cp -v'  +tmp_sgp4_file+' '+smei_sgp4_file
	ENDCASE 

	IF silent LE 1 THEN print, tmp
	spawn, tmp

	CASE hostname EQ 'smei.ucsd.edu' OF
	0: tmp = 'scp -q '+tmp_sgp4_file+' soft@smei:'+smei_sgp4_file
	1: tmp = 'cp -v ' +tmp_sgp4_file+' '+smei_sgp4_file
	ENDCASE

	IF silent LE 1 THEN print, tmp
	spawn, tmp

	spawn, 'rm '+tmp_sgp4_file

ENDIF

RETURN  &  END