Package: astronomical-almanac / 5.6-7

deltat.c.patch Patch series | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
Description: Update delta T tables.
Author: Stephen Moshier
Forwarded: patch is from upstream
--- astronomical-almanac-5.6/deltat.c	2003-05-04 11:47:38.000000000 -0400
+++ new/deltat.c	2020-03-29 11:40:20.124499191 -0400
@@ -1,32 +1,18 @@
 /* DeltaT = Ephemeris Time - Universal Time
  *
- * The tabulated values of deltaT, in hundredths of a second,
- * were taken from The Astronomical Almanac, page K8.  The program
- * adjusts for a value of secular tidal acceleration ndot. It is -25.8
- * arcsec per century squared for JPL's DE403 ephemeris.
- * ELP2000 and DE200 use the value -23.8946.
- *
- * The tabulated range is 1620.0 through 2003.0.  Bessel's interpolation
- * formula is implemented to obtain fourth order interpolated values at
- * intermediate times.
- *
- * Updated deltaT predictions can be obtained from this network archive:
+ * Tabulated values of deltaT, in hundredths of a second, are
+ * from The Astronomical Almanac and current IERS reports.
+ * Updated deltaT predictions can be obtained from this network archive,
  *    http://maia.usno.navy.mil
- * Currently (as of 2002) available series are
- *    tai-utc.dat  Changes by 1 whenever there is a leap second
- *    finals.all   EOP including UT1-UTC, always less than 1 second
- * from which deltaT = 32.184 + (tai-utc) - (UT1-UTC)
- *
- * For dates earlier than the tabulated range, the program
- * calculates approximate formulae of Stephenson and Morrison
- * or K. M. Borkowski.  These approximations have an estimated
- * error of 15 minutes at 1500 B.C.  They are not adjusted for small
- * improvements in the current estimate of ndot because the formulas
- * were derived from studies of ancient eclipses and other historical
- * information, whose interpretation depends only partly on ndot.
+ * and appended to the table.
+ *
+ * A table of values for the pre-telescopic period was taken from
+ * Morrison and Stephenson (2004).  The overall tabulated range is
+ * -1000.0 through 2019.0.  Values at intermediate times are interpolated
+ * from the tables.
  *
- * A quadratic extrapolation formula, that agrees in value and slope with
- * current data, predicts future values of deltaT.
+ * For dates earlier and later than the tabulated range, the program
+ * calculates a polynomial extrapolation formula.
  *
  * Input Y is the Julian epoch expressed in Julian years.  Y can be
  * found from the Julian date JD by
@@ -38,45 +24,60 @@
  *
  * References:
  *
+ * Morrison, L. V., and F. R. Stephenson, Historical values of the Earth's
+ * clock error deltat T and the calculation of eclipses. Journal for the
+ * History of Astronomy 35, 327-336 (2004)
+ *
  * Stephenson, F. R., and L. V. Morrison, "Long-term changes
  * in the rotation of the Earth: 700 B.C. to A.D. 1980,"
  * Philosophical Transactions of the Royal Society of London
  * Series A 313, 47-70 (1984)
  *
- * Borkowski, K. M., "ELP2000-85 and the Dynamical Time
- * - Universal Time relation," Astronomy and Astrophysics
- * 205, L8-L10 (1988)
- * Borkowski's formula is derived from eclipses going back to 2137 BC
- * and uses lunar position based on tidal coefficient of -23.9 arcsec/cy^2.
- *
  * Chapront-Touze, Michelle, and Jean Chapront, _Lunar Tables
  * and Programs from 4000 B.C. to A.D. 8000_, Willmann-Bell 1991
- * Their table agrees with the one here, but the entries are
- * rounded to the nearest whole second.
  *
  * Stephenson, F. R., and M. A. Houlden, _Atlas of Historical
  * Eclipse Maps_, Cambridge U. Press (1986)
+ *
  */
 
 #include "kep.h"
 
 /* If the following number (read from the file aa.ini)
- * is nonzero, then the program will return it
+ * is nonzero, then the program will return it for delta T
  * and not calculate anything.
  */
 double dtgiven = 0.0;
 extern double dtgiven;
-
+double deltat_value;
 
 #define DEMO 0
 #define TABSTART 1620
-#define TABEND 2013
+#define TABEND 2019
 #define TABSIZ (TABEND - TABSTART + 1)
 
-/* Note, Stephenson and Morrison's table starts at the year 1630.
- * The Chapronts' table does not agree with the Almanac prior to 1630.
- * The actual accuracy decreases rapidly prior to 1780.
- */
+/* Morrison and Stephenson (2004)
+   This table covers -1000 through 1700 in 100-year steps.
+   Values are in whole seconds.
+   Estimated standard error at -1000 is 640 seconds; at 1600, 20 seconds.
+   The first value in the table has been adjusted 28 sec for
+   continuity with their long-term quadratic extrapolation formula.
+   The last value in this table agrees with the AA table at 1700,
+   so there is no discontinuity at either endpoint.  */
+#define MS_SIZ 28
+short m_s[MS_SIZ] = {
+/* -1000 to -100 */
+ 25428, 23700, 22000, 21000, 19040, 17190, 15530, 14080, 12790, 11640,
+/* 0 to 900 */
+ 10580, 9600, 8640, 7680, 6700, 5710, 4740, 3810, 2960, 2200,
+/* 1000 to 1700 */
+ 1570, 1090, 740, 490, 320, 200, 120, 9,
+};
+
+
+/*    Entries prior to 1955 in the following table are from
+      the 1984 Astronomical Almanac and assume ndot = -26.0.
+      For dates prior to 1700, the above table is used instead of this one. */
 short dt[TABSIZ] = {
 /* 1620.0 thru 1659.0 */
 12400, 11900, 11500, 11000, 10600, 10200, 9800, 9500, 9100, 8800,
@@ -124,13 +125,11 @@
  2915, 2957, 2997, 3036, 3072, 3107, 3135, 3168, 3218, 3268,
  3315, 3359, 3400, 3447, 3503, 3573, 3654, 3743, 3829, 3920,
  4018, 4117, 4223, 4337, 4449, 4548, 4646, 4752, 4853, 4959,
-/* 1980.0 thru 2003.0 */
+/* 1980.0 thru 2019.0 */
  5054, 5138, 5217, 5296, 5379, 5434, 5487, 5532, 5582, 5630,
  5686, 5757, 5831, 5912, 5998, 6078, 6163, 6230, 6297, 6347,
- 6383, 6409, 6430, 6447,
-/* Extrapolated values */
-                         6456, 6600, 6700, 6800, 6900, 7000,
- 7100, 7200, 7300, 7400,
+ 6383, 6409, 6430, 6447, 6457, 6469, 6485, 6515, 6546, 6578,
+ 6607, 6632, 6660, 6691, 6728, 6764, 6810, 6859, 6897, 6922,
 };
 
 
@@ -143,54 +142,86 @@
 int i, iy, k;
 
 if( dtgiven != 0.0 )
-	return( dtgiven );
+{
+    deltat_value = dtgiven;
+    return dtgiven;
+}
 
 if( Y > TABEND )
 	{
-#if 0
-/* Morrison, L. V. and F. R. Stephenson, "Sun and Planetary System"
- * vol 96,73 eds. W. Fricke, G. Teleki, Reidel, Dordrecht (1982)
- */
-	B = 0.01*(Y-1800.0) - 0.1;
-	ans = -15.0 + 32.5*B*B;
-	return(ans);
-#else
-/* Extrapolate forward by a second-degree curve that agrees with
-   the most recent data in value and slope, and vaguely fits
-   over the past century.  This idea communicated by Paul Muller,
-   who says NASA used to do something like it.  */
-	B = Y - TABEND;
-	/* slope */
-	p = dt[TABSIZ-1] - dt[TABSIZ-2];
-	/* square term */
-	ans = (dt[TABSIZ - 101] - (dt[TABSIZ - 1] - 100.0 * p)) * 1e-4;
-        ans = 0.01 * (dt[TABSIZ-1] + p * B + ans * B * B);
-	if( prtflg )
-	  printf("[extrapolated deltaT] ");
-	return(ans);
-#endif
-	}
+    /* Extrapolate future values beyond the lookup table.  */
+    if (Y > (TABEND + 100.0))
+    {
+        /* Morrison & Stephenson (2004) long-term curve fit.  */
+        B = 0.01 * (Y - 1820.0);
+        ans = 32.0 * B * B - 20.0;
+    }
+    else
+    {
+        double a, b, c, d, m0, m1;
+
+        /* Cubic interpolation between last tabulated value
+           and long-term curve evaluated at 100 years later.  */
+
+        /* Last tabulated delta T value. */
+        a = 0.01 * dt[TABSIZ-1];
+        /* Approximate slope in past 10 years. */
+        b = 0.001 * (dt[TABSIZ-1] - dt[TABSIZ - 11]);
+
+        /* Long-term curve 100 years hence. */
+        B = 0.01 * (TABEND + 100.0 - 1820.0);
+        m0 = 32.0 * B*B - 20.0;
+        /* Its slope. */
+        m1 = 0.64 * B;
+
+        /* Solve for remaining coefficients of an interpolation polynomial
+           that agrees in value and slope at both ends of the 100-year
+           interval. */
+        d = 2.0e-6 * (50.0 * (m1 + b) - m0 + a);
+        c = 1.0e-4 * (m0 - a - 100.0 * b - 1.0e6 * d);
+
+        /* Note, the polynomial coefficients do not depend on Y.
+           A given tabulation and long-term formula
+           determine the polynomial.
+           Thus, for the IERS table ending at 2011.0, the coefficients are
+           a = 66.32
+           b = 0.223
+           c = 0.03231376
+           d = -0.0001607784
+        */
+
+        /* Compute polynomial value at desired time. */
+        p = Y - TABEND;
+        ans = a + p * (b  + p * (c + p * d));
+    }
+    return ans;
+        }
 
-if( Y < TABSTART )
-	{
-	if( Y >= 948.0 )
-		{
-/* Stephenson and Morrison, stated domain is 948 to 1600:
- * 25.5(centuries from 1800)^2 - 1.9159(centuries from 1955)^2
- */
-		B = 0.01*(Y - 2000.0);
-		ans = (23.58 * B + 100.3)*B + 101.6;
-		}
-	else
-		{
-/* Borkowski */
-		B = 0.01*(Y - 2000.0)  +  3.75;
-		ans = 35.0 * B * B  +  40.;
-		}
-	return(ans);
-	}
 
-/* Besselian interpolation from tabulated values.
+/* Use Morrison and Stephenson (2004) prior to the year 1700.  */
+if( Y < 1700.0 )
+{
+    if (Y <= -1000.0)
+    {
+        /* Morrison and Stephenson long-term fit.  */
+        B = 0.01 * (Y - 1820.0);
+        ans = 32.0 * B * B - 20.0;
+    }
+    else
+    {
+        /* Morrison and Stephenson recommend linear interpolation
+           between tabulations. */
+        iy = Y;
+        iy = (iy + 1000) / 100;  /* Integer index into the table. */
+        B = -1000 + 100 * iy;    /* Starting year of tabulated interval.  */
+        p = m_s[iy];
+        ans = p + 0.01 * (Y - B) * (m_s[iy + 1] - p);
+    }
+    return(ans);
+}
+
+/* Besselian interpolation between tabulated values
+ * in the telescopic era.
  * See AA page K11.
  */
 
@@ -264,6 +295,12 @@
 #endif
 
 done:
+
+ans *= 0.01;
+
+
+#if 0 /* ndot = -26.0 assumed; no correction. */
+
 /* Astronomical Almanac table is corrected by adding the expression
  *     -0.000091 (ndot + 26)(year-1955)^2  seconds
  * to entries prior to 1955 (AA page K8), where ndot is the secular
@@ -272,7 +309,6 @@
  * Entries after 1955 are referred to atomic time standards and
  * are not affected by errors in Lunar or planetary theory.
  */
-ans *= 0.01;
 if( Y < 1955.0 )
 	{
 	B = (Y - 1955.0);
@@ -282,6 +318,9 @@
 	ans += -0.000091 * (-23.8946 + 26.0) * B * B;
 #endif
 	}
+
+#endif /* 0 */
+
 return( ans );
 }
 
@@ -309,12 +348,14 @@
 		break;
 	case 1:
 		TDT = JD;
-		UT = TDT - deltat(T)/86400.0;
+                deltat_value = deltat(T);
+		UT = TDT - deltat_value/86400.0;

 		jtocal(UT); /* display the other date */
 		break;
 	case 2:
 		UT = JD;
-		TDT = UT + deltat(T)/86400.0;
+                deltat_value = deltat(T);
+		TDT = UT + deltat_value/86400.0;

 		jtocal(TDT);
 		break;
 	}