def METHOD_NAME(self, x, y, theta_E, gamma, e1, e2, center_x=0, center_y=0):
    phi_G, q = param_util.ellipticity2phi_q(e1, e2)
    gamma, q = self._param_bounds(gamma, q)
    phi_E_new = theta_E * q
    x_shift = x - center_x
    y_shift = y - center_y
    E = phi_E_new / (((3 - gamma) / 2.0) ** (1.0 / (1 - gamma)) * np.sqrt(q))
    # E = phi_E
    eta = float(-gamma + 3)
    cos_phi = np.cos(phi_G)
    sin_phi = np.sin(phi_G)
    xt1 = cos_phi * x_shift + sin_phi * y_shift
    xt2 = -sin_phi * x_shift + cos_phi * y_shift
    xt2difq2 = xt2 / (q * q)
    P2 = xt1 * xt1 + xt2 * xt2difq2
    if isinstance(P2, int) or isinstance(P2, float):
        a = max(0.000001, P2)
    else:
        a = np.empty_like(P2)
        p2 = P2[P2 > 0]  # in the SIS regime
        a[P2 == 0] = 0.000001
        a[P2 > 0] = p2
    fac = 1.0 / eta * (a / (E * E)) ** (eta / 2 - 1) * 2
    f_x_prim = fac * xt1
    f_y_prim = fac * xt2difq2
    f_x = cos_phi * f_x_prim - sin_phi * f_y_prim
    f_y = sin_phi * f_x_prim + cos_phi * f_y_prim
    return f_x, f_y