CompMech03-IVPs_project_commit1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Initial Value Problems - Project\n",
    "\n",
    "![Initial condition of firework with FBD and sum of momentum](firework.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are going to end this module with a __bang__ by looking at the flight path of a firework. Shown above is the initial condition of a firework, the _Freedom Flyer_ in (a), its final height where it detonates in (b), the applied forces in the __Free Body Diagram (FBD)__ in (c), and the __momentum__ of the firework $m\\mathbf{v}$ and the propellent $dm \\mathbf{u}$ in (d). \n",
    "\n",
    "The resulting equation of motion is that the acceleration is proportional to the speed of the propellent and the mass rate change $\\frac{dm}{dt}$ as such\n",
    "\n",
    "$$\\begin{equation}\n",
    "m\\frac{dv}{dt} = u\\frac{dm}{dt} -mg - cv^2.~~~~~~~~(1)\n",
    "\\end{equation}$$\n",
    "\n",
    "If we assume that the acceleration and the propellent momentum are much greater than the forces of gravity and drag, then the equation is simplified to the conservation of momentum. A further simplification is that the speed of the propellant is constant, $u=constant$, then the equation can be integrated to obtain an analytical rocket equation solution of [Tsiolkovsky](https://www.math24.net/rocket-motion/) [1,2], \n",
    "\n",
    "$$\\begin{equation}\n",
    "m\\frac{dv}{dt} = u\\frac{dm}{dt}~~~~~(2.a)\n",
    "\\end{equation}$$\n",
    "\n",
    "$$\\begin{equation}\n",
    "\\frac{m_{f}}{m_{0}}=e^{-\\Delta v / u},~~~~~(2.b) \n",
    "\\end{equation}$$\n",
    "\n",
    "where $m_f$ and $m_0$ are the mass at beginning and end of flight, $u$ is the speed of the propellent, and $\\Delta v=v_{final}-v_{initial}$ is the change in speed of the rocket from beginning to end of flight. Equation 2.b only relates the final velocity to the change in mass and propellent speed. When you integrate Eqn 2.a, you will have to compare the velocity as a function of mass loss. \n",
    "\n",
    "Your first objective is to integrate a numerical model that converges to equation (2.b), the Tsiolkovsky equation. Next, you will add drag and gravity and compare the results _between equations (1) and (2)_. Finally, you will vary the mass change rate to achieve the desired detonation height. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Create a `simplerocket` function that returns the velocity, $v$, the acceleration, $a$, and the mass rate change $\\frac{dm}{dt}$, as a function of the $state = [position,~velocity,~mass] = [y,~v,~m]$ using eqn (2.a). Where the mass rate change $\\frac{dm}{dt}$ and the propellent speed $u$ are constants. The average velocity of gun powder propellent used in firework rockets is $u=250$ m/s [3,4]. \n",
    "\n",
    "$\\frac{d~state}{dt} = f(state)$\n",
    "\n",
    "$\\left[\\begin{array}{c} v\\\\a\\\\ \\frac{dm}{dt} \\end{array}\\right] = \\left[\\begin{array}{c} v\\\\ \\frac{u}{m}\\frac{dm}{dt} \\\\ \\frac{dm}{dt} \\end{array}\\right]$\n",
    "\n",
    "Use [two integration methods](../notebooks/03_Get_Oscillations.ipynb) to integrate the `simplerocket` function, one explicit method and one implicit method. Demonstrate that the solutions converge to equation (2.b) the Tsiolkovsky equation. Use an initial state of y=0 m, v=0 m/s, and m=0.25 kg. \n",
    "\n",
    "Integrate the function until mass, $m_{f}=0.05~kg$, using a mass rate change of $\\frac{dm}{dt}=0.05$ kg/s. \n",
    "\n",
    "_Hint: your integrated solution will have a current mass that you can use to create $\\frac{m_{f}}{m_{0}}$ by dividing state[2]/(initial mass), then your plot of velocity(t) vs mass(t)/mass(0) should match Tsiolkovsky's_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "def heun_step(state,rhs,dt,etol=0.000001,maxiters = 100):\n",
    "    '''Update a state to the next time increment using the implicit Heun's method.\n",
    "    \n",
    "    Arguments\n",
    "    ---------\n",
    "    state : array of dependent variables\n",
    "    rhs   : function that computes the RHS of the DiffEq\n",
    "    dt    : float, time increment\n",
    "    etol  : tolerance in error for each time step corrector\n",
    "    maxiters: maximum number of iterations each time step can take\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    next_state : array, updated after one time increment'''\n",
    "    e=1\n",
    "    eps=np.finfo('float64').eps\n",
    "    next_state = state + rhs(state)*dt\n",
    "    ################### New iterative correction #########################\n",
    "    for n in range(0,maxiters):\n",
    "        next_state_old = next_state\n",
    "        next_state = state + (rhs(state)+rhs(next_state))/2*dt\n",
    "        e=np.sum(np.abs(next_state-next_state_old)/np.abs(next_state+eps))\n",
    "        if e<etol:\n",
    "            break\n",
    "    ############### end of iterative correction #########################\n",
    "    return next_state\n",
    "\n",
    "def rk2_step(state, rhs, dt):\n",
    "    '''Update a state to the next time increment using modified Euler's method.\n",
    "    \n",
    "    Arguments\n",
    "    ---------\n",
    "    state : array of dependent variables\n",
    "    rhs   : function that computes the RHS of the DiffEq\n",
    "    dt    : float, time increment\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    next_state : array, updated after one time increment'''\n",
    "    \n",
    "    mid_state = state + rhs(state) * dt*0.5    \n",
    "    next_state = state + rhs(mid_state)*dt\n",
    " \n",
    "    return next_state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.rcParams.update({'font.size': 22})\n",
    "plt.rcParams['lines.linewidth'] = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simplerocket(state,dmdt=0.05, u=250):\n",
    "    '''Computes the right-hand side of the differential equation\n",
    "    for the acceleration of a rocket, without drag or gravity, in SI units.\n",
    "    \n",
    "    Arguments\n",
    "    ----------    \n",
    "    state : array of three dependent variables [y v m]^T\n",
    "    dmdt : mass rate change of rocket in kilograms/s default set to 0.05 kg/s\n",
    "    u    : speed of propellent expelled (default is 250 m/s)\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    derivs: array of three derivatives [v (u/m*dmdt) -dmdt]^T\n",
    "    '''    \n",
    "    dstate = np.zeros(np.shape(state))\n",
    "    dstate[0] = state[1]\n",
    "    \n",
    "    dstate[1] = (u/state[2])*dmdt\n",
    "    dstate[2] = -dmdt\n",
    "    \n",
    "    \n",
    "    return dstate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The number of time steps is 10000.\n",
      "The time increment is 0.1\n",
      "40\n"
     ]
    }
   ],
   "source": [
    "dt = 0.1 \n",
    "T = 1000\n",
    "N = round(T/dt)\n",
    "\n",
    "print('The number of time steps is {}.'.format( N ))\n",
    "print('The time increment is {}'.format( dt ))\n",
    "\n",
    "# time array\n",
    "t = np.linspace(0, T, N)\n",
    "\n",
    "x0 = 0    # initial position\n",
    "v0 = 0    # initial velocity\n",
    "m0 = 0.25    # initial mass\n",
    "\n",
    "\n",
    "#initialize solution array\n",
    "num_heun = np.zeros([N,3])\n",
    "num_rk2 = np.zeros([N,3])\n",
    "\n",
    "#Set intial conditions\n",
    "num_heun[0,0] = x0\n",
    "num_heun[0,1] = v0\n",
    "num_heun[0,2] = m0\n",
    "num_rk2[0,0] = x0\n",
    "num_rk2[0,1] = v0\n",
    "num_rk2[0,2] = m0\n",
    "\n",
    "for i in range(N-1):\n",
    "    if num_heun[i][2] < .05:\n",
    "        end = i\n",
    "        print(i)\n",
    "        break\n",
    "    num_heun[i+1] = heun_step(num_heun[i], simplerocket, dt)\n",
    "    num_rk2[i+1] = rk2_step(num_rk2[i], simplerocket, dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAAEaCAYAAACimQj6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3gU1dvG8e/Z9JBQQu+hSRBEqoCABAJGEKR3kAiIAooNC2KhiAVEBaTITyBEeJEi0kEQ6YIUBaV3BKTXBNL3vH/sJqRns9kwKc/nuvZad8qZZ0FyZ2bOnKO01gghhBA5gcnoAoQQQghbSWgJIYTIMSS0hBBC5BgSWkIIIXIMCS0hhBA5hrPRBeRmRYoU0b6+vkaXIYQQOcq+ffuua62LprROQisL+fr6snfvXqPLEEKIHEUpdS61dXJ5UAghRI6RbUJLKVVVKfWaUmqeUuqoUsqslNJKqS6ZbLeXUmqbUuqOUipMKbVXKTVUKZXmd7d3PyGEEFknO10eHAy85sgGlVJTgSFABLARiAYCgG+BAKVUV611rKP2E0IIkbWy01nDQWAC0B2oDGzJTGNKqc5YgucyUFNr3VZr3RGoAhwBOgKvOGo/IYQQWS/bhJbW+nut9Tta60Va61MOaHKE9f1drfWJBMe5guWsDuC9FC732bufEEKILJYrf/AqpcoAdYEoYHHS9VrrLcBFoATQMLP7CSGEeDhyZWgBta3vh7TW4alssyfJtpnZz6HO/Xecl2c25dx/x7PqEEIIkSPl1tCqYH1Pta8/8G+SbTOzn8Ns2LWQAWs6ssPtNu+t6k5E5P2sOIwQQuRIuTW0vKzv99LYJsz67u2A/eIppQZZu8fvvXbtWrqFJnXu6mGuuFj+Wg66xTDyh44ZbkMIIXKr7NTl3ZGU9T2jM1zau188rfVMYCZAvXr1MtzOwOdGc2rOAVaZLH1R1rv8x5Qlb/Jql6/sLUlYxcTEcPPmTe7cuUNMTIzR5QiR6zk7O1OgQAF8fHxwdnZM3OTW0Aq1vnulsU3cutAEy+zdz6HG9l3MxVlP8pd7BABzQ9fjt2shrRp2z6pD5npms5nz58/j5uZGuXLlcHV1RSmV/o5CCLtorYmKiuLGjRucP3+e8uXLYzJl/uJebr08eNb6Xj6Nbcom2TYz+zmUs7ML4zsvpUy05UQt0qT44uAYzlw8mlWHzPVu3bqFs7MzJUuWxM3NTQJLiCymlMLNzY2SJUvi7OzMrVu3HNJubg2tv6zv1ZVSHqlsUz/JtpnZz+FKFCnLyHrj8TSbAbjiYmLE6p7cj0jrdptITVhYGAULFpSwEuIhU0pRsGBB7t1zzM+uXBlaWuvzwJ+AK9A16XqlVDOgDJZRL3Zmdr+s0qRWG14q0i3+8yG3GEb80D6rD5srRURE4OnpaXQZQuRJnp6ehIen9hRRxuTo0FJKfWYdXPezFFbHLftCKVU5wT7FgGnWj59rrc0O2i9L9G/3Mc+Zq8R//s31ChN/HJzGHiIlZrPZIdfThRAZZzKZMJsd8yMz2/wrVkrVUUrtinsBdayrPk2yPKGSQFXreyJa6yXAdCyjV/yjlFqplFoKnAAeBZZhGQDXIftlpbHPL6Z+5IOzhHkR21i26buHWUKuIJcGhTCGI//tZZvQAvIDDRK84p6DqpJkuc201kOA3lgu+TUDAoGTWAa87ZzaSO327pdVTE5OTOi2HN8oy+cYpfjqzGT+PpE0w4UQIndTWtv9SJJIR7169bQjZy7ef2w7r+x4iTtOlt81KkbB7K4bKFywhMOOkVsdOXKEatWqGV2GEHlWRv4NKqX2aa3rpbQuO51piXTUqtqE18sPwdn6i8ZpV3h7UXtiYqINrkwIIR4OCa0cpkvAUHq6Nor/vMftPh+GdDawIpEb+Pr6opTi7NmzRpeSqs2bN6OUwt/fP9k6pZRD7psEBQWhlCI4ODjTbSXk7+9vU7txxx81apRDj5+bSGjlQO/0+h8BUcXjP69yOsOUn94ysCIhcrfg4GCUUgQFBRldSp4noZVDff78CmpGuMZ/Dr77Cyu3zjawIpGTbdy4kSNHjlC6dGmjS7HLkSNHOHLkSKbb+eyzzzhy5AgdO8pA1dmVhFYO5e7myYROD4Z6ijIpJpyYyN8nsvyZZ5ELVapUCT8/P1xcXIwuxS5+fn74+fllup2SJUvi5+dHgQIFHFCVyAoSWjlYqaLlGdPwG/LHWh7au+Vs4v3Ng7h8/bzBlYmcJrV7WnH3YjZv3syOHTt45plnKFSoEAUKFCAwMJD9+/fHbxsSEkL9+vXx8vLCx8eHPn36cPny5WTHSnip7fr16wwePJgyZcrg7u5OpUqV+OCDD7h/P2PzyKV1Tys6OpqZM2fSvHlzfHx84gdNbtu2LfPnz0+0bUr3tHx9fXnhhRcAmDt3bvyxjLhceOTIEQYMGECFChVwd3enUKFCtGzZkhUrVqS4fXr3+lL7e0+4fMOGDQQEBFCgQAE8PT1p2LBhqsd7GCS0crj6NVryRrmX43sUnnOF4Us7yOSRwqFWrlxJs2bNuHXrFoGBgZQqVYr169fTrFkzTpw4wfDhwxk4cCAFCxYkMDAQV1dX5s+fT8uWLYmKikqxzVu3btGgQQMWLVpEgwYNCAwM5Nq1a4wbN46AgIAMB1dqx2jWrBkvvfQSO3fupHbt2nTq1IkKFSqwY8cORo4cmW4bXbp0oXHjxoDljLRfv37xryZNmmS6Rlv9+OOP1KpVi9mzZ5MvXz7atm1LzZo12bZtG+3bt+ejjz5y+DFnzZpFYGAgYWFhtGnTBj8/P/744w86dOjAkiVLHH48W+TWqUnylC4tX+XfhceYE7EFgANuUbwb8hyTXvzV4MpyBt/3Vhtdgt3Ofv7sQznO119/zaJFi+jSpQtgGRarT58+LFiwgE6dOnHt2jX279/Po48+CsDNmzdp1KgRhw4dYuHChfTt2zdZmytWrKBx48bs27ePggULAnDlyhVatWrFrl27GDVqFOPHj89U3UFBQezcuZNGjRqxZMkSSpUqFb8uIiKCTZs2pdvGl19+SXBwMDt27KBJkyYO71loi7///pt+/frh6urKsmXLaN26dfy6Q4cO0bp1a8aOHUvz5s1p3ry5w447fvx41qxZwzPPPBO/7JNPPuHDDz9kxIgR8f8/PExyppVLvNn9W1rHPphR5TfXK4ybF2RcQSJX6dGjR6IfUCaTiXfeeQeAgwcPMmbMmPjAAvDx8eHll18GSDUYlFJMnz49PrAAihcvzqRJkwCYMWMGERERdte8f/9+VqxYgZeXF8uXL08UWADu7u6Jfvg/DC+88EKiy4tJX3Pnzk1xv3HjxhEVFcX48eOT1Vy9enW++soySey33zp2hLlXX301UWABvPPOOxQoUICTJ0/y77//OvR4tpAzrVzk0+d/5vqspuxxt0wBsChmL6VXjSWo7YcGVyZyuqQ/uAAqV66c5voqVSwDPf/3338ptlmzZk0ee+yxZMubN29O6dKluXjxIvv27Yu/NJdR69atA6B9+/YULVrUrjYcrXHjxon+3JLavn07p06dSrTMbDazbt06lFKpntk0a9YMgJ07HdsRq23btsmWubq6UrFiRf766y/+++8/ypUr59BjpkdCKxdxdnbhq16r6L+gBSfcNGalmHbtR0r8XpFnnuxtdHnZ1sO6xJaTlSlTJtkyLy8vm9andrZUoUKFVI/n6+vLxYsXuXDhQkZLjXfu3DkAh/QqdJSBAwem2XkjKCgoWWjduHGDu3fvAlCsWLE027927Vqma0wotUDKnz8/kPrfbVaS0MplCnoX4fPAeQzZ0JMrLibCTSY+O/IpJQqXp1bVh3fTWOQu6U3rklXTvsjI/BAbaxmf28nJiT59+ji07fSmC8mO0/lIaOVCj5Svyai6X/DOX+8S6mTiprOJ97cOZmahlZQp5mt0eUIApDlkVNy6pPehMqJ8ecs93mPHjtndRnZQpEgRPDw8CA8P59tvv010hpseFxcXoqOjCQsLS7ZfdHQ0ly5dcnS5WS77xahwiCa12/JGucG4WLvCn3eFt5Z15N79UIMrE8LiwIEDHDx4MNnyLVu2cPHiRby8vKhbt67d7QcGBgKwfPlyrl+/bnc7YLmPAxATE5Opduzh7OxMy5YtATLczTxuhJOjR48mW7d+/XpDvk9mSWjlYl1bvkJ/r6fjPx92i+GNea0xxz7U6cCESJHWmiFDhnDnzp34ZdeuXeO1114DYNCgQXh4eNjdfu3atWnXrh2hoaF07Ngx2VlFREQEa9eutamtuB/+jhgqyh4fffQRLi4uvPbaa/z4448knVLKbDazcePG+M4ncQICAgAYM2ZMouflDh06xKuvvpr1hWcBCa1c7pUuX9FBV43/vNPtDu/NbW9gRUJYPPfcc1y8eJFKlSrRpUsXOnToQOXKlTlw4AD169dnzJgxmT5GcHAw9evXZ/v27VSsWJFWrVrRq1cv/P39KVmyJIMHD7apnYYNG1KiRAn+/PNP6tWrR79+/Rg4cCBz5szJdI22qFevHiEhIURHR9OzZ08qVKhAmzZt6N69O40bN6Z48eK0bNmSzZs3J9pvxIgR5M+fn5UrV1K1alW6dOlCkyZNqFOnDo0bN46/hJqTSGjlAaP7LuSpKJ/4z2udzvH5/P4GViQEFCpUiF27dtGxY0d27tzJ2rVrKVy4MO+//z6bNm0iX758mT6Gj48P27ZtY8qUKdSpU4fdu3ezdOlSzpw5Q9OmTfn8889tasfNzY1169bx7LPPcubMGebNm8esWbPYsmVLpmu0VY8ePfjnn38YNmwYnp6ebNmyhVWrVnH58mXq1KnDpEmTGDZsWKJ9KlWqxI4dO3juuee4ffs2q1ev5s6dO0yYMIGQkJCHVrsjyczFWcjRMxdnRkTkfV4Mfor97pEAKK15pVBHBrUfa3BlD4fMXJx9BAcH88ILL9CvXz9DRpcQxpCZi0WGuLt58k33VVSJtHQh1kox8+bPLNv0ncGVCSGE7SS08pDCBUvwZesFlLJOZxJpUow/M5ntf60yuDIhhLCNhFYeU7FsdcY1mkwh63QmoU4mPvrzXf4+scvgyoQQIn0SWnlQveoteL/qe3hYn4a/5mzivc0DOfffcYMrE3lBUFAQWmu5nyXsIqGVRz3TuC+vl3ohfh6u866KN1d15cbt5JP2CSFEdiGhlYf1ChzOi/lbo6zBddzNzGsLn+V+xD2DKxNCiJRJaOVxQzpNoKfLg56lB9yjeD3kaWJiog2sSgghUiahJRjRO5i2sb7xn3e63eWtOTLckxAi+5HQEgCM67eMZpGF4z//5nqFkSGdDaxICCGSk9ASAJicnPgq6BcaRDwYOmeV6RSf/PC8gVUJIURiEloinqurG5P6bqBmhGv8soXmv/hq4SsGViWEEA9IaIlE8nl68233tVSNfPC/xtzwzcxc/oGBVQkhhIWElkimUMFifNNhORWs0++YlWL6rWWErLVtRGwhhMgqEloiRWWK+TIxcH78OIUxSjH58jwWb/zW4MqEEHmZhJZIVZVyNZnw1EyKxViGe4o0Kb78dzort842uDKR1wQHB6OUIigoKNHys2fPopTC19c308fw9/dHKZVsIsXM8vX1RSmV6OXu7k65cuXo1q1bmnNyBQUFoZRi1KhRKa6/c+cOTZo0QSlFtWrVuHDhAgBXr15l7ty59OjRgxo1auDt7U2+fPmoUaMGb7/9Npcv59yRbyS0RJpqPvIknzeYTGFrcN03mfjs5EQ27FpocGVCPByjRo1KMzhsFRgYSL9+/ejXrx+BgYEALF68GH9/f77++usMt3ft2jWaN2/Ojh07qFOnDtu2baNMmTIAvPnmmwQFBbF48WKcnJxo3bo1/v7+XL16lS+//JLq1auzb9++TH0fozgbXYDI/urXCGBs1GeM2D+CO04mQp1MjD48Bne3fDSt3dbo8kQeVrp0aY4cOYKLi0um2woJCeH+/fuUK1fOAZUl99577+Hv7x//OTo6mtdee43p06fz3nvv0bVr1/jQSc/58+dp1aoVx44do2nTpqxatYr8+fPHr/fx8WH06NEMGDCA0qVLxy8PCwvjxRdf5Mcff6Rbt24cO3YMZ+ecFQNypiVs0rTOc3z86Ad4W6c0ueNk4oM/32Xn378YXJnIy1xcXPDz86NSpUqZbqtcuXL4+fnh6enpgMrS5+LiwsSJE/H29iYqKor169fbtN+JEydo2rQpx44do3Xr1vzyyy+JAgtg8uTJfPTRR4kCC8DLy4tZs2bh7e3N6dOn2blzp8O+z8MioSVs1qphT0ZUfgNP65QmN51NvL/nTfYc2mhwZcIR7t27x/jx46lfvz758+fHw8OD6tWrM2rUKMLCwhJte/z4cfLnz4+zszNbt25N1tbhw4fJly8fLi4uiX4wJrzUdubMGfr06UPx4sVxd3enevXqTJw4kZiYGJtrTu+e1r179/jyyy9p1KgRBQsWxMPDg4oVK9K1a1fWrFmTaNuU7mkppRg9ejQAo0ePTnRfKrOXCwE8PDx45JFHALhy5Uq62//99980bdqUc+fO0a1bN5YvX46Hh0eGjunp6UnVqlUB4u+B5SQ567xQGK7dUwOJjI7gi3+nE2Eycd3ZxHu7hjHBaQZ1/JoaXZ6w04ULFwgMDOTw4cMULVqURo0a4e7uzp49exg9ejQ///wzmzdvplChQgA88sgjfPfdd/Tq1YtevXqxf/9+ihQpAsD9+/fp2rUr9+/fZ/z48TRq1CjZ8c6cOUO9evVwd3fH39+fu3fvsmnTJoYPH8727dv56aefMJky9zv1uXPnCAwM5NixY3h5edGkSRMKFCjA+fPnWbt2LdeuXaNNmzZpttGvXz/279/PgQMHePzxx6lVq1b8uoT/nRl37twBoHjx4mlut2vXLtq0acOtW7d48cUXmTFjhl1/RjExMZw9exaAkiVLZnh/o0loiQzrEvAKseujmHBxNpEmxVVnE+9uH8xEp++pWaWh0eVl3KgCRldgv1F3Mt2E1ppu3bpx+PBhXnnlFb744ov4S2Th4eEMGjSIefPm8cYbbySauLFnz55s3ryZmTNn8vzzz7N69WqUUgwdOpTDhw/Tpk0bhg8fnuIxQ0JC6Ny5M/PmzcPd3R2wXPZq3rw5y5YtY8aMGQwZMsTu72Q2m+nYsSPHjh2jffv2zJkzJz5wAUJDQ9m9e3e67QQHBzNq1CgOHDhAhw4dHHJ2ldChQ4c4c+YMLi4uPP3006lut337dr788kvu3bvH8OHDmTBhgt3H/P7777l+/TolSpTgySeftLsdo8jlQWGX7k+/yRsl+uJqtjzHddlFMXzLQA6d2mNwZSKj1q1bx86dO2nYsCGTJk1KdE/Hw8ODGTNmUKxYMebPn8+tW7cS7Ttp0iRq1qzJ2rVrmTBhAiEhIQQHB1OmTBlCQkJQSqV4TA8PD6ZNmxYfWABVqlRh7NixAHb1pktoxYoV/PXXX/j6+rJgwYJEgQXg7e1NQEBApo6RGbdu3WLt2rV06tQJs9nMpEmT0uyEsXHjRu7du8cTTzyRqcD6559/ePvttwEYP348rq6u6eyR/UhoCbv1bv0uw4r1xMU6ieQlF8Wbm16Q4Mph4u7tdO7cOcXLTfny5aNevXrExMSwZ0/iv1t3d3cWL16Ml5cXI0eOZPDgwTg5ObFgwQIKFy6crK04Tz/9NMWKFUu2vHfv3phMJk6ePMnFixft/k7r1q2Lby+j93yySvPmzePvh/n4+NCmTRvOnTvH2rVrGTx4cJr7NmrUCBcXF3bv3s2wYcPsOv6FCxdo164dYWFhDBw4kL59+9rVjtHk8qDIlH7PjiR2ZTRTbiwhRin+swbXV8yheqX6RpdnGwdcYsvJTp8+DcDbb78d/1t4aq5du5Zs2SOPPML48eMZMmQIMTExjB49miZNmqTZToUKFVJc7urqSsmSJbl48SIXLlxI1vvNVufOnQPAz8/Prv2zQmBgICVKlEBrzeXLl9m6dSsRERE8//zz7Nixg8qVK6e679NPP83w4cPp0aMHU6ZMQSnFpEmTbD725cuXCQgIiO/AMWPGDEd8JUNIaIlM699uFOYVZqbeXJoouL40z+KxKg2MLk+kI9Y62WezZs3SHVmifPnyKe7/448/xn/es2cPWutULw3aKrP7ZzdJn9O6dOkSgYGB/PPPP/Tu3Ztdu3al+Z07derEjz/+SI8ePZg8eTKATcF19epVWrRowfHjx2nfvj3z58/Hyckp09/HKBJawiEGPjcGVpAouN7aMoCJSHBld2XLlgWga9euDB06NMP7jxo1iq1bt/LEE09w7949Vq1axVdffcVbb72V6j5xvdeSioqK4tKlSwCUKlUqw7XEiQvXY8eO2d1GVitZsiSLFi2iZs2a7N69m/nz59OnT58090kaXEopvvnmm1S3v3btGi1atODIkSM8++yzLFq0KMc9TJyU3NMSDjPwuTEM9emMc4J7XG9tGcDfJ3YZXJlIS+vWrQHLkEIZtXHjRj799FMKFizIwoULWbhwIZ6enowYMSLN3nnr169P8VLjggULMJvNVKpUyebRIVISN0zSvHnziIiIsLsdIL6zQkaeH7OVn59ffC/JUaNG2XSMuOBycXFh0qRJvP766ylud/36dVq0aMGhQ4cIDAzkp59+ypEdL5KS0BIONfC50bxauEuS4BrIn0e3GVyZSE2HDh2oW7cuW7Zs4eWXX+bmzZvJtjl9+jRTp05NtOzKlSv07t0bs9nMrFmz8PX1pXr16kyePJno6Gi6d+/O7du3Uzzm/fv3eeWVV4iMjIxfdurUKT788EMAXnvttUx9p/bt21OrVi3Onj1L796945+FihMaGsrGjbY9FB93X+3IkSOZqik1I0eOxNvbm1OnTvHDDz/YtE/S4HrjjTcSrb958yYBAQEcPHiQVq1asWzZMtzc3LKi/IdPa52tXkAvYBtwBwgD9gJDAVMG2vAHtI2vckn2DU5n+6O21lG3bl2dV81aMUrXmlNd1wiuoWsE19Atvn9U//H3r4bVc/jwYcOOnROcP39eP/bYYxrQ3t7eukmTJrpHjx66ZcuW+pFHHtGALl68ePz2sbGxOiAgQAN66NChydrr1auXBnSnTp0SLf/44481oPv27at9fHx06dKldbdu3XSbNm20u7u7BnS7du10bGxsov3mzJmjAd2vX79Ey8+cOaMBXb58+WQ1nD59WleuXDn+O7Vu3Vr36NFDN27cWOfLl083a9Ys0fbNmjXTgN60aVOi5ZcuXdKenp4a0E2bNtVBQUF6wIABevny5en/wWqty5cvn2K7CY0ePVoDumLFijo6Ojp+eb9+/TSgP/744xT3++mnn7SLi4sG9Ouvvx6/vGPHjhrQSindvXt33a9fvxRf27Zts+k7OEJG/g0Ce3VqP99TW2HEC5hqDYZwYBXwM3DXumwp4GRjO37W8Entddja5klAJdk3LrS2p7LvZ7Z+n7wcWlprHbxqrK4z+0FwNfv+Ub39rzWG1CKhlb7w8HD97bff6qeeekoXKlRIu7i46BIlSui6devq4cOH6x07dsRvO2bMGA3oWrVq6YiIiGRthYaG6ipVqmhAT5kyJX55XGh9/PHH+tSpU7p79+66aNGi2tXVVfv5+enx48frqKioZO3ZE1paa3337l09btw4XadOHe3l5aU9PDx0hQoVdPfu3fW6desSbZtaaGmt9W+//ab9/f11gQIFtFIqzSBJypbQCg0N1cWLF9eA/v777+OXpxdaWqccXHHfJb3XnDlzbPoOjpDrQgvobP2DvARUSbC8eIKQec1Bxzpkbe/9FNbFhVZQZo+T10NLa63nrf1C100QXE2/f1Rv2rP0odchoZU9JAwtkbc4KrSy0z2tEdb3d7XWJ+IWaq2vAHFP3r2nlMpUzUqpRsCjQCwwNzNtifT1fuYd3ikzCA/rILu3nE188PcHrN+5wODKhBA5UbYILaVUGaAuEAUk68Kktd4CXARKAJkd3K6/9X2d1tr+R+6Fzbq1Gsb7FV6NHx3+jpOJUUc+YcWWWQZXJoTIabJLh/3a1vdDWuvwVLbZA5S2bvu7PQdRSnkC3a0f0/uJ2VwpVRPwAq5guce1QWtttufYeV0H/5dxdXLnkxMTCLVOJPnJ6a+4H3mXHk+/kX4DQghB9gmtuDFdzqWxzb9JtrVHV8AbuIqlo0dank9h2WGlVA+t9T+ZqCHPatM0CBcXd8YcHsttJxPhJhMTLs7i/spQ+rf7yOjyxEMwatQoh4+ULvKWbHF5EMvZDMC9NLaJm4XOOxPHibs0GKK1jk5lm/3AMKC6ta5SQFvgAJZ7Yb8qpewbEE3QqmEPPqs9niIxlhPWKJNiyo1FTP0p7THvhBACsk9oxQ24pbPsAEpVBp6yfpyd2nZa62+01lO01oe11ve01pe01quBJ4BdQDEedBpJ6TiDlFJ7lVJ7U3riX0CT2s8y8cnplIy2/HXHKMV3oWuZsOBlgysTQmR32SW0Qq3vXmlsE7cuNI1t0hJ3lrVTa53hR9u11lHAZ9aPqU53qrWeqbWup7WuV7RoUTvKzBvqVHuKyS3mUjbKElxaKUKidjBqbg+DKxNCZGfZJbTOWt+TDyH9QNkk29pMKeXEg3tUmemydtT6LpcHHcCvYl2mtllCpcgHI1v/xCGGz2qN2TryuBBCJJRdQusv63t1pVRqM7bVT7JtRgRiCZp7wEI79o8TN6tdWJpbCZtVKO3H9I4rqRb5YKqEX5wvMHRWcyIi7xtYmRAiO7IrtJRSJZVSrZVSQ5RS71pfQ5RSbZRSJTPantb6PPAn4Iqlh1/S4zUDygCXgZ12lDzA+r5Qa52ZwOlmfZepeR2oZNHyzOz5K3UiHky9vt3tFi8H+3MnLPngrUKIvCtDoaWUaqeU2gVcwNJlfArwqfU1BVgJXFBK7VRKtc1gLXH3i76wdpqIO2YxYJr14+cJn5NSSn2mlDqqlPqMVCilimDp/QfpXBpUStVSSrW1Xk5MuNxZKfUmll6FAF/b9I2EzQp6F+G7oC00iSwQv2yfeziD/i+AS9fPG1iZECI7sTm0lFLfAMuw9KKLwdIFfCWwwPpaiaW7eOgQdyEAACAASURBVAzQAFiulLL5h7vWegkwHcuoF/8opVYqpZYCJ7B0NV8GfJtkt5JAVet7avpiOYM7qrVO76FkX+v3uGoN3sVKqXVYnh+baN3mXa31L7Z+L2E7dzdPpg7YTGD0g1uGh91iePnnNhw/t9/AyoQQ2YVNoaWU6oHlLOMqEAQU1FrX0Vp30Fr3sb46aK3rAgWBF6zbDlNKdU+14SS01kOA3lguFTbDci/qJPAK0Flrbc/d+Res76l2c0/gADAJOAaUA9pZ67gPzAGe0FqPt6MGYSOTkzNfDlhLZ6rHLzvtCkM39GHn3/K7ghB5nbIMqJvORkptBeoBj2mtT9nUsFJVsITAHq11s0xVmUPVq1dP79271+gycqyvFg4hJHwrscrSu7BArJn3qrxN26ZBGW7ryJEjVKtWzcEVCiFslZF/g0qpfVrreimts/XyYE1go62BBWAdqf1X4HFb9xEioTe7T+P1It1xM1t+sbrjZGL0yS8JXjXO4MrEw+br64tSirNnzz6U4ymlUEqlv6F46GwNLScgMt2tkou27iuEXYLafsjHld6gQKyl/02ESfHN9QVMXDjU4MpEThUUFIRSiuDgYKNLEXawNbSOAgHWnng2UUoVBQJ48ECuEHZp99QAJtT9kuLRluCKVYrgiK28P7ujPIQsssSRI0c4ciTDA+eIh8DW0AoGCgC/KaUap7exUqoJ8BuWwW3n2F2dEFaNHm/NtJbzqRj1YNlKp5MMmdWMe/ftHdlLiJT5+fnh5+dndBkiBbaG1jQsXcFrAFuVUueUUkuUUhOVUmOUUqOt/71EKXUO2IJllPSVWutpaTUshK0e8a3Fdx3X8FiES/yyHW53GDivGRevpTWrjUjPH3/8wdtvv029evUoXrw4rq6ulCpVii5durBr165k248aNQqlFKNGjeLKlSu89NJLlClTBjc3NypUqMB7771HREREsv1CQ0OZOXMmHTp0oHLlynh6euLl5UXt2rUZN24c4eGpTaeXmNaaKlWqoJRKsb44nTp1QinFtGnTOHv2LEop5s61TFj+wgsvxN+7Snq5MK17WtHR0cycOZPmzZvj4+ODm5sb5cqVo23btsyfP9+m+kUmaK1temEJuLewTIhoTvCKtb4SLrts3dZka/u58VW3bl0tHO9eeJh++bsmukZwjfjXszNr6P1Ht6e6z+HDhx9ihTlPQECAdnJy0jVr1tRt27bVnTt31jVq1NCAdnJy0osWLUq0/ccff6wB3b9/f126dGldqlQp3aVLF/30009rT09PDeh27dolO862bds0oIsVK6abNm2qu3fvrlu2bKm9vb01oJ944gkdHh6ebL/y5ctrQJ85cyZ+2ddff60B3bdv3xS/04ULF7Szs7P29vbWd+/e1deuXdP9+vXTlSpV0oBu3Lix7tevX/xr27Zt8ftimXEiWZs3b97UjRo10oB2c3PTLVq00D169NBPPfWULliwoC5fvryNf+J5T0b+DQJ7dSo/V23q8p6QUsoENMIyg3AFLKOvKyyjr5/F8ozVTi0z/EqX9yxkjo3lo5BuLDcdj19WOMbMiGofEPhkz2Tbp9Xd9rG5j2VZnVntn36OmY903bp11K5dm+LFiydavnLlSjp37oy3tzfnz5/H09MTsJxpjR49GoCBAwcydepUXF1dAcuf9RNPPEFYWBjbt2+nceMHdxQuXLjA8ePH8ff3x2R6cKHn9u3b9OzZk3Xr1vH555/z7rvvJqrD19eXc+fOcebMGXx9fQG4c+cOZcqUITo6mgsXLlCkSOJb7h999BFjx45l6NChfPvtg3EJgoKCmDt3LnPmzCEoKCjFP4+4s6ykPx/bt2/PihUraNSoEUuWLKFUqVLx6yIiIti0aROtW7dO+Q85j3vYXd7jaa3NWusdWutvtdZvaa1f0loPsv73FOu6PB9YImuZnJz45IWfeNEzAGfrD5YbziY+PPYJs1eONri6nOeZZ55JFlgA7dq1o2vXrty8eZNNmzYlW1+2bFkmT54cH1gA1apVo2/fvgBs3Lgx0fZlypShRYsWiQILoGDBgkyePBmAJUuW2FRzgQIF6Nu3L5GRkcyenXjsgOjoaP73v/8BMGTIEJvaS8/+/ftZsWIFXl5eLF++PFFgAbi7u0tgPQTORhcgRGYM6/oNJTZM4uvzMwlzMhFuMvHNjcVcCDnOB71DMDnJExe2un79OqtWreLgwYPcvn2bmJgYAA4ePAjA8ePHefbZZxPt06JFCzw8kk/MENeJ4b///ku2TmvNjh072Lp1KxcuXCA8PDzhbQiOHz+ebJ/UvPLKK0yfPp0ZM2YwfPjw+DBcunQply9fxt/fn0cffdTm9tKybt06wHK2JXPlGUdCS+R43Vq9Rsk/KzD6zxFccTGhlWKx/psrswIY33c1+Tzypbm/oy6x5WTfffcdb775Jvfvpz4dzN27d5MtK1euXIrb5s+fHyBZZ4wrV67QqVMnfv899WFAUzpOah599FFatmzJr7/+yrp162jTxjI/67Rplv5fQ4c67nm+c+csnX2kV6Gx7J5PSylVWSn1P6XUSaXUfaVUbCqvGEcWLERKmtZ5ju+eXkjVyAf/S291u8GAkCacv2zzQC550t69exk8eDDR0dFMmDCBo0ePEhYWhtlsRmvNiBEjgOT3d4Bkl/nSM3DgQH7//XcaN27Mhg0buHr1KlFRUWitiYy0Z/wCePXVV4EHQXXo0CG2bt1KqVKl6NChg11tiuzL3vm06mHpcNEfqAi4Y+mMkdIru0w0KXK5SuVqMLv3JhpFescvO+Qew4ur2hMZlbz7tbBYsmQJWmuGDRvG8OHDqVq1Kvny5YvvjHDy5EmHHOfevXusWbMGJycnVq1aRcuWLSlatCguLi6ZOk7btm2pUKECa9eu5ezZs0ydOhWAQYMG4ezsuItJ5ctbJlY/duyYw9oUGWdvoIzH0mtwEVAH8NZam1J7OaxaIdKRP58PMwZso725Svyyiy6KW1E3uXnnioGVZV83b1om2ixbtmyyddeuXWPDhg0OOc6dO3cwm814e3tTsGDBZOvtfcbJZDIxZMgQzGYzEyZMYN68eTg7OzNo0KAUt4/rNBJ3z85WgYGBACxfvpzr16/bVavIPHsDpQFwRGvdU2u9X2t9z5FFCZEZlp6FSxni3Tp+sF0zcCn6OpdunE3xMldeFnePJiQkhLCwBxN7h4aG0r9/f27fvu2Q4xQvXpxChQpx+/Zt/u///i/RunXr1vHVV1/Z3faAAQPw9PRk2rRphIaG0rFjR0qWTHmavdKlLfO1ZXSYptq1a9OuXbv49i9dupRofUREBGvXrrXvCwib2Rta4VimHREi2xrcaTyjKr+JT8yDJzBu6nv8e/04sTJmYbwXXniBsmXL8ueff1KxYkU6depEx44d8fX1Ze/evfTv398hx3FycmLkyJEA9O7dmyeffJJevXrRoEEDWrduzZtvvml324UKFaJPnz7xn9PqgNG+fXtMJhPffPMNgYGBDBgwIP5eW3qCg4OpX78+27dvp2LFirRq1YpevXrh7+9PyZIlGTx4sN3fQdjG3tDajeVelhDZWtum/ZnWbA7OCU6uwlQMZ28eJSLKtiGDcrtChQqxd+9eBg0ahJeXF6tXr2bv3r106tSJP//8M8XLhvZ66623WLJkCQ0bNuTQoUOsWrUKJycn5s2bx7hxmZtyplWrVgBUr16dZs1Sn8KvVq1aLFy4kPr16/P7778ze/ZsZs2aZVNXex8fH7Zt28aUKVOoU6cOu3fvZunSpZw5c4amTZvy+eefZ+o7iPRleEQMAOugub8BPbXWSx1eVS4hI2JkH4cOHSRfMSfuqQdnXc5aU9KjFPm9fAysTDhKx44dWbZsGdOmTZMznmzIUSNi2NW1Rmu9QynVA/ifUqoj8AtwAcutg5S232rPcYRwFJPJifJF/PjvxmluY+lJGKMUFyIuUST6PsUKlTG4QpEZ+/btY8WKFRQuXJjnn3/e6HJEFspMf1BX4D7Qy/pKjc7kcYRwCKUUpYtUwu32f1yLuYkZhQauxd4h8lo4pQtXyvBzR8JYAwcOJCwsjDVr1mA2mxkzZgz58qX9MLnI2ewKE6VUZ2A+lntiN7AMlBuW1j5CZBdFCpbC7Z4nl+5fJNo6+8RdFUX0jaOULlARN1d3YwsUNps1axYmk4ny5cvz0UcfOWycQZF92XsG9D6WB4eHADNlgFyR03jnK4ibqwcXbp8iXFnu64Yrzdk7JynhUYoCcp8rR5DHF/Iee6+F+AE7tNYzJLBETuXq4kaFIn4UwC1+WYxSXIy4xJWb5w2sTAiRGntD6w6WjhdC5GhKmShTpDLFnQtiwvJbuwaum+/y77Vj8jyXENmMvaG1HqivUpuPWogcpkjB0pT1LItLgqtNoSqGMzePcj9CbtcKkV3YG1ojAW/gS6WU9AwUOUJ69z+8PAtQoVAVPBMMlxmp4N+ws9y4czmryxMi13LkvUd7A2cAsAZ4HeiolPqN1J/T0lrrsXYeRwiHcHJyIjo6OtEMuylxcXbFt4gfl26c4RaWETNiUVyOvkH49TBKFa6ISUm3eCEyIjo6GicHTchq74gYZiyX/tO6PBi3Xmut8+T0sTIiRvZx6dIlXFxcKFKkiM373Lp7lSuRV4lNcBXcQyvpFi9EBl2/fp3o6OhUBzFOyuEjYgBjAOlrKnIMHx8f/v33X8Ayq66Liwvp3ZItlL8YHpFeXLx7hgjrpnHd4ou7F6egt0y5LkRqtNZER0dz9+5dbt26leos1xll7zBOoxxydCEeEjc3N8qVK8fNmzc5e/ZshnoFag137l0nnAfzL13kCh64UiBf4XTDT4i8ysnJCW9vb8qVK4ebm1v6O9jAptBSSvUBVmmtHTOxjhAGcHNzo2TJkjZfokhq9spR/O/aYsKcHtzTqhXhxpi286lQuqqjyhRCpMHWO8ohwBWl1Dql1EtKqRJZWZQQ2VH/dqOY2mgaVSIfnFntd49kwLpO/PTbDAMrEyLvsDW03sIyh1ZLYDpwQSm1TSn1plKqQpZVJ0Q2U6daM+b13cHT0aXil11zNjH232/5MLgr0THRBlYnRO5nU2hprb/WWjcFSgGDgY3AE8CXwEml1H6l1IdKqRpZV6oQ2YOnhzcTB/7CW4W7kj/W8pRHrFIsU0fpN6sBx87uN7hCIXIvu7q8AyilCgDtgE7A04Anlh6Fp4GfgGVa610OqjNHki7vud8/J3YyetPLHHN78IhioVgz/Yv1IqjtSAMrEyLnSqvLu92hleQAHkBrLAHWBiiIJcAuAT8Dy4BNeW1wXQmtvCEi4h4fz+/GGud/Ey1vEVWcMT0XyYjxQmRQWqHlkEf7tdbhWuulWus+QDEsAfY94AQMxTJWofzaKXIld/d8fDFgNSNLBFE45sHvZb+5XqH3gmZs3rfcwOqEyF1sCi2l1FdKqR62bKu1jtFa/6K1fgnLPbBmwCQsZ11C5Fo9At8i+Jkl1Il4MFrGOVcY/vdIPps/ELOMGC9Eptl6pvU6lvtWACilYpVSs9LbSVts01q/qbX+3t4ihcgpfEtXY87AXfR1eQI3s+XSe6RJ8X8xfxD0fUNOnj9kcIVC5Gy2hlYs4JLgsyLtcQeFyLNMTk6802sWX9f8hPJRD5b/5R5B/w3dmL1Sxo8Wwl62htZVoJbMnyWE7ZrW7cD/9dxMQHTx+GW3nEx8fXMRQ2Y24/KNiwZWJ0TOZFPvQaXUfKAncA44A/gDl4GjNhxDa60DMlFjjiW9B0WceWs/5/v/fuCG84PfE0tEawZXGkan5oMMrEyI7CfTXd6VUmWApUCKjaRDpiYRAjh/6TijVjzPbvd78ctMWvN0bHk+7rkAL8/8BlYnRPbhsOe0lFK+QDlgM7AO+MKW/bTWW2w+SC4ioSWS0mYzU5e+xfy76xMNvFshCl6rNZqA+p0MrE6I7MHhDxdbJ4EM1lr3z2xxuZmElkjNPyd28ummIRx0ezDdibPWtNF+fNArBA83TwOrE8JYWfFwcQXgbftLEiJve6xKI37ov5seplrxXeNjlGKF6Rg95zbkt70/G1yhENmTXaGltT6ntb7h6GKEyEucnV0Y2fcHvqn9eaLpTk65aYYf/JAPgrsSHnnfwAqFyH4yNfagUqoslhEvSgHuqWymtdZ58sEUuTwobBUZFc5nC4JYaT5ElOlBgFWKVLxedyz+ddsbWJ0QD1dW3NNyBr4FBvLgIeOkz3Bp67IM9R5USvXCMv1JTSxjFx4F5gDTMzLgrlIqGOiXxibHtNZ+WVmHhJbIqO1/ruCrvR9wwu3Bv0sXrXmWRxnRcw6ebvkMrE6IhyOt0HK2s81RwCAgBlgDnADC7GwrnlJqKjAEiMAyZ1c0EIAlIAOUUl211hkdwG0HcDKF5amOhZhFdQiRriZ1nqN+jVZ8uqAfK/VhopUiWimWcYQDIQ0ZUn0EzzzZy+gyhTCMvWda5wAfoLHW+m+HFKJUZ2AJloeWn9Jan7AuLw5sAqoBr2utJ9nYXjCWM60XtNbBRtQhZ1oiM7btW85Xf37ISdcH/0adtCYgpiwju4bgU6CogdUJkXWyovdgMWCLowLLaoT1/d24oADQWl/BcpkO4D2llEOmU8kBdYg8rmnd9vzY7w+6qOrxPQxjlWK9ywV6LG7Oj+tt+v1NiFzF3h+8/wKRjirCOuJGXSAKWJx0vfXh5ItACaCho46bXesQIo6bqwcfP/8j0+t/zaORD67mX3JRjLv0PYNnNuP85dMGVijEw2VvaP0INFNKeTmojtrW90Na6/BUttmTZFtbNbfOBzZTKTVWKRWYxllSVtYhhN3q12jFggF7CHJrjFfsg35A291u0md1O2b8PFLm6xJ5gr2h9SlwDFitlHrEAXVUsL6fS2ObuLnMK6SxTUqeB94AXgQ+wDL81D9Kqccech1CZIrJyZm3esxgTvMQ6kV4xC+/6Wxi6t0VPP99A/Yd2WZghUJkPXsfLo7EMilkIeCQUuqkUmqzUuq3FF4bbWgy7oztXhrbxPVO9LaxzP3AMKC6tf1SQFvgAPAo8KtSqrSj61BKDVJK7VVK7b127ZqNpQphO78KdZkz6A+GFXwOn5gHZ10H3CN5addgPprbnXvhaf0vLETOZVdoKaWKANuxBIITUBF4CsuUJSm90m3S+m7/k85JaK2/0VpP0Vof1lrf01pf0lqvBp4AdmHpTDIiyW6ZrkNrPVNrXU9rXa9oUendJbKIUrzYfhzz263EP6oISj+YJflnDtNtXkOW/Dbd4CKFcDx7Lw9+DjwOHMdy6a0d0DyVVwsb2gu1vqd1jyxuXWga26RLax0FfGb92MaoOoRwhDLFKjLlxU18WvENKiSYJflfVxh9fhovz3yK0xdtmfZOiJzB3oeLn8XycG5DrfUdB9Rx1vpePo1tyibZNjPi/hUnvTz4sOsQwiHaPjWAlk/04MvFg1gRs59wk+X30R1ut3j+l850zNec1zp/jbOzi8GVCpE59p5peQO/OyiwAP6yvldXSnmksk39JNtmRmHre9JRPB52HUI4jLt7Pj7oO5/ZTWZRP0FHjTtOJoIjttBtdj3W7phvYIVCZJ69oXUE2ztEpEtrfR74E3AFuiZdr5RqBpTBMkrFTgccspv1fU/ChQbUIYTD1ajSkNkv7Wa4TxeKJ+ioccLNzLsnPmPoTH/OXDxmYIVC2M/e0JoK+Duou3ucuPtMXyilKsctVEoVA6ZZP36ecLBapdRnSqmjSqnPErSDUqqWUqqtUsopyXJnpdSbWHoVAnztiDqEyI76tfuYRZ030ia2HC7WjhpaKba63aD3L534bP5AIiMjDK5SiIyxt8t7MPANsFkpNcA6kkSmaK2XANOxjDbxj1JqpVJqKZbBeB8FlmEZsDahkkBV63tCvsBK4KpSaqdSarFSah2W568mWrd5V2v9i4PqECJb8ilYgi/6r2ZG3YnUinSNXx7qZOL/Yv6ga0h9lm6aaWCFQmSMvQPmZuTRe621trnDh3VKkKHAYzyYEmQ2KUwJkmBQ3Lla66AEyysAr2Hp3l4eyz0sDVwAtgFTtdb7HFVHamTAXJGdaLOZuavH8n9XFnHJJfHvq40iC/BGq2+pVqGWQdUJ8UBWzKeVoUtjWus8ObishJbIju6G3WDC4hdZy3EiE0w46WE204pqDO/yHYW8C6fRghBZy+GjvGutTRl5Za58IYQj5fcqzNgXljL7yRnUj/SMXx5uMrHCdIwui55iypK3iY2RsQxF9iOBIkQeVbNqE2YP+oORJYIon+DB5KvOJmbeW0e32XVYvnWOcQUKkQIJLSHyuB6Bb/FTv90879KQgglGkD/uZuaDM1/x4ndPsv/YLgMrFOIBCS0hBG6uHrzd638sfm4Nz8SUju8iD7DLPZSBvw/k3dntuXLzPwOrFEJCSwiRQIki5ZkwYB3/qz8p0agakSbFGqfTdP25FZ/NH8i9CBlFXhhDQksIkUzd6gHMfmk3Y8oMolKC+123nC3Pd3We14DpP4+UzhrioZPQEkKkqmPAqywJ2sdAj2YUTTAk1EUXxbS7K+g6uw6LN05LowUhHEtCSwiRJmcXV17r9i1Lu22hM9Xwik08nuGYC9Pp+119Nu1bYWCVIq+w6+FiYRt5uFjkRuf/O843a15lk/NFotWDh5NNWtMoyofBzSfweJUGBlYocjqHj4ghbCOhJXKzv49tY+qW99jpegedILxczZqmsaUZ0moij5SvYWCFIqeS0DKIhJbICzbtXsSsvz7ngHt0ouUeZjPNzBUY1mYyZUtWNKg6kRNJaBlEQkvkJYvWT2TB2bmcdEv8M8Ur1kxz5cerz02mZOGkk4ULkZyElkEktEReY46NZe6aMSy9vJSzronXFYg1E+D0OK91mIRPgaLGFChyBAktg0hoibwqJjqK71eOZNnNtVx0UYnWFY4x08q9Aa92+Jr8+QoYVKHIziS0DCKhJfK6yKhwZix7m1V3N3E5yRxeRWLMtHRvyNDnvqSgdyGDKhTZkYSWQSS0hLC4f/8u3y57g7Xhu7junDi8CseYaeFal1eemyiXDQUgoWUYCS0hErsTep0py4axPvoAt5wSh1ehGDPNXR7nlee+pmjB4gZVKLIDCS2DSGgJkbLbd64ybeVbrI/8kxtJzrwKxJpp7lSDIW2/kt6GeZSElkEktIRI292wG0xfPpxfInZzLUl45Y8100xVY8izEylTrLxBFQojSGgZREJLCNvcu3+b6cvfZt3937mSJLzymc00MVdgYMA4/HwfN6hC8TBJaBlEQkuIjAmPCGPG8rdZG7qVS0l6G7qaNQ2ji9K74Qc8WTPAoArFwyChZRAJLSHsExl5n5kr3mX1nU3JnvMyaU3dSC+61BhGm8a9DKpQZCUJLYNIaAmROdFRkYSsHcPqKys54Zb8Z1WNCBfa+valZ8thmJycDKhQZAUJLYNIaAnhGNps5qffJrPsVEiygXkBKkYqAou2ZcCzH+Pm6mZAhcKRJLQMIqElhONt/GMBPx6YxG7XMMwq8aXDktGap9zqMfDZTynhU8qgCkVmSWgZREJLiKzz15FNBG8fww7na0SaEoeXV6yZhuay9G3yIXX8GhtUobCXhJZBJLSEyHpn/v2Hmb+OYAtnCE0yyoaT1tSO9KJd1Rfp2Kw/KsmZmcieJLQMIqElxMNz685l/rd6BL/d35OsxyFA5UgTLYoEMqDNKDzdPQ2oUNhKQssgElpCPHyxMVHM/+Uz1l5cxkG3mGTri8Zomjg/xsBnPqNccd+HX6BIl4SWQSS0hDDW5t2LWLh/Mn+43iY6yaVBD7OmXnRROtUeRsv6HQ2qUKREQssgElpCZA8n/93PnI0fsVWf4naS+14AfpFOPFX4aYKe+RDvfN4GVCgSktAyiISWENlLaNhNvl/1PhtDt3PONfl9r0IxZhqqSvRq8gG1HnnCgAoFSGgZRkJLiOzJHBvL8s3TWH1qPvtcw4hJcunQSWsej8xHK98e9AwYhpOzjLbxMEloGURCS4js7/iZvfyw+RO2m08km1UZoHS05km3ugQ9PVo6bjwkEloGkdASIueIiLzHvHWf8OuVtRxyi0223t2sqRftQ6uqQXRo0k/GOsxCEloGkdASImf6ff9KFu35hp3Ol7lvSn72VS4KnnCrTe8W71O5jJ8BFeZuEloGkdASIme7cesic9Z+zOZ7u1LsuOGiNbUivWlevis9Wr6Ki7OLAVXmPhJaBpHQEiJ3MMfGsm7HHFYfDWGPyw3CUzj7Kh6tecLJj25N3qVWlfoGVJl7SGgZREJLiNzn+o2LhGwYy467v3M8hTm+TFrzWKQ7TYo/S99n3iGfez4DqszZJLQMIqElRO62dc8Slh2YwS7TpWSD9YLlua+6uiyBNV4ksEEnGbDXRhJaBpHQEiJvCA27ybz149h6/bcUxzsEKB8Fddxq0q3JcGpUrP2QK8xZJLQMIqElRN7z1+GNLNr1NTv1GW6k8NyXSWuqR7rxhI8/vVq+Q7FCxQ2oMnuT0DKIhJYQeVdU5H2W/PYNm8+v5E+X0GQTVQJ4ms3UivahqW8XujUfjKuLqwGVZj8SWgaR0BJCAFy98S//t/Fzdt3+PcUHlwGKxpipTQWerfUyzes8m6fvf0loGURCSwiR1MHj21my8xv+iD7KhRQmqwSoGKmo5f4YHRq8Su2qDR9yhcaT0DKIhJYQIjXabGbDrh9YcziEvabL3Emh9yHAI5FOPO5Ri46NX+exirUecpXGyFGhpZTqBQwGagJOwFFgDjBda222sQ0X4CmgDdAYKA8UBq4BO4FvtdabU9k3GOiXRvPHtNY2jdsioSWEsEVEeCgLN05k639r+Mv1frIJKwGU1lSNcqF2vnp0bvoGVcs9akClD0eOCS2l1FRgCBABbASigQDAG/gZ6Kq1TvmCcOJ2WgIbrB8vA/uAe8CjQA3r8rFaV89c5QAAEJxJREFU649S2DcYS2jtAE6m0PwlrfUIW76PhJYQIqP+u3KKRVsmsuf2Lg65RhGbQoA5aY1fpBu18zega7O3qFiqkgGVZp0cEVpKqc7AEiwh85TW+oR1eXFgE1ANeF1rPcmGtlpgCb9JWuttSdZ1B+ZjOYtrobXelGR9MJbQekFrHZyZ7yShJYTIjHMXD7Fk6zfsCd3DYdcYdAoB5qw1j0a6U6dgE7o2f5NyxcoZUKlj5ZTQ2gvUBfpprUOSrGsGbMYSaKVtvUyYxrG+BwYAs7XWA5KsC0ZCSwiRzZw6t58l279m3739HHFL+Uegq1nzaJQntQo2omOTV6lYuvJDrtIx0got54ddTEqUUmWwBFYUsDjpeq31FqXURaA00BD4PZOH/Mv6XiaT7QghxENRqXwt3i0/F4AjJ3eydOcU9oUf5ESC8Q+jTIr97uHsj/iNHzZspGqkKzXy1eLZ+i9Rp2oDo0p3qGwRWkDcmCaHtNbhqWyzB0to1SbzoVXF+n4pjW2aK6VqAl7AFWA7sCGzZ3lCCJFZ1So3YmTlRgD8fWQzy3ZP48+oI5xK8GxyrFIcdo/mcOweFu3aQ8UtiuouVQmoGUTz2q0xpTBSfU6QXUKrgvX9XBrb/JtkW7sopUoAQdaPP6Wx6fMpLDuslOqhtf4nMzUIIYSj1KzmT81q/gDs/Wcta/6aw98RRzmWZAT6026a0xxl5cH3KPXXu1SnPI2rdKHdk31wdck584Bll9Dysr7fS2ObMOu7t70HUUo5A/OAAsBGrfXKFDbbj6W34UYsIZofqAOMAx4HflVK1dFaX0zlGIOAQQDlyuX8G6JCiJyj3mOtqfdYawBOnN7Dsl3TORD2F4dco4lJ0InjPxfFf/zLhrNfMfnkl1SPLU79Mm3o4v8y3p5eqTWfLWSLjhhKqZHAJ8A8rXXfVLYZB7wPzNRav2TnceI6YJwHntBaX87Avq7AFiz31KZqrV9Jbx/piCGEyA4uXz3Nz9sm8//t3XmQVeWZx/HvT8AlmsQYRMboKC5AXMZRY+KCwkRTQ6JSLmgSdaJmc0FlKlNMKlWZGcr8EcdlRiOiZWkkE2uCUcG4m2gCZtxRSgjEdYILA+4IItJCP/PH+7YcL/fevn3pvkv371N16u173vfc896nu+7T55z3vGfeWw+zcMjqsg+xBNhmfSd7fbgt+372EI45+PvssdPIBvc0afmBGMCqXFZL8V11q6q0qUjSFaSEtRw4sicJCyAiOiT9FPgN6aZlM7O2MHzYbpxz4uUArFr1Jrc9OI1Hlz7AgsFvs6IwE8d7gzbj8UEreXz1fdxw/73s3jGYUUP2ZMzok/j7g05k8OBBzfoIH2mVI60JpGQwPyIOqNBmFnA8cH5ETOvh+18G/IA0I8a4iFhcZz9HAs8CHRGxRXftfaRlZq2sY+373P3QdTz4v7ezQMt4rcyjVLpst66T0euHst/2R3DsYWex87C+G3zd8vdpSdqZNNCiA9i23AhCSa+QhqiPiYiHevDeFwNTgLdIR1hPb0I/DyGNXHw7Ij7bXXsnLTNrF7F+PXPn3cTv/3wzizte5LnNO8vezAzphuaRazdn1JZ7MXbvUxi3/3gGVZg7sR4tn7QgdZI04KHXbi6WdBHwQ+AdUsKa380m3b3ffwL/CNwXEeO7a++kZWbt6pVli7nz4WuZ/9Zj/GnwSlZVSUo7fBiM6hzO/jseyYQx32PYtkM3ad/tkrQmkm4sXg4cHhEv5PXDSNM47UXJNE75GtPxwOzS+QAl/QT4MbACOCoinqyhD39LOpq7pzjHYR51eAFwCbAZMD4i7uvu/Zy0zKw/+HDtGn772H/x0It3sHjdSx+7H6zUFp3ByI4tGb3133DMl87jgD3LXvGpqh0GYhARt0i6mjTD+0JJ97NhwtxPAbcBpdey/goYlcuP5GtkP84vXwDOr/BAtWci4qLC611JE/O+Lek54FXSEPt9gR2BTuCHtSQsM7P+YsgWW3H0EWdx9BFp4PYLf3mCux6/ngXvPsWiIatZXRiNuHYzsXDLtSxc/wQdT/2qrqRVTcskLYCIOFfS/wCTgLFseDTJz+nBo0mA7Qo/fyEv5cwFiknraeAK4Iukx5nsDwQped1AGure7RGbmVl/tseIg5g84iAAPlizkrsfvp5Hl9zLos6lvLz5hgOEr37hzF7fd8ucHuyPfHrQzAaaPz07h3uf/AUvrHmNq799FxXOclXVFqcHzcys/e0zahz7jBrXZ+/fnjMmmpnZgOSkZWZmbcNJy8zM2oaTlpmZtQ0nLTMzaxtOWmZm1jactMzMrG345uI+JOkN0tOP6zEUeLMXuzMQOYabxvHbdI5hfXaJiO3LVThptShJ8yrdEW61cQw3jeO36RzD3ufTg2Zm1jactMzMrG04abWua5vdgX7AMdw0jt+mcwx7ma9pmZlZ2/CRlpmZtQ0nLTMzaxtOWi1G0imS/ijpXUnvSZonaZKkAfG7kjRK0mRJN0p6RlKnpJA0sYZt64pdf4q5pCGSjpR0maRHJS2T1CFpqaRbJI3rZnvHUDpf0q8l/VnSW5I+lPSGpPslnaYqTzV0/BogIry0yAJcBQSwBrgTmA2szOtmAYOa3ccGxODy/HlLl4l9Ebv+FnPgqELMluXPdBOwsLD+QsewagxfBTqAp4A7gJnAI0Bn/ky3AZs5fk36/TS7A17yLwJOLHzR7FlYvwOwONdNbnY/GxCH7wIXAycDuwNzukta9cauP8Yc+DJwC3B4mbqvA+vy5/o7x7BiDMcAW5dZvzewPH+mMx2/Jv1+mt0BL/kXAfPyH+i3ytSNLfxhb/QfXn9eakxadcVuIMYcuC5/rusdw7ri9y/5M/2349ecxedLW4CknYADSackbi6tj4i5wFJgOHBwY3vX2uqN3QCO+fxc7tS1wjHskXW5/KBrhePXWE5arWH/XC6KiDUV2jxR0taSemM3UGO+Zy6XFdY5hjWQNAI4O7+8o1Dl+DXQ4GZ3wAAYkctqM8K/XNLWknpjN+BiLmk4cEZ+eWuhyjEsQ9KZpFN0Q0hHpoeS/tH/aUTMLjR1/BrISas1bJPL1VXavJfLT/ZxX9pNvbEbUDGXNBi4Efg08EBEFI8UHMPyDgNOL7xeR7qm9R8l7Ry/BvLpwdbQdd+H59TquXpjN9Bifg1wJPAKcFpJnWNYRkR8NyIEfII0cvByYCrwqKQdC00dvwZy0moNq3K5TZU2XXWrqrQZiOqN3YCJuaQrgO+QhmsfGRHLS5o4hlVExJqIWBwRU4AfAfsB0wpNHL8GctJqDUtyuUuVNjuXtLVkSS57Grt6t2srki4DLgDeICWs58s0W5JLx7B7N+TyWElD8s9Lcun4NYCTVmvoGoa8t6StKrQ5qKStJfXGrt/HXNLFwA+At4CvRMTiCk0dw9qtIF3bGgxsl9c5fg3kpNUCIuIV0pQxmwMnldZLGksavbScNJ2MZfXGrr/HXNJFwBTgHVLCerpSW8ewR44gJawVwJvg+DVcs+9u9pIWYCIb7n7fo7B+GLCIATqdC7XNiFFX7PprzIGf5L6/AxxY4zaOYer34cCpwBZl6g4DXsyf6VLHrzmLHwLZQiRNB84h3W1/P/AhacTXp0iTdE6MiPXN62Hfk3QAML2wai/ScN/ngbe7VkbEwSXb1RW7/hZzSROA3+SX80hffOU8ExEXlWw74GMo6QzSdasVpKOg5aS/v91Jf4sAdwEnRckNwY5fgzQ7a3r5+AKcAjxEmuV5NfAkMIkBMvcYMI7ys7x/bOnN2PWnmJNuHu42fsAcx7Ds5xgBXAj8gXR7wBpSMllCmoj4uL6IQ3+JXyMWH2mZmVnb8EAMMzNrG05aZmbWNpy0zMysbThpmZlZ23DSMjOztuGkZWZmbcNJy8zM2oaTllmbknSGpJA0o0n7n5r3P7UZ+7eByUnLzHqVpHE5mc1pdl+s/xnc7A6YWduaBswkz3Zu1ghOWmZWl4h4EycsazCfHjTrRZJG51NjrxeebFvaZpCk5bnd3oX1W0v6Z0lPSFopaY2kRfnaUbVHslfqy6GSbs376sjlLZIOrrKNJJ0s6Z78GTokLZX0gKTzStpudE0rnxL8Q345Ntd3LXPy+z+fX1frx6zc5tyefm7r35y0zHpRRDwDPAZsD3ytQrPxwA7AvIhYBCBpJ+Bx4N9Jj19/BPgt8Bng34CHJH2m1n5IOgf4I3AC8DJphvKXgRPze32vzDabkx6FcRPwFeC5vN0zwD7AlTXs+l7gvvzza8AvCsu9kWbovirXl01Ikj4HHAusAn5Zwz5tIGn2NPNevPS3BTib9PiPWRXqf53rJ+XXAh7O664EPlFouxXpizuAGSXvc0aF9fuRnsm0nvTcp2LdN/L6DmCfkror8vs9C4wuqRsETChZNzW3n1qyfhzVH3/yaVJC+gAYWqb+wrz9tGb/Lr203uIjLbPeN5P0hXyMpKHFiny0NIGUNH6VV48HDgEeJT2l9v2u9pEeNHg28Dpwao1HWxeQrlfPjIibixURMZN09DQEmFzo1zDSgwg7gRMiHTEWt1sfEbfXsO9uRcS7pES8BfDtYl0+pdp1FDgdsxJOWma9LCJWkE6zDSE93K/oG6Qv69sjoutJzF2nEW+NiM4y77ea9BTiwcBBNXRhbC5nVKj/eS7HFdZ9Off3kcinLPvYtFyeLan4PXQCMJx0lLa4Af2wNuOkZdY3ZuTy9JL1p5fUA+yWy0tKBi58tLAhsW1fw74/l8u/VKh/saQdpOtokK5f9bmckO4nPSl4fKGq6zrXVRttZIaHvJv1ld8BrwIHSNo3IhZKGgV8CVhOGrDQZVAu55Ie617NSz3oQ6XHkqsH79GXrgSOIiWqu/NIyiOA/yMdqZptxEnLrA9ERKekXwI/Ig2Y+KdcAtwYEesLzV/J5c0R0RtHGEuB3UlHcC+WqR9RaNelKxmO6oX91+pO0tHgVyXtCkzK66+NiHUN7Ie1EZ8eNOs7M3J5ah5OflrJ+i735PKkXtrv3Fx+q0L9mbmcU1j3e9KIw0MlfX4T99+Ry6r/FOfrd9NJ30NTSPFZB1y7ifu3fsxJy6yPRMRzpKHsOwCXADtRuDer4DbgSdLNuNdI2q70vSTtJmlS6foKfkb68v+mpONL3uck4GRSgvpZoa+vA9eQvhNulTSyZLtBko6tcf9dR3B7SOrubM71wPukU4SfBGZHxLIa92MDkE8PmvWtGcChpGHoXa8/Jp9KPA64GzgLOEXS06RrYkOBvwZGkm7W7fb0YUQ8LWkyaYTeLEmPkU4T7gF8kTSs/byIWFiy6RTSacWvAYskPZL7MAzYN5fdXg+LiJckzQf2BxZIehJYCzwbEZeUtH1H0o3A9/MqD8CwqnykZda3bgLW5J+L92Z9TES8Skoo5wHzgb1Js1fsQ7oR91LScPCaRMR04HBgNuka1snArsAsYExEbHQKLiLWkmai+AfgwbzvicBoYAEbrjnV4gTSTdTbAd8EvgMcXaHt73K5KCLmVmhjBoAiKg0wMjPre5JmA8cB50bE1c3uj7U2Jy0zaxpJB5LmXHwH2CXfSG1Wka9pmVnDSboO2IZ0/Wwz4F+dsKwWPtIys4bLs3x0ku4Pmx4Rlza5S9YmnLTMzKxtePSgmZm1DSctMzNrG05aZmbWNpy0zMysbThpmZlZ2/h/IdM9zRG+Hp4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(num_heun[:end,1],num_heun[:end,2]/m0,'-',label='implicit Heun')\n",
    "plt.plot(num_rk2[:end,1],num_rk2[:end,2]/m0,'-',label='explicit RK2')\n",
    "plt.xlabel(\"velocity\")\n",
    "plt.ylabel(\"mf/m0\")\n",
    "plt.plot(250*np.log(0.25/(.25-.05*t[:end])), (.25-.05*t[:end])/.25, label = \"analytic\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. You should have a converged solution for integrating `simplerocket`. Now, create a more relastic function, `rocket` that incorporates gravity and drag and returns the velocity, $v$, the acceleration, $a$, and the mass rate change $\\frac{dm}{dt}$, as a function of the $state = [position,~velocity,~mass] = [y,~v,~m]$ using eqn (1). Where the mass rate change $\\frac{dm}{dt}$ and the propellent speed $u$ are constants. The average velocity of gun powder propellent used in firework rockets is $u=250$ m/s [3,4]. \n",
    "\n",
    "$\\frac{d~state}{dt} = f(state)$\n",
    "\n",
    "$\\left[\\begin{array}{c} v\\\\a\\\\ \\frac{dm}{dt} \\end{array}\\right] = \n",
    "\\left[\\begin{array}{c} v\\\\ \\frac{u}{m}\\frac{dm}{dt}-g-\\frac{c}{m}v^2 \\\\ \\frac{dm}{dt} \\end{array}\\right]$\n",
    "\n",
    "Use [two integration methods](../notebooks/03_Get_Oscillations.ipynb) to integrate the `rocket` function, one explicit method and one implicit method. Demonstrate that the solutions converge to equation (2.b) the Tsiolkovsky equation. Use an initial state of y=0 m, v=0 m/s, and m=0.25 kg. \n",
    "\n",
    "Integrate the function until mass, $m_{f}=0.05~kg$, using a mass rate change of $\\frac{dm}{dt}=0.05$ kg/s, . \n",
    "\n",
    "Compare solutions between the `simplerocket` and `rocket` integration, what is the height reached when the mass reaches $m_{f} = 0.05~kg?$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rocket(state,dmdt=0.05, u=250,c=0.18e-3):\n",
    "    '''Computes the right-hand side of the differential equation\n",
    "    for the acceleration of a rocket, with drag, in SI units.\n",
    "    \n",
    "    Arguments\n",
    "    ----------    \n",
    "    state : array of three dependent variables [y v m]^T\n",
    "    dmdt : mass rate change of rocket in kilograms/s default set to 0.05 kg/s\n",
    "    u    : speed of propellent expelled (default is 250 m/s)\n",
    "    c : drag constant for a rocket set to 0.18e-3 kg/m\n",
    "    Returns\n",
    "    -------\n",
    "    derivs: array of three derivatives [v (u/m*dmdt-g-c/mv^2) -dmdt]^T\n",
    "    '''\n",
    "    dstate = np.zeros(np.shape(state))\n",
    "    dstate[0] = state[1]\n",
    "    \n",
    "    dstate[1] = (u/state[2])*dmdt - 9.8 - (c/state[2])*state[1]**2\n",
    "    \n",
    "    dstate[2] = -dmdt\n",
    "    \n",
    "    return dstate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The number of time steps is 10000.\n",
      "The time increment is 0.1\n",
      "40\n"
     ]
    }
   ],
   "source": [
    "dt = 0.1 \n",
    "T = 1000\n",
    "N = round(T/dt)\n",
    "\n",
    "print('The number of time steps is {}.'.format( N ))\n",
    "print('The time increment is {}'.format( dt ))\n",
    "\n",
    "# time array\n",
    "t = np.linspace(0, T, N)\n",
    "\n",
    "x0 = 0    # initial position\n",
    "v0 = 0    # initial velocity\n",
    "m0 = 0.25    # initial mass\n",
    "\n",
    "\n",
    "#initialize solution array\n",
    "num_heun = np.zeros([N,3])\n",
    "num_rk2 = np.zeros([N,3])\n",
    "\n",
    "#Set intial conditions\n",
    "num_heun[0,0] = x0\n",
    "num_heun[0,1] = v0\n",
    "num_heun[0,2] = m0\n",
    "num_rk2[0,0] = x0\n",
    "num_rk2[0,1] = v0\n",
    "num_rk2[0,2] = m0\n",
    "\n",
    "for i in range(N-1):\n",
    "    if num_heun[i][2] < .05:\n",
    "        end = i\n",
    "        print(i)\n",
    "        break\n",
    "    num_heun[i+1] = heun_step(num_heun[i], rocket, dt)\n",
    "    num_rk2[i+1] = rk2_step(num_rk2[i], rocket, dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAAEaCAYAAACimQj6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd1zV1f/A8ddhbxQHiqiIC0XNgYOvmrhCHLlya+BIf440y4ZZpvWt/JpZappZKppmOcotVuZOc28cKQ5UEAQREGTc8/vjwkVkXeCyz/PxuI/bPZ8z3teEt5/P53zOEVJKFEVRFKUkMCrqABRFURRFXyppKYqiKCWGSlqKoihKiaGSlqIoilJiqKSlKIqilBgmRR1AaVaxYkXp4uJS1GEoiqKUKCdPngyXUlbK7JhKWgXIxcWFEydOFHUYiqIoJYoQ4lZWx9TlQUVRFKXEKDZJSwhRXwgxRQixRghxWQihEUJIIcQr+ex3qBDioBAiSggRI4Q4IYSYKITI9rvntZ2iKIpScIrT5cHxwBRDdiiEWAxMAOKBPUAi0Bn4BugshBggpUw2VDtFURSlYBWns4YLwBfAIKAOsD8/nQkh+qNNPCFAEyllTyllX6AuEAj0BSYZqp2iKIpS8IpN0pJS/iClfEdKuV5Ked0AXU5PeX9XSnntmXFC0Z7VAbyXyeW+vLZTFEVRClip/MUrhHAGWgAJwIbnj0sp9wN3gSpAm/y2UxRFUQpHqUxaQLOU94tSyrgs6hx/rm5+2hnU7dDrjPvBi8vX1HR5RVGUZ5XWpFUr5T3Luf7A7efq5qedwWz/cymjt73M36YPmbV3LIlJCQUxjKIoSolUWpOWTcp7bDZ1YlLebQ3QTkcIMTZlevyJsLCwHAN9XqRRPCGm2v8tF80Tmb5hVK77UBRFKa2K05R3QxIp77nd4TKv7XSklMuAZQAeHh657mdEpzcI/H4328yCAdidcBaPU6sZ3PzVvIakpEhKSiIiIoKoqCiSkpKKOhxFKfVMTEywt7fHwcEBExPDpJvSmrSiU95tsqmTeiz6mbK8tjOo9wf/RNiathy1Mgbgi3PzaOrSCjcHt4IastTTaDTcuXMHc3NzatSogZmZGUKInBsqipInUkoSEhJ4+PAhd+7coWbNmhgZ5f/iXmm9PHgz5b1mNnWqP1c3P+0Mysa2PGPcPsAlIRGABCEZv3MMkfGRBTVkqRcZGYmJiQlVq1bF3NxcJSxFKWBCCMzNzalatSomJiZERhrm91dpTVqnU97dhRCWWdRp+Vzd/LQzuNYdh/LaEzesNRoAwpOjmLpnComaxIIcttSKiYmhXLlyKlkpSiETQlCuXDliY7ObKqC/Upm0pJR3gFOAGTDg+eNCiA6AM9pVL47kt11B6TD8O94PS5t5fzL8NPOOzyvoYUul+Ph4rKysijoMRSmTrKysiIvL6imi3CnRSUsI8XnK4rqfZ3I4tex/Qog6z7SpDCxJ+ThHSqkxUDuDs6/giHPjD5kY+UhX9tPln9h0dVNBD13qaDQag1xPVxQl94yMjNBoDPMrs9j8FAshmgshjqa+gOYphz57rvxZVYH6Ke/pSCk3At+iXb3ivBBimxDiV+Aa0BDYjHYBXIO0KyjNfUbSMqkpXWOf6Mr+e/S/nAo9VVghlBrq0qCiFA1D/uwVm6QF2AGtn3mlPgdV97lyvUkpJwDD0F7y6wB4A/+iXfC2f1Yrtee1XUFxGfEt08ISqP9U+6Bxkkxi6r6p3I+5X5hhKIqiFLlik7SklPuklCKn13Nt/FLK/bLp9ycpZVsppZ2U0lpK2UJKuTiny3t5bVcQKjg6c9/jIxY+CKN8sjZfRsRHMGXvFJ4kPsmhtaIoSulRbJKWkr3mPqMINfdk/oNwTKT2meXAiEA+OPwBmsLPo4qiKEVCJa0SQhgZUX3EUmrHmzH9YdrzDn/c+oPFZxYXYWRKaeDi4oIQgps3bxZ1KFnat28fQgi8vLwyHBNCGOS+iZ+fH0II/P39893Xs7y8vPTqN3X8WbNmGXT80kQlrRKkYpXq3Gg1m4HRMQyLSluQY9m5Zey4saMII1OU0s3f3x8hBH5+fkUdSpmnklYJ06L7aE7admRaRCRtn6Q99zDz8EzOhZ0rwsiUkmzPnj0EBgZSrVq1og4lTwIDAwkMDMx3P59//jmBgYH07dvXAFEpBUElrRKotu9SHlGOLx6EUyt1qSdNApP/mqxmFCp5Urt2bdzc3DA1NS3qUPLEzc0NN7f8r81ZtWpV3NzcsLe3N0BUSkFQSasEKlexCsFtP8dWSr4JDcM+ZUbhw/iHTPxrIjEJMTn0oCjpZXVPK/VezL59+zh8+DDdunWjfPny2Nvb4+3tzZkzZ3R1V69eTcuWLbGxscHBwYHhw4cTEhKSYaxnL7WFh4czfvx4nJ2dsbCwoHbt2nzwwQc8eZK7WbHZ3dNKTExk2bJldOzYEQcHB92iyT179mTt2rXp6mZ2T8vFxYWRI0cCsGrVKt1YRXG5MDAwkNGjR1OrVi0sLCwoX748Xbp0YevWrZnWz+leX1b/358t/+OPP+jcuTP29vZYWVnRpk2bLMcrDCpplVBNuw7lWLnu1EhK4qtnZhRei7zG2wfeJkmjtt5QDGfbtm106NCByMhIvL29cXJy4vfff6dDhw5cu3aNadOmMWbMGMqVK4e3tzdmZmasXbuWLl26kJCQ+UamkZGRtG7dmvXr19O6dWu8vb0JCwvj008/pXPnzrlOXFmN0aFDB8aNG8eRI0do1qwZ/fr1o1atWhw+fJgZM2bk2Mcrr7xC27ZtAe0Zqa+vr+7Vrl27fMeor59//pmmTZuyYsUKrK2t6dmzJ02aNOHgwYP07t2bmTNnGnzM5cuX4+3tTUxMDN27d8fNzY1//vmHPn36sHHjRoOPp4/SujVJmdBg5GLufd2alvEP+Cg8gg8rVQDg0N1DzD0+l/dbv1/EEZYMLu+V3EksN+f0KJRxvvrqK9avX88rr7wCaJfFGj58OOvWraNfv36EhYVx5swZGjZsCEBERASenp5cvHiRX375hREjRmToc+vWrbRt25aTJ09Srlw5AEJDQ+natStHjx5l1qxZzJ07N19x+/n5ceTIETw9Pdm4cSNOTk66Y/Hx8ezduzfHPubNm4e/vz+HDx+mXbt2Bp9ZqI9z587h6+uLmZkZmzdvxsfHR3fs4sWL+Pj48Mknn9CxY0c6duxosHHnzp3Lzp076datm67sv//9Lx9++CHTp0/X/X0oTOpMqwSztXfg0UsL0UhBn5hYXnsUpTu27vI61gauzaa1ouhv8ODB6X5BGRkZ8c477wBw4cIFPv74Y13CAnBwcOD//u//ALJMDEIIvv32W13CAnB0dGTBggUALF26lPj4+DzHfObMGbZu3YqNjQ1btmxJl7AALCws0v3yLwwjR45Md3nx+deqVasybffpp5+SkJDA3LlzM8Ts7u7O/PnzAfjmG8OuMPf666+nS1gA77zzDvb29vz777/cvn3boOPpQ51plXANPX04cm4YnvfXMCkyiiATM/600e6qMvf4XJxtnOlQvUMRR6mUdM//4gKoU6dOtsfr1q0LwL179zLts0mTJjRu3DhDeceOHalWrRp3797l5MmTuktzuRUQEABA7969qVSpUp76MLS2bdum+3N73qFDh7h+/Xq6Mo1GQ0BAAEKILM9sOnTQ/owfOWLYzSd69uyZoczMzAxXV1dOnz7NvXv3qFGjhkHHzIlKWqVAc98vuDH3EK6am3weHs5t0+pcNQeN1PD2gbdZ2W0l7hXcizrMYquwLrGVZM7OzhnKbGxs9Dqe1dlSrVq1shzPxcWFu3fvEhwcnNtQdW7dugVgkFmFhjJmzJhsJ2/4+fllSFoPHz7k8ePHAFSuXDnb/sPCwvId47OySkh2dnZA1v9vC5JKWqWAuYUVov8ynq7vgQWJfBcSzCvVa/PQ6ClxSXFM2jOJtd3X4mTjlHNnipKJnLZ1KahtX9TK/JCcMjvY2NiY4cOHG7TvnLYLKY7b+aikVUrUcm/NUbeptLkyl4oaDSvuBTGkhitPNPGEx4Uz4c8JrO6+Gjszu6IOVVEAsl0yKvXY8/ehcqNmzZoAXLlyJc99FAcVK1bE0tKSuLg4vvnmm3RnuDkxNTUlMTGRmJiYDO0SExO5f7/kPddZ/NKokmetB03nnEVLAFwTk/jv3QhMjbQPi16Pus7UvVNJTE4syhAVRefs2bNcuHAhQ/n+/fu5e/cuNjY2tGjRIs/9e3t7A7BlyxbCw8Pz3A9o7+MAJCUV/qMkJiYmdOnSBSDX08xTVzi5fPlyhmO///57kXyf/FJJqxQRRkY4+a0gEu3ZVNf4cEZHpp1ZHQs5xsy/ZyJTnulSlKIkpWTChAlERaXNeg0LC2PKlCkAjB07FktLyzz336xZM3r16kV0dDR9+/bNcFYRHx/Prl279Oor9Ze/IZaKyouZM2diamrKlClT+PnnnzP8DGs0Gvbs2aObfJKqc+fOAHz88cfpnpe7ePEir7/+esEHXgDU5cFSpmKVGpxp/wXlD44DYGLEaUIr9ua3xNMAbL+xnSrWVZjSfEpRhqkovPzyy1y4cIHatWvj5eVFUlISe/fu5fHjx7Rs2ZKPP/4432P4+/vTrVs3Dh06hKurK+3ataNSpUrcu3ePs2fPYm9vr9fK9m3atKFKlSqcOnUKDw8P3N3dMTU1pW3btrrVMgqSh4cHq1evZtSoUQwZMoT33nuPhg0bYmtrS3BwMFevXiU8PJx333033UzO6dOns2HDBrZt20b9+vVp0aIFISEhHD9+nIEDB6LRaHQTVkoKdaZVCjXtPJh/KvbTfZ5+ZSfeldKmvf9w/gfWXV5XFKEpik758uU5evQoffv25ciRI+zatYsKFSrw/vvvs3fvXqytrfM9hoODAwcPHmTRokU0b96cY8eO8euvvxIUFET79u2ZM2eOXv2Ym5sTEBBAjx49CAoKYs2aNSxfvpz9+/fnO0Z9DR48mPPnzzN58mSsrKzYv38/27dvJyQkhObNm7NgwQImT56crk3t2rU5fPgwL7/8Mo8ePWLHjh1ERUXxxRdfsHr16kKL3ZCEulRUcDw8POSJEyeKZOz4JzHcn/cfamm0/4q6blSDeW3acuj+YQAEgvle8+lSs0uRxFfYAgMDadCgQVGHoaA9+xk5ciS+vr5FsrqEUjRy8zMohDgppfTI7Jg60yqlLKxsMBqwgjipvYFcW3ObIVdiaFKxCQASybsH3uVU6KmiDFNRFCVXVNIqxWo28OBco/d0n198uI2xRu2oaaedCpygSWDSX5O4/uh6Vl0oiqIUKypplXKt+k/llM2Lus/N/vmE2W7v4GDhAEB0QjTj/hhHSGzGLSQURVGKG5W0SjlhZETtUSsIQbv2mh2x2G56i4Uvfo2ViRUAoU9CGfvHWB7FPyrKUJUyws/PDymlup+l5IlKWmWAvUMlonp+R6I0BqBe0lWebP6Orzt+jYmR9qmHoKggJu6ZyJPE/O9hpCiKUlBU0ioj6nt05mTdtOmwnqHrsLx0i8/bf45Au77bufBzvLn/TRI1atUMRVGKJ5W0ypBWQz7krGUr3WeXQ2/R1LgO01tP15UdvnuYDw9/iEZmv5CmoihKUVBJqwwxMjam5qjVPEA7CaMcMTxaNZz+tfoyrsk4Xb0dN3Yw59gctdyToijFjkpaZUy5SlWJ8FlKktT+r3dLCuTU8slMbDqRAfUG6Oqtu7yOxWcWF1WYiqIomVJJqwxya+3NiTppi2W2efALp3evYkbrGXRzSVu37Ltz37H6Yslc6kVRlNJJJa0yqvWwWZy2+o/uc70j73HvxiU+a/cZ7aq105V/ceILfrv2W1GEqCiKkoFKWmWUMDLC9bUfuSccAbARcST+NJzkp0+Z7zWf5pWb6+rOOjKLP279UVShKoqi6KikVYbZl6/Ikz4reCq1G0W6am5y/rvRWBiZs6jzItwc3ADQSA3vHHiHg8EHizJcRVEUlbTKujovtONMo7Qp7y2jAji2aT52ZnYs7bIUFzsXAJI0SUzdN5XjIceLKFJFURSVtBS06xMeL+ej+9zswmdcOfEXFSwr8P1L3+Nk7QTA0+SnTNwzkbNhZ4sqVKWM8vf3RwiBn59fuvKbN28ihMDFxSXfY3h5eSGEYN++ffnu61kuLi4IIdK9LCwsqFGjBgMHDsx2Ty4/Pz+EEMyaNSvT41FRUbRr1w4hBA0aNCA4OBiABw8esGrVKgYPHkyjRo2wtbXF2tqaRo0a8fbbbxMSUnLXGlVJS0EYGdF47A9cN3YFwEwkU277GCIe3KWKdRV+8P6BypaVAYhLimP8H+MJfFg0244rSmGbNWtWtolDX97e3vj6+uLr64u3tzcAGzZswMvLi6+++irX/YWFhdGxY0cOHz5M8+bNOXjwIM7OzgC8+eab+Pn5sWHDBoyNjfHx8cHLy4sHDx4wb9483N3dOXnyZL6+T1FRSUsBtPtvWQz7icdod4t15CH3lg8lKTGB6rbV+d77+7SV4RO1K8OrLU2UolatWjUCAwPZs2dPvvtavXo1gYGBtGrVKufKefDee+/h7++Pv78/W7Zs4fr164wfP153LPUsSR937tyhffv2nD59mvbt27N3714qVqyoO+7g4MDs2bO5ffs2Z8+eZf369ezYsYMbN24wePBgIiIiGDhwIElJSQb/ngVNJS1Fp5prA4Je/Fr3udHTMxxfMRUAV3tXlnVdhq2ZLQCRTyMZvXs0QVFBRRKrogCYmpri5uZG7dq1891XjRo1cHNzw8rKygCR5czU1JQvv/wSW1tbEhIS+P333/Vqd+3aNdq3b8+VK1fw8fFh9+7d2NnZpauzcOFCZs6cSbVq1dKV29jYsHz5cmxtbblx4wZHjhwx2PcpLCppKem80GkgR6qP0X32vL+GkzuXA1DfoT5LuyzVbWnyMP4hY3aP4fbj20USq2JYsbGxzJ07l5YtW2JnZ4elpSXu7u7MmjWLmJiYdHWvXr2KnZ0dJiYmHDhwIENfly5dwtraGlNT03S/GJ+91BYUFMTw4cNxdHTEwsICd3d3vvzyy1z96z+ne1qxsbHMmzcPT09PypUrh6WlJa6urgwYMICdO3emq5vZPS0hBLNnzwZg9uzZ6e5L5fdyIYClpSX16tUDIDQ0NMf6586do3379ty6dYuBAweyZcsWLC0tczWmlZUV9evXB8jV2V1xoZKWkkFrv7npFtZt8M90gi7+A0CTSk1Y0mUJlibaH5QHcQ8Y/ftogqNL3l9+JU1wcDCtWrXi3Xff5datW3h6evLSSy8RGRnJ7Nmzadu2LZGRkbr69erV47vvviM5OZmhQ4cSHh6uO/bkyRMGDBjAkydP+Oyzz/D09MwwXlBQEB4eHuzduxcvLy86duzI9evXmTZtGgMGDECjyf+Czbdu3aJFixa8/fbbXLhwAU9PT3r37k3VqlXZtWsXc+fOzbEPX19fXnjhBQBeeOEF3T0pX19fmjZtmu8YQTuZAsDR0THbekePHsXLy4vQ0FBee+011q1bh6mpaa7HS0pK4ubNmwBUrVo11+2LmklRB6AUP0bGxriMXUfwwnY4y/tYiaeYbRxBVJWD2FdwpIVjCxZ3XsyEPycQnxxPSGwIY34fw0rvlVS1KXk/BMyyL+oI8m5WVL67kFIycOBALl26xKRJk/jf//6nu0QWFxfH2LFjWbNmDVOnTk23ceOQIUPYt28fy5Yt49VXX2XHjh0IIZg4cSKXLl2ie/fuTJs2LdMxV69eTf/+/VmzZg0WFhaA9rJXx44d2bx5M0uXLmXChAl5/k4ajYa+ffty5coVevfuzcqVKylfvrzueHR0NMeOHcuxH39/f2bNmsXZs2fp06ePQc6unnXx4kWCgoIwNTXlpZdeyrLeoUOHmDdvHrGxsUybNo0vvvgiz2P+8MMPhIeHU6VKFf7zn//k3KCYUWdaSqbsy1ckaeAankhzAKrJUG59P5TklEs3Lau0ZGGnhZgZmQFwN+Yuo3aPIiS25E6lLasCAgI4cuQIbdq0YcGCBenu6VhaWrJ06VIqV67M2rVr051tASxYsIAmTZqwa9cuvvjiC1avXo2/vz/Ozs6sXr0aIUSmY1paWrJkyRJdwgKoW7cun3zyCUCeZtM9a+vWrZw+fRoXFxfWrVuXLmEB2Nra0rlz53yNkR+RkZHs2rWLfv36odFoWLBggW7mX2b27NlDbGwsrVq1ylfCOn/+PG+//TYAc+fOxczMLM99FRWVtJQsuTTw4LJn2iWUJvEnOLZ8qu6zp5MnCzotwNRIe4kiOCZYJa4SKPXeTv/+/TEyyvgrwdraGg8PD5KSkjh+PP3D5RYWFmzYsAEbGxtmzJjB+PHjMTY2Zt26dVSoUCHLMV966SUqV66coXzYsGEYGRnx77//cvfu3Tx/p4CAAF1/ub3nU1A6duyoux/m4OBA9+7duXXrFrt27dLNIsyKp6cnpqamHDt2jMmTJ2dbNyvBwcH06tWLmJgYxowZw4gRI/LUT1FTlweVbDXv5seRO6fxvOsPgOf91ZzY1giPXtr9t9pVa8dXXl/xxr43SNIkcSf6DqN2j2KF9wqqWFcpusBzwwCX2EqyGzduAPD222/r/hWelbCwsAxl9erVY+7cuUyYMIGkpCRmz55Nu3btMmmdplatWpmWm5mZUbVqVe7evUtwcHCG2W/6unXrFgBubm55al8QvL29qVKlClJKQkJCOHDgAPHx8bz66qscPnyYOnXqZNn2pZdeYtq0aQwePJhFixYhhGDBggV6jx0SEkLnzp11EziWLl1qiK9UJFTSUnLUauSXnP0ykBfitJMxGp2YwbXqDanbtD0AHap34MsOX/LW/rd0iWtkwEhWdltZchJXGZacnAxAhw4dclxZombNmpm2//nnn3Wfjx8/jpQyy0uD+spv++Lmvffew8vLS/f5/v37eHt7c/78eYYNG8bRo0ez/c79+vXj559/ZvDgwSxcuBBAr8T14MEDOnXqxNWrV+nduzdr167F2Ng439+nqKikpeTI2MSEWuPWcWvhi9TUBGMhErHb7Et4lX1UrFIDgE41OjG/w3ze3P8mSZokgmOCGRkwkhXeK0rm5IwypHr16gAMGDCAiRMn5rr9rFmzOHDgAK1atSI2Npbt27czf/583nrrrSzbpM5ee15CQgL3798HwMnJKdexpEpNrleuXMlzHwWtatWqrF+/niZNmnDs2DHWrl3L8OHDs23zfOISQvD1119nWT8sLIxOnToRGBhIjx49WL9+PSYmJfvXvrqnpejFrlwFjIb8xGO0N+kdeUj48kE8jX+iq9OxRke+8voKEyPtD0VwTDAjd4/kfsz9IolZ0Y+Pj3bdyQ0bNuS67Z49e/jss88oV64cv/zyC7/88gtWVlZMnz4929l5v//+e6aXGtetW4dGo6F27drZTkzISeoySWvWrCE+Pj7P/QC6yQoFsXqEm5ubbpbkrFmz9BojNXGZmpqyYMEC3njjjUzrhYeH06lTJy5evIi3tzebNm0qkRMvnqeSlqK36nVf4GaHRSRL7SUMt8RLnF06GvnMMzVe1b3SJa67MXcZuXukeo6rGOvTpw8tWrRg//79/N///R8REREZ6ty4cYPFixenKwsNDWXYsGFoNBqWL1+Oi4sL7u7uLFy4kMTERAYNGsSjR48yHfPJkydMmjSJp0+f6squX7/Ohx9+CMCUKVPy9Z169+5N06ZNuXnzJsOGDdM9C5UqOjpa76WfUu+rBQYWzHqbM2bMwNbWluvXr/Pjjz/q1eb5xDV16tR0xyMiIujcuTMXLlyga9eubN68GXNz84IIv/BJKYvVCxgKHASigBjgBDARMMpFH16A1PNV47m2/jnUv6xvHC1atJCl0ZEfZ0r5kZ3udWTNrAx19t3eJ5utbiYb+TeSjfwbyc7rO8ubUTeLIFqtS5cuFdnYJcGdO3dk48aNJSBtbW1lu3bt5ODBg2WXLl1kvXr1JCAdHR119ZOTk2Xnzp0lICdOnJihv6FDh0pA9uvXL135Rx99JAE5YsQI6eDgIKtVqyYHDhwou3fvLi0sLCQge/XqJZOTk9O1W7lypQSkr69vuvKgoCAJyJo1a2aI4caNG7JOnTq67+Tj4yMHDx4s27ZtK62trWWHDh3S1e/QoYME5N69e9OV379/X1pZWUlAtm/fXvr5+cnRo0fLLVu25PwHK6WsWbNmpv0+a/bs2RKQrq6uMjExUVfu6+srAfnRRx9l2m7Tpk3S1NRUAvKNN97Qlfft21cCUgghBw0aJH19fTN9HTx4UK/vYAi5+RkETsisfr9ndaAoXsDilMQQB2wHfgMep5T9Chjr2Y9bSvLJ6nUppc9/AfFc29SkdSiLtp/r+31Ka9LSJCfLY/MH6JJW8kx7eWbPugz1Dtw5IJuvbq5LXB1/6SivR14vgohV0tJHXFyc/Oabb+SLL74oy5cvL01NTWWVKlVkixYt5LRp0+Thw4d1dT/++GMJyKZNm8r4+PgMfUVHR8u6detKQC5atEhXnpq0PvroI3n9+nU5aNAgWalSJWlmZibd3Nzk3LlzZUJCQob+8pK0pJTy8ePH8tNPP5XNmzeXNjY20tLSUtaqVUsOGjRIBgQEpKubVdKSUsq//vpLenl5SXt7eymEyDaRPE+fpBUdHS0dHR0lIH/44QddeU5JS8rME1fqd8nptXLlSr2+gyGUuqQF9E/5g7wP1H2m3PGZJDPFQGNdTOnv/UyOpSYtv/yOU1qTlpRSxsfFysD/ttElruiZjjLo0vEM9Y7cOyI9fvTQJa4Xf35RXn54udDjVUmreHg2aSlli6GSVnG6p5W6fe67UsprqYVSylAg9cm794QQ+YpZCOEJNASSgVX56assM7ewotKYDYRQCQAbEYfZ+qFEPEj/QGibqm34tsu3ukV2I+IjGP37aC4+vFjoMSuKUvIVi6QlhHAGWgAJQIYpTFLK/cBdoArQJp/DjUp5D5BS5v2Re4UKjs48eWWtbqknJxlK6PcD0s0oBPCo4sF3Xb/DxtQGgKinUby2+zXOPDhT6DErilKyFYukBTRLeb8opYzLos7x5+rmmro8ZX8AACAASURBVBDCChiU8nF5DtU7CiHmCyGWCSE+EUJ45/csrzRybdSaq+2+RpMyo7BB4kXOL3k13YxCgKaVm/LDSz9gZ6bd9yc6MZqxf4zlyL2St5+PoihFp7j8Ek5d0+VWNnVSN23KfP0X/QwAbIEHaCd6ZOdVYCrwGvABEACcF0I0zsf4pVLTrkM5VidtPTSPx39wdNX0DPXcK7qzwnuFbgfkuKQ4Ju6ZyF+3/yq0WJWiNWvWLKSUBl8tXSk7ikvSskl5j82mTuoudLb5GCf10uBqKWViFnXOAJMB95S4nICewFm098L+FELkbUG0Uqz1sFkcK99T99nz1lJO7vghQ736DvXx7+aPo5V276BETSJv7nuTHTd2FFqsiqKUXMUlaaUuuCULbAAh6gAvpnxckVU9KeXXUspFUspLUspYKeV9KeUOoBVwFKhM2qSRzMYZK4Q4IYQ4kdkT/6WVMDKi2fgVXDBP2xiv0bH3uHwi4wOctexrsdpnNTVstUtAJctkph+czoaruV+RQVGUsqW4JK3olHebbOqkHovOpk52Us+yjkgpc/1ou5QyAfg85WP3bOotk1J6SCk9KlWqlIcwSy5TM3Oqj9vIbSPtiai5SKTydj/u3sj4x+1k44R/N3/qlNOubC2RfHzkY5afX5766IGiKEoGxSVp3Ux5z7iEdJrqz9XVmxDCGO09Ksh5AkZ2Lqe8q8uDWbB3qITx8A1EplzFdeAxyWv6E/UwNEPdSlaVWOm9kkYVGunKvj71NV+e+FIlLkVRMlVcktbplHd3IURWO7a1fK5ubnijTTSxwC95aJ8qdVe7mGxrlXHVXN0J9VnBU6ndHLKG5i7B3/XPMBUeoJxFOb5/6XtaVWmlK1t1aRUfHv6QJI3hFyhVFKVky1PSEkJUFUL4CCEmCCHeTXlNEEJ0F0Lkeh8KKeUd4BRghnaG3/PjdQCcgRAgL3OkR6e8/yKlzE/CGZjyfjzbWgpurV/iQqs5us/uCec5v2REhqnwADZmNizpsoTONdK2P99yfQtT900lPil/K3QrilK65CppCSF6CSGOAsFop4wvAj5LeS0CtgHBQogjQoieWfeUqdT7Rf9LmTSROmZlYEnKxzlSSs0zxz4XQlwWQnxOFoQQFdHO/oMcLg0KIZoKIXqmXE58ttxECPEm2lmFAF/p9Y3KuBY9xnDE9dmp8H9ydEXmeyyZG5szr8M8+tXtpyvbd2cf4/8cT3RCXm9jKopS2uidtIQQXwOb0c6iS0I7BXwbsC7ltQ3tdPEkoDWwRQih9y93KeVG4Fu0q16cF0JsE0L8ClxDO9V8M/DNc82qAvVT3rMyAu0Z3GUp5d85hOGS8j0epCTeDUKIALTPj32ZUuddKeVufb9XWddm+Gz+cXhZ99kzeAXHNmW+aZ2JkQmzPGcxqtEoXdmJ0BOMDBhJ2JOyMxNTUZSs6ZW0hBCD0Z5lPAD8gHJSyuZSyj5SyuEprz5SyhZAOWBkSt3JQohBWXb8HCnlBGAY2kuFHdDei/oXmAT0l1Im6//VdEamvGc5zf0ZZ4EFwBWgBtArJY4nwEqglZRybh5iKLOEkREtxi/nnEVLXVnzc7M5uzfz6e1CCKa2mMpbLdLOyK5EXmH4zuEERQUVeLyKohRvQp9ZWkKIA4AH0FhKeV2vjoWoizYJHJdSdshXlCWUh4eHPHHiRFGHUSzEPI4kZEFn6iRr//o8kebc7buJuk3bZ9lmy79b+Ojvj0hO+beKvbk933T6hqaVm2bZJiuBgYE0aNAgb8EripJvufkZFEKclFJ6ZHZM38uDTYA9+iYsgJSV2v8EXtC3jVJ62diVp9zo33SrwluJp5TfPJx7QZezbNO7Tm8WdVqEpYl2QmnU0yhe+/019t3ZVxghK8WIi4sLQghu3rxZKOMJIRBC5FxRKXT6Ji1j4GmOtTJKTGmrKFR0qsnTweuJwlr7mUck/diPR+EhWbZp79w+3XqF8cnxTNk7hU1XNxVKzErp4+fnhxACf3//og5FyQN9k9ZloHPKTDy9CCEqAZ1JeyBXUajp1px7PivTPcMV8l1f4mKzniHYqGIjfvT5EWcbZwA0UsOsI7P45vQ36iFkpUAEBgYSGJjrhXOUQqBv0vIH7IG/hBBtc6oshGgH/IV2cduVeY5OKZUatPbmYpsvdNuZuCVe4sriASQlJmTZpoZdDX7s/iMNHNKuiX937jtmHJpBYnJWax8rSt64ubnh5uZW1GEomdA3aS1BOxW8EXBACHFLCLFRCPGlEOJjIcTslP/eKIS4BexHu0r6Ninlkuw6Vsqm5j4jOVZ/mu5z0ydHOLVkZKYPH6eqaFmRld1W0tYp7d9N225s4//+/D8eJzwu0HhLu3/++Ye3334bDw8PHB0dMTMzw8nJiVdeeYWjR49mqD9r1iyEEMyaNYvQ0FDGjRuHs7Mz5ubm1KpVi/fee4/4+IwPhkdHR7Ns2TL69OlDnTp1sLKywsbGhmbNmvHpp58SF5fVdnrpSSmpW7cuQohM40vVr18/hBAsWbKEmzdvIoRg1SrthuUjR47U3bt6/nJhdve0EhMTWbZsGR07dsTBwQFzc3Nq1KhBz549Wbt2rV7xK/kgpdTrhTbBvQWEAppnXskpr2fLQlLqGunbf2l8tWjRQirZ+3vpJCk/stO9jnw/Ncc2CckJ8qPDH8lG/o10r96/9ZbB0cFZtrl06ZIhwy51OnfuLI2NjWWTJk1kz549Zf/+/WWjRo0kII2NjeX69evT1f/oo48kIEeNGiWrVasmnZyc5CuvvCJfeuklaWVlJQHZq1evDOMcPHhQArJy5cqyffv2ctCgQbJLly7S1tZWArJVq1YyLi4uQ7uaNWtKQAYFBenKvvrqKwnIESNGZPqdgoODpYmJibS1tZWPHz+WYWFh0tfXV9auXVsCsm3bttLX11f3OnjwoK4t2h0nMvQZEREhPT09JSDNzc1lp06d5ODBg+WLL74oy5UrJ2vWrKnnn3jZk5ufQeCEzOL3ql5T3p+VsnuvJ9odhGuhXX1doF19/SbaZ6yOyGdWriir1JT3nEmNhhMLBtMyKu157X8azqD1wHeybyclyy8sZ8GpBbqyChYVWNx5Me4V3TPUz266beNVJXdfz/O+5w3ST0BAAM2aNcPR0TFd+bZt2+jfvz+2trbcuXMHKysrQHumNXv2bADGjBnD4sWLMTMzA7R/1q1atSImJoZDhw7Rtm3amXFwcDBXr17Fy8sLI6O0Cz2PHj1iyJAhBAQEMGfOHN599910cbi4uHDr1i2CgoJwcXEBICoqCmdnZxITEwkODqZixfS33GfOnMknn3zCxIkT+eabtHUJ/Pz8WLVqFStXrsTPzy/TP4/Us6znfz/27t2brVu34unpycaNG3FyctIdi4+PZ+/evfj4+GT+h1zGFfaUdx0ppUZKeVhK+Y2U8i0p5Tgp5diU/16UcqzMJyxFP8LIiKYTf+TsMw8ft7z4Gad3ZX8rVAjBmMZjmNN+DqZG2kkdD+MfMnL3SLUTch5069YtQ8IC6NWrFwMGDCAiIoK9e/dmOF69enUWLlyoS1gADRo0YMSIEQDs2ZN+PzVnZ2c6deqULmEBlCtXjoULFwKwceNGvWK2t7dnxIgRPH36lBUr0q8dkJiYyPfffw/AhAkT9OovJ2fOnGHr1q3Y2NiwZcuWdAkLwMLCQiWsQmBS1AEoiqmZOXUnbeLq112ol3QVIyFxPzqNC7YVaNTu5Wzb9nDtgaOVI1P2TuFxwmPikuJ4Y+8bTG0xFT93P/WsTS6Eh4ezfft2Lly4wKNHj0hK0q6yf+HCBQCuXr1Kjx490rXp1KkTlpYZN2ZIncRw7969DMeklBw+fJgDBw4QHBxMXFzcs7chuHr1qt4xT5o0iW+//ZalS5cybdo0XTL89ddfCQkJwcvLi4YNG+rdX3YCAgIA7dlWWdsrrzhRSUspFqxs7Kk4djN3vu1CdXkPM5FErT9e45pN+WxXzQDwqOLBj91/ZOKfEwmOCUYimX9yPkFRQXzY5kNMjU2zbW+oS2wl2Xfffcebb77JkycZt49J9fhxxskuNWrUyLSunZ0dQIbJGKGhofTr14+//856GdDMxslKw4YN6dKlC3/++ScBAQF0767dn3XJEu38r4kTJ+rdV05u3boFoGYVFrE876clhKgjhPheCPGvEOKJECI5i5faFEnRi0Plapj4beYB2geJrUU8FTYP5c61szm2dbV35aceP9G8cnNd2W///sa4P8cR9TSqwGIuDU6cOMH48eNJTEzkiy++4PLly8TExKDRaJBSMn36dCDj/R0gw2W+nIwZM4a///6btm3b8scff/DgwQMSEhKQUvL0aV7WL4DXX38dSEtUFy9e5MCBAzg5OdGnT5889akUX3ndT8sD7YSLUYArYIF2MkZmr+Ky0aRSAlStWZ8ngzboVs1w4DGma/vz4G7Oi+WWtyjP9y99z8u10y4pHg85zrCdw9SGktnYuHEjUkomT57MtGnTqF+/PtbW1rpLq//++69BxomNjWXnzp0YGxuzfft2unTpQqVKlTA1Nc3XOD179qRWrVrs2rWLmzdvsnjxYgDGjh2LiYnhLibVrKndWP3KlSsG61PJvbwmlLloZw2uB5oDtlJKo6xeBotWKRNcGnhwv+ePPJHmAFQhjCfLX852uadUZsZm/Lftf5nSfIqu7NbjW4THhROToDaczkxERASgnVTxvLCwMP744w+DjBMVFYVGo8HW1pZy5cplOJ7XZ5yMjIyYMGECGo2GL774gjVr1mBiYsLYsWMzrZ86aST1np2+vL29AdiyZQvh4eF5ilXJv7wmlNZAoJRyiJTyjJQy1pBBKYqbR2f+9VpCotQuXemiuc2Dpb2IeRyZY9vUmYXzveZjYWwBaJd+uvX4Fg/jHqqln56Teo9m9erVxMSkJfbo6GhGjRrFo0ePDDKOo6Mj5cuX59GjR/z000/pjgUEBDB//vw89z169GisrKxYsmQJ0dHR9O3bl6pVM99mr1q1agC5XqapWbNm9OrVS9f//fv30x2Pj49n165defsCit7ymrTi0G47oigFpknHVzjbco5uuad6SVe5ubgP8XH6/Rupa82u+Hfzp7JlZV1ZSGwI92LvoVFPZeiMHDmS6tWrc+rUKVxdXenXrx99+/bFxcWFEydOMGrUqJw70YOxsTEzZswAYNiwYfznP/9h6NChtG7dGh8fH9588808912+fHmGDx+u+5zdBIzevXtjZGTE119/jbe3N6NHj9bda8uJv78/LVu25NChQ7i6utK1a1eGDh2Kl5cXVatWZfz48Xn+Dop+8pq0jqG9l6UoBcqj51iOu7+v+9zo6RkCF2W/TuGz3Cu6s67nOsyM054jehT/iJuPb5KoUWsWgvYX/okTJxg7diw2Njbs2LGDEydO0K9fP06dOpXpZcO8euutt9i4cSNt2rTh4sWLbN++HWNjY9asWcOnn36ar767du0KgLu7Ox06ZL2FX9OmTfnll19o2bIlf//9NytWrGD58uV6TbV3cHDg4MGDLFq0iObNm3Ps2DF+/fVXgoKCaN++PXPmzMnXd1BylusVMQBSFs39CxgipfzV4FGVEmpFDMM56v8+bW4u1n0+bt+NFpN/wshYv51vLl26hF11u3QzCU2MTKhhV0O3X5dSsvXt25fNmzezZMkSdcZTDBXZihgAUsrDwGBgmRDiRyHEcCGElxDixcxeeRlDUZ7V+tX/crTKMN3nllEBHFs6LtsFdp8lhKCaTTUcrdNWfUjSJBEUFaSmxJcCJ0+eZOvWrVSoUIFXX321qMNRClB+5oOaAU+AoSmvrMh8jqMoCCMjWo/9hmOLomgVuR2ANmEbOLLCFs8xX+nXhxBUtKyIubE5wdHBaKT2OaTg6GDikuJwtHJUK2iUMGPGjCEmJoadO3ei0Wj4+OOPsba2LuqwlAKUp2QihOgPrEV7pvYQ7UK5aj6xUqCEkREtJq7i1Nf9aB6zHwDP4BUcXW1Dm1c/0bsfWzNbXO1duR19m4Rk7b2xh3EPiU+Kx9nWGRMj9W+skmL58uUYGRlRs2ZNZs6cabB1BpXiK68/ne+jfXB4ArBMLZCrFBZjExMavb6es1/34oW4YwC0ubGQf36xofWgd3NoncbcxBxXe1fuxtwlOkG7a3JsYiw3Ht2gul11dZ+rhFCPL5Q9eZ096AYcllIuVQlLKWxm5hbUf/03Lpo10ZW1DvyM41tyt9+osZEx1W2rU8kqbfHTRE0iQVFBRMbn/DyYoiiFL69JKwoINmQgipIbFlY21Jy0lSsm9XVlzU+9z6kA/1z1I4SgslVlatjVwEhofxyklNyLuce9GPU8l6IUN3lNWr8DLYW6a60UIRu78lQZv40bRi4AGAtJoyNvcvavn3Pdl62ZLa7lXDE3MdeVRcZHcjPqpu6+l6IoRS+vSWsGYAvME0Kou9ZKkbGv4Ijd2O3cEdoN+cxEMm77J3H+wJYMdXO6/2FubE4tu1rYmdvpyuKS4rgRdUN330tRlNwz5L3HvCac0cBO4A2grxDiL7SXCzO7liKllPpP7VKUXKpYpTqho7dzb7kPTjIUc5FInT1juGRmTsM23QDtEkKJiYnpdtjNjLGRMc42zkSYRBAaG4pEkqxJ5vbj21SyqkQly0pqWryi5FJiYiLGei4EkJO8roihQfv8VXY/vanHpZTSMNGWMGpFjMJ17+YVjP19cOQhALHSgrsv/0y9Fh25f/8+pqamVKxYUe/+niQ+4U70nXTbmlibWqtp8YqSS+Hh4SQmJma5iPHzslsRI68/eR+jTUqKUmw4udTnzvAthK/pSUUeYS3iqbJtGP+abKC6mwe3b98GtLvqmpqa5njGZGVqRe1ytQmODiY2UbtIb2xiLNcfXae6bXWsTK0K/DspSkklpSQxMZHHjx8TGRmZ5S7XuZWnMy1FP+pMq2jcDDyB/S99KY922/ZH2BA54Fec6jQlIiKC6OhokpOT9e5PSkl0YnS6/bgEAlszW6zNrBHZXnBQlLLL2NgYW1tbHBwcMDc3z7lBiuzOtPRKWkKI4cB2KaVhNtYpI1TSKjrXz/1NpV9fwQ7tGVIEdkQP2kzNBi3y3OeB4ANMPzidxwmPdWVtq7Xls3af4WDhkO+YFUXRMsSCuauBUCFEgBBinBCiiuHCUxTDq93kP4S8vI5oqV3ZwoHHWP/SjzvX8r4N3IvOL7Kh1wZeqPSCruzw3cMM2DqA4yHH8x2zoig50zdpvYV2D60uwLdAsBDioBDiTSFErQKLTlHyoV7zDtztuYZYqd29uCKPMF/bh7s3Lua5TycbJ1Z2W8nIRiN1ZQ/iHjDm9zF8e/ZbkjX6X3ZUFCX3cnVPSwhRGegL9AO8AFO0EzLOA5uA36SUFwwfZsmkLg8WD5eO7MIlwBcr8RSAECqh8duBk0v9HFpm72DwQWYcmkHk07Qln1pXac1n7T+jslXlbFoqipIdg+2nJaV8IKX8TkrpDVQGXgW2AHWA2cBZIcQ1IcQcIUSb/AauKIbQ0NOHG12XEy9NAahCGKzqScjtnHeqzU575/Zs6LWBFo5p98n+CfmH/lv7s/f23nz1rShK5gwye1AIYQn4oD0D6w6UQ3sGdh/4DdgM7C1ri+uqM63i5fz+X6n311jMRSIA94QjxiN34Fijbr76TdIksfTsUpadW4Z85kmQQfUHMc1jGhYmFvnqX1HKGoPvXPw8KWWclPJXKeVwtGdgPsAPgDEwEe1ahTMMMZai5FXjDv247LWUBKl9PNFJhpK0sgehwdfz1a+JkQmTmk1iuffydJcFf7nyC4O3D+ZKxJV89a8oShq9kpYQYr4QYrA+daWUSVLK3VLKcYAT0AFYgPasS1GK1AsdXyGwwxJd4qomQ0lc3oOwu0H57rtllZZs6rWJLjW66MquR11n6I6hrA1cq/Z+UhQD0Pc5LQ3gL6UclfI5OeXz6AKOr0RTlweLrzN//kzDgxMwE9rZfneEExav7aKSk0u++5ZSsunaJv537H/EJ8fryl90fpHZ/5lNRUv9l5JSlLLIEJcHk9HOFNT1SfbrDipKsda0y2AutltEQsqymNXlPeK/9yH83s189y2E4JV6r/BLr19wc3DTlR8IPkC/Lf3Yc2tPvsdQlLJK36T1AGiq9s9SSpNmXYdx4T8LSXwmccV970OYARIXgKu9K2u7r+XVhq/qyiKfRvLGvjf44NAH6ZaFUhRFP/peHlwLDAFuAUFon9EKAS7rMYaUUnbOR4wllro8WDKc2r2axn+/gWkBXCpMdeTeET44/AEPnjzQlTlZO/Fpu0/xqJLpVRBFKbMMsfagM/ArkJefLrU1iVLsZZq4xuykUjXDLfgS9TSKz/75jJ1BO3VlAoGvuy+vN3sdM+Ps9/pSlLIi30nrmY5cgBrAPiAA+J8+7aSU+/UepBRRSatkObX7Rxr/PSVd4jIfs5PKBkxcAAFBAXxy9JN0C+/WKVeHz9t/nu4emKKUVQZLWs90mG42oZI5lbRKnlO719D478m6xBUsqmI6egeOzrUNOk5obCgz/57J3/f+1pWZCBPGNhnLmCZjMDUyzaa1opRuBfFwcS3g7byHpCjFU3Pv4Zx/ZnKGs7xP0nIfQm9fM+g4jtaOLO2ylPdbv4+FsXbFjCSZxJKzSxi6Y6h6IFlRspCnpCWlvCWlfGjoYBSlOGjuPZwLz0yHryZDSV7Zg/u3DJtIhBAMcRvChl4baFqpqa78csRlBm8fzLdnviVRk2jQMRWlpMvX2oNCiOpoV7xwArJaYE1KKT/J8yAlmLo8WLJpH0CeiJlIAuC+qIR8dTtOtQx/3ylZk8yawDUsOr2Ip8lPdeUNHBrwSdtPqO+QvxXpFaUkKYh7WibAN8AY0h4yfv4ZLplSlqvZg0KIocB4oAnatQsvAyuBb3Oz4K4Qwh/wzabKFSlllr99DBGHSlol39m/1uO2f4Jukd0QKpH86jaquTYokPGCooL48PCHnA1L26zSxCjlXlejMZgaq3tdSulXEEnrv8D7QBKwE7gGZPmkpJRytp79LgYmAPHAHiAR6AzYol0tfoCUUq9d9p5JWoeBfzOpcl9KOb0g41BJq3Q4u28TbnvH6RLXAxx4Onwr1es0LpDxkjXJ/HjpRxadXkSCJkFXXrd8XWZ7zqZxpYIZV1GKi4JIWrcAB6CtlPJcPuNL7bM/sBHtQ8svSimvpZQ7AnuBBsAbUsoFevbnjzZpjZRS+hdFHCpplR7nDmym3p4xWKQkrjDK82TIZmrWb5pDy7y7EXWDDw9/yLmwtB8xI2HEsAbDmNR0ElamVgU2tqIUpYKYPVgZ2G+ohJUi9azn3dREASClDEV7mQ7gPSGEQbZTKQFxKMVIkxf78G/XFTyR5gBUIhKbdS9z89LxAhvT1d6V1d1Wa/fkSplhqJEafrz0I/229uPvu3/n0IOilD55/cV7G3iaYy09pay40QJIADY8fzzl4eS7QBWgwHZELi5xKMVTo3YvE9RtlS5xVSAK+/X9uHHhaIGNaWxkjK+7L7/2/pU2VdP+yt2Nucu4P8cx49AMHsU/KrDxFaW4yWvS+hnoIISwMVAczVLeL0op47Koc/y5uvrqmLIf2DIhxCdCCO9szpIKMg6lFHD39OFWj7XESEsAyvOYChv78e+ZgwU6bnXb6izruoxP2n6CnZmdrnzr9a303tKbXUG71H5dSpmQ16T1GXAF2CGEqGeAOFLXybmVTZ3bz9XV16vAVOA14AO0y0+dF0Jkdje7IONQSokGrboS3OsnHkvtPSV7Yqm8eRDXTu4t0HGFEPSp04ctfbbg7eKtK4+Ij+CdA+8w/s/x3Hl8p0BjUJSilteHi58CLwHlgYtCiH+FEPuEEH9l8tJn86DUM7bYbOqkzk601TPMM8BkwD2lfyegJ3AWaAj8KYSoZug4hBBjhRAnhBAnwsLC9AxVKWncPDoR2nc9UVgDYEcsVbcO4fKx3QU+dkXLiszrMI+FHRdS2bKyrvzwvcP03dqX785+R0JyQjY9KErJlaekJYSoCBxCmxCMAVfgRbRblmT2yrHLlHeDXd+QUn4tpVwkpbwkpYyVUt6XUu4AWgFH0U4meX7Ke77jkFIuk1J6SCk9KlWqlNdulBKgbtP2hPXfRGTKv19sRBw1dozg0uFthTJ+xxod2dxnM0PchiBS/uo+TX7KN2e+of/W/hy7f6xQ4lCUwpTXy4NzgBeAq2gvvfUCOmbx6qRHf9Ep79ndI0s9Fp1NnRxJKROAz1M+di+qOJTSoU5jTx4N3Ew45QCwEk9x/X0kF/b/Wijj25rZ8n7r91nXYx0NHNIeeL75+Cajfx/N+wff52GcWnFNKT1M8tiuB3AfaCOljDJAHDdT3mtmU6f6c3XzI3XzyucvDxZ2HEopUKuhB7eGbEGzri+VicBCJFLvr9c4m/SUFzoPKZQY3Cu6s67HOn6+8jOLTi8iNlF7hXvbjW3sD97PGy3eoH/d/hipJzWUEi6vf4Ntgb8NlLAATqe8uwshLLOo0/K5uvlRIeX9+VU8CjsOpZSoWb8pCSO2E0JFAMxEEg0PTOR0wKpCi8HYyJhhDYaxtc/WdBM1Hic85uMjHzNsxzDOh50vtHgUpSDkNWkFov+EiBxJKe8ApwAzYMDzx4UQHQBntKtUHDHAkANT3tM9GVoEcSiliHNtdzR+O7krHAEwFck0PvIGJ7cvK9Q4KltVZl6HeXzb5VucbZx15RceXmDozqHMPDxTXTJUSqy8Jq3FgJeBprunSr3P9D8hRJ3UQiFEZWBJysc5zy5WK4T4XAhxWQjx+TP9IIRoKoToKYQwfq7cRAjxJtpZhQBfGSIORUnl5FIfkzEB3BZOAJgIDc2Ov8OJ3xYWeiztqrXjt96/Ma7JOMyMzHTlv/37G71+68XawLUkaZIKPS5FyY88b00ihJiD9hmoD4HdUsrgZhjYKAAAH2lJREFUfAcjxBK0SyXFA3+StlCtHbAZeOXZhWqfWV9wlZTS75nyPmgXto1AO1kkGO2ZYWO0U981wHQp5VxDxJEVtfZg2RV+/zbR3/eglua2ruxYw/dpNfDdIonnTvQd5h6fy747+9KV1y1fl+mtptOySsvMGypKESiIBXP1Wmk9hZRS6j3hI2VLkIloE0zqliAryGRLkGySVi1gCtrp7TXR3sOSaJPXQWCxlPKkoeLIikpaZVvEg3tEfNeTOsnXdWX/1H2L1sNmFllMB4MP8r/j/+PW4/TPz/u4+PCmx5tUsa5SRJEpSpqCSFq5ujQmpSyTU5ZU0lKiIsIIWdKD+klpux4fdRlPG785RRZTQnICqy+tZtm5ZcQlpa1WZmFswahGo/B191UryCtFyuCrvEspjXLzyl/4ilJy2TtUwun1AC6ZNtKVtbn5LUe/n4LUFM1tUTNjM8Y0HsPWPlvxcfHRlccnx7Pk7BJ6be7Ftuvb0KjbtkoxlOd7WkrO1JmWkupJTBTXF/Wm8dO0JyX+cRxEq3FLEUZF+++64yHHmXt8LpcjLqcrd6/gzrut3qVZZbU2tFK4DH55UNGPSlrKs+LjYrm8sC9N4/7RlR2r0BuPCSsxMjbOpmXBS9Yks+X6FhaeWsjD+PTT4b1dvJnaYirVbJ5/Fl9RCkZBbAKpKEouWVha0/CNrZy0flFX1urhFk4tHEJyUmIRRqZ9MLlf3X7s6LeD1xq/lm6K/O6bu3n5t5f5+uTXxCQ8/zy+ohQulbQUpRCZmVvwwhubOG7XVVfmEbWbs1/3JzEhvggj07I2tWZy88ls7buVbi7ddOUJmgSWX1hO91+7szZwLYnJRZtklbJLXR4sQOryoJKV5ORkTnzjS+vItBXhz1h50mDyJswtrIswsvTOPDjD3ONzOR+efvknZxtnJjefjLeLt1rPUDE4dU+riKikpWRHk6zh2NKxtAnboCs7b/H/7d15dFRV2u/x70MSpgBKZBLCLIMBREAFghEZFKSNLWlARRFR2wmntpf3XfZ979t29+3VXrv7VRTRlxalVRQHgooiKKiIMskMQpBBCFOQWQghQ9Vz/zgnoVIZyHyqUs9nrbOOtfeuql07MT/OqX326cslj3xMg9gqWyWt0vzqZ/5P83lx7YscyDxQqC7hogSe6PcE/S/u71HvTG1koeURCy1zPur3s/zVx0k8cG5h3a11e9L24Xk0ahLnYc+KyvZlMzttNv/a9C9OZhdeK3tQm0H8ru/v6BbXzaPemdrEQssjFlqmLFSV5TOfInHPywVlP0Z3pdVDn9IkrkUpz/TGLzm/MGPTDGZtnUW2L7ugXBBu7HQjD/d5mNaNWnvYQxPuLLQ8YqFlymP5W39i4I7/Lni8K6ojTe//lKYtQnOqeUZmBtPWT+OjnR8VuhA5pk4MY7qO4be9fkvzhnb3blN+FloesdAy5bXy3Wfpv/WvBY/T68TT8N5Pada6g3edOo/tx7czZe0UluxbUqi8flR9but+G5N6TqJp/aYe9c6EIwstj1homYpYNXcq/db/J1Hi/L+5X1oRNWkerdpV5Z2Aqt7qjNU8v/Z5NhzeUKg8NiaWCQkTuDPhThrXDZ0JJiZ0WWh5xELLVNTqT2fQe9WTxLg3VMigGb4JH9Gmc8/zPNNbqsrS/UuZum4qW49tLVTXpG4TJvWcxPju421BXlMqCy2PWGiZylj3xdv0+PYR6opzo8bDNOXMram0797X456dn1/9LE5fzEvrXmLnyZ2F6uLqx3Fvr3sZ120c9aLqedRDE8ostDxioWUqa+OSuXT58j4aSA4Ax2nC8ZR36XRZosc9Kxuf38dnuz9j2vpp7D21t1Bd8wbNmdRzEmO6jqFBdAOPemhCkYWWRyy0TFXYvOwzOi68i1hxlnn6hVgOJc+iS78hHves7HL9uXy842Ne2fgKGZkZheri6scxqcckxnUbZ6cNDWCh5RkLLVNVtq7+kjaf3EETMgHI1Pqk3zCTSwfccJ5nhpYcXw7v//g+MzbN4HDW4UJ1Tes15c4ed3Jb99uIjQmdpaxMzbPQ8oiFlqlKOzYu46LUcTTlFABZWpddw1+lR9KvPe5Z+WX7skndnsqMTTM4dOZQobomdZswIWEC4y8dT5O6TTzqofGShZZHLLRMVftpyxoav5dCM04AkK0xbLtmKpcNu9XjnlVMji+Hj3Z+xKsbXy2yrmHjmMbcnnA7d1x6BxfUu8CjHhovWGh5xELLVIf07RupO2s0rTgCQK5GsXngP+kzcpLHPau4XH8un+z8hOkbp7Pv9L5CdQ2jGzK261gmJEygZWxLj3poapKFlkcstEx1ObB7G/rvZNqoc2rNp8L6K56hX/IDHvescvL8ecz/aT7TN05nzy97CtVF14kmuVMyd/W8i04XdPKoh6YmWGh5xELLVKdD+3dx9tVk2qtzZOJXYU3vp7ky5XGPe1Z5Pr+PhbsXMn3j9CLXeQnC0HZDubvn3VzW/DKPemiqk4WWRyy0THU7nJHOqek30sl/7qhk1aV/4Kpb/sPDXlUdv/r5Zt83zNg0g/WH1xepv6LlFdzT6x4GtR6EiHjQQ1MdLLQ8YqFlasLxIxkceXkUXXznjkhWdnmC/rf/0cNeVb21h9by2ubXiizMC9CtaTfu7nk313e4nug60R70zlQlCy2PWGiZmnLy2BEOTvsV3fPSCspWdHiIAXf9zcNeVY/tx7fz+ubXmf/TfHzqK1TXOrY14y8dT0qXFFucN4xZaHnEQsvUpFMnj7F3ajIJuZsLyla0vYf+k/6B1KnjYc+qx4HTB3hjyxukbk8lKy+rUF1sTCyjLxnN+EvH07ZxW496aCrKQssjFlqmpp3J/IWdL9xEr+x1BWUrW0/gqntfqJXBBXDi7AneSXuHt9Pe5kT2iUJ1+ZM2JiRMoG+Lvva9V5iw0PKIhZbxwtmsTLZN+TW9z35fULaq5TiuvP9/am1wAZzNO8u8XfN4a8tb7Dq5q0h9wkUJTEiYwIj2I4iJivGgh6asLLQ8YqFlvHI26wxbX/gNfbKWFZR932w0Vzw0A6kT5WHPqp+qsuzAMt7c8ibfHfiuSH2LBi24tfutjO06lgvrX+hBD835WGh5xELLeCknO5tNL4yhX+Y3BWXfx91Iv8lvUCeqdgdXvp0ndvLmljf5ZNcnZPuyC9XVj6rPqE6juKXbLSRclOBRD01xLLQ8YqFlvJabm8P6F27jylOLCspWX3gDfR9+izrRkTM1/NjZY7y/7X1mb5vNkawjReova34Zt3a7les7XG83pgwBFloesdAyocCXl8eaF2/nqpMLCsrWXHA9lz/yDlERFFzgLNC7YPcC3tzyJmnH0orUN63XlJQuKYztNpY2jdp40EMDFlqesdAyocLn8/H9i3cy4MQnBWVrmgzj8kffJSo68iYlqCobDm/gnbR3+HzP5+T58wrVC8Lg+MHc0v0WElsnUkdq7wSWUGSh5RELLRNK/D4fq16axIBjHxWUrWt8Lb0eeY/oupF7Suxo1lHm7pjLe9ve42DmwSL1bRu35ZZut3DzJTfbLVJqiIWWRyy0TKjx+/ysnHYvA4/OKShbF5tEz8fmEBPBwQXOIr3f7PuG2dtms+zAsiL19aLqMaLDCFK6pNg1X9XMQssjFlomFKnfz4qX72Pg4fcLytbHJtHDgqvAnl/28O62d/lwx4ecyjlVpL5Dkw6kdEkhuXMyzRo086CHtZuFlkcstEyoUr+f5f/zEImH3ikoWx97tRtc9T3sWWjJysvis58+Y3babLYe21qkPlqiubbttaR0SSGxdSJRtfwauJpioeURCy0TypzgmkziobcLytbHDqLHY6kWXEFUlS1HtzBn+xzm/zSfzNzMIm1aNmzJ6C6jufmSm23mYSVZaHnEQsuEOguu8juTe4bP93xO6vZU1v28rki9IAxsPZCULikMaTuEulF1PehleLPQ8oiFlgkHTnA9TOKhWQVlFlxls+vELlK3p/Lxzo85nn28SH3Tek0Z1WkUyZ2TSYhLsMkbZWSh5RELLRMuiguudbFJ9HpsTkRPhy+rXF8uX+39itTtqSw7sAyl6N/VThd0IrlzMjd2upFWsa086GX4sNDyiIWWCSfFBdfaxkPo/dj7EXkBckUdOH2AD3d8yNwdc8nIzChSLwhXtbqK5M7JDG8/nNiYWA96GdostDxioWXCjfr9LH/lQRJ/nl1QtrbJMHo/+l7ELflUWT6/j1UZq5i3cx6L0hcVuVElOIv2Dms/jOROyQy4eIDNPnRZaHnEQsuEI/X7WTHtPgYeOXcd1+oLRtD3kbcjapHdqnQm9wyL0xfz8c6PWXlwZbGnD5s3aM6ojs73X93iunnQy9BhoeURCy0Trvw+P6um3c2Ao3MLylY3HeWsDh8htzWpLhmZGcz/aT7zds5jx4kdxbbp2rQrozqOYkSHEcQ3jq/hHnrPQssjFlomnPl9PlZNncSA4+fWKlzVbDRXPvRarb4Dck1RVdKOpTFv1zw+3fUpx84eK7Zdr2a9GNFhBCM6jIiYCRxhFVoiMh54ELgMiALSgNeBl1XVX8bXiAGuAUYBg4D2wEXAYWA5MFVVvy7huTOBiaW8/DZV7V6WflhomXDn9/n4/sUJ9D/xaUHZqpa3cuX9L1twVaE8fx7LDyxn3s55fLn3yyI3rMzXp0WfggCrzctHhU1oichLwEPAWWAxkAsMAxoDc4Gxquorw+sMB75wH2YAa4BMIAHo6Zb/RVX/q5jnzsQJre+A4o7dD6rqU2X5PBZapjbw5eWxZsotXBVwI8mV8XfT/97nPOxV7XU65zRf7f2KBbsXsGz/MvI0r0gbQbii1RWM7DCS69pfR9P6TT3oafUJi9ASkd8AH+CEzDWqut0tbwl8BVwKPK6qU8rwWkNxwm+Kqi4NqrsFmIVzFDdUVb8Kqp+JE1qTVHVmZT6ThZapLXJzc9j4fAr9Ms/977Sy40P0n/g3D3tV+53MPsni9MUs+GkBqzJW4Svm3+xREkX/i/szssNIhrYbWitunxIuobUa6AdMVNU3guoGA1/jBFqbsp4mLOW9XgXuAV5T1XuC6mZioWVMEdnZWWx97tdcfnZlQdnKrr+n//giJyxMNTh29hiL9ixiwe4FrM5YXewMxOg60SS2TmR4u+EMbjuYuPpxHvS08kI+tEQkHtgL5AAXqmqRCxpEZB/QBhikqkVvdlO+95sMTAU+V9URQXUzsdAyplhnszLZ/tyv6JVzbs297y/7E1emPO5hryLP4TOH+XzP5yzcvbDY9Q8B6kgd+rTow7B2wxjabmhYLeIbDqGVDHwMrFPVviW0mQvcDDysqi9V8v2eBx4D/q2qdwXVzcQJrTeAo0Aj4BDwLfBFeY7yLLRMbZR56iR7XriBhNwfAPCpsHHgc/QZOcnjnkWmjMwMFu5eyMLdC9l0ZFOJ7brHdWdo26EMbTeUrk27hvQ6iKWFVqhcKdjR3e8ppU16UNsKEZFWwF3uwzmlNL2zmLItInKrqpb8m2FMLRfb+ALaTJ7H9hevo4tvJ1Gi9Fj+ezY3aEzPwWO87l7EaRXbiok9JjKxx0T2ndrH4vTFfJn+Jet+XlfoFGLasTTSjqUxbcM04hvFM7TdUIa1G0bv5r3DaiWOUDnS+gPwV2CWqt5RQpu/An8Apqvq/RV8n2hgAc6MxMWqOryYNo8DPpzZi3uAJkBft3+9gZ+Bvqq6v4T3uA+4D6Bdu3b99uwpLYeNCV9Hf97P6Zevp73uAyBL67LnV7PoftX1HvfMABzJOsKSvUtYnL6YFQdXkOvPLbZdXP04hrQdwtB2Q+l/cX/qRXm/QHI4nB7838D/Bd5S1QkltKmK0MqfgLEXuEpVi65mWfJz6wJLgAHAS6r68PmeY6cHTW13cO9OmDGCizkMwCkacDgllU6XJXrcMxPodM5pvj3wLV/u+ZJv9n9T7E0sARpGN+TqNlczuO1gBrUexEUNLqrhnjrC4fTgKXffqJQ2+XWnSmlTIhGZghNYGcCw8gQWgKrmiMjfgI9wLlo2JuJd3LYz6bencmRWMs04QWOyyEm9jQONPqN1pwSvu2dcjeo2YmSHkYzsMJIcXw6rMlaxOH0xX6V/xdGzRwvanclzbnD5+Z7PEYSezXqS1CaJpPgkEi5KoI54f0F5qBxp3YQTBqVNxEgFRgOPqOrUcr7+P4EncFbEuFZVt1Swn12BbUCOqp73GNqOtEyk2LFpBS3mpNAE51/w+6UV9R9YxEUt23rcM1Man9/HpiObWJy+mMXpi9l7am+JbePqx3F1m6tJik8isXUiTeo2qbZ+hcPpwbY4Ey1Km/K+F4gHrlbV78rx2s8CT+LMBBymqhsq0c+BwDLgmKqe97jZQstEkh+WL6DzgjuoL853JzujOtPysUU0ahKe1wpFGlVl+4ntfLPvG5buW8r6w+vxlzBZOkqi6N28N0nxSVwTfw1dLuxSpbMRQz60wOkkzoSHKru4WESeAf4DOI4TWMVf0FD2Pj4HPA4sVNWR52tvoWUizdovZtH728lEifN3ZXO9PnT93WfUrd/A456Z8jqZfZLlB5azdP9Svt3/bYkL+gK0bNiSpPgkktokMeDiATSMaVip9w6X0BoDvI8TTEmqusMtb4GzjFMCQcs4ud8xjQbmBq8HKCJ/Af4TOAEMV9U1ZejD5ThHc58FrnHozjp8FPg7UAcYqaoLz/d6FlomEq344HkGbP5jweO1jYdw+eNz7JYmYcyvfrYc3VJwFLb56OYS28bUieGKlleQFJ/EkLZDKnRrlbAILQARmYazwvtZYBHnFsxtAnwIjAkKk5k4FwIXukg44DsygNXADyW8ZZqqPhPwvJtxFuY9BvwI7MNZrLcX0BrwA0+p6rNl+TwWWiZSLXv9KRL3TCt4vKr5GK588F+2MnwtcSTrCMsOLGPpvqV8d+A7TuUUPz/u8b6Pc0+ve4qtK004zB4EQFUfEpFvgcnAYM7dmuQ1ynFrEiDwJPoV7lacJcAzAY83AFOAq3BuZ9IHUJzweh1nqvt5j9iMiXQDJ/6V5dN+ZuCRDwDQRi097pGpSs0aNOOmzjdxU+ebyPPnsfHwRucobP9Sfjz+Y0G7a+KvqfL3DqkjrdrGjrRMJPP5fKx9fhx0TLK1CSNIRmYG3+7/lo2HN/KnxD9VaIJG2JwerG0stEykU7/fTgmacisttOy3yRhTbSywTFWz3yhjjDFhw0LLGGNM2LDQMsYYEzYstIwxxoQNCy1jjDFhw0LLGGNM2LDrtKqRiBzGuftxRTQDjlRhdyKRjWHl2PhVno1hxbRX1ebFVVhohSgRWV3SxXWmbGwMK8fGr/JsDKuenR40xhgTNiy0jDHGhA0LrdA13esO1AI2hpVj41d5NoZVzL7TMsYYEzbsSMsYY0zYsNAyxhgTNiy0QoyIjBeRpSJyUkROi8hqEZksIhHxsxKRbiLymIi8JSJpIuIXERWRMWV4boXGrjaNuYjEiMgwEfmniKwQkYMikiMi+0XkAxG59jzPtzEUeURE3hORrSJyVERyReSwiCwSkTuklLsa2vjVAFW1LUQ24CVAgSzgE2Au8ItblgpEed3HGhiD593PG7yNqY6xq21jDgwPGLOD7md6F9gUUP5nG8NSx3AfkAOsBeYBs4HlgN/9TB8CdWz8PPr5eN0B29wfBPwm4A9Nl4DylsAWt+4xr/tZA+NwL/AsMA7oDHx9vtCq6NjVxjEHhgIfAEnF1N0C5Lmfa4iNYYljeDUQW0x5DyDD/UyTbPw8+vl43QHb3B8ErHZ/Qe8spm5wwC92kX/h1eatjKFVobGLxDEHXnU/1wwbwwqN3/9xP9PbNn7ebHa+NASISDzQD+eUxPvB9aq6BNgPtAIG1GzvQltFxy6Cx3ydu4/PL7AxLJc8d382v8DGr2ZZaIWGPu7+B1XNKqHN90FtjaOiYxepY97F3R8MKLMxLAMR6Qg84D6cF1Bl41eDor3ugAGgo7svbUX49KC2xlHRsYu4MReRVsBd7sM5AVU2hsUQkUk4p+hicI5ME3H+of83VZ0b0NTGrwZZaIWGRu4+s5Q2p91942ruS7ip6NhF1JiLSDTwFnABsFhVA48UbAyLNwiYGPA4D+c7rf8OamfjV4Ps9GBoyL/uw9bUKr+Kjl2kjfkrwDBgL3BHUJ2NYTFU9V5VFaAhzszB54GngRUi0jqgqY1fDbLQCg2n3H2jUtrk150qpU0kqujYRcyYi8gU4B6c6drDVDUjqImNYSlUNUtVt6jqk8BTQG9gakATG78aZKEVGna7+/altGkb1NY4drv78o5dRZ8XVkTkn8CjwGGcwNpeTLPd7t7G8Pxed/fJIhLj/vdud2/jVwMstEJD/jTkHiLSoIQ2Vwa1NY6Kjl2tH3MReRZ4AjgKXKeqW0poamNYdidwvtuKBuLcMhu/GmShFQJUdS/OkjF1gbHB9SIyGGf2UgbOcjLGVdGxq+1jLiLPAE8Cx3ECa0NJbW0My+UanMA6ARwBG78a5/XVzbY5GzCGc1e/XxJQ3gL4gQhdzoWyrYhRobGrrWMO/MXt+3GgXxmfY2Po9DsJuB2oV0zdIGCn+5n+YePnzWY3gQwhIjINeBDnavtFQC7OjK8mOIt0jlFVn3c9rH4i0heYFlCUgDPddztwLL9QVQcEPa9CY1fbxlxEbgI+ch+uxvnDV5w0VX0m6LkRP4YichfO91YncI6CMnB+/zrj/C4CfAqM1aALgm38aojXqWlb4Q0YD3yHs8pzJrAGmEyErD0GXEvxq7wX2qpy7GrTmONcPHze8QO+tjEs9nN0BP4MfIVzeUAWTpjsxlmI+ObqGIfaMn41sdmRljHGmLBhEzGMMcaEDQstY4wxYcNCyxhjTNiw0DLGGBM2LLSMMcaEDQstY4wxYcNCyxhjTNiw0DImTInIXSKiIjLTo/d/2n3/p714fxOZLLSMMVVKRK51w+xrr/tiap9orztgjAlbU4HZuKudG1MTLLSMMRWiqkewwDI1zE4PGlOFRKS7e2rs54A72wa3iRKRDLddj4DyWBH5XyLyvYj8IiJZIvKD+91RabdkL6kviSIyx32vHHf/gYgMKOU5IiLjROQz9zPkiMh+EVksIg8HtS3ynZZ7SvAr9+Fgtz5/+9p9/e3u49L6keq2eai8n9vUbhZaxlQhVU0DVgLNgVElNBsJtARWq+oPACISD6wC/h/O7deXA58DTYE/At+JSNOy9kNEHgSWAilAOs4K5enAb9zX+m0xz6mLcyuMd4HrgB/d56UBPYEXy/DWC4CF7n8fAv4dsC1QZ4Xul9z6YgNJRNoAycAp4M0yvKeJJF4vM2+bbbVtAx7Auf1Hagn177n1k93HAixzy14EGga0bYDzh1uBmUGvc1cJ5b1x7snkw7nvU2DdrW55DtAzqG6K+3rbgO5BdVHATUFlT7vtnw4qv5bSb39yAU4gnQWaFVP/Z/f5U73+WdoWepsdaRlT9Wbj/EG+UUSaBVa4R0s34YTGO27xSGAgsALnLrVn8turc6PBB4CfgdvLeLT1KM731bNV9f3AClWdjXP0FAM8FtCvFjg3IvQDKeocMQY+z6eqH5fhvc9LVU/iBHE94O7AOveUav5R4DSMCWKhZUwVU9UTOKfZYnBu7hfoVpw/1h+rav6dmPNPI85RVX8xr5eJcxfiaODKMnRhsLufWUL9a+7+2oCyoW5/l6t7yrKaTXX3D4hI4N+hFKAVzlHalhrohwkzFlrGVI+Z7n5iUPnEoHqATu7+70ETFwo2zgVb8zK8dxt3/1MJ9TuD2oHzPRo4319VOzeQFuHcKXhkQFX+91wvFXmSMdiUd2OqyxfAPqCviPRS1U0i0g3oD2TgTFjIF+Xul+Dc1r00e8rRh5JuSy7leI3q9CIwHCeo5rszKa8BDuAcqRpThIWWMdVAVf0i8ibwFM6Eid+7e4C3VNUX0Hyvu39fVaviCGM/0BnnCG5nMfUdA9rlyw/DblXw/mX1Cc7R4A0i0gGY7JZPV9W8GuyHCSN2etCY6jPT3d/uTie/I6g832fufmwVve8Sd39nCfWT3P3XAWVf4sw4TBSRSyv5/jnuvtR/FLvf303D+Tv0JM745AHTK/n+phaz0DKmmqjqjzhT2VsCfwfiCbg2K8CHwBqci3FfEZG44NcSkU4iMjm4vAQv4Pzxv01ERge9zlhgHE5AvRDQ15+BV3D+JswRka5Bz4sSkeQyvn/+EdwlInK+szkzgDM4pwgbA3NV9WAZ38dEIDs9aEz1mgkk4kxDz39ciHsq8WZgPnA/MF5ENuB8J9YMaAd0xblY97ynD1V1g4g8hjNDL1VEVuKcJrwEuApnWvvDqrop6KlP4pxWHAX8ICLL3T60AHq5+/N+H6aqe0RkHdAH2Cgia4BsYJuq/j2o7XEReQu4zy2yCRimVHakZUz1ehfIcv878NqsQlR1H06gPAysA3rgrF7RE+dC3H/gTAcvE1WdBiQBc3G+wxoHdABSgatVtcgpOFXNxlmJYgLwjfveY4DuwEbOfedUFik4F1HHAbcB9wC/KqHtF+7+B1VdUkIbYwAQ1ZImGBljTPUTkbnAzcBDqvqy1/0xoc1CyxjjGRHph7Pm4nGgvXshtTElsu+0jDE1TkReBRrhfH9WB/gvCyxTFnakZYypce4qH36c68Omqeo/PO6SCRMWWsYYY8KGzR40xhgTNiy0jDHGhA0LLWOMMWHDQssYY0zYsNAyxhgTNv4/DXsa24hYEe4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(num_heun[:end,1],num_heun[:end,2]/m0,'-',label='implicit Heun')\n",
    "plt.plot(num_rk2[:end,1],num_rk2[:end,2]/m0,'-',label='explicit RK2')\n",
    "plt.xlabel(\"velocity\")\n",
    "plt.ylabel(\"mf/m0\")\n",
    "plt.plot(250*np.log(0.25/(.25-.05*t[:end])), (.25-.05*t[:end])/.25, label = \"analytic\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "425.49491546291245"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "num_heun[end,0] #position when mass is 0.05kg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The root of the function is 425.49491546291245 m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Solve for the mass change rate that results in detonation at a height of 300 meters. Create a function `f_dm` that returns the final height of the firework when it reaches $m_{f}=0.05~kg$. The inputs should be \n",
    "\n",
    "$f_{m}= f_{m}(\\frac{dm}{dt},~parameters)$\n",
    "\n",
    "where $\\frac{dm}{dt}$ is the variable we are using to find a root and $parameters$ are the known values, `m0=0.25, c=0.18e-3, u=250`. When $f_{m}(\\frac{dm}{dt}) = 0$, we have found the correct root. \n",
    "\n",
    "Plot the height as a function of time and use a star to denote detonation at the correct height with a `'*'`-marker\n",
    "\n",
    "Approach the solution in two steps, use the incremental search [`incsearch`](../notebooks/04_Getting_to_the_root.ipynb) with 5-10 sub-intervals _we want to limit the number of times we call the function_. Then, use the modified secant method to find the true root of the function.\n",
    "\n",
    "a. Use the incremental search to find the two closest mass change rates within the interval $\\frac{dm}{dt}=0.05-0.4~kg/s.$\n",
    "\n",
    "b. Use the modified secant method to find the root of the function $f_{m}$.\n",
    "\n",
    "c. Plot your solution for the height as a function of time and indicate the detonation with a `*`-marker."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def heun_step2(state,dmdt,rhs,dt,etol=0.000001,maxiters = 100):\n",
    "    '''Update a state to the next time increment using the implicit Heun's method.\n",
    "    \n",
    "    Arguments\n",
    "    ---------\n",
    "    state : array of dependent variables\n",
    "    rhs   : function that computes the RHS of the DiffEq\n",
    "    dt    : float, time increment\n",
    "    etol  : tolerance in error for each time step corrector\n",
    "    maxiters: maximum number of iterations each time step can take\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    next_state : array, updated after one time increment'''\n",
    "    e=1\n",
    "    eps=np.finfo('float64').eps\n",
    "    next_state = state + rhs(state,dmdt)*dt\n",
    "    ################### New iterative correction #########################\n",
    "    for n in range(0,maxiters):\n",
    "        next_state_old = next_state\n",
    "        next_state = state + (rhs(state,dmdt)+rhs(next_state,dmdt))/2*dt\n",
    "        e=np.sum(np.abs(next_state-next_state_old)/np.abs(next_state+eps))\n",
    "        if e<etol:\n",
    "            break\n",
    "    ############### end of iterative correction #########################\n",
    "    return next_state\n",
    "\n",
    "def bisect(func,xl,xu,es=0.0001,maxit=50):\n",
    "    '''bisect: root location zeroes\n",
    "    root,fx,ea,iter=bisect(func,xl,xu,es,maxit,p1,p2,...):\n",
    "    uses bisection method to find the root of func\n",
    "    arguments:\n",
    "    ------\n",
    "    func = name of function\n",
    "    xl, xu = lower and upper guesses\n",
    "    es = desired relative error (default = 0.0001 )\n",
    "    maxit = maximum allowable iterations (default = 50)\n",
    "    p1,p2,... = additional parameters used by func\n",
    "    returns:\n",
    "    -------\n",
    "    root = real root\n",
    "    and a list of [fx, ea, iter]\n",
    "    fx = function value at root\n",
    "    ea = approximate relative error ( )\n",
    "    iter = number of iterations'''\n",
    "    xr = xl\n",
    "    ea = 100\n",
    "    for iter in range(0,maxit):\n",
    "        xrold = xr\n",
    "        xr = (xl + xu)/2\n",
    "        if xr != 0:\n",
    "            ea = abs((xr - xrold)/xr) * 100\n",
    "        else:\n",
    "            ea = abs((xr - xrold)/1) * 100\n",
    "        test = func(xl)*func(xr)\n",
    "        if test < 0:\n",
    "            xu = xr;\n",
    "        elif test > 0:\n",
    "            xl = xr;\n",
    "        else:\n",
    "            ea = 0;\n",
    "        if ea <= es:\n",
    "            break\n",
    "\n",
    "    root = xr\n",
    "    fx = func(xr);\n",
    "    return root,[fx,ea,iter]\n",
    "def mod_secant(func,dx,x0,es=0.0001,maxit=50):\n",
    "    '''mod_secant: Modified secant root location zeroes\n",
    "    root,[fx,ea,iter]=mod_secant(func,dfunc,xr,es,maxit,p1,p2,...):\n",
    "    uses modified secant method to find the root of func\n",
    "    arguments:\n",
    "    ----------\n",
    "    func = name of function\n",
    "    dx = perturbation fraction\n",
    "    xr = initial guess\n",
    "    es = desired relative error (default = 0.0001 )\n",
    "    maxit = maximum allowable iterations (default = 50)\n",
    "    p1,p2,... = additional parameters used by function\n",
    "    returns:\n",
    "    --------\n",
    "    root = real root\n",
    "    fx = func evaluated at root\n",
    "    ea = approximate relative error ( )\n",
    "    iter = number of iterations'''\n",
    "\n",
    "    iter = 0;\n",
    "    xr=x0\n",
    "    for iter in range(0,maxit):\n",
    "        xrold = xr;\n",
    "        dfunc=(func(xr+dx)-func(xr))/dx;\n",
    "        xr = xr - func(xr)/dfunc;\n",
    "        if xr != 0:\n",
    "            ea = abs((xr - xrold)/xr) * 100;\n",
    "        else:\n",
    "            ea = abs((xr - xrold)/1) * 100;\n",
    "        if ea <= es:\n",
    "            break\n",
    "    return xr,[func(xr),ea,iter]\n",
    "\n",
    "def incsearch(func,xmin,xmax,ns=50):\n",
    "    '''incsearch: incremental search root locator\n",
    "    xb = incsearch(func,xmin,xmax,ns):\n",
    "      finds brackets of x that contain sign changes\n",
    "      of a function on an interval\n",
    "    arguments:\n",
    "    ---------\n",
    "    func = name of function\n",
    "    xmin, xmax = endpoints of interval\n",
    "    ns = number of subintervals (default = 50)\n",
    "    returns:\n",
    "    ---------\n",
    "    xb(k,1) is the lower bound of the kth sign change\n",
    "    xb(k,2) is the upper bound of the kth sign change\n",
    "    If no brackets found, xb = [].'''\n",
    "    x = np.linspace(xmin,xmax,ns)\n",
    "    f = func(x)\n",
    "    sign_f = np.sign(f)\n",
    "    delta_sign_f = sign_f[1:]-sign_f[0:-1]\n",
    "    i_zeros = np.nonzero(delta_sign_f!=0)\n",
    "    nb = len(i_zeros[0])\n",
    "    xb = np.block([[ x[i_zeros[0]+1]],[x[i_zeros[0]] ]] )\n",
    "\n",
    "    \n",
    "    if nb==0:\n",
    "      print('no brackets found\\n')\n",
    "      print('check interval or increase ns\\n')\n",
    "    else:\n",
    "      print('number of brackets:  {}\\n'.format(nb))\n",
    "    return xb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_m(dmdt,m0=0.25, c=0.18e-3, u=250):\n",
    "    ''' define a function f_m(dmdt) that returns \n",
    "    height_desired-height_predicted[-1]\n",
    "    here, the time span is based upon the value of dmdt\n",
    "    \n",
    "    arguments:\n",
    "    ---------\n",
    "    dmdt: the unknown mass change rate\n",
    "    m0: the known initial mass\n",
    "    c: the known drag in kg/m\n",
    "    u: the known speed of the propellent\n",
    "    \n",
    "    returns:\n",
    "    --------\n",
    "    error: the difference between height_desired and height_predicted[-1]\n",
    "        when f_m(dmdt)= 0, the correct mass change rate was chosen\n",
    "    '''\n",
    "    dt = 0.002 \n",
    "    T = 1000\n",
    "    N = round(T/dt)\n",
    "\n",
    "    # time array\n",
    "    t = np.linspace(0, T, N)\n",
    "\n",
    "    x0 = 0    # initial position\n",
    "    v0 = 0    # initial velocity\n",
    "    m0 = 0.25    # initial mass\n",
    "\n",
    "\n",
    "    #initialize solution array\n",
    "    num_heun = np.zeros([N,3])\n",
    "    num_rk2 = np.zeros([N,3])\n",
    "\n",
    "    #Set intial conditions\n",
    "    num_heun[0,0] = x0\n",
    "    num_heun[0,1] = v0\n",
    "    num_heun[0,2] = m0\n",
    "    num_rk2[0,0] = x0\n",
    "    num_rk2[0,1] = v0\n",
    "    num_rk2[0,2] = m0\n",
    "\n",
    "    for i in range(N-1):\n",
    "        if num_heun[i][2] < .05:\n",
    "            end = i\n",
    "            break\n",
    "        num_heun[i+1] = heun_step2(num_heun[i],dmdt, rocket, dt)\n",
    "        #num_rk2[i+1] = rk2_step(num_rk2[i], rocket, dt)\n",
    "\n",
    "\n",
    "    error = num_heun[end,0] - 300\n",
    "    return error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(.05,.15,40)\n",
    "f_m_vectorized = np.vectorize(f_m)\n",
    "y = f_m_vectorized(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x14eaa820>]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbAAAAECCAYAAACfTIJUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxU9b3/8dc3C4QlYUtCCGFfguxI2BQFwaVWRCxLKVir1atWbb3Wa2u32/6qvbZ622vrWkUvdbvuooLWuoHsEvZ9XwMhYcsC2fP9/THDEEL2TObMmXk/H488jmedT745zptz5jvfY6y1iIiIuE2E0wWIiIg0hAJMRERcSQEmIiKupAATERFXUoCJiIgrRTldgJvEx8fb7t27O12GiIhrrF69+pi1NqEpjq0Aq4fu3buTnp7udBkiIq5hjNnfVMfWLUQREXElBZiIiLiSAkxERFxJASYiIq6kABMREVdSgImIiCspwJqYtZZ5azN4dUWT9SQVEQlL+h5YEzp1ppj731zHV9uziYmOYGzveLrHt3K6LBGRkKArsCbUqnkUmblFABSWlPPQexvQ89dERPxDAdaEoiMjeGzqYCKMZ37FnhO8seqgs0WJiIQIBVgTG5TShn+7vKdv/r8WbCUzp9DBikREQoMCLADuv7Iv3Tu0BCCvqJRfz9ukW4kiIo2kAAuAmOhI/jh1sG/+861HWbDxiIMViYi4nwIsQEb37MCsUV1987/9YDMnTxc7WJGIiLspwALooWv7kRQXA8Dx08U8PH+LwxWJiLiXAiyA4mKieWTKQN/8e2szWLg9y8GKRETcSwEWYFf278jkIcm++V+9v4n8olIHKxIRcScFmAN+e31/2rWMBiDjVAGP/3ObwxWJiLiPAswBHVo357fXD/DNv7xiP6v2nXCwIhER91GAOeSGoclckZoAgLXw83c3UFhS5nBVIiLuoQBziDGGP9w4iFbNIgHYk32aP36iW4kiInWlAHNQctsW/HpSf9/83GX7WLLzmIMViYi4R8ADzBiTaoy5zxjzqjFmmzGm3BhjjTHT6rDvLGPMYmNMjjEm3xiTboy5xxhT4+/R0P0CYeaILkzsl+ibf/Cd9eQUlDhYkYiIOzjxBv4j4AlgNpAKmLrsZIx5GngNSAMWA58BfYGngHeMMZH+3C9QjDE8OnWQr1fikZxCfvvBJidLEhFxBScCbBPwOPBdoDewqLYdjDFTgbuBTGCwtXaStfZGoA+wFbgRuNdf+wVaYmwMj35nkG9+3rrDzN9w2MGKRESCX8ADzFo7x1r7M2vtW9ba3XXc7Rfe6c+ttTsrHOsonis6gIequCXY0P0C7lsDOzH14hTf/K/e38TRXD12RUSkOo6/cdfGGJMCDAeKgbcrr7fWLgIygCRgdGP3c9JvJ/enc9sWAOQUlPDgO3qCs4hIdYI+wIBh3ulma21BNdusqrRtY/ZzTFxMNP89fQjG+6ng1zuyeXXlAWeLEhEJUm4IsB7e6f4atjn7Lt+jwrKG7ueoMb06cPvYc+X814Kt7D122sGKRESCkxsCrLV3WtO7eL53GuuH/c5jjLnD2+0+PTs7u8ZC/eWBq1NJ7egpqaCkjPvfXEdpWXlAXltExC3cEGBnu9nX98Oghu53Hmvt89baNGttWkJCQmMOVWcx0ZH85btDiI70/ArrDp7i6a/q2t9FRCQ8uCHA8rzT1jVsc3ZdXoVlDd0vKAxIbsO/X9nXN//XL3awbLdG6RAROcsNAbbPO+1WwzZdKm3bmP2Cxl3jejGye3sAyi385P/WkpmjrvUiIuCOAFvrnQ4wxrSoZpsRlbZtzH5BIzLC8OSsYcS3bg7Asfxi7n19DSX6PExEJPgDzFp7EFgDNAOmV15vjBkHpOAZbWN5Y/cLNh3jYnjye8OI8H6il77/pEatFxHBBQHm9ah3+idjTO+zC40xicAz3tk/WmsrX5o0dL+gMqZXB372rX6++ReX7GXBhiMOViQi4jwT6JEejDEXcy48APrj6ca+E/A9lthaO7rSfs/gGf6pEPgcKAEmAnHAPGCatfaCJ0I2dL+qpKWl2fT09Dr9nv5mreWOV1bz2ZajALRqFsmHPx5Lr4Sa+qiIiDjLGLPaWpvWJMd2IMDGA1/Vtp219oJR6o0xs4B7gEFAJLANeAl4tqarqIbuV5mTAQae4aUmP7WE/cfPANC3Y2vm3XMpLZtFOVaTiEhNQirA3MzpAAPYcjiXG59ZSlGpJ3dvGJrME98dijF1eiqNiEhANWWAueUzMPHqnxzHI1MG+uY/WHeYV1fUNFqWiEhoUoC50PS0Lswc0cU3//v5W1h74KSDFYmIBJ4CzKV+N3kAA5LjACgps/zo1TVk5xU5XJWISOAowFwqJjqSZ2cPJy7G04EjM7eQe/QlZxEJIwowF+vaoSV/+94w3/PDvtl7gj8s2OpsUSIiAaIAc7nxqYk8eE2qb37usn28s/qQgxWJiASGAiwE/GhcL749KMk3/8v3N7LxUI6DFYmIND0FWAgwxvD4tCH07egZlaO4tJw7X0nnWL46dYhI6FKAhYhWzaP4+/fTiPV26jicU8g9r6lTh4iELgVYCOkR34q/zTzXqWPl3hM8+rFGrheR0KQACzFX9EvkgavOPcn5paV7eX+tOnWISOhRgIWgu8f35poBHX3zD72rTh0iEnoUYCEoIsLw5xlD6Z3o6dRRVFrObf9YxZGcAocrExHxHwVYiGrdPIrnv39upI6svCJum5vO6aJShysTEfEPBVgI65nQmuduGk5UhKdXx5Yjudz3xlrKyvUIHRFxPwVYiLukd/x5j1/5fGsWj36s4aZExP0UYGFg5siu3Hl5T9/8nCV7eW2lniEmIu6mAAsTP/9WP67uf65n4n9+sJnFO7MdrEhEpHEUYGEiIsLwxMyhDOzseYZYWbnl7tfWsPNonsOViYg0jAIsjLRsFsWLPxhBUlwMAHmFpfzwH6s4rjETRcSFFGBhpmNcDHN+kEbLZpEAHDxRwB2vrKawpMzhykRE6kcBFoYGdm5z3piJq/ef1MC/IuI6CrAwdWX/jvz6uv6++S+2ZfEfb6+nXN8RExGXUICFsdvG9uDu8b188x+sO8zvPtqMtQoxEQl+CrAw9+A1qcwe1dU3//Ly/fzlsx0OViQiUjcKsDBnjOH3Nwzk+iHJvmVPfrmLOYv3OFiViEjtFGBCZIThz9OHMD41wbfskQVbeWvVQQerEhGpmQJMAGgWFcGzs4czsnt737KH3tvAJxuPOFiViEj1FGDi06JZJHNuSaN/J89oHeUW7ntjnYacEpGgpACT88TFRPPybSPpGd8KgOKycu54eTWr9p1wuDIRkfMpwOQC8a2b88rto0hu4xlyqqCkjFv/dxUbDp1yuDIRkXMUYFKlzm1b8Orto4hv3RyA/KJSbn7pG7Zl5jpcmYiIhwJMqtUzoTWv3j6Sti2jATh1poSb5qxkd3a+w5WJiCjApBb9kuJ4+YcjiW0eBcCx/GJmv7CSgyfOOFyZiIQ7BZjUanBKW/731hG0iPaMYJ+ZW8isOSs4klPgcGUiEs4UYFInad3bM+cHaTSL8pwyB08UMHvOSrLz9CwxEXGGAkzq7NLe8Tx308VERXiew7In+zTff3Elp84UO1yZiIQjBZjUy4R+HfnrzGF4M4xtmXnMnrNST3UWkYBTgEm9XTe4E49PG+Kb33w4lxl/X05mTqGDVYlIuFGASYNMHZ7CY1MH+57qvDv7NNP/vowDx9U7UUQCQwEmDTZjRBf+NnOY7zOxgycKmPbcMnYezXO4MhEJBwowaZTrhyTz/M3Dae7tnZiVV8SMvy9n46EchysTkVCnAJNGm9CvI3NvHUmrZp7viZ08U8KsF1bwzV4NACwiTUcBJn4xplcHXvu30bRp4Rl2Kq+olJtfWsmiHXoUi4g0DQWY+M3QLm15887RvgGAC0vKuf0fq/RQTBFpEgow8at+SXG8fdcYOrdtAUBJmeWe19fw1qqDDlcmIqFGASZ+1yO+FW/dNcb3UMxyCz97dwMvfL3H4cpEJJQowKRJdG7bgrfuGkP/TnG+ZX/4eCuPf7oNa62DlYlIqFCASZOJb92cN+4czcju7X3Lnv5qN7+et4mycoWYiDSOAkyaVFxMNP/44Ugm9Ev0LXtt5QH+/c11FJeWO1iZiLidAkyaXItmkfz9+8O5YWiyb9lH6w9zxyvpFBSXOViZiLiZAkwCIjoygv+ZMZSbx3TzLVu4PZubX1pJzpkSBysTEbdSgEnAREQY/t/kAfxkQm/fslX7TnLjs0vZd+y0g5WJiBspwCSgjDH89OpUfjOpv2/ZnuzTTHlmKSv3HHewMhFxGwWYOOK2sT14atYw3yDAp86UcNOLK3l39SGHKxMRt1CAiWMmDU7mjTvODT1VUmZ54O31PP7pNsrVzV5EaqEAE0cN69qOefdcQr+kWN+yp7/azb3/t0Y9FEWkRgowcVxKu5a8fdcYrkhN8C37eGMmM59fTlZeoYOViUgwU4BJUIiNieaFm9O45ZLuvmXrD+Uw5amlbD2S61xhIhK0FGASNKIiI/jd5AE8fMMAIiMMAIdzCpn+3HIWbs9yuDoRCTYKMAk63x/TnZduGUHr5lEA5BeVcts/0nllxX6HKxORYBI2AWaMmWWMWWyMyTHG5Btj0o0x9xhjwqYN3GRc3wTe/dElvueKlZVbfjNvE4/M36KBgEUECJMAM8Y8DbwGpAGLgc+AvsBTwDvGmEgHy5NqpCbF8v49lzA4pY1v2Zwle7nr1dWcKS51sDIRCQYhH2DGmKnA3UAmMNhaO8laeyPQB9gK3Ajc62CJUoPE2BjevGMM1wzo6Fv22ZajfPfvK8jKVQ9FkXAW8gEG/MI7/bm1dufZhdbao8CPvLMP6VZi8GrRLJJnZg/n3y7r4Vu2MSOHKU+rh6JIOAvpN21jTAowHCgG3q683lq7CMgAkoDRga1O6iMywvCr6/rzyJSB5/VQ/M4zyzT8lEiYCukAA4Z5p5uttQXVbLOq0rYSxG4a3e28HooFJWU88PZ6Hnx7vT4XEwkzoR5gZ+851dT/+kClbSXIne2h2DOhlW/Z26sPccNTS9l5NM/BykQkkEI9wFp7pzU9bCrfO42taqUx5g5vl/v07OxsvxYnDZeaFMtH947lxmGdfct2ZuUz+amlvJ1+0MHKRCRQQj3AjHfa4C8OWWuft9amWWvTEhISat9BAqZV8yj+MmMIj00d7HssS0FJGQ++s4EH3tItRZFQF+oBdvZ+Uusatjm7TveeXMgYw4wRXfjw3rH0qnBL8d01h5j81FJ26JaiSMgK9QDb5512q2GbLpW2FRdKTYrlw3vH8p2Lz91S3JWVz5SnlzJ/w2EHKxORphLqAbbWOx1gjGlRzTYjKm0rLuW5pTiUx6cNJibac2qfKS7j3tfX8l8fb6W0rNzhCkXEn0I6wKy1B4E1QDNgeuX1xphxQAqeUTqWB7Y6aSrT07rwwT1j6d6hpW/Z81/v4eaXvuF4fpGDlYmIP4V0gHk96p3+yRjT++xCY0wi8Ix39o/WWv3zPISkJsXywb1jmdgv0bds2e7jTH5qKRsOnXKwMhHxl5APMGvtO8CzeEbb2GiM+cgY8x6wE+gPzMMzqK+EmDYtPA/JvP/Kvhhvf9SMUwVMe245b6mrvYjrhXyAAVhr7wZm47mdOA64BtiFZxDfqdbaMgfLkyYUEWG478o+vPiDNGJjPKN3FJeW87N3NvCr9zdSVKo/vYhbGWv1bKW6SktLs+np6U6XIQ2079hp7nxlNdsrdK1P7RjL49MHMzilrYOViYQuY8xqa21aUxw7LK7ARAC6x7fivbsvYdLgTr5l24/mceMzy/jTP7dRWKKrMRE3UYBJWGnVPIonvzeM398wgBbRnueYlpVbnl24m0lPLmHNgZMOVygidaUAk7BjjOHmMd359N8vZ3TP9r7lu7LymfbsMv6wYIuuxkRcQAEmYatrh5a8fvtoHp4ykFbNPFdj5RZeWLyXa/+6mFX7TjhcoYjURAEmYS0iwvD90d349P7LuaxPvG/53mOnmfH35TwyX1djIsFKASYCpLRrycs/HMkfvzOIWO/DMq2FOUv2MunJJWw8lONwhSJSmQJMxMsYw8yRXfnXTy9nXN9zj87ZlZXPjc8s5YnPd1Ci8RRFgoYCTKSSTm1aMPfWEfzhxoG09H42VlpueeLznUx9dhm7svSIFpFgoAATqYIxhtmjuvHJfZeR1q2db/mGQzlc97clvLhkL+XlGgRAxEkKMJEadOvQijfvHMND1/ajWaTnf5ei0nIenr+FWXNWsP/4aYcrFAlfCjCRWkRGGO4a14sPf3wp/TvF+Zav2HOCa574mjmL91CmqzGRgFOAidRRv6Q45t1zKfde0ZsI7+j2hSXlPLJgK995dhnbM/XZmEggKcBE6qFZVAT/cU0q8+65lH5Jsb7l6w+eYtKTi/mfz3ZQXKqeiiKBoAATaYDBKW356MdjeeCqvr7PxkrKLH/9YieTnlzMuoN6aKZIU1OAiTRQdGQEP57YhwU/GcvFXc89jmXH0Xy+88xSfv/RFvIKSxysUCS0KcBEGqlPx1jevusSfnt9f98I9+UWXlq6l4l/XsSH6w+j5+6J+J8CTMQPIiMMt17ag39VGlMxK6+In/zfWmbPWcmurHwHKxQJPQowET/q0t4zpuJfZw4lIba5b/my3ce59q9f86d/buNMcamDFYqEDgWYiJ8ZY7hhaGe+eGAct17a3dflvqTM8+DMq/7yNf/clKnbiiKNpAATaSJxMdH89voBzP/xZQyvMBxVxqkC7np1NTe/9A1bj+Q6WKGIuynARJpY/+Q43r5zDI9NG0z7Vs18yxfvPMa3/7aYn72znqO5hQ5WKOJOCjCRAIiIMMxI68KXD4xj9qiuvtuK1sJb6YcY//hC/uezHfp8TKQejO7D111aWppNT093ugwJAdsz83j0k60s3J593vLE2OY8cHVfpg3vQuTZlBNxMWPMamttWlMcW1dgIg5ITYpl7q0jeeW2kecNSZWVV8TP393IdX9bzPLdxx2sUCT4KcBEHHRZnwQW/OQyHps2mMQK3e63ZebxvRdW8NO31nEsv8jBCkWClwJMxGGR3s/HFj44nvuv7Ot7CjTAe2symPjnRby2cr8eoClSiQJMJEi0bBbFfVf24csHxnPd4E6+5TkFJfzq/U1MfW4Zmw/nOFihSHBRgIkEmaQ2MTw962Lm3jqCru1b+pavPXCK659cwsPzt5BfpN6KIgowkSA1PjWRf91/OT+e0JvoSE+PxHILLy7Zy8Q/L+TVFfspKi1zuEoR56gbfT2oG704ZVdWPr+Zt4nle87vmdipTQx3X9GbGWkpNI+KrGZvEec0ZTd6BVg9KMDESdZaPlh3mEcWbL2gZ2KnNjHcPb4XM0Z0UZBJUFGABQkFmASDguIyXlu5n+cW7eZYfvF565LiYrj7il7MSOtCTLSCTJynAAsSCjAJJueCbM8FV2RJcTH89Oq+TL04RSN6iKMUYEFCASbBqKC4jNe/OcBzi3aTnXd+kPVLiuUX376Iy/vEY4yCTAJPARYkFGASzApLynh95QGeWbj7giuysb3j+cW3+zEguY1D1Um4UoAFCQWYuMHpolJeWLyH57/ew5nic93sjYEbh3bmgWtS6dy2hYMVSjhRgAUJBZi4SVZuIU98sZM3Vx2krMIwVM2iIrh5dDfuuLwniXExDlYo4UABFiQUYOJGu7Ly+OMn2/l869HzljeLiuC7aV24c1xPUtq1rGZvkcZRgAUJBZi42Yo9x3n0462sP3T+eIpREYYpwzrzo/G96JXQ2qHqJFQpwIKEAkzczlrL51uzeOqrXaw/eOq8dcbAtwd24u4reqmzh/iNAixIKMAkVFhrWbrrOE99tZMVe05csP6K1ATuGteLkT3aq/u9NIoCLEgowCQUrd5/gqe+3MVX27MvWDesa1vuGteLqy7qSIS+EC0NoAALEgowCWWbMnJ4ZuEuPtmUSeW3hV4JrbhzXC+mDO1Msyg9xELqTgEWJBRgEg72ZOfzwuI9vLs6g+Ky8vPWJcXFcNvYHswc2YXYmGiHKhQ3UYAFCQWYhJOs3EJeXLqX11ccIK/SAzRjm0cxc2QXbrm0h74ULTVSgAUJBZiEo9zCEl5bcYCXlu69YKzFyAjDdYM6cftlPRic0tahCiWYKcCChAJMwllhSRnvr83ghcV72JN9+oL1I3u05/axPbhSHT6kAgVYkFCAiUB5uWXhjixe+HrvBU+IBugR34ofjOnG1OEp+pxMFGDBQgEmcr5NGTnMWbyH+RuOUFp+/ntJq2aRTE/rws1jutFTI3yELQVYkFCAiVTtSE4Bc5ft4/WVB8grLL1g/bi+CdxyaXfG9UnQ7cUwowALEgowkZqdLirlvbUZ/GPZPnZl5V+wvkd8K24a3Y1pF6fQpqVuL4YDBViQUICJ1M3ZoarmLtvLF9uyLvhidPOoCCYNTmbWqK5c3LWthqsKYQqwIKEAE6m/A8fP8MqKfbyx6mCVtxf7JcUye1RXbhjWmTh1+gg5CrAgoQATabgzxaXMW3uY11buZ/Ph3AvWt4iOZPKQZL4/phsDO2s0/FChAAsSCjCRxrPWsuFQDq+vPMCH6w9TUFJ2wTYju7fnh2O7c1X/JCLV6cPVFGBBQgEm4l+5hSXMW5vBaysOsP1o3gXrU9q14JZLujNjRBfdXnQpBViQUICJNA1rLWsOnOTl5ftZUMV3ylo3j2J6Wgq3XNKdbh1aOVSlNIQCLEgowESaXmZOIa+s2MdrKw9w6kzJeeuMgbG945k8JJlrBibpqswFFGBBQgEmEjgFxZ6xF19aurfK75Q1i4pgQmoik4cmM6FfIjHRkQ5UKbVRgAUJBZhI4FlrWbzzGC8u2cuiHRc+NRo8txivHtCRyUOSGds7nqhIPXQzWCjAgoQCTMRZh08VMH/DYT5Yd7jKrvgAibHNmTo8henDUzQGYxBQgAUJBZhI8NiVlc+H6w/z0frD7D124eNdAIZ3a8f04SlcN7iTRsZ3SMgEmDGmFTAFGOH9GQa0ABZYayfVYf9U4DfABKADkAl8DPzeWnvE3/tVpgATCT7WWjZm5PDBOs+V2bH8ogu2aREdybWDkpiR1oVRPdpr6KoACqUAGwqsrWJVrQFmjBkHfIIn8NYAO4EhQD8gGxhrrd3hr/2qogATCW4lZeUs3J7N2+kH+XJb1gXd8QF6JrRi9igNKBwooRRgvYBfAunAajxXYM9RS4B5r9x2AUnAj621T1VY99/AA3jCKc1W+IUaul91FGAi7nEsv4h5azN4Z/UhtmVe+CXpmOgIrh+czOzR3RiS0kZXZU0kZALsghc35hbgf6k9wO4FngQWWmuvqLQuEtgO9AKus9Z+3Nj9qqMAE3Gfs7cY30o/yLy1h8kvunBA4YGd45g9qhuThyTTqnmUA1WGrqYMMLf0NZ3inb5aeYW1tgx4o9J2jd1PREKEMYbBKW15ZMogVv5yIo9+ZxADkuPO22ZTRi6/eG8jaY98zp2vpPPemkOcOlPsUMVSV275p8Yw73RVNetXVdqusfuJSAhq1TyK743syswRXVh38BSvrTzAR+sPU1RaDkBBSRmfbj7Kp5uPEhlhGNWjPVf378hVA5Lo3LaFw9VLZUEfYMaYOKC9d3Z/NZsd8E57NHY/EQl9xhiGdW3HsK7t+PV1F/HO6kO8sergeSN+lJVblu0+zrLdx/ndR1sY2DmObw1IYtLgZLrHazzGYBD0AQZU/CZi1V/2gLNnXawf9juPMeYO4A6Arl27Vl+liLhS25bNuP2yntx+WU92ZeXz2ZajfLo5k3UHT5233aaMXDZl5PLf/9rBwM5xXDcomUmDO9GlfUuHKpc6d+IwxjwGTG7Aa0y01mZUc8xbqKUThzGmM3DIOxttrb3gE1hjTB9gB1BsrW3emP1qok4cIuHjaG4hn205yr+2HGX57mOUlFX9XjmkS1smDerEdYM7kazbjBdoyk4c9bkCSwZSG/Aajf2iRcX+r62AnCq2aV3Ftg3dT0SEjnEx3DS6GzeN7kZuYQlfbs1i/oYjfL0jm+Kyct926w+eYv3BU/zh462M6N6OG4elcN2gTvqOWQDUOcCstTcBNzVhLdW9bq4x5gSez7O6ARuq2KyLd7qvsfuJiFQWFxPNlGGdmTKsM7mFJXy2+SjzNxxm8c5j531ZetW+k6zad5LffbiZiRclcuOwzoxPTaRZlFs6fLuLGz4DA8/oHRPxDD9VVRCNrLCdP/YTEalSXEw0U4enMHV4CqfOFPOvzUf5aMNhlu0+Tpk3zIrLyvlkUyafbMqkXctoJg1OZsqwzlzcta2+MO1HbgmwD/AE0WzgxYorvF9Inumdfd9P+4mI1Kpty2bMGNGFGSO6kJ1XxEfrD/P+2gw2Zpz7xOLkmRJeWbGfV1bsJ7lNDNcMTOLagZ0Y3q0dkREKs8Zwy0gcrfGMYZgE3GutfbrCuseB/8BzFTW80lBSDdqvOurEISJ1sSsrj/fWZPDBusNknCqocpv41s25ZkBHrh3YiVE92xMdos8wC6mhpIwx7wOdvLMJQE/gFJ5hnc562Fq7oNJ+FQflXc25QXkvAo7hGZS34jEatV9VFGAiUh/l5ZZv9p3g/TUZ/HNzJjkFJVVu17ZlNBP7dWRCv0TG9o4PqQ4goRZg+/B0qqjJrdbauVXsmwr8J57bgu2Ao3gei/L/6vA4lXrvV5kCTEQaqqSsnBV7jvPJpkz+tTmTY/lVD1UVYeDiru0Y1zeBcakJDExuQ4SLbzWGVIC5mQJMRPyhrNySvu8En2zK5NPNmRzJKax22w6tmnF53wTGpyYwvm+i667OFGBBQgEmIv5WXm7ZkJHDl9uyWLQjmw2HTlHd23JkhGFE93ZceVFHrryooyuGtFKABQkFmIg0teP5RSzZdYxF27NZtCOb46erHxW/d2JrJl6UyFUXdWRY1+Ds1agACxIKMBEJpPJyy+bDuXy1PYsvtmWxvtL4jBXFxUQxumcHLu0dzyW9OtA7sXVQfOdMARYkFGAi4qSs3EK+3JbF51uzWLIrm8KS8mq3TYhtzstL6JwAAAgfSURBVCW9OnBpr3jG9Org2KDDCrAgoQATkWBRUFzGst3H+HzrUb7clsXR3KIat+/eoSWX903g8j4JjOnVIWBPnlaABQkFmIgEI2stu7NPs3z3MZbuOs7yPcer/c4ZQHSkIa1be0+g9Y2nf6e4JrvdqAALEgowEXGDsnLL1iO5LPMG2jd7T1BQUlbt9gmxzRnfN4GJFyUytk8Crf14daYACxIKMBFxo6LSMlbvO8miHZ6ejdsyq3+CVHSkYVSPDkzol8iEfomN7qqvAAsSCjARCQVZuYV8vfMYi3Zks2RnNifPVH+7sWdCKyakJjLhokTG9OxQ71uNCrAgoQATkVBTVm7ZcOgUX27L4outWWw5klvldr0SWvHFA+PrffxgeSKziIiEmMgIw7Cu7RjWtR0PXJ3KkZwCvtqWzZfbjrJk1zFfV/0J/RIdrvRCCjAREfHp1KYFs0Z1ZdaorhSWlLF893G+3JbFtYM61b5zgCnARESkSjHRkVzRL5ErgvDqCyA0n6AmIiIhTwEmIiKupAATERFXUoCJiIgrKcBERMSVFGAiIuJKCjAREXElDSVVD8aYbGB/A3ePB475sZxQp/aqH7VX/ai96qcx7dXNWpvgz2LOUoAFiDEmvanGAwtFaq/6UXvVj9qrfoK1vXQLUUREXEkBJiIirqQAC5znnS7AZdRe9aP2qh+1V/0EZXvpMzAREXElXYGJiIgrKcBERMSVwjrAjDGzjDGLjTE5xph8Y0y6MeYeY0yD2qW+xzPGzDXG2Bp+tgWyfn//fv48njFmfC1tVfGna6V9G9XODeWv9jLGpBpj7jPGvGqM2WaMKffWPa0p63Dr+dWQ9jLGRBtjJhpj/myMWWGMOWKMKTbGZBhj3jHGjK9h37A8v4Lh/StsH2hpjHkauBsoBL4ASoCJwFPARGPMdGttWYCOtxTYVcXyI4GqvzZB0F6ZwD9qOORI4CJgN3Cwmm3q3c4N5ef2+hFwXyDrcPn51ZD2Ggd85v3vTGA1cBroD0wFphpjHrbW/mcNxwi788vLufcva23Y/eA5Ia23gftUWN4R2OJdd19THw+Y6113i5P1u6W9ajnmZu9+v/RXOwdRe90OPAbMAHoBC73HmNZE56Xbz696txcwAXgHuKyKdd8FSr3HuELnV+N+b3/W3+SNHYw/QLq3kW6uYt24Co0b0ZTHa8QJ4Nf63dJeNRxvjHefUqBzFesD/QbTpH+ferzBNPS8dPX51dD2quUYc7zHeFHnV+N+b3/W3+SNHWw/QIq3gYqAFtVsc8i7zSVNebyGnAD+rt9N7VXDMV/wbj+/mvUBe4MJxN+nLm8wDa3D7edXQ9urDse4x3uMT3V+Nfz39nf94diJY5h3utlaW1DNNqsqbdvUx7vCGPMXY8zzxpiHjTHX1PBBpr/rr00wtpePMaYlnls8AC/Wsnl92rmhAv338Xcdbj+/mkof77Smz7PC6fyqyLH3r3DsxNHDO61pVPkDlbZt6uPdXMWyLcaYmdbajU3wevURjO1V0XQgFsgC5teybX3auaEC/ffxdx1uP7/8zhiTBNzinX23hk3D6fyqyLH3r3C8AmvtnZ6uYZt87zS2iY+3DvgJMMB7nGRgErAeT++nz40xnf34eg0RTO1VlR96py9ba0uq2aYh7dxQgf77+LsOt59ffmWMiQJeBdoAX1hrP6pis3A8vyAI3r/C8QrMeKfW6eNZa5+otOg0sMAY8xmwCBgN/AK41x+v10BB014XHMiY3sDl3tmXqtuuge3c4LLOvqwfjtUYDa3D7eeXvz2Hp3v3QeCmqjYI0/MrKN6/wvEKLM87bV3DNmfX5dWwTVMdD2ttMfCod/bbTf16tQjm9jp79bXcWru1Dq99nlrauaEC/ffxdx1uP7/8xhjzV+A2PN8Lm2itzazP/iF+flUrkO9f4Rhg+7zTbjVs06XStoE83llnv8Ve+RK8qV6vOv5+Pb8czxgTybl777V13qhJde3cUPu800D9ffxdR0P3a6hAv16dGGP+jOf2WDae8NrZwEOF6vlVm4C8f4VjgK31TgcYY1pUs82IStsG8nhndfBO8ystb6rXq06wttc1eP7nOA28WYfXrU517dxQgf77+LsOt59fjWaMeQz4KXAcuMpau6URhwvV86s2AXn/CrsAs9YeBNYAzfD0YDuPMWYcnu8qZALLA328CmZ4p6sqLmzC16tSELfXbd7pm9baxrw5VNnODRXov4+/63D7+dVYxpg/Ag8CJ/GE1/pGHjIkz686CMz7V1N/8S4Yf4BpnPu2d+8KyxM5NyTRfZX2eRTPZfGjfjreUDw9diIrLY/C86+/Mu9+1/jj9dzeXpX2j8fzRchav+zYmHYOlvaq4vgLqdsXTRvUzm4/vxrRXg97tzsJDK9j7WF5fjXm9/bn+eWXBnXjD/CMt6EKgI+A94Ac77L3q/jDzPWum+un403xrjuO518abwP/BDK8y8uAn/mrfre3V6V97/dut7UOdTeqnYOhvYCLgRUVfnK92+6ouNyf7ezm86sh7QVM9m5j8Vw1zK3m5yGdX8Hz/uXXRnXbDzALz0jKuXg+S1mNZ8iYC8bgqukEaODxegBPAMu8f/RC7x9zJ54u4bX+C7A+r+f29qq03wbvsR+sQ82Nbmen2wsYz7k312p/mqCdXXl+NaS98HxRudZ9gIU6v4Ln/ct4DyQiIuIqYdeJQ0REQoMCTEREXEkBJiIirqQAExERV1KAiYiIKynARETElRRgIiLiSgowERFxJQWYiIi40v8H4CQ2/hcv5JMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.07917657470703125, [-0.1290890616291449, 7.708739166148406e-05, 13])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bisect(f_m,.079,.08)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of brackets:  1\n",
      "\n",
      "Wall time: 1.92 s\n",
      "number of brackets:  1\n",
      "\n",
      "Wall time: 478 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[0.08888889],\n",
       "       [0.05      ]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time xl,xu = incsearch(f_m_vectorized, 0.05, .4, ns=50)\n",
    "\n",
    "%time incsearch(f_m_vectorized, 0.05, .4, ns=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The best estimate for the dmdt is [0.0791766] kg/s\n",
      "We reached a relative error of [6.88280036e-05] with 16 iterations\n"
     ]
    }
   ],
   "source": [
    "dmdt,out=bisect(f_m,xl,xu)\n",
    "print('The best estimate for the dmdt is {} kg/s'.format(dmdt))\n",
    "print('We reached a relative error of {} with {} iterations'.format(out[1],out[2]))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The best estimate for the dmdt is [0.07938595] kg/s\n",
      "We reached a relative error of [0.27532563] with 49 iterations\n"
     ]
    }
   ],
   "source": [
    "dmdt,out=mod_secant(f_m,xl,xu)\n",
    "print('The best estimate for the dmdt is {} kg/s'.format(dmdt))\n",
    "print('We reached a relative error of {} with {} iterations'.format(out[1],out[2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "dt = 0.002 \n",
    "T = 1000\n",
    "N = round(T/dt)\n",
    "\n",
    "# time array\n",
    "t = np.linspace(0, T, N)\n",
    "\n",
    "x0 = 0    # initial position\n",
    "v0 = 0    # initial velocity\n",
    "m0 = 0.25    # initial mass\n",
    "\n",
    "\n",
    "#initialize solution array\n",
    "num_heun = np.zeros([N,3])\n",
    "num_rk2 = np.zeros([N,3])\n",
    "\n",
    "#Set intial conditions\n",
    "num_heun[0,0] = x0\n",
    "num_heun[0,1] = v0\n",
    "num_heun[0,2] = m0\n",
    "num_rk2[0,0] = x0\n",
    "num_rk2[0,1] = v0\n",
    "num_rk2[0,2] = m0\n",
    "\n",
    "for i in range(N-1):\n",
    "    if num_heun[i][2] < .05:\n",
    "        end = i\n",
    "        break\n",
    "    num_heun[i+1] = heun_step2(num_heun[i],dmdt, rocket, dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaYAAAEaCAYAAABaefMNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3xV9f3H8dcnO5AQCFtGWCpDkBGGgiJFxV1RcKKAqC3ODn9WbbVWa13V1jprqwLOKm5rcVCggKhsAQFBdlgJMwkJCcn398e5mSQhuRn3Jvf9fDzyODnnfL/nfkLI/dxzzvd8vuacQ0REJFiEBToAERGR4pSYREQkqCgxiYhIUFFiEhGRoKLEJCIiQSUi0AHUdy1atHCdOnUKdBgiIvXK4sWL05xzLcvap8RUTZ06dWLRokWBDkNEpF4xs83l7dOlPBERCSpKTCIiElSUmEREJKgoMYmISFAJmsRkZrea2dtmttrM9phZrpmlmtmXZjbOzKyCvleZ2VwzO2BmGWa2yMxuNrMKfz5/+4mISO0JplF5vwFaASuBr4BMIAn4CTASGGNmlzjn8ot3MrNngZuAbGAmkOtr/www0szGOufySr+Yv/1ERKR2BVNiugJY6pzLLL7RzHrhJY6fAuOBV4rtuxQvuewETnfOrfNtbw3MAkYDtwBPlTqmX/1ERKT2Bc0lK+fcvNJJybd9FfCsb/WsUrvv9i1/U5BcfH12AZN9q3eVcWnO334iIqHruOPADNq1q9WXqS9vvEd8y+yCDWbWHhgA5ADvlO7gnJsDpABtgCHV7SciEpLy82HsWJg8GXbs8LZt3+6tjx0LtTCnXzBdyiuTmXUGfu5b/bjYrn6+5SrnXFY53RcC7Xxtv6pmvxrhnCM9PZ2DBw9y6NAh8vJ0G0ukNoWHh9OoUSOaNGlCfHw8FYyjkrIMHgy+6jZzO/Wl3/a1xOVkwQsvFO3/9tsafcmgS0xmNhEYDkQC7YFT8c7sHnbOvV+saWffstyyFsCWUm2r06/anHPs3r2bzMxMEhMTadOmDeHh4fpDEaklzjny8vLIyMggLS2NrKwsWrVqpb+5qvD9W61tkcSkS++jdcZenvzkSQamfF9if00Kxkt5Q/EGOVwFnO7bdi/wQKl2cb7lUfelisnwLeNroF8hM7vRN7R8UWpqagWHKSk9PZ3MzEySkpJo2rQpERER+gMRqUVmRkREBE2bNiUpKYnMzEzS09MDHVb9Mn8+h8MjuP3CO8iJiGJr0zY8OPJ6XLH9NS3oEpNz7nrnnAGNgF7AX4H7ga/N7LhiTQve0at6gdPffsVjfNE5l+ycS27ZssziuGU6ePAgiYmJhIeH+/vSIuKn8PBwEhMTOXjwYKBDqV+GDuXJ065hTSvvAlJ07mGe+PdfCt9IGTq0xl8y6BJTAedclnPue+fc/+GNojsZ7xmjAgUfe+KO6lykYF/xj0j+9qu2Q4cOERdX0cuKSG2Ki4vj0KFDgQ6jXvk6oSMvDhpduH7XnCkcv2drUYNaGPwQtImplIJnly40s0jf95t8y6QK+nUo1bY6/aotLy9PZ0siARQeHq4BR1VwMDuXX4+cjPM9OXPaxiWMX/yJt7NgVN4339T46wbd4Idy7McbMh4BJAK7gKW+fb3MLLacEXYDfculxbb5269G6J6SSODo769qfv/hKlL2e0/pJMRG8vjitwjDec8zPfdcrb1ufTljOh0vKe0H0gCcc1uBJUAUMLZ0BzMbjjeqbyewoGC7v/1ERELJJ99t5/2lKYXrD1/SmzbrV3mX7lJSKuhZfUGRmMzsNDO72syiy9g3FHjJt/pSqfp1D/uWj5pZt2J9WgEF6fyR0vX1qtFPRKTB23kgm9++v7Jw/ZL+7Tivd9s6e/1guZTXFe8+0jNmtgTvbCXet72nr82/8YaNF3LOTTez5/HKCK0wsy8pKsbaBPiAkgMmqtVPRKShy8933PHOcg5k5QLQrmks91/Uq05jCIozJmAO8CCwDDgBuAQ4G2gMvAuMds5dUNb9IOfcTcDVeJfnhgOjgPV4RVgvLa9CuL/9pHZ16tQJM2PTpk2BDqVcs2fPxsw444wzjtpnZjVyH2PChAmYGVOmTKn2sYo744wzKnXcgte///77a/T1JfhN+WoT89anAd6zs3+5vC9NYiKP0atmBcUZk3NuI3BfNfq/AbxRV/1EAmXKlClMnDiR8ePH13jSEvlhVzqPzFhTuP6z07syqHNinccRFIlJpMDMmTPJzc2lXS1XL64tq1evrpHjPPzww9x11120bVt31/UltGXn5nH7W8vIOeLdWu/Ztgm/OuuEgMSixCRBpWvXroEOoVq6d+9eI8dp27atkpLUqUdnrGH1Dq8qRlREGH+9oi9REYG52xMs95hEgPLvMRXcG5k9ezbz58/nnHPOoVmzZiQkJDBq1CiWLVtW2HbatGkMHDiQuLg4EhMTGTduHDt37jzqtaZMmYKZMWHCBNLS0pg8eTLt27cnJiaGrl278rvf/a7KVQIquseUm5vLiy++yIgRI0hMTCQ6OpqOHTtywQUX8Prrr5doW9Y9pk6dOjFx4kQApk6dWvhaBT9DXVq9ejWTJk2ic+fOxMTE0KxZM84880w++uijMtsf695beb/34tu/+OILRo4cSUJCAo0aNWLIkCHlvp5Uzaw1u3ll/qbC9d+e14MTWpdZKrROKDFJvfLxxx8zfPhw9u3bx6hRozjuuOP4/PPPGT58OOvWreOOO+7g+uuvp2nTpowaNYqoqChef/11zjzzTHJycso85r59+xg8eDBvv/02gwcPZtSoUaSmpvLQQw8xcuTIGilhs2/fPoYPH87PfvYzFixYQL9+/bjkkkvo3Lkz8+fP57e//e0xjzFmzBiG+uqSde3alfHjxxd+DRs2rNoxVtZbb71F3759efnll2ncuDEXXHABffr0Ye7cufz0pz/lvvv8vl1crpdeeolRo0aRkZHBeeedR/fu3fnmm2+4+OKLmT59eo2/XijZfTCbO95ZXrh+Zo9WXHtKRYVxap8u5QWRTnf9O9Ah+G3TI+fXyev85S9/4e2332bMmDEA5OfnM27cON58800uueQSUlNTWbZsGT17ek8Z7N27l1NOOYVVq1bxr3/9i2uuueaoY3700UcMHTqUxYsX07RpUwB27drFWWedxddff83999/PY489Vq24J0yYwIIFCzjllFOYPn06xx1XVI84OzubWbNmHfMYf/7zn5kyZQrz589n2LBhARn88N133zF+/HiioqL44IMPOPfccwv3rVq1inPPPZcHH3yQESNGMGLEiBp73ccee4xPP/2Uc845p3DbH//4R+69917uvvvuwv8PUjX5+Y5fv7OcPZneh7ZW8dE8NubkgFfI0BmT1CtXXHFFiTehsLAw7rzzTgBWrlzJAw88UJiUABITE/n5z715Jst78zcznn/++cKkBNC6dWueeuopAF544QWys7PL7FsZy5Yt46OPPiIuLo4PP/ywRFICiImJKfEGXxcmTpxY4lJg6a+pU6eW2e+hhx4iJyeHxx577KiYe/XqxZNPPgnAM8/U7GOAt956a4mkBHDnnXeSkJDA+vXr2bJlSzk9pSL/nLeBuetKDg1PbBwV4Kh0xiT1TOk3J4Bu3bpVuP/4448HYPv27WUes0+fPvTu3fuo7SNGjKBdu3akpKSwePHiwstoVTVjxgwAfvrTn1KVaVJq09ChQ0v8u5U2b948fvzxxxLb8vPzmTFjBmZW7hnK8OHDAViwoGareV1wwQVHbYuKiqJLly4sXbqU7du307Fjxxp9zYZuxbYDPP7Z2sL1n53elaHdWgQwoiJKTEGkri6H1Wft27c/alvxqUQq2l/eWU/nzuVPVNypUydSUlLYtm1bVUMttHmzN1lyTY3YqwnXX399hQMmJkyYcFRi2rNnT+FcRq1atarw+FWZQLMyyks6TZo0Acr/3UrZMg8f4ba3lpKb501ZcXL7BH59dmCGhpdFiUnqlbCwiq8+H2u/vwJ9zT0YFEwXER4ezrhx42r02Pn5FZelrK3fa6j6/Uer2JjmTeLdOCqcv13Zj8jw4Pk3VmKSkFdR+aOCfaXvC1VFUpI3wmnt2rXHaBncWrRoQWxsLFlZWTzzzDNVmvQyMjKS3NxcMjIyjuqXm5vLjh07ajpcKcdHy7czfXHRFYAHLz6JpOaNAxjR0YInRYoEyPLly1m5cuVR2+fMmUNKSgpxcXEMGDDA7+OPGjUKgA8//JC0tDS/jwPefRWAI0eOVOs4/oiIiODMM88EqPIQ7YJKHmvWrDlq3+effx6QnycUbd17iN++t6Jw/eK+x3FJ/6MvfweaEpOEPOccN910EwcOHCjclpqayu233w7AjTfeSGxsrN/H79evHxdeeCHp6emMHj36qLOD7Oxs/vOf/1TqWAVv8DVV+qiq7rvvPiIjI7n99tt56623cKWm1c7Pz2fmzJmFAz4KjBw5EoAHHnigxPNkq1at4tZbb639wIXcvHxue2sp6Ye9DwEdExvx4MUnBTiqsulSnoS8iy66iJUrV9K1a1fOOOMMjhw5wqxZszh48CADBw7kgQceqPZrTJkyhXPOOYd58+bRpUsXhg0bRsuWLdm+fTvLly8nISGhUhXVhwwZQps2bViyZAnJycn06tWLyMhIhg4dWlgVojYlJyczbdo0rrvuOq688kruuusuevbsSXx8PNu2beOHH34gLS2N3/zmNyVGSN5999288847fPzxx5x44okMGDCAnTt3snDhQi677DLy8/MLB4lI7Xj8s7Us3bIfgPAw46kr+hJfx1XDK0tnTBLymjVrxtdff83o0aNZsGAB//nPf2jevDn33HMPs2bNonHj6l9/T0xMZO7cuTz99NP079+fb7/9lvfee4+NGzdy2mmn8cgjj1TqONHR0cyYMYPzzz+fjRs38tprr/HSSy8xZ86casdYWVdccQUrVqzgtttuo1GjRsyZM4dPPvmEnTt30r9/f5566iluu+22En26du3K/Pnzueiii9i/fz///ve/OXDgAI8//jjTpk2rs9hD1czVu3jxfxsK1+8cdSL9OjYLYEQVs9Kn4lI1ycnJbtGiRZVqu3r1anr06FHLEUllaQqJ0BRqf4cp+7M4/29z2X/Im/hvxIkteWn8QMLCAjvS1MwWO+eSy9qnMyYRkQYqNy+fW99YUpiU2ibE8MRlfQOelI5FiUlEpIH68+drWVLsvtLTV/YLipJDx6LEJCLSAP13zS7+PqfovtIdZ59Icqe6n43WHxqVJyFrwoQJdT6PkUhd2L4/i1+/XTSVxRkntuRnp3cJYERVozMmEZEGJDcvn9veXMo+332lNk1ieLIe3FcqTolJRKQBefKLH1i0eR/gu690Vf24r1ScEpOISAMxa+1unp9dVBX+V2edwMB6cl+pOCWmOqbnxkQCpyH//W3bd4hf/mtZ4frpJ7Rk8vCuAYzIf0pMdSgiIqJEnTARqVs5OTlERDS8MV/ZuXnc9HrR80qtm0Tzl8tOrlf3lYpTYqpDCQkJ7Nmzp0F/ahMJVs459uzZQ0JCQqBDqXEPfPI9323zihBHhBnPXd2f5nHRAY7Kf0pMdSgxMZHDhw+zbds20tPTycvLU5ISqUXOOfLy8khPT2fbtm0cPnyYxMT6d8+lItMXb+ONb7YUrv/u/B4MSKrfP2PDO6cNYhERESQlJbFv3z727dvH9u3bjzlzp4hUT1hYGLGxsTRu3JhmzZo1qNlwv99+kN++XzS/0kUnH8f4UzsFLqAaosRUx8LCwmjevDnNmzcPdCgiUo8dyMpl8uuLOXzE+3B7fKs4Hr6kN2b1875ScQ3no4OISIjIz3f8+u1lbN5zCIDGUeG8cM0AGkc3jHMNJSYRkXrm+Tk/8uXq3YXrj489ma4t4wIYUc1SYhIRqUfmr0/jic/XFq7fcFpnzuvdNoAR1TwlJhGRemL7/ixufXMp+b7BvIM6JXLnOd0DG1QtUGISEakHsnPz+Nmri9mb6T2k3zI+mmeu6kdkeMN7G294P5GISAPjnOOe91ewIqXoIdpnr+pPqyYxAY6sdigxiYgEuVfmb+K9JSmF67+/sCeDOtfvh2grosQkIhLEvlqfxkOfri5cvzy5A+OGJAUwotqnxCQiEqS27j3EzW8sIc832qFvh6Y8cHGvBvEQbUWUmEREglBWTh43vrq4cCbalvHRvDBuANER4QGOrPYpMYmIBBnnHHe++x2rdxwEIDLceGFcf9okNMzBDqUpMYmIBJl/zN3Ax8u3F64/8NOT6n3F8KpQYhIRCSL/+yGVR/6zpnD96sEduXJQxwBGVPeUmEREgsSG1IwSlR2Sk5rx+wt7BTaoAFBiEhEJAgcO5XL91EUcyPIGO7RpEsNz4/oTFRF6b9Oh9xOLiASZ3Lx8bnpjMRvSMgGIjgjjxWsH0Co+NAY7lKbEJCISYA98/D3z1+8pXH/ispPp075pACMKLCUmEZEAmrZgE69+vblw/ZdnnsAFfY4LXEBBQIlJRCRA5q5L5Q8ff1+4fuHJx3HbyG4BjCg4KDGJiATA+t0Z3PR6Ubmhk9sn8PiYPg2+3FBl+DVBvJmFA5cBI4HjgPLu0Dnn3MhKHC8SOB04DxgKJAHNgVRgAfCMc252Bf2vAiYDfYBwYA3wCvC8cy6/pvuJiFTH/kM5XD91IenZRwBvBN6L1yYTE9nwyw1VRpUTk5k1Az4H+gPHSu2ukocdDnzh+34nsBjIBHoClwKXmtmDzrn7yojnWeAmIBuYCeTiJcxngJFmNtY5l1dT/UREqiM3L5+bXl/Cpj2HAIiJDOOf45Np3UDnVvKHP2dMDwEDgK14b+JrgIPVjCMfeBd4yjk3t/gOM7sceB2418xmOedmFdt3KV5y2Qmc7pxb59veGpgFjAZuAZ4qdUy/+omIVIdzjvs+XMlXPxaNwPvLZX05qV1CAKMKPuZcZU9qfB3MtgGxQC/n3M5aiero1/wnMAl42Tk3qdj2RXhJcrxzblqpPsOB2XjJp13xS3P+9itLcnKyW7Rokf8/nIiEjOdn/8ijM4rKDd1x9gnc8pPjAxhR4JjZYudccln7/Bn80AKYV1dJyWepb9m+YIOZtcdLLjnAO6U7OOfmAClAG2BIdfuJiFTHx8u3l0hKo/u14+YRGoFXFn8S03bgSE0HcgwFHyl2FNvWz7dc5ZzLKqffwlJtq9NPRMQvCzft5dfvLC9cH9IlkUcv1Qi88viTmN4FTjez2JoOpixm1gaYUOy1C3T2LTdTvi2l2lann4hIlW1IzeCGaYvIOeLdFejasjF/H5cckjXwKsuff5k/4J01/cvMWtVwPCWYWQTwGpAAzHTOfVxsd5xvmVnBITJ8y/ga6Fc8rhvNbJGZLUpNTa3gMCISyvZkHGbilIXs981C2yIuiikTB5HQKDLAkQU3f0bl/Q1YjzdybZ2ZLcY7wyhrkIArPljBDy/gDeHeCowrta/gHLhqozf871fIOfci8CJ4gx/8PY6INFzZuXncMG0Rm0sMCx9Ih8RGAY4s+PmTmCZQ9KYeD5xRQVuHN5quyszsKV/fncDIMgZbpPuWcZSvYF96sW3+9hMRqZT8fMev3l7Gki37ATCDp67oR98OoVuYtSr8SUwTazyKUszsCeA2vMoPIwueMyplk2+ZVMGhOpRqW51+IiKV8uiMNXy6ouiz9L3n92RUrzYBjKh+qXJics5NrY1ACpjZY8CvgD3AWc6578tpWjCEvJeZxZYzwm5gqbbV6Scickwvz9vI3/+3oXB9wqmduG6YxlFVRVANCzGzR4D/A/bhJaXl5bV1zm0FlgBRwNgyjjUc77mnnXj19qrVT0TkWD5avp0HPin6LH1mj9bce0HPAEZUP1UrMZlZlJmdYmZjfF+nmFmUn8d6EPgNsB8vKVXmbOVh3/JRMyt8Us03WvA53+ojZVRv8LefiEiZ5q9P49dvLytc79+xKU9f2Y/wMD2rVFX+VhePBO4HbuboIdUZZvY08AfnXG4lj3cR8Dvf6nrg1nIePFvjnHukYMU5N93MnserEL7CzL6kqBhrE+ADvHp+JfjbT0SkLCtTDvCzVxeTm+eNC+vWKo6XJwwkNkrVwv3hT3XxcOAT4Ey8odc7gA2+7zsDbYG7gYFmdl4lK3QnFvs+2fdVljnAI8U3OOduMrN5eElyOEXTV7xMBdNX+NtPRKS4LXsOMeGVhWQcLprCYtp1g2jayK+LR4J/RVwnA88CPwC/cM7NKLV/FPBX4ATgZufcCzUUa1BSEVeR0JWWcZgxz39VOIVFk5gIpk8+lRNal/lsvhRT00Vcr8WrmjCydFICcM59hnc2dQgY78fxRUSCXsbhI0x8ZWFhUoqOCOOlCQOVlGqAP4mpJzDLOZdSXgPfvlm+tiIiDUrOkXwmv7aYFSkHAAgzePrKfgzslHiMnlIZ/iSmSLyzoWM55GsrItJg5PmqOsxdl1a47aHRvTlbD9DWGH8S02bgtIqGhfv2DaPiCt4iIvWKc47fvr+CT74rmoHnl2eewJWDOgYwqobHn8T0Ed7Iu6lmdlThJzNLwBvZ1hb4sHrhiYgEB+ccf/p0NW8t3Fq4bfwpSdw2UpP91TR/nmN6DLgSuAw418w+BjbiFWztAlyI92zTNl9bEZF675n/rucfczcWrl/Svx2/v7CXJvurBf7UyttjZj8B3sB73uhqiqqNF/yGFgJXOef21kiUIiIBNGX+Rp744ofC9bN7tuaxS/sQpqoOtcKvyg/OufXAIDMbhvdwaju8pLQNmOOcm1dzIYqIBM70xdu4/+Oi+nfDurXg6av6EREeVKVGGxS/ElMBXwJSEhKRBmnGyh3cOb2olnT/jk158doBREeo1FBtUsoXESnD3HWp3PbmMvJ9Nyq6t4nnlQmDaBRVrc/zUglKTCIipSz4cQ83TFtETp5XMrNzi8a8OmkwCY30aGZdOGbqN7M8vMENPZ1zP/jWK8s55/TxQkTqjYWb9jJp6kKyc72k1DYhhteuH0zL+OgARxY6KpM0jKLRdpT6vjJ9RUTqhSVb9jHh5W85lON9/m4VH80bNwyhXdPYAEcWWo6ZmJxzYRWti4g0BN9t28/4l74l05eUWsR5Salzi8YBjiz0KMmISMhbmXKAcf/8hnTfnErNG0fxxg2D6dYqLsCRhaYqJyYzu9bMTq1EuyFmdq1/YYmI1I3VOw4y7qVvOJjtJaWmjSJ57frBmr4igPw5Y5oCXF+JdpOAV/w4vohInVi3K51x//yG/YdyAW+iv9cmDaZH2yYBjiy01ealPA18EJGgtX53Olf+4xv2ZOYAEB8dwauTBnNSu4QARya1mZjaAxm1eHwREb+s3ZnOFS9+TVrGYQDioiOYOmkQJ3c4asIECYBKPWNUxr2ibhXcP4oAegAj8Yq5iogEje+3H+Tqf37NPt/lu8ZR4bwycSD9OzYLcGRSoLIPv06hqII4wFDfV3kMyAf+7F9YIiI1b8W2A4x76RsOZHlJKS46gqnXDWRAkqZEDyaVTUzTKEpM44EfgfnltM0BUoAPnXPLy2kjIlKnlm7Zx7Uvf0u6b/RdfIx3T6mvLt8FnUolJufchILvzWw8MM85d11tBSUiUpMWb97L+JcXknG42JBwDXQIWv7UseuMBjWISD3xzYY9XDdlYWFFh8TGUbw2aTA9j9OQ8GDlzwy2m2sjEBGRmvbV+jQmTV1EVm5BmaEoXr9+CCe20cOzwawy1cU7+r5Ncc7lFVuvFOfcFr8iExGphi++38XNbywh54hXJbygIKvKDAW/ypwxbcIbYdcT+MG37ipoX5yr5GuIiNSYD5am8Ot3lpPnm+WvTZMY3rxRBVnri8okjS14CSa31LqISNB5dcEm7vtoFc73LpXUvBGvTRpMh8RGAY1LKq8y0150qmhdRCQYOOd4bvaPPP7Z2sJtJ7aO59VJg2jVJCaAkUlV6TKbiNR7zjkembGGv8/ZULitb4emTJk4kKaNogIYmfhDiUlE6rW8fMfvPljJm98WjbMa2q05L16TTONovcXVR/7Mx9TczAaZWYtS29uZ2WtmtsLMPjazfjUXpojI0XKO5POLfy0rkZTO6tmal8YPVFKqx/z5zd0N/BLoB6QBmFk0MA/oiFcnrxcwzMz6OOe21lCsIiKF0rNzmfzaEuatTyvcdkn/djx2aR8iwjU5d33mz29vBLDBOfddsW1XAEnALOBM4G9AAnBLtSMUESlld3o2l//96xJJafwpSfx5zMlKSg2AP7/BdnhFXIs7H28I+Q3Ouf86534BbADOqWZ8IiIlbEjN4JLnvuL7HQcLt91x9gncf1EvwsI0P2lD4M+lvGb4LuEVcwqw1jm3sdi2pXhzMomI1IilW/Zx3ZSFhXMphYcZD1/Sm8uSOwQ4MqlJ/iSmLKBw4IOvRFE74KVS7XIAjdMUkRoxc7VXYig71ysxFBsZzrNX9+Mn3VsHODKpaf4kpu/xBja0cM6lAVfjXcb7X6l2HYBd1YxPRIR/LdzCPe+vLCwxlNg4ipcnDNRcSg2UP4lpGvAcsMjMluDdX0oHPixoYGYxQH9gTk0EKSKhyTnHX774gb/9d33htg6JsUydOIguLVWMtaHyJzG9CAwBrsUbHp4OTHLOHSzW5iKgEUpMIuKn7Nw87pz+HR8t3164rddxTXhl4kBaxavEUEPmz3xM+cAEM7sPaAWscc6VnjjwB2A08HX1QxSRULMn4zA/e3UxizbvK9x2+gktefaqfsTHRAYwMqkLfj8a7Ztnqcy5lpxzy4Bl/h5bRELX+t0ZXDdlIVv2HircNm5IR+6/sJeeUQoR1a7ZYWZt8UblgTeZ4I7qHlNEQtNXP6bx81cXczD7CABm8NvzejBpWGfM9IxSqPA7MZnZDcAdQLdS29cBf3bO/bOasYlICHln0Vbufm8FR3wj72Ijw3nqir6c3atNgCOTuuZXYjKzKcA1eHXxHLDd931b4ATg72Y21Dk3sYbiFJEGKi/f8fhna3lhTlFBmVbx0bw0fiC92ycEMDIJFH+qi1+JNyIvFbgJaOSc6+Cca483Em8ysBu41syuqMlgRaRhOZidy/VTF5ZISt3bxPPBzUOVlEKYP2dMN0Eb9C4AABaBSURBVOBVdfiJc+774jucc4fxzpbmAkuAG4G3qh2liDQ4G1IzuGHaIn5MzSzc9pPurfjblf2I05QVIc2f335fYE7ppFScc+57M5sNDPI3MBFpuGav3c2tby4l3TfIAeCmM7ry67NPJFyFWEOeP2MvGwF7KtFuLxBb2YOa2YlmdrtvssE1ZpZvZs7MxlSi71VmNtfMDphZhpktMrObzazCn8/ffiLiH+cc//jfBq6bsrAwKUVHhPHUFX2585zuSkoC+HfGlAIMMjNzzrmyGpg3rnMg3qCIypoM3F7VYMzsWbx7XdnATCAXr6r5M8BIMxvrnMurqX4i4p/s3DzueW8F7y1NKdzWNiGGF69J1v0kKcGfM4PPgM7A42YWXnqn72zjUaALMKMKx10JPA5cjjcE/ZjljMzsUrzkshPo45y7wDk3GjgeWI1XfeKoyQr97Sci/tm69xBjX1hQIikNSGrGh7dokIMczco56Sm/gzfNxTK8GWo3A68DG/GGjXcBrsRLXPuBvv5Ore67RzUcGOucm15Om0XAAGC8c25aqX3Dgdl4yaedr5RStfqVJTk52S1atKhKP5tIKJm9dje/+Ncy9vvmUAK4PLkDD1zci+iIoz7bSogws8XOueSy9vlTK2+LmZ0LvAN0Au4p/XrAVuAyf5NSZZhZe7zkkuOLpXScc8wsBa8qxRDgq+r0E5Gqyc93PDNrPX/58gcKPv9GhBn3XtCTa09JUiUHKZdfYzKdc9+Y2fHAWLyzmnZ4CWkb3iW4d3xDx2tTP99ylXMuq5w2C32x9aMowfjbT0Qq6cChXH759jL+u2Z34bbWTaJ57ur+DEhKDGBkUh9Up4jrYeA131cgdPYtN1fQpqDIbOdi2/ztJyKVsDLlAJNfX8zWvUWf+4Z0SeTpK/vTMj46gJFJfVGfn2IrmCUss4I2BdNxxNdAv0JmdiPew8N07Nix4ihFQoRzjncWbePeD1dy+EjRrdmfDe/C/519oiqDS6VVp4hrFHApcAbQnqKaebOBd+vgUl7BBeqqjd7wv18h59yLeBMmkpyc7PdxRBqKjMNH+N37K/hgWdETInHREfx57Mmcc5KKsErV+FvE9VTgDaADRW/0BSYBD5vZ1c65edWMryLpvmVF8ysX7Esvts3ffiJShlXbD3DLG0vZmFZ0EeKE1nG8MG6Apj8Xv1Q5MZlZL+BzvAoQG4A3gU2+3Z0oeg5phpkNds6tqpFIj1bwmkkVtOlQqm11+olIMc45Xv16M3/8ZDU5eUWX7i5Lbs/9F/WiUVR9vlMggeTP/5wH8JLSw8C9pZ/zMbPf+9rcA/wBOGZJIT8t9S17mVlsOSPsBpZqW51+IuJzICuX30z/jhmrdhZuaxwVzkOje3Nxv3YV9BQ5Nn/uRg4H1jrnflvWw6fOuXzn3O+AtXj3n2qF7xmpJUAU3rD1EnwPyrbHe1B2QXX7iYhnyZZ9nPfU3BJJqWfbJnx86zAlJakR/iSmWLw39mNZAsT4cfyqeNi3fNTMCmfSNbNWwHO+1UfKSKD+9hMJWUfy8nnqy3WMfWEBKfuLLjSMPyWJ9246VfeTpMb4cylvLd5MtcfSFlhX2YOaWX+KkgJAT9/yT2Z2R8FG59yQYt9PN7Pn8QrArjCzLykqxtoE+ACvKGsJ/vYTCVWb92Tyi38tY+mW/YXb4mMieHxMH845qTJvByKV509iegF4zjd1+vyyGpjZUOB0qlYItQkwuIztx1fUyTl3k5nNA27Gu8wYDqwBXgaeL++sx99+IqGk4NmkP3y8isycomL7yUnN+MvlfemQ2CiA0UlDVeUirgBm9iTeTLbPUVTEFbxReVfjVe7+h3Pu1zUTZvBSEVdpqPZm5nD3e9/x2apdhdsiwoxfnnUCPx/eVXMnSbVUVMTVn+ri1ZmjyDnnGtQYUiUmaYjm/JDKHe8sJzW96Dn5Li0b89Tl/TRNhdSIGq0uztEP1NZVXxGpZenZufzp0zW8+e2WEtuvGZLEPef1IDZK01RI7fNn2gsVvBJpgOauS+Wud1eUGHHXIi6Kx8eczIjurQIYmYSaBnVZTUSqzjtLWs2b35acPm1Ur9Y8NLo3LeJUEVzqlhKTSAj73w+p3PXud2w/kF24rVmjSB746Ulc0KetJvOTgFBiEglBB7NzebiMs6RzerXhwYtP0rxJElBKTCIhxDnHf1bu5P6PVrG72Ig7nSVJMFFiEgkRKfuzuO+DlcwsNt05wLkneWdJupckwUKJSaSBy8t3TPlqE098vpZDxao3tIyP5v4Le3Fe7zY6S5KgosQk0oCtTDnA3e+tYEXKgRLbrx7ckTvP6U5CbGSAIhMpnxKTSAN0ICuXv375A9MWbCYvv6i6y/Gt4nj4kt4kd0oMYHQiFVNiEmlA8vMd7y7ZxqMz1pCWkVO4PSoijNt+0o0bT+9KVISekZfgpsQk0kCsTDnAfR+uZEmxqSkATu3anD9efJLmS5J6Q4lJpJ7bl5nDnz9fyxvfbqF4TeY2TWL43QU9OL+3hoBL/aLEJFJPHcnL582FW3ni87XsP5RbuD0y3Lj+tC7cMqIbjaP1Jy71j/7XitQzzjlmr03loU9Xs353Rol9p5/Qkvsv7KnLdlKvKTGJ1CPfbz/Inz5dzbz1aSW2t28Wy30X9OSsnq112U7qPSUmkXpg98Fs/vz5Wt5ZvK3EfaTGUeHcNKIbk4Z1JiZScyVJw6DEJBLEMg4f4Z9zN/Di/zaUqNoQZnDFoI788swTVHBVGhwlJpEglJ2bx2tfb+a52T+yNzOnxL7hJ7TknvN6cGKb+ABFJ1K7lJhEgsiRvHymL97GUzPXsaPYHEkAJ7aO557zezD8hJYBik6kbigxiQSB/HzHpyt38OTnP7AhLbPEvnZNY/nFmcczul87IsJVtUEaPiUmkQDKz3d8/v0u/jZzHd/vOFhiX4u4KG4Z0Y0rB3ckOkIDGyR0KDGJBEB+vjdh39P/Xceanekl9sXHRPDz4V2ZcGonPSArIUn/60XqUF6+45PvtvPMf9ezrtTDsTGRYUw4tTM/H96Fpo2iAhShSOApMYnUgdy8fD5e7iWk0veQGkWFc80pSdxwWhfNIiuCEpNIrco8fIS3Fm7l5XkbSdmfVWJfXHQE409NYtKwLiQ21hmSSAElJpFasDs9m6lfbeK1r7dwICu3xL74mAgmDu3MdUM76ZKdSBmUmERq0I+pGfzjfxt4b0kKOXn5JfYlNo5iwqmdGH9qJ01pLlIBJSaRanLOMW99GlO/2sTMNbtL1LIDSGreiOtP68KY/u2JjdKwb5FjUWIS8VPG4SO8t2QbU7/axI+pmUftP7lDU35+ehfO7tWG8DBV/BapLCUmkSrakJrBtAWbmb54GxmHjxy1f2T3Vtx4ehcGdU7UFBQiflBiEqmE3Lx8Zq7exRvfbuV/P6QetT8uOoIxA9pzzSlJdNUkfSLVosQkUoGNaZm8tXAL7y7eRlpGzlH7u7RszPhTOnHpgPbEqUqDSI3QX5JIKdm5ecxYuZO3Fm7h6w17j9pv5l2uG39qJ4Z1a6HLdSI1TIlJBG9k3dKt+3l/SQofLd9+1LNHAK2bRHNZcgcuS+5Ah8RGAYhSJDQoMUlI25SWyftLU/hgWQqb9xw6an94mDHixFZcOagDw09oqWknROqAEpOEnL2ZOXzy3XbeX5rC0i37y2zTITGWKwZ2ZMyA9rRuElPHEYqENiUmCQl7Mg7z2apd/GflDr76cQ95+e6oNvExEZzfuy0/7duOwZ0TCdOzRyIBocQkDVZq+mE+W7WTT1fs4OsNeygjFxERZpxxYisu6d+On3RvRUykKjOIBJoSkzQom9IymblmN198v5NvN+4tMxkB9O/YlNH92nF+n+NU2VskyCgxSb12JC+fxZv3MXPNbmau3lVmaSDwhngnJzXj3JPacm7vNrRNiK3jSEWkspSYpN5JyzjM/PVp/HfNbmavTS1zaDd4yWhQp0TO692Wc05qo0EMIvWEEpMEvezcPL7duJd569OYuy6N1TsOlts2JjKMYd1aMrJHK0b2aEWreCUjkfpGiUmCTm5ePqu2H+SrH9OYty6NRZv3kXMkv9z2bRNi+El3LxGd2rWFBjCI1HNKTBJwh3KOsGzLfr7dtJeFm/ayZPN+snLzym0fEWb069iU0473zox6tm2iskAiDYgSk9Qp5xwp+7P4btsBlm3dz7cb97Iy5QBHyhs+53N8qziGdmvBace3YHCX5iqYKtKA6a9batXezByWb9vP8q37+W7bAb7btr/MKt2ltWsay+DOiZzarQXDurWgTYLuFYmEipBPTGZ2FTAZ6AOEA2uAV4DnnXPl39iQEvLyHZv2ZLJmRzprdh5k9Y50Vu84SMr+rEr1P7F1PAM7N2Ngp0QGdkrkuKYazi0SqkI6MZnZs8BNQDYwE8gFRgLPACPNbKxzrvybHSHoSF4+2/ZlsSEtgw2pmfywK501O9NZuzOdwxUMUCguPjqC3u0T6NO+KQOSmpGc1IxmeshVRHxCNjGZ2aV4SWkncLpzbp1ve2tgFjAauAV4KmBBBsiRvHx2HMgmZX8Wm/dksiEtkw2pmWxIzWDL3kPk5lV8P6i4qIgwerZtQt8OTenTPoGTOzSlc/PGqkMnIuUK2cQE3O1b/qYgKQE453aZ2WRgNnCXmT3dkC7p5eU79mbmkJp+mN3p2ew6mM22fVmk7Mvylvuz2Hkwu8wip8fSukk03ds0oXvbeHr4ll1axBEVoakiRKTyQjIxmVl7YACQA7xTer9zbo6ZpQDtgCHAV3UbYeXk5Tsyso9wICu3xNf+rJzC7/dm5JCacdiXiA6zJ+NwufXjKqt1k2g6t2hMl5ZxdG0ZR4828ZzYJp7mcdE184OJSEgLycQE9PMtVznnyrs7vxAvMfWjhhPTwexcnvhsLXnOkZfvfR3Jd+QXLIttP3wkn6ycPA7l5JGdm0eW7+tQTl6FD51WV8v4aNo3i6V9s0Z0adGYLi0b06VFHJ1bNtZQbRGpVaH6DtPZt9xcQZstpdrWmMO5+UxdUNFL166mjSJpFR9Ny/hoWsXH0K5pLO2bxdKuWSztmsZyXNNYVU8QkYAJ1cQU51uWXYrak+FbxpfeYWY3AjcCdOzYscovHl6DN/7joiNIiI086qtpo0iaxEbSrFFUYRJqGR9N87gooiOUdEQkeIVqYirIDH7dbXHOvQi8CJCcnFzlYzSKCuf3F/YkIswICzPCzQgPK/kVEWaEmREdGU5sZDiNosKJiQwnNqpoPToiTKV4RKTBCdXElO5bxlXQpmBfegVt/BITGc7EoTV+hVBEpEEI1XG8m3zLpAradCjVVkRE6kCoJqalvmUvMyuv9s3AUm1FRKQOhGRics5tBZYAUcDY0vvNbDjQHq8qxIK6jU5EJLSFZGLyedi3fNTMuhVsNLNWwHO+1UcaUtUHEZH6IFQHP+Ccm25mz+NVFl9hZl9SVMS1CfABXjFXERGpQyGbmACcczeZ2TzgZmA4RdNevIymvRARCYiQTkwAzrk3gDcCHYeIiHjMuWpW9AxxZpZKxaWNKtICSKvBcCS46PfbsOn3Wz1JzrmWZe1QYgogM1vknEsOdBxSO/T7bdj0+609oTwqT0REgpASk4iIBBUlpsB6MdABSK3S77dh0++3lugek4iIBBWdMYmISFBRYhIRkaCixBQAZnaVmc01swNmlmFmi8zsZjPT76MeM7MTzex2M3vNzNaYWb6ZOTMbE+jYpHrMLNLMRprZE2b2tZntMLMcM0sxs+lmdkagY2xIdI+pjpnZs8BNQDYwk6L6fPHA+8BY51xe4CIUf5nZX4Hby9g11jk3va7jkZpjZmcCX/hWdwKLgUygJ3CSb/uDzrn7AhBeg6NP6HXIzC7FS0o7gT7OuQucc6OB44HVwGjglgCGKNWzEngcuBzoBswJbDhSg/KBd4HTnXNtfX+7lzvnegNXAHnAvWY2IqBRNhA6Y6pDZrYIGACMd85NK7VvODAbL2m1UwHZ+s/MZuMVB9YZUwNnZv8EJgEvO+cmBTqe+k5nTHXEzNrjJaUc4J3S+51zc4AUoA0wpG6jE5FqKpjpun1Ao2gglJjqTj/fcpVzLqucNgtLtRWR+uF433JHQKNoIJSY6k5n37KiSuRbSrUVkSBnZm2ACb7VdwMYSoOhxFR34nzLzAraZPiW8bUci4jUADOLAF4DEoCZzrmPAxxSg6DEVHfMt9RoE5GG4wW8xz22AuMCHEuDocRUd9J9y7gK2hTsS6+gjYgEATN7Cm8k3k5gpHNuZ4BDajCUmOrOJt8yqYI2HUq1FZEgZGZPALcBqXhJaV2AQ2pQlJjqTsFw0l5mFltOm4Gl2opIkDGzx4BfAXuAs5xz3wc4pAZHiamOOOe2AkuAKGBs6f2+B2zb410WWFC30YlIZZjZI8D/AfvwktLyAIfUICkx1a2HfctHzaxbwUYzawU851t9RFUfRIKPmT0I/AbYj5eUdGWjlqgkUR0zs+eAyXhFXL+kqIhrE+ADYIyKuNZPZtafog8Y4BX4jAfWAXsLNjrnVNmjnjGzi4APfauLgFXlNF3jnHukbqJquJSYAsDMrgJuBnoD4cAa4GXgeZ0t1V++qQ9mHaudc86O1UaCi5lNAF6pRNM5zrkzajeahk+JSUREgoruMYmISFBRYhIRkaCixCQiIkFFiUlERIKKEpOIiAQVJSYREQkqSkwiIhJUlJhEgpiZnWFmzsxmBzoWkbqixCQSYGa2yZd8OgU6FpFgEBHoAESkQt8CPYBDgQ5EpK4oMYkEMefcIbxaiiIhQ5fyRALEzCaYmaNoVuONvkt6BV+dyrvH5NvnfJcBw8zsV2a2ysyyzGybmT1pZo18bZuZ2V99bQ+b2Toz+1UFcZmZXWFmn5tZmq/PFjP7hy43Sl3QGZNI4KwHpgJjgMbAu0BGsf0ZZXUqwxvABcBs3zFPB34J9DCzq4Gv8abfmAck+vY/YWYxzrk/FT+QmUUCbwGXAFl4UzzsAk4CrgcuNbOznXOLqvrDilSWqouLBJiZbcI7a+rsnNtUat8ZeFNplJhOwXfmstG3uhb4iXNuu29fB2Ap0BxYiXcp8BrnXLZv//nAJ0A60MZ3ubDguI/gTYb3P+Bq59y2YvtuAZ4GfgS6O+eO1MCPL3IUXcoTqf9uK0hKAM65rcBrvtUkYHJBUvLt/zfwHd5ZVHLBdjNLBG7DO1MbWzwp+fo9A/wb6AqcWzs/iogSk0h9lwv8t4zt633LRc65tDL2r/Mtjyu2bQQQi3d2truc15vjW55S1UBFKkv3mETqt53lXFIruD+1rYx9xffHFNvWxbc83zcooyItKxmfSJUpMYnUb/nV3F9cuG+5Fm/AREW+qcJxRapEiUlECmz1LVc45yYEMhAJbbrHJBJ4Ob5loD8ofol3z+pMM2sa4FgkhCkxiQReim/ZI5BBOOd2Ac8CTYGPzKx76Ta+h3WvN7PWdR6ghIxAf0ITEXgfOAN43cw+B/b7tv8mALHciTdS7zJgpZktw3teKgbogJc8o3zLXQGIT0KAEpNI4D0DNAGuxqvgEO3b/se6DsQ5lwtcbmavA9cBg4A+eA/j7sCrMvEh3kO2IrVClR9ERCSo6B6TiIgEFSUmEREJKkpMIiISVJSYREQkqCgxiYhIUFFiEhGRoKLEJCIiQUWJSUREgooSk4iIBJX/B+lbSxIqETBSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(t[:end],num_heun[:end,0],'-',label='implicit Heun')\n",
    "plt.xlabel(\"time\")\n",
    "plt.ylabel(\"position\")\n",
    "plt.scatter([t[end]],[num_heun[end,0]], marker=(6,2,0), c='r')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "1. Math 24 _Rocket Motion_. <https://www.math24.net/rocket-motion/\\>\n",
    "\n",
    "2. Kasdin and Paley. _Engineering Dynamics_. [ch 6-Linear Momentum of a Multiparticle System pp234-235](https://www.jstor.org/stable/j.ctvcm4ggj.9) Princeton University Press \n",
    "\n",
    "3. <https://en.wikipedia.org/wiki/Specific_impulse>\n",
    "\n",
    "4. <https://www.apogeerockets.com/Rocket_Motors/Estes_Motors/13mm_Motors/Estes_13mm_1_4A3-3T>"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}