KHolidays

lunarphase.cpp
1 /*
2  This file is part of the kholidays library.
3 
4  SPDX-FileCopyrightText: 2004, 2007, 2009 Allen Winter <[email protected]>
5 
6  SPDX-License-Identifier: LGPL-2.0-or-later
7 */
8 
9 #include "lunarphase.h"
10 #include <config-kholidays.h>
11 
12 #include <QCoreApplication>
13 #include <QDateTime>
14 
15 using namespace KHolidays;
16 
18 {
19  return phaseName(phaseAtDate(date));
20 }
21 
23 {
24  switch (phase) {
25  case NewMoon:
26  return (QCoreApplication::translate("LunarPhase", "New Moon"));
27  case FullMoon:
28  return (QCoreApplication::translate("LunarPhase", "Full Moon"));
29  case FirstQuarter:
30  return (QCoreApplication::translate("LunarPhase", "First Quarter Moon"));
31  case LastQuarter:
32  return (QCoreApplication::translate("LunarPhase", "Last Quarter Moon"));
33  case None:
34  return QString();
35  }
36  return QString();
37 }
38 
39 static double phaseAngle(int64_t unixmsec);
40 
42 {
43  Phase retPhase = None;
44 
45  const QTime midnight(0, 0, 0);
46  const QDateTime todayStart(date, midnight, Qt::UTC);
47  const double startAngle = phaseAngle(todayStart.toMSecsSinceEpoch());
48  const QDateTime todayEnd(date.addDays(1), midnight, Qt::UTC);
49  const double endAngle = phaseAngle(todayEnd.toMSecsSinceEpoch());
50 
51  if (startAngle > endAngle)
52  retPhase = NewMoon;
53  else if (startAngle < 90.0 && endAngle > 90.0)
54  retPhase = FirstQuarter;
55  else if (startAngle < 180.0 && endAngle > 180.0)
56  retPhase = FullMoon;
57  else if (startAngle < 270.0 && endAngle > 270.0)
58  retPhase = LastQuarter;
59 
60  return retPhase;
61 }
62 
63 /*
64  Algorithm adapted from Moontool by John Walker.
65  Moontool is in the public domain: "Do what thou wilt shall be the whole of the law".
66  Unnecessary calculations have been removed.
67  See https://fourmilab.ch/moontool .
68 */
69 
70 #include <cmath>
71 
72 constexpr int64_t epoch = 315446400000; // "1980 January 0.0", a.k.a. midnight on 1979-12-31, in ms Unix time
73 constexpr double elonge = 278.833540; // ecliptic longitude of sun at epoch
74 constexpr double elongp = 282.596403; // ecliptic longitude of sun at perigee
75 constexpr double earthEcc = 0.016718; // eccentricity of earth orbit
76 static const double ecPrefactor = sqrt((1 + earthEcc) / (1 - earthEcc));
77 
78 constexpr double mmlong = 64.975464; // mean longitude of moon at epoch
79 constexpr double mmlongp = 349.383063; // mean longitude of moon at perigee
80 
81 constexpr double PI = 3.14159265358979323846;
82 
83 static double fixAngle(double degrees)
84 {
85  return degrees - floor(degrees / 360.0) * 360.0;
86 }
87 
88 static constexpr double radToDeg(double rad)
89 {
90  return rad / PI * 180.0;
91 }
92 
93 static constexpr double degToRad(double deg)
94 {
95  return deg / 180.0 * PI;
96 }
97 
98 constexpr double epsilon = 1e-6;
99 
100 static double kepler(double m, double ecc)
101 {
102  double mrad = degToRad(m);
103  double e = mrad;
104  double delta;
105  do {
106  delta = e - ecc * sin(e) - mrad;
107  e -= delta / (1 - ecc * cos(e));
108  } while (abs(delta) > epsilon);
109  return e;
110 }
111 
112 static double phaseAngle(int64_t unixmsec)
113 {
114  int64_t msec = unixmsec - epoch;
115 
116  double sunMeanAnomaly = fixAngle(msec * (360.0 / 365.2422 / 86400000.0) + elonge - elongp);
117  double trueAnomaly = 2 * radToDeg(atan(ecPrefactor * tan(kepler(sunMeanAnomaly, earthEcc) / 2)));
118  double geocentricEclipticLong = fixAngle(trueAnomaly + elongp);
119 
120  double moonMeanLong = fixAngle(msec * (13.1763966 / 86400000.0) + mmlong);
121  double moonMeanAnomaly = fixAngle(moonMeanLong - msec * (0.1114041 / 86400000.0) - mmlongp);
122  double evection = 1.2739 * sin(degToRad(2 * (moonMeanLong - geocentricEclipticLong) - moonMeanAnomaly));
123  double annualEquation = 0.1858 * sin(degToRad(sunMeanAnomaly));
124  double a3 = 0.37 * sin(degToRad(sunMeanAnomaly));
125  double correctedAnomaly = moonMeanAnomaly + evection - annualEquation - a3;
126  double centreCorrection = 6.2886 * sin(degToRad(correctedAnomaly));
127  double a4 = 0.214 * sin(degToRad(2 * correctedAnomaly));
128  double correctedLong = moonMeanLong + evection + centreCorrection - annualEquation + a4;
129  double variation = 0.6583 * sin(degToRad(2 * (correctedLong - geocentricEclipticLong)));
130  double trueLong = correctedLong + variation;
131 
132  return fixAngle(trueLong - geocentricEclipticLong);
133 }
Last quarter of moon phase.
Definition: lunarphase.h:52
qint64 toMSecsSinceEpoch() const const
QString translate(const char *context, const char *sourceText, const char *disambiguation, int n)
First quarter of moon phase.
Definition: lunarphase.h:51
static QString phaseName(Phase phase)
Return the string representation of phase.
Definition: lunarphase.cpp:22
Indication for error.
Definition: lunarphase.h:54
static QString phaseNameAtDate(const QDate &date)
Return the lunar phase as a text string for the specified date.
Definition: lunarphase.cpp:17
Phase
Phases of the moon, in traditional English notation.
Definition: lunarphase.h:49
static Phase phaseAtDate(const QDate &date)
Return the lunar phase for the specified Gregorian date.
Definition: lunarphase.cpp:41
QDate addDays(qint64 ndays) const const
This file is part of the KDE documentation.
Documentation copyright © 1996-2021 The KDE developers.
Generated on Tue May 11 2021 22:56:53 by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.