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