/* DeltaT = Ephemeris Time - Universal Time
 *
 * The tabulated values of deltaT, in hundredths of a second,
 * were taken from The Astronomical Almanac, page K8.
 * The tabulated range is 1620.0 thru 1992.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 Morrison and Stephenson
 * or Stephenson and Houlden.
 *
 * 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.
 */

/* 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;
extern double dtgiven;


#define DEBUG 0
#define TABSTART 1620.0
#define TABEND 1992.0
#define TABSIZ 373

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, 9600, 9100, 8600, 8000,
7500, 7000, 6600, 6300, 6000, 5800, 5700, 5600, 5600, 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 1988.0 */
5054, 5138, 5217, 5296, 5379, 5434, 5487, 5532, 5582,
/* Extrapolated values 1989.0 thru 1991.0 */
5630, 5670, 5720,
/* Wild guess for 1992 */
5770
};



double deltat(Y)
double Y;
{
double ans, p, B;
int d[6];
int i, iy, k;
double floor();

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 = (Y-1800.0)/100.0  - 0.1;
	ans = -15.0 + 32.5*B*B;
	return(ans);
	}

if( Y < TABSTART )
	{
/* F. R. Stephenson and M. A. Houlden, _Atlas of Historical
 * Eclipse Maps_, Cambridge U. Press (1986)
 * They estimate the uncertainty at 1500 B.C. = 15 minutes.
 */
	if( Y < 948.0 )
		{
		B = (Y - 948.0)/100.0;
		ans = (46.5 * B - 405.) * B + 1830.;
		}
	else
		{ /* stated domain is 948 to 1600 */
		B = (Y - 1850.0)/100.0;
		ans = 22.5 * B * B;
		}
	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]);
#if DEBUG
printf( "B %.4lf, ans %.4lf\n", B, ans );
#endif
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];
#if DEBUG
printf( "B %.4lf, ans %.4lf\n", B*(p-0.5), ans );
#endif
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]);
#if DEBUG
printf( "B %.4lf, ans %.4lf\n", B, ans );
#endif

done:
return(ans/100.0);
}




/* 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.
 */

#include "kep.h"

update()
{
double T;
double deltat();

/* Convert Julian date to Julian years re J2000.0.
 */
T = 2000.0 + (JD - J2000)/365.25;

switch( jdflag )
	{
	case 0:
		TDT = JD;
		UT = JD;
		break;
	case 1:
		TDT = JD;
		UT = TDT - deltat(T)/86400.0;
		jtocal(UT); /* display the other date */
		break;
	case 2:
		UT = JD;
		TDT = UT + deltat(T)/86400.0;
		jtocal(TDT);
		break;
	}
jtocal(JD);
}





/* Exercise program.
 */
#if DEBUG
main()
{
char s[20];
double ans, y;
double deltat();

loop:
printf( "year ? " );
gets(s);
sscanf( s, "%lf", &y );
ans = deltat(y);
printf( "%.4lf\n", ans );
goto loop;
}
#endif