Kstars

ksplanet.cpp
1 /*
2  SPDX-FileCopyrightText: 2001 Jason Harris <[email protected]>
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 
21 KSPlanet::OrbitDataManager KSPlanet::odm;
22 
24 {
25  //EMPTY
26 }
27 
28 bool 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 
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 
114 KSPlanet::KSPlanet(const QString &s, const QString &imfile, const QColor &c, double pSize)
115  : KSPlanetBase(s, imfile, c, pSize)
116 {
117 }
118 
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 
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
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
187 {
188  OrbitDataColl odc;
189  return odm.loadData(odc, untranslatedName());
190 }
191 
192 void 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 
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 
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 
346 void 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 
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 }
OrbitDataColl contains three groups of six QVectors.
Definition: ksplanet.h:121
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:380
double julianMillenia() const
Definition: ksnumbers.h:94
bool loadData() override
Preload the data used by findPosition.
Definition: ksplanet.cpp:186
virtual void setD(const double &x)
Sets floating-point value of angle, in degrees.
Definition: dms.h:179
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
OrbitDataManager places the OrbitDataColl objects for all planets in a QDict indexed by the planets' ...
Definition: ksplanet.h:139
void setDec0(dms d)
Sets Dec0, the catalog Declination.
Definition: skypoint.h:119
void SinCos(double &s, double &c) const
Get the sine and cosine together.
Definition: cachingdms.h:175
void append(const T &value)
OrbitDataManager()
Constructor.
Definition: ksplanet.cpp:23
void SinCos(double &s, double &c) const
Compute Sine and Cosine of the angle simultaneously.
Definition: dms.h:439
virtual QString name(void) const
Definition: skyobject.h:145
double rsun() const
Definition: ksplanetbase.h:130
double rearth() const
Definition: ksplanetbase.h:139
KSPlanet * clone() const override
Create copy of object.
Definition: ksplanet.cpp:150
void setRA0(dms r)
Sets RA0, the catalog Right Ascension.
Definition: skypoint.h:94
qint64 UID
Type for Unique object IDenticator.
Definition: skyobject.h:49
The ecliptic position of a planet (Longitude, Latitude, and distance from Sun).
Definition: ksplanetbase.h:25
static const UID invalidUID
Invalid UID.
Definition: skyobject.h:58
bool findGeocentricPosition(const KSNumbers *num, const KSPlanetBase *Earth=nullptr) override
Calculate the geocentric RA, Dec coordinates of the Planet.
Definition: ksplanet.cpp:265
QString i18n(const char *text, const TYPE &arg...)
QString & setNum(short n, int base)
Store several time-dependent astronomical quantities.
Definition: ksnumbers.h:42
const CachingDms & dec() const
Definition: skypoint.h:269
void setMag(float m)
Set the object's sorting magnitude.
Definition: skyobject.h:403
void setRearth(double r)
Set the distance from Earth, in AU.
Definition: ksplanetbase.h:145
SkipEmptyParts
const dms & ecLong() const
Definition: ksplanetbase.h:91
bool loadData(OrbitDataColl &odc, const QString &n)
Load orbital data for a planet from disk.
Definition: ksplanet.cpp:56
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
Provides necessary information about objects in the solar system.
Definition: ksplanet.h:32
UID solarsysUID(UID type) const
Compute high 32-bits of UID.
Definition: ksplanetbase.h:223
An angle, stored as degrees, but expressible in many ways.
Definition: dms.h:37
virtual void setRadians(const double &Rad)
Set angle according to the argument, in radians.
Definition: dms.h:328
QString toLower() const const
SkyObject::UID getUID() const override
Return UID for object.
Definition: ksplanet.cpp:402
const CachingDms & ra() const
Definition: skypoint.h:263
const CachingDms * obliquity() const
Definition: ksnumbers.h:56
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:700
const double & Degrees() const
Definition: dms.h:141
void setRsun(double r)
Set the solar distance in AU.
Definition: ksplanetbase.h:136
double julianCenturies() const
Definition: ksnumbers.h:88
const dms reduce() const
return the equivalent angle between 0 and 360 degrees.
Definition: dms.cpp:251
void EclipticToEquatorial(const CachingDms *Obliquity)
Convert Ecliptic longitude/latitude to Right Ascension/Declination.
const dms & ecLat() const
Definition: ksplanetbase.h:94
KSPlanet(const QString &s="unnamed", const QString &image_file=QString(), const QColor &c=Qt::white, double pSize=0)
Constructor.
Definition: ksplanet.cpp:114
static const UID UID_SOL_BIGOBJ
Big object.
Definition: ksplanetbase.h:216
void findPA(const KSNumbers *num)
Determine the position angle of the planet for a given date (used internally by findPosition() )
Provides necessary information about objects in the solar system.
Definition: ksplanetbase.h:49
This file is part of the KDE documentation.
Documentation copyright © 1996-2023 The KDE developers.
Generated on Mon May 8 2023 03:57:32 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.