1 /*
2 OpenUniverse 1.0
3 Copyright (C) 2000 Raul Alonso <amil@las.es>
4
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation; either version 2 of the License, or
8 (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software
17 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
18 */
19
20
21 /*
22 * Astronomical library
23 * Written February 1999 by Larry Manley
24 *
25 * Used in openuniverse with permission, thanks Larry :-)
26 *
27 * Based on information in "Astronomical Algorithms", by Jean Meeus, 1991
28 *
29 *
30 *
31 * Here you'll find the custom functions for the nine planets and the Moon.
32 * They return longitude, latitude and distance from the Sun in AU (the
33 * MoonPos and every custom satellite funcition must return the distance
34 * from its host planet in planet radii)
35 */
36
37 #include <stdio.h>
38 #include <math.h>
39
40 #ifndef M_PI
41 #define M_PI 3.14159265358979323846
42 #endif
43 #define DEG2RAD(x) ((x)*M_PI/180.0)
44 #define RAD2DEG(x) ((x)*180.0/M_PI)
45
46 #include "vsop87.dat" /* Hardcoded VSOP87 data */
47
48 /* Not used elsewhere yet ... but it will be */
JulianDay(int Y,int M,int D)49 double JulianDay(int Y, int M, int D)
50 {
51 double JD;
52 int A, B;
53
54 if (M < 3) {
55 Y -= 1;
56 M += 12;
57 }
58 A = (int) floor((double) Y / 100.0);
59 B = 2 - A + (int) floor((double) A / 4.0);
60 JD =
61 floor(365.25 * ((double) Y + 4716.0)) +
62 floor(30.6001 * ((double) M + 1)) + D + B - 1524.5;
63
64 return JD;
65 }
66
67 /* Planetary positioning algorithms
68 * Abridged VSOP87 theory
69 * Longitude and latitude calculated in radians
70 * Radius calculated in AU
71 */
72
ReadTerms(double T,double * data)73 double ReadTerms(double T, double *data)
74 {
75 double sum = 0.0;
76 int i, j;
77
78 i = (int) data[0];
79 for (j = 0; j < i; j++)
80 sum +=
81 (data[j * 3 + 1] *
82 cos(data[j * 3 + 2] + (data[j * 3 + 3] * T)));
83 return (sum / 1.0E8);
84 }
85
86
MercuryPos(double JD,double * Long,double * Lat,double * Rad)87 void MercuryPos(double JD, double *Long, double *Lat, double *Rad)
88 {
89 double t, L0, L1, L2, L3, L4, L5;
90 double B0, B1, B2, B3, B4;
91 double R0, R1, R2, R3;
92
93 t = (JD - 2451545.0) / 365250.0; /* Calculate tau */
94
95 L0 = ReadTerms(t, mercuryL0);
96 L1 = ReadTerms(t, mercuryL1);
97 L2 = ReadTerms(t, mercuryL2);
98 L3 = ReadTerms(t, mercuryL3);
99 L4 = ReadTerms(t, mercuryL4);
100 L5 = ReadTerms(t, mercuryL5);
101 *Long = L0 + t * (L1 + t * (L2 + t * (L3 + t * (L4 + t * L5))));
102
103 B0 = ReadTerms(t, mercuryB0);
104 B1 = ReadTerms(t, mercuryB1);
105 B2 = ReadTerms(t, mercuryB2);
106 B3 = ReadTerms(t, mercuryB3);
107 B4 = ReadTerms(t, mercuryB4);
108 *Lat = B0 + t * (B1 + t * (B2 + t * (B3 + t * B4)));
109
110 R0 = ReadTerms(t, mercuryR0);
111 R1 = ReadTerms(t, mercuryR1);
112 R2 = ReadTerms(t, mercuryR2);
113 R3 = ReadTerms(t, mercuryR3);
114 *Rad = R0 + t * (R1 + t * (R2 + t * R3));
115 }
116
VenusPos(double JD,double * Long,double * Lat,double * Rad)117 void VenusPos(double JD, double *Long, double *Lat, double *Rad)
118 {
119 double t, L0, L1, L2, L3, L4, L5;
120 double B0, B1, B2, B3, B4;
121 double R0, R1, R2, R3, R4;
122
123 t = (JD - 2451545.0) / 365250.0; /* Calculate tau */
124
125 L0 = ReadTerms(t, venusL0);
126 L1 = ReadTerms(t, venusL1);
127 L2 = ReadTerms(t, venusL2);
128 L3 = ReadTerms(t, venusL3);
129 L4 = ReadTerms(t, venusL4);
130 L5 = ReadTerms(t, venusL5);
131 *Long = L0 + t * (L1 + t * (L2 + t * (L3 + t * (L4 + t * L5))));
132
133 B0 = ReadTerms(t, venusB0);
134 B1 = ReadTerms(t, venusB1);
135 B2 = ReadTerms(t, venusB2);
136 B3 = ReadTerms(t, venusB3);
137 B4 = ReadTerms(t, venusB4);
138 *Lat = B0 + t * (B1 + t * (B2 + t * (B3 + t * B4)));
139
140 R0 = ReadTerms(t, venusR0);
141 R1 = ReadTerms(t, venusR1);
142 R2 = ReadTerms(t, venusR2);
143 R3 = ReadTerms(t, venusR3);
144 R4 = ReadTerms(t, venusR4);
145 *Rad = R0 + t * (R1 + t * (R2 + t * (R3 + t * R4)));
146 }
147
EarthPos(double JD,double * Long,double * Lat,double * Rad)148 void EarthPos(double JD, double *Long, double *Lat, double *Rad)
149 {
150 double t, L0, L1, L2, L3, L4, L5;
151 double B0, B1;
152 double R0, R1, R2, R3, R4;
153
154 t = (JD - 2451545.0) / 365250.0; /* Calculate tau */
155
156 L0 = ReadTerms(t, earthL0);
157 L1 = ReadTerms(t, earthL1);
158 L2 = ReadTerms(t, earthL2);
159 L3 = ReadTerms(t, earthL3);
160 L4 = ReadTerms(t, earthL4);
161 L5 = ReadTerms(t, earthL5);
162 *Long = L0 + t * (L1 + t * (L2 + t * (L3 + t * (L4 + t * L5))));
163
164 B0 = ReadTerms(t, earthB0);
165 B1 = ReadTerms(t, earthB1);
166 *Lat = B0 + t * B1;
167
168 R0 = ReadTerms(t, earthR0);
169 R1 = ReadTerms(t, earthR1);
170 R2 = ReadTerms(t, earthR2);
171 R3 = ReadTerms(t, earthR3);
172 R4 = ReadTerms(t, earthR4);
173 *Rad = R0 + t * (R1 + t * (R2 + t * (R3 + t * R4)));
174 }
175
MarsPos(double JD,double * Long,double * Lat,double * Rad)176 void MarsPos(double JD, double *Long, double *Lat, double *Rad)
177 {
178 double t, L0, L1, L2, L3, L4, L5;
179 double B0, B1, B2, B3, B4;
180 double R0, R1, R2, R3, R4;
181
182 t = (JD - 2451545.0) / 365250.0; /* Calculate tau */
183
184 L0 = ReadTerms(t, marsL0);
185 L1 = ReadTerms(t, marsL1);
186 L2 = ReadTerms(t, marsL2);
187 L3 = ReadTerms(t, marsL3);
188 L4 = ReadTerms(t, marsL4);
189 L5 = ReadTerms(t, marsL5);
190 *Long = L0 + t * (L1 + t * (L2 + t * (L3 + t * (L4 + t * L5))));
191
192 B0 = ReadTerms(t, marsB0);
193 B1 = ReadTerms(t, marsB1);
194 B2 = ReadTerms(t, marsB2);
195 B3 = ReadTerms(t, marsB3);
196 B4 = ReadTerms(t, marsB4);
197 *Lat = B0 + t * (B1 + t * (B2 + t * (B3 + t * B4)));
198
199 R0 = ReadTerms(t, marsR0);
200 R1 = ReadTerms(t, marsR1);
201 R2 = ReadTerms(t, marsR2);
202 R3 = ReadTerms(t, marsR3);
203 R4 = ReadTerms(t, marsR4);
204 *Rad = R0 + t * (R1 + t * (R2 + t * (R3 + t * R4)));
205 }
206
JupiterPos(double JD,double * Long,double * Lat,double * Rad)207 void JupiterPos(double JD, double *Long, double *Lat, double *Rad)
208 {
209 double t, L0, L1, L2, L3, L4, L5;
210 double B0, B1, B2, B3, B4, B5;
211 double R0, R1, R2, R3, R4, R5;
212
213 t = (JD - 2451545.0) / 365250.0; /* Calculate tau */
214
215 L0 = ReadTerms(t, jupiterL0);
216 L1 = ReadTerms(t, jupiterL1);
217 L2 = ReadTerms(t, jupiterL2);
218 L3 = ReadTerms(t, jupiterL3);
219 L4 = ReadTerms(t, jupiterL4);
220 L5 = ReadTerms(t, jupiterL5);
221 *Long = L0 + t * (L1 + t * (L2 + t * (L3 + t * (L4 + t * L5))));
222
223 B0 = ReadTerms(t, jupiterB0);
224 B1 = ReadTerms(t, jupiterB1);
225 B2 = ReadTerms(t, jupiterB2);
226 B3 = ReadTerms(t, jupiterB3);
227 B4 = ReadTerms(t, jupiterB4);
228 B5 = ReadTerms(t, jupiterB5);
229 *Lat = B0 + t * (B1 + t * (B2 + t * (B3 + t * (B4 + t * B5))));
230
231 R0 = ReadTerms(t, jupiterR0);
232 R1 = ReadTerms(t, jupiterR1);
233 R2 = ReadTerms(t, jupiterR2);
234 R3 = ReadTerms(t, jupiterR3);
235 R4 = ReadTerms(t, jupiterR4);
236 R5 = ReadTerms(t, jupiterR5);
237 *Rad = R0 + t * (R1 + t * (R2 + t * (R3 + t * (R4 + t * R5))));
238 }
239
SaturnPos(double JD,double * Long,double * Lat,double * Rad)240 void SaturnPos(double JD, double *Long, double *Lat, double *Rad)
241 {
242 double t, L0, L1, L2, L3, L4, L5;
243 double B0, B1, B2, B3, B4, B5;
244 double R0, R1, R2, R3, R4, R5;
245
246 t = (JD - 2451545.0) / 365250.0; /* Calculate tau */
247
248 L0 = ReadTerms(t, saturnL0);
249 L1 = ReadTerms(t, saturnL1);
250 L2 = ReadTerms(t, saturnL2);
251 L3 = ReadTerms(t, saturnL3);
252 L4 = ReadTerms(t, saturnL4);
253 L5 = ReadTerms(t, saturnL5);
254 *Long = L0 + t * (L1 + t * (L2 + t * (L3 + t * (L4 + t * L5))));
255
256 B0 = ReadTerms(t, saturnB0);
257 B1 = ReadTerms(t, saturnB1);
258 B2 = ReadTerms(t, saturnB2);
259 B3 = ReadTerms(t, saturnB3);
260 B4 = ReadTerms(t, saturnB4);
261 B5 = ReadTerms(t, saturnB5);
262 *Lat = B0 + t * (B1 + t * (B2 + t * (B3 + t * (B4 + t * B5))));
263
264 R0 = ReadTerms(t, saturnR0);
265 R1 = ReadTerms(t, saturnR1);
266 R2 = ReadTerms(t, saturnR2);
267 R3 = ReadTerms(t, saturnR3);
268 R4 = ReadTerms(t, saturnR4);
269 R5 = ReadTerms(t, saturnR5);
270 *Rad = R0 + t * (R1 + t * (R2 + t * (R3 + t * (R4 + t * R5))));
271 }
272
UranusPos(double JD,double * Long,double * Lat,double * Rad)273 void UranusPos(double JD, double *Long, double *Lat, double *Rad)
274 {
275 double t, L0, L1, L2, L3, L4;
276 double B0, B1, B2, B3, B4;
277 double R0, R1, R2, R3, R4;
278
279 t = (JD - 2451545.0) / 365250.0; /* Calculate tau */
280
281 L0 = ReadTerms(t, uranusL0);
282 L1 = ReadTerms(t, uranusL1);
283 L2 = ReadTerms(t, uranusL2);
284 L3 = ReadTerms(t, uranusL3);
285 L4 = ReadTerms(t, uranusL4);
286 *Long = L0 + t * (L1 + t * (L2 + t * (L3 + t * L4)));
287
288 B0 = ReadTerms(t, uranusB0);
289 B1 = ReadTerms(t, uranusB1);
290 B2 = ReadTerms(t, uranusB2);
291 B3 = ReadTerms(t, uranusB3);
292 B4 = ReadTerms(t, uranusB4);
293 *Lat = B0 + t * (B1 + t * (B2 + t * (B3 + t * B4)));
294
295 R0 = ReadTerms(t, uranusR0);
296 R1 = ReadTerms(t, uranusR1);
297 R2 = ReadTerms(t, uranusR2);
298 R3 = ReadTerms(t, uranusR3);
299 R4 = ReadTerms(t, uranusR4);
300 *Rad = R0 + t * (R1 + t * (R2 + t * (R3 + t * R4)));
301 }
302
NeptunePos(double JD,double * Long,double * Lat,double * Rad)303 void NeptunePos(double JD, double *Long, double *Lat, double *Rad)
304 {
305 double t, L0, L1, L2, L3, L4;
306 double B0, B1, B2, B3, B4;
307 double R0, R1, R2, R3;
308
309 t = (JD - 2451545.0) / 365250.0; /* Calculate tau */
310
311 L0 = ReadTerms(t, neptuneL0);
312 L1 = ReadTerms(t, neptuneL1);
313 L2 = ReadTerms(t, neptuneL2);
314 L3 = ReadTerms(t, neptuneL3);
315 L4 = ReadTerms(t, neptuneL4);
316 *Long = L0 + t * (L1 + t * (L2 + t * (L3 + t * L4)));
317
318 B0 = ReadTerms(t, neptuneB0);
319 B1 = ReadTerms(t, neptuneB1);
320 B2 = ReadTerms(t, neptuneB2);
321 B3 = ReadTerms(t, neptuneB3);
322 B4 = ReadTerms(t, neptuneB4);
323 *Lat = B0 + t * (B1 + t * (B2 + t * (B3 + t * B4)));
324
325 R0 = ReadTerms(t, neptuneR0);
326 R1 = ReadTerms(t, neptuneR1);
327 R2 = ReadTerms(t, neptuneR2);
328 R3 = ReadTerms(t, neptuneR3);
329 *Rad = R0 + t * (R1 + t * (R2 + t * R3));
330 }
331
PlutoPos(double JD,double * Long,double * Lat,double * Rad)332 void PlutoPos(double JD, double *Long, double *Lat, double *Rad)
333 {
334 double j, s, p, a, T, sumL, sumB, sumR, sa, ca;
335 double j2, j3, s2, p2, p3, p4, p5, p6;
336
337 sumL = sumB = sumR = 0.0;
338
339 T = (JD - 2451545.0) / 36525.0;
340
341 j = 34.35 + (3034.9057 * T);
342 j = DEG2RAD(j);
343 j2 = 2.0 * j;
344 j3 = 3.0 * j;
345
346 s = 50.08 + (1222.1138 * T);
347 s = DEG2RAD(s);
348 s2 = 2.0 * s;
349
350 p = 238.96 + (144.9600 * T);
351 p = DEG2RAD(p);
352 p2 = 2.0 * p;
353 p3 = 3.0 * p;
354 p4 = 4.0 * p;
355 p5 = 5.0 * p;
356 p6 = 6.0 * p;
357
358 a = p;
359 sa = sin(a);
360 ca = cos(a);
361 sumL += (-19798886 * sa) + (19848454 * ca);
362 sumB += (-5453098 * sa) + (-14974876 * ca);
363 sumR += (66867334 * sa) + (68955876 * ca);
364
365 a = p2;
366 sa = sin(a);
367 ca = cos(a);
368 sumL += (897499 * sa) + (-4955707 * ca);
369 sumB += (3527363 * sa) + (1672673 * ca);
370 sumR += (-11826086 * sa) + (-333765 * ca);
371
372 a = p3;
373 sa = sin(a);
374 ca = cos(a);
375 sumL += (610820 * sa) + (1210521 * ca);
376 sumB += (-1050939 * sa) + (327763 * ca);
377 sumR += (1593657 * sa) + (-1439953 * ca);
378
379 a = p4;
380 sa = sin(a);
381 ca = cos(a);
382 sumL += (-341639 * sa) + (-189719 * ca);
383 sumB += (178691 * sa) + (-291925 * ca);
384 sumR += (-18948 * sa) + (482443 * ca);
385
386 a = p5;
387 sa = sin(a);
388 ca = cos(a);
389 sumL += (129027 * sa) + (-34863 * ca);
390 sumB += (18763 * sa) + (100448 * ca);
391 sumR += (-66634 * sa) + (-85576 * ca);
392
393 a = p6;
394 sa = sin(a);
395 ca = cos(a);
396 sumL += (-38215 * sa) + (31061 * ca);
397 sumB += (-30594 * sa) + (-25838 * ca);
398 sumR += (30841 * sa) + (-5765 * ca);
399
400 a = s - p;
401 sa = sin(a);
402 ca = cos(a);
403 sumL += (20349 * sa) + (-9886 * ca);
404 sumB += (4965 * sa) + (11263 * ca);
405 sumR += (-6140 * sa) + (22254 * ca);
406
407 a = s;
408 sa = sin(a);
409 ca = cos(a);
410 sumL += (-4045 * sa) + (-4904 * ca);
411 sumB += (310 * sa) + (-132 * ca);
412 sumR += (4434 * sa) + (4443 * ca);
413
414 a = s + p;
415 sa = sin(a);
416 ca = cos(a);
417 sumL += (-5885 * sa) + (-3238 * ca);
418 sumB += (2036 * sa) + (-947 * ca);
419 sumR += (-1518 * sa) + (641 * ca);
420
421 a = s + p2;
422 sa = sin(a);
423 ca = cos(a);
424 sumL += (-3812 * sa) + (3011 * ca);
425 sumB += (-2 * sa) + (-674 * ca);
426 sumR += (-5 * sa) + (792 * ca);
427
428 a = s + p3;
429 sa = sin(a);
430 ca = cos(a);
431 sumL += (-601 * sa) + (3468 * ca);
432 sumB += (-329 * sa) + (-563 * ca);
433 sumR += (518 * sa) + (518 * ca);
434
435 a = s2 - p2;
436 sa = sin(a);
437 ca = cos(a);
438 sumL += (1237 * sa) + (463 * ca);
439 sumB += (-64 * sa) + (39 * ca);
440 sumR += (-13 * sa) + (-221 * ca);
441
442 a = s2 - p;
443 sa = sin(a);
444 ca = cos(a);
445 sumL += (1086 * sa) + (-911 * ca);
446 sumB += (-94 * sa) + (210 * ca);
447 sumR += (837 * sa) + (-494 * ca);
448
449 a = s2;
450 sa = sin(a);
451 ca = cos(a);
452 sumL += (595 * sa) + (-1229 * ca);
453 sumB += (-8 * sa) + (-160 * ca);
454 sumR += (-281 * sa) + (616 * ca);
455
456 a = j - s;
457 sa = sin(a);
458 ca = cos(a);
459 sumL += (2484 * sa) + (-485 * ca);
460 sumB += (-177 * sa) + (259 * ca);
461 sumR += (260 * sa) + (-395 * ca);
462
463 a = j - s + p;
464 sa = sin(a);
465 ca = cos(a);
466 sumL += (839 * sa) + (-1414 * ca);
467 sumB += (17 * sa) + (234 * ca);
468 sumR += (-191 * sa) + (-396 * ca);
469
470 a = j - p3;
471 sa = sin(a);
472 ca = cos(a);
473 sumL += (-964 * sa) + (1059 * ca);
474 sumB += (582 * sa) + (-285 * ca);
475 sumR += (-3218 * sa) + (370 * ca);
476
477 a = j - p2;
478 sa = sin(a);
479 ca = cos(a);
480 sumL += (-2303 * sa) + (-1038 * ca);
481 sumB += (-298 * sa) + (692 * ca);
482 sumR += (8019 * sa) + (-7869 * ca);
483
484 a = j - p;
485 sa = sin(a);
486 ca = cos(a);
487 sumL += (7049 * sa) + (747 * ca);
488 sumB += (157 * sa) + (201 * ca);
489 sumR += (105 * sa) + (45637 * ca);
490
491 a = j;
492 sa = sin(a);
493 ca = cos(a);
494 sumL += (1179 * sa) + (-358 * ca);
495 sumB += (304 * sa) + (825 * ca);
496 sumR += (8623 * sa) + (8444 * ca);
497
498 a = j + p;
499 sa = sin(a);
500 ca = cos(a);
501 sumL += (393 * sa) + (-63 * ca);
502 sumB += (-124 * sa) + (-29 * ca);
503 sumR += (-896 * sa) + (-801 * ca);
504
505 a = j + p2;
506 sa = sin(a);
507 ca = cos(a);
508 sumL += (111 * sa) + (-268 * ca);
509 sumB += (15 * sa) + (8 * ca);
510 sumR += (208 * sa) + (-122 * ca);
511
512 a = j + p3;
513 sa = sin(a);
514 ca = cos(a);
515 sumL += (-52 * sa) + (-154 * ca);
516 sumB += (7 * sa) + (15 * ca);
517 sumR += (-133 * sa) + (65 * ca);
518
519 a = j + p4;
520 sa = sin(a);
521 ca = cos(a);
522 sumL += (-78 * sa) + (-30 * ca);
523 sumB += (2 * sa) + (2 * ca);
524 sumR += (-16 * sa) + ca;
525
526 a = j + s - p3;
527 sa = sin(a);
528 ca = cos(a);
529 sumL += (-34 * sa) + (-26 * ca);
530 sumB += (4 * sa) + (2 * ca);
531 sumR += (-22 * sa) + (7 * ca);
532
533 a = j + s - p2;
534 sa = sin(a);
535 ca = cos(a);
536 sumL += (-43 * sa) + ca;
537 sumB += 3 * sa;
538 sumR += (-8 * sa) + (16 * ca);
539
540 a = j + s - p;
541 sa = sin(a);
542 ca = cos(a);
543 sumL += (-15 * sa) + (21 * ca);
544 sumB += sa - ca;
545 sumR += (2 * sa) + (9 * ca);
546
547 a = j + s;
548 sa = sin(a);
549 ca = cos(a);
550 sumL += (-1 * sa) + (15 * ca);
551 sumB += -2 * ca;
552 sumR += (12 * sa) + (5 * ca);
553
554 a = j + s + p;
555 sa = sin(a);
556 ca = cos(a);
557 sumL += (4 * sa) + (7 * ca);
558 sumB += sa;
559 sumR += sa - 3.0 * ca;
560
561 a = j + s + p3;
562 sa = sin(a);
563 ca = cos(a);
564 sumL += sa + (5.0 * ca);
565 sumB += sa - ca;
566 sumR += sa;
567
568 a = j2 - p6;
569 sa = sin(a);
570 ca = cos(a);
571 sumL += (8 * sa) + (3 * ca);
572 sumB += (-2 * sa) + (-3 * ca);
573 sumR += (9 * sa) + (5 * ca);
574
575 a = j2 - p5;
576 sa = sin(a);
577 ca = cos(a);
578 sumL += (-3.0 * sa) + (6 * ca);
579 sumB += sa + (2 * ca);
580 sumR += (2 * sa) - ca;
581
582 a = j2 - p4;
583 sa = sin(a);
584 ca = cos(a);
585 sumL += (6 * sa) + (-13 * ca);
586 sumB += (-8 * sa) + (2 * ca);
587 sumR += (14 * sa) + (10 * ca);
588
589 a = j2 - p3;
590 sa = sin(a);
591 ca = cos(a);
592 sumL += (10 * sa) + (22 * ca);
593 sumB += (10 * sa) + (-7 * ca);
594 sumR += (-65 * sa) + (12 * ca);
595
596 a = j2 - p2;
597 sa = sin(a);
598 ca = cos(a);
599 sumL += (-57 * sa) + (-32 * ca);
600 sumB += 21 * ca;
601 sumR += (126 * sa) + (-233 * ca);
602
603 a = j2 - p;
604 sa = sin(a);
605 ca = cos(a);
606 sumL += (157 * sa) + (-46 * ca);
607 sumB += (8 * sa) + (5 * ca);
608 sumR += (270 * sa) + (1068 * ca);
609
610 a = j2;
611 sa = sin(a);
612 ca = cos(a);
613 sumL += (12 * sa) + (-18 * ca);
614 sumB += (13 * sa) + (16 * ca);
615 sumR += (254 * sa) + (155 * ca);
616
617 a = j2 + p;
618 sa = sin(a);
619 ca = cos(a);
620 sumL += (-4 * sa) + (8 * ca);
621 sumB += (-2 * sa) + (-3 * ca);
622 sumR += (-26 * sa) + (-2 * ca);
623
624 a = j2 + p2;
625 sa = sin(a);
626 ca = cos(a);
627 sumL += -5 * sa;
628 sumR += 7 * sa;
629
630 a = j2 + p3;
631 sa = sin(a);
632 ca = cos(a);
633 sumL += (3 * sa) + (4 * ca);
634 sumB += ca;
635 sumR += (-11 * sa) + (4 * ca);
636
637 a = j3 - p2;
638 sa = sin(a);
639 ca = cos(a);
640 sumL += (-1 * sa) + (-1 * ca);
641 sumB += ca;
642 sumR += (4 * sa) + (-14 * ca);
643
644 a = j3 - p;
645 sa = sin(a);
646 ca = cos(a);
647 sumL += (6 * sa) + (-3 * ca);
648 sumR += (18 * sa) + (35 * ca);
649
650 a = j3;
651 sa = sin(a);
652 ca = cos(a);
653 sumL += (-1 * sa) + (-2 * ca);
654 sumB += ca;
655 sumR += (13 * sa) + (3 * ca);
656
657 *Long = 238.956785 + (144.96 * T) + (sumL / 1E6);
658 while (*Long > 360.0)
659 *Long -= 360.0;
660 while (*Long < 0.0)
661 *Long += 360.0;
662
663 *Lat = -3.98202 + (sumB / 1E6);
664 while (*Lat > 360.0)
665 *Lat -= 360.0;
666 while (*Lat < 0.0)
667 *Lat += 360.0;
668
669 *Rad = 40.7247248 + (sumR / 1E7);
670 *Long = DEG2RAD(*Long);
671 *Lat = DEG2RAD(*Lat);
672 }
673
674
MoonPos(double JD,double * Long,double * Lat,double * Rad)675 void MoonPos(double JD, double *Long, double *Lat, double *Rad)
676 {
677 double Lp, Lpr, D, Dr, M, Mr, Mp, Mpr, F, Fr, A1, A2, A3, E;
678 double T, T2, T3, T4, sumL, sumB, sumR;
679
680 sumL = sumB = sumR = 0.0;
681 T = (JD - 2451545.0) / 36525.0;
682 T2 = T * T;
683 T3 = T * T2;
684 T4 = T3 * T;
685
686 /* Compute moon's mean longitude <Lp>, in degrees */
687 Lp = 218.3164591 + (481267.88134236 * T) - (0.0013268 * T2) +
688 (T3 / 538841.0) - (T4 / 65194000);
689 while (Lp > 360.0)
690 Lp -= 360.0;
691 while (Lp < 0.0)
692 Lp += 360.0;
693 Lpr = DEG2RAD(Lp); /* Create radians version for trig functions below */
694
695 /* Mean elongation of the moon <D>, in degrees */
696 D = 297.8502042 + (445267.1115168 * T) - (0.0016300 * T2) +
697 (T3 / 545868.0) - (T4 / 113065000.0);
698 while (D > 360.0)
699 D -= 360.0;
700 while (D < 0.0)
701 D += 360.0;
702 Dr = DEG2RAD(D);
703
704 /* Sun's mean anomaly <M>, in degrees */
705 M = 357.5291092 + (35999.0502909 * T) - (0.0001536 * T2) +
706 (T3 / 24490000.0);
707 while (M > 360.0)
708 M -= 360.0;
709 while (M < 0.0)
710 M += 360.0;
711 Mr = DEG2RAD(M);
712
713 /* Moon's mean anomaly <Mp>, in degrees */
714 Mp = 134.9634114 + (477198.8676313 * T) + (0.0089970 * T2) +
715 (T3 / 69699.0) - (T4 / 14712000.0);
716 while (Mp > 360.0)
717 Mp -= 360.0;
718 while (Mp < 0.0)
719 Mp += 360.0;
720 Mpr = DEG2RAD(Mp);
721
722 /* Moon's argument of latitude <F> */
723 F = 93.2720993 + (483202.0175273 * T) - (0.0034029 * T2) -
724 (T3 / 3526000.0) + (T4 / 863310000.0);
725 while (F > 360.0)
726 F -= 360.0;
727 while (F < 0.0)
728 F += 360.0;
729 Fr = DEG2RAD(F);
730
731 A1 = 119.75 + (131.849 * T);
732 while (A1 > 360.0)
733 A1 -= 360.0;
734 while (A1 < 0.0)
735 A1 += 360.0;
736
737 A2 = 53.09 + (479264.290 * T);
738 while (A2 > 360.0)
739 A2 -= 360.0;
740 while (A2 < 0.0)
741 A2 += 360.0;
742
743 A3 = 313.45 + (481266.484 * T);
744 while (A3 > 360.0)
745 A3 -= 360.0;
746 while (A3 < 0.0)
747 A3 += 360.0;
748
749 E = 1.0 - (0.002516 * T) - (0.0000074 * T2);
750
751 A1 = DEG2RAD(A1);
752 A2 = DEG2RAD(A2);
753 A3 = DEG2RAD(A3);
754
755 sumL += 6288774 * sin(Mpr);
756 sumL += 1274027 * sin(2.0 * Dr - Mpr);
757 sumL += 658314 * sin(2.0 * Dr);
758 sumL += 213618 * sin(2.0 * Mpr);
759 sumL += (-185116 * sin(Mr)) * E;
760 sumL += -114332 * sin(2.0 * Fr);
761 sumL += 58793 * sin(2.0 * Dr - 2.0 * Mpr);
762 sumL += (57066 * sin(2.0 * Dr - Mr - Mpr)) * E;
763 sumL += 53322 * sin(2.0 * Dr + Mpr);
764 sumL += (45758 * sin(2.0 * Dr - Mr)) * E;
765 sumL += (-40923 * sin(Mr - Mpr)) * E;
766 sumL += -34720 * sin(Dr);
767 sumL += (-30383 * sin(Mr + Mpr)) * E;
768 sumL += 15327 * sin(2.0 * Dr - 2.0 * Fr);
769 sumL += -12528 * sin(Mpr + 2.0 * Fr);
770 sumL += 10980 * sin(Mpr - 2.0 * Fr);
771 sumL += 10675 * sin(4.0 * Dr - Mpr);
772 sumL += 10034 * sin(3.0 * Mpr);
773 sumL += 8548 * sin(4.0 * Dr - 2.0 * Mpr);
774 sumL += (-7888 * sin(2.0 * Dr + Mr - Mpr)) * E;
775 sumL += (-6766 * sin(2.0 * Dr + Mr)) * E;
776 sumL += -5163 * sin(Dr - Mpr);
777 sumL += (4987 * sin(Dr + Mr)) * E;
778 sumL += (-4036 * sin(2.0 * Dr - Mr + Mpr)) * E;
779 sumL += 3994 * sin(2.0 * Dr + 2.0 * Mpr);
780 sumL += 3861 * sin(4.0 * Dr);
781 sumL += 3665 * sin(2.0 * Dr - 3.0 * Mpr);
782 sumL += (-2689 * sin(Mr - 2.0 * Mpr)) * E;
783 sumL += -2602 * sin(2.0 * Dr - Mpr + 2.0 * Fr);
784 sumL += (2390 * sin(2.0 * Dr - Mr - 2.0 * Mpr)) * E;
785 sumL += -2348 * sin(Dr + Mpr);
786 sumL += (2236 * sin(2.0 * Dr - 2.0 * Mr)) * E * E;
787 sumL += (-2120 * sin(Mr + 2.0 * Mpr)) * E;
788 sumL += (-2069 * sin(2.0 * Mr)) * E * E;
789 sumL += (2048 * sin(2.0 * Dr - 2.0 * Mr - Mpr)) * E * E;
790 sumL += -1773 * sin(2.0 * Dr + Mpr - 2.0 * Fr);
791 sumL += -1595 * sin(2.0 * Dr + 2.0 * Fr);
792 sumL += (1215 * sin(4.0 * Dr - Mr - Mpr)) * E;
793 sumL += -1110 * sin(2.0 * Mpr + 2.0 * Fr);
794 sumL += -892 * sin(3.0 * Dr - Mpr);
795 sumL += (-810 * sin(2.0 * Dr + Mr + Mpr)) * E;
796 sumL += (759 * sin(4.0 * Dr - Mr - 2.0 * Mpr)) * E;
797 sumL += (-713 * sin(2.0 * Mr - Mpr)) * E * E;
798 sumL += (-700 * sin(2.0 * Dr + 2.0 * Mr - Mpr)) * E * E;
799 sumL += (691 * sin(2.0 * Dr + Mr - 2.0 * Mpr)) * E;
800 sumL += (596 * sin(2.0 * Dr - Mr - 2.0 * F)) * E;
801 sumL += 549 * sin(4.0 * Dr + Mpr);
802 sumL += 537 * sin(4.0 * Mpr);
803 sumL += (520 * sin(4.0 * Dr - Mr)) * E;
804 sumL += -487 * sin(Dr - 2.0 * Mpr);
805 sumL += (-399 * sin(2.0 * Dr + Mr - 2.0 * Fr)) * E;
806 sumL += -381 * sin(2.0 * Mpr - 2.0 * Fr);
807 sumL += (351 * sin(Dr + Mr + Mpr)) * E;
808 sumL += -340 * sin(3.0 * Dr - 2.0 * Mpr);
809 sumL += 330 * sin(4.0 * Dr - 3.0 * Mpr);
810 sumL += (327 * sin(2.0 * Dr - Mr + 2.0 * Mpr)) * E;
811 sumL += (-323 * sin(2.0 * Mr + Mpr)) * E * E;
812 sumL += (299 * sin(Dr + Mr - Mpr)) * E;
813 sumL += 294 * sin(2.0 * Dr + 3.0 * Mpr);
814 sumL += 3958 * sin(A1);
815 sumL += 1962 * sin(Lpr - Fr);
816 sumL += 318 * sin(A2);
817 *Long = Lp + (sumL / 1E6);
818 while (*Long > 360.0)
819 *Long -= 360.0;
820 while (*Long < 0.0)
821 *Long += 360.0;
822
823 sumB += 5128122 * sin(Fr);
824 sumB += 280602 * sin(Mpr + Fr);
825 sumB += 277693 * sin(Mpr - Fr);
826 sumB += 173237 * sin(2.0 * Dr - Fr);
827 sumB += 55413 * sin(2.0 * Dr - Mpr + Fr);
828 sumB += 46271 * sin(2.0 * Dr - Mpr - Fr);
829 sumB += 32573 * sin(2.0 * Dr + Fr);
830 sumB += 17198 * sin(2.0 * Mpr + Fr);
831 sumB += 9266 * sin(2.0 * Dr + Mpr - Fr);
832 sumB += 8822 * sin(2.0 * Mpr - Fr);
833 sumB += (8216 * sin(2.0 * Dr - Mr - Fr)) * E;
834 sumB += 4324 * sin(2.0 * Dr - 2.0 * Mpr - Fr);
835 sumB += 4200 * sin(2.0 * Dr + Mpr + Fr);
836 sumB += (-3359 * sin(2.0 * Dr + Mr - Fr)) * E;
837 sumB += (2463 * sin(2.0 * Dr - Mr - Mpr + Fr)) * E;
838 sumB += (2211 * sin(2.0 * Dr - Mr + Fr)) * E;
839 sumB += (2065 * sin(2.0 * Dr - Mr - Mpr - Fr)) * E;
840 sumB += (-1870 * sin(Mr - Mpr - Fr)) * E;
841 sumB += 1828 * sin(4.0 * Dr - Mpr - Fr);
842 sumB += (-1794 * sin(Mr + Fr)) * E;
843 sumB += -1749 * sin(3.0 * Fr);
844 sumB += (-1565 * sin(Mr - Mpr - Fr)) * E;
845 sumB += -1491 * sin(Dr - Fr);
846 sumB += (-1475 * sin(Mr + Mpr + Fr)) * E;
847 sumB += (-1410 * sin(Mr + Mpr - Fr)) * E;
848 sumB += (-1344 * sin(Mr - Fr)) * E;
849 sumB += -1335 * sin(Dr - Fr);
850 sumB += 1107 * sin(3.0 * Mpr + Fr);
851 sumB += 1021 * sin(4.0 * Dr - Fr);
852 sumB += 833 * sin(4.0 * Dr - Mpr + Fr);
853 sumB += 777 * sin(Mpr - 3.0 * Fr);
854 sumB += 671 * sin(4.0 * Dr - 2.0 * Mpr + Fr);
855 sumB += 607 * sin(2.0 * Dr - 3.0 * Fr);
856 sumB += 596 * sin(2.0 * Dr + 2.0 * Mpr - Fr);
857 sumB += (491 * sin(2.0 * Dr - Mr + Mpr - Fr)) * E;
858 sumB += -451 * sin(2.0 * Dr - 2.0 * Mpr + Fr);
859 sumB += 439 * sin(3.0 * Mpr - Fr);
860 sumB += 422 * sin(2.0 * Dr + 2.0 * Mpr + Fr);
861 sumB += 421 * sin(2.0 * Dr - 3.0 * Mpr - Fr);
862 sumB += (-366 * sin(2.0 * Dr + Mr - Mpr + Fr)) * E;
863 sumB += (-351 * sin(2.0 * Dr + Mr + Fr)) * E;
864 sumB += 331 * sin(4.0 * Dr + Fr);
865 sumB += (315 * sin(2.0 * Dr - Mr + Mpr + Fr)) * E;
866 sumB += (302 * sin(2.0 * Dr - 2.0 * Mr - Fr)) * E * E;
867 sumB += -283 * sin(Mpr + 3.0 * Fr);
868 sumB += (-229 * sin(2.0 * Dr + Mr + Mpr - Fr)) * E;
869 sumB += (223 * sin(Dr + Mr - Fr)) * E;
870 sumB += (223 * sin(Dr + Mr + Fr)) * E;
871 sumB += (-220 * sin(Mr - 2.0 * Mpr - Fr)) * E;
872 sumB += (-220 * sin(2.0 * Dr + Mr - Mpr - Fr)) * E;
873 sumB += -185 * sin(Dr + Mpr + Fr);
874 sumB += (181 * sin(2.0 * Dr - Mr - 2.0 * Mpr - Fr)) * E;
875 sumB += (-177 * sin(Mr + 2.0 * Mpr + Fr)) * E;
876 sumB += 176 * sin(4.0 * Dr - 2.0 * Mpr - Fr);
877 sumB += (166 * sin(4.0 * Dr - Mr - Mpr - Fr)) * E;
878 sumB += -164 * sin(Dr + Mpr - Fr);
879 sumB += 132 * sin(4.0 * Dr + Mpr - Fr);
880 sumB += -119 * sin(Dr - Mpr - Fr);
881 sumB += (115 * sin(4.0 * Dr - Mr - Fr)) * E;
882 sumB += (107 * sin(2.0 * Dr - 2.0 * Mr + Fr)) * E;
883 sumB += (-2235 * sin(Lpr));
884 sumB += 382 * sin(A3);
885 sumB += 175 * sin(A1 - Fr);
886 sumB += 175 * sin(A1 + Fr);
887 sumB += 127 * sin(Lpr - Mpr);
888 sumB += -115 * sin(Lpr + Mpr);
889 *Lat = sumB / 1E6;
890 while (*Lat > 360.0)
891 *Lat -= 360.0;
892 while (*Lat < -360.0)
893 *Lat += 360.0;
894
895 sumR += -20905355 * cos(Mpr);
896 sumR += -3699111 * cos(2.0 * Dr - Mpr);
897 sumR += -2955968 * cos(2.0 * Dr);
898 sumR += -569925 * cos(2.0 * Mpr);
899 sumR += (48888 * cos(Mr)) * E;
900 sumR += -3149 * cos(2.0 * Fr);
901 sumR += 246158 * cos(2.0 * Dr - 2.0 * Mpr);
902 sumR += (-152138 * cos(2.0 * Dr - Mr - Mpr)) * E;
903 sumR += -170733 * cos(2.0 * Dr + Mpr);
904 sumR += (-204586 * cos(2.0 * Dr - Mr)) * E;
905 sumR += (-129620 * cos(Mr - Mpr)) * E;
906 sumR += 108743 * cos(Dr);
907 sumR += (104755 * cos(Mr + Mpr)) * E;
908 sumR += 10321 * cos(2.0 * Dr - 2.0 * Fr);
909 sumR += 79661 * cos(Mpr - 2.0 * Fr);
910 sumR += -34782 * cos(4.0 * Dr - Mpr);
911 sumR += -23210 * cos(3.0 * Mpr);
912 sumR += -21636 * cos(4.0 * Dr - 2.0 * Mpr);
913 sumR += (24208 * cos(2.0 * Dr + Mr - Mpr)) * E;
914 sumR += (30824 * cos(2.0 * Dr + Mr)) * E;
915 sumR += -8379 * cos(Dr - Mpr);
916 sumR += (-16675 * cos(Dr + Mr)) * E;
917 sumR += (-12831 * cos(2.0 * Dr - Mr + Mpr)) * E;
918 sumR += -10445 * cos(2.0 * Dr + 2.0 * Mpr);
919 sumR += -11650 * cos(4.0 * Dr);
920 sumR += 14403 * cos(2.0 * Dr - 3.0 * Mpr);
921 sumR += (-7003 * cos(Mr - 2.0 * Mpr)) * E;
922 sumR += (10056 * cos(2.0 * Dr - Mr - 2.0 * Mpr)) * E;
923 sumR += 6322 * cos(Dr + Mpr);
924 sumR += (-9884 * cos(2.0 * Dr - 2.0 * Mr)) * E * E;
925 sumR += (5751 * cos(Mr + 2.0 * Mpr)) * E;
926 sumR += (-4950 * cos(2.0 * Dr - 2.0 * Mr - Mpr)) * E * E;
927 sumR += 4130 * cos(2.0 * Dr + Mpr - 2.0 * Fr);
928 sumR += (-3958 * cos(4.0 * Dr - Mr - Mpr)) * E;
929 sumR += 3258 * cos(3.0 * Dr - Mpr);
930 sumR += (2616 * cos(2.0 * Dr + Mr + Mpr)) * E;
931 sumR += (-1897 * cos(4.0 * Dr - Mr - 2.0 * Mpr)) * E;
932 sumR += (-2117 * cos(2.0 * Mr - Mpr)) * E * E;
933 sumR += (2354 * cos(2.0 * Dr + 2.0 * Mr - Mpr)) * E * E;
934 sumR += -1423 * cos(4.0 * Dr + Mpr);
935 sumR += -1117 * cos(4.0 * Mpr);
936 sumR += (-1571 * cos(4.0 * Dr - Mr)) * E;
937 sumR += -1739 * cos(Dr - 2.0 * Mpr);
938 sumR += -4421 * cos(2.0 * Mpr - 2.0 * Fr);
939 sumR += (1165 * cos(2.0 * Mr + Mpr)) * E * E;
940 sumR += 8752 * cos(2.0 * Dr - Mpr - 2.0 * Fr);
941 *Rad = 385000.56 + (sumR / 1000.0); /* Center of Earth to center of Moon in kilometers */
942 *Rad /= 6378.0; /* Return distance in Earth radii */
943 *Long = DEG2RAD(*Long);
944 *Lat = DEG2RAD(*Lat);
945 }
946