Modified DeltaT (#2340)

* Modified DeltaT

Current value of default DeltaT (by Espenak & Meeus 2006) in Stellarium is increasingly different from observed values published by IERS. New formulae are developed to match observed and predicted values.

* Separate original & modified

* SUG update

* small adjustment & comments

7 files changed, 149 insertions, 15 deletions

diff --git a/guide/ch_astronomical_concepts.tex b/guide/ch_astronomical_concepts.tex
index c02a16fda3..dc1dad35dd 100644
--- a/guide/ch_astronomical_concepts.tex
+++ b/guide/ch_astronomical_concepts.tex
@@ -593,10 +593,11 @@ temporal offset which amounts already to more than a minute in the 21st
 century. There have been many attempts to properly model $\Delta T$,
 and Stellarium offers several models you can choose from in the
 configuration dialog (see
-section~\ref{sec:gui:configuration:time}). The default, ``Espenak and
-Meeus (2006)'', is a widely accepted standard. But if you are a
-researcher and want to experiment with alternative models, you will
-hopefully like this feature. you can even specify your own data for
+section~\ref{sec:gui:configuration:time}). ``Espenak and
+Meeus (2006)'', is a widely accepted standard. The default, ``Modified Espenak and
+Meeus (2006)'', is created to make $\Delta T$ closely in line with observed and predicted
+values between the years 2005-2050. But if you are a researcher and want to experiment
+with alternative models, you will hopefully like this feature. you can even specify your own data for
 $a$, $b$, $c$, $y$ and the secular term for lunar acceleration $n$
 (actually $\dot{n}=dn/dt$ in units of
 $\mathrm{arcseconds}/\mathrm{century}^2$) if you can model $\Delta T$
@@ -746,8 +747,7 @@ The following list describes sources and a few details about the models for $\De
   values for 2016 and later deviate notably from current estimates, and should not be used for dates after 2015).
   Recommended range of usage: between years -720.0 and 2016.0 $\dot{n}=-25.82''/\cy^2$.
   
-\item[Espenak \& Meeus (2006).] This solution\footnote{This solution
-    is used by default.} by F. Espenak and J. Meeus, based on 
+\item[Espenak \& Meeus (2006).] This solution by F. Espenak and J. Meeus, based on 
 	\citet{2004JHA....35..327M} and a polynomial fit
   through tabulated values for 1600-2000, is used for the NASA Eclipse
   Web Site\footnote{\url{https://eclipse.gsfc.nasa.gov/eclipse.html}}
@@ -756,6 +756,11 @@ The following list describes sources and a few details about the models for $\De
   range of usage: between years -1999 and
   3000. $\dot{n}=-25.858''/\cy^2$.
 
+\item[Modified Espenak \& Meeus (2006).] This solution\footnote{This solution
+  is used by default.} is a modification from Espenak \& Meeus (2006). Two formulae are developed to make $\Delta T$
+  closely in line with observed and predicted values by IERS and also connected with original values.
+  Valid range of usage: between years -1999 and 3000. $\dot{n}=-25.858''/\cy^2$.
+
 \item[Reijs (2006).] From the Length of Day (LOD; as determined by
   \citet{2004JHA....35..327M}), Victor Reijs
   derived a $\Delta T$ formula by using a Simplex optimization with a
diff --git a/guide/ch_interface.tex b/guide/ch_interface.tex
index 1c1d14a256..d7a9c58e82 100644
--- a/guide/ch_interface.tex
+++ b/guide/ch_interface.tex
@@ -290,7 +290,7 @@ when the program starts:
 
 \noindent The lowest field allows selection of the correction model for the time
 correction $\Delta T$ (see section~\ref{sec:Concepts:DeltaT}). Default
-is ``Espenak and Meeus (2006)''. Please use other values only if you
+is ``Modified Espenak and Meeus (2006)''. Please use other values only if you
 know what you are doing.
 
 \subsection{The Tools Tab}
diff --git a/src/core/StelCore.cpp b/src/core/StelCore.cpp
index 21e50cffb5..b9932a72e0 100644
--- a/src/core/StelCore.cpp
+++ b/src/core/StelCore.cpp
@@ -207,7 +207,7 @@ void StelCore::init()
 
 	// Delta-T stuff
 	// Define default algorithm for time correction (Delta T)
-	QString tmpDT = conf->value("navigation/time_correction_algorithm", "EspenakMeeus").toString();
+	QString tmpDT = conf->value("navigation/time_correction_algorithm", "EspenakMeeusModified").toString();
 	setCurrentDeltaTAlgorithmKey(tmpDT);
 
 	// Define variables of custom equation for calculation of Delta T
@@ -2199,6 +2199,13 @@ void StelCore::setCurrentDeltaTAlgorithm(DeltaTAlgorithm algorithm)
 			deltaTstart	= -1999;
 			deltaTfinish	= 3000;
 			break;
+		case EspenakMeeusModified:
+			// Espenak & Meeus (2006) algorithm (with modified formulae) for DeltaT
+			deltaTnDot = -25.858; // n.dot = -25.858 "/cy/cy
+			deltaTfunc = StelUtils::getDeltaTByEspenakMeeusModified;
+			deltaTstart	= -1999;
+			deltaTfinish	= 3000;
+			break;
 		case EspenakMeeusZeroMoonAccel:
 			// This is a trying area. Something is wrong with DeltaT, maybe ndot is still not applied correctly.
 			// Espenak & Meeus (2006) algorithm for DeltaT
@@ -2372,7 +2379,10 @@ QString StelCore::getCurrentDeltaTAlgorithmDescription(void) const
 			description = q_("From the Length of Day (LOD; as determined by Stephenson & Morrison (%2)), Victor Reijs derived a %1T formula by using a Simplex optimisation with a cosine and square function. This is based on a possible periodicy described by Stephenson (%2). See for more info %3here%4.").arg(QChar(0x0394)).arg("<a href='http://adsabs.harvard.edu/abs/2004JHA....35..327M'>2004</a>").arg("<a href='http://www.iol.ie/~geniet/eng/DeltaTeval.htm'>").arg("</a>").append(getCurrentDeltaTAlgorithmValidRangeDescription(jd, &marker));
 			break;
 		case EspenakMeeus: // GENERAL SOLUTION
-			description = q_("This solution by F. Espenak and J. Meeus, based on Morrison & Stephenson (2004) and a polynomial fit through tabulated values for 1600-2000, is used for the %1NASA Eclipse Web Site%2 and in their <em>Five Millennium Canon of Solar Eclipses: -1900 to +3000</em> (2006). This formula is also used in the solar, lunar and planetary ephemeris program SOLEX.").arg("<a href='http://eclipse.gsfc.nasa.gov/eclipse.html'>").arg("</a>").append(getCurrentDeltaTAlgorithmValidRangeDescription(jd, &marker)).append(" <em>").append(q_("Used by default.")).append("</em>");
+			description = q_("This solution by F. Espenak and J. Meeus, based on Morrison & Stephenson (2004) and a polynomial fit through tabulated values for 1600-2000, is used for the %1NASA Eclipse Web Site%2 and in their <em>Five Millennium Canon of Solar Eclipses: -1900 to +3000</em> (2006). This formula is also used in the solar, lunar and planetary ephemeris program SOLEX.").arg("<a href='http://eclipse.gsfc.nasa.gov/eclipse.html'>").arg("</a>").append(getCurrentDeltaTAlgorithmValidRangeDescription(jd, &marker));
+			break;
+		case EspenakMeeusModified: // MODIFIED SOLUTION
+			description = q_("This solution is modified from F. Espenak and J. Meeus, based on Morrison & Stephenson (2004) and a polynomial fit through tabulated values for 1600-2000, is used for the %1NASA Eclipse Web Site%2 and in their <em>Five Millennium Canon of Solar Eclipses: -1900 to +3000</em> (2006). This formula is also used in the solar, lunar and planetary ephemeris program SOLEX. Formula for 2005-2050 is modified to match observed values and near-term predictions.").arg("<a href='http://eclipse.gsfc.nasa.gov/eclipse.html'>").arg("</a>").append(getCurrentDeltaTAlgorithmValidRangeDescription(jd, &marker)).append(" <em>").append(q_("Used by default.")).append("</em>");
 			break;
 		case EspenakMeeusZeroMoonAccel: // PATCHED SOLUTION. Experimental, it may not make sense to keep it in V1.0.
 			description = QString("%1 %2").arg(q_("PATCHED VERSION WITHOUT ADDITIONAL LUNAR ACCELERATION.")).arg(q_("This solution by F. Espenak and J. Meeus, based on Morrison & Stephenson (2004) and a polynomial fit through tabulated values for 1600-2000, is used for the %1NASA Eclipse Web Site%2 and in their <em>Five Millennium Canon of Solar Eclipses: -1900 to +3000</em> (2006). This formula is also used in the solar, lunar and planetary ephemeris program SOLEX.").arg("<a href='http://eclipse.gsfc.nasa.gov/eclipse.html'>").arg("</a>").append(getCurrentDeltaTAlgorithmValidRangeDescription(jd, &marker)).append(" <em>").append("</em>"));
@@ -2440,7 +2450,8 @@ QString StelCore::getCurrentDeltaTAlgorithmValidRangeDescription(const double JD
 		case ReingoldDershowitz:     // and
 		case MorrisonStephenson2004: // and
 		case Reijs:                  // and
-		case EspenakMeeus: // the default, range stated in the Canon, p. 14.  ... and
+		case EspenakMeeus:           // range stated in the Canon, p. 14.  ... and
+		case EspenakMeeusModified:   // the default, range stated in the Canon, p. 14.  ... and
 		case EspenakMeeusZeroMoonAccel: // and
 		case StephensonMorrisonHohenkerk2016: // and
 		case Henriksson2017:
diff --git a/src/core/StelCore.hpp b/src/core/StelCore.hpp
index a46c2eadf6..858453e2bb 100644
--- a/src/core/StelCore.hpp
+++ b/src/core/StelCore.hpp
@@ -140,7 +140,8 @@ public:
 		ReingoldDershowitz,                 //!< Reingold & Dershowitz (2002, 2007) algorithm for DeltaT
 		MorrisonStephenson2004,	//!< Morrison & Stephenson (2004, 2005) algorithm for DeltaT
 		Reijs,					//!< Reijs (2006) algorithm for DeltaT
-		EspenakMeeus,				//!< Espenak & Meeus (2006) algorithm for DeltaT (Recommended, default)
+		EspenakMeeus,				//!< Espenak & Meeus (2006) algorithm for DeltaT
+		EspenakMeeusModified,				//!< Espenak & Meeus (2006) algorithm with modified formulae for DeltaT (Recommended, default)
 		EspenakMeeusZeroMoonAccel,	//!< Espenak & Meeus (2006) algorithm for DeltaT (but without additional Lunar acceleration. FOR TESTING ONLY, NONPUBLIC)
 		Banjevic,					//!< Banjevic (2006) algorithm for DeltaT
 		IslamSadiqQureshi,			//!< Islam, Sadiq & Qureshi (2008 + revisited 2013) algorithm for DeltaT (6 polynomials)
diff --git a/src/core/StelUtils.cpp b/src/core/StelUtils.cpp
index 16dacaef47..4d0b75b707 100644
--- a/src/core/StelUtils.cpp
+++ b/src/core/StelUtils.cpp
@@ -1375,6 +1375,113 @@ double getDeltaTByEspenakMeeus(const double jDay)
 	return r;
 }
 
+// Implementation of algorithm by Espenak & Meeus (2006) with modified formulae for DeltaT computation
+double getDeltaTByEspenakMeeusModified(const double jDay)
+{
+	int year, month, day;	
+	getDateFromJulianDay(jDay, &year, &month, &day);
+
+	// Note: the method here is adapted from
+	// "Five Millennium Canon of Solar Eclipses" [Espenak and Meeus, 2006]
+	// A summary is described here:
+	// http://eclipse.gsfc.nasa.gov/SEhelp/deltatpoly2004.html
+
+	double y = yearFraction(year, month, day);
+
+	// set the default value for Delta T
+	double u = (y-1820)/100.;
+	double r = (-20 + 32 * u * u);
+
+	if (y < -500)
+	{
+		// values are equal to defaults!
+	}
+	else if (y < 500)
+	{
+		u = y/100;
+		r = (((((0.0090316521*u + 0.022174192)*u - 0.1798452)*u - 5.952053)*u + 33.78311)*u -1014.41)*u +10583.6;
+	}
+	else if (y < 1600)
+	{
+		u = (y-1000)/100;
+		r = (((((0.0083572073*u - 0.005050998)*u -0.8503463)*u +0.319781)*u + 71.23472)*u -556.01)*u + 1574.2;
+	}
+	else if (y < 1700)
+	{
+		double t = y - 1600;
+		r = ((t/7129.0 - 0.01532)*t - 0.9808)*t +120.0;
+	}
+	else if (y < 1800)
+	{
+		double t = y - 1700;
+		r = (((-t/1174000.0 + 0.00013336)*t - 0.0059285)*t + 0.1603)*t +8.83;
+	}
+	else if (y < 1860)
+	{
+		double t = y - 1800;
+		r = ((((((.000000000875*t -.0000001699)*t + 0.0000121272)*t - 0.00037436)*t + 0.0041116)*t + 0.0068612)*t - 0.332447)*t +13.72;
+	}
+	else if (y < 1900)
+	{
+		double t = y - 1860;
+		r = ((((t/233174.0 -0.0004473624)*t + 0.01680668)*t - 0.251754)*t + 0.5737)*t + 7.62;
+	}
+	else if (y < 1920)
+	{
+		double t = y - 1900;
+		r = (((-0.000197*t + 0.0061966)*t - 0.0598939)*t + 1.494119)*t -2.79;
+	}
+	else if (y < 1941)
+	{
+		double t = y - 1920;
+		r = ((0.0020936*t - 0.076100)*t+ 0.84493)*t +21.20;
+	}
+	else if (y < 1961)
+	{
+		double t = y - 1950;
+		r = ((t/2547.0 -1.0/233.0)*t + 0.407)*t +29.07;
+	}
+	else if (y < 1986)
+	{
+		double t = y - 1975;
+		r = ((-t/718.0 -1/260.0)*t + 1.067)*t +45.45;
+	}
+	else if (y < 2005)
+	{
+		double t = y - 2000;
+		r = ((((0.00002373599*t + 0.000651814)*t + 0.0017275)*t - 0.060374)*t + 0.3345)*t +63.86;
+	}
+	// WB: Polynomial data fit from IERS's DeltaT values during 2005-2021, including predicted values until 2032
+	// Data: https://cddis.nasa.gov/archive/products/iers/deltat.data & https://cddis.nasa.gov/archive/products/iers/deltat.preds
+	// Last updated: 2022 Mar 11
+	else if (y < 2032)
+	{
+		double t = y - 2000;
+		r = (-0.00331233402*t + 0.404229283)*t + 62.48;
+	}
+	// WB: Formula to create reasonable curve between final predicted year and 2050
+	// 93 is the predicted deltaT for 2050 (Espenak/Meeus)
+	// It looks more likely that this value is too high
+	// But we will follow it for now until we have a better model
+	// Small corrections (<0.2 sec) related to secular acceleration of the Moon are neglected
+	// Last updated: 2022 Mar 11
+	else if (y < 2050)
+	{
+		double finalPredictedYear = 2032.;
+		double finalPredictedDeltaT = 72.07;
+		double t = y - finalPredictedYear;
+		double diff = 2050.-finalPredictedYear;
+		r = finalPredictedDeltaT + (93.-finalPredictedDeltaT) * t * t/(diff*diff);
+	}
+	else if (y < 2150)
+	{
+		// r has been precomputed before, just add the term patching the discontinuity
+		r -= 0.5628*(2150.0-y);
+	}
+
+	return r;
+}
+
 // Implementation of algorithm by Schoch (1931) for DeltaT computation
 double getDeltaTBySchoch(const double jDay)
 {
diff --git a/src/core/StelUtils.hpp b/src/core/StelUtils.hpp
index 04cd75ea92..4820ac5413 100644
--- a/src/core/StelUtils.hpp
+++ b/src/core/StelUtils.hpp
@@ -560,6 +560,14 @@ namespace StelUtils
 	double getDeltaTByEspenakMeeus(const double jDay);
 
 	//! Get Delta-T estimation for a given date.
+	//! Note that this method is recommended for the year range:
+	//! -1999 to +3000. It gives details for -500...+2150.
+	//! Implementation of algorithm by Espenak & Meeus (2006) with modified formulae for DeltaT computation
+	//! @param jDay the date and time expressed as a Julian day
+	//! @return Delta-T in seconds
+	double getDeltaTByEspenakMeeusModified(const double jDay);
+
+	//! Get Delta-T estimation for a given date.
 	//! Implementation of algorithm by Schoch (1931) for DeltaT computation,
 	//! outdated but may be useful for science-historical purposes.
 	//! Source: Schoch, C. (1931). Die sekulare Accelaration des Mondes und der Sonne.
diff --git a/src/gui/ConfigurationDialog.cpp b/src/gui/ConfigurationDialog.cpp
index e769125278..0db13385eb 100644
--- a/src/gui/ConfigurationDialog.cpp
+++ b/src/gui/ConfigurationDialog.cpp
@@ -283,8 +283,8 @@ void ConfigurationDialog::createDialogContent()
 	idx = ui->deltaTAlgorithmComboBox->findData(core->getCurrentDeltaTAlgorithmKey(), Qt::UserRole, Qt::MatchCaseSensitive);
 	if (idx==-1)
 	{
-		// Use Espenak & Meeus (2006) as default
-		idx = ui->deltaTAlgorithmComboBox->findData(QVariant("EspenakMeeus"), Qt::UserRole, Qt::MatchCaseSensitive);
+		// Use Modified Espenak & Meeus (2006) as default
+		idx = ui->deltaTAlgorithmComboBox->findData(QVariant("EspenakMeeusModified"), Qt::UserRole, Qt::MatchCaseSensitive);
 	}
 	ui->deltaTAlgorithmComboBox->setCurrentIndex(idx);
 	connect(ui->deltaTAlgorithmComboBox, SIGNAL(currentIndexChanged(int)), this, SLOT(setDeltaTAlgorithm(int)));
@@ -1865,12 +1865,14 @@ void ConfigurationDialog::populateDeltaTAlgorithmsList()
 	algorithms->addItem(q_("Meeus & Simons (2000)"), "MeeusSimons");
 	algorithms->addItem(q_("Morrison & Stephenson (2004, 2005)"), "MorrisonStephenson2004");
 	algorithms->addItem(q_("Stephenson, Morrison & Hohenkerk (2016, 2021)"), "StephensonMorrisonHohenkerk2016");
-	// Espenak & Meeus (2006) used by default
-	algorithms->addItem(q_("Espenak & Meeus (2006)").append(" *"), "EspenakMeeus");
+	// Espenak & Meeus (2006)
+	algorithms->addItem(q_("Espenak & Meeus (2006)"), "EspenakMeeus");
 	// GZ: I want to try out some things. Something is still wrong with eclipses, see lp:1275092.
 	#ifndef NDEBUG
 	algorithms->addItem(q_("Espenak & Meeus (2006) no extra moon acceleration"), "EspenakMeeusZeroMoonAccel");
 	#endif
+	// Modified Espenak & Meeus (2006) used by default
+	algorithms->addItem(q_("Modified Espenak & Meeus (2006)").append(" *"), "EspenakMeeusModified");
 	algorithms->addItem(q_("Reijs (2006)"), "Reijs");
 	algorithms->addItem(q_("Banjevic (2006)"), "Banjevic");
 	algorithms->addItem(q_("Montenbruck & Pfleger (2000)"), "MontenbruckPfleger");