Permalink
Cannot retrieve contributors at this time
Name already in use
A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
Project-3/CompMech03-IVPs_project.ipynb
Go to fileThis commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
701 lines (701 sloc)
124 KB
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "# Initial Value Problems - Project\n", | |
| "\n", | |
| "" | |
| ] | |
| }, | |
| { | |
| "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": [ | |
| "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": 2, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "#IMPLICIT\n", | |
| "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" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 3, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "#EXPLICIT\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": 4, | |
| "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 = np.array([state[1], u/state[2]*dmdt, -dmdt])\n", | |
| " return dstate" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 5, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "y0 = 0\n", | |
| "v0 = 0\n", | |
| "m0 = .25\n", | |
| "\n", | |
| "dm = .05\n", | |
| "dmdt = .05\n", | |
| "N = 25\n", | |
| "\n", | |
| "t = np.linspace(0, 4, N)\n", | |
| "dt = t[1] -t[0]\n", | |
| "u = 250\n", | |
| "mf = .05\n", | |
| "m = np.linspace(mf, m0, N)\n", | |
| "\n", | |
| "num_sol_imp = np.zeros([N,3])\n", | |
| "num_sol_exp = np.zeros([N,3])\n", | |
| "\n", | |
| "num_sol_imp[0,0] = y0\n", | |
| "num_sol_imp[0,1] = v0\n", | |
| "num_sol_imp[0,2] = m0\n", | |
| "\n", | |
| "num_sol_exp[0,0] = y0\n", | |
| "num_sol_exp[0,1] = v0\n", | |
| "num_sol_exp[0,2] = m0\n", | |
| "\n", | |
| "for i in range(N-1):\n", | |
| " num_sol_imp[i+1] = heun_step(num_sol_imp[i], simplerocket, dt)\n", | |
| " num_sol_exp[i+1] = rk2_step(num_sol_exp[i], simplerocket, dt)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 6, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "<matplotlib.legend.Legend at 0x7f74a4215128>" | |
| ] | |
| }, | |
| "execution_count": 6, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEaCAYAAABEsMO+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3hUZfbA8e9JJaEkhIQakoDURVEBFUSlWVBBULABLljAgr2irj8riG1ld0EREBBBXUQF0XVRaWJdwYpEQCEJvbdAQtr5/XEnIWUyySSTTMr5PM88k7n3fe+cDCQn997zvq+oKsYYY0xVE+DvAIwxxhh3LEEZY4ypkixBGWOMqZIsQRljjKmSLEEZY4ypkoL8HUB1Fx0drQkJCf4OwxhjqpU1a9bsVdUYT20sQZVTQkICq1ev9ncYxhhTrYhIcklt7BKfMcaYKskSlDHGmCrJEpQxxpgqqcomKBGZICLqetzvod0wEVklIodEJFVEVovIWBHx+L2VtZ8xxpjKUSV/GYvIGcCDgMeJAkVkCjAP6AasAj4D2gGTgQUiEujLfsYYYypPlUtQIhIKzAZ2AYs8tBsC3AbsBDqr6gBVvRxoCyQClwO3+6qfMcaYylXlEhTwFPAX4BbgkId2D7ueH1LVjbkbVXUXcKvr5Tg3l+zK2s935s0jp2UcKgHkxMbBvHkV9lbGGFNdVakEJSJnAfcBb6nqYg/tYoGuQAbwbuH9qroS2AY0BbqXt59PzZtHzo2jCdi6BUEJ2LYFHTPGkpQxxhRSZRKUiNQB3gD2A3eV0Px01/NvqppWTJvvC7UtTz/fefRRAo4XfGs5dgwefbRC3s4YY6qrqjSTxHigPXCNqu4toW0r17OnkcgphdqWp5/vpKS43awpKUiFvKEpC1XlyJEjHD58mGPHjpGdne3vkIypkoKCgoiIiCAqKoqgIN+mlCqRoETkbOBuYKGq/rsUXeq5no96aJPqeq7vg34FiMgYYAxAXFych0O5ERcHyUXzY0a9MEJUEbE05W+qyu7duzl69ChRUVE0bdqUwMBA+7cxphBVJSMjg3379rFlyxbi4+MJCPDdhTm/X+ITkTBgFnAYp7quVN1cz96uV1/WfgWo6jRV7aaq3WJiPM51WNT48RAeXnBbCARenc7P05eVJyzjI0eOHOHo0aPEx8cTGRlJUFCQJSdj3BARQkNDadasGUFBQRw4cMCnx/d7ggIm4IxBuldVd5SyzxHXcz0PbXL3Hcm3raz9fGf4cJg2DeLjQYScqCAOXV6fWctu5JP7/sfB5IMV8ram9A4fPkxUVBSBgTYczpjSEBEiIyM5etTTxSnvVYUEdTmQA4wUkRX5H0B/V5tbXdtmuF4nuZ7jPRy3ZaG25ennW8OHQ1IS5OSQvXEjb/36JNs3tSAjNYPFNy1GtVwneKacjh07Rr16nv6GMcYUFh4eTlpacbVnZVMVEhQ4cfRy82ji2t/a9bqb6/WPrudOrkuE7pxRqG15+lWY4KgEBswYigQ4l5A2fb6JNdPWVMZbm2JkZ2fb2ZMxXgoICCAnJ8e3x/Tp0cpAVRNUVdw9cMrOAR5wbTvN1WcL8AMQAlxZ+Jgi0guIxZkt4pt871WmfhWtZY+WdL/3xLCrb55ZwMFNeyrr7Y0bds/JGO9UxM+M3xNUOTzren5ORNrkbhSRxsArrpcTVbVwSi9rvwrV56k+NGofxann/cjo/5tE0itj7FKfMaZWq7YJSlUXAK/izPrwq4gsFpH3gY04UyUtxJn81Sf9KlpwWDDXTKnD4JsXUafucU7p/CHrpr9V2WEYY0yVUW0TFICq3gYMx7ls1wu4CPgDZ7LXIarqdnRlWftVtOg+N3Ag9WQAAoNyiDl4Pwc37fJHKMYY43dVOkGp6ijXvacXPbR5S1V7qmoDVa2rql1VdUpJl+jK2q9CBQRS/4p3ycwMBqBx7E5Spo5Gc+xSn6laEhISEBFWrFjh71DcSkpKQkRISEgosi839qSkpHK9xxNPPIGI8MQTT5TrOKZ4VTpB1UZB0R042uyxvNcnn/IRidPn+DEiY4w3VqxYgYjQu3dvf4dS7VmCqoIiL3iUA6mdAQgIVGKOPMiBP0s7htkY48nSpUtJTEykRYsW5TrO7bffTmJiIrffbsvHVRRLUFWRBFB/yLtkZoQAENN8N1tfu8ku9RnjAyeddBIdOnQgODi4XMeJjo6mQ4cOREdH+ygyU5glqCoqqFE7jrZ4Iu91p86fsG7aTP8FZEwJRo0ahYgwe/ZsfvvtN4YMGUJMTAz16tXjnHPOYfny5XltP/roI3r16kVERAQNGjTgsssuY+PGjUWOmf9y2dGjRxk3bhytW7cmNDSUli1bcscdd7Bv3z6v4vR0D0pVmT9/PhdffDGNGzcmJCSEFi1a0K9fPyZPLljc6+4eVO/evenTpw8AK1euRETyHnbJz3uWoKqwyPMfYl+qsyxVwDdKh4fHoAEBkJBgCxyaKmv16tWceeaZbNiwgX79+tG+fXu++uorLrroIlatWsW//vUvBg0ahKpy0UUXERUVxeLFiznvvPOKTTYZGRl5SeLkk09m4MCBpKenM3nyZHr06MGuXeWvds3IyGDw4MFcffXVfPbZZ7Rr146hQ4fSoUMH1q5dyx133FHiMfr3789FF10EQJMmTRg5cmTeo3///iX0NoVVieU2TDEkgIih75I1rj2Bs7MJzHQVGCYnw5gxztfDh/svvlrqSXnS3yGU2eP6eIW/x5QpU3jppZe4995787Y99NBDPP/889x0003s3LmTFStWcO655wKQnp7OhRdeyKpVq3jllVd47LHHihzzm2++oV27dqxfvz7v3tGRI0e4/PLLWbp0KXfccQfz588vV9wPPPAAH374Ie3atWPRokV06NAhb192djYff/xxiccYN24c3bt3Z8mSJXTo0IHZs2eXK6bazs6gqrigqJOQDxsimYV22Cq8porq0aNHgeQEzi9ugA0bNjB27Ni85ARQp04d7rnnHoAClwELe+mllwoUNtSvX5+pU6cSGBjIe++9x5YtW8oc8+7du3n11VcJCAjg/fffL5CcAAIDA7nsssvKfHxTNpagqoHA3cVcYy9mdV5j/MndpayGDRvSqFGjYve3bdsWgO3bt7s9ZmRkJAMGDCiyvU2bNnTv3p2cnBy++OKLMse8bNkyMjMz6dGjB506dSrzcYxv2SW+6qCYVXhzmrewvzD8oDIuk1VnsbGxbrfXq1ePffv2ud2fu7xJenq6277uBtzm3/fVV1+xdetW74N1SXb9fBU+czL+Zb/fqgM3q/BqCGw9vRE5Wf6b+MIYd0pa8tuXS4LnZzPQ1zyWoKqDfKvwqgjaCOQmiLv2Z3594Rl/R2dMhfM0LVHuvubNm5f5+PHxzhqm69evL/MxjO9ZgqouXKvwSnYW+589C3o6m9vHTGDTwiX+jc2YCnbw4EH+85//FNm+adMmvv32W0SE8847r8zH79u3L8HBwXz99dckJiaWJ1RCQpwB9llZWeU6jrEEVf1IAFEjPib1iDN6vU74ceptHMWBjWW//m5MdXDfffexY8eJKb9SU1O59dZbyc7O5vLLLycuLq7Mx27cuDG33HILOTk5DBkyhA0bNhTYn52dzeLFi0t1rNxKwz/++MOSVDlZkUQ1JGGNCLxgIVmrehMUnEXjFjvZMHso9R79guDwEH+HZ4zP9ejRg+zsbNq1a0ffvn0JCQlh5cqV7Nmzh5NOOokpU6aU+z1eeOEF/vzzT/7zn//QqVMnevToQWxsLLt37+bXX39l9+7dpVpEND4+ntNPP50ff/yRzp0707VrV0JDQ2nfvj0PPPBAueOsTbw6gxKRSBEZLCJPishUEXlHRF51vR4kIpEVFagpKKx1T460mJD3ul2n71g38U4/RmRMxQkJCWHZsmXcfPPN/PLLL3z44YeEhIQwduxYvv32W5o2bVru9wgNDWXx4sW8+eabnHfeeaxdu5YFCxbw+++/07lzZ6+S4Pvvv89VV13F/v37efvtt3n99ddLNdDXFCQl/UUgIoHAFcBtwLlAbqlM/pIZzff8Bc7S6R/4a+G/ytStWzddvXq1395/z+zLiAlxLj1kZwWwMWsmHW4Y6bd4aoLExEQ6duzo7zAMzlx8ffr0oVevXlV27Slzgjc/OyKyRlW7eWrj8RKfiFwLTARicRLSXuBbYB2wHzgMNAAa4SyX3h3ojbNK7RYRGaeq75QqWlMm0cPf4eCMk4mM2ExgUA7ND9/Jru+70uSMk/0dmjHGlEuxCUpEvsJJOHuAfwBvqOrPJR1QRE4DRgHXAvNE5A5V7embcE1hEhxO+JD/kL64C3XC0mgQdZgfXn6CBlPeJKxhmL/DM8aYMvN0D+ok4F4gTlXvLU1yAlDVn1T1bqAlcJ/rOKYChTTuwPGTZ5BxPJiFUwez8u1TWDhyoa0fZYyp1jwmKFX9h6pmlOXAqpqhqpOA1mULzXgj4qxhbK6/lJ9XnQbAhsUb+PK5L/0clTHl07t3b1TV7j/VUsUmKFU96os3UNVjvjiOKVn7oefS/d7uea+X/205m5dt9mNExhhTdjZQt4Y5f+L5xJ3jDFjsLGtofllHW+TQGFMtlTpBiUiUiHQRkahC25uJyGwR+VFEPhCRU30fpimtwOBAhv57KL1i/8eggA8JPXocUT2xyKElKWNMNeHNGdTDwPc4xQ8AiEgI8CVwHXAqMAhYLiIt3B7BVIr6zetz7vGvbZFDY0y15k2C6gNsLlTNdzXQClgJ9AemAJHA7T6L0JRJ4N5D7nfYIofGmGrCmwQVC/xRaNsAnNkjblLVT1X1DmAzcLGP4jNlVczEmZkNG1ZyIMYYUzbeJKiGODNJ5NcD2KCqm/Jt+5F8lwGNn7hZ5JAQ0MFHSFr8iX9iMsYYL3iToNJwpjQCQERa4pxVfVWo3XEgtPyhmXIptMhhTpTATRDSL5PIzcPY9b9Sjbs2xhi/8SZB/Q6ck6+KbxgnJofNLxbY5YPYTHnlLnKYk0P6NyvJ6OYsxREZfRCWX8rBTbaGlDGm6vImQb0J1AX+JyLzgaeAVGBRbgMRCQW6ALZuchUT3u5c0jrOJifHmYS+ScttJE36K2kH0vwcmTHGuOdNgnoVeAtn6qKhQAYwWlXzl4sNxEliK30WofGZiO7Xsj/qGQA2/tyGT6Z359+X/5us47bqp/FOQkICIlLiw19TFOW+f2G9e/f2SVyzZ89GRBg1alS5jmM8K/WKuqqaA4wQkceAJsA6VT1cqNkm4EqK3pcyVUT0JY+Q/FZd3nlpHznZgSSvTGbR9Yu4Yu4VSEDRH2hjPLnooos8Lhboi4UEq5OkpCRatWpFfHw8SUlJ/g6n2vO03MY1wMeqeiT/dlXdjFNKXoSq/gD84NMIjc/FD7uLPlu+ZOm4pQCsfXstkQmR9JvQz8+Rmepm3Lhx9O7d299hlNqcOXM4duwYccUMwyityy+/nO7duxMREeGjyIw7ns6g3gKOi8hS4APgQ1XdUzlhmYrW88GeHEw6yJqpawA49u0/WD9jPe1vus3PkRlTccqbmHJFRERYcqoEnu5BPQT8hDPodhqwXUSWi8idIuKbf2XjNyLCJf+6hLaXtqHPlUsZeNNi4rmXpA8X+zs0U8OoKhdffDEiwpgxY4rsz8nJoV+/fogIt99+YhKapKQkRISEhASysrKYOHEiHTt2pE6dOjRp0oSRI0eS4uXMKCXdg1qyZAlXXHEFzZs3JyQkhKZNm9KzZ0+ee+450tJOFBS5uwc1atQoWrVqBUBycnKBe3EJCQlexWkcnpbbeEFVe+CUjd+JU07eE5gEbBaR1SLysIh0qJxQja8FBAUw9I0+dOn3KwB1fjxO3KhBNvu58SkR4c0336RFixZMnz6dt99+u8D+p556imXLlnH66afz0ksvuT3G1VdfzeOPP05cXByDBw8mJCSEOXPmcMYZZ7B+ffmLhlWVW2+9lf79+/PBBx/QokULhgwZwqmnnsqWLVsYN24cu3Z5Hj1zzjnnMGTIEADq1q3LyJEj8x5Dhw4td4y1UYlFEqq6A2eOvSki0hBnQtgrgPNxSsqfEZENwHvAQlVdXYHxGh8LadScrH6fkPncuQTPySIgw7UKb+7s5+CMpzIF/fIErH2ydG1PGg1nTSu47bsx8Of00vU/+XHo/ETBbSsGwvaPStf/zNegTdEzl8oUHR3NO++8Q+/evbn55pvp1q0bbdu2Zfny5Tz99NPUr1+f+fPnExpadIx/cnIyaWlp/Pjjj/zlL38BICMjgxtvvJG5c+dy3XXX8b///a9c8U2aNImpU6fSpEkTFi5cSPfuJ9ZVy10wsWEJ04TddNNNnH/++bz33ntER0cze/bscsVkvFwPSlUPqOpsVb0MiAGuAd4FmgOPAN+JSLKIvCwivcRdnaepcsLbdCfgo0hn4EB+Nvu5KUGfPn2KLTGPjIws0Pacc87hqaee4siRI1x11VWkpKQwbNgwcnJymD59Om3atCn2fR577LG85AQQEhLC5MmTiYiI4Pvvv+err8peOJyVlcWECRMA59Jd/uQEzhlgnz597J6TH5S6zLww14q784H5rmU3LsA5sxoI3IVzWfBx4BkfxGkqWOCufW63a0oK9leGKY6nMvPwwnNBAg8//DBffPEFS5YsoXPnzhw6dIibb76Zq6++2uP7jBgxosi2iIgIBgwYwLx581ixYgU9e/Ys0/ewevVq9u7dS2xsLP379y/TMUzFKHOCyk9VM4CPgY9FJADoBVwO7PbF8U0liItzLusVkt0gkOO79lO3SZSbTrVY5yeKXnbzxlnTil7280bvqlHM4m2Zee79qNatW3Po0CH+8pe/MGnSJI99IiMji5yN5cotPti6tezTdiW7/t+3b9++zMcwFcPnS76rao6qLlfVO1W1HD+BplIVM/t50PAs9k7rRVpx60sZ46WFCxeSmpoKOIll27Zt5T6m3U2omcqUoESkqWv597OLe/g6UFPB8s1+jgjZ0fXgJqAnNI/9nU+ue97m7TPltnbtWu666y5CQkIYPnw4hw8f5uqrryYjo/AN0BMOHjzIoUPu/0DKna2hefPmZY4pPj4ewCfVgMa3vEpQInKliCQC23CWf19VzKPwDOemOnDNfk5ODoG7D7Or+w1kpAcz74Xh/PrfEOZeNJf0Q+n+jtJUU0ePHuWqq64iLS2N5557jjlz5tCnTx/WrFnDAw884LHvPDdDHg4dOsRHHzmVjOWZzaJr165ER0ezdetWlixZUubjgFO8AU7hhSm/Uico19RH7wDtgcPAL8DXxTy+8XmkpnKJ0GTEDDawgOREZ/Dh9u+3M+/ieRw/ctzPwZnqaOzYsSQmJnLZZZdx9913ExAQwLx582jcuDH//Oc/WbhwYbF9n3rqKRITE/NeZ2Zmctddd3Ho0CG6du3KOeecU+a4goODefjhhwG4/vrri5Ss55aZF3cWl19MTAwhISHs2rWLAwcOlDkm4/CmSOIR1/NdwKuqan8i1HQinHzDZaRnNOfjWz8GYOs3W/lgyKtc8f7NhNQL83OAxt8mTpzocbzPsGHDuPDCC5kzZw5vvPEGLVu2ZNasWXn7mzVrxptvvkn//v254YYbOP300/MuueWKi4uja9eunHbaafTt25eIiAi++eYbUlJSiI6OZs6cOeX+Pu655x4SExOZMWMG3bt3p1u3brRp04b9+/ezbt06tmzZwubNm0ssNQ8ODubSSy/lgw8+4PTTT6dnz56EhYURHR3NxIkTyx1nraOqpXrgrKi7qrTta8uja9euWht896/v9Ame0Ekxd+qBSRG68fGempGa7u+wKsS6dev8HUKVFx8frzgLlnp8vPzyy5qYmKh169bVoKAg/fLLL90eb9y4cQpo9+7dNSMjQ1VVN2/erIDGx8drZmamPv3009quXTsNDQ3VmJgYHTFihG7evNnt8XLfv7BevXopoMuXL3fbb/HixXrppZdqTEyMBgcHa5MmTfScc87R559/XtPS0vLazZo1SwEdOXJkkWPs3btXb7zxRo2NjdWgoKC876E28OZnB1itJfx+FaddyURkG7BSVYf5Lj1Wf926ddPVq2vH5Bmr//ERbWU4EdHOKisbN/Qi4cElBIcXHf1fnSUmJtKxY0d/h1Hr2dIV1Y83PzsiskZVu3lq402RxKfAGV60NzVMtzsvIa3OeXmv27ZbSdLzl5CVlunHqIwxNZU3CepxoIGIPCcigRUVkKnCJICmoxex88iJ0fZt9yxDmzWwCWaNMT7nzYq6KSJyDrAIuMK1TtRWIKeY9hN8E6KpUiSApqM/Yudr/Wm69nOYAcEZrtJzm2DWGONDpU5QrolfbwfaAoHASTg3Ios0dW23BFVTBQTSZMwnZMXUIyijUMl57gSzlqBMOSUkJFDae+SmZvKmzHwccAeQhTPv3h9AakUEZao+CQwi8KD70f82wawxxhe8SVA3AseAc1T1pwqKx1QjUswEs8fDw8jalUq9JvX8EJUxpqbwpkiiBfCFJSeTx80EszlBwsdHL2TWObM4sNlG0htjys6bBLUN5wzKGEehCWazGkWwSAexltPY/8d+Zvacye61tuKKMaZsvElQ/wZ6iUjdigrGVEP5JpgN2nuQjgueJTDUGYWQuiOVWb1eZ9uq7/0bozGmWvImQT0NbAA+FJGTKigeU811GNyBEf8dQUj9EEA57+IPifytDykfL/J3aMaYasabIokPcSr4+gCJIrKJ4sdBqape5IP4TDWU0DuBkctHsvH50fS45FsAAndczaZ3XqX1Ndf7OTpjTHXhTYI6v1C/dq6HOzZ4oZZr3rU54Q/eRvpPy6kTlkad8OPEpY9m4/T9tB19n7/DM8ZUA94kqAsqLApTI0V2vZTUeks5uqI/desfJigkm5MCH2D9pL20u2uCLdNtjPHIm6mOllZkIKZmqte+B2lhX3Pooz5ERO4hIFBpv3EimVGTCTp01BlLNX68zTxhjCnCqyXfjSmLsLhO1Lnye/bvi4OvcObvO5iKqJ6Yv88mma0WRMTrx6hRo8r0Xt26dUNE8MVyNqmpqYgI9eoVHTweHR2NiLB3795yv48/rV27FhHh5JNP9ncoPuPNJT5jyiw0Jp7AUWvIjm9GYEahxZht/r5qY+TIkUW27dy5kyVLllC3bl2GDh1aZH95lmM3tVuxCUpEvgDGqerXZT24iPQEnlXV80psbGq8oPrR6MFst/ts/r7qwd3y7itWrGDJkiVER0d7XP7dW++99x5paWkkJCT47JimevF0ia8DsEpEPhORq0WkVMumikioiFwrIp8DX1B8pZ+phSQuzu32w0SwaemmSo7GVGXx8fF06NCBOnXq+DsU4yeeElRb4F9AL+AtYJeIfCwifxORISLSW0S6uJ6HiMhjIvIfYDcwFzgX+AeWoEx+bubvyyCYz7Uv8/rPY820NX4KzM/mzXMWfKzhCz++8cYb9OrVi4YNGxIcHExMTAynnnoqd955JykpKQXaeroHdfz4cf7+97/TrVs36tevT3h4OCeffDKPPfYYhw4d8kms2dnZ3HbbbYgInTt3ZuvWrQX2//TTTwwbNowWLVoQEhJC48aNGThwIMuWLStyrFNOOQURYenS4mvNbrnlFkSEJ598Mm/b0aNHefrppzn11FOpW7cuoaGhNG/enJ49e/L444+TlZVV7PHy27dvH2effTYiwrXXXsvx48d58MEHERHuv//+Yvu99dZbiAh9+/Yt1fv4nKp6fACtgZeB/TiDcrM9PHKAvcDzQEJJx64Jj65du6rx0ty5qvHxqiKa2bSFfhQ5TJ/gCX2CJ3RWh1GaOHG0Zmdm+S28devWVe4bzp2rGh6uCice4eHO9mpg+fLlCmh8fLzHdvfdd58CGhISon369NFrr71W+/fvr+3atVNAFy9eXKB9165dFdDvv/++wPYjR45o9+7dFdD69evrZZddpkOHDtWYmBgFtE2bNrply5YifQCtW7dukbgaNWqkgO7Zsydv29GjR3XgwIEKaL9+/fTQoUMF+rz99tsaHBysgJ566ql67bXX6tlnn60iooBOnDixQPsXXnhBAR0xYoTbzyY9PV0jIyNVRHTTpk2qqpqZmZn3fUZFRemll16q1157rfbp00ebNWumgB45ciTvGL/++qsC2qlTpwLH/uOPP7Rt27YK6AMPPKA5OTmqqpqUlKSBgYEaFRWlaWlpbuPq2bOnArpgwQK3+wvz5mcHWK0l5Z+SGuQ1hDCgPzAR+C/wA/AnsAb4BBiPM5g3tLTHrAkPS1Dld2jrIX2ty2v6jyZ36NGpYarz0I1P9tb0g4f9Ek+lJ6j4+ILJKfdRwi/8qqI0CergwYMaFBSkUVFRunnz5iL7161bpykpKQW2FZegbr311rzEsHPnzrztR44c0YsvvlgBPf/88wv08SZB7dq1S88444y8hJKRkVGg/ebNmzUsLEwBnTp1aoF9//nPfzQkJERFRL/44ou87Tt37tSgoCANDw/Xw4eL/r/+97//rYD27t07b9vHH3+sgPbs2bNIAsnOztYVK1YUiM1dgvruu+80JiZGAwICdPLkyUXed/DgwQrorFmziuz7+eefFdAWLVpoZmZmkf3u+C1B2cMSVEU6nnpc/3ymt+o88h5bn22rBzf+WemxVHqCEnGfoEQqN44yKk2C2rRpU94v29Jyl6AOHDigoaGhCujXX39dpM+WLVvy9v/0009520uboNavX6+tW7dWQB955BG3cT344IMK6IUXXuh2/2233aaADh48uMD2AQMGKKCvv/56kT65iXX27Nl522bOnKmAPvroo27fp7DCCWrRokUaHh6uYWFhunDhQrd9li5dqoCeccYZRfaNGTNGAX3yySdL9f6qlqCq3MMSlO/kZKTp9skXFEhS+ydF646Vyys1DjuD8k5pElR2drY2adJERUQfeeQR3bhxY4nHdZegPvnkEwW0Xbt2xfbLTQSTJk3K21aaBLVw4UJt1KiRBgYG6muvvVbs8c866ywF9K233nK7f/Xq1QpoZGRkge0LFixQQM8999wC23fs2KGBgYFar149TU1Nzdv+yy+/qIhoZGSkTp8+vcAlSHfyJ1LcUqUAACAASURBVKgpU6ZoQECAxsTE6LfffuuxX6dOnYp8zocOHdK6detqcHCwbt++3WP//HydoGygrqkyJLgOzW5bwvbsO/K2NYzZS+SGi9n89nQ/RlbB3BSOEB7ubK8hAgICePPNN2nYsCETJkygbdu2NG3alMsvv5zXXnuN1NTUUh1n27ZtALRq1arYNieddFKBtqU1dOhQ9u3bxz//+U/GjBlT5hhy3//gwYMcO3ZiCb2BAwfSqFEjvvzySzZtOlGxOnfuXLKzsxk6dCh1655YzeiUU07h2WefJTU1ldGjRxMTE0Pbtm0ZNWoUixYtIifH3TzdsH79esaOHZtXlHHWWWd5/L7vuMP5eXvllVfytr3xxhscPXqUK664gmbNmnnsX5EsQZmqRYTm1/2T3Y2mknE8GIA64enEfzGGzKgGaE2sciu08CPx8c7rGjZw+YILLiA5OZm33nqLm2++mZiYGBYtWsQtt9xC27ZtSUxMLPEYzh/eeJzHMbeNt/76178CMGHCBNavX1+uGNwJCQlh2LBhqCpvvPFG3vbcr93NuPHQQw+RlJTE5MmTueaaa0hPT+eNN95g8ODB9OzZk7S0tCJ94uLi6NWrF9nZ2dx5550lJv8RI0YQGRnJO++8w4EDzirYr776KgC33XabV9+jz5V0imUPu8TnLwd/XqZHXo1QvQ3VkEKXvyqwyq3SL/FVc6Wt4nNny5YtOmjQIAX0ggsuKLCvrJf4cqvvvL3Et2fPHv2///s/BbRx48b6yy+/uD1+SZf41qxZ4/YSX/59CQkJmpOTk3c5sFWrVnnVdSX5/vvv8yofx48fn7c9/yW+Y8eO6UUXXaSA9ujRQw8ePOjxmPfee68C+tJLL+Xdlzr55JNLFU9+donP1BoRnfsQNHgN2W8HQkahnbnTI5lqLTY2Nm/cz88//1xi++7duxMaGsqGDRv47rvviuzfvn07n332GQC9e/f2Op4nn3ySiRMnsnv3bvr06cOaNUXH5fXq1QuAOXPmuD3GrFmzin3/Ll260LlzZ5KSkli5cmXe2dPIkSNLfUbWrVu3vDOb4j6zsLAwPvzwQwYPHsw333xD37592bdvX7HHHDt2LAEBAUydOpUpU6bkbfO7kjKYPewMyt9yiqlyy6mgKjc7g/JOac6g1q9fr7Nnzy4wbifX008/7baSrLgy81tuuUUB7dKli+7evTtve2pqal6BRHnKzFVV//Wvf6mIaERERJFqwfxl5jNmzCiwb8mSJRoaGlqkzDy/v//97wrosGHDNDo6usDYp/w++eQT/fTTTzUrq+CYwIyMDO3fv3/euKZc7srMMzMz9ZprrlFATznllAJl+YXlfnaANmjQwO2/VUmsiq+KPSxBVYJiqtwOEKGf3PWJZmX4dlCvJSjvlCZBrVq1SgENDQ3V7t276zXXXKNDhw7VDh065G3//PPPC/Qp7UDdQYMG6dChQ7Vx48aKjwbqqqrOmDFDAwICtF69erp8+fIC+9566y0NCgpSQE877TQdNmyY9uzZM2+g7rPPPlvsZ7F79+68Qb4UGvuUX27ibtiwofbr10+HDRumgwYN0iZNmiigLVu2LFBhV9xA3ezsbL3++usV0Pbt2+vWrVvdvt+nn36aF9Ptt99ebPyeWIKqYg9LUJXAzUwLxwnWBVyhz4Y/pOse7q2Hk5J99naWoLxTmgS1f/9+ffHFF3XgwIHaunVrrVu3rtavX187duyoY8eO1fXr1xfpU1yCUlVNS0vTF198Ubt06aJ169bVOnXqaMeOHfXRRx/VAwcOFGlflgSleiIRhYWF6X//+98C+3744Qe95pprtGnTphoUFKSNGjXSAQMGFEm07uTed6PQ2Kf8EhMT9W9/+5v26tVLY2NjNTQ0VKOjo7VLly46YcIE3bdvX4H2xSUoVdWcnBwdO3asAtq6dWu3g6XT09PzxpCV9WfA1wlKnHYlE5GbgHmqWrRspBbr1q2b+mK9GlOCefOce04pKeTExvJd40F8+kMUV9/zbzp0Xc/BvVEc6/gmzS+4pNxvlZiYSMeOHX0QtDHVx7x58xgxYgR9+/b1OGegJ9787IjIGlXt5qmNN0US04CtIvKSiLT1op8x5Td8OCQlQU4OASkpdP/+nwx5uREdujrlwJHR+4nZNoiNUx6htH90GWMcx48fZ8KECQDce++9fo7mBG8S1EdAA+AeIFFE/isiA8XbwQDG+ICIcPJdd7IrcjLH052VYIJDsmjb8Fk2jT+f9P0H/RyhMVXf1KlTGTVqFJ07d2bdunX06dOHSy+91N9h5Sl1glLVy3BmNp+IM2P5hcBCYLOIjBORmIoJ0ZjiNblkLJnnfsWBfc3ztp3UehnHb29NdrNmNX75CmPK4/PPP+eNN95g3759DB8+nPnz5/s7pAK8GgelqltU9RGgJXAd8C0QhzOTeYqIvCkiPbwNQkSCRaSf6/LhtyKyQ0QyRGSbiCwQkd4l9B8mIqtE5JCIpIrIahEZKyIev7+y9jNVS722XWlwUyLb9/Z2NnwFEe8dIHDnTqekIjkZxoyxJGVMIQsWLEBV2bt3L3PnziU6OtrfIRVQpl/EqpqpqvNUtSdwOvA6kAUMA74UkTUickNpV+HFWRTxc+BeIB5nCY8PcNagGgIsF5Gn3HUUkSnAPKAbsAr4DGeRxMnAAhEJ9GU/UzUFhjWg+R3L2C6PoP/GBvYaUwOU+0xBVX8GngRmAeJ6nA5MB5JE5MZSHCYHeA84T1WbqeoAVb1aVU8BrsFZDPExEemTv5OIDAFuA3YCnV39LsdZDTgRuBy4vfCblbWfqeJEaH7teNjv/raoFlqx1RhTtZUrQYnI+SLyPrAZGAukAzOBa4H/AI2BaSJyp6fjqOoyVR2qqqvc7Ps3MNv1ckSh3Q+7nh9S1Y35+uwCbnW9HOfmkl1Z+5lqQOLi3G4/pA34+sWv0Ryr8jOmOvD6F7CIRIjI3SLyO7AEGAxsBx4BYlX1JlX9t6oOBM4GjgIeE1Qp/Oh6js0XRyzQFedizruFO6jqSmAb0BToXt5+phpxs3xFBsEspR+fPfAZb18ykyPJns+mrFTdGO9UxM9MqROUiHQRkRk4v7xfwrlfsxLnHlFrVX1OVffn76Oq3wEf4xRSlEfuuKsd+bad7nr+zcPg4e8LtS1PP1NdFFq+Irt5LF+2+itr6QxA22avoR91JuWDN912DwoKIiOj8E0sY4wnmZmZBAb69ta9N2dQq4EbXF/PwLl301dVP1BV9ytnOY4CQWUNUESaAqNcL9/Ltyt3tbBkD91z/0zOv7JYWfuZ6iTfwN7AbVvotf41znn4HNp1Wc8ZF3xPg4aHaHn0r/zx3NVkph4r0DUiIoJ9+/bZWZQxXjh8+DD169f36TG9SVBJwAM4l/FuVtW1pew3Ggj2NjAAEQkC5gIRwFJVXZxvdz3X81EPh8hdqSv/p1bWfvnjGuMqSV+9Z88eD4cxVUVgcCD9JvSj7/gLSDvqXP6TAGjTcj77X/0Le78/sXRDVFQUx48fZ+vWrRw5coTs7GxLVsa4oapkZGSwd+9eDhw4QFRUlE+P782ZzUlahp9SV59sb/u5TAX6AVsoWiCRW6rlbUxl7ZdHVafhTP1Et27d7DdXNdKk//WkbT+PXfMH06Sx8zdWkxbJHP/lPHa+eglNlv5A0JYtxLduzYGXX+bAqaeyffv2YpfXNqa2CwwMpH79+sTFxREaWtqRRaXjTYJaIiL/VdW/e2okIvcAF6vqheUJTET+AdyIUwreT1V3FmpyxPVcj+Ll7juSb1tZ+5kaIqz5SdS58ye2zrmPZgGTCQzKJvSHDJrOW5g3firgzz9pdM01NKqBS68bU114c4nvfODkUrT7C85ZT5mJyEs4lX97cJLTRjfNklzP8R4O1bJQ2/L0MzWIBAQSO2oSh09ZwsH9TWA+NrjXmCqmIsb5hOAMvC0TEXkeZ0aJfcAFqrqumKa5peedRCSsmDZnFGpbnn6mBmp4ej/qXb8e3et+vw3uNcZ/fJqgXDObd8WZTLYs/SfiFGIcwElOPxfXVlW3AD/gJMQr3RyrF864qZ3AN+XtZ2quoLoRSLz7E+rDEkHiB4mVHJExBkpIUCLyae7DtenC/NsKPZbhlGd3BL7wNhAReRp4CDiIk5xKc/byrOv5ORFpk+9YjYFXXC8nuimDL2s/U1MVM7j385y+zL9iPu8Pe4dju3b5KThjaiePK+qKSP5f0MqJCjhPfgEGqaqncUaF3+cyYJHr5Wrgt2Ka/q6qEwv1fQVneqJ0nAlnM3HugTXAWQ5kqKoWqSIsa7/CbEXdGiTfqr2Z0c34NLM3qw+2A6DXFcvp0udnDjR9kfihpZle0hjjSWlW1C0pQeUWOwjwKc7URi8W0zwD2Kaqm8oQ6CicyWZLslJVe7vpPwxnLsBTgEDgd5w5AV/1dBZU1n75WYKqudIOpLHkniXsXPFfRj81ncAg57/EpuQLaXrDHMKbNPFzhMZUX+VOUIUOtgr4uPAZTG1nCarm27ZwOlF77ias7okZJw7tb8ihZi8RN+R6P0ZmTPVVmgTlzYq651pyMrVRi8GjYeBaduw6M29bRNQB4o7fwK6rTiGnZUtbudeYCmDLSRhTCmFNW9Hs7m/ZEfZi3lRJfAVNFq0lYOtWW7nXmApQ7CU+EXnE9eWrqnog3+tSUdUJ5Q2uOrBLfLVP2s7NHHz7GppN+J/7ARXx8c5EtcaYYpXrHpSrgk+Bjqq6Id/rEt8XZwq+WrFkuiWoWkoVDQxA3PxEKIJmZhEQZBcojClOaRKUp7n4JuAkpL2FXhtjRJC4eOeyXiGHaMD8s2Yw4LUBNO/W3A/BGVMzlLqKz7hnZ1C12Lx5zj2nYyeq+zIIZjEDWUtnJFAZ8uRe2tz2DKENI/0YqDFVj0+r+IwxhRRauVfj4th05Th+r9MFgNPPW0OnVlNIe7MNKQte93OwxlQ/dgZVTnYGZQrb/+d+lt7zDgMG3kdY3fS87UnJPYm6chYN2rT1Y3TGVA0+PYMSkVtFJENELvXQZoCrzU3eBGpMTRJ1UhRDF97MvroPkn7sxIT5CfFfUWfiX8hsFIHauCljSuTNJb4rgP3AJx7afOJqM7Q8QRlT3UlAILHDnkQHrGPrrl7Oxq8g5M0sgvcfRmzclDEl8iZBdQB+9TRHnWty1V9xFi00ptYLa5pA7D0r2N14LjlvB7hdFDHnoYf8EpsxVZ03CSoGKM16A7uBxmULx5iaqfH5w5GDxYw53LaNr1/6muyMEifPN6ZW8SZBHeLEUuietABSyxaOMTWXxMW53X6YBnx2/2e82vlV/ljyRyVHZUzV5U2C+hHoLiInFdfAte9s4KfyBmZMjeNmUcSswCA+53wA9q3fx7z+81h0zSscTPzVHxEaU6V4k6BmA8HAQhEpUifrWpl2Ic66SrN9EZwxNUqhcVPExxMwcybN/34/oQ1C85q1aTKF8G+78uekG8k4fNiPARvjX96sByXAYuASIAv4EmeBP4D2wLk4Uyf9V1Uv8X2oVZONgzK+kLorlaWPLOXgt+8z8tE5edsP7Ysi4/dLif5oJbJlC8TFOWdiw4f7MVpjys+nCxa6DhgCvAyMpug8flnAdOBeVT3uZazVliUo40u7v1hM4I+30Shmq7PhK2AGBav/wsOdMzFLUqYa83mCynfgpkA/IN61KRlYqqo7vT5YNWcJyviaZmWy9a3HiT4+ibBH0mxJD1MjVViCMidYgjIVJX3PNkIbxyJu9imQceAAoZE2Ca2pnmyyWGOqsToxLZD4eLf7JBrev/hR1kxbQ05WsWPnjanWvE5QItJeRKaIyG8ictD1+E1EJotIh4oI0phay01pOiGwrVtTNnzbmI9u/ohXO7/Kho83YFdDTE3jVYISkVE4Y5xuAToCDVyPjsBtwE8iMtLHMRpTe7lZ0mPfDZfx6c4TP2Z7E/fy9oC3+fyvD7Lnm+V+DNYY3/JmNvMzcKr0QoAPgAE4iekvwKXAezjjpKa72hpjfGH4cKcgIicHSU6m0auLGPHV0/Qd35eQ+iEABIdm0L37VBr92Y/NEy/k2N9fdGZLt1nTTTXmzTio+cAQYISqvl1Mm2uBecC7qnq1z6KswqxIwvjT0d1HWfHkCurufIneQ1xnT1+BzgCx0nRThfm0ik9EtgNbVfXMEtp9B8SparNSR1qNWYIyVcGBH1aSvvxOmjX7Be7CStNNlefrKr5GwIZStNsIRHlxXGNMOTXs0otm9/3Mrug5qLvkBGhyMtnptWYMvakBvElQB4BiJ4rNp7WrrTGmkjW58DpnOiQ3JBoOv9aS9XPeIyfbStNN1edNgvoaOFNEBhXXQEQGAt1xJmgxxviBTJjgtjSdqyAk9Cjv37yGqZ2nkvh+opWmmyrNmwT1d9fzuyIyU0R6iUiciLR0ff06sADIydfWGFPZCpemt4xlz/A+pJ9ehy8/PIeM9FD2rNvD/CHzmX7GdP5css4SlamSSp2gVPVL4G5AgJHAMmAzkOT6+nrX8e5WVTuDMsaf8pemp2whZuYyZPAmwrrfl1eaDrBjzQ7SPr6KA9e2JKtJYytLN1WKVwN1VfVfwJnAXCAFZwbzbNfXc4AzVXWyr4M0xpRfaHQzzvu/C7lr0130uL8HQXWCaBq/g5OzfyPqg20E7d4DqpCcDGPGWJIyfmeTxZaTlZmb6urwtsOkzHyATv+Yhuwruj+7SWMCd+6q/MBMrWCTxRpjitWgRQNOfuw12O9uvnQI3LWblIlnsOuLJZUcmTEOS1DG1HJSTFk60RAXt5omW/uz5blu7PxiWeUGZmq9wqvi5hGRaeU4rqrqzeXob4ypLOPHO/ecjh3L25QTHEDAVSfGSrVsuYaZo+cRkrCdXo/3ouXZLf0Rqallik1QwE3lOK4ClqCMqQ5y5+d79FFISYG4OALGj2d/xxakLX+YFs2+ZdPaVmzZEAcb/uTPT/+k9fmtubjHfqLnvJzXh/Hjba4/41OeEtToSovCGONfw4cXSS5RAF2+4cBPX7B51Q9IwGE0xymqCv98IZHLFkFOttM4t/Iv91jG+IBV8ZWTVfGZ2mLfhn2sGr+KX+b+wp05fyeSQ0XaaFwckpzsh+hMdWNVfMYYn2nUrhGD3xjM2N/HEuEmOQGQkkLyWy+Rk5FZucGZGqlMCUpE6olIbxG5UkTO8nVQxpiqq1HbRkh8vNt9Eg3x3M+hqbFsfv3/yMpXeGGMt7xd8r2+q7pvL7AUeId8xRAicquIpIiIxzWjjDHV3PjxRSakVdeEtAANo3fTKuxp0mY354/J4zh+2Jb5MN7zZsn3cGAFTnXfYeAznHn58vsUiAUu91F8xpiqqNCEtMTHk/HSyyR3GE76sTp5zepHHmLLsu+YFD+JZX9bRvrUmbYUvSk1T1V8hd0HnA68DYxR1aMiUmBRGVX9U0Q2An19GKMxpioqVPkXCsQDGQdeJnn+40Qzl6Cg4/zv0zNJP5rOgfFTCGAx4Lo/ZZV/pgTeXOK7CtgB3KiqRz20SwZalCsqY0y1FdIwhvibXyH0uu0kBbxGePNYAPqxlBAKFU8cO+aMvzLGDW8S1EnA/1Q1vYR2e4HosodkjKkJgsLr0f76UYxNHMvQ+UOLrfzT5GS2vj8Nzc6u5AhNVedNgsrEOYsvSSyQWrZwjDE1TUBgAJ2u7ORxKfrY9Js5MDmWpJmPkZl6pJIjNFWVNwlqA3C6iBSbpEQkEjgVWFvewIwxNYu7pejzV/5Fxewkoc4zZMxrzqZ/3sDRrTbgt7bzJkG9BzQBJnho8wxQD3i3PEEZY2ogN5V/6c9OJLntFRxPP7HKb936qbSOnkXIZ2359cLRZDdvaVV/tVSppzoSkbrAaqAd8CVOwpoELMcZD3Ul0A/4DThDVWvFwAeb6siY8ju+fxc7FzxNVPZc6kc496r2vBtFxMIjBQorNDwcmTbNqv5qgNJMdeTVXHwi0hInMXXDmbFcXM+4vv4JGKSqW8oUcTVkCcoY38nJOM62914mbNcrhD20j7oZRWeiyGrUgJw/kwmJiPRDhMZXfJ6g8h14AHAJ0BoIBLYAnwDvqWqOp741jSUoY3xPc3IgKAgp5vdT2mvh7Dw6gKhL/kZE+1MqOTrjC6VJUN4M1M2jqh8BH5UpKmOMKYEEBDhVf+5mRo+GsHrHaFVvPjnfv8uWRWcQ1Pleml54pdPP1BjF/muKyAIRuVhECk9nZIwxFc/NfH85IUEcu6Ru3uuAAKVl7P9otv8a9v0jnqRZ/0fW6zNsOqUawtOfG1fgnCVtEZHxItK2kmIyxhi3VX8BM2cT9vpBtoW/wo7tnQo0j26ylYTVTxN422jnzEv1xHRKlqSqpWLvQYnIFOBqnIU1cxt9CcwE3lVVm0cfuwdljD/t/+ELjqx4luZRSwkOyYS7cOayKcQWUqx6yrVgoaqOBZrjJKlPgRzgXJwEtVNEpovI2T6M1xhjvBLV5Tzi7/2E7EuSSDp8O+omOQGQksL6yU9ybMfWSo3PlI8346CaAyOBvwLtXZsVZ4aJmcCbqrqzIoKsyuwMypiqQ+PjkZSUItsPB9ajwZxUMjOC2L77XOqceS+Ne12K3WL3H58u+a6q21X1WVXtCPQEXgeO4CSriUCKiCwSkUEiEliewI0xpizcTaeUSTC7ezQGIDgki/jY5TTZPpC9LyeQNPMxMg4e9EeophTKVJOpqt+o6migGc5Z1Qqc8VADgPeBbb4K0BhjSs1NYYXMnEHQTTeyd1fLAk1jmqaQUOcZchY0Y/ewrmQ3a2qVf1VMmQbquj2QyAXAXCAGUFWtFWdRdonPmGpCld1ffkT6dy/TLHoVwSFZzvavgBlARr6mYeHIdJtSqSL59BJfMW9QT0RuFJFVwH9xkhM4M0sYY0zVIULjcwcSd/8ysvpvJunI7RzcFwPzKZCcACTtGMduvYddv+zyS6jGUdapjvoA1+OMlQrDmYfvOPAhTsHEp+qrU7Mqzs6gjKm+NCfbNaWSm33AUzxB614RnDvqIM0H301IpM3/5ys+nepIRFrh3G8aCcThJCVwJoidBcxV1QNljNUYYyqdBARCXLzbKZUOEQFA47D/khCyhPT3JpK0vy/hPe4mpucFVgFYCTxe4hORcBEZKSLLgY3AY0A8cBCYAnRR1S6q+i9LTsaYasndlEqhYaw/ayQBwUKXvmsAqBN2nIQWn9A45SKOjIgmq1EEakUVFarYMygReR1njae6OGdLOcDnOJfwPlDVjOL6GmNMtZFbCPHoo5CSAnFxBIwfz1nDh3PyroPs+WAXwfveJbLRHqfdV9Bgwf4T962Sk8m58QYkR5HrRvjlW6ipPE11lLtsxmZgNjBLVW0YdiF2D8qYmk9zcti19F0yfn6FlhO/QPYVbZPdMIDV979J++GXEBlv96pKUq71oETkTWCmqi6viOBqCktQxtQuGhDgdp0qBZ4OeAzVQFr1bcVp159Gxys6EhwWXPlBVgPlnYvvOktOxhhTkMTFud1+PDwUzQkEhc1LN/PBiA9YdOFNJL14CbtWLkbnzrVlQLxkq3sZY4w33BRVaFgYe68fx0kXnXSivhk49eyvSWj+CU3mXwY3XGfLgHipTCvqGmNMreWmqELGjyd2+HBGAIe3HubnOT+zccFKTjrlT6ftfJDMQsc5dgx9+GHEZqsols+mOqqt7B6UMcYdp7BiARk/v0bLB5bhbtSUAinPn09I5xtpev6VSGCtmCEOqISpjowxxrgnAQE0veAq4u5fCrGx7ttEQ3yLz4nefh2vdHiBpY8uZe/vxS1qVftYgjLGmAomEycWvW8VDFzlfJ24uiN7/zjOlxO+ZErHKUw/czo/TvmQ9Mn/qNWFFXYPyhhjKpqb+1Y88wy7W9Ynbc001q0peIa1/fvtNKz3DKFfpRQYEMyYMQWPV8PZPahysntQxpjyys7IZuMnG/llzi9s+GgDQhqP1JuA7C/aNiumEQFbtxMQElL5gfqQTyeLNcYYUzECQwLpMKgDHQZ1IG1/Gn+8uwRucd82aM8+js5qxJ7UPoR1uZHGvQYiATXzbk3N/K6MMaaaCosK45SbByPx8e4bREPd+qkkNFtMkx2DOTS5CSse+4Q96/Y496hq0D0rS1DGGFMVuZtlPTiQ9AF1Cmzbv7M+K5/5H190uoXM626oUYOBLUEZY0xVNHw4TJsG8fEgAvHxBMx6g5Bph9kRNZfk7eeTfqwOa78+BYB+LCW48CITx46R89CDfgjeN6xIopysSMIY4y9Zx1L5c8mf/PL2Boa+e1Wxg4F3vtCB4w0HEXPRLdSNTajkKN0r12zmpnQsQRljqgKNi0e2pBTdEQ38w/kyJzuAXTs7khF9OY0vvoWwpi0qNcb8bCaJUhCRYSKySkQOiUiqiKwWkbEiUus/G2NM9SHPTnBzz0rQoSdeBwTm0KzFb8SHPkPIp3Fse+E0Um5+gpy4uCpZWFGrfwmLyBRgHtANWAV8BrQDJgMLRKT2TIxljKne3N6zepNj4zaSfPQ+du1oU6B5YFAO9dZsoum0CQRs2ZJXWKGjq05hRa29xCciQ4AFwE7gPFXd6NreBFgOdATuVtV/eDqOXeIzxlQXR/5MZN/nUwk/upjGTTeTdlMoYWnHi7TLigxm26RHaHzR6Aq7DGj3oDwQkdVAV2Ckqs4ptK8XsAInebVQ1ZzijmMJyhhTHR3+/RfqdzwNoZgcMA+yswLYtasTmdGXEX3hTdRtkeCz97d7UMUQkVic5JQBvFt4v6quBLYBTYHulRudMcZUvAYdOiPx7lcHJtp5CgzKoXmLX4kPHU/Y8tbseOEvJM14iPRJf6+UAcG1MkEBp7uef1PVtGLafF+orTHG1CzuVgeuU4c9Ay9g947WBbYHBCjNWiSS8OPzhD50X6UMCK6tCaqV6znZQ5vces1WHtoYY0z15aawQmbMIGbmpzS+708On7mWpNS7BAAvOgAAD/1JREFU2bmjLXk3OuaDFBoPzLFjzkztPlZbJ4ut53o+6qFNquu5fuEdIjIGGAMQF1fMKbIxxlQHw4cXu3xHgzadaNDmZeBlUjdvZO/S6cTvfcHtgGBS3IzBKqfaegaV+/mWqUJEVaepajdV7RYTE+PDsIwxpmqq16otCTc9X/wkthXwx3ptTVBHXM/1PLTJ3XfEQxtjjKld3Ny3Ijzc2e5jtTVBJbmei/lTAICWhdoaY4xxc9+KadMqZJXf2noP6kfXcycRCSumku+MQm2NMcaAx/tWvlQrz6BUdQvwAxACXFl4v2ugbizOQN1vKjc6Y4wxUEsTlMuzrufnRCRvkioRaQy84no50dMsEsYYYypObb3Eh6ouEJFXgVuBX0XkcyAT6Ac0ABbiTBprjDHGD2ptggJQ1dtE5EtgLNALCAR+B2YCr9rZkzHG+E+tTlAAqvoW8Ja/4zDGGFNQrZ3N3FdEZA+ep0zyJBrY68NwagP7zLxjn5d37PPyTnk+r3hV9TjTgSUoPxKR1SVNN28Kss/MO/Z5ecc+L+9U9OdVm6v4jDHGVGGWoIwxxlRJlqD8a5q/A6iG7DPzjn1e3rHPyzsV+nnZPShjjDFVkp1BGWOMqZIsQRljjKmSLEH5iIgME5FVInJIRFJFZLWIjBWRUn/GIhIsIv1E5CUR+VZEdohIhohsE5EFItK7Ar+FSueLz8zDsSeIiLoe9/9/e+cebVdV3eHvR4AEwlNeEYMKEvBVCgRiDeWpqKgpaCLVFpNUynAkjooPoMFaFLQYxBdiMCDgBUzURsqz9TFKDGLByM1lxCigCUloCKRIgAABg5LZP+baZrM559xz79knZ5975zfGGitnzbXXY+acPe96zjLa22nK1pekHSSdI+luSU9KelbSKkkLJB1Vdvu3NmXqS9JYSZdK+q2k5yT9QdJySXMlHdCO9m8tJB0s6UxJ35F0v6TN6XczpcVyW9e/mUVoMQBzcO+8zwG3AjcAT6W0/wBGNFnOW9MzBjySyvo+sCyXfkGn+1slndUp+0jgT8DmVN5Zne5v1fQF7A8sT8//H3AT8O/AL4HngU93us9V0RdwGPBEenYNfk/njcBDKe1pYGKn+9yCrr6We7/kw5RO67/jyun2AEzOGZRxufR9gHuT7MwmyzoB+AFwdA3Z36aXrgHHd7rfVdFZjbJHAr8B1qYfRdcbqLL1BYwGVmR/8ADbFeR7AAd1ut8V0ted6Zkr8roCtgOuSrKlne53C/r6R+CLwKnAa4BFrRioUt+JnVZOtwegNyl8ag3Zsbn/qG1KqOvKVN5Vne53VXUGXJSenwT0DBEDVaq+cFczBlzT6b5VXV/AKLaMKMbUkO+bk+/Y6b6XpL9WDVRp+o81qBaQNBYYj0+JLCjKzex2/C/5McBflVBl5t13bAlldYR26kzSm4BPAvPN7JbWW9t5ytaXpO2BM9LH2eW1tBq04fv1Aj5zAaAa8uyczkZ8OmtYU7b+w0C1xmEp/o3VdhsPcHchbyuMS/EjJZTVKdqiM0mjgGuAx4EzB9+8ylG2vsbjU3hrzOw+SRPThpLLJZ0v6c2tNrjDlKovM/sjcFv6eL6k7TJZ+vfn08erLA0Rhjml6n/Yu9tokf1T3Og28/8t5B0UksYA09PH61spq8O0S2f/BhwMvN/MhtJt1GXr6y9SvFxSDzCtID9P0vXABxu8YKpMO75fM4Ef4SPPkyT1pvQjgd2BS4CzB9jOoUqp+g8D1Ro7pXhjgzzPpHjnwVYiaVvgO8CuwG1dPn1Vus4kTQQ+BtxoZt9voW1VpGx9vSzFx+AOOr8EzAXWp7TL8EXup4APDbSxFaD075eZrUzfsWuBk3jxFHsv8LM00gpK1n9M8bVGNifd7qH9XNwV/RrgtDbX1W5K1ZmkHYBv4y/UmWWUWTHK/o5lv/lt8Wmps83sATN70sxuBk5JdU3r0vM9pf8mk3H6NXAgcDLuA2kvXFe7A9dLOq+s+rqcUvUfBqo1nk7xTg3yZLKnG+Spi6RLgNOBdcBbzGzdYMqpEGXr7ELgIOATZtbNa3P1KFtf+TzfKgrNrBdYgr8bjmuivKpRqr4k7YafedoZeIeZ3Wxm683sMTO7CXgHvjniXyWNa1TWMKFU/YeBao3VKX5Vgzz7FfI2jaQvAx8Ffo8bp+UDLaOCrE5xWTp7D34gd5qkRfmAvzwAZqS0KwfR3k6zOsVl6SufZ1WdPFn6mCbKqxqrU1yWvt6Fj5Z+YWYri0IzWwEsxkekxzXbyCHM6hSXov9Yg2qNbNv3GyTtUGdR+chC3qaQ9EXgE/jawIlmdu/gm1kp2qGzbfDzFfU4IIXdmiyvSpStr77cv/fA//gpsmeKn6khqzpl6+uVKd7QIM+TKX5ZgzzDhVL1HyOoFjCzNfgPfnvgfUW5pGPxBdV1wF3NlitpNr4r6AncOC0tpcEVoGydmdmrzUy1Ar7tHODslHZoeT3ZOrRBX2vxv/jB1zWL5e0OHJ4+9hblVacNv8mHUzw+v8U8V952+NZ9qD8iHTaUrv9On1ru9gBMYcvJ6ANz6XvjV+685FoP/CT//cAXapT3ufTME8D4TvevG3TWoJ4ehsZNEmV/xyax5Q6+Q3Ppo4DvJVkvyV9ct4Uy9ZWe2Zie+QYwMicbCXwzyR4Hdu1030vS3yL6uUmin+/XgPVft55OK2MoBHxrruGLpbfglyFuSGk3ULgYMffi7Cmk/w1brk25O+WrFWZ1us9V0Vk/dQwJA9UOfQEXJ/km4GepjLUp7SFyd6h1YyhTX/hZsewezLXAzanMh1PaH4BTOt3nFnR1OPCLXMgudf1dPn2A368B6b9eiDWoEjCzmZJ+DnwEXwsZgf91cTXwTTPb3GRR+TnsI1Koxe10+TU1JepsWFC2vszsbEl3Av+En+jfET9A+RVgtpnVWpvqGsrUl5ldI2kZftbuaOBtSbQWvyz2K9bda8S7AG+qkT7oXYll6T9cvgdBEASVJDZJBEEQBJUkDFQQBEFQScJABUEQBJUkDFQQBEFQScJABUEQBJUkDFQQBEFQScJABUEQBJUkDFQQ1EHSakmWwhf6yTsvl3fRVmpiW5E0QtIySQ9KGlmQTZB0h6TnJD0q6TJJoxuUs0TSA8l/V608+0p6VtIP2tGXoDsJAxUEzTFV0ohaAkm74G4/hhozgDcCnzWzTVmipFcAC/HbB/4bv3NtBrCgTjkfxa/TmWF13Mib2cMkb76SjiurA0F3EwYqCPqnF9gXOLGO/P3ADvj9iUMCSTsB5+M3dF9bEJ8DjAb+wcwm4cbnNuAkSUcWytkPuACYZ2Y/6afa2fjdgF9qvQfBUCAMVBD0T0+Kp9eRTwdeAK7bCm3ZWkzD74bsMbMXCrLD8QtSvwuQ5Fcn2ZsLeS8Fnsd9mzXEzB7DLxYdL+mowTc9GCqEgQqC/lkM3AucnFyA/xlJB+Mv5R/jU101kfRWSXMkLZW0XtKmtLZzjaTX1XlmlKRZkvokPZOeeUTSXZI+L2lUIf8ESQskrZX0R0kbJK2QNF/SCQPs88wUF0dP4I4ONxQu/Hw8xX9uk6T3Aifj/rgebbLezIfXzIa5gmFBGKggaI4e/OX7gUL69BR/u5/n5wKn424b7gD+Cx9ZTAV6Jf11PrOkbYD/xP3uHIDfYH89bij3A/6FnIdgSScCP8d98TyKuzRYiPsVmwKc2mQ/kTQOeD2wwsxW18iyGthLUv72/demeFUqY2fg67grj/50k+en+Gj03fXW/IJhRKd9kUSIUNWAv4gNd3syBjcui3PyEbjLhfW4B9HMUduiGmWdAuxWSBPw4fTMveQcBALHpPQlwOgazx0F7JhLW5jyf6BG3XswAOeXwBmprGvryGck+XWp7MOANbgfob1Snq/j60mvHYTe70nlT+j0dyBCZ0OMoIKgCcxsHfAjYEJuSu5t+OaJ+Wb2fD/P32hmTxbSzMwuB+4EXoePWjL2SfEdZraxxnP/Y2bP1sj/wxp1rzezJY17+CIOTfF9deTfwp3YnQY8hrv4HgucY2a/l3QE7gfoQjO7P3tI0sgmR0WZb6XDBtDmYAgSBioImqcnxdMLcQ9NIGmspA9L+qqkqyT1SOrBR2cAB+Wy9+FTXadLmilpn2J5BX6Z4vmSjmpxemzvFK+vJTSzP+FO6M4ALge+DEw0s7mp3iuA5SSnmpLeLmkpvrFik6QfSzqwQf3ZelZ/fQ6GOOFRNwia52b8pf1BSRfjGwCWNTM6kXQ+8Cka/+Z2yf5hZg9I+ji+5XoOMEfSSny0dRNwg714d925+MjnpBQ2SlqCT/1dZ2Yrm+8mu6b4qXoZ0ojxyhTyfCy143gz25S2nd+Kj4pOBXbH19UWSnqDmT1do/is3t1qyIJhRIyggqBJ0kt5PvByfOF/JE1sAJA0GTgPeA4fdbwGXz+SmYm0XRtfW8rXdynwKnzNZx6+5nUafiC2Nx0QzvKuA8YDb8FHLn34QdrPAr+V9KEBdDWbitylYa4Ckl6Jn5262sxuT8mfxI3yZDNbYGZXAP+Mb/T4uzpFZfU+MZD6g6FHGKggGBg9KX43vmliXhPPvC/FnzKzK81spb34RoW6011mts7M5prZaWb2anx0sizFswp5N5vZQjM718yOwTcwzMINxJy8QeuHbEv4Hk3mz/gGsBE/yJvxl8BjZrYil3ZXTlaLrN5mt6YHQ5QwUEEwAMysD9/OvR5YYM2d78m2Y68pCtKGi6Y3A5jZUuCS9LHeCz7Lu9HMLgIewrfIH9xkNX0pfn3DXDkkTQEmAR83s8dzoo34LRt5sjv7rE5xWb19deTBMCEMVBAMEDM72sz2NLN6U1RFsp1sZ0jaPkuUtDd+MPUl61KSTpD0TknbFtJHAO9MHx/MpZ+VrhUqlnMEPiW5GTdUzfDTFBdvhahJGpldAvzEzOYXxL8CRkvKn8OanuJ7apS1E37/31OEgRr2xCaJIGg/X8MP5L4LWCFpMT6qOBYfVd2In5PKcwjwVWCDpD78lood8XWllwPrgIty+T8NXCzpPnx7+CZ8nWci/ofobDOre9NFHjNbJelXwCGS9jezVf08ciG++WFGDdls4O+BeZKmpnwTgZVsWXvLczy+1narvfSKpWCYESOoIGgzaQfd4cD38I0Qk/BzT1fgo5QNNR67Bd9w0IevUU0GjsYN02eAQ8zswVz+j+Cjsc34S/49wCtSOW83s3MH2OzLUjy1USZJE3DDdEGtnYJm9jv8vNhifAPHG/FNHscVz3clphXqD4YxMqs3DRwEwXAl+XZ6EJ9qG7c1RjOS9sSnIX9tZke0u76g+sQIKgiCl5BGN58B9qefUVSJzMK37p+1leoLKk6MoIIgqEnakHEPfnD3IMs5LWxDXfsCK4AfmtnkdtUTdBdhoIIgCIJKElN8QRAEQSUJAxUEQRBUkjBQQRAEQSUJAxUEQRBUkjBQQRAEQSUJAxUEQRBUkv8H+WCNNp3r010AAAAASUVORK5CYII=\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "plt.plot(num_sol_imp[:,2]/m0, num_sol_imp[:,1], label = \"Implicit\", color = \"purple\")\n", | |
| "plt.plot(num_sol_exp[:,2]/m0, num_sol_exp[:,1], \"--\", label = \"Explicit\", color = \"orange\")\n", | |
| "plt.plot(m/m0, -np.log(m/m0)*250, 'o',label = \"Tsiolkovsky\", color = \"red\") \n", | |
| "plt.xlabel(\"Mass (%)\")\n", | |
| "plt.ylabel(\"Velocity (m/s)\")\n", | |
| "plt.legend(loc=\"best\")" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "As the mass decreases the velocity increases. This makes sense for two reasons. For one, less mass means it can accelerate faster. Secondly, as the fuel is used the rocket accelerates" | |
| ] | |
| }, | |
| { | |
| "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": 7, | |
| "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 = np.array([state[1], u/state[2]*dmdt-9.81-c/state[2]*state[1]**2, -dmdt])\n", | |
| " return dstate" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 8, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "y0 = 0\n", | |
| "v0 = 0\n", | |
| "m0 = .25\n", | |
| "\n", | |
| "dm = .05\n", | |
| "dmdt = .05\n", | |
| "N = 25\n", | |
| "\n", | |
| "t = np.linspace(0, 4, N)\n", | |
| "dt = t[1] -t[0]\n", | |
| "u = 250\n", | |
| "mf = .05\n", | |
| "m = np.linspace(mf, m0, N)\n", | |
| "\n", | |
| "num_sol_imp2 = np.zeros([N,3])\n", | |
| "num_sol_exp2 = np.zeros([N,3])\n", | |
| "\n", | |
| "num_sol_imp2[0,0] = y0\n", | |
| "num_sol_imp2[0,1] = v0\n", | |
| "num_sol_imp2[0,2] = m0\n", | |
| "\n", | |
| "num_sol_exp2[0,0] = y0\n", | |
| "num_sol_exp2[0,1] = v0\n", | |
| "num_sol_exp2[0,2] = m0\n", | |
| "\n", | |
| "for i in range(N-1):\n", | |
| " num_sol_imp2[i+1] = heun_step(num_sol_imp2[i], rocket, dt)\n", | |
| " num_sol_exp2[i+1] = rk2_step(num_sol_exp2[i], rocket, dt)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 9, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "<matplotlib.legend.Legend at 0x7f749bc15ef0>" | |
| ] | |
| }, | |
| "execution_count": 9, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEaCAYAAABEsMO+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3RVxfbA8e+kQgqEkEYISYjSBEEEpCpNBZUqKAj4wPJAxF55+nxP5SdiRwVBUSkC+rBRRAVEQFRUmiKCgEIglCQQQiChJCT798dJQspNuSV9f9a663rPmZkz3GWyc+bM7DEiglJKKVXZuFV0B5RSSilbNEAppZSqlDRAKaWUqpQ0QCmllKqUNEAppZSqlDwqugNVXVBQkERHR1d0N5RSqkrZvHnzMREJLq6MBignRUdHs2nTporuhlJKVSnGmP0lldEhPqWUUpWSBiillFKVUqUNUMaYycYYyX49Uky5EcaY9caYFGNMqjFmkzFmgjGm2H+bo/WUUkqVj0r5y9gY0wF4DCg2D5MxZjqwAGgPrAdWAU2BacAnxhh3V9ZTSilVfipdgDLGeANzgARgSTHlhgB3A/FAaxHpJyKDgSbATmAwcI+r6imllCpflS5AAc8ClwB3ASnFlPtX9vvjIrIn56CIJADjsz9OtDFk52g911mwAKKjwc3Nel+woMwupZRSVVWlClDGmI7Aw8BCEVlWTLkIoB2QDnxc8LyIrAMOAWFAJ2frudSCBTB2LOzfDyLW+9ixGqSUUqqAShOgjDG1gLnAceD+Eoq3zX7/Q0TOFFFmY4GyztRznSefhNOn8x87fdo6rpRSKldlWqj7HNAMGC4ix0oo2zj7vbiFXgcKlHWmnuscOGDfcaWUqqEqRYAyxnQBHgAWi8j/SlHFL/s9rZgyqdnv/i6ol48xZiwwFiAyMrKYpmyIjLSG9WwdV5WGiHDq1ClOnjzJ6dOnyczMrOguKVUpeXh4ULduXQIDA/HwcG1IqfAAZYypDcwGTmLNritVtex3e7cDdrRePiLyDvAOQPv27e1r67nnrGdOeYf5fHys46pSEBESExNJS0sjMDCQsLAw3N3dMcaUXFmpGkRESE9PJykpibi4OKKionBzc92To8rwDGoy1hqkh0TkSCnrnMp+9yumTM65U3mOOVrPdUaOhHfegagoMMZ6f+cd67iqFE6dOkVaWhpRUVEEBATg4eGhwUkpG4wxeHt706BBAzw8PEhOTnZp+xV+B4W17igLGG2MGV3gXPPs9/HGmH7AXyJyJxCbfTyqmHYbZb/H5jnmaD3XGjlSA1IldvLkSQIDA3F31/XaSpWGMYaAgACSk5OpX7++y9qtDAEKrDu57sWcj8l+BWR/3pr93tIYU7uIGXkdCpR1pp6qQU6fPk1YWFhFd0OpKsXHx4fDhw+7tM0KH+ITkWgRMbZeWNPOAR7NPnZZdp04YAvgBdxUsE1jTHcgAitbxIY813KonqpZMjMz9e5JKTu5ubmRlZXl2jZd2lr5ej77/QVjzMU5B40xIcBb2R+niEjBb8zReqoG0WdOStmnLH5mKssQn91E5BNjzAys9ES/G2O+ATKA3kAdYDFW8leX1FNKKVW+qmyAAhCRu40x3wMTsJ5huQN/Au8DM4q6C3K0nlJKqfJTqQOUiIwBxpRQZiGw0IG2HaqnlFKqfFTlZ1BKqQoSHR2NMYa1a9dWdFdsio2NxRhDdHR0oXM5fY+NjXXqGk8//TTGGJ5++mmn2lFF0wCllFIutHbtWowx9OjRo6K7UuVV6iE+pZRytdWrV5ORkUHDhg2daueee+5h+PDhBAUFuahnqiANUEqpGuWiiy5ySTtBQUEanMqYDvEppVxizJgxGGOYM2cOf/zxB0OGDCE4OBg/Pz+6devGmjVrcst+8cUXdO/enbp161KnTh0GDBjAnj17CrWZd7gsLS2NiRMnEhMTg7e3N40aNeLee+8lKSnJrn4W9wxKRFi0aBHXXXcdISEheHl50bBhQ3r37s20aflXn9h6BtWjRw969uwJwLp16zDG5L50yM9+GqCqCt0mXlURmzZt4oorrmD37t307t2bZs2a8cMPP9CnTx/Wr1/Pm2++ycCBAxER+vTpQ2BgIMuWLeOqq64qMtikp6fnBolWrVrRv39/zp49y7Rp0+jcuTMJCQlO9zs9PZ1BgwYxbNgwVq1aRdOmTRk6dCjNmzdn+/bt3HvvvSW20bdvX/r06QNAaGgoo0ePzn317dvX6T7WNDrEVxXkbBOfs0VHzjbxoElnK8Az5pmK7oLD/iv/LfNrTJ8+nVdeeYWHHnoo99jjjz/Oiy++yJ133kl8fDxr167lyiuvBODs2bNce+21rF+/nrfeeounnnqqUJsbNmygadOm7Nq1K/fZ0alTpxg8eDCrV6/m3nvvZdGiRU71+9FHH2Xp0qU0bdqUJUuW0Lx589xzmZmZLF++vMQ2Jk6cSKdOnVixYgXNmzdnzpw5TvWpptM7qKpAt4lXVUjnzp3zBSewfnED7N69mwkTJuQGJ4BatWrx4IMPAuQbBizolVdeyTexwd/fn5kzZ+Lu7s6nn35KXFycw31OTExkxowZuLm58dlnn+ULTgDu7u4MGDDA4faVYzRAVQW6TbyqQmwNZdWrVy93GwZb55s0aQJQZDbsgIAA+vXrV+j4xRdfTKdOncjKyuK7775zuM/ffvstGRkZdO7cmZYtWzrcjnItHeKrCnSb+EqlPIbJqrKIiAibx/38/EhKSrJ53s/P2if07NmzNuvaWnCb99wPP/zAwYMH7e9stv3ZP18F75xUxdI7qKrgueesbeHz0m3iVSVV0pbfrtwSPC/NQF/9aICqCnSbeFXDFZeWKOdceHi4w+1HRVmbbO/atcvhNpTraYCqKkaOhNhYyMqy3jU4qRrkxIkTfPnll4WO7927l59++gljDFdddZXD7ffq1QtPT09+/PFHdu7c6UxX8fLyAuD8+fNOtaM0QCmlqoiHH36YI0eO5H5OTU1l/PjxZGZmMnjwYCKdeCYbEhLCXXfdRVZWFkOGDGH37t35zmdmZrJs2bJStZUz0/Cvv/7SIOUknSShlKr0OnfuTGZmJk2bNqVXr154eXmxbt06jh49ykUXXcT06dOdvsZLL73E33//zZdffknLli3p3LkzERERJCYm8vvvv5OYmIiIlNhOVFQUbdu2ZevWrbRu3Zp27drh7e1Ns2bNePTRR53uZ01i1x2UMSbAGDPIGPOMMWamMeYjY8yM7M8DjTEBZdVRpVTN5eXlxbfffsu4cePYtm0bS5cuxcvLiwkTJvDTTz8RFhbm9DW8vb1ZtmwZH3zwAVdddRXbt2/nk08+4c8//6R169Z2BcHPPvuMm2++mePHj/Phhx/y3nvvlWqhr8rPlPQXgTHGHbgRuBu4EsiZKpN3yozkef8OeAv4XEQyXdrbSqh9+/ayadOmiu6GcqGdO3fSokWLiu6GwsrF17NnT7p3715p955SF9jzs2OM2Swi7YsrU+wQnzHmFmAKEIEVkI4BPwE7gOPASaAOUB+4BOgE9MDaRj3OGDNRRD4qVW+VUkqpPIoMUMaYH7ACzlHgdWCuiPxWUoPGmMuwtmm/BVhgjLlXRLq6prtKKaVqiuKeQV0EPAREishDpQlOACLyq4g8ADQCHs5uRymllLJLcUN8F4lImqMNi0g6MNUY846jbSilarYePXqUauacqp6KvINyJjgVaOd0yaWUUkqp/HShbnWmmxwqpaqwUgcoY0ygMeZyY0xggeMNjDFzjDFbjTGfG2PauL6bym45mxzu3w8iFzY51CCllKoi7LmD+hewEWvyAwDGGC/ge+BWoA0wEFhjjGloswVVfnSTQ6VUFWdPgOoJ7Cswm28Y0BhYB/QFpgMBwD0u66FyjG5yqJSq4uwJUBHAXwWO9cPKHnGniKwUkXuBfcB1LuqfclRRiTN1k0OlVBVhT4Cqh5VJIq/OwG4R2Zvn2FbyDAOqCqKbHCqlqjh7AtQZrJRGABhjGmHdVf1QoNw5wNv5rimn6CaHSqkqzp7tNv4EuhljAkXkODCCC8lh84oAElzUP+WMkSM1ICmlqix77qA+AHyBX4wxi4BngVRgSU4BY4w3cDmg+yYrpZRyij0BagawEIgBhgLpwD9FJCVPmf5YQWydy3qolKp0oqOjMcaU+KqoLTJyrl9Qjx49XNKvOXPmYIxhzJgxTrWjilfqIT4RyQJGGWOeAkKBHSJyskCxvcBNFH4upZSqhvr06VPsZoGu2EiwKomNjaVx48ZERUURGxtb0d2p8orbbmM4sFxETuU9LiL7sKaSFyIiW4AtLu2hUqrSmjhxIj169KjobpTavHnzOH36NJFOLrcYPHgwnTp1om7dui7qmbKluDuohcA5Y8xq4HNgqYgcLZ9uKaWU6zkbmHLUrVtXg1M5KO4Z1OPAr1iLbt8BDhtj1hhj7jPG6GpPpVSpiAjXXXcdxhjGjh1b6HxWVha9e/fGGMM991xIQhMbG4sxhujoaM6fP8+UKVNo0aIFtWrVIjQ0lNGjR3PAzswoJT2DWrFiBTfeeCPh4eF4eXkRFhZG165deeGFFzhz5kxuOVvPoMaMGUPjxo0B2L9/f75ncdHR0Xb1U1mK227jJRHpjDVt/D6s6eRdganAPmPMJmPMv4wxzcunq6rMafZzVQaMMXzwwQc0bNiQWbNm8eGHH+Y7/+yzz/Ltt9/Stm1bXnnlFZttDBs2jP/+979ERkYyaNAgvLy8mDdvHh06dGDXLucnDYsI48ePp2/fvnz++ec0bNiQIUOG0KZNG+Li4pg4cSIJCcWvnunWrRtDhgwBwNfXl9GjR+e+hg4d6nQfa6ISJ0mIyBGsHHvTjTH1sBLC3ghcjTWl/P+MMbuBT4HFIrKpDPurykpO9vOcBLM52c9B11LZsu1p2P5M6cpe9E/oWGDfzp/Hwt+zSle/1X+h9dP5j63tD4e/KF39K96GiwvfuZSnoKAgPvroI3r06MG4ceNo3749TZo0Yc2aNUyaNAl/f38WLVqEt3fhNf779+/nzJkzbN26lUsuuQSA9PR07rjjDubPn8+tt97KL7/84lT/pk6dysyZMwkNDWXx4sV06tQp95yIsHbtWurVq1dsG3feeSdXX301n376KUFBQcyZM8epPik794MSkWQRmSMiA4BgYDjwMRAOPAH8bIzZb4x5zRjT3dia56kqJ81+rhzQs2fPIqeYBwQE5CvbrVs3nn32WU6dOsXNN9/MgQMHGDFiBFlZWcyaNYuLL764yOs89dRTucEJwMvLi2nTplG3bl02btzIDz84PnH4/PnzTJ48GbCG7vIGJ7DuAHv27KnPnCqAPZkk8snecXcRsCh7241rsO6s+gP3Yw0L/hf4Pxf0U5U1zX6uHFDcNHOfgrkggX/961989913rFixgtatW5OSksK4ceMYNmxYsdcZNWpUoWN169alX79+LFiwgLVr19K1a1eH/g2bNm3i2LFjRERE0LdvX4faUGXD4QCVl4ikA8uB5cYYN6A7MBhIdEX7qhxERlrDeraOq8JaP1142M0eHd8pPOxnjx7LHK/rQvZOM895HhUTE0NKSgqXXHIJU6dOLbZOQEBAobuxHDmTDw4ePFjqPhS0P/v/+2bNmjnchiobLt/yXUSyRGSNiNwnIk78BKpypdnPVTlZvHgxqampgBVYDh065HSb+jShenIoQBljwrK3f+9S1MvVHVVlTLOfq3Kwfft27r//fry8vBg5ciQnT55k2LBhpKenF1nnxIkTpKSk2DyXk60hPDzc4T5FRUUBuGQ2oHItuwKUMeYmY8xO4BDW9u/ri3gVzHCuqoKRIyE2FrKyrHcNTsqF0tLSuPnmmzlz5gwvvPAC8+bNo2fPnmzevJlHH3202LoLbCx5SElJ4YsvrJmMzmSzaNeuHUFBQRw8eJAVK1Y43A5YkzfAmnihnFfqAJWd+ugjoBlwEtgG/FjEa4PLe6qUqtImTJjAzp07GTBgAA888ABubm4sWLCAkJAQ3njjDRYvXlxk3WeffZadO3fmfs7IyOD+++8nJSWFdu3a0a1bN4f75enpyb/+9S8AbrvttkJT1nOmmRd1F5dXcHAwXl5eJCQkkJyc7HCflMWeSRJPZL/fD8wQEf0TQakabsqUKcWu9xkxYgTXXnst8+bNY+7cuTRq1IjZs2fnnm/QoAEffPABffv25fbbb6dt27a5Q245IiMjadeuHZdddhm9evWibt26bNiwgQMHDhAUFMS8efOc/nc8+OCD7Ny5k3fffZdOnTrRvn17Lr74Yo4fP86OHTuIi4tj3759JU419/T05IYbbuDzzz+nbdu2dO3aldq1axMUFMSUKVOc7meNIyKlemHtqLu+tOVryqtdu3aiqpcdO3ZUdBcqvaioKMHasLTY12uvvSY7d+4UX19f8fDwkO+//95mexMnThRAOnXqJOnp6SIism/fPgEkKipKMjIyZNKkSdK0aVPx9vaW4OBgGTVqlOzbt89meznXL6h79+4CyJo1a2zWW7Zsmdxwww0SHBwsnp6eEhoaKt26dZMXX3xRzpw5k1tu9uzZAsjo0aMLtXHs2DG54447JCIiQjw8PHL/DTWBPT87wCYp4ferscqVzBhzCFgnIiNcFx6rvvbt28umTZo8ozrZuXMnLVq0qOhu1Hi6dUXVY8/PjjFms4i0L66MPZMkVgId7CivlFJKOcyeAPVfoI4x5gVjjHtZdUhVMZpgVilVRuzZUfeAMaYbsAS4MXufqINAVhHlJ7umi6rS0gSzSqkyVOoAlZ349R6gCeAOXIT1ILJQ0ezjGqCqu+ISzGqAUk6Kjo6mtM/IVfVkzzTzicC9wHmsvHt/Aall0SlVRWiCWaVUGbInQN0BnAa6icivZdQfVZVoglmlVBmyZ5JEQ+A7DU4qlyaYVUqVIXsC1CGsOyilLJpgVilVhuwZ4vsfMNYY4yvWZoVKWcFIA5JSqgzYcwc1CdgNLDXGXFRG/VFKKaUA++6glmLN4OsJ7DTG7KXodVAiIn1c0D+llFI1lD0B6uoC9Zpmv2zRxQtKKaWcYk+AuqbMeqGUUkoVYE+qo9Vl2RFVgyxYYGWbOHDAWjP13HM60UIpVYhdW74r5bSc/H3794PIhfx9mmS2SjDG2P0aM2aMQ9dq3749xhhcsZ1Namoqxhj8/PwKnQsKCsIYw7Fjx5y+TkXavn07xhhatWpV0V1xGXuG+JRynubvq9JGjx5d6Fh8fDwrVqzA19eXoUOHFjrvzHbsqmYrMkAZY74DJorIj442bozpCjwvIlc52oaqZjR/X5Vma3v3tWvXsmLFCoKCgord/t1en376KWfOnCE6OtplbaqqpbghvubAemPMKmPMMGOMd2kaNMZ4G2NuMcZ8A3xH0TP9VE1UVJ4+zd+nCoiKiqJ58+bUqlWroruiKkhxAaoJ8CbQHVgIJBhjlhtj/m2MGWKM6WGMuTz7fYgx5iljzJdAIjAfuBJ4HQ1QKi/N32dbDdn4ce7cuXTv3p169erh6elJcHAwbdq04b777uNAgbvo4p5BnTt3jldffZX27dvj7++Pj48PrVq14qmnniIlJcUlfc3MzOTuu+/GGEPr1q05ePBgvvO//vorI0aMoGHDhnh5eRESEkL//v359ttvC7V16aWXYoxh9eqi55rdddddGGN45plnco+lpaUxadIk2rRpg6+vL97e3oSHh9O1a1f++9//cv78+VL9W5KSkujSpQvGGG655RbOnTvHY489hjGGRx55pMh6CxcuxBhDr169SnUdlxORYl9ADPAacBxrUW5mMa8s4BjwIhBdUtvV4dWuXTtRdpo/XyQqSsQY633+/IruUT47duwo3wvOny/i4yNiTRuxXj4+le57KcqaNWsEkKioqGLLPfzwwwKIl5eX9OzZU2655Rbp27evNG3aVABZtmxZvvLt2rUTQDZu3Jjv+KlTp6RTp04CiL+/vwwYMECGDh0qwcHBAsjFF18scXFxheoA4uvrW6hf9evXF0COHj2aeywtLU369+8vgPTu3VtSUlLy1fnwww/F09NTAGnTpo3ccsst0qVLFzHGCCBTpkzJV/6ll14SQEaNGmXzuzl79qwEBASIMUb27t0rIiIZGRm5/87AwEC54YYb5JZbbpGePXtKgwYNBJBTp07ltvH7778LIC1btszX9l9//SVNmjQRQB599FHJysoSEZHY2Fhxd3eXwMBAOXPmjM1+de3aVQD55JNPbJ4vyJ6fHWCTlBR/SiqQWxBqA32BKcDXwBbgb2Az8BXwHNZiXu/StlkdXhqgqp9yD1BRUfmDU86rhF/4lUVpAtSJEyfEw8NDAgMDZd++fYXO79ixQw4cOJDvWFEBavz48bmBIT4+Pvf4qVOn5LrrrhNArr766nx17AlQCQkJ0qFDh9yAkp6enq/8vn37pHbt2gLIzJkz85378ssvxcvLS4wx8t133+Uej4+PFw8PD/Hx8ZGTJ08W6sP//vc/AaRHjx65x5YvXy6AdO3atVAAyczMlLVr1+brm60A9fPPP0twcLC4ubnJtGnTCl130KBBAsjs2bMLnfvtt98EkIYNG0pGRkah87ZUWIDSlwaomqLcA5QxtgOUMeXbDweVJkDt3bs395dtadkKUMnJyeLt7S2A/Pjjj4XqxMXF5Z7/9ddfc4+XNkDt2rVLYmJiBJAnnnjCZr8ee+wxAeTaa6+1ef7uu+8WQAYNGpTveL9+/QSQ9957r1CdnMA6Z86c3GPvv/++APLkk0/avE5BBQPUkiVLxMfHR2rXri2LFy+2WWf16tUCSIcOHQqdGzt2rADyzDPPlOr6IhqgKt1LA1T1o3dQ9ilNgMrMzJTQ0FAxxsgTTzwhe/bsKbFdWwHqq6++EkCaNm1aZL2cQDB16tTcY6UJUIsXL5b69euLu7u7vP3220W237FjRwFk4cKFNs9v2rRJAAkICMh3/JNPPhFArrzyynzHjxw5Iu7u7uLn5yepqam5x7dt2ybGGAkICJBZs2blG4K0JW+Amj59uri5uUlwcLD89NNPxdZr2bJloe85JSVFfH19xdPTUw4fPlxs/bxcHaB0oa6qGqrzJIIaMHHEzc2NDz74gHr16jF58mSaNGlCWFgYgwcP5u233yY1NbVU7Rw6dAiAxo0bF1nmoosuyle2tIYOHUpSUhJvvPEGY8eOdbgPOdc/ceIEp/Os+evfvz/169fn+++/Z+/evbnH58+fT2ZmJkOHDsXX1zf3+KWXXsrzzz9Pamoq//znPwkODqZJkyaMGTOGJUuWkJVlK0837Nq1iwkTJuROyujYsWOx/+57770XgLfeeiv32Ny5c0lLS+PGG2+kQYMGxdYvSxqgVOVX3bNP1JCNH6+55hr279/PwoULGTduHMHBwSxZsoS77rqLJk2asHPnzhLbsP7wtjJalFTGXv/4xz8AmDx5Mrt27XKqD7Z4eXkxYsQIRIS5c+fmHs/5b1sZNx5//HFiY2OZNm0aw4cP5+zZs8ydO5dBgwbRtWtXzpw5U6hOZGQk3bt3JzMzk/vuu6/E4D9q1CgCAgL46KOPSE5OBmDGjBkA3H333Xb9G12upFssfZXNEN+GqRtk41vfSdqhfQ7Vr1HKeQis3If4qrjSzuKzJS4uTgYOHCiAXHPNNfnOOTrElzP7zt4hvqNHj8p//vMfASQkJES2bdtms/2Shvg2b95sc4gv77no6GjJysrKHQ5s3Lhx7uy6kmzcuDF35uNzzz2XezzvEN/p06elT58+Akjnzp3lxIkTxbb50EMPCSCvvPJK7nOpVq1alao/eekQXzWQmZ7Jd5O+Y//8qdRafRGHXr6cg/97ifOpJyu6a5WTZp+otiIiInLX/fz2228llu/UqRPe3t7s3r2bn3/+udD5w4cPs2rVKgB69Ohhd3+eeeYZpkyZQmJiIj179mTz5s2FynTv3h2AefPm2Wxj9uzZRV7/8ssvp3Xr1sTGxrJu3brcu6fRo0eX+o6sffv2uXc2RX1ntWvXZunSpQwaNIgNGzbQq1cvkpKSimxzwoQJuLm5MXPmTKZPn557rKJpgKoAf634izNJZ2jddRtu7lk0DN9KROZjZC4KJu6VXsR//QGSmVnR3aw8NPtElbd7927mzp1rc7hp2bJlgJU5oiQBAQHcdtttgDX8dPTo0dxzaWlpjBs3jrNnz3L11VfTpk0bh/r6+OOP8+abb3L8+HF69+7Nhg0b8p0fP348tWvX5uuvv+a9997Ld27lypXMmjULYwwPPfSQzfZzhvJmzZrFhx9+iDEmd3gxr6+//ppVq1aRWeB3QUZGBl9//TVQ/Hfm5eXFxx9/zPDhw9myZQs9e/YkISHBZtmYmBiuv/569uzZw2effUadOnUYNWpUkW2Xm5JusfTl+iG+tKNp8vObP8jhF5qKLKDQK+VNf3mt0cuy6vFVkvB7gt3tVzvlvJBVh/jsU5ohvvXr1wsg3t7e0qlTJxk+fLgMHTpUmjdvnnv8m2++yVentAt1Bw4cKEOHDpWQkBCXLdQVEXn33XfFzc1N/Pz8ZM2aNfnOLVy4UDw8PASQyy67TEaMGCFdu3bNXaj7/PPPF/ldJCYm5i7ypcDap7wmTZokgNSrV0969+4tI0aMkIEDB0poaKgA0qhRo3wz7IpaqJuZmSm33XabANKsWTM5ePCgzeutXLkyt0/33HNPkf0vjk4zr2QvZ6eZH9+2Sfa+frskTQ3JDVDf39JFnubp3NeMNjNk62vz5ORfvzt1rSqtHLNPaICyT2kC1PHjx+Xll1+W/v37S0xMjPj6+oq/v7+0aNFCJkyYILt27SpUp6gAJSJy5swZefnll+Xyyy8XX19fqVWrlrRo0UKefPJJSU5OLlTekQAlciEQ1a5dW77++ut857Zs2SLDhw+XsLAw8fDwkPr160u/fv0KBVpbcp67UWDtU147d+6Uf//739K9e3eJiIgQb29vCQoKkssvv1wmT54sSUlJ+coXFaBERLKysmTChAkCSExMjM3F0mfPns1dQ+boz4CrA5SxypXMGHMnsEBECk8bqcHat28vrtivRrKySFz/NWkb32HNnBgO/lE33/mRj3/Axa3/JmsinAIAACAASURBVP5IMzLCRhI24B48/es5fV1V2M6dO2nRokVFd0OpcrVgwQJGjRpFr169is0ZWBx7fnaMMZtFpH1xZex5BvUOcNAY84oxpokd9VQpGDc3QrtfT8wjixmz5SWGLx1Oy5tb4u7tjl/AKWJaWesmwhrsopH5D1kfh3Hg1Ws5uu5zpIj1EEopVRrnzp1j8uTJAEU+O6sI9mxY+AVwHfAgcH/2dhrTgS+ktLdhqlTcvdxp1r8Zzfo342zKWfZ+vpwjh1vRIHw7bm7WV+1dK53IsFVwaBUnpoeQ4jmYoOsewjdKk8crpUpn5syZ/PTTT2zYsIHdu3fTs2dPbrjhhoruVq5S30GJyACszOZTsDKWXwssBvYZYyYaY4Id7YQxxtMY0zv77uwnY8wRY0y6MeaQMeYTY0yPEuqPMMasN8akGGNSjTGbjDETjDHF/vscrVeeatWtxSVjhtDwsW2kdt7J3uN3cTwxNF+ZgPqJRNV5m9rrmxP3Ykd2LfmTzIwaPguwOmeeUMpFvvnmG+bOnUtSUhIjR45k0aJFFd2lfEr9DCpfJWM8gZuBu4HOWA/70oFPgLdEZEMx1W21dzWwKvtjPFaG9DTgEqBV9vFJIvIfG3WnZ/fjLLAayAB6A/7A58BNIlLot7Wj9Qpy1TMoe0hWFvGrl3Bu60wa1F+Hd+1zuee2/XApn781BN9QX1rf2pq2t7Ul+BKH/3aomnIyT+TdWt7Hp9TZGfQZlFKOcfUzKIcCVIGLtAEmALcAOQnFfsUa/lsgIueKqpunjV5YweJ1EVlf4NwwYAHgDvQSkTV5zg3BCorxwFUisif7eCiwBmgBPCAirxdo06F6tlREgMrrXHIyR5ZMxytxAeERfzJv8j/Y90dMvjJ9x/9Bg45NCBv4IF4BQRXU03IUHW2lQyooKgpiY0usrgFKKcdU5CQJm0TkN+AZYDZgsl9tgVlArDHmjlK08a2IDC0YnLLP/Q+Yk/2x4Mqxf2W/P54TZLLrJADjsz9OtDFk52i9Sse7Xj2ix/yb8Md2knzpZsL7j8Q/3D/3vKd3Om0uX0qk52T4rCEHXutD0oavrNVE1ZVmnlCqWnDqF7Ax5mpjzGfAPqy7qLPA+1h3U18CIcA7xpj7nOzn1uz3iDzXjgDaYQ0tflywgoisAw4BYUAnZ+tVBfUuvZyrn7+WB/Y/wIjlI7hk6CW06rKDWj7WTaxXrXQiQ1dSf9/1JL0RyYE5E0lPTqzgXpcBzTyhVLVgd4AyxtQ1xjxgjPkTWAEMAg4DTwARInKniPxPRPoDXbCeJTkboHKmtR/Jc6xt9vsfxazN2ligrDP1qgw3DzeaXN+Emz6+id7z32RvygMkJeRPmV8/+CCRXi9glkQQ9+rVHP1+afW5q3LB9hU6MVUp+5TFz0ypA5Qx5nJjzLtYdxevAE2BdcAQIEZEXhCR43nriMjPwHLA4T9djTFhwJjsj5/mOZWzGYuNhw25csZ08m7c4mi9Ksk3PIKY8a8ReP9BjgR9yv64nmSc88w97+mVQaOw1QQfGMjuf3dh09ubOHeyxMeGlZuT21d4eHiQnp5exp1UqnrJyMjA3d3dpW3asw4qZybAaeBd4E0R2V6Keml2XieXMcYDmA/UBVaLyLI8p/3ytF+UnMyU/nmOOVqvSjNubjS49ka49kZOxx/m4JLX8D/9IUGhFzZ12/F9KL99t5yVD62k1YhWdBjfgQaXV9xmZU4ZOdLh/ZTq1q1LUlISDRo0sHvPH6VqqpMnT+Lv79pfmfYM8cUCj2IN440rZXAC+CfgWWIp22ZiTf2Oo/AEiZzfHPbeVzpa70IDxozNXjO1KW825arCJyycxuNeov79B4gPXULswas5lezPHz+1BCDjdAZb393KO+3eZu8z7Ymb9xQZp5IruNflJzAwkHPnznHw4EFOnTpFZmamDvkpZYOIkJ6ezrFjx0hOTiYwMNCl7dtzZ3ORIxkjsuvYvWrUGPM6cAfWVPDeIhJfoMip7Hc/ipZz7lSeY47WyyUi72ClfqJ9+/ZV9jeXcXMjrPcA6D2AM0mnuDprB1ve2ULidmviRPQl+4hpshnYzLlFL3Hg1LXU6fkYAW26VWzHy8qCBfDkk3gcOEBUTAzJr71Gcps2HD58uMjttZWq6dzd3fH39ycyMhJvb2+Xtm1PgFphjPlaRF4trpAx5kHgOhG51tFOGWNewZpYcRQrOO2xUSw2+724TWQaFSjrTL1qrXZ9fzre25Er7rmCuB/j2DRjE00DLzzy8659jsjay+CPZSSsbkpGxB00GHAP7rV8imm1CimwuNft77+pP3w49avh1utKVRX2ZDPPAuaIyO0llJsF3C4iDj0tM8a8iDWUmIQVnGxuGWmMaYQ1mSEdCLA1I88YE4c1Nb2biPzgTL2iVPRC3bKUdvggCYtfIiDjIwKDC09HP53qy9Fz/QnsOxH/Jo5tDldpOLm4Vylln3JZqGuDF+DQeIgxZgpWcEoGrikqOAGISBywJft6N9loqztWkIkHNjhbrybyDY8g5u7XqXfPEQ7V+YADcR3JPH/hfxkfvzSi6n+E389t+fmB8fy96m8kq4qOeOriXqUqHZcGKGNNeWqHlUzW3rqTgMeBE1jBaWsJVQCez35/wRhzcZ62QoC3sj9OEZGCAdPRejWScXejYb9RRD7+E2ldd7D36G2cTA7IPS/Ajx/UZv6185neYjo/vf4TZ0+crbgOO0IX9ypV6RQ7xGeMWZnn49VYC3J3FFHcA2tBbTjwiYgMK3UnjBkALMn+uAn4o4iif4rIlAJ138JKT3QW+IYLSV/rYGVbH1pEsliH6hVUnYf4ipN57hyHlszCfd87nIzPYNHU4fnOBzVKZcADW/Dp8ij1O/WtoF7awckEs0op+zidLDb7uVMO4cIU7eJsAwaKSHELYQteZwxWLr+SrBORHjbqj8BKtXQpVlLZP7FSLs0o7i7I0Xp51dQAldexP+PZOGMrv835LXeRb+9hq+g2wHp8dzShMWdDbyd88IO41/atyK4WL3sWHwcOWHdOzz2nwUmpMuKKANU75z+BlVipjV4uong6cEhE9jrQ1ypLA9QF6anpbFuwjc0zf+TWuybi459//snpVB+Onu1P4HVP4N+kdQX1UilVGbh0uw1jzHpgecEhtppOA1RhIkL86qWkb36NhqE/4OF1Pv/5LMORhLa4XXIvoX1uxbi5Nj2KUqryc+ksPhG5UoOTKg1jDA2uHkjU42s522sPe5PGknK83oXzbkJ4gy2EJd9GyvRwfp35FempVTT3ne7cq1SZqfT7HamqzS8ymph738ZvbDwHPKZxKK5VvvPpp4Ul43/m1Yav8tX9X5G0O6mCeuqAnIkV+/dbmeD377c+a5BSyiWKHOIzxjyR/Z8zRCQ5z+dSEZHJznauKtAhPvsd3/ozJ1Y/T3jACr75qDebV3fId77bbWm06B9JWL+xuHl6VVAvS0EX9yrlMKeeQWXP4BOghYjszvO5xOtipeCrEQ8WNEA57tzxJH5bsJ1fpv2e585JGPvc2zSIjufUiQCOuw0jZNAT1A6rhOuR3Nxs76FlDGjuPqWKVZoAVVwuvslYAelYgc9KuYR3YH2uuLc7HSZcxd7Ve9k4bSOn/1xFg2grL7B/wAn8eZvzK94l7thV1O76GEGVaU1VZKTtOyhd3KuUS5R6Fp+yTe+gXOvErp0cXz6JUN+l+PoX3rLraEJjzjUYS4NB91V8olpd3KuUwyoqF59SDgto1oKYhxbiNSKe2PRnSTgcne98cOg+IrL+xdm5ofz95nhSE1JtN1QenNy5VylVPL2DcpLeQZUtESFx7Rec/eVlGob9gIfnhexTG77sxDeLrqflzS254t4riOgYUYE9VUrZw6V3UMaY8caYdGPMDcWU6Zdd5k57OqpUUYwxhPbsT9Tj6zjbYzd7j93BqRN1kSzY+E0HsjKy+H3B77zX6T1mXTGLvz+Yxfk0m/tMVg66bkqpUrNniO9G4DjwVTFlvsouM9SZTilli19UDDH3vYvP7fHsdXsf35j8e1Cd3PMn0XIXGQvD2P/GME799XsF9bQIum5KKbvYk+ooDiub+DUllFsFNBORGjGVSYf4KtbhzYfZOG0jv3/4O1f2+4buN67LPSdZhsPx7XC/9H5CrxmBcavgR666bkqpXK6eJBEMJJSiXCIQYke7SjksvF04A2cP5MG4Bwnv3pGUpMDcc8ZNaBi+ibCkWzkxLZwDsx8n/UQFZqrQTRGVsos9ASoFaFSKcg2BCpxapWoi32Bfmtz7Av7j4onzms7hgy3zna8XlECk94vIpw3Z/9pATvyxpfw7qZsiKmUXewLUVqCTMeaiogpkn+sC/Opsx5RyhJuXJ42G3k34Y9tJavo9+w7349zZC+mSvGufIyp0KV+MfJWFNyxkz1d7ym+b+uees9ZJ5eXjYx1XShViT4CaA3gCi40xTQqezN46fTHWxn9zXNE5pZxRv31XGj+yDOkfx75TD3P8aCgASUcC+XvbRez5cg8Lr1/ItGbT+Hnqes4mHi7bDum6KaXsYs8kCQMsA64HzgPfY+1AC9AMuBIrddLXInK967taOekkiapDMrM4/OUH7Fn2O+ve9c+XuOvSLtvo/89lHEnqid9VEwns0LPiOlqQ7vSrqiGXbliY3aAX8BrwTwrn8TsPzAIeEpFzdva1ytIAVTUd/+s4G2ds5Nf3f+XsibPc/vS7NGpyMPd8QnxTMhr+kwaDJuDuXbviOqrplFQ15fIAlafhMKA3EJV9aD+wWkTi7W6sitMAVbWlp6Xzx4INRKSNIjj0YKHzqSf9OXb+RoL7P4lvVKGR7bKnU9NVNVVmAUpdoAGqepCsLBJWLyZ9yys0bPAT7h75t8vIPO/GkcSOeLZ9kJCeQ8pvTZVu6aGqKU0Wq1QpGTc3wq65kcjHf+B0t53sOzqGUyl1cs+7e2QREb6BevtHMrvzG2x5dwsZpzPKvmM6NV3VYHYHKGNMM2PMdGPMH8aYE9mvP4wx04wxzcuik0qVJ/+YpjS+fzY+tyWwn5c5cqhZ7rnf1rch7pcUlv1zGa82fJUVD6/g+J6jZdcZnZquajB7J0mMAWYAXlg75xaUDowTkbku6V0VoEN8NcOxn9aQuv5FVrzdhPi/6+c7N3Dc54RenEnWReNp0O9O129Tr7P4VDXk0mdQxpgOwI9Yd12fA+8Df2MFqsbA7VgJZTOBriKy0fGuVx0aoGqWM8ln+HX2r2ycvpHkvcn4+Kfx4Juv5m4DcvJEAMe5ieABE/GNiKm4jmpQU5WcqwPUImAIMEpEPiyizC3AAuBjERlmZ3+rJA1QNZNkCX+t+IuEpW/SteubmAKD5Znn3Tmc2Amvdg8Q0v3G8k1Uq1PTVRXg6gB1GDgoIleUUO5nIFJEGpS6p1WYBiiVsuM3jq+YQojPMpvb1CcdbUhqnVtpMPgRvALq22jBxXRquqoCXD2Lrz6wuxTl9gCBJZZSqpqoe0kbGj/4Id6jEtif+RwJh/Onq6wffIgo7ynsf/5Kvrr/K479eaxsO6RZ01U1YU+ASgaKTBSbR0x2WaVqFA9fX6JufYLQR/4iMWoFsYf6kH7OM/f85tVt+OWNX5jeYjrzes9jxyc7yEw/7/qO6NR0VU3YE6B+BK4wxgwsqoAxpj/QCfjB2Y4pVZWFdL2W6Ee/JvOGA+w79RAH/m7G7q0XMlHs+3Yfnw7/iKTXo4h9YwQnd29z3cV1arqqJux5BtUNWIc1S28+MBfYh5VyMwb4BzAKK5t5dxGpEUFKn0Gp0pAsYd+3+9j41kZ2LdmFZAktOuzg5gcWZZ83HD7SBtP0LsKuvx03T88SWiyBzuJTlVxZJIu9F3gV23deBit4PSgi0+zpaFWmAUrZ6+TBk2yetZkGpx+iedvCW6edPBHAcRlCcL+J+EZdXD6d0oCmylmZ5OIzxrQFHgCuAsKxAtMhrLur10Vkq2PdrZo0QClHZZ47x+Gls3CLfYeGDX8vfP68G0cSO+DeagJh147AuLmXTUd0WrqqAJosthxogFKucGL7ZpJXvUio73J8/ApPVV/zxXBqdbiPNv9og099HxstOEGnpasKoAGqHGiAUq50Pi2Nw4un4XXkfcLCrVUdGekevHrPQ5xN88Hd251Lhl5Cu7HtiOzWyDULgDVjuqoAms1cqSrGw9eXyJGPE/bILo41WUvskf789n0HzqZZd02Z5zL5fcHvrBg9mRPTGxD79r2cPuLk+iadlq4qqYK74uYyxrzjRLsiIuOcqK9UjRfUoTtBHboTnpqO6fg7W97ZwuFNhwFo12sz9eonUo9pnF85g7jETni2mUDo1cPsv6t67jnbz6B0WrqqYEUO8RljnLm3FxEpoye6lYsO8anydGTLEbbM+oUrW/2DOvVOFjp/IimYE57DCBnwCD7hUTZaKIIjs/h05p9yglPPoIwxdzhzcRF5z5n6VYUGKFUR0k8kc2TxVGolfUBog32Fzmeed+dwQkc8Lh1P2DW3YNxd/PeizvxTTtJJEuVAA5SqaEc3fEPaD1MJq/cNtWqfK3R+3utP0nhQHy4bfRn+4f6uuajO/FNO0gBVDjRAqcoiPSWZI4tfp9axDwhtsBeA/X9GMmfS7QAYd0OT65vQ7o4WXHxdC9y8vB2/mM78U04qTYAqcpJECQ37Ae2BYOCAiPzsSDtKKdfxqluPqNFPA09z7OdvSf3+NX79vk7ueckUdi/bTajbu4Qf2cixczcQ0PshAlp1sP9ikZG276B05p9yIbum+xhj/LNn9x0DVgMfAePynB9vjDlgjCl2zyilVNkK6tiL6IeXccPyOQyeP5joHtHWCZNF2x5b8KtziujgjwjYdgXxLzfnwILJZKQWnnRRJE1Iq8pBqQOUMcYHWAvcCZwEVmGlOcprJRABDHZR/5RSTvCs7Unrka0ZvWY09+y+h97/boKnd/4huLDwXUSaJ8n8KJT9r/Tj6I8rbA/f5TVypDUhIirKGtaLiirdBIkFC6znV25u1vuCBU79+1T1Zk8286eAZ4APgbEikpY9FX2OiNyep9wu4ISIdCyLDlc2+gxKVTWZ585xZPn78Nd7hDfYipt74WdGSYkNOeUzlJBBT+MTEuCaC+vMP5WHqzNJ3AwcAe4QkcLJwi7YDzS0o12lVDly9/Ym4sbxRDy2ibSuO9iXPJYTSUH5ytQPOUSdtHm8EvEGi4YsYvcXu8k67+TkhyefzB+cwPr85JPOtauqLXsmSVwErBCRsyWUOwYElVBGKVUJ+Mc0w3/C20jmWxxZ+T8ytr9Ng6ANeHpn8Nt3l5GVIez8bCc7P9uJX5gf3cZ503RgW+q17Wr/xXQremUnewJUBlCaeakRQKpj3VFKVQTj7k6D60bAdSM4eyyBw0vfJCG1LnDhjic1PpWwrDeot/MAiWtiOBMwnNAB91ErKLR0F9GZf8pO9gzx7QbaGmOKDFLGmACgDbDd2Y4ppSpGraBQom7/P4avepS7d9xNl0e74BfmR2BoElHNrbudkLC9RNWajPuyRhx4qTuHl88lKyOj+IZ15p+ykz0B6lMgFJhcTJn/A/yAj53plFKqcghuEcw1L17Dg3EP0v/t64g72JHM8xfSJnl6ZxDZ8DvCU8Zw+t0gYqcO4/iW72035sjMP531V6PZM4vPF9gENAW+xwpYU4E1WOuhbgJ6A38AHUSkcM6Vakhn8amaJu3Qfo5+MRW/058QFHrQZpmD+1twyO99Wg2/FN9gX8cupLP+qjWXpzoyxjTCCkztAcFaB5XTgAF+BQaKSJxDPa6CNECpmkpEOLZhFWk/zSDYdxW+/hcm925a3Y7l7/fHzcONJjc0oc0/2tDkhiZ4eNvx2Fvz/VVrZZaLzxjTD7geiAHcgTjgK+BTEalRibg0QCkFmefOcuTLOchfc2kQspG5z43m4J78kx8G3PU19S6uj0/7cQRfeUPJ+1Zpvr9qrcxy8YnIF8AXDvVKKVXtuHvXImLwXcBdnEk4TJszh2DeNg5usIYAvX3OcGnHjXh4ZcKh5SRPCyHF3EBArwkEtGxnu1Gd9VfjFfknjDHmE2PMdcaYgumMlFKqSLVDw2l/Vwfu+PEO7tl1D1f++0raXnvACk7Z6gUlEl1/NgG/tSfhlYuJffdR0g7F5m/I0Vl/OrGi2ihpR13Byh4xFyul0Z5y7FuVoEN8SpVMMjOJ/+Zj0n9/l7DA9XjXSi9UJvO8G/EJrTnfcBThg+7D08fT/l17dWJFleHsjrrTgWFAIBcmQnwPvA98LCKnbVasYTRAKWWfjJMpHFk+C/eDCwkL/Q13j/zPk/7c1IzPZ42m+eDmtB7Vmsa9GuPmUcoVMTqxospwepKEMcYLGAjcDlyNNSFCgDTgf8BsEfnRZT2ugjRAKeW404cPkPjldGqnfEZog78A+Pj1m9jxS8vcMr6hvvQen0LDnh0J7nZ98ZMrdGJFleF0slgRSReRj0XkOiAKeBIro4QfcAew3hiz0xjzqDEmzFUdV0rVDD7hkUTf+QKhD+8h5bIt7Dt+J8fSOuUrc+bYSZqGvkrIwf6cfCuEfa+PJGnTOtsNFjWBQidWVEmOTjPvjHVXdRNQB+uuKhNrqvn7wBcikll0C9WH3kEp5VoiwpEtR9g2fxvbP9xOg7CtjHxsYaFySYnhnPS4nno977owE1CfQVUZZbYOKs8FagNDgTFAjzynjopIjbij0gClVNnJOp/FoZXLkT+mEhL4A7Vq205QkxgfzWnfAdS/Zjz+GzbbN7EC7J+MoZxW5gGqwMWuAeYDwYCIiHsJVaoFDVBKlY/zaanEfzUH2beQsKBNeHoXTk4btzuC1Wsn0Wp4Ky4Zegk+QT42WipA77oqRHncQflhzfQbA3ThwhbwB0Qk2uGGqxANUEqVv/QTx4n/8h3cDy0iLHQb7h7WE4Wv5vXllxXWMyzjbojpHUO7mw1R11+FT4NGthvTmX8VoswySRhjegK3ATcCtbEC0zlgKdYzqJWOtKuUUqXhFRBI5IiJwETOxB8i8euZeB37jJ2/tMotI5nC3yv/ou8106j1TTKH41uSETSI4L7/zB+sdCPFSsuebOaNgdHZr0gu3C39CswG5otIcll0sjLTOyilKo/U+FT++PgPtn+4nYMbDhLSKJ7xU2bmK5N53o2EhJZkBA8muM8/8encTe+gKoDT08yNMT7GmNHGmDXAHuAprOnmJ4DpwOUicrmIvFkTg5NSqnLxC/Oj470duePHO3jgwAN0e6wtifEx+cq4e2QR3vB3oryexXtVFMc7Z5Hl5Zm/oZJSKmk6pXJRXCaJ97Cmkfti3S1lAauxhvA+F5HCuUpqIL2DUqryO7lnO8e/fQef1OWENNhbuMAPkLnQDbcTWaQHhiGTnqPW3bfbbkwnVbiEs6mOcpZd7wPmYGWNsL07WQ2mAUqpquXk7uxgdfoLQsL25R5fv6Qb3y66GgDjZojsFkm7m92J6tOOOhdfcqEBnVThEs5OklgAvC8ia1zbLaWUqjh1mraiTtM3gDdI2bWN5DWz8Dm9nD9+vpBeSbKE/d/tp2f396nzywESlzbmdK2+BFx1OwE6qaLcFPkMSkRu1eCklKrO6jZrTfRdbxLy0F5u+e4F+kztQ+SVkWDAt+4pIptaQSckbB/RATMI2NaBzIAifm2WlE5Jn1vZzaFp5kopVd3UjaxLp/s70en+TqQmpBK79EuOxLciNGRHvozr7rdkwrtAnqfwWd7emEmTKHLzvILPrfbvtz6DPrcqhssySdRU+gxKqertTPwhEle+h3vCUsJCfsPD8zz8ACwCjgFBcLafJ9O//j+a9GtO80HNiekdg0etPH//63OrQso11VFNpQFKqZrjXPIxElfMhbhPCa2/Ga/sjRe3rr2MpbMG5Zbz9PWk08h0LrrSi5BrxlC7QYRuA1JAmWWSUEqpmsi7XhCNhj8MPExG6kkOrfqQ839/yl9/RuUrl5GWQWitj4ly20HWyqc4X9cdjxPnCzeo24AUq5TbVCqllMrL068ODQePI+qRlQxZ8za3/3A7XR7tQmCTQNw9M2jSZg8Abu5ZeIw4D17562d5eyOTJhV/kRo+sUKH+JykQ3xKqbxEhKQ/Ykn97gV8zq68sNaqwHMrboa0S/04dqoj6Rf9h6hrO+HllyeKVfMFwfoMqhxogFJKFSc1djfHvp2DZ9KXhIZsx8Mz/16u5zPceXHcY2SJD9E9o2narylN+zUhoEfbaj2xQgNUOdAApZQqrXPJxzj6zUKyYpdQv87P+Pqn8ddvF7HgxVvzlYu59C9G/T7f9rT1ajKxwulksUoppVzHu14QETfdR+Sjq/G54wSJ4Us4UecBQluH5ivXrO1uTJDtNrLCw4u+QDV7ZqWz+JRSqgIYDw9CegwgpAe0fxJSDqSwe/ludi/bTUzrveBHoQXBeIHpf4jEV2M47dkDn8uGENylL8bdvVouBtYhPifpEJ9SytXSU46T+M1CvD+dSf0VO3A7LrkTK+iav+zXHw3hTMBIbvhqAl7HjhRurJI+s9J1UEopVQV51Q0kYsg9MOQeJDOToz9/Q9rWT6h9di3BWX/j5nbhxuLPDY1IObaNQdgITlClk9hqgFJKqUrMuLsT3KUPwV36AFbqpWNrFiIHv4TTcaQcCwAghboEkFKo/vm6Xhye9yyBVw7Hr3HTcu27s2r8EJ8xZgQwHmgNuAN/Ym1hP0NESpwqo0N8SqmKIlnCkS1H2PPlHrLmzafr33PxIuNCAS/gTnKHBZMSwzklnfFq2p+QHjfi4etvPbt68knrTisy0tpJuByeWek08xIYY6YDdwNnsXYLzgB6A/7A58BNIpJZdAsaoJRSlce5mbNxe+oRPI4dRwINbsOl0DOrHBnpHiQuvoSwL3fhnnHuwolyWgysAaoYxpghwCdAPHCViOzJPh4KrAFaAA+IyOvFtaMBSilVGUlmrwARGAAAEShJREFUJkkb15C65TO8T68jOHhXoUXCZ+/0ptaZc4XqZjUMx+3goTLtnwaoYhhjNgHtgNEiMq/Aue7AWqzg1bC4oT4NUEqpqiA9JZmjaz8h4+8vqOvxE/WCEpGRFLmH1fGpYZzMvAKPxn0Jvmoo3vWDXdofDVBFMMZEAHFYKwwCROSMjTIHgYZAVxH5sai2NEAppaqilF3bqNWlD97H4wufDALyjB1lZbpxNLExZzy7UKtFP0IOpeH29DNOPbfSaeZFa5v9/oet4JRtI1aAagsUGaCUUqoqqtusNbzxcqGEtFmeBrkR3Llw8+LmnkVog7+Bv2HxB0jeBcRluCC4pqY6apz9biMTY66cxQONiymjlFJV18iR1oSIqCgrx19UFG6zPyDrpWQO15lNbMJNHE2IJN9DjkVg0gu0c/q0NRPQxWrqHZRf9ntaMWVSs9/9C54wxowFxgJE6oZjSqmqbOTIQnc+nkB4vzHAGABOHz5A0vpPyIpbSeSxFbafW5XBguCaegeV8/069ABORN4RkfYi0j442LUPDpVSqrLxCY+k0bCHiHrka0xUlO1CZfDHek0NUKey3/2KKZNz7lQxZZRSqmZ57jlrrVRePj7WcRerqQEqNvu9iD8FAGhUoKxSSikbz63KamFvTX0GtTX7vaUxpnYRM/k6FCirlFIKbD63Kgs18g5KROKALViZqm4qeD57oW4E1kLdDeXbO6WUUlBDA1S257PfXzDGXJxz0BgTAryV/XFKaRLGKqWUcr2aOsSHiHxijJmBlcn8d2PMN1xIFlsHWAxMq8AuKqVUjVZjAxSAiNxtjPkemAB058J2G+9Tyu02lFJKlY0aHaAARGQhsLCi+6GUUiq/Gpks1pWMMUcpPmVScYKAYy7sTk2g35l99Puyj35f9nHm+4oSkWIzHWiAqkDGmE0lZfNV+el3Zh/9vuyj35d9yvr7qsmz+JRSSlViGqCUUkpVShqgKtY7Fd2BKki/M/vo92Uf/b7sU6bflz6D+v/2zj3a6rLM45+viKB4zRsZVppoV1NRmnBU1KysGG0gp4sB5bha4Jrsog52sXQaw8zKDENTO2owNeh4nemykjBLJQ/HRZRaIOAgypioqGhY8swfz7vj58+999nn7N8++7fPeT5rvetlv8/7ey8Pe/+e816fIAiCoJTECCoIgiAoJWGggiAIglISBqogJH1Y0h2SNkh6VlK3pNMkNaxjScMlHSvpIkl3S3pU0guS1kq6TtLEFnZhwClCZ3XKPl+SpXBGEe1tN0XrS9K2ks6SdI+kpyQ9J2mVpAWSDi+6/QNNkfqSNEbSJZL+IOl5SX+WtFzSXEn7tqL9A4WkAySdLukHkh6QtDn9bqY0WW7z+jezCE0GYA7unfd54FbgBuDplPZfwLAGy3lHesaAR1NZPwKWZdLPa3d/y6SzGmUfBvwV2JzKO6Pd/S2bvoB9gOXp+f8DbgL+E/gN8ALwhXb3uSz6Ag4GnkzPrsHv6bwReDilPQNMaHefm9DVtzLvl2yY0m79t105nR6AyRmDMjaTvidwX5Kd3mBZxwDXAUdUkf1TeukacHS7+10WnVUpewTwe2Bt+lF0vIEqWl/AKGBF5Q8eYHhOviuwf7v7XSJ93ZmeuTyrK2A4cGWSLW13v5vQ1z8DXwNOAl4HLGrGQBX6Tmy3cjo9AN1J4VOryI7K/EdtVUBdV6Tyrmx3v8uqM+CC9PwkoGuQGKhC9YW7mjHg6nb3rez6AkayZUQxuop8r4x8u3b3vSD9NWugCtN/rEE1gaQxwDh8SmRBXm5mt+N/yY8G/q6AKivefccUUFZbaKXOJL0N+Cww38xuab617adofUnaBjg1fZxdXEvLQQu+Xy/iMxcAqiKvnNPZiE9nDWmK1n8YqOY4OMW/t+pu4wHuyeVthrEpfrSAstpFS3QmaSRwNfAEcHr/m1c6itbXOHwKb42Z3S9pQtpQcpmkcyW9vdkGt5lC9WVmfwFuSx/PlTS8Ikv//kr6eKWlIcIQp1D9D3l3G02yT4rr3Wb+v7m8/ULSaGB6+nh9M2W1mVbp7N+BA4APmtlguo26aH29JcXLJXUB03LycyRdD3y0zgumzLTi+zUT+Ak+8jxeUndKPwzYBbgYOLOP7RysFKr/MFDNsX2KN9bJ82yKd+hvJZK2Bn4A7ATc1uHTV4XrTNIE4FPAjWb2oybaVkaK1tcrUnwk7qDz68BcYH1KuxRf5H4a+HhfG1sCCv9+mdnK9B27Bjiel06xdwO/TCOtoGD9xxRfc1TmpFs9tJ+Lu6JfA5zc4rpaTaE6k7Qt8H38hTqziDJLRtHfscpvfmt8WupMM3vQzJ4ys5uBE1Nd0zr0fE/hv8lknH4H7AecgPtA2h3X1S7A9ZLOKaq+DqdQ/YeBao5nUrx9nTwV2TN18tRE0sXAKcA64FgzW9efckpE0To7H9gf+IyZdfLaXC2K1lc2z/fyQjPrBpbg74aJDZRXNgrVl6Sd8TNPOwDvNrObzWy9mT1uZjcB78Y3R3xR0th6ZQ0RCtV/GKjmWJ3i19TJs3cub8NIugj4JPAn3Dgt72sZJWR1iovS2fvxA7nTJC3KBvzlATAjpV3Rj/a2m9UpLkpf2TyrauSppI9uoLyysTrFRenrvfho6W4zW5kXmtkKYDE+Ip3YaCMHMatTXIj+Yw2qOSrbvt8kadsai8qH5fI2hKSvAZ/B1waOM7P7+t/MUtEKnW2Fn6+oxb4p7NxgeWWiaH31ZP69K/7HT57dUvxsFVnZKVpfr07xhjp5nkrxK+rkGSoUqv8YQTWBma3Bf/DbAB/IyyUdhS+orgPuarRcSbPxXUFP4sZpaSENLgFF68zMXmtmqhbwbecAZ6a0g4rrycDQAn2txf/iB1/XzJe3C3BI+tidl5edFvwmH0nxuOwW80x5w/Gt+1B7RDpkKFz/7T613OkBmMKWk9H7ZdL3wK/cedm1HvhJ/geAr1Yp79/SM08C49rdv07QWZ16uhgcN0kU/R2bxJY7+A7KpI8Efphk3SR/cZ0WitRXemZjeuY7wIiMbATw3SR7Atip3X0vSH+L6OUmiV6+X33Wf8162q2MwRDwrbmGL5begl+GuCGl3UDuYsTMi7Mrl/4PbLk25Z6Ur1qY1e4+l0VnvdQxKAxUK/QFXJjkm4BfpjLWprSHydyh1omhSH3hZ8Uq92CuBW5OZT6S0v4MnNjuPjehq0OAuzOhcqnrH7Ppffx+9Un/tUKsQRWAmc2U9CvgNHwtZBj+18VVwHfNbHODRWXnsA9NoRq30+HX1BSosyFB0foyszMl3Qn8C36ifzv8AOU3gNlmVm1tqmMoUl9mdrWkZfhZuyOAdybRWvyy2G9YZ68R7wi8rUp6v3clFqX/cPkeBEEQlJLYJBEEQRCUkjBQQRAEQSkJAxUEQRCUkjBQQRAEQSkJAxUEQRCUkjBQQRAEQSkJAxUEQRCUkjBQQVADSaslWQpf7SXvvEzeRQPUxJYiaZikZZIekjQiJxsv6Q5Jz0t6TNKlkkbVKWeJpAeT/65qefaS9Jyk61rRl6AzCQMVBI0xVdKwagJJO+JuPwYbM4A3A182s02VREmvAhbitw/8HL9zbQawoEY5n8Sv05lhNdzIm9kjJG++kiYW1YGgswkDFQS90w3sBRxXQ/5BYFv8/sRBgaTtgXPxG7qvyYnPAkYBHzOzSbjxuQ04XtJhuXL2Bs4D5pnZz3qpdjZ+N+DXm+9BMBgIAxUEvdOV4uk15NOBF4FrB6AtA8U0/G7ILjN7MSc7BL8g9T8AkvyqJHt7Lu8lwAu4b7O6mNnj+MWi4yQd3v+mB4OFMFBB0DuLgfuAE5IL8L8h6QD8pfxTfKqrKpLeIWmOpKWS1kvalNZ2rpb0hhrPjJQ0S1KPpGfTM49KukvSVySNzOUfL2mBpLWS/iJpg6QVkuZLOqaPfZ6Z4vzoCdzR4YbchZ9PpPhvbZL0j8AJuD+uxxqst+LDa2bdXMGQIAxUEDRGF/7y/VAufXqKv9/L83OBU3C3DXcA/4OPLKYC3ZL+PptZ0lbAf+N+d/bFb7C/HjeUewOfJ+MhWNJxwK9wXzyP4S4NFuJ+xaYAJzXYTySNBd4IrDCz1VWyrAZ2l5S9ff/1KV6VytgB+DbuyqM33WT5BT4afV+tNb9gCNFuXyQRIpQ14C9iw92ejMaNy+KMfBjucmE97kG04qhtUZWyTgR2zqUJ+ER65j4yDgKBI1P6EmBUlecOB7bLpC1M+T9Upe5d6YPzS+DUVNY1NeQzkvzaVPbBwBrcj9DuKc+38fWk1/dD7/em8se3+zsQob0hRlBB0ABmtg74CTA+MyX3TnzzxHwze6GX5280s6dyaWZmlwF3Am/ARy0V9kzxHWa2scpzvzaz56rk/3GVuteb2ZL6PXwJB6X4/hry7+FO7E4GHsddfI8BzjKzP0k6FPcDdL6ZPVB5SNKIBkdFFd9KB/ehzcEgJAxUEDROV4qn5+IuGkDSGEmfkPRNSVdK6pLUhY/OAPbPZO/Bp7pOkTRT0p758nL8JsXzJR3e5PTYHileX01oZn/FndCdClwGXARMMLO5qd7LgeUkp5qS3iVpKb6xYpOkn0rar079lfWs3vocDHLCo24QNM7N+Ev7o5IuxDcALGtkdCLpXOBz1P/N7Vj5h5k9KOnT+JbrOcAcSSvx0dZNwA320t11Z+Mjn+NT2ChpCT71d62ZrWy8m+yU4qdrZUgjxitSyPKp1I6jzWxT2nZ+Kz4qOgnYBV9XWyjpTWb2TJXiK/XuXEUWDCFiBBUEDZJeyvOBV+IL/yNoYAOApMnAOcDz+Kjjdfj6kcxMpO3a+NpStr5LgNfgaz7z8DWvk/EDsd3pgHAl7zpgHHAsPnLpwQ/Sfhn4g6SP96GrlanIHevmyiHp1fjZqavM7PaU/FncKE82swVmdjnwr/hGjw/XKKpS75N9qT8YfISBCoK+0ZXi9+GbJuY18MwHUvw5M7vCzFbaS29UqDndZWbrzGyumZ1sZq/FRyfLUjwrl3ezmS00s7PN7Eh8A8Ms3EDMyRq0XqhsCd+1wfwVvgNsxA/yVngr8LiZrcik3ZWRVaNSb6Nb04NBShioIOgDZtaDb+deDyywxs73VLZjr8kL0oaLhjcDmNlS4OL0sdYLvpJ3o5ldADyMb5E/oMFqelL8xrq5MkiaAkwCPm1mT2REG/FbNrJU7uyzGsVV6u2pIQ+GCGGggqCPmNkRZrabmdWaospT2cl2qqRtKomS9sAPpr5sXUrSMZLeI2nrXPow4D3p40OZ9DPStUL5cg7FpyQ344aqEX6R4vytEFVJI7OLgZ+Z2fyc+LfAKEnZc1jTU3xvlbK2x+//e5owUEOe2CQRBK3nW/iB3PcCKyQtxkcVR+Gjqhvxc1JZDgS+CWyQ1IPfUrEdvq70SmAdcEEm/xeACyXdj28P34Sv80zA/xCdbWY1b7rIYmarJP0WOFDSPma2qpdHzsc3P8yoIpsNfASYJ2lqyjcBWMmWtbcsR+Nrbbfay69YCoYYMYIKghaTdtAdAvwQ3wgxCT/3dDk+StlQ5bFb8A0HPfga1WTgCNwwfQk40MweyuQ/DR+NbcZf8u8HXpXKeZeZnd3HZl+a4qn1Mkkajxum86rtFDSzP+LnxRbjGzjejG/ymJg/35WYlqs/GMLIrNY0cBAEQ5Xk2+khfKpt7ECMZiTthk9D/s7MDm11fUH5iRFUEAQvI41uvgTsQy+jqAKZhW/dP2OA6gtKToyggiCoStqQcS9+cHd/yzgtbEFdewErgB+b2eRW1RN0FmGggiAIglISU3xBEARBKQkDFQRBEJSSMFBBEARBKQkDFQRBEJSSMFBBEARBKQkDFQRBEJSS/wcZp8WYakG5aAAAAABJRU5ErkJggg==\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "plt.plot(num_sol_imp2[:,2]/m0, num_sol_imp2[:,1], label = \"Implicit\", color = \"purple\")\n", | |
| "plt.plot(num_sol_exp2[:,2]/m0, num_sol_exp2[:,1], \"--\", label = \"Explicit\", color = \"orange\")\n", | |
| "plt.plot(m/m0, -np.log(m/m0)*250, 'o',label = \"Tsiolkovsky\", color = \"red\") \n", | |
| "plt.xlabel(\"Mass (%)\")\n", | |
| "plt.ylabel(\"Velocity (m/s)\")\n", | |
| "plt.legend(loc=\"best\")" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "The velocity will not converge because Tsiolkovsky method doesn't take drag or gravity into account. " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 10, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "<matplotlib.legend.Legend at 0x7f749bbe0710>" | |
| ] | |
| }, | |
| "execution_count": 10, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEaCAYAAABEsMO+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3hUxfrA8e+kh4QQEjoBAgGE0EmQTkJoAb0oHWkBwQLYsaHeCzb0+hO5eFXUi4ooSpGONEF6EUIE6dKJID0QIIS0+f1xNpts+ia7abyf59lnc+bMzBlWs2/OOXPeUVprhBBCiOLGoagHIIQQQmRFApQQQohiSQKUEEKIYkkClBBCiGJJApQQQohiyamoB1DSVahQQfv7+xf1MIQQokTZs2fPFa11xZzqSIAqIH9/fyIjI4t6GEIIUaIopc7kVkcu8QkhhCiWil2AUkq5K6VeVkrtVkpdV0rFKaVOKaUWKKXaZ9NmiFJqi1LqhlLqllIqUik1XimV478vv+2EEELYX7G6xKeUqg2sBeoCl4BNwF3AH3gI2Adsy9DmU2AcEA+sBxKBLsAnQBel1ACtdXIWx8pXOyGEEIWj2AQopZQH8AsQALwNvK21Tky33xfwzdCmH0aQuQB00lofM5VXBjYAfYCngOm2aCeEEKLwFKdLWW9gBKfZWut/pQ9OAFrrq1rrPzO0mWh6fyU1yJjqXgTGmjZfzeKSXX7bCSGEKCTF4gtYKeUCPGbafD+PbfyAICABWJBxv9Z6E3AOqAK0KWg7IYQQhatYBCiMgOELRGutDyul2imlpiilvlBKvamUaptFmxam94Na6zvZ9Ls7Q92CtLOp+BvxrHlhDVeOXrHXIYQQokQrLvegmpjejymlZgERGfb/Sym1EBieLqjUNr3nNJf+bIa6BWlnM8fXHGfpyKXcunCLywcvM3T1UJRS9jiUEEKUWMXlDMrH9N4JGAF8iDGTrzzG7L1zQD/g03RtPE3vt3Po95bpvawN2pkppR43TUmPvHz5cg7dZM2joge3LxmHP7H2BIcXHba6DyGEKO2KyxlUaqB0AmZqrV9Kt2+ZUuo8sAuIUEq9o7U+CaSecli74mJ+25lprb8EvgQIDg62up+qLasS9GQQkZ8ZGSjWPL+GuuF1cfFwye+QSq2kpCSuXbvGjRs3SEpKKurhCCEycHJyoly5cvj4+ODkZNuQUlwC1M10P/8v406tdaRSag8QDIQCJ9O18cxYP53Ufen7z287mwp7J4xD8w8RdyWO2OhYtry7hS5TutjrcCVSSkoK0dHRuLq6UrNmTVxcXORSqBDFiNaahIQErl69SnR0NLVq1cLBwXYX5orLJb7T6X4+lU2d1PIqGdrUyqHfGln0n992NuVe3p2uH3Q1b2//cLtMmMggJiYGJycnqlatiqurqwQnIYoZpRSurq5UrVoVJycnYmJibNp/cQlQUel+9s2mTgXTe+r9od9N742UUu7ZtGmVoW5B2tlc84jm+LX1AyAlMYXVz6xG63xfeSx1bt26hbe3twQmIYo5pRTe3t7cvp3TrX3rFYsApbU+B/xm2sx0nUspVR5oadqMNLWJxghsLsCALNqEAH4Y2SJ2pDtWvtrZg3JQ9Pq0F8rB+AKWCROW4uPjKVOmTFEPQwiRB2XKlOHOneye3MmfYhGgTN41vf9LKdU8tVAp5QbMAMoBe7AMGu+Z3v+tlKqbrk0l4DPT5vta65QMx8pvO5ur2qIqwWODzdtrnltDwu0Eex+2REhJSbHp9WwhhP04ODiQkmLbr8xi89uvtV6OMb28EvCbUmqzUmoxcAIYhDHV/BGd7hqY1vonjOBVBdivlFqulFoEHAMCgSUYyV8zHitf7eyl89udKVPROFOI/SuWze9sLqxDF3tyeU+IksEev6vFJkABmKaX98XIWN4E6AXEAR8BLdLnzUvXZhwwFOOyXQjQAziOkey1X3YZyfPbzh7cy7vT7YNu5u0dU3fIhAkhxD2vuEwzN9NaLwYWW9nmB+CHfBwrX+3sodmIZkT9L4ro7dGkJKaw6ulVDFszTM4ghBD3rGJ1BnUvyzhh4uQvJzm8UCZMCCHuXRKgipEqzasQPC7dhInnZcKEyNmRI0d48sknue+++yhTpgzu7u7UrFmTdu3aMWHCBH755ZdMbfz9/VFKcfr06cIfcB5t3LgRpRShoaF2P9bkyZNRSlm8HB0d8fX1pWPHjkyfPp2EhKL9PezQoQNKKbZu3Vqk4yhsEqCKmbC3w2TChMiTefPm0axZM7744gtu375NaGgoffv2pUGDBvz555989NFHTJw4MfeOBAABAQFEREQQERHBoEGDCAgIYOvWrTz33HOEhITYfAp1cTRz5kyUUowZM6aohwIUw3tQ9zo3bze6/V83lo5cChgTJpqPbE6F+yrk0lLcSy5cuMCjjz5KQkICH330Ec888wyOjo7m/SkpKWzdujXLv7jXr19PYmIi1atXL8whF3sdOnRg1qxZFmW//vor4eHh7Ny5k//+97+8/PLLRTO4e5ScQRVDzYY3o0Y7I9tS6oQJyTAh0luxYgVxcXG0bduW559/3iI4gfFMSqdOnXjttdcytQ0ICKBBgwY4OzsX1nBLrLCwMCIijNV/VqxYUcSjufdIgCqGZMKEyM2lS5cAqFSpktVts7sHFRoailKKjRs3sm3bNsLDwylfvjzlypWjR48e7N2711x39uzZtGrVCk9PT3x8fBg2bBgXLlzIdKxZs2ahlGLkyJFcuXKFsWPH4ufnh5ubGwEBAbzxxhvExcVZ/W+4evUqb7zxBk2aNMHT0xMPDw9atmzJtGnTSExMtLq/nDRr1gyAixcvZltn+fLlhIeH4+vri4uLCzVr1mTUqFEcPXo02zYJCQl8/vnnhIaG4uPjg6urK7Vq1eIf//gHc+fOzfP4pk6dioODA9WqVbP4bwRGurD333+f4OBgvLy8cHd3p3Hjxrz11luZ0hL5+fnx2GPGwuZfffWVxT25orrkJwGqmKrSvAqtxrcyb695fg0Jt2TChDDUrFkTMC7XHThwwKZ9L1++nJCQEGJiYujRowfVqlVj7dq1hISEcOzYMV588UXGjBmDt7c3PXr0wMXFhTlz5tC1a9dsJxPExMTQunVr5s+fT+vWrenRoweXL1/m3XffpUuXLlYFqf3799O0aVPeffddrl+/TmhoKCEhIZw5c4YXXniBnj172nRSw40bNwCoXLlylvtfeuklevfuzS+//ELjxo3p168fnp6ezJo1ixYtWrB69epMba5evUrHjh0ZO3Ysv/32Gy1btqRv3774+/uzZcsW3njjjVzHlZKSwjPPPMOLL75Iw4YN2blzJ82bm5PwcPbsWVq1asXEiROJjo6mXbt2dO/enatXrzJp0iQ6dOhg/rcBDBw4kHbt2gFQr1498/24iIgI2rdvb9VnZityD6oY6/xWZw7OO8jtS7fNEya6vt8194b3gDfVm0U9hHybpCcVuI+HHnqIatWqcf78eVq0aEH37t0JCQkhKCiI4OBgypUrl+++p02bxvz58+nfvz9gfBEOGzaMH3/8kb59+3L58mX27t1LYGAgANeuXaNt27YcPHiQefPmMXz48Ex9Llu2jPbt27Nnzx68vb0B44ykW7du7Ny5k8mTJ/PBBx/kOrY7d+7w0EMPcf78eaZMmcJLL71kXoPo2rVrDBo0iHXr1jFlyhQmT56c788gldaan3/+GYAHH3wwy3/Xhx9+iKenJ6tXr7b4In/vvfd47bXXGDJkCH/++ScVKqTdRx4xYgS7du2iQ4cOzJ8/n6pVq5r3xcfHs2HDhhzHdefOHYYMGcKSJUsICQlhyZIl5s81ddz9+/fnyJEjPPvss7z33nu4uxu5sePi4hgzZgw//vgjEyZMYObMmQB89NFHzJw5k+3bt9OpUydzeVGSM6hiLHXCRKodH+3gyhHJMCGgbNmyrFu3juDgYJKSkli5ciWvvPIKXbt2xcfHh/bt2zNv3rx89T148GBzcALjflbq5IADBw7w1ltvmYMTgI+PD08++SRAtl+sSilmzJhh8SVauXJlpk+fDsDnn39OfHx8rmObNWsWp06dYuDAgUycONFigTwfHx++/fZbnJ2d+fTTTwt03zYhIYFDhw4xYsQIduzYQVhYGE899VSmelOnTgXghRdeyHSWMXHiRIKCgoiJieGrr74yl0dGRrJy5Uq8vLxYsmSJRXACcHNzo2fPntmO7cqVK4SFhbFkyRIeeeQR1q5da/G5gnG/bPfu3bRv355p06aZgxMYSV2//PJLKlSowOzZs4mNjc37B1PIJEAVc02HN6VGe5kwITJr2LAhu3fvZtu2bbz22mt06dKF8uXLk5KSwvbt2xk8eDAjR460ut/w8PBMZXXr1s1xf7169QA4f/58ln02bdqUJk2aZCrv3Lkz1atX5+bNm+zZsyfXsa1cuRKAAQMyLUQAQLVq1ahXrx5Xrlzh2LFMmdFy9O2335rvubi6utKoUSO+//57Hn30UX755ZdMmfUTEhLYscPIXZ3d5zxq1CjAeK4rVeolvz59+uDrm93qQlk7fvw4bdu2ZefOnbz00kvMmTMHF5fMK3Gnfk79+/fPMhuNp6cnLVu2JDExkcjISKvGUJjkEl8xp5QxYeLLll+iUzQn153k4PyDNB7UuKiHVqRscZmstGjXrp353kFKSgo7d+7kzTffZO3atXz77bc88MAD2X6hZ8XPzy9TmaenZ572Z3cWVLt27WyP5+/vz7lz5/jrr79yHdvJkyeB7ANUepcvX6Z+/fq51ksVEBBAhw4dALh58yaRkZGcPXuWr7/+mmbNmvHMM89k6j8xMRFHR0dq1KiRVZcEBAQAcO7cOXPZmTNnAGjQoEGex5bqscceIykpiaeffjrHS6Kpn9Pzzz/P888/n2Ofly9ftnochUUCVAlQpVkVWj3Vil0f7wJg1dOrqNOlDmUqyFpJwpKDgwPt2rVj5cqV3H///URFRbFkyRKrAlRuS5zYawmUvOSdTE42cjg/8MADFvd0smLt2UnG56CSk5N5/fXX+fe//82ECRPo1KmTxSSE1CsZqWddWbH11Y5hw4Yxe/ZsvvnmG/r27Zttpo3Uzyk0NJRatXJaPDxtwk1xJAGqhOj8VmeOLDpC7F+xxF2OY9XTq+j3Y7+iHpYophwdHQkLCyMqKqpY/IWcU1ql1H3VqlXLtZ8aNWpw9OhRxo4dywMPPGCj0WXN0dGR9957jx07drB582ZefPFF1q1bZ95fqVIlnJ2dSUxM5OzZs1meJZ46dQrA4qHo1ICR0xT07IwePZqwsDBGjRpFr169WLx4MT169MhUL/WMbvDgwTzxxBNWH6e4kHtQJYRbOTce/DJtFtGBuQdk9d17WF7+Mj979iyQ9SW5wrZv374sp8Nv2rSJc+fO4enpSVBQUK79pE4eWLBggc3HmBWlFNOmTUMpxfr16y3uJbm4uNC2bVvAeC4sK6lnZOnPdFIDyuLFi7l27ZrVYxo+fDhz584lKSmJ3r17s3Tp0kx18vs5pd7PSkpKsnpc9iABqgSp17MezUelXWL4eezPxF21/iFHUfJ99tlnjBo1il27dmXal5SUxP/+9z9++uknAAYNGlTYw8tEa824ceMsnru5fPkyzz77LACPP/64xUyz7Dz++OPUqFGDb7/9lkmTJmX5/NSBAwf45ptvbDb2li1bmj/DSZMs732+8MILgDFFe+fOnRb7PvjgA3bv3o23tzejR482l7dq1YqePXty48YN+vTpk+kB4Pj4+CyfnUqvf//+LF68GAcHB/r378/8+fMt9vfr14/mzZuzfv16xo8fT0xMTKY+Tpw4wWeffWZRlnqmd/hwMfnjV2strwK8goKCdGG6E3NHT602VU9msp7MZL1w6MJCPX5hOnToUFEPodiaNm2aBjSgq1SposPDw/WQIUN0eHi4rlatmnnfyy+/nKltrVq1NKBPnTplUR4SEqIBvWHDhiyPmdpnVjZs2KABHRISYlH+zTffaED37t1b16lTR/v6+up+/frphx56SHt5eWlAt2rVSt+6dStP/Wmt9R9//KFr1qypAe3j46NDQ0P14MGDdefOnbW/v78GdOvWrbP97DKaNGmSBnRERES2dY4fP66dnZ01oNevX2+x78UXX9SAdnR01KGhofqRRx7RgYGBGtDu7u565cqVmfq7fPmyDgoK0oB2c3PTXbt21Y888ogOCQnR3t7eOiAgwKJ++/btNaC3bNliUb5u3TpdpkwZ7ejoqGfNmmWx78yZM7pRo0Ya0GXLltUdOnTQgwcP1l27dtX16tXTgK5evbpFmzt37uhKlSppQAcHB+uIiAg9evToTH1nx5rfWSBS5/L9WuRf8CX9VdgBSmutj644ag5Qk5msDy85XOhjKAwSoLIXGxurFy9erMePH69btWqlq1evrp2dnXWZMmV0/fr1dURERKYvs1RFEaAiIiL0pUuX9JgxY3S1atW0i4uLrl27tn7ttdcyBaec+kt1/fp1PWXKFN26dWvt5eWlXVxcdPXq1XWbNm30P//5T71v374s22UlLwFKa63HjRunAd2hQ4dM+5YuXaq7d++uy5cvr52dnbWfn5+OiIjQhw9n/7sZHx+vp0+frtu2bau9vLy0q6urrlmzpu7du7eeN2+eRd3sApTWWm/ZskV7eXlppZSeMWOGxb64uDj98ccf644dO5rHVrVqVR0cHKxfeuklvWPHjkz9RUVF6Z49e2ofHx/t4OCgAT169OgcP5tUtg5Qyqgn8is4OFgXxXMESyKWsG/2PgA8q3gy7uA43H1yv0RSkhw+fJiGDRsW9TBEAcyaNYtRo0YRERGRKVO4KH2s+Z1VSu3RWgfnVEfuQZVQPab1wLOK8ezJrQu3WP1czteshRCipJEAVUK5+7jz4Bdps/r++O4P/lzxZxGOSAghbEsCVAl2X+/7aDI0LX3MiidWcCem9K/6KYS4N0iAKuHCp4fjUdkDgJvnb7L2hbVFPCIh0owcORKttdx/EvkiAaqEK+NbhgdmpD1Rv3fWXo6tsi5JphBCFEcSoEqBhn0a0nhwWvLY5Y8tJ/5G7ksXCCFEcSYBqpTo+d+elKloJI+9ee4mayfIpT4hRMkmAaqUKFPB8lLf71/9zvE1x4twREIIUTASoEqRwH6BBA5IW+l0+Zjl3I29W4QjEkKI/JMAVcr0+qSXeZ2o2L9iWfuSXOoTQpRMEqBKGY9KHvT6tJd5O+rLKE6uO1mEIxJCiPyRAFUKBQ4IpGHftHxYy0Yvk0t9QogSRwJUKaSUotdnvczJY2+cvcGKJ1cgiYGFECWJBKhSyrOyp8WsvgM/HmDft/uKcERCCGEdCVClWKOBjWgxuoV5e+X4lVw5eqUIRyRsyd/fH6WUxcvNzY2aNWsycOBANm3aVGRjO378OEop6tatW2RjECWf1QFKKVVRKdVNKTVMKfWUUmqoabuCPQYoCiZ8ejgVGhj/aRLjElk4eCFJd5OKeFTClnr06EFERAQRERH06NEDgAULFhAaGsq0adOKeHT216FDB5RSbN26taiHImzMKS+VlFJ+wBPAQ0CjHOodBJYAX2qt/7LJCEWBuHi40G9uP2a2nkny3WQu7L3AulfWEf6f8KIemrCRV199ldDQUPN2YmIizz77LDNmzODVV19lwIAB+Pn5Fd0AhcinHM+glFIBSqn5wEngdaAxcB3YAfwM/Gh632kqbwy8AZxUSs1TStWx49hFHlVpVoXuH3Y3b/82/TdZO6oUc3Z2ZurUqZQtW5aEhATWrpVn4UTJlG2AUkp9ABwE+gN7gaeBhlprX611B611b631MNN7e621LxAIPAv8AQwADpn6EUWs1fhW3Nf7PvP2kpFLuHn+ZhGOSNiTu7s79evXB+DixYuZ9t+6dYu3336bpk2b4uHhgaenJy1atOD999/nzp3s1xQ7c+YMzz33HA0bNsTDwwMvLy8CAwMZP348hw4dytPYrl+/TmhoKEop+vXrR3y8ZWLjgwcP8uijj+Lv74+rqyvly5enW7du/Pzzzxb11q1bh1KKbdu2AdCxY0eL+3Fyya/ky+kS3wRgEfCm1vpAXjrTWh8BjgD/VUo1ASYBLwAvF3SgomCUUvT+ujefN/ucm+ducufqHRYNW8TwX4bj4ChzZUqjGzduAFC5cmWL8kuXLhEWFsbBgwfx8fEhPDyclJQUNmzYwMSJE1mwYAHr1q2jfPnyFu1WrVrFoEGDuHnzJtWrVyc83LhMfPLkST7//HOqVq1KYGAgOTlz5gy9evXi0KFDPPPMM0ybNg0Hh7T//+bMmcOoUaNITEykSZMmBAcHc+nSJTZv3sy6deuYPHkykyZNAqBatWpERESwatUqLl26RM+ePalUqZK5r4z/blECaa2zfAEtsttnzctW/RTXV1BQkC5JTm08pSeryXoyxmvTO5uKekjZOnToUM4V9k3Seg55e+18LHP7nY/lvf2+SZnbb3gw7+2PfWGLj8RCrVq1NKA3bNiQad+BAwe0o6OjdnZ21tHR0Rb7+vTpowEdGhqqr1+/bi6/evWqbtOmjQb0sGHDLNqcPHlSe3p6akBPmTJFJyUlWew/ffq03rNnj3n72LFjGtABAQHmsqioKF21alWtlNIffvhhpjFHRUVpZ2dnXbZsWb1mzRqLffv379fVq1fXSim9efNmi33t27fXgN6yZUs2n5QoLLn+zqYDROpcvl+z/dNZa/27jQKgTfoRtuEf4k+nNzqZtzdO2kj09ugiHJGwpZiYGFatWkXfvn1JSUlh+vTpFhMkTp48yZIlS3B0dOTLL7+kXLly5n0+Pj588cUXKKX48ccfOX/+vHnf1KlTuXXrFkOHDmXixIk4OjpaHLdWrVq0bNky23GtWrWKTp06ce3aNebNm8eECRMy1XnnnXdITExk6tSpdO/e3WJf48aN+fDDD9Fa88knn1j9uYiSSa7t3INC/hVCjfY1ANDJmoVDFhJ/XRY4LKk6d+5svu/i4+NDr169OHPmDKtWrWLs2LEWdTdv3ozWmvbt21OvXr1MfTVt2pSgoCCSk5PZsmWLuXz16tUAjBkzxurxzZw5k969e+Pi4sK6desYMGBApjrJycmsXbvWfF8qKyEhIQDs2LHD6jGIkilP08xF6eLg5EDfOX35ovkXxF+P58aZGyx/fDn95/VHKVXUw8u7ppONV361/tJ45Vfo8vy3taEePXpQpUoVtNZcuHCBzZs3Ex8fz4gRI9i2bZvFw7Lnzp0DoHbt2tn2FxAQQGRkpLkuwNmzZwFo0KCBVWM7c+YMjz32GA4ODqxZs4bg4OAs6126dIlbt24B4Ovrm2Ofly9ftmoMouSyKkAppcoD44DOQDXALZuqWmsdUMCxCTvyruVN7696M7/ffAAOLThEVLcogh4LKuKRCWtlfA7q77//pkePHuzfv5+hQ4eyc+dO8x8e2pSPMac/RFLr2EKVKlVo2LAhv/zyC88++yyrVq3Cy8srU73k5GQAnJycGDp0aI59Zry8KEqvPAcopVRdYBNQBcjtz2zJSloCNOzbkKAng9jz+R4AVj+7mprta1IxsGIRj0wURNWqVZk/fz5NmzZl165dzJkzh2HDhgGY70edPJn9EiynTp0CoHr16uaymjVrcuLECY4ePUqVKlXyPBZXV1eWL1/OgAEDWL58OV26dGHNmjX4+PhY1KtUqRKurq4kJCQwY8YM3N3d83wMUXpZcw9qKlAV2Ar0BZoAtbN5yQO6JUSPj3pQsZERkJLuJPHT4J9IvJNYxKMSBdWgQQPGjRsHwOTJk0lKMtJbderUyfzs0IkTJzK1O3DgAJGRkTg6OtKxY0dzeWoKpZkzZ1o9FldXVxYuXMjAgQOJjIykc+fOXLp0yaKOi4sLYWFhaK1ZuHChVf27uLgAmP+NovSwJkCFAqeBblrrJVrrg1rrM9m97DJaYXPO7s70n9sfJzfjZPrS/kusfVEyD5QGr7/+OmXLluXEiRN89913ANSpU4eHHnqI5ORknnjiCWJjY831Y2JieOKJJ9Ba88gjj1CtWjXzvgkTJuDh4cH333/PBx98YL4kl+rMmTNERUVlOxZnZ2d++OEHIiIi+OOPPwgJCbGYJQgwadIknJycePrpp1mwYEGmS43JycmsW7cuU2aM1DO9w4cPW/HpiBIht3noqS+MVEZz81r/XnkV6DmopLtaH5mudXJi/vuwkd0zdpufjZrMZH1oUd6fZ7AXa56puBfl9BxUqjfffFMDuk6dOjox0fj/7OLFizowMFAD2tfXV/ft21f36dNHe3t7a0C3aNFCX7t2LVNfy5cv1x4eHhrQfn5+ul+/frpv3766RYsW2sHBQb/99tvmulk9B6W11ikpKXrs2LHmfadPn7bY//3332s3NzcNaH9/f92zZ089YMAA3bZtW+3r66sB/frrr1u0WbhwoQa0q6ur7t27tx49erQePXq0PnbsmLUfqSggWz8HZU2A2ghszGv9e+WV7wCVnKj15n7GQ5xbBmidnJC/fmwkJSVFz+s3zxyg3iv3nr7y55UiHZMEqJzlJUDdvHlTV65cWQN65syZ5vLY2Fj95ptv6saNG2t3d3ft7u6umzVrpqdMmaJv376dbX8nTpzQ48aN0wEBAdrV1VV7eXnpwMBA/dRTT+nDhw+b62UXoFK98MILGtA1a9bMFEiOHTumn3rqKd2gQQPt7u6uy5Qpo+vUqaN79OihP/74Y33+/PlM/X3yySe6adOm2t3dXWPcA5cHd4uArQOUMurlTin1IEam8k5a6+22OHsrDYKDg3VkZKT1DU/Ogp2j0rb9HoL288DR1WZjs9admDt80fwLbpw1UuRUDKzI6J2jcS1bNGM6fPgwDRs2zL2iEKJYsOZ3Vim1R2ud9XMHJnm+B6W1XgE8D/yslHpbKdVBKeWvlKqZ1Suv/eYw+ClKKW16vZhDvSFKqS1KqRtKqVtKqUil1HilVG6Z2vPVzmZqR0D9Z9K2/1oKmx+GpOwTddqbe3l3Bi4ciKOrMY338qHLLB25lLz+ESOEELZk7Zfx78BF4DWMKecngFNZvLKfw5oHSqlWGAlmc/xmVEp9CswBgoEtwI6tlGEAACAASURBVC9AfeAT4CelVJYPTOS3nU0pBUH/gYbp8uj+vRo2PQhJt+1++OxUC67GP778h3n78KLDbH1PskILIQpfngOUUioUWIfxRa6Aa8DZbF75Tu6mlHIFZmEEwqU51OuH8dDwBaCp1vpBrXUfoB5wGOgDPGWrdnahFDR/HxpPSiu7+CtsCIfE2Ozb2VmzEc24/5n7zdu/vvErx1YeK7LxCCHuTdacQb0NuAAfAD5a64pa69rZvQowprcw1pV6EriRQ72JpvdXtNbmb0+t9UUgNQHZq1lcsstvO/tQykjX02xKWtnlrfBrN0iIKZQhZKX7h92pFVLL2NCwcMhCrh67WmTjEULce6z5Em4O7NFav6q1vm6PwSilWmOsQ/WD1jrbRGemJeiDgARgQcb9WutNwDmMrBdtCtquUDSaCC2npW1f3QXrwyD+SqEOI5WjsyMD5g/Aq4aRlubujbvMe3ged2/eLZLxCCHuPdYEqDuA3a7zKKXcgG8xLh0+m0v1Fqb3g1rr7GYV7M5QtyDtCkeD56DVZ+kKFDgUXd4xj0oeDFo0SCZNCCGKhDUBagvQyF4DAd4F7gOe1lrndtqQegkxp4wVZzPULUi7wlNvLLT+GrybQee14FI+9zZ2JJMmhBBFxZoA9U8gQCmV29mN1ZRS7YDngCVa63l5aOJpes9putst03tZG7SzoJR63DQtPdIuqf8DRkH4bnCrYPu+86EoJ03I2ZoQJYM9fletWW4jGPgG+Egp1R9YA/wFpGRVWWs9Oy+dKqXcTf3GYsyuy1Oz1MPksX5B21nQWn8JfAnGg7oF6StbDs6Zy87MB5+WULZu5n121v3D7lzcd5Ezm86YJ008tvsxfOvlvHZPQTg6OpKYmGhOBiqEKL4SExNtvhSKNQFqFsYXuwLaA+1yqZ+nAAVMwZi6/qjW+u88trlpevfMoU7qvpvpyvLbruidmQfbh4BbZQhbB+UCC/XwqZMmvgz+ktjoWPOkCXtmmihbtiyxsbFUqFA8ziSFENmLjY2lbNlsLzzlizUBajb2WeepD8ZZWIRSKiLDvtTlO8eaUi0d11qPwciqDlArh35rmN5PpyvLb7uidfca/PYY6BS48zf80hFCV0GF+3Nva0Opkya+7vA1yXeTzZMmBvw0wC4r8fr4+JhXcvXy8sLZ2blkrfgrRCmntSYxMZHY2FhiYmKoWbPASYQs5DkXn70opU6Tc8BIb5/WurlSqgbGZIYEwDurGXlKqWjAD+igtd5mKstXu5zkOxeftS5uMmWZMN0ic/KATkugSlf7HzuDfbP3sSRiiXk77N0wOr7WMYcW+Xf37l2uXbvGzZs3My3xIIQoeo6OjpQtWxYfHx9cXfN+NSUvufiKPEDlRCk1C4gAXtJaf5hh3x6gJRCR8X6XUioEI/v6BaC61jqloO2yU2gBCuDqbtjYE+6aHph1cIH2P0KNvoVz/HRWPbuKXR/vMjYUDFkxhHq96hX6OIQQJZNNk8UWQ++Z3v9tWo4eAKVUJSD1YaL3swgy+W1X9HxbQdct4G5aijslAbYOgBNfFfpQJNOEEMLesg1QpqnUBZqSoZRyVEo9XpA+sqO1/gmYgZH1Yb9SarlSahHGw8SBGEuDfGKrdsVGuYbQfRuUNZ2t6BT4bQwc+r9CHUZWmSZ+6PUDty8XXaJbIUTpktMZ1OfAIaVUhGkqeJ4ppdyVUiMxkq/OKMD4cqS1HgcMBaKAEKAHcBwj2Ws/rXWWNy3y267Y8KgF3bZC+XTJLva+DPv+WbjDME2aSF0u/trxa8ztPZfEO4mFOg4hROmU7T0opdQgjMSwfhgPry4A1gM7tNans6hfG2gLdAX6YUzXjsa4f5Qp711pUaj3oDJKuAGbe8OlzYCC9nOh1sBCH8bhxYeZ32++eY5ngz4NGLBgAA6OJfkKshDCngo8ScKUH+8FjAdoq5E2zfwuRs68WMAL8MXIdA7Gc1J/AZ8C07XW8QX4NxR7RRqgwFjgcNsgqP4g1LXL1dQ8+e3j31j97GrzdutnWxP+n/AiG48Qoniz2Sw+072ovsDDQCegehbVooENGPdwlhXLSQZ2UOQBCkBrY9mOIrZmwhp2frTTvN1jWg/aPFe4SeGFECVDXgJUnh7UNd2TWWB6oZSqAFQCygHXgUtaa5nCVVSyCk7xVyDqBQiaBq72S0eUXvf/686NMzc4vPAwAGteWINXDS8C+xVu1gshROlgTSYJM1O28aJZqEjkLvGm8bzUtUiI2WNkRS+T1UmvbSkHRZ/v+nDr71tEb48GDYuHLaZs1bLUaFcj9w6EECIduYtdGl1YD9f2GD/fOAS/tIcbRwrl0M7uzgxeNhifej4AJMUn8WPvH7n6p5xgCyGsIwGqNKrxMLSbA8p0gnz7DPzSDi4VzjpOZXzLMHTVUMpULAPAnat3mNNzDrcvyTNSQoi8kwBVWvk/AiHLwNEIEiTEwK9d4ezCQjm8T4APQ1YMwcndCJIxJ2P4sfePJMbJM1JCiLyRAFWaVesJXTeBWyVjO+WukRrpyPRCOXz1+6vTf25/lIMxiePcb+dYOGQhKcn3xARPIUQBSYAq7XyDofuOtNRIaIh6Dva8YKRJsrP7et9H+Mdpz0MdXXqU1c+tlpVyhRC5kgB1L/CsA922Q4W2aWVHp8HRjwvl8PePv5+2L6Yde/cnu9k5bWcOLYQQQgLUvcOtAoStB78+xnaFdlD3iUI7fLd/d6PRwEbm7bUT1nJwwcFCO74QouTJc4BSSo1QSuW2zDtKqTZKqREFG5awCyd36LAAmr5jTKBwsioHcIEoB8XD3z5MzQ5pK24uHr6Yk+tPFtoYhBAlizVnULOAMXmoNxr4Jl+jEfbn4AiNX8+cXUJriPvLrod2cnNi8NLB+N5nHDv5bjJze8/l7Lazdj2uEKJkssclvqJPCiesd+QjWNEQzq+x62HcfdwZtnoYXn7GOlKJcYnM6TmH85Hn7XpcIUTJY48Albo8hygpzi6A31+EpFuw6QE4Yd8TYG9/b0asH4FHZQ8AEm4m8F3377j4x0W7HlcIUbLkmIsvi3tJdXO4v+QENAS6ALttMDZRWLwCoUwNiIsGnQy/PWpkn2gyyW5Z0n3r+zJi3Qhmhc7iztU7xMfE81237xi5aSQVGlSwyzGFECVLbutBpZC2BpRK93O2TYAUYKDWepFNRljMFYvlNmwh7jxs7AXX96WV1RwEbb6x62SK83vOMztsNndj7wJQtlpZRm0ZRfk65e12TCFE0bPFgoWzSAtKERjLom/LpnoCcA5YqrXel02dUqfUBCiAxFjYMgAurE0r82kFIUvBvardDhu9PZrvun9H4m0jDZK3vzcjN4+kXI1ydjumEKJo2WzBQlNnKcAsrfWjthhcaVGqAhRASiLseQ6OfZZWVsYPOi0DnxZ2O+ypDaf4odcPJMUnAeBTz4dRm0fhWcXTbscUQhSdvAQoayZJ1AZeKtiQRLHn4AytPoXgT0A5GmVxf8EvHeDcCrsdtnbn2gxcNBAHZ+N/yWvHrvFdt++IuxJnt2MKIYq3PAcorfUZWTX3HlJ/PISuBOfUy2zarpf5AOr1rEf/ef1RjsbEjEsHLvF9j++Jvx5v1+MKIYqnPF/iMzdQyg0IBqoBbtnV01rPLtjQSoZSd4kvoxtHYHNvaDYFavYvlEPu/3E/i4YuMt/99Gvrx/C1w3HxdCmU4wsh7M+m96BMHT4P/Avwyq2u1toxzx2XYKU+QAEk3wVH18zlOgWUfdI5/v717ywbvcy87R/qz5CVQ3B2d7bL8YQQhSsvASrH56AydPYoMNW0eRg4AsTmf3iixMgqON04DNsGQ9vvoHxTmx+yxaMtSIxLZNXTqwA4vfE08/vOZ9CSQTi55vl/WyFECWbNb/ozGBddhmutf7DTeERJcPcqbHoQbp00lpJv9wP49bb5Ye5/6n4S4xJZ98o6AI6vPs5PA3+i//z+EqSEuAdYc32mPrBdgpMg9gjEXzZ+TroNmx+GQ/9nJJy1sfYvtydkUoh5++iyo8ztPVeWjhfiHmBNgIoDJO20gIrtjVV6PfxNBRr2vgw7RkCS7aeFh0wKof0r7c3bJ9ae4Pvw783ZJ4QQpZM1AWo70NheAxEljHcj6LELKnZIKzv9PfzSHm6dsumhlFJ0ea8LoW+FmsvObjnL7C6zibsqz0kJUVpZE6DeBBoopSLsNRhRwrhVhLB1EDA6rSxmL6wOsvmyHUopQv4ZQvep3c1l5yPP823ot9y6IMnzhSiNsp1mrpTqlEVxT+Bl4CfgZ4xLfilZtddab7bRGIu1e2KaeW60hhP/g8injFRJACjj2alGr9r8cHu+3MOKJ1eYn5PyqefDiHUjKFdTcvcJUVIUdJr5RrLOXq6A/qZXdnQufYvSRCmo+zh4N4Ut/eDOeUBDSoJdDhf0eBAuni4sHrEYnay5duwa33T8hhHrR+BT18cuxxRCFL6cgshmcl9eQ4g0FdpAeBRsGwjO3tD4DbsdqsmQJjiXceanQT+RnJDMjbM3+KbjNwxfN5xKjSrZ7bhCiMJjdaojYUku8WUhJdHIPuHsmbncwbaZIE6sPcHch+eSdMfIgu7u686wNcOoFlTNpscRQtiWrbOZC5E3Ds5ZB6dfu8O+1yEl2WaHCugewLA1w3Apa+Tpu3P1DrPDZnN2mzwRIURJJwFKFI7fX4ZLG+HgFNj0ANy9ZrOua3WsxYj1I3D3MVb+vRt7l++7f8/JdSdtdgwhROGzZsHCrGb1ZSUBuKK1Pp7vUZUgcokvD1ISYVNv+Ht1WplHbei0CMo3t9lhLu6/yHfdvuP2xdsAOLo4MuCnAdz3j/tsdgwhhG3YY0Vda25YxQLfAv/UWt+0ol2JIgEqj1KSYf8kOPhuWpmjGwR9DAFjjJmANnD1z6vM7jqb2Ggjj7GDkwO9v+5Ns+HNbNK/EMI2bH0PajOwA2OauQKuA38Ae4EYUxnAb8BJwBN4GtiilCpj3dBFqePgCM3egY6LwKmsUZYcD7seh+1DIdE2ifF96/syassoygeUByAlKYUlI5aw6a1NyIQgIUoWawJUuOn9ENBLa+2rtW6htQ7SWlfAeIj3IMZZVhOgHkZ6pCYYmdCFgBp9jBRJ5RqllZ35EVYFwbXfbXII71rejNoyikpN0qabb5y0kWWPLiM5wXYTNIQQ9mVNgHoDI9iEaa1XZ9yptV4DdMPI1/cvrfVpYAhwF+hX8KGKUqNcAyNIBYxJK7t1HNa2gYubbHKIslXLMmrLKOp0q2Mu2ztrL3N6zZEl5IUoIawJUIOADVrrS9lV0FpfBDYAA03b0UAUxlIdQqRxKgOt/wft5oCTaUp6uUCo0Npmh3Ar58aQn4fQ/NG0iRin1p/i6w5fc/3MdZsdRwhhH9YEKD+Ms6Hc3AWqp9uOBrJYklUIwH+IkX2iUidoP8+YOGFDjs6O9J7Zm87vdDaXXT54ma/afMX5PedteiwhhG1ZE6CuAJ2UUu7ZVTDt6wRcTVdcHmNChRBZ86oHXTeBV4YTba3hr2UFXghRKUWn1zvR5/s+OLo4AnDrwi1mdZrFnyv+LFDfQgj7sSZALQcqA/OVUjUy7jSVzQMqAcvS7WqAMatPCOscmwGbHzKeobp7Nff6uWg6tCnD1g7DrbxxlpYYl8jch+ay+7PdBe5bCGF71gSoScAZ4AHguFJqo1LqW6XULKXUBuA48CDGEhyTAJRSQUBNYJ1thy1KvesHIep54+fzK2BVc7i8rcDd+of4M3r7aLz9vQHQKZqV41ey9sW16BSZhi5EcZLnAKW1vgy0wzg7csK4lDccGAGEmMpWAB1MddFa7wGctdb/svG4RWlXth7UfyptO+4vWBcCB94tcC6/Cg0qMHrnaKrfn3ardMfUHSwYuIDEO4k5tBRCFKZ8ZTNXStXECFCpv+HngS2mqeX3FMkkYWd/LYedIyEhXe6+Cu2g3XfgWSfbZnmRGJfIoqGLOLLkiLnMr40fg5cOxqOSR4H6FkLkzKapjkTWJEAVgtvRsP0Ry0t8Tp4Q9B+o82iB0iSlJKew9sW1/Paf38xl5euUZ9CSQVRuUrkgoxZC5ECW2xClg0cN6LIRmr4DyrTGZtIt+G0MbH64QBMoHBwdCJ8WTvj0cHOyrpiTMXzV5isOzD1Q8LELIfIt2wCllKppejlm2M7Ty5pBKKWclVJdlFJTlVI7lVJ/K6USlFLnlFI/KaVCc2k/RCm1RSl1Qyl1SykVqZQar5TKMQDnt50oAg5O0Ph16LETvBqklV//wyaLILZ+pjWDFg/C2cPoKzEukYWPLGTNhDWkJKUUuH8hhPWyvcRnyl6eAgRqrf+0Mpu51lrntJx8xmN1BX4xbV4A9gC3gUCM1EkAb2c12UIp9SkwDogH1gOJQBegLLAYGKC1znRXPb/tMpJLfEUgKQ72vgrHPjXOrCp1tFnXlw5eYl6feVw7lnbPyz/Un/7z+st9KSFsqED3oJRSpzECUpjW+lS67TzRWte2YqBhGMFiutZ6S4Z9g4A5gKNpLBvS7esH/IQR1DpprY+ZyitjpFxqCDyntZ6eoc98tcuKBKgidPM4lK2b9/I8ir8Rz5IRSzi67Ki5zMvPi4ELB1rM/BNC5F+pmSShlJoJjAa+1lqPTlceCQQBEVrr2RnahAAbMYJQda11SkHbZUUCVDETvQi2DoDAidD4X+Dokq9udIpmy5QtbPjXBvOfZY4ujvT6tBctx7S04YCFuDeVpkkSqesw+KUWKKX8MIJMArAgYwOt9SbgHFAFaFPQdqIEuHPBWF9KpxgLI65tCzcO56sr5aDo9EYnhvw8BDdvI/NEckIyyx9bzvLHl5N0N8mWIxdCZKGkBKh6pve/05W1ML0f1Frfyabd7gx1C9JOFHc6Gbybpm3HRMHqlnD4Q0jJX0Cp17Mej0U+RuWmaVPOo/4XxaxOs4j9yzaLLAohsmZ1gFJK1VVK/Z9SaqtS6qhS6oN0+9oopR5XSnnbaoBKqSrASNPmwnS7Uu9xncmh+dkMdQvSThR3ZapD2DpoMRUcTJf2kuPh95dgbTu4vj9f3foE+PDo9kdpMqSJuezcrnN80fILTm88bYOBCyGyYlWAUkqNBg4AEzDSHtUFKqSrUhGYAfSxxeCUUk7A90A5YL3Wenm63aZFhLidQxe3TO9lbdAu/bgeN01Jj7x8+XIO3YhCpxyg4QsQvgfKp60DxbXdsKol/DEJkvOyaowlFw8X+nzfhx7/6YFyNB6Yirscx+yus9kxbYcsJy+EHeQ5QCml2gNfYEzLfglojfnRRrPVQCzQ20bj+xxj6nc0MCzjkEzv1n4z5Ledmdb6S611sNY6uGLFivntRtiTd2Nj1d5m76adTekkOPCWcdnvxiGru1RK0ebZNkT8GmGecq6TNWtfWMvCwQtlpV4hbMyaM6iXMb7Ue2qtp2qtM61RoLVOBI5iTNMuEKXUdIyZexeALlrrCxmq3DS9e5K91H0305Xlt50oaRycodFr0HOfkb8vVfxFcM3/Hxa1OtXi8ajH8WtjnrPDwfkH+bz555zdejaHlkIIa1gToNoCu7TWO3KpFw1Uzf+QQCk1FXgGuIwRnI5lUe206b1WDl2lrlt1Ol1ZftuJkqpcA+i2BYL+C04eEPQxuBXszNeruhcRGyMIHps2S/bGmRvMCpnFxskbJfuEEDZgTYAqB/yVh3ouGEtv5Itp0sULGKvydtNaZ3ctJnXqeaMcVvltlaFuQdqJkkw5wH1PwT+OQ61HMu8//QMkxFjVpZOrEw989gADfhpgXgRRp2g2vbmJWSGziDllXX9CCEvWBKhL5G1W230YzxFZTSn1Psb9rRiM4LQvu7pa62ggCiMgDsiirxCM56YuADsK2k6UEu5VMmc/v7QVtg+DFYEQvdjqLgP7BfLkvifxD/U3l0Vvj+aL5l+w/4f8zRwUQlgXoLYBLZVS2T75q5TqBtTHyMRgFaXU28ArwHWM4JSXs5f3TO//VkqZc9sopSoBn5k2388iG0R+24nSJiURdo0BNMRfgC19YUt/Y4FEK5SrUY7h64YTNiUMByfj1+pu7F0WDV3E4hGLuRtr/cxBIe51eU51pJRqDWzHODsag7GMexIwS2v9qFKqE0bOvMpAkNY6z386KqV6A0tNm5HAwWyqHtFav5+h7WfAWIzZhetIS/rqBSwB+meTLDZf7TKSVEelQPQi2D3eCFCpnDyg8SS471mr0yWd23WOhUMWEnMi7RJf+Trl6ftDX/xa++XQUoh7h81z8SmlJgD/hzGbLxbjy/wGxpd7BYwp3C9orf9j5UBHAt/koeomrXVoFu2HAOOBJhhJZY8AXwMzcjoLym+79CRAlRIJMRD1Ipz82rLcqyEEfwJVwqzq7u7Nu6x6ehX7vk27Sq0cFaFvhtLh1Q44OJaUJC5C2IddksUqpcKBN4FgLJ+D2g/8U2u9zNqBlmQSoEqZi5sgcjzcyHASX2swtPjQyFZhhQNzD7DiiRUWl/hqdapFn+/7UK5GOVuMWIgSya7ZzJVSvhiTJhyBaK31+Xx1VMJJgCqFUhLh6H9h/yRj5d5UNfpCx4XZt8vG9dPXWTRsEdHbos1lbt5u9PykJ02GNEEVYMl6IUqqUrPcRnEmAaoUiztv5PE784ORjaLXAfCql3u7LKQkpbD53c1sfmszOiXtd67eA/V4YMYDcjYl7jkSoAqBBKh7wMWNEHsY6o21LE+Oh7vXoEy1PHd1dttZFg9bzPXT181lLmVd6PZBN4IeD0I5yNmUuDcUdEXdEQU5eMaFAEsrCVD3sP1vw+EPoMmbcN/TRmqlPLh78y7rX1vP7k93W2SErBVSi3/87x/41vO104CFKD4KGqBSKFhCVcf8ti1JJEDdo26dgp8DjbMoAK8G0PwDqP5g5geBs3F261mWjVnG1aNXzWVObk6EvhVK2+fbmp+nEqI0KmiA2kj2ASoEuIgxLTtLWuvOeRtmySYB6h51LcrIPhGbYcXeSqHQ8kPwCcpTN0nxSWx6axPbPtiGTk77dasaVJXeX/WmSrMqNhy0EMWH3e5Bmc6uZmmtH83v4EoLCVD3sOQEODodDrwNSRkS3/sPM5b68KiZp67+/v1vlo1exoXf0x4WdnByoP2r7en0RiecXPOd3lKIYikvAUquIQiRX44uEPgS9D4O9caBSndV+/T3sLw+7J0ICTdy7apqi6qM+W0MXd7rgqOr0U9KUgpb3tnCFy2+IHpHdC49CFH6SIASoqDcKkGrT41p6H4PpZWn3IVD78PZeXnqxtHZkQ6vduDJfU9Ss0PamdeVw1f4uv3XrHp2FfE3ZFFEce+QACWErZRrAJ2WQJeN4GO6cuHVAOpYdyW8wn0VGLlpJD0/6YmLZ+pqwLDr4118Uv8Tfv/md4tnqYQorSRACWFrlUOgx2/Qbg4E/xccMtw/itkHVzMtSG1BOSjuH38/Yw+MJaBHgLn89qXbLHt0GV+1/Yq/frMu47oQJY0EKCHsQTmA/xCo0tWyXGvYPQ7W3G8s63H9QI7deNfyZuiqofT7sR9efl7m8nO7zvFVm69YOmopty7cyqEHIUouCVBCFKboRXBlu+nnhbCyKWx7BG5k+8QGSikaD27M+CPj6fhGR/MkCoC9s/by3/r/ZfvU7SQn5Lo6jBAlSk7PQXXKod1GYDXwfnYVtNabCzSyEkKmmQurxB6Dva/AXxlW7lUOUGsINP5Xrvn+Yk7GsHbCWo4ssQxqvvf5Ev6fcOqG182mpRDFR1FmktBa63viwQ0JUCJfru2BPybD+RWW5coRao+Axm+AZ50cuzix9gSrn13NlSNXLMrv630f3T/qjk+Aj40HLYTtFDRAnaZgqY5q57dtSSIBShTIlV3Gsh5/r7YsV05w/xcQkPMMwOTEZHZ9sotNkzdZrDnl6OJI2wlt6fhax7SZgEIUI5LNvBBIgBI2cXm7EagurDO2lSM8cDjPy3vcuniL9a+tZ+83ey3+rPSo7EGnNzoR9HgQji73RHpMUUJIgCoEEqCETV3aDH/8CzxrQ5tvLPfFX4Hk2+BRK9vm53afY9XTqzj32zmLcm9/b0Imh9B0WFNZbl4UCxKgCoEEKGFzWhtZKBzdLMv3ToTD/we1HoHAl8G7SdbNUzT7vtvHhjc2EPtXrMW+Cg0rEPZOGA36NJCVfEWRkgBVCCRAiUKRGAtLakJiurx+VXtC4CtQqVOWS3wkxSexe8Zutk7ZStyVOIt91YKrETYljDpd60igEkVCAlQhkAAlCsXNE/DbGLi0MfM+39bGGVX1h8Ah832mu7F32fmfnWz/cDsJNxMs9vmH+hM2JYwabWvYaeBCZE0CVCGQACUK1ZVdxiq+0YvINMm2bH1o+KIxTd3RNVPTuCtxbP33VnZ/spuk+CSLffX/UZ+wd8Oo3KSyHQcvRBoJUIVAApQoErF/wpGpcHIWpFieFVGpE3TdlH3Tc7FsfnszUTOjLBZJREGTR5rQ8Y2OVGxY0T7jFsJEAlQhkAAlitSdC8aiicdmpN2fuv9LqPtYrk2vHb/Gxkkb2f/j/kwnYw0ebkD7V9vj19rPDoMWQgJUoZAAJYqFxFg4/iWcmQvdtlrOAExNUFv9Qaganuk+1cU/LvLr67/y54o/M3XrH+pP+1faE9AjQCZTCJuSAFUIJECJYkXrzDP6Lm6E9Z2Nnz1qQ72xRoYKV1+LatE7otn63lb+XJ45UFVuVpkOr3YgsH8gDk7yHJUoOAlQhUAClCj2tgyA6J8syxzdoNZgqDcefC2/Iy4duMS2D7ax/4f9lveogPJ1ytP2xbY0H9kcZ3dne49clGISoAqBBChR7N06Ccc+hxNfQcK1zPt9W0P98VBzgMWlwetnrrPjox1E/S+KpDuWs/48KnnQ+rnWtBrbCjdvt4w9CpErCVCFQAKUKDGS7hj3qP78BGKiMu93rQBtv4Nq4RbFcVfi2PXJDq9ydAAAGWNJREFULnb9dxd3rt2x2OdS1oWgx4NoNb4V5WuXt+foRSkjAaoQSIASJY7WcPU3+PNTODs/3TR1Bb1Pgqd/ls0SbiUQNTOKHVN3ZEqhhIJ6PesRPC6YuuF1Jd+fyJUEqEIgAUqUaPGX4MRM4xJg+eYQssxy/40jcOh9qDPSlFLJgeSEZPb/uJ9t/97GlcNXMnXp7e9N0JNBtHi0BR4VPQrn3yFKHAlQhUAClCgVUpKM+1NulSzL974Kh/5t/OzhD7UjoE4EeNZGp2iOrTzG7k93c3z18UxdOro40mhgI4LHBePXxk+mqQsLEqAKgQQoUWqlJMPSmnDnfOZ9lUKgziio0Q+cPbl24hqRn0ey9+u9me5TAVRpXoXgccE0GdIEFw9ZQFFIgCoUEqBEqaW1MZnixDdw5gdIiMlcx8nDCFI1B0GVriQmKA4tOMTuz3ZnWpMKwNXLlWYjmxH0WBCVGlfK3J+4Z0iAKgQSoMQ9IfkunFtu5P77exXolMx1HjgI5QLNm+f3nCdyRiT7f9ifaZo6GGdVTYc3pfEjjSlbtawdBy+KIwlQhUAClLjn3PkbTn0PJ7+B2MNGWblG8MCBDPUuQMzv3HFtx77vDrP7s91cO5b5OSzloKjTtQ5NhzelwcMNcPGUS4D3AglQhUAClLhnaQ3XIo2p6h61of44y/2HP4LfJ4BLefDrg67Rn1OH6hA1cz9Hlh4h+W5ypi6dPZxp2KchTYc3pXaX2jJdvRSTAFUIJEAJkY01bYznrdJz8QG/h0nweZBDWyuz77tjnN54OsvmnlU8aTykMU2HNaVK8yoyC7CUkQBVCCRACZEFnQK/vwxnF0Dc2azrOLpB5TDuuIexb1N9or6N5vKhy1lWrdioIoH9A2nwcAMqN6sswaoUkABVCCRACZEDreHqLuMy4NkFEBeddb1u29EV2nBh7wX++O4P9v+wn9sXb2dZtVytcjR4uAENHm5AzQ41Jbt6CSUBqhBIgBIij1JTLEUvgvM/w41DRrlrBehzwWKdqpRbF7m9ejRHdgX8f3tnHh53Ve7xzzeTNG2apk3bJF1SoC1lFWRfZVEKCoqCgNtFQL3KQ/W640Wv1+264C7KrkJB4eoFtKJel0exIAqVUi4UEGhpC21om5IuadM0aZr3/nHO0Ol0ZjJpJplJ5/08z3lO57znd5a3k98755z3nMODt9XS3prZcWLU+FEccO4BHHTeQcw8ayZVNX7C+nDBDdQQ4AbKcfaQLcuh5bfQ2wUHf2JX2bLb4eFLATBEpw7nxedm8dhvx7Hs/5ro2b67IaocVcnMs2Zy0HkHccCbDqBmYs1Q9MLZQ9xADQFuoBxnEHjw7WFaMANGFes3HcizC6bw7IIptCxtZkdP5S55VCGmnTSNGWfNYMbsGUw9dqpPBZYYbqCGADdQjjMIbHwSWu4NI6y2hzNvDI6sbDufX980m3VPZXawgHAtyPTXTmf67OnMmD2DiQdNdEeLIuMGaghwA+U4g0xXG6y9D9b+JYT2Z3aVn3IPTHsrbUvaePZXz/LMvGeY1Xwr3Z3VrFrazEvLptC9rXqXR8ZMGcOM2TOCwTpjBmOm+EkWQ40bqCHADZTjDDGdq2Ht/GCsWufDWQ9B9YSd8h1d2F3jUO82AKxXtLY00LJ0KquWNtPyfDPrVjVgtnPKr+GQBqbPns6+p+5L8wnN1E2tG9o+lSFuoIYAN1COU2K0PgB/Oi1nlq7OEby0PKxf3f+L03ZzuqibVkfzCc2vhMlHTaZyZGWW0pw9IR8D5Rp3HGfvou5gOP4WePmh4Na+6cnd1rCqR3Uz/ZAVTJm5jvm/PAvYKa+bsJEpTU+z+oEmnr67HqyCiqoKJh85maknTKX5hGamnTiNsfuO9XWsQcZHUAPER1COU+Js3xLODGxbAC8vCE4XnauDbPIb2H78vbz4txdZft9yWh5uoaF6HudcMg+A7m1VtK5qpHVlE2tfjGFlI9s6ahjdNJrm45tpOqKJpsOaaDyskfH7j/fzA/PEp/jyQNK7gCuAw4EE8AxwK3CDWQ7XoYgbKMcZZpjB1lXBYFXVweSzdhX/40No6XU5i9jUVkfrykaeePDVPPnQYa+kV46qpOGQhmCwDm98xXDVNtUOSleGMz7F1weSrgPmANuAPwPbgTOAa4EzJF1kZrsfuew4zvBFgtHTQsgkrj8UJp0JG5+AbWsz5hk7oZ2xE9p5admUXdJ7Ons4aMZPmVDXxvqHx7N43gTa1oxnW08ztfvNpPHwJhoOaWD8zPHUz6ynrrnOR1w5KFsDJekCgnFaA5xqZktiehPwF+B84EPANUVrpOM4Q8+sK0IA6FwLmxbDxsXBYG14AjY9FU6/APZ785mcuN9xtC5uZe3itWxZvYWZhz3P1Jkv7VZsV+cI1q8Zz4an6ln7wFiWrK/j2ccOhdrp1M+op35mfTBc8d/10+vL/m6ssp3ik7QQOBq41MxuT5OdBswnGK+puab6fIrPccqM3h7YvCQYrIknwuh9XhFtXddB9Z+mkLD2vIq67SuXsOLpGbukXXzV7XRvG8Hm9XVs2z6R3hFTqKidTGJMI5X1kxjZ2MzoSROpnVRL7aRaahpqhuUozKf4siCpmWCcuoG70uVmdr+kFmAqcALw96FtoeM4JUtFJYw9OIQ0aibWwOzfBwO2eQlsWYq1L8Han6Nix+bd8rev33W/VaJqOzMPW9ZnE7avrmTrczXc8NnL6eyoZXTjaGon1TKuOcHBR/4D1UwkMWosFaPqSIyuo3L0OCrHjKWqrp7q+glUjx1D9diRJEYk+qyrmJSlgQKOjPFTZtaZJc8jBAN1JG6gHMfJBwkaTgwhmQTIDLrWBaPV8WK4dmTrKi7/55fYsLKLDc9vYMOyDXSuejqvaqqqexhb3U5XZzXWa2xZs4Uta7bQs66Vw9/5w90f2BpDXFLr3SHa19dx3VVXUl1XTXVdNVWjq5g8bTnHvuZejARGAmJsVIZ/KwFUQkWCrV2TWb7m7Rxz+TE0vqpxYHrLQrkaqOkxfiFHnuQta9Nz5HEcx+kbCUY2htBw8ivJI4CmcdB0WFNI6DkCNhwMW1fR27GS7tZl9KxfAdvWoe3rSbCBEYl2Kip66O6qZsTYMXS27fyNXVOb7ff2rlQkjIpELz3beujZ1kNHa7h7a1zlC0yZ9kwfTwdWLZ3KI9fOYNbZs9xAFZikz2fmG9ECW2K82yFdkj4AfABgn332SRc7juPsGZWjoOEkACqAkYdkyGMGPVsY0b2RT71nGju6d9DR2sHm1ZvpbnmcdZvasa71aMdWKqwD0UlCnVQmOkkkuqiq6iJRuWO38wkBKhJ97qx5hd7esO5VUTV461/laqCS27/3yEPEzG4GbobgJFGoRjmO4/SJBFVjQgASIxLUNddR11wHx04Fzum7jB3dTNjRyWfm1LBt0za62rvYvnU7dLbS2nk+vT3dWE8P1rMd2xFDz3astwfi5+6pY3nDNacx8cCJg9bVcjVQydXKXLvnkrLdVzYdx3GGM4kRKDGCqhFQVVPFmMnJiaLJwKuL2bJdGH6+iYVhRYz3zZEnuYtvRY48juM4ziBRrgbqsRgfKmlUljzHpuV1HMdxhpCyNFBmthJYRHCiuShdHjfqNhM26j40tK1zHMdxoEwNVORrMf66pP2TiZIagevjx6vzOTDWcRzHKTzl6iSBmd0t6QbCSeaLJf2JnYfF1gHzCIfGOo7jOEWgbA0UgJnNkfQg8EHgNHZet3ELeV634TiO4wwOZW2gAMzsTuDOYrfDcRzH2ZWyPc28UEhaR+4jk3IxEXi5gM0pB1xn/cP11T9cX/1jIPra18wacmVwA1VEJC3s67h5Z1dcZ/3D9dU/XF/9Y7D1Vc5efI7jOE4J4wbKcRzHKUncQBWXm4vdgGGI66x/uL76h+urfwyqvnwNynEcxylJfATlOI7jlCRuoBzHcZySxA1UgZD0Lkl/lbRJ0hZJCyV9UFLeOpZUJekMSd+W9LCk1ZK6JbVIulvS6YPYhSGnEDrLUfZXJVkMnyxEe4tNofUlaZSkT0l6RNJGSVslLZd0l6ST+y6htCmkviQ1S/qBpGcldUraJmmJpBslzRiM9g8Vkg6U9BFJP5X0jKTe+Hdz4QDLHbj+zczDAANwHeF23k7gN8AvgfaY9gsgkWc5s+MzBqyOZf0cWJyS/qVi97eUdJal7GOBHqA3lvfJYve31PQFTAeWxOfXAr8C/gf4B9ANfLbYfS4VfQFHAhvisysJ53TOA1bFtM3AScXu8wB09b2U90tquLDY+i+6coZ7AC5IMSizUtKbgKej7CN5lvU64G7glAyyt8eXrgGvLXa/S0VnGcquBp4CWuIfxbA3UIXWFzAaWJr8wQNUpcknAAcUu98lpK+/x2duTtUVUAX8OMoeL3a/B6CvfwW+AbwNmAnMH4iBKug7sdjKGe4BWBgVfkkG2Wkp/1EVBajrR7G8Hxe736WqM+Dr8flzgbl7iYEqqL4IV80YcFux+1bq+gJGsnNEMSmDfEqKvKbYfS+Q/gZqoAqmf1+DGgCSmoGjCVMid6XLzex+wi/5ScAJBagyebtvcwHKKgqDqTNJxwOfAO40s18PvLXFp9D6kjQCeH/8eHXhWloaDML3awdh5gJAGeTJfTodhOmssqbQ+ncDNTCOjPFTZpbty/lIWt6BMCvGqwtQVrEYFJ1JGgncBqwHPrLnzSs5Cq2vowlTeCvN7J+STooOJTdJ+qKkEwfa4CJTUH2Z2Xbgz/HjFyVVJWXx31+OH39scYhQ5hRU/2V/3cYAmR7jXKeZv5iWd4+QNAm4LH68ZyBlFZnB0tlXgAOBd5jZ3nQadaH1dViMl0iaC1yaJv+cpHuAd+d4wZQyg/H9mgP8njDyPFvSwph+LFAPXANc2c927q0UVP9uoAZGbYw7cuTZEuMxe1qJpErgp8BY4M/DfPqq4DqTdBLwUWCemf18AG0rRQqtr/ExPpVwQee3gBuBtph2PWGRux14b38bWwIU/PtlZsvid+x24Gx2nWJfCDwQR1pOgfXvU3wDIzknPdhD+xsJV9GvBC4e5LoGm4LqTNIo4FbCC3VOIcosMQr9HUv+zVcSpqWuNLPnzWyjmd0LnBfrunSY7u8p+N9kNE5PAvsDbyHcgdRA0FU9cI+kzxWqvmFOQfXvBmpgbI5xbY48SdnmHHmyIuka4H3AGuAMM1uzJ+WUEIXW2VeBA4CPm9lwXpvLRqH1lZrnh+lCM1sIPEp4N5yeR3mlRkH1JWkcYc/TGOANZnavmbWZ2ctm9ivgDQTniP+UNCtXWWVCQfXvBmpgrIjxvjnyTEvLmzeSvg18GFhHME5L+ltGCbIixoXS2fmEDbmXSpqfGggvD4ArYtqP9qC9xWZFjAulr9Q8y7PkSaZPyqO8UmNFjAulrzcSRksPm9mydKGZLQUWEEakp+fbyL2YFTEuiP59DWpgJN2+D5U0Ksui8rFpefNC0jeAjxPWBs40s6f3vJklxWDorIKwvyIbM2IYl2d5pUSh9bUo5d8TCD9+0pkY4y0ZZKVOofW1T4w35cizMcbjc+QpFwqqfx9BDQAzW0n4gx8BXJQul3QaYUF1DfBQvuVKuprgFbSBYJweL0iDS4BC68zM9jMzZQoEt3OAK2PaEYXrydAwCPpqIfzih7CumV5ePXBU/LgwXV7qDMLf5EsxPjrVxTylvCqC6z5kH5GWDQXXf7F3LQ/3AFzIzp3R+6ekNxKO3NntWA/CTv5ngK9lKO+/4jMbgKOL3b/hoLMc9cxl7zhJotDfsXPZeQbfESnpI4GfRdlC4n1xwy0UUl/xmY74zLVAdYqsGrghytYDY4vd9wLpbz59nCTRx/er3/rPWk+xlbE3BIJrrhEWS39NOAxxU0z7JWkHI6a8OOempb+ZncemPBLzZQpXFbvPpaKzPurYKwzUYOgL+GaUdwEPxDJaYtoqUs5QG46hkPoi7BVLnoPZAtwby3wppm0Dzit2nwegq6OAh1NC8lDX51LT+/n96pf+swVfgyoAZjZH0oPABwlrIQnCr4tbgBvMrDfPolLnsI+JIRP3M8yPqSmgzsqCQuvLzK6U9Hfg3wg7+msIGyi/A1xtZpnWpoYNhdSXmd0maTFhr90pwFlR1EI4LPY7NrzXiOuA4zOk77FXYqH071e+O47jOCWJO0k4juM4JYkbKMdxHKckcQPlOI7jlCRuoBzHcZySxA2U4ziOU5K4gXIcx3FKEjdQjuM4TkniBspxsiBphSSL4Wt95L0jJe/8IWrioCIpIWmxpBckVafJjpP0V0mdklolXS9pdI5yHpX0fLy/K1OeKZK2Srp7MPriDE/cQDlOflwiKZFJIKmOcO3H3sYVwKuAL5hZVzJR0lTgPsLpA38inLl2BXBXlnI+TDhO5wrLco28mb1EvM1X0umF6oAzvHED5Th9sxCYApyZRf4OYBTh/MS9Akm1wBcJJ3Tfnib+FDAaeI+ZnUswPn8GzpZ0bFo504AvAXeY2R/7qPZqwtmA3xp4D5y9ATdQjtM3c2N8WRb5ZcAO4CdD0Jah4lLC2ZBzzWxHmuwowgGp/w0Q5bdE2YlpeX8AdBPuNsuJmb1MOFj0aEkn73nTnb0FN1CO0zcLgKeBt8QrwF9B0oGEl/IfCFNdGZE0W9J1kh6X1CapK67t3Cbp4CzPjJR0laRFkrbEZ1ZLekjSlyWNTMt/nKS7JLVI2i5pk6Slku6U9Lp+9nlOjNNHTxAuOtyUduDn+hi/0iZJbwXeQriPqzXPepN3eM3JmcspC9xAOU5+zCW8fN+Zln5ZjG/t4/kbgfcRrm34K/C/hJHFJcBCSa9JzSypAvgt4d6dGYQT7O8hGMppwH+QckOwpDOBBwl38bQSrjS4j3Cv2IXA2/LsJ5JmAYcAS81sRYYsK4AGSamn7x8U4+WxjDHA9wlXefSlm1T+QhiNvinbmp9TRhT7LhIPHko1EF7ERrj2ZBLBuCxIkScIVy60EW4QTV7UNj9DWecB49LSBFwen3malAsCgVNj+qPA6AzPnQzUpKTdF/O/M0PdE+jH5ZfA+2NZt2eRXxHlP4llHwmsJNwj1BDzfJ+wnnTQHuj9sVj+ccX+DngobvARlOPkgZmtAX4PHJcyJXcWwXniTjPr7uP5eWa2MS3NzOwm4O/AwYRRS5KmGP/VzDoyPPc3M9uaIf/vMtTdZmaP5u7hLhwR439mkf+QcIndxcDLhCu+m4FPmdk6SccQ7gH6qpk9k3xIUnWeo6Lk3UpH9qPNzl6IGyjHyZ+5Mb4sLZ5LHkhqlnS5pO9K+rGkuZLmEkZnAAekZF9EmOp6n6Q5kprSy0vjHzG+U9LJA5wea4xxWyahmfUQLqF7P3AT8G3gJDO7MdZ7M7CEeKmmpNdLepzgWNEl6Q+S9s9Rf3I9q68+O3s5fqOu4+TPvYSX9rslfZPgALA4n9GJpC8CnyH331xd8h9m9rykjxFcrq8DrpO0jDDa+hXwS9vVu+7ThJHP2TF0SHqUMPX3EzNbln83GRvj9mwZ4ojxRzGk8tHYjteaWVd0O/8NYVT0NqCesK52n6RDzWxzhuKT9Y7LIHPKCB9BOU6exJfyncBkwsJ/NXk4AEi6APgc0EkYdcwkrB/JzER01yasLaXW9wNgX8Kazx2ENa+LCRtiF8YNwsm8a4CjgTMII5dFhI20XwCelfTefnQ1ORVZlzNXGpL2IeydusXM7o/JnyAY5QvM7C4zuxn4d4Kjx7uyFJWsd0N/6nf2PtxAOU7/mBvjNxGcJu7I45mLYvwZM/uRmS2zXU9UyDrdZWZrzOxGM7vYzPYjjE4Wx/iqtLy9ZnafmX3azE4lODBcRTAQ16UatD5IuoRPyDN/kmuBDsJG3iSvBl42s6UpaQ+lyDKRrDdf13RnL8UNlOP0AzNbRHDnbgPusvz29yTdsVemC6LDRd7OAGb2OHBN/JjtBZ/M22FmXwdWEVzkD8yzmkUxPiRnrhQkXQicC3zMzNaniDoIp2ykkjyzz7IUl6x3URa5Uya4gXKcfmJmp5jZRDPLNkWVTtKT7f2SRiQTJTUSNqbuti4l6XWSzpFUmZaeAM6JH19ISf9kPFYovZxjCFOSvQRDlQ9/iXH6qRAZiSOza4A/mtmdaeIngNGSUvdhXRbjxzKUVUs4/68dN1BljztJOM7g8z3Chtw3AkslLSCMKk4jjKrmEfZJpXI48F1gk6RFhFMqagjrSpOBNcDXU/J/FvimpH8S3MO7COs8JxF+iF5tZllPukjFzJZLegI4XNJ0M1vexyNfJTg/XJFBdjXwL8Adki6J+U4ClrFz7S2V1xLW2n5jux+x5JQZPoJynEEmetAdBfyM4AhxLmHf082EUcqmDI/9muBwsIiwRnUBcArBMH0eONzMXkjJ/0HCaKyX8JI/H5gay3m9mX26n82+PsaX5Mok6TiCYfpSJk9BM3uOsF9sAcGB41UEJ4/T0/d3RS5Nq98pY2SWbRrYcZxyJd7t9AJhqm3WUIxmJE0kTEM+aWbHDHZ9TunjIyjHcXYjjm4+D0ynj1FUAbmK4Lr/ySGqzylxfATlOE5GokPGY4SNuwdYyqWFg1DXFGAp8Dszu2Cw6nGGF26gHMdxnJLEp/gcx3GcksQNlOM4jlOSuIFyHMdxShI3UI7jOE5J4gbKcRzHKUncQDmO4zglyf8DZ/k3HrT/Om0AAAAASUVORK5CYII=\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "plt.plot(num_sol_imp[:,2]/m0, num_sol_imp[:,0], label = \"Simple Rocket\", color = \"purple\")\n", | |
| "plt.plot(num_sol_imp2[:,2]/m0, num_sol_imp2[:,0], \"--\", label = \"Rocket\", color = \"orange\")\n", | |
| "plt.xlabel(\"Mass (%)\")\n", | |
| "plt.ylabel(\"Height (m)\")\n", | |
| "plt.legend(loc=\"best\")" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 11, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "The max height of the simple rocket is [598.47543949] m\n", | |
| "The max height of the rocket is [425.43300118] m\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "print(\"The max height of the simple rocket is\", num_sol_imp[:,0][-1:], \"m\")\n", | |
| "print(\"The max height of the rocket is\", num_sol_imp2[:,0][-1:], \"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": 12, | |
| "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", | |
| " N = 1000\n", | |
| " time = (.25-.05)/dmdt\n", | |
| " t = np.linspace(0, time, N)\n", | |
| " dt = t[1]- t[0]\n", | |
| " \n", | |
| " num_sol = np.zeros([N,3])\n", | |
| " \n", | |
| " num_sol[0,0] = 0 #y0\n", | |
| " num_sol[0,1] = 0 #v0\n", | |
| " num_sol[0,2] = .25 #m0\n", | |
| "\n", | |
| " for i in range(N-1):\n", | |
| " num_sol[i+1] = heun_step(num_sol[i], lambda state: rocket(state, dmdt=dmdt, u=250), dt)\n", | |
| " \n", | |
| " error = 300 - num_sol[:,0][-1] \n", | |
| " return error" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 13, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "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 = np.zeros(ns)\n", | |
| " for i in range(ns):\n", | |
| " f[i] = func(x[i])\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": 14, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "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]" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 15, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "number of brackets: 1\n", | |
| "\n" | |
| ] | |
| }, | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "array([[0.08571429],\n", | |
| " [0.07857143]])" | |
| ] | |
| }, | |
| "execution_count": 15, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| } | |
| ], | |
| "source": [ | |
| "incsearch(f_m, 0.05, .4)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 16, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "0.07912598052406028\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "dmdt, out = mod_secant(f_m, 0.00001, .08)\n", | |
| "print(dmdt)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 17, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "y0 = 0\n", | |
| "v0 = 0\n", | |
| "m0 = .25\n", | |
| "\n", | |
| "dmdt = 0.07912598052406028\n", | |
| "tf = (.25-.05)/dmdt\n", | |
| "t = np.linspace(0, tf, 1000)\n", | |
| "dt = t[1] -t[0]\n", | |
| "N = int(tf/dt)\n", | |
| "\n", | |
| "\n", | |
| "num_sol_imp3 = np.zeros([N,3])\n", | |
| "\n", | |
| "num_sol_imp3[0,0] = y0\n", | |
| "num_sol_imp3[0,1] = v0\n", | |
| "num_sol_imp3[0,2] = m0\n", | |
| "\n", | |
| "for i in range(N-1):\n", | |
| " num_sol_imp3[i+1] = heun_step(num_sol_imp3[i], lambda state: rocket(state, dmdt=dmdt, u= 250), dt)\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 18, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "Text(0, 0.5, 'height (m)')" | |
| ] | |
| }, | |
| "execution_count": 18, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaYAAAEaCAYAAABaefMNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd5xV1dX/8c+aygAz9CaDVJViodoVFU1iNLFiixFbbFge/T1G04uJJVGjsSUmiiZoEoVonqjRWBC7VAEREKQjvQ5l+vr9cc4UxpnL3Jk7c+/c+32/XvM6nLP3PneN4Kw55+yztrk7IiIiiSIt3gGIiIhUp8QkIiIJRYlJREQSihKTiIgkFCUmERFJKBnxDqCl69y5s/fp0yfeYYiItCgzZ87c5O5damtTYmqkPn36MGPGjHiHISLSopjZirradCtPREQSSsIkJjO7wcyeM7MFZrbZzErMbKOZvWFmF5uZRRh7kZm9a2bbzWynmc0ws/FmFvH7a+g4ERFpOol0K+82oCvwKfABsAvoDZwEjAHONbOz3b28+iAzewS4DigE3gRKwv4PA2PMbKy7l9X8sIaOExGRppVIiekCYLa776p+0MyGECSOM4BxwIRqbecQJJd1wPHuvjg83g2YApwFXA88WOOcDRonIiJNL2FuWbn7ezWTUnh8PvBIuHtKjeYfhNvbKpJLOGY9cG24e3stt+YaOk5ERJpYS/nBWxpuCysOmFk+MAIoBp6vOcDdpwJrgO7AkY0dJyIi1RSsgwmnQsH6mJ864ROTmfUFrgl3/12taVi4ne/ue+oYPr1G38aMExGRClN/Ays/gqn3xPzUifSMCQAzuwwYDWQC+cDRBAn0Lnd/oVrXvuG2zrnwwMoafRszTkREftUVSouq9mc8EXxlZMOPN8TkIxLxiukYgkkOFwHHh8d+AvyyRr+24fYrz6Wq2Rluc2MwrpKZXRVOLZ+xcePGCKcREUkuc85+h3+VHc0ezwLAM3LgkLFw07yYfUbCJSZ3v9LdDWgNDAEeAH4OfGRm+1XrWvFeU7QrHTZ0XPUYH3f3ke4+skuXWitqiIgknV1Fpdz48loKPIdsSigmC8qKIDsPcrvF7HMSLjFVcPc97v6Zu99KMIvuMIJ3jCoUhNu2XxlcpaKtoNqxho4TEUlpd7z0GSs276az7eA5TmHbd/6DjbgMdsZ2AkTCPWOqwwTgXuBbZpbp7iXA8rCtd4RxvcLt8mrHGjpORCRlvf7Zev4+fRUA15TczO/OP4yuB+TDASNj/lkJe8VUwzaCKeMZQMfw2OxwO8TMcuoYN6pG38aMExFJSRsLirh98tzK/dMO7cGZQ3s22ee1lMR0PEFS2gZsAnD3VcAsIAsYW3OAmY0mmNW3Dviw4nhDx4mIpCJ357bJc9m8qxiA7nmt+PWZBxOhfGmjJURiMrPjzOw7ZpZdS9sxwBPh7hM16tfdFW7vMbMB1cZ0BR4Nd++uWV+vEeNERFLKs9NW8tbCqmng9449jPats5r0MxPlGVN/gudID5vZLIKrldzw+OCwz8sE08YrufskM3uMoIzQPDN7g6pirHnAi+w9YaJR40REUsnSjTv51UsLKvcvO6YPxx7Quck/N1ES01TgDuA44ECCl2qNIEFNBia6+4u1DXT368zsPWA8wYu56cBC4Engsbqueho6TkQkFZSUlXPzc3PYUxLcpDqga1tu+8bAZvnshEhM7r4M+Gkjxj8LPNtc40REkt0jU5YwZ9U2ADLTjd+dP5RWmenN8tkJ8YxJREQSx+yVW3norSWV+zefciAH92zXbJ+vxCQiIpV2F5dyy3NzKCsPiuMc3qcjVx/fv1ljUGISEZFKd7y0gGWbglKibbMzuO+8w0hPa7qp4bVRYhIREQBe/XQtf5u2snL/598eQq+OrZs9DiUmERHhy217uG1yVYXw0w7pwTnDm666QyRKTCIiKa6s3Ln5H5+wfU8JAD3b53DnWYc0aXWHSJSYRERS3KNTlvDxsi0ApBk8cMFQ2rXOjFs8SkwiIils5oqtPPDm4sr9G046gFF9OkYY0fSUmEREUtSOwhJu+vvsyqnhI3t34IaTBuxjVNNTYhIRSUHuzo9e+JTVW/cAkNsqgwcuGEpGevzTQvwjEBGRZjd51hr+PefLyv27zz6U/A7NPzW8NkpMIiIpZtmmXfz0X59W7p8/shenHdojjhHtTYlJRCSFFJeWc+PfZrO7OKga3q9LG3727cH7GNW8lJhERFLIvf9dxLw12wHISk/j9xcMo3VWQiw0UUmJSUQkRby1cD2Pv7O0cv/73zioWauG15cSk4hICvhy2x5ueW5O5f6JB3Xh8mP6xjGiuikxiYgkuZKycm7422y27Q5KDnXPa8V95w0lrZmrhteXEpOISJK777+fM3PFVgDS04yHLhpGxzZZcY6qbkpMIiJJbMrCDfxh6heV+//vawfGveTQvigxiYgkqbXb93DLc59U7o8+sAvXNPNqtA2hxCQikoRKy4L3lbaGz5W65WVz/3mHJexzpeqUmEREktD9r3/O9OXBc6U0g4cuHE6nttlxjqp+lJhERJLM24s28Ojb1Z8rHcThfRP7uVJ1SkwiIklk3fbCvd5XOu6Azlw7OvGfK1WnxCQikiRKysq5/tlZbNlVDEDX3Gx+d37ivq9UFyUmEZEkcdcrC5mxouq50u8vHEbnFvJcqTolJhGRJPDS3C958v1llfu3fn0gR/brFMeIGk6JSUSkhVuyoYDbJs2t3D9lcDeuGd0vjhE1jhKTiEgLtquolGsmzmJXuL5Sn06tue+8wzBrWc+VqlNiEhFpodyd2ybPZcmGnQC0ykzjsYtHkNcqM86RNY4Sk4hICzXh/eW8NHdt5f6dZx3CoB55cYwoNpSYRERaoBnLt3DnKwsq979zxP6cPTw/jhHFjhKTiEgLs7GgiPHPzqK03AE4LL8dP/3W4DhHFTtKTCIiLUhFcdb1O4oA6NA6k0cvHkF2RnqcI4sdJSYRkRbk7v8s5MOlmwEwgwcuGEbP9jlxjiq2Murb0czygJOBk4BhQDegPbAV2ADMAqYAb7j7jtiHKiKS2l6YvZo/v1f1Eu3/jDmQ0Qd2iWNETWOficnMDgZuAC4CWgM1J8d3BPoDRwHXAbvN7BngYXf/NLbhioikpk/XbOf2yfMq908Z3I0bThoQx4iaTp2Jycy6AXcC4whu+a0HXgI+BD4DtgA7gDygEzCYIDmdAFwFXGlmTwE/cvf1TfYdiIgkuU07i7jqLzMoKi0HoH+XNi1m0b+GiHTFtBhoA7wAPAm86u7lEfq/DjxoZunAqcDlwGXAuQS3/EREJEolZeWMf2YWX24vBCA3O4M/XTKS3Bb+Em0kkSY/TAEOc/dz3f2VfSSlSu5e5u4vufvZwFDg7RjEKSKSkn798gI+XrYFCCY7PHjhUPp1aRvnqJpWnVdM7n5GY0/u7vOAMxt7HhGRVPT8jFU89cHyyv3/d8qBnDSwW/wCaiaaLi4ikoDmrNrGj16smj926sHdGX9ick52qEmJSUQkwWwsKOLqv86kOJzscFC3XO4d27Irhkej3u8xVTCzVsBIYD+gVV393P0vjYhLRCQlFZWWce3EmazbEUx2yGuVweOXjKBNdtQ/rlusqL5TM7sduB3IrUd3JSYRkSi4Oz964dO9lkd/6KLh9O7UJs6RNa9oKj/cTPBeE8BcgunkO5siKBGRVPT4O0uZNHN15f7tpw5MysoO+xLNFdO1QAlwprv/p4niERFJSa9/tp67X11YuT92RD7fO67lLo/eGNFMftgfeEdJSUQkthas3cH//H02Hqxiwag+HfjVWQenzGSHmqJJTGsJyhCJiEiMbNpZxJVPz2BXcRkA+R1y+EOSLWMRrWgS07+AY80sK9ZBmFmmmY0xs/vM7CMzW2tmxWa2xswmmdkJ+xh/kZm9a2bbzWynmc0ws/FmFvH7a+g4EZFYKCot4+q/zmTNtj0AtMlK54lxo+jUNjvOkcVXND+Afw7sAv5iZh1jHMdo4A3gFqA3MJOgRt8W4Bxgipn9sraBZvYI8AzBFPZ3CWr2HQg8DEwKa/fFbJyISCy4Oz/45zxmhjPwzOChi4ZxUPf6THpObvWe/ODu28zsCGAqsNTMZgCrgdpq6Lm7XxFFHOXAZOBBd3+3eoOZnU+QQH5iZlPcfUq1tnMIltpYBxzv7ovD490Iav2dBVwPPFjjnA0aJyISK398Zyn/nLWmcv+Hpw5KiXJD9WFe8bRtXx3NcoBJwDf46ppMNbm7x+yKw8z+DFwBPFk94YXJcQQwruYLvWY2mqCA7DqgZ/UitA0dV5uRI0f6jBkzGv7NiUjKeW3+Oq6ZOLNyssN5I/O555xDU2qyg5nNdPeRtbVFM138VwTLWWwGJgJLaL73mGaH2/yKA2aWT5BcioHnaw5w96lmtgboCRwJfNCYcSIisTBn1TZuqjYD7/A+HfnVmYekVFLal2gS03kEy6gPdfc1++ocYweE27XVjg0Lt/PdfU8d46YTJJhhVCWYho4TEWmUVVt2c8XT0yksCW7E9O7UmscuHk5WhuZbVRfNf41OBO8xNWtSMrPuwKXh7uRqTX3D7YoIw1fW6NuYcSIiDbZ9dwmXPTWdTTuLAWjfOpMJl2oGXm2iSUxLgWadqWZmGQS3DdsBb7r7v6s1V6yUtSvCKSpuNVaf5tLQcdXjuiqcWj5j48aNEU4jIgLFpeVcPXEGSzYEP1qy0tN4/Lsjk37Bv4aKJjE9CZwQzlxrLn8AxgCrgItrtFXckK3f7I3Gj6vk7o+7+0h3H9mlS+rVsRKR+nN3bp88l4+WVtUn+O3YQzm8b6zfukke0SSm3wEvE7xTNKapX0I1swcJZuKtA8a4+7oaXQrCbaRfOSraCqoda+g4EZGoPfDGYv45u+oJyK1fP4gzhvaMY0SJL5rJD1+E297Af4ESM1tH3e8x9W9oUGZ2H3AjsJEgKS2updvyavHUpVeNvo0ZJyISlUkzV/Pgm1U/vi4Y1YvrTmjwj8aUEU1i6lNjP4ugsGttGnybzMx+Q1ABYjNwirt/VkfXiinkQ8wsp44ZdqNq9G3MOBGRent/ySZunzy3cv+4Azpzx5mpW5g1GtEkpiafoWZmdwO3EkxLP8Xd59TV191XmdksYDgwlhoLE4YvyuYT3Ar8sLHjRETqa+G6HVwzcSal5cHv6AO75/Lod4aTma5p4fURTUmiSNOrG83M7gBuA7YRJKX6XK3cRfCS7D1m9oG7LwnP1RV4NOxzdy3VGxo6TkQkotVbdzPuyWkUFJYC0DU3mycvHUVuq8w4R9ZyJMQi8mb2beDH4e4S4IY6LncXuvvdFTvuPsnMHiNYxHCemb1BsJjhGCAPeJGgKOteGjpORCSSrbuKueTJaazfUQRA2+wMJlw2iv3a58Q5spalzsRkZpnuXtLYD6jnearPmxwZftVmKnB39QPufp2ZvQeMJ6hSng4sJJje/lhdVz0NHSciUpvdxaVc/vR0lm4MXpHMSk/j8UtGMGS/dnGOrOWps4irmS0Dfgb81etb6XXv8QaMA37m7klbQUFFXEWkpKycq/86k7cWbgDCJSwuHMbph+4X58gSV6QirpGexO0EJgCLzexHZlbXDLyaH9bbzH4CLCa4+tgRbcAiIi2Fu/PDf86rTEoAPzt9sJJSI0R6xnQYwZpFPwfuAH5pZksIZqotIJjOvYPgmUwnYDBwFNCfoLrCZoI1jf7QRLGLiMTdvf9dxPMzV1fuX3dCfy49JmlvEjWLOhNT+IzlYTObQFBE9TpgEEGl79pu7VXMVphLMLNtorvvjmm0IiIJ5Kn3l/HIlC8q98eOyOfWrx8Ux4iSwz5n5bn7LuAR4BEzGwCcAAwFuhIUV90GbABmAVPcfXlTBSsikihemvslv3ip6v3/kwZ25a6zta5SLEQ1XTx832dJE8UiItIiTFm0gZv/8UnlYn/D9m/PIxcNJ0Mv0MaE/iuKiERh2rItXDtxJiVlQVbq36UNT44bRU5Ws64KlNSUmERE6mne6u1c/lTVCrQ92+cw8coj6NAmK86RJRclJhGReliyoYBxE6axsygoNdQlN5tnrjyCHu1U1SHWlJhERPZh1ZbdfOfPH7NlV7AserucTP56xeH06dwmzpElJyUmEZEINuwo5OInPq6sf9c6K52nLhvFwO55cY4seSkxiYjUYdvuYr77xDRWbA5eyczKSOPP40YybP8OcY4suSkxiYjUYkdhCeOenMai9QUApKcZj1w0nKP7d45zZMlPiUlEpIadRaVcNmE6c1ZvB4KirPefdxinDO4W58hSQ70Tk5ktNbN76tHvLjP7Yl/9REQS0e7iUi6fMJ2ZK7ZWHrvjjIM5Y2jPOEaVWqK5YuoDdKlHv85hXxGRFmVPcRlXPDWDacu3VB77+bcGc/GRveMYVeppilt5OUBpE5xXRKTJFJaUcdVfZ/Dh0s2Vx3582iBVCo+DmCYmM2sHHAOsi+V5RUSaUlFpGddOnMm7izdVHvv+Nw7iyuP6xTGq1BWxiKuZLa1x6FwzOyHCubqF2ycaH5qISNMrLi1n/DOzmbJoY+Wxm08+kOtOGBDHqFLbvqqL96n2Zwfahl91KQZeBG5rXFgiIk2vpKycm/4+mzcWrK88dsNJA7jp5APiGJXsKzFV3Fw1YCkwCbi1jr7FwEZ31/MlEUl4xaXl3Pi32bw6v+rJw9Wj+3HLKQfGMSqBfSQmd19R8Wczexp4t/oxEZGWqKi0jPHP7H2ldMWxfbn9GwO10F8CqPdCge5+WVMGIiLSHApLgokO1Z8pXXlsX3502iAlpQQR1Qq2IiItWWFJGd/7y4y9Zt9dM7o/t33jICWlBBJVYjKzPGA8MAbYD2hVR1d39/6NjE1EJGb2FJdx5V+m8/6SqveUbjhpALeccqCSUoKpd2Iys17Au0AvgskQkXhjghIRiaVdRaVc/tR0Pl5WVdHh5pMP1Oy7BBXNFdOdwP7ALOAeYCGwoymCEhGJlaAg6zSmL6+qfXfr1w9i/Il6TylRRZOYvkZQ0eFEdy9oonhERGJm2+5ixk2YzpxV2yqP/eDUgVw9Wk8aElk0iSkPeEVJSURagg07CvnuE1XrKUFQ+05lhhJfNIlpOZDZRHGIiMTMqi27ufiJjytXnjWDX55xMN9VlfAWIZoirhOB0WbWqamCERFprCUbChj7hw8rk1J6mvG784YqKbUg0SSme4BpwCtmNriJ4hERabBP12znvD9+xLodhQBkZaTxx4tHcOYwLfLXktR5K8/M3qrlcCYwCphrZiuBlUB5Lf3c3cfEJkQRkX2btmwLVzw1nYKioFxnm6x0/jRuJEf37xznyCRakZ4xnRChLY2g8nifOtr1HpOINJspizZw7cSZFJYEvye3y8nkqctGMWz/DnGOTBoiUmI6sdmiEBFpoH/OWs33J82ltDz4fbhLbjYTrziCg7rnxjkyaag6E5O7T23OQEREouHu/PGdpdz9n4WVx/I75PDMlUfQu1ObOEYmjaUiriLS4pSXO3e8/BkT3l9eeWxg91yevvxwuuXVVcJTWgolJhFpUYpKy7jluTm8PHdt5bEj+nbk8UtG0i5Hr1omg2iKuNY2S682xcAmYCbwN3dft4/+IiL1sqOwhKv/MpMPl1ZVCP/mId25/7yhtMpMj2NkEkvRXDGdEG6duquLV2+7EPi1mV3v7k82LDwRkcCGHYWMmzCdBWurakePO6o3P/3WENLTtGxFMokmMZ0IfBu4GfgY+BuwguA9pj4EiehI4AGCq6WTgEuBP5rZAnf/MGZRi0hKWby+gEsnTGfNtj2Vx77/jYO4dnR/raWUhKJJTGXADcCN7v5wLe0Pmdl1wIPASe5+hZm9BzwB3AQoMYlI1N5fsolrJs6koDB4cTY9zbj77EMYO7JXnCOTpmLu9XsX1sxeA7q6+7B99JsNbHD3r4f7XwBZ7p6U/4pGjhzpM2bMiHcYIknpuemr+OEL8yrfUWqdlc4jFw3nxIFd4xyZNJaZzXT3kbW1RVMrbxQwvx795od9K3wGdInic0QkxZWXO795dSHfn1z14my3vGyeu/ooJaUUEM2tvCyCFWz3ZX/2Xh5jD1AUTVAikroKS8r4f8/N4eV5VdPBB/XI48lLR9KjXU4cI5PmEs0V01zgaDM7pa4OZnYycEzYt0IvYGPDwhORVLJpZxEX/umjvZLSSQO78vw1RykppZBorpjuA54H/m1mT1M1K8+B3gSz8saFfe8HMLN2wDBgcqwCFpHktHh9AZc/PZ1VW6pm3o07qjc/OX0wGenR/A4tLV29E5O7TzazHwO/BK4Mv6ozgqnjP3P3ikTUFfgt8EoMYhWRJPXGZ+v5n398ws5wyQoz+Onpg7nsmL5xjkziIaqSRO5+p5m9ClwPHA9UrL71JfAO8Ii7z6jWfzHwkxjFKiJJxt159O0vuPe/i6iYIJyTmc7vLxzGKYO7xTc4iZuoa+W5+yzg8iaIRURSyJ7iMr4/eS7/nvNl5bGe7XP40yUjGbxfXhwjk3hLmBu3ZnaQmd1kZhPNbKGZlZuZm9m59Rh7kZm9a2bbzWynmc0ws/FmFvH7a+g4EWmcL7ftYewfP9grKR3etyP/d/0xSkqSUNXFryWoEBEVM3sEuA4oBN4ESoAxwMPAGDMb6+5lsRonIo0zc8UWrv7rLDbtrHqL5DtH7M/PvjWErAz9TigREpOZPUkw4+6H7r4+3K8vd/croozlU4KJEjMIau09AYyONMDMziFILuuA48NnWphZN2AKcBbB87AHYzFORBrnH9NX8pMX51NcFiyBnpFm/OzbQ/jukb3jHJkkkjpLEplZOUFiGuTun4f79eXu3qga9Gb2NkFiGuvuk+roMwMYAYxz97/UaBsNvE2QfHq6e3ljx9VGJYlE9q2wpIyf/998/j59VeWxDq0zeeziERzZr1McI5N4iVSSKNKtvMvC7doa+wnBzPIJkksxwftVe3H3qWa2hmDm4JHAB40ZJyINs3rrbq6dOIt5a7ZXHhvYPZc/XTKSXh1bxzEySVR1JiZ3fzrSfgKoKCY739331NFnOkGCGUZVgmnoOBGJ0tTPN3LT32ezbXdJ5bEzhu7HXWcfQuusRHrELYmkJf/LqHjzbkWEPitr9G3MOBGpp/Jy5+EpS/jdG59Xvp+UkWb85PTBXHJUb62hJBE1KDGFpYZGEVQNX+Hu8biqaBtud0XoszPc5sZgXCUzuwq4CmD//etT11YkdWzfXcLNz33CWws3VB7rlpfNo98ZzojeHeMYmbQUUc3NNLN24ey8DcBrwESqlSYys+vM7EszOzK2YdYeTrit34JSjR9Xyd0fd/eR7j6ySxet6CFSYe7qbXzr4ff2SkpH9uvISzccp6Qk9VbvxGRmbQhmq10KbAX+Q9UP+QqvAt2BM2MTXkQF4bZthD4VbQXVjjV0nIjUwd3587tLOeexD1i5ZXfl8auP78fEK46gS252HKOTliaaW3n/CxxGcJV0jbvvrjmF3N2XmtnnwEkxjLEuy8NtpBcgKlbNXV7tWEPHiUgttu4q5tZJc3hjQdVVUm52Br8591BOPaRHHCOTliqaxDSWoFjr99w90sJ/K4EhjYqqfmaH2yFmllPHDLtRNfo2ZpyI1DB9+RZu/Nts1m4vrDx2aH47Hr5wOPt30lRwaZhonjH1A6bvIykBbAKa/I05d18FzCJYWXdszfbwRdl8ghdlP2zsOBGpUl7uPDJlCRc8/tFeSemKY/sy6ZqjlZSkUaJJTCVAq3r0y6dqVltTuyvc3mNmAyoOmllX4NFw9+5aqjc0dJxIytuwo5BxE6bx29cWUVYezCFq3zqTP18ykp+cPlj17qTRormVtwgYZmat3L2wtg5m1oHgOdSsaAMxs+FUJQWAweH2TjP734qD7n5ktT9PMrPHCArAzjOzN6gqxpoHvEhQlHUvDR0nkupe/XQdP/jnXLZWe2F2ZO8O/P7CYezXXkufS2xEk5gmAXeHX/9TR587CWa0PdeAWPKAI2o5fkCkQe5+nZm9B4wnqK2XDiwEngQeq+uqp6HjRFLRrqJSfvnvz/jHjKpad2Zw3Qn9ufnkA7X0ucRUnUVcv9LRrDVBqZ6BBM9e/gncSzCF/HmC5zWjgXnA4e5e3ATxJhwVcZVkN3PFVm557hNWbK6aBt49rxX3n3cYRw/oHMfIpCVraBHXvYTTw79GkISOBo4Km0aHX0awXMWZqZKURJJZSVk5D721hIffWkx5td9fTz+0B78+8xDatc6MX3CS1KIqSeTua4CjzewbwDcJZuqlA6sIXrh90et7CSYiCWvZpl3c/I9P+GTVtspjudkZ3HHmwZwxdD/VupMm1aBaee7+KkGVBxFJImXlzoT3l/Hb1xZRVFr1mPXwvh25/7zDyO+gaeDS9FpydXERiaFlm3Zx6/NzmLFia+WxzHTjllMO4qrj+5GepqskaR4NrS6eTvASbZ3vNbn7yrraRCRx1HWVNLB7LveOPYyDe7aLY3SSiqJKTGZ2BPBL4DggUlVGj/bcItL8lm7cyfcnzd3rKikjzbjuxAFcf+IAvSwrcVHv5GFmxwBvUJWQtgI7miIoEWlapWXlPPXB8q9cJQ3qkcdvzz1UV0kSV9Fc1fyCICn9Cfixu29smpBEpCl9umY7t/9zLp+uqfq9MiPNGH/iAMbrKkkSQDSJ6XBggbtf3VTBiEjT2VVUyv2vf86E95ft9V7SoB553Dv2UIbsp6skSQzRJCYD5jZVICLSdN5csJ6f/ms+a7ZVrfKSlZHGTWMO4HvH9dNVkiSUaBLTPILVaUWkhdiwo5Cf/3s+r8xbt9fxYwZ04tdnHkKfzm3iFJlI3aJJTA8Cz5jZUHf/pKkCEpHGKy0r568freD+/35OQVFp5fGObbL48WmDOGtYT1VvkIQVTa28f5jZYOB1M/sp8LLeVRJJPB8t3czP/jWfResL9jp+7oh8fvjNQXRskxWnyETqp87EZGZlEcY9DDwc4Tcud3e9xyTSjNZtL+TOVxbwf3O+3Ot4v85t+NWZB6sSuLQYkZJHY67zdY9ApJkUl5bzxHvLeOitxewurvp9snVWOjecdACXH9uH7FE5n+IAABFKSURBVIz0OEYoEp06E5O7a5qOSAJzd95etJE7XvqMpZt27dX27cP244ffHET3dnVWDRNJWLrdJtICffblDn79yme8v2TzXscP6pbLz789hKP6d4pTZCKNp8Qk0oKs31HIva8tYtKs1VRf+Sw3O4NbvnYg3z2yt5Y5lxZPiUmkBdhdXMofpy7l8XeWsqek6jlSeppx4eG9+J+TD6Rz20h1lUVaDiUmkQRWWlbOpJmruf/1z9lQULRX20kDu/KDUwdyQLfcOEUn0jSUmEQSUHm58/K8tfzu9c+/MrFhUI88fnzaII7R9G9JUkpMIgnE3ZmyaAO/fe1zFqzde1WZrrnZ/O/XD+Kc4flaTVaSmhKTSIL4aOlmfvvaImZWW7QPILdVBteM7s9lx/ShdZb+l5Xkp3/lInE2Y/kWHnxzMe8u3rTX8ZzMdC47pg9XH9+fdq0z4xSdSPNTYhKJk4+Wbub3by7mgy/2fhcpM9246PD9GX/SALrm6gVZST1KTCLNyN15f0mQkKYt37JXW5rB2cPzuWnMAfTq2DpOEYrEnxKTSDNwd97+fCO/f3Mxs1du26stPc04Y+h+jD9xAP27tI1ThCKJQ4lJpAmVlJXz0twvefydZV+ZZZeRZpwzPJ/rTuxP705asE+kghKTSBPYWVTK36et5Mn3lvHl9sK92rLS0zhvVD7XjO5PfgfdshOpSYlJJIbW7yhkwvvLeebjFRQUlu7VlpOZzvmjenHN6P6q+i0SgRKTSAzMXb2Npz5Yzr/nfElJme/V1rltFuOO6sPFR/amg1aPFdknJSaRBiouLeeVeWt56oPlfLJq21fa+3Vuw5XH9ePs4T1plamF+kTqS4lJJErrdxTyzMcrefbjlWzaWfSV9hG9O3DV8f04ZVA30lQ6SCRqSkwi9VBe7ny4dDPPTlvJa5+uo7R879t1WelpnH5oDy45ug9De7WPU5QiyUGJSSSC9TsKmTRzNf+YvoqVW3Z/pb1bXjYXH9GbCw7fny65Wg9JJBaUmERqKC0r5+1FG/n79FVMWbSBshpXRwCH9+nIuKP78LUh3cjUirEiMaXEJBJauG4HL8xaw4ufrGH9jq8+O8prlcHZw/M5f1QvBvXIi0OEIqlBiUlS2rrthfzrkzW8MHsNC9cV1NrnyH4dufDw/fn6kO6aXSfSDJSYJOUUFJbw6qfrePGTNXzwxWb8q3fq6Nw2m3NHBFdHfTurXJBIc1JikpSwo7CENxes55V565j6+UaKS8u/0ic7I42vDenO2cN6cuwBnfXsSCROlJgkaW3fXcLrC9bzyry1vLd4E8VlX01GZnB0/06cNSyfrw/pRm4rLcgnEm9KTJJU1m7fw1sLN/Df+et5f8mmr7xvVGFQjzzOHLofZwztqbp1IglGiUlatPJyZ87qbby1cANvLtjAZzWWlqju4J55fPOQHpx6cA89NxJJYEpM0uJs313CB19s4q2FG5iyaAObdhbX2fewXu355sHdOfXgHuzfSUtMiLQESkyS8IpKy5i1YhvvLdnIe0s2M2/1Nuq4Q0dmunFkv06cNLArpwzupvWORFogJSZJOGXlzoK1O/ho6WbeXbyJacu2sKekrM7+ndpkceLArowZ2JXjDuxC22z9sxZpyfR/sMRdYUkZc1ZtY/ryLUxbvpVZK7ays6i0zv5pBofkt+e4AZ0ZM6grh+W3VxVvkSSixCTNyt1Zt6OQOau280mYjOat3l7rVO7qendqzbEDOnPsgM4c3b8z7VprWrdIslJikia1fXcJc1ZvY+7qbXyyajtzV29jQ8FX69DV1C0vm1F9OnJMmIx6ddSzIpFUkfKJycwuAq4FDgXSgYXABOAxd4/8a7xUKi93Vm7ZzYK1O1iwroCFa3ewcF1BrUtF1KZflzaM6t2RUX07cnifjvTqmIOZbs+JpKKUTkxm9ghwHVAIvAmUAGOAh4ExZjbW3et+6p6CysqdNVv38MWmnSzduIslGwpYsLaAResKIk5QqK51VjoH92zH0F7tGb5/e0b26UjntlrLSEQCKZuYzOwcgqS0Djje3ReHx7sBU4CzgOuBB+MWZJyUlpWzvqCI1Vt2s2LLbpZt2sXSjUEiWrF59z6fB1WXmW4M6pHHofntODS/PUN7tad/l7aka7KCiNQhZRMT8INwe1tFUgJw9/Vmdi3wNnC7mT2UTLf0ysudLbuL2VhQxIaCItZvL2T1tj2s2bqH1Vt3s2bbHtZuL6x1cbx96dw2m0E9chnYPZeB3fMY2COXAV3bkp2hpSJEpP5SMjGZWT4wAigGnq/Z7u5TzWwN0BM4EvigeSOsn7JyZ2dhKdv3lFR+bdtTvNf+1l1VSWhjQRGbdxU3KOlU1yU3m76d29C/Sxv6d2nLwO55HNQ9V0uLi0hMpGRiAoaF2/nuvqeOPtMJEtMwYpyYdhSWcN9riyhzp6y84gvKyssp83AbHi8sKWdPSRl7issoLCkL/hzuF9WydEOsdG6bTX6HHPI75NCvcxv6dWlL385t6NulDXmqwC0iTShVE1PfcLsiQp+VNfrGTFFJOU9/GOmjm1b71pl0aZtNl9xsuuZm07NDDvkdWtOzfQ49O+TQs32OVmoVkbhJ1cTUNtzuitBnZ7jNrdlgZlcBVwHsv//+UX94LB/8t83OoF1O5le/WgfbiiTUNa8VXXKz6dw2S898RCShpWpiqsgMDXrY4u6PA48DjBw5MupztM5K52ffGkx6mgVfZlV/rnYsLc3IyUwnJyudnMx0WoV/bh1uszPS9K6PiCSdVE1MBeG2bYQ+FW0FEfo0SKvMdC47JuZ3CEVEkkJavAOIk+XhtneEPr1q9BURkWaQqolpdrgdYmY5dfQZVaOviIg0g5RMTO6+CpgFZAFja7ab2Wggn6AqxIfNG52ISGpLycQUuivc3mNmAyoOmllX4NFw9+5kqvogItISpOrkB9x9kpk9RlBZfJ6ZvUFVEdc84EWCYq4iItKMUjYxAbj7dWb2HjAeGE3VshdPomUvRETiIqUTE4C7Pws8G+84REQkYO6NK+iZ6sxsI5FLG0XSGdgUw3Ak8ejvOLnp77fhert7l9oalJjiyMxmuPvIeMchTUd/x8lNf79NI5Vn5YmISAJSYhIRkYSixBRfj8c7AGly+jtObvr7bQJ6xiQiIglFV0wiIpJQlJhERCShKDHFgZldZGbvmtl2M9tpZjPMbLyZ6e+jBTOzg8zsJjObaGYLzazczNzMzo13bNJ4ZpZpZmPM7D4z+8jM1ppZsZmtMbNJZnZCvGNMFnrG1MzM7BHgOqAQeJOq+ny5wAvAWHcvi1+E0lBm9gBwUy1NY919UnPHI7FlZicDr4e764CZwC5gMHBwePwOd/9pHMJLKvoNvRmZ2TkESWkdcKi7n+7uZwEHAAuAs4Dr4xiiNM6nwG+B84EBwNT4hiMxVg5MBo539x7h/7/nu/shwAVAGfATMzsxrlEmAV0xNSMzmwGMAMa5+19qtI0G3iZIWj1VQLblM7O3CYoD64opBZjZn4ErgCfd/Yp4x9OS6YqpmZhZPkFSKgaer9nu7lOBNUB34MjmjU5EYqBitev8uEaRBJSYms+wcDvf3ffU0Wd6jb4i0nIcEG7XxjWKJKDE1Hz6httIlchX1ugrIi2AmXUHLg13J8cxlKSgxNR82obbXRH67Ay3uU0ci4jEiJllABOBdsCb7v7vOIfU4ikxNR8Lt5ptIpJc/kDwyscq4OI4x5IUlJiaT0G4bRuhT0VbQYQ+IpIgzOxBgpl464Ax7r4uziElBSWm5rM83PaO0KdXjb4ikqDM7D7gRmAjQVJaHOeQkoYSU/OpmEo6xMxy6ugzqkZfEUlAZvYb4BZgM3CKu38W55CSihJTM3H3VcAsIAsYW7M9fME2n+CWwIfNG52I1JeZ3Q3cCmwlSEpz4hxS0lFial53hdt7zGxAxUEz6wo8Gu7eraoPIonJzO4AbgO2ESQl3d1oAipJ1MzM7FHgWoIirm9QVcQ1D3gROFdFXFsmMxtO1S8YEBT3zAUWA1sqDrq7Knu0QGb2beBf4e4MYH4dXRe6+93NE1VyUmKKAzO7CBgPHAKkAwuBJ4HHdLXUcoXLHkzZVz93t331kcRjZpcCE+rRdaq7n9C00SQ3JSYREUkoesYkIiIJRYlJREQSihKTiIgkFCUmERFJKEpMIiKSUJSYREQkoSgxiYhIQlFiEmkiZnaCmbmZvR3vWBrLzG4Lv5dvNOIcw82s3MzujWVsknyUmEQayMyWhz+s+8Q7lqZkZj2AHwHvuPurDT2Pu88C/gncaGYHxCo+ST5KTCJNZxowCLgk3oE00i8Iav79IkbnyqSqoLHIV6gkkUgDmdlygoUf+7r78vhG0zTMrBOwGvgSGOAx+IFhZtOBYUA/d1/Z2PNJ8tEVk0iUzOxSM3OqViNeFt7S8+q39up6xmRmfcLjy80szcxuMbP5ZrbHzFab2f1m1jrs28HMHgj7FpnZYjO7JUJsZmYXmNl/zWxTOGalmf2pgbccLwdaAX+pLSmZWXszuzOMf3e17+FtM/tBHed8mqB48dUNiEdSQEa8AxBpgZYQ/HA9F2gDTAZ2VmvfWdugOjwLnA68HZ73eOBmYJCZfQf4iOA22ntAx7D9PjNr5e53Vj+RmWUCfwfOBvYQLM2wHjgYuBI4x8y+5u4zoojvzHD7Rs2GMHm+T7C8x4awzy6gR3jsSGq/ZVdxrjMInl2J7M3d9aUvfTXgC1gOONCnjvYTwva3axzvEx53giVP9qvW1gvYFLbNA54HWlVrPy1s2wG0rnHeu8O2qUB+jbbrw7YlQEY9v7/WQHH41aqW9kvCc75U85wEV0Qn1XFeI1ifyoFu8f571FfifelWnkh83ejuX1bsuPsqYGK42xu41t0Lq7W/DMwluIoaWXHczDoCNxJcrY1199XVP8TdHwZeBvoDp9YztiEEExWWVY+hmm7h9g13L63xeWXu/lZtJ3V3BxaEu0PrGYukECUmkfgpAWr74b0k3M5w9021tC8Ot/tVO3YikEOwSN2GOj5varg9qp7xdQ23m+tonxZubzOzi82sfT3PC1Ur+naL2EtSkp4xicTPuppXGqGKZ1Sra2mr3t6q2rF+4fa0cGJGJF3qGV+7cLujtkZ3n2pmvwH+F/gr4Ga2kOB52GR3fy3CuSvOGU0ykxShxCQSP+WNbK8uPdwuIpgwEcnH9TzntnCbV1cHd7/NzP5AMJHhWOAY4HvA98zsv8BpdSTfinNurWcskkKUmESSw6pwO8/dL43ROStuCXaK1MndlwEPhF+Y2bHA34CvEUw3f7yWYRXnrOu2o6QwPWMSabjicJsIv+C9QfDM6uQon/VEMh8oAvqaWU59B7n7e8BT4e5hNdvNzICB4e7sRsYoSUiJSaTh1oTbQXGNAnD39cAjBM9s/s/MBtbsE76se6WZ1WvCgbvvIbjtlwmMqOV8Z5nZ8WaWVuN4DnByuLuillMPBDoA8yNM1JAUlgi/6Ym0VC8QvKv0TPg8peKZzG3uXtdMtqb0fYKZeucBn5rZJ8AygkkSvQgSaFa4XV/Pc75I8FLvyQSTGqobDdwEbDSz2cBGggkTRxO8DLwQ+GMt56xIWv+qZwySYpSYRBruYYKH+N8hqN6QHR7/FXVPsW4y7l4CnG9mzxA82zkcOBQoANYSVJn4F/BFFKd9Cvg1cImZ/SJ8B6l6WyHBpIeDgc4EyXkJwTOmJ9y9oJZzjgPKqD1piaiIq4hEFs66uxoYU9dLs1Gc6xCCF4Qnu/u5sYhPko8Sk4hEZGbdgc+B2e4+upHnmgR8Gxji7ov31V9SkyY/iEhE7r6O4Pbk8Y1dwZagwOxDSkoSia6YREQkoeiKSUREEooSk4iIJBQlJhERSShKTCIiklCUmEREJKEoMYmISEL5/0RKpGc5/NzXAAAAAElFTkSuQmCC\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "plt.plot(t[:N], num_sol_imp3[:,0], label=\"Rocket\")\n", | |
| "plt.plot(t[-1],300, '*')\n", | |
| "plt.xlabel(\"time (s)\")\n", | |
| "plt.ylabel(\"height (m)\")" | |
| ] | |
| }, | |
| { | |
| "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 | |
| } |