{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adiabatic flame temperatures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Flame temperatures depend on\n",
    "* stoichiometry\n",
    "* products\n",
    "* fuel\n",
    "* oxidizer\n",
    "* enthalpy\n",
    "\n",
    "These properties are affected by\n",
    "* heat losses (radiation)\n",
    "* differential diffusion\n",
    "* flame strain, mixing fields\n",
    "* residence time\n",
    "* soot formation\n",
    "\n",
    "**For convenience here**\n",
    "* use products of complete combustion, \n",
    "* at constant pressure (1 atm), \n",
    "* with reactants at 298.15 K (25 C)\n",
    "* adiabatic\n",
    "\n",
    "**Two kinds of energy: internal, and sensible**\n",
    "* Internal is energy stored in bonds\n",
    "* Sensible is thermal energy\n",
    "* Combustion reactants have more energy stored in chemical bonds (high internal energy) than combustion products (low internal energy). \n",
    "    * The difference is sensible energy that results in high flame temperatures.\n",
    "* That is, **internal energy is released as sensible energy**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Enthalpy, energy, C$_p$\n",
    "\n",
    "* Considering **ideal gases**\n",
    "* Total enthalpy is the sum of species enthalpies times the amount of the species.\n",
    "$$H = \\sum_k n_kh_k$$\n",
    "$$h_k = h_{f,k}^o + \\int_{T_r}^T c_{p,k}dT$$\n",
    "\n",
    "     * $n_k$ is moles of species $k$, \n",
    "     * $H$ is total enthalpy (J), \n",
    "     * $h_k$ is enthalpy per mole (J/kmol).\n",
    "     * $c_{p,k}$ is heat capacity of species $k$ (J/kmol*K)\n",
    "     \n",
    "* If we divide the first equation by the total number of moles we have\n",
    "$$ h = \\sum_k x_kh_k$$\n",
    "* Insert the second equation into this:\n",
    "$$ h = \\sum_k h_{f,k}^o + \\sum_k\\int_{T_r}^Tc_{p,k}dT$$\n",
    "$$ h = h_{f}^o + \\int_{T_r}^Tc_pdT$$\n",
    "* Here, $h_{f}^o = \\sum_kh_{f,k}^o$ and $c_p = \\sum_kc_{p,k}$\n",
    "* All of this can be done on a mass basis too, just replace $x_k$ with $y_k$ and $n$ for $m$. Then $h_k$ is J/kg instead of J/kmol, and $c_p$ is J/kg*K instead of J/kmol*K. \n",
    "    * To convert $h_k$ from per mole to per mass, divide by the molecular weight. Same for $c_p$.\n",
    "* $h_{f,k}^o$ is the heat of formation of species $k$ at 1 atm and 298.15 K. This is the enthalpy of formation from elements in their standard states.\n",
    "\n",
    "### C$_p$\n",
    "* Heat capacities are **temperature dependent**\n",
    "    * Usually have a polynomial curve fit to represent $c_p(T)$.\n",
    "    * Hence, enthalpy, entropy, etc. are temperature dependent.\n",
    "    * $c_p$ are normalized by $R$, so the units on $c_p$ are whatever the units on the chosen $R$ are.\n",
    "    * $c_p$ increases with temperature.\n",
    "\n",
    "### Heat of reaction\n",
    "* $\\Delta H_{rxn} = H_{prod}-H_{react}$\n",
    "    * **reactants and products at the SAME temperature**\n",
    "* $\\Delta H_{rxn} < 0$ since combustion reactions are exothermic.\n",
    "* $\\Delta H_{comb} = -\\Delta H_{rxn}$, so heats of combustion are positive.\n",
    "* Usually reported per mole of fuel, or per kg of fuel (be careful).\n",
    "\n",
    "### Heating values: LHV, HHV\n",
    "* As a basis for the amount of heat that can be extracted from combustion, the higher heating value (HHV) and lower heating value LHV are used.\n",
    "* The HHV treats product $H_2O$ as a liquid.\n",
    "* The LHV treats product $H_2O$ as a vapor.\n",
    "* These use products of complete combustion at 298.15 K and 1 atm.\n",
    "* **See Turns Table B1.**\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## T$_{ad}$ for CH$_4$\n",
    "\n",
    "$$CH_4 + 2O_2 + 7.52N_2 \\rightarrow CO_2 + 2H_2O + 7.52 N_2$$\n",
    "\n",
    "$H_R(298.15\\, K) = H_P(T_{ad})$\n",
    "* This is one equation in one unknown.\n",
    "* The reactant and product compositions are known. The reactant temperature is known, so the reactant enthalpy is known.\n",
    "* Solve $T_{ad}$.\n",
    "\n",
    "$$H_R = H_P(T_{ad})$$\n",
    "$$H_P(T_{ad}) = 1\\mbox{ mole}\\cdot h_{CO2}(T_{ad}) + 2\\cdot h_{H2O}(T_{ad}) + 7.52\\cdot h_{N2}(T_{ad})$$\n",
    "\n",
    "* Turns Table A.13 gives $h_{k}(T)$ as polynomial expressions:\n",
    "$$\\frac{h_{k}(T)}{RT} = a_1 + \\frac{a_2}{2}T + \\frac{a_3}{3}T^2 + \\frac{a_4}{4}T^3 + \\frac{a_5}{5}T^4 + \\frac{a_6}{T}$$\n",
    "    * where a different set of coefficients $a$ are used for each species $k$\n",
    "* We insert these expressions for $h_k(T)$ into the equation for $H_P(T_{ad})$ and solve for $T_{ad}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set some data and make a ```thermo``` class for convenience"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "import yaml\n",
    "import numpy as np\n",
    "from scipy.optimize import fsolve\n",
    "\n",
    "########################################################################\n",
    "# Thermo data from Cantera's gri30.cti file\n",
    "\n",
    "thermoData = '''\n",
    "Tlow:  200\n",
    "Tmid:  1000\n",
    "Thigh: 3500\n",
    "\n",
    "CH4:\n",
    "    a_lo: [  5.149876130E+00,  -1.367097880E-02, 4.918005990E-05,  -4.847430260E-08,   1.666939560E-11,-1.024664760E+04,  -4.641303760E+00]\n",
    "    a_hi: [  7.485149500E-02,   1.339094670E-02, -5.732858090E-06,   1.222925350E-09,  -1.018152300E-13,-9.468344590E+03,   1.843731800E+01]\n",
    "O2:\n",
    "    a_lo: [  3.782456360E+00,  -2.996734160E-03, 9.847302010E-06,  -9.681295090E-09,   3.243728370E-12,-1.063943560E+03,   3.657675730E+00]\n",
    "    a_hi: [  3.282537840E+00,   1.483087540E-03, -7.579666690E-07,   2.094705550E-10,  -2.167177940E-14,-1.088457720E+03,   5.453231290E+00]\n",
    "N2:\n",
    "    a_lo: [  3.298677000E+00,   1.408240400E-03, -3.963222000E-06,   5.641515000E-09,  -2.444854000E-12,-1.020899900E+03,   3.950372000E+00]\n",
    "    a_hi: [  2.926640000E+00,   1.487976800E-03, -5.684760000E-07,   1.009703800E-10,  -6.753351000E-15,-9.227977000E+02,   5.980528000E+00]\n",
    "CO2:\n",
    "    a_lo: [  2.356773520E+00,   8.984596770E-03, -7.123562690E-06,   2.459190220E-09,  -1.436995480E-13,-4.837196970E+04,   9.901052220E+00]\n",
    "    a_hi: [  3.857460290E+00,   4.414370260E-03, -2.214814040E-06,   5.234901880E-10,  -4.720841640E-14,-4.875916600E+04,   2.271638060E+00]\n",
    "H2O:\n",
    "    a_lo: [  4.198640560E+00,  -2.036434100E-03, 6.520402110E-06,  -5.487970620E-09,   1.771978170E-12,-3.029372670E+04,  -8.490322080E-01]\n",
    "    a_hi: [  3.033992490E+00,   2.176918040E-03, -1.640725180E-07,  -9.704198700E-11,   1.682009920E-14,-3.000429710E+04,   4.966770100E+00]\n",
    "'''\n",
    "\n",
    "########################################################################\n",
    "\n",
    "class thermo():\n",
    "\n",
    "    def __init__(self, species) :\n",
    "        \"\"\"\n",
    "        species: input, string name of species in thermoData.yaml\n",
    "        M: input (species molecular weight, kg/kmol)\n",
    "        Call as: h2o = thermo(\"H2O\"), then get enthalpy as h2o.h_mole(298.15), etc.\n",
    "        \"\"\"\n",
    "\n",
    "        s = self\n",
    "\n",
    "        #tData = open(\"thermoData.yaml\"); tData = yaml.load(tData)\n",
    "        tData = yaml.load(thermoData)\n",
    "\n",
    "        s.a_lo = tData[species][\"a_lo\"]\n",
    "        s.a_hi = tData[species][\"a_hi\"]\n",
    "\n",
    "        s.T_lo = tData[\"Tlow\"]\n",
    "        s.T_mid = tData[\"Tmid\"]\n",
    "        s.T_hi = tData[\"Thigh\"]\n",
    "\n",
    "    #--------------------------------------------------------\n",
    "\n",
    "    def h_mole(self,T) :\n",
    "        \"\"\"\n",
    "        T: input, (K)\n",
    "        return enthalpy in units of J/kmol\n",
    "        \"\"\"\n",
    "        s = self\n",
    "        \n",
    "        if T<=s.T_mid and T>=s.T_lo :\n",
    "            a = s.a_lo\n",
    "        elif T>s.T_mid and T<=s.T_hi :\n",
    "            a = s.a_hi\n",
    "        else :\n",
    "            raise ValueError(\"ERROR: temperature is out of range\")\n",
    "\n",
    "        hRT = a[0] + a[1]/2.0*T + a[2]/3.0*T*T + a[3]/4.0*T**3.0 + a[4]/5.0*T**4.0 + a[5]/T\n",
    "\n",
    "        return hRT * 8314.46 * T "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate T$_{ad}$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Adiabatic flame temperature = 2325.59812975 K\n"
     ]
    }
   ],
   "source": [
    "ch4 = thermo(\"CH4\")\n",
    "o2  = thermo(\"O2\")\n",
    "n2  = thermo(\"N2\")\n",
    "co2 = thermo(\"CO2\")\n",
    "h2o = thermo(\"H2O\")\n",
    "\n",
    "Tr = 298.15\n",
    "Hr = ch4.h_mole(Tr)  + 2*o2.h_mole(Tr)   + 7.52*n2.h_mole(Tr)\n",
    "\n",
    "def F(Tad):\n",
    "    Hp = co2.h_mole(Tad) + 2*h2o.h_mole(Tad) + 7.52*n2.h_mole(Tad)\n",
    "    return Hr - Hp\n",
    "\n",
    "Tad = fsolve(F, Tad_guess)[0]\n",
    "\n",
    "print(\"Adiabatic flame temperature =\", Tad, \"K\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Do the same thing using Cantera instead"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tad = 2325.5981297600447 K\n"
     ]
    }
   ],
   "source": [
    "import cantera as ct\n",
    "gas = ct.Solution(\"gri30.yaml\")\n",
    "P = 101325   # Pa\n",
    "\n",
    "#----------- Stoichiometric reactants\n",
    "Tr = 298.15\n",
    "gas.TPX = Tr, P, \"CH4:1, O2:2, N2:7.52\"\n",
    "hR = gas.enthalpy_mass\n",
    "\n",
    "#----------- Stoichiometric products\n",
    "gas.HPX = hR, P, \"CO2:1, H2O:2, N2:7.52\"       # note, enthalpy is reactant enthalpy\n",
    "yst = gas.Y\n",
    "print(\"Tad =\", gas.T, \"K\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute T$_{ad}$ over all $\\xi$ using Cantera\n",
    "\n",
    "* Species mass fractions are peicewise linear.\n",
    "    * That is, the lean compositions are linear between $\\xi=0$ and stoichiometric, and the rich compositions are linear between $\\xi=\\xi_{st}$ and $\\xi=1$.\n",
    "    * So, find the mass fractions at $\\xi=0$, $\\xi=\\xi_{st}$, and $\\xi=1$.\n",
    "* Lean:\n",
    "    * $(y-y_0)/(\\xi-0) = (y_{st} - y_0)/(\\xi_{st}-0)$\n",
    "    * $\\rightarrow y = y_0 + (y_{st} - y_0)\\xi/\\xi_{st}.$\n",
    "* Rich:\n",
    "    * $(y-y_{st})/(\\xi-\\xi_{st}) = (y_1 - y_{st})/(1-\\xi_{st})$\n",
    "    * $\\rightarrow y = y_{st} + (y_1 - y_{st})(\\xi-\\xi_{st})/(1-\\xi_{st})$\n",
    "    \n",
    "* Enthalpy is a conserved scalar, so \n",
    "    * $h = h_0\\cdot(1-\\xi) + h_1\\cdot(\\xi)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAESCAYAAACy36FdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4VNXWwOHfSoF0EkhICL2GLiWUIF1RUbx4ERvoxQofFhQVFbtexY6ooKJeRRGUK1iuCAqi9IAG6b23EBJ6QkuZ/f1xJmGIkAIzOZNkvc8zT5I5+5xZM2JW9t7r7C3GGJRSSim7+dgdgFJKKQWakJRSSnkJTUhKKaW8giYkpZRSXkETklJKKa+gCUkppZRX0ISklFLKK2hCUkop5RU0ISmllPIKfnYH4E0iIyNNnTp17A5DKaVKlWXLlh0wxkRd7HU0IbmoU6cOSUlJdoehlFKliojsdMd1dMhOKaWUV9CEpJRSyitoQlJKKeUVNCEppZTyCpqQlFJKeQVNSEoppbyCJiSllFJeQROSUkopr6AJSSmllFfQhKSUUsoraEJSSinlFTQhKaXKje7du3P//ffbHcZZbr/9dvr06VNmXudi6OKq6u92/wmJ78HVb0HIRS/gq5QqwDvvvIMxpsy8zsXQhKTOlnkCvhkEx/ZChRC47n27I1LK62VmZlKhQoULOrdSpUpujubCXudi3oO76JCdOtvcV6xkVK0VrPwKUjfYHZFSbpWdnc2DDz5IREQEERERjBgxAofDkXe8Tp06vPnmm2edk3+or06dOjz//PPceeedhIeHM3DgQHbs2IGIMG3aNHr16kVQUBBNmzZl9uzZBcaTfyite/fu3HvvvTz55JNERkZStWpVHn300bNizO/gwYPccsst1KhRg8DAQJo1a8Znn31W6OsMHTqURx99lKioKC699NKCP7gSUGIJSUSeFxGT75HiclycbZJF5KSIzBWRZvmuESEiE0XkqPMxUUTC87VpISLznNfYKyLPioiU1Pss1ZKXQ+JYaDMIbv0W/IPh95fsjkopt5o0aRIOh4PExETGjx/PRx99xJgxY4p9ndGjR9O4cWOSkpIYNWpU3vNPPfUUw4YNY+XKlbRr146bb76ZjIyMYsfo5+fH4sWLGTt2LGPGjGHKlCnnbX/q1CnatGnD9OnTWbt2LQ8++CBDhgxhzpw5Bb7Ol19+iTGGBQsW8MUXXxQrRk8o6SG7jUB3l59zXL5/DHgEuN3Z7llgtojEGWPSnW0mA7WAq5w/fwJMBK4FEJEwYDYwH2gHNAY+A44Db7n93ZQlOVnwwwMQXBV6vQiB4dDpAZg7CvYsgxpt7Y5QebkXflzLuuRjJfqaTWPDeO7aZoU3dFGtWjXeffddRITGjRuzadMmRo8ezcMPP1ys63Tr1o3HHnss7+cdO3YAMHz4cK699loARo0axRdffMGKFSvo3Llzka/dtGlTXnzxRQAaNWrExx9/zJw5c7jlllvO2b569eqMGDEi7+fBgwfz22+/8dVXX3HZZZed93Xq1q3LW295z6/Gkh6yyzbGpLg80sDqHQEPAa8aY6YZY9YAg4BQYICzTROsRDTYGJNojEkEhgB9RCTOef2BQBAwyBizxhgzFXgNeFh7SYVY/C7sXw3XvGUlI4CEeyEoEua8YG9sSrlRx44dcf11kJCQwN69ezl2rHjJND4+/pzPt2zZMu/72NhYAFJTU4t1bddr5F6noGvk5OTw8ssv07JlS6pUqUJISAjffvstu3btKvB12rb1rj80S7qHVE9EkoHTwFLgSWPMNqAuEAPMym1ojDkpIvOBTsB4IAHIABa7XG8RVu+nE1avKgFYYIw56dLmF+DfQB1gu2feVil3YAvMfQ2a9oUmLmWhFUOhyyPwy0jYNhfqdbcpQFUaFLen4q18fHz+Vo2WlZX1t3bBwcHnPN/f3z/v+9zEV9D8T2HXyL1OQdd48803eeutt3jnnXdo0aIFISEhPPnkk4UmwvO9B7uUZA9pKdZw3FXAPVgJaLGIVHF+D7A/3zn7XY7FAGnG5V+K8/vUfG3OdQ1c2pxFRAaLSJKIJKWlpRX3PZV+xsCMR8AvAHq/8ffj8XdCWA2Y86LVVqlSbunSpWclnCVLlhAbG0tYWBgAUVFR7Nu3L+/4qVOn2LDBu4t7Fi5cyLXXXsttt91Gq1atqF+/Pps2bbI7rGIrsYRkjJlpjPmvMWaVMeZXoI/z9QeVVAzniesjY0y8MSY+Kqoc3nOz7ger93PZMxAa/ffj/gHQ/QnYuww2TC/x8JRyt+TkZB566CE2btzI1KlTeeONNxg+fHje8Z49ezJp0iTmzp3L2rVrufPOO8nOzrYx4sI1atSIOXPmsHDhQjZs2MD999/P9u2lb0DItrJvY0wGsBZoCORW2+X/jRjtciwFiHKdC3J+XzVfm3NdA5c2KtfpDPjlSYhpYfWEzueSWyCyEcz5Nzhyzt9OqVJg4MCB5OTk0KFDB+655x7uuuuusxLSyJEj6dmzJ3379uWKK66gc+fOtG7d2saIC/f000/Tvn17evfuTdeuXQkODmbgwIF2h1VsYteduyISgDWn8wHWHE8y8J4xZpTL8VRghDFmvLOoYR1wqTFmsbNNJ6x5pMbGmI0iMhSriKGqMeaUs82TwH1ADVPIm42PjzdJSUkeeLdeavZzsGgM3DkLanUouO3a760bZvuOg9a3lkx8SqlSQUSWGWPOXeVRDCV5H9KbItJNROqKSAdgKhAMfO5MFGOAx0Wkn4g0ByZgFTFMBjDGrAd+BsaLSIKIJGAVO0w3xmx0vsxk4AQwQUSai0g/4AlgdGHJqNxJ22Tdc9Tq1sKTEVgFD9Xbwm8vQeZxz8enlCp3SnLIrgbwFVY13LdYlXYdjTE7ncdfB94GxgFJQDXgCpd7kMAqAV+JVTn3i/P723IPGmOOAr2AWOc1xmHdfzTaY++qNDIGZo6ACsFw+fNFO0cErhwF6ftg8VhPRqeUKqdsG7LzRuVmyC53+O3qN6H9PcU7d8ptsGUODPsLQs9ZuKiUKmdK3ZCd8hLZp2H2s1C1WcGFDOfT6wXIybSG7pRSyo00IZU3S8fDkZ1w5cvg41v88yvXgw5DYPmXkLLa/fEppcotTUjlyfGDMP9NaHgF1O9x4dfp+qi1vNCsp/VmWaWU22hCKk/mvQqZGdDr3xd3ncAI6Pa4dUPt5oKX1ldKqaLShFRepG2CP/8DbW+Hqo0v/nrxd0Hl+lYvKce772JXSpUOmpDKi9nPgn8QdB/pnuv5VYAr/g0HNkLSf9xzTaVUuaYJqTzYPh82zYSuj0CIG9fri7sa6vWA316GjHK4MK0qdfLv/OoN8u/kWp5pQirrjIFfn7dW7O4w1L3XFoHer0PWcd0zSakL9M477/Dll19e9HUcDgfvvfcerVq1IjAwkLCwMHr27MnMmTPPajd37lz69u1LtWrVCAoKomXLlnz66acX/fruoAmprNs4w1qpu/sT1srd7hbVCDreC8snWjvLKlUOZWZmXvC5lSpVIjw8/KJjGDBgAM888wxDhgxh7dq1JCYm0q5dO/r06cMHH3yQ127x4sW0aNGCqVOnsmbNGoYOHcrgwYOZPHnyRcdw0Ywx+nA+2rZta8qUnGxjxnU05t02xmRnee51Th0z5o1GxozvbkxOjudeR6mL1K1bNzNkyBAzbNgwEx4ebsLDw82jjz5qclz+3dauXdu88cYbfzvvvvvuO6vNc889Z+644w5TqVIl079/f7N9+3YDmKlTp5rLL7/cBAYGmiZNmphZs2YVGNOgQYPMNddcc9ZrDR061IwcOdJUqVLFREVFmUceeeSsGPObMmWKAcx33333t2MPPvigqVChgtm9e/d5z7/hhhtMv379CoyzIECSccPvYO0hlWVrpkHqOujxFPh6cHPgiqFWgUPyX7Di4ocelPKkSZMm4XA4SExMZPz48Xz00UeMGTOm2NcZPXo0jRs3JikpiVGjRuU9/9RTTzFs2DBWrlxJu3btuPnmm8nIyCh2jH5+fixevJixY8cyZswYpkyZUmD7hg0bct111/3t2IgRI8jMzGTatGnnPf/YsWNEREQUK0ZPKOktzFVJycmC31+29jpq+vd/pG7X4gZI+tSar2pyrXWvkipfZj5R8qt3xLSA3q8W65Rq1arx7rvvIiI0btyYTZs2MXr0aB5++OFiXadbt2489thjeT/v2LEDgOHDh3PttdcCMGrUKL744gtWrFhB586di3ztpk2b8uKLLwLW5nsff/wxc+bM4ZZbbjln+02bNtGkSZNzHqtevTphYWFs3LjxnMenT5/OnDlzWLRoUZHj8xTtIZVVyyfC4R3Q81nwKYH/zLkFDicPWxv5KeWlOnbsiMs+nyQkJLB3716OHTtWrOvEx597LdGWLVvmfR8bGwtAampqsa7teo3c6xT3GkWxaNEiBgwYwLvvvkv79u3dfv3i0h5SWZR1Eua9DjU7QsNeJfe61VpC+yGw9EO45Gaoaf8/cFWCitlT8VY+Pj6YfEtiZWVl/a1dcHDwOc/39/fP+z438TkcjmLF4HqN3OsUdI1GjRqxfv36cx7LTbaNGjU66/mFCxdy9dVX8+KLLzJ0qJsrcC+Q9pDKoj//Y+1bdNmzVs+lJPV8CsJi4ceHrGFDpbzM0qVLz0o4S5YsITY2lrCwMACioqLYt29f3vFTp06xYcOGEo+zOAYMGMDmzZv5/vvv/3bs9ddfp0KFCvTv3z/vufnz59O7d2+ef/55HnrooZIMtUCakMqarJOw6B2o2w3qXFryr18xFK5+A1LXWjvSKuVlkpOTeeihh9i4cSNTp07ljTfeYPjw4XnHe/bsyaRJk5g7dy5r167lzjvvJDvbu5fHuvHGG+nfvz+33347H374Idu3b2fdunU88cQTeUURNWrUAKz7kHr37s3//d//MWDAAFJSUkhJSSEtzf6b23XIrqxZ9jkcT4VuE+yLofE10LgPzH3NKqioXNe+WJTKZ+DAgeTk5NChQwdEhLvuuuushDRy5Eh27NhB3759CQkJ4amnniI5OdnGiAsnInz99deMHTuWDz74gOHDh+Pn50d8fDzTp0+nd+/eeW0nTJjAiRMnePPNN3nzzTfznq9du3ZeYYZddMdYF6V+x9isU/BuK2vR0zt+sjeWo3thXAdrHunWaSU/dKiUKjG6Y6z6uxVfWnNH3UbYHQlUqg6XPQNb51j3QymlVCE0IZUV2Zmw4G2o2cGaP/IG7e6G2DYw83Frc0CllCqAJqSyYuVXcGwPdH3Me4bHfHyh71g4dRRmekGvTSnl1TQhlQU5WbDgLas30uAyu6M5W3Qza3fZNdNg3f/sjkYp5cU0IZUFq7+BIzuhmxf1jlx1fghiWsJPD+vQnVLqvDQhlXYOh3XfUXRzaHSV3dGcm68/XPcBnDwCMx8rvL1SqlzShFTabZ4FaRug0zDv7B3limlu9eDWTIX1P9odjVLKC2lCKu0Wv2vtBtu8n92RFK7zcGt15unD4cQhu6NRSnkZTUil2Z4k2LkIEu6zhsW8Xd7Q3WErKelN2UopF5qQSrNF70BAJWjzL7sjKbqYFtB9JKz7HlZ+bXc0SikvogmptDq41ZqLaXc3VAyxO5ri6TwcanWCGSOsPZuUUgpNSKVX4lhrCKz9ELsjKT4fX+g3HsQHvh0MOd69krJSqmRoQiqNMtJg+SS45BYIjbY7mgsTXguueQt2L4WFo+2ORinlBWxLSCIyUkSMiIx1eU5E5HkRSRaRkyIyV0Sa5TsvQkQmishR52OiiITna9NCROY5r7FXRJ4V8eaa6GL682PIyYROD9gdycVpeQO0uAHmvmoVaCilyjVbEpKIdAQGA6vyHXoMeAR4AGgHpAKzRSTUpc1koA1wlfPRBpjocu0wYDaw33mNB4ERwMOeeC8lLvs0JH0Kja6EyIZ2R3Pxrn7T2mF22t1wOt3uaJRSNirxhCQilYBJwJ3AYZfnBXgIeNUYM80YswYYBIQCA5xtmmAlocHGmERjTCIwBOgjInHOSw0EgoBBxpg1xpipwGvAw2Wil7TmWzieBh3+z+5I3CMwHPp9ZC199NMjWgquVDlmRw/pI2CqMeb3fM/XBWKAWblPGGNOAvOBTs6nEoAMYLHLeYuA4/naLHCem+sXIBao4563YBNjYOmHEBkH9brbHY371O5klYKvmgLLJxbeXilVJpVoQhKRe4AGwNPnOBzj/Lo/3/P7XY7FAGnGZZtb5/ep+dqc6xqur+Ea02ARSRKRJG/YU75Au/+AfSugwxDvXiboQnR5xEqyM0bA/rV2R6OUskGJJSTnkNooYIAxJqukXrcwxpiPjDHxxpj4qKgou8Mp2NIPoWIluORmuyNxPx9f6PexdaPvfwfB6Qy7I1JKlbCS7CElAJHAWhHJFpFsoBtwr/P73H0J8tcxRwMpzu9TgCjXuSDn91XztTnXNXBpU/oc3QvrfoA2t0GFYLuj8YyQqnD9J3Boq7VVhc4nKVWulGRC+h5oAbRyeSQBXzu/34SVMHrlniAiAUAXzswZJQIhWMktVwIQnK9NF+e5uXoBycAOd76hEpX0H8BA+8F2R+JZdbvqfJJS5VSJJSRjzBFn1VveA6sY4ZDzZwOMAR4XkX4i0hyYgFXEMNl5jfXAz8B4EUkQkQRgPDDdGLPR+VKTgRPABBFpLiL9gCeA0a5zT6VK1ilYNgHiroaI2nZH43mu80nJK+yORilVQrxtpYbXgbeBcVi9p2rAFcYY1xtUBgArsSrnfnF+f1vuQWPMUaweUazzGuOAt4DSuxzA2u/gxMGy3zvK5eML/T6BoEiYcqvuMqtUOSGltdPgCfHx8SYpyQtXDPjPFVZCuj+p7FXXFWTvX/DpVVCrI9z6Lfj62R2RUuocRGSZMSb+Yq/jbT0kld/+tdZ6b23vKF/JCKB6G+gzGrbPgznP2x2NUsrD9E9Ob5f0GfhWhFYD7I7EHq1vheTlsPg9iG0Nza+3OyKllIdoD8mbZR63qs2aXQdBle2Oxj5XvgI1O8IP90PKGrujUUp5iCYkb7ZmGpw+Zg3XlWd+FeDGz6FiGHw9QIsclCqjNCF5s6TPIKqJNalf3oXGwE1fQnqKVXmXfdruiJRSbqYJyVslr4DkvyC+HBYznE/NdnDd+7BrMUwfris5KFXGaFGDt1r2GfgFQsub7I7Eu7ToDwc2wbzXILIRdH7I7oiUUm6iCckbZR6H1VOheT9rvyB1tu4j4cBm+PV5qFIfmlxrd0RKKTfQITtvtO4HyMyA1rcV3rY8ErGG7qq3gW8Hw76VdkeklHIDTUjeaPkkqFxfixkK4h8IN38FgZVh8k1wdI/dESmlLpImJG9zaDvsXGjdCKvFDAULjYYBU6whzi/7w8nDdkeklLoImpC8zYrJgJTNTfg8IaY53DwJDm6BrwdaK6MrpUolTUjexOGAlV9B/R5QqYbd0ZQedbvCPz+EnYvgu8HgyLE7IqXUBdCE5E22z4Oju6HVQLsjKX1a9IcrR1kFIT+P1HuUlCqFtOzbm6yYBAGVoHEfuyMpnRLug2PJkDgWwmL1HiWlShlNSN7i5BFY/6PVO/IPKLy9Orde/7aWF/r1OQiOgtba21SqtNCE5C3WfgfZp/QX6MXy8bHuUTpxEP53P1QItlZLV0p5PZ1D8harpkBkHMS2sTuS0s+volV5V6M9TLsbNs+2OyKlVBFoQvIGR3bBrkRoeaPee+QuFYJh4H8huqm1OviOhXZHpJQqhCYkb7B6qvW1RX974yhrAirBrd9CeG2YfDPsXWZ3REqpAmhC8garv4GaHSCijt2RlD3BkfCv760dd7+8HvavszsipdR5aEKyW8oaSF0HLW6wO5KyKywWBv0P/ALgi39A6ga7I1JKnYMmJLut/gbEF5r90+5IyraIOjDoR+uz/ryPJiWlvJAmJDs5HNb8UYPLrKEl5VmRDeH26ZqUlPJSmpDstHsJHNsDLW60O5LyQ5OSUl6rSAlJRJqIyIsiMk9EdopIqoisFZGJIjJARCp6OtAyadV/wT8I4nrbHUn5oklJKa9UYEISkTYi8iuwHLgUWAy8CTwJfA4Y4GUgWUQe18RUDDlZsO57aHwNVAyxO5ryJ39S0uo7pWxX2NJB3wFvADcYY867+5mIJADDgUexEpQqzPZ51oZyzfrZHUn5lZuUJvSBCVdb9yxV15UylLJLYUN2DY0xYwtKRgDGmERjzI1YyUsVxbofoEIo1O9pdyTlW2RDuPNnqBgKn/8Ddi62OyKlyq0CE5IxJlNE2hV2ERF5Jbe9uwIr03KyYP10iLtKV/b2BpXrwp2/QFg1mNgPNv9qd0RKlUtFKWqYISJx5zsoIi8BuvFMcexYACcPQVNdhdprhMXC7TMgsgF8dbPVg1VKlaiiJKTpwCwRqZ7/gIg8D4wACl1mQETuE5FVInLM+UgUkWtcjouIPC8iySJyUkTmikizfNeIcFb2HXU+JopIeL42LZzVgCdFZK+IPCviZSuWrv0eKoRY9x8p7xESBYOmQ2xr+OZ2WDHZ7oiUKleKkpDuwqqymy0ilXOfFJFngJHATcaY6UW4zh7gcaANEA/8BnwvIi2dxx8DHgEeANoBqc7XDHW5xmTn+Vc5H22AiS4xhQGzgf3OazyIlTAfLkJ8JSMnGzZMh0ZXgn+g3dGo/ALD4bbvoE4X+H4oLH7P7oiUKjcKTUjGGAdwE1aCmCkiwSLyJPAsMNAY831RXsgY84MxZqYxZosxZpMx5ikgHUhw9mAeAl41xkwzxqwBBgGhwACw7oXCSkKDnUUUicAQoI/LkOJAIAgYZIxZY4yZCrwGPOw1vaSdC63N43S4zntVDIEB/7X+G816Gn5+0lpVQynlUUW6MdYYcxq4FqtMfAXwPHCr8xd+sYmIr4jcDIRg3dtUF4gBZrm85klgPtDJ+VQCkOFsn2sRcDxfmwXOc3P9AsQCdS4kVrdb+711M2yDy+2ORBXEPwD6fwbth8CScfDt3ZB92u6olCrTCt3CXERcb5T5BHgb+AHIcj1mjPm2CNdqASQCAVjJ5Z/GmNUikptQ9uc7ZT+QO3cVA6QZY4zLaxoRSXUey22z5xzXyD22/RwxDQYGA9SqVauwt3BxHDmw/kdruK5CkGdfS108Hx/o/ZpVfffr85CRau1EG1DJ7siUKpMKTUjAuXpB1zsfuQzgW4RrbQRaAZWA/sDnItK9COd5jDHmI+AjgPj4eFNI84uzczGcOABN+3r0ZZQbiUDn4RBaDX64Dz67GgZOtZKUUsqtijKH5FOER1GSEcaYTOcc0jJjzEis4b/hQIqzSXS+U6JdjqUAUa5zQc7vq+Zrc65r4NLGPht+At+K0KCX3ZGo4rrkZhgwBQ5th/9cAanr7Y5IqTLH7tW+fYCKWENpKUDeb2oRCQC6cGbOKBFrzinB5fwEIDhfmy7Oc3P1ApKBHe4PvxiMgY0/Qb3uunZdadXgcrjjJ8g5bSWlLXoDrVLuVNjiqp2LeiERCXHOEZ3v+Ksi0kVE6jjvFXoF6A5Mcs4LjQEeF5F+ItIcmIA1zzQZwBizHvgZGC8iCc7188YD040xG50vMxk4AUwQkebOOa4ngNGuc0+22L8WjuyCxlfbGoa6SLGt4e45EF4LJt0Af3xsd0RKlRmF9ZA+EZE5InKL8x6fvxGRliLyOrAFuKSAa8UAX2LNI83Buk+otzFmpvP461gFE+OAJKAacIUxJt3lGgOAlViVc784v78t96Ax5ihWjyjWeY1xwFvA6ELep+dtnAEINNKtJkq98JrW+ncNr4AZj8KMx6z7y5RSF0UK6jiIiB/WvT4PAA2ArVjDX6eACCAOq2LuW+BlY0ypXsM/Pj7eJCUleebi47uBrz/crcM8ZYYjB2Y9Y5WFN+gF/T+FgHP+3aZUmSYiy4wx8Rd7ncIWV802xowzxjQGOgIfYBUi7AR+Be4GqhtjBpb2ZORRR/fCvhUQp8N1ZYqPL1w1Cq4ZDVt/g0+vhMM77Y5KqVKrKGXfABhjkrCGwVRxbZxhfW18TcHtVOnU7i6oXA/+Owg+7gE3TIC6Xe2OSqlSx+4qu/Jh4wyoXB8iG9kdifKU+j3gnt8gKBK+uA4S37cqK5VSRaYJydNOHYPtC6zqOi9ZTk95SGQDa44wrjf8MhK++z/IOln4eUopQBOS5235FRxZEKfDdeVCQBjcOBG6PwmrvoZPr4Iju+2OSqlSQROSp22cCUFVoGZ7uyNRJcXHB7o/Drd8DYe2wUfdYcdCu6NSyusVmpBE5NN8exKponLkWD2kBr2siixVvsT1tuaVAiPg83/AwjG6jYVSBShKD2kQoDvJXYjk5dZW5Q117bpyK7KhlZSaXAu/PgdfD4CTh+2OSimvVJSEpDPxF2rzLBAfqN/T7kiUnQLCrFLw3m9YPebxXWHvMrujUsrrFHUOSetXL8Tm2VA9HoIqF95WlW0i0GEw3PmLVQ7+6VXWOnhaGq5UnqImpBQRySno4dEoS6OMNEj+y1rvTKlcNdrCkPlQr4e1Dt7UO+F0euHnKVUOFHWlhsHAEU8GUuZsnWN9bahblat8gipbFXiLxsBv/7aWlbr+E6je1u7IlLJVURPSj8aYVI9GUtZsngXBURBT0ALoqtzy8YEuD0OtjjDtHmt/pZ7PQKdh1jGlyqGi/MvXQe7icuRYi2026KW/XFTBaneCoQutdQ5/fQ4m9oVj++yOSilbaJWdJ+xdZpX26nCdKorACLjhc/jHe7AnCT7oBBtm2B2VUiWu0IRkjPHR4bpi0nJvVVwi0OZfMHgeVKoOX98CPz2ia+GpckXHkzxh629Qo531l69SxRHVyNoiveN98Ocn1saOe/+yOyqlSoQmJHc7edhaoaFeD7sjUaWVX0Vr479bp8HpY/DJ5fD7K5CTZXdkSnmUJiR3274AjAPqdbc7ElXaNbgc7k2E5tfDvFetxJS6we6olPIYTUjutm0uVAiBGhe9vbxS1rDv9R9bRQ9HdlnLDi0eq4u0qjJJE5K7bZtLcngb5m3V+4iVGzW7Du5dYhXKzHoKPu8Dh3fYHZVSbqUJyZ2O7IJDW/l4b20GffoHb8/ehMOht3EpNwmNhluVwXj2AAAgAElEQVS+gr7jYN8qeL8TLP1Ie0uqzNCE5E7b5gKw0NGC5tXDeGfOZh74ajknM3WpP+UmItD6Vrh3MdTqADNHwISr4cBmuyNT6qJpQnKnbXM56luFo8H1+PH+zjx5dWNmrNnHDeMXs++o3k+i3Ci8Ftz6LVz3AaSuhw8uhQVvaSWeKtU0IbmLw4HZNo/5Oc3o0qgqIsLgrvX55F/xbE87Tt+xi0jaccjuKFVZIgKtBsB9f0DcVTDnRfi4J+xbaXdkSl0QTUjusn8NcuIAv2c2pWujyLynL2sSzbf3XkpgBV9u/mgJnyzYhtE9cJQ7hUbDjV/AjRMhPQU+6gG/vgBZp+yOTKli0YTkLs75o8WmOZ0bRJ51KC4mlB8f6MxlTary0k/rGfrlXxw7pUMrys2a/gPu/wMuuQUWjrbWxNv6u91RKVVkmpDcZfs89vjWJCq2LlVCKv7tcFiAPx/e2panrm7C7PX7+cd7C1mbfNSGQFWZFhgB142D274HDEy8DqbdDen77Y5MqUJpQnKHnGzMrkTmZsbRpWHkeZuJCPd0rcfXgztyMiuHf45bzMfzt2lpuHK/+j1gaCJ0ewLW/QBj21lbpju04lN5L01I7rBvJZJ5nCU5Tf42XHcu7epUZuaDXekeF8XLM9Zz26dLSTmq4/3KzfwDoMdIKzHFtrK2TP/kckheYXdkSp2TJiR32LkQgOU+zWhTu2grfFcOrsD429ryar8W/LXzCFe9M5+Zq3VjNuUBkQ3gXz9Av0/g6G74uAfMfAJOHbM7MqXOUmIJSURGisifInJMRNJE5EcRaZ6vjYjI8yKSLCInRWSuiDTL1yZCRCaKyFHnY6KIhOdr00JE5jmvsVdEnhURz200uGMRu32qU6dOXQL8fYt8mohwc/ta/DSsM7UqBzF00l8M+2o5h45neixUVU6JQMsb4P4kiL8Tln4IY+Nh5de60oPyGiXZQ+oOvA90AnoC2cCvIlLZpc1jwCPAA0A7IBWYLSKhLm0mA22Aq5yPNsDE3IMiEgbMBvY7r/EgMAJ42BNvCkcOjp2LWZAZR6f6hQ/XnUu9qBCmDe3EQ5c3ZOaaffQaPY/pq5K1PFy5X2A4XPMW3DMHKtWA74bAp1fqnkvKK5RYQjLGXGmM+cwYs8YYsxq4DYgCLgWrdwQ8BLxqjJlmjFkDDAJCgQHONk2wktBgY0yiMSYRGAL0EZE450sNBIKAQc7Xmgq8BjzskV5Syip8MtNZ4mjKpUWYPzoff18fHrq8ET8+0JnqEYHcP3k5QyYuI/WYzi0pD6jeFu76Ffq+by3S+nFP+N8DkJFmd2SqHLNzDinU+fqHnT/XBWKAWbkNjDEngflYvSqABCADWOxynUXA8XxtFjjPzfULEAvUces7ANixCIA1/s1pHht20ZdrHBPGt0M7MbJ3Y+ZtSuPy0fP4cslOcrQST7mbjw+0HggPJEHCfbBiMrzXFpZ8oEsQKVvYmZDeAVYAic6fY5xf898wsd/lWAyQZlzGspzfp+Zrc65ruL5GHhEZLCJJIpKUlnYBfx3uXMQeqUa9eg3x83XPx+nn68OQbvWZ+WAXmsaG8fT3a/jn+4tYuVu3tFAeEFAJrnzZqsarEQ8/PwEfdtabalWJsyUhichooDNwvTHG1hsjjDEfGWPijTHxUVFRxTvZkYNjxyIWZDXm0gZV3B5bvagQvrqnI+/c3IqUo6e47v1FPPndag5r0YPyhKhG1rbpt3wN2aesm2q/ukVXElclpsQTkoi8DdwC9DTGbHM5lOL8Gp3vlGiXYylAlOtckPP7qvnanOsarq/hHvvX4nP6KEsdTS5q/qggIkLfVtWZ80g37ry0LlP+3E3Pt+YyaelOsnO0Okq5mQjE9YZ7l8Jlz8H2BfB+R5gxAo4fsDs6VcaVaEISkXc4k4w25Du8HSth9HJpHwB04cycUSIQgjVPlCsBCM7Xpovz3Fy9gGRgh1veSK6d1vzRlqBLaFg1xK2Xzi80wJ9n+jTlp2GdaRgdylPfraH3Owv4fUOqVuMp9/MPgC4Pw7Dl0GYQ/PkfeLc1LByji7YqjynJ+5DGAXdgVcwdFpEY5yME8uaCxgCPi0g/5z1KE7CKGCY726wHfgbGi0iCiCQA44HpxpiNzpeaDJwAJohIcxHpBzwBjDZu/s1tdi5iL1Wp3yAOT97m5KpxTBhTBnfkw1vbkpXj4I4Jf3Lbf/7QdfGUZ4REQZ/RcG8i1O4Evz5nLUO0eiroH0LKzUqyh3QvVmXdHGCfy+NRlzavA28D44AkoBpwhTEm3aXNAGAlVuXcL87vb8s9aIw5itUjinVeYxzwFjDare/GGHJ2LGFJThwd67l//qggIsJVzWOYNbwbz13blDXJR+nz3kIe/WYle4/oRoDKA6LiYMAU+Nf/ILASTLsLPrkMdiYWfq5SRSQ63HNGfHy8SUpKKlrjQ9vg3dY8mXUXdz/0IvWiPDtkV5CjJ7MY9/sWJizaAcDN7WtyX48GRIcFFHyiUhfC4YBVX8Ocf0N6MjTuAz2fgaqN7Y5M2URElhlj4i/2OrqW3YXa/QcAWwOaUTcy2NZQKgX68+TVTfh9RHeub1udyUt30fX133lp+joOZJy2NTZVBvn4WDvVPrAMejwN2+fDBwnw/b1wZJfd0alSTBPShdq9lAyCiKp7SYnNHxWmenggr/RryW+PdKdPy1g+XbSdrq//zms/b9D18ZT7VQiCbiNg2AroeK81r/ReW/h5pFbkqQuiQ3YuijNkl/leRxJT/dnReyKDOtXxbGAXaGtaBu/8upkfVyUT4OfLze1rck+XesSGB9odmiqLju6Bua/CikngHwSdHrBWgKgYWvi5qlTTITs7nTqK/8ENLHM0on3dyoW3t0n9qBDevaU1s4d35eoW1ZiYuJNub/zOiG9WsjUtw+7wVFlTqQb0HWvdw1S/J8x9Bd65xFqKKFuHjlXhNCFdiD1JCIb1/k2Ii/b+v/4aVA3lrRsvYe6I7gxoX4v/rUzm8tHzGPrlMl2OSLlfVCO4aSLc/RtEN7OWInqvLfz1ha6RpwqkCelC7F5KDj5UqN0OHx/vmD8qihoRQbzQtzmLnujJvd3rs3DLAfqOW8T1Hyzmp1X7dOUH5V412sKgH+G27yE4ylpNfGy8tYhrTrbd0SkvpHNILoo6h5T56bVs3rGTRZd/x+Cu9UsgMs9IP5XFN0l7mLB4B7sOnaB6eCCDOtXmpvhaVArytzs8VZYYA5t+gbmjYN9KqFwfuj0OLfqDT9E3tVTeyV1zSJqQXBQpITlyyB5Vg8mnLqXlkE9oVTO84PalQI7DMGf9fj5btIPEbQcJ9Pelf9sa/CuhNg1LwZCkKkWMgY0z4PdXYP9qiGxkJaZm/9TEVIppUYNd9q/FL/sEq30a08wN+x95A18f4YpmMXw1uCMzhnWhT8tqTPlzN73ens+NHyby/fK9nMqydVF2VVaIQONrYMh8uPELEF9r1YcPOsHa73Q79XJOE1Jx7V4KQHb1dvi7af8jb9I0Now3briExJE9Gdm7Manpp3hoygoSXpnDS9PXaXWecg8fH2jaF4Yuhv6fWT2nb2639mFa9z9NTOWUDtm5KMqQXeZ/7+Lw2l+Z0mU2wy5vVEKR2cfhMCRuO8jkpbv4ZW0K2Q5Dx3qVuaV9La5sFkOAvw6zKDdw5MCab2Heq3BwC0Q1ga6P6lBeKaFzSB5QlIR04s2WzD9albBBX9PJQ3sgeavU9FN8k7SHr/7YxZ7DJwmt6EefS6pxfZsatK0d4TUrVqhSLCfbGrpb8CakbbCKH7o8DC1vAl8ttPFWmpA8oNCEdOIQvF6X17Nv5r5nxhFc0a/kgvMiDodhybaDTP1rDzNXp3AyK4e6kcH0a12dfm1rUF1XglAXy+GADdNh/huQsgoq1YLOD0KrW629mpRX0YTkAYUmpC2/wpfXMzLkZV559P6SC8yLZZzOZubqfUz7aw9Lth1CBBLqVeG61tW5slkMlQL1r1p1EYyBzbNh/uuw508IiYFLh0Hb26GCvYsaqzM0IXlAYQnJMfd1+H0UL7ecyTPXJ5y3XXm1+9AJvv1rL98u38POgyeo4OtD10ZRXHtJNS5vEl1ue5TKDYyxVhWf/wbsWABBkdY6ee3uhoCyUe1ammlC8oDCElL6Z9eTsn09a/45i3+2rlGCkZUuxhhW7TnKjyuTmb5qHynHThHg78NlTaK5tmUs3eOitBhCXbhdS2D+m7BlNgRUgvZDoMMQCC5fc7reRBOSBxSYkIzh5Cv1+elkM9o99DW1q+hwQVE4HIaknYf538q9zFidwqHjmYRU9KNX02iubBZDt0ZRBFbQ5KQuQPJyKzFtmA5+gdDmNki4HyJq2x1ZuaMJyQMKTEhHdsOY5rwqd/H4s29pRdkFyM5xsHjrQX5cmcysdfs5ejKLAH8fujWK4spmMVzWOFqXLFLFl7YRFr0Lq6aAcUDz6+HSByGmud2RlRvuSkg6qF9Ue5cBkFOtjSajC+TnnFPq2iiKUTkO/th+iF/Wpjgf+/HzERLqV+GKZjFc0TRat2BXRRMVB9eNgx5PwpL3YdkEWP1faNALOj8EtS+1VohQXk97SC4K6iGdnPE0Pks/YEK3hQzp2aSEIyvbHA7Dyj1H+GXtfmatTWHbgeMAtKoZzmWNq9KjcVWaxYbpHwKqaE4ehj8/gSUfwokDUKMdXPoQxF1trRCh3E6H7DygoIR0+P0r2JlygMw7fvXqTflKO2MMW1Iz+GVtCrPXp7JqzxGMgeiwivSIs5JT5waRWrGnCpd1EpZ/CYvfgyM7rYVcL30QWtwIfhXsjq5M0YTkAedNSI4cMl+qwddZnbnx2a+1QqwEpaWfZt6mNH7fkMr8TWmkn86mgq8PHepVpkdcVXo2rkqdSC0wUQXIyYZ138PCMdYK46GxkHAvtBmkJeNuognJA86bkFLXw/sdeTvkYYY/+lzJB6YAyMpxkLTjML9vTOW3DalsSbUWeq1TJYjODSPp0jCKhPpVCAvQwgh1DsbA1jlWYtqxACqGQdtB0OH/rO3X1QXTooYS5NizDB/At0Zbu0Mp1/x9fUioX4WE+lV48uom7Dp4gt83Wj2n7/7ay5dLduHrI7SqGU7nBpF0aRjJJTXDy+Sq7OoCiECDy63H3r8gcSwkvg9LPoBm/aDT/VDtErujLNe0h+TifD2kI9OG47dqMjP7/MEN7fQeB2+Ume1g+a7DLNxygAWbD7BqzxEcBkIr+tGxfhW6NIykc4NI6kYGa3GEOuPILqv44a/PITMD6naFhAespKUFEEWmQ3YecL6EdODdHmw/kEHYvb8RF6M7qJYGR09ksXjrARZsOcCCzWnsPnQSsIojOtarQsd6VUioV4XaVYI0QSk4ecRKSks+hPRkiGpsLU3U4kZdzLUINCF5wDkTksPB6Zeq801OV2557mt8ffSXV2m08+BxFm45wJJth1iy7SBp6acBiAkLoGO9ynlJShNUOZedaW1/kfgepKyG4ChraaJ2d0GQVteej84hlZRD26joOMGxSk01GZVitasEU7tKMAM71MYYw7YDx0ncepAl2w6ycMtBvl+RDJydoNrXraxDfOWNXwW45CZoeSNsnweLx8LvL8GCt6DVAKvXVKW+3VGWWZqQCpG9dzl+gF+NVnaHotxERKgfFUL9qBBu7WglqK1px0nclpugDuQlqCrBFWhbO4L4OhG0rV2ZFtUrUcFP5xbKPBGo1916pK63CiCWT4SkT6HxNdDpAajZQVeAcDNNSIU4si2JUONHtYaakMoqEaFB1RAaVA3htrwElcGfOw6TtOMwSTsPMWvdfgAq+vlwSc1w4mtH0K5OZdrUitD198q6qk2g7zjo+Sz88ZG1CsSG6VA93kpMTa7VbdbdpETnkESkK/Ao0BaIBe4wxkxwOS7Ac8BgIAJYCtxnjFnr0iYCeBf4h/Op/wEPGGOOuLRpAYwF2gOHgPHAv00hb/Zcc0gp711Jatp+woYt0hswy7HU9FP8tfOwlaR2Hmbt3qNkO6x/TnHRobStE0HbWhG0qhVO3SrB+OjwbtmVeRxWTLZ6TYd3QHhtayiv1UCoGGJ3dLYolUUNInI10Bn4C/gCuDdfQnoceBq4HdgIPOtsH2eMSXe2mQnUAu52nvYJsM0Yc63zeBiwCZgPvAg0Bj4DnjfGvFVQfH9LSMZw/KXa/JwTT7/npulcgspzMjOHFbuPsGznIf7ccZi/dh0m/VQ2AGEBflxSM5zWNcNpXSuCS2qGUzlYl6opcxw5sOEna2miPX9YezPF32kVQYRVszu6ElUqE9JZLyySAdyfm5CcvaNkYKwx5mXnc4FAKvCoMWa8iDQB1gGdjTGLnG06AwuAxsaYjSIyFHgNiDbGnHS2eRoYCtQoqJf0t4Tk3HJiQsQD3P7gS+79AFSZkuOwhvlW7DrC8t1HWL7rMJv2p+PsRFG7ShCtaobnPZrGhlHRT4d5yozdf1iJaf2P4OMHLW6wbrSNbmZ3ZCWiLFbZ1QVigFm5TxhjTorIfKAT1rBbApABLHY5bxFw3Nlmo7PNgtxk5PQL8G+gDrC9qAFl7llBBcCvut69rQrm6yM0ig6lUXQoN7arCcDx09ms3nuUFbuPsGLXEZZuO8QPzmKJCr4+NIkNo3XNcFrWqESL6pWoFxWilZylVc32cNNEOLTNWvlh+ZewcjLU72ltGli/pxZAFIE3JaQY59f9+Z7fD1R3aZPm2ssxxhgRSXU5PwbYc45r5B47KyGJyGCsOStq1ap11kkHtvxBtBGiG1504lflUHBFv7z7m3KlHD3Fit2HWe7sSU35czcTFu8AIKiCL81iw2he3UpQLWtUom6kJqlSpXI9uPoN6D7Sqsj74yP4sh9UbWb1mJr315XGC+BNCckWxpiPgI/AGrJzPZazdyVbTSzNasec81yliiumUgBXVarGVc2tOYbsHAdb046zeu9R1uw9yuq9R/nqj118luUANEmVWkGVoeujVhXe6m+s+5m+Hwq/vgAdhkD8HRAYYXeUXsebElKK82s0sMvl+WiXYylAlIhIbi/JOfdUNV+b6HzXjnY5VmRhh9exyrcJV1fSpUOUZ/j5+hAXE0pcTCj921orThc3STWNDaN+VIguIuuN/CpC61utCrytc6x5pjkvwPw3oc1t0HEoRNSxO0qv4U0JaTtWwugF/AkgIgFAF2CEs00iEII1T5Q7j5QABLv8nAi8JiIBxphTzud6YRVM7ChyNMcPUik7jSOV/qnVdapEXUiSquDrQ6OYEJpWC6NJtTDra2yYbsXhLVxXGk9ZbfWY/vzEGtJrci10GgY1dGqgRBOSiIQADZw/+gC1RKQVcMgYs0tExgBPisgGrNLtp7GKGCYDGGPWi8jPwHjn3A9YxQ7TjTEbnT9PxrqXaYKIvAQ0Ap4AXijsPiRX2ftW4wdITPOLeMdKucf5ktS2A8dZv+8Y65KPsW7fMeasT+W/SWemUGtWDqRptTCaVqtEk2qhNI0No3p4oP6RZaeYFtBvPFz+HCwdD0mfwbofoGZHa54p7upye6NtSd+H1B34/RyHPjfG3O5yY+wQzr4xdo3LNSKA9zj7xtj7z3Fj7DisG2MPAx8CLxbnxtj9s94mevHzzLxyPr0TtMpOlQ7GGNLST7PWmaTW77MS1fYDx8n91x8W4EfTWCtJNY0No0m1UBpUDdEydLucTreq8hLfh6O7rMKIjvdaw3wVguyOrkhK/X1I3sg1Ie349E5Cd87i0NB1NIzRbY5V6XYiM5sNKeln9aY27EvnZFYOYJWt16kSROOYMOJirPL1xjGh1KocpKtOlJScbFj/P2ueKfkvCKxsrTLefjCEVLU7ugJpQvIA14S0941O7E6H+OcW4qeTxaoMynEYdhy0hvw2paSzISWdjfvT2XXoRF5vKtDfl4bRIcRFh+YNGcbFhBIVUlGH/TzFGNiVaM0zbZwBvv7Q8ibrfqaqje2O7pzK4o2x3sPhoMqJrfwVdIUmI1Vm+fqcWfWclmeeP5GZzeb9GWx0JqlN+9P5fWMa3yw7MzdVObjC35JUo+hQQirqr5SLJgK1O1mPA1tgyThr7bzlEyHuGujycJktgNB/PedgjuwkwJwis4p3/jWilCcFVbDW4rukZvhZzx/MOM3G/elsTHE+9qfzTdJujmfm5LWpERFIw6ohNIwOzfvaoGqIJqoLFdkA+rwNPZ6GP8ZbRRAbf7K2Wu/8sLU9Rhnqqeq/knM4vH0FlYGgGi3sDkUpr1ElpCKdQirSqX5k3nMOh2HvkZN5CWpDSjpbUjNYtPUgmdmOvHaxlQJokJukqobQMDqEBlGhunVHUQVXgR5PWjfaLptgDedNvA5i21g9prhrwKf0j+ZoQjqH3IQU07C13aEo5dV8fISalYOoWTmIy5ueuR89x2HYfegEm1Mz2Jyazpb9GWxOzWDy0l15hRQAVUMr0jA6hIZVrZ5Ubq9KV0c/j4qhVlJqdw+s/AoWjYEpt0JkHHQeDi36W3NOpZQWNbjILWrYNLY/FdNWUuXJ9TrUoJQb5faoNqems9mZpDanZrBlf/pZQ39VgitYCcqZrBpGWxsoajFFPjnZsO57WDAaUtdCpVpw6TBrdQj/wBILQ6vsPCA3Ie17uSVbHTF0fmZW4ScppS6aMYZ9R09ZCWq/Ney3OTWDTfvT8/aZAggN8KN+lJWcrIKMYOpXDaFW5aDyvXSSMbDpF1g4GnYvhZBoq8fU9vYSSUxaZecp2aepmrWbvyK62B2JUuWGiBAbHkhseCDdGkXlPZ97o29uctqWdpytaRks2JzGVJeqPz8foXaVICtJ5UtW5WL5JBGIuwoaXQk7F8HcV+HnJ2DhGGuOqc0g8Pf+NTk1IeVzet96KuLAVG1qdyhKlXsiQtWwAKqGBXBpg8izjqWfyspLUFtSM9ialsHWtOP8tiE1b3t5gKjQijSICqF+1eC8Mvf6VUOoFhZQ9m76FYE6neH26bB9Acx9BWY+5pKY/mUt+OqlNCHlk7ZtOTWAkJotC22rlLJPaID/OcvTs3Ic7D50gq3OZLXVmaz+tyKZYy7Df4H+vtSLck1S1vd1I4MJ8C8DyyjV7WIlp+3zrcQ041FY+LaVmFrf5pWJSeeQXMTHx5tJw3tQb9N/2DFkCw1jK9sdklLKTYwxHMjIdPakMtiaejzv+71HTuatTiECNSOCrCG/fEOAlYMrlM6iCmNg21wrMe1eahU/9Hza2mrdDeXiWtTgAfHx8WbSv6rjd3ATtZ5dW74nSZUqR05m5rD9wN+H/7alZXDa5X6q8CD/M/NTLsN/NSMCS8eqLsbA1t/g1+chZRVEt4Bez0P9yy7qBltNSB4QHx9vfvhnJjtMNS59+he7w1FK2Sy3TH1LWkbefJU1BHicAxmn89r5+wp1qgSfNfRXPyqEelHBhHpjUYXDAWu/hTkvwpGdULcb9HoBYi/s3kutsvMIQ9WsvayNuNTuQJRSXsD1xt8ecWcfO3oii60HziSorWkZbEpNZ/b6/eS4FFVEh1U805tyVv7VjwqhWqUA+4b/fHysm2ibXAtJn8K81+Gj7tD8erj8BQivaUtYmpBcOLJO4082EtXQ7lCUUl6uUpA/bWpF0KZWxFnPZ2Y72HXoxN/mqr5fsfese6qCKuQrqnD2rupUKcGiCr+K1jbqrQbAonchcSxsmGEVPnQaVuKl4jpk56JF4/pm9c0HSOzxNQndetsdjlKqDDHGkJZx+qxiiq1px9maahVV5LK1qOLILpj1tLWDbXhtuHIUNL6m0PklnUPygGb1Y83a246z46611KlZw+5wlFLlxMnMHLYdOJOgbC+q2DYPZj4Oaeuhfk/o/Ya18vh5aELygKa1qph5dwRT+bmd+Ja1G+aUUqVOblHFVteiCmeySkv3cFFFTjYk/Qd+exmyT0G3EdDpQfD7+8K3mpA8oHmNEPPV/11Ci6cX2R2KUkoV6OjJLLalZfztBuCdB0+ctVJFdFhF6kVexEoV6fvh58dh7XdQtRn8492/bRCoVXYe4O/IJD2krt1hKKVUoSoF+tO6VgSt8xVVZOU4iypcqv+2pmXww4rks4oqirxSRWg03DABWtwIPz0Cn1wOHYZAz2egYohb35MmJBe+5OCoohV2SqnSy9/X58zW9C7Ot1LFX7sO8+Oq5LyVKnwEajuH/6yNFK2v9eteQfB9S2HOC9bOtZt+hus+hNoJbotdE1I+QbFN7A5BKaXcTkSICq1IVGhFOtarctax3JUqtjhXqtiSam0BMm9TKlk5Z4b/qocH0qDqTfSMa8v1u14i+LPeZLa/z30x6hzSGfGxvuaHhaupXk9X+lZKqawcBzsPnshLUptTcxNWBn7Zx3nKbxID/H5DXjimRQ3uJiLpwEa74/ASkcABu4PwEvpZnKGfxRn6WZwRZ4wJvdiL6JDd2Ta6I8uXBSKSpJ+FRT+LM/SzOEM/izNEJMkd1ykFy9MqpZQqDzQhKaWU8gqakM72kd0BeBH9LM7Qz+IM/SzO0M/iDLd8FlrUoJRSyitoD0kppZRX0ISklFLKK5SrhCQi94rIdhE5JSLLRKRLIe27OdudEpFtIvJ/JRWrpxXnsxCRfiIyS0TSRCRdRJaKyD9KMl5PKu6/C5fzOotItois8XSMJeUC/h+pICIvOs85LSK7RGRYScXrSRfwWQwQkRUickJEUkTkSxGJKal4PUVEuorI/0Rkr4gYEbm9COe0EJF5InLSed6zUoSNnMpNQhKRm4B3gFFAa2AxMFNEap2nfV1ghrNda+AV4D0Rub5kIvac4n4WQDfgN+AaZ/sZwHdF/cXtzS7gs8g9LwL4Apjj8SBLyAV+Fl8DVwGDgTjgBmCVh0P1uAv4fXEpMBH4HGgGXAc0BSaVSMCeFQKsAR4EThbSFhEJA2YD+4F2zvNGAA8X+krGmHLxAJYCH+d7bjPwynnavwZsztAgSEoAAAUOSURBVPfcJ0Ci3e+lpD+L81zjD+Atu9+LXZ8F8C3wHPA8sMbu92HHZwFcARwFIu2O3Qs+i0eBnfmeuwPIsPu9uPlzyQBuL6TNUOAYEOjy3NPAXpyFdOd7lIsekohUANoCs/IdmgV0Os9pCedo/wsQLyIXuOOV/S7wsziXUOCwu+Kyw4V+FiJyLxD9/+3dX4gVZRjH8e9jRn8wghJ0RdzyIiha6lbs382C9AeijbKI9MKg0i76A1EXlRl0IYEWRV2EixBUehdBBtGWyWJlUBEVJWu0m1t5E2WZqzxdvO/SdNpRdjgz7zvn/D5wOH9mZveZhznzzPvOe2aAZ+qLrlkVc3Ez8AnwkJlNmtl3Zva8mXX3ngQNq5iLfcCAmd1kwWJgLaE3od+sAva6e7E1tQdYBlx0qgX7oiARrjl1BqEJWfQzUNbHu7Rk/oXx77VVlVz8h5ltBJYTuijabN65MLMhQsvoLnc/WW94jaqyXawErgKuAEaATYTuu9F6QmzMvHPh7uOEAvQacBz4FTBgXX1hZqts3zk7rVS/FCTpkngObStwp7v/kDqeJpnZWcAbwCPuPpE6ngwsAJywLex39z2EojRiZkvShtYsM7sMeAHYQmhdrSHsfF9JGVfb9MvFVY8AJwndLEVLgOmSZaZL5j9Bu6/wWyUXAJjZrYQT+Xe7+1v1hNeo+eZiALgU2GFmO+JnCwAzsxPA9e7e2c3TFlW2i8PAlLv/Vvjs6/i8gv8fJbdFlVw8Bnzs7lvj+y/M7Ciw18wed/fJekLNUtm+c3Zaqb5oIbn7ceAAMNwxaZgwemYu4yXzf+ruM92NsDkVc4GZ3Uboolvv7rvri7A5FXIxBQwBVxYeLwPfx9el+ctdxe1iH7Cs45zRJfG5ta3nirk4l1DEimbf98V+tmAcuNrMzi58Ngz8BBw65ZKpR200ODrkdkLf7gbCUe52woiRwTh9J7CzMP/FwFFgW5x/Q1x+JPW6JMjFWmCGMHxzaeFxQep1aToXcyz/FL0zym6+28Ui4EdgF2Go82rC8OBdqdclQS7Wx+/IfYRza6sJAz4OpF6XLuRiEf8egP0JPBFfr4jTnwXeK8x/PqEl9DpwOXALYdTdw6f9X6lXtuHE3k+o0H8TjoCuKUwbA8Y65r8W+CzOPwHcm3odUuQivvc5HmNNx506F3Ms2zMFqUouCL89ejfuqKaAF4HzUq9Holw8AHwVc3GYMMBheer16EIeriv5/o/G6aPAoY5lhoAPgWMxF09ymiHf7q6Lq4qISB76rW9TREQypYIkIiJZUEESEZEsqCCJiEgWVJBERCQLKkgiIpIFFSSRzJnZWLwxWvExmjoukW7T75BEMmdmq4DdwPvAm8C3wC/u3urbf4h0UkESyZyZvQN87u6Ppo5FpE7qshPJ3yDhQqYiPa1fbj8h0mYfAZvj7Qz2u/sfqQMSqYO67EQyZ2bnAK8Cd8SPpt19IGFIIrVQC0kkf1uAG4CngQ8IV08W6TlqIYlkzMxWAgeBG9397dTxiNRJgxpE8nZhfJ5IGoVIA9RCEsmYmZ0JfAMcATbH15MebrMt0lPUQhLJmLvPEO5c/CXwEqEg/W5m9yQNTKQGaiGJtIyZrQOec/fFqWMR6Sa1kERaxMwWEs4r/ZU6FpFuU0ESaZdB4EFgY+pARLpNXXYiIpIFtZBERCQLKkgiIpIFFSQREcmCCpKIiGRBBUlERLKggiQiIln4BxtWR3559ACUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x113c2a588>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "#########################################################\n",
    "\n",
    "def Tad_all_mixf(oxidizer, fuel, stoicReac, stoicProd, plotLabel):\n",
    "    '''\n",
    "    Plot Tad versus mixture fraction\n",
    "    oxidizer,  input, dictionary of species names and moles\n",
    "    fuel,      input, dictionary of species names and moles\n",
    "    stoicReac, input, dictionary of species names and moles\n",
    "    stoicProd, input, dictionary of species names and moles\n",
    "    plotLabel, input, string for the plot label\n",
    "    '''\n",
    "    \n",
    "    species = [\"CH4\", \"O2\", \"N2\", \"CO2\", \"H2O\"]\n",
    "\n",
    "    #--------- Get the state at ξ = 0\n",
    "    gas.TPX = 298.15, P, oxidizer\n",
    "    y0 = gas.Y\n",
    "    h0 = gas.enthalpy_mass \n",
    "    \n",
    "    #--------- Get the state at ξ = 1\n",
    "    gas.TPX = 298.15, P, fuel\n",
    "    y1 = gas.Y\n",
    "    h1 = gas.enthalpy_mass\n",
    "    \n",
    "    #----------- Stoichiometric products\n",
    "    gas.X = stoicProd\n",
    "    yst = gas.Y\n",
    "    \n",
    "    #--------- Get Tad at all ξ\n",
    "    \n",
    "    ξ_st = gas.molecular_weights[gas.species_index(\"CH4\")] * stoicReac[\"CH4\"]/ (          \n",
    "           gas.molecular_weights[gas.species_index(\"CH4\")] * stoicReac[\"CH4\"] + \n",
    "           gas.molecular_weights[gas.species_index(\"O2\")]  * stoicReac[\"O2\"] +\n",
    "           gas.molecular_weights[gas.species_index(\"N2\")]  * stoicReac[\"N2\"])\n",
    "    \n",
    "    n = 1000\n",
    "    ξ = np.linspace(0.0,1,n)\n",
    "    T = np.zeros(n)\n",
    "    \n",
    "    for i in range(n):\n",
    "        if ξ[i] <= ξ_st:\n",
    "            y = y0 + (yst-y0)*ξ[i]/ξ_st\n",
    "        else:\n",
    "            y = yst + (y1-yst)*(ξ[i]-ξ_st)/(1-ξ_st)\n",
    "        h = h0*(1-ξ[i]) + h1*ξ[i]\n",
    "        gas.HPY = h, P, y\n",
    "        T[i] = gas.T\n",
    "        \n",
    "    #--------- Plot results\n",
    "        \n",
    "    plt.rc('font', size=14)\n",
    "    plt.plot(ξ,T, label=plotLabel)\n",
    "    plt.xlabel('ξ')\n",
    "    plt.ylabel('T (K)')\n",
    "    plt.xlim([0,1])\n",
    "    plt.ylim([298.15,5500])\n",
    "    plt.legend(frameon=False);\n",
    "    \n",
    "#########################################################\n",
    "\n",
    "#---- Do combustion in air\n",
    "\n",
    "oxidizer  = {\"O2\":1, \"N2\":3.76}\n",
    "fuel      = {\"CH4\":1}\n",
    "stoicReac = {\"CH4\":1, \"O2\":2, \"N2\":7.52}\n",
    "stoicProd = {\"CO2\":1, \"H2O\":2, \"N2\":7.52}\n",
    "\n",
    "Tad_all_mixf(oxidizer, fuel, stoicReac, stoicProd, 'burn in air')\n",
    "\n",
    "#---- Do combustion in o2\n",
    "\n",
    "oxidizer  = {\"O2\":1, \"N2\":0}\n",
    "fuel      = {\"CH4\":1}\n",
    "stoicReac = {\"CH4\":1, \"O2\":2, \"N2\":0}\n",
    "stoicProd = {\"CO2\":1, \"H2O\":2, \"N2\":0}\n",
    "\n",
    "Tad_all_mixf(oxidizer, fuel, stoicReac, stoicProd, 'burn in O2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}