bar

def delta_t_from_nu(nu, ecc, k=1.0, q=1.0, delta=1e-2):
    """Time elapsed since periapsis for given true anomaly.
    Parameters
    ----------
    nu : float
        True anomaly.
    ecc : float
        Eccentricity.
    k : float
        Gravitational parameter.
    q : float
        Periapsis distance.
    delta : float
        Parameter that controls the size of the near parabolic region.
    Returns
    -------
    delta_t : float
        Time elapsed since periapsis.
    """
    assert -np.pi <= nu < np.pi
    if ecc < 1 - delta:
        # Strong elliptic
        E = nu_to_E(nu, ecc)  # (-pi, pi]
        M = E_to_M(E, ecc)  # (-pi, pi]
        n = np.sqrt(k * (1 - ecc) ** 3 / q**3)
    elif 1 - delta <= ecc < 1:
        E = nu_to_E(nu, ecc)  # (-pi, pi]
        if delta <= 1 - ecc * np.cos(E):
            # Strong elliptic
            M = E_to_M(E, ecc)  # (-pi, pi]
            n = np.sqrt(k * (1 - ecc) ** 3 / q**3)
        else:
            # Near parabolic
            D = nu_to_D(nu)  # (-∞, ∞)
            # If |nu| is far from pi this result is bounded
            # because the near parabolic region shrinks in its vicinity,
            # otherwise the eccentricity is very close to 1
            # and we are really far away
            M = D_to_M_near_parabolic(D, ecc)
            n = np.sqrt(k / (2 * q**3))
    elif ecc == 1:
        # Parabolic
        D = nu_to_D(nu)  # (-∞, ∞)
        M = D_to_M(D)  # (-∞, ∞)
        n = np.sqrt(k / (2 * q**3))
    elif 1 + ecc * np.cos(nu) < 0:
        # Unfeasible region
        return np.nan
    elif 1 < ecc <= 1 + delta:
        # NOTE: Do we need to wrap nu here?
        # For hyperbolic orbits, it should anyway be in
        # (-arccos(-1 / ecc), +arccos(-1 / ecc))
        F = nu_to_F(nu, ecc)  # (-∞, ∞)
        if delta <= ecc * np.cosh(F) - 1:
            # Strong hyperbolic
            M = F_to_M(F, ecc)  # (-∞, ∞)
            n = np.sqrt(k * (ecc - 1) ** 3 / q**3)
        else:
            # Near parabolic
            D = nu_to_D(nu)  # (-∞, ∞)
            M = D_to_M_near_parabolic(D, ecc)  # (-∞, ∞)
            n = np.sqrt(k / (2 * q**3))
    elif 1 + delta < ecc:
        # Strong hyperbolic
        F = nu_to_F(nu, ecc)  # (-∞, ∞)
        M = F_to_M(F, ecc)  # (-∞, ∞)
        n = np.sqrt(k * (ecc - 1) ** 3 / q**3)
    else:
        raise RuntimeError

    return M / n