code: def hessian(self, x, y, coeffs, beta, center_x=0, center_y=0):
          
          shapelets = self._createShapelet(coeffs)
          r, phi = param_util.cart2polar(x, y, center=np.array([center_x, center_y]))
          kappa_shapelets=self._kappaShapelets(shapelets, beta)
          gamma1_shapelets, gamma2_shapelets=self._gammaShapelets(shapelets, beta)
          kappa_value=self._shapeletOutput(r, phi, beta, kappa_shapelets)
          gamma1_value=self._shapeletOutput(r, phi, beta, gamma1_shapelets)
          gamma2_value=self._shapeletOutput(r, phi, beta, gamma2_shapelets)
          f_xx = kappa_value + gamma1_value
          f_xy = gamma2_value
          f_yy = kappa_value - gamma1_value
          return f_xx, f_yy, f_xy
docs: returns Hessian matrix of function d^2f/dx^2, d^f/dy^2, d^2/dxdy