;+
; NAME:
;	sgp4_orbit_axis
; PURPOSE:
;	Gets ECI coordinates for Coriolis orbital axis (i.e. the normal
;	to the orbital plane, close to the sunward direction).
; CATEGORY:
;	camera/idl/toolbox
; CALLING SEQUENCE:
	FUNCTION sgp4_orbit_axis, tt, $
		degrees		= degrees	, $
		to_sphere	= to_sphere	, $
		precess		= precess
; INPUTS:
;	tt		array; type: time structure
; OPTIONAL INPUT PARAMETERS:
;	/to_sphere	if set return spherical (RA,dec) coordinates
;	/degrees	if set, angles are in degrees (default: radians)
;	/precess	(passed to sgp4_eph)
;				precess to J2000 coordinates
; OUTPUTS:
;	rr		array[3,*]; type: double
;				coordinates of a unit vector along the normal
;				to the SMEI orbital plane in rectangular (/to_sphere
;				NOT set) or spherical (RA,dec) coordinates
;				(/to_sphere) SET. 
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	TimeSystem, TimeUnit, TimeSet, gridgen, sgp4_eph, vectorproduct
;	InitVar, IsTime, TimeOp
; PROCEDURE:
;	Viewed from the Sun SMEI circles the Earth close to the terminator
;	in anti-clockwise direction (going north towards the east (dusk
;	terminator) and south in the west (dawn terminator).
;	The normal, defined as the cross produce of position and velocity
;	vector points towards the Sun.
; MODIFICATION HISTORY:
;	APR-2009, Paul Hick (UCSD/CASS)
;	MAR-2011, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Added more documentation.
;-

InitVar, to_sphere, /key
InitVar, tt, TimeSystem(/silent)
umin = TimeUnit(/min)

CASE 1 OF
IsTime(tt		): uu = tt
IsTime(tt,/orbit): uu = TimeSet(smei=tt)
ELSE			 : uu = TimeSet(tt)
ENDCASE

sztt = size(uu,/dim)
nu = n_elements(uu)
uu = reform(uu,nu,/overwrite)

; Set up 5 points spread out over one orbit.
; These are averaged to get the final result

nt = 5
dt = TimeSet(/diff,gridgen(nt,range=[-50,50]),umin)
uu = replicate(!TheTime,nu,nt)
FOR i=0,n_elements(tt)-1 DO uu[i,*] = TimeOp(/add,tt[i],dt)

uu = reform(uu,nu*nt,/overwrite)
rr = sgp4_eph(uu,/km,precess=precess)

vv = rr[3:5,*]							; Velocity vector
rr = rr[0:2,*]							; Position vector

rr /= [1,1,1]#sqrt(total(rr*rr,1))		; Convert to unit vectors
vv /= [1,1,1]#sqrt(total(vv*vv,1))

; Normal to orbital plane (points towards Sun)

rr = vectorproduct(rr,vv)

; Average over all nt points along one orbit

rr = reform(rr,[3,nu,nt])
rr = total(rr,3)/nt
rr /= [1,1,1]#sqrt(total(rr*rr,1))		; Convert to unit vectors

IF to_sphere THEN BEGIN
	rr = cv_coord(from_rect=rr, /to_sphere, degrees=degrees)
	rr[0,*] = AngleRange(rr[0,*], degrees=degrees)
ENDIF

IF nu NE 1 THEN rr = reform(rr,[3,sztt],/overwrite)

RETURN, rr  &  END