def hessian(self, x, y, Rs, theta_Rs, r_core, center_x=0, center_y=0):

        """
        :param x: x coordinate
        :param y: y coordinate
        :param Rs: scale radius
        :param rho0: central core density
        :param r_core: core radius
        :param center_x:
        :param center_y:
        :return:
        """

        if Rs < 0.0001:
            Rs = 0.0001
        x_ = x - center_x
        y_ = y - center_y
        R = np.sqrt(x_ ** 2 + y_ ** 2)

        rho0 = self._alpha2rho0(theta_Rs=theta_Rs, Rs=Rs, r_core=r_core)

        kappa = self.density_2d(x_, y_, Rs, rho0, r_core)

        gamma1, gamma2 = self.cBurkGamma(R, Rs, rho0, r_core, x_, y_)
        f_xx = kappa + gamma1
        f_yy = kappa - gamma1
        f_xy = gamma2
        return f_xx, f_yy, f_xy