File indexing completed on 2024-04-28 07:36:49

 // SPDX-License-Identifier: LGPL-2.1-or-later
 //
 // SPDX-FileCopyrightText: 2014 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          26-MAR-2014   
   ========================================================================= */
 
 #include "astr2lib.h"
 #include <cmath>
 using namespace std;
 
 #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()   // return last calculated planet velocity
 {
  return vp;
 }
 
 void Plan200::state (Vec3& rs, Vec3& vs)
 {
  /* 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          25-MAR-2014   
   ========================================================================= */
 
 
 /*--------------------- 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;
 
  }
 
 /*-------------------- Mimas, Enceladus, Dione ---------------------*/
 Vec3 PosSMimas (double t)
 {
  /* Ecliptic coordinates (in A.U.) of Mimas with respect to Saturn 
     for Equinox of Date.
     t is the time in Julian centuries since J2000.
   ===================================================================*/
  const double ie = 28.0771*M_PI/180.0;
  const double oe = 168.8263*M_PI/180.0;
 
  double d, dt, tt, a, n, ecc, gam, th1, l1, p, m;
  Vec3 rs;
  Mat3 mt1; 
 
  d = 36525.0*t + 40452.0;
  dt = d / 365.25;
  tt = 5.0616*(((36525*t + 18262.577000)/365.25)+1950.0-1866.06);
  tt *= M_PI / 180.0;
  a = 0.00124171;
  n = 381.994516;
  ecc = 0.01986;
  gam = 0.0274017; // 1.570;
  th1 = 49.4 - 365.025*dt;
  l1 = 128.839 + n*d - 43.415*sin(tt) - 0.714*sin(3.0*tt) - 0.020*sin(5.0*tt);
  p = 107.0 + 365.560*dt;
 
  // in-orbit-plane position vector
  m = l1 - p; // mean anomaly
  m = m * M_PI / 180.0;
  m = eccanom (m, ecc);  // eccentric anomaly
  rs.assign (a*(cos(m) - ecc), a*sqrt(1.0 - ecc*ecc)*sin(m), 0);
 
  // convert to ecliptic coordinates
  m = p - th1;
  m = m * M_PI / 180.0;
  mt1 = zrot (-m);
  rs = mxvct (mt1, rs);
  mt1 = xrot (-gam);
  rs = mxvct (mt1, rs);
  m = th1*M_PI/180.0 - oe;
  mt1 = zrot (-m);
  rs = mxvct (mt1, rs);
  mt1 = xrot (-ie);
  rs = mxvct (mt1, rs);
  mt1 = zrot (-oe);
  rs = mxvct (mt1, rs);
 
  // convert to ecliptic coordinates of date
  mt1 = pmatecl (-0.500002096, t);
  rs = mxvct (mt1, rs);
 
  return rs;
 }
 
 Vec3 PosSEnceladus (double t)
 {
  /* Ecliptic coordinates (in A.U.) of Enceladus with respect to Saturn 
     for Equinox of Date.
     t is the time in Julian centuries since J2000.
   ===================================================================*/
  const double ie = 28.0771*M_PI/180.0;
  const double oe = 168.8263*M_PI/180.0;
 
  double d, dt, /*tt,*/ a, n, ecc, gam, th1, l1, p, m;
  Vec3 rs;
  Mat3 mt1;
 
  d = 36525.0*t + 40452.0;
  dt = d / 365.25;
  a = 0.00158935;
  n = 262.7319052;
  ecc = 0.00532;
  gam = 0.000628319; //  0.036;
  th1 = 145.0 - 152.7*dt;
  l1 = 200.155 + n*d + 0.256333 * sin((59.4+32.73*dt)*M_PI/180.0)
       + 0.217333 * sin((119.2+93.18*dt)*M_PI/180.0); 
  p = 312.7 + 123.42*dt;
 
  // in-orbit-plane position vector
  m = l1 - p; // mean anomaly
  m = m * M_PI / 180.0;
  m = eccanom (m, ecc);  // eccentric anomaly
  rs.assign (a*(cos(m) - ecc), a*sqrt(1.0 - ecc*ecc)*sin(m), 0);
 
  // convert to ecliptic coordinates
  m = p - th1;
  m = m * M_PI / 180.0;
  mt1 = zrot (-m);
  rs = mxvct (mt1, rs);
  mt1 = xrot (-gam);
  rs = mxvct (mt1, rs);
  m = th1*M_PI/180.0 - oe;
  mt1 = zrot (-m);
  rs = mxvct (mt1, rs);
  mt1 = xrot (-ie);
  rs = mxvct (mt1, rs);
  mt1 = zrot (-oe);
  rs = mxvct (mt1, rs);
 
  // convert to ecliptic coordinates of date
  mt1 = pmatecl (-0.500002096, t);
  rs = mxvct (mt1, rs);
 
  return rs;
 }
  
 Vec3 PosSDione (double t)
 {
  /* Ecliptic coordinates (in A.U.) of Dione with respect to Saturn 
     for Equinox of Date.
     t is the time in Julian centuries since J2000.
   ===================================================================*/
 
  const double ie = 28.0771*M_PI/180.0;
  const double oe = 168.8263*M_PI/180.0;
 
  double d, dt, /*tt,*/ a, n, ecc, gam, th1, l1, p, m;
  Vec3 rs;
  Mat3 mt1;
 
  d = 36525.0*t + 40452.0;
  dt = d / 365.25;
  a = 0.00252413;
  n = 131.534920026;
  ecc = 0.001715;
  gam = 0.0005044; // 0.0289;
  th1 = 228.0 - 30.6197 * dt;
  l1 = 255.1183 + n*d - 0.0146667 * sin((59.4+32.73*dt)*M_PI/180.0)
       - 0.0125 * sin((119.2+93.18*dt)*M_PI/180.0);
  p = 173.6 + 30.8381 * dt;
 
  // in-orbit-plane position vector
  m = l1 - p; // mean anomaly
  m = m * M_PI / 180.0;
  m = eccanom (m, ecc);  // eccentric anomaly
  rs.assign (a*(cos(m) - ecc), a*sqrt(1.0 - ecc*ecc)*sin(m), 0);
 
  // convert to ecliptic coordinates
  m = p - th1;
  m = m * M_PI / 180.0;
  mt1 = zrot (-m);
  rs = mxvct (mt1, rs);
  mt1 = xrot (-gam);
  rs = mxvct (mt1, rs);
  m = th1*M_PI/180.0 - oe;
  mt1 = zrot (-m);
  rs = mxvct (mt1, rs);
  mt1 = xrot (-ie);
  rs = mxvct (mt1, rs);
  mt1 = zrot (-oe);
  rs = mxvct (mt1, rs);
 
  // convert to ecliptic coordinates 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);
 
   }