?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
moved kdeaccessibility kdeaddons kdeadmin kdeartwork kdebindings kdeedu kdegames...
[kdeedu.git] / kstars / kstars / ksnumbers.cpp
blob1e4e8e9c32e502d3f4e39d5d2660fc203e4db0be
1 /***************************************************************************
2 ksnumbers.cpp - description
3 -------------------
4 begin : Sun Jan 13 2002
5 copyright : (C) 2002-2005 by Jason Harris
6 email : kstars@30doradus.org
7 copyright : (C) 2004-2005 by Pablo de Vicente
8 email : p.devicente@wanadoo.es
9 ***************************************************************************/
10
11 /***************************************************************************
12 * *
13 * This program is free software; you can redistribute it and/or modify *
14 * it under the terms of the GNU General Public License as published by *
15 * the Free Software Foundation; either version 2 of the License, or *
16 * (at your option) any later version. *
17 * *
18 ***************************************************************************/
19
20 #include "ksnumbers.h"
21
22 // 63 elements
23 const int KSNumbers::arguments[NUTTERMS][5] = {
24 { 0, 0, 0, 0, 1},
25 {-2, 0, 0, 2, 2},
26 { 0, 0, 0, 2, 2},
27 { 0, 0, 0, 0, 2},
28 { 0, 1, 0, 0, 0},
29 { 0, 0, 1, 0, 0},
30 {-2, 1, 0, 2, 2},
31 { 0, 0, 0, 2, 1},
32 { 0, 0, 1, 2, 2},
33 {-2,-1, 0, 2, 2},
34 {-2, 0, 1, 0, 0},
35 {-2, 0, 0, 2, 1},
36 { 0, 0,-1, 2, 2},
37 { 2, 0, 0, 0, 0},
38 { 0, 0, 1, 0, 1},
39 { 2, 0,-1, 2, 2},
40 { 0, 0,-1, 0, 1},
41 { 0, 0, 1, 2, 1},
42 {-2, 0, 2, 0, 0},
43 { 0, 0,-2, 2, 1},
44 { 2, 0, 0, 2, 2},
45 { 0, 0, 2, 2, 2},
46 { 0, 0, 2, 0, 0},
47 {-2, 0, 1, 2, 2},
48 { 0, 0, 0, 2, 0},
49 {-2, 0, 0, 2, 0},
50 { 0, 0,-1, 2, 1},
51 { 0, 2, 0, 0, 0},
52 { 2, 0,-1, 0, 1},
53 {-2, 2, 0, 2, 2},
54 { 0, 1, 0, 0, 1},
55 {-2, 0, 1, 0, 1},
56 { 0,-1, 0, 0, 1},
57 { 0, 0, 2,-2, 0},
58 { 2, 0,-1, 2, 1},
59 { 2, 0, 1, 2, 2},
60 { 0, 1, 0, 2, 2},
61 {-2, 1, 1, 0, 0},
62 { 0,-1, 0, 2, 2},
63 { 2, 0, 0, 2, 1},
64 { 2, 0, 1, 0, 0},
65 {-2, 0, 2, 2, 2},
66 {-2, 0, 1, 2, 1},
67 { 2, 0,-2, 0, 1},
68 { 2, 0, 0, 0, 1},
69 { 0,-1, 1, 0, 0},
70 {-2,-1, 0, 2, 1},
71 {-2, 0, 0, 0, 1},
72 { 0, 0, 2, 2, 1},
73 {-2, 0, 2, 0, 1},
74 {-2, 1, 0, 2, 1},
75 { 0, 0, 1,-2, 0},
76 {-1, 0, 1, 0, 0},
77 {-2, 1, 0, 0, 0},
78 { 1, 0, 0, 0, 0},
79 { 0, 0, 1, 2, 0},
80 { 0, 0,-2, 2, 2},
81 {-1,-1, 1, 0, 0},
82 { 0, 1, 1, 0, 0},
83 { 0,-1, 1, 2, 2},
84 { 2,-1,-1, 2, 2},
85 { 0, 0, 3, 2, 2},
86 { 2,-1, 0, 2, 2}
87 };
88
89 const int KSNumbers::amp[NUTTERMS][4] = {
90 {-171996,-1742, 92025, 89},
91 { -13187, -16, 5736,-31},
92 { -2274, -2, 977, -5},
93 { 2062, 2, -895, 5},
94 { 1426, -34, 54, -1},
95 { 712, 1, -7, 0},
96 { -517, 12, 224, -6},
97 { -386, -4, 200, 0},
98 { -301, 0, 129, -1},
99 { 217, -5, -95, 3},
100 { -158, 0, 0, 0},
101 { 129, 1, -70, 0},
102 { 123, 0, -53, 0},
103 { 63, 0, 0, 0},
104 { 63, 1, -33, 0},
105 { -59, 0, 26, 0},
106 { -58, -1, 32, 0},
107 { -51, 0, 27, 0},
108 { 48, 0, 0, 0},
109 { 46, 0, -24, 0},
110 { -38, 0, 16, 0},
111 { -31, 0, 13, 0},
112 { 29, 0, 0, 0},
113 { 29, 0, -12, 0},
114 { 26, 0, 0, 0},
115 { -22, 0, 0, 0},
116 { 21, 0, -10, 0},
117 { 17, -1, 0, 0},
118 { 16, 0, -8, 0},
119 { -16, 1, 7, 0},
120 { -15, 0, 9, 0},
121 { -13, 0, 7, 0},
122 { -12, 0, 6, 0},
123 { 11, 0, 0, 0},
124 { -10, 0, 5, 0},
125 { -8, 0, 3, 0},
126 { 7, 0, -3, 0},
127 { -7, 0, 0, 0},
128 { -7, 0, 3, 0},
129 { -7, 0, 3, 0},
130 { 6, 0, 0, 0},
131 { 6, 0, -3, 0},
132 { 6, 0, -3, 0},
133 { -6, 0, 3, 0},
134 { -6, 0, 3, 0},
135 { 5, 0, 0, 0},
136 { -5, 0, 3, 0},
137 { -5, 0, 3, 0},
138 { -5, 0, 3, 0},
139 { 4, 0, 0, 0},
140 { 4, 0, 0, 0},
141 { 4, 0, 0, 0},
142 { -4, 0, 0, 0},
143 { -4, 0, 0, 0},
144 { -4, 0, 0, 0},
145 { 3, 0, 0, 0},
146 { -3, 0, 0, 0},
147 { -3, 0, 0, 0},
148 { -3, 0, 0, 0},
149 { -3, 0, 0, 0},
150 { -3, 0, 0, 0},
151 { -3, 0, 0, 0},
152 { -3, 0, 0, 0}
153 };
154
155
156 KSNumbers::KSNumbers( long double jd ){
157 K.setD( 20.49552 / 3600. ); //set the constant of aberration
158 updateValues( jd );
159 }
160
161 KSNumbers::~KSNumbers(){
162 }
163
164 void KSNumbers::updateValues( long double jd ) {
165 dms arg;
166 double args, argc;
167
168 days = jd;
169
170 //Julian Centuries since J2000.0
171 T = ( jd - J2000 ) / 36525.;
172
173 // Julian Millenia since J2000.0
174 jm = T / 10.0;
175
176 double T2 = T*T;
177 double T3 = T2*T;
178
179 //Sun's Mean Longitude
180 L.setD( 280.46645 + 36000.76983*T + 0.0003032*T2 );
181
182 //Mean elongation of the Moon from the Sun
183 D.setD( 297.85036 + 445267.111480*T - 0.0019142*T2 + T3/189474.);
184
185 //Sun's Mean Anomaly
186 M.setD( 357.52910 + 35999.05030*T - 0.0001559*T2 - 0.00000048*T3);
187
188 //Moon's Mean Anomaly
189 MM.setD( 134.96298 + 477198.867398*T + 0.0086972*T2 + T3/56250.0 );
190
191 //Moon's Mean Longitude
192 LM.setD( 218.3164591 + 481267.88134236*T - 0.0013268*T2 + T3/538841. - T*T*T*T/6519400.);
193
194 //Moon's argument of latitude
195 F.setD( 93.27191 + 483202.017538*T - 0.0036825*T2 + T3/327270.);
196
197 //Longitude of Moon's Ascending Node
198 O.setD( 125.04452 - 1934.136261*T + 0.0020708*T2 + T3/450000.0 );
199
200 //Earth's orbital eccentricity
201 e = 0.016708617 - 0.000042037*T - 0.0000001236*T2;
202
203 double C = ( 1.914600 - 0.004817*T - 0.000014*T2 ) * sin( M.radians() )
204 + ( 0.019993 - 0.000101*T ) * sin( 2.0* M.radians() )
205 + 0.000290 * sin( 3.0* M.radians() );
206
207 //Sun's True Longitude
208 L0.setD( L.Degrees() + C );
209
210 //Sun's True Anomaly
211 M0.setD( M.Degrees() + C );
212
213 //Obliquity of the Ecliptic
214 double U = T/100.0;
215 double dObliq = -4680.93*U - 1.55*U*U + 1999.25*U*U*U
216 - 51.38*U*U*U*U - 249.67*U*U*U*U*U
217 - 39.05*U*U*U*U*U*U + 7.12*U*U*U*U*U*U*U
218 + 27.87*U*U*U*U*U*U*U*U + 5.79*U*U*U*U*U*U*U*U*U
219 + 2.45*U*U*U*U*U*U*U*U*U*U;
220 Obliquity.setD( 23.43929111 + dObliq/3600.0);
221
222 //Nutation parameters
223 dms L2, M2, O2;
224 double sin2L, cos2L, sin2M, cos2M;
225 double sinO, cosO, sin2O, cos2O;
226
227 O2.setD( 2.0*O.Degrees() );
228 L2.setD( 2.0*L.Degrees() ); //twice mean ecl. long. of Sun
229 M2.setD( 2.0*LM.Degrees() ); //twice mean ecl. long. of Moon
230
231 O.SinCos( sinO, cosO );
232 O2.SinCos( sin2O, cos2O );
233 L2.SinCos( sin2L, cos2L );
234 M2.SinCos( sin2M, cos2M );
235
236 // deltaEcLong = ( -17.2*sinO - 1.32*sin2L - 0.23*sin2M + 0.21*sin2O)/3600.0; //Ecl. long. correction
237 // deltaObliquity = ( 9.2*cosO + 0.57*cos2L + 0.10*cos2M - 0.09*cos2O)/3600.0; //Obliq. correction
238
239 deltaEcLong = 0.;
240 deltaObliquity = 0.;
241
242 for (unsigned int i=0; i < NUTTERMS; i++) {
243
244 arg.setD ( arguments[i][0]*D.Degrees() + arguments[i][1]*M.Degrees() +
245 arguments[i][2]*MM.Degrees() + arguments[i][3]*F.Degrees() + arguments[i][4]*O.Degrees() );
246 arg.SinCos( args, argc );
247
248 deltaEcLong += (amp[i][0] + amp[i][1]/10. * T ) * args * 1e-4 ;
249 deltaObliquity += (amp[i][2] + amp[i][3]/10. * T ) * argc * 1e-4 ;
250 }
251
252 deltaEcLong/= 3600.0;
253 deltaObliquity /= 3600.0;
254
255 //Compute Precession Matrices:
256 XP.setD( 0.6406161*T + 0.0000839*T2 + 0.0000050*T3 );
257 YP.setD( 0.5567530*T - 0.0001185*T2 - 0.0000116*T3 );
258 ZP.setD( 0.6406161*T + 0.0003041*T2 + 0.0000051*T3 );
259
260 XP.SinCos( SX, CX );
261 YP.SinCos( SY, CY );
262 ZP.SinCos( SZ, CZ );
263
264 //P1 is used to precess from any epoch to J2000
265 P1[0][0] = CX*CY*CZ - SX*SZ;
266 P1[1][0] = CX*CY*SZ + SX*CZ;
267 P1[2][0] = CX*SY;
268 P1[0][1] = -1.0*SX*CY*CZ - CX*SZ;
269 P1[1][1] = -1.0*SX*CY*SZ + CX*CZ;
270 P1[2][1] = -1.0*SX*SY;
271 P1[0][2] = -1.0*SY*CZ;
272 P1[1][2] = -1.0*SY*SZ;
273 P1[2][2] = CY;
274
275 //P2 is used to precess from J2000 to any other epoch (it is the transpose of P1)
276 P2[0][0] = CX*CY*CZ - SX*SZ;
277 P2[1][0] = -1.0*SX*CY*CZ - CX*SZ;
278 P2[2][0] = -1.0*SY*CZ;
279 P2[0][1] = CX*CY*SZ + SX*CZ;
280 P2[1][1] = -1.0*SX*CY*SZ + CX*CZ;
281 P2[2][1] = -1.0*SY*SZ;
282 P2[0][2] = CX*SY;
283 P2[1][2] = -1.0*SX*SY;
284 P2[2][2] = CY;
285
286 //Compute Precession Matrices from B1950 to 1984 using Newcomb formulae
287
288 XB.setD( 0.217697 );
289 YB.setD( 0.189274 );
290 ZB.setD( 0.217722 );
291
292 XB.SinCos( SXB, CXB );
293 YB.SinCos( SYB, CYB );
294 ZB.SinCos( SZB, CZB );
295
296 //P1B is used to precess from 1984 to B1950:
297
298 P1B[0][0] = CXB*CYB*CZB - SXB*SZB;
299 P1B[1][0] = CXB*CYB*SZB + SXB*CZB;
300 P1B[2][0] = CXB*SYB;
301 P1B[0][1] = -1.0*SXB*CYB*CZB - CXB*SZB;
302 P1B[1][1] = -1.0*SXB*CYB*SZB + CXB*CZB;
303 P1B[2][1] = -1.0*SXB*SYB;
304 P1B[0][2] = -1.0*SYB*CZB;
305 P1B[1][2] = -1.0*SYB*SZB;
306 P1B[2][2] = CYB;
307
308 //P2 is used to precess from B1950 to 1984 (it is the transpose of P1)
309 P2B[0][0] = CXB*CYB*CZB - SXB*SZB;
310 P2B[1][0] = -1.0*SXB*CYB*CZB - CXB*SZB;
311 P2B[2][0] = -1.0*SYB*CZB;
312 P2B[0][1] = CXB*CYB*SZB + SXB*CZB;
313 P2B[1][1] = -1.0*SXB*CYB*SZB + CXB*CZB;
314 P2B[2][1] = -1.0*SYB*SZB;
315 P2B[0][2] = CXB*SYB;
316 P2B[1][2] = -1.0*SXB*SYB;
317 P2B[2][2] = CYB;
318
319
320 // Mean longitudes for the planets. radians
321 //
322
323 // TODO Pasar a grados
324 double LVenus = 3.1761467+1021.3285546*T; // Venus
325 double LMars = 1.7534703+ 628.3075849*T; // Mars
326 double LEarth = 6.2034809+ 334.0612431*T; // Earth
327 double LJupiter = 0.5995465+ 52.9690965*T; // Jupiter
328 double LSaturn = 0.8740168+ 21.3299095*T; // Saturn
329 double LNeptune = 5.3118863+ 3.8133036*T; // Neptune
330 double LUranus = 5.4812939+ 7.4781599*T; // Uranus
331
332 double LMRad = 3.8103444+8399.6847337*T; // Moon
333 double DRad = 5.1984667+7771.3771486*T;
334 double MMRad = 2.3555559+8328.6914289*T; // Moon
335 double FRad = 1.6279052+8433.4661601*T;
336
337 /** Contibutions to the velocity of the Earth referred to the barycenter of the solar system
338 in the J2000 equatorial system
339 Velocities 10^{-8} AU/day
340 Ron & Vondrak method
341 **/
342
343 double vondrak[36][7] = {
344 {LMars, -1719914-2*T, -25, 25-13*T,1578089+156*T, 10+32*T,684185-358*T},
345 {2*LMars, 6434+141*T,28007-107*T,25697-95*T, -5904-130*T,11141-48*T, -2559-55*T},
346 {LJupiter, 715, 0, 6, -657, -15, -282},
347 {LMRad, 715, 0, 0, -656, 0, -285},
348 {3*LMars, 486-5*T, -236-4*T, -216-4*T, -446+5*T, -94, -193},
349 {LSaturn, 159, 0, 2, -147, -6, -61},
350 {FRad, 0, 0, 0, 26, 0, -59},
351 {LMRad+MMRad, 39, 0, 0, -36, 0, -16},
352 {2*LJupiter, 33, -10, -9, -30, -5, -13},
353 {2*LMars-LJupiter, 31, 1, 1, -28, 0, -12},
354 {3*LMars-8*LEarth+3*LJupiter, 8, -28, 25, 8, 11, 3},
355 {5*LMars-8*LEarth+3*LJupiter, 8, -28, -25, -8, -11, -3},
356 {2*LVenus-LMars, 21, 0, 0, -19, 0, -8},
357 {LVenus, -19, 0, 0, 17, 0, 8},
358 {LNeptune, 17, 0, 0, -16, 0, -7},
359 {LMars-2*LJupiter, 16, 0, 0, 15, 1, 7},
360 {LUranus, 16, 0, 1, -15, -3, -6},
361 {LMars+LJupiter, 11, -1, -1, -10, -1, -5},
362 {2*LVenus-2*LMars, 0, -11, -10, 0, -4, 0},
363 {LMars-LJupiter, -11, -2, -2, 9, -1, 4},
364 {4*LMars, -7, -8, -8, 6, -3, 3},
365 {3*LMars-2*LJupiter, -10, 0, 0, 9, 0, 4},
366 {LVenus-2*LMars, -9, 0, 0, -9, 0, -4},
367 {2*LVenus-3*LMars, -9, 0, 0, -8, 0, -4},
368 {2*LSaturn, 0, -9, -8, 0, -3, 0},
369 {2*LVenus-4*LMars, 0, -9, 8, 0, 3, 0},
370 {3*LMars-2*LEarth, 8, 0, 0, -8, 0, -3},
371 {LMRad+2*DRad-MMRad, 8, 0, 0, -7, 0, -3},
372 {8*LVenus-12*LMars, -4, -7, -6, 4, -3, 2},
373 {8*LVenus-14*LMars, -4, -7, 6, -4, 3, -2},
374 {2*LEarth, -6, -5, -4, 5, -2, 2},
375 {3*LVenus-4*LMars, -1, -1, -2, -7, 1, -4},
376 {2*LMars-2*LJupiter, 4, -6, -5, -4, -2, -2},
377 {3*LVenus-3*LMars, 0, -7, -6, 0, -3, 0},
378 {2*LMars-2*LEarth, 5, -5, -4, -5, -2, -2},
379 {LMRad-2*DRad, 5, 0, 0, -5, 0, -2}
380 };
381
382 dms anglev;
383 double sa, ca;
384 // Vearth X component
385 vearth[0] = 0.;
386 // Vearth Y component
387 vearth[1] = 0.;
388 // Vearth Z component
389 vearth[2] = 0.;
390
391 for (unsigned int i=0; i<36; i++) {
392 anglev.setRadians(vondrak[i][0]);
393 anglev.SinCos(sa,ca);
394 for (unsigned int j=0; j<3; j++) {
395 vearth[j] += vondrak[i][2*j+1]*sa +vondrak[i][2*j+2]*ca;
396 }
397 }
398
399 const double UA2km = 1.49597870/86400.; // 10^{-8}*UA/dia -> km/s
400
401 for (unsigned int j=0; j<3; j++) {
402 vearth[j] = vearth[j] * UA2km;
403 }
404 }
Experimental kdeedu stuff