def METHOD_NAME(theta_a, theta_b, r_ab, alpha):
    theta_ab = theta_a - theta_b
    return alpha * theta_a + np.sum(
        np.where(
            np.abs(theta_ab) < EPSILON,
            theta_ab / 3,
            1 / np.tanh(theta_ab) - 1 / theta_ab,
        )
        + r_ab
    )