Kepler's equation without misfires
A planetarium program shows comets in the sky. To learn where a comet is on its orbit at a given moment, it has to solve Kepler's equation
$$E - e \sin E = M$$for the eccentric anomaly E. Here e is the eccentricity of the orbit (0 ≤ e < 1), and M is the mean anomaly, the time elapsed since perihelion as a fraction of the period, times 2π. The solution in [0, 2π) is unique. The planetarium solves the equation for every comet many times a second and must never get it wrong: a comet that jumps to the wrong place is noticed at once.
Input
The first line holds the number of queries q (1 ≤ q ≤ 40,000). Each of the next q lines holds two real numbers e and M (0 ≤ e ≤ 0.9999, 0 ≤ M < 2π).
Output
For each query, E in [0, 2π). An answer is accepted with an absolute or relative error of 10⁻⁹.
Examples
Example 1
4 0.0000000000 1.0000000000 0.0167000000 2.0000000000 0.9679000000 3.3602000000 0.9920000000 0.4084070450
1.000000000000 2.015078752618 3.252791912317 1.382957944904
Limits
- Time · Python
- 5000 ms
- Time · Go
- 1000 ms
- Memory
- 256 MB
- Judging
- 6 tests
Walkthrough
Open it step by step: each step gives away a little more. You can jump straight to the solution, but a problem solved with one hint stays with you longer.
1Hint 1
The standard recipe is Newton's method from E₀ = M. Try it on the fourth line of the sample, e = 0.992, M = 0.13π, and print the first ten iterates.
2Hint 2
The function f(E) = E − e sin E − M is increasing (its derivative 1 − e cos E ≥ 1 − e > 0), f(0) ≤ 0 and f(2π) > 0. Which method cannot miss under such conditions?
3The idea
The root always lies in [0, 2π], with different signs at the ends, so bisection is guaranteed to find it, in 43 steps to 10⁻¹². Newton is faster, but with e around 0.99 and small M the derivative at the start is almost zero, and the tangent throws the iterate hundreds of radians away. A reliable solution combines the two: keep an interval [lo, hi] with different signs at its ends, take a Newton step, and if it lands outside the interval, take a bisection step instead.
4Walkthrough
At each step compute f(E) and move one of the ends by its sign: negative means the root is to the right, lo = E; positive means hi = E. The interval only shrinks and always contains the root. The next iterate is the Newton step if it lands strictly inside (lo, hi); otherwise the midpoint of the interval. Stop when the step is below 10⁻¹⁵ or the interval has become shorter than that.
Start from π when e > 0.8, as celestial mechanics references advise: in our checks Newton from π never failed. But there is no need to rely on it: the bisection backup makes the solution reliable from any start.
The third test is real comets from the JPL database on January 1 of each year from 2000 to 2060. Plain Newton from M goes wrong on it eleven times, for comets Hale–Bopp and NEOWISE, whose eccentricities are 0.995 and 0.9992.
4+Solution
5Complexity
The safeguarded Newton almost always takes 4–12 steps, and where the tangent throws it out of the interval, bisection steps come to the rescue, each halving the interval. In total O(q) with a small constant.
6Common mistakes
- Newton from E₀ = M without a backup. It converges almost always, but with e ≥ 0.99 and M up to 0.5 it sometimes does not converge at all. The fourth and fifth tests are made of such cases.
- A fixed five Newton steps. Enough for circular orbits, not for almost parabolic ones.
- Kepler's iteration E ← M + e sin E. It always converges, but with e = 0.9999 the error shrinks only 1.0001 times per step: reaching 10⁻⁹ takes hundreds of thousands of steps per query.
- An answer outside [0, 2π). E + 2π is not a solution: the left-hand side of the equation grows by 2π too.