astr2lib.cpp

Go to the documentation of this file.

//
// This file is part of the Marble Virtual Globe.
//
// This program is free software licensed under the GNU LGPL. You can
// find a copy of this license in LICENSE.txt in the top directory of
// the source code.
//
// Copyright 2013 Gerhard Holtkamp
//

/* =========================================================================
  Procedures for calculating the positions of planets.
  The procedures in this unit are taken from Montenbruck, Pfleger
  "Astronomie mit dem Personal Computer", Springer Verlag, 1989
  and modified correspondingly.
  The library Astrolib has to be included.

  Copyright :Gerhard HOLTKAMP          21-SEP-2004   18:00
  ========================================================================= */

#include "astr2lib.h"
#include <cmath>
using namespace std;

#include "attlib.h"
#include "astrolib.h"

extern double frac (double f);
extern double atan21 (double y, double x);


/*-------------------- Class Plan200 --------------------------------------*/
 /*
     Ecliptic coordinates (in A.U.) and velocity (in A.U./day)
     of the Planets for Equinox of Date given in Julian centuries since J2000.
     ======================================================================
  */
Plan200::Plan200 ()
  { }

void Plan200::addthe (double c1, double s1, double c2, double s2,
                            double& cc, double& ss)
 {
    cc=c1*c2-s1*s2;
    ss=s1*c2+c1*s2;
 }

void Plan200::term (int i1, int i, int it, double dlc, double dls, double drc,
              double drs, double dbc, double dbs)
 {
     if (it == 0) addthe (c3[i1],s3[i1],c[i],s[i],u,v);
     else
      {
        u=u*tt;
        v=v*tt;
      }
     dl = dl + dlc*u + dls*v;
     dr = dr + drc*u + drs*v;
     db = db + dbc*u + dbs*v;
 }

Vec3 Plan200::velocity() const   // return last calculated planet velocity
 {
  return vp;
 }

void Plan200::state (Vec3& rs, Vec3& vs) const
 {
  /* State vector rs (position) and vs (velocity) of the Sun in
      ecliptic of date coordinates for last calculated planet
    */
  rs = rp;
  vs = vp;
 }

void Plan200::posvel ()
 {
  /* auxiliary program to calculate position and velocity
      vectors rp, vp.
      to be called by the various planet procedures.
    */

    double  cl, sl, cb, sb;

    cl = cos(l); sl = sin(l); cb = cos(b); sb = sin(b);
    rp[0] = r*cl*cb; rp[1] = r*sl*cb; rp[2] = r*sb;

    // velocity
    vp[0] = (dr*cl*cb - dl*r*sl*cb - db*r*cl*sb) * 1E-4;
    vp[1] = (dr*sl*cb + dl*r*cl*cb - db*r*sl*sb) * 1E-4;
    vp[2] = (dr*sb                 + db*r*cb) * 1E-4;

 }

/*--------------------- Mercury ------------------------------------------*/

Vec3 Plan200::Mercury (double t)   // position of Mercury at time t
 {
 /*
     Ecliptic coordinates rp (in A.U.) and velocity vp (in A.U./day)
     of Mercury for Equinox of Date.
     t is the time in Julian centuries since J2000.
  */

     const double arc = 206264.81;    // 3600*180/pi = ''/r
     const double p2 = M_PI * 2.0;

     int i;

     tt = t;
     dl = 0.0; dr = 0.0; db = 0.0;
     m1 = p2 * frac(0.4855407 + 415.2014314*t);
     m2 = p2 * frac(0.1394222 + 162.5490444*t);
     m3 = p2 * frac(0.9937861 + 99.9978139*t);
     m5 = p2 * frac(0.0558417 + 8.4298417*t);
     m6 = p2 * frac(0.8823333 + 3.3943333*t);
     c3[1] = 1.0; s3[1]= 0.0;
     c3[2] = cos(m1); s3[2] = sin(m1); c3[0] = c3[2]; s3[0] = -s3[2];
     for (i=3; i<11; i++)
        addthe(c3[i-1],s3[i-1],c3[2],s3[2],c3[i],s3[i]);

     // perturbations by Venus
     c[5] = 1.0; s[5] = 0.0; c[4] = cos(m2); s[4] = -sin(m2);
     for (i=4; i>0; i--)
        addthe(c[i],s[i],c[4],s[4],c[i-1],s[i-1]);
     term(2, 5, 0, 259.74, 84547.39, -78342.34,0.01,11683.22,21203.79);
     term(2, 5, 1, 2.3, 5.04, -7.52, 0.02, 138.55, -71.01);
     term(2, 5, 2, 0.01, -0.01, 0.01, 0.01, -0.19, -0.54);
     term(3, 5, 0, -549.71, 10394.44, -7955.45, 0.0, 2390.29, 4306.79);
     term(3, 5, 1, -4.77, 8.97, -1.53, 0.0, 28.49, -14.18);
     term(3, 5, 2, 0.0, 0.0, 0.0, 0.0, -0.04, -0.11);
     term(4, 5, 0, -234.04, 1748.74, -1212.86, 0.0, 535.41, 984.33);
     term(4, 5, 1, -2.03, 3.48, -0.35, 0.0, 6.56, -2.91);
     term(5, 5, 0, -77.64, 332.63, -219.23, 0.0, 124.4, 237.03);
     term(5, 5, 1, -0.7, 1.1, -0.08, 0.0, 1.59, -0.59);
     term(6, 5, 0, -23.59, 67.28, -43.54, 0.0, 29.44, 58.77);
     term(6, 5, 1, -0.23, 0.32, -0.02, 0.0, 0.39, -0.11);
     term(7, 5, 0, -6.86, 14.06, -9.18, 0.0, 7.03, 14.84);
     term(7, 5, 1, -0.07, 0.09, -0.01, 0.0, 0.1, -0.02);
     term(8, 5, 0, -1.94, 2.98, -2.02, 0.0, 1.69, 3.8);
     term(9, 5, 0, -0.54, 0.63, -0.46, 0.0, 0.41, 0.98);
     term(10, 5, 0, -0.15, 0.13, -0.11, 0.0, 0.1, 0.25);
     term(0, 3, 0, -0.17, -0.06, -0.05, 0.14, -0.06, -0.07);
     term(1, 4, 0, 0.24, -0.16, -0.11, -0.16, 0.04, -0.01);
     term(1, 3, 0, -0.68, -0.25, -0.26, 0.73, -0.16, -0.18);
     term(1, 0, 0, 0.37, 0.08, 0.06, -0.28, 0.13, 0.12);
     term(2, 4, 0, 0.58, -0.41, 0.26, 0.36, 0.01, -0.01);
     term(2, 3, 0, -3.51, -1.23, 0.23, -0.63, -0.05, -0.06);
     term(2, 2, 0, 0.08, 0.53, -0.11, 0.04, 0.02, -0.09);
     term(2, 0, 0, 1.44, 0.31, 0.3, -1.39, 0.34, 0.29);
     term(3, 4, 0, 0.15, -0.11, 0.09, 0.12, 0.02, -0.04);
     term(3, 3, 0, -1.99, -0.68, 0.65, -1.91, -0.2, 0.03);
     term(3, 2, 0, -0.34, -1.28, 0.97, -0.26, 0.03, 0.03);
     term(3, 1, 0, -0.33, 0.35, -0.13, -0.13, -0.01, 0.0);
     term(3, 0, 0, 7.19, 1.56, -0.05, 0.12, 0.06, 0.05);
     term(4, 3, 0, -0.52, -0.18, 0.13, -0.39, -0.16, 0.03);
     term(4, 2, 0, -0.11, -0.42, 0.36, -0.1, -0.05, -0.05);
     term(4, 1, 0, -0.19, 0.22, -0.23, -0.2, -0.01, 0.02);
     term(4, 0, 0, 2.77, 0.49, -0.45, 2.56, 0.4, -0.12);
     term(5, 0, 0, 0.67, 0.12, -0.09, 0.47, 0.24, -0.08);
     term(6, 0, 0, 0.18, 0.03, -0.02, 0.12, 0.09, -0.03);

     // perturbations by Earth
     c[4] = cos(m3); s[4] = -sin(m3);
     for (i=4; i>1; i--)
        addthe(c[i],s[i],c[4],s[4],c[i-1],s[i-1]);
     term(1, 1, 0, -0.11, -0.07, -0.08, 0.11, -0.02, -0.04);
     term(2, 4, 0, 0.1, -0.2, 0.15, 0.07, 0.0, 0.0);
     term(2, 3, 0, -0.35, 0.28, -0.13, -0.17, -0.01, 0.0);
     term(2, 1, 0, -0.67, -0.45, 0.0, 0.01, -0.01, -0.01);
     term(3, 3, 0, -0.2, 0.16, -0.16, -0.2, -0.01, 0.02);
     term(3, 2, 0, 0.13, -0.02, 0.02, 0.14, 0.01, 0.0);
     term(3, 1, 0, -0.33, -0.18, 0.17, -0.31, -0.04, 0.0);

     // perturbations by Jupiter
     c[4] = cos(m5); s[4] = -sin(m5);
     for (i=4; i>2; i--)
        addthe(c[i],s[i],c[4],s[4],c[i-1],s[i-1]);
     term(0, 4, 0, -0.08, 0.16, 0.15, 0.08, -0.04, 0.01);
     term(0, 3, 0, 0.1, -0.06, -0.07, -0.12, 0.07, -0.01);
     term(1, 4, 0, -0.31, 0.48, -0.02, 0.13, -0.03, -0.02);
     term(1, 3, 0, 0.42, -0.26, -0.38, -0.5, 0.2, -0.03);
     term(2, 4, 0, -0.7, 0.01, -0.02, -0.63, 0.0, 0.03);
     term(2, 3, 0, 2.61, -1.97, 1.74, 2.32, 0.01, 0.01);
     term(2, 2, 0, 0.32, -0.15, 0.13, 0.28, 0.0, 0.0);
     term(3, 4, 0, -0.18, 0.01, 0.0, -0.13, -0.03, 0.03);
     term(3, 3, 0, 0.75, -0.56, 0.45, 0.6, 0.08, -0.17);
     term(4, 3, 0, 0.2, -0.15, 0.1, 0.14, 0.04, -0.08);

     // perturbations by Saturn
     c[3] = cos(2*m6); s[3] = -sin(2*m6);
     term(2, 3, 0, -0.19, 0.33, 0.0, 0.0, 0.0, 0.0);

     dl = dl + (2.8 + 3.2*t);
     l = p2 * frac(0.2151379 + m1/p2 + ((5601.7+1.1*t)*t+dl)/1296.0e3);
     b = (-2522.15+(-30.18+0.04*t)*t+db) / arc;
     r = 0.3952829 +0.0000016*t + dr*1.0e-6;

     dl = 714.0+292.66*cos(m1)+71.96*cos(2*m1)+18.16*cos(3*m1)
              +4.61*cos(4*m1)+3.81*sin(2*m1)+2.43*sin(3*m1)+1.08*sin(4*m1);
     dr = 55.94*sin(m1)+11.36*sin(2*m1)+2.6*sin(3*m1);
     db = 73.4*cos(m1)+29.82*cos(2*m1)+10.22*cos(3*m1)+3.28*cos(4*m1)
              -40.44*sin(m1)-16.55*sin(2*m1)-5.56*sin(3*m1)-1.72*sin(4*m1);

     posvel();

    return rp;
 }


/*--------------------- Venus ------------------------------------------*/

Vec3 Plan200::Venus (double t)   // position of Venus at time t
 {
 /*
     Ecliptic coordinates rp (in A.U.) and velocity vp (in A.U./day)
     of Venus for Equinox of Date.
     t is the time in Julian centuries since J2000.
  */

     const double arc = 206264.81;    // 3600*180/pi = ''/r
     const double p2 = M_PI * 2.0;

     int i;

     tt = t;
     dl = 0.0; dr = 0.0; db = 0.0;
     m1 = p2 * frac(0.4861431+415.2018375*t);
     m2 = p2 * frac(0.1400197+162.5494552*t);
     m3 = p2 * frac(0.9944153+99.9982208*t);
     m4 = p2 * frac(0.0556297+53.1674631*t);
     m5 = p2 * frac(0.0567028+8.4305083*t);
     m6 = p2 * frac(0.8830539+3.3947206*t);

     c3[0] = 1.0; s3[0] = 0.0; c3[1] = cos(m2); s3[1] = sin(m2);
     for (i=2; i<9; i++)
        addthe(c3[i-1],s3[i-1],c3[1],s3[1],c3[i],s3[i]);

     // perturbations by Mercury
     c[8] = 1.0; s[8] = 0.0; c[7] = cos(m1); s[7] = -sin(m1);
     addthe(c[7],s[7],c[7],s[7],c[6],s[6]);
     term(1, 7, 0, 0.0, 0.0, 0.06, -0.09, 0.01, 0.0);
     term(2, 7, 0, 0.25, -0.09, -0.09, -0.27, 0.0, 0.0);
     term(4, 6, 0, -0.07, -0.08, -0.14, 0.14, -0.01, -0.01);
     term(5, 6, 0, -0.35, 0.08, 0.02, 0.09, 0.0, 0.0);

     // perturbations by Earth
     c[7] = cos(m3); s[7] = -sin(m3);
     for (i=7; i>0; i--)
        addthe(c[i],s[i],c[7],s[7],c[i-1],s[i-1]);
     term(1, 8, 0, 2.37, 2793.23, -4899.07, 0.11, 9995.27, 7027.22);
     term(1, 8, 1, 0.1, -19.65, 34.4, 0.22, 64.95, -86.1);
     term(1, 8, 2, 0.06, 0.04, -0.07, 0.11, -0.55, -0.07);
     term(2, 8, 0, -170.42, 73.13, -16.59, 0.0, 67.71, 47.56);
     term(2, 8, 1, 0.93, 2.91, 0.23, 0.0, -0.03, -0.92);
     term(3, 8, 0, -2.31, 0.9, -0.08, 0.0, 0.04, 2.09);
     term(1, 7, 0, -2.38, -4.27, 3.27, -1.82, 0.0, 0.0);
     term(1, 6, 0, 0.09, 0.0, -0.08, 0.05, -0.02, -0.25);
     term(2, 6, 0, -9.57, -5.93, 8.57, -13.83, -0.01, -0.01);
     term(2, 5, 0, -2.47, -2.4, 0.83, -0.95, 0.16, 0.24);
     term(3, 6, 0, -0.09, -0.05, 0.08, -0.13, -0.28, 0.12);
     term(3, 5, 0, 7.12, 0.32, -0.62, 13.76, -0.07, 0.01);
     term(3, 4, 0, -0.65, -0.17, 0.18, -0.73, 0.1, 0.05);
     term(3, 3, 0, -1.08, -0.95, -0.17, 0.22, -0.03, -0.03);
     term(4, 5, 0, 0.06, 0.0, -0.01, 0.08, 0.14, -0.18);
     term(4, 4, 0, 0.93, -0.46, 1.06, 2.13, -0.01, 0.01);
     term(4, 3, 0, -1.53, 0.38, -0.64, -2.54, 0.27, 0.0);
     term(4, 2, 0, -0.17, -0.05, 0.03, -0.11, 0.02, 0.0);
     term(5, 3, 0, 0.18, -0.28, 0.71, 0.47, -0.02, 0.04);
     term(5, 2, 0, 0.15, -0.14, 0.3, 0.31, -0.04, 0.03);
     term(5, 1, 0, -0.08, 0.02, -0.03, -0.11, 0.01, 0.0);
     term(5, 0, 0, -0.23, 0.0, 0.01, -0.04, 0.0, 0.0);
     term(6, 2, 0, 0.01, -0.14, 0.39, 0.04, 0.0, -0.01);
     term(6, 1, 0, 0.02, -0.05, 0.12, 0.04, -0.01, 0.01);
     term(6, 0, 0, 0.1, -0.1, 0.19, 0.19, -0.02, 0.02);
     term(7, 1, 0, -0.03, -0.06, 0.18, -0.08, 0.0, 0.0);
     term(8, 0, 0, -0.03, -0.02, 0.06, -0.08, 0.0, 0.0);

     // perturbations by Mars
     c[7] = cos(m4); s[7] = -sin(m4);
     for (i=7; i>5; i--)
        addthe(c[i],s[i],c[7],s[7],c[i-1],s[i-1]);
     term(1, 5, 0, -0.65, 1.02, -0.04, -0.02, -0.02, 0.0);
     term(2, 6, 0, -0.05, 0.04, -0.09, -0.1, 0.0, 0.0);
     term(2, 5, 0, -0.5, 0.45, -0.79, -0.89, 0.01, 0.03);

     // perturbations by Jupiter
     c[7] = cos(m5); s[7] = -sin(m5);
     for (i=7; i>5; i--)
        addthe(c[i],s[i],c[7],s[7],c[i-1],s[i-1]);
     term(0, 7, 0, -0.05, 1.56, 0.16, 0.04, -0.08, -0.04);
     term(1, 7, 0, -2.62, 1.4, -2.35, -4.4, 0.02, 0.03);
     term(1, 6, 0, -0.47, -0.08, 0.12, -0.76, 0.04, -0.18);
     term(2, 6, 0, -0.73, -0.51, 1.27, -1.82, -0.01, 0.01);
     term(2, 5, 0, -0.14, -0.1, 0.25, -0.34, 0.0, 0.0);
     term(3, 5, 0, -0.01, 0.04, -0.11, -0.02, 0.0, 0.0);

     // perturbations by Saturn
     c[7] = cos(m6); s[7] = -sin(m6);
     term(0, 7, 0, 0.0, 0.21, 0.0, 0.0, 0.0, -0.01);
     term(1, 7, 0, -0.11, -0.14, 0.24, -0.2, 0.01, 0.0);

     dl = dl+2.74*sin(p2*(0.0764+0.4174*t))+0.27*sin(p2*(0.9201+0.3307*t));
     dl = dl+(1.9+1.8*t);

     l = p2 * frac(0.3654783 + m2/p2 + ((5071.2+1.1*t)*t+dl)/1296.0e3);
     b = (-67.7+(0.04+0.01*t)*t+db)/arc;
     r = 0.7233482-0.0000002*t+dr*1.0e-6;

     dl = 280.0+3.79*cos(m2);
     dr = 1.37*sin(m2);
     db = 9.54*cos(m2)-13.57*sin(m2);

     posvel();

    return rp;
 }

/*--------------------- Mars ------------------------------------------*/

Vec3 Plan200::Mars (double t)   // position of Mars at time t
 {
 /*
     Ecliptic coordinates rp (in A.U.) and velocity vp (in A.U./day)
     of Mars for Equinox of Date.
     t is the time in Julian centuries since J2000.
  */

     const double arc = 206264.81;    // 3600*180/pi = ''/r
     const double p2 = M_PI * 2.0;

     int i;

     tt = t;
     dl = 0.0; dr = 0.0; db = 0.0;
     m2 = p2 * frac(0.1382208+162.5482542*t);
     m3 = p2 * frac(0.9926208+99.9970236*t);
     m4 = p2 * frac(0.0538553+53.1662736*t);
     m5 = p2 * frac(0.0548944+8.4290611*t);
     m6 = p2 * frac(0.8811167+3.3935250*t);
     c3[2] = 1.0; s3[2] = 0.0; c3[3] = cos(m4); s3[3] = sin(m4);
     for (i=4; i<19; i++)
        addthe(c3[i-1],s3[i-1],c3[3],s3[3],c3[i],s3[i]);
     c3[0] = c3[4];
     s3[0] = -s3[4];
     c3[1] = c3[3];
     s3[1] = -s3[3];

     // perturbations by Venus
     c[9] = 1.0; s[9] = 0.0; c[8] = cos(m2); s[8] = -sin(m2);
     addthe (c[8],s[8],c[8],s[8],c[7],s[7]);
     term(2, 8, 0, -0.01, -0.03, 0.1, -0.04, 0.0, 0.0);
     term(3, 8, 0, 0.05, 0.1, -2.08, 0.75, 0.0, 0.0);
     term(4, 8, 0, -0.25, -0.57, -2.58, 1.18, 0.05, -0.04);
     term(4, 7, 0, 0.02, 0.02, 0.13, -0.14, 0.0, 0.0);
     term(5, 8, 0, 3.41, 5.38, 1.87, -1.15, 0.01, -0.01);
     term(5, 7, 0, 0.02, 0.02, 0.11, -0.13, 0.0, 0.0);
     term(6, 8, 0, 0.32, 0.49, -1.88, 1.21, -0.07, 0.07);
     term(6, 7, 0, 0.03, 0.03, 0.12, -0.14, 0.0, 0.0);
     term(7, 8, 0, 0.04, 0.06, -0.17, 0.11, -0.01, 0.01);
     term(7, 7, 0, 0.11, 0.09, 0.35, -0.43, -0.01, 0.01);
     term(8, 7, 0, -0.36, -0.28, -0.2, 0.25, 0.0, 0.0);
     term(9, 7, 0, -0.03, -0.03, 0.11, -0.13, 0.0, 0.0);

     // perturbations by Earth
     c[8] = cos(m3); s[8] = -sin(m3);
     for (i=8; i>0; i--)
        addthe(c[i],s[i],c[8],s[8],c[i-1],s[i-1]);
     term(3, 9, 0, -5.32, 38481.97, -141856.04, 0.4, -6321.67, 1876.89);
     term(3, 9, 1, -1.12, 37.98, -138.67, -2.93, 37.28, 117.48);
     term(3, 9, 2, -0.32, -0.03, 0.12, -1.19, 1.04, -0.4);
     term(4, 9, 0, 28.28, 2285.8, -6608.37, 0.0, -589.35, 174.81);
     term(4, 9, 1, 1.64, 3.37, -12.93, 0.0, 2.89, 11.1);
     term(4, 9, 2, 0.0, 0.0, 0.0, 0.0, 0.1, -0.03);
     term(5, 9, 0, 5.31, 189.29, -461.81, 0.0, -61.98, 18.53);
     term(5, 9, 1, 0.31, 0.35, -1.36, 0.0, 0.25, 1.19);
     term(6, 9, 0, 0.81, 17.96, -38.26, 0.0, -6.88, 2.08);
     term(6, 9, 1, 0.05, 0.04, -0.15, 0.0, 0.02, 0.14);
     term(7, 9, 0, 0.11, 1.83, -3.48, 0.0, -0.79, 0.24);
     term(8, 9, 0, 0.02, 0.2, -0.34, 0.0, -0.09, 0.3);
     term(1, 8, 0, 0.09, 0.06, 0.14, -0.22, 0.02, -0.02);
     term(2, 8, 0, 0.72, 0.49, 1.55, -2.31, 0.12, -0.1);
     term(3, 8, 0, 7.0, 4.92, 13.93, -20.48, 0.08, -0.13);
     term(4, 8, 0, 13.08, 4.89, -4.53, 10.01, -0.05, 0.13);
     term(4, 7, 0, 0.14, 0.05, -0.48, -2.66, 0.01, 0.14);
     term(5, 8, 0, 1.38, 0.56, -2.0, 4.85, -0.01, 0.19);
     term(5, 7, 0, -6.85, 2.68, 8.38, 21.42, 0.0, 0.03);
     term(5, 6, 0, -0.08, 0.2, 1.2, 0.46, 0.0, 0.0);
     term(6, 8, 0, 0.16, 0.07, -0.19, 0.47, -0.01, 0.05);
     term(6, 7, 0, -4.41, 2.14, -3.33, -7.21, -0.07, -0.09);
     term(6, 6, 0, -0.12, 0.33, 2.22, 0.72, -0.03, -0.02);
     term(6, 5, 0, -0.04, -0.06, -0.36, 0.23, 0.0, 0.0);
     term(7, 7, 0, -0.44, 0.21, -0.7, -1.46, -0.06, -0.07);
     term(7, 6, 0, 0.48, -2.6, -7.25, -1.37, 0.0, 0.0);
     term(7, 5, 0, -0.09, -0.12, -0.66, 0.5, 0.0, 0.0);
     term(7, 4, 0, 0.03, 0.0, 0.01, -0.17, 0.0, 0.0);
     term(8, 7, 0, -0.05, 0.03, -0.07, -0.15, -0.01, -0.01);
     term(8, 6, 0, 0.1, -0.96, 2.36, 0.3, 0.04, 0.0);
     term(8, 5, 0, -0.17, -0.2, -1.09, 0.94, 0.02, -0.02);
     term(8, 4, 0, 0.05, 0.0, 0.0, -0.3, 0.0, 0.0);
     term(9, 6, 0, 0.01, -0.1, 0.32, 0.04, 0.02, 0.0);
     term(9, 5, 0, 0.86, 0.77, 1.86, -2.01, 0.01, -0.01);
     term(9, 4, 0, 0.09, -0.01, -0.05, -0.44, 0.0, 0.0);
     term(9, 3, 0, -0.01, 0.02, 0.1, 0.08, 0.0, 0.0);
     term(10, 5, 0, 0.2, 0.16, -0.53, 0.64, -0.01, 0.02);
     term(10, 4, 0, 0.17, -0.03, -0.14, -0.84, 0.0, 0.01);
     term(10, 3, 0, -0.02, 0.03, 0.16, 0.09, 0.0, 0.0);
     term(11, 4, 0, -0.55, 0.15, 0.3, 1.1, 0.0, 0.0);
     term(11, 3, 0, -0.02, 0.04, 0.2, 0.1, 0.0, 0.0);
     term(12, 4, 0, -0.09, 0.03, -0.1, -0.33, 0.0, -0.01);
     term(12, 3, 0, -0.05, 0.11, 0.48, 0.21, -0.01, 0.0);
     term(13, 3, 0, 0.1, -0.35, -0.52, -0.15, 0.0, 0.0);
     term(13, 2, 0, -0.01, -0.02, -0.1, 0.07, 0.0, 0.0);
     term(14, 3, 0, 0.01, -0.04, 0.18, 0.4, 0.01, 0.0);
     term(14, 2, 0, -0.05, -0.07, -0.29, 0.2, 0.01, 0.0);
     term(15, 2, 0, 0.23, 0.27, 0.25, -0.21, 0.0, 0.0);
     term(16, 2, 0, 0.02, 0.03, -0.1, 0.09, 0.0, 0.0);
     term(16, 1, 0, 0.05, 0.01, 0.03, -0.23, 0.0, 0.3);
     term(17, 1, 0, -1.53, 0.27, 0.06, 0.42, 0.0, 0.0);
     term(18, 1, 0, -0.14, 0.02, -0.1, -0.55, -0.01, -0.02);
     term(18, 0, 0, 0.03, -0.06, -0.25, -0.11, 0.0, 0.0);

     // perturbation by Jupiter
     c[8] = cos(m5); s[8] = -sin(m5);
     for (i=8; i>4; i--)
        addthe(c[i],s[i],c[8],s[8],c[i-1],s[i-1]);
     term(0, 8, 0, 0.05, 0.03, 0.08, -0.14, 0.01, -0.01);
     term(1, 8, 0, 0.39, 0.27, 0.92, -1.5, -0.03, -0.06);
     term(1, 7, 0, -0.16, 0.03, 0.13, 0.67, -0.01, 0.06);
     term(1, 6, 0, -0.02, 0.01, 0.05, 0.09, 0.0, 0.01);
     term(2, 8, 0, 3.56, 1.13, -5.41, -7.18, -0.25, -0.24);
     term(2, 7, 0, -1.44, 0.25, 1.24, 7.96, 0.02, 0.31);
     term(2, 6, 0, -0.21, 0.11, 0.55, 1.04, 0.01, 0.05);
     term(2, 5, 0, -0.02, 0.02, 0.11, 0.11, 0.0, 0.01);
     term(3, 8, 0, 16.67, -19.15, 61.0, 53.36, -0.06, -0.07);
     term(3, 7, 0, -21.64, 3.18, -7.77, -54.64, -0.31, 0.5);
     term(3, 6, 0, -2.82, 1.45, -2.53, -5.73, 0.01, 0.07);
     term(3, 5, 0, -0.31, 0.28, -0.34, -0.51, 0.0, 0.0);
     term(4, 8, 0, 2.15, -2.29, 7.04, 6.94, 0.33, 0.19);
     term(4, 7, 0, -15.69, 3.31, -15.7, -73.17, -0.17, -0.25);
     term(4, 6, 0, -1.73, 1.95, -9.19, -7.2, 0.02, -0.03);
     term(4, 5, 0, -0.01, 0.33, -1.42, 0.08, 0.01, -0.01);
     term(4, 4, 0, 0.03, 0.03, -0.13, 0.12, 0.0, 0.0);
     term(5, 8, 0, 0.26, -0.28, 0.73, 0.71, 0.08, 0.04);
     term(5, 7, 0, -2.06, 0.46, -1.61, -6.72, -0.13, -0.25);
     term(5, 6, 0, -1.28, -0.27, 2.21, -6.9, -0.04, -0.02);
     term(5, 5, 0, -0.22, 0.08, -0.44, -1.25, 0.0, 0.01);
     term(5, 4, 0, -0.02, 0.03, -0.15, -0.08, 0.0, 0.0);
     term(6, 8, 0, 0.03, -0.03, 0.08, 0.08, 0.01, 0.01);
     term(6, 7, 0, -0.26, 0.06, -0.17, -0.7, -0.03, -0.05);
     term(6, 6, 0, -0.2, -0.05, 0.22, -0.79, -0.01, -0.02);
     term(6, 5, 0, -0.11, -0.14, 0.93, -0.6, 0.0, 0.0);
     term(6, 4, 0, -0.04, -0.02, 0.09, -0.23, 0.0, 0.0);
     term(7, 5, 0, -0.02, -0.03, 0.13, -0.09, 0.0, 0.0);
     term(7, 4, 0, 0.0, -0.03, 0.21, 0.01, 0.0, 0.0);

     // perturbations by Saturn
     c[8] = cos(m6); s[8] = -sin(m6);
     for (i=8; i>5; i--)
        addthe(c[i],s[i],c[8],s[8],c[i-1],s[i-1]);
     term(1, 8, 0, 0.03, 0.13, 0.48, -0.13, 0.02, 0.0);
     term(2, 8, 0, 0.27, 0.84, 0.4, -0.43, 0.01, -0.01);
     term(2, 7, 0, 0.12, -0.04, -0.33, -0.55, -0.01, -0.02);
     term(2, 6, 0, 0.02, -0.01, -0.07, -0.08, 0.0, 0.0);
     term(3, 8, 0, 1.12, 0.76, -2.66, 3.91, -0.01, 0.01);
     term(3, 7, 0, 1.49, -0.95, 3.07, 4.83, 0.04, -0.05);
     term(3, 6, 0, 0.21, -0.18, 0.55, 0.64, 0.0, 0.0);
     term(4, 8, 0, 0.12, 0.1, -0.29, 0.34, -0.01, 0.02);
     term(4, 7, 0, 0.51, -0.36, 1.61, 2.25, 0.03, 0.01);
     term(4, 6, 0, 0.1, -0.1, 0.5, 0.43, 0.0, 0.0);
     term(4, 5, 0, 0.01, -0.02, 0.11,  0.05, 0.0, 0.0);
     term(5, 7, 0, 0.07, -0.05, 0.16, 0.22, 0.01, 0.01);

     dl = dl+52.49*sin(p2*(0.1868+0.0549*t))+0.61*sin(p2*(0.922+0.3307*t))
              +0.32*sin(p2*(0.4731+2.1485*t))+0.28*sin(p2*(0.9467+1.1133*t));
     dl = dl + (0.14+0.87*t-0.11*t*t);
     l = p2 * frac(0.9334591+m4/p2+((6615.5+1.1*t)*t+dl)/1296.0e3);
     r = 1.5303352 + 0.0000131*t + dr*1.0e-6;
     b = (596.32+(-2.92-0.1*t)*t+db) / arc;

     dl = 91.50 + 17.07*cos(m4) + 2.03*cos(2*m4);
     dr = 12.98*sin(m4) + 1.21*cos(2*m4);
     db = 0.83*cos(m4) + 2.8*sin(m4);

     posvel();

    return rp;
 }

/*--------------------- Jupiter ------------------------------------------*/

Vec3 Plan200::Jupiter (double t)   // position of Jupiter at time t
 {
 /*
     Ecliptic coordinates rp (in A.U.) and velocity vp (in A.U./day)
     of Jupiter for Equinox of Date.
     t is the time in Julian centuries since J2000.
  */

     const double arc = 206264.81;    // 3600*180/pi = ''/r
     const double p2 = M_PI * 2.0;

     int i;

     tt = t;
     dl = 0.0; dr = 0.0; db = 0.0;
     m5 = p2 * frac(0.0565314+8.4302963*t);
     m6 = p2 * frac(0.8829867+3.3947688*t);
     m7 = p2 * frac(0.3969537+1.1902586*t);
     c3[1] = 1.0; s3[1] = 0.0;
     c3[2] = cos(m5); s3[2] = sin(m5); c3[0] = c3[2]; s3[0] = -s3[2];
     for (i=3; i<7; i++)
        addthe(c3[i-1],s3[i-1],c3[2],s3[2],c3[i],s3[i]);

     // perturbations by Saturn
     c[10] = 1.0; s[10] = 0.0; c[9] = cos(m6); s[9] = -sin(m6);
     for (i=9; i>0; i--)
        addthe(c[i],s[i],c[9],s[9],c[i-1],s[i-1]);
     term(0, 9, 0, -0.2, 1.4, 2.0, 0.6, 0.1, -0.2);
     term(1, 9, 0, 9.4, 8.9, 3.9, -8.3, -0.4, -1.4);
     term(1, 8, 0, 5.6, -3.0, -5.4, -5.7, -2.0, 0.0);
     term(1, 7, 0, -4.0, -0.1, 0.0, 5.5, 0.0, 0.0);
     term(1, 5, 0, 3.3, -1.6, -1.6, -3.1, -0.5, -1.2);
     term(2, 10, 0, -113.1, 19998.6, -25208.2, -142.2, -4670.7, 288.9);
     term(2, 10, 1, -76.1, 66.9, -84.2, -95.8, 21.6, 29.4);
     term(2, 10, 2, -0.5, -0.3, 0.4, -0.7, 0.1, -0.1);
     term(2, 9, 0, 78.8, -14.5, 11.5, 64.4, -0.2, 0.2);
     term(2, 8, 0, -2.0, -132.4, 28.8, 4.3, -1.7, 0.4);
     term(2, 8, 1, -1.1, -0.7, 0.2, -0.3, 0.0, 0.0);
     term(2, 7, 0, -7.5, -6.8, -0.4, -1.1, 0.6, -0.9);
     term(2, 6, 0, 0.7, 0.7, 0.6, -1.1, 0.0, -0.2);
     term(2, 5, 0, 51.5, -26.0, -32.5, -64.4, -4.9, -12.4);
     term(2, 5, 1, -1.2, -2.2, -2.7, 1.5, -0.4, 0.3);
     term(3, 10, 0, -3.4, 632.0, -610.6, -6.5, -226.8, 12.7);
     term(3, 10, 1, -4.2, 3.8, -4.1, -4.5, 0.2, 0.6);
     term(3, 9, 0, 5.3, -0.7, 0.7, 6.1, 0.2, 1.1);
     term(3, 8, 0, -76.4, -185.1, 260.2, -108.0, 1.6, 0.0);
     term(3, 7, 0, 66.7, 47.8, -51.4, 69.8, 0.9, 0.3);
     term(3, 7, 1, 0.6, -1.0, 1.0, 0.6, 0.0, 0.0);
     term(3, 6, 0, 17.0, 1.4, -1.8, 9.6, 0.0, -0.1);
     term(3, 5, 0, 1066.2, -518.3, -1.3, -23.9, 1.8, -0.3);
     term(3, 5, 1, -25.4, -40.3, -0.9, 0.3, 0.0, 0.0);
     term(3, 5, 2, -0.7, 0.5, 0.0, 0.0, 0.0, 0.0);
     term(4, 10, 0, -0.1, 28.0, -22.1, -0.2, -12.5, 0.7);
     term(4, 8, 0, -5.0, -11.5, 11.7, -5.4, 2.1, -1.0);
     term(4, 7, 0, 16.9, -6.4, 13.4, 26.9, -0.5, 0.8);
     term(4, 6, 0, 7.2, -13.3, 20.9, 10.5, 0.1, -0.1);
     term(4, 5, 0, 68.5, 134.3, -166.9, 86.5, 7.1, 15.2);
     term(4, 5, 1, 3.5, -2.7, 3.4, 4.3, 0.5, -0.4);
     term(4, 4, 0, 0.6, 1.0, -0.9, 0.5, 0.0, 0.0);
     term(4, 3, 0, -1.1, 1.7, -0.4, -0.2, 0.0, 0.0);
     term(5, 10, 0, 0.0, 1.4, -1.0, 0.0, -0.6, 0.0);
     term(5, 8, 0, -0.3, -0.7, 0.4, -0.2, 0.2, -0.1);
     term(5, 7, 0, 1.1, -0.6, 0.9, 1.2, 0.1, 0.2);
     term(5, 6, 0, 3.2, 1.7, -4.1, 5.8, 0.2, 0.1);
     term(5, 5, 0, 6.7, 8.7, -9.3, 8.7, -1.1, 1.6);
     term(5, 4, 0, 1.5, -0.3, 0.6, 2.4, 0.0, 0.0);
     term(5, 3, 0, -1.9, 2.3, -3.2, -2.7, 0.0, -0.1);
     term(5, 2, 0, 0.4, -1.8, 1.9, 0.5, 0.0, 0.0);
     term(5, 1, 0, -0.2, -0.5, 0.3, -0.1, 0.0, 0.0);
     term(5, 0, 0, -8.6, -6.8, -0.4, 0.1, 0.0, 0.0);
     term(5, 0, 1, -0.5, 0.6, 0.0, 0.0, 0.0, 0.0);
     term(6, 5, 0, -0.1, 1.5, -2.5, -0.8, -0.1, 0.1);
     term(6, 4, 0, 0.1, 0.8, -1.6, 0.1, 0.0, 0.0);
     term(6, 1, 0, -0.5, -0.1, 0.1, -0.8, 0.0, 0.0);
     term(6, 0, 0, 2.5, -2.2, 2.8, 3.1, 0.1, -0.2);

     // perturbations by Uranus
     c[9] = cos(m7); s[9] = -sin(m7);
     addthe(c[9],s[9],c[9],s[9],c[8],s[8]);
     term(2, 9, 0, 0.4, 0.9, 0.0, 0.0, 0.0, 0.0);
     term(2, 8, 0, 0.4, 0.4, -0.4, 0.3, 0.0, 0.0);

     // perturbations by Saturn and Uranus
     double phi, x, y;

     phi = (2*m5-6*m6+3*m7); x = cos(phi); y = sin(phi);
     dl = dl-0.8*x+8.5*y; dr = dr-0.1*x;
     addthe(x,y,c3[2],s3[2],x,y);
     dl = dl+0.4*x+0.5*y; dr = dr-0.7*x+0.5*y; db = db-0.1*x;

     l = p2 * frac(0.0388910 + m5/p2 + ((5025.2+0.8*t)*t+dl)/1296.0e3);
     b = (227.3 - 0.3*t + db) /  arc;
     r = 5.208873 + 0.000041*t + dr*1.0e-5;

     dl = 14.5+1.41*cos(m5); dr = 3.66*sin(m5); db = 0.33*sin(m5);

     posvel();

    return rp;
 }


/*--------------------- Saturn ------------------------------------------*/

Vec3 Plan200::Saturn (double t)   // position of Saturn at time t
 {
 /*
     Ecliptic coordinates rp (in A.U.) and velocity vp (in A.U./day)
     of Saturn for Equinox of Date.
     t is the time in Julian centuries since J2000.
  */

     const double arc = 206264.81;    // 3600*180/pi = ''/r
     const double p2 = M_PI * 2.0;

     int i;

     tt = t;
     dl = 0.0; dr = 0.0; db = 0.0;
     m5 = p2 * frac(0.0565314+8.4302963*t);
     m6 = p2 * frac(0.8829867+3.3947688*t);
     m7 = p2 * frac(0.3969537+1.1902586*t);
     m8 = p2 * frac(0.7208473+0.6068623*t);
     c3[0] = 1.0; s3[0] = 0.0; c3[1] = cos(m6); s3[1] = sin(m6);
     for (i=2; i<12; i++)
        addthe(c3[i-1],s3[i-1],c3[1],s3[1],c3[i],s3[i]);

     // perturbations by Jupiter
     c[6] = 1.0; s[6] = 0.0; c[7] = cos(m5); s[7] = sin(m5);
     for (i=6; i>0; i--)
        addthe(c[i],s[i],c[7],-s[7],c[i-1],s[i-1]);
     term(0, 5, 0, 12.0, -1.4, -13.9, 6.4, 1.2, -1.8);
     term(0, 4, 0, 0.0, -0.2, -0.9, 1.0, 0.0, -0.1);
     term(1, 7, 0, 0.9, 0.4, -1.8, 1.9, 0.2, 0.2);
     term(1, 6, 0, -348.3, 22907.7, -52915.5, -752.2, -3266.5, 8314.4);
     term(1, 6, 1, -225.2, -146.2, 337.7, -521.3, 79.6, 17.4);
     term(1, 6, 2, 1.3, -1.4, 3.2, 2.9, 0.1, -0.4);
     term(1, 5, 0, -1.0, -30.7, 108.6, -815.0, -3.6, -9.3);
     term(1, 4, 0, -2.0, -2.7, -2.1, -11.9, -0.1, -0.4);
     term(2, 7, 0, 0.1, 0.2, -1.0, 0.3, 0.0, 0.0);
     term(2, 6, 0, 44.2, 724.0, -1464.3, -34.7, -188.7, 459.1);
     term(2, 6, 1, -17.0, -11.3, 18.9, -28.9, 1.0, -3.7);
     term(2, 5, 0, -3.5, -426.6, -546.5, -26.5, -1.6, -2.7);
     term(2, 5, 1, 3.5, -2.2, -2.6, -4.3, 0.0, 0.0);
     term(2, 4, 0, 10.5, -30.9, -130.5, -52.3, -1.9, 0.2);
     term(2, 3, 0, -0.2, -0.4, -1.2, -0.1, -0.1, 0.0);
     term(3, 6, 0, 6.5, 30.5, -61.1, 0.4, -11.6, 28.1);
     term(3, 6, 1, -1.2, -0.7, 1.1, -1.8, -0.2, -0.6);
     term(3, 5, 0, 29.0, -40.2, 98.2, 45.3, 3.2, -9.4);
     term(3, 5, 1, 0.6, 0.6, -1.0, 1.3, 0.0, 0.0);
     term(3, 4, 0, -27.0, -21.1, -68.5, 8.1, -19.8, 5.4);
     term(3, 4, 1, 0.9, -0.5, -0.4, -2.0, -0.1, -0.8);
     term(3, 3, 0, -5.4, -4.1, -19.1, 26.2, -0.1, -0.1);
     term(4, 6, 0, 0.6, 1.4, -3.0, -0.2, -0.6, 1.6);
     term(4, 5, 0, 1.5, -2.5, 12.4, 4.7, 1.0, -1.1);
     term(4, 4, 0, -821.9, -9.6, -26.0, 1873.6, -70.5, -4.4);
     term(4, 4, 1, 4.1, -21.9, -50.3, -9.9, 0.7, -3.0);
     term(4, 3, 0, -2.0, -4.7, -19.3, 8.2, -0.1, -0.3);
     term(4, 2, 0, -1.5, 1.3, 6.5, 7.3, 0.0, 0.0);
     term(5, 4, 0, -2627.6, -1277.3, 117.4, -344.1, -13.8, -4.3);
     term(5, 4, 1, 63.0, -98.6, 12.7, 6.7, 0.1, -0.2);
     term(5, 4, 2, 1.7, 1.2, -0.2, 0.3, 0.0, 0.0);
     term(5, 3, 0, 0.4, -3.6, -11.3, -1.6, 0.0, -0.3);
     term(5, 2, 0, -1.4, 0.3, 1.5, 6.3, -0.1, 0.0);
     term(5, 1, 0, 0.3, 0.6, 3.0, -1.7, 0.0, 0.0);
     term(6, 4, 0, -146.7, -73.7, 166.4, -334.3, -43.6, -46.7);
     term(6, 4, 1, 5.2, -6.8, 15.1, 11.4, 1.7, -1.0);
     term(6, 3, 0, 1.5, -2.9, -2.2, -1.3, 0.1, -0.1);
     term(6, 2, 0, -0.7, -0.2, -0.7, 2.8, 0.0, 0.0);
     term(6, 1, 0, 0.0, 0.5, 2.5, -0.1, 0.0, 0.0);
     term(6, 0, 0, 0.3, -0.1, -0.3, -1.2, 0.0, 0.0);
     term(7, 4, 0, -9.6, -3.9, 9.6, -18.6, -4.7, -5.3);
     term(7, 4, 1, 0.4, -0.5, 1.0, 0.9, 0.3, -0.1);
     term(7, 3, 0, 3.0, 5.3, 7.5, -3.5, 0.0, 0.0);
     term(7, 2, 0, 0.2, 0.4, 1.6, -1.3, 0.0, 0.0);
     term(7, 1, 0, -0.1, 0.2, 1.0, 0.5, 0.0, 0.0);
     term(7, 0, 0, 0.2, 0.0, 0.2, -1.0, 0.0, 0.0);
     term(8, 4, 0, -0.7, -0.2, 0.6, -1.2, -0.4, -0.4);
     term(8, 3, 0, 0.5, 1.0, -2.0, 1.5, 0.1, 0.2);
     term(8, 2, 0, 0.4, 1.3, 3.6, -0.9, 0.0, -0.1);
     term(9, 2, 0, 4.0, -8.7, -19.9, -9.9, 0.2, -0.4);
     term(9, 2, 1, 0.5, 0.3, 0.8, -1.8, 0.0, 0.0);
     term(10, 2, 0, 21.3, -16.8, 3.3, 3.3, 0.2, -0.2);
     term(10, 2, 1, 1.0, 1.7, -0.4, 0.4, 0.0, 0.0);
     term(11, 2, 0, 1.6, -1.3, 3.0, 3.7, 0.8, -0.2);

     // perturbations by Uranus
     c[5] = cos(m7); s[5] = -sin(m7);
     for (i=5; i>1; i--)
        addthe(c[i],s[i],c[5],s[5],c[i-1],s[i-1]);
     term(0, 5, 0, 1.0, 0.7, 0.4, -1.5, 0.1, 0.0);
     term(0, 4, 0, 0.0, -0.4, -1.1, 0.1, -0.1, -0.1);
     term(0, 3, 0, -0.9, -1.2, -2.7, 2.1, -0.5, -0.3);
     term(1, 5, 0, 7.8, -1.5, 2.3, 12.7, 0.0, 0.0);
     term(1, 4, 0, -1.1, -8.1, 5.2, -0.3, -0.3, -0.3);
     term(1, 3, 0, -16.4, -21.0, -2.1, 0.0, 0.4, 0.0);
     term(2, 5, 0, 0.6, -0.1, 0.1, 1.2, 0.1, 0.0);
     term(2, 4, 0, -4.9, -11.7, 31.5, -13.3, 0.0, -0.2);
     term(2, 3, 0, 19.1, 10.0, -22.1, 42.1, 0.1, -1.1);
     term(2, 2, 0, 0.9, -0.1, 0.1, 1.4, 0.0, 0.0);
     term(3, 4, 0, -0.4, -0.9, 1.7, -0.8, 0.0, -0.3);
     term(3, 3, 0, 2.3, 0.0, 1.0, 5.7, 0.3, 0.3);
     term(3, 2, 0, 0.3, -0.7, 2.0, 0.7, 0.0, 0.0);
     term(3, 1, 0, -0.1, -0.4, 1.1, -0.3, 0.0, 0.0);

     // perturbations by Neptune
     c[5] = cos(m8); s[5] = -sin(m8);
     addthe(c[5],s[5],c[5],s[5],c[4],s[4]);
     term(1, 5, 0, -1.3, -1.2, 2.3, -2.5, 0.0, 0.0);
     term(1, 4, 0, 1.0, -0.1, 0.1, 1.4, 0.0, 0.0);
     term(2, 4, 0, 1.1, -0.1, 0.2, 3.3, 0.0, 0.0);

     // perturbations by Jupiter and Uranus
     double phi, x, y;

     phi = (-2*m5+5*m6-3*m7); x = cos(phi); y = sin(phi);
     dl = dl-0.8*x-0.1*y; dr = dr-0.2*x+1.8*y; db = db+0.3*x+0.5*y;
     addthe(x,y,c3[1],s3[1],x,y);
     dl = dl+(2.4-0.7*t)*x+(27.8-0.4*t)*y;
     dr = dr+2.1*x-0.2*y;
     addthe(x,y,c3[1],s3[1],x,y);
     dl = dl+0.1*x+1.6*y; dr = dr-3.6*x+0.3*y; db = db-0.2*x+0.6*y;

     l = p2 * frac(0.2561136+m6/p2+((5018.6+t*1.9)*t+dl)/1296.0e3);
     b = (175.1 - 10.2*t + db) / arc;
     r = 9.557584 - 0.000186*t + dr*1.0e-5;

     dl = 5.84+0.65*cos(m6); dr = 3.09*sin(m6); db = 0.24*cos(m6);

     posvel();

    return rp;
 }


/*--------------------- Uranus ------------------------------------------*/

Vec3 Plan200::Uranus (double t)   // position of Uranus at time t
 {
 /*
     Ecliptic coordinates rp (in A.U.) and velocity vp (in A.U./day)
     of Uranus for Equinox of Date.
     t is the time in Julian centuries since J2000.
  */

     const double arc = 206264.81;    // 3600*180/pi = ''/r
     const double p2 = M_PI * 2.0;

     int i;

     tt = t;
     dl = 0.0; dr = 0.0; db = 0.0;
     m5 = p2 * frac(0.0564472+8.4302889*t);
     m6 = p2 * frac(0.8829611+3.3947583*t);
     m7 = p2 * frac(0.3967117+1.1902849*t);
     m8 = p2 * frac(0.7216833+0.6068528*t);
     c3[2] = 1.0; s3[2] = 0.0; c3[3] = cos(m7); s3[3] = sin(m7);
     for (i=4; i<10; i++)
        addthe(c3[i-1],s3[i-1],c3[3],s3[3],c3[i],s3[i]);
     c3[1] = c3[3];
     s3[1] = -s3[3];
     c3[0] = c3[4];
     s3[0] = -s3[4];

     // perturbations by Jupiter
     c[8] = 1.0; s[8] = 0.0; c[7] = cos(m5); s[7] = -sin(m5);
     addthe(c[7],s[7],c[7],s[7],c[6],s[6]);
     term(1, 7, 0, 0.0, 0.0, -0.1, 1.7, -0.1, 0.0);
     term(2, 7, 0, 0.5, -1.2, 18.9, 9.1, -0.9, 0.1);
     term(3, 7, 0, -21.2, 48.7, -455.5, -198.8, 0.0, 0.0);
     term(3, 6, 0, -0.5, 1.2, -10.9, -4.8, 0.0, 0.0);
     term(4, 7, 0, -1.3, 3.2, -23.2, -11.1, 0.3, 0.1);
     term(4, 6, 0, -0.2, 0.2, 1.1, 1.5, 0.0, 0.0);
     term(5, 7, 0, 0.0, 0.2, -1.8, 0.4, 0.0, 0.0);

     // perturbations by Saturn
     c[7] = cos(m6); s[7] = -sin(m6);
     for (i=7; i>4; i--)
        addthe(c[i],s[i],c[7],s[7],c[i-1],s[i-1]);
     term(2, 7, 0, 1.4, -0.5, -6.4, 9.0, -0.4, -0.8);
     term(3, 7, 0, -18.6, -12.6, 36.7, -336.8, 1.0, 0.3);
     term(3, 6, 0, -0.7, -0.3, 0.5, -7.5, 0.1, 0.0);
     term(4, 7, 0, 20.0, -141.6, -587.1, -107.0, 3.1, -0.8);
     term(4, 7, 1, 1.0, 1.4, 5.8, -4.0, 0.0, 0.0);
     term(4, 6, 0, 1.6, -3.8, -35.6, -16.0, 0.0, 0.0);
     term(5, 7, 0, 75.3, -100.9, 128.9, 77.5, -0.8, 0.1);
     term(5, 7, 1, 0.2, 1.8, -1.9, 0.3, 0.0, 0.0);
     term(5, 6, 0, 2.3, -1.3, -9.5, -17.9, 0.0, 0.1);
     term(5, 5, 0, -0.7, -0.5, -4.9, 6.8, 0.0, 0.0);
     term(6, 7, 0, 3.4, -5.0, 21.6, 14.3, -0.8, -0.5);
     term(6, 6, 0, 1.9, 0.1, 1.2, -12.1, 0.0, 0.0);
     term(6, 5, 0, -0.1, -0.4, -3.9, 1.2, 0.0, 0.0);
     term(6, 4, 0, -0.2, 0.1, 1.6, 1.8, 0.0, 0.0);
     term(7, 7, 0, 0.2, -0.3, 1.0, 0.6, -0.1, 0.0);
     term(7, 6, 0, -2.2, -2.2, -7.7, 8.5, 0.0, 0.0);
     term(7, 5, 0, 0.1, -0.2, -1.4, -0.4, 0.0, 0.0);
     term(7, 4, 0, -0.1, 0.0, 0.1, 1.2, 0.0, 0.0);
     term(8, 6, 0, -0.2, -0.6, 1.4, -0.7, 0.0, 0.0);

     // perturbations by Neptune
     c[7] = cos(m8); s[7] = -sin(m8);
     for (i=7; i>0; i--)
        addthe(c[i],s[i],c[7],s[7],c[i-1],s[i-1]);
     term(3, 8, 0, -78.1, 19518.1, -90718.2, -334.7, 2759.5, -311.9);
     term(3, 8, 1, -81.6, 107.7, -497.4, -379.5, -2.8, -43.7);
     term(3, 8, 2, -6.6, -3.1, 14.4, -30.6, -0.4, -0.5);
     term(3, 8, 3, 0.0, -0.5, 2.4, 0.0, 0.0, 0.0);
     term(4, 8, 0, -2.4, 586.1, -2145.2, -15.3, 130.6, -14.3);
     term(4, 8, 1, -4.5, 6.6, -24.2, -17.8, 0.7, -1.6);
     term(4, 8, 2, -0.4, 0.0, 0.1, -1.4, 0.0, 0.0);
     term(5, 8, 0, 0.0, 24.5, -76.2, -0.6, 7.0, -0.7);
     term(5, 8, 1, -0.2, 0.4, -1.4, -0.8, 0.1, -0.1);
     term(6, 8, 0, 0.0, 1.1, -3.0, 0.1, 0.4, 0.0);
     term(1, 7, 0, -0.2, 0.2, 0.7, 0.7, -0.1, 0.0);
     term(2, 7, 0, -2.8, 2.5, 8.7, 10.5, -0.4, -0.1);
     term(3, 7, 0, -28.4, 20.3, -51.4, -72.0, 0.0, 0.0);
     term(3, 6, 0, -0.6, -0.1, 4.2, -14.6, 0.2, 0.4);
     term(3, 5, 0, 0.2, 0.5, 3.4, -1.6, -0.1, 0.1);
     term(4, 7, 0, -1.8, 1.3, -5.5, -7.7, 0.0, 0.3);
     term(4, 6, 0, 29.4, 10.2, -29.0, 83.2, 0.0, 0.0);
     term(4, 5, 0, 8.8, 17.8, -41.9, 21.5, -0.1, -0.3);
     term(4, 4, 0, 0.0, 0.1, -2.1, -0.9, 0.1, 0.0);
     term(5, 6, 0, 1.5, 0.5, -1.7, 5.1, 0.1, -0.2);
     term(5, 5, 0, 4.4, 14.6, -84.3, 25.2, 0.1, -0.1);
     term(5, 4, 0, 2.4, -4.5, 12.0, 6.2, 0.0, 0.0);
     term(5, 3, 0, 2.9, -0.9, 2.1, 6.2, 0.0, 0.0);
     term(6, 5, 0, 0.3, 1.0, -4.0, 1.1, 0.1, -0.1);
     term(6, 4, 0, 2.1, -2.7, 17.9, 14.0, 0.0, 0.0);
     term(6, 3, 0, 3.0, -0.4, 2.3, 17.6, -0.1, -0.1);
     term(6, 2, 0, -0.6, -0.5, 1.1, -1.6, 0.0, 0.0);
     term(7, 4, 0, 0.2, -0.2, 1.0, 0.8, 0.0, 0.0);
     term(7, 3, 0, -0.9, -0.1, 0.6, -7.1, 0.0, 0.0);
     term(7, 2, 0, -0.5, -0.6, 3.8, -3.6, 0.0, 0.0);
     term(7, 1, 0, 0.0, -0.5, 3.0, 0.1, 0.0, 0.0);
     term(8, 2, 0, 0.2, 0.3, -2.7, 1.6, 0.0, 0.0);
     term(8, 1, 0, -0.1, 0.2, -2.0, -0.4, 0.0, 0.0);
     term(9, 1, 0, 0.1, -0.2, 1.3, 0.5, 0.0, 0.0);
     term(9, 0, 0, 0.1, 0.0, 0.4, 0.9, 0.0, 0.0);

     // perturbations by Jupiter and Saturn
     c[7] = cos(m6); s[7] = -sin(m6);
     c[4] = cos(-4*m6+2*m5); s[4] = sin(-4*m6+2*m5);
     for (i=4; i>2; i--)
        addthe(c[i],s[i],c[7],s[7],c[i-1],s[i-1]);
     term(0, 4, 0, -0.7, 0.4, -1.5, -2.5, 0.0, 0.0);
     term(1, 4, 0, -0.1, -0.1, -2.2, 1.0, 0.0, 0.0);
     term(3, 3, 0, 0.1, -0.4, 1.4, 0.2, 0.0, 0.0);
     term(3, 2, 0, 0.4, 0.5, -0.8, -0.8, 0.0, 0.0);
     term(4, 2, 0, 5.7, 6.3, 28.5, -25.5, 0.0, 0.0);
     term(4, 2, 1, 0.1, -0.2, -1.1, -0.6, 0.0, 0.0);
     term(5, 2, 0, -1.4, 29.2, -11.4, 1.1, 0.0, 0.0);
     term(5, 2, 1, 0.8, -0.4, 0.2, 0.3, 0.0, 0.0);
     term(6, 2, 0, 0.0, 1.3, -6.0, -0.1, 0.0, 0.0);

     l = p2 * frac(0.4734843+m7/p2+((5082.3+34.2*t)*t+dl)/1296.0E3);
     b = (-130.61+(-0.54+0.04*t)*t+db)/arc;
     r = 19.211991+(-0.000333-0.000005*t)*t+dr*1.0e-5;

     dl = 2.05+0.19*cos(m7); dr = 1.86*sin(m7); db = -0.03*sin(m7);

     posvel();

    return rp;
 }


/*--------------------- Neptune ------------------------------------------*/

Vec3 Plan200::Neptune (double t)   // position of Neptune at time t
 {
 /*
     Ecliptic coordinates rp (in A.U.) and velocity vp (in A.U./day)
     of Neptune for Equinox of Date.
     t is the time in Julian centuries since J2000.
  */

     const double arc = 206264.81;    // 3600*180/pi = ''/r
     const double p2 = M_PI * 2.0;

     int i;

     tt = t;
     dl = 0.0; dr = 0.0; db = 0.0;
     m5 = p2 * frac(0.0563867+8.4298907*t);
     m6 = p2 * frac(0.8825086+3.3957748*t);
     m7 = p2 * frac(0.3965358+1.1902851*t);
     m8 = p2 * frac(0.7214906+0.6068526*t);
     c3[0] = 1.0; s3[0] = 0.0; c3[1] = cos(m8); s3[1] = sin(m8);
     for (i=2; i<7; i++)
        addthe(c3[i-1],s3[i-1],c3[1],s3[1],c3[i],s3[i]);

     // perturbations by Jupiter
     c[6] = 1.0; s[6] = 0.0; c[5] = cos(m5); s[5] = -sin(m5);
     addthe(c[5],s[5],c[5],s[5],c[4],s[4]);
     term(0, 5, 0, 0.1, 0.1, -3.0, 1.8, -0.3, -0.3);
     term(1, 6, 0, 0.0, 0.0, -15.9, 9.0, 0.0, 0.0);
     term(1, 5, 0, -17.6, -29.3, 416.1, -250.0, 0.0, 0.0);
     term(1, 4, 0, -0.4, -0.7, 10.4, -6.2, 0.0, 0.0);
     term(2, 5, 0, -0.2, -0.4, 2.4, -1.4, 0.4, -0.3);

     // perturbations by Saturn
     c[6] = 1.0; s[6] = 0.0; c[5] = cos(m6); s[5] = -sin(m6);
     addthe(c[5],s[5],c[5],s[5],c[4],s[4]);
     term(0, 5, 0, -0.1, 0.0, 0.2, -1.8, -0.1, -0.5);
     term(1, 6, 0, 0.0, 0.0, -8.3, -10.4, 0.0, 0.0);
     term(1, 5, 0, 13.6, -12.7, 187.5, 201.1, 0.0, 0.0);
     term(1, 4, 0, 0.4, -0.4, 4.5, 4.5, 0.0, 0.0);
     term(2, 5, 0, 0.4, -0.1, 1.7, -3.2, 0.2, 0.2);
     term(2, 4, 0, -0.1, 0.0, -0.2, 2.7, 0.0, 0.0);

     // perturbations by Uranus
     c[6] = 1.0; s[6] = 0.0; c[5] = cos(m7); s[5] = -sin(m7);
     for (i=5; i>0; i--)
        addthe(c[i],s[i],c[5],s[5],c[i-1],s[i-1]);
     term(1, 6, 0, 32.3, 3549.5, -25880.2, 235.8, -6360.5, 374.0);
     term(1, 6, 1, 31.2, 34.4, -251.4, 227.4, 34.9, 29.3);
     term(1, 6, 2, -1.4, 3.9, -28.6, -10.1, 0.0, -0.9);
     term(2, 6, 0, 6.1, 68.0, -111.4, 2.0, -54.7, 3.7);
     term(2, 6, 1, 0.8, -0.2, -2.1, 2.0, -0.2, 0.8);
     term(3, 6, 0, 0.1, 1.0, -0.7, 0.0, -0.8, 0.1);
     term(0, 5, 0, -0.1, -0.3, -3.6, 0.0, 0.0, 0.0);
     term(1, 6, 0, 0.0, 0.0, 5.5, -6.9, 0.1, 0.0);
     term(1, 5, 0, -2.2, -1.6, -116.3, 163.6, 0.0, -0.1);
     term(1, 4, 0, 0.2, 0.1, -1.2, 0.4, 0.0, -0.1);
     term(2, 5, 0, 4.2, -1.1, -4.4, -34.6, -0.2, 0.1);
     term(2, 4, 0, 8.6, -2.9, -33.4, -97.0, 0.2, 0.1);
     term(3, 5, 0, 0.1, -0.2, 2.1, -1.2, 0.0, 0.1);
     term(3, 4, 0, -4.6, 9.3, 38.2, 19.8, 0.1, 0.1);
     term(3, 3, 0, -0.5, 1.7, 23.5, 7.0, 0.0, 0.0);
     term(4, 4, 0, 0.2, 0.8, 3.3, -1.5, -0.2, -0.1);
     term(4, 3, 0, 0.9, 1.7, 17.9, -9.1, -0.1, 0.0);
     term(4, 2, 0, -0.4, -0.4, -6.2, 4.8, 0.0, 0.0);
     term(5, 3, 0, -1.6, -0.5, -2.2, 7.0, 0.0, 0.0);
     term(5, 2, 0, -0.4, -0.1, -0.7, 5.5, 0.0, 0.0);
     term(5, 1, 0, 0.2, 0.0, 0.0, -3.5, 0.0, 0.0);
     term(6, 2, 0, -0.3, 0.2, 2.1, 2.7, 0.0, 0.0);
     term(6, 1, 0, 0.1, -0.1, -1.4, -1.4, 0.0, 0.0);
     term(6, 0, 0, -0.1, 0.1, 1.4, 0.7, 0.0, 0.0);

     l = p2 * frac(0.1254046+m8/p2+((4982.8-21.3*t)*t+dl)/1296.0e3);
     b = (54.77+(0.26+0.06*t)*t+db)/arc;
     r = 30.072984+(0.001234+0.000003*t)*t+dr*1.0e-5;

     dl = 1.04+0.02*cos(m8); dr = 0.27*sin(m8); db = 0.03*sin(m8);

     posvel();

    return rp;
 }

/*--------------------- Pluto ------------------------------------------*/

Vec3 Plan200::Pluto (double t)   // position of Pluto at time t
 {
 /*
     Ecliptic coordinates rp (in A.U.) and velocity vp (in A.U./day)
     of Pluto for Equinox of Date.
     t is the time in Julian centuries since J2000.
  */

     const double arc = 206264.81;    // 3600*180/pi = ''/r
     const double p2 = M_PI * 2.0;

     int i;

     tt = t;
     dl = 0.0; dr = 0.0; db = 0.0;
     m5 = p2 * frac(0.0565314+8.4302963*t);
     m6 = p2 * frac(0.8829867+3.3947688*t);
     m1 = p2 * frac(0.0385795+0.4026667*t);
     c3[0] = 1.0; s3[0] = 0.0; c3[1] = cos(m1); s3[1] = sin(m1);
     for (i=2; i<7; i++)
        addthe(c3[i-1],s3[i-1],c3[1],s3[1],c3[i],s3[i]);

     // Kepler terms and perturbations by Jupiter
     c[3] = 1.0; s[3] = 0.0; c[4] = cos(m5); s[4] = sin(m5);
     for (i=3; i>1; i--)
        addthe(c[i],s[i],c[4],-s[4],c[i-1],s[i-1]);
     addthe(c[4],s[4],c[4],s[4],c[5],s[5]);
     term(1, 3, 0, 0.06,100924.08, -960396.0,15965.1, 51987.68,-24288.76);
     term(2, 3, 0, 3274.74,17835.12, -118252.2,3632.4, 12687.49,-6049.72);
     term(3, 3, 0, 1543.52,4631.99, -21446.6,1167.0, 3504.0,-1853.1);
     term(4, 3, 0, 688.99,1227.08, -4823.4,213.5, 1048.19,-648.26);
     term(5, 3, 0, 242.27, 415.93, -1075.4, 140.6, 302.33, -209.76);
     term(6, 3, 0, 138.41, 110.91, -308.8, -55.3, 109.52, -93.82);
     term(3, 2, 0, -0.99, 5.06, -25.6, 19.8, 1.26, -1.96);
     term(2, 2, 0, 7.15, 5.61, -96.7, 57.2, 1.64, -2.16);
     term(1, 2, 0, 10.79, 23.13, -390.4, 236.4, -0.33, 0.86);
     term(0, 4, 0, -0.23, 4.43, 102.8, 63.2, 3.15, 0.34);
     term(1, 4, 0, -1.1, -0.92, 11.8, -2.3, 0.43, 0.14);
     term(2, 4, 0, 0.62, 0.84, 2.3, 0.7, 0.05, -0.04);
     term(3, 4, 0, -0.38, -0.45, 1.2, -0.8, 0.04, 0.05);
     term(4, 4, 0, 0.17, 0.25, 0.0, 0.2, -0.01, -0.01);
     term(3, 1, 0, 0.06, 0.07, -0.6, 0.3, 0.03, -0.03);
     term(2, 1, 0, 0.13, 0.2, -2.2, 1.5, 0.03, -0.07);
     term(1, 1, 0, 0.32, 0.49, -9.4, 5.7, -0.01, 0.03);
     term(0, 1, 0, -0.04, -0.07, 2.6, -1.5, 0.07, -0.02);

     // perturbations by Saturn
     c[4] = cos(m6); s[4] = sin(m6);
     for (i=3; i>1; i--)
        addthe(c[i],s[i],c[4],-s[4],c[i-1],s[i-1]);
     term(1, 2, 0, -29.47, 75.97, -106.4, -204.9, -40.71, -17.55);
     term(0, 4, 0, -13.88, 18.2, 42.6, -46.1, 1.13, 0.43);
     term(1, 4, 0, 5.81, -23.48, 15.0, -6.8, -7.48, 3.07);
     term(2, 4, 0, -10.27, 14.16, -7.9, 0.4, 2.43, -0.09);
     term(3, 4, 0, 6.86, -10.66, 7.3, -0.3, -2.25, 0.69);
     term(2, 1, 0, 4.23, 2.0, 0.0, -2.2, -0.24, 0.12);
     term(1, 1, 0, -5.04, -0.83, -9.2, -3.1, 0.79, -0.24);
     term(0, 1, 0, 4.25, 2.48, -5.9, -3.3, 0.58, 0.02);

     // perturbations by Jupiter and Saturn
     double x, y;
     y = (m5-m6); x = cos(y); y = sin(y);
     dl = dl-9.11*x+0.12*y; dr = dr-3.4*x-3.3*y; db = db+0.81*x+0.78*y;
     addthe(x,y,c3[1],s3[1],x,y);
     dl = dl+5.92*x+0.25*y; dr = dr+2.3*x-3.8*y; db = db-0.67*x-0.51*y;

     l = p2 * frac(0.6232469+m1/p2+dl/1296.0e3);
     b = -6.8232495189e-2 + db/arc;
     r = 40.7247248+dr*1.0e-5;

     dl = 0.69+0.34*cos(m1)+0.12*cos(2*m1)+0.05*cos(3*m1);
     dr = 6.66*sin(m1)+1.64*sin(2*m1);
     db = -0.08*cos(m1)-0.17*sin(m1)-0.09*sin(2*m1);

     // position
     double cl, sl, cb, sb;
     Mat3 mx;

     cl = cos(l); sl = sin(l); cb = cos(b); sb = sin(b);
     rp[0] = r*cl*cb; rp[1] = r*sl*cb; rp[2] = r*sb;
     mx = pmatecl (-0.5000021, t);   // 1950.0 -> Equinox of Date
     rp = mxvct (mx, rp);

     // velocity
     vp[0] = (dr*cl*cb - dl*r*sl*cb - db*r*cl*sb) * 1E-4;
     vp[1] = (dr*sl*cb + dl*r*cl*cb - db*r*sl*sb) * 1E-4;
     vp[2] = (dr*sb                 + db*r*cb) * 1E-4;

    return rp;
 }

/***************************************************************************/
/* =========================================================================
  Procedures for calculating the positions of various moons in the
  solar system. The algorithms for the procedures of this unit are
  taken from the Explanatory Supplement to the Astronomical Almanac
  (edited by Kenneth Seidelmann, University Science Books, Mill Valley,
  California, 1992. A number of errors in chapter 6 of this book have
  been identified and corrected in the code used here.

  Copyright :Gerhard HOLTKAMP          21-MAY-2000   18:00
  ========================================================================= */


/*--------------------- MarPhobos -----------------------------------------*/
void MarPhobos (double t, Vec3& rs, Vec3& vs)
 {
    /*
     Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
     of the Mars moon Phobos with respect to Mars for Equinox of Date.
     t is the time in Julian centuries since J2000.
     ====================================================================*/

     const double p2 = M_PI * 2.0;
     const double ax =6.26974E-5;
     const double nm = 1128.844556;
     const double ee = 0.0150;
     const double gam = 1.10*M_PI/180.0;

     double thet, ll, pp, na, ja;
     double al, dy, yr;
     Mat3 mt1;

     // time from 11-NOV-1971
     dy = t * 36525.0 + 10278.5;
     yr = dy / 365.25;

     // normalized angles
     thet = p2 * frac ((327.9 - 0.43533*dy) / 360.0);
     ll = p2 * frac ((232.41 + nm*dy + 0.00124*yr*yr) / 360.0);
     pp = p2 * frac ((278.96 + 0.43526*dy) / 360.0);
     na = p2 * frac ((47.39 - 0.0014*yr) / 360.0);
     ja = p2 * frac ((37.27 + 0.0008*yr) / 360.0);

     // in-orbit-plane position vector
     al = ll - pp; // mean anomaly
     al = eccanom (al, ee); // eccentric anomaly
     rs[0] = ax * (cos(al) - ee);
     rs[1] = ax * sqrt(1.0 - ee*ee) * sin(al);
     rs[2] = 0.0;

     // in-orbit-plane velocity vector
     dy = cos(al);
     yr = 1.0 - ee*dy;
     yr = 1.235083014E-3 / yr;
     vs.assign(-yr*sin(al),yr*sqrt(1.0-ee*ee)*dy,0.0);

     // convert to equatorial coordinates
     al = pp - thet - na;
     mt1 = zrot (-al);
     rs =   mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);
     mt1 = xrot (-gam);
     rs =   mxvct (mt1, rs);
     vs = mxvct (mt1, vs);
     mt1 = zrot (-thet);
     rs =   mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);
     mt1 = xrot (-ja);
     rs =   mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);
     mt1 = zrot (-na);
     rs =   mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);

     // convert to ecliptic coordinates of date
     mt1 = pmatequ (-0.500002096, t);
     rs =   mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);
     rs =   equecl (t, rs);
     vs =   equecl (t, vs);

 }

/*--------------------- MarDeimos ------------------------------------------*/
void MarDeimos (double t, Vec3& rs, Vec3& vs)
 {
    /*
     Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
     of the Mars moon Deimos with respect to Mars for Equinox of Date.
     t is the time in Julian centuries since J2000.
     ====================================================================== */

     const double p2 = M_PI * 2.0;
     const double ax =1.56828E-4;
     const double nm = 285.161888;
     const double ee = 0.0004;
     const double gam = 1.79*M_PI/180.0;

     double thet, ll, pp, na, ja;
     double al, dy, yr;
     Mat3 mt1;

     // time from 11-NOV-1971
     dy = t * 36525.0 + 10278.5;
     yr = dy / 365.25;

     // normalized angles
     thet = p2 * frac ((240.38 - 0.01801*dy) / 360.0);
     ll = p2 * frac ((196.55 - 0.01801*dy) / 360.0);
     ll = p2 * frac ((28.96 + nm*dy - 0.27*sin(ll)) / 360.0);
     pp = p2 * frac ((111.7 + 0.01798*dy) / 360.0);
     na = p2 * frac ((46.37 - 0.0014*yr) / 360.0);
     ja = p2 * frac ((36.62 + 0.0008*yr) / 360.0);

     // in-orbit- plane position vector
     al = ll - pp; // mean anomaly
     al = eccanom (al, ee); // eccentric anomaly
     rs[0] = ax * (cos(al) - ee);
     rs[1] = ax * sqrt(1.0 - ee*ee) * sin(al);
     rs[2] = 0.0;

     // in-orbit-plane velocity vector
     dy = cos(al);
     yr = 1.0 - ee*dy;
     yr = 7.8046400669E-4 / yr;
     vs.assign(-yr*sin(al),yr*sqrt(1.0-ee*ee)*dy,0.0);

     // convert to equatorial coordinates
     al = pp - thet - na;
     mt1 = zrot (-al);
     rs =   mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);
     mt1 = xrot (-gam);
     rs =   mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);
     mt1 = zrot (-thet);
     rs =   mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);
     mt1 = xrot (-ja);
     rs =   mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);
     mt1 = zrot (-na);
     rs =   mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);

     // convert to ecliptic coordinates of date
     mt1 = pmatequ (-0.500002096, t);
     rs =   mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);
     rs =   equecl (t, rs);
     vs =   equecl (t, vs);

 }

/*--------------------- PosJIo ------------------------------------------*/
Vec3 PosJIo (double t)
 {
  /* ---------------------------------------------------------------------- }
     Ecliptic coordinates (in A.U.) rs
     of Io with respect to Jupiter for Equinox of Date.
     t is the time in Julian centuries since J2000.
     ======================================================================*/

     const double pi2 = M_PI * 2.0;
     const double ax = 0.002819347;
     const double omj = 99.99754*M_PI/180.0;
     const double jj = 1.30691*M_PI/180.0;
     const double ii = 3.10401*M_PI/180.0;

     double l1, l2, phl, p1, p2, p3, p4;
     double pij, th1, th2, psi, gg;
     double xi, nu, ze;
     double d;  // : EXTENDED;
     Vec3 rs;
     Mat3 mt1;

     // general jupiter moon data
     d = t*36525.0 + 8544.5;
     l1 = pi2 * frac((106.07859 + 203.4889553630643*d)/360.0);
     l2 = pi2 * frac((175.733787 + 101.3747245566245*d)/360.0);
     phl = pi2* frac((184.415351 + 0.17356902*d)/360.0);
     p1 = pi2 * frac((82.380231 + 0.16102275*d)/360.0);
     p2 = pi2 * frac((128.960393 + 0.04645644*d)/360.0);
     p3 = pi2 * frac((187.550171 + 0.00712408*d)/360.0);
     p4 = pi2 * frac((335.309254 + 0.00183939*d)/360.0);
     pij = 0.23510274429; // 13.470395
     th1 = pi2* frac((308.365749 - 0.1328061*d)/360.0);
     th2 = pi2* frac((100.438938 - 0.03261535*d)/360.0);
     psi = pi2* frac((316.500101 - 0.00000248*d)/360.0);
     gg = pi2 * frac((30.2380210168 + 0.08309256178969453*d)/360.0);

     // io specific lines
     xi = -41279.0 * cos(2.0*(l1-l2))*1E-7;
     nu = -5596.0*sin(p3-p4) - 2198.0*sin(p1+p2-2.0*(pij+gg))
                +1321.0*sin(phl) - 1157.0*sin(l1-2.0*l2+p4)
            -1940.0*sin(l1-2.0*l2+p3) + 791.0*sin(l1-2.0*l2+p2)
                + 82363.0*sin(2.0*(l1-l2));
     nu *= 1E-7;
     ze = 7038.0*sin(l1-th1+nu) + 1835.0*sin(l1-th2+nu);
     ze *= 1E-7;

     // convert to rectangular coordinates in Jovian equatorial

     rs.assign (ax * (1+xi) * cos(l1-psi+nu),
                    ax * (1+xi) * sin(l1-psi+nu),
                    ax * ze);

     // convert into ecliptic of 1950.0
     mt1 = xrot (-ii);
     rs =   mxvct (mt1, rs);
     mt1 = zrot (-(psi-omj));
     rs =   mxvct (mt1, rs);
     mt1 = xrot (-jj);
     rs =   mxvct (mt1, rs);
     mt1 = zrot (-omj);
     rs =   mxvct (mt1, rs);

     // convert into ecliptic of date
     mt1 = pmatecl (-0.500002096, t);
     rs =   mxvct (mt1, rs);

    return rs;

 }

/*--------------------- PosEuropa ------------------------------------------*/
Vec3 PosEuropa (double t)
 {
  /*
     Ecliptic coordinates (in A.U.) rs
     of Europa with respect to Jupiter for Equinox of Date.
     t is the time in Julian centuries since J2000.
     =====================================================================*/

     const double pi2 = M_PI * 2.0;
     const double ax = 0.004485872;
     const double omj = 99.99754*M_PI/180.0;
     const double jj = 1.30691*M_PI/180.0;
     const double ii = 3.10401*M_PI/180.0;

     double l1, l2, l3, phl, p2, p3, p4;
     double pij, th2, th3, psi, gg;
     double xi, nu, ze;
     double d; // : EXTENDED;
     Vec3 rs;
     Mat3 mt1;


     // general jupiter moon data
     d = t*36525.0 + 8544.5;
     l1 = pi2 * frac((106.07859 + 203.4889553630643*d)/360.0);
     l2 = pi2 * frac((175.733787 + 101.3747245566245*d)/360.0);
     l3 = pi2 * frac((120.5613855 + 50.31760915340462*d)/360.0);
     phl = pi2* frac((184.415351 + 0.17356902*d)/360.0);
     // p1 = pi2 * frac((82.380231 + 0.16102275*d)/360.0);
     p2 = pi2 * frac((128.960393 + 0.04645644*d)/360.0);
     p3 = pi2 * frac((187.550171 + 0.00712408*d)/360.0);
     p4 = pi2 * frac((335.309254 + 0.00183939*d)/360.0);
     pij = 0.23510274429; // 13.470395
     th2 = pi2* frac((100.438938 - 0.03261535*d)/360.0);
     th3 = pi2* frac((118.908928 - 0.00717678*d)/360.0);
     psi = pi2* frac((316.500101 - 0.00000248*d)/360.0);
     gg = pi2 * frac((30.2380210168 + 0.08309256178969453*d)/360.0);

     // Europa specific lines
     xi = -3187*cos(l2-p3) - 1738*cos(l2-p4)
                +93748*cos(l1-l2);
     xi *= 1E-7;
     nu = -1158*sin(2*(psi-pij)) + 1715*sin(-2*(pij+gg)+th3+psi)
                -1846*sin(gg) + 2397*sin(p3-p4) - 3172*sin(phl)
                -1993*sin(l2-l3) + 1844*sin(l2-p2)
                +6394*sin(l2-p3) + 3451*sin(l2-p4)
                +4159*sin(l1-2*l2+p4) + 7571*sin(l1-2*l2+p3)
                -1491*sin(l1-2*l2+p3) - 185640*sin(l1-l2)
                -803*sin(l1-l3) + 915*sin(2*(l1-l2));
     nu *= 1E-7;
     ze = 81575*sin(l2-th2+nu) + 4512*sin(l2-th3+nu)
                -3286*sin(l2-psi+nu);
     ze *= 1E-7;

     // convert to rectangular coordinates in Jovian equatorial

     rs.assign (ax * (1+xi) * cos(l2-psi+nu),
                    ax * (1+xi) * sin(l2-psi+nu),
                    ax * ze);

     // convert into ecliptic of 1950.0

     mt1 = xrot (-ii);
     rs =   mxvct (mt1, rs);
     mt1 = zrot (-(psi-omj));
     rs =   mxvct (mt1, rs);
     mt1 = xrot (-jj);
     rs =   mxvct (mt1, rs);
     mt1 = zrot (-omj);
     rs =   mxvct (mt1, rs);

     // convert into ecliptic of date
     mt1 = pmatecl (-0.500002096, t);
     rs =   mxvct (mt1, rs);

    return rs;

  }

/*--------------------- PosGanymede --------------------------------------*/
Vec3 PosGanymede (double t)
 {
  /*
     Ecliptic coordinates (in A.U.)
     of Ganymede with respect to Jupiter for Equinox of Date.
     t is the time in Julian centuries since J2000.
     =====================================================================*/

     const double pi2 = M_PI * 2.0;
     const double ax = 0.007155352;
     const double omj = 99.99754*M_PI/180.0;
     const double jj = 1.30691*M_PI/180.0;
     const double ii = 3.10401*M_PI/180.0;

     double l1, l2, l3, l4, phl, p1, p2, p3, p4;
     double pij, th2, th3, th4, psi, gg;
     double xi, nu, ze;
     double d;  // : EXTENDED;
     Vec3 rs;
     Mat3 mt1;

     // general jupiter moon data
     d = t*36525.0 + 8544.5;
     l1 = pi2 * frac((106.07859 + 203.4889553630643*d)/360.0);
     l2 = pi2 * frac((175.733787 + 101.3747245566245*d)/360.0);
     l3 = pi2 * frac((120.5613855 + 50.31760915340462*d)/360.0);
     l4 = pi2 * frac((84.455823 + 21.57107087517961*d)/360.0);
     phl = pi2* frac((184.415351 + 0.17356902*d)/360.0);
     p1 = pi2 * frac((82.380231 + 0.16102275*d)/360.0);
     p2 = pi2 * frac((128.960393 + 0.04645644*d)/360.0);
     p3 = pi2 * frac((187.550171 + 0.00712408*d)/360.0);
     p4 = pi2 * frac((335.309254 + 0.00183939*d)/360.0);
     pij = 0.23510274429; // 13.470395
     th2 = pi2* frac((100.438938 - 0.03261535*d)/360.0);
     th3 = pi2* frac((118.908928 - 0.00717678*d)/360.0);
     th4 = pi2* frac((322.746564 - 0.00176018*d)/360.0);
     psi = pi2* frac((316.500101 - 0.00000248*d)/360.0);
     gg = pi2 * frac((30.2380210168 + 0.08309256178969453*d)/360.0);

     // Ganymede specific lines
     xi = -14691*cos(l3-p3) - 1758*cos(2*(l3-l4))
                +6333*cos(l2-l3);
     xi *= 1E-7;
     nu = -1488*sin(2*(pij+psi)) + 411*sin(-th3+psi)
                +346*sin(-th4+psi) - 2338*sin(gg)
                +6558*sin(p3-p4) + 523*sin(p1+p3-2*(pij+gg))
                +341*sin(phl) - 943*sin(l3-l4)
                +29387*sin(l3-p3) + 15800*sin(l3-p4)
                +3218*sin(2*(l3-l4)) + 226*sin(3*l3-2*l4)
                -12038*sin(l2-l3) - 662*sin(l1-2*l2+p4)
                -1246*sin(l1-2*l2+p3) + 699*sin(l1-2*l2+p2)
                +217*sin(l1-l3);
     nu *= 1E-7;
     ze = -2793*sin(l3-th2+nu) + 32387*sin(l3-th3+nu)
                +6871*sin(l3-th4+nu) - 16876*sin(l3-psi+nu);
     ze *= 1E-7;

     // convert to rectangular coordinates in Jovian equatorial

     rs.assign (ax * (1+xi) * cos(l3-psi+nu),
                    ax * (1+xi) * sin(l3-psi+nu),
                    ax * ze);

     // convert into ecliptic of 1950.0

     mt1 = xrot (-ii);
     rs =   mxvct (mt1, rs);
     mt1 = zrot (-(psi-omj));
     rs =   mxvct (mt1, rs);
     mt1 = xrot (-jj);
     rs =   mxvct (mt1, rs);
     mt1 = zrot (-omj);
     rs =   mxvct (mt1, rs);

     // convert into ecliptic of date
     mt1 = pmatecl (-0.500002096, t);
     rs =   mxvct (mt1, rs);

    return rs;

 }

/*--------------------- PosCallisto --------------------------------------*/
Vec3 PosCallisto (double t)
 {
  /*
     Ecliptic coordinates (in A.U.)
     of Callisto with respect to Jupiter for Equinox of Date.
     t is the time in Julian centuries since J2000.
     =====================================================================*/

     const double pi2 = M_PI * 2.0;
     const double ax = 0.012585436;
     const double omj = 99.99754*M_PI/180.0;
     const double jj = 1.30691*M_PI/180.0;
     const double ii = 3.10401*M_PI/180.0;

     double l3, l4, p3, p4;
     double pij, th3, th4, psi, gp, gg, ph2;
     double xi, nu, ze;
     double d;  // : EXTENDED;
     Vec3 rs;
     Mat3 mt1;

     // general jupiter moon data
     d = t*36525.0 + 8544.5;
     l3 = pi2 * frac((120.5613855 + 50.31760915340462*d)/360.0);
     l4 = pi2 * frac((84.455823 + 21.57107087517961*d)/360.0);
     // phl = pi2* frac((184.415351 + 0.17356902*d)/360.0);
     p3 = pi2 * frac((187.550171 + 0.00712408*d)/360.0);
     p4 = pi2 * frac((335.309254 + 0.00183939*d)/360.0);
     pij = 0.23510274429; // 13.470395
     th3 = pi2* frac((118.908928 - 0.00717678*d)/360.0);
     th4 = pi2* frac((322.746564 - 0.00176018*d)/360.0);
     psi = pi2* frac((316.500101 - 0.00000248*d)/360.0);
     gp = pi2 * frac((31.9785280244 + 0.033459733896*d)/360.0);
     gg = pi2 * frac((30.2380210168 + 0.08309256178969453*d)/360.0);
     ph2 = 0.91009489942; // 52.1445966929

     // Callisto specific lines
     xi = 1656*cos(l4-p3) - 73328*cos(l4-p4)
                +182*cos(l4-pij) - 541*cos(l4+p4-2*(pij+gg))
                -269*cos(2*(l4-p4)) + 974*cos(l3-l4);
     xi *= 1E-7;
     nu = -407*sin(2*(psi-p4)) + 309*sin(-2*p4+th4+psi)
                -4840*sin(2*(psi-pij)) + 2074*sin(-th4+psi)
                -5605*sin(gg) - 204*sin(2*gg)
                -495*sin(5*gp-2*gg+ph2) + 234*sin(p4-pij)
                -6112*sin(p3-p4) - 3318*sin(l4-p3)
                +145573*sin(l4-p4) + 178*sin(l4-pij-gg)
                -363*sin(l4-pij) + 1085*sin(l4+p4-2*(pij+gg))
                +672*sin(2*(l4-p4)) + 218*sin(2*(l4-pij-gg))
                +167*sin(2*l4-th4-psi) - 142*sin(2*(l4-psi))
                +148*sin(l3-2*l4+p4) - 390*sin(l3-l4)
                -195*sin(2*(l3-l4)) + 185*sin(3*l3-7*l4+4*p4);
     nu *= 1E-7;
     ze = 773*sin(l4-2*(pij+gg)+psi+nu)
                -5075*sin(l4-th3+nu) + 44300*sin(l4-th4+nu)
                -76493*sin(l4-psi+nu);
     ze *= 1E-7;

     rs.assign(ax * (1+xi) * cos(l4-psi+nu),
                  ax * (1+xi) * sin(l4-psi+nu),
                  ax * ze);

     // convert into ecliptic of 1950.0

     mt1 = xrot (-ii);
     rs =   mxvct (mt1, rs);
     mt1 = zrot (-(psi-omj));
     rs =   mxvct (mt1, rs);
     mt1 = xrot (-jj);
     rs =   mxvct (mt1, rs);
     mt1 = zrot (-omj);
     rs =   mxvct (mt1, rs);

     // convert into ecliptic of date
     mt1 = pmatecl (-0.500002096, t);
     rs =   mxvct (mt1, rs);

    return rs;

 }


/*--------------------- JupIo --------------------------------------*/
void JupIo (double t, Vec3& rs, Vec3& vs)
 {
  /*
     Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
     of Io with respect to Jupiter for Equinox of Date.
     t is the time in Julian centuries since J2000.
     ====================================================================*/

     const double dlt = 3.472222E-4;  // 30 sec

     double td;

     td = dlt / 36525;
     rs = PosJIo (t-td);
     vs = PosJIo (t+td);
     vs -= rs;
     td = 0.5/dlt;
     vs *= td;
     rs = PosJIo (t);

  }

/*--------------------- JupEuropa --------------------------------------*/
void JupEuropa (double t, Vec3& rs, Vec3& vs)
 {
  /*
     Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
     of Europa with respect to Jupiter for Equinox of Date.
     t is the time in Julian centuries since J2000.
     ====================================================================*/

     const double dlt = 3.472222E-4;  // 30 sec

     double td;

     td = dlt / 36525;
     rs = PosEuropa (t-td);
     vs = PosEuropa (t+td);
     vs -= rs;
     td = 0.5/dlt;
     vs *= td;
     rs = PosEuropa (t);

  }

/*--------------------- JupGanymede ------------------------------------*/
void JupGanymede (double t, Vec3& rs, Vec3& vs)
 {
  /*
     Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
     of Ganimede with respect to Jupiter for Equinox of Date.
     t is the time in Julian centuries since J2000.
     ====================================================================*/

     const double dlt = 3.472222E-4;  // 30 sec

     double td;

     td = dlt / 36525;
     rs = PosGanymede (t-td);
     vs = PosGanymede (t+td);
     vs -= rs;
     td = 0.5/dlt;
     vs *= td;
     rs = PosGanymede (t);

  }

/*--------------------- JupCallisto ------------------------------------*/
void JupCallisto (double t, Vec3& rs, Vec3& vs)
 {
  /*
     Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
     of Callisto with respect to Jupiter for Equinox of Date.
     t is the time in Julian centuries since J2000.
     ====================================================================*/

     const double dlt = 3.472222E-4;  // 30 sec

     double td;

     td = dlt / 36525;
     rs = PosCallisto (t-td);
     vs = PosCallisto (t+td);
     vs -= rs;
     td = 0.5/dlt;
     vs *= td;
     rs = PosCallisto (t);

  }


/*----------------------- SatRhea ------------------------------------*/
void SatRhea (double t, Vec3& rs, Vec3& vs)
 {
  /*
     Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
     of Rhea with respect to Saturn for Equinox of Date.
     t is the time in Julian centuries since J2000.
     ====================================================================*/

     const double pi2 = M_PI * 2.0;
     const double ie = 28.0771*M_PI/180.0;
     const double oe = 168.8263*M_PI/180.0;
     const double kap = 57.29578*M_PI/180.0;
     const double ax = 0.003524;
     const double sg0 = 5.7682811893E-3;  // sin(0.3305)

     double d, td, tt, pp, ot, nt;
     double lam, ii, omg, omb, ee;
     Mat3 mt1;

     // various times
     tt = t*36525.0;
     d = tt + 40452.0;
     td = 0.5219*(tt+40177)/365.25;
     tt = d / 365.25;

     // Titan perturbations
     pp = pi2 * frac((305.0 + 10.2077*tt)/360.0);
     ot = pi2 * frac((276.49 + td)/360.0);
     nt = pi2 * frac((44.5 - td)/360.0);

     // osculating orbit elements
     ee = 0.00021*sin(pp) + 0.001*sin(ot);
     omg = 0.00021*cos(pp) + 0.001*cos(ot);
     omb = atan2(ee,omg);
     ii = sin(omb);
     if (fabs(ii) > 0.5) ee = fabs(ee/ii);
     else ee = fabs(omg/cos(omb));
     pp = pi2 * frac((356.87 - 10.2077*tt)/360.0);
     ot = sin(pp);
     lam = pi2 * frac((359.4727 + 79.69004007*d)/360.0);
     lam = lam + kap*sg0*tan(0.5*ie)*ot;
     ii = ie - 7.9412480966E-4 + kap*sg0*cos(pp)
                + 3.5081117965E-4*cos(nt);
     omg = oe - 1.3613568166E-4
                 + (kap*sg0*ot + 3.5081117965E-4*sin(nt))/sin(ie);

     // in-orbit-plane position vector
     pp = lam - omb; // mean anomaly
     pp = eccanom (pp, ee);  // eccentric anomaly
     rs.assign (ax * (cos(pp) - ee), ax * sqrt(1.0 - ee*ee) * sin(pp), 0);

     // in-orbit-plane velocity vector
     d = cos(pp);
     td = 1.0 - ee*d;
     td = 4.9000413575E-3 / td;
     vs.assign (-td*sin(pp), td*sqrt(1.0-ee*ee)*d, 0);

     // convert to ecliptic coordinates
     pp = omb - omg;
     mt1 = zrot (-pp);
     rs =   mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);
     mt1 = xrot (-ii);
     rs =   mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);
     mt1 = zrot (-omg);
     rs =   mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);

     // convert to ecliptic coordinates of date
     mt1 = pmatecl (-0.500002096, t);
     rs =   mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);

  }

/*----------------------- SatTitan ------------------------------------*/
void SatTitan (double t, Vec3& rs, Vec3& vs)
 {
  /*
     Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
     of Titan with respect to Saturn for Equinox of Date.
     t is the time in Julian centuries since J2000.
     ===================================================================*/

     const double pi2 = M_PI * 2.0;
     const double ie = 28.0771*M_PI/180.0;
     const double oe = 168.8263*M_PI/180.0;
     const double kap = 57.29578*M_PI/180.0;
     const double ax = 0.00816765;
     const double sg0 = 5.7682811893E-3;  // sin(0.3305)

     double d, tt, ts, gg, ls, is, os, lms, psi;
     double lam, ii, omg, omb, ee, s4;
     Mat3 mt1;

     // various times
     ts = t + 1.0;
     tt = t*36525.0;
     d = tt+40177;
     tt = d / 365.25;

     // solar perturbations
     ls = pi2 * frac((175.4762 + 1221.5515*ts)/360.0);
     is = pi2 * frac((2.4891 + 0.002435*ts)/360.0);
     os = pi2 * frac((113.35 - 0.2597*ts)/360.0);
     lms = pi2* frac((267.2635 + 1222.1136*ts)/360.0);

     gg = pi2 * frac((41.28 - 0.5219*tt)/360.0);
     ii = ie - 1.0828022679E-2 + kap*5.2185107774E-3*cos(gg);
     s4 = sin(gg);
     omg = oe - 2.4748768793E-3 + kap*5.2185107774E-3*s4/sin(ie);
     omb = pi2 * frac((275.837 + 0.5219*tt)/360.0);

     psi = sin(is)*sin(omg-os);
     gg = cos(omg-os);
     ts = cos(is)*sin(ii) - sin(is)*cos(ii)*gg;
     psi = atan2(psi,ts);
     ts = cos(is)*sin(ii)*gg - sin(is)*cos(ii);
     gg = sin(ii)*sin(omg-os);
     gg = atan2(gg,ts);
     lms = lms - gg - os;
     gg = omb - omg - psi;

     // osculating orbit elements
     ts = 2*(lms-gg);
     ee = 0.028815 - 0.000184*cos(2*gg) + 0.000073*cos(ts);
     omb = omb + kap*(0.00630*sin(2*gg) + 0.0025*sin(ts));
     gg = 2*lms + psi;
     lam = pi2 * frac((261.3121 + 22.57697385*d)/360.0)
                 +kap*(sg0*tan(0.5*ie)*s4 - 0.000176*sin(ls)
                 -0.000215*sin(2*lms) + 0.000057*sin(gg));
     ii = ii + 0.000232*kap*cos(gg);
     omg = omg + 0.000503*kap*sin(gg);

     // in-orbit-plane position vector
     gg = lam - omb; // mean anomaly
     if (gg < 0) gg += pi2;
     gg = eccanom (gg, ee); // eccentric anomaly
     rs.assign (ax * (cos(gg) - ee), ax * sqrt(1.0 - ee*ee) * sin(gg), 0);

     // in-orbit-plane velocity vector
     d = cos(gg);
     ts = 1.0 - ee*d;
     ts = 3.2183196288E-3 / ts;
     vs.assign (-ts*sin(gg), ts*sqrt(1.0-ee*ee)*d, 0);

     // convert to ecliptic coordinates
     gg = omb - omg;
     mt1 = zrot (-gg);
     rs =   mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);
     mt1 = xrot (-ii);
     rs =   mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);
     mt1 = zrot (-omg);
     rs =   mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);

     // convert to ecliptic coordinates of date
     mt1 = pmatecl (-0.500002096, t);
     rs = mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);

  }


/*----------------------- NepTriton ------------------------------------*/
void NepTriton (double t, Vec3& rs, Vec3& vs)
 {
  /*
     Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
     of Triton with respect to Neptune for Equinox of Date.
     t is the time in Julian centuries since J2000.
     ====================================================================*/

     const double pi2 = M_PI * 2.0;
     const double ax = 0.00237142;
     const double gam = 158.996*M_PI/180.0;

     double d, nn, u, jj, ap;
     Mat3 mt1;

     // days since epoch
     d = t*36525 + 18262.5;

     // pole of fixed plane in 1950.0
     nn = pi2 * frac((359.28 + 54.308*t)/360.0);
     rs.assign (1.0,
                    pi2 * frac((298.72+2.58*sin(nn)-0.04*sin(2*nn))/360.0),
                    pi2 * frac((42.63-1.9*cos(nn)+0.01*cos(2*nn))/360.0));
     rs =   polcar (rs);
     mt1 = pmatequ (0.0, -0.500002096);
     rs =   mxvct (mt1, rs);
     rs =   carpol (rs);
     jj = M_PI/2.0;
     ap = rs[1] + jj;   // R.A. of ascending node
     jj = jj - rs[2];   // inclination

     // in-orbit-plane position vector

     u = pi2 * frac((200.913 + 61.2588532*d)/360.0);
     rs.assign (ax * cos(u), ax * sin(u), 0);

     // velocity is perpendicular to radius vector for e=0
     u = 1.0 / abs(rs);
     vs = u * rs;
     mt1 = zrot (-M_PI/2.0);
     vs =   mxvct (mt1, vs);
     u = 2.53538612E-3;  // mean orbital speed
     vs *= u;

     // convert to equatorial coordinates
     mt1 = xrot (-gam);
     rs =   mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);
     nn = pi2 * frac((151.401 + 0.57806*d/365.25)/360.0);
     mt1 = zrot (-nn);
     rs =   mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);
     mt1 = xrot (-jj);
     rs =   mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);
     mt1 = zrot (-ap);
     rs = mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);

     // convert to ecliptic coordinates of date
     mt1 = pmatequ (-0.500002096, t);
     rs =   mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);
     rs =   equecl (t, rs);
     vs =   equecl (t, vs);

  }


/*----------------------- PluCharon ------------------------------------*/
void PluCharon (double t, Vec3& rs, Vec3& vs)
 {
  /*
     Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
     of Charon with respect to Pluto for Equinox of Date.
     t is the time in Julian centuries since J2000.
     ====================================================================*/

     const double pi2 = M_PI * 2.0;
     const double ax = 0.00013102;
     const double jj = 94.3*M_PI/180.0;
     const double nn = 223.7*M_PI/180.0;

     double d, u;
     Mat3 mt1;

     // days since epoch
     d = t*36525 + 6544.5;

     // in-orbit-plane position vector
     u = pi2 * frac((78.6 + 56.3625*d)/360.0);
     rs.assign (ax * cos(u), ax * sin(u), 0);

     // velocity is perpendicular to radius vector for e=0
     d = 1.0 / abs(rs);
     vs = d * rs;
     mt1 = zrot (-M_PI/2.0);
     vs =  mxvct (mt1, vs);
     d = 1.288837578E-4;  // mean orbital speed
     vs *= d;

     // convert to equatorial coordinates
     mt1 = xrot (-jj);
     rs =   mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);
     mt1 = zrot (-nn);
     rs =   mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);

     // convert to ecliptic coordinates of date
     mt1 = pmatequ (-0.500002096, t);
     rs =   mxvct (mt1, rs);
     vs =   mxvct (mt1, vs);
     rs =   equecl (t, rs);
     vs =   equecl (t, vs);

  }