docs: def old_projection(model, data, theta=None, chain=None, n=100, extents=None, 
      uncertainties=True, title=None, fig=None, figsize=None):
      """
      Project the maximum likelihood values and sampled posterior points as 
      spectra.
  
      :param model:
          The model employed.
  
      :type model:
          :class:`sick.models.Model`
  
      :param data:
          The observed spectra.
  
      :type data:
          iterable of :class:`sick.specutils.Spectrum1D` objects
  
      :param theta: [optional]
          The optimised model parameters given the data. Either theta
          or chain should be given.
  
      :type theta:
          dict
  
      :param chain: [optional]
          The chain of sampled parameters.
  
      :type chain:
          :class:`numpy.ndarray`    
  
      :param extents: [optional]
          The wavelength extents to plot for each channel in the form of 
          [(min_chan_1, max_chan_1), ..., (min_chan_N, max_chan_N)]
      
      :type extents:
          tuple or None
  
      :param uncertainties: [optional]
          Show uncertainty of the data points.
  
      :type uncertainties:
          bool
  
      :param title: [optional]
          Title to set for the top axes.
  
      :type title:
          str
  
      :param fig: [optional]
          Overplot onto the provided figure object.
  
      :type fig:
          :class:`matplotlib.Figure` or None
      
      :param figsize: [optional]
          The figure size (x-dimension, y-dimension) in inches.
  
      :type figsize:
          tuple or None
  
      :raises ValueError:
          If a ``fig`` is provided with the incorrect number of axes.
  
      :raise TypeError:
          If the ``data`` are not provided in the correct type.
  
      :returns:
          The projection figure.
  
      :rtype:
          :class:`maplotlib.Figure`
      """
      if not isinstance(data, (tuple, list)) or \
      any([not isinstance(each, specutils.Spectrum1D) for each in data]):
          raise TypeError("Data must be a list-type of Spectrum1D objects.")
  
      K = len(data)
  
      factor = 2.0
      lbdim = 0.5 * factor
      trdim = 0.2 * factor
      whspace = 0.10
      width = np.max([len(each.disp) for each in data])/500.
      height = factor*K + factor * (K - 1.) * whspace
      dimy = lbdim + height + trdim
      dimx = lbdim + width + trdim
  
      if figsize is None:
          figsize = (dimx, dimy)
      if fig is None:
          fig, axes = plt.subplots(K, 1, figsize=figsize)
  
      else:
          try:
              axes = np.array(fig.axes).reshape((K, 1))
          except:
              raise ValueError("Provided figure has {0} axes, but data has "
                  "parameters K={1}".format(len(fig.axes), K))
  
      if chain is not None:
  
          flat_chain = chain.reshape(-1, len(model.parameters))
          map_theta = np.mean(flat_chain, axis=0)
  
          try:
              map_fluxes = model(data=data, **model._dictify_theta(map_theta))
          except:
              logger.warn("Could not draw MAP fluxes from posterior")
  
          if n > 0:
              # Draw samples from sampler.chain and compute spectra for them
              sampled_fluxes = []
              n_samples = len(flat_chain)
  
              for i in range(n):
                  sampled_theta = dict(zip(model.parameters,
                      flat_chain[np.random.randint(0, n_samples)]))
                  try:
                      sampler_flux = model(data=data, **sampled_theta)
                  except:
                      logger.warn("Could not draw sample flux from posterior")
                      continue
                  else:
                      sampled_fluxes.append(sampler_flux)
          
      elif theta is not None:
  
          sampled_fluxes = []
          map_fluxes = model(data=data, **model._dictify_theta(theta))
          
      else:
          raise ValueError("either theta or chain should be given")
  
  
      if len(data) == 1:
          axes = [axes]
  
      for k, (map_flux, observed_spectrum) in enumerate(zip(map_fluxes, data)):
  
          ax = axes[k]
  
          # Draw the random samples from the chain
          if n > 0:
              for sampled_flux in sampled_fluxes:
                  ax.plot(observed_spectrum.disp, sampled_flux[k], color="r",
                      alpha=0.5, zorder=90)
  
          # Draw the ML spectra
          ax.plot(observed_spectrum.disp, map_flux, color="r", lw=2, zorder=100)
  
          # Plot the data
          if uncertainties:
              ax.fill_between(observed_spectrum.disp,
                  observed_spectrum.flux - observed_spectrum.variance**0.5,
                  observed_spectrum.flux + observed_spectrum.variance**0.5,
                  facecolor="#cccccc", edgecolor="#666666", zorder=-1)
          ax.plot(observed_spectrum.disp, observed_spectrum.flux, color="k",
              zorder=10)
  
          # By default only show common overlap between the model and data
          if extents is None:
              finite_data = np.isfinite(observed_spectrum.flux)
              finite_model = np.isfinite(map_flux)
              finite_points = (finite_model, finite_data)
  
              x_extent = [
                  np.max([observed_spectrum.disp[s][0]  for s in finite_points]),
                  np.min([observed_spectrum.disp[s][-1] for s in finite_points]),
              ]
  
              indices = observed_spectrum.disp.searchsorted(x_extent)
              finite_flux = observed_spectrum.flux[indices[0]:indices[1]]
  
              if len(finite_flux) > 0:
                  #y_extent = [
                  #    0.9 * np.min(finite_flux[np.isfinite(finite_flux)]),
                  #    1.1 * np.max(finite_flux[np.isfinite(finite_flux)])
                  #]
                  ax.set_ylim([0.9, 1.1] * np.percentile(finite_flux[np.isfinite(finite_flux)], [0.5, 99.5]))
  
              ax.set_xlim(x_extent)
  
          else:
              ax.set_xlim(extents[k][0])
              ax.set_ylim(extents[k][1])
  
          # Labels and ticks
          if not (k < K - 1):
              ax.set_xlabel("Wavelength, $\lambda$ ($\AA$)")
  
          ax.set_ylabel("Flux, $F_\lambda$")
          ax.yaxis.set_label_coords(-0.05, 0.5)
  
          ax.xaxis.set_major_locator(MaxNLocator(5))
          ax.yaxis.set_major_locator(MaxNLocator(5))
          [l.set_rotation(45) for l in ax.get_yticklabels()]
  
  
      if title is not None and isinstance(title, (str, unicode)):
          axes[0].set_title(title)
  
      return fig