Kstars

ksasteroid.cpp
1 /*
2  SPDX-FileCopyrightText: 2001 Jason Harris <[email protected]>
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 
22 KSAsteroid::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 
150 void 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 
160 void KSAsteroid::setPerihelion(double perihelion)
161 {
162  q = perihelion;
163 }
164 
165 void KSAsteroid::setEarthMOID(double earth_moid)
166 {
167  EarthMOID = earth_moid;
168 }
169 
170 void KSAsteroid::setAlbedo(float albedo)
171 {
172  Albedo = albedo;
173 }
174 
175 void KSAsteroid::setDiameter(float diam)
176 {
177  Diameter = diam;
178 }
179 
181 {
182  Dimensions = dim;
183 }
184 
185 void KSAsteroid::setNEO(bool neo)
186 {
187  NEO = neo;
188 }
189 
191 {
192  OrbitClass = orbit_class;
193 }
194 
196 {
197  OrbitID = orbit_id;
198 }
199 
200 void KSAsteroid::setPeriod(float per)
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 
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 
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 
263 void KSAsteroid::setRotationPeriod(float rot_per)
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 }
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 setOrbitID(QString orbit_id)
Sets the asteroid's orbit solution ID.
Definition: ksasteroid.cpp:195
void setDimensions(QString dim)
Sets the asteroid's dimensions.
Definition: ksasteroid.cpp:180
void setOrbitClass(QString orbit_class)
Sets the asteroid's orbit class.
Definition: ksasteroid.cpp:190
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 setAlbedo(float albedo)
Sets the asteroid's albedo.
Definition: ksasteroid.cpp:170
void SinCos(double &s, double &c) const
Compute Sine and Cosine of the angle simultaneously.
Definition: dms.h:444
A subclass of KSPlanetBase that implements asteroids.
Definition: ksasteroid.h:41
double rsun() const
Definition: ksplanetbase.h:130
float mag() const
Definition: skyobject.h:206
KSAsteroid * clone() const override
Create copy of object.
Definition: ksasteroid.cpp:32
KCALENDARCORE_EXPORT QDataStream & operator>>(QDataStream &in, const KCalendarCore::Alarm::Ptr &)
double rearth() const
Definition: ksplanetbase.h:139
void setRA0(dms r)
Sets RA0, the catalog Right Ascension.
Definition: skypoint.h:94
void setEarthMOID(double earth_moid)
Sets the asteroid's earth minimum orbit intersection distance.
Definition: ksasteroid.cpp:165
void setDiameter(float diam)
Sets the asteroid's diameter.
Definition: ksasteroid.cpp:175
qint64 UID
Type for Unique object IDenticator.
Definition: skyobject.h:49
bool toCalculate()
toCalculate
Definition: ksasteroid.cpp:205
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
bool loadData() override
This is inherited from KSPlanetBase.
Definition: ksasteroid.cpp:269
static const UID UID_SOL_ASTEROID
Asteroids.
Definition: ksplanetbase.h:218
void setRearth(double r)
Set the distance from Earth, in AU.
Definition: ksplanetbase.h:145
const dms & ecLong() const
Definition: ksplanetbase.h:91
SkyObject::UID getUID() const override
Return UID for object.
Definition: ksasteroid.cpp:274
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.
Definition: ksasteroid.cpp:22
void setNEO(bool neo)
Sets if the comet is a near earth object.
Definition: ksasteroid.cpp:185
void setRotationPeriod(float rot_per)
Sets the asteroid's rotation period.
Definition: ksasteroid.cpp:263
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:333
const CachingDms & ra() const
Definition: skypoint.h:263
const CachingDms * obliquity() const
Definition: ksnumbers.h:56
double radians() const
Express the angle in radians.
Definition: dms.h:325
QDebug operator<<(QDebug d, const QCPVector2D &vec)
Definition: qcustomplot.h:446
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
bool findGeocentricPosition(const KSNumbers *num, const KSPlanetBase *Earth=nullptr) override
Calculate the geocentric RA, Dec coordinates of the Asteroid.
Definition: ksasteroid.cpp:39
void setRsun(double r)
Set the solar distance in AU.
Definition: ksplanetbase.h:136
void setPerihelion(double perihelion)
Sets the asteroid's perihelion distance.
Definition: ksasteroid.cpp:160
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 EclipticToEquatorial(const CachingDms *Obliquity)
Convert Ecliptic longitude/latitude to Right Ascension/Declination.
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
void setPeriod(float per)
Sets the asteroid's period.
Definition: ksasteroid.cpp:200
This file is part of the KDE documentation.
Documentation copyright © 1996-2023 The KDE developers.
Generated on Sun Oct 1 2023 04:02:40 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.