{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scipy and the Vapor-Liquid Saturation Pressure from the van der Waals Equation of State\n",
    "In this notebook we present numerical evaluation of density and residual properties for given values of pressure and temperature using the van der Waals Equation of State. Then, using these calculations we present a numerical evaluation of the saturation pressure for liquid-vapor equilibra at a given temperature.\n",
    "\n",
    "> **Here you will meet the following topics:**\n",
    ">- Thermodynamics: *equation of state*, *calculation of*: *density*, *residual properties*, *saturation point pressure*\n",
    ">- Python: `scipy`: `roots`, `optimize`: `bisect`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Supose you have to solve the following problem in a thermodynamics course:\n",
    ">\"Calculate the saturation pressure for pure hexane at the temperature of 220 K according to the van der Waals equation of state\"\n",
    "\n",
    "Here, we will show how to implement a function to calculate density of a pure fluid at given T and P.\n",
    "Then, how to calculate the residual properties. And at last, how to calculate the saturation pressure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use numpy again"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will also use the Scipy package., more specifically, we will use the  optimize module.\n",
    "The following syntax for importing packages is slightly different from the previous one. Nonetheless, both of them work in the same way. Here, we import only a subset of the resources of scipy, the subset identified as optimize."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy import optimize as opt\n",
    "#from <-- keyword \"from\"\n",
    "     #scipy <-- name of the package\n",
    "           #import <-- keyword \"import\"\n",
    "                  #optimize <-- name of the desired module\n",
    "                           #as <-- keyword \"as\" to introduce an alias\n",
    "                              #opt <-- alias name"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the *van der Waals* equation of state (Reid et al., 1987, see Eq 3-6.1 and Table 3-5, line 1)\n",
    "\n",
    "$P=\\frac{RT}{\\bar{V}-b}-\\frac{a}{{\\bar{V}}^{2}}$\n",
    "\n",
    "$a=\\frac{27}{64}\\frac{{{Tc}^{2}}{R^{2}}}{Pc}$\n",
    "\n",
    "\n",
    "$b=\\frac{{R}{Tc}}{{8}{Pc}}$\n",
    "\n",
    "where $R$ is the gas constant \n",
    "\n",
    "$R = 8.3144598$ $\\mathrm{{J} {mol^{−1}} {K^{−1}}}$\n",
    "\n",
    "and critical properties for hexane are \n",
    "\n",
    "$Tc = 507.5$ $\\mathrm{K}$\n",
    "\n",
    "\n",
    "$Pc = 30.1$ $\\times {10^5}$ $\\mathrm{Pa}$\n",
    "\n",
    "So let's code this information, as in notebook 2:\n",
    "\n",
    "First, the universal gas constant ($R$) is assigned to the variable named `R`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#gas constant\n",
    "R = 8.3144598 #J.mol^−1.K^−1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we assign the values of the critical properties of hexane ($Tc$ and $Pc$) to variables `Tc` and `Pc`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#pure component critical point of hexane\n",
    "# Tc (K)\n",
    "Tc = 507.5 #K\n",
    "\n",
    "# Pc (Pa)\n",
    "Pc = 30.1*(10**5) #Pa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we create variables `a` and `b` for the EoS parameters $a$ and $b$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.4954993637433445\n"
     ]
    }
   ],
   "source": [
    "a = 27/64*((Tc**2)*(R**2)/Pc)\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0001752320742732558\n"
     ]
    }
   ],
   "source": [
    "b = (R*Tc)/(8*Pc)\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we are interested in calculating $\\bar{V}$ for a given substance, therefore we will calculate $\\bar{V}$ given $a$ and $B$, at some condition of $P$ and $T$. For that purpose we convert the explicit $P(T,\\bar{V})$ representation of the EoS to an implicit representation. We can then find a repreesentation that corresponds to a 3$^{\\mathrm{rd}}$ degree polynomial in $\\bar{V}$, with coefficients depending on $P$, $T$, $a$ and $b$. Therefore, for a given substance (represented by a pair of $a$ and $b$), at given $P$ and $T$, the density can be obtained by solving this polynomial for $\\bar{V}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "${P}{{\\bar{V}}^3}−({P}{b}+{R}{T}){{\\bar{V}}^2}+{a}{\\bar{V}}−{a}{b}={0}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Let's define the temperature and pressure at which we want to perform the density calculation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#calculate density at what temperature and pressure?\n",
    "T = 220. #K\n",
    "P = 1*(10**5) #Pa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, let's define the coefficients of the polynomials:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "c3 = P #coefficient for v^3\n",
    "c2 = -(P*b+R*T) #coefficient v^2\n",
    "c1 = a  #coefficient v^1\n",
    "c0 = -a*b  #coefficient v^0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The volume is finally obtained by finding the roots of the polynomial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.01701555  0.00124509  0.00020641]\n"
     ]
    }
   ],
   "source": [
    "v=np.roots([c3,c2,c1,c0]) #numpy  roots finds the roots of a polynomial\n",
    "#coefficients must be provided in a list form, constructed using the square brackets\n",
    "#result is the list of real and imaginary roots found.\n",
    "print(v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, after solving a 3$^{\\mathrm{rd}}$ degree polynomial, we may obtain three real valued roots. If three roots are found in the interval $(b,\\infty )$, the minimum  and maximum corresponds to liquid and vapor aggregation states, respectively, while the intermediate root does not correspond to a phisycally relevant state.\n",
    "\n",
    "According to the principle of minimum Gibbs energy, at given $P$, $T$ (and quantity of matter), the state with least Gibbs energy is the *stable equilibrium state*. Therefore we must compare the Gibbs energy of these *states* to reveal which corresponds to the *stable equilibrium state*. While we are not able to calculate absolute values of Gibbs energy using an equation of state, we are able to calculate relative values of Gibbs energy for both vapor and liquid states with a common reference associated with the ideal gas state. The relative scale of Gibbs energy we use here is the residual Gibbs energy, since both roots are at the same temperature and pressure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.000206407191559\n",
      "0.00124508874657\n",
      "0.0170155476961\n"
     ]
    }
   ],
   "source": [
    "[v_Liq, v_intermediate, v_Vap] = np.sort(v)\n",
    "print(v_Liq)\n",
    "print(v_intermediate)\n",
    "print(v_Vap)\n",
    "v=np.array([v_Liq, v_intermediate, v_Vap]) #numpy.array is a special type of list with which numpy methods (e.g. numpy. exp, numpy.log) can work"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The residual Gibbs energy\n",
    "For the van der Waals EoS, the expression for the residual Gibbs energy ($\\bar{G}^{Res}$) is:\n",
    "\n",
    "$${\\bar{G}^{Res}} = P\\bar{V}-RT-\\frac{a}{\\bar{V}} -RTln \\left (\\frac{P(\\bar{V}-b)}{RT} \\right )$$\n",
    "\n",
    "For the origin and discussion of residual thermodynamics properties, see the references below:\n",
    "> `Smith, van Ness, Abbott, 2005, 7th-pt-BR, Ch. (6.2, 6.3)\n",
    "> O'Connell and Haile, 2005, Ch. 4\n",
    "> Tester and Modell, 1997, Ch. 9.`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "220.0 K 100000 Pa [-2238.4396586   1483.96852119  -123.05374806]\n"
     ]
    }
   ],
   "source": [
    "Residual_G = (P*v)-(R*T)-(a/v) - T*(R*np.log(P*(v-b)/(R*T))) \n",
    "print(T, \"K\", P, \"Pa\", Residual_G)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here you can note that the lowest Gibbs energy yielded correspond to the liquid phase. It means that, at this pressure and temperature, the liquid phase is the stable equilibrium state."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculation of saturation pressure\n",
    "\n",
    "Now, the saturation pressure is the pressure $P$ that makes the difference in residual Gibbs energy between trial liquid and vapor phase equals zero at a given temperature. We cannot explicitate $P$ as a function of the difference in residual Gibbs energy, so we define a function $\\Delta G^{Res}$ of $P$\n",
    "and we will solve it numerically, i.e. trying various values of $P$ according to a updating methodology and convergence criteria.\n",
    "\n",
    "To automate that task we implement a function called **`def Delta_Residual_G(P):`**, which calculates the roots of volume, as before, and evaluates the difference between residual Gibbs energy in each of these states."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-2115.38591053\n"
     ]
    }
   ],
   "source": [
    "#define equilibrium criteria:\n",
    "#Residual_G_Liquid is equal to Residual_G_Vapor :. \n",
    "#mathematically:\n",
    "#Delta_Residual_G = zero :.\n",
    "#in practical implementation for a numerical solution:\n",
    "#Delta_Residual_G < Tolerance\n",
    "def Delta_Residual_G(T,P):\n",
    "    c3 = P #coefficient for v^3\n",
    "    c2 = -(P*b+R*T) #coefficient v^2\n",
    "    c1 = a  #coefficient v^1\n",
    "    c0 = -a*b  #coefficient v^0\n",
    "    \n",
    "    V = np.roots([c3,c2,c1,c0])\n",
    "    \n",
    "    #note that near the critical point, the numerical methods underlying the roots function may yield results with a complex part\n",
    "    #here we discard any volume root for which the complex part is significant (higher then a given tolerance of 1e9i)\n",
    "    V = V[ V.imag < 1e9 ].real #filter out complex roots\n",
    "    V = V[ np.where( V > b ) ] #filter out non physical roots\n",
    "    \n",
    "    Vl = min(V)\n",
    "    Vv = max(V)\n",
    "    #print(P, V) #if you uncomment this line you will see every iteration in the output section of this cell\n",
    "    return (P*(Vl-Vv) + (-a)*(1/(Vl)-1/(Vv))+T*(R*np.log((Vv-b)/(Vl-b))))\n",
    "\n",
    "#test for T=220 and close to atmospheric pressure\n",
    "print(Delta_Residual_G(220,100000))\n",
    "#negative value means G_liq is lower"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What follows is similar to the work of Topliss et al.,:\n",
    "> Topliss, Robert J., Dimitrios Dimitrelis, and John M. Prausnitz. 1988. “Computational Aspects of a Non-Cubic Equation of State for Phase-Equilibrium Calculations. Effect of Density-Dependent Mixing Rules.” Computers & Chemical Engineering 12 (5): 483–89. doi:10.1016/0098-1354(88)85067-1.\n",
    "\n",
    "To calculate saturation pressure for a given $T$, first we isolate the range of pressure where fluid phase equilibria can occur:\n",
    "That is the range of pressure where there are 3 valid roots of volume (keep in mind the isotherms presented in notebook 2). The maximum and minimum pressure that presents 3 valid roots for volume we will call `P_upper_plateau` and `P_lower_plateau` respectively. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> given\n",
    "> \n",
    "> $T$,\n",
    "> \n",
    "> find all possible\n",
    "> \n",
    "> $\\bar{V}^*$,\n",
    "> \n",
    "> so that\n",
    "> \n",
    "> $\\left(\\frac{\\partial P}{  \\partial \\bar{V}}\\right)_{T}=0$\n",
    "> \n",
    "> in\n",
    "> \n",
    "> $\\bar{V}^* \\in (b,\\infty)$.\n",
    "> \n",
    "> Then, find the corresponding\n",
    "> \n",
    "> $P^* = P(T,\\bar{V}^*)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Analytical development of the partial derivative yields:\n",
    "\n",
    "$\\left(\\frac{\\partial P}{  \\partial \\bar{V}}\\right)_{T}=-\\frac{RT}{(\\bar{V}-b)^2} + \\frac{2a}{\\bar{V}^3}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore in the stationary points, the following equality is true:\n",
    "\n",
    "$-\\frac{RT}{(\\bar{V}-b)^2} + \\frac{2a}{\\bar{V}^3}=0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It represents another 3rd order polynomial in $\\bar{V}$ (not the same from the problem of finding $\\bar{V}$ given $P$ and $T$), as follows:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$-{RT}{\\bar{V}^3} + {2a}{(\\bar{V}-b)^2}=0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$-{RT}{\\bar{V}^3} + {2a}{(\\bar{V})^2}+{2a}{(b)^2}+{2a}{(-2\\bar{V}b)}=0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Vli: 0.000251661242908\n",
      "Vvi: 0.00233425727852\n",
      "Psupi: 388842.000342\n",
      "Pinfi: 1e-09\n"
     ]
    }
   ],
   "source": [
    "#calculate each coeffcient of the polynomial-in-V representation of dP/dV=0\n",
    "def Plimits(T):\n",
    "    rho = R*T\n",
    "    mu = 2*(-a)\n",
    "    nu = - b*(a*(-4))\n",
    "    omega=- 2*(b**2)*(a*(1))\n",
    "    #and use the method \"roots\" from package \"numpy\" to calculate the roots of this polynomial and assing the solutions to the variable plateaus as an array\n",
    "    plateaus = np.roots([rho,mu,nu,omega])\n",
    "#    print(\"plateaus:\", plateaus)\n",
    "\n",
    "    #filter out values of plateu corresponding to non-physical volume (less than b or complex)\n",
    "    plateaus = plateaus[ np.where( plateaus > b ) ] #using the numpy function 'where\n",
    "\n",
    "#    print(\"plateaus:\", plateaus)\n",
    "    Vli = min(plateaus)\n",
    "    Vvi = max(plateaus)\n",
    "\n",
    "    #find the corresponding values of pressure for each filtered volume\n",
    "    Psup = R*T/(Vvi-b)-a/(Vvi*Vvi)\n",
    "    Psupi=Psup*.999\n",
    "\n",
    "    Pinf = R*T/(Vli-b)-a/(Vli*Vli)\n",
    "\n",
    "    #note that the algorithm cannot work with pressure less than or equal to zero, therefore fix minimium pressure to slightly above zero if it was negative\n",
    "\n",
    "    if Pinf > 0:\n",
    "        Pinfi = Pinf*1.001\n",
    "    else:\n",
    "        Pinfi = 0.000000001\n",
    "        \n",
    "    return [Vli,Vvi,Psupi,Pinfi]\n",
    "\n",
    "[Vli,Vvi,Psupi,Pinfi]=Plimits(220)\n",
    "\n",
    "print(\"Vli:\",Vli)\n",
    "print(\"Vvi:\",Vvi)\n",
    "print(\"Psupi:\",Psupi)\n",
    "print(\"Pinfi:\",Pinfi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we must find the solution $P$. Thus, we ask the package optmize to use the bissect method and solve the equilibrium criteria (equality of residual Gibbs free energy for both phases):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-2115.38591053\n",
      "-2115.38591053\n"
     ]
    }
   ],
   "source": [
    "#Here, bisect is a numerical method that works on a function of 1 variable\n",
    "#However note that we wrote Delta_Residual_G as a function of 2 variables: T and P\n",
    "#So, to implement a handle of a function of one variable (P) that the method can wok with\n",
    "#and that corresponds to the previosuly defined d_G_Res(T,P) at one given T\n",
    "#we do the following\n",
    "def Delta_Residual_G_at_220K(P): #P appears here as the 1 variable inut argument\n",
    "    return Delta_Residual_G(220,P) #Temperature apppears hard-coded here as 220 K\n",
    "\n",
    "print(Delta_Residual_G(220,100000))\n",
    "print(Delta_Residual_G_at_220K(100000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saturation pressure =  29756.11279369196  Pa ( 0.2975611279369196 bar )\n"
     ]
    }
   ],
   "source": [
    "P_sat = opt.bisect(Delta_Residual_G_at_220K,\n",
    "                   Pinfi,\n",
    "                   Psupi,\n",
    "                   xtol=1.0e-9,\n",
    "                   rtol=1.0e-9,\n",
    "                   maxiter=100,\n",
    "                   full_output=True,\n",
    "                   disp=True)\n",
    "print(\"Saturation pressure = \",P_sat[0],\" Pa (\", P_sat[0]/10**5, \"bar )\" )bbb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the first element (index 0) of the returning array of `optmize.bisect` we read the value of the saturation pressure\n",
    "in the second element we read additional information about the execution of the numerical method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      converged: True\n",
      "           flag: 'converged'\n",
      " function_calls: 36\n",
      "     iterations: 34\n",
      "           root: 29756.11279369196\n"
     ]
    }
   ],
   "source": [
    "print(P_sat[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Grid, Calc and Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "29756.112793691944\n"
     ]
    }
   ],
   "source": [
    "def calc_Psat(T_in):\n",
    "    givenT=T_in\n",
    "#Here, bisect is a numerical method that works on a callback function of 1 variable\n",
    "#it is called a callback function because we call bisectt, and tell it to make a call back to our function as needed to perform its purpose, i.e. of finding the solution to the function.\n",
    "#However note that we wrote the function Delta_Residual_G as a function of 2 variables: T and P\n",
    "#So, to implement a handle of a function of one variable (P) that the method can work with as its callback function\n",
    "#and which results corresponds to the previosuly defined Delta_Residual_G(T,P) at one given T\n",
    "#we do the following\n",
    "    def Delta_Residual_G_givenT(P): #P appears here as the 1 variable inut argument\n",
    "        return Delta_Residual_G(givenT,P) #T apppears here and is obained from the surrounding scope, see below the tests:\n",
    "    \n",
    "    [_,_,Pinfi,Psupi]=Plimits(givenT)\n",
    "    ans = opt.bisect(Delta_Residual_G_givenT, Pinfi, Psupi, xtol=1.0e-9, rtol=1.0e-9, maxiter=100, full_output=True, disp=True)\n",
    "    return ans[0]\n",
    "\n",
    "#test\n",
    "print(calc_Psat(220))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "T=100\n",
    "grid_T=np.linspace(100,Tc,100)\n",
    "grid_P=np.zeros(100)\n",
    "grid_P[0]=calc_Psat(grid_T[0])\n",
    "for i in range(1,100):\n",
    "    grid_P[i]=calc_Psat(grid_T[i])\n",
    "    #print(i,grid_T[i],grid_P[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#compare with Antoine fits:\n",
    "def PAntoine_nist(T): #http://webbook.nist.gov/cgi/cbook.cgi?ID=C110543&Mask=4\n",
    "    if T < 177.70:\n",
    "        return 0\n",
    "    elif T < 264.93:\n",
    "        A=3.45604\n",
    "        B=1044.038\n",
    "        C=-53.893\n",
    "        return (10**( A - (B / (T + C)) )   )* 10**5\n",
    "    elif T < 286.18:\n",
    "        return 0\n",
    "    elif T<342.69:\n",
    "        A=4.00266\n",
    "        B=1171.53\n",
    "        C=-48.784\n",
    "        return (10**( A - (B / (T + C)) )   )* 10**5\n",
    "    else:\n",
    "        return 0\n",
    "    \n",
    "def PAntoine_ddbst(T): #http://ddbonline.ddbst.com/AntoineCalculation/AntoineCalculationCGI.exe\n",
    "    A=7.01051\n",
    "    B=1246.33\n",
    "    C=232.988\n",
    "    return (10**( A - (B / (T - 273 + C)) )   )* 101325/760.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkcAAAGHCAYAAAC6Q/beAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzsnXd4FOXWwH+HEHoHAUEQUECQImDDXsGKveD1KpaL7Vqw\nl3sVuV67gl2v2LBwbdcPRBQBQaQrIi30jjQRCKQRkpzvj3cWh3WTbMJuZjc5v+fZZ7MzZ+Y9mZ2Z\nPXPaK6qKYRiGYRiG4agUtAKGYRiGYRiJhBlHhmEYhmEYPsw4MgzDMAzD8GHGkWEYhmEYhg8zjgzD\nMAzDMHyYcWQYhmEYhuHDjCPDMAzDMAwfZhwZhmEYhmH4MOPIMAzDMAzDhxlHRkwRkRNFpEBETgha\nF6P0iMhA73tsELQuRvlHRO4VkbSwZQUi8nBQOpUVoWutlNsOF5GPY62TYcZR3BGRziLymYisEpFs\nEVknIt+KyN9Lub+bROTqWOsZYz1sTprkR7Hv0SgDRKQ2cC/wZNC6BMS+XGtPAReJSOcY6mNgxlFc\nEZFjgB+BzsB/gFuAN4F84LZS7vZmIHDjiEL0UNXvgeqqOqnsVTIMIwm5DkgB/hu0IsmGqv4C/ATc\nFbQu5Y3KQStQznkI2A4crqo7/StEpFEwKkVGRKqpak4s9qWqubHYT1kiIjVUNStoPUqCiFQFctVm\njzaSm37AyGS8byQInwADReTmZLuHJTLmOYovbYAF4YYRgKpu8X8WkWtEZLyIbBKRHBFZICI3hsms\nBA4FTvLi8QUi8p23LmLcWkT6eXItfctWichIEeklIj+KSDbQP0Z6RMw5EpFLROQnEckSkd9E5H0R\naRYm866I7BSRZiLyf97fm0XkGRGR4g627/86XURme2HMBSJyQZjc1SEdReRVEdkErPWtbyYib4vI\nRu8YzBeRayKMd6u3LlNEtnrH8nLf+loiMkREVnr72eSFVA8L0/ntCPueGDqmYcf1MhF5TETWAZlA\nbW99XW+sNd5YS708jmKPWzHU976XbSKy3Tsu1SLoe6Xv+/3dy4U4wLf+EG/du2HbHScieSLyhG9Z\nHxEZJSK/ev/LMhH5h4hUCtt2oojMFZEOIjLB+x7Wicg9EfSrIiKPesclxztOT4lIlX08PsY+ICKt\ngC7AuCjlu4nI1yKS7t0fxonIURHkuojI9945t1ZEHhJ3b9vrXljIGE1E5B1vuxwRWe/dj1qGyZ3p\njbHD02emiPT1rT9ORD4RkdW+c+75SNdPIXoUeU35GAvUAk6PZr9GdJjnKL6sBo4WkUNVdUExsjcC\n84ERQB5wLvCqiIiqvubJ3A68DOwEHgME2OStKyxuHWm5AocAHwFv4EJ+i2OkR2j/exCRfsDbwAzg\nfqAJcAdwjIh0U9Udvu0qAWOA6ThX8WnAncAyT9eiUKAdzj3/OvAucA3wqYj0VtXxYfKvApuBR4Ga\nnq6NPT3zgReBLcCZwFsiUltVX/Tk/ga8gHtqGwJUw93kj+KP8MAbwIXAS8BCoCFwHNAB+CXSsQr7\nXyLxT2AX8AxQFcgVkerAJGB/7/9eCxwDPAE0xR2/0iC4/28F7nvrDlyP+64f2CMk8hAwCPd/vwns\nhwsbfx/6flV1kYj8E3haRD5T1VEiUgP3HaUB/sTbfrhz6zkgAzjF239t4D6fnAINgK+B/3njXww8\nKSJzVXWMp58AX3rH5A1gES7UPQBoi/uOjGA4Bvc9/lycoIh0xJ3n6bj8pDzgBmCiiJygqj96cs2A\nCbhr+N9AFu68zSW63J7/4a7RF3H38MY4w6MlsMYbox/wFu5e+TguQtAN6A0M9/ZzCVAdd5/5HTgS\nuBVoDlxWzP9a7DXlE08DsoFjcfdtIxaoqr3i9ML9sOcCu4EpuAv6dKByBNmqEZZ9DSwNWzYP+C6C\n7CNAfoTlV+NuEi19y1Z6y06Lgx4nevs+wftcGdiIMwaq+OTOAgqAR3zL3vG2fTBsn7OAmVEc79D/\ndZ5vWW3gV+CnsGNSAEwEJGwfQ4F1QL2w5R8BW0PHB/gCmFuMPtuAF6PQ+e0Iyyf4j693XAuApf7j\n6K37B7ADaBO2/HHv/GteinP3EW+8/4Qt/xzY7Pvc0ju/7wuT6+iNfb9vmeB+3NbjjJqXcYZetyjO\nwddwBlNq2DHKB67wLUv19v+Jb9mVno49w/bZ39v+6JIeH3vF5oUzAPKBGhHWFQAP+z5/gTMCDvQt\na4ozlib4lr2IM5w6+5bVwz3o7HUvjDBmXW/cO4uQqeONOSX8WgyTi3Qe3+fpdoBv2V737pJcU751\ni4BRQX+f5ellYbU4oqrjgJ44a74LcA/OK/KriJwbJrsr9LeI1BGRhrgfkjbiqjlizUpPv3CdY63H\n4bgnr1fVl1OgqqNxF/TZEbYJ9xD9gAtRRsN6Vd3z9KQupDkM6OZ5hfasAt5U787i40KclyFFRBqG\nXsC3uBtsd09uO3CAiBxehC7bgaNEZP8odY+Gd/XPuRkX445RepjO43HGaWnbKiiRv4uGIlLL+3wR\nzuj5NGzszThD7uQ9O3PHuh8uBPA1zkv5uKrO3mvQvc/BWt7+JgM1cB5PPxmq+pFv293ATPY+Xy7G\nee6WhOk4wdP9ZIygaAjkaTG5Ml5I9XTgC1VdHVquqhtxDy7H+c7J3sA0VZ3nk9sOfBiFPtk4A+Qk\nEalXiMzpuHP4yQjX4h7CzuMa3jk3Decd71aEDlFfUz62AQmVx5rsWFgtzqjqLOBiEakMdAUuwLnz\nPxWRw1R1EYCIHIsL7xyN+xHYswvc08yf8pb2kZWRFsZBjwO9bZdEWLcI5wr2k6Oqv4ct2wbUj3K8\nZRGWhcZuhbvBhFjlFxKR/XAGUH+cuz4cxRl64EpoTwVmisgynPH0kapO9cnfiwsbrRWRWcBoYJiq\nRjz2UbIqwrK2uDDRb8XoXBrWhH3e5r3Xx4W8Dsbd7CMdd8X90PyxQHWFiDyKCwvOw4Vl98ILn/wb\n9yNQJ2x/dcPE10UYdxvueIRoizOq4nF8jLJhP9z9KNJ9ZCHuHGzh/X0gMDWCXKRzdC9UNVdE7gOe\nBTaJyHRgFO66DaUOHOS9F5kqISItgH/hUhP8969I57GfEl1ToeGw1hsxxYyjMkJV83DhoVkishQX\nQroE+JeItMElJC7EGU5rcRfA2bjcnGg8fIVdGCmFLM8OXxAjPfaV/DIYI0T4MQj9fx8A7xWyzVwA\ndTk07YFzgDNwHqebReRRVX3Uk/lURCbhDOJewN3AfSJygXr5MBT9veVFoXNI77E4gy1SAnakH5Ro\nKez7CI1TCReGOMN7DycjwrLeuP+7Gc5zsMdgFZG6OE/ldly4cAWQA/TAhaXDz8Hi9AvpOA93Tkc6\nPmsjLDPKht+ByiJSU1Uzg1YGQFVfEJGRwPm4c3UQ8ICInKyqc6LZh+fpGod72HoCl9OZics3eo+i\n76Wluabqs2/XuRGGGUfB8JP3Hgq39AGqAOeq6q8hIRE5NcK2hf2YbvO2qaN7J+u1KoFe58ZAj3BW\n436Q2uNyfPy099bHkoMjLGvvva8qZtvfcJ6xFFX9rhhZVDUb+BTnBayMy4l4SESeCLnbvafN14HX\nxbVvmI1r8RAyjrbhbqDhHAgsL04Hj+VALVWdEKV8LFmO+35XqWqxT+biKh9PxR2DB3FhO3814Um4\nG/15qjrFt91BlJ7lQJeAjo9RNIu899a45ObC+A2XWN0+wroOOCMiZOSuJvJ9oG20Snne3cHAYO/c\nm4MrELmKP875TjjjPRKdvfH+qqp7wnkicloUw5f0mkrBec0sGTuGWM5RHBGRkwpZFcqzCd0YQh6C\nPd+H9wTdL8K2mUT+MQ1dUHvyS0SkJu5ijpbQU/i+6BHOTzjPwI0ikurb75m4m9qoEugXDc3EV7ov\nInWAvwKzVXVz4ZuBqhbgEo4vEpFDw9eLrzeVhE2r4XkGF+K+g1QRqeSN7ZfZgksWrupbvBxX0bjn\nQUVEzsHd7KLlE6CniPSKoHNd7+YZL/6Hl1gfaaX/OIlIa+Bp4DNVfRLnSTtPRK70bZKPO4b+c7AK\nruloafkElx/2twj6VfOq5oxgmIb7vovK3Qtdm9/izhd/W5ImQF/gB1UNeVTG4K6HLj65BsAVxSkj\nItXF9Q/zsxL30BRa/q33+YEIsiH+dC/1uIPiHyyjvqY8OuKqZadEEDdKiXmO4stL3o33C5whVAWX\nY3Mp7onjXU/uW1x1wigReQNXYRUqmW4ats9ZOEPjIVxMerP3RPwtLj/kbRF5BndxXYMzTKL9oY2F\nHuALXahqnhfDfxuYJCLDvX3d5h2DIVHqFi1LgKEicoSn93W4nJLwbt6F9f+5H+e9mCEib+LKZBvg\nwjqn8EfS47cishF3Q9qEu0HdgqsYyfSMynUi8hnuqTMDl8h5OHuX1g/FJQyPEZFPcPkMVxJFfoSP\nZ3Dex1Hi+gjNwrUm6IIL97XCVdrhrb8KaKWq4flEJcbLIfoH8Lhn/Pwf7oejDS4s8QbwvCf+Nu7p\n/2Zv2/+IyEXACyIyzkuunYrzpg0TkRe97a5k3/Ip3sddc6+JyMm47ywFZ5xfggt5FltKbsQeVV0p\nIvNxlb3vFiP+D09uioi8ijNA+uPuq/f65J7GnTPjROQl3IPc9TiPUn2KPpfaAeO9azEN9+B6Ie4e\nMtzTeaeIDMCV2P8oIh/hztmuuNkBrsHd75cDz4nrTbQDl2hd7ANlCa8pcOdvJlH2ijKiJOhyufL8\nwp20b+IS99Jx+SKLce7aRmGyZ+NCLpm4i+ounMcmvAy/MTASl5ORz97l3ofhflyycU87txG5lH8F\nMKIQnfdJD8JK+X3yF+O8SFk4F/l7wP5hMu8A6RF0egRX0VLc8V7p6XQarnVAlnfsLwiTCx2T7oXs\npxGuHHgVLt/lV5zheK1P5npctdNmb5wluNyCWt76VFyOzM/eMdrh/d0/wnh34AzbLOB7XCXLBGC8\nTyZ0XC8sROcauOTmxd73vwlXWXYHLkwYkvsUZ6jVKeZYPuKN16CQY9cybPn5nu47vNcCXB+og731\ntxLWZsFbfgDuh+VL37KjcQZMBi5U8rj3ne51XnnHaE4E3d8BloctS8F5quZ6x3kLrqrtodB3Zq9g\nXt45mk5Y6bv3ff8zbFlXXGFDOs5gGAscGWGfXXBh/Czv2nrAdw7uV4QuDbxrf4F3Hm/F3VP/dN3h\n7pU/eOfpNpwX7FLf+vY4L1a6dz2+hgvF5QNX+eQi3t+Ku6Z8ctNwVayBf5fl6SXewTWMpEdc5+55\nqtonaF0SFc/b9a6q3h+0LoYBe0Lfy4F7VfWdOI4zBPgbzhguFz984rrt/4TrFTavOHkjeiznyDAq\nCF6JfDVc2MEwEgJ1BSTP4PrAxQQJm6LD6xN0JS43qVwYRh73AZ+aYRR7AjeORORGEZkjbm6adBGZ\nKiJnhMkMEje/TZaIjBWRg8PWVxWRV0Rki7j5dj6TvRv+ISL1ReRDb4xtIjLUS1j2y7QQka/EzdG0\nUUSelj/P59RFRCaJm7drtUSex+kkEZklbj6dJSISnu9iGGWOqqapaj1V3Rq0LobhR1WfVtWOMdzl\nNBEZLCL9ReRhXB5ebVzfoXKDqvZV1b7FSxolJXDjCJdTcB+u83AP4DtghIh0APCSef+OS7w7EpcL\nM0b2njByCC7+exGuWqsZrurIz0e4BMxTPdkT8HX/9Yyg0bgk9aNxuRX9cD0uQjK1cTHklZ6+9+Bm\nQ77eJ9MKV4E1HhcffwGXIGyTAsYfxRqhGYYBX+HmRHwed59eBZyhvvYQhlEUCZlzJCK/A3er6jsi\nsh54RlUHe+vq4JLbrlbVT7zPvwGXq+oXnkx7XFn10ao60zO0FgA91JuqQER64y6gA1R1o1daPhKX\nJLzFk7kBl1S7n7qqq5twTx5N1ZVuI2428fNCTz0i8hRwpqr6y0iHA3VV9aw4HjbDMAzDMGJAIniO\n9uD1hrkcV3kz1StjbIrzwgB74tMzcHOWgSuNrhwmsxhXoRCSORrYpnvP4TQO52U4yiczL2QYeYzB\ntXk/1CczKWQY+WTae6XbIZnwksoxPl0MwzAMw0hgEqLPkYh0wpUjVsOVZ16gqotFpCfOgNkUtom/\n704TIFf37godLtOUvefUQlXzRWRrmEykcULr5njv4R1R/TLpReynjohUVd9khH68hMHe/FE+bhhG\n6aiG6+00Rv88T1/CYNe8YcSUmF73CWEc4RpmdcV5aS7GNYAr7UziQVBYQ8GS0JvoZo02DCM6/oLL\nNUxU7Jo3jNgTk+s+IYwjL0wV8sjMFpEjgdtxJceC8w75vTFNcI0KATYCVeTPc4o18daFZMKr11Jw\nDb/8MkeEqdbEty703iSCjEYhs6Mwr5HHKoBu3bpRu3btvVb07t2bM844I9I2pWLAgAEMHjw4ZvsL\nagwbx8b55ptvGDNmDPPnz6dTp05kZcFPP+3Euz2sitvAsWEVwAcffECHDh3iOlB5+b5tnIo9Tk4O\njBgBw4bBxo1Qs/kt7H74OppXq8JNeXnce911EKPrPiGMowhUwnVLXek1rTsVbzZ0LwH7KOAVT3YW\nrsX7qbhpOkIJ2S1xoTq893oi0s2Xd3QqzvCa4ZN5UEQa+fKOeuFCZWk+mcdEJEVV830yi1U13Sdz\nZtj/08unS2HkAAwdOpTu3bsXI7pv1K1bt1yMYePYON27d+fBBx+kT58+jBw5kmHD4Oqrf8YVviZ8\nqCoHoEOHDnY92jg2ThHs2AGvvgqDB8OWLXDJ5QWk/n0FH9xdib7HHcN/2rVjydy5IfGYXPeBJ2SL\nyOMicryIHCginbzqrxOBDzyRIcA/RORcEekMDAPW4c1A7HmL3gKe9/oL9cDN4TRFVWd6MotwSdFv\nisgRInIs8BIwXN18TuCmh0gD3vd6GfXGVaa9rKq7PZmPgFzc/GUdReQy3BQdz/n+pdeBNiLylIi0\nF5GbcaFC/1w4hmHEgRkzoFWroLUwDCMWbN0KAwfCgQfCI4/A+efDtIW5bLx7Dv/N/ZVONWrwYYcO\n1Kocez9PIniOGuPNs4Xz0swFeqnqd+Cag4mbvPUN3KR9P+BK5XN9+xiAm6/mM9zMyd/gJgH1cwXw\nMq6SrMCTvT20UlULxM2G/hpuLp1M3ESIj/hkdoib+fwVXMv2LcBAVX3LJ7NKRM7GzZ92G86Qu05V\nbVJAw4gzM2ZAp06walXQmhiGUVq2bHFeopdegrw8uOEGuPtu2FB7BxcsWMDuggK+69qVZ6pXRyQW\nKb9/JnDjSFWvj0JmIDCwiPW7cJMK3lqEzHZc+/iixlkLnFOMzHycZ6somUl4fn3DMMqG7GyYMweq\nXPKra8NqGEZS8dtv8Nxz8PLL7vMtt8Bdd0HjxjBs40b6z17MYbVq8XmnTjSvWpVn4qhL4MaRUfb0\n7Rv/bvNlMYaNY+P4x5k92z1lbmoY3tXDKI/ft41TfsbZsgWefdYZRSJw221w553QqBHkFRRw17IV\nPL9uHdc0bcqrbdtSLSWlVOOUhITskF0REZHuwKxZs2aVScKcYZQ3hgyBBx6Aai+/zfbrrwPXEf/n\noPUqDLvmjYrOtm3OU/TCC6C6t1EEsH33bi5PS2Pstm0MPvhgbm3evNAw2s8//0yPHj0gRte9eY4M\nwygXzJgBXXrmMTM/r3hhwzACY+dOZxA9+yzk5sKtt8I99/xhFAEsy8ri3Pnz2ZibyzddunB6gwZl\nqqMZR4ZhlAtmzIBjr8phZtCKGIYRkZwceO01ePxxV55/443O29u06d5yk7Zv54L582mYmsqM7t1p\nV6NGmesaeCm/YRjGvrJpE6xcCU0Oyw5aFcMwwsjLg7ffhnbtnIfovPNg6VLnPQo3jIZt3Mhpc+bQ\ntVYtpgdkGIEZR4ZhlANmeK1cq7fJoVolu60ZRiKg6jpad+kC110HRx8NaWkwdCi0bBkuqzy8ciVX\nL1rEVU2aMKZLFxqkpgajOBZWMwyjHDB9Ouy/P2ytlk3zqlVZHrRCJWDAgAHUrVuXvn37llk1kWHE\nm6lTnZdo6lQ49VR4/33oUUiDm9yCAq5bvJgPNm3iidatua9ly6j7Fw0fPpzhw4eTnp5evHAJMOPI\nMIykZ/p091S6MieH5lWqJJVxNHjwYKtWM8oNy5bB/ffD55/DYYfBmDHQq1fh8tt37+bCBQuYkp7O\nfzt25LLGjQsXjkDoocJXrRYTzP9sGEZSk58PM2c642hFdjYHVKsWtEqGUeHYuhXuuAM6dnTX47Bh\nMGtW0YbRupwcjv/lF37JyGBc164lNoziiRlHhmEkNQsWQGYmHHmU7vEcGYZRNuze7RKrDz7YJV0/\n+igsXgx//SsUlf63IDOTnrNnsyMvjyndunF8vXplp3QUmHFkGEZSM306pKRAsy67yFWledWqQatk\nGBWC0aOhc2fXuPGSS1wF2gMPQPXqRW83eft2jps9mwaVKzOte3c61KxZNgqXADOODMNIaqZPd9Uw\nGyUHgAPMODKMuLJ4MZx1Fpx9NjRvDrNnwxtvQJMmxW/75ZYtnD53Ll1r1mRSt240S9Dr1YwjwzCS\nmlAy9ops1+NofwurGUZc2LED7r4bOnWCRYvgiy9g3Dj3cBIN72/cyAXz53NmgwZ806ULdSsnbk2Y\nGUeGYSQt27bBwoVw1FGwIieHZlWq7JmU0jCM2KDqSvHbtXMdrh991PUrOv98N1FsNLy4bh1XLVpE\nv6ZN+aRjx4S/Ts04MgwjaQk1f+zZ03mO2hSX7GAYRomYNw9OOAGuugpOPNF5jB58EKItClVV/rVq\nFbcvW8bdLVrwZvv2VE6CRq2Jr6FhGEYhTJsGDRtC27bOc9TGyvgNIybs3Al33QXdusFvv7nw2ccf\nQ4sW0e9DVbl3xQoeXrWKf7duzdNt2kTd3DFozDgyDCNpmTbNeY1EYHl2Nq3NODKMfULV5RJ17OhC\naIMGwdy5rst1SShQ5e9Ll/Ls2rW8ePDBPHjggUljGIEZR4ZhJCn5+S4Zu2dP2JmXx+bduznYwmqG\nUWrWrHGTwl54oetunZbmQmglrXHIV+X6xYt5bf16hrZvz60HHBAfheOIGUeGYSQlaWnO9d+zpwup\nARxkxpFhlJj8fNfIsWNH19X6889h5Eho1ark+8orKKDfokW8t3Ej73fowHX77x9zfcsCM44Mw0hK\npk1zHXiPOMKF1CAxjCMRaSUi34nIAhGZIyLBK2UYhTBvHhxzDAwYAP36uYeOCy+MvgrNT15BAX9d\ntIjhmzYxvGNH/hJN46MExYwjwzCSkmnTXH+VWrVgWXY2tVJS2C81NWi1AN4F/qGqhwInAruCVccw\n/syuXfDww9C9O2RkwOTJ8PLLULdu6faXV1DAlQsX8tlvv/HxoYdyaQLNk1YaErcDk2EYRhFMmwan\nnOL+Xp6dzcHVqwee8CkiHYFcVZ0KoKrbA1XIMCIwcyZce63rdP3gg+61L42qQ4bR51u28HHHjly4\n336xUzYgzHNkGEbSsXWru7H37Ok+L8/O5qDEqFRrC2SKyEgR+UlEHghaIcMIkZMD99/vrptq1Vx+\n0aOP7pthlK/K1YsWlSvDCMw4MgwjCZk2zb0fc4x7X56Ts8/5RiJyvGfU/CoiBSLSJ4LMLSKyUkSy\nRWS6iBwRJlIZOA64ETgGOF1ESlgEbRix56efoEcPGDwYHnvsjzkJ94V8Va5dtIiPN29meIcO5cYw\nAjOODMNIQqZOdZNctmkDuQUFrImBcQTUBH4BbgY0fKWIXAY8BzwCdAPmAGNEpJFP7FfgJ1Vdr6q5\nwGjgsH1VzDBKS26uyy06+mioXt15ix54APZ1WjNV5aYlS/hg0ybe79CBi5M8xygcM44Mw0g6pk51\nXiMRWJ2TQwHsc48jVf1GVR9W1RFApOSlAcAbqjpMVRfhvENZwLU+mR+BxiJSV0QqAScAC/dJMcMo\nJQsWOKPoiSecgTRtmps0dl9RVQYsW8abGzbwVvv29E3iqrTCMOPIMIykYvdul1AaCqktK4MyfhFJ\nBXoA40PLVFWBcUBP37J84EHgB5wXaomqjo6bYoYRgYICGDLEhdFyclwI7eGHIVbFnA+vWsULv/7K\nq23b0i9J+xgVh1WrGYaRVMyZA1lZcOyx7vPy7GxSRThgX7JKi6cRkAJsClu+CWjvX6CqY4Ax0e54\nwIAB1A2rn+7bty99+/YtnaZGhWb9erj6ajcX2u23O69RLJ8bnl2zhsdWr+aZNm24qXnz2O24BAwf\nPpzhw4fvtSw9PT2mY5hxZBhGUjF1qpvOoHt393l5Tg6tq1UjJYnmbfIzePBguof+GcPYB0aMcCX6\nVavCmDHQq1ds9//m+vXcs2IF/zjwQO5u2TK2Oy8BkR4efv75Z3r06BGzMSysZhhGUjF1Khx++B/l\nx6EeR3FmC5APhCdXNAE2xntwwyiK7Gy4+WY4/3w44QQ3UWysDaPPNm/mxiVLuKVZMwaVZl6RJMOM\nI8MwkoqpU/8IqYHX4yjOxpGq7gZmAXvK8sV1nDwVmLov+x4wYAB9+vT5U5jAMKJhwQI48kh45x14\n7TX43/+gUaPitysJ47dt4y8LF3JZ48a82LZt4M1W/QwfPpw+ffowYMCAmO7XwmqGYSQNa9e6VygZ\nu0CVFTk59I+BcSQiNYGD+aNSrY2IdAW2qupa4HngXRGZBczEVa/VwE0XUmosrGaUBlV4+2249VbX\n0uLHH2NTiRbOzzt3cv78+Zxcrx7vHnIIlRLIMII/QmyxDquZcWQYRtIwZYp7D3XG/nXXLnIKCmgb\nG8/R4cAEXI8jxfU0AngPuFZVP/F6Gg3ChdN+AXqr6m+xGNwwoiUjA266CT74AK6/Hl54AWrUiP04\ny7OzOXPuXDrUqMFnhx5KlUoVJ9hkxpFhGEnDlCnQtq1rAAmw1Cvjj0XOkap+TzGpBqr6KvDqPg/m\nI1StZhVqRjQsWAAXX+w8qB9+CFdcEZ9xfsvN5Yy5c6lbuTJfde5MrX3tGhknQpVrVq1mGEaFZfJk\nOO64Pz4vy84mBWiVGPOqlQoLqxnR8sEHcMMNLow2axa0b1/8NqUhKz+fc+bNY0deHtO6d2e/KlXi\nM1AMiFfiWvIlAAAgAElEQVRYreL4yAzDSGp27HBVOP5k7KXZ2RxYrVqFcvcbFY9du1w12l//6rxG\nM2bEzzDKV6VvWhoLMjP5qnNn2sS/EjQhMc+RYRhJwfTprvNvuOcoRvlGhpGQrFvnDKLZs+H116F/\nfzdtTry4c9kyRv3+O1927szhderEb6AEx4wjwzCSgilTXIlyu3Z/LFualcXJ9esHp1QMsJwjozB+\n+MEZRlWquL+PPDK+4724bh0v/vorr7Vty1kNG8Z3sBhhOUeGYVRoJk92IbXQU3OBKstzcvhbknuO\nLOfICEfV9Sy6/XZ3zn/yCcR70vuvfv+dAcuWcXeLFtwY0LQgpcFyjgzDqLDs3u3Cav58oxiX8RtG\nQrBrlwud3XKLe40dG3/DaG5GBpenpdGnUSOeatMmvoMlCYEbRyLygIjMFJEdIrJJRL4QkXZhMu+I\nSEHYa3SYTFUReUVEtojIThH5TEQah8nUF5EPRSRdRLaJyFCv8ZtfpoWIfCUimSKyUUSeFpFKYTJd\nRGSSiGSLyGoRuSfC/3WSiMwSkRwRWSIiV+/70TKMiklosll/vlEsy/gNIxHYvBlOOw2GDXMdr4cM\ngdTU+I65KTeXc+fNo2316nzQoUPCNXkMisCNI+B44CXgKOA0IBX4VkTC73hf4xqvNfVe4cH5IcDZ\nwEXACUAz4PMwmY+ADriW/2d7cm+EVnpG0GhcuPFo4GqgH67pW0imNm7G7ZVAd+AeYKCIXO+TaQWM\nAsYDXYEXgKEicnoUx8MwjDAmT3ZzqfmjT6Ey/tZJXMZvGCHmzXM5RUuXwsSJ0K9f/MfcVVDABfPn\ns1uVkZ06UTMlJf6DJgmB5xyp6ln+zyLSD9gM9AAm+1btKqwTrYjUAa4FLvcauSEi1wALReRIVZ0p\nIh2A3kAPVZ3tydwKfCUid6vqRm/9IcDJqroFmCci/wSeFJGBqpoHXIkz4K7zPi8UkW7AncBQT6Wb\ngBWqeq/3ebGIHIebbmBsKQ+VYVRYfvgBjjrqj8lmwXmOWlWrRmqSl/FbQrYxahT07QsHHwwjR0KL\nFvEfU1Xpv3gxP+/cyaRu3TggSR8y4pWQnYh3lXq41v1bw5af5IXdFonIqyLSwLeuB87QGx9aoKqL\ngTWAN9EARwPbQoaRxzhvrKN8MvM8wyjEGKAucKhPZpJnGPll2otIXZ/MuDD9x/h0MQwjSlSdceQP\nqYHzHJWHkNrgwYMZOXKkGUYVEFU39UefPi6c9sMPZWMYAQxet45hmzbx9iGHcGQSl+z37duXkSNH\nMnjw4JjuN6GMI2+W6yHAZFVN8636GrgKOAW4FzgRGC1/TA3cFMhV1R1hu9zkrQvJbPavVNV8nBHm\nl9kUYR/ESKaOiFTFMIyoWboUfvsNjj8+bHlWFm3jMaGUYZQBeXlw221wxx1w113w+edQq1bZjP3t\n1q3cs3w597ZowRWhuXiMvQg8rBbGq0BH4Fj/QlX9xPdxgYjMA5YDJ+EmikwELIvNMOLADz9ApUpw\nzDF/LCsvZfxGxSQzEy6/HL7+2jV2vOGGsht7eXY2l6el0atBAx63yrRCSRjjSEReBs4CjlfVDUXJ\nqupKEdkCHIwzjjYCVUSkTpj3qIm3Du89vHotBWgQJnNE2HBNfOtC7+GmdhNceK44mR2ququo/y2U\nf+DHchGMiswPP0CXLuD3/IfK+H8dPZo+o0btJR/r3APDiCWbNsG558LChS7X6Iwzym7szPx8Lpg/\nn4apqXzUoQMpVplWKAlhHHmG0XnAiaq6Jgr5A4CGQMiImgXk4arQvvBk2gMtgWmezDSgnoh08+Ud\nnYrz+MzwyTwoIo18eUe9gHQgzSfzmIikeGG5kMxiVU33yZwZpnYvny6FYg3hDGNvJk+Gs87ae9kS\nr4y//5VXcnD//nuti3UzuHhjCdkVh2XLoHdv15Zi0iTo1q3sxlZV/rZ4MSuys5nevTv1490joIwo\ntx2yReRVXFl+HyBTREIel3RVzfH6ED2CK8vfiPMWPQUswSU5o6o7ROQt4HkR2QbsBF4EpqjqTE9m\nkYiMAd4UkZuAKrgWAsO9SjWAb3FG0Psich+wP/Av4GVV3e3JfAQ8DLwtIk8BnYHbgNt9/9brwC3e\n+rdxRtjFOM+YYRhRsmEDLF/+53yjJVlZVBahVZJW2PixB6KKwU8/OSO/QQOYNg1atSrb8V9Yt47h\nmzfzcceOdCqr5KYyoDx3yL4RqANMBNb7Xpd66/OBLsAIYDHwJvAjcILPYAFXJj8K+My3r4vCxroC\nWISrJBsFTAL2RHtVtQA4xxtzKjAMeBdnnIVkduC8QK2An4BngIGq+pZPZhWuj9JpwC+ebtepangF\nm2EYRfDDD+49vFJtSXY2B1WrRuUkL+M3KgZjx8JJJ8FBB7k5AsvaMJq8fTv3rFjBnQccwKXxbrdd\nTgjcc6SqRd7dVDUHKDYq6+Xy3Oq9CpPZjutTVNR+1uIMpKJk5uMq5oqSmYRrMWAYRin54Qf3g7L/\n/nsvX5KVRTurVDOSgE8+gSuvdKX6n34KNWsWv00s2ZSby6VpafSsU4cnLQE7auyxyzCMhGXSJDgx\nwmPIkuxs2lmlmpHgvPmmq0q79FIYMaLsDaN8Va5IS6NAlY87dkz6hqlliR0pwzASkq1b3ZQKJ5yw\n9/LcggJWZmeb58hIaJ599o8JZIcNi/8caZEYtGoVE7dv578dO7J/VWuxVxLMODIMIyGZMsV1EA43\njlbm5JAP5jkyEhJVGDgQ7rkHHnoIXnzR9ekqa77dupV/rV7NoNatOal+/bJXIMkJPOfIMAwjEpMm\nwQEH/Dl5dUlWFkC58RxZKX/5QRXuuw+eeQaefNL9HQQbdu3iyoUL6VW/Pg+0bBmMEmVEuS3lNwzD\niMSkSc5rFN6nbkl2NjUrVWL/KlWCUSzGWCl/+UAVbr8dXnrJeYtuLbQ0KL7kq3LlwoVUFmFYhw5U\nKueNHstzKb9hGMZeZGTArFl/DqnBH5VqUs5v+kbyUFAAN9/sDKM33gjOMAJ4YvVqJmzfzgcdOtC4\nnDxABIF5jgzDSDimTYP8/MjG0eKsLMs3MhKGkGH0n//AW2/BtdcGp8vU9HQGrlrFQwceyCmWZ7RP\nmOfIMIyEY9IkaNQIDjnkz+uWWKWakSAUFLhqtEQwjLbv3s0VaWkcVacOjxx4YHCKlBPMc2QYRsLx\n/feR84125uWxITfXPEdG4KjCbbe5MNpbb8E11wSpi3LT0qVsy8tj4mGHWef4GGBH0DCMhCI7G2bM\ncNMthLPUm3DWPEdGkKjCnXfCK684r1GQhhHAB5s28d/Nm3mjXTta2YNDTDDjyDCMhGL6dMjNjdwZ\ne3GojN9+AIwAeeghGDLEGUfXXx+sLquys7ll6VL+2qQJlzdpUvwGRlSYcWQYRkIxcaKbubxTpz+v\nW5SVRZPUVOoF0W7YMIDHH4cnnoDnnnOJ2EGSr8pVixbRoHJlXmrbNlhlyhmWc2QYRkIRyjeKlDax\nOCuLQ8pZSM2aQCYPL7/svEaPPurCakHz3Nq1TE5PZ+Jhh1G3csX8ObcmkIZhlHtyclxY7amnIq9f\nlJXFUXXqlK1SccaaQCYHH3zg+hfddRf8859BawNzMzL458qV3NOiBSfUqxe0OoFhTSANwyj3zJgB\nu3ZFzjcqUGVJdna58xwZic+XX0K/fnDddW5qkKD7j+YWFHDVwoW0q1GDQa1bB6tMOcU8R4ZhJAwT\nJ0L9+tCly5/Xrd21i+yCAtqbcWSUIVOmwKWXwnnnweuvB28YAQxatYq0rCxmdu9OVSvbjwt2VA3D\nSBi+/x6OP77wfCPAPEdGmbFgAZxzDhx1FHz4ISRCWs9PO3bw5Jo1/PPAAzmsdu2g1Sm3mHFkGEZC\nkJMDU6dG7m8ELt+oqggHVqtWpnoZFZN16+CMM+DAA2HECEiE025XQQFXL1pE11q1uL9ly6DVKdck\ngB1sGIbhErF37YKTT468flFWFm1r1CAlEeIaRrkmPR3OPNN5MEePhrp1g9bI8eiqVSzNzmZWjx6k\nWjgtrphxZBhGQhDqbxQp3wjKZxm/kXjk5sKFFzrP0ZQp0KxZ0Bo5ft65k6fXrOGRVq3oXKtW0OqU\ne8z0NAwjIZgwwVWpFfZAvMiMIyPOqEL//jB5sguldewYtEaO3QUFXLd4MZ1q1rRwWhlhxpFhGIGT\nne3CaoXlG+3My2N9bi7tbdoQI4489hi89x68/bZrRJooPL12LfMyMnj7kEMsnFZG2FE2DCNwpk51\n4YzC8o2sUs2IN//9Lzz8sOt+/Ze/BK3NHyzOymLQqlXc3aIF3a06rcww48gwjMCZMAEaNYJDD428\nflFowlkzjow4MGOGa/J45ZWJ0f06RIEq/RcvpmW1ajzSqlXQ6lQoLCHbMIzAmTjRhdQKixgszMqi\neZUq1EmERjMxxuZWC5a1a12Dxx494M03E6PJY4i3NmxgUno647t2pXpKStDqJCQ2t5phGOWSjAz3\n5P7CC4XLLMzKokPNmmWnVBlic6sFR1YWnH8+VK0KX3yRGL2MQmzKzeXeFSu4pmlTTqlfP2h1Ehab\nW80wjHLJ5MmQlwennlq4zMKsLDpYSM2IIapurrRFi2DkSGjcOGiN9ubOZctIAZ456KCgVamQmHFk\nGEagfPed6yXTrl3k9bkFBSzLzqajGUdGDHnmGZeE/d570LVr0NrszditW/lo82aePeggGqamBq1O\nhcSMI8MwAmX8eDjllMJzPZZlZ5OnWm7DakbZM3YsPPAAPPggXHxx0NrsTU5+PjcvXcqJdetyddOm\nQatTYTHjyDCMwNi6FWbPLj6kBlhYzYgJK1fC5ZdD794waFDQ2vyZZ9auZVVODq+2a4ckUnZ4BcOM\nI8MwAuP7713uR2H9jQAWZmbSoHJl9rPwgrGP5OTARRdBvXrw4YeQaAVgK7OzeXzNGu484AA6mqc0\nUKxazTCMwPjuOzjoIDfzeWGkZWXRsWZNe4o29plbb4WFC1039kQsALtt2TIapabyz6IuCKNMMOPI\nMIzAGD++6JAauLDaEdYZ2NhH3nkHhg5174mWgA0wassWRv3+O58deii1ymE/r2TDwmqGYQTChg3u\nKf6UUwqXKVBlsZXxG/vIvHlw881w/fWuE3aikZOfzx3LlnFa/fpc2KhR0OoYmOfIMIyAGD/evRdl\nHK3OySG7oMCMI6PU7NwJl1ziWkW8+GLQ2kTm+XXrWL1rF1927mzh4wTBjCPDMAJh3DgX3thvv8Jl\n0rxKNUtONUqDKtx0E/z6K8yaBdWrB63Rn1mXk8O/V6/m1ubNrV1FAmFhNcMwyhxVZxyddlrRcgsz\nM6lZqRItqlYtG8WMcsV777mqtDfeKLzJaNDcv2IFNVNSbGLZBMM8R4ZhlDmLF7un+eKMozRvTrVk\nCjWIyCpgO6DAVlUtJuXciAeLF8Mtt8A118AVVwStTWSmp6fz4ebNvNmuHXUtCTuhsG/DMIwyZ9w4\nSE2F448vWm5BZiaHJl++UQHQU1Wzg1akorJrl2v02KIFvPRS0NpEpkCVO5Yt47Batbhm//2DVscI\nw4wjwzDKnHHjoGdPKCrFQlVJy8ri4qKSkhITwVIWAuWhhyAtDWbMKPocC5LhmzczY+dOJh52GClJ\n5BmtKNgFbBhGmZKXBxMmFB9SW7NrFxn5+RyaqL9uhaPAJBGZISIJGtApv4wbB889B088AYcdFrQ2\nkcnOz+eBFSu4sFEjTqxXL2h1jAgEbhyJyAMiMlNEdojIJhH5QkT+lDonIoNEZL2IZInIWBE5OGx9\nVRF5RUS2iMhOEflMRBqHydQXkQ9FJF1EtonIUBGpGSbTQkS+EpFMEdkoIk+LSKUwmS4iMklEskVk\ntYjcE0Hfk0RklojkiMgSEbl6346UYZQPfvoJduwovvnjgsxMgDIzjkTkeBEZKSK/ikiBiPSJIHOL\niKz0rv3pInJEhF0dq6o9gPOAB0WkU9yVNwA3V9/VVzvD+447gtamcAavW8fG3FyeatMmaFWMQgjc\nOAKOB14CjgJOA1KBb0VkT9GliNwH/B3oDxwJZAJjRKSKbz9DgLOBi4ATgGbA52FjfQR0AE71ZE8A\n3vCNUwkYjQs3Hg1cDfQDBvlkagNjgJVAd+AeYKCIXO+TaQWMAsYDXYEXgKEicnqJjoxhlEPGjoU6\ndeDII4uWW5CZSa2UFFqWXaVaTeAX4Gac92cvROQy4DngEaAbMAd3H9qra5+qbvDeN+LuJ93jq7YR\n4uabITsb3n0XKiXCr1sENuXm8sSaNdzSvDkHJ18+XYUh8JwjVT3L/1lE+gGbgR7AZG/x7cC/VHWU\nJ3MVsAk4H/hEROoA1wKXq+r3nsw1wEIROVJVZ4pIB6A30ENVZ3sytwJficjd3o2sN3AIcLKqbgHm\nicg/gSdFZKCq5gFX4gy467zPC0WkG3AnMNTT9yZghare631eLCLHAQOAsTE6dIaRlHz7rWv8WFxx\nzoLMTDrWqFFmlWqq+g3wDYBEHnQA8IaqDvNkbsQ9ZF0LPO0tqwFUUtUMEakFnAJ8XAbqV3iGD4eP\nP3bvzZsHrU3hPDp6NKn16tn8aQlOItrW9fBKYAFEpDXQFOeFAUBVdwAzgJ7eosNxhp5fZjGwxidz\nNLAtZBh5jPPGOsonM88zjEKMAeoCh/pkJnmGkV+mvYjU9cmMC/u/xvh0MYwKyY4dbtLPXr2Kl12Q\nlZUw+UYikop7YPPfYxR3nfuv6ybAZBGZDUwF3lXVWWWpa0Xk11+d1+jyy90rUVm8dCn/qV2bh1as\noEFqatDqGEUQuOfIj/e0NgSYrKpp3uKmOANmU5j4Jm8duBtSrmc0FSbTFOeR2oOq5ovI1jCZSOOE\n1s3x3lcUIZNexH7qiEhVVd2FYVRAJk50CdmnFxNgLlAlLTOTyxs3Llqw7GgEpBD5um4f+qCqK4ES\npQEPGDCAunXr7rWsb9++9O3bt3SaVjBUoX9/1/36lVeC1qZoHpw8mea1anHLJZcErUpSM3z4cIYP\nH77XsvT09JiOkVDGEfAq0BE4NmhFSoHVYhpGMYwdC61bw0EHFS23OieHrIKCZOxxVGIGDx5M9+6W\nllRa3nsPRo+GL7+EBg2C1qZwpv34I/9r3ZphGzZQrXbtoNVJaiI9PPz888/06NEjZmMkjHEkIi8D\nZwHHhxIaPTbiDI8m7P3U1gSY7ZOpIiJ1wrxHTbx1IZnw6rUUoEGYTHj1SRPfutB7kwgyGoXMjuK8\nRvYUaZRnxo51XqPi0oiirVQriydIjy1APpGv641/FjfKgnXrXFXa1VfDOecErU3hqCr3LVxIVxH+\nkshxP2MPCWEceYbRecCJqrrGv05VV4rIRlyF2VxPvg4uTyjkRJ0F5HkyX3gy7YGWwDRPZhpQT0S6\n+fKOTsUZXjN8Mg+KSCNf3lEvXKgszSfzmIikqGq+T2axqqb7ZM4M+zd7+XQpFHuKNMora9a4KR0e\ne6x42QVZWdRJSeGAYirVyuIJEkBVd4vILNw9YyTsSQM4Fdinud5DD0T2EFQyVOHGG6FGDRg8OGht\niuabCRP4oWVLRu/cSSXLNYopoQekmD8UqWqgL1wobRuupL+J71XNJ3Mv8DtwLtAZ+D9gKVAlbD8r\ngZNwiZNTgB/CxhoN/ITzDh0LLAbe962vhMsr+hrogqte24SrlAvJ1AHWA+/hQoCXARm46rWQTCtg\nJ/AULh/hZiAXOK2I49Ad0FmzZqlhlEfefFO1UiXV338vXvavaWl6dCmvhVmzZinOk9tdS3Yvqolr\nvXEYbgqQO7zPLbz1lwJZwFW4qtY3vPvSfiUZR+2ajwkffaQKqiNGBK1J0eTn52vXDz/U4996Swvy\n84NWp9xS2uu+sFcieI5uxP1DE8OWXwMMA1DVp70S2Tdw1Ww/AGeqaq5PfgDO7f0ZUBVXkntL2D6v\nAF7GVZgUeLK3h1aqaoGInAO8hqs0yQTexfU1CcnsEJFeOK/VTzh3+0BVfcsns0pEzgYGA7cB63DG\nU3gFm2FUGMaMcb2NoskLmZ+ZSY+yz8s4HJiAux8prqcRuAeha1X1E6+n0SDcA9wvQG9V/a2sFa3o\n/PYb3HYbXHYZ9PlTq87E4pPRo5nTrBmTCwqQRG2+ZPyJwI0jVY3qbFHVgcDAItbvAm71XoXJbMf1\nKSpqnLVAkdFrVZ0PnFiMzCScB8swKjx5eW5ah9tvj0K2oIC0zEz6NW1avHAMUdcjrcj7kaq+ivNS\nGwEyYAAUFMCL+xTQjD95eXk8nJnJORs2cOzf/ha0OkYJCNw4Mgyj/PPjj7B9O/TuXbzssuxsdqnS\nOUF6HMUbyzkqGWPHwocfui7YidPpITLDvvySpU2a8EnLlkGrUm6JV86RGUeGYcSdMWOgXj04ItJM\nZGHM8yrVOlUQ48iKMKInK8slYZ9yClx1VdDaFE3url0Mys/n4iVLOKx//6DVKbeEHirKbSm/YRjl\nlzFj3GSgxU0ZAs44apKayn5VqhQvbFQoHnvMdcP+5pvi20EEzVsjR7KmYUNG2+SySYllhxmGEVe2\nbYOZM6MLqYFLxu5cq1Z8lTKSjrQ0eOYZePBBaNs2aG2KJicnh3+npHDF0qV0NK9gUmKeI8Mw4sq4\ncS55NlrjaF5mJuc2bBhfpRIIyzkqHlW45RZo1QruvbdY8cAZOnIkGxo25OHiWsEb+4zlHBmGkZR8\n8w107AgtWhQvm5mfz/Ls7AqTjA2WcxQNH33k5uUbMwaqVQtam6LJycnhidRU/rJsGe1OPTVodco9\n8co5srCaYRhxQ9UZR2eG94svhLTMTJSKk4xtFM+OHXD33XDxxdCrV9DaFM/QESPYWKcO/+jZM2hV\njH3AjCPDMOLG3Lmwfn30xtG8zEyE4udUMyoO//qXM5Cee6542aDJycnhiSpVnNeoS5eg1TH2AQur\nGYYRN77+GmrWhOOOi05+XmYmB1WvTo2UlPgqlkBYzlHhLFwIQ4bAwIGQDK2C3vnySzY0bMhDrVsH\nrUqFwXKODMNIOr7+2vWkKWb+2D3Mz8ysUPlGYDlHhaEKd9zhjKK77gpam+LZnZvLk5UqcdmyZbS3\nvkZlhuUcGYaRVKSnw5Qp0YfUVJU5GRkVzjgyIvPVV/Dtt/D884mfhA3wwahRrGnYkIfM0C0XmHFk\nGEZcGDcO8vOjN4425uby2+7ddLUeRxWe3FznLTr11MSfWBYgLz+fx/PyuGDJEjodfnjQ6hgxwMJq\nhmHEha+/hkMOcb1pomFORgaAGUcGr7wCy5bBZ58lfidsgM++/ppljRvzcbNmQatixAjzHBmGEXNU\nYfRoOOus6LeZk5lJrZQUWidDDMWIG7//DoMGwd/+Bp07B61N8WhBAU+mp9N72TK6R1t5YCQ85jky\nDCPm/PILbNgAZ58d/TZzMzLoUrMmlZLBVRBDrFptbx57zIVjBw0KWpPo+Pq775jTvDlDdu8OWpUK\niVWrGYaRNIweDbVrR1/CDy6sdkK9evFTKkGxarU/WL7chdQGDoTGjYPWJjqeWL+eowsKOPGqq4JW\npUJi1WqGYSQNX30Fp58OVapEJ5+Tn8+irCzLN6rgPPCAM4ruuCNoTaJjyowZTG7Zkgfq1UMq2c9p\necI8R4ZhxJQtW2D6dHjzzei3ScvKIh/oamX8FZaZM+HTT+Gdd6BGjaC1iY6n09LoULUq51x6adCq\nGDHGTF3DMGLKmDEuIbtEydgZGQg2p1pFRRXuvx86dYK//jVobaJj4aJFjGzdmruBSpXNz1DesG/U\nMIyYMno0dO8O++8f/TZzMjI4qHp1atmPTIVk7FiYMAFGjoRkmTnmualT2b9ePf5y3nlBq2LEAbsT\nGYYRM/LyXH+jv/+9ZNvNycyssCG1il6tVlDgco2OPRbOOSdobaJjw4YNvN+iBYPWrqVqBT1vEwWr\nVjMMI+GZOhW2bYNzz41+m9C0IQMOOCB+iiUwFb1a7fPP4eefYdKk5Gj4CPDyuHFUadiQG0rSq8KI\nC1atZhhGwjNqFDRpAiW5R63ZtYtteXl0s0q1Ckd+Pjz8MJxxBhx/fNDaREdmVhav16vH39asoV6T\nJkGrY8QJ8xwZhhEzvvzSNX4sSVXz7J07AehWu3actDISlY8+gkWL4P33g9YkeoZ9/TXb69XjtmSx\n5oxSYZ4jwzBiwrJl7oeuJCE1gNkZGeyXmkqzaJsiGeWC3btds8fzz4dkmau1oKCAwfn5XLhsGa0O\nPTRodYw4Yp4jwzBiwldfuaaPp51Wsu1mZ2TQrVYtJFkSToyY8N57sGIF/N//Ba1J9Hw1cSJLGzfm\nvYYNg1bFiDPmOTIMIyZ8+SWccgqUNHUoZBwZFYfdu+Hf/4ZLLkmOyWVDDFm3jqNWrqTnKacErYoR\nZ0rsORKRqsBRwIFADeA3YLaqroyxboZhJAnbt8P338OQISXbbktuLut27bJ8owrGsGGwapXra5Qs\nzF+4kO9atuSjjRuTp6zOKDVRG0cicixwO3AukAqkA9lAA6CqiKwA/gO8rqo746CrYRgJyjffuB5H\nffqUbLvZGRkA5jmqQIS8RhddlFxeo5dmzKBZ7dpcXNKkOiMpico4EpGRQHfgI6AX8JOqZvvWtwGO\nB/oCd4rIVao6Ng76GoaRgIwc6bpit2hRsu1mZ2RQKyWFg6tXj49iSUBFawL5wQewcmVy5Rpt3b6d\n9/ffnwfXrCHVmj4mFEE3gfwKuEhVd0daqaorgBXAeyLSESjBxAGGYSQzu3e7KUMGDCj5trMzMuha\nsyaVKnCYoiI1gczPhyeecBVqXboErU30DB0zhoL69bmhpNUGRtwJtAmkqr5RmGEUQTZNVcfvm1qG\nYSQLkyZBejqUZoqp2Tt3Wr5RBeLzz2HpUnjwwaA1iZ78ggJerVyZy5cvZ7/WrYNWxygjrFrNMIx9\nYuRIF07r2rVk22Xk5bEkO9vyjSoIqvD443D66XDEEUFrEz1fff89qxs25O+HHBK0KkYZUppqtRRg\nAPIR3Z0AACAASURBVHAp0BLYq3ObqjaIjWqGYSQ6qs446tOn5AU8v2RkoEB3M44qBKNHw5w5MHFi\n0JqUjFfWruXI/HwO79cvaFWMMqQ0nqNHgDuBj4G6wPPA/4ACYGDMNDMMI+GZM8eVZJcmpDYrI4Oq\nIhxqCa4VgieegGOOgRNOCFqT6FmyciXftmzJzampVr5fwShNh+y/AH9T1a9EZCAwXFWXi8hc4Gjg\nxVgqaBhG4vJ//wd168JJJ5V821k7d9K1Vi1SSzIRm5GUTJniXiNHJpeN8drkyTSsU4fLzjknaFWM\nMqY0d6WmwDzv7wyc9whgFHB2LJQyDCM5+OILOOccSE0t+bazdu6khyVjVwieeQY6dHCTEicLWTk5\nvNugAdetX0+1evWCVscoY0pjHK3jj1L95bi+RwBHALtioZRhGInPihUwdy5ccEHJt83Iy2NhVpYZ\nRxWARYtgxAi45x5IJifhJ+PHs71mTW448sigVTECoDSn6hfAqd7fLwH/EpGlwDDg7VgpZhhGYjNi\nBFStCr17l3zbUDL24WYclXuefRaaNYMrrghak5Lx+vbt9FqyhDYx7J1jJA8lzjlS1ft9f38sIquB\nY4ClqvplLJUzDCNx+eILV5ZdmmKzUDJ2xxo1Yq+YkTBs2gTvvw+DBjlDOln4ZfFiZjRvzv82bQpa\nFSMgSuw5EpGqIrLndqiq01X1+X0xjETkeBEZKSK/ikiBiPQJW/+Ot9z/Gh1Br1dEZIuI7BSRz0Sk\ncZhMfRH5UETSRWSbiAwVkZphMi1E5CsRyRSRjSLytIhUCpPpIiKTRCRbRFaLyD0R/qeTRGSWiOSI\nyBIRubq0x8cwEo3Nm12CbWlCagA/WTJ2heDVV6FyZejfP2hNSsYbM2fS7PffOSeZkqSMmBL1nUlE\n9vMMkgwgXUSmi8jBMdKjJvALcDOghch8DTTBJYQ3xc3j5mcILiH8IuAEoBnweZjMR0AHXFjwbE/u\njdBKzwgajfOoHQ1cDfQDBvlkagNjgJW4+ebuAQaKyPU+mVa4BPXxQFfgBWCoiJxe5FEwjCRhxAj3\nXto5OGft3GkhNY8BAwbQp08fhg8fHrQqMSU7G157Da69FurXD1qb6MnIyeGDRo24btMmUq0HV8Iz\nfPhw+vTpw4DSzF9UFKoa1QuXT7QBuB/XBHIRMCHa7UswTgHQJ2zZO8D/itimDi4Z/ALfsvbevo70\nPnfwPnfzyfQG8oCm3uczgd1AI5/MDcA2oLL3+SZgS+izt+wJIM33+SlgbpiOw4HRRfwP3QGdNWuW\nGkaic8YZqiefXLptd+7erTJhgr61fn1slfKYNWuW4h6yumuM70+xfJX3a/7NN1VFVJcuDVqTkjH0\nq69Uxo/XVeX0eymvxPq6L4lP+3Sgn6o+qaqDgXOB40WkrCLJJ4nIJhFZJCKvioi/E3cPnLdnz5xu\nqroYWAP09BYdDWxT1dm+7cbhDuZRPpl5qrrFJzMG167gUJ/MJFXNC5NpLyJ1fTLjwvQf49PFMJKW\n7dth/Hi48MLSbT/7/9u78/ioynOB47+HsC9h30UBUdxQNlFEQEVxxaJVEbFasdd6tdbSq93srVZb\na1sr2uuCdV8QVFxQRFBAQRBkR3ZRdpBASDIJ2ZN57h/vRIcYQoacmTMneb6fz3zCzLw573NIzuQ5\n7xoZjG0z1WouVRg/3q2c3sOr/oUEeTYri+GbNnFMLdkM2FQsluSoE7Cq7ImqbsK11nQ85Hd450Pg\nBuA84DfAUGC6yHfLiXUAilQ1u9z3pUXeKyuzN/pNVS0FMsqVKT8CLy3qveqWSU1gMmlMXEybBsXF\nRz7eaElODo3q1OFkG4xdY82eDevWwa9+5XcksVnzzTcs6tSJ/7LutFov1tlqpRU8j/t6p6r6RtTT\ntSKyGrfG0jnAJ/Guv4oCtO6rMUfurbdg4EDo3PnIvn9xdjb9mjWjrg3GrrH+7//glFNg6FC/I4nN\ns198QdvGjRlhA7FrvViSIwG+EpHoAdNNgRUiEi57QROw8ayqbhGRdKAHLjnaA9QXkdRyrUftI+8R\n+Vp+9loK0KpcmfL7RbePeq/sa/sKymgVymSraqULZY4bN47mzZsf9Nro0aMZPbr8+HNjEi83F2bM\ngAceOPJjLM7J4Yo2bTyJZ9KkST8YyBwKhTw5tjkyW7bA++/DhAnB2iqkoKSEV1JTuXnrVurbiti1\nXizJ0U1xiyJGInIU0Bo3QBxgGW5g9TDcIpWISE/gaGBhpMxCoIWI9IkadzQMl/R9EVXmDyLSJmrc\n0XAgBKyLKvMXEUmJdMuVldmoqqGoMheXC3t4VCyHNH78ePpaX7dJUtOnQ0HBkY832ldUxJaCAk73\naLxRRTcOy5cvp58t3OebJ590++2NGeN3JLGZOn8+GU2bMrZXL79DMUmgysmRqr4UryAiaw314Puu\nqe4ichpuPFAGcC9uWv6eSLm/A1/hBjmjqtki8hzwiIhkAjm4DXAXqOriSJkNIjITeEZE/huoj1vh\ne5KqlrX4fIRLgl4Rkd/ixlM9ADyuqsWRMq8BfwKeF5G/A72AXwJ3Rp3SBOD2yPvP45Kwq4BLvPj/\nMsYvb74JfftC9+5H9v1Lc3IAGJCa6mFUJlnk5cFzz8HNN0OTJocvn0xe2LmTswoKOOHmm/0OxSSB\nKnX6Rw18jpf+wApcC5AC/wKWA3/GjWs6FZgKbASeAZYAQ6ISFnDLC0wDpgCfArtxax5Fuw63BMGs\nSNl5uKn6AKhqGLgsUufnuC1RXsQlZ2VlsnGtQF2BpcA/gftU9bmoMltx6yidj1u/aRxws6qWn8Fm\nTGDk5cEHH8DVVx/5MRbn5NC6bl26NWzoXWAmabz+upvNeNttfkcSm53p6XzUqZPrHglSX6CJm6q2\nHK0Vkftxaw0VHaqQiBwH/BrYpqoPVTUIVZ1L5YnaRVU4RiFwR+RxqDJZwPWHOc4OXIJUWZk1uBlz\nlZWZh1tiwJgaYfp0lyBVKznKzmZAairxv98yfnjqKbfX3pG2LPrl5blzadikCddcYOv0GqeqydEd\nuK6sJ0XkY1yLyW6gAGgJnAScjVsL6HHgKe9DNcb46c03oU8fOPbYI/t+VWVJTg7/3amTt4GZpLBs\nGSxZ8v3q6UGhqrwAXLVpE6kXHfY+3NQSVUqOVHU20F9EzgZGAWOAY4BGuNWiV+C6oCaqamacYjXG\n+CQvz61vdM89R36MbQUF7CsutvFGNdRTT0GXLhC0WfAL1q/n69at+U84fPjCptaIaZ0jVZ0PzI9T\nLMaYJPXhhx50qUUGY3s1U80kj6wseO01+MMfICXF72hi8/LKlRyTksLQyyodTWFqGVuFzRhzWG+8\nAb17w3HHHfkxFmZn071hQ9rVr+9dYElKRBqJyFYR+YffsSTCxIlQVORmqQVJfkkJb7RowU/S0qgT\ntOl1Jq4sOTLGVOrAAbeo36hR1TvOouxszqw9XWr3UIV1zWoCVXjmGRgxAjomYjMpD73/+eeEGjfm\nJ6ee6ncoJslYcmSMqdS0aZCfD9dcc+THKAyHWZ6Tw8BakByJSA+gJ25PyBpv2TJYtQp+9jO/I4nd\nyzt2MHDTJo4fMsTvUEySseTIGFOp11+H00+v3vTsFTk5FKnWlpajh4HfU0v2W3zmGbfPXtAmeqUV\nFDCjQwduKCkB2+fPlGO/EcaYQ8rOdoOxq9ultjA7m4Z16nBqEu92LiKDReQ9EdklImERubyCMreL\nyBYRyReRRSJyern3L8dtJfR12UuJiN0vBw64gdhjxwZvIPakffuoU7cu11x1ld+hmCRU5eRIROqI\nyG9EZIGILBGRh0SkUTyDM8b4a+pUKCysXpcauPFG/Zs1o35y36E3wa1ofxtupf6DiMgo3Or99wJ9\ngFXATBGJ3kX3TOBaEdmMa0H6mYj8Md6B++XNN91mxGPH+h1J7Cbu3culrVvTqm1bv0MxSSiWT6p7\ngAdx+5btwu0l9kQ8gjLGJIfJk2HQILd+TXUEYTC2qs5Q1T+p6lQqbvEZBzytqi+r6gbgViAPGBt1\njD+o6jGq2h24C3hGVf+SiPj98MILMGwYdO3qdySx+Sovj6U5OYxp187vUEySiiU5ugG4TVUvUtWR\nwAhgjIgk9a2gMebIpKfDRx/BdddV7zi7CwvZXliY9MlRZUSkHm47oNllr6mq4vZpHOhXXH76+mv4\n7DP46U/9jiR2E9PSSE1J4bLWrf0OxSSpWBaBPJqo2ReqOktEFOgE7PQ6MGOMv6ZMcdO0q7PwI7hW\nIyDoM9XaAClAWrnX03Az035AVV+qyoHHjRtH8+bND3pt9OjRjB49+gjCTJyXXoLUVLjiCr8jiY2q\nMjEtjR+3bUvDoA2UMgBMmjSJSZMmHfRaKBTytI5YkqO6uL3UohUD9bwLxxiTLF57DS64AKo7JOPz\nUIguDRrQqUEDbwKrYcaPH0/fvn39DiMm4bBLjkaNgsaN/Y4mNotzcvimoICn27f3OxRzhCq6eVi+\nfDn9+nm313ssyZEAL4pIYdRrDYEJIpJb9oKqXulVcMYYf2zf7rpMXn65+sdakJ3NoHItIwGUDpQC\n5f+itgf2JD4cf82ZAzt2wE03+R1J7CalpdGxfn3OadHC71BMEoslOaqoifhVrwIxxiSP11+Hhg1h\n5MjqHSe/tJRlOTlcH/C7dFUtFpFlwDDgPQARkcjzf1fn2GXdakHoSivz8stw/PFw5pl+RxKbUlXe\n2LePUe3akSI1epWFWqOsi823bjVVDeA9gjHmSEycCJddBtXdI3ZJTg7FqgwKwHgjEWkC9OD7mWrd\nReQ0IENVdwCP4FrPlwGLcbPXGgMvVqfeoHWr5ebC22/D734HQcsv5mVl8W1REdfaLLUao+ymws9u\nNWNMLbB6tdsO4v77q3+sBaEQzVJS6JXEiz9G6Q98glvjSHFrGoFrNR+rqm9E1jS6H9edthK4UFX3\n+RGsX9591yVIY8b4HUnsJu/dS9eGDRlQ3azf1HiWHBljDjJxIrRq5c12EAtCIc5MTQ1EF4aqzuUw\ny5uo6pPAk4mJKDm9+iqcfTZ06+Z3JLEpDoeZsm8ft3TqhATg99H4y5IjY8x3wmGXHI0aBfXrV/NY\nqnyenc2vjjrKm+BqqCCNOdqzx6199WQA08NZmZlklJRYl1oN4/uYI2NMzTdvHuzcCddfX/1jrc/L\nI7OkpCbMVIurII05mjwZ6tat/nYyfpi8dy89GzXi1CZN/A7FeCheY45sdWtjzHdefdV1lwz0YM3n\nBaEQKcAZNr6jxnjtNbjkEmjZ0u9IYlMYDjM1PZ1r2rWzLjVTJZYcGWMAyM93G4lef703s5AWhEKc\n1rQpTetaA3VN8M03sGQJJHnPX4VmZWYSKi3lattk1lSRfWoZYwA3Cyk7G264wZvjzQuFGNmmzeEL\n1nJBGXM0eTI0aeKWeAiaNyNdaqdYl1qNY2OOjDFx9dJLMGgQ9OhR/WNtLyhga0EBQ2280WEFZczR\n5Mnwox8Fb7uQonCYqfv384vOna1LrQayMUfGmLjZvRs+/hhuvNGb483LygLgbEuOaoQ1a9wjiRu2\nDmlWZiZZJSXWpWZiYsmRMYaJE6FePbj6am+ONy8U4uTGjWlT3fUATFKYPNkNwh4+3O9IYjdl3z6O\nb9SIXtalZmJg3WrG1HKqrktt5Ejwai/OeVlZnBe0KU0+SfYxR6pur70rr6z+2leJVhyZpXarLfxY\nY9mYI2NMXCxdCmvXwsMPe3O8PYWFbMzP589BW0LZJ8k+5mjVKvj6a3j8cb8jid28UIiMkhJ+bF1q\nNZaNOTLGxMXzz0PnznDBBd4c77PIHdxgG29UI0yZ4rrUzjvP70hi9/a+fRzToAF9grG3n0kilhwZ\nU4vl58OkSW4gdkqKN8ecFwrRo1EjOjVo4M0BjW9U3dpXI0e6MWlBElblnfR0rmzb1rrUTMwsOTKm\nFnvnHQiF4KabvDvmp1lZDLFWoxph9Wr46ivvBuon0hfZ2XxbVMSVttaWOQKWHBlTiz3/PAwe7M3a\nRgB7i4pYk5trg7FriClT3CD9YcP8jiR2b6en075ePQZaom6OgA3INqaW2roVZs92CZJXPo2sb3Su\nV9PeaoFknq02ZQpcfnnwZqmpKm/v28eP2rQhxbrUajSbrWaM8dRzz0GzZt7usD4nM5OeNt4oJsk6\nW23DBli/Hh56yO9IYrc2N5fNBQVcYV1qNZ7NVjPGeKakBF54Aa67zu2X5ZVPbH2jGuOdd9xWIV7N\nYkykd9PTaZaSwrn2u2iOkCVHxtRCM2bArl3wX//l3TF3FRbyVX6+danVEO++CxddBI0a+R1J7Kbu\n389FrVrRoI79iTNHxn5zjKmFnn0WevcGL3tzPsnMBOAcS44Cb9cuWLwYrrjC70hit6uwkKU5OfzI\nutRMNVhyZEwt8+23MG2aazXycqzqnKwsejVpQtugjd41P/Duu1C3Llx6qd+RxO699HRSgEtatfI7\nFBNglhwZU8s895ybfXTddd4dU1WZk5lpXWo1xLvvwrnnupWxg2ZqejpDW7SgZdBWrTRJxWarGVOL\nlJbCM8/A6NHebTIL8HV+PtsKC7kgiH9NfZZsU/mzsuDTT+Gxx/yOJHY5JSXMycri4WOP9TsUkyA2\nld8YU20zZsD27XDrrd4e9+PMTOqKMNRajmKWbFP5Z8xwsxlHjPA7kth9lJlJsSojWrf2OxSTIDV6\nKr+IDBaR90Rkl4iEReTyCsrcLyK7RSRPRD4WkR7l3m8gIk+ISLqI5IjIFBFpV65MSxGZKCIhEckU\nkWdFpEm5Ml1E5AMRyRWRPSLyDxGpU67MqSIyT0TyRWSbiNxdQbzniMgyESkQka9E5Mbq/S8ZU30T\nJrhB2P37e3vcjzMzGZiaSrO6dr8VdO+/7wbrd+nidySxm7Z/Pyc3bky3IE6xM0klKZIjoAmwErgN\n0PJvishvgV8AtwADgFxgpohEj/x8FLgU+DEwBOgEvFXuUK8BJwLDImWHAE9H1VMHmI5rUTsTuBH4\nKXB/VJlmwExgC9AXuBu4T0R+FlWmKzANmA2cBjwGPCsiAVwxxNQUO3bA9Omu1cjLgdgl4TBzMjOt\nS60GKC52vyNBbDUKq/LB/v1cZq1GxgNJcZunqjOAGQBS8fbJdwIPqOq0SJkbgDRgJPCGiKQCY4Fr\nVXVupMxNwHoRGaCqi0XkROBCoJ+qroiUuQP4QETuUtU9kfdPAM5V1XRgtYj8L/CQiNynqiXA9UA9\n4ObI8/Ui0gf4NfBsJN7/Bjar6m8izzeKyNnAOOBjb/7XjInN00+7BR+9HtayJCeH7NJSS45qgAUL\n3JijICZHS3Jy2FdcbMmR8USytBwdkoh0AzrgWmEAUNVs4AtgYOSl/rhEL7rMRmB7VJkzgcyyxChi\nFq6l6oyoMqsjiVGZmUBz4OSoMvMiiVF0mZ4i0jyqzKxypzIzKhZjEqqw0A3EvvFGaNrU22N/lJlJ\n85QU+jdr5u2BTcK9/z506AAeDt1ImPfT02lVt65tNGs8kfTJES4xUlxLUbS0yHsA7YGiSNJ0qDId\ngL3Rb6pqKZBRrkxF9eBRmVQRsU2nTMK99Rbs3Qu33eb9sT/OyOC8li2pa6sRB97777tWoyD+KKft\n38/5LZqxOm0Vqj8YnWFMTAJ4CSQt2/rZJK0nnoDzzoMTT/T2uKGSEr7IybEutRrgq69g0ya47DK/\nI4ndjoICVuXm0jJ3A32f7su+vH1+h2QCLinGHB3GHlzi0Z6DW2PaAyuiytQXkdRyrUftI++VlSk/\ney0FaFWuzOnl6m8f9V7Z1/YVlNEqlMlW1cIfnuL3ytY8iZYs65+YYFq5Ej7/3LUeeW12ZiYlqlzk\n02rEZWucRPN6vZN4S5Z1jj780C0OOmyYbyEcsQ8zMkgBtmx5k0FHD6Jdk3aH/R5TM8RrnSNUNake\nQBi4vNxru4FxUc9TgXzg6qjnhcAVUWV6Ro41IPL8BKAU6BNVZjhQAnSIPL8IKAbaRJW5BcgE6kWe\n3wqkAylRZR4E1kU9fwhYVe4cXgOmV3LefQFdtmyZGuOlm29W7dxZtbjY+2P/bMMG7blokfcHroZl\ny5Yp7malrybBZ9qhHsl2zQ8frnrBBX5HcWRGrl6tA5cu0QYPNNCHFzzsdzjGB15f90nRrSYiTUTk\nNBHpHXmpe+R52UobjwJ/FJERItILeBnYCUyF7wZoPwc8EllfqB/wPLBAVRdHymzADYp+RkROF5FB\nwP8Bk9TNVAP4CFgHvBJZy+hC4AHgcVUtjpR5DSgCnheRk0RkFPBL4F9RpzQhcg5/F5GeInIbcBXw\niHf/a8YcXno6TJwIt9/u9srykqry4f79XGyzgwIvLw/mzoVLLvE7ktgVhcPMysykWziNwtJCRp4w\n0u+QTA2QLN1q/YFPcFmf8n2i8RIwVlX/ISKNcWsStQA+Ay5W1aKoY4zDtQxNARrglga4vVw91wGP\n42aShSNl7yx7U1XDInIZ8BTwOW49pReBe6PKZIvIcOAJYCmuFek+VX0uqsxWEbkUGI9LnHbipv6X\nn8FmTFw98wyouk1mvbYmN5ddRUVcbBt8Bt4nn7gZjRdf7HcksZsfCnGgtJTQtx9zSrtTOLaVbR1i\nqi8pkiN1axNV2oqlqvcB91XyfiFwR+RxqDJZuHWKKqtnB1DpkERVXQMMPUyZeUAAJ8SamqK4GJ58\nEsaMgTZtvD/+hxkZNK5ThyE2dTrwpk+H7t3h+OP9jiR20/fvp2P9+ixY/zK3nx6H6ZimVkqKbjVj\njPfeeQd27oRf/jI+x/8wI4NzW7SgYUpKfCowCaHqBmNffLG3K6cnyocZGfSpX0RWQSY/6vkjv8Mx\nNYQlR8bUUI8+CkOHwmmneX/s7JIS5odCNt6oBvjqK9iyJZhdatsKCliXl0edzKV0ataJfp2ssd54\nIym61Ywx3lq40D2mTo3P8T/KyKBElUtsvFHgzZzppvCfc47fkcRuRmQK/9qvXuay4y6jjtj9vvGG\n/SYZUwONHw89esRvQb/3bffzGuOjj2DwYLfvXtDMzMjgtMb12JK+mhE9A7ghnElalhwZU8Ns3eoW\nfBw3Lj7bQJSqMj0jgxHxGOVtEqqw0M1UGz7c70hiVxIOMzszkzYFW2hUtxHDugVw9UqTtKxbzZga\n5t//hubN3Saz8bAoO5v04mIut/FGgbdwoVvj6IIL/I4kdl/k5JBdWsreXR9wfvfzaVTPWjGNd6zl\nyJgaJCvLrW10663x6yZ5Pz2dtvXqMSA1NT4VmIT56CNo2zY+g/bjbWZGBi3rprDqmylcdnwAN4Qz\nSc1ajoypQSZMcOsbxWv6PrjxRpe2bk1KEOd9JyE/91b76CPXahSP7td4+ygjg551DrBIS7j0uEv9\nDsf4JF57q1lyZEwNUVAAjz0GN9wAHTrEp47N+fmsy8vjgW7d4lNBLTR+/Hj69u2b8Hr37YPly+GO\nQy6bm7wyiotZkpPDgANL6N2hN51TO/sdkvFJ2U3F8uXL6dfPu6UcAni/YIypyKuvQloa/M//xK+O\nqenpNBBheMuW8avEJMTs2W4ByCCON5qVmUkY2LDpVS7pEcAN4UzSs+TImBqgtBQefhhGjoSePeNX\nz9vp6Qxv1YqmXu9iaxJu9mw48UTo1MnvSGI3KzOTrvWErOxNXHq8dakZ71lyZEwN8M47sHEj/Pa3\n8asjraiIBaEQV9gU/hphzhwYFtDZ77MyM2lbtJ2WDVtyRucz/A7H1EB2+2dMwKnCgw+6P3RnxPHv\nxNT0dOoAI2wKf+Bt3QqbNwczOdqcn8+WggLq7P6YC3tcSEod29vPeM9ajowJuJkzYcUK+MMf4lvP\n2/v2MaRFC9rUrx/fikzczZnjNpkdOtTvSGI3KzOTFOCbLW/ZeCMTN5YcGRNwDz7oWozOPTd+dWQV\nFzM7K4srrUutRpgzB/r2hSCOq5+VmUn3lEIozWP4sQFc2tsEgnWrGRNgc+fCZ5+5DWbjuezQtP37\nKVFlpCVHgafqBmPfcIPfkcQurMrszEw6HVhDv479aN+0vd8hmRrKWo6MCbAHHnCrG4+I856bb+7b\nx8DUVI5q2DC+FZm427AB9uyB887zO5LYrTxwgIySErZveZeLelzkdzimBrPkyJiAWrDAtQD86U/x\nbTXKKi5mRkYG17RtG79KTMLMmQP16sHZZ/sdSexmZWbSSCA7fREXHnuh3+GYGsySI2MC6oEH4JRT\n3NpG8fTe/v0UqXJ1u3bxrcgkxCefwIAB8dt7L55mZ2ZylGaQWr8xZx51pt/hmBrMkiNjAmjRIjdL\n7Y9/jP++WK/v3cvZzZvTuUGD+FZk4k4V5s0L5iy1onCY+aEQxRmLOa/bedRLqed3SKYGs+TImAD6\n3/+Fk0+Gq6+Obz2ZxcV8lJnJKOtSqxHWr3d7qgUxOVqcnU1eOMz2re9Zl5qJO5utZkzAzJsHs2bB\nlCnxbzV6Nz2dUlV+bMlRlYhIc2AWkIL7fP23qj7rb1TfmzsXUlLgrLP8jiR2n2Rl0ViUvJyNNoXf\nxJ0lR8YEiKprNerdG664Iv71vbZ3L+e0aEFH61KrqmxgsKoWiEgjYK2IvKWqmX4HBi456t8fmjb1\nO5LYfZKVRbvi3aS07Eb3lt39DsfUcNatZkyAzJ7tWo7uvz/+rUa7CwuZnZnJ9e1tLZmqUqcg8rRR\n5Gsc5xJWnapLjoLYpVZQWsrnoRAH9s7ngu4X+B2OqQUsOTImIFTh9793q2Ffdln865u0dy/1RaxL\nLUYi0lxEVgLbgX+qaobfMQF8/bVb3yiIydHC7GwKVUnf/RHndz/f73BMLWDJkTEB8fbbsHQpPPRQ\nfNc1KvNqWhoj2rShed3a0fsuIoNF5D0R2SUiYRG5vIIyt4vIFhHJF5FFInJ6+TKqGlLV3kA3YIyI\nJEV2OXeua20M4vpGc7KyaCqlkLuV87oFcPVKEziWHBkTACUlbtr+8OFwzjnxr2/NgQOsPHCg6/V6\noQAAIABJREFUtnWpNQFWArcBWv5NERkF/Au4F+gDrAJmikiFe6qo6r5ImcHxCjgWc+dCnz6Qmup3\nJLH7NCuLFgXb6N+pLy0bBXBDOBM4lhwZEwAvveS2fXjwwcTUN3HvXlrWrcvFrVolpsIkoKozVPVP\nqjqViscJjQOeVtWXVXUDcCuQB4wtKyAi7USkaeTfzYEhwMb4R3948+fD4KRI02KTX1rK4uxssvZ8\nyvndrEvNJEbtaC83JsDy8twWIaNGQb9+8a+vVJVX9uxhVLt21I/3qO+AEJF6QD/gu/RUVVVEZgED\no4oeA/xHXL+nAI+p6trKjj1u3DiaN29+0GujR49m9OjRHkUPu3bB1q3B7FJblJ1NkSpF+xZw/uD/\n+B2OSQKTJk1i0qRJB70WCoU8rcOSI2OS3PjxbuG+RLUafZyRwa6iIsZ26JCYCoOhDW7torRyr6cB\nPcueqOoSXJdblY0fP56+fftWO8DKLFjgvg4aFNdq4mJuVhaNKaGkYDdndQngAk3GcxXdPCxfvpx+\nHt49WnJkTBLbuxf+/nf4xS+ge4KWdnl+zx5OadKE/s2aJaZCE3fz58Oxx0IQ8925oRDNCrZy8tFn\n0aheo8N/gzEesOTImCT25z+7GUb33JOY+vYXFzM1PZ2/de+OJGJKXHCkA6VA+RHq7YE91TlwWbea\n111p0ebPD2aXWmE4zKLsbHTPXM7rbrPUzA+VdbFZt5oxtcSaNTBhAvzzn9C6dWLqfC0tjTDUtllq\nh6WqxSKyDBgGvAcgLnscBvy7OseOd7daTg6sWgX//d9xqyJulmRnUxAOw/7FDDv/dr/DMUmo7KbC\nutWMqQVU4de/dl1pv/hFoupUnvv2Wy5r3Zp29esnptIkIiJNgB58P1Otu4icBmSo6g7gEeDFSJK0\nGDd7rTHwog/hVtmiRRAOB7PlaG4oRANKqFu0h34dEzAbwZgIS46MSULTp8PHH8PUqZCoPGVxTg6r\ncnP5W6IGNyWf/sAnuDWOFLemEcBLwFhVfSOyptH9uO60lcCFkfWMktb8+a7l8YQT/I4kdnOzsmiW\nv5UBx5xNvZR6fodjahFLjoxJMoWFMG4cDBsGI0Ykrt4Ju3fTtWFDhteitY2iqepcDrP2m6o+CTzp\nZb3xHnO0YAGcdVZiVlX3Ukk4zOehEEV753Pucef6HY5JUjbmyJhaYvx42LwZ3nkncX/QMouLmbx3\nL3865hhSgvZXNODiOeaotBS++CJxA/q9tCo3l9xwGDKWcW7XO/wOxyQpG3NkTC2wcyf85S/wy1/C\nyScnrt6X09IoUWVsx46Jq9TE3dq1cOAADBx4+LLJZn4oRF3CNC76lt4devsdjqllLDkyJoncdRc0\nbQr33pu4OlWVCbt3c2WbNrSvhQOxa7JFiyAlBfr39zuS2H2WlUXTgp2cffRAUuqk+B2OqWUsOTIm\nSXz8Mbz+uttHrdxuEnE1KzOTDXl5PHXccYmr1HwnnmOOFi6EU0+FJk08PWzcqSrzQyEO7FvI0GOH\n+h2OSWLxGnMUiI2TROReEQmXe6wrV+Z+EdktInki8rGI9Cj3fgMReUJE0kUkR0SmiEi7cmVaishE\nEQmJSKaIPBuZ3htdpouIfCAiuSKyR0T+ISJ1ypU5VUTmiUi+iGwTkbu9/j8xNUtBAdx2G5xzDvzk\nJ4mt+7GdOzm1SROGtmiR2IoN4MYcvffee3EZjL1oEZx5pueHjbtv8vNJKy6mJGsFQ4+x5Mgc2ujR\no3nvvfcYP368p8cNRHIUsQY3fbZD5PHdqh0i8lvgF8AtwAAgF5gpItF9BI8ClwI/xu2U3Ql4q1wd\nrwEn4hZ2uzRS7umoeuoA03EtbmcCNwI/xU3tLSvTDJgJbAH6AncD94nIz6px7qaGe+gh2LYNnnwy\nsbOKNuXl8UFGBncedZStiF3DZGTAhg3BHW8kKE3zt9GnY0xb1RnjiSB1q5VUsp7IncADqjoNQERu\nwG0IORJ4Q0RSgbHAtZHpuojITcB6ERmgqotF5ETgQqCfqq6IlLkD+EBE7lLVPZH3TwDOVdV0YLWI\n/C/wkIjcp6olwPVAPeDmyPP1ItIH+DXwrPf/LSboNmyAv/0N7r4bTjwxsXU/vmsXberV47p27Q5f\n2ATK4sXuaxBbjj4LhWhStI+zO/embp0g/ZkyNUWQWo6OE5FdIvKNiLwqIl0ARKQbriVpdllBVc0G\nvgDK7pn64xLB6DIbge1RZc4EMssSo4hZuMXgzogqszqSGJWZCTQHTo4qMy+SGEWX6SkiCRxJYoIg\nHIaf/QyOOQb++MfE1h0qKeGFPXv4eceONEyxAa81zaJFbvHHHj0OXzbZzA+FKMxYYl1qxjdBSckX\n4bqvNgIdgfuAeSJyCi4xUlxLUbS0yHvguuOKIknTocp0APZGv6mqpSKSUa5MRfWUvbcq8nVzJWW8\nHTVmAm3CBLdI36efQqMEbzj+9O7dFIbD3N65c2IrNgeJ14DshQtdq1HQekv3FRXxVX4+ZK5g8Fl/\n9jsck+Rq9SKQqjoz6ukaEVkMbAOuATb4E5Ux1bN9O/zud3DLLTA0wTfIheEwj+7cyQ0dOtCxQYPE\nVm4OEo9FIMNht/jjXXd5etiEWJjt7mEb5G7i9M6n+xyNSXa2CGQUVQ2JyFe4TSI/xW0U2Z6DW3Xa\nA2VdZHuA+iKSWq71qH3kvbIy5WevpQCtypUpf7W2j3qv7Gv5Lc3LlzmksrvIaPHaVsD4R9V1pzVv\nDv/4R+Lrf2XPHvYUFXFXly6Jr9xDZXeN0by+gwyiTZsgFIIBA/yOJHafh0I0KD3AmW27Uz/F1t0y\n/ghkciQiTXGJ0UuqukVE9uBmmH0ZeT8VN07oici3LANKImXeiZTpCRwNLIyUWQi0EJE+UeOOhuES\nry+iyvxBRNpEjTsajusqWxdV5i8ikqKqpVFlNqrqYT+147mVgEkezzzj1jWaMSOxaxoBhFX5544d\njGzThp6NGye2co9VdOPg9R1kEC1Z4r4GcfHHhdnZaGg1Q44e7HcophYLxIBsEfmniAwRkWNE5Cxc\nglMMTI4UeRT4o4iMEJFewMvATmAqfDdA+zngERE5R0T6Ac8DC1R1caTMBtzA6WdE5HQRGQT8HzAp\nMlMN4CNcEvRKZC2jC4EHgMdVtThS5jWgCHheRE4SkVHAL/l+h29Ty23bBv/zP67l6MILE1//2/v2\n8VV+Pr8JeKuRObQlS9xA7KDtIVwcDvNFdoiizBUMtuTI+CgoLUdH4ZKO1sA+YD5wpqruB1DVf4hI\nY9yaRC2Az4CLVbUo6hjjgFJgCtAAmAHcXq6e64DHcbPUwpGyd5a9qaphEbkMeAr4HLee0ovAvVFl\nskVkOK7VaimQDtynqs9V+3/BBF5pKdxwg/uj9fDDia8/rMr927ZxfsuWnJnoJitToXgMyF68GE4P\n4HCdlQcOUKgg2es486gArkFgEq62D8g+7CeGqt6Hm8V2qPcLgTsij0OVycKtU1RZPTuAyw5TZg1g\nc1DND/zrX/DZZ252mh+5ydT0dFbn5vKkbRWSNLzuSi8uhpUrYdQozw6ZMJ9nZ1NHSzitaROaNWjm\ndzgmAOI1IDsQ3WrG1AQrVri1jO6+G4YMSXz9Gmk1OrdFC862rUJqrDVr3HY0QWw5+jwUol7eFgZb\nq5HxWSBajowJugMH4Npr4eST4f77D18+Hqamp7PywAE+7d3bnwBMQixZAikp0CeAu27Mz8qkMGMZ\ng/ud73coppaz5MiYBPjlL2HXLli2DPxYVqhUlXu2bOG8Fi1sg9kabvFiOOUUCNpExJ0FBewuLoHQ\nWgZ1uffw32BMHFlyZEycvfYavPCCe/Ts6U8Mr6alsS4vjxdPOMGfAEzCLFkSzPWNFkUWf+wiB+jY\nrKPP0ZjazsYcGRNH69e7FbDHjIEbb/QnhsJwmD9t2cJVbdtyemqqP0GYhMjLg7Vrgzne6IucHOoX\nZzKk48mHL2xMnFnLkTFxkpsLV13lNpWdMMG/Pa6e2rWLXYWF/KVbN38CMJXycir/qlVuuYggroG5\nMJRFcdaXDOo6yO9QTIDU6qn8xgSNqmsx2rbNdXM0bepPHPuLi/nztm38rGPHwK+GXVN5OZV/+XKo\nV8+NOQqS4nCYpTk5aPZazuoy0u9wTIDY3mrGBMj48W6s0eTJcOKJ/sVx75YthFV5wFqNaoVly6BX\nL38G/VfHmtxcChUa5W/jlHYBy+xMjWRjjozx2Jw5bi2ju+/2dyG+tbm5TNi9mz917Urb+raBZ22w\nfDkEcWvGL7KzES3ljOYtSamT4nc4xlhyZIyXvv4arr4azjsPHnzQvzhUlTs3baJrw4b8onNn/wIx\nCVNQ4AZjBzE5WpSdTZ28bQw+KoDT7EyNZMmRMR7JyoIRI6BNG3jjDajrY6f1pL17mZ2VxePHHUeD\nOnaZ1warV0NJSTCTo8+y9lOa9aXtp2aShn1qGuOB4mLXhZaWBu+/Dy1b+hdLVnExv/76a65q25aL\nWrf2LxCTUMuXu5WxTz3V70hik1lczObCEshZb8mRSRo2INuYalKFn/8cPvkEZsyA44/3N557tmwh\nNxzm0R49/A3EJNTy5XDSSdCokd+RxGZpTg4AXevk06pRK5+jMcax5MiYarr/frf69SuvuLFGfpqb\nlcWTu3fzWI8edA7alKVayqt1joI6GHtxTg51SvMZ0q6736GYALJ1joxJQhMmwH33wV//Ctdf728s\nuaWl3LxhA2c3b26DsAPEi3WOiorgyy/hJz/xKKgEWhTKJJy9noFdrEvNxC5e6xzZmCNjjtDrr8Nt\nt7lNZX//e7+jgXs2b2ZXURHP9+xJHb+W4za+WLfOJUhBbDlaFMqCnA2c0fkMv0Mx5jvWcmTMEZg+\n3d2ljxnjFnz0Oxf5OCODx3btYvyxx3KcrYRd66xc6b6edpq/ccRqd2Eh6aVQP28Lvdr38jscY75j\nLUfGxOijj+DKK+HSS+H558HvmfL7ioq4ccMGhrdsyS+POsrfYIwvVq2CY4+FZs38jiQ2SyKDsU9r\n3JC6dexe3SQPS46MicHs2TByJJx/vutWq1fP33hUlZs3bqRYlRdPOMG602qpVaugd2+/o4jdkuxs\n6hSHGNLhBL9DMeYglhwZU0UzZsBll8GQITBlCiTDjhz/2rGD9/fv5/mePelos9NqJVXXrRa0LjWA\n+VnphLPXcUZnWxnbJBdLjoypgqlT4Uc/ggsucP9u2NDviNy0/d9t3szvjj6aEW3a+B2O8cnOnZCZ\nGbzkSFVZnnPADcY+ygZjm+RiyZExh/HSS/DjH8Pll7sWo2RooNlZUMCotWsZ0qIFD3Tt6nc4xker\nVrmvQetW21xQQI7WoUVxGl1Su/gdjjEHseTImENQhX/9C376Uxg7FiZPTo6utAMlJYxYs4YGdeow\n+aSTqOv3iHDjq5UroUUL6BKw/KJsZewBqc0RGytnkox9qhpTgdJSt37RXXfBH/4ATz/t9q3yW6kq\nY9av5+v8fKb16kW7ZMjWjK/KBmMHLb9Ymp2NFO7l7A6n+B2KMT9gcyeNKScnx61fNH26S4puucXv\niBxV5c5Nm5i2fz/v9epFr6ZN/Q7JJIFVq+CSS/yOInbzM/ehORsZcKwNxjbJx5IjY6Js3uzGFu3Y\nAe+/Dxdf7HdE37t/2zae2L2b/xx/PJe2bu13OCYJHDgAX38dzMHYX+blQ85X9O90m9/hGPMDlhwZ\nEzFjhmsxatkSFi2CE0/0O6LvPbZzJ/dt3cpfu3Xjvzp18jsc46HqbDy7erUbGxe0wdhbCgrI0zp0\nIJvWjS3RN0fONp41Jk5KS+GBB+D++11L0SuvQKtWfkf1vUd37GDcN9/wmy5d+P3RR/sdjvFYdTae\n/fJLNxbupJM8DirOlkUGY5+RmkQXmgmkeG08a8mRqdV27YLrr4d58+DPf4Z77vF/O5BoD2/fzt2b\nN/ObLl14qHt3m9VjDrJ6NRx/fHIsLxGLpTnZSGE6gzoELKsztYYlR6bWeust+PnP3YKOc+bA0KF+\nR/Q9VeW3mzfzzx07+MPRR/OXbt0sMTI/sGYN9Argfq2fZaShORs5vecgv0MxpkJJdI9sTGJkZsIN\nN8BVV7mtQFauTK7EqDAc5sYNG/jnjh082qMHf7UWI1MBVddydErAZsKrKqvzCuDAV/TteGTdicbE\nm7UcmVpD1a1wfccdUFDgxhaNGZNc68PsLSriijVrWJaTw2snnsjo9u39DskkqT17ICMjeC1H2wsL\nOaApdJY8Uhuk+h2OMRWyliNTK3zzjds09ppr4KyzYN06N9YomRKjhaEQ/ZctY3NBAXP79LHEyFRq\nzRr3NWgtR8sjg7FPT23hcyTGHJolR6ZGy8mBP/4RTj7Z/TF55x14+21IptnwYVUe3r6dIStXcnSD\nBizp25czUu2O2lRu9Wpo1Ai6d/c7ktgsyc6GogwGd0iitTKMKce61UyNVFwMzz8P994LoRDcfTf8\n/vfQuLHfkR1sa34+N23cyKdZWdzdpQt/7daNesk0Xc4krTVrXNIftF+XBZlpkPMV/Y4f6HcoxhyS\nJUemRikpgddec9Pyt2xxY4oefDD5NuUsLS3liXfe4Z6WLWnVuDGzTzuN81q29DssEyBBHIwN8GVu\nHuR+Q5+OtjK2SV4Bu+cwpmIFBfCf/0DPnnDjjW6Q6qpVbtB1siVG87/4ggGTJ/OrVq24fto0vjz9\ndEuMTEzCYVi7NniDsfcWFZGldemgB2wwtklqlhyZQPv2W9dKdPTRcOut0LcvrFgB776bfH84NmzY\nwI9feIHB+fkQDrPw9dd56rnnaF7XGnBNbDZvhvz84LUcrTpwAIC+zSwxMsnNPpVN4JSWwuzZ8Oyz\nboB1/fpw001w551w3HF+R/dDa9as4W+LFzOpa1eOataMl/fsYcy111InOxtefNHv8EwAlc1US7Yb\ngMNZnpMDpXkMbnes36EYUylLjkwglC14N2kSTJwIO3a4jWEfecQt6Ni8ud8RHqy0tJSZn3zCv3fu\nZGbXrhzVvDlP7NrF2CuvpEGTJq5QgwZQVOROLpnWFDBHTESOAl4B2gHFwF9UdYrX9axZ4/b/69DB\n6yPH12cZ38KBb+jXtY/foRhTKUuOTNIKh2HJEpg61U2/37jR/UG4+mrXUjRgQPLlFJs2bWLiwoW8\n0KQJ21u3pk9KCq+kpXHN5ZdTv1GjgwuXbYhVWOj2MDE1QQlwp6p+KSLtgWUi8oGq5ntZybp1brPZ\nZPv9P5yVBw7AgU306Tja71CMqZSNOYojEbldRLaISL6ILBKR0/2OCWDSpElJW8eOHa6n6Sc/cXfF\nZ57pBloPGgTTprkxRhMmwBlnuD8MiTgXOHQ9Gg6zZvVq/vrqq/R/9VWO37WLR9q04fz9+/kiJYVl\nY8Zw/ahRP0yM4PuEqLDwsPV4zeqJD1Xdo6pfRv6dBqQDnm89X5YcVVUy/BzySkvZXZJCy5IM2jRu\nE7d6vGT11N56LDmKExEZBfwLuBfoA6wCZopI9T4VPJAsyVFhISxdCk884ZKhrl3dwOqxY2H9erj5\nZvjsM7dNwnPPwaWXuvFFsdbjhbJ6VJWtW7bw8rvvMva55zj63XfptX8/D7VqRbf8fF7PyCDt7LN5\n7pZbGDB4MFLZIjTRLUfl6ok3qyf+RKQfUEdVd3l53NJS14oatORodW4uKkKvJhXcKHhYj5esntpb\nj3Wrxc844GlVfRlARG4FLgXGAv/wM7BEKy52aw5t3OiSnjVr4Msv3VTkkhKoVw9694YrroDBg90m\nsK1b+x01lJSWsvnrr1mzeTMbd+zg8uefZ3GrVqS1aAEtWtArN5dr0tMZ3qgRQ4cMoWHTprFVUEFy\nZPwjIoOBu4F+QEdgpKq+V67M7cBdQAfcDc8dqrqkgmO1Al4CbvY6zq1b3dIVJwZsgenlOTkQLmFw\nmyRbW8OYClhyFAciUg/3Aftg2WuqqiIyC6gxy8KqQl4epKdDWpp77N4Nu3a5ne6HDXNJ0fbt7m4X\noFkzN/14wAD4+c+hTx+XGCV6yI2qklVSwr7iYr4tKmL3xo3s+PprdhQWslmEzY0a8U3r1hTXqweN\nGlG/Xj2OUeVn6ekMqFuXs/v1o9U551QviLLkqKCg2udjPNEEWAk8B7xd/s2o1uBbgMW4G6CZInK8\nqqZHlasPvAM8qKpfeB3kunXuaywtR8ng88w0yNtG/x6n+h2KMYdlyVF8tAFSgLRyr6cBPSv7xulr\nc1lfmnPwi1ruqZb7WvZvdV9Vv39NwxBWCJe6Ac7hMGzeW8L4aTmUlHDQo6jYTZ4qKnRfCwqgoNCt\np5JetI0DRfnk50FuPuQdgNxc930HEWjRAoo0i46nLeakYXBOe2jXHjp2UFq0dGU0Evw6lLU73fMw\n7mspYVTdnmOlGqZUlRINRx5KSThMkZZSrMrqrD3csuQdCsJhCsOlFITDFISV/HCYQoXOLXqQUq8p\nOSUlhEpLCZWUsL+4mIySEkr04P/YZi1bcnRmJt3y87kwO5vjSko4sV07TjruOG7p0IH3bva4EcBa\njpKKqs4AZgCIVDjUuaqtwS8Bs1X1tXjEuX69u8no3DkeR4+fJaFMtzJ2h/P9DsWYw7LkKHk0BPjf\njTMhf318a8rfxa93H6KvViKRNASOcJ22LICdIRafsNK9UArsjjy8FirkmRX7XPZHCYRLQEsQLYbS\nIpazgT7tTqRZSl3apaTQvW5dmqekkFq3Li0jjzb16tFGhKbNmlU4N3p3ZiahUIjly5d7G/uWLe7r\nqlUuG4X41FOBmlzP+vXfXT+etUdWtTVYRAYBVwNfisgVuPuUn6jq2goO27BcvFXy2WdubN6KFVX/\nHr9/3mFVvlm7ifp71pDeI539sj8u9XjN6glOPV5f96Ll7p5N9UU+SPOAH0ePWRCRF4HmqnpFBd9z\nHTAxYUEaU/ONOdLWGxEJEzXmSEQ6AruAgdFdZSLyd2CIqsbcXW7XvDFxccTXfTRrOYoDVS0WkWXA\nMKDsw1Uiz/99iG+bCYwBtgI2CMWYI9cQ6Iq7ppKZXfPGeMfT696So/h5BHgxkiSVDd5sDLxYUWFV\n3Q/EZYyCMbXQ5x4fLx3XQdy+3OvtgT1HckC75o3xnGfXva1zFCeq+gZuyu/9wArgVOBCVd3na2DG\nmJipajFQ1hoMHNQa7HUiZozxmY05MsYYQESaAD1w0xKWA78GPgEyVHWHiFyDa/m9le9bg68CTrCb\nHmNqFkuOjDEGEJGhuGSo/IfiS6o6NlLmNuA3uO60lbhFIJcmNFBjTNxZcmSMMcYYE01V7RGnBzAY\nN1ttFxAGLq+gzP24FYDygI+BHuXebwA8gRsQmgNMAdrFUg/wQuT16Mf0WOoBfo/rSsjGLWb5DnC8\n1+dTlXo8Op9bcds/hCKPz4GL4vCzqbQeL87lEL97v4sc6xGvz+lw9Xj087m3gmOsi+e5BOm6P1wd\nXv1eYdd9YK577JqP+VwqvY7j/UFRmx/ARZEf5o9wM13Kf4D9FsgALgNOAd4FvgHqR5V5CjfVdyhu\nA9vPgc9irOcF4AOgLdAu8mherkyl9QDTgZ8AJwK9gGmR8o28PJ8q1uPF+Vwa+X87FjfO5C9AIXCi\nxz+bw9VT7XOp4PfudGAzbiLAI17+fKpYjxc/n3uBL8sdo1W8ziVo130V6vDk9wq77gNx3WPXvOfX\nfFw/JOxx0C9GRXd3u4FxUc9TgXzgmqjnhcAVUWV6Ro41IIZ6XgDeriS2I6mnTeT9s+N8PhXV4/n5\nRMrsB26K17kcoh5PzwVoCmwEzsONn4n+APPsnA5TT7XPCfdBubySY8Tt5+PlgwRc94eoI17XiF33\nSXbdY9d8XK55m8rvExHphtvZe3bZa6qaDXzB99sR9MetRRVdZiOwndg3sD1HRNJEZIOIPBnZNbxM\nvyOopwVu4GpGnM/noHricT4iUkdErsWtQ/V5vM6lfD3xOBdck/L7qjqnXN1en1OF9Xh8TseJyC4R\n+UZEXhWRLnE6l4RJcOxeX/Ng130yXvd2zcd+Lodli0D6pwPu4q9oc9qyDb7aA0WRX4JDlamKD4G3\ngC24Zt6/AdNFZKC69LpDLPVE1nd5FJivquvidT6HqMez8xGRU4CFuJVVc3B3HBtFZKCX53Koerw8\nl0g91wK9cR8S5Xn28zlMPV6d0yLgp7g71Y7AfcC8yP9lIq8dryUqdk+vebDrPtZzScR1b9d8/K55\nS45qAXULUpZZKyKrcX215+CaR2P1JHASMKj60cVej4fnswE4DWiOW6/mZREZUp2AY6lHVTd4dS4i\nchTuD8r56hYsjIuq1OPFOalq9BYAa0RkMbANuAb3/2kqEYdrHuy6j1Vcr3u75uPLutX8swe32Fxl\n2xHsAeqLSGolZWKmqltwo/l7xFqPiDwOXAKco6rfRr3l6flUUo9n56OqJaq6WVVXqOo9uNkld3p9\nLpXU49m54Jqt2wLLRaRYRIpxgxLvFJEi3N2TF+dUaT2Ru36vzin6GCHgq8gxfLt2POBL7NX9Gdh1\nn5TXvV3zcbxuLDnySeSXZw8Hb0eQCpzB9/3Sy4CScmV6AkfjmmuPSOROoDVQ9uFTpXoiH1w/As5V\n1e3xOp/K6vHyfCpQB2iQgJ9NHdyUUy/PZRZulk9v3N3qacBS4FXgNFXd7NE5Ha4e9fCcoo/RFPch\nudvPa6e6/Iq9Oj8Du+6T9rq3az6e143GOILbHlV/AE1wv0i9cSPmfxV53iXy/m9wMxhG4H753gU2\ncfDUxCdx/bjn4DL4Bfxw2ugh64m894/IL9ExkV+cpcB6oF5V64m8n4lbX6V91KNhVJlqn8/h6vHw\nfB6M1HEMblro33AX1nke/2wOWY9X51LJ71/5GSWenFNl9Xj48/knMCRyjLNwa5qkAa3jeS5Bue4r\nq8PL3yvsug/UdY9d855d874lDrXhgWt6DOPWIYl+PB9V5j6+X9RqJhUvavV/fL+o1Zs4oYd+AAAD\nMklEQVT8cMGxQ9aDGww4A5d1F+DWqHgKaBtLPYc4filwQ7njVOt8DlePh+fzbOR78yPH+ojIB6TH\nP5tD1uPVuVTy+zeHHy4IV+1zqqweD38+k4Cdkf+37bjd67vF+1yCct1XVoeXv1eHqMOu+yS97rFr\n3rNr3rYPMcYYY4yJYmOOjDHGGGOiWHJkjDHGGBPFkiNjjDHGmCiWHBljjDHGRLHkyBhjjDEmiiVH\nxhhjjDFRLDkyxhhjjIliyZExxhhjTBRLjowxxhhjolhyZEwViUhPEflWRJrE8D0XisiKeMZljIkf\nu+5rJ0uOTK0iImERKY18Lf8oFZE/VfLtDwKPqWpu5FhDI9+XGnX8TiKyWkQ+FZFmqjoTKBKRMXE+\nNWPMIdh1b2JlyZGpbToAHSNffwWEcLt/l73+cEXfJCJHA5cCL5V7S6PKHAt8htt48UJVzYm89RJw\np3enYIyJkV33JiaWHJlaRVX3lj1wH5CqqvuiXs87xLdeDaxS1W8relNETsV9QC4ArlDVwqi33wf6\ni0g3D0/FGFNFdt2bWFlyZEzVDAaWHuK9QcCnwJuqeoOqhqPfVNUdQFrkGMaY4LDrvpay5MiYqjkG\n2F3B6wK8A7ynqpU1oe+OHMMYExx23ddSlhwZUzWNgIJDvPcucIWInF3J9+cDjT2PyhgTT3bd11KW\nHBlTNelAywpeV+DnwOvAhyJyqCb0VsC+OMVmjIkPu+5rqbp+B2BMQKwATjrEe6qqt4hIGPhARC5T\n1Xllb4pIA+DYyDGMMcFh130tZS1HxlTNTGCgiEi51797rqq3Aq8A00VkaFSZgbim+YVxj9IY4yW7\n7mspS46MqZoPgRLg/HKv60FPVG8HXgCmRX1QXgtMVNVDjV0wxiQnu+5rKVHVw5cyxiAitwEjVPXi\nGL6nNbAB6K+q2+IWnDEmLuy6r51szJExVfc00FxEmpRtJVAFXYHb7APSmMCy674WspYjY4wxxpgo\nNubIGGOMMSaKJUfGGGOMMVEsOTLGGGOMiWLJkTHGGGNMFEuOjDHGGGOiWHJkjDHGGBPFkiNjjDHG\nmCiWHBljjDHGRLHkyBhjjDEmyv8DsOmLti1qNdYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fd1dc59c4e0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "#activate a figure 1\n",
    "plt.figure(1)\n",
    "plt.subplot(121) #activate a subplot\n",
    "plt.title('Saturation pressure, hexane')\n",
    "\n",
    "#plot the calculated Psat x T\n",
    "\n",
    "#filter not-a-number and zeros\n",
    "grid_P_real = ~np.isnan(grid_P)\n",
    "filtered_grid_T = grid_T[ np.where( grid_P_real ) ]\n",
    "filtered_grid_P = grid_P[ np.where( grid_P_real ) ]\n",
    "\n",
    "#plot\n",
    "\n",
    "plt.plot(filtered_grid_T,filtered_grid_P)\n",
    "\n",
    "#calculate Antoine correlations\n",
    "\n",
    "## nist\n",
    "grid_TAntoine_nist_1 = np.linspace(177.7,264.93,100)\n",
    "grid_PAntoine_nist_1 = np.array([ PAntoine_nist(grid_TAntoine_nist_1[i]) for i in range(100) ])\n",
    "\n",
    "grid_TAntoine_nist_2 = np.linspace(286.18,342.69,100)\n",
    "grid_PAntoine_nist_2 = np.array([ PAntoine_nist(grid_TAntoine_nist_2[i]) for i in range(100) ])\n",
    "\n",
    "# plot\n",
    "plt.plot(grid_TAntoine_nist_1,grid_PAntoine_nist_1)\n",
    "plt.plot(grid_TAntoine_nist_2,grid_PAntoine_nist_2)\n",
    "\n",
    "## ddbst\n",
    "grid_TAntoine_ddbst = np.linspace(100,Tc,100)\n",
    "grid_PAntoine_ddbst = np.array([ PAntoine_ddbst(grid_TAntoine_ddbst[i]) for i in range(100) ])\n",
    "\n",
    "# plot\n",
    "\n",
    "plt.plot(grid_TAntoine_ddbst,grid_PAntoine_ddbst)\n",
    "\n",
    "# aditional options\n",
    "\n",
    "## y and x axis limits\n",
    "plt.xlabel('T (K)')\n",
    "plt.ylabel('P (Pa)')\n",
    "plt.xlim([100,Tc])\n",
    "plt.ylim([0,Pc])\n",
    "\n",
    "#activate a subplot\n",
    "plt.figure(1)\n",
    "plt.subplot(122) #activate a subplot\n",
    "\n",
    "plt.title('(log scale)')\n",
    "\n",
    "# same plot, now in log scale\n",
    "plt.semilogy(filtered_grid_T,filtered_grid_P)\n",
    "plt.semilogy(grid_TAntoine_nist_1,grid_PAntoine_nist_1)\n",
    "plt.semilogy(grid_TAntoine_nist_2,grid_PAntoine_nist_2)\n",
    "plt.semilogy(grid_TAntoine_ddbst,grid_PAntoine_ddbst)\n",
    "plt.xlabel('T (K)')\n",
    "plt.xlim([100,Tc])\n",
    "plt.ylim([100,Pc])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "In this notebook, a simple cubic equation of state was used to calculate the volume, the residual Gibbs energy and the saturation pressure of a pure fluid. Here, the saturation pressure was obtained by finding a pressure that would lead to equal residual Gibbs energy between two phases. In order to achieve that, this notebook introduced key concepts of phase equilibrium calculation (equality of molar Gibbs energy, for instance) and major tools for applied thermodynamics, such as root finding numerical methods available in scipy. After studying this notebook, the student should be able to solve simple fluid phase equilibrium problems and use numerical methods already implemented in Python libraries. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## External references\n",
    "\n",
    "### On Python\n",
    "* Guidelines for importing functions from Scipy: [in Scipy / API](https://docs.scipy.org/doc/scipy/reference/api.html#guidelines-for-importing-functions-from-scipy)\n",
    "* More opinion on Scipy, NumPy and matplotlib [by josephcslater](http://josephcslater.github.io/scipy-numpy-matplotlib-pylab.html)\n",
    "* Scipy, Optimize: [in Scipy / optimize](http://docs.scipy.org/doc/scipy/reference/optimize.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Credits\n",
    "* Initially developed in python 2.7 by Guilherme Carneiro Queiroz da Silva\n",
    "* Adapted to ipynb with python 3.5 by Iuri Soter Viana Segtovich\n",
    "* Review - Fernando de Azevedo Medeiros"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}