code: def hessian(self, x, y, theta_E, r_trunc, center_x=0, center_y=0):
          
          x_shift = x - center_x
          y_shift = y - center_y
          dphi_dr = self._dphi_dr(x_shift, y_shift, theta_E, r_trunc)
          d2phi_dr2 = self._d2phi_dr2(x_shift, y_shift, theta_E, r_trunc)
          dr_dx, dr_dy = self._dr_dx(x, y)
          d2r_dx2, d2r_dy2, d2r_dxy = self._d2r_dx2(x_shift, y_shift)
          f_xx = d2r_dx2*dphi_dr + dr_dx**2*d2phi_dr2
          f_yy = d2r_dy2*dphi_dr + dr_dy**2*d2phi_dr2
          f_xy = d2r_dxy*dphi_dr + dr_dx*dr_dy*d2phi_dr2
          return f_xx, f_yy, f_xy
docs: returns Hessian matrix of function d^2f/dx^2, d^f/dy^2, d^2/dxdy