Kstars

ksnumbers.cpp
1/*
2 SPDX-FileCopyrightText: 2002-2005 Jason Harris <kstars@30doradus.org>
3 SPDX-FileCopyrightText: 2004-2005 Pablo de Vicente <p.devicente@wanadoo.es>
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
13const 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
29const 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
93KSNumbers::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
107 computeConstantValues();
108 updateValues(jd);
109}
110
111void KSNumbers::computeConstantValues()
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
155void 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}
CachingDms::setD
void setD(const double &x) override
Sets the angle in degrees supplied as a double.
Definition cachingdms.h:59
KSNumbers::updateValues
void updateValues(long double jd)
update all values for the date given as an argument.
Definition ksnumbers.cpp:155
KSNumbers::dObliq
double dObliq() const
Definition ksnumbers.h:81
KSNumbers::KSNumbers
KSNumbers(long double jd)
Constructor.
Definition ksnumbers.cpp:93
KSNumbers::computeConstantValues
void computeConstantValues()
compute constant values that need to be computed only once per instance of the application
Definition ksnumbers.cpp:111
dms
An angle, stored as degrees, but expressible in many ways.
Definition dms.h:38
dms::SinCos
void SinCos(double &s, double &c) const
Compute Sine and Cosine of the angle simultaneously.
Definition dms.h:447
dms::radians
double radians() const
Express the angle in radians.
Definition dms.h:325
dms::setRadians
virtual void setRadians(const double &Rad)
Set angle according to the argument, in radians.
Definition dms.h:333
dms::setD
virtual void setD(const double &x)
Sets floating-point value of angle, in degrees.
Definition dms.h:179
dms::Degrees
const double & Degrees() const
Definition dms.h:141
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Mon Nov 4 2024 16:38:43 by doxygen 1.12.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.