Notes_Equilibrium_Example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import all the typical toolboxes...\n",
    "from scipy.optimize import fsolve\n",
    "from scipy.optimize import least_squares\n",
    "from scipy.interpolate import interp1d\n",
    "from numpy import*\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import cantera as ct"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simple equilibrium example:\n",
    "\n",
    "Determine the equilbrium concentrations for the following reaction at 1800K and 1 atm.\n",
    "\n",
    "$$ \\frac{1}{2} N_2 + \\frac{1}{2}O_2 \\leftrightarrow NO$$\n",
    "\n",
    "Note that for this reaction, we'll need to look up the $K_{p,NO}$ from our JANAF Tables -- for this you can use the table in the Notes folder on Teams or you can refer to the JANAF database kept up by NIST, here: https://janaf.nist.gov/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_N2 = 0.7878\n",
      "x_O2 = 0.2078\n",
      "ppmv_NO = 4419 ppmv\n"
     ]
    }
   ],
   "source": [
    "# Set up the knowns\n",
    "n_o2i = 0.21\n",
    "n_n2i = 0.79\n",
    "lnkp_no = -4.517 #note here we're using the value from the NIST site, which is slighty different from the handout\n",
    "\n",
    "# Define our equations to solve\n",
    "def eqn(vars):\n",
    "    n2, o2, no = vars\n",
    "    eq1 = 2*n_o2i - (2*o2 + no)\n",
    "    eq2 = 2*n_n2i - (2*n2 + no)\n",
    "    eq3 = exp(lnkp_no) - (no / (n2**0.5 * o2**0.5))\n",
    "    return [eq1, eq2, eq3]\n",
    "\n",
    "# Now call the fsolve function to solve our system of equatins\n",
    "n2, o2, no = fsolve(eqn, (1, 1, 1,))\n",
    "\n",
    "# Calculating mole fractions\n",
    "ntot = n2 + o2 + no #total moles\n",
    "x_n2 = n2/ntot\n",
    "x_o2 = o2/ntot\n",
    "ppmv_no = no/ntot * 1e6\n",
    "\n",
    "print(\"x_N2 = %1.4f\" % x_n2)\n",
    "print(\"x_O2 = %1.4f\" % x_o2)\n",
    "print(\"ppmv_NO = %1.f ppmv\" % ppmv_no)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Checking answers with Cantera:\n",
    "\n",
    "Now let's check our answers with cantera:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Let's look at the equilibrium gas properties:\n",
      "\n",
      "  gri30:\n",
      "\n",
      "       temperature   1800 K\n",
      "          pressure   1.0132e+05 Pa\n",
      "           density   0.19532 kg/m^3\n",
      "  mean mol. weight   28.85 kg/kmol\n",
      "   phase of matter   gas\n",
      "\n",
      "                          1 kg             1 kmol     \n",
      "                     ---------------   ---------------\n",
      "          enthalpy        1.7309e+06        4.9937e+07  J\n",
      "   internal energy        1.2122e+06        3.4971e+07  J\n",
      "           entropy            8881.2        2.5622e+05  J/K\n",
      "    Gibbs function       -1.4255e+07       -4.1126e+08  J\n",
      " heat capacity c_p            1246.4             35960  J/K\n",
      " heat capacity c_v            958.24             27645  J/K\n",
      "\n",
      "                      mass frac. Y      mole frac. X     chem. pot. / RT\n",
      "                     ---------------   ---------------   ---------------\n",
      "                 O        3.0749e-05        5.5448e-05           -14.983\n",
      "                O2           0.23056           0.20788           -29.966\n",
      "                 N         1.584e-11        3.2626e-11            -13.41\n",
      "                NO          0.004319         0.0041526           -28.393\n",
      "               NO2        1.6606e-05        1.0413e-05           -43.376\n",
      "               N2O        3.6645e-07         2.402e-07           -41.803\n",
      "                N2           0.76507            0.7879           -26.819\n",
      "     [  +46 minor]                 0                 0  \n",
      "\n",
      "%=============== Now, picking out the main ones...\n",
      "The mole fraction of N2 is 0.7879 \n",
      "The mole fraction of O2 is 0.2079 \n",
      "The ppmv of NO is 4153 ppmv \n"
     ]
    }
   ],
   "source": [
    "# Cantera is already imported, so we just need to create a gas object\n",
    "gas1 = ct.Solution('gri30.yaml')\n",
    "#if we wanted to use a different mechanism\n",
    "# gas1 = ct.Solution('burkeh2o2.yaml')\n",
    "\n",
    "# Let's set what our gas is made of\n",
    "gas1.X = {'N2':0.79, 'O2':0.21}\n",
    "\n",
    "# Now set the temperature and pressure...\n",
    "gas1.TP = 1800, 101325\n",
    "\n",
    "# Finally, calculate the equilbrium gas concentrations based on the temp and pressure. Play around with different\n",
    "# temps and pressures and see how it changes the concentrations below.\n",
    "gas1.equilibrate('TP')\n",
    "print('Let\\'s look at the equilibrium gas properties:')\n",
    "gas1() #show the gas properties\n",
    "\n",
    "print('%=============== Now, picking out the main ones...')\n",
    "# We can list specific elements\n",
    "mf_n2 = gas1['N2'].X\n",
    "print(\"The mole fraction of N2 is %1.4f \" % mf_n2)\n",
    "\n",
    "mf_o2 = gas1['O2'].X\n",
    "print(\"The mole fraction of O2 is %1.4f \" % mf_o2)\n",
    "\n",
    "ppm_no = gas1['NO'].X * 1e6\n",
    "print(\"The ppmv of NO is %1.0f ppmv \" % ppm_no)\n",
    "\n",
    "# check out here for list of gas properties: \n",
    "# https://cantera.org/documentation/docs-2.4/sphinx/html/cython/thermo.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Equilibrium example with a hydrocarbon:\n",
    "\n",
    "Now let's try one with more possibilities\n",
    "\n",
    "$$ CH_4 + 2O_2 \\rightarrow Products\\left(CO_2, H_2O, CO, OH, O_2, H_2\\right) $$\n",
    "\n",
    "Let's assume this system reaches equilibrium at 2000K and 1 atm.\n",
    "\n",
    "Note though, that since we have more produts than in the last example, we'll need some more $K_p$ values. Here, we'll need them for $OH$, $H_2O$, and $CO_2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_CO2 = 0.3249\n",
      "x_H2O = 0.6591\n",
      "ppmv_OH = 2238 ppmv\n"
     ]
    }
   ],
   "source": [
    "# Set up the knowns\n",
    "n_ch4 = 1\n",
    "n_o2 = 2\n",
    "# For these Kp's we're using the values from the handout\n",
    "lnkp_oh = -0.4499\n",
    "lnkp_h2o = 8.1521\n",
    "lnkp_co2 = 6.6301\n",
    "\n",
    "T = 2000 #K\n",
    "P = 1 #atm\n",
    "\n",
    "# Define our equations to solve\n",
    "def eqn(vars):\n",
    "    co2, h2o, co, oh, o2, h2, nt  = vars\n",
    "    eq1 = n_ch4 - (co2 + co) #carbon balance\n",
    "    eq2 = 2*n_o2 - (2*co2 + h2o + oh + 2*o2 + co) #oxygen balance\n",
    "    eq3 = 4*n_ch4 - (2*h2o + oh + 2*h2) #hydrogen balance\n",
    "    eq4 = exp(lnkp_oh)  - (oh  / (o2**0.5 * h2**0.5))*1 #Kp_OH\n",
    "    eq5 = exp(lnkp_h2o) - (h2o / (h2 * o2**0.5))*(1/nt)**(-.5) #Kp_H2O\n",
    "    eq6 = exp(lnkp_co2) - (co2 / (co * o2**0.5))*(1/nt)**(-.5) #Kp_CO2\n",
    "    eq7 = nt - (co2 + h2o + co + oh + o2 + h2) #total moles\n",
    "    return [eq1, eq2, eq3, eq4, eq5, eq6, eq7]\n",
    "\n",
    "# Now call the fsolve function to solve our system of equatins\n",
    "# also note that you'll need a pretty good guess here...\n",
    "co2, h2o, co, oh, o2, h2, nt = fsolve(eqn, (.3, .7, .001, .001, .001, .001, 1))\n",
    "\n",
    "# Calculating mole fractions\n",
    "x_co2 = co2/nt\n",
    "x_h2o = h2o/nt\n",
    "ppmv_oh = oh/nt * 1e6\n",
    "\n",
    "print(\"x_CO2 = %1.4f\" % x_co2)\n",
    "print(\"x_H2O = %1.4f\" % x_h2o)\n",
    "print(\"ppmv_OH = %1.f ppmv\" % ppmv_oh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Checking answers with Cantera:\n",
    "\n",
    "Now let's check our answers with cantera:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "%===============\n",
      "The mole fraction of CO2 is 0.3250 \n",
      "The mole fraction of H2O is 0.6592 \n",
      "The ppmv of OH is 1981 ppmv \n"
     ]
    }
   ],
   "source": [
    "# Now we'll create a gas object\n",
    "gas2 = ct.Solution('gri30.yaml')\n",
    "\n",
    "#Let's set what our gas is made of\n",
    "gas2.X = {'CH4':1/3, 'O2':2/3} #remember these shoudl be in mole fractions\n",
    "\n",
    "#Now set the temperature and pressure...\n",
    "gas2.TP = 2000, 101325\n",
    "\n",
    "#Finally, calculate the equilbrium gas concentrations based on the temp and pressure. Play around with different\n",
    "# temps and pressures and see how it changes the concentrations below.\n",
    "gas2.equilibrate('TP')\n",
    "#gas1() #show the gas properties\n",
    "\n",
    "print(\"%===============\") #just a print command to separate the outputs\n",
    "# We can list specific elements\n",
    "mf_co2 = gas2['CO2'].X\n",
    "print(\"The mole fraction of CO2 is %1.4f \" % mf_co2)\n",
    "\n",
    "mf_h2o = gas2['H2O'].X\n",
    "print(\"The mole fraction of H2O is %1.4f \" % mf_h2o)\n",
    "\n",
    "ppm_oh = gas2['OH'].X * 1e6\n",
    "print(\"The ppmv of OH is %1.0f ppmv \" % ppm_oh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Same example but with importing our $K_p$ table:\n",
    "\n",
    "Another way we can try the previous example is by importing the relevant JANAF tabels from a *.csv file usign a toolbox calld \"pandas,\" which will make looking at, and utilizing tables much easier.\n",
    "\n",
    "Doing things this way will allow us to create code that will be more flexible and enable us to solve for a wider range of temperatures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>T(K)</th>\n",
       "      <th>CO+1/2O eq CO2</th>\n",
       "      <th>1/2H2 eq H</th>\n",
       "      <th>H2+1/2O eq H2O</th>\n",
       "      <th>1/2N2 eq N</th>\n",
       "      <th>1/2O2+1/2N2 eq NO</th>\n",
       "      <th>O2+1/2N2 eq NO2</th>\n",
       "      <th>1/2O2 eq O</th>\n",
       "      <th>1/2H2+1/2O2 eq OH</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>200.00</td>\n",
       "      <td>159.6811</td>\n",
       "      <td>-125.0847</td>\n",
       "      <td>139.9748</td>\n",
       "      <td>-277.2662</td>\n",
       "      <td>-53.3900</td>\n",
       "      <td>-27.9571</td>\n",
       "      <td>-142.7264</td>\n",
       "      <td>-20.5011</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>298.15</td>\n",
       "      <td>103.7531</td>\n",
       "      <td>-82.0006</td>\n",
       "      <td>92.2077</td>\n",
       "      <td>-183.7606</td>\n",
       "      <td>-35.3306</td>\n",
       "      <td>-21.1037</td>\n",
       "      <td>-93.4813</td>\n",
       "      <td>-13.1346</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>400.00</td>\n",
       "      <td>74.6618</td>\n",
       "      <td>-59.5794</td>\n",
       "      <td>67.3216</td>\n",
       "      <td>-135.1804</td>\n",
       "      <td>-25.9536</td>\n",
       "      <td>-17.6196</td>\n",
       "      <td>-67.8546</td>\n",
       "      <td>-9.3031</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>600.00</td>\n",
       "      <td>46.2378</td>\n",
       "      <td>-37.6127</td>\n",
       "      <td>42.8974</td>\n",
       "      <td>-87.6896</td>\n",
       "      <td>-16.7997</td>\n",
       "      <td>-14.2742</td>\n",
       "      <td>-42.7592</td>\n",
       "      <td>-5.5674</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>800.00</td>\n",
       "      <td>32.0310</td>\n",
       "      <td>-26.5668</td>\n",
       "      <td>30.5926</td>\n",
       "      <td>-63.8855</td>\n",
       "      <td>-12.2212</td>\n",
       "      <td>-12.6124</td>\n",
       "      <td>-30.1586</td>\n",
       "      <td>-3.7106</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1000.00</td>\n",
       "      <td>23.5225</td>\n",
       "      <td>-19.9031</td>\n",
       "      <td>23.1622</td>\n",
       "      <td>-49.5719</td>\n",
       "      <td>-9.4731</td>\n",
       "      <td>-11.6103</td>\n",
       "      <td>-22.5729</td>\n",
       "      <td>-2.6070</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>1200.00</td>\n",
       "      <td>17.8648</td>\n",
       "      <td>-15.4375</td>\n",
       "      <td>18.1826</td>\n",
       "      <td>-40.0116</td>\n",
       "      <td>-7.6402</td>\n",
       "      <td>-10.9354</td>\n",
       "      <td>-17.5017</td>\n",
       "      <td>-1.8784</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>1400.00</td>\n",
       "      <td>13.8360</td>\n",
       "      <td>-12.2327</td>\n",
       "      <td>14.6115</td>\n",
       "      <td>-33.1715</td>\n",
       "      <td>-6.3305</td>\n",
       "      <td>-10.4473</td>\n",
       "      <td>-13.8705</td>\n",
       "      <td>-1.3634</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>1600.00</td>\n",
       "      <td>10.8245</td>\n",
       "      <td>-9.8185</td>\n",
       "      <td>11.9246</td>\n",
       "      <td>-28.0337</td>\n",
       "      <td>-5.3480</td>\n",
       "      <td>-10.0761</td>\n",
       "      <td>-11.1413</td>\n",
       "      <td>-0.9802</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>1800.00</td>\n",
       "      <td>8.4906</td>\n",
       "      <td>-7.9333</td>\n",
       "      <td>9.8304</td>\n",
       "      <td>-24.0323</td>\n",
       "      <td>-4.5838</td>\n",
       "      <td>-9.7835</td>\n",
       "      <td>-9.0144</td>\n",
       "      <td>-0.6848</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>2000.00</td>\n",
       "      <td>6.6301</td>\n",
       "      <td>-6.4198</td>\n",
       "      <td>8.1521</td>\n",
       "      <td>-20.8271</td>\n",
       "      <td>-3.9724</td>\n",
       "      <td>-9.5461</td>\n",
       "      <td>-7.3100</td>\n",
       "      <td>-0.4499</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>2200.00</td>\n",
       "      <td>5.1140</td>\n",
       "      <td>-5.1776</td>\n",
       "      <td>6.7774</td>\n",
       "      <td>-18.2015</td>\n",
       "      <td>-3.4723</td>\n",
       "      <td>-9.3492</td>\n",
       "      <td>-5.9135</td>\n",
       "      <td>-0.2593</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>2400.00</td>\n",
       "      <td>3.8550</td>\n",
       "      <td>-4.1394</td>\n",
       "      <td>5.6310</td>\n",
       "      <td>-16.0113</td>\n",
       "      <td>-3.0560</td>\n",
       "      <td>-9.1829</td>\n",
       "      <td>-4.7482</td>\n",
       "      <td>-0.1015</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>2600.00</td>\n",
       "      <td>2.7937</td>\n",
       "      <td>-3.2584</td>\n",
       "      <td>4.6604</td>\n",
       "      <td>-14.1558</td>\n",
       "      <td>-2.7042</td>\n",
       "      <td>-9.0406</td>\n",
       "      <td>-3.7608</td>\n",
       "      <td>0.0309</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>2800.00</td>\n",
       "      <td>1.8874</td>\n",
       "      <td>-2.5015</td>\n",
       "      <td>3.8283</td>\n",
       "      <td>-12.5638</td>\n",
       "      <td>-2.4030</td>\n",
       "      <td>-8.9172</td>\n",
       "      <td>-2.9137</td>\n",
       "      <td>0.1435</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>3000.00</td>\n",
       "      <td>1.1050</td>\n",
       "      <td>-1.8441</td>\n",
       "      <td>3.1066</td>\n",
       "      <td>-11.1827</td>\n",
       "      <td>-2.1423</td>\n",
       "      <td>-8.8088</td>\n",
       "      <td>-2.1789</td>\n",
       "      <td>0.2402</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>3200.00</td>\n",
       "      <td>0.4228</td>\n",
       "      <td>-1.2678</td>\n",
       "      <td>2.4753</td>\n",
       "      <td>-9.9732</td>\n",
       "      <td>-1.9148</td>\n",
       "      <td>-8.7125</td>\n",
       "      <td>-1.5356</td>\n",
       "      <td>0.3242</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>3400.00</td>\n",
       "      <td>-0.1768</td>\n",
       "      <td>-0.7585</td>\n",
       "      <td>1.9181</td>\n",
       "      <td>-8.9046</td>\n",
       "      <td>-1.7147</td>\n",
       "      <td>-8.6266</td>\n",
       "      <td>-0.9675</td>\n",
       "      <td>0.3977</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>3600.00</td>\n",
       "      <td>-0.7080</td>\n",
       "      <td>-0.3049</td>\n",
       "      <td>1.4228</td>\n",
       "      <td>-7.9538</td>\n",
       "      <td>-1.5372</td>\n",
       "      <td>-8.5490</td>\n",
       "      <td>-0.4624</td>\n",
       "      <td>0.4621</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>3800.00</td>\n",
       "      <td>-1.1819</td>\n",
       "      <td>0.1015</td>\n",
       "      <td>0.9793</td>\n",
       "      <td>-7.1023</td>\n",
       "      <td>-1.3788</td>\n",
       "      <td>-8.4786</td>\n",
       "      <td>-0.0101</td>\n",
       "      <td>0.5192</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>4000.00</td>\n",
       "      <td>-1.6065</td>\n",
       "      <td>0.4677</td>\n",
       "      <td>0.5800</td>\n",
       "      <td>-6.3349</td>\n",
       "      <td>-1.2367</td>\n",
       "      <td>-8.4139</td>\n",
       "      <td>0.3972</td>\n",
       "      <td>0.5701</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       T(K)  CO+1/2O eq CO2  1/2H2 eq H  H2+1/2O eq H2O  1/2N2 eq N  \\\n",
       "0    200.00        159.6811   -125.0847        139.9748   -277.2662   \n",
       "1    298.15        103.7531    -82.0006         92.2077   -183.7606   \n",
       "2    400.00         74.6618    -59.5794         67.3216   -135.1804   \n",
       "3    600.00         46.2378    -37.6127         42.8974    -87.6896   \n",
       "4    800.00         32.0310    -26.5668         30.5926    -63.8855   \n",
       "5   1000.00         23.5225    -19.9031         23.1622    -49.5719   \n",
       "6   1200.00         17.8648    -15.4375         18.1826    -40.0116   \n",
       "7   1400.00         13.8360    -12.2327         14.6115    -33.1715   \n",
       "8   1600.00         10.8245     -9.8185         11.9246    -28.0337   \n",
       "9   1800.00          8.4906     -7.9333          9.8304    -24.0323   \n",
       "10  2000.00          6.6301     -6.4198          8.1521    -20.8271   \n",
       "11  2200.00          5.1140     -5.1776          6.7774    -18.2015   \n",
       "12  2400.00          3.8550     -4.1394          5.6310    -16.0113   \n",
       "13  2600.00          2.7937     -3.2584          4.6604    -14.1558   \n",
       "14  2800.00          1.8874     -2.5015          3.8283    -12.5638   \n",
       "15  3000.00          1.1050     -1.8441          3.1066    -11.1827   \n",
       "16  3200.00          0.4228     -1.2678          2.4753     -9.9732   \n",
       "17  3400.00         -0.1768     -0.7585          1.9181     -8.9046   \n",
       "18  3600.00         -0.7080     -0.3049          1.4228     -7.9538   \n",
       "19  3800.00         -1.1819      0.1015          0.9793     -7.1023   \n",
       "20  4000.00         -1.6065      0.4677          0.5800     -6.3349   \n",
       "\n",
       "    1/2O2+1/2N2 eq NO  O2+1/2N2 eq NO2  1/2O2 eq O  1/2H2+1/2O2 eq OH  \n",
       "0            -53.3900         -27.9571   -142.7264           -20.5011  \n",
       "1            -35.3306         -21.1037    -93.4813           -13.1346  \n",
       "2            -25.9536         -17.6196    -67.8546            -9.3031  \n",
       "3            -16.7997         -14.2742    -42.7592            -5.5674  \n",
       "4            -12.2212         -12.6124    -30.1586            -3.7106  \n",
       "5             -9.4731         -11.6103    -22.5729            -2.6070  \n",
       "6             -7.6402         -10.9354    -17.5017            -1.8784  \n",
       "7             -6.3305         -10.4473    -13.8705            -1.3634  \n",
       "8             -5.3480         -10.0761    -11.1413            -0.9802  \n",
       "9             -4.5838          -9.7835     -9.0144            -0.6848  \n",
       "10            -3.9724          -9.5461     -7.3100            -0.4499  \n",
       "11            -3.4723          -9.3492     -5.9135            -0.2593  \n",
       "12            -3.0560          -9.1829     -4.7482            -0.1015  \n",
       "13            -2.7042          -9.0406     -3.7608             0.0309  \n",
       "14            -2.4030          -8.9172     -2.9137             0.1435  \n",
       "15            -2.1423          -8.8088     -2.1789             0.2402  \n",
       "16            -1.9148          -8.7125     -1.5356             0.3242  \n",
       "17            -1.7147          -8.6266     -0.9675             0.3977  \n",
       "18            -1.5372          -8.5490     -0.4624             0.4621  \n",
       "19            -1.3788          -8.4786     -0.0101             0.5192  \n",
       "20            -1.2367          -8.4139      0.3972             0.5701  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The easiest way to use this is to have the *.csv file in the folder our jupyter notebook is in:\n",
    "data = pd.read_csv(\"lnkp_v_temp.csv\")\n",
    "data #show the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABDZElEQVR4nO3de3xU9Z3/8deZXCYXkpEAyWQkYFS8YBAFBEEr94sV0LqteC1u/bXeawreqHXB7i6pdhdtyxZb2xUrS+m6itXWKmmVIAUFwSgBC4IRAyQGIUwSyI3M9/fHkIHJBZIwkzOTvJ+Px3nMmXPOTD5fTtt593u+53ssY4xBREREJII47C5AREREpDkFFBEREYk4CigiIiIScRRQREREJOIooIiIiEjEUUARERGRiKOAIiIiIhFHAUVEREQiTqzdBXSGz+dj3759pKSkYFmW3eWIiIhIOxhjqKqqwuPx4HCcvI8kKgPKvn37yMrKsrsMERER6YSSkhL69+9/0mOiMqCkpKQA/gampqbaXI2IiIi0R2VlJVlZWYHf8ZOJyoDSdFknNTVVAUVERCTKtGd4hgbJioiISMRRQBEREZGIo4AiIiIiEUcBRURERCKOAoqIiIhEHAUUERERiTgKKCIiIhJxFFBEREQk4kTlRG3h0ugzbCg+SHlVLekpCYzMTiPGoWf9iIiIdLUO96CsWbOGGTNm4PF4sCyLV199tcUxn3zyCTNnzsTlcpGSksLll1/OF198EdhfV1fH/fffT9++fUlOTmbmzJns2bPntBpyut4sKuXKJ9/mpufe44EVhdz03Htc+eTbvFlUamtdIiIiPVGHA8rhw4cZOnQoixcvbnX/rl27uPLKK7ngggtYvXo1H330EY8//jgJCQmBY3Jzc1m5ciUrVqxg7dq1VFdXM336dBobGzvfktPwZlEpdy/bTKm3Nmh7mbeWu5dtVkgRERHpYpYxxnT6w5bFypUrue666wLbbrzxRuLi4njxxRdb/YzX66Vfv368+OKLzJo1Czj+dOI33niDqVOnnvLvVlZW4nK58Hq9p/0snkaf4con324RTppYgNuVwNpHJuhyj4iIyGnoyO93SAfJ+nw+/vznP3PeeecxdepU0tPTGTVqVNBloE2bNtHQ0MCUKVMC2zweDzk5Oaxbt67V762rq6OysjJoCZUNxQcD4cTNAR6J/T2Pxi4P7DdAqbeWDcUHQ/Y3RURE5ORCGlDKy8uprq7mJz/5CdOmTWPVqlV84xvf4Prrr6egoACAsrIy4uPj6d27d9BnMzIyKCsra/V78/LycLlcgSUrKyt0NVcd7zlJtmq5O/Z1bo7520mPExERkfAKeQ8KwLXXXssPfvADLrnkEh599FGmT5/Os88+e9LPGmPafPzyvHnz8Hq9gaWkpCRkNaenHB8b86Xxh6ZUq4ZEats8TkRERMIrpAGlb9++xMbGMnjw4KDtF154YeAuHrfbTX19PRUVFUHHlJeXk5GR0er3Op1OUlNTg5ZQGZmdRqYrAQuoJpHDxglAunUI8I9ByXT5bzkWERGRrhHSgBIfH89ll13G9u3bg7bv2LGDgQMHAjB8+HDi4uLIz88P7C8tLaWoqIgxY8aEspx2iXFYzJ/hD1QWVqAXJYMKmvpz5s8YrAGyIiIiXajDE7VVV1ezc+fOwPvi4mIKCwtJS0tjwIABPPTQQ8yaNYurrrqK8ePH8+abb/L666+zevVqAFwuF3fccQdz586lT58+pKWl8eCDDzJkyBAmTZoUsoZ1xLScTJbcOownXt/G/pozOJsyMqwK3KkJzJ8xmGk5mbbUJSIi0lN1+Dbj1atXM378+BbbZ8+ezdKlSwH47//+b/Ly8tizZw/nn38+TzzxBNdee23g2NraWh566CGWL19OTU0NEydO5Je//GW7B7+G8jbjEzX6DBW/u5W+n/+Jz4f/kKxrHlbPiYiISIh05Pf7tOZBsUu4AgoAbz0G6xfD6Ptg6r+H9rtFRER6MNvmQekWUtz+1+ov7a1DRESkB1NAaS7l2HiTqtbnZBEREZHwU0BprqkHpUrP3xEREbGLAkpzvZoCinpQRERE7KKA0lzKscni6quhrsreWkRERHooBZTmnCkQn+JfVy+KiIiILRRQWpOiyzwiIiJ2UkBpjQKKiIiIrRRQWhO41Vh38oiIiNhBAaU1TQNl1YMiIiJiCwWU1qgHRURExFYKKK3RGBQRERFbKaC0pqkHpVoBRURExA4KKK05sQcl+h72LCIiEvUUUFrTNN19wxGoq7S3FhERkR5IAaU18UngdPnXNQ5FRESkyymgtEVPNRYREbGNAkpbAgHlS3vrEBER6YEUUNqiuVBERERso4DSFs2FIiIiYhsFlLaoB0VERMQ2Ciht0fN4REREbKOA0hb1oIiIiNhGAaUtTWNQqr/UbLIiIiJdTAGlLU2zyR6thdpDtpYiIiLS0yigtCUuARJ7+9c1DkVERKRLKaCcjMahiIiI2EIB5WR66U4eEREROyignIx6UERERGyhgHIyeh6PiIiILTocUNasWcOMGTPweDxYlsWrr77a5rF33nknlmXxzDPPBG2vq6vj/vvvp2/fviQnJzNz5kz27NnT0VLCTz0oIiIituhwQDl8+DBDhw5l8eLFJz3u1Vdf5f3338fj8bTYl5uby8qVK1mxYgVr166lurqa6dOn09jY2NFywkvP4xEREbFFbEc/cPXVV3P11Vef9Ji9e/dy33338dZbb3HNNdcE7fN6vfz2t7/lxRdfZNKkSQAsW7aMrKws/vrXvzJ16tSOlhQ+gR4UBRQREZGuFPIxKD6fj9tuu42HHnqIiy66qMX+TZs20dDQwJQpUwLbPB4POTk5rFu3rtXvrKuro7KyMmjpEoHn8ZRqNlkREZEuFPKA8uSTTxIbG8v3v//9VveXlZURHx9P7969g7ZnZGRQVtZ6T0VeXh4ulyuwZGVlhbrs1jXdZuxrgCMHu+ZvioiISGgDyqZNm/jZz37G0qVLsSyrQ581xrT5mXnz5uH1egNLSUlJKMo9tVgnJPXxr1frMo+IiEhXCWlAeffddykvL2fAgAHExsYSGxvL7t27mTt3LmeddRYAbreb+vp6Kioqgj5bXl5ORkZGq9/rdDpJTU0NWrqM7uQRERHpciENKLfddhsff/wxhYWFgcXj8fDQQw/x1ltvATB8+HDi4uLIz88PfK60tJSioiLGjBkTynJCQ3fyiIiIdLkO38VTXV3Nzp07A++Li4spLCwkLS2NAQMG0KdPn6Dj4+LicLvdnH/++QC4XC7uuOMO5s6dS58+fUhLS+PBBx9kyJAhgbt6IkogoKgHRUREpKt0OKB88MEHjB8/PvB+zpw5AMyePZulS5e26zuefvppYmNjueGGG6ipqWHixIksXbqUmJiYjpYTfr3UgyIiItLVLGOi7/7ZyspKXC4XXq83/ONRNjwHbzwIF0yHG/8nvH9LRESkG+vI77eexXMqmqxNRESkyymgnIoCioiISJdTQDmVpkGy1WXg89lbi4iISA+hgHIqvdIBC3xH4cgBu6sRERHpERRQTiUmDpL7+td1q7GIiEiXUEBpj8Blni/trUNERKSHUEBpD013LyIi0qUUUNpD092LiIh0KQWU9lAPioiISJdSQGmPXseesqweFBERkS6hgNIe6kERERHpUgoo7REYg6K7eERERLqCAkp7NPWgVH8JvkZ7axEREekBFFDaI7kfWA4wjXD4K7urERER6fYUUNojJhaS0/3rGociIiISdgoo7ZWiO3lERES6igJKe+lOHhERkS6jgNJeeh6PiIhIl1FAaS/1oIiIiHQZBZT20vN4REREuowCSnupB0VERKTLKKC0l57HIyIi0mUUUNorMJtsOTQetbcWERGRbk4Bpb2S+4IVAxg4vN/uakRERLo1BZT2csSccJlH41BERETCSQGlI3Qnj4iISJdQQOkI3ckjIiLSJRRQOkLP4xEREekSCigdEbiTRwFFREQknBRQOkJjUERERLqEAkpHaAyKiIhIl+hwQFmzZg0zZszA4/FgWRavvvpqYF9DQwOPPPIIQ4YMITk5GY/Hw7e//W327dsX9B11dXXcf//99O3bl+TkZGbOnMmePXtOuzFhpx4UERGRLtHhgHL48GGGDh3K4sWLW+w7cuQImzdv5vHHH2fz5s288sor7Nixg5kzZwYdl5uby8qVK1mxYgVr166lurqa6dOn09jY2PmWdIWmHpTD+6Gxwd5aREREujHLGGM6/WHLYuXKlVx33XVtHrNx40ZGjhzJ7t27GTBgAF6vl379+vHiiy8ya9YsAPbt20dWVhZvvPEGU6dOPeXfraysxOVy4fV6SU1N7Wz5Hefzwb/1A99R+MFWcPXvur8tIiIS5Try+x32MSherxfLsjjjjDMA2LRpEw0NDUyZMiVwjMfjIScnh3Xr1rX6HXV1dVRWVgYttnA4oFfTZZ4v7alBRESkBwhrQKmtreXRRx/l5ptvDiSlsrIy4uPj6d27d9CxGRkZlJW1PrYjLy8Pl8sVWLKyssJZ9skFxqFooKyIiEi4hC2gNDQ0cOONN+Lz+fjlL395yuONMViW1eq+efPm4fV6A0tJSUmoy20/BRQREZGwC0tAaWho4IYbbqC4uJj8/Pyg60xut5v6+noqKiqCPlNeXk5GRkar3+d0OklNTQ1abBO41Vh38oiIiIRLyANKUzj59NNP+etf/0qfPn2C9g8fPpy4uDjy8/MD20pLSykqKmLMmDGhLif0dKuxiIhI2MV29APV1dXs3Lkz8L64uJjCwkLS0tLweDx885vfZPPmzfzpT3+isbExMK4kLS2N+Ph4XC4Xd9xxB3PnzqVPnz6kpaXx4IMPMmTIECZNmhS6loWLLvGIiIiEXYcDygcffMD48eMD7+fMmQPA7NmzWbBgAa+99hoAl1xySdDn3nnnHcaNGwfA008/TWxsLDfccAM1NTVMnDiRpUuXEhMT08lmdKGmgFKtu3hERETC5bTmQbGLbfOgAHy5FZaMgaQ+8PBnXfu3RUREolhEzYPS7TQNkj1yAI7W2VuLiIhIN6WA0lGJvSEm3r+uyzwiIiJhoYDSUZalO3lERETCTAGlM3rpTh4REZFwUkDpjBQ9j0dERCScFFA6IzCbrHpQREREwkEBpTM0BkVERCSsFFA6Qz0oIiIiYaWA0hnqQREREQkrBZTO0PN4REREwkoBpTOaAkrtIWiotbUUERGR7kgBpTMSzoDYBP96tS7ziIiIhJoCSmdoNlkREZGwUkDpLN3JIyIiEjYKKJ2lHhQREZGwUUDprF4KKCIiIuGigNJZ6kEREREJGwWUztIYFBERkbBRQOks9aCIiIiEjQJKZwV6UBRQREREQk0BpbOaelDqvFB/xN5aREREuhkFlM5ypkBckn9ds8mKiIiElAJKZ1kW5lgvyrubt7B+1wEafcbmokRERLqHWLsLiFZvFpWSUZHApcAf3t7In3yxZLoSmD9jMNNyMu0uT0REJKqpB6UT3iwq5e5lmyk56gIgw6oAoMxby93LNvNmkW49FhEROR0KKB3U6DM88fo2DPCl6Q1A+rGA0nSB54nXt+lyj4iIyGlQQOmgDcUHKfXWAscDSlMPCvhDSqm3lg3FB+0oT0REpFtQQOmg8qra4+vmDAAyOHTS40RERKRjFFA6KD0lIbBeYtIBOMex76THiYiISMcooHTQyOw0Ml0JWMAnZgCNxiLdOkQ6/ss8FpDpSmBkdpqtdYqIiEQzBZQOinFYzJ8xGIA6nHxq+gNwseMzrGPHzJ8xmBiH1cY3iIiIyKl0OKCsWbOGGTNm4PF4sCyLV199NWi/MYYFCxbg8XhITExk3LhxbN26NeiYuro67r//fvr27UtycjIzZ85kz549p9WQrjQtJ5Mltw7D7Upgiy8bgCGOYtyuBJbcOkzzoIiIiJymDgeUw4cPM3ToUBYvXtzq/qeeeopFixaxePFiNm7ciNvtZvLkyVRVVQWOyc3NZeXKlaxYsYK1a9dSXV3N9OnTaWxs7HxLuti0nEzWPjKBEaMnAPDtgRWsfWSCwomIiEgIWMaYTk/YYVkWK1eu5LrrrgP8vScej4fc3FweeeQRwN9bkpGRwZNPPsmdd96J1+ulX79+vPjii8yaNQuAffv2kZWVxRtvvMHUqVNP+XcrKytxuVx4vV5SU1M7W35olGyA306GXhkwdztYurQjIiLSmo78fod0DEpxcTFlZWVMmTIlsM3pdDJ27FjWrVsHwKZNm2hoaAg6xuPxkJOTEzimubq6OiorK4OWiJGRA1YMVH8JVZpBVkREJBRCGlDKyvxP9c3IyAjanpGREdhXVlZGfHw8vXv3bvOY5vLy8nC5XIElKysrlGWfnvgk6HeBf31foa2liIiIdBdhuYvHanaZwxjTYltzJztm3rx5eL3ewFJSUhKyWkPCc4n/tbTQzipERES6jZAGFLfbDdCiJ6S8vDzQq+J2u6mvr6eioqLNY5pzOp2kpqYGLRHFc6n/dd+H9tYhIiLSTYQ0oGRnZ+N2u8nPzw9sq6+vp6CggDFjxgAwfPhw4uLigo4pLS2lqKgocEzUybzE/7qvEDo/5lhERESOie3oB6qrq9m5c2fgfXFxMYWFhaSlpTFgwAByc3NZuHAhgwYNYtCgQSxcuJCkpCRuvvlmAFwuF3fccQdz586lT58+pKWl8eCDDzJkyBAmTZoUupZ1JfexgbKHy/0DZVM9dlckIiIS1TocUD744APGjx8feD9nzhwAZs+ezdKlS3n44YepqanhnnvuoaKiglGjRrFq1SpSUlICn3n66aeJjY3lhhtuoKamhokTJ7J06VJiYmJC0CQbxCX6B8qWb/Vf5lFAEREROS2nNQ+KXSJqHpQmr94LhcvgqodhwmN2VyMiIhJxbJsHpUfTnTwiIiIho4ASKoGBsh9qoKyIiMhpUkAJlcBA2f1Quc/uakRERKKaAkqoxCVC+oX+dV3mEREROS0KKKF04nwoIiIi0mkKKKHUNFBWM8qKiIicFgWUUGqa8r60UANlRUREToMCSihlXKSBsiIiIiGggBJKJw6U1WUeERGRTlNACTVN2CYiInLaFFBCTXfyiIiInDYFlFBrGiirGWVFREQ6TQEl1DIuAkcsHPkKKvfaXY2IiEhUUkAJtbhE6Nc0ULbQ1lJERESilQJKOHiG+l91J4+IiEinKKCEw4kTtomIiEiHKaCEQ2bTQNlCDZQVERHpBAWUcDhxoKx3j93ViIiIRB0FlHCISzg+UFaXeURERDpMASVcAk82LrSzChERkaikgBIumvJeRESk0xRQwiVTM8qKiIh0lgJKuAQGyh7QQFkREZEOUkAJl7gESNdAWRERkc5QQAmnwJONNaOsiIhIRyighJPnhAnbREREpN0UUMLpxDt5NFBWRESk3RRQwin9xIGyJXZXIyIiEjUUUMIpLgHSB/vXdZlHRESk3RRQwk0TtomIiHRYyAPK0aNH+dGPfkR2djaJiYmcffbZ/PjHP8bn8wWOMcawYMECPB4PiYmJjBs3jq1bt4a6lMigO3lEREQ6LOQB5cknn+TZZ59l8eLFfPLJJzz11FP89Kc/5Re/+EXgmKeeeopFixaxePFiNm7ciNvtZvLkyVRVVYW6HPud+EweDZQVERFpl5AHlPXr13PttddyzTXXcNZZZ/HNb36TKVOm8MEHHwD+3pNnnnmGxx57jOuvv56cnBxeeOEFjhw5wvLly0Ndjv0ycsARBzUHNVBWRESknUIeUK688kr+9re/sWPHDgA++ugj1q5dy9e//nUAiouLKSsrY8qUKYHPOJ1Oxo4dy7p161r9zrq6OiorK4OWqBHrPD6jrC7ziIiItEvIA8ojjzzCTTfdxAUXXEBcXByXXnopubm53HTTTQCUlZUBkJGREfS5jIyMwL7m8vLycLlcgSUrKyvUZYfXiZd5RERE5JRCHlD+8Ic/sGzZMpYvX87mzZt54YUX+I//+A9eeOGFoOMsywp6b4xpsa3JvHnz8Hq9gaWkJMoulTTNKKs7eURERNolNtRf+NBDD/Hoo49y4403AjBkyBB2795NXl4es2fPxu12A/6elMzMzMDnysvLW/SqNHE6nTidzlCX2nUCd/IU+gfKthHERERExC/kPShHjhzB4Qj+2piYmMBtxtnZ2bjdbvLz8wP76+vrKSgoYMyYMaEuJzJkXHR8oOyhL+yuRkREJOKFvAdlxowZ/Pu//zsDBgzgoosu4sMPP2TRokV85zvfAfyXdnJzc1m4cCGDBg1i0KBBLFy4kKSkJG6++eZQlxMZYp2QMRhKP/Jf5uk90O6KREREIlrIA8ovfvELHn/8ce655x7Ky8vxeDzceeed/Mu//EvgmIcffpiamhruueceKioqGDVqFKtWrSIlJSXU5USOzEv8AWVfIQy+1u5qREREIpplTPTNHlZZWYnL5cLr9ZKammp3Oe3zwX/Dn34AZ4+Hb79qdzUiIiJdriO/33oWT1c5NlC2Yc+H/PHDPazfdYBGX9RlQxERkS4R8ks80rq3vurDBGKJqz/ET//3r+wx/ch0JTB/xmCm5WSe+gtERER6EPWgdIE3i0q56/dF/MPXH4Ah1mcAlHlruXvZZt4sKrWzPBERkYijgBJmjT7DE69vwwBbfNkADHEUA9B0geeJ17fpco+IiMgJFFDCbEPxQUq9tQBsMWcDcKm1M7DfAKXeWjYUH7SjPBERkYikgBJm5VW1gfW/+3IAGOX4hP7W/jaPExER6ekUUMIsPSUhsP6FyWBt40U4LMOsmHfaPE5ERKSnU0AJs5HZaWS6Emh6+s7yxokAzIpZTSxHsYBMVwIjs9PsKlFERCTiKKCEWYzDYv6MwQBYQL5vBPtNKunWISY5NgMwf8ZgYhx6gKCIiEgTBZQuMC0nkyW3DsPtSqCBWF5qHAfA7c7VLLl1mOZBERERaUYTtXWRaTmZTB7sZkPxQarL+mDyX+dyUwgeDY4VERFpTj0oXSjGYTH6nD5MvuJyrHMm+DduWmprTSIiIpFIAcUuI/7Z/1r4P3C03t5aREREIowCil3Omwa93HB4P/zjT3ZXIyIiElEUUOwSEwfDbvOvb3re3lpEREQijAKKnYZ9G7CgeA0c2GV3NSIiIhFDAcVOZwyAQZP96+pFERERCVBAsdvwpsGyy+Fonb21iIiIRAgFFLsNmgKpZ8KRA/DJ63ZXIyIiEhEUUOwWE3tsLArwgS7ziIiIgAJKZLj0NrAcsHst7N9udzUiIiK2U0CJBK4z/fOigGaWFRERQQElcpw4WLahxt5aREREbKaAEinOnQiuLKg9BNv+aHc1IiIitlJAiRSOGBg227+uwbIiItLDKaBEkktvBSsGSt6D8k/srkZERMQ2CiiRJDUTzr/av65eFBER6cEUUCLNiO/4Xz9aAfVH7K1FRETEJgookebs8dD7LKjzwtZX7K5GRETEFgookcbh0GBZERHp8RRQItGlt4IjFvZ+AGVb7K5GRESky4UloOzdu5dbb72VPn36kJSUxCWXXMKmTZsC+40xLFiwAI/HQ2JiIuPGjWPr1q3hKCU69UqHC6b719WLIiIiPVDIA0pFRQVXXHEFcXFx/OUvf2Hbtm3853/+J2eccUbgmKeeeopFixaxePFiNm7ciNvtZvLkyVRVVYW6nOg14tjMsh//L9RV21uLiIhIF7OMMSaUX/joo4/y97//nXfffbfV/cYYPB4Pubm5PPLIIwDU1dWRkZHBk08+yZ133nnKv1FZWYnL5cLr9ZKamhrK8iOHzweLh8PBz2DGz2H4bLsrEhEROS0d+f0OeQ/Ka6+9xogRI/jWt75Feno6l156Kc8991xgf3FxMWVlZUyZMiWwzel0MnbsWNatW9fqd9bV1VFZWRm0dHsOBwy/3b++SZd5RESkZwl5QPnss89YsmQJgwYN4q233uKuu+7i+9//Pr/73e8AKCsrAyAjIyPocxkZGYF9zeXl5eFyuQJLVlZWqMuOTJfcAjHxsO9D2FdodzUiIiJdJuQBxefzMWzYMBYuXMill17KnXfeyXe/+12WLFkSdJxlWUHvjTEttjWZN28eXq83sJSUlIS67MiU3BcunOlfVy+KiIj0ICEPKJmZmQwePDho24UXXsgXX3wBgNvtBmjRW1JeXt6iV6WJ0+kkNTU1aOkxAoNlX4LaHnBpS0REhDAElCuuuILt27cHbduxYwcDBw4EIDs7G7fbTX5+fmB/fX09BQUFjBkzJtTlRL+BV0Df86DhMGx5ye5qREREukTIA8oPfvAD3nvvPRYuXMjOnTtZvnw5v/71r7n33nsB/6Wd3NxcFi5cyMqVKykqKuL2228nKSmJm2++OdTlRD/LCh4sG9qbrkRERCJSyAPKZZddxsqVK/n9739PTk4O//qv/8ozzzzDLbfcEjjm4YcfJjc3l3vuuYcRI0awd+9eVq1aRUpKSqjL6R6G3gQxTv+ssns3212NiIhI2IV8HpSu0CPmQWnule/Bx3/wT4N/7X/ZXY2IiEiH2ToPioTJ8GODZYtegVqvvbWIiIiEmQJKtBhwOfS7ABqO+Ke/FxER6cYUUKKFZcGI7/jXP9BgWRER6d4UUKLJxbMgNhHKt0LJBrurERERCRsFlGiSeAbkXO9f18yyIiLSjSmgRJumwbJbV0JNhb21iIiIhIkCSrTpPwIycuBoLXy0wu5qREREwkIBJdqcOLOsBsuKiEg3pYASjS6+AeKS4Kvt8MV6u6sREREJOQWUaJTggpx/8q9/oMGyIiLS/SigRKsRxwbLbvsjHD5gby0iIiIhpoASrTzDIHMoNNbBR8vtrkZERCSkFFCilWUdv+V401INlhURkW5FASWaDfkmxPeCAzvh83ftrkZERCRkFFCimTMFhnzLv67BsiIi0o0ooES7psGyn7wO1fvtrUVERCREFFCiXeZQ/4BZXwMU/o/d1YiIiISEAkp3cKwXpfb9/+aPH5awftcBGn0aNCsiItEr1u4C5PTlO65gNEn0qtrN/770P/zdN4RMVwLzZwxmWk6m3eWJiIh0mHpQotybRaV8b8U/ePnoFQDcHPM3AMq8tdy9bDNvFpXaWZ6IiEinKKBEsUaf4YnXt2GA5Y0TAZji2EQ/DtF0geeJ17fpco+IiEQdBZQotqH4IKXeWgC2mwFs8g0izmrk27GrADBAqbeWDcUHbaxSRESk4xRQolh5VW3Q++ePTgPgnpg/coVjS5vHiYiIRDoFlCiWnpIQ9P5Pvsv5v8ariLEMi+N+QZb1ZavHiYiIRDoFlCg2MjuNTFcCVmCLxWMN36HQdza9rWqei1tEdqr/OBERkWiigBLFYhwW82cMBgiElDriuav+B+w3Li5wlLA8/XfEWG1/h4iISCRSQIly03IyWXLrMNyu45dxyujDj+IfxmfFkbnnTVi7yMYKRUREOk4TtXUD03IymTzYzYbig5RX1ZKeksDI7K/j2JwEf/oB/O1fIWMInDfF7lJFRETaRQGlm4hxWIw+p0/wxhHfgdKPYNNSePn/wXffhr7n2lKfiIhIR+gST3d39VOQNQrqvLDiZqittLsiERGRU1JA6e5inXDD7yAlE77aDivvAp/P7qpEREROKuwBJS8vD8uyyM3NDWwzxrBgwQI8Hg+JiYmMGzeOrVu3hruUnivFDbOWQUw8bP8zrHnK7opEREROKqwBZePGjfz617/m4osvDtr+1FNPsWjRIhYvXszGjRtxu91MnjyZqqqqcJbTs/UfAdOf9q+vzoN//NneekRERE4ibAGlurqaW265heeee47evXsHthtjeOaZZ3jssce4/vrrycnJ4YUXXuDIkSMsX748XOUIwKW3wsg7/euvfA/K/2FvPSIiIm0IW0C59957ueaaa5g0aVLQ9uLiYsrKypgy5fgtr06nk7Fjx7Ju3bpWv6uuro7KysqgRTpp6r/DwCuhvto/aLbmkN0ViYiItBCWgLJixQo2b95MXl5ei31lZWUAZGRkBG3PyMgI7GsuLy8Pl8sVWLKyskJfdE8REwc3vACuLDi4y3/7sa/R7qpERESChDyglJSU8MADD7Bs2TISEtp+SJ1lBc+/boxpsa3JvHnz8Hq9gaWkpCSkNfc4yX39g2ZjE2BnPrzz73ZXJCIiEiTkAWXTpk2Ul5czfPhwYmNjiY2NpaCggJ///OfExsYGek6a95aUl5e36FVp4nQ6SU1NDVrkNHkugZm/8K+/+5+wdaWt5YiIiJwo5AFl4sSJbNmyhcLCwsAyYsQIbrnlFgoLCzn77LNxu93k5+cHPlNfX09BQQFjxowJdTlyMhffAKPv86+/eg+UFdlbj4iIyDEhn+o+JSWFnJycoG3Jycn06dMnsD03N5eFCxcyaNAgBg0axMKFC0lKSuLmm28OdTlyKpOegC+3wmfv+AfNfm81JKXZXZWIiPRwtjyL5+GHH6ampoZ77rmHiooKRo0axapVq0hJSbGjnJ4tJha++d/w3Hio+Bz+75/hlpf920VERGxiGWOM3UV0VGVlJS6XC6/Xq/EoofLlVvjNJGg44r/sM1UDZ0VEJLQ68vutZ/GIX8ZFcN0S//r6xfDx/9pbj4iI9GgKKHLcRdfB1+b611+7H/YV2lmNiIj0YAooEmz8YzBoChythRW3QPV+uysSEZEeSAFFgjli4PrnIO0cqNwDL90OjQ12VyUiIj2MAoq0lHgG3PR7iE+B3Wvhrcdo9BnW7zrAHwv3sn7XARp9UTe2WkREoojuJZXW9Tsfrv81rLgJNvyKvM2x/Kb6isDuTFcC82cMZlpOpo1FiohId6UeFGnbBV/n08H3A/BQw6+4xNoZ2FXmreXuZZt5s6jUrupERKQbU0CRNjX6DLN3juWtxhE4raM8G/80/agAoOkCzxOvb9PlHhERCTkFFGnThuKD7KusZ07D3ezwnYnbquDZ+GdwUg/4Q0qpt5YNxQftLVRERLodBRRpU3lVLQCHSeR7DXOoNEkMd3zKq/H/wtnWvhbHiYiIhIoCirQpPSUhsP65yeT/1c/lK5PKhY4veD3+Mb7heLfFcSIiIqGggCJtGpmdRqYrAevY+w3mQq6uy2Nd42CSrTqejl/C4qTnGHmm09Y6RUSk+1FAkTbFOCzmzxgMEAgp++nNrQ0/ZFHDN2k0FtN97xDzmwn+hw2KiIiEiAKKnNS0nEyW3DoMt+v4ZRwfDl7qdTMfjHsBUjLhq+3w3ATYtBSi7+HYIiISgSxjou8XpSOPa5bQaPQZNhQfpLyqlvSUBEZmpxHjsODwV7DyLtiZ7z8w559g+jOQoPMiIiLBOvL7rYAip8/ng/W/gL/9GHxHoXc2fGspeC6xuzIREYkgHfn91iUeOX0OB1zxAPzzX8CVBRXF8NvJ8P6vdMlHREQ6RQFFQidrJNz1LlwwHRrr4S8Pwx9uhZoKuysTEZEoo4AioZXYG2Ytg6ufgph4+Mef4NmroGSj3ZWJiEgUUUCR0LMsGHUn3LHKPx7F+wU8Pw3+/jP/eBUREZFTUECR8PFcCneu8d/Z4zsK+f8Cy2/w3/kjIiJyEgooEl4JqfBPv4UZP4fYBP/tyM9eCZ+vtbsyERGJYAooEn6WBcNnw3ffhr7nQ1UpvDADCp4CXyPgn2dl/a4D/LFwL+t3HaDRp7t/RER6Ms2DIl2r/jC88TAULvO/z76Kdwb/Gz/861eUeo8/FTnTlcD8GYOZlpNpU6EiIhJqmgdFIld8Mlz3X/CNX0FcMhSvYcifpnNOVfBdPmXeWu5etpk3i0ptKlREROykgCL2GHojjd9dzQ7rLPpalfwu7ic8GPsHYvBf8mnq1nvi9W263CMi0gMpoIhtNlT1YUbNAn53dDIOy3Bf7B/5Q/y/crG1C/CHlFJvLRuKD9pbqIiIdDkFFLFNeVUtdcTzL0f/mbvrH6DSJDLCsYPXnI/zfNyTDLN2BI4TEZGeJdbuAqTnSk9JCKz/xTeKLfVn80DMy3wjZi3jYz5ifMxHvNuYQ+8jPwLOtK9QERHpcupBEduMzE4j05WAdez9HtOPh47exfj6/+T3R8fTYGL4WkwROatuhKXToXiNHj4oItJDhDyg5OXlcdlll5GSkkJ6ejrXXXcd27dvDzrGGMOCBQvweDwkJiYybtw4tm7dGupSJMLFOCzmzxgMEAgpACUmgx8e/S7j6xbxxdk3giMOPn/XP3fK81fDrrcVVEREurmQB5SCggLuvfde3nvvPfLz8zl69ChTpkzh8OHDgWOeeuopFi1axOLFi9m4cSNut5vJkydTVVUV6nIkwk3LyWTJrcNwuxKCtrtdCfzo1qkM+Pav4IFCGPk9iHHCF+vhxW/AbyfDp/kKKiIi3VTYJ2rbv38/6enpFBQUcNVVV2GMwePxkJubyyOPPAJAXV0dGRkZPPnkk9x5552n/E5N1Nb9NPoMG4oPUl5VS3pKAiOz04hxWMEHVZbCup/DB8/D0Rr/Ns+lcNXDcP7V/hlrRUQkYkXURG1erxeAtLQ0AIqLiykrK2PKlCmBY5xOJ2PHjmXdunWtfkddXR2VlZVBi3QvMQ6L0ef04dpLzmT0OX1ahhOA1EyYlge5H8OY70NcEuz7EFbcBL/6Gmx7TU9LFhHpJsIaUIwxzJkzhyuvvJKcnBwAysrKAMjIyAg6NiMjI7Cvuby8PFwuV2DJysoKZ9kS6Xqlw5R/hdwtcOUciO8FZVvgf2+DZ6+AolcCz/hpomf9iIhEl7DeZnzffffx8ccfs3ZtyyfXWs26440xLbY1mTdvHnPmzAm8r6ysVEgRSO4Lk+bDmPvhvSXw/rNQvg3+75/9DyW86kG46Hre/GQ/T7y+Tc/6ERGJImHrQbn//vt57bXXeOedd+jfv39gu9vtBmjRW1JeXt6iV6WJ0+kkNTU1aBEJSEqDCY/5e1TG/RASXPDVdnjluxx+ehj5y5+h3Hs46CN61o+ISGQLeUAxxnDffffxyiuv8Pbbb5OdnR20Pzs7G7fbTX5+fmBbfX09BQUFjBkzJtTlSE+SeAaMewRyi2DC45jE3iRX7+Y/45/l7fi53BDzDnEcBfSsHxGRSBfygHLvvfeybNkyli9fTkpKCmVlZZSVlVFT47/rwrIscnNzWbhwIStXrqSoqIjbb7+dpKQkbr755lCXIz1RQipc9SAbri0gr+EmvjKpDHSU81Tcc7zjnMMPYl/iHGuvnvUjIhLBQn6bcVvjSJ5//nluv/12wN/L8sQTT/CrX/2KiooKRo0axX/9138FBtKeim4zlvb4Y+FeHlhRSCK13BzzNnfG/ol061Bg/1bfQF5rHMPwa+5gypjL7CtURKSH6Mjvd9jnQQkHBRRpj/W7DnDTc+8F3jupZ6pjIzNj1jHW8TFx1gl3+gwYDTn/BBd9wz/4VkREQk4BRQT/rcVXPvk2Zd5amv+H/Ayq+HrMBv4p/n2Gma1YTUdYMXDOeMj5Jlxwjf9ykYiIhIQCisgxbxaVcveyzQBBIaXpQuSSW4cxbYDxz51S9H/+id+axCbAeVP9YWXQFIgLno6/SbtmwRUREQUUkRO9WVTa/nlQDuyCLf/nDytf7Ti+3ZkKF0yHId+E7LEQE9vx7xYR6eEUUESa6XAvhzH+2Wm3vOTvXancc3xfcj8YfB3vJY/npjcNhuDvCeqdUUgREQlQQBEJJZ8PSt7396psXQlHDgR27TF9eb1xNK81juETM4CmeGLhfyLz2kcm6HKPiMgxCigi4dLYAJ8VUL7+f0ja9Qa9rOOXdnb4zuS1xjG87buUT8wADA5+/93LGX1OHxsLFhGJHAooImH2x8K9PLxiAxMcHzIzZh0THIU4rYbA/oOmF+t9g+k/bBpDv3Yt9DkH2pgjSESkp1BAEQmz5nOspHCEKY4PuDrmfS53fBLUs+I/wAPZV8HZY/2vrv6IiPQ0CigiYXayOVZiOcrF1mdMSdrBnVklWCUboLEu+KC0c44HlrO+dtLJ4XQbs4h0FwooIl2gXXOs5GRCQ41/kG3xGvisAPZtBuML/rKMIf7Akn0VDBwTmCBOtzGLSHeigCLSRToVIGq9sHvd8cBSvjV4vxUDZw5jV6/hPP5xHzb5BlFH/PHdx151G7OIRBsFFJEudNqXYKr3w+drjgeWiuKg3XUmjg985/F330Ws913Ex+ZsfMToNmYRiToKKCLR7NAX7Hz/DT5e+zpXOIrIOOEJzAC1Jo5tZiBbfNlc8bVJnDv0Suh3QWB2247SGBcR6Sod+f3u3P+iiUj4nDGArRkzmNOQBRjOsfYxxrGVMY6tjHZs4wzrMMOsnQxz7IT1+bAeiE0Edw5kXgKeS8BzKfQ9/5ShRWNcRCRSqQdFJAI1v425iYWPs6wvGWIVk+MoZtaZX+Gq2Ab1VS2/pCm0eC49HlxOCC1Ng3yb/w+AxriISLioB0Ukyo3MTiPTldDiNmaDg2KTyecmk40pE7jjzgmAgYOf+Z/EXFoI+wqh9CN/aNmz0b80iU0E9xB8mUPZ+KGT86z+7DRn0kjMCX/DH1KeeH0bkwe7dblHRGyhHhSRCNXu25hb4/PBwV3HwkrhsfDyEdRXtzi0xsQHxrQUmWy2+LL5zHhoIPa0p+rX+BYROZEGyYp0EyEdI3JCaNn50bsc+PR9LrI+bznrLXDUONhtMkjyXEjmORdDv/Oh73n+JaF9/53T+BYRaU4BRaQbCUcvRNMYFwsfZ1ul5FjFDHH4l8HWblKsmrY/nJIJfQf5x7P0O//4eoo78LwhjW8RkdYooIjISZ1sqn4wuKlgZMp+npmYjOPADti/Hb76FKrL2v5SZyr0HYSv73ksKYqhsCadneZMvjDpQWNcLAjZHC66hCQSXTRIVkROKsZhMX/GYO5ethmL5mNcLL4kja9fOwlH816OmkP+oPLVDvhqO+zf4V+vKIa6Sti7CcfeTdwLNE1+W29i+Ny42WnOZKfxsMvn4bNKD5v/MZDLLuz8U551CUmke1MPikgPFrIf+aN1cGAXfLWDT7Z8wPatmzjX2sc51j4Srfq2PxefAr0HwhkDW3+NT26zbl1CEok+usQjIu0W6sskJ87hYuHDwwHOdezjXGsv51j7OMexj2yrjPRmM+S2Kqlvi+DS6BrIjS/to7AyhYZWOoF1CUkkcukSj4i0W4zDOq1biZsLnsPFwV76sdfXjwKGAicEiDmjiaksgYrdcGg3VHx+7PXY+1ovHPnKv+zddLxe4CXA57QoozclJp09ph8lph8lvnT/qzedDbv2M3pQeqfboUtIIvZSD4qIhNxpzeHSpOZQcGA59AVU7KaybCdxlSUnv3QE+KwYHCluSPX47zxq7TXVA3GJbdYf7ktI6qGRnkaXeETEduHqgfBfQlpPXyrJssrJsvbT/9hr03uPdYA4q7F9X5hwxgmhJRNfioen1lWy/UgKX5relJk0DpJCUzwJ1SUk9dBIT6SAIiIRIRw9BCe/RdofIM5MjaPgnsHEVJdB5T6oKg1+bVpvONKuv1lnYv1hhbRAaJk86hLOyj4XktOhVzok94PE3u26K6kremjUOyORSAFFRLq1kFxCMsY/zqVZaCku/pRduz7FbVWQYR2kn1XZ/sIcsf6g0rQ0BZem1+R+NCb1Y+bSHfyjMj5ofpgT23C6PTTqnZFIpYAiIt1eeC8hHX+SdDwNpFuHyOAgbqsCt3WQDKuC68+x6GsOwuH9cLjcH3Y6wGcsKujFV8bFV8bFAVID6/tx8Z0pl3Hhuef4e2iS+7Y6VqY1Gj8jkUwBRUR6BLsuIbXaw3G0Dg5/5Q8r1cdCy+H9LdZrvWXE1R4kxurg//TGJkJSGiSm+V+D1vtAYhqNCb353v99xo7qeA6ZFKpI5Hg0ia7xMwpA3ZMCiojIaQjJJaQ2rN91gFueW0dvquhreelreelDJf2OrffFSx+rkpH9jtLraIU/2DSe/I6ltjSYGA6RzCGTwkFSOGR6UWF68bWh53Gmp//xgBN47Q0JLoh1tvmdXdFDowDUfUVNQPnlL3/JT3/6U0pLS7nooot45pln+NrXvnbKzymgiEi4hetHssM9NMZAXRXUHIQjB4+/BtYPBNYPHfiSw4fK6U01SVZdp2skNtEfVBLP8L8m+F99CS5e2FzBnlonlSRTaZKpJAmvSQ6875V6Bu8+Oum0xs90hwAkrYuKgPKHP/yB2267jV/+8pdcccUV/OpXv+I3v/kN27ZtY8CAASf9rAKKiHSFcP2/7HD10Jw4fsZJPWdQTZpVxRlWNb2povex15uH9MITX9Mi4FBb2ayijms0FsaZSmxy76Bw0zzs4Ew5YUkFZwqN8SlM+MUmdlf6Wv3uUA0g1h1U9omKgDJq1CiGDRvGkiVLAtsuvPBCrrvuOvLy8k76WQUUEYl24fh/8Z0eP9PE5/M/9LH2kH/Qb82x12Pvt39ewnvbPsNlHSaVI81eD5NgNXSq7ubqTCzVJFJtEqkmkSqSqDq2Xm0SmTD0XM50pweCzfHXY0tTAGp2qarp3+fEf/MO/fu0Q3e4PBXO74/4gFJfX09SUhIvvfQS3/jGNwLbH3jgAQoLCykoKAg6vq6ujrq6492VlZWVZGVlKaCISFQLxw9BuMfPnHiHU3NO6knlCM/dcC6X9LOaBZxDx94f8l+yalpqK4+v11d1qq42OeLA2Qvi/UuVcVL45VGOkMBhEjhsjr+euO2eqUO5cKDH/7DK+KbPH1t3ONr8c93h8lS4vz/in8Xz1Vdf0djYSEZGRtD2jIwMysrKWhyfl5fHE0880VXliYh0iVA/BwlgWk4mS24d1uJHxh2CH5ng5yy1VE88sa5UhlwyEjoRtNbv3M/3frOaXtTQy6ohhSOkWDUt3n8rx4UnocHf2xMIOycGnWr/F/oaoKbCvwApwNdaTj3T0jsn2ReXFBxcnP7wYuKSadxeyROxsdQQTw1OjhgnR0igxjipwclbr25jcvJoYpzJEJcM8UnHvy8m/pST/LUVgMq8tdy9bPNpB6Bwf39H2fqwQKvZyTDGtNgGMG/ePObMmRN439SDIiIiLU3LyWTyYHfIe2diHBbzZwzm7mWbsWi9h2b+jMGd/jsjz+5LL1eaPwA1ffkJf6TpEsx9N0w4eQDyNZ4QVg4fW6r4xxdlLFn1Eb2sWpKoDbwmU0ty0zZqGZIeQy+rDuqqA5/FHBsX03DEvxzeH/QnLeAaOPmv6lHghTb2WTHHwsoJoSXwPhlfXBK12w7xo9hYjuCk1sRTeywI1Zo4anGS/2oRkxMv9weg2AT/3DlxSRCX4H89SQhq9BmeeH1bq8HTHGvfE69vY/Jgd5eNp7EloPTt25eYmJgWvSXl5eUtelUAnE4nTmfbt72JiEiwcPTOQHh7aEIWgBwx/kG5iWcEbR50lmHD+n6nHqNzb7MAZAwcrfWHlWahx7+tmo+L9/LGpl0kWvUkUkcStSRadSThX5rWz0qFXlY9NByG+iP+Xh4A03js+1q/zOUAroNTB6AXT/YPYx0LLInHl1h/eKk6GsuPj9RQExdPjXFSSxyfmv4sa5zsLw8o9dayofhgWP5z1RpbAkp8fDzDhw8nPz8/aAxKfn4+1157rR0liYhIO4Wrh6bpuyMuAFnW8R/05L6tfvfh1AM8u6Ht8TlNfv+Ny4N/4Bsb/D0y9cd6ZuoPN3s9Ag2HKfq8lLc+/CwQgBKpI8GqJ5F6Eqg/tl7HWS4HqTEN0FADDbX+7zBND840/mDUcLhFXWcAk5td/lrTOCQQUJqUV7U+wDgcbLvEM2fOHG677TZGjBjB6NGj+fWvf80XX3zBXXfdZVdJIiLSTuHqoYHoDECnGp/T1DszMjsteEdMHMS4/HcenURV7wP84oN2BKDrLm95XppCUEPN8eVoTdD7T/fu57ertx0LP/U4rXr2mH4tvj89JeGUNYSKbQFl1qxZHDhwgB//+MeUlpaSk5PDG2+8wcCBA+0qSUREIkS0BaCwj8/pbACCdoWgsy8wFHxw6lvUW/3+MNFU9yIiIiESztt0w3kLeVd8P0TBPCinSwFFREQiVTgnOutJ86AooIiIiESRnjKTrK3zoIiIiEjHhHN8Tld8f3u1PWeviIiIiE0UUERERCTiKKCIiIhIxFFAERERkYijgCIiIiIRRwFFREREIo4CioiIiEQcBRQRERGJOAooIiIiEnGicibZptn5Kysrba5ERERE2qvpd7s9T9mJyoBSVVUFQFZWls2ViIiISEdVVVXhcrlOekxUPizQ5/Oxb98+UlJSsKzQPSAp0lRWVpKVlUVJSUmPeChiT2qv2tp99aT2qq3dV7jaa4yhqqoKj8eDw3HyUSZR2YPicDjo37+/3WV0mdTU1B7xX4gmPam9amv31ZPaq7Z2X+Fo76l6TppokKyIiIhEHAUUERERiTgKKBHM6XQyf/58nE6n3aV0iZ7UXrW1++pJ7VVbu69IaG9UDpIVERGR7k09KCIiIhJxFFBEREQk4iigiIiISMRRQBEREZGIo4DSxRYsWIBlWUGL2+0O7DfGsGDBAjweD4mJiYwbN46tW7cGfUddXR33338/ffv2JTk5mZkzZ7Jnz56ubkqr1qxZw4wZM/B4PFiWxauvvhq0P1Ttq6io4LbbbsPlcuFyubjttts4dOhQmFsX7FRtvf3221uc68svvzzomGhpa15eHpdddhkpKSmkp6dz3XXXsX379qBjusu5bU9bu8u5XbJkCRdffHFgMq7Ro0fzl7/8JbC/u5zTJqdqb3c5r63Jy8vDsixyc3MD2yL+/BrpUvPnzzcXXXSRKS0tDSzl5eWB/T/5yU9MSkqKefnll82WLVvMrFmzTGZmpqmsrAwcc9ddd5kzzzzT5Ofnm82bN5vx48eboUOHmqNHj9rRpCBvvPGGeeyxx8zLL79sALNy5cqg/aFq37Rp00xOTo5Zt26dWbduncnJyTHTp0/vqmYaY07d1tmzZ5tp06YFnesDBw4EHRMtbZ06dap5/vnnTVFRkSksLDTXXHONGTBggKmurg4c013ObXva2l3O7WuvvWb+/Oc/m+3bt5vt27ebH/7whyYuLs4UFRUZY7rPOW1ve7vLeW1uw4YN5qyzzjIXX3yxeeCBBwLbI/38KqB0sfnz55uhQ4e2us/n8xm3221+8pOfBLbV1tYal8tlnn32WWOMMYcOHTJxcXFmxYoVgWP27t1rHA6HefPNN8Nae0c1/9EOVfu2bdtmAPPee+8Fjlm/fr0BzD/+8Y8wt6p1bQWUa6+9ts3PRGtbjTGmvLzcAKagoMAY073PbfO2GtO9z23v3r3Nb37zm259Tk/U1F5juud5raqqMoMGDTL5+flm7NixgYASDedXl3hs8Omnn+LxeMjOzubGG2/ks88+A6C4uJiysjKmTJkSONbpdDJ27FjWrVsHwKZNm2hoaAg6xuPxkJOTEzgmUoWqfevXr8flcjFq1KjAMZdffjkulyvi/g1Wr15Neno65513Ht/97ncpLy8P7Ivmtnq9XgDS0tKA7n1um7e1SXc7t42NjaxYsYLDhw8zevTobn1OoWV7m3S383rvvfdyzTXXMGnSpKDt0XB+o/JhgdFs1KhR/O53v+O8887jyy+/5N/+7d8YM2YMW7dupaysDICMjIygz2RkZLB7924AysrKiI+Pp3fv3i2Oafp8pApV+8rKykhPT2/x/enp6RH1b3D11VfzrW99i4EDB1JcXMzjjz/OhAkT2LRpE06nM2rbaoxhzpw5XHnlleTk5ADd99y21lboXud2y5YtjB49mtraWnr16sXKlSsZPHhw4Melu53TttoL3eu8AqxYsYLNmzezcePGFvui4b+zCihd7Oqrrw6sDxkyhNGjR3POOefwwgsvBAZjWZYV9BljTIttzbXnmEgRiva1dnyk/RvMmjUrsJ6Tk8OIESMYOHAgf/7zn7n++uvb/Fykt/W+++7j448/Zu3atS32dbdz21Zbu9O5Pf/88yksLOTQoUO8/PLLzJ49m4KCgjZrjPZz2lZ7Bw8e3K3Oa0lJCQ888ACrVq0iISGhzeMi+fzqEo/NkpOTGTJkCJ9++mngbp7mqbO8vDyQct1uN/X19VRUVLR5TKQKVfvcbjdffvlli+/fv39/RP8bZGZmMnDgQD799FMgOtt6//3389prr/HOO+/Qv3//wPbueG7bamtrovncxsfHc+655zJixAjy8vIYOnQoP/vZz7rlOYW229uaaD6vmzZtory8nOHDhxMbG0tsbCwFBQX8/Oc/JzY2NlBLJJ9fBRSb1dXV8cknn5CZmUl2djZut5v8/PzA/vr6egoKChgzZgwAw4cPJy4uLuiY0tJSioqKAsdEqlC1b/To0Xi9XjZs2BA45v3338fr9Ub0v8GBAwcoKSkhMzMTiK62GmO47777eOWVV3j77bfJzs4O2t+dzu2p2tqaaD63zRljqKur61bn9GSa2tuaaD6vEydOZMuWLRQWFgaWESNGcMstt1BYWMjZZ58d+ef3tIbYSofNnTvXrF692nz22WfmvffeM9OnTzcpKSnm888/N8b4b/tyuVzmlVdeMVu2bDE33XRTq7d99e/f3/z1r381mzdvNhMmTIiY24yrqqrMhx9+aD788EMDmEWLFpkPP/zQ7N692xgTuvZNmzbNXHzxxWb9+vVm/fr1ZsiQIV1+G9/J2lpVVWXmzp1r1q1bZ4qLi80777xjRo8ebc4888yobOvdd99tXC6XWb16ddAtmEeOHAkc013O7ana2p3O7bx588yaNWtMcXGx+fjjj80Pf/hD43A4zKpVq4wx3eectqe93em8tuXEu3iMifzzq4DSxZruM4+LizMej8dcf/31ZuvWrYH9Pp/PzJ8/37jdbuN0Os1VV11ltmzZEvQdNTU15r777jNpaWkmMTHRTJ8+3XzxxRdd3ZRWvfPOOwZoscyePdsYE7r2HThwwNxyyy0mJSXFpKSkmFtuucVUVFR0USv9TtbWI0eOmClTpph+/fqZuLg4M2DAADN79uwW7YiWtrbWTsA8//zzgWO6y7k9VVu707n9zne+YwYOHGji4+NNv379zMSJEwPhxJjuc06bnKy93em8tqV5QIn082sZY8zp9cGIiIiIhJbGoIiIiEjEUUARERGRiKOAIiIiIhFHAUVEREQijgKKiIiIRBwFFBEREYk4CigiIiIScRRQREREJOIooIiIiEjEUUARERGRiKOAIiIiIhFHAUVEREQizv8HypHdZLmg2CIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Example of how we can use the data from our data table\n",
    "TT = data['T(K)']\n",
    "fkp_co2 = data['CO+1/2O eq CO2']\n",
    "\n",
    "lnkp_co2 = interp1d(TT,fkp_co2) #setup a function that will interpolate within the bounds of our table above\n",
    "# we can see how well the fxn interpolates\n",
    "Tnew = arange(200,4000)\n",
    "kpnew = lnkp_co2(Tnew)\n",
    "plt.plot(TT,fkp_co2,'o', Tnew,kpnew,'-')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now let's setup some Kp fxns:\n",
    "# import our required Kp eqns...\n",
    "fkp_h2o = data['H2+1/2O eq H2O']\n",
    "fkp_oh = data['1/2H2+1/2O2 eq OH']\n",
    "\n",
    "lnkp_h2o = interp1d(TT,fkp_h2o)\n",
    "lnkp_oh = interp1d(TT,fkp_oh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lnkp_oh(T) -0.4499\n",
      "lnkp_h2o(T) 8.1521\n",
      "lnkp_co2(T) 6.6301\n",
      "x_CO2 = 0.3249\n",
      "x_H2O = 0.6591\n",
      "ppmv_OH = 2238 ppmv\n"
     ]
    }
   ],
   "source": [
    "# Finally, let's return to our CH4+O2 example above:\n",
    "# Set up the knowns\n",
    "n_ch4 = 1\n",
    "n_o2 = 2\n",
    "\n",
    "T = 2000 #K\n",
    "P = 1 #atm\n",
    "\n",
    "#check we're pulling the correct Kp's\n",
    "print(\"lnkp_oh(T) %1.4f\" % lnkp_oh(T))\n",
    "print(\"lnkp_h2o(T) %1.4f\" % lnkp_h2o(T))\n",
    "print(\"lnkp_co2(T) %1.4f\" % lnkp_co2(T))\n",
    "\n",
    "# Define our equations to solve\n",
    "def eqn(vars):\n",
    "    co2, h2o, co, oh, o2, h2, nt  = vars\n",
    "    eq1 = n_ch4 - (co2 + co) #carbon balance\n",
    "    eq2 = 2*n_o2 - (2*co2 + h2o + oh + 2*o2 + co) #oxygen balance\n",
    "    eq3 = 4*n_ch4 - (2*h2o + oh + 2*h2) #hydrogen balance\n",
    "    eq4 = exp(lnkp_oh(T))  - (oh  / (o2**0.5 * h2**0.5))*1 #Kp_OH\n",
    "    eq5 = exp(lnkp_h2o(T)) - (h2o / (h2 * o2**0.5))*(1/nt)**(-.5) #Kp_H2O\n",
    "    eq6 = exp(lnkp_co2(T)) - (co2 / (co * o2**0.5))*(1/nt)**(-.5) #Kp_CO2\n",
    "    eq7 = nt - (co2 + h2o + co + oh + o2 + h2) #total moles\n",
    "    return [eq1, eq2, eq3, eq4, eq5, eq6, eq7]\n",
    "\n",
    "# Now call the fsolve function to solve our system of equatins\n",
    "# also note that you'll need a pretty good guess here...\n",
    "co2, h2o, co, oh, o2, h2, nt = fsolve(eqn, (.3, .7, .001, .001, .001, .001, 1))\n",
    "\n",
    "# Calculating mole fractions\n",
    "x_co2 = co2/nt\n",
    "x_h2o = h2o/nt\n",
    "ppmv_oh = oh/nt * 1e6\n",
    "\n",
    "print(\"x_CO2 = %1.4f\" % x_co2)\n",
    "print(\"x_H2O = %1.4f\" % x_h2o)\n",
    "print(\"ppmv_OH = %1.f ppmv\" % ppmv_oh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.9"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "position": {
    "height": "502.667px",
    "left": "1070px",
    "right": "20px",
    "top": "120px",
    "width": "350px"
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}