| [1457] | 1 | /* DeltaT = Ephemeris Time - Universal Time
|
|---|
| 2 | *
|
|---|
| 3 | * original code by Stephen Moshier <moshier@world.std.com>,
|
|---|
| 4 | * adapted to xephem by Michael Sternberg <sternberg@physik.tu-chemnitz.de>
|
|---|
| 5 | * smoothed transitions and updated by Neal McBurnett <nealmcb@bell-labs.com>
|
|---|
| 6 | *
|
|---|
| 7 | **********************************************************************
|
|---|
| 8 | *
|
|---|
| 9 | * The tabulated values of deltaT, in hundredths of a second,
|
|---|
| 10 | * were taken from The Astronomical Almanac, page K8. The program
|
|---|
| 11 | * adjusts for a value of secular tidal acceleration ndot = -25.8
|
|---|
| 12 | * arcsec per century squared, the value used in JPL's DE403 ephemeris.
|
|---|
| 13 | * ELP2000 (and DE200) used the value -23.8946.
|
|---|
| 14 | *
|
|---|
| 15 | * The tabulated range is 1620.0 through 1998.0. Bessel's interpolation
|
|---|
| 16 | * formula is implemented to obtain fourth order interpolated values at
|
|---|
| 17 | * intermediate times.
|
|---|
| 18 | *
|
|---|
| 19 | * For dates earlier than the tabulated range, the program
|
|---|
| 20 | * calculates approximate formulae of Stephenson and Morrison
|
|---|
| 21 | * or K. M. Borkowski. These approximations have an estimated
|
|---|
| 22 | * error of 15 minutes at 1500 B.C. They are not adjusted for small
|
|---|
| 23 | * improvements in the current estimate of ndot because the formulas
|
|---|
| 24 | * were derived from studies of ancient eclipses and other historical
|
|---|
| 25 | * information, whose interpretation depends only partly on ndot.
|
|---|
| 26 | *
|
|---|
| 27 | * For future vaues of deltaT, the function smoothly transitions with
|
|---|
| 28 | * a linear segment back to Stephenson & Morrison's quadratic formula
|
|---|
| 29 | * in 2130.
|
|---|
| 30 | *
|
|---|
| [2551] | 31 | * Input is mj (modified julian date from MJD0 on). [stern]
|
|---|
| 32 | * Note that xephem uses a different epoch for this "mj" than the
|
|---|
| [1457] | 33 | * normal value of JD=240000.5.
|
|---|
| 34 | * See AA page B4.
|
|---|
| 35 | *
|
|---|
| [2551] | 36 | * Output double deltat(mj) is ET-UT1 in seconds.
|
|---|
| [1457] | 37 | *
|
|---|
| 38 | *
|
|---|
| 39 | * References:
|
|---|
| 40 | *
|
|---|
| 41 | * Stephenson, F. R., and L. V. Morrison, "Long-term changes
|
|---|
| 42 | * in the rotation of the Earth: 700 B.C. to A.D. 1980,"
|
|---|
| 43 | * Philosophical Transactions of the Royal Society of London
|
|---|
| 44 | * Series A 313, 47-70 (1984)
|
|---|
| 45 | *
|
|---|
| 46 | * Borkowski, K. M., "ELP2000-85 and the Dynamical Time
|
|---|
| 47 | * - Universal Time relation," Astronomy and Astrophysics
|
|---|
| 48 | * 205, L8-L10 (1988)
|
|---|
| 49 | * Borkowski's formula is derived from eclipses going back to 2137 BC
|
|---|
| 50 | * and uses lunar position based on tidal coefficient of -23.9 arcsec/cy^2.
|
|---|
| 51 | *
|
|---|
| 52 | * Chapront-Touze, Michelle, and Jean Chapront, _Lunar Tables
|
|---|
| 53 | * and Programs from 4000 B.C. to A.D. 8000_, Willmann-Bell 1991
|
|---|
| 54 | * Their table agrees with the one here, but the entries are
|
|---|
| 55 | * rounded to the nearest whole second.
|
|---|
| 56 | *
|
|---|
| 57 | * Stephenson, F. R., and M. A. Houlden, _Atlas of Historical
|
|---|
| 58 | * Eclipse Maps_, Cambridge U. Press (1986)
|
|---|
| 59 | *
|
|---|
| 60 | * from obsolete extrapolation code [stern]:
|
|---|
| 61 | * Morrison, L. V. and F. R. Stephenson, "Sun and Planetary System"
|
|---|
| 62 | * vol 96,73 eds. W. Fricke, G. Teleki, Reidel, Dordrecht (1982)
|
|---|
| 63 | *
|
|---|
| 64 | **********************************************************************
|
|---|
| 65 | *
|
|---|
| 66 | * changes by stern:
|
|---|
| 67 | * - adopted #include's for xephem
|
|---|
| 68 | * - made dt[] static
|
|---|
| [2551] | 69 | * - made mj the time argument [was: year Y].
|
|---|
| [1457] | 70 | * - updated observed and extrapolated data from tables at
|
|---|
| 71 | * ftp://maia.usno.navy.mil/ser7/ -- data deviated by up to 0.8 s
|
|---|
| 72 | * - removed references to "extern double dtgiven"
|
|---|
| 73 | * - removed DEMO #define and its references
|
|---|
| 74 | * - replaced treatment after TABEND by linear extrapolation instead
|
|---|
| 75 | * of second order version
|
|---|
| [2551] | 76 | * - installed lastmj cache (made ans static)
|
|---|
| [1457] | 77 | *
|
|---|
| 78 | * - no changes to table interpolation scheme and past extrapolations */
|
|---|
| 79 |
|
|---|
| [2551] | 80 | #include <math.h>
|
|---|
| 81 |
|
|---|
| [1457] | 82 | #include "astro.h"
|
|---|
| 83 |
|
|---|
| 84 | #define TABSTART 1620.0
|
|---|
| 85 | #define TABEND 2006.0
|
|---|
| 86 | #define TABSIZ 387
|
|---|
| 87 |
|
|---|
| 88 | /* Note, Stephenson and Morrison's table starts at the year 1630.
|
|---|
| 89 | * The Chapronts' table does not agree with the Almanac prior to 1630.
|
|---|
| 90 | * The actual accuracy decreases rapidly prior to 1780.
|
|---|
| 91 | */
|
|---|
| 92 | static short dt[TABSIZ] = {
|
|---|
| 93 | /* 1620.0 thru 1659.0 */
|
|---|
| 94 | 12400, 11900, 11500, 11000, 10600, 10200, 9800, 9500, 9100, 8800,
|
|---|
| 95 | 8500, 8200, 7900, 7700, 7400, 7200, 7000, 6700, 6500, 6300,
|
|---|
| 96 | 6200, 6000, 5800, 5700, 5500, 5400, 5300, 5100, 5000, 4900,
|
|---|
| 97 | 4800, 4700, 4600, 4500, 4400, 4300, 4200, 4100, 4000, 3800,
|
|---|
| 98 | /* 1660.0 thru 1699.0 */
|
|---|
| 99 | 3700, 3600, 3500, 3400, 3300, 3200, 3100, 3000, 2800, 2700,
|
|---|
| 100 | 2600, 2500, 2400, 2300, 2200, 2100, 2000, 1900, 1800, 1700,
|
|---|
| 101 | 1600, 1500, 1400, 1400, 1300, 1200, 1200, 1100, 1100, 1000,
|
|---|
| 102 | 1000, 1000, 900, 900, 900, 900, 900, 900, 900, 900,
|
|---|
| 103 | /* 1700.0 thru 1739.0 */
|
|---|
| 104 | 900, 900, 900, 900, 900, 900, 900, 900, 1000, 1000,
|
|---|
| 105 | 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1100, 1100, 1100,
|
|---|
| 106 | 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100,
|
|---|
| 107 | 1100, 1100, 1100, 1100, 1200, 1200, 1200, 1200, 1200, 1200,
|
|---|
| 108 | /* 1740.0 thru 1779.0 */
|
|---|
| 109 | 1200, 1200, 1200, 1200, 1300, 1300, 1300, 1300, 1300, 1300,
|
|---|
| 110 | 1300, 1400, 1400, 1400, 1400, 1400, 1400, 1400, 1500, 1500,
|
|---|
| 111 | 1500, 1500, 1500, 1500, 1500, 1600, 1600, 1600, 1600, 1600,
|
|---|
| 112 | 1600, 1600, 1600, 1600, 1600, 1700, 1700, 1700, 1700, 1700,
|
|---|
| 113 | /* 1780.0 thru 1799.0 */
|
|---|
| 114 | 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700,
|
|---|
| 115 | 1700, 1700, 1600, 1600, 1600, 1600, 1500, 1500, 1400, 1400,
|
|---|
| 116 | /* 1800.0 thru 1819.0 */
|
|---|
| 117 | 1370, 1340, 1310, 1290, 1270, 1260, 1250, 1250, 1250, 1250,
|
|---|
| 118 | 1250, 1250, 1250, 1250, 1250, 1250, 1250, 1240, 1230, 1220,
|
|---|
| 119 | /* 1820.0 thru 1859.0 */
|
|---|
| 120 | 1200, 1170, 1140, 1110, 1060, 1020, 960, 910, 860, 800,
|
|---|
| 121 | 750, 700, 660, 630, 600, 580, 570, 560, 560, 560,
|
|---|
| 122 | 570, 580, 590, 610, 620, 630, 650, 660, 680, 690,
|
|---|
| 123 | 710, 720, 730, 740, 750, 760, 770, 770, 780, 780,
|
|---|
| 124 | /* 1860.0 thru 1899.0 */
|
|---|
| 125 | 788, 782, 754, 697, 640, 602, 541, 410, 292, 182,
|
|---|
| 126 | 161, 10, -102, -128, -269, -324, -364, -454, -471, -511,
|
|---|
| 127 | -540, -542, -520, -546, -546, -579, -563, -564, -580, -566,
|
|---|
| 128 | -587, -601, -619, -664, -644, -647, -609, -576, -466, -374,
|
|---|
| 129 | /* 1900.0 thru 1939.0 */
|
|---|
| 130 | -272, -154, -2, 124, 264, 386, 537, 614, 775, 913,
|
|---|
| 131 | 1046, 1153, 1336, 1465, 1601, 1720, 1824, 1906, 2025, 2095,
|
|---|
| 132 | 2116, 2225, 2241, 2303, 2349, 2362, 2386, 2449, 2434, 2408,
|
|---|
| 133 | 2402, 2400, 2387, 2395, 2386, 2393, 2373, 2392, 2396, 2402,
|
|---|
| 134 | /* 1940.0 thru 1979.0 */
|
|---|
| 135 | 2433, 2483, 2530, 2570, 2624, 2677, 2728, 2778, 2825, 2871,
|
|---|
| 136 | 2915, 2957, 2997, 3036, 3072, 3107, 3135, 3168, 3218, 3268,
|
|---|
| 137 | 3315, 3359, 3400, 3447, 3503, 3573, 3654, 3743, 3829, 3920,
|
|---|
| 138 | 4018, 4117, 4223, 4337, 4449, 4548, 4646, 4752, 4853, 4959,
|
|---|
| 139 | /* 1980.0 thru 1995.0 */
|
|---|
| 140 | 5054, 5138, 5217, 5296, 5379, 5434, 5487, 5532, 5582, 5630,
|
|---|
| 141 | 5686, 5757, 5831, 5912, 5998, 6078,
|
|---|
| 142 | /* new USNO data (stern) */
|
|---|
| 143 | 6163, 6230,
|
|---|
| 144 | /* 1999 USNO data 1998.0 thru 2000.0 (McBurnett) */
|
|---|
| 145 | 6297, 6347, 6383,
|
|---|
| 146 | /* 1999 extrapolation (McBurnett), 2001.0 thru 2006.0
|
|---|
| 147 | Ramp up to 1.6 s/yr to transition smoothly to Stevenson formula
|
|---|
| 148 | in 2130.0 */
|
|---|
| 149 | 6440, 6510, 6600, 6750, 6900, 7060
|
|---|
| 150 |
|
|---|
| 151 | /* original 1997 USNO extrapolation (stern), 1998.0 thru 2004.0
|
|---|
| 152 | 6296, 6420,
|
|---|
| 153 | 6510, 6600, 6700, 6800, 6900 */ /* 7000, 7100, 7200, 7300, 7400, */
|
|---|
| 154 |
|
|---|
| 155 | /* Extrapolated values (USNO) (original Moshier) [1996.0 thru 2005.0]
|
|---|
| 156 | 6183, 6280, 6377, 6475,
|
|---|
| 157 | 6572, 6667, 6765, 6861, 6957
|
|---|
| 158 | */
|
|---|
| 159 | };
|
|---|
| 160 |
|
|---|
| 161 | /* calculate DeltaT = ET - UT1 in seconds. Describes the irregularities
|
|---|
| 162 | * of the Earth rotation rate in the ET time scale.
|
|---|
| 163 | */
|
|---|
| [2551] | 164 | double deltat(double mj)
|
|---|
| [1457] | 165 | {
|
|---|
| 166 | double Y;
|
|---|
| 167 | double p, B;
|
|---|
| 168 | int d[6];
|
|---|
| 169 | int i, iy, k;
|
|---|
| 170 | static double ans;
|
|---|
| [2551] | 171 | static double lastmj = -10000;
|
|---|
| [1457] | 172 |
|
|---|
| [2551] | 173 | if (mj == lastmj) {
|
|---|
| [1457] | 174 | return(ans);
|
|---|
| 175 | }
|
|---|
| [2551] | 176 | lastmj = mj;
|
|---|
| [1457] | 177 |
|
|---|
| [2551] | 178 | Y = 2000.0 + (mj - J2000)/365.25;
|
|---|
| [1457] | 179 |
|
|---|
| 180 | if( Y > TABEND && Y < 2130.0 ) {
|
|---|
| 181 | /* linear interpolation from table end; stern */
|
|---|
| 182 | B = Y - TABEND;
|
|---|
| 183 | ans = dt[TABSIZ-1] + B * (dt[TABSIZ-1] - dt[TABSIZ-2]);
|
|---|
| 184 | ans *= 0.01;
|
|---|
| 185 | return(ans);
|
|---|
| 186 | }
|
|---|
| 187 |
|
|---|
| 188 | if( Y < TABSTART || Y >= 2130.0 ) {
|
|---|
| 189 | if( (Y >= 948.0 - 15.0 && Y < TABSTART) || Y >= 2130.0 ) {
|
|---|
| 190 | /* Stephenson and Morrison, stated domain is 948 to 1600:
|
|---|
| 191 | * 25.5(centuries from 1800)^2 - 1.9159(centuries from 1955)^2
|
|---|
| 192 | * Here we offset by -15 y to minimize the discontinuity,
|
|---|
| 193 | * thus we use it from 933.0 to 1620.0,
|
|---|
| 194 | * and from the end of the table to 2130.0.
|
|---|
| 195 | * f(1620.0) = 60.955200, slope -0.079 s/y
|
|---|
| 196 | * f(2004.0) = 105.649728, slope 1.02 s/y
|
|---|
| 197 | * f(2048.0) = 155.176, slope 1.23 s/y
|
|---|
| 198 | * f(2084.0) = 202.49, slope 1.4 s/y
|
|---|
| 199 | * f(2130.0) = 272, slope .1616
|
|---|
| 200 | * f(2150.0) = 305, slope .17
|
|---|
| 201 | */
|
|---|
| 202 | B = 0.01*(Y - 2000.0);
|
|---|
| 203 | ans = (23.58 * B + 100.3)*B + 101.6;
|
|---|
| 204 | } else {
|
|---|
| 205 | /* Borkowski */
|
|---|
| 206 | /* f(2004.0) = 542.7435, slope 2.65 s/y */
|
|---|
| 207 | B = 0.01*(Y - 2000.0) + 3.75;
|
|---|
| 208 | ans = 35.0 * B * B + 40.;
|
|---|
| 209 | }
|
|---|
| 210 | return(ans);
|
|---|
| 211 | }
|
|---|
| 212 |
|
|---|
| 213 | /* Besselian interpolation from tabulated values.
|
|---|
| 214 | * See AA page K11.
|
|---|
| 215 | */
|
|---|
| 216 |
|
|---|
| 217 | /* value for 1620.1 is 121.96 or so, not 124.0 */
|
|---|
| 218 |
|
|---|
| 219 | /* Index into the table.
|
|---|
| 220 | */
|
|---|
| 221 | p = floor(Y);
|
|---|
| 222 | iy = (int)(p - TABSTART);
|
|---|
| 223 | /* Zeroth order estimate is value at start of year
|
|---|
| 224 | */
|
|---|
| 225 | ans = dt[iy];
|
|---|
| 226 | k = iy + 1;
|
|---|
| 227 | if( k >= TABSIZ )
|
|---|
| 228 | goto done; /* No data, can't go on. */
|
|---|
| 229 |
|
|---|
| 230 | /* The fraction of tabulation interval
|
|---|
| 231 | */
|
|---|
| 232 | p = Y - p;
|
|---|
| 233 |
|
|---|
| 234 | /* First order interpolated value
|
|---|
| 235 | */
|
|---|
| 236 | ans += p*(dt[k] - dt[iy]);
|
|---|
| 237 | if( (iy-1 < 0) || (iy+2 >= TABSIZ) )
|
|---|
| 238 | goto done; /* can't do second differences */
|
|---|
| 239 |
|
|---|
| 240 | /* Make table of first differences
|
|---|
| 241 | */
|
|---|
| 242 | k = iy - 2;
|
|---|
| 243 | for( i=0; i<5; i++ ) {
|
|---|
| 244 | if( (k < 0) || (k+1 >= TABSIZ) )
|
|---|
| 245 | d[i] = 0;
|
|---|
| 246 | else d[i] = dt[k+1] - dt[k];
|
|---|
| 247 | k += 1;
|
|---|
| 248 | }
|
|---|
| 249 |
|
|---|
| 250 | /* Compute second differences
|
|---|
| 251 | */
|
|---|
| 252 | for( i=0; i<4; i++ )
|
|---|
| 253 | d[i] = d[i+1] - d[i];
|
|---|
| 254 | B = 0.25*p*(p-1.0);
|
|---|
| 255 | ans += B*(d[1] + d[2]);
|
|---|
| 256 | if( iy+2 >= TABSIZ )
|
|---|
| 257 | goto done;
|
|---|
| 258 |
|
|---|
| 259 | /* Compute third differences
|
|---|
| 260 | */
|
|---|
| 261 | for( i=0; i<3; i++ )
|
|---|
| 262 | d[i] = d[i+1] - d[i];
|
|---|
| 263 | B = 2.0*B/3.0;
|
|---|
| 264 | ans += (p-0.5)*B*d[1];
|
|---|
| 265 | if( (iy-2 < 0) || (iy+3 > TABSIZ) )
|
|---|
| 266 | goto done;
|
|---|
| 267 |
|
|---|
| 268 | /* Compute fourth differences
|
|---|
| 269 | */
|
|---|
| 270 | for( i=0; i<2; i++ )
|
|---|
| 271 | d[i] = d[i+1] - d[i];
|
|---|
| 272 | B = 0.125*B*(p+1.0)*(p-2.0);
|
|---|
| 273 | ans += B*(d[0] + d[1]);
|
|---|
| 274 |
|
|---|
| 275 | done:
|
|---|
| 276 | /* Astronomical Almanac table is corrected by adding the expression
|
|---|
| 277 | * -0.000091 (ndot + 26)(year-1955)^2 seconds
|
|---|
| 278 | * to entries prior to 1955 (AA page K8), where ndot is the secular
|
|---|
| 279 | * tidal term in the mean motion of the Moon.
|
|---|
| 280 | *
|
|---|
| 281 | * Entries after 1955 are referred to atomic time standards and
|
|---|
| 282 | * are not affected by errors in Lunar or planetary theory.
|
|---|
| 283 | */
|
|---|
| 284 | ans *= 0.01;
|
|---|
| 285 | if( Y < 1955.0 ) {
|
|---|
| 286 | B = (Y - 1955.0);
|
|---|
| 287 | ans += -0.000091 * (-25.8 + 26.0) * B * B;
|
|---|
| 288 | }
|
|---|
| 289 | return( ans );
|
|---|
| 290 | }
|
|---|
| 291 |
|
|---|
| 292 |
|
|---|
| 293 | #ifdef TEST_DT
|
|---|
| 294 | main()
|
|---|
| 295 | {
|
|---|
| 296 | double ans, y;
|
|---|
| 297 |
|
|---|
| 298 | while (scanf("%lf", &y) == 1) {
|
|---|
| 299 | ans = deltat((y - 2000.0)*365.25 + J2000);
|
|---|
| 300 | printf("%.4lf %.4lf\n", y, ans);
|
|---|
| 301 | }
|
|---|
| 302 | }
|
|---|
| 303 | #endif
|
|---|
| 304 |
|
|---|
| 305 | /* For RCS Only -- Do Not Edit */
|
|---|
| [2818] | 306 | static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: deltat.c,v $ $Date: 2005-08-21 10:02:37 $ $Revision: 1.5 $ $Name: not supported by cvs2svn $"};
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
