Solving Kepler's equation by fixed point iteration

The position of a planet in its elliptical orbit as a function of time can be found through solving Kepler's equation. Let t be time elapsed since perihelion (closest approach to the Sun), and T the period of the orbit - the amount of time it takes the planet to complete one revolution. Let

M = 2 pi (t/T)

so that M varies from 0 to 2 pi as an orbit is traversed. Then

M = E - e sin(E)

where E is the angle of a point in theKepler circle constructed on the orbit. This is called Kepler's equation, and is derived from Kepler's Second Law, which states that the radial area swept out by a planet is proportional to time. The quantity M is called the mean anomaly in the literature since it measures the position of a fictitious planet moving uniformly with respect to time, and E is called the eccentric anomaly. They agree if the orbit is a circle. The position is given in terms of E as (a cos(E), b sin(E)) where a and b are the semi-major and semi-minor axes of the orbit. The quantity b is related to a and e by the formula

b = a sqrt(1 - e^2) .

Of course, Kepler's equation will tell us easily what t is if we are given E , but going in the opposite direction involves solving a transcendental equation for E. We can rewrite it as

E = M + e sin(E)

which means that E is a fixed point of the function E -> M + e sin(E). To find the fixed point of a function f(x) is simple if the root X we are looking for is stable - that is to say if |f'(X)| is less than 1. This is always true for Kepler's equation if the condition e < 1 is valid, which in fact always holds for elliptical orbit. The convergence rate will decresae, however, as e approaches 1. The following very simple applet solves Kepler's equation for ellipses by fixed point iteration. Set e and M (press `carriage return' in a window to enter the data), and E will be set to M. `Iterate' changes E to M + e sin(E). Enough iterations will converge sooner or later.