HW2

{
 "metadata": {
  "name": "",
  "signature": "sha256:101a92bab1e51ebbcfae72e79a03dc466e53ce36a04ab6fab03dd96b2e2f2770"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## A535 (ISM) Homework #2\n",
      "\n",
      "Tom Rice\n",
      "\n",
      "9 February 2015\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### 1\\. \n",
      "*An 06V star with $T_{eff}$ = 43,600 K ionizes an H II region. Assume an ionization-bounded nebula of pure hydrogen, with uniform density $n_H = 10cm^{\u22123}$ and temperature T=10,000 K. You can find the necessary data in the slides for Jan. 14 and Jan. 23.*\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### (a)\n",
      "*What is the (Str\u00f6mgren) radius of the nebula?*"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "A Str\u00f6mgren radius is given by the expression\n",
      "$$ R_s = \\left( \\frac{3 Q}{4 \\pi n_H^2 \\alpha_B} \\right)^{1/3}$$ \n",
      "\n",
      "where \n",
      "\n",
      "* $Q = 10^{49.34}$ photons/sec is the stellar radiation coefficient (given by table from 1/14 slides, attributed to Osterbrock & Ferland),\n",
      "* $n_H = 10$ cm$^{-3}$ as given above\n",
      "* $\\alpha_B = 2.59 \\times 10^{-13}$ cm$^3$ s${-1}$ is the hydrogen recombination coefficient for a 10,000 K gas (also from slides / Osterbrock & Ferland)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from __future__ import division\n",
      "import numpy as np\n",
      "import astropy.units as u\n",
      "import astropy.constants as c\n",
      "\n",
      "Q = 10**49.34 / u.s\n",
      "n_H = 10 * u.cm**(-3)\n",
      "alpha_B = 2.59e-13 * u.cm**3 / u.s\n",
      "\n",
      "R_s = (3 * Q / (4 * np.pi * n_H**2 * alpha_B))**(1/3)\n",
      "\n",
      "print \"The Stromgren radius is {0:.2f}\".format(R_s.to(u.pc))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The Stromgren radius is 19.00 pc\n"
       ]
      }
     ],
     "prompt_number": 234
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### (b)\n",
      "*Compute the typical recombination timescale in this nebula.*"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The recombination time is given by\n",
      "$$ t_r = \\frac{1}{n_e \\alpha} $$\n",
      "with \n",
      "\n",
      "* $\\alpha \\approx 2.59 \\times 10^{-13}$ as in (a), and \n",
      "* $n_e$ = 10 cm$^{-3}$ assuming complete H ionization\n",
      "\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n_e = n_H\n",
      "t_r = 1/(n_e * alpha_B)\n",
      "\n",
      "print \"The recombination time is {0:.2e}, or {1:.1f}\".format(t_r.to(u.s), t_r.to(u.yr))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The recombination time is 3.86e+11 s, or 12234.8 yr\n"
       ]
      }
     ],
     "prompt_number": 235
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### (c)\n",
      "*<b>Estimate</b> the timescale for a neutral hydrogen atom to become photoionized at a distance of half the Str\u00f6mgren radius of the nebula.*"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The photo-ionization volumetric rate is given as\n",
      "\n",
      "$$ volumetric rate_{PI} = n_H \\int_{\\nu_0}^{\\infty} \\frac{L_\\nu}{4 \\pi r^2} \\frac{a_\\nu}{h \\nu} d \\nu $$\n",
      "\n",
      "which can be re-written as\n",
      "\n",
      "$$ n_H \\frac{Q}{4 \\pi r^2} a_\\nu $$\n",
      "\n",
      "with $Q = \\int_{\\nu_0}^{\\infty} \\frac{L_\\nu}{h \\nu} d \\nu$.\n",
      "\n",
      "Here, $a_\\nu$ represents the ionization cross-section of neutral hydrogen, given by $6 \\times 10^{-18}$ cm$^{-2}$.\n",
      "\n",
      "The rate per atom factors out the $n_H$ term, leaving\n",
      "\n",
      "$$ rate_{PI} = \\frac{Q}{4 \\pi r^2} a_\\nu$$\n",
      "\n",
      "and the timescale for photoionization is just the inverse of the rate,\n",
      "\n",
      "$$ t_{PI} = \\frac{1}{rate_{PI}} $$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a_nu = 6e-18 * u.cm**2\n",
      "rate_PI = Q / (4 * np.pi * (R_s/2)**2)* a_nu\n",
      "\n",
      "t_PI = 1/rate_PI\n",
      "\n",
      "print \"The photoionization timescale is {0:.1e}, or {1:.2f}\".format(t_PI.to(u.s), t_PI.to(u.year))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The photoionization timescale is 8.2e+07 s, or 2.61 yr\n"
       ]
      }
     ],
     "prompt_number": 236
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### (d)\n",
      "Use the results of (b) and (c) to estimate the ratio of neutral hydrogen to protons at the half-radius point."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If the ratio can be expressed as a fraction of the two timescales\n",
      "\n",
      "$$ \\frac{n_{H^0}}{n_p} = \\left( \\frac{t_r}{ t_{PI} } \\right)^{-1} $$\n",
      "\n",
      "then it is given by the following:\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "H0_proton_ratio = t_PI / t_r\n",
      "\n",
      "print \"The ratio between neutral hydrogens and free protons\" \n",
      "print \"is {0:.2e}\".format(H0_proton_ratio)\n",
      "\n",
      "print \"So there is about one neutral atom per {0:.0f} free protons.\".format(1/H0_proton_ratio)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The ratio between neutral hydrogens and free protons\n",
        "is 2.13e-04\n",
        "So there is about one neutral atom per 4691 free protons.\n"
       ]
      }
     ],
     "prompt_number": 237
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### (e)\n",
      "*Use the results of (d) to estimate the optical depth to the star from $R_s/2$ at the head of the Lyman continuum. Is this reasonably consistent with your calculation in (c)?*"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Lyman photons are absorbed by neutral hydrogens, not hydrogen ions.\n",
      "\n",
      "Optical depth $\\tau$ of this emission line is given by\n",
      "\n",
      "$$ \\tau_\\nu(r) = \\int_0^r n_{H^0} a_\\nu ds$$\n",
      "\n",
      "Assuming a constant neutral-hydrogen density through the relevant path length,\n",
      "\n",
      "$$ \\tau = n_{H^0} a_\\nu d $$\n",
      "\n",
      "where $d$ is the path length."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "tau = n_H * H0_proton_ratio * a_nu * R_s/2\n",
      "\n",
      "print (\"The optical depth for Lyman photons \\n\"\n",
      "\"from the star to half the Stromgren radius\")\n",
      "print \"is {0}\".format(tau)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The optical depth for Lyman photons \n",
        "from the star to half the Stromgren radius\n",
        "is 0.375\n"
       ]
      }
     ],
     "prompt_number": 238
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This seems consistent with the idea that about half the photons should be absorbed within half the Str\u00f6mgren sphere, within an order of magnitude!"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### (f)\n",
      "*Compute the H$\\alpha$ luminosity emanating from this nebula.*\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The H$\\alpha$ emission is given by\n",
      "$$ L(\\textrm{H} \\alpha) = \\int_0^{R_s} 4 \\pi j_{H \\alpha} dV$$\n",
      "\n",
      "which can be rewritten as \n",
      "$$ \\frac{4 \\pi}{3} R_s^3 (4 \\pi j_{H \\alpha} ) $$.\n",
      "\n",
      "The table \"Table 4.4\" in notes from 1/24 suggests that at nearly all relevant densities, the strength of the H$\\beta$ line can be expressed as \n",
      "\n",
      "$$ \\frac{4 \\pi j_{H \\beta} }{ n_e n_p }\\approx 1.2 \\times 10^{-25} \\textrm{erg cm}^3 \\textrm{s}^{-1}$$\n",
      "\n",
      "and at $n_e \\sim 10$ the ratio of H$\\alpha$ relative to H$\\beta$ can be written as\n",
      "\n",
      "$$ \\frac{j_{H \\alpha} }{ j_{H \\beta}} \\approx 2.9 $$.\n",
      "\n",
      "In this formulation, $$j_{H \\alpha} = \\frac{1.2 \\cdot 2.9 \\cdot n_e n_p}{4 \\pi }\\times 10^{-25} \\textrm{erg cm}^3 \\textrm{s}^{-1} $$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "j_Ha = 1.2 * 2.9 * n_H**2 / (4*np.pi) * 10**-25 * u.erg * u.cm**3 / u.s\n",
      "\n",
      "L_Ha = 4/3*np.pi * R_s**3 * (4*np.pi * j_Ha)\n",
      "\n",
      "print \"The Stromgren sphere's H-alpha luminosity\"\n",
      "print \"is {0:.2e}, or {1:.1f} solar luminosities\".format(L_Ha.to(u.erg/u.s), L_Ha.to(u.solLum).value)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The Stromgren sphere's H-alpha luminosity\n",
        "is 2.94e+37 erg / s, or 7643.1 solar luminosities\n"
       ]
      }
     ],
     "prompt_number": 239
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### 2\\. \n",
      "Compute the thermal equilibrium of a model H II region. Assume that only H, N, and O are present. Assume that the O abundance is $7 \\times 10^{-4}$ relative to H, and that 80% of the oxygen is O$^+$ and 20% is O$^{++}$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### (a)\n",
      "*Compute the net nebular heating rates per $n_e n_p$ for a 40,000 K star using the attached Table 3.1 for $\\tau_o = 0$ at each of the available gas temperatures with recombination cooling values in Table 3.2 (use Case B).*"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "I assume that the \"net nebular heating rate\" is equivalent to a \"net effective heating rate\" (G-L$_R$) as described by Pogge on pg. 22.\n",
      "\n",
      "For these parameters, $T_i$ (the mean input energy of photoelectrons) is 2.48 $\\times 10^4$ K. Values of $\\alpha_B$ come from Table 2.1 in the slides of 1/14.\n",
      "\n",
      "T       $\\alpha_B / 10^{-13}$   $\\beta_B / 10^{-13}$\n",
      "\n",
      "5000    4.54       3.24\n",
      "\n",
      "10000   2.59       1.73\n",
      "\n",
      "20000   1.43       0.875"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def energy_gain_per_n2(T_i, alpha):\n",
      "    G = 3/2 * c.k_B * T_i * alpha\n",
      "    return G.to(u.erg / u.s * u.cm**3)\n",
      "\n",
      "def energy_loss_per_n2(T, beta):\n",
      "    L_R = c.k_B * T * beta\n",
      "    return L_R.to(u.erg / u.s * u.cm**3)\n",
      "\n",
      "def net_effective_heating(T_i, T, alpha, beta):\n",
      "    G = energy_gain_per_n2(T_i, alpha)\n",
      "    L_R = energy_loss_per_n2(T, beta)\n",
      "    return G-L_R\n",
      "\n",
      "temperature_list = [5000, 10000, 20000]\n",
      "alpha_dict = {}\n",
      "alpha_dict[5000] = 4.54e-13 * u.cm**3/u.s\n",
      "alpha_dict[10000] = 2.59e-13 * u.cm**3/u.s\n",
      "alpha_dict[20000] = 1.43e-13 * u.cm**3/u.s\n",
      "beta_dict = {}\n",
      "beta_dict[5000] = 3.24e-13* u.cm**3/u.s\n",
      "beta_dict[10000] = 1.73e-13* u.cm**3/u.s\n",
      "beta_dict[20000] = 8.75e-14* u.cm**3/u.s\n",
      "\n",
      "T_i = 2.48e4 * u.K\n",
      "\n",
      "print \"Net effective heating at each temperature (G - L_R):\"\n",
      "for T in temperature_list:\n",
      "    heating = net_effective_heating(T_i, T*u.K, alpha_dict[T], beta_dict[T])\n",
      "    print \"{1:.2e}    at    T={0} K\".format(T, heating)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Net effective heating at each temperature (G - L_R):\n",
        "2.11e-24 cm3 erg / s    at    T=5000 K\n",
        "1.09e-24 cm3 erg / s    at    T=10000 K\n",
        "4.93e-25 cm3 erg / s    at    T=20000 K\n"
       ]
      }
     ],
     "prompt_number": 240
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### (b)\n",
      "Using Table 3.6 for the collision strength = $\\Upsilon$ averaged over the Maxwellian velocity distribution of the electrons, compute the cooling rate for the [O III] $^3$P - $^1$D lines.\n",
      "The collision strength in the table refers to the *total* excitation of the 5007 and 4959 \u00c5 lines (sum of both) from the ground state; ignore the slight difference in excitation energy between these transitions. Assume the low-density limit for these lines. Evaluate the [O III] cooling rate at T=6000, 8000, 10,000 K, assuming that the $^3\\textrm{P}_0 - ^3\\textrm{P}_1$, etc. transitions are collisionally de-excited and thus can be neglected.\n",
      "\n",
      "\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The energy losses from line cooling look like this:\n",
      "\n",
      "$$ L_{line} = n_{ion} n_e q_{12} h \\nu_{line} $$.\n",
      "\n",
      "where \n",
      "\n",
      "$$ q_{21} = \\left( \\frac{2\\pi}{kT}\\right)^{1/2} \\frac{\\hbar^2}{m_e^{3/2}} \\frac{\\gamma_{12}}{g_2} $$\n",
      "$$ = \\frac{8.629 \\times 10^{-6}}{T^{1/2}} \\frac{\\gamma_{12}}{g_2} $$\n",
      "\n",
      "$$ q_{12} = \\frac{g_2}{g_1} q_{21} e^{- h \\nu / k T} $$\n",
      "\n",
      "The abundance of the oxygen 2+ ion goes as\n",
      "\n",
      "$$ n_{O^{++}} = 7 \\times 10^{-4} * 0.2 n_p $$\n",
      "\n",
      "because the O/H ratio is $7\\times 10^{-4}$ and the fraction of oxygen in the relevant ionization state is 20%.\n",
      "\n",
      "So finally,\n",
      "\n",
      "$$ \\frac{L_{21}}{n_p n_e} = \\frac{n_{O^{++}}}{n_p} q_{12} h \\nu_{21} $$\n",
      "\n",
      "with $q_{12}$ defined above.\n",
      "\n",
      "[O III], also known as O$^{++}$, has a collision strength $\\Upsilon = 2.29$ for the $^3$P - $^1$D line.\n",
      "\n",
      "The $^3$P and $^1$D states both have a statistical weight $g = 5$:\n",
      "\n",
      "For $^3$P: \n",
      "$ 2S+1 = 3 => S=1$. The shell P has $L = 1$. $J = S+L = 2$, hence $g = 2J+1 = 5$.\n",
      "\n",
      "For $^1$D: $2S+1 = 1 => S=0$. The shell D has $L=2$. $J=S+L=2$, hence $g = 2J+1 = 5$.\n",
      "\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gamma_OIII = 2.29\n",
      "\n",
      "proton_mass = c.m_p\n",
      "g_OIII = 5 # see above J's, S's, L's.\n",
      "lambda_OIII = 4983 * u.angstrom\n",
      "OIII_abundance = 7e-4 * 0.2\n",
      "\n",
      "def collisional_cooling_per_n2(T, gamma, stat_weight_1, stat_weight_2, wavelength, abundance):\n",
      "#    q21_per_n2 = ( (2*np.pi/(c.k_B*T))**(1/2) * \n",
      "#                    c.hbar**2/mass**(3/2) *\n",
      "#                    gamma/stat_weight_2 )\n",
      "    q21_per_n2 = ( 8.629e-6 / (T/u.K)**(1/2) * \n",
      "                  gamma/stat_weight_2) * u.cm**3 / u.s\n",
      "\n",
      "    photon_energy = c.h * c.c / wavelength\n",
      "\n",
      "    q12_per_n2 = (stat_weight_2/stat_weight_1 * q21_per_n2 * \n",
      "                  np.exp(-photon_energy/(c.k_B*T)) )\n",
      "\n",
      "    cooling_rate = abundance * q12_per_n2 * photon_energy\n",
      "    return cooling_rate.to(u.erg / u.s * u.cm**3)\n",
      "\n",
      "temperatures = [6000, 8000, 10000]\n",
      "\n",
      "def oIII_cooling(T):\n",
      "    return collisional_cooling_per_n2(T*u.K, gamma_OIII, g_OIII, g_OIII, lambda_OIII, OIII_abundance)\n",
      "\n",
      "print \"Collisional line cooling rates for [OIII]:\"\n",
      "\n",
      "for T in temperatures:\n",
      "    cooling = oIII_cooling(T)\n",
      "    \n",
      "    print \"{0:.2e}     at     T={1} K\".format(cooling, T)    \n",
      "\n",
      "#print collisional_cooling_per_n2(10000*u.K, oxygen_ion_mass, gamma_OIII, g_OIII)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Collisional line cooling rates for [OIII]:\n",
        "2.31e-25 cm3 erg / s     at     T=6000 K\n",
        "6.68e-25 cm3 erg / s     at     T=8000 K\n",
        "1.23e-24 cm3 erg / s     at     T=10000 K\n"
       ]
      }
     ],
     "prompt_number": 241
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### (c)\n",
      "Compute the cooling rate for the [OII] $^4$S-$^2$D 3727, 3729 \u00c5 lines using Table 3.7 analogously at the same temperatures."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Here, we use $\\lambda = 3728$ \u00c5, $\\Upsilon$ = 1.34."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "lambda_OII = 3728 * u.angstrom\n",
      "gamma_OII = 1.34\n",
      "g_OII_upper = 4 # from 4S \n",
      "g_OII_lower = 6 # from 2D\n",
      "\n",
      "OII_abundance = 7e-4 * 0.8\n",
      "\n",
      "def oII_cooling(T):\n",
      "    return collisional_cooling_per_n2(T*u.K, gamma_OII, g_OII_lower, g_OII_upper, lambda_OII, OII_abundance)\n",
      "\n",
      "print \"Collisional line cooling rates for [OII]:\"\n",
      "\n",
      "for T in temperatures:\n",
      "    cooling = oII_cooling(T)\n",
      "    \n",
      "    print \"{0:.2e}     at     T={1} K\".format(cooling, T)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Collisional line cooling rates for [OII]:\n",
        "1.19e-25 cm3 erg / s     at     T=6000 K\n",
        "5.16e-25 cm3 erg / s     at     T=8000 K\n",
        "1.21e-24 cm3 erg / s     at     T=10000 K\n"
       ]
      }
     ],
     "prompt_number": 242
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 242
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### (d)\n",
      "Compute the free-free radiative losses at the same temperatures assuming $g_{ff} = 1.3$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The free-free radiation is given by\n",
      "\n",
      "$$ L_{ff} = 1.42 \\times 10^{-27} Z^2 T^{1/2} g_{ff} n_e n_p $$.\n",
      "\n",
      "The bulk of free-free emission should come from hydrogen atoms, so heavier ions are here neglected.\n",
      "\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def free_free_cooling_rate(T, g):\n",
      "    \n",
      "    L_ff = 1.42e-27 * (T/u.K)**(1/2) * g\n",
      "    return L_ff * u.erg * u.cm**3 / u.s\n",
      "\n",
      "def ff_cooling(T):\n",
      "    return free_free_cooling_rate(T*u.K, 1.3)\n",
      "\n",
      "print \"Free-free cooling rates:\"\n",
      "for T in temperatures:\n",
      "    cooling = free_free_cooling_rate(T*u.K, 1.3)\n",
      "    \n",
      "    print \"{0:.2e}     at     T={1} K\".format(cooling, T)\n",
      "    "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Free-free cooling rates:\n",
        "1.43e-25 cm3 erg / s     at     T=6000 K\n",
        "1.65e-25 cm3 erg / s     at     T=8000 K\n",
        "1.85e-25 cm3 erg / s     at     T=10000 K\n"
       ]
      }
     ],
     "prompt_number": 243
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### (e)\n",
      "\n",
      "Compute the total cooling rates at the three temperatures by estimating that the [NII] cooling is equal to the [O III] cooling rate computed above. Compare your results with the heating rates to estimate the resulting nebular temperature in thermal equilibrium.\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$ L = L_{[OIII]} + L_{[OII]} + L_{ff} + L_{[NII]}$$\n",
      "$$ = 2L_{[OIII]} + L_{[OII]} + L_{ff} $$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "print \"Total cooling = 2*OIII_cooling + OII_cooling + ff_cooling\"\n",
      "for T in temperatures:\n",
      "    \n",
      "    total_cooling = 2*oIII_cooling(T) + oII_cooling(T) + ff_cooling(T)\n",
      "    \n",
      "    print \"{0:.2e}     at     T={1} K\".format(total_cooling, T)\n",
      "\n",
      "print \"\"\n",
      "print \"Net effective heating at each temperature (G - L_R) from above:\"\n",
      "for T in temperature_list:\n",
      "    heating = net_effective_heating(T_i, T*u.K, alpha_dict[T], beta_dict[T])\n",
      "    print \"{1:.2e}    at    T={0} K\".format(T, heating)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Total cooling = 2*OIII_cooling + OII_cooling + ff_cooling\n",
        "7.25e-25 cm3 erg / s     at     T=6000 K\n",
        "2.02e-24 cm3 erg / s     at     T=8000 K\n",
        "3.85e-24 cm3 erg / s     at     T=10000 K\n",
        "\n",
        "Net effective heating at each temperature (G - L_R) from above:\n",
        "2.11e-24 cm3 erg / s    at    T=5000 K\n",
        "1.09e-24 cm3 erg / s    at    T=10000 K\n",
        "4.93e-25 cm3 erg / s    at    T=20000 K\n"
       ]
      }
     ],
     "prompt_number": 244
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Based on the above output, the cooling and heating curves seem to balance near $ 1 \\times 10^{-24}$, which is reached somewhere between 6000 and 10000 K \u2013 it is hard to pinpoint exactly where since the heating values are populated sparsely.\n",
      "\n",
      "However, I can make a plot and"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import matplotlib.pyplot as plt\n",
      "\n",
      "fig = plt.figure()\n",
      "\n",
      "temp_array = np.linspace(5000, 20000, 50)\n",
      "\n",
      "def total_cooling(T):\n",
      "    return 2*oIII_cooling(T) + oII_cooling(T) + ff_cooling(T)\n",
      "\n",
      "cooling_array = total_cooling(temp_array)\n",
      "\n",
      "plt.plot(temp_array, cooling_array, 'b-', label=\"Total cooling curve\")\n",
      "plt.semilogy()\n",
      "plt.xlim(5000, 20000)\n",
      "plt.xlabel(\"Nebula temperature (K)\")\n",
      "plt.ylabel(\"Heating or cooling rate / $n_p n_e$ (erg cm$^3$ s$^{-1}$)\" )\n",
      "\n",
      "heating_list = [2.11e-24, 1.09e-24, 4.93e-25]\n",
      "plt.plot(temperature_list, heating_list, 'ro-', label=\"Effective heating (crudely interpolated)\")\n",
      "plt.legend()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 245,
       "text": [
        "<matplotlib.legend.Legend at 0x109259090>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAESCAYAAAC7NAEnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XdYFFf3B/DvglhBBFEUdhHEAgiKAtaomxjFLtiiYCcm\nmp9izSuxvGBL1MQkkLwpaqLGFmskiiARxRIVRFFjIQgRaRakCEhf7u+PCSsLCyy6w+zC+TzPPLCz\nszNnFtjDvXPmXhFjjIEQQgipYzpCB0AIIaRhogRECCFEEJSACCGECIISECGEEEFQAiKEECIISkCE\nEEIEQQmIEEKIICgBEUIIEUQjoQN4XYGBgQgKCkJ2dja8vLyQl5en8Hjo0KFCh0gIIaQaIm0fCSEr\nKwvLly/Hjh07lD4mhBCimQTvgpszZw5MTU3h4OCgsD4kJAQ2Njbo3LkzNm/eXOXrN2zYgAULFlT5\nmBBCiGYSPAHNnj0bISEhCutkMhkWLFiAkJAQ3Lt3DwcOHMD9+/exZ88eLFmyBKmpqWCMYcWKFRgx\nYgQcHR0rPSaEEKLZBL8GNHDgQCQkJCisi4yMRKdOnWBpaQkAmDJlCgIDA+Hj44Pp06cDAAICAhAW\nFobs7GzExcWhqKhI4fGHH35Yx2dCCCGkNgRPQMqkpKRAIpHIH4vFYkRERChs4+3tDW9vb4V1Cxcu\nrJP4CCGEvDmNTEAikUir9ksIIfUdH/Vqgl8DUsbc3BxJSUnyx0lJSRCLxWrZN2NM7Yuvry8v++Vz\noZgp3voQs7bFq60x80UjE5CzszMePHiAhIQEFBUV4eDBgxg7dqzQYRFCCFGjWieggoICFBYWqi2A\nqVOnon///oiNjYVEIsHOnTvRqFEjfPvtt3B1dYWdnR3ee+892NraquV4fn5+CA8PV8u+CCGkPgsP\nD4efnx9v+6/xRtTS0lIcP34cBw4cwOXLl1FaWgrGGHR1ddGvXz94enrCzc1NK66viEQiXpqT4eHh\nkEqlat8vnyhm/mlbvID2xaxt8QLaGTNfn501JqBBgwZh4MCBGDt2LBwdHdGkSRMAQGFhIaKjo/H7\n77/j0qVLuHDhgtqDUze+3kRCCKnPBEtAhYWF8qTzJttogvqQgLShpUkI0V7KPiP5+uyssQxblcSi\nDcmnjJ+fH6RSqdY1gcvT9iRKCNFMFf/BDQ8P5/WaudYPRlob9aUFpO3nQAjRTFV9vgjSAiobc63M\nuXPnMG3aNLUHQQghpOGptgz72rVrWLhwIXbs2IEdO3bg1KlTdRUXIUrp6Ojgn3/+qbPjSaVS/PTT\nTwCAffv2wdXVtc6OTUh9V20LaNy4cejTpw/atWsHAHj27FmdBMWn+nANSBPp6+vL+49fvnyJpk2b\nQldXFwCwbds2TJ06tdJrwsPDMX36dIVRLzSNSCSSn5enpyc8PT0FjoiQusP3NaAaixDKkg8AtG3b\nlrdA6gqfN1U1ZLm5ufLvrays8NNPP+Gdd94RMKKGp6y7nColibqU/bO+du1aXvav8kgIV65c4SUA\nUr8VFhZi8eLFMDc3h7m5OZYsWYKioiK8fPkSI0aMQGpqKgwMDNCyZUs8efIEkZGR6NevH4yMjGBm\nZoaFCxeiuLhYpWNlZGRg9uzZMDc3h7GxMdzd3eXPbd++HZ07d0br1q0xbtw4PH78WP7c5cuX4eLi\nglatWqF3795V/q7v2rULAwcOlD/W0dHBjz/+iC5dusDIyEhhIsTS0lIsW7YMbdq0QceOHfHtt99C\nR0cHpaWlSvedlJSE8ePHo23btjAxMZGP7O7n5yefggQAEhISFPYjlUqxevVqDBgwAC1atMDnn38O\nFxcXhX1/9dVXGDdunPznsXz5cnTo0AHt2rXD/PnzUVBQoNL7S4i6qZyAsrOz+YyD1FMbN25EZGQk\nbt26hVu3biEyMhIbNmxAixYtEBISAjMzM+Tk5CA7Oxvt2rVDo0aN4O/vj/T0dFy5cgVhYWH47rvv\nVDrW9OnTUVBQgHv37uHZs2dYunQpAODs2bNYuXIlDh8+jMePH6NDhw6YMmUKAC5pjRo1CosXL0ZG\nRgaWLl2KUaNGITMzU6VjBgUFISoqCrdv38ahQ4dw+vRpAFy3Y0hICG7duoUbN27g+PHjVbZMZDIZ\nRo8eDSsrKzx69AgpKSnyLktVWjN79+7Fjh07kJubi3nz5uHvv/9GXFyc/Pn9+/fLuw59fHwQFxeH\nW7duIS4uDikpKVi3bp1K50qI2jEVhYSEqLqpxqrF6WosbTgHS0tLFhYWxhhjzNramgUHB8ufO336\nNLO0tGSMMXbu3DkmFour3ddXX33F3N3d5Y9FIhGLj4+vtF1qairT0dFhWVlZlZ6bM2cOW7Fihfxx\nbm4u09PTYwkJCeyXX35hffr0Udi+X79+bNeuXYwxxqRSKfvpp58YY4zt3LmTvfXWWwqx/Pnnn/LH\nkydPZps3b2aMMfb222+zbdu2yZ87c+YME4lETCaTVYrv8uXLrE2bNkqf8/X1ZdOmTZM/fvjwocJ+\npFIp8/X1VXjNtGnT2Lp16xhjjMXGxjIDAwOWn5/PSktLWYsWLRTev8uXLzMrK6tKxyUNU1WfL3x9\n7mjkaNh8qu+DkYpE6lnUJTU1FR06dJA/trCwQGpqapXbx8bGYvTo0Wjfvj0MDQ2xatUqpKen13ic\npKQkGBsbw9DQsNJzZa2eMi1atEDr1q2RkpKCx48fw8LCQmH7Dh06VBtjeeWvkTZv3lx+Lezx48eV\nJlWsLvYOHTpAR+f1/hzLHwcAPDw8cODAAQBc68fd3R1NmzZFWloa8vLy4OTkBCMjIxgZGWHEiBF4\n/vz5ax2X1H98D0aq8m+8g4MDb0HUpbIquPqKMfUs6mJmZqYw5XpiYiLMzMwAKO9emj9/Puzs7BAX\nF4cXL15g48aNVV43KU8ikSAjIwMvXryoMYaXL18iPT0dYrEYZmZmePTokcL2jx49grm5uYpnqFz7\n9u0rzWlVXeyJiYmQyWSVntPX10deXp788ZMnTyptU/F9fPfdd5GWloZbt27h119/hYeHBwDAxMQE\nzZo1w71795CZmYnMzExkZWVR9zqpklQq1YwEdOnSJfkv6vr16+Hu7o4bN27wFhipH6ZOnYoNGzbg\n+fPneP78OdatWye/qG5qaor09HSFD8Dc3FwYGBigefPmiImJwffff6/Scdq3b48RI0bgo48+QlZW\nFoqLi+UD5E6dOhU7d+7ErVu3UFhYiJUrV6Jv376wsLDAiBEjEBsbiwMHDqCkpAQHDx5ETEwMRo8e\nXetzZeUm75o8eTL8/f2RmpqKrKwsbN68ucrrOX369EH79u3h4+ODvLw8FBQU4PLlywAAR0dHXLhw\nAUlJSXjx4gU+++wzpcctT09PD5MmTcLy5cuRmZmJoUOHAuCKJubOnYvFixcjLS0NAJCSkoLQ0NBa\nnysh6qByAlq/fj1atmyJS5cuISwsDF5eXpg/fz6fsZF6YPXq1XB2dkb37t3RvXt3ODs7Y/Xq1QAA\nGxsbTJ06FR07doSxsTGePHmCL774Avv370fLli3xwQcfYMqUKQof3NVdlN+zZw/09PRgY2MDU1NT\nBAQEAACGDBmC9evXY8KECTAzM8PDhw/x66+/AgBat26NkydPYuvWrTAxMcEXX3yBkydPwtjYuNL+\ny98TpCyW8s/PnTsXw4YNQ/fu3eHk5IRRo0ZBV1dXaTebjo4OTpw4gbi4OFhYWEAikeDQoUMAuNbM\ne++9h+7du8PFxQVjxoxRetyKPDw8EBYWhkmTJikcc/PmzejUqRP69u0LQ0NDDB06FLGxsVW+p4Tw\nSeWx4BwdHXHz5k34+PjAwcEBnp6e6NmzJ6Kjo/mOUW3qwzhq9eEcGqLg4GDMnz9foSuQEE1T12PB\nqdwCMjc3xwcffICDBw9i1KhRKCgoUKlvnpCGqKCgAKdOnUJJSQlSUlKwdu1ajB8/XuiwCNEoKieg\nQ4cOwdXVFaGhoWjVqhUyMzPx+eef8xkbL+p7FRzRDIwx+Pn5wdjYGL169UK3bt3ofhuidQSfkrs+\nqQ/dV/XhHAghmklju+AIIYQQdaIERAghRBCUgAghhAiCEhAhhBBBUAIihBAiCJUTUGxsLAoLC/mM\npU5QGbZwVq9ejTZt2sjHgvvtt98gkUhgYGCAW7duqe04Fy9ehI2Njdr2V6bi3Dx84uscyvz4449Y\nsmQJb/svr/y05jV53SnX58+fjw0bNtT6dUKYNWsW1qxZUyfHKj+HVWFhIWxtbWs1+CzfZdjVjrH9\nySefMC8vL/bDDz+w5cuXs1WrVvEyJHddqeF0tYImn0OHDh1Ys2bNmL6+vnxZuHAhY4yxR48esWbN\nmrHnz5/Lt+/YsSP7/fff3/i4VU3RoG5+fn4KUyOoU12dA2OMFRYWMolEwlJTU+vkeOWntKhJXbwP\nM2fOZKtXr+b1GNWZNWsWW7NmjUrbDh48mO3YseO1j1VxCpEtW7awZcuWVbl9VZ8vfH3uVDsl97Bh\nw9C5c2fk5ORg2rRpNPgoqZZIJMLJkyeVTsWdmJiI1q1bo3Xr1gC4GzUTExNhZ2enlmOzOrg3iu9j\n1MU5AEBgYCBsbW3Rvn37Wr+2tLT0taeNqC9kMhl0dXXfaB+q/qzVPb361KlT0bNnT3z22WfQ09NT\n675fR7W/SSYmJrh27RpsbGzwv//9743fdMKfC0FBWO3qCj+pFKtdXXEhKEiQfShz5swZDBs2TD79\ntoeHB1q2bAmZTIYePXqgc+fOALi5gyZMmIC2bduiY8eO+Oabb+T7KC0txaeffopOnTqhZcuWcHFx\nQXJyMgYNGgQA6NGjBwwMDHD48GGEh4fL58jZvHkzJk2apBDPokWLsGjRIgDAixcv4OXlBTMzM4jF\nYqxZs6bKIaZEIhGKioowc+ZMtGzZEvb29rh+/br8+erir26q8ZrOAQAsLS2xdetW9OjRA61atcKU\nKVMUusS3bNkiP4cdO3ZU25UVHByMwYMHK6y7dOkS+vfvDyMjI1hYWOCXX34BwHUXzZ8/HyNHjoS+\nvj7OnTtXqUut4lTlf/zxB2xsbNCqVSssXLhQYZRwAPj5559hZ2cHY2NjDB8+HImJiZVivHbtGtq1\na6fwumPHjsHR0VHpOZXv1goPD4dYLMaXX34JU1NTmJmZYdeuXQC4mWr379+PLVu2wMDAQD5VeXU/\nOz8/P0ycOBHTp0+HoaEhdu3aJV83ZcoUtGzZEk5OTrh9+7b8Nffv34dUKoWRkRHs7e1x4sQJpXFn\nZmZi9OjRaNu2LYyNjTFmzBikpKQAAFatWoWLFy9iwYIFMDAwgLe3NwAgJiYGQ4cORevWrWFjY4PD\nhw/L95eeno6xY8fC0NAQffr0QXx8vMLxxGIxjIyMqpx2vs6p2lSKjIxkbm5uzNHRkdnb2zN7e3vm\n4ODAR6uMN7U4XY2l7BzOnzzJVlpbK0zps9Lamp0/eVLl/apjH5aWluzMmTNKnwsPD680+2n57haZ\nTMZ69erF1q9fz4qLi9k///zDOnbsyE6fPs0Y47oOHBwcWGxsLGOMsVu3brH09PRK+2FMcabVhIQE\n1rx5c5aTk8MYY6ykpIS1b9+eRUREMMYYc3NzY/PmzWN5eXns2bNnrHfv3uzHH39Ueg6+vr6sadOm\nLDg4mJWWlrJPPvmE9e3bV6X4r1+/ziIiIphMJmMJCQnM1taWff3110rfi4rnUPbe9unThz1+/Jhl\nZGQwW1tb9sMPPzDGGAsODmbt2rVj9+7dY3l5eczT05Pp6OhU2ZXl4uLCjhw5In+ckJDADAwM2K+/\n/spKSkpYeno6u3nzJmOM664yNDRkly9fZowxVlBQUKlLrXw3T1paGjMwMGBHjx5lJSUl7KuvvmKN\nGjWSb3/8+HHWqVMnFhMTw2QyGduwYQPr37+/0vfBzs5OYTZdNzc39uWXXyo9p/LdWufOnWONGjVi\nvr6+rKSkhJ06dYo1b95cPltuxS6wmn52vr6+TE9PjwUGBjLGGMvPz5evKzvPL774gllZWbGSkhJW\nVFTErK2t2WeffcaKi4vZ2bNnmYGBAfv777/lxy/rAkxPT2fHjh1j+fn5LCcnh02aNIm5ubnJY6v4\nXufm5jKxWMx27drFZDIZi46OZiYmJuzevXuMMcbee+899t5777G8vDx2584dZm5uzgYOHKjwXo0d\nO5YFBAQofR+r+ozk67NT5b127tyZBQYGsvj4ePbw4UP5ok3qawJaNWyY0nnlVru6qrxfdeyjQ4cO\nTF9fn7Vq1Uq+lPVfK5t+u/yHzdWrV5mFhYXC859++imbPXs2Y4yxLl26VHm9qKYP77feeov98ssv\njDHGQkNDmbW1NWOMsSdPnrAmTZqw/Px8+bb79+9nb7/9ttLj+Pr6sqFDh8of3717lzVr1kyl+Cuq\naapxZQlo37598sf/+c9/2Lx58xhjjM2ePZutXLlS/lxcXFy111I6d+4s/3Ati3P8+PFKt501axab\nOXOmwrrqEtDu3btZv379FLYXi8Xy7YcPH67wWplMxpo3b84SExMrvQ+bNm1inp6ejDHug7p58+bs\nyZMnVcZZ9qF+7tw51qxZM4Upztu2bSv/p6P8tozV/LPz9fVlgwcPVnje19dX4TxLS0tZ+/bt2cWL\nF9mFCxdYu3btFLafOnUq8/PzU3r88qKjo5mRkZH8sVQqVbgG9Ouvv1ZKKB988AFbu3YtKykpYXp6\nevJExxhjK1euVLgGxBhjnp6e8inbK6rrBFTtNaDyTExMMHbsWD4aYeQNNaqiOlH39GmV59eu6hdB\nt6BA5ThEIhECAwOVXgOqyaNHj5CamgojIyP5OplMJu+eSk5OhrW1da33C7yaonr69OnYv38/PD09\n5ccsLi5WuBZSWlpaaYru8kxNTeXfN2/eXD4qfE3xx8bGYunSpbh+/Try8vJQUlICZ2fnWp1H+em/\nmzVrhsePHwPgpv/u3bu3/Lnqpv8GACMjI4VJAJOTk9GxY8cqt6845Xd1UlNTKx2//OsfPXqERYsW\nYdmyZQrbpKSkVDqOp6cnunXrhry8PBw6dAiDBg1SeP+r07p1a4VrVeWnS6+opp8doPw9Lb9OJBJB\nLBbLp3KveC5VTfOel5eHJUuW4PTp08jMzATATcrIGJNf/yl/HejRo0eIiIhQiLWkpAQzZszA8+fP\nUVJSonBsZb/LOTk5Cq8XksoJyM/PD15eXnj33XfRuHFjANwbQ0PMC6+kSROl62WurkBIiGr7cHUF\nlMyMKWva9I1iU5WFhQWsrKyqnBxNIpEgLi7utYoWJk6ciGXLliElJQXHjx/H1atX5fts0qQJ0tPT\nVbqwXt0FYYlEUm388+fPh5OTEw4ePIgWLVrg66+/xtGjR2t9LsrUZvpvAOjevbtCnBKJBJGRkSof\nr0WLFnj58qX8cflpws3MzBAYGCh/zBhTiMfCwgJr1qzB1KlTazyOWCxG3759cezYMezduxcfffRR\ntduresG+4nY1/e5VnIiwTPnzKi0tRXJyMszNzeXnXD6JPHr0SKGsvmz91q1bERsbi8jISLRt2xY3\nb95Er1695K9VFuvgwYOVzmIrk8nQqFEjJCYmomvXrgCg9Pra/fv3sXz5cqXnWtdULmfZvXs3bt26\nhZCQEJw8eRInT56s8sIaqVvDvL2xqkLrYKW1NYYuXFin+wBev5Krd+/eMDAwwJYtW5Cfnw+ZTIY7\nd+4gKioKAPD+++9jzZo1iIuLA2MMt2/fRkZGBgCuVVLxYmt5bdq0gVQqxaxZs9CxY0f5H2f79u0x\nbNgwLF26FDk5OSgtLUV8fLx8Ku/anFtN8dc01XhN51BdPJMnT8bOnTsRExODvLw8rF+/vtrXjRw5\nEufPn5c/9vT0xJkzZ3D48GGUlJQgPT1dfl+WsnN2dHTEsWPHkJ+fj7i4OIWChJEjR+Lu3bv47bff\nUFJSgoCAAIUENW/ePHz66ae4d+8eAK4IpPxF9IpmzJiBzZs3486dO9X+s8sqFDpUx9TUVKFAo6af\nXVX7vX79uvw8v/76azRt2hR9+/ZF79690bx5c2zZsgXFxcUIDw/HyZMnMWXKlEqx5ubmolmzZjA0\nNERGRgbWrl1bKdbyvxejR49GbGws9u7di+LiYhQXF+PatWuIiYmBrq4uxo8fDz8/P+Tn5+PevXvY\nvXu3QhJLSUlBRkYG+vbtq9J7xTeVE1BUVBSuXbuG3bt3Y+fOnfKFCG/QqFFw9ffHGldX+A0ejDWu\nrhju749Bo0bV6T4AYMyYMTAwMJAvEyZMkD9X3VTSOjo6OHnyJG7evImOHTuiTZs2+OCDD+RdRUuX\nLsXkyZMxbNgwGBoaYu7cuSj4t3vQz88PM2fOhJGREY4cOaL0P8eyKao9PDwU1v/yyy8oKiqSV2VN\nmjRJ4QOzYrxVnYOurm618dc01bgq51BVLMOHD4e3tzfefvttdOnSBf369QMANKmiZTx69GjExMTI\nu/AkEglOnTqFrVu3onXr1ujZs6e8oktZHEuWLEHjxo1hamqK2bNnY9q0afJtTExMcPjwYfj4+MDE\nxARxcXF466235K91c3PDihUrMGXKFBgaGsLBwQGnT5+u9H6WGT9+PBITE+Hu7o6m1bTGa5ouvTwv\nLy/cu3cPRkZGGD9+fI2/e1X93MeNG4eDBw/C2NgY+/btw7Fjx6Crq4vGjRvjxIkTCA4ORps2bbBg\nwQLs2bMHXbp0qbS/xYsXIz8/HyYmJujfvz9GjBihcKxFixbhyJEjMDY2xuLFi6Gvr4/Q0FD8+uuv\nMDc3R/v27fHJJ5+gqKgIAPDtt98iNzcX7dq1w5w5czBnzhyFuPfv349Zs2ZpRAk2ANWvLM2aNYvd\nuXNHrReg6hoA5uvry86dOyd0KK+tFj8y0kDdu3eP6erqKlyEr2jbtm1s8eLFdRjV6+vUqRMLCwsT\nOgwFfN6UzJeCggJmY2PD0tLSqtym4ufLuXPnmK+vL2+fOypPSGdjY4P4+HhYWVnJ/7MSiUQKte+a\nrj5M5lYfzoGo32+//YaRI0ciLy8PM2fORKNGjXDs2DGhw3pjx44dg4+PT5XXZ4Ti5+eH+Ph47Nmz\nR+hQ1KquJ6RTuQghRMWL2YSQurdt2zbMnj0burq6kEql+O6774QO6Y1JpVLExMRo5Id8TV2kRDU0\nJbeWqQ/nQAjRTBo7JfeMGTPkdeoAkJGRUekCFyGEEKIqlRPQ7du3FW5eMjY2psFJCSGEvDaVExBj\nTH7fBcC1gGQyGS9BEUIIqf9ULkJYtmwZ+vXrh8mTJ4MxhsOHD2PVqlV8xkYIIaQeq1URwt27d3H2\n7FmIRCK88847apvLpa7Uhwv4VHlDCOFTXRYhUBUcIYSQagleBUcIIYSoEyUgQgghgqAERAghRBAq\nV8Ft3bpVoR9QJBLB0NAQTk5OVc7Tron8/PwglUohlUqFDoUQQjRaeHg4wsPDedu/ykUIHh4eiIqK\nwpgxY8AYQ1BQEBwcHPDo0SNMnDgRK1as4C1IdaEiBEIIqT3Bq+AGDhyI4OBg6OvrA+AmUho5ciRC\nQkLg5OSE+/fvqz04daMERAghtSd4FVxaWpp8Km4A0NPTw9OnT9G8efNqJ4oihBBClFH5GpCnpyf6\n9OkDNzc3MMZw4sQJeHh44OXLl1p3QyohhBDhqdQFxxhDUlISnj59ij///BMikQgDBgyAs7NzXcSo\nNtQFRwghtSfoNSDGGBwcHHDnzh21B1CXKAERQkjtCXoNSCQSwcnJCZGRkWoPgBBCSMOkchVc165d\nERcXhw4dOqBFixbci0Ui3L59m9cA1YlaQIQQUnuCl2EnJCQoDcTS0lLtQfGFEhAhhNSe4GXYFhYW\nuHjxInbv3g1LS0vo6Ojg2bNnag+IEEJIw6ByC2jevHnQ0dHB2bNnERMTg4yMDAwbNgxRUVF8x6g2\n1AIihJCayWRARgaQlgY8ewa8/TY/n50q3wcUERGB6Oho9OzZEwBgbGyM4uJitQdECCFE/V6+BJ4+\n5RLKs2evvi+/rizhZGQALVsCbdsCbdrwF5PKCahx48aQyWTyx2lpadDRocG0CSFEKIWFXAJ5/Bh4\n8kRxefr01fLkCdeqMTXllrZtucXUFOjYEejbl/u+TRtufevWgJ7eq+PwNRGzyglo4cKFcHd3x7Nn\nz7By5UocOXIEGzZs4CcqQghpwHJzgdRUbnn8uPLy5An3NTeXSxzt2nFL+/bcY3t7YMiQV8+ZmgIG\nBvwlktdVqym579+/j7CwMADAkCFDYGtry1tgfKBrQIQQIRUXc4kjJeXVUpZoyn9fUgKYmXEJpbrF\n2Bioi44owcqwGWMQ1ZA2VdlGnQIDAxEUFITs7Gx4eXlBIpHA398f6enpcHV1hZeXl9LXUQIihPAl\nP59LIsnJ3JKU9Or75GTuufR0rovL3LzyYmb2ajE01KzWimAJaPDgwRg9ejTGjRuHLl26KDz3999/\n4/jx4wgKCsKFCxfUHlxNsrKysHz5cuzYsQMAUFpaiilTpuDQoUNKt6cERAh5HTIZ13JJTOSWpKTK\nX7OzueQhkQBiceXF3JzrCmuk8oUPzcHXZ2eNb0VoaCj27duH//u//8OdO3dgYGAAxhhyc3Nhb28P\nT09PnDlz5rUOPmfOHAQFBaFt27b466+/5OtDQkKwePFiyGQyvP/++1VOdrdhwwYsWLAAAHDixAl8\n9913mDt37mvFQghpuAoLuUSSkMAtiYnAo0fckpjIdYu1bg1YWHCLRAJYWwNvv819L5FwLRuqy6qd\nWl0DkslkeP78OQDAxMQEurq6b3TwixcvQl9fHzNmzJAnIJlMhq5du+LMmTMwNzeHi4sLDhw4gKio\nKNy4cQMff/wx2rdvDx8fHwwbNgxDhgxR2Oe4ceMQGBio9HjUAiKkYSou5lop//wDPHz4KtGULc+f\nc60US0ugQ4dXi4UF91UsBpo0EfQUBCVYC6g8XV1dmJqaqu3gAwcOlA/xUyYyMhKdOnWSD/EzZcoU\nBAYGwsfHB9OnTwcABAQEICwsDNnZ2YiLi4ONjQ2OHTuGgoICvP3222qLjxCiHRjj7l355x8gPl7x\n68OHXPfgeoGsAAAgAElEQVRZu3aAlRVXdmxlBQwfziUcS0uu6+wN/58mr0HjeiNTUlIgkUjkj8Vi\nMSIiIhS28fb2hre3t8K6wYMH10l8hBBhMMaVH8fFccuDB9zX+HhuYYzrFrO25pJM797A1KlcspFI\ngHITOhMNoXEJiO9qOj8/P/n3UqkUUqmU1+MRQlTHGNcdFhvLLQ8evPoaHw+0aAF06vRqcXfnvnbs\nyJUka1LlmDYLDw9HeHg478fRuARkbm6OpKQk+eOkpCSIxWK17b98AiKECKOggGu9xMQAf//96uuD\nB9zzXbsCXbpwy6RJQOfOXKJp2VLYuBuKiv+cr127lpfjaFwCcnZ2xoMHD5CQkAAzMzMcPHgQBw4c\nEDosQshryMgA7t9XXP7+m7snxsqKSzQ2NsA77wDz53MJp3Vrask0FG+cgE6dOgVra2vk5eXJBypV\n1dSpU3H+/Hmkp6dDIpFg3bp1mD17Nr799lu4urpCJpPBy8tLrSMu+Pn5UdcbIWr2/Dlw9+6r5d49\nLtnk5XEJxs4OsLUFBg3ivlpZKY41RjQT311xtSrDVmbXrl1wcHCASCRCr1691BUXL6gMm5A3k5MD\n3LkD/PUX97Us4eTnA926vVrs7LjF3JxaM/WBRpRhK9OnTx9s2bIF77zzjsYnIEKIamQy7uL/7dtc\nsvnrL+77Z8+4FoyDAzfg5ahRXMKhRENex2u3gK5cuQKxWKxQMq3pqAVESGXZ2VxyuXULuHmT+3r3\nLjfYZffuXLIp+2ptTffLNEQa0QLasGEDHjx4gEaNGmHo0KGIjIzEokWL1B4Un+gaEGnI0tKAGzeA\n6Gju640b3E2a9vZAjx5Az57ArFlcsqGKM6JR14B+++03uLu748WLFzh16hQMDAwwevRo3oJTN2oB\nkYbk2TMgKopbrl/nkk1ODpdkevXilp49uUo0atWQ6gg2GnZ5v/32G8RiMVxcXNQeSF2gBETqq6ws\n4Nq1Vwnn2jUu2Tg7c4uTE5dwrKzoWg2pPY1IQIsXLwYAxMfHo2nTphg8eLB8NGptQAmI1AfFxdw1\nm4iIV0tKCteacXHhFmdn7noNJRuiDhqRgC5evAiRSIS33noL+fn5uHv3LpydndUeFF8oARFt9OQJ\ncOUKcPky9zU6mmvJ9OnzaunWTTvnmSHaQSMSkLYTiUTw9fWlIgSisWQyruT58uVXS1YW0K/fq8XF\nhQoESN0oK0JYu3atZiagkJAQPHnyBB4eHmis4cPNUguIaJqCAiAyErh0Cbh4kWvhtGsHvPUW0L8/\nl3C6dqWJzoiwNLYFdOXKFRgZGeHixYsaPxspJSAitNxcrlUTHg5cuMDdd2NnBwwcyCWdt94C2rQR\nOkpCFGnEfUDKhISEoFWrVmjZsiUyMjJgbGysjrgIqRfKJ5zwcK54oFcvQCoF/PyAvn0BfX1hYyRE\nKG/cAoqKioKFhQXCwsIQFBSEvXv3qis2taMWEOFbURFXlXbmDBAWxrVwyhKOVMolnObNhY6SkNrR\nuC44ZUPxFBQUoGnTpmoLTt2oCIGoG2Nc0cCZM9xy6RI3pcCQIcC77wIDBlDCIdpLo4oQKg7F8/Tp\nU60aiodaQEQd0tKAP/4ATp/mFn19YOhQLuFIpdx8NoTUJxpxDahbt25YvXq1fCgea2trtQdEiKYp\nKQGuXgVCQriEExvLJZrhwwFfX246aEJI7dFQPIQokZ7OJZygIC7pWFhwCWf4cK40WsPvOCBErTTi\nGhANxUPqK8a4CdZOnuSSzu3bwNtvA6NHAyNHcvPdENJQaUQX3IQJEyoNxaNtaDoGUkYm40qkjx/n\nFpkMGDMGWLMGGDwY0OB6GkLqhEZNx6DtqAVE8vO5arXjx4ETJ7iWjZsbt3TvToN3EqKMRnTBaTtK\nQA3Ty5dAcDBw+DB3XadXLy7hjBsHWFoKHR0hmk+wBHT58mX069cPonrwryEloIbj5UvuWs7hw0Bo\nKDdi9MSJgLs7DXVDSG0JloDmzZuHiIgIdOnSBSNGjMDw4cPRrl07tQdSFygB1W8FBVxLZ/9+Lun0\n6wdMmsS1dExMhI6OEO0leBfc/fv3ERwcjNDQUGRlZeGdd97B8OHDMWDAAOhqyXy+lIDqH5kMOHeO\nSzrHjwOOjoCHBzB+PEDDEhKiHoInoPLy8vJw7tw5BAcH48qVK7h+/braA+MDJaD6gTHg+nVg717g\n4EGukMDDA3jvPSqXJoQPGpWAtBUlIO32+DGXdHbt4rrbpk3jEk/XrkJHRkj9phH3AdUHdB+Qdiko\nAAIDgd27ucnaxo8HfviBmzenHtTFEKLR6D4gNaIWkPa4dQvYvh04cIArm545k6tga9FC6MgIaXj4\n+uxUeaLf0tJS7NmzB+vWrQMAJCYmIjIyUu0BkYYrNxf46SeuZHr0aK5yLTqaG3l62jRKPoTUNyq3\ngObNmwcdHR2cPXsWMTExyMjIwLBhwxAVFcV3jGpDLSDNdOMGsG0bcOgQMGgQMHcuN+inlhRXElLv\nCX4NKCIiAtHR0ejZsycAwNjYGMXFxWoPiDQMhYXAkSPAN98AT55wSefOHcDMTOjICCF1ReUE1Lhx\nY8hkMvnjtLQ06Oio3INHCAAgORn48Ufu+k737sDKlcCoUdTaIaQhUjmDLFy4EO7u7nj27BlWrlyJ\nAQMG4JNPPuEzNlJPMAZcuMCNStC9O5CVBYSHc6MVjB1LyYeQhqpWVXD3799HWFgYAGDIkCGwtbXl\nLTA+0DWgulVczHWzbd0K5OQACxcCM2YALVsKHRkhpDYEvxF1xYoV2Lx5c43rNBkloLqRnQ3s2AH4\n+3OjTS9fznWzUY8tIdpJ8DLs0NDQSutOnTql1mDqgp+fH683VjVkycnAf/4DWFkBkZFc6+f8eW6S\nN0o+hGif8PBw+Pn58bb/GltA33//Pb777jvEx8fD2tpavj4nJwcDBgzAvn37eAtO3agFxI/4eGDT\nJuDoUe6G0UWLaJ4dQuoTwbrgXrx4gczMTPj4+GDz5s3yIAwMDNC6dWu1B8QnSkDqdfcu8Nln3CRv\n//d/gLc3oGW/EoQQFQh+DQgAMjMz8eDBAxQUFMjXDRo0SO1B8YUSkHrcuAFs3AhcugQsWQJ89BEV\nFhBSnwl+I+r27dsREBCA5ORkODo64urVq+jXrx/Onj2r9qCIZrpxA/jvf4GbN4GPPwb27AGaNxc6\nKkKItlL50rC/vz8iIyPRoUMHnDt3DtHR0TA0NOQzNqIh7twBJkzgigmGD+eu+SxaRMmHEPJmVE5A\nTZs2RbNmzQAABQUFsLGxwd9//81bYER4Dx4Anp7AkCFA//7c4wULgCZNhI6MEFIfqNwFJ5FIkJmZ\nCTc3NwwdOhRGRkawpFKneikxEVi7Fvj9d2DxYm7+HQMDoaMihNQ3rzUfUHh4OLKzszF8+HA0btyY\nj7h4IRKJsGrYMAzz9sagUaOEDkfjvHjBVbVt3w7Mnw8sWwYYGQkdFSFEaIIWITDGkJycDIlEAgBa\nPZvohtBQrIqPBwBKQv8qKuIGCN2wgZuH5/ZtwNxc6KgIIfWdyteARowYwWccdWpjfDz++OYbocMQ\nHGPczaPdugGnTgFnznATwlHyIYTUBZUSkEgkgpOTU72aAVX3zh3gl1+Av/4CSkqEDqfORUUBAwcC\n69YB330HBAcDDg5CR0UIaUhUvgbUtWtXxMXFoUOHDmjx79zIIpEIt2/f5jVAdRKJRCg72TVdumB9\nz57cnM/JyYCdHdCz56ule/d6WWf8/Dk3B8+JE9zNpDNn0nQIhJDqCX4j6unTp9V+cCH4AXhgZoYP\nv/ySG6IZ4OYKuH2bS0ZRUdxV+JgYbkCz8kmpZ0/A2FjA6F9fSQl3nWftWq60OiYGoNu4CCHVCQ8P\n53Xw5teqgtNWIpEIq11dMXThwpoLEIqKgHv3uKRUtty6xZWFVUxKYjEgEtXNSbyGS5e4+3eMjLgp\nsO3thY6IEKJNNGIsOG33xm9iaSk3DED5pBQdDchklZNS586C9209fcrNxRMeDnzxBTB5skbnSUKI\nhqIEpAa8vImMAY8fc4no5s1XSenpU+46kqPjq6Rkbw80bare41cR0s6dgI8PMGsWN36bvj7vhyWE\n1FOCJ6CtW7cqBCESiWBoaAgnJyc4OjqqPTA+1Olo2FlZXJdd+ZbSgwdcy6h8S8nRUa0XYx48AD78\nkJuVdPt27hCEEPImBE9AHh4eiIqKwpgxY8AYQ1BQEBwcHPDo0SNMnDgRK1asUHtw6ib4dAwFBdzI\nnuWT0l9/Aaamlbvw2rev1a6Li7lutq1buSo3b2+gkcolJoQQUjXBE9DAgQMRHBwM/X/7cnJzczFy\n5EiEhITAyckJ9+/fV3tw6iZ4AlJGJgNiYytfV9LTq5yUOnZUOrf1tWvA++9zOev777kpsQkhRF0E\nL8NOS0tTGPdNT08PT58+RfPmzdG0Dq5r1Fu6uoCtLbd4eHDrGAOSkl4lo717uYHZXrwAevSQJ6Ri\n+55Ye9gO23fp4csvuZdTkQEhRFuonIA8PT3Rp08fuLm5gTGGEydOwMPDAy9fvoSdnR2fMTY8IhFg\nYcEt48a9Wp+eLi90yDz8BzLmbcGa4gT42tlA72xPIPPfllKPHlR1QAjReLWqgrt27RouX74MABgw\nYACcnZ15C4wPGtkFV0slJa+u9Xz+OTBz4kuI7vyl2H139y4gkVTuwmvTRujwCSFaSPBrQAUFBTh6\n9CgSEhJQ8u/YaSKRCP/973/VHhRftD0BPXjADZ3TtClXZt2hQxUbFhdzQx2UT0o3b3KtoopJqUMH\n6rcjhFRL8ATk6uqKVq1awcnJCbrlbrBctmyZ2oPii7YmIMa4SeHWrOHu6VmwQGktQs07efiwcrFD\nQYHivUo9ewJdu1IJHSFETvAEZG9vjzt37qg9gLqkjQkoIwPw8uJmKd2/n8sNavX0qeINtNHRQEoK\nd9Ns+aTk4AD8OyU7IaRhETwBffDBB1iwYAG6d++u9iDqirYloD//5Crb3N2BzZuBJk3q6MA5OZVv\nov37b64MvOzm2bLERFOmElLvCZ6AbG1tERcXBysrKzT595NQK6dj0IIEJJNxCcffnxvNYOxYoSMC\nUFiofHDW1q0rX1cyN6frSoTUI4InoISEBKXrLS0t1RgOv7QhAT15Akybxg3GvX8/N9C2xiotBeLi\nKl9XYkz54Ky1vnBFCNEEgieg+kDTE9AffwAzZgBz53LFBlpZB8AYkJpaOSk9f84Nzlo+KXXrVof9\nioSQ1yVYAhowYAD+/PNP6OvrQ1ShW0UkEiE7O1vtQdUkMDAQQUFByM7OhpeXF4YOHYqXL19CKpXC\nz88Po6qY60dTExBjwKZN3Fw9e/cC77wjdEQ8yMqqXOwQHw906aKYlHr0AFq2VHjphaAghAYEoFFh\nIUqaNMEwb++a53MihKgNtYCUyMrKwvLly7Fjxw74+vrCwMAAtra2WpWAcnK4e3tSU4EjRzS8y03d\n8vOVD85qZiYvdLhQVITTu3djY7ku4FXW1nD196ckREgdqZcJaM6cOQgKCkLbtm3x119/ydeHhIRg\n8eLFkMlkeP/996scaXv58uWYNm0a0tLSkJGRgYKCApiYmGhNAoqJ4SrcBg0CAgKoNwoAN9RDucFZ\nV//8MzZkZlbabI2NDdZv3MgNVySRAG3bUuEDITwRLAEp63orH9SbdMFdvHgR+vr6mDFjhjwByWQy\ndO3aFWfOnIG5uTlcXFxw4MABREVF4caNG/j444/Rvn17+Pj4YNiwYRgyZAhWr16Nly9f4t69e2jW\nrBl+++03pTFrUgI6fhz44APg00+5kayJcn5SKfzOn6+83swMfi4u3A1SSUlcU1Is5pJRWVIq+1r2\nfYWuPUKIagQbDTs3N1ftBy0zcODAStV1kZGR6NSpk7y6bsqUKQgMDISPjw+mT58OAAgICEBYWBiy\ns7MRFxeHDRs2AAB2796NNm3aVJkwNYFMBvj6Ar/8Apw8CfTuLXREmq2kimahzMGBy+Jl8vKA5ORX\nCSkpCYiMBI4efbVOV7dyUir/VSymZighdahWdVa3bt3ChQsXIBKJMHDgQPTo0UPtAaWkpEAikcgf\ni8ViREREKGzj7e0Nb2/vSq+dOXNmjfv38/OTfy+VSiGVSl871trKyeFuLM3JAaKiuF4jUr1h3t5Y\nFR+PjfHx8nUrra0xfOFCxQ2bN+cKGrp0Ub4jxrhCiLJkVPY1NPTV96mp3I21ylpPZV9NTblERkg9\nFh4ejvDwcN6Po3IC8vf3x/bt2zF+/HgwxjBt2jTMnTtXaSJ4E3y3XsonoLqUlASMGQO4uADHjnHz\nzZGalRUarPnmG+gWFEDWtCmGL1xY+wIEkYhLLkZGXKWdMjIZNzRR+QSVmAhcvvzq+8xMrkiiqlaU\nRMIdQ4Nb4YTUpOI/52vXruXlOConoB07diAiIgItWrQAAPj4+KBv375qT0Dm5uZISkqSP05KSoJY\ny0vDoqIANzdg8WJuXjn6bKqdQaNG1U3Fm64ul1zMzIA+fZRvU1jIdfWVT1K3b3P9qWXrZLLqW1ES\nCY2rRwhq2QWnU+5Odh2e7mp3dnbGgwcPkJCQADMzMxw8eBAHDhzg5Vh14dgx4MMPuSF13NyEjoa8\nsSZNAGtrbqlKdnblVlR4+Kvvk5O5qTGqKpaQSLgkqJV3IhOiOpV/w2fPno0+ffrIu+COHz+OOXPm\nvNHBp06divPnzyM9PR0SiQTr1q3D7Nmz8e2338LV1RUymQxeXl6wtbV9o+OU5+fnVyfXfhjjJoz7\n5hvg9GmgVy9eD0c0ScuW3CgP3bopf760FEhLe1UsUZaooqJerXv2jLveVF1Vn4kJNacJr/i+FlSr\n+4CuX7+OP//8EwBXwdazZ0/eAuNDXZVhFxcD8+cD168DJ040sJtLiXoUF3NFEeWr+ioWUOTlvUpI\nVSUqAwOhz4TUA4KVYZenq6srLxLgqwtO2+XmAhMnckUGFy9yPS2E1JqeHjdbbZXT3gJ4+bJycrpy\nBTh06NW6Jk2qb0WZmwONG9fdeRFSjsotoIpVcMePH+elCo5PfLeAnj8HRo3i5nL78UfqwicCY4yb\n0bC6VtTjx1xXXnWtKFNTGsm8gRN8KB4HBwdcvXpVXgX38uVL9O3bV2EIHU0nEong6+vLyzWgR48A\nV1dg/Hhg40bqmidaQibjklD5pFQxUb14wbWUqqvqMzSkX/p6qOwa0Nq1a4VPQJGRkWj2b/lofn4+\nevfurXUJiI838c4dYMQIYPlyYNEite+eEGHl578qPVfWikpM5LarrqpPIgGaNhX2PMhrE/waEB9V\ncPXBn39yrZ6vvuJGOSCk3mnWjJtQsHNn5c8zxrWSKraiwsJerUtJ4VpJ1V2Pat+eRploYGpdBXfp\n0iX5UDwNvQruxAlgzhxuDh9XV7XtlpD6p7SUKy2v2Hoq//3z51wSqu4mXmNj6uoTgODXgOoDdV4D\n2r8fWLoU+P13GlCUELUoKuJaStUVTRQWVt+KkkiAf69TkzenMdeAZsyYAX9/fxgZGQEAMjMzsWzZ\nMvz8889qD4ov6sriO3cCq1dzN5ja26shMEKIanJyKienit83b159VZ+5OQ3GWEuCt4AcHR1x8+bN\nGtdpMnW8iT/+CGzYAJw5A3TtqqbACCHqwRjXlaes9VT2/dOnQJs21Vf1tWlDpeflCF6EwBhDRkYG\njI2NAQAZGRmQyWRqD0iT+ftzxQbh4dUPBUYIEYhIxCWPNm2qHv+qpIQbZaJ8UoqLA86efbVO2QSH\nFRMVTXD4xlROQMuWLUO/fv0wefJkMMZw+PBhrFq1is/YNMqWLcC2bcD589XfnE4I0XCNGnEJxMKi\n6m0qTnCYmAhcu8aNLlx+gkNlrSea4FBltSpCuHv3Ls6ePQuRSIR33nkHdnZ2fMamdq9ThMAYsH49\nV3QQFsZ1HxNCGjjGuLmhqiqWSEriCiqMjbV6gkONKUKoD2rbj8kYV2zw++/cNR9TUx6DI4TUL2UT\nHFZX1ZeZyZWeV1fVpwETHApehFAf1PZNXLOGSz5hYdxwWYQQolbKJjismKjKJjisqqpPIuEq/3gk\naBECYwzJycmQSCRqD0BTrV/PdfeGh1PyIYTwRJUJDstGmSifnMLDXyWpihMcKktUrznB4YWgIIQG\nBLz++dVApRYQYwwODg64c+cOb4HUBVWz+KZNwK5d3M+4XTvewyKEkNdXfoLDqlpRaWm1nuDwQlAQ\nTi9ahI3x8RABwrWARCIRnJycEBkZid71/Lb/L78Efv6Zkg8hREvo6HDJxdQUcHZWvk3FCQ4TE4GY\nGCA09FXCqjDBYejFi9j4zz+8hq5ym+zq1avYu3cvOnToIJ+SQSQS4fbt27wFx4fqpuT+5hvgf//j\nko+ZWZ2HRggh/KjNBIf/JqlG4eEIBxDOY1gqFyEkJCRwL/i3iVb2MktLS14C40N1XXA//MB1vYWH\nA1p0SoQQwovVrq7YEBoKALx1wak81oSlpSWysrLw+++/48SJE3jx4oVWJZ/q/PwzN4lcWBglH0II\nAYBh3t5YxfOQLyonIH9/f0ybNg1paWl4+vQppk2bhgAeqyPqyqFDXLl1WBgNr0MIIWUGjRoFV39/\nrOFxrpkGNyV3+dMNDgZmzQL++APo3l24uAghRJMJPhgpAOiUGx1WR8tHir14EZgxg7vRlJIPIYTU\nvQY5JfeNG8CECdz4bv36CR0NIYQ0TA1uSu7/+z9f7N8vxU8/SeHuLnREhBCiuWgwUjUSiUSQSBjW\nrwdmzhQ6GkII0Q58XQPS7gs5r2H5cko+hBCiCRpcC6gBnS4hhKgFtYAIIYTUKypXwW3dulUhC4pE\nIhgaGsLJyQmOjo68BUgIIaR+UrkLzsPDA1FRURgzZgwYYwgKCoKDgwMePXqEiRMnYsWKFXzH+sao\nC44QQmpP8BlRBw4ciODgYOjr6wMAcnNzMXLkSISEhMDJyQn3799Xe3DqRgmIEEJqT/BrQGlpaWjc\nuLH8sZ6eHp4+fYrmzZujadOmag+MEEJI/abyNSBPT0/06dMHbm5uYIzhxIkT8PDwwMuXL2FnZ8dn\njGpV3XxAhBBCXim7EZUvtSrDvnbtGi5fvgwAGDBgAJyrmn1PQ1EXHCGE1J7g14AKCgpw9OhRJCQk\noKSkRB7Uf//7X7UHxRdKQIQQUnuCj4Y9btw4tGrVCk5OTnTNhxBCyBtTuQVkb2+PO3fu8B0Pr6gF\nRAghtSd4FVz//v1x+/ZttQdACCGkYVK5BWRra4u4uDhYWVmhSZMm3ItFIq1KStQCIoSQ2hO8CCEh\nIUHpektLSzWGwy9KQIQQUnuCJ6D6gBIQIYTUnmDXgAYMGAAA0NfXh4GBgcLSsmVLtQdECCGkYaAW\nECGEkGoJXgWnbLRrbRgBmxBCiGZSOQGFhoZWWnfq1Cm1BkMIIaThqHEkhO+//x7fffcd4uPj4eDg\nIF+fk5Mjvz6kTWgwUkIIUY3gg5G+ePECmZmZ8PHxwebNm+X9gAYGBmjdujVvgfGBrgERQkjtaUQZ\ndmZmJh48eICCggL5ukGDBqk9KL5QAiKEkNoTfDDS7du3IyAgAMnJyXB0dMTVq1fRr18/nD17Vu1B\nEUIIqf9ULkLw9/dHZGQkOnTogHPnziE6OhqGhoZ8xkYIIaQeUzkBNW3aFM2aNQPAzQ1kY2ODv//+\nm7fACCGE1G8qd8FJJBJkZmbCzc0NQ4cOhZGRkVaNA0cIIUSzvNZICOHh4cjOzsbw4cPRuHFjPuLi\nBRUhEEJI7QleBVdaWop9+/bh4cOH+O9//4vExEQ8efIEvXv3VntQfKEERAghtSd4Apo3bx50dHRw\n9uxZxMTEICMjA8OGDUNUVJTag+ILJSBCCKk9wcuwIyIiEB0djZ49ewIAjI2NUVxcrPaACCGENAwq\nV8E1btwYMplM/jgtLQ06Oiq/nBBCCFGgcgZZuHAh3N3d8ezZM6xcuRIDBgzAJ598wmdshBBC6rFa\nVcHdv38fZ8+eBWMMQ4YMga2tLZ+xqR1dAyKEkNoTrAhh4cKFVQYhEokQEBCg9qD4QgmIEEJqT7Ai\nBCcnJ/nBfX19sW7dOgAAYwwikUjtARFCCGkYatUF17NnT0RHR/MZD6+oBUQIIbUneBm2JgkMDERQ\nUBCys7Ph5eUFPT09rFmzBvb29pgyZQoGDx4sdIiEEEJqoJV11OPGjcO2bdvwww8/4ODBg9DR0YGB\ngQEKCwshFouFDo8QQogKakxA+vr6MDAwgIGBAf766y/59wYGBmjZsuUbHXzOnDkwNTVVmOobAEJC\nQmBjY4POnTtj8+bNVb5+w4YNWLBgAQYOHIhTp05h06ZN8PX1faOYXgefU9byhWLmn7bFC2hfzNoW\nL6CdMfOlxgSUm5uLnJwc5OTkoKSkRP59Tk4OsrOz3+jgs2fPRkhIiMI6mUyGBQsWICQkBPfu3cOB\nAwdw//597NmzB0uWLEFqaioYY1ixYgVGjBgBR0dHeTFEq1atUFhY+EYxvQ5t/IWimPmnbfEC2hez\ntsULaGfMfBH0GtDAgQORkJCgsC4yMhKdOnWST/UwZcoUBAYGwsfHB9OnTwcABAQEICwsDNnZ2YiL\ni0Pbtm1x+vRpZGVlKZSNE0II0VwaV4SQkpICiUQifywWixEREaGwjbe3N7y9vRXWubu710l8hBBC\n1OO15gNSp4SEBIwZMwZ//fUXAODo0aMICQnB9u3bAQB79+5FREQEvvnmmzc+Ft23RAghr6dBlGGb\nm5sjKSlJ/jgpKUltlW10DxAhhGgOjSvDdnZ2xoMHD5CQkICioiIcPHgQY8eOFTosQgghaiZoApo6\ndSr69++P2NhYSCQS7Ny5E40aNcK3334LV1dX2NnZ4b333tO6QU8JIYSogBGlMjMz2YQJE5iNjQ2z\ntbVlV69eZenp6ezdd99lnTt3ZkOHDmWZmZny7T/99FPWqVMn1rVrV3b69Gn5+qioKGZvb886derE\nvGFEtqYAAA3QSURBVL29eY35008/ZXZ2dsze3p5NnTqVFRQUaFzMs2fPZm3btmX29vbydeqMsaCg\ngE2ePJl16tSJ9enThyUkJKg93uXLlzMbGxvWvXt35u7uzrKysjQm3qpiLvPFF18wkUjE0tPTNSbm\nquINCAhgNjY2rFu3buw///mPxsRbVcwRERHMxcWFOTo6MmdnZxYZGalRMScmJjKpVMrs7OxYt27d\nmL+/P2NM2L8/SkBVmDFjBvvpp58YY4wVFxezrKws9vHHH7PNmzczxhjbtGkTW7FiBWOMsbt377Ie\nPXqwoqIi9vDhQ2Ztbc1KS0sZY4y5uLiwiIgIxhhjI0aMYMHBwbzE+/DhQ2ZlZcUKCgoYY4xNnjyZ\n7dq1S+NivnDhArtx44bCH646Y/zf//7H5s+fzxhj7Ndff2Xvvfee2uMNDQ1lMpmMMcbYihUrNCre\nqmJmjPsAcnV1ZZaWlvIEpAkxK4v37Nmz7N1332VFRUWMMcaePXumMfFWFfPgwYNZSEgIY4yxU6dO\nMalUqlExP378mEVHRzPGGMvJyWFdunRh9+7dE/TvjxKQEllZWczKyqrS+q5du7InT54wxrgfZteu\nXRlj3H8JmzZtkm/n6urKrly5wlJTU5mNjY18/YEDB9iHH37IS8zp6emsS5cuLCMjgxUXF7PRo0ez\n0NBQjYz54cOHCn+46ozR1dWVXb16lTHG/eNgYmKi9njLO3bsGPP09NSoeKuKeeLEiezWrVsKCUhT\nYq4Y76RJk1hYWFil7TQlXmUxT5kyhR08eJAxxtj+/fs18veivHHjxrE//vhD0L8/jStC0AQPHz5E\nmzZtMHv2bPTq1Qtz587Fy5cv8fTpU5iamgIATE1N8fTpUwBAamqqQqWeWCxGSkpKpfXm5uZISUnh\nJWZjY2MsW7YMFhYWMDMzQ6tWrTB06FCNjrmMOmMsfx9Zo0aNYGhoiIyMDN5i//nnnzFy5EiNjzcw\nMBBisRjdu3dXWK+pMT948AAXLlxA3759IZVKERUVpdHxAsCmTZvkf4Mff/wxPvvsM42NOSEhAdHR\n0ejTp4+gf3+UgJQoKSnBjRs38NFHH+HGjRto0aIFNm3apLCNSCTSqPuK4uPj8fXXXyMhIQGpqanI\nzc3F3r17FbbRtJiV0YYYy2zcuBGNGzeGh4eH0KFUKy8vD59++inWrl0rX8c0/JaEkpISZGZm4urV\nq/j8888xefJkoUOqkZeXFwICApCYmIivvvoKc+bMETokpXJzczFhwgT4+/vDwMBA4bm6/vujBKSE\nWCyGWCyGi4sLAGDixIm4ceMG2rVrhydPngAAHj9+jLZt2wKofO9ScnIyxGIxzM3NkZycrLDe3Nyc\nl5ijoqLQv39/tG7dGo0aNcL48eNx5coVjY65jKmp6RvHWPYfmbm5ORITEwFwH2IvXryAsbGx2mPe\ntWsXTp06hX379snXaWq88fHxSEhIQI8ePWBlZYXk5GQ4OTnh6dOnGhuzWCzG+PHjAQAuLi7Q0dHB\n8+fPNTZegBtGrGxElokTJyIyMlJ+fE2Jubi4GBMmTMD06dPh5uYGQNi/P0pASrRr1w4SiQSxsbEA\ngDNnzqBbt24YM2YMdu/eDQDYvXu3/Ac4duxY/PrrrygqKsLDhw/x4MED9O7dG+3atUPLli0REREB\nxhj27Nkjf4262djY4OrVq8jPzwdjDGfOnIGdnZ1Gx1xm7NixbxzjuHHjKu3ryJEjGDJkiNrjDQkJ\nweeff47AwEA0bdpU4Tw0MV4HBwc8ffoUDx8+xMOHDyEWi3Hjxg2YmppqbMxubm44e/YsACA2NhZF\nRUUwMTHR2HgBoFOnTjh//jwA4OzZs+jSpYv8+JoQM2MMXl5esLOzw+LFi+XrBf37e/NLWfXTzZs3\nmbOzs0KpbXp6OhsyZIjScsWNGzcya2tr1rVrV3klDGOvyhWtra3ZwoULeY158+bN8jLsGTNmsKKi\nIo2LecqUKax9+/ZMT0+PicVi9vPPP6s1xoKCAjZp0iR5GejDhw/VGu9PP/3EOnXqxCwsLJijoyNz\ndHSUV/1oQrzlY27cuLH8PS7PyspKoQxb6JiVxVtUVMSmTZvG7O3tWa9evdi5c+c0Jt7yMZf/Pb52\n7Rrr3bs369GjB+vbty+7ceOGRsV88eJFJhKJWI8ePeS/u8HBwYL+/Qk+FhwhhJCGibrgCCGECIIS\nECGEEEFQAiKEECIISkCEEEIEQQmIEEKIICgBEUIIEQQlIKLxdHR0sHz5cvnjL774QmFYGWX8/Pyw\ndevWWh1HKpXi+vXrKm+/e/duPH78uFbHEML58+dx5coV3vZfWFiIwYMHgzGGhIQEODg4yJ/bvn07\nnJ2dkZWVhaVLl+LixYu8xUG0DyUgovEaN26M3377Denp6QCg0lhVrzOeVW3Hwdq1axdSU1NrfRw+\nyGSyKp87d+4cLl++XKv9lZSUqLztvn37MHr06Erv3Z49e/Dtt98iNDQUrVq1wvz58/H555/XKg5S\nv1ECIhpPT08PH3zwAb766qtKz6WlpWHixIno3bs3evfurfBBe+vWLfTv3x9dunTBjh07AADh4eEY\nM2aMfJsFCxbIhw4p76OPPoKLiwvs7e3h5+dX6fkjR44gKioKnp6e6NWrFwoKCnD9+nVIpVI4Oztj\n+PDh8vG1pFIpli5dChcXF9ja2uLatWtwd3dHly5dsGbNGgDc6MQ2NjaYNm0a7OzsMGnSJOTn5wNA\ntftdsmQJXFxc4O/vj5MnT6Jv377o1asXhg4dimfPniEhIQE//vgjvvrqK/Tq1QuXLl3CrFmzcPTo\nUfm56Ovry9+bgQMHYty4cbC3t0dpaSk+/vhj9O7dGz169MC2bduU/nwOHDggH4qlzKFDh7B582b8\n8ccf8rHAOnfujISEBGRlZSndD2mA3nh8B0J4pq+vz7Kzs5mlpSV78eIF++KLL5ifnx9jjLGpU6ey\nS5cuMcYYe/ToEbO1tWWMMebr68t69OjBCgoK2PPnz5lEImGpqans3LlzbPTo0fJ9L1iwgO3evZsx\nxphUKmXXr19njDGWkZHBGGOspKSESaVSdvv27Upxld++qKiI9evXjz1//pwxxk3GNWfOHPl2Pj4+\njDHG/P3/v727C2X/i+MA/p6nxTzGFvKw3BDzsIVSXChPF8iNKSREkZQiKxdy4ZKL/0Ik5Uqplaw8\nJIULQoqSyVKTlVkeYkOMdv4X/jtt/+2P381v/H+f19X3e3a+n+9n34tz+n7P6Zy/WExMDLu8vGQv\nLy8sLi6O3d7eMoPBwAQCAdva2mKMMdbc3MyGhobY6+vrh3E7Ojp4Ps5LqExOTrLu7m7GGGMDAwNs\neHiY/9bY2Mg0Go3L82WMsbW1NSYSifgulhMTE2xwcJAx9r7ESnZ2ttvSKm9vbyw6OpqfGwwGFhwc\nzCQSCbu4uHB7Zg0NDWxxcdGtnPyZ/LzdARLyFSEhIWhoaIBarUZgYCAvX11dxfHxMT+3Wq14fHyE\nQCBAVVUVhEIhhEIhCgsLsbu7i/Dw8C/db3Z2FpOTk3h7e4PJZIJOp3MZ23Bg/6xkdXJygqOjIxQV\nFQF4/yQWGxvL61VWVgIAZDIZZDIZ338lKSkJRqMRoaGhiI+PR15eHgCgvr4earUaZWVlH8atqanh\nx0ajEUqlEpeXl7DZbEhKSnLL8zO5ublITEwEAKysrODw8BAajQYAYLFYcHp6CqlUyutfX1+7Lekv\nkUgQGRmJ2dlZl0UvASA2NhZnZ2dfyoX8/1EHRH6Mrq4uKBQKNDU18TLGGHZ2dhAQEPDp9T4+PvDz\n84Pdbudljs9czgwGA4aHh7G3t4ewsDA0NTXh+fnZY0zHuAdjDGlpaf851iIUCnkOjmPHuWO8xXkM\nhTEGgUDwaVyRSMSPOzs70dPTg/LycmxsbHj8dAjA5RnY7XbYbDaP8QBgZGQExcXFHuM45+osKCgI\nCwsLKCgogEQicdkvyfG/CAFoDIj8IBEREVAqlZiamuKNWElJCdRqNa9zcHAA4L2hm5+fx8vLC25u\nbrC+vo6cnBwkJCRAp9PBZrPh7u6OL/nvzGKxQCQSITQ0FGazGUtLSx4bzZCQEFgsFgBAcnIyrq6u\nsL29DeB93xWdTvdL/+/8/JxfPzMzg4KCgk/jOjf+FouFvx1NT0+75Gm1Wvm5VCrls/20Wi1eX189\n5lNaWoqxsTHeQer1ejw9PbnUiYqKwsPDg9u1YrEYy8vL6Ovrw8rKCi83mUwub1Dkz0YdEPn2nBv/\n7u5uXF9f83O1Wo29vT1kZmYiLS2ND5QLBAJkZGSgsLAQeXl56O/v5/s8KZVKyGQy1NTUQKFQuN0v\nMzMTcrkcKSkpqKurQ35+vse8Ghsb0dbWBoVCAbvdDo1GA5VKhaysLMjlco9Tnz+aaZecnIzR0VGk\npqbi/v4e7e3t8Pf3/zCuc6yBgQFUV1cjOzsbYrGY/1ZRUYG5uTnI5XJsbm6itbUVGxsbyMrKwvb2\nNp+E8O94LS0tSE1NhUKhQHp6Otrb291mx/n6+kImk+Hk5MQthlQqhVarRXNzM99Se39/n39mJIS2\nYyDkGzg7O0NFRQUODw+9ncovm56ehtlshkql+rCeXq9HT08PtFrtb8qMfHf0BkTIN/FTx0Zqa2ux\nsLDw6USH8fFx9Pb2/qasyE9Ab0CEEEK8gt6ACCGEeAV1QIQQQryCOiBCCCFeQR0QIYQQr6AOiBBC\niFdQB0QIIcQr/gaTbpjdLeR2UgAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1092aac90>"
       ]
      }
     ],
     "prompt_number": 245
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "By reading the above plot carefully, it's clear that a heating-cooling balance is reached in this nebula at a temperature around 7500 K."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### (f) \n",
      "For your estimated temperature, estimate the ratio of the sum of the fluxes of the [OIII] lines to H$\\beta$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "From the slides, \n",
      "$$ 4 \\pi j_{H \\beta} / n_e n_p \\approx 1.7 \\times 10^{-25} \\textrm{erg cm}^3 s^{-1}$$\n",
      "(interpolating between T=5000 and T=10,000 K).\n",
      "\n",
      "Taking a straight ratio between this value and the energy loss from [O III] lines at T=7500 K:\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "flux_ratio = oIII_cooling(7500)/(1.7e-25 * u.erg * u.cm**3 /u.s)\n",
      "print \"The flux ratio between the [OIII] lines and H-beta\"\n",
      "print \"is {0:.2f}\".format(flux_ratio)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The flux ratio between the [OIII] lines and H-beta\n",
        "is 3.19\n"
       ]
      }
     ],
     "prompt_number": 246
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### (g)\n",
      "Suppose the O abundance were reduced by a factor of two (lower than solar metallicity) in this nebula. Assume the same ratio of O$^+$ to O$^{++}$. Qualitatively, would the [O III]/H$\\beta$ ratio decrease by a factor of two? Why or why not?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Here is how the flux lines depend on oxygen abundance and temperature:\n",
      "\n",
      "$$ F_{[OIII]} \\propto n_{O^{++}} \\frac{1}{T^{1/2}} e^{-h\\nu/ kT} $$\n",
      "$$F_{H \\beta} \\propto 1/T $$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "I think the ratio would *not* decrease by a factor of two, because even though the oxygen abundance decreases by that amount, the corresponding nebular temperature has to increase to make up for the less efficient cooling, and this drives the [O III] flux up higher than simply \"half\", while the H$\\beta$ flux would go down, per the above relations."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 246
    }
   ],
   "metadata": {}
  }
 ]
}