File indexing completed on 2024-04-21 03:44:40
0001 /* 0002 SPDX-FileCopyrightText: 2001 Jason Harris <jharris@30doradus.org> 0003 0004 SPDX-License-Identifier: GPL-2.0-or-later 0005 */ 0006 0007 #include "ksasteroid.h" 0008 0009 #include "dms.h" 0010 #include "ksnumbers.h" 0011 #include "Options.h" 0012 #ifdef KSTARS_LITE 0013 #include "skymaplite.h" 0014 #else 0015 #include "skymap.h" 0016 #endif 0017 0018 #include <qdebug.h> 0019 0020 #include <typeinfo> 0021 0022 KSAsteroid::KSAsteroid(int _catN, const QString &s, const QString &imfile, long double _JD, double _a, double _e, 0023 dms _i, dms _w, dms _Node, dms _M, double _H, double _G) 0024 : KSPlanetBase(s, imfile), catN(_catN), JD(_JD), a(_a), e(_e), i(_i), w(_w), M(_M), N(_Node), H(_H), G(_G) 0025 { 0026 setType(SkyObject::ASTEROID); 0027 setMag(H); 0028 //Compute the orbital Period from Kepler's 3rd law: 0029 P = 365.2568984 * pow(a, 1.5); //period in days 0030 } 0031 0032 KSAsteroid *KSAsteroid::clone() const 0033 { 0034 Q_ASSERT(typeid(this) == 0035 typeid(static_cast<const KSAsteroid *>(this))); // Ensure we are not slicing a derived class 0036 return new KSAsteroid(*this); 0037 } 0038 0039 bool KSAsteroid::findGeocentricPosition(const KSNumbers *num, const KSPlanetBase *Earth) 0040 { 0041 // TODO: (Valentin) TOP LEVEL CONTROL OF CALCULATION FOR ALL OBJECTS 0042 if(!toCalculate()){ 0043 return false; 0044 } 0045 0046 //determine the mean anomaly for the desired date. This is the mean anomaly for the 0047 //ephemeis epoch, plus the number of days between the desired date and ephemeris epoch, 0048 //times the asteroid's mean daily motion (360/P): 0049 0050 // All elements are in the heliocentric ecliptic J2000 reference frame. 0051 // Mean anomaly is supplied at the Epoch (which is JD here) 0052 0053 dms m = dms(double(M.Degrees() + (lastPrecessJD - JD) * 360.0 / P)).reduce(); 0054 double sinm, cosm; 0055 m.SinCos(sinm, cosm); 0056 0057 //compute eccentric anomaly: 0058 double E = m.Degrees() + e * 180.0 / dms::PI * sinm * (1.0 + e * cosm); 0059 0060 if (e > 0.005) //need more accurate approximation, iterate... 0061 { 0062 double E0; 0063 int iter(0); 0064 do 0065 { 0066 E0 = E; 0067 iter++; 0068 E = E0 - 0069 (E0 - e * 180.0 / dms::PI * sin(E0 * dms::DegToRad) - m.Degrees()) / (1 - e * cos(E0 * dms::DegToRad)); 0070 } while (fabs(E - E0) > 0.001 && iter < 1000); 0071 } 0072 0073 // Assert that the solution of the Kepler equation E = M + e sin E is accurate to about 0.1 arcsecond 0074 //Q_ASSERT( fabs( E - ( m.Degrees() + ( e * 180.0 / dms::PI ) * sin( E * dms::DegToRad ) ) ) < 0.10/3600.0 ); 0075 0076 double sinE, cosE; 0077 dms E1(E); 0078 E1.SinCos(sinE, cosE); 0079 0080 double xv = a * (cosE - e); 0081 double yv = a * sqrt(1.0 - e * e) * sinE; 0082 0083 //v is the true anomaly; r is the distance from the Sun 0084 double v = atan2(yv, xv) / dms::DegToRad; 0085 double r = sqrt(xv * xv + yv * yv); 0086 0087 //vw is the sum of the true anomaly and the argument of perihelion 0088 dms vw(v + w.Degrees()); 0089 double sinN, cosN, sinvw, cosvw, sini, cosi; 0090 0091 N.SinCos(sinN, cosN); 0092 vw.SinCos(sinvw, cosvw); 0093 i.SinCos(sini, cosi); 0094 0095 //xh, yh, zh are the heliocentric cartesian coords with the ecliptic plane congruent with zh=0. 0096 double xh = r * (cosN * cosvw - sinN * sinvw * cosi); 0097 double yh = r * (sinN * cosvw + cosN * sinvw * cosi); 0098 double zh = r * (sinvw * sini); 0099 0100 //the spherical ecliptic coordinates: 0101 double ELongRad = atan2(yh, xh); 0102 double ELatRad = atan2(zh, r); 0103 0104 helEcPos.longitude.setRadians(ELongRad); 0105 helEcPos.longitude.reduceToRange(dms::ZERO_TO_2PI); 0106 helEcPos.latitude.setRadians(ELatRad); 0107 setRsun(r); 0108 0109 if (Earth) 0110 { 0111 //xe, ye, ze are the Earth's heliocentric cartesian coords 0112 double cosBe, sinBe, cosLe, sinLe; 0113 Earth->ecLong().SinCos(sinLe, cosLe); 0114 Earth->ecLat().SinCos(sinBe, cosBe); 0115 0116 double xe = Earth->rsun() * cosBe * cosLe; 0117 double ye = Earth->rsun() * cosBe * sinLe; 0118 double ze = Earth->rsun() * sinBe; 0119 0120 //convert to geocentric ecliptic coordinates by subtracting Earth's coords: 0121 xh -= xe; 0122 yh -= ye; 0123 zh -= ze; 0124 } 0125 0126 //the spherical geocentricecliptic coordinates: 0127 ELongRad = atan2(yh, xh); 0128 double rr = sqrt(xh * xh + yh * yh); 0129 ELatRad = atan2(zh, rr); 0130 0131 ep.longitude.setRadians(ELongRad); 0132 ep.longitude.reduceToRange(dms::ZERO_TO_2PI); 0133 ep.latitude.setRadians(ELatRad); 0134 if (Earth) 0135 setRearth(Earth); 0136 0137 EclipticToEquatorial(num->obliquity()); 0138 0139 // JM 2017-09-10: The calculations above produce J2000 RESULTS 0140 // So we have to precess as well 0141 setRA0(ra()); 0142 setDec0(dec()); 0143 apparentCoord(J2000, lastPrecessJD); 0144 //nutate(num); 0145 //aberrate(num); 0146 0147 return true; 0148 } 0149 0150 void KSAsteroid::findMagnitude(const KSNumbers *) 0151 { 0152 double param = 5.0 * log10(rsun() * rearth()); 0153 double phase_rad = phase().radians(); 0154 double phi1 = exp(-3.33 * pow(tan(phase_rad / 2.0), 0.63)); 0155 double phi2 = exp(-0.187 * pow(tan(phase_rad / 2.0), 1.22)); 0156 0157 setMag(H + param - 2.5 * log10((1.0 - G) * phi1 + G * phi2)); 0158 } 0159 0160 void KSAsteroid::setPerihelion(double perihelion) 0161 { 0162 q = perihelion; 0163 } 0164 0165 void KSAsteroid::setEarthMOID(double earth_moid) 0166 { 0167 EarthMOID = earth_moid; 0168 } 0169 0170 void KSAsteroid::setAlbedo(float albedo) 0171 { 0172 Albedo = albedo; 0173 } 0174 0175 void KSAsteroid::setDiameter(float diam) 0176 { 0177 Diameter = diam; 0178 } 0179 0180 void KSAsteroid::setDimensions(QString dim) 0181 { 0182 Dimensions = dim; 0183 } 0184 0185 void KSAsteroid::setNEO(bool neo) 0186 { 0187 NEO = neo; 0188 } 0189 0190 void KSAsteroid::setOrbitClass(QString orbit_class) 0191 { 0192 OrbitClass = orbit_class; 0193 } 0194 0195 void KSAsteroid::setOrbitID(QString orbit_id) 0196 { 0197 OrbitID = orbit_id; 0198 } 0199 0200 void KSAsteroid::setPeriod(float per) 0201 { 0202 Period = per; 0203 } 0204 0205 bool KSAsteroid::toCalculate() 0206 { 0207 // Filter by magnitude, but draw focused asteroids anyway :) 0208 return ((mag() < Options::magLimitAsteroid())|| (std::isnan(mag()) != 0) || 0209 #ifdef KSTARS_LITE 0210 SkyMapLite::Instance()->focusObject() == this 0211 #else 0212 SkyMap::Instance()->focusObject() == this 0213 #endif 0214 ); 0215 } 0216 0217 QDataStream &operator<<(QDataStream &out, const KSAsteroid &asteroid) 0218 { 0219 out << asteroid.Name << asteroid.OrbitClass << asteroid.Dimensions << asteroid.OrbitID 0220 << asteroid.catN << static_cast<double>(asteroid.JD) << asteroid.a 0221 << asteroid.e << asteroid.i << asteroid.w << asteroid.N 0222 << asteroid.M << asteroid.H << asteroid.G << asteroid.q 0223 << asteroid.NEO << asteroid.Diameter 0224 << asteroid.Albedo << asteroid.RotationPeriod 0225 << asteroid.Period << asteroid.EarthMOID; 0226 return out; 0227 } 0228 0229 QDataStream &operator>>(QDataStream &in, KSAsteroid *&asteroid) 0230 { 0231 QString name, orbit_id, orbit_class, dimensions; 0232 double q, a, e, H, G, earth_moid; 0233 dms i, w, N, M; 0234 double JD; 0235 float diameter, albedo, rot_period, period; 0236 bool neo; 0237 int catN; 0238 0239 in >> name; 0240 in >> orbit_class; 0241 in >> dimensions; 0242 in >> orbit_id; 0243 0244 in >> catN >> JD >> a >> e >> i >> w >> N >> M >> H >> G >> q >> neo >> diameter 0245 >> albedo >> rot_period >> period >> earth_moid; 0246 0247 asteroid = new KSAsteroid(catN, name, QString(), JD, a, e, i, w, N, M, H, G); 0248 asteroid->setPerihelion(q); 0249 asteroid->setOrbitID(orbit_id); 0250 asteroid->setNEO(neo); 0251 asteroid->setDiameter(diameter); 0252 asteroid->setDimensions(dimensions); 0253 asteroid->setAlbedo(albedo); 0254 asteroid->setRotationPeriod(rot_period); 0255 asteroid->setPeriod(period); 0256 asteroid->setEarthMOID(earth_moid); 0257 asteroid->setOrbitClass(orbit_class); 0258 asteroid->setPhysicalSize(diameter); 0259 0260 return in; 0261 } 0262 0263 void KSAsteroid::setRotationPeriod(float rot_per) 0264 { 0265 RotationPeriod = rot_per; 0266 } 0267 0268 //Unused virtual function from KSPlanetBase 0269 bool KSAsteroid::loadData() 0270 { 0271 return false; 0272 } 0273 0274 SkyObject::UID KSAsteroid::getUID() const 0275 { 0276 return solarsysUID(UID_SOL_ASTEROID) | catN; 0277 }