{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explicit Trapezoidal Method\n",
    "\n",
    "Our goal is to implement the explicit trapezoidal method and compare its error to that of the forward Euler method.\n",
    "\n",
    "The explicit trapezoidal incorporates the forward Euler method ($\\widehat{y_{FE}}$) and is defined\n",
    "\n",
    "The method is\n",
    "\n",
    "\\begin{equation*}\n",
    "\\widehat{y}(t + \\Delta t) = \\widehat{y}(t) + \\Delta t \\cdot \\frac{f(t,\\widehat{y}(t))+f(t+\\Delta t, \\widehat{y_{FE}}))}{2}\n",
    "\\end{equation*}\n",
    "\n",
    "where \n",
    "\n",
    "\\begin{equation*}\n",
    "\\widehat{y_{FE}}(t + \\Delta t) = \\widehat{y}(t) + \\Delta t \\cdot f(t,\\widehat{y}(t))\n",
    "\\end{equation*}\n",
    "\n",
    "We first need the degradation model and the forward Euler method (for comparison)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.integrate\n",
    "\n",
    "def degradation_model( t, y, params ):\n",
    "    y1 = y[0]\n",
    "    y2 = y[1]\n",
    "    dydt = np.zeros( (2,) )\n",
    "    dydt[0] = -params[0] * y1\n",
    "    dydt[1] = -params[1] * y2\n",
    "    return dydt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "def forwardEuler( fn, t, y0, params ):\n",
    "    \"\"\" [t,y] = forwardEuler( fn, t, y0, params )\n",
    "    Given a function handle fn for dydt, find the value of\n",
    "    y at timepoints t, given y(0) = y0.\n",
    "    fn takes t, y, and params as input.\n",
    "    t is 1-dimensional (NT,)\n",
    "    y is 2-dimensional (NY,NT)\n",
    "    \"\"\"\n",
    "    y = np.zeros( (y0.size,t.size) )\n",
    "    y[:,0] = y0;\n",
    "    for i in range(1,t.size):\n",
    "        h = t[i]-t[i-1];\n",
    "        yprime = fn(t[i-1],y[:,i-1],params);\n",
    "        y[:,i] = y[:,i-1] + h*yprime;\n",
    "    \n",
    "    return (t,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implement Explicit Trapezoidal Method\n",
    "\n",
    "Implement the method so that is uses the same interface as Forward Euler.\n",
    "\n",
    "So you don't need to scan up, I retype the methor here:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\widehat{y}(t + \\Delta t) = \\widehat{y}(t) + \\Delta t \\cdot \\frac{f(t,\\widehat{y}(t))+f(t+\\Delta t, \\widehat{y_{FE}}))}{2}\n",
    "\\end{equation*}\n",
    "\n",
    "where \n",
    "\n",
    "\\begin{equation*}\n",
    "\\widehat{y_{FE}}(t + \\Delta t) = \\widehat{y}(t) + \\Delta t \\cdot f(t,\\widehat{y}(t))\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "def explicitTrapezoidal( fn, t, y0, params ):\n",
    "    \"\"\"[t,y] = explicitTrapezoidal( fn, t, y0, params )\n",
    "    Given a function handle fn for y', find the value of\n",
    "    y at timepoints t, assuming y(0) = y0.\n",
    "    fn takes t, y, and params as input.\n",
    "    t is 1-dimensional (NT,)\n",
    "    y is 2-dimensional (NY,NT)\n",
    "    \"\"\"\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run and plot result from Explicit Trapezoidal method.\n",
    "\n",
    "Run both the forward Euler and explicit trapezoidal methods on the degradation model with degradation rates (2,5) and initial conditions (1,1) with time steps from 0 to 4 with a step of 0.1.\n",
    "\n",
    "Then plot the true solution (solid lines), forward Euler's solution (dashed lines), and explicit trapezoidal's solution (dotted lines).\n",
    "\n",
    "Then plot the difference over time between the true solution and forward Euler's solution (dashed lines) and then the difference over time between the true solution and explicit trapezoidal's solution (dotted lines)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Error vs Step Size\n",
    "\n",
    "Determine how the error changes with the time step for both forward Euler and Explicit Trapezoidal. Plot two figures - one with the simulations and one with the errors. Compute the error using one method: Compute the Euclidean distance between the true solution and the approximate solution (2-norm of the difference) at each time point, then take the mean across time.\n",
    "\n",
    "Make an array of time steps (min=1e-4,max=0.2,numsteps=10). For each time step, plot each solution in a separate subfigure, and keep track of the error in an array. Then plot the errors in a second figure, with each error computation in a separate subfigure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'np' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-2-6cffbe5b90f0>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m# range from smallest timestep that doesn't take way too long to run to largest time step with qualitatively correct answer.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mdts\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinspace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1e-4\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0.3\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m \u001b[0merr1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mdts\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msize\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0merr2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mdts\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msize\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'np' is not defined"
     ]
    }
   ],
   "source": [
    "# range from smallest timestep that doesn't take way too long to run to largest time step with qualitatively correct answer.\n",
    "dts = np.linspace(1e-4,0.3,10)\n",
    "err1 = np.zeros( (dts.size,) )\n",
    "err2 = np.zeros( (dts.size,) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bonus Code\n",
    "This is the code Stephanie used to generate the illustration figure for her lecture notes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbcAAAEGCAYAAADmAds7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAABmsUlEQVR4nO3ddVxX1//A8dcBbEQUdTaimKRgze4WC2dO3abOmTB7xpybmzVn6+as2d31dbM7sVCxsANbDPL9++ODnx8qyIdQwvN8PO7Dz7333HPPAeTNufeEEhE0TdM0LSUxS+wCaJqmaVpC08FN0zRNS3F0cNM0TdNSHB3cNE3TtBRHBzdN0zQtxbFIrBtnzZpV8ufPn1i31zRNS5aOHj16X0SyJXY5krpEC2758+fnyJEjiXV7TdO0ZEkpdTWxy5Ac6MeSmqZpWoqjg5umaZqW4ujgpmmapqU4OrhpmqZpKY4ObpqmaVqKE2NwU0rNUkrdU0qdjua8UkpNVEpdVEqdVEq5JXwxNU3TNM10prTc5gB13nO+LlAoYusMTIt/sWJPRNi+fTuDBg0iPDw8MYqgaZqmJRExBjcR2QU8fE+SRsA/YnAAsFZK5UyoAr5t/u5NVBj7NXce3wfg7t27jB49miJFitCjRw/y5s2LUupD3T7FmTt3LoUKFaJQoULMnTvX5Ot27dqFm5sbFhYWLF++PNp0R48excnJCXt7e3r27MnrJZYePnxIzZo1KVSoEDVr1uTRo0fxroumadprypT13JRS+YH1IuIYxbn1wEgR2ROx/x/QX0TeGaGtlOqMoXVHvnz53K9ejf1YxBrjv+Nu5j1ISFrqZm7IuC9+QimFnZ0dBQoU4LPPPiN79uzv/Pt6S506dazvmVI9fPiQkiVLcuTIEZRSuLu7c/ToUTJnzhzjtf7+/jx9+pSxY8fi4eGBp6dnlOlKly7NxIkTKVOmDPXq1aNnz57UrVuXfv36kSVLFgYMGMDIkSN59OgRo0aNSugqalqKo5Q6KiIlE7scSd1HnaFERP4C/gIoWbJknFZJ/ddrGl6zp7Dl1WY2By6j0MDihKx+QVqVFg8PDywtLbl79y53797l5MmT3Lt3j3v37nH37l0CAgKwtLR8I+BFFQRff7ayskoxrcChQ4eSJUsWvLy8ABg0aBCpU6emZs2aZMmSBYCaNWuyefNmWrVqFWN+r6dOMzOLvvF/+/Ztnj59StmyZQFo164dq1evpm7duqxZs4YdO3YA0L59e6pUqaKDm6ZpCSYhgttNIG+k/TwRxz6Y8V91w/dKY1rMmIZZviOk632D8P3C4B8H8/Own+nXr1+U14WHh/Po0aM3At7rf48fP/7GsXv37hEcHPxOwIsuIGbLlg0Li0SbzSxGX3/9NU2bNsXLy4vw8HAWL15M586dyZv3/791efLk4eZNw7euRYsWnD9//p18vv/+e9q1a2fSPW/evEmePHmizP/u3bvkzGl4ep0jRw7u3r0b57ppmqa9LSF+G68FuiulFgNlgCcicjsB8n2v4na58fl5OF0nLuPfgAukcd9KDqcwNltspsydMpTKUeqda8zMzLCxscHGxoZixYrFeI8XL14QEBDwRsC7e/cu169f58iRI28ce/jwIZkyZYqxNfj6c4YMGT5qqzB//vzY2Nhw/Phx7t69S4kSJTA3NyckJCTK9EuWLPloZVNKpZgWsqZpSUOMwU0ptQioAmRVSt0AfgRSAYjIdGAjUA+4CLwAvvpQhX2bubkZf3q3YPux83w7JzeSy5+rWTbx9ZavKZ+7PL1K9KKYTcxBLDrp06fH1tYWW1vbGNOGhYXx8OHDN1qDr4PfoUOH3giEr1sp7wt+kf+1sbHB3Nw8zvV4rWPHjsyZM4c7d+7w9ddf8+TJE+OjQYAbN25QpUoVIGFabrlz5+bGjRtv5J87d27AUPfbt2+TM2dObt++Tfbs2eNeMU3TtLeY1KHkQyhZsqQk5KoAL4NCaD1qKcdfZsAizSas7Y4TxEvq5q9L9xLdyWeVL8HulRACAwPfeTz69qPS18ceP35MlixZTHpPmD17dtKnTx/lPYODg3FyciIkJIQLFy7w5MkT3N3dOXbsGABubm4cPXrU+A7OFB06dKBBgwYmdyjp0aMH9erVo2/fvtjY2Bg7lDx8+JDRo0fH/gupaZ8Y3aHENCkmuL02/9+jDNlwgfBUQo5Mqwi3vUJoeChNCzXlW5dvyZ4++bUQQkNDuX//frTB7+3P1atXZ8OGDVHm1aVLF6ytrRk5ciQAs2bN4tdffwUMnUy++sq0hvfhw4dp0qQJjx49Im3atOTIkYMzZ84A4Orqio+PDwBHjhyhQ4cOvHz5krp16zJp0iSUUjx48IAvvviCa9euYWtry9KlS2MVVDXtU6WDm2lSXHADuPP4Oc1Hr+J6eGZ4cIyqNW7iE3QQCzML2hRrw9dOX2OV2uqD3DuxiQjh4eFRPsYMDw/Hzc2NZcuWUahQoUQonaZp8aWDm2lS5NySOawzsGtEG751ywiZHPh3X2XynqlN5dyVmXl6JnVX1GXW6Vm8DH2Z2EVNcEqpKAObr68v9vb2VK9eXQc2TdNSvBTZcovM9/p9Wk36H0/MMmHuf5Dh3zqxN2w7u2/uJnu67HRx7UJj+8akMkv1wcuiaZoWX7rlZpoU2XKLrHjerBz+tSX18psTaluKfosfo3bZMaPGDHJZ5mL4/uE0Xt2YzVc2Ey56TkpN07SUIMUHN4DUFmZM7VKHWW2cscyYifWvivC11zL62w5kUrVJpDZPTd9dfWm5viX7bu4jsVqzH9rEiRMpVqwYbdq0YfXq1QwfPvy96aOb/3H9+vUMHTo0zuWYM2cO2bJlw9XVFVdXV/7+++8o0+l5KTVNi6tPIri9Vt3ZlkM/N8Etq/CiYDVqjdrIrqVnWFp/Kb9W+JWnwU/59t9v6fi/jpwMOJnYxU1wU6dOZevWrSxYsIDRo0fTtWvX96YfOXIk1atX58KFC1SvXt3Yw7J+/fqsW7eOFy9exLksLVq0wMfHBx8fHzp27Bhlmu+++44ZM2Zw4cIFLly4wObNm99bLk3TtNc+qeAGkCldKlb1a8SwOvlJ+1kB/r6RjfKtvHBJ5cLaxmsZUHoAFx9fpM3GNnht9+LS40uJXeQE0aVLFy5fvkzdunUZNWoUadKkIWvWrDx79gw7OzvjTCVPnz417q9Zs4b27dsDhvkfV69eDRg6rVSpUoX169d/sPJGnpdSKWWclxKItlyapmmvfXLB7bUOVRzYObA2tplSca9gPcp6T2fOvKW0LtqajU030s21GwduH6Dp2qYM2TuE24EffEaxD2r69OnkypWL7du3kz17dtzcDGvKZsyYkSpVqhjHxS1evJimTZuSKlWq987/WLJkSXbv3g3AggULjI8YI2/RDewGWLFiBc7Oznh6enL9+vV3zut5KTVNi49PNrgB5M2SgR1DG/N1yWyktv+c4QdDqN22Gy8ev6CLSxc2Nd1E22Jt2Xh5I/VX1Wf04dE8fPW+pe2Sh9u3b5MtWzbjfseOHZk9ezYAs2fPjnIg99vzP2bPnp1bt24B0KZNG+MjxshbdOu8NWzYEH9/f06ePEnNmjWNrbC40PNSapoWlU86uAGYmymGepZm5XflsMmcifN56lKi/Y+sWrOOzGkz07dUX9Y3WU+DAg1YcHYB9VbWY9qJaTwPeZ7YRY+zdOnS8erVK+N++fLl8ff3Z8eOHYSFheHoaFi27/X8j8A78z++evWKdOnSATG33AYNGmQ8BmBjY0OaNGkAQ2A9evToO2U0ZV7KqMqlaZoGOrgZueW3Yd+PHtS2z4iFc326rrxEq2+9ePbsGTktczK8/HBWeazi85yfM9VnKvVW1mPB2QUEhwUndtFjrVixYly8ePGNY+3ataN169ZvtNo8PDyMq3PPnTuXRo0aGc/5+fkZg2BMLbcRI0YYjwHGwASwdu3aKFdoyJkzJ1ZWVhw4cAAR4Z9//jHe/33l0jRNAwzTNSXG5u7uLknVuuPXpFD/VZLXe7nYVmsjO3bseOP8iXsn5OvNX4vjHEeptayWrLm4RkLDQhOptKaztbWVgIAAef78uRQvXlzCw8ON527fvi1p06aVR48eGY/dv39fqlWrJvb29lK9enV58OCB8Vz9+vXl5MmTcSrHgAEDpHjx4uLs7CxVqlSRs2fPGs+5uLgYPx8+fFgcHBykQIEC0q1bN2N531cuTUvpgCOSSL+3k9OW4mcoiau7T1/xzV87OX0/lBcXDtA0z0t+//UnMmTIABj+KNh/ez/jj47n7MOz2Fvb07NET6rkrZIs3gH16tWLhg0bUqNGDQCWL1/OmjVrmDdvXozX3r17l9atW/Pff/996GJqmvYWPUOJaXRwe4/wcOHP7ecZ878LBD9/jMXRRfz98/fGgAAQLuFsvbqVSccncfXpVVyyueDl5kXJHEn7Z+/u3bscPHgQDw8PevTowaZNm9i4cSOFCxeO8drDhw+TKlUq4zs0TdM+Hh3cTKODmwnO3n5Kx5l7uRkYztMja2mUX/hj7Gisra2NaULCQ1hzcQ3TfKZx7+U9KuSuQC+3XhTNUjTxCq5pWoqjg5tpdIcSExTLacV//WvxZek8WJX0YIu44FC+FmvXrjWmSWWWCs/CnmxouoHv3b/nZMBJmq9rTr9d/bj29Foill7TNO3To1tusbTTL4BeC4/w+EUwD7fPpp59OiZNnPjGuDGAp8FPmXN6DvN85xEaHkqzws341vlbsqXPFk3OmqZpMdMtN9PollssVS6cjf/6VKN68Zxkqd6JHeKIQ8lyLFq06I0Jl61SW9HTrScbm26kWeFmrPBbQb2V9ZhwbAJPg58mYg00TdNSPt1yiyMRYeGhawxfe4bgl8+5u/4PahTNyrRp04yDjSO79vQak30ms+nKJqxSW9HRqSOtirYirUXaRCi9pmnJlW65mUa33OJIKUWbMrZs9KqEQ/4cZG86iIPhBSnuXIIZM2a8s2xOPqt8jK40mmUNl+GczZlxR8dRf2V9lvstJzQ8NJFqoWmaljLpllsCCA4NZ/y/fkzbcQnzV4+4vmQ4FYrnZcaMGRQoUCDKaw7fOcyEYxM4EXCC/Fb56V6iOzVta2Km9N8bmqZFT7fcTKN/kyaA1BZm9KtTlEWdy5Lts5zkbj+OU+G5cXJ2Zvz48YSFhb1zTakcpZhXdx4Tq07EwsyCPjv70GpDK/bdSrmLpWqapn0suuWWwJ68DGHw6tOsO3GL9M9v4ffPYEo7FGTmzJlRzqEIEBYexoYrG5hyfAq3nt+iTI4y9HLrhVM2p49cek3TkjrdcjONbrklsEzpUjGxpSt/tHDBLEs+CnabyaXQLLi6ujJixAjjoqCRmZuZ41HQg3VN1jGg9AAuPL5A642t8d7uzeXHlxOhFpqmacmbbrl9QNcfvsB7iQ9Hrj7CJvAyx/8eiEuxQsyaNYsSJUpEe93zkOf84/sPc8/M5WXoSxoVbMR3Lt+R0zLnRyy9pmlJkW65mUa33D6gvFnSs7hzWXrXLMxjq4K49FvMXclIqVKl+OGHH95YUy2yDKky8J3Ld2xsupE2xdqw/vJ6GqxqwJjDY3j06tFHroWmaVryY1LLTSlVB5gAmAN/i8jIt87bArOAbMBDoK2I3Hgno0g+hZZbZMevPcJriQ/XHrwg3/Nz7Jran6KFCzFz5kzKlSv33mtvB95m6omprL20lnQW6ejg0IF2xduRPlX6j1R6TdOSCt1yM02MLTellDkwBagLFAdaKaWKv5VsLPCPiDgDw4HfErqgyV2JfJnZ2LMiX5TMy9UMRak0fC0vLTJSoUIFevXqRWBgYLTX5rTMyc/lf2alx0rK5izLFJ8p1F1ZN9kulqppmvahxdhyU0p9DgwTkdoR+wMBROS3SGnOAHVE5LoyLGb2RESs3pfvp9Zyi2zz6dsMWHmKoJAwCgWeZN34AeTPn58ZM2a8sZxOdE4EnGDCsQkcvnOY3Ja56ebajXp29TA3M/8Ipdc0LTHplptpTHnnlhu4Hmn/RsSxyE4ATSM+NwEyKqVs3s5IKdVZKXVEKXUkICAgLuVNEeo45mRzr0q422bhZBpHmo7bTCrLzNSsWZNvvvmGx48fv/d6l2wuzKw1kz9r/IlVait+2PMDnus82XF9hx4jp2maRsJ1KOkDVFZKHQcqAzeBd0Yui8hfIlJSREq+PYv+pyZHprT883VpBtcvxqn74Vi1GEX7/iOZO3cuxYsXZ8WKFe8NVEopyuUux+IGixlTeQzBYcH02NaDdpvacfTu0Y9YE03TtKTHlOB2E8gbaT9PxDEjEbklIk1FpAQwKOLY44QqZEplZqboWLEAq7uVJ0uG1OzAkW//3Er2nLnx9PSkQYMGXLly5f15KDPq5K/D6sarGVJ2CDcDb9Jhcwe6/tuV8w/Pf6SaaJqmJS2mvHOzAPyA6hiC2mGgtYiciZQmK/BQRMKVUiOAMBEZ+r58P+V3blF5FRLGyE3nmLPPn8KfWVI69AwThvcnLCyMIUOG0Lt3b1KnTh1jPi9DX7Lo3CL+PvU3gcGB1LWrS3fX7uS1yhvjtZqmJX36nZtpYmy5iUgo0B3YApwFlorIGaXUcKWUR0SyKsB5pZQf8Bkw4gOVN8VKm8qcYR4OzP6qFA+fh7D0iR0/Ld5F3Xr1+OGHH3B1dWXXrl0x5pPOIh1fO37Npqab+MbpG7Zd24bHag9+OfAL91/e/wg10TRNS3x6hpIk6H5gEANWnOTfs/eoWCgrdbPcZ9D33fH396dDhw6MHj36nZW/o3PvxT3+PPEnKy6sILV5atoWa8tXjl+RMXXGD1wLTdM+BN1yM40ObknU68VQf17vS7pU5vzUoAgHlv/JmDFjsLKyYvTo0Xz11VeYmZnWJ+jq06tMOT6FTf6byJQmEx0dO9KyaEu9WKqmJTM6uJlGT7+VRL1eDHV9j4rkzpyOnktPIyVbcuDwMRwcHOjYsSOVKlXi1KlTJuVna2XL6MqjWdpgKY5ZHfn96O/UX1WfFX4r9GKpmqalOLrllgwEh4Yzbqsff+66RH6bDPzxhQs+29bQp08fnjx5wvfff8/QoUPJkCGDyXkevnOY8cfGczLgJPmt8tOjRA9q2tbEMAZf07SkSrfcTKNbbslAagszBtQtysKOZQkKCcNz+n6e21bA9+w52rdvz+jRoylevDhr1641Oc9SOUoxv+58JlSdgLkyp/fO3rTc0JL9t/Z/wJpomqZ9HLrllsw8eRHCoNWnWH/yNqXzZ2FcCxeu+h6nS5cunD59mkaNGjFx4kTy5ctncp5h4WGsv7yeKT5TuP38NmVylsHLzQvHrI4fsCaapsWFbrmZRrfckplM6VMxqVUJxn3hgu/tp9SdsJv7GfJz7NgxRo8ezdatWylWrBhjx46NcmHUqJibmdPIvhHrm6ynf6n++D30o9WGVny/43suP0kei6VWqVKFxPhjKSgoiBo1auDq6sqSJUs++v3fp0OHDixfvhyAjh074uvrm6jluXr1Km5ubri6uuLg4MD06dOjTNeiRQtcXV1xdXUlf/78uLq6AuDv70+6dOmM57p06fIRS68lOyKSKJu7u7to8XPtwXNpOnWv2PZfLz0XHZMnL4PF399fPDw8BBBHR0fZs2dPrPMNDA6UqcenSun5pcV5rrMM3TtUbgfe/gA1SDiVK1eWw4cPf/T77t+/X6pXrx7n60NCQhKwNG9q3769LFu27IPlH1tBQUHy6tUrERF59uyZ2Nrays2bN997zffffy8//fSTiIhcuXJFHBwcPng5kzrgiCTS7+3ktOmWWzKWN0t6lnQuy/c1C7P+5G3qjt/N3fCMrFmzhtWrV/PkyRMqVKhAp06dePDggcn5ZkiVge9cv2NTs020LtqadZfWUX9lfcYeHsvjV48/XIVM8Pz5c+rXr4+LiwuOjo5RtpYWLVqEk5MTjo6O9O/f33jc0tISb29vHBwcqF69Oq8n77506RJ16tTB3d2dihUrcu7cuXfyfPjwIY0bN8bZ2ZmyZcty8uRJ7t27R9u2bTl8+DCurq5cunTpjWsOHz6Ms7Mzrq6u9O3bF0dHw2PeOXPm4OHhQbVq1ahevTrPnz/n66+/pnTp0pQoUYI1a9YAEBYWRt++fSlVqhTOzs78+eefAOzYsYMqVarg6elJ0aJFadOmTYwTZkdu2VpaWjJo0CBcXFwoW7Ysd+/eBSAgIIBmzZpRqlQpSpUqxd69e036npgqderUpEmTBjC0eMPDw9+bXkRYunQprVq1StByaJ+IxIqquuWWsI5efSiVRm8TuwHrZfTmsxIcGibPnj2Tvn37irm5uWTNmlXmzJkj4eHhsc775rObMmj3IHGe6yxlF5SV6T7T5Xnw8w9Qi5gtX75cOnbsaNx//PixiPx/y+3mzZuSN29euXfvnoSEhEjVqlVl1apVIiICyPz580VE5KeffpJu3bqJiEi1atXEz89PREQOHDggVatWfee+3bt3l2HDhomIyH///ScuLi4iIrJ9+3apX79+lGV1cHCQffv2iYhI//79ja2O2bNnS+7cueXBgwciIjJw4ECZN2+eiIg8evRIChUqJIGBgfLnn3/Kzz//LCIir169End3d7l8+bJs375drKys5Pr16xIWFiZly5aV3bt3v3P/yC23yC1bQNauXSsiIn379jXeo1WrVsZ8rl69KkWLFn0nz3PnzomLi0uU26NHj6L8OkR27do1cXJyknTp0snkyZPfm3bnzp0S+ffElStXJH369OLq6iqVKlWSXbt2xXi/lAjdcjNp08EtBXn2KkT6LvMR2/7rxWPSbrkcECgiIidOnJBy5coJIJUrVxZfX9845X/h4QXp+V9PcZzjKJUWV5IFvgskODQ4IasQo/Pnz4utra3069fvjV9ur395r169Wr788kvj8b///lu8vb1FRMTMzMz4GPDSpUvi4uIiz549k7Rp077xSzqqX+qurq5y6dIl436ePHnkyZMn0Qa3R48eSb58+Yz7J06ceCO4dejQwXjO3d1dHBwcjPfPmzev+Pr6SrNmzaRQoULG4/nz55ctW7bI9u3bpUaNGsbru3TpYgyOkUUX3FKnTm38I2fx4sXyzTffiIhItmzZ3vg65MqVS549exb1NyKebt68KaVKlZI7d+5Em6ZLly4yduxY4/6rV6/k/v37IiJy5MgR4/fgU6ODm2mbfiyZglimsWC0pwvT2rjh/+AF9SbsZvGhazg5ObF7925mzJjByZMncXFxYdCgQbx48SJW+dtntmdCtQnMrzefApkK8Nuh32i4uiHrLq0jLPydFY4+iMKFC3Ps2DGcnJwYPHgww4cPj3NeSinCw8OxtrbGx8fHuJ09ezYBSxy1yGMSRYQVK1YY73/t2jWKFSuGiDBp0iTj8StXrlCrVi0A4+M9AHNzc0JDTR+InypVKuN4xsjXhoeHc+DAAeP9bt68iaWl5RvXnj9/3tih4+3t7XUIDx48aDz39jCVXLly4ejoyO7du6MsY2hoKCtXrqRFixbGY2nSpMHGxrBMpLu7OwULFsTPz8/kemufFh3cUqC6TjnZ7FURN1trBqw8xbfzjvL4ZSgdO3bk/PnztG7dml9//RVHR0c2btwY6/xdsrkwq/YspteYblwstfn65uy8vtPwOOADunXrFunTp6dt27b07duXY8eOvXG+dOnS7Ny5k/v37xMWFsaiRYuoXLkyYPjl/br34MKFC6lQoQJWVlbY2dmxbNkywBBoTpw48c59K1asyIIFCwDDO6+sWbNiZRX9YvPW1tZkzJiRgwcPArB48eJo09auXZtJkyYZv3bHjx83Hp82bZqx16ufnx/Pnz+P+YsUR7Vq1WLSpEnGfR8fn3fSFClS5I0/BCJv1tbWb6QtU6aM8ZyHhwc3btzg5cuXADx69Ig9e/ZQpEiRKMvy77//UrRoUfLkyWM8FhAQQFiY4Y+oy5cvc+HCBQoUKBDPWmsplQ5uKVTOTOmY93UZBtUrxo7zAdQZv4tdfgFky5aNOXPmsH37dtKkSUP9+vXx9PTkxo0bscpfKUX53OUNi6VWGkNQaBDdt3Wn/eb2HLt7LOYM4ujUqVOULl0aV1dXfvrpJwYPHvzG+Zw5czJy5EiqVq2Ki4sL7u7uNGrUCDC0lg4dOoSjoyPbtm1j6FDDqkwLFixg5syZuLi44ODgYOzQEdmwYcM4evQozs7ODBgwgLlz58ZY1pkzZ9KpUydcXV15/vw5mTJlijLdkCFDCAkJwdnZGQcHB4YMGQIYuu8XL14cNzc3HB0d+fbbb2PVQoutiRMncuTIEZydnSlevHi0XfXj6uzZs5QpUwYXFxcqV65Mnz59cHJyAgx1jTyUY/Hixe90JNm1a5exg46npyfTp08nS5YsCVpGLeXQg7g/AWduPaHXYh8u3gvk6/J29KtThLSpzAkODmbs2LH8/PPPmJmZ8cMPP9C7d2/Spo39ZMoh4SGsurCKaSemcf/lfSrlqUTPEj0pkiXqv8wTg6WlJYGBgR/tfoGBgcbHeiNHjuT27dtMmDDho91fS5n0IG7T6OD2iXgVEsZvG88yd/9ViubIyISWJSiSw7Dsjb+/P3369GHFihXY2dkxbtw4GjVqFKd5Jl+GvmTB2QXMOj2LwOBA6hWoRzfXbuTNmPiLpX7s4LZkyRJ+++03QkNDsbW1Zc6cOSYvVaRp0dHBzTQ6uH1itp+7R9/lJ3j6KpQBdYrSoVx+zMwMQWzbtm307NmTM2fOUKNGDSZMmEDx4sXjdJ8nQU+YfXo2C84uIFRC8Szkybcu35I1XdaErI6mfXJ0cDONDm6foPuBQfRffpL/zt2jUuFsjPV0JruV4VFkaGgo06dPZ8iQITx79ozu3bszbNiwdzoLmEovlqppCUsHN9Po4PaJEhHmH7zGiA2+pE9twcimTtRyyGE8f//+fYYMGcKff/6JjY0NI0aM4JtvvsHc3DxO97v69CqTj09ms/9mMqXJRCenTrQo0kIvlqppsaSDm2l0cPvEXbwXSK/Fxzlz6ymtSudjSINipE9tYTzv4+NDz5492b17NyVKlGDixIlUqFAhzvfzfeDLxGMT2XtrL5+l/4yurl3xKOiBhZlFzBdrmqaDm4n0UIBPnH12S1Z1Lc+3lQuw+PA1Gkzcw8kbj43nXV1d2blzJ4sXLyYgIICKFSvSunXrWA8deK24TXGm15zOzFoz+Sz9Z/y470earGnC1qtbP/gYOU3TPh265aYZ7bt0n95LTxDwLAjvmoXpUrkg5mb/32Py+fPnjBo1itGjR2Nubh6voQNgeDS67fo2Jh6byOUnl3GwccDL3YuyOcsmVJU0LcXRLTfT6OCmveHJixB+WH2KDSdvU9ouC3+0cCW3dbo30iTk0AEwLJa67vI6pvhM4c7zO5TNWRYvNy8csjokRJU0LUXRwc00Orhp7xARVh67ydA1pzEzU/zS2JFGrrnfSZeQQwcAgsKCWHp+KTNOzuBR0CNq2takR4ke2GWyi091NC1F0cHNNDq4adG69uAFXkuOc+zaY5qUyM1PjRywSpvqjTShoaFMmzaNoUOHJsjQAYDA4ED+8f2HuWfmEhQWRGP7xnRx6UKODDlivljTUjgd3Eyjg5v2XqFh4UzZfomJ2y6Qwyot41u6Uir/u/P5JfTQAYAHLx/w96m/WXJ+CQpF62Kt+cbxG6zTWsejRpqWvOngZhod3DSTHLv2CK/FPtx49IJuVe3pWb0Qqczf7Wyb0EMHAG4G3mSqz1TWXVpHhlQZ+MrxK9oWa0v6VOnjla+mJUc6uJlGDwXQTOKWLzMbe1WkmVseJm27iOf0/Vy5/+7yKwk9dAAgt2VuRlQYwQqPFZTKUYpJxydRb2U9Fp1bREhYSHyqpWlaCmVScFNK1VFKnVdKXVRKDYjifD6l1Hal1HGl1EmlVL2EL6qW2CzTWDCmuQtT27jhf/859SfuZsnha++MT1NK0aJFC86dO8eQIUNYuXIlRYoUYcSIEbx69SrO9y+UuRATq01kXt155M+Un18P/krD1Q1Zf3k94RIe3+ppmpaCxPhYUillDvgBNYEbwGGglYj4RkrzF3BcRKYppYoDG0Uk//vy1Y8lk7fbT17Se+kJ9l16QG2HzxjZ1JnMGVJHmfbKlSv06dOHlStXYmdnx++//07jxo3jPHQADD06997ay4RjEzj38ByFMxeml1svKuauGK98NS2p048lTWNKy600cFFELotIMLAYaPRWGgFeL0ucCbiVcEXUkqKcmdIx/5sy/FCvKNvO3aPOhF3suXA/yrR2dnasWLGCf//9l/Tp09O0aVNq1KhhXHE6LpRSVMhdgSUNljCq4ihehLyg23/d6LC5A8fvxT1fTdNSBlOCW27geqT9GxHHIhsGtFVK3QA2Aj2iykgp1VkpdUQpdSQgICAOxdWSEjMzRedKBVndrTwZ06ai7cyD/Lzel1chYVGmr169Oj4+PkycOJETJ07g7u5O+/btuX79epTpTSqDMqNegXqsbbyWwWUGc+3ZNdptakeP/3rg98gvzvlqmpa8JVSHklbAHBHJA9QD5iml3slbRP4SkZIiUlIv2phyOOTKxLruFWj3uS0z91yh8ZS9+N19FmVaCwsLevTowcWLF+nbty9LliyhcOHCDBw4kCdPnsS5DKnMU9GiaAs2NNlAL7deHL17FM+1ngzcPZAbz+LemUXTtOTJlOB2E4i8jHKeiGORfQMsBRCR/UBaQK9K+QlJl9qc4Y0cmdWhJPcDg2gwaQ9z9l6JdjJka2trRo0axfnz5/H09GTkyJHY29szefJkQkLi3gMyfar0dHTqyKZmm+jg2IGtV7fScHVDfj34K/dfRv3YVNO0lMeU4HYYKKSUslNKpQZaAmvfSnMNqA6glCqGIbjp546foGpFP2OzVyUq2Gdl2DpfOsw+zL1n0feQtLW1Zd68eRw5cgQnJyd69OiBg4MDK1eujNcqAZnSZOJ79+/Z0GQDTeybsPT8UuqtrMek45N4Fhx1q1LTtJTDpEHcEV37xwPmwCwRGaGUGg4cEZG1ET0kZwCWGDqX9BOR/70vT91bMmV7vRjqL+t9yZDGglHNnKlZ/LMYr9m4cSP9+vXD19eX8uXLM3bsWMqWjf8qAf5P/JnsM5kt/luMi6W2LNqSNOZp4p23pn1MurekafQMJdoHdfHeM3ot9uHMrae0LpOPwfXfXAw1KqGhocyaNYuhQ4dy9+5dmjdvzm+//UbBggXjXZ4zD84w8dhE9t3ax2fpP6ObazcaFmyoF0vVkg0d3Eyjg5v2wQWHhvP71vP8tesydlkzMKFFCZzyZIrxusDAQMaOHcuYMWMICQmha9euDBkyBBsbm3iX6eDtg0w4NoFT909hl8mOniV6Uj1fdT1GTkvydHAzjQ5u2kez79J9vl9ygvuBQXxfqzDfVnpzMdTo3Lp1ix9//JFZs2aRMWNGBg0aRI8ePeK8SOprIsK2a9uYeNywWKqjjSNe7l6UyVkmXvlq2oekg5tpdHDTPqrHL4IZtOo0G07dpoxdFsZFsRhqdE6fPk3//v3ZuHEjtra2jBgxglatWmFmFr8RLaHhoay7tI6pJ6Zy5/kdPs/5Ob3ce+FgoxdL1ZIeHdxMo4Ob9tGJCCuO3eTHiMVQRzRxwsMll8nX//fff/Tp0wcfHx/c3d0ZM2YMVatWjXe5gsKCWHJuCTNOzeBx0GNq2daie4nuerFULUnRwc00Orhpiebqg+d4LfHh+LXHNI1YDDXjW4uhRic8PJwFCxYwaNAgrl+/ToMGDRg1alS8VgJ/LTA4kLm+c5l7Zi7BYcF6sVQtSdHBzTQ6uGmJKjQsnMnbLzJp20VyZkrL+BaulIxiMdTovHz5kgkTJvDbb78RGBhIx44d+emnn8iRI/6B6MHLB8w4NYMl55dgrsxpXbQ13zh9Q6Y0MXeG0bQPRQc30+jgpiUJR68+wnuJYTHU7lXt6RHNYqjRCQgI4Oeff2batGmkSZOGfv360bt3bzJkyBDvskVeLNUylSVfOX5Fm2Jt9GKpWqLQwc00OrhpSUZgUCjD1p5h+dEbuOS1ZkILV/JnjV1w8vPzY+DAgaxcuZKcOXMyfPhwvvrqK8zNzeNdPr9Hfkw6Pokd13dgk9aGLi5daFaoGanMTXuUqmkJQQc30+jgpiU5G07e5odVpwgJC2dYQweal8wT6/Fne/fupU+fPhw4cABHR0dGjRpF3bp1E2Qcm889H/44+gfH7h0jj2UeupfoTl27upi9O1e4piU4HdxMo/83aklOfeecbPaqiEsea/qtOMl384/x6HlwrPIoX748+/btY+nSpbx48YL69etTuXJldu/eHe/yuWZ3ZU6dOUytPpUMqTIwYPcAvlj3Bbtu7IrXfJiapiUcHdy0JClnpnQs6FiGgXWL8t+5u+9dDDU6SimaN2/O2bNnmTJlChcuXKBSpUrUrVuXo0ePxqt8Sikq5qnI0oZLGVVxFM9DnhsXS/W55xOvvDVNiz/9WFJL8k7ffILXEh8u3gukU0U7+tQuQhqL2L9De/HiBZMnT2bUqFE8fPiQZs2aMXz48AQZPhASFsKKCyuYfmI6D149oEreKvQs0ZNCmQvFO29Ni0w/ljSNDm5asvAyOIxfN55l3oGrFMtpxYSWrhT+LGOc8nry5Anjxo1j3LhxvHjxgrZt2zJs2DDs7OI/WPtFyAsWnF3ArNOzeB7ynIYFG9LVtSu5Ld9evF7T4kYHN9Po4KYlK/+dvUu/5ScJDAplYN2itC+XP86dRAICAhg1ahRTpkwhLCyMjh07MnjwYHLlMn22lOg8fvWYWadnsfDcQsIkjC8Kf0En505kTafX8NXiRwc30+h3blqyUr2YYTHUcgVtGLbOl6/mvH8x1PfJli0bY8eO5eLFi3zzzTfMmDGDggUL0q9fPx48eBCvclqnteb7kt+zvsl6GhVsxJLzS6i3sh6Tj08mMDgwXnnH144dO2jQoEGsrpkzZw63bt0yOb2I0LNnT+zt7XF2dubYsWNRpgsODqZz584ULlyYokWLsmLFCgDGjRtH8eLFcXZ2pnr16ly9ejVW5dU0Hdy0ZCdbxjTM6lCKnxs5sP/SA+qO382/vnfjnF/u3LmZNm0a586dw9PTk7Fjx2JnZ8dPP/3E06dP41XWHBlyMKzcMFY1WkXF3BX58+Sf1F1Zl3/O/ENQWFC88v6YYhvcNm3axIULF7hw4QJ//fUX3333XZTpRowYQfbs2fHz88PX15fKlSsDUKJECY4cOcLJkyfx9PSkX79+CVIP7RMiIomyubu7i6bFl9+dp1J3/C6x7b9eflh5Ul4EhcY7z1OnTkmTJk0EEBsbGxkzZoy8ePEiAUorcjrgtHTa0kkc5zhKjWU1ZKXfSgkJC0mQvKMyZMgQ+eOPP4z7P/zwg4wfP14qV64szZo1kyJFikjr1q0lPDxcRER++uknKVmypDg4OEinTp0kPDxcli1bJhkyZJDChQuLi4uLSV+Lzp07y8KFC437hQsXllu3br2TLk+ePBIYGPjevI4dOyblypUzscYpH3BEEun3dnLadHDTkr1XIaEyYoOv2PZfL1XHbpdTNx4nSL6HDh2SWrVqCSC5cuWSadOmSVBQUILkvf/Wfmm5rqU4znEUj1Ue8q//v8YAk5CuXLkiJUqUEBGRsLAwKVCggCxfvlysrKzk+vXrEhYWJmXLlpXdu3eLiMiDBw+M17Zt21bWrl0rIiKVK1eWw4cPG895eXmJi4vLO9tvv/0mIiL169c35ikiUq1atTeuFxF59OiR5MmTR7y9vaVEiRLi6ekpd+7ceacO3bp1k59//jmBviLJnw5upm36saSW7KWxMOeHesVY2LEML4LCaDJ1L9N2XCIsPH6dpUqVKsWWLVvYsWMHdnZ2fPfddxQtWpR58+YRFhYWr7zL5izLwvoL+aPKHwiC1w4v2mxsw6Hbh+KV79vy58+PjY0Nx48f53//+x8lSpTAxsaG0qVLkydPHszMzHB1dcXf3x+A7du3U6ZMGZycnNi2bRtnzpyJMt8//vgDHx+fd7YBAwaYXLbQ0FBu3LhBuXLlOHbsGJ9//jl9+vR5I838+fM5cuQIffv2jfPXQPs06eCmpRjl7LOy2asiNYt/xqjN52g94wC3Hr+Md76vZzbZsGEDmTJlol27djg7O7Ny5UrD4484UkpRw7YGKz1WMrzccO69uMc3//uGb7d+y5kHUQeVuOjYsSNz5sxh9uzZfP311wCkSZPGeN7c3JzQ0FBevXpF165dWb58OadOnaJTp068ehV1Zx1vb29cXV3f2UaOHAkY3mNev37dmP7GjRvkzv3mcAgbGxvSp09P06ZNAWjevPkbHU/+/fdfRowYwdq1a98or6aZQgc3LUWxTp+aKa3dGOPpzOmbT6gzfhfrTpjeESI6Sinq1avH0aNHWbp0KWFhYTRr1ozSpUuzZcuWeAU5CzMLmhRqwoamG+hTsg9nHpyh5fqW9NnZB/8n/vEue5MmTdi8eTOHDx+mdu3a0aZ7HciyZs1KYGAgy5cvN57LmDEjz549M+7H1HLz8PDgn3/+QUQ4cOAAmTJlImfOnG/cTylFw4YN2bFjB2BYhPb1gPrjx4/z7bffsnbtWrJnzx7vr4H26dHBTUtxlFI0L5mXjb0qUjC7JT0WHef7JT48exUS77zNzMxo3rw5p0+fZtasWQQEBFCnTh2qVKnCnj174pV3GvM0tHdoz6amm/jW+Vt23dhF4zWN+Wn/T9x9HvfeoKlTp6Zq1ap88cUX710dwdramk6dOuHo6Ejt2rUpVaqU8VyHDh3o0qULrq6uvHwZc2u4Xr16FChQAHt7ezp16sTUqVON51xdXY2fR40axbBhw3B2dmbevHn8/vvvAPTt25fAwECaN2+Oq6srHh4ecai59inTg7i1FC00LJxJ2y4yadsFcmdOx/gWrrjbmr4YakyCgoKYMWMGv/zyC3fv3qVevXr88ssvlChRIt553395nxknZ7DUb6lhsdRirfnGMfaLpYaHh+Pm5sayZcsoVEhPB5bc6UHcpvnkW25VqlQhMYJsUFAQNWrUwNXVlSVLlphUpjlz5tC9e/cEuX/+/Pm5f98wEXG5cuUSJM/42LFjB5kyZTK+uxk+fHiU6a5cuUKZMmWwt7enRYsWBAcbVgsICgqiRYsW2NvbU6ZMGWMHCQtzM7xrFmZZl88BaD59P+O2+hEaFp4g5U6TJg3du3fn0qVLjBw5kv379+Pm5mZs3cVH1nRZGVhmIOsar6OWbS3mnJ5D3RV1+fvU37wIeWFSHr6+vtjb21O9enUd2LRPyicf3BLL8ePHAfDx8aFFixaJWpZ9+/Yl6v1fq1ixovHdzdChQ6NM079/f7y9vbl48SKZM2dm5syZAMycOZPMmTNz8eJFvL296d+//xvXudtmYWPPijQukZuJ/13Ac/p+/O8/T7CyZ8iQgf79+3P58mWGDBnC5s2bcXJywtPTk5MnT8Yr7zwZ8/BrxV9Z1nAZbp+5MeHYBOqvqs+Sc0sICX//o9bixYtz+fJl4+M+TftUfDLB7fnz59SvXx8XFxccHR3faS0BLFq0CCcnJxwdHd/45WhpaYm3tzcODg5Ur16dgIAAAC5dukSdOnVwd3enYsWKnDt37p08Hz58SOPGjXF2dqZs2bKcPHmSe/fu0bZtWw4fPoyrqyuXLl1657p58+bh6uqKo6Mjhw692z28Q4cOb7zwt7S0NH4eM2YMpUqVwtnZmR9//DHGr83ra3fs2EGVKlXw9PSkaNGitGnTxthR4ujRo1SuXBl3d3dq167N7du3Y8w3oYkI27Ztw9PTE4D27duzevVqANasWUP79u0B8PT05L///nunk0fGtKkY94Urk1uX4HJAIPUm7mbpkevx6gzyNmtra4YPH46/vz+DBw9m69atuLi40KRJE+MfNHFVJEsRJlefzD91/yFfxnz8cvAXGq1uxMbLGwmXhGmJalqKkVgD7D72IO7ly5dLx44djfuPHxsG+r4enHrz5k3Jmzev3Lt3T0JCQqRq1aqyatUqEREBZP78+SJimMGhW7duImIYmOrn5yciIgcOHJCqVau+c9/u3bvLsGHDRETkv//+ExcXFxER2b59u9SvXz/KslauXNlY1p07d4qDg4OIiMyePdt47/bt28uyZcuM12TIkEFERLZs2WKcWSIsLEzq168vO3fufOcetra2EhAQ8Ma127dvj3Jwb3BwsHz++edy7949ERFZvHixfPXVV+/kOX/+/CgH9jZr1izKeka2fft2yZIlizg7O0udOnXk9OnT76QJCAiQggULGvevXbtm/No4ODjI9evXjecKFChgrF9Ubj56IS3+3Ce2/ddLl3lH5NHzhBmc/baHDx/Kjz/+KJkyZRJAGjZs+M5g5rgIDw+Xndd3StM1TcVxjqM0W9NMdl3f9UEGgsfV/v37Zfbs2YldjBQHPYjbpM3ClAColKoDTADMgb9FZORb5/8Aqkbspgeyi4h1QgXghODk5ETv3r3p378/DRo0oGLFim+cP3z4MFWqVCFbtmwAtGnThl27dtG4cWPMzMyMjw7btm1L06ZNCQwMZN++fTRv3tyYR1DQu3MF7tmzxzgZbLVq1Xjw4IFJ8xW2atUKgEqVKvH06VMeP35sUj3/97//GQfrAgQGBhoX6TTF68G9gHFwr7W1NadPn6ZmzZoAhIWFvdOtGwxfszZt2ph0n7e5ublx9epVLC0t2bhxI40bN+bChQtxyssUuazTsaBjWWbsvszv/zvP8fGP+f0LF8rbJ+ys/ZkzZ2bYsGF4e3szceJE/vjjD0qVKkW9evUYOnQoZcqUiVO+Sikq5alEhdwV2HRlE5OPT6brf11x/8wdLzcvXLO7Jmg94iJnzpw0aNCAChUqYG9vn9jF0T4xMT6WVEqZA1OAukBxoJVS6o3VHUXEW0RcRcQVmASs/ABljZfChQtz7NgxnJycGDx4cLQdFkyhlCI8PBxra+s3xvicPXs2wcr79jIub+9bWFgQHm54FBUeHm7sWCEiDBw40Fim1zPemyqqwb0igoODgzHPU6dO8b///e+daxcsWBDlwN7XjxEjmzJlivH8rVu3sLKyMj4erVevHiEhIcYOL6/Z2Njw+PFjQkNDgTcHBkceNBwaGsqTJ0+wsbF5b13NzRRdKhdkVdfypE9jTpu/DzJigy9BofGbfSQqmTJlYsiQIfj7+/Prr79y8OBBypYtS+3ateP1ztNMmVG/QH3WNl7LD2V+wP+JP19u+pIe23pw4dGH++PAFLa2tgwcOJCuXbsm6KNfTTOFKe/cSgMXReSyiAQDi4FG70nfCliUEIVLSLdu3SJ9+vS0bduWvn37vrMER+nSpdm5cyf3798nLCyMRYsWGWcoDw8PN77fWrhwIRUqVMDKygo7OzuWLVsGGILKiRMn3rlvxYoVWbBgAWB4p5U1a1asrKxiLO/rd4J79uwhU6ZMZMr0Zvfv/Pnzc/ToUQDWrl1LSIihY0Ht2rWZNWsWgYGGZVVu3rzJvXv3TPsiRaNIkSIEBASwf/9+AEJCQqKclqlNmzZRDuyN/G7wtW7duhnP58qVizt37hh/AR46dIjw8PB3gpNSiqpVqxrzmzt3Lo0aGX4UPTw8mDt3LgDLly+nWrVqJq/z5pg7Ext6VKRt2XzM2H2FxlP2ceHus5gvjAMrKysGDhzIlStXGDVqFMePH6d8+fLUqFGD3bt3xznfVOapaFW0FRubbqRHiR4cuXOEZmubMWjPIG4G3kzAGsROr169uHfvHosWJblfCVpKF9NzS8ATw6PI1/tfApOjSWsL3AbMoznfGTgCHMmXL1+0z5Q/hM2bN4uTk5O4uLhIyZIlje89Ik8Iu3DhQnF0dBQHBwfp16+f8doMGTKIt7e3ODg4SNWqVY3vni5fviy1a9cWZ2dnKVasmPz000/v3PfBgwfSqFEjcXJykjJlysiJEydEJOZ3br169RJXV1dxcHCQgwcPisib79zu3LkjZcqUEWdnZ+nXr5/xvZmIyPjx48XR0VEcHR2lbNmycvHixXfuEd07t8hl6tatm/GdyfHjx6VixYri7OwsxYsXl7/++iumL3msTJo0SYoXLy7Ozs5SpkwZ2bt3r/Fc3bp15ebNmyIicunSJSlVqpQULFhQPD095dWrVyIi8vLlS/H09JSCBQtKqVKl5NKlS3Eqx9Yzd8Rt+P+k8KCNMmfvlQ/+DiswMFDGjh0r2bNnF0CqVKki27dvj/d9H718JGMPjxW3f9ykxD8lZOTBkXL/xf0EKnXsHDhwQHLkyCEPHz5MlPunNOh3biZtMQ7iVkp5AnVEpGPE/pdAGRF5Z8CVUqo/kEdEesQUVJPTIG5LS0tjS0hL+QKeBdF3+Ql2nA+gapFsjPZ0IVvGDzu34YsXL/jrr78YNWoUd+7coWLFivz444+xaoFG5c7zO0w/MZ1VF1eR1jwt7Rza0b54eyxTW8Z8cQLq1q0boaGh/Pnnnx/1vimRHsRtGlMeS94E8kbazxNxLCotSYKPJDUtNrJlTMPsDqX4ycOBfZceUGf8Lv47G/fpr0yRPn16vLy8uHz5MpMmTeLy5cvUqFGDChUqxGvuysiLpZbPXZ7pJ6ZTb2U95vnO+6iLpf7666+sX7+evXv3frR7ap+4mJp2gAVwGbADUgMnAIco0hUF/ImY0iumTa/npiUH5+88lToRi6EOWpUwi6Ga4uXLlzJ16lTJmzevAFKmTBnZsGFDvB9Xngo4JR23dBTHOY5Sc1lNWXVhlYSGfZw6LV26VBwcHBJsTbxPFfqxpElbjC03EQkFugNbgLPAUhE5o5QarpSKPJtpS2BxxBdf01KEwp9lZHW3cnSuVID5B67RYNJuTt988sHvmzZtWr777jsuXLjAn3/+yZ07d6hfvz6lS5dm3bp1cW7JOWZ1ZEatGfxV8y+ypM3CkL1DaLa2Gf9de3fQe0Lz9PQkX758jBs37oPeR9NAT5ysaSbbc+E+vZf58PB5ML1rFaFzxQKYmcX9fVhsBAcHM2/ePEaMGMGVK1coUaIEQ4cOxcPDAzOzuE00JCL8e+1fJh6biP9Tf5yzOuPl7kWpHKVivjiOrly5QqlSpTh06BAFChT4YPdJyfQ7N9Po4KZpsfDoeTA/rDrFptN3+LyADb9/4UIu63Qf7f4hISEsWLCAESNGcPHiRZydnRk6dChNmjSJc5ALDQ9l7aW1TPGZwr0X9yifqzw93XpS3KZ4zBfHwejRo9m2bRubNm2KV2eZT5UObqb5ZOaW1LSEkDlDaqa2cWO0pzMnbjymzvhdrD8Z/8VQTZUqVSo6dOjA2bNnmTdvHkFBQXh6euLi4sLChQuNA9xjw8LMgqaFmrKhiWGx1NMPTtNifQv67uzL1adXE7wO3t7e3Lp1i6VLlyZ43pr2mm65aVoc+d9/jtcSH3yuP6aZWx6GeRQnY9pUH7UMYWFhLF26lF9++QVfX1/y589Pnz59+Oqrr0ifPn2c8nwW/Iw5Z+Ywz3cewWHBNC3UlC4uXciePuFWxN6/fz/NmjXD19cXa2vrBMv3U6BbbqbRwU3T4iEkYjHUycbFUEvgbpv5o5cjPDycDRs2MHLkSPbt20fWrFnp1asX3bp1I3PmuJXn/sv7/HXyL5b5LcNCWdC6WGu+dvw61oulRqdLly4opZg2bVqC5Pep0MHNNDq4aVoCOOL/EK8lPtx+8oruVe3pUc0eC/PEeeq/Z88eRo4cyYYNG7C0tKRz5854e3sbJ8SOrevPrjPVZyobLm/AMrUlXzt+TZtibUhnEb93jY8ePcLBwYEVK1bw+eefxyuvT4kObqbRwU3TEsizVyH8uPYMK4/dpEQ+a8a3cMXWJkOilefUqVOMHj2aRYsWYWZmRtu2benXrx9FixaNU37nH55n4vGJ7Lqxi2zpstHFpQtNCjUhlVncH8UuXryYX3/9laNHj5Iq1cd9pJtc6eBmGh3cNC2BrTtxi0GrThEWLgzzcMDTPU+i9gr09/dn3Lhx/P3337x69YrGjRvTv3//OC+3c/TuUSYcm8Dxe8fJlzEf3Ut0p3b+2pip2LdURYQ6depQvXp1+vXrF6fyfGp0cDONDm6a9gHcevwS7yU+HLzykHpOOfi1iRPW6VMnapkCAgKYNGkSkydP5tGjR1SpUoX+/ftTu3btWAdfEWH3zd2MPzaeC48uUDRLUXq59aJ8rvKxzuvSpUuUKVOGI0eOkD9//lhd+ynSwc00Orhp2gcSFi78tcuwGGpWyzSM+8KFcgm8GGpcBAYGMmPGDMaNG8eNGzdwcXGhf//+NG/eHAsLk9YvNgoLD2OTv2Gx1JuBNyn5WUl6ufWK9WKpv/32G3v27GH9+vV67FsMdHAzjR7npmkfiLmZ4rsqkRZDnXmQ3zae/SCLocaGpaUl3t7eXLp0idmzZxMUFETr1q0pXLgw06ZN4+XLlybnZW5mToMCDVjXeB0DSw/k8pPLfLnpS3pu68nFRxdNzqd3795cvXo1yrX/NC0udMtN0z6Cl8Fh/LLBlwUHr1E8pxUTW7linz1jYhcLMAwjWLduHSNHjuTAgQNkz56dXr160bVr11iPQXsR8oJ5vvOYc2YOL0Jf0LBAQ7q6diWXZa4Yr927dy9ffPEFvr6+7yzOq/0/3XIzjQ5umvYR/et7l34rTvI8KJTB9YvRtqxtknkMJyLs3r2bkSNHsmnTJiwtLenSpQve3t7kyhVzcIrs0atHzDw1k0XnFiEILYq0oJNzJ7KkzfLe6zp37kzq1KmZPHlyfKqSoungZhod3DTtI7v37BV9l51kp18A1YpmZ1Qz5w++GGpsnThxgtGjR7N48WIsLCz48ssv6du3L0WKFIlVPnee32HaiWmsvriatOZp6eDQgXYO7ciQKuohEg8fPsTBwYE1a9ZQunTphKhKiqODm2l0cNO0RCAi/LP/KiM2nsUqrQWjPZ2pVvSzxC7WO65cucLvv//OzJkzCQoKokmTJgwYMIBSpWK3csDlJ5eZfHwyW69uJXOazHR27swXRb4gtfm7PUgXLFjAmDFjOHLkCNu2bcPe3l6vIBCJDm6m0cFN0xKR391n9Fx0nHN3nvFlWVt+qFeMdKnNE7tY77h3755xGMHjx4+pWrUqffv2pXbt2rFajeBUwCkmHJvAwTsHyZkhJ91cu9GgQAPMzf6/ziJCrVq1qFOnDhcvXsTJyYmuXbt+iGolSzq4mUYHN01LZEGhYYzZfJ6/91zBPrslE1q64pAraXaoePbsGX/99Rfjxo3j1q1bFC1aFC8vL7788stYTdS8/9Z+xh8bj+8DXwpmKkhPt55UzVsVgBUrVlCkSBGqVq3Kl19+ScaMGRk+fPiHqlKyo4ObafRQAE1LZGkszBncoDjzvynDs1chNJ6ylz93XiI8POktap8xY0Z69+7NlStXmD9/PunTp6dLly7ky5ePwYMHc/v2bZPy+TzX5yyuv5jfK/9OmITRa3sv2m5qy+E7h9mwYQMeHh7UqVOH//77z+Q8NS0y3XLTtCTk0fNgBq48xeYzdyhX0LAYas5MH28x1Nh63cPyjz/+YM2aNVhYWNCyZUu8vb0pUaKESXmEhoey5uIapp6YalgsNXd5Sj0txagBo/D398fZ2Rn9u+L/6ZabaXRw07QkRkRYduQGw9adIZW5Gb82caK+c87ELlaMLl26xMSJE5k1axaBgYFUrlwZb29vGjRogLl5zO8RX4W+YvG5xcw4NYOnwU+pmasmtxfcxu+UH0ePHv0INUgedHAzjQ5umpZE+d9/Tq8lPpy4/hhP9zwM83DAMk3spsdKDI8fP2bmzJlMnDiRa9euUbBgQby8vOjQoQOWlpYxXv80+ClzTs9h/tn5hISF0LRQU751+TZBF0tNznRwM41+56ZpSVT+rBlY3uVzelazZ+WxG9SbsJujVx8ldrFiZG1tTe/evbl06RJLliwhW7Zs9OjRg7x589KvXz+uX7/+3uutUlvR060nG5tuxLOwJysvrKT+yvpMO6EXNU0ISqnNSqnHSqn1sbxOKaUmKqUuKqVOKqXcokmXWin1l1LKTyl1TinVLOJ4JaXUMaVUqFLKMyHq8j46uGlaEpbK3IzvaxVhybefEy7CF3/uZ/y/foSGhSd20WJkYWHBF198wf79+9m/fz+1atVi3Lhx2NnZ0bJlSw4dOvTe67Omy8qgsoNY23gt1W2rExIW8pFKnuKNAb6Mw3V1gUIRW2cgur82BgH3RKQwUBzYGXH8GtABWBiHe8eafiypacnE01chDFtzhpXHb+KWz5rxLUqQz8b07vdJwbVr15g0aRIzZszgyZMnlCtXDm9vbxo3bhzjigQikmSmKktMpj6WVEoNBx6KyPiI/REYgs4EpVQVoI+INIjFff8EdojIooj980AVEbn9VrrrQFEReR5NPnOA9SLyQWfJ1i03TUsmrNKmYlwLVya2KsGFe4HUnbCL5UdvkFh/oMZFvnz5GDNmDNevX2fChAncuXOH5s2bY29vz7hx43jy5Em01+rAFmuzgHYASikzoCUwP7rESqk/lFI+UWwDIpLkBiI/U74RcSxyHtYRH3+OeAS5TCmVKFPv6OCmacmMh0suNntVwjF3JvosO0H3hcd58iJ5PbLLmDEjPXv2xM/Pj1WrVmFra0vv3r3JmzcvXl5eXL58ObGLmOyJiD/wQClVAqgFHBeRB+9J7y0irlFsI2NxWwsgD7BPRNyA/cDYeFQjznRw07RkKLd1OhZ2Kku/OkXYcuYOdSbsYt+l+4ldrFgzNzencePG7Ny5kyNHjtCoUSOmTJlCoUKFaNq0Kbt3735vyzQwMJA//viD9etj1TfiU/I3hvdcX2FoyUXLhJbbTSBvpEvyRByL7AHwAlgZsb8MiLLjyYemg5umJVPmZoquVexZ1bU86VKZ0+Zvw2KowaFJv7NJVNzd3Zk3bx7+/v7079+fnTt3UqlSJUqXLs3ChQsJCfn/1un9+/f58ccfsbOzY//+/RQtWjQRS56krQLqAKWALe9LaELLbS3QLqLXZFngydvv28Twl8g6oErEoeqAbwLWx2QmdShRStUBJgDmwN9RNVOVUl8AwwABTohI6/flqTuUaFrCeREcyi8bzrLw4DUcclkxoWXSWQw1rl68eME///zD+PHjOX/+PLVr12bz5s1MnDgRb29vsmTJQvHixcmSJQupU6d+Z0uTJo3Jx2OTNnXq1CYNSv9QYjvOTSk1HXgsIgMi9ncDRQFLDC2tb0TkvYEv4joFTMYQLF8AX4nIkYhzPiLiGvHZFpgHWAMBEemuKaVKYQi2mYFXwB0RcTC1HrEVY3BTSpkDfkBNDC8QDwOtRMQ3UppCwFKgmog8UkplF5F778tXBzdNS3hbfe/Sf8VJXgSHMqhe0loMNa7Cw8PZvHkzFhYW1KpVi7t37+Ll5cWGDRuoVasWHh4epEmThuDgYIKDgwkKCjJ+jrwlxPHXx8zMzD5I4DQlbY0aNc4D7YEgIDjSdlNEwiJ/7SI6khwDmovIhY/9vUtMpgS3z4FhIlI7Yn8ggIj8FinNaMBPRP429cY6uGnah3Hv2Sv6LDvJrojFUEd7OpPVMmkthpoQbty4wfDhw1m5ciVLly6lWrVqH+W+IkJYWJgx2IWEhMQpQMY1wG7bti0QOAukfmvrLSLLXpdTKVUcWA+sEpHeH+WLk4SYEtw8gToi0jFi/0ugjIh0j5RmNYbWXXkMjy6HicjmKPLqjGHwH/ny5XO/evVqAlVD07TIwsOFufv9+W3TOazSWjDG04WqRVPm9FWXLl0iY8aMZM+eMuv3Nj39lmkSqkOJBYZR61WAVsCMSOMdjETkLxEpKSIls2XLlkC31jTtbWZmiq/K27GuewWyWqbhqzmHGbrmNK9CwmK+OJkpWLDgJxPYNNOZEtxM6f55A1grIiEicgVDK65QwhRR07S4KpIjI6u7leebCnb8s/8qDSbt4cyt6AdKaymLUqqnUuqsUmqBUqqxUmpoDOmzKKW2KqUuRPybOeJ4g4gZT+Jaji5KqVMRQwv2RDwyjSpdHaXU+Yj5KwdEOj5TKXUiYk7L5UqpGGfgNiW4HQYKKaXslFKpMYxyX/tWmtVEdP1USmUFCgN6FKamJQFpU5kzpEFx5n1TmqcvDYuh/rUraS6GqiW4rkBNEWkD9AOmxpB+APCfiBQC/ovYB9gANFRKxXW+t4Ui4hTRo3I0MO7tBBGdF6dgmMOyONAqUhD0FhEXEXHGMEdl97evf1uMwU1EQiMy2oLhJeZSETmjlBqulPKISLYFw0h4X2A70Pd9I+E1Tfv4KhbKxhavSlQrmp1fN56j7cyD3H7yMrGLpX0gEUMACgCblFL9gSARua+UyqiUuqKUShWRzirSfiNgbkQWc4HGYBy/tgMweS7KyETkaaTdDBiGjL2tNHBRRC6LSDCwOKI8xusjhiOki+b6N5j0zk1ENopIYREpKCIjIo4NFZG1EZ9FRL4XkeIR0XmxKflqmvZxZc6Qmult3RnVzInj1x5TZ/xuNp66HfOFWrIjIl2AW0BV4B6GIQGIyDMMgap+RNKWwEoRCQE+izQw+w4QeV7II0BFAKVUm2hmM4l2MmSlVDel1CUMLbeeUSR579yVSqnZEWUqCkyKqf56hhJN+8QopWhRKh8be1Ukv016ui44Rt9lJwgMCk3somkfTk4MA6pf+xvDlFxE/Dv77QsiWmuRW0j3gFwR5xZEM5tJtOu0icgUESkI9AcGx7YCIvJVxP3PAi1iSq+Dm6Z9ouyyZmD5d+XoUc2eFcduUH/ibo5dS/qLoWpx8hJI+3pHRPYC+SOWvjEXkdMRp+4qpXICRPwbeTKOtBH5xNhyU0qNeH0sirIsJuJx51ti7LwYMUh9MdAspgrr4KZpn7BU5mb0jlgMNTRMaD59PxP+vZAsFkPVYuUsYP/WsX8wLBwaudW2FsPsJ0T8uybSucLAaYi55SYig14fA+MsVq/VB6KaLSXKzosRc1naR+SjAA/gXEwV1sFN0zRK5c/CJq+KNHTOyR//+tHirwNce/AisYulJZxdQAn15lxsCzDM87go0rGRQE2l1AWgRsT+a1Ux9JqMi+5KqTMRLbnviQigSqlcSqmNEH3nRUABc5VSp4BTGB6xxjgsQa/ErWnaG9b43GTw6tOIwE8eDjR1y53s56dMSeI6Q4lSagKwTkT+jdj3BBqJyJcmXPsZhu781WNd4ESiW26apr2hkWtuNvWqSPFcVvRedoLui5LfYqhalH4F0gMopSZhaJX9bOK1+YBkNT+lbrlpmhalsHBh+s5L/LHVj2wZ0zDuC1c+L2iT2MX65Om5JU2jW26apkXJ3EzRrao9K7uWI10qc1r/fYDfNiXfxVC1T4sObpqmvZdzHmvW96xAy1L5+HPnZZpO28vFe4GJXSxNey8d3DRNi1H61Bb81tSJv7505+ajlzSYtJv5B66SWK81NC0mOrhpmmayWg452OJVidJ2NgxefZpO/xzhfmBQYhdL096hg5umpVSH/4b5zeDEYgh6lmDZZrdKy5wOpfixYXF2XbhPnfG72H7+XswXatpHpIObpqVUZhYQ4AervoUxhWD513B+E4QGxz/riMVQ13Yvb1gMdfZhfkyhi6FqyZMeCqBpKVl4ONw4BCeXwplV8PIhpMsMDk3A6QvIWwbM4vc37quQMEZvPs+svVcolN2S8S1dcciVKYEqoL1NDwUwjQ5umvapCA2GS9vg1DI4twFCX0KmfODUzBDoPotycWST7fILoM+yEzx+EULf2kX4poIdZmZ6ZpOEpoObaXRw07RPUVCgIcCdWmYIeBIG2R3AuTk4eoJ13pjziMLD58EMWHGS//nepby9Db83dyVHprQxX6iZTAc30+jgpmmfusAAwyPLU8sMjzAB8pUzBLrijSF9llhlJyIsOXydn9b5kiaVGb81caKuU86EL/cnSgc30+jgpmna/3t4BU4th1NL4b4fmKUC+xqGQFe4LqROb3JWlwMC8Vriw8kbT/iiZB5+bOhAhjQWH7DwnwYd3Eyjg5umae8SgTsnDR1RTq+AZ7chtSUUbWAIdHZVwDzmQBUSFs6Efy8wZcdF8mVJz/gWrpTIl/mDFz8l08HNNDq4aZr2fuFhcHWvIdD5roWgJ5AhGzg0BecvILc7xLAkzqErD/Fe4sOdp6/oVb0QXasUxMJcj0SKCx3cTKODm6Zppgt5BRe3GgKd3xYIC4IsBcCpuWHLWijaS5++CmHo6tOs9rmFu21mxrdwJW8W0x9zagY6uJlGBzdN0+Lm1RM4u84Q6K7sAgRyuhpacw5NwSrqTiRrfG4yeNVpBBjeyIEmJfRiqLGhg5tpdHDTNC3+nt6GMysNge62DygzyF/REOiKNYS0bw7qvv7wBb2XnuCQ/0MauuTil8aOZEqXKnHKnszo4GYaHdw0TUtYAX5werkh0D26AuZpoHBtQ6ArVAss0gBvLoaaPWMaxrVwpWwBvRhqTHRwM40ObpqmfRgicPNoxNRfK+F5AKTJBMU9DIHOtgKYmXHi+mO8lvjg/+A5XSoXxLtGYVJb6M4m0dHBzTQ6uGma9uGFhcKVHXByGZxbD8GBkDEXOBp6XD7PXJxfNp5l0aHrOOa2YkLLEhTMZpnYpU6SdHAzjQ5umqZ9XMEvwG+TIdBd3ArhoZC1CDg1Z1faKvTa8oiXIWEMrl+cNmXy6c4mb9HBzTQmBTelVB1gAmAO/C0iI9863wEYA9yMODRZRP5+X546uGmaxouH4LvaEOiu7QMgOGdJFr0sw4Q7TrgVs2dUM2dsLNMkbjmTEB3cTBNjcFNKmQN+QE3gBnAYaCUivpHSdABKikh3U2+sg5umaZEtnzuNvK/OUzr0ICrAl3Blzu4wJ7ZaVKZ206+o6GiX2EVMEnRwM40pb21LAxdF5LKIBAOLgUYftlialvTs2LGDBg0axOqaOXPmcOvWLeN+lSpVKFKkCK6urri6urJ8+XIAzM3NjcdcXV0ZOXJkdFm+QUTo2bMn9vb2ODs7c+zYsSjTvX3fe/cMK2cHBQXRokUL7O3tKVOmDP7+/rGqX0K69SSEDuO3YDvqJl63PNiWtgGfZ7zHL+ETcF9WhtMTPeH8ZggLSbQyasmHKbOY5gauR9q/AZSJIl0zpVQlDK08bxG5/nYCpVRnoDNAvnz5Yl9aTUtm5syZg6OjI7ly5TIeW7BgASVLvvmHd7p06fDx8Yl1/ps2beLChQtcuHCBgwcP8t1333Hw4MEo00Z135kzZ5I5c2YuXrzI4sWL6d+/P0uWLIl1OeIqODiYy5cv4+fnR3BwMBUqVOD48eNMm7uUCcHBLFm8mEYlc+O76W8cH2+DRVshXRbDYqvOX0Ce0vFebFVLmRJqiu51wCIRCVJKfQvMBaq9nUhE/gL+AsNjyQS6t6YluKFDh5IlSxa8vLwAGDRoENmzZycwMBBPT09Onz6Nu7s78+fPRynF8OHDWbduHS9fvqRcuXL8+eefrFixgiNHjtCmTRvSpUvH/v37E7yca9asoV27diilKFu2LI8fP+b27dvkzGnaEjNr1qxh2LBhAHh6etK9e3dEJEE7cYSHh3Pz5k38/PyM2/nz5/Hz8+PGjRvkzZuXwoULU7hwYRwdHXnw4AG3bt1i7ty51KxZE4CS3StELLb6n2FpHp+FcGRmxGKrnoZAl71YgpVZSwFE5L0b8DmwJdL+QGDge9KbA09iytfd3V00Lam6cuWKlChRQkREwsLCpECBArJ8+XKxsrKS69evS1hYmJQtW1Z2794tIiIPHjwwXtu2bVtZu3atiIhUrlxZDh8+bDxXuXJlKVy4sLi4uIiLi4vcv39fRETMzMyMx1xcXGTx4sUiIuLl5fXG8dfbb7/9JiIi9evXN5ZBRKRatWpv3C/yfR0dHcXFxUWGDx8u4eHhIiLi4OAg169fN6YrUKCABAQExP8LGOHRo0eSJUsWyZkzp1SuXFk6deokY8aMkTVr1si5c+ckKCjImPbWrVtSvHhx+eKLL974ekbp1VMRn0Ui85qKDMss8qOVyNRyIrv/EHl8/f3XJnPAEYnh96vexKSW22GgkFLKDkNvyJZA68gJlFI5ReR2xK4HcDYB4q6mJZr8+fNjY2PD8ePHuXv3LiVKlMDGxobSpUuTJ08eAFxdXfH396dChQps376d0aNH8+LFCx4+fIiDgwMNGzaMMu/YPJb8448/EqQ+CxYsIHfu3Dx79oxmzZoxb9482rVrlyB5v4+1tTW3b98mderUMaZNlSoVo0ePpl69ejG3HNNkBJeWhi3wnmGx1ZNL4d8f4d9hYFve0KIr3ijWi61qKUOMwU1EQpVS3YEtGFpls0TkjFJqOIa/INYCPZVSHkAo8BDo8AHLrGkfRceOHZkzZw537tzh66+/BiBNmv/vkm5ubk5oaCivXr2ia9euHDlyhLx58zJs2DBevXqVIGXw9vZm+/bt7xxv2bIlAwYMIHfu3Fy//v+vt2/cuEHu3LnfSf/6WMaMGWndujWHDh2iXbt2xuvz5MlDaGgoT548wcYmYafAMiWwAWTNmpX69evH/gaW2aHMt4bt4WXDYqsnl8J6L9jY1zDll5MnFKkLqdLFPn8tWTLpnZuIbAQ2vnVsaKTPAzE8rtS0FKNJkyYMHTqUkJAQFi5cyO7du6NM9zqQZc2alcDAQJYvX46npydgCCbPnj2Lcxliarl5eHgwefJkWrZsycGDB8mUKdM779tCQ0N5/PgxWbNmJSQkhPXr11OjRg3j9XPnzuXzzz9n+fLlVKtWLXkPms5SACr3g0p94fYJw/u5U8vh/AZIndEwibOTJ9hVNmmxVS350t9dTYtG6tSpqVq1KtbW1pibm0ebztramk6dOuHo6EiOHDkoVaqU8VyHDh3o0qVLjB1KXr58iaurq3G/Tp06Jg0HqFevHhs3bsTe3p706dMze/Zs4zlXV1d8fHwICgqidu3ahISEEBYWRo0aNejUqRMA33zzDV9++SX29vZkyZKFxYsXx3jPZEEpyOVq2GoOB/89cCpisdUTCyFDdnBsZliDLrdbjIutasmPnn5L06IRHh6Om5sby5Yto1Ch6Bfh1JKRkFdw4X+GQOe3BcKCIUvBSIut2id2CWOkB3GbRg8Q0bQo+Pr6Ym9vT/Xq1XVgS0lSpTWsStBiPvS5AB6TIVNu2DkKJrvDX1Vg/1R4djexS6rFk265aZqmPb0Fp1cY3tHdPmFYbNWuEji9XmzVKrFLaKRbbqbRwU3TNC2ygPOGTiinlsIjf8Niq0XqGAJdoZrGxVYTiw5uptGPJTUtGhMnTqRYsWK0adOG1atXM3z48Pemr1KlCrH5g83f35+FCxca93fs2EGmTJmM8z++7tE4bNgwcufO/cbck48fP+bUqVN06NAhTnUDw9g3Z2dnnJycKFeuHCdOnIgyXYcOHbCzszPe+/V4PDFxXstkJ1sRqDYIevrAN/+Ce3vw3wtL2sDYQrC2B1zZDeHhiV1S7T10b0lNi8bUqVP5999/yZMnD+XKlWPt2rUJmv/r4Na69f/PiVCxYkXWr1//Tlpvb2/69OnzxjFra2tu3LjBtWvX4jRXq52dHTt37iRz5sxs2rSJzp07Rzsv5ZgxY4zDG16LzbyWyZJSkLeUYav9G1zeYWjNnVoBx/4xLLbq1MzQosvhpHtcJjG65aZpUejSpQuXL1+mbt26jBo1ijRp0pA1a1aePXuGnZ0dISGGmemfPn36xv6yZcsoXbo0hQsXNo6L8/f3p2LFiri5ueHm5sa+fYZ1ywYMGMDu3btxdXWN80wkDRs2jHP3/XLlypE5c2YAypYty40bN2J1fXTzWqZI5hZQqAY0/Qv6XoRmMyGnMxyYBn9WhKllYdcYw2NMLUnQwU3TojB9+nRy5crF9u3byZ49O25uboBhUHaVKlXYsGEDAIsXL6Zp06akSpUKMAyYPnToEOPHj+enn34CIHv27GzdupVjx46xZMkSevbsCcDIkSOpWLEiPj4+eHt7AxiDnaurKyNGjDCW548//jAer1q1qvF4yZIljUF0+/btbzy6fL2VK1cuxvrOnDmTunXrRnt+0KBBODs74+3tTVBQEAA3b94kb968xjR58uTh5s2b0WWRcqRObxgI3nqJocdl/XGQLjNs+wUmuMDMWnBoBjy/n9gl/aTpx5KaFoPbt2+TLVs2437Hjh0ZPXo0jRs3Zvbs2cyYMcN4rmnTpgC4u7sb10YLCQmhe/fu+Pj4YG5ujp+fX7T3is1jSTAEztfrxVWtWjVOy+Zs376dmTNnsmfPnijP//bbb+TIkYPg4GA6d+7MqFGjGDp0aJRpPznps0Cpbwzbo6v/3+NyYx/YPAAKVjM8tixaD1JnSOzSflJ0cNO0GKRLl44nT54Y98uXL4+/vz87duwgLCwMR0dH47nXc0++nncSDK2uzz77jBMnThAeHk7atGkTrGyvXr0iXTrDfInbt283tgAjS58+Pfv27WPKlCnGQLxx40Zy5crFyZMn6dixI5s2bYp2TsnX03mlSZOGr776irFjxwKYPK/lJyOzLVT83rDdPWOY3/LUcrjQEVKlh6L1DYGuYFUwT5XYpU3xdHDTtBgUK1aM+fPnv3GsXbt2tG7dmiFDhsR4/ZMnT8iTJw9mZmbMnTuXsLAwIP7zTgL4+fkZg2tMLbdu3brRrVs34/61a9do2rQp8+bNo3DhwtFe93p9OBFh9erVxvuZMq/lJ+szB6j5E1T/Ea7tN7Tmzqwy/JveBuqONjza1D4Y/c5N02JQqVIljh8/TuQxoW3atOHRo0e0atUqxuu7du3K3LlzcXFx4dy5c2TIYHg85ezsjLm5OS4uLjF2KIn8zu31UjtgaK3FaSZ9YPjw4Tx48ICuXbvi6ur6xjI89erVMz7ubNOmDU5OTjg5OXH//n0GDx5sTFOgQAHs7e3p1KkTU6dOjVM5UjQzM8hfHhqON7yfa7nIMGmz1Sfcwv1I9CBuTTNBr169aNiwoXHs2fLly1mzZg3z5s1LtDIFBQVRuXJl9uzZg4WFfgjzqdCDuE2j/0domgl++OEH4xiuHj16sGnTJjZu3BjDVR/WtWvXGDlypA5smhYF3XLTNE1LRnTLzTT6nZumaZqW4ujgpmmapqU4OrhpmqZpKY4ObpqmaVqKo4ObpmmaluLo4KZpmqalOIk2FEApFQBcjePlWYGUMuW2rkvSk1LqAbouSVV86mIrItliTvZpS7TgFh9KqSMpZZyHrkvSk1LqAbouSVVKqktSpR9LapqmaSmODm6apmlaipNcg9tfiV2ABKTrkvSklHqArktSlZLqkiQly3dumqZpmvY+ybXlpmmapmnR0sFN0zRNS3GSXHBTStVRSp1XSl1USg2I4nwapdSSiPMHlVL5I50bGHH8vFKq9kcteBTiWhelVE2l1FGl1KmIf6t99MK/Wc44f08izudTSgUqpfp8tEJHI54/X85Kqf1KqTMR35u0H7Xwb4nHz1cqpdTciDqcVUoN/OiFf7OcMdWjklLqmFIqVCnl+da59kqpCxFb+49X6qjFtS5KKddIP1snlVItPm7JUyARSTIbYA5cAgoAqYETQPG30nQFpkd8bgksifhcPCJ9GsAuIh/zZFqXEkCuiM+OwM3kWI9I55cDy4A+yfjnywI4CbhE7Nsk45+v1sDiiM/pAX8gfxKuR37AGfgH8Ix0PAtwOeLfzBGfMyfx70l0dSkMFIr4nAu4DVgnVl1SwpbUWm6lgYsicllEgoHFQKO30jQC5kZ8Xg5UV0qpiOOLRSRIRK4AFyPySyxxrouIHBeRWxHHzwDplFJpPkqp3xWf7wlKqcbAFQz1SGzxqUst4KSInAAQkQciEvaRyh2V+NRFgAxKKQsgHRAMPP04xX5HjPUQEX8ROQmEv3VtbWCriDwUkUfAVqDOxyh0NOJcFxHxE5ELEZ9vAfcAPQtJPCS14JYbuB5p/0bEsSjTiEgo8ATDX9GmXPsxxacukTUDjolI0AcqZ0ziXA+llCXQH/jpI5TTFPH5nhQGRCm1JeKxUr+PUN73iU9dlgPPMbQOrgFjReThhy5wNOLz/zY5/p+PkVKqNIaW36UEKtcnySKxC6BFTynlAIzC0GpIjoYBf4hIYERDLjmzACoApYAXwH9KqaMi8l/iFitOSgNhGB5/ZQZ2K6X+FZHLiVssTSmVE5gHtBeRt1uqWiwktZbbTSBvpP08EceiTBPxWCUT8MDEaz+m+NQFpVQeYBXQTkQS8y+4+NSjDDBaKeUPeAE/KKW6f+Dyvk986nID2CUi90XkBbARcPvgJY5efOrSGtgsIiEicg/YCyTWPIfx+X+bHP/PR0spZQVsAAaJyIEELtsnJ6kFt8NAIaWUnVIqNYaX4GvfSrMWeN0ryhPYJiIScbxlRA8xO6AQcOgjlTsqca6LUsoaww/5ABHZ+7EKHI0410NEKopIfhHJD4wHfhWRyR+p3FGJz8/XFsBJKZU+IlBUBnw/UrmjEp+6XAOqASilMgBlgXMfpdTvMqUe0dkC1FJKZVZKZcbwhGPLByqnKeJcl4j0q4B/RGT5ByzjpyOxe7S8vQH1AD8Mz5sHRRwbDnhEfE6LoefdRQzBq0CkawdFXHceqJtc6wIMxvBOxCfSlj251eOtPIaRyL0lE+Dnqy2GjjGngdHJtS6AZcTxMxgCdN8kXo9SGFrOzzG0PM9EuvbriPpdBL5KBt+TKOsS8bMV8tb/edfErk9y3vT0W5qmaVqKk9QeS2qapmlavOngpmmapqU4OrhpmqZpKY4ObpqmaVqKo4ObpmmaluLo4KaleEopG6WUT8R2Ryl1M+JzoFJqamKXT9O0hKeHAmifFKXUMCBQRMYmdlk0TftwdMtN+2QppaoopdZHfB4WscbZbqXUVaVUU6XU6Ig1zzYrpVJFpHNXSu1UhnX2tkTMBahpWhKjg5um/b+CGKal8gDmA9tFxAl4CdSPCHCTMKzD5Q7MAkYkVmE1TYueXhVA0/7fJhEJUUqdwrDw5OaI46cwLDJZBMPisVsjVjkwx7BsjKZpSYwObpr2/4IARCRcKRUi//9COhzD/xWFYS7AzxOrgJqmmUY/ltQ0050HsimlPgdQSqWKWHNP07QkRgc3TTORiARjWDpmlFLqBIaZ28slaqE0TYuSHgqgaZqmpTi65aZpmqalODq4aZqmaSmODm6apmlaiqODm6Zpmpbi6OCmaZqmpTg6uGmapmkpjg5umqZpWorzfyxI/UogoFr7AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Illustrate one step of the explicit trapezoidal method\n",
    "\n",
    "def degradation_model1D( t, y, params ):\n",
    "    dydt = -params[0] * y;\n",
    "    return dydt\n",
    "fn = degradation_model1D\n",
    "\n",
    "# set up the usual problem, with a faster degradation rate \n",
    "# to exaggerate the differences somewhat.\n",
    "y0 = 1\n",
    "params = (5,)\n",
    "h = 0.1 # step size\n",
    "t0 = 0 # first time step\n",
    "\n",
    "# one step of trapezoidal method\n",
    "yprime = fn(t0,y0,params);\n",
    "yhatFE = y0 + h*yprime;\n",
    "yprimehatFE = fn(t0+h,yhatFE,params);\n",
    "yhat = y0 + h*(yprime+yprimehatFE)/2;\n",
    "\n",
    "# To get the curvature, compute the true solution over 10\n",
    "# steps within this single step.\n",
    "tru_t = np.linspace(0,h,10)\n",
    "tru_y = y0*np.exp(-params[0]*tru_t)\n",
    "\n",
    "# Plot the true solution, along with the interdediate estimate\n",
    "# and the slopes.\n",
    "plt.plot( tru_t, tru_y, 'k' );\n",
    "plt.plot( [t0, t0+h], [y0, yhatFE]);\n",
    "plt.plot( [t0+h-0.02,t0+h+0.02], [yhatFE-0.02*yprimehatFE,yhatFE+0.02*yprimehatFE]);\n",
    "plt.plot( [t0,t0+h],[y0,yhat]);\n",
    "plt.annotate( \"y0=%0.2f\\nf(y0)=%0.2f\"% (y0,yprime), (t0,y0), (t0+0.03,y0-0.05), \n",
    "                  xycoords='data', arrowprops={'arrowstyle':'->'});\n",
    "plt.annotate( \"y1=%0.2f\\nf(y0)=%0.2f\"% (tru_y[-1],fn(t0+h,tru_y[-1],params)), (t0+h,tru_y[-1]), (t0+h+0.03,y1-0.05), \n",
    "                  xycoords='data', arrowprops={'arrowstyle':'->'});\n",
    "plt.annotate( \"yhatFE=%0.2f\\nf(yhatFE)=%0.2f\"% (yhatFE,yprimehatFE), (t0+h,yhatFE), (t0+h-0.04,yhatFE-0.05), \n",
    "                  xycoords='data', arrowprops={'arrowstyle':'->'});\n",
    "plt.annotate( \"yhat=%0.2f\"% (yhat), (t0+h,yhat), (t0+h-0.02,yhat+0.1), \n",
    "                  xycoords='data', arrowprops={'arrowstyle':'->'});\n",
    "plt.annotate( \"slope of blue line = %0.2f\" % yprime, (0.005,0.7))\n",
    "plt.annotate( \"slope of green line = %0.2f\" % ((yprime+yprimehatFE)/2), (0.05,0.85) );\n",
    "plt.xlabel( \"Time\");"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}