Chapter 2, Section 2.8, Figure 2.15

There are a number of different functions required in order to produce the data quoted in th etext and then Fig.2.15.  First we must be able to calculate the Julian Date corresponding to any calendar date.  To do this we need to define a functions Trunc, Int and Frac that can truncate, take the integer and fractional parts of positive and negative real numbers.

Trunc[x_] := If[x > 0, Floor[x], Ceiling[x]]
Int[x_] := Floor[x]
Frac[x_] := x - Trunc[x]

Now define the function JulianDate to calculate the Julian Date for a given calendar year, month, day and UT.  This uses the algorithm given in Appendix A.3.

JulianDate[{year_Integer, month_Integer, 
            day_Integer, ut_}] :=
Module[{y, m, b},
  If[month <= 2,
     y = year-1; m = month+12,
     y = year;   m = month    ];
  If[year < 1582, b = -2];
  If[year > 1582, b = Int[y/400] - Int[y/100]];
  If[year == 1582,
     If[month < 10, b = -2];
     If[month > 10, b = Int[y/400] - Int[y/100]];
     If[month == 10,
        If[day <= 4,  b = -2];
        If[day >= 15, b = Int[y/400] - Int[y/100]];
        If[day > 4 && day < 15, Print["Missing day"]]
       ]
    ];
  N[If[y > 0, Int[365.25*y], Trunc[365.25*y-0.75]] + 
    Int[30.6001*(m+1)] + b + 1720996.5 + day + ut/24,
    12]];

The positions and orbits of the planets are required for 1993 September 25 at 5.32pm British Summer Time.  The latter is one hour ahead of Greenwich Mean Time (approximately equivalent to UT for our purposes) and so we need to make allowances.

date = JulianDate[{1993, 9, 25, 16+32./60.}]

NumberForm[date,12]

Next we define a function PrintPlanetElements that will use the data in Tables A.2 and A.3 to derive and print out details of the calculation of the coordinates and orbital elements of planet number id at Julian Date jd.  Note that id=1,2,...,9 corresponds to Mercury, Venus, ..., Pluto.

PrintPlanetElements[id_,jd_] :=
Module[{ a = {    0.38709893,     0.72333199,     1.00000011, 
                  1.52366231,     5.20336301,     9.53707032,
                 19.19126393,    30.06896348,    39.48168677},
        δa = {    0.00000066,     0.00000092,    -0.00000005,
                 -0.00007221,     0.00060737,    -0.00301530,
                  0.00152025,    -0.00125196,    -0.00076912},
         e = {    0.20563069,     0.00677323,     0.01671022,
                  0.09341233,     0.04839266,     0.05415060,
                  0.04716771,     0.00858587,     0.24880766},
        δe = {    0.00002527,    -0.00004938,    -0.00003804,
                  0.00011902,    -0.00012880,    -0.00036762,
                 -0.00019150,     0.00002514,     0.00006465},
         i = {    7.00487,        3.39471,        0.00005,
                  1.85061,        1.30530,        2.48446,
                  0.76986,        1.76917,       17.14175   },
        δi = {  -23.51,          -2.86,         -46.94,
                -25.47,          -4.15,           6.11,
                 -2.09,          -3.64,          11.07      },
         Ω = {   48.33167,       76.68069,      -11.26064,
                 49.57854,      100.55615,      113.71504,
                 74.22988,      131.72169,      110.30347   },
        δΩ = { -446.30,        -996.89,      -18228.25,
              -1020.19,        1217.17,       -1591.05,
               1681.40,        -151.25,         -37.33      },
         ϖ = {   77.45645,      131.53298,      102.94719,
                336.04084,       14.75385,       92.43194,
                170.96424,       44.97135,      224.06676   },
        δϖ = {  573.57,        -108.80,        1198.28,
               1560.78,         839.93,       -1948.89,
               1312.56,        -844.43,        -132.25      },
         λ = {  252.25084,      181.97973,      100.46435,
                355.45332,       34.40438,       49.94432,
                313.23218,      304.88003,      238.92881   },
        δλ = {261628.29,      712136.06,     1293740.63,
              217103.78,      557078.35,      513052.95,
              246547.79,      786449.21,      522747.90     },
       rev = {  415,            162,             99,
                 53,              8,              3,
                  1,              0,              0         },
       j2000 = 2451545.0, t, aa, ee, ii, nnode, llp, om, llam,
       ma, x, y, ea, eca, deca, p1, p2, p3, p, r},
t     = (jd - j2000)/36525;
aa    = a[[id]] + δa[[id]]*t;
Print["Semimajor axis (AU) : ",aa];
ee    = e[[id]] + δe[[id]]*t;
Print["Eccentricity        : ",ee];
ii    = Mod[i[[id]]    + δi[[id]]*t/3600, 360];
Print["Inclination         : ",ii];
nnode = Mod[Ω[[id]] + δΩ[[id]]*t/3600, 360];
If[nnode < 0, nnode = nnode + 360];
Print["Long. of Asc. Node  : ",nnode];
llp   = Mod[ϖ[[id]]   + δϖ[[id]]*t/3600, 360];
If[llp < 0, llp = llp + 360];
Print["Long. of Perihelion : ",llp];
ω = llp - nnode;
llam  = Mod[λ[[id]]  + 
            (δλ[[id]]/3600 + rev[[id]]*360)*t, 360];
If[llam < 0, llam = llam + 360];
Print["Mean Longitude      : ",llam];
ma = Mod[llam - llp, 360];
If[ma < 0, ma = ma + 360];
Print["Mean Anomaly        : ",ma];
rma = N[ma °];
eca = ea /. FindRoot[ea - ee Sin[ea] == rma, {ea,rma}];
deca = N[eca/°];
If[deca < 0, deca = deca + 360];
Print["Eccentric Anomaly   : ",deca];
x = aa*(Cos[eca] - ee); 
y = aa*Sqrt[1-ee^2]*Sin[eca];
r = Sqrt[x^2 + y^2];
Print["Radius (AU)         : ",r];
Print["{x, y} coordinates  : ",{x,y}];
p1 = {{Cos[ω °],    -Sin[ω °],     0             },
      {Sin[ω °],     Cos[ω °],     0             },
      {0,                  0,                  1             }};
p2 = {{1,                  0,                  0             },
      {0,                  Cos[ii °],    -Sin[ii °]},
      {0,                  Sin[ii °],     Cos[ii °]}};
p3 = {{Cos[nnode °], -Sin[nnode °],  0             },
      {Sin[nnode °],  Cos[nnode °],  0             },
      {0,                  0,                  1             }};
p = N[p3.p2.p1];
Print["Rotation Matrix P1  : ",MatrixForm[N[p1]]];
Print["Rotation Matrix P2  : ",MatrixForm[N[p2]]];
Print["Rotation Matrix P3  : ",MatrixForm[N[p3]]];
Print["Rotation Matrix P   : ",MatrixForm[p]];
Print["{X, Y, Z} coords.   : ",p.{x, y, 0}]
];

PrintPlanetElements[5,date]

PlanetXYZ[id_, jd_] :=
Module[{ a = {    0.38709893,     0.72333199,     1.00000011, 
                  1.52366231,     5.20336301,     9.53707032,
                 19.19126393,    30.06896348,    39.48168677},
        da = {    0.00000066,     0.00000092,    -0.00000005,
                 -0.00007221,     0.00060737,    -0.00301530,
                  0.00152025,    -0.00125196,    -0.00076912},
         e = {    0.20563069,     0.00677323,     0.01671022,
                  0.09341233,     0.04839266,     0.05415060,
                  0.04716771,     0.00858587,     0.24880766},
        de = {    0.00002527,    -0.00004938,    -0.00003804,
                  0.00011902,    -0.00012880,    -0.00036762,
                 -0.00019150,     0.00002514,     0.00006465},
         i = {    7.00487,        3.39471,        0.00005,
                  1.85061,        1.30530,        2.48446,
                  0.76986,        1.76917,       17.14175   },
        di = {  -23.51,          -2.86,         -46.94,
                -25.47,          -4.15,           6.11,
                 -2.09,          -3.64,          11.07      },
      node = {   48.33167,       76.68069,      -11.26064,
                 49.57854,      100.55615,      113.71504,
                 74.22988,      131.72169,      110.30347   },
     dnode = { -446.30,        -996.89,      -18228.25,
              -1020.19,        1217.17,       -1591.05,
               1681.40,        -151.25,         -37.33      },
        lp = {   77.45645,      131.53298,      102.94719,
                336.04084,       14.75385,       92.43194,
                170.96424,       44.97135,      224.06676   },
       dlp = {  573.57,        -108.80,        1198.28,
               1560.78,         839.93,       -1948.89,
               1312.56,        -844.43,        -132.25      },
       lam = {  252.25084,      181.97973,      100.46435,
                355.45332,       34.40438,       49.94432,
                313.23218,      304.88003,      238.92881   },
      dlam = {261628.29,      712136.06,     1293740.63,
              217103.78,      557078.35,      513052.95,
              246547.79,      786449.21,      522747.90     },
       rev = {  415,            162,             99,
                 53,              8,              3,
                  1,              0,              0         },
       j2000 = 2451545.0, t, aa, ee, ii, nnode, llp, om, llam,
       ma, x, y, ea, eca, p1, p2, p3, p, xx, yy, zz},
t     = (jd - j2000)/36525;
aa    = a[[id]] + da[[id]]*t;
ee    = e[[id]] + de[[id]]*t;
ii    = Mod[i[[id]]    + di[[id]]*t/3600, 360];
nnode = Mod[node[[id]] + dnode[[id]]*t/3600, 360];
If[nnode < 0, nnode = nnode + 360];
llp   = Mod[lp[[id]]   + dlp[[id]]*t/3600, 360];
If[llp < 0, llp = llp + 360];
om = llp - nnode;
llam  = Mod[lam[[id]]  + 
            (dlam[[id]]/3600 + rev[[id]]*360)*t, 360];
If[llam < 0, llam = llam + 360];
ma = Mod[llam - llp, 360];
If[ma < 0, ma = ma + 360];
rma = N[ma Degree];
eca = ea /. FindRoot[ea - ee Sin[ea] == rma, {ea,rma}];
x=aa*(Cos[eca] - ee); 
y=aa*Sqrt[1-ee^2]*Sin[eca];
p1 = {{Cos[om Degree],    -Sin[om Degree],     0             },
      {Sin[om Degree],     Cos[om Degree],     0             },
      {0,                  0,                  1             }};
p2 = {{1,                  0,                  0             },
      {0,                  Cos[ii Degree],    -Sin[ii Degree]},
      {0,                  Sin[ii Degree],     Cos[ii Degree]}};
p3 = {{Cos[nnode Degree], -Sin[nnode Degree],  0             },
      {Sin[nnode Degree],  Cos[nnode Degree],  0             },
      {0,                  0,                  1             }};
p = N[p3.p2.p1];
p.{x, y, 0}
];

PlanetElements[id_, jd_] :=
Module[{ a = {    0.38709893,     0.72333199,     1.00000011, 
                  1.52366231,     5.20336301,     9.53707032,
                 19.19126393,    30.06896348,    39.48168677},
        da = {    0.00000066,     0.00000092,    -0.00000005,
                 -0.00007221,     0.00060737,    -0.00301530,
                  0.00152025,    -0.00125196,    -0.00076912},
         e = {    0.20563069,     0.00677323,     0.01671022,
                  0.09341233,     0.04839266,     0.05415060,
                  0.04716771,     0.00858587,     0.24880766},
        de = {    0.00002527,    -0.00004938,    -0.00003804,
                  0.00011902,    -0.00012880,    -0.00036762,
                 -0.00019150,     0.00002514,     0.00006465},
         i = {    7.00487,        3.39471,        0.00005,
                  1.85061,        1.30530,        2.48446,
                  0.76986,        1.76917,       17.14175   },
        di = {  -23.51,          -2.86,         -46.94,
                -25.47,          -4.15,           6.11,
                 -2.09,          -3.64,          11.07      },
      node = {   48.33167,       76.68069,      -11.26064,
                 49.57854,      100.55615,      113.71504,
                 74.22988,      131.72169,      110.30347   },
     dnode = { -446.30,        -996.89,      -18228.25,
              -1020.19,        1217.17,       -1591.05,
               1681.40,        -151.25,         -37.33      },
        lp = {   77.45645,      131.53298,      102.94719,
                336.04084,       14.75385,       92.43194,
                170.96424,       44.97135,      224.06676   },
       dlp = {  573.57,        -108.80,        1198.28,
               1560.78,         839.93,       -1948.89,
               1312.56,        -844.43,        -132.25      },
       lam = {  252.25084,      181.97973,      100.46435,
                355.45332,       34.40438,       49.94432,
                313.23218,      304.88003,      238.92881   },
      dlam = {261628.29,      712136.06,     1293740.63,
              217103.78,      557078.35,      513052.95,
              246547.79,      786449.21,      522747.90     },
       rev = {  415,            162,             99,
                 53,              8,              3,
                  1,              0,              0         },
       j2000 = 2451545.0, t, aa, ee, ii, nnode, llp, om, llam,
       ma, x, y, ea, eca, p1, p2, p3, p, xx, yy, zz},
t     = (jd - j2000)/36525;
aa    = a[[id]] + da[[id]]*t;
ee    = e[[id]] + de[[id]]*t;
ii    = Mod[i[[id]]    + di[[id]]*t/3600, 360];
nnode = Mod[node[[id]] + dnode[[id]]*t/3600, 360];
If[nnode < 0, nnode = nnode + 360];
llp   = Mod[lp[[id]]   + dlp[[id]]*t/3600, 360];
If[llp < 0, llp = llp + 360];
om = llp - nnode;
{aa,ee,ii,om,nnode}
]

OrbitXYZ[id_, jd_, npts_] :=
Module[{elem, a, e, ii, om, lan, p1, p2, p3, p, x, y, z,
        f, r, cf, sf},
elem = PlanetElements[id, jd];
a   = elem[[1]];
e   = elem[[2]];
ii  = elem[[3]];
om  = elem[[4]];
lan = elem[[5]];
p1 = {{Cos[om Degree],    -Sin[om Degree],     0             },
      {Sin[om Degree],     Cos[om Degree],     0             },
      {0,                  0,                  1             }};
p2 = {{1,                  0,                  0             },
      {0,                  Cos[ii Degree],    -Sin[ii Degree]},
      {0,                  Sin[ii Degree],     Cos[ii Degree]}};
p3 = {{Cos[lan Degree],   -Sin[lan Degree],    0             },
      {Sin[lan Degree],    Cos[lan Degree],    0             },
      {0,                  0,                  1             }};
p = N[p3.p2.p1];
f = Table[N[t], {t,0,2Pi,2Pi/npts}];
cf = Cos[f];
sf = Sin[f];
ae = a*(1-e*e);
r = Map[ae/(1+e*Cos[#])&, f];
x = r*cf;
y = r*sf;
z = Table[0,{npts+1}];
Map[(p.#)&, Transpose[{x,y,z}]]
]

OrbitXY[id_, jd_, npts_] :=
  Line[Transpose[Take[Transpose[OrbitXYZ[id, jd, npts]],2]]];

PlotPlanets[idl_List, jd_, npts_] :=
Module[{},
  sun = Point[{0,0}];
  orbits = Map[OrbitXY[#, jd, npts]&, idl];
  positions = 
   Map[{PointSize[0.03],
        Point[Take[PlanetXYZ[#, jd], 2]]}&, idl];
  Show[Graphics[{sun, orbits, positions}],AspectRatio->1]
]

PlotInnerPlanets[jd_, npts_] :=
Module[{inner = {1,2,3,4}},
  sun = Point[{0,0}];
  orbits = Map[OrbitXY[#, jd, npts]&, inner];
  positions = 
   Map[{PointSize[0.03],
        Point[Take[PlanetXYZ[#, jd], 2]]}&, inner];
  Show[Graphics[{sun, orbits, positions}],AspectRatio->1]
]

PlotOuterPlanets[jd_, npts_] :=
Module[{outer = {5,6,7,8,9}},
  sun = Point[{0,0}];
  orbits = Map[OrbitXY[#, jd, npts]&, outer];
  positions = 
   Map[{PointSize[0.03],
        Point[Take[PlanetXYZ[#, jd], 2]]}&, outer];
  Show[Graphics[{sun, orbits, positions}],AspectRatio->1]
]

PlotInnerPlanets[date, 200]

PlotOuterPlanets[date, 200]