StelUtils.cpp   StelUtils.cpp
skipping to change at line 1115 skipping to change at line 1115
// For the standard epochs for many formulae, we use // For the standard epochs for many formulae, we use
// J2000.0=2000-jan-1.5=2451545.0, // J2000.0=2000-jan-1.5=2451545.0,
// 1900.0=1900-jan-0.5=2415020.0 // 1900.0=1900-jan-0.5=2415020.0
// 1820.0=1820-jan-0.5=2385800.0 // 1820.0=1820-jan-0.5=2385800.0
// 1810.0=1810-jan-0.5=2382148.0 // 1810.0=1810-jan-0.5=2382148.0
// 1800.0=1800-jan-0.5=2378496.0 // 1800.0=1800-jan-0.5=2378496.0
// 1735.0=1735-jan-0.5=2354755.0 // 1735.0=1735-jan-0.5=2354755.0
// 1625.0=1625-jan-0.5=2314579.0 // 1625.0=1625-jan-0.5=2314579.0
*/ */
double double // Implementation of algorithm by Espenak & Meeus (2006) for DeltaT computa
tion
double getDeltaTByEspenakMeeus(const double jDay)
{ {
int year, month, day;
getDateFromJulianDay(jDay, &year, &month, &day);
// Note: the method here is adapted from // Note: the method here is adapted from
// "Five Millennium Canon of Solar Eclipses" [Espenak and Meeus, 200 6] // "Five Millennium Canon of Solar Eclipses" [Espenak and Meeus, 200 6]
// A summary is described here: // A summary is described here:
// http://eclipse.gsfc.nasa.gov/SEhelp/deltatpoly2004.html // http://eclipse.gsfc.nasa.gov/SEhelp/deltatpoly2004.html
// GZ: I replaced the std::pow() calls by Horner's scheme with rever sed factors, it's more accurate and efficient. // GZ: I replaced the std::pow() calls by Horner's scheme with rever sed factors, it's more accurate and efficient.
// Old code left for readability, but can also be deleted. // Old code left for readability, but can also be deleted.
double y = year+((month-1)*30.5+day/31*30.5)/366;
// set the default value for Delta T // set the default value for Delta T
double u = (y-1820)/100.; double u = (y-1820)/100.;
double r = (-20 + 32 * u * u); double r = (-20 + 32 * u * u);
if (y < -500) if (y < -500)
{ {
// values are equal to defaults! // values are equal to defaults!
// u = (y-1820)/100.; // u = (y-1820)/100.;
// r = (-20 + 32 * u * u); // r = (-20 + 32 * u * u);
} }
skipping to change at line 1220 skipping to change at line 1226
else if (y < 2150) else if (y < 2150)
{ {
//r = (-20 + 32 * std::pow((y-1820)/100,2) - 0.5628 * (2150 - y)); //r = (-20 + 32 * std::pow((y-1820)/100,2) - 0.5628 * (2150 - y));
// r has been precomputed before, just add the term patching the discontinuity // r has been precomputed before, just add the term patching the discontinuity
r -= 0.5628*(2150.0-y); r -= 0.5628*(2150.0-y);
} }
return r; return r;
} }
// Implementation of algorithm by Schoch (1931) for DeltaT computation // Implementation of algorithm by Schoch (1931) for DeltaT computation
double getDeltaTBySchoch(const double jDay) double getDeltaTBySchoch(const double jDay)
{ {
double u=(jDay-2378496.0)/36525.0; // (1800-jan-0.5) double u=(jDay-2378496.0)/36525.0; // (1800-jan-0.5)
return -36.28 + 36.28*std::pow(u,2); return -36.28 + 36.28*std::pow(u,2);
} }
// Implementation of algorithm by Clemence (1948) for DeltaT computation // Implementation of algorithm by Clemence (1948) for DeltaT computation
double getDeltaTByClemence(const double jDay) double getDeltaTByClemence(const double jDay)
{ {
double u=(jDay-2415020.0)/36525.0; // (1900-jan-0.5) double u=(jDay-2415020.0)/36525.0; // (1900-jan-0.5)
return +8.72 + 26.75*u + 11.22*std::pow(u,2); return +8.72 + 26.75*u + 11.22*std::pow(u,2);
} }
// Implementation of algorithm by IAU (1952) for DeltaT computation // Implementation of algorithm by IAU (1952) for DeltaT computation
double getDeltaTByIAU(const double jDay) double getDeltaTByIAU(const double jDay)
{ {
double u=(jDay-2415020.0)/36525.0; // (1900-jan-0.5) double u=(jDay-2415020.0)/36525.0; // (1900-jan-0.5)
// TODO: Calculate Moon's longitude fluctuation // TODO: Calculate Moon's longitude fluctuation
return (29.9*/ ; return (29.950*u +72.318)*u +24.349 /* + 1.82144*b */ ;
} }
// Implementation of algorithm by Astronomical Ephemeris (1960) for DeltaT computation, also used by Mucke&Meeus, Canon of Solar Eclipses, Vienna 1983 // Implementation of algorithm by Astronomical Ephemeris (1960) for DeltaT computation, also used by Mucke&Meeus, Canon of Solar Eclipses, Vienna 1983
double getDeltaTByAstronomicalEphemeris(const double jDay) double getDeltaTByAstronomicalEphemeris(const double jDay)
{ {
double u=(jDay-2415020.0)/36525.0; // (1900-jan-0.5) double u=(jDay-2415020.0)/36525.0; // (1900-jan-0.5)
// TODO: Calculate Moon's longitude fluctuation // TODO: Calculate Moon's longitude fluctuation
// Note: also Mucke&Meeus 1983 ignore b // Note: also Mucke&Meeus 1983 ignore b
return (29.9*/ ; return (29.949*u +72.3165)*u +24.349 /* + 1.821*b*/ ;
} }
// Implementation of algorithm by Tuckerman (1962, 1964) & Goldstine (1973) for DeltaT computation // Implementation of algorithm by Tuckerman (1962, 1964) & Goldstine (1973) for DeltaT computation
double getDeltaTByTuckermanGoldstine(const double jDay) double getDeltaTByTuckermanGoldstine(const double jDay)
{ {
double u=(jDay-2415020.0)/36525.0; // (1900-jan-0.5) double u=(jDay-2415020.0)/36525.0; // (1900-jan-0.5)
return (36.79*u +35.06)*u + 4.87; return (36.79*u +35.06)*u + 4.87;
} }
// Implementation of algorithm by Muller & Stephenson (1975) for DeltaT com putation // Implementation of algorithm by Muller & Stephenson (1975) for DeltaT com putation
skipping to change at line 1316 skipping to change at line 1311
// Implementation of algorithm by Stephenson & Morrison (1984) for DeltaT c omputation // Implementation of algorithm by Stephenson & Morrison (1984) for DeltaT c omputation
double getDeltaTByStephensonMorrison1984(const double jDay) double getDeltaTByStephensonMorrison1984(const double jDay)
{ {
int year, month, day; int year, month, day;
double deltaT = 0.; double deltaT = 0.;
getDateFromJulianDay(jDay, &year, &month, &day); getDateFromJulianDay(jDay, &year, &month, &day);
double yeardec=year+((month-1)*30.5+day/31*30.5)/366; double yeardec=year+((month-1)*30.5+day/31*30.5)/366;
double u = (yeardec-1800)/100; double u = (yeardec-1800)/100;
if (-391 < year year <= 948) if (-391 < year && year <= 948)
deltaT = (44.3*u +320.0)*u +1360.0; deltaT = (44.3*u +320.0)*u +1360.0;
if (948 < year year <= 1600) if (948 < year && year <= 1600)
deltaT = 25.5*u*u; deltaT = 25.5*u*u;
return deltaT; return deltaT;
} }
// Implementation of algorithm by Stephenson & Morrison (1995) for DeltaT c omputation // Implementation of algorithm by Stephenson & Morrison (1995) for DeltaT c omputation
double getDeltaTByStephensonMorrison1995(const double jDay) double getDeltaTByStephensonMorrison1995(const double jDay)
{ {
double u=(jDay-2385800.0)/36525.0; // (1820-jan-0.5) double u=(jDay-2385800.0)/36525.0; // (1820-jan-0.5)
return -20.0 + 31.0*u*u; return -20.0 + 31.0*u*u;
skipping to change at line 1344 skipping to change at line 1339
int year, month, day; int year, month, day;
double u; double u;
double deltaT = 0.; double deltaT = 0.;
getDateFromJulianDay(jDay, &year, &month, &day); getDateFromJulianDay(jDay, &year, &month, &day);
double yeardec=year+((month-1)*30.5+day/31*30.5)/366; double yeardec=year+((month-1)*30.5+day/31*30.5)/366;
if (year <= 948) if (year <= 948)
{ {
u = (yeardec-948)/100; u = (yeardec-948)/100;
deltaT = (46.5*u -405.0)*u + 1830.0; deltaT = (46.5*u -405.0)*u + 1830.0;
} }
if (948 < year year <= 1600) if (948 < year && year <= 1600)
{ {
u = (yeardec-1850)/100; u = (yeardec-1850)/100;
deltaT = 2.5*u*u; deltaT = 22.5*u*u;
} }
return deltaT; return deltaT;
} }
// Implementation of algorithm by Espenak (1987, 1989) for DeltaT computati on // Implementation of algorithm by Espenak (1987, 1989) for DeltaT computati on
double getDeltaTByEspenak(const double jDay) double getDeltaTByEspenak(const double jDay)
{ {
double u=(jDay-2451545.0)/36525.0; // (2000-jan-1.5) double u=(jDay-2451545.0)/36525.0; // (2000-jan-1.5)
return (64.3*u +61.0)*u +67.0; return (64.3*u +61.0)*u +67.0;
skipping to change at line 1387 skipping to change at line 1382
// Implementation of algorithm by Chapront-TouzĂ© & Chapront (1991) for Delt aT computation // Implementation of algorithm by Chapront-TouzĂ© & Chapront (1991) for Delt aT computation
double getDeltaTByChaprontTouze(const double jDay) double getDeltaTByChaprontTouze(const double jDay)
{ {
int year, month, day; int year, month, day;
double deltaT = 0.; double deltaT = 0.;
getDateFromJulianDay(jDay, &year, &month, &day); getDateFromJulianDay(jDay, &year, &month, &day);
double u=(jDay-2451545.0)/36525.0; // (2000-jan-1.5) double u=(jDay-2451545.0)/36525.0; // (2000-jan-1.5)
if (-391 < year year <= 948) if (-391 < year && year <= 948)
deltaT = (42.4*u + 2177.0; deltaT = (42.4*u +495.0)*u + 2177.0;
if (948 < year year <= 1600) if (948 < year && year <= 1600)
deltaT = (23.6*u +100.0)*u + 102.0; deltaT = (23.6*u +100.0)*u + 102.0;
return deltaT; return deltaT;
} }
// Implementation of algorithm by JPL Horizons for DeltaT computation // Implementation of algorithm by JPL Horizons for DeltaT computation
double getDeltaTByJPLHorizons(const double jDay) double getDeltaTByJPLHorizons(const double jDay)
{ // GZ: TODO: FIXME! It does not make sense to have zeros after 1620 in a JPL Horizons compatible implementation! { // GZ: TODO: FIXME! It does not make sense to have zeros after 1620 in a JPL Horizons compatible implementation!
int year, month, day; int year, month, day;
double u; double u;
double deltaT = 0.; double deltaT = 0.;
getDateFromJulianDay(jDay, &year, &month, &day); getDateFromJulianDay(jDay, &year, &month, &day);
if (-2999 < year year < 948) if (-2999 < year && year < 948)
{ {
u=(jDay-2385800.0)/36525.0; // (1820-jan-1.5) u=(jDay-2385800.0)/36525.0; // (1820-jan-1.5)
deltaT = 31.0*u*u; deltaT = 31.0*u*u;
} }
if (948 < year year <= 1620) if (948 < year && year <= 1620)
{ {
u=(jDay-2451545.0)/36525.0; // (2000-jan-1.5) u=(jDay-2451545.0)/36525.0; // (2000-jan-1.5)
deltaT = (22.5*u +67.5)*u + 50.6; deltaT = (22.5*u +67.5)*u + 50.6;
} }
return deltaT; return deltaT;
} }
// Implementation of algorithm by Morrison & Stephenson (2004, 2005) for De ltaT computation // Implementation of algorithm by Morrison & Stephenson (2004, 2005) for De ltaT computation
double getDeltaTByMorrisonStephenson2004(const double jDay) double getDeltaTByMorrisonStephenson2004(const double jDay)
skipping to change at line 1544 skipping to change at line 1539
else if (year < 1770) else if (year < 1770)
{ {
u = 2.70 + ub; u = 2.70 + ub;
//deltaT = +10.2 + 11.3*u - std::pow(u,2) - 16.0*std::pow(u, 3) + 70.0*std::pow(u,4); //deltaT = +10.2 + 11.3*u - std::pow(u,2) - 16.0*std::pow(u, 3) + 70.0*std::pow(u,4);
deltaT = (((70.0*u -16.0)*u -1.0)*u +11.3)*u +10.2; deltaT = (((70.0*u -16.0)*u -1.0)*u +11.3)*u +10.2;
} }
else if (year < 1820) else if (year < 1820)
{ {
u = 2.05 + ub; u = 2.05 + ub;
//deltaT = +14.7 - 18.8*u - 22.0*std::pow(u,2) + 173.0*std:: pow(u,3) + 6.0*std::pow(u,4); //deltaT = +14.7 - 18.8*u - 22.0*std::pow(u,2) + 173.0*std:: pow(u,3) + 6.0*std::pow(u,4);
deltaT = (((6.0*u 173.0)*u -22.0)*u -18.8)*u +14.7; deltaT = (((6.0*u +173.0)*u -22.0)*u -18.8)*u +14.7;
} }
else if (year < 1870) else if (year < 1870)
{ {
u = 1.55 + ub; u = 1.55 + ub;
//deltaT = +5.7 + 12.7*u + 111.0*std::pow(u,2) - 534.0*std:: //deltaT = +5.7 + 12.7*u + 111.0*std::pow(u,2) - 534.0*std::
pow(u,3) 1654.0*std::pow(u,4); pow(u,3) - 1654.0*std::pow(u,4);
deltaT = -534.0)*u +111)*u +12.7)*u +5.7; deltaT = (((-1654.0*u -534.0)*u +111)*u +12.7)*u +5.7;
} }
else if (year < 1900) else if (year < 1900)
{ {
u = 1.15 + ub; u = 1.15 + ub;
//deltaT = -5.8 - 14.6*u + 27.0*std::pow(u,2) + 101.0*std::p ow(u,3) + 8234.0*std::pow(u,4); //deltaT = -5.8 - 14.6*u + 27.0*std::pow(u,2) + 101.0*std::p ow(u,3) + 8234.0*std::pow(u,4);
deltaT = (((8234.0*u +101.0)*u +27.0)*u - 14.6)*u -5.8; deltaT = (((8234.0*u +101.0)*u +27.0)*u - 14.6)*u -5.8;
} }
else if (year < 1940) else if (year < 1940)
{ {
u = 0.80 + ub; u = 0.80 + ub;
//deltaT = +21.4 + 67.0*u 443.0*std::pow(u,2) + 19.0*std:: //deltaT = +21.4 + 67.0*u - 443.0*std::pow(u,2) + 19.0*std::
pow(u,3) + 4441.0*std::pow(u,4); pow(u,3) + 4441.0*std::pow(u,4);
deltaT = (((4441.0*u + 19.0)*u +67.0)*u +21.4; deltaT = (((4441.0*u + 19.0)*u -443.0)*u +67.0)*u +21.4;
} }
else if (year < 1990) else if (year < 1990)
{ {
u = 0.35 + ub; u = 0.35 + ub;
//deltaT = +36.2 + 74.0*u + 189.0*std::pow(u,2) - 140.0*std: :pow(u,3) - 1883.0*std::pow(u,4); //deltaT = +36.2 + 74.0*u + 189.0*std::pow(u,2) - 140.0*std: :pow(u,3) - 1883.0*std::pow(u,4);
deltaT = (((-1883.0*u -140.0)*u +189.0)*u +74.0)*u +36.2; deltaT = (((-1883.0*u -140.0)*u +189.0)*u +74.0)*u +36.2;
} }
else if (year <= 2000) else if (year <= 2000)
{ {
u = 0.05 + ub; u = 0.05 + ub;
//deltaT = +60.8 + 82.0*u 188.0*std::pow(u,2) - 5034.0*std //deltaT = +60.8 + 82.0*u - 188.0*std::pow(u,2) - 5034.0*std
::pow(u,3); ::pow(u,3);
deltaT = ((-5034.0*u +82.0)*u +60.8; deltaT = ((-5034.0*u -188.0)*u +82.0)*u +60.8;
} }
return deltaT; return deltaT;
} }
// Implementation of algorithm by Reingold & Dershowitz (Cal. Calc. 1997, 2 001, 2007, Cal. Tab. 2002) for DeltaT computation. // Implementation of algorithm by Reingold & Dershowitz (Cal. Calc. 1997, 2 001, 2007, Cal. Tab. 2002) for DeltaT computation.
// GZ: Created as yet another multi-segment polynomial fit through the tabl e in Meeus: Astronomical Algorithms (1991). // GZ: Created as yet another multi-segment polynomial fit through the tabl e in Meeus: Astronomical Algorithms (1991).
// GZ: Note that only the Third edition (2007) adds the 1700-1799 term. // GZ: Note that only the Third edition (2007) adds the 1700-1799 term.
// GZ: More efficient reimplementation with stricter adherence to the sourc e. // GZ: More efficient reimplementation with stricter adherence to the sourc e.
double getDeltaTByReingoldDershowitz(const double jDay) double getDeltaTByReingoldDershowitz(const double jDay)
skipping to change at line 1621 skipping to change at line 1616
deltaT = (((((((-0.212591*c +0.677066)*c -0.861938)* c +0.553040)*c -0.181133)*c +0.025184)*c +0.000297)*c -0.00002) * 86400.0; deltaT = (((((((-0.212591*c +0.677066)*c -0.861938)* c +0.553040)*c -0.181133)*c +0.025184)*c +0.000297)*c -0.00002) * 86400.0;
} }
else //if (year >= 1800) else //if (year >= 1800)
{ {
deltaT = ((((((((((2.043794*c +11.636204)*c +28.3162 89)*c +38.291999)*c +31.332267)*c +15.845535)*c +4.867575)*c +0.865736)*c + 0.083563)*c +0.003844)*c -0.000009) * 86400.0; deltaT = ((((((((((2.043794*c +11.636204)*c +28.3162 89)*c +38.291999)*c +31.332267)*c +15.845535)*c +4.867575)*c +0.865736)*c + 0.083563)*c +0.003844)*c -0.000009) * 86400.0;
} }
} }
else if (year >= 1700) else if (year >= 1700)
{ // This term was added in the third edition (2007), its omission w as a fault of the authors! { // This term was added in the third edition (2007), its omission w as a fault of the authors!
double yDiff1700 = year-1700.0; double yDiff1700 = year-1700.0;
deltaT = ((-0.0000266484*yDiff1700 +0.003336121)*yDiff1700 0.005092142)*yDiff1700 + 8.118780842; deltaT = ((-0.0000266484*yDiff1700 +0.003336121)*yDiff1700 - 0.005092142)*yDiff1700 + 8.118780842;
} }
else if (year >= 1620) else if (year >= 1620)
{ {
double yDiff1600 = year-1600.0; double yDiff1600 = year-1600.0;
deltaT = (0.0219167*yDiff1600 -4.0675)*yDiff1600 +196.58333; deltaT = (0.0219167*yDiff1600 -4.0675)*yDiff1600 +196.58333;
} }
return deltaT; return deltaT;
} }
// Implementation of algorithm by Banjevic (2006) for DeltaT computation.
double getDeltaTByBanjevic(const double jDay)
{
int year, month, day;
getDateFromJulianDay(jDay, &year, &month, &day);
double u, c;
if (year>=-2020 && year<=-700)
{
u = (jDay-2378496.0)/36525.0; // 1800.0=1800-jan-0.5=2378496
.0
c = 30.86;
}
else if (year>-700 && year<=1620)
{
u = (jDay-2385800.0)/36525.0; // 1820.0=1820-jan-0.5=2385800
.0
c = 31;
}
else
{
u = 0.;
c = 0.;
}
return c*u*u;
}
// Implementation of algorithm by Islam, Sadiq & Qureshi (2008 + revisited
2013) for DeltaT computation.
{
int year, month, day;
getDateFromJulianDay(jDay, &year, &month, &day);
double deltaT = 0.0; // Return deltaT = 0 outside valid range.
double u;
const double ub=(jDay-2454101.0)/36525.0; // (2007-jan-0.5)
if (year >= 1620 && year <= 1698)
{
u = 3.48 + ub;
deltaT = (((1162.805 * u - 273.116) * u + 14.523) * u - 105.
262) * u + 38.067;
}
else if (year <= 1806)
{
u = 2.545 + ub;
deltaT = (((-71.724 * u - 39.048) * u + 7.591) * u + 13.893)
* u + 13.759;
}
else if (year <= 1872)
{
u = 1.675 + ub;
deltaT = (((-1612.55 * u - 157.977) * u + 161.524) * u - 3.6
54) * u + 5.859;
}
else if (year <= 1906)
{
u = 1.175 + ub;
deltaT = (((6250.501 * u + 1006.463) * u + 139.921) * u - 2.
732) * u - 6.203;
}
else if (year <= 1953)
{
// revised 2013 per email
u = 0.77 + ub;
deltaT = (((-390.785 * u + 901.514) * u - 88.044) * u + 8.99
7) * u + 24.006;
}
else if (year <= 2007)
{
// revised 2013 per email
u = 0.265 + ub;
deltaT = (((1314.759 * u - 296.018) * u - 101.898) * u + 88.
659) * u + 49.997;
}
return deltaT;
}
double getMoonSecularAcceleration(const double jDay, const double nd) double getMoonSecularAcceleration(const double jDay, const double nd)
{ {
int year, month, day; int year, month, day;
getDateFromJulianDay(jDay, &year, &month, &day); getDateFromJulianDay(jDay, &year, &month, &day);
double yeardec=year+((month-1)*30.5+day/31*30.5)/366.0; double yeardec=year+((month-1)*30.5+day/31*30.5)/366.0;
double t = (yeardec-1955.5)/100.0; double t = (yeardec-1955.5)/100.0;
// n.dot for secular acceleration of the Moon in ELP2000-82B // n.dot for secular acceleration of the Moon in ELP2000-82B
// have value -23.8946 "/cy/cy // have value -23.8946 "/cy/cy
return -0.91072 * (-23.8946 + std::abs(nd))*t*t; return -0.91072 * (-23.8946 + std::abs(nd))*t*t;
} }
double getDeltaTStandardError(const double jDay) double getDeltaTStandardError(const double jDay)
{ {
int year, month, day; int year, month, day;
getDateFromJulianDay(jDay, &year, &month, &day); getDateFromJulianDay(jDay, &year, &month, &day);
//double yeardec=year+((month-1)*30.5+day/31*30.5)/366; //double yeardec=year+((month-1)*30.5+day/31*30.5)/366;
double sigma = -1.; double sigma = -1.;
if (-1000 <= year year <= 1600) if (-1000 <= year && year <= 1600)
{ {
double cDiff1820= (jDay-2385800.0)/36525.0; // 1820.0=182 0-jan-0.5=2385800.0 double cDiff1820= (jDay-2385800.0)/36525.0; // 1820.0=182 0-jan-0.5=2385800.0
// sigma = std::pow((yeardec-1820.0)/100,2); // sigma(DeltaT ) = 0.8*u^2 // sigma = std::pow((yeardec-1820.0)/100,2); // sigma(DeltaT ) = 0.8*u^2
sigma = 0.8 * cDiff1820 * cDiff1820; sigma = 0.8 * cDiff1820 * cDiff1820;
} }
return sigma; return sigma;
} }
} // end of the StelUtils namespace } // end of the StelUtils namespace
End of changes. 21 change blocks.
37 lines changed or deleted 112 lines changed or added

This html diff was produced by rfcdiff 1.41. The latest version is available from http://tools.ietf.org/tools/rfcdiff/