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