Newton-Raphson iterations

The Newton-Raphson method is employed to solve the following Kepler equations for eccentric anomaly E for a given M:

 M = E - e*sin(E) if e < 1

and

M = e*sinh(E) - E if e > 1

Corrections dE are summed to E(N) each iteration to yield the E(N+1) approximation:

 E(N) = E(N) + dE,

until one of the following conditions is reached:

   abs(dE) < liborb_zero (1.0E-08 by default)

or

   N = liborb\_nrmax (100 by default)

The user can change these values in order to change accuracy or the performance. Once one change such parameters, the accuracy should be checked. The accuracy can be checked via the function M2Err(e,M). This function solves Kepler equation for E, and then uses the result back into Kepler equation to evaluate the Delta M, i.e.:

• M0 is the input mean anomaly
• E0 is the eccentric anomaly gotten from New-Raphson iteration
• deltaM = M0 - ( E0 - e*sin(E0) )