{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adaptive Step Size Method\n",
    "\n",
    "Our goal is to implement the adaptive step size forward Euler method with the explicit trapezoidal method embedded in it for estimating the error. More explanation is in the lecture notes PDF. Here is the algorithm.\n",
    "\n",
    "<ul>\n",
    "<li> Initialize $h$ to a small number (e.g 1/100th of the time span (t0 to tF) (Note 1/100th is chosen arbitrarily - you may want to start smaller)\n",
    "<li> Initialize lists for y, t, and the step-sizes.\n",
    "<li> While the final time in the list of time steps is $<$ tF\n",
    "\t<ul>\n",
    "\t<li> Compute YFE using step size $h$ with the forward Euler method (the beginning of the step is the last element of the time list and the last element of the y list)\n",
    "\t<li> Compute YET using step size $h$ with the explicit Trapezoidal method\n",
    "\t<li> Estimate the error (err) using the Euclidean distance between YET and YFE.\n",
    "\t<li> If err $\\leq$ ETOL, then append the time step to the list of time steps, append YFE to the list of solutions, and append $h$ to the list of step sizes.\n",
    "\t<li> Determine the step size for the next step.\n",
    "\t\t\\begin{equation*}\n",
    "\t\th_{new} = h \\cdot \\sqrt{0.9 \\frac{ETOL}{err}}\n",
    "\t\t\\end{equation*}\n",
    "    <li> If $h_{new} > 2h$ then $h_{new} = 2*h$ (make sure h doesn't more than double).\n",
    "\t<li> Set $h$ to $h_{new}$. This will be used either for the next step (if err was small enough) or to retry the current step (but with adjusted size) again.\n",
    "    </ul>\n",
    "</ul>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "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": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "IndentationError",
     "evalue": "unexpected indent (<ipython-input-1-b8f2f26870b6>, line 26)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-1-b8f2f26870b6>\"\u001b[0;36m, line \u001b[0;32m26\u001b[0m\n\u001b[0;31m    pass\u001b[0m\n\u001b[0m    ^\u001b[0m\n\u001b[0;31mIndentationError\u001b[0m\u001b[0;31m:\u001b[0m unexpected indent\n"
     ]
    }
   ],
   "source": [
    "def forwardEulerAdaptiveStepSize(f,t0,tFinal,y0,params,ETOL):\n",
    "    \"\"\"\n",
    "                  Runge-Kutta Method of Order 1 (Forward Euler) embedded\n",
    "                       in Runge-Kutta Method of Order 2 (Trapezoidal)\n",
    "    Performs Forward Euler method of integration\n",
    "    for O.D.E's, using an error estimate (assuming an embedded order 2 method is \n",
    "    the true solution) to control the step-size.\n",
    "    [t,y,h] = forwardEulerAdaptiveStepSize(f,h,t0,tFinal,y0)\n",
    "      f : function handle for f(t,y) = y'\n",
    "          (where f has signature y = f(t,y,params) with\n",
    "               t as a scalar\n",
    "               y is a column vector with the same length as y0\n",
    "               params is a vector with the parameters)\n",
    "      t0     : initial time\n",
    "      tFinal : final time\n",
    "      y0     : initial conditions\n",
    "      params : model parameters\n",
    "      ETOL   : maximum error tolerated. The step size will be chosen to guarantee \n",
    "               that the estimated error is smaller than ETOL.\n",
    "    Return values:\n",
    "      t      : the entire series of time steps\n",
    "      y      : the series of values of y at each time step\n",
    "               (has dimensions length(y0) by length(t))\n",
    "      steps  : the set of step-sizes taken internally\n",
    "      \"\"\"\n",
    "      pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run and plot result from adaptive method.\n",
    "\n",
    "Run both the adaptive forward Euler method on the degradation model with degradation rates (2,5) and initial conditions (1,1) with time beginning at 0 and ending at 4. Use an error control (EST) of 1e-4.\n",
    "\n",
    "Then plot the true solution (solid lines) and the adaptive forward Euler's solution (dashed lines).\n",
    "\n",
    "Then plot the difference over time between the true solution and adaptive forward Euler's solution (dashed lines).\n",
    "\n",
    "Finally, plot the step size over time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Write code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step Size vs Desired Estimated Error\n",
    "\n",
    "Determine which time steps are chosen, given the value of ETOL. Also, compare the ETOL to the actual error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "# range from smallest timestep that doesn't take way too long to run to largest time step with qualitatively correct answer.\n",
    "ETOLs = np.array([1e-5,1e-4,1e-3,1e-2])\n",
    "act_err = np.zeros( (ETOLs.size,) )\n",
    "params = (2,5)\n",
    "y0 = np.array((1,1))\n",
    "t0 = 0\n",
    "tF = 4\n",
    "all_steps = []\n",
    "\n",
    "# Write code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bonus Code\n",
    "### Examine step size changes over time\n",
    "\n",
    "Do we take lots of rejected steps?\n",
    "\n",
    "Examine each step size attempted. To do this, I have taken the code out of the function and added lists to keep even the rejected steps. Here, I use the degradation model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa8AAAEYCAYAAADrpHnMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAA43UlEQVR4nO3deXhd1XX///eSZMmWR1ke5UmesY0HbGOMDQkzBoJNEhIMJKEhhKZp0qb9pi1JmpbybfJr8s3TZqJpKRmAQAghoRYEcDAzBoNtkOcBWZ4keZBkebY13fX74xzZF6HhXulKV1f+vJ7nPrr3nH2O1kZCy/ucdfY2d0dERCSVpCU7ABERkXgpeYmISMpR8hIRkZSj5CUiIilHyUtERFJORrID6GoGDRrk+fn5yQ5DROSct3bt2gp3H9zUPiWvRvLz81mzZk2ywxAROeeZ2e7m9nWpy4ZmtsjMtplZkZnd08T+LDP7bbj/bTPLD7fnmtnLZnbczH7a6Jg5ZrYhPObHZmad1B0REekgXSZ5mVk6cD9wHTAVuNXMpjZq9gWgyt0nAP8BfC/cfhr4NvD1Jk79M+CLwMTwtSjx0Qe2rl7BWw99k62rV3TUtxAREbpQ8gLmAUXuXuzuNcDjwJJGbZYAD4XvnwSuNDNz9xPu/gZBEjvDzIYD/dx9lQdTiTwM3NQRwW9dvYL8Z5ZyUfH9jHvmFjau/GNHfBsREaFrJa8RwN6ozyXhtibbuHsdcATIbeWcJa2cEzO728zWmNma8vLyNoQOVZtfogd1pBlkWh0T/nQH7/zwVja8toz6uro2nVNERJrWlZJX0rj7A+4+193nDh7cZGFLq3KmXkENPajzNKo9g229ZjG16mWmv/Q5qv51PKvu/yLb1ryERyIJjl5E5NzTlaoNS4FRUZ9HhtuaalNiZhlAf6CylXOObOWcCXHehVexld9QtfklcqZewcwLr+L0yeO8+9qTsOFJLjj4FFnPPEHpH4eyJ+8G8i75DGOmzOmIUEREuj3rKrPKh8loO3AlQYJZDdzm7pui2vwlMN3dv2RmS4FPuPuno/b/GTDX3b8Ste0d4K+At4FngZ+4+7PNxTF37lzviFL5o4cr2fryY/Ta+gemnn6PdHN2pI/l4Jgbyf/oZxk+ZlLCv6eISCozs7XuPrfJfV0leQGY2fXAD4F04Bfu/h0zuw9Y4+4FZtYTeAS4ADgELHX34vDYXUA/IBM4DFzj7pvNbC7wK6AX8BzwVW+h0x2VvKJV7N9L0cuPkLNjGZPrtgKwpcc0jk64iYmX387AIR+6LScics5JmeTVFXRG8opWWryFva89wrA9T5Mf2UOdp7Gp1xxqpnySKZcvpU+/nE6LRUSkK1HyikNnJ69oOze9zf6Vvya/7DmGU84pz2Rz3wXYjE8x7SOfIKtndlLiEhFJBiWvOCQzeTXwSIRta17kyDuPMbFiBQM5ylF6szXnMrJn38KUi28gPaMr1dqIiCSeklccukLyilZXW8PmlU9T/d5vmXL4NfrYKcrJYceQaxg4/zYmzvoIlqYnHkSk+1HyikNXS17RTp88zuZXf4dtfJJpx1eRaXWU2DD2jghL78+bnewQRUQSRskrDl05eUU7UlXBtlcaSu8Lw9L7cRzMv5Gxl32OYaMmJDtEEZF2UfKKQ6okr2gV+/dEld5vA2Bzj/M5OvEmJl/+GXIGD09yhCIi8VPyikMqJq9opcWb2PPaIwzf8wz5kb3Uejqbs6NK7/sOSHaIIiIxUfKKQ6onrwYeibBz82oOrPw1Y/c9x7CG0vt+C0mb8SmmXvpxld6LSJem5BWH7pK8okXq69m+ZgVHVj/OpIoV5Jwpvb+c7DlLmTL/OpXei0iXo+QVh+6YvKLV1lSzZWUB1YVPMPXwa/S202Hp/bVh6f2lKr0XkS5BySsO3T15RTt14hibX32StI2/Y9qJt8PS++FB6f2ln2XM5FnJDlFEzmFKXnE4l5JXtCNVFWx/+VF6bfsDU0+vI82covTxlOffyNjLPqvSexHpdEpecThXk1e0irLdFL3yCAOLlzGpbjsAmzOnc2zCTUy64jPkDBqW5AhF5Fyg5BUHJa8PKinayN7XHyFv7zOMiZSEpfdzqZ36SaZcdgu9VXovIh1EySsOSl5N80iE4o1vc/CtXzN23/MMo4KTnsWWfgtJm/lppl36cTKzeiY7TBHpRpS84qDk1bpIfT3b1rzA0Xd+w6TKF8nhGEfozbaBVwSl9xctUum9iLSbklcclLziU1tTzeaVy6h97wmmHnmNbKvmIAMpHnotuRffzoQZC1V6LyJtouQVByWvtgtK739L2sbfh6X39ey1PEpG3sCISz/D6EmzWLt6Jas2bmf++ZOYc+HCZIcsIl2YklcclLwS48ihcra9/CjZ2586U3r/vC/gq9V/Tj1pZFLHo58cqgQmIs1qKXnpxoR0iP4DBzPvk18DvkZ52S52vPIIazYfoJYMwDhNGv/8/G7+vv8kLpkwiLQ0S3LEIpJKNPJqRCOvjrN29Upu+/1BasjAgAgGGJnpacweM4DbLxrDVVOG0iszPdmhikgXoJGXdAlzLlzIY4T3vKZNpHbQFB54rZg3d1SwqvgQq4oPkWZw/oj+fGrOKK6dNpQh/VR+LyIfppFXIxp5db5IxFlVXMnP39jJG0UVVNdFzuwbP7g3S2aO4KqpQ5kyvC9murwocq5QwUYclLySq7Y+wuvby/n1qt28XlRBbf3Z38/BfbO47vxhXDVlKBeNG0hWhi4vinRnSl5xUPLqOk7W1PHiloM8sWYvbxZVUO9ggAPZmelcNnkwV00ZyuWTh5DTOzPZ4YpIgil5xUHJq2s6fLKG5zfu5w/vlvLOrkMAZKQZdREnzWDumIFcNXUIV04ZyvjBfZIcrYgkgpJXHJS8ur79R07zzPoylhWWsqH0KBCMxE7W1AMwblBvrpo6lCvPG8KcMTlkpGuGD5FUpOQVByWv1LKz4gQFhWUsW1dKcfkJ0gxysjM5fKqW+ogzILsHV0wORmQfmTSIvj17JDtkEYmRklcclLxSk7uzqewoBevKKCgsY//R02RmpDG0bxZVJ2s5Xl1Hj3Rj/rhcrpoylCunDGFkTnaywxaRFqRM8jKzRcCPgHTgQXf/t0b7s4CHgTlAJXCLu+8K930D+AJQD/yVuy8Pt/8NcBfBff4NwOfd/XRzMSh5pb5IxFm96xDL1pXx7IZ9HD5ZS5+sdEYNzObIyVrKjgQ//vOG9eXqqUO5aspQpo/or1k+RLqYlEheZpYObAeuBkqA1cCt7r45qs2XgRnu/iUzWwp83N1vMbOpwG+AeUAesAKYBAwD3gCmuvspM3sCeNbdf9VcHEpe3UtNXYQ3ispZVljGnzYd4FRtPYP7ZDJ+SB+On65j876jRDwow79qyhCuPG8oCycM0iwfIl1AqsywMQ8ocvdiADN7HFgCbI5qswS4N3z/JPBTC55aXQI87u7VwE4zKwrPt4egj73MrBbIBso6oS/SRWRmpHHFeUO54ryhnKypY8WWgxQUlvLq9nJq650xudlMGdaPU7X1PL1uH795Zy89e6RxyYTBXDVlCFdMGcKQvprlQ6Sr6UrJawSwN+pzCXBRc23cvc7MjgC54fZVjY4d4e5vmdkPCJLYKeBP7v6nxt/YzO4G7gYYPXp0YnojXU52ZgaLZ+axeGYeh0/W8NzG/SwrLGX55v24w/kj+jFjRH/qIs7KokpWbDkAwKxRA7hqyhCumjqUyUM1y4dIV9CVklfCmVkOwahsLHAY+J2Zfcbdfx3dzt0fAB6A4LJhZ8cpnW9Adia3zhvNrfNGR5Xel/HYO3sxg4vyB/KpOSOpizivF1Xwgz9t5wd/2s7InF5cNSW4TzZv7EAyM1SGL5IMXSl5lQKjoj6PDLc11abEzDKA/gSFG80dexWw093LAczsD8AC4APJS85tw/r35K5Lx3HXpeMoLj9+pmLxhy++T49046OThvCvN51PxJ3Xtpfz+Oo9/OrNXfTNyuAjkwdz9ZShXDZ5MAOyNcuHSGfpSgUbGQQFG1cSJJ7VwG3uvimqzV8C06MKNj7h7p82s2nAY5wt2HgRmAjMBX4BXEhw2fBXwBp3/0lzcahgQ+Bs6f2ywlKeXreP/UdPk52ZzjVTh7Lo/GEAvLKtnBe3HqT8WDXpacbcMTlcPXUoV04ZythBvZPcA5HUlxLVhgBmdj3wQ4JS+V+4+3fM7D6ChFNgZj2BR4ALgEPA0qgCj28BdwJ1wNfc/blw+78At4Tb3wPuCgs7mqTkJY1FIs47uw6xrDAovT9yqpac7B5cP304N87MIys9jRe3HmTFlgNs3X8MCGbDvyosw589Ood0leGLxC1lkldXoOQlLampi/D6+0Hp/Qubg9L7vP49uXFmHotn5dGvZw9eChPZquJKauudgb0zuSy8vHjppMH0yepKV+tFui4lrzgoeUmsTtbU8cLmAxQUlvHq9nLqIh6sPzZrBItn5pHbJ5PXtlewYssBXtp6kCOnaslMT2P++FyunhJMWZU3oFeyuyHSZXVK8jKzZ4EvN8x4kaqUvKQtqk6cLb1/Z9ch3GHmyP4snjWCG2cMZ2DvTNburmLFlgOs2HKQnRUnAJg6vF94eXEI5+f15/HVe3lu4z6uO384t12kxzbk3NZZyetTwHeAh4Dvu3ttQk7cyZS8pL32HTnFM+v2sWxdKRtLj2IGF4/LZcmsPBZNG07/7B7sKD/Ois0HeHHLQdbsPkTEoW9WOseq68+c57sfn64EJue0TrtsaGZ9gG8DiwgKK86s5+7u/56wb9SBlLwkkXaUH6egsIyCdWXsrDhBZnoaH508mCWz8rjyvKH0ykzn0IkaXtl2kO8+u4WK4zVnju2blcEdC/JZMCGX2aNz6NlDU1bJuaUzp4eqAU4AWUBfopKXyLlo/OA+/M3Vk/jaVRPZUHqEgsIynl4fFHv0zkznmmnDWDwrjxtn5nG6NsI3n9pw5tgB2T342as7+OnLRWRlpHFh/kAWThjEwgm5TMvrrwpGOacl8rLhIuDfgQLgPnc/mZATdzKNvKSj1Uecd3YeomBdKc9u2M+RU7UM7J3J9dOHkZ2ZwZayo1w3PbjndfR0Le8UH+KNogre3FHB9gPHAejfqwfzxzUks0GMG9Rb01ZJt9NZ97xeB74U/VBxKlLyks5UUxfhte3lLFtXxgub93O6NhKU3s/KY8nMEUwZ/sG5FA8eO81bOypZWVTByqJKSg+fAmBYv54smJDLwvFBMhvWX5MJS+pTqXwclLwkWU5U17FiywGWFZbxWlh6P2FIH5aEz5CNyf3grB3uzu7Kk6zcUcGbRZW8uaOCqpNBndT4wb1ZOGEQC8YP4uJxufTP1grSknqUvOKg5CVdwaETNTy3cR/LCst4Z+chAGaOGsCSmXl8bMZwhvT78MgqEnG27D96ZlT2zs5DnKqtJ83g/BH9g0uM4wcxN1/FH5IalLzioOQlXU3Z4VNnZr3fVHaUNIOLx+eyZOYIrj1/GP17NT2qqqmLULj3cHC/rKiCwr2HqYs4mRlpzBmdwyUTB7FgfC7TR/QnI12z40vXo+QVByUv6cqKDjbMel/KrsqTZKancdnkwSyZNYIrpwxpcUR1vLqO1TsPsbKogjeKKs7Mw9g3K4OLxuWycEIul0wYxIQhfVT8IV2CklcclLwkFbg760uOULCujKfXlXHwWDW9M9O5Niy9XzhhED1aGU1VHK/mrR3BvbKVRZXsORQUCA/pm8WC8bksCCsZR2gKK0kSJa84KHlJqqmPOG8XV1KwLpj1/ujpOnJ7Z3L99OEsmZXH7NE5pMXwTNjeQyeD+2U7KnlrR8WZB6bHDurNgvG5LJwQFH/k9Na6ZdI5lLzioOQlqay6rp7XtlewrLCUFVsOcLo2wogBvbhxZh5LZuVx3rC+MV0SdHe2HTjGG+9X8OaOSt4uruRETT1mMC2vHwvHD2LBhEFcmJ9DdqZmyZeOoeQVByUv6S6OV9exYvMBlhWW8tr7FdRHnIlD+rBkVh6LZ45gdG52zOeqrY+wvuQwK4sqeaOogvf2VFFb7/RINy4YncMl4cwfM0YOaPVypUislLzioOQl3dGhEzU8u2EfBYVlvLMrKL2fNWoAS2blccOM4QzpG99DzSdr6li9q4o3iypYuaOCTWVHcYfemelh8UeQzCYPjW2kJ9IUJa84KHlJd1d6+BTPrAtK7zfvC0rvF4wfxOJZeSw6fxj9esb/QHPViRreKg5m/nhzR+WZJV8G9cnk4vGDWBjeMxs1MPbRnoiSVxyUvORcUnTwGAWFZSxbV8buypNkZqRxeVh6f8V5LZfet6T08KkgkYUFIOXHqgEYPTCbhRNyWTA+eMYst09WIrsj3YySVxyUvORc5O6sKzk76335sWr6ZGVwzbShLJk1goXjc9v8ILO7U3TwePh8WVD8cay6DoApw/udGZXNGzuQ3lkq/pCzlLzioOQl57r6iLOquJKCwjKe3biPY2Hp/Q0zzpbet+c+Vl19hA2lR3gznGB4ze4qauoiZKQZF4wewIJwcuFZowaQmaHij3OZklcclLxEzqquq+eVbeUUrCtjxeYDVNdFGJkTXXrfr93f43RtPWt2VYUTDFewofQIEYfszHQuzB/IJRMGsWBCLlOG9YvpeTXpPpS84qDkJdK049V1vLB5P8sKy3g9LL2fPLQvi2flsXhmXsKKMY6crGXVzoZlXyrYUR4UfwzsncnF43LPLP0yJjdblYzdnJJXHJS8RFpXebyaZzfup6CwlNW7qgC4YHQw6/0NM/IY3DdxhRj7j5zmzR0V4QTDlew/ehqAEQN6sXBCcL9sT+VJ3tl1iOvODxbxlO5BySsOSl4i8SmpOskz64PlW7aEpfcLJwxi8cw8rm1j6X1z3J3iihNBFWO4htnR03UfaHPTrBF8+fLxTBjcR5cZU5ySVxyUvETa7v0DxygInyHbcygovb9i8hCWzMrj8naU3jenPuJ84j9Xsq7kyIf29e2ZwQWjc5g9egBzxuQwa9QA+iYwkUrHU/KKg5KXSPu5O4V7D4ez3u+j4ng1fbMyuGbaMJbMymNBO0rvG3vs7T1886kNZz7/7dWTGDGgF+/uqWLt7iq2HTiGO5jB5KF9uWB0DnPGBElt7KDeum/WhSl5xUHJSySx6iPOWzsqKVhXynMb93PsdB2D+mRyw/ThLJ41gtmjB7Q7gTz29h6e27ivyXtex07Xsm7vEdburuLdPcHrWHipMSe7B7NH5zB7TA6zR+cwc1R/TTTchSh5xUHJS6TjnK4NSu+fXlfGii1nS+8Xz8xjyawRTB7Wt8NjiEScHeXHzySztburzlQ0pqcZU4b3ZU5UQhuZ00ujsyRJmeRlZouAHwHpwIPu/m+N9mcBDwNzgErgFnffFe77BvAFoB74K3dfHm4fADwInA84cKe7v9VcDEpeIp3j2Ola/rTpAAXrynijKCi9P29YX26cmdjS+1gcPlnDe3sOn0lmhXsPc7KmHoDBfbPCZBbcO5uW1z/h9+6kaSmRvMwsHdgOXA2UAKuBW919c1SbLwMz3P1LZrYU+Li732JmU4HfAPOAPGAFMMnd683sIeB1d3/QzDKBbHc/3FwcSl4ina/ieDXPbQgqFtfsDkrv54zJYfHMYNb7QZ08B2JdfYRtB47x7p7DvBuO0HZXBitNZ6anMW1EP2afuXeWw7D+8c3KL7FJleR1MXCvu18bfv4GgLv/f1Ftlodt3jKzDGA/MBi4J7ptQztgM1AIjPMYO6rkJZJcJVUneXrdPpYVlrJ1/zHS04wF43NZMmsE104bmrSKwfJj1Wfumb27u4r1JUeorosAwTNnF4RVjXPG5DBleD+ta5YALSWvrnRncgSwN+pzCXBRc23cvc7MjgC54fZVjY4dAZwCyoFfmtlMYC3w1+5+okN6ICLtNjInm7+4bDx/cdl4th9omPW+lK//bh3ffCqNK88LSu8vm5z40vuWDO6bxbXThnHttGEA1NRF2LzvKO/urmJtmNCeWb8PgJ490pgxckDU6GyAZtBPsK6UvDpCBjAb+Kq7v21mPyIYpX07upGZ3Q3cDTB6tJ7OF+kqJg3ty9evncz/uWYS7+09TEFhGc+s38dzG/fTNyuDa88PSu8vHpe40vtYZWakMWvUAGaNGsCdjAVg35FTvLv78JlikJ+/Ucx/vRpc9MnPzf5AZePkYX1J10PUbdaVklcpMCrq88hwW1NtSsLLhv0JCjeaO7YEKHH3t8PtTxJeYozm7g8AD0Bw2bDdPRGRhDKz4A//6Bz+8YYpvBXOev/8xv08ubaEQX2y+NiM4SyelccFo9pfet9Ww/v34oYZvbhhxnAgqK7cWHq2TP+19yv4w3vBn7XemenMGj2AOaNzuGBMDrNH5dA/Ww9Rx6or3fPKICjYuJIg8awGbnP3TVFt/hKYHlWw8Ql3/7SZTQMe42zBxovAxLBg43XgLnffZmb3Ar3d/e+ai0P3vERSR1B6fzCY9X7LQWrqIowaeLb0ftLQji+9j4e7U1J1irW7q84ktC37jhIJ/wxPGNLnA5WN4wad21NcpUTBBoCZXQ/8kKBU/hfu/h0zuw9Y4+4FZtYTeAS4ADgELHX34vDYbwF3AnXA19z9uXD7LIJS+UygGPi8u1c1F4OSl0hqOna6luVh6f3KqNL7xbPyuHFG55bex+NEdR3rShqqGoNy/cMnawHo1zPjzGXGOWNymDlqAH3OoQU7UyZ5dQVKXiKpr/xYNc9u2EfBujLWhqX3c8fksHhWHtdP7/zS+3g0TD68dncV74XPnb1/8DjukGYweVi/M/M1zh6d062XhlHyioOSl0j3svfQSZ5eX0ZBYdmZ0vtLwlnvr0li6X08jpyqpXDv2WfOCvcc5lh1MMVVbu/MD8zXOGPkAHpldo+HqJW84qDkJdJ9bdt/jIJ1pSwrLKOk6hRZGWlcOWUIi2eO4LLJg1Nm5oz6iPP+wWMfqGzcWRE8AZSRZkzNi3qIekwOef17puToTMkrDkpeIt2fu/PunsM8va6MZ9aXUXG8hr49M1g0bRhLZo3g4vG5KVfGXnm8+gNTXK0vOcKp2mCKq2H9ejJ7zIAzpfrT8vqRldH1E7WSVxyUvETOLXX1Ed7cUUnBuqD0/nh1HYP7ZnHD9OEsmZXHrCSW3rdHbX2ErfuOnUlm7+6poqTqFBA8ozZ9RP8zlxpnj85hSL+uN8WVklcclLxEzl2na+t5eWtQev/i1qD0fvTA7LD0Po+JXaz0Pl4Hj56OSmaH2VByhJr6YIqrkTm9mDMmhxOn69hbdZJbLxrNny0Ym9R4lbzioOQlIgBHT9eyfOP+M6X3EYcpw/uxeGYeN84czsicrll6H4/quno2lR09Uwjy+vZyjlXXn9k/tG8WCycOYlpef6YO78fUvH7079V5BS5KXnFQ8hKRxsqPVfPH9WUUrCvj3T2HAbgwP5j1/vrpw7vNvIWf/fnbvP5+xZnPOdk96JGexsFj1We2jRrYi6nD+51JaNNG9GNYv44pCFHyioOSl4i0ZE9lUHq/rLCU7QeOk55mXDqxofR+WEo/RPzY23v45lMbznz+7senc9tFoyk/Vs2msiNs3neUTWVH2VJ2lJ2VJ2hIHwN7Z4YJLRidTcvrR9XJWt7ZeYj543KZMyanTfEoecVByUtEYrV1/1GWFQbPkJUeDkrvr5o6lMUz87hs8uCUqOhr7LG39/Dcxn1cd/5wbruo+YnKj1fXsW1/kMw2lR5l876jbNt/7Mw9NAADsnqk8ehd89uUwJS84qDkJSLxCkrvq87Mel95Iii9v+78oPR+/rjUK71vi9r6CEUHj/OTF9/nuY37cSDd4G+vmcxfXj4h7vOlynpeIiIpycyYM2Ygc8YM5Nsfm8rKHcGs989u2M8Ta0oY3DeY9X7JrBHMHNk/JUvvY9EjPY0pw/vxhUvH8dK2g9TWReiRkcb8cbkJ/14aeTWikZeIJMrp2npe2nqQgsIyXtp6kJr6CGNyz5beTxjSl7W7q1hVXNmue0NdUSL6pcuGcVDyEpGOcORULcs37efpqNL7/NxsSqpOEXEnM6Pt94a6K102FBFJsv69evDpuaP49NxRHDx2mj+u38cDrxVTFy7mVVsXYVVxpZJXjDp33WwREWFI3558fuFYfnrbbLIy0kgzOuzeUHelkZeISJLMGZPDY1+c3y3veXU0JS8RkSSaMyZHSasNdNlQRERSjqoNGzGzcmB3O04xCKhotVXqUb9Si/qVWtSvpo1x98FN7VDySjAzW9NcaWcqU79Si/qVWtSv+OmyoYiIpBwlLxERSTlKXon3QLID6CDqV2pRv1KL+hUn3fMSEZGUo5GXiIikHCUvERFJOUpecTCzRWa2zcyKzOyeJvZnmdlvw/1vm1l+1L5vhNu3mdm1nRp4K9raLzO72szWmtmG8OsVnR58C9rz8wr3jzaz42b29U4LOgbt/D2cYWZvmdmm8OfWs1ODb0E7fg97mNlDYX+2mNk3Oj34ZsTQp4+Y2btmVmdmNzfad4eZvR++7ui8qFvX1n6Z2ayo37/1ZnZLm4Nwd71ieAHpwA5gHJAJrAOmNmrzZeC/wvdLgd+G76eG7bOAseF50pPdpwT06wIgL3x/PlCa7P4kol9R+58Efgd8Pdn9SdDPKwNYD8wMP+d2k9/D24DHw/fZwC4gP0X6lA/MAB4Gbo7aPhAoDr/mhO9zkt2nBPRrEjAxfJ8H7AMGtCUOjbxiNw8ocvdid68BHgeWNGqzBHgofP8kcKUFS6YuIfifq9rddwJF4fm6gjb3y93fc/eycPsmoJeZZXVK1K1rz88LM7sJ2EnQr66kPf26Bljv7usA3L3S3es7Ke7WtKdfDvQ2swygF1ADHO2csFvUap/cfZe7rwcijY69FnjB3Q+5exXwArCoM4KOQZv75e7b3f398H0ZcBBocgaN1ih5xW4EsDfqc0m4rck27l4HHCH4120sxyZLe/oV7ZPAu+5e3UFxxqvN/TKzPsA/AP/SCXHGqz0/r0mAm9ny8JLO33dCvLFqT7+eBE4Q/Ct+D/ADdz/U0QHHoD3/36f634xWmdk8gpHbjrYEoVnlpd3MbBrwPYJ/2XcH9wL/4e7Hw4FYd5EBXAJcCJwEXgxXqn0xuWG12zygnuAyVA7wupmtcPfi5IYlzTGz4cAjwB3u3njUGRONvGJXCoyK+jwy3NZkm/ASRn+gMsZjk6U9/cLMRgJPAZ9z9zb9C6qDtKdfFwHfN7NdwNeAb5rZVzo43li1p18lwGvuXuHuJ4FngdkdHnFs2tOv24Dn3b3W3Q8CK4GuME9ge/6/T/W/Gc0ys37AH4Fvufuqtgah5BW71cBEMxtrZpkEN4wLGrUpABqqgm4GXvLgzmQBsDSslhoLTATe6aS4W9PmfpnZAIJfwnvcfWVnBRyjNvfL3S9193x3zwd+CHzX3X/aSXG3pj2/h8uB6WaWHf7x/yiwuZPibk17+rUHuALAzHoD84GtnRJ1y2LpU3OWA9eYWY6Z5RBc1VjeQXHGq839Cts/BTzs7k+2K4pkV66k0gu4HthOcI32W+G2+4DF4fueBNVpRQTJaVzUsd8Kj9sGXJfsviSiX8A/EtxrKIx6DUl2fxLx84o6x710oWrDBPwefoagCGUj8P1k9yVBv4d9wu2bCJLx3yW7L3H06UKCEfEJglHkpqhj7wz7WgR8Ptl9SUS/wt+/2kZ/M2a1JQZNDyUiIilHlw1FRCTlKHmJiEjKUfISEZGUo+QlIiIpR8lLRERSjpKXiIikHCUvERFJOUpeIiKScpS8REQk5Sh5iYhIylHyEhGRlKP1vBoZNGiQ5+fnJzsMEZFz3tq1ayvcvcmVlpW8GsnPz2fNmjXJDkNE5JxnZrub29elLhua2SIz22ZmRWZ2TxP7PxIuX15nZjc32neHmb0fvu6I2j7HzDaE5/yxdbOlcUVEzkVdJnmZWTpwP3AdMBW41cymNmq2B/gz4LFGxw4E/plgBdx5wD+HC7gB/Az4IsECkBOBRR3UBdburuL+l4tYu7uqo76FiIjQtS4bzgOK3L0YwMweB5YQtdKru+8K90UaHXst8IK7Hwr3vwAsMrNXgH4eLjVtZg8DNwHPJTr4tburuP3BVdTURcjMSOPRu+YzZ0xO6weKiEjcuszICxgB7I36XBJua8+xI8L3LZ7TzO42szVmtqa8vDyuoBusKq6kpi5CxKG2LsKq4so2nUdERFrXlZJX0rj7A+4+193nDh7cZGFLq+aPyyUzI410gx4Zacwfl5vgKEVEpEFXumxYCoyK+jwy3BbrsZc1OvaVcPvINp4zLnPG5PDoXfNZVVzJ/HG5umQoItKBulLyWg1MNLOxBAlmKXBbjMcuB74bVaRxDfANdz9kZkfNbD7wNvA54CcJjvuMOWNylLRERDpBl7ls6O51wFcIEtEW4Al332Rm95nZYgAzu9DMSoBPAf9tZpvCYw8B/5cgAa4G7mso3gC+DDwIFAE76IBiDRER6Vzm7smOoUuZO3eu6yFlEZHkM7O17j63qX1dZuQlIiISKyUvERFJOUpeIiKScpS8REQk5Sh5iYhIylHyEhGRlNPqQ8pmNphgVvb86PbufmcLx1wCTHT3X4bH93H3ne0PV0REJLYZNpYBrwMrgPrWGpvZPwNzgcnAL4EewK+BhW0PU0RE5KxYkle2u/9DHOf8OHAB8C6Au5eZWd+2BCciItKUWO55PWNm18dxzhoPpu1wADPr3abIREREmhFL8vprggR2Opzk9piZHW2h/RNm9t/AADP7IsHlxv9JRLAiIiIQw2VDd4/rkp+7/8DMrgaOEtz3+id3f6GN8YmIiHxILNWGBtwOjHX3/2tmo4Dh7v5Oc8e4+wtm9nbD+c1sYNQs7yIiIu0Sy2XD/wQu5uzaWseB+5trbGZ/bmb7gfXAGmBt+FVERCQhYqk2vMjdZ5vZewDuXmVmmS20/zpwvrtXJCRCERGRRmIZedWaWTpnqwcHA5EW2u8ATiYgNhERkSbFMvL6MfAUMMTMvgPcDPxjC+2/AbwZ3vOqbtjo7n/V0jcxs0XAj4B04EF3/7dG+7OAh4E5QCVwi7vvMrPbgb+LajoDmO3uhWb2CjAcOBXuu8bdD7bSXxER6eJiSV5PEty3uhIw4CbgQAvt/xt4CdhAyyO0M8KR3f3A1UAJsNrMCtx9c1SzLwBV7j7BzJYC3yNIYI8Cj4bnmQ78r7sXRh13u7vrnpuISDcSS/L6A3CTu28FMLPhwAsEI6Cm9HD3v40zjnlAkbsXh9/jcWAJEJ28lgD3hu+fBH5qZhY+EN3gVuDxOL+3iIikmFjuef0vwYPH6WaWDywnuDTYnOfM7G4zG25mAxterXyPEcDeqM8l4bYm27h7HXAEyG3U5hbgN422/dLMCs3s22HZ/4eE8a4xszXl5eWthCoiIskWy0PK/xNWF/4vwczyf+7ub7ZwyK3h1+gE58C4NsYYEzO7CDjp7hujNt/u7qXh3Iq/Bz5LcN/sA9z9AeABgLlz53rj/SIi0rU0m7zMLPrSnwGjgUJgvpnNd/d/b+o4dx/bhjhKgVFRn0eG25pqU2JmGUB/gsKNBktpNOpy99Lw6zEze4zg8uSHkpeIiKSWlkZejaeF+kMz2wEwsyvc/SUz+0RT+939D01tD60GJprZWIIktZSzD0U3KADuAN4iqHh8qeF+l5mlAZ8GLo2KJwMY4O4VZtYD+BjBPIsiIpLimk1e7v4vcZ7rIwRVhjc2dTrOJr+mvledmX2F4H5aOvALd99kZvcBa9y9APg58IiZFQGHCBJc9Pfe21DwEcoCloeJKx1NECwi0m3EupLy3wPTgJ4N2939ikZNM8Ptn29LIO7+LPBso23/FPX+NPCpZo59BZjfaNsJmq+IFBGRFBZLteGjwFZgLPAvwC6Cy3yNLUpcWCIiIs2L5TmvXHf/uZn9tbu/CrxqZk0lr3QzyyEo7vgQzSovIiKJ0lK14Wh33wPUhpv2mdkNQBnQ1HNb5xHMxNFU8urwUnkRETl3tDTy+l9gNvCvZtYf+D/AT4B+wN800X6zu1+Q8AhFREQaaSl5GYC7PxN+PgJc3uERiYiItKKl5DXCzH7c3M4mZon/UWJCEhERaVlLyesUwT2smLj7r9odjYiISAxaSl6V7v5Qp0UiIiISo5ae86ppywnNbGEs20RERNqq2eTl7vOb29eKn8S4TUREpE1ieUg5JmZ2MbAAGNxoRvp+BHMLioiIJETCkhfB3IZ9wnNGzzx/lGAWeBERkYSIKXmZ2SXARHf/ZThRbx933xndJmrqqF+5++4OiFVERASIbVb5fwbmApOBXwI9gF8DzRVhnDSz/0frs9CLiIi0SSyzyn8cWAycAHD3MppZkDIU6yz0IiIibRJL8qoJVyxuWLW4dyvtc93950Ctu7/q7ncCGnWJiEjCxJK8njCz/wYGmNkXCVYkfrCF9h+Yhd7MLqDpWeg/xMwWmdk2Mysys3ua2J9lZr8N979tZvnh9nwzO2VmheHrv6KOmWNmG8JjfmxmTS7ZIiIiqaPVe17u/gMzu5qganAy8E/u/kILh8Q6C/0HmFk6cD9wNVACrDazAnffHNXsC0CVu08ws6XA94Bbwn073H1WE6f+GfBF4G2ClZoXAc+1Fo+IiHRdsRRsfBv4VXTCMrO73f2Bptq3Yxb6eUCRuxeH3+NxYAkQnbyWAPeG758EftrSSMrMhgP93H1V+Plh4CaUvEREUloslw2/CjxvZtGJ6EvNNTazwWb2TTN7wMx+0fCK4fuMAPZGfS4JtzXZxt3rCBJkbrhvrJm9Z2avmtmlUe1LWjknZna3ma0xszXl5eUxhCoiIskUy3NepQQjnt+Z2ZPu/v9oerXkBsuA1wnujdW3P8SY7ANGu3ulmc0B/tfMpsV6cDiKfABg7ty53kExiohIgsT0kLK77zGzjwI/M7PfAb1aaJ7t7v/QhlhKgVFRn0eG25pqU2JmGUB/gtnvHagOY11rZjuASWH7ka2cU0REUkwslw3XALj7aXf/PPAKwVRQzXnGzK5vQyyrgYlmNtbMMoGlQEGjNgXAHeH7m4GX3N3DS5XpAGY2DpgIFLv7PuComc0P7419jmBkKCIiKSyWasMvNvp8P0FV4AeY2TGCZ8EM+KaZVROUzVtwmPdr5fvUmdlXgOUEE/n+wt03mdl9wBp3LwB+DjxiZkXAIYIEB/AR4D4zqwUiwJfc/VC478vArwhGi8+hYg0RkZRnwRW3JnaYPeHunzazDYQPKEdz9xkdHVwyzJ0719esWZPsMEREznlmttbd5za1r6WR11+HXz8W5zd70d2vbG2biIhIWzWbvML7RTTMEG9muQSX5/a4+9rG7c2sJ9AbGGRmOZytSOxHE+XpIiIibdVswYaZPWNm54fvhwMbgTsJ7jl9rYlD/pyguOM8YG3Uaxnw08SGLSIi57KWLhuOdfeN4fvPAy+4++fMrC+wEvhhdGN3/xHwIzP7qrv/pEOiFRERoeVS+dqo91cSzAuIux8jqOhrkhKXiIh0tJZGXnvN7KsEUyrNBp4HMLNeBAtSioiIJEVLI68vEKyG/GfALe5+ONw+n2BFZRERkaRoqdrwIE1MwOvuLwMvt3RSM5sB5Eef393/0OYoRUREosQ0t2E8whnkZwCbOHtvzAElLxERSYiEJy9gvrtP7YDzioiIALFNzBuvt8xMyUtERDpMq8nLzMaZ2dNmVmFmB81sWThze3MeJkhg28xsvZltMLP1iQtZRKT91u6u4v6Xi1i7uyrZoUgbxHLZ8DGCWeQ/Hn5eCvwGuKiZ9j8HPgtsoIXnwUREkmXt7ipuf3AVNXURMjPSePSu+cwZk5PssCQOsVw2zHb3R9y9Lnz9GujZQvtydy9w953uvrvhlaB4RUTabVVxJTV1ESIOtXURVhVXJjskiVMsI6/nzOwe4HGCqsFbgGfNbCBA1LpZDd4zs8eApwlXNw7bqdpQpJ3W7q5iVXEl88flaqTQDvPH5ZKZkUZtXYQeGWnMH5eb7JAkTs2u53WmgdnOFna7u3/g/peZNfUAs7v7nW2Ir9NpPS/pqnSpK7H0D4Gur63reQHg7mPj+Wbu/vl42jcws0XAjwhWUX7Q3f+t0f4sgmKQOUAlwawfu8zsauDfgEygBvg7d38pPOYVYDhwKjzNNeHD1x1j7zuw63XIvxRGzeuwb9Od6A9I7Jq61KX/Zm03Z0yO/vulsFaTl5llA38LjHb3u81sIjDZ3Z9ppv1I4CfAwnDT68Bfu3tJC98jnaAo5GqCuRRXm1mBu2+OavYFoMrdJ5jZUuB7BJcwK4Ab3b0sXMJlOR9cP+x2d+/4odTed+ChxVBfA+mZcEeBElgrNJKIjy51iZwVS8HGLwlGNAvCz6XAv7bSvgDIC19P0/pciPOAIncvdvcagvtrSxq1WQI8FL5/ErjSzMzd33P3snD7JqBXOErrXLteDxKX1wdfd73e6SGkGt00j8+cMTk8etd8/vaayUr0cs6LJXmNd/fvEy6R4u4nObtKclMGu/svo6oTfwUMbuV7jAD2Rn0u4cOrL59p4+51wBGg8T89Pwm86+7VUdt+aWaFZvZtM2sybjO728zWmNma8vLyVkJtRv6lwYjL0oOv+Ze27TznkIaRRLqhkUSM5ozJ4S8vn6DEJee8WKoNa8JlUBzAzMYTVUXYhEoz+wzBs2AAtxLco+pQZjaN4FLiNVGbb3f30nABzd8TPH/2cONj3f0B4AEICjbaFMCoecGlQt3zilnDSEL3vEQkXrEkr3sJ1vIaZWaPEtzLaqko406Ce17/QZDw3mylPQSXIkdFfR4ZbmuqTYmZZQD9CZNieJ/tKeBz7r6j4QB3Lw2/HgvL9+fRRPJKmFHzlLTipJvmItIWsVQb/snM1hKs42UExRcVTbUNCy++6+6L44xjNTDRzMYSJKmlwG2N2hQAdwBvATcDL7m7m9kA4I/APe6+MiqWDGCAu1eYWQ/gY8CKOOMSEZEuKJa5DV9090p3/6O7PxMmgxebauvu9cAYM8uMJ4jwHtZXCCoFtwBPuPsmM7vPzBoS4c+BXDMrIqh+vCfc/hVgAvBP4b2tQjMbAmQBy8N5FQsJkuL/xBOXiIh0Tc0+pGxmPYFsgoUnL+NskUY/4Hl3P6+Z4x4GphCMlE40bHf3f09Y1B1IDymLiHQNbX1I+c+BrxGUu6/lbPI6Cvy0heN2hK80oG+8wYqIiLSm2eTl7j8CfmRmX3X3n7R2IjN7xN0/CxwOjxUREekQzd7zMrMLzWxYQ+Iys8+Fa3n9uGFS3kbmmFkecKeZ5ZjZwOhXR3VARETOPS0VbPw3wcwamNlHCOYPfJjg4eAHmmj/X8CLwHkElxmjX7qJJCIiCdPSPa/0qOVObgEecPffA783s8LGjd39x8CPzexn7v4XiQ9VREQk0NLIKz18VgrgSuClqH0t3StT4hIRkQ7V0sjrN8CrZlZBsKTI6wBmNoHg0qGIiEhStDSC+k74MPJw4E9+9oGwNOCrnRGciIhIU1qcHsrdVzWxbXvHhSMiItK6WJZEERER6VKUvEREJOUoeYmISMKt3V3F/S8XsXZ3VYecP5b1vERERGK2dncVtz+4ipq6CJkZaTx61/yEr9unkZeIiCTUquJKauoiRBxq6yKsKq5M+PdQ8hIRkYSaPy6XzIw00g16ZKQxf1xuwr+HLhuKiEhCzRmTw6N3zWdVcSXzx+Um/JIhKHmJiEgHmDMmp0OSVgNdNhQRkZRjZ2d9EgAzKwd2t+MUg4CKBIXTVamP3YP62D105z6OcffBTe1Q8kowM1vj7nOTHUdHUh+7B/WxezgX+tgUXTYUEZGUo+QlIiIpR8kr8R5IdgCdQH3sHtTH7uFc6OOH6J6XiIikHI28REQk5Sh5iYhIylHyioOZLTKzbWZWZGb3NLE/y8x+G+5/28zyo/Z9I9y+zcyu7dTA49DWPprZ1Wa21sw2hF+v6PTgY9Sen2O4f7SZHTezr3da0HFq5+/qDDN7y8w2hT/Pnp0afIza8bvaw8weCvu2xcy+0enBxyiGPn7EzN41szozu7nRvjvM7P3wdUfnRd1J3F2vGF5AOrADGAdkAuuAqY3afBn4r/D9UuC34fupYfssYGx4nvRk9ynBfbwAyAvfnw+UJrs/ie5j1P4ngd8BX092fzrg55gBrAdmhp9zu+Hv6m3A4+H7bGAXkJ/sPrWxj/nADOBh4Oao7QOB4vBrTvg+J9l9SuRLI6/YzQOK3L3Y3WuAx4EljdosAR4K3z8JXGlmFm5/3N2r3X0nUBSer6tpcx/d/T13Lwu3bwJ6mVlWp0Qdn/b8HDGzm4CdBH3sqtrTx2uA9e6+DsDdK929vpPijkd7+uhAbzPLAHoBNcDRzgk7Lq320d13uft6INLo2GuBF9z9kLtXAS8Aizoj6M6i5BW7EcDeqM8l4bYm27h7HXCE4F+usRzbFbSnj9E+Cbzr7tUdFGd7tLmPZtYH+AfgXzohzvZoz89xEuBmtjy8HPX3nRBvW7Snj08CJ4B9wB7gB+5+qKMDboP2/N1Ilb85baZZ5SWhzGwa8D2Cf8F3N/cC/+Hux8OBWHeUAVwCXAicBF40s7Xu/mJyw0qoeUA9kEdwSe11M1vh7sXJDUvioZFX7EqBUVGfR4bbmmwTXpLoD1TGeGxX0J4+YmYjgaeAz7n7jg6Ptm3a08eLgO+b2S7ga8A3zewrHRxvW7SnjyXAa+5e4e4ngWeB2R0ecfza08fbgOfdvdbdDwIrga44N2B7/m6kyt+cNlPyit1qYKKZjTWzTIIbwAWN2hQADVU9NwMveXD3tABYGlY/jQUmAu90UtzxaHMfzWwA8EfgHndf2VkBt0Gb++jul7p7vrvnAz8EvuvuP+2kuOPRnt/V5cB0M8sO/+B/FNjcSXHHoz193ANcAWBmvYH5wNZOiTo+sfSxOcuBa8wsx8xyCK6ELO+gOJMj2RUjqfQCrge2E1QAfSvcdh+wOHzfk6AKrYggOY2LOvZb4XHbgOuS3ZdE9xH4R4L7CIVRryHJ7k+if45R57iXLlptmIDf1c8QFKRsBL6f7L50wO9qn3D7JoLE/HfJ7ks7+nghwWj5BMGoclPUsXeGfS8CPp/sviT6pemhREQk5eiyoYiIpBwlLxERSTlKXiIiknKUvEREJOUoeYmISMpR8hJJAWaWa2aF4Wu/mZWG74+b2X8mOz6RzqZSeZEUY2b3Asfd/QfJjkUkWTTyEklhZnaZmT0Tvr83XKfqdTPbbWafMLPvh+tWPW9mPcJ2c8zsVQvWXVtuZsOT2wuR+Cl5iXQv4wmmPloM/Bp42d2nA6eAG8IE9hOCtZ/mAL8AvpOsYEXaSrPKi3Qvz7l7rZltIFjM8Plw+waChQsnEywW+kI4M346wdIgIilFyUuke6kGcPeImdX62ZvaEYL/341g/ruLkxWgSCLosqHIuWUbMNjMLgYwsx7hGmwiKUXJS+Qc4sFy8jcD3zOzdQSz/y9IalAibaBSeRERSTkaeYmISMpR8hIRkZSj5CUiIilHyUtERFKOkpeIiKQcJS8REUk5Sl4iIpJy/n9Vi7sABd8HSQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f = lambda t, y, params: -params[0]*y[0]\n",
    "params = (5,)\n",
    "y0 = np.array([.01,])\n",
    "t0 = 0\n",
    "tFinal = 0.12\n",
    "h = 0.1\n",
    "ETOL = 1e-4\n",
    "t = [t0];\n",
    "steps = [];\n",
    "# start with y as a list of nparrays.\n",
    "y = [np.array(y0)]; # list of 1st order solutions\n",
    "attempted_t = [t0,[]]\n",
    "attempted_y = [np.array(y0),[]]\n",
    "while t[-1] < tFinal:\n",
    "    Yi = y[-1];\n",
    "    YFE = Yi + h * f(t[0],Yi,params);\n",
    "    YET = Yi + h/2 * f(t[-1],Yi,params) + h/2 * f(t[-1]+h,YFE,params);\n",
    "    attempted_t[-1].append( t[-1]+h )\n",
    "    attempted_y[-1].append( YFE )\n",
    "    err = numpy.linalg.norm(YET-YFE,2);\n",
    "    hnew = h * np.sqrt(0.9*ETOL/err);\n",
    "    # Make sure h doesn't more than double\n",
    "    if hnew / h >= 2:\n",
    "        hnew = 2*h;\n",
    "    if err <= ETOL:\n",
    "        # The step size was small enough.  Record this\n",
    "        # info and determine the next step size.\n",
    "        t.append( t[-1] + h );\n",
    "        steps.append( h );\n",
    "        y.append( YFE )\n",
    "        h = hnew;\n",
    "        if t[-1] + h >= tFinal:\n",
    "            # If we are going to step out of bounds, then\n",
    "            # restrict the size of the last step.\n",
    "            h = tFinal - t[-1];\n",
    "        attempted_t.append( [] )\n",
    "        attempted_y.append( [] )\n",
    "    else:\n",
    "        h = hnew;\n",
    "\n",
    "# The first step taken from any time point will be blue and the second orange.\n",
    "plt.subplot( 2, 1, 1 )\n",
    "for i in range(len(attempted_t)-1):\n",
    "    plt.gca().set_prop_cycle(None) # this resets the color cycle\n",
    "    for j in range(len(attempted_t[i+1])):\n",
    "        plt.plot( [t[i], attempted_t[i+1][j]], [y[i], attempted_y[i+1][j]], '.-' )\n",
    "plt.ylabel( \"Y\")\n",
    "plt.subplot( 2, 1, 2 )\n",
    "for i in range(len(attempted_t)-1):\n",
    "    plt.gca().set_prop_cycle(None) # this resets the color cycle\n",
    "    for j in range(len(attempted_t[i+1])):\n",
    "        plt.plot( t[i], attempted_t[i+1][j]-t[i], '.' )\n",
    "plt.xlabel( 'Time')\n",
    "plt.ylabel( \"Step Size Taken \\nfrom that Time\");\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluate rejected steps for the Goodwin oscillator\n",
    "\n",
    "Here I repeat the code from above, but use it on the Goodwin oscillator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaIAAAEYCAYAAAAeWvJ8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/4klEQVR4nO3deXyU5bn/8c81k4R9SQg7BIgCIogsERAVVxSXirZuoK1LlerRHtuepT3n/LTW9pxjN09taxe07guKK66IYCuCIAmgrMGwBBLCFsIaIMtcvz+eSTIzzJbJMpOZ6/16+co869x5HOabe3nuR1QVY4wxJl5c8S6AMcaY1GZBZIwxJq4siIwxxsSVBZExxpi4siAyxhgTV2nxLkBLy87O1sGDB8e7GMYYkxIKCgr2qWrPxhyT9EE0ePBg8vPz410MY4xJCSJS3NhjrGnOGGNMXCV9jcgYY3zNLZzLx9s/5pKcSwDqX18//Pqot10//Pp4/gpJR5J9ZoW8vDy1pjljkltgSNQtXzjwQvZW7mXRjkVM6DOBw1WHeWfLO0HP0b9jf0orS4NuG9ptKF8f/Lp++ebTbqayppKiA0VclXsVM0fMDFmGVAsuESlQ1bxGHWNBZIxpi+YWzmVB8QI6pXfi4+0f16/v07EPuyp3Nfp8gqDE9n3oxk0ttfXLp3Y9laJDRfXLD0x8ABFJiWCyIArCgsiYtq2uZnHxwIvZWLGR5WXLceFi2+FtQfd34cKDJ+i2/p37U3okeK3nyiFX8t7W94JvG3wl721r2Na3U1/KjpbVL6e70qn2VEf5G8H1w66nc3pn1pev57LBlyVVMMUSRNZHZIxJKHXBc/6A89l8YDOvbnoVgKU7l0Z1/OVDLg8ZKN8d9V2AkP1A43uPD72tj/+2h5c9XH/eSwdd6veelw26jPnF8+uXu6Z35VD1oYbfcdPc+tfLdy3nnS3vkNkuk22HtnHd0Ov49shv+12LpK9FWY3IGBNvdc1saZLG4p2LG338lUOupOJERVSDDpqzzOH6hALL4BtcgTWzwFpcn4596N+5PwV7CurXPTjpwTYRRtY0F4QFkTGJaW7hXD4q/ohaTy0rdq8Iuk+4/p7A8El04YIpt2suWw5tqV/ultGNI9VHqNWGfqfMdpl85/TvUHSgiL3H9jJt8LSE/L0tiIKwIDImcdSFT01tDfl7Iv+7fHDSg0DoprS2LFwwPTjpQVSVny//ef26rPZZ7D++3+8cFw64kNvPuJ31+9bzaemnCXFd2nQQichTwFXAHlUdFWT7BcDbwFbvqjdU9eHA/QJZEBkTP74DDYoOFPFy4csRj2lrNZ3mEqw/KHDdrR/cyso9K8Oe5+7Rd/NPY/6J1za9FpfQbutBNAU4AjwXJoj+VVWvasx5LYiMiY+5hXP9/soPJ1XDp7ECr+m/5v0rb2x6w69ZD6BzWmeO1BypX27N/qU2PWpOVT8VkcHxLocxpmle2fgKL298mV1H/ft2sjtks+/YvvplC5/Gq7tOvjWdjmkd/cLpqtyrWFK6BGoajvu/gv+j4kQFldWVrC9fz9RBUxPqmidMjQjAG0TvhqkRvQ6UADtxakfrQpxnFjALICcnZ3xxcaPn4DPGNMLcwrnM3zafo9VHWVu+Nug+gf09ifRF2NYFG7HnG06B9z0BzBg+g385619o527XrGVp001zEDGIugIeVT0iIlcAj6nq0EjntKY5Y1rW02ue5tGVjwbdltstlz6d+ljwxEFgON3+4e3k7/b/LuyQ1oHJ/SYDUH6snLG9xtK1XVfyeucxpteYmN63TTfNRaKqh3xevy8ifxKRbFXdF+44Y0zL+NUXv+Kj4o/YXbk75D63jLjFAihOrh9+vd+1v2LIFX5BNPO0mdRqLe9veZ/D1YcBWL13NQDt3O148tInYw6jxmozQSQifYDdqqoiMgHnERblcS6WMSmj7i/ssb3GMn/bfIoOFAXdz/p+ElOw/iWA7Ye283nZ53771nhqyN+dn3pBJCIvAxcA2SJSAvwUSAdQ1b8A1wH3iEgNcAy4SROpXdGYJObb5xBsqh1rgmsbAmtJAFMHTfULIkFId6WT17tRrWtNkjBBpKozImz/I/DHViqOMQYngN7b+h6b9m/yWz+oyyCKDzcMArImuLbLt6Y0PHN4k/uIYpEwQWSMSSxPfPkEv1/9+6Dbbh15K2Aj4JJFsJpSa7IgMsbUq5uCp3N6ZxZtX+S3LVjzmwWQaQ4WRMYY4OS79nt26MneY3vrl635zbQUC6IQqmuqmPfhU+zcsZnR487j/IlXx7tIxjS7upFwp/c4nVc2vuK3bWjmUO458x5rfjMtLqFuaG0Jsd7Q+uJrj7Jr7iLU+/Dg2jN6M27KZRx31VBQspzzhl7MJeOnN3+BjWklgTWgwKeMtpXn35jEktQ3tLa23Tu2oCiCAEr6mj0s3PEkH07chccFb371CZOL3+LbE+9kYp+JvFn0pv3laNqMak81z61/zm/d2F5jmTZ4mn2OTauzIAph/ORLWfbFX3F5FI8LTr/jejZvegePaxcIKLDkcD5LPvavbS3duZSFxQu5e8zdrTr80ZhozC2cy4sbXuRg1UG/CUiB+getWQCZ1mZBFML5E6+GH8CaVZ9xxthzOX/i1XTI6s7SLx/AI4pLhdtzZrK73SE+Kv6IE7Un6o9dUraEz8s+Z+aImVxz6jV8uedLFu5YaH9lmrh6eu3TPFrQMCfcDcNu4LSs06wGZOLO+oga6eOCt1n89UK/PqLGPHflnjPv4Z4z70FEmq1MxoQyt3AuC4oX0L1dd+ex3D6Pnp7cbzJ/nfrXOJbOJKM2P/t2S2it2bfnFs7ljaI32Fi+EUVJd6XzyHmP8PtVv2fLQf+HVvXq0IsJfScwse9Eig8Vs758vf1Fappd4B9IWe2y2H+i4VHTNhjBtAQbrBBHdW3rq/esJn93fv0UGfuP7z/poVU1nhqWlC7h3S3v1q9funMpX1d8zQ/G/4CO6R3j8SuYJFJZXckTa57wWzc8azhTB021pjiTcKxG1AqCPYveox5mvjuTdfv9n+2X5kpjXK9xnN3vbA4eP0hhRWHCPU3RJK65hXN5fsPz7Du2j8NVh/22WQ3ItAZrmgsiEYIolMCmk1tG3EK6K50lO5ewqcJ/kslvnvpNfjj+h3Rv372VS2nair+t+Ru/W/m7+uXbR97OwC4DrQZkWpUFURCJHEQQvLYEBH2aoktcjMoexbn9zuXc/ueyYf8GFm630XiprG4wQteMrizYvgCPeuq32WAEEw/WR9QGhbpvI/BpineMuoN27nZ8VvoZf/7yz/zpyz/Vb1u6cynHao7xnZHfaZUym8QQWKPObp/NvuMN9wZdknNJPIplTKNZjSiBhaot7T++nzvn38nXB772239cr3GcN+A8pgyYwtDuQ22IeBI7Wn2Ua9++lrKjZfXrzu57tg1GMHFnTXNBtOUgCifwr+Fz+59L+bFyNuzfAEDfTn2ZMmAKUwZMYcehHXxa+ql9OSUBj8fDvYvu5YuyL6jyVPlts8EIJhFY01wKCfX8+d1Hd/NZ6Wd8WvIp8zbP45XChhmVl+5cysETB7lz9J1xKbNpmpLDJcz6aBY7juyoX3flkCupOFFhf2SYNs1qREmsqraKb7//bdbvX++3/rSs05gyYArnDzifUdmjcIkrTiU00Xhh/Qs8t/459lbuxYPHBiSYhBZLjci+gZJYhjuD64Zd57fukpxL6JjWkSfXPMnN79/Mha9eyANLHmBh8UIqqyt5YcHj/OTXM3lhweNxKrXx9cjyR/jlil9SdrSMGq3hggEX+G23AQkmGSRM05yIPAVcBexR1VFBtgvwGHAFUAncpqorW7eUbU+oJryDJw7yWeln/GPHP1hYvJC3it6iV0V7pi3rRbZC6cr3eAG4Zeq9cSx96np0xaO8uflNDpw44Lf+eO1xHpz0oA1IMEklYZrmRGQKcAR4LkQQXQF8HyeIJgKPqerESOdN5aa5aFV7qlm9ZzVPPvEAQ9e6vc9gguMZtXQYewpjJ1/CuXlXkJ6WEeeSJre5hXOZv20+eyv3suXQlqD72IAEk+ja9GAFVf1URAaH2WU6TkgpsExEuotIX1UtC3OMiUK6K52z+pxF4aRvULr+PVweUIGqTm4yvtjKyuVP8nnGbBiSxVGpIr1W6DdqFK6Rfevn1DNN88rGV/jF8l8E3ZbbLZc+nfpYDcgkrYQJoij0B3b4LJd4150URCIyC5gFkJOT0yqFSwa3TL2XF4C1q5cwasw53DL1XvZWlLHo09coKljOoV1l9DzYHoCDXy+j6IsjvDB0No9e+xd2fJHP2lWLGTX2PKZfPSu+v0gbMrdwLnM3zaXoQFHIfW4ZcYsFkElqCdM0B+CtEb0bomnuXeARVf3Mu7wQ+LGqhm13s6a55vOLh26h3YYKXIjPY9RhV/fj9DyUgXgEj0s5cckQrv7GXVZTCqHuRmVRYUnZkqD72LBs01a16aa5KJQCA32WB3jXmVZyxrgpbCx8CzzgcSkrhlfgcrsYsKc9Lo8gCOKBms+K+N/S+7jlin/Gs/MABfkLSRuUbeFE+IcoWhOcSVVtKYjmAfeJyBycwQoHrX+oddU1ua1dtZiew4cxzdtHVLzsCza+/BYuDyCQVitMWNedZcVPkHkkg66AZ9MRfl52L5edOo29hZtSsgnvgc8e8HsGVSBrgjOpKmGa5kTkZeACIBvYDfwUSAdQ1b94h2//EZiGM3z79kjNcmBNc63l7XmzKchfyJoOJWztf5TuRzMYX5hJ37K0+ia8/Z1P0LUyHbe3Ca/rNyZw98yHeHve7KTtX3p146vMKZzDrqO7OFx9+KTt1gRnko3NNReEBVHr8n1CbUNNSVCBndnH6L+3Ay5vMNW4POzr66HnTjeiTnPfaTOuYXflborW5HPqGXnMuumBOP9GsXnyqyd5fv3zfo/m9tUtoxv3j7vfwsckHQuiICyI4quuppQ2KJs+nfpQ8fbn9cFUlnUMlwr99ncAQFHKMo/R+0CH+mCqGN0Zz5ETCR1KdYMPuqd3Z+OBjRytPsruyt1hj7H7gUyyimsQicj7wD+p6rZmOWEzsSBKLLPn/JyC/IXsyjrOzpxaLjh4GtmflePyOLWk3Zkn6FvhDBHf3/kE3Y5m4FLwCBw9IxM9UV0favHua5pbOJcX1r8Q8ubTQNYMZ1JBvIPoeuC/gWeBX6lqdbOcuIksiBKPb/MdwM/+di+99qWxJ7uGib3OQj7ehMsj7M46Tp/97RGE/Z2r6HY0HZdCRecqulam4/L2Ne0+tzsH253g3FMvpHPf3i0y/U1drefigRdTcriEeUXzKK8qj+rY3G65NhDBpIy4N82JSGfgAZwBBc8D9dMEq+qjzfZGjWBBlPh8g2lMrzHMnvNzitbk4+rcjsyvjtSHUu/97XEhlGUdqw+oPd1P8OGkXXgEUPym8Z2aM5WjNUfrJwb9ePvHZLbLrK+VBK4L3Da/eD4ej4eiA0VUnKho1O9ktR+TqhIhiDKAnwAzgVfwD6KfNdsbNYIFUdtWF0qdu3Sn3Zd7cXmE/V2q6H7UqRGtOeUgq4YdhLogitNDaXO75TIia4SFj0l5cb2hVUSmAY/i3O8zTlUrm+vcJnXNuukBuMl5XTfwocOgbHC3o2jtSg52rHaGg4uCgrpbt3zW7GZM0zVnH9Fi4G5VXdcsJ2wmViNKXnVNeuXbi9lUtt6vjyizXSbvbX2vWd+vLnTg5MdqGGMccW+aS0QWRKmrboBBrH1Egftb6BgTmQVREBZExhjTeuxR4cYYY9ocCyJjjDFxZUFkjDEmriyIjDHGxJUFkTHGmLiyIDLGGBNXFkTGGGPiyoLIGGNMXFkQGWOMiSsLImOMMXGVMEEkItNEpFBEikTkJ0G23yYie0Vktfe/O+NRTmOMMc2r2R4D0RQi4gYeB6YCJcAKEZmnqusDdn1FVe9r9QIaY4xpMYlSI5oAFKnqFlWtAuYA0+NcJmOMMa0gIWpEQH9gh89yCTAxyH7fEpEpwCbgh6q6I8g+xpiWkP8MbHgbRkyHvNv8lyH467zb4lNW06YkShBF4x3gZVU9ISLfA54FLgq2o4jMAmYB5OTktF4JjUl04cJk3Hdg6WOw4T0YdimcOhWW/RF2r4cOPaB4sXOOzYtg1fNQmt+wXCfw9VevwNG9MHQqjLkZVj4PO1fByOngckPhh3D6Nc7+oYIssMwm6STE84hE5GzgIVW9zLv8HwCq+r8h9ncD+1W1W6Rz2/OITErx/dIefyss+R2smQsDz4aqI/DVnIZ9O/WCo3viVdLoZA2D/ZsalvuMAa2Bkd+EtPZQtODkILPwiqs2+2A8EUnDaW67GCgFVgAzfR87LiJ9VbXM+/pa4MeqOinSuS2ITNLx/ZLVWlj1IgycCJV7YM1rDfuJ29keSuD2tPZQc7xhOb0TVB8NfuywabDpw8aVO9z5mlOPoVD+dcPy4ClQdRTOvAkmzgrdpGiB1SzabBABiMgVwO8AN/CUqv63iDwM5KvqPBH5X+BqoAbYD9yjqhsjndeCyLR5dV+cQy+FfUWQ/2R0x7nbQe2JhuXMXKjY0rB8xg2w5tXGLVfui76PqGN2+PNFI5ZjwgkXzpPuBU8t7F4LZ3gfC2+1rEZr00HUUiyITJvh+0VXXek0o6nArtXRHZ+VC/vDBM1Vjzk/w325Rlpuyu8U7QCHUGUMDLZQwoWXOwNqqxr/e2QPh32FDcuX/AwyOkPhexZUASyIgrAgMgmt7ssroytseCvy/l36weGdwbdFEzRtXbQj9UKFV7iQ6joADpXEVq4u/eFwacPy5H8GBMpWO/1Zyfj/IgQLoiAsiEzCyX8G1r0BCGz9e+OODQwb39dJ/OXWJI2plb17f/BzBAZYp95wdLfPDi7AE7oMHbLg2P6G5fG3ObXd3WthzAw467tJE1QWREFYEJm48v1y2bfJGc5cua9x5wjsmzEtJ9xAhsBtvqEVGFSBNddIA0fSOkJNZcPyqG85/xV9DHvWw+ibIO/2k8uYgJ8HC6IgLIhMXOQ/A8sed8InWr6BAwn9ZWOIXNMKF1TdB8OBbQ3L0fRddcyGdl39B5xMvBvG3gJffwzbPvX/vMQpsCyIgrAgMq0m/xnnnp2a4w03e0ZitZ3k1ZSgmvpzWP+2/+eoUy9nEEvVkfDvm30adOgOO5Y1rLvgP2Dy9+GruS0eThZEQTQliHZ9/jmlqzfTf8wp9Dn77GYumWnz6r5YxA2lBf59AJH0PA0m3mPhk8qiGbHoG1Z1/YO+6865Hwo/8K95Z3Rx/hjyVId//5zJcMqF0LUfFC+BXetg8LnQqQcMPg8GTojp17IgCiLWINr6wXw+eNuFIrioZdK4ffQfN4zqXZvRnavpcOal9JjoM8NQ3Yeoz2ho37VJ/yNNglvxFHz2OzhYHN3+daED1txmGidY81osgTXm27B9if/wfld66LBKaw+3vhPTd5gFURCxBtGnv3yKNVsHAVK/rnf6Rq7JeoA0qUEVaqU9ae3bO/scr/A/gbsdXP+McxPiqhdCfwEleMdjyqv7/9OuO+xZC9XH4eD2yMdZk5tpTbEG1pk3wfPXwPbP/c8nbrjov+C8f2l0UWIJorY06WmrGjZlBBu2HaBW03BRy3nnHSFjxxpcxxpGvlRLZ9I8x53pQwLVnoA5M/zXbV4EK5+BgZOgc09nMsm1rzVsq62GiXc17G8h1boC2/CXPObfMRwNa3Iz8ZB328mfucB1da8Dv1NG3xQQRC5n8MTg81qosCezGlEYgX1E5csX0fX9G3FTQy1pHLriFad5LvAvDQBXBky6G9bP8x8dk9YeXGmhOxzbdYVO2aAKFVsb1o/8Foy4Ejr1dEbPdOoJG+bBxnctqJoq/xlY9if/O+ej1fM0pznWaj+mLWvGrgVrmguiuUfNlS9fxPG1n9B+1IXR9REFqw7n3ebUopb8Af7hM8H48Cugew4c2ePUkI4fiL5gmbnQZyR07OEEVcceTqB17AEb33fu8B59A0yYBSIhT7PrvZco/WoL/Ufn0ufKmdG/f6LzvdO+ch+ceimUfQlFH0FleePOZc1uxoRkQRREQgzfDtfEFmpbYIBN/QUMvcR5tsvRffD3//UfKdO+O3Tp43ypVpaDhrjL2+V27grvmAUdMr0/s6BjFru2HOLNgql4cOOmhgsn7mDA1Gm0y8wirVMnAArfWcTO1Zvpl9uB4Tm76wM31PrWUP7B3zi+7u+0H3kBPXq6/ZvX1r7u1EC3LAp/knCs1mNM1CyIgkiIIIpVpAALVtMC8Hic2lRlObwxC3aubNivWw7kToHK/c5/x+p+VlBweDrLjtyC7wCNOm6q6JG2lT01QwFBUPqkryczrQxPhyw27T8DD05/2uiO79G//QbSJ9zC9k2VHCw/Qd8ehxg6dSJbNlWxe91W+vU6wojplzgh9uIL7Fy3k34j+zH85lv8loH6156OvdlZsIFeWYeoTcukbEc1PdOLGJM+Bxc1KOASDVL6RrARbsY0iQVREG06iCKJZjBDuMDy5fGw651neOuDftTixoWHvKFf06FnL04creJ4ZQ0VO8opPj4GJ6iUdnIYt9TQ3V3KzuqROPNtNejh3kx57Sn1y5nubRys7Y+HdACEGnq4N7Ovdlj9Pn3T1lBWc4Z3/2IO1A5EcZHpLuZg7QA8pPmdZ1yn15nU+YX61kYlWIxGwQYZGNMsbNRcqgk2UibYPhA5sFwu+ky/g2vSfPuI7vPbpfCdRZS8dwIPbmckYdenGN7lCwqH/Ibdi6vra0QTO79Avw6bWV11PeUHPdRNCNnFvY+K2oHes3nolbaZDu6D7KvV+n3c4sGJExcdXIepqBVAfF7jfe2EXmnVCGrJwKU1eBDcEmI+r7p+nbo+ohaeRuel5dv5YG0ZPTplUH60istH9WXmRHtsvTHBWI3INEpj+ogKVx/hE5/gGt3xXb6qvKp++cLusyH3fD5ZObx+3Zj+q1hdOhYPbrLcxRyoHRD2tYtazu86m65ZabQ/c9rJfURNDJpggQLwwdoyurVPY8/hE7RPd1N+tIqu3uWeXTJYtqXipHON7t8Vj0K3DmkcPFbDlaP7oh5lcdE+rhrdDxHhg7Vlfu8R+HrmxJz6Mlm4mURkTXNBWBDFV+E7i9iZv45+vY4w/OqLKVx9xG+5MX1EdO1XfyzdB7JzU4XfeWIVGDaXjuhN6cFjvLWqlF2HTkQ+QSsa0acLG3Ydrl+enJvFweM1XDW6D7efk8sbBSV8sG5XyNCyEDMtzYIoCAsiAyeHzWWn92bXoeO8XlBCWYKFTXPp1j6Ng8dr6pdzszuyZV/DowZuyBuAoqwsPsDUEb2ZMTGH977cWV9Du3nSIL/ggtA1s8BtJnVZEAVhQZRaAr8c3/myFLdb+OzrRt4r1ETXjOlH+dGq+uDr0SmDt1aHeLJqI1x+Rm8+WLM76Lb26S6OVzcM209zCTWe2P99uwVqQxw+MLMDOyqOBd02Nqcbh47VMG1kb1wuYfHX+7h8VF9uO2cwr63YwYfrd4dtfvR9HS7wAvczicGCKAgLorbt7V//jf1Ll5E1eRJ9urVny8eLOTB0FP26d6ByxQr2545gX9eebCw7xLH+g1h9IMxTMmN0aq/OjOrX9aQ+It+Q8f0Z6osxXH9TLF/MgeF2zZh+jVoelNWR4v0NNaSOGW4qqxoGe2S4hapQSdRKhvToyNbyyoj7Tc7NYn9lNZeM6MWNZ+XwWsEOPt9SzhVn9CVNhPne8ItUk4t2WyzLqaLNB5GITAMeA9zAk6r6SMD2dsBzwHigHLhRVbeFO2dTgmjVB59S+vcl9L/gHMZePgWApQuWs2TlZs4ZdwqTp06s37fuQzeyb1e6dEhnUm4Pxg/K9NsW7AMZ67Zotjd2P4C3XvyQnX9fQr8LzuGam6edtH357Bcp/ccS3GPGU5J3fv3vGW59XXhMmHUzNfPf91sePyjTL2wu/+dv8+bDf+TYqlWc6NSFSesWk+apYWuXPgw8upc0Ty1bu/Qh58he0rSWDZmD+Pfz/olalzvs7xVJqLBJ9C+Rpnw5Avznm2vqzxUpuHyF2xZKuzQXJ2qa/w+FxghXy+vRMZ3yyuCzUQc2awb21Y0e0JWvSg7VL08YlMkXxQ0DVu6ekkv50RNs3HWYi07rzTVj+/Pqiu0s37qfy0b24Ya8gby2cgeLNuzhilF9uOXswYhI0H+7zb3Od32w76/GatNBJCJuYBMwFSgBVgAzVHW9zz7/BIxW1btF5CbgWlW9Mdx5Yw2iv8+eQ/ajD+NC8SBsHzCM4m59eTT3MqpdabjUQ057Ja17dw4dq2LPYf+nK7pdMP3M/lRUnuCTwoZHQ196em/G5mTiEsjfVs6CDXvrt119Zj8m5mbhEmHF1nLeWNXwD/36vAGce2o2IoIAS4v28fKKHfXbb56Yw/nDetZvFwGXCH/ftIdnlzY8ruC2yYO48LTeJ91rIwJffvgp5/35IdyeWmpdLtZdfjODxp6OejygsCt/FcM+fgO31lIrbt4ecg7bs/oz2X2QcSs+Is1TQ624+HhgHru79GSYHmbChs9wqQcVF8t7Dydvz9eke5yh1sVdelPRoStj92yqL8/arMGcsX8bAF/1yOWM8i2I9/Wo8i24gC975DKqfCtulDlDL+K5EdNQlws8tSCuoFMYBTaVtaWwaUmNDa5I28I1QcYSXrEcE1iz8xWulhcupOLBuVuvQccM548t398tu3MGHo+y3ydAc7I6IgLFPjXJ4X264AK/8BwzsDun9urMtn1HyS/2H+XZPs3Fi3dNiimM2noQnQ08pKqXeZf/A0BV/9dnn/nefT4XkTRgF9BTw/wSsQbR27P+jaGfvlv/YTjUrjNv557LnGEXoeICVdI9tdSkpRHq3V0CTWiib3U3bFrIbes/CHlDqG8wBK6vCwlfX/bI5cxyZ/ZqBb7MPoUz9m3BjaLAng7dKe7ah7N2b0QAD7Cy1zDG7dmEC9jcpQ8Dj+7D7a0FOa9r2Obzen3WYB445y5qxE2a1nLFoE6Ud+hmQRNHsTR9hXodzTGRmih9Nde2b5zZl3e+LKtfnjqiFws27KlfnnJqNp8WNfwB2q9be3YePF6/HBiWgf17OVkdUVW/fri+3dpT6/H4/dGb3TmDGo9ywCeIunVIR1U55DNQpZM3xI76vGe7NBc9OmWw78iJk8LZLfCjS4dz74WnBr0e4bT1G1r7Azt8lkuAiaH2UdUaETkI9AD2+e4kIrOAWQA5ObF9AeVc+w2qlszHXVtDrTuNdo/8lm+kteONj3ZSLW7StZanL+3H5KkTeWn5dr8mDgHapbt48c5JFO467Lft4ekjuSFvILUeZc4X2/n5exvqt/3H5adx7dj+eBTeWFXCrz5smA36R1OHcdXovngUVJV3vyrjsYVf12+/98JTuHxUX1RBUe9P+HBtGX/5R8OjDO6eksvUkb39fte6IM1/v5rqjQu8NSI36667i5HnjUfEhbhcdFr4KdVvbq/f/vSIy1k5cBT3dtxL9VvbSfPUUuNy8+fR01k65Cwe7raLE889jttTQ60rjcKsHEbs3456l18YdQXnnJJN1aub6/epHDKM6n1bcXtqGHC0nGVnTSPj0AGyJk/ihE8fkcfbR3Rg6CjuGdSeLZtLmTzuFGbcdKHf72YB1PpmTszxu+6Br8Nti/WYCUN6+IWX7zLQ7NtmTszh7NzYm0EvPb23X7BNG9nHb/nu80856ZjvXzT0pHU/mjr8pHU/nnbaSev+68rTT1r302+MrC9n4PdXepqLSbk9aC2JVCO6Dpimqnd6l78NTFTV+3z2Wevdp8S7vNm7z75g5wTrI2rMfhD/PqLp//bdk5aNSQaxDG6wPqJWlmhNc8YYYxovliAKbNaPpxXAUBEZIiIZwE3AvIB95gG3el9fBywKF0LGGGMSX8LUiABE5ArgdzjDt59S1f8WkYeBfFWdJyLtgeeBscB+4CZVDfssZxHZCxSH2yeCbAL6oFKUXQeHXQeHXQeHXYcGdddikKr2bMyBCRVEiUhE8htbzUxGdh0cdh0cdh0cdh0aNOVaJFLTnDHGmBRkQWSMMSauLIgimx3vAiQIuw4Ouw4Ouw4Ouw4NYr4W1kdkjDEmrqxGZIwxJq4siIwxxsSVBVEIIjJNRApFpEhEfhLv8sSLiGwTkTUislpEUmqKChF5SkT2eKeWqluXJSILRORr78/Y5kFpQ0Jch4dEpNT7uVjtvQcwqYnIQBH5RETWi8g6Ebnfuz6lPhNhrkPMnwnrIwoimkdSpAoR2QbkhZvPL1mJyBTgCPCcqo7yrvsVsF9VH/H+gZKpqj+OZzlbWojr8BBwRFV/E8+ytSYR6Qv0VdWVItIFKACuAW4jhT4TYa7DDcT4mbAaUXATgCJV3aKqVcAcYHqcy2Ramap+ijODh6/pwLPe18/i/ANMaiGuQ8pR1TJVXel9fRjYgPNEgJT6TIS5DjGzIAou2CMpmnSh2zAFPhKRAu/jNVJdb1WtexDNLqB3uJ2T3H0i8pW36S6pm6MCichgnKnGlpPCn4mA6wAxfiYsiEwk56rqOOBy4F5vM40BvBPupmrb9p+BU4AxQBnw27iWphWJSGfgdeAHqnrId1sqfSaCXIeYPxMWRMGVAgN9lgd416UcVS31/twDvInTbJnKdnvbyOvayvdE2D8pqepuVa1VVQ/wBCnyuRCRdJwv3xdV9Q3v6pT7TAS7Dk35TFgQBRfNIymSnoh08nZGIiKdgEuBteGPSnq+jyK5FXg7jmWJm7ovXq9rSYHPhYgI8Ddgg6o+6rMppT4Toa5DUz4TNmouhGCPpIhviVqfiOTi1ILAeaz8S6l0HUTkZeACnOntdwM/Bd4CXgVycB4vcoOqJnVHfojrcAFOE4wC24Dv+fSTJCURORdYDKwBPN7V/4nTP5Iyn4kw12EGMX4mLIiMMcbElTXNGWOMiSsLImOMMXFlQWSMMSauLIiMMcbElQWRMcaYuLIgMsYYE1cWRMYYY+LKgsgYY0xcWRAZY4yJKwsiY4wxcWVBZIwxJq7S4l2Alpadna2DBw+OdzGMMSYlFBQU7FPVno05JumDaPDgweTn58e7GMYYkxJEpLixxyRM05z30bJ7RCToMyzE8XsRKfI+inZca5fRGGNM80uYIAKeAaaF2X45MNT73yycx9K2qILiCh7/pIiC4oqWfitjjElZCdM0p6qfisjgMLtMB57zPhN+mYh0F5G+LfUwroLiCm5+chlVNR4y0ly8eOckxg/KbIm3MsaYlJZINaJI+gM7fJZLvOtOIiKzRCRfRPL37t0b05st21JOVY0Hj0J1jYdlW8pjOo8xxpjw2lIQRU1VZ6tqnqrm9ezZqMEb9Sbl9iAjzYVbID3NxaTcHs1cSmOMMZBATXNRKAUG+iwP8K5rEeMHZfLinZNYtqWcSbk9rFnOGGNaSFsKonnAfSIyB5gIHGyp/qE64wdlWgAZY0wLS5ggEpGXgQuAbBEpAX4KpAOo6l+A94ErgCKgErg9PiU1xhjTnBImiFR1RoTtCtzbSsUxxhjTSpJysIIxxpi2w4LIGGNMXFkQGWOMiSsLImOMMXFlQWSMMSauLIiMMcbEVcTh2yLSE7gLGOy7v6reEeaYc4Ghqvq09/jOqrq16cU1xhiTbKK5j+htYDHwMVAbaWcR+SmQBwwHnsa5KfUF4JzYi2mMMSZZRRNEHVX1x40457XAWGAlgKruFJEusRTOGGNM8oumj+hdEbmiEees8s6CoAAi0immkhljjEkJ0QTR/ThhdFxEDonIYRE5FGb/V0Xkr0B3EbkLp0nvieYorDHGmOQTsWlOVRvVrKaqvxGRqcAhnH6iB1V1QYzlM8YYk+SiGTUnwM3AEFX9uYgMBPqq6hehjlHVBSKyvO78IpKlqvubq9CmcQqKK+y5SsaYhBVN09yfgLOBmd7lI8DjoXYWke+JyC7gKyAfKPD+jEhEpolIoYgUichPgmzPEZFPRGSViHzVyL6rlFRQXMHNTy7jtx8VcvOTyygoroh3kYwxxk80QTRRVe8FjgOoagWQEWb/fwVGqepgVc1V1SGqmhvpTUTEjRNwlwOnAzNE5PSA3f4f8KqqjgVuwglJE8ayLeVU1XjwKFTXeFi2pTzeRTLGGD/RBFG1NyTqRsH1BDxh9t+M8+C6xpoAFKnqFlWtAuYA0wP2UaCr93U3YGcM75NSJuX2ICPNhVsgPc3FpNwe8S6SMcb4ieY+ot8DbwK9ROS/getwaiah/Aew1NtHdKJupar+c4T36Q/s8FkuwXkkuK+HgI9E5PtAJ+CSYCcSkVnALICcnJwIb5vcxg/K5MU7J1kfkTEmYUUTRK/h9PNcDAhwDbA7zP5/BRYBawhfc4rFDOAZVf2tiJwNPC8io1TV731UdTYwGyAvL0+buQxtzvhBmVEHkA1sMMa0tmiC6A3gGlXdCCAifYEFwPgQ+6er6o9iKEspMNBneYB3na/vAtMAVPVzEWkPZAN7Yng/E6BuYENVjYeMNBcv3jnJwsgY0+Ki6SN6C+cmVbeIDAbm4zS/hfKBiMwSkb4iklX3XxTvswIYKiJDRCQDZzDCvIB9tuPUzBCREUB7YG8U5zZRsIENxph4iOaG1ie8wfAWzgzc31PVpWEOmeH96RtWCoQdOaeqNSJyH07QuYGnVHWdiDwM5KvqPOBfgCdE5Ifec97mnU7INIO6gQ3VNR4b2GCMaTUS6ntcRHyb1wT4Ds69QasAVPXRFi9dM8jLy9P8/KhuYzJYH5ExpmlEpEBV8xpzTLgaUeDUPm+EWF/35hep6iIR+Waw7ar6RrD1JrE0ZmCDMcY0h5BBpKo/a+S5puCMlvtGsNPREGSmGVkNxhjT1kX7hNZ/B0biDA4AQFUvCtg1w7v+9uYsoAnNRrkZY5JBNKPmXgQ2AkOAnwHbcEa4BZrWfMUy0bBRbsaYZBDNfUQ9VPVvInK/qv4D+IeIBAsit4hk4gxsOInNvt38bJSbMSYZhAwiEclR1e1AtXdVmYhciTO/W7D7gk7DmYEhWBBFHL5tGs+m7zHGJINwNaK3gHHAL0SkG849PH/AmXT0h0H2X++dFdu0IhvlZoxp68IFkQCo6rve5YPAhS1eImOMMSklXBD1F5Hfh9oYZDbtx5qnSMYYY1JJuCA6htPnExVVfabJpTHGGJNywgVRuao+22olMcYYk5LC3UdUFcsJReScaNaZyAqKK3j8kyIKiiviXRRjjGkx4ab4mRTjOf+AM9ou0joThs2aYIxJFdHc0BoV7xNTJwM9A2bu7orzWAfTCMFmTbAgMsYko2im+IlWBtAZJ9y6+Px3CLgu0sEiMk1ECkWkSER+EmKfG0RkvYisE5GXmrHsCadu1gS3kFCzJlhzoTGmuYV8HpHfTiLnAkNV9WnvJKidVXVriH0HqWpxowoh4gY2AVOBEpy57Gao6nqffYYCrwIXqWqFiPRS1YiPCG/LzyOKy8zaO76AbYth8HkwcMJJ5bHmQmNMOM39PKK6k/4UyAOGA08D6cALQKgBCJUi8msiz9btawJQpKpbvO85B5gOrPfZ5y7gcVWt8J4vYgi1da0+a8KOL+DZq6G2CtwZcOs8vzCy5kJjTEuIpmnuWuBq4CiAqu4kxMPxvKKdrdtXf2CHz3KJd52vYcAwEVkiIstEJORs3yIyS0TyRSR/7969Ed7a1Nu22AkhrXV+blvstzlRmwuNMW1bNIMVqlRVRUQBRKRThP2jna27sdKAocAFwADgUxE5Q1UPBO6oqrOB2eA0zTXDe6eGwec5NaG6GtHg8/w22ySrxpiWEE0QvSoifwW6i8hdwB3Ak2H2j3a2bl+lwECf5QHedb5KgOWqWg1sFZFNOMHUHCEXNwn1hNWBE5zmuBB9RGCTrBpjml/EIFLV34jIVJzRb8OBB1V1QZhDop2t29cKYKiIDMEJoJuAmQH7vAXMAJ4WkWycprotkcqfyBKy83/ghKABZIwxLSWawQoPAM/4ho+IzPI2f50kltm6VbVGRO4D5uPcc/SUqq4TkYeBfFWd5912qYisB2qBf1PVNv1IUuv8N8aY6Jrmvg/cJCL3qeon3nV34+2DCeQd3n0XMNj3/Kp6R7g3UdX3gfcD1j3o81qBH3n/Swr2hFVjjIkuiEpxhlLPFZHXVPXXhHgcuNfbwGLgY5yaiwnBOv+NMSbKKX5UdbuInA/8WUTmAh3C7N5RVX/cLKVLEuEGJFjnvzEm1UUTRPkAqnocuF1E7gXGh9n/XRG5wtvUlvISckCCMcYkkIg3tKrqXQHLj6tqbuB+InJYRA4B9+OE0TEROeSzPiUFG5BgjDGmQcgakYi8qqo3iMga4KSbQlV1dMByuNkWUlYqDUhIqHuijDFtRrimufu9P69qzAlFZKGqXhxpXapIugEJISZFtSZIY0yswj0Yr8z7sxhARHoAU4DtqloQuL+ItAc6AdkikknDyLqunDxvXNJJiQEJYSZFtXuijDGxCtc09y7wE1VdKyJ9gZU4AxdOEZHZqvq7gEO+B/wA6AcU0BBEh4A/NnO5E0rK1AaCTYrqDaJUaoI0xjSvcE1zQ1R1rff17cACVf2OiHQBlgC/891ZVR8DHhOR76vqH1qktAkqZWoDYSZFTbomSGNMqwkXRNU+ry8GngBQ1cMi4gl1UFKFUJj+EN8v3JSpDUSYFDVpmiCNMa0qXBDtEJHv48x6PQ74EEBEOuA8HC+5hegPCdUMlzK1AZsU1Zg2L9FGuIYLou8CDwOXADf6PPdnEs6TWpNbiP6QUM1wVhswxrQFidinHW7U3B6cyU0D138CfHLyEQ1EZDQnT3r6RsyljIcQ/SEp0wxnjElKidinHdVcc40hIk8Bo4F1QF1fkgIRg8j7+O/HcB4F8aSqPhJiv28BrwFnqWp+c5T7JCH6Q1KqGa6xgvSpJVoTgDGppi30aYvzdIVmPKHIelU9PYbj3MAmYCpOv9QKYIaqrg/YrwvwHpAB3BcpiPLy8jQ/v2WyyvgI0qdW4BmacE0AxqSSUM1wLfkHoogUqGpeY46JONdcDD4XkUYHETABKFLVLapaBczBefxEoJ8DvwSON6GMprkF6VOzefaMaV0FxRU8/kkRBcUVQOi5LscPyuTeC09NmD8Mo3lCay5Oc9nZOE1tnwM/VNVQj+l+DieMdgEncG5s1cC56YLoD+zwWS4BJgaUZRwwUFXfE5F/i1R204qC9KlN8iReE4AxySpY7ScRm+GCiaaP6CXgceBa7/JNwMsEhISPvwHfBtbQ0EfUZCLiAh4Fboti31nALICcnJzmKoIJJ0if2niw/jRjWkmw2s+9F57aJv4NRhNEHVX1eZ/lFyLURvaq6rwYylIKDPRZHuBdV6cLMAr4u4gA9AHmicjVgf1Eqjob76PM8/LymrcTzIQW5B4jG9ZuTMuIdhBCW/g3GE0QfSAiP8Hps1HgRuB9EckCUNX9AfuvEpGXgHdwmubw7hdp1NwKYKiIDMEJoJuAmT7HHwSy65ZF5O/Av7bYqDnTPGwknTFNFvhvJtlurI8miG7w/vxewPqbcIIp8CF5HXAC6FKfdRGHb6tqjYjcB8zHGb79lKquE5GHgfwYa1kmnmwknTFNFix0ku3G+ohBpKpDGnNCVb091sJ4Hy/+fsC6B0Pse0Gs72NaSbCRdDVZCXcznTGJJLD2Eyx02soghGhFM2quI/AjIEdVZ4nIUGC4qr4bYv8BwB+Ac7yrFgP3q2pJM5XZtBU2ks6YiHyDB4hq5FtbboYLJuINrSLyCs7zhb6jqqO8wbRUVceE2H8Bzki7ugEOtwA3q+rUZit1I9gNrXEWRR+R9RmZVBXY7PatcQN4+YvteBTcAj+6dDj3Xnhqm/o3EssNrdH0EZ2iqjeKyAwAVa0U77C1EHqqqu+kqM+IyA8aUyiTRIKNpHN9zfi0xeA6j4Ji6zMyqcU3VAKb3RTa7Mi3pogmiKq8j35QABE5BZ/RcEGUi8gtOPcaAcwA7JZ64wgYwLBs1BzrM2oBVutMHOGa3h68aqRf8Hxr3AC+NW5Ayv2/iiaIHsJ5FtFAEXkRp+8n3ICEO3D6iP4PJ7yWRtjfpJKAAQyTXOvISMu1PqMmivRl9/C76/xqnUDKfdm1lnD/L741boDfH14VlVVB+3pS7f9JNKPmPhKRApznEAnOwIN9wfb1Tlz6P6p6dfMW0ySNgAEM48eexVt9D1GxfhGZp1/EafbXe6MF62fw/bL7YG2Z3/LrK0t4Y2WJNYfGKFxtE8IHT7Cmt2RvdotGNKPmFqrqxTgzXgeu86OqtSIySEQyvBOXGuMvcCog4LT5tzjBtOMpCnidm+cdsi/JCBrTz3D5qL6s2La/flkg6ESYFv7BNaa2GSl4UrXpLZKQQSQi7YGOQLaIZOLUhgC64kxQGsoWYImIzAOO1q1U1UebXlyTFHwHMCz+rV9T3bK1m6iq6WV9RmEE1oCi6WcY3qeL35fp6ytL6vfP7Jhx0vkqKqtS6osyMGyireEE1jajDZ5Uua7RClcj+h7wA6AfzvDtuiA6BPwxzHGbvf+5cOaHMya0gKa6SaOGMWFLPuN1HQUykkm5k62pjvA1oGj6GQKbf3z39z1fVbWHB99ei0c16fqTogmbNJeACDW10dVwAmubFjyxCfeo8MeAx0Tk+6r6h0gnEpHnVfXbwAHvscZEFtBUNx54KeN/vMH0Nqv25KV8U12kGlAs/QyB+9edT0TwqIbtTwKCfqHH6/9LqIDx7b/J7JhR34QWNmxqFVCU6Gs4vrVNC57YhGuaOwvYURdCIvId4FtAMfBQkMlOx4tIP+AOEXmOhhoUEHRyVGMcAU11Lk814AFPNcvWbmJk7QEmujbwRe0IXl+ZOu3rdV+iOw8ci6oGFCvfu/TrvrBD9Sf5BlPgF3pgk164Gkik4Ij2GAhdm/Htv3H5BmyYsHF7z1FbG30NxwYbNF24prm/ApcAiMgU4BHg+8AYnEcsXBew/1+AhTiToPo25UHwyVGNOVlAU93F/Wu5Y9v/kE4NX3iGc3v+MKo9JH3tyLcWlOYS0tyu+i/Hlhhp5Xu+cP1JfsHk84Ue2KTnGwK+ARHqdazHhKvN+PbfoIrLJQgaMWzg5Fpesn7OEkW4IHL71GJuBGar6uvA6yKyOnBnVf098HsR+bOq3tP8RTUpIaCp7rRti1FXDaIeVulwRnk2JXXtKFgtqNaj3DhhIP27d2iV3zVcfxI0BJPvF3pgk55fJ75vQIR6HeMx4Wozgf03vjU2CB82yfJ5aivCBpGIpKlqDXAx3ieeRjrOQsg0WcC0QOJuB7VVXJz2FXe436+vHT1WAJls4NeLRnL1Vde22ZFeofoxfGtB3xo3IG6/V6RgCtak5xsCvgER6nWsx0SqzQTrv/H9vUxiCDnpqYj8F3AFsA/IAcapqorIqcCzqnpO0ANjLYjINOAxnGcRPamqjwRs/xFwJ1AD7AXuUNXiSOe1SU+TQN3EqQdL0IJnEPXwWs0UrnQvqw+lW2v/H7UqbW74sW8TnG8/hlvgxgk5rVYLag6RbvRsiT6itnBdUk0sk56GnX1bRCYBfYGPVPWod90woLOqrmxKYQPexw1sAqYCJThPa52hqut99rkQWO6ddPUe4AJVvTHSuS2IkojPPHWlmk1v3UuaePhD9XT+r/YGPAguwOUSv76KRAwl3ya4utmW68quqqQneR+YSV7NPvu2qi4Lsm5TYwsWhQlAkapuARCROcB0oD6IVPUTn/2X4TxewqQSn/6j/h164Pngx3hqqzk7vYjHRaiuxa+vwrcDPc0lXJ83kG/GoYkr8K/811eW8FpBSX3nu28TXKIGpzEtKZpJT1tDf2CHz3IJMDHM/t8FPgi1UURm4e3TysnJaY7ymUTh03/k6n06bFtM3uDzeNEz9KS+Cr9QqlVeWr6d11eWROy0bi4FxRUnhQ4i9Z3s0PoDEYxJRIkSRFHzPmIiDzg/1D6qOhtniDl5eXnhn/xn2i6fUBoPJw0/rgulE9XOF3/gMOPAYcHX5w1kZL9ujQop34EGFZVVfj993xv8R3yBc39DvAciGJMIEiWISoGBPssDvOv8iMglwH8B56tquGcimRQWeE9MXa3kpGHGvvfB1CovLt8OgEsIG1KhwkZwws4l1A888Asdt//or3g1FRoTSrym04r4qPBWKYRIGs5ghYtxAmgFMFNV1/nsMxZ4DZimql9He24brGDg5CHSvkOEfZvK6tTdjR1sfWDYeIL8E/IdeOAbOmAjvkxiCpxKKtbBMi31qPAWp6o1InIfMB9n+PZTqrpORB4G8lV1HvBroDMw1/uk8u323CMTrUgzB7xWUEJNjQcP/jWiwJCqe+17tz6qePAPqXAj9iyATCIKnEy3NWe+T4ggAlDV94H3A9Y96PP6klYvlElKgTdojh+UWX9TZF2zW6iQChU2gX1EVuMxbc2k3B4nTabbWhKiaa4lWdOcaapQAxIsbEyyaY4+ojbbNGdMIrPZlU2qiNdn3dXq72iMMcb4SPqmORHZi/MMpVhl48y3l+rsOjSwa9HAroXDrkOD4araqKdzJ33TnKr2bMrxIpLf2PbOZGTXoYFdiwZ2LRx2HRqISKM75a1pzhhjTFxZEBljjIkrC6LIZse7AAnCrkMDuxYN7Fo47Do0aPS1SPrBCsYYYxKb1YiMMcbElQWRMcaYuLIgCkFEpolIoYgUichP4l2eeBKRbSKyRkRWxzI0sy0TkadEZI+IrPVZlyUiC0Tka+/PpJ92IcR1eEhESr2fi9UickU8y9haRGSgiHwiIutFZJ2I3O9dn1KfizDXodGfC+sjCkJE3DiPpZiK87TYFcAMVV0f9sAkJSLbgDxVTbkb9kRkCnAEeE5VR3nX/QrYr6qPeP9IyVTVH8eznC0txHV4CDiiqr+JZ9lam4j0Bfqq6koR6QIUANcAt5FCn4sw1+EGGvm5sBpRcBOAIlXdoqpVwBxgepzLZOJAVT8F9gesng486339LM4/vqQW4jqkJFUtU9WV3teHgQ1Af1LscxHmOjSaBVFw/YEdPsslxHiBk4QCH4lIgYjMindhEkBvVS3zvt4F9I5nYeLsPhH5ytt0l9RNUcGIyGBgLLCcFP5cBFwHaOTnwoLIRONcVR0HXA7c622mMYA6bdup2r79Z+AUYAxQBvw2rqVpZSLSGXgd+IGqHvLdlkqfiyDXodGfCwui4EqBgT7LA7zrUpKqlnp/7gHexGm6TGW7ve3jde3ke+JcnrhQ1d2qWquqHuAJUuhzISLpOF++L6rqG97VKfe5CHYdYvlcWBAFtwIYKiJDRCQDuAmYF+cyxYWIdPJ2RCIinYBLgbXhj0p684Bbva9vBd6OY1nipu5L1+taUuRzISIC/A3YoKqP+mxKqc9FqOsQy+fCRs2F4B1y+DvADTylqv8d3xLFh4jk4tSCwJmt/aVUuhYi8jJwAc40/7uBnwJvAa8COTiPGLlBVZO6Iz/EdbgAp/lFgW3A93z6SJKWiJwLLAbWAB7v6v/E6R9Jmc9FmOswg0Z+LiyIjDHGxJU1zRljjIkrCyJjjDFxZUFkjDEmriyIjDHGxJUFkTHGmLiyIDKmFYhID5/ZiHf5zE58RET+FO/yGRNPNnzbmFaWqrNWGxOK1YiMiSMRuUBE3vW+fkhEnhWRxSJSLCLfFJFfeZ8F9aF3OhVEZLyI/MM7Ce38gDvZjWlzLIiMSSynABcBVwMvAJ+o6hnAMeBKbxj9AbhOVccDTwEpM9OFSU5p8S6AMcbPB6paLSJrcKaX+tC7fg0wGBgOjAIWOFN94caZ4diYNsuCyJjEcgJAVT0iUq0NnbgenH+vAqxT1bPjVUBjmps1zRnTthQCPUXkbHCm4ReRkXEukzFNYkFkTBvifXT9dcAvReRLYDUwOa6FMqaJbPi2McaYuLIakTHGmLiyIDLGGBNXFkTGGGPiyoLIGGNMXFkQGWOMiSsLImOMMXFlQWSMMSau/j8GOl1h6pY6GwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "def goodwin_oscillator( t, statevars, params ):\n",
    "    X = statevars[0];\n",
    "    Y = statevars[1];\n",
    "    Z = statevars[2];\n",
    "    alpha_1 = params[0]\n",
    "    K = params[1]\n",
    "    n = params[2]\n",
    "    gamma_1 = params[3]\n",
    "    alpha_2 = params[4]\n",
    "    gamma_2 = params[5]\n",
    "    alpha_3 = params[6]\n",
    "    gamma_3 = params[7]\n",
    "        \n",
    "    dydt = np.zeros( (3,) )\n",
    "    dydt[0] = alpha_1*K**n/(K**n+Z**n) - gamma_1*X\n",
    "    dydt[1] = alpha_2*X - gamma_2*Y;\n",
    "    dydt[2] = alpha_3*Y - gamma_3*Z;\n",
    "\n",
    "    return dydt\n",
    "\n",
    "f = goodwin_oscillator\n",
    "alpha_1 = 1\n",
    "K = 0.61\n",
    "n = 11\n",
    "gamma_1 = 0.15\n",
    "alpha_2 = 1\n",
    "gamma_2 = 0.15\n",
    "alpha_3 = 1\n",
    "gamma_3 = 0.25\n",
    "params = (alpha_1,K,n,gamma_1,alpha_2,gamma_2,alpha_3,gamma_3)\n",
    "\n",
    "y0 = np.array([0.02051344, 0.27934056, 1.31619671])\n",
    "\n",
    "t0 = 0\n",
    "tFinal = 24\n",
    "h = 1\n",
    "ETOL = 1e-3\n",
    "t = [t0];\n",
    "steps = [];\n",
    "# start with y as a list of nparrays.\n",
    "y = [np.array(y0)]; # list of 1st order solutions\n",
    "attempted_t = [t0,[]]\n",
    "attempted_y = [np.array(y0),[]]\n",
    "while t[-1] < tFinal:\n",
    "    Yi = y[-1];\n",
    "    YFE = Yi + h * f(t[0],Yi,params);\n",
    "    YET = Yi + h/2 * f(t[-1],Yi,params) + h/2 * f(t[-1]+h,YFE,params);\n",
    "    attempted_t[-1].append( t[-1]+h )\n",
    "    attempted_y[-1].append( YFE )\n",
    "    err = numpy.linalg.norm(YET-YFE,2);\n",
    "    hnew = h * np.sqrt(0.9*ETOL/err);\n",
    "    # Make sure h doesn't more than double\n",
    "    if hnew / h >= 2:\n",
    "        hnew = 2*h;\n",
    "    if err <= ETOL:\n",
    "        # The step size was small enough.  Record this\n",
    "        # info and determine the next step size.\n",
    "        t.append( t[-1] + h );\n",
    "        steps.append( h );\n",
    "        y.append( YFE )\n",
    "        h = hnew;\n",
    "        if t[-1] + h >= tFinal:\n",
    "            # If we are going to step out of bounds, then\n",
    "            # restrict the size of the last step.\n",
    "            h = tFinal - t[-1];\n",
    "        attempted_t.append( [] )\n",
    "        attempted_y.append( [] )\n",
    "    else:\n",
    "        h = hnew;\n",
    "\n",
    "# The first step taken from any time point will be blue (X), orange (Y), or green (Z)\n",
    "# The second step taken will be red(X), purple (Y), or brown (Z)\n",
    "plt.subplot( 2, 1, 1 )\n",
    "for i in range(len(attempted_t)-1):\n",
    "    plt.gca().set_prop_cycle(None) # this resets the color cycle\n",
    "    for j in range(len(attempted_t[i+1])):\n",
    "        plt.plot( [t[i], attempted_t[i+1][j]], [y[i], attempted_y[i+1][j]], '.-' )\n",
    "plt.ylabel( \"Y\")\n",
    "plt.subplot( 2, 1, 2 )\n",
    "for i in range(len(attempted_t)-1):\n",
    "    plt.gca().set_prop_cycle(None) # this resets the color cycle\n",
    "    for j in range(len(attempted_t[i+1])):\n",
    "        plt.plot( t[i], attempted_t[i+1][j]-t[i], '.' )\n",
    "plt.xlabel( 'Time')\n",
    "plt.ylabel( \"Step Size Taken \\nfrom that Time\");\n",
    "plt.tight_layout()"
   ]
  }
 ],
 "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
}