;; This code computes the times of Vernal Equinox from the period 1951
;; to 2049.  It is written in the IDL language, and requires several
;; programs from the Markwardt IDL library.  It also requires the JPL
;; planetary ephemeris.

;; The times of Equinox are computed in Julian Days (UTC) and placed
;; in the JDVERN variable, and saved to a file called
;; "vernequinox.dat".  The differences between successive equinoxes
;; are also saved.

;; Load the JPL ephemeris over the entire time range of interest.
jdminmax = [jday(1951d,1,1d),jday(2049d,11,1d)]
if n_elements(info) EQ 0 then begin
    jplephread, '~/xtecal/clock/JPLEPH.405', info, raw, $
      [jdminmax(0)-30,jdminmax(1)+30]
endif

;; Initialize JDVERN, and define the initial time range to search for
;; in 1951.
jdvern = [0d]
jdminmax = jday(1951d,3,21d)+[-2d,2d]

;; Search for 99 equinoxes, starting in 1951
for i = 0, 98 do begin

iter = 0L

;; Iterate to find a solution.  The tolerance for completing the
;; solution is +/- 10 seconds (and no more than 10 iterations, but I
;; never did more than 3 iterations).
while (jdminmax(1)-jdminmax(0))*86400d GT 10d AND iter LT 10 do begin

    ;; Compute the Sun's position with respect to the earth using the
    ;; JPL ephemeris, referred to the TDT timescale.  The positions
    ;; are referred to the earth mean equator of epoch J2000.0.
    tt = rvspan(jdminmax, 100)
    jplephinterp, info, raw, tt, xs, ys, zs, object='SUN', center='EARTH'

    ;; Compute the earth precession and nutation.  This is important
    ;; because the equinox is referred to the true ecliptic *of date*,
    ;; which is defined with respect to the true earth equator of
    ;; date.  To convert from J2000, we need to precess and nutate.
    hprnutang, tt, zeta, theta, z, dpsi, deps, $
      true_obliq=eps, mean_obliq=eps0,$
      jd_ut1=ut1
    
    ;; Compose the precession and nutation rotations as a series of
    ;; quaternions
    qq = qtinv(qteuler(['z','y','z',      'x','z','x'], $
                       -zeta, +theta, -z, +eps0, -dpsi, -deps))

    ;; Apply the rotations to the Sun positions
    rs = transpose([[xs],[ys],[zs]])
    rrs = sqrt(total(rs^2,1))
    rse = qtvrot(rs, qq)

    ;; Compute the ecliptic longitude of the sun, ATAN(Y,X), and find
    ;; where it crosses zero.
    lse = atan(rse(1,*), rse(0,*))*180d/!dpi
    ams = (lse) MOD 360d
    wh = where(ams GT 180, ct) ;; Handle phase wraps
    if ct GT 0 then ams(wh) = ams(wh) - 360 
    wh = where(ams(1:*)*ams LT 0 AND abs(ams(1:*)-ams) LT 50, ct)

    if ct EQ 0 then message, 'ERROR: could not find root!'
    ;; Zoom in on the new bracketing range
    jdminmax = [tt(wh),tt(wh+1)]
    utminmax = [ut1(wh), ut1(wh+1)]
    iter = iter + 1
endwhile

;; Account for light travel time - we won't see the sun pass zero
;; longitude until 8 minutes later.  Ignore the aberration terms
jdvern1 = avg(utminmax)+avg(rrs)/(info.c/1000d)/86400d
jdvern = [jdvern, jdvern1]
jdminmax = jdvern1 + 365.25d + [-2d, +2d]

;; Print some status
daycnv, jdvern1, yr, mn, dy, hr
mi = floor((hr-floor(hr))*3600d)
se = mi MOD 60
mi = mi / 60
hr = floor(hr)
print, jdvern1, yr, mn, dy, hr, mi, se, $
  format='(D15.5," ",I4.4,"/",I2.2,"/",I2.2," ",I2.2,":",I2.2,":",I2.2)'
endfor

;; We are done!  Remainder of the code converts to calendar date and
;; saves it.
jdvern = jdvern(1:*)
jddiff = jdvern(1:*)-jdvern

daycnv, jdvern, yr, mn, dy, hr
mi = floor((hr-floor(hr))*3600d)
se = mi MOD 60
mi = mi / 60
hr = floor(hr)
jddiff1 = [0d, jddiff]

openw, 50, 'vernequinox.dat'
printf, 50, ';     Epochs of Vernal Equinox'
printf, 50, ';  Julian Date   Calendar Date UTC     Diff'
for i = 0, n_elements(jdvern)-1 do $
  printf, 50, jdvern(i), yr(i), mn(i), dy(i), hr(i), mi(i), se(i), jddiff1(i),$
  format='(D15.5," ",I4.4,"/",I2.2,"/",I2.2," ",I2.2,":",I2.2,":",I2.2," ",D10.5)'
close, 50

; plot, rebin(yr(1:*)+0d,14), rebin(jddiff-avg(jddiff),14)*86400d, charsize=1.75, charthick=2, ythick=2, xthick=2, thick=2, xtitle='Epoch Year', ytitle='Length of Equinoctial Year Minus Average', yrange=[-500,500]
; oplot, yr, (jddiff-avg(jddiff))*86400d



end