;+
; NAME:
;	smei_coriolis
; PURPOSE:
;	Returns time at which Coriolis was in specified orbit
; CALLING SEQUENCE:
	FUNCTION smei_coriolis, tt_or_orbit , $
		orbital_period	= orbital_period, $
		number			= number		, $
		fraction		= fraction		, $
		doublex 		= doublex
; INPUTS:
;	tt_or_orbit			array; type: time structure or any numerical type
;							if input are times then output is the 
;							orbit number; if input is numerical then
;							output is the orbit time.
; OPTIONAL INPUTS:
;	/orbital_period 	return orbital period
;							if set then the orbital period at the specified
;							time or orbit number is returned
;							(if tt_or_orbit is not specified then
;							tt_or_orbit=0 is assumed)
;	/number 			these three keywords are used only if a time
;	/fraction			structure is specified as input (the keywords
;	/doublex			are passed to href=smei_orbit_get=.
; OUTPUTS:
;	smei_coriolis		array; type: numerical or time structure
;							orbit number or orbit time (depending on input
;							tt_or_orbit
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar, TimeSet, IsTime, TimeOp, TimeGet, timeposn, IsType
; PROCEDURE:
;	Orbit n0(=0) starts at t0 = 2002/12/31 23:56:41.900.
;	The orbital period at this time was p0=0.0705612268519d0 days (6096.490 s),
;	and has been decreasing since with dp=0.084 msec per orbit.
;	(0.42s reduction in SMEI period per 365 days).
;
;	The effective orbital period for orbit n is p0-(n-n0)*dp.
;	The start times for each orbit becomes:
;
;	t(n0+0) = t0
;	t(n0+1) = t(n0+0)+p0-1*dp
;	t(n0+2) = t(n0+1)+p0-2*dp
;	...
;	t(n0+n) = t(n0+n-1)+p0-n*dp = t0+n*p0-0.5*n*(n+1)*dp
; MODIFICATION HISTORY:
;	JUN-2004, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

InitVar, orbital_period , /key

new = getenv('ORBIT_DEF') EQ ''

CASE new OF

0: BEGIN

	n0 = 1.0d0

	; Orbit zero starts at 2002/12/31 23:56:41.900
	; (86201900 = 86400000-198100)

	t0 = TimeSet(yr=2003, doy=0, msec=86201900)
	p0 = 0.0705612268519d0				; Nominal orbital period in days
	dp = 0.084d0						; Correction per orbit in msec

	InitVar, tt_or_orbit, n0

	CASE 1 OF
	IsType(tt_or_orbit,/string) : ttorb = timeposn(tt_or_orbit,/extract,part='name')
	IsTime(tt_or_orbit) 		: ttorb = tt_or_orbit
	IsType(tt_or_orbit,/struct) : ttorb = smei_orbit_get(tt_or_orbit,/number)+smei_orbit_get(tt_or_orbit,/fraction)
	ELSE						: ttorb = tt_or_orbit
	ENDCASE

	CASE 1 OF

	orbital_period: begin

		CASE IsTime(ttorb) OF
		0: rtn = ttorb					; Orbit number
		1: rtn = smei_coriolis(ttorb)	; Convert time to orbit number
		ENDCASE

		rtn = TimeSet(/diff, day=p0-((rtn-n0)*dp)/86400000.0d0)

	END

	IsTime(ttorb): BEGIN							; Time structure input

		tt = TimeOp(ttorb, t0, /subtract) 			; Time difference 

		; First guess at orbit number using the nominal orbital period
		; This will underestimate the number of orbits elapsed since
		; the orbital period decreases over time

		dt = TimeGet(tt, /day, /diff, /full)
		orbit = n0+dt/TimeGet(smei_coriolis(/orbital_period),/day,/diff,/full,/scalar); Underestimate # orbits elapsed
		orbit = n0+dt/TimeGet(smei_coriolis((n0+orbit)/2,/orbital_period),/day,/diff,/full)

		; The time tt at which the underestimated number of orbits
		; are elapsed is earlier than the input time ttorb

		tt = smei_coriolis(orbit)					; Exact time for 'orbit'
		tt = TimeOp(ttorb,tt,/subtract)				; Difference in days
		tt = TimeGet(tt, /day, /diff, /full)		; Positive number

		; Correct the approximate number of elapsed orbits
		; using the instantaneous orbital period

		rtn = orbit+tt/TimeGet(smei_coriolis(orbit,/orbital_period),/day,/diff,/full)
		rtn = smei_orbit_set(rtn)

	END

	ELSE: BEGIN 									; Numerical input: must be orbit number

		n  = ttorb-n0
		dt =  n*p0									; Days since t0
		ds = -n*(n+1.0d0)*0.5d0*dp					; Additional milli-seconds since t0

		rtn = TimeOp(t0, TimeSet(day=dt, msec=ds, /diff), /add)

	END

	ENDCASE

END

1: BEGIN

	common smei_coriolis_save, orbit0, norbit, torbit

	IF IsType(orbit0,/undefined) THEN BEGIN
		IF NOT txt_read(filepath(root=who_am_i(/dir),subdir='smei','smei_sgp4_orbits.txt'),torbit) THEN	$
			message, 'create new database with smei_sgp4_orbits.pro'

		orbit0 = long(strmid(torbit[0],0,6));'   286  2003_021_023458127  01:41:36.484'
		torbit = TimeSet(strmid(torbit,8,18))
		norbit = n_elements(torbit)
	ENDIF

	InitVar, tt_or_orbit, orbit0

	; For string input try to extract a time

	CASE IsType(tt_or_orbit, /string) OF
	0: ttorb = tt_or_orbit
	1: ttorb = timeposn(tt_or_orbit,/extract,part='name')
	ENDCASE

	CASE 1 OF

	orbital_period: begin

		CASE 1 OF
		IsTime(ttorb)		 : rtn = smei_coriolis(ttorb)
		IsType(ttorb,/struct): rtn = ttorb
		ELSE		  		 : rtn = smei_orbit_set(ttorb)
		ENDCASE

		n   = smei_orbit_get(rtn,/number)-orbit0
		rtn = TimeOp(/subtract,torbit[n+1],torbit[n])

	END

	IsTime(ttorb): BEGIN				; Time structure input

		rtn = TimeLInterpol(nothing, torbit, ttorb, ilow=ilow, fraction=fraction_) 
		rtn = smei_orbit_set(orbit0+ilow,fraction_)

		InitVar, number  , /key
		InitVar, fraction, /key
		InitVar, doublex , /key

		CASE 1 OF
		number	: rtn = smei_orbit_get(rtn, /number  )
		fraction: rtn = smei_orbit_get(rtn, /fraction)
		doublex : rtn = smei_orbit_get(rtn, /doublex )
		ELSE	: ; Nothing to do
		ENDCASE

	END

	ELSE: BEGIN 						; Numerical input: must be orbit number

		CASE IsType(ttorb,/structure) OF
		0: BEGIN
			wu = ttorb-orbit0
			n  = long(wu) > 0 < (norbit-2)
			wu = wu-n
		END
		1: BEGIN
			wu = smei_orbit_get(ttorb,/number)-orbit0
			n  = wu > 0 < (norbit-2)
			wu = smei_orbit_get(ttorb,/fraction)+(wu-n)
		END
		ENDCASE

		rtn = TimeOp(/meant,torbit[n],torbit[n+1],wu=wu)

	END

	ENDCASE

END

ENDCASE

RETURN, rtn  &  END