#include <cmath>
#include <cstdio>
#include <cctype>
#include <sstream>
#include <string>
using namespace std;
#include "findFile.h"
#include "Options.h"
#include "Planet.h"
#include "xpUtil.h"
#include "libephemeris/ephemerisWrapper.h"
body
Planet::parseBodyName(char *name)
{
body return_body;
char *lowercase = name;
char *ptr = name;
while (*ptr != '\0') *ptr++ = tolower(*name++);
if (strncmp(lowercase, "above", 2) == 0)
return_body = ABOVE_ORBIT;
else if (strncmp(lowercase, body_string[ARIEL], 2) == 0)
return_body = ARIEL;
else if (strncmp(lowercase, "below", 1) == 0)
return_body = BELOW_ORBIT;
else if (strncmp(lowercase, body_string[CALLISTO], 2) == 0)
return_body = CALLISTO;
else if (strncmp(lowercase, body_string[CHARON], 2) == 0)
return_body = CHARON;
else if (strncmp(lowercase, "default", 3) == 0)
return_body = DEFAULT;
else if (strncmp(lowercase, body_string[DEIMOS], 3) == 0)
return_body = DEIMOS;
else if (strncmp(lowercase, body_string[DIONE], 2) == 0)
return_body = DIONE;
else if (strncmp(lowercase, body_string[EARTH], 2) == 0)
return_body = EARTH;
else if (strncmp(lowercase, body_string[ENCELADUS], 2) == 0)
return_body = ENCELADUS;
else if (strncmp(lowercase, body_string[EUROPA], 2) == 0)
return_body = EUROPA;
else if (strncmp(lowercase, body_string[GANYMEDE], 1) == 0)
return_body = GANYMEDE;
else if (strncmp(lowercase, body_string[HYPERION], 1) == 0)
return_body = HYPERION;
else if (strncmp(lowercase, body_string[IAPETUS], 2) == 0)
return_body = IAPETUS;
else if (strncmp(lowercase, body_string[IO], 2) == 0)
return_body = IO;
else if (strncmp(lowercase, body_string[JUPITER], 1) == 0)
return_body = JUPITER;
else if (strncmp(lowercase, "major", 3) == 0)
return_body = MAJOR_PLANET;
else if (strncmp(lowercase, body_string[MARS], 3) == 0)
return_body = MARS;
else if (strncmp(lowercase, body_string[MERCURY], 2) == 0)
return_body = MERCURY;
else if (strncmp(lowercase, body_string[MIMAS], 3) == 0)
return_body = MIMAS;
else if (strncmp(lowercase, body_string[MIRANDA], 3) == 0)
return_body = MIRANDA;
else if (strncmp(lowercase, body_string[MOON], 2) == 0)
return_body = MOON;
else if (strncmp(lowercase, "naif", 2) == 0)
return_body = NAIF;
else if (strncmp(lowercase, body_string[NEPTUNE], 3) == 0)
return_body = NEPTUNE;
else if (strncmp(lowercase, body_string[NEREID], 3) == 0)
return_body = NEREID;
else if (strncmp(lowercase, "norad", 2) == 0)
return_body = NORAD;
else if (strncmp(lowercase, body_string[OBERON], 1) == 0)
return_body = OBERON;
else if (strncmp(lowercase, "path", 2) == 0)
return_body = ALONG_PATH;
else if (strncmp(lowercase, body_string[PHOBOS], 4) == 0)
return_body = PHOBOS;
else if (strncmp(lowercase, body_string[PHOEBE], 4) == 0)
return_body = PHOEBE;
else if (strncmp(lowercase, body_string[PLUTO], 2) == 0)
return_body = PLUTO;
else if (strncmp(lowercase, "random", 2) == 0)
return_body = RANDOM_BODY;
else if (strncmp(lowercase, body_string[RHEA], 2) == 0)
return_body = RHEA;
else if (strncmp(lowercase, body_string[SATURN], 2) == 0)
return_body = SATURN;
else if (strncmp(lowercase, body_string[SUN], 2) == 0)
return_body = SUN;
else if (strncmp(lowercase, "system", 2) == 0)
return_body = SAME_SYSTEM;
else if (strncmp(lowercase, body_string[TETHYS], 2) == 0)
return_body = TETHYS;
else if (strncmp(lowercase, body_string[TITANIA], 6) == 0)
return_body = TITANIA;
else if (strncmp(lowercase, body_string[TITAN], 5) == 0)
return_body = TITAN;
else if (strncmp(lowercase, body_string[TRITON], 2) == 0)
return_body = TRITON;
else if (strncmp(lowercase, body_string[URANUS], 2) == 0)
return_body = URANUS;
else if (strncmp(lowercase, body_string[UMBRIEL], 2) == 0)
return_body = UMBRIEL;
else if (strncmp(lowercase, body_string[VENUS], 1) == 0)
return_body = VENUS;
else
{
xpWarn("parseBody: Unknown body specified\n",
__FILE__, __LINE__);
return_body = UNKNOWN_BODY;
}
return(return_body);
}
/*
Rotational parameters are from Seidelmann et al. (2002), Celestial
Mechanics 82, 83--110.
*/
Planet::Planet(const double jd, const body this_body)
: index_(this_body),
julianDay_(jd),
needRotationMatrix_(true),
period_(0),
needShadowCoeffs_(true)
{
Options *options = Options::getInstance();
if (options->UniversalTime())
julianDay_ += delT(julianDay_) / 86400;
d2000_ = (julianDay_ - 2451545.0);
T2000_ = d2000_ / 36525;
switch (index_)
{
case SUN:
primary_ = SUN;
alpha0_ = 286.13;
delta0_ = 63.87;
nullMeridian0_ = 84.10;
wdot_ = 14.1844;
flipped_ = 1;
period_ = 0;
radiusEq_ = 696000;
flattening_ = 0;
break;
case MERCURY:
primary_ = SUN;
alpha0_ = 281.01 - 0.033 * T2000_;
delta0_ = 61.45 - 0.005 * T2000_;
nullMeridian0_ = 329.548;
wdot_ = 6.1385025;
flipped_ = -1;
period_ = 87.969;
radiusEq_ = 2439;
flattening_ = 0;
break;
case VENUS:
primary_ = SUN;
alpha0_ = 272.76;
delta0_ = 67.16;
nullMeridian0_ = 160.20;
wdot_ = -1.4813688;
flipped_ = 1;
period_ = 224.701;
radiusEq_ = 6051.9;
flattening_ = 0;
break;
case EARTH:
primary_ = SUN;
alpha0_ = 0 - 0.641 * T2000_;
delta0_ = 90 - .557 * T2000_;
nullMeridian0_ = 190.147;
wdot_ = 360.9856235;
flipped_ = 1;
period_ = 365.256;
// WGS84 ellipsoid
radiusEq_ = 6378.137;
flattening_ = 1/298.257223563;
break;
case MOON:
{
primary_ = EARTH;
const double E1 = (125.045 - 0.0529921 * d2000_) * deg_to_rad;
const double E2 = (250.089 - 0.1059842 * d2000_) * deg_to_rad;
const double E3 = (260.008 + 13.0120009 * d2000_) * deg_to_rad;
const double E4 = (176.625 + 13.3407154 * d2000_) * deg_to_rad;
const double E5 = (357.529 + 0.9856003 * d2000_) * deg_to_rad;
const double E6 = (311.589 + 26.4057084 * d2000_) * deg_to_rad;
const double E7 = (134.963 + 13.0649930 * d2000_) * deg_to_rad;
const double E8 = (276.617 + 0.3287146 * d2000_) * deg_to_rad;
const double E9 = ( 34.226 + 1.7484877 * d2000_) * deg_to_rad;
const double E10 = ( 15.134 - 0.1589763 * d2000_) * deg_to_rad;
const double E11 = (119.743 + 0.0036096 * d2000_) * deg_to_rad;
const double E12 = (239.961 + 0.1643573 * d2000_) * deg_to_rad;
const double E13 = ( 25.053 + 12.9590088 * d2000_) * deg_to_rad;
alpha0_ = 269.9949 + (0.0031 * T2000_ - 3.8787 * sin(E1)
- 0.1204 * sin(E2) + 0.0700 * sin(E3)
- 0.0172 * sin(E4) + 0.0072 * sin(E6)
- 0.0052 * sin(E10) + 0.0043 * sin(E13));
delta0_ = 66.5392 + (0.0130 * T2000_ + 1.5419 * cos(E1)
+ 0.0239 * cos(E2) - 0.0278 * cos(E3)
+ 0.0068 * cos(E4) - 0.0029 * cos(E6)
+ 0.0009 * cos(E7) + 0.0008 * cos(E10)
- 0.0009 * cos(E13));
nullMeridian0_ = 38.3213 + (d2000_ * (13.17635815 - 1.4E-12 * d2000_)
+ 3.5610 * sin(E1)
+ 0.1208 * sin(E2)
- 0.0642 * sin(E3)
+ 0.0158 * sin(E4)
+ 0.0252 * sin(E5)
- 0.0066 * sin(E6)
- 0.0047 * sin(E7)
- 0.0046 * sin(E8)
+ 0.0028 * sin(E9)
+ 0.0052 * sin(E10)
+ 0.0040 * sin(E11)
+ 0.0019 * sin(E12)
- 0.0044 * sin(E13));
wdot_ = 0;
flipped_ = 1;
period_ = 27.321661;
radiusEq_ = 1737.5;
flattening_ = 0.002;
}
break;
case MARS:
case PHOBOS:
case DEIMOS:
{
flipped_ = -1;
const double M1 = (169.51 - 0.4357640 * d2000_) * deg_to_rad;
const double M2 = (192.93 + 1128.4096700 * d2000_) * deg_to_rad;
const double M3 = (53.47 - 0.0181510 * d2000_) * deg_to_rad;
switch (index_)
{
case MARS:
primary_ = SUN;
alpha0_ = 317.68143 - 0.1061 * T2000_;
delta0_ = 52.88650 - 0.0609 * T2000_;
nullMeridian0_ = 176.630;
wdot_ = 350.89198226;
period_ = 686.980;
radiusEq_ = 3397;
flattening_ = 0.0065;
break;
case PHOBOS:
primary_ = MARS;
alpha0_ = 317.68 - 0.108 * T2000_ + 1.79 * sin(M1);
delta0_ = 52.90 - 0.061 * T2000_ - 1.08 * cos(M1);
nullMeridian0_ = (35.06 + 8.864 * T2000_ * T2000_ - 1.42 * sin(M1)
- 0.78 * sin(M2));
wdot_ = 1128.8445850;
period_ = 0.31891023;
radiusEq_ = 10; // non-spherical
flattening_ = 0;
break;
case DEIMOS:
primary_ = MARS;
alpha0_ = 316.65 - 0.108 * T2000_ + 2.98 * sin(M3);
delta0_ = 53.52 - 0.061 * T2000_ - 1.78 * cos(M3);
nullMeridian0_ = (79.41 - 0.520 * T2000_ * T2000_
- 2.58 * sin(M3) + 0.19 * cos(M3));
wdot_ = 285.1618970;
period_ = 1.2624407;
radiusEq_ = 6; // non-spherical
flattening_ = 0;
break;
default:
break;
}
}
break;
case JUPITER:
case IO:
case EUROPA:
case GANYMEDE:
case CALLISTO:
{
flipped_ = -1;
const double J3 = (283.90 + 4850.7 * T2000_) * deg_to_rad;
const double J4 = (355.80 + 1191.3 * T2000_) * deg_to_rad;
const double J5 = (119.90 + 262.1 * T2000_) * deg_to_rad;
const double J6 = (229.80 + 64.3 * T2000_) * deg_to_rad;
const double J7 = (352.25 + 2382.6 * T2000_) * deg_to_rad;
const double J8 = (113.35 + 6070.0 * T2000_) * deg_to_rad;
switch (index_)
{
case JUPITER:
primary_ = SUN;
alpha0_ = 268.05 - 0.009 * T2000_;
delta0_ = 64.49 + 0.003 * T2000_;
// System III (magnetic field rotation)
nullMeridian0_ = 284.95;
wdot_ = 870.5366420;
if (options->GRSSet())
{
// System II
nullMeridian0_ = 43.3;
wdot_ = 870.270;
}
period_ = 4332.589;
radiusEq_ = 71492;
flattening_ = 0.06487;
break;
case IO:
primary_ = JUPITER;
alpha0_ = (268.05 - 0.009 * T2000_ + 0.094 * sin(J3)
+ 0.024 * sin(J4));
delta0_ = (64.50 + 0.003 * T2000_ + 0.040 * cos(J3)
+ 0.011 * cos(J4));
nullMeridian0_ = 200.39 - 0.085 * sin(J3) - 0.022 * sin(J4);
wdot_ = 203.4889538;
period_ = 1.769137786;
radiusEq_ = 1821.6;
flattening_ = 0;
break;
case EUROPA:
primary_ = JUPITER;
alpha0_ = (268.08 - 0.009 * T2000_ + 1.086 * sin(J4)
+ 0.060 * sin(J5) + 0.015 * sin(J6)
+ 0.009 * sin(J7));
delta0_ = (64.51 + 0.003 * T2000_ + 0.468 * cos(J4)
+ 0.026 * cos(J5) + 0.007 * cos(J6)
+ 0.002 * cos(J7));
nullMeridian0_ = (35.67 - 0.980 * sin(J4) - 0.054 * sin(J5)
- 0.014 * sin(J6) - 0.008 * sin(J7));
wdot_ = 101.3747235;
period_ = 3.551181041;
radiusEq_ = 1560.8;
flattening_ = 0;
break;
case GANYMEDE:
primary_ = JUPITER;
alpha0_ = (268.20 - 0.009 * T2000_ - 0.037 * sin(J4)
+ 0.431 * sin(J5) + 0.091 * sin(J6));
delta0_ = (64.57 + 0.003 * T2000_ - 0.016 * cos(J4)
+ 0.186 * cos(J5) + 0.039 * cos(J6));
nullMeridian0_ = (44.04 + 0.033 * sin(J4) - 0.389 * sin(J5)
- 0.082 * sin(J6));
wdot_ = 50.3176081;
period_ = 7.15455296;
radiusEq_ = 2631.2;
flattening_ = 0;
break;
case CALLISTO:
primary_ = JUPITER;
alpha0_ = (268.72 - 0.009 * T2000_ - 0.068 * sin(J5)
+ 0.590 * sin(J6) + 0.010 * sin(J8));
delta0_ = (64.83 + 0.003 * T2000_ - 0.029 * cos(J5)
+ 0.254 * cos(J6) - 0.004 * cos(J8));
nullMeridian0_ = (259.73 + 0.061 * sin(J5) - 0.533 * sin(J6)
- 0.009 * sin(J8));
wdot_ = 21.5710715;
period_ = 16.6890184;
radiusEq_ = 2410.3;
flattening_ = 0;
break;
default:
break;
}
}
break;
case SATURN:
case MIMAS:
case ENCELADUS:
case TETHYS:
case DIONE:
case RHEA:
case TITAN:
case HYPERION:
case IAPETUS:
case PHOEBE:
{
flipped_ = -1;
const double S3 = (177.40 - 36505.5 * T2000_) * deg_to_rad;
const double S4 = (300.00 - 7225.9 * T2000_) * deg_to_rad;
const double S5 = (316.45 + 506.2 * T2000_) * deg_to_rad;
const double S6 = (345.20 - 1016.3 * T2000_) * deg_to_rad;
const double S7 = ( 29.80 - 52.1 * T2000_) * deg_to_rad;
switch (index_)
{
case SATURN:
primary_ = SUN;
alpha0_ = 40.589 - 0.036 * T2000_;
delta0_ = 83.537 - 0.004 * T2000_;
nullMeridian0_ = 38.90;
wdot_ = 810.7939024;
period_ = 10759.22;
radiusEq_ = 60268;
flattening_ = 0.09796;
break;
case MIMAS:
primary_ = SATURN;
alpha0_ = 40.66 - 0.036 * T2000_ + 13.56 * sin(S3);
delta0_ = 83.52 - 0.004 * T2000_ - 1.53 * cos(S3);
nullMeridian0_ = 337.46 - 13.48 * sin(S3) - 44.85 * sin(S5);
wdot_ = 381.994550;
period_ = 0.942421813;
radiusEq_ = 198.6;
flattening_ = 0;
break;
case ENCELADUS:
primary_ = SATURN;
alpha0_ = 40.66 - 0.036 * T2000_;
delta0_ = 83.52 - 0.004 * T2000_;
nullMeridian0_ = 2.82;
wdot_ = 262.7318996;
period_ = 1.370217855;
radiusEq_ = 249.4;
flattening_ = 0;
break;
case TETHYS:
primary_ = SATURN;
alpha0_ = 40.66 - 0.036 * T2000_ + 9.66 * sin(S4);
delta0_ = 83.52 - 0.004 * T2000_ - 1.09 * cos(S4);
nullMeridian0_ = 10.45 - 9.60 * sin(S4) + 2.23 * sin(S5);
wdot_ = 190.6979085;
period_ = 1.887802160;
radiusEq_ = 529.8;
flattening_ = 0;
break;
case DIONE:
primary_ = SATURN;
alpha0_ = 40.66 - 0.036 * T2000_;
delta0_ = 83.52 - 0.004 * T2000_;
nullMeridian0_ = 357.00;
wdot_ = 131.5349316;
period_ = 2.736914742;
radiusEq_ = 559;
flattening_ = 0;
break;
case RHEA:
primary_ = SATURN;
alpha0_ = 40.38 - 0.036 * T2000_ + 3.10 * sin(S6);
delta0_ = 83.55 - 0.004 * T2000_ - 0.35 * cos(S6);
nullMeridian0_ = 235.16 - 3.08 * sin(S6);
wdot_ = 79.6900478;
period_ = 4.517500436;
radiusEq_ = 764;
flattening_ = 0;
break;
case TITAN:
primary_ = SATURN;
alpha0_ = 36.41 - 0.036 * T2000_ + 2.66 * sin(S7);
delta0_ = 83.94 - 0.004 * T2000_ - 0.30 * cos(S7);
nullMeridian0_ = 189.64 - 2.64 * sin(S7);
wdot_ = 22.5769768;
period_ = 15.94542068;
radiusEq_ = 2575;
flattening_ = 0;
break;
case HYPERION:
primary_ = SATURN;
alpha0_ = 40.589 - 0.036 * T2000_;
delta0_ = 83.537 - 0.004 * T2000_;
nullMeridian0_ = 0;
wdot_ = 0;
period_ = 21.2766088;
radiusEq_ = 141.5; // non-spherical
flattening_ = 0;
break;
case IAPETUS:
primary_ = SATURN;
alpha0_ = 318.16 - 3.949 * T2000_;
delta0_ = 75.03 - 1.143 * T2000_;
nullMeridian0_ = 350.20;
wdot_ = 4.5379572;
period_ = 79.3301825;
radiusEq_ = 718;
flattening_ = 0;
break;
case PHOEBE:
primary_ = SATURN;
alpha0_ = 355.00;
delta0_ = 68.70;
nullMeridian0_ = 304.70;
wdot_ = 930.8338720;
period_ = 550.48;
radiusEq_ = 110;
flattening_ = 0;
break;
default:
break;
}
}
break;
case URANUS:
case MIRANDA:
case ARIEL:
case UMBRIEL:
case TITANIA:
case OBERON:
{
flipped_ = 1;
const double U11 = (141.69 + 41887.66 * T2000_) * deg_to_rad;
const double U12 = (316.41 + 2863.96 * T2000_) * deg_to_rad;
const double U13 = (304.01 - 51.94 * T2000_) * deg_to_rad;
const double U14 = (308.71 - 93.17 * T2000_) * deg_to_rad;
const double U15 = (340.82 - 75.32 * T2000_) * deg_to_rad;
const double U16 = (259.14 - 504.81 * T2000_) * deg_to_rad;
switch (index_)
{
case URANUS:
primary_ = SUN;
alpha0_ = 257.311;
delta0_ = -15.175;
nullMeridian0_ = 203.81;
wdot_ = -501.1600928;
period_ = 30685.4;
radiusEq_ = 25559;
flattening_ = 0.02293;
break;
case MIRANDA:
primary_ = URANUS;
alpha0_ = 257.42 + 4.41 * sin(U11) - 0.04 * sin(2*U11);
delta0_ = -15.08 + 4.25 * cos(U11) - 0.02 * cos(2*U11);
nullMeridian0_ = (30.70 - 1.27 * sin(U12) + 0.15 * sin(2*U12)
+ 1.15 * sin(U11) - 0.09 * sin(2*U11));
wdot_ = -254.6906892;
period_ = 1.41347925;
radiusEq_ = 235.8;
flattening_ = 0;
break;
case ARIEL:
primary_ = URANUS;
alpha0_ = 257.43 + 0.29 * sin(U13);
delta0_ = -15.10 + 0.28 * cos(U13);
nullMeridian0_ = 156.22 + 0.05 * sin(U12) + 0.08 * sin(U13);
wdot_ = -142.8356681;
period_ = 2.52037935;
radiusEq_ = 578.9;
flattening_ = 0;
break;
case UMBRIEL:
primary_ = URANUS;
alpha0_ = 257.43 + 0.21 * sin(U14);
delta0_ = -15.10 + 0.20 * cos(U14);
nullMeridian0_ = 108.05 - 0.09 * sin(U12) + 0.06 * sin(U14);
wdot_ = -86.8688923;
period_ = 4.1441772;
radiusEq_ = 584.7;
flattening_ = 0;
break;
case TITANIA:
primary_ = URANUS;
alpha0_ = 257.43 + 0.29 * sin(U15);
delta0_ = -15.10 + 0.28 * cos(U15);
nullMeridian0_ = 77.74 + 0.08 * sin(U15);
wdot_ = -41.3514316;
period_ = 8.7058717;
radiusEq_ = 788.9;
flattening_ = 0;
break;
case OBERON:
primary_ = URANUS;
alpha0_ = 257.43 + 0.16 * sin(U16);
delta0_ = -15.10 + 0.16 * cos(U16);
nullMeridian0_ = 6.77 + 0.04 * sin(U16);
wdot_ = -26.7394932;
period_ = 13.4632389;
radiusEq_ = 761.4;
flattening_ = 0;
break;
default:
break;
}
}
break;
case NEPTUNE:
case TRITON:
case NEREID:
{
flipped_ = -1;
switch (index_)
{
case NEPTUNE:
{
primary_ = SUN;
double N = (357.85 + 52.316 * T2000_) * deg_to_rad;
alpha0_ = 299.36 + 0.70 * sin(N);
delta0_ = 43.46 - 0.51 * cos(N);
nullMeridian0_ = 253.18 - 0.48 * sin(N);
wdot_ = 536.3128492;
period_ = 60189.0;
radiusEq_ = 24764;
flattening_ = 0.0171;
}
break;
case TRITON:
{
primary_ = NEPTUNE;
const double N7 = (177.85 + 52.316 * T2000_) * deg_to_rad;
alpha0_ = (299.36 - 32.35 * sin(N7) - 6.28 * sin(2*N7)
- 2.08 * sin(3*N7) - 0.74 * sin(4*N7)
- 0.28 * sin(5*N7) - 0.11 * sin(6*N7)
- 0.07 * sin(7*N7) - 0.02 * sin(8*N7)
- 0.01 * sin(9*N7));
delta0_ = (43.46 + 22.55 * cos(N7) + 2.10 * cos(2*N7)
+ 0.55 * cos(3*N7) + 0.16 * cos(4*N7)
+ 0.05 * cos(5*N7) + 0.02 * cos(6*N7)
+ 0.01 * cos(7*N7));
nullMeridian0_ = (296.53 + 22.25 * sin(N7) + 6.73 * sin(2*N7)
+ 2.05 * sin(3*N7) + 0.74 * sin(4*N7)
+ 0.28 * sin(5*N7) + 0.11 * sin(6*N7)
+ 0.05 * sin(7*N7) + 0.02 * sin(8*N7)
+ 0.01 * sin(9*N7));
wdot_ = -61.2572637;
period_ = 5.8768541;
radiusEq_ = 1353;
flattening_ = 0;
}
break;
case NEREID:
{
primary_ = NEPTUNE;
alpha0_ = 299.36;
delta0_ = 43.46;
nullMeridian0_ = 0;
wdot_ = 0;
period_ = 360.13619;
radiusEq_ = 170;
flattening_ = 0;
}
break;
default:
break;
}
}
break;
case PLUTO:
case CHARON:
{
flipped_ = 1;
alpha0_ = 313.02;
delta0_ = 9.09;
wdot_ = -56.3623195;
switch (index_)
{
case PLUTO:
primary_ = SUN;
nullMeridian0_ = 236.77;
period_ = 90465.0;
radiusEq_ = 1151;
flattening_ = 0;
break;
case CHARON:
primary_ = PLUTO;
nullMeridian0_ = 56.77;
period_ = 6.38723;
radiusEq_ = 593;
flattening_ = 0;
break;
default:
break;
}
}
break;
default:
xpExit("Planet: Unknown body specified.\n", __FILE__, __LINE__);
}
alpha0_ *= deg_to_rad;
delta0_ *= deg_to_rad;
radiusEq_ /= AU_to_km;
nullMeridian_ = fmod(nullMeridian0_ + wdot_ * d2000_, 360.0);
nullMeridian_ *= deg_to_rad;
radiusPol_ = radiusEq_ * (1 - flattening_);
omf2_ = (1 - flattening_) * (1 - flattening_);
}
Planet::~Planet()
{
}
// Return the distance from the surface to the center, in units of
// equatorial radius. Input latitude is assumed to be planetographic.
double
Planet::Radius(const double lat) const
{
double returnValue = 1;
if (index_ == JUPITER ||
index_ == SATURN)
{
double tmpLat = lat;
double tmpLon = 0;
PlanetographicToPlanetocentric(tmpLat, tmpLon);
const double ReSinLat = radiusEq_ * sin(tmpLat);
const double RpCosLat = radiusPol_ * cos(tmpLat);
// from http://mathworld.wolfram.com/Ellipse.html
returnValue = radiusPol_ / sqrt(RpCosLat * RpCosLat
+ ReSinLat * ReSinLat);
}
return(returnValue);
}
// Compute heliocentric equatorial coordinates of the planet
void
Planet::calcHeliocentricEquatorial()
{
calcHeliocentricEquatorial(true);
}
void
Planet::calcHeliocentricEquatorial(const bool relativeToSun)
{
GetHeliocentricXYZ(index_, primary_, julianDay_,
relativeToSun, X_, Y_, Z_);
}
// Return rectangular coordinates in heliocentric equatorial frame
void
Planet::getPosition(double &X, double &Y, double &Z) const
{
X = X_;
Y = Y_;
Z = Z_;
}
void
Planet::CreateRotationMatrix()
{
/*
These equations are from "Astronomy on the Personal Computer,
3rd Edition", by Montenbruck and Pfleger, Chapter 7.
*/
const double cosW = cos(nullMeridian_);
const double sinW = sin(nullMeridian_);
const double cosA = cos(alpha0_);
const double sinA = sin(alpha0_);
const double cosD = cos(delta0_);
const double sinD = sin(delta0_);
rot_[0][0] = -cosW * sinA - sinW * sinD * cosA;
rot_[0][1] = cosW * cosA - sinW * sinD * sinA;
rot_[0][2] = sinW * cosD;
rot_[1][0] = sinW * sinA - cosW * sinD * cosA;
rot_[1][1] = -sinW * cosA - cosW * sinD * sinA;
rot_[1][2] = cosW * cosD;
rot_[2][0] = cosD * cosA;
rot_[2][1] = cosD * sinA;
rot_[2][2] = sinD;
invertMatrix(rot_, invRot_);
needRotationMatrix_ = false;
}
void
Planet::PlanetocentricToXYZ(double &X, double &Y, double &Z,
const double lat, const double lon,
const double rad)
{
if (needRotationMatrix_) CreateRotationMatrix();
double r[3];
r[0] = cos(lat) * cos(lon);
r[1] = cos(lat) * sin(lon);
r[2] = sin(lat);
double newrad = rad * radiusEq_;
X = dot(invRot_[0], r) * newrad;
Y = dot(invRot_[1], r) * newrad;
Z = dot(invRot_[2], r) * newrad;
X += X_;
Y += Y_;
Z += Z_;
}
void
Planet::PlanetographicToXYZ(double &X, double &Y, double &Z,
double lat, double lon, const double rad)
{
PlanetographicToPlanetocentric(lat, lon);
PlanetocentricToXYZ(X, Y, Z, lat, lon, rad);
}
void
Planet::XYZToPlanetocentric(const double X, const double Y, const double Z,
double &lat, double &lon)
{
double rad = 0;
XYZToPlanetocentric(X, Y, Z, lat, lon, rad);
}
void
Planet::XYZToPlanetocentric(const double X, const double Y, const double Z,
double &lat, double &lon, double &rad)
{
if (needRotationMatrix_) CreateRotationMatrix();
const double r[3] = { X - X_, Y - Y_, Z - Z_ };
const double sx = dot(rot_[0], r);
const double sy = dot(rot_[1], r);
const double sz = dot(rot_[2], r);
rad = sqrt(sx * sx + sy*sy);
if (rad > 0)
lat = atan(sz/rad);
else
lat = 0;
if (cos(lat) > 1e-5)
lon = atan2(sy, sx);
else
lon = 0;
if (lon < 0) lon += TWO_PI;
rad = sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2]);
rad /= radiusEq_;
}
void
Planet::XYZToPlanetographic(const double X, const double Y, const double Z,
double &lat, double &lon)
{
double rad = 0;
XYZToPlanetographic(X, Y, Z, lat, lon, rad);
}
void
Planet::XYZToPlanetographic(const double X, const double Y, const double Z,
double &lat, double &lon, double &rad)
{
XYZToPlanetocentric(X, Y, Z, lat, lon, rad);
PlanetocentricToPlanetographic(lat, lon);
}
void
Planet::XYZToPlanetaryXYZ(const double X, const double Y, const double Z,
double &pX, double &pY, double &pZ)
{
double lat, lon, rad;
XYZToPlanetocentric(X, Y, Z, lat, lon, rad);
pX = rad * cos(lat) * cos(lon);
pY = rad * cos(lat) * sin(lon);
pZ = rad * sin(lat);
}
void
Planet::PlanetaryXYZToXYZ(const double pX, const double pY, const double pZ,
double &X, double &Y, double &Z)
{
const double lat = atan(pZ / sqrt(pX * pX + pY * pY));
double lon;
if (cos(lat) > 1e-5)
lon = atan2(pY, pX);
else
lon = 0;
const double rad = sqrt(pX * pX + pY * pY + pZ * pZ);
PlanetocentricToXYZ(X, Y, Z, lat, lon, rad);
}
void
Planet::PlanetocentricToPlanetographic(double &lat, double &lon) const
{
if (flattening_ > 0)
lat = atan(tan(lat) / omf2_);
if (index_ == EARTH || index_ == SUN) lon *= -1;
if (wdot_ > 0) lon *= -1;
if (lon < 0) lon += TWO_PI;
}
void
Planet::PlanetographicToPlanetocentric(double &lat, double &lon) const
{
if (flattening_ > 0)
lat = atan(omf2_ * tan(lat));
if (index_ == EARTH || index_ == SUN) lon *= -1;
if (wdot_ > 0) lon *= -1;
}
void
Planet::ComputeShadowCoeffs()
{
XYZToPlanetaryXYZ(0, 0, 0, sunX_, sunY_, sunZ_);
ellipseCoeffC_ = sunZ_ * sunZ_ / omf2_;
ellipseCoeffC_ += sunY_ * sunY_;
ellipseCoeffC_ += sunX_ * sunX_;
ellipseCoeffC_ -= 1;
needShadowCoeffs_ = false;
}
// x, y, z must be in planetary XYZ frame
bool
Planet::IsInMyShadow(const double x, const double y, const double z)
{
if (needShadowCoeffs_) ComputeShadowCoeffs();
double ellipseCoeffA = (z - sunZ_) * (z - sunZ_) / omf2_;
ellipseCoeffA += (y - sunY_) * (y - sunY_);
ellipseCoeffA += (x - sunX_) * (x - sunX_);
double ellipseCoeffB = (z - sunZ_) * sunZ_ / omf2_;
ellipseCoeffB += (y - sunY_) * sunY_;
ellipseCoeffB += (x - sunX_) * sunX_;
const double determinant = (ellipseCoeffB * ellipseCoeffB
- ellipseCoeffA * ellipseCoeffC_);
return(determinant > 0);
}
void
Planet::getOrbitalNorth(double &X, double &Y, double &Z) const
{
if (index_ == SUN)
{
X = cos(delta0_) * cos(alpha0_);
Y = cos(delta0_) * sin(alpha0_);
Z = sin(delta0_);
return;
}
// cross product of position and velocity vectors points to the
// orbital north pole
double pos[3] = { X_, Y_, Z_ };
double pX0, pY0, pZ0;
double pX1, pY1, pZ1;
GetHeliocentricXYZ(index_, primary_, julianDay_ - 0.5,
true, pX0, pY0, pZ0);
GetHeliocentricXYZ(index_, primary_, julianDay_ + 0.5,
true, pX1, pY1, pZ1);
double vel[3] = { pX1 - pX0, pY1 - pY0, pZ1 - pZ0 };
double north[3];
cross(pos, vel, north);
double mag = sqrt(dot(north, north));
X = north[0]/mag;
Y = north[1]/mag;
Z = north[2]/mag;
}
void
Planet::getGalacticNorth(double &X, double &Y, double &Z) const
{
// B1950.0 coordinates of the galactic north pole
const double GN_LAT = 27.4 * deg_to_rad;
const double GN_LON = 192.25 * deg_to_rad;
X = cos(GN_LAT) * cos(GN_LON);
Y = cos(GN_LAT) * sin(GN_LON);
Z = sin(GN_LAT);
}
void
Planet::getBodyNorth(double &X, double &Y, double &Z) const
{
// get direction of the rotational axis
X = cos(alpha0_) * cos(delta0_);
Y = sin(alpha0_) * cos(delta0_);
Z = sin(delta0_);
}
syntax highlighted by Code2HTML, v. 0.9.1