File indexing completed on 2024-04-28 11:25:19
0001 /* 0002 SPDX-FileCopyrightText: 2001 Jason Harris <kstars@30doradus.org> 0003 0004 SPDX-License-Identifier: GPL-2.0-or-later 0005 */ 0006 0007 #include "ksmoon.h" 0008 0009 #include "ksnumbers.h" 0010 #include "ksutils.h" 0011 #include "kssun.h" 0012 #include "kstarsdata.h" 0013 #ifndef KSTARS_LITE 0014 #include "kspopupmenu.h" 0015 #endif 0016 #include "skycomponents/skymapcomposite.h" 0017 #include "skycomponents/solarsystemcomposite.h" 0018 #include "texturemanager.h" 0019 0020 #include <QFile> 0021 #include <QTextStream> 0022 0023 #include <cstdlib> 0024 #include <cmath> 0025 #if defined(_MSC_VER) 0026 #include <float.h> 0027 #endif 0028 0029 #include <typeinfo> 0030 0031 // Qt version calming 0032 #include <qtskipemptyparts.h> 0033 0034 using namespace std; 0035 0036 namespace 0037 { 0038 // Convert degrees to radians and put it into [0,2*pi] range 0039 double degToRad(double x) 0040 { 0041 return dms::DegToRad * KSUtils::reduceAngle(x, 0.0, 360.0); 0042 } 0043 0044 /* 0045 * Data used to calculate moon magnitude. 0046 * 0047 * Formula and data were obtained from SkyChart v3.x Beta 0048 * 0049 */ 0050 // intensities in Table 1 of M. Minnaert (1961), 0051 // Phase Frac. Phase Frac. Phase Frac. 0052 // angle ill. Mag angle ill. Mag angle ill. Mag 0053 // 0 1.00 -12.7 60 0.75 -11.0 120 0.25 -8.7 0054 // 10 0.99 -12.4 70 0.67 -10.8 130 0.18 -8.2 0055 // 20 0.97 -12.1 80 0.59 -10.4 140 0.12 -7.6 0056 // 30 0.93 -11.8 90 0.50 -10.0 150 0.07 -6.7 0057 // 40 0.88 -11.5 100 0.41 -9.6 160 0.03 -3.4 0058 // 50 0.82 -11.2 110 0.33 -9.2 0059 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, 0060 -9.6, -9.2, -8.7, -8.2, -7.6, -6.7, -3.4, 0, 0 0061 }; 0062 } 0063 0064 KSMoon::KSMoon() : KSPlanetBase(i18n("Moon"), QString(), QColor("white"), 3474.8 /*diameter in km*/) 0065 { 0066 instance_count++; 0067 //Reset object type 0068 setType(SkyObject::MOON); 0069 } 0070 0071 KSMoon::KSMoon(const KSMoon &o) : KSPlanetBase(o) 0072 { 0073 instance_count++; 0074 } 0075 0076 KSMoon *KSMoon::clone() const 0077 { 0078 Q_ASSERT(typeid(this) == typeid(static_cast<const KSMoon *>(this))); // Ensure we are not slicing a derived class 0079 return new KSMoon(*this); 0080 } 0081 0082 KSMoon::~KSMoon() 0083 { 0084 instance_count--; 0085 if (instance_count <= 0) 0086 { 0087 LRData.clear(); 0088 BData.clear(); 0089 data_loaded = false; 0090 } 0091 } 0092 0093 bool KSMoon::data_loaded = false; 0094 int KSMoon::instance_count = 0; 0095 QList<KSMoon::MoonLRData> KSMoon::LRData; 0096 QList<KSMoon::MoonBData> KSMoon::BData; 0097 0098 bool KSMoon::loadData() 0099 { 0100 if (data_loaded) 0101 return true; 0102 0103 QStringList fields; 0104 QFile f; 0105 0106 if (KSUtils::openDataFile(f, "moonLR.dat")) 0107 { 0108 QTextStream stream(&f); 0109 while (!stream.atEnd()) 0110 { 0111 fields = stream.readLine().split(' ', Qt::SkipEmptyParts); 0112 0113 if (fields.size() == 6) 0114 { 0115 LRData.append(MoonLRData()); 0116 LRData.last().nd = fields[0].toInt(); 0117 LRData.last().nm = fields[1].toInt(); 0118 LRData.last().nm1 = fields[2].toInt(); 0119 LRData.last().nf = fields[3].toInt(); 0120 LRData.last().Li = fields[4].toDouble(); 0121 LRData.last().Ri = fields[5].toDouble(); 0122 } 0123 } 0124 f.close(); 0125 } 0126 else 0127 return false; 0128 0129 if (KSUtils::openDataFile(f, "moonB.dat")) 0130 { 0131 QTextStream stream(&f); 0132 while (!stream.atEnd()) 0133 { 0134 fields = stream.readLine().split(' ', Qt::SkipEmptyParts); 0135 0136 if (fields.size() == 5) 0137 { 0138 BData.append(MoonBData()); 0139 BData.last().nd = fields[0].toInt(); 0140 BData.last().nm = fields[1].toInt(); 0141 BData.last().nm1 = fields[2].toInt(); 0142 BData.last().nf = fields[3].toInt(); 0143 BData.last().Bi = fields[4].toDouble(); 0144 } 0145 } 0146 f.close(); 0147 } 0148 0149 data_loaded = true; 0150 return true; 0151 } 0152 0153 bool KSMoon::findGeocentricPosition(const KSNumbers *num, const KSPlanetBase *) 0154 { 0155 //Algorithms in this subroutine are taken from Chapter 45 of "Astronomical Algorithms" 0156 //by Jean Meeus (1991, Willmann-Bell, Inc. ISBN 0-943396-35-2. https://www.willbell.com/math/mc1.htm) 0157 //updated to Jean Messus (1998, Willmann-Bell, http://www.naughter.com/aa.html ) 0158 0159 double T, L, D, M, M1, F, A1, A2, A3; 0160 double sumL, sumR, sumB; 0161 0162 //Julian centuries since J2000 0163 T = num->julianCenturies(); 0164 0165 double Et = 1.0 - 0.002516 * T - 0.0000074 * T * T; 0166 0167 //Moon's mean longitude 0168 L = degToRad(218.3164477 + 481267.88123421 * T - 0.0015786 * T * T + T * T * T / 538841.0 - 0169 T * T * T * T / 65194000.0); 0170 //Moon's mean elongation 0171 D = degToRad(297.8501921 + 445267.1114034 * T - 0.0018819 * T * T + T * T * T / 545868.0 - 0172 T * T * T * T / 113065000.0); 0173 //Sun's mean anomaly 0174 M = degToRad(357.5291092 + 35999.0502909 * T - 0.0001536 * T * T + T * T * T / 24490000.0); 0175 //Moon's mean anomaly 0176 M1 = degToRad(134.9633964 + 477198.8675055 * T + 0.0087414 * T * T + T * T * T / 69699.0 - 0177 T * T * T * T / 14712000.0); 0178 //Moon's argument of latitude (angle from ascending node) 0179 F = degToRad(93.2720950 + 483202.0175233 * T - 0.0036539 * T * T - T * T * T / 3526000.0 + 0180 T * T * T * T / 863310000.0); 0181 0182 A1 = degToRad(119.75 + 131.849 * T); 0183 A2 = degToRad(53.09 + 479264.290 * T); 0184 A3 = degToRad(313.45 + 481266.484 * T); 0185 0186 //Calculate the series expansions stored in moonLR.txt and moonB.txt. 0187 // 0188 sumL = 0.0; 0189 sumR = 0.0; 0190 0191 if (!loadData()) 0192 return false; 0193 0194 for (const auto &mlrd : LRData) 0195 { 0196 double E = 1.0; 0197 0198 if (mlrd.nm) //if M != 0, include changing eccentricity of Earth's orbit 0199 { 0200 E = Et; 0201 if (abs(mlrd.nm) == 2) 0202 E = E * E; //use E^2 0203 } 0204 sumL += E * mlrd.Li * sin(mlrd.nd * D + mlrd.nm * M + mlrd.nm1 * M1 + mlrd.nf * F); 0205 sumR += E * mlrd.Ri * cos(mlrd.nd * D + mlrd.nm * M + mlrd.nm1 * M1 + mlrd.nf * F); 0206 } 0207 0208 sumB = 0.0; 0209 for (const auto &mbd : BData) 0210 { 0211 double E = 1.0; 0212 0213 if (mbd.nm) //if M != 0, include changing eccentricity of Earth's orbit 0214 { 0215 E = Et; 0216 if (abs(mbd.nm) == 2) 0217 E = E * E; //use E^2 0218 } 0219 sumB += E * mbd.Bi * sin(mbd.nd * D + mbd.nm * M + mbd.nm1 * M1 + mbd.nf * F); 0220 } 0221 0222 //Additive terms for sumL and sumB 0223 sumL += (3958.0 * sin(A1) + 1962.0 * sin(L - F) + 318.0 * sin(A2)); 0224 sumB += (-2235.0 * sin(L) + 382.0 * sin(A3) + 175.0 * sin(A1 - F) + 175.0 * sin(A1 + F) + 127.0 * sin(L - M1) - 0225 115.0 * sin(L + M1)); 0226 0227 //Geocentric coordinates 0228 setEcLong(dms(sumL / 1000000.0 + L * 180.0 / dms::PI)); //convert radians to degrees 0229 setEcLat(dms(sumB / 1000000.0)); 0230 Rearth = (385000.56 + sumR / 1000.0) / AU_KM; //distance from Earth, in AU 0231 0232 EclipticToEquatorial(num->obliquity()); 0233 0234 //Determine position angle 0235 findPA(num); 0236 0237 return true; 0238 } 0239 0240 void KSMoon::findMagnitude(const KSNumbers *) 0241 { 0242 // This block of code to compute Moon magnitude (and the 0243 // relevant data put into ksplanetbase.h) was taken from 0244 // SkyChart v3 Beta 0245 double phd = phase().Degrees(); 0246 0247 if (std::isnan(phd)) // Avoid nanny phases. 0248 { 0249 findPhase(nullptr); 0250 phd = phase().Degrees(); 0251 if (std::isnan(phd)) 0252 return; 0253 } 0254 int p = floor(phd); 0255 if (p > 180) 0256 p = p - 360; 0257 int i = p / 10; 0258 int k = p % 10; 0259 int j = (i + 1 > 18) ? 18 : i + 1; 0260 i = 18 - abs(i); 0261 j = 18 - abs(j); 0262 if (i >= 0 && j >= 0 && i <= 18 && j <= 18) 0263 { 0264 setMag(MagArray[i] + (MagArray[j] - MagArray[i]) * k / 10); 0265 } 0266 } 0267 0268 void KSMoon::findPhase(const KSSun *Sun) 0269 { 0270 if (Sun == nullptr) 0271 { 0272 if (defaultSun == nullptr) 0273 defaultSun = KStarsData::Instance()->skyComposite()->solarSystemComposite()->sun(); 0274 Sun = defaultSun; 0275 } 0276 0277 Q_ASSERT(Sun != nullptr); 0278 0279 // This is an approximation justified by the small Earth-Moon distance in relation 0280 // to the great Earth-Sun distance 0281 Phase = (ecLong() - Sun->ecLong()).Degrees(); // Phase is obviously in degrees 0282 double DegPhase = dms(Phase).reduce().Degrees(); 0283 iPhase = int(0.1 * DegPhase + 0.5) % 36; // iPhase must be in [0,36) range 0284 0285 m_image = TextureManager::getImage(QString("moon%1").arg(iPhase, 2, 10, QChar('0'))); 0286 } 0287 0288 QString KSMoon::phaseName() const 0289 { 0290 double f = illum(); 0291 double p = abs(dms(Phase).reduce().Degrees()); 0292 0293 //First, handle the major phases 0294 if (f > 0.99) 0295 return i18nc("moon phase, 100 percent illuminated", "Full moon"); 0296 if (f < 0.01) 0297 return i18nc("moon phase, 0 percent illuminated", "New moon"); 0298 if (fabs(f - 0.50) < 0.06) 0299 { 0300 if (p < 180.0) 0301 return i18nc("moon phase, half-illuminated and growing", "First quarter"); 0302 else 0303 return i18nc("moon phase, half-illuminated and shrinking", "Third quarter"); 0304 } 0305 0306 //Next, handle the more general cases 0307 if (p < 90.0) 0308 return i18nc("moon phase between new moon and 1st quarter", "Waxing crescent"); 0309 else if (p < 180.0) 0310 return i18nc("moon phase between 1st quarter and full moon", "Waxing gibbous"); 0311 else if (p < 270.0) 0312 return i18nc("moon phase between full moon and 3rd quarter", "Waning gibbous"); 0313 else if (p < 360.0) 0314 return i18nc("moon phase between 3rd quarter and new moon", "Waning crescent"); 0315 0316 else 0317 return i18n("unknown"); 0318 } 0319 0320 void KSMoon::initPopupMenu(KSPopupMenu *pmenu) 0321 { 0322 #ifdef KSTARS_LITE 0323 Q_UNUSED(pmenu) 0324 #else 0325 pmenu->createMoonMenu(this); 0326 #endif 0327 } 0328 0329 SkyObject::UID KSMoon::getUID() const 0330 { 0331 return solarsysUID(UID_SOL_BIGOBJ) | 10; 0332 }