;; This code computes the times of full moon 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 full moon are computed in Julian Days (UTC) and placed ;; in the JDFULL variable, and saved to a file called "fullmoon.dat". ;; The differences between successive full moons 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. jdfull = [0d] jdminmax = jday(1951d,1,23d)+[-2d,2d] ;; Search for 1200 full moons, which takes us almost up to the year ;; 2050 (more like 2048) for i = 0, 1200 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 tt = rvspan(jdminmax, 100) ;; Compute the Sun/Moon 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. jplephinterp, info, raw, tt, xs, ys, zs, object='SUN', center='EARTH' jplephinterp, info, raw, tt, xm, ym, zm, object='MOON', 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'], $ -zeta, +theta, -z, +eps0, -dpsi)) ;; Apply the rotations to the Sun/Moon positions rs = transpose([[xs],[ys],[zs]]) rm = transpose([[xm],[ym],[zm]]) rse = qtvrot(rs, qq) rme = qtvrot(rm, qq) ;; Compute the ecliptic longitude of the sun/moon, ATAN(Y,X), and ;; find where LAMBDA_MOON - LAMBDA_SUN + 180 crosses zero lse = atan(rse(1,*), rse(0,*))*180d/!dpi lme = atan(rme(1,*), rme(0,*))*180d/!dpi ams = (lme-lse+180) 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 ;; Print some status jdfull1 = avg(utminmax) jdfull = [jdfull, jdfull1] jdminmax = jdfull1 + 29.6d + [-2d, +2d] daycnv, jdfull1, yr, mn, dy, hr mi = floor((hr-floor(hr))*3600d) se = mi MOD 60 mi = mi / 60 hr = floor(hr) print, jdfull1, 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. jdfull = jdfull(1:*) jddiff = jdfull(1:*)-jdfull daycnv, jdfull, 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, 'fullmoon.dat' printf, 50, '; Epochs of Full Moon' printf, 50, '; Julian Date Calendar Date UTC Diff' for i = 0, n_elements(jdfull)-1 do $ printf, 50, jdfull(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 end