Permalink
Cannot retrieve contributors at this time
Name already in use
A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
compmech-project03/CompMech03-IVPs_project (1).ipynb
Go to fileThis commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
529 lines (529 sloc)
103 KB
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "# Initial Value Problems - Project\n", | |
| "\n", | |
| "" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "We are going to end this module with a __bang__ by looking at the flight path of a firework. Shown above is the initial condition of a firework, the _Freedom Flyer_ in (a), its final height where it detonates in (b), the applied forces in the __Free Body Diagram (FBD)__ in (c), and the __momentum__ of the firework $m\\mathbf{v}$ and the propellent $dm \\mathbf{u}$ in (d). \n", | |
| "\n", | |
| "The resulting equation of motion is that the acceleration is proportional to the speed of the propellent and the mass rate change $\\frac{dm}{dt}$ as such\n", | |
| "\n", | |
| "$$\\begin{equation}\n", | |
| "m\\frac{dv}{dt} = u\\frac{dm}{dt} -mg - cv^2.~~~~~~~~(1)\n", | |
| "\\end{equation}$$\n", | |
| "\n", | |
| "If we assume that the acceleration and the propellent momentum are much greater than the forces of gravity and drag, then the equation is simplified to the conservation of momentum. A further simplification is that the speed of the propellant is constant, $u=constant$, then the equation can be integrated to obtain an analytical rocket equation solution of [Tsiolkovsky](https://www.math24.net/rocket-motion/) [1,2], \n", | |
| "\n", | |
| "$$\\begin{equation}\n", | |
| "m\\frac{dv}{dt} = u\\frac{dm}{dt}~~~~~(2.a)\n", | |
| "\\end{equation}$$\n", | |
| "\n", | |
| "$$\\begin{equation}\n", | |
| "\\frac{m_{f}}{m_{0}}=e^{-\\Delta v / u},~~~~~(2.b) \n", | |
| "\\end{equation}$$\n", | |
| "\n", | |
| "where $m_f$ and $m_0$ are the mass at beginning and end of flight, $u$ is the speed of the propellent, and $\\Delta v=v_{final}-v_{initial}$ is the change in speed of the rocket from beginning to end of flight. Equation 2.b only relates the final velocity to the change in mass and propellent speed. When you integrate Eqn 2.a, you will have to compare the velocity as a function of mass loss. \n", | |
| "\n", | |
| "Your first objective is to integrate a numerical model that converges to equation (2.b), the Tsiolkovsky equation. Next, you will add drag and gravity and compare the results _between equations (1) and (2)_. Finally, you will vary the mass change rate to achieve the desired detonation height. \n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "1. Create a `simplerocket` function that returns the velocity, $v$, the acceleration, $a$, and the mass rate change $\\frac{dm}{dt}$, as a function of the $state = [position,~velocity,~mass] = [y,~v,~m]$ using eqn (2.a). Where the mass rate change $\\frac{dm}{dt}$ and the propellent speed $u$ are constants. The average velocity of gun powder propellent used in firework rockets is $u=250$ m/s [3,4]. \n", | |
| "\n", | |
| "$\\frac{d~state}{dt} = f(state)$\n", | |
| "\n", | |
| "$\\left[\\begin{array}{c} v\\\\a\\\\ \\frac{dm}{dt} \\end{array}\\right] = \\left[\\begin{array}{c} v\\\\ \\frac{u}{m}\\frac{dm}{dt} \\\\ \\frac{dm}{dt} \\end{array}\\right]$\n", | |
| "\n", | |
| "Use [two integration methods](../notebooks/03_Get_Oscillations.ipynb) to integrate the `simplerocket` function, one explicit method and one implicit method. Demonstrate that the solutions converge to equation (2.b) the Tsiolkovsky equation. Use an initial state of y=0 m, v=0 m/s, and m=0.25 kg. \n", | |
| "\n", | |
| "Integrate the function until mass, $m_{f}=0.05~kg$, using a mass rate change of $\\frac{dm}{dt}=0.05$ kg/s. \n", | |
| "\n", | |
| "_Hint: your integrated solution will have a current mass that you can use to create $\\frac{m_{f}}{m_{0}}$ by dividing state[2]/(initial mass), then your plot of velocity(t) vs mass(t)/mass(0) should match Tsiolkovsky's_" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 2, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "import numpy as np\n", | |
| "import matplotlib.pyplot as plt\n", | |
| "%matplotlib inline\n", | |
| "plt.rcParams.update({'font.size': 22})\n", | |
| "plt.rcParams['lines.linewidth'] = 3" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 3, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "def simplerocket(state,dmdt=0.05, u=250):\n", | |
| " '''Computes the right-hand side of the differential equation\n", | |
| " for the acceleration of a rocket, without drag or gravity, in SI units.\n", | |
| " \n", | |
| " Arguments\n", | |
| " ---------- \n", | |
| " state : array of three dependent variables [y v m]^T\n", | |
| " dmdt : mass rate change of rocket in kilograms/s default set to 0.05 kg/s\n", | |
| " u : speed of propellent expelled (default is 250 m/s)\n", | |
| " \n", | |
| " Returns\n", | |
| " -------\n", | |
| " derivs: array of three derivatives [v (u/m*dmdt-g-c/mv^2) -dmdt]^T\n", | |
| " '''\n", | |
| " dstate=np.zeros(np.shape(state))\n", | |
| " dstate=np.array([state[1],((u/(state[2])*dmdt)),-dmdt])\n", | |
| " return dstate" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 4, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "def rk2_step(state, rhs, dt): ##explicit\n", | |
| " mid_state = state + rhs(state) * dt*0.5 \n", | |
| " next_state = state + rhs(mid_state)*dt\n", | |
| " return next_state" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 5, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "def heun_step(state,rhs,dt,etol=0.000001,maxiters = 100):##implicit\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": 6, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "<matplotlib.legend.Legend at 0x7f7ee47fef28>" | |
| ] | |
| }, | |
| "execution_count": 6, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEaCAYAAABEsMO+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdeXxU1fn48c/JTPY9ZGMHRQFtKRSqUKiAUMSCIKCg1BYERNnUWvwJtVpc2mLdK5sisiiKUZBNBVwQEdeIFvUbxIWd7Ptkz8z5/XFnkkkyCZlkskzyvF+vvMa599w7Jwh5cs59znOU1hohhBCitfFp6Q4IIYQQrkiAEkII0SpJgBJCCNEqSYASQgjRKkmAEkII0SqZW7oD3i46Olr36NGjpbshhBBe5csvv8zQWsfU1UYCVCP16NGDxMTElu6GEEJ4FaXUyfO1kSk+IYQQrZIEKCGEEK1Sqw1QSql/KaW0/WtxHe2mK6UOKqVylVIWpVSiUmqBUqrO762h1wkhhGgerfKHsVLqN8D/A+qsw6SUWglsBgYBB4F3gIuBFcDrSimTJ68TQgjRfFpdgFJK+QMbgFRgRx3tpgDzgRSgn9Z6vNZ6EnARkARMAhZ66johhBDNq9UFKOBB4BLgNiC3jnZL7a/3aK1/cBzUWqcC8+xvl7iYsmvodZ5zJAGe/AUsizBejyQ02UcJIYS3alUBSil1OfBX4GWt9a462nUBBgKlwGvVz2utDwBngXhgcGOv86gjCbDrdvLzzvDvqAhOFyTDtltg44Qm+TghhPBWrSZAKaUCgI1AFnDHeZoPsL9+p7UuqqXNF9XaNuY6z3nvQT42aSZ26cjL4aH8s0Ok8aDt+AEJUkII4aTVBCjgn0BvYJHWOuM8bXvaX+ta6HWqWtvGXOc5uWcIt9nINBl5GIeCAtkTHGScO35ApvuEEMKuVVSSUEr9FrgT2K61frUel4TYXwvqaGOxv4Z64LoqlFJzgbkA3bp1q+NWLoR34dLc09yYZ2FzuPERyztE8tuiYsJtNtg+D/pNde+eolXTWpOfn09eXh6FhYVYrdaW7pIQHmM2mwkPDycqKgqz2bMhpcUDlFIqEFgP5GFk19XrMvuru9sBN/S6KrTWzwHPAQwaNMi9e426H7bdwsLsHN4JDiTNbCbLZOLRqAgezsgCWznsvgvGP9GYLopWQmtNWloaBQUFREVFER8fj8lkQil1/ouFaOW01pSWlpKZmcnp06fp3r07Pj6em5hrDVN8/8JYg3SX1jq5ntfk219D6mjjOJfvdKyh13lOv6ng40eI1tybmV1xeEdoCB8HBBhvEtc1yUeL5pefn09BQQHdu3cnIiICs9kswUm0GUop/P396dixI2azmezs7PNf5IbWEKAmATZghlLqA+cvYKy9zTz7seft70/YX7vXcd+u1do25jrPunYlAFcWFnGVpXK28cHoKAodP7we69NkHy+aT15eHlFRUZhMsvZbtF1KKSIiIigoqOvpiftaQ4ACox/DXXzF2c9fYH8/yP7+K/vrpfYpQld+U61tY67zrH5ToedwAJZkZhNufyZx1tfMM5HhRhtLsmT1tQGFhYWEhNQ1YBeibQgKCqKoqLbk6IZp8QClte6htVauvjDSzgHuth/rb7/mNHAY8AOur35PpdRwoAtGtYhPnD6rQdc1iRk7AYi22bgnq3JYvDkslK/9/Yw3ktXn9axWq4yeRLvg4+ODzWbz7D09erfm9W/76yNKqV6Og0qpWGCV/e1yrXX1P7GGXud5g2YDMN5SyNBC4zcPrRTLoqModbTZPs/1tcJryDMn0R40xd9zrw1QWuvXgdUYVR++UUrtUkptA37AKJW0HaP4q0euaxLjn4CQjijg/swsAu2/ffzk58cax1SfI6tPCCHaGa8NUABa6/nAHzGm7YYDVwE/YhR7naK1drngpKHXNYnFRwFFp3Ird2blVBx+ITyMI46pPsnqE0K0Q606QGmtZ9qfPT1WR5uXtdZDtdZhWutgrfVArfXK803RNfS6JjH5OQBuyLcwqKgYAKtS3BvdgSLJ6hNCtFOtOkC1G/asPh/g4YxMguxTfSf8fHk6MsJoI1l9oh2w2Wx069YNpRSxsbGUlZW1dJcA2LBhA0opZs6c2Syft2zZMpRSLFu2rFk+r7qZM2eilGLDhg0t8vkOEqBaC3tWX+dyK/c4LeDdHB7KpwH+xhvJ6hNt3L59+zh9+jQA6enp7NpV66YGXuvEiRMopejRo0dLd6XVkwDVmtiz+iZZChheWLme4O8xHcjzsU/1vXFrS/RMiGbxwgsvANC5c+cq79ubhQsXkpSUxMKF7XvvVAlQrYlTVt+yjEwi7At4U81mHomKNNpom0z1iTYpKyuLnTt3opRiy5YtmEwm9uzZw7lz51q6a80uOjqaPn36EB0d3dJdaVESoFobe1ZftNXGfRlZFYd3hobwTpC9+MXxA5J6Ltqcl156iZKSEkaMGMGwYcMYM2YMVquVTZs2uWyvlKpYe/Pqq68yZMgQQkJCCA0NZdSoUXz00Ucur/vss8+4++67GTRoEHFxcfj5+dGpUyeuu+46Pv3003r3d9OmTSilGDt2bK1tvvnmG5RSdO7cmfLycmbOnEnPnsZOPidPnqz4HqpP+Z3vGVRSUhJz586lV69eBAYGEhkZSb9+/Vi8eDEnT1bdTWjr1q3MmjWLSy+9lIiICAICAujVqxcLFiyomE5trSRAtUb2rL4xhUWMc6rV94/oDiQ7qhIkrpPnUaJNWb9+PUBFIsLNN99c5Xht7r//fqZPn46fnx/jxo2jS5cuvP/++4waNYpPPqlZEObee+/lySefpKysjMsuu4wJEybQoUMHtm7dyrBhw3jttRqbbbt0ww03EBsby759+/jxxx9dtlm50qi7OXfuXMxmM8OGDWPKlCkABAcHM2PGjIqv6667rl6fu2nTJvr378/atWvRWjN+/HiGDx+OzWbj8ccfZ//+/VXaT5s2jYSEBIKDgxk9ejS///3vKSkpYdWqVfz617/m2LFj9frcltDi220IF/pNhQ8fh4yjLM3M4nCAP8lmM/kmH5bEdmBdcprxP072jvJ6PZa82dJdaLATy8d57F5fffUVX3/9NaGhoRU/qCdOnEiHDh04duwYH330EcOGDXN57cqVK/n8888ZOHAgYGQC3nbbbaxdu5b777+fd955p0r7xYsXs3nzZuLi4qoc37VrF1OmTOG2225j3LhxBAUF1dlnPz8/5s6dy8MPP8yaNWt47LGqq2Hy8vLYvHkzZrOZW265BYA5c+YwevRotm7dSnR0tNtZcl988QWzZ89Ga83zzz/PrFmzqlRwSEpKqnHNyy+/zPjx46t8P+Xl5TzwwAM8/PDD3HHHHbz99ttu9aO5yAiqtVr4GaAIt2keScvARxvbTh0OCGBtRJjRxlYuz6NEm7BunbEYferUqRU/SP38/Jg+fTpQd7LEAw88UBGcwKgJ9/DDDwNw8ODBGqnqY8eOrRGcAK655hquv/56srKyaoxCajNv3jzMZjPr16+nuLi4yrmNGzdisViYNGkSnTp1qtf9zuef//wn5eXlLF68mNmzZ9coL9S3b1/69u1b5Zjzn6mD2WzmoYceolOnTuzbt4/8/KbZXaixZATVmk1+DrbdwoCSUm7LyWWVfU3UmohwLi8q4dclJZWp5zKSEl6qpKSEV155Baic1nO4+eabeeaZZ3jttdf473//67Iy/Pjx42sci42NJTIykuzsbDIzM4mPj69yPiMjg927d/Ptt9+Sk5NDeXk5AN9++y0Ax44dY9y4848QO3XqxOTJk0lISGDLli1V1kmtXr0agAULFpz3PvVhtVp59913AWMk5o5jx46xZ88efvzxRywWS0VR1/Lycmw2Gz/++CMDBgzwSD89SQJUa9ZvKnz1Ehw/wNycPD4LCODLwABsSnFPbAdeP5tMuE0bqecSoLySJ6fJvNUbb7xBVlYWF110EUOHDq1ybsCAAfTv35+vv/6ahIQEZs2aVeP6bt26ubxvWFgY2dnZNUY2zz77LHfddReFhYW19ikvL6/e/b/99ttJSEhg1apVFQFq//79JCUlcemllzJ8+PB636suGRkZFBQUYDab6dWr1/kvwAhA8+fP5/nnn0fr2jf/duf7bU4yxdfazdgJ0X0wAcvTMwmzp56nmM08EN3B2Lte26QUkvBajum73Nxchg0bVuMrNTW1Srvq3NliPDExkXnz5lFWVsajjz7K0aNHK0YUWmuWLl0KUOcP8+qGDh3KgAED+OKLL0hMTAQqkyPmz59f7/s0haeffpq1a9fSsWNHtmzZwqlTpyguLkZrjdaaIUOGAO59v81JApQ3WPgZAPFWKw86pZ6/ExzEllD7lIeUQhJe6PTp07z33nsApKWlcejQoRpfycnJABw6dKjRGWevv/46Wmtuv/12Fi9eTO/evQkODq54llNbNt75LFq0CIBVq1Zx7tw5duzYQWhoKH/6058a1V9n0dHRBAUFUV5ezk8//VSvaxwZic8++yzTpk2ja9eu+Pv7V5xv6PfbXCRAeQt7lYlRhUVMy6t8oPmfDpF84ycbHArvtH79emw2G6NGjar4rd7V1/XXG/uLNrayRFaW8Qte165da5xLT0+vkfFXXzfeeCPR0dFs2bKF5cuXU15ezp///GdCQ0NrtPWz/3t1PPeqL5PJxOjRowF4/vnn63VNXd/vO++8Q3p6ult9aG4SoLzF+Ccg2pjGW5yVQ58SY0vDcqX4a1w0OY5pDimFJLyE1pqNG41Ns8830nCc37RpE1Zrw3fD6dOnT8V9LBZLxfH8/HxmzZpFTk5ObZfWKSAggDlz5lBUVMQzzzwD1D69FxMTg5+fH6mpqWRnZ7tsU5t7770Xk8nEY4895jJF/ejRoxw9erTiveP7Xb16dZXdbn/66Sduu+02tz67JUiA8iYLP4OQjgRozRNpGYRajb9wyWYzf4vpgA3keZTwGvv37+fnn38mKCiIyZMn19l27NixREdHk5yc3Kg1OzfffDNdu3bl8OHDXHDBBUyePJlJkybRo0cPEhMTXSZh1Nf8+fMx2RfSjxgxgksuucRlO19fX8aNG0d5eTkDBgzgj3/8I3PmzGHJkiXn/YzLLruM5557ruJ76dWrF9OmTePaa6/ll7/8JX379q1SDWPp0qX4+vry7LPP0rdvX2644QbGjBnDJZdcQteuXfntb3/b4O+3ObgVoJRSEUqpa5VSDyil1iiltiilVtvfT1RKRTRVR4XdYuO3o67l5TyckVlx+GBQIM+H29dHWZJhxeUt0Tsh6s1RIWLixIkup8Kc+fr6csMNNwCNm+aLjIwkMTGRuXPnEhISwptvvkliYiKTJ0/m8OHDLqfC6qtr164VI5bzpZavXbuW2bNnY7VaSUhIYN26dWzZsqVenzNr1iwOHz7MzJkzKSsrY/v27Xz44YeYTCbuvvturrzyyoq2Q4YM4fPPP2fcuHHk5uayY8cOzpw5w7333svevXvx9fVt8PfbHNT5sjeUUiZgMjAf+B3gWBnmvEJMO71+CKwC3mjWnWlbyKBBg7Qjc6fZ7L6rYpfdJyIjWG9fuOujNc+mpDG4uMTeudnG1KBoMUlJSTUWToq26X//+x/9+/enU6dOnDx5ErO5/a3icefvu1LqS631oLra1DmCUkrdCPwMbMHYGj0LeBN4FFiCEbSWAI8BbwHZwAjgVeAnpdQN9eqpcI/T86jbs3MYaN+F11gfFU2qc70+IUSzuP/++wFjXVR7DE5NodY/RaXUIWAwkA48DWzUWv/vfDdUSvUHZgI3ApuVUou01kPrvkq4beFn8FgfzJZkHk3P4PpOHck0m8gymbgrNpoXUlLx18C/u8HSUy3dWyHapJ07d7Jjxw6++eYbvvjiC3r06NHu93DypLpGUBcCdwHdtNZ31Sc4AWitv9Za3wl0Bf5qv49oCvatOWKsNv6TXlmv70iAPw93iDLmXUtyJWlCiCZy+PBhXnjhBY4ePcrYsWPZs2cPwcHBLd2tNqPOAKW1flprXdqQG2utS7XWTwEXNKxrol7sW3NcVlzC4qzKFNntoSG8HCaLeIVoSsuWLUNrTV5eHm+//Ta9e/du6S61KbUGKK11QW3n3KG1rr3glWi8flOhp1Hr66a8fCbkV67teDQqks8C7KvGZZNDIYSXkXVQbcGMnRVbxd+fmcUvSowsPqtSLI6N5oxZNjkUQnifegcopVSUUurXSqmoasc7KqU2KKW+Ukq9oZT6lee7Kc5r8VHwD8dfw1OpGUSXGxn+OSYTd8TGUOjYN0YqTQghvIQ7I6ilwBcYyQ8AKKX8gI+APwG/AiYC+5VSnT3ZSVFP9my9OKuVJ9PS8bUnTRzz9+PvzpUm/u16ewIhhGhN3AlQI4Hj1bL5pgE9gQPAWGAlEAFInmVLsReV7V9Syt+rVT5fERluvJHMPiGEF3AnQHUBqtdmH49RPWKO1nqf1noRcBy42kP9E+5yWsQ72VLA9NzKyudrI8LZHmJPgZXMPiFEK+dOgIoEMqodGwIc01r/7HTsK5ymAUULWPhZRZC6OyubYYVFFaceiI7iC8nsE0J4AXcCVBHQwfFGKdUVY1R1qFq7EsAf0bLslc/NwGNpGVzstD3HnbHRHPe1FxGRzD4hRCvlToA6CgxzyuKbTmVxWGddgFQP9E00lr3yebDWrExNr8jsyzOZWBAXQ7ZjD6ltc1uqh0IIUSt3AtSLQDDwuVIqAXgQsAA7HA2UUv7Ar4HvPdlJ0QiT1wLGdvErUtMJsG9adtrXlzvioilRAFoy+4QQrY47AWo18DJG6aLrgFLgFq11rlObazCC2AGP9VA0Tr+pFZl9l5aWsjw9E2VPP/8qIIC/RdvTz0tyJUiJFtGjRw+UUnzwwQct3RWXTpw4gVKKHj161Djn6PuJEyca9RnLli1DKcWyZcsadZ+2pt4BSmtt01rfhFH89bdAZ6119YcXPwPXAxs910XRaOOfqCiHNKqwiL861ezbFxLMI1GRlYVlJUgJ0Wp88MEHKKUYMWJES3elRdQaoJRSNyilamxzqbU+rrX+VGud5+LcYa31Vq11iqc7KhrJXg4J4M95+VXSz18OD2WdYzdeWSMlRL299957JCUl0blz42oTLFy4kKSkJNmqo5q6RlAvA2lKqd1KqdlKqZjm6pRoIvZySAq4JyubqyyV9YCfjoqoukZKtowX4rwuvPBC+vTp0+it06Ojo+nTpw/R0dEe6lnbUFeAugf4GmPR7XPAOaXUfqXU7UopmQfyVktPgX84PsC/0jO5zL4bL8Cy6Cg+DAww3mQclSDlrY4kwJO/gGURxqsXLiOYOXMmSik2bNjAd999x5QpU4iJiSEkJIRhw4axf//+ira7d+9m+PDhhIeHExYWxoQJE/jhhx9q3NN5uqygoIAlS5ZwwQUX4O/vT9euXVm0aBGZmZlu9bOuZ1BaaxISErj66quJjY3Fz8+Pzp07M2rUKFasWFGlratnUCNGjGDkyJEAHDhwAKVUxVd7mfKra7uNR7XWQzDSxm/HSCcfCjwFHFdKJSqlliqlZD7I29iDlB/wVGo6ve1rpBzVz4/4+xntJEh5nyMJsOt2yD0NaON11+1eGaQAEhMTueyyyzh27BijRo2id+/eHDp0iKuuuoqDBw/yzDPPMHHiRLTWXHXVVURFRbFr1y6uuOKKWoNNaWlpRZD4xS9+wTXXXENxcTErVqxgyJAhpKY2fpVMaWkp1157LdOmTeOdd97h4osv5rrrrqNPnz58++23LFq06Lz3GDt2LFdddRUAcXFxzJgxo+Jr7Nixje6jN6h1y3cHrXUyRo29lUqpSIyCsJOB0Rgp5Q8rpY4BW4HtWuvEJuyv8JSlp2BZOKFaszo1jT91jOesr5kiHx8WxMWwITmVC8vKjSC1cYLxDEt43rLwpv+MsiLYdovx5UnLcs/fppFWrlzJ448/zl13VVY8ueeee/jPf/7DnDlzSElJ4YMPPuB3v/sdAMXFxYwZM4aDBw+yatUq7rvvvhr3/OSTT7j44ov5/vvvK54d5efnM2nSJN577z0WLVpEQkLjAvrdd9/Nzp07ufjii9mxYwd9+lT+Hm+1WnnzzTfPe48lS5YwePBg9u7dS58+fdiwYUOj+uSN3NoPSmudrbXeoLWeAMQANwCvAZ2AvwGfKaVOKqWeVEoNV8qxx4NolexrpGKsNtakpBFprdyiY258LKcd+0hJSSTRQoYMGVIlOIHxgxvg2LFjLFiwoCI4AQQEBPCXv/wFoMo0YHWPP/54lcSG0NBQ1qxZg8lkYuvWrZw+fbrBfU5LS2P16tX4+Piwbdu2KsEJwGQyMWGC1MGsjwZvWKi1LtBaJ2itb8AIVtcAG4BA4A7gfeBeT3RSNBGnNVI9ystZlZJOkH0hb5rZzC3xcaSYnDY7lCAlmpmrqazIyEg6dOhQ6/mLLroIgHPnzrm8Z0REBOPHj69xvFevXgwePBibzcaHH1YvkFN/77//PmVlZQwZMoRLL720wfcR9Zjiqw+tdSnwJvCmUsoHGA5MAtI8cX/RhMY/YbwmruMXpaWsSE1nXlwMJT4+nPU1c0t8LBuSU+lgsxlBqttgI7AJz/D0NJnjGVRZZYFgfAPhmv965f+3Ll26uDweEhJCZmamy/MhISGAMd3niqsFt87nDh06xJkzZ9zvrN3JkycBaoychPs8vuW7fUHvfq317Vrr5zx9f9EExj9RMZL6TXEJT6RlYLZXmzjh58ut8bHk+thna6VuX+vWb6oRjMK7Asp49dLgBODjU/ePqPOdbyh5OtE6NGgEpZSKx3juFFBbG631xw3tlGgB45+AzB/h+AGuKCrmkbQM7o6NxqYU3/v7MT8uludS0gjW2niwP3mt1/7Qa/P6TZX/N3WoqyyR41ynTp0afP/u3bsD8P33UpK0sdz69UMpdb1SKgk4i7H9+8Favho+gStazoydFftIjSks4kGnHXmPBPhze1wMxY7fLLfdIs+khFfKycnhrbfeqnH8559/5tNPP0UpxRVXXNHg+1955ZX4+vry8ccfk5SU1Jiu4udnLPkoLy9v1H28Vb0DlFLqBmAL0BvIA44AH9fy9YnHeyqah9NmhxMtBfzNKUh9HhjAHbHRlUFK9pISXuqvf/0rycnJFe8tFgvz5s3DarUyadIkunVreC2C2NhYbrvtNmw2G1OmTOHYsWNVzlutVnbt2lWvezkyDX/88cd2GaTcmeL7m/31DmC11rr9/Wm1Fws/M+rxWZK5Md9CgY8PT0dFAPBxUCB3xEbzdFoGAVobz6RkOkl4kSFDhmC1Wrn44ou58sor8fPz48CBA6Snp3PhhReycuXKRn/Go48+yk8//cRbb73FpZdeypAhQ+jSpQtpaWl88803pKWloe3PeevSvXt3BgwYwFdffUW/fv0YOHAg/v7+9O7dm7vvvrvR/Wzt3Jniuwj4WGv9jASndsBetw9gTm4e87MrK6A7gpQxkrI/k5KRlPASfn5+vP/++9x6660cOXKEnTt34ufnx4IFC/j000+Jj49v9Gf4+/uza9cuXnzxRa644gq+/fZbXn/9dY4ePUq/fv3cCoLbtm1j6tSpZGVl8corr7Bu3bp6LfRtC1R9ojiAUuoscEBrPb1pu+RdBg0apBMT23DxjH93MyqcA6sjwlgVGVFx6reFRZUjKZDECReSkpLo27dvS3dDYNTiGzlyJMOHD2+1e095O3f+viulvtRaD6qrjTsjqH3Ab9xoL9qCpacqtumYl1PXSApJQRdCeJQ7AeofQJhS6hGllKmpOiRaocVH6xmkZLpPCOE59U6S0FqfUkoNA3YAk5VS7wFnwNgx3EX7f3mmi6JVWHy0YrpvXo6xV6Vjuu/joEAWxMXwTGo6QVobKeinPq2sUiGEEA1Q7wBlL/y6ECNZwoSx9burB1j2X6WRANXWLD1Va5D6PDCAufGxrEpNI8ymjRR0kCAlWo0RI0bUK3NOtB7upJkvARYB5Rh1934ELE3RKdGKVQtSvpqKFPT/BfgzOz6OZ1PSiHLU7gMJUkKIBnEnQM0GCoFhWuuvm6g/whs4Bak5uXkE2Wz8OzoKgKP+ftzcMY7nUtKIs1qNIJX5o+wnJYRwmztJEp2BDyU4CaBiV16A6fkWHkzPxMc+ffKzny8zO8Zy1nk/qY2y/40Qwj3uBKizGCMoIQxOKeiTLAU8kp5ZUQX9jK8vMzrGcdzXPkiXICWEcJM7AepVYLhSKripOiO8kFMK+tiCQp5MzcDPZgSpVLOZmR3j+M5e8JLjB2DF5S3VUyGEl3EnQD0EHAN2KqUubKL+CG/kFKRGFBWxMjWNQPvOvFkmE7M6xvJxgH1nloyjRp0/IYQ4D3cC1E6MDL6RQJJS6qhS6l2l1D4XX3ubprui1Vp8tKIK+uDiEp5LSSPMagWg0MeHBfExvB0cZLS1JBtJFkIIUQd3AtRo4Hf2/zYDFwNX2o+7+hLtjdNWHf1LStmUnEqcfYuAcqX4f7HRbA4ztuOmJFeClBCiTu6kmf++yXoh2o6FnxnPmTKOcmFZOS+dS+W2+Bh+sj+HWt4higyTiduzc1ElufBQPNyX0sKdFkK0Ru6UOnqvKTsi2pCFnxkZe8cPEG+1sjE5jYVxMXwd4A/A8xHhZJpM3J+RhdlaJFvICyFccmvLdyHqbcZOGDQbgHCbjedS0hheWFRx+o3QEG6Pi6FAtpBv93r06IFS6rxfLbVFhuPzqxsxYoRH+rVhwwaUUsycObNR92mL3JniE8I9jhJHiesI1JqnUtN5IDqK7aHGc6iDQYHM7BjHitR0qTohuOqqq+rcLNATGwl6kxMnTtCzZ0+6d+/OiRMnWro7LaLWAKWU+hBYorX+uKE3V0oNBf6ttb6iofcQXm78E9BtMGybixnNgxlZxJZbeS7SqEJx1N+P6Z3iWJWaTu/Sssq1Ugs/a+GOi+a2ZMkSRowY0dLdqLdNmzZRWFhIt26NS/aZNGkSgwcPJjw83EM9azvqmuLrAxxUSr2jlJqmlPKvzw2VUv5KqRuVUu8CH2Jk+4n2rN9UWJYDpkAUsCgnlwedqk6kmc38uWMcBwNlrckGeCAAACAASURBVJTwHt26daNPnz4EBQU16j7h4eH06dOHjh07eqhnbUddAeoi4BlgOPAykKqUelMp9Xel1BSl1Ail1K/tr1OUUvcppd4C0oCXMFLSn0YClHC4LwVMgYBRGml1ShqhVmNBb6GPD4viYkiwT/9hSTYy/IRworXm6quvRinF3Lk1d3C22WyMGjUKpRQLFy6sOH7ixAmUUvTo0YPy8nKWL19O3759CQgIIC4ujhkzZnDq1Cm3+nK+Z1B79+5l8uTJdOrUCT8/P+Lj4xk6dCiPPPIIRUWVz2NdPYOaOXMmPXv2BODkyZNVnsX16NHDrX56s1qn+LTWucCdSqn/YmyzMQO4Ghhbx/0UkAU8DqzSWp/wXFdFm3BfijE6siQzuLiEF5NTmB8XyzlfM1aleCg6ilO+Zu7KysFHMvxENUopXnzxRfr378/atWsZOXIkN954Y8X5Bx98kPfff58BAwbw+OOPu7zHtGnT2L17NyNGjOBXv/oVhw4dYtOmTezZs4cPP/yQ3r17N6qPWmvmz5/PmjVrABg0aBDDhw8nKyuLpKQklixZwrRp0+oMNMOGDcNisbB161aCg4O57rrrKs5FR0c3qn/e5LxJElrrn4G/KKX+hjGaGgH0B2KBcCAHY9R0GNgPHNRalzRVh0UbsPholbVSm5NTWBQXw7f+xizyxvAwTprNLE/PJNixQ+9XL7XJ5IlfbvxlS3ehwb6Z8U2LfG50dDRbtmxhxIgR3HrrrQwaNIiLLrqI/fv389BDDxEaGkpCQgL+/jWfSpw8eZKioiK++uorLrnkEgBKS0uZPXs2L730En/605/4/PPPG9W/p556ijVr1hAXF8f27dsZPHhwxTmtNR988AGRkZF13mPOnDmMHj2arVu3Eh0dzYYNGxrVJ29V7zRzrXWR1nqP1nqJ1nqs1vrXWusLtdYDtdZXa63v1Vq/K8FJ1ItT1Yloq40XktO4sqCyWP4HwUHc1CmO085bdshzqTZt5MiRtaaYR0REVGk7bNgwHnzwQfLz85k6dSqnTp1i+vTp2Gw21q5dS69evWr9nPvuu68iOAH4+fmxYsUKwsPD+eKLLzh06FCDv4fy8nL+9S9jM/ENGzZUCU5gjABHjhwpCRH1JGnmouU4LegN1Jon0jJ4OjKC9RFhAPzo58f0TvE8kZbBb4pLKmv4LXXvWYHwDnWlmbtKRFi6dCkffvghe/fupV+/fuTm5nLrrbcybdq0Oj/npptuqnEsPDyc8ePHs3nzZj744AOGDh3aoO8hMTGRjIwMunTpwtixdT0NEfUhAUq0rBk74UgCbJuLCc1d2Tn0KitjWXQUZUqRYzIxNz6WJZnZTMu3GDX82tBzqZaaJmuN3E0zdzyPuuCCC8jNzeWSSy7hqaeeqvOaiIiIGqMxB8czoTNnztS7D9WdPHkSoNHPsYRBKkmIlueUhg4wwVLA+uRUosuNaujlSvFwdBQPdYikzHGNVJ4QwPbt27FYLIARWM6ePdvoe7qqGiFahgQo0Xrcl1Kxr9SvSkp55VwKl5RUPtJMCAtlbnws2T72v7aJ62SX3nbs22+/5Y477sDPz48//vGP5OXlMW3aNEpLS2u9Jicnh9zcXJfnHNUaOnXq1OA+de/eHYDvv/++wfcQlSRAidZl8VHoORygotDs1ZaCitOJgQHc0Cme7/x8jQOSPNEuFRQUMHXqVIqKinjkkUfYtGkTI0eO5Msvv+Tuu++u89rNmzfXOJabm8vu3bsBGlXNYuDAgURHR3PmzBn27m3ctnh+9h0Ayu1b1rRHEqBE6+NUaDZAax5Jz+SOrByUvfLEOV8zf+4YzxshwUZ7WdTb7ixYsICkpCQmTJjAnXfeiY+PD5s3byY2Npb//ve/bN++vdZrH3zwQZKSkirel5WVcccdd5Cbm8vAgQMZNmxYg/vl6+vL0qVLAbj55ptrpKw70sxrG8U5i4mJwc/Pj9TUVLKzsxvcJ28mSRKidXKq4afQzMnNo1dpGX+L6UC+yYdSH8X9MR34n78/S7Oy8Hcs6h00u7JIrfAqy5cvr3O9z/Tp0xkzZgybNm1i48aNdO3alfXr11ec79ixIy+++CJjx45l1qxZDBgwoGLKzaFbt24MHDiQ/v37c+WVVxIeHs4nn3zCqVOniI6OZtOmTY3+Pv7yl7+QlJTE888/z+DBgxk0aBC9evUiKyuL//u//+P06dMcP378vKnmvr6+jBs3jjfeeIMBAwYwdOhQAgMDiY6OZvny5Y3upzeQACVar35Tja9/d4OSXEYUFfHKuRTujIvmR/v0x9awEI76+/JkagYdHRXRj75lTBUKr3K+KbH+/fvTrVs35s+fj9ls5pVXXiEqKqpKmzFjxnDPPfewfPlybrjhBj788EN8fX0rziulSEhIYPny5bz44oucPHmSsLAwbrrpJh566CGPlBFSSrF27VomTpzImjVr+Pzzz/n666+JiorioosuYtGiRfWuzL527VqioqLYu3cvCQkJlJeX071793YToJS2T5uct6FSc4DNWuui8zZuRwYNGqQTExNbuhttn708EkChUjwQHcVbjik+IMJq5T9pGQwptidVmAJbxU69SUlJ9O3bt6W70e7J1hXNw52/70qpL7XWg+pq484zqOeAM0qpx5VSF7lxnRCN55Q8EaQ1y9MzWZKZVVERPcdk4rb4WJ4PD8MG4JjyO5LQcn0WQjSKOwFqNxAG/AVIUkrtUUpdo2TRgGguM3YaC3RRKOCPeRZeSE4lxp7lZFOKp6MiWBQXU5mKvu0Wo+6fEMLruFOLbwJwAbAcyADGANuB40qpJUqpmIZ2Qinlq5QaZR+dfaqUSlZKlSqlziqlXldKjTjP9dOVUgeVUrlKKYtSKlEptUApVef319DrRAuqtqh3QEkpCedS+HVxcUWTD4MCub5zPIcdxUIzjhrPsYQQXsWtH8Ra69Na678BXYE/AZ8C3YB/AqeUUi8qpYY0oB/DgXeBu4DuwJfAGxhbd0wB9iulHnR1oVJqJbAZGAQcBN7B2INqBfC6UsrkyetEK3FfSpVis88npzEzJ6/idKrZzKyOTlN+jhJJUn2iXerRowdaa3n+5GUaNFLQWpdprTdrrYcCA4B1QDkwHfhIKfWlUmpWfXfhBWzAVuAKrXVHrfV4rfU0rfUvgRsAK3CfUmqk80VKqSnAfCAF6Ge/bhLGZotJwCRgIdU09DrRyiz8rGLKzxf4a3YOK1LSCLcaJZKs9im/+XExZDpXn5CFvUJ4hUZPZWmt/wc8AKzH2LBQYQSttcAJpdTsetzjfa31dVrrgy7OvQpssL+tXoZ4qf31Hq31D07XpALz7G+XuJiya+h1orWpNuU3vKiY18+mMMBpyu+QfcrviwD770uWZEmgEMILNOoHsFJqtFJqG3AcWAAUAy8ANwJvYWxq+JxS6vZG9vMr+2sXp8/uAgwESoHXql+gtT4AnAXigcGNvU60ck51/OKtVtYlpzE7p3K1frrZzJz4WNZEhGF1HNx2i9TyE6IVcztAKaXClVJ3KqWOAnuBa4FzwN+ALlrrOVrrV7XW1wC/BQqAxgYoR1p7stOxAfbX7+pYm/VFtbaNuU60dk6p6L7Andm5rE5JI9I+5WdTipWREcyOjyXZ1HwbIdZ3raEQ3qwp/p7XO0AppX6tlHoeY3TxOEZCwQGMJIYLtNaPaK2znK/RWn8GvImRSNEgSql4YKb97VanUz3tryfruNyxs11Pp2MNvU54A6dUdIBhRcW8djaFgUWVU35fBgYwpXNH9gTbN8Frwik/s9lcZ3VtIdqKsrIyTCbP5pa5M4JKBGbZ//t5jOSCK7XWb2itbXVcV0ADSyoppczAS0A48J7WepfT6RCn+9fGYn8N9cB1wls4nkvZs/zirFaeT0ljXnYuPvbf8vJNPtwdG8290VFYHEv5mmDNVHh4OJmZmTKKEm1eXl4eoaGe/ZHpToA6AdyNMY13q9b623pedwvGjEtDrAFGAaepmSDhWCDs7r/8hl5XeQOl5trXTCWmp6c39DaiqS38rKIquhmYn5PLxuRUOpdVbl+wMzSE6zvH87W/UduPjKMeHU1FRUVRUlLCmTNnyM/Px2q1SrASbYbWmtLSUjIyMsjOzq5RG7Gx3BnZXKgb8C/Lfo31vA2rUUo9DczGSAUfpbWuXlgt3/4aQu0c5/KdjjX0ugpa6+cwSj8xaNAg+WnTmjmqor9xK2gb/UtKef1sMv/qEMWuUKOW3xlfX2Z2jOPWnFxuyckz/lFsuwW+esmYMmwEs9lM9+7dyc7OJjs7m3PnzmGz1TXhIIR3MZlMhIaG0q1bN/z967uyqH7cCVB7lVJ7tNZ17mWglPoLcLXWekxDO6WUehwjsSIdIzj94KLZCftrdxfnHLpWa9uY64S3clRFX3E5ZBwlRGv+lZHJsKIiHu4QRb7JB6tSrIqM4OPAQP6VnkHXcmtlAkUjK6P7+PjQoUMHOnTo4KFvSIj2wZ0pvtHAL+rR7hKMabkGUUr9B6OiRCbwe631/9XS1JF6fqlSKrCWNr+p1rYx1wlv5zTlB/CHgkJeP5dcpUzS1wH+TOnckYTQEGMO2JFAIRUohGh2TbEQ1Q+jMoTblFLLMZ5zZWMEp//V1lZrfRo4bP+8613cazjGuqkU4JPGXifaiPFPwLLcijVTncqtvJCcxu1ZORWV0Yt8fHgoOopb42NIcWQlJa6Ten5CNDOPBih7ZfOBGMVk3b32IeAeIAcjONVn9PJv++sjSqleTveKBVbZ3y53kWXY0OtEW+G0ZsoE3JKbx0vnUrmgtKyiySeBgUzq0pHtIcHGaErq+QnRrOrcsFAptc/p7WiMBbm1TbmZMRbUdgJe11pPq3cnlJoA7LC/TQS+q6XpUa11la0klVKrMMoTFWMUnC3DmGIMw6i2fp3WukaSRkOvq042LPRyRxIqEigAShSsiIhgY3go2mknmREFhfwjM4toq/13Fv9wWHrK1R2FEPVQnw0LzxegnEcQmsoU7bocASZqretaCFv9c2Zi1PI7nwNa6xEurp+OUWrplxi/EB/FKLm0uq5RUEOvcyYBqo2wJ1A4fOXvx99jOnDKabvwcKuVv2dmM7agsPK6QbONaUMhhFs8EaAcyQ4K2IdR2uixWpqXAme11j83oK9eSwJUG7L7LuNZk12hUjwVGcEr4VUXH46xFLA0M5voinRxBZOfMzIFhRD10ugAVe1mB4E3q0+xtXcSoNqgaqOpzwL8uS+mA8nmylUZ4VYr92RlM95SWDmt0HN4o9dNCdFe1CdAubOj7u8kOIl2wWmfKYDLi0vYdiaZyfmWiia5JhN/i4lmXlwM58xOhWeXRcg2HkJ4iOx3JIQr1er5hWjNAxlZPJuSVqVU0qGgQK7t3JHNYSH2tRVatvEQwkNqneJTSv3N/p+rtdbZTu/rRWv9r8Z2zhvIFF87cCQBts3FUb6xUCmeiQxnc1jVTL/+xSU8kJHJBRUBTJ5NCVGbRj2DsmfwaaCv1vqY0/vzfi5GCT7P1l1vpSRAtSPVnk197e/HsugofvLzqzjmqzW35uQyKyevskJydB9j2lAIUaGxAephjID0pNY6y+l9vWit73Ons95KAlQ7U200VQqsjQjn+Ygwyp1GU71KS7k/I4sBJU57QUlKuhAVPJrFJ1yTANVOVRtNHfP15R8xUXxbrZrzlHwLd2blECEp6UJU4dEsPiGEk2qZfheXlfHSuVQWZ2YT6LSdxtbQECZ06cgOR7kkSaIQot4kQAnRUNUy/UzAjLx8dpxJZqRTtYlsk4m/x3RgVnwsP/va11IdPyB1/YQ4j3oHKKXUPKVUqVJqXB1txtvbzPFM94TwAo7RlDL+OXW0WvlvWgZPp6YTX16Zkp4YGMCUzh35b2Q4xY7nVYnrjD2nhBA1uDOCmgxkAW/X0eZte5vrGtMpIbxOv6nwj+yKCukAVxYWseNMMjfn5GGyP+stV4q1EeFM6hzPwcAAo6FjzymZ9hOiCncCVB/gm7qKqNqrf3+DsWmhEO3PjJ3GaMrHmMoL0pq7snN49WwK/YtLKpqd8fVlfnwsi2KjOV2lEkW4VKIQws6dABUDpNajXRoQ27DuCNEG9JsK92dWSaLoXVbGxuRUlqVnEmat3MXlg+Agru3ciWciwilyTPttuwWWRUmgEu2eOwEqF+haj3adAct5WwnR1lVLovABplgK2HkmmUlOdf1KfRTPRYYzoUtH9gYF2rP9rJLtJ9o9dwLUV8BgpdSFtTWwn/st8HVjOyZEm1EtiaKDzcaDGVlsPpfCL0oqp/1SzGYWx8VwS3wsPzr2oZJsP9GOuROgNgC+wHal1EXVT9q3Tt+OkW27wROdE6LNcCRRDJpdeaiklM3nUnkgPZMop2m/zwIDuK5zPI9ERZDn45TtJ8+nRDvjzn5QCtgF/AEoBz7C2IEWoDfwO4xt3/dorf/g+a62TlJJQjTIxgnG6Mguz0exKiKCLWEhWJ1KJkVZrSzKzmFSfgEVxS2ltp9oAzxe6kgp5Qc8CdyCEYyclQNrgbu01iXVr22rJECJRvl3NyjJrXh7zNeX5R0i+cKRgm53UWkpd2dmM8QpE1Bq+wlv1mS1+JRS8cAooLv90EngPa11its383ISoESjVdtqXgN7g4N4LCqCVHPV3wOvKCzir1nZTlt6IIFKeCUpFtsMJEAJj6lWgLZIKTaGh/JCeBhFPpWPi01aMzXPwrycXCKd6v4xea0UoRVeQ4rFCuFNHNl+9kW+gVpzW04eu88kc22+BWX/ZdKqFK+EhzKuSyc2hoVS5rhe1k+JNsbtEZRSqjdwOzACY80TwFlgP7BCa320lkvbJBlBiSZRbd8pgCQ/Xx6Nqvl8qmtZGX/NyuHKwiIq0iskkUK0ck2RJDETWA34QeW/BSelwK1a641u9NOrSYASTcrF86n9QYE8ERXBSV/fKk0HFBfzl6ycqpsk9hxulF8SopXxaIBSSv0G+BhjWvAN4AXgJ4xA1ROYhVFQ1goM1Vp/0fCuew8JUKJZPNbHKCprVwa8GhbK6ogw8kymKk1HFBRyZ3YOF0oihWjFPB2gEoApwE1a61dqaXMjsBl4TWs9zc3+eiUJUKLZHEmA7fPAVhl4cn18WBMRxpaw0CpbzvtozURLAfOzc4l3WgQsgUq0Fp4OUOeAM1rry87T7jOgm9a6Y7176sUkQIlm5+L51BmziRWREbwVHIR2ClT+NhvT8yzMzs0l3OZoL9vOi5bn6Sy+DsCxerT7AYhy475CCHc4itA6lU3qUm5leXomCedSGFpYVHG8xMeH9RFhXN2lMy+Eh9o3SrRvOy81/kQr506AygZqLRTr5AJ7WyFEUxr/BCzLrbJJYp/SMtakpvN8ciqXOhWizTf58GRUJOO6dGRrSHBlanriOlgWIanpolVyJ0B9DFymlJpYWwOl1DXAYOBQYzsmhKinGTtrBKrLi0t45Vwqj6Wm062sIhyRZjazLKYDE7t0ZFdwEMbTKS1rqESr5M4zqGHAAYwsvZeAjcBxjInwC4A/AzdhVDMfrrVuF0FKnkGJVsdFxt+20BBWR4STaa6a8XdBaRkLsnMYXVjk9NuqCSavkWdUokk1xTqoRcATuB55KYzg9Ret9Qp3OurNJECJVulIArxxK+jKUkiFSvFyWCjrw0NrpKb3KSllQXYuw4ucFvuGdITF7WrdvWhGTVKLTyk1ALgTuALohBGYzmKMrp7WWn/VsO56JwlQolWrttAXIF8pXgwPY1N4KAU+VX/X/GVxCQuzcxlSXFwZqHz84NqVMqISHiXFYpuBBCjhFVwEqhwfH9aHh/JyWCjF1QLVwKJiFubkMsh5ew8ZUQkPkgDVDCRACa/iIlBlmHxYFx5GQmgopT5VK5hdXlTMrTm5/MY5UMmISniABKhmIAFKeCUXgSrFZGJtRBjbQkOqVKUAY0Q1LyeXy4pLpCCt8IhGBSil1HON+Gyttb61Edd7DQlQwqtV23oejKoUz0aEsyskuMr282AUpL0tJ48hRfKMSjROYwOUzeWJ+tFaa9P5m3k/CVCiTXAxojptNrEuIpwdIcE1RlT9iku4NSeX3zkHKklPF25obICa7fJEPWmt152/lfeTACXaFBcjqnNmE+vCXU/9XVJSwm05eYxw3otKav2JepBnUM1AApRok2p5RvVCeBhbQ0NqJFP0LinllpxcRhcWUWXqRKqni1pIgGoGEqBEm+ZiRJVmMrE+PJTXQkMoqZae3qO0jJtz8xhvKcDP+YQEKlFNkwUopVQIMAiIAU5prdttKo8EKNEuuAhUGSYfNoaF8WpYCEXVAlVseTl/zs3nunwLwc4/YyRQCbumKHUUCjyOUXfPsd/0Rq31LPv5ecBS4Dqt9ecN6rWXkQAl2hUXgSrTx4fN4aFsCQ0l31Q1UIVZrdyYZ2F6Xj5RNqe8K9mKvt3z6H5QSqkg4ANgDpAHvAOoas32AV2ASW71VAjhHVxUTu9gs3F7di77Tp/lrqxsossrd/DNM5l4NjKcq7p2YnlUJMmOGoDHDxj7Ua24vLm/A+FF3Nlu46/AAOAVoKfWemz1BlrrnzA2LLzSM90TQrRKjkDltGliiNbcnJvPnjNnuT8jk65O23wU20dZf+jaiXujo/jJ12ycyDhqBKoHY2SrD1GDO9ttfIOxU+6FWuti+zEbsMExxWc/tg+4RGvdpQn62+rIFJ8QuMz6swLvBAexLjyMo/5+NS4ZUVDIn/PyGeRcnULWUrUbnt7y/ULgc0dwqkMGEO3GfYUQ3s6xu+/ktaCMHysmYGxBIQnnUliTksagoqo/Oj4IDmJWxzhu6BTHW8FB9l1+rbIdvajgToAqA/zr0a4LYGlYd4QQXq3fVPhHdpVApYChRcWsT0njxXMpjCgorHLJ//n7c09sNH/o2omNYaFYHIuBE9cZgWrjhGb+JkRr4U6AOgYMUErVGqSUUhHAr4BvG9sxIYQXcwSqQVUL0vQvKeWZtAx2njnH9Xn5+Dtl9qWYzTzWIZLR3TrzaFREzYQKCVTtjjsBaisQB/yrjjYPAyHAa43plBCijXBM/VXL/OtZVs79mdnsO32O+dk5RFkrM/8KfHzYFB7G1V078f9iOvCdn31FiyNQyfRfu+FOkkQwkAhcDHyEEbCeAvYDW4DrgVHAd8BvtNYltdyqTZEkCSHc5GItVbFS7A4JYlNYGMcdAcnJoKJi/pSXz3AppdRmNMVC3a4YgWkQoDGmlx03UMDXwESt9ekG9dgLSYASooFcBCob8FFgABvDw/g8MKDGJZ3LyrkxL59JFgthNqefXbI3lddpylJH44E/ABdgJOucBt4GtmqtG7NNh9eRACVEI7lIUQf4Pz9fNoWHsTc4qEYV9UCbjQmWAqbn5XNBWXnlCdmbymtIsdhmIAFKCA+pJVClmEy8EhbC1tAQck01t5kbUlTEH3Pz+V1RcdWH6jL916o1dj+o14F1wB4tUaxWEqCEaAIrLjeqTDgpUoq3goPYHB7KD341F/52LSvjxjwL1+ZbCNUy/dfaeWJHXQ0kAxsxKkb84PFeejkJUEI0IRfPqTSQGODP5rBQ9gcFYqs2/RdkszExv4Ab8/Pp6Tz9BzKqakUaG6BWAtMwyhs5Gn0EvAC8prUudHlhOyMBSohmUMv031mziVdDQ3k9NKRGJXUwpv+m5VkYXliE2fmEVFNvcY1+BqWU8gMmArOA0RgJERooAF4F1mutP/ZYj72QBCghmtGRBNixEKxVV7EUKsWbIUG8HBbKjy6m/2LLy7ku38KU/AJindZcSe2/luPRJAmlVCdgBsZeUL3thzVGhYkXgBe11ikN7653kgAlRAupZfrvM/v034GgQHS16T+T1lxZWMS0vHwuq1KkFhlVNbOmTDMfgjGquh4Iw/h7YcVINX8B2K21ttZ+h7ZDApQQLcxFoAI4ZzbxeqiR/ZflIvuvR2kZU/MtTLBYCHdeUyWp6s2iydPMlVKBwHXATGCE06l0rXV8g2/sRSRACdFK1DL9Vwa8GxzEq6EhfOli8W+AzcbVBYVMy7NwaWlp1ZOSVNFkmnUdlFLq98BLQAygtdY1f2VpgyRACdEK1TKq+tHXl1fDQtgVEkyBT82kiktLSpiWZ+GqgkKCJFW9STXHCCoEI9NvJvBbKreAP6W17tHgG3sRCVBCtGK1ZP85kipeDQ3lexebKQbbbFxtKeC6/AIuKS2t+qxKRlUe0ZTPoEYCNwOTgUCMwFQC7MR4BrWvvSzulQAlhJeoJanif/5+JISGsickiLJqSRUAfUpKmZxvYVxBQdX6f/KsqlE8ncXXEyOLbwbQjcrR0tfAeuAlrXV2w7vrnSRACeFlahlVZfv4sCMkmK2hIZxwUVE9wGZjTEEhk/ML+HVJiYyqGskT66CCMDL1ZgK/wwhKCsgGXgbWaa2/9lSHvZEEKCG8VC1JFRo47O/P1tAQ9gUHUuLiWVXP0jKm5Fu4xlJAlK1afWwJVvXS2EoS6zCCUzBGULIB72FM4b2htS51eWE7IwFKiDagllFVro/irWBjVOXqWZVZa64sKGRKfgGDi6sVq5V1VXXyRC0+gOPABoyqEWc82sM2QAKUEG1IHaOq7/z82BoazFshwRS6GFV1LC9nQn4BEy0FdC2XGoDn09gA9SLwgtZ6f1N0rq2QACVEG1VLqnqhUuwNDuL10BCOBPi7vHRgUTETLQU109WltFIF2Q+qGUiAEqKNO5IA2+eBrbzGqR98fdkWGsyukGCXe1UF2mz8vqCQay0FDCwukSlAJxKgmoEEKCHakVqeVZUCB4IC2REawkeBAVhdpKt3LitnosVIrOhSXq0SXDucApQA1QwkQAnRDtXyrAog3eTDm8HBbA8N5icXldUBLrNPAY6uPgXYjtZWSYBqBhKghGjnahlVORIrtocG81ZwsMv9qoLsa6vGWwr4TfUpwDZeXkkCVDOQACWES1gDFQAAEF5JREFUqFBLYkWJgv1BQewICebjwIAauwADxJWX8weLEawuLiurerINPq+SANUMJEAJIWqoYwow1WRid0gQ20NcV6wA6F1SynhLAVcXFBJnbZvPqyRANQMJUEKIOtUxBfiNvx+7g4PZExJEtossQKU1lxWXcI39eVVw9Z/XXhysJEA1AwlQQoh6W3E5ZBytcbgM+DgwgN0hwewPcl1eKcBmY2RhEeMtBQwpKqbK2MsLkyskQDUDCVBCCLfVsbbKohTvBgexOySYzwP8a2xbDxBltTLWUsi4ggJ+WVJtOxAvSa6QANUMJEAJIRqllilAMJ5XvRUSxO7gYI65qAUIxvqqqwsKGGsp5OKyMq8JVhKgmoEEKCGEx9SSBQjwva8vb4YE82ZIEGlms8s2F5SWMbaggKsthfSoXg+wlWUCSoBqBhKghBAeV0cWoBVIDPDnzZBg3g0Kcrm+CqBvSWnFyKpj9UzAVhCsJEA1AwlQQogmVcfzqlLgUFAgbwcH8UFQIEUukisABhQXM9ZSyJiCQqKr71/VQsFKAlQzkAAlhGg2dQSrQqX40B6sDgYFuty+3kdrLisu5mpLIaMKCwm3tVzaugSoelBKTQfmAf0AE3AUYwv71VprW13XggQoIUQLqSO5Il8p3gsOYk9wEJ/WUrzWrDVDior5fUEhVxYWEd7MIysJUOehlFoJzAeKMXYLLgNGAaHAG8D1Wmtr7XeQACWEaAXqSK7I8vHhneAg3g4O4nAtaetmrbm8qJgxBYWMLCwishmClQSoOiilpgCvAynAFVrrH+zH44D9QF/gTq3103XdRwKUEKJVqSNYpZhM7AsO4u2QIL71d73Zosk+Dfj7gkJGFRQR1UTBSgJUHZRSicBAYIbWelO1c8OBDzCCV+e6pvokQAkhWq1aKlcAnDGbeDcoiH3BQXxTy87APlrzm+ISe7DybIKFBKhaKKW6AKcxkmAitNZFLtqcAToDQ7XWH9d2LwlQQgivUMfI6pzZxDtBQbwTHMT/aglWSmsGFpcwpqCQ0YWFxFidgtXktW6XWapPgHK92qvtG2B//c5VcLL7AiNADQBqDVBCCOEVnEc61YJVp3IrM/LymZGXT4rJxLvBQewLDuSrgICKNlopEgMDSAwM4N86kgEl9mBVUETctluMRh6uBdheA1RP++vJOtqcqtZWCCHaBkewcrEgON5q5aa8fG7KyyfVZOLd4EDeCaqaYKGV4nBAAIcDAsjzyWFeTh6896AEKA8Jsb8W1NHGYn8NrX5CKTUXmAvQrVs3z/ZMiP/f3v0Hy1nVdxx/fxLzgxhEkAjjBGj4OaOtBTG10moQageqVECkTLHFWutMw7S0KtSplBmqYMDq6BQhUqHQUsYf/BLbwrRCozL8MFctjQg1/IiECB0QEu4NUGry7R/nrHmy7N3s3j27++zdz2vmmZN9znPPc55v9u73Pj/2HLNBef2pO5JKi2S1z7ZtnP7sFKc/O8WTc+dwW75n9d2FC34+6eLbt+aLUFseK969cU1QjecsZ3QDLiIuBy6HdA+qVKfMzIZmF8lqybbtnDY5xWmTUzw1Zw63v3wR6xbM5+DG7L97LC3epXFNUJO5XNxmm0bdZJttzMxmn2qygpc8Dbj39u2cOjnFqdVPx2PPK96NcU1QG3J5QJtt9mva1sxsPFWn7HjJo+tz4eTVfZkscVwT1Pdz+TpJu03zJN/ypm3NzGyA80u1Hvp2louIjcD3gPnAe5rr8xd1l5K+qHvXYHtnZmYwpgkq+2QuL5J0cGOlpFcDl+aXqzoZMNbMzMob10t8RMR1ki4jjWS+TtI32DFY7CuAm4BLhthFM7OxNrYJCiAiVkq6AzgTWMGO6TaupMPpNszMrD/GOkEBRMS1wLXD7oeZme1sLAeLLUnSk7QfMqmdvYGnCnZnHDhm3XG8uuN4daeXeB0QEUvabeAENUSSJnY1mq/tzDHrjuPVHcerO/2O1zg/xWdmZjXmBGVmZrXkBDVclw+7AyPIMeuO49Udx6s7fY2X70GZmVkt+QzKzMxqyQnKzMxqyQmqEEm/K+nbkrZImpI0IelMSR3HWNI8ScdK+rSkuyU9LulFSZskXSfp6D4ewsCViFmbti+UFHn5SIn+DlvpeEnaTdI5ktZK2izpOUmPSPqqpF8r3f9BKxkvSUsl/a2k/5b0vKQXJK2XtFrSgf3o/6BIOkzSWZKukfSApO359+aUHtvtPf4R4aXHBfg8aXbe54F/Bm4Ens3rbgDmdtjOb+SfCeDx3NaXgXWV9X897OOtU8ymaXs58DNge27vI8M+3rrFC1gGrM8//z/A14CvAN8BXgTOHfYx1yVewBHAM/lnN5LG6bwJeCyvmwSOGvYx9xCrz1Y+X6rLKcOO/9CDM+oL8O5KQjmksn4f4Ie57qwO2zoGuA54S4u638kfugG8bdjHXZeYtWh7AXAfsCn/Uox8giodL+DlwIONP3iAeU31rwIOHfZx1yhed+afubwaK2AecEWuu3fYx91DvD4AXAycChwErOklQRX9TBx2cEZ9ASZywH+/Rd2Kyn/UnAL7+mJu74phH3ddYwZclH/+BOCqWZKgisaLNNVMAFcP+9jqHi9gITvOKPZtUf+aSv2iYR97ofj1mqCKxd/3oHogaSlwJOmSyFeb6yPim6S/5PcFfrXALhuz+y4t0NZQ9DNmkt4EfBi4NiK+3ntvh690vCTNB/4ov1xVrqf10If31zbSlQsAtahvfE9nK+ly1lgrHX8nqN4ckcv7ovW08QBrm7btxSG5fLxAW8PSl5hJWghcDTwNnDXz7tVO6XgdSbqEtzEi7pd0VH6g5AuSzpf05l47PGRF4xUR/wfcll+eL2leoy7/+xP55RWRTxHGXNH4j/10Gz1alst2o5k/2rTtjEjaF3hffnl9L20NWb9idgFwGHBaRMym0ahLx+uXcrle0lXAGU3150m6Hvi9Nh8wddaP99dK4FbSmefxkiby+uXAnsDngLO77OdsVTT+TlC9WZzLrW22mcrl7jPdiaSXAdcAewC3jfjlq+Ixk3QU8GfATRHx5R76Vkel47VXLt9KmqDzb4DVwE/zuktJN7mfBd7fbWdroPj7KyIezu+xfwCOZ+dL7BPAt/KZlhWOvy/x9aZxTbrfp/arSVPRbwTe2+d99VvRmEnaDfh70gfqyhJt1kzp91jjd/5lpMtSZ0fEQxGxOSJuBk7M+zpjRL/fU/x3MienHwAHA+8izYG0hBSrPYHrJZ1Xan8jrmj8naB6M5nLxW22adRNttlmWpI+B/wh8ARwbEQ8MZN2aqR0zC4EDgU+FBGjfG9uOqXjVd3m75orI2IC+C7ps+HoDtqrm6LxkvRK0needgeOi4ibI+KnEfFURHwNOI70cMRfSTqkXVtjomj8naB6syGXB7TZZr+mbTsm6dPAnwJPkpLT+m7bqKENuSwVs5NIX8g9Q9Ka6kL68AD447zuizPo77BtyGWpeFW3eWSabRrr9+2gvbrZkMtS8XoH6Wzp7oh4uLkyIh4E7iGdkR7daSdnsQ25LBJ/34PqTeOx79dJ2m2am8rLm7btiKSLgQ+R7g28PSJ+OPNu1ko/YjaH9P2K6RyYl1d22F6dlI7X9yr/fhXpj59me+dyqkVd3ZWO1/653NJmm8253KvNNuOiaPx9BtWDiNhI+oWfD7ynuV7SCtIN1SeAuzptV9Iq0lNBz5CS071FOlwDpWMWEb8QEWq1kB47Bzg7rzu83JEMRh/itYn0Fz+k+5rN7e0JvCG/nGiur7s+/E7+JJdHVh8xr7Q3j/ToPkx/Rjo2isd/2N9aHvUFOIUd34w+uLL+1aQhd14yrAfpm/wPAJ9s0d7H8888Axw57OMbhZi12c9VzI6RJEq/x05gxxh8h1fWLwS+lOsmyPPFjdpSMl75Z7bmn7kEWFCpWwBcluueBvYY9rEXit8adjGSxC7eX13Hf9r9DDsYs2EhPZobpJulXycNhrglr7uRpoERKx+cVzWt/212DJuyNm/XavnosI+5LjHbxT5mRYLqR7yAT+X6/wW+ldvYlNc9RmUMtVFcSsaL9F2xxjiYm4Cbc5s/yeteAE4c9jH3EKs3AHdXlsagrj+qru/y/dVV/KdbfA+qgIhYKekO4EzSvZC5pL8urgQui4jtHTZVvYb9xry08k1GfJiagjEbC6XjFRFnS7oT+BPSN/oXkb5A+RlgVUS0ujc1MkrGKyKulrSO9F27twC/mas2kQaL/UyM9j3iVwBvarF+xk8lloq/p3w3M7Na8kMSZmZWS05QZmZWS05QZmZWS05QZmZWS05QZmZWS05QZmZWS05QZmZWS05QZrOApCMkhaQrB7S/4yT9m6SnJT0n6QeSPiZpwSD2b+PBCcpsdjg5lzf2e0eSzgFuAY4hDQz6L6Rx1j4BrJG0qN99sPHgkSTMZgFJ95GmhlgSES/0cT9vBL5DGmPtmIi4J69fTEpUbwU+GxF/3q8+2PjwGZTZiJN0KPBa4JZ+Jqfso6RpvS9qJCeAiJgC/oA0eeTKPBOtWU+coMwGIN8fivzv90makLRV0hOSrpC0JNctlHS+pB9JekHSo5IuaDUXUcW7c3lDZX9H532uyW1+XNKDkp6X9LCkcyXNzdvul/uwKe9znaT3tjiG+cDx+eU/NddHmnH2LtJcQL81gzCZ7cQJymyAJF0EfIE0f9CtpOkH3g98I18mu400wvh9wO2kWW//Evh8m2ZPAl4kXWJrNh/4d9Ko0v9FmutnH9K8Y5dIOog0tcvbgG/nf/8i8I+STm9q6zDSqOdPR8RD0/RlbS6PaNNfs454ug2zwTqDNEng/fDzGWzvAl6fy83AsojYkusPJ33of0DSBRHx42pjkpaSpmW5JSImW+zvzcAdTW3+cm7zg6SpEL4EfDgituX6M0mT853PzmdKy3L5aJvja9Qta7ONWUd8BmU2WOc1khNARDwDrM4vXwt8sJFIcv1/Av9Kuu+zokV7J+e6G1rUQbon1NzmvbnNOaQzonMaySlrnOEdJGn/yvrFudza5vimcrl7m23MOuIEZTZYt7ZY92Auf1xNXhXrc/maFnUnA9tIs7y2Ml2bjX3eHhEvVisi4mfAIy32qcYm0+zLrCgnKLPBeqzFuqk2ddX6hdWVkvYGfh24o80MuLtqs5t9Ni4hLmZ6jbpWlxvNuuIEZTZAu5jquttp7k8kTaU93eW9TtrsZp8bcrl/m232a9rWbMacoMxG10m5vGlA+3uA9AXdvfLTf638Si6/P5gu2WzmBGU2giTtDhwLTEREu6fqisn3qm7JL5sfQUfSgaSnBqd75N2sK05QZqPpncACBjD2XpNVpIck/kJS42ypMdTRlaTPlEsjYvOA+2WzkBOU2WhqDA7b7v5TcRGxljTc0SLgzjyi+VeAh0iPwd8DfGyQfbLZywnKbMRIWggcB9wfEQ8Mev8RcTFpyKP/AJYDJwBPAecCKyLiuUH3yWYnj2ZuNmIkvYv0YMSFEeGzFZu1PNSR2eh5njQM0TXD7ohZP/kMyszMasn3oMzMrJacoMzMrJacoMzMrJacoMzMrJacoMzMrJacoMzMrJacoMzMrJb+H71y9g1Z3aloAAAAAElFTkSuQmCC\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "##Tsiolkovksi equation\n", | |
| "m0=0.25#kg\n", | |
| "mf=0.05\n", | |
| "dmdt=0.05\n", | |
| "y0=0\n", | |
| "v0=0\n", | |
| "u=250#m/s\n", | |
| "g=9.81 #m/s^2\n", | |
| "N=1000\n", | |
| "t1=(m0-mf)/dmdt #first launch time\n", | |
| "t=np.linspace(0,t1,N)\n", | |
| "dt=t[1]-t[0]\n", | |
| "m=np.linspace(m0,mf)\n", | |
| "v=-u*np.log(m/m0)\n", | |
| "\n", | |
| "#Numerical sol arrays\n", | |
| "numsol_heun=np.zeros([N,3])\n", | |
| "numsol_rk2=np.zeros([N,3])\n", | |
| "\n", | |
| "#Initial Conditions\n", | |
| "numsol_heun[0,0]=y0\n", | |
| "numsol_heun[0,1]=v0\n", | |
| "numsol_heun[0,2]=m0\n", | |
| "numsol_rk2[0,0]=y0\n", | |
| "numsol_rk2[0,1]=v0\n", | |
| "numsol_rk2[0,2]=m0\n", | |
| "\n", | |
| "for i in range (N-1):\n", | |
| " numsol_heun[i+1]=heun_step(numsol_heun[i],simplerocket,dt)\n", | |
| " numsol_rk2[i+1]=rk2_step(numsol_rk2[i],simplerocket,dt)\n", | |
| "\n", | |
| "plt.plot(m/m0,v,label='Analytical')\n", | |
| "plt.plot(numsol_heun[:,2]/m0,numsol_heun[:,1],'-o',label='Implicit')\n", | |
| "plt.plot(numsol_rk2[:,2]/m0,numsol_rk2[:,1],'-',label='Explicit')\n", | |
| "plt.xlabel('m/m0')\n", | |
| "plt.ylabel('Velocity (m/s)')\n", | |
| "plt.legend(loc='best')" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "2. You should have a converged solution for integrating `simplerocket`. Now, create a more relastic function, `rocket` that incorporates gravity and drag and returns the velocity, $v$, the acceleration, $a$, and the mass rate change $\\frac{dm}{dt}$, as a function of the $state = [position,~velocity,~mass] = [y,~v,~m]$ using eqn (1). Where the mass rate change $\\frac{dm}{dt}$ and the propellent speed $u$ are constants. The average velocity of gun powder propellent used in firework rockets is $u=250$ m/s [3,4]. \n", | |
| "\n", | |
| "$\\frac{d~state}{dt} = f(state)$\n", | |
| "\n", | |
| "$\\left[\\begin{array}{c} v\\\\a\\\\ \\frac{dm}{dt} \\end{array}\\right] = \n", | |
| "\\left[\\begin{array}{c} v\\\\ \\frac{u}{m}\\frac{dm}{dt}-g-\\frac{c}{m}v^2 \\\\ \\frac{dm}{dt} \\end{array}\\right]$\n", | |
| "\n", | |
| "Use [two integration methods](../notebooks/03_Get_Oscillations.ipynb) to integrate the `rocket` function, one explicit method and one implicit method. Demonstrate that the solutions converge to equation (2.b) the Tsiolkovsky equation. Use an initial state of y=0 m, v=0 m/s, and m=0.25 kg. \n", | |
| "\n", | |
| "Integrate the function until mass, $m_{f}=0.05~kg$, using a mass rate change of $\\frac{dm}{dt}=0.05$ kg/s, . \n", | |
| "\n", | |
| "Compare solutions between the `simplerocket` and `rocket` integration, what is the height reached when the mass reaches $m_{f} = 0.05~kg?$\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 7, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "def rocket(state,dmdt=0.05, u=250,c=0.18e-3):\n", | |
| " '''Computes the right-hand side of the differential equation\n", | |
| " for the acceleration of a rocket, with drag, in SI units.\n", | |
| " \n", | |
| " Arguments\n", | |
| " ---------- \n", | |
| " state : array of three dependent variables [y v m]^T\n", | |
| " dmdt : mass rate change of rocket in kilograms/s default set to 0.05 kg/s\n", | |
| " u : speed of propellent expelled (default is 250 m/s)\n", | |
| " c : drag constant for a rocket set to 0.18e-3 kg/m\n", | |
| " Returns\n", | |
| " -------\n", | |
| " derivs: array of three derivatives [v (u/m*dmdt-g-c/mv^2) -dmdt]^T\n", | |
| " '''\n", | |
| " dstate = np.zeros(np.shape(state))\n", | |
| " dstate = np.array([state[1], (u/state[2])*dmdt-g-c*state[1]**2/state[2],-dmdt])\n", | |
| " return dstate" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 8, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "<matplotlib.legend.Legend at 0x7f7ebfabddd8>" | |
| ] | |
| }, | |
| "execution_count": 8, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEaCAYAAABEsMO+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOy9eXiU1fn//zozmcm+MlkISQh7WEQgqERAAljAQlFAAf3YgrgLaGvxJ9Rqqfq12BZbq6gVrYCiGAURxIrKGnYiWEATwhrCkn3fMzPn98czmez7Hs7ruuYa5jnLcyYk855zn3sRUkoUCoVCoeho6Np7AQqFQqFQ1IQSKIVCoVB0SJRAKRQKhaJDogRKoVAoFB0SJVAKhUKh6JA4tPcCOjsmk0mGhoa29zIUCoWiU/HDDz+kSSl96+qjBKqZhIaGEhMT097LUCgUik6FECKhvj7KxKdQKBSKDokSKIVCoVB0SDqsQAkhXhFCSNtjSR397hNCRAshsoUQeUKIGCHEQiFEne+tqeMUCoVC0TZ0yA9jIcRNwP8H1JmHSQixClgPjASige+A/sCbwOdCCH1LjlMoFApF29HhBEoI4QisAZKBL+voNwt4AkgChkopp0kpZwD9gFhgBrCopcYpFAqFom3pcAIFvAgMAh4Dsuvot8z2/KyU8kzZRSllMvC47eXSGkx2TR3XcpyIgn8MgeVe2vOJqFa7lUKhUHRWOpRACSFuAX4PfCyl3FpHvyAgHCgBPqvaLqXcA1wBAoBRzR3XopyIgq1PYkm9TNIP7pRcvgqbHoa101vldgqFQtFZ6TACJYRwAtYCGcBT9XQfbnv+SUpZWEufo1X6Nmdcy7HjRfISrZz/2o/MM24k/eCJlMCFPUqkFAqFogIdRqCA/wcMABZLKdPq6dvL9lxXoNelKn2bM67lyL6M3mjFXKz96POTnMi55KS1XdijzH0KhUJho0NkkhBC3Ar8Ftgspfy0AUPcbM/5dfTJsz27t8C4SgghHgEeAQgJCaljqhrwDMKZRLz75ZMZry0n+ZgnbgHF6B0lbH4chs5u3JyKDo2UktzcXHJycigoKMBisbT3khSKFsPBwQFPT098fHxwcGhZSWl3gRJCOAMfADlo3nUNGmZ7bmw54KaOq4SU8l3gXYCRI0c2bq6JL8Cmh/G9IZfcRGfMhXosxXqSf/Qk8JYssJrhq6dh2mvNWaKigyClJCUlhfz8fHx8fAgICECv1yOEqH+wQtHBkVJSUlJCeno6iYmJ9OzZE52u5QxzHcHE9wpaDNLTUsprDRyTa3t2q6NPWVtuhWtNHddyDJ0NOiN6gyQgvNxJMfuCC3lJjtqLmPdb5daKtic3N5f8/Hx69uyJl5cXDg4OSpwUXQYhBI6OjnTv3h0HBwcyMzNbdP6OIFAzACswTwixu+IDmGLr87jt2nu21xdtzz3rmDe4St/mjGtZ7loFgHtQEe7B5b4aSUc9sZbaPrz+HtZqt1e0HTk5Ofj4+KDXq9hvRddFCIGXlxf5+XWdnjSejiBQoK1jXA0Pf1t7b9vrkbbXx23Pg20mwpq4qUrf5oxrWYbOhl7jAAgIz0ZvtAJQmu9Ayknb0VfeNeXV1wUoKCjAza2uDbtC0TVwcXGhsLA25+im0e4CJaUMlVKKmh5obucAz9iuDbONSQSOAUbgnqpzCiHGAUFo2SIOVrhXk8a1CvO2AODgZMV/eLmpLzPelYI0g/ZCefV1eiwWi9o9Ka4LdDodVqu1Zeds0dnalr/Ynl8VQvQtuyiE8APesr1cIaWs+hNr6riWZ+SDAHiEFuIaUFS2Eq4d8cJa5ui1+fEahyo6D+rMSXE90Bq/551WoKSUnwNvo2V9OCmE2CqE2AScQUuVtBkt+WuLjGsVpr0Gbt0RArrflI1w0DSxJMdA2k82U1+ZV59CoVBcZ3RagQKQUj4B/B+a2W4cMBk4i5bsdZaUssaAk6aOaxWWxAECg6sFv6E59svpsW4UpttMfcqrT6FQXId0aIGSUs63nT39vY4+H0spR0spPaSUrlLKcCnlqvpMdE0d1yrMfBcA734FuPgW2xYouHrIG6tZefUpFIrrkw4tUNcNNq8+IaD7LVnoykx9uQ6knFBefYrrB6vVSkhICEII/Pz8KC0tbe8lAbBmzRqEEMyfP79N7rd8+XKEECxfvrxN7leV+fPnI4RgzZo17XL/MpRAdRRsXn1GNwv+w8tNfZnxbuQnGbUXyqtP0cX59ttvSUxMBCA1NZWtW2statBpuXjxIkIIQkND23spHR4lUB0Jm1efZ+8C3AKL7JevHvbGUmIz9X3xaHusTKFoE/7zn/8A0KNHj0qvrzcWLVpEbGwsixZd37VTlUB1JCp59WWhN2q+GuZCPcnHPLU+0qpMfYouSUZGBlu2bEEIwYYNG9Dr9XzzzTdcvXq1vZfW5phMJsLCwjCZTO29lHZFCVRHw+bV5+BsJeCmCrn6LrqQk1ihLIdyPVd0MT766COKi4uJjIxkzJgxTJo0CYvFwrp162rsL4Swx958+umnRERE4Obmhru7OxMnTmTfvn01jjt8+DDPPPMMI0eOxN/fH6PRSGBgIHfffTeHDh1q8HrXrVuHEIIpU6bU2ufkyZMIIejRowdms5n58+fTq5dWySchIcH+Hqqa/Oo7g4qNjeWRRx6hb9++ODs74+3tzdChQ1myZAkJCZWrCW3cuJEFCxYwePBgvLy8cHJyom/fvixcuNBuTu2oKIHqiNi8+jyCi/DoWWC/fO2IF6X5tqwEMe+r8yhFl+KDDz4AsDsiPPDAA5Wu18YLL7zAfffdh9FoZOrUqQQFBbFz504mTpzIwYPVE8I899xz/OMf/6C0tJSbb76Z6dOn061bNzZu3MiYMWP47LNqxbZrZO7cufj5+fHtt99y9uzZGvusWqXl3XzkkUdwcHBgzJgxzJo1CwBXV1fmzZtnf9x9990Nuu+6desYNmwYq1evRkrJtGnTGDduHFarlZUrV7Jr165K/efMmUNUVBSurq7cfvvt/OIXv6C4uJi33nqLESNGEB8f36D7tgftXm5DUQNDZ8PelZAWR0B4NgWpRswFDlhLdVw55EXP8ekIHap2VBcgdOm29l5Ck7m4YmqLzXX8+HF+/PFH3N3d7R/Ud955J926dSM+Pp59+/YxZsyYGseuWrWKI0eOEB4eDmiegI899hirV6/mhRde4LvvvqvUf8mSJaxfvx5/f/9K17du3cqsWbN47LHHmDp1Ki4uLnWu2Wg08sgjj/Dyyy/zzjvv8Pe/V46GycnJYf369Tg4OPDwww8D8NBDD3H77bezceNGTCZTo73kjh49yoMPPoiUkvfee48FCxZUyuAQGxtbbczHH3/MtGnTKr0fs9nMn//8Z15++WWeeuop/vvf/zZqHW2F2kF1VBYdBgR6o6RHRBYIrexUYaojaT/bko9azeo8StEleP99LRh99uzZ9g9So9HIfffdB9TtLPHnP//ZLk6g5YR7+eWXAYiOjq7mqj5lypRq4gTwq1/9invuuYeMjIxqu5DaePzxx3FwcOCDDz6gqKioUtvatWvJy8tjxowZBAYGNmi++vh//+//YTabWbJkCQ8++GC19EIDBw5k4MCBla5V/JmW4eDgwEsvvURgYCDffvstubmtU12ouagdVEdm5ruw6WFcfEswDc4l7ZQHAGk/uePqX4KLb0m567naSSk6KcXFxXzyySdAuVmvjAceeIA33niDzz77jH/96181ZoafNm1atWt+fn54e3uTmZlJeno6AQEBldrT0tL46quvOHXqFFlZWZjNZgBOnToFQHx8PFOn1r9DDAwMZObMmURFRbFhw4ZKcVJvv/02AAsXLqx3noZgsVj4/vvvAW0n1hji4+P55ptvOHv2LHl5efakrmazGavVytmzZxk+fHiLrLMlUQLVkRk6G45/BBf2YBqUR0GyIwWpjiAFVw560XtKKnqj1FzPlUB1SlrSTNZZ+eKLL8jIyKBfv36MHj26Utvw4cMZNmwYP/74I1FRUSxYsKDa+JCQkBrn9fDwIDMzs9rO5t///jdPP/00BQUFNY4DzTzXUJ588kmioqJ466237AK1a9cuYmNjGTx4MOPGjWvwXHWRlpZGfn4+Dg4O9O3bt/4BaAL0xBNP8N577yFl7cW/G/N+2xJl4uvozNsCpjCEDgJHZaKz1Y4yFzhw7YgXUqK5nqtUSIpOSpn5Ljs7mzFjxlR7JCcnV+pXlcaUGI+JieHxxx+ntLSUv/3tb8TFxdl3FFJKli1bBlDnh3lVRo8ezfDhwzl69CgxMTFAuXPEE0880eB5WoPXX3+d1atX0717dzZs2MClS5coKipCSomUkoiICKBx77ctUQLVGVh0GACDq5XuN2fZL+dedibzrM22rFIhKTohiYmJ7NixA4CUlBT2799f7XHt2jUA9u/f32yPs88//xwpJU8++SRLlixhwIABuLq62s9yavPGq4/FixcD8NZbb3H16lW+/PJL3N3d+fWvf92s9VbEZDLh4uKC2Wzm3LlzDRpT5pH473//mzlz5hAcHIyjo6O9vanvt61QAtVZKKsdFVSEV9/yssrJxz3Ls56rVEiKTsYHH3yA1Wpl4sSJ9m/1NT3uuUerL9rczBIZGRkABAcHV2tLTU2t5vHXUO69915MJhMbNmxgxYoVmM1mfvOb3+Du7l6tr9GopS4rO/dqKHq9nttvvx2A9957r0Fj6nq/3333HampqY1aQ1ujBKqzMO01MGlmPP9h2Th62TyTrILL+70xF6tUSIrOhZSStWu1otn17TTK2tetW4fF0vRqOGFhYfZ58vLy7Ndzc3NZsGABWVlZtQ2tEycnJx566CEKCwt54403gNrNe76+vhiNRpKTk8nMzGzUfZ577jn0ej1///vfa3RRj4uLIy4uzv667P2+/fbblardnjt3jscee6xR924PlEB1JhYdBrfu6BwgaHQGOkP5edTVQ97qPErRqdi1axfnz5/HxcWFmTNn1tl3ypQpmEwmrl271qyYnQceeIDg4GCOHTtG7969mTlzJjNmzCA0NJSYmJganTAayhNPPIFerwXSR0ZGMmjQoBr7GQwGpk6ditlsZvjw4fzf//0fDz30EEuXLq33HjfffDPvvvuu/b307duXOXPmcNddd3HDDTcwcODAStkwli1bhsFg4N///jcDBw5k7ty5TJo0iUGDBhEcHMytt97a5PfbFjRKoIQQXkKIu4QQfxZCvCOE2CCEeNv2+k4hhFdrLVRhY4n27cjobiHwlvJve/nXnEgvi4/KuwZv3tIeq1MoGkxZhog777yzRlNYRQwGA3PnzgWaZ+bz9vYmJiaGRx55BDc3N7Zt20ZMTAwzZ87k2LFjNZrCGkpwcLB9x1Kfa/nq1at58MEHsVgsREVF8f7777Nhw4YG3WfBggUcO3aM+fPnU1payubNm9m7dy96vZ5nnnmGCRMm2PtGRERw5MgRpk6dSnZ2Nl9++SWXL1/mueeeY/v27RgMhia/37ZA1Oe9IYTQAzOBJ4CxQFlkWMUIMVnheS/wFvBFm1ambSdGjhwpyzx32oyvnrZX2U3+0YOMOJswCUnIuHRcA0psi3tQMw0q2o3Y2NhqgZOKrsn//vc/hg0bRmBgIAkJCTg4XH9RPI35fRdC/CClHFlXnzp3UEKIe4HzwAa00ugZwDbgb8BSNNFaCvwd+BrIBCKBT4FzQoi5DVqponFUOI/yG5pTqQrvlYPelBbY/ltVqXiFos144YUXAC0u6noUp9ag1p+iEGI/MApIBV4H1kop/1ffhEKIYcB84F5gvRBisZRydN2jFI1m0WH4exgi7xqBt2ZyYbsvliI9lmI9l/f70HNCGjo98JcQWHapvVerUHRJtmzZwpdffsnJkyc5evQooaGh130Np5akrh1UH+BpIERK+XRDxAlASvmjlPK3QDDwe9s8itbAVprD4GylR0SmPV9fUbqRpBhbEG9xtnKaUChaiWPHjvGf//yHuLg4pkyZwjfffIOrq2t7L6vLUKdASSlfl1KWNGViKWWJlPKfQO+mLU3RIGylOVz9S/AfVp6uJPuCC5lnbH8oKohXoWgVli9fjpSSnJwc/vvf/zJgwID2XlKXolaBklLm19bWGKSUtSe8UjSfobOhl5bry7t/Pp6h5T/u5OMe5CdrQYGqyKFCoehsqDiorsC8LfZS8QE3ZeHkY9v0SsGV/d6U5KkihwqFovPRYIESQvgIIUYIIXyqXO8uhFgjhDguhPhCCHFjyy9TUS9L4sDRE50egsZkoHfSPPwtJXou7/PBalaZJhQKReeiMTuoZcBRNOcHAIQQRmAf8GvgRuBOYJcQokdLLlLRQGzeegYXK0FjMhA6zWmiOMvA1cMVMp//pebyBAqFQtGRaIxAjQcuVPHmmwP0AvYAU4BVgBeg/CzbC1tSWRdTKQHh2fbLuYnOpJ60Resrzz6FQtEJaIxABQFVc7NPQ8se8ZCU8lsp5WLgAnBHC61P0VgqBPF69SnAu195Qsz0n93JOu+svVCefQqFooPTGIHyBtKqXIsA4qWU5ytcO04FM6CiHVh0uDzz+fAcXLuXVxS9FuOlPPsUCkWnoDECVQh0K3shhAhG21Xtr9KvGHBE0b7YMp8LHfS4NbNKeQ4finOUZ59CoejYNEag4oAxFbz47qM8OWxFgoDkFlibornYMp/rDZLgsel2zz5riY7Evd0wF9v++zc90l4rVCgUilppjEB9CLgCR4QQUcCLQB7wZVkHIYQjMAI43ZKLVDSDmasBrVx88G0ZCL1WQ6o0z4HL0d5YLQBSefYpFIoOR2ME6m3gY7TURXcDJcDDUsrsCn1+hSZie1pshYrmMXS23bPP2aeUHhFZlFVHKUxzLC90WJytRErRLoSGhiKEYPfu3e29lBq5ePEiQghCQ0OrtZWt/eLFi826x/LlyxFCsHz58mbN09VosEBJKa1SyvvRkr/eCvSQUlY9vDgP3AOsbbklKprNtNfs6ZDcg4rwq5CzLzfRmeRjHkqkFIoOyO7duxFCEBkZ2d5LaRdqFSghxFwhRLUyl1LKC1LKQ1LKnBrajkkpN0opk1p6oYpmYkuHBOAzIL+S+3nmGTfSY21FD1WMlELRYHbs2EFsbCw9ejQvN8GiRYuIjY1VpTqqUNcO6mMgRQjxlRDiQSGEb1stStFK2NIhCQH+I3JwDy60N6We8KgcI6VKxisU9dKnTx/CwsKaXTrdZDIRFhaGyWRqoZV1DeoSqGeBH9GCbt8FrgohdgkhnhRCKDtQZ2XZJbtIBY7KxMWv2N507agXeVdtEQJpcUqkOisnouAfQ2C5l/bcCcMI5s+fjxCCNWvW8NNPPzFr1ix8fX1xc3NjzJgx7Nq1y973q6++Yty4cXh6euLh4cH06dM5c+ZMtTkrmsvy8/NZunQpvXv3xtHRkeDgYBYvXkx6enqj1lnXGZSUkqioKO644w78/PwwGo306NGDiRMn8uabb1bqW9MZVGRkJOPHjwdgz549CCHsj+vF5FdXuY2/SSkj0NzGn0RzJx8N/BO4IISIEUIsE0Ioe1BnwyZSZYll7TFSUnB5vzeF6bZvg0qkOh8nomDrk5CdCEjteeuTnVKkAGJiYrj55puJj49n4sSJDBgwgP379zN58mSio6N54403uPPOO5FSMnnyZHx8fNi6dSu33XZbrWJTUlJiF4khQ4bwq1/9iqKiIt58800iIiJITm5+lExJSQl33XUXc+bM4bvvvqN///7cfffdhIWFcerUKRYvXlzvHFOmTGHy5MkA+Pv7M2/ePPtjypQpzV5jZ6DWku9lSCmvoeXYWyWE8EZLCDsTuB3NpfxlIUQ8sBHYLKWMacX1KlqKZZdguSd6oyR4XDoJ35sozXdAWnQk7vGh58R0HD3Nmkitna6dYSlanuWerX+P0kLY9LD2aEmWZ9ffp5msWrWKlStX8vTT5RlPnn32Wf7617/y0EMPkZSUxO7duxk7diwARUVFTJo0iejoaN566y2ef/75anMePHiQ/v37c/r0afvZUW5uLjNmzGDHjh0sXryYqKjmCfozzzzDli1b6N+/P19++SVhYeXf4y0WC9u2bat3jqVLlzJq1Ci2b99OWFgYa9asadaaOiONqgclpcyUUq6RUk4HfIG5wGdAIPAH4LAQIkEI8Q8hxDghhGj5JStajLIYKWcrwePS0TuWl+i4tLtbeR0plRJJ0U5ERERUEifQPrgB4uPjWbhwoV2cAJycnPjd734HUMkMWJWVK1dWcmxwd3fnnXfeQa/Xs3HjRhITE5u85pSUFN5++210Oh2bNm2qJE4Aer2e6dNVHsyG0OSChVLKfClllJRyLppY/QpYAzgDTwE7gedaYpGKVqJCjJSjh4Xg2zLQOWiBvOZCPZd2daO0wPYrEvO+EilFm1OTKcvb25tu3brV2t6vXz8Arl69WuOcXl5eTJs2rdr1vn37MmrUKKxWK3v3Vk2Q03B27txJaWkpERERDB48uMnzKBpg4msIUsoSYBuwTQihA8YBM4CUlphf0YpMe017jnkf526lBI3NIHFvN6RFUJrvwKVd3eg5MR0HJ6smUiGjNGFTtAwtbSYrO4MqLffQxOAMv/pXp/x/CwoKqvG6m5sb6enpNba7uWkhE0VFRdXagBoDbiu27d+/n8uXLzd+sTYSEhIAqu2cFI2nxUu+2wJ6d0kpn5RSvtvS8ytagWmv2XdSrv4lBI3OAFuxw5JcA5d2d8NSYrPWqrx9HZuhszUx8gwGhPbcScUJQKer+yOqvvamok4nOgZN2kEJIQLQzp2causjpTzQ1EUp2oFpr0H6WbiwB7fAYnpEZHLlgDdIQXGWgcQ93QiOTEdvkNrB/szVnfZDr8szdLb6v6mDutISlbUFBgY2ef6ePXsCcPq0SknaXBr19UMIcY8QIha4glb+PbqWR9MNuIr2Y94Wex0pj+Aiut+cZW8qTDdyOdoHq9l2YdPD6kxK0SnJysri66+/rnb9/PnzHDp0CCEEt912W5PnnzBhAgaDgQMHDhAbG9ucpWI0arXbzGZzPT27Jg0WKCHEXGADMADIAU4AB2p5HGzxlSrahgrFDr16FeIfXi5SBSmOXN5XQaRULSlFJ+X3v/89165ds7/Oy8vj8ccfx2KxMGPGDEJCmp6LwM/Pj8ceewyr1cqsWbOIj4+v1G6xWNi6dWuD5irzNDx79ux1KVKNMfH9wfb8FPC2lPL6+2ldLyw6rOXjy7uGT78CrKU6Uk94AJCf5MTlfT4EjclA54B2JqXMSYpOREREBBaLhf79+zNhwgSMRiN79uwhNTWVPn36sGrVqmbf429/+xvnzp3j66+/ZvDgwURERBAUFERKSgonT54kJSUFKWW98/Ts2ZPhw4dz/Phxhg4dSnh4OI6OjgwYMIBnnnmm2evs6DTGxNcPOCClfEOJ03WALW8fgGlQHqYh5bmBy0RK20nZzqTUTkrRSTAajezcuZNHH32UEydOsGXLFoxGIwsXLuTQoUMEBAQ0+x6Ojo5s3bqVDz/8kNtuu41Tp07x+eefExcXx9ChQxslgps2bWL27NlkZGTwySef8P777zco0LcrIBqi4gBCiCvAHinlfa27pM7FyJEjZUxMF06e8ZcQLcM5kHrKjbRTHvYm14Ci8p0UKMeJGoiNjWXgwIHtvQwFWi6+8ePHM27cuA5be6qz05jfdyHED1LKkXX1acwO6lvgpkb0V3QFll2yl+nwHVLXTgrlgq5QKFqUxgjUnwAPIcSrQgh9ay1I0QFZEtdAkVLmPoVC0XI02ElCSnlJCDEG+BKYKYTYAVwGrLX0f6VllqjoECyJs5v7fIdoxQ7LzH35SU4k7u1G8NgMdAapuaBfOlSepUKhUCiaQIMFypb4dRGas4QerfR7TQdYwnZdCVRXY9mlWkWqIMWRS7u7aUlnjVJzQQclUooOQ2RkZIM85xQdh8a4mS8FFgNmtLx7Z4G8Okcouh5VRErosLugF6YbSdhpIiSyQu4+UCKlUCiaRGME6kGgABgjpfyxldaj6AxUECnToDx0DlaSj3kBUJxlIGFnN0Ii0zG42EQq/ayqJ6VQKBpNY5wkegB7lTgpAHtVXgCf/gV0vzkThC3BbI6BhB2myvWk1qr6NwqFonE0RqCuoO2gFAqNCi7oXr0L6RFRLlKl+Q4k7DBRnKNESqFQNI3GCNSnwDghhGtrLUbRCanggu4RogXuClupDnOhnoQdJgozDFrfC3vgzVvaa6UKhaKT0RiBegmIB7YIIfq00noUnZEKIuXeo5jg29IRei36wFKs59LObuQlOWp90+K0PH8KhUJRD40RqC1oHnzjgVghRJwQ4nshxLc1PLa3znIVHZYlcfYs6K4BJYSMT0dn1ETKataRuNeH7ARb+bC8a5qThUKhUNRBY7z4bq8yrr/tURMq2OB6ZNFhzYSXFoeLqZTQiWlc2t0Nc6EerIKrB32wFGfj0z9fy+/3lxDtHEuhUChqoDEC9YtWW4Wi61BBpBw9zYTensqlPd0oydHOoZKPeWIu1OE7NBdRnA0vBcDzSe28aIVC0RFpTKqjHa25EEUXYtFhzWPvwh4MrlZCJ6aRGN2NwjStOmh6rDvmIh3db8pGUKhKyCsUihppVMl3haLBzNsCIx8EQO8oCYlMxy2wyN6cfcGVxGgfLKVCu6BKyF+3hIaGIoTo0CUwLl68iBCC0NDQam1l67948WKz7rF8+XKEECxfvrxZ81Tl+PHjvPLKK0ycOJHQ0FAcHR3x8fFh/PjxfPDBB1itNaZT7RAogVK0HtNes4uUzkESNCYDz17loXT515xI2GGitMD2axjzvoqVUiiqsHv3boQQREZGNnqs2WxmxIgRPPfccxw5coTevXszc+ZMBg8eTHR0NAsWLOCOO+6gqKio/snagVoFSgixVwhxa3MmF0KMFkLsbc4cik7OtNc08x0CoYPuN2fRbVCuvbk4y8DF73wpyrRZm1WslKKTsWPHDmJjY+nRo0ez5lm0aBGxsbEsWrSohVamER4eTlRUFGlpaezcuZNPPvmE6Ohojh8/Tvfu3fn222/5y1/+0qL3bCnq2kGFAdFCiO+EEHOEEL71WcoAACAASURBVI4NmVAI4SiEuFcI8T2wl9o9/RTXC0Nnw/Is0DsjBPgNza2UGqksoDfvqoqVUnQ++vTpQ1hYGAaDoVnzmEwmwsLCMJlMLbQycHBwICYmhnvuuQdHx8of4TfccAN//etfAfjoo49a7J4tSV0C1Q94AxgHfAwkCyG2CSH+KISYJYSIFEKMsD3PEkI8L4T4GkgBPgLGAq+jBEpRxvNJoHcGtNRIIePS0RkqxEpF+5B51kXrm3dN8/BTNJpt57cx6fNJDF07lEmfT2Lb+W3tvaQmMX/+fIQQrFmzhp9++olZs2bh6+uLm5sbY8aMYdeuXfa+X331FePGjcPT0xMPDw+mT5/OmTNnqs1Z0VyWn5/P0qVL6d27N46OjgQHB7N48WLS09Mbtc66zqCklERFRXHHHXfg5+eH0WikR48eTJw4kTfffLNS35rOoCIjIxk/fjwAe/bsQQhhfzTF5FeV4cOHA3D58uVmz9Ua1OrFJ6XMBn4rhPgXWpmNecAdwJQ65hNABrASeEtKebHllqroEjyfpO2O8q7hGlBC6O1pJO7xobTAAaQgKcaLklwH/IblICzKw6+xbDu/jeUHllNk0c4UruVfY/mB5QBM7T21HVfWdGJiYli4cCG9e/dm4sSJnDlzhv379zN58mR27NjBjz/+yG9/+1tGjx7N5MmTOXLkCFu3buXo0aOcOnWKbt26VZuzpKSEiRMncurUKSZMmMCIESPYs2cPb775Jtu3byc6Ohp/f/9mrbukpIR77rmHLVu2oNfrGTVqFCEhISQnJ3Pq1Cl27txZrzlvypQpODk5sX37dvz9/ZkypfzjNyys+VaGMhHv3r17s+dqDep1M5dSngd+J4T4A9puKhIYBvgBnkAW2q7pGLALiJZSFrfWghVdgCVxlWOlfpFGYrQPRRmaG3rGaTdK8vQEjspCX1ah9/hHXbJkxw1rb2j1exRZilgavZSl0UtbdN6T80626Hy1sWrVKlauXMnTT5d7eT777LP89a9/5aGHHiIpKYndu3czduxYAIqKipg0aRLR0dG89dZbPP/889XmPHjwIP379+f06dP2s6Pc3FxmzJjBjh07WLx4MVFRUc1a9zPPPMOWLVvo378/X375ZSVBsVgsbNtW/8526dKljBo1iu3btxMWFsaaNWuataaKSCntJr5Zs2a12LwtSYO9+KSUhVLKb6SUS6WUU6SUI6SUfaSU4VLKO6SUz0kpv1fipGgQiw7bUyM5OFvpOSEdtx6F9ua8K84kfF+lZIc6l7ouiYiIqCROoH1wA8THx7Nw4UK7OAE4OTnxu9/9DqCSGbAqK1eurOTY4O7uzjvvvINer2fjxo0kJiY2ec0pKSm8/fbb6HQ6Nm3aVG23o9frmT69fT1W//znP3Pw4EH8/f1ZtmxZu66lNpSbuaL9WHQYeo0DbG7oozPxCSsv0lycbeDitybyk7Wdlcrhd31S0axVhre3t910V1N7v379ALh69WqNc3p5eTFt2rRq1/v27cuoUaOwWq3s3dt0B+SdO3dSWlpKREQEgwcPbvI8rcW6det48cUXMRqNfPLJJy3qmNGSNCbVkULR8szbAieiYNMjCJ3Ef1gOjp6lJB31QloFlhI9l3Z3I2BENt79CrQcfl3oXKqlzWRVz6AAnPROLL91eac9gwoKCqrxupubG+np6TW2u7m5AdQa31NTwG3Ftv379zfLcSAhIQFomXOiluazzz5jwYIF6PV6NmzYYHfC6IgogWoHtp3fxgv7X6DEWmK/NipgFKsnr27HVbUjQ2drj5cCwFKIV69CHN3NJO7zwVKk15wnfvCiKNtAwIhshA7tXOrSIS3OSmGnTIReP/Y6SflJBLgG8NSIpzqtOAHodHUbeuprbypCiFaZtz3ZtGkT9913H1JKPvroI2bMmNHeS6oTJVBtzLbz2/hD9B+44eo4LMLMuW4/UmjM5VDSIfuB+ZwBc/jjqD+280rbgQoefs6mUnpNSuVytA9FmZqJL+usKyXZDvQYk4mDo1XLPJF+tks6TzSHqb2ndmpBagvqSktU1hYYGNjk+Xv27AnA6dOnmzxHS7N582bmzp2L1Wrlgw8+YO7cue29pHpRZ1BtzOvHXgerYPjlXzDm4t38+ocXmfrzY/RPuQmDWQuk+/T0p9yw9gZuWHsDLx96uZ1X3MYsibOfSxlcrPScmIZHSHl6pIJURy5ur1KlVzlPKBpJVlYWX3/9dbXr58+f59ChQwghuO2225o8/4QJEzAYDBw4cIDY2NjmLBWjUfuCZjabmzzH1q1bmT17Nmazmffee4/f/OY3zVpTW6EEqo1Jyk8iOCsMJ4srADp0BGcPZMK5+5kX8zIT439DcOZAhNT+a8rEasS6EZ024LLRVEg0q3OAwIgsfIfmUFZmrLTAgYTvTWSd14J+VVCvoin8/ve/59q1a/bXeXl5PP7441gsFmbMmEFISNMdcvz8/HjsscewWq3MmjWL+Pj4Su0Wi4WtW7c2aK4yT8OzZ882SaS+/vpr7r77bsxmM++++y4PPPBAo+doL5SJr40JcA0gqfQCe3p/Sr/UcAJz+9rbHKSRfunh9EsPp8CQw1nTMU77HiHd5QqllNpjWa6L86ppr0HIKM15QkhMg/Jw9Czl6iFvrKU6pFVw7Yg3hWlG/MOz0ZWV7Rj5oDqXUtRLREQEFouF/v37M2HCBIxGI3v27CE1NZU+ffqwatWqZt/jb3/7G+fOnePrr79m8ODBREREEBQUREpKCidPniQlJQUp66/t2rNnT4YPH87x48cZOnQo4eHhODo6MmDAAJ555pk6x6akpDBz5kxKSkoICgpi37597Nu3r8a+LRlj1VIogWpjnhrxFH+I/gOx/geI9T+Aa7EXfdNG0C8tHFNBuTeSS6kHQ69FMvRaJBnO14j3PcIZ0w/kO2ZfP+dVZc4TfwmB4mzcexRr51L7fCjO1kx8WeddKcoyEDQ6E4OrRTuXivtaMxUqFLVgNBrZtm0bf/rTn9i4cSNXr17F19eXhQsXsnz58hZxu3Z0dGTr1q18/PHHfPDBBxw/fpxDhw7h5+fH0KFDG+WgsGnTJp599ln27NnDJ598gsViYdy4cfUKVEFBAcXFWmjq5cuXWbt2ba19O6JAiYYoOIAQ4iFgvZSysN7O1xEjR46UMTExjRqz7fw2lkUvQ1L5Z++TH0j/1JH0SxuJa6lntXESK1c8zxDve5TzPv/DrC/3AtSh45Wxr3Tdw3Gb8wSA1Sy4dtSTnAQXe7PeaKHHrZm4Bth+JnrnDlGpNzY2loEDB7b3MhQ2du/ezfjx4xk3blyHrj/VWWnM77sQ4gcp5ci6+jTmDOpd4LIQYqUQol8jximqMLX3VE7MO8GogFGVrme4XuVQ6BY+Cv8TXw18i3jTUUp15Yk5BDqCsgcw4ax2XjX+7H10z+kDEqxYWRq9lBvW3sDD2x9u67fU+lRwntA5SAJHZeE/ItueEd1SoufSnm6k/eyGlEBZHr8TzUtXo1Ao2o/G7KC2oCWL1QNW4HtgFfCVbOgkXZCm7KCqUlNcVBkOFiO9M26kf+pN9Mjuh6jhO0WWUwpxfoeJ9z1CgTGnUluXMwHagnrLHCYKUo1c2e+NuUhv7+IWWET3W7I0V3TQUiotOtwOi1U7qI6G2kG1Lu22g5JSTgd6AyuANGASsBm4IIRYKoTwbehcNSzUIISYaNudHRJCXBNClAghrgghPhdCRNYz/j4hRLQQIlsIkSeEiBFCLBRC1Pn+mjqupZnaeyo//PoHTs47yZwBcyq1mfUlxPse5atBb/HRiOUcCvmSDOdrlfp4Ffkx6tKvuP+HPzMl7mFC04eis1b2Auwyu6oKtaUAXHxL6DU5FWff8p1m3lUnLnzjS0GqLUVSWpxKkaRQdEIavIOqNEgIAzAbeAKIQPs6WwJ8jlZm42Aj57sd+M72Mgn4AcgHBgFDbNdfklK+UMPYVbZ1FAE7gFJgIuAOfAHcI6W0tNS4qrTEDqo2Ht7+MIeSDlVvkOCbH0JYyi30TQvH0eJcrUuBIZd401Hi/A6R5ZJcqa3LeAHaMqIDSCuknPAgI86tvF1IfG/IpdvAPOxJAdrYy0/toBTXEy29g2qSQFW5yY3AQuBeoOzU+kc089/6hmQ3F0JMQBOL16WU0VXa5gDr0UyLE6SUuyq0zUITxSTgNinlGdt1f7TSHwOB30opX68yZ5PG1URrClRFahMrB4uBXhk3EpYyih45NR8NJrldIM7vEOdMxynVl/93GISBl8a81LkdK6qY/HKvOHLtsDeWkvJNsGtAEYGjsnBwspn83Lq3mZefEijF9USHEyjbjXoAzwIVq29JtDpRf5RSvt/M+d8DHgT+I6V8sML1GCAcmCelXFdlzDhgN5oI9ZBSWps7ribaSqDKqOu8yqOoGwNSbmFA6i24lXhVay/VFXOu23F+9j9AiluCVl7SRqc/q7Ll8QMozddx5aA3hWnlJa4dnCwE3pqJq1+Fn1sbJJxVAqW4nuhQAmUzzT0BTEPb4RSjlYf/Drgf+KWt6++klP9qxn0WAm8C30opJ9uuBQGJaKZFr5rc34UQl4EewGgp5YHmjKuNthaoirx86GU+Pf1ptetCCoKywghLuYXQzBvQy+rhbmkul/nZ/wBnTDGUOpTvqjq1u3oFV3RphdST7qTHupe3C4lpcC6mQXnYTxl7jWvVXH5KoBTXE+3pZl42qacQ4rdCiDhgO3AXcBX4AxAkpXxISvmplPJXwK1oZ0lPNvY+VSizXVX0Dhhue/6pjtiso1X6Nmdch+OPo/7IyXknOTnvJH08+tivSyFJ9I7luwFr+DD8BfaHbiLduXJdHFNBELddmM1vfniJcefm4JsXDHRyd/UKruhCB3435hI8Lh29o+0oUQrSTnlwaVc3SvNVIUSFoqPTYIESQoywmdquACuB/sAeYBbQW0r5qpQyo+IYKeVhYBvQZBcqIUQAMN/2cmOFpl6254Q6hl+q0rc54zo0m2ds5uS8k6wYuwIHUb5jKjLkc7L7Hj678VU2DXmNON9DlOrKzVwGqyMDU25l1sklzDqxhIHJERgsmmmsLGNFp8oDOG+LZrqz2S/duhfTa3IqLhW8/ApSHTn/jS85l5y0C3nXVMyUQtEBacwOKgZYYPv3e8BQKeUEKeUX9ZzT5NPElEpCCAfgI8AT2CGlrJhdscxdK7+OKcrKs1aw8zR5XKdgau+pHP/NcU7OO1k5EFhAinsCu/t+wofhLxAd+nm1XZVvfjDjzs/l1zEvMvb8bEx5WuqlUllq31V1iuzqZa7otpLyBhcrIePTMQ3OtQf2Wkt1XDngw9VDXlhKbYdxmx7WPAMVCkWHoDECdRF4Bs2M96iU8lQDxz0MGBq7MBvvoLl+J6KdaVWk7Ii/sYdoTR1XPoEQj9hipmJSU1ObOk2rs3ryavuuyqgz2q+XOBTyU/doPrvxVb4Y8g9O+x7BLMp3VUarE4OTR3P3yWeYeeL3hCWPwsGije9U2dUXHbZnRRc68L0hl54T0zC4lmeEzr7oosVMpdl+RdPi1G5KoeggNEag+kgpV0opsxpzA6lRbzxRVYQQr6N57iUBE6WUVROr5dqe3aidsrbcCteaOs6OlPJdKeVIKeVIX98mxye3GRUDgavuqpLdL7Kr73o+DP8T+0M3keFc+cfslx9C5Pl7+c0PLzH2/D345HcHOtGuatprmsnP5hXhYiql15RUPEPLa0yV5juQsMNE6ik37LaATQ/D2untsGCFQlFGYwRquxDi6fo6CSF+J4T4thlrQgixEs2xIhVNnM7U0O2i7blnHVMFV+nbnHFdgtp2VcWGAk5230PUjX9h8+DXiTfFYBal9najxYnByWOYfWIp008tpnf6sM6TrWLobPhTpt3kpzdoufwCIzLRGWyKZHOgSNhpoiRPOVAoFB2BxuTiswJrpJQL6um3GlggpdTX1a+O8X9FMyWmo4nT/2rpF4zmzFCXu3giEASMkVLub8642mhPN/OWojZ3dcdSFwak3syg5FvxKvKv1p5vyOJn/wP87H+AQmPlzWaHjav66mmtJIeN0nw9Vw55UZhaHjMlHKz4D8vBq09BszNQKDdzxfVEu7uZNwAjWjLZRiOEWIEmTpnAL2oTJwApZSJwzHa/e2qYaxyayCQBB5s7ritTm7t6saGAE4G72TDsFbYMeoNz3Y5jpdxa61rqxU2Xf8n9x5YzMf43+OeG2k/2OuyuatprsDxbyyYBGFwt9Byfju8NOXYHCmnWkRTjReJuH0rzbX8iMe+rfH6tRGhoKEKIeh/tldy17P5ViYyMbJF1rVmzBiEE8+fPb9Y8XZEWLVgotP/FcLRkso0d+xJaNoosNHE63oBhfwE+A14VQhyQUp61zeUHvGXrs6IGL8OmjuvybJ6xGaiyqxJw1fMsVz3P4lrsycCUWxmYHGGvWaWXDvZKwGkulzkVEM1Z0w+Y9aV2V/UOl/9vSZx2xnRhD0IHpsF5uHYv5uohL0pyNIeJ/GQnzn/jh//wbDx7FSKKs1XV3lZk8uTJBAQE1NpeV1tX5OLFi/Tq1YuePXty8eLF9l5Ou1Cnia/KWdLtaAG5P9fS3QEtoDYQ+FxKOaeWfjXdZzrwpe1lDPBTLV3jpJQrqox9C3gcLenr95QnffVAy7Z+dy3JYps0ripdwcRXF7WlVtJZ9fTKGMqQpNvontu72rgifT5xfof5KWAfuU7pldo6lPnvRBR88Shl3hFWC6Se9CAjzpWKuaDcAovoflMWDs627yyOnrDsUg0TVkaZ+OonNDSUhIQEdu3aRWRkZHsvpxplu6eqn5WXLl2ioKCAkJAQXFxcahraILKzs7l27Rqenp50797dfr0zClSbpjqynTuVIamUva1WTgB3SinrCoStep/5wAcN6LpHShlZw/j70BLW3oCWcikO+A/wdl27oKaOq0hXF6iK1Jawtlt+D4YkjaFv2kgMVmOlNomVS16xnArYS6LXabsZDTpYVvUKmdFBqzN19bAXpXnlRga90UrAyCw8QorKx9Wzm1ICVT+dVaBaGyVQ9Z9B/cL2mIQmTtsrXKv6GAf0lVIOa4w4AUgp10gpRQMekbWM/1hKOVpK6SGldJVShkspV9UnMk0dd71S0QOwYraKdNcr7OnzKR+N+BMHem4m27HcwivQ0TNrMFPjHufeH59j6NVIjGatPEiHylRRIWYKtDpTvaek4t0vz37NUqIF917e7425qMLZ1HIvFTfVRkgpueOOOxBC8Mgjj1Rrt1qtTJw4ESEEixaV566+ePEiQghCQ0Mxm82sWLGCgQMH4uTkhL+/P/PmzePSpfp3xBWp7wxq+/btzJw5k8DAQIxGIwEBAYwePZpXX32VwsJy36yazqDmz59Pr15aIpuEhIRKZ3GhoaGNWmdnps4zKCnljrJ/CyH2o+1gdtQxRHEdMLX3VHsy2Yq7Ks2pYhcnuu8mJGsgQ5LGEpI1yD7Os8iXWxNmcFPiLzljiuFUwD4yXK/aY6qWRi9t313VtNe0h203pXOQBITn4B5UxNXDXpgLtD+X3ERnCpId8R+ejUdoIUJILW7q+EetmnhWoe1mPvzwQ4YNG8bq1asZP3489957r739xRdfZOfOnQwfPpyVK1fWOMecOXP46quviIyM5MYbb2T//v2sW7eOb775hr179zJgwIBmrVFKyRNPPME777wDwMiRIxk3bhwZGRnExsaydOlS5syZU6fQjBkzhry8PDZu3Iirqyt33323vc1kMjVrfZ2JBjtJSCnHtuZCFJ2TMjGp7FQhueT9M5e8f8aj0MTg5DGEpdyCo0Wz0xusjgxKGc2glNFcdT/LTwHRXPA5gVVn7RhOFYsOV6oz5epfQu87Ukk+7kH2eVdA201dPexNdoIz3W/KxuBq0eKmlnvBzHcbXMYjNqzzmv8GxsW2y31NJhMbNmwgMjKSRx99lJEjR9KvXz927drFSy+9hLu7O1FRUTg6OlYbm5CQQGFhIcePH2fQIO3LU0lJCQ8++CAfffQRv/71rzly5Eiz1vfPf/6Td955B39/fzZv3syoUeXB8VJKdu/ejbe3d51zPPTQQ9x+++1s3LgRk8nEmjVrmrWmzkqbljZXdF1qc1XPcU7jYOhmPgz/E3t6byDN5UqlcYG5ffnFmQf4v2PLCU+cjEuJB9ABzH9V8vnpDZLAm7MJHpdeKVVSfpIT5/7rS0a8K9oRhVRZKJrI+PHja3Ux9/KqXN9szJgxvPjii+Tm5jJ79mwuXbrEfffdh9VqZfXq1fTt27fW+zz//PN2cQIwGo28+eabeHp6cvToUfbvrzP0sU7MZjOvvPIKoJnuKooTaDvA8ePH4+np2eR7XE/UuoMSQvzB9s+3pZSZFV43CCnlK81amaLTUuaqXtH8Z9aXEOt/kFi/gwTk9mZI0lh6ZdyI3hbP7VrqyU2Xf8mIK5M57/MjPwVEk+R+gVLa2fxXZTfl1r2Y3lNSSTnpTma85uknzTqSj3mSc8mZ7jdl4ehpLt9Nzfi+bdfbianLzbwmL7lly5axd+9etm/fztChQ8nOzubRRx9lzpy6HYjvv79qWk/w9PRk2rRprF+/nt27dzN69OgmvYeYmBjS0tIICgpiypQpTZpDUU5dJr6X0Tz3PkcLnC17XR/C1k8J1HVOzeY/SPI4T5LHeVxKPBiUfCsDk2+tEFOlrzemqs1L1Q+drT3KzqYMkoAROXiEFHLtSHncVGGakQvbfTENzqVbWB5CL6EgHVJiwa9mU157mck6IkuXLm2UF1/ZeVTv3r3Jzs5m0KBB/POf/6xzjJeXV7XdWBllZ0KXL19u8BqqkpCg+Yc19xxLoVGXQL2CJjRpVV4rFI3ij6P+aI97uuuLuziXcw6AAmMOMcHfcKzHd7aYqrF0zy03D5oKgog8fy8RCXcS63eInwP2keOUbneq+EP0H9q2+m+V3ZSLqZRek1NJ/9mdtJ/dQAqkVZB60oOcS84EjMzWxpmL4OpxcDGBV3Cdt1A0js2bN5OXp3laXr58mStXrtCnT596RtVNTVkjFO1DrQIlpfxjXa8ViqZQU6YKq87COdNxzpmO1xhT5WhxYdi1Cdx4LdIWUxVNolccVmFte/Nf1d2UXivj4R6s7aaKMrQ1F2cbSNhhwjJbh7TakqkXpGkPr57g4tP6a+3inDp1iqeeegqj0cg999zD+vXrmTNnDgcOHMBoNNY4Jisri+zs7BrPgMpijQIDA5u8pp49tRzUp0+fbvIcinKUk4SiXShzqqiaVb0spurD8Bc40POLWmKqHuPe4+0cU7XocKXKvU5eZkJvT8NvWDZCXx5GZzXrKM52wFyswx7nmZUAaTUl6Fc0lPz8fGbPnk1hYSGvvvoq69atY/z48fzwww8888wzdY5dv359tWvZ2dl89dVXAM0KFg4PD8dkMnH58mW2b9/e5HkAu8iazeZ6enZdlEAp2pWKtaoqev+VOBRyInA3nwx/mW1h75Dg9ROyQg5iz2ItpurXP7zIbefm0C1f+9bbpnWqqnj6CR10C8unzy9TcetRHogppaA0X09JrgNWi818VJKnmf2yElt3jV2UhQsXEhsby/Tp0/ntb3+LTqdj/fr1+Pn58a9//YvNmzfXOvbFF18kNrb87K+0tJSnnnqK7OxswsPDGTNmTJPXZTAYWLZsGQAPPPBANZf1Mjfz7Ozseufy9fXFaDSSnJxMZmZmk9fUmWlwHJQQ4nHgdWCGlLLGr6hCiGnAJuAJKeV7LbNExfVCzYlqJYnesSR6x9YSU2VkUMqtDEq51RZTtY8LPv/DqrPy6elP+fT0p61v/is7m7Ll9DO4Wggem0nu5UIuVEjtZDULirMdcHCy4OBs1Up5FKRBUTYEDGm99XUSVqxYUWe8z3333cekSZNYt24da9euJTg4mA8+KM+Q1r17dz788EOmTJnCggULGD58uN3kVkZISAjh4eEMGzaMCRMm4OnpycGDB7l06RImk4l169Y1+3387ne/IzY2lvfee49Ro0YxcuRI+vbtS0ZGBj///DOJiYlcuHChXldzg8HA1KlT+eKLLxg+fDijR4/G2dkZk8nEihUr6hzbVWhMNvOZQAbw3zr6/NfW525ACZSiSZQ5VVRNVFsWU3U0eBt908IZkjQWU0GQfVxgbl8Cc/uSb8jmZ/8DxPofoMCY0zbef2VnU7YM6QDuQUU4OFtwcLKWp0YCzEV6LCU6DC4W9EYJ1lJtN2V0A1O/ll9bJ6E+k9iwYcMICQnhiSeewMHBgU8++QQfn8pneZMmTeLZZ59lxYoVzJ07l71792IwGOztQgiioqJYsWIFH374IQkJCXh4eHD//ffz0ksvtUgaISEEq1ev5s477+Sdd97hyJEj/Pjjj/j4+NCvXz8WL17c4Mzsq1evxsfHh+3btxMVFYXZbKZnz57XjUA1pmBhIlo28V/U0+87YICU8roonnM9JYttT2pMVCshILcXg5PG0jtjmD2mqgyLsHDB53+cCogmyf18pVTHrZpR/UQUbH4crGZiJ0cxsKcfVrOgtECP1VzZQ0xvsOLgYkFXcenKiaLF6YyJVzsjbZrNvMpkRWhlNKpHuVXutx6YJaV0atDEnRwlUG1LbdV/XUo8GJh8K4MqxFRVpGpMVRl9PPrYTYstzokoYvM9GdjTDwApwVKsw1yop+qfXSWzXxlKqFoMJVBtQ3sKVDLaDmpcPf12A0OklNdFRkMlUO1D7XWqdPTKGMrgpLEE5lZPd1OsL9DqVPnvI8e53EOwtcx/sbGxDOyGFgtlQ1qhtEAz81VE6CQGFws6gywXquvc7NdSKIFqG9pToL4BxgODpJTnaunTB4gF9kopb2/QxJ0cJVDtT227qm75gQxOGku/tHAM1sqJQyVWEr3iOBUQzSWv2Ep1qlrS/Gf/gy3I0NzLK2A1a959ds8+GzoHicHVgk5fmx8DQQAAIABJREFU4W9TBfk2CyVQbUN7CtRc4GO0arczpZRnqrT3Bb4ABgHzpJQfNWjiTo4SqI7DtvPb+OO+P2KWleNGjGZnwlJuYXDSGDyLfauNy3ZM46eAfZz2PUyxocB+vSXMf9X+YLMSNc89G+VmPx1SVhYqB0crDs4WLci3DGX2U3Rg2lOgBLAV+CVgBvahVaAFGACMRfMK/EZK+csGTdoFUALV8ajN/IcUBGcNYEjSbYRkDURUCQM0ixLO+P7AqYBo0l0rZ11v6q6q1j/YtDNaLFTZ0qxgLtRjLq5i9hMSB2crescK51MOTrXm9lMo2pN2EyjbhEbgH8DDVHdRNwOrgaellMUNnrSTowSqY1Ob+c+jqBuDkrSYKieLa7X2JLcL/BQQzbluP2LVWezXGxtTVe8f7LUTIMvnr83bT6eXOLhY0BuU2U/RcWlXgaowcQAwESiLgksAdkgpkxo9WSdHCVTnoLZdlYPFQN+0cAYnj8E3v/qHfYEhl1i/A/zsf4B8xyz79Yaa/xr0B1uD2c9aqgmVtNbklm5V51OKDkmHEChFOUqgOh+1xVT554UyJGksvdOHoZeVDQRWLFz0OcWpgGiuepxpcExVY/5gSYmt7O0nwVykuaVXRZ1PKToiSqA6GEqgOi+1mf+cS9wZmBLBoOTRuJVUrx2U4ZzETwHRxPsepVRfbs2uyfzXKIECm7ffJSpWtrFawVyTW7rQ4qf0Tip+StExaHeBEkIMAJ4EIoEetstXgF3Am1LKuFqGdkmUQHV+ajP/CakjNGMIQ5LG0iOnf7VxJboi4n2P8FPAPjJdku3XdejsdaoaLVBl1OaWXsP5VI3xU8qRQtEOtLeTxHzgbcBIJSOHnRLgUSnl2gZP2slRAtW1qNH8B3gX+DM4aQz9U2/GaK2eJOWKRzynAqK56HMKKcqzrr99w9uMGdH07NiNOZ/SOUgcnKs4UqhAX0Ub0p5u5jcBB9BKdHwB/Ac4hyZUvYAFaAllLcBoKeXRBk3cyVEC1TWpLabKYHakf9pNDEkai3dh9YSfecZMLVGt30EKjbn8c9A/CegVgEAQ6B6Il2PN5cbrJemUllTWRl1pk5QjhaK9aE+BigJmAfdLKT+ppc+9wHrgMynlnAZN3MlRAtX1qc2pIjCnH0OSxhKacQO6KjFVFmHmfLcfuWv6Lfj2rpz1y9XgSqhnaOMXUsP5VG3xUwB6o+ZIUSkRrRIqRSvSngJ1Fbgspby5nn6HgRApZfcGTdzJUQJ1/VCb+c+12ItBtkS1zmb3Sm033e9NUM9gCg35FDvkIyuIS5N3VTWdT1k0oarqSAG2RLROVuXxp2h1WlqgGlNRtxsQ34B+ZwD1m6/ocqyevNpept5BlLuh5zv+/+2deZxcVZn3v0+tvabXENLZOukskCgkJsxEQQmgAQzKMkAcBgU3+ExQwfFV4RU3YJjg6AiOSzSCIYAKimQcHMm8IhFZIoQlJCEJ2TpLhwC9d1fXXuf9497qVFdXVXelqruru5/v53M/N3XPueee+6S7fn3Oec7ztPPi9P/hgcXf5MnZ6zlWdqDPfa6Yh/JgFTW+KZQFq3DGrPxEBkNTVxM7mndwtPvo4DtSUg11i6zRkI3DCZ6yKN4JERzuWJ/qkYCTYIfbDqdkX2w/WDAZfevr6xGRAY9NmzaNWB/jfUhm2bJleenbunXrEBGuvfbanNpJZvfu3Xz/+9/nwgsvZPbs2RQVFVFRUcF73/te7r77bkKh0MCNjCDZJCxsAxoGrAWz7LqKMiZZMWtFb9TzxFFVzBFlz8SX2DPxJWq7p7Lg2PtZwgW99wlCcaSM4kgZYUcQv7uboMsPGNoCbbQF2vA6vcyu6h+FPSWV06wjIWySw2XwlkeJhmNEeo4HojUGwvZUoKsoIXRST7N1FMCI6vzzz8+YyG+wSf7GErkGuT3vvPNoamqiqKiIJUuWcMYZZ/DWW2/x/PPPs3nzZtavX8+f/vSnfokfC4VsBOo54GIRudgY81+pKojIR4ClWE4UijLmie97St5T1Vx2hL/M/hWXeZfS7WmnKFKKK3Y8s6s75sUd9BILxQi4fATc3UQlQjAaZEfzDgCqiqqoK6sbuBNxL70EoXK6DY4JkX4efyZmfY4EHFaMP48tVO0HrWMEhermm29m2bJlI/LsE2X9+vX09PQwfXpu+VkvvfRSli5dOmAa+GyZN28et912G1deeSVlZWW91xsbG7nooot45ZVX+OIXv8j99xem43U2a1BnAX/B8tJ7ELgfOIC1YjsL+ARwNeAEzjbGPDsUHS40dA1KSeaSxy5hX6eVkSbuxQfgjnopjpThiZSk3qPhDOB3dxNy+vtcz9qpIq3HX/+I6eI0uIuT9lDBsApVfX09Bw8e5KmnnipYgYpP7w13YIOhTBPyzDPP8P73v5+ioiI6OjrweDw5tzlia1DGmGeAm7Dcyq8B/owlUI32vz9pt3fTeBEnRUnFhks3sO2abayc19eRNewM0ultobWkCZ+ng6gkpQWJFlERqKW6p46S0AQcdgp7X9jHjuYd7GjeQXuwnQE5+V2WwNhYESdieCsjVnikBCEyUSHU7SLU6SIakr5rVMe2n9D7DzXGGC688EJEhOuuu65feSwW47zzzkNE+NznPtd7vbGxERGhvr6eSCTC6tWrOfXUUykqKmLSpElcc801HDp0KKu+DLQGtXHjRi677DLq6urweDycfPLJnHnmmdx11134/cf/EEm1BnXttdcyc+ZMAA4ePNhnPa6+vj6rfqZi0aJFAAQCAVpaWnJubyjIZooPY8x/ikhcqD4A1GEJVhPW6OoeY8wree+looxCbl16Kzt37mRy+WSOdh/t/es7JjF63J30uDvxRIsoDpfhiRb33uc0TkrDFZSEKwg5/fjdXb0hlZq6mmjqahp4VFVSbR0JG31FwF0cw+WNEQk4iAaP76GK2UIVj5rucBkkFrYcKRConD7ia1RxRIQHHniAhQsXsnbtWs455xz+8R//sbf8tttu489//jOLFi3ie9/7Xso2Vq5cyeOPP86yZcs4/fTTefbZZ1m/fj1PPPEETz/9NPPmzcupj8YYVq1axZo1awBYsmQJZ599Nq2trezcuZObb76ZlStXZhSas846i+7ubh599FFKS0u5/PLLe8tqa3NPWL5nj5XSz+PxjIk1KABsAbpmCPqiKGOSSm9lryt5Y0cjvrCvt+x/v3kw3W2DYH+OPcuG1n5Xblhz7jA+vy+1tbX8+te/ZtmyZVx//fUsWbKEOXPm8NRTT3H77bdTXl7OI488gtfr7XfvwYMH8fv9vPLKK8yfPx+AUCjEpz/9aR588EE+/vGP88ILL+TUv7vvvps1a9YwadIkNmzYwNKlS3vLjDFs2rSJqqqqjG185jOf4YMf/CCPPvootbW1rFu3Lqc+JbN69WoALrroopR2KgSycTNXFCVH6ivqWVC7gCnlU1K6LY9nzjnnnLQu5pWV/feKnXXWWdx22210dXVx5ZVXcujQIa666ipisRhr165l9uz03pBf//rXe8UJrFHED3/4QyoqKnjxxRd59tkTX6WIRCLceeedgDV1lyhOYI0AzznnnLw7RGTDunXrePjhhykpKentayGS9QhKUZTciY+qnmIUp1B7e2deA9JmcjMvKSlJef2WW27h6aefZuPGjZx22ml0dHRw/fXXs3Jl5kA2V199db9rFRUVXHTRRTz00ENs2rSJM888M/uXALZs2UJzczNTp07lggsuGPiGYebJJ5/k+uuvR0T46U9/mvN05lCSVqBE5Gc5tGuMMdfncL+ijAsSp8mOdh+lLdB/C6EgeCMlFIXLcMf6e1oZDCGXH7+ru0/6DxiEB2BS6nnIHJUiHpDW4TJW7qo8rlGdiJt5fD1q1qxZdHR0MH/+fO6+++6M91RWVqYckQG9a0JHjhzJqh+JHDxoTdsW4hf/M888w8UXX0woFOIHP/hBSqEuJDKNoD6TQ7sGUIFSlCyoK6vr3feUuFZlMNZeKZcPV8xNUbiMokgpYjurxwXMGykhKhE7rFI3MYn1egCmDasU30OV4EwRj0oRi0aJBJxEE+L8xSJCqMtlCVVR3D3djOg+qg0bNtDdbYnskSNHaGpqoqFhMDEF0jMWp1+fe+45PvzhD+Pz+bjrrrv4/Oc/P9JdGpBMAvXZYeuFoih9iI96kkdVEUeYbm8bPm+7PaoqxR07vsDtNC7KQhWUhuIegN2EnYHesEpNXU2pNwCnikrhBE9plFhRGqGKe/0l7qMaZqHavn07N954Ix6PhyuuuIKHHnqIlStX8txzz6Xd19Pe3k5HR0fKNaD4XqO6ukFskE7DjBmWi//u3btPuI18s3nzZi688EK6urq44447+MpXvjLSXRoUaQXKGHPvcHZEUZT+xEdV7cH2Pq7qiaMqZ8xNcaQUb7i0N6q6AN5oMd5oMVGJEHBbdWMS7Q2rBCmiVaQbUdlCFQ30jZzexz09lVANYfR0n8/HlVdeid/v5/vf/z5f+MIXOHr0KE899RRf/vKXueeee9Le+9BDD7Fq1ao+1zo6Onj88ccBctowvHjxYmprazly5AgbN27k/PPPP+G24iIbiUQGqJmeF154gfPPP5/Ozk6+9a1v8bWvfe2E2xpu1ItPUUYBld5K5tfMZ0HtAqqK+ronRx1huj3ttJQepdPbQtjRdx3KaVyUhiqo7qljQqAWT/R4wsW2QBs7mnfQ2NGY9MBpVkBaz/HwOA4nuEujeCvDuLx9A9LGhSrY4SISSAhK29M8ZEFpb7jhBnbu3MlHP/pRbrrpJhwOBw899BAnnXQSP/jBD9iwYUPae2+77TZ27tzZ+zkcDnPjjTfS0dHB4sWLOeusE08y6Xa7ueWWWwD45Cc/2c9lPe5m3tHRMWBbEydOxOPx8NZbb9HWln2I05deeonly5fT2dnJ17/+db75zW9m3cZIckJefCJSBiwBJgKHjDF/y2uvFEVJS7pRFRiCrh6Crh6cMRdFEWutymFSjaqiBF0+Am4fUYn0rlVB0qgqRZw/hwMcpVFcxdbUXyRw/O/cPrH+vDGcRUlBaTOMqFavXp1xr89VV13F8uXLASsG3v3338+0adP4xS9+0Vtn8uTJPPDAA1xwwQV86lOfYtGiRb1TbnGmT5/O4sWLWbhwIeeeey4VFRU8//zzHDp0iNraWtavXz+I/4XMfPGLX2Tnzp38/Oc/Z+nSpSxZsoTZs2fT2trK66+/zuHDhzlw4MCAruZut5sVK1bw2GOPsWjRIs4880yKi4upra3t3ceUiQ996EN0dHRQWVnJoUOH0kZL/+53v5uXzb/5JiuBEpFy4HtYcffikS/vB/5ml/8zcAtwuTEmt51uiqJkJHEDcPJaVdQRwedpx+ex16oipX1GTk7jpCQ8gZLwBELOIAGXj5CrB5MusnoKoRIHuEuiuIqi/SJTmJhY0dMDTpzehHxUcaFKkYp+48aNGd934cKFLF++nF27drFq1SpcLhe/+tWv+kVBWL58OV/96ldZvXo1H/vYx3j66adxu48H6hURHnnkEVavXs0DDzzAwYMHmTBhAldffTW33357XsIIiQhr167l4osvZs2aNbzwwgu8+uqrVFdXM2fOHD7/+c8POjr72rVrqa6uZuPGjTzyyCNEIhFmzJgxKIGKj7ra29szBoT91re+VZAClU2w2BLgr8AioBl4GVgOrDPGfMqu04CVD+ouY8wtQ9LjAkODxSqZyCZ4Zr7Y27aXYDTY77o1qiqlKFLWO6pKpHddy+0j4uibJ6ifu3oK93RjIBpw2FN8/b3grBFVUoZfV1Fe91JlYigDryoW+Q4Wm80I6ktY4vQr4DpjjE9E+kxEG2P2icgeYORioCjKOCc+6kk9qurA5+nAEy2mKFyKJ1rcG1k9MV9VxBEm4LKyACe6q4M9BZjCmUIEK4VHUcyKnh5w9Kb5AIgEHUSCDjsVfQyHM/97qZSxRTYCdSXwJvBpY0wgQ72DwPwM5YqiDAPp16og5PQTcvpxGCfeSAnF4VKc5vg0mCvmpixUSWmokpDTT8Dt600DEp8CLHWXUl9Zb60pJQuVnRQxFhYi/uOJEwGiIQfRkAOn2xIzh2vk91IphUk2AtUAbBxAnMCa/iu8yUxFGadkWquKSRS/uwu/uwtXzENxuAxvpCRhE/Bxx4qYRAm4egi4uq3RWMKoqtRdSn3dIuhptUQGS6icHoPDHSEWsYUqkiBUYQfRsAOH0+AsSpE8cQhd1JXRQTYCFQYGE/J2KtA9YC1FUYaddNEqACKOEF3eVrq9bSk3ATuMk5JwOSXhcsKOIAG3j6CzByOm7xTghEnUOYr6CpXb4HRHiEaEqN8SpjixqBDzOYn403j+pXCoOBHq6+uHPeGgkhvZCNQbwCIR8Rpj+q/AAiJSCZwOaE4oRSlw4k4PmTcBxx0rSnsTKMLxlPVlVBF0+Qm6fISc1uRKW6CNNkA8XuocRVQGju/3cboMzvLUYZT6ev7FcHpth4pQt7VOlSehUkYP2QjUo8Cd9vGlNHXuAMqA3+TYL0VRholM+ar6OlYUUWQnV0x0rCiKlFAUKUmYAvQRdYSt8EoxP00ej+VYEfD1C6NkiqNEgg6iCZ5/xkDE9gZ0emK47HWqXqECnf4bJ2QjUP+JlajwJhFZgiVYADNE5LPAFcB5wA7g53ntpaIow0KiK3myu3rIGSDkDCDGYY+qSnAlRFdPnAKMOEIEXD29XoDxURUeD1NwURnqAey9VMWWCEVtLz+TwqGib3BaBrXxVxn9DFqgbLfy5VjC9H4gHgtkmX0I8CpwcbopQEUZjxhjRmV07Li7enuwnaNdRzHYU4ASS3CscOO1xSpxCtAV81AW8vR6AQZdPoKuAGBoIkKTx4MAdZEIlbFYf8+/QF+HinhwWnEYa53Km7Txdxj3UympGYr1vawiSRhjDgN/JyIXAR8GZgFO4DDwR+BRY0wsQxOKMq5wuVyEQqGCTak9GDJNAUYcYSJ2xApPtMiKWBEpTu0FGIwRdPX0bgQ2QJPLRRNQGotRH4n0ev45PbbnX8DRJy9VfJ0q7LfWqVxFup+qUAiHwzidzoErZsGgI0koqdFIEkom3n77bSKRCJMnTx6Vo6h0pNpbFUeMA6+9ETjRCzCR5AjricTFKk4sCtFg31BKiThctlDFp//i6PTfsNLc3Ew4HGby5MmDqj+YSBJpBUpEfgvcCzxhVMXSogKlZCISiXD48GFcLheVlZWUlJTgcDjGlFilywQMViR1b9iKBeg0qSdswo6gvV7Vg+kbnIaqaJS6qCVgxsSFytFn428ccRhrmtBjT//F0em/IcMYQzgcprOzk7a2NqZPnz7o2YJcBSqGlRn3TayAsOuMMXuy6/7YRwVKGYhYLEZbWxs+nw+/308sNnZnwTuCHX2mABNxxly4Y15cMQ+SIlafwRB1hIk4Q0RsL8BEqqIxiu0VhFhUiEUEE02VMcjgcFmHOJK+37xlUKzTf/nE6XRSXl5OdXV1VlPZucbi+wmwEqgDbgZuFpFngPuA3xhjegbdE0UZxzgcDmpqaqipqRnprgwbf9j/B77x7DcIxUL9ypxRN/Vt72LeO3/H1PZ5OOi/bhF2BGms3sae2pc4UrGLmCNGfH7PbQy3N7eywtdDqNtJ295S2veVEAv3F6vSSQEq5/RQXhfoO6qaeTZc8/u8va8yNGRcgxIRD3Ax8Cngg1gOEQbwAQ8DvzDGPDcM/SxYdASlKJm5Y/MdPLz74ZRlReFSGloWMad5MSd3zUpZx+/qZn/Nq+yp3cKx8gOW50WSWF3Y4aejsZi2PaUEO9z92nAVR6ls8FHZ0IO7OHEE64TL1sBpV+b6mkqW5DTFl6KxOqx9UJ8A5tmXDVaEifuAB4wxx068u6MTFShFGTyf3fhZNh/bnLKsPFDN7ObFzGleTLU/9UJ7p7eFvbUvs6d2C20l9teN/R221O/nZ8ea6XnLQ+ueUrqbioCkqUQxlE8JUDXHR8lJob5OFTqqGlbyKlBJDb8Xa1R1BTABS6iiWK7m9wGPG2Oi6VsYO6hAKcqJccljl7Cvc1//AgM1PVOY07yY2c3voSxU1b8O0FxyhL21L7Ov5hW6ilpJdPFb6vfz4/1ttO0roX1fCdFg/2lET3mYqtk9VMzswelJ+B50eOCSH+moaogZMoFKeEAxcDlwLdZm3TjvGGMGly5ylKMCpSi5k3ZkZYS6zgZmNy+moWUh3mhJyvvfKmtkX80r7Kt5BZ+3o49Yva/bz3de6qZ9byk97/RfxBdnjAnT/VTN6aG4Oty3cMmn4aL/yOndlNQMuUAlPexDwIPARMAYY/K7Y6tAUYFSlPyRybnCEXMyvf1U5jQvYUbrAlzGk6IFeLN8P/tqXmZ/zVZ6PJ19xGrFkQA3/S1Ex4FiYpH+ThVF1SGqZvuYMD1gxf+LU3sKfO5vub+g0stwjKDKsDz9rgXex/EJ30PGmPoTbngUoQKlKENDJucKd8TLzNbTmN3yHqZ0zMOZ4u9hQ4yjE/axr+YV9tdsJeDu7hUrb8jwiVdDXPSyIdje36nC4YoxYYafyoYeiqrCfdeqdFSVF4ZyDeoc4JPAZUAxljAFgd9jrUH973jZ3KsCpShDTybnCm+4hJmtp9HQsogpHXNSuq3HiHK0Yi97a17mQPVrBN09llgZw5yj8E8vhlnwRt8U9b3tV4apbPBRMcOva1V5JN9efDOxvPiuAaZzfLT0KvAL4EFjTOrt5GMYFShFGV4Gcluf1bKQhpZF1HU2IPSfxotKlKaK3eyteZnG6u2EXH4whvIew7JthhWvxqhuSxGpwhljwrQAlbN6KJ4Y0lFVjuQsUCJSguWpdy1WBHOxjzbgl8C9xphX89Xh0YgKlKKMHJlGViWhCcxqOZ2GlvcwOc0eq6hEaKp4g/3VW2ms3kbA7YNYjFOOwHmvGs7cFcOVIqySZ0KYylk9VNT7cRUlRQZRsRoUuYY6uhdLnEqxRCkGPIk1hfeYMab/KuY4RAVKUQqDTGJVGqygoWURDS2LmNRdn7JOjChHJ+zjQM1WDlS/Ro+nk1J/jLN2GM7baqh/O8VNDmtfVWVDD6WTgrqvKgvyEYsP4ACwDitqxJG89nAMoAKlKIVFJk9AsDYEN7QsZFbLIk7yTU9ZxxDjWHkjB6q3sr9mK92eVhreNJy71XDW64biFE27SiJUzvRTUd+DpzxpG6iOqvqRq0A9ANxnjHlqKDo3VlCBUpTCZTBiNbP1NGa1np421BLA26WH2F/zKvurtxJ0vsP7Xo9x3lbD3KOp65dMDFIxs4cJ0wI43InfsRpaKc6w7oMar6hAKcroIW30Cqw1q5mtpzGr5XQmd87GkcLBAqC5pKl3ZFXW/Sbnbo3xge2G8kD/uuKKMWFqgIqZPRpaKQkVqGFABUpRRieZxKooXMrM1tOY2XI6UzrnptxnBdBe9DaNVds4XPEaU48dYNn2GAv3GZwpvlbdpREqZlqOFZ4ynQJUgRoGVKAUZfSTycHCEymmvu1dzGw5nWntp+Ay/Tf2AvhdXTRW7eCdktdoOLSLZduCTGtO/bySkxKmAF3jc2+VCtQwoAKlKGOLTGLljnqZ3jafWa2nM71tftqU9mFHkCMVu/A5tjGvcRtn7uimLMUUoMMVo3xagIr6FFOAYzy8kgrUMKACpShjl0xi5Yy5mNIxl/rWd1Pf9i5KwhNS1osR463y/cTCr3HK/tdYvOcdkhP9gpWzqmJGDxPq/RRVRvoWjsH1KhWoYUAFSlHGB5nECiNM6p7RK1ZV/vTJHNq9R3EGX2PBnq00NB1G6P8d7K0MU1Hfw4TpftwlY3MjsArUMKACpSjjj4xiBVT4T2Jm67uob3s3k7rqU4ZcAgg6OnAFdjC3cRtT3tqNKxrsU24wlE4KUVHfQ/nUAE530vf1KBYrFahhQAVKUcY3mWIDAhSHypnRtoD6tncztX1eWieLGGGcwb3MPLKdSe9sozjQ0qc86IK36yO8t66DspODSKLmjULnChWoYUAFSlGUOH/Y/wdufeZWIiaSstwV9TCt/RTq297N9Lb5FEfK0rblDL1J3bEd1LZso6JzPw5zfKqvowT2z4ly6UmtFFWHR6VzhQrUMKACpShKOjLttRIjnNQ9gxltC5jeNp/anqlp23FEe6hteZ2alh3UtO7AE/b1lr1VAZtPhcWTO/mgp3vUiJUK1DCgAqUoymAYaCqwNFhpiVX7fKa2z02bMRgTY0JnI7Ut26lpfZ2y7iO9jhZHauDZUwXfrCD/HkrahFVgnoAqUMOACpSiKNkykFg5o26mdM6xR1cLKA9Vpa3rDnVS07qT6tbXqW7bhSfcDcD+SfDsfOHFU+CmUBsrfD3HbyoAsVKBGgZUoBRFyYWBAtpioLpnMjPaFzCjbQEnddWnjROIiVHedYia1tepaX2d8q6DOEyMXVOtkdXmU4RTHQHWvpUwuhohsVKBGgZUoBRFyScDubAXhUuZ1n6qfZyS0dHCFe6hum0X1bZgucMdbJ9hidUL8wRfkbCys4tbW9utG4bRbV0FahCIyFXAPwOnAU5gF1YK+58YY2KZ7gUVKEVRho6BR1fCRN/UXsGalGl0BZR2N/VOB5Z17WN7fZTN8+DFOYKvWHAbw+3NrdZ04BCPrFSgBkBEfgSsAgJY2YLDwHlAOfAYcIUxJpq+BRUoRVGGj4FGV55IMVM75vYKVlmoMm1dRzRIVfteqtp2UdG+i/0T32TzqcKW2dBVYrkCLvX7renAIRArFagMiMg/AL8FjgEfMMbssa9PAp4CTgVuMsbck6kdFShFUUaCgRwtMFDlP5lp7acyvf1UJnc24DSutNU9oU6q2nZT2baLt0p38+KcDl6Ye1yswBaskgV5ESsVqAyIyBZgMXCNMWZ9UtnZwCYs8ZqSaapPBUpRlEJgoNGVK+qhrnN2r2BVBCZmbK+k5xiBx+T2AAAN/klEQVSVbbvpdu1m67S9bD4lQEdJ3zorpZxbr3n+hPqrApUGEZkKHAZCQKUxxp+izhFgCnCmMea5dG2pQCmKUmgMOLrCSnc/tWMeUzrmMq19Lt5oemcLyzvwILHoLvZOfIO/zjtI64QoGAPGsHrmZaxYdntWfVSBSoOIfAT4PfCKMeY9aeo8BlwCfM4Y86N0balAKYpS6Aw0usIItb4pTO2YS33LPE7yzcJBmo3CWOtXruBejpW9wfMN2zhc/U7WIjUYgUo/ITm2mWmfD2aocyiprqIoyqhk7flre//dzzPQGBBoLjtCc9kRXp3yZ5wxFyd3zaLh7bk0tMzFG5tGYnTamNNLqGQB1bEFnL3L8OBZf+Ge/Y9lPYoaiPEqUPGxrC9DnW77XJ5cICLXAdcBTJ8+Pb89UxRFGUJWzFrBilkrej/3mQ60Z9SijghNFW/QVPEGT88Bb7iEOW/P4bQjc6kKziPqPr5+9VL9bgCOpfduP2HGq0DF3VJOaH7TGPMz4GdgTfHlq1OKoijDza1Lb+XWpbf2fu6dDkxY/gm6e9g+ZSvbp2wFoK6tiqV753GSbwY7pxwD4OQBd41mz3gVqC77nGFVsLesK0MdRVGUMUXidCDYgvXmZnr/nhfhaFUbvztjM2Cva8Vi3Djrsrz3ZbwKVKN9npGhzrSkuoqiKOOO5PWrW5++mUjC5JPDGO48AS++wTBeBeoV+7xARIpTuZkDZyTVVRRFGdckr18NNUOwrFX4GGMOAy8DHuCK5HJ7o+5UrI26J7YLTVEURcmJcSlQNv9mn+8SkdnxiyJyEvBj++PqwQSMVRRFUfLPeJ3iwxjzWxH5CVYk820i8ieOB4udAGwAfjiCXVQURRnXjFuBAjDGrBKRZ4AbgLM5nm7jPgaZbkNRFEUZGsa1QAEYY34J/HKk+6EoiqL0ZVzG4ssnIvIOmUMmZaIWaB6wlpKI2iw71F7ZofbKjlzsNcMYkzGkugrUCCIiWwYKlqj0RW2WHWqv7FB7ZcdQ22s8e/EpiqIoBYwKlKIoilKQqECNLD8b6Q6MQtRm2aH2yg61V3YMqb10DUpRFEUpSHQEpSiKohQkKlCKoihKQaIClSdE5CoR+auIdIhIt4hsEZEbRGTQNhYRt4icJyLfE5HNIvKmiIREpElEfisiy4bwFYadfNgsQ9t3ioixj/+Tj/6ONPm2l4gUi8hXRORFEWkXkR4ROSAivxGRM/Pd/+Emn/YSkaki8p8isltE/CISEJE9IrJGRGYNRf+HCxGZJyI3isiDIrJLRGL2783lObabu/2NMXrkeAA/wsrm5QceBx4DOu1rvwOcg2zng/Y9BnjTbuthYFvC9dtG+n0LyWZp2j4DiAAxu73/M9LvW2j2AmYCe+z73wL+C3gEeAEIAbeO9DsXir2ARUCbfe9hrDidG4Aj9rUu4H0j/c452OruhO+XxOPykbb/iBtntB/APyQIypyE65OA1+2yGwfZ1rnAb4H3pyhbaX/pGuCckX7vQrFZira9wA6gyf6lGPUClW97AaXA3vgfPIA7qbwGmDvS711A9nrOvudnibYC3MC9dtnWkX7vHOz1GeA7wJVAA7ApF4HK63fiSBtntB/AFtvgn0hRdnbCf5QjD8/6ud3evSP93oVqM+Au+/6PAOvGiEDl1V5YqWYMcP9Iv1uh2wso4viI4uQU5XUJ5SUj/e55sl+uApU3++saVA6IyFRgMdaUyG+Sy40xf8H6S/5kYGkeHhnP7js1D22NCENpMxH5e+BLwC+NMf+de29HnnzbS0Q8wGftj6vz19PCYAh+vqJYMxcAkqI8vk/HhzWdNa7Jt/1VoHJjkX3eYVKnjQd4MaluLsyxz2/moa2RYkhsJiJFwP1AK3DjiXev4Mi3vRZjTeEdNsbsFJH32Q4lPxWRb4vIe3Pt8AiTV3sZY8LAk/bHb4uIO15m//sO++O9xh4ijHPyav9xn24jR2ba50zRzA8l1T0hRORk4Fr746O5tDXCDJXN/hWYB3zMGDOWolHn217vts97RGQdcE1S+TdE5FHg4xm+YAqZofj5WgU8gTXyvFBEttjXzwCqgHuAL2fZz7FKXu2vApUbZfbZl6FOt30uP9GHiIgLeBCoAJ4c5dNXebeZiLwPuAnYYIx5OIe+FSL5tle1ff4AVoLO7wJrgBb72o+xFrk7gU9l29kCIO8/X8aY/fbP2HrgQvpOsW8BnrZHWkqe7a9TfLkRn5Me6qH9GqxU9IeBq4f4WUNNXm0mIsXAL7C+UFflo80CI98/Y/HfeRfWtNSXjTH7jDHtxpjfA5fYz7pmlO7vyfvvpC1O24HZwMVYOZAmYtmqCnhURL6Rr+eNcvJqfxWo3Oiyz2UZ6sTLujLUSYuI3AN8GjgGnGeMOXYi7RQQ+bbZncBc4F+MMaN5bS4d+bZXYp21yYXGmC3AS1jfDcsG0V6hkVd7iUgl1p6ncuACY8zvjTEtxphmY8x/ARdgOUd8XUTmZGprnJBX+6tA5UajfZ6Roc60pLqDRkS+B3wBeAdLnPZk20YB0mif82WzS7E25F4jIpsSD6wvD4B/tq/9/AT6O9I02ud82SuxzoE0deLXTx5Ee4VGo33Ol71WYI2WNhtj9icXGmP2An/DGpEuG2wnxzCN9jkv9tc1qNyIu30vEJHiNIvKZyTVHRQi8h3gX7DWBj5kjHn9xLtZUAyFzRxY+yvSMcs+KgfZXiGRb3u9nPDvGqw/fpKptc/dKcoKnXzba7p97shQp90+V2eoM17Iq/11BJUDxpjDWL/wHuCK5HIRORtrQfUY8Pxg2xWR1VheQW1Y4rQ1Lx0uAPJtM2NMvTFGUh1YbucAX7avLczfmwwPQ2CvJqy/+MFa10xurwp4j/1xS3J5oTMEv5NH7fPiRBfzhPbcWK77kH5EOm7Iu/1HetfyaD+Ayzm+M3p2wvWTsELu9AvrgbWTfxfwbynau92+pw1YPNLvNxpsluE56xgbkSTy/TP2EY7H4FuYcL0I+LVdtgU7X9xoO/JpL/sen33PDwFvQpkX+Ild1gpUjPS758l+mxggksQAP19Z2z/tc0baGGPhwHLNNViLpf+NFQyxw772GEmBERO+ONclXf8ox8OmvGjXS3XcPNLvXCg2G+AZY0KghsJewL/b5UHgabuNJvvaERJiqI3GI5/2wtorFo+D2QT83m7zqH0tAFwy0u+cg63eA2xOOOJBXd9IvJ7lz1dW9k936BpUHjDGrBKRZ4AbsNZCnFh/XdwH/MQYExtkU4lz2EvsIxV/YZSHqcmjzcYF+baXMebLIvIc8HmsHf0lWBso/wNYbYxJtTY1asinvYwx94vINqy9du8HlttFTVjBYv/DjO414gnA36e4fsJeifmyv6Z8VxRFUQoSdZJQFEVRChIVKEVRFKUgUYFSFEVRChIVKEVRFKUgUYFSFEVRChIVKEVRFKUgUYFSFEVRChIVKEUZA4jIIhExInLfMD3vAhH5XxFpFZEeEdkuIl8TEe9wPF8ZH6hAKcrY4DL7/NhQP0hEvgL8ETgXKzDoH7DirN0BbBKRkqHugzI+0EgSijIGEJEdWKkhJhpjAkP4nCXAC1gx1s41xvzNvl6GJVQfAO42xnxxqPqgjB90BKUooxwRmQvMB/44lOJkczNWWu+74uIEYIzpBj6JlTxylZ2JVlFyQgVKUYYBe33I2P++VkS2iIhPRI6JyL0iMtEuKxKRb4vIGyISEJFDIvKvqXIRJfAP9vl3Cc9bZj9zk93m7SKyV0T8IrJfRG4VEaddd5rdhyb7mdtE5OoU7+ABLrQ/PpRcbqyMs89j5QL68AmYSVH6oAKlKMOIiNwF/BQrf9ATWOkHPgX8yZ4mexIrwvgO4M9YWW//L/CjDM1eCoSwptiS8QD/Dyuq9GtYuX4mYeUd+6GINGCldjkH+Kv973cBD4jIPyW1NQ8r6nmrMWZfmr68aJ8XZeivogwKTbehKMPLNVhJAndCbwbb54HT7HM7MNMY02GXL8T60v+MiPyrMeZgYmMiMhUrLcsfjTFdKZ73XuCZpDZPt9u8DisVwq+BLxljonb5DVjJ+b5N35HSTPt8KMP7xctmZqijKINCR1CKMrx8Iy5OAMaYNmCN/XE+cF1cSOzyV4H/wVr3OTtFe5fZZb9LUQbWmlBym1vtNh1YI6KvxMXJJj7CaxCR6QnXy+yzL8P7ddvn8gx1FGVQqEApyvDyRIpre+3zwUTxSmCPfa5LUXYZEMXK8pqKdG3Gn/lnY0woscAYEwEOpHimxKukeZai5BUVKEUZXo6kuNadoSyxvCjxoojUAmcBz2TIgDtQm9k8Mz6FWEZ64mWpphsVJStUoBRlGBkg1XW2ae4vwUqlnW56bzBtZvPMRvs8PUOdaUl1FeWEUYFSlNHLpfZ5wzA9bxfWBt1q2/svFX9nn18Zni4pYxkVKEUZhYhIOXAesMUYk8mrLm/Ya1V/tD8mu6AjIrOwvAbTubwrSlaoQCnK6OQiwMswxN5LYjWWk8RXRSQ+WoqHOroP6zvlx8aY9mHulzIGUYFSlNFJPDhspvWnvGOMeREr3FEJ8Jwd0fwRYB+WG/zfgK8NZ5+UsYsKlKKMMkSkCLgA2GmM2TXczzfGfAcr5NFTwBnAR4Bm4FbgbGNMz3D3SRmbaDRzRRlliMjFWI4RdxpjdLSijFk01JGijD78WGGIHhzpjijKUKIjKEVRFKUg0TUoRVEUpSBRgVIURVEKEhUoRVEUpSBRgVIURVEKEhUoRVEUpSBRgVIURVEKEhUoRVEUpSD5/yJaQmduz1kPAAAAAElFTkSuQmCC\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "numsol_heun2=np.zeros([N,3])\n", | |
| "numsol_rk22=np.zeros([N,3])\n", | |
| "\n", | |
| "numsol_heun2[0,0]=y0\n", | |
| "numsol_heun2[0,1]=v0\n", | |
| "numsol_heun2[0,2]=m0\n", | |
| "numsol_rk22[0,0]=y0\n", | |
| "numsol_rk22[0,1]=v0\n", | |
| "numsol_rk22[0,2]=m0\n", | |
| "\n", | |
| "for i in range (N-1):\n", | |
| " numsol_heun2[i+1]=heun_step(numsol_heun2[i],rocket,dt)\n", | |
| " numsol_rk22[i+1]=rk2_step(numsol_rk22[i],rocket,dt)\n", | |
| " \n", | |
| "plt.plot(m/m0, v, label='Analytical')\n", | |
| "plt.plot(numsol_heun[:,2]/m0,numsol_heun[:,1],'-o',label='Implicit')\n", | |
| "plt.plot(numsol_heun2[:,2]/m0,numsol_heun2[:,1],'-o', label='Implicit 2')\n", | |
| "plt.plot(numsol_rk2[:,2]/m0,numsol_rk2[:,1],'-',label='Explicit')\n", | |
| "plt.plot(numsol_rk22[:,2]/m0,numsol_rk22[:,1],'-', label='Explicit 2')\n", | |
| "plt.xlabel('m/m0')\n", | |
| "plt.ylabel('Velocity (m/s)')\n", | |
| "plt.legend(loc='best')" | |
| ] | |
| }, | |
| { | |
| "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": 9, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "def f_m(dmdt,m0=0.25, c=0.18e-3, u=250):\n", | |
| " ''' define a function f_m(dmdt) that returns \n", | |
| " height_desired-height_predicted[-1]\n", | |
| " here, the time span is based upon the value of dmdt\n", | |
| " \n", | |
| " arguments:\n", | |
| " ---------\n", | |
| " dmdt: the unknown mass change rate\n", | |
| " m0: the known initial mass\n", | |
| " c: the known drag in kg/m\n", | |
| " u: the known speed of the propellent\n", | |
| " \n", | |
| " returns:\n", | |
| " --------\n", | |
| " error: the difference between height_desired and height_predicted[-1]\n", | |
| " when f_m(dmdt)= 0, the correct mass change rate was chosen\n", | |
| " '''\n", | |
| " h = 300\n", | |
| " m0 = 0.25\n", | |
| " mf = 0.05\n", | |
| " N = 100\n", | |
| " time = (m0-mf)/dmdt\n", | |
| " t = np.linspace(0, time, N)\n", | |
| " dt = t[1] - t[0]\n", | |
| " num_h = np.zeros([N,3])\n", | |
| " num_h[0,0] = y0\n", | |
| " num_h[0,1] = v0\n", | |
| " num_h[0,2] = m0\n", | |
| " for i in range(N-1):\n", | |
| " num_h[i+1] = heun_step(num_h[i], lambda state: rocket(state, dmdt=dmdt),dt)\n", | |
| " calc_h = num_h[:,0]\n", | |
| " error = h - calc_h[-1]\n", | |
| " return error" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 25, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "def incsearch(func,xmin,xmax,ns=50):\n", | |
| " x = np.linspace(xmin,xmax,ns)\n", | |
| " f = np.empty(ns)\n", | |
| " for i in range(ns):\n", | |
| " f[i] = func(x[i])\n", | |
| " sign_f = np.sign(f)\n", | |
| " delta_sign_f = sign_f[1:]-sign_f[0:-1]\n", | |
| " i_zeros = np.nonzero(delta_sign_f!=0)\n", | |
| " nb = len(i_zeros[0])\n", | |
| " xb = np.block([[ x[i_zeros[0]+1]],[x[i_zeros[0]] ]] )\n", | |
| " if nb==0:\n", | |
| " print('no brackets found\\n')\n", | |
| " print('check interval or increase ns\\n')\n", | |
| " else:\n", | |
| " print('number of brackets: {}\\n'.format(nb))\n", | |
| " return xb" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 26, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "def mod_secant(func,dx,x0,es=0.0001,maxit=50):\n", | |
| " iter = 0;\n", | |
| " xr=x0\n", | |
| " for iter in range(0,maxit):\n", | |
| " xrold = xr;\n", | |
| " dfunc=(func(xr+dx)-func(xr))/dx;\n", | |
| " xr = xr - func(xr)/dfunc;\n", | |
| " if xr != 0:\n", | |
| " ea = abs((xr - xrold)/xr) * 100;\n", | |
| " else:\n", | |
| " ea = abs((xr - xrold)/1) * 100;\n", | |
| " if ea <= es:\n", | |
| " break\n", | |
| " return xr,[func(xr),ea,iter]" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 27, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "number of brackets: 1\n", | |
| "\n", | |
| "3a. The two closest mass change rates within that interval are [0.08571429] and [0.07857143]\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "rate=incsearch(f_m,0.05,0.4)\n", | |
| "print('3a. The two closest mass change rates within that interval are',rate[0],'and',rate[1])" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 28, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "The root of f_m is 0.07912698052690133\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "root= mod_secant(f_m,dx=0.0001,x0=0.0857)\n", | |
| "print('The root of f_m is', root[0])" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 35, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "<matplotlib.legend.Legend at 0x7f7ebf7a6748>" | |
| ] | |
| }, | |
| "execution_count": 35, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaYAAAEaCAYAAABaefMNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3xUZfb48c9JIZXQewIJoIKAqIBio1oQXVZX0V1cKRbssrruiisrqLsKru6uinVXRBfZ/X3FighWQBEVQ5GOIC30TkIKKXN+f9zJpDAzZJJJZpKc9+s1r8u993nunCHl5D73KaKqGGOMMeEiItQBGGOMMaVZYjLGGBNWLDEZY4wJK5aYjDHGhBVLTMYYY8JKVKgDqO2aN2+uqampoQ7DGGNqlaVLlx5Q1RbezlliqqLU1FTS09NDHYYxxtQqIrLN1zlryjPGGBNWwiYxicg9IvJ/IrJORA6KSIGI7BeRz0XktyIifuqOEJGvReSoiBwTkXQRuUtE/H6+ytYzxhhTfcKpKe9BoCWwGlgMZAMdgEHAYOBaEfmVqrpKVxKRF4A7gTzgC6DAXX4qMFhEhqtqUfk3q2w9Y4wx1SucEtOvgeWqml36oIh0w0kcvwRGAa+XOncNTnLZA/RT1Y3u462A+cDVwN3As+WuWal6xhhjql/YNFmp6qLyScl9fA3wgnv3knKnH3JvHyxOLu46e4E73LvjvTTNVbaeMcaYalZbfvEWurd5xQdEJBnoBeQDb5evoKoLgZ1Aa6BvVesZY4wpJWsPvH45ZO0N+qXDPjGJSBpwu3t3dqlTZ7m3a1Q110f1H8qVrUo9Y4wxxRY+Bdu/g4VTgn7pcHrGBICIjAH6A9FAMnA+TgJ9UlXfK1U0zb312Rce2F6ubFXqGWOC5MUFm2gQGcF1fVJIio0OdTgmEH9pCYXHS/bTX3NeUTEwYV9Q3iIc75guwOnkMALo5z72Z+CxcuUS3dsTnkuVcsy9bRiEeh4iMtbdtTx9//79fi5jjCkvM6+AF77cxF/mrOO8J75g8/5jJ69kwse4ldB9OETFOftRcdBjOIxbFbS3CLvEpKq3qKoA8UA34J/AJOA7EWlbqmjxuKZAVzqsbL3SMb6qqr1VtXeLFl5n1DDG+PB2+g6y852RGG0bx5HWPCHEEZmANGwNMQ2h6DhExTrbmCRo2CpobxF2iamYquaq6lpV/QNOL7qeOGOMimW5t4knVC5RfC6r1LHK1jPGVFGRS5m+eItnf/QFqfgZO2/CVfY+6DUGbvnc2R4LbgeIsHvG5MPrwNPAL0QkWlULgK3ucx381Etxb7eWOlbZesaYKvps7R4yDjl9jhrHR/Ors5JDHJGplF+/VfLvK/8e9MuH7R1TOUdwuoxHAU3dx5a7t91EJM5HvT7lylalnjGmiv79dcnd0g3ntieuQWQIozHhqrYkpn44SekIcABAVTOAZUADYHj5CiLSH6dX3x7g2+Ljla1njKmaFRlHSN92GIDoSGHkeamhDciErbBITCJykYjcICIxXs5dALzm3n2t3Px1T7q3U0Skc6k6LYEX3buTy8+vV4V6xphKem1Ryd3SL85oS6uk2BBGY8JZuDxj6oTzHGmqiCzDuVtp6D5+urvMHJxu4x6qOktEXsKZRmiViHxOyWSsScD7lO0wUaV6xpjK2Xkkl49X7fbs33ShDRE0voVLYloIPA5cBJyKM6hWcBLUO8AMVX3fW0VVvVNEFgF34QzMjQTWA9OAl3zd9VS2njEmcG8u3kqRyxmh0bdjU7q3axTiiEw4C4vEpKpbgEeqUH8mMLOm6hljKi77eCEzl2z37N9yYccQRmNqg7B4xmSMqbveTs8gK8+ZhzmteQKDurQMcUQm3FliMsZUmyKXMu2brZ79my5IJSLCBtQa/ywxGWOqzWdr97L9UA4AjeKiuaaXDag1J2eJyRhTbaaV6iI+4tz2xDcIi8faJsxZYjLGVIuVO46wZOshAKIihFE2oNZUkCUmY0y1KD390JVntKF1IxtQayrGEpMxJugyDuUwp9SA2puti7gJgCUmY0zQvbZoi2dA7Xkdm9Ej2QbUmoqzxGSMCarD2fn8vx8yPPu3D+gUwmhMbWSJyRgTVP/5bhu5Bc5cy11aN6TfKc1DHJGpbSwxGWOCJje/iOmLt3r2b+/fyVaoNQGzxGSMCZpZSzM4lJ0PQLvGcVxxRpsQR2RqI0tMxpigKCxy8a9SXcRvvjCN6Ej7FWMCZ981xpigmLdmj2f6ocbx0fz6nJQQR2RqK0tMxpgqU1VeWbjZsz+ybwebfshUmiUmY0yVLf75IKt2HgUgJiqCkeenhjYgU6tZYjLGVNnLC3/2/Ht472SaJ8aEMBpT21liMsZUyZpdR/l64wEAIgRuvcimHzJVY4nJGFMlL5d6tnR5jzZ0aJYQwmhMXWCJyRhTaZv3H2POyl2e/dv72fRDpuosMRljKu2lBT/jnquVfqe2sMlaTVBYYjLGVMqOwzm8t3ynZ//ugZ1DGI2pSywxGWMq5dWvNlPovl06J60p56Q1DXFEpq6wxGSMCdi+rDz+V2ppC7tbMsFkickYE7DXvt5CfqELgDOSG3GRLW1hgsgSkzEmIIez85nx3TbP/l0DO9vSFiaoLDEZYwLy+uKtZOc7CwGe2iqRS7q2CnFEpq6xxGSMqbCsvAKmf1OytMVdAzsTEWF3Sya4LDEZYypsxnfbycwrBCC1WTxXntE2xBGZusgSkzGmQvIKinhtUcn0Q3cM6ESk3S2ZamCJyRhTITO/386BY86y6W0bxXL1WckhjsjUVZaYjDEnlVdQxEullra4rX8nGkTZrw9TPew7yxhzUm99v539WccBaJ0Uy/V9bNl0U30sMRlj/MrNL+KlBSV3S3cO7ERsdGQIIzJ1nSUmY4xfb32/jQPHSu6Wruttd0umekVVtKCIJAEXA4OAs4BWQGPgMLAPWAbMBz5X1czgh2qMqWm5+UVlFgK8y+6WTA04aWISke7APcAIIB4o3z+0KdAJOA+4E8gRkbeAqaq6OrjhGmNqUum7pTaNYrnOni2ZGuAzMYlIK+AJYBROk99e4CPgW2AtcAjIBJKAZsDpOMlpADAWuEVEpgMPq+reavsExphqkZNfyMsLSz9b6kxMlN0tmern745pI5AAvAdMA+apqstP+c+AZ0UkErgcuAkYA1yL0+RnjKlF3vqu7Lil63rbuCVTM/x1fpgP9FTVa1X145MkJQ9VLVLVj1T1V8CZwIIgxGmMqUF2t2RCyecdk6r+sqoXV9VVwFVVvY4xpmb959ttHMwufbdkz5ZMzbHu4saYMrKPF/LKV6V64g3qbLM8mBpl323GmDKmL97KIffdUrvGcQzvZXdLpmZVeBxTMRGJBXoDbYFYX+VU9c0qxGWMCYGjOQVlni3dY3dLJgQCSkwiMh4YDzSsQHFLTMbUMi9/9TNZ7vWWOjZP4Npe1hPP1LxAZn64D2dcE8BKnO7kx6ojKGNMzduXlcfrpVanvf/SU4mKtLslU/MCuWO6AygArlLVudUUjzEmRF74chN5Bc6okNPbJDG0e5sQR2Tqq0D+HGoPfGVJyZi6J+NQDjOXbPfs/+Gy04iw1WlNiASSmHbjTENkjKlj/vn5RgqKFIDeHZow4LQWIY7I1GeBJKYPgAtFpEGwgxCRaBEZLCLPiMh3IrJbRPJFZKeIzBKRASepP0JEvhaRoyJyTETSReQuEfH7+Spbz5i6ZOPeLN5bvsOz/8chXRCxuyUTOoH8Ap4EZANvikjTIMfRH/gcuB/oACzFmaPvEHANMF9EHvNWUUReAN7C6cL+Nc6cfacCU4FZ7rn7glbPmLrm75/9hMu5WaL/qS04Jy3YP97GBKbCnR9U9YiInAssBDaLSDqwA/A2h56q6s0BxOEC3gGeVdWvS58QketxEsifRWS+qs4vde4anKU29gD9VHWj+3grnLn+rgbuBp4td81K1TOmrlm54whzV+/x7D9w6WkhjMYYh6hqxQqKxAGzgCGcuCZTeaqqQbvjEJF/AzcD00onPHdy7AWMKj+gV0T640wguwdoV3oS2srW86Z3796anp5e+Q9nTIioKje+toRFmw4AMLRHa168oVeIozL1hYgsVdXe3s4F0l38LzjLWRwEZgCbqLlxTMvdW89oPxFJxkku+cDb5Suo6kIR2Qm0A/oCi6tSz5i65quNBzxJKTJCuP8Su1sy4SGQxHQdzjLqZ6rqzmqKx5dT3NvdpY6d5d6uUdVcH/V+wEkwZ1GSYCpbz5g6o8ilPPnxOs/+9X1S6NwyMYQRGVMikM4PzXDGMdVoUhKR1sBo9+47pU6lubfb/FQvHpiRVupYZesZU2e8t3wn6/dkARDfIJLfXXzKSWoYU3MCSUybgRrtqSYiUTjNho2AL1R1dqnTxX/eZfu5RHFTY+m5/Spbr3RcY91dy9P379/v5zLGhJ+8giKe+XSDZ//WizrSsqHP+ZiNqXGBJKZpwAB3z7Wa8jIwGMgAflvuXHEHjIr13qh6PQ9VfVVVe6tq7xYtbCCiqV2mfbOF3UfzAGieGMPYfh1DHJExZQWSmP4BzMEZUzS4ugehisizOD3x9gCDVXVPuSJZ7q2/hvHic1mljlW2njG13qHsfF6aX7KsxX2XnEJCTMCr3xhTrQL5jiz+bu4AfAoUiMgefI9j6lTZoETkGeBeYD9OUtropdjWUvH4UrzC2dZSxypbz5ha7/kvN5J13L2sRYsErrcl000YCiQxpZbbb4Azsas3lW4mE5GncGaAOAhcoqprfRQt7kLeTUTifPSw61OubFXqGVOrbTuYzYzvSvr8jB/SxZa1MGEpkMRU7T3URGQy8AecbumXqOqPvsqqaoaILAPOBoZTbmFC90DZZJymwG+rWs+Y2u6pTzZ4Jmrtk9qES06vycfFxlRcIFMS+eteXWUi8jjwIHAEJylV5G7lSZxBslNEZLGqbnJfqyXworvMZC+zN1S2njG10tJth5mzsmQY4ENDu9pErSZshcVTTxEZBkxw724C7vHxQ7NeVScX76jqLBF5CWcRw1Ui8jnOYoaDgSTgfZxJWcuobD1jaiOXS3ls9hrP/hU92nB2+yYhjMgY/3wmJhGJVtWCqr5BBa9Tejrj3u6XNwuByaUPqOqdIrIIuAtnlvJIYD1O9/aXfN31VLaeMbXNe8t38uOOowA0iIpg/OVdQhyRMf75u2P6SUQmAv/Ris70Woo4tzyjgImc5PmUqk4Hpgf6HqXqzwRm1lQ9Y2qL7OOFTJm33rN/60VppDSND2FExpycvy45x4DXgY0i8rCI+OqBV4aIdBCRPwMbce4+MqsepjGmMl5a8DP7so4D0KJhDHcM6BziiIw5OX93TD1x1iyaBDwOPCYim3B6qq3D6c6difNMphlwOnAe0AlndoWDOGsavVxNsRtj/NhxOIdXv97s2f/jZaeRaINpTS3g87vU/Yxlqoi8jjOJ6p1AV5yZvr017RX3VliJ07NthqrmBDVaY0yFPTl3PfmFzqPSM5Ibcc3ZySepYUx4OOmfT6qaDbwAvCAinYEBwJlAS5zJVY8A+4BlwHxV3VpdwRpjKmbJlkNluoc/cuXpRERY93BTOwR0X+8e77OpmmIxxgSBy6U89lFJ9/Bf9GxL79SmfmoYE15sPhJj6phZS3eweqfT5yjGuoebWsgSkzF1yNGcgjLdw2/r15F2jeNCGJExgbPEZEwd8sxnGziYnQ9A20ax3D6g0pP8GxMylpiMqSNW7zxaZvbwR35xOvENrHu4qX0sMRlTB7hcyp8/WI3LPZCj36ktuKxb69AGZUwlWWIypg6YtWwHy7cfASA6Upj0i9Nt9nBTa1liMqaWO5pTwOS5JR0exvbrSMcWiSGMyJiqscRkTC339KcbOOTu8NCucRx3DbT58EztVuHEJCKbRWRKBco9KSI/Vy0sY0xFrNpxlBnfl3R4+POV1uHB1H6B3DGlAi0qUK65u6wxphq5XMojH66meFGa/qe24LJutly6qf2qoykvDiishusaY0qZuWS7p8NDg8gIJg3rZh0eTJ0Q1MQkIo2AC4A9wbyuMaasvZl5TCnV4eH2/h1Ja54QwoiMCR6/jdEisrncoWtFZICfa7Vyb1+remjGGF8mfbiGrONOw0TH5gncaR0eTB1ysqekqaX+rUCi++VLPvA+8GDVwjLG+PLZ2r3MXV3SKPHXq3sQGx0ZwoiMCa6TJaY091aAzcAs4A8+yuYD+1XVni8ZU02OHS/kkQ9We/av653MeZ2ahTAiY4LPb2JSVU8/VBF5A/i69DFjTM16+pMN7D6aB0CzhAb8aWjXEEdkTPBVeMCDqo6pzkCMMf6tyDjCG99u9ew/8ovTaRzfIGTxGFNdbOYHY2qBgiIXD727yjNmqd+pLRjWs21ogzKmmgQ0RFxEkoC7gMFAWyDWR1FVVVsIxpgg+ffXW1i321mVNjY6gr9e1d3GLJk6q8KJSURSgK+BFJzOEP5oVYIyxpTYtC+Lf3z+k2f/votPJaVpfAgjMqZ6BXLH9ATQHlgGTAHWA5nVEZQxxlHkUh54eyX5hS4AerRrxM0Xpp2kljG1WyCJ6VKcGR0GqmpWNcVjjCnltUWbWZFRss7S34afQVSkPRo2dVsg3+FJwLeWlIypGT/vP8Yzn5Y04d076BS6tE4KYUTG1IxAEtNWILqa4jDGlFLkUv44ayXH3U143domcfsA609k6odAEtMMoL+I2DBzY6rZ699sYem2wwBERQh/u7Yn0daEZ+qJQL7TpwBLgI9F5PRqiseYem/LgWye/nSDZ//uQZ05va014Zn6w2fnBxH50svhaKAPsFJEtgPbAZeXcqqqg4MTojH1R5FL+cPbP5JX4PxYdW2TxJ0DbOZwU7/465U3wM+5CJyZx1N9nLdxTMZUwssLfya9TBPeGTSIsiY8U7/4S0wDaywKYwyrdx7lH5+V9MIbN/gUurdrFMKIjAkNn4lJVRfWZCDG1Gd5BUXc9/9WUOhyGhvOat+YO6wXnqmnrI3AmDAwZd56Nu47BkB8g0j+cd2ZNpDW1Fv2nW9MiC3aeIDXv9nq2f/zlaeT2jwhdAEZE2KBTOLqrZeeN/nAAWAp8F9V3XOS8sbUW0dzCnjg7R89+4O7tOTXfVJCGJExoRfIXHkD3FvF9+zipc/9BviriNytqtMqF54xdZeqMuGD1ezJLFmRdvI1Z9hyFqbeCyQxDQSGAfcB3wP/BbbhjGNKxUlEfYF/4twtDQJGA6+IyDpV/TZoURtTB7y9dAezf9zl2X/yVz1o0TAmhBEZEx4CSUxFwD3Avao61cv550XkTuBZYJCq3iwii4DXgHGAJSZj3Dbty2LiB2s8+7/uk8Kl3VqHMCJjwkcgnR/+DKzxkZQAUNUXgdXABPf+68AW4IKqBGlMXZJXUMTdM5eTW1AEQOeWiUz8RbcQR2VM+AgkMfUB1py0lFOmT6n9tUCLQIIypi7765x1rN/jrB4TExXB1BFnEdcgMsRRGRM+AklMDXBWsD2Z9pRdHiMXOB5IUMbUVfNW7+E/323z7P/5ytNtjSVjygkkMa0EzheRS3wVEJGLcZrtVpY6nALsr1x4xtQdOw7n8MdZJV3DL+/emhvOrcjfesbUL4Ekpmfc5WeLyCsiMkBE0kQkVUT6i8jLwGx32b8DiEgj4CycXnzG1FsFRS7G/W8FmXmFALRrHMfkX1nXcGO8qXCvPFV9R0QmAI8Bt7hfpQlO1/GJqvqO+1hL4G/Ax0GI1Zhaa8rc9Z6F/yIjhOd+cyaN4m1BaGO8CaS7OKr6hIjMA+4G+gHt3Kd2AV8BL6hqeqnyG3F68xlTb328ajf/XrTFs3//JafSq0PTEEZkTHgLKDEBqOoy4KZqiMWYOufn/cf446ySR64Xd23JHf1t1nBj/AmbSVxF5DQRGSciM0RkvYi4RERF5NoK1B0hIl+LyFEROSYi6SJyl4j4/XyVrWdMReTkF3LHjKUcO+48V2rfNJ5nhp9JRIQ9VzLGn4DvmKrRHTgzRARERF4A7gTygC+AAmAwMBUYLCLDVbUoWPWMqQhV5aF3V/HTXmcpi5ioCF767dn2XMmYCvCZmERkGs6krH9S1b3u/YpSVb05wFhW43SUSMeZa+81oL+/CiJyDU5y2QP0cz/TQkRaAfOBq3Gehz0bjHrGVNSM77bxwYqSefAev6o73draarTGVISoqvcTIi6cxNRVVX9y71eUqmqVhrKLyAKcxDRcVWf5KJMO9AJGqeqb5c71BxbgJJ92quqqaj1vevfurenp6f6KmDpCVcnKyiIzM5OcnByKirzfUOcXuth/7DjFP1oJMZE0iW9Qg5EaU3OioqJo1KgRTZs2JSqq4o1wIrJUVXt7vaafemPc293l9sOCiCTjJJd84O3y51V1oYjsxOk52BdYXJV6pn5TVfbt20d2djZNmzaldevWREZGnjAOKb/QxaZ9x2jZxPl7Ji46kk4tEu25kqmTVJX8/HwOHjxIRkYGHTp0ICKi6o/ofSYmVX3D334YOMu9XaOquT7K/ICTYM6iJMFUtp6px7KyssjOzqZDhw5ERnpvDHC5lG0Hsyl0OUkpKkLo0CzekpKps0SEmJgY2rRpw44dOzh8+DDNmjWr8nVrc++zNPd2m58y28uVrUo9U49lZmbStGlTn0lJVdlxONczY7ggtG+aQIMom5zV1H0iQuPGjcnOzg7K9SqVmESkkYhcLCK/EZHzgxJJ4BLdW3//E8fc24ZBqOchImPdXcvT9++3aQDrg5ycHBITE32e33/sOEdy8z37bRvHkhgbTp1ejale8fHx5Ob6aoQKTECJyZ2QpgH7gE+AGZSamkhE7hSRXSLSNyjRnSQc99Z7743g1/NQ1VdVtbeq9m7Rwlb0qA+Kiop83i1l5haw52ieZ79pQgOaJlhnB1O/RERE4HIF0kfOz7UqWlBEEnB6q40GDgNzKfklX2we0Bq4KijR+Zfl3vr+M7bkXFapY5WtZ+o5bxOu5hUUkXEox7Of0CCKto3jbHJWU+8E83s+kDumB4CeOHdJHVX1yvIFVHUz8BMwKDjh+bXVve3gp0xKubJVqWdMGQVFLrYeyKbI3S+8QWQE7ZvFE2FJyZgqCSQxDceZrPVWVc3xU247JZO7Vqfl7m03EYnzUaZPubJVqWeMR3EPvPwip+kiQpweeNGRtbk/kTHhIZCfoo7AD6p6stVoDwBV7y94EqqaASzDWVl3ePnz7oGyyTgDZb+taj1jiqkqGYdzyMkv7oHnzIMX18A6OxgTDIEkpgIgtgLlkinp1VbdnnRvp4hI5+KDItISeNG9O9nL7A2VrWcMu4/mcTS3wLPftnEcSXE2B15qaioiwoIFC6r9vSZNmoSI+H1ddVXZR93F8W3durXM8dGjRyMiTJ8+PaixTZo0KSjXq48C+RNvA3CWiMSqap63AiLSBOc51LJAAxGRsylJCgCnu7dPiMgDxQdVtW+pf88SkZdwJoBdJSKfUzIZaxLwPs6krGVUtp4xB44d58CxkkaDFokxNEuMCWFE9VunTp248MILvZ47++yzaziaqpk+fTpjxoxh1KhRQUuStVUgiWkWMNn9+p2PMk/g9Gj7v0rEkgSc6+X4Kf4qqeqdIrIIuAtnbr1IYD0wDXjJ111PZeuZ+iszt4DdR0rGaTSKi6Z1o4o0IpjqcuGFF1b5l/iTTz7J+PHjadOmTXCCMlUWSGKaCowC7hGR3sC77uOpInIHzvOa/sAqnJnBA6KqCzix+3lF684EZtZUPVP/HC90sf1QjmfwW3yDKFKaxFu38DqgTZs2lpTCTIWfMbl74l0KfA+cj7NEBTjJaCowAKcJ7wpVzfd2DWNqo7W7Mjl47Diu4m7hURE2B14ASj/D2bRpEyNGjKBVq1bExMTQpUsXpkyZErSBmVWNr7yCggKmTJlC165diY2NpXXr1owcOZLt27dX6FnS3r17ue2220hOTiYmJoa0tDTGjx9PXl7ZpyGpqamMGePMk/3GG2+UeVY2evToIH7a2iGgbkSquhM4X0SGAENxeupFAhk4A27fV1/raBhTC207mM2o15fwxMCmAERFRJDWPMG6hVfCihUrGDduHM2bN2fgwIHs3buXRYsWMX78eHbs2MHzzz8f6hDLKCoqYtiwYcybN4+4uDgGDRpEw4YN+fLLL+nVqxdXXnnCUM4yMjIy6NWrF6rK+eefT2ZmJosWLWLKlCmsXbuWDz/80FP22muv5bvvvuObb7454bmZr2dodVml+req6jycWR6MqbP2ZeZx42tL2J/ldHaIjBDSmicQ42Vi1tTxc2o6vKDZOvmKGnmfZ599lokTJ/LII494lkb46quvGDhwIC+++CJ//OMfSUlJOclVas7zzz/PvHnz6NChA/PnzyctzZnT+fjx44wePfqkz7amTZvGLbfcwgsvvECDBs4UVevWreOcc85h9uzZfPPNN1xwwQUAPP3000yfPp1vvvkmKM/Najv7s88YL47mFjBy2hK2u6cbEoHUZgnENbDZwiurT58+TJw4scx6Pf369eOyyy7D5XIxf/78gK9Zvtmr9KuqnnvuOQD+8pe/eJISQExMDFOnTiUhIcFv/ZSUFJ577jlPUgLo2rUrN954IwBffPFFlWOsqyp1xyQikTiDaH12SVLV7b7OGRPOso8XcvP0H1i/x5kqMTJCaJbQgIQYG0BbFUOHDvWaMLp06cLcuXPZtWuXl1r++esuXhUZGRls2bKFyMhIrr/++hPON2vWjEsuuYT333/f5zUGDRpEXNyJk8t06dIFoFKft74I6CdNRM4FHgMuAvwN3tBAr21MOMjNL+LmN34gfdthz7Ep15xBbLT/+XxrqjmsNmvfvr3X40lJSQAndAioiOpq9tq5cyfg9NiLjvY+eLpDB3/TbVbP560vKpw8ROQC4HNKEtJhILM6gjImFPIKirj1zXS+23zIc+yRK0/n2l7JrFu3LoSR1Q3BWHK7pvlrEjzZ56mNnzdcBHJX8yhOUvoXMEFVbYU8U2ccLyzi9hlLWbTpgOfYQ5d34aYLbRHj+qht27aA09xWUFDg9a6p/NRGJngCSennAOtU9TZLSqYuKShycffM5cBX6nIAAB1SSURBVCzYUPJt/ftLTuW2/p1CGJUJpfbt29OhQweKiop4++23Tzh/6NAhPvvss6C+Z3EnicLCwqBetzYKJDEJsLK6AjEmFI4XFnH3zGV8tnav59i9gzpzz2C/M2GZeuCee+4B4OGHH2bbtm2e4/n5+dx7770cOxbcuarbtXNWC7Jm48Ca8lbhrE5rTJ2QV+A035W+U7qtf0fuu+TUEEZlwsW4ceP49NNP+fTTT+natSuDBg0iISGBxYsXk5uby8iRI3nzzTfLdAevir59+9K6dWuWLVtG79696datG9HR0VxwwQWeWSHqi0DumJ4FLhKRM6srGGNqSk5+ITdN/6FMUhrbryPjh3Sx+e8MAFFRUcyePZsnnniC9u3b89lnn7FgwQL69etHenq657lT8+bNg/J+MTExzJs3jyuuuIItW7YwY8YMXnvtNRYuXBiU69cmEsgMQiLyKHAn8Agwx8YqQe/evTU9PT3UYZgAZOUVcNP0H/hha0mX8HsHdea+S071mZTWrVtH165daypEE+YKCwvp3r07GzZsID09nV69eoU6pLAQyM+JiCxV1d7ezvlsyhORIj/XnApM9fOXpaqqjWMyYedoTgEjX1/CjxlHPMf+cNlp3DWws59apr5asWKFp0mtWE5ODg888AAbNmyge/fulpSqgb/kUZX2DGsLMWFnz9E8Rk77np/2ljy0nnBFV265qGMIozLh7O6772bNmjX07NmTNm3asH//fn788UcOHDhA48aN6/2cdtXFZ2JSVRsdZuqMTfuOMWraEnaWWujv8au6c2Nf/6P3Tf02duxYZs6cyapVq/j+++8BZw686667jj/84Q+kpqaGNsA6yprbTJ23fPthxkz/gSM5BQBERQhPD+/JVWe1C3FkJtyNHDmSkSNHhjqMescSk6nT5m/Yx50zlpFb4DwyjW8QyUu/7UX/U1uEODJjjC+WmEyd9b8l25nw/moKXU7P06YJDZg2ug9npjQOcWTGGH8sMZk6x+VSpsxbzytfbfYca9c4jv/cfA4dWySGMDJjTEVYYjJ1Sk5+Ib/73wo+LTXF0Oltknh9TB9aJflcPswYE0YsMZk6Y8/RPG558wdW7yxZjeXirq149tdn2iJ/xtQi9tNq6oTl2w9z+4yl7M087jl2y4VpPDS0K5ERNqzOmNrEEpOp9f67ZDsTP1hDfpELcJZCf/yX3RlxrvcVRI0x4c0Sk6m1jhcWMenDNfx3SYbnWFJsFC/ccDYXnWLdwY2prSwxmVpp99Fc7pixjBWl5rzr0rohr9zYiw7NEkIYmTGmqiwxmVpn4U/7uf//reBgdr7n2LCebZl8TQ/iG9i3tDG1nc2HZ2qNgiIXT85dx6hpSzxJKTJC+POVp/Psr8+0pBRiqampiIjnFRERQVJSEikpKVx88cU89NBDrFxpi2BXxoIFCxARBgwYEOpQaoQlJlMrZBzKYfjL3/LKwpJBsy0axjDj5nO5+cI0W9wvjFx22WWMGjWKkSNHMmTIEDp37kx6ejqTJ0+mZ8+eDBs2jD179gTlvSZNmoSIMGnSpKBcL1SKk/rWrVtDHUpYsD8xTdibs3I3499dSVZeoedYv1Nb8PfretI8MSaEkRlvxo8ff8Jf9i6Xi9mzZ3P//fcze/Zs+vfvz+LFi2nWrFlogqxlzjnnHNatW0d8fHyoQ6kRdsdkwtbRnALG/W85d81c5klKURHCQ5d3YfroPpaUapGIiAh++ctfkp6eTufOnfnpp5/4/e9/H+qwao34+Hi6dOlC+/b1YwiEJSYTluZv2Mel/1zIByt2eY4lN4nj7dvP47b+nYioL4Nms/bA65dD1t6Tl60FmjRpwj//+U8AZsyYcUKT3sGDB5kwYQI9evQgMTGRhIQEzj77bP7xj39QUFBQpqyI8OijjwLw6KOPlnm+Vb5pb9u2bdx555107NiRmJgYmjRpwsCBA5k5c6bXOEs3Ee7du5fbbruN5ORkYmJiSEtLY/z48eTl5Z1QLysri1dffZWrrrqKzp07Ex8fT2JiImeddRZ//etfyc3NLVN++vTpiAjbtm0DIC0trcznKG7aO9kzpjVr1jBy5EhSUlKIiYmhefPmDB06lLlz53otP3r0aESE6dOns2nTJkaMGEGrVq2IiYmhS5cuTJkyBZfL5bVuTbCmPBNWjh0v5K9z1vHfJdvLHL/m7GQmDjudpNhoHzXrqIVPwfbvYOEUuPLvoY4mKIYOHUrTpk05dOgQ8+fP5ze/+Q0Aq1atYsiQIezatYvk5GQGDBiAy+Xi+++/5/7772fOnDl8/PHHNGjQAIBRo0axYsUKfvzxR3r27MmZZ57peY/S//7+++8ZMmQIR44cIS0tjauvvpqDBw+ycOFCFixYwLx583jjjTe8PqfMyMigV69eqCrnn38+mZmZLFq0iClTprB27Vo+/PDDMuV//PFHbrvtNlq2bMlpp51G7969OXjwIN9//z0TJkzgww8/ZOHChcTGOvM2du7cmVGjRjFr1iyys7O55pprSEwsmWi49L99+fDDD7nuuus4fvw43bp146KLLmLHjh188sknzJ07lwkTJvD44497rbtixQrGjRtH8+bNGThwIHv37mXRokWMHz+eHTt28Pzzz5/0/auFqtqrCq9evXqpCY7P1uzR8574XDs8+JHn1evxT/WT1btDHZquXbu2Zt/w8RaqE5NOfD3eombjCECHDh0U0Pnz55+07ODBgxXQhx9+WFVVc3JyNC0tTQF94okntKCgwFP24MGDevHFFyugEydOLHOdiRMnej1eLDc3V1NSUhTQ3/3ud1pYWOg5t2rVKm3ZsqUC+vLLL3u9LqC33HKLHj9+3HNu7dq1mpiYqIAuWrSoTL2MjAz94osvtKioqMzxw4cP65AhQxTQyZMnnxBn8f/dli1bvH6O+fPnK6D9+/cvc3z37t2alJSkgD7zzDMn1ImPj1dA582bV+bcqFGjPJ9v4sSJZeJduHChRkREaEREhG7fvt1rPL4E8nMCpKuP36vWlGdCbm9mHnfMWMotb6az62hJ88jQHq359L7+XNqtdQijC5FxK6H7cIiKc/aj4qDHcBi3KrRxBUmLFs7MHAcPHgScJq0tW7Zw3XXX8dBDDxEVVdKY07RpU9544w2io6N54YUXcH6nVczbb79NRkYGHTp04KmnniIyMtJzrnv37p4mv6efftpr/ZSUFJ577jnPXRpA165dufHGGwH44osvypRPTk5m0KBBRESU/dXauHFjnnvuOQBmzZpV4fhP5l//+heZmZmcf/753H///WXODRgwgLvvvhvw/fn69OnDxIkTy8Tbr18/LrvsMlwuF/Pnzw9arIGwpjwTMkUu5a3vt/HUvA0cO17S465pQgMm/uJ0hvVsW3+7gTdsDTENoeg4RMU625gkaNgq1JEFRfHzi+JfiB9//DEAw4cP91q+bdu2nHLKKaxdu5aNGzdy6qmnVuh9Fi5cCMANN9xAdPSJzcBjxozhrrvuYtOmTezcuZN27dqVOT9o0CDi4uJOqNelSxcAdu3adcI5VeWbb77hq6++YseOHeTm5nruBAB++umnCsVeEcWfb/To0V7P33TTTTz11FMsWrSIoqKiMokZnGZVbz9jXbp0Ye7cuV4/X02wxGRC4oeth3hs9lpW7Txa5vjwXsn8aWhXmiQ08FGzHsneB73GQO8xkP46HKsbHSAADhw4ADh3QwCbNzvj03wlptL2799f4cS0c+dOwOlU4E1sbCxt27Zl586dXhOTr15wSUlJACd0gNi7dy+/+tWvWLx4sc+YMjMzfZ4L1Mk+X1paGhEREeTl5XHw4EFatmxZ5nygn6+mWGIyNWrnkVwmz13P7B/L/iXWsXkCf726B+d1snEtHr9+q+TfdaTjAzh3FMuXLwegR48eABQVFQFwxRVX0Lx5c7/1Axn7VHyX4u/O21/TYPkmuZO55ZZbWLx4MRdccAGTJk2iZ8+eNG7cmOjoaPLz84mJCe4Qh4p8Pn8C/Xw1xRKTqRE5+YW8snAzr3z1M3kFJd1QY6IiuL1/J+4Y0InY6Eg/VzB1xZw5czh8+DBRUVGe7s8pKSls2LCBO+64gyuuuCJo75WcnAyU3JGVl5eXx+7duwFOuFsKVHZ2Nh9//DGRkZF89NFHNG7cuMz5TZs2Ven63iQnJ7NhwwY2b97M4MGDTzi/detWXC4XsbGxnrvT2iA806WpM44XFjH9my30e2oBz36xsUxSuuKMNnzx+/7cd8mplpTqicOHD3PfffcBTnfv4qalyy+/HHA6KwSiuFNCYWGh1/P9+/cH4L///a/XMm+88QaqSufOnaucmI4ePYrL5aJhw4YnJCWAt956y0stx8k+hy/Fn+/NN9/0ev71118H4MILLyzToSTcWWIy1aKwyMXb6RkMenohk2av5cCxkpVlu7VN4v+N7csLI84muUn9mGKlvnO5XHz44Yf06dOHTZs20aVLF/72t795zo8dO5aUlBTeeOMNJk6cSE5OzgnXWL16tecXbbHiZLJu3Tqv7zt8+HBSUlLYsmULDz30UJlBo2vXrmXixIkAPPDAA1X+jK1ataJJkyYcOXLkhIG78+bN4+9/990ce7LP4cutt95Kw4YNWbRokafXX7GvvvrKMw6pts2yUXtSqKkVCotczF65i6lfbuLn/dllzrVpFMvvLj6Fa3ul2HLnddjkyZOZPn064DSV7d+/n2XLlnHkiLN21lVXXcUrr7xCkyZNPHUSExOZM2cOV155JY899hhTp07ljDPOoHXr1uzdu5ctW7awdetWzj33XMaMGeOpd9lllxEfH8+7775Lv3796NSpE5GRkQwbNoxhw4YRGxvL//3f/3H55Zfz9NNP895779GnTx8OHTrEggULyM/P58Ybb2Ts2LFV/tyRkZE8/PDDPPDAA9xwww1MnTqV1NRUfv75Z5YsWcKf/vQnnnjiCa91r776ahYsWMANN9zApZde6rnjmjJlit9naq1bt+Y///kP119/PePGjePf//433bt3Z9euXXz99de4XC4mTJjAkCFDqvz5apIlJhMUeQVFvL10B69+9TMZh8pOu9I0oQF3DezMDee2tya7euCTTz4BnAfyCQkJNG7cmF69enHOOecwYsQIunfv7rVejx49WLlyJS+++CIffPABy5YtIy8vjxYtWpCSksKNN97ItddeW6ZO69at+eijj3jsscdYvnw5ixYtQlVJTk5m2LBhAPTt25cVK1YwefJk5s2bx7vvvktcXBx9+/Zl7NixjBgxImjDEn7/+9+TmprK008/zZo1a1i9ejXdu3dnxowZ3HDDDT4T0913301mZiZvvfUWH330EcePOy0MEyZMOGlnj+I5CKdMmcKXX37JrFmzaNiwIZdeein33HMPQ4cODcpnq0kSyGA1c6LevXtrenp6qMMImaM5Bfz3h+28tmgL+7OOlzmXGBPF2H4duenCNBJjavffQOvWraNr166hDsOYsBbIz4mILFXV3t7O1e7fFiZkNuzJYvrirby/fCe5BUVlzjWJj2bMBWmMOi+VRvH1bG47Y0yVWWIyFZZf6OKLdXt549utfLf50AnnWyfFcmu/jvzmnBRbTdYYU2n228Oc1JpdR5m1dAcfrNjFIfeS5qV1bZPEmPNT+eVZbYmJsmdIxpiqscRkvNp1JJe5q/fwztIdrN194hQqkRHCkG6tGXV+Kn1Sm9TfOe2MMUFnicl47DySy9xVu5mzajfLtx/xWqZd4ziuObsdvzm3PW0anTi5pTHGVFW9T0wiMgK4AzgDiATWA68DL6lq6JZwrAFFLmVFxmEWbNjP/A37WL3T++SSMVERDOnemuG9Uji/U7P6s3qsMSYk6nViEpEXgDuBPOALoAAYDEwFBovIcFUt8nOJWkVV2XwgmyVbDrH454N89dN+juYWeC0bGSGc36kZQ3u0YWiPNjSKs951qmpNlsb4EMyhR/U2MYnINThJaQ/QT1U3uo+3AuYDVwN3A8+GLMgqyisoYt3uTJZvP8IPWw/xw9ZDHDh2YueFYlERwnmdmnFFjzZc2q01TW3pCY/IyEiKiopq1XxjxtQkl8sVtNnK6/NP2UPu7YPFSQlAVfeKyB3AAmC8iDwf7k16qsqezDy27M9m/Z4sVu86ypqdmWzaf4wil/+/YlolxTDg1JYMOK0FF5zSnKRYuzPyJj4+nmPHjnmdnNMYAzk5OV4XVayMepmYRCQZ6AXkAydMZ6yqC0VkJ9AO6Av4XvWrBhS5lAPHjrPnaB57MvPYm5nHnqN5bDuUw5b92Ww5kH3CIFdfGsVF0ye1KeekNeHCzi3o2qahNU9VQFJSEgcOHKBhw4YnrAJqTH2nqhw5coSEhISgXK9eJibgLPd2jarm+ijzA05iOosgJ6ajuQX847OfKHS5KChUClwuCouUQpeL/EIXWXmFZOYVkplbQGZeAceOF1KZ5lsRSGueQPe2jeiT2oRz0ppxSstE67xQCQ0bNiQ3N5dt27bRtGlTEhMTiYyMtKRu6jVVJT8/n4MHD1JYWFhmYt6qqK+JqXgd4m1+ymwvVzZojhcWMX3x1qBes1FcNB1bJNCpRSLd2ybRrV0jurZJqvVz1IULEaFly5ZkZWWRmZnJvn37PKuuGlOfRUVF0ahRI1q2bGnPmKoo0b3N9lPmmHvbsPwJERkLjAVo3759wG8eXYkvXtOEBrRKiqV1UgytG8XSKimWto3j6NQigY7NE2liHRWqnYiQlJREUlJSqEMxpk6rr4mpuP2lUv0bVfVV4FVwZhcPtH58TCSPXHk60ZFCVGQE0ZERzr8jnG1ibBRJsdE0iosmKTaaxNgoW7/IGFNv1NfElOXeJvopU3wuy0+ZSomJiuSmC4PeQmiMMXVCfV1afat728FPmZRyZY0xxtSA+pqYlru33UTEV8f7PuXKGmOMqQH1MjGpagawDGgADC9/XkT6A8k4s0J8W7PRGWNM/VYvE5Pbk+7tFBHpXHxQRFoCL7p3J4f7rA/GGFPX1NfOD6jqLBF5CWdm8VUi8jklk7gmAe/jTOZqjDGmBtXbxASgqneKyCLgLqA/JcteTKMeLHthjDHhqF4nJgBVnQnMDHUcxhhjHBLMNTTqIxHZj/+pjfxpDhwIYjgm/NjXuG6zr2/ldVDVFt5OWGIKIRFJV9XeoY7DVB/7Gtdt9vWtHvW5V54xxpgwZInJGGNMWLHEFFqvhjoAU+3sa1y32de3GtgzJmOMMWHF7piMMcaEFUtMxhhjwoolphAQkREi8rWIHBWRYyKSLiJ3iYh9PWoxETlNRMaJyAwRWS8iLhFREbk21LGZqhORaBEZLCLPiMh3IrJbRPJFZKeIzBKRAaGOsa6wZ0w1TEReAO4E8oAvKJmfryHwHjBcVYtCF6GpLBH5JzDOy6nhqjqrpuMxwSUiFwOfuXf3AEuBbOB0oLv7+OOq+kgIwqtT7C/0GiQi1+AkpT3AGap6papeDZwCrAOuBu4OYYimalYDfwOuBzoDC0MbjgkyF/AO0E9V27h/fq9X1R7Ar4Ei4M8iMjCkUdYBdsdUg0QkHegFjFLVN8ud6w8swEla7WwC2dpPRBbgTA5sd0z1gIj8G7gZmKaqN4c6ntrM7phqiIgk4ySlfODt8udVdSGwE2gN9K3Z6IwxQVC82nVySKOoAywx1Zyz3Ns1qprro8wP5coaY2qPU9zb3SGNog6wxFRz0txbfzORby9X1hhTC4hIa2C0e/edEIZSJ1hiqjmJ7m22nzLH3NuG1RyLMSZIRCQKmAE0Ar5Q1dkhDqnWs8RUc8S9td4mxtQtL+MM+cgAfhviWOoES0w1J8u9TfRTpvhclp8yxpgwISLP4vTE2wMMVtU9IQ6pTrDEVHO2urcd/JRJKVfWGBOmROQZ4F5gP05S2hjikOoMS0w1p7graTcRifNRpk+5ssaYMCQiTwH3AweBS1R1bYhDqlMsMdUQVc0AlgENgOHlz7sH2CbjNAl8W7PRGWMqSkQmA38ADuMkpR9DHFKdY4mpZj3p3k4Rkc7FB0WkJfCie3eyzfpgTHgSkceBB4EjOEnJWjeqgU1JVMNE5EXgDpxJXD+nZBLXJOB94FqbxLV2EpGzKfkDA5zJPRsCG4FDxQdV1Wb2qIVEZBjwgXs3HVjjo+h6VZ1cM1HVTZaYQkBERgB3AT2ASGA9MA14ye6Wai/3sgfzT1ZOVeVkZUz4EZHRwOsVKLpQVQdUbzR1myUmY4wxYcWeMRljjAkrlpiMMcaEFUtMxhhjwoolJmOMMWHFEpMxxpiwYonJGGNMWLHEZIwxJqxYYjKmmojIABFREVkQ6liqSkQedH+WIVW4xtki4hKRp4MZm6l7LDEZU0kistX9yzo11LFUJxFpAzwMfKWq8yp7HVVdBrwL3CsipwQrPlP3WGIypvosAboCI0MdSBU9ijPn36NBulY0JRMaG3MCm5LImEoSka04Cz+mqerW0EZTPUSkGbAD2AV01iD8whCRH4CzgI6qur2q1zN1j90xGRMgERktIkrJasRb3E16Wrppz9czJhFJdR/fKiIRInK/iKwRkVwR2SEifxeReHfZJiLyT3fZ4yKyUUTu9xObiMivReRTETngrrNdRP5VySbHm4BY4E1vSUlEGovIE+74c0p9hgUi8pCPa76BM3nxbZWIx9QDUaEOwJhaaBPOL9drgQTgHeBYqfPHvFXyYSZwJbDAfd1+wH1AVxG5AfgOpxltEdDUff4ZEYlV1SdKX0hEooH/Ab8CcnGWZtgLdAduAa4RkUtVNT2A+K5ybz8vf8KdPL/BWd5jn7tMNtDGfawv3pvsiq/1S5xnV8aUpar2spe9KvECtgIKpPo4P8B9fkG546nu44qz5EnbUudSgAPuc6uAt4HYUuevcJ/LBOLLXXey+9xCILncubvd5zYBURX8fPFAvvsV6+X8SPc1Pyp/TZw7okE+ris461Mp0CrUX0d7hd/LmvKMCa17VXVX8Y6qZgAz3LsdgDtUNa/U+TnASpy7qN7Fx0WkKXAvzt3acFXdUfpNVHUqMAfoBFxewdi64XRU2FI6hlJaubefq2phufcrUtUvvV1UVRVY5949s4KxmHrEEpMxoVMAePvlvcm9TVfVA17Ob3Rv25Y6NhCIw1mkbp+P91vo3p5XwfhaurcHfZxf4t4+KCK/FZHGFbwulKzo28pvKVMv2TMmY0JnT/k7DbfiZ1Q7vJwrfT621LGO7u0V7o4Z/rSoYHyN3NtMbydVdaGIPAU8APwHUBFZj/M87B1V/cTPtYuvGUgyM/WEJSZjQsdVxfOlRbq3G3A6TPjzfQWvecS9TfJVQFUfFJGXcToyXAhcANwK3CoinwJX+Ei+xdc8XMFYTD1iicmYuiHDvV2lqqODdM3iJsFm/gqp6hbgn+4XInIh8F/gUpzu5q96qVZ8TV/NjqYes2dMxlRevnsbDn/gfY7zzOriAJ/1+LMGOA6kiUhcRSup6iJgunu3Z/nzIiJAF/fu8irGaOogS0zGVN5O97ZrSKMAVHUv8ALOM5sPRaRL+TLuwbq3iEiFOhyoai5Os1800MvL9a4WkX4iElHueBxwsXt3m5dLdwGaAGv8dNQw9Vg4/KVnTG31Hs5Ypbfcz1OKn8k8qKq+erJVpz/i9NS7DlgtIiuALTidJFJwEmgD93ZvBa/5Ps6g3otxOjWU1h8YB+wXkeXAfpwOE+fjDAZeD7zi5ZrFSeuDCsZg6hlLTMZU3lSch/g34MzeEOM+/hd8d7GuNqpaAFwvIm/hPNs5BzgDyAJ248wy8QHwcwCXnQ78FRgpIo+6xyCVPpeH0+mhO9AcJzlvwnnG9JqqZnm55iigCO9JyxibxNUY45+7191twGBfg2YDuFYPnAHC76jqtcGIz9Q9lpiMMX6JSGvgJ2C5qvav4rVmAcOAbqq68WTlTf1knR+MMX6p6h6c5sl+VV3BFmeC2ectKRl/7I7JGGNMWLE7JmOMMWHFEpMxxpiwYonJGGNMWLHEZIwxJqxYYjLGGBNWLDEZY4wJK/8flzWio0maZNIAAAAASUVORK5CYII=\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "yf=300\n", | |
| "dmdt=root[0]\n", | |
| "t2=(m0-mf)/dmdt\n", | |
| "t=np.linspace(0,t2,N)\n", | |
| "N2=int(t2/dt)\n", | |
| "numsol_h=np.zeros([N2,3])\n", | |
| "numsol_h[0,0]=y0\n", | |
| "numsol_h[0,1]=v0\n", | |
| "numsol_h[0,2]=m0\n", | |
| "for i in range(N2-1):\n", | |
| " numsol_h[i+1] = heun_step(numsol_h[i],lambda state:rocket(state,dmdt=dmdt),dt)\n", | |
| " heun_h=numsol_h[:,0]\n", | |
| "plt.plot(t[:N2],heun_h, label='In Flight')\n", | |
| "plt.plot(t[-1], yf, '*', label='Detonation')\n", | |
| "plt.xlabel('time (s)')\n", | |
| "plt.ylabel('height (m)')\n", | |
| "plt.legend(loc='best')" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## References\n", | |
| "\n", | |
| "1. Math 24 _Rocket Motion_. <https://www.math24.net/rocket-motion/\\>\n", | |
| "\n", | |
| "2. Kasdin and Paley. _Engineering Dynamics_. [ch 6-Linear Momentum of a Multiparticle System pp234-235](https://www.jstor.org/stable/j.ctvcm4ggj.9) Princeton University Press \n", | |
| "\n", | |
| "3. <https://en.wikipedia.org/wiki/Specific_impulse>\n", | |
| "\n", | |
| "4. <https://www.apogeerockets.com/Rocket_Motors/Estes_Motors/13mm_Motors/Estes_13mm_1_4A3-3T>" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": null, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [] | |
| } | |
| ], | |
| "metadata": { | |
| "kernelspec": { | |
| "display_name": "Python 3", | |
| "language": "python", | |
| "name": "python3" | |
| }, | |
| "language_info": { | |
| "codemirror_mode": { | |
| "name": "ipython", | |
| "version": 3 | |
| }, | |
| "file_extension": ".py", | |
| "mimetype": "text/x-python", | |
| "name": "python", | |
| "nbconvert_exporter": "python", | |
| "pygments_lexer": "ipython3", | |
| "version": "3.7.3" | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 4 | |
| } |