def solve(f, M0, e, h):
    """ Newton-Raphson solver for eccentricity anomaly
    from https://stackoverflow.com/questions/20659456/python-implementing-a-numerical-equation-solver-newton-raphson
    Inputs: 
        f (function): function to solve (transcendental ecc. anomaly function)
        M0 (float): mean anomaly
        e (float): eccentricity
        h (float): termination criteria for solver
    Returns: nextE (float): converged solution for eccentric anomaly
        Written by Logan Pearce, 2019
    """
    import numpy as np
    if M0 / (1.0 - e) - np.sqrt(6.0 * (1 - e) / e) <= 0:
        E0 = M0 / (1.0 - e)
    else:
        E0 = (6.0 * M0 / e) ** (1.0 / 3.0)
    lastE = E0
    nextE = lastE + 10 * h
    number = 0
    while abs(lastE - nextE) > h and number < 1001:
        new = f(nextE, e, M0)
        lastE = nextE
        nextE = lastE - new / (1.0 - e * np.cos(lastE))
        number = number + 1
        if number >= 1000:
            nextE = float('NaN')
    return nextE
