How to use the orbitize.kepler._calc_ecc_anom function in orbitize

def test_analytical_ecc_anom_solver(use_c = False):
    """
    Test orbitize.kepler._calc_ecc_anom() in the analytical solver regime (e > 0.95) by comparing the mean anomaly computed from
    _calc_ecc_anom() output vs the input mean anomaly
    """
    mean_anoms = np.linspace(0,2.0*np.pi,1000)
    eccs = np.linspace(0.95,0.999999,100)
    for ee in eccs:
        ecc_anoms = kepler._calc_ecc_anom(mean_anoms, ee, tolerance=1e-9, use_c=use_c)
        calc_mm = (ecc_anoms - ee*np.sin(ecc_anoms)) % (2*np.pi) # plug solutions into Kepler's equation
        for meas, truth in zip(calc_mm, mean_anoms):
            assert angle_diff(meas, truth) == pytest.approx(0.0, abs=threshold)

def profile_iterative_ecc_anom_solver(n_orbits = 1000, use_c = True):
    """
    Test orbitize.kepler._calc_ecc_anom() in the iterative solver regime (e < 0.95) by comparing the mean anomaly computed from
    _calc_ecc_anom() output vs the input mean anomaly
    """

    mean_anoms=np.linspace(0, 2.0*np.pi,n_orbits)
    eccs=np.linspace(0,0.9499999, n_orbits)
    for ee in eccs:
        ecc_anoms = kepler._calc_ecc_anom(mean_anoms, ee, tolerance=1e-9, use_c = use_c)

def test_iterative_ecc_anom_solver(use_c = False):
    """
    Test orbitize.kepler._calc_ecc_anom() in the iterative solver regime (e < 0.95) by comparing the mean anomaly computed from
    _calc_ecc_anom() output vs the input mean anomaly
    """
    mean_anoms = np.linspace(0,2.0*np.pi,100)
    eccs = np.linspace(0,0.9499999,100)
    for ee in eccs:
        ecc_anoms = kepler._calc_ecc_anom(mean_anoms, ee, tolerance=1e-9, use_c=use_c)
        calc_ma = (ecc_anoms - ee*np.sin(ecc_anoms)) % (2*np.pi) # plug solutions into Kepler's equation
        for meas, truth in zip(calc_ma, mean_anoms):
            assert angle_diff(meas, truth) == pytest.approx(0.0, abs=threshold)

def profile_mikkola_ecc_anom_solver(n_orbits = 1000, use_c = True):
    """
    Test orbitize.kepler._calc_ecc_anom() in the iterative solver regime (e < 0.95) by comparing the mean anomaly computed from
    _calc_ecc_anom() output vs the input mean anomaly
    """
    mean_anoms=np.linspace(0, 2.0*np.pi,n_orbits)
    eccs=np.linspace(.95,0.999999, n_orbits)
    for ee in eccs:
        ecc_anoms = kepler._calc_ecc_anom(mean_anoms, ee, use_c = use_c)

#first thing we have to do is figure out the average angular velocity in rad/day:
    P = (np.sqrt(sma**3/(m0 + m1)))*365.25 #returning period in days
    n = 2*np.pi/P

    #time of periastron in days:
    tperi = P*tau

    #Mean anomaly as a function of time in radians for both times:
    M = n*(time - tperi)
    M = np.mod(M,2*np.pi)

    ecc_arr = np.ones(len(M))*ecc

    #eccentric anomalies:
    E = orbitize.kepler._calc_ecc_anom(M,ecc_arr)

    nu = 2*np.arctan(np.sqrt((1 + ecc)/(1 - ecc))*np.tan(E/2))

    u = argp + nu

    #getting semi-major axis in AU from the period:

    #dx/dt,dy/dt are arrays of length time
    r = 2*np.pi*sma/(P*np.sqrt(1 - ecc**2))

    dxdt = r*(np.cos(inc)*np.cos(lan)*(ecc*np.cos(argp) + np.cos(u)) - \
              np.sin(lan)*(ecc*np.sin(argp) + np.sin(u)))
    dydt = -r*(np.cos(inc)*np.sin(lan)*(ecc*np.cos(argp) + np.cos(u)) + \
              np.cos(lan)*(ecc*np.sin(argp) + np.sin(u)))

    #getting velocity of star: