docs: def rebin_spec(spec, wavnew, oversamp=100, plot=False):
      """
      Rebin a spectrum to a new wavelength array while preserving
      the total flux
  
      Parameters
      ----------
      spec: array-like
          The wavelength and flux to be binned
      wavenew: array-like
          The new wavelength array
  
      Returns
      -------
      np.ndarray
          The rebinned flux
  
      """
      wave, flux = spec
      nlam = len(wave)
      x0 = np.arange(nlam, dtype=float)
      x0int = np.arange((nlam-1.) * oversamp + 1., dtype=float)/oversamp
      w0int = np.interp(x0int, x0, wave)
      spec0int = np.interp(w0int, wave, flux)/oversamp
  
      # Set up the bin edges for down-binning
      maxdiffw1 = np.diff(wavnew).max()
      w1bins = np.concatenate(([wavnew[0]-maxdiffw1],
                               .5*(wavnew[1::]+wavnew[0: -1]),
                               [wavnew[-1]+maxdiffw1]))
  
      # Bin down the interpolated spectrum:
      w1bins = np.sort(w1bins)
      nbins = len(w1bins)-1
      specnew = np.zeros(nbins)
      inds2 = [[w0int.searchsorted(w1bins[ii], side='left'),
                w0int.searchsorted(w1bins[ii+1], side='left')]
               for ii in range(nbins)]
  
      for ii in range(nbins):
          specnew[ii] = np.sum(spec0int[inds2[ii][0]: inds2[ii][1]])
  
      return specnew