Skip to content
Permalink
master
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
{
"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",
"where $m_f$ and $m_0$ are the mass at beginning and end of flight, $u$ is the speed of the propellent, and $\\Delta v=v_{final}-v_{initial}$ is the change in speed of the rocket from beginning to end of flight. Equation 2.b only relates the final velocity to the change in mass and propellent speed. When you integrate Eqn 2.a, you will have to compare the velocity as a function of mass loss. \n",
"\n",
"Your first objective is to integrate a numerical model that converges to equation (2.b), the Tsiolkovsky equation. Next, you will add drag and gravity and compare the results _between equations (1) and (2)_. Finally, you will vary the mass change rate to achieve the desired detonation height. \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"1. Create a `simplerocket` function that returns the velocity, $v$, the acceleration, $a$, and the mass rate change $\\frac{dm}{dt}$, as a function of the $state = [position,~velocity,~mass] = [y,~v,~m]$ using eqn (2.a). Where the mass rate change $\\frac{dm}{dt}$ and the propellent speed $u$ are constants. The average velocity of gun powder propellent used in firework rockets is $u=250$ m/s [3,4]. \n",
"\n",
"$\\frac{d~state}{dt} = f(state)$\n",
"\n",
"$\\left[\\begin{array}{c} v\\\\a\\\\ \\frac{dm}{dt} \\end{array}\\right] = \\left[\\begin{array}{c} v\\\\ \\frac{u}{m}\\frac{dm}{dt} \\\\ \\frac{dm}{dt} \\end{array}\\right]$\n",
"\n",
"Use [two integration methods](../notebooks/03_Get_Oscillations.ipynb) to integrate the `simplerocket` function, one explicit method and one implicit method. Demonstrate that the solutions converge to equation (2.b) the Tsiolkovsky equation. Use an initial state of y=0 m, v=0 m/s, and m=0.25 kg. \n",
"\n",
"Integrate the function until mass, $m_{f}=0.05~kg$, using a mass rate change of $\\frac{dm}{dt}=0.05$ kg/s. \n",
"\n",
"_Hint: your integrated solution will have a current mass that you can use to create $\\frac{m_{f}}{m_{0}}$ by dividing state[2]/(initial mass), then your plot of velocity(t) vs mass(t)/mass(0) should match Tsiolkovsky's_"
]
},
{
"cell_type": "code",
"execution_count": 24,
"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": 25,
"metadata": {},
"outputs": [],
"source": [
"def simplerocket(state,dmdt=0.05, u=250):\n",
" \n",
" \n",
" dstate = np.array([state[1], (u*dmdt)/state[2], -dmdt])\n",
" return dstate"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"def rk2_step(state, rhs, dt):\n",
" \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"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
"def heun_step(state,rhs,dt,etol=0.000001,maxiters = 100):\n",
" \n",
" e=1\n",
" eps=np.finfo('float64').eps\n",
" next_state = state + rhs(state)*dt\n",
" \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",
" \n",
" return next_state"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The height of detonation: 599.26 m.\n"
]
}
],
"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(\"The height of detonation:\", height1, 'm.') \n",
" \n"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAbcAAAE0CAYAAACvs32dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOydd3hURdfAfycFUkiQUIJAIIgiTQQTmqD0DiKIEBABBRu+vCCvKPoCoqCoWF4sH0VFUHpTehEEpISmICpVEUNvoYWWNt8f9+5mk+xudtM2ZX7Pc5/dOzNn7rl3795zZ+bMGVFKodFoNBpNQcLL0wpoNBqNRpPdaOOm0Wg0mgKHNm4ajUajKXBo46bRaDSaAoc2bhqNRqMpcGjjptFoNJoChzZueQwRmSYiSkSSRKSCG3IlRSTelJ2Zjfq0M+tUIlI2u+rNKiLiZ6NXlKf1yW+IyPM210+JyDcuyu1JI9cwp3X1NHaulbL5j14SkV9E5CMRucfTuloQkXdNHQ96WhdPoY1b3sPykPECershFwX4mt+/zVaN8iEiMtf8c6/2tC75hG4iEuisgIjUAurkkj75AS/gDqAu8BLwm4j0zckDisgZ874ekZPHKQho45b32AT8Y35/0g05S9nTwLps1UhT0LkKBAJdMyhneXBfzVl18jQtgCBzKwnUB/4HJANFga9E5H7PqaexoI1bHkMZIWMs3Yq1XPmjmN0hDczdWUqppJzSL6+glLqllBJzm+tpffI5C81Phy9TImLbk7AgxzXKu9xUSsWZW6xSapdS6iVgpJnvAwz1oH4AKKVGmP+Nap7WxVNo45Y3sR3/cKX1ZlvGpbETjcYGSzd2SxG500GZlkB54DKwLFe0yl98BCSY35t4UhGNgTZueRCl1GFgh7nb23xrdkYf83OvUuo3ewVExFdEnhGRtSJy1nQ+OSciq0Skp4hIVnQWkebmONdxEbktIrEiEi0iw0UkwAV5fxEZLCLrzHGF2yJyWkR2iMhbIlI9TXm7DiWWwX+gp5nU1o4jwGSz7GZzf4ML+n1llj3hwu9hkdlnyixxoexss+yfdvI6i8j35rFvi8g1ETkqIhtFZJSIVHVFHyccBHYB3sATDspYXqDmA7edVWbeay1E5H+ms8VlEUkQkQsi8pOIDM3onhCR+81rflhEbojILfPe2iUiH4tIUwdyOX2t7KKUug38be6WcXJeYSLygogsE5EYU8cbIvKniEwXkQgHcnPN+zrUTBpv5762/R9k6FAiIt4i8pSI/GA+C+LNZ8OK7HgmeByllN7y4AYMApS5tXZSrolNuWEOylQG/rApZ2/7HvC3I9vOpkxZO/lewP9lUPdRoKqTc4gEYjKoY3saGT+bvCib9OczqEcBk82y/cz9ZKCyE/0CMMaZFPC2G7/hcFMmHijppFwx4LpZ9s00eV+4cD4fZOL+sr1OZYHB5vdf7ZQNBOLM/CZp7omGdsq/6oLOvwPlHejWF0jMQH63HbncuFbpztem3BGzzDEnZW5moF8Sdv7HwFwXzs32f/CumXbQgR4lgW0Z1LcaKObu9corm2655V3mYjwUwXnXpKXVlgTMTpspIiHARqAGcB5jPKAaEGJ+jsJ4E+8CfJoJPccAL5jfNwDNgdJAVWC0WXdlYI2IBNvRryqwHgjDeMCPxfDIK4nRDdbO1Ouyi/p8iTHYv8jcX0eKA4BlG2LmLcAwWoJh6BzR3ZQD+NpFPcD4PZIxvFh7OCnXFcOAQsp4KyLSERhoU9dDQAWMt/cIoBewGLjlhk6OmIthUGqLSO00ed0wDNzfwFYX6roFLAWeAhoClTBaM3WAEcAZoCY252pBREoDkzFakTsx7svKQCmgFtARmARcSCOXm9cqHSLih3GeAPudFD0MvAe0wbgGpYC7gPbAEoyXxQ9EpHkauX4Y9+A5c/8N0t/XC3EBs+dhEdDITJqM4fFZEuNF09JN3RaY4UqdeRJPW1e9Od4w/owKuAYE2MkvAsSaZVY6qGOqmX8JuMtBmY6kvK3dlybPYcsNwyAlmHlrAR87dXezkX/HTv5GM+86EOHkWvik2bfbcrPJt7zprs7gGk8xy/0NiIMyG8wymzLxG64zZbc6KbMa+63Tz8306By4t1K13My0Zeb+hDRlfzDT37JzTzhsyTg5dkXznlZAozR5j5vpt4FgN+rMrWtl93yB123KdMzCsSaadaxxkH/GzB+RQT0OW24Y04Ysuo52IP+JTZk22X1Nc2PTLbe8jcU5pBj23bQ7ASXSlLUiIneQ4r49Sil11N5BlFIrMLoowL25df0wvMMABiulEu3UvRjD8AEMsO3HN1sITc3dcUqpnx0dyF7d2cSX5mc40CxtpoiEk6KjO602C5bWyYMiUtlO/aFAqzRlLViu7clMHDczWO4h6ziviJTDcH+HbJo/qZSKwXhhAGidJttyzlcxDKCr5Na18heRYuZWQkQiReRj4C0z/13z/5RZLC2lpiJSJGuqOsTSwo0B3nFQZgTGC7Ft+XyFNm55mxWkdL/Y65q0pF3F6NJIy0MYc28AfrL5U6bbgF/NcpFu6GfxCtunlDrkpJzFdbwMRnelhVY236e7cdxsQym1C7A44fS3U6Q/RrflNTLnAr8IY5wF7DtrRGF0wSUC89Lk7TE/HxWRf0kGk6yzgWXAFaAchnckGN3eXhityiOuVmTeV/82nRVOmw4hVucHoLNZ9N40onvNz1LAFHHsvZmW3LpWP2LcC9cwek12YXT1XwdaKqVey6gCEWkoIlNF5HcRuSoiyTbXxfKCV5SUbs5sQ0S8MbqKAb539NKolLoBrDR3H8puPXIDbdzyMEqpBFIeeK3EJvyVOZbWwdxdqJS6mVae1A+OX0n5U9rbLONmpd1Q0ZUxBjCcWdLKAFQxP88ppU67cdzs5ivzs7uIWMbWMFuZlrG4+Uqp6+5WrJS6RsqLRx87RSxpa5RS59PkTQP2YRi/T4ELIvKjGN6jrUTEl2xEKXWLFANuafFbXqBcbrWJ4dn6B0YXWysMh5WiDooXT6PDAYwxIIBngJOmx+UnItLd7I2wR65eKzsEA/8TkVLOConIBCAa49xqYoyVOfJKLO4gPSuUwhg/Bdf/t2VzsBWZY2jjlvexdBV5YwyKW+iJMeZmWyYtmflz+LlR1mII4jIoZ9u9FGTzPdhOvieYiTHGE4Ax5mOhGUZ3JWSuS9KCxTDcKyLWlrHpTBOZpowVpVQ88DDG+MlZjN+mOYYT0A/AGdO9PTsf3JZ7qauINMFw4ognfavSLqYu32OMq13F6K57CGN89g5SnB8WmyI+dqp5EWOcaz/Gg78uhjfnAuCs6TKfyt0+F69VI2UGD8D4fzXG6GEBuA8nLwEi0g942dzdgPF/roHxQmm5LrY9J/auTVax/f9l9n+bL9DGLY+jlNqJMQ8JUndNWr7/A/zkQNxy8yYDRVVKRA9nmzsRDSw3f7EMytnmX7Pz3aN/HKXURYwHMhgefqT5fkgp5YqXoCPWkuLlZtt6s3y/huFdaE+3K2ZX153A/cBzwByM7sMQDOORnRP3twDHMN7uLeM/K81r5AqtSel67qKUekMptUUpdcI8lzilVBxO7hmlVLJSaopSqiZGS78XxnSTGIwXun7AVrM73VYuV6+VUuqqUmob8AgpBq6diDjyjH3R/NyA0YU5Vyl1QCl1wea65HQLyfb/l9n/bb5AG7f8geVtsK6I1BCRKqS48X6rTPcmO1gcSLyAtO7d2cEx87NGBuVq2pEBsExYLuPG2EpOYemabCIiVczuycfMtKy02izOMJYQYVHmuAekOO8sctCtbFuHUkrtU0pNVUr1xnBztxjEKBHJljBL5r1kcWy5y/x0x5HEElj5jFJqo5NytVzUJ8Y0Ai9iTAmwhLm6G/vdvLl2rWyOlww8izHuBjDW5je2xRJKb56T/+x92ambHS6Qoqer/9vTZss4X6GNW/7gW4zWFxgttj5p8hyxAcNRAVK3SLKLLeZnbXG+3Ed38/OcMqKvWLAN8OxsnllmsIRCsveQscc6UgJW98fo9g3AmD+YHW/7FoMRijF+2oiUMUe3lygy3/LftUmq7qhsJrA930vAcjdkLWNrDq+7iLTEcFpxC9OIjCfFQcelc87ha2U5xikM93kwWq49bfNN71NLl6ize9KuwbbB3fs6FcqIOxtt7nZxYIQREX+MKUKQ8j/PV2jjlg9QSh3HWC0AjLd9yx9gRxpjkVbuAikPqudEpLOjsmBMHTBd013lG1KM5yf2/igi0gVjXhSkuN1b9NuHMc8NYKSI1HWim7vjD5ZuNJceouabtKWF1g942vy+OjucXUyvTEv3ch9SPCdPkuIWnwoXWhhVbL672m2YIaZXZFUMI3C/m2/tlt6C0qYBT4WIlMBJsAARuSuDcbHygL/53XrOnrpWafiQlO67122nvZiG2fLy1MWesIg8T8aeiW7d1w6w9FJUwnD5t8c7pEwz+iILx/IcOTWBTm/Zu2G0JtKGxxnkglwI8BcpoX2mYAy8lzHzqmK0rL7GcADolEY+o/Bbb9nk/4AxJ6wkxsPkv6SEG/obO5NyzeNbQltdw4h4Uhvjj3Unhkv6R8CqNHIZTeJ+0sxLxvBMK40xQO8DeDm4VmHmNbK9xt2y8TccadYZhxEtRgHvOym/HcPL9XWMh96d5nWphuF+brluxwBfN3VJN4nbDVmHk7gxvPEsobpOYhjxMFP3nsAhjNbHIexMssdoYZ0CPsbwBg7HcESpjDFt4qAplwDU8MC1cjppHXjb0b1DysRqhfFiaIkKcj+GwU8idZg8e6HNvjLzzpr/jWCb+1rsHMveJG4vUgITJAOfYfznQkydptvosDi77v/c3jyugJ0L7w+8gjF/5DJwA+PBuABo7ECmN7AZY+A4DtiNMXhr9yGWVTkPXRfb+IOWCA4hLsqGYQRiTmsc7W1t0sjmRmzJehgPQmd1uBRbMs31chSvcrITXVbblDuPmw/CDH6HcPNhYqtLbSflt7vwe50FIjOhS44YNzO/H+lfEixbAoajh90IMqQ2AI62eGCAh65VRsYthBRD+nOavCCMeXyO9NsLPJjBtY0gJSpQ2k3HlrQ9R08rkOaCVyYl+OhZjPlB8zFizMUDI+3IWMLu3MQYG/jO5uZaDHg7OFam5Dx8fb61ufHceqPCMELdMeLPxZjnfRvjLXk9hotyuuDBZGDcbMo1x3AXP2HWewmjb384dkKH2ZEvZuqwGaPrJR7D4EVjtObuTVPeqXEzy4RhhB/7EyOmoCvG7XGbch/nwG/4k039+1z4PzxnXtffMJwBEjAmD2/DaBmXyKQeOWbczDJNMSYBx5r3QwxGzMcGZr4j4xaCEYdzCsbL5inzXrhmXoNP0t4LuXytMgw3BoyzKd/Bzn0+DqPlavmf7MZ4offDaGlmdG0bYzyvLNfGbeNmlvHGGIv/AeNFLh7jubsSo5VsNxxdftnEPEmPY0YU+BWjO2ssMFYZk5gt+SUxIqsftkl7DONhfQZ4WJkRFMxxow0YYwZDlVIT0xwrU3Kago+IdCJlvbLaysESQhqNJm+Tl4zbeIzBzW+UUv1clNmN0Uzvp5T6Jk1eUwxnhTMYS2skZ1VOU/ARkUUYwZ53K6XqeVofjUaTOfKEcTNDu5zC6AeuoYwQPBnJVACOYzSl71B25gmJyAkM76rGyphsmWk5TcHHDGx8GGNw/lmlVP70EtNoNDkS3iUzRGAYtuNKqQMi8iBGxPuSGC2o1Uqp6DQyFrfxP+wZKJNdGEaqLilR7zMrpymAmPOPvDG6w7/A+E+cIZsi4Gs0Gs+QV4ybZVb+ERGZTvoJvaPN7qInbQySZfmQf3BMTJqyWZHTFExmk2bCLcZ4a44saqnRaHKHvGLcQszPhzHeoj/AiAx+0Uz7P4xQSFdJmVxriXvmLFK7JbaibezCzMpZEZFnMcLtEBgYGFGtWrZG89HkIiVKlODSpUt4eXnh5+dH2bJlKVGixNzIyMi5GUtrNJrM8vPPP19QSrmzColb5BXjZomU4gN8qZQabpO3VEROYUwH6Cci45Sx6KZl9r+7g4aZlbOilJqK4WJOZGSk2r17d2ar0mg0mkKJiDjrPcsyeSX8lm3E6XSD+Eqp3RiL+HmRslqyKxHpLXn2ItG7K6fRaDSafEJeMW7HbL7/7aCMJd2yYKdFppKTesPs1J9ZOY1Go9HkE/KKcfvF5ntJB2UsK9xaxsMsy8rXNCNY26NemrJZkdNoNBpNPiFPGDel1EmM2IdgBANNhRlJ/AFzd7cpcxzDKBYh9erJFpmmGOs4nSFliYdMy2k0Go0m/5AnjJvJ2+bnaBGxLHiIiPgBkzCWdP+Z1AZnvPn5nojcbSNTBsPDEuBdO1FGMiun0Wg0mnxAXvGWRCm1TEQ+wAieu0NEdmBMBaiPsXbRSaCXsgmpopRaKCKTgBeA30RkHUawVMtSEN9jLOeQ9liZktNoNJr8zJK9J4m9Hs+jdcpTIrCIp9XJUfJE+C1bRKQrMBgjOkgAxoTqpRgtqfMOZHpjLFVzH8Y8uYPANGCSs9ZXZuVscWUqwNWrVzl37hwJCQlOy2k0uYGvry9lypQhODjY06pocpHkpCT6fvIBW87WoIi3D9P616PJPaUyFswhRORnpVRkTtWfZ1puFpRS32Es5+COzGyMSBPuHitTcu5w9epVzp49S/ny5fH398dmcV6NJtdRSnHz5k1OnjwJoA1cIeK3U9c4dSWBRiVmgypGnYqtPa1SjpLnjFtB49y5c5QvX56AgABPq6LRICIEBARQvnx5Tp06pY1bIWL+7uMcvVWXo7fq0r1uWYoVLdiP/7zkUFIgSUhIwN/f0YwDjcYz+Pv7627yQsTN+CSW7j1l3e9Rv+CHzdXGLRfQXZGavIa+JwsXa/44w7XbiQCElwygXngJD2uU82jjptFoNAWc3Vs/oiSXAHg8MqxQvNwU7E5XjUajKeT8sn8TS4LWUuqeVTx0NZhu96/ztEq5gm65aTLF9OnTERHi4uIcllmwYAGPPPII5cuXp1ixYkRERDBnzpxc1FKj0czd/j4AF3y8uRCQyJ0hzmLGFxy0cdPkGB999BHFihXj448/ZunSpTRv3pzevXvz6aefelo1jaZQcPvWDe66fIhSiUkAPFy2vYc1yj10t6TGbZKSkoiPj8+w3LJlyyhVKmWSaIsWLTh16hQfffQRgwcPzkkVNRoN8Nu6mTx/9TwDrsKSgDJ0inrN0yrlGrrlpsmQ/v37ExkZyffff0/NmjXx8/Pj7Nmz6cpNmDABPz8/li5dCpDKsFmoW7cu586dy3GdNRoNBP46HQBfoHzpbvgVLTzzbXXLTeMSx44d45VXXmH06NGEhoayY8eOVPljx45l/PjxLFmyhLZt2zqsZ9u2bdSoUSOn1dVoCj1//7GD6gl/AJCgvLmn7SAPa5S7aOPmAcJHrPC0Chx7t6Nb5S9evMi6deuoU8dYsMESvgng9ddf59NPP2XlypU0a9bMYR3r169nyZIlTJs2LVM6azQa1znz4+dYpmrvC3qIiHLO1mcueGjjpnGJ8uXLWw2bLcOGDWP+/PmsWbOGBx980KH8sWPH6N27N126dKF///45qKlGozl78SSj/XfTJuQOelyNw6/Rs55WKdfRxk3jEqGhoXbTFy1aREREBPXr13coGxsbS/v27alYsSIzZ87MKRU1Go3JtNX/5YyvN98UD2ajfxDLGrbztEq5jjZuHsDdLsG8gKOIBsuXL6dTp0707duXmTNn4uWV2kfpxo0bdOrUifj4eFasWEFgYGBuqKvRFFqSkhW/XfsNzJC2jQPq4eXt7VmlPID2ltRkifvuu49Vq1axfPlynn/++VR5iYmJPP744xw5coRVq1ZRpkwZD2mp0RQe1h04S/Sx0dx94mHuu16UZzq842mVPIJuuWmyTP369Vm+fDnt2rUjODiYDz74AIBBgwaxcuVKJk6cSGxsLNu3b7fK1K1bl6JFi3pKZY2mwPLV5r9R+LDnWgdeiBhM6RJ3elolj6CNmyZbePjhh1m8eDFdunQhKCiIN954g7Vr1wIwZMiQdOX//vtvwsPDc1lLjaZg8+vxy+w8FguAj5fQr1G4ZxXyINq4aTJk+vTp6dL69++fzuuxXbt23L5927p/7NixnFVMo9Gk4sstf1u/d76/HGWL+3lQG8+ix9w0Go2mALDvyHYSLvanhf8KfEhkQJOCvyCpM3TLTaPRaAoA0zaNYmtQMgRtpkvcr9Qq38XTKnkU3XLTaDSafE7MyUNs9T5l3a8d3sGD2uQNtHHTaDSafM7Jlf/Ht6fP0ibuOlVuC0+2G+FplTyO7pbUaDSafMyVSxeodWIeQZLAh+cvsr3+a4Vy0nZadMtNo9Fo8jEHlnxEkNwEIMarPPXaPOVhjfIG2rhpNBpNPuVG3BXuPfatdf9s7UF4++gOOdDGTaPRaPItXyx+iT0BCSjgNKWp0+EZT6uUZ3DLxIvIHUAzoC4QCtwBXALOAb8Am5RSl7NZR41Go9Gk4fylUyxI3MGXoaWpdfs2TwZH0aGIDmlnIcOWm4h4i8jjIrIBuAAsAkYCzwI9gOfM/cXABRH5UUS6i4ge0SwgjBkzhlKlStnN69+/P5GRkbmsUWrCw8PtRlFxlWPHjiEiiAhbt25Nlz927FhExO1wYfHx8YwZM4a9e/faPd7y5cszrXNGREZG6nXzCjifLh3CFW/jEX7ey5eHOg/zsEZ5C6ctNxHpBbwLVAAEw7htB/YDscBVIBgoCdQAGmK07JoCx0VkhFJqbk4pr9FkJ8WKFWPOnDk0btw4Vfq8efMoVqyY2/XFx8fz5ptvEh4ebnehV40ms1y5mcDGfzpQN/g8fxS/QIegpgQF3uFptfIUDo2biGzFMFbngYnADKXUrxlVKCJ1gP5AL2CWiAxWSjV2LqXReJ7OnTuzcOFCJk6ciLfpSv3bb79x4MABevToQXR0tIc11GgMvtp8lJgbZYi5MZz7E84y6N9RnlYpz+GsW7IKMAyoqJQa5ophA1BK7VVKDQXCgP+Y9WgKETExMURFRRESEkJAQABt27bl0KFD1vyNGzciIvz++++p5Jo1a0b37t2t+5Yuzx9++IHatWsTGBhIkyZN+OOPP5wef+nSpURERBAYGEiJEiVo0KABmzZtylDvLl26cO3aNTZs2GBNmzt3Lk2aNKF8+fLpysfGxvLcc88RGhqKn58fDz74IDt27LDmBwUFAfDUU09Zuz1tg0nfuHGD5557juLFi1OhQgXeeOMNkpOTUx3jxx9/pEGDBvj5+REaGsqgQYOIi4tLVeb333+ncePG+Pn5Ub16dZYuXZrhuWryL7HX4/nKJkDy063b4lc0wIMa5U2cGjel1ESlVHxmKlZKxSul/gfclTnVNHmNxMTEdJtSKlWZ2NhYmjRpwqFDh5g8eTLz58/n+vXrtGrVips3b7p9zJiYGIYPH85///tf5syZw7lz5+jRo0eq4x47dsw6vvTXX3/RvXt3WrRowbJly5g1axadOnUiNjY2w2MFBgbSqVMn5syZY02bO3cuvXr1Slf29u3btGrVih9++IEJEybw/fffU7p0aVq1asWZM2cAwzABjBw5kujoaKKjo7nzzpS1tV555RWKFSvGwoUL6dOnD2+99RYLFy605u/fv5927dpRqlQpFi1axJtvvsns2bNTvQDcvHmTtm3bEhcXx+zZsxk5ciRDhw4lJibGxSusyW9M2fQX1+OTAKgaWoxOtct5WKO8icNuSaXU9ew4gFLqRnbUU+DYMB42veta2Qf6wSOfpE5b+m/4ZYZr8k1HQPPX3NMvDRcvXsTX19duXkREhPX7xx9/zPXr19m7dy8hISEANG7cmPDwcKZNm8aLL77o1nFjY2PZunUr99xzDwDJycl07dqVQ4cOUa1atXTl9+zZQ1BQEBMmTLCmdejgepy9qKgoBgwYwKRJk9i7dy8xMTF0796dd99N/VvNnDmT33//nT/++MOqW6tWrbj33nv58MMPmTBhAvXq1QOgSpUqNGzYMN2xHn74YT788EMAWrduzerVq1m8eDE9evQA4K233qJSpUosXbrU2k0aEhJCz549iY6OplGjRnz99decO3eOHTt2UKFCBcBwsGnSpInL56zJP/z+5w5OH3yNUOnDWRXKsNZV8fYST6uVJ9Hz3DQuUbx4cXbt2pVu69SpU6py69ato3Xr1gQHB1tbd0FBQURERLB79263jxseHm41HgA1atQA4MSJE3bL33fffVy5coV+/fqxdu1arl937x2tQ4cOJCUlsWbNGubOnUvLli3teoquW7eOiIgIKleubD1PgKZNm7p8nm3atEm1X6NGjVTntXPnTrp27Wo1bACPPfYYPj4+bNmyxVomIiLCatjAeJkoU6aM6yetyTd89uN/2FDyPMWrTKB32YW0rVnW0yrlWVye5yYiIUA4cEwpFWuTficwHrgfOAaMcXV8TpN/8PHxsevyX7JkSU6fPm3dv3DhAtu3b2fevHnpyrZs2dLt495xR2oPsCJFigBw69Ytu+XvvfdelixZwrvvvkuHDh3w9fWla9euTJw4kdKlS2d4vKJFi/Loo48ye/ZsNm/ezLhx4+yWs5ynvdZslSquDTPbOzfb8zp9+jShoaGpynh7e1OyZElrN+uZM2fsGjJt3AoeP2yfx7YilwHhtK8PtWvej4hutTnCnUncr2E4mDyAMQ0AESkCbMEweoJh4JqKyH1KqZPZq2oBo/lrWesqfOST9F2VeYCQkBAeeeQRRo0alS7P4mDh52esDhwfn3o4NzY21uF8Onfo2LEjHTt25MqVK6xYsYKhQ4cyePBg5s51bVZKVFQUnTp1shpGe4SEhBAZGcmkSZPS5RUtmj0Tae+8807OnTuXKi0pKYmLFy9au3zLli3LwYMH08mmldPkb1RyMmV+/IzXil5iUoniVEoMJKr1UE+rladxx7g1B/5O0yrrCVQGNmLMh+sMvAj8C8MYagoZLVu2ZP78+dSsWRN/f3+7ZSxdaAcOHOCBBx4A4Pjx4xw6dIiqVatmmy7Fixend+/ebNq0yS03/tatW/PYY49RrVo1ihcvbrdMy5YtWRRyG8oAACAASURBVLt2LRUrVnTYSsqolZkRDRo04LvvvuOdd96xdk0uXryYxMRE65havXr1mDVrFidOnLBe161bt2rjVsDYt3E+98fv4/54aHftJke7zfe0Snked4xbBWBvmrROgAIGKqWOAmtFpAPQHm3cCiXDhg1j5syZtGjRgsGDB1O+fHnOnj3Lpk2baNKkCb169aJChQrUq1ePUaNGERAQQHJyMu+88461NZIVpkyZQnR0NO3ataNcuXIcOXKEBQsW0LdvX5fr8PHxYf585w+Pvn37MnnyZJo1a8bLL7/MXXfdxcWLF9m5cydly5blpZdeokiRIlSuXJn58+dTq1Yt/Pz8qF27tst6jBw5krp16/Loo4/ywgsvcOLECV599VXatm1Lo0aNAGOawbhx4+jYsSNjxozh5s2bjBo1KltawJq8QWJCPHdsSekeP1yyMw1qt/KgRvkDdxxKSmBEKLGlEXDYNGwW9mDMcdMUQkqVKsX27dupVq0aL730Em3atOGVV17hypUrqR7ss2fPpmLFivTp04fXX3+d0aNHc++992b5+LVr1+b8+fMMGzaMNm3aMG7cOJ555hnee++9LNdti5+fHxs2bKB169a88cYbtGnThiFDhnDkyBHq169vLTd58mQuXLhAq1atqFevHqdOnXJSa2pq1qzJqlWrOHfuHN26dWPkyJH06tUr1XSBgIAA1qxZQ2BgIFFRUbz55pt8+OGHVKpUKVvPV+M5fv7uf1RKPg5AnPLn7h5ve1ij/IGknafksKDIZSBaKdXe3A8D/gGmKaUG2pSbBXRRSrkfrygfEhkZqZx5xx04cIDq1avnokYajWvoezPv88+pw3y5qCNDLl2gVHIy0ZVfpFG/dzytVrYgIj8rpXIsMK07LbeDQBPTaxKgN0aX5E9pylUAzmaDbhqNRlOombDiOb4P9qNzWDlmFCtL3cdf97RK+QZ3jNu3QCCwU0TmA28BccASSwERKYrhTXnIbg0ajUajcYlNv6zlJ9/zAMR5eZFQNQq/gELRIZYtuONQMgljjK03Rkit68AzSqkrNmU6YxjAjAP5aTQajcYuycmKj6MDqXypBbdDf6RUclGe7jja02rlK1w2bkqpZKCPiIzCWKh0v1LqappiR4HHgfSLYmk0Go3GJRb+coK9xy8DbQm82Zz3nroLL2+9RKY7OFvyJgpYoZS6ZpuulPob+NuejFLqF4wVuTUajUaTCa7cSOC9VSkT8/s/XJ37704fR1XjHGdjbrOBcyKyXEQGiEjGsYs0Go1GkyXeX76Hi9eN6D3livvxYvO7PaxR/sSZcXsVY9J2e2AqcEpENojIv0WkYq5op9FoNIWIBes+48ilfjT1WwMoRneuQUARd1wjNBYcGjel1ASlVCMM1/5/Y7j8Nwb+B/wtIrtF5DUR0e1ljUajySJX4mKZdmwyB/x82Bf+I/0qzaZdrTszFtTYJcOpAEqp00qpz5VSLTEcSZ4GVgA1gLeBP0TkgIiME5Ecm5Cn0Wg0BZnv5r3EWbORVkQperQZ5lmF8jlureemlLqklJqulHoEKA1EAQuAcsDrwA4R+UdEPhaRpqLXY9BoNJoM+Xv/Lp6IWcbCk6d54NYtuhZtQM0quq2QFTK9WKlS6rpSar5SKgrD0HUGpgP+wBDgR+C/2aGkxnOISIbbxo0bM6xnxIgRqRbUdIVbt24hInz55ZfWtIYNG9KnTx93TyPbsaebRpMZkhITiV/8L3wlibsSEnn1Qkle7jHF02rle7JlpFIpFY/RVblCRLyApkBXQK+7kc+xXSrm5s2btGjRgpEjR9KxY0drumV1bGe8+OKL9OzZM0d01GjyMztnv0mjRMP1P155E9jtc3x80i+Cq3GPbHfDMSd7bzA3TT6nYcOG1u9xcXGAsdK0bborhIWFERamF4vQaGyJ/mU5VY+mLHj7c/gzNKoe4UGNCg6Z6pYUkbIi8oCIPOhoy25FNXmbixcv0r9/f+688078/PyoVKkSL774ojXfXrfkkSNH6Ny5M0FBQQQHB9O1a1f+/ttufACHXLp0ifr16xMZGUlsbCxgGOFBgwZRpkwZ/P39adCgARs2pLxrvfrqq1SqVIm0K2IsXLgQLy8vjh83lhdZtGgRdevWJSAggJCQEBo1asS2bdsc6rJnzx5Kly7NwIEDiY+Pp3Tp0naX2mnQoAG9e/d26zw1BY8bt67z7s+vERVWmmi/ohz2qUq9PmM9rVaBwS3jJiKPi8gB4CSwC9jsYEu7UoCmgDN48GB2797NJ598wpo1axg3blw642GLpYvz6NGjTJs2ja+++or9+/fTrFkzrly54lDOlvPnz9O8eXN8fHxYv369dbHTfv36MWvWLMaMGcOiRYsoU6YMbdu2ZefOnQBERUURExPD9u3bU9U3f/58HnzwQcLCwti/fz9RUVG0b9+eFStW8O2339KuXTsuXbpkV5edO3fSokULevbsyRdffEGRIkXo06cP06dPT1XuwIED7Ny5k6eeesqlc9QUXMbOeYKjReCMjw+DQ8tws8sEfHyLeFqtgoNSyqUNwzMyCUgGLmEsSurIuG12td78vkVERChn7N+/327653s+V7Wm11K1ptdSn+/5PF3++zvft+ZP/316uvw3tr5hzZ9/aH66/OGbhlvzl/+13KmOrnLt2jUFqK+//jpdXpUqVdTUqVMdyr766quqfPny1v2PP/5Y+fr6qpiYGGvaX3/9pby9vdVHH32klFLq5s2bClBffPGFtUyDBg3UE088oU6dOqVq1KihmjZtqq5du2bN37NnjwLU3LlzrWmJiYnq7rvvVo888og1rWrVqmrIkCHW/bi4OBUQEKA+/fRTpZRS3377rSpXrpzD87HVbfPmzSooKEi9/PLLqcr89ttvClDbtm2zpg0fPlyFhYWppKQkh3XnFo7uTU3Os+/4ZfX4OwNU469qqFrTa6k3v+ntaZVyHWC3ysFnszstN8tCQkOA0kqpukqphxxtWTW6IvKOiChze9lJud4isllErohInDm5/EXTscVZ/ZmS09inTp06jB8/nsmTJ/Pnn39mWH7nzp00bNgw1TjcXXfdRb169diyZYtT2ZMnT9K0aVPCwsJYtWoVxYqlLAOyc+dOvL296datmzXN29ub7t27p6q3Z8+eLFiwgOTkZACWLVvGrVu36N69O2Cs6H369GkGDhzIunXruHHjhl1dNmzYQNu2bRkyZAgTJkxIlVerVi3q169vbb0lJSUxc+ZM+vXrh5eXvs0KK7cSkhg2fy87r3Tl1l//5qHrFRjRc5qn1SpwuPMPuwfYppT6VCmVmFMKAYhIPeAVjMVQnZX7HJgFRGK0GH8AqgKfAQtFxG4Y7czKaRwzdepU2rVrx+jRo7nnnnuoVq0aixcvdlj+9OnThIaGpksPDQ21jp05Yt++fRw5coR+/frh7++frt4SJUrg65va2yw0NDRVl2JUVBSnTp2yGrx58+bRrFkzypYtCxjGbfHixRw4cIC2bdtSqlQp+vbtm0631atX4+XlxZNPPmlX1wEDBjBv3jxu3rzJ6tWrOXPmDP3793d6fpqCzdsrDnDknOGcdcO7Eq/2XEiRIkU9rFUBxNUmHsY42+ycbEaaxykK/GEe7zsMA/eynXKPmXmngXts0kOB/WbekOySc7RltlsyP+KsW9JCcnKy2rNnj+rRo4fy9vZWR44cUUql75bs1auXevjhh9PJN2zYUHXr1k0p5bxbcuTIkapIkSJq9erVqeSnTJmivL29VXx8fKr0ESNGqJCQkFRptWrVUoMGDVJXrlxRfn5+asqUKXbP6dKlS+qbb75RISEhql+/fql0++yzz1TLli1VWFiY+ueff9LJXr16VQUGBqpZs2ap7t27q6ZNmzq4crlPQbo38wurfzulKr263LrN2p7+nikskIe6JdcC9TJnQt3iLYzQXs8DzjwLXjM/X1VKHbEkKqXOAi+YuyPsdDNmVk7jAiJCnTp1ePfdd0lKSuLw4cN2yzVo0IDo6GhOnjxpTTt27Bi7du2iSZMmGR5n7NixvPDCC3Tr1o3Nmzdb0+vXr09SUhLfffedNS0pKYlFixalqzcqKoqFCxeyePFiEhMTeeyxx+we64477uDJJ5+kU6dO7N+/P1Ve0aJFWbJkCWFhYbRq1YozZ86kyg8KCuLxxx9n4sSJLF26VDuSFGL2HdnO3J/aUNd3FwDta5WlV309PSbHcNUKAhWBs8B7gHdOWFqgAZAIzDL3p2On5YYRzFkBtwF/B3WdMMs8mFU5Z5tuuRnUr19fffTRR2rNmjVq9erVqkuXLio4OFidOXNGKZW+5Xbjxg1VoUIFVatWLbVgwQI1f/58Va1aNVWxYkV15coVpZTzlptSRitxwIABKjg4WO3evdtaplu3bqp48eJq0qRJauXKlapz587K19dX7dy5M5XOR44cUYC68847Vbt27VLlTZw4UQ0YMEDNnTtXbdq0SU2ZMkUFBwerV1991a5uly9fVnXr1lX33XefunjxYqq6Nm/erAAVFBSk4uLi3LrmOUlBujfzOrdv31KPT7lf1ZpeSzX6qoZ66f0+6vL1+IwFCzDkcMvNnZW4Y0SkCbAE6CYi601DkOyg/Duu1g0gIn7ADCAWw2nFGXXNzz+UUjcdlNkFlDfLWiYnZVZOkwGNGjXiq6++4tixY/j6+vLAAw+wZs0au+NqAP7+/vz444+89NJL9O/fHxGhZcuWfPzxxwQHB7t0TBFh6tSpXL9+nbZt2/LTTz9Ro0YNZsyYwfDhwxk1ahTXrl3j/vvvZ/Xq1dSrl7rj4e677yYiIoKff/6Z8ePHp8qrU6cOq1atYujQoVy6dIly5crxr3/9izFjxtjVpXjx4qxdu5amTZvSvn171q1bR1BQEABNmjShZMmSPProowQGBrp0bpqCxdxZ/+JwkURAuO4lPBT5IMUDdBSSHMVVKwgIMBFIwDBoyRhTA9JuyUCSu1YW+BCjxdTTJm069ltu/zbTv3NS30SzzAdZlXO2FaaWmyZz/PzzzwpQW7Zs8bQqqdD3Zu7wR/QqlTi6uPrlnVKq9ZfV1GvTHvW0SnkC8krLDRgBDMboNlwB/AnEuWNIHWFGNBkKfK+UmueCiMX3+7qTMhbdgrJBLhUi8izwLEDFinrdVo19zp8/z+HDh3nttdeIiIigcePGnlZJk8tcOBNDmdXP4S2KurfjeetUEHVemeVptQoF7hi3AcANoIlSam92KSAi/sDXwFVgkKti5qfTqQLZKJcKpdRUjNXJiYyMzFJdmoLLokWLGDRoEDVr1mTWLP1AK2wkJsRzblpvanAZgEsEsy6xBQ2LBnhYs8KBOx6B5YGfstOwmbyDMcdsmFLqtIsy18zPYk7KWPKu2aRlVk6jcZvnn3+e5ORkfvvtN2rXru1pdTS5zKTpTxCe8DsAyUo40eIzEsXPw1oVHtwxbicxWm7ZTVeMcbp+IrLRdgPamWVeMNMsi2cdMz8rOanX4mN7zCYts3IajUbjMlOXjGJqkYM8US6UYz4+7Kj8PPc93MXTahUq3OmWnAc8KyKBSilnY1aZwbIGnCPuMrc7zP095mdNEfFX9j0f66UpmxU5jUajcYm9f+5hWuxi8PLizyJFeLvUXUx58m1Pq1XocKflNhY4DCwVkSrZpYBSKlwpJfY2jKkBAMPNtDqmzHHgF6AI8HjaOkWkKcactjOAdbXNzMppNBqNK1y9lcDLS69S+dx9FElWhCYkM/rReXh564h+uY07LbelGJ6SzYEDInIUx/PclFKqbTbo54zxwALgPRHZppT6E0BEygD/Z5Z5VxmLp2aHnEaj0TgkKVkxZM4e/jwXx588SfX43QxtWZWwsnd5WrVCiTvGrVUauarmZo8c9yBUSi0UkUkYIbN+E5F1GHPwWgLBwPcYgZCzRU6j0Wic8f7qg2w4dN66/3zngbStU96DGhVu3DFurXNMi0yilBokIluAFzHG7LyBg8A0YJKj1ldm5TQajcYes5Z/yqwtoYARgeaFZlXoog2bR3F5zE0ptd6dLTuUU0r1N8faPnBSZrZSqrFSKlgpFaiUilBKfZ6RgcqsnMZ1Tp8+TYcOHShevDgiwsaNGzNd1xdffME999yDn58fERERrF/v2i22detWGjRogL+/P5UrV+aTTz5JV0ZE0m0NGzbMtK6awsWyn6bx0fkp1K/4FhW9Y2hVvQzD29zrabUKPe603DQat3j77bf59ddfmTNnDiEhIdSoUSNT9cydO5fnn3+eMWPG0KRJE77++ms6derErl27qFWrlkO5P//8k7Zt29KpUyfGjx/Pzp07GTZsGAEBAQwcODBV2f/85z/WhUoBa1xIjcYZh47+wgdHPiTex4tdgdCg/FQ+7jEALy/JWFiTo2jjpskxDh48SIMGDejQoUOW6nnjjTfo168fo0aNAqBp06bs2bOHd999l5kzZzqUmzBhAuXKlWPmzJn4+PjQokULYmJiePPNNxkwYAAiKQ+g8PBw3VrTuMXVyxfxmdmfrsWv8dUdxQlOSua5ph8Q5F/E06ppcNItKSI/mTEfM42INBaRn7JSh8bz9O/fn8jISFasWEGNGjUICAigY8eOxMbG8ueff9K8eXMCAwOJjIxk3759gNHVt379er777jtEhPDw8Ewd++jRoxw+fJgePXpY07y8vHj88cdZtWqVU9lVq1bRrVs3fHxS3uGioqI4ceIEv//+e6b00WgA4m/f4p9Jj1El+R+GXrrC6+cvMbjCQOrVbOlp1TQmzsbcqgGbReQHEekpIi6tgy4iRUWkl+mF+BOOPSo1+YiYmBhGjx7NuHHjmDp1Ktu2bePZZ58lKirKuuhnYmIiUVFRKKWIjo6mbt26NG/enOjoaOvioUopEhMTM9wsHDx4EIBq1aql0qd69erExsZy/vx57HH9+nWOHz9uV862XgtjxozBx8eHUqVK8fTTTxMbG5u1C6YpsKjkZH79v37cdzslzsPdd79GVJuXPKiVJi3OuiXvAd7ECGbcArgmIlsxJjcfAC5iBDsOBkpirJ7dCGiMEZ8xEWP5mDE5pLsmF4mNjSU6OpoqVYz5+/v27WPChAnMmDGDvn37Aobh6tixIwcPHqRhw4YEBwcTEhKSqrtvxowZLq1GbayIAZcuXQKM1bBtKVGihDW/dOnS6eQvX76coZyFfv360blzZ0qXLs3u3bsZO3Ysv/76Kzt37sRbT77VpGHe188RdWW1dT+60vM0evRFD2qksYdD46aUugIMFZFPMJa66Qe0JyXeoz0EY7HRD4H/U0odyz5VCxYbN27MkvdgZmnWrBnNmjVzWy48PNxq2MBY6BOgRYsW6dJOnjxpbSGlpXPnzuzatcvt49uOj0GK8UubnpGcvfTp06dbvz/88MNUr16dDh06sGzZMh599FG3ddUUXD6aN4ivfbZztkRxBl+6wu4SHWnYb3zGgppcJ0OHEqXUUeAlEXkdY05YM6AOUAYoDlwGzmGEtdoAbFZK3c4phQsKmTUyniJtC6hIkSLp0i1pt27dclhPSEgIxYsXd/m4lpbW5cuXU8k5apml1ddSzoKjlqAt7dq1o1ixYvzyyy/auGmszF77Kd/c/AlE+PKO4tzwKs/Lz32NeLkTxVCTW7jsLWkGGV5tbhpNpnC3W9IyZnbw4EEqVUpZzOHgwYOEhITY7ZIECAwMJCwsLN3YmqMxPFssrbqMWoWawsOOoxcZv/lOat5ZlD8C46kcD0/3XohvEZdcETQeQE8F0OQq7nZL3nXXXVStWpUFCxbQtq0RrjQ5OZkFCxbQvn17p7Lt27fnu+++Y9y4cdaxs3nz5hEWFuZ0ftzq1auJi4sjIiLCZT01BZffT15h4IzdXEvwZ3fMf2lW/nP+0/F9QkvqCCR5GW3cNLlKyZIlKVmypFsyY8aMoU+fPoSHh9O4cWNmzJjBkSNHmD17trXMpk2baNmyJevXr6dpU2P1pOHDhzNr1iyefPJJnnnmGXbt2sWUKVOYNGmStVU2depUdu/eTatWrShVqhS//PIL48aNo379+nTs2DH7TlyTL/nrfBz9pu3k2m3Dg7dEUDCjei+lUslAD2umyQht3DR5nl69ehEXF8d7773H2LFjqVmzJsuXL0/V+lJKkZSUZO3OBMPBZfXq1QwbNoz27dtTtmxZPvzww1TRSapUqcKMGTNYtGgRV69epWzZsvTt25exY8dqT8lCzr4j25n+/VQuXn8MgGA/H74dUF8btnyC2D4MNO4TGRmpdu/e7TD/wIEDDj0HNRpPou9Nxxz+Zy9DfujDGR9ocaYyK68PYubAhkRUKpGleseMGcOYMWOyR8l8joj8rJSKzKn6dctNo9FobLh49gTvr+zDiQCj63pT6N98WlVl2bBpchftw6rRaDQmsedOcm1Ke965cIrw+AS8leKZEo/QqnHW4qNqch/dctNoNBrg0vnTXJncgcrJMQB8cfocy2s9z8BH3/GwZprMoI2bRqMp9Fy5eJbYyR2oknwMgCQlnKg7noGdn/OsYppM43K3pIgMFBH/nFSmoKKddjR5DX1PpnD8zFGGLmiNeBkttmQl/PLAO0Rqw5avcWfMbSpwQkQ+FJF7ckqhgoavry83b970tBoaTSpu3ryJr6+vp9XwOH/G/M6gZV3Y7a94umwoR3x82V1nLPW6DPK0apos4o5xW46xAsBLwAERWS0inUXHKHJKmTJlOHnyJDdu3NBvyxqPo5Tixo0bnDx5kjJlynhaHY9y5sotxixcwjmfZAAu+niz4u4e1O862MOaabIDd2JLPiIiYcALwNNAG6A1cFxEJgNfKaXsL65ViAkODgbg1KlTJCQkeFgbjcboTQgNDbXem4WR47E3eOLLHcTE1qFm3BUuhK2gV2ALhvb41NOqabIJtxxKlFLHgddF5A2gB8Zab42At4E3RGQhxlI30dmuaT4mODi4UD9INJq8xNHzcfT5cgenrhirVxy61Yz3q3WlW6NGHtZMk51kap6bUipBKTVLKdUYqAt8hbE4aW9gi4j8LCJPu7p6t0aj0eQGP/28hP9MmWk1bEW8vZjcJ0IbtgJIlidxK6V+xVix+2uMxUoFw+B9ARwTkQFZPYZGo9Fkle83TGHEr69TpPREwryO4+/rzbT+9WhVI9TTqmlygCwZNxFpJSKLgb+BF4FbwDSgF7ASY0HTqSLy76wqqtFoNJll0+oveOfYJ1zz9mK/nw/3VPicb56OpMk9pTytmiaHcNu4iUhxERkqIgeBNcCjwCngdaCCUmqgUmqeUqoz8CBwHdDGTaPReISdiz6mSfRw/hNrrMpeIimZHnVGUK+yNmwFGZcdSkTkAQwHkijAH6P7cSPwKbBEKZWcVkYptUNEVgCPZYu2Go1G4yIqOZntM16j0T+TQaDntTguSDD1235GvZotPa2eJodxx1vSsq7LDeBL4FOl1O8uyF138zgajUaTJeLjb7NryjM0vrjEmvandxV69l5CqbJhHtRMk1u40y15DBiO0fX4nIuGDeAZQIdC0Gg0ucKlK+d59uuHWCmbsYRN+K1oXcoOWa8NWyHCnRZVFZWJEBumTJK7chqNRuMu/5w9wZClHfjLT/GzXzFKJSXxUGIktf81myJF/TytniYXcafltkZEhmVUSEReEpG1WdBJo9Fo3Gb/qav0/OoAQTeKW9OOB1an7pD52rAVQtwxbq2AWi6UqwHo0VqNRpNr/HjwLI9P3saZq/FsOz2Mmtd9ecK3Ph8NXIO3jx7yL4zkxK9eBEjnOanRaDTZTXJSEjO2/sXYVUdINgdNAooG8EKLVTStpidnF2ay1biZKwREABeys16NRqNJy41b13n1204UvXYbpV4DvKhQwp9p/etRNTTI0+ppPIxT42Zn7KyNk/E0H+AeoBywMBt002g0GrscO3GA4St7cbBoEpSER+Oncix4BF/0jaRUMR3SVpNxy62VzXeFYbjKZSCzD3glK0ppNBqNI47s2UTgkqe4o6wXRjwJuFXiInOeqo9fET3rSGOQkXFrbX4KsBYj3NYHDsrGAyeVUkezSTeNRqNJxa7vP6P2njEUlQQmnPPiiXKhRPrezxsDZ+Pl7e1p9TR5CKfGTSm13vJdRLYCm2zTNBqNJjeIv32LPV/9iwbnFhiv2oBXsj9vVn6DyJa9PKucJk/izkrcD+WkIhqNRmOPo8d/5/WVT/J03Glr2jGvivj0nkPk3a7MTtIURvQEEI1Gk2dZFb2Q9/e/wQU/L0YWKUn4qbPEFWlI1ee+pVhwCU+rp8nDODRuIvK6+XWSUuqSzb5LKKXeyZJmGo2m0JKcrPhi81E+WxtPuXCjH/KWCIvKt+bVJ2fq8TVNhjhruY3D8JBcCFyy2c8IMctp46bRaNzm8o14Xl7wK+sOnAOKk3iyJ6XD5tL/zifo28Gtd2xNIcaZcXsHw0hdSLOv0Wg0OcLOfdG8tPImJy/ftKYVL9OCjx95gUql9eKiGtdxaNyUUiOd7Ws0Gk12kZyUxPg5T7MsYTc1bjflJB0BeOahyrzSrhq+3u6EwdVotEOJRqPxMLHnTvLhgu4sLRYHXl6cKr+BWidrMbh7N9rWLOtp9TT5FP06pNFoPMZvmxaT/H8P8vylQwQmG/HWiyofxvWorQ2bJku4bNxE5AURiReRjk7KdDLLDMwe9TQaTUHk9q0bbJ/0PPdteIpSXCYsMYn/XoilVUJZvum1iTrVm3haRU0+x51uyW5ALLDKSZlVZpnuwJdZ0Euj0RRQtuxdSeza13nkxt/WtIsUp1LkR3zcvLsHNdMUJNwxbtWA35RSDtdqU0olichvGAuWajQajZXkpCQ+mP8C829to1jJZBrf8qJkcjK/+tWjfP+vqV02zNMqagoQ7hi30sAmF8qdA3SoLo1GY+Xc1VuMXLiJw2ort328uO3lzVslS9Ir5Gka9ByBeOnhf0324o5xuwK4aPzFDAAAIABJREFU8mpVHojLnDoajaYgoZRi6a+nGL3kD67cTKROseZcDNtExXh4tP57NKz3mKdV1BRQ3DFue4DmIlJFKfWXvQIiUgV4EPgpO5TTaDT5l7MXzjFm9XFW/X7GmrY3rj39lRdDo0ZzR5CelK3JOdwxbtOBNsD3ItJNKXXENlNE7ga+A7zNshqNppDyxZKRzD//HSHHHwEMz8fyd/gz4fHaPFjFocO1RpNtuGPc5gF9gA7AHyKyBTho5t2LMc7mA6xWSs3MVi01Gk2+4MrFs7w/vwdLA2LB14uQcksI+rs2nerX5L8da1CsqI4bockd3FnPTYlIN+Bj4BmgmblZSAQmAcOyUT+NRpMPUMnJ/LJ6BpV2juGJItdZ6V+WRBHO+gjvtfemQ9PanlZRU8hw6zVKKRUPvCgiY4GWQCUz6x9gvVLqjENhjUZTIDl78iinZr1IxI1tAJSKh2cvX2FPQCVGPfItYWUre1hDTWEkU30EphGblV1KiIgv8DBGl2djDKNZEjgPRAOfKaU2OpHvDbwA1MYY8zsIfI2xFp3DeXmZldNoNJCYmMAH81/gwX/W8PDtq9b0c4RQv/pbvND2SQ9qpyns5JUO8KbAD+b3M8DPwHWMyeCPAY+JyFil1Oi0giLyOTAIuAWsBxIwWpWfAS1F5HGlVFJ2yWk0Gtj++0Y+3DaUg0WT2FI6kPonr+GnFDtKPkq1Jz8i4o6SnlZRU8hx27iJyL3AvzHG28qbySeBDRgtrIMORJ2RDCwCJiqlNqc5Xk+MVuIoEdmglNpgk/cYhoE6Azxs8eAUkVBTn67Av4CJaerMlJxGU9iJT0xmyqa/mLPpIMmVEwAv/vH15f+Cy9Gu4Yc0aNTe0ypqNICbqwKISH9gL/A8UB0INrfqGMZir4j0c1cJpdSPSqnuaQ2bmTePlKkFfdJkv2Z+vmo7NUEpdRajuxFghIikPc/Mymk0hZZtf12g/cSf+PCHw5yKD+Ou8zXxUYr2SeEMeGYLNbRh0+QhXG65iUg94AsMg/gdMA34CxCgMvA0RnDlL0Rkv1JqVzbqucf8rGCjTwUgAogHFqQVUEptEpGTGK3LhsC2rMhpNIWVw8f2Mv+H6Xz5V+tU6XH+g/jffTdoGtHBQ5ppNI5xp1tyOIZh66OUmpMm7yCwSkR6YXQhvgz0zB4VAbjH/Dxtk1bX/PxDKXUT++zCMFJ1STFSmZXTaAoVt27f4MMFL7AsYTfe3orqPuU4kFiTwCLeDGtzL/0fDMfbSzytpkZjF3eMWxPgZzuGzYpSao6IDMXwfMwWRKQs0N/cXWSTZfEv/seJeEyaslmR02gKDQd3ryd59X/YWvYW1319Abg/dCZ3l5zOyE41CA3287CGGo1z3DFuJYEfXSh3BKiTOXVSIyI+wEygOMY8umU22cXMz+tOqrAEcA7KBjlbvZ4FngWoWLGik2o0mvxF7LmT/Dn3VerHGn+11y76MahsGconKJpU7U/vtg94WEONxjXcMW6XgCoulLvLLJsdTMZwzz9OemcSS3+IcrPOzMpZUUpNBaYCREZGZroejSavcP3GNeYveo3H/lpEfW5Y0yNvJPOciqBfr/8RFHiHBzXUaNzDHeO2DegiIl2UUkvsFRCRzhhOGN9lVTERmQgMwHDXb2kn+sk187MYjrHkXbNJy6ycRlMgmbHibeacnsMZH3jIN4HgBCN9T8CDhPb4H/8Kv9ezCmo0mcAd4/YR0AVYICIzgRnA3xgtoLuAvvx/e3ceHkWVLn78+2ZPyA5hFwYMsiurKAKigIriKCCjMy6go86o12Vm3OZeZ+aH4zbe69xxX1HcdZRFxBEZRQQFhCD7GtaEQNhMIGRP9/n9URXS5CZtluquTvJ+nqeeSupUn37rpLvfVHWdc6yzK6+9b4OJyFNYfekOYyW2zBp222Ovu9ZQVqly/rk9Ptsa+jilmpUdhwp4+NNNFJV9RE6cdUHjydbJPHQwkh9H/IWBY65xOUKlGq4+Ayd/a98s8ndgqr34EsAD3GOM+a6hAYnIk1iDLx8FxhljNteya2X3gL4iElvLnY9Dq+3bmMcp1SwcKyrnH19t5+3le6nwGrrHXEXYzz4g1hg6txpG2wdepUtMnNthKtUo9R04+Vl7qpvKOyI7YiW1HOAbrBFGGpwQROQJrC4HeViJbZ2fWLJF5AdgEDAFeKtaXedj9YvLxRqfslGPU6qpKyop5PU5/8mb28dxuLjqFv7dpQO5zpvL1LG/Jb1LPxcjVMo59R5+y05e9R6F5KfYMw08AORjJba6JMnHsTpi/01Elhljdth1tQVesPd5ooZBkBv6OKWaHOP18tq8P/HxkU/YHymMiszls+IbABjWLZU/X96Hvh11AlHVvITEwMki8nPgIfvXHcCdIjV2Dt1qjHmi8hdjzMci8iLWkFkbRORLqgZATgTmYg2EfIqGPk6ppmbzigWEf/kXDiYdZH+i1bMlO209/cKKueOy87ikX3tqea8p1aSFRHIDUn1+HmIvNfkGeMJ3gzHmdvtS6R1YswtUTl3zOn6mrmno45RqCvZu/YEf5/0XA+051n6bF8an8a0IMzA4+ix+/x9jSYjXW/tV81VrchORVxpRrzHG/KYeO8+kanDkhjzZe8B7wXqcUqFq2+41vL7wPqbvz6CrVHXBTPSE8Zuy/oz7+V85rUMPPzUo1Tz4O3O7uRH1GqDOyU0p1TjHist56sMb+NxsoCRGGJjYimsKrIF2MhLH0nHSY9yk/dVUC+Ivud0StCiUUg1SWFrBzGV7ePmbnZzZqoiSttb3Zy+mJJFefjrJ4x9jyFkjXI5SqeCrNbkZY2YEMxClVN0VFxfzbsYBXly8k6OFZQCsKr2WDqnTSfCGcXXnqQyaei9h4eEuR6qUO0LlhhKlVB0UFhXw3NzfsaRwOZ49N3LU0/tkWYfUNH7T/VEmnzeeqKhoF6NUyn0NSm4iEo91R2MakGWM+d7RqJRSp/BUVPDD/Jd46cALrIwLh6gwxqR9yNbc/0en5FjuGpPO5EGdiQjXyeOVgnomNxFJAJ7CGkcy0t78JvC9XX4b8EfgKmPMSgfjVKpFqigvY81nr9J+3XMMNfvJi4tlZVwaAFvii3j04rZcNXIQ0RF6+VEpX3VObiISByzGmp36CPADcFG13RYCzwMTAU1uSjVQUUkhH897mDFb5jDUHDy5fWxRMUOKy+kaN5C7pzxDSlKai1EqFbrqc+b2B6zE9j5wqzGmUERO6ehsjNkpIpnAhQ7GqFSLUVJWzjOz7mBh4XccihBGRBy1xs4BjhPHpi7X88zkB0lISvVfkVItXH2S2y+AA8CvjTElfvbbC/RpVFRKtTAl5R7+mZHNi4t30iUpg4OtrO/OXklO4o+HS9jc9Xr6Tryfc5NbuxypUk1DfZLb6cAXP5HYwLpk2abhISnVchSVlvP+qn28/M1ODhWUApBUNh5azSPJ4yU6/mzCr3+RcxNTXI5UqaalPsmtHKjL/cWdgRMNC0eplmHv/u28/MUD7C/JZMXe6ZQQc7LsSPhopkaVcOP4+2id3N7FKJVquuqT3LYDA0Uk2hhTWtMOIpIMnIVO8qlUjXKzMtk+/wn+M2YFx8LDIE4YHf8JC05cTduEaH57/un88uwuxEaNcztUpZq0+iS3WcBj9vKHWvZ5BIjHmitNKWXbsyWDwwueZED+l4wSDxe1TuEjewoaUjby8JiH+cWQ04iJ1Fv6lXJCfZLbs1iTlN4jIkOwkh1AVxG5BWtW6zHAJuA1R6NUqoma+/XLHFv3PlPz1/AzsOatB6YdK2B1dAJjU8dxyzWPEBMd52KUSjU/dU5u9q3/F2EltZFA5Wiso+1FgLXAFbVdtlSqJfB4DfOWL+TdTfezLdpLp/gKrs2verNtiupPxbl3M/f8yUiYjiiiVCDUa4QSY0w2cLaITAAuBbpjTfKZDXwOzNJJPlVLVVhawcer9/H6d7vJ/fEYaT0qgDByIiNY2CqODmYgsRf+gb5DxrgdqlLNXoPGljTGzAfmOxyLUk3S+swVfLVuP29sbMXxkgp7ayu653dkQ+oBhpWlkDzuSQYOGO9qnEq1JP5m4v4YmAEsMMaY2vZTqqX697L3+GDjs6yOKmB0XhzHS/5ysiwxJoI+3R/g/n4J9Esf5mKUSrVM/s7cJmGNEXlARN4EZhpjMoMTllKhqaK8jPVfvkvcD68gkbtZ2S4NEDKSTtDlcDZhKT25aUQ3Jg/qTKtonVFKKbf4e/e9CFwNdAQeBB4UkW+B14GPjDFFQYhPqZCwP3cXWV++SZed7zHIHAKgRzl0Kq8gJzKCThXRXHdxGpeNHE1YmLgcrVLK30zcd4jI74ArgJuAsVTdJfmsiHwIvGGMWRaUSJVywcLl7zNnwwusiszjvf25dDTlJ8s8Jpyri7vRud9Uxp1ztYtRKqWq83vdxBhThtUh+yMR6YjVz+0GoCfwa+AmEdmOdTb3tjEmN8DxKhVwpRUe/rXhAG8t30us9+9sSCgBhA8T4/nT0TzySGBr5yn0uPQebuzY1e1wlVI1qE8/t/3A48DjInIu1tncFKxE9wTwqIh8jpXo5htjPAGIV6mAyc7N5YN1+XywMpujhWUA9Is7HxK+AGBjVCIr+9/DmeNv5ty4eDdDVUr9hIZ2BVgOLBeRu4CrgGlYHbkn2MthQEd8VSHP6/Hwz6+eZsHuDzhBARt3P3LKIMbbSy9gfOlGxvSawoTzphEWrsNjKdUUNOp2LmNMMfA28LaIjAPeAdLsRamQlXf4ANsWvkLirg95qrNQEhMGRDCq1b9YWDiJDkkxXHdOV64eehpt4i9zO1ylVD01KrmJSDzWHZXTgOGcHDmP7MaFpZTzvB4P65bPw/P925x5/BvOEavD9WWFqcxKsC4zJqdm8dLEwYzt3ZaIcB0aS6mmqkHJTUQuAG7E6gsXi5XUSoF5WN+5LXQqQKUaa2f2Jt75+q8sK9vI1QX53FRQUPVvGHDFsTJ+jG7LlCG/Y+TACe4FqpRyTJ2Tm4h0w7pbcirQhaqPh7XAG8A7xpg8xyNUqgG8XsOynUd5f2UWefv+zob2myBSmJUQz7RjBYQB2yPOIL/3tfS9aBrPJCS7HbJSykF+k5uIxGHdETkNq4+b2Ese8B4wwxizNsAxKlVnu7O38PmOcD5YlUPWj9Y4AzEygZS0NZwID+NoeDjz0y6l76i7OaP/OS5Hq5QKFH9jS87ASmytsBKaF/gS67LjHLsPnFKuKyst4a0Fj7Lk4AI2RhXTfdckssqqxnMsMQkMLu5L9zZtmXrxQ7RO1ht5lWru/J253WivdwMzsUYj2RfwiJSqox0blnNk6ev0PLSAjHbRrImzvv7tnPJvMg4OIzEmgkmDOvPLs7vQs73e8ahUS+Ivub0LvG6M+TpYwSj1U3bnbGXNN69z1s5/k+7ZRbq9/coTcXwXFwtAUUwZf5/Sn0vP7ERMpPZLU6ol8je25PXBDESp2lR4vMxf/hmfbHqEtdGFDCgpZZLn0Cn79CmKY4KnGxOH3sXZ/ce6FKlSKlTonBwqZG05cJw5a3KY/UMOlORQkV6IR4SM2BiyI8JJKw9jU9IoooZcT5/hl/N4hL6clVIW/TRQIWV95gpmLfsH2UfGsOhwR5+S0xhWFMXmVuX0KBW+S7+V8Rffy+CUNq7FqpQKXZrclOuO5R1h26J3+DD3bf7dqhQjwoWefODOk/ukJUQzoM1vubNHR0ZoR2ul1E/Q5KZcUVpSxOYlszHrPqTvieWcLeVsTUxgYXwKAEcT9xJXAGP7dGTiwE6M7NGGiHD9Lk0pVTea3FTQVFSU8+nSGXy14wOSS3J45Kh9U4g91s34wkL+NzWZ3mVRnJ06hpduuJB4+w5IpZSqD01uKqCMMWzaf5z56w+QsfEztrWdCVEQHRHFAz8KCcYAsCP8dI50v5LZIybTrWtfd4NWSjV5mtxUQCzOmMuGA6nM3uZl95FCe2sv+ifBnmgoDQvjk7h2nJ48no4jbyC916CTfdaUUqqxNLkpx2RlruO9b59gacVGsqJgVE4vdh+fdso+rQv60lVyGZN+DZdf+2siIiLdCVYp1axpclONkrNrK1lL36Zt1r843bOLmJQkspKTAKhI3ArHoVVUOOP6tGPCmR0ZecYlREfoqCFKqcDS5KbqbX3mCuYuf4a4vC3cm7eLTj5llxQWMSM5iWivISIikZd/2Y/z+3TWYbCUUkGlyU3Vyd6jhXy+MZdV6+eyMnEmAKkJHu7Jq3oRlZpIiiKGcmfCACZdeBdtkju4Fq9SqmXT5KZq5PV6+XbtfDYd6sS8bUVsOXAcACGdrnFejkaE8WN4OMtj4kiU/pT1upJeo69hUFIqg1yOXSmlNLmpk4zXy471y3hv9VOsMJnsixSGZQ1jS+HEqn2IoMuJNnSMKWBI8jl0u/E+Orfr5mLUSin1f2lya+G8Hg/bVy8if/Usuhz6ih7mEN42qexLiAcgNmEdFE4kKjyMkT3acEm/9ozt9QUp8TEuR66UUrXT5NYCFZUUMvebl1iR9Rm9CrK5/fip08eMKSxidkI8sV4vETGxPHPNAC7o1ZaEGL1tXynVNGhyayGOl5SzeNthvtx8kENZr7Kp/VKIgr0J4dx+3Gc/4oiIGs7dKX2ZPPouUhJ11H2lVNOjya0ZW799GYtWzWR98fUs3V1Aucca6ipGRhHfdgmlYcKuqEjWRCZTHj+CmDMn0mv4BIZHxzDc5diVUqoxNLk1I8brZeeGZRzMmMP/ehayLdra3i87lXLPyJP7lZgEhhemkhAdyfCfXU7va24lJjrOpaiVUsp5mtyauMKiArZ8/xmeTV/Q9cgS0jlCOvB+uzS2YY2o3y7heygdSf9OSYzt3Y5xfdrRu8OliIi7wSulVIBocmuCjpwoZc6iZ1lx4GM2RRVwZ14+1x4/cco+FxQW8V1sDH1Ko+nRZQjLb7mQDkk6fYxSqmXQ5NYEeL2GjfuPsWjrIb7edpj1+/I5N2UVG9oVAmF8HRd3MrkdoxWZicPp0PMiPh82nvZtTnM3eKWUcoEmtxB16EgOs5c+x9rD3+ItP8bC7Mc4OasnkHn8fGi3GYCciCi+bTuFlLMmcsbQsQyJinYpaqWUCg2a3EKE8XrJ2r6WAxmfkJC1iNSKLTzftQNEQ0SUoXf4NrZ4egEQHiZ0O20gA8LPZXiv8Ywa+HPCwnVgYqWUqtTik5uI/Aq4DTgTCAe2Am8ALxpjvIF87vyCI8xb+goZOV9yR242PT2H6OpT3qe0lM3R0VSI0Cv5e3p3HcOFvdoyMj2NpLhI4NxAhqeUUk1Wi05uIvI8cDtQAnwFlANjgOeAMSIyxRjjcer5jDFsP3iCJdsPsyTzMEXlt7I91gtRcGnMCXoWnrr/+Sfi6B7RhfN6TOTiX/6KSL3cqJRSddJik5uITMZKbLnAKGNMpr29HfA1MBH4D+DpxjzP3v3b+WzFa+TnhTPv4Fhyj5ecLBuV1h5i9wOwNDaGEScM2+OHUpE+ju7DruD2jl1rq1YppZQfLTa5AX+01w9UJjYAY8xBEbkNWAw8KCLP1ufyZEV5GTt+WEz+hgWsPf41L6SWY0QYXgy5x0ecsu/hE4M5LTGH3tKJQb0uJ+qWmxkUrQMSK6VUY7XI5CYinYHBQBnwUfVyY8w3IpIDdALOAZb5q2/ttm/5Yc1sBmRvJ71wNb0oAiApKpLnxZqwc0uch9bkUx7ThpE90hh1RhtG9riQjsmPOXtwSimlWmZyAwba603GmOJa9lmFldwG4ie5ZR7ZxPUrbiPB42VJ4b5TGrRnWTlpFRWkeKLoGXE6V15/FoN69SEiPMyZo1BKKVWjlprcKmfX3Otnn6xq+9ao8nplQXgYG6KjGFhaxkFaszflHMJ7jOODoRfRNq1TY+NVSilVDy01ucXb60I/+1SOZ5VQvUBEbgVutX8t3Tht40aAQSf3OA7sBt5vbJxNTRvgiNtBhAhtiyraFlXaTJ8+XdvC0jOQlbfU5FY51IdpyIONMa8ArwCISIYxZohTgTVl2hZVtC2qaFtU0baoIiIZgay/pX75U2Cv4/3sU1lW4GcfpZRSIailJrc99tpfR7LKEYf3+NlHKaVUCGqpyW2Nve4rIrXNAzO02r61ecWZkJoFbYsq2hZVtC2qaFtUCWhbiDEN+tqpyROR1Vj3gEw1xrxVrex8rE7cuUCnQI8xqZRSylkt9cwN4HF7/TcRSa/cKCJtgRfsX5/QxKaUUk1Piz1zAxCRF7BmBCgBvqRq4OREYC5wlZMDJyullAqOlnzmhjHmduBaYDMwDvg5EI3VReBKrMGTG0xEfiUiS0XkmIicEJEMEblDREK+3Z2MXUQ6i8izIrJNRIpFpEREMkXkJRHpHoj4neT031FEYkXkfhFZJSL5IlIkIrtF5CMROc/p+J0UyNe0iDwmIsZe7nUi3kByoi1EJFJExojIUyKyQkQOiEiZiOSIyMciMjqAh+CYALxHGl+fMabFL8A/sBJa9eWqRtT5vF1HMTAfmIPVu9sAs4Fwt487GLFjDV+WZz82G+uMeC6wz95WAAx3+5iD9XfEGvEm0378QeAT4J/ASqyxTh9y+5iD1RbV6h4KVGAN+mOAe90+3mC0BTDW5/PmgF3Xh8AGn+0Pu328wXxdONa2bjdMKCzAzcCTwC+A07FuJmlwcgMm+7xYe/hsb4d1lmiAu90+7mDEjjUup8G6MyrSZ3skMMMuW+f2cQepLVoBOyo/sHzbwy5vDZzh9nEHoy2q1R0NbAJy7A+ykE5uTrYFcCHwMTCyhrKrsRK+AS5w+7iD8bpwtG3dbpxQXBxIbhn242+ooex8nz9emNvHGsjYgRiq/vtsX0N5R5/yOLePPdB/R6ybmAzwptvH5nZbVHv83+zHXw7MbALJLWjvb+A1u74Zbh93MNrC0c8ftxsnFJfGJDegs/3YUiC2ln0qL8mF1OU4p2PHOjsrt/fvUEN5B7vsBPbNTaGyBKAtorDGVzRAb7ePz822qPa4YVhnJ+/av4d0cgv2+xu4w67rC7ePPdBt4XR9IX9jQxNU1+l0fPcNFY7GbowpB76yf50uIpGVZfbPj9i/zjD2KzeEOP13HIx12THbGLNFRIbbN1C8LCLTReTcxgYcQAF5TYtIDPAm8CNwd8PDC6pgv7972OsDDtTlNKfbwtH6WurAyYHk2HQ6LghE7LcDC4BbgPE+g6UOBVKAp4H76hlnMDjdFv3tdaaIzASmViv/s4jMAq7388Z2S6Be049ijQx/jTGmqYyUH7T3t4i0B6bZv85qTF0B4nRbOFqfJjfnNWo6HZc5HrsxZpeIDAfeAsZjXXqolAEssc/wQo3TbZFqr0cB4cD/AC8BR+1tL2B9mX4cuKm+wQaY468L+zVxDzDXGPNhI2ILtqC8v0UkAngHSAK+MsZ82tC6AsjptnC0Pr0s6bxGTafjMsdjtz/ENgLpwBVYc3ulYfUjTAFmicifnXo+BzndFpXvtQisy7D3GWN2GmPyjTHzsNrDAFNDsO+fo21hj+f6BlYiv92JOoMoWO/vl7AGlMgGrgvwczWU023haH2a3JzXlKfTcTR2EUnG6tOWAFxijJlnjDlqjDlijPkEuASrL8ufRKSHv7pc4PTf0XefV6sXGmMygNVY78nRdagvmJxui8eAM4DfG2NC8bskfwL+/haRp4FfY41tO8YYk9uQeoIgUO8RR+rT5Oa8Pfa6q599QnU6nT322qnYL8M6S1thjNlVvdAYswP4HutsZnRdgwySPfbaqbbw3Wd3LftUbm9fh/qCaY+9dqotJmJ11p4qIot9F6x/eABus7e91oB4A2mPvQ7I+1tEngLuAg5jJbbM+tYRRHvstdPvEUfq0+/cnHfKdDq13BxQ1+l0gs3p2LvY62N+9sm316l+9nGD023xg8/PrbE+vKprY69P1FDmpkC8psOw+i3Vpru9JNexvmAJ2PtbRJ4Efo/1Pew4Y8zmhocZFE63haP16Zmbw4wx2VgfZFHAlOrl9nQ6nbEuOSwPbnT+BSD2/fZ6sG83AJ/6IrFukYfaz2Zc4XRbGGNysM5SwfoupXp9KVhTMIF1o03ICEBb/MwYIzUtWF0DAO6ztw1w7kgaL1DvbxF5Auuu4TysxLbOkYADKACvC2fb1u2OgKG4UIdO3FijTWwFHq+h7CqqetKn+2xvizXMkCF0h9+qd+y1tYX9mEL7Mc8B0T5l0cCLdtmPQJLbxx7ItrDLLqdqTMkBPttjgA/ssgxCrEN7INrCz/PMJIQ7cQfodfFX+zF5wGC3j8/ltnDss9P1xgmFBes/5hU+S+Ugndt9t1d7TOWbcGYtdb5A1eCfn2IN+HnM3jaH0B44uV6x+2sLrP5clePj5QDz7Dr329tKgCvdPuZgtIVd/t9UjcKwxK4jx962D5/x9EJtcbotanmOyseEbHJzsi2wZiIx9rLK3q+m5UG3jzlYr4v61ldrXG43TCgsWDczmJ9a6vMHsvf5FfAdVrIsxLob7g5CcEzJxsRehxfrIKx+bruxklkJsBNr3Lw+bh9rMNvC3mcisAjrP/VSrFkCngLS3D7WYLeFn8eEdHJzqi2wOmn/5GcPsNjt4w3m68KJz84WPVmpUkqp5klvKFFKKdXsaHJTSinV7GhyU0op1exoclNKKdXsaHJTSinV7GhyU0op1exoclNKKdXsaHJTqhkRkYEiYkTkdbdjUcpNmtyUal4m2es5rkahlMt0hBKlmhER2YQ11VCaMabE7XiUcoueuSnVTIjIGUAf4HNNbKql0+SmVBDZ34cZ++dpIpIhIoUikisiM0QkzS6LEZHpIrJdREpEJEtEHq1pXjwfk+31bJ/nG20/52K7zr+KyA4RKRaRXSLykIiE2/ueZseQYz/nBhG5LlBtoVQg6WVJpYKoMrEBTwL3AN8ABcBwoD2wHjgP+ALobZdHY81aHQcnuY2BAAACjElEQVS8aoy5tZa6VwJnAW2MMQX2ttHA11iTO3qAvljzFcYCo+w6XwL+B2sU9iJgJdAJGGFXfZ0x5l0HDl+poNHkplQQ+SS3g8AFxpgt9vYUrATUE9gI5AMTjDHH7PIBWPN9hQPdjDF7q9XbGcjCuiR5mc/20VjJDeDbanWe5VPnNmAh8AdjjMcuvwNrktmdxph0B5tBqYDTy5JKuePPlYkNwBiTh3UGBdb3ZrdWJiG7fC3wL0CwzuKqm2SXza6hDMBbQ53r7DrDsM7g7q9MbLaXsWZJP11EuvhWJiI3iMg2+/LmThG5vS4HrVSwaHJTyh0Lati2w17v9U18PjLtdccayiZhXXacV8vz1VZn5XMuMsaU+RYYYyqwJpg95TlFZALWRLN/AhKBW4DHReTaWp5bqaDT5KaUO/bVsO2EnzLf8hjfjSLSBuv7sW+NMYfr8XwNfc67gU+NMf80xpQbYxZhJbu7a6lDqaDT5KaUC4wxXj/F/spqciXW92a1XZKsS531ec7BwPfVtq0EBoqIfqaokKAvRKWavon2em6Qni8JyKu2LQ+IwPruTinXaXJTqgkTkQRgDJBhjMkK0tMeA1KqbUsBKrC6EijlOk1uSjVtE7D6wQVzLMnVwLBq284G1vzE5ValgkaTm1JNW+VAyf6+b3Pa08DlIjJFRCJF5ALgZnu7UiFBk5tSTZSIxACXAFuMMVuD9bzGmPlYt/8/ijW6ygzgjzqKiQolOkKJUk2UiFyBdRPJY8aY/3I7HqVCSYTbASilGqwYmA6843YgSoUaPXNTSinV7Oh3bkoppZodTW5KKaWaHU1uSimlmh1NbkoppZodTW5KKaWaHU1uSimlmh1NbkoppZqd/w+NPh5K2HI9twAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"v = lambda mf: u*np.log(m0/mf)\n",
"\n",
"plt.plot(rk2[:,2]/m0,rk2[:,1],label='rk2',ls='-')\n",
"plt.plot(heun[:,2]/m0,heun[:,1],label='Heun\\'s Method',ls='--')\n",
"plt.plot(heun[:,2]/m0,v(heun[:,2]),label='Tsiolkovsky', 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",
"Use [two integration methods](../notebooks/03_Get_Oscillations.ipynb) to integrate the `rocket` function, one explicit method and one implicit method. Demonstrate that the solutions converge to equation (2.b) the Tsiolkovsky equation. Use an initial state of y=0 m, v=0 m/s, and m=0.25 kg. \n",
"\n",
"Integrate the function until mass, $m_{f}=0.05~kg$, using a mass rate change of $\\frac{dm}{dt}=0.05$ kg/s, . \n",
"\n",
"Compare solutions between the `simplerocket` and `rocket` integration, what is the height reached when the mass reaches $m_{f} = 0.05~kg?$\n"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [],
"source": [
"def rocket(state,dmdt=0.05, u=250,c=0.18e-3):\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": 31,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The height of detonation: 426.36 m.\n"
]
}
],
"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(\"The height of detonation:\", height2, 'm.')\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAbcAAAE0CAYAAACvs32dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOydd3wU1fbAvycFEgIBEpqGSABBmigkNEHpRQERRIoPAbvi44m8h+2hovAsj4f+lOcDsWEB6SiIFEFAwChdlI4Q6TW00EKS8/tjZpNNspvsprAp9/v5zGd35p4zc2Z2Zs/ce889V1QVg8FgMBiKEn6+NsBgMBgMhrzGODeDwWAwFDmMczMYDAZDkcM4N4PBYDAUOYxzMxgMBkORwzg3g8FgMBQ5jHMrYIjIJyKiIpIsIlW90AsXkURb98s8tKeLvU8VkSp5td/cIiJBTnb187U9hQ0RecLp+qmIfO6h3qYMes3z21Zf4+JaqdMzelpENorI2yJSy9e2OhCRN20bd/jaFl9hnFvBw/En4wfc74VePyDQ/v5FnlpUCBGRafbDvcjXthQSeolISFYCItIAuPUa2VMY8APKAY2AZ4DfRGRgfh5QRI7a9/Xz+XmcooBxbgWPlcCf9vcHvNBzyB4BluapRYaizjkgBOiZjZzjj/tc/ppToGkHlLGXcKAp8H9AClAS+FhEbvGdeQYHxrkVMNRKGeNoVmzgyYNiN4c0s1enqGpyftlXUFDVy6oq9jLN1/YUcmbZn25fpkTEuSVhZr5bVHC5pKoJ9hKvqutU9RlgpF0eAAzzoX0AqOrz9rNRx9e2+Arj3Aomzv0fntTenGU86jsxGJxwNGO3F5Hr3Mi0ByKAM8D8a2JV4eJt4Kr9vZUvDTFYGOdWAFHVXcAv9ur99ltzVgywPzer6m+uBEQkUEQeFZElInLMDj45LiILRaSviEhubBaRtnY/1wERuSIi8SISKyIjRKSUB/rBIjJURJba/QpXROSIiPwiIq+JSN0M8i4DShyd/0Bfe1NnF4EAE23ZVfb6cg/s+9iWPejB7+HQ2WLrfOOB7FRbdo+Lsu4i8rV97Csicl5E9orIChF5SURqe2JPFuwA1gH+wF/cyDheoGYAV7LamX2vtROR/7ODLc6IyFUROSkiP4rIsOzuCRG5xb7mu0Tkoohctu+tdSLyjoi0dqOX39fKJap6Bdhnr1bK4rwiReRJEZkvIvttGy+KyB4RmSwi0W70ptn3dWV70xsu7mvn5yDbgBIR8ReRB0Xke/u/INH+b1iQF/8JPkdVzVIAF2AIoPbSMQu5Vk5yw93IVAe2Osm5Wr4Ggl3odnGSqeKi3A/4Xzb73gvUzuIcYoD92ezj5ww6QU5l/Zy2P5HNfhSYaMsOstdTgOpZ2FcKq59JgX958RuOsHUSgfAs5EoDF2zZVzOUfejB+fwnB/eX83WqAgy1v//qQjYESLDLW2W4J5q7kH/OA5t/ByLc2DYQSMpGf70LvWtxrTKdr5PcblsmLguZS9nYl4yL5xiY5sG5OT8Hb9rbdrixIxz4KZv9LQJKe3u9Cspiam4Fl2lYf4qQddOko9aWDEzNWCgiYcAKoB5wAqs/oA4QZn++hPUm3gMYnwM7RwFP2t+XA22BikBt4GV739WBxSIS6sK+2sAyIBLrD340VkReOFYzWBfbrjMe2vMRVmf/bHt9KWkBAI7labtsJpbTEixH547eth7Apx7aAdbvkYIVxdonC7meWA4U0vpbEZGuwCNO+7odqIr19h4N9AfmAJe9sMkd07AcSkMRaZihrBeWg9sHrPFgX5eBecCDQHOgGlZt5lbgeeAoUB+nc3UgIhWBiVi1yLVY92V1oALQAOgKTABOZtC7ltcqEyIShHWeANuyEN0FvAV0wroGFYAawJ3AN1gvi/8RkbYZ9AZh3YPH7fVXyHxfz8ID7JaH2UALe9NErIjPcKwXTUczdWfgM0/2WSDxtXc1i/sF62FU4DxQykV5CSDelvnOzT4m2eWngRpuZLqS9rZ2c4YytzU3LId01S5bAgS42HcvJ/3XXZSvsMsuANFZXIuADOsua25O5Y433UXZXOMPbLl9gLiRWW7LrMzBb7jU1l2ThcwiXNdO37e3x+bDvZWu5mZvm2+vj80g+729/TUX94TbmkwWx77BvqcVaJGh7D57+xUg1It9Xqtr5fJ8gRedZLrm4ljv2vtY7Kb8qF3+fDb7cVtzwxo25LD1ZTf67znJdMrra3otFlNzK9g4gkNK4zpMuxtQPoNsKiJSjrTw7ZdUda+rg6jqAqwmCvBubN0grOgwgKGqmuRi33OwHB/Aw87t+HYNobW9OkZVN7g7kKt95xEf2Z9RQJuMhSISRZqN3tTaHDhqJ7eJSHUX+68MdMgg68BxbQ/l4Lg5wXEPpfbzisj1WOHvkEfjJ1V1P9YLA0DHDMWOcz6H5QA95Vpdq2ARKW0v5UUkRkTeAV6zy9+0n6ec4qgptRaRErkz1S2OGu5+4HU3Ms9jvRA7yxcqjHMr2CwgrfnFVdOkY9s5rCaNjNyONfYG4EenhzLTAvxqy8V4YZ8jKmyLqu7MQs4ROl4Jq7nSQQen75O9OG6eoarrAEcQzmAXIoOxmi3Pk7MQ+NlY/SzgOlijH1YTXBIwPUPZJvvzHhH5q2QzyDoPmA+cBa7Hio4Eq9nbD6tWudvTHdn31d/sYIUjdkBIavAD0N0WvSmD6mb7swLwgbiP3szItbpWP2DdC+exWk3WYTX1XwDaq+oL2e1ARJqLyCQR+V1EzolIitN1cbzglSStmTPPEBF/rKZigK/dvTSq6kXgO3v19ry241pgnFsBRlWvkvaH10Gc0l/ZfWl32auzVPVSRn3S/3H8StpD6Wpx9JtV9MJET/oYwApmyagDUNP+PK6qR7w4bl7zsf3ZW0QcfWvYtUxHX9wMVb3g7Y5V9TxpLx4DXIg4ti1W1RMZyj4BtmA5v/HASRH5Qazo0Q4iEkgeoqqXSXPgjhq/4wXK41qbWJGtW7Ga2DpgBayUdCNeNoMN27H6gAAeBQ7ZEZfviUhvuzXCFdf0WrkgFPg/EamQlZCIjAVisc6tPlZfmbuoxLJutueGClj9p+D5c1slH2uR+YZxbgUfR1ORP1anuIO+WH1uzjIZycnDEeSFrMMRJGQj59y8VMbpe6iLcl/wJVYfTymsPh8HbbCaKyFnTZIOHI7hJhFJrRnbwTQxGWRSUdVE4A6s/pNjWL9NW6wgoO+Bo3Z4e17+cTvupZ4i0goriCORzLVKl9i2fI3Vr3YOq7nudqz+2XKkBT/MsVUCXOzmKax+rm1Yf/yNsKI5ZwLH7JD5dOH21/BatVA7eQDW89USq4UF4GayeAkQkUHAP+zV5VjPcz2sF0rHdXFuOXF1bXKL8/OX0+e2UGCcWwFHVddijUOC9E2Tju9/Aj+6UXfcvClASU3L6JHV4k1GA8fNXzobOefy8y6++/TBUdVTWH/IYEX4keH7TlX1JErQHUtIi3Jzrr05vp/Hii50ZdtZu6nrOuAW4HHgK6zmwzAs55GXA/dXA3FYb/eO/p/v7GvkCR1Ja3ruoaqvqOpqVT1on0uCqiaQxT2jqimq+oGq1seq6ffHGm6yH+uFbhCwxm5Od9a7ptdKVc+p6k/A3aQ5uC4i4i4y9in7czlWE+Y0Vd2uqiedrkt+15Ccn7+cPreFAuPcCgeOt8FGIlJPRGqSFsb7hdrhTS5wBJD4ARnDu/OCOPuzXjZy9V3oADgGLFfyom8lv3A0TbYSkZp28+S99rbc1NocwTCOFGH97H4PSAveme2mWdl5H6qqW1R1kqrejxXm7nCI/UQkT9Is2feSI7Clhv3pTSCJI7HyUVVdkYVcAw/t2W87gaewhgQ40lzdiOtm3mt2rZyOlwI8htXvBjDa6Td2xpFKb3oWz+zNeWmbC06SZqenz+0Ru2ZcqDDOrXDwBVbtC6wa24AMZe5YjhWoAOlrJHnFavuzoWQ93Udv+/O4WtlXHDgneM5qnFlOcKRCcvUn44qlpCWsHozV7FsKa/xgXrztOxxGZaz+0xak9Tl6PUWR/Zb/ptOmuu5kc4Dz+Z4GvvVC19G35va6i0h7rKAVr7CdyBukBeh4dM75fK0cxziMFT4PVs21r3O5HX3qaBLN6p506bCd8Pa+TodaeWdj7dUebpwwIhKMNUQI0p7zQoVxboUAVT2ANVsAWG/7jgfglwzOIqPeSdL+qB4Xke7uZMEaOmCHpnvK56Q5z/dcPSgi0gNrXBSkhd077NuCNc4NYKSINMrCNm/7HxzNaB79idpv0o4a2iDgIfv7orwIdrGjMh3NywNIi5w8RFpYfDo8qGHUdPruabNhtthRkbWxnMAtXr61O1oLKtoOPB0iUp4skgWISI1s+sUigGD7e+o5++paZWAcac13LzoPe7Eds+PlqYcrZRF5guwjE726r93gaKWohhXy74rXSRtm9GEujuU78msAnVnydsGqTWRMjzPEA70w4A/SUvt8gNXxXskuq41Vs/oUKwCgWwb97NJvveZU/j3WmLBwrD+Tf5KWbmgfLgbl2sd3pLY6j5XxpCHWg3UdVkj628DCDHrZDeJ+wC5LwYpMq4jVQR8A+Lm5VpH2NXK+xr3y8Dccae8zAStbjAL/zkL+Z6wo1xex/vSus69LHazwc8d1iwMCvbQl0yBuL3TdDuLGisZzpOo6hOXEI23b+wI7sWofO3ExyB6rhnUYeAcrGjgKKxClOtawiR223lWgng+uVZaD1oF/ubt3SBtYrVgvho6sILdgOfxk0qfJc5Xa7GO77Jj9bIQ63dfi4liuBnH7kZaYIAX4L9YzF2bbNNnJhjl5df9f68XnBri48MHAs1jjR84AF7H+GGcCLd3o3A+swuo4TgDWY3XeuvwTy62ej66Lc/5BRwaHMA91I7ESMWd0jq6WThl0r0VuySZYf4RZ7cOj3JIZrpe7fJUTs7BlkZPcCbz8I8zmd4iy/0ycbWmYhfzPHvxex4CYHNiSL87NLh9E5pcEx3IVK9DDZQYZ0jsAd0si8LCPrlV2zi2MNEe6IUNZGaxxfO7s2wzcls21jSYtK1DGxeSWdD5HXxuQ4YJXJy356DGs8UEzsHLMJQIjXeg40u5cwuobmOt0c80B/N0cK0d6Pr4+XzjdeF69UWE5od5Y+ef22+d9BesteRlWiHKm5MFk49yc5NpihYsftPd7GqttfwQuUoe50C9t27AKq+klEcvhxWLV5m7KIJ+lc7NlIrHSj+3ByinoiXO7z0nunXz4DX902v8WD56Hx+3r+htWMMBVrMHDP2HVjMvn0I58c262TGusQcDx9v2wHyvnYzO73J1zC8PKw/kB1svmYfteOG9fg/cy3gvX+Fplm24MGOMkf5eL+3wMVs3V8Zysx3qhD8KqaWZ3bVti/V85ro3Xzs2W8cfqi/8e60UuEet/9zusWrLLdHSFZRH7JH2OnVHgV6zmrNHAaLUGMTvKw7Eyq+9y2nYv1p/1UeAOtTMo2P1Gy7H6DIap6rsZjpUjPUPRR0S6kTZfWUN1M4WQwWAo2BQk5/YGVufm56o6yEOd9VjV9EGq+nmGstZYwQpHsabWSMmtnqHoIyKzsZI9r1fVJr62x2Aw5IwC4dzs1C6HsdqB66mVgic7narAAayqdDl1MU5IRA5iRVe1VGuwZY71DEUfO7HxLqzO+cdUtXBGiRkMhnxJ75ITorEc2wFV3S4it2FlvA/HqkEtUtXYDDqOsPGtrhyUzTosJ9WItKz3OdUzFEHs8Uf+WM3hH2I9E0fJowz4BoPBNxQU5+YYlb9bRCaTeUDvy3Zz0QNODskxfcifuGd/Btnc6BmKJlPJMOAWq781Xya1NBgM14aC4tzC7M87sN6i/4OVGfyUve1/WKmQzpE2uNaR9yyrTO2O3IrOuQtzqpeKiDyGlW6HkJCQ6Dp18jSbj+EaUr58eU6fPo2fnx9BQUFUqVKF8uXLT4uJiZmWvbbBYMgpGzZsOKmq3sxC4hUFxbk5MqUEAB+p6ginsnkichhrOMAgERmj1qSbjtH/3nYa5lQvFVWdhBViTkxMjK5fvz6nuzIYDIZiiYhk1XqWawpK+i3njNOZOvFVdT3WJH5+pM2W7ElGekeZq0z03uoZDAaDoZBQUJxbnNP3fW5kHNsdE3Y6dKplsd9IF/vPqZ7BYDAYCgkFxbltdPoe7kbGMcOtoz/MMa18fTuDtSuaZJDNjZ7BYDAYCgkFwrmp6iGs3IdgJQNNh51JvLG9ut7WOYDlFEuQfvZkh05rrHmcjpI2xUOO9QwGg8FQeCgQzs3mX/bnyyLimPAQEQkCJmBN6b6B9A7nDfvzLRG50UmnElaEJcCbLrKM5FTPYDAYDIWAghItiarOF5H/YCXP/UVEfsEaCtAUa+6iQ0B/dUqpoqqzRGQC8CTwm4gsxUqW6pgK4mus6RwyHitHegaDwWAoHBQY5wagqiNE5CdgKFZ2kFJYA6rfxqpJnXChM0REVmNNVdMaa5zcDuATYIK72ldO9XLCuXPnOH78OFevXs1e2GDIZwIDA6lUqRKhoaG+NsVwDVFVdp7eyU3lb8JpHtUiS4HILVmYyW6c27lz5zh27BgREREEBwcXi5vKUHBRVS5dusShQ4eoXLmycXDFjF2nd7H15FZKlyhNx2odfWqLiGxQ1Zj82n+BqrkVRY4fP05ERASlSpXytSkGAyJCqVKliIiI4PDhw8a5FTNql69N7fK1SSkG4QQFKaCkSHL16lWCg92NODAYfENwcLBpJi/G+EnR/+sv+mdYADBNkYaChrknDUUd49wMBoOhiLP5+GaSUpJ8bcY1xTg3g8FgKMIcPH+QBxY+QKdZnXhv43sUlyBC49wMOWLy5MmICAkJCW5lZs6cyd13301ERASlS5cmOjqar7766hpaaTAYZu+eDcCJSyfYEb+j2DRJG+dmyDfefvttSpcuzTvvvMO8efNo27Yt999/P+PHj/e1aQZDsaF8yfJUCLZS895b614fW3PtMEMBDF6TnJxMYmJitnLz58+nQoUKqevt2rXj8OHDvP322wwdOjQ/TTQYDDYD6w+kf93+rDywktaRrX1tzjXD1NwM2TJ48GBiYmL4+uuvqV+/PkFBQRw7diyT3NixYwkKCmLevHkA6Rybg0aNGnH8+PF8t9lgMKQR6BdIh2odCPArPvWZ4nOmhlwRFxfHs88+y8svv0zlypX55Zdf0pWPHj2aN954g2+++YbOnTu73c9PP/1EvXr18ttcg8FQzDHOzQdEPb/A1yYQ92ZXr+RPnTrF0qVLufVWa8KGQ4cOpZa9+OKLjB8/nu+++442bdq43ceyZcv45ptv+OSTT3Jks8Fg8JwUTSkWg7XdUXzP3OAVERERqY7NmeHDh/O///2PxYsXZ+nY4uLiuP/+++nRoweDBw/OP0MNBgMXrl6gy+wujF03lj/P/elrc3yCcW4Gj6hcubLL7bNnzyY6OpqmTZu61Y2Pj+fOO+/khhtu4Msvv8wvEw0Gg803e77hyIUjfL7tc4YtH1ZsxrY5Y5olfYC3TYIFAXdjY7799lu6devGwIED+fLLL/HzS/++dPHiRbp160ZiYiILFiwgJCTkWphrMBRrfjr8U+r3+2rfV2zGtjljam6GXHHzzTezcOFCvv32W5544ol0ZUlJSdx3333s3r2bhQsXUqlSJR9ZaTAUL95r9x7/a/8/2kW2454b7/G1OT7B1NwMuaZp06Z8++23dOnShdDQUP7zn/8AMGTIEL777jveffdd4uPj+fnnn1N1GjVqRMmSJX1lssFQpPETP26veju3V73d16b4DOPcDHnCHXfcwZw5c+jRowdlypThlVdeYcmSJQA8/fTTmeT37dtHVFTUNbbSYDAUF8xM3Lkku5m4t2/fTt26da+hRQaDZ5h70+BL8nsmbtPnZjAYDEWAoxeO8uKqF9l+aruvTSkQGOdmMBgMRYDPtn7G/L3z6fNtH8atH+drc3yOcW4Gg8FQyElITEid2gagSZUmPrSmYGCcm8FgMBRySpcozed3fk6nap2oF16P2yOKb5SkAxMtaTAYDEWAOmF1GNdmHInJicVy0HZGTM3NYDAYihAl/Ev42oQCgXFuBoPBYChyGOdmMBgMhZQFexfww/4fimVi5Ozwqs9NRMoBbYBGQGWgHHAaOA5sBFaq6pk8ttFgMBgMGbh49SJvrX2L01dO0yC8AePajOP60tf72qwCQ7Y1NxHxF5H7RGQ5cBKYDYwEHgP6AI/b63OAkyLyg4j0FhH/fLTbcA0ZNWoUFSpUcFk2ePBgYmLyLcmAR0RFRTF58uQc68fFxSEiiAhr1qzJVD569GhExOt0YYmJiYwaNYrNmze7PN63336bY5uzIyYmxsybV8SZumMqp6+cBiD+cjwVgl0/o8WVLGtuItIfeBOoCgiWc/sZ2AbEA+eAUCAcqAc0x6rZtQYOiMjzqjotv4w3GPKS0qVL89VXX9GyZct026dPn07p0qW93l9iYiKvvvoqUVFRLid6NRhyQ88bexJ/OZ7pO6bzxC1PmECSDLh1biKyBstZnQDeBT5T1V+z26GI3AoMBvoDU0RkqKq2zFrLYPA93bt3Z9asWbz77rv4+1sND7/99hvbt2+nT58+xMbG+thCgyGN8OBwnm3yLIPqDSI8ONzX5hQ4smqWrAkMB25Q1eGeODYAVd2sqsOASODv9n4MxYj9+/fTr18/wsLCKFWqFJ07d2bnzp2p5StWrEBE+P3339PptWnTht69e6euO5o8v//+exo2bEhISAitWrVi69atWR5/3rx5REdHExISQvny5WnWrBkrV67M1u4ePXpw/vx5li9fnrpt2rRptGrVioiIiEzy8fHxPP7441SuXJmgoCBuu+02fvnll9TyMmXKAPDggw+mNnvGxcWlll+8eJHHH3+csmXLUrVqVV555RVSUlLSHeOHH36gWbNmBAUFUblyZYYMGUJCQkI6md9//52WLVsSFBRE3bp1mTdvXrbnaig6VA6pTICfGbKckSydm6q+q6qJOdmxqiaq6v8BNXJmmqGgkZSUlGnJGKUVHx9Pq1at2LlzJxMnTmTGjBlcuHCBDh06cOnSJa+PuX//fkaMGME///lPvvrqK44fP06fPn3SHTcuLi61f+mPP/6gd+/etGvXjvnz5zNlyhS6detGfHx8tscKCQmhW7dufPXVV6nbpk2bRv/+/TPJXrlyhQ4dOvD9998zduxYvv76aypWrEiHDh04evQoYDkmgJEjRxIbG0tsbCzXXXdd6j6effZZSpcuzaxZsxgwYACvvfYas2bNSi3ftm0bXbp0oUKFCsyePZtXX32VqVOnpnsBuHTpEp07dyYhIYGpU6cycuRIhg0bxv79+z28wgZD0cStu1fVC3lxAFW9mBf7KXIsfwNWvumZbONBcPd76bfN+xts/Mwz/dbPQ9sXvLMvA6dOnSIwMNBlWXR0dOr3d955hwsXLrB582bCwsIAaNmyJVFRUXzyySc89dRTXh03Pj6eNWvWUKtWLQBSUlLo2bMnO3fupE6dOpnkN23aRJkyZRg7dmzqtrvuusvj4/Xr14+HH36YCRMmsHnzZvbv30/v3r158830v9WXX37J77//ztatW1Nt69ChAzfddBPjxo1j7NixNGli5ferWbMmzZs3z3SsO+64g3HjrAS3HTt2ZNGiRcyZM4c+ffoA8Nprr1GtWjXmzZuX2kwaFhZG3759iY2NpUWLFnz66accP36cX375hapVqwJWgE2rVq08PmdD4eHohaMs3LeQ++veT0l/M9lvVphxbgaPKFu2LOvWrcu0dOvWLZ3c0qVL6dixI6Ghoam1uzJlyhAdHU1W8965IyoqKtV5ANSrVw+AgwcPupS/+eabOXv2LIMGDWLJkiVcuODdO9pdd91FcnIyixcvZtq0abRv395lpOjSpUuJjo6mevXqqecJ0Lp1a4/Ps1OnTunW69Wrl+681q5dS8+ePVMdG8C9995LQEAAq1evTpWJjo5OdWxgvUxUqlTJ85M2FBrGbxrP2xve5u65d7PmUObIXkMaHjfUikgYEAXEqWq80/brgDeAW4A4YJSn/XOGwkNAQIDLkP/w8HCOHDmSun7y5El+/vlnpk+fnkm2ffv2Xh+3XLly6dZLlLAiwi5fvuxS/qabbuKbb77hzTff5K677iIwMJCePXvy7rvvUrFixWyPV7JkSe655x6mTp3KqlWrGDNmjEs5x3m6qs3WrOlZN7Orc3M+ryNHjlC5cuV0Mv7+/oSHh6c2sx49etSlIzPOreixM34n8/+YD8DhC4dNP1s2eHN1XsAKMGmMNQwAESkBrMZyeoLl4FqLyM2qeihvTS1itH0hd02Fd7+XuamyABAWFsbdd9/NSy+9lKnMEWARFBQEWKHyzsTHx7sdT+cNXbt2pWvXrpw9e5YFCxYwbNgwhg4dyrRpno1K6devH926dUt1jK4ICwsjJiaGCRMmZCorWTJvmouuu+46jh8/nm5bcnIyp06dSm3yrVKlCjt27Mikm1HPUPipUa4GLzR7gQmbJ9CwYkOaXdfM1yYVaLxxbm2BfRlqZX2B6sAKrPFw3YGngL9iOUNDMaN9+/bMmDGD+vXrExwc7FLG0YS2fft2GjduDMCBAwfYuXMntWvXzjNbypYty/3338/KlSu9CuPv2LEj9957L3Xq1KFs2bIuZdq3b8+SJUu44YYb3NaSsqtlZkezZs2YO3cur7/+emrT5Jw5c0hKSkrtU2vSpAlTpkzh4MGDqdd1zZo1xrkVQQL9Aulfpz/danTj4lUTypAd3ji3qsDmDNu6AQo8oqp7gSUichdwJ8a5FUuGDx/Ol19+Sbt27Rg6dCgREREcO3aMlStX0qpVK/r370/VqlVp0qQJL730EqVKlSIlJYXXX389tTaSGz744ANiY2Pp0qUL119/Pbt372bmzJkMHDjQ430EBAQwY8aMLGUGDhzIxIkTadOmDf/4xz+oUaMGp06dYu3atVSpUoVnnnmGEiVKUL16dWbMmEGDBg0ICgqiYcOGHtsxcuRIGjVqxD333MOTTz7JwYMHee655+jcuTMtWrQArGEGY8aMoVzgP84AACAASURBVGvXrowaNYpLly7x0ksv5UkN2FAwKVOiDGVKlPG1GQUebwJKymNlKHGmBbDLdmwONmGNcTMUQypUqMDPP/9MnTp1eOaZZ+jUqRPPPvssZ8+eTffHPnXqVG644QYGDBjAiy++yMsvv8xNN92U6+M3bNiQEydOMHz4cDp16sSYMWN49NFHeeutt3K9b2eCgoJYvnw5HTt25JVXXqFTp048/fTT7N69m6ZNm6bKTZw4kZMnT9KhQweaNGnC4cOHPT5G/fr1WbhwIcePH6dXr16MHDmS/v37pxsuUKpUKRYvXkxISAj9+vXj1VdfZdy4cVSrVi1Pz9dgKGyIp9mkReQMEKuqd9rrkcCfwCeq+oiT3BSgh6p6n6+oEBITE6NZRcdt376dunXrXkOLDAbPMPdmwefM5TP8d/N/eeKWJ4pc7kgR2aCq+ZaY1pua2w6glR01CXA/VpPkjxnkqgLH8sA2g8FgKNaM3zSe6Tun031ud+b9YTLPeIM3zu0LIARYKyIzgNeABOAbh4CIlMSKptzpcg8Gg8Fg8IgD5w4wc9dMABKuJlA6sFg0huUZ3ji3CcBUrHRavYFE4FFVPesk0x3LAWafyM9gMBgMbokMjWRih4lEhUbRKqIVbSPb+tqkQoXH0ZKqmgIMEJGXsCYq3aaq5zKI7QXuA8zQeYPBYMglt0Xcxpy753Au8Rwi4mtzChVZTXnTD1igquedt6vqPmCfKx1V3Yg1I7fBYDAY8oBA/0AzpU0OyKpZcipwXES+FZGHRST73EUGg8FgyBWeRrAbsiYr5/Yc1qDtO4FJwGERWS4ifxORG66JdQaDwVCMWHd0HQMWDmDX6V2+NqXQ49a5qepYVW2BFdr/N6yQ/5bA/wH7RGS9iLwgIpnnHTEYDAaDV1xJvsJrsa+x5cQW+s7vy4K9C3xtUqEm22hJVT2iqu+ranusQJKHgAVAPeBfwFYR2S4iY0Qk3wbkGQwGQ1Fmy4ktHEqw8s2XDChJTGXzd5obvJrPTVVPq+pkVb0bqAj0A2YC1wMvAr+IyJ8i8o6ItBYT3mMwGAwe0aRKE2Z1n0XjSo0Z1ngYlUMqZ69kcEuOJytV1QuqOkNV+2E5uu7AZCAYeBr4AfhnXhhp8B0iku2yYsWKbPfz/PPPp5tQ0xMuX76MiPDRRx+lbmvevDkDBgzw9jTyHFe2GQy5pUa5Gnza5VP63NTH16YUevJktjtVTcRqqlwgIn5Aa6AnYObdKOQ4TxVz6dIl2rVrx8iRI+natWvqdsfs2Fnx1FNP0bdv33yx0WAoSvhJjuscBifyfCpXe7D3cnsxFHKaN2+e+j0hIQGwZpp23u4JkZGRREaaySIMBmcOnj9I2ZJlzRQ2+UCOXhFEpIqINBaR29wteW2ooWBz6tQpBg8ezHXXXUdQUBDVqlXjqaeeSi131Sy5e/duunfvTpkyZQgNDaVnz57s2+cyP4BbTp8+TdOmTYmJiSE+Ph6wnPCQIUOoVKkSwcHBNGvWjOXL0961nnvuOapVq5ZpPNGsWbPw8/PjwIEDAMyePZtGjRpRqlQpwsLCaNGiBT/99JNbWzZt2kTFihV55JFHSExMpGLFii6n2mnWrBn333+/V+dpKHpcTbnK8BXD6TWvF7GHPZ9M1+AZXjk3EblPRLYDh4B1wCo3S8aZAgxFnKFDh7J+/Xree+89Fi9ezJgxY7IcjOpo4ty7dy+ffPIJH3/8Mdu2baNNmzacPXvWrZ4zJ06coG3btgQEBLBs2bLUyU4HDRrElClTGDVqFLNnz6ZSpUp07tyZtWvXAtCvXz/279/Pzz//nG5/M2bM4LbbbiMyMpJt27bRr18/7rzzThYsWMAXX3xBly5dOH36tEtb1q5dS7t27ejbty8ffvghJUqUYMCAAUyePDmd3Pbt21m7di0PPvigR+doKLp8uOVDtsdv5+iFowz9YSinLp3ytUlFC1X1aMGKjEwGUoDTWJOSunNuqzzdb2FfoqOjNSu2bdvmcvv7m97XBpMbaIPJDfT9Te9nKv/32n+nlk/+fXKm8lfWvJJaPmPnjEzlI1aOSC3/9o9vs7TRU86fP6+Afvrpp5nKatasqZMmTXKr+9xzz2lERETq+jvvvKOBgYG6f//+1G1//PGH+vv769tvv62qqpcuXVJAP/zww1SZZs2a6V/+8hc9fPiw1qtXT1u3bq3nz59PLd+0aZMCOm3atNRtSUlJeuONN+rdd9+duq127dr69NNPp64nJCRoqVKldPz48aqq+sUXX+j111/v9nycbVu1apWWKVNG//GPf6ST+e233xTQn376KXXbiBEjNDIyUpOTk93u+1rh7t40XBsW71usrb5q5fYZL+oA6zUf/5u9qbm9aH8+DVRU1Uaqeru7JbdOV0ReFxG1l39kIXe/iKwSkbMikmAPLn/KDmzJav850jO45tZbb+WNN95g4sSJ7NmzJ1v5tWvX0rx583T9cDVq1KBJkyasXr06S91Dhw7RunVrIiMjWbhwIaVLp00FsnbtWvz9/enVq1fqNn9/f3r37p1uv3379mXmzJmkpKQAMH/+fC5fvkzv3r0Ba0bvI0eO8Mgjj7B06VIuXrzo0pbly5fTuXNnnn76acaOHZuurEGDBjRt2jS19pacnMyXX37JoEGD8PMzt1lxp1NUJ+b2mMujNz/KgLq+jwAuanjzhNUCflLV8aqalF8GAYhIE+BZrMlQs5J7H5gCxGDVGL8HagP/BWaJiH9e6hncM2nSJLp06cLLL79MrVq1qFOnDnPmzHErf+TIESpXzjyOp3Llyql9Z+7YsmULu3fvZtCgQQQHB2fab/ny5QkMDMy0X+cmxX79+nH48OFUhzd9+nTatGlDlSpVAMu5zZkzh+3bt9O5c2cqVKjAwIEDM9m2aNEi/Pz8eOCBB1za+vDDDzN9+nQuXbrEokWLOHr0KIMHD87y/AzFhwrBFfhb47/h72f+cvIcT6t4WP1sU/OzGmkfpySw1T7eXCwH9w8XcvfaZUeAWk7bKwPb7LKn80rP3ZLTZsnCSFbNkg5SUlJ006ZN2qdPH/X399fdu3erauZmyf79++sdd9yRSb958+baq1cvVc26WXLkyJFaokQJXbRoUTr9Dz74QP39/TUxMTHd9ueff17DwsLSbWvQoIEOGTJEz549q0FBQfrBBx+4PKfTp0/r559/rmFhYTpo0KB0tv33v//V9u3ba2RkpP7555+ZdM+dO6chISE6ZcoU7d27t7Zu3drNlbv2FKV701D4oAA1Sy4BmuTMhXrFa1ipvZ4AsooseMH+fE5Vdzs2quox4El79XkXzYw51TN4gIhw66238uabb5KcnMyuXa4TwDZr1ozY2FgOHTqUui0uLo5169bRqlWrbI8zevRonnzySXr16sWqVatStzdt2pTk5GTmzp2bui05OZnZs2dn2m+/fv2YNWsWc+bMISkpiXvvvdflscqVK8cDDzxAt27d2LZtW7qykiVL8s033xAZGUmHDh04evRouvIyZcpw33338e677zJv3jwTSFKMcQSOHEk44mtTigeeekHgBuAY8Bbgnx+eFmgGJAFT7PXJuKi5YSVzVuAKEOxmXwdtmdtyq5fVYmpuFk2bNtW3335bFy9erIsWLdIePXpoaGioHj16VFUz19wuXryoVatW1QYNGujMmTN1xowZWqdOHb3hhhv07Nmzqpp1zU3VqiU+/PDDGhoaquvXr0+V6dWrl5YtW1YnTJig3333nXbv3l0DAwN17dq16WzevXu3Anrddddply5d0pW9++67+vDDD+u0adN05cqV+sEHH2hoaKg+99xzLm07c+aMNmrUSG+++WY9depUun2tWrVKAS1TpowmJCR4dc3zk6J0bxZ0kpKTdNDCQdpgcgNtMbWFrjyw0tcm+RzyuebmzUzc+0WkFfAN0EtEltmOIMWN/Oue7htARIKAz4B4rKCVrGhkf25V1UtuZNYBEbasY3BSTvUM2dCiRQs+/vhj4uLiCAwMpHHjxixevNhlvxpAcHAwP/zwA8888wyDBw9GRGjfvj3vvPMOoaGhHh1TRJg0aRIXLlygc+fO/Pjjj9SrV4/PPvuMESNG8NJLL3H+/HluueUWFi1aRJMm6RsebrzxRqKjo9mwYQNvvPFGurJbb72VhQsXMmzYME6fPs3111/PX//6V0aNGuXSlrJly7JkyRJat27NnXfeydKlSylTxhqY26pVK8LDw7nnnnsICQnx6NwMRYvNJzaz+fhmAC5cvUDpwNLZaBhyjadeEBDgXeAqlkNLwRoakHFJAZK99bLAOKwaU1+nbZNxXXP7m719bhb7e9eW+U9u9bJailPNzZAzNmzYoICuXr3a16akw9yb15aNxzZqx5kd9X+b/udrUwoEFJSaG/A8MBSr2XABsAdI8MaRusPOaDIM+FpVp3ug4njtuZCFjMM257w2OdVLh4g8BjwGcMMNZt5Wg2tOnDjBrl27eOGFF4iOjqZly5a+NsngQxpVasSsu2dRKqCUr00pFnjj3B4GLgKtVHVzXhkgIsHAp8A5YIinavant/Ox51QvHao6CWt2cmJiYsyc8AaXzJ49myFDhlC/fn2mTJnia3MMBYDQEp41uRtyjzcRgRHAj3np2GxexxpjNlxVPQ0jOm9/ZtVw7Sg777Qtp3oGg9c88cQTpKSk8Ntvv9GwYUNfm2O4xiyKW8TFq64H/xvyH29qboewam55TU+sfrpBIjIoQ1kd+/NJEekG7FHVR4A4e3u1LPbrSH0R57Qtp3oGg8HgMcv+XMaIlSO4sdyNvNPmHaLKRvnapGKHNzW36UBrEcmPcC/HHHAZF0eoXQ173THv+ib7s77drOmKJhlkc6NnMBgMHnHy0klGrhkJwJ4ze3h/8/s+tqh44o1zGw3sAuaJSM28MkBVo1RVXC1YQwMARtjbbrV1DgAbgRLAfRn3KSKtsca0HQVS55LIqZ7BYDB4SoXgCoxoMoISfiWIKB3ByOYjfW1SscSbZsl5WJGSbYHtIrIX9+PcVFU754F9WfEGMBN4S0R+UtU9ACJSCfifLfOmWpOn5oWewWAweESvWr2oG1YXgLIly/rYmuKJN86tQwa92vbiinyPIFTVWSIyAStl1m8ishRrDF57IBT4GisRcp7oGQwGgzfUDa/raxOKNd44t475ZkUOUdUhIrIaeAqrT84f2AF8AkxwV/vKqZ7BYDC4YtupbdQJq4OfSUlbYPD4l1DVZd4seWGcqg62+9r+k4XMVFVtqaqhqhqiqtGq+n52DiqnegbPOXLkCHfddRdly5ZFRFixYkWO9/Xhhx9Sq1YtgoKCiI6OZtkyz26xNWvW0KxZM4KDg6levTrvvfdeJhkRybQ0b948x7YaihdbTmzhge8e4OnlT3Phalb5IQzXEvOaYcg3/vWvf/Hrr7/y1VdfERsbS+PGjXO0n2nTpvHEE08wcOBAFi5cSP369enWrRu///57lnp79uyhc+fOVK9enQULFvD4448zfPhwPvroo0yyf//734mNjU1dPv744xzZaihexF+O5+nlT5OYksiKAyt4ec3LvjbJ4CA/c3sVh8XklnRP+/bttWfPnrneT+3atfXBBx9MXU9OTtYGDRqkzg7gjscee0xr1aqlV69eTd325JNPatWqVTUlJSV1G6Djx4/PtZ2FjeJ8b+YVySnJ+s76d7TB5Aba8quWuv/sfl+bVGjAV/O5iciPds7HHCMiLUXkx9zsw+B7Bg8eTExMDAsWLKBevXqUKlWKrl27Eh8fz549e2jbti0hISHExMSwZcsWwGrqW7ZsGXPnzkVEiIqKytGx9+7dy65du+jTp0/qNj8/P+677z4WLlyYpe7ChQvp1asXAQFpXcv9+vXj4MGD2db6DAZP8BM/hkUP46XmLzGu9TgiQyOzVzJcE7JqlqwDrBKR70Wkr4iU9GSHIlJSRPrbUYg/4j6i0lCI2L9/Py+//DJjxoxh0qRJ/PTTTzz22GP069cvddLPpKQk+vXrh6oSGxtLo0aNaNu2LbGxsamTh6oqSUlJ2S4OduzYAUCdOnXS2VO3bl3i4+M5ceKES3svXLjAgQMHXOo579fBqFGjCAgIoEKFCjz00EPEx8fn7oIZihV9bupDs+ua+doMgxNZRUvWAl7FSmbcDjgvImuwBjdvB05hJTsOBcKxZs9uAbTEys+YhDV9zKh8st1wDYmPjyc2NpaaNa3x+1u2bGHs2LF89tlnDBw4ELAcV9euXdmxYwfNmzcnNDSUsLCwdMEZn332mUezUVutFnD69GnAmg3bmfLly6eWV6xYMZP+mTNnstVzMGjQILp3707FihVZv349o0eP5tdff2Xt2rX4+/tna6uheLH11Fbqh9f3tRmGbHDr3FT1LDBMRN7DmupmEHAn0CWL/QnWZKPjgP+palzemVq0WLFiRa6iB3NKmzZtaNOmjdd6UVFRqY4NrIk+Adq1a5dp26FDh1JrSBnp3r0769at8/r4IpJu3eH8Mm7PTs/V9smTJ6d+v+OOO6hbty533XUX8+fP55577vHaVkPRZe7uubz808s8cvMjDG001IT+F2CyHeemqnuBZ0TkRawxYW2AW4FKQFngDHAcK63VcmCVql7JL4OLCjl1Mr4iYw2oRIkSmbY7tl2+fNntfsLCwihb1vOMDY6a1pkzZ9LpuauZZbTXIefAXU3QmS5dulC6dGk2btxonJshlQ3HNvBq7KsAfPTbR5QrWY5B9TPmejcUFDwexK2ql4BF9mIw5AhvmyUdfWY7duygWrW0yRx27NhBWFiYyyZJgJCQECIjIzP1rbnrw3PGUavLrlZoKF7UD6/P7RG3s+LgCuqE1aF37d6+NsmQBd5kKDEYco23zZI1atSgdu3azJw5k86drXSlKSkpzJw5kzvvvDNL3TvvvJO5c+cyZsyY1L6z6dOnExkZSYMGDdzqLVq0iISEBKKjoz2201D0CQoI4u22bzN+43gG1BtASGB+TJBiyCuMczNcU8LDwwkPD/dKZ9SoUQwYMICoqChatmzJZ599xu7du5k6dWqqzMqVK2nfvj3Lli2jdevWAIwYMYIpU6bwwAMP8Oijj7Ju3To++OADJkyYkFormzRpEuvXr6dDhw5UqFCBjRs3MmbMGJo2bUrXrl3z7sQNRYJAv0CGxwz3tRkGDzDOzVDg6d+/PwkJCbz11luMHj2a+vXr8+2336arfakqycnJqc2ZYAW4LFq0iOHDh3PnnXdSpUoVxo0bxyOPPJIqU7NmTT777DNmz57NuXPnqFKlCgMHDmT06NEmUrKYc/TCUeLOxdH8OpOKrTAizn8GBu+JiYnR9evXuy3fvn2728hBg8GXmHvTPScvneTBRQ9yMOEgb97+Jp2j8nsGr+KHiGxQ1ZjsJXOGiWM1GAyGDLwa+ypx5+JISknin6v/yYmLrpMFGAouxrkZDAZDBl5q/hJRoVH4iz9v3P4GFUu5jso1FFxMn5vBYDBkoFKpSnzS+RO2ntpKm8g2vjbHkAOMczMYDAYXVCxVkTal2vjaDEMO8bhZUkQeEZHg/DSmqGKCdgwFDXNPpnH2yln+uuyv7D2719emGPIQb/rcJgEHRWSciNTKL4OKGoGBgVy6dMnXZhgM6bh06RKBgYG+NsPnnLp0iocXP8zKgyt5aNFD/HHmD1+bZMgjvHFu32LNAPAMsF1EFolIdzE5irKkUqVKHDp0iIsXL5q3ZYPPUVUuXrzIoUOHqFSpkq/N8TlHLx7lYMJBAE5dPsXWU1t9bJEhr/BqnJuIRAJPAg9hJU5W4AAwEfhYVYtdvGx249wAzp07x/Hjx7l69eo1sspgcE9gYCCVKlUiNDTU16YUCDYc28BTy57ihaYv0OPGHr42p9iQ3+PccjSIW0QCgT5Yc721wHJyicAsrKluYvPSyIKMJ87NYDAUbE5fPk35oPK+NqNYUSAHcavqVVWdoqotgUbAx1iTk94PrBaRDSLykKezdxsMBsO1YO+ZvZxLPJdpu3FsRY9cD+JW1V+xZuz+FGuyUsFyeB8CcSLycG6PYTAYDLll0/FNPLDwAYYsHcLFqxd9bY4hn8mVcxORDiIyB9gHPAVcBj4B+gPfYfXLTRKRv+XWUIPBYMgppy6d4vHvH+dc4jl+PfErz6963tcmGfIZr52biJQVkWEisgNYDNwDHAZeBKqq6iOqOl1VuwO3ARcA49wMBoPPCA8O5+/RfwcgLCiMx2953McWGfIbjzOUiEhjrACSfkAwVvPjCmA88I2qpmTUUdVfRGQBcG+eWGswGAw5pG+dviRpEndE3EFkaKSvzTHkM96k33KEBF4EPgLGq+rvHuhd8PI4BoPBkCuSU5JJIYVAv/QD1f9S9y8+sshwrfGmWTIOGIHV9Pi4h44N4FHApEIwGAzXhEtJl/j7yr8z6qdRJnFCMcabGlVNzcGdYuske6tnMBgM3nLx6kUeWvxQaqaRCsEVeCb6GR9bZfAF3tTcFovI8OyEROQZEVmSC5sMBoMhRwQHBNOgQoPU9aSUJFN7K6Z4U3PrABz0QK4e0D5n5hgMBkPOERGeb/o8Ry8c5faI2+lbp6+vTTL4iPwI9CgBZIqcNBgMhrzGUStzzt8e4BfA+HbjMTndize5zlDijD1DQDRwMi/3azAYDBm5mnKVf/3yLz7+/eNMZcaxGbKsubnoO+uURX9aAFALuB4rgbLBYDDkC+cSzzF02VA2Ht8IQGSZSDpHdfaxVYaCRHbNkh2cviuW47o+G50twLO5McpgMBiyolRAKUr4l0hdX3VwlXFuhnRk59w62p8CLMFKt/UfN7KJwCFVNXO1GwyGfCXAL4Cxd4zlL9/9hV61evFQg4d8bZKhgJGlc1PVZY7vIrIGWOm8zWAwGK4FqpqpH61cUDlm3z2boIAgH1llKMh4HFCiqrer6pv5aYzBYDBkJP5yPI8ueZQlcZm7+41jM7jD5Hw0GAwFlt2nd/PE0ic4fvE4W05uIapsFLXL1/a1WYZCgFvnJiIv2l8nqOppp3WPUNXXc2WZwWAo9lxf+npKBZQC4HLSZTYd22Scm8EjxF1qGhFJwYqQrKuqu5zWs90nVkpJ/7wzs+ASExOj69evz17QYDDkiJ3xO3lq2VO80uIVbq96u6/NMeQRIrJBVWPya/9ZNUu+juXMTmZYNxgMhnzh7JWzlC1ZNt22m8JuYmGvhQT6m8lFDJ7j1rmp6sis1g0GgyGvUFWm7pjK+E3j+aTzJ9QLr5eu3Di23JGclMS+I8fZeCyZzQfOoApv9LrZ12blKyagxGAw+Jzxm8bz4W8fAvD3FX9nRvcZlClRxsdWFV5OHj3A/l9XcOXPtYSe/JVqV3axPrkZzyc9BkBICX/G3NMAf7+im6bMODeDweBzetbqydQdU7lw9QJlSpRJ/TRkj6akcHDvVo78thz+jOW6s5uJ1MNUcBYSuMXvj9TVC4nJ7DmewE1Viu419ti5iciTwLtAT1Vd4EamGzAHGKKqH+WNiQaDoagTWSaSfzb7J9tObWNY9DBK+pf0tUkFluQUZdvhc6yNi+f09pUMPvQKkZwhMhu9ML8LdK4TRsNqFbmlajluCCt1Tez1FW6jJTMJinwP3Axcr6oup7QREX/gELBZVbvkmZUFGBMtaTB4x96zezl9+TTRlaN9bUqhQFNS+HPnRo78upRPLrfj57jTnL+cBEAVTvFz0NBMOokawB8lanM2vBElqjXh+vqtqFy1JuKXpxPB5ApfRktmpA7wmzvHBqCqySLyG9aEpQaDwZCKqjJr9yz+vfbfhASGMPvu2YQHh/varALJ4bidHNywEP8/fyTq3AaiOEMUMPpKGc5rVKrcUcI5kFKRsn4X2BfcgItVmlKuzh1Ub9iKusEhvjK/QOCNc6sIrPRA7jhgBqMYDIZ0XEy6yMTNE7mcfJnLyZf51y//4u02b/varALB+bPx7Pl5Pok7l1L19C9E6DGX06/c5reVbclRVA4tSdPq4TSNKs/lCguIqHEjt/gXi6HFHuONczsL2TbrAkQACTkzx2AwFFVCAkN4teWrPLn0SW4sdyNP3PKEr03yGarK9iPnWbHrODeufYW2FxbSSJLdyp8lhD9CGnPbTc3p17w1NSuGmAlZs8Eb57YJaCsiNVX1D1cCIlITuA34MS+MMxgMhZeklCQC/NL/xbSKaMXYO8bSJrJNsUt6fPb0STZv3cGCo2VYuesEx85dAWB4QACdAtI7totakj3BN3MxoiXhN3ekRoMWNA4wwe3e4M3Vmgx0Ar4WkV6qutu5UERuBOYC/raswWAopiz7cxlj14/lg44fUC20WrqyLtWLRawZmpJC3I4NHF03l3IHV1ArcTtltTozEkenk1uRfCt/C/iaPf41OVHldso26MKNjdvSsGTxcv55jTfRkgLMB+4CkoDVwA67+CasfrYAYJGq3pX3phZMTLSkwZCe9ze/z8RfJwJwa8VbmdxlMv5+xaM/KCnxCjvWLSHh1/lEnlhBhB5LV56iQsyVCcQTStngQG6vVYG2tcJoHelHhSo3+Mhq31BgoiVVVUWkF/AO8CjQxl4cJAETgOF5aJ/BYChktI1sy0dbPiJJkziccJjDCYeJDPWku75wcv7iJXatmIru+I5a52JpwAW3srv8qjOsWRnqN76NW6qWJcC/4ITmFzW8asRV1UTgKREZDbQHHO0NfwLLVPVoHttnMBgKGfXC6/HYLY9x4NwBnmv6XKZEyEWBw2cusWz7MZZsO8bavSdYHTCKinIuk1yCBrOrTFOSa3WhRvO7mTbhI0b17HbtDS6G5KiH0nZiU/LKCBEJBO7AavJsieU0w4ETQCzwX1VdkYX+/cCTQEOsPr8dwKdYc9G5HZeXUz2DwQApmsKsXbOoUbYGMVXSty490fCJIhfN9+fOzRyOnc53p6rwxYkbnUqEH6QxfQNWAHCMcOIqtKbUvAdcpAAAIABJREFUzd2p3awLjYOKdiaQgkpBCb9pDXxvfz8KbAAuYA0Gvxe4V0RGq+rLGRVF5H1gCHAZWAZcxapV/hdoLyL3qWqmGNuc6hkMBjh4/iD/XP1PNh7fSLXQaszqPitd9GNRcGyaksK+bes49ssMqhxaQvWU/VQDTic35QuGpZPdWL4zN5SrRcWYe6h5821ULkCZQIorXjs3EbkJ+BtWf1uEvfkQsByrhrXDjWpWpACzgXdVdVWG4/XFqiW+JCLLVXW5U9m9WA7qKHCHI4JTRCrb9vQE/oqVE5Pc6hkMBosS/iXYdXoXAH+e+5OpO6byUIOHfGxV7tGUFPb8upqTa2dS9ehSauhhamSQaeu3mVD/K9xSI4KO9SrTvm5lIsoF+8Reg3u8cm4iMhgraKQE1ozbDkKBusAjIvK4qn7mzX5V9QfgBzdl00WkI/AwMADL+Th4wf58znlogqoesxM9rwCeF5HxGZoZc6pnMBiASqUqMTxmOK///DoPNniQv9T9i69NyjEpKcpvu/dw8YdxRB1bSi1OUMuF3CUtwY7STUmuczer23YktHTRzahfFPBmVoAmwIeAH9Z4tk+AP7CcXHXgIaAX8KGIbFPVdXlo5yb7s6qTPVWBaCARmJlRQVVXisghrNplc+Cn3OgZDMWVk5dOsufMHppf1zzd9ntr3UtM5Riql63uI8tyTnKKsnZfPAt/P8LirUe5cO40G0rOpKQkpZO7oEHsCL0NqXc3dW7vRaPSRS84pqjiTc1tBJZjG6CqX2Uo2wEsFJH+WE2I/wD65o2JAKkvUkectjWyP7eq6iU3euuwnFQj0pxUTvUMhmJFUkoS03dO57+b/ouf+DHvnnnpEh37iV+hcmyaksKeLWs4FTuFCaduZWWC8/CEUqxOuZn2/ps4Rwg7y7YisME91GnVg+hinoC4sOKNc2sFbHDh2FJR1a9EZBhW5GOeICJVgMH26mynIsdT9WcW6vszyOZGz2AoViSlJDFl+xQSrlqpYt/Z8A5jWo3xsVXes3/3rxz68QsiDi6glh6mFrAz6SQreTBVJiykBLuqPUKFG4Kp26IrTUx2kEKPN84tHDf9YhnYDdyaM3PSIyIBwJdAWaxxdPOdikvbn+5HTKYlcHZuHM+pnrNdjwGPAdxwQ/HKKmAoPgQFBPFC0xcYsmwIUaFRdK3R1dcmeczxQ/vYu/xzwvfNo1byHjI+pV39f+H9Eo/QqWEEdzW4jqbVw8yA6iKGN87tNFDTA7katmxeMBErPP8AVjCJM46AFs/yh+VeLxVVnQRMAiv9Vk73YzAUFK4mX2Xj8Y00u65Zuu23V72d/7T+D20j21LCv4SPrPOMM6dPsfOHzwnZOZd6V7ZQSTI/mhc0iG3lWlPy1j6sadWewMBAH1hquBZ449x+AnqISA9V/caVgIh0xwrCmJtbw0TkXawIyaNAexfZT87bn6Vxj6PsvNO2nOoZDEWSHw/+yNh1Yzlw/gCzus/ixvI3pivvHNXZR5Zlz8XEJJZuP868zYfYt+s3lgWOsgqcYrkTNYDfQ5qjDXpTv819NCmV1aNvKCp449zeBnoAM0XkS+AzYB9WDagGMBCrdpViy+YYERmHNZbuBJZj2+1CLM7+rOaizIGjxzjOaVtO9QyGIoeq8vm2z4k7FwfAv9f9mw86flCgB2FfTbzCttVfM+dQOWbuVi4mOnItVGazf01u9fuDZBW2B93CpZt6UbvtX2hcvoJPbTZce7xJnLzaDhZ5GxhkL84IkAwMU9U1OTVIRP6NlXz5FNBRVbe5EXUMD6gvIsFuIh+bZJDNjZ7BUOQQEZ5t8iz3zb+PUgGlaBnREkURCpZzS0lOZse67zm/7itqn1rGLZxnydU+XEy+J53c4nJ9uFwliRvbPECD67N6fzUUdbxNnDxeRFYDjojI67Gc2iFgJVaGkRw7BBF5E2vIwWksx/ZrFrYcEJGNQGP4//buPC6r6078+OcLKKAogoALKsomLlHc465xqcYtm5mkSZO0nWYmSdN00rS/tJ1mXplpG7uN7XTSdtKYPbExcYmaxCzGLS5R4xI3VFAEBdwABWR/zu+Pe1keBCLwbMD3/Xo9ryv3nHvuuYfn4eu9z1lYDLxWp6ypWOPicrDmp2zRcUq1duWOcjZnbmZmv5lOd2aJYYk8N+k5xvUa59TV39uMw0HakS+4uOMNBmRvYDCXnNIX+W/n+cpFxEWGsCg5moXDe9M/ovV0elHu1eTpt+zgVfeurcXslQb+H5CPFdhuJEg+hzUQ+zcissMYk2qXFQX8xc6zpJ5ZRpp7nFKt0saMjSz9cilnrp7h+RnPM6WP82idW2N9ZwnGc6eOkLHlNXplvk+8I5P4evKcpzu5vafxwbxxDOoT4dOPUZV3+MTEySKyEPh3+8dU4PEG3qwpxpglVT8YY94Vkb9izex/SEQ+pWYC5K7AGqyJkJ009zilWqudWTs5c9Ua2rn0y6VM7D3RpxYQvXC1hPVfZfPV3m38Mf/x6klra8snhOPdZ9Jl9D0kjZ1ND3/fqb/yPT4R3IDwWv8ebb/qswVYUnuHMeZR+1HpY1irC1QtXfMSjSxd09zjlGqN/nX4v7IubR3+4s/CuIU4jAN/vBscruRd4tPjuaw6nMvOtMs4DAhhPBUYQR+xHkEWmUCOhU6mQ/JiBk28jXE6uFrdoAaDm4i80IJyjTHmX5qQ+RXglRac7C3gLU8dp5SvunjtIstTlvNI8iN08KsZwxURHMHS6UsZ0n2IVxcPLblWyNEt7yCH32VI4S52VXyH7ZXTqtMNfqxxTGFSl2wqhtzJ4Kl3M1rnc1TN0Nid2z+3oFwD3HBwU0q13KtHXuX5A89TXFFMVKco7km6xyl9Qu8JXqlXRXkZR7evo3T/2wzK38pIsTsoCyzy2847ldMQgXEDwlmUHM3cobPo1sm3B4wr39dYcPuex2qhlGoxQSiusALHXw/+lYVxC+nUwTurQBuHg+N7N3Jl91skXPqUYVytqqSTHh3L+MWMeG4d0Y9eobommnKdBoObMWaZJyuilLpxxpjregjePfBuXj3yKt2CuvHkqCcJDvB8sDieU8DnO7Yx99ATJJkL9ebJlN6c7TOP6Mn3k5CYXO/aaUq1lK90KFFK3YDyynJWp67mzWNv8uocK5BVCQoI4uU5LxMdEu3RnpCZuddYezCLtQeyOH6+gEDKuDuwwOku7QLhnOoxm+4330f88En09dNJipV7NSu4iUgIVo/GSCDDGPOFS2ullKrXE5ueYNu5bQC8dOQlnhz1pFN6v66eWaXi8vmzpG56ndC0tSwpms9mR81CIKV05GPHaGb57yMlbDqdRt/DoHFziQrQ/0srz2nSu01EugB/wJpHsqor1qvAF3b6I8BPgbuMMbtdWE+lFHBb/G3Vwe2T9E94fMTjTr0i3angSi4pm5YTmLKKwcX7GCfWaJmF/uHVwS2ogx8zB/Wge9LvCBwcw9gg73znp9QNBzcR6QRsxlqd+hKwD5hdJ9vHwPPA7YAGN6WaqdxRzr7z1y9BMzNmJuN7jWd87/Hck3SP2wNbSXERR7eugkMrGFywkzFSbiXUeuR4i99+ZiaGMX9EDDMH9yAkUO/QlPc15V34I6zAthx42BhTJCJOA52NMWkichK4xYV1VKrdMMbw7sl3efGrF8kuymbNojXEdoutTvcTP16Y3ZIhqF+v0mHY99UhzObnSMrfzEiuWQl1ejoe6zCEqwm3kTDtPl6Mqm9OEaW8pynB7W4gG/iuMaakkXxngMEtqpVS7ZSIsDlzM1lFWQC8cOgFlkxe8jVHtZwxhoNnr/DegXOs/yobU3CBLwI/xL/Ogp9p/gO4GLOA/tMeYFA/7eeofFdTglsc8NHXBDawHlnq4klKNdO/DPsXtp7dSlhgGIlhiW4915mUfWR9/jr/d3k4m/Mia6WEst0xlCn+hzgnPciInkfvifcTN2gUcW6tkVKu0ZTgVg4E3kC+PkBh86qjVPuQX5LP8pTlHLp0iOdnPO80Zm1Y5DB+N/V3TIme4pZB2DmZqaRvfo3I9HXEVZ4iBjhYMZ/NfLM6T0RIIMfjfkCPhAgSR04jWrvuq1amKcHtBDBCRAKNMaX1ZRCRbsBwdJFPpRpUUlHCgjULyC/NB2BPzh7G9hrrlGdO/zkuPWfuhXOc3PQGXVLXMrj8MD3rpC/038Ff/O/nG0N7syg5mptjwwnw14CmWq+mBLeVwK/t148ayPNLIARrrTSlVD2CAoKYFTOLd05YH5P30t67Lri5gtV1/y0CU1Y7dd2vrdR04EiX8chNd7F72gyCAnVOR9U2NCW4/RlrkdIfishorGAHECMi38Na1XoGcAR40aW1VKqV2n9hP4VlhUzuM9lp/0NDHuLwpcM8OORBvtH/Gy47X0l5JZ+lXGDtgSwyjn/JBwE/txJq9XSsNMKR4JGUJt1B0rR7GdnNd1bfVspVbji42V3/Z2MFtcnAJDtpmv0S4ACwqKHHlkq1F1mFWTy97Wn2X9hPdEg063uvJ8Cv5uPWr2s/VixY4ZJzlZeVcnT7OlZlhbHyRAWFpRV2SjTH/PoyyC8TgJQOg7kSv4j4afczrEcfl5xbKV/VpNGWxphMYKyIzAduBWKxFvnMBD4EVuoin0pBeFA46VfSAThXeI6P0z/m1thbXVa+o7KSlD2fULBnOYmXP2M4V1lXfh+FlfOc8n0QchdXelYSM/V+kmIGuuz8Svm6Zk0lYIxZD6x3cV2UapVyinII9A8kLCisel9QQBCLBy7mpcMvMT92PoO7t3zop3E4SP1qB5d2vUlszkcM5rJT+kL/HbxYOY8BEZ1ZMLw3C4f3Ij5qXgOlKdW2NbYS97vAMmCDMcY0lE+p9upE3gn+/tXf+eTMJ3x76Ld5YuQTTunfGvQtFicupmfnun0Tm+bM8QNkff4Gfc6+T4LJqneJmAuEU9TrZtbNu5mhfcOvWw5HqfamsTu3O7DmiMwWkVeBV4wxJz1TLaV8X8bVDDakbwDgnRPv8PCwh53WUKu9HE1Tnc27xvtfZXN870b+u+DHxNSTJ48unOg+gy6j7yFp7Gyi/D23zI1Svq6x4PZX4J+A3sDTwNMi8jnwEvCOMeaaB+qnlE8oLCskpGOI077pfacTHRLNucJzJIYlcrn4Mn26NL+jRk5mKhvSSllzJJ8DmdYYOD968XRgN6LE+rnIBHG02xQ6Jt/N4IkLGdfxRuZVUKr9aWwl7sdE5N+ARcB3gJnU9JL8s4i8DbxsjNnhkZoq5QXHLh9jecpyPjz9IW/Oe9NpOix/P3+eufkZwoPDSQpPalb558+d4vSWN+l26n2SKo5xoOxRDjgmVac78GONmczokHzM0DsZMnUxYzqFNFKiUgq+pkOJMaYMa0D2OyLSG2uc2wPAQOC7wHdE5ATW3dzrxpgcN9dXKY968dCLfHzmYwDeTnmbX4z/hVP6hOgJTS7zYlY6aVveIvTUegaVH6FHrbR5/rtY45hEgJ8wKSGCeTf14htDZ9M1yDNrtinVVjRlnFsW8BzwnIiMx7qbW4wV6JYAvxKRD7EC3XpjTKUb6quU25Q7yq9bH+3epHurg1valTSMMc3qrHEpJ4O0LW/RJW09SaWHiZTr+2hVGD8iOnfkt7fcxOyhPenWSWcLUaq5mjsUYCewU0R+ANwFPIQ1kHu+/boI101fp5TPMcbwRc4X/CPlH+QU5bB83nKn4DWqxyi+PfTbzOg3g2ERw5oU2C4VlvLh4RzO7FnPzy79nIiqgFZntpCjQckUJywkYeq9jIjsxQhXXZxS7ViLlsw1xhQDrwOvi8gs4A0g0n4p5fMKygt4fOPjlFRaKzkdvHiQ5Kjk6nQR4clRT95webkXzvFRWinrj1xgZ9plHAY60ZOnAgMIwlrFutIIx4KGUxS/gISp93KTLvSplMu1KLiJSAhWj8qHgAnU/J80s2XVUsr1jDGUO8rp6F/zuK9rx67Mi53HypPWVKm7c3Y7BbcbcTErnVPb/kFI2gcklX7F+vKn2e64qTr9GkFsdIykf3AxhfELiJtyL0N79nXNRSml6tWs4CYi04FvY42FC8YKaqXAWqzv3D52VQWVaqncklzeS32PlSdXckfCHXxn6Hec0u9NupeO/h25Z+A9xHaLvaEys9KPk/H5crqd+Yik8qM1jyoE5vntYrvjJkRgTP9w5g/rxZjBq4kK7ezaC1NKNeiGg5uIDMDqLfkg0I+au7QDwMvAG8aYPJfXUKkW+vzc5/z3l/8NwMoTK3loyEP4Sc1aZQPDB/KzcT/72nIyThwga8cKIjI3EF+ZRu968jiMkNC5mP+YMpi5Q3vRMzTIVZehlGqCRoObiHTC6hH5ENYYN7FfecBbwDJjzAE311GpG1ZQVkCXjl2c9s2KmcWSL5ZQUF5AbkkumQWZxHStb84PZ8YYjmUXsOFwNhX73+InxX+kXz35KowfKUHDKIq9lbjJ9zCmdwxjXHQ9SqnmaWxuyWVYga0zVkBzAJ9iPXZcbY+BU8onfJbxGatPrmZ71nbev/19eoX0qk4LDgjmsRGP0blDZ2bHzKZTh04NlmMcDk4c3M6q7O5sOHKeM5etiXh6EM9Pat2ElRl/jnUaRVnCfOIn383QyF4NlKiU8obG7ty+bW9PA69gzUZy1u01UqoZVhxfwfas7QCsTl3No8mPOqXfN+i+Bo+tKC/j+O5PKDiwmgEXP2Mgl9lSuoQzpuY+7Tzh7DRDCAwJxzFwPolTFjNcF/lUymc1FtzeBF4yxmzyVGWU+jp5JXnkluQS1y3Oaf9t8bdVB7f0q+lfW861wisc376GiqPvk3BlO0ModEqf67+blIp+hAQGcEtSFHOH9mR44hY6BepMIUq1Bo3NLfktT1ZEqcZkXM1g6ZdL2Xx2M8mRybw852Wn9On9pvPI8EeYHzuffl3r+2YMLuVkcmr7u3RM3cCga18yQsrrzXeFzgzq1ZVlt4xmYnwEQR10tn2lWpsWjXNTylOCA4LZlLmJSlPJ3vN7ySzIpG+XmrFigf6B1z2KBDh1sZCPj57nk6PnmXLu7zwRsMpKqDPRyAXCOd19Cp2G30bSzbcyW2fbV6pV0+CmfEpOUQ4fnP6AuxLvomvHrtX7IztFMjF6IlvPbmVYxDCulF5xCm5VHJWVnNi3iczDO1iSO4W0i0XVaSUyqia4Aaf9YsjpNYOI0bcTP3wSUX5+15WnlGqdNLgpn/HLXb9kxfEVGAyhHUO5M/FOp/QfjPgBPxr1o+sGWpcUF3F85/uUHl5LbO42ksgnCfh5SSwQVp3vGP3ZGTQJ6TuOvuPvYkDsYAZ44LqUUp6nwU35jH5d+mGwJhd+//T71wW3geEDq/994dxpTu9YReDpT0gs2sdwKb2uvJn++1jlN4spCZHMGtyDGYN6EN55vnsvQinlEzS4KY9xGAcHLhxg/an1ADwz/hmn9LkD5rJ031LG9hzL/FjnIORwGA6ezSdv05/pl7GG+Mo0oqoS63x/lktXUrtN4raRs/nFzbMI7qgdQpRqbzS4KY85feU0D254ELA6gPzbqH9zmk0kslMkm+/eTGhgKAAFxWVsS73MxmMX2HLiApcKy/hlwAFuCUi7ruwMv2iyoqYROmIRiaNmMDZA39pKtWf6F0C5xcm8k8R0jXGagT+uWxxJ4Umk5KZQWlnKlrNbrrtDu5qZwbEvVhOSsZGzxR15tOyHTukbHSO5n42UGX+OBw2jKGYmfcbeRr/4ofVOjaWUap80uCmXWnF8BctTlpOan8r/TP8fpveb7pR+98C7OXb5GPNj55MclUxpyTVO7vmEwsMfEH1hK31NFlV9IOOkI4GUUYoVICNCOtIzYTb7u/Ul/uYF3BQa7uGrU0q1FhrclEtlFWaRmp8KwIb0DdcFt8WJizmbephzm9Zx6MzTJFw7wNB6OoMABEsZt0dlE3XTTGYkRXFTdCh+fgKMdfdlKKVaOQ1uqslyinL4KP0jAB4c8qBT2pwBc1h2eBlB/kEEBVgzDV8rq2Bn2mW2nLjI9uNZrC26nz5VAa1OZ5BrJpDjnUdRHjebAeNvZ0nv/u6+HKVUG6TBTTXJ8dzj3LXuLgDCg8K5b9B9BPjVvI0Ghg3kj1OX0qsgkKsHP+U/9/+dN7L6UFbpqM6zo8MQZvnvq/45U3qTFTGB4CG3kjhuDiOCdVFPpVTLaHBT9TLGcOrKKQaEDnBa2DMhLIGo4CguFF8gtySXPTl7GN97PFfyLpG2ax0VJz7hprydRJELwOmK6ZRVfs+p7M0ylshOgZTGTKfPmPn0jR3C9XONKKVU82lwU9d55fArrDy5kvSr6bw+93WSo5Kr0/zEj9n9Z5Oad5JkYrj24T9IyXmS+LIURorjurKm+h+ECkNSz65MTYxkamIko/rPITBAx54ppdxHg5u6zumrp6uXjfn0zKckRyVjjCH1QiHbTl7CHA5nafZ6QqS45qA6351dpROpIWOoHHALu26ZRs+wEI/VXymlNLi1Q+WOcvbm7GVjxkb6denHA0MecEqf0W8Gq06uIsg/iKy0/TyVtpetaflcKLA6gUTQmaeDip2OcRghtUMCuT0n023YHOJHTGNkh44opZQ3aHBrh7af287jnz0OQGxobHVwK7yaR+qej5Djn/CzIsPt104SZE5we+kYLpiE6uMvEcoxRz/C/IrJDBuLX/wtDBgzl8SoaK9cj1JK1aXBrQ3LLsxm7/m9LIhb4LR/XK9xBPkHUVJZwqkrp1j10iMMyj5IfFkKyVJ5XTkT/Q6zvzKBrkEBTIiLYFJCBJ37bKBH71701GVilFI+SINbG+QwDh748AEOXjwIQHJkMn279qWswsFXZ/PZmXaZYYWhxJRdZs61QkaWvGW9Eep8b1ZmAjgZOIRhScN5b+JEhkaH4u8n151PKaV8jQa3Vq68spwKU0FwQHD1Pj/xo1tgt+qf//f935F17T72pudRXG7dmf2rfxxPd9hzXXlp/rFcjBxP50EziR89iyGduzDE/ZehlFIupcGtldqdvZsVJ1aw/dx2vj/i+9w36D4qystIO/g5uUc3kZC3j+2hhnHFJQzL/pS387/hdPwux2AAzkovssJG4xc3ldgxtxIXFU2cNy5IKaVcSINbK5VRkFE9Bdaa3csYuvYV4q4dYqDdPX+4CN/NhxBjLf65lCtcIpR+4Z24OTacCQOGcL7HfPr0iaOP165CKaXcQ4ObjyquKGbr2a1sO7uNksoSfj/19xSXVbI/I4/d6bnsPwUEWnlLyrIYXJxNh1pfh3Wyg1oOkWR2G8V/jYpj2E3DiO4WfP3JlFKqjdHg5qOKyot4astTAPgZ4dEDy/j4fC8qHMbOEcxPu5dzS9kl4svLq/uCnKc7GaGjMP0n0yd5Nr0HJNHTK1eglFLe0+6Dm4h8E3gEGAb4AynAy8BfjTHXzyflQiUVJezJ2cPWs1v5/ojvU3zhImcPbKQyfQdRefsZ3LOUo4GBOMTguLaBCsdDTseHXUmkS0AFX3ZNxhEzkegRs+ndfxA9tHu+Uqqda9fBTUSeBx4FSoCNQDkwA/hfYIaILDbGXD/wywUqHYaH1j3AkavHAIjd9hr3FuU43WXdfbUz2QElTCkuJvtaBBuAxB4hjOkfztgB4Yzp8xo9I7rpnZlSStXRboObiNyJFdhygCnGmJP2/h7AJuB24PvAn1pynvySfHZl76KzdMVRlsCXZ/LYn5HH/ox8hnctgAgr3+HgMihyPnZhQQmnO8SRGzGayHEz2D9uFmGddUorpZT6Ou02uAE/tbf/ryqwARhjzovII8Bm4GkR+XNTH08ah4OMk1+x8sCLvFK8BSMwqcDw4dnfOOW7WjiM2K6nmVxczIyia1wzgZwKGkRB1Bi6Jk5mwIhpJIaEtvAylVKq/WmXwU1E+gCjgDLgnbrpxpgtInIOiAZuBnY0Vl7q+RR2H3yfAVkXCT6/j/7FR4ihkLkdO/BydC8AjnRy0IUCCuhSfVy+/2T+lH+Esl6j6DRhKh2GTWRox0DXXahSSrVT7TK4ASPs7RFjTHEDefZgBbcRNBLcTl46wu0bFtOl0sHWjLNODTqwrJyoigp6VlQy4VoJp6Mu0DluKCNjwhjZL4w+YcGIzHfNFSmllKrWXoPbAHt7ppE8GXXy1qvqeWWBvx+HAjsyorQMgDy6cCZ4CP8RkEz3pCkMGD6JTvqIUSmlPKK9BreqlTOLGslTaG+71E0QkYeBh+0fSw8/dPgwwEinXFeBc8DHwG9bUNVWJQK45O1K+AhtixraFjUinn32WW0Ly0B3Ft5eg1vVmGfTaK4GGGNeAF4AEJG9xpjRrqpYa6ZtUUPbooa2RQ1tixoisted5bfX0b4F9jakkTxVaQWN5FFKKeWD2mtwS7e3MY3k6Vsnr1JKqVaivQa3/fZ2iIg0NJPwmDp5G/KCa6rUJmhb1NC2qKFtUUPbooZb20KMadbXTq2eiHyJ1QfkQWPMa3XSpmIN4s4Bot09x6RSSinXaq93bgDP2dvfiEh81U4RiQL+Yv+4RAObUkq1Pu32zg1ARP6CtSJACfApNRMndwXWAHe5a+JkpZRS7tOe79wwxjwK3AccBWYBC7GWADXAbViTJzebiHxTRLaJyBURKRSRvSLymIj4fLu7su4i0kdE/iwix0WkWERKROSkiPxNRGLdUX9XcvXvUUSCReQnIrJHRPJF5JqInBaRd0Rkoqvr70rufE+LyK9FxNivp1xRX3dyRVuISAcRmSEifxCRXSKSLSJlInJORN4VkWluvASXccNnpOXlGWPa/Qv4I1ZAq/u6qwVlPm+XUQysB1Zjjew2wCrA39vX7Ym6Y01flmcfm4l1R7wGOGvvKwAmePuaPfV7xJrx5qR9/HngPWAFsBtrrtN/9/Y1e6ot6pQ9BqjAmvTHAE95+3o90RbAzFp/b7Ltst4GDtXa/5/evl5Pvi9c1rbebhi2YJGyAAAJVklEQVRfeAH/jDWNyN1AHFZnkmYHN+DOWm/WhFr7e2DdJRrgCW9ftyfqjjUvp8HqGdWh1v4OwDI77aC3r9tDbdEZSK36g1W7Pez07kCit6/bE21Rp+xA4AjWlD6rfT24ubItgFuAd4HJ9aT9E1bAN8B0b1+3J94XLm1bbzeOL75cENz22sc/UE/a1Fq/PD9vX6s76w4EUfO/z571pPeuld7J29fu7t8jVicmA7zq7WvzdlvUOf439vELgFdaQXDz2OcbeNEub5m3r9sTbeHSvz/ebhxffLUkuAF97GNLgeAG8lQ9kvOpx3GurjvW3Vm5nb9XPem97LRC7M5NvvJyQ1t0xJpf0QCDvH193myLOseNw7o7edP+2aeDm6c/38Bjdlkfefva3d0Wri7P5zs2tEI3upxO7by+wqV1N8aUAxvtH58VkQ5Vafa/f2n/uMzY71wf4urf4yisx46ZxphjIjLB7kDxfyLyrIiMb2mF3cgt72kRCQJeBXKBJ5pfPY/y9Oc7wd5mu6AsV3N1W7i0vPY6cbI7uWw5HS9wR90fBTYA3wPm1posdQwQBvwJ+HET6+kJrm6Lm+ztSRF5BXiwTvozIrIS+FYjH2xvcdd7+ldYM8PfY4xpLTPle+zzLSI9gYfsH1e2pCw3cXVbuLQ8DW6u16LldLzM5XU3xpwSkQnAa8BcrEcPVfYCW+07PF/j6rYIt7dTAH/g98DfgMv2vr9gfZl+FfhOUyvrZi5/X9jviR8Ca4wxb7egbp7mkc+3iAQAbwChwEZjzLrmluVGrm4Ll5anjyVdr0XL6XiZy+tu/xE7DMQDi7DW9orEGkcYBqwUkWdcdT4XcnVbVH3WArAew/7YGJNmjMk3xqzFag8DPOiDY/9c2hb2fK4vYwXyR11Rpgd56vP9N6wJJTKB+918ruZydVu4tDwNbq7XmpfTcWndRaQb1pi2LsAcY8xaY8xlY8wlY8x7wByssSy/EJGExsryAlf/Hmvn+XvdRGPMXuBLrM/ktBsoz5Nc3Ra/BhKBJ40xvvhdUmPc/vkWkT8B38Wa23aGMSanOeV4gLs+Iy4pT4Ob66Xb25hG8vjqcjrp9tZVdZ+HdZe2yxhzqm6iMSYV+ALrbmbajVbSQ9Ltravaonae0w3kqdrf8wbK86R0e+uqtrgda7D2gyKyufYL6z88AI/Y+15sRn3dKd3euuXzLSJ/AH4AXMQKbCebWoYHpdtbV39GXFKefufmek7L6TTQOeBGl9PxNFfXvZ+9vdJInnx7G95IHm9wdVvsq/Xv7lh/vOqKsLeF9aR5kzve035Y45YaEmu/ut1geZ7its+3iPwWeBLre9hZxpijza+mR7i6LVxant65uZgxJhPrD1lHYHHddHs5nT5Yjxx2erZ2jXND3bPs7ajawwBqldcBq4s8NHw34xWubgtjzDmsu1SwvkupW14Y1hJMYHW08RluaIv+xhip74U1NADgx/a+ZNddScu56/MtIkuweg3nYQW2gy6psBu54X3h2rb19kBAX3xxA4O4sWabSAGeqyftLmpG0sfX2h+FNc2QwXen32py3RtqC/uYIvuY/wUCa6UFAn+103KBUG9fuzvbwk5bQM2cksm19gcB/7DT9uJjA9rd0RaNnOcVfHgQt5veF/9lH5MHjPL29Xm5LVz2t9PrjeMLL6z/Me+q9aqapPNE7f11jqn6EL7SQJl/oWbyz3VYE35esfetxrcnTm5S3RtrC6zxXFXz450D1tplZtn7SoDbvH3NnmgLO/131MzCsNUu45y97yy15tPztZer26KBc1Qd47PBzZVtgbUSibFfe+x89b2e9vY1e+p90dTyGqyXtxvGF15YnRnM172a8guy83wT2I4VLIuwesM9hg/OKdmSut/Am3Uk1ji301jBrARIw5o3b7C3r9WTbWHnuR34DOt/6qVYqwT8AYj09rV6ui0aOcang5ur2gJrkPbX/u0BNnv7ej35vnDF3852vVipUkqptkk7lCillGpzNLgppZRqczS4KaWUanM0uCmllGpzNLgppZRqczS4KaWUanM0uCmllGpzNLgp1YaIyAgRMSLykrfropQ3aXBTqm25w96u9motlPIynaFEqTZERI5gLTUUaYwp8XZ9lPIWvXNTqo0QkURgMPChBjbV3mlwU8qD7O/DjP3vh0Rkr4gUiUiOiCwTkUg7LUhEnhWREyJSIiIZIvKr+tbFq+VOe7uq1vmm2efcbJf5XyKSKiLFInJKRP5dRPztvH3tOpyzz3lIRO53V1so5U76WFIpD6oKbMBvgR8CW4ACYALQE/gKmAh8BAyy0wOxVq3uBPzdGPNwA2XvBoYDEcaYAnvfNGAT1uKOlcAQrPUKg4Epdpl/A36PNQv7NWA3EA1Msou+3xjzpgsuXymP0eCmlAfVCm7ngenGmGP2/jCsADQQOAzkA/ONMVfs9GSs9b78gQHGmDN1yu0DZGA9kpxXa/80rOAG8HmdMofXKvM48DHwI2NMpZ3+GNYis2nGmHgXNoNSbqePJZXyjmeqAhuAMSYP6w4KrO/NHq4KQnb6AeADQLDu4uq6w05bVU8agKOeMg/aZfph3cH9pCqw2f4Pa5X0OBHpV7swEXlARI7bjzfTROTRG7lopTxFg5tS3rGhnn2p9vZM7cBXy0l727uetDuwHjuubeB8DZVZdc7PjDFltROMMRVYC8w6nVNE5mMtNPsLoCvwPeA5EbmvgXMr5XEa3JTyjrP17CtsJK12elDtnSISgfX92OfGmItNOF9zz/kEsM4Ys8IYU26M+Qwr2D3RQBlKeZwGN6W8wBjjaCS5sbT63Ib1vVlDjyRvpMymnHMU8EWdfbuBESKif1OUT9A3olKt3+32do2HzhcK5NXZlwcEYH13p5TXaXBTqhUTkS7ADGCvMSbDQ6e9AoTV2RcGVGANJVDK6zS4KdW6zccaB+fJuSS/BMbV2TcW2P81j1uV8hgNbkq1blUTJTf2fZur/QlYICKLRaSDiEwH/tner5RP0OCmVCslIkHAHOCYMSbFU+c1xqzH6v7/K6zZVZYBP9VZTJQv0RlKlGqlRGQRVieSXxtjfu7t+ijlSwK8XQGlVLMVA88Cb3i7Ikr5Gr1zU0op1ebod25KKaXaHA1uSiml2hwNbkoppdocDW5KKaXaHA1uSiml2hwNbkoppdocDW5KKaXanP8PTdIo9gjLVccAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(rk2[:,2]/m0,rk2[:,1],label='rk2',ls='-')\n",
"plt.plot(heun[:,2]/m0,heun[:,1],label='Heun\\'s Method',ls='--')\n",
"plt.plot(heun[:,2]/m0,v(heun[:,2]),label='Tsiolkovsky', 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_dm` that returns the final height of the firework when it reaches $m_{f}=0.05~kg$. The inputs should be \n",
"\n",
"$f_{m}= f_{m}(\\frac{dm}{dt},~parameters)$\n",
"\n",
"where $\\frac{dm}{dt}$ is the variable we are using to find a root and $parameters$ are the known values, `m0=0.25, c=0.18e-3, u=250`. When $f_{m}(\\frac{dm}{dt}) = 0$, we have found the correct root. \n",
"\n",
"Plot the height as a function of time and use a star to denote detonation at the correct height with a `'*'`-marker\n",
"\n",
"Approach the solution in two steps, use the incremental search [`incsearch`](../notebooks/04_Getting_to_the_root.ipynb) with 5-10 sub-intervals _we want to limit the number of times we call the function_. Then, use the modified secant method to find the true root of the function.\n",
"\n",
"a. Use the incremental search to find the two closest mass change rates within the interval $\\frac{dm}{dt}=0.05-0.4~kg/s.$\n",
"\n",
"b. Use the modified secant method to find the root of the function $f_{m}$.\n",
"\n",
"c. Plot your solution for the height as a function of time and indicate the detonation with a `*`-marker."
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [],
"source": [
"def f_m(dmdt,m0=0.25, c=0.18e-3, u=250):\n",
"\n",
" y0 = 0 \n",
" v0 = 0 \n",
" time_2 = (0.25-0.05)/dmdt\n",
" t2= np.linspace(0,time_2,10000)\n",
" dt=t2[1]-t2[0]\n",
" \n",
" N=int(time_2/dt)\n",
" height_desired = 300\n",
"\n",
" \n",
" height2 = np.zeros([N,3])\n",
"\n",
" \n",
" height2[0,0] = y0\n",
" height2[0,1] = v0\n",
" height2[0,2] = m0\n",
"\n",
" for i in range(N-1):\n",
" height2[i+1] = rk2_step(height2[i],lambda state: rocket(state,dmdt=dmdt, u=250),dt)\n",
" height_predicted = height2[:,0]\n",
" error = height_desired - height_predicted[-1]\n",
" \n",
" return error"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [],
"source": [
"def incsearch(func,xmin,xmax,ns=50):\n",
" \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('There are no brackets found\\n')\n",
" print('Please check interval or increase ns\\n')\n",
" else:\n",
" print('These are the number of brackets: {}\\n'.format(nb))\n",
" return xb\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [],
"source": [
"def mod_secant(func,dx,x0,es=0.0001,maxit=50):\n",
" \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": 36,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"These are the number of brackets: 1\n",
"\n",
"The Lower is: 0.05 kg/s\n",
"The Upper is: 0.0889 kg/s\n",
"\n",
"The actual root is: 0.0791 kg/s\n"
]
}
],
"source": [
"\n",
"interval = incsearch(f_m,0.05,0.4,ns=10)\n",
"root = mod_secant(f_m,0.00001,0.1)[0]\n",
"print('The Lower is:', round(interval[1,0],4),'kg/s')\n",
"print('The Upper is:', round(interval[0,0],4),'kg/s')\n",
"print()\n",
"print('The actual root is:', round(root,4),'kg/s')"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"3c)\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAbQAAAE0CAYAAABEhMaeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXwdVf3/8dcnSdN03ze674WW0g0oBSkIWFEW2RQBBUFRQAX1B4LLV8SFzV0QRQVkU9mUTdkECpS1LZTS0kJbuu970jTN9vn9MZPkJs29yU3uzd3ez8fjPubOmTMzn3ub5pNz5swZc3dEREQyXV6qAxAREUkEJTQREckKSmgiIpIVlNBERCQrKKGJiEhWUEITEZGsoIQmGcPMLjAzD1/HJLp+a+JJ9LHb4vjpoMG/UUteKyOOdVdY9mLqPpGkkhKaSBYzs2sb/uIXyVZKaCKSSvcCXaK8ro+oNz5KnYPaMlhJbwWpDkAkU7n7XcBdKQ4jo7l7JVDS2DYzK49YLXX3RutFHOsC4IKEBScZRy00ERHJCkpokvPMrIOZXWFms81si5mVm9lGM/u3mX0qxn5NDtows3Zm9h0ze9vM9pjZdjN72cy+EG5v9jUuM+ttZjeb2YdmVmZm28zsv2Z2bCN1jwnj+lFYNLSRARUvNnXO8FifitjniCbqDjWz6rDuV5pz/ESJNSik4b+VmR1gZr81s+VmttfMVpnZn8zsgIh98szsQjN71cx2mFmJmb1iZic3I5aCcN+nwp+l8vBn6xkzO9fMLKEfXgB1OUqOM7OJwGPA0Aab+gGnAqea2Z3AV9y9Ks5jdwWeAQ6PKO4IHAUcZWbHASubeayDwmMNjChuD3wSmGVmF7j73fHEF4dngM1AX+A84LUYdc8BDNgHPJikeFrFzCYRfKY+EcVDgIsJvssjga3APwl+BiIdCTxmZhe5+x1Rjj+E4GfqkAabegMnhK9zzewsd9/T2s8jddRCk5wV/uJ5gSCZrQa+CowCegITgJuAKuBLwI9bcIq/UpfM7gKmAL3C5d+A84EvNPNYjwMVwBeBwQS/jE8D1hAkkFvNrFdE/ZepP7BiNfsPqDixOScOr3P9I1z9rJm1i1H93HD5pLvvbM7xU+ARYBfwWWAAMAi4Aqgk+Fm4HriR4Pu5FhhH8O/2ceD98Bi/bfB9A7V/xDxPkMy2Ad8BDgR6AGOAa4C94bH/lIwPl9PcXS+9MuJFcMHfw9eJQOcmXl+NqH9MI8d7LNy2FugT5ZwXhXX2AQdEi6eR/Y6IOPfvoxz7tog6K5v4vGuBvo3UmRJR52uNbL822vHj/O6nRZzn5Ch1JkXU+UwC/r2vjTjesGbUvyus+2IT3+VqoHcjdX4abq8EqoEzG6kzhuCPHAe+2sj234XbdgNjo8R5fEQs01LxfylbX2qhSab6D1DcxOuP0XY2sxHASeHqFe6+JUrVO4DlQCFwVhzxfTFclgLfj1LnaoK/1pvjOnff3LDQ3ecD74arh8YRX1zcfS6wJFw9N0q188LldoJ/n3R1nbtvbaS8phWaD8xx94caVnD3D4C3w9XIrmTMrBPBH0A151ja2Mnd/TmCVhxE/y6lBZTQJFcdR9BVVw28amadG3sBnYAF4T7T4jj+jHD5orvvbqyCu+8CZjfzeP+Nsa3mF2f/Zh6rpe4Nl6eYWZfIDWaWB5wdrj7o7uWkr6ejlC9vRp3Ieg2/7xkE10gBZkf7mQp/rmr+CInnZ0qaoIQmmepYd7dYL4JrX9GMDZd5wDpit/ROD+v2ofmGhctG/0qPsKSJ7TXWx9hWGi47xqiTCPcRdJN1AM5osO1Y6gas3Et6a/S7dPfI1vKGGPvX1OvQoHxsxPs3if0zdUVYL56fKWmCEprkqm4t2KcojrqdwmVTo9hi3ixcw5s3wjKpQ8HdfSXwSrh6XoPNNV1nHwFzkhlHazXzu2zJ953snylpgobtS66qSSRb3T0ZfyXvAbpSl9ii6ZyEcyfTvcDHgGPN7AB3X29mRdS12O7zcORDDor846SLNzGziSSeWmiSq1aEy95mNigJx18VLsc0UW9sE9vTzQMEIz7zgM+HZScTJG9I/+7GZFoR8X5SyqLIYUpokquejXgf61pbS9V0ux3bcABFjfCepWOScO5IFeEyPxEH8+DesifD1fMaLOdGG9mXI2YTJHtIzs+UNEEJTXKSuy+hbmj5NWY2I1Z9M+trZj3iOMU94bIj8JModa5n/4EFibYtXPYxs0RdYqj5bJPM7GMEs5VEluekcDTrX8PVC8ys4cCZesysq5kNSH5kuUMJTXLZ1wimdOoAvGBmvzSz6eGcib3M7EAzO8fM/k7QhTiyuQd291cJZqQAuNzM/mpmh5hZTzObZGZ3AJdSv5sqGeaFy/bAdeEchu3CuQZb2mr7D8G9ZgB3E9yjFzmbSC77HsHI1jzgwfDffaaZ9TOzHmY22szOMLO/EMzycmRKo80ySmiSs9x9DTCTYDqjQuDbBPMUbiGYy28xwVD1swkSQkXjR4rqQuCtiPfvELSY3ibokrqbumtOlS39HLG4+1vAq+HqNQS3KJQTfJb/tfCY5QTX0qDu9oRnG7vxO9eE9xYeS9DlbAT/7i8CGwn+CPgAeIjgBuyuBP8WkiBKaJLTwq7HiQTzKj5OcI9SOcG1kDXAU8A3gcHuviDacaIcexfBRMRXEtxIuxfYSZBgLnT386kb5Vjc6g8T3aeAmwkSdGkTdZurYfdiLg8GqcfdNxCMBD2NYILm1UAZwc/VBoL5Q68GRrv7Y6mKMxtZ7o6wFUk9M3sUOAV4wt2bfCyJiESnFppIioRTINU8y2xerLoi0jQlNJEkCeftK4xR5SaCx7hA3TUpEWkhJTSR5JkALDGz75rZ5HCEYz8zO97MHgMuCev93d0XpzBOkayga2giSWJm04n9dGcIbsY9NRxAIiKtoISWIr179/Zhw4alOgxJoqqqKrZv387u3bvZu3cvlZWVVFdXk5+fT8eOHenZsyc9e/bELKlzCotklXnz5kWdf1WTE6fIsGHDmDt3bqrDEBHJKGa2Kto2XUMTEZGsoIQmIiJZQQlNRESyghKaiIhkBSU0ERHJCkpoIiLSdoo3wp0nQvGmhB9aCU1ERNrE26t3sOLhH+GrXsNn35Dw4+s+NBERSb6f9mVy5b669bl3BK+C9vCDxDxKTy00ERFJvsvf5fl2M9nrwXzdVflFcPBZcPnChJ1CCU1ERJKupLA36/cW0J4KyrwdedXl0L4rdOmXsHOoy1FERJLu3bU76W27ubfqOF7rcQq3jV0AJYkdGKKEJiIiSff26p3cXPEtAD4/fAicdH7Cz6EuRxERSbq3V++ofT95SPeknEMJTUREksrdmb96Z+36lCE9knIeJTQREUmq1dtL2b6nHICuRQWM6N0pKedRQhMRkaR6O6J1NmlID/LykvNQWyU0ERFJqrmrtte+nzw4OdfPQAlNRESS7PUVdQnt8OE9k3aetEloZjbWzC43s3vNbImZVZuZm9mZzdj3HDN72cx2mVmJmc01s8vMLObna+v9RERyzdaSfSzbXAJAu3xjcpIGhEB63Yd2CXB5vDuZ2a3ApUAZ8D+gAjgOuAU4zszOcveqVO8nIpKL3vyornU2aXB3OhTmJ+1c6dSieA+4GfgcMAqY3dQOZnYGQXLZCEx095Pc/TRgNPA+cBrw9VTvJyKSq15fsa32/eHDeyX1XGmT0Nz9L+5+lbs/4O7Lm7nbNeHyu+7+YcSxNhG0+ACubqQrsK33ExHJSW9EXj8bkbzrZ5BGCS1eZjYImAqUAw823O7us4F1QH9geqr2ExHJVdv3lLN0UzEABXnG1KHJu34GGZzQgMnhcpG7741S560GdVOxn4hITnojortx4qBudCxM7rCNTE5ow8Plqhh1Vjeom4r9RERy0ksfbq19f8TI5F4/g8xOaJ3D5Z4YdUrCZZcU7lfLzC4Oh/jP3bJlS4zDiIhkNnfnpQ/qfs8dPbpP0s+ZyQmtZu4UT/P9arn77e4+zd2n9emT/H9cEZFU+WjrHtbtDK7OdCrMT+r9ZzUyOaEVh8vOMerUbCuOKGvr/UREck5k6+yIkb0pLEh+usnkhLYyXA6NUWdwg7qp2E9EJOe8HHH9bOaY3m1yzkxOaG+Hy/Fm1iFKnUMb1E3FfiIiOaW8sprXIkY4fqwNrp9BBic0d18DzAcKgbMabjezmcAgglk9XkvVfiIiueatldspLQ9mABzcswPDkvT8s4YyNqGFrg+XN5rZqJpCM+sL/CFcvcHdq1O8n4hIznh28aba98eN69dm502byYnNbAp1SQHgoHD5czP7fzWF7j494v1DZnYbwbRTC83sOeomC+4K/Jtg0uB62no/EZFc4e78b0lEQjuwb5udO20SGkFCOLyR8tGxdnL3S83sFeAyYCaQDywB7gBui9Zaauv9RERywQebSlizPRiu37l9QdInJI6UNgnN3V+k7l6vePe9H7g/3fcTEcl2z71f1zqbObZPmwzXr5Hp19BERCSNRF4/O+HAtrt+BkpoIiKSIBt3lfHOmp0A5OcZx4xt2xmRlNBERCQhnly4ofb99BE96d6xsE3Pr4QmIiIJ8eS762vfnzTxgDY/vxKaiIi02rqde5m/OuhuLMgzPjm+f5vHoIQmIiKtFtk6O3JUb3p0atvuRlBCExGRBHjy3brrZ5+eOCAlMSihiYhIq6zeVsqCtbsAaJdvzDqo7bsbQQlNRERa6fGI7saPje5Dt47tUhKHEpqIiLSYu/PQvLW166cc0vajG2sooYmISIvNW7WDj7buAaBL+wJmpWB0Yw0lNBERabEH59a1zk46ZAAdCvNTFosSmoiItEhpeSVPRFw/O3Pq4BRGo4QmIiIt9NR7G9kTPpl6RJ9OTBnSPaXxKKGJiEiLRHY3njl1EGYtegJYwiihiYhI3JZtLuG1FdsAyDM4ffKgFEcU4wGfZnZHgs7h7n5Rgo4lIiJp4N7XV9W+P/7AfvTvVpTCaAKxnlh9AeC08CnSERxQQhMRyRJ79lXycMS9Z188YljqgokQK6EBzAH+2orjfxmY0Yr9RUQkzfzr7XUU76sEgsEgR47qleKIAk0ltGXu/reWHtzMjkEJTUQka7g797xW1934helDUz4YpEasQSG7gdJWHn9veBwREckCr6/YztJNxQB0LMznjKmpHwxSI2oLzd1bfUOBu18KXNra44iISHr488srat+fNnkgXYtSMxFxYzRsX0REmmXpxmKeX7IZADO46KjhKY6oPiU0ERFplttfqmudfeKgfozo0zmF0exPCU1ERJq0YddeHn1nXe36V2eOTGE0jYsroZnZKDP7s5ktM7NSM6uK8qpMVsAiItL27njlIyqrHYDDhvVkypAeKY5of00N269lZtOA54FONH2zdXqM4RQRkVbbWrKPe19fXbt+8dEjUhhNdPG00G4COgMPAFOALu6eF+2VlGhFRKTN/fmlFeytCGbVH9e/Cx8f1zfFETWu2S004HDgfXf/fLKCERGR9LK1ZB93R9xIfcXxY8jLS89OuHhaUnuBBckKRERE0s+fZi+vbZ0dNKArs8b3S3FE0cWT0N4E0rPjVEREEm5zcRn3RMyqf/nxo9NmmqvGxJPQfgZMNrPTkxWMiIikj18/+yFlFdUAjD+gK584KH1bZxDHNTR3n2NmZwN/NrPTgKeBtUB1lPovJSZEERFpax9sKuafb9WNbLzqk+PSunUG8Q0KASgkmLD4nPAVjbfg2CIikiZu+O8SwtvO+Njo3swc0ye1ATVDPPehnQHcR9BNuQ1YCZQkJywREUmVV5dtrTdn4zUnHpjiiJonnlbU9whumL4UuN3dG+1qFBGRzFVd7fzsP+/Xrp8xZRAHHdA1hRE1XzwJbRwwx93/mKxgREQktR6ct4ZF64PHWBa1y+M7nxiT4oiaL55RjrsIBoGIiEgW2rGnnBv+u6R2/ctHjWBAtw4pjCg+8SS0Z4BDLd2HuYiISIvc9PQSdpRWADCwewcuO3ZUiiOKTzwJ7ftAF+AXZqYRjCIiWWT+6h384601tevXnjKeDoX5KYwofvEkpouA/wBXAKeZ2fNEvw/N3f0nCYhPRESSrKra+eG/38PDYfrHH9iXE9L8JurGxJPQriW4v8yAYcCFjdSp2e6AEpqISAa445WPageCtC/I40cnj09xRC0TT0K7jiBRiYhIlli+pYRfPLO0dv2bx41mcM+OKYyo5eKZ+uraJMYhIiJtrKraueqhd9lXWTdfY7o+vLM59CBOEZEcdeecj5i3agcABXnGL846hHb5mZsWMjdyERFpseVbSrj56bquxm98fDQHDsiMGUGiiZrQzOx0M5vSmoOb2RQ9bkZEJL3sq6zim39/u7ar8aABXbn02JEpjqr1YrXQHgK+3srjfwN4sJXHEBGRBLr5qaW1oxoLC/L45Wczu6uxRuZ/AhERabYXl27mL698VLv+vRPHZXxXY42mRjl+MryBuqXGtWJfERFJoM3FZfy/BxfUrh83ri/nzxiWuoASrKmE1j98tYbuXRMRSbGKqmq+ft/bbC0pB6Bvl/bcdObEtH8KdTxiJbRj2ywKERFJquv/s4Q3V24HIM/g15+bRK/O7VMcVWJFTWjuPrstAxERkeR49J113DGn7rrZdz4xliNH9U5hRMmhQSEiIllsycbdXP3wwtr1WeP7cekxmT9EvzFKaCIiWWr7nnK+es889lZUATCiTyd+cdYhWXXdLJISmohIFiqrqOLiu+eyalspAB0L8/nTeVPpUtQuxZEljxKaiEiWcQ8mHZ4bztNo4SCQ0f26pDiy5FJCExHJMr9+9gMeW7C+dv37nzqQWeNbewdW+lNCExHJIg/OXcPvnl9Wu37e9CFcdNTwFEbUdpTQRESyxDOLNnL1I3UjGmeO6cO1J4/P2kEgDTU7oZnZF81sRjPqTTezL7YuLBERicery7fy9b+/TVV1MDnTgQO6css5kynIgkmHmyueT3oX8OVm1LsIuLNF0YiISNzeXbuTr/xtLuXh42CG9erI3RceltUjGhuTjNSdG21bEZE0sHRjMRfc+RZ7yoN7zfp1bc89Fx1Ony7ZNa1VcyQjoQ0CSpJwXBERibB0YzHn/Pl1tu8JJhzu3rEd9150OIN7dkxxZKkRc7b9Rq6FjYpxfawAOBA4DngrAbGJiEgUSzbu5pw/v1GbzLq0L+CuLx2W9feaxdLU42Puov7jX44MX9EYUA38onVhiYhINO9v2M25f6mfzO6+6DAmDe6e4shSq6mEdjd1Ce18YDkwJ0rdcmAd8Ki7L4hSR0REWmHx+t2c99f9k9nkIT1SHFnqxUxo7n5BzXszOx94xd0vTHZQIiKyv7dWbufCu96iuKwSUDJrqKkWWqThaLCHiEhKvLBkM5fcN4+yimBofpeiAu656PCc72aM1OyE5u6rkhmIiIg07tF31vGdBxZQGd403btze+6+8DAOOqBriiNLL/G00AAwsyJgGnAAUBStnrvf3Yq4REQEuGvOR/z4icV4OJphUI8O3HvR4Qzr3Sm1gaWhuBKamX0L+D+gOX8WKKGJiLRQVbXzsyff5445H9WWjenXmXsuOpx+XaO2JXJasxOamV0I/DJcfR9YAuxORlAiIrmstLySK/7xDs8s3lRbNnlId+684FC6dyxMYWTpLZ4W2jcJhvB/wd3vT1I8IiI5bXNxGV/+21zeXburtuzECf351Wcn0aEwP4WRpb94pr4aA7yabsnMzO4yM4/xWhJj33PM7GUz22VmJWY218wuM7OY30tL9xMRieWDTcWcduur9ZLZxUeP4NZzpiiZNUM8LbRSYHWyAkmAOcCyRso3NFbZzG4FLgXKgP8BFQTTdt0CHGdmZ7l7VaL2ExGJ5ZlFG/nWP9+pnWQ4z+DHp07gC9OHpjiyzBFPQnsVmJCsQBLgL+5+V3MqmtkZBElpI3C0u38YlvcDXgBOA74O/DYR+4mIRFNd7fz++WX8+rkPass6FuZz6zlTOHZc3xRGlnni6SL7MTAunDEk010TLr9bk5QA3H0TcEm4enUjXYgt3U9EZD979lVy6X3z6yWzQT068PAlM5TMWiBqC83Mjm6k+FfAHWb2KeBJgi7I6sb2d/eXEhJhgpnZIGAqwdyTDzbc7u6zzWwdMBCYTtAybfF+IiKNWb2tlIvvmcuSjcW1ZTNG9uKWc6bQs5NGMrZErC7HF6k/034NA84MX9F4E8dOhmPNbCLQGdgEvAI86+4NE+7kcLnI3fdGOdZbBIlpMnWJqaX7iYjUM2fZVi67fz47Sytqy7505DC+/6kDKchXB09LxUo6L9F4QktXjT2nbbGZne3uCyPKhofLWFN51Qx+GR5R1tL9REQAcHdum72cXzy9lHAWKwrz8/jpaRP47LTBqQ0uC0RNaO5+TBvG0RrvAPMIRhyuIpjFZArwM+AQ4Dkzm+Lu68L6ncPlnhjHrJmEOfJJeS3dr5aZXQxcDDBkyJAYhxGRbLO7rILvPLCAZyNulu7bpT1//MJUpmi2/IRo627BhHP33zQo2gM8aWbPArMJrmddQzD6EIIuU4i/9dnS/Wq5++3A7QDTpk3LpNaviLTCko27+do981i5rbS27NBhPbj1nCn01TRWCZPxCS0ady83s+uBR4FPRWyquQLbef+9atVsK44oa+l+IpLDHn1nHVc/vJC9FXW3p1501HCuPnEc7XS9LKHimcuxsVGPjSkHtrp7Yzc5t7WaWUIGRpStDJex7las6cxeGVHW0v1EJAeVV1bz8/+8z12vrqwt61iYz41nTOTkQw5IXWBZLJ4W2ovE0d1mZruBvwE/dPdUtVh6hcvIB5O+HS7Hm1mHKCMWD21QtzX7iUiO2birjEvvm8f81Ttry0b06cSfzpvK6H6NXmKXBIinvfsS8BrBtSQDdgLvEgzK2EHdNaY3gBUE3W/fAF42s46JCjhOnw2Xb9UUuPsaYD5QCJzVcAczmwkMIpgN5LXW7iciueW15ds46fcv10tmnxzfn0cvO1LJLMniSWifDJeLgU+5ey93n+zuU929N3AisIigFXcwMJrgXqyDCWbqTzgzm2RmJ5lZfoPyAjP7dsR5f91g1+vD5Y1mNipiv77AH8LVGxq5h62l+4lIlnN3bn9pOef99Q22lpQDwXyM15w4jtvOm0KXonYpjjD7xdPl+APCROXumxtudPenzewd4APg/9z9e2Z2DrAUOAO4IREBNzAM+Bew3cw+ANYSDJk/mOCJ2tUE01Q93SDWh8zsNoLpqhaa2XPUTTLcFfg3wWTDDT9ji/YTkexWXFbBVQ+9y3/f21hb1rtzIb/7/GRmjOydwshySzwJ7XPAC40lsxruvsnMXiDo6vueu68xs/kECSYZFhBMBHwYwWCNyQQtxLXAncCt7j4vSqyXmtkrwGXATCCfYBDJHcBt0VpZLd1PRLLTovW7uOy++fWG5E8Z0p0/nDuV/t00JL8txZPQBhFcQ2rKPuqPKlwDTIsnqOZy94+AK1qx//1A3M93a+l+IpI93J3731zNjx9fTHll3d+x5x8xlO9/+iAKCzQkv63Fk9C2AkfHGOGHmXUAjga2RRT3IBhAIiKSFUr2VfK9Rxby2IL1tWWdCvP5+ekHc+qkgTH2lGSK50+Ix4F+wANmtt+kY2HZP4G+wGMRm8YRjHoUEcl472/YzSm/f6VeMhvXvwuPfeMoJbMUi6eF9iOCkYyfBpaZ2WsEcyc6wfWrGUC7sOxHAGY2FRgC3J3AmEVE2py788+31vCjxxaxL6KL8exDB3PtKeMpapcfY29pC81OaO6+xcxmALcBJxN0LdarAjwBXOLuW8J95plZO3evQkQkQ5WWV/KDf73HI2+vqy3r0C6fn58+gdMmD0phZBIprrkc3X0D8BkzG0KQ0Gra1+uBl919ZSP7KJmJSMb6YFMxl943n2Wb6yYcGtOvM384dwqj+upG6XTSosmJ3X01cG+CYxERSRvuzt/fXMN1TyyirKKui/GsqYO47tQJdChUF2O6ydrZ9kVEWmpnaTlXP7yQpxbV3Shd1C6Pn5w6gbP0IM60FTWhhd2KAOvcvSpivVnCVpyISEZ586PtXPGPt1m/q6y2bEy/zvz+81MY219djOksVgttJcHUUQcRTGe1kubPtu9NHFtEJK1UVlXzu+eXccvzH1Id8ZvuvOlD+MGnD9IoxgwQK+msJkhMFQ3WRUSyytodpVzxj3eYu2pHbVn3ju248YyJzBrfP4WRSTyiJjR3HxZrXUQkGzzx7nqueWQhxWWVtWXTR/Tk15+bxIBuHVIYmcRL3YIikpP27KvkuscX88+5a2rL8vOMbx0/mkuOGUV+nsXYW9KREpqI5Jy5K7fz7QcWsHp73Qz5g3p04LdnT2bq0B4pjExaI+6EFj7Y8qvAEUAf4FF3vyrcNh2YCDzg7pqQWETSSnllNb957gP+OHt5vYEfpxxyAD89bQJd9RDOjBZXQjOzi4BbgcKwyIHIp9f1IZgaq4LgeWQiImlh6cZivvXPd1i8YXdtWZeiAq47dTyfmTQQM3UxZrpmJzQzOxL4E1ACfB94CXijQbWngN3AKSihiUgaqK52/vrKR9z89FLKq+pm/Jgxshe/OOsQDuiugR/ZIp4W2lUELbIT3f01YL+/aNy9wsyWAgcmLEIRkRZau6OU7zywgDc+2l5b1r4gj+9+chwXzBhGngZ+ZJV4EtoRwJs1ySyGNSihiUgKuTsPzlvLdY8vpmRf3XD8CQO78uvPTmJ0P834kY3iSWjdgLXNqFcY53FFRBJm/c69XPPIQmZ/sKW2LM/gsmNH8Y2Pj6awIJ7nGksmiSfxbAaGN6PeWGBdk7VERBKoZnb8n//n/XqtsmG9OvKrz01iyhANx8928SS0OcCZZjbN3ec2VsHMTgDGAH9JRHAiIs2xZnspVz/yLnOWbastM4MvzRjOlbPG6lEvOSKehPZr4CzgETP7MvBc5EYzOxq4A6gEfp+wCEVEoqiudu59YxU3/HcJpeV1zxIe0bsTN505kWnDeqYwOmlrzU5o7v6GmV0F3Az8l2B4vhM8wfrTBPejGfBtd1+YjGBFRGqs3LqH7z78bn+FeZwAABVOSURBVL0RjHkGX/nYCL51whjNjp+D4hq84e6/NLNFwI+BaQQJrHu4eSHwQ3d/LLEhiojUKa+s5vaXlvO755dRXll3X9novp256cyJTNa1spwV92hEd38KeMrMehEMEskH1rj7+kQHJyISae7K7XzvXwv5YFNJbVl+nvG1mSP45nGjaV+gVlkua/HwenffBmxrsqKISCvt2lvBjU8t4f43VtcrnzCwKzecPpEJA7ulKDJJJ7pfTETSlrvz5MIN/PjxxWwp3ldb3rEwn2+fMIYLZgyjIF/3lUkgakIzsy+25sDufndr9heR3LZmeyn/9+h7vLB0S73y48b15brPTGCg5mCUBmK10O4iGMXYUkpoIhK3sooq/jR7BX94cRn7IgZ99O3SnmtPGc+JE/prZnxpVKyE9hLRE9pMYBOwJOERiUjOem7xJn78xCLWbN9bW2YG5x4+hKs+OU7PK5OYoiY0dz8m2jYzqwb+6+4XJiMoEcktK7fu4bonFvP8ks31yicM7Mp1p07QtFXSLBoUIiIps7e8ij+8uIw/zV5R71ll3Tq048pZY/n8YUPI1yNepJmU0ESkzbk7/1m4kZ//533W7azfvXj2oYO5ctY4enYqTGGEkomU0ESkTb2zZic/eWIx81btqFd+yKBuXHfqBA4Z3D3KniKxKaGJSJtYt3MvNz+1hH+/U39SoR4d2/HdT47js9MG6wnS0ipKaCKSVHv2VfLH2cu5/aUV9Ybht8s3zj9iGN/4+Gi6ddToRWk9JTQRSYqqaueheWv4xTMf1JvlA2DW+H5cfeKBDO/dKUXRSTaKNVPI0U3s2z9WHXd/qcVRiUjGcneeXbyJm59eyoebS+ptmzCwKz/49EFMH9ErRdFJNovVQnuR6DdWOzArfEXbrtafSI55fcU2bnxqCW+v3lmvvF/X9lw5axynTx6o62SSNLGSzmpaN/WViOSIRet3cdNTS5n9Qf15FzsV5vOVo0dw8dEj6Fiov3EluWLNFDKsDeMQkQy0cusefvXsBzy2oP7IxcL8PM6bPpTLjh1Jr87tUxSd5Br9ySQicVu9rZRbX1jGw/PXUlld15GTZ3D6lEFccfxoBvXomMIIJRcpoYlIs63eVsotL3zIw/PXUVVd/4rECQf148pZYxnTr0uKopNcp4QmIk2Klcimj+jJlbPGMXWoJhCW1FJCE5GoVm3bwy3PL+ORtxtPZJcfN4YjRmoIvqQHJTQR2c/7G3bzx9nLeeLdDfslsiNG9OLy40frXjJJO0poIgIEN0S/+dF2bpu9nBeXbtlvuxKZpDslNJEcV13tPPf+Jm6bvXy/G6IBjhzVi29+fDSHK5FJmlNCE8lR+yqreOyd9fzppRUsazBFlRmcOKE/X5s5komD9DgXyQxKaCI5ZnNxGfe9vpr73ljF1pLyetsK8/M4Y+pAvvKxEYzo0zlFEYq0jBKaSI54b90u7pjzEU8s2EB5VXW9bZ3bF3Du9CFceORw+nUtSlGEIq2jhCaSxSqrqnl28SbunLOSN1du32/7gG5FfPGIYZxz+BC6ddAzySSzKaGJZKHNu8v451tr+Mdba1i3c+9+26cM6c6FRw1n1vj+tMvPS0GEIomnhCaSJaqrnTnLt3L/G6t5dvGmenMsAhTkGZ+eOIAvHTmcSYM10EOyjxKaSIbbVrKPB+et5e9vrmbVttL9tvfo2I5zDh/CF6YPo383XR+T7KWEJpKBqqqdV5Zt5cG5a3hm0ab9BnkAHDasJ+dOH8Ks8f0papefgihF2pYSmkgGWba5hIfnr+WR+WvZtHvfftu7FhVw+pRBnHv4EEZr1nvJMUpoImluV2kFj7+7nofmreWdNfvP5AEwaXB3zj18CCdNPIAOhWqNSW5SQhNJQ2UVVby4dAuPL1jPs+9vorxy/y7F3p0LOXXSQM6cOogDB3RNQZQi6UUJTSRNVFRVM2fZVh5bsJ5nFm2iZF/lfnXa5RvHjevHmVMHMXNsHw25F4mghCaSQtXVzpsrt/PYgvX8d+EGdpRWNFrv4IHdOGPKQE6ZNJCenQrbOEqRzKCEJtLGyiureXX5Vp5etIlnF29ia8n+gzsAhvXqyMmHHMDJhxzAGA3wEGmSEppIG9izr5IXl27h6UUbeWHJZoob6U6EYCqqkyYO4JRDBjJhYFfMrI0jFclcSmgiSbJh115eXLqF5xZv4uVlWxsd2AHB4I4TJwzglEkHMHVID/LylMREWkIJTSRBKquqmb96Jy8s3cwLSzazZGNx1LqDe3Zg1kH9mTWhP1OG9CBfSUyk1ZTQRFph0+4yXvpgCy8u3cJLH26huKzxrkSAcf27MGt8f2aN78+BA7qoO1EkwZTQROKws7Sc11dsY86ybby6fCvLt+yJWrcwP4/DhvfkmLF9OOGgfgzt1akNIxXJPUpoIjGUllfy1sodvLpsK3OWb2XR+t24R68/oFsRx47ry7Fj+zJjZC86tdd/MZG2ov9tIhG2FO9j3qrtvLVyB3NX7WDRul37PYYlUmFBHlOH9GDm2D4cO7YvY/p1VleiSIoooUnOqq52VmwtYe7KHby1cgfzVm1nZSOPX4mUZ3DwoO4cObIXR47qzdShPTSTvUiaUEKTnODurN9VxrtrdrJg7S4WrtvJu2t3xRzEUWNMv87MGNmbI0f15vARPela1K4NIhaReCmhSdapSV5LNuzm3bW7eHftThau28XWkvIm9y3Mz+OQwd2YOrQnhw7rwdShPejeUVNNiWQCJTTJaHv2VbJ0UzFLNhSzZONulmwo5v2Nu5vV8oLgac5Th/Zg2rAggU0Y2I32BepCFMlESmgtZGbnAJcAE4F8YAlwJ3Cbuzc+JYS0iLuzfU85K7buYfnmElZs3cOKLSV8uLmEVU1c84rUuX0BBw/sxsRB3Zg4qDsTB3VjUI8OGsQhkiWU0FrAzG4FLgXKgP8BFcBxwC3AcWZ2lrtXpTDEjOPu7N5byZodpazZXsrKbaUs31LCii1BAtsZZRb6aLoWFTBuQFcOGtCVQwYHCWx4r06aVkokiymhxcnMziBIZhuBo939w7C8H/ACcBrwdeC3KQsyDbk7O0or2LS7jHU79oaJay9rd5SyZsde1m4vjTphbyz5ecbIPp0Y178r4wZ0YVz/Lozr35UB3YrU8hLJMUpo8bsmXH63JpkBuPsmM7sEeBG42sx+nwtdj+WV1ewoLWf7nuC1ubiMTbv3sWl3WfgK3m/evY/yqpZ/HR0L8xneuxMj+3RmRJ9OjOjTmZF9OjGqb2dd8xIRQAktLmY2CJgKlAMPNtzu7rPNbB0wEJgOvNq2EcbP3Sktr6JkXyXFZZUUl1VQsq+SkrJKisOykrB8R2kFO0rL2bannB3hqyWtqmg6tMtnUI8ODO7ZkSE9OzKiT10C699VLS4RiU0JLT6Tw+Uid98bpc5bBAltMglOaMVlFdz89FKq3amqDm4MrnKn2j18H5QF253KamdfZRVlFdXsq6xiX0U1ZTXLiir2VVazL8ojTZKhS1EB/boWMaBbEYN6dGRwzw7BMkxivToVKmmJSIspocVneLhcFaPO6gZ1E6asopq7X4t16raXn2f06NiOHh0L6dGpkD5d2tOvSxH9u7WnX9ci+nYpon+3Ivp2aa95DUUkqfQbJj6dw2X0KdahJFx2abjBzC4GLgYYMmRI3CdP1jOz2hfk0aWoHV2KCuhSVEDn9sGrYVlN0urZKUhgvTq1p0tRgUYOikhaUEKLT81v7hjzrUfn7rcDtwNMmzYt7mN0LMzn2pMPIi/PyDMjP8/IN8MsSHb5YXmwDfLz8ihql0f7gvz9lu0L8ihql09hQZ4eLikiWUEJLT41jyDuHKNOzbbojytuoaJ2+VxwZMJ7MkVEskJeqgPIMCvD5dAYdQY3qCsiIm1ACS0+b4fL8WbWIUqdQxvUFRGRNqCEFgd3XwPMBwqBsxpuN7OZwCCCWURea9voRERymxJa/K4Plzea2aiaQjPrC/whXL0hF2YJERFJJxoUEid3f8jMbiOYaX+hmT1H3eTEXYF/E0xSLCIibUgJrQXc/VIzewW4DJhJ3eNj7kCPjxERSQkltBZy9/uB+1Mdh4iIBMy9RfcISyuZ2RZiT6EVS29gawLDyQX6zuKj7ys++r7i05rva6i792lsgxJaBjKzue4+LdVxZBJ9Z/HR9xUffV/xSdb3pVGOIiKSFZTQREQkKyihZabbUx1ABtJ3Fh99X/HR9xWfpHxfuoYmIiJZQS00ERHJCkpoIiKSFZTQMoiZnWNmL5vZLjMrMbO5ZnaZmenfMYKZjTWzy83sXjNbYmbVZuZmdmaqY0s3ZtbOzI4zs1+a2etmtsHMys1snZk9ZGbHpDrGdGNm3zCzB8zsfTPbZmYVZrbFzJ4zs/PMTE/MbYKZ/Tz8P+lm9v8SdlxdQ8sMZnYrcClQBvyPuvkjuwD/As5y96rURZg+zOw3wOWNbDrL3R9q63jSmZkdDzwbrm4E5gF7gIOACWH5T9z9/1IQXloys7VAX+A9YB3B9zUUOJzgqfaPAqdrCrzGmdmhBE8jySP4vq50918k4tj6yz4DmNkZBMlsIzDR3U9y99OA0cD7wGnA11MYYrp5D7gZ+BwwCpid2nDSWjXwMHC0uw8If7Y+5+4HA2cDVcAPzezYlEaZXs4Gerj7FHc/2d3PdvcjgIOBTcCpwPkpjTBNmVl74C6C7+nRRB9fCS0zXBMuv+vuH9YUuvsmgln/Aa5W12PA3f/i7le5+wPuvjzV8aQzd3/e3c9095cb2fZPgl8+AOe1aWBpzN1fcfc9jZQvAm4NV09o26gyxnUErf+vAbsSfXD9AkxzZjYImAqUAw823O7uswm6PfoD09s2OskBNU9eH5TSKDJHZbgsS2kUacjMDge+A9zv7o8n4xxKaOlvcrhc5O57o9R5q0FdkUQZHS43pDSKDGBmwwlaHgBJ+YWdqcysCPgbsJ3Gr28nhB4fk/6Gh8tYM/OvblBXpNXMrD9wQbj6cApDSUtm9iWC5yG2I2jBziBoJFzv7v9KZWxp6GfAWOBsd0/aUwmU0NJf53C5X599hJJw2SXJsUiOMLMC4F6gG/C/ZHURZbgjqT/4oxL4IfCr1ISTnsxsBnAF8O/wumzSqMsx/dXc06L7K6Qt/ZHgtpA1aEBIo9z9y+5uQEdgPPAb4FrgdTM7IJWxpQsz6wDcCewmGKmdVEpo6a84XHaOUadmW3GMOiLNYma/BS4iuE3kOHffmOKQ0pq773X3xe5+JcGI5EOAW1IcVrr4OTAG+La7J/06rLoc09/KcDk0Rp3BDeqKtIiZ/RL4JrCFIJl92MQuUt+dwC+Ak82snbtXpDqgFDuN4F7H882s4b1548LlJWZ2ErDM3b/cmpMpoaW/mmHT482sQ5SRjoc2qCsSNzO7Cfg2sA04wd0XpzikTLST4FpaAdCT4AbiXJdHMHgmmhHhq3siTiRpzN3XAPOBQuCshtvNbCbBCKuNBNPJiMTNzG4ArgR2ECSzBSkOKVMdTZDMdgJJG82XKdx9mLtbYy+CYfwQTH1l7j6ptedTQssM14fLG81sVE2hmfUF/hCu3qC546QlzOwnwHcJfgmf4O5q6UdhZh8zs3PDKZwabjsS+Gu4+lfNrdr2NDlxhjCzPxBMc1UGPEfd5MRdgX8DZ+o/UMDMplCX6CGYaqcL8CHBjZ0AuHvOz6xiZqdQN6feXGBRlKpL3P2GtokqfZnZBQTXyXYS9JxsJPjZGknwcwbwJMFE2NEmQhDAzO4iuO0hYZMT6xpahnD3S83sFeAygv7ofGAJcAdwm1pn9XQlmPm8odGNlOW6nhHvp4WvxswGcj6hEXwPPwE+RjB6bwbBrTUbCW4+v9fd/5268HKbWmgiIpIVdA1NRESyghKaiIhkBSU0ERHJCkpoIiKSFZTQREQkKyihiYhIVlBCExGRrKCEJpJEZuYteN0V7ntMuP5iaj9Fy5jZBY18tmg3bjf3mDsb+65EQDOFiCTb3xop6w/MIngK+UONbH8lqRG1veXUfabWTth7P8EDNUcRPDFapJYSmkgSufsFDcvM7BiChLa1se0R3gQOBEqTEVsbeqWJz9ls7n4p1M6pqIQm9SihiaQpdy8lmK9TRJpB19BE0lS0a2hmNiwsX2lmeWb2bTNbZGZ7zWytmf3KzDqGdXuY2W/CuvvM7EMz+3aMc5qZnW1mz5jZ1nCf1Wb2ZzMbloTPWGRmV5vZfDMrCc+3wcxeM7OfmllRos8p2UstNJHMdj9wEvAisIzgAZPfAg40s3OB1wkeb/IKwcz6RwO/NLMid/955IHMrB3wD+B0YC/B42Q2AROALwNnmNkn3H1uIgI3szyCR618HNhFMJP9LqAfMBb4PnALwUz2Ik1SQhPJXEMJno83xt3XA5jZYOBt4JMECWIB8AV3Lwu3fxp4ArjazH4TdmvW+AlBMnsJONfd19ZsMLOvA78H/mFm49y9MgHxH0WQzOYDR7v7nojzGcGjWXYn4DySI9TlKJLZvlmTzADcfQ1wb7g6FLikJpmF258E3iVotdUOoTeznsA3gRKCh1PWJrNwv1sIWlMjgRMTFHu/cPlyZDILz+fuPqdBwhWJSQlNJHNVAM83Ur4sXM5198aGyX8YLg+IKDsW6ADMdvfNUc43O1weEW+gUcwHqoCLzOxSM+vX1A4isSihiWSujVG6/krC5dpGtkVujxxwMSJcfjraDd/ATWGdPq0LO+Duywmu9xUCtwIbzWy5md1jZmeaWX4iziO5Q9fQRDJXdSu3R6pJHksJBpLE8kYcx43J3X9vZg8CnyG4pnYUcF74esfMZrq7rqNJsyihiQjAmnC5MFE3QTeXu28E/hi+MLNDgHuAScDVwPfaMh7JXOpyFBGA5wiuyR1vZt1TGYi7LwB+G64ekspYJLMooYkI7r6J4DpWd+AxMxvXsE54k/aXEzV4w8w+bmafMrOCBuX5wKfC1VWJOJfkBnU5ikiNqwhGPn4WeM/M3gE+Ihg8MphgXsnCcLkpAeebCPwa2GVm84ENBBMPHw4MILih+sYEnEdyhBKaiADg7hXA58zsPuBC4DCCpFNMkGzuBx4lmD0/ER4naBEeTTB7/gyCEZirCa6n3ebuWxJ0LskB5u6pjkFEslA4I/6dwN8SPdAkmceWzKUWmogk21ERD+K81t1XtvRAZvYH6p6HJlKPEpqIJNvI8AXBZMMrW3Gsc4BurQ1IspO6HEVEJCto2L6IiGQFJTQREckKSmgiIpIVlNBERCQrKKGJiEhWUEITEZGs8P8BPZGqqJOY5CsAAAAASUVORK5CYII=\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],'*',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>"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}