File indexing completed on 2022-12-06 18:58:43

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 }