code: def fit_harmonic_oscillator(orbit, omega0=[1., 1, 1], minimize_kwargs=None):
      r
      omega0 = np.atleast_1d(omega0)
  
      pot = orbit.hamiltonian.potential
      if pot is None:
          raise ValueError("The orbit object must have an associated potential")
  
      w = np.squeeze(orbit.w(pot.units))
      if w.ndim > 2:
          raise ValueError("Input orbit object must be a single orbit.")
  
      def f(omega, w):
          potential = HarmonicOscillatorPotential(omega=omega, units=pot.units)
          H = (potential.value(w[:3]).decompose(pot.units).value +
               0.5*np.sum(w[3:]**2, axis=0))
          return np.sum(np.squeeze(H - np.mean(H))**2)
  
      if minimize_kwargs is None:
          minimize_kwargs = dict()
      minimize_kwargs['x0'] = omega0
      minimize_kwargs['method'] = minimize_kwargs.get('method', 'Nelder-Mead')
      res = minimize(f, args=(w,), **minimize_kwargs)
  
      if not res.success:
          raise ValueError("Failed to fit toy potential to orbit.")
  
      best_omega = np.abs(res.x)
      return HarmonicOscillatorPotential(omega=best_omega, units=pot.units)
docs: r"""
      Fit the toy harmonic oscillator potential to the sum of the energy
      residuals relative to the mean energy by minimizing the function
  
      .. math::
  
          f(\boldsymbol{\omega}) = \sum_i (\frac{1}{2}v_i^2 + \Phi_{\rm sho}(x_i\,|\,\boldsymbol{\omega}) - <E>)^2
  
      TODO: This should fail if the Hamiltonian associated with the orbit has
            a frame other than StaticFrame
  
      Parameters
      ----------
      orbit : `~gala.dynamics.Orbit`
      omega0 : array_like (optional)
          Initial frequency guess.
      minimize_kwargs : dict (optional)
          Keyword arguments to pass through to `scipy.optimize.minimize`.
  
      Returns
      -------
      omegas : float
          Best-fit harmonic oscillator frequencies.