docs: def cross_correlate(observed, template_dispersion, template_fluxes,
      rebin="template", wavelength_range=None, continuum_degree=-1,
      apodize=0.10, rescale=False, z_limits=None, full_output=False):
      """
      Cross-correlate the observed spectrum against template fluxes.
      """
  
      template_fluxes = np.atleast_2d(template_fluxes)
      template_dispersion = np.array(template_dispersion)
      if template_dispersion.shape[0] != template_fluxes.shape[1]:
          raise ValueError("template dispersion must have size (N_pixels,) and "\
              "template fluxes must have size (N_models, N_pixels)")
  
      try:
          continuum_degree = int(continuum_degree)
      except (TypeError, ValueError):
          raise TypeError("continuum order must be an integer-like object")
  
      if not (1 > apodize >= 0):
          raise ValueError("apodize fraction must be between 0 and 1")
  
      if rebin.lower() == "template":
          # Put the template fluxes onto the observed dispersion map.
          dispersion = observed.disp
          template_flux = template_fluxes \
              * resample(template_dispersion, observed.disp)
          observed_flux = observed.flux.copy()
          observed_ivar = observed.ivariance.copy()
  
      elif rebin.lower() == "observed":
          # Put the observed fluxes onto the template dispersion map.
          dispersion = template_dispersion
          template_flux = template_fluxes
          mat = resample(observed.disp, template_dispersion)
          observed_flux = observed.flux.copy() * mat
          observed_ivar = observed.ivariance.copy() * mat
          
      else:
          raise ValueError("rebin must be either `template` or `observed`")
  
      if wavelength_range is not None:
          if not isinstance(wavelength_range, (tuple, list, np.ndarray)) \
          or len(wavelength_range) != 2:
              raise TypeError("wavelength range must be a two length tuple")
  
          indices = np.clip(
              template_dispersion.searchsorted(wavelength_range) + [0, 1],
              0, template_dispersion.size)
          template_dispersion = template_dispersion[indices[0]:indices[1]]
          template_flux = template_flux[indices[0]:indices[1]]
          observed_flux = observed_flux[indices[0]:indices[1]]
  
      N_templates = template_fluxes.shape[0]
  
      # Clip out the non-finite edges, if they exist.
      _ = np.where(np.isfinite(template_flux[0] * observed_flux))[0]
      l_idx, u_idx = _[0], _[-1]
  
      dispersion = dispersion[l_idx:u_idx]
      observed_flux = observed_flux[l_idx:u_idx]
      observed_ivar = observed_ivar[l_idx:u_idx]
      template_flux = template_flux[:, l_idx:u_idx]
  
      # Ensure an even number of points.
      N = u_idx - l_idx
      N = N - 1 if N % 2 > 0 else N
      dispersion = dispersion[:N]
      observed_flux = observed_flux[:N]
      observed_ivar = observed_ivar[:N]
      template_flux = template_flux[:, :N]
  
      finite = np.isfinite(observed_flux * observed_ivar) * (observed_flux > 1e-3)
      
      # Continuum.
      if continuum_degree >= 0:
          coefficients = np.polyfit(dispersion[finite], observed_flux[finite],
              continuum_degree)
          observed_flux /= np.polyval(coefficients, dispersion)
  
      # Interpolate over non-finite pixels.
      observed_flux[~finite] = np.interp(dispersion[~finite], dispersion[finite],
          observed_flux[finite], left=1, right=1)
  
      # Scale the flux level to that the template intensities
      if rescale:
          observed_flux = (observed_flux * template_flux.ptp()) \
              + template_flux.min()
  
      # Apodize edges.
      edge_buffer = apodize * (dispersion[-1] - dispersion[0])
      low_w_indices = np.nonzero(dispersion < dispersion[0] + edge_buffer)[0]
      high_w_indices = np.nonzero(dispersion > dispersion[-1] - edge_buffer)[0]
  
      apod_curve = np.ones(N, dtype='d')
      apod_curve[low_w_indices] = (1.0 + np.cos(np.pi*(
          1.0 - (dispersion[low_w_indices] - dispersion[0])/edge_buffer)))/2.
      apod_curve[high_w_indices] = (1.0 + np.cos(np.pi*(
          1.0 - (dispersion[-1] - dispersion[high_w_indices])/edge_buffer)))/2.
  
      apod_observed_flux = observed_flux * apod_curve
      apod_template_flux = template_flux * apod_curve
  
      fft_observed_flux = np.fft.fft(apod_observed_flux)
      fft_template_flux = np.fft.fft(apod_template_flux)
      template_flux_corr = (fft_observed_flux * fft_template_flux.conjugate())
      template_flux_corr /= np.sqrt(
          np.inner(apod_observed_flux, apod_observed_flux))
  
      z_array = np.array(dispersion.copy())/dispersion[N/2] - 1.0
  
      z = np.ones(N_templates) * np.nan
      z_err = np.ones(N_templates) * np.nan
      R = np.ones(N_templates) * np.nan
      for i in range(N_templates):
  
          denominator = np.sqrt(
              np.inner(apod_template_flux[i, :], apod_template_flux[i, :]))
          flux_correlation = template_flux_corr[i, :] / denominator 
          correlation = np.fft.ifft(flux_correlation).real
  
          # Reflect about zero
          ccf = np.zeros(N)
          ccf[:N/2] = correlation[N/2:]
          ccf[N/2:] = correlation[:N/2]
  
          # Get height and redshift of best peak
          h = ccf.max()
  
          # Scale the CCF
          #ccf -= ccf.min()
          #ccf *= (h/ccf.max())
          
          R[i] = h
          z[i] = z_array[ccf.argmax()]
          try:
              z_err[i] = (np.ptp(z_array[np.where(ccf >= 0.5*h)])/2.35482)**2
          except ValueError:
              continue
  
      # Limits:
      if z_limits is not None:
          R[~((z_limits[1] >= z) * (z >= z_limits[0]))] = np.nan
  
      # Re-measure the velocity at the best peak.
      index = np.nanargmax(R)
  
      denominator = np.sqrt(
          np.inner(apod_template_flux[index, :], apod_template_flux[index, :]))
  
      flux_correlation = template_flux_corr[index, :] / denominator
      correlation = np.fft.ifft(flux_correlation).real
  
      # FReflect about zero.
      ccf = np.zeros(N)
      ccf[:N/2] = correlation[N/2:]
      ccf[N/2:] = correlation[:N/2]
  
      h = ccf.max()
      ccf -= ccf.min()
      ccf *= h/ccf.max()
  
      # Fit +/- 5 pixels
      idx = np.argmax(ccf) - 3, np.argmax(ccf) + 3
  
      coeffs = np.polyfit(z_array[idx[0]:idx[1]], ccf[idx[0]:idx[1]], 2)
      x_i = np.linspace(z_array[idx[0]], z_array[idx[1]], 1000)
      y_i = np.polyval(coeffs, x_i)
  
      # Update the value
      z[index] = x_i[y_i.argmax()]
  
      """
      # Fit the profile peak.
      fwhm = np.ptp(z_array[np.where(ccf >= 0.5)])
      p0 = np.array([z[index], fwhm/2.355, 1.0]) # mu, sigma, peak.
      use = (z_array > (p0[0] - 3 * p0[1])) * ((p0[0] + 3*p0[1]) > z_array)
      x, y = z_array[use], ccf[use]
  
      # Fit profile:
      f = lambda x, mu, sigma, peak: peak * np.exp(-(mu - x)**2/(2*sigma**2)) 
      
      import matplotlib.pyplot as plt
      plt.plot(z_array, ccf, c='k')
      plt.plot(x, y, c='b')
      plt.plot(x, f(x, *p0), c='r')
  
      p1 = curve_fit(f, x, y, p0=p0)
      plt.plot(x, f(x, *p1[0]), c='m')
  
      plt.gca().set_ylim(0,1)
      raise a
  
      p1 = curve_fit(f, x, y, p0=p0)
      """
  
      return (z * c, z_err * c, R)