File indexing completed on 2024-04-21 14:47:05
0001 /* 0002 SPDX-FileCopyrightText: 2002-2005 Jason Harris <kstars@30doradus.org> 0003 SPDX-FileCopyrightText: 2004-2005 Pablo de Vicente <p.devicente@wanadoo.es> 0004 0005 SPDX-License-Identifier: GPL-2.0-or-later 0006 */ 0007 0008 #include "ksnumbers.h" 0009 0010 #include "kstarsdatetime.h" //for J2000 define 0011 0012 // 63 elements 0013 const int KSNumbers::arguments[NUTTERMS][5] = { 0014 { 0, 0, 0, 0, 1 }, { -2, 0, 0, 2, 2 }, { 0, 0, 0, 2, 2 }, { 0, 0, 0, 0, 2 }, { 0, 1, 0, 0, 0 }, 0015 { 0, 0, 1, 0, 0 }, { -2, 1, 0, 2, 2 }, { 0, 0, 0, 2, 1 }, { 0, 0, 1, 2, 2 }, { -2, -1, 0, 2, 2 }, 0016 { -2, 0, 1, 0, 0 }, { -2, 0, 0, 2, 1 }, { 0, 0, -1, 2, 2 }, { 2, 0, 0, 0, 0 }, { 0, 0, 1, 0, 1 }, 0017 { 2, 0, -1, 2, 2 }, { 0, 0, -1, 0, 1 }, { 0, 0, 1, 2, 1 }, { -2, 0, 2, 0, 0 }, { 0, 0, -2, 2, 1 }, 0018 { 2, 0, 0, 2, 2 }, { 0, 0, 2, 2, 2 }, { 0, 0, 2, 0, 0 }, { -2, 0, 1, 2, 2 }, { 0, 0, 0, 2, 0 }, 0019 { -2, 0, 0, 2, 0 }, { 0, 0, -1, 2, 1 }, { 0, 2, 0, 0, 0 }, { 2, 0, -1, 0, 1 }, { -2, 2, 0, 2, 2 }, 0020 { 0, 1, 0, 0, 1 }, { -2, 0, 1, 0, 1 }, { 0, -1, 0, 0, 1 }, { 0, 0, 2, -2, 0 }, { 2, 0, -1, 2, 1 }, 0021 { 2, 0, 1, 2, 2 }, { 0, 1, 0, 2, 2 }, { -2, 1, 1, 0, 0 }, { 0, -1, 0, 2, 2 }, { 2, 0, 0, 2, 1 }, 0022 { 2, 0, 1, 0, 0 }, { -2, 0, 2, 2, 2 }, { -2, 0, 1, 2, 1 }, { 2, 0, -2, 0, 1 }, { 2, 0, 0, 0, 1 }, 0023 { 0, -1, 1, 0, 0 }, { -2, -1, 0, 2, 1 }, { -2, 0, 0, 0, 1 }, { 0, 0, 2, 2, 1 }, { -2, 0, 2, 0, 1 }, 0024 { -2, 1, 0, 2, 1 }, { 0, 0, 1, -2, 0 }, { -1, 0, 1, 0, 0 }, { -2, 1, 0, 0, 0 }, { 1, 0, 0, 0, 0 }, 0025 { 0, 0, 1, 2, 0 }, { 0, 0, -2, 2, 2 }, { -1, -1, 1, 0, 0 }, { 0, 1, 1, 0, 0 }, { 0, -1, 1, 2, 2 }, 0026 { 2, -1, -1, 2, 2 }, { 0, 0, 3, 2, 2 }, { 2, -1, 0, 2, 2 } 0027 }; 0028 0029 const int KSNumbers::amp[NUTTERMS][4] = { { -171996, -1742, 92025, 89 }, 0030 { -13187, -16, 5736, -31 }, 0031 { -2274, -2, 977, -5 }, 0032 { 2062, 2, -895, 5 }, 0033 { 1426, -34, 54, -1 }, 0034 { 712, 1, -7, 0 }, 0035 { -517, 12, 224, -6 }, 0036 { -386, -4, 200, 0 }, 0037 { -301, 0, 129, -1 }, 0038 { 217, -5, -95, 3 }, 0039 { -158, 0, 0, 0 }, 0040 { 129, 1, -70, 0 }, 0041 { 123, 0, -53, 0 }, 0042 { 63, 0, 0, 0 }, 0043 { 63, 1, -33, 0 }, 0044 { -59, 0, 26, 0 }, 0045 { -58, -1, 32, 0 }, 0046 { -51, 0, 27, 0 }, 0047 { 48, 0, 0, 0 }, 0048 { 46, 0, -24, 0 }, 0049 { -38, 0, 16, 0 }, 0050 { -31, 0, 13, 0 }, 0051 { 29, 0, 0, 0 }, 0052 { 29, 0, -12, 0 }, 0053 { 26, 0, 0, 0 }, 0054 { -22, 0, 0, 0 }, 0055 { 21, 0, -10, 0 }, 0056 { 17, -1, 0, 0 }, 0057 { 16, 0, -8, 0 }, 0058 { -16, 1, 7, 0 }, 0059 { -15, 0, 9, 0 }, 0060 { -13, 0, 7, 0 }, 0061 { -12, 0, 6, 0 }, 0062 { 11, 0, 0, 0 }, 0063 { -10, 0, 5, 0 }, 0064 { -8, 0, 3, 0 }, 0065 { 7, 0, -3, 0 }, 0066 { -7, 0, 0, 0 }, 0067 { -7, 0, 3, 0 }, 0068 { -7, 0, 3, 0 }, 0069 { 6, 0, 0, 0 }, 0070 { 6, 0, -3, 0 }, 0071 { 6, 0, -3, 0 }, 0072 { -6, 0, 3, 0 }, 0073 { -6, 0, 3, 0 }, 0074 { 5, 0, 0, 0 }, 0075 { -5, 0, 3, 0 }, 0076 { -5, 0, 3, 0 }, 0077 { -5, 0, 3, 0 }, 0078 { 4, 0, 0, 0 }, 0079 { 4, 0, 0, 0 }, 0080 { 4, 0, 0, 0 }, 0081 { -4, 0, 0, 0 }, 0082 { -4, 0, 0, 0 }, 0083 { -4, 0, 0, 0 }, 0084 { 3, 0, 0, 0 }, 0085 { -3, 0, 0, 0 }, 0086 { -3, 0, 0, 0 }, 0087 { -3, 0, 0, 0 }, 0088 { -3, 0, 0, 0 }, 0089 { -3, 0, 0, 0 }, 0090 { -3, 0, 0, 0 }, 0091 { -3, 0, 0, 0 } }; 0092 0093 KSNumbers::KSNumbers(long double jd) 0094 { 0095 K.setD(20.49552 / 3600.); //set the constant of aberration 0096 0097 // ecliptic longitude of earth's perihelion, source: https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html; FIXME: We should correct this, as it changes with time. See the commit log for an order of magnitude estimate of the error. 0098 // FIXME: FIXME above seems to have been addressed? What is deltaEcLong? -- asimha 0099 P.setD(102.94719); 0100 0101 // JM 2017-09-04: Yes the longitude of the perihelion changes with time. The value below is the one for 04-09-2017 obtained from JPL Horizons system. 0102 // But some Earth orbital elements MUST be updated periodically like any other solar system body. 0103 // This is an important FIXME! 0104 // deltaEcLong is used in nutation calculations. Can P be obtained computationally? 0105 //P.setD(104.8089842092676); 0106 0107 computeConstantValues(); 0108 updateValues(jd); 0109 } 0110 0111 void KSNumbers::computeConstantValues() 0112 { 0113 // Compute those numbers that need to be computed only 0114 // once. 0115 // 0116 // Ideally, these should be computed at compile-time. When we 0117 // upgrade to C++11, we can make this function static and 0118 // constexpr (maybe?) But even otherwise, the overhead is very 0119 // negligible. 0120 0121 //Compute Precession Matrices from B1950 to 1984 using Newcomb formulae 0122 XB.setD(0.217697); 0123 YB.setD(0.189274); 0124 ZB.setD(0.217722); 0125 0126 XB.SinCos(SXB, CXB); 0127 YB.SinCos(SYB, CYB); 0128 ZB.SinCos(SZB, CZB); 0129 0130 //P1B is used to precess from 1984 to B1950: 0131 0132 P1B(0, 0) = CXB * CYB * CZB - SXB * SZB; 0133 P1B(1, 0) = CXB * CYB * SZB + SXB * CZB; 0134 P1B(2, 0) = CXB * SYB; 0135 P1B(0, 1) = -1.0 * SXB * CYB * CZB - CXB * SZB; 0136 P1B(1, 1) = -1.0 * SXB * CYB * SZB + CXB * CZB; 0137 P1B(2, 1) = -1.0 * SXB * SYB; 0138 P1B(0, 2) = -1.0 * SYB * CZB; 0139 P1B(1, 2) = -1.0 * SYB * SZB; 0140 P1B(2, 2) = CYB; 0141 0142 //P2 is used to precess from B1950 to 1984 (it is the transpose of P1) 0143 // FIXME: This can be optimized by taking the transpose of P1 instead of recomputing it from scratch 0144 P2B(0, 0) = CXB * CYB * CZB - SXB * SZB; 0145 P2B(1, 0) = -1.0 * SXB * CYB * CZB - CXB * SZB; 0146 P2B(2, 0) = -1.0 * SYB * CZB; 0147 P2B(0, 1) = CXB * CYB * SZB + SXB * CZB; 0148 P2B(1, 1) = -1.0 * SXB * CYB * SZB + CXB * CZB; 0149 P2B(2, 1) = -1.0 * SYB * SZB; 0150 P2B(0, 2) = CXB * SYB; 0151 P2B(1, 2) = -1.0 * SXB * SYB; 0152 P2B(2, 2) = CYB; 0153 } 0154 0155 void KSNumbers::updateValues(long double jd) 0156 { 0157 dms arg; 0158 double args, argc; 0159 0160 days = jd; 0161 0162 // FIXME: What is the source for these algorithms / polynomials / numbers? -- asimha 0163 0164 //Julian Centuries since J2000.0 0165 T = (jd - J2000) / 36525.; 0166 0167 // Julian Millenia since J2000.0 0168 jm = T / 10.0; 0169 0170 double T2 = T * T; 0171 double T3 = T2 * T; 0172 0173 //Sun's Mean Longitude 0174 L.setD(280.46645 + 36000.76983 * T + 0.0003032 * T2); 0175 0176 //Mean elongation of the Moon from the Sun 0177 D.setD(297.85036 + 445267.111480 * T - 0.0019142 * T2 + T3 / 189474.); 0178 0179 //Sun's Mean Anomaly 0180 M.setD(357.52910 + 35999.05030 * T - 0.0001559 * T2 - 0.00000048 * T3); 0181 0182 //Moon's Mean Anomaly 0183 MM.setD(134.96298 + 477198.867398 * T + 0.0086972 * T2 + T3 / 56250.0); 0184 0185 //Moon's Mean Longitude 0186 LM.setD(218.3164591 + 481267.88134236 * T - 0.0013268 * T2 + T3 / 538841. - T * T * T * T / 6519400.); 0187 0188 //Moon's argument of latitude 0189 F.setD(93.27191 + 483202.017538 * T - 0.0036825 * T2 + T3 / 327270.); 0190 0191 //Longitude of Moon's Ascending Node 0192 O.setD(125.04452 - 1934.136261 * T + 0.0020708 * T2 + T3 / 450000.0); 0193 0194 //Earth's orbital eccentricity 0195 e = 0.016708617 - 0.000042037 * T - 0.0000001236 * T2; 0196 0197 double C = (1.914600 - 0.004817 * T - 0.000014 * T2) * sin(M.radians()) + 0198 (0.019993 - 0.000101 * T) * sin(2.0 * M.radians()) + 0.000290 * sin(3.0 * M.radians()); 0199 0200 //Sun's True Longitude 0201 L0.setD(L.Degrees() + C); 0202 0203 //Sun's True Anomaly 0204 M0.setD(M.Degrees() + C); 0205 0206 //Obliquity of the Ecliptic 0207 // FIXME: This can be optimized by using the fact that U^3 = U^2 * U, so we can reduce the number of multiplications 0208 double U = T / 100.0; 0209 double dObliq = -4680.93 * U - 1.55 * U * U + 1999.25 * U * U * U - 51.38 * U * U * U * U - 0210 249.67 * U * U * U * U * U - 39.05 * U * U * U * U * U * U + 7.12 * U * U * U * U * U * U * U + 0211 27.87 * U * U * U * U * U * U * U * U + 5.79 * U * U * U * U * U * U * U * U * U + 0212 2.45 * U * U * U * U * U * U * U * U * U * U; 0213 Obliquity.setD(23.43929111 + dObliq / 3600.0); 0214 0215 //Nutation parameters 0216 dms L2, M2, O2; 0217 double sin2L, cos2L, sin2M, cos2M; 0218 double sinO, cosO, sin2O, cos2O; 0219 0220 O2.setD(2.0 * O.Degrees()); 0221 L2.setD(2.0 * L.Degrees()); //twice mean ecl. long. of Sun 0222 M2.setD(2.0 * LM.Degrees()); //twice mean ecl. long. of Moon 0223 0224 O.SinCos(sinO, cosO); 0225 O2.SinCos(sin2O, cos2O); 0226 L2.SinCos(sin2L, cos2L); 0227 M2.SinCos(sin2M, cos2M); 0228 0229 // deltaEcLong = ( -17.2*sinO - 1.32*sin2L - 0.23*sin2M + 0.21*sin2O)/3600.0; //Ecl. long. correction 0230 // deltaObliquity = ( 9.2*cosO + 0.57*cos2L + 0.10*cos2M - 0.09*cos2O)/3600.0; //Obliq. correction 0231 0232 deltaEcLong = 0.; 0233 deltaObliquity = 0.; 0234 0235 for (unsigned int i = 0; i < NUTTERMS; i++) 0236 { 0237 arg.setD(arguments[i][0] * D.Degrees() + arguments[i][1] * M.Degrees() + arguments[i][2] * MM.Degrees() + 0238 arguments[i][3] * F.Degrees() + arguments[i][4] * O.Degrees()); 0239 arg.SinCos(args, argc); 0240 0241 deltaEcLong += (amp[i][0] + amp[i][1] / 10. * T) * args * 1e-4; 0242 deltaObliquity += (amp[i][2] + amp[i][3] / 10. * T) * argc * 1e-4; 0243 } 0244 0245 deltaEcLong /= 3600.0; 0246 deltaObliquity /= 3600.0; 0247 0248 //Compute Precession Matrices: 0249 XP.setD(0.6406161 * T + 0.0000839 * T2 + 0.0000050 * T3); 0250 YP.setD(0.5567530 * T - 0.0001185 * T2 - 0.0000116 * T3); 0251 ZP.setD(0.6406161 * T + 0.0003041 * T2 + 0.0000051 * T3); 0252 0253 XP.SinCos(SX, CX); 0254 YP.SinCos(SY, CY); 0255 ZP.SinCos(SZ, CZ); 0256 0257 //P1 is used to precess from any epoch to J2000 0258 // Note: P1 is a rotation matrix, and P2 is its transpose = inverse (also a rotation matrix) 0259 P1(0, 0) = CX * CY * CZ - SX * SZ; 0260 P1(1, 0) = CX * CY * SZ + SX * CZ; 0261 P1(2, 0) = CX * SY; 0262 P1(0, 1) = -1.0 * SX * CY * CZ - CX * SZ; 0263 P1(1, 1) = -1.0 * SX * CY * SZ + CX * CZ; 0264 P1(2, 1) = -1.0 * SX * SY; 0265 P1(0, 2) = -1.0 * SY * CZ; 0266 P1(1, 2) = -1.0 * SY * SZ; 0267 P1(2, 2) = CY; 0268 0269 //P2 is used to precess from J2000 to any other epoch (it is the transpose of P1) 0270 // FIXME: More optimization -- just use P1(j, i) instead of P2(i, j) in code 0271 P2(0, 0) = P1(0, 0); 0272 P2(1, 0) = P1(0, 1); 0273 P2(2, 0) = P1(0, 2); 0274 P2(0, 1) = P1(1, 0); 0275 P2(1, 1) = P1(1, 1); 0276 P2(2, 1) = P1(1, 2); 0277 P2(0, 2) = P1(2, 0); 0278 P2(1, 2) = P1(2, 1); 0279 P2(2, 2) = P1(2, 2); 0280 0281 // Mean longitudes for the planets. radians 0282 // 0283 0284 // TODO Pasar a grados [Google Translate says "Jump to Degrees". --asimha] 0285 double LVenus = 3.1761467 + 1021.3285546 * T; // Venus 0286 double LMars = 1.7534703 + 628.3075849 * T; // Mars 0287 double LEarth = 6.2034809 + 334.0612431 * T; // Earth 0288 double LJupiter = 0.5995465 + 52.9690965 * T; // Jupiter 0289 double LSaturn = 0.8740168 + 21.3299095 * T; // Saturn 0290 double LNeptune = 5.3118863 + 3.8133036 * T; // Neptune 0291 double LUranus = 5.4812939 + 7.4781599 * T; // Uranus 0292 0293 double LMRad = 3.8103444 + 8399.6847337 * T; // Moon 0294 double DRad = 5.1984667 + 7771.3771486 * T; 0295 double MMRad = 2.3555559 + 8328.6914289 * T; // Moon 0296 double FRad = 1.6279052 + 8433.4661601 * T; 0297 0298 /** Contributions to the velocity of the Earth referred to the barycenter of the solar system 0299 in the J2000 equatorial system 0300 Velocities 10^{-8} AU/day 0301 Ron & Vondrak method 0302 **/ 0303 0304 double vondrak[36][7] = { 0305 { LMars, -1719914 - 2 * T, -25, 25 - 13 * T, 1578089 + 156 * T, 10 + 32 * T, 684185 - 358 * T }, 0306 { 2 * LMars, 6434 + 141 * T, 28007 - 107 * T, 25697 - 95 * T, -5904 - 130 * T, 11141 - 48 * T, -2559 - 55 * T }, 0307 { LJupiter, 715, 0, 6, -657, -15, -282 }, 0308 { LMRad, 715, 0, 0, -656, 0, -285 }, 0309 { 3 * LMars, 486 - 5 * T, -236 - 4 * T, -216 - 4 * T, -446 + 5 * T, -94, -193 }, 0310 { LSaturn, 159, 0, 2, -147, -6, -61 }, 0311 { FRad, 0, 0, 0, 26, 0, -59 }, 0312 { LMRad + MMRad, 39, 0, 0, -36, 0, -16 }, 0313 { 2 * LJupiter, 33, -10, -9, -30, -5, -13 }, 0314 { 2 * LMars - LJupiter, 31, 1, 1, -28, 0, -12 }, 0315 { 3 * LMars - 8 * LEarth + 3 * LJupiter, 8, -28, 25, 8, 11, 3 }, 0316 { 5 * LMars - 8 * LEarth + 3 * LJupiter, 8, -28, -25, -8, -11, -3 }, 0317 { 2 * LVenus - LMars, 21, 0, 0, -19, 0, -8 }, 0318 { LVenus, -19, 0, 0, 17, 0, 8 }, 0319 { LNeptune, 17, 0, 0, -16, 0, -7 }, 0320 { LMars - 2 * LJupiter, 16, 0, 0, 15, 1, 7 }, 0321 { LUranus, 16, 0, 1, -15, -3, -6 }, 0322 { LMars + LJupiter, 11, -1, -1, -10, -1, -5 }, 0323 { 2 * LVenus - 2 * LMars, 0, -11, -10, 0, -4, 0 }, 0324 { LMars - LJupiter, -11, -2, -2, 9, -1, 4 }, 0325 { 4 * LMars, -7, -8, -8, 6, -3, 3 }, 0326 { 3 * LMars - 2 * LJupiter, -10, 0, 0, 9, 0, 4 }, 0327 { LVenus - 2 * LMars, -9, 0, 0, -9, 0, -4 }, 0328 { 2 * LVenus - 3 * LMars, -9, 0, 0, -8, 0, -4 }, 0329 { 2 * LSaturn, 0, -9, -8, 0, -3, 0 }, 0330 { 2 * LVenus - 4 * LMars, 0, -9, 8, 0, 3, 0 }, 0331 { 3 * LMars - 2 * LEarth, 8, 0, 0, -8, 0, -3 }, 0332 { LMRad + 2 * DRad - MMRad, 8, 0, 0, -7, 0, -3 }, 0333 { 8 * LVenus - 12 * LMars, -4, -7, -6, 4, -3, 2 }, 0334 { 8 * LVenus - 14 * LMars, -4, -7, 6, -4, 3, -2 }, 0335 { 2 * LEarth, -6, -5, -4, 5, -2, 2 }, 0336 { 3 * LVenus - 4 * LMars, -1, -1, -2, -7, 1, -4 }, 0337 { 2 * LMars - 2 * LJupiter, 4, -6, -5, -4, -2, -2 }, 0338 { 3 * LVenus - 3 * LMars, 0, -7, -6, 0, -3, 0 }, 0339 { 2 * LMars - 2 * LEarth, 5, -5, -4, -5, -2, -2 }, 0340 { LMRad - 2 * DRad, 5, 0, 0, -5, 0, -2 } 0341 }; 0342 0343 dms anglev; 0344 double sa, ca; 0345 // Vearth X component 0346 vearth[0] = 0.; 0347 // Vearth Y component 0348 vearth[1] = 0.; 0349 // Vearth Z component 0350 vearth[2] = 0.; 0351 0352 for (auto &item : vondrak) 0353 { 0354 anglev.setRadians(item[0]); 0355 anglev.SinCos(sa, ca); 0356 for (unsigned int j = 0; j < 3; j++) 0357 { 0358 vearth[j] += item[2 * j + 1] * sa + item[2 * j + 2] * ca; 0359 } 0360 } 0361 0362 const double UA2km = 1.49597870 / 86400.; // 10^{-8}*UA/dia -> km/s 0363 0364 for (double &item : vearth) 0365 { 0366 item *= UA2km; 0367 } 0368 }