physics-based-models-1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Physics-based models (1): what are they?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The exponential growth of physics-based models articles\n",
    "\n",
    "Let's fit an exponential curve to the data on number of articles about battery physics-based models per year. The data, which we obtain with a manual search in Google Scholar with the prompt `\"physics-based\" model \"lithium\" battery` is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "years = np.arange(2000, 2023)\n",
    "articles = np.array([43, 56, 65, 54, 59, 92, 85, 99, 109, 119, 176, 215, 262, 315, 366, 459, 510, 617, 770, 866, 1090, 1520, 1900])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now use the `curve_fit` method from `scipy.optimize` to fit an exponential curve to the data. First we need to define the curve we want to fit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def curve(x, a, b):\n",
    "    return a * np.exp(b * (x - 2000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And note we shifted the data by `2000` so the problem is better conditioned. Now we can use `curve_fit` to determine the value of the parameters `a` and `b`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fit: 18.02 * exp(0.21 * (x - 2000))\n"
     ]
    }
   ],
   "source": [
    "from scipy.optimize import curve_fit\n",
    "\n",
    "parameters, _ = curve_fit(curve, years, articles)\n",
    "fit = curve(years, *parameters)\n",
    "print(f'Fit: {parameters[0]:.2f} * exp({parameters[1]:.2f} * (x - 2000))')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We finally plot the data and the fit, to produce the figure shown in the blog post."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABoOklEQVR4nO3deXhM1/8H8PfMZJdNyEoQW4h9KaK0NCQUXWiraqt9bYpW0QXVxR4aS+lm+VZbtVRrF4ldqKZiCZIgxBaxJBMh68z5/aG5v4wkZJjJnZm8X88zT829Z+58rqvxdu495yiEEAJEREREZPaUchdARERERIbBYEdERERkIRjsiIiIiCwEgx0RERGRhWCwIyIiIrIQDHZEREREFoLBjoiIiMhCMNgRERERWQgruQswB1qtFtevX4eTkxMUCoXc5RAREVE5IoTAvXv34OPjA6Xy8X1yDHalcP36dfj6+spdBhEREZVjV65cQdWqVR/bhsGuFJycnAA8/A11dnaWuRoiIiIqTzIyMuDr6yvlkcdhsCuFgtuvzs7ODHZEREQki9I8DsbBE0REREQWgsGOiIiIyEIw2BERERFZCD5jZ0AajQZ5eXlyl0FEFsTa2hoqlUruMojITDDYGYAQAikpKUhPT5e7FCKyQK6urvDy8uI8mkT0RAx2BlAQ6jw8PODg4MAfvkRkEEIIPHjwAKmpqQAAb29vmSsiIlPHYPeMNBqNFOoqVaokdzlEZGHs7e0BAKmpqfDw8OBtWSJ6LA6eeEYFz9Q5ODjIXAkRWaqCny98hpeInoTBzkB4+5WIjIU/X4iotBjsqFgdOnTAuHHj5C6D/pObm4vatWvj8OHDcpdiEJZ2PqVx5swZVK1aFffv39fZ3qZNG2zYsEGmqojI0jDYlWPvvvsuFApFkdf58+exceNGfPHFF1LbGjVqYOHChXp/R3Z2Nt599100atQIVlZWeO2114ptt2bNGjRp0gQODg7w9vbG4MGDcefOncceOzk5Gd26dYODgwM8PDwwceJE5OfnS/s3btyIzp07w93dHc7OzggMDMTOnTtLVff06dOxd+/e0p6m0S1btgx+fn5o27attO3u3bvo27cvnJ2d4erqiiFDhiAzM7PEY9y9exfvvfce/P39YW9vj2rVqiE0NBRqtVqnXWhoKFq0aAFbW1s0bdq01DVeunQJ77777lOfjzGcOHECffr0ga+vL+zt7VG/fn188803Rdrt3bsXzZs3h62tLWrXro2VK1cWabNkyRLUqFEDdnZ2aN26Nf7++2+d/dnZ2RgzZgwqVaoER0dH9OrVCzdv3pT2BwQEoE2bNggLC9P53KefforJkydDq9Ua5qSJqFxjsCvnunTpghs3bui8/Pz84ObmVqrFhp9Eo9HA3t4eoaGh6NSpU7FtDh06hAEDBmDIkCGIi4vDunXr8Pfff2PYsGGPPW63bt2Qm5uLw4cPY9WqVVi5ciWmTp0qtdm/fz86d+6Mbdu2ISYmBh07dkSPHj1w/PjxYo+Zl5eH+fPn6zzHlJqaiuXLlz/l2RuGEAKLFy/GkCFDdLb37dsXcXFxiIiIwJYtW7B//34MHz68xONcv34d169fx7x583D69GmsXLkSO3bsKHJcABg8eDB69+5dqvrWrFmDCxcu6NS7ZMkSpKWl6XU+xhATEwMPDw/8/PPPiIuLwyeffIIpU6Zg8eLFUpukpCR069YNHTt2RGxsLMaNG4ehQ4fq/CNg7dq1mDBhAqZNm4Z///0XTZo0QUhIiDRaFQDGjx+PzZs3Y926ddi3bx+uX7+Onj176tQzaNAgfPvttzr/AOnatSvu3buH7du3G/F3gojKDUFPpFarBQChVquL7MvKyhJnzpwRWVlZMlT2bAYOHCheffXVYve9+OKL4v3335d+DUDnZcjvmzt3rqhZs6bOtvDwcFGlSpUSj7Vt2zahVCpFSkqKtO3bb78Vzs7OIicnp8TPBQQEiM8//7zYfXl5eWLevHmiXbt24o033hChoaGiVatWYsuWLSUeLy0tTQwZMkRUrlxZODk5iY4dO4rY2FghhBCpqanC09NTfPXVV1L7Q4cOCWtra7F7924hhBDTpk0TTZo0EcuWLRNVq1YV9vb24s033xTp6enSZ44dOyaUSqXIyMiQtp05c0YAEMeOHZO2bd++XSgUCnHt2rUS633U77//LmxsbEReXl6RfQW1PcmhQ4fE888/LyZNmiSCgoJEcHCwmDJlirh//36x7Ys7n1WrVokKFSqIhIQEaduoUaOEv79/icd5WqNHjxYdO3aU3n/00UeiQYMGOm169+4tQkJCpPetWrUSY8aMkd5rNBrh4+MjZs6cKYQQIj09XVhbW4t169ZJbc6ePSsAiOjoaGlbTk6OsLW1la5/gUGDBol+/fqVWLM5/5whomf3uBzyKPbY0RNt3LgRVatWxYwZM6RevQIKhaLY21b6CAwMxJUrV7Bt2zYIIXDz5k2sX78eL7/8comfiY6ORqNGjeDp6SltCwkJQUZGBuLi4or9jFarxb179+Dm5lbsfisrK3zwwQcIDw/Htm3bsGvXLuzatQvdunUrsY4333wTqamp2L59O2JiYtC8eXMEBQXh7t27cHd3x08//YTp06fjn3/+wb1799C/f3+MHTsWQUFB0jHOnz+P33//HZs3b8aOHTtw/PhxjB49Wtp/4MAB1K1bV6cHNTo6Gq6urmjZsqW0rVOnTlAqlTh69GiJ9T5KrVbD2dkZVlZPP/NR27ZtsWfPHkRHR2Pv3r0YN24cvv766xJHihd3PgMGDMDLL7+Mvn37Ij8/H1u3bsUPP/yANWvWGHzEuVqt1vkzEB0dXaQ3OSQkBNHR0QAePg8YExOj00apVKJTp05Sm5iYGOTl5em0qVevHqpVqya1AQAbGxs0bdoUBw4c0Pm+Vq1aFdlGRPQ0OI+dET36kHRhKpUKdnZ2pWqrVCqluawe17ZChQp617hlyxY4OjpK77t27Yp169bptHFzc4NKpYKTkxO8vLx09vn7+8PFxUXv7y3s+eefx5o1a9C7d29kZ2cjPz8fPXr0wJIlS0r8TEpKik6oAyC9T0lJKfYz8+bNQ2ZmJt56661i92s0GixevFgKlT4+PujSpQumTZuGLl26FGl/8OBB/P3330hNTYWtra30HZs2bcL69esxfPhwvPzyyxg2bBj69u2Lli1bokKFCpg5c6bOcbKzs7F69WpUqVIFALBo0SJ069YN8+fPh5eXFy5fvgwfH58i5+/h4aGzzcrKCm5ubiWe/6Nu376NL7744rG3b0vj6NGjmDhxItq2bQtra2ssXLgQ0dHR+Pjjj3X+jBco7nwAYPny5WjcuDFCQ0OxceNGTJ8+HS1atHim2h51+PBhrF27Flu3bpW2lfRnKSMjA1lZWUhLS4NGoym2zblz56Rj2NjYwNXVtUibR6+Hj48PLl++XGTblStXoNVqoVTy39tE9PT4E8SIHB0dS3z16tVLp62Hh0eJbbt27arTtkaNGsW2exoFzxUVvMLDw/X6/Llz5/D6668/1XcXOHPmDN5//31MnToVMTEx2LFjBy5duoSRI0c+03EL++WXX/D555/j999/LxKICmi1WuTl5SEyMhINGjTA66+/jk2bNuHSpUvFtj9x4gQyMzOlh+ULXklJSTrPnM2bNw/5+flYt24d1qxZI4XAAtWqVZNCHfCwB1Or1SI+Ph4AkJWVVWxAehYZGRno1q0bAgICMH369Gc6VmJiIlasWIGRI0eiatWq2LFjBzw9PfHgwYNi25d0PhUrVsSPP/6Ib7/9FrVq1cLkyZNL/M7k5GSd3/Ovv/76iXWePn0ar776KqZNm4bg4ODSn6CB2dvbF/m9sbe3h1arRU5OjkxVEZGlYI9dOVehQgXUrl1b1hpmzpyJ559/HhMnTgQANG7cGBUqVED79u3x5ZdfFruMkpeXV5FRiQUjEB/tVfztt98wdOhQrFu3rsQBHMDDxdY//PBDnW2enp4lBszMzEx4e3sXO3q2cM/NhQsXcP36dWi1Wly6dAmNGjUqsYbiVK5cGadOndLZ5uXlpfPgPgDk5+fj7t27Rc7/Uffu3UOXLl3g5OSEP/74A9bW1nrV86h+/foBgBSAFQoFxowZU2L74s6nwP79+6FSqXDjxg3cv3+/xAE8Pj4+iI2Nld6XdHu9wJkzZxAUFIThw4fj008/1dnn5eWlM3oVePhnydnZGfb29lCpVFCpVMW2Kfi99vLyQm5uLtLT03WufeE2Be7evYtatWoV2VahQgWdnnkiMn1paWmwsrIyyGBDQ2GPnRFlZmaW+Hp03qrU1NQS2z46Wu7SpUvFtjMmGxsbaDQaoxz7wYMHRW4/FSybJIQo9jOBgYE4deqUTriJiIiAs7MzAgICpG2//vorBg0ahF9//fWxz8o9avr06ejQocNj2zRv3hwpKSmwsrJC7dq1dV6VK1cG8PD5rH79+qF379744osvMHTo0CKBLDk5GdevX5feHzlyBEqlEv7+/gCAZs2a4dy5czq/F4GBgUhPT0dMTIy0LSoqClqtFq1bty6x5oyMDAQHB8PGxgZ//fWXQXsCa9SoUarnLYs7H+DhbdLZs2dj8+bNcHR0xNixY0s8xqO/548LdnFxcejYsSMGDhyIr776qsj+wMBAREZG6myLiIhAYGAggId/9lu0aKHTRqvVIjIyUmrTokULWFtb67SJj49HcnKy1KbA6dOn0axZsyduIyLTt3TpUlSqVElnRgbZGXkgh0Uo76NihRCic+fO4pVXXhFXr14Vt27dkrb7+/uLjRs3PvZ74uLixPHjx0WPHj1Ehw4dxPHjx8Xx48el/StWrBBWVlZi6dKl4sKFC+LgwYOiZcuWolWrVlKbjRs3Cn9/f+l9fn6+aNiwoQgODhaxsbFix44dwt3dXUyZMkVqs2bNGmFlZSWWLFkibty4Ib0Kjzh9FlqtVrRr1040adJE7Ny5UyQlJYlDhw6Jjz/+WBqt+uGHH4oaNWoItVotNBqNaNeunejWrZt0jGnTpokKFSqITp06idjYWLF//35Rt25d8fbbb0ttbt++LaytrcWpU6d0vr9Lly6iWbNm4ujRo+LgwYOiTp06ok+fPtL+q1evCn9/f3H06FEhxMM/x61btxaNGjUS58+f1/k9yc/Plz6XmJgojh8/LkaMGCHq1q0rXa/HjTbWR3Hnk5GRIWrWrCkmTJgghBDi5MmTwtbWVmeU6dM4deqUcHd3F/369dM539TUVKnNxYsXhYODg5g4caI4e/asWLJkiVCpVGLHjh1Sm99++03Y2tqKlStXijNnzojhw4cLV1dXnVHZI0eOFNWqVRNRUVHin3/+EYGBgSIwMFCnnqSkJKFQKMSlS5d0tr/44otixowZJZ6HOf+cIbJkgYGBAoBYtmyZUb9Hn1GxDHalwGAnRHR0tGjcuLGwtbXVme4EgFixYsVjv6d69epFpkt59N8U4eHhIiAgQNjb2wtvb2/Rt29fcfXqVWn/ihUrinzm0qVLomvXrsLe3l5UrlxZfPDBBzrTdhQ3TQsAMXDgwMf/xughIyNDvPfee8LHx0dYW1sLX19f0bdvX5GcnCz27NkjrKysxIEDB6T2SUlJwtnZWSxdulQI8f9TiixdulT4+PgIOzs78cYbb4i7d+/qfM9bb70lJk+erLPtzp07ok+fPsLR0VE4OzuLQYMGiXv37ul8FwCxZ88eIYQQe/bsKfb3A4BISkp64u9b4TbP6tHzGTRokGjUqJHIzs6Wts2fP1+4ubnp/DnQ17Rp04o9l+rVq+u027Nnj2jatKmwsbERNWvWLPbP9KJFi0S1atWEjY2NaNWqlThy5IjO/qysLDF69GhRsWJF4eDgIF5//XVx48YNnTZff/21zjQqQjwM4NbW1uLKlSslnoc5/5whslS3bt0SCoVCABDJyclG/S59gp1CiBLudZEkIyMDLi4u0tQQhWVnZyMpKQl+fn4Gf8CdLN/06dOxadMmnefFinPy5El07twZFy5ceOqBMqbE0s6nNHJzc1GnTh388ssveP7556XtkyZNQlpaGr777rsSP8ufM0Sm5+eff0b//v3RuHFjnDhxwqjf9bgc8ig+Y0dkBho3bozZs2cjKSlJ7lIMwtLOpzSSk5Px8ccf64Q64OGI+MLL9xGReSiYNkmf57fLAkfFEpmJ0q7Dai4s7XyepGCgx6M++OADGaohomeRn5+PHTt2ADC9YMceOyIZTZ8+/Ym3YYmIyLRER0cjPT0dbm5uaNOmjdzl6GCPHREREZEe/P39sXjxYmRlZUnTc5kKBjsiIiIiPXh4eDx2InY58VYsERERkYVgsCMiIiIqpe3bt2P58uU6KwaZEgY7IiIiolJavHgxRo4cidWrV8tdSrEY7IiIiIhK4cGDB4iKigIAdO/eXeZqisdgR0RERFQKe/bsQXZ2NqpVq4YGDRrIXU6xGOzKsXfffRcKhQIKhQLW1tbw9PRE586d8dNPP0Gr1Zb6OCtXroSrq6vxCiUiIjIBhVebUCgUMldTPAa7cq5Lly64ceMGLl26hO3bt6Njx454//330b17d+Tn58tdHhERkUkQQpjsMmKFMdiZgAURCQiPTCx2X3hkIhZEJBjtu21tbeHl5YUqVaqgefPm+Pjjj/Hnn39i+/btWLlyJQAgLCwMjRo1QoUKFeDr64vRo0cjMzMTALB3714MGjQIarVa6v2bPn06AOB///sfWrZsCScnJ3h5eeGdd95Bamqq0c6FiIjIWOLi4pCcnAw7Ozt07NhR7nJKxGBnAlRKBcKKCXfhkYkIi0iASlm23b0vvfQSmjRpgo0bNwIAlEolwsPDERcXh1WrViEqKgofffQRAKBt27ZYuHAhnJ2dcePGDdy4cQMffvghACAvLw9ffPEFTpw4gU2bNuHSpUvlbn1QIiKyDMePH4dKpcJLL70EBwcHucspEVeeMAGhQXUAAGH/9cyFBtWRQt2EznWl/WWpXr16OHnyJABg3Lhx0vYaNWrgyy+/xMiRI7F06VLY2NjAxcUFCoUCXl5eOscYPHiw9OuaNWsiPDwczz33HDIzM+Ho6Fgm50FERGQI/fv3R/fu3XHnzh25S3ksBjsTUTjcLY46j1yNVrZQBzx8lqDgwdDdu3dj5syZOHfuHDIyMpCfn4/s7Gw8ePDgsf9qiYmJwfTp03HixAmkpaVJAzKSk5MREBBQJudBRERkKBUrVkTFihXlLuOxeCvWhIQG1YGNSolcjRY2KqVsoQ4Azp49Cz8/P1y6dAndu3dH48aNsWHDBsTExGDJkiUAgNzc3BI/f//+fYSEhMDZ2Rlr1qzBsWPH8Mcffzzxc0RERKZGn5ki5MZgZ0LCIxOlUJer0ZY4oMLYoqKicOrUKfTq1QsxMTHQarWYP38+2rRpg7p16xZZRsXGxgYajUZn27lz53Dnzh3MmjUL7du3R7169ThwgoiIzNKgQYPQtm1baXJiU8ZbsSbi0WfqCt4DMGrPXU5ODlJSUqDRaHDz5k3s2LEDM2fORPfu3TFgwACcPn0aeXl5WLRoEXr06IFDhw5h2bJlOseoUaMGMjMzERkZiSZNmsDBwQHVqlWDjY0NFi1ahJEjR+L06dP44osvjHYeRERExqDRaLB161bcuXMHVlamH5vYY2cCihsoERpUBxM61y12tKwh7dixA97e3qhRowa6dOmCPXv2IDw8HH/++SdUKhWaNGmCsLAwzJ49Gw0bNsSaNWswc+ZMnWO0bdsWI0eORO/eveHu7o45c+bA3d0dK1euxLp16xAQEIBZs2Zh3rx5RjsPIiIiYzh69Cju3LkDV1dXtG3bVu5ynkghhBByF2HqMjIy4OLiArVaDWdnZ5192dnZSEpKgp+fH+zs7J7q+Av+m9KkuJ658MhEaLQC4zvXfapjE5H5M8TPGSJ6Op988gm+/vpr9O7dG7/99pssNTwuhzzK9PsUy4HHhTY5B1AQERGVdwWrTXTv3l3mSkqHt2KJiIiIinH16lWcOHECCoUCXbp0kbucUmGwIyIiIirGtm3bAABt2rRB5cqVZa6mdHgrloiIiKgYtWrVQs+ePdG+fXu5Syk1BjsiIiKiYgQFBSEoKEjuMvTCW7EGwsHFRGQs/PlCRKXFYPeMrK2tAQAPHjyQuRIislQFP18Kft4QkfHt2bMHCQkJcpehN96KfUYqlQqurq7SclkODg5QKBQyV0VElkAIgQcPHiA1NRWurq5QqVRyl0RULgghMGTIECQlJWH79u1mMyIWYLAzCC8vLwDgWqhEZBSurq7SzxkiMr5z584hKSkJNjY2aNeundzl6IXBzgAUCgW8vb3h4eGBvLw8ucshIgtibW3NnjqiMlYwKXGHDh3g6OgoczX6kTXYzZw5Exs3bsS5c+dgb2+Ptm3bYvbs2fD395faZGdn44MPPsBvv/2GnJwchISEYOnSpfD09JTaJCcnY9SoUdizZw8cHR0xcOBAzJw5U2ex3r1792LChAmIi4uDr68vPv30U7z77rsGPR+VSsUfwERERGauINh169ZN5kr0J+vgiX379mHMmDE4cuQIIiIikJeXh+DgYNy/f19qM378eGzevBnr1q3Dvn37cP36dfTs2VPar9Fo0K1bN+Tm5uLw4cNYtWoVVq5cialTp0ptkpKS0K1bN3Ts2BGxsbEYN24chg4dip07d5bp+RIREZFpU6vVOHjwIADzDHYKYULj6G/dugUPDw/s27cPL7zwAtRqNdzd3fHLL7/gjTfeAPDwvnf9+vURHR2NNm3aYPv27ejevTuuX78u9eItW7YMkyZNwq1bt2BjY4NJkyZh69atOH36tPRdb7/9NtLT07Fjx44n1qXP4rtERERkvtatW4e33noL/v7+OHfunNzlANAvh5jUdCdqtRoA4ObmBgCIiYlBXl4eOnXqJLWpV68eqlWrhujoaABAdHQ0GjVqpHNrNiQkBBkZGYiLi5PaFD5GQZuCYxAREREBwPbt2wGYZ28dYEKDJ7RaLcaNG4fnn38eDRs2BACkpKTAxsYGrq6uOm09PT2RkpIitSkc6gr2F+x7XJuMjAxkZWXB3t5eZ19OTg5ycnKk9xkZGc9+gkRERGTyFi1ahNdffx21atWSu5SnYjLBbsyYMTh9+rR0X1tOM2fOxOeffy53GURERFTGKlSogB49eshdxlMziVuxY8eOxZYtW7Bnzx5UrVpV2u7l5YXc3Fykp6frtL9586Y0p5OXlxdu3rxZZH/Bvse1cXZ2LtJbBwBTpkyBWq2WXleuXHnmcyQiIiIyNlmDnRACY8eOxR9//IGoqCj4+fnp7G/RogWsra0RGRkpbYuPj0dycjICAwMBAIGBgTh16pTO5MARERFwdnZGQECA1KbwMQraFBzjUba2tnB2dtZ5ERERkWXr168fPv74Y+lRLnMk66jY0aNH45dffsGff/6pM3edi4uL1JM2atQobNu2DStXroSzszPee+89AMDhw4cBPJzupGnTpvDx8cGcOXOQkpKC/v37Y+jQofj6668BPJzupGHDhhgzZgwGDx6MqKgohIaGYuvWrQgJCXlinRwVS0REZNlu3LgBHx8fAMU/my8nvXKIkBGAYl8rVqyQ2mRlZYnRo0eLihUrCgcHB/H666+LGzdu6Bzn0qVLomvXrsLe3l5UrlxZfPDBByIvL0+nzZ49e0TTpk2FjY2NqFmzps53PIlarRYAhFqtfpbTJSIiIhP1ww8/CADiueeek7uUIvTJISY1j52pYo8dERGRZevZsyf++OMPTJ8+HdOmTZO7HB1mO48dERERUVnLyclBREQEAPOdv64Agx0RERGVawcOHEBmZiY8PT3RvHlzuct5Jgx2REREVK5t3boVwMPeOqXSvKOReVdPRERE9Izs7OxQqVIls78NC8g83Ym54OAJIiIiy6bRaCCEgJWVySzKJdEnh5he9URERERlTKVSyV2CQfBWLBEREVm0BREJCI9MLHbf9HVHsSAivowrMh4GOyIiIrJoKqUCYcWEu7nbTmNlzG2Ehc3H7du3ZarOsHgrloiIiCxaaFAdAEBYRIL0PjwyEUv2X0b6gZ/hnvoPKleuLGeJBsNgR0RERBavcLhbHHUeuRotat4/gz2Hf8O7778vc3WGw1uxREREVC6EBtWBjUqJXI0WNioFzq4PA2D+q00UxmBHRERE5UJ4ZOJ/oU6JXI1AVs0OcHR0xAsvvCB3aQbDYEdEREQWLzwyEWERCZjQuS4SvuqK5tbX4Nq+Hxr2nghbW1u5yzMYBjsiIiKyaIVDXcGzdtd2/Yj0Az/jRuUWJU6FYo44eIKIiIgsmkYrdEIdAMycORNbtmxBpTbe0GgtZxEuLilWClxSjIiIiOSiTw7hrVgiIiIiC8FgR0REROVGXl4eJk6ciIiICGg0GrnLMTgGOyIiIio3Dh48iHnz5qFv375QKBRyl2NwDHZERERUbmzduhUA0LVrVyiVlheDLO+MiIiIiEpQEOwsabWJwhjsiIiIqFy4ePEizp07B5VKheDgYLnLMQoGOyIiIioXtmzZAgBo3749XF1d5S3GSBjsiIiIqFz49ddfAQCvvvqqzJUYD4MdERERWbx79+7h1q1bUKlU6NOnj9zlGA2XFCMiIiKL5+TkhMTERJw9exaenp5yl2M07LEjIiKickGhUCAgIEDuMoyKwY6IiIgs2p07d5CTkyN3GWWCwY6IiIgs2qeffgpvb2+sWrVK7lKMjsGOiIiILFZubi7Wrl2LtLQ0+Pj4yF2O0THYERERkcXatm2bFOpeeuklucsxOgY7IiIislj/+9//AADvvPMOVCqVzNUYH4MdERERWaS0tDRptYn+/fvLXE3ZYLAjIiIii7Ru3Trk5uaiUaNGaNy4sdzllAkGOyIiIrJIP//8M4Dy01sHcOUJIiIislDLly/H//73P7zzzjtyl1JmFEIIIXcRpi4jIwMuLi5Qq9VwdnaWuxwiIiIqR/TJIbwVS0RERGQhGOyIiIjIosTExKB3797YsWOH3KWUOT5jR0RERBZl9erV+P3336FSqdClSxe5yylT7LEjIiIii5GXl4dff/0VQPkaDVuAwY6IiIgsxq5du3Dr1i14eHigc+fOcpdT5hjsiIiIyGIULCHWp08fWFmVvyfOGOyIiIjIImRkZODPP/8EAPTr10/mauTBYEdEREQWYcOGDcjOzka9evXQokULucuRRfnroyQiIiKL5OTkhMaNG+Ott96CQqGQuxxZcOWJUuDKE0REROYjPz/fop6v48oTREREVG5ZUqjTF4MdERERmTUhBDZu3IjMzEy5S5Edgx0RERGZtRMnTqBXr16oXr06cnJy5C5HVgx2REREZNYK5q576aWXYGtrK3M18mKwIyIiIrOl0Wjwyy+/ACifS4g9isGOiIiIzFZkZCRSUlJQqVIldOnSRe5yZMdgR0RERGar4DZs7969YWNjI3M18mOwIyIiIrOUmZmJjRs3Aii/S4g9isGOiIiIzNLBgwfx4MED1KpVC23atJG7HJNQfmfwIyIiIrPWpUsXXL58GZcvXy63S4g9isGOiIiIzFa1atVQrVo1ucswGbwVS0RERGYnPz9f7hJMEoMdERERmZ22bduiW7duSExMlLsUk8JbsURERGRW4uLicOzYMVhZWcHNzU3uckwKe+yIiIjIrPz8888AgJdffhmVKlWSuRrTwmBHREREZkOr1WLNmjUAuIRYcRjsiIiIyGzs27cPV65cgYuLC7p37y53OSaHwY6IiIjMRsFt2DfffBN2dnYyV2N6GOyIiIjILGRlZWH9+vUAeBu2JBwVS0RERGZj9uzZ2L17N9q1ayd3KSZJIYQQchdh6jIyMuDi4gK1Wg1nZ2e5yyEiIqJyRJ8cwluxRERERBaCt2KJiIjIZCyISIBKqUBoUB2d7Rs3bsTG+CwENGiAj19pKk9xZoA9dkRERGQyVEoFwiISEB6pu1TYZ78dxkG1K06dPClTZeaBPXZERERkMgp66sIiEqT3U9dG437NjlAfXINvNi6UsTrTJ2uP3f79+9GjRw/4+PhAoVBg06ZNOvvfffddKBQKnVeXLl102ty9exd9+/aFs7MzXF1dMWTIEGRmZuq0OXnyJNq3bw87Ozv4+vpizpw5xj41IiIiekqhQXUwoXNdhEUkoO4n27H6+F2kH/gZbV3U8PDwkLs8kyZrsLt//z6aNGmCJUuWlNimS5cuuHHjhvT69ddfdfb37dsXcXFxiIiIwJYtW7B//34MHz5c2p+RkYHg4GBUr14dMTExmDt3LqZPn47vvvvOaOdFREREzyY0qA5sVErkarSAJg/qw7+hX79+cpdl8p75VqxGo8GpU6dQvXp1VKxYUa/Pdu3aFV27dn1sG1tbW3h5eRW77+zZs9ixYweOHTuGli1bAgAWLVqEl19+GfPmzYOPjw/WrFmD3Nxc/PTTT7CxsUGDBg0QGxuLsLAwnQBIREREpiM8MhG5Gi2slEA+rOHeYQBeffVVucsyeXr32I0bNw4//vgjgIeh7sUXX0Tz5s3h6+uLvXv3Gro+7N27Fx4eHvD398eoUaNw584daV90dDRcXV2lUAcAnTp1glKpxNGjR6U2L7zwAmxsbKQ2ISEhiI+PR1pamsHrJSIiomcTHpmIsIgETOhcF53ubkb6gZ/h0Pot/BB9Te7STJ7ewW79+vVo0qQJAGDz5s1ISkrCuXPnMH78eHzyyScGLa5Lly5YvXo1IiMjMXv2bOzbtw9du3aFRqMBAKSkpBS5125lZQU3NzekpKRIbTw9PXXaFLwvaPOonJwcZGRk6LyIiIjI+AqHutCgOnjw4AHu/70er9WyKna0LOnS+1bs7du3pVuj27Ztw5tvvom6deti8ODB+Oabbwxa3Ntvvy39ulGjRmjcuDFq1aqFvXv3IigoyKDfVdjMmTPx+eefG+34REREVDyNVkihDgBWr16NefPmoVKlSqi59yI0Wi6Y9Th699h5enrizJkz0Gg02LFjBzp37gwAePDgAVQqlcELLKxmzZqoXLkyzp8/DwDw8vJCamqqTpv8/HzcvXtXCp9eXl64efOmTpuC9yU9uzdlyhSo1WrpdeXKFUOfChERERVjfKFQV8DDwwMqlQqhQXUwvnNdmSozD3oHu0GDBuGtt95Cw4YNoVAo0KlTJwDA0aNHUa9ePYMXWNjVq1dx584deHt7AwACAwORnp6OmJgYqU1UVBS0Wi1at24ttdm/fz/y8vKkNhEREfD39y9xsIetrS2cnZ11XkRERFR21Go1Ll26JHcZZkfvYDd9+nT88MMPGD58OA4dOgRbW1sAgEqlwuTJk/U6VmZmJmJjYxEbGwsASEpKQmxsLJKTk5GZmYmJEyfiyJEjuHTpEiIjI/Hqq6+idu3aCAkJAQDUr18fXbp0wbBhw/D333/j0KFDGDt2LN5++234+PgAAN555x3Y2NhgyJAhiIuLw9q1a/HNN99gwoQJ+p46ERERlZElS5bAz88Po0aNkrsU8yKeQVZW1rN8XOzZs0cAKPIaOHCgePDggQgODhbu7u7C2tpaVK9eXQwbNkykpKToHOPOnTuiT58+wtHRUTg7O4tBgwaJe/fu6bQ5ceKEaNeunbC1tRVVqlQRs2bN0qtOtVotAAi1Wv1M50tERERPlp2dLby8vAQAsXr1arnLkZ0+OUQhhNDrKUSNRoOvv/4ay5Ytw82bN5GQkICaNWvis88+Q40aNTBkyBBDZ0/ZZWRkwMXFBWq1mrdliYiIjGz16tUYOHAgfHx8kJSUpDNlWXmkTw7R+1bsV199hZUrV2LOnDk6v9ENGzbEDz/8oH+1RERERP8RQiAsLAwA8N5775X7UKcvvYPd6tWr8d1336Fv3746o2CbNGmCc+fOGbQ4IiIiKl+ioqJw4sQJODg4cIWop6B3sLt27Rpq165dZLtWq9UZeUpERESkr4LeukGDBsHNzU3masyP3sEuICAABw4cKLJ9/fr1aNasmUGKIiIiovInPT0dx44dg0KhwLhx4+QuxyzpvfLE1KlTMXDgQFy7dg1arRYbN25EfHw8Vq9ejS1bthijRiIiIioHXF1dcenSJezbt6/Yu4P0ZHqPigWAAwcOYMaMGThx4gQyMzPRvHlzTJ06FcHBwcaoUXYcFUtERERy0SeHPFWwK28Y7IiIiIzrwoUL8PPzg1Kp91NiFs+o050QERERGVJ2djbatm2LBg0a4OLFi3KXY9ZK9YxdxYoVoVAoSnXAu3fvPlNBREREVL6sWbMGqampsLW1ha+vr9zlmLVSBbuFCxcauQwiIiIqjwpPSBwaGgpra2uZKzJvfMauFPiMHRERkXHs3LkTXbp0gaOjI65cuQJXV1e5SzI5Rn3Gbtu2bdi5c2eR7bt27cL27dv1PRwRERGVYwW9dUOHDmWoMwC9g93kyZOh0WiKbNdqtZg8ebJBiiIiIiLLd/r0aezatQtKpRKhoaFyl2MR9A52iYmJCAgIKLK9Xr16OH/+vEGKIiIiIsu3du1aAEDPnj3h5+cnczWWQe+VJ1xcXHDx4kXUqFFDZ/v58+dRoUIFQ9VFREREFm7GjBno0KEDPDw85C7FYujdY/fqq69i3LhxuHDhgrTt/Pnz+OCDD/DKK68YtDgiIiKyXAqFAkFBQWjUqJHcpVgMvYPdnDlzUKFCBdSrVw9+fn7w8/ND/fr1UalSJcybN88YNRIREZEFyc7ORmZmptxlWCS9g52LiwsOHz6MrVu3YvTo0fjggw8QGRmJqKgojmYhIiKiJ1q9ejWqVq0qjYglw9H7GTvgYddpcHAwgoODDV0PERERWTCtVouwsDCo1WqoVCq5y7E4pQp24eHhGD58OOzs7BAeHv7YthyuTERERCXZvn074uPj4ezsjMGDB8tdjsUp1coTfn5++Oeff1CpUqXHDkdWKBQWuXgvV54gIiIyjKCgIERFReHDDz/E3Llz5S7HLOiTQ7ikWCkw2BERET272NhYNGvWDCqVChcvXkS1atXkLsksGHVJsRkzZuDBgwdFtmdlZWHGjBn6Ho6IiIjKiYLBEm+99RZDnZHo3WOnUqlw48aNIpMJ3rlzBx4eHsUuN2bu2GNHRET0bNLT0+Hl5YWcnBwcO3YMLVu2lLsks6FPDtF7VKwQAgqFosj2EydOwM3NTd/DERERUTng6uqKuLg4/PXXXwx1RlTqYFexYkUoFAooFArUrVtXJ9xpNBpkZmZi5MiRRimSiIiIzF+tWrUwfvx4ucuwaKUOdgsXLoQQAoMHD8bnn38OFxcXaZ+NjQ1q1KiBwMBAoxRJRERE5is7Oxt2dnZyl1EulDrYDRw4EPn5+VAoFHjppZfg6+trzLqIiIjIAmg0GjRt2hQNGjRAeHg4qlSpIndJFk2vZ+ysrKwwatQonD171lj1EBERkQXZsmUL4uPjkZqayqVHy4De0520atUKx48fN0YtREREZGHmz58PABg5ciQqVKggczWWT+9RsaNHj8YHH3yAq1evokWLFkUuUuPGjQ1WHBEREZmvY8eO4cCBA7C2tsbYsWPlLqdc0DvYvf322wB014RVKBTSNCiWOI8dERER6W/BggUAHmYHHx8fmaspH/QOdklJScaog4iIiCzIlStX8PvvvwMApzgpQ3oHu+rVqxujDiIiIrIg3377LTQaDTp27IhmzZrJXU65oXewK3DmzBkkJycjNzdXZ/srr7zyzEURERGRefv444/h7e2Nhg0byl1KuaL3WrEXL17E66+/jlOnTknP1gGQVqKwxGfsuFYsERERyUWfHKL3dCfvv/8+/Pz8kJqaCgcHB8TFxWH//v1o2bIl9u7d+7Q1ExERkQXQarXQs8+IDEjvYBcdHY0ZM2agcuXKUCqVUCqVaNeuHWbOnKkzUpaIiIjKn40bN6Jp06ZYt26d3KWUS3oHO41GAycnJwBA5cqVcf36dQAPB1XEx8cbtjoiIiIyK2FhYTh58iROnToldynlkt6DJxo2bIgTJ07Az88PrVu3xpw5c2BjY4PvvvsONWvWNEaNREREZAaio6MRHR0NGxsbjB49Wu5yyiW9g92nn36K+/fvAwBmzJiB7t27o3379qhUqRLWrl1r8AKJiIjIPISFhQEA+vbtCy8vL5mrKZ/0HhVbnLt376JixYrSyFhLw1GxREREj5eUlITatWtDq9Xi1KlTnObEgPTJIU89j11hbm5uhjgMERERmanw8HBotVoEBwcz1MlI78ETRERERIXdunULy5YtAwBMmDBB5mrKN4P02BEREVH55e7ujr/++gvr169HcHCw3OWUawZ5xs7S8Rk7IiIikovBV55o3rw50tLSADwcCfvgwYNnr5KIiIjMXkE+INNQqmB39uxZaYqTzz//HJmZmUYtioiIiExfREQEqlWrhjlz5shdCv2nVM/YNW3aFIMGDUK7du0ghMC8efPg6OhYbNupU6catEAiIiIyPRqNBhMnTkRmZqa0ChXJr1TP2MXHx2PatGm4cOEC/v33XwQEBMDKqmgmVCgU+Pfff41SqJz4jB0REZGuVatW4d1334WLiwsuXLiASpUqyV2SxdInh+g9eEKpVCIlJQUeHh7PVKQ5YbAjIiL6f1lZWahbty6uXr2K2bNn46OPPpK7JItm1AmKtVrtUxdGRERE5m/hwoW4evUqqlWrhtDQULnLoUKeah67CxcuYOHChTh79iwAICAgAO+//z5q1apl0OKIiIjItNy6dQszZ84EAHz11Vews7OTuSIqTO+VJ3bu3ImAgAD8/fffaNy4MRo3boyjR4+iQYMGiIiIMEaNREREZCL279+PrKwsNG/eHO+8847c5dAj9H7GrlmzZggJCcGsWbN0tk+ePBm7du3i4AkiIiILd/78eWRmZqJp06Zyl1IuGHXwhJ2dHU6dOoU6derobE9ISEDjxo2RnZ2tf8UmjsGOiIiI5GLwlScKc3d3R2xsbJHtsbGx5WqkLBERUXly4sQJHD9+XO4y6An0HjwxbNgwDB8+HBcvXkTbtm0BAIcOHcLs2bMxYcIEgxdIRERE8hJCYNSoUThy5Ah++OEHDB48GACwICIBKqUCoUF1inwmPDIRGq3A+M51y7rcck3vYPfZZ5/ByckJ8+fPx5QpUwAAPj4+mD59Ooc8ExERWaANGzYgOjoaDg4O6Nq1q7RdpVQgLCIBAHTCXXhkIsIiEjCBoa7M6R3sFAoFxo8fj/Hjx+PevXsAACcnJ4MXRkRERPLLzc3F5MmTAQAffvghvL29pX0FYa5wuCsc6orrySPjeqp57Aow0BEREVm25cuX48KFC/D09MTEiROL7C8c7hZHnUeuRstQJyO9B08QERFR+aBWq/H5558DAD7//HM4OjoW2y40qA5sVErkarSwUSkZ6mTEYEdERETFmjVrFu7cuYN69ephyJAhJbYLj0yUQl2uRovwyMQyrJIKe6ZbsURERGS5atasicqVK2P27Nmwsio+Mjz6TF3BewDsuZOBXj12eXl5CAoKQmIikzgREZGlGzZsGJKSktCjR49i9xc3UCI0qA4mdK6LsIgE9tzJQK8eO2tra5w8edJYtRAREZGJKem5OgDQaEWxAyUK3mu0ei1uRQag95Ji48ePh62tbZG1Yi0ZlxQjIqLyQgiBESNGIDg4GL169YJCoZC7pHJPnxyi9zN2+fn5+Omnn7B79260aNECFSpU0NkfFham7yGJiIjIROzcuRPff/89Vq1ahdatW8PX11fukkgPege706dPo3nz5gCAhIQEnX1M9UREROZLo9FIc9WNHTuWoc4M6R3s9uzZY4w6iIiISGarVq3C6dOn4erqik8++UTucugpPPU8dufPn8fOnTuRlZUF4OE9eSIiIjJP9+/fx2effQYA+PTTT+Hm5iZzRfQ09A52d+7cQVBQEOrWrYuXX34ZN27cAAAMGTIEH3zwgcELJCIiIuNbsGABrl+/jho1amDs2LFyl0NPSe9gN378eFhbWyM5ORkODg7S9t69e2PHjh16HWv//v3o0aMHfHx8oFAosGnTJp39QghMnToV3t7esLe3R6dOnYrMoXf37l307dsXzs7OcHV1xZAhQ5CZmanT5uTJk2jfvj3s7Ozg6+uLOXPm6HfSREREFiwtLQ2zZ88GAHz99dewtbWVuSJ6WnoHu127dmH27NmoWrWqzvY6derg8uXLeh3r/v37aNKkCZYsWVLs/jlz5iA8PBzLli3D0aNHUaFCBYSEhCA7O1tq07dvX8TFxSEiIgJbtmzB/v37MXz4cGl/RkYGgoODUb16dcTExGDu3LmYPn06vvvuO71qJSIislQVK1bEhg0bMHDgQPTu3VvucuhZCD05OjqKhIQE6dcXLlwQQghx7Ngx4ebmpu/hJADEH3/8Ib3XarXCy8tLzJ07V9qWnp4ubG1txa+//iqEEOLMmTMCgDh27JjUZvv27UKhUIhr164JIYRYunSpqFixosjJyZHaTJo0Sfj7+5e6NrVaLQAItVr9tKdHRERE9FT0ySF699i1b98eq1evlt4rFApotVrMmTMHHTt2NFjgTEpKQkpKCjp16iRtc3FxQevWrREdHQ0AiI6OhqurK1q2bCm16dSpE5RKJY4ePSq1eeGFF2BjYyO1CQkJQXx8PNLS0or97pycHGRkZOi8iIiILBH/jrMsege7OXPm4LvvvkPXrl2Rm5uLjz76CA0bNsT+/ful+/OGkJKSAgDw9PTU2e7p6SntS0lJgYeHh85+KysruLm56bQp7hiFv+NRM2fOhIuLi/TiPD5ERGSJDhw4wGfPLYzewa5hw4ZISEhAu3bt8Oqrr+L+/fvo2bMnjh8/jlq1ahmjxjI3ZcoUqNVq6XXlyhW5SyIiIjIoIQQmTpyIjIwMXLx4Ue5yyED0nqAYeHhL1NgTF3p5eQEAbt68CW9vb2n7zZs30bRpU6lNamqqzufy8/Nx9+5d6fNeXl64efOmTpuC9wVtHmVra8sRQUREZNHWrVsnDUycPn263OWQgTzVBMVpaWmYN28ehgwZgiFDhmD+/Pm4e/euQQvz8/ODl5cXIiMjpW0ZGRk4evQoAgMDAQCBgYFIT09HTEyM1CYqKgparRatW7eW2uzfvx95eXlSm4iICPj7+6NixYoGrZmIiMgc5OTkYMqUKQCAjz76qMSODjI/ege7/fv3o0aNGggPD0daWhrS0tIQHh4OPz8/7N+/X69jZWZmIjY2FrGxsQAeDpiIjY1FcnIyFAoFxo0bhy+//BJ//fUXTp06hQEDBsDHxwevvfYaAKB+/fro0qULhg0bhr///huHDh3C2LFj8fbbb8PHxwcA8M4778DGxgZDhgxBXFwc1q5di2+++QYTJkzQ99SJiIgswsKFC3Hx4kV4e3tzcQFLo++Q24YNG4phw4aJ/Px8aVt+fr4YPny4aNiwoV7H2rNnjwBQ5DVw4EAhxMMpTz777DPh6ekpbG1tRVBQkIiPj9c5xp07d0SfPn2Eo6OjcHZ2FoMGDRL37t3TaXPixAnRrl07YWtrK6pUqSJmzZqlV52c7oSIiCzF6dOnhY2NjQAgVq5cKXc5VAr65BCFEPot8mpvb4/Y2Fj4+/vrbI+Pj0fTpk2ltWMtSUZGBlxcXKBWq+Hs7Cx3OURERE9t9erVGDp0KIKDg7F582YoFAq5S6In0CeH6D14onnz5jh79myRYHf27Fk0adJE38MRERFRGRowYACaNGkCDw8PhjoLVKpgd/LkSenXoaGheP/993H+/Hm0adMGAHDkyBEsWbIEs2bNMk6VREREZDDsiLFcpboVq1QqoVAo8KSmCoUCGo3GYMWZCt6KJSIic5aTk4P+/ftj0qRJaNGihdzlkJ4Mfis2KSnJIIURERFR2Zs2bRrWrVuHw4cP48KFC5yr1YKVKthVr17d2HUQERGRERw+fBhz584FACxevJihzsI91coT169fx8GDB5GamgqtVquzLzQ01CCFERER0bO5f/8+Bg4cCK1WiwEDBkjzwJLl0jvYrVy5EiNGjICNjQ0qVaqkM6JGoVAw2BEREZmIKVOm4Pz586hSpQq++eYbucuhMqD3PHa+vr4YOXIkpkyZAqXyqVYkMzscPEFEROZmz549eOmllwAAO3fuRHBwsMwV0dPSJ4foncwePHiAt99+u9yEOiIiInO0cuVKAMCIESMY6soRvdPZkCFDsG7dOmPUQkRERAayYsUKfPfdd9LACSof9L4Vq9Fo0L17d2RlZaFRo0awtrbW2R8WFmbQAk0Bb8USERGRXIy6pNjMmTOxc+dOaUmxRwdPEBERkTzu3r2LOXPm4NNPP4Wjo2OR/QsiEqBSKhAaVKfIvvDIRGi0AuM71y2LUslI9A528+fPx08//YR3333XCOUQERHR0woNDcWaNWtw6tQpbN26tch+lVKBsIiEh20LhbvwyESERSRgAkOd2dM72Nna2uL55583Ri1ERET0lDZs2IA1a9ZAqVRi2rRpxbYpCHOFw13hUFdcTx6ZF72D3fvvv49FixYhPDzcGPUQERGRnlJTUzFy5EgAD+eua9WqVYltC4e7xVHnkavRMtRZEL2D3d9//42oqChs2bIFDRo0KDJ4YuPGjQYrjoiIiB5PCIERI0bg9u3baNy4MaZOnfrEz4QG1ZFCnY1KyVBnQfQOdq6urujZs6cxaiEiIiI9rVmzBps2bYK1tTVWr14NGxubJ34mPDJRCnW5Gi3CIxMZ7iyE3sFuxYoVxqiDiIiI9JSfn4/p06cDAKZNm4YmTZo88TOPPlNX8B4Aw50F0DvYERERkWmwsrLCgQMHsGDBAkyaNOmJ7YsbKFHcgAoyX3oHOz8/v8fOV3fx4sVnKoiIiIhKz9vbG3PmzClVW41WFDtQouC9RqvXmgVkgvQOduPGjdN5n5eXh+PHj2PHjh2YOHGioeoiIiKiEiQlJeHEiRN47bXX9Prc4yYfZk+dZXiq6U6Ks2TJEvzzzz/PXBARERGVTKvVYtCgQdi3bx9mz56Njz76SO6SyIQoDXWgrl27YsOGDYY6HBERERVj0aJF2LdvHypUqIBevXrJXQ6ZGIMFu/Xr18PNzc1QhyMiIqJHxMfHY/LkyQCAefPmoVatWjJXRKZG71uxzZo10xk8IYRASkoKbt26haVLlxq0OCIiInooPz8fAwcORHZ2Njp37owRI0bIXRKZIL2D3aMPaiqVSri7u6NDhw6oV6+eoeoiIiKiQubOnYujR4/CxcUFP/7442NnqKDySyGE4NjmJ8jIyICLiwvUajWcnZ3lLoeIiMqZa9euoWbNmsjNzcWqVaswYMAAuUuiMqRPDuEExURERCZoQUQCVEoFQoPqoEqVKtiwYQM2b96M/v37IzwyERqteOz0JVQ+lXrwhFKphEqleuzLyoo5kYiIyBBUSgXCIhIQHpkIAOjevTuWL1+ORVHnEfZf6CN6VKmT2B9//FHivujoaISHh0Or1RqkKCIiovIuNKgOYmJidJb6Km5JMKLCnukZu4Jh15s3b0bfvn0xY8YMVK9e3ZD1mQQ+Y0dERGVt48aN6NWrF7yCBsG2ZS/YqJTI1WgZ6sohfXLIU81jd/36dQwbNgyNGjVCfn4+YmNjsWrVKosMdURERGXt33//Rf/+/QEAbwY4SqHORqVkqKPH0ivYqdVqTJo0CbVr10ZcXBwiIyOxefNmNGzY0Fj1ERERlSvXr1/HK6+8ggcPHiAkJAQ1e4yRQl2uRis9c0dUnFI/YzdnzhzMnj0bXl5e+PXXX/Hqq68asy4iIqJy58GDB3j11Vdx7do11K9fHy+9NxcLI89Lt18LnrEDwJ47Klapn7FTKpWwt7dHp06doFKpSmy3ceNGgxVnKviMHRERGZtWq0WfPn3w+++/o1KlShj77VasjLld5Jk6DqAof4wyj92AAQM4yzUREZGR3Lt3D5cuXYK1tTU2btyImJyKmNDZrUh4K3iv0XJ9ASqKK0+UAnvsiIioLGRlZeHIkSPo2LGj3KWQCTH6qFgiIiIyjLS0NOnX9vb2DHX0TBjsiIiIZHLlyhUEBATgs88+4yT/ZBAMdkRERDLIzMzEK6+8gpSUFPz555/IysqSuySyAAx2REREZUyr1aJ///6IjY2Fh4cHNm/ejAoVKshdFlkABjsiIqIy9sknn2DTpk2wsbHBH3/8wZWbyGAY7IiIiMrQ6tWrMWvWLADAjz/+iLZt28pcEVkSBjsiIqIycuXKFQwbNgwA8PHHH6Nfv34yV0SWptQTFBMREdGz8fX1xY8//ojt27fjiy++kLscskCcoLgUOEExERERyYUTFBMREZkIjUaDKVOmICUlRe5SqBxgsCMiIjKijz76CLNmzUKHDh2Qn58vdzlk4RjsiIiIjOSHH35AWFgYAGDGjBmwsuKj7WRcDHZERERGsHfvXowaNQoA8Pnnn+Ott96SuSIqDxjsiIiIDOz8+fPo1asX8vPz8fbbb+Ozzz6TuyQqJxjsiIiIDCg9PR09evTA3bt30apVK/z0009QKBRyl0XlBIMdERGRAanVagBA1apVsWnTJtjb28tcEZUnfIqTiIjoGS2ISIBKqUBoUB1Ur14d0dHRSElJgbe3N8IjE6HRCozvXFfuMqkcYI8dERHRM8rNyUZYRALCIxMBAK6urqhXrx7CIxMR9l/oIyoL7LEjIiJ6Brdv38baqe8i3d4fYXi49mtoUB0p1E3oXBehQXVkrpLKCwY7IiKip3Tz5k106tQJp0+fhofHNfQdMwZhEQlYHHUeuRotQx2VOd6KJSIiegrXrl3Diy++iNOnT8PHxwf79u3DV++0g41KiVyNFjYqJUMdlTkGOyIiIj1dvnwZL7zwAuLj41GtWjXs379feqauINTlarTSM3dEZYW3YomIiPRw9+5dvPDCC0hOTkbNmjURFRWF6tWrF3mmruA9APbcUZlhsCMiItJDxYoV0adPH/zxxx+IiopClSpVih0oUfBfhjsqSwohhJC7CFOXkZEBFxcXqNVqODs7y10OERHJTAgh/d0A6M5j9yjOY0fPSp8cwmBXCgx2RETl27///otZs2Zh1apVXEmCypw+OYS3YomIiB7j6NGjCAkJgVqtRvXq1TF37ly5SyIqEUfFEhERleDAgQPo1KkT1Go1nn/+eXz22Wdyl0T0WAx2RERExYiMjESXLl2QmZmJjh07YseOHXwch0wegx0REdEjtm/fjm7duuHBgwcICQnB1q1b4ejoKHdZRE/EYEdERFTIgwcPMGjQIOTk5OCVV17Bn3/+yQETZDYY7IiIiApxcHDAX3/9hUGDBmH9+vWwtbWVuySiUuOoWCIiIjxcUcLNzQ0A0KpVK7Rq1Urmioj0xx47IiIqdxZEJOis4/r999+jVq1a+OeffxAemYgF/60WQWRuTDrYTZ8+HQqFQudVr149aX92djbGjBmDSpUqwdHREb169cLNmzd1jpGcnIxu3brBwcEBHh4emDhxIvLz88v6VIiIyISolAqE/RfuFi9ejOHDhyM9PR2fr/8bYf+tIkFkjkz+VmyDBg2we/du6b2V1f+XPH78eGzduhXr1q2Di4sLxo4di549e+LQoUMAAI1Gg27dusHLywuHDx/GjRs3MGDAAFhbW+Prr78u83MhIiLTEBpUB1qtFmERCUg/8PDvjK4fhuMUquus90pkbkx6SbHp06dj06ZNiI2NLbJPrVbD3d0dv/zyC9544w0AwLlz51C/fn1ER0ejTZs22L59O7p3747r16/D09MTALBs2TJMmjQJt27dgo2NTanq4JJiRESW5c6dO3jnnXdwNNMNru37QQkttFAy1JFJ0ieHmPStWABITEyEj48Patasib59+yI5ORkAEBMTg7y8PHTq1ElqW69ePVSrVg3R0dEAgOjoaDRq1EgKdQAQEhKCjIwMxMXFle2JEBGRSbh69SpatGiBXbt2Iff4n1ApBLRQwkalZKgjs2fSwa5169ZYuXIlduzYgW+//RZJSUlo37497t27h5SUFNjY2MDV1VXnM56enkhJSQEApKSk6IS6gv0F+0qSk5ODjIwMnRcREVkGHx8fNGzYELVq1cK473dCIxSwUSmRq9HqDKggMkcm/Yxd165dpV83btwYrVu3RvXq1fH7778bdbLImTNn4vPPPzfa8YmIqGzl5ORAq9XC3t4eSqUSP//8M747dAVLDyRLt1/DIxMR9t9oWPbckbky6R67R7m6uqJu3bo4f/48vLy8kJubi/T0dJ02N2/ehJeXFwDAy8uryCjZgvcFbYozZcoUqNVq6XXlyhXDnggREZWZq1ev4sUXX8TIkSNR8Fj56phbOqEOeBjmJnSuK42WJTJHZhXsMjMzceHCBXh7e6NFixawtrZGZGSktD8+Ph7JyckIDAwEAAQGBuLUqVNITU2V2kRERMDZ2RkBAQElfo+trS2cnZ11XkREZH727NmD5s2b4+jRo/jrr7+kf6hrtKLYgRIF4U6jNdlxhUSPZdKjYj/88EP06NED1atXx/Xr1zFt2jTExsbizJkzcHd3x6hRo7Bt2zasXLkSzs7OeO+99wAAhw8fBvBwupOmTZvCx8cHc+bMQUpKCvr374+hQ4fqNd0JR8USEclnwX/zyhV3ezQ8MhEarcD4znV1tgshMH/+fEyaNAlarRZNmzbFhg0bULNmzbIqm8hg9MkhJv2M3dWrV9GnTx/cuXMH7u7uaNeuHY4cOQJ3d3cAwIIFC6BUKtGrVy/k5OQgJCQES5culT6vUqmwZcsWjBo1CoGBgahQoQIGDhyIGTNmyHVKRESkp4LJhAHdZ98Knomb8Eiou3fvHgYPHoz169cDAAYMGIBvv/0WDg4OZVc0kUxMusfOVLDHjohIXoVDXOGBDo/eThVC4MUXX8SBAwdgbW2Nb775BiNHjoRCwZUkyHzpk0MY7EqBwY6ISH4FYa5gapKSJhPevXs3Bg8ejLVr10rPXBOZM4uaoJiIiAh4eBu2INQVnkw4Pz8fp0+fltp16tQJCQkJDHVULjHYERGRWQiPTJRCXcFkwrdu3UKXLl3w/PPPIzHx/6cosbOzk7FSIvkw2BERkckr/ExdwlddpfnmWgz8FJGRkcjPz0dCQoLcZRLJzqRHxRIRERU3UML+4l7cO7wPTm37oKaVFf6aORoNGjSQuVIi+THYERGRSSs8mXB2djbGjh2LH3/8EQBQp04dvDZ8BEMd0X8Y7IiIyKQVnnw4PDwcP/74IxQKBb788ktMnjwZSiWfKiIqwGBHRERmY9y4cThw4ADee+89BAcHy10OkcnhP3OIiMhk7du3D7169UJeXh4AwMbGBps3b2aoIyoBgx0REZmctLQ0DB06FB06dMDGjRsRHh4ud0lEZoG3YomIyGQIIbBu3TqEhobi5s2bAIARI0ZgyJAhMldGZB4Y7IiIyCRcuXIFo0ePxpYtWwAA9erVw/fff4927drJXBmR+eCtWCIiMpgFEQkIj0wsdl94ZCIWRJQ8iXBBqLO2tsa0adMQGxvLUEekJ/bYERGRwaiUCoT9F94KJhMGdCcZLsn8+fORlZWF8PBwBAQEGL1WIkvEYEdERAZTEOYKh7viVo7Izs7GF198gaysLISFhQEA6tati927d8tTOJGFUAghhNxFmLqMjAy4uLhArVbD2dlZ7nKIiExeQZizUSmRq9HqhLq9e/di+PDhSExMhEKhwJkzZ1CvXj2ZKyYyXfrkED5jR0REBhcaVEcKdTYqJUKD6khTmHTs2BGJiYnw9vbGhg0bGOqIDIi3YomIyODCIxOlUJer0WLYNxuxeeZoaQqTkSNHYtasWXBxcZG5UiLLwh47IiIyqMLP1CV81RWj2vki4oYtsmt1RP369XHgwAF8++23DHVERsBgR0REBlMQ6t7r4IexHWsBACZ1b4yOle/DtX0/DF2wgVOYEBkRgx0RERlMbl4enndOw7xBQVi/fr20fcWHb2FC57pQKFUyVkdk+fiMHRFRObYgIgEqpUJnzrkC4ZGJ0GgFxj9m7rkCeXl5WLVqFcJnzMCVK1cAACtWrMBbb70ltSnuO4jIsNhjR0RUjhVMKPzoahEFt1RVSsVjP6/RaPDzzz+jfv36GDZsGK5cuQIfHx98++23+PPPP41ZOhEVgz12RETlWGknFC7JgAED8MsvvwAA3N3d8fHHH2PEiBGwt7c3buFEVCwGOyKicq5wuFscdb7IhMKFCSGQn58Pa2trAMDAgQOxfft2fPTRRxg7diwcHR3LtHYi0sVbsUREVOyEwoUJIbB7924EBgZi5syZ0vbOnTvj8uXLmDx5MkMdkQlgsCMioiITChd+5u7gwYPo2LEjOnfujKNHj2LZsmXIy8sDACgUCjg5OclVNhE9grdiiYjKuUefqSt4f+3aNZz4ZRZ27twJALCxscHIkSMxZcoU6VYsEZkWBjsionKsuIESoUF1EB0djbVngPR7FWFlZYXBgwfj008/ha+vr8wVE9Hj8FYsEZGZWFDMtCQFwiMTseC/ka360GiFFOq0Wq20/Yu32+Le4V/RpFkzxMfHY/ny5Qx1RGaAwY6IyEw865xzxRnXqQ4CxGX07NkT7777rrS9du3aOLs+DPsWf4SaNWs+a+lEVEZ4K5aIyEw865xzhanVaqxatQpLly5FfHw8AECpVOKLL75A9erVAQAeHh4GPgMiMjYGOyIiM6LPnHPFOXv2LL755hv8/PPPuH//PgDAyckJAwYMwOjRo6VQR0TmicGOiMjMhAbVkUJdcXPOPU5UVBSWL18OAAgICMCYMWPQv39/TllCZCEY7IiIzExxc84VF+6uXbuG7777Dg0aNMBbb70FAOjfvz8OHz6MYcOG4cUXX4RCof9zeURkuhjsiIjMSElzzgEPe/KEENi3bx+WLFmCP/74AxqNBk2bNsWbb74JhUIBZ2dnrFmzRuazICJjYbAjIjITJc05Bzx85u7vv/9G7JqZiIuLkz7Tvn17jBkzBkII9s4RlQMMdkRERrDgv+lHirtFGh6ZCI1WYHznunods/Ccc4WFBtXBr7/+ii0H46GOi0OFChXQr18/jBkzBo0aNXqm8yAi88J57IiIjMAYc86N/y/UZWVlYd26dbh27Zq075PXmsPrdgy++eYbXLt2DcuWLWOoIyqH2GNHRGQEhpxzDgBycnKwY8cO/P777/jrr7+QmZmJzz77DDNmzAAAvPzyy+jWrRtvtxKVcwx2RERG8qxzzuXn52PXrl1Yu3YtNm3ahIyMDGlftWrV4OLiIr1XKnkDhogAhRBCyF2EqcvIyICLiwvUajWcnZ3lLoeIzEzdT7ZL05MkfNW11J/LycmBp6cn1Go1AKBKlSp488030bt3b7Ru3Zq9c0TlhD45hD12RERGVJo55zQaDfbt24fff/8dJ06cwOHDh6FQKGBra4vhw4cjKysLvXv3Rtu2bdkzR0SPxWBHRGQkj5tzbmzHWjh06BDWrl2L9evX4+bNm9LnYmJi0LJlSwDAnDlzZKmdiMwTgx0REQw/PcmT5pybO3curu78XmpfsWJF9OzZE71790bTpk2f7WSIqNxinz4REQw/PUnBnHOjX/TDoUOHcOHCBQAPw91rtaxwL/M+XFxcMHDgQGzbtg03b97EDz/8gM6dO8PKiv/mJqKnw58eREQw7PQkFy9ehP2FKOzatQuf9YxERkYGxo0bhwULFgAA5g0KQpBXLoKDf4Ktra3hT4aIyi0GOyKi/zzL9CRZWVn48MMPsWvXLpw/f15nn5ubm857Kysr9OjRw3CFExH9h8GOiKiQ0KA6UqizUSmLDXUajQb//vsvkpOT0atXLwCAnZ0dNm/ejCtXrsDKygpt27ZFcHAwgoOD0bx5c6hUqrI+FSIqhxjsiIgKKWl6kitXriAiIgI7d+7E7t27cffuXbi5ueG1116DSqWCQqHAzJkz4eTkhI4dO8LJyUnuUyGicojBjojoPyVNT7Jo0SKc/3ORTlsXFxd06NAB6enpqFSpEgCgb9++cpRNRCRhsCMiAjBn6yksPZCMZlZX8U6T1gAe3paNjIzEiXohcL19G/U0SQgJCUFwcDBatWrF0atEZHL4U4mIyh2tVotz584hOjpael2r2BQQWmw6/BsGtfLGK6+8AgCYM7AjwnaeRY1JkzClRxN5CyciegIGOyIyO/pOJpyeng4AcHV1BQD8/PPPGDhw4COfPIPatWvjlf794eHhIW319/fHcn9/Q58CEZFRMNgRkdkpmEwYgE64K3gmrl9TV/z44wGpN+7s2bMICwvDuHHjAABt2rSBg4MDWrVqhcDAQAQGBqJNmzZwd3eX43SIiAyGwY6IjMrQS3UBxU8m/OUfMfjhaAqyj63HV7NXFvlM4bnl6tSpA7VazWfkiMji8KcaEekwdBB7Uu/ahFIeS6PR4Pz58zh58qT0atKgK8IiIM07l3HoF6Qd/AX29vZ47rnndHrjPD09pWMpFAqGOiKySPzJRlRGjNFzZQyGCmIFnmWpLo1GgxEjRuDkyZM4ffo0srKydPa3vnkTNkHTpHnnlo7pgToLJ6Bx48awtrbWq04iIkvAYEdURgwdmIzFkGumFnfMwkt1jX7RD3FxcTq9cK6urlizZg0AQKVSYffu3bh8+TIAwN7eHo0aNULjxo3RuHFjXHYMwPr4bGky4VT3Fni7hf71ERFZCgY7ojJijMBkrF7AZ1kztThCCJ2lupRCi5UTe2NiXBxyc3N12laqVAlCCCgUCgDAl19+CTs7OzRu3Bi1atWSluYKj0zE+mImEy5cPxFRecNgR1SGDB2YjNkLWJo1UwsTQuDatWtITExEQkKCzn+zsrIw4ceIQkt1ARft/ZGbexyOjo5o1KgRmjRpIvXEFdavX78i31VcIC4uOBMRlTcMdkRlTN/A9KRjAYbtBSxQ3Jqp771UG7du3UJCQgKuXLmCPn36SO27dOmCXbt2FXssl7Zv69QUunw7/kI/jB8/Dp/1fA5KpVKv2jRaUez5FbzXaIWeZ0tEZBkUQgj+BHyCjIwMuLi4QK1Ww9nZWe5yzJ65DCIwloLgVRCYnjWAGeOYBcd7rZYVXK9FIzLFBldcmyDr73VI3bNKanf//n04ODgAAIYPH44VK1bAz88PdevWRZ06dVC3bl2cU/nhz4sajO9UB+93qlvkOwxx/kRElkyfHMIeO3oiU53+wtiMEUBLWmQeeLZbh/r0AmZmZiI5ORnJycm4fPlykV+/t2wbFu1NwoTOdXHy19n45scfATzsdXNt3w8uOTlwvRaNOnXqID09XQp2c+fOxZIlS4qMRl0QkYAJtYr+PrJ3jYjI8Bjs6IlMafqLkhgjhBn6vI35XNj/3zZVIFejxUer9uA5h9tSYJs1a5YUwMaNG4cf/wtrxUlLT5dq/N/1FwFA6oGLyXaFy4ufYWLXgCKfc3FxKfZ4j/t9Z08dEZFhMdjRE5Xl9BemNIjA0Of9LM+FaTQa3Lp1Czdu3EBKSgpu3LiBt99+Gw4ODlJNitNbcDniJ9i37Inf0Q/fHdgI9eHfAABjxoyB/3/rnVavXh0uLi6oXr06qlWrVuS/TZs2hb29PQCgf//+6N+/v1TH63qdMRERlTU+Y1cK5vSMnTGfXzPGs2F1P9ku3T5M+KrrMx2rpNuchnrezJDnXSArKws3btyQXi+//LLUs7ZkyRL8+OOPuHHjBlJTU6HVanU+e+bMGURct0JYRAKaWV3Dpq9GAACUSiWqhAyFsvErqKY+ieddMzBmzBhUrVoVwMOQWDBlCBERmT4+Y1eOmdL0F09S3KhLQ40QNUQvYOHjlua8tVot0tLScPv2bZ3XW2+9BScnJwDA8uXLsXLlSty+fRupqanIyMjQOcaZM2dQv359AMCdO3dw/PhxaZ9SqYSHhwe8vb3h5eUFhUIh9QL2qBmAD7seRNWqVVGlShVYWVn9F+TrFAnyDHVERJaLwc7ClPX0F8/yTJjcgwgeJysrC3fv3kV6ejpW/J2CXI0WKoVArkaL7pOX4H+T3kHFihUBAF999RUWLlyIu3fvFulVA4A2bdqgQYMGAICbN2/iyJEjOvvt7Ozg7e0Nb29vaDQaaXvv3r3x3HPPwcvLC97e3nB3dy8SyurV+/9f+/n5Ffm9ICKi8oXBzgIZo+fKkEGsLAYRWP83iGBhRDzGdX74bNnRo0fx77//Ij09XXqlpaVJv/7rr7/g5eUFAPjkk0+wYMECaSRo+oGfoT78G1zavo3T7fth3vY4fPVOOwAPe+pu374t1eDi4oLKlStLr8KLzb/xxhto3LixtM/b2xvOzs7SKguF+fv7S8/FERERlQaDnYUy5G1TQwexJw0iSFer8c8//+DevXs6r4yMDNy7dw9TpkyRnkNbuHAh1q1b93C/7/NAo+64d/hX3D2wBi5t38ZC9INS+fD8165diwULFpRY1507d6Rg5+rqCtfn+8ClXV8o47aidk4iKoeEoHJlFW7nnceak7Xh6fmwx3LIkCF4/fXXUblyZbi5ucHGxqbE7wgICEBAQNERpURERIbAYGehDHnbVKMVGN+pDoa1raqz/eXqClyua4uExERsVJ/CgwcPcP/+fdy/fx8PHjzAlClT/n9dz/BwREZG6uwP/+/X9+/fx7Vr1+Dk5CQFped++qnEekaOHCkFu8uXL+Pw4cMPe9YadZd61gAgI3otxr43VgqfzZo1w+uvv/4wtP33qlixovTrgsEFAPDxxx/DKfACVEolQmd1K/b3t2AUq4+PD3x8fJ7q95aIiMiQytWo2CVLlmDu3LlISUlBkyZNsGjRIrRq1eqJnzOFUbFarVZn2aUbN27g3r17yMnJQU5ODnJzc6Vfb72kxeZL/98rNiJ8E3Zet0Yzq2vwz7+I7OxsZGVlITs7G0IIrF69Wjrue++9h6ioKGRlZUmv7Oxs5OTkSHUU3DZ84403sGHDhhJrvnfvHhwdHQEAgwYNwsqVK0tse+3aNSkcffLJJ1i9ejWcnJyKfU2fPh3u7u4AgOPHj+PixYvYc8sB9na2GNjSA05OTnB1dYWTkxOUSmW5WM2CiIgslz45pNwEu7Vr12LAgAFYtmwZWrduLd3Ci4+Ph4eHx2M/W1bB7s0338SpU6ekgFb45eHhgWvXrkltAwMDizyED/z/6gCFb3V27twZxx5U1nlWrIBKpUJ+fr70/vXXX8emTZtKrDErKwt2dnYAgPHjx2Pr1q2oUKECHBwcUKFCBZ1fh4WFScFuz549OH/+fLHtKlSoAD8/P51n0YiIiOghBrtitG7dGs899xwWL14M4GHPk6+vL9577z1Mnjz5sZ8tq2D33HPP4Z9//il2X6VKlXQe0A8JCcGRI0dga2ur89IEdIGdjQ1i18yU2s6dOxdxcXG45toYKmtrNLe+AXt7e9jZ2cHe3h5jxoyReuFiY2ORlpYGe3t7nTYFv3Z0dCz2QX8iIiIyDga7R+Tm5sLBwQHr16/Ha6+9Jm0fOHAg0tPT8eeffz7282UV7I4dO4asrKwiYc3W1hZ2dnaoVKmS0b6biIiITBMnKH7E7du3odFo4OnpqbPd09MT586dK9K+4PZngUcnkTWW5557rky+h4iIiCyT8slNyp+ZM2fCxcVFevn6+spdEhEREdETlYtgV7lyZahUKty8eVNn+82bN6V5ywqbMmUK1Gq19Lpy5UpZlUpERET01MpFsLOxsUGLFi0QGRkpbdNqtYiMjERgYGCR9ra2tnB2dtZ5EREREZm6cvGMHQBMmDABAwcORMuWLdGqVSssXLgQ9+/fx6BBg+QujYiIiMggyk2w6927N27duoWpU6ciJSUFTZs2xY4dO4oMqCAiIiIyV+ViupNnZQorTxAREVH5pE8OKRfP2BERERGVBwx2RERERBaCwY6IiIjIQjDYEREREVkIBjsiIiIiC8FgR0RERGQhGOyIiIiILASDHREREZGFYLAjIiIishDlZkmxZ1GwOEdGRobMlRAREVF5U5A/SrNYGINdKdy7dw8A4OvrK3MlREREVF7du3cPLi4uj23DtWJLQavV4vr163BycoJCoZC7nFLLyMiAr68vrly5wjVuTQyvjWnj9TFdvDami9fGeIQQuHfvHnx8fKBUPv4pOvbYlYJSqUTVqlXlLuOpOTs7838yE8VrY9p4fUwXr43p4rUxjif11BXg4AkiIiIiC8FgR0RERGQhGOwsmK2tLaZNmwZbW1u5S6FH8NqYNl4f08VrY7p4bUwDB08QERERWQj22BERERFZCAY7IiIiIgvBYEdERERkIRjsTNjMmTPx3HPPwcnJCR4eHnjttdcQHx+v0yY7OxtjxoxBpUqV4OjoiF69euHmzZs6bZKTk9GtWzc4ODjAw8MDEydORH5+vk6bvXv3onnz5rC1tUXt2rWxcuVKY5+e2Sur67N3714oFIoir5SUlDI5T3NkqGsTGhqKFi1awNbWFk2bNi32u06ePIn27dvDzs4Ovr6+mDNnjrFOyyKU1bW5dOlSsf/fHDlyxJinZ/YMcX1OnDiBPn36wNfXF/b29qhfvz6++eabIt/Fv3eMg8HOhO3btw9jxozBkSNHEBERgby8PAQHB+P+/ftSm/Hjx2Pz5s1Yt24d9u3bh+vXr6Nnz57Sfo1Gg27duiE3NxeHDx/GqlWrsHLlSkydOlVqk5SUhG7duqFjx46IjY3FuHHjMHToUOzcubNMz9fclNX1KRAfH48bN25ILw8PjzI5T3NkiGtTYPDgwejdu3ex35ORkYHg4GBUr14dMTExmDt3LqZPn47vvvvOaOdm7srq2hTYvXu3zv83LVq0MPg5WRJDXJ+YmBh4eHjg559/RlxcHD755BNMmTIFixcvltrw7x0jEmQ2UlNTBQCxb98+IYQQ6enpwtraWqxbt05qc/bsWQFAREdHCyGE2LZtm1AqlSIlJUVq8+233wpnZ2eRk5MjhBDio48+Eg0aNND5rt69e4uQkBBjn5JFMdb12bNnjwAg0tLSyu5kLMzTXJvCpk2bJpo0aVJk+9KlS0XFihWlayWEEJMmTRL+/v6GPwkLZaxrk5SUJACI48ePG6v0cuFZr0+B0aNHi44dO0rv+feO8bDHzoyo1WoAgJubG4CH/yrKy8tDp06dpDb16tVDtWrVEB0dDQCIjo5Go0aN4OnpKbUJCQlBRkYG4uLipDaFj1HQpuAYVDrGuj4FmjZtCm9vb3Tu3BmHDh0y9ulYlKe5NqURHR2NF154ATY2NtK2kJAQxMfHIy0tzUDVWzZjXZsCr7zyCjw8PNCuXTv89ddfhim6HDHU9VGr1dIxAP69Y0wMdmZCq9Vi3LhxeP7559GwYUMAQEpKCmxsbODq6qrT1tPTU3r+KiUlRSc0FOwv2Pe4NhkZGcjKyjLG6VgcY14fb29vLFu2DBs2bMCGDRvg6+uLDh064N9//zXyWVmGp702pVGa60clM+a1cXR0xPz587Fu3Tps3boV7dq1w2uvvcZwpwdDXZ/Dhw9j7dq1GD58uLSNf+8Yj5XcBVDpjBkzBqdPn8bBgwflLoWKYczr4+/vD39/f+l927ZtceHCBSxYsAD/+9//DP59lob/75guY16bypUrY8KECdL75557DtevX8fcuXPxyiuvGPz7LJEhrs/p06fx6quvYtq0aQgODjZgdVQS9tiZgbFjx2LLli3Ys2cPqlatKm338vJCbm4u0tPTddrfvHkTXl5eUptHR5MVvH9SG2dnZ9jb2xv6dCyOsa9PcVq1aoXz588b6Aws17Ncm9J42utHxr82xWndujX/vyklQ1yfM2fOICgoCMOHD8enn36qs49/7xgPg50JE0Jg7Nix+OOPPxAVFQU/Pz+d/S1atIC1tTUiIyOlbfHx8UhOTkZgYCAAIDAwEKdOnUJqaqrUJiIiAs7OzggICJDaFD5GQZuCY1Dxyur6FCc2Nhbe3t4GPiPLYYhrUxqBgYHYv38/8vLypG0RERHw9/dHxYoVn/1ELFBZXZvi8P+bJzPU9YmLi0PHjh0xcOBAfPXVV0W+h3/vGJGsQzfosUaNGiVcXFzE3r17xY0bN6TXgwcPpDYjR44U1apVE1FRUeKff/4RgYGBIjAwUNqfn58vGjZsKIKDg0VsbKzYsWOHcHd3F1OmTJHaXLx4UTg4OIiJEyeKs2fPiiVLlgiVSiV27NhRpudrbsrq+ixYsEBs2rRJJCYmilOnTon3339fKJVKsXv37jI9X3NiiGsjhBCJiYni+PHjYsSIEaJu3bri+PHj4vjx49Io2PT0dOHp6Sn69+8vTp8+LX777Tfh4OAgli9fXqbna07K6tqsXLlS/PLLL+Ls2bPi7Nmz4quvvhJKpVL89NNPZXq+5sYQ1+fUqVPC3d1d9OvXT+cYqampUhv+vWM8DHYmDECxrxUrVkhtsrKyxOjRo0XFihWFg4ODeP3118WNGzd0jnPp0iXRtWtXYW9vLypXriw++OADkZeXp9Nmz549omnTpsLGxkbUrFlT5zuoeGV1fWbPni1q1aol7OzshJubm+jQoYOIiooqq9M0S4a6Ni+++GKxx0lKSpLanDhxQrRr107Y2tqKKlWqiFmzZpXRWZqnsro2K1euFPXr1xcODg7C2dlZtGrVSmeKDiqeIa7PtGnTij1G9erVdb6Lf+8Yh0IIIYzSFUhEREREZYrP2BERERFZCAY7IiIiIgvBYEdERERkIRjsiIiIiCwEgx0RERGRhWCwIyIiIrIQDHZEREREFoLBjoiIiMhCMNgRERERWQgGOyKiUhJCoFOnTggJCSmyb+nSpXB1dcXVq1dlqIyI6CEGOyKiUlIoFFixYgWOHj2K5cuXS9uTkpLw0UcfYdGiRahatapBvzMvL8+gxyMiy8ZgR0SkB19fX3zzzTf48MMPkZSUBCEEhgwZguDgYDRr1gxdu3aFo6MjPD090b9/f9y+fVv67I4dO9CuXTu4urqiUqVK6N69Oy5cuCDtv3TpEhQKBdauXYsXX3wRdnZ2WLNmjRynSURmSiGEEHIXQURkbl577TWo1Wr07NkTX3zxBeLi4tCgQQMMHToUAwYMQFZWFiZNmoT8/HxERUUBADZs2ACFQoHGjRsjMzMTU6dOxaVLlxAbGwulUolLly7Bz88PNWrUwPz589GsWTPY2dnB29tb5rMlInPBYEdE9BRSU1PRoEED3L17Fxs2bMDp06dx4MAB7Ny5U2pz9epV+Pr6Ij4+HnXr1i1yjNu3b8Pd3R2nTp1Cw4YNpWC3cOFCvP/++2V5OkRkIXgrlojoKXh4eGDEiBGoX78+XnvtNZw4cQJ79uyBo6Oj9KpXrx4ASLdbExMT0adPH9SsWRPOzs6oUaMGACA5OVnn2C1btizTcyEiy2EldwFERObKysoKVlYPf4xmZmaiR48emD17dpF2BbdSe/TogerVq+P777+Hj48PtFotGjZsiNzcXJ32FSpUMH7xRGSRGOyIiAygefPm2LBhA2rUqCGFvcLu3LmD+Ph4fP/992jfvj0A4ODBg2VdJhFZON6KJSIygDFjxuDu3bvo06cPjh07hgsXLmDnzp0YNGgQNBoNKlasiEqVKuG7777D+fPnERUVhQkTJshdNhFZGAY7IiID8PHxwaFDh6DRaBAcHIxGjRph3LhxcHV1hVKphFKpxG+//YaYmBg0bNgQ48ePx9y5c+Uum4gsDEfFEhEREVkI9tgRERERWQgGOyIiIiILwWBHREREZCEY7IiIiIgsBIMdERERkYVgsCMiIiKyEAx2RERERBaCwY6IiIjIQjDYEREREVkIBjsiIiIiC8FgR0RERGQhGOyIiIiILMT/AY0Fh5Gt6c2QAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(np.arange(1999, 2024), curve(np.arange(1999, 2024), *parameters), 'k--', label=f'Fit: {parameters[0]:.2f} * exp({parameters[1]:.2f} * (x - 2000))')\n",
    "plt.plot(years, articles, 'x', label='Data')\n",
    "plt.xlabel('Year')\n",
    "plt.ylabel('Number of articles')\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "\n",
    "# uncomment below to save the figure\n",
    "# plt.savefig(\"articles.png\", dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dropping a ball from the top of a tower\n",
    "\n",
    "We can use PyBaMM to simulate the simple physics-based model used to illustrate the blog post: dropping a ball from the top of a tower. Mathematically, if $x$ is the height of the ball, we can write the model as\n",
    "$$\n",
    "    x''(t) = -g, \\quad x(0) = x_0, \\quad \\text{ and } \\quad x'(0) = 0,\n",
    "$$\n",
    "where $g$ is the acceleration of gravity (in absolute value) and $x_0$ is the initial position of the ball (here $'$ denotes the time derivative). For more details on how the model is defined and solved, you can check the PyBaMM examples [Creating a simple ODE model](https://docs.pybamm.org/en/latest/source/examples/notebooks/creating_models/1-an-ode-model.html) and [ODE solver](https://docs.pybamm.org/en/latest/source/examples/notebooks/solvers/ode-solver.html).\n",
    "\n",
    "We first define the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m A new release of pip is available: \u001b[0m\u001b[31;49m23.1.2\u001b[0m\u001b[39;49m -> \u001b[0m\u001b[32;49m23.2\u001b[0m\n",
      "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m To update, run: \u001b[0m\u001b[32;49mpip install --upgrade pip\u001b[0m\n",
      "Note: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "%pip install pybamm==23.5 -q\n",
    "import pybamm\n",
    "\n",
    "# We define the model as a function so we can generate multiple instances of it\n",
    "def create_model():\n",
    "    model = pybamm.BaseModel()\n",
    "\n",
    "    # Define variables to solve for\n",
    "    x = pybamm.Variable(\"Position [m]\")\n",
    "    v = pybamm.Variable(\"Velocity [m.s-1]\")\n",
    "\n",
    "    # Define parameters\n",
    "    g = pybamm.Parameter(\"Gravity acceleration [m.s-2]\")\n",
    "    x0 = pybamm.Parameter(\"Initial position [m]\")\n",
    "\n",
    "    # Define equations and initial conditions\n",
    "    model.rhs = {x: v, v: -g}\n",
    "    model.initial_conditions = {x: x0, v: 0}\n",
    "\n",
    "    # Define event (collision with ground)\n",
    "    model.events.append(pybamm.Event(\"Collision with ground\", x))\n",
    "\n",
    "    # Define output variables (to plot)\n",
    "    model.variables = {\n",
    "        \"Position [m]\": x,\n",
    "        \"Velocity [m.s-1]\": v,\n",
    "        \"Acceleration [m.s-2]\": -g,\n",
    "    }\n",
    "\n",
    "    return model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can generate, process, and solve the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = create_model()\n",
    "\n",
    "param = pybamm.ParameterValues(\n",
    "    {\n",
    "        \"Gravity acceleration [m.s-2]\": 9.81,\n",
    "        \"Initial position [m]\": 100,\n",
    "    }\n",
    ")\n",
    "\n",
    "param.process_model(model)\n",
    "disc = pybamm.Discretisation()  # use the default discretisation\n",
    "disc.process_model(model, remove_independent_variables_from_rhs=False)\n",
    "solver = pybamm.ScipySolver()\n",
    "t_eval = np.linspace(0, 5, 100)\n",
    "solution = solver.solve(model, t_eval)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use PyBaMM's plotting functionality to easy plot the solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d0be606bdb6a4b8b8feb29b07bfc8982",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='t', max=4.5152327782515345, step=0.04515232778251534…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<pybamm.plotting.quick_plot.QuickPlot at 0x7f62b951a520>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solution.plot([\"Position [m]\", \"Velocity [m.s-1]\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can easily change the parameters, for example to simulate what would happen if we repeated the experiment in the moon, where gravity is $1.62 \\mathrm{m} \\; \\mathrm{s}^{-2}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_moon = create_model()\n",
    "\n",
    "param_moon = pybamm.ParameterValues(\n",
    "    {\n",
    "        \"Gravity acceleration [m.s-2]\": 1.62,\n",
    "        \"Initial position [m]\": 100,\n",
    "    }\n",
    ")\n",
    "\n",
    "param_moon.process_model(model_moon)\n",
    "disc = pybamm.Discretisation()  # use the default discretisation\n",
    "disc.process_model(model_moon, remove_independent_variables_from_rhs=False)\n",
    "solver = pybamm.ScipySolver()\n",
    "t_eval = np.linspace(0, 15, 100)\n",
    "solution_moon = solver.solve(model_moon, t_eval)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can equally easily plot both solutions together:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ecdbe6fae084469f8c323ac5710a2691",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='t', max=11.11110217382865, step=0.1111110217382865),…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot = pybamm.dynamic_plot(\n",
    "    [solution, solution_moon], \n",
    "    output_variables=[\"Position [m]\", \"Velocity [m.s-1]\"],\n",
    "    labels=[\"Earth\", \"Moon\"],\n",
    ")\n",
    "\n",
    "# uncomment below to save the figure\n",
    "# plot.fig.savefig(\"earth_moon.png\", dpi=300)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "env",
   "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.17"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}