Chapter 1, Section 1.3, Figure 1.2

Define a function to calculate the radius of an ellipse of semi-major axis a and eccentricity e at an angular coordinate θ using the polar equation of an ellipse (to be derived in Sect.2.3):

r[a_,e_,θ_]:=a(1-e^2)/(1+e Cos[θ])

Now define a function that will draw an ellipse of semi-major axis a and eccentricity e using Mathematica's ParametricPlot command:

ellipse[a_,e_]:=
ParametricPlot[{r[a,e,θ] Cos[θ], r[a,e,θ] Sin[θ]},
               {θ,0,2Pi},
               AspectRatio->Automatic,Axes->None]

Draw an ellipse with a=1 and e=0.5

ell=ellipse[1,0.5];

The foci of the ellipse are at (-1,0) (the empty focus) and (0,0) (the Sun) in the heliocentric polar coordinate system we are using.  Define a graphics object to be used to produce a large point at each focus.

focus1=Graphics[{PointSize[0.05],Point[{-1,0}]}]

focus2=Graphics[{PointSize[0.05],Point[{0,0}]}]

Place the planet on the ellipse at a position of, say, θ=57.78° by finding the corresponding radius and hence the cartesian coordinates x and y.

radius = r[1, 0.5, 57.78°]

{x,y}=radius*{Cos[57.78°],Sin[57.78°]}

Define graphics objects for the planet at this point and for the line joining the Sun and planet.

planet=Graphics[{PointSize[0.02],Point[{x,y}]}]

radial=Graphics[{Line[{{0,0},{x,y}}]}]

Show the ellipse, the foci, the planet and the radial line.

Show[ell,focus1,focus2,planet,radial]

In order to illustrate Kepler's second law we define a function to solve Kepler's equation M=E-eSin[E] in order to find the eccentric anomaly E for given values of the mean anomaly M and the eccentricity e (see Sect.2.4).  Note that we are taking E=M as a first approximation in the FindRoot command.

ksolve[Ma_, e_]:= FindRoot[Ea - e Sin[Ea] == Ma, {Ea,Ma}]

The mean anomaly is a quantity that is a linear function of time so we choose three equal intervals in M to illustrate Kepler's second law: (i) -0.2->0.3, (ii) 2.4->2.9 and (iii) 4.6->5.1 radians.

Mlist={{-0.2,0.3},{2.4,2.9},{4.6,5.1}}

Now solve Kepler's equation to find the corresponding eccentric anomaly E for each of these values. (Note the use of the RuleDelayed (:>) command.)

Elist=
Mlist /. 
{start_, end_} :> {Ea /. ksolve[start,0.5], Ea /. ksolve[end,0.5]}

Calculate the (x,y) coordinates of the points corresponding to each value of the eccentric anomaly using the equations x=Cos[E]-e and y=(Sqrt[1 - (e^2)])Sin[E].

points = Elist /. 
{start_, end_} :> {{Cos[start]-0.5, Sqrt[0.75]*Sin[start]},
                   {Cos[end]-0.5,   Sqrt[0.75]*Sin[end]}}

Create a graphics object with the lines drawn from the Sun to each of the six points along the ellipse.

plines=Graphics[Map[Line[{{0,0},#}]&,Flatten[points,1]]]

Display the ellipse, the Sun and the three ranges of θ covered in equal time intervals (i.e. equal intervals of M).

Show[ell,focus2,plines]