def hessian(self, x, y, center_x = 0, center_y = 0, **kwargs):
        """
        returns Hessian matrix of function d^2f/dx^2, d^f/dy^2, d^2/dxdy
        (un-normalized!!!) interpolated from the numerical deflection table
        """

        diff = 1e-6
        alpha_ra, alpha_dec = self.derivatives(x, y, center_x = center_x, center_y = center_y,
                                               **kwargs)

        alpha_ra_dx, alpha_dec_dx = self.derivatives(x + diff, y, center_x = center_x, center_y = center_y,
                                               **kwargs)
        alpha_ra_dy, alpha_dec_dy = self.derivatives(x, y + diff, center_x = center_x, center_y = center_y,
                                               **kwargs)

        dalpha_rara = (alpha_ra_dx - alpha_ra) / diff
        dalpha_radec = (alpha_ra_dy - alpha_ra) / diff
        dalpha_decdec = (alpha_dec_dy - alpha_dec) / diff

        f_xx = dalpha_rara
        f_yy = dalpha_decdec
        f_xy = dalpha_radec

        return f_xx, f_yy, f_xy