Kstars

ksnumbers.cpp
1 /*
2  SPDX-FileCopyrightText: 2002-2005 Jason Harris <[email protected]>
3  SPDX-FileCopyrightText: 2004-2005 Pablo de Vicente <[email protected]>
4 
5  SPDX-License-Identifier: GPL-2.0-or-later
6 */
7 
8 #include "ksnumbers.h"
9 
10 #include "kstarsdatetime.h" //for J2000 define
11 
12 // 63 elements
13 const int KSNumbers::arguments[NUTTERMS][5] = {
14  { 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 },
15  { 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 },
16  { -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 },
17  { 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 },
18  { 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 },
19  { -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 },
20  { 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 },
21  { 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 },
22  { 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 },
23  { 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 },
24  { -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 },
25  { 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 },
26  { 2, -1, -1, 2, 2 }, { 0, 0, 3, 2, 2 }, { 2, -1, 0, 2, 2 }
27 };
28 
29 const int KSNumbers::amp[NUTTERMS][4] = { { -171996, -1742, 92025, 89 },
30  { -13187, -16, 5736, -31 },
31  { -2274, -2, 977, -5 },
32  { 2062, 2, -895, 5 },
33  { 1426, -34, 54, -1 },
34  { 712, 1, -7, 0 },
35  { -517, 12, 224, -6 },
36  { -386, -4, 200, 0 },
37  { -301, 0, 129, -1 },
38  { 217, -5, -95, 3 },
39  { -158, 0, 0, 0 },
40  { 129, 1, -70, 0 },
41  { 123, 0, -53, 0 },
42  { 63, 0, 0, 0 },
43  { 63, 1, -33, 0 },
44  { -59, 0, 26, 0 },
45  { -58, -1, 32, 0 },
46  { -51, 0, 27, 0 },
47  { 48, 0, 0, 0 },
48  { 46, 0, -24, 0 },
49  { -38, 0, 16, 0 },
50  { -31, 0, 13, 0 },
51  { 29, 0, 0, 0 },
52  { 29, 0, -12, 0 },
53  { 26, 0, 0, 0 },
54  { -22, 0, 0, 0 },
55  { 21, 0, -10, 0 },
56  { 17, -1, 0, 0 },
57  { 16, 0, -8, 0 },
58  { -16, 1, 7, 0 },
59  { -15, 0, 9, 0 },
60  { -13, 0, 7, 0 },
61  { -12, 0, 6, 0 },
62  { 11, 0, 0, 0 },
63  { -10, 0, 5, 0 },
64  { -8, 0, 3, 0 },
65  { 7, 0, -3, 0 },
66  { -7, 0, 0, 0 },
67  { -7, 0, 3, 0 },
68  { -7, 0, 3, 0 },
69  { 6, 0, 0, 0 },
70  { 6, 0, -3, 0 },
71  { 6, 0, -3, 0 },
72  { -6, 0, 3, 0 },
73  { -6, 0, 3, 0 },
74  { 5, 0, 0, 0 },
75  { -5, 0, 3, 0 },
76  { -5, 0, 3, 0 },
77  { -5, 0, 3, 0 },
78  { 4, 0, 0, 0 },
79  { 4, 0, 0, 0 },
80  { 4, 0, 0, 0 },
81  { -4, 0, 0, 0 },
82  { -4, 0, 0, 0 },
83  { -4, 0, 0, 0 },
84  { 3, 0, 0, 0 },
85  { -3, 0, 0, 0 },
86  { -3, 0, 0, 0 },
87  { -3, 0, 0, 0 },
88  { -3, 0, 0, 0 },
89  { -3, 0, 0, 0 },
90  { -3, 0, 0, 0 },
91  { -3, 0, 0, 0 } };
92 
93 KSNumbers::KSNumbers(long double jd)
94 {
95  K.setD(20.49552 / 3600.); //set the constant of aberration
96 
97  // 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.
98  // FIXME: FIXME above seems to have been addressed? What is deltaEcLong? -- asimha
99  P.setD(102.94719);
100 
101  // 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.
102  // But some Earth orbital elements MUST be updated periodically like any other solar system body.
103  // This is an important FIXME!
104  // deltaEcLong is used in nutation calculations. Can P be obtained computationally?
105  //P.setD(104.8089842092676);
106 
108  updateValues(jd);
109 }
110 
112 {
113  // Compute those numbers that need to be computed only
114  // once.
115  //
116  // Ideally, these should be computed at compile-time. When we
117  // upgrade to C++11, we can make this function static and
118  // constexpr (maybe?) But even otherwise, the overhead is very
119  // negligible.
120 
121  //Compute Precession Matrices from B1950 to 1984 using Newcomb formulae
122  XB.setD(0.217697);
123  YB.setD(0.189274);
124  ZB.setD(0.217722);
125 
126  XB.SinCos(SXB, CXB);
127  YB.SinCos(SYB, CYB);
128  ZB.SinCos(SZB, CZB);
129 
130  //P1B is used to precess from 1984 to B1950:
131 
132  P1B(0, 0) = CXB * CYB * CZB - SXB * SZB;
133  P1B(1, 0) = CXB * CYB * SZB + SXB * CZB;
134  P1B(2, 0) = CXB * SYB;
135  P1B(0, 1) = -1.0 * SXB * CYB * CZB - CXB * SZB;
136  P1B(1, 1) = -1.0 * SXB * CYB * SZB + CXB * CZB;
137  P1B(2, 1) = -1.0 * SXB * SYB;
138  P1B(0, 2) = -1.0 * SYB * CZB;
139  P1B(1, 2) = -1.0 * SYB * SZB;
140  P1B(2, 2) = CYB;
141 
142  //P2 is used to precess from B1950 to 1984 (it is the transpose of P1)
143  // FIXME: This can be optimized by taking the transpose of P1 instead of recomputing it from scratch
144  P2B(0, 0) = CXB * CYB * CZB - SXB * SZB;
145  P2B(1, 0) = -1.0 * SXB * CYB * CZB - CXB * SZB;
146  P2B(2, 0) = -1.0 * SYB * CZB;
147  P2B(0, 1) = CXB * CYB * SZB + SXB * CZB;
148  P2B(1, 1) = -1.0 * SXB * CYB * SZB + CXB * CZB;
149  P2B(2, 1) = -1.0 * SYB * SZB;
150  P2B(0, 2) = CXB * SYB;
151  P2B(1, 2) = -1.0 * SXB * SYB;
152  P2B(2, 2) = CYB;
153 }
154 
155 void KSNumbers::updateValues(long double jd)
156 {
157  dms arg;
158  double args, argc;
159 
160  days = jd;
161 
162  // FIXME: What is the source for these algorithms / polynomials / numbers? -- asimha
163 
164  //Julian Centuries since J2000.0
165  T = (jd - J2000) / 36525.;
166 
167  // Julian Millenia since J2000.0
168  jm = T / 10.0;
169 
170  double T2 = T * T;
171  double T3 = T2 * T;
172 
173  //Sun's Mean Longitude
174  L.setD(280.46645 + 36000.76983 * T + 0.0003032 * T2);
175 
176  //Mean elongation of the Moon from the Sun
177  D.setD(297.85036 + 445267.111480 * T - 0.0019142 * T2 + T3 / 189474.);
178 
179  //Sun's Mean Anomaly
180  M.setD(357.52910 + 35999.05030 * T - 0.0001559 * T2 - 0.00000048 * T3);
181 
182  //Moon's Mean Anomaly
183  MM.setD(134.96298 + 477198.867398 * T + 0.0086972 * T2 + T3 / 56250.0);
184 
185  //Moon's Mean Longitude
186  LM.setD(218.3164591 + 481267.88134236 * T - 0.0013268 * T2 + T3 / 538841. - T * T * T * T / 6519400.);
187 
188  //Moon's argument of latitude
189  F.setD(93.27191 + 483202.017538 * T - 0.0036825 * T2 + T3 / 327270.);
190 
191  //Longitude of Moon's Ascending Node
192  O.setD(125.04452 - 1934.136261 * T + 0.0020708 * T2 + T3 / 450000.0);
193 
194  //Earth's orbital eccentricity
195  e = 0.016708617 - 0.000042037 * T - 0.0000001236 * T2;
196 
197  double C = (1.914600 - 0.004817 * T - 0.000014 * T2) * sin(M.radians()) +
198  (0.019993 - 0.000101 * T) * sin(2.0 * M.radians()) + 0.000290 * sin(3.0 * M.radians());
199 
200  //Sun's True Longitude
201  L0.setD(L.Degrees() + C);
202 
203  //Sun's True Anomaly
204  M0.setD(M.Degrees() + C);
205 
206  //Obliquity of the Ecliptic
207  // FIXME: This can be optimized by using the fact that U^3 = U^2 * U, so we can reduce the number of multiplications
208  double U = T / 100.0;
209  double dObliq = -4680.93 * U - 1.55 * U * U + 1999.25 * U * U * U - 51.38 * U * U * U * U -
210  249.67 * U * U * U * U * U - 39.05 * U * U * U * U * U * U + 7.12 * U * U * U * U * U * U * U +
211  27.87 * U * U * U * U * U * U * U * U + 5.79 * U * U * U * U * U * U * U * U * U +
212  2.45 * U * U * U * U * U * U * U * U * U * U;
213  Obliquity.setD(23.43929111 + dObliq / 3600.0);
214 
215  //Nutation parameters
216  dms L2, M2, O2;
217  double sin2L, cos2L, sin2M, cos2M;
218  double sinO, cosO, sin2O, cos2O;
219 
220  O2.setD(2.0 * O.Degrees());
221  L2.setD(2.0 * L.Degrees()); //twice mean ecl. long. of Sun
222  M2.setD(2.0 * LM.Degrees()); //twice mean ecl. long. of Moon
223 
224  O.SinCos(sinO, cosO);
225  O2.SinCos(sin2O, cos2O);
226  L2.SinCos(sin2L, cos2L);
227  M2.SinCos(sin2M, cos2M);
228 
229  // deltaEcLong = ( -17.2*sinO - 1.32*sin2L - 0.23*sin2M + 0.21*sin2O)/3600.0; //Ecl. long. correction
230  // deltaObliquity = ( 9.2*cosO + 0.57*cos2L + 0.10*cos2M - 0.09*cos2O)/3600.0; //Obliq. correction
231 
232  deltaEcLong = 0.;
233  deltaObliquity = 0.;
234 
235  for (unsigned int i = 0; i < NUTTERMS; i++)
236  {
237  arg.setD(arguments[i][0] * D.Degrees() + arguments[i][1] * M.Degrees() + arguments[i][2] * MM.Degrees() +
238  arguments[i][3] * F.Degrees() + arguments[i][4] * O.Degrees());
239  arg.SinCos(args, argc);
240 
241  deltaEcLong += (amp[i][0] + amp[i][1] / 10. * T) * args * 1e-4;
242  deltaObliquity += (amp[i][2] + amp[i][3] / 10. * T) * argc * 1e-4;
243  }
244 
245  deltaEcLong /= 3600.0;
246  deltaObliquity /= 3600.0;
247 
248  //Compute Precession Matrices:
249  XP.setD(0.6406161 * T + 0.0000839 * T2 + 0.0000050 * T3);
250  YP.setD(0.5567530 * T - 0.0001185 * T2 - 0.0000116 * T3);
251  ZP.setD(0.6406161 * T + 0.0003041 * T2 + 0.0000051 * T3);
252 
253  XP.SinCos(SX, CX);
254  YP.SinCos(SY, CY);
255  ZP.SinCos(SZ, CZ);
256 
257  //P1 is used to precess from any epoch to J2000
258  // Note: P1 is a rotation matrix, and P2 is its transpose = inverse (also a rotation matrix)
259  P1(0, 0) = CX * CY * CZ - SX * SZ;
260  P1(1, 0) = CX * CY * SZ + SX * CZ;
261  P1(2, 0) = CX * SY;
262  P1(0, 1) = -1.0 * SX * CY * CZ - CX * SZ;
263  P1(1, 1) = -1.0 * SX * CY * SZ + CX * CZ;
264  P1(2, 1) = -1.0 * SX * SY;
265  P1(0, 2) = -1.0 * SY * CZ;
266  P1(1, 2) = -1.0 * SY * SZ;
267  P1(2, 2) = CY;
268 
269  //P2 is used to precess from J2000 to any other epoch (it is the transpose of P1)
270  // FIXME: More optimization -- just use P1(j, i) instead of P2(i, j) in code
271  P2(0, 0) = P1(0, 0);
272  P2(1, 0) = P1(0, 1);
273  P2(2, 0) = P1(0, 2);
274  P2(0, 1) = P1(1, 0);
275  P2(1, 1) = P1(1, 1);
276  P2(2, 1) = P1(1, 2);
277  P2(0, 2) = P1(2, 0);
278  P2(1, 2) = P1(2, 1);
279  P2(2, 2) = P1(2, 2);
280 
281  // Mean longitudes for the planets. radians
282  //
283 
284  // TODO Pasar a grados [Google Translate says "Jump to Degrees". --asimha]
285  double LVenus = 3.1761467 + 1021.3285546 * T; // Venus
286  double LMars = 1.7534703 + 628.3075849 * T; // Mars
287  double LEarth = 6.2034809 + 334.0612431 * T; // Earth
288  double LJupiter = 0.5995465 + 52.9690965 * T; // Jupiter
289  double LSaturn = 0.8740168 + 21.3299095 * T; // Saturn
290  double LNeptune = 5.3118863 + 3.8133036 * T; // Neptune
291  double LUranus = 5.4812939 + 7.4781599 * T; // Uranus
292 
293  double LMRad = 3.8103444 + 8399.6847337 * T; // Moon
294  double DRad = 5.1984667 + 7771.3771486 * T;
295  double MMRad = 2.3555559 + 8328.6914289 * T; // Moon
296  double FRad = 1.6279052 + 8433.4661601 * T;
297 
298  /** Contributions to the velocity of the Earth referred to the barycenter of the solar system
299  in the J2000 equatorial system
300  Velocities 10^{-8} AU/day
301  Ron & Vondrak method
302  **/
303 
304  double vondrak[36][7] = {
305  { LMars, -1719914 - 2 * T, -25, 25 - 13 * T, 1578089 + 156 * T, 10 + 32 * T, 684185 - 358 * T },
306  { 2 * LMars, 6434 + 141 * T, 28007 - 107 * T, 25697 - 95 * T, -5904 - 130 * T, 11141 - 48 * T, -2559 - 55 * T },
307  { LJupiter, 715, 0, 6, -657, -15, -282 },
308  { LMRad, 715, 0, 0, -656, 0, -285 },
309  { 3 * LMars, 486 - 5 * T, -236 - 4 * T, -216 - 4 * T, -446 + 5 * T, -94, -193 },
310  { LSaturn, 159, 0, 2, -147, -6, -61 },
311  { FRad, 0, 0, 0, 26, 0, -59 },
312  { LMRad + MMRad, 39, 0, 0, -36, 0, -16 },
313  { 2 * LJupiter, 33, -10, -9, -30, -5, -13 },
314  { 2 * LMars - LJupiter, 31, 1, 1, -28, 0, -12 },
315  { 3 * LMars - 8 * LEarth + 3 * LJupiter, 8, -28, 25, 8, 11, 3 },
316  { 5 * LMars - 8 * LEarth + 3 * LJupiter, 8, -28, -25, -8, -11, -3 },
317  { 2 * LVenus - LMars, 21, 0, 0, -19, 0, -8 },
318  { LVenus, -19, 0, 0, 17, 0, 8 },
319  { LNeptune, 17, 0, 0, -16, 0, -7 },
320  { LMars - 2 * LJupiter, 16, 0, 0, 15, 1, 7 },
321  { LUranus, 16, 0, 1, -15, -3, -6 },
322  { LMars + LJupiter, 11, -1, -1, -10, -1, -5 },
323  { 2 * LVenus - 2 * LMars, 0, -11, -10, 0, -4, 0 },
324  { LMars - LJupiter, -11, -2, -2, 9, -1, 4 },
325  { 4 * LMars, -7, -8, -8, 6, -3, 3 },
326  { 3 * LMars - 2 * LJupiter, -10, 0, 0, 9, 0, 4 },
327  { LVenus - 2 * LMars, -9, 0, 0, -9, 0, -4 },
328  { 2 * LVenus - 3 * LMars, -9, 0, 0, -8, 0, -4 },
329  { 2 * LSaturn, 0, -9, -8, 0, -3, 0 },
330  { 2 * LVenus - 4 * LMars, 0, -9, 8, 0, 3, 0 },
331  { 3 * LMars - 2 * LEarth, 8, 0, 0, -8, 0, -3 },
332  { LMRad + 2 * DRad - MMRad, 8, 0, 0, -7, 0, -3 },
333  { 8 * LVenus - 12 * LMars, -4, -7, -6, 4, -3, 2 },
334  { 8 * LVenus - 14 * LMars, -4, -7, 6, -4, 3, -2 },
335  { 2 * LEarth, -6, -5, -4, 5, -2, 2 },
336  { 3 * LVenus - 4 * LMars, -1, -1, -2, -7, 1, -4 },
337  { 2 * LMars - 2 * LJupiter, 4, -6, -5, -4, -2, -2 },
338  { 3 * LVenus - 3 * LMars, 0, -7, -6, 0, -3, 0 },
339  { 2 * LMars - 2 * LEarth, 5, -5, -4, -5, -2, -2 },
340  { LMRad - 2 * DRad, 5, 0, 0, -5, 0, -2 }
341  };
342 
343  dms anglev;
344  double sa, ca;
345  // Vearth X component
346  vearth[0] = 0.;
347  // Vearth Y component
348  vearth[1] = 0.;
349  // Vearth Z component
350  vearth[2] = 0.;
351 
352  for (auto &item : vondrak)
353  {
354  anglev.setRadians(item[0]);
355  anglev.SinCos(sa, ca);
356  for (unsigned int j = 0; j < 3; j++)
357  {
358  vearth[j] += item[2 * j + 1] * sa + item[2 * j + 2] * ca;
359  }
360  }
361 
362  const double UA2km = 1.49597870 / 86400.; // 10^{-8}*UA/dia -> km/s
363 
364  for (double &item : vearth)
365  {
366  item *= UA2km;
367  }
368 }
void setD(const double &x) override
Sets the angle in degrees supplied as a double.
Definition: cachingdms.h:59
virtual void setD(const double &x)
Sets floating-point value of angle, in degrees.
Definition: dms.h:179
void SinCos(double &s, double &c) const
Compute Sine and Cosine of the angle simultaneously.
Definition: dms.h:439
void computeConstantValues()
compute constant values that need to be computed only once per instance of the application
Definition: ksnumbers.cpp:111
KSNumbers(long double jd)
Constructor.
Definition: ksnumbers.cpp:93
void updateValues(long double jd)
update all values for the date given as an argument.
Definition: ksnumbers.cpp:155
double dObliq() const
Definition: ksnumbers.h:81
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:328
double radians() const
Express the angle in radians.
Definition: dms.h:320
const double & Degrees() const
Definition: dms.h:141
This file is part of the KDE documentation.
Documentation copyright © 1996-2022 The KDE developers.
Generated on Fri Aug 12 2022 04:00:55 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.