Kstars

ksasteroid.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 "ksasteroid.h"
8
9#include "dms.h"
10#include "ksnumbers.h"
11#include "Options.h"
12#ifdef KSTARS_LITE
13#include "skymaplite.h"
14#else
15#include "skymap.h"
16#endif
17
18#include <qdebug.h>
19
20#include <typeinfo>
21
22KSAsteroid::KSAsteroid(int _catN, const QString &s, const QString &imfile, long double _JD, double _a, double _e,
23 dms _i, dms _w, dms _Node, dms _M, double _H, double _G)
24 : KSPlanetBase(s, imfile), catN(_catN), JD(_JD), a(_a), e(_e), i(_i), w(_w), M(_M), N(_Node), H(_H), G(_G)
25{
26 setType(SkyObject::ASTEROID);
27 setMag(H);
28 //Compute the orbital Period from Kepler's 3rd law:
29 P = 365.2568984 * pow(a, 1.5); //period in days
30}
31
33{
34 Q_ASSERT(typeid(this) ==
35 typeid(static_cast<const KSAsteroid *>(this))); // Ensure we are not slicing a derived class
36 return new KSAsteroid(*this);
37}
38
40{
41 // TODO: (Valentin) TOP LEVEL CONTROL OF CALCULATION FOR ALL OBJECTS
42 if(!toCalculate()){
43 return false;
44 }
45
46 //determine the mean anomaly for the desired date. This is the mean anomaly for the
47 //ephemeis epoch, plus the number of days between the desired date and ephemeris epoch,
48 //times the asteroid's mean daily motion (360/P):
49
50 // All elements are in the heliocentric ecliptic J2000 reference frame.
51 // Mean anomaly is supplied at the Epoch (which is JD here)
52
53 dms m = dms(double(M.Degrees() + (lastPrecessJD - JD) * 360.0 / P)).reduce();
54 double sinm, cosm;
55 m.SinCos(sinm, cosm);
56
57 //compute eccentric anomaly:
58 double E = m.Degrees() + e * 180.0 / dms::PI * sinm * (1.0 + e * cosm);
59
60 if (e > 0.005) //need more accurate approximation, iterate...
61 {
62 double E0;
63 int iter(0);
64 do
65 {
66 E0 = E;
67 iter++;
68 E = E0 -
69 (E0 - e * 180.0 / dms::PI * sin(E0 * dms::DegToRad) - m.Degrees()) / (1 - e * cos(E0 * dms::DegToRad));
70 } while (fabs(E - E0) > 0.001 && iter < 1000);
71 }
72
73 // Assert that the solution of the Kepler equation E = M + e sin E is accurate to about 0.1 arcsecond
74 //Q_ASSERT( fabs( E - ( m.Degrees() + ( e * 180.0 / dms::PI ) * sin( E * dms::DegToRad ) ) ) < 0.10/3600.0 );
75
76 double sinE, cosE;
77 dms E1(E);
78 E1.SinCos(sinE, cosE);
79
80 double xv = a * (cosE - e);
81 double yv = a * sqrt(1.0 - e * e) * sinE;
82
83 //v is the true anomaly; r is the distance from the Sun
84 double v = atan2(yv, xv) / dms::DegToRad;
85 double r = sqrt(xv * xv + yv * yv);
86
87 //vw is the sum of the true anomaly and the argument of perihelion
88 dms vw(v + w.Degrees());
89 double sinN, cosN, sinvw, cosvw, sini, cosi;
90
91 N.SinCos(sinN, cosN);
92 vw.SinCos(sinvw, cosvw);
93 i.SinCos(sini, cosi);
94
95 //xh, yh, zh are the heliocentric cartesian coords with the ecliptic plane congruent with zh=0.
96 double xh = r * (cosN * cosvw - sinN * sinvw * cosi);
97 double yh = r * (sinN * cosvw + cosN * sinvw * cosi);
98 double zh = r * (sinvw * sini);
99
100 //the spherical ecliptic coordinates:
101 double ELongRad = atan2(yh, xh);
102 double ELatRad = atan2(zh, r);
103
104 helEcPos.longitude.setRadians(ELongRad);
105 helEcPos.longitude.reduceToRange(dms::ZERO_TO_2PI);
106 helEcPos.latitude.setRadians(ELatRad);
107 setRsun(r);
108
109 if (Earth)
110 {
111 //xe, ye, ze are the Earth's heliocentric cartesian coords
112 double cosBe, sinBe, cosLe, sinLe;
113 Earth->ecLong().SinCos(sinLe, cosLe);
114 Earth->ecLat().SinCos(sinBe, cosBe);
115
116 double xe = Earth->rsun() * cosBe * cosLe;
117 double ye = Earth->rsun() * cosBe * sinLe;
118 double ze = Earth->rsun() * sinBe;
119
120 //convert to geocentric ecliptic coordinates by subtracting Earth's coords:
121 xh -= xe;
122 yh -= ye;
123 zh -= ze;
124 }
125
126 //the spherical geocentricecliptic coordinates:
127 ELongRad = atan2(yh, xh);
128 double rr = sqrt(xh * xh + yh * yh);
129 ELatRad = atan2(zh, rr);
130
131 ep.longitude.setRadians(ELongRad);
132 ep.longitude.reduceToRange(dms::ZERO_TO_2PI);
133 ep.latitude.setRadians(ELatRad);
134 if (Earth)
135 setRearth(Earth);
136
138
139 // JM 2017-09-10: The calculations above produce J2000 RESULTS
140 // So we have to precess as well
141 setRA0(ra());
142 setDec0(dec());
143 apparentCoord(J2000, lastPrecessJD);
144 //nutate(num);
145 //aberrate(num);
146
147 return true;
148}
149
150void KSAsteroid::findMagnitude(const KSNumbers *)
151{
152 double param = 5.0 * log10(rsun() * rearth());
153 double phase_rad = phase().radians();
154 double phi1 = exp(-3.33 * pow(tan(phase_rad / 2.0), 0.63));
155 double phi2 = exp(-0.187 * pow(tan(phase_rad / 2.0), 1.22));
156
157 setMag(H + param - 2.5 * log10((1.0 - G) * phi1 + G * phi2));
158}
159
160void KSAsteroid::setPerihelion(double perihelion)
161{
162 q = perihelion;
163}
164
165void KSAsteroid::setEarthMOID(double earth_moid)
166{
167 EarthMOID = earth_moid;
168}
169
170void KSAsteroid::setAlbedo(float albedo)
171{
172 Albedo = albedo;
173}
174
176{
177 Diameter = diam;
178}
179
181{
182 Dimensions = dim;
183}
184
185void KSAsteroid::setNEO(bool neo)
186{
187 NEO = neo;
188}
189
191{
192 OrbitClass = orbit_class;
193}
194
196{
197 OrbitID = orbit_id;
198}
199
201{
202 Period = per;
203}
204
206{
207 // Filter by magnitude, but draw focused asteroids anyway :)
208 return ((mag() < Options::magLimitAsteroid())|| (std::isnan(mag()) != 0) ||
209#ifdef KSTARS_LITE
210 SkyMapLite::Instance()->focusObject() == this
211#else
212 SkyMap::Instance()->focusObject() == this
213#endif
214 );
215}
216
217QDataStream &operator<<(QDataStream &out, const KSAsteroid &asteroid)
218{
219 out << asteroid.Name << asteroid.OrbitClass << asteroid.Dimensions << asteroid.OrbitID
220 << asteroid.catN << static_cast<double>(asteroid.JD) << asteroid.a
221 << asteroid.e << asteroid.i << asteroid.w << asteroid.N
222 << asteroid.M << asteroid.H << asteroid.G << asteroid.q
223 << asteroid.NEO << asteroid.Diameter
224 << asteroid.Albedo << asteroid.RotationPeriod
225 << asteroid.Period << asteroid.EarthMOID;
226 return out;
227}
228
229QDataStream &operator>>(QDataStream &in, KSAsteroid *&asteroid)
230{
231 QString name, orbit_id, orbit_class, dimensions;
232 double q, a, e, H, G, earth_moid;
233 dms i, w, N, M;
234 double JD;
235 float diameter, albedo, rot_period, period;
236 bool neo;
237 int catN;
238
239 in >> name;
240 in >> orbit_class;
241 in >> dimensions;
242 in >> orbit_id;
243
244 in >> catN >> JD >> a >> e >> i >> w >> N >> M >> H >> G >> q >> neo >> diameter
245 >> albedo >> rot_period >> period >> earth_moid;
246
247 asteroid = new KSAsteroid(catN, name, QString(), JD, a, e, i, w, N, M, H, G);
248 asteroid->setPerihelion(q);
249 asteroid->setOrbitID(orbit_id);
250 asteroid->setNEO(neo);
251 asteroid->setDiameter(diameter);
252 asteroid->setDimensions(dimensions);
253 asteroid->setAlbedo(albedo);
254 asteroid->setRotationPeriod(rot_period);
255 asteroid->setPeriod(period);
256 asteroid->setEarthMOID(earth_moid);
257 asteroid->setOrbitClass(orbit_class);
258 asteroid->setPhysicalSize(diameter);
259
260 return in;
261}
262
264{
265 RotationPeriod = rot_per;
266}
267
268//Unused virtual function from KSPlanetBase
270{
271 return false;
272}
273
275{
276 return solarsysUID(UID_SOL_ASTEROID) | catN;
277}
A subclass of KSPlanetBase that implements asteroids.
Definition ksasteroid.h:42
SkyObject::UID getUID() const override
Return UID for object.
void setPerihelion(double perihelion)
Sets the asteroid's perihelion distance.
KSAsteroid * clone() const override
Create copy of object.
void setPeriod(float per)
Sets the asteroid's period.
void setDiameter(float diam)
Sets the asteroid's diameter.
void setDimensions(QString dim)
Sets the asteroid's dimensions.
void setEarthMOID(double earth_moid)
Sets the asteroid's earth minimum orbit intersection distance.
void setNEO(bool neo)
Sets if the comet is a near earth object.
void setOrbitID(QString orbit_id)
Sets the asteroid's orbit solution ID.
bool toCalculate()
toCalculate
void setAlbedo(float albedo)
Sets the asteroid's albedo.
bool findGeocentricPosition(const KSNumbers *num, const KSPlanetBase *Earth=nullptr) override
Calculate the geocentric RA, Dec coordinates of the Asteroid.
bool loadData() override
This is inherited from KSPlanetBase.
KSAsteroid(int catN, const QString &s, const QString &image_file, long double JD, double a, double e, dms i, dms w, dms N, dms M, double H, double G)
Constructor.
void setRotationPeriod(float rot_per)
Sets the asteroid's rotation period.
void setOrbitClass(QString orbit_class)
Sets the asteroid's orbit class.
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_ASTEROID
Asteroids.
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:403
qint64 UID
Type for Unique object IDenticator.
Definition skyobject.h:49
float mag() const
Definition skyobject.h:206
void setType(int t)
Set the object's type identifier to the argument.
Definition skyobject.h:195
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
double radians() const
Express the angle in radians.
Definition dms.h:325
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
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Wed Nov 6 2024 12:13:16 by doxygen 1.12.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.