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

0001 /*
0002     This file is part of the kholidays library.
0003 
0004     SPDX-FileCopyrightText: 2004, 2007, 2009 Allen Winter <winter@kde.org>
0005 
0006     SPDX-License-Identifier: LGPL-2.0-or-later
0007 */
0008 
0009 #include "lunarphase.h"
0010 
0011 #include <QCoreApplication>
0012 #include <QDateTime>
0013 
0014 using namespace KHolidays;
0015 
0016 QString LunarPhase::phaseNameAtDate(const QDate &date)
0017 {
0018     return phaseName(phaseAtDate(date));
0019 }
0020 
0021 QString LunarPhase::phaseName(LunarPhase::Phase phase)
0022 {
0023     switch (phase) {
0024     case NewMoon:
0025         return (QCoreApplication::translate("LunarPhase", "New Moon"));
0026     case FullMoon:
0027         return (QCoreApplication::translate("LunarPhase", "Full Moon"));
0028     case FirstQuarter:
0029         return (QCoreApplication::translate("LunarPhase", "First Quarter Moon"));
0030     case LastQuarter:
0031         return (QCoreApplication::translate("LunarPhase", "Last Quarter Moon"));
0032     case None:
0033         return QString();
0034     case WaxingCrescent:
0035         return (QCoreApplication::translate("LunarPhase", "Waxing Crescent"));
0036     case WaxingGibbous:
0037         return (QCoreApplication::translate("LunarPhase", "Waxing Gibbous"));
0038     case WaningGibbous:
0039         return (QCoreApplication::translate("LunarPhase", "Waning Gibbous"));
0040     case WaningCrescent:
0041         return (QCoreApplication::translate("LunarPhase", "Waning Crescent"));
0042     }
0043     return QString();
0044 }
0045 
0046 static double phaseAngle(int64_t unixmsec);
0047 
0048 LunarPhase::Phase LunarPhase::phaseAtDate(const QDate &date)
0049 {
0050     Phase retPhase = None;
0051 
0052     const QTime midnight(0, 0, 0);
0053     const QDateTime todayStart(date, midnight, Qt::UTC);
0054     const double startAngle = phaseAngle(todayStart.toMSecsSinceEpoch());
0055     const QDateTime todayEnd(date.addDays(1), midnight, Qt::UTC);
0056     const double endAngle = phaseAngle(todayEnd.toMSecsSinceEpoch());
0057 
0058     if (startAngle > endAngle) {
0059         retPhase = NewMoon;
0060     } else if (startAngle < 90.0 && endAngle > 90.0) {
0061         retPhase = FirstQuarter;
0062     } else if (startAngle < 180.0 && endAngle > 180.0) {
0063         retPhase = FullMoon;
0064     } else if (startAngle < 270.0 && endAngle > 270.0) {
0065         retPhase = LastQuarter;
0066     } else if (endAngle < 90.0) {
0067         retPhase = WaxingCrescent;
0068     } else if (endAngle < 180.0) {
0069         retPhase = WaxingGibbous;
0070     } else if (endAngle < 270.0) {
0071         retPhase = WaningGibbous;
0072     } else if (endAngle < 360.0) {
0073         retPhase = WaningCrescent;
0074     }
0075 
0076     return retPhase;
0077 }
0078 
0079 /*
0080     Algorithm adapted from Moontool by John Walker.
0081     Moontool is in the public domain: "Do what thou wilt shall be the whole of the law".
0082     Unnecessary calculations have been removed.
0083     See https://fourmilab.ch/moontool .
0084 */
0085 
0086 #include <cmath>
0087 
0088 constexpr int64_t epoch = 315446400000; // "1980 January 0.0", a.k.a. midnight on 1979-12-31, in ms Unix time
0089 constexpr double elonge = 278.833540; // ecliptic longitude of sun at epoch
0090 constexpr double elongp = 282.596403; // ecliptic longitude of sun at perigee
0091 constexpr double earthEcc = 0.016718; // eccentricity of earth orbit
0092 static const double ecPrefactor = sqrt((1 + earthEcc) / (1 - earthEcc));
0093 
0094 constexpr double mmlong = 64.975464; // mean longitude of moon at epoch
0095 constexpr double mmlongp = 349.383063; // mean longitude of moon at perigee
0096 
0097 constexpr double PI = 3.14159265358979323846;
0098 
0099 static double fixAngle(double degrees)
0100 {
0101     return degrees - floor(degrees / 360.0) * 360.0;
0102 }
0103 
0104 static constexpr double radToDeg(double rad)
0105 {
0106     return rad / PI * 180.0;
0107 }
0108 
0109 static constexpr double degToRad(double deg)
0110 {
0111     return deg / 180.0 * PI;
0112 }
0113 
0114 constexpr double epsilon = 1e-6;
0115 
0116 static double kepler(double m, double ecc)
0117 {
0118     double mrad = degToRad(m);
0119     double e = mrad;
0120     double delta;
0121     do {
0122         delta = e - ecc * sin(e) - mrad;
0123         e -= delta / (1 - ecc * cos(e));
0124     } while (abs(delta) > epsilon);
0125     return e;
0126 }
0127 
0128 static double phaseAngle(int64_t unixmsec)
0129 {
0130     int64_t msec = unixmsec - epoch;
0131 
0132     double sunMeanAnomaly = fixAngle(msec * (360.0 / 365.2422 / 86400000.0) + elonge - elongp);
0133     double trueAnomaly = 2 * radToDeg(atan(ecPrefactor * tan(kepler(sunMeanAnomaly, earthEcc) / 2)));
0134     double geocentricEclipticLong = fixAngle(trueAnomaly + elongp);
0135 
0136     double moonMeanLong = fixAngle(msec * (13.1763966 / 86400000.0) + mmlong);
0137     double moonMeanAnomaly = fixAngle(moonMeanLong - msec * (0.1114041 / 86400000.0) - mmlongp);
0138     double evection = 1.2739 * sin(degToRad(2 * (moonMeanLong - geocentricEclipticLong) - moonMeanAnomaly));
0139     double annualEquation = 0.1858 * sin(degToRad(sunMeanAnomaly));
0140     double a3 = 0.37 * sin(degToRad(sunMeanAnomaly));
0141     double correctedAnomaly = moonMeanAnomaly + evection - annualEquation - a3;
0142     double centreCorrection = 6.2886 * sin(degToRad(correctedAnomaly));
0143     double a4 = 0.214 * sin(degToRad(2 * correctedAnomaly));
0144     double correctedLong = moonMeanLong + evection + centreCorrection - annualEquation + a4;
0145     double variation = 0.6583 * sin(degToRad(2 * (correctedLong - geocentricEclipticLong)));
0146     double trueLong = correctedLong + variation;
0147 
0148     return fixAngle(trueLong - geocentricEclipticLong);
0149 }