Kepler’s Equation

Introduction

Since at least the time of Hipparchus and Eudoxus the ancient world believed that the Sun, moon, planets, and stars moved around the earth in circular orbits. Aristotle’s Physics put the heavenly bodies on perfect crystal spheres. This theory was further formalized in the second century by Ptolemy in his Almagest which served the basic needs of astronomers for the next 1,500 years. Over that time, Aristotelian physics became an article of faith not to be questioned.

But there were unexplained irregularities that never quite fit the theory. The ancients knew that the seasons were of different lengths. For example, Winter is about 89 days in length while Summer is about 94 days. Why did the Sun sometimes speed up and slow down like that? They also noticed that Mars would move East across the night sky for a few years, slow down, and then reverse course and move West for a few months before looping back to its original course. (What we now call apparent retrograde motion.)

Copernicus took a step forward by putting the Sun at the center of the universe but he still retained and tweaked Ptolemy’s epicycles. Tycho Brahe proposed a hybrid model in which the planets orbited the Sun but the Sun and Moon orbited the Earth. Yet he too retained the circular orbits. In his Epitome of Copernican Astronomy, Johannes Kepler summed up the quandary nicely:

The ancients wished it to be the office of the astronomer to bring forward such causes of this apparent irregularity as would bear witness that the true movement of the planet or spheres is most regular, most equal, and most constant, and also of the most simple figure, that is, exactly circular. And they judged that you should not listen to him who laid down that there was actually any irregularity at all in the real movements of these bodies (Book IV, 571).

So the problem came down to metaphysics. The circle was a symbol of perfection. And the heavens were the perfect creation of God. Who would suggest otherwise and earn the wrath of all Christendom? And yet Kepler ended up doing exactly that, finally proposing the ellipse as the figure the planets describe as they orbit the Sun. Kepler worked out a method (Kepler’s Equation) that determines where an orbiting body like the Earth is in its orbit around the Sun. That is the subject I’ll explore here.

The Problem

The version of Kepler’s Equation (KE) I’ll walk through works for any elliptical orbit in the interval [0, 1). That is to say where 0 <= e < 1 with 0 being a circle and 1 being a parabola. In the figure below we have an elliptical orbit of planet P. The eccentricity of the ellipse is greatly exaggerated for the purposes of illustration; Earth’s orbital eccentricity e is actually about 0.016 making it closer to a circle.

Of the two foci, the Sun S is shown at one focus and the other is empty. Planet P now at time t is moving counterclockwise and went through apsis A (perihelion) at time τ. What we want to know is where is P now at time t? To get there we have to explore three kinds of angles in orbital mechanics: the mean anomaly, the eccentric anomaly, and the true anomaly.

Properties of an Ellipse

Every ellipse has a center, two foci, and a major and minor axis. In the figure of the ellipse below I’ve labeled AB as the major axis and CD the minor axis. I’ve also labeled semi-major axis a and semi-minor axis b.

Earlier I said the eccentricity of the ellipse is going to be in the interval [0, 1). Its eccentricity is related to a and b by:

$$e=\sqrt{1-\left\lgroup b \over a \right\rgroup^2}$$

Mean Anomaly

Kepler’s second law tells us that “planets sweep out equal areas in equal times.” That is to say P does not move uniformly in its orbit. It speeds up at perihelion when it is closest to the gravitational pull of S and slows down when it is furthest away at aphelion. This irregularity is the main reason the problem is hard to solve.

But what if we set aside the elliptical motion for now and created a fictitious planet P’ moving at a constant speed in a uniform circular orbit? It would be easier to calculate that angle and then use it later to solve the real problem. In this diagram P’ moves in a circle around center C while the real planet P moves in an ellipse around focus S:

This mean (average) movement of P’ between any two points in time always sweeps out an equal area for any duration of time. The angle created is called the mean anomaly. In the figure above mean anomaly angle M is the interior angle of the area CP’A. Let T represent the time it takes for P’ to complete one orbit around S. For example, if P is the Earth then T = 365.256 days.  So the radius vector CP’ sweeps counterclockwise from A to 2π (0 to 360 degrees) once in time T. At a given moment t we can determine M by:

$$M = {2 \pi t  \over T}$$

We will need M later to solve KE. First, let’s look at a second angle called the eccentric anomaly that Kepler created as a stepping stone to get to the true anomaly.

Eccentric Anomaly

Here’s Kepler’s trick for creating a new angle he called the eccentric anomaly. As with angle M describe a circle whose diameter is equal to the ellipse’s major axis AB as in the figure below. In other words, the center of circle C is also the center of the ellipse. Draw a line segment perpendicular to AB that runs from R on up through P to the circle at point Q. Angle QCA is eccentric anomaly E. Angle PSA is the true anomaly ν.  Angle E forms a relationship with angle ν as we shall see next.

True Anomaly

Note in the figure below that a is also the radius of the described circle. Angle ν is the true anomaly but I’ve now added the radius vector to the diagram. The radius vector (hypotenuse r in PSR) is the distance center-to-center measured in astronomical units (AU) between the Sun S and the planet P.

Now we see why Kepler described the circle with a diameter equal to the major axis of the ellipse:

$${b \over a} = {PR \over QR}$$

This section goes through the relationship between ν and E. If you’re not up for the trigonometry feel free to skip down to the last function (85). Otherwise following Smart (pp. 112-113) \( PR = r \sin \nu \) and \(QR = CQ \sin E\) = \( a \sin E\) so:

$$r \sin \nu = b \sin E \tag{61}$$

And \(SR = r \cos \nu \) so \(SR = CR – CS = a \cos E – ae\) so:

$$r \cos \nu = a (\cos E – e) \tag{62}$$

$$r = a(1 – e \cos E)$$

$$2r \sin^2 {\nu \over 2} = r (1 – e \cos E)$$

$$=a(1 – e \cos E)-a(\cos E – e) \tag{63}$$

$$2r \sin^2 = a (1 + e) (1 – \cos E) \tag{64}$$

$$2r \cos^2 = a (1 – e) (1 + \cos E) \tag{65}$$

Dividing (65) into (64):

$$ \tan^2 {\nu \over 2} = {1+e \over 1-e} \cdot {1- \cos E \over 1 + \cos E}$$

$$ \tan {\nu \over 2} = \sqrt{1+e \over 1-e} \tan {E \over 2}$$

Inverting so we can solve for ν:

$$\nu = 2 \arctan {\left\lgroup \sqrt{1+e \over 1-e} \tan {E \over 2}\right\rgroup}$$

Smart (p. 119) applies the formula for the logarithmic series to derive:

$$\nu = E + (e + 0.25 e^3) \sin E + 0.25e^2 \sin 2E + 0.0833e^3 \sin 3E \tag{85}$$

Function (85) can easily be translated to programming code and solved if we know e and E. That’s where KE comes in.

Kepler’s Equation

Colwell (p. 3-4) walks through the derivation of KE. Area PSA = area PSR + area PRA. Then:

$$= {1 \over 2} a (\cos E – e)(r \sin \nu) + {b \over a} Area \ QRA$$

$$= {1 \over 2} ab (\cos E – e)( \sin E) + {b \over a} \left\lgroup {1 \over 2} a^2 E = {1 \over 2} a^2 \sin E \cos E \right\rgroup$$

$$= {1 \over 2} a^2 \sqrt {1 – e^2} (E – e \sin E)$$

Finally:

$$M = E – e \sin E$$

Which I will put in the form:

$$E = M + e \sin E$$

Where M is the mean anomaly angle, e is the ellipse’s eccentricity, and E is the eccentric anomaly angle. Because of \(\sin E\) this is a transcendental equation. Therefore, we have to use an iteration method (or a series expansion) to solve for E. I am going to show two iteration methods, the first of which would have been familiar to Kepler. Once we have E we can get coordinates P(x,y) from:

$$x = a (\cos E – e)$$

$$y = b \sin E$$

Solving Kepler’s Equation

Now I’m going to convert the functions needed above into the C# programming language since I’m fluent in that. As with my previous posts on the subject I’m going to rely heavily on Meeus for the core algorithms as well as bit and pieces of functions from various Wiki articles. To keep it simple I’m just going to calculate the Earth’s orbit; however, the public source code available in my GitHub repository is structured such that you could add Mars or any other planet to the solution.

As detailed in a previous post I have a Moment struct that accepts any DateTime value and can return the Julian Day (JD) or the Time T in Julian centuries of 36525 ephemeris days from the J2000.0 epoch.

A second struct is named Earth which accepts a Moment value and calculates the necessary orbital elements for the Earth, its mean anomaly M and eccentricity e in particular. Following Meeus, to find M for a given moment we need the four coefficients for both Earth’s mean longitude L and the perihelion’s mean longitude P which I’ve implemented as a simple Dictionary:

/// <summary>
/// Orbital element mean longitude L coefficients from Meeus Table 31.A
/// </summary>
public Dictionary<string, double> L => new Dictionary<string, double>
{
    { "a0", 100.466457 },
    { "a1", 36000.7698278 },
    { "a2", 0.00030322 },
    { "a3", 0.000000020 }
};

/// <summary>
/// Orbital element perihilion longitude π (ϖ) coefficients from Meeus Table 31.A 
/// </summary>
public Dictionary<string, double> P => new Dictionary<string, double>
{
    { "a0", 102.937348 },
    { "a1", 1.7195366 },
    { "a2", 0.00045688 },
    { "a3", -0.000000018 }
};

/// <summary>
/// Mean anomaly of the Earth from Meeus p.210 where M = L - π or 
/// Colwell p. 3 has M = 2πt / T. Note that Colwell's use of π means 
/// literally 3.1415... and NOT the longitude of the perihilion as 
/// Meeus uses it.
/// </summary>
public double MeanAnomaly
{
    get
    {
        var T = Moment.TimeT;

        var L0 = L["a0"] +
            (L["a1"] * T) +
            (L["a2"] * Math.Pow(T, 2)) +
            (L["a3"] * Math.Pow(T, 3));

        var P0 = P["a0"] + 
            (P["a1"] * T) + 
            (P["a2"] * Math.Pow(T, 2)) + 
            (P["a3"] * Math.Pow(T, 3));

        return L0 - P0;
    }
}

I use the same implementation for eccentricity e:

/// <summary>
/// Orbital element eccentricity e coefficients from Meeus Table 31.A 
/// </summary>
public Dictionary<string, double> e => new Dictionary<string, double>
{
    { "a0", 0.01670863 },
    { "a1", -0.000042037 },
    { "a2", -0.0000001267 },
    { "a3", -0.00000000014 }
};

/// <summary>
/// Eccentricity of the Earth's orbit at a given moment in time. 
/// </summary>
/// <param name="m">Moment</param>
/// <returns></returns>
public double Eccentricity
{
    get
    {
        var T = Moment.TimeT;
        return e["a0"] + 
            (e["a1"] * T) + 
            (e["a2"] * (T * T)) +
            (e["a3"] * (T * T * T));
    }
}

Now that I have M and e I can solve KE. Meeus provides four methods. I implemented the first two and put some instrumentation around them to get a sense of their performance. The first method is a simple fixed-point iteration, which is what Kepler himself did to solve the equation. Let \(E_0 = M\) be the default value for E in the first calculation. Then take the output and pass it back into the equation a second time for greater accuracy:

$$E_1 = M + e \sin E_0$$

$$E_2 = M + e \sin E_1$$

$$E_3 = M + e \sin E_2$$

and so on until the desired accuracy is reached. Kepler stopped at \(E_2\) because that was good enough in his day. But in my implementation of this first method I found it needed 120 iterations (taking 0.0928 ms) to achieve the desired accuracy:

/// <summary>
/// Meeus First Method (p. 196)
/// </summary>
/// <returns></returns>
private static double FirstMethod(double e, double M)
{
    double e0 = e.ToDegrees();
    double E = M;

    var stopwatch = new Stopwatch();
    stopwatch.Start();

    int n = 120;
    for(int i = 1; i < n; i++)
    {
        E = M + e0 * Math.Sin(E.ToRadians());
    }

    stopwatch.Stop();
    var ts = stopwatch.Elapsed;
    Console.WriteLine(n + " iterations " + "(" + ts.TotalMilliseconds + "ms)");
    return E;
}

The second method is much faster and produces the desired level of accuracy after only 5 iterations and 0.0059 ms:

/// <summary>
/// Meeus Second Method (p. 199)
/// </summary>
/// <returns></returns>
private static double SecondMethod(double e, double M)
{
    double e0 = e.ToDegrees();
    double E0 = M;
    double E1 = 0;

    var stopwatch = new Stopwatch();
    stopwatch.Start();

    int n = 5;
    for(int i = 1; i < n; i++)
    {
        E1 = E0 + (M + e0 * Math.Sin(E0.ToRadians()) - E0) / 
                  (1 - e * Math.Cos(E0.ToRadians()));
        E0 = E1;
    }

    stopwatch.Stop();
    var ts = stopwatch.Elapsed;
    Console.WriteLine(n + " iterations:  " + "(" + ts.TotalMilliseconds + "ms)");
    return E1;
}

Given e and E we can now calculate the true anomaly ν using equation (85) above:

/// <summary>
/// True anomaly expressed as a series in terms of e and E
/// via Smart pp. 118-119
/// </summary>
/// <param name="e"></param>
/// <param name="E"></param>
/// <returns></returns>
private static double GetTrueAnomaly2(double e, double E)
{
    return E + (e + 0.25 * Math.Pow(e, 3)) * 
        Math.Sin(E.ToRadians()) + (0.25 * Math.Pow(e, 2) * 
        Math.Sin(2 * E.ToRadians())) + (0.083333 * Math.Pow(e, 3) * 
        Math.Sin(3 * E.ToRadians()));
}

And of course now that we have ν we can get the radius vector r:

/// <summary>
/// Radius vector from Meeus (25.5) on p. 164 and also see section
/// on Wiki entry for "True Anomaly" entitled "Radius from true anomaly"
/// </summary>
/// <param name="e">eccentricity</param>
/// <param name="v">true anomaly</param>
/// <returns></returns>
private static double GetRadiusVector(double e, double v)
{
    return 1.000001018 * (1 - Math.Pow(e, 2)) / 1 + e * Math.Cos(v.ToRadians());
}

Output

I ran the program for 9 PM this evening as of this writing (2019-04-07 21:00 UTC):

Now to check this output against JPL’s Horizons system:

Here’s the comparison between the three highlighted values against my output:

	Eccentricity e	Mean anomaly M (deg)	True anomaly v (deg)
My Output	0.01670	92.58	93.55
Horizons	0.01648	90.13	92.02

I’m going to assume that JPL is the gold standard! But because I’m following Meeus (Chap. 31) my value for M accords with the VSOP87 theory while JPL uses ICRF. A quick check with Mathematica shows M = 92.15. In any case, this isn’t bad for “armchair” computing on a personal laptop.

Sources

I am a software developer by trade and not a mathematician. So I was absolutely dependent on all three of these sources for my understanding of KE. Most of the trigonometry above is copied exactly from Smart and I retained his formula numbers for easy reference. If there is no formula number then you can certainly blame me rather than Smart for any errors. My task was largely to organize and make sense of the material for my own benefit and then use that understanding to implement the software solution.

Colwell, Peter. Solving Kepler’s Equation Over Three Centuries. Willmann-Bell. 1993.

Meeus, Jean. Astronomical Algorithms. 2nd Ed. Willmann-Bell. 1998.

Smart, W.M. Textbook on Spherical Astronomy. 6th Ed. Cambridge University Press. 1977.

Shout Outs

A thank you to Ketil Wright for pointing out equal sign typos in equations (63) and (66) and a place where it should have read “focus” (singular) instead of “foci” (plural).