{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Natural gas hydrates - Phase Equilibria\n",
    "Liquid(Water)-Solid(Hydrate-sI)-Gas(Methane) Equilibria\n",
    "> cf. Segtovich et al., 2016, supplementary material "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Overview"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Equilibrium criteria\n",
    "\n",
    "Equality of pressure, temperature and chemical potential among all phases for each component \n",
    "\n",
    "$$\\mu_w^L=\\mu_w^H$$\n",
    "$$\\mu_w^G=\\mu_w^H$$\n",
    "$$\\mu_g^L=\\mu_g^H$$\n",
    "$$\\mu_g^G=\\mu_g^H$$\n",
    "\n",
    "Note that the solubility of water in the gas phase is known to be small ( circa 100 ppm ) and does not significantly affect the chemical potential of other components in the gas phase\n",
    "Similarly, the solubility of guest components in the liquid water phase is known to be small ( even lower than 1 ppm ) and does not significantly affect the chemical potential of water in that liquid water phase.\n",
    "\n",
    "Therefore we are able to split the problem in two sections with .\n",
    "\n",
    "First we can calculate phase equilibria using equality of chemical componets present in significant ammount in each phase\n",
    "\n",
    "$$\\mu_w^L=\\mu_w^H$$\n",
    "\n",
    "$$\\mu_g^G=\\mu_g^H$$\n",
    "\n",
    "Then we may calculate the solubility of guest components in water and of water in gas phase at equilibrium condition, if we want to do so."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Models\n",
    "\n",
    "Regarding models for said chemical potentials in each phase,\n",
    "\n",
    "We know an equation of state for the gas phase, from which we are able to calculate fugacity as function of temperature, of either density or pressure, and of composition.\n",
    "\n",
    "$$f_g^G \\leftarrow f_g^G(T,P,x_g^G)$$\n",
    "\n",
    "We know an equation of state for the adsorption of guest components in an hypothesized empty lattice phase, from which we are able to calculate chemical potential difference for water between filled and empty lattice conditions, as well as composition of guests in the lattice, as function of temperature, and fugacity of the guests in the hydrate phase.\n",
    "\n",
    "$$\\Delta\\mu_w^{H-EL} \\leftarrow \\Delta\\mu_w^{H-EL}(T,f_g^H)$$\n",
    "\n",
    "Regarding the liquid water phase and the hypothesized empty lattice phase, we can calculate a difference in the cehmical potential of water between these phases from classical thermodynamics integration as a function of temperature and pressure, given a initial condition of T0, P0 and a reference chemical potential difference.\n",
    "\n",
    "$$\\Delta\\mu_w^{Lw-EL} \\leftarrow \\Delta\\mu_w^{Lw-EL}(T,P)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Devised algorithm\n",
    "\n",
    "Having hold of these models we can devise an algorithm to achieve equality of chemical potential in the following \n",
    "\n",
    "* $\\mu_g^G=\\mu_g^H$ is met by means o meeting $f_g^G=f_g^H$\n",
    "\n",
    "* $\\mu_w^{Lw}=\\mu_g^H$ is met by means of meeting $\\mu_w^{Lw}-\\mu_w^{EL}=\\mu_g^H-\\mu_w^{EL}$ which we name $\\Delta\\mu_w^{Lw-EL}=\\Delta\\mu_g^{H-EL}$\n",
    "\n",
    "Then,\n",
    "\n",
    "Equality of guest fugacity is met by **direct substitution** of gas fugacity into hydrate model:\n",
    "\n",
    "$$f_g^H \\impliedby f_g^G(P,T,x_g^G)$$\n",
    "\n",
    "making the functionality of $\\Delta\\mu_w^{H-EL}(T,P,f_g^H)$ become $\\Delta\\mu_w^{H-EL}(T,P,f_g^H(T,P,x_g^G))$ therefore $\\Delta\\mu_w^{H-EL}(T,P,x_g^G)$.\n",
    "\n",
    "For given $P$ and $x_g^G$ , the equality of the chemical potential differences for water is met by numerical solution of this equality for $T$.\n",
    "\n",
    "$T^{SOLUTION} \\leftarrow SOLVE \\left({\\Delta\\mu_w^{H-EL}(T<^{SOLVE},P<^{GIVEN},x_g^G<^{GIVEN})=\\Delta\\mu_w^{Lw-EL}(T<^{SOLVE},P<^{GIVEN})}\\right)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Implementation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Implement model for calculation of fugacity in gas phase\n",
    "e.g. Peng-Robinson [REF](DOI) (Poling, Abreu Lecture Notes, PR original)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def PR_fug_G(T,P,x): #Peng-Robinson\n",
    "    R = 8.31446 #m^3*Pa/(mol*K)  \n",
    "    \n",
    "    Ncomp=1 #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<INPUT\n",
    "\n",
    "    k=np.zeros([Ncomp,Ncomp])\n",
    "    Tc=np.ones([Ncomp])\n",
    "    Pc=np.ones([Ncomp])\n",
    "    w=np.ones([Ncomp])\n",
    "    \n",
    "    sigma=1+np.sqrt(2)\n",
    "    epsilon=1-np.sqrt(2)\n",
    "    \n",
    "    #e.g.METHANE\n",
    "    Tc[0]=-82.59+273#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<INPUT\n",
    "    Pc[0]=45.99*1e5#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<INPUT\n",
    "    w[0]=0.010#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<INPUT\n",
    "    \n",
    "    ai_crit = (0.45724*((R)**2)*(Tc[:])**2)/(Pc[:])\n",
    "    bi_crit = (0.07780*(R)*(Tc[:]))/(Pc[:])\n",
    "    \n",
    "    PRk = 0.37464 + 1.54226*w[:]-0.26992*(w[:])**2\n",
    "    \n",
    "    alphaAi= (1. + PRk*(1.-np.sqrt(T/Tc[:])))**2\n",
    "    \n",
    "    ThetaAi=ai_crit[:]*alphaAi[:]\n",
    "    \n",
    "    b_m=sum(x*bi_crit)\n",
    "    \n",
    "    Theta_m=0\n",
    "    for i in range(Ncomp):\n",
    "        for j in range(Ncomp):\n",
    "            Theta_m += x[i]*x[j]*np.sqrt(ThetaAi[i]*ThetaAi[j])*(1.- k[i,j])\n",
    "            \n",
    "    Theta_m=sum(sum(np.array( #sums over all elements of the following temporary 2d array\n",
    "        [x[i]*x[j]*np.sqrt(ThetaAi[i]*ThetaAi[j])*(1.- k[i,j])\n",
    "        for (i,j) in np.ndindex(Ncomp,Ncomp)]).reshape((Ncomp,Ncomp))        \n",
    "    ))            \n",
    "    \n",
    "    c0 = -(b_m**3)*sigma*epsilon + (-R*T*sigma*epsilon*(b_m**2)-b_m*Theta_m)/P\n",
    "    c1 = (b_m**2)*(sigma*epsilon-epsilon-sigma) + ((R*T)*(-sigma*b_m-epsilon*b_m) + Theta_m)/P\n",
    "    c2 = epsilon*b_m+sigma*b_m-b_m-R*T/P\n",
    "    c3 = 1\n",
    "        \n",
    "    V = np.roots(np.array([c3,c2,c1,c0]))\n",
    "\n",
    "    \n",
    "    V = V[ np.isreal(V) ].real #filter out complex roots\n",
    "    V = V[ np.where( V > b_m ) ] #filter out non physical roots\n",
    "\n",
    "    V = max(V) #GASphase-only EOS implementation  \n",
    "    \n",
    "#    print(\"debug V=\",V*1000)\n",
    "    \n",
    "    dbidni=bi_crit[:]\n",
    "\n",
    "    dthetaidni=np.zeros(Ncomp)\n",
    "    for i in range(Ncomp):\n",
    "        sum1 = 0.\n",
    "        for j in range(Ncomp):\n",
    "            sum1 += x[j]*np.sqrt(ThetaAi[j])*(1.-k[i,j])\n",
    "\n",
    "        dthetaidni[i]=np.sqrt(ThetaAi[i])*sum1\n",
    "\n",
    "    qsi = (1./(b_m*(epsilon-sigma)))*np.log((V+epsilon*b_m)/(V+sigma*b_m))\n",
    "    \n",
    "    lnPhi = (((dbidni[:]/b_m)*((P*V)/(R*T)-1) \n",
    "                -np.log(P*(V-b_m)/(R*T)) \n",
    "                -(Theta_m/(R*T))*qsi)*\n",
    "                ((2*dthetaidni[:]/Theta_m)- \n",
    "                (dbidni[:]/b_m)))\n",
    "\n",
    "    fug = np.exp(lnPhi)*x*P\n",
    "\n",
    "    return fug, V\n",
    "\n",
    "def test_PR():\n",
    "    R = 8.31446 #m^3*Pa/(mol*K)  \n",
    "    tryT=273\n",
    "    tryP=1e5\n",
    "    tryx=np.array([1.])\n",
    "    ans=PR_fug_G(tryT,tryP,tryx)\n",
    "    print(\"Z\",tryP*ans[1]/(R*tryT))\n",
    "    print(\"f\",ans[0])\n",
    "\n",
    "    tryT=273\n",
    "    tryP=1e9\n",
    "    tryx=np.array([1.])\n",
    "    ans=PR_fug_G(tryT,tryP,tryx)\n",
    "    print(\"Z\",tryP*ans[1]/(R*tryT))\n",
    "    print(\"f\",ans[0])\n",
    "\n",
    "    tryT=298\n",
    "    tryP=1e5\n",
    "    tryx=np.array([1.])\n",
    "    ans=PR_fug_G(tryT,tryP,tryx)\n",
    "    print(\"Z\",tryP*ans[1]/(R*tryT))\n",
    "    print(\"f\",ans[0])\n",
    "    \n",
    "    tryT=298\n",
    "    tryP=1e9\n",
    "    tryx=np.array([1.])\n",
    "    ans=PR_fug_G(tryT,tryP,tryx)\n",
    "    print(\"Z\",tryP*ans[1]/(R*tryT))\n",
    "    print(\"f\",ans[0])\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Z 0.997061903565\n",
      "f [ 99706.49713971]\n",
      "Z 12.685200677\n",
      "f [  1.71815606e+13]\n",
      "Z 0.997781120382\n",
      "f [ 99778.1888341]\n",
      "Z 11.7011817254\n",
      "f [  8.32866536e+12]\n"
     ]
    }
   ],
   "source": [
    "#test\n",
    "test_PR()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The model for $\\Delta\\mu_w^{H-EL}$, van der Waals and Platteeuw, 1959\n",
    "langmuir constants as function of temperature\n",
    "$$C_{i,j}=\\frac{A_{i,j}}{T}exp\\frac{B_{i,j}}{T}$$\n",
    "\n",
    "composition as occupancy of cavities as function of langmuir constants and fugacities\n",
    "$$\\Theta_{i,j}=\\frac{C_{i,j}f_{i}}{1+\\sum_k\\left[{C_{k,j}f_{k}}\\right]}$$\n",
    "\n",
    "chemical potential difference as function of occupancy\n",
    "$$\\Delta\\mu_w^{H-EL}=RT\\left( \\sum_j\\left[{\\nu_j ln\\left(1- \\sum_i\\left[{\\Theta_{i,j}}\\right]\\right)}\\right] \\right)$$\n",
    "\n",
    "where A and B are empirical parameters presented by Parrish and Prausnitz (...), Munck (...), Sloan (...), Segtovich, 2016."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def vdwnp_H_EL(T,fug):\n",
    "    \n",
    "#parameters\n",
    "\n",
    "    R = 8.31446 #!m^3*Pa/(mol*K)  \n",
    "\n",
    "    Nguest=1    #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<INPUT components\n",
    "    Ncav=2 #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<INPUT structure\n",
    "\n",
    "    UCnw = 46 #sI #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<INPUT structure\n",
    "    ni = np.array([2/UCnw, 6/UCnw]) #sI #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<INPUT structure\n",
    "\n",
    "    Nguest = 1\n",
    "    A_pSWP = np.zeros([Nguest,Ncav])\n",
    "    B_pSWP = np.zeros([Nguest,Ncav])\n",
    "    \n",
    "    #methane\n",
    "    A_pSWP[0,0] = .7228e-8 #Pa^-1.K #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<INPUT components,structure\n",
    "    B_pSWP[0,0] = 3187. #K #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<INPUT components,structure\n",
    "    A_pSWP[0,1] = 23.35e-8 #Pa^-1.K #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<INPUT components,structure\n",
    "    B_pSWP[0,1] = 2653. #K #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<INPU2*T components,structure\n",
    "\n",
    "#LANGMUIR CONSTANTS\n",
    "    Clang=np.zeros([Nguest,Ncav])\n",
    "    for i in range(Nguest):\n",
    "        for j in range(Ncav):\n",
    "            Clang[i,j] = (A_pSWP[i,j]/T)*np.exp(B_pSWP[i,j]/T) #Pa\n",
    "\n",
    "#OCCUPANCY\n",
    "    occ=np.zeros([Nguest,Ncav])\n",
    "    for j in range(Ncav):\n",
    "        sum1 = 0.      \n",
    "        for i in range(Nguest):\n",
    "            sum1 += Clang[i,j]*fug[i]\n",
    "\n",
    "        for i in range(Nguest):\n",
    "            occ[i,j] = Clang[i,j]*fug[i]/(1.+sum1)         \n",
    "#vdwnp\n",
    "    Dmu_H_EL = 0.\n",
    "    for j in range(Ncav):\n",
    "        sum1 = 0.\n",
    "        for i in range(Nguest):\n",
    "            sum1 += occ[i,j]\n",
    "            \n",
    "        Dmu_H_EL += ni[j]*np.log(1.-sum1)\n",
    "\n",
    "    return Dmu_H_EL\n",
    "\n",
    "# Here there is a standard test for vwnp\n",
    "def test_vwnp():\n",
    "    tryT=273\n",
    "    tryELV=1e-5\n",
    "    tryf=np.array([100000])\n",
    "    print(vdwnp_H_EL(tryT,tryf))\n",
    "    \n",
    "    tryT=273\n",
    "    tryf=np.array([1e9])\n",
    "    print(vdwnp_H_EL(tryT,tryf))\n",
    "    \n",
    "    tryT=298\n",
    "    tryELV=1e-5\n",
    "    tryf=np.array([100000])\n",
    "    print(vdwnp_H_EL(tryT,tryf))\n",
    "    \n",
    "    tryT=298\n",
    "    tryf=np.array([1e9])\n",
    "    print(vdwnp_H_EL(tryT,tryf))\n",
    "   \n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.127099254882\n",
      "-1.59686925346\n",
      "-0.0637573948653\n",
      "-1.43275104482\n"
     ]
    }
   ],
   "source": [
    "#test\n",
    "test_vwnp()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The model for $\\Delta\\mu_w^{EL-Lw}$, Saito (...)\n",
    "Later reviewed by Parrish and Prausnitz (...), Munck (...), Sloan (...), Segtovich, 2016.\n",
    "$$\\Delta\\mu_w^{EL-Lw}=RT\\left( \\frac{\\Delta\\mu_w^{EL-Lw}|_{T0,P0}}{RT_0}-\\int_{T0}^{T} \\frac{\\Delta h_w^{EL-Lw}|_{P0}}{RT^2}+ \\int_{P0}^{P} \\frac{\\Delta v_w^{EL-Lw}|_{T}}{RT} \\right)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def saito_EL_Lw(T,P):\n",
    "    \n",
    "    R = 8.31446 #!m^3*Pa/(mol*K)  \n",
    "    \n",
    "    P0=612.62 #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<INPUT\n",
    "    T0=273.18 #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<INPUT\n",
    "    \n",
    "    Dmu00 = 1264.00 #EL-LW  #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<INPUT structure\n",
    "    Dh00  = -4858.00 #EL-LW   #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<INPUT structure\n",
    "    Dcp00 = -39.16 #EL-LW?   #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<INPUT structure\n",
    "\n",
    "    \n",
    "    #-S_DH_EL_LW_RT2_dT (T0->T, P)\n",
    "    GenVantHoffFactor = (1./R) * (((1./T)-(1/T0))*(Dh00-Dcp00*T0) - (Dcp00)*np.log(T*(1/T0)))\n",
    "       \n",
    "    #VLW\n",
    "    V_Lw = (18.)*1e-6 # ( ( 1g/cm3 / 18g/mol ) ^-1 ) * 10^-6 m3/cm3 = ___ m3/mol\n",
    "\n",
    "    #VEL\n",
    "    \n",
    "    D_V_EL_Lw = 4.6e-6 #*1.2   #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<INPUT structure\n",
    "    # Try multiplying this parameter for 1.2 and see the effect in the phase diagram\n",
    "    \n",
    "    V_EL = V_Lw+D_V_EL_Lw\n",
    "    \n",
    "    #Poyinting_EL_Lw (T, P0->P)\n",
    "    Poynting_EL_Lw = D_V_EL_Lw*(P-P0)/(R*T)\n",
    "    \n",
    "    Dmu_EL_Lw = ( Dmu00/(R*T0) + GenVantHoffFactor + Poynting_EL_Lw ) #EL-PW\n",
    "    \n",
    "    return Dmu_EL_Lw\n",
    "\n",
    "# Here there is a standar test\n",
    "def test_saito():\n",
    "    tryT=273\n",
    "    tryP=1e5\n",
    "    print(saito_EL_Lw(tryT,tryP))\n",
    "\n",
    "    tryT=273\n",
    "    tryP=1e9\n",
    "    print(saito_EL_Lw(tryT,tryP))\n",
    "    \n",
    "    tryT=298\n",
    "    tryP=1e5\n",
    "    print(saito_EL_Lw(tryT,tryP))\n",
    "\n",
    "    tryT=198\n",
    "    tryP=1e9\n",
    "    print(saito_EL_Lw(tryT,tryP))  \n",
    "    \n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.555290889116\n",
      "2.58165602505\n",
      "0.752126561468\n",
      "2.81098609671\n"
     ]
    }
   ],
   "source": [
    "#test\n",
    "test_saito()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAE0CAYAAADQYm9sAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcnfPd//HXO7JW5G5paotkqJKlQrWTWptBiK1oVDVE\nbxrCzR1d9BetqoTW3U25g6KISLVZai1aSm4mKJpULRGJpdkEibFFyDbi8/vjukZOxizXzJyZOefM\n+/l4nEfO9b2W8/3OZOYz310RgZmZWRad2jsDZmZWPBw0zMwsMwcNMzPLzEHDzMwyc9AwM7PMHDTM\nzCwzBw0zM8vMQcPMzDJz0DBLSfq+pP+RVFI/F5LOknRFe+fDSkNJ/XCYtVDPiDgvIj4EkLS3pOsk\n/UTSOZIOlHS5pK9KekjSK5LOT6+9ID2eJenIpnxoPp4n6VBJv5Q0WtJ5kr5bcy4ifgtUNekrYVaP\nzu2dAbNCJGk48P+AIyNibZp2AHBaRJwtaUvgiIj4GUBEXCRpN+CuiLi7KZ8VEXe15HmSjgBOBr4R\n6bpAaa3psoj4XlPyYtYYBw2zjQQgqQtwHXBUTcAAiIgHJVW2U97qJGkz4DLg0Nh0IbnLgSWSrouI\n50jLZtZSbp4y26jml+6+QO+IeKqOa25t6YdImitpz5Y+J7UPsFVELMxNjIgPgGeAEXn6HDPANQ2z\numwLvFnXiYi4PudwV0nj0vcCdsn4/J8Az9eR3pznbQ+8V8+5t9PzZnnjoGH2ca8An67rhKQtI+Kt\n9PD5iPhVzrkv1XPPbyLinJrjiLijns9t8HmSdgZ2S193R8S/gFeBLet53jbAc/WcM2sWN0+ZbVTT\n7v8osFzSF+u45tgmPVDaCdi9pRlLfZUkoF0G/CBNexRYJam81uf2BL4A3J6nzzYDHDTMcgV81B9w\nGnCppE/WnJS0LbCqic8sA5bmJkg6Lu3AblrmIi6LiNlAH2BRTl7/G7iw1uXfA26IiHk1tzf188zq\n4uYpszpExP2S3iUJHItIgsUbEfEHSYeTDHHdSdK5EfFLST8EhgBbSaqKiHsk7QXMBk6q9fiLgYXA\nEwBZn5dz/zHpM2ryepuk9yX9CniJJKhURcRP8/xlMUPe7tUsIemCiLgoj887HlgLnA38NCIq0/Te\nwA5pn0RTn/lVoBLYJiJebMJ9eS2bdVxunjJrJRExA5gLdAV65Jzam2Q4bJNI+hrJyKtbgW/kI49m\nTeXmKbON8j4BLp0/sX+t5GVpX0RTn3U7ze/Y9uQ+ywvXNMw2WtUWCxY2p1mqJSSdRT1DiM2ayn0a\nZmaWmWsaZmaWmYOGmZll5qBhZmaZOWiYmVlmDhpmZpaZg4aZmWXmoGHWiiQ9K+krxfLctpZbDkmL\nJB1YO732uaY+1/LLQaMDkvRVSQ9JekXS+WnaBenxLElH1nHPcZIelrRM0gV5+OxXJV0saeuWlCXj\nZ54m6eV0x7yDctLHS3pL0th67vswXdqcWvfclPWzI+LzEfFQ83Pf+s9t6i/kfKqvHC0tX2t93c3L\niHRIEXGXpC2BIyLiZ2naRZJ2A+6KiLvruOdmSZ9I72n2wne1PvvHzX1OEz/zOkkfAidFxP/lnLoa\nWBoRk+u7tYnprU7SZhGxob0+v9D569P6XNOwjuJWYIikPjlpRwM3N3BPo+s1pX+lnyPpaUlvS5om\nqWut8wdKGifp5lr3TpT0v+n7cyW9JOndtGnlmFrPGCfpaeA9SZvVrh1kuD83j9Nr8ijp90Bf4K70\n3h9Qi6RtJd0i6XVJ/86tmUn6gqQnJK1MnztN0kU55zeprUmaXOt8nbWcetKHSJon6U1JN+SUoUVf\nH2saBw2rk6TBkk6V9DVJt7XRZ46VtE7SKUrWgKqStJ+kEyW9Jqn2wn+ZRcQ7wD3ACTnJvSKivv21\nm+I44BBgR5Jd+k6u45rpwGGSNgdQsr7VccAf0/MvAftGRC+SDZX+UKvp7pvAYcAn6/lLurH7c/M4\nuCaPEfEtkk2ijoyIXhFxSe5DJQm4C3iSZO/0g4DvSDpYUheSBRSnkGw5ezMf39kwn7WyE4CDgc+S\n7J9+fs65ln59LCMHjY5t1/QvtHGSzgV2zTl3GnBvurJqfXta51VEXAHcB8yJiPNINj56DLgN+F5E\nPNzCj5gKnAgg6bNA5v0oGjExIlakgekuYI/aF0TEUuBfwNfSpIOA9yNiTnr+1ohYkb6/Oc3bkFqf\n8WpErKsrAxnvbyiP9dWqyoFPR8TFEbEhIhYD1wMjgb2AzhFxeXruVmBOxuc2xxXp1+Adkk2oRuac\na+nXxzJyn0bH9nxE/KrmQNKXcs7dDjwh6SGSPakbJGkc0L12MslfmlMiYkkj9/ePiAXAn4FjJHUH\n/knyy7UrcG8ePvNu4HpJg4ADgGsaKdYGoEuttC5Ada20FTnvV5P8RV6XaSS/6P6Q/jv1o0xL3yLZ\norUsTdqcTVemXdZQRjPcnzWPtfUDtpf0Vs1Hkfyx+TCwHcme5bka/D63UO7XYEn6+XWd+5gMXx/L\nyEHD6vMCMBA4HLiusdE1ucGnmfYDaoLG3SRNJv8DjAEeTf+6bNFnRsQ6SbcCo4DXM+xpsZTkl8zz\nOWk71jpuipuBSyRtT1Lj2AtAUl/gWuCAiHgsTXuSTf9Kr7eZJ+P9DWmoCellYGFE7Fr7hJIhrdvX\nSu5L0hRUYzXwiZzjbdJnNscOOe/7Aa/mHLfm18dyFEXzlKQdJV0v6U/p8Sck3Sjpd5JOaOx+a5ax\nwLsRcRMwkY1/meb9B01SN+BLABFRBawjaZt+CvgCH//LviWmAWeR1GIaMwM4X9L2SgwDjgRuac4H\nR8QbwCxgMskv4prgsznwIfCGpE6STgE+34RHt/T+5cBO9ZybTbLPyDhJ3dNO5kFprfQx4IO0L6qz\npBF8vMnnSeCENF+HAkObkK/azkq/F1sC55H0E2XR0q+P5SiKoBERiyLi1JykEcDNEXE6cFQ7Zato\nSTqcpCP0y2lfBpJ+SPIDf7Kkw0h+cY+WdCLQk+SXyinpPefl4bP3TEe0jAf+zqZ/jd4CPJ6+ryRp\nCsmXB4AnM/aPXAQ8CjwCvAX8AjghIp7Luaaxjt7a56eSNLn98aMLIuYDvyEp83JgUPqZDX3GR2nN\nvD/XL4CfKJmz8v1NPiTiQ5JAuQewCHgduI5kEEE1yc/iKcCbJJ3tt9Z69ndJfkbfJmmSq73zYGR4\nX3M8laTP6yWSPomLGyhfU74+1gTtsgmTpEkk/xFXRMTgnPRDgf8lCWaTIuKXte77U0R8I/0F99eI\neEbSHyPixLbMv5nVTdJk4OWIaPYEUCts7VXTmAwMz01IhyBemaYPAkZK6l/rvpqmkWVAn1ppZmbW\nytolaETEIyTV1VxDgBcjYkla7Z1OMvkKSVtKuhrYI21OuRX4uqTfkgwfNLPC4P2jS1whjZ7ank1H\nVSwj7VSLiLeA/6p1/bcbe6Ak/wc2aweSftLeebCWiYg6W3GKoiO8JSKiIF/jx49v9zy4fC5fRyxf\nKZctX+VrSCEFjVdIxnjX6MPHJw412YQJE6isrGzpY8zMSl5lZSUTJkxo8Jr2DBpi007sOcDOkvql\nC5F9E7izpR8yYcIEKioqWvoYM7OSV1FRUZhBQ9JUkvHvu0haKumUSBYZG0syDnseMD2S8dUlqdQD\nmctX3Eq5fKVcNmj98rXLPI22IinGjx9PRUVFyf9HMTNrqcrKSiorK7nwwguJejrCSz5olHL5zMxa\ng6R6g0YhdYS3CneEm5llk6Uj3DUNMzPbRIeuaZiZWf6UfNBw85SZWTZunnLzlJlZk7l5yszM8sJB\nw8zMMiv5oOE+DTOzbNyn4T4NM7Mmc5+GmZnlhYOGmZllVvJBw30aZmbZuE/DfRpmZk3mPg0zM8sL\nBw0zM8vMQcPMzDJz0DAzs8xKPmh49JSZWTYePeXRU2ZmTebRU2ZmlhcOGmZmlpmDhpmZZeagYWZm\nmTlomJlZZg4aZmaWWckHDc/TMDPLxvM0PE/DzKzJPE/DzMzywkHDzMwyc9AwM7PMHDTMzCwzBw0z\nM8vMQcPMzDIr2qAhaYCkGZJ+K+nY9s6PmVlHULRBAzgMuDwizgK+1d6ZMTPrCAomaEiaJGmFpGdq\npR8qaYGkFySdm3PqJuCbkn4FbNmmmTUz66AKZka4pP2A94DfR8TgNK0T8AJwEPAqMAf4ZkQsyLmv\nE3BrRHytjmd6RriZWRM1NCO8c1tnpj4R8YikfrWShwAvRsQSAEnTgaOBBem15wGfAH7dppk1MysC\nGzbAvHnw97/DFlvAqFEtf2bBBI16bA+8nHO8jCSQkAaS0xt7QO7iWxUVFVRUVOQ1g2ZmheLdd+Ef\n/4BHH00CxT/+AdtsA/vsA1/7WFvMRpWVlZkXdi2Y5imAtPZwV07z1LHA8IgYkx6PAoZExNkZn+fm\nKbMSUlVVxeLFiykrK6N37971pnUUS5fCI48kAeLvf4eXXoI990yCxD77wN57Q3O+JEXRPFWPV4C+\nOcd90rTMJkyY4BqGWQmYNm0Go0efSdeuZaxfv5hJk64C+FjayJHHt3NOW8eGDTB3bhIcagLF2rWw\n336w775w0klJwOjatfmfkaXGUWg1jTKSmsZu6fFmwPMkHeGvAbOBkRExP+PzXNMwKwFVVVX069ef\nNWseBAYDz9CjxwFEfMjatbM2SVuyZEFJ1DjWroXZs+Hhh5Mg8dhjSVPTvvsmgWK//WDnnUF11gda\npihqGpKmAhXAVpKWAuMjYrKkscB9JMODJ2UNGDVc0zArfosXL6Zr1zLWrBmcpgymU6c+wDqSgJGk\ndenSj8WLFxdl0HjnnY21iIcfhiefhIEDYf/9YcwYmDIFPvOZ1s1D0dU08s01DbPSkK+aRiH1f7z+\nehIcHnooeb30EpSXJ0Fi//1hr72gZ8/2yVtR1DTMzOrTu3dvJk26itGjD6BLl35UVy/J6dPYNK2+\nYFBXn0hb9n8sWwazZm0MEq+9ljQ1feUrcNVV8MUvtqw/oq2UfE1j/Pjxbp4yKxHNHT1VX02lNfs/\nlixJgsSsWVBZCStXJgFi6NDk38GDYbPNWuWjm62meerCCy+st6ZR8kGjlMtnZtnMmTOHgw8+g5Ur\nn/gorVevPZk583eUl5czf/58Zs+ezZAhQxgwYECzPmPJEnjwwSRAzJoFq1cnAaLmNXAgdCqYhZsa\n5uYpM+vQysqSJil4hpqaRnX1EsrKyhg79rtceeW1wA7Ay/z3f5/GFVdMbPSZL7+cBIgHH0xeq1fD\nAQdARQWcey707986I5vaW8nXNNw8ZWawsU8jt/9jjz0GM3DgF4HHqQkmsBfPPffEx2ocy5cnweGB\nB5J/V65MAkRNoBgwoPiDhJun3DxlZjlq939cccUVnH32FSTrotb4HJdffjYnnjiWWbOSIPHAA/Dq\nq0kz04EHJoFi0KDiaW5qKjdPmZnVYeuttyZZ0q6m2Woh8F9cdtmJ/PjHyeimAw+E3/8e9tij8Dqu\n20PJBw1P7jMz+PiQ22uvvZottzwYaRwRbwPvA68hvc/EiWL48OIYAptPntzn5ikzI2mW6tu3P2vX\nPgrsCqwEYODAzSkr+zf33fdjNttsHvAqkydfW7LrV2XVUPOUg4aZlayqKpg5E6ZNq+LuuzcQsc1H\n53r2HMYDD/yc8vLygpopXgjcp2FmHcL69cleEvfdB3/7W7I0x9ChsPfe3bnvvoNZt+5aakZJbdjw\nJGVlZUAy49zBIhsHDTMrWhFJYLj33iRQPPQQ7LorHHIIXHZZsn5T0i+xBa+++mWuvHIvkh0WljF6\n9GkOFM1Q8s1TnqdhVlpWrUrmSdx7b/Jatw6GD09ew4bBVlt9/J6Ny4jcCmwOvE+PHseWzDLq+eJ5\nGu7TMCt6Eck+13/9axIk5syBIUPg0EOT1+c/3/ikusaWEbFNuSPczIrKqlXwf/+XBIp77oHOneGw\nw2CffVayww4vMXBg3ybVENpjwcJi1lDQKNH5jGZWTCLg+efh0kvhoINgu+3gt79N1m+6/35YuBD2\n338GY8bsxNFHj6Ffv/5MmzYj8/Nrllbv0eMAevXakx49DmhwGXWrn2saZtYu1q1LVoP9y1+S15o1\ncMQRyeuggzbdgChfNQUPrc2mQw+59Yxws8KxYkXS5HTXXUnz06BBSZC45RbYfff6+ybq2u61OVu7\nemhtwzwj3DUNs3YVAXPnJkHirrtgwQI4+GA48kg4/HDI+vvbfRJtq0PXNMysbVVXJ/Ml7rwzeQEc\ndRT87GfJjnW113PK0mRU33avDhhtzzUNM2uxd99NhsPecUfy7847J4Hi6KMbHhLb1H273SfRNjzk\n1szybvnypCZxxx3wyCOw335JoPjqV2H77Ru/301OhcvNU2aWFwsXwu23w223wXPPJZPr/vM/Yfp0\n6NUruaaqqoo5cxqvDeSrc9valudpmFm9amZj//Sn8IUvwN57J/MpfvKTpKYxbRocf/zGgDFt2gz6\n9evPwQef0ehcik337YbcfbutcLl5ysw2EQFPPw0335wMhV2zBkaMSF777lv/7nXNaW6qa9/ujr6X\nRSHo0M1Tnqdh1rgIeOKJJEjcckty/PWvw003QXl542s7QfOam0aOPJ5hww5053aB8DwN1zTM6hUB\n//oX/OlPyatLFzjuuCRY7LHHpoEiy6gld2yXDq89ZWZAEiieegp+9KNkWOzxxyeLAd5xR9JXcfHF\nSd9FbsDI2k/h9Z06Btc0zDqABQuSEU7TpydrPh1/fPKqXaOorTm1B8+lKH4duk/DrKNasmRjoFix\nIgkSU6Yke1FIyS/3f/6z4V/uzemn8PpOpc3NU2Yl5I034OqrYf/94YtfTOZVXHYZvPxy8u+Xv5wE\njKxNTh4Wa7W5ecqsyK1enczM/sMf4OGHk4UATzgh2f609jpP0PQmJw+L7XjcPGVWYjZsSPaiuOmm\npBN7yBAYNQquuOIN3nhjEWVlZXTtmp8mJw+LtVxunjIrIvPnJyOfysrgnHNgt92S5Tz+9jfo3HkG\ngwbt2ipNTr1796a8vNwBw4q3eUrSDsDlwJvAixHxyzqucfOUFb0330w6s6dMgWXL4MQT4VvfSgJG\nDTc5WT6VavPUbsDNETFV0rT2zoxZPn3wQVJ7mDw52SP78MPhootg2LBkXkVtbnKytlIwNQ1Jk4Aj\ngRURMTgn/VDgf0ma0ibV1CgkbQncAnwI3BQRU+p4pmsaVlSefz4JFDfdBDvsAKecAgcd9AZvv73I\ns7GtzRTLjPDJwPDcBEmdgCvT9EHASEn909OnABdExDCSYGNWlN5/H268MRkmO3QofPhhUrt4/HHo\n1WsGgwc33k/h2djWVgqmpgEgqR9wV01NQ9JewPiIOCw9/iEQEfFLSYOACcAbwKqIGFfH81zTsIJU\ns0Dg9dcn6z7tuy+MHg1HHJGsAQWejW3tp5j7NLYHXs45XgYMAYiIecBxjT1gwoQJH733arfW3lau\nhKlT4brr4O234dRTYe7cune682xsaytZVretUeg1jWOB4RExJj0eBQyJiLMzPs81DWtXVVVVLFq0\nmHfe2ZkZMz7FbbclndljxsBBB0GnBhqI3U9h7aWYaxqvAH1zjvukaZl5Pw1rL5Mn38rppz/Mhg2n\nE7GS449fyoIFu7P11tnur+mnGD36gE2GxjpgWGspuv00JJWR1DR2S483A54HDgJeA2YDIyNifsbn\nuaZhbW7ePLj00jXccMNakrEm/0FLagnup7C2VhSjpyRNBR4FdpG0VNIpEbEBGAvcB8wDpmcNGDUm\nTJiQua3OrLmqq5PtUQ84IGl+kt5iiy1OIAkYkNsf0VSejW1tpbKycpN+4LoUVE0j31zTsNa2fDlc\ney1cffUGtt32fc48E771rV6sXOn+CCteRVHTMCsmjz+eLOcxYADMmvVv3n57KAsXHsDZZ+/IrbfO\n8LwJK1klX9MYP368O8ItL9avT5qgLr20muXLN3D66R8wcuRadt9913prFO6PsGJS0xF+4YUX1lvT\nKPmgUcrls7ZRVQXXXJNsbvSpT63gpZe+T/fuL1BdvZDzzjuHSy65lZUrn/jo+l699mTmzN9RXl7e\njrk2az43T5k1w7x5cNppsMsuydapU6e+xaJFA1m//lzefXcOa9Y8yMUX/9o721mHUvJBw6OnrCki\nktVlhw9PRkH17ZssInj99bD55v+ma9cykmYogMF07boj5513jvsurCTkbfSUpM4kS3bsnSZtDmwA\nVpP8iTU1Ita2JLOtwc1TltW6dcnyHpdemuyh/f3vw8iR0K3bxmsamqENuO/CSkZDzVONBg1J5cD+\nwP0RMbeO858FjgCejohZechv3jhoWGPefjvpr5g4cQNlZe/x/e8Hxx33SVTnj4s3L7KOoaVBY7e6\ngkUd1+0ELIuI9c3LZv559JTVZ8kSuOwy+P3v4fOfX8Ts2SfRvfsa1q9f3Ggg8IgoK1WtOnoq3W71\nqxFxVUsy2Zpc07Dann4afv1ruOce+Pa3YdSoN9l77108Cc8sR95GT0k6Q9Jjkp4CpgG75yODZq0p\nAior4bDDktduu8HChUnwWL9+4cc6t5u73IdZR9DUVW5XR8TekkZExG3pVqxmBenDD+HOO+EXv4C3\n3oJx4+COOzbt3C4rK8sZMpvUNDxk1qx+TR1yu6ekbsC7kv6bZPXZguYhtx1PdTVMmQIDBnzAj370\nPmPGrGT+/GTDo9yAAd4m1SxX3hcslLRLRLyQvh8H/Dsibm1JJluT+zQ6ltWrYdIkuOQS2GKLFbz4\n4n/RvfsSqqvduW3WFC0aPZU+YCtgR5JhtdV5zl+rcdDoGFauhKuugokTYa+94Iwz3mbEiJ3duW3W\nTC3qCJc0AngRuAN4RtI2ec6fWbNUVcH558NOO33Iww+/wc03v8Udd8BWW73kzm2zVpKlT2MosG1E\n9AEOATyTydrVK6/A974Hu+4Kjz/+Eu+//yUefXQ4w4d/jmnTZtTq3AZ3bpvlT5agMS8i1gFExMvA\nspoTkvq0VsbMalu4EE4/PRkyK0Fl5Zs8+uiXWbfuRlaufII1ax5k9OgzAdy5bdZKsgy5HSjpKznH\nAyQdkr4/Hhid/2zlz4QJEzwjvMjNnw8//zncffeHHHPMcv7+964MGPBp5sxJ5lisWfPxZqiRI49n\n2LAD3blt1gQ1M8IbkmUZkWeAfwJ1dYp8MSIG15FeENwRXtyeegouvhgeeggqKp7hzjuPplu3LT9a\n6mPYsAO9papZK2jp2lPDImJmU88VAgeN4vTYY0mwePJJOOccGDGiioED6w4OM2c+4AUEzfKsoaDR\nYPNUOpHvyfrO5wYMSTukfR5mmdXMj+jXr4y5c3tz8cVJ38UPfwi33ALdu8OcOYvdDGVWIBoMGhGx\nTtLBkrYA7oiINbWvkfRJ4BvAc4CDhmU2bdoMvv3tM+nU6STWrBnFNtt04+c/78UJJ0CXLhuva2yp\nj969eztYmLWRLM1TWwDfAwLoC3wAdGHjJkzLgOsjYmXrZrXp3DxVuJYvr6Jv3x9RXX0F0ANYQvfu\nX2Lp0ufqDADex8Ks7bS0T+MaYCXQJ30dFhGr857LVuCgUXjWrYObboKLLlrLq68uYMOGPT4616vX\nnsyc+TvKy8vrvNdLfZi1jZYujT43Is6NiBNJhtj6zztrsvfeg4sueo++fdczdep6rrxyLV26HERT\nJuD17t2b8vJyBwyzdpQlaHy093dELAdWtV528s+r3LavN9+E8eNhu+3WcuGF/8fq1afw+OPb8v77\nf+OGGzwBz6yQ5GWVW0kvAfcC/0pfn61Z2VbSZyLi9bzkthW4ear95G6nevjha7jllqGsW3c9tYfM\nAm5yMiswLW2euhG4G9gBuBi4It297xLgkrzl0krC3Llw0kmw557QtSs8+yx85zvP0r37BupaQNBN\nTmbFpdFlRCLiZ+nbe2vSJO0EfBkY00r5siJQe47FJZcks7i/8x244gr45CeT67p08e54ZqWiqdu9\nAhARC4GFkl7Jc36sSCRzLM5GGs3atd9g2227cdFFvT62nSps3B1v9OgDNhky69qFWfFp0s59xcZ9\nGq3j3/9+g/79L+eDD84HugKL6N69nKVL5zcYCDxk1qw4NHsZEbNcixbB5ZfDDTd8EmlvkoABsCNd\nu/b9qI+iPp65bVb8snSEWwdSVVXFnDlzqKqqAiAiWWV2xAgoL0+W95g1ayWdO4/CmxyZdTxFW9OQ\ntB9wIkkZBkTEfu2cpaJXs1RH165lrFv3GieffDv/+MeXef99+O53k5ncm28OsJX7KMw6qKLv05B0\nNPCZiLiujnPu08ioqqoq3ZviYWAgUE2nTrP4wx/KOf74/6BTHXVS91GYlaaWztNoE5ImSVqRbvqU\nm36opAWSXpB0bh23ngBMbZtclqYIuO22N/ngg+kkAQOgCz17jmPnnV+oM2CAl/Uw64gKJmgAk4Hh\nuQmSOgFXpumDgJGS+uec3wF4JyLeb8uMlor33oNrroHBg+E3v9kZmAk8m551P4WZfVzBBI2IeAR4\nu1byEODFiFgSEdXAdODonPOjSYKNNaB25/azz8JZZ0HfvnDffTBxIjz/fGemTNmTHj2Gei0oM6tX\noXeEb8+mGzstIwkkAETEhMYekLv4VkVFBRUVFXnLXDHYuA/F51izZjA77XQxq1b15tRT4ZlnoE+f\njdd6FzyzjqmysjLzwq4F1REuqR9wV0QMTo+PBYZHxJj0eBQwJCLOzvi8Dt0RXlVVxQ47HMm6dXcD\nvYFVdO06lkWLfs122zkgmFndiqIjvB6vkOwWWKNPmpZZR1wafe1amDYNDjusG+vX/5UkYABsQffu\nz/DKK4vbMXdmVqjysjR6W5JURlLT2C093gx4HjgIeA2YDYyMiPkZn9ehahrz5sF118Ef/wh77AHH\nH/8uY8fuwtq191F7SXI3PZlZfYpiGRFJU4EKYCtJS4HxETFZ0ljgPpJa0aSsAaOjWLUKZsyASZNg\n6VI45RSYPRt23BGgF5tvPtGT8MwsbwqqppFvkmL8+PEl1wEeAY8+mgSK22+HoUNh9Gg47DDoXMef\nAZ6EZ2YcIeYjAAAMuElEQVRZ1HSIX3jhhfXWNEo+aJRS+V59NdkJ78Ybk+PRo5MNj7bZpl2zZWYl\npiiap1rLhAkTiqamUVeNYN06uPNOmDwZHnsMvv51uOEG2HtvUJ3fUjOz5sky9NY1jQKx6WKBiznv\nvOksX34wM2YkM7ZPOSVZaTZZMNDMrPU0VNNw0CgAH18scC3Sq4wbtzVnnLE5XsnDzNpSMc/TaLFC\nn6fx3ntw1VWrqK7+CxsXC+xOz55f59hjn3PAMLM2U3TzNPKtEGoadfVTbNgADz6Y7E/x5z9Defk6\nZs06k+rq7wK74fkUZtaeOnRNoz1NmzaDfv36c/DBZ9CvX39++ct7GTcO+vWDceOSCXgLFsD993dj\nypRD6NGjwosFmllBc02jlWzsp3gEGACsQXqTsWO3ZMyYTzBoUN33eD6FmbU3D7lt4yG377wDEye+\nx/r195IEDIAe9Ox5NKNGXcOgQeV13te7d28HCzNrNx5y20o1jbpqBGvXwj33JOs+3X8/7LffOmbO\nHMP69T/A/RRmVkw6dE0j3zadT7GUc875E8uXH8Dtt8Puu8MJJySLBn7qU92YNu1wRo+u8LpPZlYy\nXNNogqqqKvr27c/atY8DnwOqkeZzwQU7ceqpPTfZ0Cj3HvdTmFkx6dA1jXz1aSxYAL/5zXrWr3+S\njVt8dGGLLU7miCN+R58+7qcws+LmPo2MNY36agNLl8L06cmGRitWwFFHrebGG49i3bpL8f4UZlaq\nOnRNozG5fRTr1y/m0ksnsWHDMUybltQuRoyASy+Fr3wFNtvsEwwdepr3pzCzDqtD1zQ2zqV4CBgE\nrAI+5Nhju3Hyyd055BDo2rXu+9xPYWalyjWNOqxaBddc8y7V1beRBAyALdhii30499yJlJfX3UcB\n7qcws46r5JcRqaqq+uj9mjVwyy1w3HHQpw889FAfpBnAs+kVz/DBB89T5lUCzczqVPI1jW233ZFj\nj/0fOnc+m7/8Bb70JfjmN+Gaa2CrrboxbdpQRo8e6j4KM+vwPHpKCqimU6fZXHzxbpx88hZ1bo3q\nPgozs4069CZMEPTqtSczZ/6uwX4KMzNLdPCl0Z+hunqJ+ynMzPKg5IOG96YwM8ufkm+eev311x0w\nzMyaoEP3aZRy+czMWkOH7tOYMGFCo0PIzMwsGXI7YcKEBq9xTcPMzDbRoWsaZmaWPw4aZmaWmYOG\nmZll5qBhZmaZOWiYmVlmDhpmZpaZg4aZmWVWtPtpSBLwU6AXMCcibmrnLJmZlbxirmkcDfQB1gPL\n2jkvZmYdQsEEDUmTJK2Q9Eyt9EMlLZD0gqRzc07tCvw9In4AnNmmmTUz66AKJmgAk4HhuQmSOgFX\npumDgJGS+qenlwFvp+83tFUmzcw6soIJGhHxCBuDQI0hwIsRsSQiqoHpJM1SALcBh0qaCMxqu5ya\nmXVchd4Rvj3wcs7xMpJAQkSsAU5t7AG5KzZWVFRQUVGR1wyamRW7ysrKzKuBF9Qqt5L6AXdFxOD0\n+FhgeESMSY9HAUMi4uyMz/Mqt2ZmTVTMq9y+AvTNOe6TpmXm/TTMzLIpuv00JJWR1DR2S483A54H\nDgJeA2YDIyNifsbnuaZhZtZERVHTkDQVeBTYRdJSSadExAZgLHAfMA+YnjVg1HBNw8wsm6KraeSb\naxpmZk1XFDUNMzMrfCUfNNw8ZWaWjZun3DxlZtZkbp4yM7O8KPmg4eYpM7Ns3Dzl5ikzsyZz85SZ\nmeVFyQcNN0+ZmWXj5ik3T5mZNZmbp8zMLC8cNMzMLDMHDTMzy6zkg4Y7ws3MsnFHuDvCzcyazB3h\nZmaWFw4aZmaWmYOGmZllVvJBwx3hZmbZuCPcHeFmZk3mjnAzM8sLBw0zM8vMQcPMzDJz0DAzs8wc\nNMzMLDMHDTMzy6zkg4bnaZiZZeN5Gp6nYWbWZJ6nYWZmeeGgYWZmmTlomJlZZg4aZmaWmYOGmZll\n5qBhZmaZFW3QkDRU0kOSrpb0lfbOj5lZR1C0QQMIYBXQDVjWznkxM+sQCiZoSJokaYWkZ2qlHypp\ngaQXJJ1bkx4RD0XEEcAPgYvaOr8tVeqz1F2+4lbK5SvlskHrl69gggYwGRiemyCpE3Blmj4IGCmp\nf6373gG6tkkO88j/cYuby1e8Srls0Prl69yqT2+CiHhEUr9ayUOAFyNiCYCk6cDRwAJJXyMJJv9B\nEljMzKyVFUzQqMf2wMs5x8tIAgkRcTtwe3tkysysoyqoBQvTmsZdETE4PT4WGB4RY9LjUcCQiDg7\n4/MKp3BmZkWkvgULC72m8QrQN+e4T5qWSX2FNjOz5imkjnAApa8ac4CdJfWT1BX4JnBnu+TMzMwK\nJ2hImgo8CuwiaamkUyJiAzAWuA+YB0yPiPntmU8zs46sYIJGRJwQEdtFRLeI6BsRk9P0eyJi14j4\nXET8or3zmYWkPpIekDRP0lxJZ6fpu0t6TNKTkmZL+lLOPT+S9KKk+ZIOab/cN66R8j0q6WlJf5bU\nM+eeYipfN0n/SL9PcyWNT9M/Jek+Sc9L+puk/8i5pxTK93VJz0raIGnPWveUQvl+leb/KUm3SuqV\nc09RlK+Bsl2U/tw9KeleSdvk3JPfskWEX3l+AdsAe6TvewILgAHA34BD0vTDgAfT9wOBJ0n6mMqA\nl0gHKRTiq4HyzQb2S9NPBi4qxvKlef5E+u9mwOMko/Z+CYxL088FflFi5dsV+BzwALBnzrUDSqR8\nw4BOafovgJ8X4/evnrL1zDk/Fri6tcpWMDWNUhIRyyPiqfT9eyS/VLcDPiSZVwLwSTZ26h9F0vT2\nQUQsBl4kHVpciOop3/bA5yLikfSymcCx6fuiKh9ARKxO33Yj+YELkjlCU9L0KcAx6fuSKF9EPB8R\nL7JpvyIk5S6F8s2MiA/T9MdJBtZAkX3/6inbezmXbE7yuwZaoWwOGq1MUhmwB/AP4HvAJZKWAr8C\nfpReVns+yitpWsHLKd/jwDxJR6WnvsHGH8qiK5+kTpKeBJYD90fEHGDriFgBSeAEPpNeXirlq08p\nlu/bwF/T90VVvvrKJuln6e+WE4AL0svzXjYHjVaUtunfAnwn/Uvgv9L3fUkCyA3tmb+WqqN8o4Gz\nJM0h+WtnfXvmryUi4sOI+AJJ4BsiaRBJbWOTy9o+Z/lRq3xfljSwvfOUTw2VT9KPgeqImNZuGWyB\n+soWEeenv1v+SNJE1SocNFqJpM4kv1Bviog/p8n/GRF3AETELUB5mv4KsEPO7U2aj9Ie6ipf2rwx\nPCLKgenAv9PLi658NSLiXaASOBRYIWlrgLSj8fX0smIv34Mk5atPyZRP0snA4SR/jdcoyvI18L2b\nCoxI3+e9bA4arecG4LmImJiT9oqkoQCSDiJpX4Rk7sk3JXWVtCOwM0mnciH7WPkk9U7/7QScD1yT\nniqq8kn6dM3IKEk9gIOB+STlODm97D+Bmj8GSqF8C2pflvO+JMon6VDg/wFHRcS6nFuKpnwNlG3n\nnMuOYeP3M+9lK/QZ4UVJ0r7AicDctO0xgPOA04DLJW0GrAXGAETEc5L+BDwHVANnRjr0oRA1UL5d\nJJ2VHt8WETdC8ZUP2BaYkga/TsCMiPirpMeBP0n6NrCEpN+mlMp3DHAF8GngbklPRcRhJVS+F0lW\nxL5fEsDjEXFmkZWvvrLdImkXkg7wJcAZ0Dr/Nwtq7SkzMytsbp4yM7PMHDTMzCwzBw0zM8vMQcPM\nzDJz0DAzs8wcNMzMLDMHDTMzy8xBw8zMMnPQMGtDkrq1dx7MWsJBwyyPJH1O0j2Sxki6X9L1kk6X\n9ISkI4CekoZJekbStZK2S++rlHSNpK0k7SapvJGPMmsXXnvKLL++QLIgXrWkEcCvIuIFSWuAXhHx\nJjBT0j+BP0TEq5L6Ar+JiLvSZ7wp6fuSnoyID9qpHGZ1ck3DLL9eiIjq9P0uEfFC+r4MuCPnOkGy\n7ziwT07AqHE/cFxrZtSsORw0zPKoZhvcdKnql3JObRkRa2pd3oNk35FFdTxnLrB3a+XTrLkcNMxa\nxxA23begex3XHA78FPjfNsmRWR44aJi1jiEk+8LXqKv/8NaIeARYKOmkOs5v3io5M2sBBw2z1lHO\npjWNDQ1c+0PgAkmfqJXe0D1m7cJBwyyPJO0u6QfAYGBEzRa4wOqca4YBXwRGpUNug6Qm8sd0S05q\n32NWKLxzn1kbkHQOMCki3sl4/WeBoRFxQ+vmzKxpXNMwaxvXk+4pntERwNRWyotZszlomLWBiFgJ\nPCdph8aulbQT8HRErG39nJk1jZunzAqMpK4Rsb6982FWFwcNMzPLzM1TZmaWmYOGmZll5qBhZmaZ\nOWiYmVlmDhpmZpaZg4aZmWX2/wE+kz/VZK7vAgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x99190b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#define a script to calculate the phase diagram\n",
    "def PhaseDiagram():\n",
    "    from scipy.optimize import bisect as SOLVE\n",
    "\n",
    "    #specified gas composition\n",
    "    x=np.array([1.])\n",
    "\n",
    "    #definition of residue function for numerical solution along the PxT phase diagram\n",
    "    def RES(T,P):\n",
    "        Dmu_EL_Lw=saito_EL_Lw(T,P)\n",
    "        fug_G=PR_fug_G(T,P,x)[0]\n",
    "        Dmu_H_EL=vdwnp_H_EL(T,fug_G)\n",
    "        return Dmu_EL_Lw+Dmu_H_EL\n",
    "    \n",
    "    #define a grid\n",
    "    Tgrid = np.zeros(100)\n",
    "    Pgrid = np.logspace(5,9,100)\n",
    "    \n",
    "    #seek solution for T for each P in the grid\n",
    "    for i in range(100):\n",
    "\n",
    "        #call the solver imported from scipy.optimize\n",
    "        ans=SOLVE(\n",
    "            #define a lambda function for thew numerical method\n",
    "            #so that the numerical method can treat the RES of T and P as a function of 1 variable and tehrefore provide solutions for T for each given P\n",
    "            lambda T, P=Pgrid[i]: RES(T,P)\n",
    "            #set lower and upper search limits and tolerances\n",
    "            , 150, 400, xtol=1e-9, rtol=1e-9, maxiter=100, full_output=True, disp=True)\n",
    "\n",
    "    #record solution after each point in the grid\n",
    "        Tgrid[i]=ans[0]\n",
    "\n",
    "    #plot:\n",
    "    \n",
    "    #import plotting tools\n",
    "    import matplotlib.pyplot as plt\n",
    "    %matplotlib inline\n",
    "    \n",
    "    #activate figure\n",
    "    #plot calculated data\n",
    "    plt.figure(1)\n",
    "    plt.semilogy(Tgrid,Pgrid)\n",
    "    \n",
    "    #set plot limits\n",
    "    plt.xlim([273,333])\n",
    "    plt.ylim([1e5,1e10])\n",
    "\n",
    "    #Let's also plot some experimental data along our calculations\n",
    "    \n",
    "    #retrieve data from a text file provided with this notebook\n",
    "    table = np.loadtxt('support-files/C1-sI-hydrates-HLwGequilibria-expData.txt', #data from NIST\n",
    "                       dtype='float', \n",
    "                       comments='#', \n",
    "                       converters=None,\n",
    "                       skiprows=0,\n",
    "                       usecols=None,\n",
    "                       unpack=False,\n",
    "                       ndmin=0)\n",
    "    \n",
    "    #split data into arrays for T and P\n",
    "    T_eq = table[:,0]\n",
    "    P_eq = table[:,1]    \n",
    "    \n",
    "    #plot exp data as scatter\n",
    "    plt.scatter(T_eq,P_eq)\n",
    "    \n",
    "    #format\n",
    "    plt.title(r'$[\\mathrm{CH}_4; \\mathrm{H}_2\\mathrm{O}]$'+'\\n'+'$\\mathrm{HsI-Lw-V}$'+' Univariant equilibria')\n",
    "    plt.xlabel(r'$T\\mathrm{(K)}$')\n",
    "    plt.ylabel(r'$P\\mathrm{(Pa)}$')\n",
    "\n",
    "    #show plot!\n",
    "    plt.show()    \n",
    "    return\n",
    "\n",
    "#call the phase diagram script\n",
    "_=PhaseDiagram()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Coming soon:\n",
    ">- preserve this notebook for beginner lectures\n",
    ">- present improvements in intermediate notebooks\n",
    "\n",
    "* add ice below 273 K\n",
    "* add Langmuir parameters for\n",
    "* add acurate correlation for Lw and H volumes and its Poynting integrals\n",
    "* add activity coefficient model for water soluble inhibitors\n",
    "* add ethanol' inhibitor-guest dual-role\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "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": 1
}