1 #include <cmath>
2 #include <cstdio>
3 #include <cctype>
4 #include <cstring>
5 #include <sstream>
6 #include <string>
7 using namespace std;
8
9 #include "findFile.h"
10 #include "Options.h"
11 #include "Planet.h"
12 #include "xpUtil.h"
13
14 #include "libephemeris/ephemerisWrapper.h"
15
16 body
parseBodyName(char * name)17 Planet::parseBodyName(char *name)
18 {
19 body return_body;
20 char *lowercase = name;
21 char *ptr = name;
22 while (*ptr != '\0') *ptr++ = tolower(*name++);
23
24 if (strncmp(lowercase, "above", 2) == 0)
25 return_body = ABOVE_ORBIT;
26 else if (strncmp(lowercase, body_string[ARIEL], 2) == 0)
27 return_body = ARIEL;
28 else if (strncmp(lowercase, "below", 1) == 0)
29 return_body = BELOW_ORBIT;
30 else if (strncmp(lowercase, body_string[CALLISTO], 2) == 0)
31 return_body = CALLISTO;
32 else if (strncmp(lowercase, body_string[CHARON], 2) == 0)
33 return_body = CHARON;
34 else if (strncmp(lowercase, "default", 3) == 0)
35 return_body = DEFAULT;
36 else if (strncmp(lowercase, body_string[DEIMOS], 3) == 0)
37 return_body = DEIMOS;
38 else if (strncmp(lowercase, body_string[DIONE], 2) == 0)
39 return_body = DIONE;
40 else if (strncmp(lowercase, body_string[EARTH], 2) == 0)
41 return_body = EARTH;
42 else if (strncmp(lowercase, body_string[ENCELADUS], 2) == 0)
43 return_body = ENCELADUS;
44 else if (strncmp(lowercase, body_string[EUROPA], 2) == 0)
45 return_body = EUROPA;
46 else if (strncmp(lowercase, body_string[GANYMEDE], 1) == 0)
47 return_body = GANYMEDE;
48 else if (strncmp(lowercase, body_string[HYPERION], 1) == 0)
49 return_body = HYPERION;
50 else if (strncmp(lowercase, body_string[IAPETUS], 2) == 0)
51 return_body = IAPETUS;
52 else if (strncmp(lowercase, body_string[IO], 2) == 0)
53 return_body = IO;
54 else if (strncmp(lowercase, body_string[JUPITER], 1) == 0)
55 return_body = JUPITER;
56 else if (strncmp(lowercase, "major", 3) == 0)
57 return_body = MAJOR_PLANET;
58 else if (strncmp(lowercase, body_string[MARS], 3) == 0)
59 return_body = MARS;
60 else if (strncmp(lowercase, body_string[MERCURY], 2) == 0)
61 return_body = MERCURY;
62 else if (strncmp(lowercase, body_string[MIMAS], 3) == 0)
63 return_body = MIMAS;
64 else if (strncmp(lowercase, body_string[MIRANDA], 3) == 0)
65 return_body = MIRANDA;
66 else if (strncmp(lowercase, body_string[MOON], 2) == 0)
67 return_body = MOON;
68 else if (strncmp(lowercase, "naif", 2) == 0)
69 return_body = NAIF;
70 else if (strncmp(lowercase, body_string[NEPTUNE], 3) == 0)
71 return_body = NEPTUNE;
72 else if (strncmp(lowercase, body_string[NEREID], 3) == 0)
73 return_body = NEREID;
74 else if (strncmp(lowercase, "norad", 2) == 0)
75 return_body = NORAD;
76 else if (strncmp(lowercase, body_string[OBERON], 1) == 0)
77 return_body = OBERON;
78 else if (strncmp(lowercase, "path", 2) == 0)
79 return_body = ALONG_PATH;
80 else if (strncmp(lowercase, body_string[PHOBOS], 4) == 0)
81 return_body = PHOBOS;
82 else if (strncmp(lowercase, body_string[PHOEBE], 4) == 0)
83 return_body = PHOEBE;
84 else if (strncmp(lowercase, body_string[PLUTO], 2) == 0)
85 return_body = PLUTO;
86 else if (strncmp(lowercase, "random", 2) == 0)
87 return_body = RANDOM_BODY;
88 else if (strncmp(lowercase, body_string[RHEA], 2) == 0)
89 return_body = RHEA;
90 else if (strncmp(lowercase, body_string[SATURN], 2) == 0)
91 return_body = SATURN;
92 else if (strncmp(lowercase, body_string[SUN], 2) == 0)
93 return_body = SUN;
94 else if (strncmp(lowercase, "system", 2) == 0)
95 return_body = SAME_SYSTEM;
96 else if (strncmp(lowercase, body_string[TETHYS], 2) == 0)
97 return_body = TETHYS;
98 else if (strncmp(lowercase, body_string[TITANIA], 6) == 0)
99 return_body = TITANIA;
100 else if (strncmp(lowercase, body_string[TITAN], 5) == 0)
101 return_body = TITAN;
102 else if (strncmp(lowercase, body_string[TRITON], 2) == 0)
103 return_body = TRITON;
104 else if (strncmp(lowercase, body_string[URANUS], 2) == 0)
105 return_body = URANUS;
106 else if (strncmp(lowercase, body_string[UMBRIEL], 2) == 0)
107 return_body = UMBRIEL;
108 else if (strncmp(lowercase, body_string[VENUS], 1) == 0)
109 return_body = VENUS;
110 else
111 {
112 xpWarn("parseBody: Unknown body specified\n",
113 __FILE__, __LINE__);
114 return_body = UNKNOWN_BODY;
115 }
116 return(return_body);
117 }
118
119 /*
120 Rotational parameters are from Seidelmann et al. (2002), Celestial
121 Mechanics 82, 83--110.
122 */
Planet(const double jd,const body this_body)123 Planet::Planet(const double jd, const body this_body)
124 : index_(this_body),
125 julianDay_(jd),
126 needRotationMatrix_(true),
127 period_(0),
128 needShadowCoeffs_(true)
129 {
130 Options *options = Options::getInstance();
131 if (options->UniversalTime())
132 julianDay_ += delT(julianDay_) / 86400;
133
134 d2000_ = (julianDay_ - 2451545.0);
135 T2000_ = d2000_ / 36525;
136
137 switch (index_)
138 {
139 case SUN:
140 primary_ = SUN;
141
142 alpha0_ = 286.13;
143 delta0_ = 63.87;
144
145 nullMeridian0_ = 84.10;
146 wdot_ = 14.1844;
147
148 flipped_ = 1;
149
150 period_ = 0;
151 radiusEq_ = 696000;
152 flattening_ = 0;
153
154 break;
155 case MERCURY:
156 primary_ = SUN;
157
158 alpha0_ = 281.01 - 0.033 * T2000_;
159 delta0_ = 61.45 - 0.005 * T2000_;
160
161 nullMeridian0_ = 329.548;
162 wdot_ = 6.1385025;
163
164 flipped_ = -1;
165
166 period_ = 87.969;
167
168 radiusEq_ = 2439;
169 flattening_ = 0;
170
171 break;
172 case VENUS:
173 primary_ = SUN;
174
175 alpha0_ = 272.76;
176 delta0_ = 67.16;
177
178 nullMeridian0_ = 160.20;
179 wdot_ = -1.4813688;
180
181 flipped_ = 1;
182
183 period_ = 224.701;
184
185 radiusEq_ = 6051.9;
186 flattening_ = 0;
187
188 break;
189 case EARTH:
190 primary_ = SUN;
191
192 alpha0_ = 0 - 0.641 * T2000_;
193 delta0_ = 90 - .557 * T2000_;
194
195 nullMeridian0_ = 190.147;
196 wdot_ = 360.9856235;
197
198 flipped_ = 1;
199
200 period_ = 365.256;
201
202 // WGS84 ellipsoid
203 radiusEq_ = 6378.137;
204 flattening_ = 1/298.257223563;
205
206 break;
207 case MOON:
208 {
209 primary_ = EARTH;
210
211 const double E1 = (125.045 - 0.0529921 * d2000_) * deg_to_rad;
212 const double E2 = (250.089 - 0.1059842 * d2000_) * deg_to_rad;
213 const double E3 = (260.008 + 13.0120009 * d2000_) * deg_to_rad;
214 const double E4 = (176.625 + 13.3407154 * d2000_) * deg_to_rad;
215 const double E5 = (357.529 + 0.9856003 * d2000_) * deg_to_rad;
216 const double E6 = (311.589 + 26.4057084 * d2000_) * deg_to_rad;
217 const double E7 = (134.963 + 13.0649930 * d2000_) * deg_to_rad;
218 const double E8 = (276.617 + 0.3287146 * d2000_) * deg_to_rad;
219 const double E9 = ( 34.226 + 1.7484877 * d2000_) * deg_to_rad;
220 const double E10 = ( 15.134 - 0.1589763 * d2000_) * deg_to_rad;
221 const double E11 = (119.743 + 0.0036096 * d2000_) * deg_to_rad;
222 const double E12 = (239.961 + 0.1643573 * d2000_) * deg_to_rad;
223 const double E13 = ( 25.053 + 12.9590088 * d2000_) * deg_to_rad;
224
225 alpha0_ = 269.9949 + (0.0031 * T2000_ - 3.8787 * sin(E1)
226 - 0.1204 * sin(E2) + 0.0700 * sin(E3)
227 - 0.0172 * sin(E4) + 0.0072 * sin(E6)
228 - 0.0052 * sin(E10) + 0.0043 * sin(E13));
229 delta0_ = 66.5392 + (0.0130 * T2000_ + 1.5419 * cos(E1)
230 + 0.0239 * cos(E2) - 0.0278 * cos(E3)
231 + 0.0068 * cos(E4) - 0.0029 * cos(E6)
232 + 0.0009 * cos(E7) + 0.0008 * cos(E10)
233 - 0.0009 * cos(E13));
234 nullMeridian0_ = 38.3213 + (d2000_ * (13.17635815 - 1.4E-12 * d2000_)
235 + 3.5610 * sin(E1)
236 + 0.1208 * sin(E2)
237 - 0.0642 * sin(E3)
238 + 0.0158 * sin(E4)
239 + 0.0252 * sin(E5)
240 - 0.0066 * sin(E6)
241 - 0.0047 * sin(E7)
242 - 0.0046 * sin(E8)
243 + 0.0028 * sin(E9)
244 + 0.0052 * sin(E10)
245 + 0.0040 * sin(E11)
246 + 0.0019 * sin(E12)
247 - 0.0044 * sin(E13));
248 wdot_ = 0;
249
250 flipped_ = 1;
251
252 period_ = 27.321661;
253
254 radiusEq_ = 1737.5;
255 flattening_ = 0.002;
256
257 }
258 break;
259 case MARS:
260 case PHOBOS:
261 case DEIMOS:
262 {
263 flipped_ = -1;
264
265 const double M1 = (169.51 - 0.4357640 * d2000_) * deg_to_rad;
266 const double M2 = (192.93 + 1128.4096700 * d2000_) * deg_to_rad;
267 const double M3 = (53.47 - 0.0181510 * d2000_) * deg_to_rad;
268
269 switch (index_)
270 {
271 case MARS:
272 primary_ = SUN;
273
274 alpha0_ = 317.68143 - 0.1061 * T2000_;
275 delta0_ = 52.88650 - 0.0609 * T2000_;
276
277 nullMeridian0_ = 176.630;
278 wdot_ = 350.89198226;
279
280 period_ = 686.980;
281
282 radiusEq_ = 3397;
283 flattening_ = 0.0065;
284
285 break;
286 case PHOBOS:
287 primary_ = MARS;
288
289 alpha0_ = 317.68 - 0.108 * T2000_ + 1.79 * sin(M1);
290 delta0_ = 52.90 - 0.061 * T2000_ - 1.08 * cos(M1);
291
292 nullMeridian0_ = (35.06 + 8.864 * T2000_ * T2000_ - 1.42 * sin(M1)
293 - 0.78 * sin(M2));
294 wdot_ = 1128.8445850;
295
296 period_ = 0.31891023;
297
298 radiusEq_ = 10; // non-spherical
299 flattening_ = 0;
300
301 break;
302 case DEIMOS:
303 primary_ = MARS;
304
305 alpha0_ = 316.65 - 0.108 * T2000_ + 2.98 * sin(M3);
306 delta0_ = 53.52 - 0.061 * T2000_ - 1.78 * cos(M3);
307
308 nullMeridian0_ = (79.41 - 0.520 * T2000_ * T2000_
309 - 2.58 * sin(M3) + 0.19 * cos(M3));
310 wdot_ = 285.1618970;
311
312 period_ = 1.2624407;
313 radiusEq_ = 6; // non-spherical
314 flattening_ = 0;
315
316 break;
317 default:
318 break;
319 }
320 }
321 break;
322 case JUPITER:
323 case IO:
324 case EUROPA:
325 case GANYMEDE:
326 case CALLISTO:
327 {
328 flipped_ = -1;
329
330 const double J3 = (283.90 + 4850.7 * T2000_) * deg_to_rad;
331 const double J4 = (355.80 + 1191.3 * T2000_) * deg_to_rad;
332 const double J5 = (119.90 + 262.1 * T2000_) * deg_to_rad;
333 const double J6 = (229.80 + 64.3 * T2000_) * deg_to_rad;
334 const double J7 = (352.25 + 2382.6 * T2000_) * deg_to_rad;
335 const double J8 = (113.35 + 6070.0 * T2000_) * deg_to_rad;
336
337 switch (index_)
338 {
339 case JUPITER:
340 primary_ = SUN;
341
342 alpha0_ = 268.05 - 0.009 * T2000_;
343 delta0_ = 64.49 + 0.003 * T2000_;
344
345 // System III (magnetic field rotation)
346 nullMeridian0_ = 284.95;
347 wdot_ = 870.5366420;
348
349 if (options->GRSSet())
350 {
351 // System II
352 nullMeridian0_ = 43.3;
353 wdot_ = 870.270;
354 }
355
356 period_ = 4332.589;
357
358 radiusEq_ = 71492;
359 flattening_ = 0.06487;
360
361 break;
362 case IO:
363 primary_ = JUPITER;
364
365 alpha0_ = (268.05 - 0.009 * T2000_ + 0.094 * sin(J3)
366 + 0.024 * sin(J4));
367 delta0_ = (64.50 + 0.003 * T2000_ + 0.040 * cos(J3)
368 + 0.011 * cos(J4));
369 nullMeridian0_ = 200.39 - 0.085 * sin(J3) - 0.022 * sin(J4);
370 wdot_ = 203.4889538;
371
372 period_ = 1.769137786;
373
374 radiusEq_ = 1821.6;
375 flattening_ = 0;
376
377 break;
378 case EUROPA:
379 primary_ = JUPITER;
380
381 alpha0_ = (268.08 - 0.009 * T2000_ + 1.086 * sin(J4)
382 + 0.060 * sin(J5) + 0.015 * sin(J6)
383 + 0.009 * sin(J7));
384 delta0_ = (64.51 + 0.003 * T2000_ + 0.468 * cos(J4)
385 + 0.026 * cos(J5) + 0.007 * cos(J6)
386 + 0.002 * cos(J7));
387 nullMeridian0_ = (35.67 - 0.980 * sin(J4) - 0.054 * sin(J5)
388 - 0.014 * sin(J6) - 0.008 * sin(J7));
389 wdot_ = 101.3747235;
390
391 period_ = 3.551181041;
392
393 radiusEq_ = 1560.8;
394 flattening_ = 0;
395
396 break;
397 case GANYMEDE:
398 primary_ = JUPITER;
399
400 alpha0_ = (268.20 - 0.009 * T2000_ - 0.037 * sin(J4)
401 + 0.431 * sin(J5) + 0.091 * sin(J6));
402 delta0_ = (64.57 + 0.003 * T2000_ - 0.016 * cos(J4)
403 + 0.186 * cos(J5) + 0.039 * cos(J6));
404 nullMeridian0_ = (44.04 + 0.033 * sin(J4) - 0.389 * sin(J5)
405 - 0.082 * sin(J6));
406 wdot_ = 50.3176081;
407
408 period_ = 7.15455296;
409
410 radiusEq_ = 2631.2;
411 flattening_ = 0;
412
413 break;
414 case CALLISTO:
415 primary_ = JUPITER;
416
417 alpha0_ = (268.72 - 0.009 * T2000_ - 0.068 * sin(J5)
418 + 0.590 * sin(J6) + 0.010 * sin(J8));
419 delta0_ = (64.83 + 0.003 * T2000_ - 0.029 * cos(J5)
420 + 0.254 * cos(J6) - 0.004 * cos(J8));
421 nullMeridian0_ = (259.73 + 0.061 * sin(J5) - 0.533 * sin(J6)
422 - 0.009 * sin(J8));
423 wdot_ = 21.5710715;
424
425 period_ = 16.6890184;
426
427 radiusEq_ = 2410.3;
428 flattening_ = 0;
429
430 break;
431 default:
432 break;
433 }
434 }
435 break;
436 case SATURN:
437 case MIMAS:
438 case ENCELADUS:
439 case TETHYS:
440 case DIONE:
441 case RHEA:
442 case TITAN:
443 case HYPERION:
444 case IAPETUS:
445 case PHOEBE:
446 {
447 flipped_ = -1;
448
449 const double S3 = (177.40 - 36505.5 * T2000_) * deg_to_rad;
450 const double S4 = (300.00 - 7225.9 * T2000_) * deg_to_rad;
451 const double S5 = (316.45 + 506.2 * T2000_) * deg_to_rad;
452 const double S6 = (345.20 - 1016.3 * T2000_) * deg_to_rad;
453 const double S7 = ( 29.80 - 52.1 * T2000_) * deg_to_rad;
454 switch (index_)
455 {
456 case SATURN:
457 primary_ = SUN;
458
459 alpha0_ = 40.589 - 0.036 * T2000_;
460 delta0_ = 83.537 - 0.004 * T2000_;
461 nullMeridian0_ = 38.90;
462 wdot_ = 810.7939024;
463
464 period_ = 10759.22;
465
466 radiusEq_ = 60268;
467 flattening_ = 0.09796;
468
469 break;
470 case MIMAS:
471 primary_ = SATURN;
472
473 alpha0_ = 40.66 - 0.036 * T2000_ + 13.56 * sin(S3);
474 delta0_ = 83.52 - 0.004 * T2000_ - 1.53 * cos(S3);
475 nullMeridian0_ = 337.46 - 13.48 * sin(S3) - 44.85 * sin(S5);
476 wdot_ = 381.994550;
477
478 period_ = 0.942421813;
479
480 radiusEq_ = 198.6;
481 flattening_ = 0;
482
483 break;
484 case ENCELADUS:
485 primary_ = SATURN;
486
487 alpha0_ = 40.66 - 0.036 * T2000_;
488 delta0_ = 83.52 - 0.004 * T2000_;
489
490 nullMeridian0_ = 2.82;
491 wdot_ = 262.7318996;
492
493 period_ = 1.370217855;
494
495 radiusEq_ = 249.4;
496 flattening_ = 0;
497
498 break;
499 case TETHYS:
500 primary_ = SATURN;
501
502 alpha0_ = 40.66 - 0.036 * T2000_ + 9.66 * sin(S4);
503 delta0_ = 83.52 - 0.004 * T2000_ - 1.09 * cos(S4);
504
505 nullMeridian0_ = 10.45 - 9.60 * sin(S4) + 2.23 * sin(S5);
506 wdot_ = 190.6979085;
507
508 period_ = 1.887802160;
509
510 radiusEq_ = 529.8;
511 flattening_ = 0;
512
513 break;
514 case DIONE:
515 primary_ = SATURN;
516
517 alpha0_ = 40.66 - 0.036 * T2000_;
518 delta0_ = 83.52 - 0.004 * T2000_;
519
520 nullMeridian0_ = 357.00;
521 wdot_ = 131.5349316;
522
523 period_ = 2.736914742;
524
525 radiusEq_ = 559;
526 flattening_ = 0;
527
528 break;
529 case RHEA:
530 primary_ = SATURN;
531
532 alpha0_ = 40.38 - 0.036 * T2000_ + 3.10 * sin(S6);
533 delta0_ = 83.55 - 0.004 * T2000_ - 0.35 * cos(S6);
534 nullMeridian0_ = 235.16 - 3.08 * sin(S6);
535 wdot_ = 79.6900478;
536
537 period_ = 4.517500436;
538
539 radiusEq_ = 764;
540 flattening_ = 0;
541
542 break;
543 case TITAN:
544 primary_ = SATURN;
545
546 alpha0_ = 36.41 - 0.036 * T2000_ + 2.66 * sin(S7);
547 delta0_ = 83.94 - 0.004 * T2000_ - 0.30 * cos(S7);
548 nullMeridian0_ = 189.64 - 2.64 * sin(S7);
549 wdot_ = 22.5769768;
550
551 period_ = 15.94542068;
552
553 radiusEq_ = 2575;
554 flattening_ = 0;
555
556 break;
557 case HYPERION:
558 primary_ = SATURN;
559
560 alpha0_ = 40.589 - 0.036 * T2000_;
561 delta0_ = 83.537 - 0.004 * T2000_;
562
563 nullMeridian0_ = 0;
564 wdot_ = 0;
565
566 period_ = 21.2766088;
567
568 radiusEq_ = 141.5; // non-spherical
569 flattening_ = 0;
570
571 break;
572 case IAPETUS:
573 primary_ = SATURN;
574
575 alpha0_ = 318.16 - 3.949 * T2000_;
576 delta0_ = 75.03 - 1.143 * T2000_;
577 nullMeridian0_ = 350.20;
578 wdot_ = 4.5379572;
579
580 period_ = 79.3301825;
581
582 radiusEq_ = 718;
583 flattening_ = 0;
584
585 break;
586 case PHOEBE:
587 primary_ = SATURN;
588
589 alpha0_ = 355.00;
590 delta0_ = 68.70;
591
592 nullMeridian0_ = 304.70;
593 wdot_ = 930.8338720;
594
595 period_ = 550.48;
596
597 radiusEq_ = 110;
598 flattening_ = 0;
599
600 break;
601 default:
602 break;
603 }
604 }
605 break;
606 case URANUS:
607 case MIRANDA:
608 case ARIEL:
609 case UMBRIEL:
610 case TITANIA:
611 case OBERON:
612 {
613 flipped_ = 1;
614
615 const double U11 = (141.69 + 41887.66 * T2000_) * deg_to_rad;
616 const double U12 = (316.41 + 2863.96 * T2000_) * deg_to_rad;
617 const double U13 = (304.01 - 51.94 * T2000_) * deg_to_rad;
618 const double U14 = (308.71 - 93.17 * T2000_) * deg_to_rad;
619 const double U15 = (340.82 - 75.32 * T2000_) * deg_to_rad;
620 const double U16 = (259.14 - 504.81 * T2000_) * deg_to_rad;
621
622 switch (index_)
623 {
624 case URANUS:
625 primary_ = SUN;
626
627 alpha0_ = 257.311;
628 delta0_ = -15.175;
629 nullMeridian0_ = 203.81;
630 wdot_ = -501.1600928;
631
632 period_ = 30685.4;
633
634 radiusEq_ = 25559;
635 flattening_ = 0.02293;
636
637 break;
638 case MIRANDA:
639 primary_ = URANUS;
640
641 alpha0_ = 257.42 + 4.41 * sin(U11) - 0.04 * sin(2*U11);
642 delta0_ = -15.08 + 4.25 * cos(U11) - 0.02 * cos(2*U11);
643 nullMeridian0_ = (30.70 - 1.27 * sin(U12) + 0.15 * sin(2*U12)
644 + 1.15 * sin(U11) - 0.09 * sin(2*U11));
645 wdot_ = -254.6906892;
646
647 period_ = 1.41347925;
648
649 radiusEq_ = 235.8;
650 flattening_ = 0;
651
652 break;
653 case ARIEL:
654 primary_ = URANUS;
655
656 alpha0_ = 257.43 + 0.29 * sin(U13);
657 delta0_ = -15.10 + 0.28 * cos(U13);
658 nullMeridian0_ = 156.22 + 0.05 * sin(U12) + 0.08 * sin(U13);
659 wdot_ = -142.8356681;
660
661 period_ = 2.52037935;
662
663 radiusEq_ = 578.9;
664 flattening_ = 0;
665
666 break;
667 case UMBRIEL:
668 primary_ = URANUS;
669
670 alpha0_ = 257.43 + 0.21 * sin(U14);
671 delta0_ = -15.10 + 0.20 * cos(U14);
672 nullMeridian0_ = 108.05 - 0.09 * sin(U12) + 0.06 * sin(U14);
673 wdot_ = -86.8688923;
674
675 period_ = 4.1441772;
676
677 radiusEq_ = 584.7;
678 flattening_ = 0;
679
680 break;
681 case TITANIA:
682 primary_ = URANUS;
683
684 alpha0_ = 257.43 + 0.29 * sin(U15);
685 delta0_ = -15.10 + 0.28 * cos(U15);
686 nullMeridian0_ = 77.74 + 0.08 * sin(U15);
687 wdot_ = -41.3514316;
688
689 period_ = 8.7058717;
690
691 radiusEq_ = 788.9;
692 flattening_ = 0;
693
694 break;
695 case OBERON:
696 primary_ = URANUS;
697
698 alpha0_ = 257.43 + 0.16 * sin(U16);
699 delta0_ = -15.10 + 0.16 * cos(U16);
700 nullMeridian0_ = 6.77 + 0.04 * sin(U16);
701 wdot_ = -26.7394932;
702
703 period_ = 13.4632389;
704
705 radiusEq_ = 761.4;
706 flattening_ = 0;
707
708 break;
709 default:
710 break;
711 }
712 }
713 break;
714 case NEPTUNE:
715 case TRITON:
716 case NEREID:
717 {
718 flipped_ = -1;
719
720 switch (index_)
721 {
722 case NEPTUNE:
723 {
724 primary_ = SUN;
725
726 double N = (357.85 + 52.316 * T2000_) * deg_to_rad;
727 alpha0_ = 299.36 + 0.70 * sin(N);
728 delta0_ = 43.46 - 0.51 * cos(N);
729
730 nullMeridian0_ = 253.18 - 0.48 * sin(N);
731 wdot_ = 536.3128492;
732
733 period_ = 60189.0;
734
735 radiusEq_ = 24764;
736 flattening_ = 0.0171;
737
738 }
739 break;
740 case TRITON:
741 {
742 primary_ = NEPTUNE;
743
744 const double N7 = (177.85 + 52.316 * T2000_) * deg_to_rad;
745 alpha0_ = (299.36 - 32.35 * sin(N7) - 6.28 * sin(2*N7)
746 - 2.08 * sin(3*N7) - 0.74 * sin(4*N7)
747 - 0.28 * sin(5*N7) - 0.11 * sin(6*N7)
748 - 0.07 * sin(7*N7) - 0.02 * sin(8*N7)
749 - 0.01 * sin(9*N7));
750 delta0_ = (43.46 + 22.55 * cos(N7) + 2.10 * cos(2*N7)
751 + 0.55 * cos(3*N7) + 0.16 * cos(4*N7)
752 + 0.05 * cos(5*N7) + 0.02 * cos(6*N7)
753 + 0.01 * cos(7*N7));
754
755 nullMeridian0_ = (296.53 + 22.25 * sin(N7) + 6.73 * sin(2*N7)
756 + 2.05 * sin(3*N7) + 0.74 * sin(4*N7)
757 + 0.28 * sin(5*N7) + 0.11 * sin(6*N7)
758 + 0.05 * sin(7*N7) + 0.02 * sin(8*N7)
759 + 0.01 * sin(9*N7));
760
761 wdot_ = -61.2572637;
762
763 period_ = 5.8768541;
764
765 radiusEq_ = 1353;
766 flattening_ = 0;
767 }
768 break;
769 case NEREID:
770 {
771 primary_ = NEPTUNE;
772
773 alpha0_ = 299.36;
774 delta0_ = 43.46;
775
776 nullMeridian0_ = 0;
777 wdot_ = 0;
778
779 period_ = 360.13619;
780
781 radiusEq_ = 170;
782 flattening_ = 0;
783 }
784 break;
785 default:
786 break;
787 }
788 }
789 break;
790 case PLUTO:
791 case CHARON:
792 {
793 flipped_ = 1;
794
795 alpha0_ = 313.02;
796 delta0_ = 9.09;
797
798 wdot_ = -56.3623195;
799
800 switch (index_)
801 {
802 case PLUTO:
803 primary_ = SUN;
804
805 nullMeridian0_ = 236.77;
806
807 period_ = 90465.0;
808
809 radiusEq_ = 1151;
810 flattening_ = 0;
811 break;
812 case CHARON:
813 primary_ = PLUTO;
814
815 nullMeridian0_ = 56.77;
816
817 period_ = 6.38723;
818
819 radiusEq_ = 593;
820 flattening_ = 0;
821 break;
822 default:
823 break;
824 }
825 }
826 break;
827 default:
828 xpExit("Planet: Unknown body specified.\n", __FILE__, __LINE__);
829 }
830
831 alpha0_ *= deg_to_rad;
832 delta0_ *= deg_to_rad;
833
834 radiusEq_ /= AU_to_km;
835
836 nullMeridian_ = fmod(nullMeridian0_ + wdot_ * d2000_, 360.0);
837 nullMeridian_ *= deg_to_rad;
838
839 radiusPol_ = radiusEq_ * (1 - flattening_);
840 omf2_ = (1 - flattening_) * (1 - flattening_);
841 }
842
~Planet()843 Planet::~Planet()
844 {
845 }
846
847 // Return the distance from the surface to the center, in units of
848 // equatorial radius. Input latitude is assumed to be planetographic.
849 double
Radius(const double lat) const850 Planet::Radius(const double lat) const
851 {
852 double returnValue = 1;
853 if (index_ == JUPITER || index_ == SATURN)
854 {
855 double tmpLat = lat;
856 double tmpLon = 0;
857 PlanetographicToPlanetocentric(tmpLat, tmpLon);
858
859 const double ReSinLat = radiusEq_ * sin(tmpLat);
860 const double RpCosLat = radiusPol_ * cos(tmpLat);
861
862 // from http://mathworld.wolfram.com/Ellipse.html
863 returnValue = radiusPol_ / sqrt(RpCosLat * RpCosLat
864 + ReSinLat * ReSinLat);
865 }
866
867 return(returnValue);
868 }
869
870 // Compute heliocentric equatorial coordinates of the planet
871 void
calcHeliocentricEquatorial()872 Planet::calcHeliocentricEquatorial()
873 {
874 calcHeliocentricEquatorial(true);
875 }
876
877 void
calcHeliocentricEquatorial(const bool relativeToSun)878 Planet::calcHeliocentricEquatorial(const bool relativeToSun)
879 {
880 GetHeliocentricXYZ(index_, primary_, julianDay_,
881 relativeToSun, X_, Y_, Z_);
882 }
883
884 // Return rectangular coordinates in heliocentric equatorial frame
885 void
getPosition(double & X,double & Y,double & Z) const886 Planet::getPosition(double &X, double &Y, double &Z) const
887 {
888 X = X_;
889 Y = Y_;
890 Z = Z_;
891 }
892
893 double
Illumination(const double oX,const double oY,const double oZ)894 Planet::Illumination(const double oX, const double oY, const double oZ)
895 {
896 return(50 * (ndot(X_, Y_, Z_, X_ - oX, Y_ - oY, Z_ - oZ) + 1));
897 }
898
899 void
CreateRotationMatrix()900 Planet::CreateRotationMatrix()
901 {
902 /*
903 These equations are from "Astronomy on the Personal Computer,
904 3rd Edition", by Montenbruck and Pfleger, Chapter 7.
905 */
906 const double cosW = cos(nullMeridian_);
907 const double sinW = sin(nullMeridian_);
908
909 const double cosA = cos(alpha0_);
910 const double sinA = sin(alpha0_);
911
912 const double cosD = cos(delta0_);
913 const double sinD = sin(delta0_);
914
915 rot_[0][0] = -cosW * sinA - sinW * sinD * cosA;
916 rot_[0][1] = cosW * cosA - sinW * sinD * sinA;
917 rot_[0][2] = sinW * cosD;
918
919 rot_[1][0] = sinW * sinA - cosW * sinD * cosA;
920 rot_[1][1] = -sinW * cosA - cosW * sinD * sinA;
921 rot_[1][2] = cosW * cosD;
922
923 rot_[2][0] = cosD * cosA;
924 rot_[2][1] = cosD * sinA;
925 rot_[2][2] = sinD;
926
927 invertMatrix(rot_, invRot_);
928
929 needRotationMatrix_ = false;
930 }
931
932 void
PlanetocentricToXYZ(double & X,double & Y,double & Z,const double lat,const double lon,const double rad)933 Planet::PlanetocentricToXYZ(double &X, double &Y, double &Z,
934 const double lat, const double lon,
935 const double rad)
936 {
937 if (needRotationMatrix_) CreateRotationMatrix();
938
939 double r[3];
940 r[0] = cos(lat) * cos(lon);
941 r[1] = cos(lat) * sin(lon);
942 r[2] = sin(lat);
943
944 double newrad = rad * radiusEq_;
945 X = dot(invRot_[0], r) * newrad;
946 Y = dot(invRot_[1], r) * newrad;
947 Z = dot(invRot_[2], r) * newrad;
948
949 X += X_;
950 Y += Y_;
951 Z += Z_;
952 }
953
954 void
PlanetographicToXYZ(double & X,double & Y,double & Z,double lat,double lon,const double rad)955 Planet::PlanetographicToXYZ(double &X, double &Y, double &Z,
956 double lat, double lon, const double rad)
957 {
958 PlanetographicToPlanetocentric(lat, lon);
959 PlanetocentricToXYZ(X, Y, Z, lat, lon, rad);
960 }
961
962 void
XYZToPlanetocentric(const double X,const double Y,const double Z,double & lat,double & lon)963 Planet::XYZToPlanetocentric(const double X, const double Y, const double Z,
964 double &lat, double &lon)
965 {
966 double rad = 0;
967 XYZToPlanetocentric(X, Y, Z, lat, lon, rad);
968 }
969
970 void
XYZToPlanetocentric(const double X,const double Y,const double Z,double & lat,double & lon,double & rad)971 Planet::XYZToPlanetocentric(const double X, const double Y, const double Z,
972 double &lat, double &lon, double &rad)
973 {
974 if (needRotationMatrix_) CreateRotationMatrix();
975
976 const double r[3] = { X - X_, Y - Y_, Z - Z_ };
977
978 const double sx = dot(rot_[0], r);
979 const double sy = dot(rot_[1], r);
980 const double sz = dot(rot_[2], r);
981
982 rad = sqrt(sx * sx + sy*sy);
983 if (rad > 0)
984 lat = atan(sz/rad);
985 else
986 lat = 0;
987
988 if (cos(lat) > 0)
989 lon = atan2(sy, sx);
990 else
991 lon = 0;
992
993 if (lon < 0) lon += TWO_PI;
994
995 rad = sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2]);
996 rad /= radiusEq_;
997 }
998
999 void
XYZToPlanetographic(const double X,const double Y,const double Z,double & lat,double & lon)1000 Planet::XYZToPlanetographic(const double X, const double Y, const double Z,
1001 double &lat, double &lon)
1002 {
1003 double rad = 0;
1004 XYZToPlanetographic(X, Y, Z, lat, lon, rad);
1005 }
1006
1007 void
XYZToPlanetographic(const double X,const double Y,const double Z,double & lat,double & lon,double & rad)1008 Planet::XYZToPlanetographic(const double X, const double Y, const double Z,
1009 double &lat, double &lon, double &rad)
1010 {
1011 XYZToPlanetocentric(X, Y, Z, lat, lon, rad);
1012 PlanetocentricToPlanetographic(lat, lon);
1013 }
1014
1015 void
XYZToPlanetaryXYZ(const double X,const double Y,const double Z,double & pX,double & pY,double & pZ)1016 Planet::XYZToPlanetaryXYZ(const double X, const double Y, const double Z,
1017 double &pX, double &pY, double &pZ)
1018 {
1019 double lat, lon, rad;
1020 XYZToPlanetocentric(X, Y, Z, lat, lon, rad);
1021
1022 pX = rad * cos(lat) * cos(lon);
1023 pY = rad * cos(lat) * sin(lon);
1024 pZ = rad * sin(lat);
1025 }
1026
1027 void
PlanetaryXYZToXYZ(const double pX,const double pY,const double pZ,double & X,double & Y,double & Z)1028 Planet::PlanetaryXYZToXYZ(const double pX, const double pY, const double pZ,
1029 double &X, double &Y, double &Z)
1030 {
1031 const double lat = atan(pZ / sqrt(pX * pX + pY * pY));
1032 double lon;
1033 if (cos(lat) > 1e-5)
1034 lon = atan2(pY, pX);
1035 else
1036 lon = 0;
1037
1038 const double rad = sqrt(pX * pX + pY * pY + pZ * pZ);
1039
1040 PlanetocentricToXYZ(X, Y, Z, lat, lon, rad);
1041 }
1042
1043 void
PlanetocentricToPlanetographic(double & lat,double & lon) const1044 Planet::PlanetocentricToPlanetographic(double &lat, double &lon) const
1045 {
1046 if (flattening_ > 0)
1047 lat = atan(tan(lat) / omf2_);
1048
1049 if (index_ == EARTH || index_ == SUN) lon *= -1;
1050 if (wdot_ > 0) lon *= -1;
1051 if (lon < 0) lon += TWO_PI;
1052 }
1053
1054 void
PlanetographicToPlanetocentric(double & lat,double & lon) const1055 Planet::PlanetographicToPlanetocentric(double &lat, double &lon) const
1056 {
1057 if (flattening_ > 0)
1058 lat = atan(omf2_ * tan(lat));
1059
1060 if (index_ == EARTH || index_ == SUN) lon *= -1;
1061 if (wdot_ > 0) lon *= -1;
1062 }
1063
1064 void
ComputeShadowCoeffs()1065 Planet::ComputeShadowCoeffs()
1066 {
1067 XYZToPlanetaryXYZ(0, 0, 0, sunX_, sunY_, sunZ_);
1068
1069 ellipseCoeffC_ = sunZ_ * sunZ_ / omf2_;
1070 ellipseCoeffC_ += sunY_ * sunY_;
1071 ellipseCoeffC_ += sunX_ * sunX_;
1072 ellipseCoeffC_ -= 1;
1073
1074 needShadowCoeffs_ = false;
1075 }
1076
1077 // x, y, z must be in planetary XYZ frame
1078 bool
IsInMyShadow(const double x,const double y,const double z)1079 Planet::IsInMyShadow(const double x, const double y, const double z)
1080 {
1081 if (needShadowCoeffs_) ComputeShadowCoeffs();
1082
1083 double ellipseCoeffA = (z - sunZ_) * (z - sunZ_) / omf2_;
1084 ellipseCoeffA += (y - sunY_) * (y - sunY_);
1085 ellipseCoeffA += (x - sunX_) * (x - sunX_);
1086
1087 double ellipseCoeffB = (z - sunZ_) * sunZ_ / omf2_;
1088 ellipseCoeffB += (y - sunY_) * sunY_;
1089 ellipseCoeffB += (x - sunX_) * sunX_;
1090
1091 const double determinant = (ellipseCoeffB * ellipseCoeffB
1092 - ellipseCoeffA * ellipseCoeffC_);
1093
1094 return(determinant > 0);
1095 }
1096
1097 void
getOrbitalNorth(double & X,double & Y,double & Z) const1098 Planet::getOrbitalNorth(double &X, double &Y, double &Z) const
1099 {
1100 if (index_ == SUN)
1101 {
1102 X = cos(delta0_) * cos(alpha0_);
1103 Y = cos(delta0_) * sin(alpha0_);
1104 Z = sin(delta0_);
1105 return;
1106 }
1107
1108 // cross product of position and velocity vectors points to the
1109 // orbital north pole
1110 double pos[3] = { X_, Y_, Z_ };
1111
1112 double pX0, pY0, pZ0;
1113 double pX1, pY1, pZ1;
1114
1115 GetHeliocentricXYZ(index_, primary_, julianDay_ - 0.5,
1116 true, pX0, pY0, pZ0);
1117 GetHeliocentricXYZ(index_, primary_, julianDay_ + 0.5,
1118 true, pX1, pY1, pZ1);
1119
1120 double vel[3] = { pX1 - pX0, pY1 - pY0, pZ1 - pZ0 };
1121
1122 double north[3];
1123 cross(pos, vel, north);
1124
1125 double mag = sqrt(dot(north, north));
1126
1127 X = north[0]/mag;
1128 Y = north[1]/mag;
1129 Z = north[2]/mag;
1130 }
1131
1132 void
getBodyNorth(double & X,double & Y,double & Z) const1133 Planet::getBodyNorth(double &X, double &Y, double &Z) const
1134 {
1135 // get direction of the rotational axis
1136 X = cos(alpha0_) * cos(delta0_);
1137 Y = sin(alpha0_) * cos(delta0_);
1138 Z = sin(delta0_);
1139 }
1140