;+
; NAME:
;	smei_normal
; CALLING SEQUENCE:
	PRO smei_normal, _extra=_extra
; INCLUDE:
	@compile_opt.pro
;-

yr0  = 1900
doy0 = 0.5d0
t0   = TimeSet(yr=yr0,doy=doy0)

; ra0 = RA of fictitious sun on 2000 jan 1.0 in degrees

yr1  = 2000
doy1 = 1.5d0
t1   = TimeSet(yr=yr1,doy=doy1)

dt  = TimeOp(t1,t0,TimeUnit(/year))	; Julian years since t0
ra0 = [(18+24*(dt-floor(dt))) mod 24,38,45,836+5.42d0*dt+0.00929d0*dt*dt]
ra0 = AngleUnits(from_time=ra0, /to_degrees)
ra0 = ra0[0]

t_start = TimeSet('2003/03/01')
t_stop  = TimeSet('2011/03/01')

ndays = TimeOp(t_stop,t_start, TimeUnit(/day))+1
tt = TimeOp(/add,t_start, TimeSet(/diff,gridgen(ndays),TimeUnit(/day)))

; RA=18h38m45.836s + 8,640,184.542s TE + 0.0929s TE^2
; where TE is the number of Julian centuries (36525 days)
; since 1900 0.5 January.

nt  = n_elements(tt)
dra = fltarr(nt)
dec = fltarr(nt)

openw, /get_lun, iu,  'orbit_normal_data.txt'
FOR i=0,nt-1 DO BEGIN
	ra = ra_fictitious_sun(tt[i],/degrees,/scalar)
	pp = sgp4_orbit_axis(tt[i],/to_sphere,/degrees,/precess)

	dra[i] = AngleRange(pp[0]-ra,/pi,/degrees)
	dec[i] = pp[1]

	printf, iu, format='(A,3F12.4)',	$
		TimeGet(tt[i], /ymd,/scalar, upto=TimeUnit(/day)),	$
		ra, dra[i], dec[i]
ENDFOR
close, iu

whatis, dra, dra
plotcurve, tt, dra,	$
	/ynozero	, $
	ytitle	 = "RA-RA(fictitious Sun)"	, $
	charsize = 1.5
	_extra	 = _extra

write_png, 'orbit_normal_plot.png', tvrd(/true)

return

yr0 = 2007
doy0 = 83.5
t0 = timeset(yr=yr0,doy=doy0)

yr = 2007
doy = 83.5
nstep = 6
tt = timeset(yr=yr,doy=doy+gridgen(401,range=[-150,150]*nstep))

jd0 = timeget(t0,/njd,/scalar)
jd  = timeget(tt,/njd)

ra = 360.0/365.25*(jd-jd0)

pp = sgp4_orbit_axis( tt, /to_sphere, /degrees, /precess )
print, pp
plotcurve, tt, AngleRange(pp[0,*]-ra,/pi,/degrees),	$
	/ynozero	, $
	ytitle	= "RA-360/365.25*(t-[2007,doy 083.5])"	, $
	_extra	= _extra

RETURN  &  END