Index: trunk/SophyaExt/XephemAstroLib/deltat.c
===================================================================
--- trunk/SophyaExt/XephemAstroLib/deltat.c	(revision 3654)
+++ trunk/SophyaExt/XephemAstroLib/deltat.c	(revision 4017)
@@ -1,41 +1,32 @@
 /* DeltaT = Ephemeris Time - Universal Time
  *
- * original code by Stephen Moshier <moshier@world.std.com>,
- * adapted to xephem by Michael Sternberg <sternberg@physik.tu-chemnitz.de>
- * smoothed transitions and updated by Neal McBurnett <nealmcb@bell-labs.com>
- *
- **********************************************************************
- *
- * 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 = -25.8
- * arcsec per century squared, the value used in JPL's DE403 ephemeris.
- * ELP2000 (and DE200) used the value -23.8946.
- *
- * The tabulated range is 1620.0 through 1998.0.  Bessel's interpolation
- * formula is implemented to obtain fourth order interpolated values at
- * intermediate times.
- *
- * 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.
- *
- * For future vaues of deltaT, the function smoothly transitions with
- * a linear segment back to Stephenson & Morrison's quadratic formula
- * in 2130.
- *
- * Input is mj (modified julian date from MJD0 on). [stern]
- *  Note that xephem uses a different epoch for this "mj" than the
- *  normal value of JD=240000.5.
- * See AA page B4.
- *
- * Output double deltat(mj) is ET-UT1 in seconds.
+ * Adapted 2011/4/14 from Stephen Moshier <moshier@world.std.com>,
+ * cosmetic changes only.
+ *
+ * Compile as follows to create stand-alone test program:
+ *   cc -DTEST_MAIN deltat.c libastro.a
+ *
+ * Tabulated values of deltaT, in hundredths of a second, are
+ * from The Astronomical Almanac and current IERS reports.
+ * A table of values for the pre-telescopic period was taken from
+ * Morrison and Stephenson (2004).  The overall tabulated range is
+ * -1000.0 through 2011.0.  Values at intermediate times are interpolated
+ * from the tables.
+ *
+ * For dates earlier and later than the tabulated range, the program
+ * calculates a polynomial extrapolation formula.
+ *
+ * Updated deltaT predictions can be obtained from this network archive,
+ *    http://maia.usno.navy.mil
+ * then appended to the dt[] table and update TABEND.
+ *
+ * Input is XEphem's MJD, output is ET-UT in seconds.
  *
  *
  * 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
@@ -44,37 +35,11 @@
  * 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)
  *
- * from obsolete extrapolation code [stern]:
- * Morrison, L. V. and F. R. Stephenson, "Sun and Planetary System"
- * vol 96,73 eds. W. Fricke, G. Teleki, Reidel, Dordrecht (1982)
- *
- **********************************************************************
- *
- * changes by stern:
- *   - adopted #include's for xephem
- *   - made dt[] static
- *   - made mj the time argument [was: year Y].
- *   - updated observed and extrapolated data from tables at
- *	ftp://maia.usno.navy.mil/ser7/ -- data deviated by up to 0.8 s
- *   - removed references to "extern double dtgiven"
- *   - removed DEMO #define and its references
- *   - replaced treatment after TABEND by linear extrapolation instead
- *	of second order version
- *   - installed lastmj cache (made ans static)
- *
- *   - no changes to table interpolation scheme and past extrapolations */
+ */
 
 #include <math.h>
@@ -82,13 +47,35 @@
 #include "astro.h"
 
-#define TABSTART 1620.0
-#define TABEND 2017.0
-#define TABSIZ 398
-
-/* 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] = {
+#define TABSTART 1620
+#define TABEND 2011
+#define TABSIZ (TABEND - TABSTART + 1)
+
+/* 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,
@@ -96,4 +83,5 @@
     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,
@@ -101,4 +89,5 @@
     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,
@@ -106,4 +95,5 @@
     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,
@@ -111,10 +101,13 @@
     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,
@@ -122,4 +115,5 @@
     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,
@@ -127,4 +121,5 @@
     -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,
@@ -132,97 +127,119 @@
     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 1995.0 */
-     5054, 5138, 5217, 5296, 5379, 5434, 5487, 5532, 5582, 5630,
-     5686, 5757, 5831, 5912, 5998, 6078,
-    /* new USNO data (stern) */
-     6163, 6230,
-    /* 1999 USNO data 1998.0 thru 2000.0 (McBurnett) */
-     6297, 6347, 6383, 
-
-     6409, 6430, 6447, 6507, 6546, 6546, 6546,
-
-    /* 2008+ added 5/19/2009 from http://maia.usno.navy.mil/ser7/deltat.preds */
-
-    6546,
-    6570,
-    6650,
-    6710,
-    6800,
-    6800,
-    6900,
-    6900,
-    7000,
-    7000,
+    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 2011.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, 6457, 6469, 6485, 6515, 6546, 6578,
+    6607, 6632,
 };
 
-/* calculate  DeltaT = ET - UT1 in seconds.  Describes the irregularities
+
+/* Given MJD return DeltaT = ET - UT1 in seconds.  Describes the irregularities
  * of the Earth rotation rate in the ET time scale.
  */
-double deltat(double mj)
+double
+deltat(double mj)
 {
-	double Y;
-	double p, B;
+	static double ans, lastmj;
+	double Y, p, B;
 	int d[6];
 	int i, iy, k;
-	static double ans;
-	static double lastmj = -10000;
-
-	if (mj == lastmj) {
-	    return(ans);
+
+	if (mj == lastmj)
+	    return (ans);
+	lastmj = mj;
+
+	mjd_year (mj, &Y);
+
+	if( Y > TABEND ) {
+	    /* 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);
 	}
-	lastmj = mj;
-
-	Y = 2000.0 + (mj - J2000)/365.25;
-
-	if( Y > TABEND) {
-	    /* linear interpolation from table end; stern */
-	    B = Y - TABEND;
-	    ans = dt[TABSIZ-1] + B * (dt[TABSIZ-1]  - dt[TABSIZ-11])/10;
-	    ans *= 0.01;
-	    return(ans);
+
+
+	/* 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);
 	}
 
-	if( Y < TABSTART) {
-	    if( Y >= 948.0 - 15.0 ) {
-		/* Stephenson and Morrison, stated domain is 948 to 1600:
-		 * 25.5(centuries from 1800)^2 - 1.9159(centuries from 1955)^2
-		 * Here we offset by -15 y to minimize the discontinuity,
-		 * thus we use it from 933.0 to 1620.0,
-		 * and from the end of the table to 2130.0.
-		 * f(1620.0) = 60.955200, slope -0.079 s/y
-		 * f(2004.0) = 105.649728, slope 1.02 s/y
-		 * f(2048.0) = 155.176, slope 1.23 s/y
-		 * f(2084.0) = 202.49, slope 1.4 s/y
-		 * f(2130.0) = 272, slope .1616
-	         * f(2150.0) = 305, slope .17
-		 */
-		B = 0.01*(Y - 2000.0);
-		ans = (23.58 * B + 100.3)*B + 101.6;
-	    } else {
-		/* Borkowski */
-	        /* f(2004.0) = 542.7435, slope 2.65 s/y */
-		B = 0.01*(Y - 2000.0)  +  3.75;
-		ans = 35.0 * B * B  +  40.;
-	    }
-	    return(ans);
-	}
-
-	/* Besselian interpolation from tabulated values.
+	/* Besselian interpolation between tabulated values
+	 * in the telescopic era.
 	 * See AA page K11.
 	 */
 
-	/* value for 1620.1 is 121.96 or so, not 124.0 */
-
-	/* Index into the table.
-	 */
+	/* Index into the table.  */
 	p = floor(Y);
-	iy = (int)(p - TABSTART);
-	/* Zeroth order estimate is value at start of year
-	 */
+	iy = (int) (p - TABSTART);
+	/* Zeroth order estimate is value at start of year */
 	ans = dt[iy];
 	k = iy + 1;
@@ -230,44 +247,39 @@
 	    goto done; /* No data, can't go on. */
 
-	/* The fraction of tabulation interval
-	 */
+	/* The fraction of tabulation interval */
 	p = Y - p;
 
-	/* First order interpolated value
-	 */
+	/* 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
-	 */
+	/* Make table of first differences */
 	k = iy - 2;
-	for( i=0; i<5; i++ ) {
+	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;
+	    else
+		d[i] = dt[k+1] - dt[k];
+	    k += 1;
 	}
 
-	/* Compute second differences
-	 */
+	/* 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( iy+2 >= TABSIZ )
+	if (iy+2 >= TABSIZ)
 	    goto done;
 
-	/* Compute third differences
-	 */
+	/* 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( (iy-2 < 0) || (iy+3 > TABSIZ) )
+	if ((iy-2 < 0) || (iy+3 > TABSIZ) )
 	    goto done;
 
-	/* Compute fourth differences
-	 */
+	/* Compute fourth differences */
 	for( i=0; i<2; i++ )
 	    d[i] = d[i+1] - d[i];
@@ -275,5 +287,10 @@
 	ans += B*(d[0] + d[1]);
 
-	done:
+    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
@@ -284,25 +301,34 @@
 	 * are not affected by errors in Lunar or planetary theory.
 	 */
-	ans *= 0.01;
-	if( Y < 1955.0 ) {
-	    B = (Y - 1955.0);
-	    ans += -0.000091 * (-25.8 + 26.0) * B * B;
-	}
+	if( Y < 1955.0 )
+		{
+		B = (Y - 1955.0);
+	#if 1
+		ans += -0.000091 * (-25.8 + 26.0) * B * B;
+	#else
+		ans += -0.000091 * (-23.8946 + 26.0) * B * B;
+	#endif
+		}
+
+#endif /* 0 */
+
 	return( ans );
 }
 
 
-#ifdef TEST_DT
-main()
+#ifdef TEST_MAIN
+
+/* Exercise program.
+ */
+#include <stdio.h>
+#include <stdlib.h>
+
+int main(int ac, char *av[])
 {
-	double ans, y;
-
-	while (scanf("%lf", &y) == 1) {
-		ans = deltat((y - 2000.0)*365.25 + J2000);
-		printf("%.4lf %.4lf\n", y, ans);
-	}
+	double ans, mj, y = atof(av[1]);
+	year_mjd (y, &mj);
+	ans = deltat(mj);
+	printf( "%.4lf\n", ans );
+	return (0);
 }
 #endif
-
-/* For RCS Only -- Do Not Edit */
-static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: deltat.c,v $ $Date: 2009-07-16 10:34:37 $ $Revision: 1.8 $ $Name: not supported by cvs2svn $"};
