def METHOD_NAME(self, x, y, theta_E, r_trunc, center_x=0, center_y=0):
    x_shift = x - center_x
    y_shift = y - center_y
    r = np.sqrt(x_shift * x_shift + y_shift * y_shift)
    if isinstance(r, int) or isinstance(r, float):
        if r < r_trunc:
            f_ = theta_E * r
        elif r < 2 * r_trunc:
            f_ = theta_E * r_trunc + 1.0 / 2 * theta_E * (3 - r / r_trunc) * (
                r - r_trunc
            )
        else:
            f_ = 3.0 / 2 * theta_E * r_trunc
    else:
        f_ = np.zeros_like(r)
        f_[r < r_trunc] = theta_E * r[r < r_trunc]
        r_ = r[(r < 2 * r_trunc) & (r > r_trunc)]
        f_[
            (r < 2 * r_trunc) & (r > r_trunc)
        ] = theta_E * r_trunc + 1.0 / 2 * theta_E * (3 - r_ / r_trunc) * (
            r_ - r_trunc
        )
        f_[r > 2 * r_trunc] = 3.0 / 2 * theta_E * r_trunc
    return f_