Reference: Carroll, Bradley W. & Ostlie, Dale A. (2007), *An Introduction to Modern Astrophysics*, 2nd Edition; Pearson Education – Chapter 2, Problem 2.15.

We can use Kepler’s laws to generate numerical solutions for orbits on a computer. Carroll & Ostlie provide a program called *Orbit* (in Fortran or C++) that produces textual output that can be input to a spreadsheet or graphing program. However, this is a bit primitive, so I decided to implement a similar routine in Maple, since Maple has built-in plotting features.

The idea is to simulate a complete orbit of a planet by calculating its distance from the focus of the ellipse (essentially, the distance from the Sun) as a function of the angle from perihelion. The equation of an ellipse is:

Given the planet’s period we can calculate its semimajor axis from Kepler’s third law:

so the equation 1 becomes

To use this formula in a numerical solution, we need to know how much changes for a given time increment . We can use Kepler’s second law in the form

which gives the rate at which area is swept out in terms of the constant total angular momentum and reduced mass . The angular momentum is

The area increment is given in terms of by

so we get

The numerical solution divides the period up into equal time intervals , and starts with at perihelion and . We then use the value of to calculate , add this to the current value of and calculate the next value of from 3. We then use the new value of to get the next and so on until we’ve covered a complete period.

The Maple code for doing this is:

with(plots): G := 0.6673e-10; AU := 0.14959787066e12; M__sun := 0.19891e31; yr := 31558145.0; rad2deg := 180/Pi; secsYear := 365.25*(3600*24); orbit := proc (M__strsun, a__AU, e, n) local M__star, a, P, dt, t, theta, LoM, r, i, dtheta; M__star := M__strsun*M__sun; a := a__AU*AU; P := sqrt(4*Pi^2*a^3/(G*M__star)); dt := P/(n-1); t := Array(0 .. n, (i) -> i*dt); r := Array(0 .. n, datatype = float); theta := Array(0 .. n, datatype = float); theta[0] := 0.; LoM := sqrt(G*M__star*a*(1-e^2)); for i from 0 to n-1 do r[i] := a*(1-e^2)/(1+e*cos(theta[i])); dtheta := LoM*dt/r[i]^2; theta[i+1] := evalf(theta[i]+dtheta); if 0 < i and 1.0 < r[i]/AU and r[i-1]/AU < 1.0 then print("Passes 1 AU at t = ", evalf(t[i]/secsYear)) end if end do; r[n] := a*(1-e^2)/(1+e*cos(theta[n])); theta := theta*rad2deg; r := r/AU; print(polarplot(r, theta, angularunit = degrees)); print(plot(t/secsYear, r, labels = ["t (years)", "r (AU)"])) end proc;

The code pretty much just implements the algorithm given above. The Maple procedure *orbit* on line 8 takes as its arguments the mass of the central star in solar masses, the semimajor axis of the planet in AU, the eccentricity and the number of time steps . It then converts these quantities into SI units using the conversion factors given at the start, creates arrays for the time , the radius and the angle , calculates (as LoM), and then enters a for loop to calculate and for each time increment. The ‘if’ statement finds the time at which the planet crosses from to and prints this out.

After the loop, we calculate the final value of , convert and to degrees and AU, respectively, and print out a couple of plots. The first plot is a polar plot of as a function of , so it shows the elliptical orbit. The second plot graphs as a function of (the latter in years) for one complete period.

For Halley’s comet, and from which we find . If we input these values (along with M\_\_strsun = 1 since the central star is the Sun), and choose time increments, we get the time at which the comet first crosses after perihelion as or around 39 days. The plots are as follows: