{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Emission Line Galaxy Classification\n",
    "\n",
    "[Chris Richardson, Elon University](https://facstaff.elon.edu/crichardson17/) \n",
    "\n",
    "<i>Funding for this tutorial made possible through the [CATL Scholars Program](https://www.elon.edu/u/academics/catl/sotl/catl-scholars-program/)</i>\n",
    "\n",
    "<i>This tutorial is intended for students to complete in groups. A version with solutions is available to instructors upon <a href=\"mailto:crichardson17@elon.edu\">request</a>.</i>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Objectives\n",
    "\n",
    "- To explore classifying emission line galaxies according to the radiation source excitating gaseous clouds within each galaxy\n",
    "- To understand how emission line ratios can help constrain excitation mechanisms\n",
    "- To explore how models can highlight the strengths and weakness of various diagnostic diagrams\n",
    "- To provide a foundation to explore your own questions about excitation mechanisms in galaxies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learning Outcomes\n",
    "\n",
    "By the end of this tutorial you should be able to:\n",
    "\n",
    "- Use emission line ratios to explore their diagnostic potential for finding AGN\n",
    "- Use diagnostic diagrams to demarcate AGN and SF boundaries\n",
    "- Use photoionization models to test the validity of diagnostic diagrams"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prerequisites\n",
    "\n",
    "- Computational\n",
    "    - Basic Python skills (e.g., for loops, conditional statements, defining functions)\n",
    "    - Basics of <code>matplotlib</code> visualization and SQL queries (e.g., see [here](https://escip.io/notebooks/sdss_tutorial_1.html))\n",
    "    - Using the <code>[astropy](https://docs.astropy.org/en/stable/units/index.html)</code> package for fundamental constants and unit conversion\n",
    "\n",
    "- Astrophysics\n",
    "    - Introductory knowledge of the formation of [spectral lines](https://openstax.org/books/astronomy-2e/pages/5-5-formation-of-spectral-lines)\n",
    "    - Basic understanding of the physical picture present in [star forming galaxies](http://skyserver.sdss.org/edr/en/proj/advanced/galaxies/spectra.asp) and [active galaxies](https://imagine.gsfc.nasa.gov/science/objects/active_galaxies1.html) (i.e., AGN)\n",
    "    - Advanced conceptual understanding of using emission lines to probe spectral energy distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import relevant libraries\n",
    "\n",
    "import pandas as pd # data analysis \n",
    "import numpy as np # more data analysis\n",
    "import seaborn as sns # plotting library\n",
    "from matplotlib import pyplot as plt # another plotting library\n",
    "from astroquery.sdss import SDSS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make notebook have a clean appearance\n",
    "\n",
    "pd.set_option('display.max_colwidth', None)\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Choosing a sample of galaxies\n",
    "\n",
    "[<code>GalSpecLine</code>](http://skyserver.sdss.org/dr16/en/help/browser/browser.aspx#&&history=description+galSpecLine+U) and [<code>emissionLinesPort</code>](http://skyserver.sdss.org/dr16/en/help/browser/browser.aspx#&&history=description+emissionLinesPort+U) are examples of two schema available to search for emission line galaxies. The most common emission lines are featured in both schema, but there are several lines not cross listed between the two catalogs. The methods used for measuring the emission line strengths differ between the catalogs, each with their unique merits. We will use [GalSpecLine](http://skyserver.sdss.org/dr16/en/help/browser/browser.aspx#&&history=description+galSpecLine+U) for this tutorial, however <code>emissionLinesPort</code> could be used in a follow up project. \n",
    "\n",
    "While you can include as many or as little emission lines as you'd like in your query, you should always include the emission line flux errors and the galaxy redshifts. The errors on the emission line measurements enable you to improve the quality of your sample if only specific emission lines interest you. It is better to impose these restrictions after your query as it can dramatically reduce the sample size. The redshift of each galaxy is important because certain emission lines will shift into and outside the SDSS window of observation.\n",
    "\n",
    "- At what $z$ will the shortest wavelength emission line shift into the observation window?\n",
    "- At what $z$ will the longest wavelength emission line shift outside the observation window?\n",
    "\n",
    "The SDSS aperture is fixed in size and is typically aligned with the nuclear center of each galaxy. \n",
    "\n",
    "- How might the size of the SDSS aperture tie into your ability to discern galaxy wide properties at certain redshifts?\n",
    "- What sort of bias does this introduce into your results?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Your answers to these questions:</b>\n",
    "<br>\n",
    "<br>\n",
    "<br>\n",
    "<br>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>plate</th>\n",
       "      <th>fiberid</th>\n",
       "      <th>mjd</th>\n",
       "      <th>z</th>\n",
       "      <th>zwarning</th>\n",
       "      <th>oii_3726_flux</th>\n",
       "      <th>oii_3729_flux</th>\n",
       "      <th>neiii_3869_flux</th>\n",
       "      <th>h_delta_flux</th>\n",
       "      <th>h_gamma_flux</th>\n",
       "      <th>...</th>\n",
       "      <th>oiii_5007_flux_err</th>\n",
       "      <th>hei_5876_flux_err</th>\n",
       "      <th>oi_6300_flux_err</th>\n",
       "      <th>nii_6548_flux_err</th>\n",
       "      <th>h_alpha_flux_err</th>\n",
       "      <th>nii_6584_flux_err</th>\n",
       "      <th>sii_6717_flux_err</th>\n",
       "      <th>sii_6731_flux_err</th>\n",
       "      <th>ariii7135_flux_err</th>\n",
       "      <th>bptclass</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>266</td>\n",
       "      <td>4</td>\n",
       "      <td>51630</td>\n",
       "      <td>0.064656</td>\n",
       "      <td>0</td>\n",
       "      <td>39.057380</td>\n",
       "      <td>27.096570</td>\n",
       "      <td>2.468194</td>\n",
       "      <td>6.669699</td>\n",
       "      <td>12.437470</td>\n",
       "      <td>...</td>\n",
       "      <td>3.302969</td>\n",
       "      <td>2.745098</td>\n",
       "      <td>2.520172</td>\n",
       "      <td>1.046200</td>\n",
       "      <td>3.448143</td>\n",
       "      <td>3.155578</td>\n",
       "      <td>2.919328</td>\n",
       "      <td>2.860605</td>\n",
       "      <td>4.992754</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>266</td>\n",
       "      <td>120</td>\n",
       "      <td>51630</td>\n",
       "      <td>0.069341</td>\n",
       "      <td>0</td>\n",
       "      <td>-1.330543</td>\n",
       "      <td>7.919315</td>\n",
       "      <td>-0.333586</td>\n",
       "      <td>-1.529107</td>\n",
       "      <td>10.477300</td>\n",
       "      <td>...</td>\n",
       "      <td>2.545394</td>\n",
       "      <td>1.778672</td>\n",
       "      <td>1.680469</td>\n",
       "      <td>0.703648</td>\n",
       "      <td>2.371513</td>\n",
       "      <td>2.122362</td>\n",
       "      <td>2.055846</td>\n",
       "      <td>1.841887</td>\n",
       "      <td>3.126379</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>266</td>\n",
       "      <td>146</td>\n",
       "      <td>51630</td>\n",
       "      <td>0.052176</td>\n",
       "      <td>0</td>\n",
       "      <td>145.169700</td>\n",
       "      <td>184.163200</td>\n",
       "      <td>26.046470</td>\n",
       "      <td>67.026570</td>\n",
       "      <td>114.044300</td>\n",
       "      <td>...</td>\n",
       "      <td>4.271465</td>\n",
       "      <td>3.850465</td>\n",
       "      <td>3.498755</td>\n",
       "      <td>2.040700</td>\n",
       "      <td>10.994270</td>\n",
       "      <td>6.155217</td>\n",
       "      <td>4.846725</td>\n",
       "      <td>4.428813</td>\n",
       "      <td>3.586765</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>266</td>\n",
       "      <td>533</td>\n",
       "      <td>51630</td>\n",
       "      <td>0.085351</td>\n",
       "      <td>0</td>\n",
       "      <td>31.840760</td>\n",
       "      <td>28.041630</td>\n",
       "      <td>2.403825</td>\n",
       "      <td>12.097370</td>\n",
       "      <td>18.920390</td>\n",
       "      <td>...</td>\n",
       "      <td>2.793131</td>\n",
       "      <td>1.910691</td>\n",
       "      <td>1.961116</td>\n",
       "      <td>0.822471</td>\n",
       "      <td>3.265500</td>\n",
       "      <td>2.480760</td>\n",
       "      <td>2.806020</td>\n",
       "      <td>2.475492</td>\n",
       "      <td>2.315306</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>266</td>\n",
       "      <td>630</td>\n",
       "      <td>51630</td>\n",
       "      <td>0.073687</td>\n",
       "      <td>0</td>\n",
       "      <td>12.878990</td>\n",
       "      <td>21.359510</td>\n",
       "      <td>3.436489</td>\n",
       "      <td>0.142290</td>\n",
       "      <td>5.390950</td>\n",
       "      <td>...</td>\n",
       "      <td>1.988206</td>\n",
       "      <td>1.574178</td>\n",
       "      <td>1.350679</td>\n",
       "      <td>0.568662</td>\n",
       "      <td>2.131675</td>\n",
       "      <td>1.715214</td>\n",
       "      <td>1.609602</td>\n",
       "      <td>1.671389</td>\n",
       "      <td>1.910809</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9995</th>\n",
       "      <td>1207</td>\n",
       "      <td>463</td>\n",
       "      <td>52672</td>\n",
       "      <td>0.055911</td>\n",
       "      <td>0</td>\n",
       "      <td>-0.273259</td>\n",
       "      <td>21.205140</td>\n",
       "      <td>0.332328</td>\n",
       "      <td>6.859152</td>\n",
       "      <td>4.029799</td>\n",
       "      <td>...</td>\n",
       "      <td>2.268165</td>\n",
       "      <td>2.104548</td>\n",
       "      <td>1.901967</td>\n",
       "      <td>0.799355</td>\n",
       "      <td>2.466196</td>\n",
       "      <td>2.411038</td>\n",
       "      <td>2.136849</td>\n",
       "      <td>2.025468</td>\n",
       "      <td>2.065840</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9996</th>\n",
       "      <td>1207</td>\n",
       "      <td>512</td>\n",
       "      <td>52672</td>\n",
       "      <td>0.040428</td>\n",
       "      <td>0</td>\n",
       "      <td>168.423800</td>\n",
       "      <td>208.261600</td>\n",
       "      <td>16.263230</td>\n",
       "      <td>33.077510</td>\n",
       "      <td>53.475750</td>\n",
       "      <td>...</td>\n",
       "      <td>3.645575</td>\n",
       "      <td>2.059486</td>\n",
       "      <td>1.995668</td>\n",
       "      <td>0.826124</td>\n",
       "      <td>5.905511</td>\n",
       "      <td>2.491776</td>\n",
       "      <td>2.573609</td>\n",
       "      <td>2.407138</td>\n",
       "      <td>1.849287</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9997</th>\n",
       "      <td>1208</td>\n",
       "      <td>202</td>\n",
       "      <td>52672</td>\n",
       "      <td>0.051229</td>\n",
       "      <td>0</td>\n",
       "      <td>75.635860</td>\n",
       "      <td>81.785600</td>\n",
       "      <td>4.387515</td>\n",
       "      <td>17.040500</td>\n",
       "      <td>28.414330</td>\n",
       "      <td>...</td>\n",
       "      <td>3.219081</td>\n",
       "      <td>2.728835</td>\n",
       "      <td>2.423338</td>\n",
       "      <td>1.162954</td>\n",
       "      <td>4.954031</td>\n",
       "      <td>3.507735</td>\n",
       "      <td>3.237844</td>\n",
       "      <td>3.007984</td>\n",
       "      <td>2.462693</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9998</th>\n",
       "      <td>1208</td>\n",
       "      <td>490</td>\n",
       "      <td>52672</td>\n",
       "      <td>0.093831</td>\n",
       "      <td>0</td>\n",
       "      <td>37.472160</td>\n",
       "      <td>44.705480</td>\n",
       "      <td>2.661830</td>\n",
       "      <td>16.227150</td>\n",
       "      <td>30.224240</td>\n",
       "      <td>...</td>\n",
       "      <td>2.448881</td>\n",
       "      <td>1.810508</td>\n",
       "      <td>1.965327</td>\n",
       "      <td>0.843513</td>\n",
       "      <td>4.036311</td>\n",
       "      <td>2.544227</td>\n",
       "      <td>2.512585</td>\n",
       "      <td>2.451001</td>\n",
       "      <td>2.777650</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9999</th>\n",
       "      <td>1208</td>\n",
       "      <td>593</td>\n",
       "      <td>52672</td>\n",
       "      <td>0.071304</td>\n",
       "      <td>0</td>\n",
       "      <td>71.641330</td>\n",
       "      <td>85.087980</td>\n",
       "      <td>5.179615</td>\n",
       "      <td>23.398060</td>\n",
       "      <td>36.458060</td>\n",
       "      <td>...</td>\n",
       "      <td>2.943066</td>\n",
       "      <td>2.206493</td>\n",
       "      <td>2.083630</td>\n",
       "      <td>1.008152</td>\n",
       "      <td>4.212832</td>\n",
       "      <td>3.040816</td>\n",
       "      <td>2.661236</td>\n",
       "      <td>2.490553</td>\n",
       "      <td>2.830445</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>10000 rows × 40 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "      plate  fiberid    mjd         z  zwarning  oii_3726_flux  oii_3729_flux  \\\n",
       "0       266        4  51630  0.064656         0      39.057380      27.096570   \n",
       "1       266      120  51630  0.069341         0      -1.330543       7.919315   \n",
       "2       266      146  51630  0.052176         0     145.169700     184.163200   \n",
       "3       266      533  51630  0.085351         0      31.840760      28.041630   \n",
       "4       266      630  51630  0.073687         0      12.878990      21.359510   \n",
       "...     ...      ...    ...       ...       ...            ...            ...   \n",
       "9995   1207      463  52672  0.055911         0      -0.273259      21.205140   \n",
       "9996   1207      512  52672  0.040428         0     168.423800     208.261600   \n",
       "9997   1208      202  52672  0.051229         0      75.635860      81.785600   \n",
       "9998   1208      490  52672  0.093831         0      37.472160      44.705480   \n",
       "9999   1208      593  52672  0.071304         0      71.641330      85.087980   \n",
       "\n",
       "      neiii_3869_flux  h_delta_flux  h_gamma_flux  ...  oiii_5007_flux_err  \\\n",
       "0            2.468194      6.669699     12.437470  ...            3.302969   \n",
       "1           -0.333586     -1.529107     10.477300  ...            2.545394   \n",
       "2           26.046470     67.026570    114.044300  ...            4.271465   \n",
       "3            2.403825     12.097370     18.920390  ...            2.793131   \n",
       "4            3.436489      0.142290      5.390950  ...            1.988206   \n",
       "...               ...           ...           ...  ...                 ...   \n",
       "9995         0.332328      6.859152      4.029799  ...            2.268165   \n",
       "9996        16.263230     33.077510     53.475750  ...            3.645575   \n",
       "9997         4.387515     17.040500     28.414330  ...            3.219081   \n",
       "9998         2.661830     16.227150     30.224240  ...            2.448881   \n",
       "9999         5.179615     23.398060     36.458060  ...            2.943066   \n",
       "\n",
       "      hei_5876_flux_err  oi_6300_flux_err  nii_6548_flux_err  \\\n",
       "0              2.745098          2.520172           1.046200   \n",
       "1              1.778672          1.680469           0.703648   \n",
       "2              3.850465          3.498755           2.040700   \n",
       "3              1.910691          1.961116           0.822471   \n",
       "4              1.574178          1.350679           0.568662   \n",
       "...                 ...               ...                ...   \n",
       "9995           2.104548          1.901967           0.799355   \n",
       "9996           2.059486          1.995668           0.826124   \n",
       "9997           2.728835          2.423338           1.162954   \n",
       "9998           1.810508          1.965327           0.843513   \n",
       "9999           2.206493          2.083630           1.008152   \n",
       "\n",
       "      h_alpha_flux_err  nii_6584_flux_err  sii_6717_flux_err  \\\n",
       "0             3.448143           3.155578           2.919328   \n",
       "1             2.371513           2.122362           2.055846   \n",
       "2            10.994270           6.155217           4.846725   \n",
       "3             3.265500           2.480760           2.806020   \n",
       "4             2.131675           1.715214           1.609602   \n",
       "...                ...                ...                ...   \n",
       "9995          2.466196           2.411038           2.136849   \n",
       "9996          5.905511           2.491776           2.573609   \n",
       "9997          4.954031           3.507735           3.237844   \n",
       "9998          4.036311           2.544227           2.512585   \n",
       "9999          4.212832           3.040816           2.661236   \n",
       "\n",
       "      sii_6731_flux_err  ariii7135_flux_err  bptclass  \n",
       "0              2.860605            4.992754         3  \n",
       "1              1.841887            3.126379         2  \n",
       "2              4.428813            3.586765         1  \n",
       "3              2.475492            2.315306         2  \n",
       "4              1.671389            1.910809         2  \n",
       "...                 ...                 ...       ...  \n",
       "9995           2.025468            2.065840         2  \n",
       "9996           2.407138            1.849287         1  \n",
       "9997           3.007984            2.462693         1  \n",
       "9998           2.451001            2.777650         1  \n",
       "9999           2.490553            2.830445         1  \n",
       "\n",
       "[10000 rows x 40 columns]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create a query to generate a sample of galaxies\n",
    "\n",
    "query = 'SELECT TOP 10000 s.plate, s.fiberid, s.mjd, s.z, s.zwarning, '\n",
    "query += 'g.oii_3726_flux, g.oii_3729_flux, ' #-- Gather possible emission lines of interest\n",
    "query += 'g.neiii_3869_flux, g.h_delta_flux, '\n",
    "query += 'g.h_gamma_flux, g.oiii_4363_flux, '\n",
    "query += 'g.h_beta_flux, g.oiii_4959_flux, '\n",
    "query += 'g.oiii_5007_flux, g.hei_5876_flux, '\n",
    "query += 'g.oi_6300_flux, g.nii_6548_flux, '\n",
    "query += 'g.h_alpha_flux, g.nii_6584_flux, '\n",
    "query += 'g.sii_6717_flux, g.sii_6731_flux, '\n",
    "query += 'g.ariii7135_flux, '\n",
    "query += 'g.oii_3726_flux_err, g.oii_3729_flux_err, ' #-- Gather emission line errors\n",
    "query += 'g.neiii_3869_flux_err, g.h_delta_flux_err, '\n",
    "query += 'g.h_gamma_flux_err, g.oiii_4363_flux_err, '\n",
    "query += 'g.h_beta_flux_err, g.oiii_4959_flux_err, '\n",
    "query += 'g.oiii_5007_flux_err, g.hei_5876_flux_err, '\n",
    "query += 'g.oi_6300_flux_err, g.nii_6548_flux_err, '\n",
    "query += 'g.h_alpha_flux_err, g.nii_6584_flux_err, '\n",
    "query += 'g.sii_6717_flux_err, g.sii_6731_flux_err, '\n",
    "query += 'g.ariii7135_flux_err, '\n",
    "query += 'e.bptclass '\n",
    "query += 'FROM GalSpecLine AS g, SpecObj AS s, GalSpecExtra AS e ' \n",
    "query += 'WHERE '\n",
    "query += 's.specobjid = g.specobjid '\n",
    "query += 'AND s.specobjid = e.specobjid '\n",
    "query += 'AND z BETWEEN 0.04 and 0.1 '#-- Lower limit needed for global parameters w/ small aperture\n",
    "query += 'AND h_alpha_flux > h_alpha_flux_err*5 ' #-- Make sure Balmer lines have adequate S/N\n",
    "query += 'AND h_beta_flux > h_beta_flux_err*5 '\n",
    "query += 'AND h_alpha_flux_err > 0 ' #-- Consistency check\n",
    "query += 'AND h_beta_flux_err > 0 '\n",
    "query += 'AND sigma_balmer < 509.55 ' #-- Exclude the broad line emitting galaxies\n",
    "query += 'AND s.class = \\'GALAXY\\' '  #-- Looking for galaxies, not stars\n",
    "query += 'AND s.zwarning = 0 '\n",
    "#print (query) # useful for debugging, remove first # to uncomment\n",
    "\n",
    "# send query to SDSS\n",
    "gals = SDSS.query_sql(query,data_release=16)\n",
    "\n",
    "gals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using emission lines to find active galaxies\n",
    "\n",
    "In the lecture-tutorial on emission line ratio diagnostics, you discovered how to find emission line ratios sensitive to the excitation source. In summary, ratios featuring emission lines with relatively different ionization potentials effectively probe the shape of the ionizing radiation field and therefore the excitation source (e.g., starlight, AGN).\n",
    "\n",
    "As a trial run, let's explore how well the emission line ratio [O II] $\\lambda$3727 / [Ne III] $\\lambda$3869 separates star forming (SF) galaxies and AGN. N.B. [O II] $\\lambda$3727 is actually the sum of the [O II] $\\lambda$3726 and [O II] $\\lambda$3729 doublet. \n",
    "\n",
    "- What are relavent ionization potentials for [O II] and [Ne III] line emission? \n",
    "- Would you expect AGN to display higher or lower [O II]/[Ne III] ratios compared to SF galaxies?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Your answers to these questions:</b>\n",
    "<br>\n",
    "<br>\n",
    "<br>\n",
    "<br>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAU50lEQVR4nO3df7DddX3n8eebJBIxsZDchE2JmHQFxR9bxDtikelmy8qm1m5AQVKppgUns2uLtNqV4DgYZwcn45a1dCi73BGbdAoBNsiPyVgrExrbikKTGGtCpAFFuHCXhCCNkVIM+94/zjfx5HKTnHxyzvd7bu7zMcOc7+/vO9/5kFc+35+RmUiSdKSOa7oASdL4ZIBIkooYIJKkIgaIJKmIASJJKjK56QKOxsDAQM6bN6/pMiRpXNm4ceOzmTnraLczrgNk3rx5bNiwoekyJGlciYgfdWM7nsKSJBUxQCRJRQwQSVKRcX0NRJJK/exnP2N4eJgXX3yx6VJ6ZurUqcydO5cpU6b0ZPsGiKQJaXh4mOnTpzNv3jwioulyui4z2bVrF8PDw8yfP78n+/AUlqQJ6cUXX2TmzJnHZHgARAQzZ87saQ/LAJE0YR2r4bFPr/98BogkqYjXQCQJWL68me1de+213HrrrUyaNInjjjuOm266iauuuoqRkRFe/epXA/CZz3yGiy66qLsFdoEBIqlj3f5Ltu7t95tvfetbrF27lk2bNnH88cfz7LPP8tJLLwFwyy23MDg42HCFh2aASFJDRkZGGBgY4PjjjwdgYGCg4YqOjNdAJKkh559/Pk8++SSnn346H/vYx/jGN76xf96ll17KmWeeyZlnnsmuXbsarPLgehYgEfHliNgREVvaps2IiPsiYnv1e1LbvKsj4tGIeCQi/lOv6pKkfjFt2jQ2btzI0NAQs2bN4pJLLmHlypVA6xTW5s2b2bx5MzNnzmy20IPoZQ9kJbBw1LRlwLrMPA1YV40TEW8GFgNvqda5MSIm9bA2SeoLkyZNYsGCBXzuc5/jhhtu4M4772y6pI71LEAy82+B50ZNXgSsqoZXARe0Tb8tM/81M38IPAq8s1e1SVI/eOSRR9i+ffv+8c2bN/P617++wYqOTN0X0U/OzBGAzByJiNnV9FOAb7ctN1xNe4WIWAosBTj11FN7WKqkiaSJO8D27NnDFVdcwfPPP8/kyZN5wxvewNDQUF/esjuWfrkLa6zHJXOsBTNzCBgCGBwcHHMZSRoP3vGOd/DAAw+8Yvr69evrL6ZA3XdhPRMRcwCq3x3V9GHgdW3LzQWerrk2SdIRqDtA7gWWVMNLgHvapi+OiOMjYj5wGvBQzbVJko5Az05hRcRqYAEwEBHDwGeBFcAdEXE58ARwMUBmbo2IO4CHgb3A72Xmy72qTZJ09HoWIJn5WweZdd5Blr8WuLZX9UiSussn0SVJRQwQSVKRfrmNV5Ka1dD73O+66y7e//73s23bNt70pjcB8NBDD/GpT32Kp556iunTpzNnzhxWrFjB2972NpYvX84XvvAFHn/8cWbPbj1KN23aNPbs2dPd+jtgD0SSGrR69WrOPfdcbrvtNgCeeeYZPvjBD/L5z3+e7du3s2nTJq6++moee+yx/esMDAxw3XXXNVXyfgaIJDVkz549fPOb3+Tmm2/eHyA33HADS5Ys4Zxzztm/3LnnnssFF1ywf/yyyy7j9ttv57nnRr8tql4GiCQ15O6772bhwoWcfvrpzJgxg02bNrF161bOOuusQ643bdo0LrvsMq6//vqaKh2bASJJDVm9ejWLFy8GYPHixaxevfoVy5x99tmcccYZXHnllQdM//jHP86qVavYvXt3LbWOxYvoktSAXbt2cf/997NlyxYigpdffpmIYMmSJWzatIlFixYB8OCDD7JmzRrWrl17wPonnngiH/rQh7jxxhubKB+wByJJjVizZg0f+chH+NGPfsTjjz/Ok08+yfz58zn//PNZuXLlAS9ZfOGFF8bcxic+8Qluuukm9u7dW1fZB7AHIklQ+/vcV69ezbJlyw6Y9oEPfIBbb72V22+/nauuuoqnnnqK2bNnMzAwwDXXXPOKbQwMDHDhhRfyxS9+sa6yDxCZ4/eN6IODg7lhw4amy5AmjF7/HVvn3+Hbtm3jjDPOqG+HDRnrzxkRGzNz8Gi37SksSVIRA0SSVMQAkTRhjedT+J3o9Z/Pi+jSMaKJb3qPZ1OnTmXXrl3MnDmTiLG+qj2+ZSa7du1i6tSpPduHASJpQpo7dy7Dw8Ps3Lmz6VJ6ZurUqcydO7dn2zdAJE1IU6ZMYf78+U2XMa55DUSSVMQAkSQVMUAkSUUMEElSEQNEklTEAJEkFTFAJElFDBBJUhEDRJJUxACRJBUxQCRJRQwQSVIRA0SSVMQAkSQVaSRAIuIPI2JrRGyJiNURMTUiZkTEfRGxvfo9qYnaJEmdqT1AIuIU4OPAYGa+FZgELAaWAesy8zRgXTUuSepTTZ3Cmgy8OiImAycATwOLgFXV/FXABc2UJknqRO0BkplPAX8MPAGMAP+cmV8HTs7MkWqZEWD2WOtHxNKI2BARG47lT1FKUr9r4hTWSbR6G/OBXwReExG/3en6mTmUmYOZOThr1qxelSlJOowmTmH9R+CHmbkzM38GfAU4B3gmIuYAVL87GqhNktShJgLkCeBdEXFCRARwHrANuBdYUi2zBLingdokSR2aXPcOM/PBiFgDbAL2At8BhoBpwB0RcTmtkLm47tokSZ2rPUAAMvOzwGdHTf5XWr0RSdI44JPokqQiBogkqYgBIkkqYoBIkooYIJKkIgaIJKmIASJJKmKASJKKGCCSpCIGiCSpiAEiSSpigEiSihggkqQiBogkqYgBIkkqYoBIkooYIJKkIgaIJKmIASJJKmKASJKKGCCSpCIGiCSpiAEiSSpigEiSihggkqQiBogkqYgBIkkqYoBIkooYIJKkIgaIJKmIASJJKtJIgETEiRGxJiK+HxHbIuJXImJGRNwXEdur35OaqE2S1JnJDe33euBrmXlRRLwKOAH4NLAuM1dExDJgGXBVQ/VJXbd8edMVSN1Vew8kIl4L/CpwM0BmvpSZzwOLgFXVYquAC+quTZLUuSZOYf0SsBP484j4TkR8KSJeA5ycmSMA1e/ssVaOiKURsSEiNuzcubO+qiVJB2giQCYDZwH/KzPfDvyU1umqjmTmUGYOZubgrFmzelWjJOkwmgiQYWA4Mx+sxtfQCpRnImIOQPW7o4HaJEkdqj1AMvP/Ak9GxBurSecBDwP3AkuqaUuAe+quTZLUuabuwroCuKW6A+sHwO/SCrM7IuJy4Ang4oZqkyR1oJEAyczNwOAYs86ruRRJUiGfRJckFTFAJElFDBBJUhEDRJJUxACRJBUxQCRJRQ55G29EXHOI2ZmZ/73L9UiSxonDPQfy0zGmnQB8FJgJGCCSNEEdMkAy87p9wxExHbgSuAy4DbjuYOtJko59h30SPSJmAJ8ALqX1nY6zMvPHvS5MktTfDncN5H8A7weGgLdl5p5aqpIk9b3D3YX1SeAXgc8AT0fE7uq/n0TE7t6XJ0nqV4e7BuJtvpKkMRkQkqQiBogkqYgBIkkq0tQXCSXpFZYvH9/bn2jsgUiSihggkqQiBogkqYgBIkkqYoBIkooYIJKkIgaIJKmIASJJKmKASJKKGCCSpCIGiCSpiAEiSSpigEiSihggkqQijQVIREyKiO9ExNpqfEZE3BcR26vfk5qqTZJ0eE32QK4EtrWNLwPWZeZpwLpqXJLUpxoJkIiYC/wG8KW2yYuAVdXwKuCCmsuSJB2BpnogfwJ8Cvh/bdNOzswRgOp39lgrRsTSiNgQERt27tzZ80IlSWOrPUAi4n3AjszcWLJ+Zg5l5mBmDs6aNavL1UmSOtXEN9HfDfzniHgvMBV4bUT8JfBMRMzJzJGImAPsaKA2SVKHau+BZObVmTk3M+cBi4H7M/O3gXuBJdViS4B76q5NktS5fnoOZAXwnojYDrynGpck9akmTmHtl5nrgfXV8C7gvCbrkSR1rp96IJKkccQAkSQVMUAkSUUMEElSEQNEklTEAJEkFTFAJElFDBBJUhEDRJJUxACRJBVp9FUmUj9ZvrzpCqTxxR6IJKmIASJJKmKASJKKeA1EzenWRQcvXkiNsAciSSpigEiSihggkqQiBogkqYgBIkkqYoBIkooYIJKkIgaIJKmIASJJKmKASJKKGCCSpCIGiCSpiAEiSSpigEiSivg6d5Xpp1eod62Wbm1HmhjsgUiSitQeIBHxuoj4m4jYFhFbI+LKavqMiLgvIrZXvyfVXZskqXNN9ED2Ap/MzDOAdwG/FxFvBpYB6zLzNGBdNS5J6lO1B0hmjmTmpmr4J8A24BRgEbCqWmwVcEHdtUmSOtfoNZCImAe8HXgQODkzR6AVMsDsBkuTJB1GY3dhRcQ04E7gDzJzd0R0ut5SYCnAqaee2rsCj1X9dPeUpHGtkR5IREyhFR63ZOZXqsnPRMScav4cYMdY62bmUGYOZubgrFmz6ilYkvQKTdyFFcDNwLbM/J9ts+4FllTDS4B76q5NktS5Jk5hvRv4MPC9iNhcTfs0sAK4IyIuB54ALm6gNklSh2oPkMz8e+BgFzzOq7MWSVI5n0SXJBUxQCRJRQwQSVIRA0SSVMQAkSQVMUAkSUUMEElSEQNEklTEAJEkFTFAJElFDBBJUhEDRJJUxACRJBVp7IuE0rFqwfrlXdnO+gXd2Y7UK/ZAJElF7IGoK9av7/0+Fizo8fa71HNQ/1q+/NjYR7+wByJJKmIPRONGHb0cSZ2zByJJKmKASJKKGCCSpCIGiCSpiAEiSSriXVjjxUS6uVzSuGAPRJJUxACRJBUxQCRJRQwQSVIRA0SSVMS7sHrNu6ckHaPsgUiSihggkqQifRcgEbEwIh6JiEcjYlnT9UiSxtZX10AiYhLwZ8B7gGHgHyLi3sx8uNnKes9vXagX/D67eqnfeiDvBB7NzB9k5kvAbcCihmuSJI0hMrPpGvaLiIuAhZn50Wr8w8DZmfn7bcssBZZWo28FttRe6JEbAJ5tuogOWGd3WWf3jIcaYfzU+cbMnH60G+mrU1hAjDHtgITLzCFgCCAiNmTmYB2FHQ3r7C7r7K7xUOd4qBHGV53d2E6/ncIaBl7XNj4XeLqhWiRJh9BvAfIPwGkRMT8iXgUsBu5tuCZJ0hj66hRWZu6NiN8H/hqYBHw5M7ceYpWheio7atbZXdbZXeOhzvFQI0ywOvvqIrokafzot1NYkqRxwgCRJBXp2wA53CtNouVPq/n/GBFndbpuzXVeWtX3jxHxQET8ctu8xyPiexGxuVu31RXWuCAi/rmqY3NEXNPpujXX+d/aatwSES9HxIxqXi3HstrXlyNiR0SM+QxSH7XNw9XZD23zcDX2S9s8XJ390jZfFxF/ExHbImJrRFw5xjLda5+Z2Xf/0bqA/hjwS8CrgO8Cbx61zHuBv6L17Mi7gAc7XbfmOs8BTqqGf31fndX448BAHxzLBcDaknXrrHPU8r8J3F/nsWzb168CZwFbDjK/8bbZYZ2Nts0Oa2y8bXZSZx+1zTnAWdXwdOCfevl3Z7/2QDp5pcki4C+y5dvAiRExp8N1a6szMx/IzB9Xo9+m9WxLnY7mePTVsRzlt4DVParlkDLzb4HnDrFIP7TNw9bZB22zk2N5MH11LEdpsm2OZOamavgnwDbglFGLda199muAnAI82TY+zCsPwsGW6WTdbjnSfV1OK/n3SeDrEbExWq9o6YVOa/yViPhuRPxVRLzlCNftho73FREnAAuBO9sm13EsO9UPbfNINdE2O9V02+xYP7XNiJgHvB14cNSsrrXPvnoOpM1hX2lyiGU6WbdbOt5XRPwHWv+Tnts2+d2Z+XREzAbui4jvV//SqbvGTcDrM3NPRLwXuBs4rcN1u+VI9vWbwDczs/1fhHUcy071Q9vsWINtsxP90DaPRF+0zYiYRivE/iAzd4+ePcYqRe2zX3sgnbzS5GDL1Pk6lI72FRH/DvgSsCgzd+2bnplPV787gLtodSFrrzEzd2fmnmr4q8CUiBjoZN0662yzmFGnCGo6lp3qh7bZkYbb5mH1Sds8Eo23zYiYQis8bsnMr4yxSPfaZx0XdgouBE0GfgDM5+cXc94yapnf4MALQQ91um7NdZ4KPAqcM2r6a4DpbcMP0HoTcRM1/ht+/lDpO4EnquPaV8eyWu4XaJ2Lfk3dx3JUHfM4+IXfxttmh3U22jY7rLHxttlJnf3SNqtj8xfAnxxima61z748hZUHeaVJRPyXav7/Br5K626CR4EXgN891LoN1nkNMBO4MSIA9mbrbZ0nA3dV0yYDt2bm1xqq8SLgv0bEXuBfgMXZalH9diwBLgS+npk/bVu9lmO5T0SspnV30EBEDAOfBaa01dl42+ywzkbbZoc1Nt42O6wT+qBtAu8GPgx8LyI2V9M+TesfC11vn77KRJJUpF+vgUiS+pwBIkkqYoBIkooYIJKkIgaIJKmIASJJKmKA6JgWEfMi4l/a7oknIuZGxD0RsT0iHouI6yPiVQdZd0s1vCAi1lbDl1Svu17btuyciPh6/Pz1419t20ZGxBVty94QEb9zBH+G34mIG6rh5RHxR9Xwyoi4qBq+JSKe2zcu1cEA0UTwWGaeCa1vIQBfAe7OzNOA04FpwLWdbiwzbwc+OmryQloPYAH8XWa+t23eDuDKsUKqWzLzUuDeXm1fGosBoonm14AXM/PPATLzZeAPgcuqN6mWWsiBb7NttxNYBywZPSMi/m1EfK16U+vfRcSbjqIGqVYGiCaatwAb2ydk622lTwBvKNlgREwC3piZDx9isRXAJ6tl2w0BV2TmO4A/Am4sqUFqQl++C0vqoWDsV1QfbHonzuaV31w4QGb+MCIeAj60f4etV26fA/yf6l1JAMcX1iDVzgDRRLMV+ED7hIh4La3XWD9WuM1fBzp5Qd7ngTXAvm9BHAc8v+/6jDTeeApLE8064ISI+AjsP/10HbAyM18o3OZ51XYPKTO/DzwMvK8a3w38MCIurmqJiPjlwhqk2hkgmlCqV4FfCFwcEduBfwJepPXK6yMWEbNoXZQf/dW3g7mWA789filweUR8l1bvqGff9Za6zde565hWfRd6bWa+tcvbXUDrovdtwNzMXNE+PTPf1839dVjTSlp/1jV171sTkz0QHeteBn6h/UHCoxURl9C6W+rHmfmX+8Kj8hLw1n0PEtYlIm4B/j2t3pRUC3sgkqQi9kAkSUUMEElSEQNEklTEAJEkFfn/q5q6cgcjxgwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Remove possible +/- inf values\n",
    "gals_flux = gals[(gals['neiii_3869_flux']>0) & (gals['oii_3726_flux']>0) & (gals['oii_3729_flux']>0) \n",
    "                 & (gals['oi_6300_flux']>0) & (gals['oiii_5007_flux']>0)]\n",
    "\n",
    "# Make sure relavent emission lines have well measured fluxes\n",
    "quality_gals = gals_flux[(gals_flux['neiii_3869_flux']>5.0*gals_flux['neiii_3869_flux_err']) &\n",
    "                (gals_flux['oii_3726_flux']>5.0*gals_flux['oii_3726_flux_err']) &\n",
    "                (gals_flux['oii_3729_flux']>5.0*gals_flux['oii_3729_flux_err']) &\n",
    "                (gals_flux['oi_6300_flux']>5.0*gals_flux['oi_6300_flux_err'])]\n",
    "\n",
    "# Split the sample into AGN and SF galaxies\n",
    "gals_sf = quality_gals[quality_gals['bptclass']==1]\n",
    "gals_agn = quality_gals[(quality_gals['bptclass']==4)]\n",
    "\n",
    "oii_doublet_sf = gals_sf['oii_3726_flux']+gals_sf['oii_3729_flux']\n",
    "oii_doublet_agn = gals_agn['oii_3726_flux']+gals_agn['oii_3729_flux']\n",
    "\n",
    "oii_neiii_sf = np.log10(oii_doublet_sf/gals_sf['neiii_3869_flux'])\n",
    "oii_neiii_agn = np.log10(oii_doublet_agn/gals_agn['neiii_3869_flux'])\n",
    "\n",
    "plt.hist(oii_neiii_sf,color='blue',alpha=0.5,label='SF')\n",
    "plt.hist(oii_neiii_agn,color='red',alpha=0.5,label='AGN')\n",
    "plt.xlim(0.0,2.0)\n",
    "#plt.title('[O II]/[Ne III] as a function of galaxy classification')\n",
    "plt.xlabel('[O II]/[Ne III]')\n",
    "plt.ylabel('N')\n",
    "plt.legend()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The histogram shows the number of galaxies in each bin as a function of galaxy classification. It is useful in telling us that the $N_{SF}$ > $N_{AGN}$ and that each distribution peaks at a differnt [O II]/[Ne III]. But, the actual location of those peaks is bit uncertain because the histogram is rather blocky.\n",
    "\n",
    "The Python package <code>matplotlib</code> will try to select a bin width that results in a visually appealing histogram. We can adjust the number of bins, and therefore the bin width, using the \"bins\" kwarg. However, histograms will always appear bumpy and depend upon the endpoints of the bins. One way around this is to use a Keneral Density Estimator (KDE) plot, which you can read about [here](https://towardsdatascience.com/histograms-vs-kdes-explained-ed62e7753f12). We can use the \"seaborn\" package to make publicaton quality KDE plots that show a smooth distribution of our data. The plot now shows the <i>density</i> of galaxies rather than the <i>number</i> of galaxies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEGCAYAAACgt3iRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA+AUlEQVR4nO3deXxU9fX4/9dJCIQAgiyCAgpqgCQsAcNOFVAprmi1FWrVtrbWWv10tbWf7+dnXap1bV0QFRXBVhatoqgguCCrIAEREEQQQQIoe5A1JDm/P94THWKWSTJ37syd83w85kFyl7knYTJn7ns5b1FVjDHGmLpI8TsAY4wxic+SiTHGmDqzZGKMMabOLJkYY4ypM0smxhhj6qye3wFEU8uWLbVDhw5+h2GMMQlj6dKlO1W1VV2fJ1DJpEOHDuTn5/sdhjHGJAwR2RSN57FmLmOMMXVmycQYY0ydWTIxxhhTZ4HqMzHGmNo4evQoBQUFHD582O9QPJOenk67du1IS0vz5PktmRhjkl5BQQFNmjShQ4cOiIjf4USdqrJr1y4KCgro2LGjJ9ewZi5jTNI7fPgwLVq0CGQiARARWrRo4emdlyUTY4yBwCaSMl7/fJ4lExFpLyKzRWSNiHwsIr+t4BgRkUdEZL2IrBCRXmH7hovI2tC+W7yK0xhjTN15eWdSDPxRVbOAfsBvRCS73DHnAZmhx3XA4wAikgo8FtqfDYyq4FyTrI4ehY8+gvfeg+3b/Y7GBNDJJ4NI9B4nnxzZde+66y5ycnLo3r07ubm5LF68mMGDB9O5c2dyc3PJzc3lv//9r7c/fC151gGvqtuAbaGvvxaRNUBbYHXYYSOA59St0LVIRJqJyIlAB2C9qm4AEJHJoWPDzzXJ5uhRePBBeOghaNQImjSBDRtg+HD417+gbVu/IzQBsXkzzJ4dvecbMqT6Y95//31ef/11li1bRoMGDdi5cydFRUUAPP/88+Tl5UUvIA/EZDSXiHQAegKLy+1qC2wO+74gtK2i7X0ree7rcHc1nBxp+jeJp6AALrwQMjLgnnugrAbboUMwcSL06wczZ0K23cCaxLRt2zZatmxJgwYNAGjZsqXPEdWM5x3wItIYeAn4naruK7+7glO0iu3f3ag6VlXzVDWvVas61yoz8WjtWujb1yWMu+76NpEANGwI114LV17p7lB27fItTGPqYtiwYWzevJlOnTpxww03MGfOnG/2XXnlld80c+2K09e4p3cmIpKGSyTPq+rLFRxSALQP+74dsBWoX8l2k2w2bYKzz4af/ATOO6/y44YPd8f++Mfw5puuodqYBNK4cWOWLl3KvHnzmD17NldccQX33HMPkBjNXF6O5hLgGWCNqv6zksOmAVeHRnX1AwpDfS1LgEwR6Sgi9YGRoWNNMvn6azj/fBgxoupEUuYXv4DPP4eXXvI+NmM8kJqayuDBg7n99tsZPXo0LyXQa9nLZq6BwFXAUBFZHnqcLyLXi8j1oWOmAxuA9cBTwA0AqloM3AjMBNYAL6jqxx7GauKNKlx9NZx6Klx+eWTnpKbCDTfA73/v+lKMSSBr165l3bp133y/fPlyTjnlFB8jqhkvR3PNp+K+j/BjFPhNJfum45KNSUaPPQZr1sAjj9SsySo3Fzp2hGeegRtv9Cw8E2zt20c2Aqsmz1ed/fv3c9NNN7F3717q1avH6aefztixY7k80g9TPrPaXCb+rF4Nf/sbPPww1K9f8/NHjoR//AN+9SvwqKidCbYvvoj9Nc844wwWLlz4ne3vvfde7IOpBSunYuLL0aNuZNbPfgbt2tXuObKzoU0bmDIlurEZYyplycTElwcegAYN4IIL6vY8F10Ejz8enZiMMdWyZi4TP9atg/vvhzFj6j60d8AAePRR1++SlRWd+EyFVq+GadNg+XLXKtm3L/zwh3DCCX5HZmLJ7kxMfFB1fRyjRrkmqrqqV8/NPRk7tu7PZSpUUABXXAGDB8OSJW4uaevWLrF07uxyeWmp31GaWLE7ExMfXnzRvTv9v/8Xvec85xy45RbXdJaaGr3nNcye7cY5DB8OEya4QgRlLrjA1ba6915YtAiee85+/cnA7kyM/w4ccHNDbrwxuu86p5wCTZvC/PnRe07DK6+4Zqw//9mNkwhPJGXat3c1OdeudTecWmExJBMklkyM/+65x43A6t49+s995pmuEKSJijfecIUG7r4bzjij6mMbNIA77oB58+Dpp2MTX9T4VIN+6tSpiAiffPLJN9s++OADBg8eTGZmJr169eKCCy5g5cqVANx2221kZGSwPWwphsaNG0f3dxEha+Yy/tq8GUaPhief9Ob5hwyBm25ynfrW1lInK1fCNde4BNGpU2TnZGTAX/8Kf/yj61vJzPQ0xOjxowY9MGnSJAYNGsTkyZO57bbb+Oqrr/jRj37ExIkTGTBgAADz58/ns88+o1u3boCrLvzggw9y7733Ri/eWrA7E+OvW26Biy/2bujPiSdCixau8d7U2p49brT19dfXvMp/hw6u/ua111pzV1X279/PggULeOaZZ5g8eTIAo0eP5pprrvkmkQAMGjSISy655Jvvf/7znzNlyhR2794d65CPYcnE+GfZMnjrLTckyEv9+rmGflMrqvDLX0Lv3m5MQ21ccgls2QKvvRbV0ALllVdeYfjw4XTq1InmzZuzbNkyPv74Y3r16lXleY0bN+bnP/85Dz/8cIwirZglE+Ofm292peUzMry9zoAB8Oqr3l4jwMaPhxUrXEKprdRUuO46+MMfXJED812TJk1i5MiRAIwcOZJJkyZ955i+ffuSlZXFb3/722O2/8///A8TJkxg377yS0bFjiUT44+334bPPnMl5r2WmQl798L69d5fK2C2bYM//Qn+8pfalUkL16cPHH88VPAemfR27drFu+++yy9+8Qs6dOjA/fffz5QpU8jJyWHZsmXfHLd48WLuvPNOCgsLjzm/WbNm/PjHP2bMmDGxDv0blkxM7Km6caXXXOMmF3otJcW10cyc6f21Auamm1y+P+20uj+XCPzoR64Gp01mPNZ///tfrr76ajZt2sTGjRvZvHkzHTt2ZNiwYYwfP/6YApAHDx6s8Dn+8Ic/8OSTT1JcXByrsI9ho7lM7L38sltv5KyzYnfNM86A6dPhNxWueGAqMGsWLF4c3WG9eXludYAZM+pefs1TMa5BP2nSJG655ZZjtl122WVMnDiRKVOm8Je//IUtW7Zwwgkn0LJlS2699dbvPEfLli259NJL+de//hW9uGtANEDDK/Ly8jQ/P9/vMExVSkogJwd++lPXMR4rhYVw1VWwc2fd22uSQHExdO3qurQGDYruc7/zDsyZA3PnRvd562LNmjVkJUENt4p+ThFZqqp1XhPYy2V7x4nIdhFZVcn+m8NWYFwlIiUi0jy0b6OIrAzts+wQJJMnu9lsffvG9rpNm7qJYxWsF2G+68knoUkTGDgw+s995pmuOOSnn0b/uY1/vOwzGQ8Mr2ynqt6vqrmqmgv8FZijquEDpYeE9tc5Y5o4UVwMt97qluOta1Xg2ujZ03X8myrt3w+33+7KoHjx35SWBsOGWQ3OoPEsmajqXCDSWTSjABvjEXQTJ0KzZlDNuHnP5Oa6eS2mSo884n5Vp5/u3TXOP98ViDxyxLtr1FSQmvwr4vXP5/toLhHJwN3BvBS2WYFZIrJURK6r5vzrRCRfRPJ37NjhZaimLoqL4bbbXL+FH3cl4DoBVq1yH71NhfbudQUaf/ITb6/Trp2rwxkvkxjT09PZtWtXYBOKqrJr1y7S09M9u0Y8jOa6CFhQrolroKpuFZETgLdE5JPQnc53qOpYYCy4DnjvwzW1MnGim2SQm+tfDA0auIWy5s93tdPNdzz8sOvOirAuYZ0MHQr/+Q9cfrn316pOu3btKCgoIMgfSNPT02lX26WwIxAPyWQk5Zq4VHVr6N/tIjIV6APE0dgPUyMlJa4R/sYb/Y7EVSZ++21LJhXYv981ccVqZOmZZ8ITT7iBdk2bxuaalUlLS6Njx47+BpHgfG3mEpGmwFnAq2HbGolIk7KvgWFAhSPCTIKYMsUNDfLzrqRM9+7w3nt+RxGXnnjC/RfF4q4EoHFjNyZi6tTYXM94y8uhwZOA94HOIlIgIteKyPUicn3YYZcCs1T1QNi21sB8EfkI+AB4Q1Xf9CpO47HSUlez/Mc/9q+vJFxWllsX3vpNjlFU5PpKQqWhYmbwYLcSo0l8njVzqeqoCI4ZjxtCHL5tA9DDm6hMzE2d6kqm9O7tdyROgwZuMY5Fi2pfAjeApkxxdySxXm+kf3/45z9h925o3jy21zbR5ftoLhNgqu6uZNSo+LgrKZOTE1/Tr32mCvffDz/4QeyvnZ7uKt28/nrsr22iy5KJ8c6bb7oaXP37+x3Jsbp3j+4qegnuvffgwAH/bh779YOXXqr+OBPfLJkY79xxh2uET4mzl1l2Nnz4oS2sEfLQQzBihH//Tf37u9xeSTFckyDi7K/cBMa8eW5pvVhWBo5U48bQtq1LKElu82bX4nfuuf7F0LQpdO7sqhSbxGXJxHjjzjvhhz90S+zFo+xsN3kxyT3+OJx9NjRs6G8cffvGz2x4UzuWTEz0LV8OH30E3/++35FULjs76Tvhi4rcWiUXX+x3JC6ZzJjhBgOYxGTJxETf3Xe7oUHxvG5It26uHH0Sv3tNm+aGA8dqkmJV2rVzI8hXrvQ7ElNblkxMdG3Y4MqVXHih35FUrXVrN1x5wwa/I/HN44/HT1UZETea7I03/I7E1JYlExNdDzzg1mNt1MjvSKom4uabvP++35H4YsMGN/7gzDP9juRbvXvbfJNEZsnERM/Ona468KWX+h1JZLp0gQUL/I7CF08/7QoAxFNLZG6u62rbt8/vSExtWDIx0TN6tBsKnCh1MXJykjKZlJTA+PFw3nl+R3Ks9HT3X2J1OBOTJRMTHYcOwWOPwWWX+R1J5DIzYf16N/07icya5fJ9PFZcz821+SaJypKJiY7nnnMVeeNhaFCk6td3CeWDD/yOJKaeftrfSYpV6dXLVlZOVJZMTN2VlvpXKbCuunRJqk74Xbvcm/XQoX5HUrHMTNi+HQoK/I7E1JQlE1N3b7zhPuX3SMCVA7KykqrfZPJkV1ixSRO/I6lYSoqrIvzOO35HYmrKkompuwcecCO44qnMfKSys10zV5JMXhw3Lv6XcenRw/pNEpElE1M3H30En3zilsxLRK1auY/Dn3/udySeW7PGFXY84wy/I6laz55uRFeS5PfA8HLZ3nEisl1EKly/XUQGi0ihiCwPPW4N2zdcRNaKyHoRucWrGE0UPPywK+6UluZ3JLVTNnlx0SK/I/Hc+PGuqGO81t4s07YtFBcnRX4PFC/vTMYD1RVrmKequaHHHQAikgo8BpwHZAOjRCTbwzhNbe3c6VY1uuACvyOpm86dA99vUloK//kPDBvmdyTVE3FDhG2+SWLxLJmo6lxgdy1O7QOsV9UNqloETAZGRDU4Ex1PPw0DB0KzZn5HUjfZ2YEf0fXee67TPR7nllSkWzfrhE80fveZ9BeRj0RkhojkhLa1BTaHHVMQ2lYhEblORPJFJH/Hjh1exmrClZTAmDFuib5El5np+n0OHfI7Es+UNXElirI7E+s3SRx+JpNlwCmq2gN4FHgltL2iIUGVvqRUdayq5qlqXqtWraIfpanYjBlw3HGuiSjRpae7j+wBXXnx4EF49dX4nVtSkbJ+kyQu6pxwfEsmqrpPVfeHvp4OpIlIS9ydSPuwQ9sBW30I0VTl0Ufjv8x8TXTpEthO+FdecS15LVr4HUnkRKB7d7f6s0kMviUTEWkj4iYmiEifUCy7gCVApoh0FJH6wEhgml9xmgps2uTmZgwZ4nck0RPgCsKJ1sRVJjsb5szxOwoTKS+HBk8C3gc6i0iBiFwrIteLyPWhQy4HVonIR8AjwEh1ioEbgZnAGuAFVf3YqzhNLTz1lHt3atDA70iiJzs7kHcmX34JixfDoEF+R1Jz3bsn/crKCaWeV0+sqqOq2T8aGF3JvunAdC/iMnVUXAzPPAN33eV3JNF10kmuA76gwK0hGxATJ7pEkp7udyQ117GjG33+5ZfQpo3f0Zjq+D2ayySaWbPcrPFTT/U7kugqm7y4eLHfkUTVhAmJ1fEeLiXF3Z3Mn+93JCYSlkxMzTz1VPzWL6+rgE1eXLXKVeDNzfU7ktqzfpPEYcnERG7HDjeTLEgd7+Gys2HhQr+jiJqyu5J4L59Sla5dbURXorBkYiI3eTL07w+NG/sdiTe6dIEVK6CoyO9I6qykxJVPifcKwdXp3BnWrYP9+/2OxFTHkomJ3IQJiTnGNFIZGdC+PSxf7nckdfbuu/G7NG9NlC2GuWSJ35GY6lgyMZFZu9bNL4n3+uV1FZA6XePGBSfvZ2VZJ3wisGRiIvP8866vJJEb4CPRpUvCN9Lv2+cWvwxKMsnJsfkmicCSiameqksmiTrGtCa6dk34O5MXX3Q3kE2b+h1JdOTkuGau0lK/IzFVsWRiqvfhh3D0aDCKOlbnpJPg8GE3eTFBPfNMcO5KAI4/3iXG1av9jsRUxZKJqd7EiW5Z3kRc472mRNzdSYIOEV63Dj79FPr29TuS6MrJSfgbxsCzZGKqpgovvABnneV3JLGTlZWw/SbPPOPmlCbqKsqV6dLFOuHjnSUTU7Vly9yn9aCVT6lKTk5CJpPiYjd6OxGW5q2pgAyyCzRLJqZqL74I3/tecjRxlUnQmXJlZdMSfW5JRTp2hK1bYc8evyMxlbFkYiqn+m0ySSZlM+U++MDvSGrkySeDWzYtNdW1PgasDmegWDIxlVuzxq352qmT35HEXnZ2QjV1bdvm1kwP0iiu8jp3TthxEUnBkomp3CuvwIABydXEVSYnJ6HK1Y4b58ZIZGT4HYl3srOtEz6eWTIxlZs61RV2TEZdu7qZcsXFfkdSrdJSGDsWLrjA70i8lZ0NS5fa5MV45eWyveNEZLuIrKpk/5UisiL0WCgiPcL2bRSRlSKyXETyvYrRVOHLL92EhR49qj82iJo2dcv7ffih35FUa9Ysd0cS9NbIZs3cf8vatX5HYiri5Z3JeGB4Ffs/B85S1e7AncDYcvuHqGququZ5FJ+pyvTp0Lt38CYs1ESCLKbx6KNw4YXJ0RppnfDxy7Nkoqpzgd1V7F+oqmUD/RYBwVl4Owheew3ykjyPd+sGs2f7HUWVNm50ndJB7ngP16mTdcLHq3jpM7kWmBH2vQKzRGSpiFxX1Ykicp2I5ItI/o4dOzwNMmkUFbkFMYJWk6OmunVzPb5x3Ej/xBNuOHB6ut+RxEZWFixa5HcUpiK+JxMRGYJLJn8J2zxQVXsB5wG/EZEzKztfVceqap6q5rVq1crjaJPEggVukajjj/c7En+1auUa6VdV2O3nu0OH4Omn4aKL/I4kdk4/HdavhwMH/I7ElOdrMhGR7sDTwAhV3VW2XVW3hv7dDkwF+vgTYZJ64w3XX2Kge/e4beqaPNnNvWjf3u9IYqd+fZdQli71OxJTnm/JREROBl4GrlLVT8O2NxKRJmVfA8OA+PxoGFRvvmn9JWV69IC33/Y7iu9QhX/+E0aM8DuS2Ovc2Zq64lG9SA4SkZeAccAMVY2oAVlEJgGDgZYiUgD8DUgDUNUngFuBFsAYccNQikMjt1oDU0Pb6gETVfXNGvxMpi62bXNreXTp4nck8aFHDxg92vWbpPjeKvyNuXNd6bBkzPmdO1vRx3gUUTIBHgd+BjwiIi8C41X1k6pOUNVR1ez/BfCLCrZvAJJ0ckMcmDXLLdMX9OV5I9Wypes7Wr4cevXyO5pv3H8/XHppXOW3mMnKgmef9TsKU15EL0VVfVtVrwR6ARuBt0ITDX8mIkk8ESGAZsyAnj39jiK+5ObGVVPXp5+6T+ZBLDUfiZNOciXjtm71OxITLuLPNSLSAvgp7m7iQ+BhXHJ5y5PITOypuiHBydh2UpWePd0dW5x48EE3STFZhgOXJ+JKqyRYUefAiyiZiMjLwDwgA7hIVS9W1SmqehPQ2MsATQx9/LF7h2rTxu9I4kturuvxPXLE70jYvh2mTEnOjvdwmZk2Ez7eRHpn8rSqZqvqP1R1G4CINACwcicB8vbbcdUvEDeaNIEOHeKi1/fRR2HwYGje3O9I/JWVZTPh402kyeTvFWzz/y/LRNfMmclb2LE6cdDUtX8/jBkDl1/uaxhxoUsXV4MzjosTJJ0qk4mItBGRM4CGItJTRHqFHoNxTV4mKIqL3Uc963yvWO/ervilj556yuX6dlbFjqZNXRXhT6ocU2piqbqhwd/Hdbq3A/4Ztv1r4H89isn4YdkyaN3a/YWa78rOhg0b4Kuv3O8pxoqK4IEH4NZbY37puJWV5Trhs7P9jsRANXcmqjpBVYcAP1XVIWGPi1X15RjFaGLh3Xetiasq9eq5+Tdv+TN48T//cWVTOnf25fJxKTMzLrqxTEh1zVw/CX3ZQUT+UP4Rg/hMrLz9tqtDZSrXq5crzR9jJSVw991wxRUxv3RcswrC8aW6DvhGoX8bA00qeJggOHrUjbO0ZFK1fv1cJ3yMl/J96SW3kmJubkwvG/cyM2HdOlc92fivyj4TVX0y9O/tsQnH+GLpUjetuGlTvyOJb61auTk4778P3/teTC6pCnfeCT/+cXKspFgT9etDx45uVNeAAX5HYyKdtHifiBwnImki8o6I7AxrAjOJbs4ctxCUqV6fPvDqqzG73BtvuM73fv1idsmE0qmTTV6MF5HOMxmmqvuAC4ECoBNws2dRmdiaPduSSaT69YNp02JyKVW44w4YOdLuSipjFYTjR6TJpKyY4/nAJFWtdG13k2BKStxfoyWTyHTuDHv3wtq1nl9q9mzYsQPOrHSdUVM2PNj4L9Jk8pqIfALkAe+ISCvgsHdhmZhZuRJatLD6HJFKSYGBA2HqVM8vdfvtbgSXrQZQufbtYdcul3SNvyItQX8L0B/IU9WjwAEgyUvNBcTcudC1q99RJJaBA90QKw8tWgSffQbnnOPpZRJeSoqbtLhkid+RmJosrZMFXCEiVwOX45bTNYluzhxLJjWVmwvr17sVKT1yxx3wox+5uZKmapmZNt8kHkQ6muvfwAPAIKB36FFltWARGSci20WkwvXbxXlERNaLyAoR6RW2b7iIrA3tuyXin8bUjCosWGD9JTVVr54bi+rR3clHH0F+Pgwf7snTB06XLtYJHw8ivTPJAwaq6g2qelPo8T/VnDMeqOrP4TwgM/S4Drc0MCKSCjwW2p8NjBIRq77jhc8/dwnF1i+pue99DyZN8uSp//53uOwyN4/CVC8ryyVfVb8jSW6RJpNVQI3ecVR1LlDVqK8RwHPqLAKaiciJQB9gvapuUNUiYDLWP+ONefPcXYmNO625M85wI7qi3NS1fj288w5cdFFUnzbQWrSAhg3d7874J9Jk0hJYLSIzRWRa2aOO124LbA77viC0rbLtFRKR60QkX0Tyd9iQjpqZO9d9rDM1l5bmOuKnTInq0/7jHy6RZNgCDzWSlWWTF/0WaTK5DbgEuBt4MOxRFxV9HNYqtldIVceqap6q5rVq1aqOISWZ+fOtv6QuhgyBf/87ak+3davrhvnBD6L2lEmjUyfrN/FbpEOD5wAbgbTQ10uAZXW8dgHQPuz7dsDWKrabaNq9G7ZsgdNO8zuSxJWb6zJAlFZoevBBOPdcK5FWG1lZlkz8Fulorl8C/wWeDG1qC7xSx2tPA64OjerqBxSG1pdfAmSKSEcRqQ+MDB1romnhQsjJsRlxdZGaGrW7k717Ydw4W5K3tjp1cjn9sE2l9k2kzVy/AQYC+wBUdR1wQlUniMgk3DrxnUWkQESuFZHrReT60CHTgQ3AeuAp4IbQcxcDNwIzgTXAC6r6cY1+KlO9BQvcmEpTN+eeC889V+fFyMeMgf79fVnEMRDS0+GUU2D5cr8jSV6RTok6oqpFEhr1IyL1qKIfA0BVR1WzX3FJqqJ903HJxnhl7lwYYYPk6uz006FRI3jvPRg6tFZPcfgwPPww3HNPdENLNl26uE54q7Dsj0jvTOaIyP8CDUXkXOBFIPZLzpnoOHrUfYSzxbOj45xz4Jlnan36v//tclLHjlGMKQl17uxuuI0/Ik0mtwA7gJXAr3B3Df/nVVDGY8uXQ9u20Lix35EEw9lnu+V8CwtrfGppKdx7L/zwhx7ElWSys214sJ8iHc1Viutwv0FVL1fVp0LNVCYRLVxo80ui6fjjIS+vVjPiX3sNGjSAHj08iCvJtGvn8vlXX/kdSXKqMpmERlrdJiI7gU+AtSKyQ0RujU14xhPz51syibbhw+HJJ6s/rpx773UjuKwIQd2lpLgBilb00R/V3Zn8DjeKq7eqtlDV5kBfYKCI/N7r4IxHFi60SsHRdsYZsH07LIt8+tXixfDFF7b4VTTZ5EX/VJdMrgZGqernZRtUdQPwk9A+k2gKCuDQIddnYqInNRXOPx8efzziU+67Dy691Kb6RFNWlrvxNrFXXTJJU9Wd5Teq6g6+XcrXJJL333d3JdauEn3Dh8MLL8C+fdUe+vnn8O67cN55MYgriWRlufElxcV+R5J8qksmRbXcZ+KVTVb0TosWriP+ueeqPfRf/3KJxAo6RleTJm7i58qVfkeSfKpLJj1EZF8Fj68BqxCYiObPt/klXrrwQhg9usrFNfbudfnm0ktjF1Yyyc62fhM/VJlMVDVVVY+r4NFEVa2ZK9EcPgyrV7vZXcYbubmujWXOnEoPeeop6NsXrMi1N7p0cUv1mNiqyRrwJtEtWwYdOriVhIw3RNzdycMPV7j76FG3y8rMe6drV7sz8YMlk2Ty/vs2vyQWhg2D2bNh8+bv7HrpJdembzeH3mnf3jUlfvml35EkF0smyWTePOt8j4WMDFdipdwwYVW4/363vrvxTkqK3Z34wZJJslB1U4NtsmJsjBjhOkfCFtiYPx927bKqtrHQubPNN4k1SybJ4osvoKQE2rTxO5LkcPLJbhXLsDXi77vP9ZWk2F+d53JyrBM+1uxlnSzKVla0yYqxM2IEPPQQqLJunZviM2yY30Elh6wsWLXKFXswseFpMhGR4SKyVkTWi8gtFey/WUSWhx6rRKRERJqH9m0UkZWhfflexpkUbLJi7PXt69q13n+ff/7TDfJKT/c7qOTQsCGceirk2ztHzHiWTEQkFXgMOA/IBkaJyDGz5VT1flXNVdVc4K/AHFXdHXbIkND+PK/iTBoLFrg7ExM7KSlw8cUcue8hJk60hS1jzZq6YsvLO5M+wHpV3aCqRcBkoKo/p1FAzReEMNU7cADWrnUlVU1sDR+OvjmLS/pspUULv4NJLjk5Vc4dNVHmZTJpC4QPtC8IbfsOEckAhgMvhW1WYJaILBWR6zyLMhnk57t1YRs08DuSpFNUvzHv6BB+1/AJv0NJOl27ujL/paV+R5IcvEwmFfX0Vlaw6CJgQbkmroGq2gvXTPYbEalw1QcRuU5E8kUkf8eOHXWLOKgWLLDJij6ZNQs+bD+CnPljkaNWGzWWmjeHZs3g44/9jiQ5eJlMCoD2Yd+3A7ZWcuxIyjVxqerW0L/bgam4ZrPvUNWxqpqnqnmtrNhRxay4oy9KStxKvl2Gd+DwCSfTat7LfoeUdLp2hblz/Y4iOXiZTJYAmSLSUUTq4xLGtPIHiUhT4Czg1bBtjUSkSdnXwDBglYexBlfZZEVLJjG3YAHUr++mm+wacCFtX6q4XpfxTteurrKN8Z5nyURVi4EbgZnAGuAFVf1YRK4XkevDDr0UmKWqB8K2tQbmi8hHwAfAG6r6plexBtratW6cpN21xZQqPP88nHWWm9pTmDOQ9C830mj9R36HllR69HAjuqpYEcBEST0vn1xVpwPTy217otz344Hx5bZtAHp4GVvSsPXeffHhh1BYCN27hzakprKr7/m0nfoon978tK+xJZM2bdyyyOvW2WBGr9kM+KCbN886333w73/D0KHHlk7Z3e8CTpjzIvX27/UtrmQj4hK69Zt4z5JJ0M2fb3cmMfbJJ7BpE/Tqdez24uOa83Xn3rSe9W9/AktSXbvCu+/6HUXwWTIJsp074auvoGNHvyNJKhMmwJAhUK+CRuRd/S/kpFces0b8GOrRw92Z2K/cW5ZMgqyshEpqqt+RJI3162HNGleWqyL7T+tBSnERTVdaffRYadfOrXC5YYPfkQSbJZMgmzvXhgTH2IQJMHiwGxJcIRF29Tmfk14dE8uwkpqIa3K0pi5vWTIJsjlzoFs3v6NIGhs2wIoV0L9/1cft6T2MFoveoF7hrtgEZujWDd5+2+8ogs2SSVAdPAirV1vZ+RgaN87dlVRXAq2k0XEU5gygzcwJMYnLuDuT2bOt38RLlkyCavFiyMy0BTRiZP16WLkSBg6M7PjdfYZz4utj7d0tRtq0cUl+zRq/IwkuSyZBNWeOrV8SQ08/7eaVVNpXUs6B03qQWnSY41Yv8jYw843cXHjnHb+jCC5LJkE1e7b1l8TI6tXw6acwYEANThJhV+/hnPjak57FZY7VsyfMnOl3FMFlySSIiopg6VJLJjEydiycey6kpdXsvD29h9Fq/lRSD37tTWDmGL16uQGOR4/6HUkwWTIJoiVL4OSToXFjvyMJvPx8+PJL6FPhAglVKz6uOftPy6XV7BeiH5j5jmbN3JyTxYv9jiSYLJkE0XvvWQmVGCgthccfh+HDaz8vdHfvYZz0+tjoBmYq1auXNXV5xZJJEL37bli5WuOV2bNdQulRh/rW+7L6kr5tAw2/WBu9wEylevWCGTP8jiKYLJkETVGRu4+vyzucqVZRkesrueACN8O61lLrsafX2bSZMS5qsZnKde3qlvjZvbv6Y03NWDIJmsWL4ZRToEkTvyMJtJdfhtat4fTT6/5ce3p/nzaznkNKiuv+ZKZK9eu7UV2zZvkdSfBYMgkaa+LyXGEhTJzo7kqi4fCJHSk+rgXHL7V6H7GQlwevveZ3FMHjaTIRkeEislZE1ovILRXsHywihSKyPPS4NdJzTSXeesvNzjKeefZZ14rYunX0nnPPGdbUFSt9+rhO+NJSvyMJFs+SiYikAo8B5wHZwCgRqaiE7TxVzQ097qjhuSbcgQNuvVibX+KZTZvcLOrvfz+6z7un51Caf/CmrcIYA23auGHC+fl+RxIsXt6Z9AHWq+oGVS0CJgMjYnBu8po71xV2zMjwO5LAGj0azj47+lN4Sho15etOeTbnJEZ694bXX/c7imDxMpm0BTaHfV8Q2lZefxH5SERmiEhZMalIz0VErhORfBHJ37FjRzTiTlwzZ1oTl4cWL4bNm2HQIG+ef88Z59BmxrPePLk5Rt++8OqrfkcRLF4mk4oGTJYvkboMOEVVewCPAq/U4Fy3UXWsquapal6rVq1qG2swzJzpehdN1BUXw6OPwkUXVbwcbzTsy+pDxpZ1NNyy3psLmG906wZffOEeJjq8TCYFQPuw79sBW8MPUNV9qro/9PV0IE1EWkZyriln61bYtg06dfI7kkB6+WU47jiPF65MrceenkNp/aatc+K11FS3iNm0aX5HEhxeJpMlQKaIdBSR+sBI4Jj/OhFpI+KmfIlIn1A8uyI515Tz5puuIdjWe4+6PXvgP/+BESPqOEExkmvlnUubWRNsqFEM9OvnPiSY6PAsmahqMXAjMBNYA7ygqh+LyPUicn3osMuBVSLyEfAIMFKdCs/1KtZAmDbNmrg88tRT7lcbzaHAlTnU9nRK0xrQdMU87y+W5Hr3djVR9+71O5Jg8HSeiapOV9VOqnqaqt4V2vaEqj4R+nq0quaoag9V7aeqC6s611SiqMgViqpN6VpTpbVrYeFCV2I+JkRcR/yb42N0weTVsKGr1WUTGKPDZsAHwYIF0L49HH+835EESmkp/OtfcN557o0nVvb0OptW818m5fDB2F00SQ0cCJMn+x1FMFgyCYJXX7W7Eg/MmgVHjrjmkFgqbtqSA6dk03L+K7G9cBIaMMBNzyos9DuSxGfJJNGpwtSpNVwz1lRn/35XFfiSSyDFh7+SPWecw4lWXsVzjRu7qVnW1FV3lkwS3cqVUFICp53mdySBMn48dO7sCjD7obDrIJqszaf+ji3+BJBEBg2CSZP8jiLxWTJJdGV3JV6PWU0iGze6+Z/nn+9fDFq/AXu7n0nrt/7tXxBJYtAg19Rla5zUjSWTRPfCC64X0USFqut0HzbM/yVh9uSdy4kznnVBGc80auTKq7z4ot+RJDZLJolszRrYudOqBEfR7Nmwa1d8dEEd6NiVlKJDNPlkid+hBN6QITDBCg/UiSWTRDZ5Mpx1lj89xAF08CCMGQOXXhonhQRE2H3GubbOSQz06eM+m23c6HckicvehRKVqus1HDzY70gCY8IEOPXU+BrLsKf3ME6YPYWUosN+hxJoaWnu7uS55/yOJHFZMklU+flw6BBkZfkdSSB8/jlMnw4XXuh3JMc6enxrDrXrRIv5Vi/da9//PjzzjJVFqy1LJolq3DjXS2yjuOpMFR580L2ZHHec39F81+68cznpjaf8DiPwOnWC9HSYM8fvSBKTJZNEdOQITJkSw4JRwTZzppuk2L+/35FUrLDb92j86VIafLnJ71ACTcR9Phs71u9IEpMlk0T08stw+uluMWtTJ4WF8Pjj8IMfxEmnewW0fgP29hpqM+Jj4Jxz4I033CBJUzOWTBLRo4/CBRf4HUUgjB4NPXvCySf7HUnVdvc5jxOnP+OqHRjPNG3qJjGOs7xdY5ZMEs2KFbBhg3cLkSeR/HxYtsxVBY53h9qeTnHjZjRfMtPvUALvoovgsccsb9eUJZNE889/uiFH8domkyAOHoT77oPLL4cGDfyOJjK7+p7HSa+O8TuMwOvSxc2Kf/NNvyNJLJZMEsmWLfDKK3DxxX5HkvCeeMLNKUmkkdV7ew6l6aoFNNhR4HcogSbilmh+4AG/I0ksniYTERkuImtFZL2I3FLB/itFZEXosVBEeoTt2ygiK0VkuYjkexlnwnjwQTfcJB7HryaQ/HyYPz/xcnJpg4auI37aE36HEnhDh7oZ8cuX+x1J4vAsmYhIKvAYcB6QDYwSkexyh30OnKWq3YE7gfKD8oaoaq6q2uLm27bBs8/CD3/odyQJ7euv4Z574IorICPD72hqbueAiznp9aeQo0V+hxJoaWluLZv77vM7ksTh5Z1JH2C9qm5Q1SJgMjAi/ABVXaiqe0LfLgLaeRhPYrvjDjerrlUrvyNJWKrwj39A165urZJEdKT1KRxucwqt3rMSt1678EKYMcNVRzDV8zKZtAU2h31fENpWmWuBGWHfKzBLRJaKyHWVnSQi14lIvojk79ixo04Bx601a9wkxZEj/Y4kob30kut2SvRR1TsHjqD9Cw9aaXqPNW7sRnbdfbffkSQGL5NJRXU+Knz1i8gQXDL5S9jmgaraC9dM9hsRObOic1V1rKrmqWpeqyB+aleFG26AK6+EZs38jiZhrVjhivj95CeuCSOR7cvuT9q+XRy3aqHfoQTeZZe5dU6++MLvSOKfl8mkAGgf9n07YGv5g0SkO/A0MEJVd5VtV9WtoX+3A1NxzWbJ5/nn3cfpSy7xO5KE9dVXcNttMGpUQFoJU1LYMegSTp5sDfpea9rUrbh5551+RxL/vEwmS4BMEekoIvWBkcC08ANE5GTgZeAqVf00bHsjEWlS9jUwDFjlYazxacsW+N3v4OabbV5JLRUWwp/+5JZ9SaRhwNXZ02c4TVfOJ+OLT/wOJfCuuMI1kX76afXHJjPPkomqFgM3AjOBNcALqvqxiFwvIteHDrsVaAGMKTcEuDUwX0Q+Aj4A3lDV5JpCdPSoG7l16aWQmel3NAlp/374859dNdizzvI7mugqbdDQ9Z1MvMfvUAKvaVP3p3jLdyY3mHCiAerEy8vL0/z8AExJUYUbb4SPPnKjuGwlxRorLIQ//hFOOsnl4yBW6k89UEiXe37K0qc+5HCbDn6HE2iHD8PPfgb//S8MHOh3NNElIkujMf3C3qXi0X33waxZ7qOQJZIa27ABrr/ezXAPaiIBKGnUlF39L+KU56xB32vp6fDLX8Kvf201uypj71Tx5sEHXVXgu+92YxNNxEpKXNv2734HZ5/tOk6DmkjK7DjrclrOe5n0LZ/5HUrgDRniRgI+/rjfkcQnSybxorTUNfA/8ogr5hiIYUexUVoKCxfCr37livP95jeQlyQ1E0oaHceOMy/j1KesQd9rIq71+W9/gwIrj/Yd9fwOwAC7d8PVV7vRW4884nr8ksihQ2747tdfu0UkwX0CbNjQ3Zw1auRKn9Sr57qTjhyBHTtg40bXrTRnjjtm6FDo1i35WgZ3nnkZXe75KU1WL+br7L5+hxNoHTu6Ufq/+IWbHR/0O9+asA54v737LlxzjVsz9pe/TPwZdRHYuxfefx8WLXKT+/fuhebNXeIo+/GPHoWiIpdoyh4lpW5f/TQ4/ng44QS3qFXXrnDiiX79NPHh+A/epHn+LJY+vsSGkXusuNjd/f7pT+5PNtFFqwPe7kz8smePmz/yxhvw+99Dn2DPyVSFJUvcisMrVrg1Izp3hgEDXIteJHcTJSXuk2Cy3XlEYk/eMJp/8CYnvvEU2y6+vvoTTK3Vqwf/+7/whz+4kV3Z5cvXJilLJrFWWuqq//71r261xKeeCnRHe2kpzJvnfuTiYvfHd8klbnRMTdkH7iqkpLDlB//DaU/+md39LuDICe2rP8fU2imnuKauyy6DDz6AJk38jsh/1swVSwsXuh68khJXbytRS9dGaPlyt8b60aNuGZbsbGtj9toJb/2H9C83suL+WXYL5zFVeOgh94Hp1VcT99dt80wSSUGBq/j7gx+4BccfeijQieTLL+H//g/uuss1Y/32t5CTY4kkFrYPHUn9vTs4edK9focSeGWju7ZscS3VAfpcXiuWTLx0+LCrENetm1to/Nln4dxzA/uuWlTkKvP+8pduMcibb4ZevRL3E1tCSq3Hpqv+j3Yv/ovmHyRXBSI/pKXB7bfD9OlWDNL+zL2gCq+95ioLvvMOjBnjajE0bOh3ZJ5ZtMiNbl62zH1KGzYM6tf3O6rkdLRZKzZecytZd/2ExmuX+h1O4DVpAvfeC+PGuTkoyXqHYh3w0bZuHdx0E6xd6/pFevf2OyJPFRS4CfubNsGIEcGqzJvIDnbsSsHlv6P7X4az8p7pfN0l2K9DvzVv7uYa//WvsHMnPPywG/WVTOzOJFr273e1tPr2hQ4dYOzYQCeSwkLXuf7rX0Pr1q6ooiWS+FLYbRAFl/+e7n8eTqvZL/gdTuA1b+6qIS1b5u7Mg7rwa2UsmdRVaanrKOjUCT78EJ580nW2B3Ty4YEDMGECXHWV62i/+WY38zygP27C29d1ABuuu4fTHv8jXe66irTdX/kdUqA1bgx//7urVt29u5slnyxsaHBtqcJbb7l309JSV6Y2Jyc21/bBnj1uwuErr7g7kHPPtfJhiSTlyCHazJzA8Utm8tWwq9h23rUcOK17YAeDxINly1zT14AB7o6lfZxO/YnW0GBLJjVVlkRuv92NCfzpT93KSwH8o1SFlSth2jRX/qRnT/ejWhJJXGl7ttPi/ddptnw2KcVH+bpTLw61zaTo+BPQtAZQWkLq4YPU27+HtH27qbdvF/X276Xewa9JOXwQKTkKIpSmNaAkowlHm7bi8AntOdQuk0PtO7P/1O4cPrFjIP8eauPwYZg40Y3HueYaV4KlXTu/ozqWJZMKeJpMDhyAyZPdHJFDh9xankOHBm5atqpbnnTuXDcQLSXFdf307h3oifrJR5X6u7bRcMs66u/+itSDX5NSXISmpKBpDShJz6Ak4ziKGzamJKMJpemNKE1rgNZLg9JSpLjo26RTuJP6u7aRvn0zDbd+RkrRYb7u1IvCnIEUdv8e+7L7UdLoOL9/Yl/t3AkvvOCWKTrnHFfhesiQ+OikT4hkIiLDgYeBVOBpVb2n3H4J7T8fOAj8VFWXRXJuRaKeTA4dcu+okye7jxbdu8OFF7p31oBMnigpgS++gNWrXZfPsmVuSG9ODuTmuk9R9iHT1ES9r3eTsekTMjatptGmNWR88QmH2p5GYbfvUdh1IPuy+yXt3cv+/a5h4+23Yft2uOgiuOACGDwYWrTwJ6a4TyYikgp8CpwLFABLgFGqujrsmPOBm3DJpC/wsKr2jeTcitQpmai6ca4rVrhiO3PmQH7+t9UIa/G/XVrqyqrv3ev6HML/LSyEffvc/kOH3O1wUZF7c1d1f2epqe6TS3q6m6LSqJF7NGniJgU2afJtifZGjdwx9eu7h4h7npIS97wHDrhr7d7tyr1v2eKG8xYUuOfq0MHVG+rUyZqxTHRJ8VEaFnxKo89XkfHFJzTa+DFSUsL+03vw9ek9OXBqdw6e3IVD7TIpbnJ80iSZLVtc8/GyZe5tp21b9zm1Z0/3YS4z01XF9npwSyIkk/7Abar6/dD3fwVQ1X+EHfMk8J6qTgp9vxYYDHSo7tyKVJhMZs2CBQvcIhgHD7qPBvv2uXf0HTtg27ZjxvAdSW/KnlaZ7Gqdw1dtunOkXiNKSvjmUVzs3pzLSqQfOXJsmfSDB+FI0XdjS6vn1uTIyPg2OaSnu4nx9eu7F0y9eu7vqOxvqbT02GsWFbmkc+TIsY9vElFp5b+bFHHXzMhwCahpUzeUsUULm1xoYi99/06abl/PcTs30GT3Jprs3nTM/iMNm3KoyQkcbtySIxnNKWp4HMUNGlNcP4PitHRK6zWgJDWN0tQ0vjz9e3x12gCffpLoKCmBzz93Tcyff+4+6O3b9+3+lBS35EKrVm75haZN3YfJsrV+GjRwj7Q090hNdY+UFPfo0sUNV65IIpSgbwtsDvu+AHf3Ud0xbSM8FwARuQ64LvTtERFZFb4/G7IaQkZ1wRYhR4tJT+NwEWz+2D2o/dh8QRXCEnUxsA/27CvheHzoZ1FcQ+JBYGf1h+/BpzhryOKMLv/ibABAKqWkcVTkUCEcKoTt6yo8egdQdgO9j0al6+hwJDZx1tTuetC8uPbni0C9tNLS1NQvv3TD8WtHFT5cdsx70reiUijQy2RS0b1q+Z+ksmMiOddtVB0LjAUQkfxoZFgviUj+Vj0a1zGCxRltFmd0iUj+pjj/W4ey96StcR2niESlo9nLZFIAhI+sbgdsjfCY+hGca4wxJk54OSRpCZApIh1FpD4wEphW7phpwNXi9AMKVXVbhOcaY4yJE57dmahqsYjcCMzEDe8dp6ofi8j1of1PANNxI7nW41rzf1bVuRFcdmz0f5KoS4QYweKMNoszuizO6IlKjIGatGiMMcYfwZh5Z4wxxleWTIwxxtRZQiQTERkuImtFZL2I3FLBfhGRR0L7V4hIr0jPjXGcV4biWyEiC0WkR9i+jSKyUkSWR2uoXh3iHCwihaFYlovIrZGeG+M4bw6LcZWIlIhI89C+mPw+RWSciGwvP78pbH+8vDarizNeXpvVxen7azOCGH1/XYau1V5EZovIGhH5WER+W8Ex0Xt9qmpcP3Ad8J8Bp+KGDH8EZJc75nxgBm5+Sj9gcaTnxjjOAcDxoa/PK4sz9P1GoGWc/D4HA6/X5txYxlnu+IuAd334fZ4J9AJWVbLf99dmhHH6/tqMMM54eG1WGWM8vC5D1zoR6BX6ugmuRJVn752JcGfSB1ivqhtUtQiYDIwod8wI4Dl1FgHNROTECM+NWZyqulBV94S+XYSbPxNrdfmdxNXvs5xRwCSPYqmUqs4FdldxSDy8NquNM05em5H8PisTs99nDWP05XUJoKrbNFQ4V1W/BtbgqouEi9rrMxGSSWUlVyI5JpJzo6Wm17oW94mgjAKzRGSpuBIxXok0zv4i8pGIzBCRslW/4vL3KSIZwHDgpbDNsfp9ViceXps15ddrM1J+vzYjEk+vSxHpAPQEFpfbFbXXZxxU069WTMqyREHE1xKRIbg/2EFhmweq6lYROQF4S0Q+CX0C8iPOZcApqrpfXGXnV4DMCM+Nlppc6yJggaqGf1qM1e+zOvHw2oyYz6/NSMTDazNScfG6FJHGuIT2O1XdV353BafU6vWZCHcmdSnLEsm50RLRtUSkO/A0MEJVd5VtV9WtoX+3A1Nxt5m+xKmq+1R1f+jr6UCaiLSM5NxYxhlmJOWaEmL4+6xOPLw2IxIHr81qxclrM1K+vy5FJA2XSJ5X1ZcrOCR6r89YdATV5YG7e9oAdOTbjqCccsdcwLGdSB9Eem6M4zwZN9t/QLntjYAmYV8vBIb7GGcbvp3Q2gf4IvS7javfZ+i4prj260Z+/D5D1+hA5R3Gvr82I4zT99dmhHH6/tqsLsY4el0K8BzwUBXHRO31GffNXOpPWRav4rwVaAGMEbdoSbG6yqetgamhbfWAiar6po9xXg78WkSKgUPASHWvsHj7fQJcCsxS1QNhp8fs9ykik3AjjFqKSAHwNyAtLEbfX5sRxun7azPCOH1/bUYQI/j8ugwZCFwFrBSR5aFt/4v74BD116eVUzHGGFNnidBnYowxJs5ZMjHGGFNnlkyMMcbUmSUTY4wxdWbJxBhjTJ1ZMjHGGFNnlkxMUhKRDiJyKGz8PSLSTkReFZF1IvKZiDwsIvUrOXdV6OvBIvJ66OsrQuW6Xw879kQRmSXflk6fHrYvR0TeFZFPQ9f8/yQ0CaHc9cKv8VMRGR36+jYR+VPo6/tF5Muy742JNUsmJpl9pqq54NZ1AF4GXlHVTKAT0Bi4K9InU9UpwC/KbR6Om/gFME9Vzw9dryEwDbhHVTsBPXBl4G+ozQ+iqjcDT1R7oDEesWRijDMUOKyqzwKoagnwe+DnoeqvtTWcYyvwlvkxrgjgrND1DgI3Ap4ukmWMVyyZGOPkAEvDN6irsPoFcHptnlBEUoHOqro6wut9BjQWkeNqcz1j/GTJxBhHqLjEdmXbI9GX764fEcnzWo0jk3AsmRjjfAzkhW8I3SG0xy1fWhvnAZUV8qvoeqcC+9WtimdMQrFkYozzDpAhIlfDN01UDwLjQ/0ZtXF26Hkr8jwwSETOCV2vIfAIcF8tr2WMryyZGAOEyphfCvxQRNYBnwKHcSW7a0xEWuE69MuvbFd2vUO4NbX/T0TWAiuBJcDo2lzPGL/F/XomxsSKqm7GLbVa3XEbga6hr98D3qvgsO8Ds6p5npW4dTGqu94311DV8cD40Ne3VXeuMbFidyYmWZUATcMnLdaViFwBjAH2qOp/VPWesN1FQNfwSYvRJCL3Az8BDlR3rDFesMWxjDHG1JndmRhjjKkzSybGGGPqzJKJMcaYOrNkYowxps7+f3i9UWYX35LOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.kdeplot(oii_neiii_sf,color='Blue',fill=True,label='SF')\n",
    "sns.kdeplot(oii_neiii_agn,color='Red',fill=True,label='AGN')\n",
    "plt.xlim(0.0,2.0)\n",
    "plt.xlabel('[O II]/[O III]')\n",
    "plt.ylabel('Density')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The KDE plot of [O II]/[Ne III] shows that no matter how we slice it, this ratio alone is not capable of separating SF and AGN excitation in galaxies. Let's assess the ability of [O III] $\\lambda$5007 /[O I] $\\lambda$6300 to act as an effective excitation diagnostic. Try generating a KDE plot for this emission line ratio in the space below and ask yourself the following questions.\n",
    "\n",
    "- What are relavent ionization potentials for [O III] and [O I] line emission? \n",
    "- Would you expect AGN to display higher or lower [O III]/[O I] ratios compared to SF galaxies?\n",
    "- Can you make a single cut to separate all galaxies with an AGN from those only experiencing SF?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Your answers to these questions:</b>\n",
    "<br>\n",
    "<br>\n",
    "<br>\n",
    "<br>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code for KDE plot\n",
    "#oiii_oi_sf = \n",
    "#oiii_oi_agn = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Diagnostic Diagrams\n",
    "\n",
    "We can combine multiple line ratios in a diagram to increase their utility in separating AGN and SF activity. This assesses the shape of the inicident radiation field (starlight or AGN) at more energies than possible with a single emission line ratio. Let's combine [O II]/[Ne III] and [O III]/[O I] into a single plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'oiii_oi_sf' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m\u001b[0m",
      "\u001b[0;31mNameError\u001b[0mTraceback (most recent call last)",
      "\u001b[0;32m<ipython-input-18-cf6c5fbac493>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mscatter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moii_neiii_sf\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0moiii_oi_sf\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcolor\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'blue'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mlabel\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'SF'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.5\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0malpha\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.75\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mscatter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moii_neiii_agn\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0moiii_oi_agn\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcolor\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'red'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mlabel\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'AGN'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.5\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0malpha\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.75\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mxlabel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'[O II]/[Ne III]'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mylabel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mr'[O III]/[O I]'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'oiii_oi_sf' is not defined"
     ]
    }
   ],
   "source": [
    "plt.scatter(oii_neiii_sf,oiii_oi_sf,color='blue',label='SF',s=0.5,alpha=0.75)\n",
    "plt.scatter(oii_neiii_agn,oiii_oi_agn,color='red',label='AGN',s=0.5,alpha=0.75)\n",
    "\n",
    "plt.xlabel('[O II]/[Ne III]')\n",
    "plt.ylabel(r'[O III]/[O I]')\n",
    "\n",
    "plt.xlim(0.2,1.8)\n",
    "plt.ylim(0.3,2.0)\n",
    "\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adding AGN/SF demarcations\n",
    "\n",
    "The [O III]/[O I] vs. [O II]/[Ne III] diagram does a reasonable job in separating AGN and SF galaxies. We can make this separation more explicit by adding a demarcation line. Since we might want to generate this plot over and over, we can make define a function that will do it for us."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_diagram():\n",
    "\n",
    "    o2ne3_demarc = np.linspace(0.2, 1.19)\n",
    "    o3o1_demarc = (0.055 / (o2ne3_demarc - 1.19)) + 1.35\n",
    "\n",
    "    plt.scatter(oii_neiii_sf,oiii_oi_sf,color='blue',label='SF',s=0.5,alpha=0.75)\n",
    "    plt.scatter(oii_neiii_agn,oiii_oi_agn,color='red',label='AGN',s=0.5,alpha=0.75)\n",
    "\n",
    "    plt.plot(o2ne3_demarc,o3o1_demarc,'k')\n",
    "\n",
    "    plt.xlabel('[O II]/[Ne III]')\n",
    "    plt.ylabel(r'[O III]/[O I]')\n",
    "\n",
    "    plt.xlim(0.2,1.8)\n",
    "    plt.ylim(0.3,2.0)\n",
    "\n",
    "    plt.legend()\n",
    "\n",
    "######    \n",
    "\n",
    "make_diagram()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Justifying diagnostic diagrams with models\n",
    "\n",
    "Up until this point, we have been using observed galaxies flagged as an AGN or star forming galaxy in the <code>GalSpecExtra</code>. While we would like to take this information for granted, it is important to note these tags are based on a set of assumptions that may or may not hold for all galaxies. We can use photoionization models to understand the limits of these AGN/SF distinctions and our newly created diagnostics diagram. \n",
    "\n",
    "Photoionization models predict the observed emission line spectrum for gas clouds illuminated by an AGN or starlight in a galaxy. The set that we will use comes from [this paper](https://iopscience.iop.org/article/10.3847/1538-4357/ac510c) and a link to the simulated emission region data has been provided on Moodle. To make the visualization easier, we can restrict the model parameter space to regions that we're interested in exploring. A brief description of the parameters is given below along with accepted values within the grid.\n",
    "\n",
    "<code>LOGU</code> - dimensionaless quantity (ionization parameter) specifying the ionization level of the gas clouds\\\n",
    "options: -4.0 to -0.5 in 0.25 increments\n",
    "\n",
    "<code>LOGZ</code> - metallicity of the gas clouds relative to solar metallicity\\\n",
    "options: [0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.7, 1.0, 1.5, 2.0]\n",
    "\n",
    "<code>AGNFRAC</code> - the fraction of AGN excitation as opposed to starlight excitation\\\n",
    "options: [0.0, 0.04, 0.08, 0.16, 0.32, 0.64, 1.0] \n",
    "\n",
    "Let's see how the emission line predictions compare to our observed galaxies on the diagnostic diagram when we select a simulation with 100% AGN excitation and Z = 0.2 Z$_{\\odot}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget --no-check-certificate 'https://drive.google.com/uc?export=download&id=1vQrGvdOf2w9VAu7F0lub5tEAi75cIMm4' -O model_file\n",
    "\n",
    "model_grid = pd.read_csv('model_file', delimiter=',')\n",
    "\n",
    "#print (model_grid.columns) # useful to check out the column names in the file\n",
    "\n",
    "# function to restrict the model parameter space\n",
    "def apply_limits(grid,f_agn,U,Z):\n",
    "    \n",
    "    grid = grid[(grid['AGNFRAC'] == f_agn) & (round(grid['LOGZ'],5) == round(np.log10(Z),5)) \\\n",
    "                & (grid['MIXING'] == ' coincident') & (grid['LOGU'] >= U[0]) & (grid['LOGU'] <= U[1])]\n",
    "    \n",
    "    return grid\n",
    "\n",
    "U = [-3.75,-2.75] # range of ionization parameter\n",
    "f_agn = 1.0\n",
    "Z = 0.2\n",
    "\n",
    "agn_models = apply_limits(model_grid,f_agn,U,Z)\n",
    "#print (agn_models)\n",
    "\n",
    "model_x = np.log10((agn_models[' \\'O__2_372603A\\'']+agn_models[' \\'O__2_372881A\\''])/agn_models[' \\'NE_3_386876A\\''])\n",
    "model_y = np.log10(agn_models[' \\'O__3_500684A\\'']/agn_models[' \\'O__1_630030A\\''])\n",
    "\n",
    "make_diagram()\n",
    "\n",
    "plt.plot(model_x,model_y,color='Green',marker='o')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The AGN model assuming 0.2 Z$_{\\odot}$ lies in the AGN part of our diagram and including a range of ionization parameter covers some of the variation in the region. It seems that the models preliminarily suggest the galaxies flagged as AGN really are AGN and the diagnostic diagram can pick them out. In the cells below, determine if this is a robust result by answering the following questions.\n",
    "\n",
    "- Do all of the 100% AGN simulations fall within the AGN part of the diagnostic diagram?\n",
    "- If so, what range of ionization parameters and metallicities do the observed AGN span? If not, what values disagree with our chosen demarcations?\n",
    "- Do all of the purely SF simulations fall within the SF part of the diagnostic diagram?\n",
    "- If so, what range of ionization parameters and metallicities do the observed SF galaxies span? If not, what values disagree with our chosen demarcations?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code to answer the questions above"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The BPT Diagram\n",
    "\n",
    "The diagnostic diagram that we have developed does have it's limitations. For example, both of the [O II] and [Ne III] emission lines are unavailable for low-$z$ galaxies because they fall outside the SDSS observing window. This means that we're unable to identify AGN activity for local galaxies. Another limitation is that some of the emission lines can be rather weak, making a reliable detection difficult. We can quantify this by looking at much our sample reduces down after applying our S/N cut."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print ('Fraction of initial SF galaxies included after S/N cut: ',\n",
    "       \"{:.3}\".format(len(gals_sf)/len(gals[gals['bptclass']==1])))\n",
    "print ('Fraction of initial AGN galaxies included after S/N cut: ',\n",
    "       \"{:.3}\".format(len(gals_agn)/len(gals[gals['bptclass']==4])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therein lies the problem: most galaxies only have a handful of reliable emission lines measurements. This is one reason why the diagnostic diagram we have been analyzing isn't a popular choice in research literature even though it does a reasonable job of distinguishing between AGN and star forming galaxies.\n",
    "\n",
    "The so-called \"gold standard\" for classifying AGN with optical spectroscopy is known as the \"BPT\" diagram, named after the authors that developed it in [this paper](https://ui.adsabs.harvard.edu/abs/1981PASP...93....5B/abstract). The BPT diagram uses emission line ratios from [O III] $\\lambda$5007, [N II] $\\lambda$6584, H$\\beta$, and H$\\alpha$, all of which are known as \"strong lines\" since they're strong enough to be reliably measured in most galaxies. Using strong lines for diagnostics drastically improves the number of galaxies we can include in a sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gals_flux = gals[(gals['oiii_5007_flux']>0) & (gals['nii_6584_flux']>0)]\n",
    "\n",
    "quality_gals = gals_flux[(gals_flux['oiii_5007_flux']>5.0*gals_flux['oiii_5007_flux_err']) &\n",
    "                (gals_flux['nii_6584_flux']>5.0*gals_flux['nii_6584_flux_err'])]\n",
    "\n",
    "gals_sf = quality_gals[quality_gals['bptclass']==1 | (quality_gals['bptclass']==2)]\n",
    "gals_agn = quality_gals[(quality_gals['bptclass']==4)]\n",
    "\n",
    "print ('Fraction of initial SF galaxies included after S/N cut: ',\n",
    "       \"{:.3}\".format(len(gals_sf)/len(gals[gals['bptclass']==1])))\n",
    "print ('Fraction of initial AGN galaxies included after S/N cut: ',\n",
    "       \"{:.3}\".format(len(gals_agn)/len(gals[gals['bptclass']==4])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The BPT diagram remains a hot topic in academic literature today. In fact, the <code>GalSpecExtra</code> schema actually uses the BPT diagram to classify galaxies as AGN or star forming. If we use the standard demarcations from literature, it is no surprise in the plot the below that all of the AGN and star forming galaxies are perfectly separated. There is even an additional class of galaxies known as \"composites\" that are thought to represent a transition state between purly star forming galaxies and those excited purely by an AGN."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_bpt_plot():\n",
    "    \n",
    "    gals_comp = quality_gals[(quality_gals['bptclass']==3)]\n",
    "\n",
    "    oiii_Hb_sf = np.log10(gals_sf['oiii_5007_flux']/gals_sf['h_beta_flux'])\n",
    "    oiii_Hb_agn = np.log10(gals_agn['oiii_5007_flux']/gals_agn['h_beta_flux'])\n",
    "    oiii_Hb_comp = np.log10(gals_comp['oiii_5007_flux']/gals_comp['h_beta_flux'])\n",
    "\n",
    "    nii_Ha_sf = np.log10(gals_sf['nii_6584_flux']/gals_sf['h_alpha_flux'])\n",
    "    nii_Ha_agn = np.log10(gals_agn['nii_6584_flux']/gals_agn['h_alpha_flux'])\n",
    "    nii_Ha_comp = np.log10(gals_comp['nii_6584_flux']/gals_comp['h_alpha_flux'])\n",
    "\n",
    "    plt.scatter(nii_Ha_sf,oiii_Hb_sf,color='blue',label='SF',s=0.5,alpha=0.75)\n",
    "    plt.scatter(nii_Ha_agn,oiii_Hb_agn,color='red',label='AGN',s=0.5,alpha=0.75)\n",
    "    plt.scatter(nii_Ha_comp,oiii_Hb_comp,color='magenta',label='Comp',s=0.5,alpha=0.75)\n",
    "\n",
    "    comp_n2ha = np.linspace(-2.5,-0.1,100)\n",
    "    comp_o3hb = 0.61/(comp_n2ha-0.05)+1.3\n",
    "    plt.plot(comp_n2ha,comp_o3hb,color='k',linestyle='--',linewidth=1.0)\n",
    "\n",
    "    agn_n2ha = np.linspace(-2.5,0.3,100)\n",
    "    agn_o2hb = 0.61/(agn_n2ha-0.47)+1.19\n",
    "    plt.plot(agn_n2ha,agn_o2hb,color='k',linewidth=1.0)\n",
    "\n",
    "    plt.xlabel(r'[N II]/H$\\alpha$')\n",
    "    plt.ylabel(r'[O III]/H$\\beta$')\n",
    "\n",
    "    plt.xlim(-1.5,0.5)\n",
    "    plt.ylim(-1.0,1.25)\n",
    "\n",
    "    plt.legend()\n",
    "\n",
    "######\n",
    "    \n",
    "make_bpt_plot()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All of the strong lines on the BPT diagram fall within the SDSS observation window for both local galaxies and those at higher redshifts. Combined with the superior statistics shown above, this makes the BPT diagram a suitable choice for tracing galaxy evolution. We can return to the photoionization models to investigate how a purely star forming model compares to the selected observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "make_bpt_plot()\n",
    "\n",
    "U = [-3.75,-2.75]\n",
    "f_agn = 0.0\n",
    "Z = 0.4\n",
    "\n",
    "agn_models = apply_limits(model_grid,f_agn,U,Z)\n",
    "\n",
    "#print (model_grid.columns) # useful to check out the column names in the file\n",
    "\n",
    "model_x = np.log10(agn_models[' \\'N__2_658345A\\'']/agn_models[' \\'H__1_656281A\\''])\n",
    "model_y = np.log10(agn_models[' \\'O__3_500684A\\'']/agn_models[' \\'H__1_486133A\\''])\n",
    "\n",
    "plt.plot(model_x,model_y,color='Cyan',marker='o')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All diagnostic diagrams have their limitations and it is important to keep the limitations in mind when interpreting results. In the space below, use the photionization models to explore the limitations of the BPT diagram in the same way that you previously explored the limitations of the [O III]/[O I] vs. [O II]/[Ne III] diagram. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code to compare photoionization models to selected observations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Follow Up Questions\n",
    "\n",
    "Developing new diagnostic diagrams to separate AGN and SF galaxies remains a trendy topic in academic research. This Lecture-Tutorial has laid the foundation for you to investigate this yourself. A sample of potential follow up questions is given below, but don't let these limit your curiousity.\n",
    "\n",
    "- What are the effects of selecting models with non-binary $f_{AGN}$?\n",
    "- We know where BPT classfied AGN lie on the [O III]/[O I] vs. [O II]/[Ne III] diagram, but what about vice-versa, and what does this tell you about using only one diagram to find AGN?\n",
    "- What other diagrams can effectively seperate AGN and SF galaxies?"
   ]
  },
  {
   "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.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}