{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running the calculation\n",
    "**Running the code under Gouy-Chapman conditions takes approximately 19 minutes (iMac with 4 Ghz i7 processor).**\n",
    "\n",
    "This is due to the Poisson-Boltzmann solver being run multiple times to minimise the difference between input and output mole fractions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyscses.defect_species import DefectSpecies\n",
    "from pyscses.set_of_sites import SetOfSites\n",
    "from pyscses.constants import boltzmann_eV\n",
    "from pyscses.calculation import Calculation, calculate_activation_energies\n",
    "from pyscses.set_up_calculation import calculate_grid_offsets\n",
    "from pyscses.grid import Grid\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "boundary_conditions = 'dirichlet'\n",
    "site_charges = False\n",
    "systems = 'gouy-chapman'\n",
    "core_models = False\n",
    "site_models = 'continuum'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 0.0005\n",
    "\n",
    "conv = 1e-8\n",
    "grid_x_min = -2.0e-8\n",
    "grid_x_max = +2.0e-8\n",
    "bulk_x_min = -2.0e-8\n",
    "bulk_x_max = -1.0e-8\n",
    "\n",
    "dielectric = 1\n",
    "\n",
    "index = 111\n",
    "\n",
    "b = 5e-9\n",
    "c = 5e-9\n",
    "\n",
    "temp = [773.15]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "valence = [ +1.0, -1.0 ]\n",
    "site_labels = ['site_1', 'site_2']\n",
    "defect_labels = ['defect_1', 'defect_2']\n",
    "mole_fractions = [ [ 0.2, 0.2 ] ]\n",
    "initial_guess = [ [ 0.2, 0.2 ] ]\n",
    "mobilities = [1.0, 1.0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = '../input_data/example_data_2_one_seg_energies.txt'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ```limits``` are the widths between the midpoint of the final site and an 'imaginary' site extending beyond the calculation unit. These are calculated to use in delta_x calculations where the width of each site is required to calculate volumes or densities. The ```laplacian_limits``` are the widths between the end sites.\n",
    "\n",
    "```\n",
    "limits, laplacian_limits = calculate_grid_offsets( data, grid_x_min, grid_x_max, 'single' )\n",
    "```\n",
    "\n",
    "If the calculation is being run over a number of different conditions, the mole fractions are looped over. \n",
    "\n",
    "```\n",
    "for m in mole_fractions:\n",
    "```\n",
    "\n",
    "If the calculation is being run over a number of different conditions, the temperatures are looped over.\n",
    "\n",
    "```\n",
    "for t in temp:\n",
    "```\n",
    "\n",
    "Next the different defect species are defined using the ```Defect_Species``` class. This collects and stores the information about each defect, the defect label, the defect valence and the defect mole fraction.\n",
    "\n",
    "```\n",
    "defect_species = { l : Defect_Species( l, v, m ) for l, v, m in zip( defect_labels, valence, m) }\n",
    "```\n",
    "\n",
    "The input data file and system information are passed into ```set_of_sites_from_input_data```, which formats and splits wach line of the input file to create ```Site``` objects for the individual sites, which are then grouped together into a ```SetOfSites``` object. \n",
    "\n",
    "```\n",
    "all_sites = SetOfSites.set_of_sites_from_input_data( data, [grid_x_min, grid_x_max], defect_species, site_charges, core_models, t )\n",
    "```\n",
    "\n",
    "The code checks whether the model being used is continuum, if so, the sites are passed into ```form_continuum_sites``` which interpolates the defect segregation energies onto a regular grid.\n",
    "\n",
    "```\n",
    "if site_models == 'continuum':\n",
    "    all_sites, limits = SetOfSites.form_continuum_sites( all_sites, grid_x_min, grid_x_max, 1000, b, c, defect_species, laplacian_limits, site_labels, defect_labels )\n",
    "```\n",
    "\n",
    "ace charge models typically consider dopant ions as either mobile or immobile. If the dopant ions are considered immobile, the model follows a Mott-Schottky approximation. If the dopant ions are considered mobile, the model follows a Gouy-Chapman approximation. The code checks whether a Mott-Schottky or Gouy-Chapman approximation is being applied and fixes the mole fractions of the defects accordingly. \n",
    "\n",
    "```\n",
    "if systems == 'mott-schottky':\n",
    "    for site in all_sites.subset( 'site_2' ):\n",
    "        site.defect_with_label('defect_2').fixed = True\n",
    "if systems == 'gouy-chapman':\n",
    "    for site in all_sites.subset( 'site_2' ):\n",
    "        site.defect_with_label('defect_2').fixed = False\n",
    "```\n",
    "\n",
    "A grid is created using the ```Grid``` class. This includes the array of $x$ coordinates, b/c (the width and height of the cell used in the atomistic simulation) and the full set of sites ```all_sites```\n",
    "\n",
    "```\n",
    "grid = Grid.grid_from_set_of_sites( all_sites, limits, laplacian_limits, b, c )\n",
    "```\n",
    "\n",
    "Following this, a Calculation object is created. This class contains the functions which calculate the space charge properties if the system.\n",
    "\n",
    "```\n",
    "c_o = Calculation( grid, bulk_x_min, bulk_x_max, alpha, conv, t, boundary_conditions, initial_guess )\n",
    "```\n",
    "\n",
    "Subgrids are created for each of the defect species.\n",
    "\n",
    "```\n",
    "c_o.form_subgrids( site_labels )\n",
    "```\n",
    "\n",
    "If the system follows Gouy-Chapman conditions, and periodic boundary conditions are being enforced, a mole fraction correction is required to ensure the output bulk mole fraction is equivalent to the input bulk mole fraction. This function runs the Poisson-Boltzmann calculation and minimises the difference between the target and output mole fractions. This is not required when running the calculation using Mott-Schottky conditions. \n",
    "\n",
    "```\n",
    "c_o.mole_fraction_correction( m, system )\n",
    "```\n",
    "\n",
    "Next, the Poisson-Boltzmann equation is solved for the given system. The potential and charge density are calculated. This uses a function ```calculation``` which takes as arguments ``` grid, conv, temp, alpha and boundary_conditions ```. The ```calculation``` function creates numpy arrays the same length as the grid including zeros for the potential (```phi```) and the charge density (```rho```). \n",
    "\n",
    " **Solving the Poisson-Boltzmann equation**\n",
    "\n",
    "1. The invertible matrix $A$ is created from the distances between each of the points on the grid. Initially the full matrix is created using ```scipy.sparse.diags```  and is then converted into a sparse tridiagonal matrix using ```scipy.sparse.csc_matrix``` and the calculation follows the finite difference approximation method defined above. \n",
    "2. The convergence is initialised as 1 and while it is larger than ```conv``` defined in the notebook the calculation will iterate. \n",
    "3. Using the ```phi``` array of zeros, the charge density is calculated at each site. Using the calculated value for the charge density at each site $\\Phi = A^{-1}b$ is solved using ```scipy.linalg.spsolve```. This returns ```predicted_phi``` at each site, which is damped using $\\alpha$ to become the new ```phi``` array. The convergence is updated using $\\frac{ \\sum_i ( ( \\Phi_{predicted} - \\Phi )^2 )}{L_{grid.x}}$ where $L_{grid.x}$ is the length of the grid.\n",
    "\n",
    "```\n",
    "c_o.solve(system)\n",
    "```\n",
    "\n",
    "Using the calculated equilibrium electrostatic potentials at each site, the defect mole fractions can be calculated. \n",
    "```\n",
    "c_o.mole_fractions()\n",
    "```\n",
    "\n",
    "The resistivity ratio can be calculated from the previously calculated defect concentrations. ```c_o.calculate_resistivity_ratio``` takes two arguments. Firstly, whether the space charge potential is positive or negative and secondly a limit for the value of the potential that the site must exceed to be considered as the space charge region.\n",
    "\n",
    "```\n",
    "c_o.calculate_resistivity_ratio( 'positive', 2e-2 )\n",
    "```\n",
    "\n",
    "If the calculation is run over a number of temperatures, and the resistivity ratios stored in a list the activation energy can be calculates.\n",
    "\n",
    "```\n",
    "c_o.calculate_activation_energies( resistivity_ratios, temp )\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- 1132.7030682563782 seconds ---\n"
     ]
    }
   ],
   "source": [
    "limits, laplacian_limits = calculate_grid_offsets( data, grid_x_min, grid_x_max, 'single' )\n",
    "\n",
    "for m in mole_fractions:\n",
    "    for t in temp:\n",
    "    \n",
    "        defect_species = { l : DefectSpecies( l, v, m, mob ) for l, v, m, mob in zip( defect_labels, valence, m, mobilities ) }\n",
    "\n",
    "        all_sites = SetOfSites.set_of_sites_from_input_data( data, [grid_x_min, grid_x_max], defect_species, site_charges, core_models, t )\n",
    "        if site_models == 'continuum':\n",
    "            all_sites, limits = SetOfSites.form_continuum_sites( all_sites, grid_x_min, grid_x_max, 1000, b, c, defect_species, laplacian_limits, site_labels, defect_labels )\n",
    "        if systems == 'mott-schottky':\n",
    "            for site in all_sites.subset( 'site_2' ):\n",
    "                site.defect_with_label('defect_2').fixed = True\n",
    "        if systems == 'gouy-chapman':\n",
    "            for site in all_sites.subset( 'site_2' ):\n",
    "                site.defect_with_label('defect_2').fixed = False\n",
    "        grid = Grid.grid_from_set_of_sites( all_sites, limits, laplacian_limits, b, c )\n",
    "        \n",
    "        c_o = Calculation( grid, bulk_x_min, bulk_x_max, alpha, conv, dielectric, t, boundary_conditions )\n",
    "        c_o.form_subgrids( site_labels )\n",
    "        if systems == 'gouy-chapman':\n",
    "            c_o.mole_fraction_correction( m, systems, initial_guess )\n",
    "        c_o.solve(systems)\n",
    "        c_o.mole_fractions()\n",
    "        c_o.calculate_resistivity_ratio( 'positive', 2e-2 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The space charge properties can then be plot, using ```matplotlib.pyplot```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEMCAYAAADTfFGvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuYXXV97/H3Z+6TyY0kAwQIJFwVlIJGRKviI1XRUwFbOEK14ikt9rS0p496erCeRz1Un1P0PFVbsZZWKloLKEpNK4IXtF7KLcjNEC4hCAQC5DYhk7nume/5Y60dNjt7JnvP7LVmL/J5Pc882Xut317rOzsz+zO/329dFBGYmZnVq22uCzAzs2JxcJiZWUMcHGZm1hAHh5mZNcTBYWZmDXFwmJlZQxwcZmbWEAeHmZk1xMFhZmYN6ZjrArKwbNmyWLly5VyXYWZWKHfeeefWiOjfV7vcgkPSGcDngHbgHyPir6rWvwH4LHAicF5EXJcuPwn4O2AhMAF8MiKunW5fK1euZO3atc3/JszMXsQkPVZPu1yGqiS1A5cDbwOOB86XdHxVs8eB9wH/UrV8CHhvRJwAnAF8VtLibCs2M7Op5NXjOAXYEBEbASRdA5wF3F9uEBG/StdNVr4wIh6qePyUpGeBfmAg+7LNzKxaXpPjhwJPVDzflC5riKRTgC7gkSbVZWZmDSrMUVWSlgNfBf5bREzWWH+RpLWS1m7ZsiX/As3M9hN5BceTwIqK54ely+oiaSHwHeAjEXFrrTYRcUVErI6I1f39+zwowMzMZiiv4LgDOEbSKkldwHnAmnpemLa/HvhK+UgrMzObO7kER0SUgIuBm4D1wNcjYp2kSyWdCSDpVZI2AecCfy9pXfry/wq8AXifpLvTr5PyqNvMzPamF+OtY1evXh0+j8NaTURw9e1PcGR/H6ceuXSuyzHbi6Q7I2L1vtq9KM8cN2tF92zayV9cfx9tgoc/+Xba2zTXJZnNSGGOqjIruv98ZCsAkwGPbBmc42rMZs7BYZaT9Zt37Xl8/1PPzWElZrPj4DDLyZM7hjj58ORqOU8ODM9xNWYz5+Awy8nmnSMcuWw+i3o72bzTwWHF5eAwy0FpYpJnnhth+aIeli/q4emdI3NdktmMOTjMcvDsrlEmA5YvToLjqQEHhxWXg8MsB5vTHsYhi3pZvrjXQ1VWaA4OsxxsHRwFYNn8bvrnd7NjaJzSxF7X6jQrBAeHWQ52Do0DsHheJ4vndQLw3EhpLksymzEHh1kOBobHADigr4sD5nUly4bG5rIksxlzcJjlYGBonI420dfVzqK0xzEwPD7HVZnNjIPDLAcDw+MsnteJJBb3JsFRHr4yKxoHh1kOdg6NsygNjMXloaphD1VZMTk4zHIwMDy2JzDKPY4B9zisoBwcZjkYGBrfExgLHRxWcA4OsxwMVAxVtbeJhT0d7PTkuBWUg8MsB4OjJRb0PH/ftAU9nezyeRxWUA4Os4xFBLtHS/R1Px8cfd3t7B51cFgxOTjMMjZamqQ0GVXB0cHuMQeHFZODwyxj5Z7F/IrgmN/dwaB7HFZQDg6zjO0enQB4YY+jq8NDVVZYDg6zjA3u6XG071nW192xJ1DMisbBYZax8lxG3wuGqto9VGWF5eAwy1g5IPaaHB8tERFzVZbZjOUWHJLOkPSgpA2SLqmx/g2SfiGpJOmcqnUXSHo4/bogr5rNmqHW5HhfdwelyWC05Js5WfHkEhyS2oHLgbcBxwPnSzq+qtnjwPuAf6l67RLgY8CrgVOAj0k6IOuazZpld40eRzlEPEFuRZRXj+MUYENEbIyIMeAa4KzKBhHxq4i4F6j+E+ytwPcjYntE7AC+D5yRR9FmzTCYToLP73phjwPwBLkVUl7BcSjwRMXzTemyrF9rNueG0l7FvIqjqspHWHmC3IroRTM5LukiSWslrd2yZctcl2O2x/D4BJ3torP9+V+33rT3MTzuHocVT17B8SSwouL5Yemypr02Iq6IiNURsbq/v3/GhZo128j4JN0d7S9Y1tvZnq5zcFjx5BUcdwDHSFolqQs4D1hT52tvAt4i6YB0Uvwt6TKzQhgpTdDT+cJftXJwDI85OKx4cgmOiCgBF5N84K8Hvh4R6yRdKulMAEmvkrQJOBf4e0nr0tduB/6SJHzuAC5Nl5kVwsj4xN49jq7kV89DVVZEHftu0hwRcQNwQ9Wyj1Y8voNkGKrWa68Ersy0QLOMjJYm9+px9JR7HA4OK6AXzeS4WasaHZ/YExRlnuOwInNwmGVsZHxyr+CYVz6qynMcVkAODrOMJXMcL/xVKz/3UJUVkYPDLGPJUVUv7HG0tYmezjYHhxWSg8MsY8lQ1d6/ar2d7Yx4qMoKyMFhlrHR0gQ9VYfjQhIcQw4OKyAHh1nGRsYn6e7cOzh6uto9VGWF5OAwy1ityXFIh6ocHFZADg6zjI3WOBwXkuBwj8OKyMFhlqGJyWBsYorJ8a52n8dhheTgMMvQaCkJhlo9ju6OdobHfetYKx4Hh1mGRtNg6Kkxx9HT2bYnWMyKxMFhlqGRaXocPZ3te4LFrEgcHGYZGkmDobvGHEdPZ5uPqrJCcnCYZagcDLVOAOzp8OG4VkwODrMM7QmOWpPjnW2MlDxUZcXj4DDL0LRDVR3tTEwGpQmHhxWLg8MsQ9Mdjlte5l6HFY2DwyxDe3ocUxyOm7TxPIcVi4PDLEPTngDo28daQTk4zDI07eR4R7nH4aEqKxYHh1mGRqY9czwJE589bkXj4DDL0HQ9jj2T4+5xWME4OMwyNFqaZnI8XTbqOQ4rGAeHWYZGxifoaBMd7VMPVY14qMoKJrfgkHSGpAclbZB0SY313ZKuTdffJmllurxT0lWS7pO0XtKH86rZbLZGpriJEzx/UqCHqqxocgkOSe3A5cDbgOOB8yUdX9XsQmBHRBwNfAa4LF1+LtAdES8HXgm8vxwqZq1upDRR8yZO8Pz1qzw5bkWTV4/jFGBDRGyMiDHgGuCsqjZnAVelj68DTpckIIA+SR1ALzAGPJdP2Wazk9xvvHaPw5PjVlR5BcehwBMVzzely2q2iYgSsBNYShIiu4HNwOPA/4uI7VkXbNYMyf3Gp+hx+MxxK6giTI6fAkwAhwCrgA9KOrK6kaSLJK2VtHbLli1512hW02jJPQ578ckrOJ4EVlQ8PyxdVrNNOiy1CNgG/A5wY0SMR8SzwM+B1dU7iIgrImJ1RKzu7+/P4Fswa9zIND2Ornb3OKyY8gqOO4BjJK2S1AWcB6yparMGuCB9fA5wc0QEyfDUmwAk9QGnAg/kUrXZLI2MT0x5VFVbm+jqaNtzrodZUeQSHOmcxcXATcB64OsRsU7SpZLOTJt9CVgqaQPwAaB8yO7lwHxJ60gC6J8i4t486jabreSoqtrBAclJgO5xWNF05LWjiLgBuKFq2UcrHo+QHHpb/brBWsvNimC6oSpI5jl8OK4VTREmx80Ka7rDcSEJDk+OW9E4OMwyNFqavsfR7aEqKyAHh1mG6ulxeHLcisbBYZah0WmuVQXJSYDucVjRODjMMjI5GYxN7Hty3MFhRePgMMvIyDT3Gy/r7vDkuBWPg8MsI8NjSXD0ThccnW2+H4cVjoPDLCMj6aT3tENVHe2MusdhBePgMMtIucexr8lxnwBoRePgMMtIedJ7+uDwHIcVj4PDLCPl4JhujsOH41oROTjMMlLuSezrqKrSZFCacK/DisPBYZaR4Tp7HIDPHrdCcXCYZeT5OY7pTwCsbGtWBA4Os4wM1zM5nl7HasQ9DisQB4dZRuo5qqq707ePteJxcJhlZM9RVV3TT45XtjUrgoaDQ1KfpKl/E8wMgOGx9KiqjunmODw5bsWzz+CQ1CbpdyR9R9KzwAPAZkn3S/q0pKOzL9OseEZKE3S2i452T47bi0s9PY6bgaOADwMHR8SKiDgQeB1wK3CZpPdkWKNZIQ2PTUw7vwHPB4evV2VF0lFHm/uBH0XEvZULI2I78E3gm5I6syjOrMhGS/sOju4OT45b8dTT43gQ+LSkX0n6lKSTqxtExHjzSzMrtuGxiWlP/oOKoSpf6NAKZJ/BERGfi4jXAKcB24ArJT0g6WOSjs28QrOCGhmf/u5/UDE57qEqK5C6j6qKiMci4rKIOBk4HzgbWJ9ZZWYFNzxeR4/Dh+NaAdUdHJI6JL1D0teA75IMYf1WZpWZFdzI+ATddQ9VucdhxVHP4bhvlnQlsAn4A+A7wFERcV5EfLveHUk6Q9KDkjZIuqTG+m5J16brb5O0smLdiZJukbRO0n2Seurdr9lcGamjx+HJcSuieo6q+jDwL8AHI2LHTHaSnjB4OfBmkgC6Q9KaiLi/otmFwI6IOFrSecBlwLskdQD/DPxuRNwjaSngyXhrefXMcbS1ia72Nt/MyQpln8EREW9qwn5OATZExEYASdcAZ5Ec6lt2FvDx9PF1wOclCXgLcG9E3JPWs60J9Zhlrp45DkiuV+Xbx1qRNDLHIUnvkfTR9Pnhkk6p8+WHAk9UPN+ULqvZJiJKwE5gKXAsEJJukvQLSX8+RX0XSVorae2WLVvq/bbMMjMyvu/zOMC3j7XiaeRaVV8AXkNyRBXALpLhp6x1kJyl/u7033dKOr26UURcERGrI2J1f39/DmWZTa+eM8fBt4+14mkkOF4dEX8MjACk8x1ddb72SWBFxfPD0mU126TzGotIzhvZBPwkIrZGxBBwA/CKBuo2y11EMDQ+QV/3voOjr6uDobFSDlWZNUcjwTGeTnIHgKR+oN7+9R3AMZJWSeoCzgPWVLVZA1yQPj4HuDkiArgJeLmkeWmgnMYL50bMWs5oaZKJyWBe176PP5nX1c7QmHscVhz1HFVV9jfA9cCBkj5J8uH+v+t5YUSUJF1MEgLtwJURsU7SpcDaiFgDfAn4qqQNwHaScCEidkj6a5LwCeCGiPhOA3Wb5a4cBH3T3IujrK+7g92j7nFYcewzOCQpEl+TdCdwOiDg7IhYX9lmuu1ExA0kw0yVyz5a8XgEOHeK1/4zySG5ZoVQDoJ53fX1OLbsGs26JLOmqafH8SNJ3wS+HREPkNyPA0ldkt5EMrz0I+DLmVVpVjDl+4331TVU1eGhKiuUeoLjDOD3gKslrQIGgB6SIafvAZ+NiLuyK9GseJ7vcex7qCqZ4/BQlRVHPScAjpAcivuF9L4by4DhiBjIujizoir3IObVcThuMsfhHocVRyOT4+X7bmzOqBazF41yj6OvzjmO4fEJJiaD9jZlXZrZrDVyOK6Z1WlPj6Oeo6rSeZBhnwRoBeHgMMvAnsNx6+hx9Kbh4nkOK4pGr1V1rKTj0osPmtkUyiFQV48jnUAf8jyHFURdcxySjie5Ym35TPEOSb9VdVl0M0uVJ7vrO3M8abPbPQ4riHp7HP8E/EVEvCwiXgb8L+Br2ZVlVmxDYyV6Otvqmuwuz3H4XA4rinrOHP/vJJc8XyjpveXFwFJJfxQRX8iyQLMi2j1Wqqu3Ac+f6+HLjlhR1POT/UqgE3hV1fKudJ2ZVRkanahrfgOenwdxj8OKop4TAH9f0iPA5yJiA4CkI4HfjIgLsy7QrIiGxibqutwIeKjKiqfeEwA/APyHpOvT52cDF2dTklnx7R4r1XW5EajscXioyoqhruCIiG9Luo/k/t+QXJ9qQ3ZlmRVbQz2O9FwPX3bEiqLuS45ExEbgixnWYvaisXu0xNK++m6Q2d3RRpvc47Di8JnjZhnYPVaq66xxAEn0dXUw6KOqrCAcHGYZ2DVSYkFP/dcQXdDTwa4RB4cVg4PDrMkigsGGg6OTQQeHFYSDw6zJRsYnKU0G87s7637N/J4Odo2OZ1iVWfM4OMyabNdIEgAeqrIXKweHWZPtSie5GwmO+d0dHqqywnBwmDVZuefQ6BzHcw4OKwgHh1mTPT9UVf8cx8Kejj2vM2t1Dg6zJisPOc2v8zyOctvR0iRjpcl9NzabYw4Osyab2VBV0tYnAVoR5BYcks6Q9KCkDZIuqbG+W9K16frbJK2sWn+4pEFJH8qrZrOZeG4GQ1Xlth6usiLIJTgktQOXA28DjgfOT29HW+lCYEdEHA18Brisav1fA9/Nulaz2Sr3Ghoaqkp7HD4k14ogrx7HKcCGiNgYEWPANcBZVW3OAq5KH18HnC5JAJLOBh4F1uVUr9mM7Rop0dfVXtdtY8sWODisQPIKjkOBJyqeb0qX1WwTESVgJ8ntaeeT3OP8/0y3A0kXSVorae2WLVuaVrhZo3aNjDc0TAWw0ENVViBFmBz/OPCZiBicrlFEXBERqyNidX9/fz6VmdUwOFraM/RUr/KwlifHrQga++meuSeBFRXPD0uX1WqzSVIHsAjYBrwaOEfSp4DFwKSkkYj4fPZlmzWu0SvjgoeqrFjyCo47gGMkrSIJiPOA36lqswa4ALgFOAe4OSICeH25gaSPA4MODWtlzw2Ps2hefTdxKisPbe0c9lCVtb5chqrSOYuLgZuA9cDXI2KdpEslnZk2+xLJnMYGknuc73XIrlkRDAyPc8C8xuY4ujra6Otqd3BYIeTV4yAibgBuqFr20YrHI8C5+9jGxzMpzqyJBobGWdzbWHAALJ7XxcCQg8NaXxEmx80KY2IyeG5knMUNDlUBLJ7XycDQWAZVmTWXg8OsiZ4bHiciCYFGLZ7XyYCHqqwAHBxmTVT+4J9RcPR2ucdhheDgMGui8gf/4t7Gh6oWzev05LgVgoPDrInKPY5FM+pxdDIwNE5yFLpZ63JwmDVRucdxwAwmxw+Y10VpMnz2uLU8B4dZE5UPp53J4bjlXooPybVW5+Awa6Lyh/7CmZzH0euzx60YHBxmTbRzeJyFPR0NXVK9rHzuh3sc1uocHGZNtGNobEYn/8Hzh/Du8CG51uIcHGZNtH33GAf0zSw4yhPq23c7OKy1OTjMmmjr4Bj982cWHEv6umgTbBscbXJVZs3l4DBroq2Doyyb3z2j17a3iSV9XWwZdI/DWpuDw6xJJieD7bvHWDrDHgfA0r5utrrHYS3OwWHWJAPD40xMxox7HADLFnR5qMpanoPDrEnKH/hLZxEcSY/DQ1XW2hwcZk2yJQ2OZbMYqlo230NV1vocHGZNUu4p9M9yqGpobIKhMV+vylqXg8OsSZoxVLWsrzvdloerrHU5OMyaZOvgKO1tmtEFDsuWLUiGubZ4uMpamIPDrEm27BplaV8XbTO4TlVZ//weAJ59zsFhrcvBYdYkm3eOsHxx76y2cfCiJDieeW6kGSWZZcLBYdYkT+8cYfnCnlltY2lfF53tYvNOB4e1LgeHWZM8vXNkT49hptraxEELe3h653CTqjJrvtyCQ9IZkh6UtEHSJTXWd0u6Nl1/m6SV6fI3S7pT0n3pv2/Kq2azeu0aGWfXaInlswwOgOWLetzjsJaWS3BIagcuB94GHA+cL+n4qmYXAjsi4mjgM8Bl6fKtwDsi4uXABcBX86jZrBFPpx/0s53jAFi+qJenPcdhLSyvHscpwIaI2BgRY8A1wFlVbc4CrkofXwecLkkRcVdEPJUuXwf0Spr5gfJmGSj3EJrZ44iIWW/LLAt5BcehwBMVzzely2q2iYgSsBNYWtXmt4FfRISPVbSWUu5xHDzLyXFIjqwaK02yw7eQtRZVmMlxSSeQDF+9f4r1F0laK2ntli1b8i3O9ntP7RxGgoOaEBzlXstTA54gt9aUV3A8CayoeH5YuqxmG0kdwCJgW/r8MOB64L0R8UitHUTEFRGxOiJW9/f3N7l8s+k9vm2I5Qt76OqY/a/UiiXzkm1uH5r1tsyykFdw3AEcI2mVpC7gPGBNVZs1JJPfAOcAN0dESFoMfAe4JCJ+nlO9Zg15bPsQhy+d15RtHbG0L9nmNgeHtaZcgiOds7gYuAlYD3w9ItZJulTSmWmzLwFLJW0APgCUD9m9GDga+Kiku9OvA/Oo26xej20b4oglfU3Z1vzuDpb2dfH49t1N2Z5Zs3XktaOIuAG4oWrZRysejwDn1njdJ4BPZF6g2QztHi2xdXC0aT0OgMOXznOPw1pWYSbHzVpVeS7iiCYGxxFLHBzWuhwcZrP02LZkSKlZQ1WQzHM8tXOY0dJE07Zp1iwODrNZ2rg1DY5lzetxrFrWR0RytJZZq3FwmM3Sw88MsnxRDwt7Zn4Dp2rHHDQfgIeeGWzaNs2axcFhNksPPr2LYw9a0NRtHtU/nzbBg8/saup2zZrBwWE2C6WJSTZsGeS4g5sbHD2d7axc1sdDTzs4rPU4OMxm4bHtQ4yVJpve4wA49sAFPOQeh7UgB4fZLDywOflgPy6L4Dh4Ab/atpuhsVLTt202Gw4Os1m4d9MAXe1tHHvw/KZv+8RDFzEZsO6p55q+bbPZcHCYzcJdTwxw/CEL6e5ob/q2Tzp8MQB3Pz7Q9G2bzYaDw2yGShOT3LdpJyetWJzJ9pfN7+awA3q5+wkHh7UWB4fZDD30zCDD4xOZBQfAr61YzF2P78hs+2Yz4eAwm6FbNm4DYPXKAzLbx6tXLeGpnSM+g9xaioPDbIZ+vmErq5b1cdgBzbvUSLVfP3oZAD/bsDWzfZg1ysFhNgNjpUlu3biN16Uf7Fk5clkfyxf18LMNvh2ytQ4Hh9kM3PboNobGJnj9MdkGhyRef8wyfvrQVl8p11qGg8NsBm64bzN9Xe284djs72//tpcvZ9doiZ897OEqaw0ODrMGjU9McuMvn+b0lx5ET2fzz9+o9utHLWNRbyf/fu/mzPdlVg8Hh1mDblr3NDuGxjn75ENy2V9XRxv/5cTl3HDfZgaGxnLZp9l0HBxmDfrKLY+xYkkvpx17YG77/N1Tj2C0NMk31m7KbZ9mU3FwmDXg9ke3c/uj27ngNStpb1Nu+33p8oW8etUS/uGnGxke8yS5zS0Hh1mdJieDy258gAMXdPPuVx+R+/4/+JbjeHbXKFf+/NHc921WycFhVqev3f44dz62gw+99Th6u7KfFK92yqolvPWEg/jcDx/mYd+nw+aQg8OsDvduGuAT/34/rz9mGee+8rA5q+MTZ7+c+d0d/MnVd/HcyPic1WH7NweH2T7c88QA773ydpbN7+az7zoJKb+5jWr9C5IaHtkyyIVfvoOdQw4Py19uwSHpDEkPStog6ZIa67slXZuuv03Syop1H06XPyjprXnVbPu30sQkX/75o7zriltY0NPB1X9wKkvnd891Wbzh2H4+866TuPuJAc68/Gfcll5s0SwvHXnsRFI7cDnwZmATcIekNRFxf0WzC4EdEXG0pPOAy4B3SToeOA84ATgE+IGkYyPCh5ZYJp7dNcKNv3yaf/r5r3h0627eeFw/nzrnRA5c0DPXpe3xmycewvJFPfzp1Xfzritu5bRj+3nPqUfw+mOW5XJSou3fcgkO4BRgQ0RsBJB0DXAWUBkcZwEfTx9fB3xeyZjAWcA1ETEKPCppQ7q9W3Kq3V5EJieDofEJhkZL7B6bYPvuUZ4aGGHzzmE2PDvIvZt28uAzu4iAEw5ZyBff80reesJBczo8NZVXHrGEH3zgNL70s41cdctj/MFX1tLd0cbJhy/mpcsXclT/fA5a2MNBC7tZOr+beZ3t9Ha1093R1pLfjxVHXsFxKPBExfNNwKunahMRJUk7gaXp8lurXntoFkUODI1xzhdfmEcR8cLntV4Y0z6taztR40VR1apmm5oFTb3fpu57n9upse+s3qs6/p9KE8Hw+NQd1QPmdXLiYYt5+8uXc8bLDubYgxZM2bZV9Ha1c/GbjuGiNxzFLRu38eMHn+Wuxwe45vYnpvxe2wS9ne10drTRJtEm0d7GnsdtbdCePmaG+TLTWJppoDkGn/eS5Qv52/NPznQfeQVH5iRdBFwEcPjhh89oG+1t4rhaHxaa9ml5/3W0mX47tX5p9lpSY8OqWriv/dRuU2Pf1W1q/nbmtO+91u/7o6K6SbtEX3cHfd3tyb9dHSya18khi3pZvriHhT2d+9xmq+rqaOO0Y/s5Lb3o4uRk8MyuEZ59bpRnd42yffcoQ2MTDI9PMDw2wdDYBKWJSSYDJiKYnAwmI5iYJP03eT4TM3vVzF9Y/UfO/m7FAb2Z7yOv4HgSWFHx/LB0Wa02myR1AIuAbXW+loi4ArgCYPXq1TP6SVrQ08nl737FTF5q1lLa2sTyRb0sX5T9h4jtf/I6quoO4BhJqyR1kUx2r6lqswa4IH18DnBzJOMPa4Dz0qOuVgHHALfnVLeZmVXJpceRzllcDNwEtANXRsQ6SZcCayNiDfAl4Kvp5Pd2knAhbfd1kon0EvDHPqLKzGzuqNbkadGtXr061q5dO9dlmJkViqQ7I2L1vtr5zHEzM2uIg8PMzBri4DAzs4Y4OMzMrCEODjMza8iL8qgqSVuAx2axiWXA1iaV00yuqzGuqzGuqzEvxrqOiIj+fTV6UQbHbElaW88haXlzXY1xXY1xXY3Zn+vyUJWZmTXEwWFmZg1xcNR2xVwXMAXX1RjX1RjX1Zj9ti7PcZiZWUPc4zAzs4Y4OABJn5b0gKR7JV0vafEU7c6Q9KCkDZIuyaGucyWtkzQpacqjJCT9StJ9ku6WlPnVHRuoK+/3a4mk70t6OP33gCnaTaTv1d2Sqi/v38x6pv3+01sFXJuuv03SyqxqabCu90naUvEe/X4ONV0p6VlJv5xivST9TVrzvZJyuXFOHXW9UdLOivfqoznVtULSjyTdn/4u/o8abbJ7zyJiv/8C3gJ0pI8vAy6r0aYdeAQ4EugC7gGOz7iulwLHAT8GVk/T7lfAshzfr33WNUfv16eAS9LHl9T6f0zXDebwHu3z+wf+CPhi+vg84NoWqet9wOfz+nlK9/kG4BXAL6dY/3bguyQ3gzwVuK1F6noj8O95vlfpfpcDr0gfLwAeqvH/mNl75h4HEBHfi4hS+vRWkrsMVjsF2BARGyNiDLgGOCvjutZHxINZ7mMm6qwr9/cr3f5V6eOrgLMz3t906vn+K+u9DjhdM73pdnPryl1E/ITkPjxTOQv4SiRuBRZLWt4Cdc2JiNgcEb9IH+8C1gOHVjXL7D1zcOzt90hSutqhwBMVzzex93/UXAnge5LuTO+93grm4v06KCI2p4+fBg6aol2QrFZjAAAGYUlEQVSPpLWSbpWUVbjU8/3vaZP+4bITWJpRPY3UBfDb6fDGdZJW1Fift1b+/XuNpHskfVfSCXnvPB3iPBm4rWpVZu9ZXvccn3OSfgAcXGPVRyLi22mbj5DcZfBrrVRXHV4XEU9KOhD4vqQH0r+U5rquppuursonERGSpjpk8Ij0/ToSuFnSfRHxSLNrLbB/A66OiFFJ7yfpFb1pjmtqVb8g+XkalPR24F9Jbm+dC0nzgW8CfxYRz+W13/0mOCLiN6ZbL+l9wG8Cp0c6QFjlSaDyL6/D0mWZ1lXnNp5M/31W0vUkwxGzCo4m1JX7+yXpGUnLI2Jz2iV/doptlN+vjZJ+TPLXWrODo57vv9xmk6QOYBGwrcl1NFxXRFTW8I8kc0dzLZOfp9mq/LCOiBskfUHSsojI/BpWkjpJQuNrEfGtGk0ye888VEVylAnw58CZETE0RbM7gGMkrZLURTKZmdkROfWS1CdpQfkxyUR/zSNAcjYX79ca4IL08QXAXj0jSQdI6k4fLwN+neR+9s1Wz/dfWe85wM1T/NGSa11V4+Bnkoyfz7U1wHvTI4VOBXZWDEvOGUkHl+elJJ1C8pmadfiT7vNLwPqI+OspmmX3nuV9NEArfgEbSMYC706/yke6HALcUNHu7SRHLzxCMmSTdV3vJBmXHAWeAW6qrovk6Jh70q91rVLXHL1fS4EfAg8DPwCWpMtXA/+YPn4tcF/6ft0HXJhhPXt9/8ClJH+gAPQA30h//m4Hjsz6Paqzrv+b/izdA/wIeEkONV0NbAbG05+tC4E/BP4wXS/g8rTm+5jmKMOc67q44r26FXhtTnW9jmRu896Kz6235/We+cxxMzNriIeqzMysIQ4OMzNriIPDzMwa4uAwM7OGODjMzFrcvi62OIPtfSq9OOL69EKIDV3qxsFhZtb6vgyc0YwNSXotyflLJwIvA14FnNbINhwcZmYtLmpcbFHSUZJuTK9R91NJL6l3cyTnEHUB3UAnyflYdXNwmM0xSR+X9KH08X/OcBuLJf1RcyuzFncF8CcR8UrgQ8AX6nlRRNxCcmLn5vTrpoho6OoA+821qsxaQTqWrIiYrLU+Il47w00vJrm/R10fHlZs6cUNXwt8o2J6onwpnd8iuRJAtScj4q2Sjia5p0759hHfl/T6iPhpvft3j8P2W+kd1N6cPv6EpL+t0ea96eXF75H01XTZByT9Mv36s4q2Uy1fqeSOe18huY7YCkkfkfSQpJ+R3BSr3Haw4jXrJf1DOon5PUm96bp/TYcn1un5y+j/FXCUkrvQfTpt9x5Jt6fL/l5Se5PfQps7bcBARJxU8fVSgIj4VkS8rMbXW9PXvhO4NSIGI2KQ5DYSr2lo73lcV8Vf/mrFL5K7u/0YeDfwHaC9av0JJNd0WpY+XwK8kuS6P33AfJLrFJ081fL0dSuBSeDU9Hm57TxgIcm1qj6UrhuseE0JOCl9/nXgPeU60n97SYJoadr+lxW1v5Tk8uid6fMvAO+d6/fcX7P6ea3+P/5P4Nz0sYBfq3M77yK5llsHyfzGD4F3NFKLh6psvxURP0mHjj4AvDEiJqqavAn4RqSXyI6I7ZJ+F7g+InYDSPoW8HqSX9xay+9Kt/VYJHdhI11+faRXYtbU9z1/NCLuTh/fSfLBAfCnkt6ZPl5Bcv+Hp6teezpJQN2RDmX0MsVl5q31Sbqa5Da1yyRtAj5G8gfP30n63yQBcA3JxRb35TqSn+37SCbKb4yIf2ukHgeH7bckvZzk3s3bIrn9ZpZ2z+A1oxWPJ4BeSW8EfgN4TUQMpfcT6anxWgFXRcSHZ7BfazERcf4Uqxo+RDf9A+n9s6nHcxy2X0rvOfE1kvsyD6b3ZKl2M3CupKXpa5YAPwXOljQvvf/JO9NlUy2v5Sdp2970XirvaKD0RcCONDReApyaLt8FLKho90PgHCV3hUTSEklHNLAfsym5x2H7HUnzgG8BH4yI9ZL+ErgMuLGyXUSsk/RJ4D8kTQB3RcT7JH2Z5P4ZkNzn4650uzWXV4uIX0i6lmRY4VmSmyvV60bgDyWtBx4kuQcEEbFN0s/TM4u/GxH/Mx3C+J6kNpL7Sfwx8FgD+zKryffjMDOzhnioyszMGuLgMDOzhjg4zMysIQ4OMzNriIPDzMwa4uAwM7OGODjMzKwhDg4zM2vI/wfIM6AWunrs+QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1163321d0>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa8AAAEMCAYAAABz4zlBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcXFWd9/HPN510ks6+EQJJCEuQTRQMCOjMMIAS0RFwxMEVHJRxxBl9XHGc58FRmZHxGZ1xRlRGGMEHRUSROKKIgivDErawBEgIBBISsnV2Olv/nj/uKVJ0qtJ1O1XddSvf9+tVr7517rn3/rpS6V+dc0+do4jAzMysSAYNdABmZmZ5OXmZmVnhOHmZmVnhOHmZmVnhOHmZmVnhOHmZmVnhOHmZmVnhOHmZmVnhOHmZmVnhDB7oAFrVxIkTY8aMGQMdhplZodx7772rImJSb/WcvBpkxowZzJ07d6DDMDMrFEmLa6nnbkMzMyscJy8zMyscJy8zMyscJy8zMyscJy8zMyscJy8zMyscJy8zMyscJy+zFrdo5UZ+v2BVxX1Pr9rE7xas7OeIzPacv6Rs1uJO+ZffAPD0F9+4y76T/++vq+4za2ZueZntJbq7Y6BDMKsbJy+zvcT6rm0DHYJZ3Th5mbW49rbsv/maTVsHOBKz+nHyMmtx40YMAaBzs5OXtQ4nL7MWN66jHYDOTS/tNozwPTArrqZLXpKGSbpb0oOSHpH0D6n8QEl3SVoo6fuS2lP50PR8Ydo/o+xcn07lj0s6vax8dipbKOnisvLc1zBrdmM7spbXmh4trw1btg9EOGZ10XTJC9gCnBIRrwBeCcyWdAJwGfCViDgE6AQuSPUvADpT+VdSPSQdAZwLHAnMBi6X1CapDfga8AbgCODtqS55r2FWBDtbXi9NXms2uhvRiqvpkldkNqanQ9IjgFOAG1L51cBZafvM9Jy0/1RJSuXXRcSWiHgKWAgcnx4LI2JRRGwFrgPOTMfkvYZZ0+toz77O2bn5pd2GPVtiZkXSdMkLILWQHgBWALcCTwJrI6LUz7EE2D9t7w88C5D2rwMmlJf3OKZa+YQ+XMOs6ZU+ZrnlZa2kKZNXROyIiFcCU8laSocNcEg1kXShpLmS5q5c6Sl3rDmUxmX0bGl56LwVWVMmr5KIWAvcDpwIjJVUms5qKrA0bS8FpgGk/WOA1eXlPY6pVr66D9foGe8VETErImZNmjSpj7+1WWOs7ZG8VqfkNXRwU/8ZMKuo6d61kiZJGpu2hwOvA+aTJbG3pmrnATel7TnpOWn/bZGNAZ4DnJtGCh4IzATuBu4BZqaRhe1kgzrmpGPyXsOs6QXZW3Xlhi0vKS8972hv6/eYzPZUM07MOwW4Oo0KHARcHxH/LelR4DpJXwDuB65M9a8EviNpIbCGLBkREY9Iuh54FNgOXBQROwAkfQi4BWgDroqIR9K5PpXnGmZFsnx9FxFBaazR8+u7Bjgis75ruuQVEfOAYyqULyK7/9WzvAs4p8q5LgUurVB+M3BzPa5hVhRd27pZ/8J2xqTvfS138rICa7puQzOrs7IO7vKEtXxdV8/dZoXh5GW2Fyklr+7uYMUGt7ysuJy8zFpcsPO7XsvXvQBkw+a37cjaXB56ZEXk5GW2F5gyehhtg8QzazYDsHh19nPa+OEDGZZZnzl5me0FBrcNYvr4Dp5elSWtp1dtAuDAiSMHMiyzPnPyMmtxpa8kHjhxBItS0npq1SbaBolp44Z7aRQrJCcvs72AlCWvp1ZtpLs7WLRqI9PHdzCkzX8CrJj8zjXbSxw+ZTRd27pZsGIjDz67jsP2HTXQIZn1mZOXWYsrdQoeN2McADc9sJSla1/guBnjX7LfrEicvMz2AgKmj+9g39HDuPzXTwJw/IHj8ap0VlROXmYtrjQeQxJvmzUVgCP3G82R+41OFQYoMLM90HRzG5pZ/ZUm4/2bU2dyxH6jmTVjPJIQbnpZMTl5me1FhrQNYvZRUwY6DLM95m5DsxbXW6+gew2tiJy8zPYC1ToHPWDDisrJy6zF9TaDhmfYsCJy8jLbG1RpYbnhZUXl5GVmZoXj5GXW4jxgw1qRk5fZXsADNqzVOHmZtbpemlYer2FF5ORlthdQlSZWtXKzZufkZWZmhePkZdbiopd+w972mzUjJy+zvUDVARv9GoVZ/Th5mbW43gZkeMCGFVHTJS9J0yTdLulRSY9I+nAqHy/pVkkL0s9xqVySvippoaR5ko4tO9d5qf4CSeeVlb9K0kPpmK8q3bXuyzXMiqDquAw3vaygmi55AduBj0XEEcAJwEWSjgAuBn4VETOBX6XnAG8AZqbHhcDXIUtEwCXAq4HjgUtKySjVeX/ZcbNTea5rmJnZwGi65BURyyLivrS9AZgP7A+cCVydql0NnJW2zwSuicydwFhJU4DTgVsjYk1EdAK3ArPTvtERcWdkM5Je0+Ncea5h1vR67TbsnzDM6qpPyUvSCElt9Q6mwnVmAMcAdwGTI2JZ2rUcmJy29weeLTtsSSrbXfmSCuX04RpmhVBtxWSvpGxFVVPykjRI0jsk/VTSCuAxYFm6L/UlSYfUOzBJI4EfAh+JiPXl+1KLqaEfGPtyDUkXSporae7KlSsbFJlZPr0OhXfTywqo1pbX7cDBwKeBfSNiWkTsA7wWuBO4TNK76hWUpCFkievaiPhRKn6+1FWXfq5I5UuBaWWHT01luyufWqG8L9d4iYi4IiJmRcSsSZMm1f4LmzVYtQEbnmDDiqrW5HVaRHw+IuZFRHepMN1P+mFE/Dnw/XoElEb+XQnMj4gvl+2aA5RGDJ4H3FRW/p40IvAEYF3q+rsFeL2kcWmgxuuBW9K+9ZJOSNd6T49z5bmGmZkNgMG1VIqIbaVtSZ+KiMt2V2cPvQZ4N/CQpAdS2d8BXwSul3QBsBh4W9p3M3AGsBDYDLw3xbNG0ueBe1K9z0XEmrT9QeDbwHDgZ+lB3muYFUHvAzbcb2jF02vyknR9+VPglcAuyateIuL3VP/2yakV6gdwUZVzXQVcVaF8LnBUhfLVea9hVmTuNbSiqqXltT4i3ld6IsnfcTIrkF4Xo3TDywqolntel/Z4/plGBGJmjVN9SZR+DsSsTnpNXhHxFICkien5mt0fYWZm1lh5vqS8y70jM2t+nmHDWlGe5OUOBrOCqj4vr/9bWzHlSV7+gGZWSL0sRukRG1ZAbnmZ7QU8w4a1mjzJ69MNi8LMzCyHmpNXRDzcyEDMrDE8YMNaUU3TQ5VImkX2Pa8D0rEim4Di6AbEZmZ1UrXbsH/DMKubXMkLuBb4BPAQ0N1LXTMzs4bIm7xWRsSchkRiZg3h6aGsFeVNXpdI+hbwK2BLqbBszS0za0JVv8/l4YZWUHmT13uBw4Ah7Ow2DMDJy6xJ+Xtc1oryJq/jIuJlDYnEzBrGAzas1eT5nhfAHZKOaEgkZmZmNcrb8joBeEDSU2T3vDxU3qzJ1dJpGBFVl00xa0Z5k9fshkRhZg1VdWJe5ysrqJqSl6RDgMkR8Yce5a8BljciMDOrj1rGa0Q4kVmx1HrP61+B9RXK16d9ZtbMqq2k7CEbVlC1Jq/JEfFQz8JUNqOuEZmZmfWi1uQ1djf7htcjEDNrjJoGbDQ8CrP6qjV5zZX0/p6Fkt4H3FvfkMys3jxgw1pNraMNPwLcKOmd7ExWs4B24OxGBGZm9VHLDBtZHWcyK46akldEPA+cJOlPgaNS8U8j4raGRWZmdeMZNqzV1DpUXpG5Hbh9d3XqGp2ZmVkFtd7zul3S30iaXl4oqV3SKZKuBs6rR0CSrpK0QtLDZWXjJd0qaUH6OS6VS9JXJS2UNE/SsWXHnJfqL5B0Xln5qyQ9lI75qtK0An25hlmr8KdOK5pak9dsYAfwPUnPSXo0TRG1AHg78K8R8e06xfRtdp3J42LgVxExk2w5lotT+RuAmelxIfB1yBIRcAnwauB4sqVcxqVjvg68v+y42X25hlmReMCGtZqakldEdEXE5RHxGuAA4FTgmIg4ICLeHxH31yugiPgtsKZH8ZnA1Wn7auCssvJrUpfmncBYSVOA04FbI2JNRHQCtwKz077REXFn6uK8pse58lzDrBBqnWHDrEjyzipPRGyLiGURsbYRAVUxOSKWpe3lwOS0vT/wbFm9Jalsd+VLKpT35RpmhVFt0l1PxmtFlTt5DbTUYmro58S+XkPShZLmSpq7cuXKBkRmZmZQnOT1fKmrLv1ckcqXAtPK6k1NZbsrn1qhvC/X2EVEXBERsyJi1qRJk3L9gmaNEjV8DquljlkzyZW80ojDcb3XrLs57BzNeB5wU1n5e9KIwBOAdanr7xbg9ZLGpXhfD9yS9q2XdEIaZfieHufKcw2zwnDnoLWavOt5TQbukXQfcBVZQqjrRzZJ3wNOBiZKWkI2avCLwPWSLgAWA29L1W8GzgAWApuB9wJExBpJnwfuSfU+FxGlQSAfJBvROBz4WXqQ9xpmReEBG9aKciWviPh7Sf+brCXzXuA/JF0PXBkRT9YjoIh4e5Vdp1aoG8BFVc5zFVmC7Vk+l52zhJSXr857DbOiqDrDhptkVlB9GW0YZKPxlgPbgXHADZL+uc6xmZmZVZSr5SXpw2T3iVYB3wI+ERHbJA0i+8LyJ+sfopntCXcJWivKe89rPPCWiFhcXhgR3ZLeVL+wzKyeqq2Y7JWUrajydhsO65m4JF0GEBHz6xaVmdVNTUPl3TqzgsmbvF5XoewN9QjEzBrIAzasxdS6JMpfkw0xP0jSvLJdo4A/NCIwMzOzamq95/Vdsu9D/RM7Z1sH2FD2/Skza0I1fc/LM2xYwdS6kvI6YB3Z8idmVjBVl0Tp1yjM6qeme16Sfp9+bpC0Pj02lJ43NkQzM7OXqrXl9dr0c1RjwzGzegt6b2F5tKEVTd6Jec+RNCpt/72kH0k6pjGhmVm9eHooazV5h8r/74jYIOm1wGnAlcA36h+WmdVNTQM2zIolb/LakX6+EbgiIn4KtNc3JDOrN8+wYa0mb/JaKumbwLnAzZKG9uEcZmZmeyRv4nkbaaHHiFhLNqP8J+oelZnVTW3TQ7nj0Iol78S8O4BhwDmSyo/9Rf1CMrN684ANazV5k9dNwFrgPmBL/cMxs3qrbYYNs2LJm7ymRsTshkRiZg3jFpa1mrz3vO6Q9PKGRGJmZlajvC2v1wLvlbSIrNtQQETE0XWPzMzqopYuQY/XsKLJm7y8dpdZAVX9npf7E62g8nYbPgP8EXBeWlE5gMl1j8rM6qamYfBueVnB5E1elwMnsnNplA3A1+oakZnVXdWh8v0bhlnd5O02fHVEHCvpfoCI6JTk6aHMzKxf5W15bZPURupkkDQJ6K57VGZWNzUN2HC/oRVM3uT1VeBGYB9JlwK/B/6x7lGZWb/weA0rqlzdhhFxraR7gVPJusvPioj5DYmsCUmaDfwb0AZ8KyK+OMAhmfWqpvEabnhZweS950VEPAY81oBYmlrqLv0a8DpgCXCPpDkR8ejARmbWu2pD4t3wsqKqKXlJ+uju9kfEl+sTTlM7HlgYEYsAJF0HnAnUNXlt39HNlu273kbs+ben5/d2aun+6e0clev03F/hmN7O4b4p60VEEJHdn+uOoLv0PLL7cW4ZFsvgNjF0cFtjr1FjvVHp58uA44A56fmfAXfXO6gmtT/wbNnzJcCr632R3y9cxfn/dU+9T9v0ekuaWR3ttk6lHLlLgq7pOrs/R+Xr9DzH7mMd3CY62gczYmgbHe2DGTl0MONGtLPfmGHsN3Y4B04cwSumjWXM8CEVIsyntgEbe65r2w4eX76BJ57fwKJVm1i1YQtrNm1l9aatbNqynS3bu+natuPFn9u7g4ig24mp5bz7hAP4/FlHNfQaNSWviPgHAEm/BY6NiA3p+WeBnzYsuoKRdCFwIcD06dP7dI6DJ43k78447CVlPT919vy/XulTac/RY7Xd99j9MZVOsWud3q+7S1HP6zboOr2do9JBu55j12P68u+zbUc3m7ZuZ/OWHWzaup1NW7bzbOdmbnm4i607dra8D9t3FG98+RTOPnZ/po7r2PVENarW9t3TVvELW3fwkwef4yfznuOup9awNfUaDGkTE0YMZcLIdsaPaGe/scMYOriNoYMHMWxI9rNtkGgblH00kIQEgyQGaedzUSrfozCtnx0xZUzDr5H3ntdkYGvZ863sPTNsLAWmlT2fmspeFBFXAFcAzJo1q0+fJ6eN7+DCPz64rzFawXV3B6s3beWJ5zfwwLNr+fXjK/iXW5/gq7ct4NzjpvPx17+MMR05W2M1fHLJuxjlju7gu3ct5su3PkHn5m0cNHEE7z7hAI6bMZ5DJ49k+vgOBrd5kXVrnLzJ6xrgbkk3pudnAVfXN6SmdQ8wU9KBZEnrXOAdAxuStZpBg8SkUUOZNGoorzlkIhf96SEs6dzMN37zJN+7+xlue2wFX3/XsRw9dWyu89ZzMcp1m7fxt9fdz2+eWMmJB03gI6fN5PgDx/vepvWrXB+NIuJS4L1AZ3q8NyL2iu95RcR24EPALcB84PqIeGRgo7K9wdRxHXzhrJfzw78+CYB3/udd3Lu4c0BiWd+1jXddeRd3PLmKL5x1FN99/6t59UETnLis3+VKXsreoUcAYyLi34DVko5vSGRNKCJujohDI+LglMjN+s0rpo3lBx84kfEj23n/NXNZuvaFmo6r14CN7u7gomvvY/6y9Xzz3a/iXScc4KRlA8YT85oVyH5jh3PV+cexdXs3H7v+gZrvVVUdsJHj2v91x9P8bsEq/uHMIznlsL3lVrc1q7zJ69URcRHQBdnEvIAn5jXrR9mI1MO5c9EafnTf0l7r12OGjeXruvjSLY9x6mH78I7j+zaS1qyePDGvWQGde9w0jp46hq/88okXh6fvTtXuvRq7/f71l0/Q3Q2fffOR7iq0ptDXiXkne2Jes4EzaJD4X6cdypLOF/jxA723vvbEivVd/PC+JZx7/DSmje/7d83M6mlPJuaFvWxiXrNmcvLLJnHo5JFce+di3jZrWtV6tSx3srs6/++uZ9jeHfzlaw7sU5xmjZB3tOFHgTOAoenxBkkXSHplI4Izs+ok8Y7jp/PgknU8vHTd7uvmLC/p7g5+MPdZ/uTQScyYOKJPcZo1Qt5uw1nAB8jm+dsf+CtgNvCfkj5Z59jMrBdnHzOVwYPET+Y9V7VOTQMSq9S595lOlq3r4uxj9u9bgGYNkjd5TSWb2/BjEfEx4FXAPsAfA+fXOTYz68WYjiGcdMhEbnl4+W6Hzfd1vMZPHnyOYUMGcdrhHhpvzSVv8toH2FL2fBswOSJe6FFuZv1k9pH78vTqzTz+/Ia6njciuO2xFbz2kEmMGJp76T+zhsqbvK4F7pJ0iaRLgD8A35U0gjqva2VmtTn18H0A+N0Tqyrur+l7XhXKnl69mSWdL/Anh07cg+jMGqPm5JWmhvo22ZIfa9PjAxHxuYjYFBHvbEyIZrY7k0cP46BJI7jjycrJK1NtJeXq/Ya/fWIlAH986KQ9Cc+sIWruC4iIkHRzRLwcmNvAmMwspxMPmsCP71/Kth3dDKnTUiR3LlrN1HHDOWCCRxla88n7Lr9P0nENicTM+uzEgyewaesOHnlu/S77ahps2HNBzQjue6aTVx0wrj4BmtVZ7rkNgf+R9KSkeZIekjSvEYGZWe1eOS1b32vekrUV9+cdbfjcui6eX7+FY6blWzfMrL/kHUJ0ekOiMLM9sv/Y4UwY0c68Jbt+WbmWmed7zrBx/zPZemHHuuVlTSrv9FCLJY0DZgLDynYtrmtUZpaLJI6eOoaHKiQvyD/Dxrwl62gfPIjDp4yuS3xm9ZYreUl6H/Bhsi8rPwCcAPwPcEr9QzOzPI6eOpbfPLGAzVu309G+Z9/Lmr9sPTP3GVm3wR9m9Zb3nflh4DhgcUT8KXAM2ZB5MxtgR+w3mu6ABc9vzH1sz57Fx5dv4LB93eqy5pU3eXVFRBeApKER8RjwsvqHZWZ5HTp5FAALVuyavPIM2FizaSsrNmzhsH1H1TM8s7rK27ewRNJY4MfArZI68f0us6Ywbdxw2gcPYkGPaaLyzrDx2PJsuP1hU5y8rHnlHbBxdtr8rKTbgTHAz+selZnlNrhtEAdNHFG55ZVjho1St+PLJjt5WfPq813diPhNPQMxsz03c/IoHni2c4/O8dSqTYwcOphJo4bWKSqz+ss72nAo8OfAjPJjI+Jz9Q3LzPpi5j4j+cmDz9G1bQfDhrQBNa6kXNa3+PTqTRwwoQP1tl6K2QDKO2DjJuBMYDuwqexhZk3ggAkdACzp3PyS8qp5qEL506s2McPzGVqTy9ttODUiZjckEjPbY9PGZ8nrmTWbOWSf7J5VTQM2Up1tO7pZ0vkCbzx6SqNCNKuLvC2vOyS9vCGRmNkem15KXqtra3n1LF7a+QLbu8MzyVvTqyl5lU3A+1qymeUfb8TEvJLOkfSIpG5Js3rs+7Skhenap5eVz05lCyVdXFZ+oKS7Uvn3JbWn8qHp+cK0f0Zfr2HWbCaMaKejvY1n1rzQp+OfXp3dBThwopOXNbdauw3f1NAodnoYeAvwzfJCSUcA5wJHAvsBv5R0aNr9NeB1wBLgHklzIuJR4DLgKxFxnaRvABcAX08/OyPiEEnnpnp/0cdrmDUVSUwf38Eza3a2vGpZEqVkSWeW9KaOG17nyMzqq6aWV0QsjojFwOeAdWXP1wOX1CuYiJgfEY9X2HUmcF1EbImIp4CFwPHpsTAiFkXEVuA64My06vMpwA3p+KuBs8rOdXXavgE4NdXPdY16/c5m9TZtfAfPrunRbVjte149+hOXrXuBwYPEPqOGVaxv1izy3vM6OiJenMswIjrJ5jdstP2BZ8ueL0ll1conAGsjYnuP8pecK+1fl+rnvYZZU5o6bjjPdm5+cfh7TUuipCrPre1i8uhhtA3yMHlrbnlHGw6SNC4lLSSNz3sOSb8E9q2w6zMRcVPOeJqKpAuBCwGmT58+wNHY3mrKmGFs3rqDDVu2M3rYkKywxgEbz619gf3GutVlzS9v8voXspWUf5CenwNcmucEEXFazmsCLAWmlT2fmsqoUr4aGCtpcGpdldcvnWuJpMFkU1yt7sM1dhERVwBXAMyaNSvPrQazupk8Oks+z6/r2pm8arRsXdeLqzKbNbNc3YYRcQ3ZgIrn0+MtEfGdRgTWwxzg3DRS8ECyxTDvBu4BZqaRhe1kAy7mRNZPcjvw1nT8eWRfsC6d67y0/VbgtlQ/1zUa/Pua9dmUMdlgi+Xru4DaBmwEQXd3sGzdC+w31oM1rPnlntswjbJryEg7SWcD/w5MAn4q6YGIOD0iHpF0fbruduCiiNiRjvkQcAvQBlwVEY+k030KuE7SF4D7gStT+ZXAdyQtBNaQJSP6eA2zprNvanktW9f1YlnVCTbKdqzatIVtO8LdhlYIe7bcap1FxI3AjVX2XUqFLsqIuBm4uUL5IrKRgj3Lu8i6O/f4GmbNaJ/R2YS6z5eSV40zbCxbm9UvtdzMmpnX+DZrMcOGtDGuY8iL3Yaw65D4neU7t59P9UstN7Nm5uRl1oImjx7G8rJuw1qs2rgVgImj2hsRklldOXmZtaB9xwzLOWADVm7YAsCEEV7Hy5qfk5dZC5o8atiLyQh2M2CjbM+qjVsY2zGE9sH+s2DNz+9SsxY0YWQ7azZtpbs7apxhI1i5YQuTRrrVZcXg5GXWgiaMHMr27mB91zZgN0uilJWv3LiFiU5eVhBOXmYtaOLIbNBFaRBGLVZt3MKkUU5eVgxOXmYtqDToYs2mrbkGbLjlZUXh5GXWgiakltfqjdmgjd7miN+8ZQebt+5wy8sKw8nLrAWVkteqTVupYbwGq1KSK3U3mjU7Jy+zFjS+o0fLq+oMG1n56k3ZvbEJTl5WEE5eZi1ocNsgxnYMYXWNAzbWbs7qje1w8rJicPIya1ETRrSzetMWooYhG50peY1z8rKCcPIya1ETRg59cah89Rk2Mms2Zd8HG9eRb/FKs4Hi5GXWosZ1DGHd5m011V27eSuDRO6Vl80GipOXWYsaO7ydtS/UNtqwc/NWxgwfwqBBvQ2qN2sOTl5mLWpMxxDWllpevUwP1blpm+93WaE4eZm1qDHDh7Blezdd27p7rdu5eStjfb/LCsTJy6xFlZLRuhe2vmTpk3Kl8s7NW93yskJx8jJrUWOHZ8lo247eb3pt2xGMG+HkZcXh5GXWovJ2A3qYvBWJk5dZixozfGcyqmU9L8+uYUXi5GXWosqTVyPqmw0kJy+zFlXebdjbDBsAo4YNbmg8ZvXk5GXWokYOHUxbji8de3YNKxInL7MWJSlXV6BbXlYkTZW8JH1J0mOS5km6UdLYsn2flrRQ0uOSTi8rn53KFkq6uKz8QEl3pfLvS2pP5UPT84Vp/4y+XsOs2Y1NyauWARuj3PKyAmmq5AXcChwVEUcDTwCfBpB0BHAucCQwG7hcUpukNuBrwBuAI4C3p7oAlwFfiYhDgE7gglR+AdCZyr+S6vX1GmZNbZRbXtaimip5RcQvImJ7enonMDVtnwlcFxFbIuIpYCFwfHosjIhFEbEVuA44U9nysKcAN6TjrwbOKjvX1Wn7BuDUVD/XNRrx+5vV26ihWUKqNsNG+ZANJy8rkqZKXj38JfCztL0/8GzZviWprFr5BGBtWSIslb/kXGn/ulQ/7zXMml6ehDSi3cnLiqPf362SfgnsW2HXZyLiplTnM8B24Nr+jG1PSboQuBBg+vTpAxyNWTbisFZeDsWKpN+TV0Sctrv9ks4H3gScGvHiSkRLgWll1aamMqqUrwbGShqcWlfl9UvnWiJpMDAm1c97jUq/2xXAFQCzZs2qYRUls8YqDcKoZcCGWZE0VbehpNnAJ4E3R8Tmsl1zgHPTSMEDgZnA3cA9wMw0srCdbMDFnJT0bgfemo4/D7ip7Fznpe23Arel+rmu0Yjf36zeRqZuw+5aVqQ0K5Bm6+T+D2AocGs2hoI7I+IDEfGIpOuBR8m6Ey+KiB0Akj4E3AK0AVdFxCPpXJ8CrpP0BeB+4MpUfiXwHUkLgTVkyYg+XsOsqY1OyWvT1h0V95caXh2zalrpAAAJt0lEQVTtbf0UkVl9NFXySsPXq+27FLi0QvnNwM0VyheRjRTsWd4FnFOPa5g1u9I9r41d23dbzyMNrWiaqtvQzOqrdM9r45bekpe/oGzF4uRl1sJK97yqtbxS93yuUYlmzcDJy6yFlboDe295OXlZsTh5mbWw0oCNDV3bKu7v2pYN5HDysqJx8jJrYSOHZveyNlRpeZVaZKOG+p6XFYuTl1kLK7Woqn3Nq3QvzC0vKxonL7MW1tv3t0rdiSOdvKxgnLzMWph6mf9pW3fWJBvX0d4f4ZjVjT9umbW4fzz75Rw2ZVTFfR88+WC6u4Nzj59Wcb9Zs1J4zrOGmDVrVsydO3egwzAzKxRJ90bErN7qudvQzMwKx8nLzMwKx8nLzMwKx8nLzMwKx8nLzMwKx8nLzMwKx8nLzMwKx8nLzMwKx19SbhBJK4HFfTx8IrCqjuHUU7PG5rjycVz5OK589iSuAyJiUm+VnLyakKS5tXzDfCA0a2yOKx/HlY/jyqc/4nK3oZmZFY6Tl5mZFY6TV3O6YqAD2I1mjc1x5eO48nFc+TQ8Lt/zMjOzwnHLy8zMCsfJqwlI+pKkxyTNk3SjpLFV6s2W9LikhZIu7oe4zpH0iKRuSVVHDkl6WtJDkh6Q1C+LmOWIrb9fs/GSbpW0IP0cV6XejvR6PSBpToNi2e3vLmmopO+n/XdJmtGIOPoQ1/mSVpa9Pu/rp7iukrRC0sNV9kvSV1Pc8yQd2yRxnSxpXdnr9X/6Ka5pkm6X9Gj6v/jhCnUa95pFhB8D/ABeDwxO25cBl1Wo0wY8CRwEtAMPAkc0OK7DgZcBvwZm7abe08DEfn7Neo1tgF6zfwYuTtsXV/q3TPs2NjiOXn934IPAN9L2ucD3++HfrZa4zgf+oz/fT+m6fwwcCzxcZf8ZwM8AAScAdzVJXCcD/z0Ar9cU4Ni0PQp4osK/ZcNeM7e8mkBE/CIitqendwJTK1Q7HlgYEYsiYitwHXBmg+OaHxGPN/IafVVjbP3+mqXzX522rwbOavD1qqnldy+P9QbgVElqgrgGRET8FlizmypnAtdE5k5grKQpTRDXgIiIZRFxX9reAMwH9u9RrWGvmZNX8/lLsk8qPe0PPFv2fAm7vlEGSgC/kHSvpAsHOpgyA/GaTY6IZWl7OTC5Sr1hkuZKulNSIxJcLb/7i3XSh6d1wIQGxJI3LoA/T91MN0ia1uCYatXM/wdPlPSgpJ9JOrK/L566nI8B7uqxq2Gv2eB6nMR6J+mXwL4Vdn0mIm5KdT4DbAeubaa4avDaiFgqaR/gVkmPpU+LzRBb3e0urvInERGSqg3nPSC9ZgcBt0l6KCKerHesBfUT4HsRsUXSX5G1Dk8Z4Jia2X1k76eNks4AfgzM7K+LSxoJ/BD4SESs76/rOnn1k4g4bXf7JZ0PvAk4NVJncQ9LgfJPoFNTWUPjqvEcS9PPFZJuJOsa2uPkVYfY+v01k/S8pCkRsSx1j6yoco7Sa7ZI0q/JPrXWM3nV8ruX6iyRNBgYA6yuYwx9iisiymP4Ftl9xGbQkPfTnipPGBFxs6TLJU2MiIbPeShpCFniujYiflShSsNeM3cbNgFJs4FPAm+OiM1Vqt0DzJR0oKR2shvsDRmlloekEZJGlbbJBp9UHBU1AAbiNZsDnJe2zwN2aSFKGidpaNqeCLwGeLTOcdTyu5fH+lbgtiofnPo1rh73RN5Mdi+lGcwB3pNG0J0ArCvrIh4wkvYt3auUdDzZ3/VGfwghXfNKYH5EfLlKtca9Zv09QsWPiqN2FpL1Cz+QHqURYPsBN5fVO4NsRM+TZF1njY7rbLI+6i3A88AtPeMiGzX2YHo80h9x1RrbAL1mE4BfAQuAXwLjU/ks4Ftp+yTgofSaPQRc0KBYdvndgc+RfUgCGAb8IL3/7gYO6qd/u97i+qf0XnoQuB04rJ/i+h6wDNiW3lsXAB8APpD2C/haivshdjMCt5/j+lDZ63UncFI/xfVasvvd88r+dp3RX6+ZZ9gwM7PCcbehmZkVjpOXmZkVjpOXmZkVjpOXmZkVjpOXmZn1qrcJgvtwvn9OE/rOT5P35pqazMnLzMxq8W1gdj1OJOkksu83Hg0cBRwH/Emeczh5mZlZr6LCBMGSDpb08zSv6e8kHVbr6ci+Z9gODAWGkH1fs2ZOXmaGpM9K+njavqOP5xgr6YP1jcya3BXA30TEq4CPA5fXclBE/A/ZF9CXpcctEZFrJhXPbWi2l0n3FhQR3ZX2R8RJfTz1WLI1wmr6A2bFlibkPQn4QdntqtK0Z28hmzWlp6URcbqkQ8jW5Cst/3SrpD+KiN/Ven23vMwGUFqJ9nVp+wuS/r1Cnfek5UEelPSdVPZRSQ+nx0fK6lYrn6Fs9eJryOaenCbpM5KekPR7soU9S3U3lh0zX9J/phvrv5A0PO37ceoqekQ7l8H5InCwstV8v5TqvUvS3ansm5La6vwS2sAZBKyNiFeWPQ4HiIgfRcRRFR6np2PPBu6MiI0RsZFsGagTc129P+bA8sMPPyo/yFbJ/TXwTuCnQFuP/UeSzQM4MT0fD7yKbJ64EcBIsnntjqlWno6bAXQDJ6TnpbodwGiy+Q0/nvZtLDtmO/DK9Px64F2lONLP4WTJcEKq/3BZ7IeTLW8yJD2/HHjPQL/mfuzR+7Xnv/EdwDlpW8ArajzPX5DN/TmY7H7Xr4A/yxOLuw3NBlBE/DZ1430UODkidvSocgrwg0jLW0TEGknvBm6MiE0Akn4E/BHZH49K5fency2ObDVbUvmNkVYxkFRttv2nIuKBtH0v2R8vgL+VdHbanka2ftTyHseeSpYk70ndSsOpskSMNT9J3wNOBiZKWgJcQvah6+uS/p4sCV1HNkFwb24ge28/RDZ44+cR8ZM88Th5mQ0gSS8HpgCrI1tKvZE29eGYLWXbO4Dhkk4GTgNOjIjNaT2yYRWOFXB1RHy6D9e1JhMRb6+yK/fw+fQh7a/2JB7f8zIbIGndqmuBM4GNaV23nm4DzpE0IR0zHvgdcJakjrSG2tmprFp5Jb9NdYen9dj+LEfoY4DOlLgOA05I5RuAUWX1fgW8VdkK20gaL+mAHNcxq8otL7MBIKkD+BHwsYiYL+nzwGXAz8vrRcQjki4FfiNpB3B/RJwv6dtka3BBtk7Y/em8Fct7ioj7JH2frItnBdkikbX6OfABSfOBx8nWkCIiVkv6Q5qB4WcR8YnUnfQLSYPI1qO6CFic41pmFXk9LzMzKxx3G5qZWeE4eZmZWeE4eZmZWeE4eZmZWeE4eZmZWeE4eZmZWeE4eZmZWeE4eZmZWeH8f4qXLHLwMScWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x105ee00f0>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEMCAYAAAAvaXplAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcVOWd7/HPr3qloVm6my00uxghRhHQGDRKxCg4icuNjjqJ0YmJGuOYGcdk8JoYr4OvJCZ3cq+5GseMGXXCuMRJDEbiEjWJcckAigt0QCQijQhIA90svf/uH+c0FF1V3VVdS3cV3/frVS+qznmq6leH7v7W8zxnMXdHREQkWqS/CxARkYFH4SAiIjEUDiIiEkPhICIiMRQOIiISQ+EgIiIxFA4iIhJD4SAiIjEUDiIiEqO4vwvoq5qaGp80aVJ/lyEikldWrlz5gbuP7K1d3obDpEmTWLFiRX+XISKSV8xsYzLtNKwkIiIxFA4iIhJD4SAiIjHyds5BRA5fbW1t1NfX09zc3N+lDFjl5eXU1tZSUlLSp+crHEQk79TX11NZWcmkSZMws/4uZ8Bxd3bs2EF9fT2TJ0/u02tkfVjJzH5qZtvM7M0E683Mbjez9Wb2upnNynZNIpLfmpubqa6uVjAkYGZUV1en1bPKxZzDvcCCHtYvBKaFtyuAH+egJhHJcwqGnqW7fbI+rOTufzCzST00OQe434Prlb5sZsPNbKy7b8l2bSIZ8cYjsH1t/HUWgZkXw4hJOS1JJF0DYc5hHLAp6nF9uCwmHMzsCoLeBRMmTMhJcSK9evQr0NEKxPum5tDZDvO/leuqRNKSV7uyuvvd7j7H3eeMHNnr0d8iudHZAZ+4Hm7eFXsrKgPv6O8KJQduvvlmfvCDH/TY5vbbb2f69Ol87nOfS+m1d+3axZ133plOeSkbCOGwGRgf9bg2XCaSJ7yX1b2sl8PGnXfeydNPP82SJUtSet7hGg5LgS+Eey2dCOzWfIOI5INbb72VI488kpNPPpm1aw/OO/3sZz/jhBNOYObMmVx55ZV0dHRw1VVXsWHDBhYuXMgPf/jDhO0A7r//fo455hiOPfZYLrnkEhYtWsTbb7/NzJkz+frXv56Tz5b1OQczewCYB9SYWT3wbaAEwN3vApYBZwHrgX3A32a7JpGMS7RniPaoKVgrV67kwQcfZNWqVbS3tzNr1ixmz55NXV0dDz30EC+88AIlJSVcffXVLFmyhLvuuosnnniC5557jpqamoTtZs+ezeLFi3nxxRepqamhoaGBxsZG3nzzTVatWpWzz5eLvZUu7mW9A1/Ndh0iWdPrsJGGlbLpfz22mjXvNWb0NWd8aCjf/sxHemzz/PPPc95551FRUQHA2WefDcAzzzzDypUrOf744wHYv38/o0aNinl+ona7d+/mggsuoKamBoCqqioaGw/9fBs2bODWW29l9+7dPPLII+l92AQGwt5KIiIFw9259NJL+c53vtOndj/60Y96fY8pU6Zwzz33cP7556dVa08UDiJpc+LvxkqwXBPSWdXbN/xsOeWUU7jsssu44YYbaG9v57HHHuPKK69k/vz5nHPOOfzDP/wDo0aNoqGhgaamJiZOnHjI8xO1O+200zjvvPO47rrrqK6upqGhgcrKSpqamnL6+RQOIiJ9MGvWLC688EKOPfZYRo0adWB4aMaMGSxevJgzzjiDzs5OSkpKuOOOO2LCIVG7E088kRtvvJFTTz2VoqIijjvuOO69915OOukkjj76aBYuXMj3v//9rH8+8zz9VjNnzhzXleBkQLh5GJz6T/DJ/xm77taxcPyX4Ix/zn1dBayuro7p06f3dxn9ZseOHdx44408/fTTfOlLX+KGG26I2y7edjKzle4+p7f3UM9BJOvy8wuYDFzV1dXcddddWX2PgXCcg4iIDDAKB5F0HBiW1YS0FBaFg4iIxFA4iGSCjpCWAqNwEEmHhoykQCkcREQkhsJBJC2akJbCpHAQEZEYCgeRTNCEtBQYhYNIOpIaMtKw0uFAlwkVEZE+0WVCRQ47mpA+HOkyoSKSgP7wH650mVAR6Z3mnfvPbxbB+29k9jXHfBQWfrfHJv15mdBHH32Uxx9/nMbGRi6//HLOOOOM9D5vHAoHkXT0NmRkoN7F4SUXlwk999xzOffcc9m5cyfXX3+9wkFEJEYv3/CzZSBcJnTx4sV89atfzcrnUziIZEQPE9JSkPrzMqG33XYbixYtYuHChcyaNSsrn0+XCRVJR3sLLB4Fp30LTrk+dv13J8IxF8JZt+W+tgJ2uF8m9Pbbb+e+++7j+OOPZ+bMmVx11VVx2+kyoSL9TUdCSw5de+21XHvttVl9Dx3nIJKOXiekDU1ISz5SOIiISAyFg0hGaEJaCovCQSQtSQwZ5elOH3J4UziIZIImpHMuX/e0zJV0t4/CQSQdmpDuF+Xl5ezYsUMBkYC7s2PHDsrLy/v8GtqVVUTyTm1tLfX19Wzfvr2/SxmwysvLqa2t7fPzFQ4iGaEJ6VwqKSlh8uTJ/V1GQdOwkkhaNCEthSkn4WBmC8xsrZmtN7NFcdZPMLPnzOxVM3vdzM7KRV0iGaMJaSkwWQ8HMysC7gAWAjOAi81sRrdm3wQedvfjgIuA3F4PT6SvNCEtBSoXPYcTgPXuvsHdW4EHgXO6tXFgaHh/GPBeDuoSEZEEcjEhPQ7YFPW4HvhYtzY3A0+Z2d8Bg4HTc1CXSAZpQloKy0CZkL4YuNfda4GzgP8ws5jazOwKM1thZiu0C5sMDJqQlsKUi3DYDIyPelwbLot2OfAwgLu/BJQDNd1fyN3vdvc57j5n5MiRWSpXpA80IS0FJhfhsByYZmaTzayUYMJ5abc27wLzAcxsOkE4qGsgA58mpKVAZT0c3L0duAZ4Eqgj2CtptZndYmZnh83+Efiymb0GPABc5jouXkSk3+TkCGl3XwYs67bspqj7a4CTclGLSHZoWEkKy0CZkBbJU711cE0T0pKXFA4imaAJaSkwCgeRdGhCWgqUwkFERGIoHEQyQsNKUlgUDiJp0YS0FCaFg4iIxFA4iKSjq1eQaG8lTUhLnlI4iIhIDIWDSEZoQloKi8JBJKtMo0qSlxQOIiISQ+EgkgmakJYCo3AQSYeOYZACpXAQyQhNSEthUTiIpEVHSEthUjiIiEgMhYNIJiSckAZNSEs+UjiIpENDRlKgFA4iIhJD4SCSFk1IS2FSOIiISAyFg0gmJJqQBjQhLflI4SCSjt6GjHoMDZGBS+EgkhEKASksCgeRtGhCWgqTwkFERGIoHEQyQXMLUmAUDiLpSGpCWsNKkn8UDiIZoZ6DFBaFg0haNCEthUnhICIiMRQOIpmgCWkpMAoHkXRoQloKVE7CwcwWmNlaM1tvZosStPlrM1tjZqvN7D9zUZdI5qjnIIWlONtvYGZFwB3Ap4B6YLmZLXX3NVFtpgE3ACe5+04zG5XtukQyQxPSUphy0XM4AVjv7hvcvRV4EDinW5svA3e4+04Ad9+Wg7pERCSBXITDOGBT1OP6cFm0I4EjzewFM3vZzBbEeyEzu8LMVpjZiu3bt2epXJE+0IS0FJiBMiFdDEwD5gEXAz8xs+HdG7n73e4+x93njBw5MsclisTRbcho+TsN/PDpdexrbQ8WaEJa8lTW5xyAzcD4qMe14bJo9cCf3L0N+IuZrSMIi+U5qE8kA4zmtg6+eO9ymprbaW7v4IaF0/u7KJE+y0XPYTkwzcwmm1kpcBGwtFubRwl6DZhZDcEw04Yc1CaSpoO9gt+v205TczulxRF+/dqWcKkmpCU/ZT0c3L0duAZ4EqgDHnb31WZ2i5mdHTZ7EthhZmuA54Cvu/uObNcmkkm/X7edyrJivnHmh9m8az/1O/f1d0kifZaLYSXcfRmwrNuym6LuO3BdeBPJP2asea+Rj4wbytypNQCseGcntf1clkhfDZQJaZH8FA4Zdbrz5/cbmTF2GNNGD6E4Yqzb2qQJaclbCgeRDPhgTyvNbZ1MH1tJSVGECVUV/OWDvf1dlkifKRxEMqBhbysAk2sGAzBl5GA2bA/DQRPSkocUDiJpCf7wN+xtA6B2RAUAU0YO4S879uI655LkKYWDSAY07GultCjCqMoyAGpHDKK1vZOOTvUaJD8lvbeSmd0Stl8FrHL3dVmrSiRfeFfPoZVxIwYRiQQ9hdFDywFo73SKNSEteSjpcHD3m8xsNDATOM/MjnD3L2evNJH80bCvnXHDBx14PHZYEA5tHU55fxUlkoZew8HM/i/w9x7YSnDA2pNZr0wkjzQ1tzGqtuzA4zHDunoOnZqQlryUzJxDE7DUzAYDmNmZZvZCdssSyRfBH/6m5nZGDj0YDjWDyyiOGG2d/VWXSHp67Tm4+zfN7G+A35lZK7AHiHs1N5HDVVuHM3LIwXCIRIxRlWW0dygdJD8lM6w0n+BiPHuBscAX3X1ttgsTyQtRQ0YjK8sOWVU9pIz2xlwXJJIZyQwr3Qh8y93nAecDD5nZaVmtSiTPOBYnHEq1K6vkrWSGlU6Luv+GmS0E/guYm83CRPJN9LASQNXgUto7XRPSkpdSPgjO3bcA87NQi0heqxpcesjj6jAcXMc5SB7q0xHS7r4/04WI5DMHhg0qOWRZ9ZAy3NHQkuQlnT5DJB3hkNGg0iKKiw79daoaXIoDbe3aY0nyj8JBJAMqSotiltUMCYaZWjo6cl2OSNoUDiIZUFFaErNseEUQDu3qOUgeUjiIpCUYVorXcxhRUYpjtOlAOMlDCgeRDBhUGrtX+IiKoDfR2qEJack/CgeRdIQT0oPLYnsOQ8tLcKBdcw6ShxQOImno9K5hpdieQyRiFEU0rCT5SeEgkoa9rUGvIF44ABRHIgoHyUsKB5E0NO1vBaCiLH44RCKmM7NKXlI4iKShqbkNgME99Bw0IS35SOEgkoaucEjUcyiKmCakJS8pHETS0Li/q+cQu7cSQFEkQpt6DpKHFA4iadjdFQ5lsUdIAxQVGd7ZSXObeg+SXxQOImk4MKzUQ88BYOe+1pzVJJIJCgeRNHQNKxUVxf9VKo4YADv3tuWsJpFMUDiIpKGxuec/+kURw4Bd6jlInlE4iKShaX/P4dB1jYed+9RzkPyicBBJw8Geg8VdH/QcXHMOkndyEg5mtsDM1prZejNb1EO7z5qZm9mcXNQlkq7G/e09ri8OJ6Q1rCT5JuvhYGZFwB3AQmAGcLGZzYjTrhL4GvCnbNckkikHeg4Wv+cQMSMSMQ0rSd7JRc/hBGC9u29w91bgQeCcOO3+Gfge0JyDmkTS1tLeQUsSxy+URkzDSpJ3chEO44BNUY/rw2UHmNksYLy7P56DekQyYlcyvQEziossubYiA0i/T0ibWQT4F+Afk2h7hZmtMLMV27dvz35xIj1o2BvdG4g/rARQUhRRz0HyTi7CYTMwPupxbbisSyVwNPA7M3sHOBFYGm9S2t3vdvc57j5n5MiRWSxZpHcf7GnB6O28SUaJeg6Sh3IRDsuBaWY22cxKgYuApV0r3X23u9e4+yR3nwS8DJzt7ityUJtIn21rbDn4IMGENEBJkeYcJP9kPRzcvR24BngSqAMedvfVZnaLmZ2d7fcXyZZtTS29NwJKIsbu/W10dOrsrJI/4p+EPsPcfRmwrNuymxK0nZeLmkTSta2pmYqSXr5fmVFSBO7ByfdqhpTlpjiRNPX7hLRIvtre1MKIwaXho8TDSqXFwRlb39+tvbQlfygcRPpoW1MLVYPjX8fhIKOsKAiOrY0KB8kfCgeRPtre1EJVRdhz6GFCurQ4+DV7X+EgeUThINJH2xqbqTowrJRYSZERMdiqYSXJIwoHkT7Y29LO3taO3oeVLAiGmiFlbFE4SB5ROIj0QdcQUVUSE9IAY4eVa1hJ8orCQaQPNjXsA2BMZRK7prozemi5JqQlrygcRPpg0879AIwaWh4sSDghHSwfM6xcw0qSVxQOIn1Qv3MfpUWRg3sr9aJ2xCCamtvZrXMsSZ5QOIj0QX3DfsaNGETEejklRtijmFg9GICNDXuzXZpIRigcRPqgfuc+akcMilrS84T0xOoKADbu2JfFqkQyR+Eg0gebdu6ndkQF9HrKbsCdCVVd4aCeg+QHhYNIipqa22jY23poz6GXCemK0mJGVZap5yB5Q+EgkqL12/YAMG3UkCSfEfQuJlUPVjhI3lA4iKTora1hOIyu7H1UKapHMbG6gg0faFhJ8oPCQSRFb21rorQ4cmAeIdDzhDTAh8dU8sGeFj7Yk9xFgkT6k8JBJEVvbdvD1JFDKIoYyU5IA0wfOxSAte83ZbE6kcxQOIik6K2tezhydLf5hh5O2d3lqDGVANRtacxGWSIZpXAQScGOPS1s3rX/QC8gOUHPoXpIGSMry6jbop6DDHwKB5EUvFa/C4CZ44cHCzy5I6S7TB87lNXv7c5GaSIZpXAQScGqd3cRMfjouGHd1vQ+rAQwe8II1m5torFZ51iSgU3hIJKCVzft4sNjhjK4rDhckvyENMDxk0bgDis37sxOgSIZonAQSVJreyevbNzJ7InDY1cm7DgcumLmhOEUR4wV7zRkvD6RTFI4iCTplXd3sre1g09MG5niMw/2HCpKizl63DBeWL8js8WJZJjCQSRJv1+3neKIMXdq9cGFKU5IA5x21Cheq9/F9iYdDCcDl8JBJAnuzjN1W5k9cQSV5SVxWiQ3IQ1w+vTRuMNzf96WuQJFMkzhIJKENVsaWbd1D58+9kPd1qQ2IQ0wfWwltSMG8djr72WuQJEMUziIJOGXr2ympMj49EfHpvjM2B6FmfHZWbX8cf0HbGrQWVplYFI4iPSiqbmNh1Zs4owZYxgxOME1o3s8fUZs7+KCObUAPPDf72agQpHMUziI9GLJn96lqbmdK06ZEruytwnpBGpHVHDW0WO5/6WN7NrXmmaFIpmncBDpwfamFu54dj2nHjmSY8fHOb7hgAQ9hx56FF87fRp7W9v5P799K70iRbJA4SCSgLvzzUffoLm9g299ekaiVsm8UNzFR46u5JITJ3LfS+/w8gYd9yADi8JBJIEf/vYtnly9lX9acBRHJH1J0O563sX1GwuOYnL1YL7ys5W8o6vEyQCicBDppq2jk1sfX8Ptz7zF+bNrufzkyb0/KcUJ6S5Dyoq557LjceCCf32JVZt2pVyvSDbkJBzMbIGZrTWz9Wa2KM7668xsjZm9bmbPmNnEXNQlEs3d+cO67Zx35wv85Pm/cMmJE/neZ4/BevrD38cJ6WiTawbz8JUfp7Qowmd//CKLf72GbY3Nab+uSDqKe2+SHjMrAu4APgXUA8vNbKm7r4lq9iowx933mdlXgNuAC7Ndmxze3J3tTS28Xr+bFRt38sSbW3hnxz7GDR/EnZ+bxVkpHdOQ+oR0tCNHV/L4tSfznWV/5p4X/sJ9L73DKdNGMu/DI5k5fgTTRg+hvKQohXpE0pP1cABOANa7+wYAM3sQOAc4EA7u/lxU+5eBz2ermDf+8Cv2bY3aO8TB4/wCe/cH3dp0/77o8dp0a+Qc2iZYH+91o9oEbx5T36FtwGPaeMyy2Dbdljm4dV8f84Q426v75z605u6fKW7NMa/rMfX1tj0PvE636vZ5Cbs7ymjsKGVnexkNbSVsaCqmfncLLe2dABRHjBMmV3H1J4/g3JnjKC3OYKc6yd7F8IpSvnf+MVz9yan8x0sbeWrNVp4JT7ERMagZUsaYIUVMHtzK6OJ9DI00U1HcyaBIB4OsnRLaKPY2IuZEwucYXfediHX9T3hYU7yfiOzJ5XsVuiNOuZDq0bVZfY9chMM4YFPU43rgYz20vxz4TbwVZnYFcAXAhAkT+lRM6/J/52NNz/XeUApam5XQNHg0zYM/RGT00VR9eC6l006AwTUpvlJvf/hT/5M4sXow3zzrKG6cuZ+dda/Q/O5KIjs3MHR/PRU7m0CXgjjsrZs6pyDCIWlm9nlgDnBqvPXufjdwN8CcOXP6NNh7xN/+hA+ao09Z4LG/v+4Yh36Ljfmy7HG+l8d8i43X5tD1cd78kEXmcerr/o3P+/DeB9r08DnjfuONbWPxPkNP9catOfYzQfeOQWy9seXFvjftLdC6F1r3BP+2NFHS+B5Vu+th10Z46z+h7qdgRTDlVDjxajji9KSHhOIUGltDsvbvhJd/DK8uwRrrqQIYNgFGTYWquTBkNFRUwaARUFYJRaVQXEZnpBQvKqUjUoJjdDp0YnQC7kanGx0HerfhzaL/7Z2n8jlin5yW9Gd2Csuk6jFZf49chMNmYHzU49pw2SHM7HTgRuBUd8/auYyHVaV6Ln4peB1t8P7rUPdreP1hWHI+TDwZzr0DRkzq+bkZmJA+8Dqv3AdP3QQtu+GIT8H8m4Kwquz9D0HXIJhmJSRTchEOy4FpZjaZIBQuAv4muoGZHQf8K7DA3XUeY8mtohIYNzu4zbsBXr0ffnsL/PhkuGhJ8Ae6V2lMSHe0wa+ugdcfhMmnwJnfgTFHp/QRRDIt67uyuns7cA3wJFAHPOzuq83sFjM7O2z2fWAI8HMzW2VmS7Ndl0hcxaVw/JfgK3+EYbWw5ALY8Lv0XrOn3kVnJzzyxSAY5v1PuORXCgYZEHIy5+Duy4Bl3ZbdFHX/9FzUIZK04RPgb5fBv58FD18KX34WqqfGaZjmhPRzt0LdUjhjMcz9u75WK5JxOkJaJJGKKrj4P4P7v7wy+JafSCqT1102vgjP/wCO+zx8/Jq+1SiSJQoHkZ5UTYEF34X65bDyp7Hrk5qQjtOmvQUe+1qwJ9LC2/oWLiJZpHAQ6c2xFwUTxc8uhpamBI1SnJBeeR98sA7+6n9D6eCMlCmSSQoHkd6YwfxvB8cgLL8n9ed37120t8AffwgT5sK0T2WmRpEMUziIJKN2DkydDy/+CNqiT4rXh+McXnsAmt6DU7+h4SQZsBQOIsk66VrY9wHUPRa7LuEf+dgj11n+bzD6ozBlXoYLFMkchYNIsiadEhwx/cp9B5elOiH93qvw/hsw+1L1GmRAUziIJCsSgVlfgHeehx1vd1uZ5IT0K/dD8SA45q+zUqJIpigcRFJx7MXBv6t/kfxzujoOHe3BAW9H/RWUD8t4aSKZpHAQScXQD8H4j8GarjO8pDAh/e5LsG8HTP9MVkoTySSFg0iqpp8dnMW1YcPBZclMSNc9BsXlwenARQY4hYNIqmaE54us+3XyE9Lu8OdfB7vDlg3JankimaBwEEnV8Akwagas/23Uwl4mpLesgsbNMP3TWS9PJBMUDiJ9ccT8YA6hbV/vbd3h7WfD52lISfKDwkGkL6bOh45WeOePybXf8HsYfTQMGZXdukQyROEg0hcTPh4cr9A1tNTTAW3t++Hdl2FyMleUExkYFA4ifVFSDpNOhh3re2/bsAE6WnS6DMkrCgeRvpp6WtSDXiakI8UwcW7WSxLJFIWDSF9N/kTybWuP1y6sklcUDiJ9NeojybedMi9bVYhkhcJBpK8iESip6KVROKykyWjJMwoHkXR0HbewvyFxm9IhwcWCRPJIcX8XIJLXzr49OGI6Uc9gzhfhyDOhqCS3dYmkSeEgko5BI+DMWxOvn/rJ3NUikkEaVhIRkRgKBxERiaFwEBGRGAoHERGJoXAQEZEYCgcREYmhcBARkRgKBxERiWGe1AXSBx4z2w5s7OPTa4APMlhOpqiu1Kiu1Kiu1BRqXRPdfWRvjfI2HNJhZivcfcCd7EZ1pUZ1pUZ1peZwr0vDSiIiEkPhICIiMQ7XcLi7vwtIQHWlRnWlRnWl5rCu67CccxARkZ4drj0HERHpwWERDmb2fTP7s5m9bma/NLPhCdotMLO1ZrbezBbloK4LzGy1mXWaWcK9D8zsHTN7w8xWmdmKAVRXrrdXlZk9bWZvhf+OSNCuI9xWq8xsaRbr6fHzm1mZmT0Urv+TmU3KVi0p1nWZmW2P2kZfylFdPzWzbWb2ZoL1Zma3h3W/bmazBkhd88xsd9T2uikHNY03s+fMbE34u/i1OG2yu73cveBvwBlAcXj/e8D34rQpAt4GpgClwGvAjCzXNR34MPA7YE4P7d4BanK4vXqtq5+2123AovD+onj/j+G6PTnYRr1+fuBq4K7w/kXAQwOkrsuA/5ern6eo9z0FmAW8mWD9WcBvCC68fSLwpwFS1zzg1zneVmOBWeH9SmBdnP/HrG6vw6Ln4O5PuXt7+PBloDZOsxOA9e6+wd1bgQeBc7JcV527r83me/RFknXlfHuFr39feP8+4Nwsv19Pkvn80fU+Asw3MxsAdfULd/8D0MPFtjkHuN8DLwPDzWzsAKgr59x9i7u/Et5vAuqAcd2aZXV7HRbh0M0XCdK2u3HApqjH9cT+Z/QXB54ys5VmdkV/FxPqj+012t23hPffB0YnaFduZivM7GUzy1aAJPP5D7QJv5zsBqqzVE8qdQF8NhyKeMTMxme5pmQN5N/Bj5vZa2b2GzP7SC7fOByOPA74U7dVWd1eBXMNaTP7LTAmzqob3f1XYZsbgXZgyUCqKwknu/tmMxsFPG1mfw6/7fR3XRnXU13RD9zdzSzRrnYTw+01BXjWzN5w97czXWseewx4wN1bzOxKgt7Naf1c00D2CsHP1B4zOwt4FJiWizc2syHAfwF/7+6NuXjPLgUTDu5+ek/rzewy4NPAfA8H7LrZDER/g6oNl2W1riRfY3P47zYz+yXB0EFa4ZCBunK+vcxsq5mNdfctYfd5W4LX6NpeG8zsdwTfujIdDsl8/q429WZWDAwDdmS4jpTrcvfoGv6NYC5nIMjKz1S6ov8ou/syM7vTzGrcPavnXTKzEoJgWOLuv4jTJKvb67AYVjKzBcA3gLPdfV+CZsuBaWY22cxKCSYQs7anS7LMbLCZVXbdJ5hcj7tXRY71x/ZaClwa3r8UiOnhmNkIMysL79cAJwFrslBLMp8/ut7zgWcTfDHJaV3dxqXPJhjPHgiWAl8I98I5EdgdNYzYb8xsTNdckZmdQPB3M6shH77fPUCdu/9LgmbZ3V65nIHvrxuwnmBsblV469qD5EPAsqh2ZxHsFfA2wfBKtus6j2CcsAXYCjzZvS6CvU5eC2+rB0pd/bS9qoFngLeA3wJKs/FYAAAEVElEQVRV4fI5wL+F9+cCb4Tb6w3g8izWE/P5gVsIvoQAlAM/D3/+/huYku1tlGRd3wl/ll4DngOOylFdDwBbgLbw5+ty4CrgqnC9AXeEdb9BD3vw5biua6K218vA3BzUdDLBXOPrUX+3zsrl9tIR0iIiEuOwGFYSEZHUKBxERCSGwkFERGIoHEREJIbCQURkAOjtBIB9eL3bwpP21YUn6Evp1C0KBxGRgeFeYEEmXsjM5hIc43MMcDRwPHBqKq+hcBARGQA8zgkAzWyqmT0RnlfteTM7KtmXIzjOphQoA0oIjllKmsJBJMvM7GYzuz68/2IfX2O4mV2d2cokD9wN/J27zwauB+5M5knu/hLBAY5bwtuT7p7SkfAFc24lkYEgHNc1d++Mt97d5/bxpYcTXB8iqT8Okv/Ck+7NBX4eNV3QdWqY/0Fw1Ht3m939TDM7guC6LF2XJ3jazD7h7s8n+/7qOUjBCq+k9anw/mIz+1GcNl8IT139mpn9R7jsOjN7M7z9fVTbRMsnWXDltfsJzns13sxuNLN1ZvZHggsndbXdE/WcOjP7SThp+JSZDQrXPRoOI6y2g6do/y4w1YIrkX0/bPd5M/vvcNm/mllRhjeh9K8IsMvdZ0bdpgO4+y/c/eg4tzPD554HvOzue9x9D8FlCj6e0rvn4twluunWHzeCK3z9Dvgc8DhQ1G39RwjOQVQTPq4CZhOcp2YwMITgnDrHJVoePm8S0AmcGD7ualsBDCU4t9L14bo9Uc9pB2aGjx8GPt9VR/jvIIKwqQ7bvxlV+3SCU2+XhI/vBL7Q39tct7R/Zrv/P78IXBDeN+DYJF/nQoLzjxUTzDc8A3wmlVo0rCQFy93/EA7zXAfMc/eObk1OA37u4amX3b3BzC4BfunuewHM7BfAJwh+MeMtfzV8rY0eXI2LcPkvPTwDsCW+jvVf3H1VeH8lwR8GgGvN7Lzw/niCawe83+258wlCaHk45DCIBKcwl/xgZg8QXJK0xszqgW8TfLH5sZl9k+CP/IMEJwDszSMEP99vEExOP+Huj6VSj8JBCpaZfZTgWrw7PLjUYjbt7cNzWqLudwCDzGwecDrwcXffF16PojzOcw24z91v6MP7ygDk7hcnWJXy7q3hF6Er06lHcw5SkMJrFiwhuM7unvCaHt09C1xgZtXhc6qA54FzzawivH7GeeGyRMvj+UPYdlB4LY7PpFD6MGBnGAxHEVw4HqCJ4ELzXZ4Bzrfg6oCYWZWZTUzhfUR6pJ6DFBwzqwB+Afyju9eZ2T8D3wOeiG7n7qvN7Fbg92bWAbzq7peZ2b0E11+A4DoRr4avG3d5d+7+ipk9RND930ZwAZ5kPQFcZWZ1wFqC6wfg7jvM7IXw6NnfuPvXw6GGp8wsQnAtgq8CG1N4L5GEdD0HERGJoWElERGJoXAQEZEYCgcREYmhcBARkRgKBxERiaFwEBGRGAoHERGJoXAQEZEY/x/6STQUBXv8HAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x118efc550>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(grid.x, c_o.phi)\n",
    "plt.xlabel( '$x$ $\\mathrm{coordinate}$' )\n",
    "plt.ylabel('$\\Phi$ $\\mathrm{( eV )}$')\n",
    "plt.show()\n",
    "\n",
    "plt.plot(grid.x, c_o.rho)\n",
    "plt.xlabel( '$x$ $\\mathrm{coordinate}$'  )\n",
    "plt.ylabel(' $\\mathrm{charge density}$ $(\\mathrm{C m}^{-1})$')\n",
    "plt.show()\n",
    "\n",
    "plt.plot(grid.x, c_o.mf[site_labels[1]], label = '$\\mathrm{defect_1}$')\n",
    "plt.plot(grid.x, c_o.mf[site_labels[0]], label = '$\\mathrm{defect_2}$')\n",
    "plt.xlabel( '$x$ $\\mathrm{coordinate}$'  )\n",
    "plt.ylabel('$x_{i}$')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}