customorbit.cpp

Go to the documentation of this file.

// customorbit.cpp
//
// Copyright (C) 2001, Chris Laurel <claurel@shatters.net>
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 2
// of the License, or (at your option) any later version.

#include <cassert>
#include <vector>
#include <fstream>
#include <iomanip>
#include <celmath/mathlib.h>
#include "astro.h"
#include "customorbit.h"
#include "vsop87.h"
#include "jpleph.h"

using namespace std;


#define TWOPI 6.28318530717958647692
#define LPEJ 0.23509484  // Longitude of perihelion of Jupiter

// These are required because the orbits of the Jovian and Saturnian
// satellites are computed in units of their parent planets' radii.
static const double JupiterRadius = 71398.0;
static const double SaturnRadius = 60330.0;

// The expressions for custom orbits are complex, so the bounding radii are
// generally must computed from mean orbital elements.  It's important that
// a sphere with the bounding radius completely enclose the orbit, so we
// multiply by this factor to make the bounding radius a bit larger than
// the apocenter distance computed from the mean elements.
static const double BoundingRadiusSlack = 1.2;

static bool jplephInitialized = false;
static JPLEphemeris* jpleph = NULL;


double gPlanetElements[8][9];
double gElements[8][23] = {
        {   /*     mercury... */

            178.179078, 415.2057519,    3.011e-4,       0.0,
            75.899697,  1.5554889,      2.947e-4,       0.0,
            .20561421,  2.046e-5,       3e-8,           0.0,
            7.002881,   1.8608e-3,      -1.83e-5,       0.0,
            47.145944,  1.1852083,      1.739e-4,       0.0,
            .3870986,   6.74,           -0.42
        },

        {   /*     venus... */

            342.767053, 162.5533664,    3.097e-4,       0.0,
            130.163833, 1.4080361,      -9.764e-4,      0.0,
            6.82069e-3, -4.774e-5,      9.1e-8,         0.0,
            3.393631,   1.0058e-3,      -1e-6,          0.0,
            75.779647,  .89985,         4.1e-4,         0.0,
            .7233316,   16.92,          -4.4
        },

        {   /*     mars... */

            293.737334, 53.17137642,    3.107e-4,       0.0,
            3.34218203e2, 1.8407584,    1.299e-4,       -1.19e-6,
            9.33129e-2, 9.2064e-5,      7.7e-8,         0.0,
            1.850333,   -6.75e-4,       1.26e-5,        0.0,
            48.786442,  .7709917,       -1.4e-6,        -5.33e-6,
            1.5236883,  9.36,           -1.52
        },

        {   /*     jupiter... */

            238.049257, 8.434172183,    3.347e-4,       -1.65e-6,
            1.2720972e1, 1.6099617,     1.05627e-3,     -3.43e-6,
            4.833475e-2, 1.6418e-4,     -4.676e-7,      -1.7e-9,
            1.308736,   -5.6961e-3,     3.9e-6,         0.0,
            99.443414,  1.01053,        3.5222e-4,      -8.51e-6,
            5.202561,   196.74,         -9.4
        },

        {   /*     saturn... */

            266.564377, 3.398638567,    3.245e-4,       -5.8e-6,
            9.1098214e1, 1.9584158,     8.2636e-4,      4.61e-6,
            5.589232e-2, -3.455e-4,     -7.28e-7,       7.4e-10,
            2.492519,   -3.9189e-3,     -1.549e-5,      4e-8,
            112.790414, .8731951,       -1.5218e-4,     -5.31e-6,
            9.554747,   165.6,          -8.88
        },

        {   /*     uranus... */

            244.19747,  1.194065406,    3.16e-4,        -6e-7,
            1.71548692e2, 1.4844328,    2.372e-4,       -6.1e-7,
            4.63444e-2, -2.658e-5,      7.7e-8,         0.0,
            .772464,    6.253e-4,       3.95e-5,        0.0,
            73.477111,  .4986678,       1.3117e-3,      0.0,
            19.21814,   65.8,           -7.19
        },

        {   /*     neptune... */

            84.457994,  .6107942056,    3.205e-4,       -6e-7,
            4.6727364e1, 1.4245744,     3.9082e-4,      -6.05e-7,
            8.99704e-3, 6.33e-6,        -2e-9,          0.0,
            1.779242,   -9.5436e-3,     -9.1e-6,        0.0,
            130.681389, 1.098935,       2.4987e-4,      -4.718e-6,
            30.10957,   62.2,           -6.87
        },

        {   /*     pluto...(osculating 1984 jan 21) */

            95.3113544, .3980332167,    0.0,            0.0,
            224.017,    0.0,            0.0,            0.0,
            .25515,     0.0,            0.0,            0.0,
            17.1329,    0.0,            0.0,            0.0,
            110.191,    0.0,            0.0,            0.0,
            39.8151,    8.2,            -1.0
        }
};


// Useful version of trig functions which operate on values in degrees instead
// of radians.
static double sinD(double theta)
{
    return sin(degToRad(theta));
}

static double cosD(double theta)
{
    return cos(degToRad(theta));
}


static double Obliquity(double t)
{
    // Parameter t represents the Julian centuries elapsed since 1900.
    // In other words, t = (jd - 2415020.0) / 36525.0

    return degToRad(2.345229444E1 - ((((-1.81E-3*t)+5.9E-3)*t+4.6845E1)*t)/3600.0);
}

static void Nutation(double t, double &deps, double& dpsi)
{
    // Parameter t represents the Julian centuries elapsed since 1900.
    // In other words, t = (jd - 2415020.0) / 36525.0

    double ls, ld;      // sun's mean longitude, moon's mean longitude
    double ms, md;      // sun's mean anomaly, moon's mean anomaly
    double nm;      // longitude of moon's ascending node
    double t2;
    double tls, tnm, tld;       // twice above
    double a, b;

    t2 = t*t;

    a = 100.0021358*t;
    b = 360.*(a-(int)a);
    ls = 279.697+.000303*t2+b;

    a = 1336.855231*t;
    b = 360.*(a-(int)a);
    ld = 270.434-.001133*t2+b;

    a = 99.99736056000026*t;
    b = 360.*(a-(int)a);
    ms = 358.476-.00015*t2+b;

    a = 13255523.59*t;
    b = 360.*(a-(int)a);
    md = 296.105+.009192*t2+b;

    a = 5.372616667*t;
    b = 360.*(a-(int)a);
    nm = 259.183+.002078*t2-b;

    //convert to radian forms for use with trig functions.
    tls = 2*degToRad(ls);
    nm = degToRad(nm);
    tnm = 2*degToRad(nm);
    ms = degToRad(ms);
    tld = 2*degToRad(ld);
    md = degToRad(md);

    // find delta psi and eps, in arcseconds.
    dpsi = (-17.2327-.01737*t)*sin(nm)+(-1.2729-.00013*t)*sin(tls)
        +.2088*sin(tnm)-.2037*sin(tld)+(.1261-.00031*t)*sin(ms)
        +.0675*sin(md)-(.0497-.00012*t)*sin(tls+ms)
        -.0342*sin(tld-nm)-.0261*sin(tld+md)+.0214*sin(tls-ms)
        -.0149*sin(tls-tld+md)+.0124*sin(tls-nm)+.0114*sin(tld-md);
    deps = (9.21+.00091*t)*cos(nm)+(.5522-.00029*t)*cos(tls)
        -.0904*cos(tnm)+.0884*cos(tld)+.0216*cos(tls+ms)
        +.0183*cos(tld-nm)+.0113*cos(tld+md)-.0093*cos(tls-ms)
        -.0066*cos(tls-nm);

    // convert to radians.
    dpsi = degToRad(dpsi/3600);
    deps = degToRad(deps/3600);
}

static void EclipticToEquatorial(double t, double fEclLat, double fEclLon,
                                 double& RA, double& dec)
{
    // Parameter t represents the Julian centuries elapsed since 1900.
    // In other words, t = (jd - 2415020.0) / 36525.0

    double seps, ceps;  // sin and cos of mean obliquity
    double sx, cx, sy, cy, ty;
    double eps;
    double deps, dpsi;

    t = (astro::J2000 - 2415020.0) / 36525.0;
    t = 0;
    eps = Obliquity(t);         // mean obliquity for date
    Nutation(t, deps, dpsi);
    eps += deps;
    seps = sin(eps);
    ceps = cos(eps);

    sy = sin(fEclLat);
    cy = cos(fEclLat);                          // always non-negative
    if (fabs(cy)<1e-20)
        cy = 1e-20;             // insure > 0
    ty = sy/cy;
    cx = cos(fEclLon);
    sx = sin(fEclLon);
    dec = asin((sy*ceps)+(cy*seps*sx));
    RA = atan(((sx*ceps)-(ty*seps))/cx);
    if (cx<0)
        RA += PI;               // account for atan quad ambiguity
    RA = pfmod(RA, TWOPI);
}


// Convert equatorial coordinates from one epoch to another.  Method is from
// Chapter 21 of Meeus's _Astronomical Algorithms_
void EpochConvert(double jdFrom, double jdTo,
                  double a0, double d0,
                  double& a, double& d)
{
    double T = (jdFrom - astro::J2000) / 36525.0;
    double t = (jdTo - jdFrom) / 36525.0;
            
    double zeta = (2306.2181 + 1.39656 * T - 0.000139 * T * T) * t +
        (0.30188 - 0.000344 * T) * t * t + 0.017998 * t * t * t;
    double z = (2306.2181 + 1.39656 * T - 0.000139 * T * T) * t +
        (1.09468 + 0.000066 * T) * t * t + 0.018203 * t * t * t;
    double theta = (2004.3109 - 0.85330 * T - 0.000217 * T * T) * t -
        (0.42665 + 0.000217 * T) * t * t - 0.041833 * t * t * t;
    zeta  = degToRad(zeta / 3600.0);
    z     = degToRad(z / 3600.0);
    theta = degToRad(theta / 3600.0);

    double A = cos(d0) * sin(a0 + zeta);
    double B = cos(theta) * cos(d0) * cos(a0 + zeta) -
        sin(theta) * sin(d0);
    double C = sin(theta) * cos(d0) * cos(a0 + zeta) +
        cos(theta) * sin(d0);

    a = atan2(A, B) + z;
    d = asin(C);
}


double meanAnomalySun(double t)
{
    double t2, a, b;

        t2 = t*t;
        a = 9.999736042e1*t;
        b = 360*(a - (int)a);
        
    return degToRad(3.5847583e2 - (1.5e-4 + 3.3e-6*t)*t2 + b);
}

void auxJSun(double t, double* x1, double* x2, double* x3, double* x4,
             double* x5, double* x6)
{
    *x1 = t/5+0.1;
    *x2 = pfmod(4.14473+5.29691e1*t, TWOPI);
    *x3 = pfmod(4.641118+2.132991e1*t, TWOPI);
    *x4 = pfmod(4.250177+7.478172*t, TWOPI);
    *x5 = 5 * *x3 - 2 * *x2;
    *x6 = 2 * *x2 - 6 * *x3 + 3 * *x4;
}

void computePlanetElements(double t, vector<int> pList)
{
    // Parameter t represents the Julian centuries elapsed since 1900.
    // In other words, t = (jd - 2415020.0) / 36525.0

    double *ep, *pp;
    double aa;
    int planet;

    for(unsigned i = 0; i < pList.size(); i++)
    {
        planet = pList[i];
        ep = gElements[planet];
        pp = gPlanetElements[planet];
        aa = ep[1]*t;
        pp[0] = ep[0] + 360*(aa-(int)aa) + (ep[3]*t + ep[2])*t*t;
        *pp = pfmod(*pp, 360.0);
        pp[1] = (ep[1]*9.856263e-3) + (ep[2] + ep[3])/36525;

        for(unsigned j = 4; j < 20; j += 4)
            pp[j/4+1] = ((ep[j+3]*t + ep[j+2])*t + ep[j+1])*t + ep[j+0];

        pp[6] = ep[20];
        pp[7] = ep[21];
        pp[8] = ep[22];
    }
}

void computePlanetCoords(int p, double map, double da, double dhl, double dl,
                         double dm, double dml, double dr, double ds,
                         double& eclLong, double& eclLat, double& distance)
{
    double s, ma, nu, ea, lp, om, lo, slo, clo, inc, spsi, y;

    s = gPlanetElements[p][3] + ds;
    ma = map + dm;
    astro::anomaly(ma, s, nu, ea);
    distance = (gPlanetElements[p][6] + da)*(1 - s*s)/(1 + s*cos(nu));
    lp = radToDeg(nu) + gPlanetElements[p][2] + radToDeg(dml - dm);
    lp = degToRad(lp);
    om = degToRad(gPlanetElements[p][5]);
    lo = lp - om;
    slo = sin(lo);
    clo = cos(lo);
    inc = degToRad(gPlanetElements[p][4]);
    distance += dr;
    spsi = slo*sin(inc);
    y = slo*cos(inc);
    eclLat = asin(spsi) + dhl;
    spsi = sin(eclLat);
    eclLong = atan(y/clo) + om + degToRad(dl);
    if (clo < 0)
        eclLong += PI;
    eclLong = pfmod(eclLong, TWOPI);
    distance *= KM_PER_AU;
}

void ComputeGalileanElements(double t,
                             double& l1, double& l2, double& l3, double& l4,
                             double& p1, double& p2, double& p3, double& p4,
                             double& w1, double& w2, double& w3, double& w4,
                             double& gamma, double& phi, double& psi,
                             double& G, double& Gp)
{
    // Parameter t is Julian days, epoch 1950.0.
    l1 = 1.8513962 + 3.551552269981*t;
    l2 = 3.0670952 + 1.769322724929*t;
    l3 = 2.1041485 + 0.87820795239*t;
    l4 = 1.473836 + 0.37648621522*t;

    p1 = 1.69451 + 2.8167146e-3*t;
    p2 = 2.702927 + 8.248962e-4*t;
    p3 = 3.28443 + 1.24396e-4*t;
    p4 = 5.851859 + 3.21e-5*t;

    w1 = 5.451267 - 2.3176901e-3*t;
    w2 = 1.753028 - 5.695121e-4*t;
    w3 = 2.080331 - 1.25263e-4*t;
    w4 = 5.630757 - 3.07063e-5*t;

    gamma = 5.7653e-3*sin(2.85674 + 1.8347e-5*t) + 6.002e-4*sin(0.60189 - 2.82274e-4*t);
    phi = 3.485014 + 3.033241e-3*t;
    psi = 5.524285 - 3.63e-8*t;
    G = 0.527745 + 1.45023893e-3*t + gamma;
    Gp = 0.5581306 + 5.83982523e-4*t;
}




class MercuryOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
    const int p = 0;  //Planet 0
    vector<int> pList;
    double t;
    double map[4];
    double dl, dr, dml, ds, dm, da, dhl;
    double eclLong, eclLat, distance;    //heliocentric longitude, latitude, distance

    dl = dr = dml = ds = dm = da = dhl = 0.0;

    // Calculate the Julian centuries elapsed since 1900
    t = (jd - 2415020.0)/36525.0;

    // Specify which planets we must compute elements for
    pList.push_back(0);
    pList.push_back(1);
    pList.push_back(3);
    computePlanetElements(t, pList);

    // Compute necessary planet mean anomalies
    map[0] = degToRad(gPlanetElements[0][0] - gPlanetElements[0][2]);
    map[1] = degToRad(gPlanetElements[1][0] - gPlanetElements[1][2]);
    map[2] = 0.0;
    map[3] = degToRad(gPlanetElements[3][0] - gPlanetElements[3][2]);

    // Compute perturbations
    dl = 2.04e-3*cos(5*map[1]-2*map[0]+2.1328e-1)+
         1.03e-3*cos(2*map[1]-map[0]-2.8046)+
         9.1e-4*cos(2*map[3]-map[0]-6.4582e-1)+
         7.8e-4*cos(5*map[1]-3*map[0]+1.7692e-1);

    dr = 7.525e-6*cos(2*map[3]-map[0]+9.25251e-1)+
         6.802e-6*cos(5*map[1]-3*map[0]-4.53642)+
         5.457e-6*cos(2*map[1]-2*map[0]-1.24246)+
         3.569e-6*cos(5*map[1]-map[0]-1.35699);

    computePlanetCoords(p, map[p], da, dhl, dl, dm, dml, dr, ds,
                        eclLong, eclLat, distance);

    // Corrections for internal coordinate system
    eclLat -= (PI/2);
    eclLong += PI;

    return Point3d(cos(eclLong) * sin(eclLat) * distance,
                   cos(eclLat) * distance,
                   -sin(eclLong) * sin(eclLat) * distance);
    };

    double getPeriod() const
    {
        return 87.9522;
    };

    double getBoundingRadius() const
    {
        return 6.98e+7 * BoundingRadiusSlack;
    };
};

class VenusOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
    const int p = 1;  //Planet 1
    vector<int> pList;
    double t;
    double map[4], mas;
    double dl, dr, dml, ds, dm, da, dhl;
    double eclLong, eclLat, distance;    //heliocentric longitude, latitude, distance

    dl = dr = dml = ds = dm = da = dhl = 0.0;

    //Calculate the Julian centuries elapsed since 1900
        t = (jd - 2415020.0)/36525.0;

    mas = meanAnomalySun(t);

    //Specify which planets we must compute elements for
    pList.push_back(1);
    pList.push_back(3);
    computePlanetElements(t, pList);

    //Compute necessary planet mean anomalies
    map[0] = 0.0;
    map[1] = degToRad(gPlanetElements[1][0] - gPlanetElements[1][2]);
    map[2] = 0.0;
    map[3] = degToRad(gPlanetElements[3][0] - gPlanetElements[3][2]);

    //Compute perturbations
    dml = degToRad(7.7e-4*sin(4.1406+t*2.6227));
    dm = dml;

        dl = 3.13e-3*cos(2*mas-2*map[1]-2.587)+
             1.98e-3*cos(3*mas-3*map[1]+4.4768e-2)+
             1.36e-3*cos(mas-map[1]-2.0788)+
             9.6e-4*cos(3*mas-2*map[1]-2.3721)+
             8.2e-4*cos(map[3]-map[1]-3.6318);

        dr = 2.2501e-5*cos(2*mas-2*map[1]-1.01592)+
             1.9045e-5*cos(3*mas-3*map[1]+1.61577)+
             6.887e-6*cos(map[3]-map[1]-2.06106)+
             5.172e-6*cos(mas-map[1]-5.08065e-1)+
             3.62e-6*cos(5*mas-4*map[1]-1.81877)+
             3.283e-6*cos(4*mas-4*map[1]+1.10851)+
             3.074e-6*cos(2*map[3]-2*map[1]-9.62846e-1);

    computePlanetCoords(p, map[p], da, dhl, dl, dm, dml, dr, ds,
                        eclLong, eclLat, distance);

    //Corrections for internal coordinate system
    eclLat -= (PI/2);
    eclLong += PI;

    return Point3d(cos(eclLong) * sin(eclLat) * distance,
                   cos(eclLat) * distance,
                   -sin(eclLong) * sin(eclLat) * distance);
    };

    double getPeriod() const
    {
        return 224.7018;
    };

    double getBoundingRadius() const
    {
        return 1.089e+8 * BoundingRadiusSlack;
    };
};

class EarthOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
        double t, t2;
        double ls, ms;    // mean longitude and mean anomaly
        double s, nu, ea; // eccentricity, true anomaly, eccentric anomaly
        double a, b, a1, b1, c1, d1, e1, h1, dl, dr;
        double eclLong, distance;

        // Calculate the Julian centuries elapsed since 1900
        t = (jd - 2415020.0)/36525.0;

        t2 = t*t;
        a = 100.0021359*t;
        b = 360.*(a-(int)a);
        ls = 279.69668+.0003025*t2+b;
        ms = meanAnomalySun(t);
        s = .016751-.0000418*t-1.26e-07*t2;
        astro::anomaly(degToRad(ms), s, nu, ea);
        a = 62.55209472000015*t;
        b = 360*(a-(int)a);
        a1 = degToRad(153.23+b);
        a = 125.1041894*t;
        b = 360*(a-(int)a);
        b1 = degToRad(216.57+b);
        a = 91.56766028*t;
        b = 360*(a-(int)a);
        c1 = degToRad(312.69+b);
        a = 1236.853095*t;
        b = 360*(a-(int)a);
        d1 = degToRad(350.74-.00144*t2+b);
        e1 = degToRad(231.19+20.2*t);
        a = 183.1353208*t;
        b = 360*(a-(int)a);
        h1 = degToRad(353.4+b);
        dl = .00134*cos(a1)+.00154*cos(b1)+.002*cos(c1)+.00179*sin(d1)+
            .00178*sin(e1);
        dr = 5.43e-06*sin(a1)+1.575e-05*sin(b1)+1.627e-05*sin(c1)+
            3.076e-05*cos(d1)+9.27e-06*sin(h1);

        eclLong = nu+degToRad(ls-ms+dl) + PI;
        eclLong = pfmod(eclLong, TWOPI);
        distance = KM_PER_AU * (1.0000002*(1-s*cos(ea))+dr);

        // Correction for internal coordinate system
        eclLong += PI;
        
        return Point3d(-cos(eclLong) * distance,
                       0,
                       sin(eclLong) * distance);
    };

    double getPeriod() const
    {
        return 365.25;
    };

    double getBoundingRadius() const
    {
        return 1.52e+8 * BoundingRadiusSlack;
    };
};

class LunarOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
        double jd19, t, t2;
        double ld, ms, md, de, f, n, hp;
        double a, sa, sn, b, sb, c, sc, e, e2, l, g, w1, w2;
        double m1, m2, m3, m4, m5, m6;
        double eclLon, eclLat, horzPar, distance;
        double RA, dec;

        // Computation requires an abbreviated Julian day:
        // epoch January 0.5, 1900.
        jd19 = jd - 2415020.0;
        t = jd19/36525;
        t2 = t*t;

        m1 = jd19/27.32158213;
        m1 = 360.0*(m1-(int)m1);
        m2 = jd19/365.2596407;
        m2 = 360.0*(m2-(int)m2);
        m3 = jd19/27.55455094;
        m3 = 360.0*(m3-(int)m3);
        m4 = jd19/29.53058868;
        m4 = 360.0*(m4-(int)m4);
        m5 = jd19/27.21222039;
        m5 = 360.0*(m5-(int)m5);
        m6 = jd19/6798.363307;
        m6 = 360.0*(m6-(int)m6);

        ld = 270.434164+m1-(.001133-.0000019*t)*t2;
        ms = 358.475833+m2-(.00015+.0000033*t)*t2;
        md = 296.104608+m3+(.009192+.0000144*t)*t2;
        de = 350.737486+m4-(.001436-.0000019*t)*t2;
        f = 11.250889+m5-(.003211+.0000003*t)*t2;
        n = 259.183275-m6+(.002078+.000022*t)*t2;

        a = degToRad(51.2+20.2*t);
        sa = sin(a);
        sn = sin(degToRad(n));
        b = 346.56+(132.87-.0091731*t)*t;
        sb = .003964*sin(degToRad(b));
        c = degToRad(n+275.05-2.3*t);
        sc = sin(c);
        ld = ld+.000233*sa+sb+.001964*sn;
        ms = ms-.001778*sa;
        md = md+.000817*sa+sb+.002541*sn;
        f = f+sb-.024691*sn-.004328*sc;
        de = de+.002011*sa+sb+.001964*sn;
        e = 1-(.002495+7.52e-06*t)*t;
        e2 = e*e;

        ld = degToRad(ld);
        ms = degToRad(ms);
        n = degToRad(n);
        de = degToRad(de);
        f = degToRad(f);
        md = degToRad(md);

        l = 6.28875*sin(md)+1.27402*sin(2*de-md)+.658309*sin(2*de)+
            .213616*sin(2*md)-e*.185596*sin(ms)-.114336*sin(2*f)+
            .058793*sin(2*(de-md))+.057212*e*sin(2*de-ms-md)+
            .05332*sin(2*de+md)+.045874*e*sin(2*de-ms)+.041024*e*sin(md-ms);
        l = l-.034718*sin(de)-e*.030465*sin(ms+md)+.015326*sin(2*(de-f))-
            .012528*sin(2*f+md)-.01098*sin(2*f-md)+.010674*sin(4*de-md)+
            .010034*sin(3*md)+.008548*sin(4*de-2*md)-e*.00791*sin(ms-md+2*de)-
            e*.006783*sin(2*de+ms);
        l = l+.005162*sin(md-de)+e*.005*sin(ms+de)+.003862*sin(4*de)+
            e*.004049*sin(md-ms+2*de)+.003996*sin(2*(md+de))+
            .003665*sin(2*de-3*md)+e*.002695*sin(2*md-ms)+
            .002602*sin(md-2*(f+de))+e*.002396*sin(2*(de-md)-ms)-
            .002349*sin(md+de);
        l = l+e2*.002249*sin(2*(de-ms))-e*.002125*sin(2*md+ms)-
            e2*.002079*sin(2*ms)+e2*.002059*sin(2*(de-ms)-md)-
            .001773*sin(md+2*(de-f))-.001595*sin(2*(f+de))+
            e*.00122*sin(4*de-ms-md)-.00111*sin(2*(md+f))+.000892*sin(md-3*de);
        l = l-e*.000811*sin(ms+md+2*de)+e*.000761*sin(4*de-ms-2*md)+
            e2*.000704*sin(md-2*(ms+de))+e*.000693*sin(ms-2*(md-de))+
            e*.000598*sin(2*(de-f)-ms)+.00055*sin(md+4*de)+.000538*sin(4*md)+
            e*.000521*sin(4*de-ms)+.000486*sin(2*md-de);
        l = l+e2*.000717*sin(md-2*ms);
        eclLon = ld+degToRad(l);
        eclLon = pfmod(eclLon, TWOPI);

        g = 5.12819*sin(f)+.280606*sin(md+f)+.277693*sin(md-f)+
            .173238*sin(2*de-f)+.055413*sin(2*de+f-md)+.046272*sin(2*de-f-md)+
            .032573*sin(2*de+f)+.017198*sin(2*md+f)+.009267*sin(2*de+md-f)+
            .008823*sin(2*md-f)+e*.008247*sin(2*de-ms-f);
        g = g+.004323*sin(2*(de-md)-f)+.0042*sin(2*de+f+md)+
            e*.003372*sin(f-ms-2*de)+e*.002472*sin(2*de+f-ms-md)+
            e*.002222*sin(2*de+f-ms)+e*.002072*sin(2*de-f-ms-md)+
            e*.001877*sin(f-ms+md)+.001828*sin(4*de-f-md)-e*.001803*sin(f+ms)-
            .00175*sin(3*f);
        g = g+e*.00157*sin(md-ms-f)-.001487*sin(f+de)-e*.001481*sin(f+ms+md)+
            e*.001417*sin(f-ms-md)+e*.00135*sin(f-ms)+.00133*sin(f-de)+
            .001106*sin(f+3*md)+.00102*sin(4*de-f)+.000833*sin(f+4*de-md)+
            .000781*sin(md-3*f)+.00067*sin(f+4*de-2*md);
        g = g+.000606*sin(2*de-3*f)+.000597*sin(2*(de+md)-f)+
            e*.000492*sin(2*de+md-ms-f)+.00045*sin(2*(md-de)-f)+
            .000439*sin(3*md-f)+.000423*sin(f+2*(de+md))+
            .000422*sin(2*de-f-3*md)-e*.000367*sin(ms+f+2*de-md)-
            e*.000353*sin(ms+f+2*de)+.000331*sin(f+4*de);
        g = g+e*.000317*sin(2*de+f-ms+md)+e2*.000306*sin(2*(de-ms)-f)-
            .000283*sin(md+3*f);
        w1 = .0004664*cos(n);
        w2 = .0000754*cos(c);
        eclLat = degToRad(g)*(1-w1-w2);

        hp = .950724+.051818*cos(md)+.009531*cos(2*de-md)+.007843*cos(2*de)+
         .002824*cos(2*md)+.000857*cos(2*de+md)+e*.000533*cos(2*de-ms)+
         e*.000401*cos(2*de-md-ms)+e*.00032*cos(md-ms)-.000271*cos(de)-
         e*.000264*cos(ms+md)-.000198*cos(2*f-md);
        hp = hp+.000173*cos(3*md)+.000167*cos(4*de-md)-e*.000111*cos(ms)+
         .000103*cos(4*de-2*md)-.000084*cos(2*md-2*de)-
         e*.000083*cos(2*de+ms)+.000079*cos(2*de+2*md)+.000072*cos(4*de)+
         e*.000064*cos(2*de-ms+md)-e*.000063*cos(2*de+ms-md)+
         e*.000041*cos(ms+de);
        hp = hp+e*.000035*cos(2*md-ms)-.000033*cos(3*md-2*de)-
         .00003*cos(md+de)-.000029*cos(2*(f-de))-e*.000029*cos(2*md+ms)+
         e2*.000026*cos(2*(de-ms))-.000023*cos(2*(f-de)+md)+
         e*.000019*cos(4*de-ms-md);
        horzPar = degToRad(hp);

        // At this point we have values of ecliptic longitude, latitude and
        // horizontal parallax (eclLong, eclLat, horzPar) in radians.

        // Now compute distance using horizontal parallax.
        distance = 6378.14 / sin(horzPar);

        // Finally convert eclLat, eclLon to RA, Dec.
        EclipticToEquatorial(t, eclLat, eclLon, RA, dec);

        // RA and Dec are referred to the equinox of date; we want to use
        // the J2000 equinox instead.  A better idea would be to directly 
        // compute the position of the Moon in this coordinate system, but
        // this was easier.
        EpochConvert(jd, astro::J2000, RA, dec, RA, dec);

        // Corrections for internal coordinate system
        dec -= (PI/2);
        RA += PI;

        return Point3d(cos(RA) * sin(dec) * distance,
                       cos(dec) * distance,
                       -sin(RA) * sin(dec) * distance);
    };

    double getPeriod() const
    {
        return 27.321661;
    };

    double getBoundingRadius() const
    {
        return 405504 * BoundingRadiusSlack;
    };
};

class MarsOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
    const int p = 2;  //Planet 2
    vector<int> pList;
    double t;
    double map[4], mas, a;
    double dl, dr, dml, ds, dm, da, dhl;
    double eclLong, eclLat, distance;    //heliocentric longitude, latitude, distance

    dl = dr = dml = ds = dm = da = dhl = 0.0;

    //Calculate the Julian centuries elapsed since 1900
        t = (jd - 2415020.0)/36525.0;

    mas = meanAnomalySun(t);

    //Specify which planets we must compute elements for
    pList.push_back(1);
    pList.push_back(2);
    pList.push_back(3);
    computePlanetElements(t, pList);

    //Compute necessary planet mean anomalies
    map[0] = 0.0;
    map[1] = degToRad(gPlanetElements[1][0] - gPlanetElements[1][2]);
    map[2] = degToRad(gPlanetElements[2][0] - gPlanetElements[2][2]);
    map[3] = degToRad(gPlanetElements[3][0] - gPlanetElements[3][2]);

    //Compute perturbations
    a = 3*map[3]-8*map[2]+4*mas;
    dml = degToRad(-1*(1.133e-2*sin(a)+9.33e-3*cos(a)));
    dm = dml;
    dl = 7.05e-3*cos(map[3]-map[2]-8.5448e-1)+
             6.07e-3*cos(2*map[3]-map[2]-3.2873)+
             4.45e-3*cos(2*map[3]-2*map[2]-3.3492)+
             3.88e-3*cos(mas-2*map[2]+3.5771e-1)+
             2.38e-3*cos(mas-map[2]+6.1256e-1)+
             2.04e-3*cos(2*mas-3*map[2]+2.7688)+
             1.77e-3*cos(3*map[2]-map[2-1]-1.0053)+
             1.36e-3*cos(2*mas-4*map[2]+2.6894)+
             1.04e-3*cos(map[3]+3.0749e-1);

    dr = 5.3227e-5*cos(map[3]-map[2]+7.17864e-1)+
             5.0989e-5*cos(2*map[3]-2*map[2]-1.77997)+
             3.8278e-5*cos(2*map[3]-map[2]-1.71617)+
             1.5996e-5*cos(mas-map[2]-9.69618e-1)+
             1.4764e-5*cos(2*mas-3*map[2]+1.19768)+
             8.966e-6*cos(map[3]-2*map[2]+7.61225e-1);
    dr += 7.914e-6*cos(3*map[3]-2*map[2]-2.43887)+
             7.004e-6*cos(2*map[3]-3*map[2]-1.79573)+
             6.62e-6*cos(mas-2*map[2]+1.97575)+
             4.93e-6*cos(3*map[3]-3*map[2]-1.33069)+
             4.693e-6*cos(3*mas-5*map[2]+3.32665)+
             4.571e-6*cos(2*mas-4*map[2]+4.27086)+
             4.409e-6*cos(3*map[3]-map[2]-2.02158);

    computePlanetCoords(p, map[p], da, dhl, dl, dm, dml, dr, ds,
                        eclLong, eclLat, distance);

    //Corrections for internal coordinate system
    eclLat -= (PI/2);
    eclLong += PI;

    return Point3d(cos(eclLong) * sin(eclLat) * distance,
                   cos(eclLat) * distance,
                   -sin(eclLong) * sin(eclLat) * distance);
    };

    double getPeriod() const
    {
        return 689.998725;
    };

    double getBoundingRadius() const
    {
        return 2.49e+8 * BoundingRadiusSlack;
    };
};

class JupiterOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
    const int p = 3;  //Planet 3
    vector<int> pList(1, p);
    double t, map;
    double dl, dr, dml, ds, dm, da, dhl, s;
    double dp;
    double x1, x2, x3, x4, x5, x6, x7;
    double sx3, cx3, s2x3, c2x3;
    double sx5, cx5, s2x5;
    double sx6;
    double sx7, cx7, s2x7, c2x7, s3x7, c3x7, s4x7, c4x7, c5x7;
    double eclLong, eclLat, distance;    //heliocentric longitude, latitude, distance

    dl = dr = dml = ds = dm = da = dhl = 0.0;

    //Calculate the Julian centuries elapsed since 1900
        t = (jd - 2415020.0)/36525.0;

    computePlanetElements(t, pList);

    map = degToRad(gPlanetElements[p][0] - gPlanetElements[p][2]);

    //Compute perturbations
    s = gPlanetElements[p][3];
    auxJSun(t, &x1, &x2, &x3, &x4, &x5, &x6);
    x7 = x3-x2;
    sx3 = sin(x3);
    cx3 = cos(x3);
    s2x3 = sin(2*x3);
    c2x3 = cos(2*x3);
    sx5 = sin(x5);
    cx5 = cos(x5);
    s2x5 = sin(2*x5);
    sx6 = sin(x6);
    sx7 = sin(x7);
    cx7 = cos(x7);
    s2x7 = sin(2*x7);
    c2x7 = cos(2*x7);
    s3x7 = sin(3*x7);
    c3x7 = cos(3*x7);
    s4x7 = sin(4*x7);
    c4x7 = cos(4*x7);
    c5x7 = cos(5*x7);
    dml = (3.31364e-1-(1.0281e-2+4.692e-3*x1)*x1)*sx5+
              (3.228e-3-(6.4436e-2-2.075e-3*x1)*x1)*cx5-
              (3.083e-3+(2.75e-4-4.89e-4*x1)*x1)*s2x5+
              2.472e-3*sx6+1.3619e-2*sx7+1.8472e-2*s2x7+6.717e-3*s3x7+
              2.775e-3*s4x7+6.417e-3*s2x7*sx3+
              (7.275e-3-1.253e-3*x1)*sx7*sx3+
              2.439e-3*s3x7*sx3-(3.5681e-2+1.208e-3*x1)*sx7*cx3;
    dml += -3.767e-3*c2x7*sx3-(3.3839e-2+1.125e-3*x1)*cx7*sx3-
              4.261e-3*s2x7*cx3+
              (1.161e-3*x1-6.333e-3)*cx7*cx3+
              2.178e-3*cx3-6.675e-3*c2x7*cx3-2.664e-3*c3x7*cx3-
              2.572e-3*sx7*s2x3-3.567e-3*s2x7*s2x3+2.094e-3*cx7*c2x3+
              3.342e-3*c2x7*c2x3;
    dml = degToRad(dml);
    ds = (3606+(130-43*x1)*x1)*sx5+(1289-580*x1)*cx5-6764*sx7*sx3-
             1110*s2x7*sx3-224*s3x7*sx3-204*sx3+(1284+116*x1)*cx7*sx3+
             188*c2x7*sx3+(1460+130*x1)*sx7*cx3+224*s2x7*cx3-817*cx3+
             6074*cx3*cx7+992*c2x7*cx3+
             508*c3x7*cx3+230*c4x7*cx3+108*c5x7*cx3;
    ds += -(956+73*x1)*sx7*s2x3+448*s2x7*s2x3+137*s3x7*s2x3+
             (108*x1-997)*cx7*s2x3+480*c2x7*s2x3+148*c3x7*s2x3+
             (99*x1-956)*sx7*c2x3+490*s2x7*c2x3+
             158*s3x7*c2x3+179*c2x3+(1024+75*x1)*cx7*c2x3-
             437*c2x7*c2x3-132*c3x7*c2x3;
    ds *= 1e-7;
    dp = (7.192e-3-3.147e-3*x1)*sx5-4.344e-3*sx3+
             (x1*(1.97e-4*x1-6.75e-4)-2.0428e-2)*cx5+
             3.4036e-2*cx7*sx3+(7.269e-3+6.72e-4*x1)*sx7*sx3+
             5.614e-3*c2x7*sx3+2.964e-3*c3x7*sx3+3.7761e-2*sx7*cx3+
             6.158e-3*s2x7*cx3-
             6.603e-3*cx7*cx3-5.356e-3*sx7*s2x3+2.722e-3*s2x7*s2x3+
             4.483e-3*cx7*s2x3-2.642e-3*c2x7*s2x3+4.403e-3*sx7*c2x3-
             2.536e-3*s2x7*c2x3+5.547e-3*cx7*c2x3-2.689e-3*c2x7*c2x3;
    dm = dml-(degToRad(dp)/s);
    da = 205*cx7-263*cx5+693*c2x7+312*c3x7+147*c4x7+299*sx7*sx3+
             181*c2x7*sx3+204*s2x7*cx3+111*s3x7*cx3-337*cx7*cx3-
             111*c2x7*cx3;
    da *= 1e-6;

    computePlanetCoords(p, map, da, dhl, dl, dm, dml, dr, ds,
                        eclLong, eclLat, distance);

    //Corrections for internal coordinate system
    eclLat -= (PI/2);
    eclLong += PI;

    return Point3d(cos(eclLong) * sin(eclLat) * distance,
                   cos(eclLat) * distance,
                   -sin(eclLong) * sin(eclLat) * distance);
    };

    double getPeriod() const
    {
        return 4332.66855;
    };

    double getBoundingRadius() const
    {
        return 8.16e+8 * BoundingRadiusSlack;
    };
};

class SaturnOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
    const int p = 4;  //Planet 4
    vector<int> pList(1, p);
    double t, map;
    double dl, dr, dml, ds, dm, da, dhl, s;
    double dp;
    double x1, x2, x3, x4, x5, x6, x7, x8;
    double sx3, cx3, s2x3, c2x3, s3x3, c3x3, s4x3, c4x3;
    double sx5, cx5, s2x5, c2x5;
    double sx6;
    double sx7, cx7, s2x7, c2x7, s3x7, c3x7, s4x7, c4x7, c5x7, s5x7;
    double s2x8, c2x8, s3x8, c3x8;
    double eclLong, eclLat, distance;    //heliocentric longitude, latitude, distance

    dl = dr = dml = ds = dm = da = dhl = 0.0;

    //Calculate the Julian centuries elapsed since 1900
        t = (jd - 2415020.0)/36525.0;

    computePlanetElements(t, pList);

    map = degToRad(gPlanetElements[p][0] - gPlanetElements[p][2]);

    //Compute perturbations
    s = gPlanetElements[p][3];
    auxJSun(t, &x1, &x2, &x3, &x4, &x5, &x6);
    x7 = x3-x2;
    sx3 = sin(x3);
    cx3 = cos(x3);
    s2x3 = sin(2*x3);
    c2x3 = cos(2*x3);
    sx5 = sin(x5);
    cx5 = cos(x5);
    s2x5 = sin(2*x5);
    sx6 = sin(x6);
    sx7 = sin(x7);
    cx7 = cos(x7);
    s2x7 = sin(2*x7);
    c2x7 = cos(2*x7);
    s3x7 = sin(3*x7);
    c3x7 = cos(3*x7);
    s4x7 = sin(4*x7);
    c4x7 = cos(4*x7);
    c5x7 = cos(5*x7);
    s3x3 = sin(3*x3);
    c3x3 = cos(3*x3);
    s4x3 = sin(4*x3);
    c4x3 = cos(4*x3);
    c2x5 = cos(2*x5);
    s5x7 = sin(5*x7);
    x8 = x4-x3;
    s2x8 = sin(2*x8);
    c2x8 = cos(2*x8);
    s3x8 = sin(3*x8);
    c3x8 = cos(3*x8);
    dml = 7.581e-3*s2x5-7.986e-3*sx6-1.48811e-1*sx7-4.0786e-2*s2x7-
              (8.14181e-1-(1.815e-2-1.6714e-2*x1)*x1)*sx5-
              (1.0497e-2-(1.60906e-1-4.1e-3*x1)*x1)*cx5-1.5208e-2*s3x7-
              6.339e-3*s4x7-6.244e-3*sx3-1.65e-2*s2x7*sx3+
              (8.931e-3+2.728e-3*x1)*sx7*sx3-5.775e-3*s3x7*sx3+
              (8.1344e-2+3.206e-3*x1)*cx7*sx3+1.5019e-2*c2x7*sx3;
    dml += (8.5581e-2+2.494e-3*x1)*sx7*cx3+1.4394e-2*c2x7*cx3+
              (2.5328e-2-3.117e-3*x1)*cx7*cx3+
              6.319e-3*c3x7*cx3+6.369e-3*sx7*s2x3+9.156e-3*s2x7*s2x3+
              7.525e-3*s3x8*s2x3-5.236e-3*cx7*c2x3-7.736e-3*c2x7*c2x3-
              7.528e-3*c3x8*c2x3;
    dml = degToRad(dml);
    ds = (-7927+(2548+91*x1)*x1)*sx5+(13381+(1226-253*x1)*x1)*cx5+
             (248-121*x1)*s2x5-(305+91*x1)*c2x5+412*s2x7+12415*sx3+
             (390-617*x1)*sx7*sx3+(165-204*x1)*s2x7*sx3+26599*cx7*sx3-
             4687*c2x7*sx3-1870*c3x7*sx3-821*c4x7*sx3-
             377*c5x7*sx3+497*c2x8*sx3+(163-611*x1)*cx3;
    ds += -12696*sx7*cx3-4200*s2x7*cx3-1503*s3x7*cx3-619*s4x7*cx3-
             268*s5x7*cx3-(282+1306*x1)*cx7*cx3+(-86+230*x1)*c2x7*cx3+
             461*s2x8*cx3-350*s2x3+(2211-286*x1)*sx7*s2x3-
             2208*s2x7*s2x3-568*s3x7*s2x3-346*s4x7*s2x3-
             (2780+222*x1)*cx7*s2x3+(2022+263*x1)*c2x7*s2x3+248*c3x7*s2x3+
             242*s3x8*s2x3+467*c3x8*s2x3-490*c2x3-(2842+279*x1)*sx7*c2x3;
    ds += (128+226*x1)*s2x7*c2x3+224*s3x7*c2x3+
             (-1594+282*x1)*cx7*c2x3+(2162-207*x1)*c2x7*c2x3+
             561*c3x7*c2x3+343*c4x7*c2x3+469*s3x8*c2x3-242*c3x8*c2x3-
             205*sx7*s3x3+262*s3x7*s3x3+208*cx7*c3x3-271*c3x7*c3x3-
             382*c3x7*s4x3-376*s3x7*c4x3;
    ds *= 1e-7;
    dp = (7.7108e-2+(7.186e-3-1.533e-3*x1)*x1)*sx5-7.075e-3*sx7+
             (4.5803e-2-(1.4766e-2+5.36e-4*x1)*x1)*cx5-7.2586e-2*cx3-
             7.5825e-2*sx7*sx3-2.4839e-2*s2x7*sx3-8.631e-3*s3x7*sx3-
             1.50383e-1*cx7*cx3+2.6897e-2*c2x7*cx3+1.0053e-2*c3x7*cx3-
             (1.3597e-2+1.719e-3*x1)*sx7*s2x3+1.1981e-2*s2x7*c2x3;
    dp += -(7.742e-3-1.517e-3*x1)*cx7*s2x3+
             (1.3586e-2-1.375e-3*x1)*c2x7*c2x3-
             (1.3667e-2-1.239e-3*x1)*sx7*c2x3+
             (1.4861e-2+1.136e-3*x1)*cx7*c2x3-
             (1.3064e-2+1.628e-3*x1)*c2x7*c2x3;
    dm = dml-(degToRad(dp)/s);
    da = 572*sx5-1590*s2x7*cx3+2933*cx5-647*s3x7*cx3+33629*cx7-
             344*s4x7*cx3-3081*c2x7+2885*cx7*cx3-1423*c3x7+
             (2172+102*x1)*c2x7*cx3-671*c4x7+296*c3x7*cx3-320*c5x7-
             267*s2x7*s2x3+1098*sx3-778*cx7*s2x3-2812*sx7*sx3;
    da += 495*c2x7*s2x3+688*s2x7*sx3+250*c3x7*s2x3-393*s3x7*sx3-
             856*sx7*c2x3-228*s4x7*sx3+441*s2x7*c2x3+2138*cx7*sx3+
             296*c2x7*c2x3-999*c2x7*sx3+211*c3x7*c2x3-642*c3x7*sx3-
             427*sx7*s3x3-325*c4x7*sx3+398*s3x7*s3x3-890*cx3+
             344*cx7*c3x3+2206*sx7*cx3-427*c3x7*c3x3;
    da *= 1e-6;
    dhl = 7.47e-4*cx7*sx3+1.069e-3*cx7*cx3+2.108e-3*s2x7*s2x3+
              1.261e-3*c2x7*s2x3+1.236e-3*s2x7*c2x3-2.075e-3*c2x7*c2x3;
    dhl = degToRad(dhl);

    computePlanetCoords(p, map, da, dhl, dl, dm, dml, dr, ds,
                        eclLong, eclLat, distance);

    //Corrections for internal coordinate system
    eclLat -= (PI/2);
    eclLong += PI;

    return Point3d(cos(eclLong) * sin(eclLat) * distance,
                   cos(eclLat) * distance,
                   -sin(eclLong) * sin(eclLat) * distance);
    };

    double getPeriod() const
    {
        return 10759.42493;
    };

    double getBoundingRadius() const
    {
        return 1.50e+9 * BoundingRadiusSlack;
    };
};

class UranusOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
    const int p = 5;  //Planet 5
    vector<int> pList(1, p);
    double t, map;
    double dl, dr, dml, ds, dm, da, dhl, s;
    double dp;
    double x1, x2, x3, x4, x5, x6;
    double x8, x9, x10, x11, x12;
    double sx4, cx4, s2x4, c2x4;
    double sx9, cx9, s2x9, c2x9;
    double sx11, cx11;
    double eclLong, eclLat, distance;    //heliocentric longitude, latitude, distance

    dl = dr = dml = ds = dm = da = dhl = 0.0;

    //Calculate the Julian centuries elapsed since 1900
        t = (jd - 2415020.0)/36525.0;

    computePlanetElements(t, pList);

    map = degToRad(gPlanetElements[p][0] - gPlanetElements[p][2]);

    //Compute perturbations
    s = gPlanetElements[p][3];
    auxJSun(t, &x1, &x2, &x3, &x4, &x5, &x6);
    x8 = pfmod(1.46205+3.81337*t, TWOPI);
    x9 = 2*x8-x4;
    sx9 = sin(x9);
    cx9 = cos(x9);
    s2x9 = sin(2*x9);
    c2x9 = cos(2*x9);
    x10 = x4-x2;
    x11 = x4-x3;
    x12 = x8-x4;
    dml = (8.64319e-1-1.583e-3*x1)*sx9+(8.2222e-2-6.833e-3*x1)*cx9+
              3.6017e-2*s2x9-3.019e-3*c2x9+8.122e-3*sin(x6);
    dml = degToRad(dml);
    dp = 1.20303e-1*sx9+6.197e-3*s2x9+(1.9472e-2-9.47e-4*x1)*cx9;
    dm = dml-(degToRad(dp)/s);
    ds = (163*x1-3349)*sx9+20981*cx9+1311*c2x9;
    ds *= 1e-7;
    da = -3.825e-3*cx9;
    dl = (1.0122e-2-9.88e-4*x1)*sin(x4+x11)+
             (-3.8581e-2+(2.031e-3-1.91e-3*x1)*x1)*cos(x4+x11)+
             (3.4964e-2-(1.038e-3-8.68e-4*x1)*x1)*cos(2*x4+x11)+
             5.594e-3*sin(x4+3*x12)-1.4808e-2*sin(x10)-
             5.794e-3*sin(x11)+2.347e-3*cos(x11)+9.872e-3*sin(x12)+
             8.803e-3*sin(2*x12)-4.308e-3*sin(3*x12);
    sx11 = sin(x11);
    cx11 = cos(x11);
    sx4 = sin(x4);
    cx4 = cos(x4);
    s2x4 = sin(2*x4);
    c2x4 = cos(2*x4);
    dhl = (4.58e-4*sx11-6.42e-4*cx11-5.17e-4*cos(4*x12))*sx4-
              (3.47e-4*sx11+8.53e-4*cx11+5.17e-4*sin(4*x11))*cx4+
              4.03e-4*(cos(2*x12)*s2x4+sin(2*x12)*c2x4);
    dhl = degToRad(dhl);

    dr = -25948+4985*cos(x10)-1230*cx4+3354*cos(x11)+904*cos(2*x12)+
             894*(cos(x12)-cos(3*x12))+(5795*cx4-1165*sx4+1388*c2x4)*sx11+
             (1351*cx4+5702*sx4+1388*s2x4)*cos(x11);
    dr *= 1e-6;

    computePlanetCoords(p, map, da, dhl, dl, dm, dml, dr, ds,
                        eclLong, eclLat, distance);

    //Corrections for internal coordinate system
    eclLat -= (PI/2);
    eclLong += PI;

    return Point3d(cos(eclLong) * sin(eclLat) * distance,
                   cos(eclLat) * distance,
                   -sin(eclLong) * sin(eclLat) * distance);
    };

    double getPeriod() const
    {
        return 30686.07698;
    };

    double getBoundingRadius() const
    {
        return 3.01e+9 * BoundingRadiusSlack;
    };
};

class NeptuneOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
    const int p = 6;  //Planet 6
    vector<int> pList(1, p);
    double t, map;
    double dl, dr, dml, ds, dm, da, dhl, s;
    double dp;
    double x1, x2, x3, x4, x5, x6;
    double x8, x9, x10, x11, x12;
    double sx8, cx8;
    double sx9, cx9, s2x9, c2x9;
    double s2x12, c2x12;
    double eclLong, eclLat, distance;    //heliocentric longitude, latitude, distance

    dl = dr = dml = ds = dm = da = dhl = 0.0;

    //Calculate the Julian centuries elapsed since 1900
        t = (jd - 2415020.0)/36525.0;

    computePlanetElements(t, pList);

    map = degToRad(gPlanetElements[p][0] - gPlanetElements[p][2]);

    //Compute perturbations
    s = gPlanetElements[p][3];
        auxJSun(t, &x1, &x2, &x3, &x4, &x5, &x6);
    x8 = pfmod(1.46205+3.81337*t, TWOPI);
    x9 = 2*x8-x4;
        sx9 = sin(x9);
        cx9 = cos(x9);
    s2x9 = sin(2*x9);
        c2x9 = cos(2*x9);
        x10 = x8-x2;
        x11 = x8-x3;
        x12 = x8-x4;
        dml = (1.089e-3*x1-5.89833e-1)*sx9+(4.658e-3*x1-5.6094e-2)*cx9-
              2.4286e-2*s2x9;
        dml = degToRad(dml);
        dp = 2.4039e-2*sx9-2.5303e-2*cx9+6.206e-3*s2x9-5.992e-3*c2x9;
        dm = dml-(degToRad(dp)/s);
        ds = 4389*sx9+1129*s2x9+4262*cx9+1089*c2x9;
        ds *= 1e-7;
        da = 8189*cx9-817*sx9+781*c2x9;
        da *= 1e-6;
        s2x12 = sin(2*x12);
        c2x12 = cos(2*x12);
        sx8 = sin(x8);
        cx8 = cos(x8);
        dl = -9.556e-3*sin(x10)-5.178e-3*sin(x11)+2.572e-3*s2x12-
             2.972e-3*c2x12*sx8-2.833e-3*s2x12*cx8;
        dhl = 3.36e-4*c2x12*sx8+3.64e-4*s2x12*cx8;
        dhl = degToRad(dhl);
        dr = -40596+4992*cos(x10)+2744*cos(x11)+2044*cos(x12)+1051*c2x12;
        dr *= 1e-6;

    computePlanetCoords(p, map, da, dhl, dl, dm, dml, dr, ds,
                        eclLong, eclLat, distance);

    //Corrections for internal coordinate system
    eclLat -= (PI/2);
    eclLong += PI;

    return Point3d(cos(eclLong) * sin(eclLat) * distance,
                   cos(eclLat) * distance,
                   -sin(eclLong) * sin(eclLat) * distance);
    };

    double getPeriod() const
    {
        return 60190.64325;
    };

    double getBoundingRadius() const
    {
        return 4.54e+9 * BoundingRadiusSlack;
    };
};

class PlutoOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
    const int p = 7;  //Planet 7
    vector<int> pList(1, p);
    double t, map;
    double dl, dr, dml, ds, dm, da, dhl;
    double eclLong, eclLat, distance;    //heliocentric longitude, latitude, distance

    dl = dr = dml = ds = dm = da = dhl = 0.0;

    //Calculate the Julian centuries elapsed since 1900
    t = (jd - 2415020.0)/36525.0;

    computePlanetElements(t, pList);

    map = degToRad(gPlanetElements[p][0] - gPlanetElements[p][2]);

    computePlanetCoords(p, map, da, dhl, dl, dm, dml, dr, ds,
                        eclLong, eclLat, distance);

    //Corrections for internal coordinate system
    eclLat -= PI / 2;
    eclLong += PI;

    return Point3d(cos(eclLong) * sin(eclLat) * distance,
                   cos(eclLat) * distance,
                   -sin(eclLong) * sin(eclLat) * distance);
    };

    double getPeriod() const
    {
        return 90779.235;
    };

    double getBoundingRadius() const
    {
        return 7.38e+9 * BoundingRadiusSlack;
    };
};


// Compute for mean anomaly M the point on the ellipse with
// semimajor axis a and eccentricity e.  This helper function assumes
// a low eccentricity; orbit.cpp has functions appropriate for solving
// Kepler's equation for larger values of e.
static Point3d ellipsePosition(double a, double e, double M)
{
    // Solve Kepler's equation--for a low eccentricity orbit, just a few
    // iterations is enough.
    double E = M;
    for (int k = 0; k < 3; k++)
        E = M + e * sin(E);

    return Point3d(a * (cos(E) - e),
                   0.0,
                   a * sqrt(1 - square(e)) * -sin(E));
}


class PhobosOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
        double epoch = 2433283.0 - 0.5; // 00:00 1 Jan 1950
        double a     = 9380.0;
        double e     = 0.0151;
        double w0    = 150.247;
        double M0    =  92.474;
        double i     =   1.075;
        double node0 = 164.931;
        double n     = 1128.8444155;
        double Pw    =  1.131;
        double Pnode =  2.262;

        double refplane_RA  = 317.724;
        double refplane_Dec =  52.924;
        double marspole_RA  = 317.681;
        double marspole_Dec =  52.886;

        double t = jd - epoch;
        t += 10.5 / 1440.0;     // light time correction?
        double T = t / 365.25;

        double dnode = 360.0 / Pnode;
        double dw    = 360.0 / Pw;
        double node = degToRad(node0 + T * dnode);
        double w    = degToRad(w0 + T * dw - T * dnode);
        double M    = degToRad(M0 + t * n  - T * dw);

        Point3d p = ellipsePosition(a, e, M);

        // Orientation of the orbital plane with respect to the Laplacian plane
        Mat3d Rorbit     = (Mat3d::yrotation(node) *
                            Mat3d::xrotation(degToRad(i)) *
                            Mat3d::yrotation(w));

        // Rotate to the Earth's equatorial plane
        double N = degToRad(refplane_RA);
        double J = degToRad(90 - refplane_Dec);
        Mat3d RLaplacian = (Mat3d::yrotation( N) *
                            Mat3d::xrotation( J) *
                            Mat3d::yrotation(-N));

        // Rotate to the Martian equatorial plane
        N = degToRad(marspole_RA);
        J = degToRad(90 - marspole_Dec);
        Mat3d RMars_eq   = (Mat3d::yrotation( N) *
                            Mat3d::xrotation(-J) *
                            Mat3d::yrotation(-N));

        return RMars_eq * (RLaplacian * (Rorbit * p));
    }
    
    double getPeriod() const
    {
        return 0.319;
    }

    double getBoundingRadius() const
    {
        return 9380 * BoundingRadiusSlack;
    }
};


class DeimosOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
        double epoch = 2433283.0 - 0.5;
        double a     = 23460.0;
        double e     = 0.0002;
        double w0    = 290.496;
        double M0    = 296.230;
        double i     = 1.793;
        double node0 = 339.600;
        double n     = 285.1618919;
        double Pw    = 26.892;
        double Pnode = 54.536;

        double refplane_RA  = 316.700;
        double refplane_Dec =  53.564;
        double marspole_RA  = 317.681;
        double marspole_Dec =  52.886;

        double t = jd - epoch;
        t += 10.5 / 1440.0;     // light time correction?
        double T = t / 365.25;

        double dnode = 360.0 / Pnode;
        double dw    = 360.0 / Pw;
        double node = degToRad(node0 + T * dnode);
        double w    = degToRad(w0 + T * dw - T * dnode);
        double M    = degToRad(M0 + t * n  - T * dw);

        Point3d p = ellipsePosition(a, e, M);

        // Orientation of the orbital plane with respect to the Laplacian plane
        Mat3d Rorbit     = (Mat3d::yrotation(node) *
                            Mat3d::xrotation(degToRad(i)) *
                            Mat3d::yrotation(w));

        // Rotate to the Earth's equatorial plane
        double N = degToRad(refplane_RA);
        double J = degToRad(90 - refplane_Dec);
        Mat3d RLaplacian = (Mat3d::yrotation( N) *
                            Mat3d::xrotation( J) *
                            Mat3d::yrotation(-N));

        // Rotate to the Martian equatorial plane
        N = degToRad(marspole_RA);
        J = degToRad(90 - marspole_Dec);
        Mat3d RMars_eq   = (Mat3d::yrotation( N) *
                            Mat3d::xrotation(-J) *
                            Mat3d::yrotation(-N));

        return RMars_eq * (RLaplacian * (Rorbit * p));
    }

#if 0
    // More accurate orbit calculation for Mars from _Explanatory
    // Supplement to the Astronomical Almanac_  There's still a bug in
    // this routine however.
    Point3d computePosition(double jd) const
    {
        double d = jd - 2441266.5;  // days since 11 Nov 1971
        double D = d / 365.25;      // years

        double a = 1.56828e-4 * KM_PER_AU;
        double n = 285.161888;
        double e = 0.0004;
        double gamma = degToRad(1.79);
        double theta = degToRad(240.38 - 0.01801 * d);

        double h = degToRad(196.55 - 0.01801 * d);
        double L = degToRad(28.96 + n * d - 0.27 * sin(h));
        double P = degToRad(111.7 + 0.01798 * d);

        // N and J give the orientation of the Laplacian plane with respect
        // to the Earth's equator and equinox.
        //    N - longitude of ascending node
        //    J - inclination
        double N = degToRad(46.37 - 0.0014 * D);
        double J = degToRad(36.62 + 0.0008 * D);

        // Compute the mean anomaly
        double M = L - P;
        
        // Solve Kepler's equation--for a low eccentricity orbit, just a few
        // iterations is enough.
        double E = M;
        for (int i = 0; i < 3; i++)
            E = M + e * sin(E);

        // Compute the position in the orbital plane (y = 0)
        double x = a * (cos(E) - e);
        double z = a * sqrt(1 - square(e)) * -sin(E);

        // Orientation of the orbital plane with respect to the Laplacian plane
        Mat3d Rorbit     = (Mat3d::yrotation(theta) *
                            Mat3d::xrotation(gamma) *
                            Mat3d::yrotation(P - theta));

        Mat3d RLaplacian = (Mat3d::yrotation( N) *
                            Mat3d::xrotation(-J) *
                            Mat3d::yrotation(-N));

        double marspole_RA  = 317.681;
        double marspole_Dec =  52.886;
        double Nm = degToRad(marspole_RA + 90);
        double Jm = degToRad(90 - marspole_Dec);
        Mat3d RMars_eq   = (Mat3d::yrotation( Nm) *
                            Mat3d::xrotation( Jm) *
                            Mat3d::yrotation(-Nm));

        // Celestia wants the position of a satellite with respect to the
        // equatorial plane of the planet it orbits.

        // to Laplacian...
        Point3d p = Rorbit * Point3d(x, 0.0, z);

        // to Earth equatorial...
        p = RLaplacian * p;

        // to Mars equatorial...
        return RMars_eq * p;
    }
#endif

    double getPeriod() const
    {
        return 1.262441;
    }

    double getBoundingRadius() const
    {
        return 23462 * BoundingRadiusSlack;
    }
};


// static const double JupAscendingNode = degToRad(20.453422);
static const double JupAscendingNode = degToRad(22.203);

class IoOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
    //Computation will yield latitude(L), longitude(B) and distance(R) relative to Jupiter
    double t;
    double l1, l2, l3, l4;
    double p1, p2, p3, p4;
    double w1, w2, w3, w4;
    double gamma, phi, psi, G, Gp;
    double sigma, L, B, R;
    double T, P;

    // Epoch for Galilean satellites is 1976.0 Aug 10
    t = jd - 2443000.5;

    ComputeGalileanElements(t,
                            l1, l2, l3, l4,
                            p1, p2, p3, p4,
                            w1, w2, w3, w4,
                            gamma, phi, psi, G, Gp);

    // Calculate periodic terms for longitude
    sigma = 0.47259*sin(2*(l1 - l2)) - 0.03478*sin(p3 - p4)
          + 0.01081*sin(l2 - 2*l3 + p3) + 7.38e-3*sin(phi)
          + 7.13e-3*sin(l2 - 2*l3 + p2) - 6.74e-3*sin(p1 + p3 - 2*LPEJ - 2*G)
          + 6.66e-3*sin(l2 - 2*l3 + p4) + 4.45e-3*sin(l1 - p3)
          - 3.54e-3*sin(l1 - l2) - 3.17e-3*sin(2*(psi - LPEJ))
          + 2.65e-3*sin(l1 - p4) - 1.86e-3*sin(G)
          + 1.62e-3*sin(p2 - p3) + 1.58e-3*sin(4*(l1 - l2))
          - 1.55e-3*sin(l1 - l3) - 1.38e-3*sin(psi + w3 - 2*LPEJ - 2*G)
          - 1.15e-3*sin(2*(l1 - 2*l2 + w2)) + 8.9e-4*sin(p2 - p4)
          + 8.5e-4*sin(l1 + p3 - 2*LPEJ - 2*G) + 8.3e-4*sin(w2 - w3)
          + 5.3e-4*sin(psi - w2);
    sigma = pfmod(sigma, 360.0);
    sigma = degToRad(sigma);
    L = l1 + sigma;

    // Calculate periodic terms for the tangent of the latitude
    B = 6.393e-4*sin(L - w1) + 1.825e-4*sin(L - w2)
      + 3.29e-5*sin(L - w3) - 3.11e-5*sin(L - psi)
      + 9.3e-6*sin(L - w4) + 7.5e-6*sin(3*L - 4*l2 - 1.9927*sigma + w2)
      + 4.6e-6*sin(L + psi - 2*LPEJ - 2*G);
    B = atan(B);

    // Calculate the periodic terms for distance
    R = -4.1339e-3*cos(2*(l1 - l2)) - 3.87e-5*cos(l1 - p3)
      - 2.14e-5*cos(l1 - p4) + 1.7e-5*cos(l1 - l2)
      - 1.31e-5*cos(4*(l1 - l2)) + 1.06e-5*cos(l1 - l3)
      - 6.6e-6*cos(l1 + p3 - 2*LPEJ - 2*G);
    R = 5.90569 * JupiterRadius * (1 + R);

    T = (jd - 2433282.423) / 36525.0;
    P = 1.3966626*T + 3.088e-4*T*T;
    L += degToRad(P);

    L += JupAscendingNode;

    // Corrections for internal coordinate system
    B -= (PI/2);
    L += PI;

    return Point3d(cos(L) * sin(B) * R,
                   cos(B) * R,
                   -sin(L) * sin(B) * R);
    };

    double getPeriod() const
    {
        return 1.769138;
    };

    double getBoundingRadius() const
    {
        return 423329 * BoundingRadiusSlack;
    };
};

class EuropaOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
    // Computation will yield latitude(L), longitude(B) and distance(R) relative to Jupiter
    double t;
    double l1, l2, l3, l4;
    double p1, p2, p3, p4;
    double w1, w2, w3, w4;
    double gamma, phi, psi, G, Gp;
    double sigma, L, B, R;
    double T, P;

    // Epoch for Galilean satellites is 1976 Aug 10
    t = jd - 2443000.5;

    ComputeGalileanElements(t,
                            l1, l2, l3, l4,
                            p1, p2, p3, p4,
                            w1, w2, w3, w4,
                            gamma, phi, psi, G, Gp);

    // Calculate periodic terms for longitude
    sigma = 1.06476*sin(2*(l2 - l3)) + 0.04256*sin(l1 - 2*l2 + p3)
          + 0.03581*sin(l2 - p3) + 0.02395*sin(l1 - 2*l2 + p4)
          + 0.01984*sin(l2 - p4) - 0.01778*sin(phi)
          + 0.01654*sin(l2 - p2) + 0.01334*sin(l2 - 2*l3 + p2)
          + 0.01294*sin(p3 - p4) - 0.01142*sin(l2 - l3)
          - 0.01057*sin(G) - 7.75e-3*sin(2*(psi - LPEJ))
          + 5.24e-3*sin(2*(l1 - l2)) - 4.6e-3*sin(l1 - l3)
          + 3.16e-3*sin(psi - 2*G + w3 - 2*LPEJ) - 2.03e-3*sin(p1 + p3 - 2*LPEJ - 2*G)
          + 1.46e-3*sin(psi - w3) - 1.45e-3*sin(2*G)
          + 1.25e-3*sin(psi - w4) - 1.15e-3*sin(l1 - 2*l3 + p3)
          - 9.4e-4*sin(2*(l2 - w2)) + 8.6e-4*sin(2*(l1 - 2*l2 + w2))
          - 8.6e-4*sin(5*Gp - 2*G + 0.9115) - 7.8e-4*sin(l2 - l4)
          - 6.4e-4*sin(3*l3 - 7*l4 + 4*p4) + 6.4e-4*sin(p1 - p4)
          - 6.3e-4*sin(l1 - 2*l3 + p4) + 5.8e-4*sin(w3 - w4)
          + 5.6e-4*sin(2*(psi - LPEJ - G)) + 5.6e-4*sin(2*(l2 - l4))
          + 5.5e-4*sin(2*(l1 - l3)) + 5.2e-4*sin(3*l3 - 7*l4 + p3 +3*p4)
          - 4.3e-4*sin(l1 - p3) + 4.1e-4*sin(5*(l2 - l3))
          + 4.1e-4*sin(p4 - LPEJ) + 3.2e-4*sin(w2 - w3)
          + 3.2e-4*sin(2*(l3 - G - LPEJ));
    sigma = pfmod(sigma, 360.0);
    sigma = degToRad(sigma);
    L = l2 + sigma;

    // Calculate periodic terms for the tangent of the latitude
    B = 8.1004e-3*sin(L - w2) + 4.512e-4*sin(L - w3)
      - 3.284e-4*sin(L - psi) + 1.160e-4*sin(L - w4)
      + 2.72e-5*sin(l1 - 2*l3 + 1.0146*sigma + w2) - 1.44e-5*sin(L - w1)
      + 1.43e-5*sin(L + psi - 2*LPEJ - 2*G) + 3.5e-6*sin(L - psi + G)
      - 2.8e-6*sin(l1 - 2*l3 + 1.0146*sigma + w3);
    B = atan(B);

    // Calculate the periodic terms for distance
    R = 9.3848e-3*cos(l1 - l2) - 3.116e-4*cos(l2 - p3)
      - 1.744e-4*cos(l2 - p4) - 1.442e-4*cos(l2 - p2)
      + 5.53e-5*cos(l2 - l3) + 5.23e-5*cos(l1 - l3)
      - 2.9e-5*cos(2*(l1 - l2)) + 1.64e-5*cos(2*(l2 - w2))
      + 1.07e-5*cos(l1 - 2*l3 + p3) - 1.02e-5*cos(l2 - p1)
      - 9.1e-6*cos(2*(l1 - l3));
    R = 9.39657 * JupiterRadius * (1 + R);

    T = (jd - 2433282.423) / 36525.0;
    P = 1.3966626*T + 3.088e-4*T*T;
    L += degToRad(P);

    L += JupAscendingNode;

    // Corrections for internal coordinate system
    B -= (PI/2);
    L += PI;

    return Point3d(cos(L) * sin(B) * R,
                   cos(B) * R,
                   -sin(L) * sin(B) * R);
    };

    double getPeriod() const
    {
        return 3.5511810791;
    };

    double getBoundingRadius() const
    {
        return 678000 * BoundingRadiusSlack;
    };
};

class GanymedeOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
    //Computation will yield latitude(L), longitude(B) and distance(R) relative to Jupiter
    double t;
    double l1, l2, l3, l4;
    double p1, p2, p3, p4;
    double w1, w2, w3, w4;
    double gamma, phi, psi, G, Gp;
    double sigma, L, B, R;
    double T, P;

    //Epoch for Galilean satellites is 1976 Aug 10
    t = jd - 2443000.5;

    ComputeGalileanElements(t,
                            l1, l2, l3, l4,
                            p1, p2, p3, p4,
                            w1, w2, w3, w4,
                            gamma, phi, psi, G, Gp);

    //Calculate periodic terms for longitude
    sigma = 0.1649*sin(l3 - p3) + 0.09081*sin(l3 - p4)
          - 0.06907*sin(l2 - l3) + 0.03784*sin(p3 - p4)
          + 0.01846*sin(2*(l3 - l4)) - 0.01340*sin(G)
          - 0.01014*sin(2*(psi - LPEJ)) + 7.04e-3*sin(l2 - 2*l3 + p3)
          - 6.2e-3*sin(l2 - 2*l3 + p2) - 5.41e-3*sin(l3 - l4)
          + 3.81e-3*sin(l2 - 2*l3 + p4) + 2.35e-3*sin(psi - w3)
          + 1.98e-3*sin(psi - w4) + 1.76e-3*sin(phi)
          + 1.3e-3*sin(3*(l3 - l4)) + 1.25e-3*sin(l1 - l3)
          - 1.19e-3*sin(5*Gp - 2*G + 0.9115) + 1.09e-3*sin(l1 - l2)
          - 1.0e-3*sin(3*l3 - 7*l4 + 4*p4) + 9.1e-4*sin(w3 - w4)
          + 8.0e-4*sin(3*l3 - 7*l4 + p3 + 3*p4) - 7.5e-4*sin(2*l2 - 3*l3 + p3)
          + 7.2e-4*sin(p1 + p3 - 2*LPEJ - 2*G) + 6.9e-4*sin(p4 - LPEJ)
          - 5.8e-4*sin(2*l3 - 3*l4 + p4) - 5.7e-4*sin(l3 - 2*l4 + p4)
          + 5.6e-4*sin(l3 + p3 - 2*LPEJ - 2*G) - 5.2e-4*sin(l2 - 2*l3 + p1)
          - 5.0e-4*sin(p2 - p3) + 4.8e-4*sin(l3 - 2*l4 + p3)
          - 4.5e-4*sin(2*l2 - 3*l3 + p4) - 4.1e-4*sin(p2 - p4)
          - 3.8e-4*sin(2*G) - 3.7e-4*sin(p3 - p4 + w3 - w4)
          - 3.2e-4*sin(3*l3 - 7*l4 + 2*p3 + 2*p4) + 3.0e-4*sin(4*(l3 - l4))
          + 2.9e-4*sin(l3 + p4 - 2*LPEJ - 2*G) - 2.8e-4*sin(w3 + psi - 2*LPEJ - 2*G)
          + 2.6e-4*sin(l3 - LPEJ - G) + 2.4e-4*sin(l2 - 3*l3 + 2*l4)
          + 2.1e-4*sin(2*(l3 - LPEJ - G)) - 2.1e-4*sin(l3 - p2)
          + 1.7e-4*sin(l3 - p3);
    sigma = pfmod(sigma, 360.0);
    sigma = degToRad(sigma);
    L = l3 + sigma;

    //Calculate periodic terms for the tangent of the latitude
    B = 3.2402e-3*sin(L - w3) - 1.6911e-3*sin(L - psi)
      + 6.847e-4*sin(L - w4) - 2.797e-4*sin(L - w2)
      + 3.21e-5*sin(L + psi - 2*LPEJ - 2*G) + 5.1e-6*sin(L - psi + G)
      - 4.5e-6*sin(L - psi - G) - 4.5e-6*sin(L + psi - 2*LPEJ)
      + 3.7e-6*sin(L + psi - 2*LPEJ - 3*G) + 3.0e-6*sin(2*l2 - 3*L + 4.03*sigma + w2)
      - 2.1e-6*sin(2*l2 - 3*L + 4.03*sigma + w3);
    B = atan(B);

    //Calculate the periodic terms for distance
    R = -1.4388e-3*cos(l3 - p3) - 7.919e-4*cos(l3 - p4)
      + 6.342e-4*cos(l2 - l3) - 1.761e-4*cos(2*(l3 - l4))
      + 2.94e-5*cos(l3 - l4) - 1.56e-5*cos(3*(l3 - l4))
      + 1.56e-5*cos(l1 - l3) - 1.53e-5*cos(l1 - l2)
      + 7.0e-6*cos(2*l2 - 3*l3 + p3) - 5.1e-6*cos(l3 + p3 - 2*LPEJ - 2*G);
    R = 14.98832 * JupiterRadius * (1 + R);

    T = (jd - 2433282.423) / 36525.0;
    P = 1.3966626*T + 3.088e-4*T*T;
    L += degToRad(P);

    L += JupAscendingNode;

    //Corrections for internal coordinate system
    B -= (PI/2);
    L += PI;

    return Point3d(cos(L) * sin(B) * R,
                   cos(B) * R,
                   -sin(L) * sin(B) * R);
    };

    double getPeriod() const
    {
        return 7.154553;
    };

    double getBoundingRadius() const
    {
        return 1070000 * BoundingRadiusSlack;
    };
};

class CallistoOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
    //Computation will yield latitude(L), longitude(B) and distance(R) relative to Jupiter
    double t;
    double l1, l2, l3, l4;
    double p1, p2, p3, p4;
    double w1, w2, w3, w4;
    double gamma, phi, psi, G, Gp;
    double sigma, L, B, R;
    double T, P;

    //Epoch for Galilean satellites is 1976 Aug 10
    t = jd - 2443000.5;

    ComputeGalileanElements(t,
                            l1, l2, l3, l4,
                            p1, p2, p3, p4,
                            w1, w2, w3, w4,
                            gamma, phi, psi, G, Gp);

    //Calculate periodic terms for longitude
    sigma =
        0.84287*sin(l4 - p4)
        + 0.03431*sin(p4 - p3)
        - 0.03305*sin(2*(psi - LPEJ))
        - 0.03211*sin(G)
        - 0.01862*sin(l4 - p3)
        + 0.01186*sin(psi - w4)
        + 6.23e-3*sin(l4 + p4 - 2*G - 2*LPEJ)
        + 3.87e-3*sin(2*(l4 - p4))
        - 2.84e-3*sin(5*Gp - 2*G + 0.9115)
        - 2.34e-3*sin(2*(psi - p4))
        - 2.23e-3*sin(l3 - l4)
        - 2.08e-3*sin(l4 - LPEJ)
        + 1.78e-3*sin(psi + w4 - 2*p4)
        + 1.34e-3*sin(p4 - LPEJ)
        + 1.25e-3*sin(2*(l4 - G - LPEJ))
        - 1.17e-3*sin(2*G)
        - 1.12e-3*sin(2*(l3 - l4))
        + 1.07e-3*sin(3*l3 - 7*l4 + 4*p4)
        + 1.02e-3*sin(l4 - G - LPEJ)
        + 9.6e-4*sin(2*l4 - psi - w4)
        + 8.7e-4*sin(2*(psi - w4))
        - 8.5e-4*sin(3*l3 - 7*l4 + p3 + 3*p4)
        + 8.5e-4*sin(l3 - 2*l4 + p4)
        - 8.1e-4*sin(2*(l4 - psi))
        + 7.1e-4*sin(l4 + p4 - 2*LPEJ - 3*G)
        + 6.1e-4*sin(l1 - l4)
        - 5.6e-4*sin(psi - w3)
        - 5.4e-4*sin(l3 - 2*l4 + p3)
        + 5.1e-4*sin(l2 - l4)
        + 4.2e-4*sin(2*(psi - G - LPEJ))
        + 3.9e-4*sin(2*(p4 - w4))
        + 3.6e-4*sin(psi + LPEJ - p4 - w4)
        + 3.5e-4*sin(2*Gp - G + 3.2877)
        - 3.5e-4*sin(l4 - p4 + 2*LPEJ - 2*psi)
        - 3.2e-4*sin(l4 + p4 - 2*LPEJ - G)
        + 3.0e-4*sin(2*Gp - 2*G + 2.6032)
        + 2.9e-4*sin(3*l3 - 7*l4 + 2*p3 + 2*p4)
        + 2.8e-4*sin(l4 - p4 + 2*psi - 2*LPEJ)
        - 2.8e-4*sin(2*(l4 - w4))
        - 2.7e-4*sin(p3 - p4 + w3 - w4)
        - 2.6e-4*sin(5*Gp - 3*G + 3.2877)
        + 2.5e-4*sin(w4 - w3)
        - 2.5e-4*sin(l2 - 3*l3 + 2*l4)
        - 2.3e-4*sin(3*(l3 - l4))
        + 2.1e-4*sin(2*l4 - 2*LPEJ - 3*G)
        - 2.1e-4*sin(2*l3 - 3*l4 + p4)
        + 1.9e-4*sin(l4 - p4 - G)
        - 1.9e-4*sin(2*l4 - p3 - p4)
        - 1.8e-4*sin(l4 - p4 + G)
        - 1.6e-4*sin(l4 + p3 - 2*LPEJ - 2*G);
    sigma = pfmod(sigma, 360.0);
    sigma = degToRad(sigma);
    L = l4 + sigma;

    //Calculate periodic terms for the tangent of the latitude
    B =
        - 7.6579e-3 * sin(L - psi)
        + 4.4134e-3 * sin(L - w4)
        - 5.112e-4  * sin(L - w3)
        + 7.73e-5   * sin(L + psi - 2*LPEJ - 2*G)
        + 1.04e-5   * sin(L - psi + G)
        - 1.02e-5   * sin(L - psi - G)
        + 8.8e-6    * sin(L + psi - 2*LPEJ - 3*G)
        - 3.8e-6    * sin(L + psi - 2*LPEJ - G);
    B = atan(B);

    //Calculate the periodic terms for distance
    R =
        - 7.3546e-3 * cos(l4 - p4)
        + 1.621e-4  * cos(l4 - p3)
        + 9.74e-5   * cos(l3 - l4)
        - 5.43e-5   * cos(l4 + p4 - 2*LPEJ - 2*G)
        - 2.71e-5   * cos(2*(l4 - p4))
        + 1.82e-5   * cos(l4 - LPEJ)
        + 1.77e-5   * cos(2*(l3 - l4))
        - 1.67e-5   * cos(2*l4 - psi - w4)
        + 1.67e-5   * cos(psi - w4)
        - 1.55e-5   * cos(2*(l4 - LPEJ - G))
        + 1.42e-5   * cos(2*(l4 - psi))
        + 1.05e-5   * cos(l1 - l4)
        + 9.2e-6    * cos(l2 - l4)
        - 8.9e-6    * cos(l4 - LPEJ -G)
        - 6.2e-6    * cos(l4 + p4 - 2*LPEJ - 3*G)
        + 4.8e-6    * cos(2*(l4 - w4));
    
    R = 26.36273 * JupiterRadius * (1 + R);

    T = (jd - 2433282.423) / 36525.0;
    P = 1.3966626*T + 3.088e-4*T*T;
    L += degToRad(P);

    L += JupAscendingNode;

    //Corrections for internal coordinate system
    B -= (PI/2);
    L += PI;

    return Point3d(cos(L) * sin(B) * R,
                   cos(B) * R,
                   -sin(L) * sin(B) * R);
    };

    double getPeriod() const
    {
        return 16.689018;
    };

    double getBoundingRadius() const
    {
        return 1890000 * BoundingRadiusSlack;
    };
};


static const double SatAscendingNode = 168.8112;
static const double SatTilt = 28.0817;
// static const double SatAscendingNode = 169.530;
// static const double SatTilt = 28.049;

// Calculations for the orbits of Mimas, Enceladus, Tethys, Dione, Rhea,
// Titan, Hyperion, and Iapetus are from Jean Meeus's Astronomical Algorithms,
// and were originally derived by Gerard Dourneau.

void ComputeSaturnianElements(double t,
                              double& t1, double& t2, double& t3,
                              double& t4, double& t5, double& t6,
                              double& t7, double& t8, double& t9,
                              double& t10, double& t11,
                              double& W0, double& W1, double& W2,
                              double& W3, double& W4, double& W5,
                              double& W6, double& W7, double& W8)
{
    t1 = t - 2411093.0;
    t2 = t1 / 365.25;
    t3 = (t - 2433282.423) / 365.25 + 1950.0;
    t4 = t - 2411368.0;
    t5 = t4 / 365.25;
    t6 = t - 2415020.0;
    t7 = t6 / 36525;
    t8 = t6 / 365.25;
    t9 = (t - 2442000.5) / 365.25;
    t10 = t - 2409786.0;
    t11 = t10 / 36525;

    W0 = 5.095 * (t3 - 1866.39);
    W1 = 74.4 + 32.39 * t2;
    W2 = 134.3 + 92.62 * t2;
    W3 = 42.0 - 0.5118 * t5;
    W4 = 276.59 + 0.5118 * t5;
    W5 = 267.2635 + 1222.1136 * t7;
    W6 = 175.4762 + 1221.5515 * t7;
    W7 = 2.4891 + 0.002435 * t7;
    W8 = 113.35 - 0.2597 * t7;
}


static Point3d SaturnMoonPosition(double lam, double gam, double Om, double r)
{
    double u = lam - Om;
    double w = Om - SatAscendingNode;

    u = degToRad(u);
    w = degToRad(w);
    gam = -degToRad(gam);
    r = r * SaturnRadius;

    // Corrections for Celestia's coordinate system
    u = -u;
    w = -w;

    double x = r * (cos(u) * cos(w) - sin(u) * sin(w) * cos(gam));
    double y = r * sin(u) * sin(gam);
    double z = r * (sin(u) * cos(w) * cos(gam) + cos(u) * sin(w));

    return Point3d(x, y, z);
}


static void OuterSaturnMoonParams(double a, double e, double i,
                                  double Om_, double M, double lam_,
                                  double& lam, double& gam,
                                  double& r, double& w)
{
    double s1 = sinD(SatTilt);
    double c1 = cosD(SatTilt);
    double e_2 = e * e;
    double e_3 = e_2 * e;
    double e_4 = e_3 * e;
    double e_5 = e_4 * e;
    double C = (2 * e - 0.25 * e_3 + 0.0520833333 * e_5) * sinD(M) +
        (1.25 * e_2 - 0.458333333 * e_4) * sinD(2 * M) +
        (1.083333333 * e_3 - 0.671875 * e_5) * sinD(3 * M) +
        1.072917 * e_4 * sinD(4 * M) + 1.142708 * e_5 * sinD(5 * M);
    double g = Om_ - SatAscendingNode;
    double a1 = sinD(i) * sinD(g);
    double a2 = c1 * sinD(i) * cosD(g) - s1 * cosD(i);
    double u = radToDeg(atan2(a1, a2));
    double h = c1 * sinD(i) - s1 * cosD(i) * cosD(g);
    double psi = radToDeg(atan2(s1 * sinD(g), h));

    C = radToDeg(C);
    lam = lam_ + C + u - g - psi;
    gam = radToDeg(asin(sqrt(square(a1) + square(a2))));
    r = a * (1 - e * e) / (1 + e * cosD(M + C));
    w = SatAscendingNode + u;
}


class MimasOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
        // Computation will yield latitude(L), longitude(B) and distance(R)
        // relative to Saturn.
        double t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11;
        double W0, W1, W2, W3, W4, W5, W6, W7, W8;

        ComputeSaturnianElements(jd,
                                 t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11,
                                 W0, W1, W2, W3, W4, W5, W6, W7, W8);

        double L = 127.64 + 381.994497 * t1 - 43.57 * sinD(W0) -
            0.720 * sinD( 3 * W0) - 0.02144 * sinD(5 * W0);
        double p = 106.1 + 365.549 * t2;
        double M = L - p;
        double C = 2.18287 * sinD(M) + 0.025988 * sinD(2 * M) +
            0.00043 * sinD(3 * M);
        double lam = L + C;
        double r = 3.06879 / (1 + 0.01905 * cosD(M + C));
        double gam = 1.563;
        double Om = 54.5 - 365.072 * t2;

        return SaturnMoonPosition(lam, gam, Om, r);
    };

    double getPeriod() const
    {
        return 0.9424218;
    };

    double getBoundingRadius() const
    {
        return 189000 * BoundingRadiusSlack;
    };
};


class EnceladusOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
        // Computation will yield latitude(L), longitude(B) and distance(R)
        // relative to Saturn.
        double t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11;
        double W0, W1, W2, W3, W4, W5, W6, W7, W8;

        ComputeSaturnianElements(jd,
                                 t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11,
                                 W0, W1, W2, W3, W4, W5, W6, W7, W8);

        double L = 200.317 + 262.7319002 * t1 + 0.25667 * sinD(W1) +
            0.20883 * sinD(W2);
        double p = 309.107 + 123.44121 * t2;
        double M = L - p;
        double C = 0.55577 * sinD(M) + 0.00168 * sinD(2 * M);
        double lam = L + C;
        double r = 3.94118 / (1 + 0.00485 * cosD(M + C));
        double gam = 0.0262;
        double Om = 348 - 151.95 * t2;

        return SaturnMoonPosition(lam, gam, Om, r);
    };

    double getPeriod() const
    {
        return 1.370218;
    };

    double getBoundingRadius() const
    {
        return 239000 * BoundingRadiusSlack;
    };
};


class TethysOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
        // Computation will yield latitude(L), longitude(B) and distance(R)
        // relative to Saturn.
        double t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11;
        double W0, W1, W2, W3, W4, W5, W6, W7, W8;

        ComputeSaturnianElements(jd,
                                 t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11,
                                 W0, W1, W2, W3, W4, W5, W6, W7, W8);

        double lam = 285.306 + 190.69791226 * t1 + 2.063 * sinD(W0) +
            0.03409 * sinD(3 * W0) + 0.001015 * sinD(5 * W0);
        double r = 4.880998;
        double gam = 1.0976;
        double Om = 111.33 - 72.2441 * t2;

        return SaturnMoonPosition(lam, gam, Om, r);
    };

    double getPeriod() const
    {
        return 1.887802;
    };

    double getBoundingRadius() const
    {
        return 295000 * BoundingRadiusSlack;
    };
};


class DioneOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
        // Computation will yield latitude(L), longitude(B) and distance(R)
        // relative to Saturn.
        double t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11;
        double W0, W1, W2, W3, W4, W5, W6, W7, W8;

        ComputeSaturnianElements(jd,
                                 t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11,
                                 W0, W1, W2, W3, W4, W5, W6, W7, W8);

        double L = 254.712 + 131.53493193 * t1 - 0.0215 * sinD(W1) -
            0.01733 * sinD(W2);
        double p = 174.8 + 30.820 * t2;
        double M = L - p;
        double C = 0.24717 * sinD(M) + 0.00033 * sinD(2 * M);
        double lam = L + C;
        double r = 6.24871 / (1 + 0.002157 * cosD(M + C));
        double gam = 0.0139;
        double Om = 232 - 30.27 * t2;
        // cout << "Dione: " << pfmod(lam, 360.0) << ',' << gam << ',' << pfmod(Om, 360.0) << ',' << r << '\n';

        return SaturnMoonPosition(lam, gam, Om, r);
    };

    double getPeriod() const
    {
        return 2.736915;
    };

    double getBoundingRadius() const
    {
        return 378000 * BoundingRadiusSlack;
    };
};


class RheaOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
        // Computation will yield latitude(L), longitude(B) and distance(R)
        // relative to Saturn.
        double t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11;
        double W0, W1, W2, W3, W4, W5, W6, W7, W8;

        ComputeSaturnianElements(jd,
                                 t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11,
                                 W0, W1, W2, W3, W4, W5, W6, W7, W8);
        double e1 = 0.05589 - 0.000346 * t7;

        double p_ = 342.7 + 10.057 * t2;
        double a1 = 0.000265 * sinD(p_) + 0.01 * sinD(W4);
        double a2 = 0.000265 * cosD(p_) + 0.01 * cosD(W4);
        double e = sqrt(square(a1) + square(a2));
        double p = radToDeg(atan2(a1, a2));
        double N = 345 - 10.057 * t2;
        double lam_ = 359.244 + 79.69004720 * t1 + 0.086754 * sinD(N);
        double i = 28.0362 + 0.346898 * cosD(N) + 0.01930 * cosD(W3);
        double Om = 168.8034 + 0.736936 * sinD(N) + 0.041 * sinD(W3);
        double a = 8.725924;

        double lam, gam, r, w;
        OuterSaturnMoonParams(a, e, i, Om, lam_ - p, lam_,
                              lam, gam, r, w);
        // cout << "Rhea (intermediate): " << e << ',' << pfmod(lam_, 360.0) << ',' << pfmod(i, 360.0) << ',' << pfmod(Om, 360.0) << '\n';
        // cout << "Rhea: " << pfmod(lam, 360.0) << ',' << gam << ',' << pfmod(w, 360.0) << ',' << r << '\n';

        return SaturnMoonPosition(lam, gam, w, r);
    };

    double getPeriod() const
    {
        return 4.517500;
    };

    double getBoundingRadius() const
    {
        return 528000 * BoundingRadiusSlack;
    };
};


class TitanOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
        // Computation will yield latitude(L), longitude(B) and distance(R)
        // relative to Saturn.
        double t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11;
        double W0, W1, W2, W3, W4, W5, W6, W7, W8;

        ComputeSaturnianElements(jd,
                                 t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11,
                                 W0, W1, W2, W3, W4, W5, W6, W7, W8);
        double e1 = 0.05589 - 0.000346 * t7;

        double L = 261.1582 + 22.57697855 * t4 + 0.074025 * sinD(W3);
        double i_ = 27.45141 + 0.295999 * cosD(W3);
        double Om_ = 168.66925 + 0.628808 * sinD(W3);
        double a1 = sinD(W7) * sinD(Om_ - W8);
        double a2 = cosD(W7) * sinD(i_) - sinD(W7) * cosD(i_) * cosD(Om_ - W8);
        double g0 = 102.8623;
        double psi = radToDeg(atan2(a1, a2));
        double s = sqrt(square(a1) + square(a2));
        double g = W4 - Om_ - psi;

        // Three successive approximations will always be enough
        double om;
        for (int n = 0; n < 3; n++)
        {
            om = W4 + 0.37515 * (sinD(2 * g) - sinD(2 * g0));
            g = om - Om_ - psi;
        }

        double e_ = 0.029092 + 0.00019048 * (cosD(2 * g) - cosD(2 * g0));
        double q = 2 * (W5 - om);
        double b1 = sinD(i_) * sinD(Om_ - W8);
        double b2 = cosD(W7) * sinD(i_) * cosD(Om_ - W8) - sinD(W7) * cosD(i_);
        double theta = radToDeg(atan2(b1, b2)) + W8;
        double e = e_ + 0.002778797 * e_ * cosD(q);
        double p = om + 0.159215 * sinD(q);
        double u = 2 * W5 - 2 * theta + psi;
        double h = 0.9375 * square(e_) * sinD(q) + 0.1875 * square(s) * sinD(2 * (W5 - theta));
        double lam_ = L - 0.254744 * (e1 * sinD(W6) + 0.75 * square(e1) * sinD(2 * W6) + h);
        double i = i_ + 0.031843 * s * cosD(u);
        double Om = Om_ + (0.031843 * s * sinD(u)) / sinD(i_);
        double a = 20.216193;

        double lam, gam, r, w;
        OuterSaturnMoonParams(a, e, i, Om, lam_ - p, lam_,
                              lam, gam, r, w);

        return SaturnMoonPosition(lam, gam, w, r);
    };

    double getPeriod() const
    {
        return 15.94544758;
    };

    double getBoundingRadius() const
    {
        return 1260000 * BoundingRadiusSlack;
    };
};


class HyperionOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
        // Computation will yield latitude(L), longitude(B) and distance(R)
        // relative to Saturn.
        double t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11;
        double W0, W1, W2, W3, W4, W5, W6, W7, W8;

        ComputeSaturnianElements(jd,
                                 t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11,
                                 W0, W1, W2, W3, W4, W5, W6, W7, W8);
        double eta = 92.39 + 0.5621071 * t6;
        double zeta = 148.19 - 19.18 * t8;
        double theta = 184.8 - 35.41 * t9;
        double theta_ = theta - 7.5;
        double as = 176 + 12.22 * t8;
        double bs = 8 + 24.44 * t8;
        double cs = bs + 5;
        double om = 68.898 - 18.67088 * t8;
        double phi = 2 * (om - W5);
        double chi = 94.9 - 2.292 * t8;
        double a = 24.50601 -
            0.08686 * cosD(eta) -
            0.00166 * cosD(zeta + eta) +
            0.00175 * cosD(zeta - eta);
        double e = 0.103458 -
            0.004099 * cosD(eta) -
            0.000167 * cosD(zeta + eta) +
            0.000235 * cosD(zeta - eta) +
            0.02303 * cosD(zeta) -
            0.00212 * cosD(2 * zeta) +
            0.000151 * cosD(3 * zeta) +
            0.00013 * sinD(phi);
        double p = om +
            0.15648 * sinD(chi) -
            0.4457 * sinD(eta) -
            0.2657 * sinD(zeta + eta) -
            0.3573 * sinD(zeta - eta) -
            12.872 * sinD(zeta) +
            1.668 * sinD(2 * zeta) -
            0.2419 * sinD(3 * zeta) - 
            0.07 * sinD(phi);
        double lam_ = 177.047 +
            16.91993829 * t6 +
            0.15648 * sinD(chi) +
            9.142 * sinD(eta) +
            0.007 * sinD(2 * eta) -
            0.014 * sinD(3 * eta) +
            0.2275 * sinD(zeta + eta) +
            0.2112 * sinD(zeta - eta) -
            0.26 * sinD(zeta) -
            0.0098 * sinD(2 * zeta) -
            0.013 * sinD(as) +
            0.017 * sinD(bs) -
            0.0303 * sinD(phi);
        double i = 27.3347 + 0.643486 * cosD(chi) + 0.315 * cosD(W3) +
            0.018 * cosD(theta) - 0.018 * cosD(cs);
        double Om = 168.6812 + 1.40136 * cosD(chi) + 0.68599 * sinD(W3) -
            0.0392 * sinD(cs) + 0.0366 * sinD(theta_);

        double lam, gam, r, w;
        OuterSaturnMoonParams(a, e, i, Om, lam_ - p, lam_,
                              lam, gam, r, w);

        return SaturnMoonPosition(lam, gam, w, r);
    };

    double getPeriod() const
    {
        return 21.276609;
    };

    double getBoundingRadius() const
    {
        return 1640000 * BoundingRadiusSlack;
    };
};


class IapetusOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
        // Computation will yield latitude(L), longitude(B) and distance(R)
        // relative to Saturn.
        double t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11;
        double W0, W1, W2, W3, W4, W5, W6, W7, W8;

        ComputeSaturnianElements(jd,
                                 t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11,
                                 W0, W1, W2, W3, W4, W5, W6, W7, W8);
        double L = 261.1582 + 22.57697855 * t4;
        double om_ = 91.796 + 0.562 * t7;
        double psi = 4.367 - 0.195 * t7;
        double theta = 146.819 - 3.198 * t7;
        double phi = 60.470 + 1.521 * t7;
        double Phi = 205.055 - 2.091 * t7;
        double e_ = 0.028298 + 0.001156 * t11;
        double om0 = 352.91 + 11.71 * t11;
        double mu = 76.3852 + 4.53795125 * t10;
        double i_ = 18.4602 - 0.9518 * t11 - 0.072 * square(t11) +
            0.0054 * cube(t11);
        double Om_ = 143.198 - 3.919 * t11 + 0.116 * square(t11) +
            0.008 * cube(t11);
        double l = mu - om0;
        double g = om0 - Om_ - psi;
        double g1 = om0 - Om_ - phi;
        double ls = W5 - om_;
        double gs = om_ - theta;
        double lT = L - W4;
        double gT = W4 - Phi;
        double u1 = 2 * (l + g - (ls + gs));
        double u2 = l + g1 - (lT + gT);
        double u3 = l + 2 * (g - (ls + gs));
        double u4 = lT + gT - g1;
        double u5 = 2 * (ls + gs);

        double a = 58.935028 + 0.004638 * cosD(u1) + 0.058222 * cosD(u2);
        double e = e_ -
            0.0014097 * cosD(g1 - gT) +
            0.0003733 * cosD(u5 - 2 * g) +
            0.0001180 * cosD(u3) +
            0.0002408 * cosD(l) +
            0.0002849 * cosD(l + u2) +
            0.0006190 * cosD(u4);
        double W = 0.08077 * sinD(g1 - gT) +
            0.02139 * sinD(u5 - 2 * g) -
            0.00676 * sinD(u3) +
            0.01380 * sinD(l) +
            0.01632 * sinD(l + u2) +
            0.03547 * sinD(u4);
        double p = om0 + W / e_;
        double lam_ = mu - 
            0.04299 * sinD(u2) -
            0.00789 * sinD(u1) -
            0.06312 * sinD(ls) -
            0.00295 * sinD(2 * ls) -
            0.02231 * sinD(u5) +
            0.00650 * sinD(u5 + psi);
        double sum = l + g1 + lT + gT + phi;
        double i = i_ +
            0.04204 * cosD(u5 + psi) +
            0.00235 * cosD(sum) +
            0.00360 * cosD(u2 + phi);
        double w_ = 0.04204 * sinD(u5 + psi) +
            0.00235 * sinD(sum) +
            0.00358 * sinD(u2 + phi);
        double Om = Om_ + w_ / sinD(i_);

        double lam, gam, r, w;
        OuterSaturnMoonParams(a, e, i, Om, lam_ - p, lam_,
                              lam, gam, r, w);

        return SaturnMoonPosition(lam, gam, w, r);
    };

    double getPeriod() const
    {
        return 79.330183;
    };

    double getBoundingRadius() const
    {
        return 3660000 * BoundingRadiusSlack;
    };
};


class PhoebeOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
        double t = jd - 2433282.5;
        double T = t / 365.25;

        double a = astro::AUtoKilometers(0.0865752f) / SaturnRadius;
        double lam_ = 277.872 - 0.6541068 * t - 90;
        double e = 0.16326;
        double pi = 280.165 - 0.19586 * T;
        double i = 173.949 - 0.020 * T;
        double Om = 245.998 - 0.41353 * T;

        double lam, gam, r, w;
        OuterSaturnMoonParams(a, e, i, Om, lam_ - pi, lam_,
                              lam, gam, r, w);

        return SaturnMoonPosition(lam, gam, w, r);
    };

    double getPeriod() const
    {
        return 548.2122790;
    };

    double getBoundingRadius() const
    {
        return 15100000 * BoundingRadiusSlack;
    };
};


class UranianSatelliteOrbit : public CachingOrbit
{
private:
    double a;
    double n;
    double L0;
    double L1;
    double *L_k, *L_theta, *L_phi;
    int LTerms;
    double *z_k, *z_theta, *z_phi;
    int zTerms;
    double *zeta_k, *zeta_theta, *zeta_phi;
    int zetaTerms;

public:
    UranianSatelliteOrbit(double _a,
                          double _n,
                          double _L0, double _L1,
                          int _LTerms, int _zTerms, int _zetaTerms,
                          double* _L_k, double* _L_theta, double* _L_phi,
                          double* _z_k, double* _z_theta, double* _z_phi,
                          double* _zeta_k, double* _zeta_theta,
                          double* _zeta_phi) :
        a(_a), n(_n), L0(_L0), L1(_L1),
        LTerms(_LTerms), zTerms(_zTerms), zetaTerms(_zetaTerms),
        L_k(_L_k), L_theta(_L_theta), L_phi(_L_phi),
        z_k(_z_k), z_theta(_z_theta), z_phi(_z_phi),
        zeta_k(_zeta_k), zeta_theta(_zeta_theta), zeta_phi(_zeta_phi)
    {
    };

public:
    double getPeriod() const
    {
        return 2 * PI / n;
    }

    double getBoundingRadius() const
    {
        // Not quite correct, but should work since e is pretty low
        // for most of the Uranian moons.
        return a * BoundingRadiusSlack;
    }

    Point3d computePosition(double jd) const
    {
        double t = jd - 2444239.5;
        int i;

        double L = L0 + L1 * t;
        for (i = 0; i < LTerms; i++)
            L += L_k[i] * sin(L_theta[i] * t + L_phi[i]);

        double a0 = 0.0;
        double a1 = 0.0;
        for (i = 0; i < zTerms; i++)
        {
            double w = z_theta[i] * t + z_phi[i];
            a0 += z_k[i] * cos(w);
            a1 += z_k[i] * sin(w);
        }

        double b0 = 0.0;
        double b1 = 0.0;
        for (i = 0; i < zetaTerms; i++)
        {
            double w = zeta_theta[i] * t + zeta_phi[i];
            b0 += zeta_k[i] * cos(w);
            b1 += zeta_k[i] * sin(w);
        }

        double e = sqrt(square(a0) + square(a1));
        double p = atan2(a1, a0);
        double gamma = 2.0 * asin(sqrt(square(b0) + square(b1)));
        double theta = atan2(b1, b0);

        L += degToRad(174.99);

        // Now that we have all the orbital elements, compute the position
        double M = L - p;

        // Iterate a few times to compute the eccentric anomaly from the
        // mean anomaly.
        double ecc = M;
        for (i = 0; i < 4; i++)
            ecc = M + e * sin(ecc);

        double x = a * (cos(ecc) - e);
        double z = a * sqrt(1 - square(e)) * -sin(ecc);

        Mat3d R = (Mat3d::yrotation(theta) *
                   Mat3d::xrotation(gamma) *
                   Mat3d::yrotation(p - theta));
        return R * Point3d(x, 0, z);
    }
};


static double uran_n[5] =
{ 4.44352267, 2.49254257, 1.51595490, 0.72166316, 0.46658054 };
static double uran_a[5] =
{ 129800, 191200, 266000, 435800, 583600 };
static double uran_L0[5] =
{ -0.23805158, 3.09804641, 2.28540169, 0.85635879, -0.91559180 };
static double uran_L1[5] =
{ 4.44519055, 2.49295252, 1.51614811, 0.72171851, 0.46669212 };
static double uran_L_k[5][3] = {
{  0.02547217, -0.00308831, -3.181e-4 },
{ -1.86050e-3, 2.1999e-4, 0 },
{ 6.6057e-4, 0, 0 },
{ 0, 0, 0 },
{ 0, 0, 0 }
};
static double uran_L_theta[5][3] = {
{ -2.18167e-4, -4.36336e-4, -6.54502e-4 },
{ -2.18167e-4, -4.36336e-4, 0 },
{ -2.18167e-4, 0, 0 },
{ 0, 0, 0 },
{ 0, 0, 0 }
};
static double uran_L_phi[5][3] = {
{ 1.32, 2.64, 3.97 },
{ 1.32, 2.64, 0 },
{ 1.32, 0, 0 },
{ 0, 0, 0 },
{ 0, 0, 0 },   
};
static double uran_z_k[5][5] = {
{ 1.31238e-3, -1.2331e-4, -1.9410e-4, 0, 0 },
{ 1.18763e-3, 8.6159e-4, 0, 0, 0 },
{ -2.2795e-4, 3.90496e-3, 3.0917e-4, 2.2192e-4, 5.4923e-4 },
{ 9.3281e-4, 1.12089e-3, 7.9343e-4, 0, 0 },
{ -7.5868e-4, 1.39734e-3, -9.8726e-4, 0, 0 }
};
static double uran_z_theta[5][5] = {
{ 1.5273e-4, 0.08606, 0.709, 0, 0 },
{ 4.727824e-5, 2.179316e-5, 0, 0, 0 },
{ 4.727824e-5, 2.179132e-5, 1.580524e-5, 2.9363068e-6, -0.01157 },
{ 1.580524e-5, 2.9363068e-6, -6.9008e-3, 0, 0 },
{ 1.580524e-5, 2.9363068e-6, -6.9008e-3, 0, 0 }
};
static double uran_z_phi[5][5] = {
{ 0.61, 0.15, 6.04, 0, 0 },
{ 2.41, 2.07, 0, 0, 0 },
{ 2.41, 2.07, 0.74, 0.43, 5.71 },
{ 0.74, 0.43, 1.82, 0, 0 },
{ 0.74, 0.43, 1.82, 0, 0 }
};
static double uran_zeta_k[5][2] = {
{ 0.03787171, 0 },
{ 3.5825e-4, 2.9008e-4 },
{ 1.11336e-3, 3.5014e-4 },
{ 6.8572e-4, 3.7832e-4 },
{ -5.9633e-4, 4.5169e-4 }
};
static double uran_zeta_theta[5][2] = {
{ -1.54449e-4, 0 },
{ -4.782474e-5, -2.156628e-5 },
{ -2.156628e-5, -1.401373e-5 },
{ -1.401373e-5, -1.9713918e-6 },
{ -1.401373e-5, -1.9713918e-6 }
};
static double uran_zeta_phi[5][2] = {
{ 5.70, 0 },
{ 0.40, 0.59 },
{ 0.59, 1.75 },
{ 1.75, 4.21 },
{ 1.75, 4.21 },
};

static UranianSatelliteOrbit* CreateUranianSatelliteOrbit(int n)
{
    assert(n >= 1 && n <= 5);
    n--;

    return new UranianSatelliteOrbit(uran_a[n], uran_n[n],
                                     uran_L0[n], uran_L1[n],
                                     3, 5, 2,
                                     uran_L_k[n], uran_L_theta[n],
                                     uran_L_phi[n], uran_z_k[n],
                                     uran_z_theta[n], uran_z_phi[n],
                                     uran_zeta_k[n], uran_zeta_theta[n],
                                     uran_zeta_phi[n]);
};


class TritonOrbit : public CachingOrbit
{
    Point3d computePosition(double jd) const
    {
        // Pole of the fixed reference plane
        //double t = jd - 2433282.5;
        double epoch = 2433282.5;
        double t = jd - epoch;
        double T = (jd - 2451545.0) / 36525;
        double N = degToRad(359.28 + 54.308 * T);
        double refplane_RA  = degToRad(298.72 +
                                       2.58 * sin(N) - 0.04 * sin(2 * N));
        double refplane_Dec = degToRad(42.63 -
                                       1.90 * cos(N) - 0.01 * cos(2 * N));
        double nepPole_RA = degToRad(299.36);
        double nepPole_Dec = degToRad(43.46);

        double a = 354800;
        double e = 0.0000;
        double n = degToRad(61.2588532);
        double L0 = degToRad(200.913);
        double L  = L0 + n * t;
        double gamma = degToRad(158.996);
        double theta = degToRad(151.401 + 0.57806 * t / 365.25);

        double E = L - theta;
        Point3d p(a * cos(E), 0.0, a * -sin(E));

        // Orientation of the orbital plane with respect to the Laplacian plane
        Mat3d Rorbit     = (Mat3d::yrotation(theta) *
                            Mat3d::xrotation(gamma));

        // Rotate to the Earth's equatorial plane
        double Nr = degToRad(refplane_RA);
        double Jr = degToRad(90 - refplane_Dec);
        Mat3d Rrefplane = (Mat3d::yrotation( Nr) *
                           Mat3d::xrotation( Jr) *
                           Mat3d::yrotation(-Nr));

        // Rotate to the Neptunian equatorial plane
        double Nn = degToRad(nepPole_RA);
        double Jn = degToRad(90 - nepPole_Dec);
        Mat3d RNept_eq   = (Mat3d::yrotation( Nn) *
                            Mat3d::xrotation(-Jn) *
                            Mat3d::yrotation(-Nn));

        return RNept_eq * (Rrefplane * (Rorbit * p));
    }
    
    double getPeriod() const
    {
        return 5.877;
    }

    double getBoundingRadius() const
    {
        return 354800 * BoundingRadiusSlack;
    }
};


class JPLEphOrbit : public CachingOrbit
{
public:
    JPLEphOrbit(const JPLEphemeris& e, JPLEphemItem item,
                double _period, double _boundingRadius) :
        ephem(e),
        body(item),
        period(_period),
        boundingRadius(_boundingRadius)
    {
    };

    double getPeriod() const
    {
        return period;
    }

    double getBoundingRadius() const
    {
        return boundingRadius;
    }

    Point3d computePosition(double jd) const
    {
        Point3d baryPos = ephem.getPlanetPosition(body, jd);

        // Convert from solar system barycentric position to heliocentric
        // position.
        Point3d sunPos = ephem.getPlanetPosition(JPLEph_Sun, jd);
        Point3d heliocentricPos = Point3d(0.0, 0.0, 0.0) + (baryPos - sunPos);

        // Rotate from the J2000 mean equator to the ecliptic
        // 23 deg 26' 21".448
        double J2000_equator = degToRad(23.4392911);
        heliocentricPos = heliocentricPos * Mat3d::xrotation(J2000_equator);

        return Point3d(heliocentricPos.x,
                       heliocentricPos.z,
                       -heliocentricPos.y);
    }

private:
    const JPLEphemeris& ephem;
    JPLEphemItem body;
    double period;
    double boundingRadius;
};


static Orbit* CreateJPLEphOrbit(JPLEphemItem body,
                                double period,
                                double boundingRadius)
{
    Orbit* o = new JPLEphOrbit(*jpleph, body, period, boundingRadius);
    return new MixedOrbit(o,
                          jpleph->getStartDate(),
                          jpleph->getEndDate(),
                          astro::SolarMass);
}


static double yearToJD(int year)
{
    return (double) astro::Date(year, 1, 1);
}


Orbit* GetCustomOrbit(const string& name)
{
    // Attempt to load JPL ephemeris data if we haven't tried already
    if (!jplephInitialized)
    {
        jplephInitialized = true;
        ifstream in("data/jpleph.dat", ios::in | ios::binary);
        if (in.good())
            jpleph = JPLEphemeris::load(in);
        if (jpleph != NULL)
        {
            clog << "Loaded DE" << jpleph->getDENumber() <<
                " ephemeris. Valid from JD" <<
                setprecision(8) <<
                jpleph->getStartDate() << " to JD" <<
                jpleph->getEndDate() << '\n';
        }
    }

    if (name == "mercury")
        return new MixedOrbit(new MercuryOrbit(), yearToJD(-4000), yearToJD(4000), astro::SolarMass);
    if (name == "venus")
        return new MixedOrbit(new VenusOrbit(), yearToJD(-4000), yearToJD(4000), astro::SolarMass);
    if (name == "earth")
        return new MixedOrbit(new EarthOrbit(), yearToJD(-4000), yearToJD(4000), astro::SolarMass);
    if (name == "moon")
        return new MixedOrbit(new LunarOrbit(), yearToJD(-2000), yearToJD(4000), astro::EarthMass + astro::LunarMass);
    if (name == "mars")
        return new MixedOrbit(new MarsOrbit(), yearToJD(-4000), yearToJD(4000), astro::SolarMass);
    if (name == "jupiter")
        return new MixedOrbit(new JupiterOrbit(), yearToJD(-4000), yearToJD(4000), astro::SolarMass);
    if (name == "saturn")
        return new MixedOrbit(new SaturnOrbit(), yearToJD(-4000), yearToJD(4000), astro::SolarMass);
    if (name == "uranus")
        return new MixedOrbit(new UranusOrbit(), yearToJD(-4000), yearToJD(4000), astro::SolarMass);
    if (name == "neptune")
        return new MixedOrbit(new NeptuneOrbit(), yearToJD(-4000), yearToJD(4000), astro::SolarMass);
    if (name == "pluto")
        return new MixedOrbit(new PlutoOrbit(), yearToJD(-4000), yearToJD(4000), astro::SolarMass);
    if (name == "mercury-jpl")
        return CreateJPLEphOrbit(JPLEph_Mercury,  0.2408  * 365.25, 6.0e7);
    if (name == "venus-jpl")
        return CreateJPLEphOrbit(JPLEph_Venus,    0.6152  * 365.25, 1.0e8);
    if (name == "earth-jpl")
        return CreateJPLEphOrbit(JPLEph_EarthMoonBary,      365.25, 1.6e8);
    if (name == "mars-jpl")
        return CreateJPLEphOrbit(JPLEph_Mars,      1.8809 * 365.25, 2.4e8);
    if (name == "jupiter-jpl")
        return CreateJPLEphOrbit(JPLEph_Jupiter,  11.86   * 365.25, 8.0e8);
    if (name == "saturn-jpl")
        return CreateJPLEphOrbit(JPLEph_Saturn,   29.4577 * 365.25, 1.5e9);
    if (name == "uranus-jpl")
        return CreateJPLEphOrbit(JPLEph_Uranus,   84.0139 * 365.25, 3.0e9);
    if (name == "neptune-jpl")
        return CreateJPLEphOrbit(JPLEph_Neptune, 164.793  * 365.25, 4.7e9);
    if (name == "pluto-jpl")
        return CreateJPLEphOrbit(JPLEph_Pluto,   248.54   * 365.25, 6.0e9);
    if (name == "moon-jpl")
        return CreateJPLEphOrbit(JPLEph_Moon,     27.321661,        5.0e5);
    if (name == "phobos")
        return new PhobosOrbit();
    if (name == "deimos")
        return new DeimosOrbit();
    if (name == "io")
        return new IoOrbit();
    if (name == "europa")
        return new EuropaOrbit();
    if (name == "ganymede")
        return new GanymedeOrbit();
    if (name == "callisto")
        return new CallistoOrbit();
    if (name == "mimas")
        return new MimasOrbit();
    if (name == "enceladus")
        return new EnceladusOrbit();
    if (name == "tethys")
        return new TethysOrbit();
    if (name == "dione")
        return new DioneOrbit();
    if (name == "rhea")
        return new RheaOrbit();
    if (name == "titan")
        return new TitanOrbit();
    if (name == "hyperion")
        return new HyperionOrbit();
    if (name == "iapetus")
        return new IapetusOrbit();
    if (name == "phoebe")
        return new PhoebeOrbit();
    if (name == "miranda")
        return CreateUranianSatelliteOrbit(1);
    if (name == "ariel")
        return CreateUranianSatelliteOrbit(2);
    if (name == "umbriel")
        return CreateUranianSatelliteOrbit(3);
    if (name == "titania")
        return CreateUranianSatelliteOrbit(4);
    if (name == "oberon")
        return CreateUranianSatelliteOrbit(5);
    if (name == "triton")
        return new TritonOrbit();
    else
        return CreateVSOP87Orbit(name);
}

---