docs: def calc_ispec(model, ph):
      """Compute isotropic spectrum `phr` of `ph` from 2D spectrum.
  
      Parameters
      ----------
      model : pyqg.Model instance
          The model object from which `ph` originates
      ph : complex array
          The field on which to compute the variance
  
      Returns
      -------
      kr : array
          isotropic wavenumber
      phr : array
          isotropic spectrum
      """
  
      if model.kk.max()>model.ll.max():
          kmax = model.ll.max()
      else:
          kmax = model.kk.max()
  
      # create radial wavenumber
      dkr = np.sqrt(model.dk**2 + model.dl**2)
      kr =  np.arange(dkr/2.,kmax+dkr,dkr)
      phr = np.zeros(kr.size)
  
      for i in range(kr.size):
          fkr =  (model.wv>=kr[i]-dkr/2) & (model.wv<=kr[i]+dkr/2)
          dth = pi / (fkr.sum()-1)
          phr[i] = ph[fkr].sum() * kr[i] * dth
  
      return kr, phr