code: def hessian(self, q, t=0.):
          
          if (self.R is not None and
                  not np.allclose(np.diag(self.R), 1., atol=1e-15, rtol=0)):
              raise NotImplementedError("Computing Hessian matrices for rotated "
                                        "potentials is currently not supported.")
          q = self._remove_units_prepare_shape(q)
          orig_shape,q = self._get_c_valid_arr(q)
          t = self._validate_prepare_time(t, q)
          ret_unit = 1 / self.units['time']**2
          hess = np.moveaxis(self._hessian(q, t=t), 0, -1)
          return hess.reshape((orig_shape[0], orig_shape[0]) + orig_shape[1:]) * ret_unit
docs: Compute the Hessian of the potential at the given position(s).
  
          Parameters
          ----------
          q : `~gala.dynamics.PhaseSpacePosition`, `~astropy.units.Quantity`, array_like
              The position to compute the value of the potential. If the
              input position object has no units (i.e. is an `~numpy.ndarray`),
              it is assumed to be in the same unit system as the potential.
  
          Returns
          -------
          hess : `~astropy.units.Quantity`
              The Hessian matrix of second derivatives of the potential. If the input
              position has shape ``q.shape``, the output energy will have shape
              ``(q.shape[0],q.shape[0]) + q.shape[1:]``. That is, an ``n_dim`` by
              ``n_dim`` array (matrix) for each position.