markley.py

import numpy as np

def esolve_markley(M, ecc):
    M = M % (2*np.pi)
    flip = M > np.pi
    M[flip] = 2*np.pi - M[flip]
    alpha = (3*np.pi + 1.6*(np.pi-np.abs(M))/(1+ecc) )/(np.pi - 6/np.pi)
    d = 3*(1 - ecc) + alpha*ecc
    r = 3*alpha*d * (d-1+ecc)*M + M**3
    q = 2*alpha*d*(1-ecc) - M**2
    w = (np.abs(r) + np.sqrt(q**3 + r**2))**(2/3)
    E = (2*r*w/(w**2 + w*q + q**2) + M) / d
    f_0 = E - ecc*np.sin(E) - M
    f_1 = 1 - ecc*np.cos(E)
    f_2 = ecc*np.sin(E)
    f_3 = 1-f_1
    d_3 = -f_0/(f_1 - 0.5*f_0*f_2/f_1)
    d_4 = -f_0/(f_1 + 0.5*d_3*f_2 + (d_3**2)*f_3/6)
    E = E -f_0/(f_1 + 0.5*d_4*f_2 + d_4**2*f_3/6 - d_4**3*f_2/24)
    E[flip] = 2*np.pi - E[flip]
    return E