Voltz_comp.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:0fa27a4c5d785c3f19055cf4baec6603cab9df64ae69dda6921e148156394ee2"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import cantera as ca\n",
      "import os\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as pl\n",
      "import scipy.optimize as opt\n",
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 51
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def get_rates_cantera(temperature,co_conc,gas,interface,dt=10000.):\n",
      "    gas.TPX = (temperature,ca.one_atm, [0.1,0.0,0.0,0.9-co_conc,co_conc])\n",
      "    interface.TPX = (temperature,ca.one_atm,[1,0,0,0])\n",
      "    interface.advance_coverages_test(dt)\n",
      "    return interface.net_production_rates[1]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 100
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "temperatures = np.linspace(600,1000,10)\n",
      "co_conc = 0.0003\n",
      "ref_folder = \"/home/vonrickenbach/testcases_mb/kelvin_cell_final/wc_1d_tests/ref_sim/finite_vol_avg_wc_only_CO/\"\n",
      "\n",
      "gas = ca.import_phases(os.path.join(ref_folder,\"cantera_input.cti\"),[\"gas\"])[0]\n",
      "interface = ca.Interface(os.path.join(ref_folder,\"cantera_input.cti\"),\"pt_surf\",[gas])\n",
      "r_const = [9.2E12,2.7,12028.0,-962.0]\n",
      "\n",
      "rates_cantera = [get_rates_cantera(temperature,co_conc,gas,interface) for temperature in temperatures]\n",
      "rates_voltz   = [get_rates_voltz(r_const,temperature,co_conc,gas,interface) for temperature in temperatures]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 115
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def get_rates_voltz(r_const,temperature,co_conc,gas,interface):\n",
      "    gas.TPX = (temperature,ca.one_atm, [0.1,0.0,0.0,0.9-co_conc,co_conc])\n",
      "    interface.TPX = (temperature,ca.one_atm,[1,0,0,0])\n",
      "    \n",
      "    k1 = r_const[0]\n",
      "    k2 = r_const[1]\n",
      "    E1 = r_const[2]\n",
      "    E2 = r_const[3]\n",
      "    \n",
      "    co_conc = gas.concentrations[4]\n",
      "    wc_thickness = 1E-5 \n",
      "    return k1*np.exp(-E1/temperature)*co_conc / (1. +k2 *np.exp(-E2/temperature) )**2 * wc_thickness"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 116
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def min_fun(r_const,temperatures,co_conc,gas,interface):\n",
      "    pen_val = 0.0\n",
      "    for idx,temperature in enumerate(temperatures):\n",
      "        pen_val = pen_val + (rates_cantera[idx]\n",
      "                            -get_rates_voltz(r_const,temperature,co_conc,gas,interface))**2\n",
      "    return pen_val"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 117
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "opt_result = opt.fmin(min_fun,r_const,args=(temperatures,co_conc,gas,interface),maxiter=100000)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Warning: Maximum number of function evaluations has been exceeded.\n"
       ]
      }
     ],
     "prompt_number": 118
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "rates_cantera = [get_rates_cantera(temperature,co_conc,gas,interface) for temperature in temperatures]\n",
      "rates_voltz   = [get_rates_voltz(opt_result,temperature,co_conc,gas,interface) for temperature in temperatures]\n",
      "\n",
      "pl.plot(temperatures,rates_cantera)\n",
      "pl.plot(temperatures,rates_voltz)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 119,
       "text": [
        "[<matplotlib.lines.Line2D at 0x44e6250>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAD/CAYAAAAuYdt6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVPX3P/AXq8LMMILIkruisgnmkpCmoKamluCSgKWm\nllqpaSqYG6aWtogtv4L8umWCmqK5fdCAMCnXTBTGVFBAMVkcRGZkG+b8/rgxNrIaywxwno8HDz/3\nvu99z7n3Q/Pi3vu+9xoQEYExxhhrQIa6LoAxxljzw+HDGGOswXH4MMYYa3AcPowxxhochw9jjLEG\nx+HDGGOswVUbPnK5HL6+vhCLxejUqRMiIiIqXTYkJAT29vaQSqWYMWMGiouLa9xPTEwMHB0dIRKJ\nMGTIEKSnp2vafvnlF3h7e6NVq1bo3Llzuc9NTU2Ft7c3RCIRnJycEBMTU6ONZ4wxphvVhs8777yD\nli1bIisrC7t27cKcOXMgk8nKLXf8+HFs2LABsbGxSEtLw82bN7Fq1aoa9ZOTk4Px48dj3bp1yM3N\nRd++fTFp0iTNumKxGDNnzsSnn35aYY3+/v7o06cP5HI51q1bhwkTJiAnJ+epdwZjjLEGQlVQKBRk\nampKN27c0MybMmUKBQUFlVvW39+fli1bppmOjY0lOzu7GvUTFhZGAwYM0LQplUoyMzOja9euaX3G\nzz//TJ06ddKad+3aNWrRogUpFArNvEGDBlFoaGhVm8YYY0yHjKsKpuvXr8PY2BgODg6aee7u7oiL\niyu3rEwmg6+vr2bazc0NmZmZyM3NRWpqapX9JCUlwd3dXdNmbm4OBwcHJCYmonv37lWGZ1JSErp0\n6QKRSKTVd1JSktZyBgYGVfbDGGOsclTHD8Op8rSbQqGAhYWF1jyJRIL8/PwKl5VKpZrpsvXy8/Or\n7aeidgsLCygUimo34MnPLVu3ohoZY4zphyrDRywW4+HDh1rz8vLyIJFIql02Ly8PgBAylfVTFjgS\niaTGn1OTGh88eFAuzMoQkV7/rFq1Suc1cJ1cZ2Otkeusu5+8wjwM3ja42u/g/6rK8OnevTtUKhWS\nk5M18xISEuDq6lpuWRcXF1y6dElrOVtbW1haWlbaj4uLi2bdhIQETZtSqURKSoqmvSouLi64efOm\n1lHSv/tmjDH2dHIe5WDo90NxMu1kvX1GleEjEokwbtw4rFy5Eo8ePUJ8fDwOHz6M119/vdyyU6ZM\nwZYtW3D16lXk5uZizZo1eOONN2rUj6+vLxITExEZGYnCwkKsXr0avXr10lzvISIUFhaipKQERISi\noiLNMO7u3bujV69eWL16NQoLCxEZGYnExESMHz++TncUY4w1BxkPMzBo2yBcuHsBXSy71NvnVDvU\n+ptvvkFBQQFsbGzw2muvITQ0FE5OTkhPT4dEIsGdO3cAACNGjMCSJUvg7e2NTp06oWvXrli9enW1\n/QCAtbU19u/fj2XLlsHKygoXLlzA7t27NeuePHkS5ubmGD16NG7fvg0zMzOMHDlS0757925cuHAB\nVlZWWLZsGfbv34/WrVvX2U5ijLHmIFmejIHbBuJqzlW42rgi/o34evusKke7AYClpSUOHDhQbn6H\nDh3KXdRfsGABFixY8FT9lBk6dCiuXr1aYZuXlxfUanWl63bs2BG//PJLpe2MMcaqdjnzMobvHI5M\nZSb6t+2PY5OPwcrMqt4+jx+vwxhjzdzp26cxePtgZCozMbTzUERPiYaVmRWysurvMzl8GGOsGfs5\n5WcM2zkMDwofwNfRF0cDjkJsKsaPPwL1OW6Lw0fPeHl56bqEGuE661ZjqLMx1AhwnU8j8mokxkSM\nwaOSR5jqPhV7J+5F/oMW8PMDXn0VqM+nlBkQUd3etqqnyp5w0Ew2lzHGqrT90nbMODQDalJjfv/5\n2DhiI346aIjZs4GsLEAkAj79FHj77fr57uTwYYyxZmbTmU1YcFwYHBY8OBjv9lyJ+fMNsGuX0O7l\nBWzdCnTuXH/fndWOdmOMMdY0EBFWxa3Cml/XAAA2jdiErjnz4eoK3LsHmJsDGzYAb78NGNbzRRk+\n8mGMsWZATWq8F/Uevjr3FQwNDPH1sK04+91U7NghtA8YAGzfDvzr+c8A+MiHMcbYf6RSqzD9p+nY\neXknTI1MsaTzbqx71RcZGUDLlsBHHwHz5gFGRg1XEx/5MMZYE1aoKoTfPj/8dO0niExEeCHjIKK+\nGQYA8PAQjnZ69Kh8/fr67uTwYYyxJiq/KB8+e3wQeysWEmNLmEUeQ9ZFD7RoAaxZAyxcWP3RDodP\nLXH4MMaak/uP7mNU+CicyzgH81I7PAo7AWT1RL9+wtGOs3PN+qmv706+yZQxxpqYu/l3MXj7YJzL\nOAfj/M549HU8THJ7Yt064Pffax489YkHHDDGWBNyM/cmhu4YhtS8W0CWM1Q7T+BZh7bYsQPo2VPX\n1T3GRz6MMdZEJGYl4rnQgULwZPSD0c5fsXpRW5w9q1/BA/CRD2OMNQknk89ixM6XUGSYC9zyhsuV\nn/DDrxL06qXryirGRz6MMdbIfXUkBt7bhwrBc+0VBLU/houn9Td4AD7yYYyxRquwEPBffRAHjScB\nJsWQpr6GqPe2wuM5E12XVi0+8mGMsUbo3Dmgi+/3OGg6ATAuRm/VO/j72x2NIngADh/GGGtUioqA\nDz4APOZ/ib89pgKGpXij63Jc+PArmLVsPF/pfNqNMcYaiT/+AKZMJcisPwRGBgMA1nt/jsBBC3Vb\n2H/QeGKSMcaaqeJiYOVK4Ln+asjaLwC8g2EIQ2x5ZUujDB6Aj3wYY0yvJSQAU6cCCVdUwMtvAs9u\nh4mhCcLHh2OC8wRdl/ef8ZEPY4zpoZIS4eGfffsCCYlFMJ/6KvDsdpibmOOw/+FGHTwAH/kwxpje\nSUwUjnYuXgRgqkC7Rb64YxoNaQspjk0+hufbP6/rEmuNj3wYY0xPqFTAxx8DffoIwdOumxzOH7+I\nO6bRsBHZ4OS0k00ieAA+8mGMMb2QmgpMmiTcvwMAr835GxddhkOWk4gO0g6Ifj0a3Vp302mNdYmP\nfBhjTMeOHxeOds6dA9q3B3b8dAunHV+ALCcRjtaO+G36b00qeAAOH8YY0xm1Gli7FnjpJUAuB0aN\nAvbEyrA0eSBSclPQ2743fp32K9pZtNN1qXWOT7sxxpgOPHgAvP46cOQIYGAABAcDI2ecx+jwl3C/\n4D4GdRyEQ36HIG0p1XWp9YLDhzHGGtjly8C4cUBKCmBpCezaBRh2O45hOydAUazA6G6j8ePEH2Fm\nYqbrUusNn3ZjjLEGtGsX4OEhBE+vXsCFC8Adm80YHT4aimIFAnoG4MCkA006eAAOH8YYaxDFxcC8\necBrrwEFBcJ9PPG/qfF/tz7AW0feQimVYunApdjpuxMmRo3jydS1wafdGGOsnt29C0ycCPz+O2Bi\nAnz1FTBleiGmH3oDuxN3w8jACN+O/hZv9nlT16U2GA4fxhirRydPCvfvZGYC7doB+/YBDj3vY/gP\nPohPj4fYVIx9E/dhhMMIXZfaoKo97SaXy+Hr6wuxWIxOnTohIiKi0mVDQkJgb28PqVSKGTNmoLi4\nuMb9xMTEwNHRESKRCEOGDEF6erpWe2BgIKytrWFtbY2goCCttqSkJHh5eaFVq1Zo37491q5dW6ON\nZ4yx+kIEhIQAQ4cKwePtLbwSwbpbCjy3eCI+PR5tJW0R/0Z8swseAABVw8/Pj/z8/EipVFJ8fDxJ\npVJKSkoqt1xUVBTZ2tqSTCaj3Nxc8vLyoqCgoBr1k52dTVKplPbt20dFRUW0ePFi8vDw0KwbGhpK\nPXr0oIyMDMrIyCBnZ2cKDQ3VtD/77LO0fPlyUqvVlJKSQvb29nTo0CGt+gBQDTaXMcZqLT+f6NVX\niYQIIlqyhKikhOj39N/J+hNrQjDI7Vs3up13W9elVqu+vjur7FGhUJCpqSnduHFDM2/KlClaoVLG\n39+fli1bppmOjY0lOzu7GvUTFhZGAwYM0LQplUoyMzOja9euERGRp6cnbd68WdO+detWrXBq0aIF\nXb16VTM9ceJEWr9+vfaGcvgwxhrAX38ROTsLoSOREO3bJ8zfl7SPWq5tSQgGjdg5gvIK83RbaA3V\n13dnlafdrl+/DmNjYzg4OGjmubu7IykpqdyyMpkM7u7ummk3NzdkZmYiNze32n6SkpK01jU3N4eD\ng4OmvaK+/13D8OHDsWPHDqhUKvz11184ffo0hg0bVu1RH2OM1aXISKBfP0AmA5ychMfljBtH+Pz3\nzzHxx4koVBXizd5v4rD/YVi0sNB1uTpV5YADhUIBCwvtHSSRSJCfn1/hslLp4ztxy9bLz8+vth+F\nQgEbGxutdgsLC632J/tWKBSa6ZCQELz44ov47LPPUFpailWrVqFPnz4VblNwcHCF8728vODl5VVh\nG2OMVUWlApYvBzZsEKYnTgS2bgVamqvw7v/m45vz3wAA1g9djyUDlsDAwECH1WqLi4tDXFxcg39u\nleEjFovx8OFDrXl5eXmQSCTVLpuXlwdACJnK+ikLJIlEUuXnVNS3WCwGADx69AhDhgzBunXrEBAQ\ngL///hsTJkyAjY0N5syZU67OysKHMcb+i+xswM8PiI0FjIyATz4BFiwAlCUK+Oz2w9EbR2FqZIrv\nfb7HJNdJui63nOr+8F69enW9fG6Vp926d+8OlUqF5ORkzbyEhAS4urqWW9bFxQWXLl3SWs7W1haW\nlpaV9uPi4qJZNyEhQdOmVCqRkpKi1f5k32U1JCUlIT8/H6+99hoMDQ3Rtm1bTJo0CceOHXuqHcEY\nY0/r3Dmgd28heGxsgJgYYOFC4G/FXQzaNghHbxyFlZkVYqbE6GXw6FR1F4X8/PzI39+flEolnTp1\niqRSKclksnLLRUVFkZ2dHclkMpLL5TR48GBaunRpjfopG+22f/9+KigooMWLF5Onp6dm3dDQUHJy\ncqKMjAy6c+cOOTs7U1hYGBERyeVyEolEFB4eTqWlpfT333+Th4eH1uAHIh5wwBirO2o1UWgokamp\nMLDA05Pozh2h7UrmFWq/sT0hGNT1i650Leeaboutpfr67qy2R7lcTj4+PiQSiahjx44UERFBRERp\naWkkFovp9u3HQwU3btxItra2ZGFhQdOnT6fi4uJq+ykTHR1Njo6OZGZmRt7e3pSWlqbVvmTJErKy\nsiIrKysKDAzUajt27Bg9++yzZGFhQXZ2dvTWW29RQUGB9oZy+DDG6sCjR0TTpj0eRv3uu0RFRULb\nzyk/k8XHFoRgkOf/eVKWIku3xdaB+vruNPin8yav7AJfM9lcxlg9uHULGD8e+PNPwMwM+O474Vlt\nALDtz21468hbUKlVmOA8Ad/7fN8kHg5aX9+d/GBRxhirgago4W2jf/4JdO0KnDkjBA8RYeUvKzH9\n0HSo1Cosfn4x9kzY0ySCpz7xs90YY6wKajWwbh2wapVwom3MGGDnTqBVK6BIVYSZh2fih8s/wNDA\nEF+99BXe7ve2rktuFDh8GGOsErm5wttGjx4V3ja6Zg3wwQeAoSGQW5CLcXvHIS41DiITEfZM2IPR\n3UfruuRGg8OHMcYqkJAgvG305k3hbaPh4cDIkULbrdxbGBU+Cn/l/AV7sT2OBBxBb/veui24keHw\nYYyxJ/zwA/DWW8JL33r3Fl6D0Lmz0HYu4xxejngZWcosuNq44mjAUXSQdtBtwY0QDzhgjLF/FBcD\n774rnGorKADeeAOIj38cPAf/Ogiv7V7IUmZhWJdhiH8jnoPnP+IjH8YYA5CRITyT7fRpwNQU+Ppr\nYOZM4VoPAHxx5gssOL4ABMIbvd5A2JiwZvG66/rC4cMYa/bi4oS3jWZlAe3bC6fZnntOaCtVl2Lh\niYX48uyXAIA13muw7IVlevVw0MaIw4cx1mwRARs3AoGBQGmp8NbRiAigTRuhXVmsxOTIyfjp2k8w\nMTTB1rFb8Zrba7otuong8GGMNUv5+cCMGcCPPwrTS5cKQ6mNjITpTEUmXo54Gefvnkerlq1wcNJB\nDO40WHcFNzEcPoyxZufGDWDsWODqVUAiAXbsAHx9H7dfzb6KUeGjkPogFZ1adcKxgGNwauOku4Kb\nIA4fxlizcuoU4OMDyOWAszNw4ADQvfvj9rjUOPju8cWDwgd4ru1zOOR3CLZiW90V3ETxUGvGWLOx\nc6dwXUcuFx6Tc+aMdvD8cPkHDN85HA8KH8DH0Qe/TP2Fg6eecPgwxpo8ImDlSmDKFKCkBJg/Hzh4\nUDjlJrQT1pxcg9cPvI4SdQne83gP+ybug7mJuW4Lb8L4tBtjrEkrLASmTxdGsRkaAl9+CbzzzuP2\nktISzDoyC9subYMBDLBp5CbM6z9PdwU3Exw+jLEmKztbuL7z+++AWAzs3Qu89NLj9rzCPEz4cQKi\nb0bDzNgMEeMjMNZxrO4KbkY4fBhjTdLVq8Do0cIL4Nq3B44cAdzcHren56Vj1K5RSMpOgq3IFof9\nD6Nf2366K7iZ4fBhjDU5MTHCG0fz8oC+fYFDhwB7+8ft5zLOYezusbinuAcnayccm3wMnVp10lm9\nzREPOGCMNSlbtgivPsjLE16JcPKkdvD838X/wwvbXsA9xT14dfLCb9N/4+DRAQ4fxliToFYDQUHC\nw0BVKmDJEuHpBeb/DFgrUhVh1pFZePPwmyguLcY7/d7B8deOw9LMUreFN1N82o0x1ug9eiQMo96/\nHzA2Br75BnjzzcftGQ8zMOHHCThz5wxaGLVA2JgwTO01VXcFMw4fxljjdu8e8MorwPnzgFQqPJF6\n2LDH7afSTmHijxORqcxEB2kHRL4aiT7P9NFdwQwAhw9jrBG7ckV4UkF6uvDCt6NHAad/HsFGRPj6\n3NdYeGIhVGoVhnQegt3jd6ONqI1ui2YA+JoPY6yRiooCBgwQgsfTU3hUTlnwFJQUYOrBqZgXNQ8q\ntQqLnl+E468d5+DRI3zkwxhrdL75Bpg7VxhkMGkSsG0bYGYmtKU+SMW4PePw570/YW5ijq2vbMUk\n10m6LZiVw+HDGGs0SkuBRYuATZuE6eXLgdWrhcfmAED0zWj47fPD/YL76GLZBQcnHURP2566K5hV\nisOHMdYoKBRAQABw+DBgYgJs3gxM/WfAGhHhs98/Q1BMENSkxksOL2HXuF08jFqPcfgwxvTenTvA\nyy8Dly4BlpbCO3gG//NSUUWxAjMOzcDepL0AgOWDliN4cDCMDI10WDGrDocPY0yv/fmnMKLt7l3A\nwUEY0Vb2Dp5keTJ89/giMSsRElMJvvf9Hj6OProtmNUIhw9jTG8dPgz4+wNKJfDCC8IRT+vWQtvR\n60cxOXIy8ory4GjtiAOTDsDR2lG3BbMa46HWjDG9QyQMKhg7Vgie118Hfv5ZCB41qbHm5Bq8HPEy\n8ory4OPog7Mzz3LwNDJ85MMY0ysqlfCm0W++EabXrAGWLQMMDIT370w5OAWHrh2CAQyw1nstlr6w\nFIYG/Hd0Y8PhwxjTGw8fCvftREUBLVoI9+/4+wttV7OvwmePD67fv45WLVshfFw4Xur2UtUdMr3F\n4cMY0wtpacLAgsREwNoa+Okn4PnnhbbIq5GYenAqFMUK9LTpiQOTDqCrVVfdFsxqhY9VGWM6d/48\n0L+/EDyOjsDZs0LwlKpL8UHMBxi/dzwUxQr4ufrh9IzTHDxNQLXhI5fL4evrC7FYjE6dOiEiIqLS\nZUNCQmBvbw+pVIoZM2aguLi4xv3ExMTA0dERIpEIQ4YMQXp6ulZ7YGAgrK2tYW1tjaCgoHKf/cUX\nX6BLly4Qi8VwdnbGjRs3qt14xpju7d8v3LOTmQkMHQqcPg106QLIC+QYHT4aH8d/DCMDI3w+/HOE\njwuHyFSk65JZXaBq+Pn5kZ+fHymVSoqPjyepVEpJSUnllouKiiJbW1uSyWSUm5tLXl5eFBQUVKN+\nsrOzSSqV0r59+6ioqIgWL15MHh4emnVDQ0OpR48elJGRQRkZGeTs7EyhoaGa9s2bN5ObmxtdvXqV\niIhu3rxJcrlcqz4AVIPNZYw1ELWaaMMGImFsG9GMGUTFxULbpb8vUedNnQnBIOtPrCn2Zqxui23G\n6uu7s8oeFQoFmZqa0o0bNzTzpkyZohUqZfz9/WnZsmWa6djYWLKzs6tRP2FhYTRgwABNm1KpJDMz\nM7p27RoREXl6etLmzZs17Vu3btWEU2lpKbVr145iY6v+5eTwYUx/FBcLYVMWPBs2CGFERLTr8i4y\nW2tGCAb1CetDaQ/SdFtsM1df351Vnna7fv06jI2N4eDgoJnn7u6OpKSkcsvKZDK4u7trpt3c3JCZ\nmYnc3Nxq+0lKStJa19zcHA4ODpr2ivoua7tz5w4yMjJw5coVdOjQAV26dEFwcDCEfcYY0ze5ucDI\nkcCWLUDLlsLL35YsAUpJhYXHF2Jy5GQUqAowrdc0nHrjFDpIO+i6ZFYPqhztplAoYGFhoTVPIpEg\nPz+/wmWlUqlmumy9/Pz8avtRKBSwsbHRarewsNBqf7JvhUIBQAgfAPj555+RmJiI3NxcDB8+HO3a\ntcPMmTPL1RkcHFzhtnp5ecHLy6vCNsZY3bh5Exg9GvjrL8DWFjh0CHjuOSBLmYVJ+yYhLjUOxobG\n+HLkl5jddzYMDAx0XXKTFxcXh7i4uAb/3CrDRywW4+HDh1rz8vLyIJFIql02Ly8PgBAylfVTFkgS\niaTKz6mob7FYDAAw++clHkuWLIGFhQUsLCwwa9YsHDt27KnChzFWv37/XXhiQU4O4OoKHDkCdOwI\nnM84j3F7x+HOwzuwE9th38R9GNBhgK7LbTaq+8N79erV9fK5VZ526969O1QqFZKTkzXzEhIS4Orq\nWm5ZFxcXXLp0SWs5W1tbWFpaVtqPi4uLZt2EhARNm1KpREpKilb7k32X1dCjRw+YmpqWq4f/YmJM\nf0REAEOGCMEzciTw229C8Gz9cyte2PYC7jy8g+fbP48/3vqDg6e5qO6ikJ+fH/n7+5NSqaRTp06R\nVColmUxWbrmoqCiys7MjmUxGcrmcBg8eTEuXLq1RP2Wj3fbv308FBQW0ePFi8vT01KwbGhpKTk5O\nlJGRQXfu3CFnZ2cKCwvTtE+ZMoXGjBlD+fn5dPv2bXJ0dKStW7dq1QcecMBYg1OridaseTywYM4c\nopISoiJVEc0+MpsQDEIwaM6ROVSkKtJ1uawC9fXdWW2PcrmcfHx8SCQSUceOHSkiIoKIiNLS0kgs\nFtPt27c1y27cuJFsbW3JwsKCpk+fTsVl4yar6KdMdHQ0OTo6kpmZGXl7e1NamvYIlyVLlpCVlRVZ\nWVlRYGCgVtvDhw/Jz8+PJBIJtW/fntasWVN+Qzl8GGtQRUVE06YJoWNgQBQSIoRRxsMMen7L84Rg\nUIs1LWjLxS26LpVVob6+Ow3+6bzJKzsN10w2lzGdevAAGD8eiI0FzM2B8HDhes9v6b9hwo8TcE9x\nD+0s2iHy1Uj0a9tP1+WyKtTXdyc/240xVqdSU4URbTKZMKLtyBGgTx/CN+e/xfyo+VCpVfDq5IU9\nE/bARmRTbX+saeLwYYzVmfPnhdddZ2YCLi7CW0dt2xZi+qE52H5pOwBggccCfPLiJzA25K+f5oz/\n32eM1YmDB4GAAKCgABg2DPjxR+BOcSLGbPZHYlYizIzN8H+v/B8CegboulSmB/ip1oyxWil76+i4\ncULwTJ8OHD1K2HntK/T9ri8SsxLRzaobTs84zcHDNPjIhzH2n6lUwIIFwNdfC9Nr1wIz5mVi7I/T\nEJUcBQCY2XsmQkaEQGwq1mGlTN9w+DDG/hOFQnjL6JEjgKkpsH07YNHnKNxC30D2o2xYmVlh88ub\nMc5pnK5LZXqIh1ozxp7a3bvCW0f//BOwsgL27C/AwUeL8f/O/z8AwNDOQ7HDZwfaWrTVcaWstnio\nNWNML1y5Igylvn0b6NoV2LgrAfP/CIAsWwYTQxN8NPQjLPRcCEMDvqTMKsfhwxirsRMngAkTgPx8\nwPN5NUat/gITfw5CcWkxHK0dET4uHM/aP6vrMlkjwKfdGGM1snkzMGcOUFoKvDL5bzwaPg3Rt04A\nAGb3nY3Ph38OcxNzHVfJ6lp9fXdy+DDGqqRWA8uWAevXC9PjPvgJJ6UzcL/gPlqbtcaWV7ZgrONY\n3RbJ6g2HTy1x+DD29AoLgalTgb17AcMWjzBwzUL8+igMADC863BsH7sd9hJ7HVfJ6hOHTy1x+DD2\ndHJyhIeB/v47YN71IqzeDMCdwmswNTLFhmEbMK//PB5U0Axw+NQShw9jNXf9OjBqFJByUw3pS59D\n2X8ZVFQC5zbOCB8XDnc7d12XyBoIh08tcfgwVjOnTgE+PoC8JAPiKVOgaBMLAHj3uXfxybBPYGZi\npuMKWUPi+3wYY/UuPBx44w2guGskTMa/CYWxHDYiG2x9ZStGdx+t6/JYE8LhwxgDEfDRR8DyDxXA\nyPeA3ltQAuAlh5ewbew22IptdV0ia2L4tBtjzVxJCTBrFrDt+Hlg/GSg9Q20MGqBz4Z/hnf6vaP5\nb4c1T3zNp5Y4fBgr78EDYPyEUsQWfQJ4rwSMVOhp0xPh48PhauOq6/KYHuBrPoyxOpWWBrw4/jZu\n9Hwd6HQSADC//3ysH7YeLY1b6rg61tRx+DDWDF24AAybvxd5w2YBZg9g3dIWO8dvx0iHkboujTUT\nfNqNsWZmd2Q+Xgufh9Ke2wEAIzqPwc7xW9FG1Ea3hTG9xNd8aonDhzFg4ednEZI+GbBKgZG6JUJe\n2oh3+8/mQQWsUhw+tcThw5qz4pJSDF7+Mc60DAYMS2EHd8TMiYCzjZOuS2N6jgccMMb+k6t30zBw\n42uQS+IBAKNavY/Id9ahhXELHVfGmjN+KiBjTdi3pyLg+o0b5JJ4GCjs8bn7zzg6/zMOHqZzfOTD\nWBOUV5iH18LfxZHbPwAmgCjdBzHvbUb/nta6Lo0xABw+jDU5v9/+HeN2TUZmUSpQbI4uNzbh9Ncz\nYWPDgwqY/uDwYayJKCgpwNpTa/HxqfUgqIG7vTFSGY7IHT1gxg+iZnqGw4exJiD6ZjTmHJmD5Nxk\ngAyA3wIrqOYGAAAgAElEQVSxqM+H2PCtKQz5yi7TQxw+jDViWcosvH/iffxw+Yd/ZrjA8GgYvv1g\nAN56S7e1MVYVDh/GGiE1qbHtz21Y/PNi5BbmwkDVEhS3CtY3FuLH3abw8tJ1hYxVjcOHsUZGli3D\n7COzcSr9FADA4OZw0OFv8Fy3rth/AWjXTscFMlYDHD6MNRIFJQX4KP4jbIjfgBJ1CcxKbVBw4AtQ\n4iS89ZYBvvwSaMG377BGgsOHsUYg+mY05hydg2R5MgDA5vZbyNq1HqZqS/y/zcDMmToukLGnVO04\nGLlcDl9fX4jFYnTq1AkRERGVLhsSEgJ7e3tIpVLMmDEDxcXFNe4nJiYGjo6OEIlEGDJkCNLT07Xa\nAwMDYW1tDWtrawQFBVX4+SdPnoShoSFWrFhR3WYx1ihkKbPw+oHX8eLOF5EsT0ZnkQssD8Qja0sY\n2llbIj6eg4c1TtWGzzvvvIOWLVsiKysLu3btwpw5cyCTycotd/z4cWzYsAGxsbFIS0vDzZs3sWrV\nqhr1k5OTg/Hjx2PdunXIzc1F3759MWnSJM26YWFh+Omnn3D58mVcvnwZhw8fRlhYmNbnl5SUYP78\n+fDw8OAn9LJGT01qbLm4BY5fO+KHyz+gpXFLvGz2MdKXXURuwgB4ewN//AH066frShn7j6gKCoWC\nTE1N6caNG5p5U6ZMoaCgoHLL+vv707JlyzTTsbGxZGdnV6N+wsLCaMCAAZo2pVJJZmZmdO3aNSIi\n8vT0pM2bN2vat27dSh4eHlqf//HHH1NgYCBNmzaNli9fXq4+AFTN5jKmF5KykuiFrS8QgkEIBg3d\nPpxemZZMABFAtGgRUUmJrqtkzUV9fXdWeeRz/fp1GBsbw8HBQTPP3d0dSUlJ5ZaVyWRwd3fXTLu5\nuSEzMxO5ubnV9pOUlKS1rrm5ORwcHDTtFfX97xrS0tKwbds2rFixgl+ZwBqtQlUhVvyyAr1Ce+FU\n+inYiGzwxQsRyA6JwqHtXSESAXv2AJ9+Chjz1VrWyFX5K6xQKGBhYaE1TyKRID8/v8JlpVKpZrps\nvfz8/Gr7USgUsLGx0Wq3sLDQan+yb4VCoZmeN28e1q5dC5FIBAMDgypPuwUHB1c438vLC158cwTT\nkScHFLzV5y0MpfWYNd4SDx4A3boBBw4ALi46LpQ1OXFxcYiLi2vwz60yfMRiMR4+fKg1Ly8vDxKJ\npNpl8/LyAAghU1k/ZYEkkUiq/JyK+haLxQCAw4cPQ6FQYOLEiQCEFx5VdfRTWfgwpgtPPqHApY0L\nvh0dhpM7B8BvpXCi7ZVXgO+/B/719xdjdaa6P7xXr15dL59b5Wm37t27Q6VSITk5WTMvISEBrq6u\n5ZZ1cXHBpUuXtJaztbWFpaVlpf24/PNnnIuLCxISEjRtSqUSKSkpWu1P9l1WQ2xsLC5cuAB7e3vY\n29tj79692LRpE3x9fZ9qRzDWkCoaUPDRkI/wi99FfDZvAMoGbK5ZIxzxcPCwJqe6i0J+fn7k7+9P\nSqWSTp06RVKplGQyWbnloqKiyM7OjmQyGcnlcho8eDAtXbq0Rv1kZ2eTVCql/fv3U0FBAS1evJg8\nPT0164aGhpKTkxNlZGTQnTt3yNnZmcLCwoiIKD8/nzIzMykzM5Pu3btHkyZNooULF1Jubq5WfeAB\nB0xPPDmgYPjO4ZR8P5kSE4m6dRMGFbRqRXTsmK4rZaz+vjur7VEul5OPjw+JRCLq2LEjRUREEBFR\nWloaicViun37tmbZjRs3kq2tLVlYWND06dOpuLi42n7KREdHk6OjI5mZmZG3tzelpaVptS9ZsoSs\nrKzIysqKAgMDK6132rRptGLFivIbyuHDdKygpICWxy4nkw9NCMEgm09tKOJKBKnVatq7l0gkEoLH\nzY0oOVnX1TImqK/vToN/Om/yygYhNJPNZXqmogEF64euh8TEEkuXAp99JiwXEABs3gyYm+uwWMb+\npb6+O3nAJmP1KFuZjYUnFmoNKAgbE4YBHQYgOxsY4QfExgpDpz//HJg7F+B7pFlzwOHDWD0oe+XB\nkuglkBfI0dK4JVYOWon3n38fpkamOH8eGD8euH0bsLUF9u4FBg3SddWMNRwOH8bq2NXsq5h1ZJbm\nlQfDuw7HN6O+QVerrgCALVuAt98GiosBDw9g3z6gbVtdVsxYw+PwYayOFKoK8dGpj7A+fj1K1CWw\nEdlg04hN8HP1g4GBAYqKgHnzgO++E5afMwfYtAkwNdVt3YzpAocPY3Ug5mYMZh+dXW5AgaWZJQDg\nzh3hNNu5c8I7d0JDgWnTdFgwYzrG4cNYLWQrs/H+ifex8/JOANoDCsrExQGvvgpkZwMdOgCRkUCf\nPjoqmDE9weHD2H+gUquw7c9tCIoJqnBAASA8GmfTJmDxYqC0FBg2DIiIAKytdVw8Y3qAw4exp0BE\nOPDXASyLXYa/cv4CUH5AAQAolcJL3nbvFqYDA4F16wAjI11UzZj+4fBhrIbiUuMQFB2EsxlnAQBd\nLLvgoyEf4VWXV7WepJ6cDPj6AomJgFgMbN8uXO9hjD3G4cNYNS7du4SlMUsRlRwFALAV2WLl4JWY\n2Xum5hRbmaNHgcmTgbw8oEcP4aGgTk66qJox/cbhw1glUuQpWPHLCkQkRgAAJKYSLBmwBO95vAex\nqVhrWbVaeAJ12Rs7fHyAHTuAJ15jxRj7B4cPY0/IVGRiza9rEPZHGFRqFUyNTPFOv3fwwQsfwNq8\n/GiBBw+A114TjnoMDIRrO4GBgGGVLyxhrHnj8GHsHw+LHuLT3z9FyOkQKEuUMDQwxLRe07DaazU6\nSDtUuM6VK8L1nZQUwMoKCA8HRoxo4MIZa4Q4fFizV6gqxLfnv8W6U+twv+A+AOCVHq/goyEfwcWm\n4vdWEwG7dgGzZgGPHgG9egn373Tu3JCVM9Z4cfiwZqtUXYqdl3diVdwqpOelAwAGdhiI9UPXa90k\n+qTkZODdd4Hjx4Xp118HwsIAM7OGqJqxpoHDhzU7RITD1w/jg5gPkJSdBADoadMTHw/9GKO6jdIa\nNv1vRUXAhg3ARx8J/7tVK+DTT4EZM/g1CIw9LQ4f1qzEp8cjKDoIv93+DQDQqVUnfOj1IQJ6BsDI\nsPI7QKOjhSdR37ghTE+ZIgSPjU1DVM1Y08Phw5qFK5lXsDRmKY7eOAoAaGPeBssHLcesPrPQwrhF\npev9/TewcOHjJxU4OQHffgsMHtwQVTPWdHH4sCYt9UEqVv6yEj9c/gEEgthUjPc938f7nu9D0kJS\n6XqlpcA33wDLlwMPHwrXc1auFIKIX4HAWO1x+LAmKVuZjbWn1uLb89+iRF0CE0MTzOk3B8teWAYb\nUdXnys6fB2bPBi5eFKbHjAG++gro1Kn+62asueDwYU1KflE+Np7eiM9OfwZFsQIGMMBrbq/hQ68P\n0dmy6nHQDx4AH3wgvGuHCGjfHvjyS2DsWB5QwFhd4/BhTUKRqgjf/fEd1vy6BtmPsgEAo7qNwsdD\nP4abrVuV65bds/P++0BWFmBsDCxYIJxmE4urXJUx9h9x+LBGTU1qhF8Jx4pfViD1QSoAwLOdJ9YP\nW49BHQdVu/5ffwmj2H75RZgeOFAYUODqWo9FM8Y4fFjjREQ4duMYPoj9AJczLwMAnNs446MhH+GV\nHq9Ueq9OmYIC4Rlsn3wClJQArVsLQ6enTuVnsjHWEDh8WKNz+vZpBEYH4lT6KQBAe4v2WO21GlPc\np1R5r06ZY8eEJxTcuiVMz5wJrF8vBBBjrGFw+LBGQ5YtwwcxH+Cnaz8BAKzMrLDshWV4u9/baGnc\nstr1b98G3ntPeAYbAPTsKZxiG1D5k3QYY/WEw4fpvaSsJHz6+6fYeXkn1KSGuYk5FnouxCLPRZC2\nlFa7fkmJMGpt1Srh9dYiEbB6NTBvHmBi0gAbwBgrh8OH6SUiwvGU4wg5E4ITKScAAMaGxpjddzZW\nDFoBO7Fdjfr5/Xfhnp0rV4TpceOATZuEYdSMMd3h8GF6paCkADsv78SmM5twNecqAMDcxBzTek3D\nQo+F6GrVtUb93L8vvNBtyxZhunNn4OuvgVGj6qtyxtjT4PBheuHv/L/xzYVv8O35bzXv1GkraYu5\nz83Fm33ehJWZVY36UauB7duBJUuEADIxEf73Bx8A5ub1uAGMsafC4cN06tK9Swg5E4KIKxEoUZcA\nAPo+0xcLPBZgovNEmBjV/KLMlSvAnDnAb8IDq+HtLTyfzdGxPipnjNUGhw9rcGpS4+j1o9h4ZiPi\nUuMAAAYwwDincVjgsQAD2g+o9j6df1MogA8/BDZuFB4IamMj/O+AAH4sDmP6isOHNRhlsRLbL23H\nF2e/wA258GIcsakYM56dgXn956GLZZen6o8IOHgQmD9fGEZtYCAc+axbB1ha1scWMMbqCocPq3d3\nHt7B1+e+xnd/fIfcwlwAQAdpB8zvPx8znp1Ro+HST7p1C5g7FzgqvJ4HvXsL9+w891xdVs4Yqy/V\nPkhELpfD19cXYrEYnTp1QkRERKXLhoSEwN7eHlKpFDNmzEBxcXGN+4mJiYGjoyNEIhGGDBmC9PR0\nrfbAwEBYW1vD2toaQUFBmvnZ2dnw9/dH27Zt0apVKwwcOBDnzp2r8Q5g9ed8xnkE7A9A5y86Y8Nv\nG5BbmAvPdp7YO2EvUualYKHnwqcOnuJi4OOPARcXIXgsLIR7eM6d4+BhrFGhavj5+ZGfnx8plUqK\nj48nqVRKSUlJ5ZaLiooiW1tbkslklJubS15eXhQUFFSjfrKzs0kqldK+ffuoqKiIFi9eTB4eHpp1\nQ0NDqUePHpSRkUEZGRnk7OxMoaGhRER08+ZNCgkJoXv37pFarabvvvuOrK2tSaFQaNUHgGqwuayW\nVKUq2i/bTwO3DiQEgxAMMlptRJN+nERnbp+pVd+//ELk5EQknHAj8vMjunu3bupmjFWsvr47q+xR\noVCQqakp3bhxQzNvypQpWqFSxt/fn5YtW6aZjo2NJTs7uxr1ExYWRgMGDNC0KZVKMjMzo2vXrhER\nkaenJ23evFnTvnXrVq1wepKFhQVdvHhRe0M5fOpVXmEehZwOoc6bOmtCR/qxlBadWERpD9L+c78q\nFVF0tBA0ZaHTrRvRiRN1WDxjrFL19d1Z5TWf69evw9jYGA4ODpp57u7uiIuLK7esTCaDr6+vZtrN\nzQ2ZmZnIzc1Fampqlf0kJSXB3d1d02Zubg4HBwckJSWhe/fukMlkWu1ubm5ISkqqsOZLly6huLhY\n67NY/Ul9kIovz36JLX9uwcOihwCArpZdMb//fEzrNa3KV1VXRq0Whkvv2QPs2wdkZgrzW7QQ7tdZ\nsgRoWf2j3BhjeqzK8FEoFLCwsNCaJ5FIkJ+fX+GyUunj8/dl6+Xn51fbj0KhgI2N9quNLSwstNqf\n7FuhUJSr4eHDh3j99dcRHBwMiaTiL73g4OAK53t5ecHLy6vCNqaNiHD6zmmEnAlB5NVIqEkNABjU\ncRAWeizEmO5javR0ae0+hes2e/YAe/cCGRmP27p2Bfz8gOnTgS5PNyCOMVaNuLi4Cg8o6luV4SMW\ni/Hw4UOteXl5eRV+sT+5bF5eHgAhZCrrpyyQJBJJlZ9TUd/iJ14xWVBQgJdffhnPP/88AgMDK92m\nysKHVa+ktAT7r+5HyJkQnMsQBnUYGxojoGcAFngsQG/73k/VHxFw6RKwe7cQOKmpj9s6dAAmTRJ+\nevfm+3UYqy/V/eG9evXqevncKsOne/fuUKlUSE5O1pzGSkhIgGsFr3l0cXHBpUuXMGHCBM1ytra2\nsLS0hKmpaYX9uLi4aNbdsWOHpi+lUomUlBSt9kuXLqFv374V1lBUVAQfHx906NABYWFh/3lnsIo9\nKHyAzX9sxlfnvsLth7cBCK8zmN13Nt7p9w6ekTzzVP0lJQmBs2cPcOPG4/nPPANMnCgc5fTvz4HD\nWJNW3UUhPz8/8vf3J6VSSadOnSKpVEoymazcclFRUWRnZ0cymYzkcjkNHjyYli5dWqN+yka77d+/\nnwoKCmjx4sXk6empWTc0NJScnJwoIyOD7ty5Q87OzhQWFkZERMXFxTRmzBjy8fEhlUpV6XaABxw8\nteT7yTT32FwSrRNpBhH0+KoHfXv+W1IWK5+qr2vXiD78kMjF5fHAAYCoTRuiOXOITp4kKi2tpw1h\njP1n9fXdWW2PcrmcfHx8SCQSUceOHSkiIoKIiNLS0kgsFtPt27c1y27cuJFsbW3JwsKCpk+fTsXF\nxdX2UyY6OpocHR3JzMyMvL29KS1Ne4TUkiVLyMrKiqysrCgwMFAzPy4ujgwMDEgkEpFYLNb8xMfH\na28oh0+NqNVqirsVR2MjxpJBsIEmdIZ9P4yOXj9KpeqaJ8StW0Tr1xM9+6x24FhaEs2cSfTzz0Ql\nJfW3LYyx2quv706Dfzpv8sqeFdZMNvepybJlCL8SjvAr4bj1QHi/tKmRKSb3nIz3PN6Dm61bjfq5\ncwf48UfhlNrZs4/nW1gAPj7CNZxhwwBT0/rYCsZYXauv705+vE4zlp6Xjt2JuxF+JRwJmQma+c9I\nnsGbvd/EnL5zYCu2rbafzExhSPSePcCpU4/nm5sDr7wiBM7IkTw8mjH2GIdPM5PzKAf7ZPsQfiUc\np9IfJ0Wrlq0w0XkiAnoG4IUOL1Q7VPr+fSAyUhg4EBcn3JsDCAEzapQwaGD0aH6HDmOsYhw+zYCy\nWIlD1w4hPDEcUclRUKlVAICWxi3xSo9XMLnnZIzoOgItjFtU2c+DB8JTpPfsAaKjAZXQDUxMhKCZ\nNEk40qnkFivGGNPg8GmiSkpLcCLlBMITw3Hwr4N4VPIIAGBkYISRDiMR4BoAH0efap9AoFAAhw4J\ngRMVJTzYEwCMjIARI4TA8fHhVxgwxp4Oh08ToiY1fkv/DeGJ4fgx6UfN66gB4Pn2zyPANQATXSbC\nRmRTRS9AQYHwxOg9e4R/CwqE+QYGwttBJ00Cxo0D2rSpz61hjDVlzW6027PPEsRiQCTCU/375Dwz\nM/24CZKIcOneZXx/MRx7r0bgrvK2pq1DSxc8L5kMdyM/tCzoDIUCmp/8fFQ6/eABUFLy+DMGDBAC\nZ8IEwN5eBxvJGNOZ+hrt1uzCRxiyXhf9/ffg+ve/LVoAjx5VHAJVBUQu3cL9ZyLwyGEXyFr2uLAH\nHYArAcJPVs//vH39+gmDBiZOBNq3r4MdxhhrlDh8aqlsB54/T1AqhS/w2vxbWKiDjRBlAS57gZ7h\nQPvTj+c/ag3ja69Cmh4AS8XzkIgNIZE8Dj6xGE81LZEIR3aMMcbhU0t1vQNVqsdHLGWh9F+CrKjo\n8VFQRYFgLHqIawYHcaEoHFeU0VCjFABgbizCS519ENAzAKMdX0QLE5M62S7GGPs3Dp9aakxPOChS\nFSEqOQq7ruzC4euHUagSDrOMDY3xksNLCOgZgJe7vwyRqUjHlTLGmjp+wkETV6ouxa9pvyI8MRz7\nZPvwoPCBpm1Qx0GY3HMyxjuNR2vz1jqskjHG6gaHjw4RES7+fRHhieHYnbgbd/Pvatp62fVCgGsA\n/Fz90F7KV/wZY00Lh08Du5t/F9E3ozU/fyv+1rR1seyCgJ4B8Hf1h3MbZx1WyRhj9YvDp57lF+Xj\nZNpJ/HzzZ0TfjIYsW6bVbi+2x0SXiQhwDcBzbZ/715Bwxhhrujh86lhJaQnOZZzThM3ZjLOaZ6kB\ngMhEBK9OXhjWZRhe7PIinNs4c+AwxpodDp9aIiJczbmK6JvR+Pnmz4hLjYOiWKFpNzIwwvPtn8ew\nLsMwrPMw9G/XH6ZG/DIbxljzxuHzH9zNv4uYmzGIviVct/n3QAEAcLR2xItdXsSwLsMwuONgSFtK\ndVQpY4zpJw6fGii7blN2dPPkdRs7sZ3myGZol6FoZ9FOR5UyxljjwOFTgbLrNtE3oxF9Kxpn7pwp\nd91mcKfBmqMblzYufN2GMcaeAocPhOs2f+X8pRkkEJcah/zifE27kYERPNt5agYJ8HUbxhirnWYb\nPn/n/42YWzGawKnouk3ZqTSvTl583YYxxupQswuf96LeQ/TNaCRlJ2nNtxXZCmHTZRiGdh7KTxVg\njLF61OweLIpg4Z+y6zbDOguB42rjytdtGGPsCfxg0TqyYtAKDOsyDB7tPPi6DWOM6UizO/JpJpvL\nGGN1or6+Ow3rtDfGGGOsBjh8GGOMNTgOH8YYYw2Ow4cxxliD4/BhjDHW4Dh8GGOMNTgOHz0TFxen\n6xJqhOusW42hzsZQI8B1NhYcPnqmsfxCcp11qzHU2RhqBLjOxoLDhzHGWIPj8GGMMdbgqg0fuVwO\nX19fiMVidOrUCREREZUuGxISAnt7e0ilUsyYMQPFxcU17icmJgaOjo4QiUQYMmQI0tPTtdoDAwNh\nbW0Na2trBAUFabWlpqbC29sbIpEITk5OiImJqdHGM8YY041qw+edd95By5YtkZWVhV27dmHOnDmQ\nyWTlljt+/Dg2bNiA2NhYpKWl4ebNm1i1alWN+snJycH48eOxbt065Obmom/fvpg0aZJm3bCwMPz0\n00+4fPkyLl++jMOHDyMsLEzT7u/vjz59+kAul2PdunWYMGECcnJyarVjGGOM1Z8qw0epVCIyMhJr\n1qyBubk5BgwYgLFjx2Lnzp3llt2xYwdmzpwJJycntGrVCitXrsT27dtr1E9kZCRcXV0xfvx4mJqa\nIjg4GAkJCbh+/bqm70WLFuGZZ57BM888g0WLFmn6vn79Ov7880+sXr0aLVq0wLhx4+Dm5ob9+/fX\n4W5ijDFWp6gKFy9eJHNzc615n3/+Ob388svllnV3d6e9e/dqpnNycsjAwIDkcnm1/cybN4/efvtt\nrfaePXtSZGQkERFJpVI6d+6cpu3ChQskkUiIiCgyMpKcnJy01p07dy7NnTtXax4A/uEf/uEf/vmP\nP3WtyiMfhUIBCwsLrXkSiQT5+fkVLiuVPn7VdNl6+fn51fZTUbuFhYVW+5N9KxSKCtueXJcxxpj+\nqfJlcmKxGA8fPtSal5eXB4lEUu2yeXl5AISQqayfssCRSCRVfk5FfYvF4kprfPDgQbkwI36PD2OM\n6Y0qj3y6d+8OlUqF5ORkzbyEhAS4urqWW9bFxQWXLl3SWs7W1haWlpaV9uPi4qJZNyEhQdOmVCqR\nkpKi1f5k32U1uLi44ObNm5ojoSf7ZowxpoeqOy/n5+dH/v7+pFQq6dSpUySVSkkmk5VbLioqiuzs\n7Egmk5FcLqfBgwfT0qVLa9RPdnY2SaVS2r9/PxUUFNDixYvJ09NTs25oaCg5OTlRRkYG3blzh5yd\nnSksLEzT7uHhQYsWLaKCggLav38/tWrVinJycmpzOpIxxlg9qjZ85HI5+fj4kEgkoo4dO1JERAQR\nEaWlpZFYLKbbt29rlt24cSPZ2tqShYUFTZ8+nYqLi6vtp0x0dDQ5OjqSmZkZeXt7U1pamlb7kiVL\nyMrKiqysrCgwMFCrLTU1lZydncnAwIAMDAzI3t6eTp06pem3R48eZG5uXmm/rVu3ptatW5frtz5E\nRESQo6MjiUQi6tq1K8XHx9OtW7fIwMCAxGKx5mft2rU6qVMkEmnVYWRkpDV4Q1/2Z1V16tP+JCK6\nffs2jRkzhqysrMjOzo7effddUqlURKQ/+7OyGvVtX8pkMvL29iapVEoODg504MABTZu+7Muq6tTl\n/vzqq6+oT58+1KJFC5o2bZpWW2323a1bt8jLy4vMzc3J0dGRoqOja1RP3Q9h0IETJ05Qx44d6ezZ\ns0REdPfuXcrIyKDs7GyysLCgffv2UVFRES1evJg8PDw064WGhlKPHj0oIyODMjIyyNnZmUJDQxu8\nzrJfSLVaXeF6DV1nGYVCQWKxWBPkZUeo+rI/K6tT3/anr68vTZs2jYqKiujevXvUs2dP+vLLL/Xq\n97OyGvVpX5aUlFC3bt0oJCSE1Go1xcbGkkgkouvXr+vVvqyqTl3uz8jISDp48CDNmTNHK3xq+9+1\nh4cHvf/++1RYWKg585SdnV1tPU0ifDw9PWnr1q3l5oeFhdGAAQM000qlkszMzOjatWua9TZv3qxp\n37p1q9ZOb6g6y34hy/4armi9hqyzzPbt26lr166aaX3bn5XVqW/7s3v37vS///1PM7148WKaNWuW\nXu3PimqcPXs2paam6s2+vHLlConFYq15w4cPpxUrVujVvqyqTn343Vy+fLlW+NRm3127do1atGhB\nCoVC0z5o0KAaBWajf7ZbaWkp/vjjD2RlZaFbt25o37495s6di8LCQiQlJcHd3V2zrLm5ORwcHJCU\nlAQAkMlkWu1ubm6atoass0zHjh3Rvn17TJ8+Hffv39fMb8g6/23Hjh2YMmWKZlqf9mdVdZbRl/05\nYsQIhIeHo6CgABkZGfjf//6Hl156qVwdutyfFdU4cuRITbu+7MsnqdVqJCYmQiaTwc3NTTNfX343\n/11nUlISDAwMAOh2f9ITI39r8991UlISunTpApFIpGl3d3evUc2NPnwyMzNRUlKC/fv3Iz4+Hpcu\nXcKff/6JtWvXQqlU/uf7hxqyzjZt2uDChQtIT0/HH3/8gfz8fEyePFmzbkPWWSYtLQ2//vorpk6d\nqpmnT/uzqjr1bX8GBwcjMTERFhYWaN++Pfr164exY8fW6v62hqrR2tpab/Zljx49YGNjg08//RQl\nJSU4ceIEfv31Vzx69Kja+/30pU592J9lAVimNv9d1+Y+y0YfPmZmZgCAuXPnwtbWFq1bt8bChQtx\n7Nixau9Tqur+oYasUyQSoXfv3jA0NISNjQ2+/vprnDhxAkqlssHrLLNz50688MIL6Nixo2aePu3P\nqurUp/1JRBgxYgQmTpyIR48eIScnB3K5HIGBgXqzP6uqUZ/2pYmJCQ4ePIijR4/C3t4eGzduxKuv\nvop27drpzb6sqM6QkBBNnfqwP5888qnNvqvpfZYVafThY2lpiXbt2lXYVpv7hxqyzsqo1WoADVtn\nmct0cBkAAAICSURBVO+//17raKKsDn3Zn1XVWRld7M+cnBz88ccfePfdd2FiYgIrKytMmzYNx44d\n05v9WVWNldHV72bPnj0RFxeHnJwcREVFISUlBf3799ebfVlRnf/73/+QkpKC5557rtLlG3J/Pnnk\no7P7LGt15UpPrFy5kvr160dZWVkkl8tp4MCBtHLlylrfP9RQdZ49e5b++usvKi0tpZycHHr11Vdp\nyJAhOqvzt99+I5FIpHURkaj292M1VJ36tD/VajU988wztGHDBlKpVJSbm0s+Pj40efJkvdmfVdWo\nT/uSiOjy5ctUUFBASqWSPv30U+rSpQsVFxdTVlaWXuzL6urU5f5UqVRUUFBAQUFB9Prrr1NhYSGp\nVCqd3WfZJMKnpKSE3n77bWrVqhXZ2dnR/PnzqaioiIhqd/9QQ9RZWFhIERER1LlzZxKJRGRvb09T\np06lzMxMndU5a9YsmjJlSoVt+rQ/K6tT3/bnmTNnaODAgdSqVSuytramSZMmUVZWFhHpz/6srEZ9\n25eLFy8mS0tLEovFNGrUKEpJSdG06cu+rKpOXe7PVatWae6FLPtZvXo1EdX+PksvLy8yMzMjR0dH\niomJqVE9BkT80DPGGGMNq9Ff82GMMdb4cPgwxhhrcBw+jDHGGhyHD2OMsQbH4cMYY6zBcfgwxhhr\ncBw+jDHGGtz/B06Vwohznx0kAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x44eb610>"
       ]
      }
     ],
     "prompt_number": 119
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "opt_result"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 120,
       "text": [
        "array([  8.22370949e+09,   4.43265816e+00,   4.62628660e+03,\n",
        "        -3.31921953e-05])"
       ]
      }
     ],
     "prompt_number": 120
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}