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?

compmech-project03/CompMech03-IVPs_project.ipynb

Go to file
 
 
Cannot retrieve contributors at this time
1008 lines (1008 sloc) 132 KB
Raw Blame
Open in GitHub Desktop
{
"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": 23,
"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": 24,
"metadata": {},
"outputs": [],
"source": [
"def euler_step(state, rhs, dt):\n",
" '''Update a state to the next time increment using Euler's method.\n",
" \n",
" Arguments\n",
" ---------\n",
" state : array of dependent variables\n",
" rhs : function that computes the RHS of the DiffEq\n",
" dt : float, time increment\n",
" \n",
" Returns\n",
" -------\n",
" next_state : array, updated after one time increment'''\n",
" \n",
" next_state = state + rhs(state) * dt\n",
" return next_state"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"def heun_step(state,rhs,dt,etol=0.000001,maxiters = 100):\n",
" '''Update a state to the next time increment using Heun's method.\n",
" \n",
" Arguments\n",
" ---------\n",
" state : array of dependent variables\n",
" rhs : function that computes the RHS of the DiffEq\n",
" dt : float, time increment\n",
" etol : tolerance in error for each time step corrector\n",
" maxiters: maximum number of iterations each time step can take\n",
" \n",
" Returns\n",
" -------\n",
" next_state : array, updated after one time increment'''\n",
" e=1\n",
" eps=np.finfo('float64').eps\n",
" next_state = state + rhs(state)*dt\n",
" ################### New iterative correction #########################\n",
" for n in range(0,maxiters):\n",
" next_state_old = next_state\n",
" next_state = state + (rhs(state)+rhs(next_state))/2*dt\n",
" e=np.sum(np.abs(next_state-next_state_old)/np.abs(next_state+eps))\n",
" if e<etol:\n",
" break\n",
" ############### end of iterative correction #########################\n",
" return next_state"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"def simplerocket(state,dmdt=0.05, u=250):\n",
" '''Computes the right-hand side of the differential equation\n",
" for the acceleration of a rocket, without drag or gravity, in SI units.\n",
" \n",
" Arguments\n",
" ---------- \n",
" state : array of three dependent variables [y v m]^T\n",
" dmdt : mass rate change of rocket in kilograms/s default set to 0.05 kg/s\n",
" u : speed of propellent expelled (default is 250 m/s)\n",
" \n",
" Returns\n",
" -------\n",
" derivs: array of three derivatives [v (u/m*dmdt-g-c/mv^2) -dmdt]^T\n",
" '''\n",
" \n",
" dstate = np.zeros(np.shape(state))\n",
" derivs = np.array([state[1],(u/state[2])*(dmdt),-dmdt])\n",
" \n",
" return derivs\n",
" return dstate"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
"#Initial conditions\n",
"y_0 = 0\n",
"v_0 = 0\n",
"m_0 = .25\n",
"\n",
"m_f = .05\n",
"dmdt = .05\n",
"N = 500\n",
"dt = 4/N"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"#Initialize array\n",
"num_sol_heun = np.zeros([N,3])\n",
"num_sol_euler = np.zeros([N,3])"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [],
"source": [
"#IC's\n",
"num_sol_heun[0,0] = y_0\n",
"num_sol_heun[0,1] = v_0\n",
"num_sol_heun[0,2] = m_0"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [],
"source": [
"num_sol_euler[0,0] = y_0\n",
"num_sol_euler[0,1] = v_0\n",
"num_sol_euler[0,2] = m_0"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [],
"source": [
"for i in range(N-1):\n",
" num_sol_heun[i+1] = heun_step(num_sol_heun[i],simplerocket,dt)\n",
"for i in range(N-1):\n",
" num_sol_euler[i+1] = euler_step(num_sol_euler[i],simplerocket,dt)\n",
"\n",
"m=num_sol_euler[:,2]"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"According to the graph, the solutions converge\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAn8AAAEaCAYAAABgnzGPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdeXhU5dnH8e9NEsIadlDCKpugkVUKYgERFBQ0ghYXUN4q2qq4YC1gW+uCS7WuuFWqYgUqKogKVsUliCwiIIiogAiJIigQ2UNCkuf9Y2YwhGQyM5lJJpnf57rmOuSc55y58zgmd57VnHOIiIiISGyoUt4BiIiIiEjZUfInIiIiEkOU/ImIiIjEECV/IiIiIjFEyZ+IiIhIDIkv7wAquoYNG7pWrVqFdO+BAweoWbNmeAOqxFRfwVF9BU91FpzS1NfKlSt3OucahTkkEQmAkr9SatWqFStWrAjp3rS0NPr37x/egCox1VdwVF/BU50FpzT1ZWbp4Y1GRAKlbl8RERGRGBK1yZ+Z3Wtmzvv6k59yl5rZIjPbY2b7zWyFmV1nZn6/t1DvExEREanIojLRMbNTgT8DfrcfMbMngRlAD2ARsABoDzwBvGZmceG8T0RERKSii7rkz8wSgWnAT8AbfsqNAK4FtgOnOOeGOucuANoBXwMXANeH6z4RERGRyiDqkj/gLqAT8Adgj59yk7zHCc65jb6TzrmfgD96v5xYRDduqPeJiIiIVHhRleCY2W+AW4CZzrm3/JRrBnQHcoBXC193zi0EtgLHAb1Ke1+4zf18K33u/5Ax7xygz/0fMvfzrZF6KxEREZGjRE3yZ2bVgBeBTODGEop39R7XOeeyiinzWaGypbkvbOZ+vpVJc9aydbfn7bfuzuLmWav569y1kXg7ERERkaNE0zp/9wAdgIudcztLKNvae/S3TlRGobKluS9sHnx3PVmH844654DpyzxvOzk1JRJvKyIiIgJESfJnZqcBNwFznXOzArillvd4wE+Z/d5j7TDcdxQzuxq4GqBJkyakpaX5edzRfC1+RZm+LIMaB7dzWtOEgJ8XS/bv3x9UXcc61VfwVGfBUX2JVEzlnvyZWXXgBWAvnlm4Ad3mPfpdCiaM9x3FOfcs8CxAjx49XDAr3Ccv+9BvAjjjmzxuu3RQacKrtLT7QnBUX8FTnQVH9SVSMUXDmL978ayxN945ty3Ae/Z5j7X8lPFd21fgXKj3hc2tZ3c4koEW5UBOniaAiIiISMREQ/J3AZAPXGFmaQVfwGBvmT96z/3b+/UW77Gln+c2L1S2NPeFTWrXZC7r1cJvmb+8rskfIiIiEhnRkPyBJ45+RbyaeK+f4P26h/frz73Hk7zdxkU5tVDZ0twXVpNTU6ieUHzVH8jJ0+xfERERiYhyT/6cc62cc1bUC8/SLwC3es918d7zPbAKqApcVPiZZtYPaIZnF4+lBd4rpPsi4b7hp/i9Pn1Zhrp/RUREJOzKPfkrhfu8x3+YWVvfSTNrDDzl/fJ+51x+mO4Lq9SuyVT1N/gPdf+KiIhI+FXY5M859xrwNJ7dONaa2VtmNgfYiGd7uLnAE+G6LxLGpFT1e13dvyIiIhJu5b7US2k45641s0+A6/CMCYwDvgGeB54urvUu1PvC7bSmCRyscdyRBZ6LosWfRUREJJyiOvlzzo0BxpRQZiYwM4Rnh3RfuE1OTWH2yh/IOlx8vjl9WQY9WtYntWtyGUYmIiIilVGF7fatTEqa/AEa/yciIiLhoeQvCqR2TWZUCWv/afyfiIiIhIOSvyhR0tp/oOVfREREpPSU/EURdf+KiIhIpCn5iyKBdv+q9U9ERERCpeQvykxOTSkxAbx51molgCIiIhISJX9RqKTxfw4Y/4oSQBEREQmekr8oVdL4v3yn8X8iIiISPCV/USq1a3KJs381/k9ERESCpeQvimn2r4iIiISbkr8oFujs38umLi2jiERERKSiU/IX5QKZ/bt4U6YSQBEREQmIkr8KIJDdPxZvytT4PxERESmRkr8KQuP/REREJByU/FUQGv8nIiIi4aDkrwKZnJpCnzb1/ZbR+D8RERHxJ6jkz8zqmlmqmd1pZs+Y2ctm9rT36/PNrG6kAhWPGWN7U8LwPxZvyuSvc9UFLCIiIseKL6mAmcUBw4Frgd8C5rtUoJjzHc3sY+Ap4HXnXF4YYxWvBy/qwk2zVvstM31ZBj1a1ie1a3IZRSUiIiIVgd/kz8wuAe4HmuFJ9nYCy4CvgExgL5AENAA6Ab2A/kA/4Hszm+icezlSwceq1K7JrEjPZPqyDL/lJsz+QsmfiIiIHKXY5M/MFuNJ5nYAjwEvOufWlPRAM+sCjAEuAWaY2TjnXJ/whCs+k1NT2LxjP4s3ZRZbJjs3n8umLmXG2N5lGJmIiIhEM3+jx9oA44EWzrnxgSR+AM651c65m4DmwC3e50gEzBjbO6AJIBr/JyIiIj5+kz/n3GPOuZxQHuycy3HOPQqcEFpoEogZY3uXuAD09GUZSgBFREQE8JP8OecOhOMNnHMHw/EcKV4gC0BPX5ahHUBERERE6/xVBqldk0vs/gW4edZqJYAiIiIxLuDkz8zqm1k3M6tf6PzxZjbNzD43s9fNrHP4w5SSBDL+zwHjX1ECKCIiEsuCafmbBHyGZyIHAGZWFfgEGA10Bs4HPjIzrS9SDgIZ/5fvtAewiIhILAsm+TsD2Fxo1u9IoDWwEBgMPAnUBa4PW4QSlEDG/x3IydMEEBERkRgVTPLXDPi20LmheHoTr3LOveecGwdsBoaEKT4JUmrXZEb1alFiOc0AFhERiU3BJH/18OzwUVBvYINz7rsC5z6nQNewlL3JqSlKAEVERKRIwSR/WXi2cQPAzJrjaQ1cXKhcNpBY+tCkNJQAioiISFGCSf6+AU4vMNv3Ujxdvh8XKtcM+CkMsUkpTU5NKXECCGgNQBERkVgSTPL3ElATWG5mrwB3AfuBN3wFzCwR6AasD2eQErr7hp9CFSu53ITZX0Q+GBERESl3wSR/TwMz8WzXdiGQA4x1zu0pUGYYngRxYdgilFJJ7ZrMw7/rQkn5X3ZuPpdNXVomMYmIiEj5CTj5c87lO+dGAW2A04Bk59wrhYp9B1wEvBi+EKW0Ursm88jILiW2AC7elKkEUEREpJIrNvkzs4vNrHbh8865zc65Zc65vUVcW+Wcm+2c2x7uQKV0Am0BVAIoIiJSuflr+ZsJ/Gxm88zsSjNrVFZBSWT4WgBLsnhTpmYAi4iIVFL+kr8JwGo8CzY/C/xoZh+Z2Q1mVvIaIhKVUrsml7gHMGgJGBERkcqq2OTPOfegc643nqVbbsCzpEsf4FFgs5mtMLNJZnZi2YQq4TJjbG8lgCIiIjGqxAkfzrltzrknnXNnAk2A3wPzgU7APcA6M/vazCabWY/IhivhEkwCqDUARUREKo/4YAo7534BpgHTzKwmcC4wHE/X8G3AJDP7AZgDzAU+ds65sEYsYTNjbG/a3Tafw/n+y908azXg6TIWESkPK1eurFelSpUr4uLiLnPONYQS56+JxBwz2++cW5Sbm/sq8FH37t2LzMGCSv4Kcs4dAF4BXjGzqsAgPIngMOBGPF3Ffwcmh/oeEnkPXtSF8a+sJt9Piu6A8a8oARSR8rFy5cqq8fHx0+rXr5/SsGHD/dWqVdtlptxPpCDnHLm5uXF79+4duGPHjkFZWVn/Xrly5T+KSgBDTv4KvWEOnq7g+WZWBegHXAD8HI7nS+T4krmbZ63GXxNtvlMLoIiUm0uSkpJSmjVrlqmkT6RoZkZCQkJegwYN9tSpU6fKt99+e9WBAweWAx8WLhvMDh8B8S4G/ZFz7gbn3LPhfr6EX6BLwDjgplmrNQlERMpUQkLCOfXr189R4icSmPj4+PxGjRoRHx9/UZHXQ3momR0HNAWqFVfGObcklGdL+UjtmsyK9EymL8sosez0ZRn0aFlfLYAiUlZOqVmz5sHyDkKkIklKStpvZr8t6lpQyZ+ZXQTcBbQvoagL9tlS/ianpgAElACqC1hEyopzrnpcXNyB8o5DpCKJj4/Pdc7VLfJaoA8xs4uBGXhmWO0BtgD7wxGgRI9AE0BNAhGRsqQuX5Hg+Pt/JpjWudu8xxuBp51zuaUJSqJXoAlgvoMJs79Q8iciIlKBBDPhox2wxDk3RYlf5Tc5NYVRvUrexS87N592t83XQtAiIiIVRDDJXybwfaQCkegTaAJ4ON/TBawEUEREJPoFk/y9B5waqUAkOgWaAPq6gEVEpOwlJyenmFn3efPm1fZXrmfPnh3MrPvjjz/eoKxik+gTTPL3dyDJzP5hZnGRCkiij7qARUREKo+AJ3w45zLM7HTgDWC4mX0A/AAUuTOsc+7e8IQo0WByagqbd+xn8aZMv+UO53sWgl6Rnnlk4oiIiIhEj2CWejHgejwTP+KANlDkjmDmPa/kr5KZMbY3l01dWmICCFoIWkREJFoF0+07ERiHJ7GbDzyGJ8Er/LoHJX6V1oyxvQPqAgbPQtDqAhYRqRg+/PDDmkOHDj2hSZMmpyQkJHSrV69e5wEDBrR99913axUuu379+qpm1j05ObnYLh4z625m3f2dnzp1ar0uXbqcWKNGja41a9bs2rt37/ZFvZ+EVzDr/F0JHAROd86tjlA8UgEE2gXs2wtYXcAiItHt73//e5O77767GUCnTp0OduvWbf+2bduqLly4sM7ChQvrPPDAA+m33HLLznC+50033dR0ypQpx3fr1m3/GWecsefrr7+uvmzZstrDhg1r//bbb68fOHCgdnWJkGCSv2TgIyV+AsF3AW/esZ8ZY3uXQWQiEmvMOKZ1qaJwjpXlHcNrr72WdNdddzVr1KjR4ZdffnnTgAEDjiRd7733Xs0RI0a0mzhxYotBgwbtO+WUU7LD9b7Tpk1rnJaW9vVvf/vbgwB5eXmMGjWq5csvv9zw9ttvbzpw4MCN4XovOVow3b5b8bT8iQDBdQEv3pTJZVOXRjgiEZHYNmzYsPa+btWiXp999tkxXap33XVXU4AnnnhiS8HED+Css846cPPNN2/Lzc21KVOmNApnrH/+85+3+hI/gLi4OP75z39uBVi5cmXt7Oxs7ekXIcG0/M0Crjazms45NcUKEPhWcOBJAP86d626gEVEIuT000/f27hx48PFXV+4cGGdXbt2Hfndv23btvgvv/yyZq1atfKGDx++t6h7zjzzzH133303K1asCOtYvBEjRuwpfC45OTk3KSkpb+/evXE//fRTXIsWLbSjWAQEk/zdDfQH3jSzq51zmyITklQ0wSSA6gIWkXCLhq7TaDFhwoTtQ4cO3Vfc9Z49e3bYtWvXkSRuw4YNVZ1z7N+/Py4hIcFv93lmZmYwOUOJ2rZtm1PU+Vq1auXt3bs3LisrK5jeSQlCMP8h3wRygTOAr83sO4pf5885584OQ3xSQQTbAjjo4TQWjO8f4ahERMSf3NxcA0/CddZZZ+32V7ZBgwYBt8Ll5eWVWCYuTvtFlJdgkr+Bhe5r730Vpaj1/6SSm5yaQo+W9Zkw+wuyc4tc+/uIjT8foN1t83nwoi5aC1BEpJyccMIJOQDx8fFu9uzZWwK9LzEx0QEcPHiwyNa5jRs3Vg1LgBIRwSR/gyIWhVQaqV2TSe2aTLvb5nPYf/7H4XwY/8rqI/eJiEjZat269eF27dplbdy4sfq8efNq++syLuj444/PTUhIcLt3747/8ccf45s2bXpUq+Drr79eJzIRSzgE3J/unPsgmFckg5bo9+BFXagSwDytfAcTZn8R+YBERKRIt99++48AV155Zes5c+YkFb5+6NAhmzFjRp3333+/pu9cYmKi69Gjx36AW2+9tWl+/q9/7b/77ru1/vGPf+gv+iimwZQSEaldk3n4d11IjC/5I5adm0+72+ZrNxARkXIwatSo3X//+99/2LVrV8KIESPatWrV6uQBAwa0HTx48AmnnHLKiY0aNeo8atSotqtWrapR8L4777xza0JCgps5c2ajdu3anTRkyJATUlJSOp5zzjkdrrjiip/L6/uRkin5k4hJ7ZrM+slDaNe4ZollD+d7dgP569y1ZRCZiIgUdMcdd/y0aNGir373u9/tzM/PZ8mSJUmLFi2qs3fv3viePXvue+ihh9KvuOKKo1b1HzRo0IG33nprQ+/evfdt3769alpaWh2AJ554YvNjjz32Y/l8JxKIYsf8mdnHwETn3JJQH25mfYD7nHN9Q32GVHwLxvfXbiAiIhG0devWgP5yXr58+frirvXu3Turd+/e6cG879lnn73/7LPP3lDUNedckUvwFHfeJ9DvRULnr+XvRGCRmS0ws5FmlhjIA80s0cwuMbP3gY8pfkawxJBgdwMZ9HBaZAMSERGJUf5m+7YD7gSuBQYA+8xsMbAU+BrYBewFkoAGQCegN9AHqIVnTcDHgDsiFLtUMJNTU9i8Y39ALYBaCkZERCQyik3+nHN7gJvM7HFgHHAFMAQY7Od5BmQCDwFPOee2hC9UqQxmjO0dcBewbxzgivRMbQknIiISJiWu8+ec+w642cxuA/rh2eKtC9AYqAPsBn4GVgEfAYucc9mRClgqvhlje/PXuWsD2g0Eft01ZGDdSEYlIiISGwJe5Nk5lwW8432JlEowu4GAJwFcWd/o3z/ysYmIiFRmWupFyk0wS8EAfJ3pNBFERESklJT8SblbML4/fdrUD6isbyKIFoQWEREJjZI/iQrBLAWjBaFFRERCp+RPosbk1JSAE0DwjANUAigiIhIcJX8SVSanpvDoyMD2BAZPAnjZ1KURjkpERKTyUPInUSfYiSCLN2VqHKCIiEiAlPxJ1ApmIojGAYqIiARGyZ9EtWAmgoDGAYqIiJQk4OTPzK4ys+qRDEakKKFMBFE3sIiISNGCafl7FvjBzB4ys3aRCkikKJNTU7j6lKoBTwRRN7CIxJLk5OQUM+te0mvevHm1S/tevmeFI24pHwFv7wbMA4YANwM3mtn7wJPAPOeci0RwIgWd1jSB2y4dxGVTl7J4U2ZA90xflsHmHfuZMbZ3hKMTESl/p59++t7GjRsfLu56cnJysdckdgSzt+95ZtYc+CPwe+AsYBDwvZk9AzznnNsRShBmlgD0Bc4B+gAtgQbADmAp8IRzLs3P/Zd64zoFiAO+AV4AnnbOFbtxbKj3SfmaMbY3f527lunLMgIqv3hTJoMeTmPB+P6RDUxEpJxNmDBh+9ChQ/eVdxwS3YKa8OGc+945dxvQHBgNLANaAPcAGWb2kpmF0sTSD3gfGI8n8VsJvA5kAiOAj8zsrqJuNLMngRlAD2ARsABoDzwBvGZmceG8T6JDsOMAtS2ciIiIR0izfZ1zh51zM5xzfYCuwHNALnAp8ImZrTSz35tZYoCPzAdmA32dc8c754Y650Y651KAi4E84G9mdkbBm8xsBHAtsB04xXvfBUA74GvgAuD6wm8W6n0SXYJdENo3DrDj3/6nJFBEYtq8efNqm1n3nj17dijq+vr166uaWffk5OSUYJ6bnZ1tDzzwQKPu3bt3SEpK6pKYmNitZcuWJ1911VXNfvzxx2N6Gx9//PEGZtZ9xIgRrbZv3x43ZsyY5snJySkJCQndBg4c2CbU70/8K/VSL865NcCdeLpLzfvqCkwFtpjZlQE840Pn3IXOuUVFXJsFTPN+OarQ5Une4wTn3MYC9/yEpzsXYKKZFf4+Q71PooxvQehA1wMEyDqcr8kgIiJhlpmZWeW0005rP2HChBYbNmyoftJJJx3s37//ntzcXHvuueea9OjRo+P69eurFnNvfI8ePTrNnTu3QceOHQ8OHDhwt7+xi1I6wUz4OIaZDcTTgjYUz5i5Q8BMPF2oo/CM4XvWzGo65x4vxVt97j02K/DezYDuQA7wauEbnHMLzWwrkAz0ApaU5j6JbjPG9mbu51uZMPsLsnMDG67pGzM4OTWoP2xFJBqNH9+URx45PqCyF1+8k//+N/2oc5dc0pKXX24Y0P0337yNhx/+8ahzAwa05aOP6gR0/4MPpvOnP+0MqGwFcvnll7datWpVrcGDB//yn//8J71Ro0Z5ALm5uYwbNy75mWeeOW706NGtly9fvr7wvWlpaXX69Omz96233tpUr149jbmPsKBbtsysjpndZGbfAO8CqcCPwG1AM+fcVc65Wc65YcBpwAHghlLG6VtaZluBc129x3XOuaxi7vusUNnS3CdRLtht4UB7A4tI5TJs2LD2xS3zUrt27S6Ret+VK1dWmz9/fr2mTZvmvPrqq5t9iR9AfHw8TzzxxNb27dtnffbZZ7WWL19+zJrB8fHx7rnnnktX4lc2Am75M7NueFr5Lgaq4+neTQOmAG8UNTvWOfepmc3HM2kjJGZ2HDDG++XsApdae49H//V2NN900NYFzoV6n1QQC8b3D2o5mMWbMjlh4nweHtmF1K7JEY5ORCRy/C31Ur169YglVm+++WYdgDPPPHNPrVq1jln+LS4ujp49e+7fsGFD9Y8//rhmz549j2p86dSp08EOHTrkRCo+OVow3b4rvMeDwL+BKc65LwO470CQ73OEmcUD04E6wAfOubcKXK5V4PnF2e89FlzUMtT7CsZ1NXA1QJMmTUhLS/PzKD9vsn9/yPfGomDqa2w7qJoTx0ff55VcGM+Mo5tmreZfC9YwoWeN0IOMIvp8BU91Fpyoqq+HH/7xmK7YYPz3v+nHdAUH48MPvw353jAqr6Vevvvuu0SAl156qdFLL73UyF/ZHTt2HJMTNGvWTIlfGQomKduCZ1Hn55xzu4O4byxwTTBBFfAMcCbwPcdO9jDvMdgFpkO97wjn3LN4djyhR48ern///iE9Jy0tjVDvjUXB1lf//jD386385fW1HMgJLAn8OtNxzyoqxZqA+nwFT3UWHNVX5ZWXF9jPzMLlTzrppIMdOnQobkgVACeffPKhwueqVaum7t4yFEzy1yaUnTy89wT3KQLM7DHgSjzLsZzpnNteqIjvL5taFM93reBfQaHeJxVQatdkUrsmB7Uo9MafD9Bq4nxG9WqhySAiUiklJibmAxw8eLDIsf+bNm0KdKk2AJo3b54D0KdPn33/+te/fih9hBJJwUz4eNfMxpdUyMxuNrP3ShETZvYQnkkiO/AkfhuLKLbFe2zp51HNC5UtzX1SgQW7JiB4JoMMejgtckGJiJSTli1bHgbIyMhIzM7OtsLX582bF9jMZa9hw4btBXjnnXfqHj6sFVqiXTDJ30Dg5ADKdcLTVRsSM3sAz04fu4BBzrmviinqW/7lJDM7ZuaQ16mFypbmPqngQpkN7GsF1JqAIlKZtG/fPqd58+bZ+/bti7vjjjuaFLz20ksv1X3hhRcaB/O8008//eDAgQN3Z2RkJJ577rltNm3alFC4THp6esJdd93VWMlh+SvVOn/FqIpn/HzQzOx+4FbgFzyJ35riyjrnvjezVUA34CLgP4We1Q/PuoDb8ewPXKr7pPJYML5/UN3A4GkFnLU8gwcv0oxgEYle//jHP4574YUXGhR3/bLLLsscPnz4XoA77rhj61VXXXXC/fffn/zmm2/Wa9GiRfbmzZurbdiwofr111+/bcqUKYGtm+g1a9aszYMHD263YMGCuieddFKdDh06HGzWrFnOvn374rZt21b1u+++q5afn8+f/vSnHQkJCSGPu5fSC2vyZ2aGZwHloBevNLO7gQnAbjyJXyCtbvfhWaj5H2a2xDn3rfdZjYGnvGXuL2IZmlDvk0picmoKk1NTGPRwGht/9jfx+1e+7eFWpGdqLKCIRKVPPvkkyd/1zp07H/Qlf7///e9/SUxM/PbBBx88fv369dXT09OrderU6eCrr7668eSTTz4UbPJXv379/CVLlqz/17/+Vf+///1vg3Xr1tVYt25djaSkpLzGjRsfvvTSS3dccMEFu2vUqKHEr5yZvzkchcbuDcSzmHNx3bDxeBZjbgq85pwbGXAQZucBb3i/XAGsK6boN865+wvd+xSeLdkOAe8Dh/F0OycBc4ELnXPHTDgJ9b7CevTo4VasWFFSsSJpplxwIlVfwawJ6FM9oQr3DT8lqlsB9fkKnuosOKWpLzNb6ZzrEUjZNWvWbOncuXOl2xFDJNLWrFnTsHPnzq0Kny+p5W9ggX87PIld0xLu+QL4c1DRQcGNWXt4X0VZCByV/DnnrjWzT4DrgH54tpn7BngeeLq41rtQ75PKx7c13J9eXUNufmB/kPr2B1YroIiIVDQlJX+DvEcD3sOznds/iymbA2x1zn0XbBDOuWnAtGDvK3D/TDx7CpfJfVL5+JaECbYVcPqyDGav/CHqWwFFRER8/CZ/zrkPfP82s8XAwoLnRCobXyvghNlfkJ0bWOOvWgFFRKQiCXipF+fcbwuPtxOpjHxLwozq1SKo+6Yvy6Dj3/7H3M+3RigyERGR0gtmnT+RmDI5NSXoBNDXCqh1AUVEJFoV2+1rZrd5//m0c+6XAl8HxDl3b6kiE4kCk1NT6NGyflD7A4PWBRQRkejlb8zfZDwzfF/Ds+iy7+uSmLeckj+pFELZHxh+XRdw0pwvNCFERESihr/k7148SdzOQl+LxKRQWwE1IURERKJJscmfc+6v/r4WiUWhtgKCloUREZHooAkfIiGYnJrCoyO7ULNqXFD3+VoBL5uqbaNFRKR8KPkTCVFq12TW3TWYR0d2ITE+uP+VFm/KpN1t87UsjIiIlLmAf2OZ2R/NLMfMzvVTZqi3zFXhCU8k+vnWBezTpn7JhQvwTQjR2oAiIlKWgmmuGA5kAv/zU+Z/3jIXliYokYpoxtje6goWEZGoF0zydyKw1jlX7J5Xzrk8YC3QqbSBiVREvq7gYBeHBk9XcKuJ87VAtIiIRFQwyV8j4KcAyv0MNA4tHJHKIdQJIeCZFTzo4bTwByUilYqZdQ/2NWLEiFahvNfJJ5/c0cy6f/zxxzVKG/eePXuqmFn3GjVqdC18rV69ep3NrPu2bdv8LUUX9T777LNqZta9Xbt2J5V3LEUJpnL3AM0DKJcM7A8tHJHKw7cszNzPtzJh9hdk5xbbaH6MjT8foNXE+Yzq1UJrA4pIkYYPH76r8Lmff/454ZNPPkmqXv97f0QAACAASURBVL16/pAhQ34pfL1Pnz76/SxBJX+fA2eYWRvn3KaiCphZG+A04ONwBCdSGZR2bcCZyzJ4eKS2iRORo82ePXtL4XPz5s2r/cknnyTVq1cvt6jroZozZ86mAwcOVOnQoUN2uJ4p5SeYbt9pQAIw18zaFb5oZm2BuUCct6yIFBBqV3A+nlnBWhpGRMpL+/btc7p27XqoRo0a2umrEggm+ZsFvA2cBKwzsw/N7Cnv6wPgK++1d51z0yMQq0iFV5q1AbU0jEjZmb4svX7Pe95PaT1xfvee97yfMn1ZenBrOVUAU6ZMaXDqqad2SEpK6hIfH9+tXr16nTt06NBpzJgxzTdu3Fi1YFl/Y/6ysrLsjjvuaHLyySd3rFmzZtfq1at3bdeu3Uk33nhj0127dgU/8LkIubm5jBo1qoWZdW/fvn2nTZs2JRS8vmTJkurDhg1r3bhx41MSEhK61a9fv/OAAQPavvnmm7ULP6t9+/adzKz7G2+8ccw1n0svvbSFmXW/5ZZbjved27t3b5Vbb731+A4dOnSqXr1616pVq3Zr3LjxKd26dTvx5ptvbnr48OGAvpft27fHde3a9UQz6z5s2LDWWVlZ9oc//KGZmXW/+uqrmxV33zPPPFPfzLr36tWrfUBv5EfAv32ccw7Pci9Pe0/1B/7gfZ3hPfc0cEFpgxKp7HxrAz46sgvxVSyoe31LwygJFImM6cvS698976uWP+/LruqAn/dlV7173lctK1MCePXVVze74YYbWq1Zs6Zmp06dDg4ZMuSXlJSUg9nZ2VVefPHFxitWrKgeyHP27NlTpXfv3h3uvPPOZlu2bKnWq1evvf3799+za9eu+Mcff/z4bt26dSycqAVr3759Vc4666y2M2bMaNS7d+99y5Yt+6ZNmzZHMq1nn322Xr9+/TrOmzevfv369XMHDx78S6tWrbLT0tLqnH/++e3/8pe/HFfweRdffPEugBdeeKFhUe+XlZVl8+bNq29mjB07dhfA4cOH6du3b/t//vOfTX/66aeqvXr12nf22Wf/csIJJxz64Ycfqj766KPHHzx4sMScat26dYm9evXquHr16prXXHPNT2+88cbm6tWru/Hjx/8cFxfHK6+80vDgwYNF/lKYOnVqI4A//OEPPwdee0ULajaNcy4HuM7M7gbOBFp6L6UDHzjntpc2IJFY4hsPeNnUpSzelBnUvb4kcEV6piaFiITR4x9sTM7OzT/qF3l2bn6Vxz/YmDyqV8vg/keNQrt27Yp7/vnnG9epUyfv008//apDhw45Ba+vWrWqWlJSUl4gz7ruuuuarVmzpmaHDh2yFixYsKF58+a54EkKzzvvvBM+/vjjOqNHj261ZMmSjaHEunXr1vjBgwe3+/LLL2ucf/75mbNmzdqSmJh4pOt5/fr1VW+88cZWubm59sADD6TfeuutO33XXnnllaRRo0a1ve+++5J/+9vf7h88ePB+gLFjx+669957k9999926v/zyS5V69eodNRtvxowZdfft2xfXs2fPfSeeeGIOwJw5c+qsWbOmZrdu3fYvWrRoQ8Hu77y8PN55551a1apV89sl/tFHH9UYMWJEuz179sTfe++9GZMmTdrhu9a+ffucM844Y/f7779f99///nf9G2644ajJPJ9++mn1VatW1WrcuPHhyy67bHcodVlQSNu7Oee2O+dmOOfu9b5mKPETCZ1vgehgWwHBMylErYAi4bNjX3bVYM5XNDt27IjLy8uzNm3aZBVO/AC6det2qG3btiX2Ye7cuTPu1VdfbQjwxBNPpPsSP4A6derkP//88+lVq1Z1S5cuTVq6dGlALYkFffHFF4m9evU68csvv6xx/fXXb587d+7mgokfwKOPPtro0KFDVfr06bO3YOIH8Lvf/W7vyJEjdzrneOihh5r4zjdv3jy3b9++ew8dOlTlxRdfrFf4fV966aUGAKNGjTqSgG3fvj0eoE+fPvsKj3uMi4vj3HPP3V84toJmzJhR55xzzumQlZVV5cUXX9xUMPHzGTdu3M8AU6dOPWa5vMcee6wxwOjRo3ckJJSqIRXQ3r4iUSO1azLf3ntOSAtEa5cQkfBpVDvxmITI3/mKpl27djkNGjTI/fzzz2uNGzcu+csvv0wM5TlpaWk1c3JyrGXLltkDBw48UPh6mzZtDvfp02cvwPvvv1/s+LqifPDBB7X69et34rZt2xIffPDB9ClTphT51+2SJUtqA4wePfqYZW8Arr766p0An3766VHvf/nll+8EmD59+lFdvxkZGfGLFy+uU6NGjfwrrrjiyFI5vXv3PmhmTJs2rfHDDz/cMJh1CO+///5Gl19+edvq1avnz58/f8OoUaOKbLk777zz9rVt2/bQl19+WaPg2MrMzMwqc+fOrR8fH+/GjRu3s6h7gxV08mdmHczsSTNbZ2a7va91ZvaEmZ0YjqBEYtnk1BS23H9u0HsFg3YJEQmHG85stzUxvspRXYGJ8VXybzizXaVoXo+Li2Pq1Kmbk5KS8p544onjUlJSTm7YsGHnQYMGtXnwwQcb7tmzJ6Dc4Pvvv68K0Lx582KXf2nVqlU2wNatW4NqNR0zZswJu3fvjp88eXLGn/70p2ITnp9++qkqQNu2bYuMoWPHjtkA+/bti9u3b9+R7+viiy/eU7du3dxVq1bV+uqrr47E9u9//7tBXl4eQ4YM+SUpKenIZ6Bnz55ZkyZN2nrw4MEqt9xyS8umTZt2btmy5ckjRoxoNX369Lp5eUX3km/evDlx0qRJLcyMd955Z8OAAQOOSZILuvrqq38CmDJlypHWv6effrphVlZWlbPPPnt3y5YtA5tVUoKgkj8zGwOsxjPJoyOQ5H11BK4FVpvZFeEITCTWhbpXMHi6gv/vnQPqChYJwaheLTP/NrRTeuPaiTkGNK6dmPO3oZ3SK8N4P58LLrhgb3p6+hfPPPPM5ksuuWRH/fr1D3/wwQd1//znP7ds27btyatWrapW0jM880DBrPjhKr4ywfItYP3II48cv2bNmmJbJgOJoSjVqlVz559/fqZzjqlTpx5p/Xv55ZcbAPzf//3fMQnnPffcs33jxo1r77333oyhQ4dmZmdnV5kzZ06D0aNHt+nevfuJ+/fvPyaI448/PufUU0/dn5eXx7hx45qXlFhfc801mbVr186bN29evR07dsQBPP/8840ArrvuulJP9PAJOPkzs1OBqUBV4HVgKJ6krxNwLjAbzzqAU71lRaSUCi4NE2wS6ND6gCKhGtWrZebyvwxcu/n+c1cu/8vAtZUp8fOpU6dO/jXXXJM5c+bMjA0bNny1cePGL84888zdO3fuTBg3blyJO3q1aNEiByAjI6PY5Cw9PT0RIDk5Oagu88cee2zrTTfdtO3nn39OGDBgwInLly8vcszgcccdlwOwcePGImP45ptvEgFq166dV7t27aNac6+66qpdAK+88kqD/Px8Fi1aVGPjxo3Vk5OTc4YMGVLkTiitW7c+PGnSpB1vvfXW5u3bt3+xcOHCr1u2bJm9Zs2ampMnT25SuHy1atXchx9+uOH000/fu3z58tr9+/dv72/5m6SkpPyRI0fuPHToUJUnn3yy4Ztvvln7u+++q9auXbus4mIKRTAtf7d6y49yzl3onHvbObfeOfeNc+5/zrmLgFF4ZhD/KVwBisivSWAo4wF96wOqO1hE/GnTps3hu++++0eAb775psQ9fPv373+gatWqLj09PfHDDz+sWfj6li1bEpYsWZIEMHDgwH3BxvPII4/8eNttt23NzMyMP/vss9svWrTomJhOO+20fQAzZsxoUNQzpk6d2gDgN7/5zTHvf/rppx9s37591o8//lj17bffrv3cc881ABg5cuTOKlUCS4/69u178Morr/wZYO3atUXWWa1atdyCBQu+HThw4O7Vq1fX7Nu3b/vt27cXmwDefPPNP1epUoVp06Y1evLJJxsDXHXVVcdMECmNYJK/04GVzrn/FlfAe+0zoG9pAxORY4W6S4jP9GUZSgJFYtwXX3yROGXKlAZFdUG+/vrrdQGaNm1aYktdw4YN8y688MKdAOPGjWvx448/HpkEsXfv3iq///3vW2ZnZ1vv3r339u7dOyuUWO+5557t99xzT8aePXvizznnnPbvv//+UUnmTTfdtKNatWr5ixYtSnrkkUeOmrwxZ86cpJdffrmRmXHLLbf8VNTzL7300p0Azz77bMM33njjqLX9CnrttdeSXn/99aTc3NyjzmdnZ9uCBQvqADRv3rzYOqtWrZp7++23Nw0dOjTzq6++qtG3b98O33//fZGTRjp16pTTr1+/Penp6Ynvvfde3Vq1auVdc801RU5oCVUwyV8DYEMA5TYClWYhTJFoU5pdQnymL8vghInqDhaJRdu2bUu44YYbWjVu3LhLly5dThw2bFjrIUOGnHDCCSec9OCDDzatWrWqmzx58g+BPOupp576oXPnzge++uqrGu3bt08ZOHBgmyFDhpzQunXrlI8++qhOixYtsl966aUtpYn3tttu2/HQQw+lHzx4MO78889vP2/evCMzdzt06JDz6KOPpsfFxbnx48e37NixY6fzzjuvdffu3TtceOGF7XJycmzixIlbfWv8FTZ27NjM+Ph499Zbb9XfvXt3/Kmnnnpkbb+Cli9fXnP48OHtGjRo0OW0005rf95557UeOHBgm+Tk5FMWLVqUdNxxx+VMmjSpyATTJyEhgblz526+8MILd23cuLF63759O2zevLnIdVt8y76AZ/xjnTp18osqF6pgfnP8ArQJoNwJ3rIiEkGl2SUEtGewSKzq1q1b1u233/5Dnz599u7atSv+gw8+qLto0aI6Zsbo0aN3fPbZZ+vOP//8gLpp69Spk79kyZL1t99++w8tW7Y8tGTJkqSPPvqobt26dXPHjRu3beXKlV8X3I0jVDfffPPOp556anN2drZddNFFbWfPnp3ku3bNNddkLly48OuhQ4dm7ty5M/7tt9+u991331Xr37//nrlz52649957i12HuGnTprn9+vXb4/u64Np+BV1yySW/3HDDDds6dOiQ9d1331V79913661cubJWkyZNciZOnLh19erVXwUyEzcuLo5Zs2ZtGT169I4tW7ZU69u3b4f169cfMxN68ODB+6pWreoAbrzxxrB2+QJYoDNxzGwOcD4w3Dn3RjFlhgFvAK8750aELcoo1qNHD7dixYqQ7k1LS6N///7hDagSU33599e5a5m+LCPk+6snVOG+4aeQ2jU5jFFVLPqMBac09WVmK51zPQIpu2bNmi2dO3cOy/pmIhXB008/Xf/aa69t3atXr31Lly4NpNe1SGvWrGnYuXPnVoXPB9Py97D3+KqZPW9m/cyshZk19/77OeA1PA0KDxf/GBGJBN/6gKFMCgEtFC0iEg2ysrLsoYceOh7gxhtv9NuVHKqAkz/n3CfATYABVwAfApuBLd5//5/3eTc55xaHPVIRCYgvCTyjeWiTQrRQtIhI2XvggQcajRgxolXHjh1P2rRpU7Xf/OY3+y6++OI9Jd8ZvKBGizvnpgA9gelABpAL5Hn//R+gp3PuiXAHKSLBu+Kkamy5/9xS7RncauJ87RssIlIGPvzww9pz5sxpsGfPnrjzzjsv8/XXX/8uUu8V8N50Ps65z/G0/IlIBZDaNZnUrskhjwn0dQdPmvNFzI8JFBGJlHfeeSdiyV5hoa0TISIVTmn2DAaNCRQRqSyU/InEmNLsGQy/jglUd7CISMVUbLevmT1biuc659w1pbhfRCLI1xU89/OtTJj9Bdm5wa8fqu5gKUvOOcyCH7sqEqv8LeXnb8zfVaV5T0DJn0iUUxIoFYGZZeXl5VWJj48P6y4HIpVZbm5uvJkVubOJv+RvbITiEZEoE84k8NZXV/PgRV2UBEo4fXHgwIHOderUKfIXmYgca+/evbWccwuKulZs8uecey5yIYlINCqYBP7l9bUcyMkL+hmH81ESKGF1+PDhtzMzM09NSkpS169IAHJzc6vs2LGD3NzcV4u6rgkfInKM1K7JrLtrMI+O7EJifGg/JnxJoBaMljD47969e9d+//33DbKysqoGui2pSCxxznH48OH4Xbt21f3222+TsrKy/g18VFTZoNf5AzCzWkAPoBGQ4Zz7NPRwRSRahaM7GDwLRs9clsHDI9USKMHr3r17zsqVK6/YtWvX5bt37x7tnGuIZ7cpESnAzPY75xZ4W/w+6t69e5F/KQWV/JlZbeAh4HIgwXv6ReBT7/U/ApOAC51zy0MNXkSiSzi6g/PxtATeNGs1dasncMd5JykRlIB17959N/C49yUipRBwf46Z1QDS8MwC3gss4Ni/vN4DmgEXhCk+EYkiBbuDQ10nEGB31mEtGC0iUk6CGcxzC9AV+C/Q2jk3uHAB59wmYCMwIDzhiUg0ClcSqAWjRUTKXjDJ3++AbcCVzrkDfsqlA+rLEYkBviSwNNvGwa/LxCgJFBGJvGCSvzbAcufcoRLK7QQahh6SiFREvm3j6tVIKLlwMXxJYLvb5isJFBGJkGCSv8NAYgDlmgFaiFMkBqV2Tebz289iy/3nMqpXi5Cfo2ViREQiJ5jkbwPQ1cyKTQDNrC7QGfiytIGJSMU2OTXlSBJYmjU5pi/LUBIoIhJGwSR/s4EmwL1+ykwGagFFrigtIrFncmoKm+8/t1QLRsOvSaDGBYqIlE4wP4mnAOuBm8xsoZnd4D3f0szGmtl7wB+BdcC/wxyniFRwqV2TWT95SKmTQN+4QLUGioiEJuBFnp1zB8zsLDwtgL8FTvde6u99GbAaON85lx3eMEWksgjHgtE+05dlMH1ZhhaNFhEJQlB/fjvnvnfO9QTOA/6FZ1HnD/Hs8jES6OGc+z7sUYpIpROutQJBi0aLiAQjpL19nXPzgHlhjkVEYpCvJRDgr3PXMn1ZRsjP8i0aXT2hCvcNP0UtgSIiRSi25c/MXjOzIWamzbNFpEz4ZghrXKCISOT4++k6HE/r3vdmdo+ZtSujmEQkxoVrcgholrCISGH+fqo+DfwCNAUmAt94Z/leYWY1yiQ6EYlpBZPA0uwcAmoNFBHxKTb5c85dhyfxG4lnYkc+nlm+zwPbzWyqmZ1WJlGKSEwrvHNIacei+FoDu9z5nloDRSTm+O1Pcc7lOOdedc4NAVoCf8Gz00ct4EpgkZl9bWa3mtlxkQ9XRGJduBaNhl9nCWsvYRGJJQH/5HTO/eicu8851xHoAzwH7AM6APcDGWb2hpmdb2alW7dBRKQEBbuES7tUjG8v4THvHFBroIhUeiH92eycW+qcGwscD1wBpAFxwFBgDqCfnCJSJnzrBfq6hEtLrYEiUtmVqs/EOZflnHvJOXcmMBjYiWenj0bhCE5EJBjhWioGfm0N1AQREalsSvXT0cxqmdmVZrYIeIdfkz7t8iEi5Sacs4RBy8WISOUSUvJnZmeY2X+A7cCzeMYA5gCvAkOA1mGLUEQkROGeJazlYkSkMgh4ezcza41nfN8VQAs48nN0NfACMN0590vYIxQRCYPJqSlMTk1h7udb+cvrazmQk1eq501flnFkK7pRvVowOTUlHGGKiESc3+TPu5jzRcAYPGv8mff1CzATeM45tzrCMYqIhE3hvYRnLMvAlfKZvkRQewqLSEXgb2/f5/B06z4P9POefh+4BDjeOTdOiZ+IVGQF1wws7XIxoG5hEakY/LX8/Z/3uBmYBrzgnPsh4hGJiJQxX2tgWloa7+9ucKQ7tzR8rYEGXKZuYRGJIv4mfMwAznTOtXHO3a3ET0RiQTiXiwFwaDs5EYku/vb2He2c+6gsgxERiRbhXi4Gfl1AWt3CIlKeAp7tKyISiwpOEAnXTGFQt7CIlB8lfyIiASqcCE6Y/QXZufmleqavW3j6sgzqVk/gjvNO0mxhEYmo0g9oERGJQZHuFtb4QBGJFCV/IiKlUHgXkSql3UbES+MDRSRS1O0rIhImvl1EIHwLSIPGB4pIeKnlT0QkAgouIB2ubmEtGyMi4aCWPxGRCCq8ndzMTzPID0NzoK9b+KZZqzVRRESCopY/EZEyMjk1he/uO/fI+MAwDQ88anxgx7/9Ty2CIuKXkj8RkXIQ7n2FfbS/sIiURN2+IiLlqPDagXe+tY5fDh4Oy7M1UUREiqLkT0QkSkRqfGDBhaSVCIqIkj8RkSgUqWVjlAiKiMb8iYhEuUiNDyy4dIzGCIrEjphv+TOzS4E/AqcAccA3wAvA08650m3aKSISRpEcHwi/tghWT6jCfcNP0dIxIpVUTCd/ZvYkcC1wCPgAOAycCTwBnGlmFznn8soxRBGRIkUyEfTNGL5p1molgiKVUMwmf2Y2Ak/itx3o65zb6D3fBPgIuAC4Hnis3IIUEQlApCaKgBJBkcoolsf8TfIeJ/gSPwDn3E94uoEBJppZLNeRiFQwhReSrhKulaQ5eg3Bjn/7H0t+DF+Xs4iUnZhMbMysGdAdyAFeLXzdObcQ2AocB/Qq2+hERMIj0ongs1/k0O62+dpRRKSCicnkD+jqPa5zzmUVU+azQmVFRCqsgolgOGcNH86H8bNWKwEUqUBiNflr7T2m+ymTUaisiEilkNo1mXV3DT6SCCbGl+5XQT7w4LvrwxOciERcrE74qOU9HvBTZr/3WLvwBTO7GrgaoEmTJqSlpYUUxP79+0O+NxapvoKj+gpeLNZZXeBfA6uz5MfDzPw6h/0hDuPbujsr5upOpKKK1eTPN/IlpPlwzrlngWcBevTo4fr37x9SEGlpaYR6byxSfQVH9RW8WK6z/sBt3n+HsnRMct3qMVt3IhVNrCZ/+7zHWn7K+K7t81NGRKTSKbh0DJS8vVwV4NazO5RJbCJSerE65m+L99jST5nmhcqKiMQk3/ZyRc0arp5QhYdHdtHafyIVSKy2/H3uPZ5kZtWLmfF7aqGyIiIxb3JqCpNTU4DY7iYXqchisuXPOfc9sAqoClxU+LqZ9QOa4dn9Y2nZRiciIiISOTGZ/Hnd5z3+w8za+k6aWWPgKe+X9zvn8ss8MhEREZEIidVuX5xzr5nZ03i2cltrZu8Dh4EzgSRgLvBEOYYoIiIiEnYxm/wBOOeuNbNPgOuAfkAc8A3wPPC0Wv1ERESksonp5A/AOTcTmFnecYiIiIiUBXMupHWOxcvMduB/mzh/GgI7wxhOZaf6Co7qK3iqs+CUpr5aOucahTMYEQmMkr9yZGYrnHM9yjuOikL1FRzVV/BUZ8FRfYlUTLE821dEREQk5ij5ExEREYkhSv7K17PlHUAFo/oKjuoreKqz4Ki+RCogjfkTERERiSFq+RMRERGJIUr+RERERGKIkr8gmNmlZrbIzPaY2X4zW2Fm15lZSPUY7PPMbJqZOT+vb0r3HYZXuOrLzDqY2Y1mNt3MvjGzfO/3e2FZxhFp5V1fFe3zBeGpMzNLMLMzzewhM1tmZtvMLMfMtprZa2bWvyziKAvlXV8V8TMmUhnF/A4fgTKzJ4FrgUPAB/y6D/ATwJlmdpFzLq+MnrcY+LaI89sCff9IC3N9/RG4MQriiJhoqS+vqP98QVjrrB+wwPvv7cBK4ADQCRgBjDCzu51zt0c4joiKlvryqhCfMZFKyzmnVwkvPD/QHJ4fTO0KnG8CfOW9dmOknwdM814bU951Usb1dRXwAPA7oA2Q5n3GhWUZRwzUV4X4fIW7zoABwGvAb4u4NhLI9T7vDH3GSl1fFeYzppdelflV7gFUhBewwvsD6/IirvUr8EO1SiSfV1F+cIa7vop4RqDJTETjqIT1VSE+X2X93xb4t/d5z+kzVur6qjCfMb30qsyvqBqPEo3MrBnQHcgBXi183Tm3ENgKHAf0KuvnRZto+f6iJY6SVJQ4o0k51Nnn3mOzco4jJNFSXyISPZT8layr97jOOZdVTJnPCpWN9PPOMLOHzexZM7vbzM6OooHl4a6vih5HSaIxzmj+fEHZ11k777HweLRo/G9XlGipr4Ki/TMmUqlpwkfJWnuP6X7KZBQqG+nnXV7Eua/M7GLn3NoAYoikcNdXRY+jJNEYZzR/vqAM68zMjgPGeL+cXV5xlFK01FdB0f4ZE6nU9JdWyWp5jwf8lNnvPdaO8PNWAzcAJ3mf0xQYCqzBM9PufTNLDiCGSAp3fVX0OEoSTXFWhM8XlFGdmVk8MB2oA3zgnHurPOIIg2ipL6g4nzGRSk0tfyUz7zFc++CF/Dzn3KOFTh0A5pvZAmAhnvE6k4DrSxVh6YS7vkIVLXGUJGrirCCfLyi7OnsGz1Io3wOjyjGO0oqW+qpInzGRSk0tfyXb5z3W8lPGd22fnzKReh7OuRzgPu+X5wRyTwSF/fur4HGUJOrjjLLPF5RBnZnZY8CVeNaxO9M5t7084giTaKmvYkXhZ0ykUlPyV7It3mNLP2WaFypbls/z8a2MX95dJlu8x3B/fxU1jpL43jva44yWzxdEuM7M7CE8XZM78CQyG8sjjjDyvXd511dJoukzJlKpKfkrmW/ZgpPMrHoxZU4tVLYsn+fTwHvc77dU5EXq+6uocZSkosQZLZ8viGCdmdkDwHhgFzDIOfdVecQRZtFSXyWJps+YSKWm5K8EzrnvgVVAVeCiwtfNrB+e9ay2A0vL+nkF/M57/MxvqQiL4PdXIeMoSUWJkyj5fEHk6szM7gduBX7Bk8isKY84wi1a6isAUfMZE6n0ynuV6YrwAi7k1xXw2xY43xhYRxFbI+EZv/INcF+YntcFz6y4uELn4/H85Z3nve/sylZfRTw/jcB2rAg6jlitr4r0+YpEnQF3e+/5BegeyThitb4q2mdML70q86vcA6goL+Ap7w+mLOAtYA6wx3vu9SJ+oE3zXpsWpueleq/twvPX+avAO3hW5nfeH5x/Lu96ikR9Ad2AZQVee71lNxQ8H444YrW+KtrnK5x1BpznPe/wtDpNK+Y1eEAP6wAACGBJREFUUZ+x0OurIn7G9NKrsr7KPYCK9AIuBRZ7f5keAFYC11HEfpj+fjmH+LzWwKPAEu8Py0PeH+IbgecJorWiotUX0L/AL5tiX+GII1brqyJ+vsJVZ3gWJS6xvoA0fcZCr6+K+hnTS6/K+DLnHCIiIiISGzThQ0RERCSGKPkTERERiSFK/kRERERiiJI/ERERkRii5E9EREQkhij5ExEREYkhSv5EREREYoiSP5ECzGyLmTnv674Sys4oUDatjEKMKDOLM7O1ZpZuZokFzrfyfp9byiCG4d73uj7S7yUiEouU/IkU73IziyvqgpklAReUcTxl4Y/AycAdzrns8gjAOTcHWAHcaWb1yyMGEZHKTMmfSNFWAE2BQcVcvxiojmdv00rBzGoBdwKbgf+Uczh3AfWB28o5DhGRSkfJn0jRpnmPY4q5PgbPRvQvlUEsZeUKPAnXNOdcXjnH8jawHbjKzGqWcywiIpWKkj+Ron0KfAWcb2Z1C14wsw5Ab+BdYFtxDzCzgWb2pJmtMbNdZpbtHUv3opl1LOaeamY20cxWmdl+7z3bzGypmU02s2qFyvc0s1fNbKuZHTazPWb2rZnNNLMBQX7P13qPQbX6mVkdM/vQO05vrplVL3DNzOxqM/vczLLMbIeZzTGzFDMb471nWuFnepPPGUAd4NIgvw8REfFDyZ9I8aYB1YBLCp0f4z2+UML9zwBXArnAIjytWTnA5cAKMzu9YGEzqwLMB+4DTgAWArPxJKHNgb8AdQuUHwR8AlwI/Ay8DnwI/OI997sAv0/MrB3QCfjWObcliPuae2M4A3gKGO6cyypQ5F/e18nAYuB9IAVPct2jhMe/7z2eH2g8IiJSsvjyDkAkir2EJxEbAzwNntmweJK3TOBN4Dw/9/8JSHPO7fadMDMDrsaTGD5rZic555z38unAAGAV0Nc5d6DQfacBews8fxKQAFzqnPtvwTc2swZAqyC+1/7e49JAbzCzzngS2uOBCc65BwpdTwXGAruBM51zq7znqwD/wFM//iwDHNDXzOKdc7mBxiYiIsVTy59IMZxz24F3gJ4FumnPwjMRZKZzLqeE++cWTPy855xz7l/AEqAjntY2nybe46KCiV+B+xY75w4WUf5/Rbz3LufcSv/f4VG6eI9fB1LYzM7C05rZEE/y+UARxW7wHh/yJX7e2PLxTOT43t97eOtuG1AbaBNIXCIiUjIlfyL+TfMexxQ6TiMAZtbMzK4xs0fM7Dkzm+Yd43act0j7AsVX4ZlEcqWZXWtmTQo/r5Dl3uNMM+tT3LI0AWrsPe4qqaCZjcHTPZ0HnOWce7mIMvF4WioBZha+7pw7jKdLuySZ3mNJdSEiIgFSt6+If2/iSYhGm9mDeMafrQ2kVc3M7sTTwuXv/7Mk3z+cc5vM7Gbgn8CTwJNm9h2eVsI3gNcLzcKdhKfFboj3dcDMVuIZ9/eSc+67wL9N6niPe/2WgmZ4xjo6YLBz7tNiyjUEEoF8im/hSw8gLl88df2WEhGRgKnlT8QPb9fuTDzj2l7Ak9CUNNEDMxsB3A5k4Rn31gao4Zwz55wBvjF6Vuj9pgAt8Sy2/P/t3U+IlVUYx/Hvj9EWFYEoEVjhkBUkGVrQtLBFBYJWoOUmhhaRq8iViwyjgoiUAqOIiIiCJgKxPxQRFJbZLihKigIrigoxJRRaTMF9Wjzn1Zd35t57bjkU3d9nc5h3zvvec2cxPDznfZ4zA0wA08BeskikHSweAa4GbgQeIzOH1wIPAd9IumuEr9psT583cFYWlrxb1r2nWwndR/S53qu4t1nPbxVzzcysgoM/s+FeLOPNZOXuTMU9W8p4f0Q8HxHfdapgV/a7MSKORMSzETEdESvI7N6hMt7XmduLiP0RsSMirgeWljmLyMzhsGCucbSMS4fM+4PMfr4BTAH7JS2bZ95xYJb8H3NRn2etqFhXs56jA2eZmVk1B39mQ5RihY/JgGZvRNQEIs2xZHO2PEvxyJoRPv9z4Mny41VD5v4eEbuAn8g2NZdXfkxTkHHFwFmcyoZuIbOXa4APJV3QmfMnWa0Lc1vlIGkxcNugz5G0hHw38iRweNi6zMysjoM/swoRsS4ilkVEbcPhr8u4VdJZzUVJ5wMvMc97gJJukLShFEu0r08AG8qPP7Suby999rrPuYbcpu6RQWCND8p4Xc3k0nZlGngBWAUckHRhZ9pTZdwuqakmblq9PAJcPORjpsjt5YP/gRNHzMz+Nxz8mS2MPcAJYCNwuJzC8TbwLXAuuW3atZqsoj1WTsyYkfQ6mT3cTB53tqs1fyfwo6SvJO0rp3ocJBsoTwC7I6LvCSRtEfE98AVwiaTJynt6wN3A02TV8kfteyNiHxkcLgE+kfSepFfIwHgbpXciuZU8n5vK+GbNeszMrI6DP7MFUCpt1wKvktmrW8i+fs+R2bUT89z2FvAwuQW7ktwWXUcGfQ8CqyOiXSF7D5lF7JEnbGwClpfnrI+IHSMu+5ky3ll7Q+k/eC+wG5gkA8B2+5qtZPHKl+W7rCd7CU4Bv5Q5x7rPLdnOO8i/05xWMWZm9vfp9OECZjbOJJ1DbiufBC5d6K1WSe+Tlcq3lyxh+3e3khm/JyJi2EkgZmY2Amf+zAzIYhEywzjJCNm/QSStknR259piSTvJwO9X8oi4rgfIBs+Pnol1mJnZac78mdkpZbv1M7Lp82URMfsPn/cyuR39KfAz2az5SvKIvFlgc0S807lnE/AasK30PTQzszPIwZ+ZLRhJG8n3/taSPfsWkef1HgAej4hD/+LyzMzGkoM/MzMzszHid/7MzMzMxoiDPzMzM7Mx4uDPzMzMbIw4+DMzMzMbIw7+zMzMzMbIX07FX2l31XeCAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"#Plot all 3 solutions\n",
"plt.plot(m,num_sol_heun[:,1], 'b', label = 'Heun');\n",
"plt.plot(m, num_sol_euler[:,1], 'r--', label = 'Euler');\n",
"plt.plot(m,-np.log(m/0.25)*250, 'o', label = 'Tsiolkovsky')\n",
"plt.xlabel('Mass (kg)')\n",
"plt.ylabel('Velocity (m/s)')\n",
"plt.grid(True);\n",
"plt.legend(loc='center left',bbox_to_anchor=(1,0.5));\n",
"\n",
"print('According to the graph, the solutions converge')"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0.00000000e+00 0.00000000e+00 3.20000000e-03 9.60512821e-03\n",
" 1.92205293e-02 3.20513644e-02 4.81028115e-02 6.73800651e-02\n",
" 8.98883364e-02 1.15632854e-01 1.44618862e-01 1.76851624e-01\n",
" 2.12336418e-01 2.51078541e-01 2.93083307e-01 3.38356046e-01\n",
" 3.86902108e-01 4.38726859e-01 4.93835682e-01 5.52233979e-01\n",
" 6.13927168e-01 6.78920688e-01 7.47219992e-01 8.18830555e-01\n",
" 8.93757868e-01 9.72007439e-01 1.05358480e+00 1.13849549e+00\n",
" 1.22674508e+00 1.31833915e+00 1.41328331e+00 1.51158317e+00\n",
" 1.61324437e+00 1.71827258e+00 1.82667347e+00 1.93845274e+00\n",
" 2.05361610e+00 2.17216929e+00 2.29411807e+00 2.41946821e+00\n",
" 2.54822551e+00 2.68039577e+00 2.81598484e+00 2.95499856e+00\n",
" 3.09744282e+00 3.24332350e+00 3.39264652e+00 3.54541782e+00\n",
" 3.70164335e+00 3.86132909e+00 4.02448103e+00 4.19110520e+00\n",
" 4.36120762e+00 4.53479437e+00 4.71187152e+00 4.89244517e+00\n",
" 5.07652145e+00 5.26410650e+00 5.45520649e+00 5.64982760e+00\n",
" 5.84797605e+00 6.04965807e+00 6.25487992e+00 6.46364786e+00\n",
" 6.67596820e+00 6.89184726e+00 7.11129139e+00 7.33430694e+00\n",
" 7.56090031e+00 7.79107790e+00 8.02484617e+00 8.26221155e+00\n",
" 8.50318054e+00 8.74775964e+00 8.99595537e+00 9.24777429e+00\n",
" 9.50322298e+00 9.76230802e+00 1.00250361e+01 1.02914137e+01\n",
" 1.05614477e+01 1.08351447e+01 1.11125114e+01 1.13935546e+01\n",
" 1.16782810e+01 1.19666975e+01 1.22588108e+01 1.25546278e+01\n",
" 1.28541554e+01 1.31574004e+01 1.34643699e+01 1.37750707e+01\n",
" 1.40895098e+01 1.44076943e+01 1.47296310e+01 1.50553272e+01\n",
" 1.53847899e+01 1.57180261e+01 1.60550431e+01 1.63958479e+01\n",
" 1.67404478e+01 1.70888500e+01 1.74410618e+01 1.77970903e+01\n",
" 1.81569429e+01 1.85206269e+01 1.88881497e+01 1.92595186e+01\n",
" 1.96347411e+01 2.00138246e+01 2.03967766e+01 2.07836046e+01\n",
" 2.11743160e+01 2.15689185e+01 2.19674196e+01 2.23698270e+01\n",
" 2.27761483e+01 2.31863912e+01 2.36005633e+01 2.40186724e+01\n",
" 2.44407263e+01 2.48667328e+01 2.52966996e+01 2.57306348e+01\n",
" 2.61685460e+01 2.66104413e+01 2.70563287e+01 2.75062160e+01\n",
" 2.79601114e+01 2.84180228e+01 2.88799584e+01 2.93459262e+01\n",
" 2.98159344e+01 3.02899912e+01 3.07681049e+01 3.12502835e+01\n",
" 3.17365355e+01 3.22268691e+01 3.27212927e+01 3.32198146e+01\n",
" 3.37224433e+01 3.42291873e+01 3.47400549e+01 3.52550548e+01\n",
" 3.57741955e+01 3.62974855e+01 3.68249336e+01 3.73565483e+01\n",
" 3.78923384e+01 3.84323126e+01 3.89764796e+01 3.95248484e+01\n",
" 4.00774276e+01 4.06342263e+01 4.11952533e+01 4.17605176e+01\n",
" 4.23300282e+01 4.29037941e+01 4.34818244e+01 4.40641282e+01\n",
" 4.46507146e+01 4.52415929e+01 4.58367723e+01 4.64362620e+01\n",
" 4.70400714e+01 4.76482097e+01 4.82606865e+01 4.88775111e+01\n",
" 4.94986930e+01 5.01242417e+01 5.07541668e+01 5.13884778e+01\n",
" 5.20271845e+01 5.26702964e+01 5.33178233e+01 5.39697750e+01\n",
" 5.46261614e+01 5.52869921e+01 5.59522772e+01 5.66220266e+01\n",
" 5.72962502e+01 5.79749582e+01 5.86581605e+01 5.93458674e+01\n",
" 6.00380889e+01 6.07348353e+01 6.14361169e+01 6.21419439e+01\n",
" 6.28523267e+01 6.35672757e+01 6.42868014e+01 6.50109142e+01\n",
" 6.57396248e+01 6.64729436e+01 6.72108814e+01 6.79534488e+01\n",
" 6.87006566e+01 6.94525155e+01 7.02090364e+01 7.09702303e+01\n",
" 7.17361080e+01 7.25066805e+01 7.32819589e+01 7.40619543e+01\n",
" 7.48466778e+01 7.56361406e+01 7.64303541e+01 7.72293294e+01\n",
" 7.80330780e+01 7.88416113e+01 7.96549408e+01 8.04730780e+01\n",
" 8.12960344e+01 8.21238218e+01 8.29564517e+01 8.37939361e+01\n",
" 8.46362866e+01 8.54835152e+01 8.63356337e+01 8.71926542e+01\n",
" 8.80545887e+01 8.89214493e+01 8.97932482e+01 9.06699976e+01\n",
" 9.15517098e+01 9.24383971e+01 9.33300719e+01 9.42267467e+01\n",
" 9.51284341e+01 9.60351465e+01 9.69468968e+01 9.78636976e+01\n",
" 9.87855616e+01 9.97125018e+01 1.00644531e+02 1.01581662e+02\n",
" 1.02523909e+02 1.03471283e+02 1.04423799e+02 1.05381470e+02\n",
" 1.06344309e+02 1.07312328e+02 1.08285543e+02 1.09263966e+02\n",
" 1.10247611e+02 1.11236492e+02 1.12230622e+02 1.13230015e+02\n",
" 1.14234685e+02 1.15244646e+02 1.16259912e+02 1.17280498e+02\n",
" 1.18306416e+02 1.19337683e+02 1.20374311e+02 1.21416315e+02\n",
" 1.22463711e+02 1.23516512e+02 1.24574732e+02 1.25638388e+02\n",
" 1.26707493e+02 1.27782063e+02 1.28862112e+02 1.29947656e+02\n",
" 1.31038709e+02 1.32135287e+02 1.33237405e+02 1.34345079e+02\n",
" 1.35458324e+02 1.36577156e+02 1.37701589e+02 1.38831641e+02\n",
" 1.39967327e+02 1.41108662e+02 1.42255663e+02 1.43408346e+02\n",
" 1.44566726e+02 1.45730822e+02 1.46900647e+02 1.48076220e+02\n",
" 1.49257557e+02 1.50444674e+02 1.51637588e+02 1.52836315e+02\n",
" 1.54040874e+02 1.55251281e+02 1.56467553e+02 1.57689707e+02\n",
" 1.58917761e+02 1.60151732e+02 1.61391638e+02 1.62637496e+02\n",
" 1.63889325e+02 1.65147141e+02 1.66410964e+02 1.67680810e+02\n",
" 1.68956699e+02 1.70238648e+02 1.71526677e+02 1.72820803e+02\n",
" 1.74121045e+02 1.75427422e+02 1.76739953e+02 1.78058657e+02\n",
" 1.79383553e+02 1.80714660e+02 1.82051997e+02 1.83395585e+02\n",
" 1.84745442e+02 1.86101588e+02 1.87464044e+02 1.88832828e+02\n",
" 1.90207962e+02 1.91589466e+02 1.92977359e+02 1.94371662e+02\n",
" 1.95772396e+02 1.97179582e+02 1.98593240e+02 2.00013392e+02\n",
" 2.01440059e+02 2.02873261e+02 2.04313021e+02 2.05759360e+02\n",
" 2.07212299e+02 2.08671861e+02 2.10138067e+02 2.11610941e+02\n",
" 2.13090503e+02 2.14576776e+02 2.16069784e+02 2.17569548e+02\n",
" 2.19076092e+02 2.20589438e+02 2.22109611e+02 2.23636633e+02\n",
" 2.25170528e+02 2.26711319e+02 2.28259030e+02 2.29813687e+02\n",
" 2.31375311e+02 2.32943929e+02 2.34519565e+02 2.36102242e+02\n",
" 2.37691987e+02 2.39288824e+02 2.40892778e+02 2.42503875e+02\n",
" 2.44122141e+02 2.45747601e+02 2.47380281e+02 2.49020208e+02\n",
" 2.50667407e+02 2.52321905e+02 2.53983730e+02 2.55652907e+02\n",
" 2.57329465e+02 2.59013430e+02 2.60704830e+02 2.62403692e+02\n",
" 2.64110046e+02 2.65823918e+02 2.67545337e+02 2.69274332e+02\n",
" 2.71010931e+02 2.72755164e+02 2.74507060e+02 2.76266648e+02\n",
" 2.78033959e+02 2.79809021e+02 2.81591865e+02 2.83382522e+02\n",
" 2.85181022e+02 2.86987396e+02 2.88801675e+02 2.90623890e+02\n",
" 2.92454074e+02 2.94292258e+02 2.96138474e+02 2.97992754e+02\n",
" 2.99855131e+02 3.01725639e+02 3.03604310e+02 3.05491178e+02\n",
" 3.07386276e+02 3.09289638e+02 3.11201299e+02 3.13121294e+02\n",
" 3.15049657e+02 3.16986423e+02 3.18931628e+02 3.20885307e+02\n",
" 3.22847498e+02 3.24818235e+02 3.26797556e+02 3.28785497e+02\n",
" 3.30782097e+02 3.32787392e+02 3.34801421e+02 3.36824222e+02\n",
" 3.38855833e+02 3.40896294e+02 3.42945644e+02 3.45003922e+02\n",
" 3.47071169e+02 3.49147425e+02 3.51232731e+02 3.53327127e+02\n",
" 3.55430657e+02 3.57543360e+02 3.59665280e+02 3.61796459e+02\n",
" 3.63936941e+02 3.66086768e+02 3.68245985e+02 3.70414636e+02\n",
" 3.72592766e+02 3.74780419e+02 3.76977642e+02 3.79184481e+02\n",
" 3.81400981e+02 3.83627190e+02 3.85863154e+02 3.88108923e+02\n",
" 3.90364544e+02 3.92630066e+02 3.94905539e+02 3.97191011e+02\n",
" 3.99486534e+02 4.01792157e+02 4.04107933e+02 4.06433913e+02\n",
" 4.08770149e+02 4.11116695e+02 4.13473603e+02 4.15840928e+02\n",
" 4.18218724e+02 4.20607047e+02 4.23005951e+02 4.25415494e+02\n",
" 4.27835732e+02 4.30266723e+02 4.32708524e+02 4.35161195e+02\n",
" 4.37624795e+02 4.40099385e+02 4.42585023e+02 4.45081773e+02\n",
" 4.47589696e+02 4.50108855e+02 4.52639314e+02 4.55181136e+02\n",
" 4.57734386e+02 4.60299131e+02 4.62875437e+02 4.65463370e+02\n",
" 4.68063000e+02 4.70674394e+02 4.73297622e+02 4.75932756e+02\n",
" 4.78579865e+02 4.81239022e+02 4.83910301e+02 4.86593775e+02\n",
" 4.89289518e+02 4.91997608e+02 4.94718120e+02 4.97451131e+02\n",
" 5.00196722e+02 5.02954970e+02 5.05725958e+02 5.08509766e+02\n",
" 5.11306477e+02 5.14116175e+02 5.16938945e+02 5.19774873e+02\n",
" 5.22624046e+02 5.25486553e+02 5.28362482e+02 5.31251924e+02\n",
" 5.34154973e+02 5.37071719e+02 5.40002259e+02 5.42946688e+02\n",
" 5.45905103e+02 5.48877602e+02 5.51864286e+02 5.54865255e+02\n",
" 5.57880613e+02 5.60910464e+02 5.63954913e+02 5.67014068e+02\n",
" 5.70088038e+02 5.73176933e+02 5.76280866e+02 5.79399950e+02\n",
" 5.82534302e+02 5.85684038e+02 5.88849278e+02 5.92030143e+02]\n"
]
}
],
"source": [
"#Calculate heights\n",
"height = num_sol_euler[:,0]\n",
"print(height)"
]
},
{
"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": 34,
"metadata": {},
"outputs": [],
"source": [
"def rocket(state,dmdt=0.05, u=250,c=0.18e-3):\n",
" '''Computes the right-hand side of the differential equation\n",
" for the acceleration of a rocket, with drag, in SI units.\n",
" \n",
" Arguments\n",
" ---------- \n",
" state : array of three dependent variables [y v m]^T\n",
" dmdt : mass rate change of rocket in kilograms/s default set to 0.05 kg/s\n",
" u : speed of propellent expelled (default is 250 m/s)\n",
" c : drag constant for a rocket set to 0.18e-3 kg/m\n",
" Returns\n",
" -------\n",
" derivs: array of three derivatives [v (u/m*dmdt-g-c/mv^2) -dmdt]^T\n",
" '''\n",
" dstate = np.zeros(np.shape(state))\n",
" derivs = np.array([state[1],(u/state[2])*(dmdt)-9.81-((c/state[2])*((state[1])**2)),-dmdt])\n",
" \n",
" return derivs\n",
" return dstate"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [],
"source": [
"#Set initial conditions\n",
"y_0 = 0\n",
"v_0 = 0\n",
"m_0 = .25\n",
"\n",
"m_f = .05\n",
"dmdt = .05\n",
"N = 500\n",
"dt = 4/N"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [],
"source": [
"#Initialize array\n",
"num_sol_heun = np.zeros([N,3])\n",
"num_sol_euler = np.zeros([N,3])"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [],
"source": [
"#IC's\n",
"num_sol_heun[0,0] = y_0\n",
"num_sol_heun[0,1] = v_0\n",
"num_sol_heun[0,2] = m_0"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [],
"source": [
"num_sol_euler[0,0] = y_0\n",
"num_sol_euler[0,1] = v_0\n",
"num_sol_euler[0,2] = m_0"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [],
"source": [
"for i in range(N-1):\n",
" num_sol_heun[i+1] = heun_step(num_sol_heun[i],rocket,dt)\n",
"for i in range(N-1):\n",
" num_sol_euler[i+1] = euler_step(num_sol_euler[i],rocket,dt)\n",
"\n",
"m=num_sol_euler[:,2]"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"According to the graph, the solutions converge\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAn8AAAEaCAYAAABgnzGPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd1zV1RvA8c9howKKigMVHIALFXGhpqhYWpqklQ1H5arMUitXu8yGP21oWZpmqeXIlVqZC3ORiol7K7hwoYKCzPP744ICMu69XGT4vF+v+7pyzznf+3gzeDjne56jtNYIIYQQQoj7g1VhByCEEEIIIe4dSf6EEEIIIe4jkvwJIYQQQtxHJPkTQgghhLiPSPInhBBCCHEfsSnsAIq7ChUqaE9PT7PG3rx5k9KlS1s2oBJMPi/TyOdlOvnMTJOfzyssLOyy1rqihUMSQhhBkr988vT0ZOfOnWaNDQkJITAw0LIBlWDyeZlGPi/TyWdmmvx8XkqpCMtGI4Qwliz7CiGEEELcR4ps8qeUmqCU0mmPN3Lp94xSapNS6rpS6oZSaqdSaqhSKte/m7njhBBCCCGKsyKZ6CilmgOjgFyPH1FKfQPMA5oBm4A1gDcwFfhNKWVtyXFCCCGEEMVdkUv+lFL2wGzgArA8l369gJeBKKCR1rqb1voxwAs4CDwGvGKpcUIIIYQQJUGRS/6AD4H6wIvA9Vz6jU17Hq21Ppr+otb6AvBS2pdjslnGNXecEEIIIUSxV6QSHKVUS+B14Bet9Ypc+lUD/IFEYFHWdq31RuAsUBlold9xlrbsv7O0+XQ9z/11kzafrmfZf2cL6q2EEEIIITIpMsmfUsoB+AmIBl7Lo7tf2vN+rXV8Dn12ZOmbn3EWs+y/s4xdspez1wxvf/ZaPCMW7ObtZXsL4u2EEEIIITIpSnX+PgZ8gKe01pfz6Fsz7Tm3OlGRWfrmZ5zFTFx9mPiklEyvaWBuqOFtxwf7FsTbCiGEEEIARST5U0q1BoYDy7TWC4wYUibt+WYufW6kPTtZYFwmSqnBwGCASpUqERISksvlMkuf8cvO3NBISsVF0bqqrdHXu5/cuHHDpM/6fiefl+nkMzONfF5CFE+FnvwppRyBH4EYDLtwjRqW9pxrKRgLjstEaz0dmA7QrFkzbUqFe/fQ9bkmgPMOpTDumc75Ca/EktMXTCOfl+nkMzONfF5CFE9F4Z6/CRhq7I3UWp83ckxs2nOZXPqkt8VmeM3ccRbz5kM+tzPQ7NxMTJENIEIIIYQoMEUh+XsMSAX6K6VCMj6ALml9Xkp77Ye0r0+lPXvkct3qWfrmZ5zFBPu582yrGrn2eWupbP4QQgghRMEoCskfGOJon82jUlp7rbSvm6V9/V/ac4O0ZePsNM/SNz/jLGp8sC+Otjl/9DcTU2T3rxBCCCEKRKEnf1prT621yu6BofQLwJtprzVJG3Ma2AXYAU9kvaZSqj1QDcMpHtsyvJdZ4wrCJz0b5do+NzRSln+FEEIIYXGFnvzlwydpz58ppeqkv6iUcgO+TfvyU611qoXGWVSwnzt2ud38hyz/CiGEEMLyim3yp7X+DZiG4TSOvUqpFUqpJcBRDMfDLQOmWmpcQXjO1y7Xdln+FUIIIYSlFXqpl/zQWr+slNoMDMVwT6A1cAiYBUzLafbO3HGW1rqqLXGlKt8u8JwdKf4shBBCCEsq0smf1vo54Lk8+vwC/GLGtc0aZ2njg31ZHHaG+KSc8825oZE083Al2M/9HkYmhBBCiJKo2C77liR5bf4Auf9PCCGEEJYhyV8REOznTp88av/J/X9CCCGEsARJ/oqIvGr/gZR/EUIIIUT+SfJXhMjyrxBCCCEKmiR/RYixy78y+yeEEEIIc0nyV8SMD/bNMwEcsWC3JIBCCCGEMIskf0VQXvf/aWDkQkkAhRBCCGE6Sf6KqLzu/0vVcv+fEEIIIUwnyV8RFeznnufuX7n/TwghhBCmkuSvCJPdv0IIIYSwNEn+ijBjd/8+O2PbPYpICCGEEMWdJH9FnDG7f7ccj5YEUAghhBBGkeSvGDDm9I8tx6Pl/j8hhBBC5EmSv2JC7v8TQgghhCVI8ldMyP1/QgghhLAESf6KkfHBvrSp7ZprH7n/TwghhBC5MSn5U0qVVUoFK6U+UEp9p5Sar5SalvZ1D6VU2YIKVBjMGxRAHrf/seV4NG8vkyVgIYQQQtzNJq8OSilroCfwMvAAoNKbMnTT6c9KqX+Ab4GlWusUC8Yq0kx8ognDF+zOtc/c0EjAMFsohBBCCJEu1+RPKfU08ClQDUOydxkIBQ4A0UAM4AyUB+oDrYBAoD1wWik1Rms9v6CCv18F+7mzMyL6doKXk7mhkTTzcCXYz/0eRSaEEEKIoi7H5E8ptQVDMncJ+Ar4SWsdntcFlVJNgOeAp4F5SqlhWus2lglXpBsf7MvJSzfYcjw6135vLd0ryZ8QQgghbsvt7rHawEightZ6pDGJH4DWerfWejhQHXg97TqiAMwbFJDnBpCbiSly/58QQgghbss1+dNaf6W1TjTnwlrrRK31l0At80ITxpg3KCDPAtBzQyMlARRCCCEEkEvyp7W+aYk30FrHWeI6ImfGFICeGxopJ4AIIYQQQur8lQTGFIAGGLFgtySAQgghxH3O6ORPKeWqlGqqlHLN8noVpdRspdR/SqmlSqnGlg9T5GV8sG+eCaAGRi6UBFAIIYS4n5ky8zcW2IFhIwcASik7YDPQF2gM9AA2KKVke2khGB/sm+f9f6lazgAWQggh7memJH8dgJNZdv32BmoCG4EuwDdAWeAVi0UoTGLM/X+yA1gIIYS4f5mS/FUDjmV5rRuG1cSBWuu/tdbDgJNAVwvFJ0xk7P1/sgNYCCGEuD+ZkvyVw3DCR0YBwBGt9YkMr/1HhqVhce8Zc/8fSAIohBBC3I9MSf7iMRzjBoBSqjqG2cAtWfolAPb5D03khySAQgghhMiOKcnfIaBtht2+z2BY8v0nS79qwAULxCbyyZgNICA1AIUQQoj7iSnJ3xygNLBdKbUQ+BC4ASxP76CUsgeaAoctGaQw3yc9G2Gl8u43evGegg9GCCGEEIXOlORvGvALhuPaHgcSgUFa6+sZ+nTHkCButFiEIl+C/dyZ/GQT8sr/EpJTeXbGtnsSkxBCCCEKj9HJn9Y6VWvdB6gNtAbctdYLs3Q7ATwB/GS5EEV+Bfu580XvJnnOAG45Hi0JoBBCCFHC5Zj8KaWeUko5ZX1da31Sax2qtY7Jpm2X1nqx1jrK0oGK/DF2BlASQCGEEKJky23m7xfgolJqpVJqgFKq4r0KShSM9BnAvGw5Hi07gIUQQogSKrfkbzSwG0PB5unAOaXUBqXUq0qpvGuIiCIp2M+dNrVd8+wnJWCEEEKIkinH5E9rPVFrHYChdMurGEq6tAG+BE4qpXYqpcYqperem1CFpcwbFCAJoBBCCHGfynPDh9b6vNb6G611J6AS8AKwCqgPfAzsV0odVEqNV0o1K9hwhaWYkgBKDUAhhBCi5LAxpbPW+iowG5itlCoNPAL0xLA0PA4Yq5Q6AywBlgH/aK21RSMWFjNvUABe41aRlJp7vxELdgOGJWMhhCgMYWFh5aysrPpbW1s/q7WuAHnuXxPivqOUuqG13pScnLwI2ODv759tDmZS8peR1vomsBBYqJSyAzpjSAS7A69hWCp+Dxhv7nuIgjfxiSaMXLib1FxSdA2MXCgJoBCicISFhdnZ2NjMdnV19a1QocINBweHK0pJ7idERlprkpOTrWNiYoIuXbrUOT4+/oewsLDPsksAzU7+srxhIoal4FVKKSugPfAYcNES1xcFJz2ZG7FgN7lN0aZqmQEUQhSap52dnX2rVasWLUmfENlTSmFra5tSvnz56y4uLlbHjh0bePPmze3A+qx9TTnhwyhpxaA3aK1f1VpPt/T1heUZWwJGA8MX7JZNIEKIe8rW1vZhV1fXREn8hDCOjY1NasWKFbGxsXki23ZzLqqUqgxUBRxy6qO13mrOtUXhCPZzZ2dENHNDI/PsOzc0kmYerjIDKIS4VxqVLl06rrCDEKI4cXZ2vqGUeiC7NpOSP6XUE8CHgHceXbWp1xaFb3ywL4BRCaAsAQsh7hWttaO1tfXNwo5DiOLExsYmWWtdNts2Yy+ilHoKmIdhh9V14BRwwxIBiqLD2ARQNoEIIe4lWfIVwjS5/T9jyuzcuLTn14BpWuvk/AQlii5jE8BUDaMX75HkTwghhChGTNnw4QVs1VpPkcSv5Bsf7EufVnmf4peQnIrXuFVSCFoIIYQoJkxJ/qKB0wUViCh6jE0Ak1INS8CSAAohhBBFnynJ399A84IKRBRNxiaA6UvAQggh7j13d3dfpZT/ypUrnXLr16JFCx+llP/XX39d/l7FJooeU5K/9wBnpdRnSinrggpIFD2yBCyEEEKUHEZv+NBaRyql2gLLgZ5KqXXAGSDbk2G11hMsE6IoCsYH+3Ly0g22HI/OtV9SqqEQ9M6I6NsbR4QQQghRdJhS6kUBr2DY+GEN1IZsTwRTaa9L8lfCzBsUwLMztuWZAIIUghZCCCGKKlOWfccAwzAkdquArzAkeFkfHyOJX4k1b1CAUUvAYCgELUvAQghRPKxfv750t27dalWqVKmRra1t03LlyjXu2LFjndWrV5fJ2vfw4cN2Sil/d3f3HJd4lFL+Sin/3F6fMWNGuSZNmtQtVaqUX+nSpf0CAgK8s3s/YVmm1PkbAMQBbbXWuwsoHlEMGLsEnH4WsCwBCyFE0fbee+9V+uijj6oB1K9fP65p06Y3zp8/b7dx40aXjRs3unz++ecRr7/++mVLvufw4cOrTpkypUrTpk1vdOjQ4frBgwcdQ0NDnbp37+79xx9/HA4KCpJTXQqIKcmfO7BBEj8Bpi8Bn7x0g3mDAu5BZEKI+41S3DW7VFxoTVhhx/Dbb785f/jhh9UqVqyYNH/+/OMdO3a8nXT9/fffpXv16uU1ZsyYGp07d45t1KhRgqXed/bs2W4hISEHH3jggTiAlJQU+vTp4zF//vwK7777btWgoKCjlnovkZkpy75nMcz8CQGYtgS85Xg0z87YVsARCSHE/a179+7e6cuq2T127Nhx15Lqhx9+WBVg6tSppzImfgAPPvjgzREjRpxPTk5WU6ZMqWjJWEeNGnU2PfEDsLa25n//+99ZgLCwMKeEhAQ506+AmDLztwAYrJQqrbWWqVgBGH8UHBgSwLeX7ZUlYCGEKCBt27aNcXNzS8qpfePGjS5Xrly5/bP//PnzNvv27StdpkyZlJ49e8ZkN6ZTp06xH330ETt37rTovXi9evW6nvU1d3f3ZGdn55SYmBjrCxcuWNeoUUNOFCsApiR/HwGBwO9KqcFa6+MFE5IobkxJAGUJWAhhaUVh6bSoGD16dFS3bt1ic2pv0aKFz5UrV24ncUeOHLHTWnPjxg1rW1vbXJfPo6OjTckZ8lSnTp3E7F4vU6ZMSkxMjHV8fLwpq5PCBKb8h/wdSAY6AAeVUifIuc6f1lo/ZIH4RDFh6gxg58khrBkZWMBRCSGEyE1ycrICQ8L14IMPXsutb/ny5Y2ehUtJScmzj7W1nBdRWExJ/oKyjPNOe2Qnu/p/ooQbH+xLMw9XRi/eQ0JytrW/bzt68SZe41Yx8YkmUgtQCCEKSa1atRIBbGxs9OLFi08ZO87e3l4DxMXFZTs7d/ToUTuLBCgKhCnJX+cCi0KUGMF+7gT7ueM1bhVJued/JKXCyIW7b48TQghxb9WsWTPJy8sr/ujRo44rV650ym3JOKMqVaok29ra6mvXrtmcO3fOpmrVqplmBZcuXepSMBELSzB6PV1rvc6UR0EGLYq+iU80wcqIfVqpGkYv3lPwAQkhhMjWu+++ew5gwIABNZcsWeKctf3WrVtq3rx5LmvXri2d/pq9vb1u1qzZDYA333yzamrqnd/2V69eXeazzz6T3+iLMLmZUhSIYD93Jj/ZBHubvP+JJSSn4jVulZwGIoQQhaBPnz7X3nvvvTNXrlyx7dWrl5enp2fDjh071unSpUutRo0a1a1YsWLjPn361Nm1a1epjOM++OCDs7a2tvqXX36p6OXl1aBr1661fH196z388MM+/fv3v1hYfx+RN0n+RIEJ9nPn8PiueLmVzrNvUqrhNJC3l+29B5EJIYTI6P3337+wadOmA08++eTl1NRUtm7d6rxp0yaXmJgYmxYtWsROmjQpon///pmq+nfu3PnmihUrjgQEBMRGRUXZhYSEuABMnTr15FdffXWucP4mwhg53vOnlPoHGKO13mruxZVSbYBPtNbtzL2GKP7WjAyU00CEEKIAnT171qjfnLdv3344p7aAgID4gICACFPe96GHHrrx0EMPHcmuTWudbQmenF5PZ+zfRZgvt5m/usAmpdQapVRvpZS9MRdUStkrpZ5WSq0F/iHnHcHiPmLqaSCdJ4cUbEBCCCHEfSq33b5ewAfAy0BHIFYptQXYBhwErgAxgDNQHqgPBABtgDIYagJ+BbxfQLGLYmZ8sC8nL90wagZQSsEIIYQQBSPH5E9rfR0YrpT6GhgG9Ae6Al1yuZ4CooFJwLda61OWC1WUBPMGBRi9BJx+H+DOiGg5Ek4IIYSwkDzr/GmtTwAjlFLjgPYYjnhrArgBLsA14CKwC9gAbNJaJxRUwKL4mzcogLeX7TXqNBC4c2pIUNmCjEoIIYS4Pxhd5FlrHQ/8lfYQIl9MOQ0EDAlgmKsiMLDgYxNCCCFKMin1IgqNKaVgAA5Ga9kIIoQQQuSTJH+i0K0ZGUib2q5G9U3fCCIFoYUQQgjzSPJXSBITCzuCosWUUjBSEFoIIYQwn9H3/AnL6tEDUk6Ux6rxVOq//yQV6rsVdkiFLn1Hr6kbQWQnsBBCCGE8mfkrBBcuwJo10OLIctotGkbZBlXZ4fYIW1+bT3x0fGGHV6jGB/vyZW/jzgQGQwL47IxtBRyVEEIIUXJI8lcIQkMBrXmWeQDYkELzS3/Q+uunSSxfmU0+A9iy8iqpeW+CLZFM3Qiy5Xi03AcohBBCGEmSv0LQowecP5vK3g5PEe7UNlObCzHUPvIn7bo7U6sWvP02HM7xJMaSzZSNIHIfoBBCCGEcSf4KScXK1ri9257GMZuIDDnBhg4fcsq2DgBz6UMq1kREwMcfQ926MMznbzY+MZUrhy8XcuT3likbQcCwDCwJoBBCCJEzo5M/pdRApZRjQQZzv6rRviYd1r+Dx60j7JuxjZg+Q3HNMuH18JEvaP/bMJzrViG0SjD/vvU7SfHJhRPwPTY+2NfkBFCWgYUQQojsmTLzNx04o5SapJTyKqiA7mfKStFwYCvGz/Hg/HlYuhQeewyq2UTxIH8DYEsyraKW03JCD6LLVGdDy9Gc+LPkrwuPD/ZlcCM7ozeCyDKwEOJ+4u7u7quU8s/rsXLlSqf8vlf6tSwRtygcppR6WQl0BUYAryml1gLfACu11roggruf2dlBcLDhER1Rmi1jplJu1Rx8Y+/sbK2UGkWl7Z/Dw58T7tSG671ewG/CkzhVKVOIkRec1lVtGfdMZ56dsY0tx6ONGjM3NJKTl24wb1BAAUcnhBCFr23btjFubm5JObW7u7vn2CbuH6ac7fuoUqo68BLwAvAg0Bk4rZT6Dpiptb5kThBKKVugHfAw0AbwAMoDl4BtwFStdUgu459Ji6sRYA0cAn4Epmmtc9wza+64e83Vw4l2v74EvMSpv49w6v0fqfvvT1ROPX+7T+PYLTB7C5dnj+b1fmfoN8ieNm1AqcKLu6DMGxTA28v2Gl0PcMvxaDpPDmHNyMCCDUwIIQrZ6NGjo7p16xZb2HGIos2kDR9a69Na63FAdaAvEArUAD4GIpVSc5RS5kyxtAfWAiMxJH5hwFIgGugFbFBKfZjdQKXUN8A8oBmwCVgDeANTgd+UUtaWHFfYPB/0JnDrJ1S4GcmOd1fwb5VgkjLk8Kt5kBk/2/PAA4aNIp99BufPlbyJWVPvA5Rj4YQQQggDs3b7aq2TtNbztNZtAD9gJpAMPANsVkqFKaVeUErZG3nJVGAx0E5rXUVr3U1r3Vtr7Qs8BaQA7yilOmQcpJTqBbwMRAGN0sY9BngBB4HHgFeyvpm544oSGwcbmn/QjZbnlnJt7xk2dv8fx+3qMYsXbvc5cgTGjIEl1YaxrUpP/v1wNSlJRWZCM99MLQidfh9gvXf+lCRQCHFfW7lypZNSyr9FixY+2bUfPnzYTinl7+7ubtIRSgkJCerzzz+v6O/v7+Ps7NzE3t6+qYeHR8OBAwdWO3fu3F2rjV9//XV5pZR/r169PKOioqyfe+656u7u7r62trZNg4KCapv79xO5y3epF611OPABhuVSlfbwA2YAp5RSA4y4xnqt9eNa603ZtC0AZqd92SdL89i059Fa66MZxlzAsJwLMEYplfXvae64Iqliw0q0//11asXvZ8K2jgweDE5pt/SW4iZ99c8ERC2l5XtdOOdYm/VBE4jaHVW4QVtIekFoY+sBAsQnpcpmECGEsLDo6Gir1q1be48ePbrGkSNHHBs0aBAXGBh4PTk5Wc2cObNSs2bN6h0+fNguh7E2zZo1q79s2bLy9erViwsKCrqW272LIn/ydbavUioIwwxaNwz3zN0CfsGwhNoHwz1805VSpbXWX+fjrf5Le66W4b2rAf5AIrAo6wCt9Ual1FnAHWgFbM3PuOJAWSlatoKWrWDyZFi8GA5MDMF5353bP6qnnKL6urdI8nuPrVV7YDt0CP6jOmFl5OxZUTVvUADL/jvL6MV7SEg2bnZTzgYWogQZObIqX3xRxai+Tz11mV9/jcj02tNPezB/fgWjxo8YcZ7Jk89leq1jxzps2OBi1PiJEyN4440SV7S1X79+nrt27SrTpUuXqz///HNExYoVUwCSk5MZNmyY+3fffVe5b9++Nbdv335XiYqQkBCXNm3axKxYseJ4uXLlSs4SVRFl8k98pZSLUmq4UuoQsBoIBs4B44BqWuuBWusFWuvuQGvgJvBqPuNMLy1zPsNrfmnP+7XWOR2IuyNL3/yMK1ZKl4Z+/eDTvY9w8o+D/OM/gmh1Z3bMlmRan1tM87ceJNLRm/VdPuPivouFGHH+mXosHMjZwEKIkqV79+7eOZV5cXJyalJQ7xsWFuawatWqclWrVk1ctGjRyfTED8DGxoapU6ee9fb2jt+xY0eZ7du331Uz2MbGRs+cOTNCEr97w+iZP6VUUwyzfE8BjhiWd0OAKcDy7HbHaq3/VUqtwrBpwyxKqcrAc2lfLs7QVDPtOfNvb5mlbwetmeE1c8cVWzW71qVm18kkXJ/AtnG/UXru9zSK2Xy73TP5OJ6rx7Bp9Spe7vUPQ4ZAp05gVUwnA9eMDDSpHMyW49HUGrOKyb2bEOznXsDRCSFEwcmt1Iujo2OBJVa///67C0CnTp2ulylT5q5dhtbW1rRo0eLGkSNHHP/555/SLVq0yDT5Ur9+/TgfH5/EgopPZGbKsu/OtOc44AdgitZ6nxHjbpr4PrcppWyAuYALsE5rvSJDc3oxu5u5XOJG2nPGopbmjssY12BgMEClSpUICQnJ5VK5vMmNG2aPNdsT1Uh44iOWb43C9se/aX1sKWW5BsAPDGDxYsNycdWq8XTrdo4uQWcoV7Fo7BY25fMa5AV2idZsOJ2Sd2cMO46GL9jN92vCGd2ilPlBFiGF8u+rmJPPzDRF6vOaPPncXUuxpvj114i7loJNsX79MbPHWlBhlXo5ceKEPcCcOXMqzpkzp2JufS9dunRXTlCtWjVJ/O4hU5KyUxiKOs/UWl8zYdwgYIgpQWXwHdAJOM3dmz3SK9iZmpmYO+42rfV0DCee0KxZMx0YGGjWdUJCQjB3bL4FAuOe4lb0VLaOXYTNwl9YeO3J283nzjkyY3pNhk0P5mqNRpR6YyhNhwagrAqvcKCpn1dgICz77yxvLd3LzUTjksCD0ZqPd1EiagIW6r+vYko+M9PI51VypaQY9z0za/8GDRrE+fj45HRLFQANGza8lfU1BwcHWe69h0xJ/mqbc5JH2hjT/hUBSqmvgAEYyrF00lpn3Z6a/ptNbsdZpLdl/C3I3HElkoNrKVp/3x++78/O/TB9Ovz8M1y7Bh3YQEP2QeQ+ePUX9o9qysUnhtJ80tOUqVg8jnkO9nMn2M/dpKLQRy/exHPMKvq0qiGbQYQQJZK9vX0qQFxcXLY3+Bw/ftzYUm0AVK9ePRGgTZs2sd9///2Z/EcoCpIpd3WtVkqNzKuTUmqEUurvfMSEUmoShk0ilzAkfkez6XYq7dkjl0tVz9I3P+NKvAYN4Kuv4OxZmD0b+rmvy9x+axcd5gwgsVI11jUbxYl1JwsnUDOYWhMQDJtBOk8OKbighBCikHh4eCQBREZG2ickJNy1pLNy5Urjdi6n6d69ewzAX3/9VTYpSSq0FHWmJH9BQEMj+tXHsFRrFqXU5xhO+rgCdNZaH8iha3r5lwZKqZymoZpn6ZufcfeNUqWgf3/of2YCRxfsYkvdAcTjcLvdVUfTKWwinkG1CXV7lNAP/y4WxaPN2Q2cPgsoNQGFECWJt7d3YvXq1RNiY2Ot33///UoZ2+bMmVP2xx9/dDPlem3bto0LCgq6FhkZaf/II4/UPn78uG3WPhEREbYffvihmySHhS9fdf5yYIfh/nmTKaU+Bd4ErmJI/MJz6qu1Pq2U2gU0BZ4Afs5yrfYY6gJGYTgfOF/j7ldeT/rh9eQPxJz6nJ3DZ+G56huqJ58CwApNq0sr4L0VfPLl59iOe5MXXgBX4+stF4o1IwNNWgYGwyzggu2RTHxCdgQLIYquzz77rPKPP/5YPqf2Z599Nrpnz54xAO+///7ZgQMH1vr000/df//993I1atRIOHnypMORI0ccX3nllfNTpkwxrm5imgULFpzs0qWL15o1a8o2aNDAxcfHJ65atWqJsbGx1ufPn7c7ceKEQ2pqKm+88cYlW1vborGT8D5l0eRPKaUwFFA2uXilUuojYDRwDUPiZ8ys2ycYCjV/ppTaqrU+lnYtN+DbtD6fZlOGxtxx9y1nT1ceWPYGqUkj2PXJnzB1Ck0vGTo8FiIAACAASURBVFb3k7Hm26tPceZNeOcdePZZGDoU/IpwlcTxwb6MD/al8+QQjl7MbeP3HenHw+2MiJZ7AYUQRdLmzZudc2tv3LhxXHry98ILL1y1t7c/NnHixCqHDx92jIiIcKhfv37cokWLjjZs2PCWqcmfq6tr6tatWw9///33rr/++mv5/fv3l9q/f38pZ2fnFDc3t6Rnnnnm0mOPPXatVKlSkvgVMpXbHo4s9+4FYSjmnNMyrA2GYsxVgd+01r2NDkKpR4HlaV/uBPbn0PWQ1vrTLGO/xXAk2y1gLZCEYdnZGVgGPK61vmvDibnjsmrWrJneuXNnXt2yVdx3yp1Zd5hTo74lYl8sfRJnZWrz5CSLXAYRP2AoLT/qjl2p/P+eUVCflyk1AdM52lrxSc9GRXoWsLj/+yoM8pmZJj+fl1IqTGvdzJi+4eHhpxo3blziTsQQoqCFh4dXaNy4sWfW1/P6iRyU4c8aQ2JXNY8xe4BRJkUHGRcKm6U9srMRyJT8aa1fVkptBoYC7TEcM3cImAVMy2n2ztxx4o5qnXyoFvYVTeNg5nyYMgV27za0vcQ0ml1fB5PXcfbL6hwMfImGXw6ksm+u5Z8KRfrRcG8sCic51bhfSNPPB5ZZQCGEEMVNXhs+Oqc9HsRQH291hteyPtoDdbTWTbTWJhXK1FrP1lorIx6BOYz/RWvdRmvtrLUurbX211p/k1cCZ+44kVmpUvDCC7BrF2zZAn16J/Ecs2+3u6eeJmj9OMo1qsbGmv3ZPWMHphcNKljBfu4cm/AwbWqbdsPi3NBI6r3zJ8v+O1tAkQkhhBCWlWvyp7Vel/ZYC2wBNmZ4Letjk9b6xL0JWxRFSkHr1jBnvi16+042PTCOK1Z3zkm3J5H2p36myeAW7CvTkpABc4i7mlCIEd9t3qAAk0vCpM8Cyo5gIYQQxYHRP+G01g9kvd9OiJxUal6DB/75GOdrp/l36M8cLNM8U7tv3HYCZ/XjZvnqjB8cyfHjhRRoNtJLwvRpVcOkcTILKIQQojgwpc6fECazdXKg5dS+1IvdzpE5/7LVqx8J2N1uP68r886M6nh5Qbdu8NdfkFpEFt3HB/uanADKLKAQQoiiLscNH0qpcWl/nKa1vprha6NorSfkKzJR4nj3aYF3nxZcOzKRHcN/oPbf05ia8gqg0BpWrTI8Xqy8jMfbRuH/RR/KVsvtFL6CNz7Yl2YeriadDwxSF1AIIUTRldtu3/EYdvj+hqHocvrXeVFp/ST5E9kq6+1G2z/GkZo4isf+TOXM9/Dnn+mtmhej3qPxb3u4/tto1jV8nmoTXsanu3ehxWvO+cBwpy7g2CV7inxZGCGEEPeP3JK/CRiSuMtZvhbCIqzsbOjaA7r2gGPH4Ntv4eD0zTS+uQcAF2LotO8rePQrtpd7iHNdupPU6gFsHawLJV5zZwGlLIwQQoiiJMfkT2v9dm5fC2FJderA5Mlw883GbBnxNe5Lp+KZeOR2e4urq2nx62oiFk7iaNDLNPryBdzq3vtz5MydBQTDUvDisDMyCyiEEKJQyYYPUaSUruJMm/nD8Ig7yJ7//c2Oqo+Sirrd7pFykqDVb+JUz51V9d8gNJRCqRk4PtiXL3s3obSdabOQ6bOAz86474+NFkIIUUgk+RNFkrK2otHrnWl+djkXtx5nU8AootWdmT5HbnHwIAQEQPPmMHs2xMff2xiD/dzZ/2EXk+sCAmw5Ho3XuFVSFkYIIcQ9Z/RPLKXUS0qpRKXUI7n06ZbWZ6BlwhMCKgfU5IGtn7Fr+a+EDprJkVJNSEXxLS8DEBYGzz8P1avDvB4LOb39/D2NL70uoKmng6RvCJHagEIIIe4lU6YregLRwJ+59Pkzrc/j+QlKiOzYONnRavoLeN/YxcGF+wh8vhb29nfa7a+cpffvz1C5ZQ3+qfY0O77cgjbyrF5LSD8dRJaChRBCFGWmJH91gb25nXurtU4B9gL18xuYEDlSigZP1GfWLDh7Fj77DDw8YAjfY0MKtiTT7ux8mo9oy8FSTVn37ExiouLuSWjpS8GmFocGw1Kw55hVUiBaCCFEgTIl+asIXDCi30XAzbxwhDBN+fIwahQcPw4Pvd2CvWUfyNReP2E3nX4ZSHKVaqz1e5Njf9+b46fN3RAChl3BnSeHWD4oIUSJopTyN/XRq1cvT3Peq2HDhvWUUv7//PNPqfzGff36dSullH+pUqX8sraVK1eusVLK//z587mVoivyduzY4aCU8vfy8mpQ2LFkx5QP9zpQ3Yh+7sAN88IRwjzW1tDyo27wUTdOLQ/n3Fvf0GT/XEph2AXiylWCdv+P1Icmsa3CIySMfo+2w5thU4DfXtLLwiz77yyjF+8hIdn4c+uOXryJ55hV9GlVQ2oDCiGy1bNnzytZX7t48aLt5s2bnR0dHVO7du16NWt7mzZt5OezMCn5+w/ooJSqrbU+nl0HpVRtoDXwjyWCE8Icnj0a49ljOrGRn7FlxI9UX/ENNZIMM35WaAIur+TBN4dxeAq89BIMGAAVKxZcPPmtDfhLaCSTe8sxcUKIzBYvXnwq62srV6502rx5s3O5cuWSs2s315IlS47fvHnTysfHJ8FS1xSFx5Rl39mALbBMKeWVtVEpVQdYBlin9RWiUDnVKEebxSOpHn+U8AmrCKvUFYDDeLOWICIjYexYwy7hgX1usW/+vgKNx9yl4FQMu4KlNIwQorB4e3sn+vn53SpVqpSc9FUCmJL8LQD+ABoA+5VS65VS36Y91gEH0tpWa63nFkCsQphFWVvReOzD+Ef9wbmQI2x8dgblK9z5p5+QAAnzFtHwaV/CnAPZ+MoibsUmFUgs+akNKKVhhLh35oZGuLb4eK1vzTGr/Ft8vNZ3bmjEvT9SqIBNmTKlfPPmzX2cnZ2b2NjYNC1XrlxjHx+f+s8991z1o0eP2mXsm9s9f/Hx8er999+v1LBhw3qlS5f2c3R09PPy8mrw2muvVb1y5YpFzuNMTk6mT58+NZRS/t7e3vWPHz9um7F969atjt27d6/p5ubWyNbWtqmrq2vjjh071vn999+dsl7L29u7vlLKf/ny5Xe1pXvmmWdqKKX8X3/99Srpr8XExFi9+eabVXx8fOo7Ojr62dnZNXVzc2vUtGnTuiNGjKialGTcz42oqChrPz+/ukop/+7du9eMj49XL774YjWllP/gwYOr5TTuu+++c1VK+bdq1Srfh90b/dNHa60xlHuZlvZSIPBi2qND2mvTgMfyG5QQBaVqey8Gz23H6dPw00+GAtEArzAVAP/YjbT/5kmulvVkzQMfcnpHVIHEkV4b8MveTbCxUnkPyCC9NIwkgUIUjLmhEa4frTzgcTE2wU4DF2MT7D5aecCjJCWAgwcPrvbqq696hoeHl65fv35c165dr/r6+sYlJCRY/fTTT247d+50NOY6169ftwoICPD54IMPqp06dcqhVatWMYGBgdevXLli8/XXX1dp2rRpvayJmqliY2OtHnzwwTrz5s2rGBAQEBsaGnqodu3atzOt6dOnl2vfvn29lStXurq6uiZ36dLlqqenZ0JISIhLjx49vN96663KGa/31FNPXQH48ccfK2T3fvHx8WrlypWuSikGDRp0BSApKYl27dp5/+9//6t64cIFu1atWsU+9NBDV2vVqnXrzJkzdl9++WWVuLi4PHOq/fv327dq1are7t27Sw8ZMuTC8uXLTzo6OuqRI0detLa2ZuHChRXi4uKy/aEwY8aMigAvvvjiReM/veyZNPWgtU7UWg8FagB9gbfTHn2BGlrroVpruR9AFHkODtCvH2zfDts3xqM8PEjmzi+oVVLP0Xnze1RqUbA1A4P93Dk24WGTC0TDnSRQSsMIYVlfrzvqnpCcmunnY0JyqtXX646WiBtvr1y5Yj1r1iw3FxeXlL179+4LDQ09smLFipP//PPP0VOnTu0LCwvb37x5c6PqYw0dOrRaeHh4aR8fn/iDBw/uXbdu3fE///zzxMmTJ/e2a9fuemRkpH3fvn09zY317NmzNq1bt/bZsGGDS48ePaI3bNhw1NXV9fbuucOHD9u99tprnsnJyerzzz+POHTo0IEVK1ac3LVr16H58+cftbW11Z988on7X3/9VSZ9zKBBg65YW1vr1atXl7169epdedC8efPKxsbGWjdv3jy2bt26iQBLlixxCQ8PL920adMb586dC9+wYcOxFStWnAwNDT1y/vz5PStXrjzs4OCQ6w+JDRs2lHrggQfqnj592n7ChAmR33333RkrK8Pbe3t7J3bo0OHa9evXrX/44Ye7fiD8+++/jrt27Srj5uaW9Oyzz14z9/NMZ9bxblrrKK31PK31hLTHPK11wUyRCFHAmrdzpMWphVwPj2Bzx3e5aFXpdpsdSbdrBh4q5ce6p38g5uIti8eQXiDa1FlAMGwKkVlAISznUmyCnSmvFzeXLl2yTklJUbVr14738fFJzNretGnTW3Xq1MlzDfPy5cvWixYtqgAwderUiOrVqyent7m4uKTOmjUrws7OTm/bts1527ZtRs0kZrRnzx77Vq1a1d23b1+pV155JWrZsmUn7e3tMyVYX375ZcVbt25ZtWnTJubNN9+8nLHtySefjOndu/dlrTWTJk26/Y29evXqye3atYu5deuW1U8//VQu6/vOmTOnPECfPn1u76aOioqyAWjTpk1s1vsera2teeSRR25kjS2jefPmuTz88MM+8fHxVj/99NPxsWPHXsraZ9iwYRcBZsyYcVe5vK+++soNoG/fvpdsbfM1kQrI2b5C3Fa+kTtt131A+RuR7Hz9V/a5tMnUXi8hnObzR+JdM4mXX4b9+y37/umzgOYUiJZTQoSwnIpO9nclRLm9Xtx4eXklli9fPvm///4rM2zYMPd9+/bZ5z3qbiEhIaUTExOVh4dHQlBQ0M2s7bVr105q06ZNDMDatWtzvL8uO+vWrSvTvn37uufPn7efOHFixJQpU7L97Xbr1q1OAH379r2r7A3A4MGDLwP8+++/md6/X79+lwHmzp2baek3MjLSZsuWLS6lSpVK7d+//+1SOQEBAXFKKWbPnu02efLkCqbUIfz0008r9uvXr46jo2PqqlWrjvTp0yfbmbtHH300tk6dOrf27dtXKuO9ldHR0VbLli1ztbGx0cOGDbuc3VhTmZz8KaV8lFLfKKX2K6WupT32K6WmKqXqWiIoIQqTtaMdzf73FA2vbebkkv/YWn8gcRh+af2ZflyIc2LaNGjYEDp0gOU/XSMpwfgafnkZH+zLqU8fMWspWE4JESL/Xu3kddbexirT/9T2Nlapr3byKhHT69bW1syYMeOks7NzytSpUyv7+vo2rFChQuPOnTvXnjhxYoXr168blRucPn3aDqB69eo53u7l6emZAHD27FmTZk2fe+65WteuXbMZP3585BtvvJFjwnPhwgU7gDp16mQbQ7169RIAYmNjrWNjY2//vZ566qnrZcuWTd61a1eZAwcO3I7thx9+KJ+SkkLXrl2vOjs73/430KJFi/ixY8eejYuLs3r99dc9qlat2tjDw6Nhr169POfOnVs2JSUl2/hOnjxpP3bs2BpKKf76668jHTt2vCtJzmjw4MEXAKZMmXJ79m/atGkV4uPjrR566KFrHh4eFtmNaFLyp5R6DtiNYZNHPcA57VEPeBnYrZTqb4nAhCgKaj7WhNb7Z5AScZYtvSbxV51hmdpDQiDmuWGcLe3Fmq6TuHgo2mLvbe5ZwWBYCn7+r5uyFCyEGfq08oh+p1v9CDcn+0QFuDnZJ77TrX5En1YelvsfvJA99thjMREREXu+++67k08//fQlV1fXpHXr1pUdNWqUR506dRru2rXLIa9rGPaBglI5366S3sdU6QWsv/jiiyrh4eE5zkwaE0N2HBwcdI8ePaK11syYMeP27N/8+fPLAzz//PN3JZwff/xx1NGjR/dOmDAhslu3btEJCQlWS5YsKd+3b9/a/v7+dW/cuHFXEFWqVEls3rz5jZSUFIYNG1Y9r8R6yJAh0U5OTikrV64sd+nSJWuAWbNmVQQYOnRovjd6pDM6+VNKNQdmAHbAUqAbhqSvPvAIsBhDHcAZaX2FKDGcapSjzW8jWXHEhw0b4PHHDaeKuHGBJ1mIZ8oJOv/1BmXqVWNDnYGEz/4PM7/nZZKxNIypSaBG6gMKYa4+rTyit78VtPfkp4+EbX8raG9JSvzSubi4pA4ZMiT6l19+iTxy5MiBo0eP7unUqdO1y5cv2w4bNizPE71q1KiRCBAZGZljchYREWEP4O7ubtKS+VdffXV2+PDh5y9evGjbsWPHutu3b8/2nsHKlSsnAhw9ejTbGA4dOmQP4OTklOLk5JRpNnfgwIFXABYuXFg+NTWVTZs2lTp69Kiju7t7YteuXbM9CaVmzZpJY8eOvbRixYqTUVFRezZu3HjQw8MjITw8vPT48eMrZe3v4OCg169ff6Rt27Yx27dvdwoMDPTOrfyNs7Nzau/evS/funXL6ptvvqnw+++/O504ccLBy8srPqeYzGHKzN+baf37aK0f11r/obU+rLU+pLX+U2v9BNAHw6khb1gqQCGKEqUgMBAWLYJTp+CzPvuIU6Vvt5cing7HZ9L4+aaEO7Vhw+BfibuW/9uE0pNAc+4HTK8PKMvBQojc1K5dO+mjjz46B3Do0KE8z/ANDAy8aWdnpyMiIuzXr19fOmv7qVOnbLdu3eoMEBQUFGtqPF988cW5cePGnY2OjrZ56KGHvDdt2nRXTK1bt44FmDdvXvnsrjFjxozyAC1btrzr/du2bRvn7e0df+7cObs//vjDaebMmeUBevfufTl9F25e2rVrFzdgwICLAHv37s32MytTpoxes2bNsaCgoGu7d+8u3a5dO++oqKgcE8ARI0ZctLKyYvbs2RW/+eYbN4CBAwfetUEkP0xJ/toCYVrrX3PqkNa2A2iX38CEKOqqVYPn5nSidPQZ/h08kyOlm2Rqb3JzKx1mPEOsaw3WtHqHEyGmHe2WHXNPCUk3NzRSkkAh7nN79uyxnzJlSvnsliCXLl1aFqBq1ap5/tZaoUKFlMcff/wywLBhw2qcO3fu9iaImJgYqxdeeMEjISFBBQQExAQEBMSbE+vHH38c9fHHH0dev37d5uGHH/Zeu3ZtpiRz+PDhlxwcHFI3bdrk/MUXX2TavLFkyRLn+fPnV1RK8frrr1/I7vrPPPPMZYDp06dXWL58eabafhn99ttvzkuXLnVOTk7O9HpCQoJas2aNC0D16tVz/MwcHBz0H3/8cbxbt27RBw4cKNWuXTuf06dPZ7tppH79+ont27e/HhERYf/333+XLVOmTMqQIUOy3dBiLlOSv/LAESP6HQVKTCFMIfJiV7YULb9/Ae/YXRz5cQuhtZ8hkTtb8SvpC3T+dzyHOrxIUBAsXgxGFoLPVn5OCUk3NzSSWmNkOViI+9H58+dtX331VU83N7cmTZo0qdu9e/eaXbt2rVWrVq0GEydOrGpnZ6fHjx9/xphrffvtt2caN25888CBA6W8vb19g4KCanft2rVWzZo1fTds2OBSo0aNhDlz5pzKT7zjxo27NGnSpIi4uDjrHj16eK9cufL2zl0fH5/EL7/8MsLa2lqPHDnSo169evUfffTRmv7+/j6PP/64V2JiohozZszZLl26ZLtkOmjQoGgbGxu9YsUK12vXrtlkrO2X0fbt20v37NnTq3z58k1at27t/eijj9YMCgqq7e7u3mjTpk3OlStXThw7dmy2CWY6W1tbli1bdvLxxx+/cvToUcd27dr5nDx5Mtu6LellX8Bw/6OLi4vldhViWvJ3FahtRL9aaX2FuL8ohfdzrWl1bB43D0Sy5aEPOW99pybsd7zIunWG+wU9POC99+DMCfOXhPNzSgjImcFC3K+aNm0a/+67755p06ZNzJUrV2zWrVtXdtOmTS5KKfr27Xtpx44d+3v06GHUMq2Li0vq1q1bD7/77rtnPDw8bm3dutV5w4YNZcuWLZs8bNiw82FhYQcznsZhrhEjRlz+9ttvTyYkJKgnnniizuLFi53T24YMGRK9cePGg926dYu+fPmyzR9//FHuxIkTDoGBgdeXLVt2ZMKECTnWIa5atWpy+/btr6d/nbG2X0ZPP/301VdfffW8j49P/IkTJxxWr15dLiwsrEylSpUSx4wZc3b37t0HjNmJa21tzYIFC0717dv30qlTpxzatWvnc/jw4bt2Qnfp0iXWzs5OA7z22msWXfIFUMbuxFFKLQF6AD211stz6NMdWA4s1Vr3sliURVizZs30zp07zRobEhJCYGCgZQMqwYrj55WakET4+BXE/PgbQefmkKzvLNfakMQx6hDh3hqH4S/RbMQDWFmbnsSle3vZXuaGmr+07GhrxSc9GxHsVyIOMTBLcfw3Vpjy83kppcK01s2M6RseHn6qcePGFqlvJkRxMG3aNNeXX365ZqtWrWK3bdtmzKprtsLDwys0btzYM+vrpsz8TU57XqSUmqWUaq+UqqGUqp7255nAbxgmFCbnfBkh7h9W9rb4fdST9md+4fgpa95+GyqnnTIZzDI8iKTd2fm0eLM9xxx9WRs8lSsnrud+0Ryk1wc0Z1MISKFoIYQoCuLj49WkSZOqALz22mu5LiWby+jkT2u9GRgOKKA/sB44CZxK+/PzadcbrrXeYvFIhSjmatSAjz6CyEjDbuGnq23K1O6dtJ+g5cNwqF2VDXUGsXvWLrPKxaQngR2qm7cpRApFCyHEvff5559X7NWrl2e9evUaHD9+3KFly5axTz31lHmzAXkw6W5xrfUUoAUwF4gEkoGUtD//DLTQWk+1dJBClCS2tob7/nqe/ppTS/9jq+8QbnBnA1tp4uhw/AeaDPBnX+mWrOs7m5gLpm+U69/AgVOfPpKvM4M9x6ySc4OFEOIeWL9+vdOSJUvKX79+3frRRx+NXrp06YmCei+Ttwpqrf/TWvfXWtfUWttrre3S/vyc1vq/gghSiJLKM7gJrfd8h3XUObb1/YZjjg0ztfvGb6fT3Of5qMYMXnwRwsNNf4/8nBkMd5aDJQkUQoiC89dff53QWoddvXo1fPny5SerVKmSnPco85hXJ0IIYVGOlZwJ+Pll6tzcw+EfNvFv7WdIwLABLB4Hfkjsy/ffQ5MmEBAAP87S3Lxq2k7h/JwZDHJPoBBClBSS/AlRlCiFz4C2tDw2j4RjZ9jS4zN+qvAG1yh3u0toKHw7YCfx5d1Z2+QNDi47bNJb5OfMYLhzT6DMBAohRPGUbXVpAKXU9HxcV2uth+RjvBD3PefaFWmzbBStNdTfBNOm3SkQPYgZVNCXCQqfBI9NIsypPTFPDqL5p70oUyHP89gJ9nMn2M+dZf+dZfTiPSQkm14/NH0mcOySPfd9iRhR8LTWKGV+KSQh7je5lfLLMfkDBubnPQFJ/oSwAKWgXTvD4+JF+Hl2KoHvbIUMq77+sRth5kaiZ77KuoZ9qfzOIHDL+9qSBIriQCkVn5KSYmVjY2PRUw6EKMmSk5NtlFLZnmySW/I3qIDiEUKYyc0N3hhlhR4Zzp5Jf5MwdQZ+Z37HhhQAXImm076voPdXhDm0Yv2TL9Lssydwrpz7Ge2WTALfXLSbiU80kSRQWNKemzdvNnZxccn2B5kQ4m4xMTFltNZrsmvLMfnTWs8suJCEEPmhbKxpNLorjO5K9P7zHBz1IzXW/ED1pJO3+/jfCoWfQ2mysBX+z/gweDC0aGGYScxJxiTwraV7uZmYYnJsSalIEigsKikp6Y/o6Ojmzs7OsvQrhBGSk5OtLl26RHJy8qLs2mXDhxDFnGuDKrRZNY5q8cfYN/lvtns8TlLa73UbaUf4LR9mzYJWraBxY/h+YgzXzuQ+gRLs587+D7vwZe8m2NuY920iPQmUgtHCAn6NiYnZe/r06fLx8fF2xh5LKsT9RGtNUlKSzZUrV8oeO3bMOT4+/gdgQ3Z9c1v2zZFSqgzQDKgIRGqt/zU/XCGEJShrKxqO6AwjOnPt8AU2DpjEmjPtIOJOn7174fSor7Ea9TnrvZ7GZfgL+A1pkeOZwpZYDgZDwehfQiOZ3FtmAoXp/P39E8PCwvpfuXKl37Vr1/pqrStgOG1KCJGBUuqG1npN2ozfBn9//2x/UzIp+VNKOQGTgH6AbdrLPwH/prW/BIwFHtdabzc3eCFE/pT1qYTL+IeZ0j6QPv/CjBkwfz7Ex6UykB9wJpaOR6fD0OkcGV6fiE4v0OCTvlRtkv0uEUssB6dimAkcvmA3ZR1tef/RBpIICqP5+/tfA75Oewgh8sHo9RylVCkgBMMu4BhgDXf/5vU3UA14zELxCSHyQSnDcu/MmXDuHPz08VlS7R0z9fFOOkDnv96gop87Wys9xubRK0i4mX1h+YzLwebWCQS4Fp8kBaOFEKKQmHIzz+uAH/ArUFNr3SVrB631ceAo0NEy4QkhLMXFBfqOq07NuAMc+mETW+u9kOlMYVuSaX1xGW0/f5SrTtVZ22w0e0Ljsr2WpZJAKRgthBD3ninJ35PAeWCA1vpmLv0iAFnLEaKIUlaKugPa0vrATGwuRfHvkFnsLds2U5/KOgrPsMU0DnCkWTP49lu4evXua6Ungfk5Ng7k/GAhhLiXTEn+agPbtda38uh3GahgfkhCiHvFoUIZWn73PL5XN3Fm7SE2txnNRevKAPzI84AiLAyGDoUqVeCDoE3s+HIzqSl330OcfmxcuVK2d7UZKz0J9Bq3SpJAIYQoIKYkf0mAvRH9qgFSiFOIYqZaJx/abv6UCnGn2T1+BdceewH7DP/HJyRAx3XjaD7iAU461GVt5884/e+5TNcI9nPnv3cf5NSnj9CnVQ2zY5EyMUIIUXBMSf6OAH5KqRwTQKVUWaAxsC+/gQkhCoeVnQ1N3urGN0uqcP48fPMN+PuDF0d4gM0A1E4+c/TwsQAAGfVJREFUQtDaMVRtVZ3trl0IGfIrMVGZ7w8cH+x7OwnMT02OuaGRkgQKIYQFmZL8LQYqARNy6TMeKANkW1FaCFG8lCsHL78MO3fC8lW2bG04mFicbrdbk0qLq6sJnP4MVKlMSJ2B7Jj0DylJd+oBjg/25eSnj+SrYDTcSQLlvkAhhMgfU74TTwEOA8OVUhuVUq+mve6hlBqklPobeAnYD/xg4TiFEIWs3sM1ab33e+yio9j+ys+EuwZmancmlsDjM2n+Rnv2l27O6FGa/fvvtAf7uXN4fNd8J4Hp9wXKbKAQQpjH6O/AaTt8HwTCgAeAL9KaAoHvgCAgHHhEa51g2TCFEEWFfblStJjSl8ZXNnBx+yk2dxnPKTuvTH22JjXn84mKhg2hWTP4+mu4dNGwSSRjEpifMjFwZzawyQd/y2ygEEIYyaRfv7XWp7XWLYBHge8xFHVej+GUj95AM631aYtHKYQoktyae9D2z7fwvHWYIz9tY0ujl7iqyvET/W/3CQuD116D1ZX7s9m9N1ve+oOEm8kWqxUIUjRaCCFMYdbZvlrrlcBKC8cihCiulMK7Xyu8+7Ui6cYXvB1ix08/w/LlkJgIZbnK43ohDucSYMJCLnxSiX2NnsXtjX70eLbx7WPe3l62l7mhkWaHkV402tHWik96NpLj44QQIhs5zvwppX5TSnVVSsnh2UIIo9mWseeRboqFCyEqCr77Dl7xXoMDd+4GqaQv0Cl8Mr59m3DEwZe1nT7hxPpTt3cIy32BQghRcHL77toTw+zeaaXUx0opr1z6CiHEXcqVgyFD4KPDTxKxci+bA97kgnWVTH18EvcRtH4ctTrVZLdTW9Y8NZOWlS2zOQRkl7AQQmSV23fVacBVoCowBjiUtsu3v1Kq1D2JTghRYng80pC2Wz+nYvxp9nz+F6G1niYOx0x9mtzYwo0FK6lWDTp3hmvh7uwY1TXfJ4eAzAYKIUS6HJM/rfVQDIlfbwwbO1Ix7PKdBUQppWYopVrfkyiFECWGla01/2/v3sOkqM48jn9fBpAREBQFFRVZRYMKIhDFGARFRcQYQJCEEEPWaDaouGt0I16SYDAgKKIRL6wXokJkEYUQESMqBG8RBZSFwXhDDIqAIgJyn3f/qGpoOzPdPdP37t/neeo5UnWq+u2TSvNyqs457a/tSZf3p2Br1/L6lY/yZote7CIY9DGZH1FZCXPnwk9/Ci1awD9//jwTdu7mnZt6pjxpNGiUsIiUtrjPU9x9h7tPc/deQCvgBoKVPhoBlwALzKzCzK41s4MzH66IFJPygxpx8l2D6bRmNpvf+YSXB93Npq69iX7TeNe2nQxc+Eu+M/I8vmp8KGfeMJHpRzTjjgEnpvxIODJKWGsJi0gpqck8f5+4+yh3bwucBjwIbAKOBUYDq8xsppl938xSm7dBREpO02Oac9rky3n2b+V89BGMGQMnnghn8xwHsR6Ag3wdZy6/m46Xn0rHU7ox4YVXuKZVM/ZNcaqYyFrCQ+ZsUW+giBS9Wv2z2d1fdfdLgUOAnwDzgDLgfOBJQL+cIlJrhx8O114LS5bAnbOPYcEpv2RN2aHfqNNq9wec9epIrhh6KjNv/SUPvPIm57doXM0Vk6feQBEpdik9M3H3re7+qLv3AM4F1gMGHJSO4ERE2vQ6mq6v3UaLbav4v7te4OVvXcJGa/LNOjuWc9aC3/DDq4ez7597c3q9DtSrk9oj4UhvoAaIiEixSenX0cwamdklZrYAmMPepE+rfIhIWlndMk648gxOq3iAfTeu4Y3h03mt1UVsYe/kA1MZSEUFPDqyJe+N6kX5kg6UV6Y+VammixGRYlKr5M/MzjCzR4A1wESCdwB3ANOAXkDrtEUoIhKjXuMGdP59P7qsnErZ+rW8fvXjvN6yL0836P+NeiuePZSnxl7NylvPp9dbC8E9pc/VdDEiUgySXt7NzFoTvN/3E+AI2DPbwhLgYeAxd9+Q9ghFROJo0KwhJ98+EG4fyNtb4OmnYepUmD0bjt32Ft/iHQDunTMC5sBTbbsz/Nwr2VavPqSwgNFjr63asxTd4C5HMLJPu7R8HxGRTIub/IWTOQ8AhhDM8WfhtgGYAjzo7ksyHKOISFIaNoSLLgq2TZvgpfG7WPSHnrRfN5e67Aagb8U8+lbMA+CGs/6DyR17ByenIRHUmsIiUgjire37IMFj3YeAbuHuucAPgUPc/UolfiKSrxo3hl43dabj2jlseW8NLw+ZyOJmPdgd9bN3y9z7WDnme4yfdRvl27dBak+F9VhYRApCvJ6/n4blh8Ak4GF3/2fGIxIRSbMmRx3IaQ9fClzKFxWfsXzkdBrPnkq7LxdQB6dPxXyWVZzGjdzC/mctpXHHVcF7LSmMFYn0BhrwIz0WFpE8Ei/5mww85O4vZisYEZFMO6BtC747eSgwlM+Xf8aK0TMonz2dmRv7wy7YMLcdG+a2Y9+2qzmk10J21a2b0iNhZ28i2LS8Hr+94Hg9FhaRnIq3tu+PlfiJSDFrdlwLTnvk53Rc/1f+uu4krr++gj59oEEDaFxRxj/G9WP8rNvY/+uNwUjhFEcLRyaQ1mNhEcmlpEf7iogUs6ZN4eyzP+OWW9qyeTP87cFNLBzbl3MqnqFPxXwAZrTtxg09L2dL/fLgpDQMEtFjYRHJNiV/IiIxGjWC865qA1dNZ+v6Lbw+9hl2/+90elT8hT4VFwFBIvirXsPYXje1KWP0WFhEsi219Y9ERIpc+YENOfnW/pz64Z/YZ+M63vz1TF45+mK6r3iLd8ZdGIwU/np7Op4Kf+OxcIcRf9VqIiKSEUr+RESSVH+/BnQacQHfefePNNqylrfHzaXZPu3p/EQTVo3pzaZFR+CVYRKo9wNFJE8p+RMRqYW65fVo/1896Lr4Lh5Z3YOlS+Hqbu1o/mJvPhvTg4GLnsX2ZIKpiawt3FqJoIikgZI/EZEUmcEJJ8CNN8LChfD+O7s5v2kbnvqfOxg3a1zaRgtH3g/UY2ERSYWSPxGRNGt5TEO6ThtGhy9eoMfTT3DP7kqmTpnEwEXPUqeyMu3TxigRFJGaUPInIpJBTVofwKkTBnPKx09w8xNjmda4jEf/Mot+i15M22Ph6ESw7U3PKBEUkbiU/ImIZMk+TRrQ6abz6Lrsfm6fM5Y5Jx5K36/3g11laRktDFpfWEQS0zx/IiI5YGV1OHZQJ+4YBHcAa9fC6MdXM/PjZeyusyOslMLiwmgiaRGpmpI/EZE80Lw5jBvWknG0ZMem7Qwb+zzPbq9kT2dgmiaSViIoIkr+RETyTP3G+3DfzecFf3DnqrsWMPOTr4I/p9gbqERQRJT8iYjkMzPuvOp07gRmLF7N9dPf5uudu/ccS0V0IggwWImgSEko+eTPzAYBvwDaA2XACuBh4F53r8xlbCIi0fqc1HLPur8zFq9mxKxlbPh6Z9quH0kEy+vVYVS/9lpjWKRIlXTyZ2YTgKHANuB5YCfQA7gb6GFmA9x9dw5DFBGpUiYTwciI4f+cukSJoEgRKtnkz8wuJEj81gCnu/u74f4WwItAX+AK4M6cBSkikoToRPDGGUuZ8vdVVKZh2hhQIihSjEp5nr/hYfmrSOIH4O6fETwGBrjOzEq5jUSkwIzs044PRvVm5ejeDO5yBHVSey3wG6LnEGx70zO88kn6HjmLSPaUZGJjZocBnYAdwLTY4+4+H1gNHAx0yW50IiLpkelEcOLbO2hz/dNaUUSkwJRk8gecFJbL3H1rNXUWxtQVESlY0Yng+IEdaFi/LC3X3VkJV09dogRQpICUavLXOiw/ilNnVUxdEZGi0Oekliy7+dw9ieA+dVP7q6ASGPvsO+kJTkQyrlQHfDQKyy1x6mwOy8axB8zsMuAygBYtWjBv3rxaBbF58+Zan1uK1F41o/aquVJss6bA/WeV88onO5lSsYPNtXyNb/WXW0uu7UQKVakmf5E3X2o1Hs7dJwITATp37uzdu3evVRDz5s2jtueWIrVXzai9aq6U26w7cH3437WZOqZl0/KSbTuRQlOqyd+msGwUp07k2KY4dUREik701DEQTB8z+bVV1f5ruQ5wbc9jsxKbiKSuVN/5WxmWreLUOTymrohISRrZpx0fjq561HB5vTqMG9hBc/+JFJBS7flbHJbHm1l5NSN+vx1TV0Sk5I3s027P+r+l/JhcpJCVZM+fu38MLALqAwNij5tZN+AwgtU/Xs1udCIiIiKZU5LJX2hUWN5qZkdHdppZc+Ce8I+j3b0y65GJiIiIZEipPvbF3Z8ws3sJlnJbamZzgZ1AD2A/YAZwdw5DFBEREUm7kk3+ANx9qJm9BFwOdAPKgBXAQ8C96vUTERGRYlPSyR+Au08BpuQ6DhEREZFsMPdazXMsITNbR/xl4uI5EFifxnCKndqrZtReNac2q5lU2quVux+UzmBEJDlK/nLIzN5w9865jqNQqL1qRu1Vc2qzmlF7iRSmUh7tKyIiIlJylPyJiIiIlBAlf7k1MdcBFBi1V82ovWpObVYzai+RAqR3/kRERERKiHr+REREREqIkj8RERGREqLkrwbMbJCZLTCzjWa22czeMLPLzaxW7VjT65nZJDPzONuK1L5heqWrvczsWDO7ysweM7MVZlYZft/+2Ywj03LdXoV2f0F62szM6plZDzO73cxeM7NPzWyHma02syfMrHs24siGXLdXId5jIsWo5Ff4SJaZTQCGAtuA59m7DvDdQA8zG+Duu7N0vZeB96rY/2myn59paW6vXwBX5UEcGZMv7RXK+/sL0tpm3YDnwv9eA7wJbAGOAy4ELjSz37n7rzMcR0blS3uFCuIeEyla7q4twUbwg+YEP0xtova3AJaHx67K9PWASeGxIblukyy318+AMcBFwFHAvPAa/bMZRwm0V0HcX+luM+BM4AmgaxXHBgK7wuudoXss5fYqmHtMm7Zi3nIeQCFswBvhD9bFVRzrFvWjWieT1yuUH850t1cV10g2mcloHEXYXgVxf2X7f1vggfB6D+oeS7m9CuYe06atmLe8eh8lH5nZYUAnYAcwLfa4u88HVgMHA12yfb18ky/fL1/iSKRQ4swnOWizxWF5WI7jqJV8aS8RyR9K/hI7KSyXufvWauosjKmb6eudYWbjzGyimf3OzHrm0Yvl6W6vQo8jkXyMM5/vL8h+m7UJy9j30fLxf7uq5Et7Rcv3e0ykqGnAR2Ktw/KjOHVWxdTN9PUurmLfcjP7gbsvTSKGTEp3exV6HInkY5z5fH9BFtvMzA4GhoR/nJ6rOFKUL+0VLd/vMZGipn9pJdYoLLfEqbM5LBtn+HpLgGHA8eF1DgXOB94iGGk318xaJhFDJqW7vQo9jkTyKc5CuL8gS21mZnWBx4AmwPPuPisXcaRBvrQXFM49JlLU1POXmIVlutbBq/X13H18zK4twNNm9hwwn+B9neHAFSlFmJp0t1dt5UscieRNnAVyf0H22uw+gqlQPgYG5zCOVOVLexXSPSZS1NTzl9imsGwUp07k2KY4dTJ1Pdx9BzAq/ON5yZyTQWn/fgUeRyJ5H2ee3V+QhTYzszuBSwjmsevh7mtyEUea5Et7VSsP7zGRoqbkL7GVYdkqTp3DY+pm83oRkZnxc/3IZGVYpvv7FWociUQ+O9/jzJf7CzLcZmZ2O8GjyXUEicy7uYgjjSKfnev2SiSf7jGRoqbkL7HItAXHm1l5NXW+HVM3m9eLaBaWm+PWyrxMfb9CjSORQokzX+4vyGCbmdkY4Grgc+Bsd1+eizjSLF/aK5F8usdEipqSvwTc/WNgEVAfGBB73My6EcxntQZ4NdvXi3JRWC6MWyvDMvj9CjKORAolTvLk/oLMtZmZjQauBTYQJDJv5SKOdMuX9kpC3txjIkUv17NMF8IG9GfvDPhHR+1vDiyjiqWRCN5fWQGMStP1OhCMiiuL2V+X4F/eu8PzehZbe1Vx/Xkkt2JFjeMo1fYqpPsrE20G/C48ZwPQKZNxlGp7Fdo9pk1bMW85D6BQNuCe8IdpKzALeBLYGO57qooftEnhsUlpul6f8NjnBP86nwbMIZiZ38Mfzv/OdTtlor2AjsBrUdtXYd1/RO9PRxyl2l6Fdn+ls82AC8L9TtDrNKma7TrdY7Vvr0K8x7RpK9Yt5wEU0gYMAl4O/zLdArwJXE4V62HG+8u5ltdrDYwHXgl/LLeFP+LvAg9Rg96KQmsvoHvUXzbVbumIo1TbqxDvr3S1GcGkxAnbC5ine6z27VWo95g2bcW4mbsjIiIiIqVBAz5ERERESoiSPxEREZESouRPREREpIQo+RMREREpIUr+REREREqIkj8RERGREqLkT0RERKSEKPkTiWJmK83Mw21UgrqTo+rOy1KIGWVmZWa21Mw+MrN9ovYfGX7PlVmIoV/4WVdk+rNEREqRkj+R6l1sZmVVHTCz/YC+WY4nG34BnAD81t235yIAd38SeAMYYWYH5CIGEZFipuRPpGpvAIcCZ1dz/AdAOcHapkXBzBoBI4APgUdyHM7NwAHA9TmOQ0Sk6Cj5E6napLAcUs3xIQQL0T+ahViy5ScECdckd9+d41hmA2uAn5lZwxzHIiJSVJT8iVTt78By4Ptm1jT6gJkdC5wKPAt8Wt0FzOwsM5tgZm+Z2edmtj18l+6PZta2mnMamNl1ZrbIzDaH53xqZq+a2UgzaxBT/2Qzm2Zmq81sp5ltNLP3zGyKmZ1Zw+88NCxr1OtnZk3M7IXwPb0ZZlYedczM7DIzW2xmW81snZk9aWbtzGxIeM6k2GuGyedkoAkwqIbfQ0RE4lDyJ1K9SUAD4Icx+4eE5cMJzr8PuATYBSwg6M3aAVwMvGFm342ubGZ1gKeBUcC/AfOB6QRJ6OHADUDTqPpnAy8B/YG1wFPAC8CGcN9FSX5PzKwNcBzwnruvrMF5h4cxnAHcA/Rz961RVe4PtxOAl4G5QDuC5LpzgsvPDcvvJxuPiIgkVjfXAYjksUcJErEhwL0QjIYlSN6+AP4MXBDn/GuAee7+ZWSHmRlwGUFiONHMjnd3Dw9/FzgTWASc7u5bYs77DvBV1PWHA/WAQe7+p+gPNrNmwJE1+K7dw/LVZE8wsxMJEtpDgF+5+5iY432AS4EvgR7uvijcXwe4laB94nkNcOB0M6vr7ruSjU1ERKqnnj+Rarj7GmAOcHLUY9pzCAaCTHH3HQnOnxGd+IX73N3vB14B2hL0tkW0CMsF0Ylf1Hkvu/vXVdR/porP/tzd34z/Db+hQ1hWJFPZzM4h6M08kCD5HFNFtWFheXsk8QtjqyQYyPFxvM8I2+5ToDFwVDJxiYhIYkr+ROKbFJZDYspJJMHMDjOzn5vZHWb2oJlNCt9xOzisckxU9UUEg0guMbOhZtYi9noxXg/LKWZ2WnXT0iSpeVh+nqiimQ0heDy9GzjH3R+vok5dgp5KgCmxx919J8Ej7US+CMtEbSEiIknSY1+R+P5MkBD92MzGErx/tjSZXjUzG0HQwxXv/2f7Rf7D3d83s/8CbgMmABPM7AOCXsKZwFMxo3CHE/TY9Qq3LWb2JsF7f4+6+wfJf02ahOVXcWvBYQTvOjpwrrv/vZp6BwL7AJVU38P3URJxReJpGreWiIgkTT1/InGEj3anELzX9jBBQpNooAdmdiHwa2ArwXtvRwH7uru5uwGRd/Qs5vP+ALQimGx5MlAGDAamEQwSiU4W1wCdgB7AaIKew1OA3wLvmNm/1+CrRh5P7xe3VjCwZE4Y9/jYkdDV8Gr2VyZxbiSeDUnUFRGRJCj5E0lsUlieTzByd3IS5wwIy+vd/QF3/yBmFOzR1Z3o7mvc/T53H+zuRxL07i0Ny+ti6la6+wvuPtzdTweahXXqEvQcJkrmItaGZbME9XYQ9H7OALoAL5jZgVXU+xzYTvAbc3g11zoyibgi8ayNW0tERJKm5E8kgXCwwksECc00d08mEYksS/YvjzzDwSMn1eDz3wLuDP94YoK6W9z9VuCfBNPUHJvkx0QGZBwXtxZ7ekMHEPRengTMM7ODY+rsJBitC/86VQ5mVg+4MN7nmNn+BO9GfgW8lyguERFJjpI/kSS4e1d3P9Ddk51weEVYXmpm9SM7zaw58EeqeA/QzM40s/PCwRLR+8uA88I/fhS1/5pwnr3Y63QmeExdSZAEJuPFsDw1mcrhtCuDgYeA44H5ZnZYTLU/hOU1ZhYZTRyZ6mUkcESCj+lC8Hh5QR6sOCIiUjSU/IlkxnhgI9AbeC9cheMvwPtAI4LHprHaE4yiXR+umDHZzJ4i6D3sR7Dc2a1R9W8EVpnZcjObHq7qsYBgAuUyYIy7V7sCSTR3/xB4GzjKzFoneU4l8DPgboJRy3+LPtfdpxMkh/sDC83sOTObQpAYDyOcO5HgUXJVzgrLmcnEIyIiyVHyJ5IB4UjbjsDjBL1X3yOY128iQe/axipOmwWMIHgEezTBY9GuBEnfb4D27h49QvZygl7ESoIVNvoCLcPr9HT34TUM+56wvDjZE8L5B68ExgCtCRLA6OlrLiUYvLIs/C49CeYS7AJ8EtZZH3vdsLdzEEE7/ctUMSIiUnu2d3EBESllZtaQ4LHyV0CbTD9qNbO5BCOV+4e9hNHHLiDo8bvd3ROtBCIiIjWgnj8RAYLBIgQ9jK2pQe9fPGZ2vJntG7OvnpndSJD4rSNYIi7WTQQTPP8+HXGIiMhe6vkTkT3Cx62LCSZ9Psbdt6d4vccIHkcvAlYTTNbcjmCJvO1AP3efHXNOX+BJYFg476GIiKSRkj8RyRgz603w3l9Hgjn76hKs1zsfuM3dl+YwPBGRkqTkT0RERKSE6J0/ERERkRKi5E9ERESkhCj5ExERESkhSv5ERERESoiSPxEREZES8v9pxhQYVj/21gAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"#Plot all 3 solutions\n",
"plt.plot(m,num_sol_heun[:,1], 'b', label = 'Heun');\n",
"plt.plot(m, num_sol_euler[:,1], 'r--', label = 'Euler');\n",
"plt.plot(m,-np.log(m/0.25)*250, 'o', label = 'Tsiolkovsky')\n",
"plt.xlabel('Mass (kg)')\n",
"plt.ylabel('Velocity (m/s)')\n",
"plt.grid(True);\n",
"plt.legend(loc='center left',bbox_to_anchor=(1,0.5));\n",
"\n",
"print('According to the graph, the solutions converge')"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0.00000000e+00 0.00000000e+00 2.57216000e-03 7.72160343e-03\n",
" 1.54534606e-02 2.57728686e-02 3.86849713e-02 5.41949195e-02\n",
" 7.23078706e-02 9.30289887e-02 1.16363445e-01 1.42316416e-01\n",
" 1.70893088e-01 2.02098649e-01 2.35938299e-01 2.72417240e-01\n",
" 3.11540684e-01 3.53313847e-01 3.97741953e-01 4.44830232e-01\n",
" 4.94583920e-01 5.47008260e-01 6.02108500e-01 6.59889897e-01\n",
" 7.20357712e-01 7.83517213e-01 8.49373675e-01 9.17932376e-01\n",
" 9.89198605e-01 1.06317765e+00 1.13987482e+00 1.21929541e+00\n",
" 1.30144473e+00 1.38632811e+00 1.47395086e+00 1.56431831e+00\n",
" 1.65743580e+00 1.75330866e+00 1.85194225e+00 1.95334191e+00\n",
" 2.05751301e+00 2.16446090e+00 2.27419096e+00 2.38670856e+00\n",
" 2.50201907e+00 2.62012789e+00 2.74104041e+00 2.86476202e+00\n",
" 2.99129813e+00 3.12065414e+00 3.25283546e+00 3.38784752e+00\n",
" 3.52569573e+00 3.66638552e+00 3.80992233e+00 3.95631159e+00\n",
" 4.10555876e+00 4.25766927e+00 4.41264858e+00 4.57050215e+00\n",
" 4.73123544e+00 4.89485391e+00 5.06136305e+00 5.23076833e+00\n",
" 5.40307522e+00 5.57828923e+00 5.75641582e+00 5.93746052e+00\n",
" 6.12142880e+00 6.30832618e+00 6.49815816e+00 6.69093026e+00\n",
" 6.88664799e+00 7.08531687e+00 7.28694244e+00 7.49153021e+00\n",
" 7.69908572e+00 7.90961451e+00 8.12312211e+00 8.33961408e+00\n",
" 8.55909596e+00 8.78157330e+00 9.00705165e+00 9.23553658e+00\n",
" 9.46703365e+00 9.70154842e+00 9.93908646e+00 1.01796533e+01\n",
" 1.04232546e+01 1.06698959e+01 1.09195828e+01 1.11723209e+01\n",
" 1.14281157e+01 1.16869728e+01 1.19488978e+01 1.22138964e+01\n",
" 1.24819742e+01 1.27531366e+01 1.30273894e+01 1.33047381e+01\n",
" 1.35851883e+01 1.38687457e+01 1.41554159e+01 1.44452045e+01\n",
" 1.47381170e+01 1.50341592e+01 1.53333366e+01 1.56356548e+01\n",
" 1.59411196e+01 1.62497364e+01 1.65615110e+01 1.68764489e+01\n",
" 1.71945558e+01 1.75158373e+01 1.78402991e+01 1.81679468e+01\n",
" 1.84987859e+01 1.88328223e+01 1.91700614e+01 1.95105089e+01\n",
" 1.98541705e+01 2.02010518e+01 2.05511585e+01 2.09044961e+01\n",
" 2.12610703e+01 2.16208867e+01 2.19839511e+01 2.23502690e+01\n",
" 2.27198460e+01 2.30926878e+01 2.34688001e+01 2.38481885e+01\n",
" 2.42308585e+01 2.46168160e+01 2.50060664e+01 2.53986154e+01\n",
" 2.57944687e+01 2.61936318e+01 2.65961105e+01 2.70019103e+01\n",
" 2.74110369e+01 2.78234959e+01 2.82392930e+01 2.86584336e+01\n",
" 2.90809236e+01 2.95067684e+01 2.99359737e+01 3.03685452e+01\n",
" 3.08044884e+01 3.12438090e+01 3.16865125e+01 3.21326045e+01\n",
" 3.25820907e+01 3.30349767e+01 3.34912681e+01 3.39509704e+01\n",
" 3.44140892e+01 3.48806301e+01 3.53505988e+01 3.58240008e+01\n",
" 3.63008416e+01 3.67811268e+01 3.72648621e+01 3.77520529e+01\n",
" 3.82427049e+01 3.87368235e+01 3.92344144e+01 3.97354830e+01\n",
" 4.02400349e+01 4.07480757e+01 4.12596109e+01 4.17746459e+01\n",
" 4.22931864e+01 4.28152378e+01 4.33408057e+01 4.38698954e+01\n",
" 4.44025127e+01 4.49386628e+01 4.54783513e+01 4.60215838e+01\n",
" 4.65683656e+01 4.71187022e+01 4.76725990e+01 4.82300616e+01\n",
" 4.87910953e+01 4.93557057e+01 4.99238980e+01 5.04956778e+01\n",
" 5.10710504e+01 5.16500213e+01 5.22325958e+01 5.28187793e+01\n",
" 5.34085773e+01 5.40019950e+01 5.45990378e+01 5.51997111e+01\n",
" 5.58040203e+01 5.64119705e+01 5.70235673e+01 5.76388158e+01\n",
" 5.82577214e+01 5.88802894e+01 5.95065250e+01 6.01364336e+01\n",
" 6.07700203e+01 6.14072905e+01 6.20482493e+01 6.26929020e+01\n",
" 6.33412538e+01 6.39933100e+01 6.46490756e+01 6.53085559e+01\n",
" 6.59717561e+01 6.66386812e+01 6.73093365e+01 6.79837270e+01\n",
" 6.86618580e+01 6.93437344e+01 7.00293614e+01 7.07187440e+01\n",
" 7.14118873e+01 7.21087964e+01 7.28094763e+01 7.35139320e+01\n",
" 7.42221685e+01 7.49341908e+01 7.56500038e+01 7.63696126e+01\n",
" 7.70930220e+01 7.78202370e+01 7.85512625e+01 7.92861034e+01\n",
" 8.00247646e+01 8.07672510e+01 8.15135672e+01 8.22637183e+01\n",
" 8.30177090e+01 8.37755442e+01 8.45372284e+01 8.53027666e+01\n",
" 8.60721635e+01 8.68454238e+01 8.76225521e+01 8.84035532e+01\n",
" 8.91884318e+01 8.99771924e+01 9.07698397e+01 9.15663783e+01\n",
" 9.23668127e+01 9.31711476e+01 9.39793875e+01 9.47915368e+01\n",
" 9.56076002e+01 9.64275820e+01 9.72514867e+01 9.80793188e+01\n",
" 9.89110826e+01 9.97467826e+01 1.00586423e+02 1.01430008e+02\n",
" 1.02277543e+02 1.03129031e+02 1.03984477e+02 1.04843884e+02\n",
" 1.05707258e+02 1.06574603e+02 1.07445922e+02 1.08321219e+02\n",
" 1.09200500e+02 1.10083767e+02 1.10971025e+02 1.11862279e+02\n",
" 1.12757531e+02 1.13656787e+02 1.14560050e+02 1.15467323e+02\n",
" 1.16378611e+02 1.17293919e+02 1.18213248e+02 1.19136604e+02\n",
" 1.20063991e+02 1.20995411e+02 1.21930869e+02 1.22870368e+02\n",
" 1.23813912e+02 1.24761505e+02 1.25713151e+02 1.26668852e+02\n",
" 1.27628613e+02 1.28592436e+02 1.29560326e+02 1.30532286e+02\n",
" 1.31508319e+02 1.32488429e+02 1.33472619e+02 1.34460893e+02\n",
" 1.35453252e+02 1.36449702e+02 1.37450245e+02 1.38454884e+02\n",
" 1.39463623e+02 1.40476464e+02 1.41493411e+02 1.42514466e+02\n",
" 1.43539633e+02 1.44568915e+02 1.45602314e+02 1.46639834e+02\n",
" 1.47681477e+02 1.48727246e+02 1.49777143e+02 1.50831172e+02\n",
" 1.51889336e+02 1.52951636e+02 1.54018076e+02 1.55088658e+02\n",
" 1.56163385e+02 1.57242259e+02 1.58325282e+02 1.59412457e+02\n",
" 1.60503787e+02 1.61599273e+02 1.62698918e+02 1.63802725e+02\n",
" 1.64910694e+02 1.66022829e+02 1.67139132e+02 1.68259604e+02\n",
" 1.69384248e+02 1.70513066e+02 1.71646059e+02 1.72783229e+02\n",
" 1.73924579e+02 1.75070110e+02 1.76219823e+02 1.77373721e+02\n",
" 1.78531805e+02 1.79694076e+02 1.80860536e+02 1.82031188e+02\n",
" 1.83206031e+02 1.84385067e+02 1.85568299e+02 1.86755726e+02\n",
" 1.87947351e+02 1.89143174e+02 1.90343197e+02 1.91547420e+02\n",
" 1.92755845e+02 1.93968473e+02 1.95185304e+02 1.96406339e+02\n",
" 1.97631580e+02 1.98861027e+02 2.00094680e+02 2.01332541e+02\n",
" 2.02574609e+02 2.03820886e+02 2.05071371e+02 2.06326066e+02\n",
" 2.07584970e+02 2.08848085e+02 2.10115409e+02 2.11386943e+02\n",
" 2.12662688e+02 2.13942643e+02 2.15226809e+02 2.16515184e+02\n",
" 2.17807770e+02 2.19104566e+02 2.20405572e+02 2.21710787e+02\n",
" 2.23020211e+02 2.24333843e+02 2.25651683e+02 2.26973731e+02\n",
" 2.28299986e+02 2.29630447e+02 2.30965112e+02 2.32303983e+02\n",
" 2.33647057e+02 2.34994333e+02 2.36345811e+02 2.37701489e+02\n",
" 2.39061366e+02 2.40425441e+02 2.41793713e+02 2.43166180e+02\n",
" 2.44542841e+02 2.45923694e+02 2.47308737e+02 2.48697970e+02\n",
" 2.50091389e+02 2.51488994e+02 2.52890782e+02 2.54296751e+02\n",
" 2.55706900e+02 2.57121226e+02 2.58539727e+02 2.59962401e+02\n",
" 2.61389246e+02 2.62820258e+02 2.64255435e+02 2.65694776e+02\n",
" 2.67138277e+02 2.68585935e+02 2.70037748e+02 2.71493712e+02\n",
" 2.72953826e+02 2.74418085e+02 2.75886487e+02 2.77359027e+02\n",
" 2.78835704e+02 2.80316514e+02 2.81801453e+02 2.83290517e+02\n",
" 2.84783704e+02 2.86281008e+02 2.87782427e+02 2.89287957e+02\n",
" 2.90797593e+02 2.92311331e+02 2.93829168e+02 2.95351099e+02\n",
" 2.96877120e+02 2.98407227e+02 2.99941414e+02 3.01479677e+02\n",
" 3.03022012e+02 3.04568414e+02 3.06118877e+02 3.07673398e+02\n",
" 3.09231971e+02 3.10794590e+02 3.12361251e+02 3.13931948e+02\n",
" 3.15506677e+02 3.17085430e+02 3.18668203e+02 3.20254991e+02\n",
" 3.21845786e+02 3.23440584e+02 3.25039378e+02 3.26642163e+02\n",
" 3.28248931e+02 3.29859678e+02 3.31474396e+02 3.33093079e+02\n",
" 3.34715720e+02 3.36342313e+02 3.37972851e+02 3.39607327e+02\n",
" 3.41245734e+02 3.42888066e+02 3.44534314e+02 3.46184472e+02\n",
" 3.47838532e+02 3.49496487e+02 3.51158329e+02 3.52824050e+02\n",
" 3.54493643e+02 3.56167100e+02 3.57844412e+02 3.59525573e+02\n",
" 3.61210572e+02 3.62899403e+02 3.64592056e+02 3.66288524e+02\n",
" 3.67988797e+02 3.69692867e+02 3.71400724e+02 3.73112361e+02\n",
" 3.74827767e+02 3.76546934e+02 3.78269853e+02 3.79996513e+02\n",
" 3.81726905e+02 3.83461020e+02 3.85198849e+02 3.86940380e+02\n",
" 3.88685604e+02 3.90434512e+02 3.92187092e+02 3.93943335e+02\n",
" 3.95703230e+02 3.97466766e+02 3.99233933e+02 4.01004720e+02\n",
" 4.02779117e+02 4.04557111e+02 4.06338693e+02 4.08123850e+02\n",
" 4.09912571e+02 4.11704846e+02 4.13500661e+02 4.15300006e+02\n",
" 4.17102869e+02 4.18909238e+02 4.20719100e+02 4.22532443e+02]\n"
]
}
],
"source": [
"#Calculate heights\n",
"height = num_sol_euler[:,0]\n",
"print(height)"
]
},
{
"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": 42,
"metadata": {},
"outputs": [],
"source": [
"def f_m(dmdt,m0=0.25, c=0.18e-3, u=250):\n",
" ''' define a function f_m(dmdt) that returns \n",
" height_desired-height_predicted[-1]\n",
" here, the time span is based upon the value of dmdt\n",
" \n",
" arguments:\n",
" ---------\n",
" dmdt: the unknown mass change rate\n",
" m0: the known initial mass\n",
" c: the known drag in kg/m\n",
" u: the known speed of the propellent\n",
" \n",
" returns:\n",
" --------\n",
" error: the difference between height_desired and height_predicted[-1]\n",
" when f_m(dmdt)= 0, the correct mass change rate was chosen\n",
" '''\n",
" \n",
" mf = 0.05\n",
" dm = m0-mf\n",
" T = dm/dmdt\n",
" g = 9.81\n",
" dt = T/100\n",
" \n",
" N = 100\n",
" t = np.linspace(0, T, N)\n",
"\n",
" #Initial conditions of position and velocity\n",
" y0 = 0\n",
" v0 = 0\n",
"\n",
" #Initialize array\n",
" num_sol = np.zeros([N,3])\n",
"\n",
" #IC's\n",
" num_sol[0,0] = y0\n",
" num_sol[0,1] = v0\n",
" num_sol[0,2] = m0\n",
"\n",
"\n",
" for i in range(N-1):\n",
" num_sol[i+1] = heun_step(num_sol[i], rocket, dt)\n",
" \n",
" error = num_sol[-1,0]-300\n",
" \n",
" return error\n"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"number of brackets: 1\n",
"\n"
]
},
{
"data": {
"text/plain": [
"array([[0.10833333],\n",
" [0.05 ]])"
]
},
"execution_count": 43,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def incsearch(func,xmin,xmax,ns=50):\n",
" '''incsearch: incremental search root locator\n",
" xb = incsearch(func,xmin,xmax,ns):\n",
" finds brackets of x that contain sign changes\n",
" of a function on an interval\n",
" arguments:\n",
" ---------\n",
" func = name of function\n",
" xmin, xmax = endpoints of interval\n",
" ns = number of subintervals (default = 50)\n",
" returns:\n",
" ---------\n",
" xb(k,1) is the lower bound of the kth sign change\n",
" xb(k,2) is the upper bound of the kth sign change\n",
" If no brackets found, xb = [].'''\n",
" x = np.linspace(xmin,xmax,ns)\n",
" f = np.zeros(ns)\n",
" for i in range(ns):\n",
" f[i] = func(x[i])\n",
" sign_f = np.sign(f)\n",
" delta_sign_f = sign_f[1:]-sign_f[0:-1]\n",
" i_zeros = np.nonzero(delta_sign_f!=0)\n",
" nb = len(i_zeros[0])\n",
" xb = np.block([[ x[i_zeros[0]+1]],[x[i_zeros[0]] ]] )\n",
"\n",
" \n",
" if nb==0:\n",
" print('no brackets found\\n')\n",
" print('check interval or increase ns\\n')\n",
" else:\n",
" print('number of brackets: {}\\n'.format(nb))\n",
" return xb\n",
"\n",
"incsearch(lambda x: f_m(x), 0.05, 0.4, ns=7)"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"number of brackets: 1\n",
"\n",
"[[0.06060606]\n",
" [0.05707071]]\n"
]
}
],
"source": [
"mn = 0.05 \n",
"mx = 0.4\n",
"\n",
"mass_change_rates = incsearch(lambda dmdt: f_m(dmdt,m0=0.25, c=0.18e-3, u=250),mn,mx,ns=100)\n",
"print(mass_change_rates)"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {},
"outputs": [],
"source": [
"def mod_secant(func,dx,x0,es=0.0001,maxit=50):\n",
" '''mod_secant: Modified secant root location zeroes\n",
" root,[fx,ea,iter]=mod_secant(func,dfunc,xr,es,maxit,p1,p2,...):\n",
" uses modified secant method to find the root of func\n",
" arguments:\n",
" ----------\n",
" func = name of function\n",
" dx = perturbation fraction\n",
" xr = initial guess\n",
" es = desired relative error (default = 0.0001 )\n",
" maxit = maximum allowable iterations (default = 50)\n",
" p1,p2,... = additional parameters used by function\n",
" returns:\n",
" --------\n",
" root = real root\n",
" fx = func evaluated at root\n",
" ea = approximate relative error ( )\n",
" iter = number of iterations'''\n",
"\n",
" iter = 0;\n",
" xr=x0\n",
" for iter in range(0,maxit):\n",
" xrold = xr;\n",
" dfunc=(func(xr+dx)-func(xr))/dx;\n",
" xr = xr - func(xr)/dfunc;\n",
" if xr != 0:\n",
" ea = abs((xr - xrold)/xr) * 100;\n",
" else:\n",
" ea = abs((xr - xrold)/1) * 100;\n",
" if ea <= es:\n",
" break\n",
" return xr,[func(xr),ea,iter]\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/opt/conda/lib/python3.7/site-packages/ipykernel_launcher.py:25: RuntimeWarning: divide by zero encountered in double_scalars\n",
"/opt/conda/lib/python3.7/site-packages/ipykernel_launcher.py:27: RuntimeWarning: invalid value encountered in double_scalars\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"The root function is inf [-300.0, nan, 28]\n"
]
}
],
"source": [
"#Find the root\n",
"n=np.arange(3,30)\n",
"\n",
"err_modsec=np.zeros(len(n))\n",
"\n",
"for i in range(0,len(n)):\n",
" \n",
" root,out = mod_secant(f_m,0.001,1,es=0,maxit=n[i])\n",
" \n",
"\n",
"print('The root function is ', root,out)"
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA5wAAAJ7CAYAAACPhJ5fAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd5gkVfm38fvZSAbJURZBQFBEWQVEAck/EBVY0pLTAisoQaKAZFBBFgmLILDgIjmDZAExu0RfBBVwCRKEJe4CG8/7x6lxmma6J/ZWz8z9ua66qrrrnOpnpnvmmu9U1TmRUkKSJEmSpJ42oOwCJEmSJEl9k4FTkiRJktQQBk5JkiRJUkMYOCVJkiRJDWHglCRJkiQ1hIFTkiRJktQQBk5J6qCIGB8RqVhu6GCf44v2E7v52g8UxxnXnePMLr2t3kaq+F50dTm+4lgtz+1e3lfUO0TEYhFxTkT8KyI+rPjefbuD/YfVeU8mF8cdHxEbNvpr6aqIWKGi5q+WXY+k/snAKUkdEBHzAJV/qG4REQt285iVf9Cu343j7N5ynO7UI/UVxc/rH4ADgBWAoT38EnMXx90JuDciLo6I6OHX6PUi4vTid9PTZdciqTwGTknqmBHkPzJbDAF2KKkW9S7/B8xbY3mhaPO7Om1Onc319gUjgU8BCdgbWJLW7+etXTjeaRX95wNWBb4DvFns3xM4pHslS1LfNKjsAiSpl9i1WD8HTAU+Uzx3/ux48ZTS+rPjddTzUkof1NpXcVZ6ZkppcgeO5Vm0jvl8sX4ipXRxDxxvWtX783fg7xHxZ+Av5H/gHwqc2QOvJUl9imc4JakdEbEMsF7x8JfFArBmRHy6nKok1TFXsX67kS+SUnoYuLd4uERELN/I15Ok3sjAKUnt24XW35fji6XlzNSubfZoRzGI0L8rnrq/jYFJhlW0/9ggPC33gAKXVjxXfYyJFfvWb+vYbdTW7kBHEbFQRJwZEc8UA7K8GhE3dWZgkogYFBF7RsSdRf9pEfF6RNwdETt15Z64iBhd1D4zIpZqp+26Fd+Pjav2LRsRZ0XEExHxXlHbyxHxWERcGBFbdba2nlBv0KDqz0hEbBwRt0fEaxExJSIej4j9I2JARZ8lIuKMiPhnRHxQvA+XFf9kaa+WJYt79B6LiLeLz8Fzxf2Mq/TA1zokIg6IiAcj4o2ImFq8BzdExBY1+kwsfiZ2L55ar+rnYVx362pD5f2Ji9ZrGBHrRcSVEfFi8fW8FRF/iojDI2Luen2L/nMW35N7ivdqakS8EhF/joiTI2LVzhYfEZtGHgQpRcRvImLeNtosFBE/jIi/RMSbxeu+EBFXRMSX2mi/WfE+HFE8tVIbv5vu7GytknonL6mVpPbtUqz/mFJ6BvIf98DXgZ0j4riUUr8ZsCciVgLuB5aoeHox4FvAlhGxfweO8UngFlovfWyxMLBxsewUEdumlKZ0oryrgTHAYPJ9fD+p03anYv0KcF9FbesBt/PRe3Yhf71LFDWPBObpRF2zVUT8ADi56unVyJeAfxHYJyK+APwaWLyizRzkf6JsFBFrppReqnH8HYBLgDmrdi1XLLtFxP4ppYu6WP8ywB3keyUrLQFsBWwVEVcCu6eUpnXlNXpQ5T/v32qrQRHyxwAHVu0aAqxZLKMjYrOUUpsD7BTv143AslW7Fi+WLwObAcM7WnhE7AhcRv55uQEYmVKaWtVmU+AqYIGq7suQfw52jIhjUkreayypTZ7hlKQ6IuLLwMrFw19W7GrZHgas24VDr8JH/5jenI8PFvN8O8d4vmi3X8Vz1cfo9pmmShExJ3Ab+Q//acBxwKeBRYBNgEeBc8gjeNY6xnzAb8jBbRL53rfPAJ8AVgSOAj4gD7bz887Ul1KaRA4qADvXqWEIsG3x8MqU0qzi+QHA5eSw+Sw5fK0ALASsBGwE/Ij235syrQecRA4JXybXvho54APsXZyhvRF4l/x9WAJYCvgeMJ08yM6P2zp4RHwD+BU5bP6eHACXLl7nq8DNwEDg5xGxSWeLj4ih5MC/KjCT/E+DVcn/jPhKxdexI3BWVfdVyJ/7K4rH1YMx7dvZejqg5ffD+8DEGm2OoTVs/hbYgPwz8+li34fkIHlXRFQHOyJfqnt/0eZ94BTgC+Tv+VLApsDPqBF42xIRB5K/T4OBi4Dt2giba5EHWVoAeIz8PV8WWBD4EvlnJYBTIqLy5+0e8vf7p8Xjf/Lx300dmp5GUh+QUnJxcXFxqbEA55Ivn50KLFTx/LzAlGLfxXX6H1+0mdjGvmHFvgSs304dDxTtxrWxb/eW47RzjPUrXm9YF2v+fsUxRraxf27gbxVt2qr3Z8W+d4GVatSwUcUxhnfyPRtR0fdzNdp8u6LN6hXPf67i+c/Phs/XxOK1Huhg+5badq/zGUnA2Db2DyEPepXIofL5ys90RbuTKz7z81XtmwN4rdh/BzCgRp2XF23+1oXvyUEVX8eoNvYHOUy3tFmtjTbjOvN9rVFH5c/n8TXarAHMKtqcU6PNkuR/ziTymfTBbbT5ZsVr/biN/fcW+94Hvlyn5kFVj1eoOO5X23iPE3BqjWMNIF8unMgDIw2t0e4nRZv/AEOq9p1e7Hu6kT9HLi4uzb14hlOSaoiIwcD2xcNfp3z2DICU0nvkMzkAI4ozf/3B7sX6zymlX1XvTPny1yNrdS7uU9ureHhiSukfbbVLKd1LPgsKrZe+dtStwDvt9G05G/NkSumxiucrbzX5Tydft1m8TxvvQcqXnt5QPBxE/v5Pqm5HDnOQA+rqVft2JN+nOIscemfVqOEHxfqzEVF92XR79i7Wj6SULqzemVJKtJ6JrWzfSEMiYp5imTciPlNcOn43OQD/jtqf+13JZxEBDkwpTa9ukFK6hXx5M8CeEa33Lxf3ZW5YPDwtpfSXWkWmlGbU+yIiYkBE/Jz8/iTg4JTS0TWab0o+qw+wZ6o6+1nhBPI/J5Yk32YgSR9h4JSk2rYgX8YHH72ctsXlxXo+8v2LfVpEfILWS3RvrNP0TnLoactXaB1B9MGKP+I/tgBPFO06fE8aQPGH8bXFw5GVf7wXX8f85PcW8gBQlf5BvpwX4LKIWLEzr90k/pRSeqfGvmcrtu/qQJvFq/ZtVKwfB6bUee/eAl4v2nbmnsLKz9h1tdqllF4jX5oK8LWOHr8bjgLeK5Z3ydOinE++tPRo4Oup9r3GLQNp/T2l9Pc6r9HymV2Ij14Kv2HF9rjOlf0RQ4vXGAXMAHZNKY2p077lvX4ZmFjnvQb4V7Hu1M+qpP7BwClJtbWMQPsW+Z6yavcAr1a17cuWJZ/NgY+OzPkRKaWZtP4BWm2liu2/0PpHfFvLQUW7RbpQa0uQrJzSpsUI8qWhidZ7/Vpqf58cLiDfV/uPiHg6Ii6KiF0iYgma38t19lXOCfpKWw3SR+cNrT5z3/L+fYH67917tL5vnXn/PknrZ6xeOAN4slhXD6Izux0FrFVnf0t9Hf16KvsAtEy18mZK6cVO1lbpEmBr8j+DvpVSqv5nS7WW93pJ2n+vP1u07crPqqQ+zsApSW2IiAVpPQv2e+AzEbF65UK+36/lLMsmEbFYGbXORpWjsk5up22t/fN34XXn6EKf39I6sE/14EEtl9k+2NYf8Cmls8kD4fyJHEpXIl+2eTnwUkTc2uRnPmd2pFHxj4H2VE9N0+j3r3JKjvY+Y++10adRTkgpRUopyF/PyuRBlWYVr399RNSaEqWlvo5+PZV9IF9BUb2/Kz5RrGfSscGFZtfPqqQ+zsApSW3bgXwPG8A3yKOvtrVsV7QZSOfvNZzdOjp1S60psyr/YG5vSpBa+yuPMW/LH/HtLMM6WPf/FPf5tdxjOqIY+ZTIc3O2nPGseYYnpXRTSmlt8nQv3wbOIJ+hGkD+PPwp6sxl2oe1vH/XdfC9i5TS8Z04fmWo6uhnrLtBrFNSSlNTSv9IKR0BHFw8vSh5gJy2tNTXmZ+Z99rY7m6w3pN8yfi8wJ0RsXY77Vve69914r3er+4RJfVLBk5JaltXLpHdpf0mpfqwYrveIEdL1nj+eVpD68o12hARA8nTPbTluYrt6gFpelrLfbfzk0Mi5HkDB5C/FzXvEWyRUno9pXRzSumwlNKq5EFzZpHPFh1Uv3ef1PL+faFBx3+B1s9Ye1P6tEwrNLFBtXTEOcAfiu1dizlqq00s1h39eir7ADxTrBcs5iftqlfJI1X/g3zW9M5i2pNaWt7rz0WE87ZL6jIDpyRViYhPkydihzwqZN3/6gOHFW1Xj4jPdeKlKkerHNiNkv93nCLs1VJ5z15bfxi39N+orX0ppbdovQ9tqzqvsxmtAwNVe5A8oiXAHnWO0W0ppaeAR4qHO1etb6szsE69Y14F/L/i4We6V2GvdHexXj4ienywnuIz1nIv4za12hWXr7bMf/u7nq6jo4oz6S33/A6kdXTeSi31rRIR9T4zI4r1JOCpiufvrdjerSt1tkgpVYfOu+qEzpb3en7q/7zX0/K7qTu/3yT1cgZOSfq4yj/qruxA+6vIZ72gc2dG36L1bE6ts4odUTm1Rc3jpJReoDV01vrD9UjyQDu1jCvWa0bEyOqdxbQntS4tJKX0LnBx8XD3iKgZKorjzdfNgXpaznJuXgSk1aqer369pSpG3mxr/5y0fo/bmlKkrxsP/LfY/kV79y3XOOPXnpbPxxoRsVeNNmNoveT9F114jR6TUvotcH/xcMeI+FRVk1/SGrzObutsYURsQetZ+IuLINty/L/TGjqPiog1atXSkTORNULnmm00vQ34Z0Xd1V9X9Wt/qo3Xb/kZWSwi/JtT6qf84ZekCsUUGi33Yj6ZUvpbe31SSi8BDxUPd2rnLGNlv/dpPZNxYER8PiLmiohBnbyE7VFaA+8JEbFsRAwpjlNdy6XF+tsRcV5ErBgRn4iINSLiQvKE8M9S23m0XuJ3aUQcGxHLR8RCEbEx8AD57Gm9kVKPJv+xOwC4NiIujoj1ImKxopZPR8Q2EfEL4EVgnQ5+H9pyJXmQlCG0TmMzCbijRvuNyQMD/SIito6IFYqalilCwT20TpXTkX9G9CnF1B97kD9vKwKPRcQhEbFq8X1aLCKGR8ToiLgP+GsXXmYs0PJzd0FEnBYRK0fEghGxVkTcSL60GWBsSumJtg8zW/2wWA+iaj7OlNLL5J8ryJ+vuyNi/eJnZvmIOIrWKVFeBE5r4/j7keeWnYs8ndCJEbFa8T1fIiI2iIgzaZ3Ls642Qufd1aGzGFRqV/IVCUsAD0fEMcXvqQUjYpFie++IuJ0cTqsHDXq4WM8LHBsRi7f8fjOASv1ISsnFxcXFpVjIf4SlYjm6E/32rei3acXzxxfPTazRb1RFv+plWEW7B4rnxtU4zq9qHGNiVbu5yZeZ1nrNszpQc0ugbKv/zOJ70V69S5AvNaxVR+XyzW6+p3dUHe/8Om1372BNp/XQ521icbwHOti+5fV3b2Nf3e959dfX1dcp9n+T1jP09ZZJXfy+LEO+dLnesX8FDKnRf1xnvq81jjGs4rWO70D7+4q2U4Glq/YNIN/vWe/rmQisXOf4awAvtXOMCVV9VqjY99U2jrk4eYqjRA60a7bRZl3ylRHtvdfTgbmq+gYwoUb7O3viZ8jFxaX5F/+7JEkfVTnwT2fOYF1L62VzHb6sNqV0IXkgmweAN2k9U9lZewDHkM92Tqb1Ut3q15tCHqX1ZPLZjanF695Hnpvv4Lb6VR3jH+QpYX5KHlhkGvkyy1uBDVJKP+/AMV4Bvka+N+xa8mAxHxbHeoV8ieKRwKdTSre0d7x2VF8+W2/+wWvIlzaOIU+L8iL5e/QB+QzOOGDtlNJRtQ7QHxTvyafI9y/+FngDmAFMIc/BehX5LOSwLh7/ReCLwIHkqwfeJP98vQLcBGyZUhqZUprWrS+kZ/2wWA8BDq/ckVKalVI6kPwPravJwXEaOeT9mfxZXzWlVG9+24fJ/+w5lPw9mUT+nrxcHONEPj4FUF2p9Uzn09S4vDblS4ZXIA+SdR/5Z306+WfiOeBG8u+fxVK+aqOybyKf1f1p8RqVc7xK6ici/y6QJEmSJKlneYZTkiRJktQQBk5JkiRJUkMYOCVJkiRJDWHglCRJkiQ1hIFTkiRJktQQnZlYXG1YeOGF07Bhw8ouQ5IkSZJK8fDDD7+RUlqkrX0Gzm4aNmwYEyZMKLsMSZIkSSpFRDxfa5+X1EqSJEmSGsLAKUmSJElqCAOnJEmSJKkhDJySJEmSpIYwcEqSJEmSGsLAKUmSJElqCAOnJEmSJKkhDJySJEmSpIYwcEqSJEmSGsLAKUmSJElqCAOnJEmSJKkhDJySJEmSpIYwcEqSJEmSGsLAKUmSJElqCAOnJEmSJKkhDJySJEmSpIYwcEqSJEmSGsLAKUmSJElqCAOnJEmSJKkhDJySJEmSpIYwcEqSJEmSGsLAKUmSJElN7vnn4bnnyq6i8wyckiRJktTEZsyAHXeE1VeHBx4ou5rOMXBKkiRJUhM74QT44x9h/vlhtdXKrqZzDJySJEmS1KQeeABOOQUGDIDx42HBBcuuqHMMnJIkSZLUhN54A3baCVKCY46B9dYru6LOM3BKkiRJUpNJCfbaC15+GdZZB449tuyKusbAKUmSJElN5vzz4ZZbYIEF4IorYNCgsivqGgOnJEmSJDWRJ56AQw/N2xddBMsuW2493WHglCRJkqQm8f77sMMOMHUqjBoFI0aUXVH3GDglSZIkqUkcfDA89RR85jNw1lllV9N9Bk5JkiRJagLXXQcXXghDh8JVV8Fcc5VdUfcZOCVJkiSpZM8/D/vsk7fPPBNWW63cenqKgVOSJEmSSjRjRp5v8+234ZvfhNGjy66o5xg4JUmSJKlEJ50Ev/89LLUUXHIJRJRdUc8xcEqSJElSSR58EE4+OYfM8eNhoYXKrqhnGTglSZIkqQSTJsHOO8OsWfCDH8D665ddUc8zcEqSJEnSbJYS7L03vPQSrL02/PCHZVfUGAZOSZIkSZrNLrgAbroJ5p8ffvUrGDSo7Ioaw8ApSZIkSbPR3/4GBx+cty+6CIYNK7WchjJwSpIkSdJs8v77sOOOMHVqvqR2223LrqixDJySJEmSNJsceig8+SSsvDKMGVN2NY3XtIEzIk6NiFQs36/TbmREPBQR70TE5IiYEBHfiYi6X1tX+0mSJElSV9xwQ753c8gQuOoqmHvusitqvKYMVxHxJeBwILXT7jzgCmA48BBwD7AicC5wXUQM7Ml+kiRJktQVzz+fL6EFOOMM+Pzny61ndmm6wBkRQ4FxwGvAzXXabQOMBl4FVkspfSOltBXwaeApYCvggJ7qJ0mSJEldMX16vm/zrbfgG9+AA/pR2mi6wAmcCKwC7Ae8U6fdUcX6iJTSv1qeTCm9BuxfPDyyjUtku9pPkiRJkjrtuOPgj3+EpZeGceMgouyKZp+mClURsSZwKPCrlNKtddotDawBTAOurd6fUnoQ+A+wOLBWd/tJkiRJUlfcdRecfjoMGJDn21xoobIrmr2aJnBGxBzAZcCbwPfaaf6FYv1kSumDGm3+WtW2O/0kSZIkqVNeeQV22SVvn3gifO1r5dZThkFlF1DhFGAlYIeU0hvttF2uWD9fp80LVW2700+SJEmSOmzmTNhpJ3j9ddhwQzjyyLIrKkdTnOGMiK8ABwE3pZSu7kCXeYr1lDptJhfreXugnyRJkiR12Kmnwv33w6KLwvjxMLCfzoNReuCMiDmBS4F3yaPHdqhbsa47bUoP9vvoQSJGFfN2Tnj99de7cyhJkiRJfcyDD8Lxx+fBga64AhZfvOyKylN64AROJc+BeUhK6ZUO9nmvWM9Tp03Lvvcqnutqv49IKV2YUhqeUhq+yCKL1C1UkiRJUv/x+uswciTMmgVHHQUbbVR2ReVqhns4twJmAbtFxG5V+1Yu1vtHxDeAZ1JKewMTi+eXrXPcZYr1xIrnutpPkiRJkuqaNQt23x1efhnWWQdOOKHsisrXDIET8pnW9ers/1SxLFA8frRYrxoRc9YYcfZLVW2700+SJEmS6jrrLPj1r2HBBeHKK2FQs6StEpV+SW1KaVhKKdpayNOkABxWPLd60edF4BFgCLBt9TEjYj1gaeBV4I8Vr9WlfpIkSZJUz5//3DoS7bhxsMwydZv3G6UHzm44rVj/KCJWaHkyIhYFzi8enp5SmtVD/SRJkiTpY95+G3bYAWbMgIMOgi23LLui5tFrT/KmlK6LiLHA/sDfIuJeYDqwITAfcBNwbk/1kyRJkqRqKcHee8PEiTB8OPzoR2VX1Fx6beAESCmNjojfAd8h3wM6EHgauAQYW+ssZVf7SZIkSVKlsWPh+uth3nnhqqtgyJCyK2oukVK3pqTs94YPH54mTJhQdhmSJEmSZrPHHoO11oKpU+Hqq2G77cquqBwR8XBKaXhb+3rzPZySJEmSVIrJk2H77XPYHDWq/4bN9hg4JUmSJKkTUoL994d//hM++1kYM6bsipqXgVOSJEmSOuGyy2D8eJhrLrjmGphzzrIral4GTkmSJEnqoCefhNGj8/Z558FnPlNuPc3OwClJkiRJHTB5Mmy7LXzwAeyyC+y2W9kVNT8DpyRJkiS1I6V8ZvOpp2CVVfJ0KBFlV9X8DJySJEmS1I5LLoFf/jLft3nttTD33GVX1DsYOCVJkiSpjieegAMOyNtjx+YznOoYA6ckSZIk1fDee/m+zQ8/hD33hF13Lbui3sXAKUmSJEltSAlGjcrzbX7uc3DOOWVX1PsYOCVJkiSpDT//OVx1FcwzT75vc665yq6o9zFwSpIkSVKVRx6B730vb194Iay0Urn19FYGTkmSJEmq8M47+b7NadNg331hxx3Lrqj3MnBKkiRJUiEl2GsveO45WH11GDOm7Ip6NwOnJEmSJBXOPReuvx7mnTfftznHHGVX1LsZOCVJkiQJ+Otf4dBD8/bFF8MKK5RbT19g4JQkSZLU7731Fmy3HUyfDgcckO/hVPcZOCVJkiT1aynBHnvAxIkwfDiccUbZFfUdBk5JkiRJ/dpZZ8HNN8P888M118DQoWVX1HcYOCVJkiT1W3/6ExxxRN6+9FJYbrly6+lrDJySJEmS+qVJk/J9mzNmwMEHw1ZblV1R32PglCRJktTvzJoFu+4KL74Ia64Jp59edkV9k4FTkiRJUr9z6qnw61/DggvC1VfDkCFlV9Q3GTglSZIk9Sv33gvHHQcRcMUVsOyyZVfUdxk4JUmSJPUbL70EO+6Yp0I59ljYbLOyK+rbDJySJEmS+oVp0/IgQW+8AZtsks9yqrEMnJIkSZL6hcMPhz/+EZZeOl9KO3Bg2RX1fQZOSZIkSX3eNdfA2WfD4MFw7bWw8MJlV9Q/GDglSZIk9WlPPw177ZW3zzwT1lqr3Hr6EwOnJEmSpD5ryhQYMQImT4YddoADDii7ov7FwClJkiSpT0oJ9t0XnnwSVl4ZLrooT4Wi2cfAKUmSJKlPuuCCPDjQ3HPD9dfDPPOUXVH/Y+CUJEmS1Of89a9w0EF5+6KLYJVVyq2nvzJwSpIkSepTJk3K921Omwbf+Q7suGPZFfVfBk5JkiRJfcasWbDLLvDCC/DlL+dRaVUeA6ckSZKkPuOUU+COO2ChhfJ8m0OHll1R/2bglCRJktQn3HMP/PCHeSTaK66AT36y7Ipk4JQkSZLU6734IowcmadCOe442HTTsisSGDglSZIk9XLTpsF228Ebb8Amm8Cxx5ZdkVoYOCVJkiT1aoccAn/6EyyzTL6UduDAsitSCwOnJEmSpF7rl7+E886DIUPguutg4YXLrkiVDJySJEmSeqXHHoNRo/L2OefkaVDUXAyckiRJknqdN9+ErbeGDz+EvfaCffYpuyK1xcApSZIkqVeZNQt23hn+/W8YPhzOPTdPhaLmY+CUJEmS1KuccALccQcstFC+b3OOOcquSLUYOCVJkiT1GrfeCieeCAMGwFVXwbLLll2R6jFwSpIkSeoVnnkGdtklb59yCmy0Ubn1qH0GTkmSJElNb8oU2GoreOedvD7iiLIrUkcYOCVJkiQ1tZTyKLT/7//BSivBuHEOEtRbGDglSZIkNbWf/QyuvBLmnhtuuAHmm6/sitRRBk5JkiRJTeuhh+D738/b48bBKquUWo46ycApSZIkqSm9/DJsuy3MmAGHHQYjRpRdkTrLwClJkiSp6UyblsPma6/B178Op55adkXqCgOnJEmSpKZzyCHwhz/A0kvn+TYHDSq7InWFgVOSJElSU/nlL+G882DIELj+elh00bIrUlcZOCVJkiQ1jcceg1Gj8vY558CXv1xuPeoeA6ckSZKkpvDmm7D11vDhh7DnnnnuTfVuBk5JkiRJpZsxA3bYAf79b1hjjXxJbUTZVam7DJySJEmSSnfUUXDPPbDIInDjjTDHHGVXpJ5g4JQkSZJUqiuvhDPOyCPRXncdLLNM2RWppxg4JUmSJJXmscdgr73y9pgxsO665dajnmXglCRJklSKN96Ab38bPvgA9tgDRo8uuyL1NAOnJEmSpNluxgzYfnt4/vk89cn55ztIUF9k4JQkSZI02x1+OPzmN7DYYnDDDQ4S1FcZOCVJkiTNVuPHw1lntQ4StNRSZVekRjFwSpIkSZptHnkE9tknb59zDnz1q+XWo8YycEqSJEmaLf773zxI0Icfwt57w777ll2RGs3AKUmSJKnhpk+H7baDF1+EtdaCc891kKD+wMApSZIkqeG+/3148EFYfHG4/noYOrTsijQ7GDglSZIkNdRll8HPfgaDB+cRaZdcsuyKNLsYOCVJkiQ1zF//2nqv5nnnwdprl1uPZi8DpyRJkqSGeO012HprmDoV9tuvdXRa9R8GTkmSJEk9bto02HZbeOklWGcdOPvssitSGQyckiRJknrcwQfDQw/l+zWvuw6GDCm7IpXBwClJkiSpR11wAZx/fg6ZN9yQR6ZV/2TglCRJktRjHnwQDjwwb190Eay5Zrn1qFwGTkmSJEk94t//hm22gRkz8rybu+5adkUqm4FTkiRJUrdNngzf+hZMmgSbbQann152RWoGBk5JkiRJ3TJrVj6b+be/wUorwZVXwsCBZVelZknf48oAACAASURBVGDglCRJktQtJ5wAN94I888Pt9wCCyxQdkVqFgZOSZIkSV127bVw4okwYABcfTWsuGLZFamZGDglSZIkdcmjj8Juu+XtM86ATTcttx41HwOnJEmSpE577bU8SNAHH8Duu8NBB5VdkZqRgVOSJElSp0yblqc/efFFWHttuOACiCi7KjUjA6ckSZKkDksJRo+G3/8ell4abrgBhg4tuyo1KwOnJEmSpA4791y4+GKYYw646SZYfPGyK1IzM3BKkiRJ6pB774WDD87bl14Ka6xRbj1qfgZOSZIkSe165hnYbjuYOROOPhp22KHsitQbGDglSZIk1fXuu/DNb8Jbb8GWW8JJJ5VdkXoLA6ckSZKkmmbOhJEj4amnYNVVYfx4GGCKUAf5UZEkSZJU01FHwe23w4ILws03w3zzlV2RehMDpyRJkqQ2XXop/OQnMGgQXHstLL982RWptzFwSpIkSfqYBx+EfffN2+edBxtsUG496p0MnJIkSZI+4tlnYeutYfp0OOggGDWq7IrUWxk4JUmSJP3P22/DN74Bb74Jm28OZ5xRdkXqzQyckiRJkgCYMQO23x6efho++1m48koYOLDsqtSbGTglSZIkAXDwwXD33bDIInDrrY5Iq+4zcEqSJEni/PPh3HNhyBC48UYYNqzsitQXGDglSZKkfu6ee+C7383bF18M66xTbj3qOwyckiRJUj/29NOw7bYwcyYcfTTsvHPZFakvaZrAGREHRsQ1EfFUREyKiOkR8XpE3BsRO0dE1Ok7MiIeioh3ImJyREyIiO9ERN2vr6v9JEmSpL5g0qQ8Iu0778A228BJJ5VdkfqaQWUXUOEIYFHg/wF/AKYAywIbABsCIyJi65TSrMpOEXEeMBr4ELgPmF60PxfYMCK2TSnNrH6xrvaTJEmS+oJp0/Jcm88+C1/8Ilx2GQzwtIt6WDMFzh2AR1NKUyqfjIhVyYHwW8BuwKUV+7Yhh8ZXgXVTSv8qnl8MuB/YCjgAOLvqmF3qJ0mSJPUFKcH++8NvfwtLLgm33AJzz112VeqLmuZ/GCml31WHzeL5J4HziocbV+0+qlgf0RIaiz6vAfsXD49s4xLZrvaTJEmSer0zz4RLLoE558xhc6mlyq5IfVVvCVQzivWHLU9ExNLAGsA04NrqDimlB4H/AIsDa3W3nyRJktQX3HILHH543r78clhjjXLrUd/W9IEzIpYD9ise3lqx6wvF+smU0gc1uv+1qm13+kmSJEm92uOPw8iR+ZLak0+GESPKrkh9XTPdwwlAROwBrAcMBpYGvkIOxqellG6saLpcsX6+zuFeqGrbnX6SJElSr/Xqq7DlljBlSp765Oijy65I/UHTBU5gHfLgQC1mAMcCP61qN0+x/th9nxUmF+t5e6CfJEmS1Cu9/34Omy++CGuvDRddBLUnHZR6TtNdUptS2julFMBcwKrAGOB44E8RsWRF05YfkdTJl+hqv9YDRIwq5uyc8Prrr3f1MJIkSVLDzZyZz2hOmACf+hTcfDPMMUfZVam/aLrA2SKl9EFK6e8ppcPIo8p+njxHZov3ivU8H+vcqmXfexXPdbVfZW0XppSGp5SGL7LIInUOI0mSJJXriCPgxhthgQXg9tvBP181OzVt4KzSMvfmlhExuNieWKyXrdNvmaq23eknSZIk9Spjx+YpUAYPzqFz5ZXLrkj9TW8JnG+T7+UcBCxYPPdosV41Iuas0e9LVW2700+SJEnqNe64Aw44IG9fdBGsv36p5aif6i2Bc11y2HwbeAMgpfQi8AgwBNi2ukNErEce5fZV4I8tz3e1nyRJktRbPP44bLcdzJoFxxwDu+3Wfh+pEZoicEbE1yJip4gY2sa+dYCLi4cXp5RmVuw+rVj/KCJWqOizKHB+8fD0lNKsqsN2tZ8kSZLU1F5+Gb7xDZg8GXbcEU48seyK1J81y7Qoy5Pv0zw3Ih4hn12ct3h+laLN7eTpUf4npXRdRIwF9gf+FhH3AtOBDYH5gJv46EBD3eonSZIkNbMpU/L0Jy+9BOusA5dc4vQnKlezBM4HgZOArwErAl8hT1/yKnA9MD6ldFNbHVNKoyPid8B3gPWAgcDTwCXA2FpnKbvaT5IkSWpGM2fCyJHwyCOw/PJw001Of6LyRUpdno5SwPDhw9OECRPKLkOSJEn93MEHw5gx8IlPwJ/+BCuuWHZF6i8i4uGU0vC29jXFPZySJEmSuu7cc3PYHDw4n9k0bKpZGDglSZKkXuz22+F738vbF18M665bbj1SJQOnJEmS1Es99hhsv32e/uSHP4Rddim7IumjDJySJElSL/Sf/+TpT6ZMgZ13zoFTajYGTkmSJKmXefdd2HzzHDq/9jX4xS+c/kTNycApSZIk9SLTp8M228ATT8BKK+VBgoYOLbsqqW0GTkmSJKmXSAn22QfuvRcWXRTuuAMWXLDsqqTaDJySJElSL3H88XDZZTDXXHl02uWWK7siqT4DpyRJktQLXHIJnHgiDBgAV18Nw4eXXZHUPgOnJEmS1OTuugtGjcrb55+fR6eVegMDpyRJktTEHn0URoyAmTPhyCNh333LrkjqOAOnJEmS1KReeAG22AImT4aRI+GUU8quSOocA6ckSZLUhN5+O8+1+corsP76+R7OAf71rl7Gj6wkSZLUZKZOha22giefhFVWgRtvdK5N9U4GTkmSJKmJpAR77QUPPABLLJHn2lxggbKrkrrGwClJkiQ1kR/8AK64AuaZJ8+1+clPll2R1HUGTkmSJKlJ/PzncNppMHAgXHcdfOELZVckdY+BU5IkSWoCt98Oo0fn7QsvhE03LbceqScYOCVJkqSS/fnPsO22MGsWHHss7Lln2RVJPcPAKUmSJJXoH//Ic21+8EEOmiecUHZFUs8xcEqSJEklefnlfOnspEk5dP785xBRdlVSzzFwSpIkSSV45x3YfHN4/nlYc024+moYNKjsqqSeZeCUJEmSZrOpU2HrreHxx2HFFeG222DuucuuSup5Bk5JkiRpNpo1C3bbDX7zG1h8cbjrLlh44bKrkhrDwClJkiTNJinBIYfky2fnnRfuuAOGDSu7KqlxDJySJEnSbHLGGXD22TB4MNx0E6y+etkVSY1l4JQkSZJmg1/+Eg4/vHV7gw3KrUeaHQyckiRJUoPdfXeeYxNgzBjYfvty65FmFwOnJEmS1EATJuQRaWfMgMMOg+99r+yKpNnHwClJkiQ1yDPP5Lk2p0yBnXeG008vuyJp9jJwSpIkSQ3w2muw2Wbw+uuwySZw8cUwwL++1c/4kZckSZJ62Lvv5jObzz4La6wB110HQ4aUXZU0+xk4JUmSpB704Yfw7W/DI4/A8svD7bfnOTel/sjAKUmSJPWQmTNhp53g/vth8cXz6LSLLVZ2VVJ5DJySJElSD0gJ9tsPbrgBFlgA7roLPvWpsquSymXglCRJknrAD34Av/gFzDkn3HYbrLZa2RVJ5TNwSpIkSd3005/CaafBwIFw7bWwzjplVyQ1h0Gd7RARiwCrA4sBCwBvAf8FHk0pvdGz5UmSJEnN7fLL4dBD8/a4cbDFFqWWIzWVDgXOiFga2Bf4FrBqnXZPAjcBF6aUXuqRCiVJkqQmdeutsOeeefuss2DnncutR2o2dQNnRCwPnAZ8u6LtW8BTwJvAu8B8wELAysBni+XIiLgROCql9FxjSpckSZLK89BDsN12eWTao4+Ggw4quyKp+dQMnBHxY+C7wBBgAnAZcG9K6R91+qwMbAzsBmwLfCsifpZSOrxHq5YkSZJK9PjjsOWWec7NUaPg5JPLrkhqTvUGDToUuBVYLaX05ZTSefXCJkBK6emU0jkppeHA54HbgEN6rlxJkiSpXM89B5ttBu+8AyNGwPnnQ0TZVUnNqd4ltcNTSo929cAppb8BIyLiC109hiRJktRMXn0VNt44rzfcEMaPzyPTSmpbzTOc3QmbjTiOJEmSVKa334ZNN81nOIcPhxtvhKFDy65Kam7OwylJkiS14/334ZvfhCeegJVWgl//Guadt+yqpOZn4JQkSZLqmDYNttkmj0q71FJw992wyCJlVyX1Dh2ah7NFRHwCGA18HVgSmKNG05RSWr6btUmSJEmlmjkzz615552w8MJw773wyU+WXZXUe3Q4cEbECsCDwOJAe+Nwpe4UJUmSJJVt1qw85cm118J888Fdd8HKK5ddldS7dOYM55nAEsBDwFnAv4DJjShKkiRJKlNKcOihcMklMOeccPvt8MUvll2V1Pt0JnCuD0wENk4pTWtINZIkSVITOPFEGDMGBg/Oo9F+9atlVyT1Tp0ZNCgBfzFsSpIkqS8bMwaOPx4GDIArr8xToUjqms4EzsfI929KkiRJfdIll8DBB+ftiy/Oo9NK6rrOBM4zgK9GxFcaVYwkSZJUlmuvhX32ydtnnw27715qOVKf0OF7OFNKt0XEwcDtEXEucBfwEjCrRvsXeqZESZIkqbHuuAN22imPTHviifDd75ZdkdQ3dGoeTuBR4DXg6GKpJXXh2JIkSdJs99BD+dLZ6dPzyLTHHFN2RVLf0Zl5ONcH7gSGFE9NwmlRJEmS1Is9/DBssQV88AHsvTf85CcQ7c04L6nDOnMW8iRy2PwxcHpK6e3GlCRJkiQ13t//nkegfe892H57uOACw6bU0zoTOFcHHk4pHdmoYiRJkqTZ4dlnYeONYdIk2HxzuPxyGDiw7Kqkvqczo9R+APyrUYVIkiRJs8MLL8CGG8LLL8N668F118GQIe33k9R5nQmcDwGrNqoQSZIkqdFefjmHzeefh7XWgltvhTnnLLsqqe/qTOA8Flg+Ir7XqGIkSZKkRvnvf2GjjeCZZ+CLX8xTocw7b9lVSX1bZ+7hHA5cCvw0IkbQ/jycl3e/PEmSJKn73nwz37P51FPw2c/C3XfDAguUXZXU93UmcI4jz68ZwDrAV9ppb+CUJElS6d55J49G+8QTsOKKcO+9sNBCZVcl9Q+dCZyXkwOnJEmS1CtMnpxHoZ0wAZZbDu67DxZbrOyqpP6jw4EzpbR7A+uQJEmSetQHH8A3vwl/+AMsswz85jew9NJlVyX1L50ZNEiSJEnqFaZOha23hvvvhyWWyGc2hw0ruyqp/zFwSpIkqU+ZPh223x7uvBMWXjjfs/npT5ddldQ/1QycETEqIgZ25+ARMTAiRnXnGJIkSVJHzZwJu+wCN98Mn/hEDpurrFJ2VVL/Ve8M5wXA3yNit4jo1HS4ETFnROwOPAWM7UZ9kiRJUofMmgV77glXX53n17zrLvj858uuSurf6gXOHYE5gEuAVyPiFxGxY0QMa6txRCwXESMj4hLgVeBiYAiwQ8+WLEmSJH1USjB6NFx+Ocw1F/z61/ClL5VdlaSao9SmlK6OiJuBQ4DRwJ7AHgARMRV4E3gXmA9YiBwuIc/T+RJwKnB2SunDhlUvSZKkfi8l+N734Oc/hznmgFtvha9+teyqJEE706IUYfHUiPgRsDXwbWBdYClgyWJp8SJwP3ATcEtKaVZDKpYkSZIKKcEhh8A558CQIXDjjbDBBmVXJalFh+bhTCnNBK4tFiJiYWBRYH7gbeC/KaVJjSpSkiRJqpYSHHYYjBkDgwfDDTfAZpuVXZWkSh0KnNVSSm8Ab/RwLZIkSVKHpARHHglnnpnD5vXXwxZblF2VpGrOwylJkqReJSX4wQ/gxz+GQYPg2mthyy3LrkpSWwyckiRJ6jVSguOOg9NOg4ED8xQo3/pW2VVJqsXAKUmSpF7jhBPg5JNz2LzqKth667IrklSPgVOSJEm9wkkn5cA5YABccQWMGFF2RZLaY+CUJElS0zv11Hwp7YABMH48bL992RVJ6ggDpyRJkpraj36UBwmKgMsugx13LLsiSR1l4JQkSVLTOuOMPP1JBIwbBzvvXHZFkjqjw4EzInaNiK90oN1aEbFr98qSJElSf3fWWXDYYTlsXnwx7OpfmFKv05kznOOAvTvQbi/g0i5VI0mSJAFjxsAhh+Ttiy6CPfYotx5JXdOIS2qjAceUJElSP3HGGXDwwXn7ggtgr73KrUdS1zUicC4NTG7AcSVJktTH/ehH+TJagAsvhH33LbceSd0zqN7ONu7FXKHO/ZmDgM8AGwJ/7YHaJEmS1I+ccgocc0y+Z/MXv4A99yy7IkndVTdwku/bTBWP1ymWWgKYBZzRvbIkSZLUn5xwAhx/fA6bl14Ku+1WdkWSekJ7gfNyWgPnbsCzwO9rtJ0G/Ae4OaX0eM+UJ0mSpL4sJfjhD+Gkk2DAgDzPplOfSH1H3cCZUtq9ZTsidgN+l1Ly4gZJkiR1W0r5EtpTT81hc/x42HHHsquS1JPaO8NZaTkcDEiSJEk9ICU46qg8SNDAgfCrX8F225VdlaSe1uHAmVJ6vpGFSJIkqX9IKY9Ee+aZMGgQXHUVbLNN2VVJaoTOnOEEICLmAIYDSwJz1GqXUrq8G3VJkiSpD0opz7F59tkweDBccw18+9tlVyWpUToVOCPiYOA4YL4ONDdwSpIk6X9Sgu9+F849F4YMgeuugy23LLsqSY3U4cAZEXsCZxYPnwKeBt5tRFGSJEnqW2bNggMOgLFjc9i88UbYfPOyq5LUaJ05w/ld8hQpu6SUftWgeiRJktTHzJoF++0HF10EQ4fCTTfBZpuVXZWk2aEzgXNF4A+GTUmSJHXUjBmwxx55ypM55oCbb4ZNNim7KkmzS2cC5/vAC40qRJIkSX3LtGkwciRcfz3MPTfcdhusv37ZVUmanToTOP8AfLZRhUiSJKnv+PBDGDECbr8d5p8f7rgD1l677KokzW4DOtH2BGDliNitUcVIkiSp95syJY8+e/vtsNBC8JvfGDal/qrmGc6IWLeNp38KXBIRmwO3ky+xndVW/5TSb3ukQkmSJPUa774LW2wBv/sdLLYY3HsvfNZr5KR+q94ltQ+QR6WtFsCIYqkltXPsjx4wYjCwLrA5sA6wLLAQ8DrwR+DclNIDdfqPBPYHVgMGkqdsuRQYm1JqMxB3p58kSZI+7s03YdNNYcIEWHppuO8+WHHFsquSVKZ6ofC3tB04G2E94J5i+1XgYWAKsAqwDbBNRJyUUjquumNEnAeMBj4E7gOmAxsC5wIbRsS2KaWZPdVPkiRJH/ff/8LGG8MTT8Byy+WwudxyZVclqWw1A2dKaf3ZWMcs4Hrg7JTSQ5U7ImJ74Arg2Ii4P6V0f8W+bcih8VVg3ZTSv4rnFwPuB7YCDgDOrjpml/pJkiTp4/7zH9hoI3j66XxG87778hlOSerMoEENk1L6TUppRHXYLPZdDYwrHu5ctfuoYn1ES2gs+rxGvlQW4MiIqP46u9pPkiRJFSZOhHXXzWHzc5+D3/7WsCmpVW8JVI8W6//9+oqIpYE1gGnAtdUdUkoPAv8BFgfW6m4/SZIkfdS//pXD5nPPwRprwP3354GCJKlFZwb2aWvU2rZMA95IKT3TtZLa9Oli/UrFc18o1k+mlD6o0e+vwFJF2z90s58kSZIKTz6ZL6N99VX4ylfg17/O821KUqUOB05qj1rbpoh4F7gMODal9F4n66o8zuLA7sXD6yt2tdyG/nyd7i9Ute1OP0mSJJFHod1sM5g0Cb7+dbjlFphnnrKrktSMOnNJ7W/JU5REsbwNPAE8BrxVPAfwZ+A5YB7gQOChiJirK8VFxCBgPDA/cF9K6daK3S2/1qbUOcTkYj1vD/STJEnq9x58EDbYIIfN//s/uP12w6ak2joTODcr1n8HNk8pLZRS+kJKaY2U0sLA/wFPks+Cfo58Gewfiu3vdrG+C8hTlbzIxwcMagm4nZ26pav9Wg8QMSoiJkTEhNdff72rh5EkSepVbrstn9l87z3Yfnu46SaYc86yq5LUzDoTOI8hh8cNUkp3Vu9MKd0FbAx8FjgupTQRGAlMJc+l2SkRcTawF3nqkg1TSq9WNWm5TLfe/9Ra9lVe0tvVfv+TUrowpTQ8pTR8kUUWqXMYSZKkvuFXv4KttoIPP4RRo+CKK2DIkLKrktTsOhM4twfuTyn9t1aDYlqR+4HtiscvAo8AK3amqIg4k3xW9HVy2PxXG80mFutl6xxqmaq23eknSZLUL40dCzvvDDNmwBFHwAUXwMCBZVclqTfoTOBcmny2sj1TySO8tngRGNrRF4mIHwOHAJOAjVNKf6/RtGWqlFUjotbFHF+qatudfpIkSf1KSnDaaTB6dOv26adDRPt9JQk6FzjfANatE9Io9q1LDostPkEeYKhdEXE6cBh5EKKNU0qP12pbcfZ0CLBtG8dajxySXyUPdtStfpIkSf1JSvls5tFH54A5diwceWTZVUnqbToTOG8FFgOuiYhlqncWz10NLArcUrFrZfKotXVFxEnAEeRwunFKqSNnF08r1j+KiBUqjrUocH7x8PSU0qwe6idJktTnzZwJ++4LP/kJDBqU79fcb7+yq5LUG0VKHRusNSIWAf5CvvdxOvns3/Pk0V6XBb4CDC6e+3JK6fWIWAP4K3BySum4Osf+JnBz8XACebTbtjydUjq9qu/5wP7Ah8C9RW0bAvMBNwEjUkoz23jNLvWrNnz48DRhwoT2mkmSJPUK06bBLrvANdfAHHPAddfBFluUXZWkZhYRD6eUhre1b1BHD1IEyK8AY4EtyZfOfqQJcBuwf0rp9aLPwxExuAPBbcGK7eHF0pYHgY8EzpTS6Ij4HfAdYD1gIPA0cAkwttZZyq72kyRJ6qvefx+22QbuvBPmnTdPg7Ju9V98ktQJHT7D+ZFOEZ8kB86WwYFeBh4qpkLpVzzDKUmS+oK334ZvfAN+/3tYeOEcOtdYo+yqJPUGPXKGs1JK6QVgfLeqkiRJUlN47TXYbDN47DFYemm45x5YeeWyq5LUF3QpcEqSJKlveO452GQTePZZWGEFuPdeWLbebOWS1Ak1A2dx2SzAf1JKMysed0hxFlSSJElN6rHH8pnN116DL34R7rgDFl207Kok9SX1znBOBGYBqwD/LB539IbP1M6xJUmSVKL774dvfQveew823BBuuAHmm6/sqiT1NfVC4Qvk4Di96rEkSZJ6seuvh5Ej8xQo220Hl18OQ4eWXZWkvqhm4EwpDav3WJIkSb3PBRfA6NGQEhxwAJx9NgwYUHZVkvoqf71IkiT1AynBCSfA/vvn7ZNOgp/9zLApqbG8z1KSJKmPmzkTDjwQxo7NAfOCC2Cffcqu6v+3d99hdlX1/sff3zQCJICETpAuXVqQAKGGrmDoRZByASkCXjriFX+A0ptUkSYiIP2KFKUISCdSxIRIh0AIhJYQTINZvz/WmZvCzOTMzDmzzznzfj3PfnbO3vvsfGexgHyy9l5LUnfQ7sAZEcsBPwLWAxYE/jeldFzp3GDg28DNKaXPKlmoJEmS2m/KFNhrL7j11vye5o03wg47FF2VpO6iXYEzIv4LuAToUzqUgAVmuGRB4DLyREPXVKJASZIkdcyECTBsWJ6Rdp554E9/go03LroqSd1J2U/tR8QGwG+AycCxwLpAzHLZfcAEYPtKFShJkqT2++AD2GSTHDYXWQQefdSwKanrtWeE8zjyiOY2KaUnASJmzpsppWkR8W9gpYpVKEmSpHZ57TXYemt4/XVYbjn4619h6aWLrkpSd9SeecnWA55pDpttGA0s2vGSJEmS1FHPPAPrr5/D5tprw+OPGzYlFac9gXNe4N0yruuDs99KkiR1ubvuyo/RjhsHW22VH6ddaKGiq5LUnbUncH4IlPP3YysA73WsHEmSJHXEb36TJwiaNAn22y+Hz/79i65KUnfXnsD5OLBWRAxq7YKI2AL4FvBwJ+uSJElSGVKCk06Cgw+GpiY4+WS46iro3bvoyiSpfYHzfPKstLdHxJYRMdN3I2Ij4GrgS+CiypUoSZKklkydCvvsA7/6FfTsCVdeCb/4BcSs6whIUkHKftcypfR0RBwHnA3cS17+JAHDIuK75PU4AzgqpfRSNYqVJElSNmEC7LQTPPAAzD033HILbLNN0VVJ0szaM8JJSulcYFtgODAPOWDOBywI/AsYllK6oNJFSpIkabr33oMNN8xhc+GF4ZFHDJuSalO7Z5NNKd0H3BcRA8iTCPUERqeUxlS6OEmSJM1sxIgcLkePhhVWgHvvddkTSbWrw8uXpJQ+Bj6uYC2SJElqw8MP55lox4/Pa23+6U8wYEDRVUlS69r1SK0kSZKKcdNNeW3N8eNhxx3z47SGTUm1rtURzoj4YWdunFK6rjPflyRJUl725Oyz4fjj8+cjjoDzzsuz0kpSrWvrkdprybPQdpSBU5IkqROmTYPDDoPf/jZ/PuccOOoolz2RVD/aCpyP0nrg3Bj4ABhV8YokSZLE+PGwyy5w//3Qty9cf31eBkWS6kmrgTOltElr5yKiCbg3pbR/NYqSJEnqzt56C777XRg5EhZaKE8OtO66RVclSe3X4VlqJUmSVHnPPAPbbQcffggrrwx33w1LLVV0VZLUMc5SK0mSVCNuuw023jiHzc03h8cfN2xKqm8GTkmSpII1z0S7884weTIccADccw/MN1/RlUlS5/hIrSRJUoGmTYMf/xiuuCJ/PvNMOPZYZ6KV1BgMnJIkSQWZdSba3/8+j3JKUqMwcEqSJBXg7bfzTLQjRjgTraTG1WrgjIiNZvPdRdq6JqX0aIerkiRJamDPPptnov3gA1hppTwT7dJLF12VJFVeWyOcDwOplXMJ2Kq0tXbe0VNJkqRZ3HQT7Ldfnhxo6FC49VYnB5LUuNoKhe/QeuCUJElSOzQ1wS9+Aaeemj8feCBccgn07l1oWZJUVa0GzpTSUl1YhyRJUsP64gvYZ5+8zmaPHnD++XD44c5EK6nx+dirJElSFb37Lmy/PTz/PMw7L/zxj7BVay8lSVKDMXBKkiRVydNPw7BhMHYsLLss/PnPsOKKRVclSV2nR9EFSJIkNaIbwPb6cAAAIABJREFUboCNN85hc9NNc/g0bErqbgyckiRJFdTUBD/7GfzgBzBlChx8MPzlLzBgQNGVSVLX85FaSZKkCvniC9h7b7jjjjw50IUXwmGHOTmQpO7LwClJklQBo0fnyYFeeCFPDnTzzbDllkVXJUnFMnBKkiR10lNP5cmBPvgAll8e7roLVlih6KokqXi+wylJktQJ11yTJwf64AMYOjSHT8OmJGUGTkmSpA6YNg2OOAL23x+mTs3vat57L8w/f9GVSVLt8JFaSZKkdho3DnbdFR5+GPr0gUsugQMOKLoqSao9Bk5JkqR2eP75/L7mO+/AoovCbbfBeusVXZUk1SYfqZUkSSrTjTfCBhvksLnuujB8uGFTktpi4JQkSZqNr76C44+HPfeESZPye5uPPAKLLVZ0ZZJU23ykVpIkqQ2ffgq77w5//Sv06gUXXACHHgoRRVcmSbXPwClJktSKESPg+9+H11+HBRaAW2/NS6BIksrjI7WSJEktuOMOGDw4h80118zvaxo2Jal9DJySJEkzaGqCk0+GHXeEiRNhjz3gscdgySWLrkyS6o+P1EqSJJV88gnsvTfccw/06AFnnglHH+37mpLUUQZOSZIk8vqaO+0Eb74J88+fl0DZcsuiq5Kk+uYjtZIkqdv73e9g/fVz2Fx7bXjuOcOmJFWCgVOSJHVbU6bkJU723RcmT4YDDvB9TUmqJB+plSRJ3dK778LOO8PTT8Mcc8DFF+fAKUmqHAOnJEnqdh56CHbfHcaNg29+E267DQYNKroqSWo8PlIrSZK6jZTgrLNgiy1y2NxiC/jHPwybklQtBk5JktQtTJiQH6E9/vi81uZJJ8G998ICCxRdmSQ1Lh+plSRJDW/kSNhhB3jlFZhnHvj972H77YuuSpIanyOckiSpod14I3znOzlsrroqDB9u2JSkrmLglCRJDWnyZDjkENhzT/jii7x/6ilYfvmiK5Ok7sNHaiVJUsN57TXYdVd4/nno0wd+/Ws46CCIKLoySepeDJySJKmh3HYb7L9/niRomWXglltgrbWKrkqSuicfqZUkSQ1h6lQ48sg8E+2ECbDjjvDcc4ZNSSqSI5ySJKnuvfUW7LYbPPMM9O4N55wDhx/uI7SSVDQDpyRJqmt/+hPssw989hksuSTcfHOelVaSVDwfqZUkSXVp2jQ47jj4/vdz2Pze9/IjtIZNSaodjnBKkqS68+67+RHaJ56Anj3h9NPh6KOhh3+VLkk1xcApSZLqyr33wg9/CB99BIsvDjfdBEOGFF2VJKkl/j2gJEmqC1On5lHMbbfNYXPLLfM6m4ZNSapdjnBKkqSa99prsMceMHx4foT21FPh+ON9hFaSap2BU5Ik1bQ//AEOPhgmTsyz0N54I6y3XtFVSZLK4d8LSpKkmjRxIuy7L+y1V/71LrvACy8YNiWpnjjCKUmSas4LL+RZaF95Bfr2hQsvhAMPhIiiK5MktYcjnJIkqWakBBddBOuum8PmKqvk9zYPOsiwKUn1yMApSZJqwscfw7BhcMQReUbagw+GZ5/NoVOSVJ98pFaSJBXukUfgBz+A996D+eaDK6+EnXYquipJUmc5wilJkgozbRr8/Oew2WY5bK6/fn5/07ApSY3BEU5JklSI117LM9A+/XR+P/Okk+AXv4Be/ulEkhqG/0mXJEldKiW4+mo48kj44gsYOBCuuw423bToyiRJleYjtZIkqct8/DHsvDMccEAOm7vtBv/8p2FTkhqVI5ySJKlL3H8/7LMPvP8+9O8Pl1ySH6l1uRNJalyOcEqSpKqaPBn++79hyy1z2NxgA3jxRdh7b8OmJDU6A6ckSaqal16C73wHLrggTwZ02mnw8MOw9NJFVyZJ6go+UitJkiquqQl+/Ws44QSYMgWWXx7+8AdYZ52iK5MkdSUDpyRJqqgxY2DfffM7mwAHHgjnnQf9+hValiSpAAZOSZJUMTfdBIceCp9+CgMGwJVXwrBhRVclSSqKgVOSJHXaxx/noHnzzfnzVlvBNdfAoosWW5ckqVhOGiRJkjrl7rth1VVz2Jx7brj8crj3XsOmJMkRTkmS1EETJuTlTq6+On/ecEO49lpYZplCy5Ik1RBHOCVJUrv97W/w7W/nsDnHHHDOOfmYYVOSNCNHOCVJUtkmTYITT4QLL8yf114brrsOVl652LokSbXJwClJksry9NOwzz7w739Dr17ws5/BT38KvXsXXZkkqVYZOCVJUpumToVTToHTT4empjyaed11eXRTkqS2GDglSVKrXngB9t0XXnwRIuCYY+DUU6Fv36IrkyTVAwOnJEn6milT4LTT4Iwz4Msv82RA116bZ6KVJKlcNTNLbUSsEBFHRsT1ETEqIpoiIkXEzmV8d8+I+HtEjI+IiRExPCIOi4g2f76Ofk+SpEb2zDP5cdnTToOvvoIjjsgjnIZNSVJ71dII5yHAke39UkRcAhwKTAYeBKYBQ4GLgaERsUtK6atKfU+SpEY1aRKcfDKce25+V/Nb34KrroIhQ4quTJJUr2ppJO9fwNnAbsBywCOz+0JE7EQOjWOBb6eUvpdS2gFYHngZ2AH4caW+J0lSo3r8cVhjDTj77Pz52GPz+5uGTUlSZ9RM4EwpXZlSOi6ldHNK6fUyv3ZiaX98SunVGe71AXnEFOCEFh6R7ej3JElqKF98AUcemR+XfeWVPAPtk0/CWWfBnHMWXZ0kqd7VbaCKiIHA2sBU4JZZz6eUHgHeAxYBBnf2e5IkNZqHHoLVVoNf/xp69Mjraj73HHznO0VXJklqFHUbOIE1S/sRKaVJrVzz7CzXduZ7kiQ1hAkT4OCDYehQePNNWH11ePbZvNzJHHMUXZ0kqZHUc+BcurR/u41r3pnl2s58T5Kkunf33bDqqvCb30Dv3nDKKTlsrulfsUqSqqCWZqltr36l/RdtXDOxtO9fge9JklS3xo7N72refHP+PGgQXHNNDp+SJFVLPY9wRmmfuuh7028QcVBpzc7h48aN6+htJEmquqYmuPJKWGmlHDbnmisve/Lkk4ZNSVL11XPg/Ly079fGNc3nPp/hWEe/939SSleklAallAYtuOCCsy1UkqQijBoFm24KBx4In30GW28NI0bAUUdBr3p+xkmSVDfqOXC+Vdov2cY1S8xybWe+J0lSXZg6NU8AtPrq8OijsNBCcOONcM89sNRSRVcnSepO6vnvN58v7VeJiDlbmXF2nVmu7cz3JEmqeY8/DgcdBCNH5s/77w9nnw3zz19sXZKk7qluRzhTSqOB54A+wC6zno+IjYGBwFjgyc5+T5KkWjZ+PBx6KAwZksPm8svndTavusqwKUkqTt0GzpLTS/szI2K55oMRsRBwaenjGSmlpgp9T5KkmnP77XlSoMsuy+9mnnQSvPhifn9TkqQiRUodnqy1oiJiLaaHPYCVycuSvAp80nwwpTR4lu9dChwCTAYeAKYBQ4F5gDuBnVNKX7Xw+3Xoe7MaNGhQGj58eNk/pyRJlfLWW3DEEXDXXfnz4MHw2986+6wkqWtFxD9SSoNaOldL73DOA6zbwvHl2/pSSunQiHgMOAzYGOgJjAKuBi5rbZSyo9+TJKloU6bkpU1OOw0mTYL+/eH00+Hgg6Fnz6KrkyRpupoZ4axXjnBKkrrSgw/CYYfBv/+dP++xRw6fiy5abF2SpO6rXkY4JUlSK95/H44+Oi9vArDCCnDJJTB0aLF1SZLUlnqfNEiSpIb25Zdw0UWw4oo5bPbtC7/8ZZ4UyLApSap1jnBKklSjnn4aDjkEni+tCv2978Gvfw1LL11sXZIklcsRTkmSaswnn8CPfgTrrZfD5je/Cf/7v3k2WsOmJKmeGDglSaoRTU1w9dX5/cwrrshrap54IowcCdtvX3R1kiS1n4/USpJUA556Cg4/HJonPt900zwp0EorFVuXJEmd4QinJEkFGjsW9t03Pz47fDgsvjj84Q95+RPDpiSp3jnCKUlSAaZOzRMAnXIKfP459OkDxxyTH6Ht16/o6iRJqgwDpyRJXey+++DII+GVV/Ln7beH886DZZctti5JkirNwClJUhd57TU46qg82yzkyYEuuAC23rrYuiRJqhbf4ZQkqcomToSf/hRWWSWHzf794Zxz4J//NGxKkhqbI5ySJFVJUxPceCMcdxyMGZOP7bsvnH46LLJIoaVJktQlDJySJFXB44/nx2efeSZ/XmcduOgiWHfdYuuSJKkr+UitJEkV9MYbsOuuMGRIDpuLLAJXX53X2TRsSpK6G0c4JUmqgPHj4Ze/hAsvzEuezDknHHts3lzmRJLUXRk4JUnqhC+/hCuugJNPho8+ysd++MMcPgcOLLY2SZKKZuCUJKkDUoJ774VjjoGXX87HNtoIzj0XBg0qtjZJkmqFgVOSpHZ66SU4+mi4//78edll4eyzYdgwiCi2NkmSaomTBkmSVKYxY+Cgg2CNNXLYnG8+OO88GDkSdtjBsClJ0qwc4ZQkaTbGj4ezzoLzz4dJk6BXL/jxj+HnP4cBA4quTpKk2mXglCSpFVOmwKWX5gmAPv44H9tpp/x5hRWKrU2SpHpg4JQkaRZffQU33AD/8z/w9tv52EYbwZlnwuDBxdYmSVI9MXBKklSSEtx3H5xwAvzzn/nYqqvCGWfAttv6jqYkSe1l4JQkCXj2WTjuOHj44fx5iSXg1FNhr72gZ89CS5MkqW4ZOCVJ3dqrr8JJJ8Ett+TP3/hG/nzYYdC3b7G1SZJU7wyckqRu6Z134LTT4Oqr8zubffvCT34Cxx+flzuRJEmdZ+CUJHUr778Pp58Ov/kNTJ0KPXrA/vvD//t/MHBg0dVJktRYDJySpG7ho4/yWpoXX5zX0oyAPfaAk092iRNJkqrFwClJamiffQbnnQfnnw8TJ+Zjw4bBKafAaqsVW5skSY3OwClJakgTJ8JFF8HZZ8Onn+ZjW2+dZ54dNKjY2iRJ6i4MnJKkhjJpElx+eX5Pc9y4fGzjjfMEQUOGFFubJEndjYFTktQQJk+GK6/MQXPMmHxs3XXhl7+EzTbL72xKkqSuZeCUJNW1SZPgiivgzDPzDLQAa6yRRzS33dagKUlSkQyckqS69J//5EdnzzoLPvggH1tjDfif/8mTAvXoUWx9kiTJwClJqjMTJ8Jll8E558CHH+Zja68NP/85bLedI5qSJNUSA6ckqS58/jlccgmce25eUxNgnXXyOpo+OitJUm0ycEqSatr48XDxxXktzU8+yccGD85Bc6utDJqSJNUyA6ckqSZ9/HFeR/PCC+Gzz/KxIUNy0Bw61KApSVI9MHBKkmrK6NF5NPOKK/LEQJDX0Tz5ZNhkE4OmJEn1xMApSaoJo0blGWevvx6mTcvHttkGTjwRNtyw2NokSVLHGDglSYV69lk44wy44w5IKS9nsvvucMIJsPrqRVcnSZI6w8ApSepyKcFDD8Hpp8ODD+ZjffrAfvvBscfCsssWW58kSaoMA6ckqcs0NcGdd+YRzWefzcf694dDDoGf/AQWXbTY+iRJUmUZOCVJVfef/8B118H558Mrr+RjCy6YQ+ahh8J88xVbnyRJqg4DpySpaj74AC65BC69NC9zArDkknDMMbD//jDXXMXWJ0mSqsvAKUmquJEj89Im118PU6bkY+usA0cfDTvtBL38v48kSd2C/8uXJFVE80RA554L996bj0XA97+fg+aQIa6hKUlSd2PglCR1ytSpcNNNeUTzxRfzsTnnhH33ze9ofutbhZYnSZIKZOCUJHXIxx/Db38LF10EY8bkYwsvDD/+MRx8MCywQLH1SZKk4hk4JUnt8sILOWTecANMnpyPrbIKHHUU7Lkn9O1bbH2SJKl2GDglSbM1bVpeP/Oii+Dvf59+fJtt4IgjYKutfD9TkiR9nYFTktSqcePyY7OXXgrvvZeP9e8P++0Hhx3m+5mSJKltBk5J0tf84x95NPOmm6Yva7Liivn9zB/+MIdOSZKk2TFwSpKAPNvs7bfnoPnEE/lYBGy3HRx+OGy+uY/NSpKk9jFwSlI398YbcMUVcM018OGH+di888J//Vd+bHaZZYqtT5Ik1S8DpyR1Q19+CX/+M1x+OfzlL9OPr7ZaDpl77QVzz11cfZIkqTEYOCWpGxk9Gq68Mm/Na2f27Qu77QY/+hEMHuxjs5IkqXIMnJLU4L76Ko9iXn453H03NDXl4yusAAcfnCcBmn/+YmuUJEmNycApSQ1qzJj8XuZvfwtvv52P9e4Nu+6ag+ZGGzmaKUmSqsvAKUkNZOrU/G7m1VfDvfdOH81cZpn8yOy++8JCCxVaoiRJ6kYMnJLUAEaMyCHz97+HcePysd69YdiwHDQ33xx69Ci2RkmS1P0YOCWpTo0fD3/8I1x1FTzzzPTjq66alzT5wQ9gwQWLq0+SJMnAKUl1JCV45JE8mnnrrTBpUj4+zzyw556w//4waJDvZkqSpNpg4JSkOvDGG/CHP8Dvfgevvz79+Kab5tHMHXaAueYqrj5JkqSWGDglqUZ98gncckt+L/Pxx6cfHzgQ9tsvTwC0zDKFlSdJkjRbBk5JqiFTpsA99+SQeffdedZZyKOXO+4Ie+8NQ4dCz57F1ilJklQOA6ckFSwleOIJuP76PAnQp5/m4z16wBZb5JC5ww7Qr1+xdUqSJLWXgVOSCvLKK/m9zOuvz+9oNlt99Rwy99gDFlusuPokSZI6y8ApSV3orbfyKOYf/wjPPz/9+GKL5WVM9t4bVlutsPIkSZIqysApSVX23ntw8805ZD799PTj88yTH5Xde2/YZBPfy5QkSY3HwClJVfDhh3mdzJtugscey+9pQp78Z/vtYffdYautoG/fYuuUJEmqJgOnJFXIJ5/A7bfnkcyHHoKmpnx8jjngu9+F3XbL+7nnLrZOSZKkrmLglKROGDsW7rwzB82//Q2+/DIf790btt02h8ztt8+Pz0qSJHU3Bk5Jaqc334Q77sgh84knpj8u27NnXsZk993zu5nf+EaxdUqSJBXNwClJZXj55Rwwb7tt5tll55gDttwSdtoJvvc9GDCguBolSZJqjYFTklqQEjz3XA6Zt98Oo0ZNP9evX34Xc8cdYZttoH//4uqUJEmqZQZOSSqZNAkefBDuugv+/GcYM2b6ufnnz+9i7rQTbL65s8tKkiSVw8ApqVt7//0cLu+6Cx54IIfOZostBsOG5ZHMjTbKEwFJkiSpfAZOSd1KSvDCCzlg3nUXDB8+8/m114bttsvbmmtCRDF1SpIkNQIDp6SGN3FiXrLknnvyaOa7704/17dvfkR2u+3ypD+LLVZcnZIkSY3GwCmp4aQEL70E992Xt8ceg2nTpp9fdNEcLrfbDoYOhbnmKq5WSZKkRmbglNQQPvkE7r8f/vKXHDLff3/6uR49YPBg2GqrHDTXWisfkyRJUnUZOCXVpa++yu9fNo9iPvMMNDVNP7/oorD11nnbfPM8y6wkSZK6loFTUl1ICV5+OS9b8uCD8PDDMH789PO9e8Mmm+SAudVWsNpqTvgjSZJUNAOnpJr19tvTA+ZDD8HYsTOfX3bZHC633ho23RT69SumTkmSJLXMwCmpZnz4YZ5Ntjlgvv76zOcXWQQ22yxP9DN0KCy5ZDF1SpIkqTwGTkmFGTMG/v53ePTRvP3rXzOfn3fe/Jhsc8BcaSUfk5UkSaonBk5JXSIleOONHCybQ+asI5h9+8KQITlcbrZZnk22l/+VkiRJqlv+UU5SVTQ1wciRMwfMMWNmvqZfP9hgA9hoo7wNGpRDpyRJkhqDgVNSRUyYAE8/DU89BU8+mfeffjrzNQMGwIYbTg+Yq6/uCKYkSVIj8496ktqtqQn+/e/p4fLJJ2HEiPzY7IwWX3x6uNxoI1hxRejRo5iaJUmS1PUMnJJm65NPYPjw6QHz6ae/PnrZuzesuSast970bYklnORHkiSpOzNwSprJ+PHw3HM5YDZvb7zx9esWW2zmcLnWWr5/KUmSpJkZOKVubOJEeP75mcPlK698/bq+fWGNNWDw4BwuBw929FKSJEmzZ+CUuomxY+HFF+GFF6bvR436+nuXffrkyXwGDZq+rbyyk/tIkiSp/fwjpNRgpk3LQfLFF2fePvzw69f26gXf/vbM4XKVVXLolCRJkjrLwCnVqZRg9Oi81uWIEfCvf+VgOWIETJ369evnmSePXM64rbaa711KkiSpegycUo1rasrBcsSI6eFy5Mi8TZzY8neWWebr4XKppXznUpIkSV3LwCnViEmT4PXX86Q9r74KL788PVh+8UXL31lwwfx+5Sqr5P3qq+dHZOeZp2trlyRJklpi4JS60JdfwttvTw+Vr7wyfXvnna9P4NNsoYWmh8rm/cor58ApSZIk1apuHzgjYk/gEODbQE9gFHANcFlKqanI2lSfPv8c3nxz5u2NN+C11/II5rRpLX+vZ8/8KOy3vgXLLw8rrJDD5UorwQILdO3PIEmSJFVCtw6cEXEJcCgwGXgQmAYMBS4GhkbELimlrwosUTVo0iR4992WQ+Wbb8LHH7f9/YEDc6hs3pZfPu+XXhp69+6an0GSJEnqCt02cEbETuSwORbYKKX0aun4wsDfgB2AHwMXFlakutyUKfDee3mSnubt3Xdn/jy7QNm3b56gZ+ml87bMMnm/7LKw3HIw99xd8qNIkiRJheu2gRM4sbQ/vjlsAqSUPoiIQ4CHgRMi4iIfra1/kybB2LHw/vt539KvR49uea3KWfXqBYsvPj1UNgfK5m2RRaBHj6r/SJIkSVLN65aBMyIGAmsDU4FbZj2fUnokIt4DFgcGA090bYWancmT80jjRx+1vI0bN3OYnDChvPv27JnD5MCBsMQSeWv+dfN+4YUNlJIkSVI5umXgBNYs7UeklCa1cs2z5MC5JgbOikspL/Uxfjx89tns959+OnOgbG2ZkNb06ZNHHhdZBBZdtOVfDxyY9z17VudnliRJkrqb7ho4ly7t327jmndmubYhpQRNTXn76qs8g+qUKdO3qVNn/tza9p//wMSJOQh+8cXsfz1hQv79Oqp37zxz6wILwIAB03894zZjmPzGNyCicu0mSZIkafa6a+DsV9q3NU42sbTvX+VaKm7YMHjiiZmDZPOvZ/3c2rqPXWHOOWHeeWG++Wa/n2++vOZkc5js398AKUmSJNW67ho4m6NKh+JWRBwEHATwzW9+s1I1Vcynn+Z3GMvVo8f0rXfv/PjpHHOUv/Xpk2debd769Zv9r/v3z9+VJEmS1Li6a+D8vLTv18Y1zec+n/VESukK4AqAQYMGFThG2LLbb4cvv8zvIs4YJmf83PzrCEcKJUmSJFVHdw2cb5X2S7ZxzRKzXFs3BgwougJJkiRJgu66uMPzpf0qETFnK9esM8u1kiRJkqR26JaBM6U0GngO6APsMuv5iNgYGAiMBZ7s2uokSZIkqTF0y8BZcnppf2ZELNd8MCIWAi4tfTwjpdTU5ZVJkiRJUgPoru9wklK6NSIuAw4BXoqIB4BpwFBgHuBO4OICS5QkSZKkutZtAydASunQiHgMOAzYGOgJjAKuBi5zdFOSJEmSOq5bB06AlNINwA1F1yFJkiRJjaY7v8MpSZIkSaoiA6ckSZIkqSoMnJIkSZKkqjBwSpIkSZKqwsApSZIkSaoKA6ckSZIkqSoMnJIkSZKkqjBwSpIkSZKqwsApSZIkSaoKA6ckSZIkqSoMnJIkSZKkqjBwSpIkSZKqwsApSZIkSaoKA6ckSZIkqSoMnJIkSZKkqjBwSpIkSZKqIlJKRddQ1yJiHPB20XW0YAHgo6KL6KZs++LY9sWy/Ytj2xfHti+ObV8c2744tdr2S6aUFmzphIGzQUXE8JTSoKLr6I5s++LY9sWy/Ytj2xfHti+ObV8c27449dj2PlIrSZIkSaoKA6ckSZIkqSoMnI3riqIL6MZs++LY9sWy/Ytj2xfHti+ObV8c2744ddf2vsMpSZIkSaoKRzglSZIkSVVh4KwDEbFnRPw9IsZHxMSIGB4Rh0VEh/75Vfp+jaxSbRUR10ZEamMbVa2fod5ExAoRcWREXB8RoyKiqdRGO3fyvvb7MlS6/e375YmI3hExNCLOjYinIuL9iJgaEe9FxK0RsUkn7m3fb0M12t5+X76IODwibo6IlyPi44iYFhHjIuKBiNgrIqKD97Xfz0al295+3zkR8asZ2uqYDt6jJvt9ryJ/c81eRFwCHApMBh4EpgFDgYuBoRGxS0rpq6Lu18iq1FaPA6+1cPz9ztTaYA4BjqzkDe337VLx9i+x77dtY+D+0q/HAv8AvgBWBnYCdoqIU1NKP2/PTe37ZalK25fY72fveGAh4F/AE+S2XxLYjNxXd46IHVNKTeXe0H5ftoq3fYn9vp0iYh3gOCABHf1Lltrt9ykltxrdyP+jS+R/QZef4fjCwMjSuSOLul8jb1Vo+2tL39m36J+t1jfgAOAsYFdgWeDhUtvtXAv/LBt9q0L72/fLa6fNgFuBDVs4txvwZakdN23HPe37xbW9/b78thoCzN3C8VXIfwGQgP3acT/7fXFtb7/v2D+HOYARwHvAHaU2PKad96jpfu9jBbXtxNL++JTSq80HU0ofkEchAE5oxzB5pe/XyGyrgqSUrkwpHZdSujml9HoFbuk/y3aoQvurDCmlh1JKO6eU/t7CuT+S/yAHsFc7bmvfL0OV2l5lSik9llL6ooXjI4BLSh+3aMct7fdlqkLbq2NOIT9RcTAwvoP3qOl+3+3/ZatVETEQWBuYCtwy6/mU0iPkvwlZBBjc1fdrZLZV4/CfpRrI86X9wHIutu9XVLvaXhX1ZWk/uZyL7fcV1a62V8dExLrA0cANKaW7OniPmu/3vsNZu9Ys7UeklCa1cs2zwOKla5/o4vs1smq21aYR8W2gH/AB8Bhwf2r/+xEqj/2+dtj3O2f50r7cd6Ds+5XT3rafkf2+gyJiafKID0Axgob8AAAKfElEQVS5fxC331dAB9t+Rvb7MkREX+B3wCd0bu6Emu/3Bs7atXRp/3Yb17wzy7Vdeb9GVs22+mELx0ZGxO4ppZfaeS/Nnv2+dtj3OygiFgH2LX28rcyv2fcroINtPyP7fZkiYj/yBE69yaPJ65OfxDs9pXRHmbex33dAhdp+Rvb78vwSWAHYPaX0USfuU/P93kdqa1e/0v5rz9bPYGJp37+A+zWyarTVC8AR5Bfx+wGLAd8DXiQ/t/9ARCze/lI1G/b74tn3OyEiegHXA/MCD7bjkSv7fid1ou3Bft8RGwD7AHsCG5WO/Q/5/bZy2e87phJtD/b7skXE+sBPgDtL74p3Rs33ewNn7WqeEjnV6P0aWcXbKqV0QUrpopTSyJTSFyml91NKdwPfAZ4iT0t+Ytt3UQfY7wtm3++0y8nT2o+mfZPW2Pc7r6Ntb7/vgJTSASmlAOYiB5YLgF8AT0XEYmXexn7fARVqe/t9mSJiTuAaYAJ5GZNO37K0r9l+b+CsXZ+X9v3auKb53OdtXFOt+zWyLmurlNJU4PTSx207cy+1yH5fo+z7sxcRFwL/RV6eYGhKaWw7vm7f74ROtn2r7Pezl1KaVAosx5LDyerkdQTLYb/vhE62fVv3td/P7FfAt4CjUkqVWJu05vu9gbN2vVXaL9nGNUvMcm1X3q+RvVXad1VbjSrtfcyk8t4q7e33tcm+34qIOJf8aNo4cuB5dTZfmdVbpb19v50q0PazY78v3zWl/XYR0buM698q7e33ndfetp8d+/10OwBNwD4R8fCMG7B16ZpDSseuLON+b5X2NdvvnTSodjVPw75KRMzZyqxT68xybVfer5F1dVsNKO0ntnmVOsJ+X9vs+y2IiLOAo4CPgS1SSiM7cBv7fgdUqO1nx35fvs/Iy3P0AuYnz3jaFvt95bS37WfHfj+zHuSJmlqzTGmbr4x71Xy/d4SzRqWURgPPAX2AXWY9HxEbk2cSGws82dX3a2QFtNWupf2zFbiXZmC/r3n2/VlExBnAscCn5MDzYkfuY99vv0q1fRns9+XbiBx4PgNmO4un/b6i2tX2ZbDfl6SUlkopRUsbeZkUgGNLx9Yo43413+8NnLWt+Xn3MyNiueaDEbEQcGnp4xkzrmsUEadHxKiIOJ2va/f9urGKtX1ErBER34uInrMc7xURR5Ef3QI4v+I/RTdhvy+Wfb8yIuJU4HjyH/C2SCnN9m+i7fuVUcm2t9+XLyI2jIgfRMQcLZzbALiq9PGqlNJXM5yz33dSpdvefl999dzvfaS2hqWUbo2Iy4BDgJci4gFgGnnWvHmAO/n6y9yLktf0WbRC9+uWKtz2SwF3AJ9ExCvAu+RpqVcjTxneBByfUvpLdX6a+hIRazH9P46Qp1IH+FVEHNN8MKU0eIZr7PcVUuH2Xwr7flkiYnvgZ6WPrwGHR0RLl45KKZ0xw2f7fidVoe2Xwn5frmXJ7wpeHBHPkUdg+peON/+3527yEh0zst93XqXbfins99VWt/3ewFnjUkqHRsRjwGHkZ717kl+8vhq4rL1/U1Hp+zWyCrbVi8CF5GnBlwTWJE9d/S75P/aXpJT+UeHy69k8wLotHF++oze037dLJdvfvl+++Wf49aDS1pJHgDNaOfc19v2yVLrt7fflewQ4FdiQPGvn+uQlHsYCtwHXp5TubO9N7fdlqXTb2+8LVsv9PlKq2SVbJEmSJEl1zHc4JUmSJElVYeCUJEmSJFWFgVOSJEmSVBUGTkmSJElSVRg4JUmSJElVYeCUJEmSJFWFgVOSJEmSVBUGTkmSZhARqQPbtaXvblL6/HCxP0XnRcTxpZ9l607cY62IaIqIcypZmySpfvQqugBJkmrM71o4tgiwFfAFcGsL5x+rakVdLCIWBU4CHk0p3dfR+6SUnouI24EjIuI3KaVXK1akJKkuREqp6BokSappEbEJ8Dfg7ZTSUm1cNxfwTeA/KaV3uqa6youIK4ADgaEppYc6ea/VgH8Ct6WUdq5EfZKk+mHglCRpNsoNnI0gIgYA7wJjgOVSBf6gEBHPAmsCy9RzEJcktZ/vcEqSVCGtvcMZEUuVjr8VET0i4qiIGBERkyLi3Yg4rzQ6SkR8IyIuKF07JSJejYij2vg9IyJ2j4i/RsRHpe+8ExG/jYilOvBj7A/0Ba5rKWxGxHwR8atS/f+Z4Wd4OCJObOWevwN6Aj/qQD2SpDpm4JQkqWvdAJwCvAn8FZgb+G/gtoiYH3ga2A14lvxu6FLAuRHx01lvFBG9ye+U3ggMAUYCfyK/a3oA8FxEDGpnfcNK+wda+P3mAh4HTgQWKF1zJ/AasDJwciv3bL7X99tZiySpzjlpkCRJXWdJYDLwrZTSGICIWAJ4HtgaeAR4Edg7pTS5dP67wJ+BEyLigpTSf2a436nAjsCjwA9SSu82n4iIHwMXATdFxIoppS9nV1wpUK4DTAP+0cIlO5OD5d3AsBnvGRE9gY1bufW/gU+BVSJi4ZTSB7OrRZLUGBzhlCSpax3RHDYBUkqjgetLH5cEDmkOm6Xzd5Mn3ekP/N9oZWk09AhgIrDLjGGz9L2LycFwWWCbMmtbBegNvDljDTNYuLR/YNYAm1L6qrUJhkqP5r5c+rhGmbVIkhqAgVOSpK4zDWgplL1W2g9PKX3Uwvnm5UQWm+HYpsCcwCMppQ9b+f0eKe3XK7O+hUr7j1s5/0xpf3xE7BUR85V5X4BPSvuF27xKktRQfKRWkqSuM7aVR1snlvbvtnBuxvN9Zzi2TGn/3YiY3UyyC5ZZ37yl/YSWTqaUHomIs4BjgN8DKSJGkd81vS2l9Jc27t18z/aEVElSnTNwSpLUdZo6eX5GPUv7fwNPzebap8u852el/TytXZBSOj4iLidPADQE2IC8ZueBEfFX4LuthOrme35aZi2SpAZg4JQkqT6NLu1fSintW6F7Nj+aO6Cti1JKbwIXlDYiYgh5ptwtycuqXNHC15rv2drjv5KkBuQ7nJIk1acHyO+Ebt7OdynbMgKYAiwdEXOW+6WU0mPAtaWPq896PiICWLH08flO1ihJqiMGTkmS6lBpaZFLyO9E/ikiVpz1moj4RkQcEBFlTdSTUppEfvy2N7B2C/fbISI2iogesxyfE9i89PHtFm69IvANYEQbExxJkhqQj9RKklS/jiPPXLsr8K+IeAF4kzy50BLASkCf0r7ctS/vBDYiB8jHZjm3MXAkMC4ingfGkScaWh+YHxgF/KaFezaH0f8tswZJUoNwhFOSpDqVUpqWUtqNPIHPn8nh8/vkANgLuAHYAXi9Hbe9FpgE/LD0KOys584EXgFWBXYBvkNe1uW/ge+klMa3cM99gK9oOYxKkhpY5LWYJUmSstIstD8ChqaUWlo3tD33Wg34J3nZlJ0rUZ8kqX4YOCVJ0kwiYhHyKObzKaWNO3mvW4HtgVVSSq9Woj5JUv3wkVpJkjSTlNJY4DRgo4jYuqP3iYi1gB2BiwybktQ9OcIpSZIkSaoKRzglSZIkSVVh4JQkSZIkVYWBU5IkSZJUFQZOSZIkSVJVGDglSZIkSVVh4JQkSZIkVYWBU5IkSZJUFf8fvikoQI+0bcIAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 1080x720 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"def rocket2(state,dmdt=0.05, u=250,c=0.18e-3):\n",
" g = 9.81\n",
" dstate = np.zeros(np.shape(state))\n",
" dstate[0] = state[1]\n",
" dstate[1] = (u/state[2])*dmdt - g - (c/state[2])*state[1]**2\n",
" dstate[2] = -dmdt\n",
" return dstate\n",
"t = np.linspace(0,(m_0-m_f)/dmdt,N)\n",
"dt = (m_0 - m_f)/(dmdt*N)\n",
"num_sol = np.zeros([N,3])\n",
"num_sol[0,0] = y_0\n",
"num_sol[0,1] = v_0\n",
"num_sol[0,2] = m_0\n",
"for i in range(N-1):\n",
" num_sol[i+1] = heun_step(num_sol[i], rocket2, dt)\n",
"\n",
"\n",
"fig = plt.figure(figsize=(15,10))\n",
"plt.plot(t,num_sol[:,0],linewidth=2, color='b')\n",
"plt.xlabel('Time (s)')\n",
"plt.ylabel('Height (m)')\n",
"plt.title('Altitude vs Time of Rocket');\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## References\n",
"\n",
"1. Math 24 _Rocket Motion_. <https://www.math24.net/rocket-motion/\\>\n",
"\n",
"2. Kasdin and Paley. _Engineering Dynamics_. [ch 6-Linear Momentum of a Multiparticle System pp234-235](https://www.jstor.org/stable/j.ctvcm4ggj.9) Princeton University Press \n",
"\n",
"3. <https://en.wikipedia.org/wiki/Specific_impulse>\n",
"\n",
"4. <https://www.apogeerockets.com/Rocket_Motors/Estes_Motors/13mm_Motors/Estes_13mm_1_4A3-3T>"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}