From cb9d9b4b16c7de668772a93899f65ee313efac91 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Henrik=20Grubbstr=C3=B6m=20=28Grubba=29?= <grubba@grubba.org> Date: Sun, 9 Sep 2012 21:06:56 +0200 Subject: [PATCH] Calendar.Event: Compensate for leap seconds. The solar event calculations now compensate for the difference between TDT (Terrestrial Dynamical Time) and UTC. Also fixes some AutoDoc mk II markup bugs. --- lib/modules/Calendar.pmod/Event.pmod | 99 +++++++++++++++++++++++++--- 1 file changed, 91 insertions(+), 8 deletions(-) diff --git a/lib/modules/Calendar.pmod/Event.pmod b/lib/modules/Calendar.pmod/Event.pmod index 5437abab20..cb0ff45659 100644 --- a/lib/modules/Calendar.pmod/Event.pmod +++ b/lib/modules/Calendar.pmod/Event.pmod @@ -876,14 +876,14 @@ class Solar(int|void event_type) inherit Day_Event; //! @array - //! @item + //! @elem array(float) 0.. //! @array - //! @item 0 amplitude - //! Days. - //! @item 1 phase - //! Radians @ year 2000. - //! @item 2 period - //! Radians/Year. + //! @elem float 0 + //! Amplitude in days. + //! @elem float 1 + //! Cosine phase in radians in year 2000. + //! @elem float 2 + //! Period in radians/year. //! @endarray //! @endarray protected constant periodic_table = ({ @@ -913,6 +913,88 @@ class Solar(int|void event_type) ({ 0.00008, 0.2697, 294.24635, }), }); + //! Polynomial terms for calculating the difference between TDT and UTC. + //! + //! The polynomials have been taken from NASA: + //! @url{http://eclipse.gsfc.nasa.gov/SEhelp/deltatpoly2004.html@} + //! + //! Each entry is an @expr{array(float)@}: + //! @array + //! @elem float 0 + //! End year for the range. + //! @elem float 1 + //! Year offset. + //! @elem float 2 + //! Year divisor. + //! @elem float 3.. + //! Polynomial factors with the highest exponent first. + //! @endarray + protected constant deltat_polynomials = ({ + ({ -500.0, -1820.0, 100.0, + -20.0, 0.0, 32.0, }), + ({ 500.0, 0.0, 100.0, + 0.0090316521, 0.022174192, -0.1798452, + -5.952053, 33.78311, -1014.41, 10583.6, }), + ({ 1600.0, -1000.0, 100.0, + 0.0083572073, -0.005050998, -0.8503463, + 0.319781, 71.23472, -556.01, 1574.2, }), + ({ 1700.0, -1600.0, 1.0, + 0.000140272128, -0.01532, -0.9808, 120.0, }), + ({ 1800.0, -1700.0, 1.0, + 0.000000851788756, 0.00013336, -0.0059285, 0.1603, 8.83, }), + ({ 1860.0, -1800.0, 1.0, + 0.000000000875, -0.0000001699, 0.0000121272, -0.00037436, + 0.0041116, 0.0068612, -0.332447, 13.72, }), + ({ 1900.0, -1860.0, 1.0, + 0.000004288643, -0.0004473624, 0.01680668, -0.251754, 0.5737, 7.62, }), + ({ 1920.0, -1900.0, 1.0, + -0.000197, 0.0061966, -0.0598939, 1.494119, -2.79, }), + ({ 1941.0, -1920.0, 1.0, + 0.0020936, -0.076100, 0.84493, 21.20, }), + ({ 1961.0, -1950.0, 1.0, + 0.0003926188, 0.0042918455, 0.407, 29.07, }), + ({ 1986.0, -1975.0, 1.0, + 0.001392758, 0.0038461538, 1.067, 45.45, }), + ({ 2005.0, -2000.0, 1.0, + 0.00002373599, 0.000651814, 0.0017275, -0.060374, 0.3345, 63.86, }), + ({ 2050.0, -2000.0, 1.0, + 0.005589, 0.32217, 62.92 }), + ({ 2150.0, -1820.0, 100.0, + 32.0-185.724, 56.28, -20.0, }), + ({ Math.inf, -1820.0, 100.0, + 32.0, 0.0, -20.0, }), + }); + + //! Terrestrial Dynamical Time difference from standard time. + //! + //! @returns + //! Difference between TDT and UTC in seconds at the + //! specified time. + int deltat(int unadjusted_utc) + { + // Approximation of the year. This ought to be good enough for + // most purposes given the uncertainty in the table values. + // 31556952 == 365.2425 * 24 * 60 * 60. + float y = 1970.0 + unadjusted_utc/31556952.0; + + array(float) polynomial; + int l, c, h = sizeof(deltat_polynomials); + do { + c = (l + h)/2; + polynomial = deltat_polynomials[c]; + if (y < polynomial[0]) h = c; + else l = c + 1; + } while (l < h); + float u = (y + polynomial[1])/polynomial[2]; + + float deltat = 0.0; + foreach(polynomial; int i; float factor) { + if (i < 3) continue; + deltat = deltat * u + factor; + } + return (int)round(deltat); + } + //! Calculate the next event. //! //! Based on Meeus Astronomical Algorithms Chapter 27. @@ -1035,7 +1117,8 @@ class Solar(int|void event_type) } // Convert into an UTC timestamp. - return (new_jd - 2440588)*86400 + (int)(offset * 86400.0); + int utc = (new_jd - 2440588)*86400 + (int)(offset * 86400.0); + return utc - deltat(utc); } Calendar.TimeRanges.TimeRange next(void|Calendar.TimeRanges.TimeRange from, -- GitLab