diff --git a/autotests/testseasons.h b/autotests/testseasons.h --- a/autotests/testseasons.h +++ b/autotests/testseasons.h @@ -28,8 +28,8 @@ { Q_OBJECT private Q_SLOTS: - void test2005(); - void test2007(); + void test_data(); + void test(); }; #endif diff --git a/autotests/testseasons.cpp b/autotests/testseasons.cpp --- a/autotests/testseasons.cpp +++ b/autotests/testseasons.cpp @@ -27,24 +27,43 @@ QTEST_MAIN(SeasonsTest) -void SeasonsTest::test2005() +Q_DECLARE_METATYPE(AstroSeasons::Season) + +void SeasonsTest::test_data() { - AstroSeasons as; - QVERIFY(as.seasonAtDate(QDate(2005, 3, 22)) == AstroSeasons::MarchEquinox); - QVERIFY(as.seasonAtDate(QDate(2005, 6, 22)) == AstroSeasons::JuneSolstice); - QVERIFY(as.seasonAtDate(QDate(2005, 9, 22)) == AstroSeasons::SeptemberEquinox); - QVERIFY(as.seasonAtDate(QDate(2005, 12, 22)) == AstroSeasons::DecemberSolstice); - QVERIFY(as.seasonAtDate(QDate(2005, 1, 22)) == AstroSeasons::None); - QVERIFY(as.seasonAtDate(QDate(2005, 12, 31)) == AstroSeasons::None); + QTest::addColumn("date"); + QTest::addColumn("season"); + + // Test data obtained from https://data.giss.nasa.gov/ar5/srvernal.html + QTest::newRow("2005-03-20") << QDate(2005, 3, 20) << AstroSeasons::MarchEquinox; + QTest::newRow("2005-06-21") << QDate(2005, 6, 21) << AstroSeasons::JuneSolstice; + QTest::newRow("2005-09-22") << QDate(2005, 9, 22) << AstroSeasons::SeptemberEquinox; + QTest::newRow("2005-12-21") << QDate(2005, 12, 21) << AstroSeasons::DecemberSolstice; + QTest::newRow("2005-01-22") << QDate(2005, 1, 22) << AstroSeasons::None; + QTest::newRow("2005-12-31") << QDate(2005, 12, 31) << AstroSeasons::None; + QTest::newRow("2007-03-21") << QDate(2007, 3, 21) << AstroSeasons::MarchEquinox; + QTest::newRow("2007-06-21") << QDate(2007, 6, 21) << AstroSeasons::JuneSolstice; + QTest::newRow("2007-09-23") << QDate(2007, 9, 23) << AstroSeasons::SeptemberEquinox; + QTest::newRow("2007-12-22") << QDate(2007, 12, 22) << AstroSeasons::DecemberSolstice; + QTest::newRow("2018-03-20") << QDate(2018, 3, 20) << AstroSeasons::MarchEquinox; + QTest::newRow("2018-06-21") << QDate(2018, 6, 21) << AstroSeasons::JuneSolstice; + QTest::newRow("2018-09-23") << QDate(2018, 9, 23) << AstroSeasons::SeptemberEquinox; + QTest::newRow("2018-12-21") << QDate(2018, 12, 21) << AstroSeasons::DecemberSolstice; + QTest::newRow("2020-03-20") << QDate(2020, 3, 20) << AstroSeasons::MarchEquinox; + QTest::newRow("2020-06-20") << QDate(2020, 6, 20) << AstroSeasons::JuneSolstice; + QTest::newRow("2020-09-22") << QDate(2020, 9, 22) << AstroSeasons::SeptemberEquinox; + QTest::newRow("2020-12-21") << QDate(2020, 12, 21) << AstroSeasons::DecemberSolstice; + QTest::newRow("2042-03-20") << QDate(2042, 3, 20) << AstroSeasons::MarchEquinox; + QTest::newRow("2042-06-21") << QDate(2042, 6, 21) << AstroSeasons::JuneSolstice; + QTest::newRow("2042-09-22") << QDate(2042, 9, 22) << AstroSeasons::SeptemberEquinox; + QTest::newRow("2042-12-21") << QDate(2042, 12, 21) << AstroSeasons::DecemberSolstice; } -void SeasonsTest::test2007() +void SeasonsTest::test() { + QFETCH(QDate, date); + QFETCH(AstroSeasons::Season, season); + AstroSeasons as; - QVERIFY(as.seasonAtDate(QDate(2007, 3, 22)) == AstroSeasons::MarchEquinox); - QVERIFY(as.seasonAtDate(QDate(2007, 6, 22)) == AstroSeasons::JuneSolstice); - QVERIFY(as.seasonAtDate(QDate(2007, 9, 22)) == AstroSeasons::SeptemberEquinox); - QVERIFY(as.seasonAtDate(QDate(2007, 12, 22)) == AstroSeasons::DecemberSolstice); - QVERIFY(as.seasonAtDate(QDate(2007, 1, 2)) == AstroSeasons::None); - QVERIFY(as.seasonAtDate(QDate(2007, 11, 30)) == AstroSeasons::None); + QCOMPARE(as.seasonAtDate(date), season); } diff --git a/src/astroseasons.h b/src/astroseasons.h --- a/src/astroseasons.h +++ b/src/astroseasons.h @@ -61,6 +61,16 @@ None }; + + /** + * Return the Gregorian date on which the season occurs in given year. + * + * @param season Season to return a date for + * @param year Year for which to return the date + * @since 5.50 + */ + static QDate seasonDate(Season season, int year); + /** Return the season for the specified Gregorian date. The enum 'None' is returned if one of the supported seasons diff --git a/src/astroseasons.cpp b/src/astroseasons.cpp --- a/src/astroseasons.cpp +++ b/src/astroseasons.cpp @@ -2,6 +2,7 @@ This file is part of the kholidays library. Copyright (c) 2004,2006-2007 Allen Winter + Copyright (c) 2018 Daniel Vrátil This library is free software; you can redistribute it and/or modify it under the terms of the GNU Library General Public @@ -21,11 +22,117 @@ #include "astroseasons.h" +#include +#include + #include #include using namespace KHolidays; +/* Correction tables and calculation algorithms below are based on the book + * Astronomical Algorithms by Jean Meeus, (c) 1991 by Willman-Bell, Inc., + * + * The correction tables are precise for years -1000 to +3000 but according + * to the book they can be used for several centuries before and after that + * period with the error still being quite small. As we only care about date + * and not the time, the precision of the algorithm is good enough even for + * greater time span, therefor no checks for the given year are performed + * anywhere in the code. + */ +namespace { + +double meanJDE(AstroSeasons::Season season, int year) +{ + if (year <= 1000) { + // Astronomical Algorithms, Jean Meeus, chapter 26, table 26.A + // mean season Julian dates for years -1000 to 1000 + const double y = year / 1000.0; + switch (season) { + case AstroSeasons::MarchEquinox: + return 1721139.29189 + (365242.13740 * y) + (0.06134 * std::pow(y, 2)) + (0.00111 * std::pow(y, 3)) - (0.00071 * std::pow(y, 4)); + case AstroSeasons::JuneSolstice: + return 1721233.25401 + (365241.72562 * y) - (0.05323 * std::pow(y, 2)) + (0.00907 * std::pow(y, 3)) + (0.00025 * std::pow(y, 4)); + case AstroSeasons::SeptemberEquinox: + return 1721325.70455 + (365242.49558 * y) - (0.11677 * std::pow(y, 2)) - (0.00297 * std::pow(y, 3)) + (0.00074 * std::pow(y, 4)); + case AstroSeasons::DecemberSolstice: + return 1721414.39987 + (365242.88257 * y) - (0.00769 * std::pow(y, 2)) - (0.00933 * std::pow(y, 3)) - (0.00006 * std::pow(y, 4)); + case AstroSeasons::None: + Q_ASSERT(false); + return 0; + } + } else { + // Astronomical Algorithms, Jean Meeus, chapter 26, table 26.B + // mean season Julian dates for years 1000 to 3000 + const double y = (year - 2000) / 1000.0; + switch (season) { + case AstroSeasons::MarchEquinox: + return 2451623.80984 + (365242.37404 * y) + (0.05169 * std::pow(y, 2)) - (0.00411 * std::pow(y, 3)) - (0.00057 * std::pow(y, 4)); + case AstroSeasons::JuneSolstice: + return 2451716.56767 + (365241.62603 * y) + (0.00325 * std::pow(y, 2)) + (0.00888 * std::pow(y, 3)) - (0.00030 * std::pow(y, 4)); + case AstroSeasons::SeptemberEquinox: + return 2451810.21715 + (365242.01767 * y) - (0.11575 * std::pow(y, 2)) + (0.00337 * std::pow(y, 3)) + (0.00078 * std::pow(y, 4)); + case AstroSeasons::DecemberSolstice: + return 2451900.05952 + (365242.74049 * y) - (0.06223 * std::pow(y, 2)) - (0.00823 * std::pow(y, 3)) + (0.00032 * std::pow(y, 4)); + case AstroSeasons::None: + Q_ASSERT(false); + return 0; + } + } + + return 0; +} + +double periodicTerms(double t) +{ + // Astronomical Algorithms, Jean Meeus, chapter 26, table 26.C + // The table gives the periodic terms in degrees, but the values are converted to radians + // at compile time so that they can be passed to std::cos() + struct Periodic { + constexpr Periodic(int a, double b_deg, double c_deg) + : a(a), b_rad(b_deg * (M_PI / 180.0)), c_rad(c_deg * (M_PI / 180.0)) + {} + + int a; + double b_rad; + double c_rad; + } periodic[] = { + {485, 324.96, 1934.136}, {203, 337.23, 32964.467}, {199, 342.08, 20.186}, {182, 27.85, 445267.112}, + {156, 73.14, 45036.886}, {136, 171.52, 22518.443}, { 77, 222.54, 65928.934}, { 74, 296.72, 3034.906}, + { 70, 243.58, 9037.513}, { 58, 119.81, 33718.147}, { 52, 297.17, 150.678}, { 50, 21.02, 2281.226}, + { 45, 247.54, 29929.562}, { 44, 325.15, 31555.956}, { 29, 60.93, 4443.417}, { 18, 155.12, 67555.328}, + { 17, 288.79, 4562.452}, { 16, 198.04, 62894.029}, { 14, 199.76, 31436.921}, { 12, 95.39, 14577.848}, + { 12, 287.11, 31931.756}, { 12, 320.81, 34777.259}, { 9, 227.73, 1222.114}, { 8, 15.45, 16859.074} + }; + + return std::accumulate(std::begin(periodic), std::end(periodic), 0.0, + [t](double s, const Periodic &p) { return s + p.a * std::cos(p.b_rad + p.c_rad * t); }); +} + +// Returns julian date of given season in given year +double seasonJD(AstroSeasons::Season season, int year) +{ + // Astronimical Algorithms, Jean Meeus, chapter 26 + const auto jde0 = meanJDE(season, year); + const auto T = (jde0 - 2451545.0) / 36525; + const auto W_deg = 35999.373 * T + 2.47; + const auto W_rad = W_deg * (M_PI / 180.0); + const auto dLambda = 1 + (0.0334 * std::cos(W_rad)) + (0.0007 * std::cos(2 * W_rad)); + const auto S = periodicTerms(T); + return jde0 + (0.00001 * S) / dLambda; +} + +} + +QDate AstroSeasons::seasonDate(Season season, int year) +{ + if (season == None) { + return {}; + } + const qint64 jd = round(seasonJD(season, year)); + return QDate::fromJulianDay(jd); +} + QString AstroSeasons::seasonNameAtDate(const QDate &date) { return seasonName(seasonAtDate(date)); @@ -50,23 +157,10 @@ AstroSeasons::Season AstroSeasons::seasonAtDate(const QDate &date) { - // see http://www.hermetic.ch/cal_sw/ve/ve.php - Season retSeason = None; - - const int year = date.year(); - //Use dumb method for now - if (date == QDate(year, 6, 22)) { - return JuneSolstice; - } - if (date == QDate(year, 12, 22)) { - return DecemberSolstice; + for (auto season : { JuneSolstice, DecemberSolstice, MarchEquinox, SeptemberEquinox }) { + if (seasonDate(season, date.year()) == date) { + return season; + } } - if (date == QDate(year, 3, 22)) { - return MarchEquinox; - } - if (date == QDate(year, 9, 22)) { - return SeptemberEquinox; - } - - return retSeason; + return None; }