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
    E = theta_E / ((3.0 - gamma) / 2.0) ** (1.0 / (1.0 - gamma))
    # E = phi_E_spp
    eta = -gamma + 3.0
    P2 = xt1**2 + xt2**2
    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
    kappa = (
        1.0
        / eta
        * (a / E**2) ** (eta / 2 - 1)
        * ((eta - 2) * (xt1**2 + xt2**2) / a + (1 + 1))
    )
    gamma1 = (
        1.0
        / eta
        * (a / E**2) ** (eta / 2 - 1)
        * ((eta / 2 - 1) * (2 * xt1**2 - 2 * xt2**2) / a)
    )
    gamma2 = (
        4 * xt1 * xt2 * (1.0 / 2 - 1 / eta) * (a / E**2) ** (eta / 2 - 2) / E**2
    )
    f_xx = kappa + gamma1
    f_yy = kappa - gamma1
    f_xy = gamma2
    return f_xx, f_xy, f_xy, f_yy