P02-The-Heun-Method.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Project 2: The Heun Method and Coupled Equations\n",
    "===========================\n",
    "\n",
    "Two Learning Goals\n",
    "------------------\n",
    "\n",
    "This project involves two related learning goals 1) Begin to learn a slightly improved method of integrating a differential equation (The Heun Method) and 2) learn how to turn higher order equations into a set of first order coupled equations along with a simple technique to manage the extra complexity.\n",
    "\n",
    "The Heun Method\n",
    "----------------\n",
    "\n",
    "While Euler's Method is quite easy and accurate enough for simple problems\n",
    "it is less successful when applied to more challenging situations. Depending\n",
    "on the nature of the function $f(s,t)$ Euler's Method can  easily\n",
    "under or overestimate the change. Getting satisfactory results using\n",
    "the Euler Method can require very small time steps, making the algorithm\n",
    "very slow. \n",
    "\n",
    "There are many ways to improve the situation and we'll see\n",
    "a number of them in this course. Our first improvement is called\n",
    "\"The Heun Method\".  The idea is to take a \"fake\" time step, re-evaluate\n",
    "the function $f(s,t)$ there, and then average the two results\n",
    "for the \"real\" timestep. This process is animated in the online slides, and I ecourage\n",
    "you to watch those if you can. (Check under Resources -> Slides on ACE).\n",
    "\n",
    "Starting at $x_1$ find the slope $f(s,t)$ and\n",
    "extrapolate out to $x_1 + \\Delta x$.\n",
    "\n",
    "<img width=\"250\" src=\"https://raw.githubusercontent.com/sspickle/sci-comp-notebooks/master/imgs/heun1.png\">\n",
    "\n",
    "Next.. evaluate the slope at the extrapolated value of s.\n",
    "\n",
    "<img width=\"250\" src=\"https://raw.githubusercontent.com/sspickle/sci-comp-notebooks/master/imgs/heun2.png\">\n",
    "\n",
    "Using this slope.. extrapolate *again*. Note that while the\n",
    "first extrapolation underestimated the change in s, the second\n",
    "overestimated. Heun's method is to use the *average* of these \n",
    "two estimates:\n",
    "\n",
    "<img width=\"250\" src=\"https://raw.githubusercontent.com/sspickle/sci-comp-notebooks/master/imgs/heun3.png\">\n",
    "\n",
    "You can see that the average is much closer to the exact\n",
    "solution than the original Euler Method result. In short\n",
    "here's how it works:\n",
    "\n",
    "Starting at $x=x_1$, take a fake Euler step,\n",
    " \n",
    "first evaluate f(s,t):  \n",
    "\n",
    "$$f_1 = f(s_1, x_1)$$\n",
    " \n",
    "next jump to the next x: \n",
    "\n",
    "$$x_2 = x_1 + \\Delta x$$\n",
    " \n",
    "estimate the new value of s: \n",
    "\n",
    "$$s_2 = s_1 + f(s_1,x_1) \\Delta x$$\n",
    "  \n",
    "Evaluate f(s,t) at the new position:\n",
    "    \n",
    "$$f_2 = f(s_2, x_2)$$\n",
    "            \n",
    "Extrapolate using the average of the 'f' values:\n",
    "    \n",
    "$$s_3 = s_1 + \\frac{1}{2}(f_1 + f_2) \\Delta x$$\n",
    "            \n",
    "Now $s_3$ is the Heun's Method result.\n",
    "\n",
    "Let's re-work our project from last time using both the Euler and Heun methods to compare."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fce1adecc70>,\n",
       " <matplotlib.lines.Line2D at 0x7fce1adecd60>]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXhU5fn/8fedhIAsIhRkMywKqCiIitZUhCAuqEXUqmjtD7BarH7VtrZVEHFDRetWrbW4Vm3d0FZFxaosEZeggIqCgoKgKDuIgiIhyf3745xgCAmZSWbmTJLP67rOxZwzz5y5zyTMnWc5z2PujoiISKwyog5ARERqFyUOERGJixKHiIjERYlDRETiosQhIiJxUeIQEZG4KHFI2jCzfDM7N3w8wszeKPPcJjPbM7rods7MdjGz583sGzN7Kjx2nZmtNbOVZtYxvIbMqGONhZmdZ2Z/TdC53My6VvO1vczsrUTEIYmjxCFJYWZLzWxz+GVZurWv7vncvam7f5bIGAHM7Orwi+135Y7/Ljx+dYynOhVoA/zE3U8zs47AH4Ee7t7W3b8Ir6G4mnHuamZ/NbMvws9ycbjfqjrnq+K9soErgJvD/c7hZ7Gp3DY00e9dnrt/AGwws8HJfi+JnRKHJNPg8MuydFsedUCV+AQYVu7Y8PB4rDoBn7h7UbjfEVjn7qtrGlz4RT4V2A8YBOwK5ALrgENrev4KDAEWuPtX5Y7vVu7n+WQS3nsbM8sKHz4KnJfM95L4KHFIyphZCzN7wczWmNnX4eM9YnzttuYOM3vIzP5uZi+a2UYze9vM9ipT9hgzWxg2G91tZq+VNoFVYhbQ2Mz2C1+/H9AoPF56zu2azsrGZGbXAFcCQ8O/xM8DXgXah/sPlfmrPSt8bb6ZjTOzN8NreGUntYdhBInoZHf/yN1L3H21u49z98nlP58yn9F14eM8M/vSzC4Pm86WmtlZO/k8jgNe28nz2ynbxFjZZ1XmuYZmdktYc1plZhPMbJdycV5mZiuBf4YvywcGmlnDWGOS5FLikFTKIPgy6ETwRbgZuKua5zoDuAZoASwCrgcIv3yfBkYDPwEWAj+L4Xz/4sdax/BwPybufhVwA/Bk+Jf4PQRfvsvD/RGVvPSXwNnA7kA28KdKyh0F/M/dN8UaUwXaAq2ADgTXd6+Z7V1J2Z4En1sy3Ah0B3oDXcN4riwXZ0uC35GRAGHNZytQWbySYkockkzPmtmGcHvW3de5+3/c/Xt330jwZd+/mud+xt3fCZuGHiX4IgI4Hpjv7v8Nn7sTWBnD+f4NnGlmDQiS0r+rGVc8/unun7j7ZmAiP15DeT8BViTg/ca6+xZ3fw14ETi9knK7ARsrOL62zM9zg5ntG8+bm5kRJIM/uPv68HfgBoLPu1QJcFUY5+YyxzeGcUkayKq6iEi1neTuU0p3zKwxcDtBO32L8HAzM8usRqdx2WTwPdA0fNweWFb6hLu7mX1Z1cnc/QszW0TwRfapuy8LvueSqrJrKG8d0K6G7/W1u39XZv9zgs+qwrJAswqOtyrTh1MdrYHGwJwyn60BZUearXH3Hyp4bTNgQw3eWxJINQ5JpT8SNDf81N13BfqFxxP5Db0C2NZvEv6VG1M/CvBIGOMjFTz3HcGXXul529YgxnhNAY41syY7KfM9ZeIjaPIpq0W513cEKhus8AFBc1KstvtsKnjvUmsJmif3c/fdwq25u5dNmDtM121mHQia8pLVfCZxUuKQVGpG8MWxwcxaAlcl4T1eBHqa2UlhR/T/UfkXWXlPAscQNBuVNxfYz8x6m1kj4OpEBBujfxHUov5jZvuYWYaZ/STs7D4+LPM+8EszyzSzQVTcBHiNmWWb2RHAz4GnKnm/yZW8vjLvA6eYWeOwg/6cigq5ewlwH3C7me0OQVIws2OrOH9/YJq7b4kjJkkiJQ5Jpb8CuxD85TkT+F+i38Dd1wKnAX8haOLpAcwGqvzScffN7j6lXNt66XOfANcS/PX/KVDhqKFkCL8wjwIWEIzW+hZ4h6Cz++2w2O+AwQTNOWcBz5Y7zUqCJqjlBH1Cv3X3BZW85fPAPrbjfTcbyt3HcUl4/HagEFgFPByevzKXEQxmmGlm3xJ8nlV1ep8FTKiijKSQaSEnqcvMLAP4EjjL3adHHU8UzCwP+Le7x9pkh5mNJLh58fdJCyy2OHoB97h7bpRxyPbUOS51Ttj08TZBs9ifCfpQZkYaVC3j7vdGHQNsu3NcSSPNqKlK6qJcYDFBk9hggtFdOzQ/iUj1qKlKRETiohqHiIjEpV70cbRq1co7d+4cdRgiIrXKnDlz1rp76/LH60Xi6Ny5M7Nnz446DBGRWsXMPq/ouJqqREQkLkocIiISFyUOERGJixKHiIjERYlDRETiknaJw8wGhct+LjKzURU839DMngyff9vMOqc+ShGR+iutEoeZZQJ/J1h2swfBimw9yhU7h2BRmq4Es3LelKx4CpYVMP718RQsK0jWW4iI1DpplTiAQ4FF7v6ZuxcCTwBDypUZQjB1MwRrSw+0JCzVVrCsgLxrrmTM1d8z8JGBSh4iIqF0SxwdKLPsJ8F02B0qKxMuY/kNwZrM2zGzkWY228xmr1mzJu5A8pfms/WTI/Fp17Bl5Z7kL82P+xwiInVRuiWOhHH3e929j7v3ad16hzvmq5TXOY+Gfe+GBpux1y8nr3Ne4oMUEamF0i1xfAXklNnfIzxWYZlwadDmBCu9JVRuTi7Tzp9I31Pm4R+eye6FWhJARATSL3HMArqZWRczywbOACaVKzMJGB4+PpVgLeKkzA2fm5PLxNt/SoMGxo03JuMdRERqn7RKHGGfxYXAy8DHwER3n29m15rZiWGxB4CfmNki4BJghyG7idSuHfzmN/Dww/DFF8l8JxGR2qFeLOTUp08fr8nsuF98AV27wsiRcNddCQxMRCSNmdkcd+9T/nha1TjSVceOMHw43H8/rFgRdTQiItFS4ojR6NFQVAS33BJ1JCIi0VLiiNGee8IvfwkTJkA1bgsREakzlDjicPnlsHkz3H571JGIiERHiSMO++wDp50WdJB//XXU0YiIREOJI05jxsDGjXDnnVFHIiISDSWOOPXqBUOGwB13wLffRh2NiEjqKXFUwxVXBE1Vd98ddSQiIqmnxFENffrAoEFw663w3XdRRyMiklpKHNV0xRWwdi3cd1/UkYiIpJYSRzUdfjjk5cFf/gI//BB1NCIiqaPEUQNjxwZTkAwd8z+tECgi9YYSRw006lpARs5MJj2wL0c+OEjJQ0TqBSWOGnjt83zoPw6+6cSWd4dqeVkRqReUOGogr3Me2XtPhw6zYMbl9O0wIOqQRESSTomjBnJzcpk2fCrDLl6Cb+jMotcOizokEZGk00JOCeAOhxwCGzbAggWQlZW0txIRSRkt5JREZnDllbB4MTz2WNTRiIgklxJHggweDL17w3XXBQs+iYjUVUocCVJa6/j0U3jyyaijERFJHiWOBBoyBHr2hHHjoLg46mhERJJDiSOBMjKCWsfChfDUU1FHIyKSHEocCXbKKbDffkGto6Qk6mhERBIvbRKHmbU0s1fN7NPw3xaVlCs2s/fDbVKq46xKRkYwh9VHH8F//hN1NCIiiZc2iQMYBUx1927A1HC/IpvdvXe4nZi68GJ36qnB+uTXXqtah4jUPemUOIYAD4ePHwZOijCWGsnMDGod8+bBs89GHY2ISGKlU+Jo4+4rwscrgTaVlGtkZrPNbKaZVZpczGxkWG72mjVrEh5sVYYOhe7dVesQkbonpYnDzKaY2bwKtiFly3kwD0plc6F0Cm+B/yXwVzPbq6JC7n6vu/dx9z6tW7dO7IXEIDMzWCVw7lx4/vmUv72ISNKkzVxVZrYQyHP3FWbWDsh3972reM1DwAvu/vTOyiV7rqrKFBXBvvtCs2YwZ05wk6CISG1RG+aqmgQMDx8PB54rX8DMWphZw/BxK+Bw4KOURRinrCwYMwbeew9G/OUpLfQkInVCOiWOG4GjzexT4KhwHzPrY2b3h2X2BWab2VxgOnCju6dt4gDYq99MrMUSHrmzM0c+PFDJQ0RqvbSZANzd1wEDKzg+Gzg3fPwW0DPFodXIG19Nh36L4bn7KfzoGPKX5pObkxt1WCIi1ZZONY46Ka9zHg0PehJafopPv5Z+HfOiDklEpEaUOJIsNyeXaWe/wmkXfIyv7MWKWaptiEjtljajqpIpqlFVZRUXBzPnmsEHHwTDdUVE0lltGFVVp2VmwtVXB3NYab0OEanNlDhS6NRToVevIIFolUARqa2UOFIoIwOuuSZYJfDf/446GhGR6lHiSLEhQ+Dgg4M5rLZujToaEZH4KXGkmFmQNJYsgX/+M+poRETip8QRgeOOg8MOg+uugy1boo5GRCQ+ShwRMAuWll22DO67L+poRETio8QRkYEDoV8/uOEG2Lw56mhERGKnxBGR0r6OFSvgH/+IOhoRkdgpcUSof/+g5nHjjbBpU9TRiIjERokjYuPGwZo1cNddUUciIhIbJY6I5eYGo6xuvhm+/TbqaEREqqbEkQauvRbWr4dT/jBDCz2JSNpT4kgDW9sUkLHvJKY+2psj/3G6koeIpDUljjSQvzQfjhwDhU3Ykv/HYF9EJE0pcaSBvM55NGy3GOv9MP7O+eydeWzUIYmIVEqJIw3k5uQyddhU/nz592RnNeDZCQdFHZKISKWyog5AArk5ueTm5MKsYITVn/4UrN0hIpJuVONIM6NGQfPmMHp01JGIiFQsbRKHmZ1mZvPNrMTMdljjtky5QWa20MwWmdmoVMaYCi1aBElj8mTIz486GhGRHaVN4gDmAacAMyorYGaZwN+B44AewJlm1iM14aXORRdBhw5w2WXgHnU0IiLbS5vE4e4fu/vCKoodCixy98/cvRB4AhiS/OhSa5ddgpsC33kH/vvfqKMREdle2iSOGHUAlpXZ/zI8tgMzG2lms81s9po1a1ISXCINGwY9esDll0NRUdTRiIj8KKWJw8ymmNm8CraE1xrc/V537+PufVq3bp3o0yddVhaMHw+ffAIPPhh1NCIiP0rpcFx3P6qGp/gKyCmzv0d4rE4aPBgOPxyuvhrOOguaNIk6IhGR2tdUNQvoZmZdzCwbOAOYFHFMSWMWrNWxYgXccUfU0YiIBNImcZjZyWb2JZALvGhmL4fH25vZZAB3LwIuBF4GPgYmuvv8qGJOhb594cQT4aabYN26qKMREQHzejDes0+fPj579uyow6i2+fODu8h//3u49daooxGR+sLM5rj7DvfVpU2NQyq3334wfHiwSuDnn0cdjYjUd0octcQ11wBWwsnnfaj1OkQkUkoctcSXFFBy6G2893JPBtxwmZKHiERGiaOWyF+aT8kR10GTlWx58SamL8mPOiQRqaeUOGqJvM55NGxciB11BSzLZevcU6MOSUTqKSWOWqJ0sadxf9iLfXpu4v6/dOO776KOSkTqIyWOWiQ3J5cx/Udz391N+fJL+Mtfoo5IROojJY5aqG9fGDo0SBxffBF1NCJS3yhx1FKltY3LLos2DhGpf5Q4aqmOHeHSS+GJJ+CNN6KORkTqEyWOWuzSS4OVAn//eygpiToaEakvlDhqsSZNgskP58yBhx+OOhoRqS/iThxm1iRc+1vSwC9/CYcdBqNHw7ffRh2NiNQHVSYOM8sws1+a2YtmthpYAKwws4/M7GYz65r8MKUyZsFaHatWwQ03RB2NiNQHsdQ4pgN7AaOBtu6e4+67A32BmcBNZvarJMYoVTj00GCN8ttvh8WLo45GROq6KtfjMLMG7r61pmWiVNvX44jF8uXQvTsccwz8979RRyMidUG11+MoTQhmdpqZNQsfjzWz/5rZQWXLSHTatw/6OZ55BqZNizoaEanL4ukcH+vuG82sLzAQeAD4R3LCkuq45BJol/MDZ56zhhmLZ0YdjojUUfEkjuLw3xOAe939RSA78SFJdb2/toB1/c9i9dLWDDx/stbsEJGkiCdxfGVm9wBDgclm1jDO10uS5S/Np7jbc7Dv0xRNH81/3nw/6pBEpA6KZThurpkZcDrwMnCsu28AWgJ/TnJ8Eoe8znlkZ2aTcfwlkLmVNyacSRVjH0RE4pYVQ5lhwN+BT4D/ARsB3H0FsCJ5oUm8StfsyF+az5pd13H7VV148kk444yoIxORuqTK4bjbCprtAxwHHAs0J7i/43/Am+5evLPXxnj+04CrgX2BQ929wvGzZraUIHkVA0UVDRUrrz4Mxy2vuDi4o3zZMliwAHbbLeqIRKS2qfZw3FLuvsDdb3f3QcCRwBvAacDbCYpxHnAKMCOGsgPcvXcsSaO+ysyEe+6BNWuCYboiIolSrc5td9/s7pPd/aJEfXm7+8fuvjAR55LAQQfB734HEyZAgQZYiUiCxJw4zKyPmT1jZu+a2QelWzKDq4QDr5jZHDMbGcH71yrXXgs5OXDeebBVt2mKSALE0jle6lGCUVQfAtVa/cHMpgBtK3hqjLs/F+Np+rr7V2a2O/CqmS1w9x2at8KkMhKgY8eO1Qm3TmjaFP72NzjppGAuq0svjToiEant4ukcf8Pd+yY5HswsH/hTZZ3j5cpeDWxy91t2Vq4+do6Xd/LJ8PLLMH8+dOkSdTQiUhvUuHMcuMrM7jezM83slNItgTFWKVwLpHS+rCbAMQSd6lKFO+8MOswvvBDd2yEiNRJP4jgb6A0MAgaH288TFYiZnWxmXwK5wItm9nJ4vL2ZTQ6LtQHeMLO5wDvAi+7+v0TFUJfl5MC4cTB5Mjz9dNTRiEhtFk9T1UJ33zvJ8SSFmqoCRUXw058GU7AvWADNm0cdkYiks0Q0Vb1lZj0SGJOkWFYW3HsvrF4No0ZFHY2I1FbxJI7DgPfNbGE4FPfDiIbjSg0cfDAMPXc5EybA7f/+KOpwRKQWimc47qCkRSEpU7CsgGfaHQ+t3+CSC1rQ44B3OLbnoVGHJSK1SCyz4xqAu39e0Va2jKS//KX5bM3YCCcPg+9aM+qPjaMOSURqmViaqqab2UVmtt1ddGaWbWZHmtnDwPDkhCeJVjr1emaHuWTl3cj7r+6vUVYiEpcqR1WZWSPg18BZQBdgA9AIyAReAe529/eSHGeNaFTV9gqWFZC/NJ++HQZwydDDWLIE5s2DthXd0y8i9VZlo6piHo4bnqQB0ArYHC7mVCsocVTu44+DyRCPOgomTQI1OopIqUQMx8Xdt7r7itqUNGTn9t0Xxo+HF16Af/4z6mhEpDbQmuHCxRdDXl4wBfvSpVFHIyLpTolDyMiAhx4KmqlGjICSas19LCL1RbUSh5k1THQgEq1OneCOO+C114J/RUQqU90ax91mdnRCI5HIjRgBgwcHS81+pJvKRaQS1V069hygo5ndaWatEhyTRMQM7rsPmjWDYcO0YqCIVKy6TVXHEtzT0RW4z8xOSmhUEpk2bYI1yufMCZadFREpL541x980swHhbjvgQXc/3t1PBgbs5KVSy/ziF3D8aau5/nrnr4+qzUpEthdPjeM84EIzmwosdPfPyjw3NrFhSZQKlhUwbd998dbz+MPI3Xn27TlRhyQiaSTmxOHu89z9F8CfgTFm9oKZHRA+922yApTUy1+az9bMb+D0U6GoIb87p436O0Rkm+r0cSwCxgFfAfpTtA7aNhFi68U0OPl8vpi/B5ddFnVUIpIuYl6Pw8ymA92AzcBH4TYiOWFJlHJzcpk6bCr5S/PJOzuPx1rD7bdD375wyilRRyciUYtnzfGDgI/dfXNyQ0o8TXJYM1u2QL9+wTrlc+ZA165RRyQiqVDjSQ7d/d3amDSk5ho2hIkTITMTTj0VNuu3QKRe01xVEpNOneBf/4K5c4NJEUWk/lLikJidcEIwHcn998Mjj0QdjYhEJW0Sh5ndbGYLzOwDM3vGzHarpNwgM1toZovMbFSq46zvrr02mIL9t7+FDz+MOhoRiULaJA7gVWB/d+8FfAKMLl/AzDKBvwPHAT2AM82sR0qjrOeysuDxx6F5czjtNNi4MeqIRCTV0iZxuPsr7l4U7s4E9qig2KHAInf/zN0LgSeAIamKUQJt2wbJ49NPg8kQi4ujjkhEUiltEkc5vwZequB4B2BZmf0vw2M7MLORZjbbzGavWbMmCSHWb3l5wb0dzz6Lbg4UqWdSmjjMbIqZzatgG1KmzBigCHi0Ju/l7ve6ex9379O6deuahi4VOOTkAg47ZTa33hrMqCsi9UPMd44ngrsftbPnzWwE8HNgoFd8Z+JXQE6Z/T3CY5JiBcsKGPjIQLbsX0zG/Ge58MJBdO5sDBoUdWQikmxp01RlZoOAS4ET3f37SorNArqZWRczywbOACalKkb5Uf7SfAqLCynJKIRfnMHue67i9NM10kqkPkibxAHcBTQDXjWz981sAoCZtTezyQBh5/mFwMvAx8BEd58fVcD12baJEC2Tho238o9/f0WzZsG9HitWRB2diCRTzHNV1Waaqyo5CpYVBBMhds4jNyeX996DI46AffaB116DJk2ijlBEaqKyuaqUOCShXngBhgyBwYPhP/8J5rcSkdqpxpMcisTi5z+HO+6A556DSy+NOhoRSYaUjqqS+uHCC4ObA2+7DfbaCy64IOqIRCSRlDgkKW67DZYsgYsugo4dg5qIiNQNaqqSpMjMhMceg4MOCtbweOWVqCMSkURR4pCkadoUrntwFi1zVjH4xBKmT486IhFJBCUOSZqCZQWcPKk/q045gK3NP+b4E4p5/fWooxKRmlLikKTZdnd541XY8KNp2noDxx8PBQVRRyYiNaHEIUmz3d3lzTfwwNNLaNsWBg0C3VYjUntpVJUkTW5OLlOHTS1zd3kfDpwG/fvD0UfD9OnQu3fUUYpIvHTnuKTc0qXQrx98/z3k58P++0cdkYhURHeOS9ro3BmmTYOGDWHgQPj446gjEpF4KHFIJLp2DZKHGRx5JCxcGHVEIhIrJQ6JzN57w9SpwZrlhx8Ob78ddUQiEgslDonUfvvBW29B8+YwYEAwu66IpDclDoncmoYFDL3tDjp13cRJJ8EDD0QdkYjsjIbjSqRK1y4vLC6kwUnX0aflp5x77m4sXw5XXBH0gYhIelGNQyJVend5sRezNetrfn7VPQwbBldeCeefH/R/iEh6UY1DIlV6d3lhcSHZmdkM7NqPMQ9B+/Zw442wciU8/jjsskvUkYpIKSUOidSOd5fnAjB+PHToABdfDEcdBc8/Dy1bRhysiAC6c1zS3NNPw69+BV26wOTJwb8ikhq6c1xqpdJFoFauhIMPDpKHiEQrbRKHmd1sZgvM7AMze8bMdquk3FIz+9DM3jczVSPqgX79YM4c6NQJTjgh6DhXp7lIdNImcQCvAvu7ey/gE2D0TsoOcPfeFVWhpG5a1aCAk2++mRNOX824cXDccbB2bdRRidRPaZM43P0Vdy8Kd2cCe0QZj6SP0ns9rn1rNNN6dWbUTYuZMSNYz/ydd6KOTqT+SZvEUc6vgZcqec6BV8xsjpmNrOwEZjbSzGab2ew1a9YkJUhJjbL3ehQWF7Jr7kTefBMyM6FvX/jHP6AejPEQSRspTRxmNsXM5lWwDSlTZgxQBDxayWn6uvtBwHHA/5lZv4oKufu97t7H3fu0bt064dciqVN2JcHszGzyOudx8MFBv8dRR8EFF8CwYcH6HiKSfGk1HNfMRgDnAQPdvcqvATO7Gtjk7rfsrJyG49Z+BcsKdrjXA6CkBK6/Hq66KlgQ6oknoEePCAMVqUPSfjiumQ0CLgVOrCxpmFkTM2tW+hg4BpiXuiglKrk5uYw+YvR2SQMgIwPGjoWXXoIVK+DAA4M7zouKKjmRiNRY2iQO4C6gGfBqONR2AoCZtTez0tH7bYA3zGwu8A7worv/L5pwJZ0ceyzMnw+DB8Po0fCzn8FHH0UdlUjdlDZTjrh710qOLweODx9/BhyQyrik9li8pYCDLs7ngAG/4I6runPggXDNNfCnP0FW2vymi9R+6VTjEKm20iG7V+aPZfyG3jz86izVPkSSRIlD6oTyQ3Y/2DSFp54KOss/+0x9HyKJpMQhdUJFQ3bNYOjQoLZRtvbx7rtRRytSuylxSJ1QOj37uAHjmDps6najr3bfHZ56Cp58EpYuhT594Ne/DkZhiUj80uo+jmTRfRxS6ptv4Lrr4I47IDsbLr8c/vAHLRQlUpG0v49DJBU++raAlieO59Gp73H00TBmDOy7L0ycqGlLRGKlxCH1RunIq7HTxzJ8xuFcemcB06bBbrsFfSFHHAGqmIpUTYlD6o3yI6/yl+YzYEAw59V998Gnn8Ihh8Dw4bBkSdTRiqQvJQ6pNyoaeQXBLLvnnhskjssuC4bwdusGI0bAwoWRhiySltQ5LvVKZZMllvXVV3DLLXDPPfDDD0Ez1uWXQ8+eKQ5WJGKVdY4rcYiUU5pcejU9itefPIS//x02bYKTToIrrgjWPhepDzSqSiQGZTvQT5vcnyH/V8DnnwfTtufnB/eAHH88vPlm1JGKREeJQ6SMijrQW7aEq6+Gzz+H8eNh1qxg5cFDDoEHH9QCUlL/KHGIlFFZBzrArrvCqFHB3ed33QWbN8M550CHDvD736sjXeoPJQ6RMnY2dUmpD9YX8G2v8dw7uYAZM2DQILj7bthnHxg4EJ5+GrZujSB4kRRR57hIHEr7QAqLC8nOzN6WXFatggceCEZiffEFtGsXDPH91a+ge/eooxapHnWOiyRARX0gAG3aBEN2P/sMnn8eevcO5sTae++gQ/3WW+HLL6ONXSRRlDhE4rCzPhAIbib8+c9h7D0FXPrUXVx4xVIgWIWwY0fo3x8mTIC1a1Mfu0iiqKlKJE5V3URYUXPWTzbn8sQT8PjjsGBBsJTt0UfDGWcEw3tbtYrgQkSqoKYqkQTJzcll9BGjK73zvKLmrO7d4corg0Wl3nsPLrkE5s8P5sVq0yZYYOr66+H99zVLr6Q/JQ6RBNtZc5ZZ0P9x003w2IwCzp/wECMuWkZhYXBX+oEHQk4OjBwJzz0X3LEukm7UVCWSBNVpzurSIJeXXoIXX4RXXoGNG4PFpvr1gwEDgv6RQw4JjomkQmVNVVlRBCNS12jf6n0AAAy1SURBVOXm5FbalAUVN2flHpHL2WfD2WdDYSG88Qbc98Rypk/JYsqU3YFgpcLc3CCJ9O8PP/0pNGqUqqsSCaRN4jCzccAQoARYDYxw9+UVlBsOXBHuXufuD6cuSpHEKG3OKq1xlB+dlZ0Nu3Qr4LlOAykcUUijH9pxVZcprJy/N6+9FkyB4g4NGwbJo18/OPTQoEbStm0klyT1SNokDuBmdx8LYGYXA1cCvy1bwMxaAlcBfQAH5pjZJHf/OtXBitRE6R3qO2vOKlsrodEKfJ//8tffjAbg66+DGslrr8GLr27i+hsa4yVBl2WHDkECKd0OPhhatkzp5UkdlzaJw92/LbPbhCAxlHcs8Kq7rwcws1eBQcDjyY9QJLGqas7aWa2kRQsYPBha9S7g7hYDscFZNFjdh/PbPcCaRV2YNQueffbHc+21V5BAevYMtv33hy5dIEPDY6Qa0iZxAJjZ9cAw4BtgQAVFOgDLyux/GR6r6FwjgZEAHTt2TGygIikQT62kpMFmiveYQZsBT/DXa36slcyZE8zm+/KMdbw6I5OJE3fb9trGjaFHjx8Tyf77B/sdOgSjv0Qqk9JRVWY2BaioBXaMuz9XptxooJG7X1Xu9X8Kj18X7o8FNrv7LTt7X42qkrqqsrmzKivToKglf+szFVvTk3nz2LatXPlj+YaNiuneLZNu3dhha9tWSaU+SYtRVe5+VIxFHwUmE/RnlPUVkFdmfw8gv8aBidRScfeVZK1nTYsXGH3i9uvgTn5/Fif/7Qq2rupK0dd7s+suv2L+/JY8//z2M/02blLMrm3X0rVLA3rv05JOndhu2313JZb6IG2aqsysm7t/Gu4OARZUUOxl4AYzaxHuHwOMTkV8IumqJn0lpeZunEJxp6l4x1fAMjlhwHeMPmI0RUXBbL+ffgqvvLOEO16azMp1nVj1YSfee3tXvtu4/VdIw4awe/vNZLdcSbdOTenVtTXt27Nta9cu2HbZJdGfgqRS2iQO4EYz25tgOO7nhCOqzKwP8Ft3P9fd14fDdmeFr7m2tKNcRCoWS62ksuSSlQV77hls7zZ+AnwseDEZlsmYAeM4v+doPv+cbVvBvOVMfKuA4lUd+GxxBlOeL6Fo64498M2aF9Fotw10aJtN15xdad06qK2U/bd0a9ECGjRI9qck8dCd4yICVO9u9/Llxr8+nrHTx1LsxWRaJtfmjWNkj9EsXw4rVsDy5TDz4y94YMZLFG9sRcZ3bchpcBCbvm7MunWVx9a0aTCkuPy2OesrVhd/Qq+OnejdeU+aN2eHrVkzjR6rrrTo4xCR9FVVk1d1ai4DuuTRqlUw+2+vXkGZla8/Ck2DmotZJucNGLetWWzdOli9GvLnz+ePz9xE0abmZG5pzQldzqXR1vasXw/r1wcTRK5aW8j69a2huAMv7+S6zJwGjQrZtZnRonk2TZsGyaRp0x+3Taxg9dbP6N42h33adaRJk2DUWenWpAks+vYD5q6fyRF79uFnex7ELrsEd+03bPhjv05VybeuUOIQkZglI7mUbRZr0ybYXtgwiZL9HsO9GLdMDhjQgNFHbN+dOf71W7li2lhKCrPJKGzJxQeM5czu5/HNN2zb5n6+lH+88Rhbf2jM11t35cAOJ9KwuBWbNsGqVbB4Maz/Zgtrvm4MW39KfsnOvhJ7Ab24o9xRsyCBZGUXsbGkE2QNxbIK6dr6O1o0bULDhmy3bSxey7rC5eS0aENOyzZkZ7PD9tWmJXz27Sfs26Yr+7TZiwYN2GH7ZP1HfLjuXfp06E3vDvuTlcUO29zVc5j33TSO3rtvQhOZEoeIJFQyk0v5Mg2zsim0QrIbr+f0/r04NGf7MuNffxzPuhL3YrBMBgxYVUECui1oXispJqOkEWMOu57zDriE779n23b/zMd5aNZESrZmY8WNOXHPoeTtMYjNm+GHH2DzZnjjs/eZueR9vCgbincha7cMWjTtzpYt8N13QU3p603fs3TtN3hxc94tdnbN2kpJUQMKC4P5yX7UBejC8zv9pHsAPbh/p2UOJuPM67hxv6sqbFqsLiUOEUm5RCSXRCag7MxsCikku4FxXM9cOpS/rXiPzjz+/cvbznPZsJHklktSBcu2MvCRi7eVeWDY1B3KjH/9ju36gEaFzXQQzD1WVAQ3TL+Fa6ffSElRBhneiD8eOorfHHgBW7eybXtg1iPcO+uflJRAhjfk/+1/Nr/YZyhFRWzbnpn/PE/Pf5aSNnN+nEgzQYlDneMiUqfF0u+QqjLx3rCZ7DJVqaxzXIlDRCSF0imRVUWJQ4lDRCQuWnNcREQSQolDRETiosQhIiJxUeIQEZG4KHGIiEhclDhERCQu9WI4rpmtIZiqvTpaAWsTGE660/XWXfXpWkHXmwid3L11+YP1InHUhJnNrmgcc12l66276tO1gq43mdRUJSIicVHiEBGRuChxVO3eqANIMV1v3VWfrhV0vUmjPg4REYmLahwiIhIXJQ4REYmLEkfIzAaZ2UIzW2Rmoyp4vqGZPRk+/7aZdU59lIkTw/VeYmYfmdkHZjbVzDpFEWciVHWtZcr9wszczGr1EM5YrtfMTg9/vvPN7LFUx5hIMfwudzSz6Wb2Xvj7fHwUcSaCmT1oZqvNbF4lz5uZ3Rl+Fh+Y2UFJCcTd6/0GZAKLgT2BbGAu0KNcmQuACeHjM4Ano447ydc7AGgcPj6/tl5vLNcalmsGzABmAn2ijjvJP9tuwHtAi3B/96jjTvL13gucHz7uASyNOu4aXG8/4CBgXiXPHw+8BBhwGPB2MuJQjSNwKLDI3T9z90LgCWBIuTJDgIfDx08DA83MUhhjIlV5ve4+3d2/D3dnAnukOMZEieVnCzAOuAn4IZXBJUEs1/sb4O/u/jWAu69OcYyJFMv1OrBr+Lg5sDyF8SWUu88A1u+kyBDgEQ/MBHYzs3aJjkOJI9ABWFZm/8vwWIVl3L0I+Ab4SUqiS7xYrrescwj+iqmNqrzWsDqf4+4vpjKwJInlZ9sd6G5mb5rZTDMblLLoEi+W670a+JWZfQlMBi5KTWiRiPf/drVkJfqEUreY2a+APkD/qGNJBjPLAG4DRkQcSiplETRX5RHUJGeYWU933xBpVMlzJvCQu99qZrnAv8xsf3cviTqw2ko1jsBXQE6Z/T3CYxWWMbMsgirvupREl3ixXC9mdhQwBjjR3bekKLZEq+pamwH7A/lmtpSgXXhSLe4gj+Vn+yUwyd23uvsS4BOCRFIbxXK95wATAdy9AGhEMCFgXRTT/+2aUuIIzAK6mVkXM8sm6PyeVK7MJGB4+PhUYJqHvVG1UJXXa2YHAvcQJI3a3Aa+02t192/cvZW7d3b3zgT9OSe6++xowq2xWH6XnyWobWBmrQiarj5LZZAJFMv1fgEMBDCzfQkSx5qURpk6k4Bh4eiqw4Bv3H1Fot9ETVUEfRZmdiHwMsEojQfdfb6ZXQvMdvdJwAMEVdxFBJ1TZ0QXcc3EeL03A02Bp8IxAF+4+4mRBV1NMV5rnRHj9b4MHGNmHwHFwJ/dvVbWnmO83j8C95nZHwg6ykfU1j/6zOxxgqTfKuyzuQpoAODuEwj6cI4HFgHfA2cnJY5a+vmJiEhE1FQlIiJxUeIQEZG4KHGIiEhclDhERCQuShwiIhIXJQ6ROJnZbmZ2QZn99mb2dJLe6yQzu3Inz/c0s4eS8d4ildFwXJE4hVPqv+Du+6fgvd4iuCFx7U7KTAF+7e5fJDseEVCNQ6Q6bgT2MrP3zexmM+tcuj6CmY0ws2fN7FUzW2pmF4Zrm7wXTijYMiy3l5n9z8zmmNnrZrZP+Tcxs+7AltKkYWanmdk8M5trZjPKFH2eWnxDqtQ+Shwi8RsFLHb33u7+5wqe3x84BTgEuB743t0PBAqAYWGZe4GL3P1g4E/A3RWc53Dg3TL7VwLHuvsBQNm7+GcDR9TgekTioilHRBJvurtvBDaa2TcENQKAD4FeZtYU+Bk/TucC0LCC87Rj+zmV3gQeMrOJwH/LHF8NtE9g/CI7pcQhknhlZxIuKbNfQvB/LgPY4O69qzjPZoJZmAFw99+a2U+BE4A5ZnZwOMdUo7CsSEqoqUokfhsJpmOvFnf/FlhiZqfBtnWiD6ig6MdA19IdM9vL3d929ysJaiKl02d3Bypcg1okGZQ4ROIU/pX/ZthRfXM1T3MWcI6ZzQXmU/FytjOAA8ssUXyzmX0YdsS/RbC+NgTrw9eF1QulltBwXJE0ZmZ3AM+7+5RKnm8IvAb0DZc0Fkk61ThE0tsNQOOdPN8RGKWkIamkGoeIiMRFNQ4REYmLEoeIiMRFiUNEROKixCEiInFR4hARkbj8f40ylbKml8xcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#\n",
    "# No data this time. We'll just compare the analytical and computer\n",
    "# solutions to the problem. First let's do the Euler Method:\n",
    "#\n",
    "\n",
    "m = 0.0038  # kg (mass of cup plus 6 paper clips. See proj. 1)\n",
    "g = 9.8     # m/s^2, gravitational field strength at Earth's surface\n",
    "b = 0.00375 # effective drag coefficient, measured experimentally\n",
    "\n",
    "s = 0.0     # s is the 'state', in this case, it's just the velocity.\n",
    "dt = 0.03   # pick a reasonable time step\n",
    "t = 0.0     # start at t=0.0\n",
    "tf = 1.0    # stop in 10 sec\n",
    "\n",
    "def derivs_drag(s,t):\n",
    "    \"\"\"\n",
    "    The derivs_drag function computes the time derivative of the \"state\"\n",
    "    for this particular problem. In this case the \"state\" is just the velocity\n",
    "    of an object falling with drag.\n",
    "    \"\"\"\n",
    "    dsdt = b*s**2/m - g\n",
    "    return dsdt\n",
    "\n",
    "def EulerStep(s, t, derivs, dt):\n",
    "    \"\"\"\n",
    "    EulerStep is a generic implimentation of the Euler Algorithm.\n",
    "    Notice that it knows *nothing* about the state, nor the physics\n",
    "    of the problem (which is handled by the 'derivs' function\n",
    "    that is passed in by the programmer).\n",
    "    \"\"\"\n",
    "    return s + derivs(s,t)*dt\n",
    "\n",
    "v0=np.sqrt(m*g/b)\n",
    "tau=v0/(2*g)\n",
    "\n",
    "def VAnalytical(t):\n",
    "    return v0*(1-np.exp(t/tau))/(1+np.exp(t/tau))\n",
    "\n",
    "tlist = [t]\n",
    "slist = [s]\n",
    "\n",
    "while t <= tf:\n",
    "    s = EulerStep(s, t, derivs_drag, dt)\n",
    "    t += dt\n",
    "    tlist.append(t)\n",
    "    slist.append(s)\n",
    "    \n",
    "plt.xlabel(\"time (s)\")\n",
    "plt.ylabel(\"$v_y$ (m/s)\")\n",
    "plt.title(\"Falling Muffin Cup (Euler)\")\n",
    "v_euler = np.array(slist)\n",
    "v_analytical = VAnalytical(np.array(tlist))\n",
    "plt.plot(tlist,v_euler,'g.',tlist,v_analytical,'b-')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMS difference per time step= 0.0267378365429191\n"
     ]
    }
   ],
   "source": [
    "#\n",
    "# Let's get a numerical measure of the \"error\" per time step\n",
    "#\n",
    "\n",
    "rms_diff_euler = np.sqrt((v_euler-v_analytical)**2).sum()/len(tlist)\n",
    "print (\"RMS difference per time step=\", rms_diff_euler)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fce1b178610>,\n",
       " <matplotlib.lines.Line2D at 0x7fce1b178580>]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXhU5dnH8e+dhIRdRUAFAlEWEVHRBnRslSioaN21vm6vWLW4vNrFYiuKC4pbRau2toq17guoVbHgwhZxATRUVBRBFGS1BBFEWUPu949z0AgJzISZOTPJ73Nd52LOmSdn7pNofnmeszzm7oiIiMQrJ+oCREQkuyg4REQkIQoOERFJiIJDREQSouAQEZGEKDhERCQhCg7JCGZWamYXhK/PNbM3q7z3rZntEV11W2dmjczsJTNbaWbPhNuGmtkyM/vSzNqHx5Abda3xMLMLzeyuqOsAMLPLzOy2qOuQH1NwSNKZ2TwzWxP+sty0tKnt/ty9qbt/nswaAczsejNzM/vNZtt/E26/Ps5dnQrsAuzs7r8ws/bA74Fu7r6ru88Pj2FjLetsbmZ3mdn88Hv5Wbjesjb728Zn5QODgdvD9aLwe5G3WbuHzWxosj+/Gg8AZ5lZ6zR8lsRJwSGpclz4y3LTsjjqgmowGzhns239w+3x6gDMdveKcL098JW7L93e4sJf5OOBvYF+QHMgBnwF9Nre/VfjBOATd1+Ugn0nzN3XAi+z5c9IIqTgkLQws53M7N9mVm5mX4ev28X5tW5mncLXD5vZvWY22sxWmdlUM+tYpe2RZjYrHDb6m5m9vmkIrAbvAo3NbO/w6/cGGobbN+3zR0NnVWsysyHAtcD/hL2BC4GxQJtw/eHN/2oPh+VuNLO3wmN4bSu9h3MIgugkd//Y3Svdfam73+juYzb//lT5Hg0NX5eY2UIzuyocOptnZmdt5ftxNPD6Vt6vlpkdZGZvm9kKM3vfzEqqvDfPzPpWWb/ezB4PX2/63vQPe1TLzOzqzXZfCvw80ZokdRQcki45wEMEf523B9YAf63lvk4HhgA7AXOAmwDCX77PAoOAnYFZwMFx7O8xfviLtn+4Hhd3vw64GRgR9qzuJ/jluzhcP7eGLz0T+CXQGsgHBtbQri/wirt/G29N1dgVaAm0JTi+4Wa2Zw1t9yH4vsXNzNoCo4GhQAuCY3nOzFolsJufAXsCfYBrzWyvKu/NBPZLpCZJLQWHpMoL4V+fK8zsBXf/yt2fc/fV7r6K4Jd971ru+3l3fyccGnoC6BFuPwb4yN3/Fb53D/BlHPt7HDjDzBoQhNLjtawrEQ+5+2x3XwOM5Idj2NzOwJIkfN417r7O3V8n+CV/Wg3tdgRWVbN9WZWf5wqC4NvkbGCMu48Je0RjgTKCn0e8hrj7Gnd/H3ifHwfFKmCHBPYlKZa37SYitXKiu4/btGJmjYE/E4zT7xRubmZmubU4aVw1DFYDTcPXbYAFm95wdzezhdvambvPN7M5BD2HT919gZklWFLCajqGzX0F7Ladn/W1u39XZf0Lgu9VtW2BZtVsb1nlHA5m9nCV9zoAvzCz46psawBMTKDGrX0/mgErE9iXpJh6HJIuvycYijjQ3ZsDh4bbk/kbegnw/XkTC377x3UeBXg0rPHRat77DmhcZb+7bkeNiRoHHGVmTbbSZjVV6iMYmqpqp82+vj1Q08UKHwBdEqxxAfCYu+9YZWni7reG7//o+1dNfduyF0EvRDKEgkPSpRnBeY0VZtYCuC4FnzEa2MfMTgxPRP8f8f+SGgEcSTBstLn3gb3NrIeZNQSuT0axcXqM4Bfzc2bW1cxyzGzn8GT3pqGg6cCZZpZrZv2ofghwiJnlm9khwLHAMzV83pgavn5rHgeOM7OjwhoahiflN4X2dOB0M2tgZsUEly8nojfBlVWSIRQcki53AY2AZcAU4JVkf4C7LwN+AfyJYIinG8FY+7o4vnaNu48Lzzls/t5s4AaCv/4/Bd7cvE2quPs6ghPknxBcrfUN8A7Bye6pYbPfAMcBK4CzgBc2282XBENQiwnOCV3k7p/U8JEvAV0tgftu3H0BwWW8VwHlBEF3BT/8frkG6BjWMAR4Mt59h0F9DPBIvF8jqWeayEnqKjPLARYCZ7l7IuPtdUZ4Wezj7h7vkB1mNoDg5sXfpqyw+Gu5DCh09z9EXYv8QCfHpU4xs6MI/hJfQ/BXrxH0cCRO7j486ho2cfe/RF2DbElDVVLXxIDPCIbEjiO4umuL4ScRqT0NVYmISELU4xARkYTUi3McLVu29KKioqjLEBHJKtOmTVvm7ls8OqZeBEdRURFlZWVRlyEiklXM7IvqtmuoSkREEqLgEBGRhCg4REQkIQoOERFJiIJDREQSknHBYWb9wqk/55jZldW8X2BmI8L3p5pZUfqrFBGpvzIqOMwsF7iXYOrNbgSzsnXbrNn5BBPTdCKYGOi2VNUzecFkbnnjFiYvmJyqjxARyToZFRxAL2COu3/u7uuBpwke11zVCfzwiOVngT6WgunaJi+YTMmQa7n6+tX0ebSPwkNEJJRpwdGWKlN/EjwSu21NbcKpLFcSzMv8I2Y2wMzKzKysvLw84UJK55Wy4dPD8AlDWPflHpTOK014HyIidVGmBUfSuPtwdy929+JWrba4Y36bSopKKPjp36HBGuzNQZQUlSS/SBGRLJRpwbEIKKyy3i7cVm2bcHrQHQhme0uqWGGMCReP5Gcnz4APz6T1+liyP0JEJCtlWnC8C3Q2s93NLB84HRi1WZtRQP/w9anABE/Rs+FjhTFG/vlA8vKMW29NxSeIiGSfjAqO8JzFpcCrwExgpLt/ZGY3mNnxYbMHgZ3NbA5wObDFJbvJtNtucMEF8MgjMH9+Kj9JRCQ71IuJnIqLi317no47fz506gQXXgh/0USWIlJPmNk0dy/efHtG9TgyVfv20L8/PPAAfPll1NWIiERLwRGnK6+EDRtg2LCoKxERiZaCI04dO8KZZ8Lf/w7LlkVdjYhIdBQcCRg0CNasgbvuiroSEZHoKDgS0K0bnHJKcIJ8xYqoqxERiYaCI0GDB8M338Bf/xp1JSIi0VBwJGi//eC44+DPf4ZVq6KuRkQk/RQctXD11bB8Odx3X9SViIikn4KjFg48EI44Irg0d82aqKsREUkvBUctDR4MS5fCP/4RdSUiIuml4KilQw8Nlttug3Xroq5GRCR9FBzbYfBgWLQITh/8smYIFJF6Q8GxHZrsORlr9y4v/KMrhz90lMJDROoFBcd2eP2LUuzQobBid9a99wtNLysi9YKCYzuUFJWQv9dY2G0aTLqKn7UribokEZGUU3Bsh1hhjAn9x3P2pXPx5R2Z96amlxWRuk8TOSWBO+y/f3BPx8cfQ25uyj5KRCRtNJFTCpnBNdfA7NkwYkTU1YiIpJaCI0lOOgm6d4ehQ2HjxqirERFJHQVHkuTkBL2OmTPh2WejrkZEJHUUHEl0yimw115w441QWRl1NSIiqaHgSKLc3KDX8dFH8PzzUVcjIpIaCo4kO+006NIFbrhBvQ4RqZsyJjjMrIWZjTWzT8N/d6qh3UYzmx4uo9Jd57bk5gbPsPrgAxiVcdWJiGy/jAkO4EpgvLt3BsaH69VZ4+49wuX49JUXvzPOgE6dgl5HPbhNRkTqmUwKjhOAR8LXjwAnRljLdsnLC2YJfO89+Pe/o65GRCS5Mik4dnH3JeHrL4FdamjX0MzKzGyKmdUYLmY2IGxXVl5envRit+Wss2D33dXrEJG6J63BYWbjzGxGNcsJVdt58ByUmn7ddghvgT8TuMvMOlbXyN2Hu3uxuxe3atUquQcShwYNgl5HWRm88kraP15EJGXSGhzu3tfdu1ezvAj818x2Awj/XVrDPhaF/34OlAL7p6n8hP3v/0KHDjBkiHodIlJ3ZNJQ1Sigf/i6P/Di5g3MbCczKwhftwR+CnyctgoTlJ8PgwbB1Klw/p1Pa6InEakTMik4bgWOMLNPgb7hOmZWbGb/CNvsBZSZ2fvAROBWd8/Y4ADo2mcK1nwhD93dnsMf6aPwEJGslxd1AZu4+1dAn2q2lwEXhK/fBvZJc2nb5e0lE+GQhTD6XtbPLqF0XimxQs3bISLZK5N6HHVSSVEJBT0fhx3n4RNuoHeHkqhLEhHZLgqOFIsVxphw3iucPOAjfFExX72v3oaIZDfNAJgmGzZA167QvDlMmxY8hl1EJJNpBsCINWgA110H06fDCy9EXY2ISO0pONLorLNgzz2DANGTc0UkWyk40ig3F66/HmbMgJEjo65GRKR2FBxpdtppsPfeQYBUVERdjYhI4hQcaZaTEzyCZNYsePLJqKsREUmcgiMCJ50EPXoET87dsCHqakREEqPgiEBOThAan30Gjz4adTUiIolRcETk2GOhZ0+48UZYvz7qakRE4qfgiIhZ0Ov44gt48MGoqxERiZ+CI0JHHQUHHww33QRr10ZdjYhIfBQcETILhqoWLYLhw6OuRkQkPgqOiB1+OJSUwJCh6xky9nbN1yEiGU/BkQFO+78ZLC/PZ8iwcvo8qsmeRCSzKTgywIpdXoKOY/E3/sC67xpSOq806pJERGqk4MgAJUUlFBx5PaxpSc7kP1BSVBJ1SSIiNVJwZIBYYYyJg4axz+Efkzf1CoryNNmTiGQuBUeGiBXG+Nf93ajYkMsNN0RdjYhIzRQcGaRTJxgwAB54AGbPjroaEZHqKTgyzDXXQMOGMHhw1JWIiFQvY4LDzH5hZh+ZWaWZbTHHbZV2/cxslpnNMbMr01ljOuy6K1x+OTzzDLz7btTViIhsKWOCA5gBnAxMqqmBmeUC9wJHA92AM8ysW3rKS5+BA6FlS7jySnCPuhoRkR/LmOBw95nuPmsbzXoBc9z9c3dfDzwNnJD66tKrefNgyGrCBBg7NupqRER+LGOCI05tgQVV1heG27ZgZgPMrMzMysrLy9NSXDJdeCEUFcEf/wiVlVFXIyLyg7QGh5mNM7MZ1SxJ7zW4+3B3L3b34latWiV79ylXUABDh8L06TBiRNTViIj8IC+dH+bufbdzF4uAwirr7cJtddIZZ8Dtt8PVV8Mpp0B+ftQViYhk31DVu0BnM9vdzPKB04FREdeUMjk5cOutMHcu3H9/1NWIiAQyJjjM7CQzWwjEgNFm9mq4vY2ZjQFw9wrgUuBVYCYw0t0/iqrmdDjqKDjssGDejlWroq5GRCSDgsPdn3f3du5e4O67uPtR4fbF7n5MlXZj3L2Lu3d095uiqzg9zIJeR3k53HFH1NWIiGRQcEjNevWCU0+F24dt5KoX79Z8HSISKQVHljj10vdYvca59eY8TfYkIpFScGSJz3NewX7yAP7uhaxb0lGTPYlIZBQcWaKkqISCvrdAwbfw2p307lASdUkiUk8pOLJErDDGhItHcMwFZVR+egRff6jJnkQkGub14Cl6xcXFXlZWFnUZSbF+PeyzT/D6ww91U6CIpI6ZTXP3LZ5Wrh5HlsnPhzvvDCZ6uvfeqKsRkfpIwZGFjjkmuDFwyJDg/g4RkXRScGQhs6DX8e23cO21UVcjIvWNgiNLdesGl1wCw4fDBx9EXY2I1CcKjix2/fWw447wu99ppkARSR8FRxZr0SI4zzFhArz4YtTViEh9oeDIchddFAxbDRwI69ZFXY2I1AcKjiyXlwd//jN89hncfXfU1YhIfaDgqAOOPBKOPTaYava//426GhGp6xQcdcQdd8DatfCr3/6XW964RU/PFZGUUXDUEV26wCnnLualEa0Y/ORzevS6iKRMwsFhZk3MLDcVxcj26Xzi09B4GZWj72bdhg169LqIpMQ2g8PMcszsTDMbbWZLgU+AJWb2sZndbmadUl+mxOPofWI0OGowLPgpudMHUFJUEnVJIlIH5cXRZiIwDhgEzHD3SgAzawEcBtxmZs+7++OpK1PiESuMUfpnOPvzL1g28W72yI/nxysikphtPlbdzBq4+4btbROluvRY9Xh88gnsuy+cdho8rjgXkVqq9WPV4wmETA6N+qhrV7jySnjiCRg3LupqRKSuifvkuJkVm9nzZvYfM/vAzD40Mz1eL0NddRV06gQXXxxcpisikiyJXFX1BPAQcApwHHBs+G9SmNkvzOwjM6s0sy26RlXazQtDa7qZ1Z/xpwQ1bAh//zvMmQM33xx1NSJSlyQSHOXuPsrd57r7F5uWJNYyAzgZmBRH28PcvUd1Y2/yg7594eyz4dZbg/MeIiLJkEhwXGdm/zCzM8zs5E1Lsgpx95nuPitZ+5PAHXdA06bBwxD16HURSYZEguOXQA+gH8EQ1abhqnRz4DUzm2ZmA2pqZGYDzKzMzMrK6/H8qq1bw5/+BK+/Do88EnU1IlIXbPNy3O8bms1y9z2368PMxgG7VvPW1e7+YtimFBjo7tWevzCztu6+yMxaA2OBy9x9q8Nb9e1y3M1VVsKhhwbDVZ98Ai1bRl2RiGSDWl+OW8XbZtZte4pw977u3r2aJe5piNx9UfjvUuB5oNf21FQf5OTA/ffDypVwxRVRVyMi2S6R4DgImG5ms6K6HDd8TlazTa+BIwlOqss27L13EBoPPwylpVFXIyLZLJGhqg7VbU/WlVVmdhLwF6AVsAKY7u5HmVkb4B/ufoyZ7UHQy4DgcSlPuvtN29p3fR+q2mT1aujeHTbmrOGC+/9C3y6HECuMRV2WiGSomoaq4nnkiPk2GsXTJkoKjh/8+bGZXH7OXljJEBr2vY3x54xXeIhItbbnHMdEM7vMzNpvtsN8MzvczB4B+ierUEmttUUvYPs8iU+6inUL9taj10UkYfEERz9gI/CUmS0OH6c+F/gUOAO4y90fTmGNkkQlRSUUHHcFNCnHn3+Y2K6HRV2SiGSZeB5yuNbd/+buPwU6AH2A/d29g7v/yt3fS3mVkjSxwhgTLnqWc699A1+6N2MeOCjqkkQkyyQ0YUP4FNwlZlaQonokDWKFMWIDYxTMgWHD4Pjj4Wc/i7oqEckWtZ1z/G9mdkRSK5G0GzYMdt8d+veHb7+NuhoRyRa1Cg53Px9ob2b3mJnuQ85STZsG93XMnQsDB0ZdjYhki1oFh5kdBewOdAIeMLMTk1qVpM0hhwShcf/98PLLUVcjItkgkYmc3jKzTZfg7Ab8092PcfeTCOYelyx1ww3BneXnnw/Ll0ddjYhkukR6HBcCl5rZeGCWu39e5b1rkluWpFPDhvDYY1BeDpdeGnU1IpLp4g4Od5/h7qcAVwBXm9m/zWy/8L1vUlWgpMf++8O118JTT8Ezz0RdjYhkstqc45gD3AgsAqYltxyJ0qBB0KtXME/5kiVRVyMimSqRcxwTzWwhQVhcBSwHzk1RXRKBvLxgsqfvvoNTz/6amyfdwuQFk6MuS0QyTCI3AP4emOnua1JVjESva1e48I9zuXvI7kxuMo+GB/bRgxBF5EcSOcfxH4VG/dDqsKdhj3H4y3exbtGeehCiiPxIbe8clzrs8D1KaHjaBdDoa3zkSHru3CfqkkQkgyg4ZAuxwhgTLnmKX908npwVnRh+fS8yd7YVEUm3hB5yKPVHrDBG7NIYnVbDH/8Y3GF+2WVRVyUimUA9DtmqgQPh2GPh97+HqVOjrkZEMoGCQ7YqJye4RLdNGzjtNPjqq6grEpGoKThkm1q0CO4m//JLOOccqKyMuiIRiZKCQ+LSsyfceSeMGQO33RZ1NSISJQWHxO2SS+D002HwYCgtjboaEYlKxgSHmd1uZp+Y2Qdm9ryZ7VhDu35mNsvM5pjZlemusz4zg+HDoXPnIEC+/DLqikQkChkTHMBYoLu77wvMBgZt3sDMcoF7gaOBbsAZZtYtrVXWc82awbPPwjffwM9PXsnQibfqeVYi9UzGBIe7v+buFeHqFKBdNc16AXPc/XN3Xw88DZyQrhol0L07/H7oHP4zeQeuvSaPPo/2UXiI1CMZExybOQ+obiLTtsCCKusLw21bMLMBZlZmZmXl5eUpKLF+a9zzGazn3/C3BrJ26v/qeVYi9Uha7xw3s3HArtW8dbW7vxi2uRqoAJ7Yns9y9+HAcIDi4mI9MCPJSopKKDj2SNZ+vQc++l4a/8/sqEsSkTRJa3C4e9+tvW9m5wLHAn3cq3060iKgsMp6u3CbpFmsMMaEX77GK73e5onfHcK1l3Sjb49g7nIRqdus+t/P6Wdm/YA7gd7uXu3YkpnlEZw470MQGO8CZ7r7R1vbd3FxsZeVlSW5Ytlk/nw48EAoKAgeS7LLLlFXJCLJYGbT3L148+2ZdI7jr0AzYKyZTTez+wDMrI2ZjQEIT55fCrwKzARGbis0JPXat4eXXoKlS+H442H16qgrEpFUypgeRyqpx5EeL7wAJ58cLCNHBs+5EpHslQ09DslyJ54Iw4bBc8/BVVdFXY2IpIrm45Ck+t3vYM6c4HlWnTrBBRdEXZGIJJuCQ5LKDO65B+bOhYsugtWNP+a7whcpKSohVhiLujwRSQINVUnS5eXBiBFQ1Pk7fnNeWwaPeEp3l4vUIQoOSYnmzeGUGx+G/O+ofORV1v23g+4uF6kjFBySMiceeAAFvzwWPJfKR8bR2Y6KuiQRSQIFh6RMrDDGxIH3ctlfR9E8pzUDzz6AL76IuioR2V4KDkmpWGGMe869gInjG7ByJRx+OCxcGHVVIrI9FBySFgccAK++CuXl0KcPLFkSdUUiUlsKDkmbXr3glVdg0SLo2zd4RImIZB8Fh6TVwQfD6NHBfR59+8JXX0VdkYgkSsEhade7N4waBbNnwxFHwKsfvsMtb9yi+zxEsoSCQyLRty88/zzM+KiSo48xBr98u24SFMkSCg6JzNFHw+nXP4cv7kHlQ+NYt2In3SQokgUUHBKpi89uR/7Zp8LyLviDb7JHZb+oSxKRbVBwSKRihTFKh17JJX95lh1y2nHpqfszdWrUVYnI1ig4JHKxwhj3DjiXd6c0oHnz4CbB0aOjrkpEaqLgkIzRqRO8/TZ07QonnAD//GfUFYlIdRQcklF22QVKS4O7y88/H4YOhXowu7FIVtFETpJxmjWDl14KZg+85hpYvBjO/MNk3lhQqgmhRDKAgkMyUn4+PPIItGkTTEN7f2k5nHITBQ1vZPw54xUeIhHSUJVkLDO49VY49tdjqfzkWCr/MYl15W11r4dIxDImOMzsdjP7xMw+MLPnzWzHGtrNM7MPzWy6mZWlu05Jv6sGNg3u9VixO5X3v0ODOSdGXZJIvZYxwQGMBbq7+77AbGDQVtoe5u493L04PaVJlGKFMUpvuoKBDz9Fl04NuOK8vbjmGti4MerKROqnjAkOd3/N3SvC1SlAuyjrkcwSK4xx+/9cwvR3mnLeecHVVkcfDcuWRV2ZSP2TMcGxmfOAl2t4z4HXzGyamQ1IY02SARo1ggcfhAcegEmTggmi3nkn6qpE6pe0BoeZjTOzGdUsJ1RpczVQATxRw25+5u4HAEcD/2dmh9bwWQPMrMzMysrLy5N+LBKtCy6At96C3Fw45BC47z54e/5kPZ5dJA3MM+juKjM7F7gQ6OPuq+Nofz3wrbsP21q74uJiLyvTefS6aPlyOPtsePllyO3xOP7zSyhoVKFLdkWSwMymVXcuOWOGqsysH/AH4PiaQsPMmphZs02vgSOBGemrUjJNixbw739Dn/MmsfH9M6n823us++xAXbIrkkIZExzAX4FmwNjwUtv7AMysjZmNCdvsArxpZu8D7wCj3f2VaMqVTJGTAzde34D8844CoPKhibz30Ll8913EhYnUURk1VJUqGqqqHyYvmMxrM9/io6fP4pmHdqNjx+BBiYdWexZMRLYl44eqRLZXrDDGdUcOZOQ/d+P114OHI5aUwG9+g3ofIkmk4JA66dBD4YMP4NJL4Z57YL/94I03gl6JrrwS2T4KDqmzmjQJQqO0NOh99O7tHPI/0xj8ym30ebSPwkOklhQcUuf17h30Pg48aRobJ19K5T0zWVt2OhM+L426NJGspOCQeqFJE7jzrg0UXNgbdvwCf+GfPHbpZbzxRtSViWQfBYfUG7HCGBOvvpWbnpjI9X+ZzXcrmnLooXDaaTB3btTViWQPXY4r9dbq1TBsWDBR1MaNcPnlMGgQzFgxmdJ5mm1QpKbLcRUcUu8tXAhXXQWPPQY7t1rPNz/7NRv3/ScFDfL06BKp13Qfh0gN2rWDRx+FqVOhSeulbHj+Pir/MoO1ZWcwfs6kqMsTyTgKDpFQr17w1JgF5J95OhR8i7/wIPee8xvuuw/WrYu6OpHMoeAQqeLg9jFKb/0NNz39KsMenklRu4ZcfDHssQfcfXdwXkSkvtM5DpGtcIcJE+DGG+H116F1azj1/C9o1ftZjup2sM5/SJ2mcxwitWAGffoEd59PmgRFXVfwt1s6MOSkX3LomVN45o3/RF2iSNopOETidMghcOLQv5PzqxgUTaDi7cs47dAD6NcPXnwRKiqirlAkPfKiLkAkm5QUlVDQ4UbWtzudBt914OyKUl4eUciJJwZXZ114YTCt7dwNuhdE6i6d4xBJ0OQFPw6Figp46SX4+99h7FjIzauErv/Ci+8jf4/JTDh3nMJDspJuAFRwSBrMng3nXzOVN1/qDGtawA7zOeSYxdx9xUH06BGcMxHJFjo5LpIGXbrAn4ZV0vAPnbBT/pecXT/i7Wd6ccABsNdeMGRIEC6guUEke6nHIZICVYezujSO8dxz8NRTfD8z4Z7dv+WzdkOp7PYUBS3K9WgTyUgaqlJwSAZYtAhGjoQ7hi9h0Se7BRt3+w+HHbGGmy/+KT17Qm5utDWKbKKhKpEM0LYt/O538Mxr8yj47T5Y36vIyV/D648fTCwGu+4K55wDI0bAihUazpLMpMtxRSIQK4wx8fLh4XBWDns2MV59FUaPDpbHHoPcXMcLK/BO35Df6Y+M/cNtHLK7hrMkehqqEskwGzcGT+q9+r63KB3bBL7sAUB+o3UcdmgBvXsH0+EWF0N+/paXB4skS01DVRnT4zCzG4ETgEpgKXCuuy+upl1/YHC4OtTdH0lflSKpl5sLBx8MNxfm0GfPg1m3YidyFxzGz/OHMWvarlx1VdCucWPotv8KpjccS2X7UvLb3cmEX41SeEjKZTBTrxkAAAvMSURBVEyPw8yau/s34etfA93c/aLN2rQAyoBiwIFpwE/c/eut7Vs9DslW1fUmysuD52a9/jo8M+a/fPnZLkFj20jrouUc07sVPXtCz56w775QUKBeidROVl1VZWaDgPbufvFm288AStz9wnD9fqDU3Z/a2v4UHFJXTV4wmcPvO5X1X/yEnMUH0tMuYc6HO1FeHryfnw97dP2WTxs+SeWu75LfZhZjfns7h+95YLSFS1bI+KEqADO7CTgHWAkcVk2TtsCCKusLw23V7WsAMACgffv2yS1UJEPECmNMuOjZKr2JnXCH+fPh3XeD5dmxy9j4/unwzgDWAX3uC+YX6d49WPbZJ/h3eaMpvLV4onolsk1p7XGY2Thg12reutrdX6zSbhDQ0N2v2+zrB4bbh4br1wBr3H3Y1j5XPQ6pzyYvmMzhD/dl/bJ25C3rQf82t7NyfntmzIBZs4KT8QDkbICdPyVn5885vfcB9D6gDZ07Q+fO0KYN5ORoyKu+yYgeh7v3jbPpE8AY4LrNti8CSqqstwNKt7swkTosVhhjwrnjqvzC/6EHvm5d8AiUm559kZETZ+Lle1K5vDMjH2rNk/f/sI9GjaBNh++Ym/NffKcG5O00kptP3pGjiveiQwdo3vyHtgqXui9jznGYWWd3/zR8fRnQ291P3axNC4IT4geEm/5DcHJ8+db2rR6HyNZNXjCZPo/2Yf3G9eTn5jP27PEUWoxPP+X7ZcyUT/lkdgV8vQdsLPjR1++4I3ToAM1aL2fyqhFUNptH3g5LueOUgRzRY2/atIFmzX54yKPCJTtk/MlxM3sO2JPgctwvgIvcfZGZFYevLwjbnQeEFyRyk7s/tK19KzhEtm1bv8w3hcu6DRvIX9uOvxw8iuZr9+GLL/h+Kft4KV8uKoB1O2zx9U2awG67QdOdv+HD1a9Q2WQJuU2XM/CIczhoz460bg2tWgXT8zZrBlMWbjtcFECplfHBkUoKDpHkiDtcviugwer2DIs9yU4Ve7NkCSxeHCxTZ85n3sL18O0usL5ZtZ/TIL+SDQ2XQOOl5DRayaFd96Fzu51p0YLvl/LKWVw/9ddU5C+jQePVvND/Ufp260le3pY1KYBqR8Gh4BBJi3jDZf3G9TSobM6T/cbSvsFPWLo0uEdl6VJ46T+TeWPmLHx1C1i7M61zupCzthVffQUbNmz98xs3hh12CJbcRquYuWoKlfnfkFuwmpP260vn3XajaVNo2jTo2SxeO4sbpvyBiryVNCjYwEOn3stPO/WgceNgX40aBUNsyQqgbAopBYeCQyRjJBIu+bn53z923h1Wr4bly2Hix+/zq5F/pGJ1E3LX78wl+wyiRc7urFzJ98v0eXOZs6Qc1jWFDU1puLEVG9Y0+uFKsjg1bLSRtbYcGqzBGqylY6u27Ny8CQ0bBsHSqBF8W7mMCQtGU5m7mpwGFZyx38ns0aotBQV8vyz87jPufPcWKnK+o0ED56YjrmO/tnuRn8+PlhnL3mPaf6fw06JeHNThJzRowI+WvLz0DOUpOBQcIlklGX+9VxdAB7WLsX49rFoF334Lb86ZzvnP/JaKtQXkbWzO1bGbaNuwC6tX8/0yYdZU3vrsPbyiAKtoTNcd96ewcRfWroU1a4JlydfL+WrVaqgogI0F5FY2ZuOGFF64mlMRXEKdU0Gzho0oyM/7PlTy8qCCNSz8dh52zK8p6PxWreZ8yYjLcUVE4hUrjG3zF9222sQKY4w/Z/wW4bKpB9CyJRQV9aDjnrdUadNli/30W1BJn0cv/z6AHjxnPLHCH7eZvGDWFiF1YNsgpNatC5a3507j9BH92bABGngT7u33AF132pf16/l+eXr6v3jqg+fwilyssoDjO59Cnw792LCB75eJn73JxM/exCsNq8yne5tD6NG6mIoKvl/eX/I5C5d8TGX+16zfuJ7SeaVJGxpTj0NEJA7pOsdR0zBdKtpsi4aqFBwikiUy5US8gkPBISKSEE0dKyIiSaHgEBGRhCg4REQkIQoOERFJiIJDREQSouAQEZGE1IvLcc2snOBR7bXREliWxHIynY637qpPxwo63mTo4O6tNt9YL4Jje5hZWXXXMddVOt66qz4dK+h4U0lDVSIikhAFh4iIJETBsW3Doy4gzXS8dVd9OlbQ8aaMznGIiEhC1OMQEZGEKDhERCQhCo6QmfUzs1lmNsfMrqzm/QIzGxG+P9XMitJfZfLEcbyXm9nHZvaBmY03sw5R1JkM2zrWKu1OMTM3s6y+hDOe4zWz08Kf70dm9mS6a0ymOP5bbm9mE83svfC/52OiqDMZzOyfZrbUzGbU8L6Z2T3h9+IDMzsgJYW4e71fgFzgM2APIB94H+i2WZtLgPvC16cDI6KuO8XHexjQOHx9cbYebzzHGrZrBkwCpgDFUded4p9tZ+A9YKdwvXXUdaf4eIcDF4evuwHzoq57O473UOAAYEYN7x8DvAwYcBAwNRV1qMcR6AXMcffP3X098DRwwmZtTgAeCV8/C/QxM0tjjcm0zeN194nuvjpcnQK0S3ONyRLPzxbgRuA2YG06i0uBeI73V8C97v41gLsvTXONyRTP8TrQPHy9A7A4jfUllbtPApZvpckJwKMemALsaGa7JbsOBUegLbCgyvrCcFu1bdy9AlgJ7JyW6pIvnuOt6nyCv2Ky0TaPNezOF7r76HQWliLx/Gy7AF3M7C0zm2Jm/dJWXfLFc7zXA2eb2UJgDHBZekqLRKL/b9dKXrJ3KHWLmZ0NFAO9o64lFcwsB7gTODfiUtIpj2C4qoSgJznJzPZx9xWRVpU6ZwAPu/sdZhYDHjOz7u5eGXVh2Uo9jsAioLDKertwW7VtzCyPoMv7VVqqS754jhcz6wtcDRzv7uvSVFuybetYmwHdgVIzm0cwLjwqi0+Qx/OzXQiMcvcN7j4XmE0QJNkonuM9HxgJ4O6TgYYEDwSsi+L6f3t7KTgC7wKdzWx3M8snOPk9arM2o4D+4etTgQkeno3KQts8XjPbH7ifIDSyeQx8q8fq7ivdvaW7F7l7EcH5nOPdvSyacrdbPP8tv0DQ28DMWhIMXX2eziKTKJ7jnQ/0ATCzvQiCozytVabPKOCc8Oqqg4CV7r4k2R+ioSqCcxZmdinwKsFVGv9094/M7AagzN1HAQ8SdHHnEJycOj26irdPnMd7O9AUeCa8BmC+ux8fWdG1FOex1hlxHu+rwJFm9jGwEbjC3bOy9xzn8f4eeMDMfkdwovzcbP2jz8yeIgj9luE5m+uABgDufh/BOZxjgDnAauCXKakjS79/IiISEQ1ViYhIQhQcIiKSEAWHiIgkRMEhIiIJUXCIiEhCFBwiCTKzHc3skirrbczs2RR91olmdu1W3t/HzB5OxWeL1ESX44okKHyk/r/dvXsaPuttghsSl22lzTjgPHefn+p6REA9DpHauBXoaGbTzex2MyvaND+CmZ1rZi+Y2Vgzm2dml4Zzm7wXPlCwRdiuo5m9YmbTzOwNM+u6+YeYWRdg3abQMLNfmNkMM3vfzCZVafoSWXxDqmQfBYdI4q4EPnP3Hu5+RTXvdwdOBnoCNwGr3X1/YDJwTthmOHCZu/8EGAj8rZr9/BT4T5X1a4Gj3H0/oOpd/GXAIdtxPCIJ0SNHRJJvoruvAlaZ2UqCHgHAh8C+ZtYUOJgfHucCUFDNfnbjx89Uegt42MxGAv+qsn0p0CaJ9YtslYJDJPmqPkm4ssp6JcH/cznACnfvsY39rCF4CjMA7n6RmR0I/ByYZmY/CZ8x1TBsK5IWGqoSSdwqgsex14q7fwPMNbNfwPfzRO9XTdOZQKdNK2bW0d2nuvu1BD2RTY/P7gJUOwe1SCooOEQSFP6V/1Z4ovr2Wu7mLOB8M3sf+Ijqp7OdBOxfZYri283sw/BE/NsE82tDMD98XZi9ULKELscVyWBmdjfwkruPq+H9AuB14GfhlMYiKaceh0hmuxlovJX32wNXKjQkndTjEBGRhKjHISIiCVFwiIhIQhQcIiKSEAWHiIgkRMEhIiIJ+X/hVl0iVa0B4AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#\n",
    "# No data this time. We'll just compare the analytical and computer\n",
    "# solutions to the problem.\n",
    "#\n",
    "\n",
    "m = 0.0038  # kg (mass of cup plus 6 paper clips. See proj. 1)\n",
    "g = 9.8     # m/s^2, gravitational field strength at Earth's surface\n",
    "b = 0.00375 # effective drag coefficient, measured experimentally\n",
    "\n",
    "s = 0.0     # s is the 'state', in this case, it's just the velocity.\n",
    "dt = 0.03   # pick a reasonable time step\n",
    "t = 0.0     # start at t=0.0\n",
    "tf = 1.0    # stop in 10 sec\n",
    "\n",
    "def derivs_drag(s,t):\n",
    "    \"\"\"\n",
    "    The derivs_drag function computes the time derivative of the \"state\"\n",
    "    for this particular problem. In this case the \"state\" is just the velocity\n",
    "    of an object falling with drag.\n",
    "    \"\"\"\n",
    "    dsdt = b*s**2/m - g\n",
    "    return dsdt\n",
    "\n",
    "def HeunStep(s, t, derivs, dt):\n",
    "    \"\"\"\n",
    "    HeunStep is a generic implimentation of the Heun Method.\n",
    "    Notice that it knows *nothing* about the state, nor the physics\n",
    "    of the problem (which is handled by the 'derivs' function\n",
    "    that is passed in by the programmer).\n",
    "    \"\"\"\n",
    "    f1=derivs(s,t)\n",
    "    f2=derivs(s+f1*dt,t+dt)\n",
    "    return s + 0.5*(f1+f2)*dt\n",
    "\n",
    "v0=np.sqrt(m*g/b)\n",
    "tau=v0/(2*g)\n",
    "\n",
    "def VAnalytical(t):\n",
    "    return v0*(1-np.exp(t/tau))/(1+np.exp(t/tau))\n",
    "\n",
    "tlist = [t]\n",
    "slist = [s]\n",
    "\n",
    "while t <= tf:\n",
    "    s = HeunStep(s, t, derivs_drag, dt)\n",
    "    t += dt\n",
    "    tlist.append(t)\n",
    "    slist.append(s)\n",
    "    \n",
    "plt.xlabel(\"time (s)\")\n",
    "plt.ylabel(\"$v_y$(m)\")\n",
    "plt.title(\"Falling Muffin Cup (Heun)\")\n",
    "v_heun = np.array(slist)\n",
    "v_analytical = VAnalytical(np.array(tlist))\n",
    "plt.plot(tlist,v_heun,'g.',tlist,v_analytical,'b-')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMS difference per time step= 0.001909920488169262\n"
     ]
    }
   ],
   "source": [
    "#\n",
    "# Let's get a numerical measure of the \"error\" per time step\n",
    "#\n",
    "\n",
    "rms_diff_heun = np.sqrt((v_heun-v_analytical)**2).sum()/len(tlist)\n",
    "print (\"RMS difference per time step=\", rms_diff_heun)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wow! That 0.2% rms error for Euler compared to 2.6% for the Euler Method, over 10x better.\n",
    "\n",
    "Heun’s method is considered a 2nd order method. That means that\n",
    "the size of the error goes like $dt^2$. Euler’s Method, in\n",
    "contrast, is 1st order. Let’s graph the error using Heun’s Method and\n",
    "Euler’S Method as we change the time step size for an example problem\n",
    "(like the simple harmonic oscillator, which we’ll study next!).\n",
    " \n",
    "Note that halving the time step using  Euler’s Method reduces the\n",
    "error in proportion (i.e., by half). However, reducing the time step\n",
    "using Heun’s method reduces the error by a factor of 4, so you get 4\n",
    "times the accuracy by doubling the number of steps. Of course, each\n",
    "step now requires more work as well, so it’s not all for free.\n",
    "\n",
    "Let's quantify that even more. The following block of code computes the velocity as a function of time for a number of different time steps using both the Euler and Heun methods comparing the output with the analytical solution along the way. Notice that the log of the rms error per time step $\\log(\\sigma)$ is linear in the log of the time step with a slope of 1 (Euler) or 2 (Heun). We'll discuss what that means exactly in class, but you should think about it yourself as well. You don't need to understand this code in detail, but you should understand what it's doing at a high level (not the nuts and bolts, but the basic idea)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(-3.5, -1.5, 'log($\\\\sigma$) =  0.996  log(dt) +  0.028')"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEcCAYAAADHiMP9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdd3yUVfb48c9NAiR0pEgJJCABIQmEEMCCIIogyYK0XWBRwIZ1RVewsQoWVhRWQdB1UVdUYll/uCtKQL8gKKiAlCBIL6EqEmogPTm/P57JkDKTTJKZzCQ579drXmbmKXPmEnPmee695xoRQSmllPLzdgBKKaV8gyYEpZRSgCYEpZRSNpoQlFJKAZoQlFJK2WhCUEopBWhCUEopZaMJQSmlFKAJQfkIY0ySMaa/m871ojHm4VLsv8EYE17M9o7GmERjTIox5iF3xOhOxphfjDHXV/cYVPlpQvAhtj+KacaYC8aY34wxC40xdQttzzTGNCl03BZjjBhjQm3PextjfjDGnDPGnDbGfG+M6VGxn8Y7jDFNgXHAv0px2GzguWK2PwasEpF6IvJaeeLzBBEJF5HVFfV+jpJ3RcdQmDu/UFRnmhB8z2ARqQtEAd2AJwttPwiMyXtijIkEaud7Xh/4EpgHXAa0Ap4FMsobmDEmwJXXSnsON5sAJIhIWimOWQL0M8Y0d7I9BPilLMF48vNWQFuqakYTgo8Skd+Ar7ASQ34fYH0DzjMeeD/f8w624z8SkRwRSRORr0XkZ2fvZYxpaYxZbIw5aYw5mP+2iO2b1+PGmJ+Bi8aYACevdTLGrDbGnLXdPhhS3DmK++zFncu2Pdp2VZRijPnUGPOJMeYF2+ZBwLeF9q9hjJlhiyPLdjUltngQkXRgEzDQQSzfAP2A+bYrtw4uxOf089pe/3+F9p9rjHnN9vMTxpj9ts+2wxgzrKRz5/92XMK/gxhj2ud7vjBfu+XFdsz23ruNMTc6aI8PgDbAF7b2eCxfXP3z/TzFGPOzMeaiMeYdY8zlxphltnOvMMY0yndOp79/Dt6/SIzFxFTS7/WTtjY+Y4x51xgTWJq2qJJERB8+8gCSgP62n4OBbcDcwtuB3UAnwB84ivUNVoBQoD5wCngP649joxLe0w/rj+EzQE2gHXAAGJjvPROB1kCQo9eAGsA+4CnbOW4AUoCOzs7h7LO7cK6awCFgkm3f4UAm8IJt+0mgR6FzvwSss71/HWAF8BnQLt8+rwGvOIltNXCX7edi4yvp89r+rVKBerbn/sCvwFW2538EWtr+XUYBF4EWxZ27FG0nQPt851qYr906AkeAlrbnocAVJf2eOvndTbK19+VYV6i/A5uxrngDgW+Aaa78/hV6D6cxFo6ppPPa9t9ua8fLgO/L0hZV7aFXCL7nf8aYFKxfyN+BaQ72ybtKuAnYCRzL2yAi54HeWP/zvwWcNMYsMcZc7uT9egBNReQ5EckUkQO240bn2+c1ETkiBW/D5H/tKqAuMNN2jm+wbluNKeEcjpR0rquAANv5skTkM2BDvuMbYv0RBMAYUw94CLjN9v4XgcXAZbbPmifFdmxJXPmsTj+viBzC+uOY983/BiBVRNbZtn8qIsdFJFdEPgH2Aj1dOXcpYnMkB6gFdDbG1BCRJBHZ78JxzswTkRMicgxYA6wXkS1iXY39Fys5gGu/f2WJ0ZXzzre142lgBpfayd1tUWloQvA9Q0WkHnA9cCXQxME+HwB/xrpf/n7hjSKyU0QmiEgwEIH1jXOOk/cLAVrabjGcNcacxfqGmT+BHHFwXP7XWgJHRCQ332uHsL4dFncOR0o6V0vgmNi+ujk49xmgXr7nfYADIrI332uNgN8KvW894Kwb4nMUU2EfcumPz59tzwEwxowz1oimvH+LCIr+Djg7t6uxFSEi+4CHgenA78aYj40xLUs6rhgn8v2c5uB53mAJV37/yhJjaX+vD2G1nyfaotLQhOCjRORbrEv62Q62HcLqXI7FuvVR3Hl22c4T4WSXI8BBEWmY71FPRGLzn8bRqfP9fBxobYzJ//vUhnxXLk7O4UhJ5/oVaGWMMfm2t87388/Y+lFsmmIlCQBsxw3D+uacXydgqxviy1Pc5/0UuN4YE2yL5UNbbCFY32IfBBqLSEOs2xqm0PHOzl1SbKnkG4AAFOhEF5EPRaQ3l25BvuTkfdy5iIorv3+uxFg4JlfOm//3pg1W+5X0PlWaJgTfNge4yRjT1cG2O4EbbLdA7IwxVxpjHrX9scEY0xrr2+g6J++xAUixdaIFGWP8jTERpnTDVNdj/bF5zFgduNcDg4GPS3EOV8/1I9Yl/YO2DtVbKHhLJQHom+/5diDaGBNljAkCXsT6H/yTvB1snYndgf9zQ3wlEpGTWP0S72L90dpp21THFttJW1y34zyRlyW2RODPtn/jm8nXTsaaa3GDMaYWkI71LT4Xx05g3ZN3B5d//0qIsXBMrpz3AWNMsDHmMmAqtt+JUrZFlaIJwYfZ/nC8j9UxVnjbfhHZ6OCwFKAXsN4YcxErEWwHHnXyHjnAH7BGMx0EkoG3gQaliDMT6w/PINvxbwDjbFcnpVLSuWzbh2MlxLPArVjf9vOG1b4PxNr++GNroxlYieIA1rfiWBHJyve2g4HVInKcErjxs36I1RFsv10kIjuAf2AlvRNAJFZnp0tciG2SbftZYCzwv3yH1wJm2o77DWhG0SHPeV4E/ma7FTPZ1ficxFya37/iYiwQk4vn/RD4Guv3Yj+QN+KqNG1RpZiCt2KVqnyMMeuBN0XkXdvzvwO/i4izfhNHx98pIts9GKbyIcaYJKyRYyu8HYsv0YktqtIxxvTFGnqbjPVNtwuwPG+7iDxVmvOJSC+3BqhUJaUJQVVGHYH/YN1zPwCMFJFfvRuSUpWf3jJSSikFaKeyUkopm0p9y6hJkyYSGhrq7TCUUqpS2bRpU7KINC38eqVOCKGhoWzc6GjkpVJKKWeMMYccva63jJRSSgGaEJRSStloQlBKKQVoQlBKKWWjCUEppRSgCUEppZSNTyQEY8wfjbX+a64xJsbb8SillK+K3xZP6JxQ/J71I3ROKPHb4t12bl+Zh7Adq6Txv7wdiFJK+ar3tr7HPV/cQ0aOVe390LlDTPxiIgBjI8eW+/w+kRDyFggpuAiWUkqpExdOsHzfcpbuXcr/2/H/kEKLw6VmpTJ15dSqkxBKwxgzEZgI0KZNGy9Ho5RS7pUruWw8vpGEvQks3buUjcetagzN6zYvkgzyHD532C3vXWEJwRizgkJruNpMFZHPXT2PiCwAFgDExMRoqValVKV3Ju0MX+//moR9CSzbu4yTqScxGK4Kvorn+z1PXFgcXZt3pd3cdhw6V7TqRJsG7vlyXGEJQUT6V9R7KaWUr4nfFs/UlVM5fO4wreu35t4e92IwLN27lB+P/EiO5HBZ0GXc3P5mYtvHMrD9QJrUblLgHDNunMHELyaSmpVqf612jdrMuHGGW2KsdLeMlFKqsonfFs/dS+4mLTsNgMPnD/PUSmthv27Nu/FE7yeIC4ujZ6ue+Pv5Oz1PXj9BXmJp06ANM26c4Zb+A/CRBXKMMcOAeUBTrAXAE0VkYEnHxcTEiFY7VUr5qj2n9pCwN4EnVz5JenZ6ke2t6rXi6F+PVnhcxphNIlJkiL9PXCGIyH+B/3o7DqWUKo/07HS+TfqWhL0JJOxLYN/pfcXufzzleAVF5hqfSAhKKVVZHT532EoAexNYeXAlqVmpBAYEckPbG3i418PEhsXS771+Hu0MdhdNCEopVYz8ncFtGrTh2X7PEtog1H4VsP337QCENgzl9qjbiQuL4/rQ6wmqEWQ/h6c7g93FJ0pXqPKrW7dumY9NS0ujb9++5OTkON0nMzOTPn36kJ2dXeb3cWb58uV07NiR9u3bM3PmTIf7zJ07l4iICMLDw5kzZ45L286ePcvIkSO58sor6dSpEz/++GOJsXiqHadPn87s2bM5e/Ysb7zxhv11b7ers/2OHDlCv3796Ny5M+Hh4cydO9e+/6uvvkp4eDgRERGMGTOG9PSi98arivht8Uz8YiKHzh1CEA6dO8SE/03g+veu55V1r9CsTjNm3zSbnQ/s5MBDB5gfO59BYYMKJAOwOoMXDF5ASIMQDIaQBiEsGLzAbZ3BbiMilfbRvXt3UZY6deqU+dj58+fLnDlzStxv+vTpsmjRojK/jyPZ2dnSrl072b9/v2RkZEiXLl3kl19+KbDPtm3bJDw8XC5evChZWVly4403yt69e0vcNm7cOHnrrbdERCQjI0POnDlTYjyeasdp06bJrFmz5ODBgxIeHl5gmyvtumrVKhk/frzLsbjSrsXtd/z4cdm0aZOIiJw/f17CwsLkl19+kaNHj0poaKikpqaKiMgf//hHeffdd12Oq7LIyc2RdUfWSYMXGwjTKfJo+nJTOZd+ztthlhmwURz8TdUrhCrmlVdeISIigoiIiALflp9//nk6duxI7969GTNmDLNnz7Zvi4+P55ZbbrE/37p1K3369KFz5874+flhjOGZZ55h6NChxMe7r5AWwIYNG2jfvj3t2rWjZs2ajB49ms8/LzhPcefOnfTq1YvatWsTEBBA3759+eyzz4rddu7cOb777jvuvPNOAGrWrEnDhg1LFVtp27JwO86YMYMOHTrQu3dvdu/eDcATTzzB/v37iYqKYsqUKQBea9fi9mvRogXR0dEA1KtXj06dOnHs2DEAsrOzSUtLIzs7m9TUVFq2bOnW2L3ldNppPt7+MeP+O47ms5tz1TtXcS7jnMN9k1OTqV+rfgVH6Hnah1CFbNq0iXfffZf169cjIvTq1Yu+ffuSnZ3N4sWL2bp1K1lZWURHR9O9e3fAumVx4MABQkNDAUhPT2fUqFG8//779OzZk6effpr09HSeffZZcnNz+emnnxy+93XXXUdKSkqR12fPnk3//s7nJB47dozWrVvbnwcHB7N+/foC+0RERDB16lROnTpFUFAQCQkJxMTEFLvt4MGDNG3alNtvv52tW7fSvXt35s6dS506dTzSloXbcdOmTXz88cckJiaSnZ1t32/mzJls376dxMTEAp/PWbuWlSvt6up+SUlJbNmyhV69elG/fn0mT55MmzZtCAoKYsCAAQwYMMCtsVcUEeHnEz/b+wJ+OPIDuZJL46DG1uSwsFge/7/HOZpSdFior3UGu4smhCpk7dq1DBs2zP5Hb/jw4axZs4bc3FxuueUWAgMDCQwMZPDgwfZjkpOTC3xzXrFiBdHR0fTs2ROALl26sHz5cowx+Pv7U7NmTVJSUqhXr16B916zZo3HPlenTp14/PHHGTBgAHXq1CEqKgp/f/9it2VnZ7N582bmzZtHr169mDRpEjNnzuT555936T1L25aF23HNmjUMGzaM2rVrAzBkyBCn71Vcu/bq1YuMjAwuXLjA6dOniYqKAuCll15i4MASp+qU24ULFxgxYgRz5syhfv36nDlzhs8//5yDBw/SsGFD/vjHP7Jo0SJuvfVWj8dSWoU7g2fcOIMhHYaw8uBK+6igYynWVU90i2ie6v0UsWGxBSaHCVIpOoPdRRNCNRcUFFSgU3D79u1ERkban2/evNl+6wAgIyODwMDAIucp6xVCq1atOHLkiP350aNHadWqVZH97rzzTvvtn6eeeorg4OBit+U9evXqBcDIkSOL7Vgtr8LtWFrO2jXv2/rq1atZuHAhCxcudOl8rrZrcftlZWUxYsQIxo4dy/DhwwHrC0Pbtm1p2rQpYCXKH374wecSQl5ncN4f8kPnDjHuv+MwGHIkh3o16zHgigHEhsUyqP0gWtRr4fA8np4ZXBbx8TB1Khw+DG3awIwZMNZd4TjqWKgsD+1UvqROnTqyadMmiYyMlIsXL8qFCxckPDxcNm/eLBs2bJBu3bpJWlqapKSkSFhYmMyaNct+bHBwsKSlpYmIyIIFC2T06NEiIrJ7927p0KGDJCcni4hIcnKydOzY0a1xZ2VlSdu2beXAgQP2Ts3t27cX2e/EiRMiInLo0CHp2LFjgQ5iZ9t69+4tu3btEhGrU3fy5MklxpPXqVyWtszfjnnHp6amyvnz56V9+/Yya9YsSU5OljZt2hR4T1fatbSdyq62q7P9cnNz5bbbbpNJkyYV2H/dunXSuXNnuXjxouTm5sq4cePktddeczmuipCWlSbNXm7msDO4/t/ryzcHvpGM7Axvh1kmixaJ1K4tApcetWtbr5cGTjqV9QqhComOjmbChAn22z133XUX3bp1A6xbFl26dOHyyy8nMjKSBg0a2I8bMGAAa9eupX///owZM4YlS5YQERFBkyZN+Oijj2jcuDEAq1atIi4uzq0xBwQEMH/+fAYOHEhOTg533HEH4eHhAMTGxvL222/TsmVLRowYwalTp6hRowavv/56gdszzrbNmzePsWPHkpmZSbt27Xj33XddjqssbZm/HaOjoxk1ahRdu3alWbNm9OjRA4DGjRtz7bXXEhERwaBBg5g1a5ZX29XZfmvXruWDDz4gMjLSfpvq73//O7GxsYwcOZLo6GgCAgLo1q0bEydOdGvsZXHo7CF7X8DKAyvtNYMKS8lMoV/bfhUcnftMnQqpqQVfS021XnfLVYKjLFFZHnqF4LqUlBQREbl48aJ0797dPqRQxPo2e+utt5Z4jmHDhsnu3bs9FmNl4awtXW3HwrRdSy8zO1NWHVwlU76eIp1f72y/Amg7p608uPRBp1cIIa+GeDv0cjGm4NVB3sOY0p0HvUKo3iZOnMiOHTtIT09n/PjxBfoFoqOj6devHzk5OfbO2sIyMzMZOnQoHTp0qKiQfZaztnSlHQvTdi3KUWfw2Mix/HbhN5btXUbCvgS+3v815zPOU8OvBn1C+nBntzuJC4ujQ+MOGGO4qvVVVa4zWASaNIGTJ4tuc9daYT5R7bSstNqpUlVL4c5ggAC/AFrXb83BswcBaFmvJbHtY4nrEMeNbW+kXq16Ts/lS53B5bFnDzz8MCxbBsZYySFP7dqwYEHpbhk5q3aqCUEp5TNav9La4bj/Wv61eKbvM8SFxdHl8i7VZv31CxfghRfglVcgKAimT4fLLoNp08o3ysiny18rpaonEWHria329YMdJQOAzJxMnrruqQqOzntE4KOPYMoUOH4cJkyAF1+E5rZFiMeP98z7akJQSlWolIwUVhxYYR8VlLcmQPcW3WlQq4HDchFVdWawI1u3wl/+AmvWQEwMLF4MV11VMe+tCUEp5RbO7tmLCLtP7bbPDv7u0Hdk5WZRv1Z9BlwxgLiwOG5ufzPN6zZ32IdQ2TuDXXX6NDz9NLz5pnVb6K234I47wK8CK85pH4JSqtwc/SGv5V+LPiF92H9mPwfOHAAgvGk4sWGxxIXFcU3ra6jhX8PhuapKZ7ArcnLg7betuQRnzsADD8Czz0KjRp57T+1UVkp5TOicUIcrghkMcR3iiAuLY1D7QYQ0DPFCdL7BUcmJtm2t20ObN0PfvjBvHuSrHOMx2qmslHKrzJxMvj/8PQl7ExwmA7CKw30x5osKjsz3xMfDxImXZhkfOmR1DOfkQHAwfPwx/OlP1pBSb9KEoJRy2a8pv7Js3zIS9lqTw1IyU6jhV4PAgEDSs4sW9wtpUH2vCPJzVHIiJwfq14ddu8DFquwe5xMJwRgzCxgMZAL7gdtF5Kx3o1Kqeijunn1Obg4bjm2wDwvd8tsWAFrVa8XoiNHEhcVxQ9sbWLJnSbXtDHbF4cOOX09J8Z1kAD6SEID/A54UkWxjzEvAk8DjXo5JqSrPUZnou5fczQ9HfuBs+lmW71vO6bTT+Bk/rml9DX+/4e/EdYgjsllkgclhvlgm2lckJUFgIKQ5qLfnrpIT7uJzncrGmGHASBEp8TdJO5WVKh9nncEATWo3YVD7QcSFxTHgigE0CvLgsJcqKC0NXnrJeohAbi5kZV3aXpaSE+5SmTqV7wA+cbbRGDMRmAjQxtfSq1KVxPmM86w4sMJpMgA4MfkEfkaXXS8tEfjsM3j0UavzeMwYePll+PZbDy5s4yYVlhCMMSuA5g42TRWRz237TAWyAacrjovIAmABWFcIHghVqSpHRNiVvMveF7Dm8Bqyc7MxGISi/xuFNAjRZFAGO3bAQw/BypXQpYuVBPr0sbaNHet7CaCwCksIIuJ8HUXAGDMB+ANwo/jafSylKqHUrFRWJ61m6Z6lJOxLIOlsEgARzSJ49OpHiQ2L5eDZg9y/9H7tDC6nc+esyWTz5kHdujB/PtxzDwT44j2YYvhEuMaYm4HHgL4iklrS/kopx6ODrgm+xn4VsCppFenZ6dSuUZsb297IE9c+waCwQQXqAvUJ6UOAX4B2BpdRbi68/z488QT8/jvcdZd1K8i25HSl4xOdysaYfUAt4JTtpXUicm9Jx2mnsqquHJWKyH/7p/1l7YkLiyM2LJY+IX0IDAj0VqhV1saN1izjdeus4nPz5lnF6CoDn+5UFpH23o5BqcrieMpxHlr2UIFkANas4EaBjVh/13rCGod5Kbqqp3DJiccft0pNvPMONGsGCxfCbbdVbBE6T/GJhKCUci4nN4f1x9bbbwUl/pbodN+z6Wc1GbiRo5IT998P/v7wyCPwzDPQoIF3Y3QnTQhK+aDk1GS+2vcVS/cu5av9X3E67TT+xp9rWl/DzBtnMnf9XH698GuR46rTugEVwVHJCYDLL4d//KPi4/E0TQhKVTBHncFjIsaQ+Fui/Spg/dH1CEKzOs0Y3GEwsWGx3NTuJvvksOAGwVoqogIccjJN49eiubhK0ISgVAVyVCpi/H/H88DSB+wrhfVo2YNpfacRGxZL95bdHc4H0FIRnpWRYV0BFF7QPk9VnROrCUGpCiIiTPl6SpHO4BzJISs3i/eGvsfAKwZyed3LXTrf2MixmgA84Msv4eGHYf9+a9TQ9u2Qnq+Qa+3a1tDSqqgK9Isr5btSs1JZumcp9y+9n7Zz2zq87w+QlpXGuK7jXE4Gyv327oW4OBg8GGrUgK+/hp9+slYzCwmxrhZCQrxXf6gi6BWCUm524MyBS5PDDq4iIyeDOjXq0L9df1IyUziddrrIMdoZ7D0XLsDf/27dIqpZE2bNsspP1Kxpba8MJSfcRROCUi5ytm5AZk4maw6tYenepSTsTWD3qd0AhF0Wxn0x99knh9UKqFWtF5H3NSLwyScweTIcO2bNJXjpJWjRwtuReY9PzFQuK52prCqKoz/kNf1rEtkskt2ndnMh8wK1/Gtxfej1xIbFMqj9IKfzAarbIvK+6OefrauAb7+Fbt2s2kPXXOPtqCqOs5nKmhCUckHIqyEcPl902St/48/d0XcTGxbLDW1voE5NH1r+ShVx5ow1meyNN6BRI6tz+K67rIlm1YlPl65QyhedvHiS5fuWk7AvwWEyAGuE0D//8M8Kjky5In/Jidat4aab4PPP4fRpuPdeeP55uOwyb0fpWzQhKGWTK7ls+XWLvS9gw7EN9slhdWrU4WLWxSLH6CLyvqlwyYnDh63aQx07wooV0LWrd+PzVTrsVFV58dviCZ0Tit+zfoTOCSV+26X1l86ln+PTXz7l9s9vp+U/WhLzVgzTV09HEKZfP52f7v6JXx/9lX8N/he1a9QucF7tDPZdzkpOpKVpMiiO9iGoKs1RZ3BgQCBDrxzKrym/8v2R78nOzaZhYENubn8zse1jGdh+IM3qNHN4Lu0M9n1ZWZeGjBZmjLWGQXWnncqqWipuEfmul3clNiyW2LBYrgq+igA/vYNa2a1YYY0e2rnT8faQEEhKqtCQfJJ2KqtqZf/p/STsTSh2EfnEe52XkVaVS1KStaj9Z59Bu3bw17/CP/9p3SLKU5VLTriLJgRVJWRkZ7Dm8Br7+sF7Tu0BIMAvgOzc7CL7a2dw1ZCWBi+/DDNnWgvUvPCClRgCAyE6uuDCNjNmVJ8Zx2WlCUH5tOLu2x89f5SEvQkk7E1gxYEVXMy6SC3/WvRr248HezzIoLBBrD+2XmcGV0Ei8L//WVcCSUnwpz/B7NnW8NI81ankhLu4lBCMMQYIFpEjHo5HKTtHpaLv/PxOPtn+CYfOHeLnEz8DVh2gcV3HERsWS7/QfgUmh7W/zFqdVTuDq45du6x+gv/7P4iIgG++gX79vB1V1eByp7IxZpuIRHo4nlLRTuWqrbgO4b4hfe2LyHdu2hnrO4uqys6fh+eeg7lzoW5d6+f77oMAvc9Rau7oVN5sjOkhIj+5MS4AjDHPA7cAucDvwAQROe7u91G+L1dy2fzrZpbuWVpsh/DqCasrLijlVbm5sGiRtbj9iRNw551WddKmTb0dWdVTmoTQCxhrjDkEXAQMICLSxQ1xzBKRpwGMMQ8BzwD3uuG8qhI4m36Wr/d/TcLeBJbtW8bvF3/HYKjpX5PMnMwi+2uHcNWWv+TE5ZdDvXrWWgW9esGSJdCjh7cjrLpKkxAGeioIETmf72kdoPJOjlAlEhF+OfmLfUTQ94e/J0dyaBTYyJocFhbLwCsG8vWBr7VDuJopXHLit9+sx8SJ1jBSP62t4FE+MzHNGDMDGAecA/qJyEkn+00EJgK0adOm+yFnq2Arryo8Oujpvk/TtHZT+6igI+et8QlRzaOIbW9NDusV3KvI5DCdHVy9hIRYVwaOXtcJZe7jlpnKxpiuwHW2p2tEZGspjl0BNHewaaqIfJ5vvyeBQBGZVtI5tVPZNzkqF5Gnbs263NTuJvuaAa3qt/JChMoXffcd9O3reJuWnHCvcncqG2MmAXcDn9leWmSMWSAi81w5XkT6u/hW8UACUGJCUL4lIzuDbw99y31f3ucwGTSr04wjjxyhpr+TQjOqWjp2DKZMgY8+stYlyMkpuk8bXWG0QpSmD+FOoJeIXAQwxrwE/Ai4lBCKY4wJE5G9tqe3ALvKe05VMY6cO2LdBtpnTQ5zlAjynLx4UpOBssvIgFdftWYXZ2dbC9eEhMBf/lKwUqmWnKg4pUkIBsifu3Nsr7nDTGNMR6xhp4fQEUY+Kzs3mx+O/GDvC9j2+zbAGvkzoesEYsNiuW/pffY+gvx0IXmVJyEBHn7YGj00dCi88gq0bWttq1VLS054S2kSwrvAemPMf23PhxbozdYAACAASURBVALvuCMIERnhjvOo8nPUiXtTu5tYtncZCfsS+GrfV5zLOEeAXwDXtbmOWTfNIi4sjiubXGmfHPZixos6Okg5tG8fPPIIfPmltVjNV1/BgAEF99GSE97jUqdyXukKoCnQ2/byGhHZ4sHYSqSdyu7lqDPYz/iRK1ZvXvO6ze0jgm664ibq16pf7Ll0dJDKc/GiNZls9mxrrYJp06zyE87WLVCeVe5RRlq6omo7k3aGDvM6kJyWXGRbg1oN+Gb8N0Q1j8LP6EBw5ToR+M9/YPJkOHoUbrsNXnoJWrTwdmTVm0+XrlAVT0TY9vs2e1/AD0d+IEccDO8AzmecJ7pFdAVHqCq7bdusq4DVqyEqCj7+GK691ttRqeL4SukKVQEuZF5g5YGV9lFBR88fBaBb82480fsJ3t78NicunihynHYGK1fklZw4dMgqN3HxIjRsaM0wvvtua0ip8m2lKX89EWsEkPJBzu7Z7zm1x34V8O2hb8nMyaRezXrcdMVNTO87nUFhg2hZryUAnZp20s5gVSbx8dYf/bwVylJSrCqkM2bAvTpmsNLQPoQqwFFncIBfAI2DGtu/8Xdq0sm+fnDvNr2dzgfQzmBVFi1aWDWHCtOSE75J+xCqsMf/7/EiE8Kyc7M5n3Ge+YPmExsWS9tGbV0619jIsZoAlMtOnIAnnnCcDMBxXSLlu0rbh3CrMSYJ7UPwqqycrEuTw/YlcCzlmMP90rPTeaDnAxUcnaoOsrJg/nyYPt26TVS/vrWATWFacqJy8Yny16pkv134jeX7lpOwN4Gv939tnxzWJ6QPjQIbcSb9TJFjtDNYecLKldbooR074OabYc4c2LixYNlq0JITlVFpEsJhYCzQTkSeM8a0wapeqh3N5eDsnn2u5PLTsZ/sVwEbj1t9JS3qtmBk55HEhsXSv11/6teq77APQTuDlbsdOgSPPgqLF0O7dvD55zB4sFWJtGNHax8tOVG5laZT+Z9YtYZuEJFOxphGwNci4rX1iyp7p7KjP+Q1/WvSs1VPdifv5mTqSfyMH1cFX0Vs+1jiOsTR9fKuDtcP1s5g5SlpaTBrFsycaT1/6ilrollgoHfjUmXnjpnKm0Uk2hizRUS62V7bKiJd3Ryryyp7Qgh5NYTD54v2uvkZP8ZEjLGvHNa4dmMvRKeqOxHrKuCRR6yRQn/8o1V6QvsFKj93jDLKMsb4Y1ve0hjTFOuKQZVCSkYKKw9ak8McJQOwFppfNHxRBUem1CW7dsGkSfD11xAebvUb3HCDt6NSnlaahPAa8F+gmW25y5HA3zwSVRUiIpcmh+1L4Nukb8nKzaJezXrUDqhNanbR9QN0EXlVkfIvah8cDJGRViKoXdvqML7/fqhRw9tRqorgckIQkXhjzCbgRqwhp0NFZKfHIvNxxd2zT89OZ3XSavsM4f1n9gPQuWlnJvWaRFyHOK5pfQ2f7vhUO4OVVxVe1P7IEevRt69VlK5ZM+/GpypWqdZU9jXe6kNw1BkcGBDI6IjRJKcms/LAStKy0wgKCOKGtjfYZwiHNgx1eC7tDFbeEhpqjR4qTGcYV23l7lT2Rd5KCKFzQjl0zvFo27YN2xIXFkdsWCzXh15PUI2gCo5OKdecOgVNmjjepovaV23u6FSu9n678BvL9i5zmgwA9j+03+GwUKV8RU4O/Otf8LdiegB1JFH1VGJCMMb0AI6IyG+25+OAEVgT0qaLyGnPhug9Obk5/HT8J3tfwKZfNwHgb/wdrh0Q0iBEk4HyaWvWWIvYb90K/fpZy1c+/7zOMFYWV64Q/gX0BzDG9AFmAn8BooAFWKONqozTaaf5at9XJOxLYPm+5SSnJuNn/Lg6+Gpm3DCDuLA4tv2+jXu+vEc7g1WlcewYPPYYfPghtG5tdRiPHGndGmrdWmcYK4srCcE/31XAKGCBiCwGFhtjEj0XmmcU7sR94YYXiGwWydK9S0nYm8CPR38kV3JpUrsJN7e/mbiwOAZcMYDLgi6zn6Nrc2u2sHYGK1+XkWENHX3+ecjOtm4TPfEE1KlzaR9d1F7lKbFT2RizHYgSkWxjzC5gooh8l7dNRCLcFowxjwKzgaYiUnRx30JK26nsaHSQwSDWXDtiWsbYF5GPaRmDv58u8aQqr2XLrMlle/fCkCHw6qtWDSKlytOp/BHwrTEmGUgD1thO2B4458YAWwMDsIroecTUlVOLrBsgCI2DGrP9/u00r9vcU2+tVIXZv98qN/HFF9Chg5UYbr7Z21GpyqDEhCAiM4wxK4EWWMXs8i4p/LD6EtzlVeAx4HM3nrOAw+cc55rTaac1GahKLzUVXnzRKkRXowa89BI8/DDUdLw4nlJFuDTsVETWOXhtj7uCMMbcAhwTka0ljdIxxkzEWt+ZNqUcG9emQRuHQ0Z13QBVGeVf1L5JE6sY3alTVn/Ayy9Dy5bejlBVNq4MO11S3HYRGeLKGxljVmCtn1DYVOAprNtFJRKRBVijm4iJiSnVrLoZN87QUhGqSihcciI52Rox9PTT8Nxz3o1NVV6udCqfBI5g9SWsx6pjZCci35YrAGMigZVA3l/pYOA40DNv7oMzZZmprKUiVFXQpo1Vc6gwLTmhXFHm0hW2ktc3AWOALsBS4CMR+cVDgSYBMZ4YZaRUZZebCwsXwp13Ot6uJSeUK5wlBL+SDhSRHBFZLiLjgauAfcBqY8yDHohTKeXEhg1w9dVWMqhVy/E+WnJClUeJCQHAGFPLGDMcWAQ8wKW1EdxOREJduTpQqrr4/XcrCfTqZd0m+uADePttq8REflpyQpWXK53K7wMRQALwrIhs93hUSimysuCNN2DaNKvzeMoUq9O4Xj1ruzFackK5lyt9CLnARdvT/DsbQESkvodiK5H2IaiqatUqqwjdL7/AwIEwdy507OjtqFRVUeaZyiLi0m0lpVT5HT4MkyfDp59C27bwv/9ZZSe0iK6qCGX+Y2+M8TPG6AWqUm6Qng4vvABXXglffmnNJfjlF7jlFk0GquKUmBCMMfWNMU8aY+YbYwYYy1+AA8CfPB+iUlVLfLy1dKWfnzVv4K9/hfBwq38gLg527rR+DtLF9lQFc6V0xQfAGeBH4C6sWcUGGCoila78tVLeVHiG8eHDVhXSVq1gxQq48UbvxqeqN1cSQjsRiQQwxrwN/Aq0EZF0j0amVBU0dWrB1cny+PtrMlDe50ofQlbeDyKSAxzVZKBU6YlYhegccVSGQqmK5kpC6GqMOW97pABd8n42xpz3dIBKVQWJidCnj/PtOsNY+QJXSlf4i0h926OeiATk+9lrcxCUqgxOnYL774fu3WH3brjrrqKdxTrDWPkKnWOglAfk5MCbb1orli1YAA8+CHv2wFtvWY+QEGs4aUiItV1nGCtf4NICOUop133/vTXLeMsWuP56eO01iIy8tF0XtVe+Sq8QlHKTX3+F226D3r3h5En45BP45puCyUApX+bKxLQPbP+d5PlwlKp8MjOtdYw7dID//McaWrprF/zpTzrLWFUurtwy6m6MaQncYat8WnjFtNMeiUypSuCrr2DSJKvDePBga5LZFVd4OyqlysaVhPAm1hKX7YBNFEwIYntdqSovb1H7w4ehRQvrsWkThIVBQgIMGuTtCJUqH1eqnb4GvGaM+aeI3FcBMSnlcwqXnDh+3HqMGgXvved8BTOlKhOXRxmJyH3GmK7AdbaXvhORnz0TllK+5amnHJecWLdOk4GqOlweZWSMeQiIB5rZHvG2qqdKVWm//GLdJnLE2etKVUalGXZ6F9BLRJ4RkWeAq4C7PROWUt539iw8/DB07WqVqnZES06oqqQ0CcEAOfme51BoxFFZGWOmG2OOGWMSbY9Yd5xXqbLIzYV337WWrHztNWuB+9df10XtVdVXmpnK7wLrjTH/tT0fCrzjxlheFZHZbjyfUqW2YYM1y3jDBrj6amv0UPfu1rZ69XRRe1W1laZT+RVjzGqgt+2l20Vki0eiUqqC/f671XH8zjtw+eXWyKFbby14q0hLTqiqrlS1jERkM7DZQ7E8aIwZB2wEHhWRMx56H6XssrPhjTfgmWfg4kV49FHr5/pax1dVQxVWy8gYs8IYs93B4xbgn8AVQBTWimz/KOY8E40xG40xG0+ePFlB0auqaPVq6NbNmmncsyf8/DPMnq3JQFVfFVbtVET6u7KfMeYt4MtizrMAWAAQExMj7olOVWX5Zxi3aWMtav/991bdoZAQ+OwzGDpU6w4p5RPlr40xLUTkV9vTYcB2b8ajqo7CM4wPHbKuCGrUgOnT4bHHii5Yo1R1VeqEYIz5MzCES8NOvxCRj8oZx8vGmCis2khJwD3lPJ9SgPNF7Zs1g2nTKj4epXxZWa4Q+orI6LwnxpjXgXIlBBG5rTzHK+WMs0Xtjx+v2DiUqgzKkhBqGWPigCNAMKAX3MrnXLgAL7zgfLvOMFaqqLKMMrofaATEApcBD7g1IqXKQQQ+/NCaZfzSS3DddbqovVKuKnVCEJFUEVkkIjNFZBHwkAfiUqrUtm6Fvn2tyWMtWsAPP8B33+mi9kq5qiydyv/J/xRr7sBLbotIqVI6fRqefhrefBMuu8z6g3/HHeDvb23XGcZKuaYsfQjnReSuvCfGmH+6MR6lXJaTA2+/bY0kOnMG7r8fnnsOGjXydmRKVU5l6UMofPd1qjsCUao0fvjBml18770QHg5btsC8eZoMlCqPsvQhHCz0/LT7wlGqeL/+CuPGwbXXwokT8NFHVgmKLl28HZlSlZ/Lt4yMMXuBbcBWIBHYKiJJHopLKeBS2YlDh6BhQ0hPt9YrePJJqzpp3brejlCpqqM0Vwj/An4DTgGDgO3GmG3GmOeMMTU8Ep2q1vLKTuRNLjt7FjIz4e9/tx6aDJRyLyPiWn04Y0yiiETlex4F3AEcBkJEpMLXV46JiZGNGzdW9NuqCtKqleMZxSEhkJRU4eEoVWUYYzaJSEzh10szyuicMaaLiPwMICKJxpi+ItLVGOOpNRJUNZSaak0qc1ZeQhe2V8ozSpMQ7gUWGWMSsfoQOgJ5ZcNqujswVf2IWKWo//pX649+7dqOC9Np2QmlPMPlPgQR2Qn0BJYDzYB9wB+MMXWAjz0TnqouduyAm26CkSOhQQNr5NCCBbqwvVIVqTSjjC4DHsFKBjuA9/Mtc1lMGTGlnDt3Dp591ppDULeu9d9774WAfL+ZurC9UhWjNKOMPgZSgC+A2sBaY0xPj0SlqrzcXFi40CpCN2cO3H477NkDDz5YMBmMHWt1IOfmWv/VZKCU55SmD6GpiLxs+/lLY8wnwIfAVe4PS1VlGzfCX/4C69bBVVfBl19CTJHxDkq5T1ZWFkePHiU9Pd3boVSowMBAgoODqVHDtZkBpUkIp40xkSKyDUBEDhhjapd0kFJ5Tp60JpO98461YtnChXDbbeBXlgIqSpXC0aNHqVevHqGhoZhqsni2iHDq1CmOHj1K27ZtXTqmNAnhfmCxMWYN1ozlcGB/6cNU1UH+he1bt4brr4clS6yFax55BJ55xuo8VqoipKenV6tkAGCMoXHjxpw8edLlY0pMCMaYD4AtWCUrbgCuBzrZXnu0TJGqKq3wwvaHD8P770NEBPznP9Cpk3fjU9VTdUoGeUr7mV25QngX6ArcZvtvA6xRRjWBwcCnpQtRVXXOFrY/f16TgVK+rMS7tyLyjYi8KiITRKQb0AF4Esibl6CUXUaG84Xtjxyp2FiUKqv4bfGEzgnF71k/QueEEr8tvtzn9Pf3Jyoqyv6YOXNmsfsvXLiQBx98sNzvWxqlXiBHRLKx+hC2AYvcFYgx5i9Y6zPnAEtF5DF3nVtVjKVLYdIk59t1hrGqDOK3xTPxi4mkZlmXuYfOHWLiFxMBGBtZ9nHPQUFBJCYmuiVGR7KzswkIKMuaZ5f4xPgOY0w/4Bagq4iEA7O9HJIqhX374A9/sB41asDjj+sMY1V5TV051Z4M8qRmpTJ1pWfWAgsNDSU5ORmAjRs3cv311xfZ5+TJk4wYMYIePXrQo0cPvv/+ewCmT5/ObbfdxrXXXsttt91W7ljKl07c5z5gpohkAIjI716OR7ngwgWrDPU//gG1asHs2db8gpo1ITJSZxiryunwOcfVE5297qq0tDSiouwFo3nyyScZNWqUS8dOmjSJRx55hN69e3P48GEGDhzIzp07AdixYwdr164lKCioXPGB7ySEDsB1xpgZQDowWUR+crSjMWYiMBGgjd6D8AoR+OQTmDwZjh2zVjCbORNatLi0jy5sryqrNg3acOhc0Y6wNg3K9/emPLeMVqxYwY4dO+zPz58/z4ULFwAYMmSIW5IBVGBCMMasAJo72DTVFsdlWLOeewD/Mca0EweLNYjIAmABWOsheC5i5cjPP8NDD8G330J0tDWM9JprvB2VUu4z48YZBfoQAGrXqM2MGz1zzzMgIIDc3FwApzOpc3NzWbduHYGBgUW21alTx22xVFgfgoj0F5EIB4/PgaPAZ2LZAOQCTSoqNlWyM2es20HdusH27fCvf8GGDZoMVNUzNnIsCwYvIKRBCAZDSIMQFgxeUK4O5eKEhoayadMmABYvXuxwnwEDBjBv3jz7c091TvtEpzLwP6AfgDGmA9Ych2SvRqQAyMmBt96CDh3gjTfgvvusInQTJ4K/v7ejU8ozxkaOJenhJHKn5ZL0cJJbkkFeH0Le44knngBg2rRpTJo0iZiYGPyd/E/12muvsXHjRrp06ULnzp158803yx2PIy4voelJxpiawL+BKCATqw/hm5KO0yU03S9/yYnLL4egIDh4EK67zipN3bWrtyNUqvR27txJp2o6K9LRZ3fHEpoeIyKZwK3ejqO6K1xy4rffwBi4/36YP9/6WSlVdfnKLSPlA556qmjJCRFrwpkmA6WqPk0ICoAVK5wvXq+L2itVPWhCqOaSkmDECGs9Y2ez3nW6h1LVgyaEaiotzVrLuFMnWL7cmkmsi9orVb35RKeyqjgi8L//wV//al0djBoFs2ZZi9iAVXZCS04oVT3pFUI1smsXDBwIw4dD3bqwahV8/PGlZAC6qL1SnpSWlkbfvn3JyclhypQphIeHM2XKFN58803ef/99wCp7ffz4cfsxo0ePZu/evRUSn14hVAPnz8Nzz8HcuVYieO01a4JZOSvlKqVK6d///jfDhw/H39+fBQsWcPr06SKT0RYuXEhERAQtW7YE4L777uPll1/mrbfe8nh8+iehCsvNhUWLrHLUJ07AnXda1UmbNvV2ZEp5z8PLHybxN/eWfohqHsWcm+eUuF98fDwffvghQ4YM4cKFC3Tv3p0nn3ySnTt3UrduXUJDQ9m4cSNjx44lKCiIH3/8keuuu44JEya4Zb2DkugtoyoiPh5CQ8HPz/rvCy9A794wfrzVF7B+vVWCQpOBUt6RmZnJgQMHCA0NZcmSJfbqp/lLYI8cOZKYmBji4+NJTEwkKCgIPz8/2rdvz9atWz0eo14hVAGFZxgfOgRPPw3168O//20lBT9N/UoBuPRN3hOSk5Np2LBhmY5t1qwZx48fp3v37m6OqiBNCFWAs0XtGzSA22+v+HiUUkUFBQU5LW9dkvT0dLeteVAc/d5YBThb1P7o0YqNQynlXKNGjcjJySkxKdSrV4+UlJQCr+3Zs4eIiAhPhgdoQqjUjh2DP//Z+XadYayUbxkwYABr164tdp8JEyZw7733EhUVRVpaGidOnCAoKIjmzR2tL+ZemhAqoYwMa8nKjh3hs89g6FCrTHV+OsNYKd/zwAMP8N577wHYl8AEmD59OpMnTwZgxIgR7N69296p/OGHH3LPPfdUSHyaECqZhARrAfsnn4T+/WHHDvjvf60RRCEhVlXSkBCrDIVOKlPKt0RHR9OvXz9ycnJcPqZhw4aMHz/eg1Fdop3KlcS+ffDII/Dll9bqZcuXW7OO8+ii9kpVDnfccUep9r+9AkeG6BWCj7t40RpFFB4Oq1fDyy/Dtm0Fk4FSSrmDXiH4KBH4z39g8mRrtNCtt8JLL4FtNrtSSrmdXiF4WeEZxvHx1hXADTfA6NHQpAmsXQsffKDJQCnlWXqF4EWOZhhPmGDVIGrYEP75T7j7bihU+0oppTxCrxC8yNEM4+xsa8jonj1w772aDJTyBkdX7uVVt27dAs8XLlzIgw8+WP4Tu5FPXCEYYz4BOtqeNgTOikiUF0OqEM7WKr54ERo3rthYlFIWR1fuEydaP1f1kXw+cYUgIqNEJMqWBBYDn3k7Jk87caLocpV5dIaxUt7j6Mo9NdV63VNOnjzJiBEj6NGjBz169OD7778HrAlrs2fPtu8XERFBUlISSUlJdOrUibvvvpvw8HAGDBhAWlpauePwiYSQxxhjgD8BH3k7Fk/JyoJXX7XmEqSnF12kRmcYK+Vdzq7cnb3uqrS0NKKiouyPZ555xr5t0qRJPPLII/z0008sXryYu+66q8Tz7d27lwceeIBffvmFhg0bsnjx4vIFiI/cMsrnOuCEiDhdL84YMxGYCNCmkn2VXrkSHnrIml18880wZw5s3KhrGCvlS9q0cVwwsrx/bvLWP8izcOFCNm7cCMCKFSvYsWOHfdv58+cLlLZwpG3btkRFWXfWu3fvTlJSUvkCpAITgjFmBeCoOtNUEfnc9vMYSrg6EJEFwAKAmJgYcWuQHnLoEDz6KCxeDO3aweefw+DBVpmJjh01ASjlS2bMKNiHAJ6/cs/NzWXdunUEBgYWeD0gIIDc3Fz78/yVUmvVqmX/2d/fv3LdMhKR/iIS4eDxOYAxJgAYDnxSUTF5WlqatZZxp05WDaLnn4dffoEhQ6xkoJTyPWPHWrXAKrI22IABA5g3b579ed6VRGhoKJs3bwZg8+bNHDx40HNB4Ft9CP2BXSJS6av4i1hXAeHhMG0a/OEPsGsX/O1vUOgLgFLKB40dC0lJ1pygpCTPX8W/9tprbNy4kS5dutC5c2fefPNNwKp8evr0acLDw5k/fz4dOnTwaBxGxDfuuhhjFgLrRORNV4+JiYmRvHtwvmL3bpg0Cb76ykoIr71mzTpWSnnPzp076dSpk7fD8ApHn90Ys0lEYgrv6zNXCCIyoTTJwNsKT1x5+2147DGrNPW6dVaH8ZYtmgyUUpWHr40yqhScTVwRgTvugBdfhGbNvBujUkqVls9cIVQmjiauiEDz5vDOO9UvGRSekl8aaWlp9O3bt9gFQzIzM+nTpw/Z2dllfh9Hjhw5Qr9+/ejcuTPh4eHMnTvX4X7Lly+nY8eOtG/fnpkzZ7p0/NmzZxk5ciRXXnklnTp14scff3R4bk+1Xd6EprNnz/LGG2/YX/dUW4LzdiosJyeHbt268Yc//MGlbaGhoURGRhIVFUVMTJG7HMqNNCGUgbNF7U+cqNg4qoJ///vfDB8+HP9iijbVrFmTG2+8kU8+ce8AtICAAP7xj3+wY8cO1q1bx+uvv15gLDhYf6AeeOABli1bxo4dO/joo4/s+xR3/KRJk7j55pvZtWsXW7du9cj9a1farnBCcLUtV69ezYQJE1yOpbh2Kmzu3LlO28PZtlWrVpGYmIiv9RlWNZoQSiEnx6pA6uek1SrZPDm3e+WVV4iIiCAiIoI5c+bYX3/++efp2LEjvXv3ZsyYMQWm4sfHx3PLLbfYn2/dupU+ffrQuXNn/Pz8MMbwzDPPMHToUOLdUWEsnxYtWhAdHQ1AvXr16NSpE8eOHSuwz4YNG2jfvj3t2rWjZs2ajB49ms8//7zY48+dO8d3333HnXfeCVh/hBs2bFhsLM7aDpy3X+G2mzFjBh06dKB3797s3r0bgCeeeIL9+/cTFRXFlClTADzSlsW1U35Hjx5l6dKlDmfiFrdNVQztQ3DR2rXwl79AYiJ07gwHDlilJ/JU95ITmzZt4t1332X9+vWICL169aJv375kZ2ezePFitm7dSlZWFtHR0XTv3h2wbl8cOHCA0NBQwJp0M2rUKN5//3169uzJ008/TXp6Os8++yy5ubn89NNPDt/7uuuuIyUlpcjrs2fPpn///i7Fn5SUxJYtW+jVq1eB148dO0br1q3tz4ODg1m/fn2xxx84cICmTZty++23s3XrVrp3787cuXOpU6dOqdquW7du9lIGhduvcNtt2rSJjz/+mMTERLKzs+37zZw5k+3btxeYIRsREeG0LcvK1XZ6+OGHefnllx3+eznbZoxhwIABGGO45557mJhXaU65nV4hlOD4cWu1suuug1OnrFXMtm+3RhXpovaXrF27lmHDhlGnTh3q1q3L8OHDWbNmDd9//z233HILgYGB1KtXj8GDB9uPSU5OLvDNecWKFURHR9OzZ08AunTpwunTpzHG4O/vT82aNR3+IVmzZg2JiYlFHq4mgwsXLjBixAjmzJlD/fr1S/3ZCx+fnZ3N5s2bue+++9iyZQt16tQp9p66s7YDnLZf4bZbs2YNw4YNo3bt2tSvX58hQ4Y4fb/i2rJXr15ERUVx1113sWTJEnvdna+++qrU7VLYl19+SbNmzexfCFzdtnbtWjZv3syyZct4/fXX+e6778odi7fk7/eZMmUK4eHhTJkyhTfffJP3338fsEpaHD9+3H7M6NGj2bvXaTUft9IrBCcyM62ho88/bxWk+9vf4IknIO9Lni5qX35BQUEFpuJv376dyMhI+/PNmzfbb8kAZGRkFJnaD+W7QsjKymLEiBGMHTuW4cOHF9neqlUrjhw5Yn9+9OhRWrVqVezxwcHBBAcH2682Ro4cWWxCKIvCbVdaztoy71v96tWrWbhwIQsXLnTpfCW1E1jJbcmSJSQkJJCens758+e59dZbWbRoUbHb8s7TrFkzhg0bxoYNG+jTp08ZP7l35e/3WbBgAadPny7SB7Rw4UIiIiJoaVsi8b777uPll1/mrbfe8nyAIlJpH927dxdPWLZMpEMHERAZMkRk/36PvE2VUadOHdm0aZNEO0hvHAAADHlJREFURkbKxYsX5cKFCxIeHi6bN2+WDRs2SLdu3SQtLU1SUlIkLCxMZs2aZT82ODhY0tLSRERkwYIFMnr0aBER2b17t3To0EGSk5NFRCQ5OVk6duzo1rhzc3Pltttuk0mTJjndJysrS9q2bSsHDhyQjIwM6dKli2zfvr3E43v37i27du0SEZFp06bJ5MmTHZ6/uLYTkWLbL3/b5Z0jNTVVzp8/L+3bt5dZs2ZJcnKytGnTpsB7utKWq1atkvHjxxe7j6vt5Oz8cXFxJW67cOGCnD9/3v7z1VdfLcuWLXM5rjw7duyw/zxpkkjfvu59FPMrVMDVV18tBw8elMGDB4ufn5907dpVPv74Y5k2bZrMmjVLPv30U6lTp4506NBBunbtKqmpqZKTkyOhoaGSlZVV6s9d+LPnATaKg7+pessonwMH4JZbYNAg6/myZVYJinbtvBtXZRAdHc2ECRPo2bMnvXr14q677qJbt2706NGDIUOG0KVLFwYNGkRkZCQNGjSwHzdgwADWrl0LwJgxY7hw4QIRERFMnDiRjz76iMa2lYJWrVpFXFycW2P+/vvv+eCDD/jmm2/st0YSEhIAiI2N5fjx4wQEBDB//nwGDhxIp06d+NOf/kR4eHiJx8+bN4+xY8fSpUsXEhMTeeqpp0rddkCx7Ze/7aKjoxk1ahRdu3Zl0KBB9OjRA4DGjRtz7bXXEhERYe9U9kRbFtdOeW1ZFidOnKB379507dqVnj17EhcXx8033+zO0CtM/n6fJUuW2Kufjho1yr7PyJEjiYmJIT4+nsTERIKCgvDz86N9+/Zs3brV80E6yhKV5VGWK4RFi0RCQkSMsf67aJHIxYsif/ubSK1aInXrirz0kkhGRqlPrZxISUkREZGLFy9K9+7dZdOmTfZtmzZtkltvvbXEcwwbNkx2797tsRh9mbP2c7XtCquObenoW3JFO3bsWIErszp16th/zrtCEBHp27ev/PTTTwWO/fOf/yxLliwp0/uW5gqhWvUhOJphfMcdVu2hU6esPoGXXwbbrTvlJhMnTmTHjh2kp6czfvz4Av0C0dHR9OvXj5ycHKfj6TMzMxk6dKjHC3v5Kmft50rbFVbd29KbytPvk56eTlBQkJsjKqpaJQRHM4wzM+H8eVizBnr39k5cVd2HH35Y7PY77rij2O01a9Zk3Lhx7gypUimu/Upqu8Kqe1t6U6NGjcjJySE9Pd1hh36eevXqFRkksWfPHiIiIjwdYvXqQ3C2BF5WliYDpZTn5e/3cWbChAnce++9REVFkZaWxokTJwgKCqJ5c0fri7lXtbpCcLY0XkhIxceilKp+HnjgAV599VX69+9fYInM6dOn238eMWIEI0aMsD9/8803ueeeeyokvmp1hTBjhjWjOL/qPsNYKVVx8vf7uKphw4aMHz/eg1FdUq0SgjeWxlNKqfzuuOMOlwcBANx+++0EBFTMzZxqdcsIdIaxUtWViGCq2WLmUsoVMavVFYJSqnoKDAzk1KlTpf4DWZmJCKdOnSp2RFNh1e4KQSlV/QQHB3P06FFOnjzp7VAqVGBgIMHBwS7vrwlBKVXl1ahRg7Zt23o7DJ+nt4yUUkoBmhCUUkrZaEJQSikFgKnMve7GmJOAkyXvS9QESHZjON5SVT4H6GfxVfpZfFN5PkuIiDQt/GKlTgjlYYzZKCIx3o6jvKrK5wD9LL5KP4tv8sRn0VtGSimlAE0ISimlbKpzQljg7QDcpKp8DtDP4qv0s/gmt3+WatuHoJRSqqDqfIWglFIqH00ISimlgGqYEIwxjxpjxBjTxMn28caYvbZHxaxKUQrGmOeNMT8bYxKNMV8bY1o62S/Htk+iMeb/t3e2MXZVVRh+3iAfYUA6SIqWKsMY0BJqoTZN5CMhKUViiAXRaASNtok/0BAkomJR/kBC0x8mxmgNSEA7iSbFolWglBKC1GCBwjClXxlKkzIJmlSjVmK19vXH3iOXO+feOfPRc++ZWU+yc8/eZ9191pp171l3731mr99UrWcZJmBLV/sEQNJaSXuyPRslzWkhd0DSULb5har1LMMEbLlG0l5Jw5K+XbWeZZD0GUmvSjomqeUjmjXxS1lbJu8X27OmAO8HNpP+me2sgvNnAvvza28+7u203k06vrvh+BZgXQu5w53WdTpsqYNPsp5XA+/Kx2uANS3kDhR99rqplLEFOAF4DegHTgIGgQs7rXuBnguADwFPA0vayNXBL+PaMlW/zLYRwveBbwKtVtI/Dmyx/RfbfwW2ANdUpVwZbP+9odpDa1u6npK2dL1PAGw/Yftorj4HlN9zuMsoactSYNj2ftv/Bn4BrKhKx7LY3m17b6f1mA5K2jIlv8yagCBpBTBie7CN2DnAwYb6G7mtq5B0j6SDwI3A91qInSLpBUnPSbquQvUmRAlbauGTJlYCj7U4Z+AJSS9K+kqFOk2WVrbU0S/tqJtfWjElv8yofAiSngTeW3BqNfAd0lC462lnh+1f214NrJZ0B/A14K4C2XNtj0jqB56SNGT7teOodiHTZEtXMJ4tWWY1cBQYaNHN5dkvc4EtkvbYfub4aNyaabKlKyhjSwlq45fjyYwKCLavKmqXtBA4DxjMOVXnAzskLbX9ZoPoCHBlQ30+ab6uUlrZUcAA8CgFN1HbI/l1v6SngUtIc4uVMg22dIVPYHxbJH0JuBZY5jyhW9DHqF/+LGkjaYhf+Y1nGmwZIa3JjTI/t1XOBD5j7fqohV9KMCW/zIopI9tDtufa7rPdRxpGLW4KBpAWnK+W1CuplzSi2Fyxum2RdH5DdQWwp0CmV9LJ+fgs4DJgVzUalqeMLdTAJ5Ce7CCtT33S9lstZHoknT56TLJlZ3ValqOMLcDzwPmSzpN0EvA5oCufZhuPuvilJFPzS6dXzjtRaHiiAFgC3N9wbiUwnMuXO61rge4Pkz6srwCbgHOa7QAuBYZITxgMAas6rfdkbamDT7KOw6S525dzWZfb5wGP5uP+7JNB4FXSNEDHdZ+MLbn+CWAfaeTZrbZcT/oBeAT4E7C5xn4Z15ap+iW2rgiCIAiAWTJlFARBEIxPBIQgCIIAiIAQBEEQZCIgBEEQBEAEhCAIgiATASEIgiAAIiAEQRAEmQgIQa2QtEzS+g7rcIuk3ZIGmtrnSLq5qe0PFek05tpBMFEiIAR1YxHpv2c7yc3Acts3NrXPyef+j+1LK9JpzLWDYKJEQAjqxiLSJoUflvRUznD1ZN6zCUkLJD2Ts33dLml4sheSdJuknbncmtvWkbY6eEzS15veci/wwazT2ix/OL/25SxkD0raJ2lA0lWStillglua5W6StD338RNJJxTo1SPpd5IGs26fbXHtMX016DGQRzkbJJ062b9RMMPo9P4cUaJMpJBGB2eT9py5OLd9C7iHtHvvDuCS3P5j4JFJXuejpH2geoDT8vVG+z1Acca9PmBnU9vhhnNHgYWkH2IvAg8AIm3s9wgpI9Ym4MT8nh8BXyy4zg3AfQ31M5qv3aqvLGfgstz+APCNTvs1SneUGCEEtUHSiaSb35XAs7ZHp452AXOBTwGDtl9qaB/M7+2X9FNJG3K9R9JDku6T1Dz1A3A5sNH2P20fBn4FXDFFE1532nn3GCnAbLVtUuDpA5aRAtHzkl7O9f6CfoaA5ZLWSLrC9t8KZNr1ddD2tny8PtsaBDMrH0Iw41kA7AYuJN0UR1lIuvl/hHeuL1wEPA4pLwSwajQgkILHBtubJP2SapLAHGk4PtZQP0b6Lgp4yPYd7TqxvU/SYtKulndL2gr8rEmssC9JfYxNVRo7XAZArCEE9WIR6Rf/CCkokDPCfYF0QzwEXJDbLwZuyvJFzOftVIP/LTj/e+A6SafmPfKvz23t+AdwelljCtgKfDpn7ULSmZLObRaSNA94y/Z6YC2wuODa7fr6gKSP5ePPA89OQedgBhEBIagTo08Y/RyYJ2mIlER8pe1DuX1Jbl8FHMgjgyLe4O3k8WO+B7Z3AA8C24E/kvIzvNQs1/SeQ8C2vNC7dqLG2d4F3EnK7fsKsAV4X4HoQmB7ngq6C7i7+drj9LUX+Kqk3UAvaa0lCCIfQjBzkHRanu9H0u3AGbbvzPX3kBaelwP3Az8Afgj8i7Qe0dV5g6eLPGX0W9sXdViVoAuJgBDMGCR9l5Qy8D/ANuA220fav2t2EQEhaEcEhCAIggCINYQgCIIgEwEhCIIgACIgBEEQBJkICEEQBAEQASEIgiDIREAIgiAIgAgIQRAEQeZ/HJw7LHVeFlsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\"\"\"\n",
    "Compare Euler and Heun methods with different time steps\n",
    "\"\"\"\n",
    "#\n",
    "# No data this time. We'll just compare the analytical and computer\n",
    "# solutions to the problem.\n",
    "#\n",
    "\n",
    "m = 0.0038  # kg (mass of cup plus 6 paper clips. See proj. 1)\n",
    "g = 9.8     # m/s^2, gravitational field strength at Earth's surface\n",
    "b = 0.00375 # effective drag coefficient, measured experimentally\n",
    "\n",
    "s = 0.0     # s is the 'state', in this case, it's just the velocity.\n",
    "dt = 0.03   # pick a reasonable time step\n",
    "t = 0.0     # start at t=0.0\n",
    "tf = 1.0    # stop in 10 sec\n",
    "\n",
    "def derivs_drag(s,t):\n",
    "    \"\"\"\n",
    "    The derivs_drag function computes the time derivative of the \"state\"\n",
    "    for this particular problem. In this case the \"state\" is just the velocity\n",
    "    of an object falling with drag.\n",
    "    \"\"\"\n",
    "    dsdt = b*s**2/m - g\n",
    "    return dsdt\n",
    "\n",
    "def EulerStep(s, t, derivs, dt):\n",
    "    \"\"\"\n",
    "    EulerStep is a generic implimentation of the Euler Algorithm.\n",
    "    Notice that it knows *nothing* about the state, nor the physics\n",
    "    of the problem (which is handled by the 'derivs' function\n",
    "    that is passed in by the programmer).\n",
    "    \"\"\"\n",
    "    return s + derivs(s,t)*dt\n",
    "\n",
    "def HeunStep(s, t, derivs, dt):\n",
    "    \"\"\"\n",
    "    HeunStep is a generic implimentation of the Heun Method.\n",
    "    Notice that it knows *nothing* about the state, nor the physics\n",
    "    of the problem (which is handled by the 'derivs' function\n",
    "    that is passed in by the programmer).\n",
    "    \"\"\"\n",
    "    f1=derivs(s,t)\n",
    "    f2=derivs(s+f1*dt,t+dt)\n",
    "    return s + 0.5*(f1+f2)*dt\n",
    "\n",
    "v0=np.sqrt(m*g/b)\n",
    "tau=v0/(2*g)\n",
    "\n",
    "def VAnalytical(t):\n",
    "    return v0*(1-np.exp(t/tau))/(1+np.exp(t/tau))\n",
    "\n",
    "\"\"\"\n",
    "OK.. ready to run!\n",
    "\"\"\"\n",
    "dt = 0.1\n",
    "dtlist = []\n",
    "err_heunlist=[]\n",
    "err_eulerlist=[]\n",
    "\n",
    "while dt > 1e-4:\n",
    "    t = 0.0\n",
    "\n",
    "    s = 0.0\n",
    "    s_euler = 0.0\n",
    "    s_exact = VAnalytical(t)\n",
    "    \n",
    "    vlst = [s]\n",
    "    tlst = [t]\n",
    "    \n",
    "    velist = [s_exact]\n",
    "    veulerlst = [s_euler]\n",
    "    \n",
    "    while t < tf:\n",
    "    \n",
    "        # do a Heun Step and an Euler Step\n",
    "    \n",
    "        s = HeunStep(s, t, derivs_drag, dt)\n",
    "        s_euler = EulerStep(s_euler, t, derivs_drag, dt)\n",
    "        t += dt\n",
    "        s_exact = VAnalytical(t)\n",
    "        \n",
    "        vlst.append(s)    \n",
    "        velist.append(s_exact)\n",
    "        veulerlst.append(s_euler)\n",
    "        tlst.append(t)\n",
    "    \n",
    "    #\n",
    "    # now compute the rms error for Heun and Euler methods\n",
    "    #\n",
    "    \n",
    "    diff_heun = np.array(vlst)-np.array(velist)\n",
    "    diff_euler = np.array(veulerlst)-np.array(velist)\n",
    "    \n",
    "    heunErr = np.sqrt((diff_heun**2).sum()/len(tlst))\n",
    "    eulerErr = np.sqrt((diff_euler**2).sum()/len(tlst))\n",
    "    err_heunlist.append(heunErr)\n",
    "    err_eulerlist.append(eulerErr)\n",
    "    dtlist.append(dt)\n",
    "    dt = dt/1.5\n",
    "\n",
    "coefs_euler = np.polyfit(np.log10(dtlist), np.log10(err_eulerlist), 1)\n",
    "coefs_heun = np.polyfit(np.log10(dtlist), np.log10(err_heunlist), 1)\n",
    "\n",
    "plt.plot(np.log10(dtlist), np.log10(err_eulerlist), \"go\", np.log10(dtlist), coefs_euler[1] + np.log10(dtlist)*coefs_euler[0], \"g\")\n",
    "plt.plot(np.log10(dtlist), np.log10(err_heunlist), \"bo\", np.log10(dtlist), coefs_heun[1] + np.log10(dtlist)*coefs_heun[0], \"b\")\n",
    "plt.xlabel(\"$log_{10}$ of time step\")\n",
    "plt.ylabel(\"$log_{10}$ of RMS error\")\n",
    "plt.title(\"RMS error log($\\sigma$) for various time steps\")\n",
    "plt.legend(['Euler', '(fit)','Heun', '(fit)'], loc=4)\n",
    "\n",
    "plt.text(-3.2, -7.5, \"log($\\sigma$) = %6.3f log(dt) + %6.3f\" % (coefs_heun[0], coefs_heun[1]))\n",
    "plt.text(-3.5, -1.5, \"log($\\sigma$) = %6.3f  log(dt) + %6.3f\" % (coefs_euler[0], coefs_euler[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Coupled Equations\n",
    "=================\n",
    "\n",
    "Suppose in addition to the velocity of the muffin cup, we also want to track the position. How could we do it? The first thing to notice is that the position and the velocity are related by yet another first order differential equation:\n",
    "\n",
    "$$\\frac{dy}{dt} = v$$\n",
    "\n",
    "So that we can introduce the position $y$ by adding a second first order equation to our system.\n",
    "\n",
    "Let's re-work the Heun method example from above to do this. Now we'll define the *state* as not just the velocity, but include both the velocity *and* the position.\n",
    "\n",
    "`s = array([y, v])`\n",
    "\n",
    "The primary changes then relate to the `derivs_drag` method and how we set initial state and extract values from the state to make graphs and comparisons. Examine the following code. Notice how the state is now an array that contains the $y$ and $v_y$ values at each time. Notice also that the `HeunStep` function is identical to the one used earlier! This means that we can use the same function to handle any number of coupled first order equations. Neat!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fce1b35a3a0>]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcAAAAEWCAYAAADxQkdBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXhU5fn/8fedhIRFXEARERA3FBRFxSVWJYoKbpUq+nOre3GptrXV1g3FrWq1VuvXVmjrvoFaERUXBCIuAcUqgiAIgiwuoIiCLCHk/v1xTjDGLDPJZM6Zmc/ruubKzDnPPHOfoHPn2c5j7o6IiEiuyYs6ABERkSgoAYqISE5SAhQRkZykBCgiIjlJCVBERHKSEqCIiOQkJUCJBTMrNbNzw+dnmtkb1c6tNLPtoouufmbWysyeM7NvzezJ8NiNZvaVmX1hZl3Da8iPOtZEmNl5ZnZn1HEAmNnFZnZr1HFIdlIClJQzs/lmtjr80q96dGpsfe6+kbt/ksoYAcxsqJm5mf22xvHfhseHJljVIGBLoL27n2BmXYE/AD3dvaO7LwivYX0j49zYzO40swXh73Ju+HrzxtTXwGcVAlcDt4Wvu4W/i4Ia5R4wsxtT/fm1+Bdwqpl1SMNnSY5RApTmckz4pV/1+CzqgOowGzi9xrEzwuOJ2gaY7e4V4euuwNfuvqSpwYUJaRywCzAA2BgoBr4G9mlq/bU4FvjI3Rc3Q91Jc/c1wIv89N9IpMmUACUtzGwzM3vezJaa2Tfh884JvtfNbIfw+QNmdo+ZvWBmK8xsspltX63s4WY2K+yO/IeZvVbVtVqHd4DWZrZL+P5dgJbh8ao6f9QlWz0mM7sOuAb4f2Hr7DxgLNApfP1AzVZU2N17g5m9GV7DK/W05k4nSKi/cPcZ7l7p7kvc/QZ3H1Pz91Ptd3Rj+LzEzBaZ2ZVhl+x8Mzu1nt/HEcBr9ZyvlZntZ2ZvmdlyM5tqZiXVzs03s0OrvR5qZo+Ez6t+N2eELdyvzOyqGtWXAkclG5NIQ5QAJV3ygPsJWktdgdXA/zWyrpOA64DNgDnATQBhEnkKuAJoD8wC9k+gvof5oYVxRvg6Ie5+LfBnYETY0h1GkEQ+C1+fWcdbTwHOAjoAhcCldZQ7FHjJ3VcmGlMtOgKbA1sTXN9wM9upjrK9CH5vCTOzrYEXgBuBdgTX8rSZbZFENQcAOwH9gGvMrEe1czOB3ZOJSSQRSoDSXEaFrYHlZjbK3b9296fdfZW7ryBIWn0bWfcz7v522OX4KNA7PH4k8KG7/zc893fgiwTqewQ42cxaECTXRxoZVzLud/fZ7r4aGMkP11BTe+DzFHzeEHdf6+6vESSrE+sotymwopbjX1X791xOkMCrnAaMcfcxYQt1LDCF4N8jUde5+2p3nwpM5ccJbwWwSRJ1iSSkoOEiIo0y0N1frXphZq2BvxGMY20WHm5rZvmNmBxSPamtAjYKn3cCFladcHc3s0UNVebuC8xsDkFL7mN3X2hmSYaUtLquoaavga2a+FnfuPv31V5/SvC7qrUs0LaW45tXG+PEzB6odm4b4AQzO6basRbAhCRirO/30Rb4Nom6RBKiFqCkyx8Iurj2dfeNgYPC46nMNJ8DG8YVLchiCY0zAg+FMT5Uy7nvgdbV6u3YhBiT9SrQ38za1FNmFdXiI+jyrG6zGu/vCtQ1KekDoHuSMS4EHnb3Tas92rj7LeH5H/3+aomvIT0IWoUiKaUEKOnSlmDcb7mZtQOubYbPeAHoZWYDwwknvybxL9sRwOEE3ZE1TQV2MbPeZtYSGJqKYBP0MEGCedrMdjazPDNrH05qqepifB84xczyzWwAtXctX2dmhWZ2IHA08GQdnzemjvfX5xHgGDPrH8bQMpx8U/XHx/vASWbWwsz6ECwbSUZfgpmgIimlBCjpcifQCvgKmAS8lOoPcPevgBOAvxB0HfYkGItam8B7V7v7q+GYXM1zs4HrCVpjHwNv1CzTXNx9LcFEmI8IZpd+B7xNMKllcljst8AxwHLgVGBUjWq+IOja/IxgzPR8d/+ojo98DtjZkli36e4LCZZPXAksJUjYl/HD98sQYPswhuuAxxKtO/yD40jgwUTfI5Io04a4kq3MLA9YBJzq7smMR2WNcDnCI+6eaFcwZjaYYBH/75otsMRjuRjo4u5/jDoWyT6aBCNZxcz6E7SMVhO0QoygxSkJcvfhUcdQxd3vjjoGyV7qApVsUwzMJehqPYZgNupPujVFRNQFKiIiOUktQBERyUk5MQaYl5fnrVq1ijoMEZGMsmrVKnf3rG0o5UQCbNWqFd9//33DBUVEZAMzy+rx86zN7CIiIvVRAhQRkZykBCgiIjlJCVBERHKSEqCIiOSk2CVAMxtgZrPMbI6ZXV7L+SIzGxGen2xm3dIfpYiIZPr3dawSoJnlA/cARxDcyf9kM+tZo9g5BBt87kCwweqtzRXPtFHDKT2vP9NGxebWiCIisRC37+vGiFUCBPYB5rj7J+5eDjxBsM1Kdcfyw9YoTwH9rBm27542ajjbD7qAA4aPZfsTz6s7CZaVwc03Bz9FRHJHbL6vGytuCXBrgr3EqiwKj9Vaxt0rgG+B9jUrMrPBZjbFzKZUVFQkHcjXLz5N4XqnAKdFRfD6J8rKoF8/GDIk+FlXElSSFJHMVFD1PRo+Blc7l7Lv66hk7Z1gwi1dhgO0adMm6Tt+tz/ieMrvG4tX5LOOfNofcfxPC5WWQnk5rF8f/CwtheLiH5epSpLl5VBYCOPG/bRMVbnSUigpqf28iEj6Vbh7n6iDaC5xS4CLgS7VXncOj9VWZpGZFQCbEOz+nVK9Bg5m2pPw1jXv8/CHp/Fgr/1/WqikJEhqVcmtpOSnZVKZJEVE4iM239eNFbcu0HeAHc1sWzMrBE4CRtcoMxo4I3w+CBjvzbSnU6+Bg/n5y/9gSov9ueWWWgoUFwfJ6oYb6k5aVUkyPz+5JCkiEm+x+r5ujNjtB2hmRwJ3AvnAfe5+k5ldD0xx99Fm1hJ4GNgDWAac5O6f1FdnmzZtvCk3w77oIhg+HObMga5dG1FBQ92b6iYVkRgys1Xu3qae8yn/vk6n2CXA5tDUBLhgAeywA5x3Htx9dwoDqy5VSVJEJEUaSoCZLm5doLHUtSuccQb861/wxRfN9CHFxXDFFXUnNXWTioiklBJggi6/HNatg9tvjyiARMYSQUsuREQSpC7QJPzyl/Df/8Knn8Lmm6cgsGSpm1RE0khdoLLBFVfA6tVw550RBaBuUhGRlFECTELPnnD88cFEmOXLo46mFuomFRFJmLpAkzR1KvTuHSz9u/rqlFSZWuomFZEUUReo/Mjuu8Mxx8Df/gYrVkQdTS3UTSoikhAlwEa46ipYtgzuvTfqSBoh0W5SEZEspy7QRjr88KA7dP58aNUqpVU3v0TuKKO7zojkvGzvAlUCbKSJE6FvX/j73+Hii1NadfQ0TigiZH8CVBdoIx10UPC49VZYuzbqaFJM44QikgOUAJvg6qth8WI46eoXKVuYRUsKNE4oIjlAXaBN8NaCMg74WQG+cnNaXrI74896meIuWdJVqDFAkZynLlCp02uflmIH3QjLt2XteydQOr806pBSp6HlFKAF9SKS0eK2I3xGKelWQmGPG1iz1bsw8UoOuGNJ1CGljybKiEiGUwuwCYq7FDP+jHGcdtE8fNn2zH8jhxKAJsqISIbTGGAKuMMeewQ3yp4xI5g7kvXUAhTJehoDlAaZwZAhMHs2jBgRdTRpUlwcJL0bblDyE5GMpBZgilRWBvcJXb8epk3LkVZgIjSbVCRjqQUoCcnLC1qBM2fCU09FHU1MVHWTDhkS/NRsURGJESXAFDr+eOjRI+gVrKyMOpoY0EQZEYkxJcAUys8PGjsffgjPPBN1NDGgO8qISIxpDDDF1q8Pdo5v2RLeey/oGs1pGgMUyVjZPgYYmwRoZu2AEUA3YD5wort/U0u59cC08OUCd/95Q3WnMwECPPwwnH560AocODBtH5u5lCRFYkkJME3M7C/AMne/xcwuBzZz9z/VUm6lu2+UTN3pToAVFcFYYNu28O67wTIJqYPWE4rEVrYnwDh10B0LPBg+fxDI2LZTQUGwa/x778Hzz0cdTcxpooyIRCROCXBLd/88fP4FsGUd5Vqa2RQzm2RmdSZJMxsclptSUVGR8mAbcuqpsO22cP31wZ1ipA6aKCMiEUlrF6iZvQp0rOXUVcCD7r5ptbLfuPtmtdSxtbsvNrPtgPFAP3efW9/nprsLtMp//gPnngtjxsARR6T94zOHxgBFYinbu0DjNAY4Cyhx98/NbCug1N13auA9DwDPu3u9S8+jSoDl5dC9O3TsGHzHayywCZQkRdIu2xNgnLpARwNnhM/PAJ6tWcDMNjOzovD55sDPgBlpizBJhYXBlnqTJ8M5dzyRXbvGp5PuKCMizSBOCfAW4DAz+xg4NHyNmfUxs3+HZXoAU8xsKjABuMXdY5sAAXbuNwnbeBH339WVQx7spyTYGJooI5JxzKydmY01s4/Dnz8Z0grLrTez98PH6HTGGJsNcd39a6BfLcenAOeGz98CeqU5tCZ56/MJcOAieOEeymeXUDq/lOIu6sJLStVEmaqlEpooI5IJLgfGVVvadjnwk6VtwGp3753e0AKxGQNsTlGNAQKULSzjkPsGsObOqVibr3jjrXXs31UJMGkaAxRJu6aMASY6r6Mxa7tTRQkwDcoWlnH7/y3jv385itGj4ZhjIgtFRCRhTUyAy6tm9puZAd9Un+lfrVwF8D5QQTCsNaopMScVoxJgeqxbBzvvDBtvHNwdJufvEdoc1EoUSSkzK+eHW08CDHf34dXOR7K0LVViMwaY7Vq0gGuvhTPOgFGj4Ljjoo4oy+iWaiLNocLd+9R10t0PreucmX1pZltV6wJdUkcdi8Ofn5hZKbAHkJYEqHZIGp16Kuy0U5AItV9gimmmqEjcxH5pmxJgGuXnw9ChMH06jBwZdTRZRrdUE4mb2C9t0xhgmlVWwm67BTtGTJ8e3DhbUkRjgCIple13glECjMDTT8OgQfDgg8G+gSIicaQEmAXilgArK2GvvWDFCpg5M5ggI2miVqJIwrI9AWoMMAJ5ecE2SXPnwkMPRR1NDtE9RUWkGiXAiBx9NOy9N9xwQzBpUdJAM0VFpBolwIiYBa3ATz8N9g2UNNBMURGpRmOAEXKHAw4IkuCcOdCyZdQR5QCNAYokLNvHAJUAIzZ+fDAcdddd8JvfRB2NiMgPlACzQJwTIMDBB8MHH5bzm0fv4vCdD9B2SSISC9meADUGGAMn/no6y5YWct3tS+n3kDbNjVxZGdx8s2aJimQ5JcAYWL7lc7D9WPz1P7L2+5aUzi+NOqTcpaUSIjlDCTAGSrqVUHT4UFi9OXllf6SkW0nUIeUuLZUQyRlKgDFQ3KWYCVfcTq9DZlAw+TK6FWgMMDJaKiGSMzQJJkbmzIEePeDcc+Gf/4w6mhympRIiQPZPglECjJlf/xqGDYMZM6B796ijEZFclu0JUF2gMTNkSLAg/uqro45E6qRZoiJZITYJ0MxOMLMPzazSzPrUU26Amc0yszlmdnk6Y0yHjh3h97+HJ5+Ed96JOhr5Cc0SFckasUmAwHTgOGBiXQXMLB+4BzgC6AmcbGY90xNe+lx6KWy+OVx+eXC7NIkRzRIVyRqxSYDuPtPdZzVQbB9gjrt/4u7lwBPAsc0fXXptvHHQwBg/HsaOjToa+RHNEhXJGrFJgAnaGlhY7fWi8NhPmNlgM5tiZlMqKirSElwqnXcedOsGf/pTsIGuxERxMYwbF+xjNW6cZomKZLCCdH6Ymb0KdKzl1FXu/mwqP8vdhwPDIZgFmsq606GoCG68EU47DUaMgJNPjjoi2aC4WIlPJAukNQG6+6FNrGIx0KXa687hsax08slw221w1VVw/PFBj5uIiKRGpnWBvgPsaGbbmlkhcBIwOuKYmk1eHtxyC8ybF6wNlAyipRIisRebhfBm9gvgbmALYDnwvrv3N7NOwL/d/ciw3JHAnUA+cJ+739RQ3Zm0EL4m92C2/fTpMHcutG0bdUTSoKqlEuXlQbNdY4WSobQQPk3c/Rl37+zuRe6+pbv3D49/VpX8wtdj3L27u2+fSPLLdGZBK3DpUvjrX6OORhKipRIiGSE2CVDqts8+MGgQ3Hb7eq589i7tFxh3WiohkhGUADPEoIveY9Vq55Y/F2jT3LjTUgmRjJDWWaDSeJ/kvYTtNQl/5zzW7n0vpfNLKe6iL9bY0lIJkdhTCzBDlHQroejQm6FoJbxyB323KYk6JBGRjKYEmCGKuxQz/oIRHHnuFCo/Poxvpql1kfG0VEIkUrFZBtGcMnkZRE3l5dCrV/B82jQtjs9YWiohGUDLICRWCgvhjjtg9my4556oo5FG01IJkcgpAWagI4+E/v3huuuC9YGSgbRUQrJcJuzxqgSYgcyCVuDKlXDNNVFHI42ipRKS/WK/x6uWQWSonj3hwguDbtALLoDddos6IkmalkpIFnP3mQBmVl+xDXu8hmWr9nid0ewBohZgRhs6FDbdFC65RDvHi0izKKjaVzV8DE5x/Qnv8doc1ALMYO3aBeOAF18Mzz4LAwdGHZGkXFlZMEGmpEStRYlChbvXN36Xtj1em4MSYIY7/3z45z/h0kvhiCOCjXQlS2iphMRcpu/xqi7QDFdQAH/7W7BV0l13RR2NpJSWSkj2i3SP16QToJm1CWfuSEwcfjgcfTTceCN8+WXU0UjKaKmEZDAz+4WZLQKKgRfM7OXweCczGwPg7hXARcDLwExgpLt/mLYYG7oTjJnlEWTlU4G9gbVAEfAV8AIwzN3nNHOcTZJNd4Kpy+zZsOuuMOD4Lym+8D5KupXoZtnZQGOAEqFsvxNMIgnwNeBV4FlgurtXhsfbAQcDpwDPuPsjzRxro+VCAgQ4efBnPPHvjuSdtw9FnWcw7vRxSoIi0mjZngAT6QI91N1vcPcPqpIfgLsvc/en3f14YETzhSiJ2nHgE9D6KypfuIu169ZROr806pBERGKrwQTo7utgw21t2obPh5jZf81sz+plJFpH9CqmRf+rYeHPyH9/MCXdSqIOSdJBu0qINErCu0GY2QfuvpuZHQDcCNwGXOPu+zZngKmQK12gAG8tKOO0Yzvx1Sdb8/HsArbcMuqIpFlpqYQ0o7h3gZpZG2CNu69vzPuTmQVa9QFHAcPd/QVAm/HEzP5dixnz+DasWV3AH/4QdTTS7LRUQnKImeWZ2Slm9oKZLQE+Aj43sxlmdpuZ7ZBMfckkwMVmNgz4f8AYMytK8v2SJjvvDJdfDo8+Cq++GnU00qy0VEJyywRge+AKoKO7d3H3DsABwCTgVjM7LdHKEpkFWhxW3AoYAExz94/NbCugl7u/0rjrSJ9c6gKtsmbNjzfObdky2nikGWmphDSTuHWBmlmLhuacJFJmQ9kEEuA/gX2B2cBLwEvu/kWC8SbMzE4AhgI9gH3cfUod5eYDKwi6ZOu9T12VXEyAELT+DjsMhgyB66+POhoRyTRxS4BVwv0FrwK2IbilpwHu7knti5PMJJidCfZs6g9sQtAUfQl4s7EDkDXq7wFUAsOASxtIgH3c/atE687VBAjwy1/CiBHwwQdB16iISKJinABnAZcB0wjyBgDu/mlS9SSaAGt8eCuCRfBHAMWJtMKSqLsUJcCUWbIkSHy77QYTJgSb6UoOUjepNEKME+Ab7n5AU+tp1G4Q7r4aGBM+0s2BV8zMCW7DNry2QuG+VYMBCgtzd7Jqhw7wl7/Ar34FDz4IZ54ZdUSSdloqIdnnWjP7NzCO4PacALj7f5OpJOEEWEufa9UHJtznmqK9ow5w98Vm1gEYa2YfufvEmoXCxDgcghZgojFmo7PPhgceCLZMOvpo2HzzqCOStKptqYQSoGS2s4CdgRb80AXqQPMkQOBRaulzTUYK9o7C3ReHP5eY2TPAPsBPEqD8IC8Phg2D3r3hssvg/vujjkjSqmqpRFULUEslJPPt7e47NbWSZNbxLXX30e4+z90/rXo0NYBkhFsxVd2OrQ1wODA9nTFkql12CZLfAw9orXTOKS4Ouj1vuEHdn5It3jKznk2tJJlZoP2Ak2lin2s99f8CuBvYAlgOvO/u/c2sE/Bvdz/SzLYDngnfUgA85u43NVR3Lk+CqW7VqmDLpPV5qzl32N0c2v1A7RYhInWK8SSYmQQL4ucR5KNmXwbxCEGf64dU63N197OT+cAoKAH+4G8Pz+T3p/fASq6j5aG3asskEalTjBPgNrUdT7ZXMpkxwJT0uUq01nQbhfV6D594JWt3ep7S+aVKgBLQUgmJOTMzD9SZ6KrKJFJfMmOAKelzlWiVdCuh6JjLoM1S/JkHKO54cNQhSRxULZUYMiT4qa2VJJ4mmNnFZta1+kEzKzSzQ8zsQeCMRCtLJgHuB7xvZrPM7AMzm2ZmHyTxfomB4i7FjD//Kc685nV8yS6M+dd+UYckcaBdJSQzDCC4DebjZvZZuAvEJ8DHBHNU7nT3BxKtLJkxwJT0uUZBY4C1O/98GD4cJk6EA5p8TwXJaFosL7WI6xggBDe9BjYHVrv78kbVkcDNsBvsT02mzzUKSoC1W7kSdt89eD51Kmy0UbTxSMQ0Big1xDkBpkIiCbAUeBp41t0XVDteSLAH0xnAhGSanemmBFi311+Hvn1h8GC4996ooxGROMn2BJjIGGBtfa7zaGSfq8TLgQcGt0gbNgxefDHqaERE0iep3SBS0ecaBbUA67dmDfTpA8uWwfTp0K5d1BFJbKmbNKfEtQVoZmMJdg2a2qR6GrkdUpG7r224ZDwoATbsvfdgn33ghBPgsceijkZiSRNlck6ME+CewF+B+cCV7v55Y+pJZhlEdf8ws8Ma+V6JoT32gGuugccfhyefjDoaiSUtlZCYcPf/ufvBwPPAS2Z2bbhPbVIalQDd/Rygq5n93cy0uU6WuOKKoBV4wQXweaP+npKsVrWrRH6+dpWQyJmZAbOAfwIXAx+b2S+TqqORXaD9gQOBPQluRPqgu49KuqI0URdo4j76KGgN7rn/Nxx17b0cvG2JbpUmP9AYYE6JcRfom8C2BPemngRMBj4CfgsUufvghOpJYiH8m8DV7j7BzM4EJrr7J+G5u9z9t0lfRZooASbnd0Pncdd122LHnEfLfR/WDbNFclSME+AuwIza1p+b2Ux375FIPcl0gZ4HXGRm44BZVckvNCSJeiTmtjj4CdjuVfzFO1m7eCdK55dGHZKIyAbu/mE9N185KtF6Ek6A7j7d3Y8n2BX+KjN73sx2D899l2g9En+HbFdCyxPPhVbf4CNHsnf7flGHJJmkrAxuvlk31M5xZnaCmX1oZpVm1qeecvPDe0u/b2ZTmvq5NRpn9ceY7BigmW0M9ADOBs5x92S2VIqEukCTV7awjPufncN9vzuN444zRowAs6ijktjTUoms0pQuUDPrQbB37DCCNXu1Jjczmw/0cfevGh1oIyWcvMxsArAjsBqYET7ObJ6wJGrFXYopvqiYHVbBn/4U3DHm4oujjkpir7alEkqAOcndZwJYjP9yTqb19gdgpruvbq5gJH4uvTS4X+gf/hAskdh336gjklirWipR1QLUUolMV1CjW3K4uw9P8Wc48IqZOTCsGeqvU6OWQWQadYE2zbJlsOee4A7/+x+0bx91RBJrWiqRNRrqAjWzV4GOtZy6yt2fDcuUUn8X6NbuvtjMOgBjgYvdfWLTo2+YEqAk5J13gj0DDz0UnnsO8hp7DyERyRipWAbRUAKsUXYosNLdb2/KZyZKX2OSkL33hjvugDFj4NZbo45GRLKBmbUxs7ZVz4HDgenp+nwlQEnYhRfCSSfB1VfrNpDSRFoqkfXM7BdmtggoBl4ws5fD453MbExYbEvgDTObCrwNvODuL6Utxrh0gZrZbcAxQDkwFzirti2XzGwAcBeQD/zb3W9pqG51gabOihVBa3D5cnj/fehYW++/SH20VCJjxPVOMKkSpxbgWGBXd98NmA1cUbOAmeUD9wBHAD2Bk82sZ1qjzHFt28JTT8F338FRx33LjRNuoWyh/oqXJGhXCYmJ2CRAd3/F3SvCl5OAzrUU2weY4+6fuHs58ARwbLpilMCuu8IfbpzD/8o24ZohBfR7qJ+SoCROu0pITMT1Li5nAyNqOb41sLDa60VArSvTzGwwMBigsLAw1fHlvNZ7P4ntvQn+5qWsafcxpQeV6obZkpji4qDbU0slJGJpTYAJrhm5CqgAHm3KZ4WLKYdDMAbYlLrkp0q6lVB09OGs+WY7/IV7aP3/ZkcdkmSS4mIlPolcWhOgux9a3/lwm6WjgX513Ol7MdCl2uvO4TFJs+IuxYw/6xVe2uctHr3kQK65sCeH9oZddok6MhGRxMRpFugA4A6gr7svraNMAcEEmX4Eie8d4BR3/7C+ujULtHktWBDcIq2oCCZPhi23jDoiEUkFzQJNn/8D2gJjw20x7oUfrxkJJ8lcBLwMzARGNpT8pPl17RrcHWbJEvj5z2HVqqgjkqygtYLSzGLTAmxOagGmx6hRcNxxwWPkSN0uTZpAawVjQS1AkQQNHAi33w5PPw1XXhl1NJLRtFZQ0iCuyyAkQ11yCcyZE9wvdIcd4Nxzo45IMpK2VZI0UBeopFxFBRxzDIwdC3c8NIPvuzxLSbcSrROU5GhbpchlexeoEqA0i+++gz33/Z658yrIO/dAiraaw7jTxykJimSQbE+AGgOUZrHxxnD8DQ9A4fdUPvgya7/chtL5pVGHJSKygRKgNJuB++5J0VlHg+dT+eCr7Gj9ow5Jso2WSkgTqAtUmlXZwjIeH/chD15yBptt0oLXXoNttok6KskKWirR7NQFKtIExV2K+fuZ5zJhXAu+/RYOOQQWLYo6KskKWiohTaQEKGmx557w8suwdGnwR/vnn0cdkWQ8baskTaQuUEmrt96Cww8PukEnTIAOHaKOSDKalko0q2zvAlUClLR77TU44ohgofyECdC+fdQRiUhtsj0BqgtU0q5vXxg9GmbPhsMOg5envc3Nr9+sXeVFJK3UApTIvPgiHDuwkooO72s1dYcAABSeSURBVGK/7E9RmzVaLC+pp27SRlMLUKSZHHEEnDT0afyz3lTe/yprl2+mxfKSWlVLJYYMCX5qvaBUowQokbrgtM4UnjYIlnXH//MG21UOiDokySZaKiH1UAKUSBV3Kab0xsu58O6n2CSvMxcN2oPJk6OOSrKGlkpIPTQGKLExZw707w9ffBFsqHvUUVFHJFlBY4CNlu1jgEqAEitffglHHglTp8Lw4XD22VFHJJK7sj0BqgtUYmXLLYM/1vv1g3POgRtvhBz4G01EIqAd4SV22raF554LdpMfMgQ++wxO+WMZry8s1ca6IpIy6gKV2HKHK66AW2+FvB6j4fhTKGpZqbWCknoaJ6yVukBFImIGt9wCR/9mLJUfHU3lvyeydunWWisoqaW1gjkrNgnQzG4zs4/M7AMze8bMNq2j3Hwzm2Zm75vZlHTHKel35aUbBWsFl29L5bC3aTFnYNQhSTbRWsFmkcR3+gAzm2Vmc8zs8nTGGJsECIwFdnX33YDZwBX1lD3Y3Xu7e5/0hCZRKu5STOlNl3HpA4/TfYcWXHZ2D4YMCb6vRJpMawWbS4Pf6WaWD9wDHAH0BE42s57pCjCWY4Bm9gtgkLufWsu5+UAfd/8q0fo0Bpg9Vq+Giy6C++4LbqT92GOw+eZRRyUZT2OAtUrVGGBd3+lmVgwMdff+4esrANz95qZ+ZiLi1AKs7mzgxTrOOfCKmb1rZoPTGJPEQKtW8J//wL/+BRMnBhvtvv121FFJxisuDmZcKfk1l7q+07cGFlZ7vSg8lhZpXQZhZq8CHWs5dZW7PxuWuQqoAB6to5oD3H2xmXUAxprZR+4+sZbPGgwMBigsLExJ/BIf554Le+wBgwbBgQfCXXfBbkeW8dqnWiohkkIFNeZaDHf34VUvUvSdHplYdYGa2ZnAeUA/d1+VQPmhwEp3v72+cuoCzV7LlsFppwVbK+X3fgQ/6kKKWlVoqYSkXg52kza1C7Sh73R1gYbMbADwR+DndSU/M2tjZm2rngOHA9PTF6XETbt28Pzz0O/siayfegqV/3iPtXP31VIJSS0tlUhaIt/pwDvAjma2rZkVAicBo9MVY2wSIPB/QFuCbs33zexeADPrZGZjwjJbAm+Y2VTgbeAFd38pmnAlLvLy4IahLSg8uz8AlfdP4L37z0SNfkkZLZVojAa/0929ArgIeBmYCYx09w/TFWCsukCbi7pAc0PZwjJemfkmHz5xKk/evxXbbx/MFj3ooKgjk4xX1QIsLw+WSowblxPdoNl+JxglQMlKEyfCWWfBvHlw8cXw5z9Dm6z931jSQmOAWUcJULLW998HM9vvvhu23x7uvx8KupVROl8zRUUSoQSYBZQAc9trrwX7Cs6b5+Ttdw9+8NUUtS7XTFGRBmR7AozTJBiRZtG3L3zwAez7i3dZX3YRlX+fyZopJzH+k9KoQ5NsU1YGN9+sWaIZQvsBSk5o0wbuuHMdB2/Rl7Uv3IqPuo+H56zkoH8EC+lFmixHJ8pkMrUAJWcUdylmwlW3cNOjExh692y+X74RBx0EJ54YTJYRaRItlcg4GgOUnLVqFdx+e7Dh7vr18PvfB5Nmpi/XRBlphCxsAWb7GKASoOS8RYvgyivh4Yeh/RblfHfAb1i/230UtSjQRBlJTpYtlVACzAJKgJKIt9+GE85exIIPO0O72dhBt3L9b7pz9cF/ijo0kUhkewLUGKBIaJ994PExCyk85SQoWomP+g/3nP5b7r0X1q6NOjoRSTW1AEVqKFtYxoR5pRTNG8hT9/Zg0iTo1An++Ef41a+gdeuoIxRJj2xvASoBitTDHcaPhxtuCBbUd+gAg875lC36PkX/nvtrfFCSl0HjhEqAWUAJUFLh9dfh0quX8/bETaHVMgr2epDH/tyXEw7cM+rQJFNk2EzRbE+AGgMUSdCBB8LAG/9J3q+Kodt4Kt66mBMP2pMBA+DZZ6GiIuoIJfa0VjBWdCcYkSSUdCuhaJsbKO98Ei2+34bTKkp5cUQXBg6Ezp3hvPPg3HNh3jqtJZRalJQELb+qFmBJSdQR5TR1gYokqWzhj5NbRQU89xz8858wdizkF1TCzv/F+9xL4XZljD/zVSVB+YHGAGNDCVAkhWbPhnOGTOaN53aE1e1gkwUceORn3HXZfvTuDWZRRyiSOCXALKAEKOlUtrCMQ/5zFGunH4VNPwWb25/1FXnstBOcfHLw6N79py1JkbhRAswCSoCSbtWTW/fWxTz9NDz+eLCUwh122nUlczvfSGXPxylqt1S3XJMfi0k3qRJgFlAClLhYvBhGjoS/Dv+cxR9tFRzc6n8cfNhq/nzBz9h7b8jPjzZGiViMlkpkewLUMgiRNNp6a7jkEnjylfkU/a4XduiV5BWu5rVH9qe4GDp2hNNPhxEjYPnyoCV58+s3U7ZQG6zmDC2VSBstgxCJQHGXYib8fnjYTZrHTm2Ml1+GF14IHg8/DPn5jnepwHf4jsId/sTYP97KgduqmzTraalE2qgLVCRm1q+HyZPhqnvfpHRsG/iiNwCFrdZy8EFF9O0LfftCnz7B96Mm02QhjQGmRWwSoJndABwLVAJLgDPd/bNayp0BXB2+vNHdH2yobiVAyURlC8vo91A/1i7fjPyFB3N04e3MercjM2YE51u3hp57LOf9ln+nsmsphZ2nMf5Xo5UEJWWUANPEzDZ29+/C578Berr7+TXKtAOmAH0AB94F9nL3b+qrWwlQMlVtrbulS2HixGBG6ZNjvuSLuVsGhW09Hbot48i+W7D33rD33rDbblBUpFaiNI4SYATM7Aqgq7tfUOP4yUCJu58Xvh4GlLr74/XVpwQo2apsYRmH3DuI8k/3Iu+zfdnbLmTOtM1YujQ4X1gI2+28ko9bPkZlx3co7DSLMb+7jUN22jfawCUjKAGmkZndBJwOfAsc7O5La5y/FGjp7jeGr4cAq9399lrqGgwMBigsLNxrrXY0lSxVs3XnDgsWwDvvBI+nxs7nkxntYO3GG96z3Xaw667Bo1ev4OeyVpN487MJaiVmijSMEyoBpvLDzF4FOtZy6ip3f7ZauSsIEt21Nd6fcAKsTi1AyWVlC8s45IFDKf+qMwVf9eaMTrfx7YKuTJ8Os2YFk24AyFsH7T8mr/0nnNR3T/ru2Ykdd4Qddww2BM7LU1dqbKRprWC2J8C0LoNw90MTLPooMAa4tsbxxUBJtdedgdImByaSxYq7FDP+zFerJa6uG86tXRvcv/Smp55l5ISZ+NKdqFy2IyPv78Bjw36oo1Ur6LTN98zL+xLfrAUFm43kz8dtSv8+PdhmG9j4h8alkmQ61LZWMOY31o6j2HSBmtmO7v5x+PxioK+7D6pRph3BxJeqHUj/RzAJZll9dasFKFK/qhmn5evLKcwvZOxp4+hixXz8MRseYyZ9zEezK+Cb7WB90Y/ev+mmsM020LbDMspWjKCy7XwKNlnCX4+/lMN670KnTtC27Q83A1eSbCK1AFMiTgnwaWAngmUQnwLnu/tiM+sTPj83LHc2cGX4tpvc/f6G6lYCFGlYQ0lpw7KMdesoXNOZu/cfzcZrevHpp2x4TJmxhC8WF8HaTX7y/jZtYKutYKP23zFt1UtUtvmc/I2Wcelhp7PfTtvToQNssQV06BAky0mLGk6SOZ1INQbYZLFJgM1JCVAkNRJOkt8X0WJVV24vfozNKnbh88/hs8+Cx+SZC5i/qBxWbgnlbWv9nBaFlaxr+Tm0XkJeq285aOde7Ni5Pe3aseGxtHIWQyf/horCr2jRehWjzniIQ3vuTUHBT2NSIm2cpiRAM7sNOAYoB+YCZ7n78lrKzQdWAOuBCnfv0/iIk4xRCVBEUinRJFm+vpwWlRvz2ICxdG2xF0uWBGsclyyB5/5XxuszZ+Gr2sGa9nTI607emi34+mtYt67+z2/dGjbZJHjkt1rBzBWTqCz8jvyiVfxi90PZcaut2Ggj2GijoKX52ZpZXD/pj1QUfEuLonXcP+gefrZDb1q3Dupq1Srouk1VIs2kZNvEBHg4MN7dK8zsVgB3/1Mt5eYDfdz9qyYF2whKgCKSdskkycL8wg3bRbnDqlWwbBlMmDGVX438ExWr2pBf3p4Le11Bu7xt+fZbNjzenz+POZ8vhbUbwbqNaLl+C9atbvXDzNcEtWy1njW2DFqsxlqsYfsttqb9xm1o2TJIkK1awcrKrxi/8AUq81eR16KCk3c/ju222JqiIjY8Fn0/lzveuZmKvO9p0cK56bBr2X3rHhQW8qPH9K/e490vJ/Gzbvuw3zZ70aIFP3oUFCTWRTxt1HC+fvFp2h9xPL0GDk763ylVXaBm9gtgkLufWsu5+SgBNh8lQJHMk4rWVG2JdL/OxZSXw4oVsHIlvDHnfc558ndUrCmiYP3GXFV8E1u37M6qVWx4jJ81mTfnvodXFGEVrdl50z3o0ro7a9bA6tXB4/NvlvH1ilVQUQTri8ivbM36dc040T6vIli6kldB25atKCos2JAcCwpgtzWlPPb5ERSyhvIWMHfksKSToJmVA9OqHRru7sOTDdXMngNGuPsjtZybB3xDcHevYY2pv7GUAEUkqzVXIq1ZrrYy+24dJNu1a4PHW/Pe5aQRZ7BuHbTwNtwz4F/svNlulJez4fHE+//l8Q+exivyscoifr7j8fTbZgDr1rHhMWHuG0yY+wZeaVhlIft1OpDeHfpQUcGGx15v/Y5fz72bAipZZ/Dmrw6nZNjLSf3uGmoBJrK228yuIrh95XFeS8Ixs63DCY8dgLHAxe4+MalAG0kJUEQkAekaA2xssq1ZZtqo4Wx/4nm0qIB1BY1uATapC9TMzgTOA/q5+6oEyg8FVjZ0c5NUUQIUEYmZVCXbKMcAzWwAcAfBmu6ldZRpA+S5+4rw+Vjgend/qTGfmXSMSoAiIlKbJibAOUAR8HV4aJK7n29mnYB/u/uRZrYd8Ex4vgB4zN1vanLgicaoBCgiIrXJ9oXweVEHICIiEgUlQBERyUlKgCIikpOUAEVEJCcpAYqISE7KiVmgZlYJrG7k2wuAihSGE3e63uyVS9cKut5UaOXuWdtQyokE2BRmNiWd23NETdebvXLpWkHXKw3L2swuIiJSHyVAERHJSUqADUvb1hwxoevNXrl0raDrlQZoDFBERHKSWoAiIpKTlABFRCQnKQGGzGyAmc0yszlmdnkt54vMbER4frKZdUt/lKmTwPX+3sxmmNkHZjbOzLaJIs5UaOhaq5U73szczDJ6Knki12tmJ4b/vh+a2WPpjjGVEvhvuauZTTCz98L/no+MIs5UMLP7zGyJmU2v47yZ2d/D38UHZrZnumPMKO6e8w8gH5gLbAcUAlOBnjXKXAjcGz4/CRgRddzNfL0HA63D5xdk6vUmcq1hubbARGAS0CfquJv533ZH4D1gs/B1h6jjbubrHQ5cED7vCcyPOu4mXO9BwJ7A9DrOHwm8CBiwHzA56pjj/FALMLAPMMfdP3H3cuAJ4NgaZY4FHgyfPwX0MzNLY4yp1OD1uvsEd18VvpwEdE5zjKmSyL8twA3ArcCadAbXDBK53l8B97j7NwDuviTNMaZSItfrwMbh802Az9IYX0q5+0RgWT1FjgUe8sAkYFMz2yo90WUeJcDA1sDCaq8XhcdqLePuFcC3QPu0RJd6iVxvdecQ/FWZiRq81rCbqIu7v5DOwJpJIv+23YHuZvammU0yswFpiy71ErneocBpZrYIGANcnJ7QIpHs/9s5rSDqACTezOw0oA/QN+pYmoOZ5QF3AGdGHEo6FRB0g5YQtOwnmlkvd18eaVTN52TgAXf/q5kVAw+b2a7uXhl1YBIttQADi4Eu1V53Do/VWsbMCgi6Ur5OS3Spl8j1YmaHAlcBP3f3tWmKLdUauta2wK5AqZnNJxg3GZ3BE2ES+bddBIx293XuPg+YTZAQM1Ei13sOMBLA3cuAlsDmaYku/RL6f1sCSoCBd4AdzWxbMyskmOQyukaZ0cAZ4fNBwHgPR50zUIPXa2Z7AMMIkl8mjxHVe63u/q27b+7u3dy9G8F458/dfUo04TZZIv8tjyJo/WFmmxN0iX6SziBTKJHrXQD0AzCzHgQJcGlao0yf0cDp4WzQ/YBv3f3zqIOKK3WBEozpmdlFwMsEs8ruc/cPzex6YIq7jwb+Q9B1ModgEPqk6CJumgSv9zZgI+DJcK7PAnf/eWRBN1KC15o1Erzel4HDzWwGsB64zN0zsjcjwev9A/AvM7uEYELMmZn6x6uZPU7wx8vm4ZjmtUALAHe/l2CM80hgDrAKOCuaSDODboUmIiI5SV2gIiKSk5QARUQkJykBiohITlICFBGRnKQEKCIiOUkJUCRJZrapmV1Y7XUnM3uqmT5roJldU8/5Xmb2QHN8tki20zIIkSSFW2E97+67puGz3iJYmP9VPWVeBc529wXNHY9INlELUCR5twDbm9n7ZnabmXWr2p/NzM40s1FmNtbM5pvZReHeiu+FN55uF5bb3sxeMrN3zex1M9u55oeYWXdgbVXyM7MTzGy6mU01s4nVij5HBt+YQSQqSoAiybscmOvuvd39slrO7wocB+wN3ASscvc9gDLg9LDMcOBid98LuBT4Ry31/Az4X7XX1wD93X13oPpdeaYABzbhekRykm6FJpJ6E9x9BbDCzL4laKEBTAN2M7ONgP354TZzAEW11LMVP75n5ZvAA2Y2EvhvteNLgE4pjF8kJygBiqRe9Z0zKqu9riT4fy4PWO7uvRuoZzXBriMAuPv5ZrYvcBTwrpntFd7Ds2VYVkSSoC5QkeStINhGqVHc/TtgnpmdABDeuX/3WorOBHaoemFm27v7ZHe/hqBlWLXtTXdgemPjEclVSoAiSQpbXW+GE1Jua2Q1pwLnmNlU4EPg2FrKTAT2sB/6SW8zs2nhhJu3gKnh8YOBbNjNXiSttAxCJMbM7C7gOXd/tY7zRcBrwAHuXpHW4EQynFqAIvH2Z6B1Pee7Apcr+YkkTy1AERHJSWoBiohITlICFBGRnKQEKCIiOUkJUEREcpISoIiI5KT/D9pRcknNGBqQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#\n",
    "# No data this time. We'll just compare the analytical and computer\n",
    "# solutions to the problem.\n",
    "#\n",
    "\n",
    "m = 0.0038  # kg (mass of cup plus 6 paper clips. See proj. 1)\n",
    "g = 9.8     # m/s^2, gravitational field strength at Earth's surface\n",
    "b = 0.00375 # effective drag coefficient, measured experimentally\n",
    "\n",
    "vi = 0.0    # intial velocity\n",
    "yi = 0.0    # initial height\n",
    "\n",
    "s = np.array([yi,vi])  # <-- now 'state' is position *and* velocity\n",
    "\n",
    "dt = 0.03   # pick a reasonable time step\n",
    "t = 0.0     # start at t=0.0\n",
    "tf = 1.0    # stop in 10 sec\n",
    "\n",
    "def derivs_drag(s,t):\n",
    "    \"\"\"\n",
    "    The derivs_drag function computes the time derivative of the \"state\"\n",
    "    for this particular problem. Now the \"state\" is both the position\n",
    "    *and* the velocity of the falling muffin tin.\n",
    "    \"\"\"\n",
    "    y=s[0]  # <- y is the first element of the state vector\n",
    "    v=s[1]  # <- v is the second element of the state vector\n",
    "\n",
    "    dvdt = b*v**2/m - g   # in this problem dv/dt = a doesn't depend on y\n",
    "    dydt = v\n",
    "    \n",
    "    return np.array([dydt, dvdt])\n",
    "\n",
    "def HeunStep(s, t, derivs, dt):\n",
    "    \"\"\"\n",
    "    HeunStep is a generic implimentation of the Heun Method.\n",
    "    Notice that it knows *nothing* about the state, nor the physics\n",
    "    of the problem (which is handled by the 'derivs' function\n",
    "    that is passed in by the programmer). Note that this function has *not*\n",
    "    changed at all.\n",
    "    \"\"\"\n",
    "    f1=derivs(s,t)\n",
    "    f2=derivs(s+f1*dt,t+dt)\n",
    "    return s + 0.5*(f1+f2)*dt\n",
    "\n",
    "v0=np.sqrt(m*g/b)  # terminal velocity\n",
    "tau=v0/(2*g)    # time constant\n",
    "\n",
    "def VAnalytical(t):\n",
    "    return v0*(1-np.exp(t/tau))/(1+np.exp(t/tau))\n",
    "\n",
    "tlist = [t]\n",
    "ylist = [yi]\n",
    "vlist = [vi]\n",
    "\n",
    "while t <= tf:\n",
    "    s = HeunStep(s, t, derivs_drag, dt)\n",
    "    t += dt\n",
    "    tlist.append(t)\n",
    "    ylist.append(s[0])  # append the y value\n",
    "    vlist.append(s[1])  # append the v value\n",
    "\n",
    "fig, ax1 = plt.subplots()\n",
    "ax1.set_xlabel(\"time (s)\")\n",
    "ax1.set_ylabel(\"$v_y$(m/s)\")\n",
    "ax1.set_title(\"Falling Muffin Cup (Heun)\")\n",
    "v_heun = np.array(vlist)\n",
    "v_analytical = VAnalytical(np.array(tlist))\n",
    "ax1.plot(tlist,v_heun,'g.',tlist,v_analytical,'b-')\n",
    "\n",
    "ax2 = ax1.twinx()\n",
    "ax2.set_ylabel(\"$y$ (m)\")\n",
    "ax2.plot(tlist,ylist,'r.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Project 2 (option 1): Spaceshot to Proxima Centauri\n",
    "============================\n",
    "\n",
    "Now we're going to combine all this into our next project. Proxima Centauri is about 4.24 light-years away. It is the nearest start to our Sun. It is dreamed that some day humans will develop the technology needed to send astronauts to visit Proxima Centauri. Let's assume we have that technology already and consider modeling how that might work. Suppose you have an advanced propulsion system capable of producing incredible quantities of thrust with very little fuel. Let's design the mission so that humans can arrive at Proxima Centauri in five years or less (according to their own biological clocks).\n",
    "\n",
    "Use the **Heun** **Method** to compute the position and velocity of the space-craft as it leaves our solar system and heads to Proxima Centauri. You should use a `state` array like the one used in the last Heun Method example. Ignore gravitational effects from our Sun and/or Proxima Centauri. Just assume that you have control of the thrust which can be directed toward or away from the Sun and turned on and off. You can make the thrust whatever you like, but understand that the human body will not survive an acceleration of much more than 4 of 5 times $g$.\n",
    "\n",
    "Think about how you can do this and what limitations you might encounter. We'll discus this more in class. This project is a good option if you've had Modern Physics (PHYS 250) *or* if you want to learn about the effects of special relativity in a concrete situation.\n",
    "\n",
    "Project 2 (option 2): Baseball with quadratic drag\n",
    "============================\n",
    "\n",
    "Suppose you're modeling a baseball in flight, but including the effect of air drag. Assume a baseball is launched with an initial speed of 80 m/s in a direction 35 degrees above the horizontal.  You may assume that the magnitude of the drag force is given by:\n",
    "\n",
    "$$ f_d = \\frac{1}{2} C_d A \\rho v^2 $$\n",
    "\n",
    "Where $C_d$ is the drag coefficient (around 1), $A$ is the cross sectional area of the ball, $\\rho$ is the density of the air and $v$ is the speed of the ball (magnitude of velocity).\n",
    "\n",
    "Use Heun's method to find the trajectory of the ball and, in particular, the maximum altitude and range (distance from the batter where the ball reaches the ground). Remember, position and velocity will both be two component vectors so the `state` vector will have four components over all. This project is a good option if you don't want to mess with relativity.\n",
    "\n",
    "Project 2 (option 3): You come up with your own problem\n",
    "==============================\n",
    "\n",
    "It's OK with me, if you don't like either of these two options, to come up with your own. However it must be at least a second order system (have a state vector with two components) and you must use the Heun method to solve it.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}