{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Runge-Kutta 4th Order Method\n",
    "\n",
    "Our goal is to implement the fourth order Runge-Kutta method commonly refered to as RK4. We will compare its accuracy to those of the forward Euler and explicit trapezoidal methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.integrate\n",
    "import numpy.linalg\n",
    "import matplotlib.pyplot as plt\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": 94,
   "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)\n",
    "\n",
    "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",
    "    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",
    "        yhatFE = y[:,i-1] + h*yprime;\n",
    "        yprimehatFE = fn(t[i],yhatFE,params);\n",
    "        y[:,i] = y[:,i-1] + h/2*(yprime+yprimehatFE);\n",
    "\n",
    "    return (t,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implement the method\n",
    "\n",
    "\\begin{eqnarray*}\n",
    "f_1 &=& f\\left(t,\\widehat{y}(t)\\right)\\\\\n",
    "f_2 &=& f\\left(t+\\frac{\\Delta t}{2}, \\widehat{y}(t) + \\frac{\\Delta t}{2} f_1\\right)\\\\\n",
    "f_3 &=& f\\left(t+\\frac{\\Delta t}{2}, \\widehat{y}(t) + \\frac{\\Delta t}{2} f_2\\right)\\\\\n",
    "f_4 &=& f\\left(t+\\Delta t, \\widehat{y}(t) + \\Delta t f_3 \\right)\\\\\n",
    "\\widehat{y}(t + \\Delta t)&=& \\widehat{y}(t) + \\frac{\\Delta t}{6} \\cdot (f_1 + 2f_2 + 2f_3 + f_4)\\\\\n",
    "\\end{eqnarray*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RK4( fn, t, y0, params ):\n",
    "    \"\"\"\n",
    "    [t,y] = RK4( 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",
    "    Uses a 4th order explicit Runge-Kutta method\n",
    "    \"\"\"\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run and plot result from RK4.\n",
    "\n",
    "Run the RK4 method on the degradation model with degradation rates (2,5) and initial conditions (1,1) with time beginning at 0 and ending at 2.\n",
    "\n",
    "Then plot the true solution (solid lines) and the forward Euler's solution (dashed lines), explicit Trapezoidal's solution (dotted lines), and RK4's solution (dash-dotted lines).\n",
    "\n",
    "Then plot the difference over time between the true solution and each of the other methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Write code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the error vs timestep size\n",
    "\n",
    "For 10 time step sizes ranging from $10^{-4}$ to 0.2, determine the error for forward Euler, explicit trapezoidal, and RK4. To compute the error, use the mean over time of the 2-norm at each time step. \n",
    "\n",
    "Plot the simulations in one figure and the error vs step size in another."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bonus Code\n",
    "### Examine substeps\n",
    "\n",
    "Examine the substeps taken to get the final slope estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnkAAAEGCAYAAAAZlYMoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAB570lEQVR4nO3dd1zU9R8H8Nfn7hjH3nucjOM4BEQRRXFructB7lVq7nKUlTO1tNJKM0c50dxajtziCDeKLAcCCoLsPY8b398fiD9UXMjm/Xw8eMTd9/P9fN/fbwjv+0zGcRwIIYQQQkjDwqvtAAghhBBCSNWjJI8QQgghpAGiJI8QQgghpAGiJI8QQgghpAGiJI8QQgghpAES1HYAb8vExIQTiUS1HQYhhNQrN27cSOc4zvQd6zATCAQbADQFNRIQUttUACIUCsXYFi1apFZUoN4leSKRCMHBwbUdBiGE1CuMsbh3rUMgEGywsLBwNTU1zeLxeLT+FiG1SKVSsbS0NGlycvIGAH0rKkOfxAghhLyppqamprmU4BFS+3g8HmdqapqD0pb1isvUYDyEEELqNx4leITUHU/+Pb40l6MkjxBCCCGkAap3Y/IIIYQ0TsnJyfyOHTu6AEB6eroaj8fjjIyMFABw69atO5qamu/Uyjh16lRrhUKBtWvXJgJAVFSUeufOnV1u3bp128TERPnud0BIzaIkjxBCSL1gYWGhvHv37m0AmDFjhpWOjo5y0aJFKWXH5XI51NTUKl3/0qVLHzdt2tTt5s2bGc2bNy+eNGmS7Zw5cxIpwSP1VbUleYyxTQB6A0jlOO6FQYGMMQZgJYCeAAoBjOY47mZ1xUMIIaThGTBggEhDQ0MVERGh5ePjk6+np6cqn/w5Ozu7HTly5L6Li0vJmjVrjNauXWsul8tZ8+bNCwICAuIEgv//GdTR0eF++OGHR5MmTbKbPn16SkFBAX/ixImZtXZzhLyj6mzJ2wJgNYCAlxzvAcD5yVcrAGuf/JcQQkgd98W+UNuo5DytqqxTbKFb+NNAz0dve15SUpL6zZs37woEAsyYMcOqojI3b97U3Ldvn1FwcPBdDQ0Nbvjw4Xbr1q0znjJlSkb5coMGDcrZtGmTyYQJE0Tnzp27W9l7IaQuqLYkj+O4C4wx0SuKfAAggOM4DsAVxpgBY8yS47ik6ojnQXoBDuzYgUEDusHG+lVhEUIIqU/69++fVb5FriLHjx/XjYiI0PL09HQFgOLiYp6ZmZmiorKTJ09OLS4uZp6enrJqCJeQGlObY/KsAZT/xJbw5L0XkjzG2HgA4wHAzs6uUhc7fSMK+je1cCT0FoodduO98X0hNXWtVF2EENLYVabFrbro6Oioyr4XCAScSvX0JWQyGQMAjuOYv79/xu+//574uvp4PB54PFp8gtR/9eKnmOO4PziO8+Y4ztvUtHK78ox73wsWdmHQKk6G2qOWODMnFLPWzcPVx1erOFpCCCG1RSQSyW7duqUNAEFBQVqJiYkaANC9e/fcI0eOGCYmJgoAICUlhR8VFaVem7ESUt1qM8lLBGBb7rXNk/eqTe9FP6DVKB14Z/wK3UI5mtzqhEfnS3f6KZAXIEeWU52XJ4QQUs1GjhyZlZWVxXdycnJbuXKlmb29fTEAtGjRonju3LmJXbp0EYvFYmnnzp3Fjx49qvxUXELqAVY6JK6aKi8dk3fkJbNrewGYgtLZta0ArOI4zud1dXp7e3Pvundt6r0rUK4ahojHbRE1+StM6uKKDYd3YV/yTvi5t8Qw12FwMHB4p2sQQkhdwhi7wXGc97vUERoa+tDT0zO9qmIihLy70NBQE09PT1FFx6pzCZWdADoCMGGMJQBYAEANADiOWwfgKEoTvGiULqEyprpieZ6ZS2sULwzE6WPR2H36ARIfJMHpihr6qn2G+49uYFj4aLiLJBguHY72Nu1rKixCCCGEkCpTnbNrh7zmOAdgcnVd/3U0Te2xbIQd7M/cheeRj6EWpYYU3bbgoTOcM5ojKu0adsr2PE3y5Eo51PjUsk8IIYSQ+qFR73jBGMOkTk5ISJTCWm8fLt1Twfr6RTx06AOoWqF9LxEAILkgGQMPD0Qfhz4YKhkKWz3bV1dMCCGEEFLL6sXs2mrFV4PNyA1Ibz0bbV3DoWxfBJHsEroKTsG9mSMA4PbpVLxX9BF239mNXn/3wtQzU3E16SqqczwjIYQQQsi7oCQPABiDac9vkN97PZrpP0ChRzIOdf4IShWHophYpF9MhPHVppj+4DeM05uOsLQwfHrqU6QVpdV25IQQQgghFWrU3bXP0/EeDIWxDcKupOL3a0mIzFFgcdoFND20ETkt+yFa1h28E/aY5rgCJu/JYaZlBgD48vyXsNKxwmDJYFhoW9TyXRBCCCGEUEveCwRN/DBpSH98388drWN/wwalDGpzvoVxVCBa/DsNLWxSkZdRDA9zDwBAiaIECk6BzZGb0X1/d8w6PwshqSHUlUsIIdVg9uzZFk5OTm5isVgqkUikgYGB2gCwaNEis7y8vCr5mxYSEqLZrFkzibq6evP58+ebV0WdhNQGSvJeYmhzUwwzf4QFsp/wT8x/SF61BXpdOkF/+7fo6XAHhhbaAIDATVHo+/hTHOh6CCOkI3Dp8SWMPDYS/0T/U7s3QAghDczp06e1T5w4YRAeHn47Kirq9tmzZ6McHBxKAGD9+vXm+fn5VfI3zczMTLFy5cr4Tz/9NKUq6iOktlCS9zJqQuh+egzFjj3xBQIQd2QuTnw4ETZr1sBk+DAAQPGDh9DS5uFO0GOcXvYAfo/74Wiv45jXeh4623UGAJyJO4O1t9YivYjWDyWEkHeRmJioZmRkpBAKhRwAWFpaKkQikXzJkiVmqampah06dBC3atVKDAAHDhzQa9asmUQqlbr26NHDIScnhwcA1tbW7hMmTLARi8VSd3d314iICI3nr2Ntba3o0KFDoZqaGnXJkHqNxuS9ipoQWsO2o+TEfIy4+hv+PZ2Dhe7L8Z2aBjQUCiROnAhrAOIvv0VorBaCjz5ExIVE9Pj0Pehr6AMAQlJDsPX2VvwZ/ie6i7pjmHQY3Izdave+CCGkKvzRyeWF91z7ZKLdjDTI8nnY2sf5heOeg9LRakIG8pIF2DnE8Zlj48/ee9XlPvzww9ylS5daiUSipn5+frlDhgzJ7NWrV/7cuXNT165da37+/PkoS0tLRVJSkuD777+3vHDhQpSenp5qzpw5FosXLzZfvnx5EgDo6+sroqKibq9evdp46tSptmfPno1+p+dASB1FLXmvw+NBvccSqHr+DM6tP/bfTMCQP68grVABi/nzwMnlyJo0Cp6P/0b/aW6wcNCHkWVpV25hbglmtpiJwx8exkDxQJyJP4PBRwZjwaUFtXxThBBS/+jr66siIiJur169Os7U1FQxatQox1WrVhk/X+7cuXPaMTExmj4+PhKJRCLdtWuXcXx8vHrZ8VGjRmUCwLhx4zJDQkJ0avIeCKlJ1JL3hng+n6C3D8B3S8KpvWsxfdU9zB7dH00PHUTqypXI2rYdgrNn0S1gK9R11MCpOBz+7RYYY/Dt74hvWn2DqV5T8ff9v2GpYwkAyC/Jx56oPejv1B8Gmga1e4OEEPK2XtXypqGjeuVxXQvF61ruKiIQCNC7d++83r1753l4eBRt27bNeNq0aRnly3AcBz8/v9zDhw8/qKgOHu//7RuMMeqSJQ0WteS9pR4SQyzTP4A/FN/gl/Xrcfh+Niy++Qb2O/6Ctk9LqFlZAQA4lQrNutiiKL8Eh369hUOrbqE4BRjpNhLd7LsBAC4+vohfbvyCbvu6YeGlhbifdb82b40QQuq00NBQjfDw8Kdj6EJCQoQ2NjYlAKCtra0sG3fXsWPHguDgYJ2y8Xa5ubm8sLCwp+cFBAQYAcDGjRsNvby8Cmr2LgipOdSS97bUNKE+9gR42wfiz7Qf8M3eNNxNHotZ7zWDlZcXAECRmYm4kSNhOXkyhi7shsgLjxF87CH2fH8dfaZ6wk5a2rvwvuh9NNFvgh13duBI7BHsv78frSxbYXXn1dAUaNbmXRJCSJ2Tm5vLnzZtml1ubi6fz+dzIpFItnXr1jgAGDVqVHr37t3F5ubmJVevXo1av379w8GDBzuUlJQwAFiwYEGih4eHDACysrL4YrFYqq6uzu3atSv2+evEx8cLWrZsKS0oKOAzxrj169eb37lzJ8LIyEhVs3dMyLth9W09N29vby44OLi2wwCKc6HaMxq82DNYrfgAIY5T8OsQL+hqqqEkLg6JM2aiODISOl27wGLefKh0DRFxIRHNutqBL+AhOTYH+mZCCHVKh4lkF2dj3/19uJ91Hz+0/wEAEJQYhGamzaCjTkNGCCHvhjF2g+M473epIzQ09KGnp2e9XirA2traPTg4+I6lpaWitmMhpCqEhoaaeHp6iio6Ri15laWpB96w3eD+nYn22Rr45W46+q25hA0jvSGyt4do9y5kbt2KtFW/IbZ3b5jP/hLNBwwAYwwqpQon/oxASZECzbvbw6OzLQw0DTDWfezT6rOLszE1cCrUeer40OlDDHUdCns9+1q8YUIIIYTUJzQm713w1cD6rITHiB+x7RMfGOfdxcjVxxF0Px1MIIDxJ5/A4eA/0HRxQf5/QWCMAQB4fB56T/WEldgQV/6JxV/zLuN20GOolP/vCTDQNMD2ntvRxa4L9kTtQZ+/+2DymcmIzX6hZ4EQQsgbSkxMDKdWPNJYUJL3rhgDGEMbe138pf0rdvDmYv7mQ9h88QE4joO6SAS7gK2w+v47AIAs9gEyAwJgZC5Er0ke6DezOXSNNXF2+10k3Mt6pmo3Yzd83+57nBp4ChM8J+BOxh2o80u7d5MLklEoL6zx2yWEEEJI/UBJXlURaEDgvwnW6kU4qLkQh4/8g6/2h0OmUILxeOBpl66dl3PoIFK+X4q4ocMgi46GlbMB+n/RAh983gy2rkYAgDuXHiMpOvtp1SZCE0xqNgmnBp6Cja4NAGDxlcXotq8bfg7+GY/zH9f47RJCCCGkbqMkryrZ+4KNOwMdA2Ps0fweeTf3YtifV5GWJ3taxPSzz2D1048oiYvDg379kbZmDSCXw0Zi9HS83o1jcTiw/CaOrg1DZtL/Z/fzefyn3491H4tWlq0QcDsAPQ70wPSz0xGaFlqjt0sIIYSQuouSvKpm7Aj2yWkIbLwwz/42Ih5n44PVQYhIzAEAMMag36cPHP49At1u3ZC+6jdkbN369HQen4dBc33Qqq8DEu5lYdeiqzi77Q4KcmTPXMbLzAs/d/wZx/ofw2i30biWfA1Xk64CABQqBWTKZ8sTQgghpHGhJK86aBsDIw/Bcsw27JvQFvpcLvzX/oeAyw9RtmSNwNgY1j+vgO36dTAaNgwAIHvwAKqiIqhp8OHdU4QRi33h3tEG966moChPXuGlLHUsMb3FdJz2P43hrsMBAKfjTuO9fe9hdchqpBam1sw9E0JIDZg9e7aFk5OTm1gslkokEmlgYKA2APj4+LhcuHBBq7quW1RUxNq0aSOWSCTSP//807AydYSGhmr4+Pi4SCQSqYODg9uQIUPsAeDSpUvC3bt361dFnHl5ebyOHTs6NWnSxM3Jyclt0qRJ1lVRL6mfaAmV6qJWuphxU/MSHNb/CZF8bQw9+CmC7qfjx4EeMNAqnUCh06EDAIBTKJAwcRI4TgXLxYuh7eMDoa462g0Sw7uX6Ol6ehd2R0HPWBPuHWzAV/t/ji4UCJ9+b61jDQ8TD/wR9gc2hm/Ee6L3MNx1ONxN3Wvq7gkhpMqdPn1a+8SJEwbh4eG3hUIhl5SUJJDJZKwmrn3p0iUtALh79+7tNz1HoVBAIPj/n9nJkyfbTZs2LWX48OHZAHDt2jUhAAQHB2sFBwdrDxo0KKcqYp05c2ZKnz598oqLi1nbtm3Fe/bs0fvoo49yq6JuUr9QS151E6hD0Ho8PIpv4LzJj7h97y56rvwP1x9mPlOMCQSwWLgAUHGIHzkKSQsXQpmfDwBPEzylUoWclEJc3BeNvxZewb2ryeBULy5m7W7qjt+6/IYj/Y5gsGQwziecx5yLc562Ita3BbAJIQQAEhMT1YyMjBRCoZADAEtLS4VIJHqhm2P9+vVGYrFY6uzs7DZx4sSnLVlaWlpen3zyia2Tk5Obr6+v+PHjxwIAiIyM1GjXrp2zm5uba4sWLVxCQkKe2XIoMTFRMGbMmCbh4eFaEolEGhkZqXHw4EFdV1dXqVgslvr7+4uKiooYULrY8sSJE62lUqnrpk2bnmnxS01NVbO3ty8pe+3j41NUXFzMli5danX48GHDslbC3Nxcnr+/v8jd3d3V1dVVun37dgMAWLVqlXGXLl0cfXx8XOzt7ZvOnDnT8vl719XVVfXp0ycPADQ1NTkPD4/CR48eqb/DYyf1GO14UVOiTwN7RqFEoIPxqq9wIccM07uKMamTE/i8/38QVRUVIW3Vb8jcuhUCMzPYbwuAuq3tM1U9up2JS39HI/1RPkxsddB1tBTG1i/fFaNAXoDE/ESIDcUolBdi2NFh6NmkJwaKB8JQs1K9DoSQeqaqd7yYd3GebXRWdJV2jzoZOhUubrv40cuO5+Tk8Fq1aiUpLi7m+fn55Q4ZMiSzV69e+UBpd+3y5csf2dnZyX19fSU3bty4Y2pqqmjXrp148uTJqSNGjMhmjLVYs2bNg4kTJ2bOmjXLMjU1VS0gICDe19dX/Mcff8S5u7vLAgMDtb/55hvrK1euRJW/9pEjR3RXrFhhfvbs2ejCwkLm4ODgfvLkyXseHh6yfv36iby8vArnz5+fam1t7T5mzJjUJUuWpDwf/8qVK43nzp1r6+XlVdClS5ecyZMnZ5iYmChXrVplHBwcrB0QEBAPAFOmTLGWSqVFkyZNykxPT+d7e3u7hoWF3d6yZYvh4sWLrcPDwyN1dHRUXl5e0o0bNz5o3759hetppaen85s1ayY9efLkPalUWlJRGVL/vWrHC2rJqylOXYGPj0Odz8NG4x3o7W6JFaeiMHzDVaTkFj8txhMKYT77S4h27YR2mzZQs7ICAHCq/y+UbCs1wkdft0S3j6XgVBw0ddQAAEp5xdsqaqtpQ2woBgBkybJgKjTFqpBV6LavGxZcWoCorKgKzyOEkLpEX19fFRERcXv16tVxpqamilGjRjmuWrXKuHyZoKAg7datW+dZWVkp1NTUMGjQoMzz58/rAACPx8PYsWMzAeDjjz/OuHbtmk5OTg4vJCREx9/f31EikUgnTZpkn5qaqvaqOEJDQzVtbGxkZXvhjh49OiMoKEi37PjIkSOzKjrvs88+ywgPD4/s379/5oULF3RbtmwpKWsBLO/cuXN6v/zyi6VEIpH6+fm5yGQyFh0drQ4Afn5+uRYWFkodHR2uV69eWefOnavwE75cLkf//v0dxo8fn0IJXuNFY/JqkoU7MO4M+ColVurboJ2zEeYfuoseK//DCn9PdJKYPS0q9PCA0MMDAKDIyEDcyFEwmTQRej17gjEGxmMQ+1jAuaU5GGPgOA4HV4ZAx1ATrfo6QN9UWGEI1jrW+OO9PxCTHYO/7vyFwzGHceD+AfzzwT9wNHCskcdACKn/XtXiVp0EAgF69+6d17t37zwPD4+ibdu2GU+bNi2jMnUxxqBUKqGrq6t4m7F2r6Orq1vxJ24AIpFI/vnnn2d8/vnnGc7Ozm7BwcEv/LLmOA779u2L9vT0fGaZhKCgIO2ynZPK30NFhg4dKnJwcCieP38+zb5rxKglr6bpWQEGtmAcB//YebjkfR7mOmoYs+U6lhy5jRLFi78bVAUF4Glp4fHMWUiYPAXylP/3ApT9A1epOFg5G+DBrTTsWHgF/+2OQlHeyz+8ORo4Yr7vfJz2P43FbRc/TfB+v/U7tkZuRW4JjdElhNQtoaGhGuHh4Rplr0NCQoQ2NjbP/KJr165dwdWrV3WTkpIECoUCe/fuNerYsWM+AKhUKmzevNkQALZs2WLs4+OTZ2RkpLKxsSkpGz+nUqlw+fLlij8lP+Hp6VmcmJioHhERoQEAAQEBxu3atct7Xfz79u3TK5soEh8fL8jOzubb29uX6OnpKfPz85/+Pe7UqVPuihUrzFVPenAuXrz4NJ6goCC9lJQUfn5+Pjt69KhBhw4d8p+/zrRp06xyc3P5GzdurJVEnNQdlOTVFk4FaJvAMGQNjlhtxsetLLAh6AEGrruEh+kFzxRVt7ODaNdOmM2ejYJLlxDbqzey9ux5ZgIFn89D6w8cMXyxLyRtLBF+PhHb5l1GUsyrJ2vpa+jjQ6cPS0PiOISlhWF58HJ03dsVS64swYOcB1V+64QQUhm5ubn8kSNHNnF0dHQTi8XSu3fvCn/44Ydntvyxt7eXL1iwILFDhw5iV1dXN09Pz4Ky2axCoVB17do1bWdnZ7cLFy7oLl26NAkAdu7cGbt582YTFxcXqbOzs9v+/fsNXhWHlpYWt27duof+/v6OYrFYyuPxMGvWrLTXxX/8+HE9FxcXNxcXF2m3bt3E3377bYKdnZ2iR48eeVFRUcKyiRfLli17rFAomEQikTo5ObnNnTv36eQRDw+Pgr59+zq6ubm59enTJ+v58XgxMTFqv/32m+X9+/c13dzcpBKJRPrzzz+bvPlTJg0JTbyoTRwHXPoNODUPsG2FM81+xYwjCVCqOHzXryk+aPbi8kYl8fFImjcffAMD2Kz89aVVZyUXIORUPNoPEkOgzkd2aiF0jTXB578+r7+dcRt/3fkLxx4cg1wlx5xWczBYMvhd7pQQUsuqeuJFfaSlpeVVWFgYUttxVNbzEzQIAWjiRd3FGNB2GuC/FUgKRZebU3F0mh8kFrr4bNctfLE3FIUlimdOUbezg92WzbD6/jsAgCw2FhmbNoNTKp8pZ2ihjc4jXCFQ50OpUOHwqlvYtegaYkJSX7uEitRYiu/8vsPJgScxqdkk+Fr5AgDC08Kx8+5OFMornMhFCCGEkDqEWvLqikfXAJUSsPeFQqnCyjP3sfpsNBxMtPHbkOaQWulVeFrqypXIWLsOmu7usPxuCTTF4hfKcByHh2HpuPx3DLKSC2HeRA9t+jvBytngrUL8OfhnbI7cDF01XfRz7ochkiGw0bWpzN0SQmoYteQR0jBRS159YOsD2Je2mAku/YKZVhH465NWyCtW4MM1F7Gt3JZo5ZlOmwarFcshT0jAgwEDkfbbanAlz064YIyhiacpBs/zQacREuRnFuPvFTeRGFXhLP+Xmt5iOrb12Ia21m3x152/0OvvXvj28reVv2dCCCGEVBtaQqWuUcqB+6eB+Eto02UBjk2bhJn7wjDvYCSCotPxw4D/b4kGlCZw+r16QbtNG6R8vxTpv/8OpqEBk/HjXqiax+dB2tYKzi3Ncf9aCqycDAAAD8PTYWKjAx1DzRfOKY8xhmZmzdDMrBmSC5Kx594e6GuUbreoVClx9MFRdLPvBk3Bq+shhBBCSPWj7tq6SF4MHJwEROwHmo+CqsdybLqSgB+O34WpjgZWDfGCt8iowlPz/wuClncL8IRCyGIfQM3SAjzhy1cDUMpV2PrNRZQUK+HZ2QbN37eHhtYr1wGt0MXEi5hwegIMNAzgL/bHIJdBMNc2f+t6CCHVg7prCWmYqLu2vlHTBPpvANrNBG5uBW/PcIz1a4J9E9pAwOdh0B9X8NuZ+1BWsG+tTjs/8IRCcAoFEiZNQmzfD1Bw5epLL8VX42HgbG84NjfFzZPx2Db3MkJOxUMhV770nIq0sWqDTe9vQnOz5tgQvgHd93fHl+e/RI6sSvbbJoQQQshboiSvruLxgC7zgb6/Aa69AcbgaWuAf6f5oddLtkQrjwkEsFj0LcBjiB89Gknz5kOZV/FanXomQnQb44aPvmkJc5EeLu2PRkZCQYVlX4YxhpYWLbGy80oc7X8UQ12HIi4vDrrqpTv93Mu8B7nyhX3ECSHkrcyePdvCycnJTSwWSyUSiTQwMFAbABYtWmSWl5dXJX/T1q5dayQWi6VisVjq5eUled3iyITUVdRdW59EnQB0LcFZuGPvjQQsOBgJoTr/hS3RylMVFSFt9Wpkbt4CgYkJ7Ldvg7qd3Ssvk56QDxOb0u0Qbxx/CBMbXdi5Gb10+5yX4TgOjDEUK4rRbV83qPHU8JHLR/AX+8NYaPz6CgghVaYhdNeePn1ae9asWbaXL1++JxQKuaSkJIFMJmMikUhubW3tHhwcfMfS0lLx+ppe7dSpU9rNmjUrNjU1Ve7Zs0dvyZIlVmFhYXer4h4IqWrUXdsQKBXAiTnA5h5g0afxkbctDk9tCzNdjVduicYTCmH+xRcQ7d4N7fbtoGZdusAyp3rp1opPEzyFXIm7l5NxZHUoDv4agpSHb7fVWVlSqM5Xx/d+30NsKMbvt35Ht33dMCdoDh7mPHyr+gghjVtiYqKakZGRQigUcgBgaWmpEIlE8iVLlpilpqaqdejQQdyqVSsxABw4cECvWbNmEqlU6tqjRw+HnJwcHgBYW1u7T5gwwUYsFkvd3d1dy7YmK69bt24FpqamSgDo1KlTQXJysvrzZQipD6glrz7JTQJ2fASkRAI9fwJafoJiuRLf/XsH267EwcNGH6sGe0Fkov3KahQZGYgbPgImkyZBr3evV7bQKRUqRP73GMFHH6AoTw6nFmZoM8AJukaVm0EbmxOLHXd24FDMIazrug7NzZsjryQPQoEQAh5N9iakulRHS96QI0Ncni/Txb5L5lj3sWkF8gLe2BNjnZ8/3tuxd/ow12EZaYVpgmmB0xzLH9vZe+e9V10/JyeH16pVK0lxcTHPz88vd8iQIZm9evXKB0qTt7KWvKSkJEGfPn0cAwMD7+vp6anmzJljIZPJ2PLly5Osra3dhw8fnvbDDz8kr1692nj//v2GZ8+ejX7ZNefPn29+7949zd27d8e9yTMipKZRS15DoWcJjDkGOHUF/p0BnJwHTQEPiz9sinXDW+BhegF6/xaEg7cSX1mNqrAQfD09PP7iCyRMmAh5UtJLy/IFPHh0ssHwRb7w7iXCo7uZT9frq8wHBAd9B8xtPRdn/M/Ay8wLALDy5kr0PNATmyM200QNQshL6evrqyIiIm6vXr06ztTUVDFq1CjHVatWvTD249y5c9oxMTGaPj4+EolEIt21a5dxfHz809a4UaNGZQLAuHHjMkNCQnRedr3Dhw/rbt++3WTlypUJ1XNHhFSvam3JY4x1B7ASAB/ABo7jlj133A7AVgAGT8p8xXHc0VfV2ahb8sooFcDx2QBYaYvek5a4hKxCfLbrFm7EZcG/hQ2+/cANWuoVt45xSiWy/voLqb/8CsbjweyLWTAYNOi14+4UJUoI1PkAgGPrw2FsrYNmXW2hrln5VrgLCRewJXILridfh1AgRB+HPhjmOgwOBg6VrpMQ8qyGMCbveZs3bzbctm2bcWBgYHT5lrwdO3bo79y50+jw4cMPnj/H2tra/cyZM/ckEkmJTCZjFhYWHllZWaHPl7t69apw4MCBjv/+++99Dw8PWc3cESFvr1Za8hhjfAC/A+gBQApgCGNM+lyxuQD2cBznBWAwgDXVFU+DwhcAPZcDPX4sTfDS7gEF6bAx1MLu8a0xpZMT9t1MQJ/fgnD7ccXj6BifD6ORI+Fw6CA0PdxRcPXqG02sKEvwFCVKMADXjzzA9nmXEX4uAUrly8f5vUp7m/bY9P4m7OuzD91F3fFP9D/YdmcbgNLWQhVXuXoJIQ1LaGioRnh4+NMxdCEhIUIbG5sSANDW1laWjbvr2LFjQXBwsE7ZeLvc3FxeWFjY0/MCAgKMAGDjxo2GXl5eLywlcP/+fXV/f3/HTZs2PaAEj9Rn1TkIygdANMdxsQDAGNsF4AMAt8uV4QCUbcqqD+BxNcbTsDBW+qVUALuGAioFMGwfBCbOmPW+C9o4GuOz3bfw4ZqL+PJ9F4xp2wR83otJnLqtLew2bQJXXLoUiyw2Fvlnz8Jo1Cgwwct/PATqfHT/1B3JsTm4/HcMLuyKQuiZR+gxwR3G1i/t/XglFyMXLGq7CNNbTIdCVTpBLjQtFKeWToGntTd8x34DXYOKZxETQhq+3Nxc/rRp0+xyc3P5fD6fE4lEsq1bt8YBwKhRo9K7d+8uNjc3L7l69WrU+vXrHw4ePNihpKSEAcCCBQsSyxK2rKwsvlgslqqrq3O7du2Kff46c+fOtczOzhZMnTrVHgAEAgEXERFxpybvlZCqUG3dtYyxgQC6cxw39snrEQBacRw3pVwZSwAnARgC0AbQleO4GxXUNR7AeACws7NrERdH41+f8eg6sHNwaaI3eAcgagsASM+XYfa+MJy5m4rmdgb4caAnnMxenYClrVqF9DVrodm0KSy/WwJNlxfGVb+A4zjERWTg1ul49JzoAXVNAWSF8krtnPG8kJSbSBw/EY73clGoASR1lML90y9hL231znUT0pg0xO7ayqjKpVYIqQvq8sSLIQC2cBxnA6AngG2MsRdi4jjuD47jvDmO8zY1Na3xIOs825bA2NOAtimw7UMgbC8AwERHAxtGeeOXQZ6ITS9Az1X/Yc25aChe0a1qMnUqrH/9BfKkJDwYMBBpq1ZBVVLyysszxiByN8GH05tDXVMAlVKF/T/ewJHVochIzH+nW/Myb47eB69Csf57JHlaw/7UbeQPGI24iRNQcPlypSZ/EEIIIY1BdSZ5iQBsy722efJeeZ8A2AMAHMddBqAJwKQaY2q4jJoAn5wEbHyAkADgyTp4jDH087LBqekd0EVihh+P30O/NZdwJ+klY/UYg1737nA4chj6vXoifc1aZG7Z+lahcBwgaWOJ5Ngc7FpyDWe23EZeZsU7c7wp9w790DvgNEyP7YdsWG/IQsMRP+ZjXOzigzO/fYWCvMx3qp8Q0jgkJiaGUyseaSyqs7tWACAKQBeUJnfXAQzlOC6yXJljAHZzHLeFMeYK4AwAa+4VQdHs2tdQyABFMaCpDxRlA+raAP//3aZHw5Mw/2AEsgvlmNzJCZM7OUFd8PJcP//iRWg1bw6eUAhZbCzULCzA09J6o1CKC+S4eTwOYWdLVx/o/0VzmNnrveasN6OSyZDw9y5E//ErLB8XI1/IkNrVE14TvoaVo0eVXIOQhoS6awlpmGqlu5bjOAWAKQBOALiD0lm0kYyxRYyxvk+KzQQwjjEWCmAngNGvSvDIGxBolCZ4KiWwYxDw10Cg+P9rz/V0t8TJ6R3Q28MSK8/cR9/VQQhLyH5pdTpt24InFIJTKJAwcRJi+36AgkuX3igUTW01tBnghGGLWsOzqy1MbEv3sU2Ny4WiRPlOt8nT0IDd4FHoePoGilbOQZrEDKIjt5DZexAiJo5G4c2b1JVLCCGkUaMdLxqykL+Aw9MAEzEwdA9gYPvM4dO3UzDnn3Ck5ckwvr0jPu/qDE01/kurK7x+HUlz56EkLg76A/rD/MsvwdfXf6uQ5CVKBHx9CQJ1Hnz6NIFLa0vwKpj1WxkJUTcR+ecKiM7fhyo3D0VO1sj9wA++w2dBQ1i5Gb+ENBTUkkdIw/SqljxK8hq62HPA7pGAmiYwdDdg5fXM4ZwiOb7/9w52Bz+Cg6k2fhzgAW+R0UurUxUXI/33NcjYtAl8I0OItm+Hur39W4WUeC8Ll/6OQerDXBhZacO3nyPsmxq/0Tp9b0JVWIicQ4cQufZHGKcUIUeHh6z3vdFiwjcws339bGFCGiJK8ghpmOry7FpS3Rw6Ap+cAPjqwD+Tn07IKKMvVMMPAz2w7RMfyOQq+K+/jG8PR6KwpOJxyTxNTZjNnAHRnt3Q7dgJajY2AEp30HhT1i6GGDi7Bd4f1xRKhQr//h6GtPi8St/iCzFqacFw8GC0ORuMvB+mI8veAKL915Dy/ocImvARiiIiX18JIaROmj17toWTk5ObWCyWSiQSaWBgoDYA+Pj4uFy4cOHNBgxXQlFREWvTpo1YIpFI//zzT8PK1BEaGqrh4+PjIpFIpA4ODm5DhgyxB4BLly4Jd+/e/XbdIq/Qrl07ZxcXF6mTk5Pb0KFD7RQKmmfSWNGO8I2BmSsw9gxQkg/weKWJHu/Z/L6dsylOTG+PH4/fxeaLD3H6Tgp+6O+BNk4VT3YWurlBuHgRAECRno64YcNhMmki9Pr2faMWOcYYnFqYoUkzE8RHZDydkHH3ShLMRXowtNB+x5sGeDwefD4YD3wwHg8iLiFy/U9wCLqPh+cGgu/phuRe3vAd/DnU1DXf+VqEkOp3+vRp7RMnThiEh4ffFgqFXFJSkkAmk1VNF8BrXLp0SQsA7t69e/t1ZcsoFAoIyi0qP3nyZLtp06alDB8+PBsArl27JgSA4OBgreDgYO1BgwZVyebdBw8ejDEyMlKpVCr06NHDcdOmTYbjx4/Pqoq6Sf1CLXmNha45YOxYur7J4WnA8W9KJ2eUo6MhwKIPmmL3+NbgM4ahG67im7/DkVcsf2XVquJi8I2M8Hj2V3g0/lPIH7/5xiV8Pg9NPEvXPpTLlAjaex87F13DuR33UJBTdbsJNWnaBr1/+xsuFy7A/OuvkJ/0CKbfb8V1P28cWzQOmcm0wDYhdV1iYqKakZGRQigUcgBgaWmpEIlEL/yCWr9+vZFYLJY6Ozu7TZw40brsfS0tLa9PPvnE1snJyc3X11f8+PFjAQBERkZqtGvXztnNzc21RYsWLiEhIc988ktMTBSMGTOmSXh4uJZEIpFGRkZqHDx4UNfV1VUqFoul/v7+oqKiIgaULrY8ceJEa6lU6rpp06ZnWvxSU1PV7O3tny486uPjU1RcXMyWLl1qdfjwYcOyVsLc3Fyev7+/yN3d3dXV1VW6fft2AwBYtWqVcZcuXRx9fHxc7O3tm86cOdOyoudkZGSkAgC5XM7kcjmrqqEwpP6hMXmNjUoFHP8KuLYekPQG+v8JqL/Yw1FUosQvp6Ow4b9YmOtp4vv+7ujk8vItxTiVCll/7UDqL7+AATCdOQOGQ4e+9Ti7wtwSBP/7AJH/PQZPwNCsqx28utlBXVi1jc4KeQmuHViL3L92wj4qBzIB8LitEzrPWIG7hoUITgmGt7k3mpk1q9LrElJbqnpM3uNv5tjK7t+v0u5RDWfnQqvvv3v0suM5OTm8Vq1aSYqLi3l+fn65Q4YMyezVq1c+UNpdu3z58kd2dnZyX19fyY0bN+6Ympoq2rVrJ548eXLqiBEjshljLdasWfNg4sSJmbNmzbJMTU1VCwgIiPf19RX/8ccfce7u7rLAwEDtb775xvrKlStR5a995MgR3RUrVpifPXs2urCwkDk4OLifPHnynoeHh6xfv34iLy+vwvnz56daW1u7jxkzJnXJkiUpz8e/cuVK47lz59p6eXkVdOnSJWfy5MkZJiYmylWrVhkHBwdrBwQExAPAlClTrKVSadGkSZMy09PT+d7e3q5hYWG3t2zZYrh48WLr8PDwSB0dHZWXl5d048aND9q3b1/4/LX8/Pycw8LCtDt06JDz999/PxC8YptKUr/RmDzyfzwe0PNH4P2lwN1/ga29gfzUF4oJ1fn4pqcr9k9sAx0NAcZsvo6Ze0KRXVjx7heMx4PRiOFwOHQIQi8vFAYHV2oihZaeOtoPccGQBa0gcjfBzRNxKMx99Y4blSFQU0ebQZ+h+6Er4G1fhYQ2DrC79AAnJn2AT46MxG83V2HcybG4lXqryq9NCKkcfX19VURExO3Vq1fHmZqaKkaNGuW4atUq4/JlgoKCtFu3bp1nZWWlUFNTw6BBgzLPnz+vA5QO4Rg7dmwmAHz88ccZ165d08nJyeGFhITo+Pv7O0okEumkSZPsU1NTX7knY2hoqKaNjY2sbC/c0aNHZwQFBemWHR85cmSFXaOfffZZRnh4eGT//v0zL1y4oNuyZUtJWQtgeefOndP75ZdfLCUSidTPz89FJpOx6OhodQDw8/PLtbCwUOro6HC9evXKOnfuXIVLBwQFBd1PTk4OLSkp4R0+fLhqFigl9Q6l9o2V76TSJVX2jwMCPgQm/AfwXlw+xcvOEEem+WF1YDTWnIvBhftpWPJhU7zvZlFhteo21rDd8Ce44tIdLmQxMcgLDITxmDFgb/FJ0sBcC++Pa4q8zGLoGpX2nPy3OwoWjvpwam4GVkXLrgCAi3c3uHh3gyIrC2cOzIaCXYYKQIlchgvH/oD7gB/B16ElWAgp71UtbtVJIBCgd+/eeb17987z8PAo2rZtm/G0adMyKlMXYwxKpRK6urqKtxlr9zq6urov3TtSJBLJP//884zPP/88w9nZ2S04OFj4fBmO47Bv375oT0/PZ8asBAUFaT//4flVH6a1tLS4Pn36ZP/9998G/fr1q3ibI9KgUUteY+baBxj9L/DeogoTvDIaAj5mvueCg5PbwkRHA59uu4EpO24iI7/iMXOMMfCEpb+3co8eQ9qKn/Hwo0EovnPnrUMsS/BKihVIjMrGyQ2R2LssGAl3q34bM4GhITr0mQR1NQ3wwCDgGOwDziK6fQckL/kOJQ8fVvk1CSFvLjQ0VCM8PFyj7HVISIjQxsbmmab+du3aFVy9elU3KSlJoFAosHfvXqOOHTvmA4BKpcLmzZsNAWDLli3GPj4+eUZGRiobG5uSsvFzKpUKly9ffiHxKs/T07M4MTFRPSIiQgMAAgICjNu1a/faJQL27dunVzZRJD4+XpCdnc23t7cv0dPTU+bn5z/9e9ypU6fcFStWmKuerIZw8eLFp/EEBQXppaSk8PPz89nRo0cNOnTo8MwG4Tk5Oby4uDg1AJDL5Th27Ji+RCIpel1spGGiJK+xs2kBOHUt/f7an0Do7pcWbWqtj0NT2mJmNzFORCaj2y8XcCj08St3ljCdOgXWK1dCnpqKBwP9kfrLr1DJ3n5ChbqmAB/NaYkuo11RlF+Cg7/ewuFVt5CbUbW/u5qZNcOf723A1ObTsLF3AN5buQ+63boia/duxHTvgfhPP0V+0EXaTYOQWpCbm8sfOXJkE0dHRzexWCy9e/eu8Icffnhmppe9vb18wYIFiR06dBC7urq6eXp6FpTNZhUKhapr165pOzs7u124cEF36dKlSQCwc+fO2M2bN5u4uLhInZ2d3fbv32/wqji0tLS4devWPfT393cUi8VSHo+HWbNmpb0u/uPHj+u5uLi4ubi4SLt16yb+9ttvE+zs7BQ9evTIi4qKEpZNvFi2bNljhULBJBKJ1MnJyW3u3LlPJ494eHgU9O3b19HNzc2tT58+Wc+Px8vNzeX16tXLSSwWS6VSqZuJiYn8iy++eG1spGGiiReklEoFbO9XunhypzlA+y+AV3QDRKXk4Yt9YQh9lI1uUnN892FTmOm9fCkSZXY2Un78CTkHDsB0xgyYjB9X6VAVciXCzyUi4nwC/L9uCU1tNahUXJXtnFHhNdPSkLVrN7J27YIyIwPqjo4wGjEc+n37vvFevoTUJloMuXR2bWFhYUhtx1FZz0/QIASgiRfkTfB4wNC9gOcQ4Ox3wMHJgOLlEx7E5ro4MLEN5vR0xYWoNHT9+Tz2Bj96aQsX38AAVt9/B7stm2E0YjiA0vF6qoKCtw5VoMaHVzc7DFvkC01tNXAqDgd+uoGgffdRnP/q5V4qS2BqCtOpU+B0NhCWy5aCp6GB5IXf4n7HTkj56SfIExOr5bqEEEJIZVFLHnkWxwHnfwDOLS3dLWP4gVeO1wOAB+kFmL0vDNceZqK92BRL+7vD2uCVQ1rAKRSI7dUbnFwOi0WLoOPXttIhlxQrELTnPu5eToKapgDN37eDR2dbqKm/Ou53wXEcikJCkBmwDXmnTgEcB90uXWA0cgSE3t5VtkUbIVWFWvIIaZho71ry9m7tBAozgDZT3qi4SsVh+9U4LDt2FwzA1z1dMdTH7pVdqIU3byJp7jyUxMZCv18/mM/+EnwDg0qHnJGYjysHY/EwLB3aBhroO60ZjKzefeeM15EnJSFrxw5k79kLZU4ONFxdYTRiBPR69QRPQ+P1FRBSAyjJI6RhoiSPvJsH/wHq2oB189cWfZRZiK8PhCMoOh2tHYyw5EN3OJm9fPkRlUyG9LVrkfHnBvANDSH6azvU7e3fKdzH97MQfj4RXUdLwRfwkJdZDB1DjWpvXVMVFSHn8GFkbdsG2f1o8I2MYDh4EAwGDYaa+csXkiakJlCSR0jDREkeqTyVCvijPZARAwzYAEh6vfYUjuOwJ/gRlvx7B0UlSoxpK8K0Ls7Q1Xz5+qLFd+4ga9duWCyYD8bjgVMqwfjv3t0qlymxff5l6JsI4dvfCZaOVbYH+EtxHIfCK1eQGbAN+efOAXw+9Lp3L+3K9fCo9usTUhFK8ghpmGjiBak8Hq90XJ6pBNg1DLiy7rWnMMYwqKUdzs7qiIEtbLAh6AE6LT+PfTcSoFJV/KFC09UVlt8uBOPxoEhPR0yPnsg+8Pc7L1XCFzC07NUEOWlFOPDTDRxdG4as5Lef7PE2GGPQ9vWF7do1cDxxHEbDhiL/7Fk8/GgQHg4ajJx//wUnr54JIoQQQkgZSvLI6+mYlS6aLOkFHJ8NHJsNqJSvPc1ERwPLBnjgn0ltYWskxKy9oei/9hJCH2W/8jxVsQwCU1MkffMNHn0yFiUJlZ+5yuPz0LS9NYYv9kWrvk2QcC8LOxddQ0Zi/utPrgLqdnYw//prOJ0/D/M5c6DIzsLjmbMQ3aUr0tetgyKz6hd1JqQhmz17toWTk5ObWCyWSiQSaWBgoDYALFq0yCwvL69K/qZt377doKz+pk2bup44cYK2vCH1EnXXkjenUgIn5wIF6UD/P165jt4Lp6o4/B2SiKXH7iKjQIaPWtjii+4uMNGpeGICp1Iha9cupC1fAQ6A2eefw3DE8HceV1eUV4Koaynw6GwDxhge3c6EWRM9aAhrZoc/TqVC/oULyNq2HQUXL4Kpq0OvT2+YzZgBgbHx6ysgpJIaQnft6dOntWfNmmV7+fLle0KhkEtKShLIZDImEonk1tbW7sHBwXcsLS0V73qdnJwcnq6urorH4+Hq1avCwYMHOzx48CCyKu6BkKpG3bWkavD4QPelQL91pQle1kMgL+XNTuUxDGhhg7OzOmBcOwfsv5mATsvPYWPQA8iVL27zyHg8GA0dCofDh6DVogWKbt2qkokTQl11eHaxBWMMskI5jq4Lw/a5lxF65hGU8pduN1llGI8H3Y4dYbdxAxyOHIZ+/34ouHSZFlQm5A0kJiaqGRkZKYRCIQcAlpaWCpFIJF+yZIlZamqqWocOHcStWrUSA8CBAwf0mjVrJpFKpa49evRwyMnJ4QGAtbW1+4QJE2zEYrHU3d3dtWxrsvL09fVVPF7pn8e8vDweLYlE6itqySOVw3HAHx2Awixg2F7ATPJWp0en5mPRkdu4EJUGZzMdLOzrhrZOJi+5FAdOJgNPUxOy6GjknT4N408+AVN7+USON5Ual4vLf8cg4W4WdI010aqvA8QtzcGqcfeM53FyeZXcCyGvUh0teQ/8P3J5voxut66ZJuPHp6kKCnhxo8c4P39cv2+fdKMRIzIUaWmCR5MmO5Y/1mTvnnuvun5OTg6vVatWkuLiYp6fn1/ukCFDMnv16pUPlCZvZS15SUlJgj59+jgGBgbe19PTU82ZM8dCJpOx5cuXJ1lbW7sPHz487YcffkhevXq18f79+w3Pnj0b/fy1AgICDBYsWGCdmZmptn///vtdu3at3sG8hFQSteSRqscY0GcloJQBG98DYs+/1elOZjrYOqYl/hzpDZlChWEbrmLCtht4lFn4QlnGGHiapVum5R4/gbRfV+KB/0coinz33hMzez30/awZ+kzzhIaWAKe33EZW8osxVCdK8Ah5M/r6+qqIiIjbq1evjjM1NVWMGjXKcdWqVS+Mczh37px2TEyMpo+Pj0QikUh37dplHB8fr152fNSoUZkAMG7cuMyQkJAKx9uNHDky+8GDB5G7du2Knj9/vnVFZQip66glj7yb7Hjgr4+AjPtA39+AZkPfuopiuRIb/ovF72djoOI4TOjgiAkdHCF8yY4VeadPI+nbb6HMzILxx2NgMnny0yTwXXAqDskPcp8us3LzZBxsJUYwtdN957oJqW0NYUze8zZv3my4bds248DAwOjyLXk7duzQ37lzp9Hhw4cfPH+OtbW1+5kzZ+5JJJISmUzGLCwsPLKyskJfdR0bGxv369evV8l4P0KqGrXkkepjYAd8fBywbwNc+xNQvv3vQE01PqZ0dsaZmR3QTWqOlWfuo+vP53EsPKnCJVR0u3aF45Ej0O/3ITL+3IDMgG1VcSdgPPY0wSsukCPkZDz2fH8dJzdGIje9qEquQQipvNDQUI3w8PCnY+hCQkKENjY2JQCgra2tLBt317Fjx4Lg4GCdsvF2ubm5vLCwsKfnBQQEGAHAxo0bDb28vF7oho2IiNBQqUrH6AYFBWmVlJQwc3NzSvBIvVMzUwpJwyY0AIbtB0ryAb4AkOUDfDVA8HZbelkZCLF6aHMMb52BhYciMfGvm2jjaIyFfd0gNn+2NY2vrw+rJUug36cvhJ6lCwzL7t+HwNIKfJ1338pMU1sNwxf7IuREHELPPELMzVQ07WCNlr2aQFObulcJqQ25ubn8adOm2eXm5vL5fD4nEolkW7dujQOAUaNGpXfv3l1sbm5ecvXq1aj169c/HDx4sENJSQkDgAULFiR6eHjIACArK4svFoul6urq3K5du2Kfv87OnTsNd+/ebSwQCDhNTU3Vtm3bYssmYhBSn1B3LalaHAf85Q8oioFB2wChYaWqUShV2HEtHitORiFfpsCI1vaY3k0MfWHFCRanUCCmVy9wJXJYfrsQOu3bv8tdPCM/S4brR2IRfSMVQ79tDW192o+W1D8Nsbu2MqpyqRVC6gLqriU1hzHA3R+Iv1I6ISPrYaWqEfB5GOkrwtlZHTG4pS22Xn6ITsvPYde1eCgr2DWDCQSwWrYMPC0tPBr/KR7Png1FVtY73kwpHUMNdBrhipHft4G2vgY4jsOx9eGIuJAIVQXLvxBCCCF1ASV5pOp5DgJG/gPkpwIbugIJNypdlZG2Or7r547DU/zgaKqNrw6E48PfL+JG3IsJnJaXF5r8fQAmkyYi59+jiO3VGyUPH1b+Pp6joVXaiigrVKAwpwTnd9zDzkXXEBuS9s7brxFCakZiYmI4teKRxoKSPFI9RH7AJ6cANS1g/8eA8t32am1qrY89n/pi5eBmSM0rxoC1lzBjzy2k5hY/U46nrg7TadPQZP8+6L7/HtTs7ACUdudWFU1tNfT/ojl6THAHY8Cx9eE48NMN5GbQ5AxCCCF1ByV5pPqYioGxZ4DBO0onYnBc6VclMcbwQTNrBM7siIkdHXEkNAmdV5zHHxdiUKJ4tttU08UFlgsWgPF4UKSlIaZ7D2Tv21dlLW6MMTg0M8XgeT7oNFwCjgO0dEuX4VKUvH5fX0IIIaS6UZJHqpeOKWDuVvr9qXnAsdmle+C+A20NAWZ3l+DE9PbwaWKE74/eRfeVF3A+Kq3C8lxJCdQsLZE0dx7iP/4YJY8evdP1y+PxeZD6WWHAly0gUOdDKVdh5+JrCAy4g/ys4tdXQAghhFQTSvJIzShrQbu2Htg1DCh59x2CmphoY9Poltg02hsqFYdRm65h7NZgxGc8u2OFmrU17LZugcXChSgOC0ds3w+QuXVrlY6jK9vbUqlUoYm7Ce5dTcb2+Vdw+e9oyArfrauaEEIIqQxK8kjNYAx4bwnQczlw/wSwuSeQl1wlVXeWmOPE9PaY3V2CSzHp6PrLefx4/C5yi/+fXDEeD4aDB8Hh3yPQ9vFBUWgYqmPTcXVNAfw+csawb1vD0csUN0/EY9vcy8hOrdmt0ghpqGbPnm3h5OTkJhaLpRKJRBoYGKgNAD4+Pi4XLlzQqq7rHjlyRPfUqVNvtQjnjBkzrObPn29eXTFVxoABA0SbN29+67WtfvzxR9PVq1cbA8CqVauMHz58+HQ9K2tra/ekpKRKrbt76dIl4e7du/XLXtfEMysqKmJt2rQRSyQS6Z9//lmpdb5CQ0M1fHx8XCQSidTBwcFtyJAh9sCL91NVOnfu7OTs7Oz2tufRYsikZvmMK90lY++Y0kRv8tXS8XrvSEPAx8SOjujf3BpLj97BmnMx2HEtHhM7OGKkr+jpFmlqFhawWbcWXEkJAEAWHY3ckydhMnYsmLr6qy7xVvRMhOj2sRuadbXDvSvJ0DcVAgAyHufDyEIbjFf1CSYhDd3p06e1T5w4YRAeHn5bKBRySUlJAplMViP/mAIDA3V1dHSU3bp1e/duiHroyy+/fDoeZvv27SbNmjUrEolE79xNERwcrBUcHKw9aNCgnHetCwAUCgUEglenNpcuXdICgLt3796ubL2TJ0+2mzZtWsrw4cOzAeDatWtCoOrvBwC2bt1qoK2tXalxTtSSR2qe+H1gzFGg68IqSfDKM9fTxK+DvXBkqh88bQyw9NhddPjpLLZfiYP8yZp2jDHwNEoXNM49eRLpq37Dg4H+KAoPr9JYAMDUThd+HzmDMYbifDn2/3gDu7+/jrjIDFp2hZC3lJiYqGZkZKQQCoUcAFhaWioqSjTWr19vJBaLpc7Ozm4TJ060LntfS0vL65NPPrF1cnJy8/X1FT9+/FgAAJGRkRrt2rVzdnNzc23RooVLSEjIM5th37t3Tz0gIMB03bp15hKJRHr8+HGde/fuqbdu3VosFoulvr6+4vv371f4KTEsLEyrWbNmEnt7+6YrVqwwKXt/3rx55k2bNnUVi8XS6dOnW1V0rpaWllfZ95s3bzYcMGCACAA2bdpk6Ozs7Obi4iL19vZ2AUqTkE8//dSmrM6ffvrJBABUKhVGjhxpJxKJmrZp00acnp7+QgaUmJgocHNzcwWAy5cvCxljLcrux9bWtmleXh6vrIVt8+bNhhEREVojR450kEgk0vz8fAYAP/74o5lUKnUVi8XS558fABQWFrKBAweKxGKx1NXVVXr48GHd4uJitnTpUqvDhw8blm9Vu3PnjtDHx8fFxsbGfcmSJWZldaxZs8bI3d3dVSKRSIcOHWqveLJqgpaWlte4ceNsXFxcpGfOnNGZNGmStaOjo5tYLJaOHz/e5vl7HTNmTJPw8HAtiUQijYyM1Dh48KCuq6urVCwWS/39/UVFRUUMKG2hnDhxorVUKnXdtGnTMy1+qampavb29iVlr318fIoqup/c3Fyev7+/yN3d3dXV1VW6fft2A6C0NbRLly6OPj4+Lvb29k1nzpxpWdHPQE5ODm/VqlXmCxcuTKro+OtQSx6pHVbNSr8AIGwPoJABzUdUWfVNrfWx9WMfXI3NwE8n7mHuPxH440IsZnQTo4+nFfhPWtJMJ02CpkSC5IXf4uGgwTAaMxqmU6aAJxRWWSxlNLQE6DjMBVcPxuLIb6GwdjFEm/6OMLPXq/JrEVLdTqz91Tb9UVyVdo+a2NoXvj/x85fOjPrwww9zly5daiUSiZr6+fnlDhkyJLNXr1755cs8fPhQbeHChdY3bty4Y2pqqmjXrp1427ZtBiNGjMguKirieXt7F2zcuPHRrFmzLL/66iurgICA+LFjx9r/8ccfce7u7rLAwEDtiRMn2l25ciWqrE4XF5eSkSNHpuno6CgXLVqUApR2nw0bNixj6tSpGb/++qvxxIkTbU+fPh3zfMx37twR3rhx405eXh7fy8tLOmDAgJybN28Ko6OjNcPCwu5wHIeuXbs6HTt2TKdHjx75z59fkWXLllmePHkyqkmTJvL09HQ+APz6668m+vr6yoiIiDtFRUWsZcuWkj59+uRevXpVKzo6WiM6OjoiISFBzd3d3W306NEZ5euztrZWyGQyXmZmJu/s2bM6bm5uhadPn9bhOC7f2NhYoaur+3T5gjFjxmStXbvWbPny5Y/at2//dByKiYmJ4vbt23eWLVtmumzZMvPdu3fHlb/GDz/8YMYYQ1RU1O2QkBDNnj17OsfExER8/fXXj4ODg7UDAgLiAWDGjBnC6OhozUuXLt3Lzs7mu7q6Nv3iiy/SIiMjNfbt22cUHBx8V0NDgxs+fLjdunXrjKdMmZJRVFTEa9WqVcGff/6ZkJyczP/0009FsbGxETweD2XPp/y9rlmzJm7FihXmZ8+ejS4sLGRdunRxOXny5D0PDw9Zv379RD/99JPp/PnzUwHA2NhYcfv27TvP/z+YPHlySs+ePcVeXl4FXbp0yZk8eXKGiYmJ8vn7mTJlinWnTp1y9+7d+zA9PZ3v7e3t2rdv31wACAsL0w4PD4/U0dFReXl5ST/44IOc8s/0yfOw/uyzz1J0dHQqtfI+teSR2sVxQMR+4NAUIHDJOy2xUpFWDsbYO8EXm0e3hLaGAJ/vvoWeK//DqdspT1vSdDt3hsO/R2AwYAAyN25CZsC2Ko2hDOMxiFtaYOjC1vD7yBkZifnYuyyYxusR8ob09fVVERERt1evXh1namqqGDVqlOOqVauMy5cJCgrSbt26dZ6VlZVCTU0NgwYNyjx//rwOAPB4PIwdOzYTAD7++OOMa9eu6eTk5PBCQkJ0/P39HSUSiXTSpEn2qampr+1iCAkJ0R4/fnwmAEycODHzxo0bOhWV69GjR7aOjg5naWmp8PX1zf3vv/+0jx8/rnfhwgU9qVQqdXNzk8bExGjevXv3hdavl/H29s4fNmyYaMWKFSZlrVmnT5/W27Nnj7FEIpF6eXm5ZmVlCW7fvq15/vx53Y8++ihTIBBAJBLJfX19815W5+nTp3WCgoJ0v/zyy6T//vtP9/Tp0zqtW7d+o8Rz6NChWQDg4+NT+OjRoxf2frx06ZLOiBEjMgDAy8ur2MrKqiQ8PLzCe37vvfeyhUIhZ2lpqTAyMpInJCQIjh8/rhsREaHl6enpKpFIpEFBQXqxsbEaAMDn8zF69OgsADA2NlZqaGioBg0aJNq6davB65Kj0NBQTRsbG1nZvsajR4/OCAoKerpZ+siRIyvcOumzzz7LCA8Pj+zfv3/mhQsXdFu2bCkpawEs79y5c3q//PKLpUQikfr5+bnIZDIWHR2tDgB+fn65FhYWSh0dHa5Xr15Z586de+Zn6NKlS8IHDx5ojBw5MvtV9/Aq1JJHahdjwKDtwL8zgAs/lW6D9sHvgKDq9odljKGTxAwdxKY4Ep6En0/ew7iAYHjZGeCL913QxtEEfF1dWC5eBP2+faDp7g4AKI6KgpqlJfi6uq+5wtvhC3jw7GwLV19LPIxIh4FZaWNI1PVk2LgYQUuv6sYGElJdXtXiVp0EAgF69+6d17t37zwPD4+ibdu2GU+bNi3j9We+iDEGpVIJXV1dxduMz3rbazz/muM4fP7550lffPHFK/cBLn9u+QRix44d8YGBgdqHDh3Sb9GihfTGjRu3OY5jK1asiB8wYEBu+TqOHDnyRpMA2rVrl3fhwgXdhIQE9WHDhmWvWLHCAgDXu3fvNxpbpqmpyQGAQCDgFArFO42T1NDQePppn8/nQ6FQMI7jmL+/f8bvv/+e+Hx5dXV1Vdl4OTU1Ndy6devOoUOH9Pbt22e4du1as/Ktsm+rfCvm80Qikfzzzz/P+PzzzzOcnZ3dgoODX+gC4jgO+/bti/b09JSVfz8oKEi7op+N8v777z+diIgILWtra3eFQsEyMzMFPj4+LteuXbv3pvFTSx6pfXw1oM8qoMt8IHwvsK0foCh5/Xlvicdj6OtphVMzOmBpf3ckZRdj6J9XMWLjVYQ+ygYAaLVsCZ6mJjilEglTpyK2dx/knTtX5bEAgLpQAHFLCwBAYW4JArfexfZ5l3H93wcoKaZdlwh5XmhoqEZ4ePjTT4AhISFCGxubZ35ZtGvXruDq1au6SUlJAoVCgb179xp17NgxHygdn1Y2s3TLli3GPj4+eUZGRiobG5uSsjFXKpUKly9ffuGPta6urjIvL+9p15+Xl1fBhg0bDIHSMYDe3t4VtngdO3bMoLCwkCUnJ/OvXLmi6+fnV9CjR4/cbdu2meTk5PAA4MGDB2qJiYkvNLoYGxvLb968qalUKnHw4MGnY8IiIyM1OnfuXPDrr78+NjQ0VMTGxqp369YtZ+3ataZlE1HCwsI0cnNzeR06dMjbt2+fkUKhQFxcnNqVK1cq/NTatWvX/P379xs1adJExufzYWBgoDh79qx+t27dXrgvHR0dZU5ODr+iel6mbdu2+du3bzcqiy0pKUndw8OjWE9PT5mfn//aXKR79+65R44cMSx7TikpKfyoqKgXPhHn5OTwMjMz+YMGDcpZt27do7t3775ySIGnp2dxYmKiekREhAYABAQEGLdr167C1s7y9u3bp1f2rOPj4wXZ2dl8e3v7kufvp1OnTrkrVqwwV6lKc8WLFy8+/dkKCgrSS0lJ4efn57OjR48adOjQ4ZlnPXv27LTU1NSwxMTE8AsXLtwViUSyt0nwAGrJI3UFY0C7mYCBPZAeBQiqrzVLjc/DEB879POyxvYrcVhzLgYf/H4R77uZY9Z7LnA21wXj82H9009ImjMXCRMmQq93b5h/8zUERkbVEpOWnjoGzW2Jqwdjce3wA4SfT4RPLxFc/azA59NnMUIAIDc3lz9t2jS73NxcPp/P50QikWzr1q3PjP2yt7eXL1iwILFDhw5ijuNY165ds8tmQAqFQtW1a9e0f/rpJytjY2P5gQMHYgFg586dsePGjbP/4YcfLBUKBevXr1+mr6/vM/sUDhgwIHvgwIGOx44dM/j111/j161bFz9y5EjRypUrLYyNjRUBAQEPK4rZ1dW1sE2bNi5ZWVmCWbNmJYlEIrlIJJJHRkZqtmzZUgIAWlpaqr/++uuBtbX1M5/uvv3228QPPvjAycjISOHp6VlYUFDAA4Dp06fbPHz4UIPjOObn55fbunXrolatWhU9fPhQw93d3ZXjOGZkZCQ/evRozIgRI7LPnDmj5+Tk1NTKykrm5eVVYTLq4uJSwnEcK0twfH1985OSktRNTU1fmNU5cuTI9KlTp9p/8cUXquDg4BfGq1Xkyy+/TB05cqS9WCyW8vl8rF+//qFQKOR69OiRt3z5ckuJRCKdOXPmSycXtGjRonju3LmJXbp0EatUKqipqXGrVq2KF4vFzyT52dnZ/N69ezuVJWCLFy9+ZYuzlpYWt27duof+/v6OSqUSnp6ehbNmzap4Zf1yjh8/rjdr1iw7DQ0NFQB8++23CXZ2dorn72fZsmWPx48fbyeRSKQqlYrZ2trKzp49Gw0AHh4eBX379nVMTk5WHzhwYMbz4/GqAqvOGX6Mse4AVgLgA9jAcdyyCsp8BGAhAA5AKMdxQ19Vp7e3NxccHFwN0ZI6J+EGwKkA25bVepm8Yjk2BT3En//ForBEgQ+9rDG9qxi2RlrgSkqQ/sefSF+/HnwdHYh27oC6SFSt8STH5uDSgWikPszDsEWtoWv0xkN1CHkpxtgNjuO836WO0NDQh56enq/sYqzLtLS0vAoLC0NqOw5CVq1aZVx+gsa7CA0NNfH09BRVdKzaWvIYY3wAvwPoBiABwHXG2CGO426XK+MM4GsAbTmOy2KMmVVcG2l0OA448TWQFAr0/wOQflBtl9LVVMNnXZ0xwtcea89FY+vlOBwOfYyhPnaY3NkJZlMmQ/e9bsjevQdqdnal4cnlYGpVu/xLGQsHffSb2RzZKYXQNdIEx3G4uC8aTTxNYC2u1LqdhBBCGqFqa8ljjPkCWMhx3PtPXn8NABzHLS1X5kcAURzHbXjTeqklrxEpSAd2DsFv+feg5/w+RvZYD8ar/q7LpJwirDoTjT3Bj6DO52FMWxE+be8Ifa3SpE6RloaHg4fA+NPxMPD3r5adM8orzC3Bnu+voyBbBnt3Y/h+6Ahj6won8hHyUtSSR0jD9KqWvOr8i2kNoHxfeMKT98oTAxAzxi4yxq486d59AWNsPGMsmDEWnJb22q5y0lBom0A18h88NLHH8vQr+G5PLyjkxdV+WUt9IZb2d8eZGR3QTWqONedi0O7HQPx+NhqFJQpwcjnUbGyQPH8B4kePQUn8O7e2v5KWnjqGL2oN336OSIrOwa4l13Bm620U5Vf95BRCCCENx0uTPMbYUcaYqJqvLwDgDKAjgCEA/mSMGTxfiOO4PziO8+Y4ztvU1LSaQyJ1CU9dGz8NOYsxOs7YLUvAtPPTUSivmXXlRCbaWDXEC0entUNLkRF+OnEP7X88hx0PS2C+YQMsFn2L4shIxPb9ABmbt1TrDhYCdT6av2+PEUt80ayrHRLuZYEvKP3nSztnEEIIqcirWvI2AzjJGJvDGKvM4KNEALblXts8ea+8BACHOI6Tcxz3AEAUSpM+Qp7i8QWYMeAA5vnMwaXHlzH632FITYmosetLrfSwcXRL7J/oC0dTbSw4FIkuKy7gjGMb2B8+DG1fXxRHRlZ7ty0AaGqroe0AJwz/1hfqmgKolCoc+OkGbp6Mg0Jeqa0NCSGENFAvTfI4jtsLoDkAPQDBjLFZjLEZZV9vUPd1AM6MsSaMMXUAgwEceq7MPyhtxQNjzASl3bexb30XpFH4yHUwfuu8CnHZMRj672BERR+r0eu3sDfCrvGtEfCxD4y01TFrbyh677yHyElzYbFkMYDSBZTTflsNrqR6u1L5aqX/dGWFCqgL1XD5QAz+mn8Fdy4lQaWilj1CCCGvH5NXAqAAgAYA3ee+XonjOAWAKQBOALgDYA/HcZGMsUWMsb5Pip0AkMEYuw3gLIAvOI6r1MrlpHFoZ9MeAW2+Awdg5H9f4NL132v0+owxtBeb4tCUtlg7rDk4jsOEv0LQb0Mw/rufhrwzZ5D+++94MGAAisLCqj0eoa46+kz1xAfTvaClp47AgDvYveQa8rOqf+wiIbVh9uzZFk5OTm5isVgqkUikgYGB2gCwaNEis7y8vCodZ37+/HktgUDQomwBZULqm5fOrn0yCeJnlLa+LeI4rk5ssEmzawkAJCffwuRjoxDDlJhn/T4GdFtRK3EolCocCEnEytP3kZhdhNYORpitnwGdtSugSE2F0ciRMJ02FTytKt3HvUIcxyHmZhqiriWj+6fu4PEYivPl0NSpnqVeSP3SEGbXnj59WnvWrFm2ly9fvicUCrmkpCSBTCZjIpFIbm1t7R4cHHzH0tKySraLUSgU8PPzE2toaKhGjx6dMWbMmAr3MCWktlV2du0cAP4cx31VVxI8QspYWDTD1oHH0ZppY+Hjk1gV/DNU3Cv3oa4WAj4PH3nbInBWByzsI0V0aj76hTAs6jcPhe/3QeaWLcjctr1GYmGMwamFGXpO9ACPxyArlOOvBVdw/I9wZKfQP2FS/yUmJqoZGRkphEIhBwCWlpYKkUgkX7JkiVlqaqpahw4dxK1atRIDwIEDB/SaNWsmkUqlrj169HAo20LM2trafcKECTZisVjq7u7uWrad1fO+//57sw8++CDLxMSE9hgk9Va17nhRHaglj5SnkBfhu0vfYt/Df9HD/j0sbj0fGppvtB93tSgsUWDXtUf4879YJOUUoweXgh4DOqJnCxEU0dFQs7IEX/e1ox2qhFymRMipeISciodKroLUzwotezeBll71bRlH6q7qaMnbu/S6y/NlHLxMM1t0F6WVFCt4B38JeWEinUsri3SPzrYZBTkywdE1YY7lj/l/3fKV+3Lm5OTwWrVqJSkuLub5+fnlDhkyJLNXr175QGnyVtaSl5SUJOjTp49jYGDgfT09PdWcOXMsZDIZW758eZK1tbX78OHD03744Yfk1atXG+/fv9+wbJupMg8ePFD76KOPHK5cuXLvo48+EvXu3TuHWvJIXVVb6+QRUu0EakLMb78U01tMx7G4kxi/qzOysx7UWjxa6gJ87NcE57/ohOX+nrhv7ohpB+6gy09ncWf8RMT06o28wLM1EouaBh8+vZtgxGJfSNtZ4XbQY2ybdxl5mTRej9RP+vr6qoiIiNurV6+OMzU1VYwaNcpx1apVxs+XO3funHZMTIymj4+PRCKRSHft2mUcHx//9NPNqFGjMgFg3LhxmSEhIS+sLD5p0iTbZcuWJfD5/Oq9IUKqGbXkkQbj+IVFmBO7B5YqhjVd18LOzq+2Q4JKxeHUnRSsOReDwrBwzArdC7vsxxB27w6beXMhMH7h71O1yU4pRExIKlp0FwEAHt/PgnkT/afr7ZGGrSGMyXve5s2bDbdt22YcGBgYXb4lb8eOHfo7d+40Onz48Auf+Kytrd3PnDlzTyKRlMhkMmZhYeGRlZUV+nyZsu+zsrIEmpqaqpUrV8aNGDEiuwZui5C3Qi15pFHo3n4+NjT/CjmMw/AzE3Ar/K/aDgk8HsP7bhb4Z1IbfPvlQOz4eBG2unZH7slTiHyvB5IjomosFgNzracJXn5WMQ7+egs7vr2K+8EptKAyqRdCQ0M1wsPDn46hCwkJEdrY2JQAgLa2trJs3F3Hjh0LgoODdcrG2+Xm5vLCwsKenhcQEGAEABs3bjT08vIqeP46iYmJ4WVfPXr0yFqxYkU8JXikPhLUdgCEVCUvj+H4y0CESacn4pMbS/G9lj7ed+xd22GBMYY2jiZo42iC8J5NsXtfJ+icOozNO+/Dv6UM49vYwc6s5sYSahtooOdED1z+OxonN0Ti1ql4+PZ3go0LrRRB6q7c3Fz+tGnT7HJzc/l8Pp8TiUSyrVu3xgHAqFGj0rt37y42NzcvuXr1atT69esfDh482KGkpIQBwIIFCxI9PDxkAJCVlcUXi8VSdXV1bteuXbQ2K2mwqLuWNEjZWQ8w7fxMhOTcx/Tmn2OM2xgwXt1quI5Ny8cfF2Jx9mIkfjy7Cnc790P7mePhamVQYzGoVByiriXj6qFYFGaXYMR3vtAx1Kyx65Oa0xC7ayujqpdaIaS2UXctaXQMDJvgzz470UPUA7/c/BWL9/SAQl63Jhw4mOpg2QAP/D3BF5yNHdr/uxnh/sMw4+fDuP4ws0Zi4PEYJK0tMezb1ug91fNpghd65hFyM4pqJAZCCCHVg5I80mBp8DWwrP0yjNWVYK/sMabs6ICC/OTaDusFVmIROhzaBf0F30JSmIJRG+dg98zvMXDNRQTerZnxcgI1PmxdjQAAeZnFuPxPDP5acAVB++6jOF9e7dcnpKYkJiaGUyseaSwoySMNGo/x8Fn/vVho1Q1XuAKM2vs+UpJDX39iDWOMwWrIR3A9cRT6HdvjA91CJOXK8PGWYPRY+R8O3kqEQlkziz3rGmli2LetIfaxQOiZR9g27zJuHH8IRYmyRq5P6jSVSqVitR0EIaTUk3+PL/3jQGPySKNx8fpqzIxYB20Aa3r9BRczz9oOqUIcx4GTy6HkC3Di8EXc3/MP1li3h7mJLj5t7wB/b1toqtXM+l0Zifm48k8MkmJyMHyxLzS1aYu0+qqKxuQdsrCwkJqamubweLz69ceDkAZGpVKxtLQ0/eTk5Nuenp59KypDSR5pVO7d/xeTr3+HPE6JFR1XwM+69tfSe5X0deuR9uuvkNvYY0ProTikMIaJjjrGtG2C4a3toS+smaSrIEcGbX0NcCoOpzZFwrmlOUQeJmCMGnXqi6pI8m7cuGEmEAg2AGgK6gkipLapAEQoFIqxLVq0SK2oACV5pNFJKUjBlMApuJ8ZhTmWneH/3i+1HdIr5V+4gKQFC6FIToas70D81qQbzjzMha6GAMNa2+NjPxHMdGtmRmxeZjEO/hqCnNQiWDrpo01/J1g41N42cuTNVUWSRwipXyjJI41SgbwAX+x6D/+pcjFGxxmff7gHPH7dXTZSmV+AtJ9/RtaOHTCdMQPJvT7C2vMxOBaeBAGfB/8WNhjf3gH2xtrVH4tShTsXk3DtyAMU5ZbAoZkpOg53gVCH9sStyyjJI6TxoSSPNFoKeTGW7v8Qe2SJeI9vgO8HHoGGZt1ulSoMCYGmVAqehgaK797FY3V9/HkrHftvJEChUuF9NwuM8LWHr4NxtXellhQrEHrmER6EpmPA7Bbg83lQKVXg8akXry6iJI+QxoeSPNKocSoVth4djxUZV9GMU8OqAUdgqGtV22G9FqdUIqZnT6gKC2Exfz6KW7XDposPset6PLIL5XA208HINiL097KGtkb1tlByKg6MxyCXKbFryTWIW5rDq5sd1IV1t2W0MaIkj5DGh5I8QgCc/G8JvondC3NdG2yQfgrLtPuAqB1g61Pbob1UUWQkkubOg+zOHei+/z4s5s6BwsAIh0IfI+DyQ0Qklo7bG9DCBiN87eFoqlOt8RTmliBoTxTuB6dCU0cN3j1FaNreGnwBtezVBZTkEdL4UJJHyBOhaaFYd+Rj/BIfAw0wML4GMOpQnU70OLkcGZs2I/3338GEQoh27YRGkybgOA4hj7IRcOkh/g1PglzJoZ2zCUb6itBZYgY+r/q6clPjcnHpQAwS72VBz0QTA770hpYejderbZTkEdL4UJJHSDnZp+dDN2gl+ABUjAde57lAu5m1HdZryWJjkb17D8xmfwnG44ErKQFTL02s0vJk2H09HtuvxCM5txjWBkIMb22PQS1tYaRdPckXx3F4dDsTsaHp6DBEDMYYctOLoGcirJbrkdejJI+QxoeSPELKe3QN3NY+UCmKUcIYTnWYgj4dFter9eDkKal4OGgQjMeNheGQIWC80u5ShVKFU7dTsPXyQ1yJzYS6gIe+nlYY6WsPDxuDao0pP0uG7fMvw8pJH779nGBqp1ut1yMvoiSPkMaHkjxCnvfoGhSx57A+NxLrMm9igPMAzGk9B2q8+rHbgzw5GUlz5qLg4kUIW7SA5eLF0HBo8kyZqJQ8BFx+iAM3E1FYokQzWwOMamOPnu6W0BBU/W4aSrkK4ecTEHzsIWQFCoh9zNGqrwO17NUgSvIIaXwoySPkJTiOw+pbq/FH2B9oY9UGKzqsgI569U5eqCocxyHnn4NIWbYMXFERTKZMgfG4sS+0SOYWy7H/RgK2XY5DbHoBjLXVMdjHFsNa2cPKoOoTMFmRAjdPxCHszCNwAEZ+14bG69UQSvIIaXwoySPkNf6+/zcWXV6EJgZNsKbLGlhoW9R2SG9MkZaG5MVLwDQ1YP3jjy8tp1JxuBiTjq2X4hB4NwUA8J7UAiPbVM+ae/lZMjy6kwnXNpYAgNiQNNhKjaCmUTN78jZGlOQR0vhQkkfIG7j8+DJmnJsBLYEWVndZDVdj19oO6a2UTcQovncPuf8ehcnkSeBpaFRY9lFmIf66Go/d1+ORVbbmnq89+jW3gU41rLmXnVqIv+ZfgZa+Onx6N4FrG0taULkaUJJHSONDSR4hb+h+1n1MOjMJObIcLO+wHO1t2td2SG8t/Y8/kfbzz1Bv0gSW3y2BVvPmLy1bLFficOhjBFyOQ3hiDnQ0BBjQ3BojfEVwMqvabuvH0dm4fCAaybG5MLTQQusPHdHE06ReTXip6yjJI6TxoSSPkLeQVpiGyWcm417WPXzj8w0GSQbVdkhvLT/oIpLnz4c8KQmGQ4fCdPp08HVevuctx3G49SgbAZfj8G9YEkqUKvg5mWCkrz26uJpX2Zp7HMfhQWg6rvwTg7wsGUZ+50v74VYhSvIIaXwoySPkLRXKC/HlhS9xPuE8RruNxvQW08Fj9at7UVVQgNRfVyJr+3aYfv45TD4d/0bnpefLsPv6I2y/EoeknNI194a1tsMgb1sY61Tc/fvWsSlVyEgsgKmdLjiOw5V/YuHS2gJGli9PRMnrUZJHSONDSR4hlaBUKfHD9R+w8+5OdLPvhu/9voemQLO2w3prRaGh0JBIwNPQQPGdO1CztATfwOC15ymUKpy+k4KAy3G4FJMBdT4P3dzMMbCFDdo7m1ZZ6152SiH2Lr0OuUwJ1zaWaNnbATqGVZNMNjaU5BHS+FCSR0glcRyH7Xe246frP8Hd1B2/df4NRppGtR1WpXBKJWJ69oSqoBAW8+ZB7/333vjc+yl5+OtqPP65lYjsQjnM9TTQz8sGA1vYVMnYvaL8Etw4Gofw8wlgPAbPzrbw7imimbhviZI8QhofSvIIeUdn4s7gq/++gonQBGu6rkET/SavP6kOKr5zB0lz5qL49m3odusG83lzoWZm9sbnyxRKBN5Jxb4bCTgXlQalioOXnQH8W9iit6cl9DTfbTHp3PQiXDkYi6SYbAz7tjUEapTkvQ1K8ghpfCjJI6QKhKWFYWrgVChUCqzstBLeFvXzbymnUCBj82ak/7YaTFMTol27Xtgt402k5hXjn5BE7A1OwP3UfGgIeHjfzQL+3jZo42jyTt25JcUKqGsKoJSrcGjVLUjbWkLsYwFWRV3EDRUleYQ0PpTkEVJFHuU9wuQzk5GQl4DFbRejl0Ov2g6p0mQPHiB79x6YffkFGI8HVUkJeOpvP9OV4ziEJ+Zgb3ACDt5KRG6xApb6mhjQ3AYDWtigiUnlJ1Pkphfh2PpwpD/Kh7GNDtr0c4St1IiWXXkJSvIIaXwoySOkCuXIcjD93HRcT76OKc2mYLzH+HqfdMhTUvDQ/yMYjx0Lw2FDwfiV6yYtlitx+k4K9t1IwIWoNKg4wNveEP7eNujpbgndSnTncioO92+k4OrBWOSmF8NGYoj3xzWFpnb92Ge4JlGSR0jjQ0keIVWsRFmCBZcW4EjsEfRz6od5vvOgxqu/SYc8JRVJ8+ai4MJ/EDZrBsvvlkDD0fGd6kzJLcaBm4nYd+MRYtIKoKnGQ4+mlvBvYYPWDsbgvWXXq1KhQuR/iYi/nYlekzzAGHvarUtKUZJHSONDSR4h1YDjOKwJXYN1oevQ2rI1fu74M3TVdWs7rErjOA65hw8j5bvvoSoshMmkiTCeMOGdWynLFlreeyMBh0MfI69YAWsDIQa0sMHA5jawM9aqVL1FeSXYsfAqnH3M4d1DBC09WlSZkjxCGh9K8gipRgejD2LhpYUQ6Yvwe5ffYaVjVdshvRNFRgaSlywBT10dVj/8UKV1F8uVOBGZjH03EhAUnQ6OA3yaGMG/RWl3rvZb7JtblF+CqwdjcftiEgRqPHi9ZwfPLraNumWPkjxCGh9K8gipZleTrmL62enQEGhgdZfVcDN2q+2Q3hknl4OpqaH43j3kHjkCk8mTwdOsusWgH2cX4e+QROy7kYAH6QXQUueXdud628BHZPTG3blZyQW48k8sYm+lQUtfHUPmt2q04/UoySOk8aEkj5AaEJMdg0mnJyFLloUf2/+IjrYdazukKpGxYQNSl6+Aur09LJcshlbLllVaP8dxuBGXhX03EnAkLAn5MgVsjYQY2NwW/Ztbw9bozbpzk2Nz8OhOJlr2Kl0OJi0+Dya2OvV+UszboCSPkMaHkjxCakh6UTqmnJmCO5l38JXPVxgiGVLbIVWJgsuXkTRvPuQJCTAcOgSmM2aCr1P1+8wWliiedudejM4AAPg6GGNACxu852b+xostZyUXYOe3V2Em0kOb/o6wcjas8ljrIkryCGl8qjXJY4x1B7ASAB/ABo7jlr2k3AAA+wC05DjulRkcJXmkPiuUF+Kr/77C2UdnMUI6AjNbzASfV/93blAVFiJt5UpkBmyD6fTpMBk/rlqvl5BV+GR2bgLiMwuhzuehvdgUfTwt0cXVHDqvGL+nUqpw90oyrh2KRUFOCUTuxmj9oSOMrd99C7a6jJI8QhqfakvyGGN8AFEAugFIAHAdwBCO424/V04XwL8A1AFMoSSPNHRKlRLLg5dj+53t6GLXBUvbLYVQIKztsKpEUXg4NFxcwFNXR/Ht2xBYWkJgWH0tZRzHIeRRNo6EJuFoeBKSc4uhLuChk4spenlYoYvE7KUTNuQlSoQFPsLNE/FQqTiMXtoGGloNd7weJXmEND7VmeT5AljIcdz7T15/DQAcxy19rtyvAE4B+ALALErySGOx/fZ2/Hj9RzQ1aYrfOv8GY6FxbYdUZTilEjE9e0KVlw+LeXOh2717tY9/U6k43IzPwpGw0oQvNU8GTTUeOkvM0MvdCp0kptBSfzHhK86XIzk2ByIPE3Ach8j/HsOphVmDm6BBSR4hjU91JnkDAXTnOG7sk9cjALTiOG5KuTLNAczhOG4AY+wcXpLkMcbGAxgPAHZ2di3i4uKqJWZCalpgfCA2RmzE+q7roaPesLoLi+9FIWnOHBRHRECnSxdYzJ8PNXOzGrm2UsUh+GEm/g1PwtHwZKTnyyBU46OLqxl6e1iio4sZNNVe7CbPSMzHriXXoCEUoEV3Edw7WUNQQbn6iJI8QhqfWkvyGGM8AIEARnMc9/BVSV551JJHGhqO4xrsLE9OoUBmwDakrVwJpq4O0e7d0HBoUqMxKFUcrj7IwL9hSTgekYyMghJoq/PRxdUcvT0s0V5s+kzCl56Qh8t/xyA+MhM6Rhpo1dcBYh+Lt96Fo66hJI+QxqfWumsZY/oAYgDkPznFAkAmgL6vSvQoySOk/imJi0PWnj0wmzULjDGoZDLwNDRqPA6FUoUrsZn4N/wxjkckI6tQDh0NAbpJSxM+P2cTaAhKE76Eu5m4dCAG+VnFGL7Yt94vpExJHiGNT3UmeQKUTrzoAiARpRMvhnIcF/mS8udALXmENHjylBQ8HOgP47GfwHD4cDB+7XSHypUqXI7JwJGwxzgRmYKcIjl0NQV4T2qB3h6WaOtkAjUeQ25GMfRNhVApVTi7/S6atreBeRO9Won5XVCSR0jjU20fTTmOUzDGpgA4gdIlVDZxHBfJGFsEIJjjuEPVdW1CSF3GoCmVImXpMuQcPQqrJUug4exc41GoPVl2pb3YFEs+VOFiTDqOhCbh5O1k7L+ZAH2hGt53M0cvDyu0MdJAXkoR4iIycPdyMhybm8GxuSly04tgLTaEhYN+jcdPCCGvQ4shE0JqHMdxyD3yL1K++w7KggKYTPgUJpMm1YmxiTKFEkH30/FvWBJO3k5BvkwBAy01dHezQA+JOdSj8xF6+hGUchXAAIGAhw+me9X5RI9a8ghpfOr3IBNCSL3EGIN+n97QbtsGKd99D3n8ozqR4AGAhqB0UkYXV3MUy5W4EJWGf8OTcDj0MXZdfwQjbXX42+hC70EhwAFKpQqJUVl1PskjhDQ+lOQRQmqNwMgI1iuWg5PLAQDFd+8i59BhmE6dAp6w9heI1lTj4z03C7znZoFiuRLn7qXhSNhjnAxLwwcQgA8GdT4P1uLGsTUaIaR+oSSPEFLrmFrpwsMFFy8ic9Mm5J0+DcvFi6HdyqeWI/s/TTU+uje1QPemFigaqMSx8w9RlFCILh3sqBWPEFIn0Zg8QkidUnDlKpLmz4c8Ph4GH30Esy9mga+rW9th1Xs0Jo+QxodX2wEQQkh52q1bweHgPzD6+GNk79uHrB07azskQgipl6i7lhBS5/CEQph/+QX0evV8urxKUUQk1KwsITAyquXoCCGkfqCWPEJInSV0cwNPXR2cUonEmTMQ26s3co78i/o2zIQQQmoDJXmEkDqP8fmw+e03qNna4vGsWUiYOAny5OTaDosQQuo0SvIIIfWCplgM0c4dMPtqNgquXEFs7z6QxcbWdliEEFJn0Zg8Qki9wfh8GI8eDd3OnZG9dy/UmzQBAKiKi8HT1Kzl6AghpG6hljxCSL2jbmcHs5kzwRiDPDkZ0d26IWPjJqhKSmo7NEIIqTMoySOE1GuMz4fQwxORa1djy/ABSLxwrrZDIoSQOoGSPEJIvSYwNYXN6t9gOG4sClRK7F79E07N+gxyatUjhDRylOQRQuo9xhjcP52E0b+sg52OAcIexeCvrz9H0v17tR0aIYTUGtrWjBDS4MQEX8HpTeuQn5kBiYkVui5eCg3Dxr2IMm1rRkjjQy15hJAGx9G7NUYvXwMXmya4m5aITeNG4N4e2h6NENK4UJJHCGmQNLS00Hv5KnwwaDQYgCP7/8I/40ahIDmptkMjhJAaQUkeIaRBc+o/EB9v2gE3C1vE5qRj6xdTEH39Sm2HRQgh1Y7G5BFCGo340ydx7sQhpMU/hJPUA51GjoVeE4faDqtG0Jg8QhofSvIIIY2KUqHA9YN7cXn3dvBVHNq064LmU6eDx2vYHRuU5BHS+DTs32qEEPIcvkCA1gOGYPD0b6DLF+D8pbPYOcIfmZHhtR0aIYRUKUryCCGNkqVvW4zcvg+tpF5IKylCwMLZCN6xFZxKVduhEUJIlaDuWkJIo5cWegsnf/8ZyTmZsJa4oduYT2Esalhj9ai7lpDGh5I8QggBwHEcIs+dxrmtf0JeWAAvcVO0m7cIfA2N2g6tSlCSR0jjQ921hBCC0q3RmnbqhhHzvoelpjZu3I/ElhEDkXAusLZDI4SQSqEkjxBCytF3dMKgrbvRtVMPFKmU2LNmBU7MnAq5TFbboRFCyFuhJI8QQp7DGIPnhMkY/et6iHQMEZHwANu/+gyJ9+7UdmiEEPLGaEweIYS8xoOb13Fq4xrkpafDxcQCXRctg6axcW2H9VZoTB4hjQ8leYQQ8gZKigpxav7XuBsfAy2lCh192sNU3whaPi2h5eVV2+G9FiV5hDQ+lOQRQshbiPlnP67+sQYeD5LAAweeugbstmyu84keJXmEND40Jo8QQt6C44cD0KX7h+BxHBgHKGUyPNq/r7bDIoSQF1CSRwghb0m3nR94mpoAjweOx0PQras49PP3KMjOqu3QCCHkKUFtB0AIIfWNlpcX7DZvQuG169Bs0RwFj2Jxef9OPIoIQ8dR4yBt3xmMsdoOkxDSyNGYPEIIqQIZiY9wcv1veHzvNkSezdFt3BTomZrVdlhP0Zg8Qhof6q4lhJAqYGxti8ELl6HzxxOQeO8OtsychJDjh8GpVLUdGiGkkaIkjxBCqgjj8eD1fm+MXv47rCVSBG5ej10Lv0JG4qPaDo0Q0ghRkkcIIVVMz9QM/b/+Ft0nTUdmQjy2fTkVV//eA6VCUduhEUIaEUryCCGkGjDG4NahC0b/vBaOLVohaFcA/pozAykPYmo7NEJII0FJHiGEVCNtA0P0mfE1+s74BoXZWfjrm+n4b8cWKEpKajs0QkgDV61JHmOsO2PsHmMsmjH2VQXHZzDGbjPGwhhjZxhj9tUZDyGE1BbnVm0wesVaSNt3xrWD+xDw5VQk3I2s7bAIIQ1YtSV5jDE+gN8B9AAgBTCEMSZ9rlgIAG+O4zwA7APwY3XFQwghtU1TRwfdJ36OAXMWQ6lQYPeC2TizaS1KigprOzRCSANUnS15PgCiOY6L5TiuBMAuAB+UL8Bx3FmO48p+u10BYFON8RBCSJ0g8vDCqOWr0bxHX9w6eRRbZk3Gg1s3ajssQkgDU51JnjWA8usGJDx572U+AXCsogOMsfGMsWDGWHBaWloVhkgIIbVDXVOITqPHY8iiH6GmoYkDSxfg2O8/oygvt7ZDI4Q0EHVi4gVjbDgAbwA/VXSc47g/OI7z5jjO29TUtGaDI4SQamQldsWIH1ahdf9BuHvxPLbMnISoK0Gob7sREULqnupM8hIB2JZ7bfPkvWcwxroCmAOgL8dxsmqMhxBC6iSBmhraDhqBYd//Al1jExz+ZRkOrfge+VmZtR0aIaQeq84k7zoAZ8ZYE8aYOoDBAA6VL8AY8wKwHqUJXmo1xkIIIXWemcgBQ5esQLuho/Hw1g1smTkREWdPUaseIaRSqi3J4zhOAWAKgBMA7gDYw3FcJGNsEWOs75NiPwHQAbCXMXaLMXboJdURQkijwOPz4fPBQIz48TeY2jXBiXUrse+7echJTa7t0Agh9Qyrb58Qvb29ueDg4NoOgxBCqh2nUiHszHFc+GszVCoV3p/wGSRt2leqLsbYDY7jvKs4REJIHSao7QAIIYRUjPF48OzWE028WiJw8zoYWr5qgQJCCHkWJXmEEFLH6ZmY4sMv5tV2GISQeqZOLKFCCCGEEEKqFiV5hBBCCCENECV5hBBCCCENECV5hBBCCCENECV5hBBCCCENECV5hBBCCCENECV5hBBCCCENECV5hBBCCCENUL3b1owxlgYgrpKnmwBIr8JwGgJ6JhWj5/IieiYvqk/PxJ7jONPaDoIQUnPqXZL3LhhjwbR347PomVSMnsuL6Jm8iJ4JIaQuo+5aQgghhJAGiJI8QgghhJAGqLEleX/UdgB1ED2TitFzeRE9kxfRMyGE1FmNakweIYQQQkhj0dha8gghhBBCGgVK8gghhBBCGqAGk+Qxxrozxu4xxqIZY19VcFyDMbb7yfGrjDFRuWNfP3n/HmPs/RoNvBpV9pkwxkSMsSLG2K0nX+tqPPhq8gbPpD1j7CZjTMEYG/jcsVGMsftPvkbVXNTV6x2fibLcz8mhmou6+r3Bc5nBGLvNGAtjjJ1hjNmXO9Ygf1YIIfUMx3H1/gsAH0AMAAcA6gBCAUifKzMJwLon3w8GsPvJ99In5TUANHlSD7+276mWn4kIQERt30MtPRMRAA8AAQAGlnvfCEDsk/8aPvnesLbvqTafyZNj+bV9D7X4XDoB0Hry/cRy/34a5M8KfdEXfdW/r4bSkucDIJrjuFiO40oA7ALwwXNlPgCw9cn3+wB0YYyxJ+/v4jhOxnHcAwDRT+qr797lmTRUr30mHMc95DguDIDquXPfB3CK47hMjuOyAJwC0L0mgq5m7/JMGrI3eS5nOY4rfPLyCgCbJ9831J8VQkg901CSPGsAj8q9TnjyXoVlOI5TAMgBYPyG59ZH7/JMAKAJYyyEMXaeMdauuoOtIe/y/7ox/5y8iiZjLJgxdoUx9mGVRla73va5fALgWCXPJYSQaiGo7QBInZQEwI7juAzGWAsA/zDG3DiOy63twEidY89xXCJjzAFAIGMsnOO4mNoOqiYxxoYD8AbQobZjIYSQ8hpKS14iANtyr22evFdhGcaYAIA+gIw3PLc+qvQzedJ1nQEAHMfdQOnYJHG1R1z93uX/dWP+OXkpjuMSn/w3FsA5AF5VGVwteqPnwhjrCmAOgL4cx8ne5lxCCKluDSXJuw7AmTHWhDGmjtJJBM/P9DsEoGyW20AAgRzHcU/eH/xkpmkTAM4ArtVQ3NWp0s+EMWbKGOMDwJMWGmeUDh6v797kmbzMCQDvMcYMGWOGAN578l59V+ln8uRZaDz53gRAWwC3qy3SmvXa58IY8wKwHqUJXmq5Qw31Z4UQUs80iO5ajuMUjLEpKP1FygewieO4SMbYIgDBHMcdArARwDbGWDSATJT+0saTcntQ+sdJAWAyx3HKWrmRKvQuzwRAewCLGGNylA62n8BxXGbN30XVepNnwhhrCeBvlM6K7MMY+5bjODeO4zIZY4tR+scfABY19mcCwBXAesaYCqUfGJdxHNcgkrw3/PfzEwAdAHufzFeK5ziub0P9WSGE1D+0rRkhhBBCSAPUULprCSGEEEJIOZTkEUIIIYQ0QJTkEUIIIYQ0QJTkEUIIIYQ0QJTkEUIIIYQ0QJTkEfISjDFjxtitJ1/JjLHEJ9/nM8bW1HZ8hBBCyKvQEiqEvAHG2EIA+RzHLa/tWAghhJA3QS15hLwlxlhHxtiRJ98vZIxtZYz9xxiLY4z1Z4z9yBgLZ4wdZ4ypPSnXgjF2njF2gzF2gjFmWbt3QQghpKGjJI+Qd+cIoDOAvgC2AzjLcZw7gCIAvZ4ker8BGMhxXAsAmwB8V1vBEkIIaRwaxLZmhNSyYxzHyRlj4SjdAuv4k/fDAYgAuABoCuDUk+2v+ACSaiFOQgghjQgleYS8OxkAcBynYozJuf8PdFWh9N8YAxDJcZxvbQVICCGk8aHuWkKq3z0ApowxXwBgjKkxxtxqOSZCCCENHCV5hFQzjuNKAAwE8ANjLBTALQBtajUoQgghDR4toUIIIYQQ0gBRSx4hhBBCSANESR4hhBBCSANESR4hhBBCSANESR4hhBBCSANESR4hhBBCSANESR4hhBBCSANESR4hhBBCSAP0PyRmLEeI5I8fAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fn = lambda t, y, params: -params[0]*y\n",
    "params = (5,)\n",
    "y0 = 1\n",
    "t0 = 0\n",
    "h = 0.2\n",
    "tru_t = np.linspace( t0, t0+h, 10 )\n",
    "tru_y = y0*np.exp(-params[0]*tru_t)\n",
    "\n",
    "f1 = fn(t0,y0,params);\n",
    "k1_t = t0+h/2\n",
    "k1_y = y0+0.5*h*f1\n",
    "f2 = fn(k1_t,k1_y,params);\n",
    "k2_t = t0+h/2\n",
    "k2_y = y0+0.5*h*f2\n",
    "f3 = fn(k2_t,k2_y,params);\n",
    "k3_t = t0+h\n",
    "k3_y = y0+h*f3\n",
    "f4 = fn(k3_t,k3_y,params);\n",
    "kf_y = y0 + h/6*(f1+2*f2+2*f3+f4);\n",
    "kf_t = t0 + h\n",
    "eps = h/10;\n",
    "\n",
    "plt.plot( tru_t, tru_y.T, 'C0' );\n",
    "plt.plot( [t0,k1_t], [y0,k1_y], 'C1--' );\n",
    "plt.plot( [k1_t-eps,k1_t+eps], [k1_y-eps*s1,k1_y+eps*s1], 'C2')\n",
    "plt.plot( [t0,k2_t], [y0,k2_y], 'C2--' );\n",
    "plt.plot( [k2_t-eps,k2_t+eps], [k2_y-eps*s2,k2_y+eps*s2], 'C3')\n",
    "plt.plot( [t0,k3_t], [y0,k3_y], 'C3--' );\n",
    "plt.plot( [k3_t-eps,k3_t+eps], [k3_y-eps*s3,k3_y+eps*s3], 'C5')\n",
    "plt.plot( [t0,kf_t], [y0,kf_y], 'C4--' );\n",
    "plt.plot( k1_t, k1_y, 'C1.')\n",
    "plt.plot( k2_t, k2_y, 'C2.')\n",
    "plt.plot( k3_t, k3_y, 'C3.')\n",
    "plt.plot( kf_t, kf_y, 'C4.')\n",
    "plt.legend( bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0., \n",
    "           labels=['True Y', 'Step 1','Slope for Step 2','Step 2','Slope for Step 3','Step 3','Slope to be used with others for Step 4','Step 4']);\n",
    "plt.xlabel( 'Time')\n",
    "plt.ylabel( 'Y');"
   ]
  }
 ],
 "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
}