;+
; NAME:
;	smei_sgp4_orbits
; PURPOSE:
;	Build list of start times for all orbits using TLE data base.
; CALLING SEQUENCE:
	PRO smei_sgp4_orbits, file, new=new, remote=remote
; INCLUDE:
	@compile_opt.pro	; On error, return to caller
; INPUTS:
;	file		existing file with start times
;				(by default the file in the SMEI tree is used)
; OPTIONAL INPUTS:
;	/new		creates a new file from scratch
;				(by default the existing file in the SMEI tree
;				is updated)
;	/remote 	copies the modified file with start times
;				to the SMEI master tree
; CALLS:
;	InitVar, txt_read, TimeUnit, TimeSet, TimeGet, TimeOp
;	sgp4_eph, who_am_i
; 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.
; 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; pphick@ucsd.edu)
;		Added check to avoid calculating start times past
;		the time of the last TLE available to sgp4_eph.
;-

InitVar, new, /key

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

IF NOT new THEN BEGIN
	InitVar, file, filepath(root=who_am_i(/dir),subdir='smei','smei_sgp4_orbits.txt')
	new = 1-txt_read(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.

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

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

tlm = tlm[last_tlm]

; Time of last TLE

tlm = TimeSet(yr=2000+long(strmid(tlm,18,2)),doy=double(strmid(tlm,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

	content = content[0:i-2]			; Drop last record
END
1: BEGIN
	norbit = 286L						; Orbit near max orbital period
	tt = TimeSet('2003_021_023458127')
	InitVar, file, filepath(root=getenv('TUB'),'smei_sgp4_orbits.txt')

	period1 = !TheTime
	dgeolat = 0.0d0

	destroyvar, content
END
ENDCASE

norbit_max = 40000L

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

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

message, /info, 'write to '+hide_env(file)
openw, /get_lun, iu, 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,tlm,usec) GE 0 THEN BEGIN
		message, /info, 'passed time of last TLE @ '+TimeGet(/_ydoy,tlm,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)

	printf, iu, format='(I6,6A,E12.3)',	$
		norbit,							$
		'  ',TimeGet(tt,/_ydoy),		$
		'  ',TimeGet(period1,format='hh:mm:ss.mss'),	$
		'  ',TimeGet(period2,format='hh:mm:ss.mss'),	$
		dgeolat

	print, format='(I6,6A,E12.3)',		$
		norbit,							$
		'  ',TimeGet(tt,/_ydoy),		$
		'  ',TimeGet(period1,format='hh:mm:ss.mss'),	$
		'  ',TimeGet(period2,format='hh:mm:ss.mss'),	$
		dgeolat

	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

	norbit += 1 							; Increment orbit counter

ENDREP UNTIL norbit GT norbit_max

free_lun, iu

InitVar, remote, /key
smei		= getenv('SMEI')
remote_smei = getenv('SMEI_MASTER')

remote = remote AND smei NE '' and remote_smei NE '' AND strpos(file,smei) EQ 0

IF remote THEN BEGIN
	
	tmp = 'scp '+file+' '+remote_smei+strmid(file,strlen(smei))
	print, tmp
	spawn, tmp

ENDIF

RETURN  &  END