def METHOD_NAME(self, x, y, theta_E, gamma, center_x=0.0, center_y=0.0):
    gamma = self._gamma_limit(gamma)
    xt1 = x - center_x
    xt2 = y - center_y
    r2 = xt1 * xt1 + xt2 * xt2
    a = np.maximum(r2, 0.000001)
    r = np.sqrt(a)
    alpha = theta_E * (r2 / theta_E**2) ** (1 - gamma / 2.0)
    fac = alpha / r
    f_x = fac * xt1
    f_y = fac * xt2
    return f_x, f_y