Learn how easy it is to sync an existing GitHub or Google Code repo to a SourceForge project! See Demo

Close

Diff of /trunk/src/libplanet/Planet.cpp [000000] .. [r1] Maximize Restore

  Switch to side-by-side view

--- a
+++ b/trunk/src/libplanet/Planet.cpp
@@ -0,0 +1,1030 @@
+#include <cmath>
+#include <cstdio>
+#include <cctype>
+#include <string>
+using namespace std;
+
+#include "Options.h"
+#include "Planet.h"
+#include "xpUtil.h"
+
+#include "libmoons/libmoons.h"
+
+extern void
+ephemeris(const body b, const double tjd, 
+          double &X_, double &Y_, double &Z_);
+
+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, body_string[NEPTUNE], 3) == 0)
+        return_body = NEPTUNE;
+    else if (strncmp(lowercase, body_string[NEREID], 3) == 0)
+        return_body = NEREID;
+    else if (strncmp(lowercase, body_string[OBERON], 1) == 0)
+        return_body = OBERON;
+    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 Davies et al. (1996), Celestial
+  Mechanics 63, 127--148.
+*/
+Planet::Planet(const double jd, const body this_body) 
+    : index_(this_body), julianDay_(jd), needRotationMatrix_(true), period_(0)
+{
+    Options *options = Options::getInstance();
+    if (options->UniversalTime())
+	julianDay_ += delT(jd) / 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;
+	radius_ = 695000;
+	flattening_ = 0;
+
+	break;
+    case MERCURY:
+	primary_ = SUN;
+
+        alpha0_ = 281.01 - 0.033 * T2000_;
+        delta0_ = 61.45 - 0.005 * T2000_;
+
+        nullMeridian0_ = 329.68;
+        wdot_ = 6.1385025;
+
+        flipped_ = -1; 
+
+	period_ = 87.969;
+
+	radius_ = 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;
+
+	radius_ = 6051;
+	flattening_ = 0;
+
+        break;
+    case EARTH:
+	primary_ = SUN;
+	
+	alpha0_ = 0 - 0.641 * T2000_;
+	delta0_ = 90 - .557 * T2000_;
+	
+	nullMeridian0_ = 190.16;
+	wdot_ = 360.9856235;
+	
+	flipped_ = 1;
+
+	period_ = 365.256;
+	
+	radius_ = 6378;
+	flattening_ = 0.0034;
+	
+	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;
+
+	radius_ = 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.681 - 0.108 * T2000_;
+	    delta0_ = 52.886 - 0.061 * T2000_;
+	    
+	    nullMeridian0_ = 176.901;
+	    wdot_ = 350.891983;
+	    
+	    period_ = 686.980;
+	    
+	    radius_ = 3393;
+	    flattening_ = 0.0052;
+	    
+	    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;
+
+	    radius_ = 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;
+	    radius_ = 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.536;
+
+	    period_ = 4332.589;
+
+	    radius_ = 71492;
+	    flattening_ = 0.0649;
+
+            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;
+
+	    radius_ = 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;
+
+	    radius_ = 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;
+
+	    radius_ = 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;
+
+	    radius_ = 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;
+
+	    radius_ = 60268;
+	    flattening_ = 0.0980;
+
+            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;
+	    
+	    radius_ = 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;
+
+	    radius_ = 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;
+
+	    radius_ = 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;
+
+	    radius_ = 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;
+
+	    radius_ = 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;
+
+	    radius_ = 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;
+
+	    radius_ = 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;
+
+	    radius_ = 718;
+	    flattening_ = 0;
+
+	    break;
+	case PHOEBE:
+	    primary_ = SATURN;
+	    
+	    alpha0_ = 355.00;
+	    delta0_ = 68.70;
+
+	    nullMeridian0_ = 304.70;
+	    wdot_ = 930.8338720;
+
+	    period_ = 550.48;
+
+	    radius_ = 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;
+	    
+	    radius_ = 25559;
+	    flattening_ = 0.0229;
+
+	    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;
+
+	    radius_ = 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;
+
+	    radius_ = 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;
+
+	    radius_ = 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;
+
+	    radius_ = 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;
+
+	    radius_ = 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;
+	    
+	    radius_ = 24764;
+	    flattening_ = 0.017;
+
+	}
+	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;
+
+	    radius_ = 1353;
+	    flattening_ = 0;
+	}
+	break;
+	case NEREID:
+	{
+	    primary_ = NEPTUNE;
+
+	    alpha0_ = 299.36;
+	    delta0_ = 43.46;
+	    
+	    nullMeridian0_ = 0;
+	    wdot_ = 0;
+
+	    period_ = 360.13619;
+
+	    radius_ = 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;
+	    
+	    radius_ = 1150;
+	    flattening_ = 0;
+	    break;
+	case CHARON:
+	    primary_ = PLUTO;
+
+	    nullMeridian0_ = 56.77;
+
+	    period_ = 6.38723;
+
+	    radius_ = 593;
+	    flattening_ = 0;
+	    break;
+	default:
+	    break;
+	}
+    }
+    break;
+    default:
+	xpExit("Planet: Unknown body specified.\n", __FILE__, __LINE__);
+    }
+
+    alpha0_ *= deg_to_rad;
+    delta0_ *= deg_to_rad;
+
+    radius_ /= AU_to_km;
+
+    nullMeridian_ = fmod(nullMeridian0_ + wdot_ * d2000_, 360.0) * deg_to_rad;
+}
+
+// Compute heliocentric equatorial coordinates of the planet
+void
+Planet::calcHeliocentricEquatorial()
+{
+    calcHeliocentricEquatorial(true);
+}
+
+void
+Planet::calcHeliocentricEquatorial(const bool relativeToSun)
+{
+    if (primary_ == SUN)
+    {
+	ephemeris(index_, julianDay_, X_, Y_, Z_);
+    }
+    else
+    {
+	switch(primary_)
+	{
+	case EARTH:
+	    moon(julianDay_, X_, Y_, Z_);
+	    break;
+	case MARS:
+	    marsat(julianDay_, index_, X_, Y_, Z_);
+	    break;
+	case JUPITER:
+	    jupsat(julianDay_, index_, X_, Y_, Z_);
+	    break;
+	case SATURN:
+	    satsat(julianDay_, index_, X_, Y_, Z_);
+	    break;
+	case URANUS:
+	    urasat(julianDay_, index_, X_, Y_, Z_);
+	    break;
+	case NEPTUNE:
+	    nepsat(julianDay_, index_, X_, Y_, Z_);
+	    break;
+	case PLUTO:
+	    plusat(julianDay_, X_, Y_, Z_);
+	    break;
+	default:
+	    break;
+	}
+	if (relativeToSun)
+	{
+	    double Prx, Pry, Prz;
+	    ephemeris(primary_, julianDay_, Prx, Pry, Prz);
+	    X_ += Prx;
+	    Y_ += Pry;
+	    Z_ += Prz;
+	}
+    }
+}
+
+// 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;
+}
+
+//  Converts planetocentric to Heliocentric Equatorial coordinates
+void
+Planet::PlanetocentricToXYZ(double &X, double &Y, double &Z,
+			    const double lat, const double lon, 
+			    const double rad)
+{
+    if (needRotationMatrix_) CreateRotationMatrix();
+
+    double newlon = lon;
+    if (index_ == EARTH || index_ == SUN) newlon *= -1;
+    if (wdot_ > 0) newlon *= -1;
+
+    double r[3];
+    r[0] = cos(lat) * cos(newlon);
+    r[1] = cos(lat) * sin(newlon);
+    r[2] = sin(lat);
+
+    double newrad = rad * radius_;
+    X = dot(invRot_[0], r) * newrad;
+    Y = dot(invRot_[1], r) * newrad;
+    Z = dot(invRot_[2], r) * newrad;
+
+    X += X_;
+    Y += Y_;
+    Z += Z_;
+}
+
+//  Converts Heliocentric Equatorial coordinates to planetocentric
+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);
+}
+
+//  Converts Heliocentric Equatorial coordinates to planetocentric
+void
+Planet::XYZToPlanetocentric(const double X, const double Y, const double Z,
+			    double &lat, double &lon, double &rad)
+{
+    if (needRotationMatrix_) CreateRotationMatrix();
+
+    double r[3] = { X - X_, Y - Y_, Z - Z_ };
+
+    double sx = dot(rot_[0], r);
+    double sy = dot(rot_[1], r);
+    double sz = dot(rot_[2], r);
+
+    lat = atan(sz/sqrt(sx * sx + sy*sy));
+    if (cos(lat) > 1e-5)
+	lon = atan2(sy, sx);
+    else
+	lon = 0;
+
+    if (index_ == EARTH || index_ == SUN) lon *= -1;
+    if (wdot_ > 0) lon *= -1;
+    if (lon < 0) lon += TWO_PI;
+
+    rad = sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2]);
+    rad /= radius_;
+}
+
+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;
+
+    ephemeris(index_, julianDay_ - 0.5, pX0, pY0, pZ0);
+    ephemeris(index_, julianDay_ + 0.5, 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_);
+}
+