acids-bases-first.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Acids and Bases I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy.optimize as spo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Definitions \n",
    "\n",
    "The reaction of acids with water has to be recorded as follows: \n",
    "\n",
    "$$ AH(aq) + H_2O(l) \\rightleftharpoons A^-(aq) + H_3O^+(l) $$\n",
    "\n",
    "where $AH/A^-$ is an acid-base couple."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using a numerical method for Acid/Bases calculations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#https://github.com/rnelsonchem/pHcalc/blob/master/pHcalc/pHcalc.py\n",
    "#https://pypi.org/project/pHcalc/\n",
    "\n",
    "class Neutral(object):\n",
    "    \"\"\"A nonreactive ion class.\n",
    "    This object defines things like K+ and Cl-, which contribute to the\n",
    "    overall charge balance, but do not have any inherent reactivity with\n",
    "    water.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    charge : int\n",
    "        The formal charge of the ion.\n",
    "    conc : float\n",
    "        The concentration of this species in solution.\n",
    "    Attributes\n",
    "    ----------\n",
    "    charge : int\n",
    "        The formal charge of the ion.\n",
    "    conc : float\n",
    "        The concentration of this species in solution.\n",
    "    \"\"\"\n",
    "    def __init__(self, charge=None, conc=None):\n",
    "        if charge == None:\n",
    "            raise ValueError(\n",
    "                \"The charge for this ion must be defined.\")\n",
    "\n",
    "        self.charge = charge \n",
    "        self.conc = conc\n",
    "\n",
    "    def alpha(self, pH):\n",
    "        '''Return the fraction of each species at a given pH.\n",
    "        Parameters\n",
    "        ----------\n",
    "        pH : int, float, or Numpy Array\n",
    "            These are the pH value(s) over which the fraction should be\n",
    "            returned.\n",
    "        Returns\n",
    "        -------\n",
    "        Numpy NDArray\n",
    "            Because this is a non-reactive ion class, this function will\n",
    "            always return a Numpy array containing just 1.0's for all pH\n",
    "            values.\n",
    "        '''\n",
    "        if isinstance(pH, (int, float)):\n",
    "            length = 1\n",
    "        else:\n",
    "            length = len(pH)\n",
    "        ones = np.ones(length).reshape(-1,1)\n",
    "        return ones\n",
    "\n",
    "\n",
    "        \n",
    "class Acid(object):\n",
    "    '''An acidic species class.\n",
    "    This object is used to calculate a number of parameters related to a weak\n",
    "    acid in an aqueous solution. \n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    Ka : None (default), float, list, Numpy Array\n",
    "        This defines the Ka values for all acidic protons in this species. It\n",
    "        can be a single Ka value (float), a list of floats, or a Numpy array\n",
    "        of floats. Either this value or pKa needs to be defined. The other\n",
    "        will then be calculated from the given values.\n",
    "    pKa : None (default), float, list, Numpy Array\n",
    "        The pKa value(s) for all the acidic protons in this species.  This\n",
    "        follows the same rules as Ka (See Ka description for more details),\n",
    "        and either this value or Ka must be defined.\n",
    "    charge : None (default), int\n",
    "        This is the charge of the fully protonated form of this acid. This\n",
    "        must be defined.\n",
    "    conc : None (default), float\n",
    "        The formal concentration of this acid in solution. This value must be\n",
    "        defined.\n",
    "    Note\n",
    "    ----\n",
    "    There is no corresponding Base object. To define a base, you must use a\n",
    "    combination of an Acid and Neutral object. See the documentation for\n",
    "    examples.\n",
    "    '''\n",
    "    def __init__(self, Ka=None, pKa=None, charge=None, conc=None):\n",
    "        # Do a couple quick checks to make sure that everything has been\n",
    "        # defined.\n",
    "        if Ka == None and pKa == None:\n",
    "            raise ValueError(\n",
    "                \"You must define either Ka or pKa values.\")\n",
    "        elif charge == None:\n",
    "            raise ValueError(\n",
    "                \"The maximum charge for this acid must be defined.\")\n",
    "\n",
    "        # Make sure both Ka and pKa are calculated. For lists of values, be\n",
    "        # sure to sort them to ensure that the most acidic species is defined\n",
    "        # first.\n",
    "        elif Ka == None:\n",
    "            if isinstance(pKa, (int, float)):\n",
    "                self.pKa = np.array( [pKa,], dtype=float)\n",
    "            else:\n",
    "                self.pKa = np.array(pKa, dtype=float)\n",
    "                self.pKa.sort()\n",
    "            self.Ka = 10**(-self.pKa) \n",
    "        elif pKa == None:\n",
    "            if isinstance(Ka, (int, float)):\n",
    "                self.Ka = np.array( [Ka,], dtype=float)\n",
    "            else:\n",
    "                self.Ka = np.array(Ka, dtype=float)\n",
    "                # Ka values must be in reverse sort order\n",
    "                self.Ka.sort()\n",
    "                self.Ka = self.Ka[::-1]\n",
    "            self.pKa = -np.log10(self.Ka)\n",
    "        # This temporary Ka array will be used to calculate alpha values. It\n",
    "        # starts with an underscore so that it won't be confusing for others.\n",
    "        self._Ka_temp = np.append(1., self.Ka)\n",
    "        \n",
    "        # Make a list of charges for each species defined by the Ka values.\n",
    "        self.charge = np.arange(charge, charge - len(self.Ka) - 1, -1)\n",
    "        # Make sure the concentrations are accessible to the object instance.\n",
    "        self.conc = conc \n",
    "\n",
    "    def alpha(self, pH):\n",
    "        '''Return the fraction of each species at a given pH.\n",
    "        Parameters\n",
    "        ----------\n",
    "        pH : int, float, or Numpy Array\n",
    "            These are the pH value(s) over which the fraction should be\n",
    "            returned.\n",
    "        Returns\n",
    "        -------\n",
    "        Numpy NDArray\n",
    "            These are the fractional concentrations at any given pH. They are\n",
    "            sorted from most acidic species to least acidic species. If a\n",
    "            NDArray of pH values is provided, then a 2D array will be\n",
    "            returned. In this case, each row represents the speciation for\n",
    "            each given pH.\n",
    "        '''\n",
    "        # If the given pH is not a list/array, be sure to convert it to one\n",
    "        # for future calcs.\n",
    "        if isinstance(pH, (int, float)):\n",
    "            pH = [pH,]\n",
    "        pH = np.array(pH, dtype=float)\n",
    "\n",
    "        # Calculate the concentration of H3O+. If multiple pH values are\n",
    "        # given, then it is best to construct a two dimensional array of\n",
    "        # concentrations.\n",
    "        h3o = 10.**(-pH)\n",
    "        if len(h3o) > 1:\n",
    "            h3o = np.repeat( h3o.reshape(-1, 1), len(self._Ka_temp), axis=1)\n",
    "\n",
    "        # These are the powers that the H3O+ concentrations will be raised.\n",
    "        power = np.arange(len(self._Ka_temp))\n",
    "        # Calculate the H3O+ concentrations raised to the powers calculated\n",
    "        # above (in reverse order).\n",
    "        h3o_pow = h3o**( power[::-1] )\n",
    "        # Calculate a cumulative product of the Ka values. The first value\n",
    "        # must be 1.0, which is why _Ka_temp is used instead of Ka.\n",
    "        Ka_prod = np.cumproduct(self._Ka_temp)\n",
    "        # Multiply the H3O**power values times the cumulative Ka product.\n",
    "        h3o_Ka = h3o_pow*Ka_prod\n",
    "\n",
    "        # Return the alpha values. The return signature will differ is the\n",
    "        # shape of the H3O array was 2-dimensional. \n",
    "        if len(h3o.shape) > 1:\n",
    "            den = h3o_Ka.sum(axis=1)\n",
    "            return h3o_Ka/den.reshape(-1,1)\n",
    "        else:\n",
    "            den = h3o_Ka.sum()\n",
    "            return h3o_Ka/den\n",
    "\n",
    "\n",
    "class System(object):\n",
    "    '''An object used to define an a system of acid and neutral species.\n",
    "    This object accepts an arbitrary number of acid and neutral species\n",
    "    objects and uses these to calculate the pH of the system. Be sure to\n",
    "    include all of the species that completely define the contents of a\n",
    "    particular solution.\n",
    "    Parameters\n",
    "    ----------\n",
    "    *species \n",
    "        These are any number of Acid and Neutral objects that you'd like to\n",
    "        use to define your system.\n",
    "    Attibutes\n",
    "    ---------\n",
    "    species : list\n",
    "        This is a list containing all of the species that you input.\n",
    "    pHsolution \n",
    "        This is the full minimization output, which is defined by the function\n",
    "        scipy.optimize.minimize. This is only available after running the\n",
    "        pHsolve method.\n",
    "    pH : float\n",
    "        The pH of this particular system. This is only calculated after\n",
    "        running the pHsolve method.\n",
    "    '''\n",
    "    def __init__(self, *species):\n",
    "        self.species = species\n",
    "\n",
    "\n",
    "    def _diff_pos_neg(self, pH):\n",
    "        '''Calculate the charge balance difference.\n",
    "        Parameters\n",
    "        ----------\n",
    "        pH : int, float, or Numpy Array\n",
    "            The pH value(s) used to calculate the different distributions of\n",
    "            positive and negative species.\n",
    "        Returns\n",
    "        -------\n",
    "        float or Numpy Array\n",
    "            The absolute value of the difference in concentration between the\n",
    "            positive and negatively charged species in the system. A float is\n",
    "            returned if an int or float is input as the pH: a Numpy array is\n",
    "            returned if an array of pH values is used as the input.\n",
    "        '''\n",
    "        twoD = True\n",
    "        if isinstance(pH, (int, float)) or pH.shape[0] == 1:\n",
    "            twoD = False\n",
    "        else:\n",
    "            pH = np.array(pH, dtype=float)\n",
    "        # Calculate the h3o and oh concentrations and sum them up.\n",
    "        h3o = 10.**(-pH)\n",
    "        oh = (10.**(-14))/h3o\n",
    "        x = (h3o - oh)\n",
    "\n",
    "        # Go through all the species that were given, and sum up their\n",
    "        # charge*concentration values into our total sum.\n",
    "        for s in self.species:\n",
    "            if twoD == False:\n",
    "                x += (s.conc*s.charge*s.alpha(pH)).sum()\n",
    "            else:\n",
    "                x += (s.conc*s.charge*s.alpha(pH)).sum(axis=1)\n",
    "        \n",
    "        # Return the absolute value so it never goes below zero.\n",
    "        return np.abs(x)\n",
    "        \n",
    "\n",
    "    def pHsolve(self, guess=7.0, guess_est=False, est_num=1500, \n",
    "                method='Nelder-Mead', tol=1e-5):\n",
    "        '''Solve the pH of the system.\n",
    "        The pH solving is done using a simple minimization algorithm which\n",
    "        minimizes the difference in the total positive and negative ion\n",
    "        concentrations in the system. The minimization algorithm can be\n",
    "        adjusted using the `method` keyword argument. The available methods\n",
    "        can be found in the documentation for the scipy.optimize.minimize\n",
    "        function.\n",
    "        \n",
    "        A good initial guess may help the minimization. It can be set manually\n",
    "        using the `guess` keyword, which defaults to 7.0. There is an\n",
    "        automated method that can be run as well if you set the `guess_est`\n",
    "        argument. This will override whatever you pass is for `guess`. The\n",
    "        `est_num` keyword sets the number of data points that you'd like to\n",
    "        use for finding the guess estimate. Too few points might start you\n",
    "        pretty far from the actual minimum; too many points is probably\n",
    "        overkill and won't help much. This may or may not speed things up.\n",
    "        Parameters\n",
    "        ----------\n",
    "        guess : float (default 7.0)\n",
    "            This is used as the initial guess of the pH for the system. \n",
    "        guess_est : bool (default False)\n",
    "            Run a simple algorithm to determine a best guess for the initial\n",
    "            pH of the solution. This may or may not slow down the calculation\n",
    "            of the pH.\n",
    "        est_num : int (default 1500)\n",
    "            The number of data points to use in the pH guess estimation.\n",
    "            Ignored unless `guess_est=True`.\n",
    "        method : str (default 'Nelder-Mead')\n",
    "            The minimization method used to find the pH. The possible values\n",
    "            for this variable are defined in the documentation for the\n",
    "            scipy.optimize.minimize function.\n",
    "        tol : float (default 1e-5)\n",
    "            The tolerance used to determine convergence of the minimization\n",
    "            function.\n",
    "        '''\n",
    "        if guess_est == True:\n",
    "            phs = np.linspace(0, 14, est_num)\n",
    "            guesses = self._diff_pos_neg(phs)\n",
    "            guess_idx = guesses.argmin()\n",
    "            guess = phs[guess_idx]\n",
    "            \n",
    "        self.pHsolution = spo.minimize(self._diff_pos_neg, guess, \n",
    "                method=method, tol=tol)\n",
    "        \n",
    "        if self.pHsolution.success == False:\n",
    "            print('Warning: Unsuccessful pH optimization!')\n",
    "            print(self.pHsolution.message)\n",
    "            \n",
    "        if len(self.pHsolution.x) == 1:\n",
    "            self.pH = self.pHsolution.x[0]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# pH of strong acid solutions\n",
    "\n",
    "First of all, HCl completely dissociates in water to give equal amounts of $H_3O^+$ and $Cl^-$. Because $H_3O^+$ is adjusted internally, all you need to define is $Cl^-$. This implies a single equivalent of $H_3O^+$ in order to balance the charge of the system.\n",
    "\n",
    "$$ HCl(g) +H_2O(l) \\longrightarrow H_3O^+(l) + Cl^-(aq)$$\n",
    "\n",
    "You just need to define the amount of Cl-. The solver will find the correct $H_3O^+$ concentration. \"M\" stands for molar, which means $1\\,mol.L^{-1}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HCl at 1e-1 M has a  pH = 7.00\n"
     ]
    }
   ],
   "source": [
    "# [HCl] = 1.0 x 1E-1 \n",
    "cl1 = Neutral(charge=-1, conc=1E-30) # [H3O+] in mol.L-1 \n",
    "s = System(cl1)\n",
    "s.pHsolve()\n",
    "print('HCl at 1e-1 M has a  pH = {:0.2f}'.format(s.pH))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HCl at 1e-8 M has a pH = 6.98\n"
     ]
    }
   ],
   "source": [
    "# [HCl] = 1.0 x 1E-8   \n",
    "\n",
    "cl2 = Neutral(charge=-1, conc=1e-8)\n",
    "s = System(cl2)\n",
    "s.pHsolve()\n",
    "print('HCl at 1e-8 M has a pH = {:0.2f}'.format(s.pH))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# pH of strong basic solutions\n",
    "\n",
    "Again for strong bases the dissociation is going to be complette :\n",
    "\n",
    "$$ KOH  \\longrightarrow HO^- + K^+$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "KOH at 0.1 M has a pH = 13.00\n"
     ]
    }
   ],
   "source": [
    "# KOH, just need to define the amount of K+.\n",
    "a = Neutral(charge=+1, conc=0.1) # to define a base another combination of is needed\n",
    "s = System(a)\n",
    "s.pHsolve()\n",
    "print('KOH at 0.1 M has a pH = {:0.2f}'.format(s.pH))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "KOH at 0.1 M has a pH = 7.00\n"
     ]
    }
   ],
   "source": [
    "# KOH, just need to define the amount of K+.\n",
    "a = Neutral(charge=+1, conc=1E-10) # to define a base another combination of is needed\n",
    "s = System(a)\n",
    "s.pHsolve()\n",
    "print('KOH at 0.1 M has a pH = {:0.2f}'.format(s.pH))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Question 1)\n",
    "Write a cell with a solution that could give you a pH of 3.3\n",
    "You need to chose a strong acid.\n",
    "\n",
    "Answer here :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "### BEGIN SOLUTION\n",
    "\n",
    "### END SOLUTION"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Question 2)\n",
    "Write a cell with a solution that could give you a pH of 8.7\n",
    "You need to chose a strong base.\n",
    "\n",
    "Answer here :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "### BEGIN SOLUTION\n",
    "\n",
    "### END SOLUTION"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  Titration of carbonic acid\n",
    "\n",
    "The following theoretical example aims at showing the double reaction : \n",
    "\n",
    "$$ H_2CO_3(aq) + HO^-(aq) \\rightleftharpoons HCO_3^{-}(aq) +H_2O(l) $$ pKa1=6.4 at 25°C\n",
    "\n",
    "$$ HCO_3^{-}(aq) + HO^-(aq) \\rightleftharpoons CO_3^{2-}(aq) +H_2O(l) $$ pKa2=10.3 at 25°C\n",
    "\n",
    "A picture is diplayed and saved on the same directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XecXXWd//HXJ9NnMpM2k94TEiAIBEJHuoAsig2RhRVcWRZ1EVRkwYplV1zXn65lURYV6aKABVBAJCAlhFCTEJJAep9JmZLpM5/fH+c74WaYlsy998yd+34+Hvcx957yPZ/vOXe+n9Pu95i7IyIi2WtI3AGIiEi8lAhERLKcEoGISJZTIhARyXJKBCIiWU6JQEQkyykRDGBm9iUzuyXFy5hsZnVmlpPK5XSx3DFm9pSZ1ZrZ99O0zPlmdlkKyn23mS1Pdrk9LK/H74WZrTGzM7oZd6uZfTt10SWfmV1kZo/GHcdgpkQQo9AAd7zazawh4fNF7v6f7n5ZmHaqmbmZ5fZzmXs1Eu6+zt2Huntbf+uzjy4HqoAyd/9CVxOY2dFm9rCZ7TKzHWa20Mw+kd4we+fuf3f32Wlc3p7vRZzC93Fmp2E3mNkdCZ/NzL5oZivD93udmd1oZgV9XY673+nuZyYzdtmbEkGMQgM81N2HAuuA9yUMu3Nfy+tvkkizKcDr3s0vGs3sOOBvwJPATGAU8Cngvfu6oNAY6bsejx8RJf2PA6VE2+804N44g5K96Z9jAOu0d/VU+LsrHDEcZ2aXmtkzZvYDM9sB3GBmM8zsb2a23cyqzOxOMxseyrsdmAz8KZRxbecjDTMbb2Z/DHvgb5rZv3SK514zuy2c0llqZvN6iP94M3vBzKrD3+PD8FuBS4BrQxxdncb4HvBrd/+uu1d55EV3/2goY4SZPWhmlWa2M7yfmLDs+Wb2H2b2DFAPTA+jZoQji2oz+4OZjUyY5/2hTrvC/AcljFtjZteY2Wth3t+YWWEYd4qZbUiYdpKZ3R9i225mP+lm/RxtZs+F5W02s5+YWX7C+Dlm9ljYFlvN7EsJ2yFxr/ufzGxtWNaXu9seCcpDubVm9qSZTQnl/NQ6naYzsz+Z2dV9KLOr+h0AfBq4yN2fc/dWd18KfBg428xO62M5l5rZ0wmfu/xehXHzzexb4f+i1sweNbPyMK7QzO4I62lXmHfM/tRt0HF3vQbAC1gDnNFp2A3AHeH9VMCB3ITxlwKtwJVALlBEtPf8HqAAqCBKID/sbjmdyyXaA/9foBA4HKgETk+IpxE4B8gBvgMs6KY+I4GdwD+F2C4Mn0eF8bcC3+5m3mKgDTi1h/U1iqhBKSba0/wt8PuE8fOJjrLmhOXnhWEbgUOAEuC+hPU7C9gd1l0ecC3wJpCfsN4WAuND3ZYBV4RxpwAbwvsc4FXgB2EZhcCJ3dThSODYEN/UUObVYVwpsBn4QiijFDimi+/FwUAdcFLY5v8vfCfO6GaZtwK1CdP/D/B0GHc0sAkYEj6XEyXRMd2U5cDMHr6zVwBru5n3SeA7ffzfuDQhxt6+V/OBt8L2LAqfbwzj/hX4U/jO5IT1Xxb3//5AeOmIIPNtcvcfe7S31eDub7r7Y+7e5O6VRA3DyX0pyMwmAScC/+7uje7+CnAL0T9dh6fd/WGPrincDhzWTXH/AKx099tDbHcDbwDv60MoI4iOVjd3N4G7b3f3+9y93t1rgf/gnfW81d2XhuW3hGG3u/sSd98NfBX4qEUXyi8AHgrrrgX4b6KG5PiE8n7k7pvcfQdRg3J4F6EdTZQsvujuu8N6fLqL6fDoCGdBiG8N8POEOpwLbHH374cyat39+S6K+QjwoLs/5e5NoU7t3a234KGE6b8MHGdmk9x9IVANnB6m+xgw39239lDWS2HvepeZ7QKuSxhXTvfbcHMYv6/68r36lbuvcPcGolNQHduphWgHYqa7t4X1X7MfMQw6SgSZb33iBzMbbWb3mNlGM6sB7qDv/3DjgR2hYe2wFpiQ8HlLwvt6oNC6vjYxPsybqHNZ3dlJ1JiN624CMys2s5+HUyI1REc+w23vu5/WdzFr4rC1RHv/5Z3jdff2MG1PdR/aRfmTiPaCW7uLPaEOs8IprS2hDv/J29tqEtGebW/Gk1CnkOC29zJP4vR1wI5QDsCvgYvD+4uJkn1PjnD34R0v4MaEcVV0vw3HAVX29l1rdWZW18uyoG/fq+620+3AI8A9ZrbJzP7LzPL6sMxBT4kgc3TXTWzn4d8Jww519zKif2brQzkQnRYYaWalCcMmE51O2VebiC4IJ+pTWe5eDzxHdOqnO18AZhOdLikjOtUBvdd1Uqd4WogarL3iNTML0+5r3dcDk7tJjp3dRLQ3e0Cow5d4O/71wIw+lLGZhDqZWTHRXm9PEqcfSnS6ZVMYdAdwnpkdBhwE/L4PMXTnb8AkMzs6cWA48jwWeNzfvmut46aJ3vTne9Xi7t9w94OJjvTOJbqInfWUCDJHJdFe8vRepislOme8y8wmAF/sNH5rd2W4+3rgWeA74cLaocAngX2+gwl4GJhlZv9oZrlmdgHR+ewH+zj/tcClFt16OArAzA4zs3vC+FKggaieI4Gv97Hci83s4NBgfhP4XTjNdS/wD2Z2ethL/ALQRLQ+9sVCosb5RjMrCevxhG6mLQVqgDozO5DorqgODwJjzexqMysws1IzO6aLMn4HnGtmJ4YLzd+k9//rcxKm/xbwfNj2uPsG4AWivef7wumV/eLuK4CfAXea2bFmlmNmc4iuzfzV3f+6H8Xu9/fKzE41s3eFo8Yaop2AdN82PSApEWSIsJf8H8Az4Xzssd1M+g3gCKJzvQ8B93ca/x3gK6GMa7qY/0KiC5ebgAeAr7v7Y/sR73aiPa4vEJ2quBY4192r+jj/s0S3GZ4GrLLorqibiRoCgB8SncOvAhYAf+ljaLcTXTDdQnQR9rNhecuJjp5+HMp8H9HtvM19LLcj7rYw70yii9UbiK4/dOUa4B+JLt7+H/CbhHJqiS5cvy/EuhI4tYvlLQU+A9xFlIB2hmX25C6ixLmD6ILpRZ3G/xp4F72fFuqLfyO6znQH0Q7KX4gu4PZ0tNetfn6vxhIlzhqiC/NPhriynrnrwTQi8jYzO4mogZwarpXIIKcjAhHZI5wWuwq4RUkgeygRiAgAFv2AbhfRHT0/jDkcSSOdGhIRyXI6IhARyXIZ0UlZeXm5T506Ne4wREQyyosvvljl7hW9TZcRiWDq1KksWrQo7jBERDKKmXX+FXaXdGpIRCTLKRGIiGS5lCUCM/ulmW0zsyUJw75nZm9Y1Kf7Axb6yRcRkfik8ojgVuDsTsMeAw5x90OBFcD1KVy+iIj0QcoSgbs/RdSXSeKwRxO6510ATHzHjCIiklZxXiP4Z+DP3Y00s8vNbJGZLaqsrExjWCIi2SWWRBCeq9pKD90bu/vN7j7P3edVVPR6G6yIiOyntP+OwMwuIepG9nRX/xYikqXa2p3dza3UNrZS19hKXVMLNXvet1Lb2EJdYysfOmIiU8tLUhpLWhOBmZ0N/DtwcuhfX0Qk47S2tVPbGDXiNY0t1DS0hL/R59rQmNc1tlLb9Pbnjka/trGF3c29PxPHDOZOGZG5icDM7gZOAcrNbAPRgzCuBwqAx6InAbLA3a9IVQwiIl1pCQ155wa8pqGlU+Pe9TR9acRL8nMoLcxjaGEuQwtyKS3MZdywQkoL9h5WWpjL0IK86G9hLqUFuXvmK87LYcgQ63VZ/ZWyRODuF3Yx+BepWp6IZJf2dqe2sZVdDc3sqm9hV0MLu+qbqW5oiT7Xt7CroZnq+haqOzXm9b005EMMyoryKCvMo6wol9KCPKaWF4fPbw9/+3Nu9Lcoj6EFUSOfk4YGPFkyoq8hERm82tqdXfXNoSFvobqjYQ+Ne3XCuMTP1Q0t9HSVcWhBLsOK8hhenMewojxmVAx9u2EvTGi8OxrzhIa9JD+HcNYiKygRiEjSuDsNLW1sr2tmZ30z23c3s3N3Mzs6vRLH7eqhQTeDssKoMR9elMew4nymjCze6/Pw0NhHDX7+noY/L0c96PSVEoGIdMvdqWtqpaqumcraJrbXNe1pwLeHBr1zI9/U2vUTLnOHGCNK8hlZnM/IknwOGlvGyJL8MCyPESX5YQ/+7ca9tDAvo06xZColApEs09G4V9Y2UVXXTFVdE1V1TeFzE5W1zXt97q5hLy3IjRrxknzGlBVy0LioYR8ZGvuOcR2vssLcrDrdkkmUCEQGifZ2Z0d9M1trGsOriS3VjWyrbaSytpnKuiaqemjchxiMLCmgfGg+FaUFTCsv2fO+fOjbr1FDo9MvBbk5MdRSUkGJQGSA69iD31rTtKeR31LTyLbwueP9ttpGWtr2PtluBqNK8ikfWkBFaQEzyksoL40a+45hHQ38yJJ8nYbJUkoEIjFrbGljS3Ujm3Y1sHFXA5t2Re83VTewaVcDm6sbu7zdsbQwlzFlhYwpK+CY6SMZU1bI2PA5Gl5IRWmBLppKr5QIRFLI3dlZ38L6HfXvaOg3VzewcVcjVXVN75ivorSA8cOLmDWmlJNmVTC2rJCxwwoZXdrxt4CSAv37SnLomyTST02tbWzc2cC6HfWs31HPuj2vBtbvqKeuqXWv6Yvycpgwoojxw4s4aFwZ44cXhVchE4YXMXZYoc6/S1opEYj0QV1TK2uqdvNWZR3rtr/d2K/fUc/mmsa97oMvyB3CpJHFTB5ZzDHTRjJpZDGTRhQxYUQRE4YXMawoT3fPyICiRCAStLa1s2FnA6tDg7+qajerK3ezqqqOrTV7n74ZU1bA5JHFHDtjFJNDo9/xKh9akJb+YUSSRYlAsk5TaxurKnezYmtteNWxqrKOdTvq97rrZnhxHtPLSzhxZgXTK0qYXl7C9IqhTBlVTGGeTt3I4KFEIINWa1s7a7bXs3JrLcsTGv3VVbtpa48a/NwhxtTyEmaOHsqZc8YyrbyEGRUlTCsfysiS/JhrIJIeSgQyKNQ3t7Jscy1LN1WzdGMNSzZVs3JrHc1t0Q+nzGDKyGJmjSnl7DljmTW2lNljSplWXkJ+rm6vlOymRCAZp6axhSUbq3l9Uw1LNlazZFMNqyrrCDv5jCjO45AJw7j0hKnMHlPKrDGlzBw9lKJ8nc4R6YoSgQxobe3Oym21vLxuFy+v28nL63bxZmXdnrt0xpYVcsiEMv7hXeOYM76MQyYMY9ywQt2VI7IPlAhkQKlpbGHRmh28uDZq9F9dv2vP06BGFOcxd/II3n/YeA6dNJw548soH1oQc8QimU+JQGJVXd/CwjU7WLBqO8+v3s7rm2po9+gi7kHjyvjwkROZO3k4cyeNYMqoYu3pi6SAEoGkVUNzGwtWb+fvK6pYsGo7y7bU4A75uUM4YvJwrjztAI6ZPpK5k0bonL5ImigRSEq5Oyu31fHUikqeXFHJ86t30NzaTkHuEI6cMoLPnTGLY6aN5LBJw3VvvkhMlAgk6Zpa23j2re08unQr85dvY3N1IwAHjB7KPx07hZNnVXD0tJFq+EUGCCUCSYq6plbmL9/GI0u3Mv+NbdQ2tVKSn8O7D6jgqtMrOGlWBeOHF8Udpoh0IWWJwMx+CZwLbHP3Q8Kw84EbgIOAo919UaqWL6nX0NzGX5dt5Q+vbOSplVU0t7YzqiSfc941jrMOGcPxM8q11y+SAVJ5RHAr8BPgtoRhS4APAT9P4XIlhdranWffquKBlzfyyJIt7G5uY2xZIRcfM4Wz5oxh3tSResqVSIZJWSJw96fMbGqnYcsA3QKYgdZu383dC9dz30sbqKxtorQgl3MPHc95c8dzzLRRavxFMtiAvUZgZpcDlwNMnjw55miyU0tbO399fSt3LVzH31dWkTPEOHX2aD58xAROPXC0TvuIDBIDNhG4+83AzQDz5s3zXiaXJNq5u5k7FqzltgVrqaxtYvywQj7/nll8dN4kxg4rjDs8EUmyAZsIJP3WVO3mF0+v5rcvrqexpZ1TZlfw8eOmcPKs0Tr1IzKIKREISzdV8+PH3+SR17eQN2QIH5g7nsvePZ1ZY0rjDk1E0iCVt4/eDZwClJvZBuDrwA7gx0AF8JCZveLuZ6UqBunZ8i21/PCvK/jzki2UFuby6VNmcMlxUxldptM/ItkklXcNXdjNqAdStUzpm9VVu/n+o8t5aPFmhubnctXpB/DPJ05jWFFe3KGJSAx0aiiLVNe38KO/reS259aQnzOEz5wyk8vePY3hxXoko0g2UyLIAq1t7dy1cB0/eGwFuxpa+NhRk/j8e2ZTUaq+/EVEiWDQe3ndTq6/fzFvbKnluOmj+Mq5BzFn/LC4wxKRAUSJYJCqbWzhvx9Zzm0L1jKmtJCfXXwkZ80Zo191i8g7KBEMQvOXb+O6+xaztbaRS46byhfOnEVpoS4Ei0jXlAgGkYbmNv7z4WXcvmAts8YM5aaLj2fu5BFxhyUiA5wSwSDx2oZdXP2bV1hVuZvLTpzGNWfNVl9AItInSgQZzt359bNr+PZDy6goLeCuy47h+JnlcYclIhlEiSCD7W5q5br7F/OnVzdxxkGj+f75hzOsWNcCRGTfKBFkqDe31XHFHS+yqrKOL541m0+dPIMh6hhORPaDEkEGenplFZ+680Xyc4Zw+yeP4QSdChKRflAiyDC/eWEdX35gCTMqhvKLS+cxcURx3CGJSIZTIsgQ7e3O9x5dzk3z3+LdB5Tz04uOoEy/DRCRJFAiyACtbe1ce99r3P/SRv7xmMl84/1zyMsZEndYIjJIKBEMcM2t7Vx1z8v8eckWPnfGLD57+kx1EyEiSaVEMIA1trRxxR0vMn95JV8992A+eeK0uEMSkUFIiWCAampt419vf5GnVlby3Q+/iwuOmhx3SCIySCkRDECtbe189u6XeXJFJf/14UP56FGT4g5JRAYxXXEcYNrbnWt++yqPLN3K1993sJKAiKScEsEA862HXuf3r2zimjNn8YkTdE1ARFJPiWAAue25NfzqmTVcevxUPnPqzLjDEZEsoUQwQDzxxjZu+ONSzjhoNF8992DdIioiaZOyRGBmvzSzbWa2JGHYSDN7zMxWhr96agqwpmo3n737ZQ4cW8b/fGwuOeo8TkTSKJVHBLcCZ3cadh3wuLsfADwePme1xpY2PnXnS+TkGP93yTxKCnQjl4ikV8oSgbs/BezoNPg84Nfh/a+BD6Rq+Znihj8uZdnmGn5wweFMGF4UdzgikoXSfY1gjLtvBgh/R3c3oZldbmaLzGxRZWVl2gJMpz++uol7XljPZ06dwamzu10VIiIpNWAvFrv7ze4+z93nVVRUxB1O0m2raeSrv1/C4ZOG87kzZsUdjohksXQngq1mNg4g/N2W5uUPCO7Olx5YTGNLG9//6GHkqidREYlRulugPwKXhPeXAH9I8/IHhIcXb+Gvy7bxxbNmM6NiaNzhiEiWS+Xto3cDzwGzzWyDmX0SuBF4j5mtBN4TPmeV+uZW/uOh1zloXBmXHj817nBERFLX6Zy7X9jNqNNTtcxM8L9PvMWm6kb+58K5OiUkIgOCWqI02lLdyM1/X8UHDh/PUVNHxh2OiAigRJBWN81/k/Z25wtnzo47FBGRPZQI0mRLdSN3L1zP+fMmMmlkcdzhiIjsoUSQJjfNf5N2dz59inoVFZGBRYkgDarrW/jNovV8cO4EHQ2IyICjRJAGv31xPY0t7Vx6wtS4QxEReQclghRrb3duX7CWo6aOYM74YXGHIyLyDkoEKbZo7U7Wbq/nH4+ZHHcoIiJdUiJIsd+/spGivBzOPHhs3KGIiHRJiSCFmlvbeXjxZs6cM0YPnBGRAUuJIIWefauKXfUtvP+w8XGHIiLSLSWCFHpyRSUFuUM4YWZ53KGIiHRLiSCFnlpRybHTR1GYlxN3KCIi3VIiSJENO+t5q3I3J80afE9XE5HBRYkgRf6+sgqAk5UIRGSAUyJIkZfW7mRUST4zKkriDkVEpEdKBCmyeGM1h04chpnFHYqISI+UCFKgobmNFVtredfE4XGHIiLSKyWCFHh9czXtDodOUN9CIjLwKRGkwGsbqgF410QlAhEZ+JQIUmDxhmrGlBUwpqww7lBERHoVSyIws6vMbImZLTWzq+OIIZVWbqtj1pjSuMMQEemTtCcCMzsE+BfgaOAw4FwzOyDdcaSKu7OmajfTy3XbqIhkhjiOCA4CFrh7vbu3Ak8CH4whjpTYvruZ2qZWpioRiEiGiCMRLAFOMrNRZlYMnANM6jyRmV1uZovMbFFlZWXag9xfq6t2AygRiEjGSHsicPdlwHeBx4C/AK8CrV1Md7O7z3P3eRUVmdNNw4ad9QBMGqGH1ItIZojlYrG7/8Ldj3D3k4AdwMo44kiFzdWNAIwbpjuGRCQzxPLYLDMb7e7bzGwy8CHguDjiSIWt1Y2UFebqiWQikjHiaq3uM7NRQAvwGXffGVMcSbe5upFxw4riDkNEpM9iSQTu/u44lpsOW2oaGavTQiKSQfTL4iTbUt3IWP2iWEQySI9HBGa2GPCuRgHu7oemJKoM5e7sqm9hREl+3KGIiPRZb6eGzg1/DXiI6J5/6UZjSzvNbe0MK8qLOxQRkT7rMRG4+9qO92bWlPhZ3qmmsQWAsiLdMSQimUPXCJKouiFKBDoiEJFM0ts1giMSPhZ1+oy7v5SSqDKUEoGIZKLezmF8P+H9FuC/w3sjuoh8WiqCylTV9UoEIpJ5ertGcCqAmRUBnwZOJEoAfwduSnl0GUZHBCKSifp6VfPXQA3wo/D5QuA24KOpCCpTdSSCskIlAhHJHH1NBLPd/bCEz0+Y2aupCCiTvX3XkBKBiGSOvt419LKZHdvxwcyOAZ5JTUiZq7qhhdKCXHKGWNyhiIj0WV+PCI4BPm5m68LnycCyjl8e6xfGkeqGFh0NiEjG6WsiODulUQwSNQ0tulAsIhmnT4lAvyjum+iIQL8qFpHMol8WJ1G1jghEJAMpESRRTUOrEoGIZBwlgiTSEYGIZCIlgiRpbm2noaVNPyYTkYyjRJAke7qXKFYiEJHMokSQJOpnSEQylRJBkqh7CRHJVEoESaIjAhHJVLEkAjP7nJktNbMlZna3mRXGEUcy1ajnURHJUGlPBGY2AfgsMM/dDwFygI+lO45k0xGBiGSquE4N5RI9+jIXKAY2xRRH0ujpZCKSqdKeCNx9I9EjL9cBm4Fqd3+083RmdrmZLTKzRZWVlekOc59VN7RQlJdDfq4uu4hIZonj1NAI4DxgGjAeKDGziztP5+43u/s8d59XUVGR7jD3WU2jflUsIpkpjt3XM4DV7l7p7i3A/cDxMcSRVOp5VEQyVRyJYB1wrJkVm5kBpwPLYogjqXbV64hARDJTHNcIngd+B7wELA4x3JzuOJKtsq6JitKCuMMQEdlnsVzZdPevu/uB7n6Iu/+TuzfFEUcyVdY0Mbo0438OISJZSLe4JEFDcxu1Ta06IhCRjKREkASVtdEBzWglAhHJQEoESbCtthGA0WU6NSQimUeJIAm26YhARDKYEkESbNzZAMD44UUxRyIisu+UCJJg7Y7dDCvK0+8IRCQjKREkwbodDUweWRx3GCIi+0WJIAnW76hn8iglAhHJTEoE/dTS1s6GnfU6IhCRjKVE0E9rqnbT0ubMGjM07lBERPaLEkE/Ld9aC8ABo0tjjkREZP8oEfTTiq11DDGYOVpHBCKSmZQI+mnpxmqmlZdQmJcTdygiIvtFiaAf3J3XNlZz2MThcYciIrLflAj6YUtNI5W1TRw6cVjcoYiI7Dclgn5YuHoHAEdOGRlzJCIi+0+JoB+efXM7ZYW5HDy+LO5QRET2mxJBPzy7qopjp48iZ4jFHYqIyH5TIthP63fUs35HA8fPGBV3KCIi/aJEsJ+eebMKgONmlMcciYhI/ygR7KdHX9/KxBFF6lpCRDKeEsF+qG1s4emVVZw9Zyxmuj4gIpkt7YnAzGab2SsJrxozuzrdcfTHE8sraW5r5+xDxsYdiohIv+Wme4Huvhw4HMDMcoCNwAPpjqM/HnptE+VDC5g7eUTcoYiI9Fvcp4ZOB95y97Uxx9Fn2+uaeHzZNj5w+HjdNioig0LcieBjwN1djTCzy81skZktqqysTHNY3fv9K5tobXfOnzcp7lBERJIitkRgZvnA+4HfdjXe3W9293nuPq+ioiK9wXXD3bn3hfUcOnEYs8fq+QMiMjjEeUTwXuAld98aYwz75Lm3trN8ay0XHTM57lBERJImzkRwId2cFhqofvnMakaV5HPe4RPiDkVEJGliSQRmVgy8B7g/juXvj9VVu3n8jW1cdOwUPYRGRAaVtN8+CuDu9UBGddLzq2dWkzvEuPhYnRYSkcEl7ruGMsKW6kbueWE9H5o7kdGlhXGHIyKSVEoEffC/89+kvd35t9Nmxh2KiEjSKRH0YtOuBu5ZuJ7z501i0sjiuMMREUk6JYJe/PSJN3F0NCAig5cSQQ/Wba/n3kXr+ei8SUwYXhR3OCIiKaFE0IMb/7KM3CFDuPK0A+IORUQkZZQIurFw9Q4eXryFK06ewdhhulNIRAYvJYIutLc733rwdcYNK+Tyk6bHHY6ISEopEXTh/pc3snhjNdeePZuifP2KWEQGNyWCTuqbW/neI29w2MRhnHeY+hQSkcFPiaCTnz25iq01TXztfQczRA+eEZEsoESQYHN1Azc/9RbnHjqOI6eMjDscEZG0UCJI8F9/WU67w3XvPTDuUERE0kaJIHhl/S4eeHkjl504jYkj1JWEiGQPJQKiR1B++8HXKR+az6dPVVcSIpJdlAiAR5ZuYdHanXz+PbMZWhDLIxpERGKT9YmgubWdG//8BrPGDOWj8ybGHY6ISNplfSK4Y8Fa1myv5/pzDiI3J+tXh4hkoaxu+arrW/jR31by7gPKOWVWRdzhiIjEIqsTwU/nv0l1QwvXv/cgzPTjMRHJTlmbCLbVNnLbc2v44OETOHh8WdzhiIjEJpZEYGbDzex3ZvaGmS0zs+PSHcPPn1xFS5tz5el61oCIZLe47pX8H+Av7v4RM8sH0voLrm21jdyxYC0fOHwC08pL0rloEZEBJ+2JwMzKgJOASwHcvRloTmcMP3/ruZYoAAAOYElEQVRyFa3tzpV6DrGISCynhqYDlcCvzOxlM7vFzN6xW25ml5vZIjNbVFlZmbSFV9U1cceCtXxw7gSm6mhARCSWRJALHAHc5O5zgd3AdZ0ncveb3X2eu8+rqEjerZ13Pb+OptZ2rjh5RtLKFBHJZHEkgg3ABnd/Pnz+HVFiSLnm1nZuX7CWk2dVMHP00HQsUkRkwEt7InD3LcB6M5sdBp0OvJ6OZT+8eDOVtU184oSp6ViciEhGiOuuoSuBO8MdQ6uAT6Rjob96dg3Ty0s46QD9ilhEpEMsicDdXwHmpXOZy7fU8ur6XXztXD2CUkQkUdb8svi+lzaQO8Q47/DxcYciIjKgZEUiaG1r5/6XNnLagaMZNbQg7nBERAaUrEgEz6/eQVVdEx86YkLcoYiIDDhZkQgeWbqFwrwhnDxrdNyhiIgMOIM+EbS3O48u3crJsyooys+JOxwRkQFn0CeC1zZWs6WmkbPmjI07FBGRAWnQJ4KnVlRiBqfO1mkhEZGuDPpE8Nxb2zl4XBkjSvLjDkVEZEAa1ImgsaWNF9ft5Njpo+IORURkwBrUieDldbtobm3nOCUCEZFuDepE8Nyq7QwxOHr6yLhDEREZsAZ1IpgwvJCPHDmRssK8uEMRERmw4up9NC0uOGoyFxw1Oe4wREQGtEF9RCAiIr1TIhARyXJKBCIiWU6JQEQkyykRiIhkOSUCEZEsp0QgIpLllAhERLKcuXvcMfTKzCqBtfswSzlQlaJwBrJsrHc21hmys97ZWGfoX72nuHtFbxNlRCLYV2a2yN3nxR1HumVjvbOxzpCd9c7GOkN66q1TQyIiWU6JQEQkyw3WRHBz3AHEJBvrnY11huysdzbWGdJQ70F5jUBERPpusB4RiIhIHykRiIhkuYxIBGZ2tpktN7M3zey6LsYXmNlvwvjnzWxqwrjrw/DlZnZWX8uMW7LrbGaTzOwJM1tmZkvN7Kr01abvUrGtw7gcM3vZzB5MfS32TYq+38PN7Hdm9kbY5selpzZ9l6J6fy58v5eY2d1mVpie2vTN/tbZzEaF/986M/tJp3mONLPFYZ4fmZntc2DuPqBfQA7wFjAdyAdeBQ7uNM2ngZ+F9x8DfhPeHxymLwCmhXJy+lLmIKzzOOCIME0psGIg1TlV9U6Y7/PAXcCDcdczHXUGfg1cFt7nA8PjrmsavuMTgNVAUZjuXuDSuOuapDqXACcCVwA/6TTPQuA4wIA/A+/d19gy4YjgaOBNd1/l7s3APcB5naY5j+iLD/A74PSQFc8D7nH3JndfDbwZyutLmXFKep3dfbO7vwTg7rXAMqJ/nIEkFdsaM5sI/ANwSxrqsK+SXmczKwNOAn4B4O7N7r4rDXXZFynZ1kSP3y0ys1ygGNiU4nrsi/2us7vvdvengcbEic1sHFDm7s95lBVuAz6wr4FlQiKYAKxP+LyBdzZge6Zx91agGhjVw7x9KTNOqajzHuFwcy7wfBJjToZU1fuHwLVAe/JD7rdU1Hk6UAn8KpwOu8XMSlIT/n5Ler3dfSPw38A6YDNQ7e6PpiT6/dOfOvdU5oZeyuxVJiSCrs53db7ntbtp9nX4QJGKOkczmQ0F7gOudvea/Y4wNZJebzM7F9jm7i/2N7gUScW2zgWOAG5y97nAbmCgXQdLxbYeQbRHPQ0YD5SY2cX9ijK5+lPn/pTZq0xIBBuASQmfJ/LOw70904RDwmHAjh7m7UuZcUpFnTGzPKIkcKe735+SyPsnFfU+AXi/ma0hOhQ/zczuSEXw+ylV3+8N7t5xxPc7osQwkKSi3mcAq9290t1bgPuB41MS/f7pT517KnNiL2X2Lu4LKH24wJILrCLK8h0XWOZ0muYz7H2B5d7wfg57X1RaRXTBptcyB2Gdjej84Q/jrl86691p3lMYeBeLU1Jn4O/A7PD+BuB7cdc1Dd/xY4ClRNcGjOhc+5Vx1zUZdU4YfynvvFj8AnAsb18sPmefY4t75fRxBZ5DdJfLW8CXw7BvAu8P7wuB3xJdNFoITE+Y98thvuUkXE3vqsyB9Ep2nYnuOHDgNeCV8NrnL0ym1btT2acwwBJBCr/fhwOLwvb+PTAi7nqmqd7fAN4AlgC3AwVx1zOJdV5DdHRQR3QkcHAYPi/U9y3gJ4QeI/blpS4mRESyXCZcIxARkRRSIhARyXJKBCIiWU6JQEQkyykRiIhkOSWCFDAzN7PvJ3y+xsxuSFLZt5rZR5JRVi/LOT/0WvlEksu9wcyuSWaZncr/Xuh98ntJLHPPOg/dNRycrLITlrHGzMr7MX9dMuPpw/Lmm1mfH6huZt80szNSGE+v9Tezq82sOOHzw2Y2PFUxZZLcuAMYpJqAD5nZd9y9Ku5gOphZjru39XHyTwKfdvekJoI0+Fegwt2bUlG4u1+WinJTYR+3d0q5+9fijgG4GrgDqAdw93PiDWfg0BFBarQSPWf0c51HdN6j79iTMbNTzOxJM7vXzFaY2Y1mdpGZLQx9jc9IKOYMM/t7mO7cMH9O2Bt+wcxeM7N/TSj3CTO7C1jcRTwXhvKXmNl3w7CvEf0A7Wed96z7GqeZTTGzx0Msj5vZ5C6WPcPM/mJmL4b6HBiGnx/iedXMnupiPgt1XRKWeUEY/kei7nqf7xiWMM/JZvZKeL1sZqU9lGNm9hMze93MHgJGJ5SzZ084cS/UzD5iZrcmbOObwnpfFZb9y3CEdWvn+iS40sxeCrEcaGZDzGylmVWEcodY1Od8uZlNM7Pnwvb+Vqfts9f2NrPPhzouMbOrw7CjwrYpNLMSi46iDjGz283svITy7jSz95tZkZndE+b5DVCUMM2ZIZaXzOy3FvVn1XmbJR5VrTGzbyTWtYvp54Tv1CthmQd0V5dO851iCc+cCNvxUjP7LFH/Q09YOMq1hKOwbtbR1LDN/i+sn0fNrKjzMgeFuH9pNxhfRL/8KyP6JeAw4BrghjDuVuAjidOGv6cAu4ieG1AAbAS+EcZdRegaIsz/F6IkfgDRLwwLgcuBr4RpCoh+VTotlLsbmNZFnOOJemqsIDo6/BvwgTBuPjCvi3n6GuefgEvC+38Gfh/e3wBcE94/DhwQ3h8D/C28X0zUmyR00Y8+8GHgMaJuBcaEOoxLXJ9dzPMn4ITwfmiob5flAB9KGD4+1PcjnddL4rKAjwC3Jmyje4h+8n8eUAO8K2yzF4HDu4hvDaE7BKI+6W8J779O1EEgwJnAfeH9H4GPh/efYe/v0Z7tDRwZ1mdJqPdSYG4Y922i3jp/Clwfhp2csK2GEfXvn0v0PIdfhuGHEu3szAPKgaeAkjDu34GvdVG/WxPWYZd17TT9j4GLwvt8osTTU10S6/9gQjk/ITyTICy3vNM6L++uXGBqqOfhYfp7gYvjbl9S8dIRQYp41LPnbcBn92G2Fzx6bkAT0c/FO7rQXUz0pexwr7u3u/tKor5LDiRqJD5uZq8QdS89iihRACz0qN/2zo4C5nvUSVcrcCdRP/bJiPM4ogfBQPRT/xMTCwh7jccDvw0x/5yoEQZ4BrjVzP6FqDHu7ETgbndvc/etwJOhLj15Bvh/Yc9weKhvd+WclDB8E1GC3Fd/8qj1WAxsdffF7t5O1MhM7Waejo4AX0yY5pfAx8P7fwZ+Fd6fANwd3t/eqZzE7X0i8IBH/dnXhWW8O4z7JvAeogb9vwDc/UlgppmNBi4kSjytROvkjjDNa0RdV0DUx83BwDNhO14CTOl+tfRY10TPAV8ys38Hprh7Qy916Y+eyl3t7q/0EmvG0zWC1Poh8BJv//NCtIcxBKJTEER7Ox0Sz2u3J3xuZ+9t1blfECfa+7zS3R9JHGFmpxDtIXZl3x9pt29xdo4x0RBgl7sf/o4J3a8ws2OIHibzipkd7u7bEybZ57jd/cZwmuccYIFFFy57Kqcvfa8kTtP5kYiJ66Tz+upuHXVM19YxjbuvN7OtZnYa0VHTRX2IMXF791THkUR7wHkh/o75bg/L+RhR8ulpeQY85u4X9rCcrryjronc/S4ze57oO/CImV1G37b7nv+voC+Pquyp3MRt10bCKbHBREcEKeTuO4gOJz+ZMHgN0aEoRKcN8vaj6PPD+eIZRA8hWQ48AnzKoq6mMbNZ1vvDSJ4HTg7nnHOI9gCf3I94uvIsUUMCUaPydOLIcMS02szOD/GamR0W3s9w9+c9usBYxd5d90J0KuICi66LVBDtrS7sKZhQ5mJ3/y7RabMDeyjnKeBjYfg44NRuit1qZgeZ2RDgg72ukf13C9He+L3+9sXfZ9h7/XbnKeADZlYcvg8fJOqZFKLrWF8lOhL8bsI8txJdWMXdlyaUcxGAmR1CdHoIYAFwgpnNDOOKzWzWftRxL2Y2HVjl7j8iOg12aC916bAWONiiZ/8OA05PGFdL9JjWzvpS7qCmI4LU+z7wbwmf/w/4g5ktJDpH3t3eek+WEzXYY4Ar3L3RzG4hOmx9KRxpVNLLI+vcfbOZXQ88QbRX9LC7/2E/4unKZ4FfmtkXQyyf6GKai4CbzOwrRAnxHqKueb8XLg4a0Tp6tdN8DxCdenqVaC/1Wnff0ks8V5vZqUR7da8Tddfb3FU5ZvYAcBrRaZ0VdJ8crwMeJHqi1BKivetU+CPRUWXikeVVwF1mdhXRMya65O4vWXSBuiNR3uLuL5vZx4HWsOedAzxrZqe5+9/cfauZLSPqtbTDTURPPOvovXZhKL/SzC4F7jazgjDtV4jWW39cAFxsZi3AFuCb7r6jq7p0qu96M7uX6NTVSiBx/M3An81ss7ufmjBPd+toaj/rkDHU+6jIAGfRXUo/cPdknA/vy/KKiZLgEe5enY5lSrx0akhkADOz64j2+K9P0/LOIOrP/8dKAtlDRwQiIllORwQiIllOiUBEJMspEYiIZDklAhGRLKdEICKS5f4/YcCRLpuWruQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "na_moles = np.linspace(1e-4, 10.e-3, 500)\n",
    "sol_volume = 1. # Liter\n",
    "carbo= Acid(pKa=[6.4,10.3], charge=0, conc=1.e-3) # pKa given at 25°C\n",
    "phs = []\n",
    "for mol in na_moles:\n",
    "     na = Neutral(charge=1, conc=mol/sol_volume)\n",
    "     system = System(carbo, na)\n",
    "     system.pHsolve(guess_est=True)\n",
    "     phs.append(system.pH)\n",
    "plt.plot(na_moles, phs)\n",
    "\n",
    "plt.title('Titration of Carbonic acid by HO- ions')\n",
    "plt.xlabel('Number of moles of sodium hydroxyde in solution')\n",
    "plt.ylabel('pH')\n",
    "plt.savefig(\"titration.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}