HW10-Q4-question.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homework 10, question 4\n",
    "First we need to set up the notebook to draw the plots inline, and import some libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy\n",
    "import scipy.integrate\n",
    "import scipy.optimize"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell will just try to make the graphs look nicer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    import seaborn as sns\n",
    "except ImportError:\n",
    "    # This block will be run if there's an ImportError, i.e you don't have seaborn installed.\n",
    "    sns = False\n",
    "    print (\"If you want to try different figure formatting, \"\n",
    "           \"type 'conda install seaborn' at an anaconda command prompt or terminal. \"\n",
    "           \"See https://stanford.edu/~mwaskom/software/seaborn/ for details\")\n",
    "    # If not using seaborn, we can still control the size of the figures this way\n",
    "    from pylab import rcParams\n",
    "    rcParams['figure.figsize'] = 4, 4\n",
    "else:\n",
    "    # This block will be run if there is no ImportError\n",
    "    sns.set_style(\"ticks\")\n",
    "    sns.set_context(\"poster\",rc={\"figure.figsize\": (4, 4)})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The system description\n",
    "The following complex gas phase reactions follow elementary rate laws\n",
    "$$\\require{mhchem}$$\n",
    "(1) $$\\ce{A + 2B -> C}$$\n",
    "(2) $$\\ce{2A + 3C -> D}$$\n",
    "and take place in a PFR. The feed is stoichiometric for reaction (1) in A and B with $F_{A0}$ = 5 mol/min. The reactor volume is 10 dm<sup>3</sup> and the total entering concentration is $C_{T0}$ = 0.2 mol/dm<sup>3</sup>. The entering pressure is 100 atm and the entering temperature is 300 K. The coolant flow rate is 50 mol/min and the entering coolant fluid has a heat capacity of $C_{P,cool}$ = 10 cal/mol/K and enters at a temperature of 325 K.\n",
    "\n",
    " * $-r_{1A} = k_{1A} C_A C_B^2$\n",
    " * $k_{1A} = 40$ dm<sup>6</sup>/mol<sup>2</sup>/min at 300 K with $E_1 = 8000$ cal/mol\n",
    " * $\\Delta H_{Rx1B}= -1500$ cal/mole B \n",
    " * $-r_{2C} = k_{2C} C_A^2 C_C^3$\n",
    " * $\\Delta H_{Rx2A} = -10,000$ cal/mole A\n",
    " * $k_{2C} = 2$ dm<sup>12</sup>/mol<sup>4</sup>/min at 300K with $E_2 = 12,000$ cal/mol\n",
    " * $C_{P,A} = 10$ cal/mol/K\n",
    " * $C_{P,B} = 12$ cal/mol/K\n",
    " * $C_{P,C} = 14$ cal/mol/K\n",
    " * $C_{P,D} = 16$ cal/mol/K\n",
    " * $Ua = 80$ cal/min/K\n",
    " * $T_{a0} = 325$ K\n",
    " * $\\dot{m} = 50$ mol/min\n",
    " * $C_{P,cool} = 10$ cal/mol/K\n",
    " \n",
    "Plot $F_A$, $F_B$, $F_C$, $F_D$, $y$, $T$, $T_a$ as a function of $V$ for:\n",
    "\n",
    " * (a) Co-current heat exchange (cooling fluid at $\\dot{m}$ enters at $T_{a0}$ at $V=0$)\n",
    " * (b) Counter-current heat exchange (cooling fluid at $\\dot{m}$ enters at $T_{a0}$ at $V=10$ dm<sup>3</sup>)\n",
    " * (c) Constant temperature heat exchange (as if $\\dot{m}=\\infty$)\n",
    " * (d) Adiabatic (no heat exchange)\n",
    " \n",
    " **Notes:**\n",
    "Neglect pressure drop.\n",
    "The description above uses Fogler's syntax where $k_{2C}$ reminds you that it is for a per-mole-of-C basis rate law, and $\\Delta H_{Rx2A}$ means the enthalpy change of reaction 1 on a per-mole-of-A basis (though it is per mole of A consumed, not produced; the forward reaction is exothermic). The tests below assume you're doing everything on a \"per mole of reaction as written\" basis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "DHr1 = -15000. * 2 # value given is per mole of B\n",
    "DHr2 = -10000. * 2 # value given is per mole of A\n",
    "Cpa = 10 # cal/mol/K\n",
    "Cpb = 12 # cal/mol/K\n",
    "Cpc = 14 # cal/mol/K\n",
    "Cpd = 16 # cal/mol/K\n",
    "Ua = 80 # cal/min/K\n",
    "Ta0 = 325 # K\n",
    "mdot = 50 # mol/min\n",
    "Cpcool = 10 # cal/mol/K\n",
    "Fa0 = 5 # mol/min\n",
    "Fb0 = 10 # mol/min\n",
    "Fc0 = 0\n",
    "Fd0 = 0\n",
    "P0 = 100 # atm\n",
    "T0 = 300 # K\n",
    "Ct0 = 0.2 # initial total concentration, mol/dm3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def k1(T):\n",
    "    \"Return the rate constant for reaction 1 at the given temperature\"\n",
    "\n",
    "    return 0####FIXME####"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Some tests\n",
    "assert k1(300) == 40.\n",
    "assert round(k1(400),1) == 1146.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def k2(T):\n",
    "    \"Return the rate constant for reaction 2 at the given temperature\"\n",
    "\n",
    "    return 0####FIXME####"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Some tests\n",
    "assert k2(300) != 2., \"Beware: The given value is per mole of C\"\n",
    "assert round(k2(300),2) == 0.67\n",
    "assert round(k2(400),2) == 102.23"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def dVector_dV(vector, V, heatexchange='co-current'):\n",
    "    \"\"\"\n",
    "    Returns the volume-derivative of a vector containing\n",
    "    the variables (Fa, Fb, Fc, Fd, y, T, Ta).\n",
    "    The parameter `heatexchange` must be one of \n",
    "    'co-current', 'counter-current', 'constant-Ta', or 'adiabatic'.\n",
    "    If not specified, the default is 'co-current'.\n",
    "    \"\"\"\n",
    "    Fa, Fb, Fc, Fd, y, T, Ta = vector\n",
    "    assert heatexchange in ['co-current', 'counter-current', 'constant-Ta', 'adiabatic']\n",
    "\n",
    "    Ft = 0####FIXME####\n",
    "    Ct = 0####FIXME####\n",
    "    Ca = 0####FIXME####\n",
    "    Cb = 0####FIXME####\n",
    "    Cc = 0####FIXME####\n",
    "    Cd = 0####FIXME####\n",
    "    r1 = 0####FIXME####\n",
    "    r2 = 0####FIXME####\n",
    "\n",
    "    ra = 0####FIXME####\n",
    "    rb = 0####FIXME####\n",
    "    rc = 0####FIXME####\n",
    "    rd = 0####FIXME####\n",
    "    \n",
    "    # Qg is \"heat generated\"\n",
    "    Qg = 0####FIXME####\n",
    "    \n",
    "    # Qr is \"heat removed\"\n",
    "    if heatexchange == 'adiabatic':\n",
    "        Qr = 0####FIXME####\n",
    "    else:\n",
    "        Qr = 0####FIXME####\n",
    "    \n",
    "    sumFiCpi = 0####FIXME####\n",
    "    \n",
    "    dFa_dV = 0####FIXME####\n",
    "    dFb_dV = 0####FIXME####\n",
    "    dFc_dV = 0####FIXME####\n",
    "    dFd_dV = 0####FIXME####\n",
    "    dy_dV = 0\n",
    "    dT_dV = 0####FIXME####\n",
    "\n",
    "    if heatexchange == 'co-current':\n",
    "        dTa_dV = 0####FIXME####\n",
    "    elif heatexchange == 'counter-current':\n",
    "        dTa_dV = 0####FIXME####\n",
    "    else:\n",
    "        dTa_dV = 0\n",
    "    return (dFa_dV,\n",
    "            dFb_dV,\n",
    "            dFc_dV,\n",
    "            dFd_dV,\n",
    "            dy_dV,\n",
    "            dT_dV,\n",
    "            dTa_dV)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Some tests\n",
    "np.testing.assert_allclose(dVector_dV((Fa0, Fb0, Fc0, Fd0, 1.0, T0, Ta0), 0), (-0.047407407407407398, -0.094814814814814796, 0.047407407407407398, 0.0, 0, 20.130718954248362, -4))\n",
    "np.testing.assert_allclose(dVector_dV((10, 10, 10, 10, 1.0, 300, 600), 0, heatexchange='counter-current'), (-0.0050004166666666678, -0.010000000000000002, 0.0049993750000000012, 2.0833333333333339e-07, 0, 46.442315705128202, 48))\n",
    "np.testing.assert_allclose(dVector_dV((10, 10, 10, 10, 1.0, 300, 600), 0, heatexchange='constant-Ta'),(-0.0050004166666666678, -0.010000000000000002, 0.0049993750000000012, 2.0833333333333339e-07, 0, 46.442315705128202, 0))\n",
    "np.testing.assert_allclose(dVector_dV((10, 10, 10, 10, 1.0, 400, 600), 0, heatexchange='adiabatic'),(-0.060448050006020868, -0.12086577476161001, 0.060410143442981207, 7.5813126079323574e-06, 0, 3.4868043224544398, 0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def make_plots(odeint_result):\n",
    "    \"\"\"\n",
    "    Make some plots from the given odeint result,\n",
    "    and print the selectivity of C to D\n",
    "    \"\"\"\n",
    "    Fa, Fb, Fc, Fd, y, T, Ta = odeint_result.T\n",
    "    plt.plot(volumes, T, label=\"T\")\n",
    "    plt.plot(volumes, Ta, label=\"Ta\")\n",
    "    plt.legend(loc=\"best\")\n",
    "    plt.show()\n",
    "\n",
    "    plt.plot(volumes, Fa, label=\"F$_A$\")\n",
    "    plt.plot(volumes, Fb, label=\"F$_B$\")\n",
    "    plt.plot(volumes, Fc, label=\"F$_C$\")\n",
    "    plt.plot(volumes, Fd, label=\"F$_D$\")\n",
    "    plt.legend(loc=\"best\")\n",
    "    plt.show()\n",
    "\n",
    "    ## Skip this, it's boring because we assume no pressure drop\n",
    "    #plt.plot(volumes, y, label=\"y\")\n",
    "    #plt.legend(loc=\"best\")\n",
    "    #plt.show()\n",
    "    \n",
    "    print(\"Maximum temperature reached is {0:.0f} K\".format(T.max()))\n",
    "    \n",
    "    selectivity = Fc[-1] / Fd[-1]\n",
    "    print(\"Exit flow rates Fc={0:.3g}, Fd={1:.3g} mol/min\".format(Fc[-1] , Fd[-1]))\n",
    "    print(\"Selectivity C/D = {0}\".format(selectivity))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Co-current heat exchange"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASIAAAEVCAYAAABEweijAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XtclGXC//HPzMBwmOGMgAgCigc8i6JC2spqim3bWj91\nK23zVPrss9pmudmaqdVDrnlYixTN6Mm2fbYyDdt2axMzPGGmqXlARRwFBZEzMxyGGe7fHyMDEyCg\nwMzA9X69eDXcc8/NNcl8uU73dckkSZIQBEGwIrm1CyAIgiCCSBAEqxNBJAiC1YkgEgTB6kQQCYJg\ndSKIBEGwuhYFUXl5Oa+++ir33XcfkZGRzJs3j/T0dItztmzZQmxsLMOGDWPu3LlkZmZaPK/X64mP\nj2fs2LFERkayePFi8vLy2u6dCIJgt2QtmUc0b948fvzxRxYvXkzfvn3Zs2cP33zzDZ999hmhoaEk\nJCSwfft2li5dSmBgIJs3byYvL48vv/wStVoNwEsvvcS3337LsmXLcHV1Zf369bi6urJr1y5kMtk9\nv5FHH30UgF27dt3ztTqKvZXZ3soLoswdoS3K69DcCWfPnuXQoUO89tprTJ8+HYCYmBg0Gg2bNm3i\n9ddfJykpiUWLFjFz5kwARowYQWxsLDt37mT27Nlcu3aN5ORkNmzYQFxcHAD9+vUjLi6OlJQUJk6c\neNdvQBAE+9ds00yj0SCTybjvvvssjkdGRnLw4EFOnTpFRUUFsbGx5ufc3d2JioriwIEDAKSlpSGT\nyRg/frz5nJCQEMLDw0lNTW2jtyIIgr1qNogCAgKQJIkbN25YHM/Ozkar1XL69GkAevbsafF8cHAw\nGo0GMIWZr68vzs7OTZ4jCELX1WwQDR48mNDQUFavXs2ZM2coKyvjH//4h7kmU1NTg1KpxMHBspWn\nUqnQarUAaLVaVCpVg2vXP0cQhK6r2T4ipVLJO++8w/PPP2/uIxo2bBjz58/nnXfeQS6XN9nZLJfX\n5VxLzmmp2s6x+i5dukRgYCBXr15FoVC0+prWoNfrAVPt0h7YW3lBlLkjVFVVcfXqVaZOndrg89zS\nDuxmgwigd+/efP7559y8eRODwUCPHj1ISEhALpfj6uqKXq/HaDRaBIBOp8PNzQ0AtVqNTqdrcN36\n59yrwMBANBoNkyZNapPrdaQJEyZYuwitYm/lBVHmjlBRUdFoy6clmg2iyspKvv76a6Kjo/H39zcf\nv3jxIn369KF3797U1NSQnZ1NSEiI+fmsrCzCwsIACA0NJT8/H71ej1KptDgnKiqq1YVuLGWvXr3K\npEmT+OijjwgICGj1NQVBuDu5ubnMnDmTbdu2WWRAazQbRA4ODqxatYpnn32W2bNnA6YA+e6775g7\ndy7Dhw/HycmJvXv3Mm/ePABKSko4duwYixYtAiA6OhqDwcC+ffvMw/cajYaMjAyeffbZuyr4z9XW\nxgICAggKCmqTawqC0HL30iXSoiCaNm0aiYmJeHt7o1KpWL9+PT4+PsyePRtXV1dmzZrFpk2bkMlk\nhISEkJiYiLu7O9OmTQNMo2NxcXGsWLGCsrIy3Nzc2LhxIxEREXZX/RQEoe21qI/ohRdeQC6X8+ab\nb6LX6xkzZgxLly7Fw8MDgCVLlqBQKEhKSqK8vJzIyEjWrl1rnlUNsGbNGuLj41m3bh2SJBETE8Py\n5cvbZFa1IPxcZZWBv32VTkWVgam/6E2wf9v0RQrto0W3eNiD7OxsJkyYQEpKimiadXHVBiOvvXeU\nHy/eAkAhl/Gr+8J48sEInJUt+tsrtEJbfPbE3fdCp2I01rDuo+PmEAIw1kjsOZDJ25+ctGLJhDsR\nQSR0KsmpmRw+nQOAo4OcyWNCcHQw/Zqn/nids5kF1iye0AQRREKn8u3xLPPjF58cyR+mD2Purwea\nj21P/omamk7RG9GpiCASOo2cfB2anFIAwoM8GD2oOwBx0aEE+5sGTjKyS9j3wzWrlVFonAgiodNI\nO5NjfjzmdggBOCjkzH94sPn7v32VjtFY06FlE+5MBJHQaTQVRACR/f0Y1rcbAAUllRw7f7NDyybc\nmQgioVMoLqvivKYQgO4+KnoGNJw39GBMmPnx12lXO6xsQvNEEAmdwvfncqmdETd6UECjE2WjBvjj\n5eYEwIn0m+QXV3RkEYU7EEEkdAp3apbVclDImTjKtIBfjQTffC86rW2FmGZqx1566SV2797d5PMy\nmYz//u//5g9/+EMHlqrjSZLEuSumZpmbqyP9Q72bPPeBUSF8mnIJgG++v8qMiX1RyMVtRtYmgsiO\n/f73v+fxxx83f7906VLCwsL4/e9/bz5Wf+mWziq/uBJdRTUAvXp43DFYuvuqGNrHl1OX8rlVVMGZ\ny/kM7dOto4oqNEEEkR0LDg4mODjY/L2LiwteXl4MGTLEiqXqeFdzS82PQ7q7N3v+L0cGc+pSPgCH\nTt8QQWQDRB+RYPdqJzEChLUgiEYN7I6DwlRrOvJTDkYx09rqRBAJdk9zo3U1IrWLI8P6+gGmYf9z\nV8T9Z9bW5ZtmB368zkdfn6eiymC1Mrg4OTBzcgTjhvewWhnsWW3TTCajxesO3TckkB9uT2o8dOoG\ng3v7tlv5hOZ1+SDatf8S1281XNi/Y1Wx67sMEUR3odpQQ3ZeGQCBvqoWrzc0ZlAACZ/KMNZIHD59\ng2emDkYuRs+spssH0aPj+9hEjejR8eFW+/n27PotLQajqY+nJc2yWmpXJUP7duNEeh5Ft2dlD+zl\n017FFJrR5YNo3PAeoiZix+p3VIcGtDyIAMYOCeREeh5g6rQWQWQ9orNasGtX6wdRYOuCaNTAAGrv\nBPn+bC6dZNVkuySCSLBr9WtErWmaAXionegfYpqFnVOgIztPbH9uLSKIOhGZTNbldkWpDSInpYIA\n79bvMjpqYN1mnN+fzW2zcgmtI4KoE9m9ezfx8fHWLkaH0VZUm++gDwlwu6tRr9H1guioCCKrEUEk\n2K3c/LppF0F+d7dvWZCfmu4+pppU+tVCSrRVbVI2oXVEEAl262ZRufmxv7frXV1DJpOZm2eShHmS\no9CxWhRENTU1vPvuu0yaNInhw4czY8YM0tLSLM7ZsmULsbGxDBs2jLlz55KZmWnxvF6vJz4+nrFj\nxxIZGcnixYvJy8tru3cidDl5hfceRACjBtatUCCaZ9bRoiDavn07f/3rX5k2bRqbN28mODiY+fPn\nk56eDkBCQgJbt25l/vz5bNy4kbKyMubMmYNWWzcKsXLlSvbs2cMLL7zAmjVruHDhAgsWLBBDpsJd\nqx9EfvcQRAPCfFA5m6bUnbx4i2qDWFi/o7UoiD7//HMefvhhnnnmGaKjo3nzzTfp1q0bO3fuRKfT\nkZSUxKJFi5g5cyaxsbG89957aLVadu7cCcC1a9dITk5m1apVTJ06lUmTJrFt2zbS09NJSUlp1zco\ndF4WTTOvuw8iB4Wc4f1MN8FWVBk4JzZh7HAtCiK9Xo9KVTc0KpfLUavVFBcXc+rUKSoqKoiNjTU/\n7+7uTlRUFAcOHAAgLS0NmUzG+PHjzeeEhIQQHh5OampqG70VoauprRHJ5TJ8PJzv6VpRA+qaZ2KH\nj47XoiCaOXMmycnJHDlyBK1WywcffMDly5d56KGHuHLlCgA9e/a0eE1wcDAajQYAjUaDr68vzs7O\nTZ4jCK0hSRI3bweRr6cLCsW9jbuM6O9vnmX9w3nRT9TRWnSv2eOPP05aWhpz5swBTCMNf/zjHxk/\nfjzbtm1DqVTi4GB5KZVKZe4j0mq1FjWq+ufk5op/dKH1SnV6KvVG4N6aZbU81E70DfbiwrUirt/S\ncSNfS6Cv+p6vK7RMi4KodhRs9erV9OrVi8OHD/P222+jVquRJKnJ2bxyed1fqZacIwgtlVdUv6Pa\npU2uOXKAPxeuFQHww7mbPHy/CKKO0mwQHT9+nBMnTvDWW28xadIkAKKiojAYDKxbt47nnnsOvV6P\n0WhEoVCYX6fT6XBzM00yU6vV6HQN1/ypf05rPProow2O6fX6Vl9HsF95hXV7krVFjQhgZIQ/H31l\nGgk+dv4mD9/fu02u21UsXLgQpVJpcWzXrl0tem2z1ZHc3FxkMhlDhw61OD5ixAgqKyuRy+VIkkR2\ndrbF81lZWYSFmXbWDA0NJT8/v0FY1D9HaL2XXnqJ/v37N/kVERFBQkKCtYvZLm620dB9fb17eODt\nbtqA8czlAquuUdXVNFsjCg0NRZIkjh8/zoMPPmg+fvLkSRQKBQ888ABvvvkme/fuZd68eQCUlJRw\n7NgxFi1aBEB0dDQGg4F9+/YRFxcHmDqwMzIyePbZZ1td6MZSNjs7mwkTJrT6WvasK28nlNcGs6p/\nTiaTMaK/P998fw2DsYaTF28RPbjxzRqFhhITEwkKCrqr1zYbRAMHDmT8+PGsXr2a4uJievfuzdGj\nR9m+fTtPPfUU/v7+zJo1i02bNiGTyQgJCSExMRF3d3emTZsGmEbH4uLiWLFiBWVlZbi5ubFx40Yi\nIiK6XHi0pa68nVB71IjANIxfuwPsD+dviiDqIC3qrH7rrbfYuHEjW7dupaSkhJCQEF555RVmzJgB\nwJIlS1AoFCQlJVFeXk5kZCRr165Fra7r7FuzZg3x8fGsW7cOSZKIiYlh+fLlXW7ZCmu4efMmGzdu\n5ODBgxQXF+Pt7c2DDz7ICy+80GC0017U1ogUchk+7vc2h6i+oX264aCQYTBK/HD+5h0HY4S206Lf\nQqVSyYsvvsiLL77Y6PMKhYIlS5awZMmSJq/h7OzMq6++yquvvnp3JRXuSk1NDfPmzUOpVPLqq6+i\nUqlITU3lvffeIzQ0lMcee8zaRWw1SZLMkxnbYg5Rfa7Ojgzs5cOpS/kUllZy5UYpvXp4tNn1hcbZ\n55/DNnT42g98fOYLKqutt/yDs6MTvx30a2J6jmzza+fm5uLl5cWqVavo3ds0CjR69Gi+++47jh07\nZpdBZDGHqA2bZbVGRgSYd4I9dj5XBFEH6PJBtCf9G3LKrLwKQCV8kb63XYIoMDCQDz/8EEmS0Gg0\naDQa0tPTKSwstNspDxb9Q200dF9f1AB/3ttzBjDNJ/rtxH5t/jMES10+iB7u/4BN1Ih+3f+Bdrv+\nxx9/zFtvvUVhYSF+fn4MHToUZ2dnu135wHIyY9sHUY9uarr7qsjJ13HhWhEl2io81E5t/nOEOl0+\niGJ6jmyXmoitOHLkCCtXruTZZ5/l8ccfx9PTE2h8Uqi9qF0eFqCbZ9vMqv65qAh/9hzIRJLgxIU8\nYkcEN/8i4a6J+ys6uVOnTuHg4MDChQvNIZSbm8ulS5fstkZUUFJpfnyvd903pf7d+GJR/fbX5WtE\nnd3gwYMxGo38z//8D5MmTeL69eskJiZiNBqprKxs/gI2qCOCaGAvX1ycHKioMnDiQh4GYw0ObTg6\nJ1gSQdSJNLad0H333cef/vQnPvzwQ3bu3ElAQAAPPfQQkiTx0UcfNbhH0B4UlNQ1zXw82qdp5ugg\nJ7KfH4dO36C80sDZzAKG9unWLj9LEEHUqezevbvR43PmzDEv4VLf4sWL27tI7aK2RuSsVODq3H6/\nwlED/Dl0+gYAx87dFEHUjkRdU7ArkiRRWGoKIh8P53ad9Twyom6xtGPnRD9RexJBJNiVUp3evLh9\nezXLanmonejX0wuAG/k6svPK2vXndWUiiAS7UlsbAvBup47q+iy3pBZrWbcXEUSCXbEYMWvDm12b\nMmpA/S2pc9r953VVIogEu9IRI2b19QxwI8DHNHs7XVNIcZnYkro9iCAS7EpHzCGqTyaTMWaQaU2i\nGkl0WrcXEUSCXenoIAIYPbB+80wEUXsQQSTYlY5umgFEhHrj5mpaFP7Hi7eo1Iu1rNuaCCLBrtTW\niOQy8HLrmDviFQo5owaa7j3TVxv58cKtDvm5XYkIIsGu1AaRp5tTm67M2JzRA+vWrk47I0bP2poI\nIsFu6KuNlJWbFnPz7qBmWa3h/bqhdDTdk/f92VwMxpoO/fmdnQgiwW7Un8zYEXOI6nNWOjCivx8A\n2opqfsrI79Cf39mJIBLshjVGzOqLGRJofnz4J9E8a0siiAS7YY0Rs/pGDfA3r0mU9lMOxhr7XFjO\nFokgEuyGtWtErs6ODO9nWgqkWFvFuSsFHV6GzkoEkWA3rB1EADGD6zXPbq9VJNw7EUSC3bDorLZC\n0wxg9KAAFHLTIkWHT+dQI5pnbaLZIPr+++/p379/k185OaZOuy1bthAbG8uwYcOYO3cumZmZFtfR\n6/XEx8czduxYIiMjWbx4MXl5Vt5PTLArln1E1qkRubkqzSs1FpZWiuZZG2l2nc2BAwfyySefWByr\nrKxk8eLFDB48mO7du5OQkMD27dtZunQpgYGBbN68mTlz5vDll1+iVqsBWLlyJd9++y3Lli3D1dWV\n9evXs2DBAnbt2iX2FhdapLZG5KxU4OJkvVWOxw3rwYkLpj+iqSevM6i3r9XK0lk0WyNSqVQMGTLE\n4uubb75BLpfz5ptvotPpSEpKYtGiRcycOZPY2Fjee+89tFotO3fuBODatWskJyezatUqpk6dyqRJ\nk9i2bRvp6emkpKS0+5sU7J8kSeY+ovZeIrY5YwZ3N4+eHTp1Q0xubAOt7iPKyMjg73//O8899xye\nnp6cPHmSiooKYmNjzee4u7sTFRXFgQMHAEhLS0MmkzF+/HjzOSEhIYSHh5Oamnrv70Lo9LQV1eYl\nYr3drdM/VEvt4sjICNPkxlKdntOXxOTGe9XqINq4cSNhYWFMnz4dAI1GA0DPnj0tzgsODjY/p9Fo\n8PX1xdnZuclzBOFObGHErL77hwWZH3/3Y7YVS9I5tCqIsrKy+Pbbb5k7d675mE6nQ6lU4uBg2WZX\nqVRotVoAtFotKpWqwfXqnyMId1K/o9q7g2/vaEzUAH+clKZ7z9LO5KCvNlq5RPatVT1+n376KR4e\nHvz61782H5Mkqcn2ulxel3MtOaelGtu3Xa/Xt/o6gv0otLEakbOTA6MHBpD643XKKw0cO3eT+4YG\nNv/CTmzhwoUolUqLY7t27WrRa1uVAikpKUycOBFHR0fzMbVajV6vx2i0/Iug0+lwc3Mzn6PT6Rpc\nr/45gnAnBTYwh+jnxkfWNc/2/ZBlxZLYvxbXiHJycrh8+TLLli2zOB4aGookSWRnZxMSEmI+npWV\nRVhYmPmc/Px89Hq9RWJmZWURFRXV6kI3lrLZ2dlMmDCh1dcS7EP9GpEtNM0AIvv54al2olhbxfH0\nmxSXVeHZQYu12aLExESCgoKaP7ERLa4RnT59GplMxtChQy2ODx8+HKVSyd69e83HSkpKOHbsGNHR\n0QBER0djMBjYt2+f+RyNRkNGRgYxMTF3VXCha7G1zmowrdz4i9u1ImONRKrotL5rLa4RXbp0CS8v\nL9zd3S2Ou7q6MmvWLDZt2oRMJiMkJITExETc3d2ZNm0aYBodi4uLY8WKFZSVleHm5sbGjRuJiIgQ\ntRihRQpL6zqrvWykRgQwISqY5NTLAKT8kMXD9/e2consU4uDqKCgoEEI1VqyZAkKhYKkpCTKy8uJ\njIxk7dq15lnVAGvWrCE+Pp5169YhSRIxMTEsX75czKoWWqS2RuShVuLoYDu3SIYFehDa3R1NTimZ\n10vQ5JQS2r3xz4nQNJkkSZ3irr3aPqKUlJS7bqcKtslorOGRF79AkiAs0J23no9t/kUd6PPvMnhv\nz1kAfnN/b+b/ZpCVS9Sx2uKzZzt/WgShCUVlVdT+ubSVEbP6fhEZZL4jf98PWWJO0V0QQSTYPMvl\nP2ynf6iWl5uzeTfYsnI9R8Qysq0mgkiwebY2q7oxk8fUTV35z9GrViyJfRJBJNg8Wxy6/7mhfbrh\n7+0KwOmMfG7cErcutYYIIsHm1W+a2WqNSC6XMWl0Xa3o6zRRK2oNEUSCzbOsEdleZ3WtiaN6Ir/d\naf3N91epEp3WLSaCSLB5tnbDa1O83Z2JGVzbaV3NdyfETOuWEkEk2LyC27OqHRQy3FyVzZxtXb8e\n18v8+IsDmXSSaXrtTgSRYNMkSSK/2BRE3u7O5qaPrYoI9aZ3kAcAmpxSzmSKxfVbQgSRYNN0lQYq\nqkx9Ld28XK1cmubJZDIeus+yViQ0TwSRYNNqa0MA3Txtt6O6vvuH98BdZWpCHj2Tw418MZTfHBFE\ngk27VVRufuxrJ0GkdFTwYIxpLa4aCXbvv2zlEtk+EUSCTbOoEXnZRxABPDQ2DKWjaU3rlGPXKKo3\nF0poSASRYNNu1Qsie6kRAXionZg02rSzTbWhhi8Oir6iOxFBJNi0W3bYR1Rr6i/CzaN8/zp0hfLK\naiuXyHaJIBJs2q0i+w0if29X7h/WAzCN/olaUdNEEAk2rbZG5OKkQOXi2MzZtmf6hD7ULkK6e/9l\ndBWiVtQYEUSCzTLWSBTeXgLE19PFLpcV7hngbt4VVldRzZ5UMYLWGBFEgs0qLqvEYDTdItHN0/Yn\nMzblsUl9qZ0Qnpx6GW252Az050QQCTbLXofufy7Iz43xI4IBU1/Rrv0ZVi6R7RFBJNgsex26b8xv\nH+hrXtc6+bvLFp3wgggiwYbZ4+0dTQn0VTMlJhQAvaGGv3113roFsjEiiASbVb/WYO81IoDHHuiH\nq7NpK8Fvj2eRkV1s5RLZDhFEgs261Un6iGp5qJ2YPqEvAJIE25PPiPWKbmtxEB05coQZM2YwdOhQ\nfvnLX/L2229TU1Njfn7Lli3ExsYybNgw5s6dS2am5eQtvV5PfHw8Y8eOJTIyksWLF5OXl9d270To\ndCz6iGx4idjWeHhcL/xuh+rZzAL2i1UcgRYG0fHjx3n66acJDw9n27ZtzJo1i3fffZctW7YAkJCQ\nwNatW5k/fz4bN26krKyMOXPmoNXWLX+wcuVK9uzZwwsvvMCaNWu4cOECCxYsEH8RhCbl326aeaqd\nzDeQ2julo4L5vxls/j5pz1m0YpJjy4Jow4YNjBs3jvj4eEaPHs3cuXOZPXs2R48eRafTkZSUxKJF\ni5g5cyaxsbG89957aLVadu7cCcC1a9dITk5m1apVTJ06lUmTJrFt2zbS09NJSUlp1zco2Cd9tZFi\nbRUAvp62u0713RgzKICREf4AFGur+OjfouO62SAqLCzkxIkT/Pa3v7U4vmTJEnbs2MGpU6eoqKgg\nNrZuP3J3d3eioqI4cOAAAGlpachkMsaPH28+JyQkhPDwcFJTU9vorQidSX4nGrr/OZlMxoJHBqN0\nMH38vjx8hfNXCq1cKutqNoguXrwIgJOTEwsXLmTIkCHExMSQkJCAJElcuXIFgJ49e1q8Ljg4GI1G\nA4BGo8HX1xdnZ+cmzxGE+nIKdObHAT4qK5akfQT4qJjxQF3H9aaPT3Tp7YdaVCOSJIlly5bRu3dv\ntm/fzhNPPEFiYiLbt29Hp9OhVCpxcHCweJ1KpTL3EWm1WlSqhr9M9c8RhPpy8+uCqLtv5wsigP8X\n24dePUwL7V+/peNvXbiJ5tDcCQaDAYBx48axdOlSAEaNGkVRURFbtmzhmWeeafJmRLm8Ludack5L\nPfroow2O6fXi/p3O5Ea9GlH3TlgjAnBQyPnjY8N5buN3GGskklMvM3pgAIN6+1q7aHdl4cKFKJWW\n2z3t2rWrRa9tNgVcXU03G44dO9bieExMDBUVFbi5uaHX6zEaLauVOp0ONzc3ANRqNTqdjp+rf44g\n1JfTBWpEAGGBHvz2gX6AqYm27qPjlOq63h/VZmtEISGm/byrqy2HGGtrSkqlEkmSyM7ONp8LkJWV\nRViYaQHx0NBQ8vPz0ev1FomZlZVFVFRUqwvdWMpmZ2czYcKEVl9LsE25t2tECrnM7m/vaM6MCX34\n8UIe5zWFFJRUsukfP/Ly3FF2t+xJYmIiQUFBd/XaZmtE4eHh+Pv789VXX1kc379/P35+fjz44IMo\nlUr27t1rfq6kpIRjx44RHR0NQHR0NAaDgX379pnP0Wg0ZGRkEBMTc1cFFzqvmhqJ3ALT7h3+3q4o\nFJ37BgCFQs4LM0eYF377/lwuyV1s3aJma0QymYznnnuOl156iVWrVjF58mQOHz5McnIyq1evRqVS\nMWvWLDZt2oRMJiMkJITExETc3d2ZNm0aYBodi4uLY8WKFZSVleHm5sbGjRuJiIgQtRihgYKSSqoN\npln7AZ24WVafn7cri2cM440PjgHw/hdnCe3uzrC+flYuWcdoNogApk6dilKpJDExkd27dxMQEMDq\n1auZPn06YJpTpFAoSEpKory8nMjISNauXYtarTZfY82aNcTHx7Nu3TokSSImJobly5fbXfVTaH85\nBXUjqYGdtKO6MTFDAnlkfDi792dQI8HaD39gwx9/0SmnL/ycTOok91jU9hGlpKTcdTtVsA1fp10l\n4dOTADz9m0E8fH9vK5eo4xiNNazansbJi7cACPZXs/YP41C7Kpt5pfW0xWevcze+BbuUU2+L5q7S\nNKulUMj505MjCfAxjVZn3dTy+vvfo+/kkx1FEAk2J6cLzCG6EzdXJSvnj8Htdi3obGYB6z46jtFY\n08wr7ZcIIsHm5OabRsxkMsw1g64myM+NV+aPNq86cOSnHDb83wmMNZ2iJ6UBEUSCTZEkydxZ7evp\ngqND51j+4270D/Hmxd+NNK91nfrjdTb9o3OGkQgiwaaUaPVUVJn6Q7pis+znRg0IYOmTI81bV397\nPJu1Hx6j2tC5+oxEEAk2pavc2tEa9w0J5IUnRpjD6PDpHFZvT6O8svMsqCaCSLAp9ecQiRpRnXHD\ne7B89ijzGkanLuXzp7cPmG+FsXciiASbci23zPy4h5/6Dmd2PaMGBrD6mWhUt3cCuZpbxvObUvkp\nI9/KJbt3IogEm6LJKTU/Du3ubsWS2KZBvX15c/H95tHEUp2elxMP8fHeC9TYcSe2CCLBply9HUTO\nSgV+Xl1z6L45wf5urH/2FwwJN61bVCPB3/6dzivbDtvtDrIiiASboa2oJr+kEoCQAHdz56zQkLtK\nyavPRDN9Qh/zsVOX8vnDun3s/f6q3e2OI4JIsBlX6zXLegaIBfOao1DI+d2DA1j19Bg83ZwAKK80\nsOnjk/wvP+RhAAAgAElEQVR5yyGu5ZY2cwXbIYJIsBlXc0X/0N0Y0d+fd5b+knHDepiPnblcwOL1\n+9m6+zQlt7dlsmUiiASbUb9GFCKCqFXcVUr+9ORI/jx7lHl7bmONxD8PXuGZN/byf1+n2/RGjiKI\nBJtxtd7QfUiACKK7ET24O5uX/pLpE/qY5xyVVxr4+38uMP/1//Dhv89TWFpp5VI2JIJIsAmSJJmH\n7j3VTuY+D6H1nJ0c+N2DA9iybAK/HBlMbZ+/rtLAJ3svMu/1/7Dh78c5f6XQZjq1W7RCoyC0t8LS\nSnS3mw4h3UVHdVvw83LluccjmT6hD5/svch3J7KpkcBglPj2eDbfHs+mZ4AbE0YGc//wIKvuqCuC\nSLAJ9ScyimZZ2wryc2PJEyOYGRfBPw9m8nXaVSqqTLvwXMst4/1/nuN/vzxHRKg3MUMCGT0woMOX\npxVBJNgE0VHd/vy9XZn38CCemNyf1B+z+SrtKhlZxYBpT7VzVwo5d6WQ7clnCPJTE9nfj6F9ujGo\nlw+uzo7tWjYRRIJNuCJu7egwLk4OTB4TyuQxoVzLLWX/iWwOnrxhsTJmdp6W7Dwte1IzkcsgNNCD\nAaHe9A3xok+wJ4G+6jadcCqCSLAJF68WAaYNFcVkxo7TM8Cd3z04gCenRKDJKeXITzkcT7/Jpaxi\navuxayTIvF5C5vUSOHQFABcnBaHdPQgNdMdNUX7P5RBBJFhdibaKG7fXIQrr4YGzUvxadjSZTEZY\noAdhgR48Mbk/xWVV/JSRz6mMW5zNLCA7T2txfkWVkfOaQs5rCqkuL7znny/+xQWru3ityPy4f08v\nK5ZEqOXp5sS44T0YN9w0W7tEW8WFa0VculbMpawirtwoobC07WZsiyASrC79al0Q9Qv1tmJJhKZ4\nqJ0YNSCAUQMCzMeKy6q4mlvKyTMZxO+7w4tboEUTGouLi+nfv3+Dr2effdZ8zpYtW4iNjWXYsGHM\nnTuXzMxMi2vo9Xri4+MZO3YskZGRLF68mLy8vHsrvdApXLhaV7XvHyJqRPbC082JoX26MSEq+J6v\n1aIaUXp6OjKZjKSkJFSquvkFnp6eACQkJLB9+3aWLl1KYGAgmzdvZs6cOXz55ZfmbadXrlzJt99+\ny7Jly3B1dWX9+vUsWLCAXbt2iW2nuzBjjWRumnmqnfD3FmsQdUUtCqILFy7g4+NDdHR0g+d0Oh1J\nSUksWrSImTNnAjBixAhiY2PZuXMns2fP5tq1ayQnJ7Nhwwbi4uIA6NevH3FxcaSkpDBx4sQ2fEuC\nPcm6WWbetaNfiJf4o9RFtTiI+vXr1+hzp06doqKigtjYWPMxd3d3oqKiOHDgALNnzyYtLQ2ZTMb4\n8ePN54SEhBAeHk5qaqoIoi4sXVPXLOsnmmU2R5IkKg1V6KrLKddXUF5t+qowVFJRXUl5dSU5N3Lu\n+ee0OIicnJx47LHHOHfuHF5eXvzud79j3rx5XLlimlfQs2dPi9cEBwezb5+pB0uj0eDr64uzs3OD\nczQazT2/CcF+XajXUd0/RHRUt7eamhpKq8ooriw1f5VWlVFSWUZplZayKi1leh1avQ6tvhydvpwa\n6c5bXevbYHnaZoOopqaGy5cv4+rqyosvvkhgYCD79+9nw4YNVFZW4ujoiFKpxMHB8lIqlQqt1jT3\nQKvVWvQt1T8nNzf3nt+EYL/O364RyWUQHuxp5dLYP71BT56ugFztLfJ0+eSXF5GvK6SgooiC8iKK\nK0ubDRZraFGNaOvWrQQGBhIcbOodj4qKQqfTsX37dhYuXNhku14urxuUa8k5LfXoo482OKbX61t9\nHcG6bhVVcP2W6Y9Vrx4euDiJ2SQtVVal5VrJDbJKbnC9NJcbZbncKMujoLyo+Re3gFLhiFqpQqV0\nRa10ReXoiqujC66OLrg4Opv/6+LgjLaghD9yhIULF6JUKi2us2vXrhb9vGb/5eVyOaNHj25wfNy4\ncXz88ce4uLig1+sxGo0oFHX7lOt0OtzcTFP11Wo1Ol3DjeDqnyN0PT9erJu+MbyfnxVLYttKKku5\nXHiVjEINmUVZaIqyKKwobvV1VI4ueLt64e3igYezO17Opv96Orvh7uSGu5MatZMKd6UapYOy+Qve\nlu2Q3eqy/FyzQZSXl8f+/ft54IEH8PKq60ysqjLNqvTw8ECSJLKzswkJCTE/n5WVRVhYGAChoaHk\n5+ej1+stEjMrK4uoqKhWF7qxlM3OzmbChAmtvpZgPScuiCBqzC1dAWduXuD8rQzS8zPI1d5q0euc\nHZzorvbD360bAepudHP1oZvKh24qb3xcvXB2aN/F5hITEwkKCrqr1zYbRHq9nldeeYWKigqeeuop\n8/GvvvqKsLAwJk2axCuvvMLevXuZN28eACUlJRw7doxFixYBEB0djcFgYN++febhe41GQ0ZGhsWk\nSKHrMNZInLpo+oC5OCm6dEe13qDnp7wLnMw5y8ncc9xsJnicHJwI8ehBT88e9PQIJMg9gB7u3fF0\ndrfb6Q/NBlFQUBC/+tWv2LRpEzKZjN69e/Pvf/+bvXv3snnzZlxcXJg1a5b5+ZCQEBITE3F3d2fa\ntGmAaXQsLi6OFStWUFZWhpubGxs3biQiIkLUYrqojKwi82LuQ8K74ejQtVYtrqiu5Ifrp0nLPsHp\n3PNUGRvv41TI5IR6BdPHJ4w+3mH09u5JgNrvrvpWbVmLegffeOMN3nnnHXbs2MGtW7fo3bs3b7/9\ntnle0JIlS1AoFCQlJVFeXk5kZCRr1641z6oGWLNmDfHx8axbtw5JkoiJiWH58uV2m+DCvTlxoe6v\n/vC+3axYko5jrDFyKvc832nS+OHGaaqNDXfVkMlk9PEOY5B/Pwb69aWPT1i7N6lsgUyyldWz71Ft\nH1FKSspdt1OFjvOntw+Yh+63vjSBQF91M6+wXwXlRaRkHiQl8xBFFSUNnlcrVUR2H8SIHoMZ7N8f\ntbJjl2m9V23x2RPjpUKH01ZUc+H2/WUBPq6dNoQuF17liwt7Scs60WDujsrRhdHBkcQEj2CgX18U\nckUTV+kaRBAJHS7tpxxqakwV8chOOFqWfiuDT89+yU830y2Oy5AxrPsAYsNiGBE4GEdF+64DbU9E\nEAkd7sCp6+bHY+ttk2zvMguv8ffTn3P65nmL4yqlKxN7jeWB3uPwU/taqXS2TQSR0KFKtFWcvD1s\n7+3uzIAwHyuX6N4VlBfx0andHLx2zOK4j6sXD/d7gNiwaJwdnZt4tQAiiIQOdqRes2zssEAUbbgT\nREczGA3882IKn537N1WGumVTfVy8+H8Dp/CL0DGi+dVCIoiEDnXgZF2zbJwdN8syCjQkHvsb10rq\n3o+rowuPDogjLnx8q26REEQQCR2osLSSny7nA+Dn5UI/O1wo32A08MnZf5Kc/h/zvvEyZMT2iuGJ\nwb/B3VncO3k3RBAJHea7E9nmvbLGDethd5NZb5TmsiktiStFWeZjwe7dWTjqSfr4hFmxZPZPBJHQ\nIWpqJP59WGP+PnbkvS+43pEOXTtG4rGPzH1BcpmcRwfE8WjEFBwU4mN0r8T/QaFD/Hgxz7yl8aDe\nPoQE2Me20oYaIztO7uSrS/vNxwLU3Vg0Zo6oBbUhEURCh/jXIY358YMx9vEB1lbp2HD4Xc7kXTAf\ni+k5kgUjZ+IihuPblAgiod3dLCzn2HnTksBebk6MGdTdyiVqXk5ZHm+kJpjXAlLI5Dw1fDqTw39h\nd31b9kAEkdDu/n34irmTevKYUJtf8uNy4VXeSE2gtMq0jK1aqeL5+55hoF9fK5es8xJBJLSrUp2e\nfx027fSikMuYPCakmVdY1+nc86w7tJXK253SPdwCePH+3xOg7hpLlViLCCKhXX3+XYZ5A8VfjgzG\n19PFyiVq2okbZ1h/aCvVNQYA+viEsWzc73Fz6pyrA9gSEURCuynV6fnnwUzAVBuaMdF2mzbHrp9i\nw+F3MdaYQnN494E8F/N0l1iUzBaIIBLaTXLqZYvaUICPbS74deLGGYsQiuoxlOei54v5QR1I/J8W\n2kV+cQXJqZcBkNtwbejMzXTWH9pqDqExwZEsHjMXhy6+UFlHs+3hC8Fu/e8/z1GlN324HxjV0yZr\nQxkFGv5yMNHcJxTVY6gIISsRQSS0ubOZBXz3o2nTPZWzA09OibByiRrKKcvjjQPvmG/ZGBowgD9G\nzxMhZCUiiIQ2ZTTWsG33T+bvn5jcHw+1bXX4FleU8D/fvUXZ7XlCfX168cJ9C8TaQVYkgkhoU7v2\nZ5B5w7RTRbC/Gw/eZ1u3c1QZ9Pzl4BbydAXA7XlC4/4LJ7F+kFWJIBLazJUbJfz967oF4//r/w3B\nQWE7v2I1Ug3vHP2Ay4VXAfBy9uDPv/iDmCdkA2znt0Swa9WGGv76fz9iMJru5Xj4/l4M7m1bC8V/\ncuYL0rJPAOCkUPLiuP+im8r+18zuDFoVRHq9nilTpvDSSy9ZHN+yZQuxsbEMGzaMuXPnkpmZ2eB1\n8fHxjB07lsjISBYvXkxeXt69l16wGf/75Vlzk6xHNzW/e3CAlUtk6fC14+w69xVgWlFx0Zg59PK2\n7dtNupJWBVFCQgJXrlxpcGzr1q3Mnz+fjRs3UlZWxpw5c9BqteZzVq5cyZ49e3jhhRdYs2YNFy5c\nYMGCBXSSTWa7vNQfs9mTWjeDeskTkTg52s7ok6Yoiy3f7zB///iQ3zAqaJgVSyT8XIsnNJ47d44P\nP/wQb29v8zGdTkdSUhKLFi1i5syZAIwYMYLY2Fh27tzJ7NmzuXbtGsnJyWzYsIG4uDgA+vXrR1xc\nHCkpKUycOLGN35LQka7mlPLWJyfN3899eCB9bWgt6rIqLW8eTKTKqAdM6wn9pv8kK5dK+LkW1YiM\nRiPLly9n/vz5+PnV7cx58uRJKioqiI2NNR9zd3cnKiqKAwcOAJCWloZMJmP8+PHmc0JCQggPDyc1\nNbWN3oZgDQUlFax+L808cfH+4T349dheVi5VnZqaGjYdSeJWeSEAYZ7B/FfUk2I9IRvUoiDatm0b\nBoOBBQsWWBzXaDQA9OzZ0+J4cHCw+TmNRoOvry/Ozs5NniPYn/LKal7dfpRbRRUAhHZ35w/Th9nU\nh/zjM1+Yd11VK1U8P3aBGKa3Uc02zS5fvszWrVvZsWMHDg6Wp+t0OpRKZYPjKpXK3Eek1WpRqRpO\n71epVOTm5t5L2QUrqao2Ev+/35s7p308nFk5fwwuTrZz6+IP10+x+3xd5/Sz0XPxEyNkNuuOvzmS\nJPHyyy8zffp0hgwZ0ujzTf0FlMvrKlstOac1Hn300QbH9Hr9XV1LaJ2qaiOvJx3l1CXT/mSuzg6s\nejraptYZyi3LI+HoB+bvfzv41wwNsK1RvM5o4cKFKJWWNc5du3a16LV3DKIdO3aQm5vLu+++i9Fo\ntBjlMhqNqNVq9Ho9RqMRhaJulESn0+HmZtpoTq1Wo9PpGly7/jmCfaisMhD/v9+b9653UipYMXc0\nod1tZ0eOKoOe9Ye2UV5tajKODBzC1IjJVi6V0Jw7BtHevXvJzc1l5MiRFsfT09P5/PPPefXVV5Ek\niezsbEJC6uZkZGVlERZmmtofGhpKfn4+er3eIi2zsrKIioq6q0I3lrLZ2dlMmDDhrq4nNK9EW8Vr\n7x3lwrUiwBRCK+eNYZANTVqUJIl3j/+dq7e3gfZX+fLfo59CLhPzdjtCYmIiQUFBd/XaO/4Lvfba\na+zcuZPPPvvM/BUaGkpsbCyfffYZU6ZMQalUsnfvXvNrSkpKOHbsGNHR0QBER0djMBjYt2+f+RyN\nRkNGRgYxMTF3VWihY+UW6Hgx4aA5hFycTCE0ONx2Qgjgm8uppGqOAuCocOT5+55BpXS1cqmElrhj\njSg0NLTBMWdnZzw9PRkwwNTmnjVrFps2bUImkxESEkJiYiLu7u5MmzYNMI2OxcXFsWLFCsrKynBz\nc2Pjxo1ERESIGowdOHXpFn/ZcYyy8moAPNVOrJw/hvBgTyuXzNLF/Eze//FT8/fPjHiCUC/72k22\nK2v1MIdMJrPofF6yZAkKhYKkpCTKy8uJjIxk7dq1qNV1NxKuWbOG+Ph41q1bhyRJxMTEsHz5cpsa\n6hUsSZJEcupl3v/nOWpqTH2D3X1UrHpmDIG+tnWTaFFFCesPbzOvsjip9/38ImyMlUsltIZM6iT3\nWdT2EaWkpNx1O1UwKSvXs+kfP3L0bN30imF9uvGn343EzdW25uHojdWs3reBS4UawLTzxurYJWK9\n6Q7UFp898a8lWDiRnsemj3+ksLTSfOw39/dmzkMDUNjQkh5wu3P6h7+bQ8jLxYPn73tGhJAdEv9i\nAgC6imo++PIc/z6iMR9TuTjyx8eG2+wW0Z+d+xffadIAcJQ7sPS+hXi72FbfldAyIoi6OEmSSDuT\nQ+KunyxqQYN6+/DcY5H4edvmqNP+K0f45Mw/zd8vjHqScJ9Q6xVIuCciiLqw67e0bPv8J06k160N\npXSQ8+SDA3h4XC/kctscTDh2/RRbj/3N/P1jgx9mXOgoK5ZIuFciiLqgUp2eT/Ze5MtDmeYVFQGG\n9vHl99OG2tyoWH0/XD9t2gxRqgFgQq+xPBIRZ+VSCfdKBFEXUlll4J+HrrAz5SK6SoP5uLe7E7Mf\nGsj4yCCbnlLxffZJNh7ZbrEZ4vwRj9l0mYWWEUHUBVRWGfgq7Sqf7btEsbbKfNzRQc7D43oxY2Jf\nXJ1tdysdSZL48uI+Pjz5GRKmGtyYINOOrAqxD1mnIIKoE9NWVPPvw1dITr1MidZydYLxkUE8OSXC\nZjuja+kNet7/8VNSMg+aj8UEj+APY+aIzRA7ERFEnVBeUTlfHMjk67SrVFQZLJ4bMyiAJyb3JyzQ\nw0qla7nskhw2HtlOVskN87GH+0/iiSG/ETeydjIiiDoJSZI4m1nAPw9e4ciZHPNtGbWiB3dnxsS+\nhAfZ/jwbvUHP7vNfk5z+Hwy396V3lDswJ3IGE3uPs3LphPYggsjO6Sqq2X8im38fvsLV3DKL5xwU\nMsZHBvPI+N70DLCdNYOaUlNTw6FrP/DxmT3mnVgBAt38+WP0fEK9xK07nZUIIjtUW/v55vtrHDx1\nA3210eJ5N1clU2JCeTAmFB8P21k5sSkGo4HDWcdJPv81WaU55uNymZwH+8QyY/CvcXZwsmIJhfYm\ngsiO5Bbo+PZ4Nt8ezyInv+Gql+HBnvwqJpRxw4Nsal+xpuTpCth/5Qgplw9SVFli8VxEt3DmDP+t\nqAV1ESKIbFyJtoqDp27w3YlszmsKGzzv4qTg/uFBTB4TQp9g29lPrCmlVVqOZv3I4awfOJd3yTwc\nX6uXV08eG/wwQwMGiPlBXYgIIhtUVq4n7accDpy8zqmM/AYdzwADe/kwMaon9w0NtKndM36uRqrh\nWvENTuWe4/iN01woyGx0h9/I7oN4qN9EBvr1FQHUBdnub3AXU1xWxdGzORw+ncOpS7cwNhI+3X1V\nxI4IJnZEEAE+DbdosgU1NTVkld4g/dZlzt+6xNlblyipLG303G6u3vwibAzjw2LEVj9dnAgiK7pZ\nWE7amRyO/JTD+SsFNJI9+Hg4M3ZoD34R2YPwIE+bqi3USDXc0hVwpSiLzKJrXC7UkFFwlQpDZZOv\n8Vf5MrLHUKKDI+njE2ZT70ewHhFEHaimRuLy9WKOns3l6JlcNDmN1xS83Z2IGRLI2KE9iAj1tvpd\n8JIkUVhRzPXSXK6X5pJVmkNW8XWuld6gorrp0AFQKhyJ6NaHIf4RDOs+gCD37iJ8hAZEELWzyioD\npy7d4tj5mxw7l0thaVWj5/l5uxIzuDvRg7vTP6Tjw6faWE1BeRF5uoLbX/nklt0iV5tHTlkeVcaW\nbWDp5qSmr08Y/Xx7E9EtnN5eIWLFRKFZ4jekHeQW6Pjh/E2OnbvJT5fzqTbUNHper0APRg8KYMyg\n7oQFurdLTUGSJHT6coorSymqLKGowvRVWFFMYUUxBeVFFJQXUVxZ2mAEqzmezu6EeAYR6hlEL++e\n9PYOpZurt6jxCK3W6YKoRmr8Q9+eqg1GzmYW8MP5PH44f5Prt7SNnuegkDO4tw+jBwYQNSCgVTec\nSpJEtbGa8uoKyqsr0OrL0VWXo9OXo9WXU1alpUyvQ1ulo0yvpbRSS0lVGSVVZeZlM+6Wr6s3gW7+\nBLr508M9gCCP7gR7BOLuZLvrFgn2pdMF0eIvX8G9myeuji64Orrg4uCEs6MTzg7OODk44axQ4uTg\nhFLhePtLiaPCAUe5Aw5yBxwVDijkChQyBQ5yBXKZ/PaX7PaNljJkMigsreTi1SLSrxaRkV1kmt0s\nA5CQu0sgM32pXBX07uFOWJAbPbq5IlMYMNRoOHIzA/0NPXpjtfmryqinylBFpaGKKoOeCkMlldWV\nVBiqqDBU3nOgNMVNqcLH1YtuKh+6qXzwU/ngp/IlQN0NP5UPSgfb2rlD6Hw6XRABVN7+MBdWFLf/\nD3MBWR9o6gYEA3ABuHAduN7+xanPTanCw9kdD2c3PJzd8XR2x9vFA09nD7xdPPF29cTHxQsnETSC\nlXW6IOrp0QNJraCiuoLy6krz3dv2SC6T367ROd+u3TmjUrrg4uiCytEFldIVtdIVtVJl/nJzUuHu\npEatVIlFwwS70aIgqq6uJiEhgS+++IKioiKGDBnCiy++aN52GmDLli188sknFBUVERkZycsvv0yv\nXr3Mz+v1etatW8e//vUvysvLGTt2LC+//DJ+fn5t+ob+NO6/LDZ5MxgNpiZOvSaP3qinyni7WWSo\nprrGQLWxmipDNTfyS7mWV8r1/FJKy6sAUxNLJrvdkXv7vx5qJ7p5ueDn5YKXmwvy2zvgypEhlyuQ\ny2QoZAoUclMTz+F2c89RYWoCKhWOOCoccZQ74uRgaiI6OShxUihxdnDCycEJR7mD6PgVuoQWBVF8\nfDxffPEFS5cupWfPnnzwwQf87ne/44svvqB79+4kJCSwfft2li5dSmBgIJs3b2bOnDl8+eWX5q2n\nV65cybfffsuyZctwdXVl/fr1LFiwgF27drXrh81B4YCbQo1bEx2rRaWV/HD+JsfP3+TkRS0VVQrA\n6/ZXHTdXR4b382NEf38i+/nh6SbuBheEttJsEGm1Wnbu3MkLL7zAb3/7WwAiIyMZPXo0ycnJPPnk\nkyQlJbFo0SJmzpwJwIgRI4iNjWXnzp3Mnj2ba9eukZyczIYNG4iLM+240K9fP+Li4khJSWHixInt\n+BYtSZJE5vUSvj9nmtdzKavpfqTwIA9GRPgzsr8/fXp6obDR7XUEwd41G0QuLi58+umn9OjRw3xM\noVAgk8nQ6/WcOnWKiooKYmNjzc+7u7sTFRXFgQMHmD17NmlpachkMsaPH28+JyQkhPDwcFJTU9s9\niKoNNfx0OZ+jZ3L4/mwu+SWNzwZWOTuYaz0jIvzwcnNu13IJgmDSbBApFAr69+8PmGoT2dnZvP32\n28hkMh5++GEOHToEQM+ePS1eFxwczL59+wDQaDT4+vri7Ozc4ByNRtMW76OB8spqjqfnkXYmhx/O\n36S8svFO6yA/NaMGBDBygD8Rod442Nj+7oLQFbRq1Oydd94hISEBmUzG4sWLCQ0N5T//+Q9KpRIH\nB8tLqVQqtFrTxD6tVotK1fBucZVKRW5u7j0Uv6EjP93gwn+uc+JCXqMzmhVyGQN7+TBqYACjBgTQ\n3dc272IXhK6kVUE0adIkxowZQ1paGu+88w56vR5nZ+cmO5vl8rraRUvOaQvv7TmLo6u3xTEXJwdG\n9PdjzKDujIjwR+1iu3t4CUJX1Kog6tu3LwAjR45Ep9ORlJTE888/j16vx2g0olDUzVvR6XS4ubkB\noFar0ekaLm1a/5zWePTRRxscq6oy3UxafXsSo9rVkWF9ujGivz8RYXVNruKCm3TANEdB6DJqWzUL\nFizAyclyNHnXrl0tukazQZSfn09qaipxcXG4utbdGxUREYFer8fDw8PcdxQSEmJ+Pisri7CwMABC\nQ0PJz89Hr9ejVCotzomKimpRQZtz9epVALKPJJqP/QR82CZXFwShOQaDoUEQtVSzQVRaWsqf//xn\nZDIZjzzyiPn4wYMH8fHxYeLEiSiVSvbu3cu8efMAKCkp4dixYyxatAiA6OhoDAYD+/btMw/fazQa\nMjIyePbZZ1td6MZSdurUqVRUVLBt2zaLmpktW7hwIQCJiYnNnGkb7K28IMrcERYsWIDBYCA5ObnB\ngFRLNRtEvXr1YvLkyaxZswa9Xk9wcDBff/01X3zxBW+88QYqlYpZs2axadMmZDIZISEhJCYm4u7u\nzrRp0wDT6FhcXBwrVqygrKwMNzc3Nm7cSEREBBMmTLirgv+cXC5HpVJZ1MpsXW3tsP5McFtmb+UF\nUeaO4OTkhJOT012HELSwj2jt2rUkJCSwbds2bt26RXh4OG+99RYPPPAAAEuWLEGhUJCUlER5eTmR\nkZGsXbvWPKsaYM2aNcTHx7Nu3TokSSImJobly5eLWxgEQWhZEDk5OfH888/z/PPPN/q8QqFgyZIl\nLFmypMlrODs78+qrr/Lqq6/eXUkFQei0xOw9QRCsTiY1tsmUIAhCBxI1IkEQrE4EkSAIVieCSBAE\nqxNBJAiC1YkgEgTB6kQQCYJgdSKIBEGwOhFEgiBYnQgiQRCsTgSRIAhW1ymC6JNPPmHy5MkMHTqU\nxx57jJMnT1q7SHdUU1PD+++/z4MPPsjw4cP51a9+xUcffWTtYrWIXq9nypQpvPTSS9YuSrOOHDnC\njBkzGDp0KL/85S95++23seU7mmpqanj33XeZNGkSw4cPZ8aMGaSlpVm7WI1KSUkhMjKywfEtW7YQ\nGxvLsGHDmDt3LpmZmS26nt0H0e7du1m1ahW/+c1vePvtt3F3d2f+/Plcv97BG823wjvvvMNf//pX\npk6dypYtW5gyZQrx8fG899571i5asxISErhy5Yq1i9Gs48eP8/TTTxMeHs62bduYNWsW7777Lps3\nb0Xo9QoAAAXRSURBVLZ20Zq0fft2/vrXvzJt2jQ2b95McHAw8+fPJz093dpFs3DixAn+9Kc/NTie\nkJDA1q1bmT9/Phs3bqSsrIw5c+aYN9G4I8nOxcbGSqtXrzZ/X11dLU2YMEF6/fXXrViqphmNRiky\nMlJ66623LI6vXr1aiomJsVKpWubs2bPSsGHDpOjoaGnZsmXWLs4dPfHEE9LChQstjq1fv1568skn\nrVSi5k2ZMsXi/6vRaJTGjx8vvfbaa1YsVZ2qqipp27Zt0qBBg6RRo0ZJw4cPNz+n1Wql4cOHS9u3\nbzcfKykpkSIjI6X333+/2WvbdY3o6tWr3Lhxw2JzRwcHB8aPH8+BAwesWLKmabVaHnnkEfOicrXC\nwsIoLCyksrLxzR+tzWg0snz5cubPn4+fn5+1i3NHhYWFnDhxwrwzca0lS5awY8cOK5WqeXq93mLb\nLblcjlqtprjYNrZ7SE1NZfv27SxbtoxZs2ZZPNfcRqvNsesg0mg05uVp6wsKCiIrK8sm+wPc3d15\n+eWXzZtW1tq3bx8BAQH3tNxme9q2bRsGg4EFCxZYuyjNunjxImBa0G/hwoUMGTKEmJgYEhISbPJ3\notbMmTNJTk7myJEjaLVaPvjgAy5fvsxDDz1k7aIBMGTIEFJSUpg5c2aDlVVrm+uNbbTakk1UW7Wd\nkK2pbXv+fPNGlUpFTU0N5eXljW7saGs+/fRTjhw5wooVK6xdlEZdvnyZrVu3smPHjgYbadqiwsJC\nJEli2bJlPPTQQ8ydO5fvv/+eLVu24OzszPz5861dxEY9/vjjpKWlMWfOHMC0F+Af//hHi63arelO\nNWGdTtfsRqt3Yvu/VXdQ+9etozZvbA979uxh1apVTJkyhZkzZ1q7OA1IksTLL7/M9OnTGTJkiLWL\n0yIGg2l78XHjxrF06VIARo0aRVFREVu2bGHevHk2uVZ67SjT6tWr6dWrF4cPH+btt9/Gzc2NJ554\nwtrFuyNJku7pc2jXQVS7OaNOp8Pbu253V51Oh0KhwMXFxVpFa5H333+ftWvXMnHiRN58801rF6dR\nO3bsIDc3l3fffRej0WjRtPn5ppq2onb/vbFjx1ocj4mJ4e9//zvZ2dkEBwdbo2hNOn78OCdOnOCt\nt95i0qRJAERFRWEwGHjzzTd55JFHbPr3Wa1WN7vR6p3YfpXhDkJCQpAkiaysLIvj2dnZhIaGWqdQ\nLbRhwwb+8pe/MHXqVDZt2mSzTZ69e/eSm5vLyJEjGThwIIMGDeLChQvs3r2bQYMGcePGDWsXsYHa\nPsPq6mqL47U1JVusDeXm5iKTyRg6dKjF8REjRlBZWWnT01HAtImqdHuj1frqb7R6J3YdRKGhoXTv\n3p29e/eaj1VXV7N//36io6OtWLI7++CDD9i2bRuzZ8/mjTfesOkm5GuvvcbOnTv57LPPzF8hISHE\nxsby2Wef2eQIWnh4OP7+/nz11VcWx/fv34+fn59N7hdW+0E+fvy4xfGTJ0+iUCgICAiwUslaZvjw\n4eaNVmvVbrTaks+ibf4ZboWnn36a119/HTc3NyIjI/nb3/5GcXExTz31lLWL1qhbt26xfv16+vXr\nx5QpUzh16pTF84MHD7apYGqsZuns7IynpycDBgzo+AK1gEwm47nnnuOll15i1apVTJ48mcOHD5Oc\nnMzq1autXbxGDRw4kPHjx7N69WqKi4vp3bs3R48eZfv27Tz11FMWewTaIldX12Y3Wr0Tuw+iJ554\nAr1ez44dO9ixYwf9+/cnKSnJJv/qgWmr7urqai5evMhjjz3W4PkjR47g6elphZK1nEwms8nmTX1T\np05FqVSSmJjI7t27CQgIYPXq1UyfPt3aRWvSW2+9xcaNG9m6dSslJSWEhITwyiuvMGPGDGsXrVE/\n/x1oyUarTV5LsuWJFYIgdAm20wYQBKHLEkEkCILViSASBMHqRBAJwv9vp44FAAAAAAb5Ww9jT0HE\nTkTATkTATkTATkTATkTALn4bFnnMbDtnAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x103837f50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARkAAAEVCAYAAADU0pFsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4lFXe+P/39JnMZBLSeyG0UEJCNRQFURRXV2AVUfQS\nFZXvrj62nyu7uvuIrrqrjxWUIliw7K4iiu66roIiIKEGEBEIhPSQXmeS6fP7Y8KEGEidyUzCeV3X\nXLnnzF0+0cyHc5/7FInT6XQiCILgJVJfByAIwsAmkowgCF4lkowgCF4lkowgCF4lkowgCF4lkowg\nCF7V7SSzdetWxo0b16581apVzJw5k/T0dO68805Onz7tkQAFQejfupVksrOz+f3vf9+ufOXKlaxZ\ns4YlS5bw8ssv09jYyB133IHBYPBIkPPnz2f+/PkeOVdf6W8x97d4QcTcFzwRb5eSjMVi4c033+T2\n229HLpe3+cxoNPLWW29x//33s2jRImbOnMn69esxGAxs3LixV8EJgtD/dSnJbN++nXXr1rFs2TJu\nvfXWNp8dPnyY5uZmZs6c6S7T6/VMnDiRHTt2eDZaQRD6nS4lmbS0NLZu3cqiRYuQSCRtPsvLywMg\nISGhTXl8fDz5+fmeiVIQhH5L3vkuEBERccHPjEYjSqWy3W2UVqv1WJuMIAj9V5eSTEecTme72s1Z\nUmn3n5Cfr5Hp5MmTxMTEUFBQgEwm6/Y5fcFisQBQXFzs40i6pr/FCyLmvmA2mykoKGDu3Lntvs+b\nNm3q0jl6nWR0Oh0WiwW73d4mARiNRgIDA3t7egBiYmLIz89n9uzZHjlfX5o1a5avQ+iW/hYviJj7\nQnNzM1qttkfH9jrJJCUl4XQ6KS4uJjEx0V1eVFREcnJyt893vuxYUFDA7NmzSb13MgT0vP+gXhVI\npC6MSF04US2v6MAI9KrAC9bGBOFiVlZWxqJFi1i7dm2b73d39DrJZGRkoFQq2bJlC3fddRcA9fX1\n7Nu3j/vvv7+3pwdw15BWLniWsKhwapvrqWmuo7qpluqmWqqaaqhqqqHSWEOlsRqz3XLe85iwUeAo\no6ChDBpay7UKDfFBMSQEx5IUHEdicByJQbEo5UqPxC8I/V1vmil6nWQCAgK49dZbefXVV5FIJCQm\nJrJ69Wr0ej033HBDb0/fjlquIjowgujA8zdGO51OGsyNlBuqXC9jJWcaKyhrrKC0sRyjtbndMUZr\nM8ercjlelesuk0gkxAVGkRySQMqgRFJCEkkaFI9SpvD47yQIA1mPkswvby0efvhhZDIZb731Fk1N\nTYwbN47nn38enU7nkSC7G1uQWk+QWs+wsMFtPjubgEoayilpKKOooZTi+jMU1ZdSb25st29RwxmK\nGs6wPX8PADKJlKTgeIaGJjMsLJkRYUMI04b02e8mCP2RpD9Mv1lcXMysWbPYunUrcXFxXrlGvamB\ngroSCupKyK8rIq+2iJLGMjr7zxMaMIjU8KGMDB/KqIhhROnCRfuOMGB44rvX69ulgSJIrSctSk9a\nVKq7zGQzk19bzOnaAk7VFHCqOo8yQ2Wb46qbatlZsJedBXsBCNEEMzpiOGlRqYyJHMEgTVCf/h6C\n4G9EkumAWq5iRHgKI8JT3GUNZgM5VafJqT7N8cpTnKopwOawuT+vaa5je8Eethe4brESg+NIjxpJ\nRvRohocNRibtH/18BMFTRJLpJr1Kx4TYNCbEpgFgsVs5VZ3Pz5U5HK3IIafqNNZzkk5BXTEFdcVs\nPv41WoWG9OhRTIwdS3r0KAIUGl/9GoLQZ0SS6SWlTMHIiKGMjBjKDaN+hcVm4XhVLj+WH+Nw2TEK\n6lp7dhqtzfxQuJ8fCvcjl8pJi0rlkrgMJsaORasM8OFvIQjeI5KMhynlStKiUkmLSuXWsa7bp0Nn\nfib7zBEOlx3DbDMDYHPYyC49QnbpEeRSOelRI5mWOJEJMWmif44woIgk42UhmmAuHzyFywdPwWK3\ncrTiBPuKD7O35BANZtcAUpvDxv7SH9lf+iMauZrJ8RnMTM5kRNgQ8aRK6PdEkulDSpmCjOjRZESP\nZsn4mzlWdYqswgPsLs52J5xmm4lteVlsy8siShfOrMHTmJF8CUFqvY+jF4SeEUnGR6RSKaMihjEq\nYhh3jFvAkfIT7CjYw97iQ+5hEWWGSj748VP+8dPnXBKXwdVDZzAsdLCo3Qj9ikgyfkAmlZEePZL0\n6JE0W03sLsrmu7xd7mEOdofd3WA8eFAC1w6fxSXx45GLx+FCPyCSjJ/RKNTMHDyFmYOnUNxwhi25\nO/k+L8s95up0bSGv7X6bD378jOuGX8GswdNQiYZiwY+JdZf8WJw+msUZN7L613/lngmLiA+KcX9W\n3VTLOwc/5rf/epzNx77GZDX5MFJBuDBRk+kHVHIlV6RMY9bgqRwpP87nx7/hx/JjADSaDXzw46d8\nfuIb5qVexewhl4mR4gPEbbfdxr59+877WVhYGDt37uzjiHpGJJl+RCKRuPvgnK4p4NNj/2VP8UHA\nlWw2HPqEf+d8y02jr+PSpMlIJaKi2t+NHz+exx57rF25QtF//iERSaafGhySyCNT76GwroSPj/7b\nnWyqm2p5Y+8G/nPyO25Pv5GREUN9HKnQG4GBgaSlpfk6jF4R/9T1cwnBsTwy9R6eu3IZYyJHuMvz\naot48ruXWLH7beqa630YoXCxE0lmgEgJSeRPMx7gj5feR5w+2l2+o2AvD/5nOVtyd3Y6N47gn+x2\ne7tXfyJulwaY9OhRjIkcwX9Pfc8/f/qCZquJJmsza/d/wA+F+1g68VYideG+DrPP7DhYwgf/PUaz\n2db5zl6iUclZdFUq0zNiu33stm3bGDVqVJsyiURCVlYWwcHBngrRq0SSGYBkUhnXDLucKfHjeffQ\nRn4o3A/A0YocHv3vMyzOWMDM5MyLoufwpm0nKak0+jgKM5u+P9WjJDNhwgT++Mc/tquF6vUXHmay\nZcsWiouLWbx4cbev5w0iyQxgwZogHsi8i2mJk3hz/4fUNNdhsplZve89ss8c4f9NvG3ATzExf8ZQ\nv6jJzJ8xpEfH6nQ6Ro4c2eX9zWYzq1atYvLkyT26njeIJHMRGB8zhhFX/4n12f90TxO6t/gQ+bVF\nPDTlblJCeraeTn8wPSO2RzWI/urjjz9m7ty5HDlyxNehuImG34uEVhnA/1xyBw9m3oVGoQagwljN\nn7b+H9vysnwcneAJRUVFaLVaEhMTqaqq8nU4biLJXGSmJEzgb7P/SPKgeMA1l80bezfw7sGNOBwO\nH0cn9MbGjRuZO3cuYWFhVFZWdn5AHxFJ5iIUpQvn6VmPMiM5013275ytvLhrLRbb+VffFHyjq43z\n33//PQ6Hg82bN7Nv3z6/qsmINpmLlFKm4P9NvI2k4Dg2HPoEh9PBvpLDPLXtVR6b/v8IVPX9wnxC\nW++9916X9jObzZw4cYJHHnkEcC1M+MILL2C323u1vKyniJrMRUwikXDNsMv5/bSlqGSu6SJyqk/z\n5HcvU2dq6ORowR/s37+fO++8E7PZ7C778MMPkclkrF+/3i9ugcUKkgIAp6rz+euO193TgMYGRvGn\nmQ8QoukfHb4E7/DEd0/UZAQAhoQmsfzyRxikdq14WdJYxvLvXqbB1NjJkYLQMZFkBLdYfRRPXv4w\noZpBAJxprOCZ7StoapmVTxB6QiQZoY3owAj+9/KHCG5ZHSGvtojnd6zCarf6ODKhvxJJRmgnShfO\n45fdj7ZlGd2fK0+yZv8HYhS30CMiyQjnlRgcx2PTf4dC6urlsD1/D58f/8bHUQn9kUgywgWNCE9h\n6cTb3O8//PEzskt/8mFEQn/ksSTjcDh48803mT17NhkZGSxYsIDdu3d76vSCj0xPmsT8kVcD4MTJ\n63veobqp1sdRCf2Jx5LMunXreOWVV7jhhht44403iI+PZ8mSJRw/ftxTlxB8ZMHo68iIdk2c1Ggx\n8mrWeuyO/jU7m+A7Hksyn332Gb/+9a+55557yMzM5IUXXiA8PJyNGzd66hKCj0glUn436XYGaVx9\naI5X5bLx6Jc+jkroLzyWZCwWC1qttvXEUik6nY66ujpPXULwIb06kAcuucs9YO/TY1+RV1vk46iE\n/sBjSWbRokVs3ryZrKwsDAYD7777Lrm5uVx77bWeuoTgYyMjhjJ3xFUAOJwO3ti7AZu4bRI64bFR\n2DfffDO7d+/mjjvuAFyD7x588EFmzJjhqUsIfuA3o65hb8khShrKKKgr5rNj/+WGUdf4OizBj3ks\nydx5552cPn2a5cuXM3jwYHbt2sWKFSsIDAzklltu6fJ55s+f367MYhFznPiLs1NE/Onb/8PpdLLp\n5/8wPXHiRbUCQl/xp2Vqly5dilKpbFO2adOmLh3rkSRz4MABsrOzee2115g9ezYAEydOxGaz8cIL\nLzBv3jw0Go0nLiX4gWFhg5kzZAZfnvwOm8PGe4c38f9NvdfXYQ1IYpnaFmVlZUgkEsaOHdumfPz4\n8axbt46SkhKGDOnabO3ny45nh5sL/uPG0deyo3AfjWYDe4sPcbQih1ERw3wd1oDjL8vUrl692rdT\nPSQlJeF0Ojlw4ECb8kOHDiGTyYiKivLEZQQ/olUGcNPo69zv38n+yC8mSBL8j0dqMqNGjWLGjBks\nX76curo6UlJS2LNnD+vWreP2229HpxNTOQ5EswZP5etT2ymsL6GgvoQfCvczPWmSr8MacM63LK0/\nTKvZVR5r+H3ttdd4+eWXWbNmDfX19SQmJvLnP/+ZBQsWeOoSgp+RSWXcOnY+z25fAcAnP3/J1IQJ\nSKX+MyRuV+F+/vnTF5is5s539hK1QsVNo69jSsKEbh8rlqk9h1Kp5LHHHjtvI5UwcI2NSmVoaDIn\nq/MobSxnV9EBpiVO9HVYbp8f/4YzjRW+DcIEXxzf0qMk09Vlanfu3MnTTz/NFVdcQXJyMiaTiezs\nbJ5++uk2nWR9QaxWIPSKRCLhhlHX8Nz21wFXbWZKwnikEv+ozfx6xJV+UZO5bsSVPTq2q8vUTps2\nDYCHHnoIudz1td63bx/79u3zeV81kWSEXkuPGkVKSCK5NQWUNJSxp/ggmfHjfR0W4FrMric1iP6m\nqKiI+Ph4d4I5ceIEVVVVTJzo+1qlf/xzI/RrrtrMr9zv/5PznQ+juTgdOHCAqKgodu7cyb/+9S9W\nrlzJmjVrfH6rBCLJCB6SET2KaF0E4BqlnS8GT/apAwcOMH/+fKZNm8a1117L1KlTeeaZZ3wdFiCS\njOAhUomUq4Ze5n7/1anvfRjNwNHVZWoPHTrE6NGj3e/DwsL4+eefvRVWt4gkI3jMZUmXoJKrANhZ\nsBeD2ejjiPq39957j1WrVnW6X01NDUql0j22yOFw8P7773P99dd7O8QuEQ2/gsdolQFcmjiJb3J3\nYLFb+S4vi+tGXOHrsAa0H3/8kY8//hiADz74AIfDQU5ODjNmzGDx4sW+Da6FSDKCR1015DK+yd0B\nwLenf+Da4bO6XOUXui8tLc0vxjZ1RNwuCR6VEBzL0JAkwLXUbX5dsW8DEnxOJBnB46Ylto5f2lmw\n14eRCP5AJBnB4zITxrtvkX4o3C9GZ1/kRJIRPC5YrSctcgQANc11HKs65eOIBF8SSUbwiqkJrd3Z\ndxacfwpJ4eIgkozgFZPi0lHIXFNE7i46IFY1uIiJJCN4RYBC41510mht5kRVro8jEnxFJBnBa8ZF\nj3FvZ5ce8WEkgi+JJCN4zbjo1hndss/85MNIBF8SSUbwmmBNEIMHJQBQ0lBGhaHKxxEJviCSjOBV\n42JaRwaL2szFSYxdErxqXPQYNh79EoCDZ37i6qEzfBtQP+JPK0j2hkgyglcNDkkgSBVIvbmRnypy\nMNssqOTKzg8UALGCpCB0SiqRkh49iu/zd2O1W/m5MoeM6NGdHygA/rOCZG+INhnB69IiU93bxyrF\nEIOLjajJCF6XGt66DrpIMt0nVpAUhE6EaUMIDwihsqmGUzX5WGwWlKJdpkvECpKC0EWp4UOpLNiD\n3WHnZE0+oyKG9cl1K3f8QOGH/8De3Nwn1zsfmUZDwi0LCZ8+tdvHdnUFSYfDwUcffURdXR2hoaHY\nbDbCw8MxGo0+n+tXJBmhT6SGD2F7wR7AdcvUV0mm5NPNmEpL++RaF2KtraX0s809SjJdWUHSarVy\n//33M2/ePBYuXOguv+GGG/jb3/7W7Wt6mkgyQp9o2y5zss+uGzvver+oycTO815t4vnnnycqKoqr\nrrqqTfmMGTNISUnx2nW7SiQZoU9EB0a6+8vkVJ3G5rAjl3q/8TJ8+tQe1SD6i5qaGv7xj3/wzTff\ntPvsrrvu8kFE7YlH2EKfkEgkpIYPBcBst5BXW+jjiAaG7OxsEhISiIqKaveZRqPxQUTtiSQj9Bnx\nKLv7OltORiKRtGsEBtd6TGVlZd4Kq1s8mmSysrJYsGABY8eO5fLLL2fFihXtWsWFi9fwsNb2gVM1\n+b4LpJ/oygqSkydPpqqqisrKSndZZWUlpaWl563d+ILH2mQOHDjA3Xffza9//WseeeQRjh49yiuv\nvIJUKuV3v/udpy4j9GMJQTHIpXJsDht5NeJ2yRN0Oh0rVqzgxRdfZPTo0QQHB6PT6bj66qt9HZqb\nx5LMSy+9xPTp03n22WcBV4atq6tjz549IskIAMhlchKDYsmtLaDcWIXBYkSn1Po6rH5vxIgR/PWv\nf/V1GBfkkSRTU1NDdnZ2u6rdww8/7InTCwNIckgCubUFAOTVFjGmZekUYeDySJtMTk4OACqViqVL\nl5KWlsaUKVNYuXKlaJMR2jg7Ux5Abk2BDyMR+opHkkxNTQ1Op5Nly5aRkpLCunXruOWWW1i1ahXr\n16/3xCWEASIlJNG9fVo8xr4oeOR2yWazATB9+nQeffRRACZNmkRtbS2rVq3irrvu6vRRnHBxiNdH\ni8bfi4xHkkxAQAAA06ZNa1M+ZcoUPvzwQ4qLi4mPj+/SuebPn9+uzGKx9D5IwS+Ixt/+aenSpSiV\nbUfOb9q0qUvHeuR2KTHRVQW2Wq1tys/WcEQtRjhXckhru8xpUZsZ8DxSkxkyZAiRkZF89dVXXHfd\nde7ybdu2ERERQVxcXJfPdb7sWFxczKxZszwRquAHUgYlsKVl+3RtIWlRqR3uL/je6tWru/U9PpdH\najISiYSHHnqIb7/9lieffJKsrCxefPFFNm/ezH333eeJSwgDyGDR+HtR8VhnvLlz56JUKlm9ejWf\nfvopUVFRLF++nBtvvNFTlxAGiHMbfwvqin0djuBlHp3q4ZprruGaa67x5CmFAUgukxMbGElBfQll\nhkoxHecAJ0ZhCz4RHxQDgNPppKSx3MfRCN4kkozgE2eTDEBRvW+nxxS8S8yMJ/hEfFC0e7uwvsSH\nkfgvsUytIPSCqMl0jVimVhB6KFwbikqmxGy3UCiSzAWJZWoFoYekEilxLbdM1U21NFl8t5qA4F0i\nyQg+0+aWqUHUZi7Ebre3e/Un4nZJ8JmEX7TLnDsHsKccPVjCtv+ewGy2efzcXaVSyZlx1XBGZcR2\n+1ixTK0g9MK5NRlvtcvs2pZLdaXRK+fuKgNmsr7P7VGS6eoytVlZWTz99NNkZmaSmppKQ0MDJ0+e\n5NFHHyUkJKRX8feWSDKCz/TFE6YpM1L8oiaTOWNI5zueR1eWqQXIzMzEZDJx//33u2s4+/bt46mn\nnuKVV17p0bU9RSQZwWcGqYPQKgMwWpq8lmRGZcT2qAbR35w5cwa1Wt3mFiotLY17773Xh1G5iIZf\nwWckEom7XabBbKDe1ODjiPqv/fv3M2HChDZllZWVyOW+r0eIJCP4VGxg6wJkpWIMU4/t37+f8ePH\ntynbu3cv48aN81FErUSSEXwqRt+aZEoa/GNZVX/S1VklDxw40K4m88UXX3DHHXd4I6xu8X1dSrio\nxeoj3dslDaImc6733nuvS/vV1dVhMBiIjW1te/roo49IT09n8uTJ3gqvy/pVktl3tJyyRgUalZwA\ntbzlpwKtWo5MJipl/VGsvnWgZGmjqMl017Fjx/j73/9OYGAgGzduxOFwUFFRQXJyMg888ICvwwP6\nWZJZ89mPKALOP5OaSilDq5aj1SjRaRRoNQoCAxQEBijRBSjd23qtkkCt62eQToVKIevj30I4V1jA\nIJQyBRa7lWJxu9RtqampPPXUU74Oo0P9Ksl0xGyxY7bYqWkwd+s4jUpGkE5FkE5FsE7FIL2aYJ2K\nEL1rO6TlNShQJWpLXiCVSIkJjCS/rpgqYw1mmwWVmCVvQOlXSeY3lw9FrQul2WyjyWSlyeT6aTTZ\naGq2YjRZMTZbsdm7vjRus9lOs7mJsuqmDveTSiA4UEVokIawYNcrPFhD+CANEYMCCB+kIVinEsu/\n9ECMPor8umKcODnTWEHSoJ7Nii/4p36VZOZkJnW6LIPT6cRstWNsttLYZMXQZKGxyUKD0dry00Kj\n0UK90UyDwUKdwUyD0UyzueNBZw4n1DSYqWkwc7Ko7rz7KBUyIkM0RIZoiQoNICpUS3SYluhQLVGh\nWhRyURM6n9jAcxp/G8+IJDPA9Ksk0xUSiQS1Uo5aKSc0SNPl40wWG3WNZuoMZmobTNQ2mqlpMFHb\nYKa6vpmaBhNVdSYamy68mqXFaqeo3EBRuaHdZ1IJRIZoiY3QEed+BZIQFUhgwMV9e3Bu4694wjTw\nDLgk01NqpZyoUDlRoR0vmWq22qmub6aqzvWqqG2moqaJ8pomKmubqahtwu5of7vmcMKZaiNnqo3s\nP9b2ixSiV5EQqScpRk9StJ7kmCDiIwMvmprPuY+xS0Xj74Ajkkw3qRQyYsJ0xITpzvu53eGkuq7Z\nlVCqjJRVGymtMlJaaaC0yojV5mh3jOs2rJJDJyvdZXKZhIQoPUPighkSF8TQhEEkReuRD8DG5+jA\nSCRIcOIUHfIGIJFkPEwmlRAREkBESABjh4a3+czucFJZ20RxhYHiCgNF5Y0UljVQWN5Ik6ntKGGb\n3cnpknpOl9Tz9R5XmVIuJSUumOGJgxiZHEJqUijBgaq++tW8RilTEKENpdxYRamhAofDgVQ68JLp\nxUokmT4kk0qIamkEnpDaeovgdDqprGsm/0wDeaX15JU2cLqknjNVbedBsdgcHMuv4Vh+DZ99nwtA\nbLiOMUPCGJMSStqQ8H6bdGL0UZQbq7DarVQ11RChC/N1SAOS0+kEhwOH1YrDasVps7l/urZbftqs\nOK026s+c6fU1RZLxAxKJhIhBAUQMCmDSyNaxPMZmK7kldZwqqiOnsI4ThbVU1bWdC7ek0kBJpYGv\nsvIBSI7RkzEsgvGpEaQmhfabdp1YfRQHz/wEQHFD2UWVZBw2Gw6TCXuzCbvJ9XKYTNjNZhwtL9e2\nxfXeYsFutuCwtL6cVqtr22rFYbHisFpwWKw4rRYc1pZE0pJYcHa9i0eV5cIPOrpKJBk/ptUoSBsS\nTtqQ1tuu6vpmjufX8nNeNUfzqskrqefcdua80gbyShvYtO0UGpWccSMimDwqiokjo9Bp/HcZjTaP\nsRvKGBcz2ofRdJ/DasVa34C1oR5rfQO2xkasDY3YDAZsBgN2oxGb0YjN2IS9qQl7UzP25iZsTc04\nrVZfh+9VIsn0M6FBGqaO1TB1rGseFmOzlaN51Rw5VcWhnEryz7TOydJstvHD4VJ+OFyKXCZh7NBw\npqfHkjkmmgC1fyWcNo+x/WwMk9Nux1xZiam8AlN5BebKSixVVZira7DU1GCtrcNmaN9twaekUqQK\nBVKlAolc4dpWyJEqlS3v5UjkcqQKBRKFAqlcjkShQCKTucrkcqQKOXJDI7yW06tQRJLp57QaBZNG\nRrlvs2oaTGQfr2D/8XIOnqhwNyjb7E4OHK/gwPEKXt94mMmjopg1MYGM4RHIpL7vpdx2NLZvkozT\n6cRcXo4h9zTG/AKaCotoLi7GVFaO0+bB6TslEmRqNbIADTJN60uqUrm21SpkajVStRqZSoVUpUSq\natlWKl3vlee8VEqkirPvXQlFIvPMmDx5cTG89mrvzuGRSAS/EaJXc8WkBK6YlIDV5uDo6Sqyjpwh\n68gZahtd47qsNgc7D5ey83ApoUFqrpqcyOxLErvVedHTAlU69CodDWYDJQ1lOJ1Orw/RcNrtGE7l\nUn/kJxp+PkbjiZxu10gkMhmK4GCUg4JRBAejCApCoQ9ErtejCNQh1wUiD9Qh12qR67TIAgKQaTRI\nLqKnZyLJDGAKuZT0YRGkD4vg3nlpHM2rZvvBEnYeKsHQ7GoHqK438eHXJ/jnlhymjY3l+ssGMzR+\nkE/ijdVH01B5EoPFSKPZgF4d6PFr2JubqT2QTdWu3dQdOozd2PFKBhKZDHVUJOqYaNRRUagjI1CF\nh6MKC0MZFopCr7+oEkZPiCRzkZBKJYxJCWNMShj3zB3N3qPlfL23gIMnKnA6XX14vj9YzPcHixk7\nNIybrhzOmJS+fcITq4/iWOVJwPWEaaSHkozT6aTx2HHKvt5C9Q+7cFzgiYlMqyVw6BC0KYPRDU4m\nICEBdUw0Uj+YJ7c/8/h/PYvFwvXXX096ejrPPfecp08veIBCLmPq2Bimjo2hrNrIV1n5fLW7AGNL\n7ebwySoOn6wibUgYt/9qJMMS+qZm88snTCMjhvbqfE6Hg+rdeyjeuAlj7ul2n8s0GoLSRhOcnk7Q\n6JFo4uJErcQLPJ5kVq5cSV5eHunp6Z4+teAFUaFaFl87ipuuHM6WvYV8+v0pKmtdfXF+PFXFI69u\nZ9rYGO68bjThg7zbZuPJJ0x1Px4hb91bNBUUtimXaTSETskkbGomQWljkCr86ynbQOTRJPPzzz/z\n3nvv+XzFOqH7NCo5100fzNWZSXx3oIh/fHPCnWx2Hi5l37FybrpiGPNmDPHa+Kk4D0wqbqmr5/Ta\nN6n+IatNeUBiAjG/vpawqVOQaXzXwH0x8liSsdvtPP744yxZsoRvvvnGU6cV+phCLmX25ERmjIvj\ny135fLTlBI1NVswWOxu+PMbOw6U8uDCD5Jggj187JCAYlUyJ2W7pUZKpO3SYnFdew1rbOt+PNjmJ\nhFsWMmgpc7ZrAAAgAElEQVTiBDGhmI947J+ktWvXYrPZ/GLFOqH3lAoZcy9LYfWyK7jqkkTOfj9P\nl9Tz8Cvf8+m2U+3WZ+4tqURKTEt/maqmGky2rk2l6nQ6Kd64iaNPPu1OMHKdjsFL72Hsi88TMmmi\nSDA+5JEkk5uby5o1a3jmmWf8YsU6wXP0WiX33ZjO8/dNJzbcNdeOze7krS+O8szbe92Pwj2lzWJv\nXZjAyulwkLf+bQre+8A9Jkc/ehTpr71E9JyrPNYpTei5XmcEp9PJE088wY033khaWlqvA5o/f367\nMosHBmkJvTMiKYRXH5nJhn//zOc7XE9q9hwt45FXvud/777kgvPrdFfsL9plBockXHBfp9PJqZWr\nqNj6rbssbsENJCxcIJKLhy1duhSlsu0Mjps2berSsb2uyWzYsIGysjIeeOAB7HY7tnO6X9vtHc+b\nK/QvKoWMu+eO4Y+LJxKgdv37VFpl5NHXdnA8v8Yj12iTZBo7nmag6O//bE0wEgmDl95N4qKbRYLx\nM72uyWzZsoWysrI2S2RKJBKOHz/OZ599xtatW4mJieny+c6XHYuLi5k1a1ZvQxU8JHNMDInRep5a\nt5uSSiMNRguPr97Fn++a3G6iru6KC2p9jF1Yf+EkU/HtdxT982P3+6EP3k/EjMt6dW3hwlavXt3p\nJP4X0usk8/TTT2P8RdfsRx55hOTkZO6//34iIiJ6ewnBD8WE6Xj+/kt55u09/JxXg8Vq56n1e3hy\nySWMGdLznsLRuggUUjlWh42iupLz7mM4ncep11e73yfetkgkGD/W69ulpKQkRo0a1ealVqsJDg5m\n5MiRoiF4ANNrlTx97xTGjXD9Q2Kx2lm+fnevbp1kUpn7lqncWIXJamrzucNq5eSrK9yjoiOumEXs\nb+b1+HqC93mlV5VEIhGPDC8SSoWMxxdPImOY6zbJbLHzzNt7Ka/peLG8jiQEtS4cX9TQ9pap6KON\nNOUXAKCJiyPl3iXib83PeSXJfPrppzz77LPeOLXgh5QKGY/fOZmRya6e3nUGM0+v302TqWePtxOC\nW5NM4Tm3TIbc0xRvbGmzk0oZ+uD9SJUX95pV/YEYDSZ4hEoh44+LJxEZEgBAQVkjr/zjYI867J1b\nkymsL3Vv57+zARyuJWXi5s8lcOiQXkYt9AWRZASPCdKp+PNdk92Pt7OOnGHL3sJOjmovIbj1aWRh\nvasmU3foMPU/HgFAFR5G/MIFHohY6AsiyQgelRCl574bWkfgr/3sCKVV3ZttbpA6CJ3S1bu4sK4E\nh8NBwfsfuj+Pv/kmMXq6HxFJRvC46RmxzBzv6lNhsth56cNsHOdZuvdCJBIJCUGu2kyjxUjRjm0Y\nTp4CXI294nF1/yKSjOAV985LI6Jl/pkTBbVs3de92yZ3u4zTSenGT93libeKHr39jUgygldoNQp+\nd85t04Yvj7ln3uuKs+0yUVU2HIWuxt+AxARCLpns2UAFrxNJRvCacSMi3Eu11BnM/OObE10+9mxN\nZmxOa3+b6GuvEX1i+iGRZASvuuv6Ue6Z9L7YcZrSyq41AscHxRDQbGdokWtOGZlWS/il070Wp+A9\nIskIXhUTpmPuZSmAa0WEj7ee7NJxGoWazCIZMle3GMIuvwyZWu2tMAUvEklG8Lr5M4egUbn6znx3\noIiKLgw5cNrtDDvhWnLXCVgzx3gzRMGLRJIRvC4wQMk1U5IAV23mk+86r83U/3QUZaNrcGRBtJI8\nhZ+tNS10mUgyQp+Ye9kQlArXo+dv9hZS02DqcP/K7Tvc28eT1JyqzvdmeIIXiSQj9IngQBVXXZII\nuNbi/tfO9outneWwWKjO2u3aVybhdJySUzX5fRGm4AUiyQh9Zt5lQ5C2PIH+Zm8hNrvjvPvVHsjG\nbnS125QPHoRVIaXMUInB3PG61YJ/EklG6DPhgzRMSG3pN9NoZs/R86+tVPl9662Sc9wI9/apmgLv\nBih4hUgyQp+6OjPRvf3frPx2n9uMRmr2HwBAHqgjamJrD19xy9Q/iSQj9KlxIyIJC3aNaTqYU0lZ\nddtboJp9B3BaXcMPQqdkMiQixf2ZSDL9k0gyQp+SSSXMnnxObWZ321ugmt273dthUzKJCYxEI3d1\nwsutzvf4qpWC94kkI/S5KycluBuAtx0ock8DYTebqc0+BLiWmdWPHoVUKnUv8FZvbqTMUOmTmIWe\nE0lG6HNhwRrSWtZnqqo3caxldYO6Q4dxmF1jlQZNnIC0ZaWLURHD3MceKT/ex9EKvSWSjOATl2W0\nzuO745Bris3qrD3ustBLJrm3x0S2PmESSab/EUlG8IlLxsQgl7numX44XIrVbKF2334ApCoVwRmt\nc9GkhCShlqsAOFqRg8N5/v41gn8SSUbwCZ1GwfgRkYBrrpkft+zGZnCNTxo0Lh2ZSuXeVy6VMbLl\nlslgMZJfW9z3AQs9JpKM4DPT089Z+uS7ne7t881+NyZiuHtb3DL1LyLJCD4zaVSUa9Ck00lA3jEA\nJDIZIRPGt9v33HaZnypEkulPRJIRfEajkjNxZCTR5mp0NtdYpaAxo5HrdO32jQ+KIUgVCMCxylNY\n7T1bnVLoeyLJCD6VOTqaYYbWlQxCznmqdC6JRMLoSNctk8Vu5URVbp/EJ/SeSDKCT01IjWRYU2uS\nGTRx4gX3TY8a5d7eXXzQq3EJniOSjOBTkqpyQi2uaTZLVGHkGy+8GsGE2DRkUtfEV3uKDuJwiEfZ\n/YFIMoJP1ezZ697O0SWw+8iZC+6rVQYwNmok4Bpi8HNl1yYlF3xLJBnBp6p2tQ6IzNEmkHXkTIeD\nIKfEtz55yio64NXYBM8QSUbwmebSUoy5rgbcel0YtUo9pVVGCssbL3jMhJg05FLXmKY9xQexO+x9\nEqvQcx5LMg6Hg7fffptrrrmGjIwMfvWrX/HBBx946vTCAFS1c5d7W5rWWkPp6JYpQKkhveWWqcFs\nELdM/YDHkszrr7/OK6+8wty5c1m1ahVz5szh2WefZf369Z66hDDAVO1o7eWb+usr3NtZP104yQBk\nnnPL9H3e7g72FPyBR5KMw+HgnXfeYcmSJdxzzz1ccskl3Hfffdx000289dZbnriEMMAYCwppKiwC\nQDd0KImpg0mK1gOQW1zf4QJwE2PT0ChcE1ntKjpAnanB+wELPeaRJGMwGJg3bx5XXnllm/Lk5GRq\namowmTpeY0e4+FSds65S2PSpAFwyOtpdtruD2oxaoWZm8hQAbA4bW3J3XnBfwfc8kmT0ej1PPPEE\nI0aMaFP+7bffEhUVhVqsYSycw+l0UrXzB9cbiYSwaa6EkTmmNcns6qBdBuDqIZchwdWn5ptT27GJ\nBmC/5bWnSx9//DFZWVksWbLEW5cQ+qn6Iz9hKisHQD8yFVVoKADJMXoiQgIAOJZXTb3BfMFzRAVG\nkBHt6gFca6pnT3G2l6MWekrujZN+/vnnPPnkk8yZM4dFixZ169j58+e3K7NYLJ4KTfAD5V9/496O\nvHKWe1sikZA5OprN23NxOGH3T2XuVSfPZ86wmWSf+QmAL45vITN+PFKJ6JXhDUuXLkWpVLYp27Rp\nU5eO9fj/kbfffpvHHnuMyy+/nBdeeMHTpxf6OWt9vXuaTblOR+iUzDafn3vLtPNwSYfnSotMJV7v\n2v90bSG7i0Rtxh95tCbz0ksvsXbtWubNm8czzzyDVNr9HHa+7FhcXMysWbPOs7fQ31R8uw2nzQZA\n+MzL2syAB5CaFEJYkJqqehM/nqyktsHEIP352/QkEgkL067nhZ2rAfj7j5uZFJuOXOaVCvpFbfXq\n1cTFxfXoWI/VZN59913Wrl3L4sWLee6553qUYISBzel0UnbOrVLU7Cvb7SOVSrg0w/XH7HC2TjJ+\nIRNi0hge5loArtxYxTe5OzrcX+h7HskElZWVvPjiiwwfPpw5c+Zw+PDhNi8xWlYAqDt4CFOp66mR\nfmQqAQnx593vsnGt/2J+f7Dj+XwlEgm3jp3nfr/x6L9pMF14WILQ9zxSr9y5cydWq5WcnBwWLlzY\n7vOsrCyCg4M9cSmhHyv++BP3dtScqy64X3KMnvhIHUXlBnIK6yitMhAT1n62vLOGh6UwKS6dvcWH\naLQYefPA33l4yt1IJBeeNkLoOx5JMvPmzWPevHmd7yhctOqPHqXhZ9c8vuroKMKmTrngvhKJhMvG\nxfH+f1xz+W4/WMLCK4dfcH+AOzNu4mj5CYzWZvYUH+SHwn1MSzz/LHtC3xINJ0KfKP6otRYT95t5\nSGSyDve/LKP1lmnrvkL3UrYXEhIQzB3jbnK/X3/gH1Qaq3sYreBJIskIXteYc5K6Q4cBUIaFET7j\nsk6PiQrVMmqwq5NeWXUT2ScqOj1meuIkJsW6FoUzWpv5245VNFvFkBZfE0lG8Cqn00n+2++638fN\nvx6pQtGlY6+dluze/vcPeZ3uL5FIuGfCLYRrXcmpsL6EV3e/JR48+JhIMoJXVe34wd0Wo4qIIOKK\nrvd3umR0NCEtfWQOHC/nTJWx02P06kCWTf8tGrnruOzSI6ze/75IND4kkozgNXaTifx3NrjfJ9+5\nuF3nu47IZVKuzkwCwOmEL3d1XpsB1xpND01Z4n66tC0vi5V73hGz6PmISDKC1xR//AmWalfja9DY\ntAuuqdSRqy9JRC5rGW29txBjc9cWdUuPHsX9kxe7xzLtLNzHCztX02Rp7nYMQu+IJCN4RcOx4xRv\n+sz1Ripl8JI7e9RvZZBezdQ015rZxmYrn33f9UXdpiVO4sHMu5C1JJrsMz/xhy1/pbih42kkBM8S\nSUbwOFtTEzkvvwot7SBxv5l3wd69XbFw9jCkUleC2rz9VIdTQPzSJfHjeGz6bwlQaAA401jBsq+f\n48ucb3E4RTtNXxBJRvAop9PJ6bXrMZe7Hjnrhg4hfuGCXp0zLiKQWRNcSarZbOfjrd2bPDw9ehTP\nXbmMWH0U4Frm9p2DH/Pkty+RV1vUq9iEzokkI3hU6ef/ovK7bQBIVSqGPfQAUnnvO5bfPHsECrnr\nz/XLXXmUVXf+pOlc0YERPHvFY1wxeJq77HhVLsu+fo439mygwlDV6xiF8xNJRvCYmn372/SJSVl6\nN5rYGI+cO3yQhmumuPrNWG0OVnx0qNNewL+kUai5Z+Ii/njp/YQFhADgxMm2/Cz+58v/ZeWed8iv\n7XhAptB9IskIHtF4IocT//ey61kzEDvveiIun+nRayy8cpi738yPp6r4and+j86THj2SV+b8LzeP\nuR613PVI3eF0sD1/D7//+hme2PIC357eJZ5EeYhIMkKvNeac5OiTT+NoWZUiZPJEEm/r3rSrXaEL\nUHLfjWPd79/+4mi3b5vOUsqVzBt5NSt/9TTzUq92L7ECkFN9mtX73uPuzb/n/3auYVteFg1mQ6/j\nv1iJKcSEXmk4foKfn/oL9ibXOkn6kakMe+iBTgdA9tTEkVHMmhjP1n1FmCx2/vLWHp6/fzoB6q4N\nVfglvTqQm9Ou5/oRs9mWn8XXp7ZT2uia5NzqsLG35BB7Sw4hQcLgkATSIlMZGTGUYaGD2yQm4cJE\nkhF6rOqHLHJefhWn1dVBLjB1BKl/ehyZRuPV6y65fgxHcqupqGmioKyRv723nz/fORmZrOcV8wCl\nhmuGXc6coTM5XnWKHQX72F2UjcHiqik5cZJbU0BuTQGfHvsKqURKQlAMKSFJpIQkkBQcT0JwLEpZ\nz5LdQCaSjNBtToeD4o2bKPzg7+4y/aiRpD7xR+QB3k0wADqNgv+9azK/X7EDo8lG9vEKXvvoEP+z\nIL1XiQZcgyxTw4eSGj6UO8fdxPHKU+wvOcyhsp/dNRxwteHk1xWTX1fM1tOtx0bpwonXxxCjjyQm\nMJLowAiidOHoVYEX7SRaIskI3WJtaOTkK69Se+Cguyzs0mkM/Z/7ujy62hMSovQsu30iT765G7vD\nybf7i2g223j01vEo5J65VZNLZYyOHM7oSNeEWVXGGn6qOMHxylMcr8ptk3TA1UfoTGMFZxor4BdT\nE6vlKsK1oa5XQAihAYMI1QwiNCCYQZpgBqn1qAfo7ZdIMkKX1ew/QO7rq7HU1LjL4hbcQMLNNyHx\nwcTx6cMieOjmcbz092wcDidZR87wpzVZ/P62Ce6nUJ4Upg1hRnImM5Jdy7gYLU3k1hSQX1dEXm0R\nhfWllDaWn3cgpslmpqi+lKL60gueXyNXE6zWE6QOJEitJ0gViF6tQ68KRK/SEajSEajUolNpCVTq\nUMoU/aJ2JJKM0ClrQwP577xHxdZv3WVynY6hD/0PIRPG+zAy16TjGrWcv727D4vNwdHT1fzPi9/x\n0M3jGD8i0qvX1ioDSItKJS0q1V1mc9gpM1RQ2lBOaWM5ZYZKyg2VVBiqqGquxem8cN+eZpuJZoOJ\nM4bOJ+gCUEjl6JRatMoA90unCCBAqUGrCECr1BCgaH1pFOrWn3I1KrmqT5KUSDLCBTntdsq+3kLh\nBx9ia2x9hKsfPYqhD9yHOiLCh9G1mjQyiqfuncJfN+yjrtFMvcHCk2/u5tKMWO68bhShQd5vJzpL\nLpURp48mTh/d7jObw05Ncx3VTTVUGWtbtmupNdVT21xPnameWlMDVnvXRppbHTbXsab6HsUqkUhQ\ny1UEyDWoFSo0cjUahQq1XI1arkItV2Gu7X1fIZFkhHacTic1e/dT8N77NBe19oCVqlQk3X4bUXOu\n8sntUUdGDQ7ltUdm8NKH2RzKqQRcE5Dv+7mMa6cN5vpLUwjSdX0uG2+QS2VEaEOJ0IZC+Pn3cTqd\nNNtMNJgN1JsaaDAbaDA10mgxun82WowYzEaMliYaLQaMlibsPRjs6XQ6abaaXFOUXiCXWESSETzJ\n6XRSu28/RR9txHDyVJvPQqdmknzHYlThYT6KrnODAtUsvzuTL3fl8f5/jmE02dwDKj/fcZrLJ8Qz\nJzOJ5JggX4d6QRKJxH17E6W7QCb6BafTidlmxmhtxmhposna3Gbb9TLRZG3G1PKz2WZyJ5iz21aH\nzSu/k0gyAg6rlcrvd1D6+Rc0FRS2+UybkkLS4tsIThvjo+i6RyqVcO20wUwbG8t7/znGlpaVDswW\nO//Zlc9/duWTEhfE9LGxTB0bQ1So1tch95pEIkGtUKNWqAkNGNTj89gcdkwtScdkM9NsNVFYVMjv\nyOpVfCLJXMRM5eWUf72F8m+2Yq1ve1+viYsjfuECwqZm+t2tUVcEB6q4f0E6N84aysZvT/Ld/iIs\nNtctRW5xPbnF9bzz75+JDdcybkQkY1JCGZkc6vNbKl+SS2XoVK6nV2cFmJS9P2+vzyD0K3aTierd\ne6jY+h31Px5p97k2JYW438wjNHNyv0wuvxQVquW+G9O5/Vcj+XZ/EVv2FpJ/psH9eUmlkZLK03yx\nw9WjLjpUy9D4YAbHBpEUoycxSk9okLpfPCr2VyLJXATsZjN1Bw9R9cMuavbsw2H+xcxyUikhEycQ\nc/216EeOHJBfqMAAJddfmsL1l6ZQVN7IzsOl7D9WxsmiOs59qnym2siZaiPbD7X2ptOoZMSG64gK\n1RIdpiV8UADhwRrCgzWEBqnRavpHfxVfEUlmgLI2NFJ74AA1e/dTm33QPUL6XKqIcCIun0nklVeg\nCgv1QZS+ER8ZyM2zh3Pz7OHUG8z8dLqao6erOZZXTf6ZBmz2tn1Zms12ThXXc6r4/I+KlQoZgwJV\nDApUEaRTERyoQq9VotcqCQxwvbQaBTqNAm3LS62UXTSJSSSZAcJpt2M4lUvtwUPUZR+k8eQp9xy7\n55IFBBB6yWTCZ1xK0JjRA+KWqDeCdCqmpsUwNc01uZbVZievtIH8M65XUXkjxRUGquou/CjXYrVT\nXtNEeU1Tl68rlYBGrUCrlqNRyVGrXD81KjlqpQy1So5a2bKtlKFSyFApZSgVrm2l4txtqXtboZCh\nUkiRy6R+k8REkumnnHY7xrx86o8epf7IURqO/uyebuGXZNoAQiZOJDRzMoPGZSBV9r4xb6BSyGUM\nSxjEsIS2T2lMFpsrkVQ3UVHbREVtM9X1zVTXm6hpMFHXaKLZ3PV1nRxO1+oLXV3ipScUcilKuRSF\nXIZCIUUhcyUjudy1rZC3vuTu9zLkMol7H0N9Za/jEEmmn7A2NGI4eZKG4ydoPH6CxpyT570FOksd\nE82gcRmETJ6EfmSqR+bZvZiplXISo1wNwRdiMtuoN1qoN5hpMFpoMJppbLLS2GTB2GTF0Ox6NZms\nNJlsNJltNDVbaTbbsHdzKtGusNocWG0OoOf9X6xNNZ3v1Anxl+eHrI2NGE/nYcg9jTH3NIZTuZjK\nyjo8RqbVEpw2mqC0NIIz0tFER/VRtMJZ6pbbnsiQgG4d53Q6sdocNJttmCz2lp82zGY7JourzGK1\nY7LYMVtd2+aWMrPVjtXmcP+0tHxutbu2z5adTThWm8MrCa0jHk0yH330EevXr6esrIzU1FSWLVtG\nenq6Jy8xoDisVppLSmkqLKSpoBBjQQFN+QWYKzufOV8RFIR+5AgCU1MJGjMabWKC12ajE7xLIpG4\n21j6oi+y3eHEarNjO5t47A7Xtt313mZvTUhnzpTw8Ledn7MjHksyn376KU8++ST33Xcfo0eP5v33\n32fJkiVs3ryZ2NhYT12m33E6ndgaGmguPUNzaSnNxSWuxFJU7KqddGEheIlcjjYpEd3QIQQOH0bg\niOGoo6L8pmFP6F9kUgkypRy60DRXrLX0+noeSzIrVqxg4cKF/Pa3vwVgypQpXH311bzzzjs8/vjj\nnrqMX3JYrZirqjGXl2OqqMBcXkHzmTJMZeWYys5gN3bjqYNSSUBiItrByegGJ6NNGYw2KbFPJ4QS\nBE/ySJIpKCigtLSUmTNbl8CQy+XMmDGDHTt2eOISPuOw2bDW1mGpqcFcXY2lutqVUCqrsFRVY66s\nxFJbCx3ME3I+ErkcTUw0mvg4AhISCEiIR5uUiDoyUtz2CAOKR5JMfn4+EomExMTENuVxcXEUFRXh\ndDr9qmrvsFqxNjRia2jAWl+Ptb4BS10d1ro6rHX1ru3aWiy1da4xPd1MIOdShoagjo5GEx2NOiYa\nTWwsAXGxqCIjxBMf4aLgkb9yg8E1oZFW23ZEq1arxeFw0NTU1O6z3nJYrdibmrA1Nbt+Go3Yja6f\nNqMBm8GIzdDys7ERW2Mj1sZGbI2GC/Yn6Qm5TocqPBxVRBiqiEjUkeGoIiLRREeiioxEprp4B9wJ\nAngoyZydUvBCtRWph3qVvvXUZ+jVepwORxdrF3IgqOUFKICQlldnJBIkMlnrSy5HKpcjkZ+7LW/b\nY7ax5XWqGchveQk94Z6msu2PczZcy5ScLXO6DvrlLq3H/+LvpaM/n46myGy/b4ef9vL4vvfL373B\nUN3rc3okyQQGBgJgNBoJCWn9BhuNRmQyGZpurMMzf/78dmXmlgF9FRY7BpkTkLS8vMzuAByA93pl\nCoI/a2quA+Dee+9F9Yta+aZNm7p0Do8kmcTERJxOJ0VFRcTHx7vLi4uLSUpK6vX5CwoKAPgm641e\nn0sQhO6z2WztkkxXeSTJJCUlER0dzZYtW5gyZQoAVquVbdu2tXni1BXny45z586lubmZtWvXIusn\nT16WLl0KwOrVq30cSdf0t3hBxNwX7r33Xmw2G5s3b0at7tkyMx57vHH33Xfzl7/8hcDAQMaNG8f7\n779PXV0dt99+e6/PLZVK0Wq17Z5e+TNlyyDEuLg4H0fSNf0tXhAx9wWVSoVKpepxggEPJplbbrkF\ni8XChg0b2LBhAyNGjOCtt97qN/8xBUHwDo921Fi8eDGLFy/25CkFQejnLu4ZiwRB8DqJszudAgRB\nELpJ1GQEQfAqkWQEQfAqkWQEQfAqkWQEQfAqkWQEQfAqkWQEQfAqkWQEQfAqkWQEQfAqkWQEQfAq\nkWQEQfAqv08yH330EVdddRVjx45l4cKFHDp0yNchdcjhcPD2229zzTXXkJGRwa9+9Ss++OADX4fV\nJRaLhTlz5vCHP/zB16F0KisriwULFjB27Fguv/xyVqxY0a1pM/uaw+HgzTffZPbs2WRkZLBgwQJ2\n797t67DOa+vWrYwbN65d+apVq5g5cybp6enceeednD59ukvn8+skc3bBuOuvv54VK1ag1+tZsmQJ\nJSUlvg7tgl5//XVeeeUV5s6dy6pVq5gzZw7PPvss69ev93VonVq5ciV5eXm+DqNTBw4c4O6772bI\nkCGsXbuWW2+9lTfffJM33vDfmRPXrVvHK6+8wg033MAbb7xBfHw8S5Ys4fjx474OrY3s7Gx+//vf\ntytfuXIla9asYcmSJbz88ss0NjZyxx13uBcR6JDTj82cOdO5fPly93ur1eqcNWuW8y9/+YsPo7ow\nu93uHDdunPO1115rU758+XLnlClTfBRV1xw9etSZnp7uzMzMdC5btszX4XTolltucS5durRN2Ysv\nvui87bbbfBRR5+bMmdPmv6vdbnfOmDHD+fTTT/swqlZms9m5du1a5+jRo52TJk1yZmRkuD8zGAzO\njIwM57p169xl9fX1znHjxjnffvvtTs/ttzWZ/rhgnMFgYN68eVx55ZVtypOTk6mpqcFkMvkoso7Z\n7XYef/xxlixZQkREhK/D6VBNTQ3Z2dncdNNNbcoffvhhNmzY4KOoOmexWNosCySVStHpdNTV1fkw\nqlbbt29n3bp1LFu2jFtvvbXNZ4cPH6a5ubnNd1Gv1zNx4sQufRf9Nsl0ZcE4f6PX63niiScYMWJE\nm/Jvv/2WqKioXk1h6E1r167FZrNx7733+jqUTuXk5ACuaSGXLl1KWloaU6ZMYeXKlX75N3HWokWL\n2Lx5M1lZWRgMBt59911yc3O59tprfR0aAGlpaWzdupVFixa1W9ro7C10QkJCm/L4+Hjy8/M7Pbff\nLmHoiwXjvOHjjz8mKyuLP/3pT74O5bxyc3NZs2YNGzZsQN4PVrSsqanB6XSybNkyrr32Wu688072\n7kwuWwsAAAOBSURBVN3LqlWrUKvVLFmyxNchntfNN9/M7t27ueOOOwDXGmUPPvggM2bM8G1gLTqq\nwRqNRpRKZbu/D61W26U2Gb/9q3L20YJx3vT555/z5JNPMmfOHBYtWuTrcNpxOp088cQT3HjjjaSl\npfk6nC6x2WwATJ8+nUcffRSASZMmUVtby6pVq7jrrrv8aknks84+jVm+fDmDBw9m165drFixgsDA\nQG655RZfh9chZwfLTHfle+i3ScaTC8b5wttvv83zzz/PFVdcwQsvvODrcM5rw4YNlJWV8eabb2K3\n29vcbtjtdr9cfiYgIACAadOmtSmfMmUKH374IcXFxW3W/vIHBw4cIDs7m9dee43Zs2cDMHHiRGw2\nGy+88ALz5s3z679nnU6HxWJp9zdhNBrd39OO+G114NwF487lqQXjvOmll17ib3/7G3PnzuXVV1/1\n29uQLVu2UFZWxoQJExg1ahSjR4/mxIkTfPrpp4wePZrS0lJfh9jO2TY6q7Xtqp5nazj+WIspKytD\nIpEwduzYNuXjx4/HZDL5dZcMcK2r5nQ6KS4ublNeVFREcnJyp8f7bZI5d8G4s84uGJeZmenDyDr2\n7rvvsnbtWhYvXsxzzz3n17d1Tz/9NBs3buSTTz5xvxITE5k5cyaffPKJXz5pGjJkCJGRkXz11Vdt\nyrdt20ZERIRfLsFz9kt64MCBNuWHDh1CJpMRFRXlo8i6JiMjA6VS2ea7WF9fz759+7r0XfTPf2Jb\neHPBOG+orKzkxRdfZPjw4cyZM4fDhw+3+XzMmDF+lXTOVyNUq9UEBwczcuTIvg+oCyQSCQ899BB/\n+MMfePLJJ7nqqqvYtWsXmzdvZvny5b4O77xGjRrFjBkzWL58OXV1daSkpLBnzx7WrVvH7bffjk6n\n83WIHQoICODWW2/l1VdfdT/xXb16NXq9nhtuuKHT4/06yfS3BeN27tyJ1WolJyeHhQsXtvs8KyuL\n4OBgH0TWdRKJxC9vOc41d+5clEolq1ev5tNPPyUqKorly5dz4403+jq0C3rttdd4+eWXWbNmDfX1\n9SQmJvLnP/+ZBQsW+Dq08/rl38DDDz+MTCbjrbfeoqmpiXHjxvH88893KUGKJVEEQfAq/6m7C4Iw\nIIkkIwiCV4kkIwiCV4kkIwiCV4kkIwiCV4kkIwiCV4kkIwiCV4kkIwiCV4kkIwiCV/3/qqoXUwE/\nFYIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x115ff6910>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum temperature reached is just under 900 K\n",
      "Exit flow rates Fc=4.59, Fd=0.003ish mol/min\n",
      "Selectivity C/D = about 1200\n"
     ]
    }
   ],
   "source": [
    "initial_values = (Fa0, Fb0, Fc0, Fd0, 1.0, T0, Ta0)\n",
    "volumes = np.linspace(0,10,100)\n",
    "result = scipy.integrate.odeint(dVector_dV, initial_values, volumes)\n",
    "make_plots(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Counter-current heat exchange\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "With cooling water exit temperature 400.0\n",
      "Actual entrance Ta is 187.696323043\n",
      "Target entrance Ta is Ta0 = 325\n",
      "[ about 500ish ]\n",
      "Trying again with Ta_exit = about 500ish\n",
      "Actual entrance Ta is 325.000000816\n",
      "Target entrance Ta is Ta0 = 325\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10afcdd50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.figure.Figure at 0x116d73d10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum temperature reached is about 1100 K!!\n",
      "Exit flow rates Fc=4.59, Fd=0.003ish again mol/min\n",
      "Selectivity C/D = about 1200\n"
     ]
    }
   ],
   "source": [
    "Ta_exit = 400. # guess the cooling fluid temperature at V=0\n",
    "initial_values = (Fa0, Fb0, Fc0, Fd0, 1.0, T0, Ta_exit)\n",
    "volumes = np.linspace(0,10,100)\n",
    "result = scipy.integrate.odeint(dVector_dV, initial_values, volumes, args=('counter-current',))\n",
    "Ta_entrance = result[-1][6] # cooling fluid temperature at V=10 dm3\n",
    "print(\"With cooling water exit temperature {0}\".format(Ta_exit))\n",
    "print(\"Actual entrance Ta is {0}\".format(Ta_entrance))\n",
    "print(\"Target entrance Ta is Ta0 = {0}\".format(Ta0))\n",
    "\n",
    "# We could guess many times until we're right\n",
    "# but instead we can use scipy.optimize.fsolve to find it for us\n",
    "def make_this_zero(Ta_exit):\n",
    "    \"We want to find the Ta_exit that will make this function return zero\"\n",
    "     \n",
    "    return 0####FIXME####\n",
    "\n",
    "fsolve_result = scipy.optimize.fsolve(make_this_zero, 400)\n",
    "print fsolve_result\n",
    "Ta_exit = fsolve_result[0]\n",
    "\n",
    "print(\"Trying again with Ta_exit = {0}\".format(Ta_exit))\n",
    "initial_values = (Fa0, Fb0, Fc0, Fd0, 1.0, T0, Ta_exit)\n",
    "result = scipy.integrate.odeint(dVector_dV, initial_values, volumes, args=('counter-current',))\n",
    "Ta_entrance = result[-1][6]\n",
    "print(\"Actual entrance Ta is {0}\".format(Ta_entrance))\n",
    "print(\"Target entrance Ta is Ta0 = {0}\".format(Ta0))\n",
    "\n",
    "make_plots(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Constant temperature heat exchange \n",
    "This is the case of having infinite cooling flowrate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.figure.Figure at 0x116738210>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.figure.Figure at 0x116ec7610>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum temperature reached is a bit over 800 K\n",
      "Exit flow rates Fc=4.45, Fd=0.002ish mol/min\n",
      "Selectivity C/D = about 1900\n"
     ]
    }
   ],
   "source": [
    "initial_values = (Fa0, Fb0, Fc0, Fd0, 1.0, T0, Ta0)\n",
    "result = scipy.integrate.odeint(dVector_dV, initial_values,\n",
    "                                volumes, args=('constant-Ta',))\n",
    "make_plots(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adiabatic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10a2f0d90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.figure.Figure at 0x116d998d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum temperature reached is over 1500 K!!\n",
      "Exit flow rates Fc=4.77, Fd=0.007ish mol/min\n",
      "Selectivity C/D = under 700\n"
     ]
    }
   ],
   "source": [
    "initial_values = (Fa0, Fb0, Fc0, Fd0, 1.0, T0, Ta0)\n",
    "result = scipy.integrate.odeint(dVector_dV, initial_values,\n",
    "                                volumes, args=('adiabatic',))\n",
    "make_plots(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Now play around with reversible reactions (with equilibrium constant being a function of temperature) and pressure drop..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}