Add eccentric anomaly determination by binary search.

[ktyldev]

Background

The eccentric anomaly E is normally calculated using Newton-Rhapson (NR). This method is straightforward to understand and implement, and usually converges on a good result in a few iterations. The accuracy of the result can be increased by running the method for an arbitrarily large number of iterations.

Problem

Real-time contexts, such as interactive simulations or games the worst-case performance of an algorithm becomes more important than the best or average case, as there may be severe constraints on per-frame processing budgets. An iterative method such as NR must therefore be capped in its maximum number of iterations according to the maximum allowed budget. We must allow for the worst case high-iteration outcome if using NR, as it uses a variable number of iterations depending on the input values of the orbit's eccentricity and the mean anomaly. This makes for inefficient processing, as most of the time NR resolves quickly, but the rest of the application must be restricted so as to allow NR to converge with reasonable accuracy in the worst case.

Solution

In Astronomical Algorithms a method is described to determine E using a binary search (BS). The BS resolves in a constant number of iterations. As a consequence, the best- and average-cases take longer to compute than NR. However, in the worst case, BS far outperforms NR in terms of accuracy per computational cost. Since BS resolves in constant time, it is also much easier to budget for, as it is highly predictable.

[ktyldev]

There is one remaining TODO regarding a sign() calculation. I would usually have extracted this to a function, but that would appear out of place given the convention of the rest of the file and library at large. It is therefore written out twice in the newly implemented function.

[kartikkumar]

@ktyldev thanks for this! I'm a bit swamped at the moment, but I'll try to schedule some time to go through this PR.