1 /*
2     SPDX-FileCopyrightText: 2001 Jason Harris <jharris@30doradus.org>
3 
4     SPDX-License-Identifier: GPL-2.0-or-later
5 */
6 
7 #include "ksplanet.h"
8 
9 #include "ksnumbers.h"
10 #include "ksutils.h"
11 #include "ksfilereader.h"
12 
13 #include <cmath>
14 #include <typeinfo>
15 
16 #include "kstars_debug.h"
17 
18 KSPlanet::OrbitDataManager KSPlanet::odm;
19 
OrbitDataManager()20 KSPlanet::OrbitDataManager::OrbitDataManager()
21 {
22     //EMPTY
23 }
24 
readOrbitData(const QString & fname,QVector<OrbitData> * vector)25 bool KSPlanet::OrbitDataManager::readOrbitData(const QString &fname, QVector<OrbitData> *vector)
26 {
27     QFile f;
28 
29     if (KSUtils::openDataFile(f, fname))
30     {
31         KSFileReader fileReader(f); // close file is included
32         QStringList fields;
33         while (fileReader.hasMoreLines())
34         {
35             fields = fileReader.readLine().split(' ', QString::SkipEmptyParts);
36 
37             if (fields.size() == 3)
38             {
39                 double A = fields[0].toDouble();
40                 double B = fields[1].toDouble();
41                 double C = fields[2].toDouble();
42                 vector->append(OrbitData(A, B, C));
43             }
44         }
45     }
46     else
47     {
48         return false;
49     }
50     return true;
51 }
52 
loadData(KSPlanet::OrbitDataColl & odc,const QString & n)53 bool KSPlanet::OrbitDataManager::loadData(KSPlanet::OrbitDataColl &odc, const QString &n)
54 {
55     QString fname, snum;
56     QFile f;
57     int nCount = 0;
58     QString nl = n.toLower();
59 
60     if (hash.contains(nl))
61     {
62         odc = hash[nl];
63         return true; //orbit data already loaded
64     }
65 
66     //Create a new OrbitDataColl
67     OrbitDataColl ret;
68 
69     //Ecliptic Longitude
70     for (int i = 0; i < 6; ++i)
71     {
72         snum.setNum(i);
73         fname = nl + ".L" + snum + ".vsop";
74         if (readOrbitData(fname, &ret.Lon[i]))
75             nCount++;
76     }
77 
78     if (nCount == 0)
79         return false;
80 
81     //Ecliptic Latitude
82     for (int i = 0; i < 6; ++i)
83     {
84         snum.setNum(i);
85         fname = nl + ".B" + snum + ".vsop";
86         if (readOrbitData(fname, &ret.Lat[i]))
87             nCount++;
88     }
89 
90     if (nCount == 0)
91         return false;
92 
93     //Heliocentric Distance
94     for (int i = 0; i < 6; ++i)
95     {
96         snum.setNum(i);
97         fname = nl + ".R" + snum + ".vsop";
98         if (readOrbitData(fname, &ret.Dst[i]))
99             nCount++;
100     }
101 
102     if (nCount == 0)
103         return false;
104 
105     hash[nl] = ret;
106     odc      = hash[nl];
107 
108     return true;
109 }
110 
KSPlanet(const QString & s,const QString & imfile,const QColor & c,double pSize)111 KSPlanet::KSPlanet(const QString &s, const QString &imfile, const QColor &c, double pSize)
112     : KSPlanetBase(s, imfile, c, pSize)
113 {
114 }
115 
KSPlanet(int n)116 KSPlanet::KSPlanet(int n) : KSPlanetBase()
117 {
118     switch (n)
119     {
120         case MERCURY:
121             KSPlanetBase::init(i18n("Mercury"), "mercury", KSPlanetBase::planetColor[KSPlanetBase::MERCURY], 4879.4);
122             break;
123         case VENUS:
124             KSPlanetBase::init(i18n("Venus"), "venus", KSPlanetBase::planetColor[KSPlanetBase::VENUS], 12103.6);
125             break;
126         case MARS:
127             KSPlanetBase::init(i18n("Mars"), "mars", KSPlanetBase::planetColor[KSPlanetBase::MARS], 6792.4);
128             break;
129         case JUPITER:
130             KSPlanetBase::init(i18n("Jupiter"), "jupiter", KSPlanetBase::planetColor[KSPlanetBase::JUPITER], 142984.);
131             break;
132         case SATURN:
133             KSPlanetBase::init(i18n("Saturn"), "saturn", KSPlanetBase::planetColor[KSPlanetBase::SATURN], 120536.);
134             break;
135         case URANUS:
136             KSPlanetBase::init(i18n("Uranus"), "uranus", KSPlanetBase::planetColor[KSPlanetBase::URANUS], 51118.);
137             break;
138         case NEPTUNE:
139             KSPlanetBase::init(i18n("Neptune"), "neptune", KSPlanetBase::planetColor[KSPlanetBase::NEPTUNE], 49572.);
140             break;
141         default:
142             qDebug() << "Error: Illegal identifier in KSPlanet constructor: " << n;
143             break;
144     }
145 }
146 
clone() const147 KSPlanet *KSPlanet::clone() const
148 {
149     Q_ASSERT(typeid(this) == typeid(static_cast<const KSPlanet *>(this))); // Ensure we are not slicing a derived class
150     return new KSPlanet(*this);
151 }
152 
153 // TODO: Get rid of this dirty hack post KDE 4.2 release
untranslatedName() const154 QString KSPlanet::untranslatedName() const
155 {
156     if (name() == i18n("Mercury"))
157         return "Mercury";
158     else if (name() == i18n("Venus"))
159         return "Venus";
160     else if (name() == i18n("Mars"))
161         return "Mars";
162     else if (name() == i18n("Jupiter"))
163         return "Jupiter";
164     else if (name() == i18n("Saturn"))
165         return "Saturn";
166     else if (name() == i18n("Uranus"))
167         return "Uranus";
168     else if (name() == i18n("Neptune"))
169         return "Neptune";
170     else if (name() == i18n("Earth"))
171         return "Earth";
172     else if (name() == i18n("Earth Shadow"))
173         return "Earth Shadow";
174     else if (name() == i18n("Moon"))
175         return "Moon";
176     else if (name() == i18n("Sun"))
177         return "Sun";
178     else
179         return name();
180 }
181 
182 //we don't need the reference to the ODC, so just give it a junk variable
loadData()183 bool KSPlanet::loadData()
184 {
185     OrbitDataColl odc;
186     return odm.loadData(odc, untranslatedName());
187 }
188 
calcEcliptic(double Tau,EclipticPosition & epret) const189 void KSPlanet::calcEcliptic(double Tau, EclipticPosition &epret) const
190 {
191     double sum[6];
192     OrbitDataColl odc;
193     double Tpow[6];
194 
195     Tpow[0] = 1.0;
196     for (int i = 1; i < 6; ++i)
197     {
198         Tpow[i] = Tpow[i - 1] * Tau;
199     }
200 
201     if (!odm.loadData(odc, untranslatedName()))
202     {
203         epret.longitude = dms(0.0);
204         epret.latitude  = dms(0.0);
205         epret.radius    = 0.0;
206         qCWarning(KSTARS) << "Could not get data for name:" << name() << "(" << untranslatedName() << ")";
207         return;
208     }
209 
210     //Ecliptic Longitude
211     for (int i = 0; i < 6; ++i)
212     {
213         sum[i] = 0.0;
214         for (int j = 0; j < odc.Lon[i].size(); ++j)
215         {
216             sum[i] += odc.Lon[i][j].A * cos(odc.Lon[i][j].B + odc.Lon[i][j].C * Tau);
217             /*
218             qDebug() << "sum[" << i <<"] =" << sum[i] <<
219                 " A = " << odc.Lon[i][j].A << " B = " << odc.Lon[i][j].B <<
220                 " C = " << odc.Lon[i][j].C;
221                 */
222         }
223         sum[i] *= Tpow[i];
224         //qDebug() << name() << " : sum[" << i << "] = " << sum[i];
225     }
226 
227     epret.longitude.setRadians(sum[0] + sum[1] + sum[2] + sum[3] + sum[4] + sum[5]);
228     epret.longitude.setD(epret.longitude.reduce().Degrees());
229 
230     //Compute Ecliptic Latitude
231     for (uint i = 0; i < 6; ++i)
232     {
233         sum[i] = 0.0;
234         for (int j = 0; j < odc.Lat[i].size(); ++j)
235         {
236             sum[i] += odc.Lat[i][j].A * cos(odc.Lat[i][j].B + odc.Lat[i][j].C * Tau);
237         }
238         sum[i] *= Tpow[i];
239     }
240 
241     epret.latitude.setRadians(sum[0] + sum[1] + sum[2] + sum[3] + sum[4] + sum[5]);
242 
243     //Compute Heliocentric Distance
244     for (uint i = 0; i < 6; ++i)
245     {
246         sum[i] = 0.0;
247         for (int j = 0; j < odc.Dst[i].size(); ++j)
248         {
249             sum[i] += odc.Dst[i][j].A * cos(odc.Dst[i][j].B + odc.Dst[i][j].C * Tau);
250         }
251         sum[i] *= Tpow[i];
252     }
253 
254     epret.radius = sum[0] + sum[1] + sum[2] + sum[3] + sum[4] + sum[5];
255 
256     /*
257     qDebug() << name() << " pre: Lat = " << epret.latitude.toDMSString() << " Long = " <<
258         epret.longitude.toDMSString() << " Dist = " << epret.radius;
259     */
260 }
261 
findGeocentricPosition(const KSNumbers * num,const KSPlanetBase * Earth)262 bool KSPlanet::findGeocentricPosition(const KSNumbers *num, const KSPlanetBase *Earth)
263 {
264     if (Earth != nullptr)
265     {
266         double sinL, sinL0, sinB, sinB0;
267         double cosL, cosL0, cosB, cosB0;
268         double x = 0.0, y = 0.0, z = 0.0;
269 
270         double olddst = -1000;
271         double dst    = 0;
272 
273         EclipticPosition trialpos;
274 
275         double jm = num->julianMillenia();
276 
277         Earth->ecLong().SinCos(sinL0, cosL0);
278         Earth->ecLat().SinCos(sinB0, cosB0);
279 
280         double eX = Earth->rsun() * cosB0 * cosL0;
281         double eY = Earth->rsun() * cosB0 * sinL0;
282         double eZ = Earth->rsun() * sinB0;
283 
284         bool once = true;
285         while (fabs(dst - olddst) > .001)
286         {
287             calcEcliptic(jm, trialpos);
288 
289             // We store the heliocentric ecliptic coordinates the first time they are computed.
290             if (once)
291             {
292                 helEcPos = trialpos;
293                 once     = false;
294             }
295 
296             olddst = dst;
297 
298             trialpos.longitude.SinCos(sinL, cosL);
299             trialpos.latitude.SinCos(sinB, cosB);
300 
301             x = trialpos.radius * cosB * cosL - eX;
302             y = trialpos.radius * cosB * sinL - eY;
303             z = trialpos.radius * sinB - eZ;
304 
305             //distance from Earth
306             dst = sqrt(x * x + y * y + z * z);
307 
308             //The light-travel time delay, in millenia
309             //0.0057755183 is the inverse speed of light,
310             //in days/AU
311             double delay = (.0057755183 * dst) / 365250.0;
312 
313             jm = num->julianMillenia() - delay;
314         }
315 
316         ep.longitude.setRadians(atan2(y, x));
317         ep.longitude.reduce();
318         ep.latitude.setRadians(atan2(z, sqrt(x * x + y * y)));
319         setRsun(trialpos.radius);
320         setRearth(dst);
321 
322         EclipticToEquatorial(num->obliquity());
323 
324         setRA0(ra());
325         setDec0(dec());
326         apparentCoord(J2000, lastPrecessJD);
327 
328         //nutate(num);
329         //aberrate(num);
330     }
331     else
332     {
333         calcEcliptic(num->julianMillenia(), ep);
334         helEcPos = ep;
335     }
336 
337     //determine the position angle
338     findPA(num);
339 
340     return true;
341 }
342 
findMagnitude(const KSNumbers * num)343 void KSPlanet::findMagnitude(const KSNumbers *num)
344 {
345     double cosDec, sinDec;
346     dec().SinCos(cosDec, sinDec);
347 
348     /* Computation of the visual magnitude (V band) of the planet.
349     * Algorithm provided by Pere Planesas (Observatorio Astronomico Nacional)
350     * It has some similarity to J. Meeus algorithm in Astronomical Algorithms, Chapter 40.
351     * */
352 
353     // Initialized to the faintest magnitude observable with the HST
354     float magnitude = 30;
355 
356     double param = 5 * log10(rsun() * rearth());
357     double phase = this->phase().Degrees();
358     double f1    = phase / 100.;
359 
360     if (name() == i18n("Mercury"))
361     {
362         if (phase > 150.)
363             f1 = 1.5;
364         magnitude = -0.36 + param + 3.8 * f1 - 2.73 * f1 * f1 + 2 * f1 * f1 * f1;
365     }
366     else if (name() == i18n("Venus"))
367     {
368         magnitude = -4.29 + param + 0.09 * f1 + 2.39 * f1 * f1 - 0.65 * f1 * f1 * f1;
369     }
370     else if (name() == i18n("Mars"))
371     {
372         magnitude = -1.52 + param + 0.016 * phase;
373     }
374     else if (name() == i18n("Jupiter"))
375     {
376         magnitude = -9.25 + param + 0.005 * phase;
377     }
378     else if (name() == i18n("Saturn"))
379     {
380         double T     = num->julianCenturies();
381         double a0    = (40.66 - 4.695 * T) * dms::PI / 180.;
382         double d0    = (83.52 + 0.403 * T) * dms::PI / 180.;
383         double sinx  = -cos(d0) * cosDec * cos(a0 - ra().radians());
384         sinx         = fabs(sinx - sin(d0) * sinDec);
385         double rings = -2.6 * sinx + 1.25 * sinx * sinx;
386         magnitude    = -8.88 + param + 0.044 * phase + rings;
387     }
388     else if (name() == i18n("Uranus"))
389     {
390         magnitude = -7.19 + param + 0.0028 * phase;
391     }
392     else if (name() == i18n("Neptune"))
393     {
394         magnitude = -6.87 + param;
395     }
396     setMag(magnitude);
397 }
398 
getUID() const399 SkyObject::UID KSPlanet::getUID() const
400 {
401     SkyObject::UID n;
402     if (name() == i18n("Mercury"))
403     {
404         n = 1;
405     }
406     else if (name() == i18n("Venus"))
407     {
408         n = 2;
409     }
410     else if (name() == i18n("Earth"))
411     {
412         n = 3;
413     }
414     else if (name() == i18n("Mars"))
415     {
416         n = 4;
417     }
418     else if (name() == i18n("Jupiter"))
419     {
420         n = 5;
421     }
422     else if (name() == i18n("Saturn"))
423     {
424         n = 6;
425     }
426     else if (name() == i18n("Uranus"))
427     {
428         n = 7;
429     }
430     else if (name() == i18n("Neptune"))
431     {
432         n = 8;
433     }
434     else
435     {
436         return SkyObject::invalidUID;
437     }
438     return solarsysUID(UID_SOL_BIGOBJ) | n;
439 }
440