This is deltat.c in view mode; [Download] [Up]
/********************************************************/
/* deltat.c */
/* */
/* computes DeltaT = Ephemeris Time - Universal Time */
/********************************************************/
/***** description
*
* $Id: deltat.c,v 1.6 1993/04/30 18:09:46 craig Exp $
*
* 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 = -23.8946
* arcsec per century squared.
*
* The tabulated range is 1620.0 through 1994.0. Bessel's interpolation
* formula is implemented to obtain fourth order interpolated values at
* intermediate times. Outside 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.
*
* Input Y is the Julian epoch expressed in Julian years. Y can be
* found from the Julian date JD by
* Y = 2000.0 + (JD - 2451545.0)/365.25.
* See AA page B4.
*
* Output double deltat(Y) is ET-UT in seconds.
*
*
* References:
*
* 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)
*/
/***** modification history
*
* $Log: deltat.c,v $
* Revision 1.6 1993/04/30 18:09:46 craig
* Added a output file option to the update () routine.
*
* Revision 1.5 1993/04/27 16:57:46 craig
* Made the jtocal printing optional using the prtflg flag.
*
* Revision 1.4 1993/04/22 19:40:23 craig
* minor changes.
*
* Revision 1.3 1993/04/21 21:07:26 craig
* Changed path for the satellite.h include.
* Changed ecnsts to pcnsts.
*
* Revision 1.2 1993/04/21 15:11:05 craig
* First working version. Ran through indent and converted to ansi.
* Added hooks to work with the satellite programs.
*
*
*/
#define TABSTART 1620.0
#define TABEND 1994.0
#define TABSIZ 375
/***** include files *****/
#include <math.h>
#include "aaproto.h"
#include "satellite.h"
/***** global variables *****/
/* from cnstinit.c */
extern int jdflag;
extern int prtflg;
extern double JD;
extern double TDT;
extern double UT;
extern struct PCONSTANTS pcnsts;
/* If the following number (read from the file aa.ini)
* is nonzero, then the program will return it
* and not calculate anything.
*/
double dtgiven = 0.0;
/***** local global variables *****/
/* 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.
*/
static short dt[TABSIZ] = {
/* 1620.0 thru 1659.0 */
12400, 11900, 11500, 11000, 10600, 10200, 9800, 9500, 9100, 8800,
8500, 8200, 7900, 7700, 7400, 7200, 7000, 6700, 6500, 6300,
6200, 6000, 5800, 5700, 5500, 5400, 5300, 5100, 5000, 4900,
4800, 4700, 4600, 4500, 4400, 4300, 4200, 4100, 4000, 3800,
/* 1660.0 thru 1699.0 */
3700, 3600, 3500, 3400, 3300, 3200, 3100, 3000, 2800, 2700,
2600, 2500, 2400, 2300, 2200, 2100, 2000, 1900, 1800, 1700,
1600, 1500, 1400, 1400, 1300, 1200, 1200, 1100, 1100, 1000,
1000, 1000, 900, 900, 900, 900, 900, 900, 900, 900,
/* 1700.0 thru 1739.0 */
900, 900, 900, 900, 900, 900, 900, 900, 1000, 1000,
1000, 1000, 1000, 1000, 1000, 1000, 1000, 1100, 1100, 1100,
1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100,
1100, 1100, 1100, 1100, 1200, 1200, 1200, 1200, 1200, 1200,
/* 1740.0 thru 1779.0 */
1200, 1200, 1200, 1200, 1300, 1300, 1300, 1300, 1300, 1300,
1300, 1400, 1400, 1400, 1400, 1400, 1400, 1400, 1500, 1500,
1500, 1500, 1500, 1500, 1500, 1600, 1600, 1600, 1600, 1600,
1600, 1600, 1600, 1600, 1600, 1700, 1700, 1700, 1700, 1700,
/* 1780.0 thru 1799.0 */
1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700,
1700, 1700, 1600, 1600, 1600, 1600, 1500, 1500, 1400, 1400,
/* 1800.0 thru 1819.0 */
1370, 1340, 1310, 1290, 1270, 1260, 1250, 1250, 1250, 1250,
1250, 1250, 1250, 1250, 1250, 1250, 1250, 1240, 1230, 1220,
/* 1820.0 thru 1859.0 */
1200, 1170, 1140, 1110, 1060, 1020, 960, 910, 860, 800,
750, 700, 660, 630, 600, 580, 570, 560, 560, 560,
570, 580, 590, 610, 620, 630, 650, 660, 680, 690,
710, 720, 730, 740, 750, 760, 770, 770, 780, 780,
/* 1860.0 thru 1899.0 */
788, 782, 754, 697, 640, 602, 541, 410, 292, 182,
161, 10, -102, -128, -269, -324, -364, -454, -471, -511,
-540, -542, -520, -546, -546, -579, -563, -564, -580, -566,
-587, -601, -619, -664, -644, -647, -609, -576, -466, -374,
/* 1900.0 thru 1939.0 */
-272, -154, -2, 124, 264, 386, 537, 614, 775, 913,
1046, 1153, 1336, 1465, 1601, 1720, 1824, 1906, 2025, 2095,
2116, 2225, 2241, 2303, 2349, 2362, 2386, 2449, 2434, 2408,
2402, 2400, 2387, 2395, 2386, 2393, 2373, 2392, 2396, 2402,
/* 1940.0 thru 1979.0 */
2433, 2483, 2530, 2570, 2624, 2677, 2728, 2778, 2825, 2871,
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 1990.0 */
5054, 5138, 5217, 5296, 5379, 5434, 5487, 5532, 5582, 5630,
5686,
/* Extrapolated values 1991.0 thru 1993.0 */
5740, 5800, 5860,
/* Wild guess for 1994 */
5930
};
/**********/
/* deltat */
/**********/
double deltat (double Y)
{
int d[6];
int i, iy, k;
double ans, p, B;
if (dtgiven != 0.0)
{
return (dtgiven);
}
if (Y > TABEND)
{
/* 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);
}
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. See AA page K11. ***/
/* Index into the table. */
p = floor (Y);
iy = p - TABSTART;
/* Zeroth order estimate is value at start of year */
ans = dt[iy];
k = iy + 1;
if (k >= TABSIZ)
{
goto done; /* No data, can't go on. */
}
/* The fraction of tabulation interval */
p = Y - p;
/* First order interpolated value */
ans += p * (dt[k] - dt[iy]);
if ((iy - 1 < 0) || (iy + 2 >= TABSIZ))
{
goto done; /* can't do second differences */
}
/* Make table of first differences */
k = iy - 2;
for (i = 0; i < 5; i++)
{
if ((k < 0) || (k + 1 >= TABSIZ))
{
d[i] = 0;
}
else
{
d[i] = dt[k + 1] - dt[k];
}
k += 1;
}
/* Compute second differences */
for (i = 0; i < 4; i++)
{
d[i] = d[i + 1] - d[i];
}
B = 0.25 * p * (p - 1.0);
ans += B * (d[1] + d[2]);
/* printf ("B %.4lf, ans %.4lf\n", B, ans); */
if (iy + 2 >= TABSIZ)
{
goto done;
}
/* Compute third differences */
for (i = 0; i < 3; i++)
{
d[i] = d[i + 1] - d[i];
}
B = 2.0 * B / 3.0;
ans += (p - 0.5) * B * d[1];
/* printf ("B %.4lf, ans %.4lf\n", B * (p - 0.5), ans); */
if ((iy - 2 < 0) || (iy + 3 > TABSIZ))
{
goto done;
}
/* Compute fourth differences */
for (i = 0; i < 2; i++)
{
d[i] = d[i + 1] - d[i];
}
B = 0.125 * B * (p + 1.0) * (p - 2.0);
ans += B * (d[0] + d[1]);
/* printf ("B %.4lf, ans %.4lf\n", B, ans); */
done:
/* Astronomical Almanac table is corrected by adding the expression
-1.915914(Julian centuries from 1955)^2 seconds to entries prior to
1955. 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 = 0.01 * (1955.0 - Y);
ans += -1.915914 * B * B;
}
return (ans);
}
/* Routine to fill in time values for TDT and UT.
* These and the input date, JD are global variables.
* jdflag is set up on reading the initialization file aa.ini.
*/
int update (FILE *out)
{
double T;
/* Convert Julian date to Julian years re J2000.0. */
T = 2000.0 + (JD - pcnsts.J2000) / (pcnsts.dapcen / 100.);
switch (jdflag)
{
case 0:
TDT = JD;
UT = JD;
break;
case 1:
TDT = JD;
UT = TDT - deltat (T) / pcnsts.secpda;
if (prtflg)
{
jtocal (out, UT); /* display the other date */
}
break;
case 2:
UT = JD;
TDT = UT + deltat (T) / pcnsts.secpda;
if (prtflg)
{
jtocal (out, TDT);
}
break;
}
if (prtflg)
{
jtocal (out, JD);
}
return (0);
}
/* Exercise program. */
/*
#include <stdio.h>
main ()
{
char s[20];
double ans, y;
loop:
printf ("year ? ");
gets (s);
sscanf (s, "%lf", &y);
ans = deltat (y);
printf ("%.4lf\n", ans);
goto loop;
}
*/
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.