Kstars

ksplanet.cpp
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// Qt version calming
19#include <qtskipemptyparts.h>
20
21KSPlanet::OrbitDataManager KSPlanet::odm;
22
27
28bool KSPlanet::OrbitDataManager::readOrbitData(const QString &fname, QVector<OrbitData> *vector)
29{
30 QFile f;
31
32 if (KSUtils::openDataFile(f, fname))
33 {
34 KSFileReader fileReader(f); // close file is included
35 QStringList fields;
36 while (fileReader.hasMoreLines())
37 {
38 fields = fileReader.readLine().split(' ', Qt::SkipEmptyParts);
39
40 if (fields.size() == 3)
41 {
42 double A = fields[0].toDouble();
43 double B = fields[1].toDouble();
44 double C = fields[2].toDouble();
45 vector->append(OrbitData(A, B, C));
46 }
47 }
48 }
49 else
50 {
51 return false;
52 }
53 return true;
54}
55
56bool KSPlanet::OrbitDataManager::loadData(KSPlanet::OrbitDataColl &odc, const QString &n)
57{
58 QString fname, snum;
59 QFile f;
60 int nCount = 0;
61 QString nl = n.toLower();
62
63 if (hash.contains(nl))
64 {
65 odc = hash[nl];
66 return true; //orbit data already loaded
67 }
68
69 //Create a new OrbitDataColl
70 OrbitDataColl ret;
71
72 //Ecliptic Longitude
73 for (int i = 0; i < 6; ++i)
74 {
75 snum.setNum(i);
76 fname = nl + ".L" + snum + ".vsop";
77 if (readOrbitData(fname, &ret.Lon[i]))
78 nCount++;
79 }
80
81 if (nCount == 0)
82 return false;
83
84 //Ecliptic Latitude
85 for (int i = 0; i < 6; ++i)
86 {
87 snum.setNum(i);
88 fname = nl + ".B" + snum + ".vsop";
89 if (readOrbitData(fname, &ret.Lat[i]))
90 nCount++;
91 }
92
93 if (nCount == 0)
94 return false;
95
96 //Heliocentric Distance
97 for (int i = 0; i < 6; ++i)
98 {
99 snum.setNum(i);
100 fname = nl + ".R" + snum + ".vsop";
101 if (readOrbitData(fname, &ret.Dst[i]))
102 nCount++;
103 }
104
105 if (nCount == 0)
106 return false;
107
108 hash[nl] = ret;
109 odc = hash[nl];
110
111 return true;
112}
113
114KSPlanet::KSPlanet(const QString &s, const QString &imfile, const QColor &c, double pSize)
115 : KSPlanetBase(s, imfile, c, pSize)
116{
117}
118
119KSPlanet::KSPlanet(int n) : KSPlanetBase()
120{
121 switch (n)
122 {
123 case MERCURY:
124 KSPlanetBase::init(i18n("Mercury"), "mercury", KSPlanetBase::planetColor[KSPlanetBase::MERCURY], 4879.4);
125 break;
126 case VENUS:
127 KSPlanetBase::init(i18n("Venus"), "venus", KSPlanetBase::planetColor[KSPlanetBase::VENUS], 12103.6);
128 break;
129 case MARS:
130 KSPlanetBase::init(i18n("Mars"), "mars", KSPlanetBase::planetColor[KSPlanetBase::MARS], 6792.4);
131 break;
132 case JUPITER:
133 KSPlanetBase::init(i18n("Jupiter"), "jupiter", KSPlanetBase::planetColor[KSPlanetBase::JUPITER], 142984.);
134 break;
135 case SATURN:
136 KSPlanetBase::init(i18n("Saturn"), "saturn", KSPlanetBase::planetColor[KSPlanetBase::SATURN], 120536.);
137 break;
138 case URANUS:
139 KSPlanetBase::init(i18n("Uranus"), "uranus", KSPlanetBase::planetColor[KSPlanetBase::URANUS], 51118.);
140 break;
141 case NEPTUNE:
142 KSPlanetBase::init(i18n("Neptune"), "neptune", KSPlanetBase::planetColor[KSPlanetBase::NEPTUNE], 49572.);
143 break;
144 default:
145 qDebug() << Q_FUNC_INFO << "Error: Illegal identifier in KSPlanet constructor: " << n;
146 break;
147 }
148}
149
150KSPlanet *KSPlanet::clone() const
151{
152 Q_ASSERT(typeid(this) == typeid(static_cast<const KSPlanet *>(this))); // Ensure we are not slicing a derived class
153 return new KSPlanet(*this);
154}
155
156// TODO: Get rid of this dirty hack post KDE 4.2 release
157QString KSPlanet::untranslatedName() const
158{
159 if (name() == i18n("Mercury"))
160 return "Mercury";
161 else if (name() == i18n("Venus"))
162 return "Venus";
163 else if (name() == i18n("Mars"))
164 return "Mars";
165 else if (name() == i18n("Jupiter"))
166 return "Jupiter";
167 else if (name() == i18n("Saturn"))
168 return "Saturn";
169 else if (name() == i18n("Uranus"))
170 return "Uranus";
171 else if (name() == i18n("Neptune"))
172 return "Neptune";
173 else if (name() == i18n("Earth"))
174 return "Earth";
175 else if (name() == i18n("Earth Shadow"))
176 return "Earth Shadow";
177 else if (name() == i18n("Moon"))
178 return "Moon";
179 else if (name() == i18n("Sun"))
180 return "Sun";
181 else
182 return name();
183}
184
185//we don't need the reference to the ODC, so just give it a junk variable
186bool KSPlanet::loadData()
187{
188 OrbitDataColl odc;
189 return odm.loadData(odc, untranslatedName());
190}
191
192void KSPlanet::calcEcliptic(double Tau, EclipticPosition &epret) const
193{
194 double sum[6];
195 OrbitDataColl odc;
196 double Tpow[6];
197
198 Tpow[0] = 1.0;
199 for (int i = 1; i < 6; ++i)
200 {
201 Tpow[i] = Tpow[i - 1] * Tau;
202 }
203
204 if (!odm.loadData(odc, untranslatedName()))
205 {
206 epret.longitude = dms(0.0);
207 epret.latitude = dms(0.0);
208 epret.radius = 0.0;
209 qCWarning(KSTARS) << "Could not get data for name:" << name() << "(" << untranslatedName() << ")";
210 return;
211 }
212
213 //Ecliptic Longitude
214 for (int i = 0; i < 6; ++i)
215 {
216 sum[i] = 0.0;
217 for (int j = 0; j < odc.Lon[i].size(); ++j)
218 {
219 sum[i] += odc.Lon[i][j].A * cos(odc.Lon[i][j].B + odc.Lon[i][j].C * Tau);
220 /*
221 qDebug() << Q_FUNC_INFO << "sum[" << i <<"] =" << sum[i] <<
222 " A = " << odc.Lon[i][j].A << " B = " << odc.Lon[i][j].B <<
223 " C = " << odc.Lon[i][j].C;
224 */
225 }
226 sum[i] *= Tpow[i];
227 //qDebug() << Q_FUNC_INFO << name() << " : sum[" << i << "] = " << sum[i];
228 }
229
230 epret.longitude.setRadians(sum[0] + sum[1] + sum[2] + sum[3] + sum[4] + sum[5]);
231 epret.longitude.setD(epret.longitude.reduce().Degrees());
232
233 //Compute Ecliptic Latitude
234 for (uint i = 0; i < 6; ++i)
235 {
236 sum[i] = 0.0;
237 for (int j = 0; j < odc.Lat[i].size(); ++j)
238 {
239 sum[i] += odc.Lat[i][j].A * cos(odc.Lat[i][j].B + odc.Lat[i][j].C * Tau);
240 }
241 sum[i] *= Tpow[i];
242 }
243
244 epret.latitude.setRadians(sum[0] + sum[1] + sum[2] + sum[3] + sum[4] + sum[5]);
245
246 //Compute Heliocentric Distance
247 for (uint i = 0; i < 6; ++i)
248 {
249 sum[i] = 0.0;
250 for (int j = 0; j < odc.Dst[i].size(); ++j)
251 {
252 sum[i] += odc.Dst[i][j].A * cos(odc.Dst[i][j].B + odc.Dst[i][j].C * Tau);
253 }
254 sum[i] *= Tpow[i];
255 }
256
257 epret.radius = sum[0] + sum[1] + sum[2] + sum[3] + sum[4] + sum[5];
258
259 /*
260 qDebug() << Q_FUNC_INFO << name() << " pre: Lat = " << epret.latitude.toDMSString() << " Long = " <<
261 epret.longitude.toDMSString() << " Dist = " << epret.radius;
262 */
263}
264
265bool KSPlanet::findGeocentricPosition(const KSNumbers *num, const KSPlanetBase *Earth)
266{
267 if (Earth != nullptr)
268 {
269 double sinL, sinL0, sinB, sinB0;
270 double cosL, cosL0, cosB, cosB0;
271 double x = 0.0, y = 0.0, z = 0.0;
272
273 double olddst = -1000;
274 double dst = 0;
275
276 EclipticPosition trialpos;
277
278 double jm = num->julianMillenia();
279
280 Earth->ecLong().SinCos(sinL0, cosL0);
281 Earth->ecLat().SinCos(sinB0, cosB0);
282
283 double eX = Earth->rsun() * cosB0 * cosL0;
284 double eY = Earth->rsun() * cosB0 * sinL0;
285 double eZ = Earth->rsun() * sinB0;
286
287 bool once = true;
288 while (fabs(dst - olddst) > .001)
289 {
290 calcEcliptic(jm, trialpos);
291
292 // We store the heliocentric ecliptic coordinates the first time they are computed.
293 if (once)
294 {
295 helEcPos = trialpos;
296 once = false;
297 }
298
299 olddst = dst;
300
301 trialpos.longitude.SinCos(sinL, cosL);
302 trialpos.latitude.SinCos(sinB, cosB);
303
304 x = trialpos.radius * cosB * cosL - eX;
305 y = trialpos.radius * cosB * sinL - eY;
306 z = trialpos.radius * sinB - eZ;
307
308 //distance from Earth
309 dst = sqrt(x * x + y * y + z * z);
310
311 //The light-travel time delay, in millenia
312 //0.0057755183 is the inverse speed of light,
313 //in days/AU
314 double delay = (.0057755183 * dst) / 365250.0;
315
316 jm = num->julianMillenia() - delay;
317 }
318
319 ep.longitude.setRadians(atan2(y, x));
320 ep.longitude.reduce();
321 ep.latitude.setRadians(atan2(z, sqrt(x * x + y * y)));
322 setRsun(trialpos.radius);
323 setRearth(dst);
324
325 EclipticToEquatorial(num->obliquity());
326
327 setRA0(ra());
328 setDec0(dec());
329 apparentCoord(J2000, lastPrecessJD);
330
331 //nutate(num);
332 //aberrate(num);
333 }
334 else
335 {
336 calcEcliptic(num->julianMillenia(), ep);
337 helEcPos = ep;
338 }
339
340 //determine the position angle
341 findPA(num);
342
343 return true;
344}
345
346void KSPlanet::findMagnitude(const KSNumbers *num)
347{
348 double cosDec, sinDec;
349 dec().SinCos(cosDec, sinDec);
350
351 /* Computation of the visual magnitude (V band) of the planet.
352 * Algorithm provided by Pere Planesas (Observatorio Astronomico Nacional)
353 * It has some similarity to J. Meeus algorithm in Astronomical Algorithms, Chapter 40.
354 * */
355
356 // Initialized to the faintest magnitude observable with the HST
357 float magnitude = 30;
358
359 double param = 5 * log10(rsun() * rearth());
360 double phase = this->phase().Degrees();
361 double f1 = phase / 100.;
362
363 if (name() == i18n("Mercury"))
364 {
365 if (phase > 150.)
366 f1 = 1.5;
367 magnitude = -0.36 + param + 3.8 * f1 - 2.73 * f1 * f1 + 2 * f1 * f1 * f1;
368 }
369 else if (name() == i18n("Venus"))
370 {
371 magnitude = -4.29 + param + 0.09 * f1 + 2.39 * f1 * f1 - 0.65 * f1 * f1 * f1;
372 }
373 else if (name() == i18n("Mars"))
374 {
375 magnitude = -1.52 + param + 0.016 * phase;
376 }
377 else if (name() == i18n("Jupiter"))
378 {
379 magnitude = -9.25 + param + 0.005 * phase;
380 }
381 else if (name() == i18n("Saturn"))
382 {
383 double T = num->julianCenturies();
384 double a0 = (40.66 - 4.695 * T) * dms::PI / 180.;
385 double d0 = (83.52 + 0.403 * T) * dms::PI / 180.;
386 double sinx = -cos(d0) * cosDec * cos(a0 - ra().radians());
387 sinx = fabs(sinx - sin(d0) * sinDec);
388 double rings = -2.6 * sinx + 1.25 * sinx * sinx;
389 magnitude = -8.88 + param + 0.044 * phase + rings;
390 }
391 else if (name() == i18n("Uranus"))
392 {
393 magnitude = -7.19 + param + 0.0028 * phase;
394 }
395 else if (name() == i18n("Neptune"))
396 {
397 magnitude = -6.87 + param;
398 }
399 setMag(magnitude);
400}
401
402SkyObject::UID KSPlanet::getUID() const
403{
404 SkyObject::UID n;
405 if (name() == i18n("Mercury"))
406 {
407 n = 1;
408 }
409 else if (name() == i18n("Venus"))
410 {
411 n = 2;
412 }
413 else if (name() == i18n("Earth"))
414 {
415 n = 3;
416 }
417 else if (name() == i18n("Mars"))
418 {
419 n = 4;
420 }
421 else if (name() == i18n("Jupiter"))
422 {
423 n = 5;
424 }
425 else if (name() == i18n("Saturn"))
426 {
427 n = 6;
428 }
429 else if (name() == i18n("Uranus"))
430 {
431 n = 7;
432 }
433 else if (name() == i18n("Neptune"))
434 {
435 n = 8;
436 }
437 else
438 {
439 return SkyObject::invalidUID;
440 }
441 return solarsysUID(UID_SOL_BIGOBJ) | n;
442}
CachingDms::SinCos
void SinCos(double &s, double &c) const
Get the sine and cosine together.
Definition cachingdms.h:175
EclipticPosition
The ecliptic position of a planet (Longitude, Latitude, and distance from Sun).
Definition ksplanetbase.h:26
KSFileReader
I totally rewrote this because the earlier scheme of reading all the lines of a file into a buffer be...
Definition ksfilereader.h:55
KSNumbers
There are several time-dependent values used in position calculations, that are not specific to an ob...
Definition ksnumbers.h:43
KSNumbers::obliquity
const CachingDms * obliquity() const
Definition ksnumbers.h:56
KSNumbers::julianMillenia
double julianMillenia() const
Definition ksnumbers.h:94
KSNumbers::julianCenturies
double julianCenturies() const
Definition ksnumbers.h:88
KSPlanetBase::solarsysUID
UID solarsysUID(UID type) const
Compute high 32-bits of UID.
Definition ksplanetbase.h:223
KSPlanetBase::UID_SOL_BIGOBJ
static const UID UID_SOL_BIGOBJ
Big object.
Definition ksplanetbase.h:216
KSPlanetBase::setRearth
void setRearth(double r)
Set the distance from Earth, in AU.
Definition ksplanetbase.h:145
KSPlanetBase::KSPlanetBase
KSPlanetBase(const QString &s=i18n("unnamed"), const QString &image_file=QString(), const QColor &c=Qt::white, double pSize=0)
Constructor.
Definition ksplanetbase.cpp:37
KSPlanetBase::findPA
void findPA(const KSNumbers *num)
Determine the position angle of the planet for a given date (used internally by findPosition() )
Definition ksplanetbase.cpp:240
KSPlanetBase::ecLat
const dms & ecLat() const
Definition ksplanetbase.h:94
KSPlanetBase::ecLong
const dms & ecLong() const
Definition ksplanetbase.h:91
KSPlanetBase::rearth
double rearth() const
Definition ksplanetbase.h:139
KSPlanetBase::EclipticToEquatorial
void EclipticToEquatorial(const CachingDms *Obliquity)
Convert Ecliptic longitude/latitude to Right Ascension/Declination.
Definition ksplanetbase.cpp:81
KSPlanetBase::rsun
double rsun() const
Definition ksplanetbase.h:130
KSPlanetBase::setRsun
void setRsun(double r)
Set the solar distance in AU.
Definition ksplanetbase.h:136
KSPlanet::OrbitDataColl
OrbitDataColl contains three groups of six QVectors.
Definition ksplanet.h:122
KSPlanet::OrbitDataManager
OrbitDataManager places the OrbitDataColl objects for all planets in a QDict indexed by the planets' ...
Definition ksplanet.h:140
KSPlanet::OrbitDataManager::loadData
bool loadData(OrbitDataColl &odc, const QString &n)
Load orbital data for a planet from disk.
Definition ksplanet.cpp:56
KSPlanet::OrbitDataManager::OrbitDataManager
OrbitDataManager()
Constructor.
Definition ksplanet.cpp:23
KSPlanet::OrbitData
This class contains doubles A,B,C which represent a single term in a planet's positional expansion su...
Definition ksplanet.h:95
KSPlanet::clone
KSPlanet * clone() const override
Create copy of object.
Definition ksplanet.cpp:150
KSPlanet::getUID
SkyObject::UID getUID() const override
Return UID for object.
Definition ksplanet.cpp:402
KSPlanet::findGeocentricPosition
bool findGeocentricPosition(const KSNumbers *num, const KSPlanetBase *Earth=nullptr) override
Calculate the geocentric RA, Dec coordinates of the Planet.
Definition ksplanet.cpp:265
KSPlanet::KSPlanet
KSPlanet(const QString &s="unnamed", const QString &image_file=QString(), const QColor &c=Qt::white, double pSize=0)
Constructor.
Definition ksplanet.cpp:114
KSPlanet::loadData
bool loadData() override
Preload the data used by findPosition.
Definition ksplanet.cpp:186
KSPlanet::untranslatedName
QString untranslatedName() const
return the untranslated name This is a dirty way to solve a lot of localization-related trouble for t...
Definition ksplanet.cpp:157
KSPlanet::calcEcliptic
virtual void calcEcliptic(double jm, EclipticPosition &ret) const
Calculate the ecliptic longitude and latitude of the planet for the given date (expressed in Julian M...
Definition ksplanet.cpp:192
SkyObject::invalidUID
static const UID invalidUID
Invalid UID.
Definition skyobject.h:66
SkyObject::setMag
void setMag(float m)
Set the object's sorting magnitude.
Definition skyobject.h:472
SkyObject::UID
qint64 UID
Type for Unique object IDenticator.
Definition skyobject.h:57
SkyObject::name
virtual QString name(void) const
Definition skyobject.h:154
SkyPoint::apparentCoord
void apparentCoord(long double jd0, long double jdf)
Computes the apparent coordinates for this SkyPoint for any epoch, accounting for the effects of prec...
Definition skypoint.cpp:720
SkyPoint::dec
const CachingDms & dec() const
Definition skypoint.h:269
SkyPoint::ra
const CachingDms & ra() const
Definition skypoint.h:263
SkyPoint::setRA0
void setRA0(dms r)
Sets RA0, the catalog Right Ascension.
Definition skypoint.h:94
SkyPoint::setDec0
void setDec0(dms d)
Sets Dec0, the catalog Declination.
Definition skypoint.h:119
dms
An angle, stored as degrees, but expressible in many ways.
Definition dms.h:38
dms::reduce
const dms reduce() const
return the equivalent angle between 0 and 360 degrees.
Definition dms.cpp:251
dms::SinCos
void SinCos(double &s, double &c) const
Compute Sine and Cosine of the angle simultaneously.
Definition dms.h:447
dms::PI
static constexpr double PI
PI is a const static member; it's public so that it can be used anywhere, as long as dms....
Definition dms.h:385
dms::setRadians
virtual void setRadians(const double &Rad)
Set angle according to the argument, in radians.
Definition dms.h:333
dms::setD
virtual void setD(const double &x)
Sets floating-point value of angle, in degrees.
Definition dms.h:179
dms::Degrees
const double & Degrees() const
Definition dms.h:141
i18n
QString i18n(const char *text, const TYPE &arg...)
QList::append
void append(QList< T > &&value)
QString::setNum
QString & setNum(double n, char format, int precision)
QString::toLower
QString toLower() const const
Qt::SkipEmptyParts
SkipEmptyParts
This file is part of the KDE documentation.
Documentation copyright © 1996-2025 The KDE developers.
Generated on Fri May 2 2025 12:02:39 by doxygen 1.13.2 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.