1 /*
2 SPDX-FileCopyrightText: 2001 Jason Harris <kstars@30doradus.org>
3
4 SPDX-License-Identifier: GPL-2.0-or-later
5 */
6
7 #include "ksmoon.h"
8
9 #include "ksnumbers.h"
10 #include "ksutils.h"
11 #include "kssun.h"
12 #include "kstarsdata.h"
13 #ifndef KSTARS_LITE
14 #include "kspopupmenu.h"
15 #endif
16 #include "skycomponents/skymapcomposite.h"
17 #include "skycomponents/solarsystemcomposite.h"
18 #include "texturemanager.h"
19
20 #include <QFile>
21 #include <QTextStream>
22
23 #include <cstdlib>
24 #include <cmath>
25 #if defined(_MSC_VER)
26 #include <float.h>
27 #endif
28
29 #include <typeinfo>
30
31 using namespace std;
32
33 namespace
34 {
35 // Convert degrees to radians and put it into [0,2*pi] range
degToRad(double x)36 double degToRad(double x)
37 {
38 return dms::DegToRad * KSUtils::reduceAngle(x, 0.0, 360.0);
39 }
40
41 /*
42 * Data used to calculate moon magnitude.
43 *
44 * Formula and data were obtained from SkyChart v3.x Beta
45 *
46 */
47 // intensities in Table 1 of M. Minnaert (1961),
48 // Phase Frac. Phase Frac. Phase Frac.
49 // angle ill. Mag angle ill. Mag angle ill. Mag
50 // 0 1.00 -12.7 60 0.75 -11.0 120 0.25 -8.7
51 // 10 0.99 -12.4 70 0.67 -10.8 130 0.18 -8.2
52 // 20 0.97 -12.1 80 0.59 -10.4 140 0.12 -7.6
53 // 30 0.93 -11.8 90 0.50 -10.0 150 0.07 -6.7
54 // 40 0.88 -11.5 100 0.41 -9.6 160 0.03 -3.4
55 // 50 0.82 -11.2 110 0.33 -9.2
56 static const double MagArray[19] = { -12.7, -12.4, -12.1, -11.8, -11.5, -11.2, -11.0, -10.8, -10.4, -10.0,
57 -9.6, -9.2, -8.7, -8.2, -7.6, -6.7, -3.4, 0, 0
58 };
59 }
60
KSMoon()61 KSMoon::KSMoon() : KSPlanetBase(i18n("Moon"), QString(), QColor("white"), 3474.8 /*diameter in km*/)
62 {
63 instance_count++;
64 //Reset object type
65 setType(SkyObject::MOON);
66 }
67
KSMoon(const KSMoon & o)68 KSMoon::KSMoon(const KSMoon &o) : KSPlanetBase(o)
69 {
70 instance_count++;
71 }
72
clone() const73 KSMoon *KSMoon::clone() const
74 {
75 Q_ASSERT(typeid(this) == typeid(static_cast<const KSMoon *>(this))); // Ensure we are not slicing a derived class
76 return new KSMoon(*this);
77 }
78
~KSMoon()79 KSMoon::~KSMoon()
80 {
81 instance_count--;
82 if (instance_count <= 0)
83 {
84 LRData.clear();
85 BData.clear();
86 data_loaded = false;
87 }
88 }
89
90 bool KSMoon::data_loaded = false;
91 int KSMoon::instance_count = 0;
92 QList<KSMoon::MoonLRData> KSMoon::LRData;
93 QList<KSMoon::MoonBData> KSMoon::BData;
94
loadData()95 bool KSMoon::loadData()
96 {
97 if (data_loaded)
98 return true;
99
100 QStringList fields;
101 QFile f;
102
103 if (KSUtils::openDataFile(f, "moonLR.dat"))
104 {
105 QTextStream stream(&f);
106 while (!stream.atEnd())
107 {
108 fields = stream.readLine().split(' ', QString::SkipEmptyParts);
109
110 if (fields.size() == 6)
111 {
112 LRData.append(MoonLRData());
113 LRData.last().nd = fields[0].toInt();
114 LRData.last().nm = fields[1].toInt();
115 LRData.last().nm1 = fields[2].toInt();
116 LRData.last().nf = fields[3].toInt();
117 LRData.last().Li = fields[4].toDouble();
118 LRData.last().Ri = fields[5].toDouble();
119 }
120 }
121 f.close();
122 }
123 else
124 return false;
125
126 if (KSUtils::openDataFile(f, "moonB.dat"))
127 {
128 QTextStream stream(&f);
129 while (!stream.atEnd())
130 {
131 fields = stream.readLine().split(' ', QString::SkipEmptyParts);
132
133 if (fields.size() == 5)
134 {
135 BData.append(MoonBData());
136 BData.last().nd = fields[0].toInt();
137 BData.last().nm = fields[1].toInt();
138 BData.last().nm1 = fields[2].toInt();
139 BData.last().nf = fields[3].toInt();
140 BData.last().Bi = fields[4].toDouble();
141 }
142 }
143 f.close();
144 }
145
146 data_loaded = true;
147 return true;
148 }
149
findGeocentricPosition(const KSNumbers * num,const KSPlanetBase *)150 bool KSMoon::findGeocentricPosition(const KSNumbers *num, const KSPlanetBase *)
151 {
152 //Algorithms in this subroutine are taken from Chapter 45 of "Astronomical Algorithms"
153 //by Jean Meeus (1991, Willmann-Bell, Inc. ISBN 0-943396-35-2. https://www.willbell.com/math/mc1.htm)
154 //updated to Jean Messus (1998, Willmann-Bell, http://www.naughter.com/aa.html )
155
156 double T, L, D, M, M1, F, A1, A2, A3;
157 double sumL, sumR, sumB;
158
159 //Julian centuries since J2000
160 T = num->julianCenturies();
161
162 double Et = 1.0 - 0.002516 * T - 0.0000074 * T * T;
163
164 //Moon's mean longitude
165 L = degToRad(218.3164477 + 481267.88123421 * T - 0.0015786 * T * T + T * T * T / 538841.0 -
166 T * T * T * T / 65194000.0);
167 //Moon's mean elongation
168 D = degToRad(297.8501921 + 445267.1114034 * T - 0.0018819 * T * T + T * T * T / 545868.0 -
169 T * T * T * T / 113065000.0);
170 //Sun's mean anomaly
171 M = degToRad(357.5291092 + 35999.0502909 * T - 0.0001536 * T * T + T * T * T / 24490000.0);
172 //Moon's mean anomaly
173 M1 = degToRad(134.9633964 + 477198.8675055 * T + 0.0087414 * T * T + T * T * T / 69699.0 -
174 T * T * T * T / 14712000.0);
175 //Moon's argument of latitude (angle from ascending node)
176 F = degToRad(93.2720950 + 483202.0175233 * T - 0.0036539 * T * T - T * T * T / 3526000.0 +
177 T * T * T * T / 863310000.0);
178
179 A1 = degToRad(119.75 + 131.849 * T);
180 A2 = degToRad(53.09 + 479264.290 * T);
181 A3 = degToRad(313.45 + 481266.484 * T);
182
183 //Calculate the series expansions stored in moonLR.txt and moonB.txt.
184 //
185 sumL = 0.0;
186 sumR = 0.0;
187
188 if (!loadData())
189 return false;
190
191 for (const auto &mlrd : LRData)
192 {
193 double E = 1.0;
194
195 if (mlrd.nm) //if M != 0, include changing eccentricity of Earth's orbit
196 {
197 E = Et;
198 if (abs(mlrd.nm) == 2)
199 E = E * E; //use E^2
200 }
201 sumL += E * mlrd.Li * sin(mlrd.nd * D + mlrd.nm * M + mlrd.nm1 * M1 + mlrd.nf * F);
202 sumR += E * mlrd.Ri * cos(mlrd.nd * D + mlrd.nm * M + mlrd.nm1 * M1 + mlrd.nf * F);
203 }
204
205 sumB = 0.0;
206 for (const auto &mbd : BData)
207 {
208 double E = 1.0;
209
210 if (mbd.nm) //if M != 0, include changing eccentricity of Earth's orbit
211 {
212 E = Et;
213 if (abs(mbd.nm) == 2)
214 E = E * E; //use E^2
215 }
216 sumB += E * mbd.Bi * sin(mbd.nd * D + mbd.nm * M + mbd.nm1 * M1 + mbd.nf * F);
217 }
218
219 //Additive terms for sumL and sumB
220 sumL += (3958.0 * sin(A1) + 1962.0 * sin(L - F) + 318.0 * sin(A2));
221 sumB += (-2235.0 * sin(L) + 382.0 * sin(A3) + 175.0 * sin(A1 - F) + 175.0 * sin(A1 + F) + 127.0 * sin(L - M1) -
222 115.0 * sin(L + M1));
223
224 //Geocentric coordinates
225 setEcLong(dms(sumL / 1000000.0 + L * 180.0 / dms::PI)); //convert radians to degrees
226 setEcLat(dms(sumB / 1000000.0));
227 Rearth = (385000.56 + sumR / 1000.0) / AU_KM; //distance from Earth, in AU
228
229 EclipticToEquatorial(num->obliquity());
230
231 //Determine position angle
232 findPA(num);
233
234 return true;
235 }
236
findMagnitude(const KSNumbers *)237 void KSMoon::findMagnitude(const KSNumbers *)
238 {
239 // This block of code to compute Moon magnitude (and the
240 // relevant data put into ksplanetbase.h) was taken from
241 // SkyChart v3 Beta
242 double phd = phase().Degrees();
243
244 if (std::isnan(phd)) // Avoid nanny phases.
245 {
246 findPhase(nullptr);
247 phd = phase().Degrees();
248 if (std::isnan(phd))
249 return;
250 }
251 int p = floor(phd);
252 if (p > 180)
253 p = p - 360;
254 int i = p / 10;
255 int k = p % 10;
256 int j = (i + 1 > 18) ? 18 : i + 1;
257 i = 18 - abs(i);
258 j = 18 - abs(j);
259 if (i >= 0 && j >= 0 && i <= 18 && j <= 18)
260 {
261 setMag(MagArray[i] + (MagArray[j] - MagArray[i]) * k / 10);
262 }
263 }
264
findPhase(const KSSun * Sun)265 void KSMoon::findPhase(const KSSun *Sun)
266 {
267 if (Sun == nullptr)
268 {
269 if (defaultSun == nullptr)
270 defaultSun = KStarsData::Instance()->skyComposite()->solarSystemComposite()->sun();
271 Sun = defaultSun;
272 }
273
274 Q_ASSERT(Sun != nullptr);
275
276 // This is an approximation justified by the small Earth-Moon distance in relation
277 // to the great Earth-Sun distance
278 Phase = (ecLong() - Sun->ecLong()).Degrees(); // Phase is obviously in degrees
279 double DegPhase = dms(Phase).reduce().Degrees();
280 iPhase = int(0.1 * DegPhase + 0.5) % 36; // iPhase must be in [0,36) range
281
282 m_image = TextureManager::getImage(QString("moon%1").arg(iPhase, 2, 10, QChar('0')));
283 }
284
phaseName() const285 QString KSMoon::phaseName() const
286 {
287 double f = illum();
288 double p = abs(dms(Phase).reduce().Degrees());
289
290 //First, handle the major phases
291 if (f > 0.99)
292 return i18nc("moon phase, 100 percent illuminated", "Full moon");
293 if (f < 0.01)
294 return i18nc("moon phase, 0 percent illuminated", "New moon");
295 if (fabs(f - 0.50) < 0.06)
296 {
297 if (p < 180.0)
298 return i18nc("moon phase, half-illuminated and growing", "First quarter");
299 else
300 return i18nc("moon phase, half-illuminated and shrinking", "Third quarter");
301 }
302
303 //Next, handle the more general cases
304 if (p < 90.0)
305 return i18nc("moon phase between new moon and 1st quarter", "Waxing crescent");
306 else if (p < 180.0)
307 return i18nc("moon phase between 1st quarter and full moon", "Waxing gibbous");
308 else if (p < 270.0)
309 return i18nc("moon phase between full moon and 3rd quarter", "Waning gibbous");
310 else if (p < 360.0)
311 return i18nc("moon phase between 3rd quarter and new moon", "Waning crescent");
312
313 else
314 return i18n("unknown");
315 }
316
initPopupMenu(KSPopupMenu * pmenu)317 void KSMoon::initPopupMenu(KSPopupMenu *pmenu)
318 {
319 #ifdef KSTARS_LITE
320 Q_UNUSED(pmenu)
321 #else
322 pmenu->createMoonMenu(this);
323 #endif
324 }
325
getUID() const326 SkyObject::UID KSMoon::getUID() const
327 {
328 return solarsysUID(UID_SOL_BIGOBJ) | 10;
329 }
330