examples/notebooks/compare-ecker-data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparing with Experimental Data\n",
    "\n",
    "In this notebook we show how to compare results generated in PyBaMM with experimental data. We compare the results of the DFN model (see the [DFN notebook](./models/DFN.ipynb)) with the experimental data from Ecker et. al. [[3]](#References). Results are compared for a constant current discharge at 1C and at 5C."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we import pybamm and any other packages required by this example, and then change our working directory to the root of the pybamm folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Note: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "%pip install pybamm -q    # install PyBaMM if it is not installed\n",
    "import pybamm\n",
    "import os\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "os.chdir(pybamm.__path__[0]+'/..')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then load the Ecker data in from the `.csv` files using `pandas`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "voltage_data_1C = pd.read_csv(\"pybamm/input/discharge_data/Ecker2015/Ecker_1C.csv\", header=None).to_numpy()\n",
    "voltage_data_5C = pd.read_csv(\"pybamm/input/discharge_data/Ecker2015/Ecker_5C.csv\", header=None).to_numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the data is Time [s] vs Voltage [V]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We load the DFN model and select the parameter set from the Ecker paper [1]. We update the C-rate an `InputParameter` so that we can re-run the same model at different C-rates without the need to rebuild the model. This is done by passing the flag `[input]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# choose DFN\n",
    "model = pybamm.lithium_ion.DFN()\n",
    "\n",
    "# pick parameters, keeping C-rate as an input to be changed for each solve\n",
    "chemistry = pybamm.parameter_sets.Ecker2015\n",
    "parameter_values = pybamm.ParameterValues(chemistry=chemistry)\n",
    "parameter_values.update({\"Current function [A]\": \"[input]\"})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this comparison we choose a fine mesh of 1 finite volume per micron in the electrodes and separator and 1 finite volume per 0.1 micron in the particles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "var = pybamm.standard_spatial_vars\n",
    "var_pts = {\n",
    "    var.x_n: int(parameter_values.evaluate(model.param.L_n / 1e-6)),\n",
    "    var.x_s: int(parameter_values.evaluate(model.param.L_s / 1e-6)),\n",
    "    var.x_p: int(parameter_values.evaluate(model.param.L_p / 1e-6)),\n",
    "    var.r_n: int(parameter_values.evaluate(model.param.R_n_typ / 1e-7)),\n",
    "    var.r_p: int(parameter_values.evaluate(model.param.R_p_typ / 1e-7)),\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create a simulation using our model, parameters and number of grid points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim =  pybamm.Simulation(model, parameter_values=parameter_values, var_pts=var_pts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then solve the model for a 1C and 5C discharge "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "C_rates = [1, 5]  # C-rates to solve for\n",
    "capacity = parameter_values[\"Nominal cell capacity [A.h]\"]\n",
    "t_evals = [\n",
    "    np.linspace(0, 3800, 100), \n",
    "    np.linspace(0, 720, 100)\n",
    "] # times to return the solution at\n",
    "solutions = [None] * len(C_rates)  # empty list that will hold solutions\n",
    "\n",
    "# loop over C-rates\n",
    "for i, C_rate in enumerate(C_rates):\n",
    "    current = C_rate * capacity\n",
    "    sim.solve(t_eval=t_evals[i], solver=pybamm.CasadiSolver(mode=\"fast\"),inputs={\"Current function [A]\": current})\n",
    "    solutions[i] = sim.solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we plot the numerical solution against the experimental data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6AAAAEYCAYAAABCw5uAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABccElEQVR4nO3deXyU1dn/8c+Vyb5DFggJ+yogglIEEVxwV1BbtdhNayuP1lZra1t92h8itU93tdZutLVqa9214oJVEfcFAdn3nbAmAUISkpDl/P6YSQgh65DZku/79ZrXzNxzz8x1M2HOXPc55zrmnENEREREREQk0KJCHYCIiIiIiIh0DUpARUREREREJCiUgIqIiIiIiEhQKAEVERERERGRoFACKiIiIiIiIkGhBFRERERERESCQgmoiIiIiIiIBIUSUJFOwsy+bWaLzKzSzB5p9FiqmT1gZtvNrNTMNvnuZ4YoXBERkbBjZm+bWYWvrSw1s3UNHssxs7+b2W4zKzGztWZ2j5klhTJmkUijBFSk89gF3As83HCjmcUC84ERwEVAKjABKALGBTlGERGRcPdt51yy7zIUwMy6Ax8BCcAE51wKcD6QDgwMWaQiESg61AGISMdwzj0PYGZjgbwGD30N6AOc45wr9W3bB/w0uBGKiIhErO8BJcBXnHO1AM65HcBtIY1KJAKpB1Sk8zsPeK1B8ikiIiLN+7mZFZrZB2Z2tm/becDzdcmniPhPCahI55cB7A51ECIiIhHgR8AAIBeYA7xkZgNRWyrSYZSAinR+RUBOqIMQEREJd865T5xzJc65Sufco8AHwCWoLRXpMEpARTq/N4ELVaVPRESk3RxgeNvSK81Mv51FTpD+E4l0EmYWbWbxgAfwmFm8mUUD/wR2AM+Z2TAzizKzDDP7XzO7JKRBi4iIhAkzSzezC+vaTzP7MjAZeA24D28V+UfNrK9v/1wzu8/MRoUwbJGIowRUpPP4CVAO3Al8xXf7J865SrzFE9YCbwCHgIVAJvBJaEIVEREJOzF4lzMrAAqB7wBXOOfWO+f2A2cAVcAnZlaCd4mzYmBjiOIViUjmnAt1DCIiIiIiItIFqAdUREREREREgkIJqIiIiIiIiASFElAREREREREJCiWgIiIiIiIiEhTRoQ6gvTIzM12/fv1CHYaIiHRyixcvLnTOZYU6jo6mdlRERIKhuXY04hLQfv36sWjRolCHISIinZyZbQt1DIGgdlRERIKhuXZUQ3BFREREREQkKJSAioiIiIiISFAoARUREREREZGgiLg5oCIinVlVVRX5+flUVFSEOpQuIz4+nry8PGJiYkIdioiInCC1o8HX3nZUCaiISBjJz88nJSWFfv36YWahDqfTc85RVFREfn4+/fv3D3U4IiJygtSOBpc/7aiG4IqIhJGKigoyMjLUaAaJmZGRkaEz5SIinYTa0eDypx1VAioiEmbUaAaX/r1FRDoXfa8HV3v/vbtsAvrqit3kv/Mo3D8SZqV7r5c/HeqwREREwp5zjsc/2cb+siOhDkVERCJMl0xAa2odC1/8M93fugOKdwDOe/3SrUpCRaTL83g8jB49mhEjRnDKKafw29/+ltraWgDefvtt0tLSGD16NKNHj+a8884DYNasWSQmJrJv377610lOTg5J/BJ4hR/9i7NfPZf0X2dTe98ItZ0iIg2oHW1Zl0xAPVHGzMRnSbRGZ26ryql5457QBCUiEiYSEhJYunQpq1at4o033mDevHncc8/R78ZJkyaxdOlSli5dyptvvlm/PTMzk9/+9rehCDnimZnHzD4zs5ebeCzOzJ4ys41m9omZ9Wvw2F2+7evM7MKgBLv8abIW/IBcKyQKR9ShfJxO4IqI1FM72rIumYACRB3a2eR2O7STyx96n9+9uYGVO4txzgU5MhGR8JGdnc2cOXN46KGHWv0+vOGGG3jqqafYv39/kKLrVG4D1jTz2DeAA865QcD9wC8BzGw4MB0YAVwE/NHMPAGPdP5sqCo/ZpNVlePmzw74W4uIRBq1o8frusuwpOX5ht8eqzS+B54o44H567n/zfXkpMVz7rBsvpT4CcNX348V7/Q+d8pMGHVNCAIXka7inpdWsXrXoQ59zeG9Url76oh2PWfAgAHU1NTUDwt67733GD16NABXX301P/7xjwHvUKEbbriB3/3ud8ec6ZWWmVkecCnwM+B7TexyOTDLd/tZ4CHzVny4HHjSOVcJbDGzjcA44KOABlyc377tIiIhonY0PHXdBHTKTO+cz4ZncWMSSL30pzw/aiKFpZUsWLuPN9fspfKzJ+lvc7C6IbvFO3Bzb8VASaiIdDmTJk3i5ZePGykKwK233sro0aO54447ghxVRHsA+CGQ0szjucAOAOdctZkVAxm+7R832C/ft+04ZjYDmAHQp0+fE4u2mRO4O2szWL5iN5ecnHNiry8i0sl19Xa06yagdYnj/Nnes7aNejUzk+O4emxvrh7bm9r7riPq0LHzRa26nMIXf8x7NWdw3kk9SImPCfYRiEgn194zrIGyefNmPB4P2dnZrFnT3ChRr/T0dL70pS/xhz/8IUjRRTYzuwzY55xbbGZnB+p9nHNzgDkAY8eOPbG5JU2cwHUYyxJO567nljMqL428bokn9BYiIh1B7Wh46roJKHiTzTb0YDY3X7R7TQG3P7WM2OgoLhjegy+d3ocJA7TwrYh0HgUFBdx00018+9vfbvN32/e+9z0+97nPUV1dHeDoOoWJwDQzuwSIB1LN7F/Oua802Gcn0BvIN7NoIA0oarC9Tp5vW2CNuga2fwyLHga8uazhuLjmLd5xA7jtyRSemjGeaE+XLTMhIlJP7ejx1Dq0RVpek5stLY/nbp7AtZ/rzbvrC/jSXz9hym/f4a2nH/KWpdf6oiISgcrLy+vLx5933nlccMEF3H333W1+fmZmJldeeSWVlZUBjLJzcM7d5ZzLc871w1tQ6K1GySfAXOA63+2rfPs43/bpviq5/YHBwMKgBL7hdeqSzzpR1eXMSnyWxdsO8OD8DUEJQ0QkHKkdbZlFWpXXsWPHukWLFgX3TZc/3eR8UaY+WN+DWlFVw8vLd7Pj7Uf4n0O/O3aJl0b7iog0Z82aNZx00kmhDqPLaerf3cwWO+fGBisG3xDcO5xzl5nZbGCRc26umcUD/wTGAPuB6c65zb7n/Bi4AagGvuucm9fa+3RIOzorncYJqO8ouGPEOzy3JJ9/f3M8EwZmnNj7iIi0k9rR0GhPOxrwHlB/1zYLK6Ou8SaQab0B8143SijjYzxcdVoet0c92eT6ohX/vVtLuoiISLOcc2875y7z3Z7pnJvru13hnLvaOTfIOTeuLvn0PfYz59xA59zQtiSfHaaZkUGk5XHPtBH0y0ji1ic/Y9Wu4qCFJCIikSEYQ3DbvbZZWBp1Ddy+EmYd9F4315vZTBn62NLdnH//uzzywRYOVVQ1/z7Ln/YO29XwXRERCVdTZnpH9zQUkwBTZpIUF81fvnoa0VHGVX/6iNdW7glNjCIiEpYCmoA2WNvsb83scjnwqO/2s8AUi/QKPs2cFS5P7ElSXDSzXlrN2Hvf5H/+uYi5y3ZRVtlgcnHdUN/iHYDzXr90q5JQEREJL62MDBrSI4UXb5nIkJ4p3PSvxfxhwUaNAhIRESDwVXAfwL+1zQob7tSh65cFWjPriyZdPJsXR01kRX4xz3+Wz6srdvPfVXuJi45i0uBMJg/J4toPZxHT8HngfZ35szV/VEREwkvjSvJ1I3h8S5tlT5nJUzO+wI+eW86v/7uOz7Yf5AcXDmVoz+Z+EoiISFcQsAS0I9c269D1ywKtlfVFT85L4+S8NP7fpcNZtO0AryzfxYJ1Bby5Zh9fidsJTfX/NjOsV0REJCw0LtbnG8ETDzzwxasZ0SuVB+dv5MIH3uXikT25dcpgTspJDWnIIiISGoHsAT2Rtc0iWxvWF42KMsb17864/t25B9hWVMbhv+aQXLH7uH2LorN54b3NnNq3GyN6pRIX7Tl2h+VPN5vwioiIBNz82ceO/IH6ETw26hpmTB7INWN78/f3t/DIB1uZt3IPZw3J4ivj+3LusGw8UZE9+0ZERNouYAmoc+4u4C44prR8c2ubfcSxa5t1OX0zkuCS2ccN3620OH7PtTzyireOU6wniuG9UhndO53RvdM54/BbZL39A6zRWWdASaiI+MXj8XDyySfX358+fTp33nlnwN5v7ty5rF69OqDv8fbbbxMbG8sZZ5wRsPfo0pobqdNge3piLN+/YCjfPHMAj3y4lX8v3MaNjy2iV1o808f14eqxeeSkJTT9OiIiEUTtaMsCPQf0OA3XNgP+DvzTzDbiW9ss2PGElSaG78ZNmcmsUddw86EKlmw7wJLtB1i2o5inPt3BIx9u5f3YWVjU8Wedne+ss4h0cgEYAZGQkMDSpUs7Jr5WVFdXM23aNKZNmxbQ93n77bdJTk5WAhooaXm+AnpNbG+8KTGG284bzLfOGcj8NXv518fbue+N9Tzw5nrOHJzF1aflcf7wHsTHeI5/PRGRjqZ2tE06sh21SOtw7JAFtDuB6ppa1u8t5aQ5fbAmFgOvxbiu9385OTeNkblpnJybRl63BI4pMqyhuyJhp10LaDeedwfepTAarVPcXsnJyZSWlh6zrbi4mHHjxjF37lyGDh3Ktddey7nnnsuNN95IcnIyN954I6+//jo9e/bkySefJCsri02bNnHLLbdQUFBAYmIif/3rXxk2bBjXX3898fHxfPbZZ0ycOJFRo0axaNEiHnroIa6//noSEhL47LPP2LdvHw8//DCPPfYYH330EaeffjqPPPIIAK+//jp33303lZWVDBw4kH/84x8kJyfTr18/rrvuOl566SWqqqp45plniI+PZ/z48Xg8HrKysvj973/PpEmTWv13b24B7UgXkHb0BP8WtxaW8dySfJ5bnM+u4gpS46O57JRefOHUXE7t041IL5AvIsGjdjT829FgrAMqARDtG4przSz7UhyTTWHpEea8u5lvPb6ESb9awOjZb3DtnI/52Sur+XTuX6idqyVfRCJaC/PuTkR5eTmjR4+uvzz11FOkpaXVN2xPPvkkBw4c4MYbbwSgrKyMsWPHsmrVKs466yzuueceAGbMmMHvf/97Fi9ezG9+8xu+9a1v1b9Hfn4+H374Iffdd99x73/gwAE++ugj7r//fqZNm8btt9/OqlWrWLFiBUuXLqWwsJB7772XN998kyVLljB27NhjXiczM5MlS5Zw880385vf/IZ+/fpx0003cfvtt7N06dLjGk3pAC0ty9KG9a37ZSbx/QuG8t6PzuWf3xjHOcOyeX5JPl/400ec/Zu3eeDN9WwrKgv6YYlIJ6d2NCTtaNCH4EoHa2bZl25T72XeqElUVNWwfm8JK3ceYsXOYlbvKubRj7ZxXdSviGpi6O6R/86i9qQvaOiTSCRow7w7fzQ3dOj888/nmWee4ZZbbmHZsmX126OiovjiF78IwFe+8hU+//nPU1payocffsjVV19dv19lZWX97auvvhqPp+nvmalTp2JmnHzyyfTo0aN+Hs2IESPYunUr+fn5rF69mokTJwJw5MgRJkyYUP/8z3/+8wCcdtppPP/8837+K0i7NVWAr5nquPX7N+KJMiYNzmLS4CxKK6uZt2I3zy/Zye/mb+CBNzcwpk86V47J5bJRveieFBvgAxKRTk/tKBD8dlQJaKRrZdmX+BgPo/LSGZWXXv+Uqppaon/adLHh6NJdDJr5GgOzkhneK5XhOan11xnJcUd31PBdkdBrx7y7jlBbW8uaNWtITEzkwIED5OU1/T5mRm1tLenp6c3OgUlKSmr2feLivN81UVFR9bfr7ldXV+PxeDj//PN54oknWny+x+Ohurq6LYcmgdJS70IrbUZyXDRXj+3N1WN7s+tgOXOX7eKFJTuZ+eIqZr+0mkmDM7liTC7nD+9BYqx+zoiIH9SOtvj8QLWjGoLbGYy6Bm5fCbMOeq9badRjPFHNDt2tTMrh2+cMok/3RD7dsp+fz1vLV/++kNPufZNxP3uT6x5eyNx/PkD1i9/R8F2RUJsy0ztXpaGYBO/2ALj//vs56aST+Pe//83Xv/51qqqqAG+D+uyzzwLw73//mzPPPJPU1FT69+/PM888A4Bz7pizvSdi/PjxfPDBB2zcuBHwDl1av359i89JSUmhpKSkQ95f2qGDehd6pSdw01kD+e/tk5l32yS+Mak/a/eUcNuTSzntp29y6xOfMX/NXo5U13ZA0CLSZagdBYLfjuqUYVfVzNDdhIvu4XujhtZvOlB2hDV7DrF61yFW7z7Emt0lnLb/90RbxbGvV1VO6aszWZE0heE5qaQlxgTpQES6sFZGQPirbu5KnYsuuoivf/3r/O1vf2PhwoWkpKQwefJk7r33Xu655x6SkpJYuHAh9957L9nZ2Tz11FMAPP7449x8883ce++9VFVVMX36dE455ZQTig0gKyuLRx55hGuvvbZ+ONK9997LkCFDmn3O1KlTueqqq3jxxRebLJ4gARKA3oWTclI5KSeVH104jE+37uc/S3cxb+Vu5i7bRXpiDBePzGHaKb0Y17+71hcVkZapHQ1JO6oquF2Zn8No3az0pivvOmNA5eMA5KTFM6xnCsNyUr3XPVMZuOdVohf8VMN2RVrQrup9YaKpan+RRlVwAyRAFSYbO1Jdy3sbCpi7bBevr9pLeVUNPVLjuPTkXkwb3YtT8tJUSVeki1A7GhrtaUfVA9qVNVUwog2smTPaLjWXR6eOY83uQ6zdfYi1e0p4f2MhVTWOaVHv84uYvxFtR7w7F++g5sXvUHz4CN1O/7J+GIiIdEbt6V04gdoCsdFRTDmpB1NO6sHhI9XMX7OPuct28a+Pt/HwB1vo3T2BqaN6MfWUXgzrmaI2R0QkhJSASvs1M3zXc/7dnDUki7OGZNVvPlJdy6aCUvo+9n0Sy48c8zKemgoOv3o3Z83LYnCPZIb2TGFID+9lcI9kspLj9CNBJAJE+llbCbC2nOxsZ7XcliTGRjP1FG+yWVxexX9X7eHl5bv5y7ub+ePbmxiUnexLRnMYkJXsxwGJiHSsrtaOKgGV9mvHGe3Y6ChOykmF8j1NvlRuVBFXjM5l/d4S5q3cwxMLj/aspifGMCQ7hS/EfMilBX8lqWIPtSm5RJ1/N6ahu9KJOed08iWIIm0qSqd0AtVyW5KWEMM1Y3tzzdjeFJVW8urKPby8bBcPzF/P/W+uZ0SvVKae0ovLRuWQ1y3xBA9CRMKF2tHgam87qgRU/NPe4bvNDNu1tDx+esVIwPvHW1BSyYZ9pazfW8L6vaX02Poi0w79jgS8vaeeknzKn7+FP76xnt19pjIwK5mBWUkMzE6mT/dEYjwq7CyRLT4+nqKiIjIyMtR4BoFzjqKiIuLj40MdStcWoLX4GspIjuOr4/vy1fF92VNcwSsrvIWLfjFvLb+Yt5bT+nZj6qgcLh3Vi6yUuNZfUETCktrR4PKnHVURIgkOfwtR3D+yycS1wJPNJVF/oqDk6IK80VFGn4xEBmQmMyAribMr3+a0jQ8SW7Yb0vIwFT2SCFBVVUV+fj4VFRWt7ywdIj4+nry8PGJijq3erSJEQdTMdz1pvb3LiwXQ9qLDvLR8Fy8t28XaPSVEGUwclMm0U3px4ciepMarqrtIJFE7GnztbUeVgErw+FNgYlY6NFFxFwxmHaS4vIrNBaVsLihjY0EpWwrK2FxYysii17nXM4dEOzrvtII4nuh5B4X9p9E3I4l+GUn0y0gkK0VzTUXkeEpAgyhI1XJbs35vCXOX7uLFZTvZsb+c2Ogozh2azeWje3HOsGziYzxBi0VEJNIpAZXI5OdZcXf/SKyJ5+2xLCZWPkhN7dG/+4QYD326J9InI5G+3ROZXLmAz216iPjDuyEtF5tyt3pORbqgYCSgZhYPvAvE4Z0W86xz7u5G+9wPnOO7mwhkO+fSfY/VACt8j213zk1r7T3Dth1t60nKE6iW21bOOZbuOMiLS3fx8vLdFJZWkhIXzQUjejJtdC/OGJihKR8iIq1QAiqRyd+z4i30nFb9v/3sPFDOlqIythcdZlvRYbbvL2Nb0WFGHXidn0Yd23NaThxz0m5jR95l9O6WSO/uCfTpnkjv7olkJccRpYXORTqlICWgBiQ550rNLAZ4H7jNOfdxM/t/BxjjnLvBd7/UOdeuUq4R3Y6GoKe0uqaWjzfv58WlO3lt1R5KKqrpnhTLxSN7ctmoXozr3x2P2gERkeNoHVCJTO1ZQ66hZooekZZHjCeKfplJ9MtMOu5hd/8tWPGxy8UkUMm1JY8wdcPp7D1UecxjsdFR5HVLqE9M87ol0rtbInndEsjrlkD3pFgN7xWRZjnvWeC6+vsxvktLZ4avBe5u4fHOLUDVclsS7YnizMGZnDk4k3uvHMk76wqYu2wXzy3J5/FPtpOZHMfFI3tyyck5SkZFRNpACaiEv/ZW3IVm1yplyswWn2bNVFzMri3gk/89j4qqGnYeLGf7/sPk7z9M/oFydhw4zI795SzLP8jkigVcFv00vayQXS6THzCdpekXkJueQG63BHLTvYlpr/QEctLi6ZEar2FcIl2cmXmAxcAg4A/OuU+a2a8v0B94q8HmeDNbBFQDv3DO/aeZ584AZgD06dOn44IPNn+q5XbgkN24aA8XjOjJBSN6UlZZzYJ1+3hl+W6eWbyDf368jczkWM4f3oMLR/TkjIGZxEbr+11EpDEloNI5BaDnFCA+xuNb+qWJEW/Ln8bN/QdW7U1686yQ/4v6G4/Fp/Cfsokszz/IgcNVxzwlyuDLiZ/wHfcEWbUFFMdk8+mg71Ay+EqyU+LpkRpHVkocaQkx6kkV6aScczXAaDNLB14ws5HOuaYmuU/HO0e0psG2vs65nWY2AHjLzFY45zY18R5zgDngHYLb8UcRJK18Rx+n8ZDd4h3e+3DCPaZJcdFcNqoXl43qVZ+M/nfVXuYu3cUTC3eQEhfN2cOyOX94D84emqVquiIiPpoDKtLQicwvakPBpLLKanYdLGdXcQW7D5aTsuEFzt/4M2Ld0aG9h10sd1Z9k7m1Zx4NwWNkJMWRmRJLRlIc3ZNimVSxgPN2/ZmUyr2UJ+SQf+odVI24irSEGFITYkiJi1bSKnICQlEF18xmAoedc79p4rHPgFuccx8289xHgJedc8+29B4R3Y629zs6BMu7VFTV8OGmQl5buYf5a/ZRVHaEGI9xev8MppyUzXkn9aB398SAvLeISDjRHFCRtvC35xTaNDQsKS6awT1SGNwjxbvhg7+CO3ZeaaId4bcZL/KlK+9gX0klBSWVFJZWUui73l92hP67X+HiI38kwVcsKbF8F3nv38mdCzbWJ65RBinxMaTER9dfp8ZHM7liAVML/0Z61T5K4nqwZPCt7Os3jeS4aJLiokmK9fiuo0mK896Oi45SMisSAGaWBVQ55w6aWQJwPvDLJvYbBnQDPmqwrRveZLXSzDKBicCvghN5iLT3O9qfIbsnKD7Gw7nDenDusB7U1Do+236AN1bv5Y01e7nnpdXc89JqBmcnc+5J2Zw7NJtT+3bTVAwR6VKUgIo05s+cU2j/0DBo9kdQTMkuxg/IaP55998EVccWS0q0I/wy7T9MPPdmDpVXc6iiiuLyKkoqqimpqOJQRTVD9s7jmsMPEY836U2t3MO4FbO4c0n+MT2ujV0R/QE/jH6KnhRRYJn8K/l6lqSdR1JsNMlx0STHH71Oiff2vtYlvqkJ0aTGe3tlk2I9SmRFjpUDPOqbBxoFPO2ce9nMZgOLnHNzfftNB550xw5bOgn4i5nV+p77C+fc6mAGHxLt+Y72Z8huBy7x4okyxvbrzth+3bnrkpPYWljG/LX7mL9mL39/bwt/eWczKfHRTB6SxdlDsjhraBbZKfF+v5+ISCTQEFyRjuLP8F1/h4e1sMwMsw42/7xm3q86JY+NX/6IsspqyiprKKusprSymsNHasjZ8RLnrLuXmNqK+v0rLY4/Jt/Ka1GTKfXtW1pZfcz6qnWmRb3PD32FmXaTyUP2JT5OnkJ6YgzpCTF0S4ylW1Is3RtcMpNjyUyOIzM5jqQ4nSeT0AjFENxg6FLtaHu+l4O8xEtJRRUfbCzkrbX7WLCugIIS74nBEb1SOXtoFmcNyWZMn3T1jopIxNIQXJFA82f4rp/Vev3qbYVme1yjS3YyrGdq08+5/y/QIPkEiHOV3B71JLff/pP6bc45Kqtr63tcSyqqiVvzHIM+eZjoGu/zcynkHvsL/05OYX7MWRSWHmHDvlIOlB2h7EjNMe9Rl7gmWBH7ojJ5OvXrbOx5CTlp8eSkxdMr3VtNOK9bgoo0iUjT2vO93N4lXk6wtzQlPoaLRuZw0cgcnHOs3n2It9cV8M66Av78zmb+sGATKXHRTBiYweQhWZw1JEtzR0WkU1ACKtKR2jt81985p8FMXNs4h8rMiI/xEB/jISslzrvx2fuh5tjkNdZVcn35Y1x/0w+P2V5RVcOBw0coKj1C1MpnGLLwaOLa0xUwo/h3/PLwEf5xeDxHamqPeW5KfDS9uyXSp3si/TKTGJCZRP+sJPpnJpGhtVhFura2fi+3Z75oB1fXNTNG9EpjRK80bjlnEIcqqvhwYxHvrC/g3fUFvL56LwD9MhKZNDiLSYMzmTAwgxRV1hWRCKQEVCTU/JlzGszE1d/eVmjXD7r4GA85aQnkpCXA0w8cl7jGU8ndic/x/+6aRVHZEXYXl7PzQHn9Wqzb9x9m/b4S5q/dS1XN0aHA3RJjvIWfspMZ2jOF4TmpDMtJJVlDe0WkofZ817W3t7SdUuNjuGhkTy4a2RPnHJsKynh/QwHvbSjkuSX5/PPjbXiijNG905k4KJNJgzMZ3VvDdUUkMugXmEikClbi6m9vK3T4UGGK84mKMrJSvOujjspLP26XmqVP4ebfg6dkF6XxPXkl65s8c+QM5i7bRckn1fX79c1IZGSvNEblpTEqL52RuanqTRDpytrzXedPdV0/h+yaGYOykxmUncz1E/tzpLqWJdsP8P6GQt7bWMhDb23gwfkbSIr1cPqADM4YmMHEQZkM7ZFCVJRGf4hI+FECKtLVBGuYMAR/juvyp/G8clv9+6VU7Gb6nt8wfeqDuJOvZndxBWt2H2L1rkOs3n2IZfkHeWXFbgDMYEh2Cqf27cZpvku/jEQN3xXpKtrzXedPdd0OGrIbGx3F+AEZjB+QwR0XDqX4cBUfbS7kg41F9UWNADKSYpkwMIMzBnqH6+r7TETChargikhg+XPW399qlH5UFS4qrWTFzmKW7ShmyfYDfLb9AIcqvD2lWSlxjOvfnfH9uzN+QAaDspP1A64LURVcaVZ7v6P8rXhe917t+A7debCcDzcW8tGmIj7YVMjeQ97qujlp8YwfkMGEgRlMGJBBXrcEfZ+JSECpCq6IhEYw57j6MSwuIzmOs4dmc3bl27BiNo58qjN7sXDAd3jmyHg+2bKfV5Z7e0mzUuK8w9sGZjJxcCa56QntOy4R6Rza+x3lz5Bd8KvnNDc9gavH9ubqsb1xzrG5sIyPNhXx0aYi3l1fwAuf7azf7/QB3pNrp/fvTp/u6iEVkeBQD6iIdB7+9jK00JvhTr6aHfvL64e4fbipiMJSb4/CwKwkJg3OYvKQTMYPyCAxVuf0OhP1gEqH8fe76UR6TptQW+vYsK+UjzcX8fHmIj7Zsp/9ZUcA6Jkaz7j+3esvg7KSNYdURE6IekBFpPPzd85pCxUtbdQ19MlIpE9GH774uT445/0B996GQt5dX8CTn27nkQ+3EuuJYlz/7r4F5LM0XFdEjvL3u8nfntNmREUZQ3umMLRnCted0Y/aWsemglI+3rKfhVv28/HmIuYu2wVAemIMY/t247S+3Rnbrxsn56YRH+Px631FRBpSAioinUcQhu6aGUN6pDCkRwrfOLM/FVU1LNp6gHfW7+PtdQXc+8oa7n1lDbnpCZw9NItzhmZzxiD1jop0af5+N53IMljQ6vzRqCjzLlPVI4Wvju+Lc47t+w+zcMt+Pt26n0VbD/DmGm9Ro1hPFCNzUxnTpxun9unGmD7p9NI0BBHxg4bgioh0YIGQ/RPu5DWbzNvr9vH+xkIOH6kh1hPF6QO6c/bQbM4dlk3/zKTAHId0KA3BlZDztyDbiT63gaLSShZvO8Di7QdYvPUAK3YWU1ldC0CP1DhG907nlN7pjO6dzsm5aVrOSkTqNdeOKgEVEfH3h1orz6us9vaOLli7j7fXF7BxXyngXYP07CFZnD0smwkDMjSsLUwpAZWw4Of6oR09f7TOkepa1uw+5KsafpDl+QfZWnS4/vEBWUmcnJvGyblpjMxNY3ivVFKVlIp0SUpARURa4s+PvHb+wNux/zAL1nmH6n64qZCKqlrioqM4fUAGZw/J4qyhWQzITNLc0TChBFQi2qx0oKnfeAazDnboWx0oO8Ky/IOsyC9m+c5iVuQXs+dQRf3jfbonMjwnlZNyUhmWk8JJPVPJ65agIkcinZyKEImItMSf5WLaWSCkd/dEvjahH1+b0I+Kqho+2bKft9ft4531Bcx+eTW8DHndEpg0OIuzhmQyYWAmaQnqORARP5zo/NF26JYU613Oamh2/bZ9JRWs2nWI1b7Lql3F/Hf1Hur6PRJjPQzOTmZwjxSG9EhmcHYKA7OSye2WgEeJqUinpgRURMRf/v7AW/408fNnc1ZxPmel5cH5M9mRdxnvrC/gnfUFvLRsF08s3I4nyhjdO50zB2UyaXAmp/ROJ8YTFZhjEZHOxd/Kux0kOyWe7KHxnNMgKS2rrGb93hLKFz/J8NX3k1q4j72FGfz8s2v4v9ozAYiNjmJAZhL9M5Pom5FE/8xE+mV4b2enxKnXVKQTCNgQXDOLB94F4vAmus865+5utM/1wK+Bnb5NDznn/tbS62rokIiEDX/mjrbhOVU1tSzdcZB31xfw7oZCVuQfpNZBclw04wdkMGlwJhMHZTIwS8N1A0lDcCXi+Tt/NNAxNfoOdNEJbJnwf3yaeh6bCsrYtK+ULYVlbN9/mOrao79T46Kj6N09kb7dE8nrlkBet6PXvdLj6Z4Uq+9EkTAS9Dmg5v0GSHLOlZpZDPA+cJtz7uMG+1wPjHXOfbutr6uGU0TCSnt/4PlRGKT4cBUfbirk3Q2FvL+xgB37vT/cctLimTgokzMHZXLGoAyyU+I74ojERwmoSAC04zuwuqaWXQcr2FLkTUa3+663FR1m54FySiqrj9k/PiaKXukJ5KYnkJMWT06a79p3v2davAoiiQRR0OeAOm9mW+q7G+O7RFbFIxGR1rR37qgfC8unJcZw8ck5XOzeg22zcfH5HI7vyVNpX+fBNWN4drH3uUN7pHDmYO9w3dP7Z5AQq+q6IhJm2vEdGO2Jok9GIn0yEpt+SnkV+QcOk3+gnF0Hy9l5oJydB8vZVVzBuj0FFJRW0rifJSnWQ4+0eHLS4umRGk/PVG9i2vA6IzlO81BFAiigc0DNzAMsBgYBf3DOfdLEbl8ws8nAeuB259xxp8XMbAYwA6BPnz4BjFhEJMBOYN5o3bA1A5IqdnNDzf1cf/mDrMq4kPc3entH//nxNv7+/hZiPVGM7deNyUOyOGtIFsN6pmhomoiEXgcWR0pLiCEtIY0RvdKafPxIdS17D1Ww51AFu4sr2FNczq6DFfXbPt5UxN6SSmpqj81SPVFGdkoc2anx9EyN8yWnCfRMi6NnakJ9oqqTfCL+CcoyLGaWDrwAfMc5t7LB9gyg1DlXaWb/A3zROXduS6+loUMiEtH8XXO0jcPWyo/U8OnW/by3oYD3NhSydk8JANkpcUweksWUYdmcOThTi8W3gYbgigSAv9+BgYhj/mxccT61KbnsPO0HrMu+mD2HKthb7E1Q9x6qYE+x99J4uC94E+C6ob11Pai90hLISfcO/+2VHk9irOp9StcV0mVYnHMHzWwBcBGwssH2oga7/Q34VTDiEREJmbofWO0tDNLGYWsJsR4mD8li8pAsAPYequCd9QW8u76A11ft4dnF+cR4jPEDMpgyLJsLR/YkJy3hRI9KRKRt/P0O7EiNRpR4SvLp88Gd9JmaCOObjqOssro+Od3tS1D3NLheufMQhaWVxz0vPTGGnLQEctPj6+en5nbzXud1SyQzWYWTpOsJZBGiLKDKl3wmAK8Dv3TOvdxgnxzn3G7f7SuBHznnxrf0ujpzKyJdkh/Fi+o1ONN/JCmH13rM4MGCMWwqKANgdO90Lh7Zk0tOzqF396bnWnVFweoBPdGq8WZ2HfAT3/Z7nXOPtvR+akelyzuR79MWVFbXsLe4kt3F5ewurmBXsXdu6q6DFd45qgfLKak4tic1LjqqvpJv7+4J9OmeSO9uid5qvxmJGq0iES0UPaA5wKO+eaBRwNPOuZfNbDawyDk3F7jVzKYB1cB+4PoAxiMiErn8XdOv0Zn+uLJdXL7jl1w+9UE25VzCayv3MG/lbn4+by0/n7eWz/XrxhVjcrn05BzSE2MDekhSrxI4t2HVeDOb17BqvM9TjavGm1l34G5gLN5Cf4vNbK5z7kBQIheJRH4Ug2uLuGjP8UWTGlVKPzz1x2zPvZT8/d6EtK6I0o4Dh1m64yDF5VXHvGZGUix9M46uhdovM5EBmcn0z0oiOU7DeyUyBWUOaEfSmVsR6bL8WdOvjWf6d+w/zNxlu3jhs51s3FdKjMe4YERPvjq+L6f3794lh4iFYg6omSXiXbbs5oaF+5pbtszMrgXOds79j+/+X4C3nXNPNPceakelywtQD+hx/JjvWlxexY79h+uXm9m+v4ythYfZVlTGruKKY/btkRrHwKxkBmYlMyg7mcHZ3uuslLgu+Z0t4Sekc0BFRKQDtHfJF2jzmf7e3RO55ZxBfOvsgazadYjnluTz3OJ8Xlm+m8HZyXxlfF+uOi2PJJ1xD4gTqBqfCzT8JZ3v29b49VVNXqSOvyNK2mv+7GPfA7z3588+9ru8wcnFtLQ80qbMZGQT3/XlR2rYtr+MrYVlbCooY3NBGZsKSvnP0p3HDO1NS4hhaI8UhvRMZmiPFIblpDKsZ4qG80rY0C8JEZHOrJ1LHpgZI3PTGJmbxg8vHMZLy3fx+MfbuHvuKh54cz1fn9if687oR1qCfsh0JOdcDTC6rmq8mY1sWDUeeAl4okHV+EeBFqvGN3r9OcAc8PaAdlzkIhEoWIWQ2nICsHEvafEO7/2GcfokxHoY1jOVYT1Tj9nunKOgpJIN+0rZuK+U9XtLWLenhBeX7jomMe3dPYGTeqYyolcaJ+elMjI3jeyU+BM+TJH2UgIqItKZ+XOm33c2PqE4n2vS8rhmykwWp53HHxZs4r431jPn3c1cd0Zf/uesgaTqjHqH8qNq/E7g7AaP5QFvBzZKkU7AnxEl7dWWE4Bt7SWt08RUDBt1Ddmp8WSnxjNxUGb9rs45dhdXsHbPIdbsLmH17kOs2XWIN9bspW4GXo/UOEblpTO6dzqj8tIYlZeuE4wScEpARUQ6s/ae6W/mbPxpUx/k4euvYdWuYv749ib++PYmnvp0B9+/YCjXjO2NJ0rzjfzVRNX484FfNtqnvmo8MA1Y47v9X+D/zKyb7/4FwF1BCFtEWtOWE4DtKYjUjt5S8I5o6ZWeQK/0BM4d1qN+e2llNat3HWLFzmJW7ixmWf5B3li9t/7xwdnJnNa3G6f27cbYvt3on5mkOaXSoVSESEREjmpjcY4V+cXc89IqFm07wEk5qcyaOpzTB2QEMdDAC+IyLKPwDqltWDV+dsOq8Wb2c7yJZ13V+Judc2t9z78B+F/fy/3MOfePlt5P7ahIELVWPK49BZH8LZ7UhgJ2xYerWL7zIEu3H2TJ9gMs2X60Im9mchzj+ndjXL/ujB+YwZDsFKJ00lHaoLl2VAmoiIgcNSsd72oejRnMOnjMFuccr6zYzc9fXcvOg+V888z+3HHhUOJjPEEINPBCUQU3GNSOioSR9lTKbcf3s1+v30BtrWNzYSmLth5g4Zb9fLJlPzsPel8jIymW8QMzOGNgBpMGZR277IxIA6qCKyIirWtH0SIz47JRvTh3WDb/9+oa/vb+Ft7bUMj9XxzN8F6px7+GiIgcqz3TJNpZVK7+ddszx9QnKsoYlJ3CoOwUpo/zVs7esf8wH28u4qNNRXywqZBXlntnBfTLSGTykCwmD87ijEEZJMYqvZCWqQdURESO8vNsOcCCdfv44bPLOXj4CD+5dDjXndEvsLEGmHpARSSs+PP97G+vaSsJsXOOzYVlvLe+gHc3FPLRpiLKq2qIjY5iwoAMzh2WzbnDsundXb2jXZl6QEVEpHX+Lk+w/GnOmT+bhVX5FMVnMfvlq9hS+BX+32XDVaBIRKQj+PP93N5e0zYWOjIzBmYlMzArmesn9qeyuoZFWw/w1tp9vLV2H3fPXcXdc1dxUk4qF47owYUjejKsZ4qKGQmgHlARETlRTZyVPxIVzx0VN1A25PM8eO0YkuIi73ynekBFJOK1t9fU30JHjWwuKOXNNXt5fdVeFm8/gHPQNyORS07O4dKTcxjRK1XJaBegHlAREQmMJuYYxdZW8H+pLzBq3Zlc85eP+MfXP6cFz0VEgq29vabtWRYGmh2uOyArmRlZycyYPJCCkkreXLOXeSv3MOfdzfzp7U30zUhk6qheXDGmF4OyU078OCWiqAdUREROTAtzjBZMX88t/15Cn+6JPDVjAmmJkbPAuXpARaTLaU8PqB9zUveXHeH1VXt4efluPtxUSK2DEb1SuWJ0LpeP6aUTlZ1Mu5dhMbMH2/C6h5xzPznR4NpDDaeISJhp5QfL+xsKueGRTxmZm8q/vnl6xFRIbGsCGq7tZXPUjopIs9qTVJ7gcN19JRW8tGw3Ly7dyfL8YjxRxllDsrjqtDymnJRNXHTnWNKrK/MnAd0GzGzlde90zp3UAfG1mRpOEZEw04YfLK+t3M23Hl/CmYOz+NvXxhIbHRWiYNuuHQloWLaXzVE7KiItakMVXMC/CrvN2LivlOeX5PP8kp3sOVRBemIMnx+Tx7XjejO4h4boRip/5oDe75x7tJUX7XbCkYmISGRrwxyji0bm8IvPj+KHzy3n9qeX8uD0MZ2pOq7aSxHpPEZd03rlc/BvXdLGfMnuoOJ8fpiWxx2XzuT9hHN46tMd/PPjrTz8wRbG9u3GteP6cOmoHOJj1CvaGbTUA/o559ynQY6nVTpzKyISuf7yziZ+Pm8tP7hwKLecMyjU4bSoHT2gYdleNkftqIh0iBNYN7otzy8sreT5Jfk8sXAHWwrL6JYYwxc/14cvn95H64tGiOba0ZbGQM0xsw1m9lMzGx7A2EREpIuYMXkAl47K4b431rN42/5Qh9NR1F6KSNcz6hpvspjWGzDvdVuTT2iygjpV5d7tQGZyHDMmD+St75/F4988nXH9uzPn3U1M/vUCvvnoIj7aVESkFVMVr2aH4DrnxpjZUGA68KyZVQFPAE8657YGKT4REekMfMOsrDif36fmMjvpKm59Ip5Xb5tEWkLkVMZtitpLEemy2jpctyltXPLFzJg4KJOJh9+iet89eEp2sntLBr9Ydw0/7XEJN5zZn6mn5KhoUQRpsQqEc26dc+4e59xw4GtAGjDfzD4ISnQiIhL56oZZFe8AHFGH8vl/tX9hXOmb3PX88k5xBlvtpYhIOzU3V7Sp7b52JLokH8PRi0Lui3+YSRVvccczy5j8qwX85Z1NHKqoCmzM0iHaVIbQzKKAbKAHkATsC2RQIiLSiTQxzMpTU87spOd5dcUenljYRBGLCKX2UkSkjabM9M75bCgmwbu9sSbakejaCu6MfYbHbhjH4OwUfj5vLRN//hY/n7eGfSUVAQxcTlSLi7GZ2STgWuAKYAXwJHC7c6448KGJiEin0Mwwq+TKPUwanMnsl1dxzrAsctISmtwvEqi9FBFppzZUUK/XTDtixflMHpLF5CFZrNxZzF/e3cxf393MIx9s5dpxfZgxeQC90iO3bemsmk1AzWwHsA1vIzrLOaezuCIi0n7NlOq3tDz+78qTmXLfO/zqtXXc/8XRwY+tA6i9FBHxUwcu+TIyN43fXzuG758/hD+9vYnihY/jFj+FsyJqUnKJPv9u/+erSodqaQjumc65M51zD6kxFRERv7UwzKp390RunNSfFz7byZLtB0IT34lTeykiEkjtGK7bLzOJXw5Zy2/jHybXCjEc0SX5HHnh2xz4+F9BClha0lIC+vXWnmxmszouFBER6ZRaKdV/89mDyEqJY/ZLqyO1IJHaSxGRQGrvki/zZxNVfeyc0VhXyeF5d3PPS6soKq0MfMzSrJbmgH7TzA618LjhLTk/q0MjEhGRzqeFYVbJcdH84MKh/PDZ5cxdtovLR+cGObgTpvZSRCTQ2rPkSzNzRntZEY99tI1nFuVz46QBfHNSf5LiWiyJIwHQUg/oX4GUFi7Jvn1EREROyFWn5jEyN5VfzFvL4SPVoQ6nvdReioiEk2aWeLG0PP773cmcOSiT+99cz1m/XsC/Pt5GdU1tkAPs2ppN+Z1z9wQzEBER6YKWP+0dKlWcz3NJvfhByeXMebc33z1vSKgja7MTaS/NLB54F4jD2yY/65y7u9E+3wO+CVQDBcANzrltvsdq8FbdBdjunJvmbywiIp3GlJne9acbLt3imzM6KDuZP3/1NJZsP8Av5q3lJ/9ZyY53HuF2e5L4w7tbrsYrHaJN64CKiIh0ON/C4t7Kho64sp38Ou7v7H7/McoqI64X1F+VwLnOuVOA0cBFZja+0T6fAWOdc6OAZ4FfNXis3Dk32ndR8ikiAm2aM3pqn248NWM8L07eye3lfyD+8C7Aedukl271tlESEBr0LCIiodHEwuJxrpLv1D7Bc0u+ydcm9AtNXEHkvFWXSn13Y3wX12ifBQ3ufgx8JTjRiYhEsDbMGTUzTln3IN5zgQ1UlVP75j1EqRc0INQDKiIiodFckYioIv7xwVZqayOyIm67mZnHzJYC+4A3nHOftLD7N4B5De7Hm9kiM/vYzK5o4T1m+PZbVFBQ0CFxi4h0Cs20RRzayeur9gQ3li6i1QTUzIaY2XwzW+m7P8rMfhL40EREpFNrpkhERUJPthSWsWBdZC2p6W976Zyrcc6NBvKAcWY2spnX/wowFvh1g819nXNjgS8BD5jZwGbeY45zbqxzbmxWVlb7DkxEpDNrpi0qsExm/HMxN/1zMfsOVXg3Ln8a7h8Js9K91xqm65e29ID+FbgLqAJwzi3HW05eRETEf80sLB574T3kpMXz9/e3hCYu/51Qe+mcOwgsAC5q/JiZnQf8GJjmnKts8JydvuvNwNvAGL+jFxHpipppizIu/xk/umgYC9btY8p97/DhC3/CNahboLmi/mtLApronFvYaFuXqQ4hIiIB0kyRiOjRX+RrE/rx4aYiVu9qaXnNsNPu9tLMssws3Xc7ATgfWNtonzHAX/Amn/sabO9mZnG+25nARGD1iR6EiEiX0kJbdPPZA3ntu5MZnpNKn89+gzWqW0BVubeegbRLW4oQFfqG9DgAM7sK2B3QqEREpGtopkjEl8b14cH5G3j4gy385upTQhCYX/xpL3OAR83Mg/ek8NPOuZfNbDawyDk3F++Q22TgGTODo8utnAT8xcxqfc/9hXNOCaiISHu1ULCof2YST9w4Hptd1PRzm5tDKs1qSwJ6CzAHGGZmO4EtqAKfiIgEgm9d0LTifD6Mz+any65m30VDyU6JD3VkbdHu9tI3TPe4YbPOuZkNbp/XzHM/BE4+kYBFRKR1UVHmnStavOP4B5uZQyrNa3UIrnNus6/xywKGOefOdM5tDXhkIiLStTRaF7Rb1V7u9cxhyctzQh1Zm6i9FBHpxJqYK1pOHGtH3H78vipW1KJWe0DN7HuN7gMUA4udc0tbeF488C4Q53ufZ51zdzfaJw54DDgNKAK+qMZaRKSLamJd0EQ7wpj1DwK3hiamdvC3vRQRkQhQN0R3/mwozqcquRf3V3+Rvy7oyS1uHd89bzDRnqijJ1Pr2rO6YkUNX6OLa8sQ3LG+y0u++5cBy4GbzOwZ59yvmnleJXCuc67UzGKA981snnPu4wb7fAM44JwbZGbTgV8CX/TrSEREJLI1M48mq7aQDXtLGNwjJcgBtZu/7aWIiESCBnNFY4DvHqnm4NxVPLRgI59sKeLBa8eQ08TJ1PpiRUpAgbZVwc0DTnXOfd859328vZXZwGTg+uae5LxKfXdjfJfGq4pfDjzqu/0sMMV8p4xFRKSLaWYezS4yeGVFRNS+86u9FBGRyJQYG82vrjqFB744mtW7DjH19+/jmitKpGJF9dqSgGbj7c2sUwX0cM6VN9p+HDPzmNlSYB/whnPuk0a75AI7AJxz1XiHKmU08TozzGyRmS0qKChoQ8giIhJxmlmL7bn0G3g1MhJQv9tLERGJXFeMyeU/t0wkJT6GXe64VMZLxYrqtSUBfRz4xMzuNrO7gQ+Af5tZEq2sN+acq3HOjcZ7VnicmY30J0jn3Bzn3Fjn3NisrCx/XkJERMJdM2uxpZ/+ZdbvLWXD3pJQR9gav9tLERGJbIN7pPCfWybycuY3Oexij30wJsF7klWANswBdc791MxeA87wbbrJObfId/vLbXkT59xBM1sAXASsbPDQTqA3kG9m0UAa3mJEIiLSFTWxFtvFhyqY9dIqXlmxm++G8TzQjmgvRUQkcqUlxHDjLXcy74lETln3IL2iinCpuXjOu1vzPxtoSw8ozrlPgSeAF4B9ZtanteeYWZaZpftuJwDnA2sb7TYXuM53+yrgLedc43miIiLShWWnxvO5ft0jYhiuP+2liIh0HlFRxqVfvo3Fn3+XYVVPcIH7AzvyLjt+xy68VEurCaiZTTOzDXgX1H7Hdz2vDa+dAywws+XAp3jngL5sZrPNbJpvn78DGWa2EfgecKc/ByEiIp3bZaNywn4Y7gm0lyIi0slcPjqXf35jHIWlR7jyjx+wbMfBow82Wve6fqmWLpKEtqUH9KfAeGC9c64/cB7wcctPAefccufcGOfcKOfcSOfcbN/2mc65ub7bFc65q51zg5xz45xzm0/gWEREpJO6aGRPzAj3arh+tZciItI5nT4gg+duPoP4GA/T53zMext8xVRbWqqlC2hLAlrlnCsCoswsyjm3AO86ZyIiIoHTYHhS9t/Gcnv20nAfhqv2UkREjjEoO5kXvjWRvhmJfOORRbyxem/zS7J0kaVa2pKAHjSzZOBd4HEz+x1QFtiwRESkS2tieNK3Sh5kWMFr4TwMV+2liIgcJysljidnjOekXqnc9K/FHE7IaXrHLrJUS1sS0MuBw8DtwGvAJqCJmbQiIiIdpInhSdG1Ffww+mk+2hy2xdLVXoqISJPSE2N5/Junc1rfbtx16AqqPfHH7tCFlmppSwI60zlX65yrds496px7EPhRoAMTEZEurJlhSL2iili6/WBwY2k7tZciItKs5LhoHv36OA4MvJLvV9xAWUIODde9rl+qpZNXyG1LAnp+E9su7uhARERE6jUzDGl/dBZLG1YSDC9qL0VEpEUJsR7mfPU0CvpdzskHf8u8L6yB21cem3x28gq5zSagZnazma0AhprZ8gaXLcDy4IUoIiJdzpSZ3uFIDcUksGTQrWwuLOPg4SOhiasJai9FRKQ94mM8/PVrYxnTpxu3PvkZC9buO/pgF6iQ21IP6L+BqcBc33Xd5TTn3FeCEJuIiHRVo67xDkdK603D4UnJY68FCLdeULWXIiLSLklx0fzj659jWM9U/udfi/lwY6H3gS5QIbelBNQDHAJuAUoaXDCz7oEPTUREurRR13iHJc06WD88aVTvdMzCLgFVeykiIu2WGh/DYzeMo19GIv/zz8Ws21PSfCXcTlQhN7qFxxYDznfbGj3mgAEBiUhERKQZyXHRDMlOCbcEVO2liIj4pVtSLP/4+jiu/MMH3PDIp7x6zl2kvfn9Y4fhdrIKuc32gDrn+jvnBvgu/Rtd1JiKiEhIjO6dztIdB3HOtb5zEKi9FBGRE5GbnsDD13+OA4eP8JWFfam8+IHjpqDUFynqBFrqAa1nZtOAyb67bzvnXg5cSCIiIs0b3SedpxbtYGvRYfpnJoU6nGOovRQREX+MzE3j99eO4cbHFnHLikH85bYVeKIaD6rBWw13/mzvnNC0PG/PaIQlp60uw2JmvwBuA1b7LreZ2f8FOjAREZGmjOmTDsDSHQdCG0gj/raXZhZvZgvNbJmZrTKze5rYJ87MnjKzjWb2iZn1a/DYXb7t68zswg48JBERCaIpJ/XgnmkjeHPNXn7933XH79BJlmhpSw/oJcBo51wtgJk9CnwG/G8gAxMREanX4Izv0LQ8ro69nKXb+3LlmLAqyuBve1kJnOucKzWzGOB9M5vnnPu4wT7fAA445waZ2XTgl8AXzWw4MB0YAfQC3jSzIc65mo49NBERCYavTujH2j0l/PmdTYzunc5FI3sefbClJVoiqBe01R5Qn/QGt9MCEIeIiEjTGp3xteId3Ov5K6kbXgh1ZE1Jb3C7Te2l8yr13Y3xXRpPcL0ceNR3+1lgipmZb/uTzrlK59wWYCMwzs/YRUQkDMycOpxTeqdzxzPL2FxQevSBTrJES1sS0J8Dn5nZI76zuYuBnwU2LBEREZ8mzvjGuUqmlzxCRVVYdfT53V6amcfMlgL7gDecc5802iUX2AHgnKsGioGMhtt98n3bGr/+DDNbZGaLCgoK2ndUIiISVHHRHv745VOJjY7ipn8t5vCRau8DnWSJlmYTUDP7g5lNdM49AYwHngeeAyY4554KVoAiItLFNXNmN4ciVu06FORgjtcR7aVzrsY5NxrIA8aZ2ciOjNE5N8c5N9Y5NzYrK6sjX1pERAIgNz2BB6ePYeO+Uu58boW38vuUmd4lWRqKwCVaWuoBXQ/8xsy2ArcDO5xzc51ze4ISmYiICDR7ZneXywiX9UA7rL10zh0EFgAXNXpoJ9AbwMyi8Q7vLWq43SfPt01ERCLcmYMz+f4FQ5m7bBdPL9rhnec59cGIX6KlpXVAf+ecmwCchbeRe9jM1prZ3WY2JGgRiohI19bMGd+/xnwlLBLQE20vzSzLzNJ9txOA84G1jXabC1znu30V8JbzLoQ6F5juq5LbHxgMLOyI4xIRkdC7+ayBTBiQwT0vrWZbUZk32bx9Jcw66L2OsOQT2jAH1Dm3zTn3S+fcGOBa4ApgTaADExERAZo941swYFpYLcVyAu1lDrDAzJYDn+KdA/qymc32rSsK8Hcgw8w2At8D7vS95yrgabzLvrwG3KIKuCIinUdUlPGba07BE2V87+ll1NQ2qlG3/Gm4fyTMSvdeR8CSLK0uw+Ib6nMx3jLvU4C3gVkBjUpERKShUdccd5Z31IFNvLpiD4cqqkiNjwlRYEf5214655YDY5rYPrPB7Qrg6mae/zNUHFBEpNPKTU/gp5eP5LtPLeXP72zilnMGeR+oqxJfV6ivbl1QCOue0ZaKEJ1vZg/jrah3I/AKMNA5N90592KwAhQREWlK96RYAA6VV4U0DrWXIiISaJeP7sVlo3K4/431rNxZ7N3Y0rqgYaylIbh3AR8CJznnpjnn/u2cKwtSXCIiIi1KjvMO4imtrA5xJGovRUQksMyMe68YSUZyLN99aimV1TURuy5oS0WIznXO/c05Fz4TbERERHzqEtCyECegai9FRCQY0hNj+cUXRrFxXylz3tkcseuCtlqESEREJBwl+RLQkoqQ94CKiIgExTlDs7n05Bx+v2AjBeN+FJHrgioBFRGRiJQSX9cDqqKvIiLSdcycOpxYTxTfWzsEF4HrgioBFRGRiJQUJkNwRUREgqlHajx3XDCE9zYU8gpnHrsuKIT9sixKQEVEJCIlx/qG4CoBFRGRLuarE/pxcm4a97y0mkMVvmrwdcuyFO8A3NFlWcIsCVUCKiIiESkpzgOoB1RERLoeT5TxsytHUlRayX2vr/dujJBlWZSAiohIRIr2RBEfExUOy7CIiIgE3ai8dK4d14d/fbyNLYVlEbMsixJQERGJWMlx0UpARUSky/rueUOIjY7i1/9dGzHLsigBFRGRiJUcF60huCIi0mVlpcRx46QBvLpiD1tHfz8ilmVRAioiIhErKS6aUq0DKiIiXdiNkweQmRzLD9cPi4hlWaJDHYCIiIi/NARXRES6uuS4aG6bMpj/9+Iq3pp8FlPqlmMJU+oBFRGRiJUcF03ZESWgIiLStU0f14f+mUn88rW11NS6UIfTIiWgIiISsTQEV0REBGI8UfzgwqGs31vKc0vCq+ptY0pARUQkYiXHR1NaWRPqMERERELu4pE9GZmbyp/f3nS0F3T503D/SJiV7r1e/nRIY4QAJqBm1tvMFpjZajNbZWa3NbHP2WZWbGZLfZfwKtEkIiJhzTsHtCrUYYiIiIScmXHTWQPZXFjG66v2eJPNl26F4h2A816/dGvIk9BA9oBWA993zg0HxgO3mNnwJvZ7zzk32neZHcB4RESkk0mOi6aiqpbqmtpQhyIiIhJyF4/MoW9GIn9+ZxNu/myoKj92h6pymB/alCtgCahzbrdzbonvdgmwBsgN1PuJiEjXkxTnLeZedkTDcEVERDxRxozJA1iWXwzFzcwFbW57kARlDqiZ9QPGAJ808fAEM1tmZvPMbEQzz59hZovMbFFBQUEgQxURkQiSHOcB0FIsIiIiPl84NY/M5DiKPFlN75CWF9yAGgl4AmpmycBzwHedc4caPbwE6OucOwX4PfCfpl7DOTfHOTfWOTc2K6uZf0gREelykuNiACiL0AS0jfUSftCgVsJKM6sxs+6+x7aa2QrfY4uCfwQiIhJu4mM83HBmP2aXX0VtdMKxD8YkwJTQlt0JaAJqZjF4k8/HnXPPN37cOXfIOVfqu/0qEGNmmYGMSUREOo8kXw9oSeQuxdJqvQTn3K/raiUAdwHvOOf2N9jlHN/jY4MWtYiIhLUvn96Xt2LO5rHM70Fab8C811MfhFHXhDS26EC9sJkZ8HdgjXPuvmb26Qnsdc45MxuHNyEuClRMIiLSuSTXzQGN0B5Q59xuYLfvdomZ1dVLWN3MU64FnghSeCIiEqHSEmL48vg+zH63mnPvWEifjMRQh1QvkD2gE4GvAuc2GDp0iZndZGY3+fa5ClhpZsuAB4HpzjkXwJhERKQTSY6P7AS0oVbqJWBmicBFeEcW1XHA62a22MxmtPDaqqUgItLF3DCxP2bGk59uD3UoxwhYD6hz7n3AWtnnIeChQMUgIiKdW1KstxkrifAEtJV6CXWmAh80Gn57pnNup5llA2+Y2Vrn3LuNn+icmwPMARg7dqxO9IqIdAE9UuM5e0gWzy3J5/sXDMUT1WJqFjRBqYIrIiISCCmdoAe0tXoJDUyn0fBb59xO3/U+4AVgXKDiFBGRyHP12Dz2Hqrk3Q0FsPxpuH8kzEr3Xi9/OiQxKQEVEZGIVbcOaGmEFiFqS70E335pwFnAiw22JZlZSt1t4AJgZWAjFhGRSHLusB50T4ply1v/gJduheIdgPNev3RrSJLQgA3BFRERCbQYTxSx0VGUHonMBJSj9RJWmNlS37b/BfoAOOf+7Nt2JfC6c66swXN7AC94c1iigX87514LRtAiIhIZYqOjuGJ0Lhcuugms/NgHq8ph/uygV8VVAioiIhEtJS46YofgtqVegm+/R4BHGm3bDJwSkMBERKTTuOZzeeQsKmz6weL84AaDhuCKiEiES4qLjtghuCIiIoE2rGcqhZ6sph9MywtuMCgBFRGRCJccF01pZU2owxAREQlb60bczmEXe+zGmASYMjPosSgBFRGRiJYcwUNwRUREgmHUxTP4Se0MDsb0AAzSesPUB4M+/xM0B1RERCJcUpyHwtIjoQ5DREQkbKUlxlA1/CrOWn8WC38yhbhoT8hiUQ+oiIhEtOT4GPWAioiItGLaKb0oLq9i8bYDIY1DCaiIiES05DgPJUpARUREWjRhYAbRUca765upiBskSkBFRCSiJcVqDqiIiEhrkuOiObVvN97bUBDSOJSAiohIREuOj+bwkRpqal2oQxEREQlrZw3JYtWuQxSWVoYsBiWgIiIS0ZLjvPX0yo6oF1RERKQlkwZnAvDBxtANw1UCKiIiEa0+AdUwXBERkRaN6JVGt8QY3lkfumG4SkBFRCSiJfkS0NIKJaAiIiIt8UQZEwdl8t6GQpwLzdQVJaAiIhLR6npAS9UDKiIi0qrJQ7IoKKlk3d6SkLy/ElAREYloyfF1Q3BrQhyJiIhI+KubB/peiJZjUQIqIiIRLSm2rge0KsSRiIiIhL+ctAQGZyfzboiWY1ECKiIiES0lvi4BVQ+oiIhIW0wanMUnW/ZTURX8tlMJqIiIRLQkVcEVERFpl8lDMjlSXcvCLfuD/t5KQEVEJKIlxXkAFSESERFpq9P7ZxDrieK9EAzDVQIqIiIRLS7aQ6wnSgmoiIhIGyXEevhc/268v7Eo6O+tBFRERCJeUpxH64CKiIi0w8jcNDbtK6WmNrjrgSoBFRGRiJcUF605oCIiIu0wIDOJIzW17DxQHtT3VQIqIiIRLzkuOiKH4JpZbzNbYGarzWyVmd3WxD5nm1mxmS31XWY2eOwiM1tnZhvN7M7gRi8iIpFsQFYyAJsKS4P6vtFBfTcREZEAiNQEFKgGvu+cW2JmKcBiM3vDObe60X7vOecua7jBzDzAH4DzgXzgUzOb28RzRUREjjMgMwmALQVlnDM0eO+rHlAREYl4yfGROQTXObfbObfEd7sEWAPktvHp44CNzrnNzrkjwJPA5YGJVEREOpvuSbGkxkezOcg9oEpARUQk4iXFRVMSgQloQ2bWDxgDfNLEwxPMbJmZzTOzEb5tucCOBvvk0/bkVUREujgzY0BWMpsLyoL6vkpARUQk4iXHRmYPaB0zSwaeA77rnDvU6OElQF/n3CnA74H/+PH6M8xskZktKigI/ppvIiISngZkJrGlUAmoiIhIu3iH4NaEOgy/mFkM3uTzcefc840fd84dcs6V+m6/CsSYWSawE+jdYNc837bjOOfmOOfGOufGZmVldfgxiIhIZBqQlcTu4goOHwneSVwloCIiEvGSfEWIaoO8ltmJMjMD/g6scc7d18w+PX37YWbj8LbdRcCnwGAz629mscB0YG5wIhcRkc6gf6a3Em4we0FVBVdERCJecpwHgMNVNSTHRVTTNhH4KrDCzJb6tv0v0AfAOfdn4CrgZjOrBsqB6c45B1Sb2beB/wIe4GHn3Kogxy8iIhFsQJa3Eu7mgjJG9EoLyntGVCstIiLSlOS4GABKK6ojKgF1zr0PWCv7PAQ81MxjrwKvBiA0ERHpAvrXLcUSxB5QDcEVEZGIl+TrAY3QtUBFRERCIj7GQ256ApsLgrcUixJQERGJeCnx3l7PSK6EKyIiEgoDsoJbCVcJqIiIRLykWG8Cqh5QERGR9umfmcTmgjK85QUCL2AJqJn1NrMFZrbazFaZ2W1N7GNm9qCZbTSz5WZ2aqDiERGRzispTgmoiIiIPwZkJlFSWU1h6ZGgvF8ge0Crge8754YD44FbzGx4o30uBgb7LjOAPwUwHhER6aQ0BFdERMQ/A7K8S7EEax5owBJQ59xu59wS3+0SYA2Q22i3y4HHnNfHQLqZ5QQqJhER6ZzUAyoiIuKfukq4m4M0DzQoc0DNrB8wBvik0UO5wI4G9/M5PknFzGaY2SIzW1RQUBCwOEVEJDIlKwEVERHxS256ArHRUUErRBTwBNTMkoHngO865w758xrOuTnOubHOubFZWVkdG6CIiES8uOgooqOM0goloCIiIu0RFWX0z0iK/CG4AGYWgzf5fNw593wTu+wEeje4n+fbJiIi0mZmRlJctOaAioiI+GFAVlLkD8E1MwP+Dqxxzt3XzG5zga/5quGOB4qdc7sDFZOIiHReyXHRlFbWhDoMERGRiNM/M4ntRYepqqkN+HtFB/C1JwJfBVaY2VLftv8F+gA45/4MvApcAmwEDgNfD2A8IiLSiXkT0KpQhyEiIhJxBmQlU13ryD9QXl+UKFACloA6594HrJV9HHBLoGIQEZGuIynOQ5l6QEVERNptQJavEm5BacAT0KBUwRUREQm05PgYSjQHVEREpN0G+JLOYFTCVQIqIiKdQnKcR0WIRERE/JCeGEv3pFg2FQQ+AQ3kHFAREZGgOWdoNkN7pIY6DBERkYj0tQl96d0tMeDvowRUREQ6havH9m59JxEREWnSd88bEpT30RBcERERERERCQoloCIiIiIiIhIUSkBFREREREQkKJSAioiIiIiISFAoARUREREREZGgUAIqIiIiIiIiQaEEVERERERERIJCCaiIiIiIiIgEhTnnQh1Du5hZAbCtg14uEyjsoNcKFR1D+OgMx6FjCA+d4Rgg8o+jr3MuK9RBdDS1oy3S8YQ3HU940/GEr1AdS5PtaMQloB3JzBY558aGOo4ToWMIH53hOHQM4aEzHAN0nuOQ5nW2z1jHE950POFNxxO+wu1YNARXREREREREgkIJqIiIiIiIiARFV09A54Q6gA6gYwgfneE4dAzhoTMcA3Se45DmdbbPWMcT3nQ84U3HE77C6li69BxQERERERERCZ6u3gMqIiIiIiIiQaIEVERERERERIKiSyagZnaRma0zs41mdmeo42mNmW01sxVmttTMFvm2dTezN8xsg++6m2+7mdmDvmNbbmanhijmh81sn5mtbLCt3TGb2XW+/TeY2XVhcAyzzGyn77NYamaXNHjsLt8xrDOzCxtsD9nfm5n1NrMFZrbazFaZ2W2+7RHzWbRwDJH2WcSb2UIzW+Y7jnt82/ub2Se+mJ4ys1jf9jjf/Y2+x/u1dnwhPIZHzGxLg89itG972P09SccI5f8lf3VUuxQuOvL7PRx05HdkODEzj5l9ZmYv++5H7PFYBP4ebYmZpZvZs2a21szWmNmESD0eMxvaoA1eamaHzOy7YXs8zrkudQE8wCZgABALLAOGhzquVmLeCmQ22vYr4E7f7TuBX/puXwLMAwwYD3wSopgnA6cCK/2NGegObPZdd/Pd7hbiY5gF3NHEvsN9f0txQH/f35gn1H9vQA5wqu92CrDeF2vEfBYtHEOkfRYGJPtuxwCf+P6Nnwam+7b/GbjZd/tbwJ99t6cDT7V0fCE+hkeAq5rYP+z+nnTpkL+DiGtHfXGfcLsUTpeO+n4Pl0tHfUeG2wX4HvBv4GXf/Yg9HiLw92grx/Mo8E3f7VggPZKPp8FxeYA9QN9wPZ6u2AM6DtjonNvsnDsCPAlcHuKY/HE53v84+K6vaLD9Mef1MZBuZjnBDs459y6wv9Hm9sZ8IfCGc26/c+4A8AZwUcCD92nmGJpzOfCkc67SObcF2Ij3by2kf2/Oud3OuSW+2yXAGiCXCPosWjiG5oTrZ+Gcc6W+uzG+iwPOBZ71bW/8WdR9Rs8CU8zMaP74Aq6FY2hO2P09SYeIyHa0g9qlsNGB3+9hoQO/I8OGmeUBlwJ/8903Ivh4mhGRf29mlob3pNTfAZxzR5xzB4nQ42lkCrDJObeNMD2erpiA5gI7GtzPp+Ufs+HAAa+b2WIzm+Hb1sM5t9t3ew/Qw3c7nI+vvTGH67F82zdc4eG6oQxEwDH4hvOMwXtWOSI/i0bHABH2WfiGYi0F9uFNujYBB51z1U3EVB+v7/FiIIMQH0fjY3DO1X0WP/N9FvebWZxvW9h+FnJCOtPnF4lt6XFO8Ps9bHTQd2Q4eQD4IVDru59BZB9PZ/k9Ct4RRAXAP3xDpP9mZklE7vE0NB14wnc7LI+nKyagkehM59ypwMXALWY2ueGDzjlHy70QYScSY/b5EzAQGA3sBn4b0mjayMySgeeA7zrnDjV8LFI+iyaOIeI+C+dcjXNuNJCHtxdpWGgjar/Gx2BmI4G78B7L5/AOq/1R6CIU8U+kfBc21hm+3+t0hu/IOmZ2GbDPObc41LF0oM70ezQa75D8PznnxgBleIeo1ouw4wHAN6d4GvBM48fC6Xi6YgK6E+jd4H6eb1vYcs7t9F3vA17A+6W8t66r3He9z7d7OB9fe2MOu2Nxzu31NZC1wF85OvQxbI/BzGLw/jh53Dn3vG9zRH0WTR1DJH4WdXzDfBYAE/AOe4luIqb6eH2PpwFFhMlxNDiGi3xDAZ1zrhL4BxH0WYhfOtPnF4ltab0O+n4POyf4HRkuJgLTzGwr3mHq5wK/I3KPpzP9HgVvj19+g1E8z+JNSCP1eOpcDCxxzu313Q/L4+mKCeinwGDzViGLxdtNPTfEMTXLzJLMLKXuNnABsBJvzHWVI68DXvTdngt8zVfdajxQ3KDrPdTaG/N/gQvMrJtveOUFvm0h02h8/JV4PwvwHsN081ax6w8MBhYS4r833/yRvwNrnHP3NXgoYj6L5o4hAj+LLDNL991OAM7HO2drAXCVb7fGn0XdZ3QV8Jbv7GVzxxeqY1jboHEzvPNLGn4WYfX3JB0iotrRVkRiWwp06Pd7WOjA78iw4Jy7yzmX55zrh/f/yFvOuS8TocfTyX6P4pzbA+wws6G+TVOA1UTo8TRwLUeH30K4Ho8Lg2pNwb7grfy0Hu/cgh+HOp5WYh2At8LgMmBVXbx45wXMBzYAbwLdfdsN+IPv2FYAY0MU9xN4h0VW4T3L9A1/YgZuwFtkZSPw9TA4hn/6YlyO9z9vToP9f+w7hnXAxeHw9waciXe4xXJgqe9ySSR9Fi0cQ6R9FqOAz3zxrgRm+rYPwJtAbsQ7ZCbOtz3ed3+j7/EBrR1fCI/hLd9nsRL4F0crWYbd35MuHfa3EDHtaIOYO6RdCpdLR36/h8OlI78jw+0CnM3RKrgReTxE6O/RVo5pNLDI9zf3H7xV2SP5eJLw9pqnNdgWlsdjviBEREREREREAqorDsEVERERERGREFACKiIiIiIiIkGhBFRERERERESCQgmoiIiIiIiIBIUSUBEREREREQkKJaAiIiIiIiISFEpARcKQmWWY2VLfZY+Z7fTdLjWzPwbg/R4xsy1mdlML+0wys9VmtrKj319ERKSjqA0VCW9aB1QkzJnZLKDUOfebAL7HI3gXyX62lf36+fYbGahYREREOoraUJHwox5QkQhiZmeb2cu+27PM7FEze8/MtpnZ583sV2a2wsxeM7MY336nmdk7ZrbYzP5rZjlteJ+rzWylmS0zs3cDfVwiIiKBpjZUJDwoARWJbAOBc4FpwL+ABc65k4Fy4FJfA/p74Crn3GnAw8DP2vC6M4ELnXOn+F5bRESks1EbKhIC0aEOQEROyDznXJWZrQA8wGu+7SuAfsBQYCTwhpnh22d3G173A+ARM3saeL6jgxYREQkDakNFQkAJqEhkqwRwztWaWZU7Oqm7Fu//bwNWOecmtOdFnXM3mdnpwKXAYjM7zTlX1JGBi4iIhJjaUJEQ0BBckc5tHZBlZhMAzCzGzEa09iQzG+ic+8Q5NxMoAHoHOE4REZFwozZUJADUAyrSiTnnjpjZVcCDZpaG9//8A8CqVp76azMbjPfs73xgWUADFRERCTNqQ0UCQ8uwiIhKyIuIiPhJbahI+2gIrogAFAM/bW0RbeAloDBoUYmIiIQ/taEi7aAeUBEREREREQkK9YCKiIiIiIhIUCgBFRERERERkaBQAioiIiIiIiJBoQRUREREREREguL/A318AGEzHnXoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 936x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4))\n",
    "\n",
    "# plot the 1C results\n",
    "t_sol = solutions[0][\"Time [s]\"].entries\n",
    "ax1.plot(t_sol, solutions[0][\"Terminal voltage [V]\"](t_sol))\n",
    "ax1.plot(voltage_data_1C[:,0], voltage_data_1C[:,1], \"o\")\n",
    "ax1.set_xlabel(\"Time [s]\")\n",
    "ax1.set_ylabel(\"Voltage [V]\")\n",
    "ax1.set_title(\"1C\")\n",
    "ax1.legend([\"DFN\", \"Experiment\"], loc=\"best\")\n",
    "\n",
    "# plot the 5C results\n",
    "t_sol = solutions[1][\"Time [s]\"].entries\n",
    "ax2.plot(t_sol, solutions[1][\"Terminal voltage [V]\"](t_sol))\n",
    "ax2.plot(voltage_data_5C[:,0], voltage_data_5C[:,1], \"o\")\n",
    "ax2.set_xlabel(\"Time [s]\")\n",
    "ax2.set_ylabel(\"Voltage [V]\")\n",
    "ax2.set_title(\"5C\")\n",
    "ax2.legend([\"DFN\", \"Experiment\"], loc=\"best\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a 1C discharge we observe an excellent agreement between the model and experiment, both in terms of the overall shape of the curve and the capacity. The agreement between model and experiment is less good at 5C, but in line with other implementations of the DFN (e.g. [[6]](#References)). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "The relevant papers for this notebook are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] Joel A. E. Andersson, Joris Gillis, Greg Horn, James B. Rawlings, and Moritz Diehl. CasADi – A software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1):1–36, 2019. doi:10.1007/s12532-018-0139-4.\n",
      "[2] Marc Doyle, Thomas F. Fuller, and John Newman. Modeling of galvanostatic charge and discharge of the lithium/polymer/insertion cell. Journal of the Electrochemical society, 140(6):1526–1533, 1993. doi:10.1149/1.2221597.\n",
      "[3] Madeleine Ecker, Stefan Käbitz, Izaro Laresgoiti, and Dirk Uwe Sauer. Parameterization of a Physico-Chemical Model of a Lithium-Ion Battery: II. Model Validation. Journal of The Electrochemical Society, 162(9):A1849–A1857, 2015. doi:10.1149/2.0541509jes.\n",
      "[4] Madeleine Ecker, Thi Kim Dung Tran, Philipp Dechent, Stefan Käbitz, Alexander Warnecke, and Dirk Uwe Sauer. Parameterization of a Physico-Chemical Model of a Lithium-Ion Battery: I. Determination of Parameters. Journal of the Electrochemical Society, 162(9):A1836–A1848, 2015. doi:10.1149/2.0551509jes.\n",
      "[5] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, and others. Array programming with NumPy. Nature, 585(7825):357–362, 2020. doi:10.1038/s41586-020-2649-2.\n",
      "[6] Giles Richardson, Ivan Korotkin, Rahifa Ranom, Michael Castle, and Jamie M. Foster. Generalised single particle models for high-rate operation of graded lithium-ion electrodes: systematic derivation and validation. Electrochimica Acta, 339:135862, 2020. doi:10.1016/j.electacta.2020.135862.\n",
      "[7] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). ECSarXiv. February, 2020. doi:10.1149/osf.io/67ckj.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "pybamm.print_citations()"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}