docs/source/examples/notebooks/creating_models/3-negative-particle-problem.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A step towards the Single Particle Model"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the [previous notebook](./2-a-pde-model.ipynb) we saw how to solve a PDE model in pybamm. Now it is time to solve a real-life battery problem! We consider the problem of spherical diffusion in the negative electrode particle within the single particle model. That is,\n",
    "$$\n",
    "  \\frac{\\partial c}{\\partial t} = \\nabla \\cdot (D \\nabla c),\n",
    "$$\n",
    "with the following boundary and initial conditions:\n",
    "$$\n",
    "  \\left.\\frac{\\partial c}{\\partial r}\\right\\vert_{r=0} = 0, \\quad \\left.\\frac{\\partial c}{\\partial r}\\right\\vert_{r=R} = -\\frac{j}{FD}, \\quad \\left.c\\right\\vert_{t=0} = c_0,\n",
    "$$\n",
    "where $c$ is the concentration, $r$ the radial coordinate, $t$ time, $R$ the particle radius, $D$ the diffusion coefficient, $j$ the interfacial current density, $F$ Faraday's constant, and $c_0$ the initial concentration. \n",
    "\n",
    "In this example we use the following parameters:\n",
    "\n",
    "| Symbol | Units              | Value                                          |\n",
    "|:-------|:-------------------|:-----------------------------------------------|\n",
    "| $R$      | m                | $10 \\times 10^{-6}$                            |\n",
    "| $D$      | m${^2}$ s$^{-1}$ | $3.9 \\times 10^{-14}$                          |\n",
    "| $j$      | A m$^{-2}$       | $1.4$                                          |\n",
    "| $F$      | C mol$^{-1}$     | $96485$                                        |\n",
    "| $c_0$    | mol m$^{-3}$     | $2.5 \\times 10^{4}$                            |\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up the model\n",
    "As before, we begin by importing the PyBaMM library into this notebook, along with any other packages we require, and start with an empty `pybamm.BaseModel`\n"
   ]
  },
  {
   "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[plot,cite]\" -q    # install PyBaMM if it is not installed\n",
    "import pybamm\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "model = pybamm.BaseModel()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then define all of the model variables and parameters. Parameters are created using the `pybamm.Parameter` class and are given informative names (with units). Later, we will provide parameter values and the `Parameter` objects will be turned into numerical values. For more information please see the [parameter values notebook](../parameterization/parameter-values.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = pybamm.Parameter(\"Particle radius [m]\")\n",
    "D = pybamm.Parameter(\"Diffusion coefficient [m2.s-1]\")\n",
    "j = pybamm.Parameter(\"Interfacial current density [A.m-2]\")\n",
    "F = pybamm.Parameter(\"Faraday constant [C.mol-1]\")\n",
    "c0 = pybamm.Parameter(\"Initial concentration [mol.m-3]\")\n",
    "\n",
    "c = pybamm.Variable(\"Concentration [mol.m-3]\", domain=\"negative particle\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we define our model equations, boundary and initial conditions, as in the previous example. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# governing equations\n",
    "N = -D * pybamm.grad(c)  # flux\n",
    "dcdt = -pybamm.div(N)\n",
    "model.rhs = {c: dcdt}\n",
    "\n",
    "# boundary conditions\n",
    "lbc = pybamm.Scalar(0)\n",
    "rbc = -j / F / D\n",
    "model.boundary_conditions = {c: {\"left\": (lbc, \"Neumann\"), \"right\": (rbc, \"Neumann\")}}\n",
    "\n",
    "# initial conditions\n",
    "model.initial_conditions = {c: c0}"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we add any variables of interest to the dictionary `model.variables`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.variables = {\n",
    "    \"Concentration [mol.m-3]\": c,\n",
    "    \"Surface concentration [mol.m-3]\": pybamm.surf(c),\n",
    "    \"Flux [mol.m-2.s-1]\": N,\n",
    "}"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the model"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to discretise and solve the model we need to provide values for all of the parameters. This is done via the `pybamm.ParameterValues` class, which accepts a dictionary of parameter names and values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "param = pybamm.ParameterValues(\n",
    "    {\n",
    "        \"Particle radius [m]\": 10e-6,\n",
    "        \"Diffusion coefficient [m2.s-1]\": 3.9e-14,\n",
    "        \"Interfacial current density [A.m-2]\": 1.4,\n",
    "        \"Faraday constant [C.mol-1]\": 96485,\n",
    "        \"Initial concentration [mol.m-3]\": 2.5e4,\n",
    "    }\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here all of the parameters are simply scalars, but they can also be functions or read in from data (see [parameter values notebook](../parameterization/parameter-values.ipynb))."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As in the previous example, we define the particle geometry. Note that in this example the definition of the geometry contains a parameter, the particle radius $R$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = pybamm.SpatialVariable(\n",
    "    \"r\", domain=[\"negative particle\"], coord_sys=\"spherical polar\"\n",
    ")\n",
    "geometry = {\"negative particle\": {r: {\"min\": pybamm.Scalar(0), \"max\": R}}}"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both the model and geometry can now be processed by the parameter class. This replaces the parameters with the values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "param.process_model(model)\n",
    "param.process_geometry(geometry)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now set up our mesh, choose a spatial method, and discretise our model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "submesh_types = {\"negative particle\": pybamm.Uniform1DSubMesh}\n",
    "var_pts = {r: 20}\n",
    "mesh = pybamm.Mesh(geometry, submesh_types, var_pts)\n",
    "\n",
    "spatial_methods = {\"negative particle\": pybamm.FiniteVolume()}\n",
    "disc = pybamm.Discretisation(mesh, spatial_methods)\n",
    "disc.process_model(model);"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model is now discretised and ready to be solved."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solving the model"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As is the previous example, we choose a solver and times at which we want the solution returned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2021-11-19 15:29:22,931 - [WARNING] processed_variable.get_spatial_scale(520): No length scale set for negative particle. Using default of 1 [m].\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6AAAAEYCAYAAABCw5uAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABw/UlEQVR4nO3dd3hVVdbH8e9KIaH3XqT3KkUQ6RbEXrGMXbGgYp1BneI7o2MvY28oWLEPdkSkiDSDIL0X6b2EnrLeP+6JRiYkN5Dk5ia/z/Oc5567T7nr5EJ21tn77G3ujoiIiIiIiEh+i4l0ACIiIiIiIlI8KAEVERERERGRAqEEVERERERERAqEElAREREREREpEEpARUREREREpEDERTqAglalShWvX79+pMMQEZEoMmPGjC3uXjXScUSK6k4REcmtw9WdxS4BrV+/PklJSZEOQ0REooiZrYp0DJGkulNERHLrcHWnuuCKiIiIiIhIgVACKiIiIiIiIgVCCaiIiIiIiIgUiGL3DKiIiIiIiBQ/KSkprFmzhv3790c6lCIlMTGROnXqEB8fH9b+SkBFRERERKTIW7NmDWXLlqV+/fqYWaTDKRLcna1bt7JmzRoaNGgQ1jHqgisiIiIiIkXe/v37qVy5spLPPGRmVK5cOVetyvmWgJpZXTMbZ2YLzGyemQ0Jyu83s7VmNitYBmQ65h4zW2pmi8zslEzlHc1sTrDtGQv+1ZhZgpm9H5RPM7P6+XU9IiIiIiIS3ZR85r3c/kzzswU0FbjT3VsAXYHBZtYy2PaUu7cPlq8Agm0XAa2A/sALZhYb7P8iMAhoEiz9g/JrgO3u3hh4CngkH69HREQkX2Vz8/aC4H26mXXKtP+lmW7ozgq2tw+2jQ9u6GZsqxaU6+atiIhETL49A+ru64H1wXqymS0AamdzyFnASHc/AKwws6VAFzNbCZRz9ykAZvYmcDbwdXDM/cHxHwHPmZm5u+f9Ff1u1Ky1TFy8hScubJefHyMiIsVPxs3bn82sLDDDzMYAc4FzgZcz7+zu7wDvAJhZG2CUu8/KtMul7p50yGf8dvPWzC4idPN2YL5cTSZ//ugXAMomxlM2Me6313LBernfykPvS8TpKSERkaKoQAYhCu6udgCmAd2Bm83sciCJUEW7nVByOjXTYWuCspRg/dBygtfVAO6eamY7gcrAlkM+fxChFlTq1at31Nezbsd+Pv55DUNPbU7VsglHfT4RERE4/M1bdx8DOXZzuhh4L4yPicjN23nrdrF190GS96ew52BajvsnxMUEiWncHxLWjPWqZROoVaEktSskUqtCSaqVTSQ2Rl3rRKTw2rFjB++++y433XRTWPs/99xzPP300yxbtozNmzdTpUoVIDTwz5AhQ/jqq68oVaoUw4cP59hjjwXgm2++YciQIaSlpXHttdcydOhQALZt28bAgQNZuXIl9evX54MPPqBixYoMHz6cu+++m7POOovXXnstyzj27dtHt27dmD9/PuvWrfstjiOV7wmomZUBPgZuc/ddZvYi8C/Ag9cngKuBrGoNz6acHLb9XuD+CvAKQKdOnY66gj2+UWUAJi/bwlnts2vUFREROTKH3LwNx0BCyWVmb5hZGqF6+IEgyYzIzdsvb+3x23pqWjq7D6SSvD+VXftTSN6fGiwpf3jd9YftKWzYtZ/k/Sns2pfKvpQ/JrFxMUaN8olBUlqSWhUSqV2hVPBakloVSlI6QYP/i0jk7NixgxdeeCHsBLR79+6cfvrp9O7d+w/lX3/9NUuWLGHJkiVMmzaNG2+8kWnTppGWlsbgwYMZM2YMderUoXPnzpx55pm0bNmShx9+mH79+jF06FAefvhhHn74YR55JPT04sCBA3nuuecOG0fJkiWZNWsW9evXP9JL/4N8/U1sZvGEKr133P0TAHffmGn7q8AXwds1QN1Mh9cB1gXldbIoz3zMGjOLA8oD2/L+Sv6ode3yVCgVzzdzNygBFRGRPHfozdsw9j8O2OvuczMVX+rua4OuvB8DlwFvEqGbt5nFxcZQoVQJKpQqccTnSN6fwvqd+1m7Yx/rduxj7fbQ67od+5m+Yhsbdu0nLf2PYVcoFU+t8iX/0HJaq0JJGlQpTeNqZUiMjz3Mp4lIUfN/n89j/rocf73mSsta5fjHGa0Ou33o0KEsW7aM9u3bc9JJJ/HYY49le74OHTpkWT5q1Cguv/xyzIyuXbuyY8cO1q9fz8qVK2ncuDENGzYE4KKLLmLUqFG0bNmSUaNGMX78eACuuOIKevfu/VsCmtm8efO46qqrOHjwIOnp6Xz88cc0adIkzJ9AePItAQ1Gqh0GLHD3JzOV1wy6GAGcQ+i5FoDPgHfN7EmgFqHBhqa7e5qZJZtZV0J3gS8Hns10zBXAFOB84Pv87kIEEBtjXNS5Hq9MXMbijck0rV42vz9SRESKiaxu3obhIg7pfuvua4PXZDN7F+hCKAGNyM3bvBbqkht/2Do4Ld3ZuGt/KDkNEtOM9TXb9zJtxVaS96f+tn+MQf0qpWleoyzNa5SjWY2yNK9RlroVSxGjrr0ikgcefvhh5s6dy6xZs0hOTqZ9+/ZZ7vfuu+/SsmXLLLcBrF27lrp1f2+3q1OnDmvXrs2yfNq0UCeajRs3UrNmTQBq1qzJpk2bsjz3Sy+9xJAhQ7j00ks5ePAgaWk5PzKRW/nZAtqd0N3WOWY2Kyi7F7g4GKHPgZXA9QDuPs/MPgDmExqEYbC7Z1zxjcBwoCShwYe+DsqHAW8FAxZtI1QBF4jrejRg5E+/cu8nc/jg+m6qnERE5Kgd7uZtDsfEABcAPTOVxQEV3H1LkNCeDnwXbI7IzduCFhtjv7VwdjrMPrv2p7Buxz6Wb97Dwg3JLFy/i3nrdvH13A1k/ERKlYilafVQMtosSE6b1yhLxdJH3norIpGXXUtlQShbtiyzZs06omOz+pVtZoctz41u3brx4IMPsmbNGs4999w8b/2E/B0FdxJZd/P5KptjHgQezKI8CWidRfl+QpVugatcJoH7BrTg7o9m89LEZdzUu3EkwhARkaLlcDdvEwj1/qkKfGlms9w9Y77snsAad1+e6TwJwOgg+YwllHy+GmyL2M3bwqZcYjzlasTTvEY5BrSp+Vv5ngOpLN6YzKINySzcEHodPW8DI39a/ds+1com0KxGWVrULEez6qHkVN14RSRcycnJ9OjRI8ttObWA1qlTh9Wrf/99tGbNGmrVqsXBgwezLAeoXr0669evp2bNmqxfv55q1aplee5LLrmE4447ji+//JJTTjmF1157jb59+x7JJR6WnsY/Cud3rMOExZt5bPQiGlUtwymtakQ6JBERiWLZ3LwF+PQwx4wnNN925rI9QMfD7B+xm7fRonRCHB3qVaRDvYq/lbk7m5MP/JaQLtiwi0Ubkhk+eSUHU9OBUKtrk2pl6FS/Ip3rV6Jz/UrUqlAyUpchIoVM2bJlSU5O/m39SFtAzzzzTJ577jkuuugipk2bRvny5alZsyZVq1ZlyZIlrFixgtq1azNy5Ejefffd344ZMWIEQ4cOZcSIEZx11qFj1oUsX76chg0bcuutt7J8+XJmz56tBLQwMTMev6Adq7fv47aRs3j72i50PKZSpMMSERGRPGZmVCuXSLVyifRsWvW38tS0dFZu3fNbYjpr9Q4+/Xktb0/9FYDaFUrSuX5FOtWvRJcGlWhctYwe2xEppipXrkz37t1p3bo1p556ao6DED3zzDM8+uijbNiwgbZt2zJgwABee+01BgwYwFdffUXjxo0pVaoUb7zxBgBxcXE899xznHLKKaSlpXH11VfTqlWoq/HQoUO58MILGTZsGPXq1ePDDz/M8jPff/993n77beLj46lRowZ///vf8/aHAFgRfOwjW506dfKkpEPn5D46m5L3M/DlqWzatZ8RV3ehU30loSIiRYmZzXD3wz1KWOTlR91ZlKWmpbNwQzI/rdzGTyu3MX3FdrbsPgCERuLtdEwoIe1cvyJtalegRFxMhCMWKR4WLFhAixYtIh1GoTJ8+HCSkpKynYYlQ/369UlKSspyHtCsfraHqzvVApoHqpVNZOSgrlz8ylQuf306z19yLH2aZ92vWkRERIq2uNgYWtcuT+va5bmqewPcnVVb9/LTym0krdzOTyu38d2C0AiUCXExtKtbgS71K9GpfkWOPaYi5RLjI3wFIlJclCxZkq+//pprr72W1157Lct99u3bR7du3UhJSSEm5uhvmKkFNA9t2rWfq4b/xMINyfzrrNZcctzRT9wtIiKRpxZQtYDmtS27D5C0chs/BQnpvHW7SEt3Ygya1yhH5/oV6dqwMic0qUJZJaQieWLBggU0b9481yPDSvbcnYULF6oFNBKqlUvkg+u7Mfjdn7n30zks27yboac2Jz5WXWtERETkd1XKJNC/dU36tw6NvrvnQCqzVu9g+optJK3axgdJaxgxZRVxMUaXBpXo27wafZtXo2HVMhGOXCR6JSYmsnXrVipXrqwkNI+4O1u3biUxMTHsY9QCmg9S09J54MsFDJ+8ki71K/HcJR2oVi78L0VERAoXtYCqBbSgpaSlM/PXHYxduJFxCzexeONuAOpXLkXf5tXp27waXRpU0vOjIrmQkpLCmjVr2L9/f6RDKVISExOpU6cO8fF/7K1xuLpTCWg+GjVrLUM/nkOZxDiev+RYujTQ4EQiItFICagS0EhbvW0v4xZt4vuFm5i8bCsHU9MpkxDHCY2r0LdFNfo0q0bVsgmRDlNE5DdKQAMFXYku2pDMjW/PYNW2vdzStzE392lMnLrkiohEFSWgSkALk70HU5m8dCtjF25i3MJNbNgVas1pV6c8fZpXo1/z6rSqVU7TvYhIRCkBDUSiEk3en8I/Rs3jk5lr6VCvAk8PbM8xlUsXaAwiInLklIAqAS2s3J3563cxbuEmxi7cxKzVO3CHqmUT6NusGn2aV+OEJlUok6BhP0SkYCkBDUSyEv38l3Xc9+kc0tKdf5zZigs61tED0CIiUUAJqBLQaLF19wEmLN7M2IWbmLh4M8n7UykRG0PXRpU5s10tTmlVXaPqikiBUAIaiHQlum7HPu74YBZTl2/jxBbVeODsNtQorwGKREQKMyWgSkCjUUpaOkkrtzNu0Sa+nrue1dv2kRAXw4ktqnNW+1r0alaVhLjYSIcpIkWUEtBAYahE09KdN35cwePfLiI+Noa/ntaCCzvVVWuoiEghpQQ08nWnHB135+dfd/DZrLV8MXs9W/ccpHzJeAa0qcFZ7WvTpX4lPTMqInlKCWigMFWiK7fsYegns5m6fBsnNK7CQ+e2oW6lUpEOS0REDqEEtPDUnXL0UtLSmbR0C5/NWsfoeRvYezCNmuUTObNdLc5sX4uWNcvppriIHDUloIHCVommpzvvTv+Vh75agAN3n9KMy7vVJ1Z3IUVECg0loIWr7pS8s/dgKt8t2MSomWuZsHgzqelOk2plOKt9Lc5qX1s3xkXkiCkBDRTWSnTtjn3c88kcJi7eTJva5fn3OW1oU6d8pMMSERGUgBbWulPy1rY9B/lqznpGzVrLTyu3A3BsvQqc3aE2p7WpSeUymmdURMKnBDRQmCtRd+eL2ev55xfz2br7AJd3q88dJzelnEarExGJKCWghbfulPyxZvtePvtlHaNmrmPRxmRiY4weTapwdvvanNSyOqU1rYuI5EAJaCAaKtFd+1N4YvQi3py6iqplEvjHGa0Y0KaGnscQEYkQJaCFv+6U/LNwwy7+O3Mdn81ay7qd+ylVIpazO9TmT8cdQ8ta5SIdnogUUkpAA9FUif6yegf3fjqHeet20atpVf55ViuOqVw60mGJiBQ7SkCjp+6U/JOe7iSt2s4HSav5/Jd1HEhNp+MxFbms6zGc2qaGpnQRkT9QAhqItko0NS2dt6au4olvF3MwLZ0bejXixl6NKFlCv+RFRAqKEtDoqjsl/+3Ye5CPZqzh7amrWLl1L5VKl+DCTnW59Lh6GrhIRAAloL+J1kp0w879/PurBXz2yzpqVyjJ305vwSmt1C1XRKQgKAGNzrpT8l96uvPjsi28NWUV3y3YiAN9mlXjT13r0atpNY3qL1KMKQENRHslOnX5Vu7/bB4LNyRzQuMq3H9mSxpXKxvpsEREijQloNFdd0rBWLdjHyOn/8p7P61mc/IB6lQsyaXHHcOFnepoBF2RYuhwdWdMDgftymFJNrPF+Re2HKprw8p8ccsJ/N+ZrZi9Zgf9n/6BB76YT/L+lEiHJiIiR8nM6prZODNbYGbzzGxIUH5B8D7dzDpl2r++me0zs1nB8lKmbR3NbI6ZLTWzZyzoMmNmCWb2flA+zczqF/iFSpFUq0JJ7ji5GT/+pS/PXdKBOhVL8sg3C+n20Pfc/v4sZqzaRnFr+BCR/5VtAgosc/dy2SxlgT1ZHZhNJfqYmS00s9lm9qmZVQjKVYmGKS42hiuOr8+4u3pzfsc6DPtxBX0en8DHM9aQnq5f7CIiUSwVuNPdWwBdgcFm1hKYC5wLTMzimGXu3j5YbshU/iIwCGgSLP2D8muA7e7eGHgKeCR/LkWKqxJxMZzethYjB3VjzO09ubhLXb6bv5HzXpzCgGcm8e60X9lzIDXSYYpIhOSUgJ4XxjkOt8/hKtExQGt3bwssBu7JdIwq0VyoXCaBh89ry39v6k6diiW588NfOP+lycxduzPSoYmIyBFw9/Xu/nOwngwsAGq7+wJ3XxTuecysJlDO3ad4qMnpTeDsYPNZwIhg/SOgX8aNXZG81qR6Wf7vrNZMvbcf/z6nDQD3fjqHrv8eyz9GzWXJxuQIRygiBS3bBNTdl+d0gsPtk00l+q27Z9z2mgrUye78qkRz1q5uBT658XgePb8tv27byxnPTWLox7PZnHwg0qGJiMgRCnr1dACm5bBrAzObaWYTzKxHUFYbWJNpnzVBWca21QBBfbwTqJxXcYtkpXRCHJccV4+vbj2Bj288nhNbVue96as56amJXD38J5JWbot0iCJSQHJ6BrS5mX1tZl+aWSMzG25mO8xsupm1CPdDsqlErwa+zvRelegRiokxLuxUl+/v6s3V3Rvw0Yw19Hl8PC9NWMaB1LRIhyciIrlgZmWAj4Hb3H1XNruuB+q5ewfgDuBdMysHZHUzNuMZjey2ZY5hkJklmVnS5s2bc3cBIodhZnQ8piJPDWzPlHv6csdJTZm1egfnvzSFC1+awrhFm/ScqEgRl1MX3FeAF4C3ge+Bb4CKwL+A58L5gMNVomZ2H6Fuuu8ERapE80C5xHj+dnpLvr29J8c1qMTDXy/k5KcmMnreBv1CFxGJAmYWT6jefMfdP8luX3c/4O5bg/UZwDKgKaGbtZl7GNUB1gXra4C6wWfFAeWB/2l+cvdX3L2Tu3eqWrXq0V2USBYql0ng1n5NmPSXPvzjjJas2b6Xq974iQHPTOKzX9aRmpYe6RBFJB/klICWdffP3f09IMXdR3rI54QS0WwdrhI1syuA04FLg261qkTzWMOqZRh2ZWfevLoLJWJjuP6tGVz62jQWrM/uRrqIiERS8BjJMGCBuz8Zxv5VzSw2WG9IaJyE5e6+Hkg2s67BOS8HRgWHfQZcEayfD3zvukMpEVSqRBxXdW/A+Lv78Nj5bTmYmsat782k35MTeGfaKvanqCeXSFGSUwIam2n90IqwRHYHHq4SNbP+wF+AM919b6ZyVaL5oGfTqnw9pAf/PKsV89fv4rRnfuDeT+ewdbeeDxURKYS6A5cBfTONCj/AzM4xszVAN+BLMxsd7N8TmG1mvxAaC+EGd8+4EXsj8BqwlNBN3YxHXoYBlc1sKaEeR0ML5MpEclAiLoYLOtVlzO29eOlPHalQMp77Pp1Lz0fH8fKEZezWyLkiRYJll6+Z2fWEWi93H1LeGLjZ3W/L5tgTgB+AOUBGH4p7gWeABGBrUDbV3W8ws/OAfxLqlpsG/CNoaSWY82w4UJJQBXqLu7uZJQJvEXq+dBtwUU4DJxXnybR37D3I098t4a2pqygVH8ut/ZpwxfH1KRGX030IEZHizQ4zmXZxUZzrTokcd2fysq28OH4Zk5ZuoVxiHFccX58rj69P5TIJkQ5PRHJwuLoz2wS0KFIlCks3JfOvLxYwYfFmGlQpzX0DWtCvRTWK4QDCIiJhUQKqulMi65fVO3hx/DJGz99AQlwMF3Wux3U9G1K7QslIhyYih3G4ujPXTV9m9nPehCSR0rhaWUZc3YU3rupMjMG1byZx+evTWay5uERERKQQale3Ai9d1pExt/fk9La1eHvqKno9Oo47PpiluURFokyuW0DNbGYwUm1U0l3cP0pJS+etKat4+rvF7DmYxiVd6nH7SU2pVDrbR3xFRIoVtYCq7pTCZe2Ofbz2w3Lem/4r+1PSOblldW7q05j2dStEOjQRCeRZCyjwZR7EI4VEfGwMV58QGnnu0uPq8e70X+n92Dhen7SCFA1/LiIiIoVQ7Qol+ccZrZg8tB+39m3M1OVbOfv5H7lm+E8a8V+kkMtVC2gwL2dcxvtMI+1FDd3Fzd7ijcn864v5/LBkCw2rlOZePR8qIqIWUNWdUsjtPpDKiMkreSkYLffs9rW5/cSm1KtcKtKhiRRbR9UCambXm9lGYDaQBMwIXqWIaVq9LG9e3YVhV3SC4PnQPw2bxvx1upsoIiIihVOZhDgG92nMD3/uw/U9G/HVnPX0e3I8/xg1l83JmnpOpDAJqwXUzJYA3dx9S/6HlL90Fzd8KWnpvDN1FU+PXcLOfSkM7FSXO05uSrWyiZEOTUSkQKkFVHWnRJcNO/fzzPdLeP+n1STExXDNCQ24rmdDyiXGRzo0kWLjaJ8BXQbszduQpLCLj43hyu4NmHBXH67u3oCPf15Dn8fG8/y4pexPSYt0eCIiIiJZqlE+kX+f04bv7uhF3+bVePb7pfR8dByvTlyuv2FEIizcFtAOwBvANOC3fgzufmv+hZY/dBf3yK3YsoeHvlrAt/M3UrtCSf7cvxlntqul50NFpMhTC6jqToluc9fu5NHRi5i4eDM1yydy24lNOO/YOsTFHsl4nCISjqNtAX0Z+B6YSuj5z4xFipEGVUrzyuWdePe64yhfMp4hI2dx7ouTmbFqe6RDExERETms1rXL8+bVXXj3uuOoXi6Rv3w8h5OfnsjXc9aT2ykJReTohNsCOtndjy+AePKd7uLmjbR05+Of1/DY6EVsTj7AGe1q8Zf+zahTUaPNiUjRk5sWUDOrFMZu6e6+4+iiKjiqO6UocXe+nb+Rx0YvYumm3bStU54/n9KcE5pUiXRoIkXK4erOcBPQB4FVwOf8sQuupmEp5vYcSOXlCct4eeJyHLj2hAbc1KcxZRLicjxWRCRa5DIB3Q+sA7J7PiHW3evlSXAFQHWnFEVp6c4nP6/h6e+WsHbHPro3rsyfT2lOu7oVIh2aSJFwtAnoiiyK3d0b5kVwBUmVaP5Yt2Mfj36zkP/OWkeVMgncdXJTLuhUl9gYPR8qItEvlwnoTHfvcLT7FCaqO6UoO5CaxjtTf+W5cUvZtucgp7auwZ0nN6NxtTKRDk0kqh1VAlqUqBLNX7NW7+CBL+aTtGo7zWuU5a+ntVSXFhGJerlMQBPdff/R7lOYqO6U4mD3gVRe+2E5r05czr6UNC497hjuPLkpFUqViHRoIlHpaAchyuqENY4uJCmK2tetwIc3dOP5S45l94FU/jRsGtcM/4llm3dHOjQRkQKRVWJ56HOh0ZR8ihQXZRLiuO3Epkz8cx/+1PUY3pm2ij6Pj+fdab+Sll68GmxE8tPRjD09LM+ikCLFzDitbU2+u6MXf+nfnGkrtnHKUxO5/7N5bN9zMNLhiYjkKzPrbmYLzGyemR1nZmOAJDNbbWbdIh2fiGSvcpkE/nlWa768tQdNqpfl3k/ncPbzP2rUf5E8oi64ku+27D7Ak2MWM3L6r5RNjOfWfk24rOsxlIjT3FsiEh1y2QV3OnANUIbQ4H1nu/skMzsWeNbdu+djqPlCdacUV+7O57PX8+8vF7Bh137OPbY2Q09tTrWyiZEOTaTQO6IuuGZWKbsl/8KVoqRKmQT+fU4bvhrSg7Z1yvOvL+ZzytMT+XbeBs29JSJFUby7z3H3KcBmd58E4O4/AyUjG5qI5IaZcWa7Woy9sxc39W7EF7+sp+/jE3h14nJS0tIjHZ5IVMqpCWoGkBS8HrroVqjkSvMa5Xjz6i68cWVnYgwGvTWDS16dxrx1OyMdmohIXspct95zyDaNZiIShUonxPHn/s0ZfXtPujSoxINfLaD/0xP5YcnmSIcmEnXUBVciIiUtnfem/8pTYxazY18KF3Ssw10nN6NaOXVpEZHCJ5ddcM8EvnP3vYeUNwLOc/dH8yPG/KS6U+SPxi7YyD+/mM+qrXs5pVV1/npaS+pWKhXpsEQKlaOehiWoUHsGb8e7+xd5GF+BUSVauOzcl8Jz3y9h+OSVxMfGcGOvRlzXsyGJ8bGRDk1E5De5SUAPc3wNd9+QlzEVJNWdIv9rf0oawyat4Lnvl5Luzg29GnFj70b6G0YkcFTTsJjZw8AQYH6wDDGzh/I2RCmOypeM577TWjLm9l70bFKVJ8Yspu/j4/nvzLWka8hzESk6vop0ACKStxLjYxncpzFj7+zFSS2r85+xS+j3xAS+mbteY1yIZCPcYUgHACe5++vu/jrQHzgt/8KS4qZ+ldK8dFlHRg7qSqUyJbjt/Vmc8+JkZqzaFunQRETygoW1k1ldMxuXaRqXIUH5BcH7dDPrlGn/k8xshpnNCV77Zto23swWmdmsYKkWlCeY2ftmttTMpplZ/Ty+VpFipVaFkjx3ybG8d11XyibGccPbP3PZsOks3ZQc6dBECqXczINRIdN6+TyOQwSArg0r89ngE3j8gnZs2LmP816cwuB3f2b1tr05HywiUni9GuZ+qcCd7t4C6AoMNrOWwFzgXGDiIftvAc5w9zbAFcBbh2y/1N3bB8umoOwaYLu7NwaeAh7J/eWIyKG6NarMF7ecwP1ntGT2mh30f/oHHvxyPsn7UyIdmkihEm4C+hAw08yGm9kIQqPg/jv/wpLiLCbGOL9jHcbd1Ztb+zVh7IKN9HtyAo98s1C/xEUkqphZRTNrC0w1s2ODuUAPy93XB9O14O7JwAKgtrsvcPdFWew/093XBW/nAYlmlpBDWGcBI4L1j4B+ZhZWC62IZC8uNoYruzdg3F29Ob9jHV6btII+j0/gs1/WqVuuSCCsBNTd3yN0J/aTYOnm7iOzOyabbkSVzGyMmS0JXitmOuaeoEvQIjM7JVN5x6B70VIzeyajolQ3oqKtVIk47jipKePu6s3pbWry4vhl9Hl8PO9O+5VUzb0lIoWcmf0LmA08AzwRLI/n4vj6QAdgWpiHnAfMdPcDmcreCLrf/i1TklkbWA3g7qnATqByFp8/yMySzCxp82ZNNSGSG5XLJPDweW0ZNbg7tSskcut7M7nuzRls2Lk/0qGJRFxuuuBWDV5jgePN7Nwc9j9cN6KhwFh3bwKMDd4TbLsIaEXoGdMXzCxjGLEXgUFAk2DpH5SrG1ExULN8SZ4c2J7Pbu5OgyqluffTOZz2zCTNvSUihd2FQCN37+3ufYKlb45HAWZWBvgYuM3dd4WxfytCdeD1mYovDbrm9giWyzJ2z+IU/9M04+6vuHsnd+9UtWrVLA4RkZy0rVOBT27qzl9Pa8GkpZs56ckJvDf9V7WGSrEW7ii4rwOvE7q7ekawnJ7dMYfrRsQfu/6MAM4O1s8CRrr7AXdfASwFuphZTaCcu0/x0P/WNw85Rt2Iiom2dSrwwfXdePHSY9mbksplw6Zz1Rt6yF9ECq25/HH8hLCYWTyh5PMdd/8kjP3rAJ8Cl7v7soxyd18bvCYD7wJdgk1rgLrBsXGExnXQiG8i+SQ2xri2R0NG39aT1rXLc88nc7jk1Wms3LIn0qGJRERcmPt1dfeWR/ohh3Qjqu7u6yGUpGaMykcoOZ2a6bA1QVlKsH5oecYxv3UjMrOMbkRbDvn8QYRaUKlXr96RXoYUAmbGqW1q0rdFNUZMXsmzY5dyytM/cOlx9bjtxKZUKl0i0iGKiGTIGD9hLvBbt1h3P/NwBwQ3UYcBC9z9yZw+wMwqAF8C97j7j5nK44AK7r4lSGhPB74LNn9GaMCiKcD5wPeu5hiRfHdM5dK8e91xvP/Tah78cgH9/zORO09qxtUnNCA2Ru0nUnyEm4BOMbOW7j4/tx9waDeibBooD9clKLuuQmF3IwJegdBk2jnFLIVfQlwsg3o24rxj6/D0d0t4Z9qvfDpzLbf0bcwVx9cnIU6TQItIxI0g1C12DhDug+vdCXWVnWNms4Kye4EE4FlCj8N8aWaz3P0U4GagMfA3M/tbsP/JwB5gdJB8xhJKPjNG4h0GvGVmSwm1fF50xFcoIrliZlzUpR69m1Xjr/+dy4NfLeCL2et45Py2NK9RLtLhiRSIcBPQEYSS0A2E7uIa4O7eNruDDtONaKOZ1QxaP2sCGcPC/9YlKFAHWBeU18miPPMxa9SNqHiqXCaBf53dmsu7HcO/v1rAv79ayFtTV3HPqS04tXUN1CNbRCJoi7s/k5sD3H0Sh58z9NMs9n8AeOAw+3c8zGfsBy7ITVwikrdqlE/k1cs78sXs9dz/2TxOf2YSN/VpzOA+jXQTXYq8cAchep3QHdn+/P785xnZHZBNN6KMrj8Er6MylV8UjGzbgNBgQ9OD7rrJZtY1OOflhxyTcS51IyrGmlQvyxtXdeHNq7tQKj6Om975mQtfnsIvq3dEOjQRKb5mmNlDZtYtYwqWnKZhEZHiw8w4o10txtzRizPa1eKZsUs4/ZlJ/Pzr9kiHJpKvLJx8zcy+D3fkvkzHnAD8wB+7Ht1L6DnQD4B6wK/ABe6+LTjmPuBqQiPo3ubuXwflnYDhQEnga+AWd3czSyQ06XYHgm5E7r48u7g6derkSUlJubkUiTKpael8kLSGJ8csYsvug5zToTZ3n9KMWhVKRjo0EYlSZjbD3Tvl8phxWRR7buvTwkB1p0j+G7dwE/d+OocNu/ZzdfcG3HlyU0qVCLezokjhc7i6M9wE9AVCI/l9zh8HUshxdL7CRpVo8ZG8P4UXxy/jtUkriDEY1KMh1/dqROkE/TIXkdw5kgS0KFHdKVIwkven8Og3i3hr6irqVirJw+e2pXvjKpEOS+SIHK7uDLcLbklCiefJhDkNi0iklU2M58/9m/P9nb04qWUNnvl+Kb0fH88HP60mLV09tUWk4KkLrohkp2xiPP86uzXvD+pKXEwMl742jb98NJud+1IiHZpIngmrBbQo0V3c4mvGqu088OV8Zv66gxY1y/G301pwvO4qikgY8qoF1Mxedffr8iKmgqS6U6Tg7U9J4+nvlvDqD8upXLoE/zq7Nae0qhHpsETCdkQtoMH8mTmdOMd9RAqDjsdU5JMbj+fZizuwa18Kl7w2jWtHJLF88+5IhyYixUQ0Jp8iEhmJ8bEMPbU5/72pO5XLJHD9WzMY/M7PbNl9IOeDRQqxbFtAzWw5cFd2xwP/dPdWeR1YftFdXIHQXcXXf1zBC+OWsT8ljT91PYYh/ZpQsXSJSIcmIoVQblpAc+pm6+4/501UBUd1p0hkpaSl88rE5fznuyWUKxnHo+e3pW/z6pEOSyRbRzQIkZm9Eca5d7r7bUcRW4FSJSqZbU4+wFPfLWbk9F8pmxjPrf2acFnXYygRF+7j0SJSHOQyAc1q9NsMGgVXRI7Yog3JDBk5k4Ubkrms6zHcO6AFJUto3lApnI5qFNyiRJWoZGXRhmQe+HI+PyzZQv3KpbhnQAtOblmd0NSzIlLcaRRc1Z0ihcX+lDQeH72I1yatoHG1Mjw9sD2ta5ePdFgi/+NoR8EVKdKa1SjLm1d34Y2rOhMXG8P1b83g4lenMnftzkiHJiJRyszizexWM/soWG42s/hIxyUi0S0xPpa/nt6St67pwq59KZzzwo+8PGEZ6RrhX6KEElCRgJnRp1k1vhnSg3+d3ZrFG3dzxnOTuOvDX9i4a3+kwxOR6PMi0BF4IVg6BmUiIketR5OqjL6tJ/2aV+ehrxdy6WvTWLdjX6TDEsmRuuCKHMau/Sk8//1S3vhxJbExxvW9GjKoZ0NKlYiLdGgiUsCOpAuumf3i7u1yKosGqjtFCi9358OkNdz/+TziYoyHzm3LaW1rRjoskaPrgmtmCWZ2iZnda2Z/z1jyPkyRwqNcYjz3DGjBd3f0om/zajz93RL6Pj6Bj2esUTcXEQlHmpk1ynhjZg2BtAjGIyJFkJlxYee6fHVrDxpULcPgd3/mzg9+IXl/SqRDE8lSuF1wRwFnAanAnkyLSJFXr3Ipnr/0WD66oRvVyyVw54e/cObzk5i2fGukQxORwu1uYJyZjTezCcD3wJ0RjklEiqj6VUrz0Q3duLVvYz6duYYBz/zAjFXbIx2WyP8Iqwuumc1199YFEE++UzciORrp6c5nv6zjkW8Wsn7nfvq3qsHQU5tTv0rpSIcmIvnoSEfBNbMEoBmhebMXuntUziCvulMkuiSt3MZt789i/c793NynMbf0bUxcrIZ+kYJ1tKPgTjazNnkck0jUiYkxzu5Qm+/v7M1dJzdl4pLNnPTUBB74Yj4796qri4j8zsxigVOA3kA/YLCZ3RHRoESkWOhUvxJfDenBWe1q8Z+xS7jg5Sms2qrOi1I4hJuAngDMMLNFZjbbzOaY2ez8DEykMCtZIpab+zZh/F29ObdDHYb9uILej49jxOSVpKSlRzo8ESkcPgeuBCoDZTMtIiL5rlxiPE8ObM+zF3dg6abdDPjPD3yYtJriNgCpFD7hdsE9Jqtyd1+V5xHlM3Ujkvwwf90uHvhyPpOXbaVh1dLcN6AFfZtXw8wiHZqI5IEjHAV3tru3za+YCpLqTpHotnbHPu54fxbTVmxjQJsa/PucNlQoVSLSYUkRd1RdcINEswJwRrBUiMbkUyS/tKxVjneuPY5hV4T+j10zIok/DZvG/HW7IhyZiETQ12Z2cqSDEBGpXaEk717Xlb/0b86Y+Rvp//QPTF66JdJhSTEV7jQsQ4B3gGrB8raZ3ZKfgYlEGzOjX4vqjL6tJ/ef0ZJ563Zx2rM/MPTj2WxK3h/p8ESk4E0FPjWzfWa2y8ySzUx3pUQkImJjjBt7N+LTm7pTKiGWS16bxsNfLyRVjw5JAQu3C+5soJu77wnelwamRGPXInUjkoKyc28Kz36/hBFTVhIfG8NNvRtxbY+GJMbHRjo0EcmlI+yCuxw4G5jjUf7QlepOkaJl38E0/vnFfN6b/ivHNajEs5d0oFrZxEiHJUXM0Y6Ca/xx8uy0oExEDqN8qXj+enpLvr29Fz2aVOHxbxfT9/HxjJq1lvT0qP5bVETCswSYm5vk08zqmtk4M1tgZvOCHkiY2QXB+3Qz63TIMfeY2dJgoMBTMpV3DAYNXGpmz1jwULqZJZjZ+0H5NDOrnzeXKyLRomSJWB46tw1PXtiOX9bs4LRnJjFV85tLAQk3AX0DmGZm95vZ/YS6FQ3Lt6hEipAGVUrz8mWdGDmoK5XKlGDIyFmc8+JkklZui3RoIpK/1gPjgwTxjowlh2NSgTvdvQXQldDULS2BucC5wMTMOwfbLgJaAf2BF4LpXwBeBAYBTYKlf1B+DbDd3RsDTwGPHOV1ikiUOvfYOowafAJlE+K45NWpvDh+mW6SS74LdxCiJ4GrgG3AduAqd386H+MSKXK6NqzMZ4NP4PEL2rFh5z7Of2kKg9/5mdXb9kY6NBHJHyuAsUAJwpyGxd3Xu/vPwXoysACo7e4L3H1RFoecBYx09wPuvgJYCnQxs5pAOXefErTAvkmoO3DGMSOC9Y+AfhmtoyJS/DSrUZbPbjmBU9vU5JFvFjLorRma21zyVVx2G82snLvvMrNKwMpgydhWyd3VhCOSCzExxvkd6zCgTQ1embiclycsZ8z8jVx1Qn0G92lMucT4SIcoInnE3f/vaI4PusZ2AKZls1ttQr2SMqwJylKC9UPLM45ZHcSYamY7Cc1V+ochMc1sEKEWVOrVq3eklyEiUaBMQhzPXdyBzsdU5MGvFnD6cz/w4qUdaV27fKRDkyIopxbQd4PXGUBSpiXjvYgcgVIl4rjtxKaMu6s3Z7avxSsTl9P7sfG8NXWVRqMTiXLBoypHtY+ZlQE+Bm5z9+xGzs2q5dKzKc/umD8WuL/i7p3cvVPVqlWzC1dEigAz48ruDXj/+m6kpjnnvjiZ96b/SpSPoSaFULYtoO5+evDaoGDCESleapRP5PEL2nHl8fX51xfz+dt/5/Lm5JXcd1oLejerFunwROTIXJvDdCtG6LnN+7PcaBZPKPl8x90/yeGz1gB1M72vA6wLyutkUZ75mDVmFgeUJ/SIjYgIx9aryJe39mDIyJnc88kcflq5jQfPbkPJEhrFX/JGuPOAjg2n7JDtr5vZJjObm6nsfTObFSwrzWxWUF4/mCctY9tLmY7RKH5S5LWuXZ6Rg7ry8mUdSUlL58o3fuLy16ezaENypEMTkdx7lT8+83noUibY538EddwwYEEw/kJOPgMuCurEBoQGG5ru7uuBZDPrGpzzcmBUpmOuCNbPB76P9mliRCRvVSpdguFXdeG2E5vw6cy1nPPCjyzfvDvSYUkRkdMzoIlAKaCKmVXk92475YBaOZx7OPAcoYEPAHD3gZnO/QSwM9P+y9y9fRbnyRjFbyrwFaFR/L4m0yh+ZnYRoVH8BmZxvEhUMDNOaVWDPs2q8eaUlTwzdgmn/mciF3Wpxx0nNaVKmYRIhygiYTjKZz+7A5cBczJu0gL3AgnAs0BV4Eszm+Xup7j7PDP7AJhPaATdwe6eMW3ajYTq4pKE6s2vg/JhwFtmtpRQy+dFRxGviBRRsTHGbSc2pUO9itw2ciZnPvcjj57flgFtakY6NIlylt1Nz2D+sdsIJZtr+T0B3QW86u7PZXvyUKvkF+7e+pByA34F+rr7kmz2qwmMc/fmwfuLgd7ufr2ZjQbud/cpQReiDUDVnO7iajJtiRbb9xzkP2OX8PbUVSTGxzK4T2Ou6l6fxHh1gREpaIebTLu4UN0pUryt27GPwe/+zMxfd3B19wYMPbU5JeLCnc1RiqvD1Z3Z/stx9/8Ez3/e5e4N3b1BsLTLKfnMQQ9go7svyVTWwMxmmtkEM+sRlNUmzFH8CLWmVs7qw8xskJklmVnS5s2bjyJskYJTsXQJ7j+zFaNv70nXhpV55JuF9HtiAp//sk4DAoiIiEiBqVWhJO8P6sZV3evz+o8ruOiVKazfuS/SYUmUCnce0GfNrLWZXWhml2csR/G5FwPvZXq/Hqjn7h2AO4B3zawceTCKXxC/RvKTqNWoahleu6IT71x7HOVKxnPLezM578XJzPx1e6RDExERkWKiRFwM/zijFc9d0oFFG5I57ZlJ/LBEDTuSe+EOQvQPQs+ePAv0AR4FzjySDwy6y54LvJ9RFkygvTVYnwEsA5oS3ih+GefUKH5SpHVvXIUvbjmBR89ry+rt+zjnhcnc+t5M1mzfG+nQRCQLZlbVzO41s1eCgfleN7PXIx2XiMjROL1tLT675QSqlCnB5a9P5z/fLSE9XT2zJHzhdt4+H+gHbHD3q4B2hAZEOBInAgvd/beutUElHRusNyQ0it9yjeIn8kexMcaFnesy7q7e3NK3MaPnbaDfExN4bPRCdh9IjXR4IvJHowjdHP0O+DLTIiIS1RpVLcN/B3fn7Pa1eeq7xVw5/Ce27TkY6bAkSoSbgO5z93QgNegauwlomN0BZvYeMAVoZmZrzOyaYNNF/LH7LUBPYLaZ/QJ8BNzg7hmtmTcCrwFLCbWMZh7Fr3Iwit8dwNAwr0Uk6pVJiOPOk5sx7q7eDGhTk+fHLaP3Y+N4b/qvpOkupEhhUcrd/+LuH7j7xxlLpIMSEckLpUrE8eSF7fj3OW2YumwrZz43iYUbspsCWSQk21Fwf9vJ7AVCw8BfBNwJ7AZmBa2hUUUj+UlRNGv1Dh74Yj5Jq7bTvEZZ7jutBT2a6HlnkbxyJKPgmtkDwGR3/yqfwiowqjtFJDu/rN7BdW8msedAKk8NbM/JrWpEOiQpBA5Xd+aYgAZdX+u4++rgfX2gnLvPzo9A85sqUSmq3J2v527goa8XsHrbPvo0q8p9p7WgcbWykQ5NJOodYQKaDJQGDgIpQbG7e7m8ji+/qe4UkZxs3LWfQW8mMXvtTu46uRk39W5EKI2Q4uqIpmGBUE0J/DfT+5XRmnyKFGVmxoA2Nfnujl7cO6A5SSu3c8rTP/D3UXP1XIZIBLh7WXePcffEYL1sNCafIiLhqF4ukfev78YZbWvx2OhF3Pb+LPanpEU6LCmEwn0GdKqZdc7XSEQkTyTExTKoZyPG392bS7rU451pv9LrsXG8OnE5B1JVEYgUJDM708weD5bTIx2PiEh+SoyP5T8XtefuU5oxatY6Br48hY279kc6LClkwk1A+wBTzGyZmc02szlmplZQkUKscpkE/nV2a74Z0oNOx1Tkwa8WcNKTE/l6zno0YLRI/jOzh4EhwPxgGRKUiYgUWWbG4D6NeeWyjizZtJszn5vE7DU7Ih2WFCLhDkJ0TFbl7r4qzyPKZ3qORYqriYs388CX81m8cTdd6lfir6e3oG2dCpEOSyQqHOEzoLOB9sEo8gTTjc1097b5EWN+Ut0pIkdiwfpdXDsiiS27D/DYBe04s12tSIckBeiInwENPODuqzIvwAN5G6KI5KeeTavy1a09+Pc5bVi+ZTdnPvcjd7w/i/U790U6NJGirEKm9fKRCkJEJBJa1CzHZzd3p12dCtz63kweH72IdE0XV+yFm4C2yvwmuIvbMe/DEZH8FBcbwyXH1WPcXb25qXcjvpiznj6Pj+fJbxex50BqpMMTKWoeAmaa2XAzGwHMAP4d4ZhERApU5TIJvH3tcVzUuS7PjVvKDW/P0N8cxVy2CaiZ3RMMI9/WzHYFSzKwCRhVIBGKSJ4rmxjPn/s3Z+wdvTipZQ2e+X4pfR4fzwdJq0nTnUmRPOHu7wFdgU+CpZu7j4xsVCIiBa9EXAwPnduGf5zRku8WbOS8FyezetveSIclEZJtAuruD7l7WeAxdy8XLGXdvbK731NAMYpIPqlbqRTPXtyBj288ntoVS/Lnj2ZzxrOTmLxsS6RDE4laZtY8eD0WqAmsAVYDtYIyEZFix8y4qnsDhl/VhbU79nHW8z8yfcW2SIclERDWIEQAZlYbOAaIyyhz94n5FFe+0UAKIllzdz6fvZ5Hvl7I2h37OLFFde4d0JyGVctEOjSRiMvNIERm9oq7DzKzcVlsdnfvm8fh5TvVnSKSl5Zt3s11I5JYvX0vD5zdmoGd60U6JMkHh6s7wx0F92HgIkLDyGdMJOjufmaeRlkAVImKZG9/Shqv/7iCF8YtY39KGpd1O4Yh/ZpQoVSJSIcmEjFHOApuorvvz6ksGqjuFJG8tnNvCje/9zM/LNnCVd3rc9+AFsTFhjs8jUSDox0F9xygmbsPcPczgiXqkk8RyVlifCw39W7MuLt6c2HnuoyYvJJej41n2KQVHExNj3R4ItFkcphlIiLFTvlS8bxxZWeu6l6fN35cyVXDf2LnvpRIhyUFINwEdDkQn5+BiEjhUrVsAv8+pw1fDelB2zrl+dcX8znl6Yl8O28D4XbdFymOzKyGmXUESppZBzM7Nlh6A6UiG52ISOERFxvDP85oxcPntmHq8q2c8/yPLN+8O9JhST6Ly3kXAPYCs8xsLHAgo9Ddb82XqESk0GheoxxvXt2F8Ys38+CXCxj01gy6NqzEX09rSevamtZQJAunAFcCdYAnM5UnA/dGIiARkcLsoi71aFClNDe+8zNnP/8jz196LD2aVI10WJJPwn0G9Iqsyt19RJ5HlM/0HIvIkUtNS+e96b/y1HdL2L73IOcdW4e7T2lG9XKJkQ5NJF8d4TOg57n7x/kVU0FS3SkiBWH1tr1c92YSSzft5qFz23BBp7qRDkmOwlENQhScoCRQz90X5XVwBUmVqMjR27kvhRfGLeWNH1cSG2Pc0KsR1/VsQKkS4XaqEIkuR5KABsedBrQCfrtL4+7/zGb/usCbQA0gHXjF3f9jZpWA94H6wErgQnffbmaXAndnOkVb4Fh3n2Vm4wlNA7Mv2Hayu28ys4TgMzoCW4GB7r4yu+tQ3SkiBSV5fwo3vv0zk5Zu4Y6TmnJL38aYWaTDkiNwVIMQmdkZwCzgm+B9ezP7LE8jFJGoUb5kPPcMaMF3d/SiT/OqPPXdYvo+PoGPZ6whPV3Ph4oAmNlLwEDgFsCACwhNZ5adVOBOd28BdAUGm1lLYCgw1t2bAGOD97j7O+7e3t3bA5cBK919VqbzXZqx3d03BWXXANvdvTHwFPDI0V+tiEjeKJsYz+tXdubcDrV5csxi7v10DqlpGgSxKAl3EKL7gS7ADoCgcmuQLxGJSNSoV7kUL1zakQ9v6Ea1cgnc+eEvnPX8j0xbvjXSoYkUBse7++WEkr3/A7oB2fYnc/f17v5zsJ4MLABqA2cBGY+9jADOzuLwi4H3wogr87k+AvqZmhdEpBApERfDExe246bejXhv+moGvTWDvQdTIx2W5JFwE9BUd995SJmaOUQEgM71K/Hfm7rz9MD2bNl9gIGvTOWGt2awcsueSIcmEkkZ833uNbNaQAq5uHlrZvWBDsA0oLq7r4dQkgpUy+KQgfxvAvqGmc0ys79lSjJrA6uDc6UCO4HKWXz+IDNLMrOkzZs3hxu2iEieMDP+3L85/zq7NeMXbeLiV6ayZfeBnA+UQi/cBHSumV0CxJpZEzN7Fs1lJiKZxMQYZ3eozfd39ubOk5oycclmTnpqAg98MV/zeklx9bmZVQAeA34m9OxmOC2UmFkZ4GPgNnffFcb+xwF73X1upuJL3b0N0CNYLsvYPYtT/M9NZXd/xd07uXunqlU1GqWIRMZlXY/hpT91ZNHGZM57cbJubhcB4SagtxAaROEA8C6hu6W35VNMIhLFSpaI5ZZ+TRh/V2/O7VCHYT+uoPdj4xgxeSUpeoZDigkziyH0zOaOYCTcY4Dm7v73MI6NJ5R8vuPunwTFG82sZrC9JrDpkMMu4pDk1t3XBq/JhOruLsGmNQRdgc0sDigPbMv1RYqIFJCTW9Xg3eu6krw/lXNfnMzMX7dHOiQ5CmEloO6+193vc/fOwfJXd9+f85EiUlxVK5fII+e35ctbetCiZjn+8dk8+j89ke8XbiTc0bdFopW7pwNPZHp/IItHWf5H0E12GLDA3TPPIfoZkDEl2hXAqEzHxBAa4GhkprI4M6sSrMcDpwNzszjX+cD3rv+UIlLIHVuvIh/feDxlEuK4+NWpjJm/MdIhyREKdxTcMUE3ooz3Fc1sdL5FJSJFRsta5Xjn2uN47fJOuMPVw5O4bNh0FqzPsVehSLT71szOy+UAP90JdZXtGzy7OcvMBgAPAyeZ2RLgpOB9hp7AGndfnqksARhtZrMJjWK/Fng12DYMqGxmS4E7CEbUFREp7BpUKc3HNx5P0+pluf6tJN6euirSIckRCGseUDOb6e4dcio7ZPvrhO64bnL31kHZ/cB1QMZoBve6+1fBtnsIDQ2fBtzq7qOD8o7AcKAk8BUwxN39SOYxA81lJhJJKWnpvDN1FU+PXcKufSlc2Kkud5zclGplE3M+WCSCjmQeUDNLBkoTmlplP6FnL93dy+VDiPlKdaeIFCZ7D6Yy+J2fGbdoM4P7NOKuk5tprtBC6KjmAQXSzaxeppMdQ86j4A4H+mdR/lSmOckyks+WhJ5faRUc84KZxQb7vwgMApoES8Y5NY+ZSJSJj43hyu4NmHBXH67q3oCPf15Dn8fG8/y4pexPSYt0eCJ5yt3LunuMu5dw93LB+6hLPkVECptSJeJ49fJOXNS5Ls+PW8adH/7CwVSNMxEtwk1A7wMmmdlbZvYWMBG4J7sD3H0i4Q9qcBYwMnhGZgWwFOgSDLRQzt2nBM+nvMnvc59pHjORKFW+VDx/O70l397eixOaVOGx0Yvo+/h4Rs1aq+dDpcgws7HhlImISO7Fxcbw0LltuOOkpnzy81quGfETyfs16n40CHcQom+AY4H3gQ+AjhldZI/AzWY228xeN7OKQdlvc5IF1gRltYP1Q8v/cEx285iB5jITKawaVCnNy5d1YuSgrlQqU4IhI2dxzguTmbFKA3JK9DKzRDOrBFQJxkyoFCz1gVoRDk9EpMgwM27t14RHz2/L5GVbGfjyVDbu0jiphV24LaAQGtBgG6FEr6WZ9TyCz3sRaAS0B9bz+wiBh5uTLLu5ysKaxww0l5lIYde1YWU+G3wCj1/QjvU793Hei1MY/O7PrN62N9KhiRyJ64EZQPPgNWMZBTwfwbhERIqkCzvV5fUrO7Ny6x7OfWEySzYmRzokyUa4o+A+AvxIqCvu3cFyV24/zN03untaMDz9q2QxJ1mgDrAuKK+TRfkfjtE8ZiLRLybGOL9jHcbd1ZvbTmzC9ws20e+JCTz09QJ2qUuNRBF3/4+7NwDucveG7t4gWNq5+3ORjk9EpCjq1bQqH1zfjQOp6Zz34mSmr1BaUFiF2wJ6NtDM3U9z9zOC5czcfljGJNqBc/jjnGQXmVmCmTUgNNjQdHdfDySbWdfg+c7L+X3uM81jJlIElSoRx20nNmXcXb05o10tXp6wnD6PjeftqatITdMAAxI93P1ZMzvezC4xs8szlkjHJSJSVLWuXZ5PbzqeKmUT+NOwaXw1Z32kQ5IshJuALgfic3NiM3sPmAI0M7M1ZnYN8KiZzQnmJesD3A7g7vMIPVs6H/gGGOzuGUNi3gi8RmhgomXA10G55jETKcJqlE/kiQvb8fnNJ9CoWhn++t+5nPqfHxi/aFOkQxMJSzBo3+PACUDnYMnVVC4iIpI7dSuV4uMbjqdN7fIMfvdnXp+0ItIhySHCnQf0Y6AdMBY4kFHu7rfmX2j5Q3OZiUQfd2f0vI089PUCVm3dS6+mVbnvtBY0rV420qFJMXGE84AuAFoWhd45qjtFJNrsT0ljyMiZjJ63kdtObMKQfk00V2gBO1zdGRfm8Z8Fi4hIgTMz+reuQd/m1XhzykqeGbuE/k9P5OIu9bj9pKZUKZMQ6RBFsjIXqEFo0D0RESlAifGxPH/JsQz9ZA5Pf7eE3ftTue+0FkpCC4GwElB3H2FmJYCmQdEid9eoICJSoErExXBtj4acd2wd/jN2CW9PXcWoWesY3KcxV3WvT2J8bKRDFMmsCjDfzKbzx95DuR5DQUREci8uNoZHz2tLmYQ4Xpu0gt0HUnnwnDbExigJjaSwElAz6w2MAFYSmv6krpld4e4T8y0yEZHDqFi6BPef2YrLuh3DQ18t4JFvFvLOtFUMPbU5p7WpqbubUljcH+kARESKu5gY4x9ntKRsYhzPfr+U3QdSeWpge+JjczMbpeSlcH/yTwAnu3svd+8JnAI8lX9hiYjkrFHVMrx2RWfeufY4yiTEcfO7MznvxcnM/HV7pEMTwd0nELpxGx+s/wT8HNGgRESKITPjzpObMfTU5nwxez03vDWD/SlpOR8o+SLcBDTe3RdlvHH3xeRyVFwRkfzSvXEVvry1B4+c14bV2/dxzguTGTJyJmt37It0aFKMmdl1wEfAy0FRbeC/EQtIRKSYu6FXIx44uzXfL9rEVW/8xO4DqZEOqVgKNwFNMrNhZtY7WF4FZuRnYCIiuREbYwzsXI9xd/Xmlr6N+WbuBvo+Pp7HRi9UBSORMhjoDuwCcPclQLWIRiQiUsz9qesxPHVhe6av3MafXpvGjr0HIx1SsRNuAnojMA+4FRhCaL7OG/IrKBGRI1UmIY47T27GuLt6M6BNTZ4ft4zej41n5PRfSUuP+tkwJLoccPff/rIxszhA/whFRCLs7A61efHSY5m/bhcXvTKVTcn7Ix1SsRJuAhoH/Mfdz3X3c4BnAA03KSKFVq0KJXlqYHv+O7g79SuXYugnczjtmR+YtGRLpEOT4mOCmd0LlDSzk4APgc8jHJOIiAAnt6rB61d2ZtXWvQx8eaoe2ylA4SagY4GSmd6XBL7L+3BERPJW+7oV+PCGbrxw6bHsOZjKn4ZN4+rhP7F0U3KkQ5OibyiwGZgDXA98Bfw1ohGJiMhvTmhShbev7cKW3Qe44MXJLN+8O9IhFQvhJqCJ7v7bNxKsl8qfkERE8paZMaBNTcbc3ot7Tm3OTyu2ccrTP/CPUXPZtkfPfki+KQm87u4XuPv5wOv88WauiIhEWMdjKjFyUFcOpKZz4ctTWLB+V6RDKvLCTUD3mNmxGW/MrCOgdmoRiSqJ8bFc36sR4+/uzSVd6vH2tF/p9dg4Xp24nAOpGo5d8px6D4mIRIFWtcrz/vXdiI+NYeDLUzSdWz4LNwG9DfjQzH4wsx+A94Gb8y0qEZF8VLlMAv86uzXfDOlBp2Mq8uBXCzj5qYl8M3c97hojRvKMeg+JiESJxtXK8MH13ahYugSXvjaNycs0ZkR+CSsBdfefgOaERsO9CWjh7pqGRUSiWpPqZXnjqi68eXUXEuJiuOHtnxn48lRmr9kR6dCkaFDvIRGRKFK3Uik+vL4bdSqW5Mo3fmLsgo2RDqlICrcFFHdPcfe57j7H3VPyMygRkYLUs2lVvrq1B/8+pw3Lt+zmzOd+5I73Z7F+p3IFOSq3kcveQ2ZW18zGmdkCM5tnZkOC8kpmNsbMlgSvFYPy+ma2z8xmBctLmc7V0czmmNlSM3vGzCwoTzCz94PyaWZWP5+uX0Qk6lQrl8j7g7rRvEZZrn9rBp//si7SIRU5YSegIiJFWVxsDJccV49xd/Xmxt6N+GLOevo8Pp4nxyxmz4HUSIcnUegIew+lAne6ewugKzDYzFoSGlF3rLs3IfRs6dBMxyxz9/bBknmO7heBQUCTYOkflF8DbHf3xsBTwCNHc50iIkVNxdIleOfa4zj2mIrcOnImI6f/GumQihQloCIimZRNjOcv/Zsz9o5enNSyBs+MXUKfx8fzQdJq0tL1fKjkWmegLdABuNjMLs9uZ3df7+4/B+vJwAKgNnAWMCLYbQRwdnbnMbOaQDl3n+KhB5vfzHRM5nN9BPTLaB0VEZGQsonxjLiqCz2bVGXoJ3N47YflkQ6pyAgrAbWQP5nZ34P39cysS/6GJiISOXUrleLZizvw8Y3HU6tCSf780WzOeHaSBiWQsJnZW8DjwAmEEtHOQKdcHF+fUOI6Daju7ushlKQC1TLt2sDMZprZBDPrEZTVBtZk2mdNUJaxbXVwrlRgJ1A5VxcnIlIMlCwRy6uXd+K0NjV54MsFPP3dYg1WmAfiwtzvBSAd6Av8E0gGPiZUmYqIFFkdj6nIpzcdz+ez1/PI1wu55NVpnNSyOvec2pyGVctEOjwp3DoBLf0I/loxszKE6tnb3H1XNg2U64F67r41GOTov2bWCsjqgIw4stuWOYZBhLrwUq9evVxegYhI0VAiLoZnLu5AqRKxPP3dEtId7jipaaTDimrhdsE9zt0HA/sB3H07UCLfohIRKUTMjDPb1WLsnb34c/9mTFm2lZOfmsg/P5/Pjr0HIx2eFF5zgRq5PcjM4gkln++4+ydB8cagW21G99pNAO5+wN23BuszgGVAU0ItnnUynbYOkDGSxhqgbnCuOKA8sO3QONz9FXfv5O6dqlatmtvLEBEpMmJjjEfOa8vATnV5ZuwSnhqzONIhRbVwE9AUM4sluENqZlUJtYiKiBQbifGx3NS7MePu6s0FneoyfPIKej02ntcnreBgqn4lyv+oAsw3s9Fm9lnGkt0BwbOYw4AF7v5kpk2fAVcE61cAo4L9qwb1M2bWkNBgQ8uDbrrJZtY1OOflGccccq7zge+PpJVWRKQ4iYkxHjq3DRd2qsN/lIQelXC74D4DfApUM7MHCVVYf823qERECrGqZRN46Nw2XHH8MTz45QL++cV83pq6intObc5JLauj8VwkcP8RHNMduAyYY2azgrJ7gYeBD8zsGuBX4IJgW0/gn2aWCqQBN7h7RmvmjcBwoCTwdbBAKMF9y8yWEmr5vOgI4hQRKXZiYoyHz22LO/xn7BLM4LYT1R03tyzcm55m1hzoR+jZkbHuviA/A8svnTp18qSkpEiHISJFhLszftFmHvhyPss276Fbw8rcd1oLWtcuH+nQJA+Z2Qx3D3sAoUzHVef38RKmu/umvI2sYKjuFBH5XXq685ePZ/PhjDXcfmJThpzYJNIhFUqHqzvDHQW3K7DW3Z939+eANWZ2XF4HKSISbcyMPs2r8c1tPfnXWa1YtDGZM56bxN0f/sLGXfsjHZ5EkJldCEwn1Fp5ITDNzM6PbFQiInK0YoJnQs/vWIenvlvMM2OXRDqkqBJuF9wXgWMzvd+TRZmISLEVHxvDZd3qc2b72rwwbilv/LiSL+es54ZejbiuR0NKloiNdIhS8O4DOme0egbjJ3xHaO5NERGJYhlJqDs8OWYxBtzSTy2h4Qh3ECLLPECBu6eTQ/JqZq+b2SYzm5up7DEzW2hms83sUzOrEJTXN7N9ZjYrWF7KdExHM5tjZkvN7JmMybLNLMHM3g/KpwXzpYmIRFT5kvHcM6AFY+7oSe9mVXlyzGL6PD6eT35eQ3q6xnkpZmIO6XK7lfDrXRERKeRiY4xHz2/LucfW5okxi3lWLaFhCbciXG5mt5pZfLAMAZbncMxwoP8hZWOA1u7eFlgM3JNp2zJ3bx8sN2Qqf5HQPGRNgiXjnNcA2929MfAU8EiY1yIiku+OqVyaFy7tyIc3dKNauQTu+OAXzn7hR6av+J/ZLqTo+iYYAfdKM7sS+JLfBwISEZEiIDbGeOz8dpzbIZSEPve9ktCchJuA3gAcD6wlNH/YcQSTUx+Ou0/kkHnF3P1bd08N3k7lj3OU/Y9grrNy7j4laIF9Ezg72HwWMCJY/wjol9E6KiJSWHSuX4n/3tSdpwa2Y3PyAS58eQo3vj2DVVv3RDo0yWfufjfwMtAWaAe84u5/jmxUIiKS12JjjMcuCCWhj3+7mOfHLY10SIVaWM+ABl2I8nqY9quB9zO9b2BmM4FdwF/d/QegNqGEN8OaoIzgdXUQX6qZ7QQqA1sO/SAzG0SQMNerVy+PL0NEJHsxMcY5HerQv1VNXvthOS9OWMZ3CzZy5fH1ublvE8qXjI90iJKHzKwxUN3df3T3T4BPgvKeZtbI3ZdFNkIREclrGUmoA4+NXgTA4D6NIxtUIRVWAmpmiYS6vLYCEjPK3f3qI/lQM7sPSAXeCYrWA/XcfauZdQT+a2atCE35cqiMh6iy2/bHQvdXgFcgNJT8kcQsInK0SpaI5ZZ+TRjYuS5PfLuY1yat4KMZa7j9pKZc3KUe8bF6PLCIeJrQ3J2H2htsO6MggxERkYIRG2M8fkE73J3HRi/CDG7qrST0UOH+tfMWUAM4BZhAqOts8pF8oJldAZwOXJoxsJG7H3D3rcH6DGAZ0JRQi2fmbrp1gHXB+hqgbnDOOKA8h3T5FREpjKqVS+SR89vyxS0n0LxGOf4+ah79n57I9ws3Eu7czFKo1Xf32YcWunsSUL/gwxERkYISG2M8cWF7zmpfi0e/WcSL49Xp5VDhJqCN3f1vwB53HwGcBrTJ7YeZWX/gL8CZ7r43U3lVM4sN1hsSGmxoubuvB5LNrGvwfOflwKjgsM+AK4L184HvXX+5iUgUaVWrPO9edxyvXd4Jd7h6eBKXDZvOwg27Ih2aHJ3EbLaVLLAoREQkImJjjCcuaMeZ7WrxyDcLlYQeItx5QFOC1x1m1hrYQA53cc3sPaA3UMXM1gD/IDTqbQIwJhgvaGow4m1P4J9mlgqkATe4e0Zr5o2ERtQtSWj0wIwRBIcBb5nZUkItn3n9jKqISL4zM05sWZ1ezary9tRVPP3dEgb85wcGdq7L7Sc1pVrZ7HIZKaR+MrPr3P3VzIVmdg0wI0IxiYhIAYqLjeHJC0PPhD7yzULM4IZejSIdVqFg4TQamtm1wMeEWj2HA2WAv7n7y/kaXT7o1KmTJyUlRToMEZEs7dh7kGe/X8qbU1ZSIjaGm/o05poTGpAYHxvp0Io1M5vh7p3C3Lc68ClwkN8Tzk5ACeAcd9+QP1HmH9WdIiJHJjUtnds/+IXPf1nHPac25/pilIQeru7MtgXUzIa4+3+ABe6+HZgINMynGEVEir0KpUrwt9Nb8qeux/Dw1wt4bPQi3p32K3/u34wz29VCs00Vfu6+ETjezPoArYPiL939+wiGJSIiERAXG8NTF4YGJnro61BL6KCexScJzUpOz4BeFbw+m9+BiIjI7xpUKc3Ll3Xiveu6UqFUPENGzuKcFyYzY5XGWosW7j7O3Z8NFiWfIiLFVFxsDE8PbM9pbWvy768W8urE5ZEOKaJyegZ0gZmtBKqZWeYR/Qxwd2+bb5GJiAjdGlXm85tP4JOZa3ls9ELOe3EKp7WtydD+zalbqVSkwxMREZEwxMXG8J+B7cHhwa8WAHBdz+LZsTTbBNTdLzazGsBo4MyCCUlERDKLiTHO71iHAW1q8PKE5bw8cRlj5m/k6u4NuKlPI8olxkc6RBEREclBXGwMT1/UHggloYnxMVzWrX5EY4qEcEbB3QzMcfdV+R2MiIgcXqkScdx+UlMu7lKPx0Yv4qUJy/gwaTW3n9SUizrXJS423Jm1REREJBLigyR0f0oafxs1jzKJcZzToU6kwypQOf614u5phKZSKVEA8YiISA5qlE/kiQvb8fnNJ9CoWhn++t+5DHjmByYs3hzp0ERERCQH8bExPH/psXRrWJm7PpzN6HlRNzj6UQn3dvkq4Ecz+5uZ3ZGx5GdgIiKSvTZ1yvP+oK689KeOHEhN54rXp3PF69NZvDE50qGJiIhINhLjY3n1ik60rl2eW96dyY9Lt0Q6pAITbgK6Dvgi2L9spkVERCLIzOjfugZjbu/FX09rwcxft9P/6Ync9+kctuw+EOnwRERE5DDKJMQx4qrONKhSmuveTGLGqu2RDqlAmLtHOoYCpcm0RaQo277nIP8Zu4S3pq6iVHwsg/s25srj65MYHxvp0KLa4SbTLi5Ud4qI5J9Nu/ZzwctT2L7nICMHdaNlrXKRDilPHK7uDKsF1MzGmdn3hy55H6aIiByNiqVLcP+ZrRh9W0+Oa1iJh79eyElPTeDL2espbjccRUREokG1com8fc1xlE6I4/LXp7F88+5Ih5Svwu2Cexdwd7D8DZgF6FaoiEgh1bhaGV67ojPvXHscpUvEMfjdnzn/pSnMWr0j0qGJiIjIIepWKsVb1xyHO/zptWms3bEv0iHlm7ASUHefkWn50d3vAI7L59hEROQodW9chS9v7cEj57Vh1da9nP38jwwZObNIV2zRzMzqBr2OFpjZPDMbEpRXMrMxZrYkeK0YlJ9kZjPMbE7w2jfTucab2SIzmxUs1YLyBDN738yWmtk0M6sfkYsVEZE/aFytDCOu7kLygVT+9No0NicXzbEcwu2CWynTUsXMTgFq5HNsIiKSB2JjjIGd6zH+7t7c0rcx38zdQN/Hx/PY6IXsPpAa6fDkj1KBO929BdAVGGxmLYGhwFh3bwKMDd4DbAHOcPc2wBXAW4ec71J3bx8sm4Kya4Dt7t4YeAp4JH8vSUREwtW6dnneuLIzG3bu57Jh09i5NyXSIeW5cLvgziDU5XYGMAW4k1AFJiIiUaJMQhx3ntyMcXf15tTWNXh+3DJ6PzaekdN/JS1dz4cWBu6+3t1/DtaTgQVAbeAsYESw2wjg7GCfme6+LiifBySaWUIOH5P5XB8B/czM8uwiRETkqHSqX4mXL+vI8s17uHL4dPYUsZvF4XbBbeDuDYPXJu5+srtPyu/gREQk79WqUJKnL+rAfwd3p37lUgz9ZA6nPfMDk5YUnznIokHQNbYDMA2o7u7rIZSkAtWyOOQ8YKa7Z+6z9UbQ/fZvmZLM2sDq4FypwE6gchafP8jMkswsafPmzXl1WSIiEoaeTavyzMXt+WX1Dga9lcT+lLRIh5Rnsk1AzayzmdXI9P5yMxtlZs+YWaX8D09ERPJL+7oV+PCGbrxw6bHsOZjKn4ZN45rhP7F0U9EefS8amFkZ4GPgNnffFcb+rQh1pb0+U/GlQdfcHsFyWcbuWZzif5rA3f0Vd+/k7p2qVq2a20sQEZGj1L91TR49vx0/Lt3KLe/NJCUtPdIh5YmcWkBfBg4CmFlP4GHgTUJ3S1/J39BERCS/mRkD2tRkzO29uOfU5kxfsY1Tnp7IP0bNZfueg5EOr1gys3hCyec77v5JULzRzGoG22sCmzLtXwf4FLjc3ZdllLv72uA1GXgX6BJsWgPUDY6NA8oD2/LzmkRE5Mic37EO/3dmK8bM38ifP5pNehF4ZCanBDTW3TMqpYHAK+7+sbv/DWicv6GJiEhBSYyP5fpejRh/d28u6VKPt6f9Sq/HxvHaD8s5kFp0uv0UdkE32WHAAnd/MtOmzwgNMkTwOirYvwLwJXCPu/+Y6TxxZlYlWI8HTgfmZnGu84HvXZPEiogUWlccX5+7Tm7KpzPX8vfP5kb9vN45JqDB3VGAfsD3mbbFZbG/iIhEscplEvjX2a35ZkgPjj2mIg98uYCTn5rIN3PXR32FFyW6E+oq2zfT9CkDCPVAOsnMlgAnBe8BbiZ0Q/hvh0y3kgCMNrPZhObuXgu8GhwzDKhsZkuBO/h9RF0RESmkBvdpzPU9G/L21F95dPSiSIdzVHJKIt8DJpjZFmAf8AOAmTUm1A1XRESKoCbVyzL8qi5MWLyZB7+czw1v/0yXBpX422ktaVOnfKTDK7KCAf4ONyJtvyz2fwB44DD7dzzMZ+wHLjiiAEVEJCLMjKGnNif5QCovjl9G2cQ4buodnR1Ss01A3f1BMxsL1AS+zdRFJwa4Jb+DExGRyOrVtCrdG/Xgg6Q1PDlmEWc8N4lzj63N3ac0o2b5kpEOT0REpNgwM/51Vmv2HEjl0W8WUTYxnsu6HhPpsHItx2607j41i7LF+ROOiIgUNnGxMVxyXD3OaFeTF8YvY9ikFXw1Zz2Dejbihl4NKVVCT2SIiIgUhNgY4/EL2rHnQCp/HzWXMgmxnNOhTqTDypWw5gEVEREpmxjPX/o3Z+wdvTipZQ2eGbuE3o+N58Ok1UViVD4REZFoEB8bw3OXHEvXBpW568PZfDtvQ6RDypV8S0DN7HUz22RmczOVVTKzMWa2JHitmGnbPWa21MwWmdkpmco7mtmcYNszGRNpm1mCmb0flE8LJuwWEZF8VrdSKZ69uAMf33g8tSqU5O6PZnPGc5OYsmxrpEMTEREpFhLjY3n1ik60rl2em9+dydTl0VMH52cL6HCg/yFlQ4Gx7t4EGBu8x8xaAhcBrYJjXjCz2OCYF4FBQJNgyTjnNcB2d28MPEVoAm4RESkgHY+pyKc3Hc8zF3dgx94ULn51KoPeTGLFlj2RDk1ERKTIK5MQx4irOlOvcimuezOJRRuSIx1SWPItAXX3ifzvxNZnASOC9RHA2ZnKR7r7AXdfASwFugSTbZdz9ynBAEhvHnJMxrk+AvpltI6KiEjBMDPObFeLsXf24s/9mzF52VZOenIC//x8Pjv2Hox0eCIiIkVahVIlGH5VZ0rGx3LF69NZt2NfpEPKUUE/A1rd3dcDBK/VgvLawOpM+60JymoH64eW/+EYd08lNC1M5XyLXEREDisxPpabejdm3F29uaBTXYZPXkGvx8bz+qQVpKSlRzo8ERGRIqtOxVIMv6oLew6kcuUb09m5NyXSIWWrsAxClFXLpWdTnt0x/3tys0FmlmRmSZs3bz7CEEVEJCdVyybw0Llt+GpID9rWKc8/v5jPKU9NZMz8jfw+k5eIiIjkpZa1yvHyZR1ZsWUP172VxP6UtEiHdFgFnYBuDLrVErxuCsrXAHUz7VcHWBeU18mi/A/HmFkcUJ7/7fILgLu/4u6d3L1T1apV8+hSRETkcJrXKMebV3fhjSs7YwbXvZnEJa9OY966nZEOTUREpEg6vnEVHr+gHdNXbOOOD2YV2hHqCzoB/Qy4Ili/AhiVqfyiYGTbBoQGG5oedNNNNrOuwfOdlx9yTMa5zge+d91eFxEpNMyMPs2r8c1tPfnXWa1YuGEXpz87iT9/9Aubdu2PdHgiIiJFzlnta3PfgBZ8NWcD//xifqHsfZRvs4eb2XtAb6CKma0B/gE8DHxgZtcAvwIXALj7PDP7AJgPpAKD3T2j3fhGQiPqlgS+DhaAYcBbZraUUMvnRfl1LSIicuTiY2O4rFt9zmxfmxfGLeWNH1fyxez13NCrEdf1aEjJErE5n0RERETCcm2PBqzfuZ/Xf1xBzfKJXN+rUaRD+gMrjFlxfurUqZMnJSVFOgwRkWJr1dY9PPLNQr6as4Ea5RL5c/9mnN2+NjExhXcgczOb4e6dIh1HpKjuFBGJLunpzi0jZ/Ll7PU8PbA9Z3eonfNBeexwdWdhGYRIRESKiWMql+aFSzvy4Q3dqFYugTs++IWzX/iR6SuyfIxfREREcikmxnjywnZ0bViJuz/6hUlLtkQ6pN8oARURkYjoXL8S/72pO08NbMfm5ANc+PIUbnx7Bqu27ol0aCIiIlEvIS6Wly/rRMMqZbjh7RmFZiBAJaAiIhIxMTHGOR3q8P2dvbnzpKZMWLyZk56cyL+/WsDOfYV7HjMREZHCrnzJeIZf3ZmyiXFc+cZPrN62N9IhKQEVEZHIK1killv6NWH8Xb05p0NtXv1hOb0fG8ebU1aSkpYe6fBERESiVs3yJRlxdRcOpKRxxRvT2b7nYETjUQIqIiKFRrVyiTxyflu+uOUEmtcox99HzaP/0xP5fuHGQjmUvIiISDRoWr0sr13RmTXb93HNiJ/Yn5KW80H5RAmoiIgUOq1qlefd647jtcs74Q5XD0/i8tens3DDrkiHJiIiEpW6NKjE0wPbM3P1Dm55byZp6ZG5sasEVERECiUz48SW1fnmtp7844yWzF6zkwH/+YF7PpnNpuT9kQ5PREQk6gxoU5N/nN6SMfM38o/P5kakd1FcgX+iiIhILpSIi+Gq7g04p0Ntnv1+KSMmr+SzWeu4qU9jrjmhAYnxsZEOUUREJGpc2b0B63ft5+UJy6lRLpGb+zYp0M9XC6iIiESFCqVK8LfTWzLmjl6c0KQKj41eRL8nJjBq1lo9HyoiIpILfzmlOed0qM3j3y7mw6TVBfrZSkBFRCSqNKhSmpcv68R713WlQql4hoycxbkvTmbGqu2RDu2omVldMxtnZgvMbJ6ZDQnKK5nZGDNbErxWzHTMPWa21MwWmdkpmco7mtmcYNszZmZBeYKZvR+UTzOz+gV+oSIiElExMcYj57XlhMZVGPrJHMYt2lRwn11gnyQiIpKHujWqzOc3n8Bj57dl7fZ9nPfiZG5+9+dCMcfZUUgF7nT3FkBXYLCZtQSGAmPdvQkwNnhPsO0ioBXQH3jBzDL6JL8IDAKaBEv/oPwaYLu7NwaeAh4piAsTEZHCpURcDC/+6ViaVS/L4Hd+ZvaaHQXyuUpARUQkasXEGBd0qsv4u3szpF8TvluwkX5PTuDhrxeSvD8l0uHlmruvd/efg/VkYAFQGzgLGBHsNgI4O1g/Cxjp7gfcfQWwFOhiZjWBcu4+xUP9k9885JiMc30E9MtoHRURkeKlbGI8w6/qTKXSJbh6+E+s2ron3z9TCaiIiES9UiXiuP2kpoy/qw9ntK3FSxOW0fux8QX+XEteCrrGdgCmAdXdfT2EklSgWrBbbSDzRa4JymoH64eW/+EYd08FdgKVs/j8QWaWZGZJmzdvzqOrEhGRwqZauURGXN2F1HTn8tens2X3gXz9PCWgIiJSZNQon8gTF7bj85tPoFG1MqzZvi/SIR0RMysDfAzc5u7ZTX6aVculZ1Oe3TF/LHB/xd07uXunqlWr5hSyiIhEsUZVyzDsis6kpKazfkf+TnWmaVhERKTIaVOnPO8P6kpqhCbZPhpmFk8o+XzH3T8JijeaWU13Xx90r80YLWINUDfT4XWAdUF5nSzKMx+zxszigPLAtny5GBERiRodj6nIuLt7kxCXv9ObqQVURESKJDMjPja6qrngWcxhwAJ3fzLTps+AK4L1K4BRmcovCka2bUBosKHpQTfdZDPrGpzz8kOOyTjX+cD3rnlsREQE8j35BLWAioiIFCbdgcuAOWY2Kyi7F3gY+MDMrgF+BS4AcPd5ZvYBMJ/QCLqD3T0tOO5GYDhQEvg6WCCU4L5lZksJtXxelM/XJCIi8hsloCIiIoWEu08i62c0Afod5pgHgQezKE8CWmdRvp8ggRURESlo0dU3SURERERERKKWElAREREREREpEEpARUREREREpEAoARUREREREZECoQRURERERERECoQSUBERERERESkQVtzmnjazzcCqPDhVFWBLHpynMCgq16LrKHyKyrXoOgqfgr6WY9y9agF+XqESRt1ZlP5tFRT9zHJHP6/c0c8r9/Qzy51wfl5Z1p3FLgHNK2aW5O6dIh1HXigq16LrKHyKyrXoOgqfonQtRYG+j9zTzyx39PPKHf28ck8/s9w5mp+XuuCKiIiIiIhIgVACKiIiIiIiIgVCCeiReyXSAeShonItuo7Cp6hci66j8ClK11IU6PvIPf3Mckc/r9zRzyv39DPLnSP+eekZUBERERERESkQagEVERERERGRAqEEVERERERERAqEEtAjYGb9zWyRmS01s6GRjicnZrbSzOaY2SwzSwrKKpnZGDNbErxWzLT/PcG1LTKzUyIY9+tmtsnM5mYqy3XcZtYxuP6lZvaMmVkhuZb7zWxt8L3MMrMBhf1azKyumY0zswVmNs/MhgTlUfW9ZHMdUfWdmFmimU03s1+C6/i/oDyqvo8criWqvpPiyKKsToykw/3ukeyZWayZzTSzLyIdSzQwswpm9pGZLQz+rXWLdEyFmZndHvx/nGtm75lZYqRjKmwO83fsYf/WyJG7a8nFAsQCy4CGQAngF6BlpOPKIeaVQJVDyh4FhgbrQ4FHgvWWwTUlAA2Ca42NUNw9gWOBuUcTNzAd6AYY8DVwaiG5lvuBu7LYt9BeC1ATODZYLwssDuKNqu8lm+uIqu8k+MwywXo8MA3oGm3fRw7XElXfSXFbiMI6McI/ryx/90Q6rsK+AHcA7wJfRDqWaFiAEcC1wXoJoEKkYyqsC1AbWAGUDN5/AFwZ6bgK20Iu/iYPZ1ELaO51AZa6+3J3PwiMBM6KcExH4ixCv6AIXs/OVD7S3Q+4+wpgKaFrLnDuPhHYdkhxruI2s5pAOXef4qH/IW9mOqbAHOZaDqfQXou7r3f3n4P1ZGABoV/eUfW9ZHMdh1NYr8PdfXfwNj5YnCj7PiDbazmcQnstxUxRqRMLxBH87in2zKwOcBrwWqRjiQZmVo5QsjAMwN0PuvuOiAZV+MUBJc0sDigFrItwPIVOLv8mz5ES0NyrDazO9H4Nhb/ycOBbM5thZoOCsuruvh5CFSJQLSgv7NeX27hrB+uHlhcWN5vZ7KBrQ0bXhai4FjOrD3Qg1FIVtd/LIdcBUfadBF3TZgGbgDHuHrXfx2GuBaLsOylmCnudUWhl8btHsvY08GcgPcJxRIuGwGbgjaDb8mtmVjrSQRVW7r4WeBz4FVgP7HT3byMbVdQ43N8aOVICmntZPUtU2Oey6e7uxwKnAoPNrGc2+0bj9cHh4y7M1/Mi0AhoT+iX3hNBeaG/FjMrA3wM3Obuu7LbNYuyQnMtWVxH1H0n7p7m7u2BOoRaAFtns3uhvQ447LVE3XdSzOjnfQRy8Tu0WDOz04FN7j4j0rFEkThCXSVfdPcOwB5C3SMlC8FNzbMIPcpRCyhtZn+KbFRFnxLQ3FsD1M30vg6FvKne3dcFr5uATwl1mdoYdFUjeN0U7F7Yry+3ca8J1g8tjzh33xj8wZ0OvMrvXZ0L9bWYWTyhP5zecfdPguKo+16yuo5o/U4Agi5W44H+ROH3kVnma4nm76SYKOx1RqFzmN+hkrXuwJlmtpJQ9+6+ZvZ2ZEMq9NYAazL1IPmIUEIqWTsRWOHum909BfgEOD7CMUWLw/2tkSMloLn3E9DEzBqYWQngIuCzCMd0WGZW2szKZqwDJwNzCcV8RbDbFcCoYP0z4CIzSzCzBkATQgN6FBa5ijvoEpBsZl2DkTAvz3RMRGX8pw2cQ+h7gUJ8LcHnDgMWuPuTmTZF1fdyuOuItu/EzKqaWYVgvSShinQhUfZ9BPFneS3R9p0UQ1FVJ0ZaNr9DJQvufo+713H3+oT+bX3v7mqdyoa7bwBWm1mzoKgfMD+CIRV2vwJdzaxU8P+zH6FnsyVnh/tbI2fhjlak5Q8jQQ0gNHLdMuC+SMeTQ6wNCY1K+AswLyNeoDIwFlgSvFbKdMx9wbUtIoKjRwLvEepyl0Lojt41RxI30InQH63LgOcAKyTX8hYwB5gd/CeuWdivBTiBUPe62cCsYBkQbd9LNtcRVd8J0BaYGcQ7F/h7UB5V30cO1xJV30lxXIiiOjHSy+F+90Q6rmhYgN5oFNxwf1btgaTg39l/gYqRjqkwL8D/Ebp5OzeocxIiHVNhW8jl3+Q5LRacVERERERERCRfqQuuiIiIiIiIFAgloCIiIiIiIlIglICKiIiIiIhIgVACKiIiIiIiIgVCCaiIiIiIiIgUCCWgIiIiIiIiUiCUgIpEKTOrbGazgmWDma0N1neb2Qv58HnDzWyFmd2QzT49zGy+mc3N688XEZHCz8zSgrporpl9aGalcnFsezMbkOn9mWY2NIdjdh9NvGHGVT+jXjOzTmb2TB6cc7yZLTKzM3N53OSj/ezcMrNGGX9fFPRnS9GkeUBFigAzux/Y7e6P5+NnDCc0CfhHOexXP9ivdX7FIiIihZOZ7Xb3MsH6O8AMd38yjOPigD8Bndz95iP5vCOINdbd08LYrz55XK+Z2XjgLndPyqtzZvEZce6emofnO+KftUhmagEVKWLMrLeZfRGs329mI8zsWzNbaWbnmtmjZjbHzL4xs/hgv45mNsHMZpjZaDOrGcbnXBDc4f7FzCbm93WJiEjU+QFobGZnmNk0M5tpZt+ZWXX4rY56xcy+Bd4E/gkMDFrbBprZlWb2XLBvdTP7NKhzfjGz4w/9MDO728x+MrPZZvZ/WQUU9BL6p5lNA7qZ2d+DY+YGsViwX8fgc6YAgzMdf2gde1embXOD1tLSZvZlcPxcMxuY0w8qaBF9yswmmtkCM+tsZp+Y2RIzeyBz/JnW/xzU57+Y2cOZzvNvM5sADDGzfsHPfY6ZvW5mCcF+K83s/8zs52Bb86C8l/3eu2qmmZXNKXaR3FICKlL0NQJOA84C3gbGuXsbYB9wWpCEPguc7+4dgdeBB8M479+BU9y9HZCrLkQiIlK0BS2apwJzgElAV3fvAIwE/pxp147AWe5+CaF65X13b+/u7x9yymeACUGdcyww75DPOxloAnQB2gMdzaxnFqGVBua6+3HuPgl4zt07B62bJYHTg/3eAG51925HcPn9gXXu3i447zdhHnfQ3XsCLwGjCCW+rYErzaxy5h3N7FTgbOC44GfyaKbNFdy9F/A8MBwYGNT7ccCNmfbb4u7HAi8CGYn0XcBgd28P9CD0t4JInlICKlL0fe3uKYT+CIjl94pwDlAfaEaoghtjZrOAvwJ1wjjvj8BwM7suOK+IiEjJoC5JAn4FhhGqU0ab2RzgbqBVpv0/c/dwkpy+hBIl3D3N3Xcesv3kYJkJ/Aw0J5SQHioN+DjT+z5B6+yc4DNamVl5QknchGCft8KIL7M5wIlm9oiZ9cgi1sP5LNPx89x9vbsfAJYDdQ/Z90TgDXffC+Du2zJty0jemwEr3H1x8H4EkDkp/yR4nUHo7wEI1e1PmtmthH4GedaFVyRDXKQDEJF8dwDA3dPNLMV/f/A7ndDvACNU0eXqLq+732BmxxFqXZ1lZu3dfWteBi4iIlFnX9B69hszexZ40t0/M7PewP2ZNu/Jo8814CF3fzmH/fZnPPdpZonAC4SeO11tofEUEoNzhTNISip/bMxJBHD3xWbWERgAPGRm37r7P8M434HgNT3Tesb7Q/9mzy7GPZn2Cefz0jLO7+4Pm9mXQexTzexEd18YRuwiYVMLqIgsAqqaWTcAM4s3s1Y5HIOZNXL3ae7+d2AL/3t3VkREBKA8sDZYvyKb/ZKBwz1zOJag+6iZxZpZuUO2jwauNrOMAZBqm1m1HOJKDF63BMedD+DuO4CdZnZCsP3Swxy/klB3YMzsWKBBsF4L2OvubwOPZ+yTx74ldL2lgs+slMU+C4H6ZtY4eH8ZMCGL/X4T1O1z3P0RQq3YzfMwZhFACahIsefuBwlVuo+Y2S/ALOB/BnfIwmPBwAVzgYnAL/kXpYiIRLH7gQ/N7AdCNywPZxzQMmMQokO2DSHUXXYOoS6jf7hR6u7fAu8CU4J9PuLwyWzGMTuAVwl1ef0v8FOmzVcBzweDEB2ui/DHQKWgy/GNQEZX1zbA9KD8PuCBLI8+Cu7+DaEuu0nB59yVxT77CV3Hh8HPJJ3Q86XZuS0YOOkXQtf9dZ4GLoKmYRGRMJmmYREREckTVgDTsOQ10zQskkfUAioi4doJ/MvMbjjcDmbWA/ic7O9wi4iIFHfbCA3kV+hHkTezRkEr68ZIxyJFg1pARUREREREpECoBVREREREREQKhBJQERERERERKRBKQEVERERERKRAKAEVERERERGRAvH/yO2f1MOhcwUAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 936x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# solve\n",
    "solver = pybamm.ScipySolver()\n",
    "t = np.linspace(0, 3600, 600)\n",
    "solution = solver.solve(model, t)\n",
    "\n",
    "# post-process, so that the solution can be called at any time t or space r\n",
    "# (using interpolation)\n",
    "c = solution[\"Concentration [mol.m-3]\"]\n",
    "c_surf = solution[\"Surface concentration [mol.m-3]\"]\n",
    "\n",
    "# plot\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4))\n",
    "\n",
    "ax1.plot(solution.t, c_surf(solution.t))\n",
    "ax1.set_xlabel(\"Time [s]\")\n",
    "ax1.set_ylabel(\"Surface concentration [mol.m-3]\")\n",
    "\n",
    "r = mesh[\"negative particle\"].nodes  # radial position\n",
    "time = 1000  # time in seconds\n",
    "ax2.plot(r * 1e6, c(t=time, r=r), label=f\"t={time}[s]\")\n",
    "ax2.set_xlabel(\"Particle radius [microns]\")\n",
    "ax2.set_ylabel(\"Concentration [mol.m-3]\")\n",
    "ax2.legend()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the [next notebook](./4-comparing-full-and-reduced-order-models.ipynb) we consider the limit of fast diffusion in the particle. This leads to a reduced-order model for the particle behaviour, which we compare with the full (Fickian diffusion) model. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "The relevant papers for this notebook are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "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] 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",
      "[3] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). Journal of Open Research Software, 9(1):14, 2021. doi:10.5334/jors.309.\n",
      "[4] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, and others. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods, 17(3):261–272, 2020. doi:10.1038/s41592-019-0686-2.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "pybamm.print_citations()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pybamm",
   "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.9.16"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  },
  "vscode": {
   "interpreter": {
    "hash": "187972e187ab8dfbecfab9e8e194ae6d08262b2d51a54fa40644e3ddb6b5f74c"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}