code: def hessian(self, x, y, theta_E, s, q):
          
          alpha_ra, alpha_dec = self.derivatives(x, y, theta_E, s, q)
          diff = self._diff
          alpha_ra_dx, alpha_dec_dx = self.derivatives(x + diff, y, theta_E, s, q)
          alpha_ra_dy, alpha_dec_dy = self.derivatives(x, y + diff, theta_E, s, q)
  
          f_xx = (alpha_ra_dx - alpha_ra) / diff
          f_xy = (alpha_ra_dy - alpha_ra) / diff
          # f_yx = (alpha_dec_dx - alpha_dec)/diff
          f_yy = (alpha_dec_dy - alpha_dec) / diff
          return f_xx, f_yy, f_xy
docs: returns Hessian matrix of function d^2f/dx^2, d^f/dy^2, d^2/dxdy