Kstars

kscomet.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 "kscomet.h"
8
9#include "ksnumbers.h"
10#include "kstarsdata.h"
11
12#include <QDir>
13
14namespace
15{
16int 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.
26int 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
39QMap<QChar, qint64> cometType;
40}
41
42KSComet::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 QRegularExpression rePerReg("^(\\d+)[PD](-([A-Z]+))?");
89 QRegularExpressionMatch rePer = rePerReg.match(_s);
90 if (rePer.hasMatch())
91 {
92 // Comet number
93 qint64 num = rePer.captured(1).toInt();
94 // Fragment number
95 qint64 fragmentN = letterDesigToN(rePer.captured(2));
96 // Set UID
97 uidPart = num << 16 | fragmentN;
98 return;
99 }
100
101 // Non periodic comet or comets with single approach
102 // Designations like C/(2006 A7)
103 QRegularExpression reProReg("^([PCXDA])/.*\\((\\d{4}) ([A-Z])(\\d+)(-([A-Z]+))?\\)");
104 QRegularExpressionMatch rePro = reProReg.match(_s);
105 if (rePro.hasMatch())
106 {
107 // Fill comet type
108 if (cometType.empty())
109 {
110 cometType.insert('P', 0);
111 cometType.insert('C', 1);
112 cometType.insert('X', 2);
113 cometType.insert('D', 3);
114 cometType.insert('A', 4);
115 }
116 qint64 type = cometType[rePro.captured(1).at(0)]; // Type of comet
117 qint64 year = rePro.captured(2).toInt(); // Year of discovery
118 qint64 halfMonth = letterToNum(rePro.captured(3).at(0));
119 qint64 nHalfMonth = rePro.captured(4).toInt();
120 qint64 fragment = letterDesigToN(rePro.captured(6));
121
122 uidPart = 1LL << 43 | type << 40 | // Bits 40-42 (3)
123 halfMonth << 33 | // Bits 33-39 (7) Hope this is enough
124 nHalfMonth << 28 | // Bits 28-32 (5)
125 year << 16 | // Bits 16-27 (12)
126 fragment; // Bits 0-15 (16)
127 return;
128 }
129 // qDebug() << Q_FUNC_INFO << "Didn't get it: " << _s;
130}
131
133{
134 Q_ASSERT(typeid(this) == typeid(static_cast<const KSComet *>(this))); // Ensure we are not slicing a derived class
135 return new KSComet(*this);
136}
137
139{
140 // Compute and store the estimated Physical size of the comet's coma, tail and nucleus
141 // References:
142 // * https://www.projectpluto.com/update7b.htm#comet_tail_formula [Project Pluto / GUIDE]
143 // * http://articles.adsabs.harvard.edu//full/1978BAICz..29..103K/0000113.000.html [Kresak, 1978a, "Passages of comets and asteroids near the earth"]
144 NuclearSize = pow(10, 2.1 - 0.2 * M1);
145 double mHelio = M1 + K1 * log10(rsun());
146 double L0, D0, L, D;
147 L0 = pow(10, -0.0075 * mHelio * mHelio - 0.19 * mHelio + 2.10);
148 D0 = pow(10, -0.0033 * mHelio * mHelio - 0.07 * mHelio + 3.25);
149 L = L0 * (1 - pow(10, -4 * rsun())) * (1 - pow(10, -2 * rsun()));
150 D = D0 * (1 - pow(10, -2 * rsun())) * (1 - pow(10, -rsun()));
151 TailSize = L * 1e6;
152 ComaSize = D * 1e3;
153 setPhysicalSize(ComaSize);
154}
155
157{
158 // All elements are in the heliocentric ecliptic J2000 reference frame.
159
160 double v(0.0), r(0.0);
161
162 // Different between lastJD and Tp (Time of periapsis (Julian Day Number))
163 long double deltaJDP = lastPrecessJD - JDp;
164
165 if (e > 0.98)
166 {
167 //Use near-parabolic approximation
168 double k = 0.01720209895; //Gauss gravitational constant
169 double a = 0.75 * (deltaJDP) * k * sqrt((1 + e) / (q * q * q));
170 double b = sqrt(1.0 + a * a);
171 double W = pow((b + a), 1.0 / 3.0) - pow((b - a), 1.0 / 3.0);
172 double c = 1.0 + 1.0 / (W * W);
173 double f = (1.0 - e) / (1.0 + e);
174 double g = f / (c * c);
175
176 double a1 = (2.0 / 3.0) + (2.0 * W * W / 5.0);
177 double a2 = (7.0 / 5.0) + (33.0 * W * W / 35.0) + (37.0 * W * W * W * W / 175.0);
178 double a3 = W * W * ((432.0 / 175.0) + (956.0 * W * W / 1125.0) + (84.0 * W * W * W * W / 1575.0));
179 double w = W * (1.0 + g * c * (a1 + a2 * g + a3 * g * g));
180
181 v = 2.0 * atan(w) / dms::DegToRad;
182 r = q * (1.0 + w * w) / (1.0 + w * w * f);
183 }
184 else
185 {
186 //Use normal ellipse method
187 //Determine Mean anomaly for desired date. deltaJDP is the difference between current JD minus JD of comet epoch.
188 // In JPL data, the Modified Julian Day is given to designate the epoch of comet data, which we convert to JD.
189 // Check http://astro.if.ufrgs.br/trigesf/position.html#17 for more details
190
191 dms m = dms(double(360.0 * (deltaJDP) / P)).reduce();
192 double sinm, cosm;
193 m.SinCos(sinm, cosm);
194
195 //compute eccentric anomaly:
196 double E = m.Degrees() + e * 180.0 / dms::PI * sinm * (1.0 + e * cosm);
197
198 if (e > 0.005) //need more accurate approximation, iterate...
199 {
200 double E0;
201 int iter(0);
202 do
203 {
204 E0 = E;
205 iter++;
206 E = E0 - (E0 - e * 180.0 / dms::PI * sin(E0 * dms::DegToRad) - m.Degrees()) /
207 (1 - e * cos(E0 * dms::DegToRad));
208 }
209 while (fabs(E - E0) > 0.001 && iter < 1000);
210 }
211
212 // Assert that the solution of the Kepler equation E = M + e sin E is accurate to about 0.1 arcsecond
213 //Q_ASSERT( fabs( E - ( m.Degrees() + ( e * 180.0 / dms::PI ) * sin( E * dms::DegToRad ) ) ) < 0.10/3600.0 );
214
215 double sinE, cosE;
216 dms E1(E);
217 E1.SinCos(sinE, cosE);
218
219 double xv = a * (cosE - e);
220 double yv = a * sqrt(1.0 - e * e) * sinE;
221
222 //v is the true anomaly in degrees
223 v = atan2(yv, xv) / dms::DegToRad;
224 // Comet-Sun Heliocentric Distance in AU
225 r = sqrt(xv * xv + yv * yv);
226 }
227
228 //Precess the longitude of the Ascending Node to the desired epoch
229 // i, w, and N are supplied in J2000 Epoch from JPL
230 // http://astro.if.ufrgs.br/trigesf/position.html#16
231 //dms n = dms(double(N.Degrees() + 3.82394E-5 * (lastPrecessJD - J2000))).reduce();
232
233 //vw is the sum of the true anomaly and the argument of perihelion
234 dms vw(v + w.Degrees());
235
236 double sinN, cosN, sinvw, cosvw, sini, cosi;
237
238 // Get sin's and cos's for:
239 // Longitude of ascending node
240 N.SinCos(sinN, cosN);
241 // sum of true anomaly and argument of perihelion
242 vw.SinCos(sinvw, cosvw);
243 // Inclination
244 i.SinCos(sini, cosi);
245
246 //xh, yh, zh are the heliocentric cartesian coords with the ecliptic plane congruent with zh=0.
247 double xh = r * (cosN * cosvw - sinN * sinvw * cosi);
248 double yh = r * (sinN * cosvw + cosN * sinvw * cosi);
249 double zh = r * (sinvw * sini);
250
251 //the spherical ecliptic coordinates:
252 double ELongRad = atan2(yh, xh);
253 double ELatRad = atan2(zh, r);
254
255 helEcPos.longitude.setRadians(ELongRad);
256 helEcPos.longitude.reduceToRange(dms::ZERO_TO_2PI);
257 helEcPos.latitude.setRadians(ELatRad);
258 setRsun(r);
259
260 //xe, ye, ze are the Earth's heliocentric cartesian coords
261 double cosBe, sinBe, cosLe, sinLe;
262 Earth->ecLong().SinCos(sinLe, cosLe);
263 Earth->ecLat().SinCos(sinBe, cosBe);
264
265 double xe = Earth->rsun() * cosBe * cosLe;
266 double ye = Earth->rsun() * cosBe * sinLe;
267 double ze = Earth->rsun() * sinBe;
268
269 //convert to geocentric ecliptic coordinates by subtracting Earth's coords:
270 xh -= xe;
271 yh -= ye;
272 zh -= ze;
273
274 //Finally, the spherical ecliptic coordinates:
275 ELongRad = atan2(yh, xh);
276 double rr = sqrt(xh * xh + yh * yh);
277 ELatRad = atan2(zh, rr);
278
279 ep.longitude.setRadians(ELongRad);
280 ep.longitude.reduceToRange(dms::ZERO_TO_2PI);
281 ep.latitude.setRadians(ELatRad);
282 setRearth(Earth);
283
284 // Now convert to geocentric equatorial coords
286
287 // JM 2017-09-10: The calculations above produce J2000 RESULTS
288 // So we have to precess as well
289 setRA0(ra());
290 setDec0(dec());
291 apparentCoord(J2000, lastPrecessJD);
293
294 return true;
295}
296
297// m = M1 + 2.5 * K1 * log10(rsun) + 5 * log10(rearth)
298void KSComet::findMagnitude(const KSNumbers *)
299{
300 setMag(M1 + 2.5 * K1 * log10(rsun()) + 5 * log10(rearth()));
301}
302
303void KSComet::setEarthMOID(double earth_moid)
304{
305 EarthMOID = earth_moid;
306}
307
308void KSComet::setAlbedo(float albedo)
309{
310 Albedo = albedo;
311}
312
313void KSComet::setDiameter(float diam)
314{
315 Diameter = diam;
316}
317
319{
320 Dimensions = dim;
321}
322
323void KSComet::setNEO(bool neo)
324{
325 NEO = neo;
326}
327
329{
330 OrbitClass = orbit_class;
331}
332
334{
335 OrbitID = orbit_id;
336}
337
338void KSComet::setPeriod(float per)
339{
340 Period = per;
341}
342
343void KSComet::setRotationPeriod(float rot_per)
344{
345 RotationPeriod = rot_per;
346}
347
348//Unused virtual function from KSPlanetBase
350{
351 return false;
352}
353
355{
356 return solarsysUID(UID_SOL_COMET) | uidPart;
357}
A subclass of KSPlanetBase that implements comets.
Definition kscomet.h:44
void setPeriod(float per)
Sets the comet's period.
Definition kscomet.cpp:338
void setDimensions(QString dim)
Sets the comet's dimensions.
Definition kscomet.cpp:318
void setAlbedo(float albedo)
Sets the comet's albedo.
Definition kscomet.cpp:308
SkyObject::UID getUID() const override
Return UID for object.
Definition kscomet.cpp:354
bool loadData() override
Unused virtual function inherited from KSPlanetBase thus it's simply empty here.
Definition kscomet.cpp:349
void setOrbitID(QString orbit_id)
Sets the comet's orbit solution ID.
Definition kscomet.cpp:333
KSComet * clone() const override
Create copy of object.
Definition kscomet.cpp:132
void setDiameter(float diam)
Sets the comet's diameter.
Definition kscomet.cpp:313
void setEarthMOID(double earth_moid)
Sets the comet's earth minimum orbit intersection distance.
Definition kscomet.cpp:303
void setRotationPeriod(float rot_per)
Sets the comet's rotation period.
Definition kscomet.cpp:343
void setOrbitClass(QString orbit_class)
Sets the comet's orbit class.
Definition kscomet.cpp:328
bool findGeocentricPosition(const KSNumbers *num, const KSPlanetBase *Earth=nullptr) override
Calculate the geocentric RA, Dec coordinates of the Comet.
Definition kscomet.cpp:156
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 findPhysicalParameters()
Estimate physical parameters of the comet such as coma size, tail length and size of the nucleus.
Definition kscomet.cpp:138
void setNEO(bool neo)
Sets if the comet is a near earth object.
Definition kscomet.cpp:323
There are several time-dependent values used in position calculations, that are not specific to an ob...
Definition ksnumbers.h:43
const CachingDms * obliquity() const
Definition ksnumbers.h:56
A subclass of TrailObject that provides additional information needed for most solar system objects.
UID solarsysUID(UID type) const
Compute high 32-bits of UID.
void setPhysicalSize(double size)
set the planet's physical size, in km.
static const UID UID_SOL_COMET
Comets.
void setRearth(double r)
Set the distance from Earth, in AU.
const dms & ecLat() const
const dms & ecLong() const
double rearth() const
void EclipticToEquatorial(const CachingDms *Obliquity)
Convert Ecliptic longitude/latitude to Right Ascension/Declination.
double rsun() const
void setRsun(double r)
Set the solar distance in AU.
void setMag(float m)
Set the object's sorting magnitude.
Definition skyobject.h:416
qint64 UID
Type for Unique object IDenticator.
Definition skyobject.h:49
int type(void) const
Definition skyobject.h:189
void setType(int t)
Set the object's type identifier to the argument.
Definition skyobject.h:196
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 CachingDms & dec() const
Definition skypoint.h:269
const CachingDms & ra() const
Definition skypoint.h:263
void setRA0(dms r)
Sets RA0, the catalog Right Ascension.
Definition skypoint.h:94
void setDec0(dms d)
Sets Dec0, the catalog Declination.
Definition skypoint.h:119
An angle, stored as degrees, but expressible in many ways.
Definition dms.h:38
void reduceToRange(enum dms::AngleRanges range)
Reduce this angle to the given range.
Definition dms.cpp:446
const dms reduce() const
return the equivalent angle between 0 and 360 degrees.
Definition dms.cpp:251
void SinCos(double &s, double &c) const
Compute Sine and Cosine of the angle simultaneously.
Definition dms.h:447
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
virtual void setRadians(const double &Rad)
Set angle according to the argument, in radians.
Definition dms.h:333
const double & Degrees() const
Definition dms.h:141
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
char toLatin1() const const
QRegularExpressionMatch match(QStringView subjectView, qsizetype offset, MatchType matchType, MatchOptions matchOptions) const const
QString captured(QStringView name) const const
bool hasMatch() const const
const QChar at(qsizetype position) const const
int toInt(bool *ok, int base) const const
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Sat Dec 21 2024 17:04:47 by doxygen 1.12.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.