Kepler’s Equation

What is it?  That was the last question in my quiz.  Kepler’s three Laws are often mentioned, at least in astronomy textbooks; his Equation is more of a mystery, with a tang to it.

Here is a sample orbit.  I’ll have to keep referring to this diagram, which took me some days to program and more to improve, and is the best diagram of the Keplerian orbit you’ll ever see.  Enlarge it any way you can to see it better.

(The orbit of my fictitious planet has a period of 5 years, and an eccentricity of 0.6 – that is the measure of how far it is from circular.  It has an inclination to the ecliptic plane, so that it has ascending and descending nodes, but that makes no difference to the shape of the ellipse in the diagram, which is from a perpendicular viewpoint.)

Kepler’s three laws govern the way a planet revolves around the Sun (or any other small body around a larger).  The first law describes the orbit’s shape (an ellipse, with the Sun at one of the two foci), and the third law is about the orbit’s size (tying size in time to size in space – “period squared equals average distance cubed”).

It’s left to the second law to tell us where the planet is: its position on its orbit at a certain time.

This would be simple if the orbit were a circle, which it never exactly is.  The planet moves faster along the inner part of its orbit, nearer to the Sun.  And Kepler’s second law is the quantification of this that he discovered: the “radius vector” – the line from the Sun to the planet – “sweeps out” equal areas in equal spans of time.  Or, as we might put it, when the radius vector is shorter the arc traveled by the planet is longer.

The radius vector is the red line from Sun to planet.  The light blue sectors are spaces swept by the radius vector in various spans of 60 days, and they are equal in area.

That’s a description; it doesn’t lead immediately to a calculation (as the third law does).  During his ten years of intense thought leading up to the publication of his Astronomia Nova in 1609, Johannes Kepler managed to translate this description into an actual procedure for finding the planet’s position.  There are several ways, but what he came up with looks so paradoxical that, if Kepler announced it to you, you might take it for a joke, a magical spell – “Abracadabra!” – rather than a serious piece of mathematics.

We measure along an orbit from the perihelion, the point nearest to the Sun.  So, supposing we know the date when the planet is at the perihelion, what we want for any other date is how far around the orbit it has reached.  This angle is called the “anomaly” (a special use of this word, whose general meaning is “a deviation from the usual”).

If the orbit were a circle, the anomaly would be simple: 90° when the planet is a quarter of the way around, 180° at the aphelion, and so on.  This is called the “mean anomaly.”  It’s found easily: if the perihelion date is 2000 Jan. 1 and the period is 10 years, the mean anomaly at 2001 Jan. 1 is 36°.

But that’s not where the planet is: being not far from perihelion, it is traveling faster, so its “true anomaly” is greater.  That’s what we want.]

In the picture, the planet is a quarter of the way into the period of its orbit, so its mean anomaly is 90°.  Its true anomaly is well ahead of this.

After (I imagine) scribbling maze after maze of curves and triangles to find trigonometrical relations, Kepler devised a way to get from the mean anomaly to the true anomaly, but it had to get there by way of a third angle: the “eccentric anomaly.”  This is – wait for it –

The angle, at the center of the ellipse (not at the Sun), from the perihelion point to a point where a perpendicular from the major axis through the planet meets the circle enclosing the orbit.

At the moment shown by the picture, an imaginary plant having the eccentric anomaly has gone around the enclosing circle as far as the imaginary planet with the mean anomaly, and farther, to where it is as far, parallel to the major axis, as the true planet is.

Kepler’s Equation finds this steppingstone from the circle to the ellipse:

E = M + e sin E

– the eccentric anomaly (E) is the mean anomaly (M) plus the eccentricity  (e) multiplied by the sine of the eccentric anomaly  (E).  Abracadabra!

Hold on!  “Eccentric anomaly” is on both sides of the equation.  How can you find the sine of the quantity you haven’t yet found?

You can; it’s called a transcendental equation, which has to be solved by iteration.  You start with a number that’s about as near as you can guess – the mean anomaly will do – and keep running the same equation over and over until the difference between each answer and the last becomes very-very small.

I learned all this from Jean Meeus’s Astronomical Algorithms, chapter on “Equation of Kepler.”  I’ll show you my Fortran subroutine that performs it.  (Until at least a year after I started the Astronomical Calendar I didn’t know a sine, a radian, an anomaly, or a subroutine from a pipsissewa.)  Ee, aom, aoe are my arbitrary names for the eccentricity, mean anomaly, and eccentric anomaly.  I’ll explain some of my Fortran habits another time.

subroutine kepler (ee,aom, aoe)
implicit double precision (a-i,k-p,r-z)
aoe=aom
eedg=ee*57.29577951d0
8001 corr=(aom+eedg*sind(aoe)-aoe) / (1-ee*cosd(aoe))
aoe=aoe+corr
if (dabs(corr).gt..0000001) go to 8001
end

It says: keep adding a correction, using a  slight rearrangement of Kepler’s equation, until the correction becomes less than 0.0000001 of a degree.

Kepler’s equation assumes that the angles are expressed not in degrees but in radians (a radian is the angle made by a wall that is the same length as its distance away, and there are 57.2957795… degrees in a radian.  Computer languages, too, tend to use radians rather than degrees.  I can think better in degrees, so I multiply the eccentricity by that number and, instead of the sin and cos functions provided, use my own sind and cosd.

Back when I was working this out, I must have started by making it run through the loop a few times, then more, to see how small the correction was getting.  Like Pi and the degree-to=radian number and some others, it can’t reach absolute zero.  It starts oscillating either side of zero, and there is a danger that the loop can become infinite.

Procedures such as finding a sine are elaborate, and they happen over and over within one performance of Kepler’s equation.  The equation is just a part – a deep-down part – in the calculation of the position of a moving body.  My program does all of this for every day around a 5-year orbit, and then many times more for drawing other curves and points in the diagram.  For a chart showing, say, the movements of a dozen comets, it gets done thousands of times more.  On earlier computers that I once used, it all might have taken hours.  Now it is done quicker than I can lift my finger off the keyboard.  Kepler had to use pen and paper.  I’m not even sure whether he could save labor by using logarithms; they were invented about the same time.