docs/source/examples/notebooks/creating_models/1-an-ode-model.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Creating a simple ODE model\n",
    "\n",
    "In this series of notebooks, we will run through the steps involved in creating a new model within pybamm. Before using pybamm we recommend following the [Getting Started](../getting_started/tutorial-1-how-to-run-a-model.ipynb) guides.\n",
    "\n",
    "In this notebook we create and solve the following simple ODE model:\n",
    "\n",
    "$$\n",
    "  \\frac{\\textrm{d}x}{\\textrm{d} t} = 4x - 2y, \\quad x(0) = 1,\n",
    "$$\n",
    "\n",
    "$$\n",
    "  \\frac{\\textrm{d}y}{\\textrm{d} t} = 3x - y, \\quad y(0) = 2.\n",
    "$$\n",
    "\n",
    "We begin by importing the PyBaMM library into this notebook, along with NumPy and Matplotlib, which we use for plotting:\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"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up the model"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first initialise the model using the `BaseModel` class. This sets up the required structure for our model. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = pybamm.BaseModel()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we define the variables in the model using the `Variable` class. In more complicated models we can give the variables more informative string names, but here we simply name the variables \"x\" and \"y\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = pybamm.Variable(\"x\")\n",
    "y = pybamm.Variable(\"y\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now use the symbols we have created for our variables to write out our governing equations. Note that the governing equations must be provied in the explicit form `d/dt = rhs` since pybamm only stores the right hand side (rhs) and assumes that the left hand side is the time derivative. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "dxdt = 4 * x - 2 * y\n",
    "dydt = 3 * x - y"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The governing equations must then be added to the dictionary `model.rhs`. The dictionary stores key and item pairs, where the key is the variable which is governed by the equation stored in the corresponding item. Note that the keys are the symbols that represent the variables and are not the variable names (e.g. the key is `x`, not the string \"x\")."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.rhs = {x: dxdt, y: dydt}"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The initial conditions are also stored in a dictionary, `model.initial_conditions`, which again uses the variable as the key"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.initial_conditions = {x: pybamm.Scalar(1), y: pybamm.Scalar(2)}"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can add any variables of interest to our model. Note that these can be things other than the variables that are solved for. For example, we may want to store the variable defined by $z=x+4y$ as a model output. Variables are added to the model using the `model.variables` dictionary as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.variables = {\"x\": x, \"y\": y, \"z\": x + 4 * y}"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the keys of this dictionary are strings (i.e. the names of the variables). The string names can be different from the variable symbol, and should in general be something informative. The model is now completely defined and is ready to be discretised and solved!"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the model"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first discretise the model using the `pybamm.Discretisation` class. Calling the method `process_model` turns the model variables into a `pybamm.StateVector` object that can be passed to a solver. Since the model is a system of ODEs we do not need to provide a mesh or worry about any spatial dependence, so we can use the default discretisation. Details on how to provide a mesh will be covered in the following notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "disc = pybamm.Discretisation()  # use the default discretisation\n",
    "disc.process_model(model);"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that the model has been discretised it is ready to be solved. Here we choose the ODE solver `pybamm.ScipySolver` and solve, returning the solution at 20 time points in the interval $t \\in [0, 1]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "solver = pybamm.ScipySolver()\n",
    "t = np.linspace(0, 1, 20)\n",
    "solution = solver.solve(model, t)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After solving, we can extract the variables from the solution. These are automatically post-processed so that the solutions can be called at any time $t$ using interpolation. The times at which the model was solved at are stored in `solution.t` and the statevectors at those times are stored in `solution.y`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_sol, y_sol = solution.t, solution.y  # get solution times and states\n",
    "x = solution[\"x\"]  # extract and process x from the solution\n",
    "y = solution[\"y\"]  # extract and process y from the solution"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then plot the numerical solution against the exact solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6AAAAEYCAYAAABCw5uAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABYyUlEQVR4nO3dd3hUVf7H8fdJMmmkQQotQOg9SEdFpFiwoVhYbFgQrGvb/a11ERFdd11FEXtDxQJ27CAWFGkB6UhvoYYWAqTn/P6YwFImkDKZO5l8Xs8zz8zce3PvZ644Z763nGOstYiIiIiIiIhUtiCnA4iIiIiIiEj1oAJUREREREREfEIFqIiIiIiIiPiEClARERERERHxCRWgIiIiIiIi4hMhTgc4kYSEBJuSkuJ0DBERkaPMmzdvp7U20ekcaidFRMRfldRW+nUBmpKSQlpamtMxREREjmKM2eB0BlA7KSIi/quktlKX4IqIiIiIiIhPqAAVERERERERn1ABKiIiIiIiIj7h1/eAioj4g/z8fNLT08nJyXE6ivhYeHg4ycnJuFwup6OIiPg1tZXVV1nbShWgIiInkZ6eTnR0NCkpKRhjnI4jPmKtZdeuXaSnp9O4cWOn44iI+DW1ldVTedpKXYIrInISOTk5xMfHq0GtZowxxMfH62i+iEgpqK2snsrTVnqlADXGvGmM2WGMWVLCfGOMGWuMWW2MWWSM6eSN7ZbKokkwph2MjHM/L5rkX+sTkSpBDWr1pP/uIj6k31hVnr4zq6ey/nf31hnQ8UD/E8w/D2he/BgOvOSl7Z7Yoknw5Z2QuQmw7ucv7yz/F5q313fkevWFKyIiIlWJN3+/VNZvLBHxO165B9RaO90Yk3KCRS4G3rHWWmCWMSbOGFPXWrvVG9sv0bRRkJ999LT8bHZPfph/rWhJkDEEBbmr9iCD+70xGAOG4mlBBlewwRUcxI1z/0mMh/Ud/PYRpgf1IjTEvdyhR4QrmIjQYPezK5jw0CBCg4OOPkpw6Av30HoPfeECpA4q3+deNMn92TPTITYZ+o0o/7pExHGbNm1iyJAhbN++HWMMw4cP56677jo8f/z48fTu3ZtGjRp5/ejzs88+S61atRgyZAjjx4/nnHPOoV69egAMHjyYxx57jObNm3t1m6VlraVfv358/vnnZGZmlmkf+dtnEal0XvxtUFBYRM78D4n8/l6CCv73+6Xgi78yc1UGfyaeR25BIbkFRe5Hvvt1QZE9vA5rj17nQ6seppaH31j7vxnBVzk9iIlwERPuIiYihNji19HhIYQEn+Bcin4PVSs5OTn06tWL3NxcCgoKuPzyy3n00Uex1mKMYeTIkYwcOfLwe2+6++67ufTSS+nVqxfPPvssw4cPJzIyEoCzzjqLjz76iJo1a3p1m6WVnZ1N//79+fHHH1m8eDG33nor+/btIzg4mIceeoi//OUvJe6jyvosxh77DVDeFbkL0K+ste08zPsKeNJa+1vx+2nAfdbaNA/LDsd9lpSGDRt23rBhQ/lDjYwDjv98RRh6hn1CkQWLdT9b93ORtRQVWeyh1xYKiorIL7SsDbuKIA//XousoUnue6WKFGQoLkxDiAgN4uPs4dS2Gccttze0Dm93+4ro8BCiwkOICQ8hOtxFVFgI0cWvo8NDCAs5SUEL4IqAi8bqS1eknJYvX07r1q0d2/7WrVvZunUrnTp1Iisri86dO/P5558TGxvLI488QqNGjWjcuDG//vorr7zyite2W1BQQKdOnZg/fz4hISH07t2b//73v3Tp0gWAX375hQkTJvDaa695bZtl8fXXX/PDDz8wZsyYMu+jsnwWT//9jTHzrLVdfPJBT6BLly42Le24plTkaKX4bZBbUMi2zBy2Zuawc38uO7Ny2bk/j537c8nIynVP25/H3oN5HMgr5LfQO0kO2nncptKLEuiZNxYAYyA8JJgwl/sAvOsExeKvOZcS5Ok320l+Y0WHh5AUHUZSdDhJMWGHX3fMnErHBY8QXKjfQ77idFtpreXAgQNERUWRn59Pz549ee655wgLC+Ott94CoG/fvsyZM4cnnnjCa9vdtWsXF1xwAbNmzQIgJSWFtLQ0EhISAHj77bdJT0/noYce8to2y+KFF16goKCAu+66i5UrV2KMoXnz5mzZsoXOnTuzfPly1q1b53EfleWzlKWt9LtecK21rwKvgrthrdDKYpOLL+U4WlBsMr/f06+sueDZZPdRtGMURtfjm5vPIL+wiPzCIvIKi8grKCInv4ic/EIO5hWSnV9ITn4h2cWvs/MLyckrJGn58V/eADG52xnzw8qT5nIFG2IjQqkZ6aJmZCgv7XyI+ILjjyDmfj+STXXOp1aNUOIiXAR5qqRLoiOIIo6qW7cudevWBSA6OprWrVuzefNm2rRpw+OPP0737t1p164dkydPBmDNmjXcfvvtZGRkEBkZyWuvvUazZs049dRTeeqpp+jduzcPPPAAQUFBPP7446SkpDBo0CC+/fZbIiIieP/992nWrBk//vgjnTp1IiQkhI8//pi0tDSuvvpqIiIimDlzJmeccQbXX389BQUFhISUvjmZMGECY8eOJS8vj+7du/Piiy8yf/58hg4dypw5cygsLKRbt25MnDiRnTt3MmLECKKjo1m9ejV9+vThxRdfJCgoiPfee4/hw4eXeR9587OIVAklXBG2a/LDDP21Ppv3ZpORlXvcnwUHGWrVCCUhKozE6DCaJkVRMzKUmHAX9X/b5XFT9YN2sXjkOYSFBOMKNqU/0zTG8282Yuvz+4192ZeTz77sAvZl55OZnc++HPfzngN57MjKZUdWLvM37mHHvlxyC4r4LfQ/BAcd/5n3fvkwn+7rSvPaUTRPiqZ2TJjuWwwQxhiioqIA95Aw+fn5GGPo2LEjERERnHrqqeTn5/PSS+47AefNm8e9997L/v37SUhIYPz48URGRtKtWzcmT55My5YtufLKK+nbty/Dhg0jKiqKYcOGMWXKFOrUqcOHH35IYmIin3zyCf37u+9EHDt2LFu2bKFPnz4kJCTw008/MWDAAM4444wyF6BPPfUUkyZNIjc3l4EDB/Loo4/y2WefMW7cOH744Qe2bdvGmWeeyfTp0/nuu+/47LPPyMzMZPPmzVxzzTU88sgjALz33nu8//77ALRo0eLw+uvVq0dSUhIZGRke95E3P8uxfNXKbgYaHPE+uXha5eo3wvMRv34jyrwqYwz0e8Tj+lznjKRNvZjyZSzhCzcoLplVd57HgdwCsnIOPfLJyilgf27x6+J5ew+6v4D3HMyjZsHxZ1MBXPu3cNYzvwDuBiUhKpTE6DASixuV/70O/9/76DBqrPgU4+1LhEWqsEe/XMqyLfu8us429WJ45KK2pVp2/fr1/PHHH3Tv3p0tW7bwyCOPcOONN9K4cWNuv/12XnrpJYYPH87LL79M8+bNmT17Nrfddhs//vgj48eP5/LLL+f555/nu+++Y/bs2YfXGxsby+LFi3nnnXe4++67+eqrr5gxYwadO3cG4PLLL2fcuHFHnTUEaNasGQsXLjy83MksX76ciRMnMmPGDFwuF7fddhvvvfceQ4YMYcCAATz88MNkZ2dzzTXX0K5dO37++WfmzJnDsmXLaNSoEf379+fTTz/l8ssvZ8aMGR7P+JZmH3njs4j4m7yCIjbuPsj6nQdYt/MA63YdYP3OA0zITPfY6UfNgh1EhYXQp2Ui9eMiqRcXTt3YCBKiQ0mMCqNmZGjJB6wXe/79YmKTiQ4vx7i5JfxmCzrrEerFRVCPiFKtxlpLVm4B0U96LpBj8nYw6qtlh99Hh4fQonY0zZOiaFE7mnb1Y2lbL4YaYR5+IuuAfKk51VYWFhbSuXNnVq9eze2330737t1ZsGABb731Ftdeey19+/bl4Ycf5pFHHuGvf/0rX3zxBYmJiUycOJGHHnqIN998k3HjxnH99ddz1113sWfPHoYNGwbAgQMH6NKlC2PGjGHUqFE8+uijjBs3jhkzZnD55ZcDcOedd/LMM8/w008/HT5rWLNmTXJzc9m1axfx8fGl+qxTpkxh1apVzJkzB2stAwYMYPr06QwcOJBPPvmEF154ge+++45HH32UOnXqADBnzhyWLFlCZGQkXbt25YILLiA1NZW1a9eSkpJy3DbmzJlDXl4eTZs29biPRo8e7ZXP4omvCtDJwB3GmA+B7kBmpd//Cf/7UvDWl4W31wcnLJJdwUHERYYSFxla+vWVUNDm1ajLcwNPYfeBPHbtzyMjK5eM4stqlm/NYuf+3KPuzThkRtgD1DfHH0HM+e4R0uucR93YCM9f0iLidfv37+eyyy7j2WefJSYmhpiYGF577TXGjx/PGWecwTXXXMP+/fv5/fffueKKKw7/XW6u+8xG27Ztufbaa7nwwguZOXMmoaH/+2658sorDz/fc889gPvS35NdTpWUlHT4Mp7SmDZtGvPmzaNr166A+96UpKQkAEaMGEHXrl0JDw9n7Nixh/+mW7duNGnS5HC+3377jcsvv5zdu3cTHR1d5n3krc8iUqlOUOwUFlk27DrAyu1ZrNi2n5Xbs/hz2z7W7zpI4RFteVyki8YJNdgbkkStgu3HbSIoNpkJN3UvXz4vHuQHvPYbyxhDTLirxKvgTFwyacPOYuX2LFbvcO+7ldv38/3SbXw4d1PxOqBZYhTtk2NpXz+W1ORY2u+eQug3d+uAvJ8LDg5mwYIF7N27l4EDB7JkyRI6dOjAc889x8iRI7nkkku4+OKLWbp0KUuWLOHss88G3IXroatozj77bD766CNuv/12Fi5ceHjdQUFB/OUvfwHgmmuu4dJLLwXcbWViYuIJcx1qX8pSgE6ZMoWOHTsC7rZt1apV9OrVi+eff5527drRo0ePw233odyH1n/ppZfy22+/Ua9ePeLi4o5b/9atW7n22mt5++23CQoK8riPvPVZPPFK5WCM+QDoDSQYY9KBRwAXgLX2ZeAb4HxgNXAQuMEb2y2V1EHe/WKojPWB94raEhqE8P6PcnFq/RL/rKjIsjc7nx1ZOe7itPhR7yfPRxBDD2zlrGemAxATHkK9uAjqxobTsFYkDWpF0rBWJI3ia9CgVgSRocf8M9MRRKnCSnum0tvy8/O57LLLuPrqqw83eodcf/31h18XFRURFxfHggULPK5n8eLFxMXFsWPHjqOmH3kJ2qHXERERJx3XKycnh4iIo89MzJ49m5tvvhmAUaNGMWDAgMPzrLVcd911/Otf/zpuXbt27WL//v3k5+eTk5NDjRo1jst25PuQkBCKiooICnKf2yntPirLZ/EWY0wD4B2gNu7OCV611j53zDJXA/cBBsgCbrXWLjx2XVINeOigMP/zvzJx1gY+zO3Bqu37yS0oAtzFUqNakbSoHc157erSJLEGjRPcj8MHsBeN9m6xCJVzUN6bv7FK+D1k+o0gISqMhKgwTmuacHiWtZaMrFyWbMlkUXomi9Mz+XXVTj6d775gb0bYQx4PyDNtlH7DeOBUW3lIXFwcffr04bvvvqNdO3f3NCNHjgTcbYi1lrZt2zJz5szj/raoqIjly5cTGRnJnj17SE5O9riNiraVn332GY8++igAr7/++lFX5FhreeCBBw63pUdKT08nKCiI7du3H9UGemorPWXbt28fF1xwAY8//jg9evQ46m+P3Edl+Sxl5a1ecK88yXwL3O6NbQUkb37hlrNBCCq+z6NWjVBa1TlixnzPRxDzo+rx3MBT2LI3h62Z2WzZm8OWvdmkrd9DVm7BUcsmRofRsLgoPbvgF85Z+wQhhcX/M+gIoshJWWsZOnQorVu35t577z3hsjExMTRu3JiPPvqIK664AmstixYtokOHDnz66afs3r2b6dOnc+GFFzJnzpzDR0YnTpzI/fffz8SJEzn11FMBaN26NatXrz687ujoaLKyso7a3sqVKw837occuuTJk379+nHxxRdzzz33kJSUxO7du8nKyqJRo0bcfPPNPPbYY6xbt4777ruPcePGAe7LhNatW0ejRo2YOHHi4fs+W7Zsydq1a2nWrFmZ9lFZPosXFQB/s9bON8ZEA/OMMVOttcuOWGYdcKa1do8x5jzc/SGU8/SUVEV7D+axZPM+Onw9guhj7tl0FeXQb8srTGl4JkNOjadF7Wha1ommWVLU8Qd6j1UZxeKh9fpr213Gz2yMISkmnL4x4fRtVRtwf/du35fLwvS91PvI8wF5m5nOb6sy6JpSi3BXcKV8FCmdjIwMXC4XcXFxZGdnM3XqVO677z6Py7Zs2ZKMjAxmzpx5+L7HlStX0rZtW8aMGUPr1q154oknuOGGG5g5cyYul4uioiI+/vhjBg8ezPvvv0/Pnj2B/7WVvXv3Bv7Xvhy6bNVay7Zt2467DHbgwIEMHDjQY75zzz2Xf/7zn1x99dVERUWxefNmXC4XtWrV4sYbb+SDDz7g7bff5plnnuHvf/87AFOnTmX37t1ERETw+eef8+abb1KzZk0KCwvJyckhPDycvLw8Bg4cyJAhQw5fNnwipf0sZaVrJwORD44ghp070uMZVWstew/ms3H3QTbsPsim3QfZuOsgG3YfYM663fwteywh5pijRPnZ7PnyYT7a25mmiVE0SYyiQc2IE3etLlKNzJgxg3fffZf27dtzyimnAPDEE09w/vnne1z+vffe49Zbb2X06NHk5+czePBg6tevz/3338+0adNo0KABd9xxB3fddRdvv/02AHv27CE1NZWwsDA++OADAM477zyuvfbaw+u9/vrrueWWWw533LNv3z4iIiIO339SGm3atGH06NGcc845FBUV4XK5eOGFF/jll19wuVxcddVVFBYWctppp/Hjjz8SFBRE165dueOOOw53QnSowb7gggv4+eefadasWZn3kTc+S1kU33aytfh1ljFmOVAfWHbEMr8f8SezcPeXIAEqJ7+QxZszmb9hDws27WXx5kzS97jb2rVhW93nwY9Rl528c2O38m3Qn4vFylLBz2yMoU5sOHVi68AUzwfkt9h4rn1jDuGuIHo0iadX80TObJlIk4Qa6tzIx7Zu3cp1111HYWEhRUVFDBo0iAsvvNDjsqGhoXz88cfceeedZGZmUlBQwN13301ISAivv/46c+bMITo6ml69ejF69GgeffRRatSowZw5cxg9ejRJSUlMnDgRcLdFr7zyCjfddBMAw4cPp3///tSrV4+ffvqJefPm0aNHjzJ1cHfOOeewfPnywweEo6KimDBhAi+//DJnnHEGPXv2pEOHDofv9QT37SqXXXYZ6enpXHPNNYfPqJ5zzjn89ttvnHXWWUyaNInp06eza9cuxo8fD7iHKTvUbh7LG5/FE68Nw1IZ1L28n/DiJbN2ZBymhKFxmuT8r5t1V7AhJb4GTROjaFE7ipZ1YmhVN5qU+BoEe+oQQZf1SiVyumv5ynZsN+tHGjhwIP/5z388jo85ZswYYmJiGDp0aKVl+/nnn/nvf//LV199ddy8rVu3MmTIEKZOnVrh7Zzos3h7GJbiYcumA+2stR576TDG/B1oZa29ycM87w1XJt5zgnbIWsuWzBzmb9jD/I17mL9xL8u2ZJJf6G4PG9SKILV+HO3qu+85PPXLMwned3yv+8Q2gHuW+PJTySElDGWTe96zzIjsw/SVO5m+MoO1Ow8A0LBWJOe2rc25bevQqWHNso0+UEUFelsZFRXF/v37Pc7r2bMnX331lcf7Le+66y4GDBhAv35lG4GjLMaPH09aWtrhK4eONH/+fMaMGcO7775b4e2c6LNU6WFYxA958aipOcHQOAv/cQ5rdu5nzY79rMk4wJoMd+cAU5Zt41CfCmEhQYcvO2pVJ5rWdWNov3sKMVPvVccAIpXgySefZOvWrR4L0Li4uKPOkPpa3bp1GTZsGPv27SMmppw9kRfz1WcxxkQBnwB3n6D47AMMBXp6mu/V4crEOzzcs1n0xZ38/OcOPs4/lXkb9rB9n7szsHBXEKnJcQzt2YRODePo2LAmidFhR6/vLM+97lfonk2pmBIu6Q1LHURfOHzZ7qbdB/l5ZQbTlm9n/O/ree3XdSREhXFOcTF6apN4QkN0hVegefrpp9m4caPHArRdu3aVWnyeTKdOnejTpw+FhYUEB1fsMnFvfRadARXfKsVg2MfKyS9k9Y79/Lktiz+37mPF9qzDPfcCJQ6GbWMbYHSkWLwg0I/qyol56wyoMcYFfAV8b619poRlUoHPgPOstScdDFrtpH+wY9piPIwTnl6UwOAar9GpYU06NYyjc6NatKobjas0t5joyp4qb19OPj/9uYMpS7fz04odHMwrpGakiwtS63LJKfXp3KgmZvFHAfPfWW1l9aYzoOK/ytERQrgrmHb1Y2lXP/ao6Tv357JiWxb1J5TcMcB1b84htfhv2yfHUi82XPdkiIjPGfcXzxvA8hMUnw2BT4FrS1N8inPyC4tYlJ7JrLW7mL1uN+P3bsZT01I/aBe/3de3fBupjvdsBpiYcBcXn1Kfi0+pT05+Ib+u2snkhVv4eF46E2Zt5IboOTxY+DKuInXMKNWLClDxPS81qglRYSQ0CytxrK9MVxI79uXw0uqdh8dFi68RSseGcXRqVJNODWvSITmOiFAPlyPoyLOIeNfpwLXAYmPMguJpDwIN4fCQZSOAeODF4gNlBeW9z1S8y1rLyu37+XVVBr+t3smcdbs5mFcIQPOkKDJDk6iZf/w4myZW/UiJW7grmLPb1ObsNrXZn1vA90u2ceY3d/+v+DxEQ7tINaACVKq+EnrqrXnRaL5L7UVOfiHLt+5j8eZMFm7K5I+Ne/hhuXsMxJAgQ+u6MXQ6oihNTv8Kc8y9PDoiKSIVYa39DY/9mh61zE3AcZ0OiTN27Mvht9U7+W3VTn5bvZMdWe7bPpok1OCyTsmc1jSebo1rER8VVjnjbErAigoL4bLOyfBlhsf5NjOdFdv20apOxe5tF/FXKkCl6jvJZb3hrmA6NqxJx4Y1wd2bNbsP5PHHxuLeCDfsZVJaOm/PdPckOTP8AeqiwaZFRAKWh6tcsltdxqx1u9wF56qdrNjuHie2ZqSL05slcEbzBHo2T6R+nIcB2CtrnE0JbCVcwbXFxtP/2V/p2DCOG05vzHnt6pTuvmGRKkIFqASGMl7WW6tGKP1a16Zfa3evdQWFRfy5LYv5G/dQ5/uS7yldsW0fLZKiq0V36iIiAclDj7W5n93BQwUL+TT/NEJDguiaUpNLOrbijOYJtKkbU7rvfN2zKWVVwhVcsec8xj9z2zBh1gbu/OAP6saGc91pKVzZtSGxkS7n8op4iQpQESAkOOh/HR3N8nxEcnOR+4hkzUgX3RrXokeTeM5onkjTRA02LcfQPcQifimvoIii7x8hPP/oq1zCbC4jIj7mkmvvplvjWoS7KjZUgUiplHDmPCp1EEOBG05L4acVO3jjt3U8+e2fPPfDKq7okswNpzemcUINR6N7hdrKaksFqMixSjgiGdH3Uf4b2oFZa3cxa+0uvl/q7nCiXmw4ZzRP5IwWCfRslkBcZKhDwcUveDi7UtF7iOfOncvQoUOZM2cOhYWFdOvWjYkTJ9KuXTsvhRYJXDv25fDTih389Ke7A6FFbPF4N25c3nZ6tUj0fUCp3k5w5jwoyBy+WmvZln28OWMdH87ZxIRZG7ggtR6392lade8TrYS2csSIEdSqVYu7774bgIceeoikpCTuuusuLwQWb9I4oCKelOKo3KbdB/l11U6mr8xgxpqdZOUUYAykJsfRq3kCvVokckqDOPd9GzrKV6WVaWyzMe08nkEntgFUYFzahx9+mJycHLKzs0lOTuaBBx4o97qkbLw1DmhlUDt5PGstK7ZnMWXpdqYu287izZkA1I0Np0+rJP65ahARB7cc/4cV/H9UxBd2ZOXw1oz1vPP7eg7kFXJOm9rc0bcZqclxTkdzvK1cv349l156KfPnz6eoqIjmzZszZ84c4uPjy7U+KRuNAypSUaW4l6dBrUiu6t6Qq7o3pKCwiIXpmfy6KoPpKzN44afVPP/jaqLDQri3zkKu3fk0IYUa56ta8DAY/Qmnl9KIESPo2rUr4eHhjB07tkLrEgk0BYVFzNuwhynLtjNl2TY27c7GGOjYII5/9G9J31ZJtKwd7b5dYtGj6rFWqqyk6HDu69+Km3s14a0Z63lrxjqmLNvOmS0SuefsFpyyZ0rVOOBdCW1lSkoK8fHx/PHHH2zfvp2OHTuq+PRTKkBFvCAkOIjOjWrSuVFN7j6rBZnZ+cxcs5Of/syg/5LbCOH4cb7stFEYf2wUpGJK6NWQCo4HuGvXLvbv309+fj45OTnUqBEA9/+IlFYJvdb+uiqDKcu2M235dvYczCc0OIjTm8VzW+9m9GudRFJ0+PHrUo+1EgDiIkO55+wW3HRGY96dtYHXf13Hmy/9m6fC3iDMuocM8usD3pXUVt50002MHz+ebdu2ceONN1ZoXVJ5dAmuSCWzI+MwHP//WRGG0Z1/56w2SXRLqUWIulj3W2W6rOjY+1rAfXblorEV+gEwYMAABg8ezLp169i6dSvjxo0r97qkbHQJrsM8/D+Va8J4sHAYn+SdRkx4CH1bJXFO2zr0apFIVJiOrUv1sz+3gKJn2hKTu+34mT66vNwf2sq8vDzat29Pfn4+q1atIjhYHYr5ii7BFfEjpoSjfLuDE5kwewNvzlhHfI1Qzmlbh/Pb16FHk3iN91WVVcLZlXfeeQeXy8VVV11FYWEhp512Gj/++CN9+/b1UmgR/1X4w6MEe+i19p9hH3HpkHvo1riWvjOl2osKC4Hc7R7n2cx0T/1uOauSrkQIDQ2lT58+xMXFqfj0YypARSpbCb3qJlz0OAtanc0vKzL4Zsk2Ji/YzAdzNhIX6eKcNrU5r31dTm+aQGiIflhVOV4eD3DIkCEMGTIEgODgYGbPnu21dYv4o537c/luyTa+WbyVCZnpnnutzd/B6c0SfB9OxF+VcMB7K/FMm7WBK7s28K+rrSph7NyioiJmzZrFRx995NX1inepABWpbCc4yhcJnNe+Lue1r0tOfiHTV2bw7ZJtfLt4G5PS0okOD+Hs1rW5qEM9ejZP0FF+EQlYO7Jy+H7JNr5evJU563ZTZKFJQg2ywuoQm+fpssKK3SsmEnA8HPAuCongk6gbefrzJUyYuYF/XtiGns0D88DNsmXLuPDCCxk4cCDNmzd3Oo6cgApQEV8oxVG+cFcw57Stwzlt65BbUMiM1Tv5dvE2pizbzqd/bCa+RigXpNblko716dggDrP4I3WiISJV2o59OXxbfKZzzvrdWAtNE2twR59mnJ9a191z7eLH1GutSGl4OOAd1G8Ed7S/guZLt/H4N8u55o3ZnNW6No9c1IYGtSKdzetlbdq0Ye3atU7HkFJQASrih8JCgunbqjZ9W9Xm8YIiflmZwecLNjNx7ibembmBG2Pm8kDBS7iKNLSLr1hr3UM4SLXizx31VVV7DuTxzZKtTF6w5XDR2aJ2FHf2bc4FqXVpUTv66D9Qr7UipefhgLcB+rerS++WSbw1Yz3P/7iKs8f8wp39mjPsjCZevbpKbWX1VNa2Ur3gilQhWTn5fLdkG72/7Uti4Y7jF9BA6pVi3bp1REdHEx8fr4a1GrHWsmvXLrKysmjcuPFR89QL7kkcM2xKzpkP8505g8kLtzB9ZQYFRZamiTW4qEM9LkytS7Ok6JOvU0S8YsvebEZOXsqUZdtpUTuKJwa2p0tKrQqvV21l9VSetlIFqEhVNDIOPAztYjH8fOVKejVPJDhIX/7ekp+fT3p6Ojk5OSdfWAJKeHg4ycnJuFyuo6arAD0BD8MrZNtQ7su/ibTos7jolHoM6FCPNnVj9CNVxEFTl23nkS+WsCUzh8FdG3D/ea2Iiwwt9/rUVlZfZW0rdQmuSFV0gp7ubnhrLnViwrm8czJXdEmmUXwNBwIGFpfLddxRPRE5XmGRJf/7Rwg/ZtiUCJPHUzW/wPW3JwjSwTERv3B2m9qc1jSe56at4o3f1jHtzx08MbA9Z7epXa71qa2U0lKXmiJVUb8R7k44juSKIOmSJ3jp6k60qhvNiz+v5synfmbwqzP57I90cvILnckqIgHNWssfG/fw6JdL6fGvaYTu3+JxubADW1R8iviZGmEhPHh+a764/XTia4Qy7J007pm4gL0H85yOJgFMZ0BFqqISOuUISR3EebiHdtmamc0n89KZlJbOPRMX8thXyxnUpQFXd28YcD3fiYjvbdx1kM/+2Mxnf6SzftdBQoOD6NMqkZwt9YjM9lCEatgUEb/Vrn4sk+/oybgfV/HCz2uYsXonTwxsz1nlPBsqciK6B1QkwBUVWX5fs4sJszYwdfl2iqylT8skru3RiDNbJOqMhEg5VNd7QDOz8/lm8VY+nZ/O3PV7ADi1STwDO9bn3HZ1iI1webwHFFcEXDRWPdeKVAFLNmfy948W8ue2LEY3Wc5V+98iaN9m9UAtZaZOiESErZnZfDB7Ix/M3URGVi4NakVwTfdGXNGlAbXWfK5hDkRKqToVoPmFRUxfmcGn8zczdfl28gqKaJpYg0s7JXNJx/rUj4s4/o+O6QVX3yciVUteQRFTJz5Pn5WjiTRHXI6rg0lSBpVagBpj+gPPAcHA69baJ4+Zfz3wFLC5eNI4a+3rJ1uvClCRypFXUMSUZdt4d+YGZq/bzaWu33ky5DVCbe7/FlIjI1KishagxpgGwDtAbdxdWL9qrX3umGUM7rb0fOAgcL21dv6J1uuVdtJDsWjbX8GSzfv4ZH46Xy7cwq4DedSqEcqADvUY2LE+qcmx6sFWJNCNaeexw0MN+SalVWm94BpjgoEXgLOBdGCuMWaytXbZMYtOtNbeUdHtiUjFhYYEcWFqPS5MrcfK7Vkkvn4Pofm5Ry+Un+3+UaoCVMQbCoC/WWvnG2OigXnGmKnHtJXnAc2LH92Bl4qfK8+xl8tmbiL/87/y1DfLeXVvF0KDgzirTRIDOyZzZotEQkPUd6FItZGZ7nGyzUxHh5+kIrzRCVE3YLW1di2AMeZD4GLg2AJURPxQi9rRkL/D4zybmU5eQSFhIcE+TiUSWKy1W4Gtxa+zjDHLgfoc3VZeDLxj3ZcmzTLGxBlj6hb/beWYNuroezUBV1EON+VNIGXgDVzQvi6xka4S/lhEAtoJhnxbumx7uYdrEfHGocz6wJH/OtOLpx3rMmPMImPMx8WXInlkjBlujEkzxqRlZGR4IZ6InFQJvVNuLorn9Cd/YtyPq8g8mO/jUCKByRiTAnQEZh8zq1TtqVfbyRLOcCQVZXBV94YqPkWqMw9DvhWFRPBejesY9k4aT377JwWFRQ6Fk6rMV9fSfAmkWGtTganA2yUtaK191VrbxVrbJTEx0UfxRKo5D42MdUVw4IwHaVc/hv9OWcnp//6RJ7/9k4ys3BJWIiInY4yJAj4B7rbW7ivPOrzaTpY0NIqGTBGR1EHuviBiGwAGYhsQNGAsf737Ia7q3pCXf1nDNW/M1u8CKTNvXIK7GTjyjGYy/+tsCABr7a4j3r4O/McL2xURb/EwrqjpN4KWqYMYDyzbso+XflnDq9PX8NaMdfylawOG92pCck2NJypSWsYYF+7i8z1r7aceFjlpe+p1/UZ4HjKl34hK3ayIVBGpg47rCyIceGJgezo1rMlDny3mgrG/8sLVneiaUsuZjFLlVLgXXGNMCLAS6Ie7oZwLXGWtXXrEMofvYTHGDATus9b2ONm61QuuiH9Zt/MAr/yyhk/mp2MtXHxKfW7t3YRmSdFORxPxqXL0gmtwX/2z21p7dwnLXADcgbsX3O7AWGtttxOtt7J6wVXnYyJSGsu37uPWCfPYtCebB89vzY2np6iHbDmssodhOR94FvcwLG9aax83xowC0qy1k40x/wIG4O4FcDdwq7X2z5OtVwWoiH/ampnNa9PX8f6cDeQWFHF+u7rcfVZzmtdWISrVQzkK0J7Ar8Bi4NBNUw8CDQGstS8XF6njgP64h2G5wVp7wkZQ7aSIOG1fTj5/n7SQKcu2c3nnZB4f2E6dFwpQyQVoZVHDKuLfdu3P5c0Z6xg/Yz0H8wsZ0KEed/ZrTtPEKKejiVSqshaglUXtpIj4g6Iiy7PTVjF22io6N6rJK9d2JiEqzOlY4rCS2koN6CUi5RYfFcb/nduKX+/ry829mjJl6XbOfuYX/jZpIRm/v+sexHpknPt50SSn44qIiEglCAoy3Ht2C164qhNLt2Ry8bgZLNtSrn7WpBpQASoiFVarRij3n9eKX+/rw42nN4bFk6jx/b3F44dZ9/OXd6oIFRERCWAXpNbl41tOo7DIctlLv/Pdkm1ORxI/pAJURLwmISqMhy9sw3/iPifS5B09Mz/b3dGJiIiIBKx29WOZfMfptKwTzS0T5vHKL2vw51v+xPdUgIqI1wVneR45wmamczCvwMdpRERExJeSYsL5cHgPLkyty7++/ZORk5dSWKQiVNxUgIqI95UwiP3monh6P/Uz78/eSEFhkcdlREREpOoLdwUzdnBHhvdqwtszN3DrhHnk5Bc6HUv8gApQEfG+fiPcg9kfyRVB3pkP0aBWJA9+tphzn53O90u36bIcERGRABUUZHjw/NY8clEbpi7fzpWvzWL/3PfVSWE1F+J0ABEJQIcGsT9mcPsmqYP4uK9lyrLt/Pu7P7n53Xl0S6nFiIva0K5+rLOZRUREpFLccHpj6saGM3XiOEK+fg3Idc841Ekh/O+3gwQ8jQMqIo4oKCxiYtomnpmykt0H8/hLlwb87ZyWJEZr3DDxfxoHVESk7HKfakPYAQ/9RMQ2gHuW+D6QVCqNAyoifiUkOIiruzfix7/3Zujpjfl4Xjp9//szr05fQ16B7g8VEREJNGEHtniekZnu2yDiKBWgIuKo2AgXD1/Yhu/v6UXXxrV44ps/OffZ6Uxbvl33h4qIiASSEjopLHG6BCQVoCLiF5omRvHm9V1564auGAND307jurfmsm7nAaejiYiIiDd46KTwoA1lcau7HAokTlABKiJ+pU/LJL6/uxcPX9CaPzbs4dxnpzNm6kp31+2LJqnnPBERkaoqdRBcNNZ9zyeGwphkXoy+k0t+rc8XCzyPIS6BR73giojfcQUHcdMZTRjQoR6Pfb2c56atYv/c93mw8CWCC3PcC6nnPBERkaonddDhdjsYuDknn7lvp3H3xAUUFFou66zLcQOdzoCKiN9Kignn+Ss78u7QbgzLn/C/4vOQ/Gz3UC8iIiJSJUWHu3j7xm6c3jSBv3+8kE/mqUOiQKcCVET83hnNE6ltd3qeqZ7zREREqrRwVzCvDemiIrSaUAEqIlWCKaGHvPyoej5OIiIiIt4WEeouQk9rGs/fP17Ip/NVhAYqFaAiUjV46Dkvm1D+sfcSnpm6UmOHioiIVHERocG8PqQrpzWN528fqQgNVCpARaRqOKbnPGIbUHjBc9B+EGOnreLC539lwaa9TqcUERGRCji2CP38D/WOG2jUC66IVB1H9JwHEAWM6QoXdajLQ58t4dIXZ3Dj6Y352zktiQgNdi6niIiIlNuhIvSG8XP420cLqREWwtltajsdS7xEZ0BFpMrr26o2U+7pxVXdG/L6b+s499npzF67y+lYIiIiUk4RocG8fl1X2tWP5fb35zNjdQmdEUqVowJURAJCdLiL0Ze058PhPTAGBr82i8e/XkZOfqHT0UQwxrxpjNlhjFlSwvxYY8yXxpiFxpilxpgbfJ1RRMTfRIWF8PYNXWkcX4Nh76Qxf+MepyOJF6gAFZGA0qNJPN/ceQZXd2/Ia7+u46Lnf2PJ5kynY4mMB/qfYP7twDJrbQegN/C0MSbUB7lERPxaXGQo7w7tRmJ0GNe/OYflW/c5HUkqSAWoiAScGmEhjL6kPW/f2I19Oflc8sIMnp+2ioJC9ZQrzrDWTgd2n2gRINoYY3Df3rwbKPBFNhERf5cUE86Eod2pERbCtW/MYd3OA05HkgpQASoiAevMFolMuftMLkity9NTV3LZyzNZk7Hf6VginowDWgNbgMXAXdZaj0dMjDHDjTFpxpi0jIwMX2YUEXFMg1qRvDu0O0XWcs3rs9mxL8fpSFJOKkBFJKDFRrp4bnBHxl3VkQ27DnDB2F+ZMGsD1lqno4kc6VxgAVAPOAUYZ4yJ8bSgtfZVa20Xa22XxMRE3yUUEXFYs6Qo3r6hG3sO5nH9W3PJysl3OpKUgwpQEakWLkytx5S7e9E1pRYPf76EWyfM50Da+zCmHYyMcz8vmuR0TKm+bgA+tW6rgXVAK4cziYj4nfbJsbx4dSdWbM/itvfmk6/ba6ocFaAiUm0kxYTz9g3deOj81kSs+ISgr+6GzE2AdT9/eaeKUHHKRqAfgDGmNtASWOtoIhERP9W7ZRL/urQ9v67aycQ3nsbqYHKV4pUC1BjT3xizwhiz2hhzv4f5YcaYicXzZxtjUryxXRGRsgoKMgzr1YR/x31OBLlHz8zPhmmjnAkmAc0Y8wEwE2hpjEk3xgw1xtxijLmleJHHgNOMMYuBacB91loNeiciUoJBXRrwUuoaLt38H4wOJlcpIRVdgTEmGHgBOBtIB+YaYyZba5cdsdhQYI+1tpkxZjDwb+AvFd22iEh5he7f4nlGZrpvg0i1YK298iTztwDn+CiOiEhA6L/9VYzJO3rioYPJqYOcCSUn5Y0zoN2A1dbatdbaPOBD4OJjlrkYeLv49cdAv+Ku5kVEnBGb7HFydmRdHwcRERGR8jAlHTTWwWS/5o0CtD6w6Yj36cXTPC5jrS0AMoF4L2xbRKR8+o0AV8RRk3II4769l/Cvb5drzFARERF/V8LB5BKni1/wu06INL6ZiPhE6iC4aCzENgAMxDYg+JKxxHS7ild+WctVGmNMRETEv3k4mJxNGPtOf9ChQFIaFb4HFNgMNDjifXLxNE/LpBtjQoBYYJenlVlrXwVeBejSpYsG6hORypM66Kh7RFzA6FOgS6NaPPDpYs4f+xtjrzyF05omOBZRRERESnCoDZ82CjLTyatRj4f3DWTDvBTe61RIWEiws/nEI2+cAZ0LNDfGNDbGhAKDgcnHLDMZuK749eXAj1ajwIuIn7qkY30m33E6sREhXPP6bF74aTVFRfrKEhER8Tupg+CeJTByL6H/t4y+V9xB2oY9PPDJYlRu+KcKF6DF93TeAXwPLAcmWWuXGmNGGWMGFC/2BhBvjFkN3AscN1SLiIg/aV47msl39OSC1Ho89f0KbnonjcyD+U7HEhERkRO4ILUu957dgk//2MxLv6xxOo544I1LcLHWfgN8c8y0EUe8zgGu8Ma2RER8pUZYCGMHn0LXlJo89tUyLnlxBq9e25nmtaOdjiYiIiIl+GvfZqzesZ//fLeC5knRnN2mttOR5Ah+1wmRiIg/McYw5NQUPhjWg6ycAga++DtTlm5zOpaIiIiUwBjDfy5PJTU5lnsmLmD1jv1OR5IjqAAVESmFLim1+PKvp9MksQbD353Hcz+s0n2hIiIifircFczL13QmLCSIm99NIytHt9H4CxWgIiKlVDc2gkk3n8qlHesz5oeV3PrePPbnFjgdS0RERDyoFxfBuKs6sX7XQf42aaEOHPsJFaAiImUQ7grm6UEd+OeFbfhh+Q4ufXEG63cecDqWiIiIeHBq03gePL81U5Zt58WfVzsdR1ABKiJSZsYYhvZszDs3dmNHVi4XvzCD31fvdDqWiIiIeHDj6Slccko9np66kp9W7HA6TrWnAlREpJxOb5bA5Nt7khQdxpA35zDz85dgTDsYGed+XjTJ6YgiIiLVnjGGf12aSus6Mdz1wR9s2KUrl5ykAlREpAIaxkfyyW2ncW+dBXT4YwRkbgKs+/nLO1WEioiI+IGI0GBeubYzxhhue28+OfmFTkeqtlSAiohUUEy4i1sL3yPS5B09Iz8bpo1yJpSIiIgcpUGtSJ4Z1IGlW/Yx+utlTseptlSAioh4gcnc7HlGZrpvg4iIiEiJ+rWuzc29mjBh1kYmL9zidJxqSQWoiIg3xCZ7nJwXVc/HQURERORE/n5uS7o0qskDnyxiTcZ+p+NUOypARUS8od8IcEUcNSmbMB7OulQ95IqIiPgRV3AQz1/VkTBXMLfrflCfUwEqIuINqYPgorEQ2wAwENuAnP7PsCDubK57a44u8xEREfEjdWMjeGZQB/7clsXIyUudjlOthDgdQEQkYKQOcj+K1QQ+6pDPsHfSuPODP9ixL4ebzmjiXD4RERE5rHfLJG7v05QXflpDjybxXNKxvtORqgWdARURqUSxES7eubEb57Wrw+ivlzP6q2UUFVmnY4mPGWPeNMbsMMYsOcEyvY0xC4wxS40xv/gyn4hIdXXPWS3omlKThz9fwsZdB52OUy2oABURqWThrmDGXdWJ605txOu/reOuiQvILdD9JtXMeKB/STONMXHAi8AAa21b4ArfxBIRqd5CgoMY85dTMAbumvgHBYVFTkcKeCpARUR8IDjIMHJAW+7r34ovF27h+jfnkpWT73Qs8RFr7XRg9wkWuQr41Fq7sXj5HT4JJiIiJNeM5ImB7flj417GTlvldJyApwJURMRHjDHc2rspzwzqwNz1u7nqtdnsPpDndCzxDy2AmsaYn40x84wxQ0pa0Bgz3BiTZoxJy8jI8GFEEZHAdVGHelzeOZmNv4wn96nWMDIOxrSDRZOcjhZwVICKiPjYpZ2SeeXazqzcnsVfXpnJtswcpyOJ80KAzsAFwLnAP40xLTwtaK191VrbxVrbJTEx0ZcZRUQC2uimy3nS9QZhB7YAFjI3wZd3qgj1MhWgIiIO6Ne6NuNv6MaWvdlc8crv6vhA0oHvrbUHrLU7gelAB4cziYhUK+G/jCac3KMn5mfDtFHOBApQKkBFRBxyatN43h/Wg6ycAi5/+XdWbMtyOpI45wugpzEmxBgTCXQHljucSUSkeslML9t0KRcVoCIiDurQII5JN58KwF9encnCTXudDSSVwhjzATATaGmMSTfGDDXG3GKMuQXAWrsc+A5YBMwBXrfWljhki4iIVILY5LJNl3JRASoi4rAWtaP5+JbTiA4P4arXZjFzzS6nI4mXWWuvtNbWtda6rLXJ1to3rLUvW2tfPmKZp6y1bay17ay1zzoYV0Skeuo3AlwRR03KNWEU9R3hUKDApAJURMQPNIyP5ONbTqNeXATXvzWH6SvVu6mIiIhPpQ6Ci8ZCbAPAcCC8Lv+XO5R39nd1OllAUQEqIuInaseEM/HmU2mcUIOb3knj5xUaClJERMSnUgfBPUtg5F4i71vOvuYDefK7P1mbsd/pZAFDBaiIiB+pVSOUD4b1oFliFMPfmcfS7153j0Om8chERER8yhjDvy9LJSwkmL99tJCCwiKnIwUEFaAiIn6mZo1Q3h/WneFxaTSe+YB7HDKNRyYiIuJztWPCGXVxW/7YuJdXpq91Ok5AUAEqIuKH4iJDuTf4QyJN3tEzNB6ZiIiITw3oUI8L2tfl2R9WsnzrPqfjVHkqQEVE/FTQvs2eZ2g8MhEREZ8xxjD6knbERrj4x8eLdCluBVWoADXG1DLGTDXGrCp+rlnCcoXGmAXFj8kV2aaISLWh8chERET8Qs0aoYy6uB2LN2fy2q/rnI5TpVX0DOj9wDRrbXNgWvF7T7KttacUPwZUcJsiItWDh/HIDtpQFre6y6FAIiIi1df57evSv20dxvywUr3iVkBFC9CLgbeLX78NXFLB9YmIyCHHjEdWFJPMSzF3cdlvyfyicUJFRER8btTFbQkPCeK+TxZRVGSdjlMlVbQArW2t3Vr8ehtQu4Tlwo0xacaYWcaYS060QmPM8OJl0zIy9ANLRKq5I8YjC7p3KTfdeh/NkqIY/k4av6/Z6XQ6ERGRaiUpJpx/XtiGuev3MGH2BqfjVEknLUCNMT8YY5Z4eFx85HLWWguUdBigkbW2C3AV8KwxpmlJ27PWvmqt7WKt7ZKYmFiWzyIiEvBiI128O7QbjeIjGTo+jbnrdzsdSUREpFq5vHMyZzRP4N/f/kn6noNOx6lyTlqAWmvPsta28/D4AthujKkLUPy8o4R1bC5+Xgv8DHT02icQEalm4qPCmHBTd+rGhnPDW3NZsGmv05FERESqDWMMTwxsjwUe/GwJ7vNwUloVvQR3MnBd8evrgC+OXcAYU9MYE1b8OgE4HVhWwe2KiFRrSdHhvD+sB7VqhDLkjdks3ZLpdCQREZFqo0GtSP5xbkumr8zg0/klDJsmHlW0AH0SONsYswo4q/g9xpguxpjXi5dpDaQZYxYCPwFPWmtVgIqIVFCd2HDeH9adqLAQrn1jDqt3qEc+ERERXxlyagodG8bx+DfL2Xswz+k4VUaFClBr7S5rbT9rbfPiS3V3F09Ps9beVPz6d2tte2tth+LnN7wRXEREILlmJO8P60GQMQx5Yzab92Y7HUlERKRaCAoyPH5JezKz8/n3d386HafKqOgZUBERcVhKQg3evrErWTkFXPvGbHbtz3U6koiISLXQpl4MN56ewgdzNpGmjgFLRQWoiEgAaFsvltev68LmPdncMH4u+3MLnI4kIiJSLdx9VgvqxYbz0GdLyC8scjqO31MBKiISILo3ieeFqzqxdMs+hr+TRm5BodORREREAl6NsBBGDmjLiu1ZvPHbOqfj+D0VoCIiAeSsNrX5z2Wp/L5mF3d9sIDCInUNLyIiUtnOaVuHs9vU5tkfVrJpt8YGPREVoCIiAeayzsn888I2fLd0Gw99tljjk/kBY8ybxpgdxpglJ1muqzGmwBhzua+yiYiId4wc0JYgYxg5eana3hNQASoiEoCG9mzMHX2a8eHcTfx3ygqn4wiMB/qfaAFjTDDwb2CKLwKJiIh31Y+L4J6zWjDtzx18v3S703H8lgpQEZEA9bdzWnBltwa88NMa3pu9wek41Zq1djpwsu4R/wp8Auyo/EQiIlIZrj89hVZ1onn0y6UczFOHgJ6oABURCVDGGB67uB19Wibyz8+XMG25jsb6K2NMfWAg8JLTWUREpPxcwUE8dkk7tmbm8NLPa5yO45dUgIqIBLCQ4CDGXdWJtvViueP9P1i4aa/TkcSzZ4H7rLUn7b/fGDPcGJNmjEnLyMio/GQiIlImXVNqcckp9Xhl+lo27lKHRMdSASoiEuBqhIXwxvVdiI8KZdKbT1PwdFsYGQdj2sGiSU7HE7cuwIfGmPXA5cCLxphLPC1orX3VWtvFWtslMTHRhxFFRKS07j+vNSFBhse+XuZ0FL+jAlREpBpIig7n49PTedi+TEhWOmAhcxN8eaeKUD9grW1srU2x1qYAHwO3WWs/dzaViIiUV53YcP7atzlTl23nl5W6WuVIKkBFRKqJOnP/QwR5R0/Mz4Zpo5wJVI0YYz4AZgItjTHpxpihxphbjDG3OJ1NREQqx409U2icUINfP3kRO0ZXHx0S4nQAERHxkcz0sk0Xr7HWXlmGZa+vxCgiIuIjYSHBjGu3msYzx2Fyiw8AH7r6CCB1kHPhHKQzoCIi1UVsctmmi4iISIW0Xf4skUZXHx1JBaiISHXRbwS4Io6adNCGMqvx7Q4FEhERCXC6+ug4KkBFRKqL1EFw0ViIbQAYbGwyb9W6h2vnNmLOut1OpxMREQk8uvroOCpARUSqk9RBcM8SGLkXc89Srhn2fzSoGcnN76axYdcBp9OJiIgEFg9XH1lXhHt6NaUCVESkGouNdPHG9V0psjD07TT25eQ7HUlERCRwHHH1kcWQXpTAgg6jqm0HRKACVESk2mucUIOXr+nM+p0HuP29+RQUFjkdSUREJHAUX31UNGIPN9V6izuXNSO3oNDpVI5RASoiIpzaNJ7Rl7Tj11U7eeyrZU7HERERCTjBQYYHz2/Npt3ZvPP7BqfjOEYFqIiIADC4W0OGndGYt2duYMKs6tswioiIVJZeLRLp3TKRsT+uYveBvJP/QQBSASoiIofdf15rerdM5NEvlzJ3vXrGFRER8bYHz2/NgdwCxk5b5XQUR6gAFRGRw4KDDM8N7khyzUhunTCfrZnZTkcSEREJKC1qRzO4W0MmzNrA2oz9TsfxORWgIiJylNgIF69e25nsvAJueXceOfnVt6MEERGRynDPWS0ICwniX9/+6XQUn1MBKiIix2leO5pn/nIKC9MzeeizJVhrnY4kIiISMBKjw7itTzOmLtvOzDW7nI7jUypARUTEo3Pb1uGufs35ZH46b/++3uk4IiIiAWVoz8bUiw3niW+WV6sDvSpARUSkRHf1a85ZrWvz2NfLq90RWhERkcoU7grmnrNbsHhzJt8u2eZ0HJ+pUAFqjLnCGLPUGFNkjOlyguX6G2NWGGNWG2Pur8g2RUTEd4KCDGP+0oGU+Ehuf38+6XsOOh1JREQkYFzaKZnmSVH89/sVFBQWOR3HJyp6BnQJcCkwvaQFjDHBwAvAeUAb4EpjTJsKbldERHwkOtzFa0O6kF9QxG3vzSe3QJ0SiYiIeENwkOHv57Zk7c4DfDwv3ek4PlGhAtRau9xau+Iki3UDVltr11pr84APgYsrsl0REfGtJolR/HdQBxalZzLqy2VOxxEREQkY57SpTceGcTz7w6pq0fO8L+4BrQ9sOuJ9evE0j4wxw40xacaYtIyMjEoPJyIipXNu2zrc3KsJ783eyGd/VI+jtCIiIpXNGMN9/VuxbV9Otej076QFqDHmB2PMEg+PSjmLaa191VrbxVrbJTExsTI2ISIi5fR/57akW+NaPPDpYv7cts/pOCIiIgGhR5N4zmyRyIs/ryEzO9/pOJXqpAWotfYsa207D48vSrmNzUCDI94nF08TEZEqJiQ4iHFXdiQ63MWtE+aTlRPYjaSIiIiv/N+5LcnMzufV6WucjlKpfHEJ7lyguTGmsTEmFBgMTPbBdkVEpBIkxYQz7sqObNx9kH98vKhajV1WXsaYN40xO4wxS0qYf7UxZpExZrEx5ndjTAdfZxQREWe1qx/LgA71ePO39ezYl+N0nEpT0WFYBhpj0oFTga+NMd8XT69njPkGwFpbANwBfA8sByZZa5dWLLaIiDipe5N4/nFuS1zLPmb/v1vDyDgY0w4WTXI6mr8aD/Q/wfx1wJnW2vbAY8CrvgglIiL+5d6zW5BfWMTzP652OkqlCanIH1trPwM+8zB9C3D+Ee+/Ab6pyLZERMS/DK85j7ywNwjLyXVPyNwEX97pfp06yLlgfshaO90Yk3KC+b8f8XYW7ttVRESkmklJqMHgbg34YM5Ghp3RhIbxkU5H8jpfXIIrIiIByEwbRZjNPXpifjZMG+VMoMAxFPi2pJnqLV5EJLD9tW9zgoIMz/+4yukolUIFqIiIlE9mCUOxlDRdTsoY0wd3AXpfScuot3gRkcBWOyaca7o34tM/NrNu5wGn43idClARESmf2BKuEi1pupyQMSYVeB242Fq7y+k8IiLinFt6N8EVbHh+WuCdBVUBKiIi5dNvBLgijpp00IaysePfHQpUdRljGgKfAtdaa1c6nUdERJyVFB3OkFNT+HzBZtZk7Hc6jlepABURkfJJHQQXjYXYBoChMCaZf4fcxrVzG2l80GMYYz4AZgItjTHpxpihxphbjDG3FC8yAogHXjTGLDDGpDkWVkRE/MLNvZoQ7gpmbICdBa1QL7giIlLNpQ463ONtMHDR+t1MeHUWD362hLGDT8EY42w+P2GtvfIk828CbvJRHBERqQLio8K47rQUXv5lDXf0aUbz2tFOR/IKnQEVERGv6ZJSi3vPbsGXC7cwce4mp+OIiIhUacPPaEKkK5hnA+gsqApQERHxqlvPbErPZgmM/HIpK7dnOR1HRESkyqpZI5Qbezbm60Vb+XPbPqfjeIUKUBER8aqgIMMzf+lAVFgId7w/n5z8QqcjiYiIVFk39WxCdFgIz04NjLOgKkBFRMTrkqLDeWbQKazcvp8nv/3T6TgiIiJVVmykixt7Nua7pdtYsjnT6TgVpgJUREQqRa8Widxwegrjf1/PLysznI4jIiJSZd3YszEx4SGM+3G101EqTAWoiIhUmvv6t6JF7Sj+/tFCdu3PdTqOiIhIlRQb4eL601L4bum2Kt+/ggpQERGpNOGuYJ4b3JHMg/nc/+lirLVORxIREamSbji9MZGhwbzwU9U+C6oCVEREKlXrujH8o39Lpi7bzocamkVERKRcatYI5doejfhy4RbW7zzgdJxyUwEqIiKV7sbTG9OzWQKjvlzG2oz9TscRERGpkoae0ZiBIb8T+0pHGBkHY9rBoklOxyoTFaAiIlLpgoIM/72iA2GuIO6euID8wiKnI4mIiFQ5Sesm86+Q16iZvx2wkLkJvryzShWhKkBFRMQn6sSG86+B7VmUnsmzP6x0Oo6IiEjVM20UofaYTv3ys2HaKGfylIMKUBER8Znz2tflis7JvPTzGuZt2ON0HBERkaolM71s0/2QClAREfGpERe1oW5sBH//aCHZeYVOxxEREak6YpPLNt0PqQAVERGfig538dQVqazbeYB/f/en03FERESqjn4jwBVx1CQbEuGeXkWoABUREZ87rWkC15+Wwvjf1/P7mp1OxxEREakaUgfBRWMhtgEWQ3pRAt80fsA9vYpQASoiIo74R/+WpMRH8o+PF7E/t8DpOCIiIlVD6iC4Zwlm5F5GN5/E/ataVal2VAWoiIg4IjI0hKcHdWDL3mwe/3q503FERESqnFt7NyUrp4AP52x0OkqpqQAVERHHdG5Ui2FnNOGDORv5ecUOp+OIiIhUKR0axNGjSS3e+G0deQVVY4xtFaAiIuKoe85uQfOkKO7/ZDGZB/OdjiMiIlKl3HxmU7Zm5vDlwi1ORykVFaAiIuKocFcwTw/qQMb+XB79cqnTcURERKqU3i0SaVk7mlemr8Fa63Sck1IBKiIijktNjuP23k359I/NTF223ek4XmeMedMYs8MYs6SE+cYYM9YYs9oYs8gY08nXGUVEpGoyxnDzmU1YuX0/P6/IcDrOSVWoADXGXGGMWWqMKTLGdDnBcuuNMYuNMQuMMWkV2aaIiASmO/o2p1WdaKZ/8gJFz7SFkXEwph0smuR0NG8YD/Q/wfzzgObFj+HASz7IJCIiAeKiDvWoFxvOy7+scTrKSVX0DOgS4FJgeimW7WOtPcVaW2KhKiIi1VdoSBCvdFjLAwUvEbQvHbCQuQm+vLPKF6HW2unA7hMscjHwjnWbBcQZY+r6Jp2IiFR1ruAgbuzZmNnrdvPHxj1OxzmhChWg1trl1toV3gojIiLVW6MF/yXS5B09MT8bpo1yJpDv1Ac2HfE+vXjacYwxw40xacaYtIwM/7/USkREfGNwt4bEhIfw6vS1Tkc5IV/dA2qBKcaYecaY4SdaUA2riEg1lpletunVkLX2VWttF2ttl8TERKfjiIiIn4gKC+HaUxvx3dJtrNt5wOk4JTppAWqM+cEYs8TD4+IybKentbYT7ntcbjfG9CppQTWsIiLVWGxy2aYHjs1AgyPeJxdPExERKbXrTkvBFRzEa7/671nQkxag1tqzrLXtPDy+KO1GrLWbi593AJ8B3cofWUREAla/EeCKOGpSUXCEe3pgmwwMKe4NtweQaa3d6nQoERGpWpKiw7msUzIfz0snIyvX6TgeVfoluMaYGsaY6EOvgXNwd14kIiJytNRBcNFYiG2AxbCVRP4deit5bS53OlmFGGM+AGYCLY0x6caYocaYW4wxtxQv8g2wFlgNvAbc5lBUERGp4oad0Zj8wiLenbXB6SgehVTkj40xA4HngUTga2PMAmvtucaYesDr1trzgdrAZ8aYQ9t731r7XQVzi4hIoEodBKmDMMCy5dt55e00In5ezd1ntXA6WblZa688yXwL3O6jOCIiEsCaJEbRr1US783awG29mxLuCnY60lEq2gvuZ9baZGttmLW2trX23OLpW4qLT6y1a621HYofba21j3sjuIiIBL5+rWtz8Sn1eOGn1azYluV0HBERkSrhxtMbs+tAHl8s8L/uBHzVC66IiEi5jLiwDVFhITzw6SKKiqzTcURERPzeqU3jaVUnmjd/W4/7Ihv/oQJURET8WnxUGA+e35r5G/fy4dxNJ/8DERGRas4Yw9CejVmxPYsZq3c5HecoKkBFRMTvXd45me6Na/Hkt8v9tlc/ERERf3JRh3okRIXyxm/+NSSLClAREfF7xhgeH9ie7PxCRn+9zOk4IiIifi/cFcw1PRrx04oMVu/Y73Scw1SAiohIldAsKYpbezfjiwVb+HVVhtNxRERE/N41PRoRGhLE+N/XOR3lMBWgIiJSZdzWuymNE2rw8OdLyMkvdDqOiIiIX0uICuOSU+rxybzN7D2Y53QcQAWoiIhUIeGuYEZf0o4Nuw7ywk+rnY4jIiLi927s2Zjs/ELen7PR6SiAClAREaliTm+WwMCO9Xn5lzWs3qGxQUVERE6kVZ0YTm8Wzzu/byC/sMjpOCpARUSk6nnogtZEhobw4GdL/G58MxEREX8ztGdjtu3L4ZvFW52OogJURESqnoSoMB44rxVz1u3mo3npTscRERHxa71bJNEkoQZvzljvdBQVoCIiUjUN6tKALo1q8uS3f5J5MN/pOCIiIn4rKMgw5NRGLNy0l0Xpe53N4ujWRUREyikoyDDq4nbsPZjH01NXOB1HRETEr13aOZnI0GDembnB0RwqQEVEpMpqUy+Ga3s0YsKsDSzdkul0HBEREb8VE+7i0k71mbxwC7sPODckiwpQERGp0u49pyU1I0MZ8cVSiorUIZGIiEhJhpyaQl5BEZPSNjmWQQWoiIhUabERLu47rxXzNuzh0z82Ox1HRETEb7WoHU2PJrV4d+YGCh06aKsCVEREqrzLOyXTsWEcT367nMxsdUgkIiJSkutOTWHz3mx++nOHI9tXASoiIlVeUJBh1IB27DqQx7M/rHQ6joiIiN86u01t6sSE8/bM9Y5sXwWoiIgEhPbJsVzVrSHvzNzAn9v2OR1HRETEL4UEB3F194b8umonazP2+3z7KkBFRCRg/N+5LYkJD2HE50uxVh0SiYiIeDK4W0NcwYZ3Z/l+SBYVoCIiEjDiIkP5R/9WzFm/my8WbHE6zlGMMf2NMSuMMauNMfd7mN/QGPOTMeYPY8wiY8z5TuQUEZHAlxgdxvnt6/JxWjoHcgt8um0VoCIiElD+0qUBHZJjeeKb5T5vVEtijAkGXgDOA9oAVxpj2hyz2MPAJGttR2Aw8KJvU4qISHUy5NQUsnIL+HyBb3uQVwEqIiIBJSjI8MiAtuzIyuXFn1c7HeeQbsBqa+1aa20e8CFw8THLWCCm+HUs4F+ncEVEJKB0ahhH23oxvPP7Bp/etqICVEREAk6nhjW5tGN9Xvt1HRt3HXQ6DkB94MhRv9OLpx1pJHCNMSYd+Ab4q6cVGWOGG2PSjDFpGRkZlZFVRESqAWMM152awortWcxdv8dn21UBKiIiAekf/VsREmR44pvlTkcprSuB8dbaZOB84F1jzHHttLX2VWttF2ttl8TERJ+HFBGRwHFRh3pEh4fw3mzfdUakAlRERAJSndhwbuvdlNDlH5PzVGsYGQdj2sGiSU7E2Qw0OOJ9cvG0Iw0FJgFYa2cC4UCCT9KJiEi1FBEazMhGS/nH8suxPmonQyp17SIiIg4aXnMehaFvEH4g1z0hcxN8eaf7deogX0aZCzQ3xjTGXXgOBq46ZpmNQD9gvDGmNe4CVNfYiohI5Vk0iYGb/0OQyXa/90E7qTOgIiISsEJ/Hk0EuUdPzM+GaaN8msNaWwDcAXwPLMfd2+1SY8woY8yA4sX+BgwzxiwEPgCutxrMVEREKtO0UQQVZB89rZLbyQqdATXGPAVcBOQBa4AbrLV7PSzXH3gOCAZet9Y+WZHtioiIlEpmetmmVyJr7Te4Oxc6ctqII14vA073dS4REanGHGgnK3oGdCrQzlqbCqwEHjh2gVKOfSYiIuJ9scllmy4iIlKdONBOVqgAtdZOKb6sCGAW7k4VjlWasc9ERES8r98IcEUcPc0V4Z4uIiJS3TnQTnrzHtAbgW89TC/N2GciIiLelzoILhoLsQ0A436+aKyvOyASERHxTw60kye9B9QY8wNQx8Osh6y1XxQv8xBQALxX0UDGmOHAcICGDRtWdHUiIlLdpQ5SwSkiIlISH7eTJy1ArbVnnWi+MeZ64EKgXwm99ZVm7LMjt/cq8CpAly5d1PufiIiIiIhIgKjQJbjFvdv+AxhgrT1YwmKHxz4zxoTiHvtsckW2KyIiIiIiIlVPRe8BHQdEA1ONMQuMMS8DGGPqGWO+gZLHPqvgdkVERERERKSKqdA4oNbaZiVM3wKcf8T748Y+ExERERERkerFm73gioiIiIiIiJRIBaiIiIiIiIj4hPHcca1/MMZkABu8tLoEYKeX1lVdaJ+VnfZZ2WmflY/2W9l5c581stYmemld5aZ20nHaZ+Wj/VZ22mdlp31Wdt7eZx7bSr8uQL3JGJNmre3idI6qRPus7LTPyk77rHy038pO++zEtH/KTvusfLTfyk77rOy0z8rOV/tMl+CKiIiIiIiIT6gAFREREREREZ+oTgXoq04HqIK0z8pO+6zstM/KR/ut7LTPTkz7p+y0z8pH+63stM/KTvus7Hyyz6rNPaAiIiIiIiLirOp0BlREREREREQcpAJUREREREREfCLgClBjTH9jzApjzGpjzP0e5ocZYyYWz59tjElxIKZfKcU+u9cYs8wYs8gYM80Y08iJnP7kZPvsiOUuM8ZYY0y17wa8NPvMGDOo+N/aUmPM+77O6G9K8f9mQ2PMT8aYP4r//zzfiZz+xBjzpjFmhzFmSQnzjTFmbPE+XWSM6eTrjE5TO1l2aifLTu1k+aitLDu1lWXjF+2ktTZgHkAwsAZoAoQCC4E2xyxzG/By8evBwESnc1eBfdYHiCx+fav22cn3WfFy0cB0YBbQxenc/r7PgObAH0DN4vdJTueuAvvsVeDW4tdtgPVO53b6AfQCOgFLSph/PvAtYIAewGynM/t4/6idrJx9pnayjPuseDm1k2Xcb2ory7XP1FYevT8cbycD7QxoN2C1tXattTYP+BC4+JhlLgbeLn79MdDPGGN8mNHfnHSfWWt/stYeLH47C0j2cUZ/U5p/ZwCPAf8GcnwZzk+VZp8NA16w1u4BsNbu8HFGf1OafWaBmOLXscAWH+bzS9ba6cDuEyxyMfCOdZsFxBlj6vomnV9QO1l2aifLTu1k+aitLDu1lWXkD+1koBWg9YFNR7xPL57mcRlrbQGQCcT7JJ1/Ks0+O9JQ3EdFqrOT7rPiyxUaWGu/9mUwP1aaf2ctgBbGmBnGmFnGmP4+S+efSrPPRgLXGGPSgW+Av/omWpVW1u+8QKN2suzUTpad2snyUVtZdmorva/S28kQb65MApsx5hqgC3Cm01n8mTEmCHgGuN7hKFVNCO5Li3rjPnsw3RjT3lq718lQfu5KYLy19mljzKnAu8aYdtbaIqeDiVRHaidLR+1khaitLDu1lX4m0M6AbgYaHPE+uXiax2WMMSG4T8Xv8kk6/1SafYYx5izgIWCAtTbXR9n81cn2WTTQDvjZGLMe9/Xzk6t5Bwul+XeWDky21uZba9cBK3E3stVVafbZUGASgLV2JhAOJPgkXdVVqu+8AKZ2suzUTpad2snyUVtZdmorva/S28lAK0DnAs2NMY2NMaG4O0+YfMwyk4Hril9fDvxoi++4raZOus+MMR2BV3A3qtX9XgM4yT6z1mZaaxOstSnW2hTc9wMNsNamORPXL5Tm/83PcR/RxRiTgPsyo7U+zOhvSrPPNgL9AIwxrXE3qhk+TVn1TAaGFPfy1wPItNZudTqUD6mdLDu1k2WndrJ81FaWndpK76v0djKgLsG11hYYY+4AvsfdK9ab1tqlxphRQJq1djLwBu5T76tx34A72LnEzivlPnsKiAI+Ku6HYqO1doBjoR1Wyn0mRyjlPvseOMcYswwoBP7PWlttz7qUcp/9DXjNGHMP7k4Wrq/mhQLGmA9w/zhLKL7f5xHABWCtfRn3/T/nA6uBg8ANziR1htrJslM7WXZqJ8tHbWXZqa0sO39oJ0013v8iIiIiIiLiQ4F2Ca6IiIiIiIj4KRWgIiIiIiIi4hMqQEVERERERMQnVICKiIiIiIiIT6gAFREREREREZ9QASoS4IwxccaY25zOISIi4q/UVor4jgpQkcAXB6hRFRERKVkcaitFfEIFqEjgexJoaoxZYIx5yukwIiIifkhtpYiPGGut0xlEpBIZY1KAr6y17ZzOIiIi4o/UVor4js6AioiIiIiIiE+oABURERERERGfUAEqEviygGinQ4iIiPgxtZUiPqICVCTAWWt3ATOMMUvUsYKIiMjx1FaK+I46IRIRERERERGf0BlQERERERER8QkVoCIiIiIiIuITKkBFRERERETEJ1SAioiIiIiIiE+oABURERERERGfUAEqIiIiIiIiPqECVERERERERHzi/wGMJFmGva32iwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 936x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "t_fine = np.linspace(0, t[-1], 1000)\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4))\n",
    "ax1.plot(t_fine, 2 * np.exp(t_fine) - np.exp(2 * t_fine), t_sol, x(t_sol), \"o\")\n",
    "ax1.set_xlabel(\"t\")\n",
    "ax1.legend([\"2*exp(t) - exp(2*t)\", \"x\"], loc=\"best\")\n",
    "\n",
    "ax2.plot(t_fine, 3 * np.exp(t_fine) - np.exp(2 * t_fine), t_sol, y(t_sol), \"o\")\n",
    "ax2.set_xlabel(\"t\")\n",
    "ax2.legend([\"3*exp(t) - exp(2*t)\", \"y\"], loc=\"best\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the [next notebook](./2-a-pde-model.ipynb) we show how to create, discretise and solve a PDE model in pybamm."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "The relevant papers for this notebook are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "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). ECSarXiv. February, 2020. doi:10.1149/osf.io/67ckj.\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": "Python 3 (ipykernel)",
   "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.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}