Fixed Stephenson & Houlden (1986) algorithm for DeltaT

- Added 2 new unit tests
- Added 2 TODO's
- Updated description for S&H solution
- Updated valid ranges
- TODO: Update related parts of SUG

5 files changed, 180 insertions, 54 deletions

diff --git a/src/core/StelCore.cpp b/src/core/StelCore.cpp
index c240da6864..689cc821af 100644
--- a/src/core/StelCore.cpp
+++ b/src/core/StelCore.cpp
@@ -2075,14 +2075,14 @@ void StelCore::setCurrentDeltaTAlgorithm(DeltaTAlgorithm algorithm)
 			deltaTnDot = -26.0; // n.dot = -26.0 "/cy/cy
 			deltaTfunc = StelUtils::getDeltaTByStephensonMorrison1984;
 			deltaTstart	= -391;
-			deltaTfinish	= 1600;
+			deltaTfinish	= 1600; // TODO: 1630..1980 are tabulated values
 			break;
 		case StephensonHoulden:
-			// Stephenson & Houlden (1986) algorithm for DeltaT. The limits are implicitly given by the tabulated values.
+			// Stephenson & Houlden (1986) algorithm for DeltaT.
 			deltaTnDot = -26.0; // n.dot = -26.0 "/cy/cy
 			deltaTfunc = StelUtils::getDeltaTByStephensonHoulden;
-			deltaTstart	= -600;
-			deltaTfinish	= 1650;
+			deltaTstart	= -1500;
+			deltaTfinish	= 1600; // TODO: 1630..1980 are tabulated values from Stephenson & Morrison (1984)
 			break;
 		case Espenak:
 			// Espenak (1987, 1989) algorithm for DeltaT
@@ -2321,7 +2321,7 @@ QString StelCore::getCurrentDeltaTAlgorithmDescription(void) const
 			description = q_("This formula was published by F. R. Stephenson and L. V. Morrison in the article <em>Long-term changes in the rotation of the earth - 700 B.C. to A.D. 1980</em> (%1).").arg("<a href='http://adsabs.harvard.edu/abs/1984RSPTA.313...47S'>1984</a>").append(getCurrentDeltaTAlgorithmValidRangeDescription(jd, &marker));
 			break;
 		case StephensonHoulden:
-			description = q_("This algorithm is used in the PC planetarium program Guide 7.").append(getCurrentDeltaTAlgorithmValidRangeDescription(jd, &marker));
+			description = q_("This algorithm was published by F. R. Stephenson and M. A. Houlden in the book <em>Atlas of Historical Eclipse Maps</em> (%1). This algorithm is used in the PC planetarium program Guide 7.").arg("<a href='https://assets.cambridge.org/97805212/67236/frontmatter/9780521267236_frontmatter.pdf'>1986</a>").append(getCurrentDeltaTAlgorithmValidRangeDescription(jd, &marker));
 			break;
 		case Espenak: // limited range, but wide availability?
 			description = q_("This algorithm was given by F. Espenak in his <em>Fifty Year Canon of Solar Eclipses: 1986-2035</em> (1987) and in his <em>Fifty Year Canon of Lunar Eclipses: 1986-2035</em> (1989).").append(getCurrentDeltaTAlgorithmValidRangeDescription(jd, &marker));
diff --git a/src/core/StelUtils.cpp b/src/core/StelUtils.cpp
index d2e8d04e51..ede71d9cc6 100644
--- a/src/core/StelUtils.cpp
+++ b/src/core/StelUtils.cpp
@@ -1205,7 +1205,7 @@ QString hoursToHmsStr(const float hours, const bool lowprecision)
 //  1820.0=1820-jan-0.5=2385800.0
 //  1810.0=1810-jan-0.5=2382148.0
 //  1800.0=1800-jan-0.5=2378496.0
-//  1735.0=1735-jan-0.5=2354755.0
+//  1735.0=1735-jan-0.5=2354756.0
 //  1625.0=1625-jan-0.5=2314579.0
 //
 // Up to V0.15.1, if the requested year was outside validity range, we returned zero or some useless value.
@@ -1381,7 +1381,7 @@ double getDeltaTByStephenson1978(const double jDay)
 // Implementation of algorithm by Stephenson (1997) for DeltaT computation
 double getDeltaTByStephenson1997(const double jDay)
 {
-	double u=(jDay-2354755.0)/36525.0; // (1735-jan-0.5)
+	double u=(jDay-2354756.0)/36525.0; // (1735-jan-0.5)
 	return -20.0 + 35.0*u*u;
 }
 
@@ -1432,9 +1432,27 @@ double getDeltaTByStephensonMorrison1995(const double jDay)
 // Implementation of algorithm by Stephenson & Houlden (1986) for DeltaT computation
 double getDeltaTByStephensonHoulden(const double jDay)
 {
+	int year, month, day;
+	double u, deltaT = 0.;
+	getDateFromJulianDay(jDay, &year, &month, &day);
+
+	if (year <= 948)
+	{
+		//u = (year-948)/100;
+		u = (yearFraction(year, month, day)-948)/100;
+		deltaT = (46.5*u -405.0)*u +1830.0;
+	}
+	if (948 < year && year <= 1600)
+	{
+		//u = (year-1850)/100;
+		u = (yearFraction(year, month, day)-1850)/100;
+		deltaT = 22.5*u*u;
+	}
+
+	return deltaT;
 	// This formula found in the cited book, page (ii), formula (1).
-	double T=(jDay-2415020.0)/36525; // centuries from J1900.0
-	return (36.79*T+35.06)*T+4.87;
+	//double T=(jDay-2415020.0)/36525; // centuries from J1900.0
+	//return (36.79*T+35.06)*T+4.87;
 }
 
 // Implementation of algorithm by Espenak (1987, 1989) for DeltaT computation
diff --git a/src/core/StelUtils.hpp b/src/core/StelUtils.hpp
index fe215b67d2..9de4dc22d2 100644
--- a/src/core/StelUtils.hpp
+++ b/src/core/StelUtils.hpp
@@ -661,6 +661,8 @@ namespace StelUtils
 
 	//! Get Delta-T estimation for a given date.
 	//! Implementation of algorithm by Stephenson & Houlden (1986) for DeltaT computation
+	//! Source: STEPHENSON F.R and HOULDEN M.A. - Atlas of Historical Eclipse Maps - Cambridge Univ.Press. (1986)
+	//! [https://assets.cambridge.org/97805212/67236/frontmatter/9780521267236_frontmatter.pdf]
 	//! @param jDay the date and time expressed as a Julian day
 	//! @return Delta-T in seconds or 0 if year > 1600
 	double getDeltaTByStephensonHoulden(const double jDay);
diff --git a/src/tests/testDeltaT.cpp b/src/tests/testDeltaT.cpp
index c19da5a8dc..97b735bf62 100644
--- a/src/tests/testDeltaT.cpp
+++ b/src/tests/testDeltaT.cpp
@@ -604,61 +604,33 @@ void TestDeltaT::testDeltaTByStephensonMorrison1995WideDates()
 	}
 }
 
-void TestDeltaT::testDeltaTByStephenson1997WideDates()
+void TestDeltaT::testDeltaTByStephenson1997()
 {
-	// test data taken from http://www.staff.science.uu.nl/~gent0113/deltat/deltat_old.htm
+	// test data taken from http://ebooks.cambridge.org/ebook.jsf?bid=CBO9780511525186
 
 	QVariantList data;
-
-	/*
-	// FIXME: WTF?
-	data << -500 << 16800;
-	data << -450 << 16000;
-	data << -400 << 15300;
-	data << -350 << 14600;
-	data << -300 << 14000;
-	data << -250 << 13400;
-	data << -200 << 12800;
-	data << -150 << 12200;
-	data << -100 << 11600;
-	data <<  -50 << 11100;
-	*/
+	//data << -500 << 16800;
+	//data << -400 << 15300;
+	//data << -300 << 14000;
+	//data << -200 << 12800;
+	//data << -100 << 11600;
 	data <<    0 << 10600;
-	data <<   50 << 10100;
 	data <<  100 <<  9600;
-	data <<  150 <<  9100;
 	data <<  200 <<  8600;
-	data <<  250 <<  8200;
 	data <<  300 <<  7700;
-	data <<  350 <<  7200;
 	data <<  400 <<  6700;
-	data <<  450 <<  6200;
 	data <<  500 <<  5700;
-	data <<  550 <<  5200;
 	data <<  600 <<  4700;
-	data <<  650 <<  4300;
 	data <<  700 <<  3800;
-	data <<  750 <<  3400;
-	/*
-	// FIXME: WTF?
-	data <<  800 <<  3000;
-	data <<  850 <<  2600;
-	data <<  900 <<  2200;
-	data <<  950 <<  1900;
-	data << 1000 <<  1600;
-	data << 1050 <<  1350;
-	data << 1100 <<  1100;
-	data << 1150 <<   900;
-	data << 1200 <<   750;
-	data << 1250 <<   600;
-	data << 1300 <<   460;
-	data << 1350 <<   360;
-	data << 1400 <<   280; // 300
-	data << 1450 <<   200; // 230
-	data << 1500 <<   150; // 180
-	data << 1550 <<   110; // 140
-	data << 1600 <<    80;  // 110
-	*/
+	//data <<  800 <<  3000;
+	//data <<  900 <<  2200;
+	//data << 1000 <<  1600;
+	//data << 1100 <<  1100;
+	//data << 1200 <<   750;
+	//data << 1300 <<   470;
+	//data << 1400 <<   300;
+	//data << 1500 <<   180;
+	//data << 1600 <<   110;
 
 	while(data.count() >= 2)
 	{
@@ -1322,6 +1294,138 @@ void TestDeltaT::testDeltaTByChaprontTouze()
 	}
 }
 
+void TestDeltaT::testDeltaTByIAUWideDates()
+{
+	// Source of test data: https://web.archive.org/web/20150226203358/http://user.online.be/felixverbelen/dt.htm
+	QVariantList data;
+	data << -2000 << 42757.0;
+	data << -1900 << 40524.0;
+	data << -1800 << 38350.0;
+	data << -1700 << 36236.0;
+	data << -1600 << 34181.0;
+	data << -1500 << 32187.0;
+	data << -1400 << 30253.0;
+	data << -1300 << 28378.0;
+	data << -1200 << 26564.0;
+	data << -1100 << 24809.0;
+	data << -1000 << 23115.0;
+	data <<  -900 << 21480.0;
+	data <<  -800 << 19905.0;
+	data <<  -700 << 18390.0;
+	data <<  -600 << 16935.0;
+	data <<  -500 << 15539.0;
+	data <<  -400 << 14204.0;
+	data <<  -300 << 12929.0;
+	data <<  -200 << 11713.0;
+	data <<  -100 << 10557.0;
+	data <<     0 <<  9462.0;
+	data <<   100 <<  8426.0;
+	data <<   200 <<  7450.0;
+	data <<   300 <<  6534.0;
+	data <<   400 <<  5678.0;
+	data <<   500 <<  4882.0;
+	data <<   600 <<  4145.0;
+	data <<   700 <<  3469.0;
+	data <<   800 <<  2852.0;
+	data <<   900 <<  2296.0;
+	data <<  1000 <<  1799.0;
+	data <<  1100 <<  1362.0;
+	data <<  1200 <<   985.0;
+	data <<  1300 <<   668.0;
+	data <<  1400 <<   411.0;
+	data <<  1500 <<   214.0;
+	data <<  1600 <<    76.0;
+	data <<  1700 <<    -1.0;
+	//data <<  1800 <<   -19.0;
+	//data <<  1900 <<    24.0;
+	data <<  2000 <<   126.0;
+
+	while(data.count() >= 2)
+	{
+		int year = data.takeFirst().toInt();
+		int yout, mout, dout;
+		double JD;
+		double expectedResult = data.takeFirst().toDouble();
+		double acceptableError = 1.0; // TODO: Increase accuracy to 0.1 seconds
+		StelUtils::getJDFromDate(&JD, year, 1, 1, 0, 0, 0);
+		double result = StelUtils::getDeltaTByIAU(JD);
+		double actualError = qAbs(qAbs(expectedResult) - qAbs(result));
+		StelUtils::getDateFromJulianDay(JD, &yout, &mout, &dout);
+		QVERIFY2(actualError <= acceptableError, QString("date=%2 year=%3 result=%4 expected=%5 error=%6 acceptable=%7")
+							.arg(QString("%1-%2-%3 00:00:00").arg(yout).arg(mout).arg(dout))
+							.arg(year)
+							.arg(result)
+							.arg(expectedResult)
+							.arg(actualError)
+							.arg(acceptableError)
+							.toUtf8());
+	}
+}
+
+void TestDeltaT::testDeltaTByStephensonHoulden()
+{
+	// Source of test data: https://web.archive.org/web/20150226203358/http://user.online.be/felixverbelen/dt.htm
+	QVariantList data;
+	data << -2000 << 54181.0;
+	data << -1900 << 51081.0;
+	data << -1800 << 48073.0;
+	data << -1700 << 45159.0;
+	data << -1600 << 42338.0;
+	data << -1500 << 39610.0; // begin of valid range
+	data << -1400 << 36975.0;
+	data << -1300 << 34433.0;
+	data << -1200 << 31984.0;
+	data << -1100 << 29627.0;
+	data << -1000 << 27364.0;
+	data <<  -900 << 25194.0;
+	data <<  -800 << 23117.0;
+	data <<  -700 << 21133.0;
+	data <<  -600 << 19242.0;
+	data <<  -500 << 17444.0;
+	data <<  -400 << 15738.0;
+	data <<  -300 << 14126.0;
+	data <<  -200 << 12607.0;
+	data <<  -100 << 11181.0;
+	data <<     0 <<  9848.0;
+	data <<   100 <<  8608.0;
+	data <<   200 <<  7461.0;
+	data <<   300 <<  6406.0;
+	data <<   400 <<  5445.0;
+	data <<   500 <<  4577.0;
+	data <<   600 <<  3802.0;
+	data <<   700 <<  3120.0;
+	data <<   800 <<  2531.0;
+	data <<   900 <<  2035.0;
+	data <<  1000 <<  1625.0;
+	data <<  1100 <<  1265.0;
+	data <<  1200 <<   950.0;
+	data <<  1300 <<   680.0;
+	data <<  1400 <<   455.0;
+	data <<  1500 <<   275.0;
+	data <<  1600 <<   140.0;
+
+	while(data.count() >= 2)
+	{
+		int year = data.takeFirst().toInt();
+		int yout, mout, dout;
+		double JD;
+		double expectedResult = data.takeFirst().toDouble();
+		double acceptableError = 1.0; // 1 sec
+		StelUtils::getJDFromDate(&JD, year, 1, 1, 0, 0, 0);
+		double result = StelUtils::getDeltaTByStephensonHoulden(JD);
+		double actualError = qAbs(qAbs(expectedResult) - qAbs(result));
+		StelUtils::getDateFromJulianDay(JD, &yout, &mout, &dout);
+		QVERIFY2(actualError <= acceptableError, QString("date=%2 year=%3 result=%4 expected=%5 error=%6 acceptable=%7")
+							.arg(QString("%1-%2-%3 00:00:00").arg(yout).arg(mout).arg(dout))
+							.arg(year)
+							.arg(result)
+							.arg(expectedResult)
+							.arg(actualError)
+							.arg(acceptableError)
+							.toUtf8());
+	}
+}
+
 void TestDeltaT::testDeltaTStandardError()
 {
 	// the test data was obtained from https://eclipse.gsfc.nasa.gov/SEhelp/uncertainty2004.html
diff --git a/src/tests/testDeltaT.hpp b/src/tests/testDeltaT.hpp
index f078d79440..07649c41b2 100644
--- a/src/tests/testDeltaT.hpp
+++ b/src/tests/testDeltaT.hpp
@@ -40,7 +40,7 @@ private slots:
 	void testDeltaTByStephensonMorrison2004WideDates();	
 	void testDeltaTByStephensonMorrisonHohenkerk2016GenericDates();
 	void testDeltaTByStephensonMorrisonHohenkerk2016SpecialDates();
-	void testDeltaTByStephenson1997WideDates();	
+	void testDeltaTByStephenson1997();
 	void testDeltaTByMeeusSimons();	
 	void testDeltaTByKhalidSultanaZaidiWideDates();
 	void testDeltaTByMontenbruckPfleger();
@@ -51,6 +51,8 @@ private slots:
 	void testDeltaTByAstronomicalEphemeris();
 	void testDeltaTByTuckermanGoldstine();
 	void testDeltaTByChaprontTouze();
+	void testDeltaTByIAUWideDates();
+	void testDeltaTByStephensonHoulden();
 	void testDeltaTStandardError();
 private:
 	QVariantList genericData;