Chapter 9, Section 9.9, Figures 9.33–34

Here we make use of many of the functions defined in, for example, Fig.2.15.nb.  First we need to define 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]];

Calculate the Julian date corresponding to 1st September 1989 when Pluto was at perihelion.

date = JulianDate[{1989, 9, 1, 0.}]

Next we define a function PrintPlanetXYZ 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.

PrintPlanetXYZ[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}]
];

Here we define a function PlanetXYZ that is identical to PrintPlanetXYZ except that it does not print out any of the results.  Again it uses the data in Tables A.2 and A.3 to derive the (X,Y,Z) coordinates of planet number id at Julian Date jd.  Note that id=1,2,...,9 corresponds to Mercury, Venus, ..., Pluto.

PlanetXYZ[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, p1, p2, p3, p, xx, yy, zz},
t     = (jd - j2000)/36525;
aa    = a[[id]] + δa[[id]]*t;
ee    = e[[id]] + δe[[id]]*t;
ii    = Mod[i[[id]]    + δi[[id]]*t/3600, 360];
nnode = Mod[Ω[[id]] + δΩ[[id]]*t/3600, 360];
If[nnode < 0, nnode = nnode + 360];
llp   = Mod[ϖ[[id]]   + δϖ[[id]]*t/3600, 360];
If[llp < 0, llp = llp + 360];
ω = llp - nnode;
llam  = Mod[λ[[id]]  + 
            (δλ[[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 °];
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[ω °],    -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];
p.{x, y, 0}
];

Now define a function PlanetElements that uses the data in Tables A.2 and A.3 to return the semi-major axis, eccentricity, inclination, argument of perihelion and longitude of ascending node for planet number id at Julian date jd.

PlanetElements[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, p1, p2, p3, p, xx, yy, zz},
t     = (jd - j2000)/36525;
aa    = a[[id]] + δa[[id]]*t;
ee    = e[[id]] + δe[[id]]*t;
ii    = Mod[i[[id]]    + δi[[id]]*t/3600, 360];
nnode = Mod[Ω[[id]] + δΩ[[id]]*t/3600, 360];
If[nnode < 0, nnode = nnode + 360];
llp   = Mod[ϖ[[id]]   + δϖ[[id]]*t/3600, 360];
If[llp < 0, llp = llp + 360];
ω = llp - nnode;
{aa,ee,ii,ω,nnode}
]

Now define a function OrbitXYZ that returns a list of npts values of (X,Y,Z) around the orbit of planet id at Julian date jd.  Note that this function makes use of PlanetElements to return the orbital elements at jd.

OrbitXYZ[id_, jd_, npts_] :=
Module[{elem, a, e, ii, ω, Ω, p1, p2, p3, p, x, y, z,
        f, r, cf, sf},
elem = PlanetElements[id, jd];
a   = elem[[1]];
e   = elem[[2]];
ii  = elem[[3]];
ω   = elem[[4]];
Ω   = elem[[5]];
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[Ω °],   -Sin[Ω °],    0             },
      {Sin[Ω °],    Cos[Ω °],    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}]]
]

In order to make drawings of the planetary orbits we need to be able to take the (X,Y) coordinates from the list obtained from OrbitXYZ and make them into a Line of joined up points.  Hence we define a function OrbitXY.

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

Now we can define a function PlotPlanets that will draw the selected planetary orbits at the given Julian date jd.

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

Finally we define a function PlotNeptunePluto that uses the previous functions to plot the orbits of Neptune and Pluto at the specified date.

PlotNeptunePluto[jd_, npts_] :=
Module[{outer = {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]
]

Now do the plot.

nepplu= PlotNeptunePluto[date, 200]

Print out the elements for Neptune at that date.

PrintPlanetXYZ[8, date]

Print out the elements for Pluto at that date.

PrintPlanetXYZ[9, date]

Use the value of Neptune's longitude of ascending node to create a graphics object for an approximate line of nodes.

nodeneptune = 
Graphics[Line[
{31*{Cos[131.726°], Sin[131.726°]},
 31*{Cos[(131.726+180)°], Sin[(131.726+180)°]}}]]

Use the value of Pluto's longitude of ascending node to create a graphics object for an approximate line of nodes.

nodepluto = 
Graphics[Line[
{43*{Cos[110.305°], Sin[110.305°]},
 43*{Cos[(110.305+180)°], Sin[(110.305+180)°]}}]]

Display the orbits, the positions and the lines of nodes (see Fig.9.33).

Show[nepplu, nodeneptune, nodepluto]

Calculate the Julian date at the start of 1989, our reference year.

y1989 = JulianDate[{1989, 1, 1, 12.}]

Define a simple function, Separation, that calculates the distance between points defined by two position vectors, v1 and v2.

Separation[v1_, v2_] :=
Sqrt[Apply[Plus, (v1-v2)^2]]

Calculate the separation of Neptune and Pluto for a period of 1000 years centred on 1989.

sep =
Table[{1989+i, 
       Separation[PlanetXYZ[8, y1989+i*365.25],
                  PlanetXYZ[9, y1989+i*365.25]]},
      {i,-500,500,2}];

Plot the results (see Fig.9.34).

ListPlot[sep, PlotJoined->True, Frame->True, 
PlotRange->{{1489, 2489}, {0,80}}, AxesOrigin->{1489,0}]