1 /* 2 * Stellarium 3 * Copyright (C) 2002 Fabien Chereau 4 * Copyright (C) 2021 Georg Zotti 5 * 6 * This program is free software; you can redistribute it and/or 7 * modify it under the terms of the GNU General Public License 8 * as published by the Free Software Foundation; either version 2 9 * of the License, or (at your option) any later version. 10 * 11 * This program is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 * GNU General Public License for more details. 15 * 16 * You should have received a copy of the GNU General Public License 17 * along with this program; if not, write to the Free Software 18 * Foundation, Inc., 51 Franklin Street, Suite 500, Boston, MA 02110-1335, USA. 19 */ 20 21 #ifndef ROTATIONELEMENTS_HPP 22 #define ROTATIONELEMENTS_HPP 23 24 //! Class used to store rotational elements, i.e. axis orientation for the planetary body. 25 //! Data are read from ssystem_(major|minor).ini in SolarSystem.cpp::loadPlanets(). 26 //! 27 //! Notes of early 2021: The original planet axis model in Stellarium seems to be an outdated model untraceable in the literature. 28 //! Stellarium used to have elements w.r.t. VSOP87 (ecliptical) frame for planets, and elements w.r.t. the host planet's equator for planet moons: 29 //! ascendingNode 30 //! obliquity 31 //! offset 32 //! period 33 //! We cannot remove these old data as long as some Moon orbits for the outer planets depend on them. 34 //! The problem is that without a reference we have no idea whether any of these elements and algorithms are correct. 35 //! 36 //! IAU standards like the Report of the IAU Working Group on Cartographic Coordinates and Rotational elements 2009 (DOI:10.1007/s10569-010-9320-4) 37 //! has axes given w.r.t. J2000 ICRF. 38 //! 39 //! Between V0.15 and 0.20.*, if the planet elements in ssystem.ini had elements "rot_pole_ra" and "rot_pole_de" given, the poles were transformed to 40 //! ecliptically-based directions (VSOP87) obliquity/ascendingNode. But these were only ra/de_J2000 poles. Some axes have precession, which need to be modelled/updated. 41 //! The new way (0.21+): We still allow using the previous 4 elements in Planet::computeTransMatrix(.,.) and Planet::getSiderealTime(), but only if WGCCRE elements are not available. 42 //! 43 //! ra=ra0+T*ra1 44 //! de=de0+T*de1 ( --> obliquity, ascendingNode) 45 //! rot_rotation_offset [degrees] =W0 46 //! rot_periode [hours] = computed from rot_pole_W1[deg/day] if that exists. 360 [deg] / _rot_ [hours] --> 360 * _rot_ / 24 [deg/hours] 47 //! 48 //! If rot_pole... values are given, then they are ICRF and transformed on the fly to VSOP87, stored in here. 49 //! 50 //! Since 0.21 we use the WGCCRE elements and orientations directly, so that axes behave properly. 51 //! In addition, the objects with more complicated element behaviour are updated with two special functions. 52 //! New keys in ssystem_*.ini, their storage in RotationElements, and their equivalents in the IAU report: 53 //! rot_pole_ra [degrees] re.ra0 [rad] constant term for alpha_0 54 //! rot_pole_de [degrees] re.ra1 [rad/ct] constant term for delta_0 55 //! rot_pole_ra1 [degrees/cy] re.de0 [rad] T factor for alpha_0 56 //! rot_pole_de1 [degrees/cy] re.de1 [rad/ct] T factor for delta_0 57 //! rot_pole_W0 [degrees] re.W0 [degrees] constant term fo W. 58 //! rot_pole_W1 [degrees/day] re.W1 [degrees/day] d factor for W 59 60 #include <QString> 61 #include "VecMath.hpp" 62 63 // epoch J2000: 12 UT on 1 Jan 2000 64 #define J2000 2451545.0 65 66 class RotationElements 67 { 68 public: 69 enum ComputationMethod { 70 Traditional, // Stellarium prior to 0.21, This is unfortunately not documented. Orbits and axes given w.r.t. parent axes. 71 WGCCRE // Orientation as described by the IAU Working Group on Cartographic Coordinates and Rotational Elements reports of 2009 or later 72 }; 73 enum PlanetCorrection // To be used as named argument for correction calculations regarding orientations of planetary axes. 74 { // See updatePlanetCorrections() 75 EarthMoon=3, 76 Mars=4, 77 Jupiter=5, 78 Saturn=6, 79 Uranus=7, 80 Neptune=8 81 }; 82 RotationElements(void)83 RotationElements(void) : period(1.), offset(0.), epoch(J2000), obliquity(0.), ascendingNode(0.), 84 method(Traditional), ra0(0.), ra1(0.), de0(0.), de1(0.), W0(0.), W1(0.), 85 currentAxisRA(0.), currentAxisDE(0.), currentAxisW(0.) {} 86 double period; // [deprecated] (sidereal) rotation period [earth days] 87 double offset; // [deprecated] rotation at epoch [degrees] 88 double epoch; // JDE (JD TT) of epoch for these elements 89 double obliquity; // [deprecated] tilt of rotation axis w.r.t. ecliptic [radians] 90 double ascendingNode; // [deprecated] long. of ascending node of equator on the ecliptic [radians] 91 // new elements for 0.21+: The 6/9 new entries after the switch are enough for many objects. More corrections can be applied where required. 92 ComputationMethod method; // The reference system in use for the respective object. WGCCRE is preferred, but we don't have all data. 93 double ra0; // [rad] RA_0 right ascension of north pole. ssystem.ini: rot_pole_ra /180*M_PI 94 double ra1; // [rad/century] rate of change in axis ra ssystem.ini: rot_pole_ra1 /180*M_PI 95 double de0; // [rad] DE_0 declination of north pole ssystem.ini: rot_pole_de /180*M_PI 96 double de1; // [rad/century] rate of change in axis de ssystem.ini: rot_pole_de1 /180*M_PI 97 // These values are only in the modern algorithms. invalid if ra0=0. 98 double W0; // [deg] mean longitude of prime meridian along equator measured from intersection with ICRS plane at epoch. 99 double W1; // [deg/d] mean longitude motion. W=W0+d*W1. 100 double currentAxisRA; // [rad] Mostly for infostring: RA=RA0+d*RA1(+corrections) 101 double currentAxisDE; // [rad] Mostly for infostring: DE=DE0+d*DE1(+corrections) 102 double currentAxisW; // [deg] Mostly for infostring: W =W0 +d*W1 (+corrections) 103 104 //! 0.21+: Axes of planets and moons require terms depending on T=(jde-J2000)/36525, described in Explanatory Supplement 2013, Tables 10.1 and 10.10-14, 105 //! updated in WGCCRE reports 2009 and 2015. 106 //! Others require frequent updates, depending on jde-J2000. (Moon etc.) 107 //! These should be updated as frequently as needed, optimally with the planet. Light time correction should be applied when needed. 108 //! best place to call update is the SolarSystem::computePlanets() 109 struct PlanetCorrections { 110 double JDE_E; // keep record of when these values are valid: Earth 111 double JDE_M; // keep record of when these values are valid: Mars 112 double JDE_J; // keep record of when these values are valid: Jupiter 113 double JDE_S; // keep record of when these values are valid: Saturn 114 double JDE_U; // keep record of when these values are valid: Uranus 115 double JDE_N; // keep record of when these values are valid: Neptune 116 117 double E1; // Earth corrections. These are from WGCCRE2009. 118 double E2; 119 double E3; 120 double E4; 121 double E5; 122 double E6; 123 double E7; 124 double E8; 125 double E9; 126 double E10; 127 double E11; 128 double E12; 129 double E13; 130 double M1; // Mars corrections. These are from WGCCRE2015. 131 double M2; 132 double M3; 133 double M4; 134 double M5; 135 double M6; 136 double M7; 137 double M8; 138 double M9; 139 double M10; 140 double Ja; // Jupiter axis terms, Table 10.1 141 double Jb; 142 double Jc; 143 double Jd; 144 double Je; 145 double Na; // Neptune axix term 146 double J1; // corrective terms for Jupiter' moons, Table 10.10 147 double J2; 148 double J3; 149 double J4; 150 double J5; 151 double J6; 152 double J7; 153 double J8; 154 double S1; // corrective terms for Saturn's moons, Table 10.12 155 double S2; 156 double S3; 157 double S4; 158 double S5; 159 double S6; 160 double U1; // for Cordelia // corrective terms for Uranus's moons, Table 10.14. 161 double U2; // for Ophelia 162 double U3; // for Bianca (not in 0.20) 163 double U4; // for Cressida 164 double U5; // for Desdemona 165 double U6; // for Juliet 166 double U7; // for Portia (not in 0.20) 167 double U8; // for Rosalind (not in 0.20) 168 double U9; // for Belinda (not in 0.20) 169 double U10;// for Puck 170 double U11; 171 double U12; 172 double U13; 173 double U14; 174 double U15; 175 double U16; 176 double N1; // corrective terms for Neptune's moons, Table 10.15 (N=Na!) 177 double N2; 178 double N3; 179 double N4; 180 double N5; 181 double N6; 182 double N7; 183 }; 184 static PlanetCorrections planetCorrections; 185 //! Update the planet corrections. planet=3(Moon), 4 (Mars), 5(Jupiter), 6(Saturn), 7(Uranus), 8(Neptune). 186 //! The values are immediately converted to radians. 187 static void updatePlanetCorrections(const double JDE, const PlanetCorrection planet); 188 189 //! Finetuning correction functions. These functions are attached to the planets in Planet::setRotationElements() 190 //! and called at the right places. This avoids ugly sequences of name query ifelses. 191 //! arguments are d=JDE-J2000, T=d/36525, J2000NPoleRA, J2000NPoleDE 192 typedef double (*axisRotFuncType)(const double, const double); 193 typedef void (*axisOriFuncType)(const double, const double, double*, double*); 194 axisRotFuncType corrW; 195 axisOriFuncType corrOri; 196 static const QMap<QString, axisRotFuncType>axisRotCorrFuncMap; 197 static const QMap<QString, axisOriFuncType>axisOriCorrFuncMap; 198 199 //! These corr* functions can be used as corrW. corrWnil(const double d,const double T)200 static double corrWnil(const double d, const double T){Q_UNUSED(d) Q_UNUSED(T) return 0;} 201 static double corrWMoon(const double d, const double T); 202 static double corrWMercury(const double d, const double T); 203 static double corrWMars(const double d, const double T); 204 //! The default W delivers SystemII longitude. 205 //! We have to shift by GRS position and texture position. 206 //! The final value will no longer be W but the rotation value required to show the GRS. 207 static double corrWJupiter(const double d, const double T); 208 static double corrWNeptune(const double d, const double T); 209 static double corrWPhobos(const double d, const double T); 210 static double corrWDeimos(const double d, const double T); 211 static double corrWIo(const double d, const double T); 212 static double corrWEuropa(const double d, const double T); 213 static double corrWGanymede(const double d, const double T); 214 static double corrWCallisto(const double d, const double T); 215 static double corrWAmalthea(const double d, const double T); 216 static double corrWThebe(const double d, const double T); 217 static double corrWMimas(const double d, const double T); 218 static double corrWTethys(const double d, const double T); 219 static double corrWRhea(const double d, const double T); 220 static double corrWJanus(const double d, const double T); 221 static double corrWEpimetheus(const double d, const double T); 222 static double corrWCordelia(const double d, const double T); 223 static double corrWOphelia(const double d, const double T); 224 static double corrWBianca(const double d, const double T); 225 static double corrWCressida(const double d, const double T); 226 static double corrWDesdemona(const double d, const double T); 227 static double corrWJuliet(const double d, const double T); 228 static double corrWPortia(const double d, const double T); 229 static double corrWRosalind(const double d, const double T); 230 static double corrWBelinda(const double d, const double T); 231 static double corrWPuck(const double d, const double T); 232 static double corrWAriel(const double d, const double T); 233 static double corrWUmbriel(const double d, const double T); 234 static double corrWTitania(const double d, const double T); 235 static double corrWOberon(const double d, const double T); 236 static double corrWMiranda(const double d, const double T); 237 static double corrWTriton(const double d, const double T); 238 static double corrWNaiad(const double d, const double T); 239 static double corrWThalassa(const double d, const double T); 240 static double corrWDespina(const double d, const double T); 241 static double corrWGalatea(const double d, const double T); 242 static double corrWLarissa(const double d, const double T); 243 static double corrWProteus(const double d, const double T); 244 245 //! These functions can be used as corrOri. corrOriNil(const double,const double,double *,double *)246 static void corrOriNil(const double, const double, double*, double*){} // Do nothing. 247 static void corrOriMoon(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 248 static void corrOriMars(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 249 static void corrOriJupiter(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 250 static void corrOriNeptune(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 251 static void corrOriPhobos(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 252 static void corrOriDeimos(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 253 static void corrOriIo(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 254 static void corrOriEuropa(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 255 static void corrOriGanymede(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 256 static void corrOriCallisto(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 257 static void corrOriAmalthea(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 258 static void corrOriThebe(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 259 static void corrOriMimas(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 260 static void corrOriTethys(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 261 static void corrOriRhea(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 262 static void corrOriJanus(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 263 static void corrOriEpimetheus(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 264 static void corrOriAriel(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 265 static void corrOriUmbriel(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 266 static void corrOriTitania(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 267 static void corrOriOberon(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 268 static void corrOriMiranda(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 269 static void corrOriCordelia(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 270 static void corrOriOphelia(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 271 static void corrOriBianca(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 272 static void corrOriCressida(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 273 static void corrOriDesdemona(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 274 static void corrOriJuliet(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 275 static void corrOriPortia(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 276 static void corrOriRosalind(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 277 static void corrOriBelinda(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 278 static void corrOriPuck(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 279 static void corrOriTriton(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 280 static void corrOriNaiad(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 281 static void corrOriThalassa(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 282 static void corrOriDespina(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 283 static void corrOriGalatea(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 284 static void corrOriLarissa(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 285 static void corrOriProteus(const double d, const double T, double* J2000NPoleRA, double* J2000NPoleDE); 286 287 static bool flagCustomGrsSettings; //!< Use custom settings for calculation of position of Great Red Spot? 288 static double customGrsJD; //!< Initial JD (epoch) for calculation of position of Great Red Spot 289 static double customGrsLongitude; //!< Longitude of Great Red Spot at customGrsJD (System II, degrees) 290 static double customGrsDrift; //!< Annual drift of Great Red Spot position (degrees) 291 292 //! Retrieve magnitude variation depending on angle Ls [radians]. 293 //! Source: A. Mallama: The magnitude and albedo of Mars. Icarus 192(2007) 404-416. 294 //! @arg albedo true to return longitudinal albedo correction, Ls=the average of sub-earth and sub-solar planetographic longitudes 295 //! false for the Orbital Longitude correction. Ls here is the apparent solar longitude along Mars' "ecliptic". (0=Martian Vernal Equinox) 296 static double getMarsMagLs(const double Ls, const bool albedo); 297 private: 298 static const QList<Vec3d> marsMagLs; 299 }; 300 301 #endif // ROTATIONELEMENTS_HPP 302