Skip to content
Permalink
eed6ab4919
Switch branches/tags

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?
Go to file
 
 
Cannot retrieve contributors at this time
579 lines (579 sloc) 107 KB
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Initial Value Problems - Project\n",
"\n",
"![Initial condition of firework with FBD and sum of momentum](../images/firework.png)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We are going to end this module with a __bang__ by looking at the flight path of a firework. Shown above is the initial condition of a firework, the _Freedom Flyer_ in (a), its final height where it detonates in (b), the applied forces in the __Free Body Diagram (FBD)__ in (c), and the __momentum__ of the firework $m\\mathbf{v}$ and the propellent $dm \\mathbf{u}$ in (d). \n",
"\n",
"The resulting equation of motion is that the acceleration is proportional to the speed of the propellent and the mass rate change $\\frac{dm}{dt}$ as such\n",
"\n",
"$$\\begin{equation}\n",
"m\\frac{dv}{dt} = u\\frac{dm}{dt} -mg - cv^2.~~~~~~~~(1)\n",
"\\end{equation}$$\n",
"\n",
"If we assume that the acceleration and the propellent momentum are much greater than the forces of gravity and drag, then the equation is simplified to the conservation of momentum. A further simplification is that the speed of the propellant is constant, $u=constant$, then the equation can be integrated to obtain an analytical rocket equation solution of [Tsiolkovsky](https://www.math24.net/rocket-motion/) [1,2], \n",
"\n",
"$$\\begin{equation}\n",
"m\\frac{dv}{dt} = u\\frac{dm}{dt}~~~~~(2.a)\n",
"\\end{equation}$$\n",
"\n",
"$$\\begin{equation}\n",
"\\frac{m_{f}}{m_{0}}=e^{-\\Delta v / u},~~~~~(2.b)\n",
"\\end{equation}$$\n",
"\n",
"$$v = uln(\\frac{m_{f}}{m_{0}}),~~~~~(2.c)$$\n",
"\n",
"where $m_0$ and $m_f$ 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 function `simplerocket` 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": [
"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",
" \n",
" dstate = np.array([state[1], u/state[2]*dmdt, -dmdt])\n",
" \n",
" return dstate"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"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\n",
"\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",
" \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": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Height of detonation: 599.26 m\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAbcAAAE0CAYAAACvs32dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOydd3hURdfAfycJEBIIEHoHUURAFBOaICBFulQh8CKgiAXlFVEUfamCgi8CL5aPYgOlVwGpgnSQojSliAVBOoQWWiCZ74+5u9kku5vdtE2Z3/PcZ/feOWfu2bt399yZOXNGlFIYDAaDwZCV8PO1AQaDwWAwpDbGuRkMBoMhy2Gcm8FgMBiyHMa5GQwGgyHLYZybwWAwGLIcxrkZDAaDIcthnFsGQ0S+FBElIjEiUsoLvYIiEm3pzkhFe5pbdSoRKZZa9aYUEQl0sCvC1/ZkNkTkRYfrp0Tkaw/19iTQq53WtvoaJ9dKOfxGL4nIzyIyXkTu87WtNkRkjGXjYV/b4iuMc8t42P5k/IBuXuhFADms99+kqkWZEBGZY/24V/nalkxCBxEJdicgIlWBh9PJnsyAH5AfqA68BhwQkR5peUIROWPd14PS8jxZAePcMh4bgb+t9097oWeTPQ2sTVWLDFmdq0Aw0D4JOdsf99W0NSdD0wjIa20FgZrA/4BYIBfwhYg85DvzDDaMc8tgKJ0yxtatWNWTH4rVHVLL2p2plIpJK/syCkqpW0opsbY5vrYnk7PAenX5MCUijj0J89PcoozLTaVUlLVFKqV2KaVeAwZb5QFAfx/aB4BSapD126jka1t8hXFuGRPH8Q9PWm+OMh6NnRgMDti6sRuLSHEXMo2BksBlYFm6WJW5GA/csd7X86UhBo1xbhkQpdRvwA5rt5v11OyO7tbrXqXUAWcCIpJDRPqIyBoROWsFn5wTkZUi0kVEJCU2i8jj1jjXCRG5LSKRIrJdRAaKSJAH+rlFpJ+IrLXGFW6LyGkR2SEi74rIAwnknQaU2Ab/gS7WoWZOAgEmW7Kbrf31Htj3hSX7jwffh01nv6WzxAPZWZbs707K2ojIt9a5b4vINRH5U0Q2iMgQEanoiT1uOAzsAvyBf7mQsT1AzQNuu6vMutcaicj/rGCLyyJyR0QuiMgmEemf1D0hIg9Z1/w3EbkhIrese2uXiEwQkQYu9NL6WjlFKXUb+MvaLeLmc5UWkZdEZJmIHLdsvCEiv4vINBEJc6E3x7qvi1qHRju5rx1/B0kGlIiIv4g8IyLfW/8F0dZ/w/LU+E/wOUops2XADegLKGtr6kaunoPcABcy5YFfHeScbd8CuZ3oNneQKeak3A/4vyTq/hOo6OYzhAPHk6jjxwQ6gQ5lEQ7HX0yiHgVMtmR7WvuxQHk39gWhx5kU8J4X3+FASycaKOhGLg9w3ZIdkaDsMw8+z4fJuL8cr1MxoJ/1fp8T2WAgyiqvl+CeqO1E/i0PbP4FKOnCth7A3ST0dzvRS49rlejzOsgdtWSOuZG5mYR9MTj5HQNzPPhsjr+DMdaxwy7sKAhsS6K+VUAeb69XRtlMyy3jMgf9pwjuuyZtrbYYYFbCQhEJBTYAlYHz6PGASkCo9ToE/STeFvg4GXYOB16y3q8HHgcKAxWBoVbd5YHVIhLixL6KwDqgNPoPfiQ6Iq8guhusuWXXZQ/t+Rw92L/Q2l9LXACAbXvVKpuPdlqCdnSu6GTpAXzloR2gv49YdBRrZzdy7dEOFOLGWxGRVsBzDnU9BpRCP72HAV2BRcAtL2xyxRy0Q6kmItUSlHVAO7i/gK0e1HULWAo8A9QGyqJbMw8Dg4AzQBUcPqsNESkMTEa3Inei78vyQCGgKtAKmARcSKCXntcqESISiP6cAAfdiP4GfAA8gb4GhYB7gBbAEvTD4oci8ngCvZ7oe/CctT+MxPf1AjzA6nlYCNSxDk1GR3wWRD9o2rqpmwHTPakzQ+Jr72o21xv6x6iAa0CQk/KcQKQls8JFHVOt8kvAPS5kWhH3tPZggjKXLTe0Q7pjla0BApzU3cFB/30n5RussutAmJtrEZBg32nLzaHc9qS7KolrPMWS+wsQFzLrLZmNyfgO11q6W93IrMJ56/RT6/j2NLi34rXcrGPLrP2xCWS/t46/6+SecNmScXPuMtY9rYA6Ccqeso7fBkK8qDO9rpXTzwu84yDTKgXnmmjVsdpF+RmrfFAS9bhsuaGnDdlsHepC/yMHmSdS+5qmx2ZabhkbW3BIHpyHabcGCiSQtSMi+YkL3x6ilPrT2UmUUsvRXRTg3dy6nujoMIB+Sqm7TupehHZ8AL0d+/GtFkIDa3eUUuonVydyVncq8bn1Wg5omLBQRMoRZ6M3rTYbttbJoyJS3kn9RYEmCWRt2K7tyWScNznY7iH7OK+IlECHv0MqzZ9USh1HPzAANE1QbPvMV9EO0FPS61rlFpE81lZARMJFZALwrlU+xvo9JRdbS6mBiORMmakusbVwjwPvu5AZhH4gdpTPVBjnlrFZTlz3i7OuSduxq+gujYQ8hp57A7DJ4UeZaAP2WXLhXthniwrbr5Q64kbOFjpeBN1daaOJw/tpXpw31VBK7QJsQTi9nIj0QndbXiN5IfAL0eMs4DxYIwLdBXcXmJugbI/12k5EXpEkJlmnAsuAK0AJdHQk6G5vP3Sr8qinFVn31b+tYIXTVkCIPfgBaGOJ3p9Ada/1WgiYIq6jNxOSXtfqB/S9cA3da7IL3dV/HWislHo7qQpEpLaITBWRX0TkqojEOlwX2wNeLuK6OVMNEfFHdxUDfOvqoVEpdQNYYe0+ltp2pAfGuWVglFJ3iPvDayIO6a+ssbSW1u4CpdTNhPrE/+PYR9yP0tlmGzcr7IWJnowxgA5mSagDUMF6PaeUOu3FeVObL6zXTiJiG1vDamXaxuLmKaWue1uxUuoacQ8e3Z2I2I6tVkqdT1D2JbAf7fw+Bi6IyA+io0ebiEgOUhGl1C3iHLitxW97gPK41SY6svVXdBdbE3TASi4X4vkS2HAIPQYE0Ac4aUVcfiQinazeCGek67VyQgjwPxEp5E5IRMYC29GfrQp6rMxVVGI+F8dTQiH0+Cl4/rstloatyDTDOLeMj62ryB89KG6jC3rMzVEmIcn5cQR6IWtzBFFJyDl2L+V1eB/ipNwXzECP8QShx3xsNER3V0LyuiRt2BzD/SJibxlbwTThCWTsKKWigfro8ZOz6O/mcXQQ0PfAGSu8PTX/uG33UnsRqYcO4ogmcavSKZYt36LH1a6iu+seQ4/P5icu+GGRpRLgpJqX0eNcB9F//NXR0ZzzgbNWyHy8cPt0vFZ1lJU8AP37qovuYQF4EDcPASLSE3jD2l2P/j1XRj9Q2q6LY8+Js2uTUhx/f8n93WYKjHPL4CildqLnIUH8rknb+7+BTS7UbTdvLJBLxWX0cLd5k9HAdvPnSULOsfyak/c+/eEopS6i/5BBR/iR4P0RpZQnUYKuWENclJtj6832/ho6utCZbVesrq7iwEPAC8BsdPdhKNp5pObE/S3AMfTTvW38Z4V1jTyhKXFdz22VUsOUUluUUv9YnyVKKRWFm3tGKRWrlJqilKqCbul3RU83OY5+oOsJbLW60x310vVaKaWuKqW2AU8S5+Cai4iryNiXrdf16C7MOUqpQ0qpCw7XJa1bSI6/v+T+bjMFxrllDmxPg9VFpLKIVCAujPcbZYU3OcEWQOIHJAzvTg2OWa+Vk5Cr4kQHwDZhuYgXYytpha1rsp6IVLC6Jztax1LSarMFw9hShEVY4x4QF7yz0EW3smMdSim1Xyk1VSnVDR3mbnOIESKSKmmWrHvJFthyj/XqTSCJLbHyGaXUBjdyVT2057jlBF5GTwmwpbm6F+fdvOl2rRzOFws8jx53Axjp8B07YkulN9fNb/bB1LTNCReIs9PT3+1pq2WcqTDOLXPwDbr1BbrF1j1BmSvWowMVIH6LJLXYYr1WE/fLfXSyXs8pnX3FhmOCZ3fzzJKDLRWSsz8ZZ6wlLmF1L3S3bxB6/mBqPO3bHEZR9PhpHeLGHL1eosh6yh/jcOgBV7LJwPHzXgK+80LXNrbm8rqLSGN00IpXWE5kNHEBOh595jS+VrZznEKHz4NuuXZxLLeiT21dou7uSacO2wFv7+t4KJ13dru129aFE0ZEcqOnCEHc7zxTYZxbJkApdQK9WgDop33bD2BHAmeRUO8CcX9UL4hIG1eyoKcOWKHpnvI1cc7zI2c/FBFpi54XBXFh9zb79qPnuQEMFpHqbmzzdvzB1o3m0Z+o9SRta6H1BJ613q9KjWAXKyrT1r3cnbjIyZPEhcXHw4MWRgWH9552GyaJFRVZEe0EHvLyqd3WW1DYcuDxEJECuEkWICL3JDEuVhLIbb23f2ZfXasEjCOu++4dx2kvlmO2PTy1daYsIi+SdGSiV/e1C2y9FGXRIf/OeJ+4aUafpeBcviOtJtCZLXU3dGsiYXqcvh7ohQJ/EJfaZwp64L2IVVYR3bL6Ch0A0DqBflLpt951KP8ePSesIPrP5D/EpRv6CyeTcq3z21JbXUNnPKmG/mEVR4ekjwdWJtBLahL301ZZLDoyrTB6gD4A8HNxrUpb18jxGndIxe9wsFVnFDpbjAL+60b+R3SU6zvoP73i1nWphA4/t123Y0AOL21JNInbC12Xk7jR0Xi2VF0n0U68tGV7F+AIuvVxBCeT7NEtrFPABHQ0cDl0IEp59LSJw5beHaCyD66V20nrwHuu7h3iJlYr9IOhLSvIQ2iHH0P8NHnOUpt9YZWdtX4bIQ73tTg5l7NJ3H7EJSaIBT5B/+ZCLZumOdiwKLXu//TefG6AkwufG3gTPX/kMnAD/cc4H6jrQqcbsBk9cBwF7EYP3jr9E0upno+ui2P+QVsGh1APdUujEzEndI7OticS6KZHbska6D9Cd3V4lFsywfVyla9yshtbVjnIncfLP8Ikvody1p+Joy3V3Mj/6MH3dRYIT4YtaeLcrPKeJH5IsG130IEeTjPIEN8BuNqigd4+ulZJObdQ4hzpTwnK8qLn8bmyby/waBLXNoy4rEAJN5Nb0vEz+tqABBe8PHHJR8+i5wfNQ+eYiwYGO9Gxpd25iR4bWOxwcy0C/F2cK1l6Pr4+3zjceF49UaGdUCd0/rnj1ue+jX5KXocOUU6UPJgknJuD3OPocPF/rHovofv2B+IkdZgT/TyWDZvRXS/RaIe3Hd2auz+BvFvnZsmURqcf+x2dU9AT5/aUg9yENPgONznUv9+D38ML1nU9gA4GuIOePLwN3TIukEw70sy5WTIN0JOAI6374Tg652Mtq9yVcwtF5+Gcgn7YPGXdC9esa/BRwnshna9VkunGgFEO8i2d3Oej0C1X2+9kN/qBPhDd0kzq2tZF/1/Zro3Xzs2S8UePxX+PfpCLRv/vrkC3kp2mo8ssm1gf0udYGQX2obuzRgIjlZ7EbCsviM6s/pvDsY7oP+szQH1lZVCwxo3Wo8cM+iulJiY4V7L0DFkfEWlN3Hpl1ZSLJYQMBkPGJiM5t9Howc2vlVI9PdTZjW6m91RKfZ2grAE6WOEMemmN2JTqGbI+IrIQnex5t1Kqhq/tMRgMySNDODcrtcspdD9wZaVT8CSlUwo4gW5K51dO5gmJyD/o6Kq6Sk+2TLaeIetjJTb+DT04/7xSKnNGiRkMhjRJ75IcwtCO7YRS6pCIPIrOeF8Q3YJapZTankDHFjb+qzMHZbEL7aSqE5f1Prl6hiyINf/IH90d/hn6N3GGVMqAbzAYfENGcW62WflHRWQaiSf0DrW6i552cEi25UP+xjXHE8imRM+QNZlFggm36PHWNFnU0mAwpA8ZxbmFWq/10U/RH6Izg1+0jv0fOhXSVeIm19rynrnL1G7LreiYuzC5enZE5Hl0uh2Cg4PDKlVK1Ww+hnSkQIECXLp0CT8/PwIDAylWrBgFChSYEx4ePidpbYPBkFx++umnC0opb1Yh8YqM4txsmVICgM+VUgMdypaKyCn0dICeIjJK6UU3bbP/vR00TK6eHaXUVHSIOeHh4Wr37t3JrcpgMBiyJSLirvcsxWSU9FuOGacTDeIrpXajF/HzI261ZE8y0tvKnGWi91bPYDAYDJmEjOLcjjm8/8uFjO24bcFOm05ZN/WWdlJ/cvUMBoPBkEnIKM7tZ4f3BV3I2Fa4tY2H2ZaVr2JlsHZGjQSyKdEzGAwGQyYhQzg3pdRJdO5D0MlA42FlEn/E2t1t6ZxAO8WcxF892abTAL2O0xnilnhItp7BYDAYMg8ZwrlZvGe9DhUR24KHiEggMAm9pPtPxHc4o63XD0TkXgedIugIS4AxTrKMJFfPYDAYDJmAjBItiVJqmYh8iE6eu0NEdqCnAtREr110EuiqHFKqKKUWiMgk4CXggIisRSdLtS0F8S16OYeE50qWnsFgMGRqZs+G8+fhX/+Cgq5GgLIGGSL9liMi0h7oh84OEoSeUL0U3ZI670KnG3qpmgfR8+QOA18Ck9y1vpKr54gnUwGuXr3KuXPnuHPnjls5gyE9yJEjB0WKFCEkJMTXphjSkdiYGOaHleWpfSfxy5kTli+HJk18Zo+I/KSUCk+r+jNMy82GUmoxejkHb3RmoTNNeHuuZOl5w9WrVzl79iwlS5Ykd+7cOCzOazCkO0opbt68ycmTJwGMg8tG+P38MwEnTtK7DBT0i+XDWrV8bVKakuGcW1bj3LlzlCxZkqCgIF+bYjAgIgQFBVGyZElOnTplnFt24ssv6RgJHSPhbs9ukNdpAqYsQ0YKKMmS3Llzh9y5Xc04MBh8Q+7cuU03eXbixg2YFddJFdD7OR8akz4Y55YOmK5IQ0bD3JPZjMWL4epV/f7ee6FePd/akw4Y52YwGAxZnHe/GsJJ28jIs89CNni4Mc7NYDAYsjCbVsxl2Lq/KOUHlSrD3X9187VJ6YJxboZkMW3aNESEqKgolzLz58/nySefpGTJkuTJk4ewsDBmz56djlYaDIZB77+u30RB5I2cBJRxl1Y362CiJQ1pxvjx4ylfvjwTJkygUKFCrFixgm7dunHhwgX69evna/MMhizP7etXOZDzpF7nJAqead3S1yalG8a5GbwmJiaG6OjoJOWWLVtGoUKF7PuNGjXi1KlTjB8/3jg3gyEdWDTjP0Q1AOpBwX1+jPjgG1+blG6YbklDkvTq1Yvw8HC+/fZbqlSpQmBgIGfPnk0kN3bsWAIDA1m6dClAPMdmo3r16pw7dy7NbTYYDPB/h7/Wb/zh1bCGBAa5W8Yya2FabgaPOHbsGG+++SZDhw6laNGi7NixI175yJEjGT16NEuWLKFZs2Yu69m2bRuVK1dOa3MNhmzPgS0L2ZJfh/8HxMBzPcf52KL0xTg3X5ARwnC9zCl68eJF1q5dy8MP6wUbbOmbAN555x0+/vhjVqxYQcOGDV3WsW7dOpYsWcKXX36ZLJMNBoPn/N/SoRCs33eIKkXxCg+7V8hiGOdm8IiSJUvaHZsjAwYMYN68eaxevZpHH33Upf6xY8fo1q0bbdu2pVevXmloqcFgOHXsN6ZMOQhVgHB4qdHrvjYp3THOzeARRYsWdXp84cKFhIWFUbNmTZe6kZGRtGjRgjJlyjBjxoy0MtFgMFi82r8L6iqwHXL9JtSfkP0CuExAiS9Qyvebl7hK1/Tdd9+xd+9eevToQWxs4lWCbty4QevWrYmOjmb58uUEBwd7fW6DweAFMTFs+nm/fffJxx7Gz9/fhwb5BuPcDCniwQcfZOXKlXz33Xe8+OKL8cru3r3LU089xdGjR1m5ciVFihTxkZUGQzZi2TJOnohlRHEoXlr4+KMFvrbIJ5huSUOKqVmzJt999x3NmzcnJCSEDz/8EIC+ffuyYsUKJk6cSGRkJD/++KNdp3r16uTKlctXJhsMWZfx4wkAhp6GoYPegtL3+Noin2CcmyFVqF+/PosWLaJt27bkzZuXYcOGsWbNGgBeffXVRPJ//fUX5cqVS2crDYYszq5dsHmzfh8QAK+84lt7fIioZIy/GOIIDw9Xu3fvdll+6NAhHnjggXS0yGDwDHNvZkG6doU5c/T77t3hm4ybkUREflJKhadV/WbMzWAwGLIAO9YtofDuObx3L9zxAwYM8LVJPsU4N4PBYMgCvDjoOS78DoN/hwerBUH16r42yacY52YwGAyZnD8O7GDvvgv2/a5N2/nQmoyBcW4Gg8GQyfly/gB4FqgMuYsIQ0Z/7WuTfI6JljQYDIZMzOWzf/PJ3W1QHOgMXxTtly0nbSfEtNwMBoMhE/N/U/tw1Zoyev/VnHTuPda3BmUQjHMzGAyGTMr1S+eYcH2tff/te3rgnyOnDy3KOBjnZjAYDJmU519szoW/FSgoe82fbr0n+tqkDINXY24ikh9oCFQHigL5gUvAOeBnYKNS6nIq22gwGAyGBJw5/gezl+2Bm0AJeK5XS3IEBvnarAxDki03EfEXkadEZD1wAVgIDAaeBzoDL1j7i4ALIvKDiHQSETOimUUYPnw4hQoVclrWq1cvwsPTLMmAR5QrV45p06YlW//YsWOICCLC1q1bE5WPHDkSEfE6XVh0dDTDhw9n7969Ts/33XffJdvmpAgPDzfr5mVxevZujrqp3/tHwb9f/8K3BmUw3LbcRKQrMAYoBQjauf0IHAQigatACFAQqAzURrfsGgAnRGSQUmpOWhlvMKQmefLkYfbs2dStWzfe8blz55InTx6v64uOjmbEiBGUK1fO6UKvBkOyuXyZ8fvP0KMs/PwPPN2qLiGhhX1tVYbCpXMTka1oZ3UemAhMV0rtS6pCEXkY6AV0BWaKSD+lVF33WgaD72nTpg0LFixg4sSJ+Fuh1AcOHODQoUN07tyZ7du3+9hCg8Fi/HiqnIviJ2B3pVJU/XyVry3KcLjrlqwADADKKKUGeOLYAJRSe5VS/YHSwOtWPYZsxPHjx4mIiCA0NJSgoCCaNWvGkSNH7OUbNmxARPjll1/i6TVs2JBOnTrZ921dnt9//z3VqlUjODiYevXq8euvv7o9/9KlSwkLCyM4OJgCBQpQq1YtNm7cmKTdbdu25dq1a6xfv95+bM6cOdSrV4+SJUsmko+MjOSFF16gaNGiBAYG8uijj7Jjxw57ed68eQF45pln7N2ex44ds5ffuHGDF154gXz58lGqVCmGDRuWaMHXH374gVq1ahEYGEjRokXp27cvUVFR8WR++eUX6tatS2BgIA888ABLly5N8rMaMjEXLsCECfbd8CEfEBjkfc9CVsetc1NKTVRKRSenYqVUtFLqf0D2XEwoC3L37t1EW8JVJSIjI6lXrx5Hjhxh8uTJzJs3j+vXr9OkSRNu3rzp9TmPHz/OwIED+c9//sPs2bM5d+4cnTt3jnfeY8eO2ceX/vjjDzp16kSjRo1YtmwZM2fOpHXr1kRGRiZ5ruDgYFq3bs3s2bPtx+bMmUPXrl0Tyd6+fZsmTZrw/fffM3bsWL799lsKFy5MkyZNOHPmDKAdE8DgwYPZvn0727dvp3jx4vY63nzzTfLkycOCBQvo3r077777LgsWxC0sefDgQZo3b06hQoVYuHAhI0aMYNasWfEeAG7evEmzZs2Iiopi1qxZDB48mP79+3P8+HEPr7Ah0/Hf/4LtAadKFejSxbf2ZFBcdksqpa6nxgmUUjdSo56sxvANwxmxcYRHsn0e6cPUNlPjHXt+2fN89vNnHukPazCM4Q2He2tiPC5evEiOHDmcloWFhdnfT5gwgevXr7N3715CQ0MBqFu3LuXKlePLL7/k5Zdf9uq8kZGRbN26lfvuuw+A2NhY2rdvz5EjR6hUqVIi+T179pA3b17Gjo2byNqyZUuPzxcREUHv3r2ZNGkSe/fu5fjx43Tq1IkxY8bEk5sxYwa//PILv/76q922Jk2acP/99zNu3DjGjh1LjRo1AKhQoQK1a9dOdK769eszbtw4AJo2bcqqVatYtGgRnTt3BuDdd9+lbNmyLF261N5NGhoaSpcuXdi+fTt16tThq6++4ty5c+zYsYNSpUoBOsCmXr16Hn9mQ+Zh18bveHX5h3ydF+69Brz7LphsJE4x89wMHpEvXz527dqVaGvdunU8ubVr19K0aVNCQkLsrbu8efMSFhaGu3XvXFGuXDm78wCoXLkyAP/8849T+QcffJArV67Qs2dP1qxZw/Xr3j2jtWzZkpiYGFavXs2cOXNo3Lix00jRtWvXEhYWRvny5e2fE6BBgwYef84nnngi3n7lypXjfa6dO3fSvn17u2MD6NixIwEBAWzZssUuExYWZndsoB8mihQp4vmHNmQaur7cje0HFff5Qd+6haB9e1+blGHxeJ6biIQC5YBjSqlIh+PFgdHAQ8AxYLin43OGzENAQIDTkP+CBQty+vRp+/6FCxf48ccfmTt3biLZxo0be33e/Pnzx9vPmVNnX7h165ZT+fvvv58lS5YwZswYWrZsSY4cOWjfvj0TJ06kcOGko8ly5cpFu3btmDVrFps3b2bUqFFO5Wyf01lrtkIFz4aZnX02x891+vRpihYtGk/G39+fggUL2rtZz5w549SRGeeW9Zj7+Rj++PWa3rkCFZ7oBCK+NSoD480k7rfRASaPoKcBICI5gS1opydoB9dARB5USp1MXVOzFsMbDk9RV+HUNlMTdVVmBEJDQ3nyyScZMmRIojJbgEVgYCCgQ+UdiYyMdDmfzhtatWpFq1atuHLlCsuXL6d///7069ePOXM8m5USERFB69at7Y7RGaGhoYSHhzNp0qREZbly5UqR/TaKFy/OuXPn4h2LiYnh4sWL9i7fYsWKcfjw4US6CfUMmRsVG8uk/WOgJbABihbNyetDE997hji8cW6PA38laJV1AcoDG9Dz4doALwOvoJ2hIZvRuHFj5s2bR5UqVcidO7dTGVsX2qFDh3jkkUcAOHHiBEeOHKFixYqpZku+fPno1q0bGzdu9CqMv2nTpnTs2JFKlSqRL18+pzKNGzdmzZo1lClTxmUrKalWZlLUqlWLxYsX8/7779u7JhctWsTdu3ftY2o1atRg5syZ/PPPP/brunXrVuPcshgrZr/LxoJXoCD4VYU59TLeg21GwwBWhTkAACAASURBVBvnVgrYm+BYa0ABzyml/gTWiEhLoAXGuWVLBgwYwIwZM2jUqBH9+vWjZMmSnD17lo0bN1KvXj26du1KqVKlqFGjBkOGDCEoKIjY2Fjef/99e2skJUyZMoXt27fTvHlzSpQowdGjR5k/fz49evTwuI6AgADmzZvnVqZHjx5MnjyZhg0b8sYbb3DPPfdw8eJFdu7cSbFixXjttdfImTMn5cuXZ968eVStWpXAwECqVavmsR2DBw+mevXqtGvXjpdeeol//vmHt956i2bNmlGnTh1ATzMYNWoUrVq1Yvjw4dy8eZMhQ4akSgvYkDG4G32LgT+NBus56wWq0LBNT98alQnwJqCkADpDiSN1gN8sx2ZjD3qOmyEbUqhQIX788UcqVarEa6+9xhNPPMGbb77JlStX4v2xz5o1izJlytC9e3feeecdhg4dyv3335/i81erVo3z588zYMAAnnjiCUaNGkWfPn344IMPUly3I4GBgaxfv56mTZsybNgwnnjiCV599VWOHj1KzZo17XKTJ0/mwoULNGnShBo1anDq1CmPz1GlShVWrlzJuXPn6NChA4MHD6Zr167xpgsEBQWxevVqgoODiYiIYMSIEYwbN46yZcum6uc1+I7PPnmGQ/l0F37e2zD8pcTj2YbESMJ5Si4FRS4D25VSLaz90sDfwJdKqecc5GYCbZVS2WJWYXh4uHIXHXfo0CEeeOCBdLTIYPAMc29mfI4e2MmDPWtxuzGQB973f4K3B6/2tVmpgoj8pJRKs8S03nRLHgbqiUioFS3ZDd0luSmBXCngbCrZZzAYDNmWdj2ac3svcBAKPe5H/7mzk9QxaLzplvwGCAZ2isg84F0gClhiExCRXOhoyiNOazAYDAaDR/wwezIH917SO7ehz/3tyR2S8nHp7II3zm0SMAudTqsTEA30UUpdcZBpg3aASSfyMxgMBoNzYmNpNHEa7xeHnPmhULkcjBpnxtq8weNuSaVULNBdRIagFyo9qJS6mkDsT+ApIPGiWAaDwWDwjOnTYccO3gb65Qrg78Xf4mfSbHmFuyVvIoDlSqlrjseVUn8BfznTUUr9jF6R22AwGAzJ4dIleOst+26e19+kSkPP86MaNO66JWcB50TkOxHpLSJmJTyDwWBIY24OfA3On9c7pUvDO+/41qBMijvn9hZ60nYLYCpwSkTWi8i/RaRMulhnMBgM2YiJ7/Uj34rpjLtHh6Lzv/9BcLCvzcqUuHRuSqmxSqk66ND+f6ND/usC/wP+EpHdIvK2iCRed8RgMBgMXnH5whkGjvuEO6fhjWPQrVFJ6NDB12ZlWpKMllRKnVZKfaqUaowOJHkWWA5UBt4DfhWRQyIySkTSbEKewWAwZGX6vd6UO7YQvQB4fbjJH5kSvFrPTSl1SSk1TSn1JFAYiADmAyWAd4AdIvK3iEwQkQYiZj0Gg8FgSIpfti5mTtlf4EWgDHRp9Qjhj5kgkpSQ7MVKlVLXlVLzlFIRaEfXBpgG5AZeBX4A/pMaRhp8h4gkuW3YsCHJegYNGhRvQU1PuHXrFiLC559/bj9Wu3Ztunfv7u3HSHWc2WYwJIeYO9H0WdCDu/5AYajTLg8z5nq+ioXBOd6k33KJUioa3VW5XET8gAZAe8Csu5HJcVwq5ubNmzRq1IjBgwfTqlUr+3Hb6tjuePnll+nSpUua2GgwZGY+/G9bfswfBUCOGPi860wCcuT0sVWZn1Rxbo5Yk73XW5shk1O7dm37+6go/QOsUKFCvOOeULp0aUqXNotFGAyOrJz3fwy5vAqsNPNDAhpRufaTvjUqi5CsbkkRKSYij4jIo6621DbUkLG5ePEivXr1onjx4gQGBlK2bFlefvlle7mzbsmjR4/Spk0b8ubNS0hICO3bt+evv5zmB3DJpUuXqFmzJuHh4URGRgLaCfft25ciRYqQO3duatWqxfr1cc9ab731FmXLliXhihgLFizAz8+PEydOALBw4UKqV69OUFAQoaGh1KlTh23btrm0Zc+ePRQuXJjnnnuO6OhoChcu7HSpnVq1atGtWzevPqch63Hj2hU6vvIKd6YCf0CNy8G8PWi5r83KMnjl3ETkKRE5BJwEdgGbXWwJVwowZHH69evH7t27+eijj1i9ejWjRo1K5DwcsXVx/vnnn3z55Zd88cUXHDx4kIYNG3LlyhWXeo6cP3+exx9/nICAANatW2df7LRnz57MnDmT4cOHs3DhQooUKUKzZs3YuXMnABERERw/fpwff/wxXn3z5s3j0UcfpXTp0hw8eJCIiAhatGjB8uXL+eabb2jevDmXLl1yasvOnTtp1KgRXbp04bPPPiNnzpx0796dadOmxZM7dOgQO3fu5JlnnvHoMxqyLi3bPMTN8wquArNh/BOTCcgZ6Guzsg5KKY82dGRkDBALXEIvSurKuW32tN7MvoWFhSl3HDx40OnxYcOGKfQ8TTVs2LBE5QMGDLCXf/jhh4nK+/TpYy+fMmVKovKuXbvay2fOnOnWRk+5du2aAtRXX32VqKxChQpq6tSpLnXfeustVbJkSfv+hAkTVI4cOdTx48ftx/744w/l7++vxo8fr5RS6ubNmwpQn332mV2mVq1a6l//+pc6deqUqly5smrQoIG6du2avXzPnj0KUHPmzLEfu3v3rrr33nvVk08+aT9WsWJF9eqrr9r3o6KiVFBQkPr444+VUkp98803qkSJEi4/j6NtmzdvVnnz5lVvvPFGPJkDBw4oQG3bts1+bODAgap06dIqJibGZd3phat705AO7N6tBpQWJbn1b/TJllV8bVG6A+xWafjf7E3LzZYD5lWgsFKqulLqMVdbMn2tHRF5X0SUtb3hRq6biGwWkSsiEmVNLn/ZCmxxV3+y9AzOefjhhxk9ejSTJ0/m999/T1J+586d1K5dO9443D333EONGjXYsmWLW92TJ0/SoEEDSpcuzcqVK8mTJ25d3J07d+Lv708Hh8mv/v7+dOrUKV69Xbp0Yf78+cTGxgKwbNkybt26RadOnQC9ovfp06d57rnnWLt2LTdu3HBqy/r162nWrBmvvvoqY8eOjVdWtWpVatasaW+9xcTEMGPGDHr27Imfn7nNsi03b0KPHow7odgrUPfBEOYv3OVrq7Ic3vzC7gO2KaU+VkrdTSuDAESkBvAmVgYaN3KfAjOBcHSL8XugIvAJsEBEnKbRTq6ewTVTp06lefPmDB06lPvuu49KlSqxaNEil/KnT5+maNGiiY4XLVrUPnbmiv3793P06FF69uxJ7ty5E9VboEABcuTIkahexy7FiIgITp06ZXd4c+fOpWHDhhQrVgzQzm3RokUcOnSIZs2aUahQIXr06JHItlWrVuHn58fTTz/t1NbevXszd+5cbt68yapVqzhz5gy9evVy+/kMWZw33oCDBwGoJsFsWfwzOQNzJ6Fk8BpPm3jocbZZadmMtM6TC/jVOt9itIN7w4lcR6vsNHCfw/GiwEGr7NXU0nO1JbdbMjPirlvSRmxsrNqzZ4/q3Lmz8vf3V0ePHlVKJe6W7Nq1q6pfv34i/dq1a6sOHToopdx3Sw4ePFjlzJlTrVq1Kp7+lClTlL+/v4qOjo53fNCgQSo0NDTesapVq6q+ffuqK1euqMDAQKfdu0opdenSJfX111+r0NBQ1bNnz3i2ffLJJ6px48aqdOnS6u+//06ke/XqVRUcHKxmzpypOnXqpBo0aODiyqU/WenezCzELlqkFMRtLu657AAZqFtyDVAjeS7UK95Fp/Z6EXAXWfC29fqWUuqo7aBS6izwkrU7yEk3Y3L1DB4gIjz88MOMGTOGmJgYfvvtN6dytWrVYvv27Zw8edJ+7NixY+zatYt69eoleZ6RI0fy0ksv0aFDBzZv3mw/XrNmTWJiYli8eLH9WExMDAsXLkxUb0REBAsWLGDRokXcvXuXjh07Oj1X/vz5efrpp2ndujUHrSduG7ly5WLJkiWULl2aJk2acObMmXjlefPm5amnnmLixIksXbrUBJJkY3b8sJT8Azowv6R1oGNH6NPHpzZlaTz1gkAZ4CzwAeCfFp4WqAXcBWZa+9Nw0nJDJ3NWwG0gt4u6/rFkHk2pnrvNtNw0NWvWVOPHj1erV69Wq1atUm3btlUhISHqzJkzSqnELbcbN26oUqVKqapVq6r58+erefPmqUqVKqkyZcqoK1euKKXct9yU0q3E3r17q5CQELV79267TIcOHVS+fPnUpEmT1IoVK1SbNm1Ujhw51M6dO+PZfPToUQWo4sWLq+bNm8crmzhxourdu7eaM2eO2rhxo5oyZYoKCQlRb731llPbLl++rKpXr64efPBBdfHixXh1bd68WQEqb968KioqyqtrnpZkpXszo3P75g2Vr6S/DvIKRPV5JFipyEhfm+VTSOOWmzcrcR8XkXrAEqCDiKyzHEGsC/n3Pa0bQEQCgelAJDpoxR3VrddflVI3XcjsAkpasrbJScnVMyRBnTp1+OKLLzh27Bg5cuTgkUceYfXq1U7H1QBy587NDz/8wGuvvUavXr0QERo3bsyECRMICQnx6JwiwtSpU7l+/TrNmjVj06ZNVK5cmenTpzNw4ECGDBnCtWvXeOihh1i1ahU1asTveLj33nsJCwvjp59+YvTo0fHKHn74YVauXEn//v25dOkSJUqU4JVXXmH48OFObcmXLx9r1qyhQYMGtGjRgrVr15I3b14A6tWrR8GCBWnXrh3BZvmSbMlLL9fjyqkYvXMbHur8PBQo4FujsjqeekFAgInAHbRDi0VPDUi4xQIx3npZYBy6xdTF4dg0nLfc/m0dX+ymvomWzIcp1XO3ZaeWmyF5/PTTTwpQW7Zs8bUp8TD3ZvqwaenHym8oimdR5EM1bFjO1yZlCMgoLTdgENAP3W24HPgdiPLGkbrCymjSH/hWKTXXAxVb7Pd1NzI22/Kmgl48ROR54HmAMmXMuq0G55w/f57ffvuNt99+m7CwMOrWretrkwzpzJk/99N586vEBgNloEbnIFb+74CvzcoWeOPcegM3gHpKqb2pZYCI5Aa+Qs/T7+upmvXqdqpAKurFQyk1Fb06OeHh4Smqy5B1WbhwIX379qVKlSrMnDnT1+YY0pm70bfo8lF9zhTQIzeFbgrN83QlMChPEpqG1MCbiMCSwKbUdGwW76PnmA1QSp32UOea9eruLrGVXXM4llw9g8FrXnzxRWJjYzlw4ADVqlXztTmGdKZDl6psCtYB36JgdvgY/EK8W/bJkHy8cW4n0S231KY9epyup4hscNyA5pbMS9Yx2+JZx6zXsm7qtaW+OOZwLLl6BoPB4DFv9WvPsm//gM+BCzDSvwlNOr3pa7OyFd50S84FnheRYKWUuzGr5GBbA84V91hbfmt/j/VaRURyK+eRjzUSyKZEz2AwGDzil/XL+e9n3+qdc1B8TS7ePrjSt0ZlQ7xpuY0EfgOWikiF1DJAKVVOKSXONvTUAICB1rGHLZ0TwM9ATuCphHWKSAP0nLYzgH21zeTqGQwGg0dcuULVlwfSsyjgDwEhsPXbH/HzT/WlMw1J4M0VX4qOlHwcOCQif+J6nptSSjVLBfvcMRqYD3wgItuUUr8DiEgR4P8smTFKL56aGnoGg8HgmpgY6NYNDh1iGtCseAAMG075Bx72tWXZEm+cW5MEehWtzRlpHkGolFogIpPQKbMOiMha9By8xkAI8C06EXKq6BkMBoNb3n4bVqyw73b9cLp2dgaf4I1za5pmViQTpVRfEdkCvIwes/MHDgNfApNctb6Sq2cwGAzOmDy4FxETp9uDAhg0yDg2H+PxmJtSap03W2oYp5TqZY21fehGZpZSqq5SKkQpFayUClNKfZqUg0qunsFzTp8+TcuWLcmXLx8iwoYNG5Jd12effcZ9991HYGAgYWFhrFvn2S22detWatWqRe7cuSlfvjwfffRRIhkRSbTVrl072bYashdfTPwPL42ZTukSsLcA0KYNvPeer83K9phRTkOa8d5777Fv3z5mz55NaGgolStXTlY9c+bM4cUXX2T48OHUq1ePr776itatW7Nr1y6qVq3qUu/333+nWbNmtG7dmtGjR7Nz504GDBhAUFAQzz33XDzZ119/3b5QKWDPC2kwuGPf9jU8/5/3IQaijkHTCv6c++YbxCxG63OMczOkGYcPH6ZWrVq0bNkyRfUMGzaMnj17MmTIEAAaNGjAnj17GDNmDDNmzHCpN3bsWEqUKMGMGTMICAigUaNGHD9+nBEjRtC7d29ExC5brlw501ozeMWVc8fpPu9JYqsDW0ACYe6EaUi+fL42zYCbbkkR2WTlfEw2IlJXRDalpA6D7+nVqxfh4eEsX76cypUrExQURKtWrYiMjOT333/n8ccfJzg4mPDwcPbv3w/orr5169axePFiRIRy5col69x//vknv/32G507d7Yf8/Pz46mnnmLlSvdzh1auXEmHDh0ICIh7houIiOCff/7hl19+SZY9BgNA9M0oOr1fnV/y34Ym4NcSxg56gUZtuvvaNIOFu7ZzJWCziHwvIl1EJJcnFYpILhHpakUhbsJ1RKUhE3H8+HGGDh3KqFGjmDp1Ktu2beP5558nIiLCvujn3bt3iYiIQCnF9u3bqV69Oo8//jjbt2+3Lx6qlOLu3btJbjYOHz4MQKVKleLZ88ADDxAZGcn58+ed2nv9+nVOnDjhVM+xXhvDhw8nICCAQoUK8eyzzxIZGZmyC2bIsqjYWJ4f8jBrC8TdI1+17MPrwyb70CpDQtx1S94HjEAnM24EXBORrejJzYeAi+hkxyFAQfTq2XWAuuj8jHfRy8cMTyPbDelIZGQk27dvp0IFPX9///79jB07lunTp9OjRw9AO65WrVpx+PBhateuTUhICKGhofG6+6ZPn+7RatR6RQy4dOkSoFfDdqSAtRbWpUuXKFy4cCL9y5cvJ6lno2fPnrRp04bChQuze/duRo4cyb59+9i5cyf+/v5J2mrIXrzwYg2ml/zDvv+uNKLHy1N9aJHBGS6dm1LqCtBfRD5CL3XTE2hBXL5HZwh6sdFxwP8ppY6lnqlZiw0bNqQoejC5NGzYkIYNG3qtV65cObtjA73QJ0CjRo0SHTt58qS9hZSQNm3asGvXLq/P7zg+BnHOL+HxpPScHZ82bZr9ff369XnggQdo2bIly5Yto127dl7basi6PNOtHtNm/wz1gEbQ+0ZFBn/wva/NMjghyYASpdSfwGsi8g56TlhD4GGgCJAPuAycQ6e1Wg9sVkrdTiuDswrJdTK+ImELKGfOnImO247dunXLZT2hoaHk82LA3dbSunz5cjw9Vy2zhPba5Gy4agk60rx5c/LkycPPP/9snJvBzifDXmTanK16ZwvcL0FMWrPHREZmUDyOlrSSDK+yNoMhWXjbLWkbMzt8+DBly8Yt5nD48GFCQ0OddkkCBAcHU7p06URja67G8ByxteqSahUashGbNtFt3FeMKgFnT0LuwsIPX+8hR2CQry0zuMBMBTCkK952S95zzz1UrFiR+fPn06yZTlcaGxvL/PnzadGihVvdFi1asHjxYkaNGmUfO5s7dy6lS5d2Oz9u1apVREVFERYW5rGdhizMzz9DmzaEXo/mzxvQsGJOJn+xlBLlTKxcRsY4N0O6UrBgQQoWLOiVzvDhw+nevTvlypWjbt26TJ8+naNHjzJr1iy7zMaNG2ncuDHr1q2jQQO9etLAgQOZOXMmTz/9NH369GHXrl1MmTKFSZMm2VtlU6dOZffu3TRp0oRChQrx888/M2rUKGrWrEmrVq1S74MbMidHjkDz5nD1KgBBRYuxc8UWqJBqC6MY0gjj3AwZnq5duxIVFcUHH3zAyJEjqVKlCt9991281pdSipiYGHt3JugAl1WrVjFgwABatGhBsWLFGDduXLzsJBUqVGD69OksXLiQq1evUqxYMXr06MHIkSNNpGQ2Z8cPS5k18F9MPB+lD+TPD2vWGMeWSRDHPwOD94SHh6vdu3e7LD906JDLyEGDwZeYe9M1+3f8QI0WjYm+Au0qw6I/cyNr10GdOimqd/jw4QwfPjx1jMzkiMhPSqnwtKrfhPkYDAaDA+eO/Urdbk2IvgTEwreHYcPHI1Ls2Azpi3FuBoPBYHH++CEafRRGVDulU1MIvNHnSR5/dqCvTTN4iRlzMxgMBuDCiSM0/t8j/JpPT9OVHvBGdGv++8kSH1tmSA6m5WYwGLI9kaf+oOn46hzIpxMQ+MXCjIp9+e8ny3xsmSG5eOzcROQ5EcmdlsZkVUzQjiGjYe7JOP46tJfy9Sqy985NAETBtKIv0O3FT31smSEleNNymwr8IyLjROS+tDIoq5EjRw5u3rzpazMMhnjcvHmTHDly+NoMn/Pr7k1Urv8IV/+KhWnAOfiyUG+e7msy/Gd2vHFu36FXAHgNOCQiq0SkjZgcRW4pUqQIJ0+e5MaNG+Zp2eBzlFLcuHGDkydPUqRIEV+b41tOnmTPsxHcumr9Lq/DM5fr0uuVz31rlyFV8Ca35JMiUhp4CXgWeAJoCpwQkcnAF0op54trZWNCQkIAOHXqFHfu3PGxNQaD7k0oWrSo/d7Mlhw7Bo0b0/3P01wuAv0uw3Od6/PZNxt9bZkhlfAqWlIpdQJ4R0SGAZ3Ra73VAd4DhonIAvRSN9tT3dJMTEhISPb+IzEYMhK//QZNmsCJEwC8EhlA00/+y/0vvOZjwwypSbKiJZVSd5RSM5VSdYHqwBfoxUm7AVtE5CcRedbT1bsNBoMhPVg682N2tKptd2zkzAmLFhnHlgVJ8VQApdQ+9IrdX6EXKxW0w/sMOCYivVN6DoPBYEgpk8a+Qbvn/s1j1y+xPz8QFATLl0ObNr42zZAGpMi5iUgTEVkE/AW8DNwCvgS6AivQC5pOFZF/p9RQg8FgSC5ff9qfvv8Zh7oFd07DYwUgdtVK3T1pyJJ47dxEJJ+I9BeRw8BqoB1wCngHKKWUek4pNVcp1QZ4FLgOGOdmMBh8wucTe/LMuYk6BA6QIPh04Lv4PVbft4YZ0hSPA0pE5BF0AEkEkBvd/bgB+BhYopSKTaijlNohIsuBjqlircFgMHiIio1l1KimDFU/6Mf4GlDwuh/zBkynUZvuvjbPkMZ4Ey1pW9flBvA58LFS6hcP9K57eR6DwWBIEdG3btJvcDhT8x60H3vkcm5WfLmTouVdr8JuyDp443SOAZ+i57Nd9kKvD/CCN0YZDAZDcrlw+jjVHqvI6Ty39aCJQJNLoSwacoC8BUv42jxDOuGNc6ugkpFiw9KJ8VbPYDAYvOXc4QOUb/gwN85aoyR5oGutskwb/Qs5c+fxrXGGdMWbgJLVIjIgKSEReU1E1qTAJoPBYPCeffso1KQ59+SKG/4Pu1OEbz44ahxbNsQb59YE8KSzujLQOHnmGAwGQzJYvhzq1cPv5Cl2nYAiJSCiYxi7d5zFP8AkiM6OpEWgR04gUeSkwWAwpDaxMTHEfDyRHK8PhFj9txOYN4TTX8zFr3lzH1tn8CWp6tysFQLCgAupWa/BYDAk5Ma1K9SsXx65e4l9yuqGKlcOvvsOvypVfGydwde4dW5Oxs6ecDOeFgDcB5QAFqSCbQaDweCU3/Zvo0aL+lw9pWPVOlWFRXlqw5IlkN2X8jEASbfcHHPTKLTjSiqWdj/wZkqMMhgMBlfsWjud9iuf5Wpg3OjH3thgYtetxS8o2IeWGTISSTm3ptarAGvQ6bY+dCEbDZxUSv2ZSrYZDAZDPKZ/+jwvnPmM2yFAJ+BzaB5ekeUrDuLn7+9r8wwZCLfOTSm1zvZeRLYCGx2PGQwGQ3oQffM6A9+tx0eBe+3/Wvn9hGkThtC2xwjfGmfIkHizEvdjaWmIwWAwOOPgz1uo2+5xLte+C1acSJUrufj26eXcW93MOjI4J8XruRkMBkNasXjy+zzY4DEun7gL3wJnocOVEmz/z5/GsRnc4rLlJiLvWG8nKaUuOex7hFLq/RRZZjAYsi+xsTBuHLWGDsY/0Jo4ewfanKjE/G8P4OdvcrEb3OPuDhmFjpBcAFxy2E8KseSMczMYDN4TGQm9esGyZZQAZgVBRB4Y2rc7Qz/4xtfWGTIJ7pzb+2gndSHBvsFgMKQJ+5bN4KF+g+Hvv+3HOt1fh8vTviBPxQd8aJkhs+HSuSmlBrvbNxgMhtQiNiaGpzo8wqJV+xlZBux/Nq+/DqNHkyeHyQ9p8A4TUGIwGHzK+eOHuL9OCIuW7odoGHoBdpYNhsWL4cMPwTg2QzIwzs1gMPiMNfNGU+2Tqvxe54ZOuQ4E5hTuTPkK2rXzrXGGTI3Hzk1EXhKRaBFp5UamtSXzXOqYZzAYsiK3r1/l9XfCaXboHc4Ex0Io0AoeDi/MqV9PUbfZU7420ZDJ8abl1gGIBFa6kVlpyXRKiVEGgyHrsmLeFKo+U5jxuX6yHytyQ1jZeSR7dp0jf6FiPrTOkFXwZrJIJeCAUsrlWm1KqRgROYBesNRgMBjsxMbE0KPro8xcvBNyAWWBPNDicmG++vcPFC3vyVrIBoNneNNyKwyc9UDuHGDWnDAYDHGcPs2ZNs2YvWIn3AWug9938FFQR5aPO2McmyHV8ca5XQFKeyBXEohKnjkGgyFLoRTMmgVVqlBi5TpGhOjDuQoK8waNo9/ABYifiWszpD7edEvuAR4XkQpKqT+cCYhIBeBRYFNqGGcwGDIvV47/Tr4Bg2DhQvuxwafhfOuHGPb5ckKLlvShdYasjjePTNOAHMC3InJfwkIRuRed2tTfkjUYDNmUt15pR2i1+/hiV5xjo2xZ+OEHJi7baxybIc3xpuU2F+gOtAR+FZEtwGGr7H7gMau+VUqpGalqpcFgyBREnvqDhh3DOPDjFQBeCoa2uaBQjz4wbhzkzetjCw3ZBY9bbkophZ4OMMk61BB40doet45NAtqnon0GgyEToGJjmf/FOERUpgAAIABJREFUACpPrMiBalfs/ywxV2HruKEwdapxbIZ0xat1I5RS0cDLIjISaIwO5gX4G1inlDqTyvYZDIYMzskju3h58pMsyX8GgtBbfbjnXB7WLdxCuUoP+dpEQzYkWYsiWU5sZmoZISI5gProLs+6aKdZEDgPbAc+UUptcKPfDXgJqIYe8zsMfIVei87lvLzk6hkMBrh7J5oeXeuyJHA3NxxG4Utc9+PTngNp12uM74wzZHsyyop/DYDvrfdngJ+A6+jJ4B2BjiIyUik1NKGiiHwK9AVuAeuAO+hW5SdAYxF5SikVk1p6BoMBNi35hjZ9n+HqqRidOqsckANevF6ZMa+vIF/RsknUYDCkLV5PMBGR+0XkUxH5VUQuW9uvIvKJiFRKph2xwEKgvlKquFKqtVKqi1LqQSACiAGGiMjjjkoi0hHtoM4A1Sy99sB9wCH0+N8rTj5DsvQMhmxPdDSMGkWeXs9y9aL17BcJhTb5san6R0z676/GsRkyBF45NxHpBexFB5E8AIRY2wNoZ7FXRHp6a4RS6gelVCel1GYnZXOJm1rQPUHx29brW0qpow46Z9HdjQCDRCTh50yunsGQfVm/Hh56CIYM4ZHLd/lXUcAPHq1biqOLTvLYk/18baHBYMebVQFqAJ+hF6ZYDLRGO7XKQCt0yysH8Jklm5rssV5LOdhTCggDooH5CRWUUhuBk0AxoHZK9QyG7Mr+H9fxVuvK0KgRHD5sPz6t8COs+OZTtm45YZIdGzIc3oy5DUQ7w+5KqdkJyg4DK0WkKzrQ5A2gS+qYCOjuQoDTDseqW6+/KqVuutDbhU4HVh3YlkI9gyFbcetGFN271mPh6n0g0KQIND0H5MkDI0cS0K8fLfz9fW2mweAUb7rd6gE/OXFsdqyyXejIx1RBRIoBvaxdh3QHlLde/3ajfjyBbEr0DIZsw4+rPqfWW4VZuHkf3AZuQbf8oCK6wJEj0L8/GMdmyMB449wKAr95IHcUHT+VYkQkAJgB5EPPo1vmUJzHer3upgpbAmfH2aPJ1XO063kR2S0iu8+fP++mGoMhc3H++CGef/MB6uzow/5Ct/TkHCBHARjY5Tlk9hwoUcK3RhoMHuBNt+QloIIHcvdYsqnBZHR4/gkSB5OI9aq8rDO5enaUUlOBqQDh4eHJrsdgyChEXYnkzTfaMKvQNq4Exx3PXQ7adq7GlElrCQkt7DP7DAZv8abltg2oKSJtXQmISBt0EMbWlBomIhOB3uhw/cZOsp9cs17z4Bpb2TWHY8nVMxiyJCPe/Beh5Qsy6YttXLkad/zJy8U4GLGZ2XP3GcdmyHR449zGW6/zReRLEWkgImVEpLT1/gtgAXrO2njX1SSNiIwD/o3OUNLYMVzfgWPWq7tJNbb15445HEuunsGQtTh0iNjmzZj4zSzuXEL3ZayCe68EsKzCEJZMOE25qvV8baXBkCw87pZUSm0Rkf5ox9XT2hwR9GTr/kqpZLfcROS/wADgItBUKXXQhahtekAVEcntIvKxRgLZlOgZDFmDS5dgxAj49FP87t7ls1DoJEBOaFOpMvNGbicwOMTXVhoMKcKrScpKqY+Bmuggj+PoBeNjrPdfAzWVUp8k1xgRGYOecnAJ7dj2ubHlBPAzet7dU07qaoCeF3cGnZ8yRXoGQ2bnxrUr9O1Z///bu/fwKKrzgePfNwQI93AHDSgiKlqriNWKraBUoAVFAaX1QrQCLaBiKSgq6oOi4AWFCqJAEAQvFPGCiiKKoKjcRFBRAREwXKSCBAkSIMn7++NMzJJfsuYyu7PJvp/nmWeSObNn3jnZ2TczO3MOWSe1hPHjITsbgJ57hRs7n8oXS5cw7411lthMhVDiHjhU9VNVTVXVFqpaVVWreD9fp6qlPtPxRhq4DcjAJbbi1DXamz/oDZaaV1cj4Anv1zGFdIJc2tcZU+5obi633tid5ObJTHrmAwY1DLnfq317WL2ax99cx2ln+/YEjzGBi4mOk0XkUmCE9+s3wE0iUtiqX6vqL12Nq+qLIjIJ12XW5yLyDvkdINfGjQz+/84kS/s6Y8qbD16bwK2LbmfZikz3byMw/Qe47+QUjrl/HPToAYUfa8aUazGR3Dj6ubizvakwS4CjxtFQ1YHeqOCDcKML5A1dM40wQ9eU9nXGlAdfLX+d21/o58ZYS8a9w9cCAl3OPoWazy6Beo0CjtKYyCkyuYnI5DLUq6r6jxKsPJ38zpFLs7HngOei9TpjYtWajxdy/ZC/sbbjHjQ5f3mVJLikZysevPN5Wp7WNrgAjYmScGdufctQrwLFTm7GmDLKyKD7lecwb/FGd4G9Hu7WL+Dqn45n1PUz7bZ+E1fCJbd+UYvCGFM6mZnwn//Aww8jdTJcYgNYDB2OT2ZszwmcddHVQUZoTCCKTG6qmhbNQIwxxXd4/z6qTJ0Go0eD17/ptP3QsAZUqSaM6NeH2+9LI8E6NzZxKlZuKDHGFEPmvh/pf0Nn5ry3igWV4aKQfrvrtTiRt/tewR8H30WVpGrBBWlMDChVchORmrg7GhsC36nqcl+jMsYcJefIYWZNHsSAsWkc3Oz66u7XGjbtApo3h7vvhtRUOiba/6vGQAkf4haRWt5dlLuBd4EXCLlxREQGiMh3InKOv2EaE5+yD2fxzMT+nDq8JtftnsrBkEEotu6AnY+Ogg0b4IYbwBKbMb8o9tEgItWBxbjRqXfjurDqVGC1t4GJwOXACn9CNCb+/Lx/H/fc3puXkt7l21rZrmsBgNZQ+XjodMopTJ+2gAZNmwcZpjExqyT/6v0bl9ieB/qr6gEROepBZ1XdJCIbgYt8jNGYuJF78CB9UjvwwoIV5PyE62LAGzK3ThYMSbqQW1Y8Q+2GKUGGaUzMK0lyuxLYCdygqllh1tsKnFqmqIyJN1lZkJZGwpgxLMzd5hIbwPtQt6swpNpF3PSvadRpZGdqxhRHSZJbS2DBryQ2cJcsG5Q+JGPih2ZmIlOnwkMPwc6dADzQ0PWgINWgyzEnMXvoe9Sqf0ywgRpTzpQkuR0BqhZjvRQgs3ThGBMfNn6+gv4392blli3s2gY1svPLbpBGfPG3k7n9oRk0SmkRXJDGlGMluVtyA9BGRIpMcCKSDJwBfFHWwIypiL778mNuHHoGJ517LosXb+HAFrijpVfYtCmMGwebN/PYc+9bYjOmDEqS3OYCjYEHwqwzCqgJzClLUMZUNOs+epXUoS1p+UI7Jtb67KhvpWdrAkyYAN9+C4MHQ/XqwQVqTAVRksuSjwOpwC0icjYu2QEcJyL9cKNadwTWAVN9jdKYcmrSw0OZ/mEaK9pk/HLXIwDtICld+HvXjoyd8DJUrxlYjMZURMVObt6t/51wSe2PQF4X4x28SYA1QHdVPeRvmMaUIzk5fPjUw/x51Aj278xx46n9FjdiINB+bx1uO28IXcaNQBJK1I+CMaaYStSlgaqmA+eISDfgL8AJuEM2HXgTmGuDfJq4lZkJ06fDuHG03rKJzCre8gzgS7isWVNu63wvv+9SltGkjDHFUar+elT1deB1n2MxplxavmgevDyHc2e9DhkZgBtO7fxGsDQdWp1am/F/fYQ/97JRpIyJlnAjcb8IpAFvqaoWtZ4x8eq5p+5j+KOjSd94kFanwIaMkMLkZKZd2ouMnt35XftugcVoTLwKd+bWA9dH5E4RmQFMV9WN0QnLmNiUfTiLV2beyWOfT+GjnfvdAzLAxi3weV04vf6JcMstkJpKq5p2k4gxQQmX3CYBvYFjgOHAcBFZCkwD5qjqz1GIz5iYsGX9Gl58aQQT9yxgS61sqAvUwd0skgHJ9RNJv/UOTh90D9hNIsYErsijUFUH4RJbb1xv/7m4uySnAd+LyBQRaReVKI0JyPOTH+Ck0+vQ4jdtGJb+hktsnsoKf+rYhBemjGFv+hH+ctNIS2zGxIiwN5So6mHcA9lzROQY3HNufYCTgRuAv4vIBlzCm6mq30c4XmMi79AhmDMHJk7k1m3L2LbNW74S6Ab1DwoDqp7PwNTHadryzCAjNcYUodj/ZqrqDlUdraqtgfNxN5vsxyW6McB3IvKqiHQXkUqRCdeYyNm7/nO4805o1gyuvRaWLWPYkfzymjsTmJqcSvodu7lv5AeW2IyJYaV9FOBj4GMRuRnoBVyHe5C7mzf9ADTxJ0RjIic3J4ex9w7gseee4YfMQ2TsProT4xt/TOTZtnXod3U//n7zKBIq2f9txpQHZRqXXlUPAjOBmSJyMTALaOhNxsSs3enreebZW5mU/ibfTD4CXkK77yQYswFISYEBA0jo25fljRoFGqsxpuTKlNxEpCbuhpPrgHa4LrjA9VhiTEzJzclh4dyxzFj+BHOrb+VwItAI1zXWarfO4tq14aXpcMklkFimw8MYE6BSHb0iciFwPe5ZuGq4pHYImIe7ueRtvwI0pqy+Wr2UoXf2Y+GyrzlyNvm9onqqnwktqc8Dw+6l218HBhKjMcZfxU5uItICd7dkKtCc/LO0NcDTwCxV3et7hMaURm4uLFoEkyfz2IoXmb/V62RnNe4aQwKck1GD/i160fumMdSsZ18RG1ORhE1uIlIdN5TNdbhn3MSb9gLPAWmquibCMRpTbBvXfkTL+YtJmJrmxkcDHkiEKVVx1xYOwNU7T2LYNfdxxgVXBhqrMSZywvUtmYZLbDVwCS0XeAd32fFl7xk4YwJ3+GAmd992LdPfmM+urYd59hi4KuRb3wbZ0OP0etQ99TQeGDPDRrg2Jg6EO3O73ptvBqYDT6vqtqJXNya61i6ZzdNvjeFZXcvu+QruRI37a8JVAMnJ0KcP9OvH3N/8JshQjTFRFi65PQtMU9X3ohWMMb9m/dplzHjmLt7M/ZA1yQchyStoA2xyP+7JqUzujCkkXHElVKsWVKjGmAAVmdxU9dpoBmJMkbKzeXfSGPr85352fJsFzci/ruBJSUmgxQXHcs+QUXTs3ieQMI0xscMe5DGx67PPYOZMmDmT+vt3seMgoMBW4EdIqgWXHzyO63//Ty66bAiVKlf5lQqNMfHCkpuJKcsXzWPkQ7dy197DnLdi8y/LzwQaHQP/2wE1miRwU8LF3HbTUyQ3Pi64YI0xMcuSmwlcxq6tzJ19D3enPc+Oz9xNuNmnFugJoEkTJnY+h6ROXej21wGBxGmMKT8suZlAHDrwE2++OJpZa2fyeo3tHEoEmgKfufIl+yGnehKVul8O11wDnTrRy7rDMsYUk31amKjJPnKYyY8OZ8LzT7NJMjh8GW406zy/ARZC/ZTK9OrwB3LHv0il5HoBRWuMKc8suZnIUoU1a2D2bOa/Mo1B639wyxOBLvxyK/9ZGdW4pnEnuqwcQuuzLggqWmNMBWHJzUTEK7PG03bdJprNfQs2bgTgUqBafTi4B8iG+p8J/2zbjqu7Dqf1ud0CjdcYU7FYcjO+2bDyLYbeO5i3P97IoT1K6ikwfePR6/SoW4XVKdW58W/X03/IGBLt9n1jTARYcjNlsnntEma/NprZe953PYZkAHtc2by8Ea1r1oTu3aF3b2Z16gRVqwYVrjEmTlhyMyW2fNE8Rj0ynLWZm0jv6PWfnewVngYsBRIhOakGh+akUbXrpdYNljEmqiy5meLZtAnmzuW/L02m93KvE8fqQAegkvu1ajZ0SWrCiQPaMuz2x2nczHrfN8YEw5KbKZTm5vLaC09wwYYdJL8yH9auBdzQ6wk1IPcA8DNU+ha6NGxE71aX0/2KEdRumBJo3MYYA5bcTAjNzeXTxc8zbPSdLP1kK4f3wh2t4P6Qm0ISgTMawNb6lbn0/HO4Z8REjj/ljMBiNsaYwlhyi3O5OdksW5DG3KVTeOnwWrbUynY3hOx15S9UgvsBqlSBTp2gZ09WdetGQoMGAUZtjDHhWXKLQz/v38eEh4cy8/WX2Jq8l/3tFariJoDWwKdAZaiSVAt97imka1eoXRuAhIDiNsaY4rLkFi/27YM334R58xi9fC6jvvXucmwAtM9frU4WdG1wPM1vOothd4ynXhP7Ds0YU/5YcqvAPn7nZSZNHc3kvbVIeu8DOHIEgMGJMCoRyAZ2Q3I69Kp7Mj3bXsNF3W+hSrWagcZtjDFlZcmtAsm7IeTVJZN5ZMpSft6ZC8AFx0HfI/nrNciG01omklCnDtd1v5KBQx8iqbolNGNMxWHJrZzL3Pcji+Y9ztvr5jIv50vSa+a4L8VqATvdOmk1oC9A27Zw6aXQvTtf/Pa3IBJc4MYYE0GW3Mqj//2PR0cPZtz8V0nfehAuBH5fYJ2TgW+gfrPKnPf7P8KCGZBi358ZY+KDJbfyIDcXVq+GN96A+fNh5UrWNFPSv/PK1/NLcqubJXQ9fDydunTlwnEDSTmhdVBRG2NMYCy5xaidW9bzyINDmffBEnYfyuTHb5TQi4g374GZ3s+VM+CfB86kR7tUzu/Sn8pJ1YMI2RhjYoYltxihubl8vXI+byx6kjd2fcgHCRnkPOkVJsD7jaH9Lu/3SpU4u+359KmXyWU9rqX7VTeRUKlSUKEbY0zMifvkJiJXAQOA3+K6AP4aeBqYpKq5kdz2j7u2M+HR4bz4znwyzttHesMcV1DXW6Ep7qaQXJjSpCrtO/eGrl3h4ouhbl1mRDI4Y4wpx+I6uYnIRGAgkAW8CxwBOgITgI4icoWq5vi2QVVYtw4WLIAFCzjh84Xs+94rOwFoePTqx51WheatGpPasw/X/mMEVE3yLRRjjKnI4ja5iUhPXGL7HrhAVTd6yxsD7wGXAzcC48uynY2fr2DC4/fQaM8u7lz+P9i+/ZeyNsfB4l9WhFqtoFPWsXRt0Ykuf7mZpvecWZZNG2NM3Irb5Abc7s1vy0tsAKq6S0QG4PLOcBF5vCSXJ7MPZ7F84dMsWPE8sz5ezuaFrpuruifCnduPXvfqn+GjenB6q8akXnoF//jXaOsdxBhjfBCXyU1EUoC2wGFgTsFyVV0iItuBY3E32X8Urr4PF85lztwnSK/+Ne9W3cG+JNyD1Cfkr7N3B+ysDk2rJLvvzDp3pm+nTvRt1sy/HTPGGAPEaXID2njzdap6sIh1VuKSWxvCJLfVaz7hD516uR71b+WXUakBaAzUgpo1Emh7cjMODR8Lf+oOifHa7MYYEx3x+inbwptvDbNO3iPSLcKsg6r3wyFgO9AcUjIr0ZmWdD6lK39Y25+mLU4pW7TGGGNKJF6TW94XWwfCrJPpzWsVLBCR/kB/79dDwBcATHMLtpFDGhtIYwPwWNmjLT8aALuDDiJGWFvks7bI12DkyJHWFs7Jkaw8XpNbXmcfGnatIqjqZGAygIisUtWz/QqsPLO2yGdtkc/aIp+1RT4RWRXJ+uN1UOX93jzcrYl5ZfvDrGOMMSYGxWty2+LNjwuzTt5tjFvCrGOMMSYGxWty+9SbnyYi1YpY53cF1i3KZH9CqhCsLfJZW+SztshnbZEvom0hqqX62qncE5FPgLOAVFV9pkBZe9xD3N8Dx0a6j0ljjDH+itczN4DR3vxBETkxb6GINAKe8H4dY4nNGGPKn7g9cwMQkSdwIwJkAe+Q33FybeAVoJevHScbY4yJing+c0NVBwJXA18CFwOX4voaUeAyXOfJpSYiV4nIByKyT0QyRWSViAwSkZhvdz9jF5EUEXlcRNaLyEERyRKRjSLypIic8Os1BMvvv6OIVBORW0VkpYhkiMjPIrJZROaIyPl+x++nSL6nReQBEVFvGupHvJHkR1uISGUR6SgiY0VkmYjsFJHDIrJdRF4UkQ4R3AXfROAYKXt9qhr3EzAOl9AKTr3KUOdEr46DwOvAy8BP3rKXgEpB73c0Ysd1X7bXe2067oz4FWCbt2w/0C7ofY7W3xHX481G7/W7gFeB/wIrcH2djgh6n6PVFgXq/h2QDeR69Q0Nen+j0RbAn0I+b3Z6dc0GPg9Zfm/Q+xvN94VvbRt0w8TCBPQFHgKuBFribiYpdXIDeoa8WVuFLG+MO0tUYHDQ+x2N2HH9ciruzqjKIcsrA2le2dqg9ztKbVED+CbvAyu0Pbzy+sBJQe93NNqiQN1VgXW4DuxejvXk5mdbABcBLwJ/LKSsNy7hK3Bh0PsdjfeFr20bdOPE4uRDclvlvb5PIWXtQ/54CUHvayRjB5LI/++zSSHlx4SUVw963yP9d8TdxKTAjKD3Lei2KPD6B73XXwJMLwfJLWrHNzDVqy8t6P2ORlv4+vkTdOPE4lSW5AakeK89BFQrYp28S3IxdTnO79hxZ2dHvPWbFlLe1CvLxLu5KVamCLRFFVz/igq0Dnr/gmyLAq87F3d28qz3e0wnt2gf38Agr64FQe97pNvC7/pi/saGcqi4w+mErhsrfI1dVY8A73q/jhSRynll3s+jvF/T1HvnxhC//45tcZcd01X1KxFp591A8ZSIjBSR88oacARF5D0tIknADOBHYHDpw4uqaB/frbz5Th/q8pvfbeFrffHacXIk+TacTgAiEftA4C2gH/DnkM5SfwfUBcYDw0oYZzT43Rane/ONIjIdSC1QfreIzAWuDXNgByVS7+n7cT3D/1VVy0tP+VE7vkWkCXCd9+vcstQVIX63ha/1WXLzX5mG0wmY77Gr6rci0g54Bvgz7tJDnlXA+94ZXqzxuy3qefMLcEPaPgI8Cezxlj2B+zL9J+DvJQ02wnx/X3jviVuAV1R1dhlii7aoHN8ikgjMAuoA76rqa6WtK4L8bgtf67PLkv4r03A6AfM9du9D7AvgRKA7bmyvhrjnCOsCc0Xkbr+25yO/2yLvWEvEXYYdpqqbVDVDVefh2kOB1Bh89s/XtvD6c30al8gH+lFnFEXr+H4S16FEOnBNhLdVWn63ha/1WXLzX3keTsfX2EUkGfdMWy2gi6rOU9U9qrpbVV8FuuCeZblLRFqFqysAfv8dQ9eZUrBQVVcBn+COyQ7FqC+a/G6LB4CTgCGqGovfJYUT8eNbRMYDN+D6tu2oqt+Xpp4oiNQx4kt9ltz8t8WbHxdmnVgdTmeLN/cr9q64s7RlqvptwUJV/QZYjjub6VDcIKNkizf3qy1C19lcxDp5y5sUo75o2uLN/WqLy3EPa6eKyOLQCfcPD8AAb9nUUsQbSVu8eUSObxEZC9wM/IBLbBtLWkcUbfHmfh8jvtRn37n576jhdIq4OaC4w+lEm9+xN/fm+8Ksk+HN64VZJwh+t8XqkJ/r4z68CmrgzTMLKQtSJN7TCbjnlopygjclF7O+aInY8S0iDwFDcN/DXqyqX5Y+zKjwuy18rc/O3Hymqum4D7IqwBUFy73hdFJwlxw+jm504UUg9h3evG3oYwAh9VXG3SIPRZ/NBMLvtlDV7bizVHDfpRSsry5uCCZwN9rEjAi0xfGqKoVNuEcDAIZ5y870b0/KLlLHt4iMwd01vBeX2Nb6EnAEReB94W/bBv0gYCxOFOMhblxvE18Dowsp60X+k/QnhixvhOtmSInd7rdKHHtRbeG95oD3mglA1ZCyqsAkr+xHoE7Q+x7JtvDKLiG/T8kzQ5YnAS94ZauIsQfaI9EWYbYznRh+iDtC74v7vNfsBdoGvX8Bt4Vvn52BN04sTLj/mJeFTHmddG4IXV7gNXkH4fQi6nyC/M4/X8N1+LnPW/Yysd1xcoliD9cWuOe58vrH2w7M8+rc4S3LAi4Lep+j0RZe+cPk98LwvlfHdm/ZNkL604u1ye+2KGIbea+J2eTmZ1vgRiJRb1rprVfYNDzofY7W+6Kk9RUZV9ANEwsT7mYG/bWpJH8gb52rgA9xyfIA7m64QcRgn5Jlib0Yb9azcM+5bcYlsyxgE67fvFOD3tdotoW3zuXAItx/6odwowSMBRoGva/Rboswr4np5OZXW+Ae0v7Vzx5gcdD7G833hR+fnXE9WKkxxpiKyW4oMcYYU+FYcjPGGFPhWHIzxhhT4VhyM8YYU+FYcjPGGFPhWHIzxhhT4VhyM8YYU+FYcjOmAhGRNiKiIjIt6FiMCZIlN2Mqlh7e/OVAozAmYNZDiTEViIisww011FBVs4KOx5ig2JmbMRWEiJwEnAq8aYnNxDtLbsZEkfd9mHo/Xyciq0TkgIh8LyJpItLQK0sSkZEiskFEskTkOxG5v7Bx8UL09OYvhWyvg7fNxV6d94nINyJyUES+FZERIlLJW7eZF8N2b5ufi8g1kWoLYyLJLksaE0V5iQ14CLgFWALsB9oBTYDPgPOBBUBrr7wqbtTq6sAUVe1fRN0rgDOABqq631vWAXgPN7hjDnAabrzCasAFXp1PAo/gemH/GVgBHAv8wav6GlV91ofdNyZqLLkZE0UhyW0XcKGqfuUtr4tLQCcDXwAZQDdV3eeVn4kb76sS0EJVtxaoNwX4DndJsmvI8g645AawtECdZ4TUuR54G/i3quZ45YNwg8xuUtUTfWwGYyLOLksaE4y78xIbgKruxZ1BgfverH9eEvLK1wDzAcGdxRXUwyt7qZAygNxC6lzr1ZmAO4O7NS+xeZ7CjZLeUkSah1YmIn1EZL13eXOTiAwszk4bEy2W3IwJxluFLPvGm28NTXwhNnrzYwop64G77DiviO0VVWfeNhep6uHQAlXNxg0we9Q2RaQbbqDZu4DaQD9gtIhcXcS2jYk6S27GBGNbIcsyw5SFlieFLhSRBrjvx5aq6g8l2F5ptzkYeE1V/6uqR1R1ES7ZDS6iDmOizpKbMQFQ1dwwxeHKCnMZ7nuzoi5JFqfOkmyzLbC8wLIVQBsRsc8UExPsjWhM+Xe5N38lSturA+wtsGwvkIj77s6YwFlyM6YcE5FaQEdglap+F6XN7gPqFlhWF8jGPUpgTOAsuRlTvnXDPQcXzb4kPwHOLbDsHODTX7ncakzUWHIzpnzL6yg53PdtfhsPXCIiV4hIZRG5EOjrLTcmJlhyM6acEpEkoAvwlap+Ha3tqurruNv/78f1rpIG3G69mJhYYj2UGFNS52wUAAAAYElEQVROiUh33E0kD6jqnUHHY0wsSQw6AGNMqR0ERgKzgg7EmFhjZ27GGGMqHPvOzRhjTIVjyc0YY0yFY8nNGGNMhWPJzRhjTIVjyc0YY0yFY8nNGGNMhWPJzRhjTIXzf3J4/NMbMzMvAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"N = 500\n",
"t = np.linspace(0,4.5, N)\n",
"dt = t[1]-t[0]\n",
"\n",
"y0 = 0\n",
"v0 = 0\n",
"m0 = 0.25\n",
"u = 250 \n",
"\n",
"rk2 = np.zeros([N,3])\n",
"heun = np.zeros([N,3])\n",
"\n",
"rk2[0,0] = y0 ; rk2[0,1] = v0 ; rk2[0,2] = m0\n",
"heun[0,0] = y0 ; heun[0,1] = v0 ; heun[0,2] = m0\n",
"\n",
"for i in range(N-1):\n",
" rk2[i+1] = rk2_step(rk2[i], simplerocket, dt)\n",
" heun[i+1] = heun_step(heun[i], simplerocket, dt)\n",
"\n",
"loc1 = np.where(heun[:,2]<=0.05)[0][0]\n",
"height1 = round(heun[loc1,0],2)\n",
"print(\"Height of detonation:\", height1, 'm') \n",
" \n",
"v = lambda mf: u*np.log(m0/mf)\n",
"\n",
"plt.plot(rk2[:,2]/m0,rk2[:,1],label='rk2',color='red',ls='-')\n",
"plt.plot(heun[:,2]/m0,heun[:,1],label='Heun\\'s Method',color='green',ls='--')\n",
"plt.plot(heun[:,2]/m0,v(heun[:,2]),label='Tsiolkovsky', color='black', ls=':')\n",
"plt.vlines(0.2,-100,400,lw=0.5,label='mf=0.05')\n",
"plt.xlabel('m/m\\u2080')\n",
"plt.ylabel('Velocity (m/s)')\n",
"plt.ylim(0,600)\n",
"plt.xlim(1,0)\n",
"plt.title('Velocity vs Mass Ratio')\n",
"plt.legend(loc='best',prop={'size':15});"
]
},
{
"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",
"$$m\\frac{dv}{dt} = u\\frac{dm}{dt} -mg - cv^2~~~~~~~~(1)$$\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",
"\n",
"**Answer:** Comparison: The simple rocket solution shows the velocity of the rocket exponentially increasing beacuse there is no drag or gravity effects. Once the mass ratio goes to zero, the line to turn horizontal and the rocket will maintain whatever velocity it was at once it burned all of its fuel. The rocket equation that incorporates drag and gravity shows a rocket increasing in velocity but only until a certain point where it then starts to level out. At this point the force of drag and gravity are equal to the thrust of the engine. I would expect that once the mass ratio reached zero, the velocity would start to exponetially decay since there is no thrust upward and only force acting downward. The height achieve by both rockets at a final mass of 0.05 kg is shown above the graphs. "
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"def rocket(state,dmdt=0.05, u=250,c=0.18e-3):\n",
" '''Computes the right-hand side of the differential equation\n",
" for the acceleration of a rocket, with drag, in SI units.\n",
" \n",
" Arguments\n",
" ---------- \n",
" state : array of three dependent variables [y v m]^T\n",
" dmdt : mass rate change of rocket in kilograms/s default set to 0.05 kg/s\n",
" u : speed of propellent expelled (default is 250 m/s)\n",
" c : drag constant for a rocket set to 0.18e-3 kg/m\n",
" Returns\n",
" -------\n",
" derivs: array of three derivatives [v (u/m*dmdt-g-c/mv^2) -dmdt]^T\n",
" '''\n",
" \n",
" g=9.81\n",
" dstate = np.array([state[1], u/state[2]*dmdt-g-c/state[2]*state[1]**2, -dmdt])\n",
" \n",
" return dstate"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Height of detonation: 426.36 m\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAbcAAAE0CAYAAACvs32dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOydZ5gUVdaA38MQZggDDDlKkCBJYIagIFlykCRhEXDBBGtiBcOSFFbYRXGR/UQwr4BIlBwEAQkjSQQREFAQJOccJpzvR1X39Mx0z3RP6gn3fZ56uqvuObdOV1f3qXvvueeKqmIwGAwGQ2Yim78NMBgMBoMhpTHOzWAwGAyZDuPcDAaDwZDpMM7NYDAYDJkO49wMBoPBkOkwzs1gMBgMmQ7j3NIZIvKpiKiIRIlIaR/0ConIPVt3Zgra09auU0WkeErVm1xEJNDFrt7+tiejISLPulw/FZH/eam3O45ew9S21d+4uVbq8hu9LCI/ishkEankb1sdiMhE28aD/rbFXxjnlv5w/MlkA/r6oNcbyGG//zJFLcqAiMgc+8e9yt+2ZBC6iUiehAREpAZQO43syQhkAwoAdYCXgZ9FpH9qnlBEztj39WupeZ7MgHFu6Y+NwB/2+yd80HPIngbWpqhFhszONSAP0DUROccf97XUNSdd0wLIZ2+FgPrAf4BoIBfwiYg86D/zDA6Mc0tnqJUyxtGtWMObH4rdHdLA3p2lqlGpZV96QVXvqKrY2xx/25PBmW+/enyYEhHXnoR5qW5R+uW2qt6wt0uqukNVXwZG2uXZgZf8aB8Aqvqa/duo6m9b/IVxbukT1/EPb1pvrjJejZ0YDC44urFbikgJDzItgVLAFWBpmliVsZgMRNjvG/vTEIOFcW7pEFU9BGyzd/vaT80J0c9+/UlVf3YnICI5ROQpEVkjImft4JNzIrJSRHqJiCTHZhFpbo9znRCRuyJySUTCRWS4iOT2Qj9IRJ4XkbX2uMJdETktIttE5C0ReSCOvNuAEsfgP9DLPtTGTSDAh7bsJnt/vRf2fWLL/unF9+HQ2WvrLPZCdrYte8RNWScR+cY+910RuS4iv4vIBhEZJSKVvbEnAQ4CO4AA4C8eZBwPUHOBuwlVZt9rLUTkP3awxRURiRCRCyLyvYi8lNg9ISIP2tf8kIjcEpE79r21Q0TeE5GmHvRS+1q5RVXvAkft3aIJfK4yIvKciCwVkeO2jbdE5IiIfC4ioR705tj3dTH70AQ397Xr7yDRgBIRCRCRJ0XkW/u/4J7937A8Jf4T/I6qmi0dbsAQQO3t0QTkGrvIDfMgUx74xUXO3fYNEORGt62LTHE35dmADxKp+3egcgKfIQw4nkgdP8TRCXQp6+1y/NlE6lHgQ1t2gL0fDZRPwL7cWONMCvzTh+9wuK1zDyiUgFxe4KYt+2acso+8+DzvJOH+cr1OxYHn7fd73MjmAW7Y5Y3j3BMN3ci/6oXN+4BSHmzrD0Qmor/TjV5aXKt4n9dF7rAtcywBmduJ2BeFm98xMMeLz+b6O5hoHzvowY5CwNZE6lsF5PX1eqWXzbTc0i9zsP4UIeGuSUerLQqYHbdQREKADUA14DzWeEBVIMR+HYX1JN4FmJoEO8cCz9nv1wPNgSJAZWC0XXd5YLWIBLuxrzKwDiiD9Qc/DisirxBWN1hb264rXtrzMdZg/wJ7fy0xAQCO7UW7bB6W0xIsR+eJHrYewGde2gHW9xGNFcX6eAJyXbEcKMSMtyIiHYDBLnU9ApTGenoPBfoAC4E7PtjkiTlYDqWWiNSKU9YNy8EdBbZ4UdcdYAnwJNAQuA+rNVMbeA04A1TH5bM6EJEiwIdYrcjtWPdleaAwUAPoAEwDLsTRS8trFQ8RCcT6nAD7ExA9BPwLaI11DQoDFYB2wGKsh8V3RKR5HL0BWPfgOXt/DPHv6/l4gd3zsAB4yD70IVbEZyGsB01HN3Ub4Atv6kyX+Nu7ms3zhvVjVOA6kNtNeU7gki2zwkMdM+zyy0AFDzIdiHlaqxmnzGPLDcshRdhla4Dsburu5qL/tpvyDXbZTSA0gWuRPc6+25abS7njSXdVItd4ui13FBAPMuttmY1J+A7X2rpbEpBZhfvW6f/Zx8NT4d6K1XKzjy219yfFkf3WPv6Wm3vCY0smgXOXte9pBR6KU9bTPn4XCPahzrS6Vm4/L/CGi0yHZJxril3Hag/lZ+zy1xKpx2PLDWvakMPW0R7033eRaZ3S1zQtNtNyS984gkPy4j5MuyNQMI6sExEpQEz49ihV/d3dSVR1OVYXBfg2t24AVnQYwPOqGumm7oVYjg9gkGs/vt1CaGrvjlfVXZ5O5K7uFOJj+7Uc0CxuoYiUI8ZGX1ptDhytk4dFpLyb+osBreLIOnBc25NJOG9ScNxDznFeESmJFf4OKTR/UlWPYz0wADwap9jxma9hOUBvSatrFSQiee2toIiEich7wFt2+UT795RUHC2lpiKSM3mmesTRwj0OvO1B5jWsB2JX+QyFcW7pm+XEdL+465p0HLuG1aURl0ew5t4AfO/yo4y3AXtsuTAf7HNEhe1V1V8TkHOEjhfF6q500Mrl/ec+nDfFUNUdgCMIZ6AbkYFY3ZbXSVoI/AKscRZwH6zRG6sLLhL4Ok7Zbvv1MRH5myQyyToFWApcBUpiRUeC1e2dDatVedjbiuz76gU7WOG0HRDiDH4AOtmiVeKo/mS/Fgami+fozbik1bX6DuteuI7Va7IDq6v/JtBSVV9PrAIRaSgiM0Rkn4hcE5Fol+vieMDLRUw3Z4ohIgFYXcUA33h6aFTVW8AKe/eRlLYjLTDOLR2jqhHE/OG1Epf0V/ZYWnt7d76q3o6rT+w/jj3E/CjdbY5xsyI+mOjNGANYwSxxdQAq2q/nVPW0D+dNaT6xX3uIiGNsDbuV6RiLm6uqN32tWFWvE/Pg0c+NiOPYalU9H6fsU2AvlvObClwQke/Eih5tJSI5SEFU9Q4xDtzR4nc8QHndahMrsvUXrC62VlgBK7k8iOePY8MBrDEggKeAk3bE5fsi0sPujXBHml4rNwQD/xGRwgkJicgkIBzrs1XHGivzFJWY38Px5FAYa/wUvP/dFk/FVmSqYZxb+sfRVRSANSjuoBfWmJurTFyS8uMI9EHW4QhuJCLn2r2Uz+V9sJtyfzATa4wnN9aYj4NmWN2VkLQuSQcOx1BFRJwtYzuYJiyOjBNVvQc0wRo/OYv13TTHCgL6Fjhjh7en5B+3417qKiKNsYI47hG/VekW25ZvsMbVrmF11z2CNT5bgJjgh4W2SnY31QzFGufaj/XHXwcrmnMecNYOmY8Vbp+G1+ohtZMHYP2+GmH1sADUJIGHABEZALxi767H+j1Xw3qgdFwX154Td9cmubj+/pL6u80QGOeWzlHV7VjzkCB216Tj/R/A9x7UHTdvNJBLYzJ6JLT5ktHAcfPnTUTOtfy6m/d+/eGo6kWsP2SwIvyI8/5XVfUmStATa4iJcnNtvTneX8eKLnRn21W7q6sE8CDwDPAVVvdhCJbzSMmJ+5uBY1hP947xnxX2NfKGR4npeu6iqmNUdbOq/ml/lhuqeoME7hlVjVbV6apaHaul3wdruslxrAe6AcAWuzvdVS9Nr5WqXlPVrUBnYhxcWxHxFBk71H5dj9WFOUdVD6jqBZfrktotJNffX1J/txkC49wyBo6nwToiUk1EKhITxvul2uFNbnAEkGQD4oZ3pwTH7NdqichVd6MD4JiwXNSHsZXUwtE12VhEKtrdk93tY8lptTmCYRwpwnrb4x4QE7yzwEO3smsdqqp7VXWGqvbFCnN3OMTeIpIiaZbse8kR2FLBfvUlkMSRWPmMqm5IQK6Gl/Yct53AUKwpAY40V/fjvps3za6Vy/migaexxt0Axrl8x644Uul9ncBvtmZK2uaGC8TY6e3v9rTdMs5QGOeWMfgSq/UFVoutX5wyT6zHClSA2C2SlGKz/VpLEl7uo4f9ek6t7CsOXBM8JzTPLCk4UiG5+5Nxx1piElYPxOr2zY01fzAlnvYdDqMY1vjpQ8SMOfq8RJH9lD/R5dADnmSTgOvnvQws80HXMbbm8bqLSEusoBWfsJ3IBGICdLz6zKl8rRznOIUVPg9Wy7WXa7kdferoEk3onnTrsF3w9b6OhVp5Z8Pt3S4enDAiEoQ1RQhifucZCuPcMgCqegJrtQCwnvYdP4BtcZxFXL0LxPxRPSMinTzJgjV1wA5N95b/EeM833f3QxGRLljzoiAm7N5h316seW4AI0WkTgK2+Tr+4OhG8+pP1H6SdrTQBgB/td+vSolgFzsq09G93I+YyMmTxITFx8KLFkZFl/fedhsmih0VWRnLCTzo41O7o7egiO3AYyEiBUkgWYCIVEhkXKwUEGS/d35mf12rOLxLTPfdG67TXmzH7Hh46uJOWUSeJfHIRJ/uaw84einuwwr5d8fbxEwz+igZ5/IfqTWBzmwpu2G1JuKmxxnihV4I8BsxqX2mYw28F7XLKmO1rD7DCgDoGEc/sfRbb7mUf4s1J6wQ1p/JP4hJN3QUN5Ny7fM7Ultdx8p4Ugvrh1UCKyR9MrAyjl5ik7ifsMuisSLTimAN0GcHsnm4VmXsa+R6jbul4Hc40q7zBla2GAX+nYD8D1hRrm9g/emVsK9LVazwc8d1Owbk8NGWeJO4fdD1OIkbKxrPkarrJJYTL2Pb3gv4Fav18StuJtljtbBOAe9hRQOXwwpEKY81beKgrRcBVPPDtUpw0jrwT0/3DjETqxXrwdCRFeRBLIcfRew0ee5Sm31il521fxvBLve1uDmXu0nc2YhJTBAN/BfrNxdi2/S5iw0LU+r+T+vN7wa4ufBBwAis+SNXgFtYf4zzgEYedPoCm7AGjm8AO7EGb93+iSVXz0/XxTX/oCODQ4iXumWwEjHHdY7uttZxdNMit2Q9rD/ChOrwKrdknOvlKV/lhwnYsspF7jw+/hEm8j2Us/9MXG2plYD8D158X2eBsCTYkirOzS4fQPyHBMcWgRXo4TaDDLEdgKftHjDIT9cqMecWQowj3RWnLB/WPD5P9v0EPJzItQ0lJitQ3M3klnT9jP42IM4FL09M8tGzWPOD5mLlmLsHjHSj40i7cxtrbGCRy821EAjwcK4k6fn5+nzpcuP59ESF5YR6YOWfO25/7rtYT8nrsEKU4yUPJhHn5iLXHCtc/E+73stYffvDcZM6zI1+XtuGTVhdL/ewHF44VmuuShz5BJ2bLVMGK/3YEaycgt44t54ucu+lwnf4vUv9e734PTxjX9efsYIBIrAmD2/FahkXTKIdqebcbJmmWJOAL9n3w3GsnI8N7HJPzi0EKw/ndKyHzVP2vXDdvgbvx70X0vhaJZpuDBjvIt/ezX0+Hqvl6vid7MR6oA/Eamkmdm0bYf1fOa6Nz87NlgnAGov/FutB7h7W/+4KrFay23R0GWUT+0P6HTujwB6s7qxxwDi1JjE7ygthZVY/5HKsO9af9RmgidoZFOxxo/VYYwYvqeqUOOdKkp4h8yMiHYlZr6yWelhCyGAwpG/Sk3ObgDW4+T9VHeClzk6sZvoAVf1fnLKmWMEKZ7CW1ohOrp4h8yMiC7CSPe9U1Xr+tsdgMCSNdOHc7NQup7D6gauplYInMZ3SwAmspnQBdTNPSET+xIquaqTWZMsk6xkyP3Zi40NYg/NPq2rGjBIzGAypkt4lKYRiObYTqnpARB7GynhfCKsFtUpVw+PoOMLGf3HnoGx2YDmpOsRkvU+qniETYs8/CsDqDv8I6zdxhhTKgG8wGPxDenFujln5h0Xkc+JP6B1tdxc94eKQHMuH/IFnjseRTY6eIXMymzgTbrHGW1NlUUuDwZA2pBfnFmK/NsF6in4HKzP4RfvYB1ipkK4RM7nWkfcsoUztjtyKrrkLk6rnRESexkq3Q548eUKrVk3RbD6GNKRgwYJcvnyZbNmyERgYSPHixSlYsOCcsLCwOYlrGwyGpLJr164LqurLKiQ+kV6cmyNTSnbgY1Ud7lK2REROYU0HGCAi49VadNMx+9/XQcOk6jlR1RlYIeaEhYXpzp07k1qVwWAwZElEJKHes2STXtJvuWacjjeIr6o7sRbxy0bMasneZKR3lLnLRO+rnsFgMBgyCOnFuR1zeX/Ug4zjuGPBTofOfQnUW8ZN/UnVMxgMBkMGIb04tx9d3hfyIONY4dYxHuZYVr66ncHaHfXiyCZHz2AwGAwZhHTh3FT1JFbuQ7CSgcbCziRe197daeucwHKKOYm9erJDpynWOk5niFniIcl6BoPBYMg4pAvnZvNP+3W0iDgWPEREAoFpWEu67yK2w5lgv/5LRO530SmKFWEJMNFNlpGk6hkMBoMhA5BeoiVR1aUi8g5W8txtIrINaypAfay1i04CfdQlpYqqzheRacBzwM8ishYrWapjKYhvsJZziHuuJOkZDAaDIWOQbpwbgKoOF5GtwPNY2UFyY02onozVkjrvRmeIiGzGWqqmKdY8uYPAp8A0T62vpOolhWvXrnHu3DkiIiISFzYYUpkcOXJQtGhRgoOD/W2KIQ1RVfbs2cODDz6IyzqqmZZ0kVsyI5PYPLdr165x9uxZSpUqRVBQUJa4qQzpF1Xl9u3bnDx5kmLFihkHl8X4+eef2bFjB/nz56d79+5+tUVEdqlqWGrVn65abpmRc+fOUapUKXLnzu1vUwwGRITcuXNTqlQpTp06ZZxbFqNmzZrUrFmT6OjMH06QngJKMiUREREEBXmacWAw+IegoCDTTZ6FyZYt8//1Z/5PmA4wXZGG9Ia5Jw2ZHePcDAaDIZMTHh5OZGSkv81IU4xzMxgMhkzM0aNHefjhhylbtiwjR44kqwQRGudmSBKff/45IsKNGzc8ysybN4/OnTtTqlQp8ubNS2hoKF999VUaWmkwGD76yMpFf/r0aXbv3p1luqSNczOkGpMnTyZv3ry89957LFmyhObNm9O3b1+mTp3qb9MMhixDkSJFKF7cyjc/ePBgP1uTdpipAAafiYqK4t69e4nKLV26lMKFCzv3W7RowalTp5g8eTLPP/98appoMBhsXn75Zf72t7+xdOlSOnXq5G9z0gzTcjMkysCBAwkLC+Obb76hevXqBAYGcvbs2XhykyZNIjAwkCVLlgDEcmwO6tSpw7lz51LdZoPBEEOOHDno1q0b2bNnnfZM1vmkhmRx7NgxRowYwejRoylWrBjbtm2LVT5u3DgmTJjA4sWLadOmjcd6tm7dSrVq1VLbXIPBkMUxzs0fpIcBXR8jpi5evMjatWupXdtasOHkyZPOsjfeeIOpU6eyYsUKmjVr5rGOdevWsXjxYj799NMkmWwwGLwnOjo6S0zW9kTW/eQGnyhVqpTTsbkybNgwPvjgA1avXp2gYzt27Bh9+/alS5cuDBw4MPUMNRgMXL9+nfLly/P3v/+dw4cP+9scv2Ccm8ErihUr5vb4ggULCA0NpX79+h51L126RLt27ShbtiwzZ85MLRMNBoPNF198wfHjx5k8eTJdu3bNMnPbXDHOzR+o+n/zEU9zY5YtW8ZPP/1E//793SZjvXXrFh07duTevXssX76cPHny+Hxug8HgG6tXr3a+f/bZZ7PM3DZXjHMzJIuaNWuycuVKli1bxrPPPhurLDIykp49e3L48GFWrlxJ0aJF/WSlwZC1WLx4MStWrOCxxx7LssMAJqDEkGzq16/PsmXLaNu2LcHBwbzzzjsADBkyhBUrVjBlyhQuXbrEDz/84NSpU6cOuXLl8pfJBkOmJlu2bLRr14527dr52xS/YZybIUVo0qQJCxcupEuXLuTLl48xY8awZs0aAF588cV48kePHqVcuXJpbKXBYMgqmJW4k0liK3EfOHCABx54IA0tMhi8w9ybBn+S2itxmzE3g8FgyAScOHGC/v37s3v3bn+bki4wzs1gMBgyAe+++y5ffvkldevWZfjw4f42x+8Y52YwGAwZnGvXrjmXtgFo3ry5H61JHxjnZjAYDBmc4OBgtmzZQs+ePQkNDc3SUZIOTLSkwWAwZAJq167N3LlzuXv3bpactB0X03IzGAyGTISZP2phnJvBYDAYMh3GuRkMBkMGZfbs2SxevDhLJkZODJ/G3ESkANAMqAMUAwoAl4FzwI/ARlW9ksI2GgwGgyEON27c4MUXX+TChQvUq1ePefPmcd999/nbrHRDoi03EQkQkZ4ish64ACwARgJPA48Dz9j7C4ELIvKdiPQQkYBUtNuQhowdO5bChQu7LRs4cCBhYamWZMArypUrx+eff55k/WPHjiEiiAhbtmyJVz5u3DhExOd0Yffu3WPs2LH89NNPbs+3bNmyJNucGGFhYVk2YW5WYerUqVy4cAGAc+fOUbx4cT9blL5IsOUmIn2AiUBpQLCc2w/AfuAScA0IBgoB1YCGWC27psAJEXlNVeeklvEGQ0qSN29evvrqKxo1ahTr+Ndff03evHl9ru/evXu8+eablCtXzu1CrwZDchg0aBDnz5/ngw8+YPTo0SaQJA4eW24isgWYCeQCpgB1VLWoqnZW1ddU9d+q+qH9+qqqdlLVIkBdYCoQBMyy6zEY0j2dOnVi/vz5REVFOY/9/PPPHDhwgI4dO/rRMoMhPkWLFmXy5Mn89ttv9O/f39/mpDsS6pasCAwDyqrqMFXd402FqvqTqr4ElAH+btdjyEIcP36c3r17ExISQu7cuWnTpg2//vqrs3zDhg2ICPv27Yul16xZM3r06OHcd3R5fvvtt9SqVYs8efLQuHFjfvnllwTPv2TJEkJDQ8mTJw8FCxakQYMGbNy4MVG7u3TpwvXr11m/fr3z2Jw5c2jcuDGlSpWKJ3/p0iWeeeYZihUrRmBgIA8//DDbtm1zlufLlw+AJ5980tnteezYMWf5rVu3eOaZZ8ifPz+lS5dmzJgx8RZ8/e6772jQoAGBgYEUK1aMIUOGcOPGjVgy+/bto1GjRgQGBvLAAw+wZMmSRD+rIfNQqlQpsmc3U5bjkqBzU9UpqnovKRWr6j1V/Q9QIWmmGdIbkZGR8ba4UVqXLl2icePG/Prrr3z44YfMnTuXmzdv0qpVK27fvu3zOY8fP87w4cP5xz/+wVdffcW5c+d4/PHHY5332LFjzvGl3377jR49etCiRQuWLl3KrFmz6NixI5cuXUr0XHny5KFjx4589dVXzmNz5syhT58+8WTv3r1Lq1at+Pbbb5k0aRLffPMNRYoUoVWrVpw5cwawHBPAyJEjCQ8PJzw8nBIlSjjrGDFiBHnz5mX+/Pn069ePt956i/nz5zvL9+/fT9u2bSlcuDALFizgzTffZPbs2bEeAG7fvk2bNm24ceMGs2fPZuTIkbz00kscP37cyytsMGROPLp7Vb2ZEidQ1VspUU9mY+yGsby58U2vZJ+q+xQzOs2IdezppU/z0Y8fedCIzZimYxjbbKyvJsbi4sWL5MiRw21ZaGio8/17773HzZs3+emnnwgJCQGgUaNGlCtXjk8//ZShQ4f6dN5Lly6xZcsWKlWqBEB0dDRdu3bl119/pWrVqvHkd+/eTb58+Zg0aZLzWPv27b0+X+/evRk0aBDTpk3jp59+4vjx4/To0YOJEyfGkps5cyb79u3jl19+cdrWqlUrqlSpwrvvvsukSZOoV68eABUrVqRhw4bxztWkSRPeffddAB599FFWrVrFwoULefzxxwF46623uO+++1iyZAkBAVZ8VkhICL169SI8PJyHHnqIzz77jHPnzrFt2zZKly4NWAE2jRs39vozGzIOf/75J1999RXPP/88gYGB/jYnXWPmuRm8In/+/OzYsSPeFncsau3atTz66KMEBwc7W3f58uUjNDSUhNa980S5cuWczgOgWrVqgPUjd0fNmjW5evUqAwYMYM2aNdy86dszWvv27YmKimL16tXMmTOHli1buo0UXbt2LaGhoZQvX975OQGaNm3q9eds3bp1rP1q1arF+lzbt2+na9euTscG0L17d7Jnz87mzZudMqGhoU7HBtbDRNGiRb3/0IYMw8iRIxkxYgRVq1Zl9erV/jYnXeN1R62IhADlgGOqesnleAlgAvAgcAwY6+34nCHjkD17drch/4UKFeL06dPO/QsXLvDDDz/w9ddfx5Nt2bKlz+ctUKBArP2cOXMCcOfOHbfyVapUYfHixUycOJH27duTI0cOunbtypQpUyhSpEii58uVKxePPfYYs2fPZtOmTYwfP96tnONzumvNVqzo3TCzu8/m+rlOnz5NsWLFYskEBARQqFAhZzfrmTNn3Doy49wyH3v27OF///sfAH/88YfHnhSDhS+jkK9jBZjUxZoGgIjkBDZjOT3BcnBNRaSmqp5MWVMzF2ObjU1WV+GMTjPidVWmB0JCQujcuTOjRo2KV+YIsHB0p9y7F3s499KlSx7n0/lChw4d6NChA1evXmX58uW89NJLPP/888yZ492slN69e9OxY0enY3RHSEgIYWFhTJs2LV5ZSoVklyhRgnPnzsU6FhUVxcWLF51dvsWLF+fgwYPxdOPqGTI+1apVY+rUqYwdO5aGDRvSokULf5uUrvHFuTUHjsZplfUCygMbsObDdQKGAn/DcoaGLEbLli2ZO3cu1atXJygoyK2MowvtwIED1K1bF7BWEf7111+pXLlyitmSP39++vbty8aNGwkPD/da79FHH6V79+5UrVqV/Pnzu5Vp2bIla9asoWzZsh5bSYm1MhOjQYMGLFq0iLffftvZNblw4UIiIyOdY2r16tVj1qxZ/Pnnn87rumXLFuPcMiE5cuRg6NCh9OvXL17ErCE+vji30sBPcY51BBQYrKq/A2tEpD3QDuPcsiTDhg1j5syZtGjRgueff55SpUpx9uxZNm7cSOPGjenTpw+lS5emXr16jBo1ity5cxMdHc3bb7/tbI0kh+nTpxMeHk7btm0pWbIkhw8fZt68eT7NA8qePTtz585NUKZ///58+OGHNGvWjFdeeYUKFSpw8eJFtm/fTvHixXn55ZfJmTMn5cuXZ+7cudSoUYPAwEBq1arltR0jR46kTp06PPbYYzz33HP8+eefvPrqq7Rp04aHHnoIsKYZjB8/ng4dOjB27NMDTHkAACAASURBVFhu377NqFGjUqQFbEif5M+f3+NDlyEGXwJKCmJlKHHlIeCQ7dgc7Maa42bIghQuXJgffviBqlWr8vLLL9O6dWtGjBjB1atXY/2xz549m7Jly9KvXz/eeOMNRo8eTZUqVZJ9/lq1anH+/HmGDRtG69atGT9+PE899RT/+te/kl23K4GBgaxfv55HH32UMWPG0Lp1a1588UUOHz5M/fr1nXIffvghFy5coFWrVtSrV49Tp055fY7q1auzcuVKzp07R7du3Rg5ciR9+vSJNV0gd+7crF69mjx58tC7d2/efPNN3n33XZNj0JDlEW+zSYvIFSBcVdvZ+2WAP4BPVXWwi9wsoIuq+p6vKAMSFhamCUXHHThwgAceeCANLTIYvMPcm+mfixcvMmrUKMaMGRMvuCijIyK7VDXVEtP60nI7CDS2oyYB+mJ1SX4fR640cDYFbDMYDIYszciRI5k2bRqVK1d2RkoavMMX5/YlkAfYLiJzgbeAG8Bih4CI5MKKpvzVbQ0Gg8Fg8IrffvuN6dOnA3Dt2jUzzuYjvji3acBsrHRaPYB7wFOqetVFphOWA0w8kZ/BYDAYPFKxYkVWrVpFlSpVaNeuHZ07d/a3SRkKr6MlVTUa6Ccio7AWKt2vqtfiiP0O9ATMSgAGg8GQTFq3bs3evXu5cuUKIuJvczIUHp2biPQGlqvqddfjqnoUOOpOR1V/xFqR22AwGAwpQM6cOU3GmSSQULfkbOCciCwTkUEiknjuIoPBYDAkC28j2A0Jk5BzexVr0nY7YAZwSkTWi8gLIlI2TawzGAyGLMTGjRt5+OGH+fnnn/1tSobHo3NT1Umq+hBWaP8LWCH/jYD/AEdFZKeIvC4i8dcdMRgMBoNP3Llzh2eeeYYffviBunXrMnv2bH+blKFJNFpSVU+r6v+pakusQJK/AsuBasA/gV9E5ICIjBeRVJuQZzAYDJmZH374gaNHrXCGoKAgmjZt6meLMjY+reemqpdV9XNV7QwUAXoD84CSwBvANhH5Q0TeE5GmYsJ7DAaDwSuaNWvGTz/9ROPGjZk4cSKlSpXyt0kZmiQvVqqqN1V1rqr2xnJ0nYDPgSDgReA74B8pYaTBf4hIotuGDRsSree1116LtaCmN9y5cwcR4eOPP3Yea9iwIf369fP1Y6Q47mwzGJLLAw88wMaNG3n22Wf9bUqGx5dVATyiqvewuiqXi0g2oCnQFTDrbmRwXJeKuX37Ni1atGDkyJF06NDBedyxOnZCDB06lF69eqWKjQZDZiJbtiS3OQwupIhzc8We7L3e3gwZnIYNGzrfO9aQqlixYqzj3lCmTBnKlDGLRRgMrhw9epSQkBCTWisVSNIjgogUF5G6IvKwpy2lDTWkby5evMjAgQMpUaIEgYGB3HfffQwdOtRZ7q5b8vDhw3Tq1Il8+fIRHBxM165dnQPq3nL58mXq169PWFgYly5dAiwnPGTIEIoWLUpQUBANGjRg/fqYZ61XX32V++67L958ovnz55MtWzZOnDgBwIIFC6hTpw65c+cmJCSEhx56iK1bt3q0Zffu3RQpUoTBgwdz7949ihQp4napnQYNGtC3b1+fPqch8xEREUH37t2pWbMma9eu9bc5mQ6fnJuI9BSRA8BJYAewycMWd6UAQybn+eefZ+fOnbz//vusXr2a8ePHJzgZ1dHF+fvvv/Ppp5/yySefsH//fpo1a8bVq1c96rly/vx5mjdvTvbs2Vm3bp1zsdMBAwYwa9Ysxo4dy4IFCyhatCht2rRh+/btAPTu3Zvjx4/zww8/xKpv7ty5PPzww5QpU4b9+/fTu3dv2rVrx/Lly/nyyy9p27Ytly9fdmvL9u3badGiBb169eKjjz4iZ86c9OvXj88//zyW3IEDB9i+fTtPPvmkV5/RkHn55z//ye7duzlx4gSdOnUyq6enNKrq1YYVGRkFRAOXsRYl9eTcNnlbb0bfQkNDNSH279/v9viYMWMUa8kgHTNmTLzyYcOGOcvfeeedeOVPPfWUs3z69Onxyvv06eMsnzVrVoI2esv169cV0M8++yxeWcWKFXXGjBkedV999VUtVaqUc/+9997THDly6PHjx53HfvvtNw0ICNDJkyerqurt27cV0I8++sgp06BBA/3LX/6ip06d0mrVqmnTpk31+vXrzvLdu3croHPmzHEei4yM1Pvvv187d+7sPFa5cmV98cUXnfs3btzQ3Llz69SpU1VV9csvv9SSJUt6/Dyutm3atEnz5cunr7zySiyZn3/+WQHdunWr89jw4cO1TJkyGhUV5bHutMLTvWlIG+bNm6eFChVSQN99911/m5PmADs1Ff+bfWm5vWG/vggUUdU6qvqIpy2JvtaJiLwtImpvryQg11dENonIVRG5YU8uH2oHtiRUf5L0DO6pXbs2EyZM4MMPP+TIkSOJym/fvp2GDRvGGoerUKEC9erVY/PmzQnqnjx5kqZNm1KmTBlWrlxJ3rwx6+Ju376dgIAAunXr5jwWEBBAjx49YtXbq1cv5s2bR3R0NABLly7lzp079OjRA7BW9D59+jSDBw9m7dq13Lp1y60t69evp02bNrz44otMmjQpVlmNGjWoX7++s/UWFRXFzJkzGTBggAkaMNCjRw/27dvHP/7xD1588UV/m5Pp8OUXVgnYqqpTVTUytQwCEJF6wAislkdCcv8HzALCsFqM3wKVgf8C80UkICX1DJ6ZMWMGbdu2ZfTo0VSqVImqVauycOFCj/KnT592u7JwsWLFnGNnnti7dy+HDx9mwIABBAUFxau3YMGC5MiRI169rl2KvXv35tSpU06H9/XXX9OsWTOKFy8OWM5t4cKFHDhwgDZt2lC4cGH69+8fz7ZVq1aRLVs2nnjiCbe2Dho0iK+//prbt2+zatUqzpw5w8CBAxP8fIasQ/HixRk/fjwBAeYvJ8XxtomHNc42OzWbkfZ5cgG/2OdbhOXgXnEj190uOw1UcjleDNhvl72YUnqetqR2S2ZEEuqWdBAdHa27d+/Wxx9/XAMCAvTw4cOqGr9bsk+fPtqkSZN4+g0bNtRu3bqpasLdkiNHjtScOXPqqlWrYulPnz5dAwIC9N69e7GOv/baaxoSEhLrWI0aNXTIkCF69epVDQwMdNu9q6p6+fJl/d///qchISE6YMCAWLb997//1ZYtW2qZMmX0jz/+iKd77do1zZMnj86aNUt79OihTZs29XDl0p7MdG8aMh6ko27JNUC9pLlQn3gLK7XXs0BCkQWv26+vquphx0FVPQs8Z+++5qabMal6Bi8QEWrXrs3EiROJiori0KFDbuUaNGhAeHg4J0+edB47duwYO3bsoHHjxomeZ9y4cTz33HN069aNTZs2OY/Xr1+fqKgoFi1a5DwWFRXFggUL4tXbu3dv5s+fz8KFC4mMjKR79+5uz1WgQAGeeOIJOnbsyP79+2OV5cqVi8WLF1OmTBlatWrFmTNnYpXny5ePnj17MmXKFJYsWWICSbIwJ06coEuXLhw/ftzfpmQNvPWCQFngLPAvICA1PC3QAIgEZtn7n+Om5YaVzFmBu0CQh7r+tGUeTq5eQptpuVnUr19fJ0+erKtXr9ZVq1Zply5dNDg4WM+cOaOq8Vtut27d0tKlS2uNGjV03rx5OnfuXK1ataqWLVtWr169qqoJt9xUrVbioEGDNDg4WHfu3OmU6datm+bPn1+nTZumK1as0E6dOmmOHDl0+/btsWw+fPiwAlqiRAlt27ZtrLIpU6booEGDdM6cObpx40adPn26BgcH66uvvurWtitXrmidOnW0Zs2aevHixVh1bdq0SQHNly+f3rhxw6drnppkpnszvRMZGalNmjRRQAsUKKDLli3zt0l+h1RuufmyEvdxEWkMLAa6icg62xFEe5B/29u6AUQkEPgCuIQVtJIQdezXX1T1tgeZHUApW9YxOSmpeoZEeOihh/jkk084duwYOXLkoG7duqxevdrtuBpYiWG/++47Xn75ZQYOHIiI0LJlS9577z2Cg4O9OqeIMGPGDG7evEmbNm34/vvvqVatGl988QXDhw9n1KhRXL9+nQcffJBVq1ZRr17sjof777+f0NBQdu3axYQJE2KV1a5dm5UrV/LSSy9x+fJlSpYsyd/+9jfGjh3r1pb8+fOzZs0amjZtSrt27Vi7di358uUDoHHjxhQqVIjHHnuMPHnyePXZDJmLrVu3smXLFgCuXbtmJm2nBd56QUCAKUAElkOLxpoaEHeLBqJ89bLAu1gtpl4uxz7HfcvtBfv4ogTqm2LLvJNcvYS2rNRyMySNXbt2KaCbN2/2tymxMPdm2rJ582YtW7asjh071t+mpAtILy034DXgeaxuw+XAEeCGL47UE3ZGk5eAb1T1ay9UHLHfNxOQcdiWLwX0YiEiTwNPA5Qta9ZtNbjn/PnzHDp0iNdff53Q0FAaNWrkb5MMfqRRo0bs2bMn1tQVQ+rhi3MbBNwCGqvqTyllgIgEAZ8B14Ah3qrZr76ux55UvVio6gys1ckJCwsza8Ib3LJgwQKGDBlC9erVmTVrlr/NMaQDChQo4G8Tsgy+RASWAr5PScdm8zbWHLNhqnraS53r9mtCj0COsusux5KqZzD4zLPPPkt0dDQ///wztWrV8rc5hjRm7ty53LyZUCeRITXxxbmdxGq5pTRdscbpBojIBtcNaGvLPGcfcyyedcx+vS+Beh2pL465HEuqnsFgMHjNokWL6NWrFw0bNvQ4HcaQuvji3L4GmopIaoR7OdaAi7s5Qu0q2Pth9v5u+7W63a3pjnpxZJOjZzAYDF5x9uxZZxaaffv2MXr0aP8alEXxxbmNAw4BS0SkYkoZoKrlVFXcbVhTAwCG28dq2zongB+BnEDPuHWKSFOsOW1nAOdqm0nVMxgMBm8pVqwYkydPJleuXJQrV45p06b526QsiS8BJUuwIiWbAwdE5Hc8z3NTVW2TAvYlxARgHvAvEdmqqkcARKQo8IEtM1GtxVNTQs9gMBi8YtCgQdStWxeAggUL+tmarIkvzq1VHL3K9uaOVI8gVNX5IjINK2XWzyKyFmsOXksgGPgGKxFyiugZDAaDL9SpUydxIUOq4YtzezTVrEgiqjpERDYDQ7HG5AKAg8CnwDRPra+k6hkMBoM7fvzxR2rXrm2WMkpHeP1NqOo6X7aUME5VB9pjbe8kIDNbVRuparCq5lHVUFX9v8QcVFL1DN5z+vRp2rdvT/78+RERNmzYkOS6PvroIypVqkRgYCChoaGsW+fdLbZlyxYaNGhAUFAQ5cuX5/33348nIyLxtoYNGybZVkPWYtu2bTz88MN07dqV69fNDKL0gnnMMKQa//znP9mzZw9fffUV4eHhzjEIX5kzZw7PPvss/fv3Z+XKlVSvXp2OHTuyb9++BPWOHDlCmzZtKF++PMuXL+eZZ55h2LBhfPzxx/Fk//73vxMeHu7cPvnkkyTZashanD9/nscee4y7d++yZMkSBg0a5G+TDA5SM7dXVthMbknPtGzZUrt27ZrseipXrqxPPvmkcz8qKkpr1KjhXB3AE08//bRWqlRJIyIinMeee+45LV26tEZHRzuPATp16tRk25nRyMr3ZkoRFRWlr732mgIaEhKiR44c8bdJGQb8tZ6biHxv53xMMiLSSES+T04dBv8zcOBAwsLCWL58OdWqVSN37tx06NCBS5cuceTIEZo3b06ePHkICwtj7969gNXVt27dOhYtWoSIUK5cuSSd+/fff+fQoUM8/vjjzmPZsmWjZ8+erFy5MkHdlStX0q1bN7Jnjxla7t27N3/++WeirT6DwRuyZcvGhAkTmDZtGvPmzaNixRSbJWVIJgl1S1YFNonItyLSS0RyeVOhiOQSkT52FOL3eI6oNGQgjh8/zujRoxk/fjwzZsxg69atPP300/Tu3du56GdkZCS9e/dGVQkPD6dOnTo0b96c8PBw5+KhqkpkZGSim4ODBw8CULVq1Vj2PPDAA1y6dInz58+7tffmzZucOHHCrZ5rvQ7Gjh1L9uzZKVy4MH/961+5dOlS8i6YIUvx7LPP0qJFC3+bYXAhoWjJSsCbWMmMWwDXRWQL1uTmA8BFrGTHwUAhrNWzHwIaYeVnjMRaPmZsKtluSEMuXbpEeHi488l07969TJo0iS+++IL+/fsDluPq0KEDBw8epGHDhgQHBxMSEhIrOOOLL77wajVqq9cCLl++DMRPOOuYO3T58mWKFCkST//KlSuJ6jkYMGAAnTp1okiRIuzcuZNx48axZ88etm/fTkBAQKK2GrIWu3btIjQ01N9mGBLBo3NT1avASyLyPtZSNwOAdsTke3SHYC02+i7wgaoeSzlTMxcbNmxIVvRgUmnWrBnNmjXzWa9cuXKxulzuv/9+gFhPq45jJ0+edLaQ4tKpUyd27Njh8/lFJNa+w/nFPZ6Ynrvjn3/+ufN9kyZNeOCBB2jfvj1Lly7lscce89lWQ+bl008/ZdCgQbz++uuMHz/ehP6nYxKd56aqvwMvi8gbWHPCmgG1gaJAfuAKcA4rrdV6YJOq3k0tgzMLSXUy/iJuCyhnzpzxjjuO3blzx2M9ISEhPq1C7GhpXblyJZaep5ZZXHsdcg48tQRdadu2LXnz5uXHH380zs3gZNOmTTz99NMATJgwgUKFCvH3v//dz1YZPOH1JG5VvQ2ssjeDIUn42i3pGDM7ePAg990Xs5jDwYMHCQkJcdslCZAnTx7KlCkTb2zN0xieK45WXWKtQkPWIiwszNmir127ttPRGdInvmQoMRiSja/dkhUqVKBy5crMmzePNm2sdKXR0dHMmzePdu3aJajbrl07Fi1axPjx451jZ19//TVlypShRo0aHvVWrVrFjRs3zLiKIRZBQUEsWLCAf/zjH7z00kvky5fP3yYZEsA4N0OaUqhQIQoVKuSTztixY+nXrx/lypWjUaNGfPHFFxw+fJjZs2c7ZTZu3EjLli1Zt24dTZs2BWD48OHMmjWLJ554gqeeeoodO3Ywffp0pk2b5myVzZgxg507d9KqVSsKFy7Mjz/+yPjx46lfvz4dOnRIuQ9uyBTkyJGDf//73/42w+AFxrkZ0j19+vThxo0b/Otf/2LcuHFUr16dZcuWxWp9qSpRUVHO7kywAlxWrVrFsGHDaNeuHcWLF+fdd99l8ODBTpmKFSvyxRdfsGDBAq5du0bx4sXp378/48aNM5GSWZwTJ05w6NAhWrZs6W9TDElAXP8MDL4TFhamO3fu9Fh+4MABj5GDBoM/MfemZ86cOUOzZs34/fffmTVrFj17xlv+0ZBMRGSXqoYlLpk0TByrwWAwxOGZZ57h119/JSIigv79+3P69Gl/m2TwEePcDAaDIQ7Tpk2jcuXKBAQEMHPmTEqUKOFvkww+YsbcDAaDIQ4lS5Zkw4YN7Ny5k06dOvnbHEMSMC03g8FgcEOJEiWMY8vAeO3cRGSwiASlpjGZFRO0Y0hvmHsyhsuXL9OpU6d4E/4NGRtfWm4zgD9F5F0RqZRaBmU2cuTIwe3bt/1thsEQi9u3b5MjRw5/m+F3zp07R/PmzVm2bBnNmjVj//79/jbJkEL44tyWYa0A8DJwQERWiUgnMTmKEqRo0aKcPHmSW7dumadlg99RVW7dusXJkycpWrSov83xOydOnOD3338H4OzZsyQ0rceQsfAlt2RnESkDPAf8FWgNPAqcEJEPgU9U1f3iWlmY4OBgAE6dOkVERISfrTEYrN6EYsWKOe/NrExoaCjLly+nQ4cOTJ061bl8kyHjk6RJ3CKSA3gca623hwAF7gHzsZa6CU9JI9MziU3iNhgM6Z8LFy5QuHBhf5uRpUiXk7hVNUJVZ6lqI6AO8AnW4qR9gc0isktE/urt6t0Gg8GQFhw4cCDeMkiAcWyZkGRPBVDVPVgrdn+GtVipYDm8j4BjIjIouecwGAyG5LJlyxYaNWpE+/btuXHjhr/NMaQyyXJuItJKRBYCR4GhwB3gU6APsAJrQdMZIvJCcg01GAyGpHLu3DnatGnD5cuXCQ8Pp1+/fv42yZDK+OzcRCS/iLwkIgeB1cBjwCngDaC0qg5W1a9VtRPwMHATMM7NYDD4jaJFizJp0iTn+1GjRvnZIkNq43W0pIjUxQog6Q0EYXU/bgCmAotVNTqujqpuE5HlQPcUsdZgMBiSyHPPPUdERAQdOnSgYsWK/jbHkMr4klvSERJ4C/gYmKqq+7zQu+njeQwGgyFZREVFER0dHW+i+gsvmE6krIIv3ZLHgOFYXY/PeOnYAJ4CTCoEg8GQJty6dYuePXsyePBgkzghC+NLi6qiJuFOsXWifNUzGAwGX7l58ybNmjVzZhopUaIEEydO9LNVBn/gS8tttYgMS0xIRF4WkTXJsMlgMBiSRO7cualfv75zPyIiwrTesii+tNxaAX96IVcNaJk0cwwGgyHpiAhTpkzh+PHjtG/fnueee87fJhn8RGoEeuQE4kVOGgwGQ0rjaJW55m/Pnj07S5YsweR0z9qk6GKl9goBocCFlKzXYDAY4hIREcHQoUPdjqkZx2ZIsOXmZuysdQLjadmBSkBJrATKBoPBkCpcuXKFzp07s2nTJgDuv/9+evbs6WerDOmJxLolW7m8VyzHVTIRnb3AiOQYZTAYDAmRN29ecuWKycu+YsUK49wMsUjMuT1qvwqwBivd1jseZO8BJ1X19xSyzWAwGNySPXt25syZQ8OGDRk8eDAjRpjnaUNsEnRuqrrO8V5EtgAbXY8ZDAZDWqCq8cbRChUqxN69ewkKCvKTVYb0jNcBJar6iKqa2ZAGgyFNOX/+PK1atWL+/PhD+caxGTxhcj4aDIZ0y759+2jbti0nT55k27ZtVKlShZo1a/rbLEMGwKNzE5E37LfTVPWyy75XqOrbybLMYDBkee677z7y5MkDWDkjN2/ebJybwSvEU2oaEYnGipB8QFUPuewnWidWSsmAlDMz/RIWFqaOPHYGgyHl2bt3Lx06dGDGjBm0a9fO3+YYUggR2aWqYalVf0Ldkm9jObMLcfYNBoMhVbh8+TIFCxaMdaxWrVr89ttv5MyZ009WGTIiHp2bqo5MaN9gMBhSClVl6tSpjBw5kg0bNlC3bt1Y5caxJY+oiHvc/GU3wbv3w7ZtoArTp/vbrFTFBJQYDAa/M3LkSN5+2xqm79mzJz/++CP58+f3s1UZl7NH97F145ds/+17tl0/yM7cV3j8F/h4iS2QNy988AEEZN7RI+PcDAaD3xk0aBBTp07l+vXrFChQgGvXrhnn5iUaHc2R3evYvPUrNh3fzOboYxwOjrAKswN2L+/2Ui5KN27AgQNQo0Zam5tmeO3cROQ5YArQVVWXe5DpCCwEhqjqxyljosFgyOxUqFCBDz74gF27djFhwgQCAwP9bVL6JSoKfvoJNm1i84/f0KPEJs7mthdiyetZ7VKebNzr0p6c9R+C+vWhQoW0sddPeIyWjCco8i1QEyipqm6XtBGRAOAk8JOqtk0xK9MxJlrSYPCNgwcPcv78eR555BF/m5Ih0Oho9v+wlI2bZ/Lslrtk2/g9XL0KwJ/BUMbNEtI5I6He9Xw8lKcKDSo0of5DPShTtQGSLUUXgkkW/oyWjEtV4GdPjg1AVaNE5GesBUsNBoPBiary0Ucf8dJLLxEcHMzevXspWrSov81Klxzbt5l1333Muj828F32E86W2cPbofbVGLnS16DcZbgcBI3uFOGRQnVpXLszYc36Epi3gJ+sTx/44tyKABu9kDsHmEcyg8EQixs3bvDmm29y+/Ztbt++zZAhQ9ym1MqKXDv/J2uXvc/qA0tZG32E3/NFWgXBseXWlYfaZ4CSJaFJE3jkETaFVqJkWHOyBZgQCld8uRpXgTJeyJUCbiTNHIPBkFnJly8fn3zyCe3ataN69eqMGTPG3yb5D1XYuxdWrmTI71P5qPgpIgOAPO7FC94Rmt8tQdWeXeDTF6BKFbATSZdOO6szFL44t91AcxGpqKq/uRMQkYrAw8D3KWGcwWDIuERGRpI9e+y/mLZt2zJnzhw6d+6c5ZIeXzn7B6e/XcgD6/fBqlVw6hQAhZpDZKnYsrnvwSO3CtGySANaNOhD7SaPE5DDzPXzBV+c2+dAa+AbEemmqoddC0XkfmAREGDLGgyGLMqiRYsYNmwYa9asoVKlSrHKevXq5Ser0haNjuaX8MUsXf8hKy78QHj+a4Segm2fxpZrdwTGN4W6V4Jom68ObcJ60/DRJ8kZlEDooyFRfHFuXwP9gPbALyKyGThol1XBGmfLDqxS1ZkpaqXBYMgwjBkzhrfeeguAgQMH8v333xOQiScLuxJx9zabV0xjyfaZLIn4OWbszJ5rtqMUnM8NRW4BBQtC69bUb9ua003CKF6hlt/szox4PRUAQERyAu8BTxHfMUYCHwHDVPVuilmYzjFTAQyG2Pz44480aNCAyMhISpYsyffff0/FihX9bVaqEXnlEgvnj2PJgW9YkfMPLgd6/k+tczGQzwv0p1a7gVCvHmTPukEg6WkqAKp6DxgqIuOAlsB9dtEfwDpVPZPC9hkMhgxG3bp1GTVqFEeOHGHKlCnxEiFnCk6cgKVLYfFiZMN3PP9CJOeC44vluwtt75Sm0/0daNP+eT74fB61xo5Nc3OzIkl6bLCd2KyUMkJEcgBNsLo8G2E5zULAeSAc+K+qbkhAvy/wHFALa8zvIPAZ1lp0HuflJVXPYDBAdHQ0M2bMoFq1ajRp0iRW2ahRoxA7mi+z8Ov2lSxYPZmwLcdovfqI83gA0PEQfGrnei59I4DO2R6gS2hfmnYYSq48rl5vXpranJVJL23ipsC39vszwC7gJtZk8O5AdxEZp6qj4yqKyP8BQ4A7wDogAqtV+V+gpYj0VNWolNIzGAxw9OhRBgwYwKZNm6hUqRJ79uyJFf2YGRybRkezb8si5q97nwXX5vNtXQAAIABJREFUtvFLfmu0pXt+K7LOlf637qeMlqJz06ep07R3usoEklXx2bmJSBXgBaAZ1pw2sFJurcdqYR30oJoQ0cACYIqqbopzvl5YrcRRIrJeVde7lHXHclBngCaOCE4RKWbb0xX4G1ZOTJKrZzAYLHLlysWePXsAOHz4MFOnTmXEiBF+tir5aHQ0u76byYIN01hwe5eVgFgAlxzOKyrBzdzZydO4BXTuDJ060bRsWZr6zWqDO3xybiIyEJgG5MT6yh0EAw8Ag0XkGVX9wpd6VfU74DsPZV+LyKPAIKxozfUuxa/br6+6Tk1Q1bN2oucNwGsiMjVON2NS9QwGA1CyZEkmTZrE0KFDGTFiBC+88IK/TUo60dGc/34lE9eMZsG9PfyRLwpyYG0uBEVA21sl6VH5MQKOj4ZCxfxirsE7fFkVoB5WNGQ2rPlsnwK/YTm58sBfgW7ARyKyX1V3pKCdu+1X52R8ESkNhAL3cNORraobReQkVuuyIbA1OXoGQ1blzJkz/PLLL7Rs2TLW8cGDB9O0aVOqVKniJ8uSQVQUbNoE8+fDokXkuniK/w6He/lii+W9Bx3vlKV7tR606/YqeQqaXJgZBV9absOxHFs/Vf0qTtlBYKWI9MHqQnwFSMmZmo5ZoKddjtWxX39R1dse9HZgOak6xDippOoZDFmKyMhIPvjgA0aNGkVAQAAHDx6Mleg4W7ZsGcqxaXQ0P66fzex1/6HP8j8I23vBWRYMPPo7LK8MBe4InSPK0/3B3rTuOjzLJyDOqPji3BoDu9w4Nieq+pWIvIQV+ZgiiEhxYKC9u8ClqLz9+kcC6sfjyCZHz2DIUkRERPD+++9z7do1AF599VU+++wzP1vlO4d2rGL28ol8dX0rh4IjIBfcLQ1he12EChfm9QINeb5KQ5p3edFkB8kE+OLcCuFhXCwOh4HaSTMnNiKSHZiJNZy7TlWXuhQ77r6bCVThSODs2tmQVD1Xu54GngYoW7ZsAtUYDBmXoKAgpk6dSvv27alcuTJ/+ctf/G2S15w8tJOvF7zF7HPr2FXgljV44hKRP7c6/GdnYbJ37wk9ekCTJjTKwhOqMyO+fJuXAW/SDFSwZVOCD7HC809gBZO44gho8T7FSvL0nKjqDGAGWBlKklqPwZBeuHfvHps3b6ZFixaxjrdr1465c+fSuXNncuXK5SfrvOPamWPMnTOa2ceWsqHAFVSAOD2Kee9Bt7sV6BP2JDLuFchlVvzOrPji3LYCXUSki6oudicgIp2wgjAWJdcwEZmCFSF5BmjpJvvJdfs1of4DR9l1l2NJ1TMYMiUrVqzg5Zdf5rfffmPPnj1Ur149VnnPnj39ZJkX3LxpZQqZPZtz21by1JBIZx5HBzkjocONEvSp0ZuOPUcSFBziH1sNaYovzm0y0AWYJyIzgS+Ao1gtoApAf6zWVbQtm2RE5F2suXTnsRzbYTdix+zX+9yUOXCsP3fM5VhS9QyGTIeqMnnyZA4dOgTAyy+/zOrVq9P1JOyIO7f4dsG/qfXtXkrPX2M5OOB+oP6fsL00ZIuG5lcL0rdCF7r1GkuBYgn93A2ZEa+dm6putoNFJgMD7M0VAaKAl1R1S1INEpF/A8OAi8Cjqrrfg6hjekB1EQnyEPlYL45scvQMhkyHiPDee+9Ru3Zt8ubNS5s2bVDVdOfcoqMi2bLiQ2ZvnsY8OcDFIGX87/CPOCPnr5y7n1OVavJ499GUqJgiQ/+GDIqviZOn2kvdOCIiS2I5tZPARqwMI0l2CCIyEWvKwWUsx7YnAVtOiMiPQF2gJ/C/OHU1xZoXdwYrP2Wy9AyGjE5ERARLliyhW7dusZxXzZo1mTlzJi1btowV6u9vNDqavZvmMXvNZL66t4sTeaMgd0z57JrwxiaQqlWhb1/o04ee99/vP4MN6Qqfw4Ns5xW31ZZs7JUGXgWuYDk2b5zkBKyJ2P8Ska2qesSuqyjwgS0z0U2WkaTqGQwZkm+++YYRI0Zw+PBhli9fTvv/b+++46uo0gaO/54UIKF3qQoG2YB0BQlGQomggBiaUlZkF1HBRV0Rfd+1oC4IyO67FgQXYVkXiZSAlFAEBBYFKQIKSO8gXQkBAiHkvH/MJLkJSUyZW5I8389nPpPMOXPmzMm998nMPXPOww+nS+/bt6+Xanarg9u/JnrhGKIvruOnsonWeEgZJqGuedmfLuWac+P7iRRrdg/42JWm8j6f6PsqIo8Ar9m/HgD+lMVtkT3GmLEpvxhj5orIJKyR/XeIyErSBkAuA3yJNRByOnndT6mCasWKFezfb311PXLkSDp16uRbE4ieOgWzZ7Nt8RSa37/L2lY2fZYKCUIfE0rfNs9wf5dn8fP3iY8v5aN85dXh2n3pHnvJzFpgrOsGY8xQ+1bpMKzZBVKmrplGNlPX5HU/pQqiN954g88++wx/f38GDhzIzZs3vR7cLp45SrFFSwj+IgZWr4bkZJoI1L4bjtld+EsmwqPX7qBv84FE9hihD1erHMsyuInIP/NRrjHGPJ2LzNOB6fk42Exgpqf2U8pXnTp1io8++ohRo0YRGJg28m/VqlWZN28eLVq0oEIF73WFT7j0C7FzRjNzZzSxpU4xKRb+4PIFhJ+BgTv82N6oMv1C+9Ct92s6nqPKk+yu3Abno1wD5Di4KaXy7+9//zuvv/46V69epUaNGgwdOjRdemRkpFfqlZR4jVXz/8bMzdOYX+wQ8cVJfbh6ZiM7uIlA27bQrx9v9+wJXgzAqnDILrg95bFaKKXyTUS4evUqAKNGjWLgwIGULFnSK3UxyclsWDaFmf+dyGyzk3PBJtPB7OIqBJP43hsUe7w/1Kx5awal8ijL4GaMmerJiiilci6zZ9GeeeYZJkyYQKVKlRg/fjzBwcFZ7O1GO3eyK/p9ul6fzpHSSRB0a5Z6lwLpVzqMvl1eof69D3m+jqpI8JUOJUqpHEhMTGTatGl88MEHrFu3jooVK6amBQUFsXbtWurUqePZziJHjkB0NMycCTt3UjcALrycPkv1K348FtCUfu2fp0X7AYifn+fqp4qkPAU3ESmF1aOxMnDMGLPR0VoppTIVFRXFkiVLABg/fjzjxo1Llx7ioYeYzx7Zxew5o4g+uZzXYuN56EBaWlASRO2Ghb8TeiXdRb+wITzQ9Tn8A4tlXaBSDsvVv08iUtruRXkeWAV8gUvHERF5VkSOiUhLZ6uplAIYNGhQ6s9z587lxo0bHjv2pXMn+GziEDq/WInq/7qbP12dy/ry8cxs5JIpKAgee4wJAz7j9GsXmfLeHtpF/VkDm/K4HF+5iUgwsAZrdurzwFbgwQzZvgImAlHAJmeqqFTRc+PGDdatW3fLFDQ9evQgMjKSyMhIhg0blq67vztcu3yRpXPfZeaPn7M4+CTXArllGpnYu+BGl84E9h0AjzwCpUtT2a21Uuq35ea25EtYgS0aGGKMuSIi6R50NsYcFJH9QPvMClBKZc8Yw5QpUxgzZgzHjh1j165dhIaGpqb7+fnx1VdfubcSN29ybNks3lozipiA/cSV4JbRQgDu/7UM/Wo9TK9BbxBYO/TWDEp5UW6CWx/gFPBHY8y1bPIdBRrkq1ZKFVEiwsKFCzl69CgAo0ePZsaMGe4/sDGwebPVKWTWLIrHn2b6S5Cc4YuLJhdL0K9SOx6Pep3aDVq7v15K5VFugtudwPLfCGxg3bKslPcqKVW0vf7668TGxlKpUiUaN27s1mPt3riYmUvG8VjsEe7+/kTq9qpAx0PwVQjUjQ+gX8n76PvQyzS47xG31kcpp+QmuN0AcjLPfE3gct6qo1TRcOHCBT766CM2bdrE4sWL0z2z1qpVK2bNmkWXLl3c8hD28d0b+WL+28w8t5rt5RLAD65XgPGumapW5e3KbXn7vgdpGTlIu+6rAic3wW0f0ExEihtjrmeWQUTKAU3QST6VylJCQgL169fnwoULAKxZs4Z27dqly9OnTx9Hj3nu2G7mzHmLL44vZV35S9ZGl44h0XfD2I2l8evZy5obLSKCVgH6GKwquHLz71gM1t2KMdnk+StQCmuuNKVUJoKCgujVq1fq79OnT3fLceLPn+TfE5+i84uVqDa1AcMuz0oLbLbiSdAzrgbvN3oZc+oUTJsGHTuCBjZVwOXmFfwh1iSlL4jIPVjBDuB2EXkKa1brDsAu4FNHa6lUAfXtt99y6dIlHnoo/TBTI0aMYNOmTYwYMcLZq7SEBIiNhehojm1cxJNP3bil675fMnSMq0i/elE82ucNylau5dzxlfIROQ5udtf/B7GCWjhwv50UYS8CbAe6Z3XbUqmi4ujRo/Tv359vv/2WOnXqsG/fPgJcroZCQkLYunWrI8e6ce0qq+b/jaZf/chtMcshPh6AhkCjM7CjqpUv7GJp+tXoTO9eb1LljoaOHFspX5Wrew/GmONASxHpCjwM1MWa5PM4sBSI0Uk+lYIqVaqwd+9eAA4fPsycOXPo27evY+Un30zi2yWTif5mEnNkN+eDDBP2wkvx6fO9dKI2p2s05PFH/8LtDds4dnylfF2ebqwbYxYDix2ui1IF0vHjxwkKCqJSpbQnYIKCgnjmmWcYN24cAwYMoEWLFvk+jklOZtuaaKJX/YMvErdxotRNcBn4P7oRvLQBqFcP+vaFxx9nYKg+XK2Kpuxm4p4LTAWWGWOM56qkVMGwY8cORo8ezdy5cxk5ciRjxqTva/XCCy/w9NNPUzOf85Tt3bSU6CVjiY7fwL4yN6AY1uKixmU/Iio2I2nTRALuaWlN/qlUEZbdlVsPrDEiT4nIv4Hpxpj9nqmWUr5v//79zJo1C4BPPvmE1157Ld0caq7T0eTa0aMwezYblk4hrO1+6xvtMumzVEwQeptQ+rZ5hvu7PIufv/ZwVCpFdu+GScBjQHXgVeBVEfkGmAbMMcZc9UD9lPIJly5dokyZ9NGle/fu1KlTh8OHD9O4cWPOnDlDnTp18nyM4z99R4Wlqyk5ZwFstGaRailwW3M4bc9iXSoRoq7XpW+LgXSMGkFgCS9MSKpUAZDdTNzDRORFoDvwB6Ajab0kPxSRWcC/jDHrPVJTpbxg27ZtfPTRR0RHR7Nx40YaNUqb38Xf35/JkydTpUoVmjZtmqfyT+7dzJwvxzD71Eo2lL/Mf+bBgB/T0v0NPLHTnwOhVel79+N06f0XgspUyO9pKVXoZXsfwxiTiPVA9hwRqY71nNsTQH3gj8AfRGQf1tXcf4wxp91cX6U86t1332XOHGtMgo8//phJkyalS3/wwYyzPv22n/dvJebLMcw+uYJvUh6qLm+tZje0g1tAAERGQp8+jIuKgrKZDMuvlMpSjkcoMcb8bIx51xgTCrTB6mwSjxXoxgLHRGSBiHQXEQ/Oca+UMzKb+PO5555L/fmnn34ir32rTh/6kYkT+tD2hXLU/LwFw6/GpAU2m38ymMqVMZ9+CmfOwJIl8OSTGtiUyoO8PgqwAdggIsOBXsCTWA9yd7WXc8BtzlRRKfcxxvD1118zceJEjh8/zqZNm9INYhweHs7IkSOJioqiVatW6dJ+09mzEBPDyhWTebDxjxgh9QothV8ytI+rQJ86XYjq8Rcq1arvzIkpVcTlq3uVMSYB+A/wHxGJBGYAle1FKZ8XFxdHt27dSEhIAGDDhg2EhYWlposI48aNy3F5547tpvzSNQTMiYHVqyE5mfuKQfEGWLNYYwW0iLjy9Ln9IXr0fI3KOtGnUo7LV3ATkVJYPSqfBMKwOiyDNWKJUj7FGENiYiLFi6fN3FSuXDn69+/Pp59aw6GuXr06XXDLiZ/3b2X+grHEnFjB2rIXWT7DmgstRalE6LYPzlUvR59anenR4y9UrXO3I+eklMpcnoKbiLQDBmE9CxeEFdSuAwuxOpd85VQFlcqvc+fOMX36dKZMmcLgwYMZOXJkuvTnnnuOEiVKMHToUEJzOKLHkZ3fELNoHPPOrGV9eXvMK5dOIR0PYT1IHR4OffoQHfUo/tVrOHhWSqns5Di4iUgdrN6SA4HapF2lbQf+BcwwxvzqeA2VyqelS5emBrQpU6YwYsQI/Fwm32zSpAkffvjhb5azZ/NS5i37P2J++Yat5azbmBm/QxMDZ2pVgPffhJ49oYYV0LSHlVKelW1wE5FgrKlsnsR6xk3s5VdgJjDVGLPdzXVUKsfi4uIom6F3Ya9evRg+fDhxcXGcPXuWgwcPUq9evd8uzBj48UeIieGz76cxsOVJa3uGKWT8k6FtXHl61nqQqEdfpdqdeXvmTSnlnOzGlpyKFdhKYgW0ZGAl1m3H+fYzcEr5hAULFjB16lSWL1/O/v37qV27dmpacHAwb7/9NmXKlKF3796ULFkyy3JMcjLbVs2g2cqdSMw8OHgQgA6lgZZp+QJvQmR8ZXrW6cIjUa9qL0elfEx2V26D7PVhYDrWaCQn3F4jpfJg0qRJLF++HIBp06YxatSodOnDhw/Pct+kxGt8s2Qy8zdOZ17STk6UuskP/4bGZ9Ly1IiH9kf8KFfuNnrU707XHq9StkrtLMtUSnlXdsHtc2CaMWa1pyqj1G85f/48Z8+epUGDBum2Dxo0KDW47du37zfLufLrWZZ/OYEFu2JYHHCYX4IMlEhLjwm1g1vp0tC1K/TsycpOnZBSpZw8HaWUm2Q3tuTvPVkRpbJz4MABXnnlFRYtWkRYWBhr1qxJl969e3fefPNNBgwYQEhISKZlnDm8k0UL32PBkWWsKHWW6wFA6Vvzlb8mSNPG8NJfoWNHKGFFPZ1ERqmCQ+fIUAVCyZIlWbBgATdv3mTt2rUcOnSIunXrpqaXKFHilluRAOzbB19+CQsW8HHx9bzdlls6hIA1H9ojfqFE3TOAiG7DdbR9pQo4DW7Kpxw/fpzo6GiGDBlCuXJpUahatWp07tyZ2NhYWrVqxYULF9IFtxTJN5PY+NU0tnwzhz/NOwF79qSmda+GFdxsd8cV59FS99D9gSG0aD8A8cvxUKtKKR+nwU35jGHDhjFp0iSMMVSoUIHBgwenSx8zZgzvvffeLQ9aX7t8ka8Xvs+X279gEfs4XTIZikHPk9ZkhCmanRF6na1IWK0wund+gbpN23ngrJRS3qDBTfmMkJCQ1FH3P//881uCW+PGjVN/PrlvC7FL/sHioytZFXyGq8WwHlpxseguePqnIOjUCbp3R7p2ZU6lSu4+DaWUD9DgpjwmOTmZ9evXM2PGDAAmT56cLv3xxx/nlVdeoV27dvz+97/PuDNs3swHsW8wPX4d21JGCMnk+7PKV4VuySE0fmkwdH8OgvX7M6WKGg1uymP27NlDeHg4YHUAGTduXLrRRKpVq8aZM2coX94a08rExSErVsDixbB0KZw9y+4usO3eW8uuf6kY3Uo0pnvrQbTuNBj/wGIeOSellG/S4KbcYufOndSrVy/dCPwNGjSgadOmbN++nWvXrrF48WL69++fbr/zh7YwfcVEFp/+L+VOXyRmVvrJQbvug8n3WiOEtL1Uga7V29IlchghzTp45LyUUgWDBjflqE8++YQPP/yQXbt2sWDBAh555JF06c8++yxbt25lwIABhIWFcf3KJdYtm8yS72cRe30H+8rYs2GXh6BSkBAAQUn2zlWq0D68EzG1b6Nj1+GUqVzTsyenlCowNLgpRx05coRdu3YBMGvWrFuC25AhQziwbRXLVr7P2JierA4+a3UGKW4vLhICYVO7erRt3dcaJaRFC4L8/OjhmVNRShVgGtxUrh0/fpzZs2cD8NJLL6VLe+yxxxg7dixBQUEEBQVZG69csWalXraMxK+W0rT3Ia4UI9POIMGJEHn1Nrre8SAPP/w81d9s7uazUUoVRhrcVK788MMPNG1qTelSuXJlnn/+eQIC0l5GTZo0YV5MDDVKX+O/W2ayrkcLwmN3QqI1iUQxoP1hWOQyiH69S4F0DgylS7M+tO0yjBKlMol6SimVCxrcVKaMMezevZvf/e536Sb2bNSoEdWrV+fnn3/m3LlzrFmzho4dO3LxzFFWLv6AZXsWs4wDnCyVDMDgJAjPMDlS1KFiULU8nWtG0KnDEO5s2t6Tp6aUKgI0uKlbTJgwgU8//ZS9e/eyfv16WrdunZrm5+dHnz592LnjR+5pWI2VK0fx5uIebCwbz00/IMOg+ctCwADSqBF07gydOzOoTRsGFc/wBZtSSjlIg5u6xZ49e9i7dy8AMTExVnAzBnbvhhUr6H7qO6a2/I6VKfEp6NYyyl6HyGs16XxHR24eG0VArds9dwJKqSJPg1sRdOPGDdauXcv8+fMJCQnhxRdfTJfeo0cPpk6dSnBwEDu3fkXSEwMIWPk1nDoFQGhJiH85fZli4J64knQu05xOrfrRquOTBBQrgVJKeYMGtyJo2bJlqV30Q0NDU4Nb/IWfWbt0Mst+WEDNbgGcaJzA8sAdbPl0B/edStu/6hVofBoulvInUu4ksl5n2nd6hsq1QzM7nFJKeZwGt0Ls2LFjrF279pZxGjt06EBQUBAJCQns3r2b555ryQ8Be/muzCWS/LG+N2uRln9lXbjvBFCuHLRvD5GRrAm/l3KhzXSaGKWUT9LgVgglJydz//33s2HDBgDCwsK48847re74mzcTvHo1v7vdnx9Lwc2GMLHCZvC/tZxiSdAmvhx1Ix+ECSOgeXPwtzKW9+QJKaVULmlwK+ASExNJSkoi2GXkez8/PypWrJj6+/g/9eWTpHLw7bdw9SoAj7WBbZG3ltf0YhAdgxsS2aQH93d6iuCyOkWMUqrg0eBWQK1evZrJkyezbNky3nnnHYYPH05S4jW+X/05a7bM5fC1deAH1IElshlWpN8/4oi1DrkUSDv/O2kfEkmHzs/q92ZKqUJBg1sBdeDAgdQhsD58/22WH/4r64LOEZ/SPb8lcA9QAk4AZ9ZaHUGoWxciImgR8QDHWzagZv1M5o9RSqkCToObj7p69SqxsbEsWbKEq1evMmvWLOuW4nffwbp1tPl6aWreA3EXOFCa9N+b2dOZ1b7sTzvuIOH9ZyGyN9SuDVh/eB1TXylVWGlw81GXLl2iT58+APj5CRsidtL6232QZM3/0gCo3B7O1QeqAGLtV/OyP+3M7bS7vS0REYOo0yjcK/VXSilvKvLBTUT6Ac8CjbGuffYA/wImGWOS3XnshIQE1qxZQ2xsLO+88w6XT+9j3drPWHd4Dd/cOAjVgFOQnGyYdPonWiel379raVhV0p/wy7WIqPUA7doOpG7jCO2er5Qq8op0cBORicBQ4BqwCrgBdAA+AjqISG9jzE23HPzmTdq1acXGbTsAmHXkY87fa886ndLx8R4gDrgLLl4C9gING0J4OISHMymsJcXvCHFL9ZRSqiArssFNRHpiBbbTwAPGmP329qrAaiAKeA54Pz/HuXDhAitXrqRq6ZJE+AfC+vWwYQN89x2lysen5jt/wkCGvh0BTaF5fEnCg0OJbNsFPv4TuHTx16GHlVIqc0U2uAH/Y69fSQlsAMaYMyLyLLAGeFVEPszt7UmTnMy+LcsZ/8EYpn3+DQBl68HF/enz9SgJqyoB9YBQa6LO+66UJ7xcY8IbdeW+9k9QsnyVPJ+gUkoVVUUyuIlITawBphKBORnTjTFrReQkUAO4D1ifXXk//bCZxTFTuBF4mA2//MCG4uf5Jcikm/4l7iRcKA4Vr6dtGxhQhVV9A2hduRnhzaNo3rYvgSWCbz2AUkqpXCmSwQ1oZq93GWMSssizGSu4NSOb4LZ1+/c0bNrSukc4EnCdRLoqUBooC3InbPULIbJhZ2jdGsLCKHn77cSI5P9slFJKpVNUg1sde300mzzHMuTNlLH7gHAdOAlYj5FRMUFofb0y977YiAeadeXetv30FqNSSnmImNRP56JDRP4XGA18bowZkEWe0cD/Av80xjydIW0IMMT+9W5gpxurW5BUAs57uxI+QtsijbZFGm2LNPWNMaXdVXhRvXJLuReYp8hujPkn8E8AEdlijLnHqYoVZNoWabQt0mhbpNG2SCMiW9xZflF92jelD36pbPKkpMVnk0cppZQPKqrB7Yi9vj2bPLUy5FVKKVVAFNXgts1eNxSRoCzy3Jshb1b+6UyVCgVtizTaFmm0LdJoW6Rxa1sUyQ4lACLyPdAcGGiM+SxDWlush7hPAzXcPcakUkopZxXVKzeAd+31OBFJHaBRRKoAH9u/jtXAppRSBU+RvXIDEJGPsWYEuAasJG3g5DLAl0Avtw2crJRSym2K8pUbxpihQH/gJyASeARrrBEDPIo1eHKeiUg/EVknInEicllEtojIMBHx+XZ3su4iUlNEPhSRvSKSICLXRGS/iEwWkbruqL+TnP47ikiQiIwUkc0iclFErorIYRGZIyJtnK6/k9z5mhaRMSJi7GWEE/V1JyfaQkQCRaSDiPxNRL4TkVMikigiJ0VkrohEuPEUHOOG90j+yzPGFPkF+AdWQMu49MpHmRPtMhKAxcB84JK9bR7g7+3z9kTdsYYv+9Xe9zjWFfGXwAl7WzwQ5u1z9tTfEWvEm/32/meABcBsYBPWWKevefucPdUWGcq+F0gCku3yRnj7fD3RFkBHl8+bU3ZZs4AdLtvf9vb5evJ14VjberthfGEBBgPjgT7AnVidSfIc3ICeLi/Wei7bq2JdJRrgeW+ftyfqjjUup8HqGRXosj0QmGqn/eDt8/ZQW5QEDqR8YLm2h51eEbjL2+ftibbIUHZxYBfWAHbzfT24OdkWQHtgLhCeSdpjWAHfAO28fd6eeF042rbebhxfXBwIblvs/Z/IJK2tyx/Pz9vn6s66AyVI++/ztkzSq7ukB3v73N39d8TqxGSAf3v73LzdFhn2H2fv3w2YXgCCm8fe38CndnlTvX3enmgLRz9/vN04vrjkJ7gBNe19rwNBWeRJuSXnU7fjnK471tXZDTt/tUzSq9lpl7E7N/nK4oa2KIY1pqDAHCxTAAAIUklEQVQBQr19ft5siwz7tcK6Ovnc/t2ng5un39/AMLus5d4+d3e3hdPl+XzHhgIop9PpuOb1FY7W3RhzA1hl//qWiASmpNk//9X+daqxX7k+xOm/Ywus247HjTG7RSTM7kDxiYi8JSKt81thN3LLa1pESgD/Bn4Bns979TzK0+/vevb6lANlOc3ptnC0vKI6cLI7OTadjhe4o+5DgWXAU8BDLoOl3guUB94HXs5lPT3B6bZoZK/3i8h0YGCG9DdEJAb4fTZvbG9x12t6NFAfeNwYU1BGyvfY+1tEbgOetH+NyU9ZbuJ0WzhangY356UMuHwlmzyX7bXbpnvII8frbow5JCJhwGfAQ1i3HlJsAf5rX+H5GqfbooK9fgDwByYAk4EL9raPsb5MvwT8IbeVdTPHXxf2a+IF4EtjzKx81M3TPPL+FpEAYAZQFlhljFmU17LcyOm2cLQ8vS3pvHxNp+Nljtfd/hDbCYQA3bHms6qM9RxheSBGRN5w6ngOcrotUt5rAVi3YV82xhw0xlw0xizEag8DDPTBZ/8cbQt7PNd/YQXyoU6U6UGeen9PxhpQ4jiQ6ZyTPsDptnC0PA1uzivI0+k4WncRKYf1TFtpoLMxZqEx5oIx5rwxZgHQGetZltdFpF52ZXmB039H1zxTMiYaY7YA32O9JyNyUJ4nOd0WY4C7gD8bY3zxu6TsuP39LSLvA3/EGtu2gzHmdF7K8QB3vUccKU+Dm/OO2Ovbs8njq9PpHLHXTtW9C9ZV2nfGmEMZE40xB4CNWFczETmtpIccsddOtYVrnsNZ5EnZflsOyvOkI/baqbaIwnpYe6CIrHFdsP7hAXjW3vZpHurrTkfstVve3yLyN2A4cA4rsO3PbRkedMReO/0ecaQ8/c7Neemm08mic0BOp9PxNKfrXttex2WT56K9rpBNHm9wui22uvxcEevDK6NK9vpyJmne5I7XtB/Wc0tZqWsv5XJYnqe47f0tIuOBP2N9DxtpjPkp79X0CKfbwtHy9MrNYcaY41gfZMWA3hnT7el0amLdctjg2dplzw11/9let3B9DMClvECsLvKQ9dWMVzjdFsaYk1hXqWB9l5KxvPJYUzCB1dHGZ7ihLe4wxkhmC9ajAQAv29uaOncm+eeu97eIjMXqNfwrVmD7wZEKu5EbXhfOtq23HwT0xYUcPMSNNdrEHuDdTNJ6kfYkfYjL9ipYwwwZfHf4rVzXPau2sPe5Yu/zEVDcJa04MMlO+wUo6+1zd2db2GndSBtTsqnL9hLAF3baFnzsgXZ3tEU2x5mODz/E7abXxTv2Pr8CLbx9fl5uC8c+O73eOL6wYP3H/J3LkjJI5z7X7Rn2SXkTTs+izI9JG/xzEdaAn3H2tvn49sDJuap7dm2B9TxXyvh4J4GFdpk/29uuAY96+5w90RZ2+nukjcLwX7uMk/a2E7iMp+dri9NtkcUxUvbx2eDmZFtgzURi7GWznS+z5VVvn7OnXhe5LS/Lenm7YXxhwerMYH5ryc0fyM7TD/gWK1heweoNNwwfHFMyP3XPwYu1OdZzboexgtk14CDWuHkNvH2unmwLO08U8DXWf+rXsWYJ+BtQ2dvn6um2yGYfnw5uTrUF1kPav/nZA6zx9vl68nXhxGdnkZ6sVCmlVOGkHUqUUkoVOhrclFJKFToa3JRSShU6GtyUUkoVOhrclFJKFToa3JRSShU6GtyUUkoVOhrclCpERKSZiBgRmebtuijlTRrclCpcetjr+V6thVJepiOUKFWIiMgurKmGKhtjrnm7Pkp5i165KVVIiMhdQANgqQY2VdRpcFPKg+zvw4z985MiskVErojIaRGZKiKV7bQSIvKWiOwTkWsickxERmc2L56LnvZ6nsvxIuxjrrHLfEdEDohIgogcEpHXRMTfzlvLrsNJ+5g7RGSAu9pCKXfS25JKeVBKYAPGAy8Aa4F4IAy4DfgRaAMsB0Lt9OJYs1YHA1OMMUOyKHsT0ASoZIyJt7dFAKuxJne8CTTEmq8wCHjALnMyMAFrFParwCagBnC/XfQAY8znDpy+Uh6jwU0pD3IJbmeAdsaY3fb28lgBqD6wE7gIdDXGxNnpTbHm+/IH6hhjjmYotyZwDOuWZBeX7RFYwQ3gmwxlNnEpcy/wFfCSMeamnT4Ma5LZg8aYEAebQSm309uSSnnHGymBDcAY8yvWFRRY35sNSQlCdvp2YAkgWFdxGfWw0+ZlkgaQnEmZP9hl+mFdwY1MCWy2T7BmSb9TRGq7FiYiT4jIXvv25kERGZqTk1bKUzS4KeUdyzLZdsBeH3UNfC722+vqmaT1wLrtuDCL42VVZsoxvzbGJLomGGOSsCaYTXdMEemKNdHs60AZ4CngXRHpn8WxlfI4DW5KeceJTLZdzibNNb2E60YRqYT1/dg3xphzuTheXo/5PLDIGDPbGHPDGPM1VrB7PosylPI4DW5KeYExJjmb5OzSMvMo1vdmWd2SzEmZuTlmC2Bjhm2bgGYiop8pyifoC1Gpgi/KXn/poeOVBX7NsO1XIADruzulvE6Dm1IFmIiUBjoAW4wxxzx02DigfIZt5YEkrEcJlPI6DW5KFWxdsZ6D8+RYkt8DrTJsawls+43brUp5jAY3pQq2lIGSs/u+zWnvA91EpLeIBIpIO2CwvV0pn6DBTakCSkRKAJ2B3caYPZ46rjFmMVb3/9FYo6tMBf5HRzFRvkRHKFGqgBKR7lidSMYYY/7i7foo5UsCvF0BpVSeJQBvATO8XRGlfI1euSmllCp09Ds3pZRShY4GN6WUUoWOBjellFKFjgY3pZRShY4GN6WUUoWOBjellFKFjgY3pZRShc7/AzRBqQI98VFdAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"N = 500\n",
"t = np.linspace(0,4.5, N)\n",
"dt = t[1]-t[0]\n",
"\n",
"y0 = 0\n",
"v0 = 0\n",
"m0 = 0.25\n",
"u = 250\n",
"g = 9.81\n",
"c = 0.18e-3\n",
"dmdt=0.05\n",
"\n",
"rk2 = np.zeros([N,3])\n",
"heun = np.zeros([N,3])\n",
"\n",
"rk2[0,0] = y0 ; rk2[0,1] = v0 ; rk2[0,2] = m0\n",
"heun[0,0] = y0 ; heun[0,1] = v0 ; heun[0,2] = m0\n",
"\n",
"for i in range(N-1):\n",
" rk2[i+1] = rk2_step(rk2[i], rocket, dt)\n",
" heun[i+1] = heun_step(heun[i], rocket, dt)\n",
"\n",
"v = lambda mf: u*np.log(m0/mf)\n",
"\n",
"loc2 = np.where(heun[:,2]<=0.05)[0][0]\n",
"height2 = round(heun[loc2,0],2)\n",
"print(\"Height of detonation:\", height2, 'm')\n",
"\n",
"plt.plot(rk2[:,2]/m0,rk2[:,1],label='rk2',color='red',ls='-')\n",
"plt.plot(heun[:,2]/m0,heun[:,1],label='Heun\\'s Method',color='green',ls='--')\n",
"plt.plot(heun[:,2]/m0,v(heun[:,2]),label='Tsiolkovsky', color='black', ls=':')\n",
"plt.vlines(0.2,-100,225,lw=0.5,label='mf=0.05')\n",
"plt.xlabel(u'm/m\\u2080')\n",
"plt.ylabel('Velocity (m/s)')\n",
"plt.ylim(0,600)\n",
"plt.xlim(1,0)\n",
"plt.title('Velocity vs Mass Ratio')\n",
"plt.legend(loc='best',prop={'size':15});"
]
},
{
"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_m` 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": 8,
"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",
" \n",
" x = np.linspace(xmin,xmax,ns)\n",
" f = np.empty(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",
" 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\n",
"\n",
"def mod_secant(func,dx,x0,es=0.0001,maxit=50):\n",
" '''mod_secant: Modified secant root location zeroes\n",
" root,[fx,ea,iter]=mod_secant(func,dfunc,xr,es,maxit,p1,p2,...):\n",
" uses modified secant method to find the root of func\n",
" arguments:\n",
" ----------\n",
" func = name of function\n",
" dx = perturbation fraction\n",
" xr = initial guess\n",
" es = desired relative error (default = 0.0001 )\n",
" maxit = maximum allowable iterations (default = 50)\n",
" p1,p2,... = additional parameters used by function\n",
" returns:\n",
" --------\n",
" root = real root\n",
" fx = func evaluated at root\n",
" ea = approximate relative error ( )\n",
" iter = number of iterations'''\n",
"\n",
" iter = 0;\n",
" xr=x0\n",
" for iter in range(0,maxit):\n",
" xrold = xr;\n",
" dfunc=(func(xr+dx)-func(xr))/dx;\n",
" xr = xr - func(xr)/dfunc;\n",
" if xr != 0:\n",
" ea = abs((xr - xrold)/xr) * 100;\n",
" else:\n",
" ea = abs((xr - xrold)/1) * 100;\n",
" if ea <= es:\n",
" break\n",
" return xr,[func(xr),ea,iter]"
]
},
{
"cell_type": "code",
"execution_count": 9,
"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",
" time = (m0-0.05)/dmdt\n",
" t = np.linspace(0,time,1000)\n",
" dt = t[1] - t[0]\n",
" N = int(time/dt)\n",
" y0 = 0\n",
" v0 = 0\n",
" \n",
" rk2 = np.zeros([N,3])\n",
" rk2[0,0] = y0 \n",
" rk2[0,1] = v0\n",
" rk2[0,2] = m0\n",
" \n",
" for i in range(N-1): \n",
" rk2[i+1] = rk2_step(rk2[i], lambda state: rocket(state,dmdt=dmdt,u=250),dt)\n",
" pred_h = rk2[:,0] \n",
" error = 300 - pred_h[-1]\n",
" \n",
" return error"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"number of brackets: 1\n",
"\n",
"3a) Lower: 0.05 kg/s\n",
" Upper: 0.0889 kg/s\n",
"\n",
"3b) Actual root: 0.0789 kg/s\n"
]
}
],
"source": [
"interval = incsearch(f_m,0.05,0.4,ns=10)\n",
"root = mod_secant(f_m,0.00001,0.1)[0]\n",
"print('3a) Lower:', round(interval[1,0],4),'kg/s')\n",
"print(' Upper:', round(interval[0,0],4),'kg/s')\n",
"print()\n",
"print('3b) Actual root:', round(root,4),'kg/s')"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"3c)\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAbQAAAE0CAYAAABEhMaeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXwdVf3/8dcnSdukbbqkO13p3gItlLJDWzZlR1YRRBH8giCg+P3Koj8F5asgbvhVQFHZQQUU2UQ2oVAoSwuUAl1oS/em+5Kmzf75/TGT5CbNTXKTe3O39/PxuI/JnDkz87m3t/lkzpxzxtwdERGRdJeT7ABERETiQQlNREQyghKaiIhkBCU0ERHJCEpoIiKSEZTQREQkIyihSdows4vMzMPXjHjXb0888T52Rxw/FTT6N2rLa3nEse4Ly15N3juSZFJCE8lgZnZT41/8IplKCU1EkukhoDDK65aIevtEqTOxI4OV1JaX7ABE0pW73wfcl+Qw0pq7VwE7m9pmZhURq7vcvcl6Ece6CLgobsFJ2tEVmoiIZAQlNMl6ZlZgZt82s5lmttHMKsys2Mz+aWYnNbNfi502zKyTmf23mb1vZqVmtsXMXjezC8Ptrb7HZWZ9zeznZvapmZWZ2WYze87Mjm6i7owwrhvDouFNdKh4taVzhsc6KWKfw1qoO9zMasK6/9Wa48dLc51CGv9bmdleZvYbM1tqZrvNbIWZ/cHM9orYJ8fMLjazN81sq5ntNLNZZnZqK2LJC/f9d/hdqgi/Wy+Y2QVmZnF98wKoyVGynJlNAp4ChjfaNAA4HTjdzO4F/svdq2M8dg/gBeCQiOKuwJHAkWZ2LLC8lceaGB5rcERxF+AE4PNmdpG7PxBLfDF4AdgA9Ae+DMxupu75gAHlwGMJiqddzGx/gvfUL6J4GHApwWd5BLAJ+BvBdyDSEcBTZnaJu98T5fjDCL5Tkxtt6gscH74uMLNz3L20ve9H6ukKTbJW+IvnFYJkthK4DBgNFAH7ArcB1cDXgB+14RR/pj6Z3QdMAfqEy/uBrwIXtvJYTwOVwFeAoQS/jM8AVhEkkDvMrE9E/ddp2LFiJXt2qDixNScO73P9NVw918w6NVP9gnD5rLtva83xk+AfwHbgXGAQMAT4NlBF8F24BfgZwedzEzCe4N/tGGBBeIzfNPq8gbo/Yv5DkMw2A/8NTAB6A2OBG4Dd4bH/kIg3l9XcXS+90uJFcMPfw9eJQPcWXpdF1J/RxPGeCretBvpFOeclYZ1yYK9o8TSx32ER5/5tlGPfFVFneQvvdzXQv4k6UyLqfKOJ7TdFO36Mn/3UiPOcGqXO/hF1vhCHf++bIo43ohX17wvrvtrCZ7kS6NtEnf8Nt1cBNcDZTdQZS/BHjgOXNbH9/8JtO4BxUeI8LiKWqcn4v5SpL12hSbr6F1DSwuv30XY2s5HAKeHqt919Y5Sq9wBLgc7AOTHE95VwuQv4fpQ61xP8td4aP3b3DY0L3f094MNw9aAY4ouJu88BFoarF0Sp9uVwuYXg3ydV/djdNzVRXnsVmgu84e6PN67g7ouB98PVyKZkzKwbwR9AtedY1NTJ3f0lgqs4iP5ZShsooUm2Opagqa4GeNPMujf1AroB88J9psZw/MPD5avuvqOpCu6+HZjZyuM918y22l+cA1t5rLZ6KFyeZmaFkRvMLAc4L1x9zN0rSF3PRylf2oo6kfUaf96HE9wjBZgZ7TsVfq9q/wiJ5TslLVBCk3R1tLtbcy+Ce1/RjAuXOcAamr/SOzOs24/WGxEum/wrPcLCFrbXWtvMtl3hsmszdeLhYYJmsgLgrEbbjqa+w8pDpLYmP0t3j7xaXtfM/rX1ChqVj4v4+R2a/059O6wXy3dKWqCEJtmqZxv2yY+hbrdw2VIvtmYHC9fy1vWwTGhXcHdfDswKV7/caHNt09lnwBuJjKO9WvlZtuXzTvR3SlqgbvuSrWoTySZ3T8RfyaVAD+oTWzTdE3DuRHoIOAo42sz2cve1ZpZP/RXbwx72fMhCkX+cFHoLM5tI/OkKTbLVsnDZ18yGJOD4K8Ll2BbqjWthe6p5lKDHZw7wpbDsVILkDanf3JhIyyJ+3j9pUWQxJTTJVi9G/Nzcvba2qm12O7pxB4pa4ZilGQk4d6TKcJkbj4N5MLbs2XD1y42Wc6L17MsSMwmSPSTmOyUtUEKTrOTuC6nvWn6DmR3eXH0z629mvWM4xYPhsitwc5Q6t7Bnx4J42xwu+5lZvG4x1L63/c3sKILZSiLLs1LYm/XP4epFZta440wDZtbDzAYlPrLsoYQm2ewbBFM6FQCvmNkvzezQcM7EPmY2wczON7O/EDQhjmrtgd39TYIZKQC+ZWZ/NrPJZlZkZvub2T3AFTRspkqEueGyC/DjcA7DTuFcg229avsXwVgzgAcIxuhFziaSzb5H0LM1B3gs/HefbmYDzKy3mY0xs7PM7E8Es7wckdRoM4wSmmQtd18FTCeYzqgz8B2CeQo3Eszl9wlBV/XzCBJCZdNHiupi4N2Inz8guGJ6n6BJ6gHq7zlVtfV9NMfd3wXeDFdvIBiiUEHwXl5u4zErCO6lQf3whBebGvidbcKxhUcTNDkbwb/7q0AxwR8Bi4HHCQZg9yD4t5A4UUKTrBY2PU4imFfxaYIxShUE90JWAf8GrgaGuvu8aMeJcuztBBMRf5dgIO1uYBtBgrnY3b9KfS/Hkna/mehOAn5OkKB3tVC3tRo3L2ZzZ5AG3H0dQU/QMwgmaF4JlBF8r9YRzB96PTDG3Z9KVpyZyLK3h61I8pnZk8BpwDPu3uJjSUQkOl2hiSRJOAVS7bPM5jZXV0RapoQmkiDhvH2dm6lyG8FjXKD+npSItJESmkji7AssNLPrzOyAsIfjADM7zsyeAi4P6/3F3T9JYpwiGUH30EQSxMwOpfmnO0MwGPf0sAOJiLSDElqS9O3b10eMGJHsMCSBqqur2bJlCzt27GD37t1UVVVRU1NDbm4uXbt2paioiKKiIswSOqewSEaZO3du1PlXNTlxkowYMYI5c+YkOwwRkbRiZiuibdM9NBERyQhKaCIikhGU0EREJCMooYmISEZQQhMRkYyghCYiIh1n3TqYPh2Ki+N+aCU0ERHpEO+v3Mqyb12Hz5qF/+hHcT++BlYnydSpU13j0EQkaxQUQFnZnuX5+bB7d6sPY2Zz3X1qU9t0hSYiIom3bBn/mXIcu/O6AFCdXwAXXACffRa3UyihiYhIwpX07svamk50qaqgLLcTORXl0KMHDBwYt3No6isREUm4D1Zto++ubTx0wIm8fdzZ3FE6J+ggEkdKaCIiknBzV2zl9jO+D8CFhw6HL3wt7udQk6OIiCTceyu31f184PDeCTmHEpqIiCRUTY3z/oqtdetKaCIikpY+3bCTkvIqAPoVdmFI74KEnEcJTUREEmpuxNXZlGG9EvZQ25RJaGY2zsy+ZWYPmdlCM6sxMzezs1ux7/lm9rqZbTeznWY2x8y+aWbNvr+O3k9EJBu989nmup8T1dwIqdXL8XLgW7HuZGZ3AFcAZcDLQCVwLPA74FgzO8fdq5O9n4hINnJ33lq2pW79kL37JOxcqXRF8RHwc+CLwGhgZks7mNlZBMmlGJjk7qe4+xnAGGABcAZwZbL3ExHJVss376J4RzDlVWGXPPbZq0fCzpUyCc3d/+Tu17r7o+6+tJW73RAur3P3TyOOtZ7gig/g+iaaAjt6PxGRrDR7aX1z48F7F5GXm7hfj2n7i9fMhgAHAhXAY423u/tMYA0wEDg0WfuJiGSz2cvqE9phoxLX3AhpnNCAA8Llx+4ebarmdxvVTcZ+IiJZKbh/Vp/QDh2phBbN3uFyRTN1Vjaqm4z9RESy0tKNpWwsKQegR34eEwYl7v4ZpHdC6x4uS5upszNcFiZxvzpmdmnYxX/Oxo0bmzmMiEj6e3PpprqfDxnZh9ycxIw/q5XOCa32k4n1CaUdvV8dd7/b3ae6+9R+/fq19TAiImlh5qL6P9yPGtM34edL54RWEi67N1OndltJRFlH7yciknXKq6p5M6KH4/Sxif8jPp0T2vJwObyZOkMb1U3GfiIiWWfu8q3srgzmmBjRpyvD+3RL+DnTOaG9Hy73MbNoM10e1KhuMvYTEck6MxfXNzdO64CrM0jjhObuq4D3gM7AOY23m9l0YAjBrB6zk7WfiEg2ikxoHdHcCGmc0EK3hMufmdno2kIz6w/cGa7e6u41Sd5PRCRrFG8vY2Fx0JWgc25Owsef1UqZyYnNbAr1SQFgYrj8qZn9T22hux8a8fPjZnYXwbRT883sJeonC+4B/JNg0uAGOno/EZFs8tKC9XU/H7R3b7p16ZhUkzIJjSAhHNJE+ZjmdnL3K8xsFvBNYDqQCywE7gHuina11NH7iYhkixc+qU9ox00Y0GHnTZmE5u6vUj/WK9Z9HwEeSfX9REQyXUlZJbMjBlQfP7HjElq630MTEZEUMnPxRiqrg/knJg7qwZDeXTvs3EpoIiISNy9GNDd25NUZKKGJiEicVFTV8MrCDXXrSmgiIpKWZi3ZyI6yKgAG9ypI6NOpm6KEJiIicfH0vHV1P58yaRBmiZ1dvzElNBERabeyyuoG989OnbxXh8eghCYiIu326qIN7CwPmhv37tutw5sbQQlNRETi4OkPk9vcCEpoIiLSTqXlVby8ILnNjaCEJiIi7fTSgvWUVQaz/o0bUMjYAYVJiUMJTURE2uXxuavrfj518qCkxaGEJiIibbZm225mLQnmbjSDM6YMSVosSmgiItJmf5+7Gg+mbuTI0X0Z3KsgabEooYmISJvU1DiPzV1Vt37u1KFJjEYJTURE2uitzzazastuAHrk53X43I2NKaGJiEibPDanvjPIFw4YTH6n3CRGo4QmIiJtsLW0gmfn1w+mTnZzIyihiYhIG/xtzioqqoKxZ/sO7pGUqa4ay4u2wczuidM53N0vidOxREQkyaprnAdnr6hb/+phI5Iy1VVjURMacBHgQHujdEAJTUQkQ7y8YD1rtgWdQXp37ZS0qa4aay6hAbwB/Lkdx/86cHg79hcRkRTzQMTV2XkHD0t6Z5BaLSW0Je5+f1sPbmYzUEITEckYSzaU1M0MkmNwwSHDkhxRveY6hewAdrXz+LvD44iISAb442uf1f183IQBDOndNYnRNBT1Cs3de7X34O5+BXBFe48jIiLJt35HGU+8v6Zu/b+mjUxiNHtSt30REWmVe2Z9RkV10FX/wOG9OWhEUZIjakgJTUREWrR9dyUPv72ybv0b00clMZqmKaGJiEiLHn57BTvLqwAY3b87x47vn+SI9hRTQjOz0Wb2RzNbYma7zKw6yqsqUQGLiEjHKi2v4s+v13cGuXTaSHJykj+QurGWuu3XMbOpwH+AbrQ82Dr13qmIiLTJA7NXsLm0AoDBvQr4wv6DkxxR02K5QrsN6A48CkwBCt09J9orIdGKiEiHKimr5A+vLa1bv/KY0XTOS81f8a2+QgMOARa4+5cSFYyIiKSW+99czrZdlQAMLSrg7AOHJDmi6GJJs7uBeYkKREREUsuOskr+GHHv7KpjxtApNzWvziC2hPYOkFqj6EREJGHueGUJ23cHV2fD+3TlzANS895ZrVgS2k+AA8zszEQFIyIiqWHVll3c+8byuvXvHD+WvBS+OoMY7qG5+xtmdh7wRzM7A3geWA3URKn/WnxCFBGRjvaLFxbVPcBz8tBenJYij4hpTiydQgA6E0xYfH74isbbcGwREUkB81Zt48kP1tatf/+kCSnxAM+WxDIO7SzgYYJmys3AcmBnYsISEZFkcHd+8q8Fdeuf32cAB++dWnM2RhPLVdT3CAZMXwHc7e5NNjWKiEj6+tf8Yt75bAsAeTnGdSeMT3JErRdLQhsPvOHuv09UMCIikjwlZZX8+JmP69a/fOhwRvbrnsSIYhNLl5XtBJ1AREQkA93+0qes31EOQN/uXbjm+LFJjig2sSS0F4CDLB3uDIqISEw+WbuD+95cXrf+g1Mm0LOgU/ICaoNYEtr3gULgF2amHowiIhmipsb5f/+cT3WNA3D4qD5p0U2/sVgS0yXAv4BvA2eY2X+IPg7N3f3mOMQnIiIJ9vDbK3hv5TYAOuUaPz5937Topt9YLAntJoLxZQaMAC5uok7tdgeU0EREUtyqLbu45bmFdeuXThvJ6P7p0xEkUiwJ7ccEiUpERDJATY3z3cfnsauiGoAx/btz9bFjkhxV28Uy9dVNCYxDREQ62ENvr+CtZcGYsxyDn58zmS55uUmOqu1Se6ZJERFJiBWbS7k1oqnxsumj2H9oryRG1H5KaCIiWaaiqoar//J+g6bGbx+Xvk2NtaImNDM708ymtOfgZjZFj5sREUktv3pxMfNWbweC6a1+eW56NzXWau4K7XHgynYe/yrgsXYeQ0RE4mTWp5v4w2tL69avPWEck4akd1NjLTU5iohkiU07y7nm0Q/wsL/6UWP68vUjRyY3qDhqqZfjCeEA6rZKn2maRUQyWFV1DVc98j4bS2rnauzML8+dTE5O+g2gjqalhDYwfLWHxq6JiCTZbc8vYvayzQCYwS/OmUz/wvwkRxVfzSW0ozssChERSZhnPlzL3a8tq1u/5rixzBjXP4kRJUbUhObuMzsyEBERib/F60u49vEP69aPm9CfK48encSIEkedQkREMtTW0goufWBO3XizEX268stz98+o+2aRlNBERDJQeVU1lz04l+WbdwFQ0CmXP1w4Ne2ecRYLJTQRkQzj7lz/9/m8s3xLXdmvzp3MuIGFSYwq8ZTQREQyzG9e/pQn3l9Tt37DieM5cb9BSYyoYyihiYhkkL/PXc3tL31at/6lg4dx6bTMGTzdHCU0EZEM8eIn67n27/U9Go8a05cfn75PWj59ui2U0EREMsDspZv55iPvUV0TzGUxfmAhd1wwhU652fNrvtXv1My+YmaHt6LeoWb2lfaFJSIirTV/9Xb+64E5VFTVADCsqCsPXHwwPfIzt0djU2JJ3fcBX29FvUuAe9sUjYiIxGTx+hK+eu877CyvAqB/YRceuuQQ+vfIrGmtWiMR16LZ0VgrIpJki9eX8KW732JLaQUAPQs68eAlhzCsT9ckR5YciUhoQ4CdCTiuiIiEapPZ5jCZde+Sx71fOyjjx5o1p9nZ9pu4Fza6mftjecAE4Fjg3TjEJiIiTWgqmd1/8cFMGdY7yZElV0uPj7mPho9/OSJ8RWNADfCL9oUlIiJNWVi8gwv++PYeyezA4dmdzKDlhPYA9Qntq8BS4I0odSuANcCT7j4vPuGJiEituSu28rV732FHWdABRMmsoWYTmrtfVPuzmX0VmOXuFyc6KBERaWjm4o1848G57K4MZs4v7JLH/ZeomTFSS1dokfZGnT1ERDrcMx+u5Zq/fUBlddBg1qdbZ+6/+GD2HdwzyZGlllYnNHdfkchARERkTw+9tYIfPPkRHt78GdyrgAcvOZiR/bonN7AUFMsVGgBmlg9MBfYCoo7cc/cH2hGXiEhWq6lxbnt+Eb+fubSubFS/bjx4ySHs1asgiZGlrpgSmpldA/wQ6NGK6kpoIiJtUFZZzX8/No9nP1xXVzZpSE/u+9rBFHXrnMTIUlurE5qZXQz8MlxdACwEdiQiqFiY2X0EPTCjWeTu46Psez5wOTAJyCV4T/cCd7l7TTPnbNN+IiIt2VJawaUPzGHOiq11ZceO78//fekAunWJuVEtq8Ty6VxN0IX/Qnd/JEHxtMcbwJImytc1UYaZ3QFcAZQBLwOVBIPCfwcca2bnuHt1vPYTEWnJ8k2lfO2+d/lsU2ld2VcOG86Np+5Dbo5mFWxJLAltLPBmiiYzgD+5+32tqWhmZxEkpWJgmrt/GpYPAF4BzgCuBH4Tj/1ERFry2uKNXPnIe3VjzMzg+ydN4JIj986a55m1VyxzOe4CViYqkA52Q7i8rjYpAbj7eoKmRIDrzazx59PW/UREmuTu/PG1ZVwUMWC6S14Od10wha8fNVLJLAaxXKG9CeybqEA6ipkNAQ4kmNnkscbb3X2mma0BBgOHErzvNu8nIhJNWWU1N/xjPk+8v6aubGCPfP5w4YFMHtoriZGlp1gS2o+AN83sq+5+f6ICaoejzWwS0B1YD8wCXmyik8YB4fJjd98d5VjvEiSmA6hPTG3dT0RkD+u27+ayB+fy4ertdWUHDu/NXV+eQv/C7HuWWTxETWhmNq2J4l8B95jZScCzBE2QTfbqc/fX4hJh6zX1FIBPzOw8d58fUbZ3uGxuoHht0+reEWVt3U9EpIE5y7fwjYfeY9PO8rqyL04dyo+/sA9d8nKTGFl6a+4K7VUazrRfy4Czw1c03sKx4+kDYC5Bj8MVBGPkpgA/ASYDL5nZFHevvaavHV5f2vhAEWqn+Ip8sFBb96tjZpcClwIMGzasmcOISKb6yzsr+eGTH9VNY5WXY/zw1IlceOhw3S9rp+aSzms0ndBSirvf3qioFHjWzF4EZhLcz7qBoPch1D9RO9b31tb96rj73cDdAFOnTk35z1ZE4qeiqoabn/mEB9+qb+Qp6taZO86fwmGj+iQxsswRNaG5+4wOjCPu3L3CzG4BngROithUEi6bmwitdltJRFlb9xORLLd+RxlXPPwecyMGS08Y1IO7LzyQoUVdkxhZZsn0YecLw+XgiLLl4XJ4M/sNbVS3PfuJSBZ7e9lmvvnI+w3ul508aRA/P3sSXTtn+q/gjpXpn2btdXzkY2/eD5f7mFlBlB6LBzWq2579RCQLuTv3vLGcn/5rAdU1wR2GHIPrThjPpdM0viwRYpnLsalej02pADa5e1PTUHW0c8Plu7UF7r7KzN4j6DhyDo0mUTaz6cAQgtlAZrd3PxHJPqXlVVz/j/k8PW9tXVlRt8787ksHcPjovkmMLLPFcoX2KjF0iDCzHcD9wA/cPSH3lMxsf4Ik8lzk/Ilmlkcw9+TVYdGvG+16C8Hg6J+Z2Zu1ydfM+gN3hnVubWIMW1v3E5Es8dmmUi57cA6L19c3DE0e2ou7Lpiix74kWCwJ7bWw/uHh+lbqx6ENB4oIEt7bQD9gBHAVMMPMDnf3XXGKOdII4Algi5ktBlYTdJnfj+B5bTUE01Q9H7mTuz9uZncRTFc138xeon6S4R7APwkmGyYe+4lIdnjh42L++9F5lJRX1ZWdf8gwbjx1osaXdYBYEtoJwEvAJ8D/uPu/Izea2eeBnxMktf2AgcCDBAnwauDWeATcyDyCiYAPJkiqB4TnX03wOJc73H1uUzu6+xVmNgv4JjCd+sfA3EMzj4Fp634ikrmqa5xfvbiIO16pfxhn57wc/vcL+3Lu1KHN7CnxZO6ta0U0s/8luOIa4+4botQZACwmSCTfM7OhwCKC6aIOamqfbDV16lSfM2dOssMQkXbaWlrB1X99n9c/3VRXNrhXAX+48ED2HdwziZFlJjOb6+5Tm9oWy6zwXwReiZbMoG7W+VcIO2O4+yrgPYJHz4iIZJQPVm3jlN/OapDMpo3txzNXHalklgSxNDkOIUhOLSmn4bivVUCT2VREJB25O/e9GXTJr53CCuCqY0bz7ePG6mGcSRJLQtsETGtmDBZmVgBMAzZHFPcGtrU9RBGR1LGjrJLrHv+Q5z4qrisrzM/jV+fuz/ETByQxMomlyfFpYADwaHhvrIGw7G9Af+CpiE3jgWXtCVJEJBV8tGY7p/52VoNktt/gnjx71VFKZikgliu0G4ETgZOBJWY2m2B2eyfoYXg40CksuxHAzA4EhtFoELKISDpxdx55ZyU/evoTKqrqOzJ/5bDhfP/kCeqSnyJandDcfaOZHQ7cBZxK0LTYoArwDHC5u28M95lrZp0iBz2LiKST0vIqvvfEfJ78oH7Wj+5d8rj1rP04ZdJeSYxMGotpLkd3Xwd8wcyGESS02s4fa4HX3X15E/somYlIWlpUXMIVD89l6cb6xyCOH1jInRdMYWS/5h68IcnQpsmJ3X0l8FCcYxERSQnuzt/eXcVNT39MWWV9E+N5Bw3lptP2Ib+TmhhTUabPti8iEpPtuyv53j/m8+z8dXVlBZ1y+ckZ+3LmlCFJjExaEjWhhc2KAGvcvTpivVXCqzgRkbQxd8UWrv7LB6zZVj8yaUz/7tx5wRTGDChMYmTSGs1doS0nmNx3IsF0Vstp/Wz73sKxRURSRnWNc9erS/j1S5/WPbsM4IJDhvH/Tp5IQWc1MaaD5pLOSoLEVNloXUQkYxRvL+Pbf3uft5ZtqSvrkZ/Hz86axIn7DUpiZBKrqAnN3Uc0ty4iku5e+mQ93318Hlt3VdaVHTSiN7efdwCD9eyytKNmQRHJOrsrqrnluQU8MHtFXVmOwVXHjOGqY0aTlxvLJEqSKpTQRCSrfLh6G9f87YMGY8sG9czn11/cn0NH9kliZNJeMSc0MxsNXAYcRvBk6ifd/dpw26HAJOBRd9eExCKSMqqqa7jr1aX85uVPqYro+PG5iQP42VmT6N2tcxKjk3iIKaGZ2SXAHUDtv7wDfSOq9COYGquS4InRIiJJt2JzKdf87QPeW1n/d3bXzrnceOpEzp06FDM97iUTtDqhmdkRwB+AncD3gdeAtxtV+zewAzgNJTQRSTJ356/vruLmZz5hV0X9LHwHDu/Nr86dzPA+3ZIYncRbLFdo1xJckZ3o7rOBPf6qcfdKM1sETIhbhCIibbCxpJwb/vEhLy3YUFeWl2Ncc/xYLps2Uh0/MlAsCe0w4J3aZNaMVSihiUgS/fujdXz/iY/YXFpRVza6f3du/+L+7Du4ZxIjk0SKJaH1BFa3ol7nGI8rIhIXm3eWc+NTH/PMh+salF90+AiuP3G8JhXOcLEkng3A3q2oNw5Y07ZwRETa5l/z1/GDfza8KhvQowu/OGcyR43pl8TIpKPEktDeAM42s6nuPqepCmZ2PDAW+FM8ghMRacmmneXc+OTHDWbHBzh36hC+f/JEehZ0SlJk0tFiSWi/Bs4B/mFmXwdeitxoZtOAe4Aq4Ldxi1BEpAnuzrPz1/HDJz9mS8RV2aCe+dxy5n7MGNc/idFJMrQ6obn722Z2LfBz4DmC7vlO8ATrkwnGoxnwHZqB3GsAABUySURBVHefn4hgRUQANpSUceOTH/PcR8UNys87aCjfO3kCPfJ1VZaNYuq84e6/NLOPgR8BUwkSWK9w83zgB+7+VHxDFBEJ1NQE48pueW4BJWVVdeV79cznlrMmMX2s7pVls5h7I7r7v4F/m1kfgk4iucAqd18b7+BERGot2VDCDf+Yz7vLtzYo/9LBQ7nhJF2VSTu617v7ZmBzHGMREdlDeVU1d76ylDtfXUJldf0cjCP6dOUnZ+zHEaP7NrO3ZBONFxORlPX2ss3c8MR8lkXMjJ+XY1w2fSRXHTNG48qkgagJzcy+0p4Du/sD7dlfRLLX1tIKfvbvhfz13VUNyvcf2otbz9qP8QN7JCkySWXNXaHdR9CLsa2U0EQkJrWdPm57fiHbIp4i3b1LHteeMI4LDhlObo5mxpemNZfQXiN6QpsOrAcWxj0iEclK81Zt44dPfsS81dsblH9u4gB+dPo+DOpZkKTIJF1ETWjuPiPaNjOrAZ5z94sTEZSIZI+tpRXc9vwi/vruSjziT+ghvQv44SkT+dw+A5MXnKQVdQoRkaSI1rzYOS+Hb0wfxRUzRqnTh8RECU1EOtyc5Vu4+ZlP9mhePHpcP246bR89eFPaRAlNRDrMqi27uPXfC3m20eNdhvQu4MZT9+G4Cf33eHCwSGspoYlIwu0sr+LOV5bwp1mfUVFVU1eu5kWJJyU0EUmY6hrnsTmr+MULi9m0s7zBtpMnDeL6E8YztKhrkqKTTKOEJiIJ8eaSTdz87AIWrNvRoHzykJ784JSJTB1RlKTIJFM1N1PItBb2HdhcHXd/rc1RiUja+mjNdm57fhGvLd7YoHxgj3yuO3Ecp08eTI4GR0sCNHeF9irRB1Y78PnwFW27rv5EssiKzaX84oXFPD2v4YM3Cjrlctn0kVw6bSRdO+vXgiROc9+ulbRv6isRyQIbSsr47ctL+Ms7K6mqqf+VkWNw1pQhfOdzYzXLh3SI5mYKGdGBcYhImikpq+Tu15bxp9c/Y3dldYNtn5s4gP/5/DjGDihMUnSSjXT9LyIxKSmr5IHZK/jj68sazPABcPDeRVx3wngOHN47SdFJNlNCE5FW2Vlexf1vLm8ykY0fWMh1J4xnxrh+GhgtSaOEJiLNai6RDSvqyjXHj1HPRUkJSmgi0qSWEtlVx4zmCwcMplNuTpIiFGlICU1EGti8s5z731zO/bNXsH33nonsymNGc4YSmaQgJTQRAWDNtt388bVl/PXdlZRV1jTYNrSogKuOHsMZU5TIJHUpoYlkucXrS/j9zKU89cHaBuPIAIb36coVM0Zx5pQhSmSS8pTQRLKQuzN3xVb+8NoyXvxk/R7bJw7qweUzRnHSfoPIVWcPSRNKaCJZpKKqhn/NX8c9b3zGh40erglw6MgiLp8xmmlj+qr7vaQdJTSRLLCltIJH3l7BA7NXsKGkfI/tn5s4gG/MGMWUYRoQLelLCU0kgy0qLuHeNz7jiffXUF7VsKNH57wczjxgMJccuTdjNEWVZAAlNJEMU1FVwwufFPPwWyuZvWzzHtv7F3bhK4cN50sHD6NP9y5JiFAkMZTQRDLEqi27+Ms7K3l0zio27azYY/t+g3tyyZF7c9J+g+icpx6LknmU0ETSWFV1Da8s2sjDb69g5uKNeKMHPuXmGJ+bOIBLjtybA4f3VkcPyWhKaCJpaMXmUv7+3hoen7OKtdvL9tg+sEc+5x08lC8eNFTPIpOsoYQmkiZ2llfxr/nreHzuat75bMse281g2ph+XHDIMI4Z3588DYSWLKOEJpLCamqctz7bzONzV/Pc/OI9HqQJ0KdbZ849aChfOmgYw/p0TUKUIqlBCU0kBS0qLuGpeWt48oO1rN66e4/tuTnGjLH9OGfqEI4ZP0CdPERQQhNJGSs2l/L0vLU8NW8ti9fvbLLO2AHdOefAoZx+wF70L8zv4AhFUpsSmkgSFW8v45kP1/L0vLXMa2IqKoCeBZ04ff+9OPvAIew3uKd6KopEoYQm0sFWbC7l+Y+Lef7j9by3cuseXe0B8jvlcOyEAZw2eS+mj+1Hfqfcjg9UJM0ooYkkmLuzYF1JmMSKWVhc0mS9TrnGtDH9OG3/vThuwgC6ddF/T5FY6H+MSAJUVNUwZ8UW/rNgA89/UsyqLXt27ADIMTh0ZB9Om7wXJ+w7kF5dO3dwpCKZQwlNJE7W7yjj1UUbeGXhRmYt2cTO8qom63XOy+Go0X35/D4DOXZCf82nKBInSmgibVRZXcO8Vdt4ddFGXlm0gY/X7ohat7BLHkeP78/n9xnI9HH96K7mRJG40/8qkVaqqXEWrS/hjSWbeHPpZt5etpnSij0HOtca3KuAY8b355gJ/Tl8VB+65Kljh0giKaGJNGPl5l28sXQTbyzZxOylm9lcuucs9rXycoyDRhRxzPj+HD2+H6P6dVcXe5EOpIQmEqqucRYW72Duiq28u3wrc5dvaXLi30iDexVwxOg+HDO+P0eM7kthfqcOilZEGlNCk6y1q6KKD1ZtY87yrcxZsZX3V2ylJEpHjlpF3Tpz2Mg+HD66D0eM6svwPl11FSaSIpTQJCtUVNWweH0J81Zv48NV2/lwzXYWry+huqaJUc0RunXO5eC9izhidF8OH9WX8QMLyclRAhNJRUpoknEqqmpYunEnH6/dwYertzFv9XYWrNtBRVVNi/sO6NGFqSOKOGh4b6aOKGL8wEI9hkUkTSihSdpydzaUlLNg3Q4WFpewMFwu2bCTqhauvCB4ftiY/t2DBDaiN1OHFzGkd4GaEEXSlBJaG5nZ+cDlwCQgF1gI3Avc5e4tXwpIq9XUOGu27Wbpxp0s21jK0o07WbpxJ4uKS9i6q7LVxxlaVMCkIb2YNLgnk4b0Yt/BPdSJQySDKKG1gZndAVwBlAEvA5XAscDvgGPN7Bx3jz5ASfZQU+OsLylj1ZbdrNqyixWbS1kaJq/PNpVS3ormwkhDehcwfmAPJg/pyX5DggRW1E3TSolkMiW0GJnZWQTJrBiY5u6fhuUDgFeAM4Argd8kLcgUVFVdw8ad5azfUc6arbtZtXUXK7fsYtWWXazeups1W3dTUR37hW33LnmMH1jI+EGFjB/YgwmDChk7oFBXXiJZSAktdjeEy+tqkxmAu683s8uBV4Hrzey32dD0uLuimi27Ktiys4LNpeVsKCln/fYy1peUUby9nPU7yli/o4xNO8tpxW2tqPp068yoft0Z2a9b3XLsgELd8xKROkpoMTCzIcCBQAXwWOPt7j7TzNYAg4FDgTc7NsLY1dQ4uyqr2VlWRUlZJSXlVZSUVdWt7yyvYkf487ZdlWwurWBraQVbwtfuyvi1rPbu2omhRV0Z2rsrQ4u6MqpfN0b2686oft00C72ItEgJLTYHhMuP3b3p54HAuwQJ7QDinNBKyir5+fOLqHGnuiZIRtXu1LiHPwdlwXanqsYpr6qmrLJmj2VZZTXlVTWt6soeL327d6Z/YT579cpnSJi0hvYuYGhRV4b0LlAzoYi0ixJabPYOlyuaqbOyUd24Kaus4YHZzZ2643XKNYq6daaoWxeKunWiX/cuDOiZz4DCfAb2zGdAj3wG9OhC/8J8OudpPJeIJI4SWmy6h8vSZursDJeFjTeY2aXApQDDhg2L+eS5CZqhoqBTLoX5eXTPz6MwvxOFXfLo3iVvj7JeXTuFyav+1b1Lnu5hiUhKUEKLTe1v7jZ1b3D3u4G7AaZOnRrzMbp2zuWmUyeSm2OYGbk5Rq4ZOTlGbg7kmJETlueYkZdj5HfKJb9TDl3y9lx26ZRD59wcTeUkIhlBCS02JeGyezN1areVNFOnTfI75XLREXFvyRQRyQi6qRGb5eFyeDN1hjaqKyIiHUAJLTbvh8t9zKwgSp2DGtUVEZEOoIQWA3dfBbwHdAbOabzdzKYDQwhmEZndsdGJiGQ3JbTY3RIuf2Zmo2sLzaw/cGe4ems2zBIiIpJK1CkkRu7+uJndRTDT/nwze4n6yYl7AP8kmKRYREQ6kBJaG7j7FWY2C/gmMJ36x8fcgx4fIyKSFEpobeTujwCPJDsOEREJmHs7pkCXNjOzjTQ/hVZz+gKb4hhONtBnFht9XrHR5xWb9nxew929X1MblNDSkJnNcfepyY4jnegzi40+r9jo84pNoj4v9XIUEZGMoIQmIiIZQQktPd2d7ADSkD6z2Ojzio0+r9gk5PPSPTQREckIukITEZGMoIQmIiIZQQktjZjZ+Wb2upltN7OdZjbHzL5pZvp3jGBm48zsW2b2kJktNLMaM3MzOzvZsaUaM+tkZsea2S/N7C0zW2dmFWa2xsweN7MZyY4x1ZjZVWb2qJktMLPNZlZpZhvN7CUz+7LpEe4tMrOfhv8n3cz+J27H1T209GBmdwBXAGXAy9TPH1kIPAGc4+7VyYswdZjZ7cC3mth0jrs/3tHxpDIzOw54MVwtBuYCpcBEYN+w/GZ3/2ESwktJZrYa6A98BKwh+LyGA4cQPNX+SeBMTYHXNDM7iOBpJDkEn9d33f0X8Ti2/rJPA2Z2FkEyKwYmufsp7n4GMAZYAJwBXJnEEFPNR8DPgS8Co4GZyQ0npdUAfwemufug8Lv1RXffDzgPqAZ+YGZHJzXK1HIe0Nvdp7j7qe5+nrsfBuwHrAdOB76a1AhTlJl1Ae4j+JyejPfxldDSww3h8jp3/7S20N3XE8z6D3C9mh4D7v4nd7/W3R9196XJjieVuft/3P1sd3+9iW1/I/jlA/DlDg0shbn7LHcvbaL8Y+COcPX4jo0qbfyY4Or/G8D2eB9cvwBTnJkNAQ4EKoDHGm9395kEzR4DgUM7NjrJArVPXh+S1CjSR1W4LEtqFCnIzA4B/ht4xN2fTsQ5lNBS3wHh8mN33x2lzruN6orEy5hwuS6pUaQBM9ub4MoDICG/sNOVmeUD9wNbaPr+dlzo8TGpb+9w2dzM/Csb1RVpNzMbCFwUrv49iaGkJDP7GsHzEDsRXMEeTnCRcIu7P5HM2FLQT4BxwHnunrCnEiihpb7u4XKPNvsIO8NlYYJjkSxhZnnAQ0BP4OVENRGluSNo2PmjCvgB8KvkhJOazOxw4NvAP8P7sgmjJsfUVzumReMrpCP9nmBYyCrUIaRJ7v51dzegK7APcDtwE/CWme2VzNhShZkVAPcCOwh6aieUElrqKwmX3ZupU7utpJk6Iq1iZr8BLiEYJnKsuxcnOaSU5u673f0Td/8uQY/kycDvkhxWqvgpMBb4jrsn/D6smhxT3/JwObyZOkMb1RVpEzP7JXA1sJEgmX3awi7S0L3AL4BTzayTu1cmO6AkO4NgrONXzazx2Lzx4fJyMzsFWOLuX2/PyZTQUl9tt+l9zKwgSk/HgxrVFYmZmd0GfAfYDBzv7p8kOaR0tI3gXloeUEQwgDjb5RB0nolmZPjqFY8TSQpz91XAe0Bn4JzG281sOkEPq2KC6WREYmZmtwLfBbYSJLN5SQ4pXU0jSGbbgIT15ksX7j7C3a2pF0E3fgimvjJ337+951NCSw+3hMufmdno2kIz6w/cGa7eqrnjpC3M7GbgOoJfwse7u670ozCzo8zsgnAKp8bbjgD+HK7+WXOrdjxNTpwmzOxOgmmuyoCXqJ+cuAfwT+Bs/QcKmNkU6hM9BFPtFAKfEgzsBMDds35mFTM7jfo59eYAH0eputDdb+2YqFKXmV1EcJ9sG0HLSTHBd2sUwfcM4FmCibCjTYQggJndRzDsIW6TE+seWppw9yvMbBbwTYL26FxgIXAPcJeuzhroQTDzeWNjmijLdkURP08NX02ZCWR9QiP4HG4GjiLovXc4wdCaYoLB5w+5+z+TF1520xWaiIhkBN1DExGRjKCEJiIiGUEJTUREMoISmoiIZAQlNBERyQhKaCIikhGU0EREJCMooYkkkJl5G173hfvOCNdfTe67aBszu6iJ9xZt4HZrj7mtqc9KBDRTiEii3d9E2UDg8wRPIX+8ie2zEhpRx1tK/Xtq74S9jxA8UHM0wROjReoooYkkkLtf1LjMzGYQJLRNTW2P8A4wAdiViNg60KwW3merufsVUDenohKaNKCEJpKi3H0XwXydItIKuocmkqKi3UMzsxFh+XIzyzGz75jZx2a228xWm9mvzKxrWLe3md0e1i03s0/N7DvNnNPM7Dwze8HMNoX7rDSzP5rZiAS8x3wzu97M3jOzneH51pnZbDP7XzPLj/c5JXPpCk0kvT0CnAK8CiwheMDkNcAEM7sAeIvg8SazCGbWnwb80szy3f2nkQcys07AX4Ezgd0Ej5NZD+wLfB04y8w+5+5z4hG4meUQPGrlGGA7wUz224EBwDjg+8DvCGayF2mREppI+hpO8Hy8se6+FsDMhgLvAycQJIh5wIXuXhZuPxl4BrjezG4PmzVr3UyQzF4DLnD31bUbzOxK4LfAX81svLtXxSH+IwmS2XvANHcvjTifETyaZUccziNZQk2OIunt6tpkBuDuq4CHwtXhwOW1ySzc/izwIcFVW10XejMrAq4GdhI8nLIumYX7/Y7gamoUcGKcYh8QLl+PTGbh+dzd32iUcEWapYQmkr4qgf80Ub4kXM5x96a6yX8aLveKKDsaKABmuvuGKOebGS4PizXQKN4DqoFLzOwKMxvQ0g4izVFCE0lfxVGa/naGy9VNbIvcHtnhYmS4PDnagG/gtrBOv/aFHXD3pQT3+zoDdwDFZrbUzB40s7PNLDce55HsoXtoIumrpp3bI9Umj0UEHUma83YMx22Wu//WzB4DvkBwT+1I4Mvh6wMzm+7uuo8mraKEJiIAq8Ll/HgNgm4tdy8Gfh++MLPJwIPA/sD1wPc6Mh5JX2pyFBGAlwjuyR1nZr2SGYi7zwN+E65OTmYskl6U0EQEd19PcB+rF/CUmY1vXCccpP31eHXeMLNjzOwkM8trVJ4LnBSurojHuSQ7qMlRRGpdS9Dz8VzgIzP7APiMoPPIUIJ5JTuHy/VxON8k4NfAdjN7D1hHMPHwIcAgggHVP4vDeSRLKKGJCADuXgl80cweBi4GDiZIOiUEyeYR4EmC2fPj4WmCK8JpBLPnH07QA3Mlwf20u9x9Y5zOJVnA3D3ZMYhIBgpnxL8XuD/eHU0SeWxJX7pCE5FEOzLiQZw3ufvyth7IzO6k/nloIg0ooYlIoo0KXxBMNry8Hcc6H+jZ3oAkM6nJUUREMoK67YuISEZQQhMRkYyghCYiIhlBCU1ERDKCEpqIiGQEJTQREckI/x+5Jc99vBbgrAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"N=500\n",
"time = (m0-0.05)/dmdt\n",
"t = np.linspace(0,time,N)\n",
"dt = t[1] - t[0]\n",
"h = 300\n",
"\n",
"height = np.zeros([N,3])\n",
"height[0,0] = y0\n",
"height[0,1] = v0\n",
"height[0,2] = m0\n",
" \n",
"for i in range(N-1):\n",
" height[i+1] = rk2_step(height[i], lambda state: rocket(state, dmdt=root, u=250), dt)\n",
" pred_h = height[:,0]\n",
"\n",
"print('3c)')\n",
"plt.plot(t, pred_h)\n",
"plt.plot(t[-1], pred_h[-1],'*',color='red',label = 'Detonation',lw=25)\n",
"plt.ylabel('Height [m]')\n",
"plt.xlabel('Time [s]')\n",
"plt.title('Height v Time');"
]
},
{
"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
}