| 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 2011.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 2011
|
|---|
| 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 1984 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 2011.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,
|
|---|
| 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
|
|---|
| 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 |
|
|---|
| 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
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
