File indexing completed on 2022-12-06 19:13:34

0001 /*
0002     This file is part of the kholidays library.
0003 
0004     SPDX-FileCopyrightText: 2004, 2006-2007 Allen Winter <winter@kde.org>
0005     SPDX-FileCopyrightText: 2018 Daniel Vrátil <dvratil@kde.org>
0006 
0007     SPDX-License-Identifier: LGPL-2.0-or-later
0008 */
0009 
0010 #include "astroseasons.h"
0011 
0012 #include <cmath>
0013 #include <numeric>
0014 
0015 #include <QCoreApplication>
0016 #include <QDate>
0017 
0018 using namespace KHolidays;
0019 
0020 /* Correction tables and calculation algorithms below are based on the book
0021  * Astronomical Algorithms by Jean Meeus, (c) 1991 by Willman-Bell, Inc.,
0022  *
0023  * The correction tables are precise for years -1000 to +3000 but according
0024  * to the book they can be used for several centuries before and after that
0025  * period with the error still being quite small. As we only care about date
0026  * and not the time, the precision of the algorithm is good enough even for
0027  * greater time span, therefore no checks for the given year are performed
0028  * anywhere in the code.
0029  */
0030 namespace
0031 {
0032 double meanJDE(AstroSeasons::Season season, int year)
0033 {
0034     if (year <= 1000) {
0035         // Astronomical Algorithms, Jean Meeus, chapter 26, table 26.A
0036         // mean season Julian dates for years -1000 to 1000
0037         const double y = year / 1000.0;
0038         switch (season) {
0039         case AstroSeasons::MarchEquinox:
0040             return 1721139.29189 + (365242.13740 * y) + (0.06134 * std::pow(y, 2)) + (0.00111 * std::pow(y, 3)) - (0.00071 * std::pow(y, 4));
0041         case AstroSeasons::JuneSolstice:
0042             return 1721233.25401 + (365241.72562 * y) - (0.05323 * std::pow(y, 2)) + (0.00907 * std::pow(y, 3)) + (0.00025 * std::pow(y, 4));
0043         case AstroSeasons::SeptemberEquinox:
0044             return 1721325.70455 + (365242.49558 * y) - (0.11677 * std::pow(y, 2)) - (0.00297 * std::pow(y, 3)) + (0.00074 * std::pow(y, 4));
0045         case AstroSeasons::DecemberSolstice:
0046             return 1721414.39987 + (365242.88257 * y) - (0.00769 * std::pow(y, 2)) - (0.00933 * std::pow(y, 3)) - (0.00006 * std::pow(y, 4));
0047         case AstroSeasons::None:
0048             Q_ASSERT(false);
0049             return 0;
0050         }
0051     } else {
0052         // Astronomical Algorithms, Jean Meeus, chapter 26, table 26.B
0053         // mean season Julian dates for years 1000 to 3000
0054         const double y = (year - 2000) / 1000.0;
0055         switch (season) {
0056         case AstroSeasons::MarchEquinox:
0057             return 2451623.80984 + (365242.37404 * y) + (0.05169 * std::pow(y, 2)) - (0.00411 * std::pow(y, 3)) - (0.00057 * std::pow(y, 4));
0058         case AstroSeasons::JuneSolstice:
0059             return 2451716.56767 + (365241.62603 * y) + (0.00325 * std::pow(y, 2)) + (0.00888 * std::pow(y, 3)) - (0.00030 * std::pow(y, 4));
0060         case AstroSeasons::SeptemberEquinox:
0061             return 2451810.21715 + (365242.01767 * y) - (0.11575 * std::pow(y, 2)) + (0.00337 * std::pow(y, 3)) + (0.00078 * std::pow(y, 4));
0062         case AstroSeasons::DecemberSolstice:
0063             return 2451900.05952 + (365242.74049 * y) - (0.06223 * std::pow(y, 2)) - (0.00823 * std::pow(y, 3)) + (0.00032 * std::pow(y, 4));
0064         case AstroSeasons::None:
0065             Q_ASSERT(false);
0066             return 0;
0067         }
0068     }
0069 
0070     return 0;
0071 }
0072 
0073 double periodicTerms(double t)
0074 {
0075     // Astronomical Algorithms, Jean Meeus, chapter 26, table 26.C
0076     // The table gives the periodic terms in degrees, but the values are converted to radians
0077     // at compile time so that they can be passed to std::cos()
0078     struct Periodic {
0079         constexpr Periodic(int a, double b_deg, double c_deg)
0080             : a(a)
0081             , b_rad(b_deg * (M_PI / 180.0))
0082             , c_rad(c_deg * (M_PI / 180.0))
0083         {
0084         }
0085 
0086         int a;
0087         double b_rad;
0088         double c_rad;
0089     }
0090     // clang-format off
0091     periodic[] = {
0092         {485, 324.96,  1934.136}, {203, 337.23, 32964.467}, {199, 342.08,    20.186}, {182,  27.85, 445267.112},
0093         {156,  73.14, 45036.886}, {136, 171.52, 22518.443}, { 77, 222.54, 65928.934}, { 74, 296.72,   3034.906},
0094         { 70, 243.58,  9037.513}, { 58, 119.81, 33718.147}, { 52, 297.17,   150.678}, { 50,  21.02,   2281.226},
0095         { 45, 247.54, 29929.562}, { 44, 325.15, 31555.956}, { 29,  60.93,  4443.417}, { 18, 155.12,  67555.328},
0096         { 17, 288.79,  4562.452}, { 16, 198.04, 62894.029}, { 14, 199.76, 31436.921}, { 12,  95.39,  14577.848},
0097         { 12, 287.11, 31931.756}, { 12, 320.81, 34777.259}, {  9, 227.73,  1222.114}, {  8,  15.45,  16859.074}
0098     };
0099     // clang-format on
0100 
0101     return std::accumulate(std::begin(periodic), std::end(periodic), 0.0, [t](double s, const Periodic &p) {
0102         return s + p.a * std::cos(p.b_rad + p.c_rad * t);
0103     });
0104 }
0105 
0106 // Returns julian date of given season in given year
0107 double seasonJD(AstroSeasons::Season season, int year)
0108 {
0109     // Astronimical Algorithms, Jean Meeus, chapter 26
0110     const auto jde0 = meanJDE(season, year);
0111     const auto T = (jde0 - 2451545.0) / 36525;
0112     const auto W_deg = 35999.373 * T + 2.47;
0113     const auto W_rad = W_deg * (M_PI / 180.0);
0114     const auto dLambda = 1 + (0.0334 * std::cos(W_rad)) + (0.0007 * std::cos(2 * W_rad));
0115     const auto S = periodicTerms(T);
0116     return jde0 + (0.00001 * S) / dLambda;
0117 }
0118 
0119 }
0120 
0121 QDate AstroSeasons::seasonDate(Season season, int year)
0122 {
0123     if (season == None) {
0124         return {};
0125     }
0126     const qint64 jd = round(seasonJD(season, year));
0127     return QDate::fromJulianDay(jd);
0128 }
0129 
0130 QString AstroSeasons::seasonNameAtDate(const QDate &date)
0131 {
0132     return seasonName(seasonAtDate(date));
0133 }
0134 
0135 QString AstroSeasons::seasonName(AstroSeasons::Season season)
0136 {
0137     switch (season) {
0138     case JuneSolstice:
0139         return QCoreApplication::translate("AstroSeasons", "June Solstice");
0140     case DecemberSolstice:
0141         return QCoreApplication::translate("AstroSeasons", "December Solstice");
0142     case MarchEquinox:
0143         return QCoreApplication::translate("AstroSeasons", "March Equinox");
0144     case SeptemberEquinox:
0145         return QCoreApplication::translate("AstroSeasons", "September Equinox");
0146     case None:
0147         return QString();
0148     }
0149     return QString();
0150 }
0151 
0152 AstroSeasons::Season AstroSeasons::seasonAtDate(const QDate &date)
0153 {
0154     for (auto season : {JuneSolstice, DecemberSolstice, MarchEquinox, SeptemberEquinox}) {
0155         if (seasonDate(season, date.year()) == date) {
0156             return season;
0157         }
0158     }
0159     return None;
0160 }