1 /* DeltaT = Ephemeris Time - Universal Time
2 *
3 * Adapted 2011/4/14 from Stephen Moshier <moshier@world.std.com>,
4 * cosmetic changes only.
5 *
6 * Compile as follows to create stand-alone test program:
7 * cc -DTEST_MAIN deltat.c libastro.a
8 *
9 * Tabulated values of deltaT, in hundredths of a second, are
10 * from The Astronomical Almanac and current IERS reports.
11 * A table of values for the pre-telescopic period was taken from
12 * Morrison and Stephenson (2004). The overall tabulated range is
13 * -1000.0 through 2018.0. Values at intermediate times are interpolated
14 * from the tables.
15 *
16 * For dates earlier and later than the tabulated range, the program
17 * calculates a polynomial extrapolation formula.
18 *
19 * Updated deltaT predictions can be obtained from this network archive,
20 * http://maia.usno.navy.mil
21 * then appended to the dt[] table and update TABEND.
22 *
23 * Input is XEphem's MJD, output is ET-UT in seconds.
24 *
25 *
26 * References:
27 *
28 * Morrison, L. V., and F. R. Stephenson, Historical values of the Earth's
29 * clock error deltat T and the calculation of eclipses. Journal for the
30 * History of Astronomy 35, 327-336 (2004)
31 *
32 * Stephenson, F. R., and L. V. Morrison, "Long-term changes
33 * in the rotation of the Earth: 700 B.C. to A.D. 1980,"
34 * Philosophical Transactions of the Royal Society of London
35 * Series A 313, 47-70 (1984)
36 *
37 * Chapront-Touze, Michelle, and Jean Chapront, _Lunar Tables
38 * and Programs from 4000 B.C. to A.D. 8000_, Willmann-Bell 1991
39 *
40 * Stephenson, F. R., and M. A. Houlden, _Atlas of Historical
41 * Eclipse Maps_, Cambridge U. Press (1986)
42 *
43 */
44
45 #include <math.h>
46
47 #include "astro.h"
48
49 #define TABSTART 1620
50 #define TABEND 2018
51 #define TABSIZ (TABEND - TABSTART + 1)
52
53 /* Morrison and Stephenson (2004)
54 * This table covers -1000 through 1700 in 100-year steps.
55 * Values are in whole seconds.
56 * Estimated standard error at -1000 is 640 seconds; at 1600, 20 seconds.
57 * The first value in the table has been adjusted 28 sec for
58 * continuity with their long-term quadratic extrapolation formula.
59 * The last value in this table agrees with the AA table at 1700,
60 * so there is no discontinuity at either endpoint.
61 */
62 #define MS_SIZ 28
63 short m_s[MS_SIZ] = {
64 /* -1000 to -100 */
65 25428, 23700, 22000, 21000, 19040, 17190, 15530, 14080, 12790, 11640,
66
67 /* 0 to 900 */
68 10580, 9600, 8640, 7680, 6700, 5710, 4740, 3810, 2960, 2200,
69
70 /* 1000 to 1700 */
71 1570, 1090, 740, 490, 320, 200, 120, 9,
72 };
73
74
75 /* Entries prior to 1955 in the following table are from
76 * the 2020 Astronomical Almanac and assume ndot = -26.0.
77 * For dates prior to 1700, the above table is used instead of this one.
78 */
79 short dt[TABSIZ] = {
80 /* 1620.0 thru 1659.0 */
81 12400, 11900, 11500, 11000, 10600, 10200, 9800, 9500, 9100, 8800,
82 8500, 8200, 7900, 7700, 7400, 7200, 7000, 6700, 6500, 6300,
83 6200, 6000, 5800, 5700, 5500, 5400, 5300, 5100, 5000, 4900,
84 4800, 4700, 4600, 4500, 4400, 4300, 4200, 4100, 4000, 3800,
85
86 /* 1660.0 thru 1699.0 */
87 3700, 3600, 3500, 3400, 3300, 3200, 3100, 3000, 2800, 2700,
88 2600, 2500, 2400, 2300, 2200, 2100, 2000, 1900, 1800, 1700,
89 1600, 1500, 1400, 1400, 1300, 1200, 1200, 1100, 1100, 1000,
90 1000, 1000, 900, 900, 900, 900, 900, 900, 900, 900,
91
92 /* 1700.0 thru 1739.0 */
93 900, 900, 900, 900, 900, 900, 900, 900, 1000, 1000,
94 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1100, 1100, 1100,
95 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100,
96 1100, 1100, 1100, 1100, 1200, 1200, 1200, 1200, 1200, 1200,
97
98 /* 1740.0 thru 1779.0 */
99 1200, 1200, 1200, 1200, 1300, 1300, 1300, 1300, 1300, 1300,
100 1300, 1400, 1400, 1400, 1400, 1400, 1400, 1400, 1500, 1500,
101 1500, 1500, 1500, 1500, 1500, 1600, 1600, 1600, 1600, 1600,
102 1600, 1600, 1600, 1600, 1600, 1700, 1700, 1700, 1700, 1700,
103
104 /* 1780.0 thru 1799.0 */
105 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700,
106 1700, 1700, 1600, 1600, 1600, 1600, 1500, 1500, 1400, 1400,
107
108 /* 1800.0 thru 1819.0 */
109 1370, 1340, 1310, 1290, 1270, 1260, 1250, 1250, 1250, 1250,
110 1250, 1250, 1250, 1250, 1250, 1250, 1250, 1240, 1230, 1220,
111
112 /* 1820.0 thru 1859.0 */
113 1200, 1170, 1140, 1110, 1060, 1020, 960, 910, 860, 800,
114 750, 700, 660, 630, 600, 580, 570, 560, 560, 560,
115 570, 580, 590, 610, 620, 630, 650, 660, 680, 690,
116 710, 720, 730, 740, 750, 760, 770, 770, 780, 780,
117
118 /* 1860.0 thru 1899.0 */
119 788, 782, 754, 697, 640, 602, 541, 410, 292, 182,
120 161, 10, -102, -128, -269, -324, -364, -454, -471, -511,
121 -540, -542, -520, -546, -546, -579, -563, -564, -580, -566,
122 -587, -601, -619, -664, -644, -647, -609, -576, -466, -374,
123
124 /* 1900.0 thru 1939.0 */
125 -272, -154, -2, 124, 264, 386, 537, 614, 775, 913,
126 1046, 1153, 1336, 1465, 1601, 1720, 1824, 1906, 2025, 2095,
127 2116, 2225, 2241, 2303, 2349, 2362, 2386, 2449, 2434, 2408,
128 2402, 2400, 2387, 2395, 2386, 2393, 2373, 2392, 2396, 2402,
129
130 /* 1940.0 thru 1979.0 */
131 2433, 2483, 2530, 2570, 2624, 2677, 2728, 2778, 2825, 2871,
132 2915, 2957, 2997, 3036, 3072, 3107, 3135, 3168, 3218, 3268,
133 3315, 3359, 3400, 3447, 3503, 3573, 3654, 3743, 3829, 3920,
134 4018, 4117, 4223, 4337, 4449, 4548, 4646, 4752, 4853, 4959,
135
136 /* 1980.0 thru 2018.0 */
137 5054, 5138, 5217, 5296, 5379, 5434, 5487, 5532, 5582, 5630,
138 5686, 5757, 5831, 5912, 5998, 6078, 6163, 6230, 6297, 6347,
139 6383, 6409, 6430, 6447, 6457, 6469, 6485, 6515, 6546, 6578,
140 6607, 6632, 6660, 6691, 6728, 6764, 6810, 6859, 6897,
141 };
142
143
144 /* Given MJD return DeltaT = ET - UT1 in seconds. Describes the irregularities
145 * of the Earth rotation rate in the ET time scale.
146 */
147 double
deltat(double mj)148 deltat(double mj)
149 {
150 static double ans, lastmj;
151 double Y, p, B;
152 int d[6];
153 int i, iy, k;
154
155 if (mj == lastmj)
156 return (ans);
157 lastmj = mj;
158
159 mjd_year (mj, &Y);
160
161 if( Y > TABEND ) {
162 /* Extrapolate future values beyond the lookup table. */
163 if (Y > (TABEND + 100.0)) {
164 /* Morrison & Stephenson (2004) long-term curve fit. */
165 B = 0.01 * (Y - 1820.0);
166 ans = 32.0 * B * B - 20.0;
167
168 } else {
169
170 double a, b, c, d, m0, m1;
171
172 /* Cubic interpolation between last tabulated value
173 * and long-term curve evaluated at 100 years later.
174 */
175
176 /* Last tabulated delta T value. */
177 a = 0.01 * dt[TABSIZ-1];
178 /* Approximate slope in past 10 years. */
179 b = 0.001 * (dt[TABSIZ-1] - dt[TABSIZ - 11]);
180
181 /* Long-term curve 100 years hence. */
182 B = 0.01 * (TABEND + 100.0 - 1820.0);
183 m0 = 32.0 * B*B - 20.0;
184 /* Its slope. */
185 m1 = 0.64 * B;
186
187 /* Solve for remaining coefficients of an interpolation polynomial
188 * that agrees in value and slope at both ends of the 100-year
189 * interval.
190 */
191 d = 2.0e-6 * (50.0 * (m1 + b) - m0 + a);
192 c = 1.0e-4 * (m0 - a - 100.0 * b - 1.0e6 * d);
193
194 /* Note, the polynomial coefficients do not depend on Y.
195 * A given tabulation and long-term formula
196 * determine the polynomial.
197 * Thus, for the IERS table ending at 2011.0, the coefficients are
198 * a = 66.32
199 * b = 0.223
200 * c = 0.03231376
201 * d = -0.0001607784
202 */
203
204 /* Compute polynomial value at desired time. */
205 p = Y - TABEND;
206 ans = a + p * (b + p * (c + p * d));
207 }
208
209 return (ans);
210 }
211
212
213 /* Use Morrison and Stephenson (2004) prior to the year 1700. */
214 if( Y < 1700.0 ) {
215 if (Y <= -1000.0) {
216 /* Morrison and Stephenson long-term fit. */
217 B = 0.01 * (Y - 1820.0);
218 ans = 32.0 * B * B - 20.0;
219
220 } else {
221
222 /* Morrison and Stephenson recommend linear interpolation
223 * between tabulations.
224 */
225 iy = Y;
226 iy = (iy + 1000) / 100; /* Integer index into the table. */
227 B = -1000 + 100 * iy; /* Starting year of tabulated interval. */
228 p = m_s[iy];
229 ans = p + 0.01 * (Y - B) * (m_s[iy + 1] - p);
230 }
231
232 return (ans);
233 }
234
235 /* Besselian interpolation between tabulated values
236 * in the telescopic era.
237 * See AA page K11.
238 */
239
240 /* Index into the table. */
241 p = floor(Y);
242 iy = (int) (p - TABSTART);
243 /* Zeroth order estimate is value at start of year */
244 ans = dt[iy];
245 k = iy + 1;
246 if( k >= TABSIZ )
247 goto done; /* No data, can't go on. */
248
249 /* The fraction of tabulation interval */
250 p = Y - p;
251
252 /* First order interpolated value */
253 ans += p*(dt[k] - dt[iy]);
254 if( (iy-1 < 0) || (iy+2 >= TABSIZ) )
255 goto done; /* can't do second differences */
256
257 /* Make table of first differences */
258 k = iy - 2;
259 for (i=0; i<5; i++) {
260 if( (k < 0) || (k+1 >= TABSIZ) )
261 d[i] = 0;
262 else
263 d[i] = dt[k+1] - dt[k];
264 k += 1;
265 }
266
267 /* Compute second differences */
268 for( i=0; i<4; i++ )
269 d[i] = d[i+1] - d[i];
270 B = 0.25*p*(p-1.0);
271 ans += B*(d[1] + d[2]);
272 if (iy+2 >= TABSIZ)
273 goto done;
274
275 /* Compute third differences */
276 for( i=0; i<3; i++ )
277 d[i] = d[i+1] - d[i];
278 B = 2.0*B/3.0;
279 ans += (p-0.5)*B*d[1];
280 if ((iy-2 < 0) || (iy+3 > TABSIZ) )
281 goto done;
282
283 /* Compute fourth differences */
284 for( i=0; i<2; i++ )
285 d[i] = d[i+1] - d[i];
286 B = 0.125*B*(p+1.0)*(p-2.0);
287 ans += B*(d[0] + d[1]);
288
289 done:
290
291 ans *= 0.01;
292
293 #if 0 /* ndot = -26.0 assumed; no correction. */
294
295 /* Astronomical Almanac table is corrected by adding the expression
296 * -0.000091 (ndot + 26)(year-1955)^2 seconds
297 * to entries prior to 1955 (AA page K8), where ndot is the secular
298 * tidal term in the mean motion of the Moon.
299 *
300 * Entries after 1955 are referred to atomic time standards and
301 * are not affected by errors in Lunar or planetary theory.
302 */
303 if( Y < 1955.0 )
304 {
305 B = (Y - 1955.0);
306 #if 1
307 ans += -0.000091 * (-25.8 + 26.0) * B * B;
308 #else
309 ans += -0.000091 * (-23.8946 + 26.0) * B * B;
310 #endif
311 }
312
313 #endif /* 0 */
314
315 return( ans );
316 }
317
318
319 #ifdef TEST_MAIN
320
321 /* Exercise program.
322 */
323 #include <stdio.h>
324 #include <stdlib.h>
325
main(int ac,char * av[])326 int main(int ac, char *av[])
327 {
328 double ans, mj, y = atof(av[1]);
329 year_mjd (y, &mj);
330 ans = deltat(mj);
331 printf( "%.4lf\n", ans );
332 return (0);
333 }
334 #endif
335