Kstars

kscomet.cpp
1 /*
2  SPDX-FileCopyrightText: 2001 Jason Harris <[email protected]>
3 
4  SPDX-License-Identifier: GPL-2.0-or-later
5 */
6 
7 #include "kscomet.h"
8 
9 #include "ksnumbers.h"
10 #include "kstarsdata.h"
11 
12 #include <QDir>
13 
14 namespace
15 {
16 int letterToNum(QChar c)
17 {
18  char l = c.toLatin1();
19  if (l < 'A' || l > 'Z' || l == 'I')
20  return 0;
21  if (l > 'I')
22  return l - 'A';
23  return l - 'A' + 1;
24 }
25 // Convert letter designation like "AZ" to number.
26 int letterDesigToN(QString s)
27 {
28  int n = 0;
29  foreach (const QChar &c, s)
30  {
31  int nl = letterToNum(c);
32  if (nl == 0)
33  return 0;
34  n = n * 25 + nl;
35  }
36  return n;
37 }
38 
39 QMap<QChar, qint64> cometType;
40 }
41 
42 KSComet::KSComet(const QString &_s, const QString &imfile, double _q, double _e, dms _i, dms _w,
43  dms _Node, double Tp, float _M1, float _M2, float _K1, float _K2)
44  : KSPlanetBase(_s, imfile), q(_q), e(_e), M1(_M1), M2(_M2), K1(_K1), K2(_K2), i(_i), w(_w), N(_Node)
45 {
46  setType(SkyObject::COMET);
47 
48  //Find the Julian Day of Perihelion from Tp
49  //Tp is a double which encodes a date like: YYYYMMDD.DDDDD (e.g., 19730521.33333
50  // int year = int(Tp / 10000.0);
51  // int month = int((int(Tp) % 10000) / 100.0);
52  // int day = int(int(Tp) % 100);
53  // double Hour = 24.0 * (Tp - int(Tp));
54  // int h = int(Hour);
55  // int m = int(60.0 * (Hour - h));
56  // int s = int(60.0 * (60.0 * (Hour - h) - m));
57 
58  //JDp = KStarsDateTime(QDate(year, month, day), QTime(h, m, s)).djd();
59 
60  // JM 2021.10.21: In the new JPL format, it's already in JD
61  JDp = Tp;
62 
63  //compute the semi-major axis, a:
64  if (e == 1)
65  a = 0;
66  else
67  a = q / (1.0 - e);
68 
69 
70  //Compute the orbital Period from Kepler's 3rd law:
71  P = 365.2568984 * pow(a, 1.5); //period in days
72 
73  //If the name contains a "/", make this name2 and make name a truncated version without the leading "P/" or "C/"
74  // 2017-10-03 Jasem: We should keep the full name
75  /*if (name().contains(QDir::separator()))
76  {
77  setLongName(name());
78  setName(name().replace("P/", " ").trimmed());
79  setName(name().remove("C/"));
80  }*/
81 
82  // Try to calculate UID for comets. It's derived from comet designation.
83  // To parse name string regular expressions are used. Not really readable.
84  // And its make strong assumptions about format of comets' names.
85  // Probably come better algorithm should be used.
86 
87  // Periodic comet. Designation like: 12P or 12P-C
88  QRegExp rePer("^(\\d+)[PD](-([A-Z]+))?");
89  if (rePer.indexIn(_s) != -1)
90  {
91  // Comet number
92  qint64 num = rePer.cap(1).toInt();
93  // Fragment number
94  qint64 fragmentN = letterDesigToN(rePer.cap(2));
95  // Set UID
96  uidPart = num << 16 | fragmentN;
97  return;
98  }
99 
100  // Non periodic comet or comets with single approach
101  // Designations like C/(2006 A7)
102  QRegExp rePro("^([PCXDA])/.*\\((\\d{4}) ([A-Z])(\\d+)(-([A-Z]+))?\\)");
103  if (rePro.indexIn(_s) != -1)
104  {
105  // Fill comet type
106  if (cometType.empty())
107  {
108  cometType.insert('P', 0);
109  cometType.insert('C', 1);
110  cometType.insert('X', 2);
111  cometType.insert('D', 3);
112  cometType.insert('A', 4);
113  }
114  qint64 type = cometType[rePro.cap(1).at(0)]; // Type of comet
115  qint64 year = rePro.cap(2).toInt(); // Year of discovery
116  qint64 halfMonth = letterToNum(rePro.cap(3).at(0));
117  qint64 nHalfMonth = rePro.cap(4).toInt();
118  qint64 fragment = letterDesigToN(rePro.cap(6));
119 
120  uidPart = 1LL << 43 | type << 40 | // Bits 40-42 (3)
121  halfMonth << 33 | // Bits 33-39 (7) Hope this is enough
122  nHalfMonth << 28 | // Bits 28-32 (5)
123  year << 16 | // Bits 16-27 (12)
124  fragment; // Bits 0-15 (16)
125  return;
126  }
127  // qDebug() << Q_FUNC_INFO << "Didn't get it: " << _s;
128 }
129 
131 {
132  Q_ASSERT(typeid(this) == typeid(static_cast<const KSComet *>(this))); // Ensure we are not slicing a derived class
133  return new KSComet(*this);
134 }
135 
137 {
138  // Compute and store the estimated Physical size of the comet's coma, tail and nucleus
139  // References:
140  // * https://www.projectpluto.com/update7b.htm#comet_tail_formula [Project Pluto / GUIDE]
141  // * http://articles.adsabs.harvard.edu//full/1978BAICz..29..103K/0000113.000.html [Kresak, 1978a, "Passages of comets and asteroids near the earth"]
142  NuclearSize = pow(10, 2.1 - 0.2 * M1);
143  double mHelio = M1 + K1 * log10(rsun());
144  double L0, D0, L, D;
145  L0 = pow(10, -0.0075 * mHelio * mHelio - 0.19 * mHelio + 2.10);
146  D0 = pow(10, -0.0033 * mHelio * mHelio - 0.07 * mHelio + 3.25);
147  L = L0 * (1 - pow(10, -4 * rsun())) * (1 - pow(10, -2 * rsun()));
148  D = D0 * (1 - pow(10, -2 * rsun())) * (1 - pow(10, -rsun()));
149  TailSize = L * 1e6;
150  ComaSize = D * 1e3;
151  setPhysicalSize(ComaSize);
152 }
153 
155 {
156  // All elements are in the heliocentric ecliptic J2000 reference frame.
157 
158  double v(0.0), r(0.0);
159 
160  // Different between lastJD and Tp (Time of periapsis (Julian Day Number))
161  long double deltaJDP = lastPrecessJD - JDp;
162 
163  if (e > 0.98)
164  {
165  //Use near-parabolic approximation
166  double k = 0.01720209895; //Gauss gravitational constant
167  double a = 0.75 * (deltaJDP) * k * sqrt((1 + e) / (q * q * q));
168  double b = sqrt(1.0 + a * a);
169  double W = pow((b + a), 1.0 / 3.0) - pow((b - a), 1.0 / 3.0);
170  double c = 1.0 + 1.0 / (W * W);
171  double f = (1.0 - e) / (1.0 + e);
172  double g = f / (c * c);
173 
174  double a1 = (2.0 / 3.0) + (2.0 * W * W / 5.0);
175  double a2 = (7.0 / 5.0) + (33.0 * W * W / 35.0) + (37.0 * W * W * W * W / 175.0);
176  double a3 = W * W * ((432.0 / 175.0) + (956.0 * W * W / 1125.0) + (84.0 * W * W * W * W / 1575.0));
177  double w = W * (1.0 + g * c * (a1 + a2 * g + a3 * g * g));
178 
179  v = 2.0 * atan(w) / dms::DegToRad;
180  r = q * (1.0 + w * w) / (1.0 + w * w * f);
181  }
182  else
183  {
184  //Use normal ellipse method
185  //Determine Mean anomaly for desired date. deltaJDP is the difference between current JD minus JD of comet epoch.
186  // In JPL data, the Modified Julian Day is given to designate the epoch of comet data, which we convert to JD.
187  // Check http://astro.if.ufrgs.br/trigesf/position.html#17 for more details
188 
189  dms m = dms(double(360.0 * (deltaJDP) / P)).reduce();
190  double sinm, cosm;
191  m.SinCos(sinm, cosm);
192 
193  //compute eccentric anomaly:
194  double E = m.Degrees() + e * 180.0 / dms::PI * sinm * (1.0 + e * cosm);
195 
196  if (e > 0.005) //need more accurate approximation, iterate...
197  {
198  double E0;
199  int iter(0);
200  do
201  {
202  E0 = E;
203  iter++;
204  E = E0 - (E0 - e * 180.0 / dms::PI * sin(E0 * dms::DegToRad) - m.Degrees()) /
205  (1 - e * cos(E0 * dms::DegToRad));
206  }
207  while (fabs(E - E0) > 0.001 && iter < 1000);
208  }
209 
210  // Assert that the solution of the Kepler equation E = M + e sin E is accurate to about 0.1 arcsecond
211  //Q_ASSERT( fabs( E - ( m.Degrees() + ( e * 180.0 / dms::PI ) * sin( E * dms::DegToRad ) ) ) < 0.10/3600.0 );
212 
213  double sinE, cosE;
214  dms E1(E);
215  E1.SinCos(sinE, cosE);
216 
217  double xv = a * (cosE - e);
218  double yv = a * sqrt(1.0 - e * e) * sinE;
219 
220  //v is the true anomaly in degrees
221  v = atan2(yv, xv) / dms::DegToRad;
222  // Comet-Sun Heliocentric Distance in AU
223  r = sqrt(xv * xv + yv * yv);
224  }
225 
226  //Precess the longitude of the Ascending Node to the desired epoch
227  // i, w, and N are supplied in J2000 Epoch from JPL
228  // http://astro.if.ufrgs.br/trigesf/position.html#16
229  //dms n = dms(double(N.Degrees() + 3.82394E-5 * (lastPrecessJD - J2000))).reduce();
230 
231  //vw is the sum of the true anomaly and the argument of perihelion
232  dms vw(v + w.Degrees());
233 
234  double sinN, cosN, sinvw, cosvw, sini, cosi;
235 
236  // Get sin's and cos's for:
237  // Longitude of ascending node
238  N.SinCos(sinN, cosN);
239  // sum of true anomaly and argument of perihelion
240  vw.SinCos(sinvw, cosvw);
241  // Inclination
242  i.SinCos(sini, cosi);
243 
244  //xh, yh, zh are the heliocentric cartesian coords with the ecliptic plane congruent with zh=0.
245  double xh = r * (cosN * cosvw - sinN * sinvw * cosi);
246  double yh = r * (sinN * cosvw + cosN * sinvw * cosi);
247  double zh = r * (sinvw * sini);
248 
249  //the spherical ecliptic coordinates:
250  double ELongRad = atan2(yh, xh);
251  double ELatRad = atan2(zh, r);
252 
253  helEcPos.longitude.setRadians(ELongRad);
254  helEcPos.longitude.reduceToRange(dms::ZERO_TO_2PI);
255  helEcPos.latitude.setRadians(ELatRad);
256  setRsun(r);
257 
258  //xe, ye, ze are the Earth's heliocentric cartesian coords
259  double cosBe, sinBe, cosLe, sinLe;
260  Earth->ecLong().SinCos(sinLe, cosLe);
261  Earth->ecLat().SinCos(sinBe, cosBe);
262 
263  double xe = Earth->rsun() * cosBe * cosLe;
264  double ye = Earth->rsun() * cosBe * sinLe;
265  double ze = Earth->rsun() * sinBe;
266 
267  //convert to geocentric ecliptic coordinates by subtracting Earth's coords:
268  xh -= xe;
269  yh -= ye;
270  zh -= ze;
271 
272  //Finally, the spherical ecliptic coordinates:
273  ELongRad = atan2(yh, xh);
274  double rr = sqrt(xh * xh + yh * yh);
275  ELatRad = atan2(zh, rr);
276 
277  ep.longitude.setRadians(ELongRad);
278  ep.longitude.reduceToRange(dms::ZERO_TO_2PI);
279  ep.latitude.setRadians(ELatRad);
280  setRearth(Earth);
281 
282  // Now convert to geocentric equatorial coords
284 
285  // JM 2017-09-10: The calculations above produce J2000 RESULTS
286  // So we have to precess as well
287  setRA0(ra());
288  setDec0(dec());
289  apparentCoord(J2000, lastPrecessJD);
291 
292  return true;
293 }
294 
295 // m = M1 + 2.5 * K1 * log10(rsun) + 5 * log10(rearth)
296 void KSComet::findMagnitude(const KSNumbers *)
297 {
298  setMag(M1 + 2.5 * K1 * log10(rsun()) + 5 * log10(rearth()));
299 }
300 
301 void KSComet::setEarthMOID(double earth_moid)
302 {
303  EarthMOID = earth_moid;
304 }
305 
306 void KSComet::setAlbedo(float albedo)
307 {
308  Albedo = albedo;
309 }
310 
311 void KSComet::setDiameter(float diam)
312 {
313  Diameter = diam;
314 }
315 
317 {
318  Dimensions = dim;
319 }
320 
321 void KSComet::setNEO(bool neo)
322 {
323  NEO = neo;
324 }
325 
327 {
328  OrbitClass = orbit_class;
329 }
330 
332 {
333  OrbitID = orbit_id;
334 }
335 
336 void KSComet::setPeriod(float per)
337 {
338  Period = per;
339 }
340 
341 void KSComet::setRotationPeriod(float rot_per)
342 {
343  RotationPeriod = rot_per;
344 }
345 
346 //Unused virtual function from KSPlanetBase
348 {
349  return false;
350 }
351 
353 {
354  return solarsysUID(UID_SOL_COMET) | uidPart;
355 }
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
void setDimensions(QString dim)
Sets the comet's dimensions.
Definition: kscomet.cpp:316
void setDec0(dms d)
Sets Dec0, the catalog Declination.
Definition: skypoint.h:119
static constexpr double DegToRad
DegToRad is a const static member equal to the number of radians in one degree (dms::PI/180....
Definition: dms.h:390
void SinCos(double &s, double &c) const
Compute Sine and Cosine of the angle simultaneously.
Definition: dms.h:444
bool findGeocentricPosition(const KSNumbers *num, const KSPlanetBase *Earth=nullptr) override
Calculate the geocentric RA, Dec coordinates of the Comet.
Definition: kscomet.cpp:154
A subclass of KSPlanetBase that implements comets.
Definition: kscomet.h:43
double rsun() const
Definition: ksplanetbase.h:130
void setDiameter(float diam)
Sets the comet's diameter.
Definition: kscomet.cpp:311
double rearth() const
Definition: ksplanetbase.h:139
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
int type(void) const
Definition: skyobject.h:188
bool loadData() override
Unused virtual function inherited from KSPlanetBase thus it's simply empty here.
Definition: kscomet.cpp:347
KSComet(const QString &s, const QString &image_file, double q, double e, dms i, dms w, dms N, double Tp, float M1, float M2, float K1, float K2)
Constructor.
Definition: kscomet.cpp:42
void setPeriod(float per)
Sets the comet's period.
Definition: kscomet.cpp:336
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
const dms & ecLong() const
Definition: ksplanetbase.h:91
int indexIn(const QString &str, int offset, QRegExp::CaretMode caretMode) const const
SkyObject::UID getUID() const override
Return UID for object.
Definition: kscomet.cpp:352
int toInt(bool *ok, int base) const const
void setEarthMOID(double earth_moid)
Sets the comet's earth minimum orbit intersection distance.
Definition: kscomet.cpp:301
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
void findPhysicalParameters()
Estimate physical parameters of the comet such as coma size, tail length and size of the nucleus.
Definition: kscomet.cpp:136
virtual void setRadians(const double &Rad)
Set angle according to the argument, in radians.
Definition: dms.h:333
void setOrbitClass(QString orbit_class)
Sets the comet's orbit class.
Definition: kscomet.cpp:326
QString cap(int nth) const const
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
const QChar at(int position) const const
void setOrbitID(QString orbit_id)
Sets the comet's orbit solution ID.
Definition: kscomet.cpp:331
void setNEO(bool neo)
Sets if the comet is a near earth object.
Definition: kscomet.cpp:321
KSComet * clone() const override
Create copy of object.
Definition: kscomet.cpp:130
void reduceToRange(enum dms::AngleRanges range)
Reduce this angle to the given range.
Definition: dms.cpp:438
const dms reduce() const
return the equivalent angle between 0 and 360 degrees.
Definition: dms.cpp:251
void setPhysicalSize(double size)
set the planet's physical size, in km.
Definition: ksplanetbase.h:197
void setRotationPeriod(float rot_per)
Sets the comet's rotation period.
Definition: kscomet.cpp:341
void EclipticToEquatorial(const CachingDms *Obliquity)
Convert Ecliptic longitude/latitude to Right Ascension/Declination.
char toLatin1() const const
static const UID UID_SOL_COMET
Comets.
Definition: ksplanetbase.h:220
void setAlbedo(float albedo)
Sets the comet's albedo.
Definition: kscomet.cpp:306
const dms & ecLat() const
Definition: ksplanetbase.h:94
void setType(int t)
Set the object's type identifier to the argument.
Definition: skyobject.h:195
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 Fri Jun 9 2023 04:02:22 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.