09b-Specutils/Specutils_analysis_solutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Specutils Analysis\n",
    "\n",
    "![Specutils: An Astropy Package for Spectroscopy](data/specutils_logo.png)\n",
    "\n",
    "\n",
    "This notebook provides an overview of some of the spectral analysis capabilities of the Specutils Astropy coordinated package.  While this notebook is intended as an interactive introduction to specutils at the time of its writing, the canonical source of information for the package is the latest version's documentation: \n",
    "\n",
    "https://specutils.readthedocs.io\n",
    "\n",
    "Note that the below assumes you have knowledge of the material in the [overview notebook](Specutils_overview.ipynb).  If this is not the case you may wish to review that notebook before proceding here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports\n",
    "\n",
    "We start with some fundamental imports for working with specutils and simple visualization of spectra:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import astropy.units as u\n",
    "\n",
    "import specutils\n",
    "from specutils import Spectrum1D, SpectralRegion\n",
    "specutils.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for plotting:\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "# for showing quantity units on axes automatically:\n",
    "from astropy.visualization import quantity_support\n",
    "quantity_support();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sample Spectrum and SNR\n",
    "\n",
    "For use below, we also load the sample SDSS spectrum downloaded in the [overview notebook](Specutils_overview.ipynb).  See that notebook if you have not yet downloaded this spectrum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from urllib.request import urlretrieve\n",
    "\n",
    "sdss_spectrum_path = 'data/sdss_spectrum.fits'\n",
    "\n",
    "url = 'https://data.sdss.org/sas/dr16/sdss/spectro/redux/26/spectra/1323/spec-1323-52797-0012.fits'\n",
    "urlretrieve(url, sdss_spectrum_path)\n",
    "\n",
    "sdss_spec = Spectrum1D.read(sdss_spectrum_path, format='SDSS-III/IV spec')\n",
    "plt.step(sdss_spec.wavelength, sdss_spec.flux);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example file already has uncertainties, but they are initially in inverse variance form.  We convert that to standard deviation form to simplify some of the operations below, and for the same region unmask the whole spectrum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sdss_spec.uncertainty.uncertainty_type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.nddata import StdDevUncertainty\n",
    "\n",
    "sdss_spec.uncertainty = StdDevUncertainty(sdss_spec.uncertainty.quantity**-0.5)\n",
    "sdss_spec.uncertainty.uncertainty_type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sdss_spec.mask[:] = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With these uncertainties, it is straightforward to use one of the fundamental quantifications of a spectrum: the whole-spectrum signal-to-noise ratio:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from specutils import analysis\n",
    "\n",
    "analysis.snr(sdss_spec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spectral Regions\n",
    "\n",
    "Most analysis required on a spectrum requires specification of a part of the spectrum - e.g., a spectral line.  Because such regions may have value independent of a particular spectrum, they are represented as objects distrinct from a given spectrum object.  Below we outline a few ways such regions are specified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ha_region = SpectralRegion((6563-50)*u.AA, (6563+50)*u.AA)\n",
    "ha_region"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regions can also be raw pixel values (although of course this is more applicable to a specific spectrum):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pixel_region = SpectralRegion(2100*u.pixel, 2600*u.pixel)\n",
    "pixel_region"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Additionally, *multiple* regions can be in the same `SpectralRegion` object. This is useful for e.g. measuring multiple spectral features in one call:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "HI_wings_region = SpectralRegion([(1.44*u.GHz, 1.43*u.GHz), (1.41*u.GHz, 1.4*u.GHz)])\n",
    "HI_wings_region"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While regions are useful for a variety of analysis steps, fundamentally they can be used to extract sub-spectra from larger spectra:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from specutils.manipulation import extract_region\n",
    "\n",
    "subspec = extract_region(sdss_spec, pixel_region)\n",
    "plt.step(subspec.wavelength, subspec.flux)\n",
    "\n",
    "analysis.snr(subspec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Line Measurements\n",
    "\n",
    "While line-fitting (detailed more below) is a good choice for high signal-to-noise spectra or when detailed kinematics are desired, more empirical measures are often used in the literature for noisier spectra or just simpler analysis procedures. Specutils provides a set of functions to provide these sorts of measurements, as well as similar summary statistics about spectral regions.  The [analysis part of the specutils documentation](https://specutils.readthedocs.io/en/latest/analysis.html) provides a full list and detailed examples of these, but here we demonstrate some example cases."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: these line measurements generally assume your spectrum is continuum-subtracted or continuum-normalized. Some spectral pipelines do this for you, but often this is not the case.  For our examples here we will do this step \"by-eye\", but for a more detailed discussion of continuum modeling, see the next section.  Based on the above plot we estimate a continuum level for the area of the SDSS spectrum around the H-alpha emission line, and use basic math to construct the continuum-normalized and continuum-subtracted spectra."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# estimate a reasonable continuum-level estimate for the h-alpha area of the spectrum\n",
    "sdss_continuum = 205*subspec.flux.unit\n",
    "\n",
    "sdss_halpha_contsub = extract_region(sdss_spec, ha_region) - sdss_continuum\n",
    "\n",
    "plt.axhline(0, c='k', ls=':')\n",
    "plt.step(sdss_halpha_contsub.wavelength, sdss_halpha_contsub.flux)\n",
    "plt.ylim(-50, 50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the continuum level identified, we can now make some measurements of the spectral lines that are apparent by eye - in particular we will focus on the H-alpha emission line. While there are techniques for identifying the line automatically (see the fitting section below), here we assume we are doing \"quick-look\" procedures where manual identification is possible. \n",
    "\n",
    "In the cell below, change the values for `LOWER` and `UPPER` to make a spectral region that just encompasses the  H-alpha line (the middle of the three lines). You may find it useful to change the values, re-run the cell, and change again to \"hone in\" on the right number."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "LOWER = 6562 * u.angstrom\n",
    "UPPER = 6575 * u.angstrom\n",
    "halpha_lines_region = SpectralRegion(LOWER, UPPER)\n",
    "\n",
    "plt.step(sdss_halpha_contsub.wavelength, sdss_halpha_contsub.flux)\n",
    "\n",
    "yl1, yl2 = plt.ylim()\n",
    "plt.fill_between([halpha_lines_region.lower, halpha_lines_region.upper], \n",
    "                 yl1, yl2, alpha=.2)\n",
    "plt.ylim(yl1, yl2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can now call a variety of analysis functions on the continuum-subtracted spectrum to estimate various properties of the line (you can see the full list of relevant analysis functions [in the analysis part of the specutils docs](https://specutils.readthedocs.io/en/stable/analysis.html#functions)):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analysis.centroid(sdss_halpha_contsub, halpha_lines_region)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analysis.fwhm(sdss_halpha_contsub, halpha_lines_region)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analysis.line_flux(sdss_halpha_contsub, halpha_lines_region)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Equivalent width, being a continuum dependent property, can either be computed directly from the spectrum if the continuum level is given, or measured on a continuum-normalized spectrum. The latter is mainly useful if the continuum is non-uniform over the line being measured."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analysis.equivalent_width(sdss_spec, sdss_continuum, regions=halpha_lines_region)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sdss_halpha_contnorm = sdss_spec / sdss_continuum\n",
    "analysis.equivalent_width(sdss_halpha_contnorm, regions=halpha_lines_region)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "Load one of the spectrum datasets you made in the overview exercises into this notebook (i.e., your own dataset, a downloaded one, or the blackbody with an artificially added spectral feature).  Make a flux or width measurement of a line in that spectrum directly.  Is anything odd?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.modeling import models\n",
    "\n",
    "wl = np.linspace(4000, 12000, 1024)*u.angstrom\n",
    "flux = models.BlackBody(temperature=5800*u.K)(wl)\n",
    "\n",
    "noise = np.random.randn(flux.size)*1e-6 * flux.unit\n",
    "bbspec = Spectrum1D(spectral_axis=wl, flux=flux+noise, uncertainty=StdDevUncertainty(1e-6*np.ones_like(flux)))\n",
    "\n",
    "gmodel = models.Gaussian1D(amplitude=bbspec.flux.mean(), mean=6000, stddev=50)\n",
    "gflux = gmodel(bbspec.spectral_axis.value)\n",
    "\n",
    "bbspec_w_gaussian = bbspec + Spectrum1D(spectral_axis=bbspec.spectral_axis, flux=gflux)\n",
    "plt.step(bbspec_w_gaussian.spectral_axis, bbspec_w_gaussian.flux)\n",
    "\n",
    "gregion = SpectralRegion(5900*u.angstrom, 6100*u.angstrom)\n",
    "plt.axvline(gregion.lower, c='k')\n",
    "plt.axvline(gregion.upper, c='k')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analysis.line_flux(bbspec_w_gaussian, gregion)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Is that right?  Lets compare to a flat chunk of the size of the line region\n",
    "200*3e-5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# It's too big! I.e., it includes the continuum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Continuum Subtraction\n",
    "\n",
    "While continuum-fitting for spectra is sometimes thought of as an \"art\" as much as a science, specutils provides the tools to do a variety of approaches to continuum-fitting, without making a specific recommendation about what is \"best\" (since it is often very data-dependent).  More details are available [in the relevant specutils doc section](https://specutils.readthedocs.io/en/latest/fitting.html#continuum-fitting), but here we outline the two basic options as it stands: an \"often good-enough\" function, and a more customizable tool that leans on the [`astropy.modeling`](http://docs.astropy.org/en/stable/modeling/index.html) models to provide its flexibility."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The \"often good-enough\" way\n",
    "\n",
    "The `fit_generic_continuum` function provides a function that is often sufficient for reasonably well-behaved continuua, particular for \"quick-look\" or similar applications where high precision is not that critical.  The function yields a continuum model, which can be evaluated at any spectral axis value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from specutils.fitting import fit_generic_continuum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "generic_continuum = fit_generic_continuum(sdss_spec)\n",
    "\n",
    "generic_continuum_evaluated = generic_continuum(sdss_spec.spectral_axis)\n",
    "\n",
    "plt.step(sdss_spec.spectral_axis, sdss_spec.flux)\n",
    "plt.plot(sdss_spec.spectral_axis, generic_continuum_evaluated)\n",
    "plt.ylim(100, 300);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Note that in some versions of astropy/specutils you may see a warning that the \"Model is linear in parameters\" upon executing the above cell. This is not a problem unless performance is a serious concern, in which case more customization is required.)\n",
    "\n",
    "With this model in hand, continuum-subtracted or continuum-normalized spectra can be produced using basic spectral manipulations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sdss_gencont_sub = sdss_spec - generic_continuum(sdss_spec.spectral_axis)\n",
    "sdss_gencont_norm = sdss_spec / generic_continuum(sdss_spec.spectral_axis)\n",
    "\n",
    "ax1, ax2 = plt.subplots(2, 1)[1]\n",
    "\n",
    "ax1.step(sdss_gencont_sub.wavelength, sdss_gencont_sub.flux)\n",
    "ax1.set_ylim(-50, 50)\n",
    "ax1.axhline(0, color='k', ls=':')  # continuum should be at flux=0\n",
    "\n",
    "ax2.step(sdss_gencont_norm.wavelength, sdss_gencont_norm.flux)\n",
    "ax2.set_ylim(0, 2)\n",
    "ax2.axhline(1, color='k', ls='--');  # continuum should be at flux=1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The customizable way\n",
    "\n",
    "The `fit_continuum` function operates similarly to `fit_generic_continuum`, but is meant for you to provide your favorite continuum model rather than being tailored to a specific continuum model. To see the list of models, see the [astropy.modeling documentation](http://docs.astropy.org/en/stable/modeling/index.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from specutils.fitting import fit_continuum\n",
    "from astropy.modeling import models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, suppose you want to use a 3rd-degree Chebyshev polynomial as your continuum model. You can use `fit_continuum` to get an object that behaves the same as for `fit_generic_continuum`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "chebdeg3_continuum = fit_continuum(sdss_spec, models.Chebyshev1D(3))\n",
    "\n",
    "generic_continuum_evaluated = generic_continuum(sdss_spec.spectral_axis)\n",
    "\n",
    "plt.step(sdss_spec.spectral_axis, sdss_spec.flux)\n",
    "plt.plot(sdss_spec.spectral_axis, chebdeg3_continuum(sdss_spec.spectral_axis))\n",
    "plt.ylim(100, 300);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This then provides total flexibility.  For example, you can also try other polynomials like higher-degree Hermite polynomials:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hermdeg7_continuum = fit_continuum(sdss_spec, models.Hermite1D(degree=7))\n",
    "hermdeg17_continuum = fit_continuum(sdss_spec, models.Hermite1D(degree=17))\n",
    "\n",
    "plt.step(sdss_spec.spectral_axis, sdss_spec.flux)\n",
    "plt.plot(sdss_spec.spectral_axis, chebdeg3_continuum(sdss_spec.spectral_axis))\n",
    "plt.plot(sdss_spec.spectral_axis, hermdeg7_continuum(sdss_spec.spectral_axis))\n",
    "plt.plot(sdss_spec.spectral_axis, hermdeg17_continuum(sdss_spec.spectral_axis))\n",
    "plt.ylim(150, 250);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This immediately demonstrates the tradeoffs in polynomial fitting: while the high-degree polynomials capture the wiggles of the spectrum better than the low, they also *over*-fit near the strong emission lines."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "Try combining the `SpectralRegion` and continuum-fitting functionality to only fit the parts of the spectrum that *are* continuum (i.e. not including emission lines).  Can you do better?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "LOWER = 6562 * u.angstrom\n",
    "UPPER = 6575 * u.angstrom\n",
    "halpha_lines_region_ex = SpectralRegion(LOWER, UPPER)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "better_hermdeg17_continuum = fit_continuum(sdss_spec, models.Hermite1D(degree=17), exclude_regions=halpha_lines_region_ex)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.step(sdss_spec.spectral_axis, sdss_spec.flux)\n",
    "plt.plot(sdss_spec.spectral_axis, hermdeg17_continuum(sdss_spec.spectral_axis))\n",
    "plt.plot(sdss_spec.spectral_axis, better_hermdeg17_continuum(sdss_spec.spectral_axis))\n",
    "plt.ylim(100, 300);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apparent that the fit is better around the Halpha line"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "Using the spectrum from the previous exercise, first subtract a continuum, then re-do your measurement.  Is it better?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "generic_continuum_bbspec = fit_generic_continuum(bbspec) # we cheat a smidge and fit on the line-less version\n",
    "bbspec_w_gaussian_contsub = bbspec_w_gaussian - generic_continuum_bbspec(bbspec_w_gaussian.spectral_axis)\n",
    "\n",
    "plt.step(bbspec_w_gaussian_contsub.spectral_axis, bbspec_w_gaussian_contsub.flux)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analysis.line_flux(bbspec_w_gaussian_contsub, gregion)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Much better number than before!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Line-Fitting\n",
    "\n",
    "In addition to the more empirical measurements described above, `specutils` provides tools for doing spectral line fitting. The approach is akin to that for continuum modeling: models from [astropy.modeling](http://docs.astropy.org/en/stable/modeling/index.html) are fit to the spectrum, and either those models can be used directly, or their parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from specutils import fitting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fitting machinery must first be given guesses for line locations. This process can be automated using functions designed to identify lines (more detail on the options is [in the docs](https://specutils.readthedocs.io/en/latest/fitting.html#line-finding)).  For data sets where these algorithms are not ideal, you may substitute your own (i.e., skip this step and start with line location guesses). \n",
    "\n",
    "Here we identify the three lines near the Halpha region in our SDSS spectrum, finding the lines above about a $\\sim 3 \\sigma$ flux threshold.  They are then output as an astropy Table:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "halpha_lines = fitting.find_lines_threshold(sdss_halpha_contsub, 3)\n",
    "\n",
    "plt.step(sdss_halpha_contsub.spectral_axis, sdss_halpha_contsub.flux, where='mid')\n",
    "for line in halpha_lines:\n",
    "    plt.axvline(line['line_center'], color='k', ls=':')\n",
    "\n",
    "halpha_lines_region"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(If you see a warning about the signal-to-noise, you can ignore it, or follow the instructions it gives to supress the warning. It is occurring because our cutout has a lot of real flux so it *could* be the case instead that we forgot to subtract the continuum.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now for each of these lines, we need to fit a model. Sometimes it is sufficient to simply create a model where the center is at the line and excise the appropriate area of the line to do a  line estimate.  This is not *too* sensitive to the size of the region, at least for well-separated lines like these.  The result is a list of models that carry with them them the details of the fit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "halpha_line_models = []\n",
    "for line in halpha_lines:\n",
    "    line_region = SpectralRegion(line['line_center']-5*u.angstrom,\n",
    "                                 line['line_center']+5*u.angstrom)\n",
    "    line_spectrum = extract_region(sdss_halpha_contsub, line_region)\n",
    "    # here's the workaround from above again\n",
    "    line_spectrum = Spectrum1D(flux=line_spectrum.flux, spectral_axis=line_spectrum.spectral_axis, uncertainty=line_spectrum.uncertainty)\n",
    "    line_estimate = models.Gaussian1D(mean=line['line_center'])\n",
    "    line_model = fitting.fit_lines(line_spectrum, line_estimate)\n",
    "    \n",
    "    halpha_line_models.append(line_model)\n",
    "    \n",
    "plt.step(sdss_halpha_contsub.spectral_axis, sdss_halpha_contsub.flux, where='mid')\n",
    "for line_model in halpha_line_models:\n",
    "    evaluated_model = line_model(sdss_halpha_contsub.spectral_axis)\n",
    "    plt.plot(sdss_halpha_contsub.spectral_axis, evaluated_model)  \n",
    "    \n",
    "halpha_line_models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more complicated models or fits it may be better to use the `estimate_line_parameters` function instead of manually creating e.g. a `Gaussian1D` model and setting the center.  An example of this pattern is given below.\n",
    "\n",
    "Note that we provided a default `Gaussian1D` model to the `estimate_line_parameters` function above.  This function makes reasonable guesses for `Gaussian1D`, `Voigt1D`, and `Lorentz1D`, the most common line profiles used for spectral lines, but may or may not work for other models.  See [the relevant docs section](https://specutils.readthedocs.io/en/latest/fitting.html#parameter-estimation) for more details.\n",
    "\n",
    "In this example we also show an example of a *joint* fit of all three lines at the same time.  While the difference may seems subtle, in cases of blended lines this typically provides much better fits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "halpha_line_estimates = []\n",
    "for line in halpha_lines:\n",
    "    line_region = SpectralRegion(line['line_center']-3*u.angstrom,\n",
    "                                 line['line_center']+3*u.angstrom)\n",
    "    line_spectrum = extract_region(sdss_halpha_contsub, line_region)\n",
    "    line_estimate = fitting.estimate_line_parameters(line_spectrum, models.Gaussian1D())\n",
    "    \n",
    "    halpha_line_estimates.append(line_estimate)\n",
    "\n",
    "# this could be done more flexibly with a for loop but we are explicit here for simplicity\n",
    "combined_model_estimate = halpha_line_estimates[0] + halpha_line_estimates[1] + halpha_line_estimates[2]\n",
    "combined_model_estimate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "combined_model = fitting.fit_lines(sdss_halpha_contsub, combined_model_estimate)\n",
    "\n",
    "plt.step(sdss_halpha_contsub.spectral_axis, sdss_halpha_contsub.flux, where='mid')\n",
    "plt.plot(sdss_halpha_contsub.spectral_axis, \n",
    "         combined_model(sdss_halpha_contsub.spectral_axis))  \n",
    "    \n",
    "combined_model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "Fit a spectral feature from your own spectrum using the fitting methods outlined above. Try the different line profile types (Gaussian, Lorentzian, or Voigt).  If you are using the blackbody spectrum (where you know the \"true\" answer for the spectral line), compare your answer to the true answer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.step(bbspec_w_gaussian.spectral_axis, bbspec_w_gaussian.flux);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "line_region = SpectralRegion(5800*u.angstrom, 6200*u.angstrom)\n",
    "line_spectrum = extract_region(bbspec_w_gaussian_contsub, line_region)\n",
    "line_estimate = fitting.estimate_line_parameters(line_spectrum, models.Gaussian1D())\n",
    "line_estimate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "line_model = fitting.fit_lines(line_spectrum, line_estimate)\n",
    "line_model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.step(line_spectrum.spectral_axis, line_spectrum.flux);\n",
    "model_wl = np.linspace(line_spectrum.spectral_axis[0], line_spectrum.spectral_axis[-1], 100)\n",
    "plt.plot(model_wl, line_model(model_wl))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bonus: we can fit the *combined* spectrum since it's a blackbody with a single line:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# amplitude of 1 for the bb is what we used to create it so that's a good estimate\n",
    "combined_estimated_model = line_estimate + models.BlackBody(temperature=5800*u.K) \n",
    "\n",
    "plt.step(bbspec_w_gaussian.spectral_axis, bbspec_w_gaussian.flux);\n",
    "model_wl = np.linspace(bbspec_w_gaussian.spectral_axis[0], bbspec_w_gaussian.spectral_axis[-1], 100)\n",
    "plt.plot(model_wl, combined_estimated_model(model_wl))\n",
    "combined_estimated_model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# notice this is using the *line spectrum* to focus on getting the line part right\n",
    "full_spec_model = fitting.fit_lines(line_spectrum, combined_estimated_model)\n",
    "full_spec_model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.step(bbspec_w_gaussian.spectral_axis, bbspec_w_gaussian.flux);\n",
    "model_wl = np.linspace(bbspec_w_gaussian.spectral_axis[0], bbspec_w_gaussian.spectral_axis[-1], 100)\n",
    "plt.plot(model_wl, full_spec_model(model_wl))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# huzzah!  the fit gives a pretty close answer given the noise\n",
    "gmodel, full_spec_model.left"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}