docs: def resample(old_wavelengths, new_wavelengths):
      """
      Resample a spectrum to a new wavelengths map while conserving total flux.
  
      :param old_wavelengths:
          The original wavelengths array.
  
      :type old_wavelengths:
          :class:`numpy.array`
  
      :param new_wavelengths:
          The new wavelengths array to resample onto.
  
      :type new_wavelengths:
          :class:`numpy.array`
      """
  
      data = []
      old_px_indices = []
      new_px_indices = []
      for i, new_wl_i in enumerate(new_wavelengths):
  
          # These indices should span just over the new wavelength pixel.
          indices = np.unique(np.clip(
              old_wavelengths.searchsorted(new_wavelengths[i:i + 2], side="left")\
                  + [-1, +1], 0, old_wavelengths.size - 1))
          N = np.ptp(indices)
  
          if N == 0:
              # 'Fake' pixel.
              data.append(np.nan)
              new_px_indices.append(i)
              old_px_indices.extend(indices)
              continue
  
          # Sanity checks.
          assert (old_wavelengths[indices[0]] <= new_wl_i \
              or indices[0] == 0)
          assert (new_wl_i <= old_wavelengths[indices[1]] \
              or indices[1] == old_wavelengths.size - 1)
  
          fractions = np.ones(N)
  
          # Edges are handled as fractions between rebinned pixels.
          _ = np.clip(i + 1, 0, new_wavelengths.size - 1)
          lhs = old_wavelengths[indices[0]:indices[0] + 2]
          rhs = old_wavelengths[indices[-1] - 1:indices[-1] + 1]
          fractions[0]  = (lhs[1] - new_wavelengths[i])/np.ptp(lhs)
          fractions[-1] = (new_wavelengths[_] - rhs[0])/np.ptp(rhs)
  
          # Being binned to a single pixel. Prevent overflow from fringe cases.
          fractions = np.clip(fractions, 0, 1)
          fractions /= fractions.sum()
  
          data.extend(fractions) 
          new_px_indices.extend([i] * N) # Mark the new pixel indices affected.
          old_px_indices.extend(np.arange(*indices)) # And the old pixel indices.
  
      return sparse.csc_matrix((data, (old_px_indices, new_px_indices)),
          shape=(old_wavelengths.size, new_wavelengths.size))