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.ipynb
Go to fileThis commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
664 lines (664 sloc)
114 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", | |
"![Initial condition of firework with FBD and sum of momentum](../images/firework.png)" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"We are going to end this module with a __bang__ by looking at the flight path of a firework. Shown above is the initial condition of a firework, the _Freedom Flyer_ in (a), its final height where it detonates in (b), the applied forces in the __Free Body Diagram (FBD)__ in (c), and the __momentum__ of the firework $m\\mathbf{v}$ and the propellent $dm \\mathbf{u}$ in (d). \n", | |
"\n", | |
"The resulting equation of motion is that the acceleration is proportional to the speed of the propellent and the mass rate change $\\frac{dm}{dt}$ as such\n", | |
"\n", | |
"$$\\begin{equation}\n", | |
"m\\frac{dv}{dt} = u\\frac{dm}{dt} -mg - cv^2.~~~~~~~~(1)\n", | |
"\\end{equation}$$\n", | |
"\n", | |
"If we assume that the acceleration and the propellent momentum are much greater than the forces of gravity and drag, then the equation is simplified to the conservation of momentum. A further simplification is that the speed of the propellant is constant, $u=constant$, then the equation can be integrated to obtain an analytical rocket equation solution of [Tsiolkovsky](https://www.math24.net/rocket-motion/) [1,2], \n", | |
"\n", | |
"$$\\begin{equation}\n", | |
"m\\frac{dv}{dt} = u\\frac{dm}{dt}~~~~~(2.a)\n", | |
"\\end{equation}$$\n", | |
"\n", | |
"$$\\begin{equation}\n", | |
"\\frac{m_{f}}{m_{0}}=e^{-\\Delta v / u},~~~~~(2.b) \n", | |
"\\end{equation}$$\n", | |
"\n", | |
"where $m_f$ and $m_0$ are the mass at beginning and end of flight, $u$ is the speed of the propellent, and $\\Delta v=v_{final}-v_{initial}$ is the change in speed of the rocket from beginning to end of flight. Equation 2.b only relates the final velocity to the change in mass and propellent speed. When you integrate Eqn 2.a, you will have to compare the velocity as a function of mass loss. \n", | |
"\n", | |
"Your first objective is to integrate a numerical model that converges to equation (2.b), the Tsiolkovsky equation. Next, you will add drag and gravity and compare the results _between equations (1) and (2)_. Finally, you will vary the mass change rate to achieve the desired detonation height. \n" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"1. Create a `simplerocket` function that returns the velocity, $v$, the acceleration, $a$, and the mass rate change $\\frac{dm}{dt}$, as a function of the $state = [position,~velocity,~mass] = [y,~v,~m]$ using eqn (2.a). Where the mass rate change $\\frac{dm}{dt}$ and the propellent speed $u$ are constants. The average velocity of gun powder propellent used in firework rockets is $u=250$ m/s [3,4]. \n", | |
"\n", | |
"$\\frac{d~state}{dt} = f(state)$\n", | |
"\n", | |
"$\\left[\\begin{array}{c} v\\\\a\\\\ \\frac{dm}{dt} \\end{array}\\right] = \\left[\\begin{array}{c} v\\\\ \\frac{u}{m}\\frac{dm}{dt} \\\\ \\frac{dm}{dt} \\end{array}\\right]$\n", | |
"\n", | |
"Use [two integration methods](../notebooks/03_Get_Oscillations.ipynb) to integrate the `simplerocket` function, one explicit method and one implicit method. Demonstrate that the solutions converge to equation (2.b) the Tsiolkovsky equation. Use an initial state of y=0 m, v=0 m/s, and m=0.25 kg. \n", | |
"\n", | |
"Integrate the function until mass, $m_{f}=0.05~kg$, using a mass rate change of $\\frac{dm}{dt}=0.05$ kg/s. \n", | |
"\n", | |
"_Hint: your integrated solution will have a current mass that you can use to create $\\frac{m_{f}}{m_{0}}$ by dividing state[2]/(initial mass), then your plot of velocity(t) vs mass(t)/mass(0) should match Tsiolkovsky's_" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 1, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"import numpy as np\n", | |
"import matplotlib.pyplot as plt\n", | |
"%matplotlib inline\n", | |
"plt.rcParams.update({'font.size': 22})\n", | |
"plt.rcParams['lines.linewidth'] = 3" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 2, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"def simplerocket(state,dmdt=0.05, u=250):\n", | |
" '''Computes the right-hand side of the differential equation\n", | |
" for the acceleration of a rocket, without drag or gravity, in SI units.\n", | |
" \n", | |
" Arguments\n", | |
" ---------- \n", | |
" state : array of three dependent variables [y v m]^T\n", | |
" dmdt : mass rate change of rocket in kilograms/s default set to 0.05 kg/s\n", | |
" u : speed of propellent expelled (default is 250 m/s)\n", | |
" \n", | |
" Returns\n", | |
" -------\n", | |
" derivs: array of three derivatives [v (u/m*dmdt) -dmdt]^T\n", | |
" '''\n", | |
" dstate = np.array([state[1], (u*dmdt)/state[2], -dmdt])\n", | |
" return dstate" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 3, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"def rk2_step(state, rhs, dt):\n", | |
" '''Update a state to the next time increment using modified Euler's method.\n", | |
" \n", | |
" Arguments\n", | |
" ---------\n", | |
" state : array of dependent variables\n", | |
" rhs : function that computes the RHS of the DiffEq\n", | |
" dt : float, time increment\n", | |
" \n", | |
" Returns\n", | |
" -------\n", | |
" next_state : array, updated after one time increment'''\n", | |
" \n", | |
" mid_state = state + rhs(state) * dt*0.5 \n", | |
" next_state = state + rhs(mid_state)*dt\n", | |
" \n", | |
" return next_state" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 4, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"def heun_step(state,rhs,dt,etol=0.000001,maxiters = 100):\n", | |
" '''Update a state to the next time increment using the implicit Heun's method.\n", | |
" \n", | |
" Arguments\n", | |
" ---------\n", | |
" state : array of dependent variables\n", | |
" rhs : function that computes the RHS of the DiffEq\n", | |
" dt : float, time increment\n", | |
" etol : tolerance in error for each time step corrector\n", | |
" maxiters: maximum number of iterations each time step can take\n", | |
" \n", | |
" Returns\n", | |
" -------\n", | |
" next_state : array, updated after one time increment'''\n", | |
" e=1\n", | |
" eps=np.finfo('float64').eps\n", | |
" next_state = state + rhs(state)*dt\n", | |
" ################### New iterative correction #########################\n", | |
" for n in range(0,maxiters):\n", | |
" next_state_old = next_state\n", | |
" next_state = state + (rhs(state)+rhs(next_state))/2*dt\n", | |
" e=np.sum(np.abs(next_state-next_state_old)/np.abs(next_state+eps))\n", | |
" if e<etol:\n", | |
" break\n", | |
" ############### end of iterative correction #########################\n", | |
" return next_state" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 75, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"y0 = 0 # initial position\n", | |
"v0 = 0 # initial velocity\n", | |
"m0 = 0.25 # initial mass\n", | |
"dmdt = 0.05\n", | |
"time = (0.25-0.05)/dmdt\n", | |
"t= np.linspace(0,time,10000)\n", | |
"dt=t[1]-t[0]\n", | |
"N=int(time/dt)\n", | |
"mf = np.linspace(0.25, 0.05, N)\n", | |
"\n", | |
"#initialize solution array\n", | |
"num_heun = np.zeros([N,3])\n", | |
"num_rk2 = np.zeros([N,3])\n", | |
"\n", | |
"#Set intial conditions\n", | |
"num_heun[0,0] = y0\n", | |
"num_heun[0,1] = v0\n", | |
"num_heun[0,2] = m0\n", | |
"\n", | |
"num_rk2[0,0] = y0\n", | |
"num_rk2[0,1] = v0\n", | |
"num_rk2[0,2] = m0\n", | |
"\n", | |
"d_m = mf/m0\n", | |
"vf = -250*np.log(d_m)\n", | |
"for i in range(N-1):\n", | |
" num_heun[i+1] = heun_step(num_heun[i], simplerocket, dt)\n", | |
"for i in range(N-1):\n", | |
" num_rk2[i+1] = rk2_step(num_rk2[i], simplerocket, dt)\n", | |
" \n", | |
"#print(num_rk2[:,1])" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 101, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAtwAAAEaCAYAAAAi8IBdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd1hUR/cH8O/ZwiId6aA0USwoAgYb2IJdUyRqbNGoUaMmlsSWRBOjSUzUvIktxhZLNHZjLD+7Yk8ES+wdolJEOtK2zO+PBUQEpCwu5XyeZ599996Ze8+CLzk7e2aGhBBgjDHGGGOMlQ+JvgNgjDHGGGOsKuOEmzHGGGOMsXLECTdjjDHGGGPliBNuxhhjjDHGyhEn3IwxxhhjjJUjmb4DqOysra2Fq6urvsNgjLFKJSws7IkQwkbfcTDG2KvACXcZubq6IjQ0VN9hMMZYpUJEEfqOgTHGXhUuKWGMMcYYY6wcccLNGGOMMcZYOaqwCTcRfUtEIvvxaRHt+hPRCSJKIqJUIgolojFEVOR7K20/xhhjjDHGSqJCJpdE9BqAyQCK3HeeiBYDWA+gGYATAA4CqAdgEYCtRCTVZT/GGGOMMcZKqsIl3ESkALAaQAyAnUW0CwYwGkA0gCZCiB5CiLcB1AVwHcDbAMbqqh9jjDHGGGOlUeESbgBfA2gIYBSApCLaTct+niKEuJ1zUAgRA+DD7JdTCygRKW0/xhhjjDHGSqxCJZVE1BzAJwA2CCF2FdGuFgA/AFkAtuQ/L4QIAfAIgD2AFmXtVx4y0p/izMpP8Oje1fK8DWOMMcYY07MKk3ATkSGANQDiAYx7SXOf7OerQoj0Qtqcy9e2LP106sa5Q3jygx/8H6zA3r9GQKNWl9etGGOMMcaYnlWYhBvANwA8AXwkhHjykrZu2c9FbZzwX762ZemnU1KZARJkcejvaIcFlsn4ZefU8roVY4wxxhjTswqx0yQRtQIwHsCfQohNxehikv38tIg2qdnPpjro9xwiGgFgBAA4OzsXcamC1fVpg3Gh3rhmEAMA2JSwF29EjUBth7olvhZjjLHChYWFWUokksFSqXSAEMIaAOk7JsZY1UBEqUKIEyqVaguAo35+foWurqf3hJuIagD4DUAytKuHFKtb9nORywbqsN9zhBDLACwDgGbNmpXqWtPe/g2XdnZBnEyCBJkE3+waiqUjTpQlLMYYY3mEhYUZyGSy1TVr1mxsbW2damhoGEfE+TZjrOyEEFCpVNLk5OSg2NjYjunp6SvCwsK+LyzprgglJd9Cuwb2RCFEVDH7pGQ/mxTRJudcSp5jpe2nc/bWtTHQ5p3c16HyeBw98lt53pIxxqqbfmZmZo1r1aoVX6NGjSxOthljukJEkMvlaisrqyQPD4/kGjVqDAfQvrD2FSHhfhuABsBgIjqW9wGgS3abD7OPrch+HZ797FLEdWvna1uWfuVi+Bsz4Z9hDN+MDGx9FA3PE9/iaUpied+WMcaqBblc3q1mzZqcaDPGypVMJtPY2NhAJpP1LrTNqwyoCBIAbYs47579sMh+fSH7uRER1ShkxZHX8rUtS79y81X3tTBe0R41oQLwGGfXTUKL0ctfxa0ZY6yqa2JsbJym7yAYY1WfmZlZKhEFFnZe7yPcQghXIQQV9IB2mUAAmJR9rGl2nwcAzgMwAPDCpwkiagugFrS7SZ7Jc69S9StPtWvVw12fz3Nf+8dswY1zh17FrRljrEoTQtSQSqUafcfBGKv6ZDKZSghRaMmy3hPuMvgu+/l7IvLIOUhEtgCWZL+cI4TI/8e2tP3KTbOeo/CvoXZgXUIC8QfHI+Upl5YwxlhZcTkJY+xVeNnfmkqbcAshtgL4BdpdIS8T0S4i2g7gNrRbw/8JYJGu+pUnkkhg028xkoQhVpqb4iMHYNamQa8yBMYYY4wxVk4qbcINAEKI0QAGQFsm0hZAZwB3AIwFECyEKHALx9L2K08OLp5Y7fkWfqppiSwJ4YDkPnafXP2qw2CMMcYYYzpWoRNuIcSQ7NrteUW02SCEaC2EMBNCGAsh/IQQi19WElLafuVpTJ9FaJApBQCoibDoxjwkJMXqKxzGGGOMMaYDFTrhrm5kMjmmtlkII402538kJ3y1uZ+eo2KMMVZVOTk5NSYiv927dxe6uzIA+Pv7exKR34IFC6xeVWyMVSWccFcwvvUD0afGsxUSjxjEYOOBH/UYEWOMMcYYKwtOuCugCb0XommGIvf1sgcr8fBxuP4CYowxxhhjpcYJdwUkkUoxvcsKmKu1pSWxMglm/TlQz1ExxhhjjLHS4IS7gqrn0hSDar6Z+/q0IgnL/vy8iB6MMcbYq3fkyBHjHj16uNvZ2TWRy+W+lpaW3h06dPDYv3//C5uA3Lx504CI/JycnBoXdj0i8iMiv6KOL1++3LJp06b1jYyMfIyNjX1atmxZr6D7MVZRcMJdgY1861u0yHw2j+X3+D8REX5VjxExxhhjz3z55Zd2QUFB9ffu3WtpY2OjDAoKSnRxcckMCQkx79atm+f8+fOtdX3P8ePHO44aNcpdLpeL9u3bJ9nZ2WWdPXvWtGfPnvUOHTpkrOv7MaYLMn0HwIr25VsbMHhXdyiEGj88jsPT9SOgnhICqYx/dYwxVhauU/e8MIpaWYTP6R6m7xi2bt1q9vXXX9eysbFRbty48W6HDh2e5pw7cOCAcXBwcN2pU6c6d+zYMaVJkyaZurrv6tWrbY8dO3Y9MDAwDQDUajUGDhzosnHjRusZM2Y4BgUF3dbVvRjTFR7hruBq2bpiar2p2PDwMbyystBQeQX/rP1M32ExxhirQnr27Fkvp2SjoMe5c+deKNf4+uuvHQFg0aJF4XmTbQDo1KnT0wkTJkSpVCpauHChjS5jnTx58qOcZBsApFIp5s2b9wgAwsLCTDMzM4veY5sxPeBh0kqgY+AgnLlzBS0jlgIA/COW4frfr6NB8856jowxxlhVEBAQkGxra6ss7HxISIh5XFxcbs4QFRUlu3LlirGJiYm6V69eyQX1ef3111NmzZqF0NBQndZWBwcHJ+U/5uTkpDIzM1MnJydLY2JipM7Ozipd3pOxsuKEu5LwH/QNrv1wCg2zLkNKAhkHx+CBywHUtnfXd2iMMVYpVYSyjIpiypQp0T169Egp7Ly/v79nXFxcbuJ869YtAyEEUlNTpXK5vMjSnPj4eJ3mGh4eHlkFHTcxMVEnJydL09PT+dt7VuFwwl1JSGUyWL23Bokr2uCssQYzrRWov7MfVg4/DYlUqu/wGGOMVSMqlYoAbZLbqVOnxKLaWllZFXu0Wa1Wv7SNlP+bxyohTrgrEbtadbDjtYmY8WQ1ACDUMA0/bBqBqf1X6jcwxhhj1Yq7u3sWAMhkMrFt27bw4vZTKBQCANLS0gochb59+7aBTgJkrILhr10qmbe7f4J2Wc/mn2zO+ht7T63TY0SMMcaqGzc3N2XdunXTExMTZbt37zZ9eQ8tBwcHlVwuF4mJibLIyMgXBv127NhhrttIGasYOOGuhGa9uwVu2RVsSiLMvzEH/0Xd0m9QjDHGqpUZM2ZEAsCwYcPctm/fbpb/fEZGBq1fv94879rYCoVCNGvWLBUAJk2a5KjRaHLb79+/3+T77793egWhM/bKccJdCVmYWmF6q59gkr31+2OZBJ/t6geVqtAJ5owxxphODRw4MPHLL798GBcXJw8ODq7r6urq1aFDB48uXbq4N2nSpL6NjY33wIEDPc6fP2+Ut9/MmTMfyeVysWHDBpu6des26tq1q3vjxo0bdOvWzXPw4MGP9fV+GCtPJUq4iciCiN4ioplEtJSINhLRL9mv3yQii/IKlD3vtUavY5j1O7mvLymyMPP3fnqMiDHGWHXz1VdfxZw4ceJanz59nmg0Gpw+fdrsxIkT5snJyTJ/f/+U+fPnRwwePDg+b5+OHTs+3bVr162WLVumREdHGxw7dswcABYtWnT/559/jtTPO2GsfJEQougGRFIAvQCMBhAIIGdB+bwLy4s8z8cBLAGwQwjx8unGlVyzZs1EaGio3u7/ycquOCB7CACQCIFpjkPxbqeJeouHMcaKg4jChBDNyvMely5dCvf29n5SnvdgjLEcly5dsvb29nYt6FyRq5QQUT8AcwDUgjbBfgLgLIBrAOIBJAMwA2AFoCGAFgDaAWgL4AERTRVCbNTJu2AF+mbgNjxY0wrXFWpoiLD44Uo0vtsGjeqU63/HGGOMMcZYMRWacBPRKWgT6FgAPwNYI4S49LILElFTAEMA9AOwnog+EkK01k24LD9DhRFmdfwNHxwZiASZBEkSwsmdI+Ex+gQUhkYvvwBjjDHGGCtXRdVw1wEwEYCzEGJicZJtABBCXBRCjAdQG8An2ddh5cjTzQcfuYyCmVqDBTFPMDL5Di4tHQaRZ/Y3Y4wxxhjTjyITbiHEz0KIArdQfRkhRJYQ4icAvPf4K9A76CPMVgSjXXo6AMA/cS/+3jRHz1ExxhhjjLFCE24hxFNd3EAIkaaL67CXa9fvK5wz75L7utmNubhyapceI2KMMcYYY7wOdxVCEgkaj1qFW7J6AIAkKbDk0qc4f/24niNjjDHGGKu+ip1wE1FNIvIlopr5jjsQ0WoiukBEO4jIW/dhsuIyrGEMi/c346xBTfRzsscJYwN8eXIM4hKj9R0aY4wxxli1VJIR7mkAzkE7GRIAQEQGAE4CGATAG8CbAI4SEW/Nqke2Tm6IaTUNsVIpACDcAJi0+U3eiZIxxhhjTA9KknC3B3A/32olfQG4AQgB0AXAYgAWAMbqLEJWKm+2G4GBhoG5r88p0vDZ2rf1GBFjjDHGWPVUkoS7FoA7+Y71gHZ3yeFCiANCiI8A3AfQVUfxsTL45N1f8HqWXe7r/5NG4MdNY/QYEWOMMcZY9VOShNsS2p0m82oJ4JYQ4l6eYxeQp+yE6dec9/6Cd4Yi9/W69BBsOviTHiNijDHGGKteSpJwp0O7hTsAgIhqQzvqfSpfu0wACrAKwVBhhPnv7IRb9mrqKiL89GA5Tl3aq9/AGGOMMcaqiZIk3DcABORZpaQ/tOUk+decqwUgRgexMR2xs3LC7HbLYaXS7jyZKpXgq3OTcf/RDT1HxhhjjDFW9ZUk4V4HwBjAP0S0GcDXAFIB7MxpQEQKAL4AbuoySFZ2Teq2wNQGn8FQIwAA0XLCkh39kZGuk/2NGGOMMcZYIUqScP8CYAO0W7W/AyALwAdCiKQ8bXpCm5SH6CxCpjNdWg3ASKteICHQNfUpZsfdxbXF/aBRq/UdGmOMsVfMycmpMRH55X0oFApfBweHxt26dXPfs2ePib5jrGiuXLmiICI/mUzmV1ibJUuW1JTJZL4SicTv66+/tn2V8bGKq9gJtxBCI4QYCKAOgFYAnIQQm/M1uwegN4A1uguR6dLwN77GZFl7zImNg0IAvqkhOLd0BIRGo+/QGGOM6UFAQEByr1694nr16hUXEBCQBAD/93//Z9mjRw/PmTNncsJYAnPmzLEZO3asGwD63//+Fz5jxozHZb2mt7d3fSLyO3z4sHFB5y9cuGBIRH5ubm6NynovVn4KTbiJ6F0iMs1/XAhxXwhxVgiRXMC580KIbUII3tawAhvQ/2f8Y9M793Xz2K04u266HiNijDGmL1OmTInetm1b+LZt28IPHz58Nzw8/MqAAQNiAeCbb76pdffuXbm+Y6wMpk2bZj9t2jRnmUwmVqxYcXfcuHFx+o6JVRxFjXBvAPCYiHYT0TAisnlVQbHyRRIJ/Ef9ijCTdgC0M1+vxK/Fj5t4vyLGGKvuFAqFWLp06QNjY2ONUqmkXbt2mek7poruww8/dJozZ46ToaGhZuPGjXeGDBmSqO+YWMVSVMI9BcBFaDexWQYgkoiOEtHHROT8SqJj5UYilcJr7B+4bNAE31pZYkFNC6xNP4Y1e77Rd2iMMcb0zMTERLi6umYAQExMzAsj3HZ2dk2IyK+w0W8/Pz9PIvLbv3+/SWHHQ0JCjDp06OBhbm7e1NDQ0Ld+/foNFyxYYFXQ9QDgwYMHsv79+zvb2to2USgUvs7Ozl7jxo1zTEtLo8LuBwAajQZLly6t2apVq7oWFhZN5XK5r6OjY+P+/fu73L5926DkP51n1Go1+vfv77x06VJ7ExMT9c6dO2/16tXrhQqAtLQ0yqmRL+xa1tbW3kTk999//8kAYOvWrWZE5Pfvv/8aA0BQUFD9vPX2hw8fNu7Ro4e7r69vIwAIDw83zHs+b4nJf//9J/vqq6/sWrduXdfJyamxQqHwNTU1berj41N/7ty51mqey1XuZIWdEELMBTCXiBwA9Mp+BAJoC+B/RHQBwDYAO4QQvL5cJaQwNIL1Bxtw9s9uAAA1ERY9/gM1jzugZ5uheo6OMcaYPqWkpEgBwM7OTqnra//111/mv/76q12dOnUy2rRpk/Tw4UPFxYsXjceNG+ealJQknT59+nO1z3fv3pUHBgbWj4qKMrCyslJ16NAhMTMzU7J8+XK7U6dOmapUKiroPpmZmdS9e3f3w4cPWxgaGmoaNWqUZmNjo7xx40aNP/74w3rv3r2We/bsudm6dev0kr4HpVKJXr16ue3evbtmzZo1Vbt27brVqlWrEl+nMM7OzspevXrFHTlyxDwxMVHWtm3bJCsrK1XOeXt7e1VgYGBKZmYmHTp0yMLExETdqVOnxDznc39vW7dutZg5c2YtBweHLGdn50wfH5+njx8/ll+8eNF48uTJLseOHTPbs2fPvfwxMN0pNOHOIYSIArAYwGIisgTwJrTJdxC0SwDOJqJb0CbffwohQssxXqZjDja1MS9oLcYcHogYuQQZEsL3d+bDwsQagb5v6Ds8xhgrP1+ZF7rSRIX3VVJYeV4+NDTU8NGjRwqZTCZ69uz5wohtWf3yyy/2P//8c/hHH32UW+e8YMECq3HjxrnOnTvX8ZNPPok1MjISOeeGDx/uEhUVZdC+ffuknTt33jM1NdUAQHh4uLxDhw717t+/b1jQfT766COnw4cPWzRv3jxl48aN911dXXOT0FmzZtnOmDGj9oABA+rcuXPnikz20pToOZ07d/Y4evSouYODQ9a+fftuNWnSJLPEP4gi+Pv7p2/bti3c29u7fmJiomz69OlRr7/++nNr+TZq1Ci2U6dOKb6+vhbW1tbKbdu2hRd0rcDAwNSQkJDrbdq0Sct7/P79+/JOnTrV27t3r+Xvv/9uMXDgQC6FKSclWRYQQogEIcRqIcQbAGwAvAtgCwBHAJ8B+JuIIojof0TUlogK/MTJKhZPNx/M8p8Pc7V2pZIkqQTTL0xD6NVj+g2MMcbYKxUbGyvdvHmz2TvvvOOh0Wgwe/bsB3Xq1NH5CHe3bt0S8ibbAPDxxx/Hubi4ZKakpEhPnTpllHP86tWrimPHjpnLZDKxbNmyiJxkGwBcXV2Vs2fPfljQPSIjI2Vr1qyxNTExUe/YseNe3mQbAKZPn/44ICAgOSIiQrF9+/YS1amr1WocPXrUHAAWL14coetkW9dee+21jPzJNgC4ubkpZ86c+RAAtm7davnqI6s+SvZxLg8hxFMAmwFsJiIDAB2hHfnuCWAcgI8BfAlgtg7iZOWsZZMumJYcja9vz0WaRII4mQRTz4zBT4pV8PJoru/wGGOMlZOePXvWy3/MwMBAbNmy5XZwcLDOR7cBoFu3bgWOpLq7u2dEREQoHj58aADgKQAcOnTIBAD8/PxSPTw8Xkj+33333aQRI0Zonj59+twg4p49e0yzsrIoMDAwxcHBQZW/HwAEBASknDx50uz06dMmffr0KfZ7lUgkaNq0aer58+dNRowY4eru7n7T29u7QifdmZmZtHPnTtOzZ8+aREdHy7KysiRCCCQlJUkB4N69ewp9x1iVlTrhzksIkQVgD4A9RCSBts77bQBlXn+SvTrdA4YgLSMV3z9YikwJIUYuwaRjw7DQYCM8nL30HR5jjOlWOZdlVBYBAQHJtra2SiEEHj9+LA8NDTXNzMykESNGuHl6et7w8vLSeSLp5uaWVdBxU1NTNQCkp6fnfkP+6NEjAwCoVatWgX0kEgkcHByy7ty581xZSU4CefjwYQsiKrJ8KDY2tkT5EBHhyJEjtzt06FD3/PnzJkFBQZ6HDh2qsEl3aGio4TvvvOMRERFRaFKdmpoqfZUxVTc6SbjzEkJoABzNfrBKpnfQWKTvTcb/Hm+AiggP5YQJ+/thSc+dqG3vru/wGGOM6diUKVOie/TokZLzOiIiQt6xY8e6t2/frtG/f3+3ixcv3pBISlSBCo1GU2RJaUmvB2iT3CLOifzH1Go1AYCbm1uGj4/P0xd7PePv71/k+YKYm5tr8ibdHTt29Dx48GCpk25NOW1Ap1arkZNsd+7cOeHTTz+N8fb2zrC0tFTLZDL8/fffNVq0aNFQiBd+hEyHSpVwE5E9tHXbBU5SAAAhxOnSBsX0671unyF9RwqWJO2ChgixMg0ure2Dmh8egrGphb7DY4wxVo5cXFyUmzZtuufv79/w8uXLxkuXLq05evTo+Lxt5HK5AIDk5GQpgBfKPHJGpXXB0dExCwCyy0wKFBUV9cK52rVrZwFAo0aN0gqbTFhWxU26FQqFAAClUkkZGRlkaGj4XHabkpIiSUxM1PkgKAD8/fffNSIiIhT29vZZe/bsuSeVPj+QfePGDS4leQVK9BGTiHoT0XUAjwCcA3CikMdxHcfJXrGRb3+H943aw0Ktxoqox+iRdhvhC99ARlqqvkNjjDFWznx8fDIGDRoUCwBz5sxxVCqfz6nt7OyyAODKlSsvDLydOXOmRmxsrM52pwwKCkoFgLCwMJOC1v3euHGjeUHlED179kyWSqXixIkT5vHx8SUfUi+mnKTb19c3NTY2Vt6xY0fPf//997kkViqVwtraWimEwOXLl19IcLdt22ZW2AhzzoebwpY+VCgUGuDZiH5+T548kQHa5R3zJ9sAsGHDhkLXPme6U+x/gET0LoCNADwBJAP4F8DpQh5ndB4pe+XG91mIr43ehVeWtmyuUdYl3P65ByfdjDFWDcyePTvK2NhY8+DBA8WSJUueS8ratm2bAgDz58+3T0hIyM0lbt26ZTB06FA3Xcbh5eWVGRgYmKxUKmnEiBHOqampuYllRESE/PPPP69VUD9XV1flgAEDYpOSkqRdunSpmz8JBoCYmBjpvHnzrCMjI8s0upw/6Q4KCnoh6W7VqlUKAMyYMcMxMzMz9z2cOXOmxrRp02oXdm0HB4dCP9wAQO3atVUSiQQxMTEGeX8XORo3bpxBRLh27ZrR4cOHjfOe++GHH2wOHDjAX12/AiX5xPdZ9vM4ADZCCB8hRGBhj3KIlelB+3dn4Iz7x7mvG2dewMlFXZGY8kSPUTHGGCtvjo6OqlGjRkUDwLx58xzyjnJPnjw5xtbWVnnp0iVjT09Pr06dOtVp3rx5PW9v70aWlpbKJk2alLgmuigrV66MsLOzUx45csTCzc2tcbdu3dw7dOjg0bBhQy8LCwuVl5dXGgAYGBg8Vwj966+/PuzcuXPCuXPnTHx9fRt5eXk16Natm3v79u09PD09G9aqVct70qRJLvHx8WWeMPiypPurr76KMjIy0uzbt8/Sw8OjUdeuXd19fX3rt2nTpkGbNm2S825qk9ebb76ZCACff/65c1BQUJ2+ffu69O3b1+XatWsGAGBqaqoJCAhIysrKosaNGzd688033fr27esybtw4R0C79F+fPn2eKJVK6ty5c/1WrVrV69mzp1udOnUaTZ061XnMmDHRZX3v7OVKknDXBXBaCLFQCFHgPwpWNbV8bxbOuIwCANyVyzDbKg6jN3TipJsxxqq46dOnx1hZWakePnyoWLRokXXOcXt7e/WJEydu9OjRI16pVNKxY8fMY2JiDMaOHRt99OjROzKZTKcz8OrWrZv1zz//XOvXr98TIsKhQ4csbt++bThkyJDHISEht+Lj43PKJp7LTwwNDcW+ffvurVu37m7btm2TYmJi5AcPHrS4ePGisRACb731Vvy6devuenp66mR1kaKSbh8fn4xDhw7daNeuXVJiYqLs6NGjFk+fPpXMnj37wR9//BFR2DWHDx+eMGvWrAfOzs4ZJ0+eNN+8ebP15s2braOionLLazZs2BAeHBwcl5WVRXv27LHcvHmz9V9//ZW7rvbvv/8e8d133/3n4eGRfuHCBePjx4+bOzg4ZG3fvv3WsGHD4gq+M9MlKu6sVCJ6BCBECNG/fEOqXJo1ayZCQ6vH5poHfpuAbzX7ESfTDgR4ZcrxS78DsDC1fklPxhh7HhGFCSGalec9Ll26FO7t7c0jA1Xc1atXFY0bN/YyNTVVJyQkXCzNCiiM6cKlS5esvb29XQs6V5J/lQcAvKaTiFil1On9/6G1pH7u6ysKJUb/wSPdjDHGypdarcbJkyeN8h+/ffu2wcCBA92EEHjnnXfiONlmFVVJJgl8CeAcEX0P4DMhhLqcYmIV2Dfvbwd+64W/JLcBAJcVSoz+ozOW9NvPI92MMcbKRWZmJgUGBjZwdHTMcnd3zzA3N1dHRkYaXLt2zSgzM5Pq1auXPnfu3Ef6jpOxwhQ74RZC/EdEAQB2AuhFRIcBPARQ4ErtQohvdRMiq2i0Sffb+EtyBwBwWZGFUX90xMLee2Bj6ajn6BhjjFU1BgYGYuzYsdEhISGmV69eNUpJSZEaGBgIDw+P9J49eyZMmzbtsZmZWfnsHMOYDhQ74SbtFk9joZ08KQVQB0BBBeCUfZwT7irsm/d3PJd0X1WoMGpLV/z81k7UsnXVb3CMMcaqFJlMhoULF/IINqu0SlJSMhXARwBUAPYAuAOAF2Suxr55fwdka/piO64BAG4pNBi9syd+6rIR7rUb6Tk6xhhjjLGKoSQJ9zAAaQAChBAXyykeVsnMHLwJhuuHYIMqDABw3wDYsr0vhvbeDRtHV/0GxxhjjDFWAZRkOq8TgOOcbLP8pg1YjfcN20AiBHqmPMWkhAfIXN4ZURE39R0aY4wxxpjelSThfgTtCDdjL5jYdzGmWPTB9NhESADUEtGQ/NYVD25f0ndojDHGGGN6VZKEexOAtkRkXF7BsMqt/1szcLPNEmQK7eZXdvcV/wEAACAASURBVIiDYn0PHDu7Vc+RMcYYY4zpT0kS7lkAbgH4i4jqlFM8rJJr+vq7uB20EmlCAQFgiZUMU65/iQ375+s7NMYYY4wxvSjJpMm/oF2hpD2A60R0D4Wvwy2EEJ11EB+rhLwC38SNGiY4cGIEtplpNwabG/UbErZHY0yvuXqOjjHGGGPs1SpJwh2Ur1+97EdBClqfm1Uj9Zu9jljJT7C5MBWxMglURFiasg+Jvz/G5wPX6Ds8xhhjjLFXpiQJd8dyi4JVSYG+PfGTqQ2mHhuOBwYEANioPo+kld0xZ8hfkEileo6QMcYYY6z8lWRr98PlGQirmprUbYFfTXZg/K53cEuhrT76P9l/SFnZDv8bvB+GCiM9R8gYY4wxVr5KMmmSsVKp7VAXK/oeQtMMRe6xk4pEDF8diJg43qmXMcb0Ta1Ww8HBoTER+dWsWdM7MzOT9B1TjgULFlgRkV9wcLDrq7jfxIkTHYnIb+LEiY6v4n75BQcHuxKR34IFC6xK0k+lUmH+/PnWLVq0qGdpaektk8l8LS0tvd3d3Rt17drVfdasWbaRkZElqWwolJOTU2Mi8rt586aBLq5XHETkR0R+r+p+usYJN3slLM1tsGxICFpnWuQeu2SYheHbu+DKnb/1GBljjLEdO3aYRUdHGwBAQkKCbOPGjeb6jqk83Lx504CI/JycnBrrOxZdSkhIkDRv3tzz008/dQkNDTV1dXXN7Ny5c2KLFi1S5HK5OHDggOWMGTNqh4SEVMilnUv7IaMyKfSTDhEdBzBVCHG6tBcnotYAvhNCtCntNVjVUUNhjCXDjmHamrewVxoOAEiTqIA/+uHum7+jTpNW+g2QMcaqqVWrVlkDgK2trfLx48fy1atXWw8ePDhR33Hpw6RJkx4PGjQo3t7eXqXvWIpr8uTJjufPnzfx8PDI2LNnz+169epl5T3/6NEj2apVq2o6Ojoq9RVjWZ0/f/6qvmMoi6JGuOsDOEFEB4moLxEpimibi4gURNSPiA4BOI7CVzJh1ZBEKsX3Q3fhfcM2MFVrsCgmFl7qONhvexuXjmzWd3iMMVbtxMTESA8fPmxBRFi7du09qVSKEydOmIeHh8v1HZs+ODg4qHx8fDIcHBwqTcL9119/1QSA77///kH+ZBsAnJycVNOnT3/ctm3bSrtjuI+PT4aPj0+GvuMoraIS7roAFgJoC2ADgBgi2kNEXxBRMBG1IyLf7OdgIppORHsBPAbwO4BAAD+DE25WgIl9F+NHt+molan9ksWYMuAVMgJ/b/5Bz5Exxlj1snz5cqusrCzy9/dP6dy5c2rr1q2T1Go1li1bVujX+3nraZcvX27ZtGnT+kZGRj7GxsY+LVu2rLd//36TgvodOXLEeOTIkbW8vLwaWFlZecvlcl9bW9smXbp0cT98+HCxyx0WLVpkRUR+gYGBdQtr888//9QgIj9bW9smSqUSwcHBrvXr128MAJGRkQY57yF/icnLarjPnz9v2K9fPxdnZ2cvQ0NDXzMzs6b16tVrOGLEiFq3bt16rqZ59erVFr1793b18PBoZGpq2lShUPg6Ozt7DRo0yPnOnTs6+0ATHx8vAwB7e/sSj2BrNBosXry4pr+/v6eZmVlThULhW7t27VLF+LLabn9/f08i8tu9e7cp8KzEZ/v27VYAMG7cONe8v5e8JSZF1XBHRUXJPvzwQyc3N7dGhoaGviYmJj7e3t7158yZY6NUvvgjyTsvICEhQTJy5MhaTk5OjQ0MDHxtbW2bDBgwwDkmJkanS6kVmnALIZKEEOOhHeleAO0GN10BzASwGcBhAOeynzcD+ApAFwBKAPMBeAohJgohknUZMKs6WrR/F3Hv7kY0bAAAUhLQ3JuHT1d0gUpVab/1YoyxSmX9+vXWADBw4MA4ABg8eHAcAGzYsMH6ZX3Hjx/vOGrUKHe5XC7at2+fZGdnl3X27FnTnj171jt06NALCfQXX3zhtHLlSjulUkne3t5PX3/99UQLCwvV/v37LTt37lx/1apVlsWJefjw4fE1a9ZUnTp1yuzKlSsFfgP/008/2QDAoEGDYuVyOVq3bp3auXPnBACoUaOGplevXnE5j+7duycU576LFi2yatGiRcONGzdaCyHQvn37RH9//xQhBC1fvtxu3759pvnirLNnzx7LGjVqaFq3bp3cunXr5KysLMnvv/9u89prrzX8999/i1U98DIODg5Z2e/ZVq1WF7ufRqPBW2+95TZ27Fi3CxcuGDdu3Phpx44dE4UQ9Pvvv9s0a9asUUhISLktJ2ZmZqbp1atXXO3atTMBwNfXNzXv78XT0zPzZde4cuWKwtfXt8HSpUvtU1NTpR06dEh87bXXUm7dulVj2rRpzm3btq2bnp5e4CTg5ORkafPmzetv2rTJumHDhmkBAQHJGRkZkg0bNth06NChni4nD790tqoQ4h6ACUT0GbSj3e0ANAVgC8AcQCK0o9rnARwFcEII8dIfEGMA4NqgGZ6MOoLbK3oBFI6JdjZIlTzC41WtMDf4T9hZOek7RMZYFfbDuR8c111b51Cctl1cuzyZ23ZuRN5jk0ImuewL3/fSxBQABjUcFDX5tcmReY8N3T/U41z0uWJNUPzE75OIIV5DnhSnbXGdOnWqxo0bN2oYGxtrBg8enAAA/fv3T/zkk09UERERiv3795t07tw5tbD+q1evtj127Nj1wMDANEC72snAgQNdNm7caD1jxgzHoKCg23nbT5w4Mbp58+b3ateu/Vy5xoYNG8wHDx5cZ+LEiS69e/dOMjU1LWgX61yGhoZi4MCBsQsWLHBYsGCBzbJlyx7mPR8fHy/ZuXOnlVQqFR999NGT7Hs/6d69e3L9+vUtLS0tVdu2bQsvyc8qJCTEaPz48S4A6Mcff4wYN27cE4nk2bjl+fPnDfP3Wbp06b2+ffs+936USiU+/fRTxwULFjiMHTvW+fjx47fz9yupoUOHxs6cObPWli1brE+dOmXWsWPHRH9//6fNmzdP8/HxycgbZ14//PCDza5du2paWVmp9u3bd7NZs2YZgHbFk+HDh9des2aNbf/+/evcuXPnSo0aNXS+qaGDg4Nq27Zt4cHBwa4PHjxQDB48+MnHH38cV5Jr9OvXzz06Otqga9euCVu3br1vZGQkAODOnTvyoKAgzzNnzph9+umnjosXL35hWbRDhw5ZtG3bNuncuXM3zM3NNQAQHh4ub9myZf1r164ZrVq1yvLDDz+M18V7LfYqJUKIdCHEPiHEVCFEFyGErxCijhDCTwjRVQjxuRDiECfbrKSs7Z3hNP4wFtf0QGr2H4ULigwM3d4FoVeP6Tc4xhirwn799VcbAOjevXt8TlJoaGgo3nzzzXgAWLFiRZEfJiZPnvwoJ9kGAKlUinnz5j0CgLCwMNP8I4TvvPNOcv5kGwD69++f1LVr14SkpCTpnj17TPOfL8iECRNipVKp2Lx5s3VaWtpz9/nll1+s09LSJJ06dUp0dXXVyVems2bNclCr1TRixIjoCRMmPMmfxPr6+mb4+vo+V2M8fPjwhPwfHuRyOX7++edIGxsb5alTp8wSEhLKvGLcjBkzYj799NNIQ0NDTWRkpMGaNWtsx4wZ49asWbNGVlZW3u+9957z/fv3XygPWbx4sR0ATJs27VFOsg0AMpkMS5cufWhvb58VGRlpsHr16mJ98/Cq7du3z+TKlStGxsbGmt9++y0iJ9kGAA8PD+UPP/zwHwCsWbPGNv+/EQAwMjLSrF27Njwn2QYAV1dX5bBhwx4DwJEjR8x0FSsvC8gqBCMTc/wwLASdlM9GtP8zAMb9PQabDv6kx8gYY6xqSk9Pp5zJdsOGDXtu5PyDDz54AgB79+61TEpKKjRXCA4OTsp/zMnJSWVmZqbOysqigupgo6KiZAsWLLAaMWJErb59+7oEBwe7BgcHu968ebMGANy8ebNYZRaurq7Kzp07JyYlJUlXrFhRM++5VatW2QDAmDFjHhfnWi+jUqlw+vRpMwAYPXp0ib5l+PfffxWzZ8+2HTJkSO3evXu75rxftVpNGo0G165dK3NZiUQiwdy5c6PCw8P//fHHHyPefvvtOA8PjwwiQmJiomzdunU2Pj4+jU6cOJFbHnL37l35w4cPFRKJBB9++OELo8qGhoaiV69e8QAQEhJSrA9Br9qRI0dMAaBDhw6JdnZ2L9TS9OnTJ9nGxkb59OlTycmTJ18ojWnUqFGas7PzCx8AGzRokAEA0dHROquz18kC6IzpgoGBAvOH78NPmz/C2rSjUBIhWSrBd49W4P6Gy5jaf6W+Q2SMVTGTX5scmb/MoyTmtp0bkb/MpCRWdV51p7R9y2rdunUWSUlJUhcXl8xOnTo9zXuudevW6fXr10+/ceNGjd9++81y/PjxBX7N7+Hh8cKKGABgYmKiTk5Olqanpz+XrM+dO9d6xowZtTMyMgpN4pOTk4s9WW3cuHExe/futVy+fLltTinCrl27TO/du2fo4eGR0b1790LLYUoiKipKlp6eLpFKpcLLy6tY3+QrlUq89957Lps2bbIWovBqjMTERJ1NzrOzs1NPmDDhyYQJE54AQGRkpGzlypU1f/jhB8ekpCTp+++/73bnzp2rABAREWEAANbW1sq8I8N51alTJxMAoqKiKuSKNY8ePZIDgKura6G/k9q1a2fGxsbK//vvPwMAz/07d3JyKrBfzoh3ZmamzgameYSbVTjj+yzEV+4TUFOl/YZHTYT1yn8wdnl7pGU8fUlvxhhjxbFmzRprAEhJSZH6+fl55n/ExcXJAGDdunWFlpVIpcXPFY8fP240ZcoUF5VKRdOnT3948eLFK0lJSRfUanWYECJszJgx0QAghCj2RLVOnTo9bdCgQdqVK1eMjh8/bgQAixcvtgGAnLIAfZk9e7bdxo0bra2trZXLli27d/v27X/T0tLOCyHChBBhTZs2fQqU7P2WlKOjo2r69OmPFy9eHA4Ad+/eNbx8+bIi+74AAKLCb1/UB4XS0GiKLM0vsWK+h0JPFlbbXh444WYV0htthmFx21Wok6f8L8TgCYasbY37j27oMTLGGKv87ty5Iz979qwZoF1S7vz58yb5H7GxsXIAOH/+vIkuVtPYuHGjpRAC77///uOvv/46xtvbO9PMzEyTk/Tcu3evVPcYOXLkYwBYuHChbXh4uPzQoUMWxsbGmpEjR5Zo8l1RHBwcVIaGhhq1Wk1Xr14tVpx//vmnJQD8/PPPER988EGCh4eHMu/Ew4iICJ2sUFIcb731Vu6KcdHR0TIAcHV1zQKA2NhYeWGreNy/f18BAA4ODsWqg5fL5QIAkpOTC8wvIyMjdboVfK1atZTAszgL8vDhQwMAcHZ2LvDbmFeFE25WYXl5NMdvfQ+jeeaz5VyvK9S4+PsbuH3xhB4jY4yxym3p0qXWGo0GLVu2TMkZcS3o0bVr14Sc9mW9Z0JCggwAateu/ULiExkZKTt58mSpJqh98MEH8RYWFqrdu3dbfvXVV/ZqtZp69eoVZ2lp+cJwqkKhEACgUqlKNKosk8nQqlWrZABYsmRJsX4WSUlJzyW2ee3YscMs5+ehCy8bOb5z505uouvi4qIEgDp16ihr1aqVqdFosHTp0hfWXM/MzKQdO3bUBIC2bdumFCcOOzu7LAC4cuVKjfznzp07ZxgdHV1gwm1gYFCq30uHDh1SAODIkSMWsbGxL3zdsm3bNrPY2Fi5kZGRJiAgQK+b/nDCzSo0S3MbLBt2Em+L+iAhMCwxCW+nx8B5x9v4ZxtPpmSMsZLSaDTYtGmTFQD069evyFHgnLW5t27daqVSlW3jRU9PzwwA2Lhxo1XeiZgJCQmSgQMHuqakpJSqltnIyEj079//SUZGhmTNmjW2ADBu3LgCy0kcHBxUcrlcxMXFyQpK0IryxRdfREmlUvz66692eTdkyXHhwgXDCxcu5C4N6O7ungEACxcutMm7NvbVq1cVH3/8sXNJ7v0y3t7eDebOnWv95MmTF97TjRs3DIYPH+6a3e5p3p0oR48eHQMA3333nWPe2FUqFUaPHl0rMjLSwNHRMWvIkCHFWqc8JzH/8ccf7ePj43N/x3fu3JG///77boWVqDg6OmYBwPXr119YWrEoXbp0SfXy8kp7+vSpZOjQoc55R+rv378vnzRpUm0AGDJkyOPC6tRfFU64WYUnkUrx9ZAt+NxmEN5L0H6rpSAl/C9/iX8WDERGOtd1M8ZYce3evdv04cOHCkNDQ83AgQOLTKSCg4OTLSwsVLGxsfItW7YUa73wwowZM+aJvb191rVr14zc3Nwad+rUqU7Hjh3ruLm5Nbl8+bJR7969S73G+IQJEx7n1JP7+/un+Pn5FbgFuEKhEO3atUtSq9XUtGnThm+88YZb3759XUaPHv3STR/at2+fNm/evHBAuyOis7OzV/fu3d2DgoLq1KtXr6Gvr2+jEydO5G7289lnn0XJZDLxxx9/2NSpU8erR48e7q1bt67r4+PTyMHBQenj46Oz/3iFh4crJk+e7OLg4ODt5eXVoGvXru7dunVzb9q0af1GjRo1vnDhgrG9vX3W2rVr7+ftN2XKlNgePXrEx8bGyps3b94wMDCwbs+ePd3c3Ny8Vq1aZWtmZqbesGHD3eKuwT1p0qTH9vb2WVeuXDHy9PT06tSpU50WLVrUa9KkiZeJiYm6sPccHBycKJFIsGrVKruAgIC6ffr0cenbt6/LwYMHX7r76B9//HHPzs5OuXv37pouLi6Nu3fv7t6hQwePRo0aed2/f9+wZcuWKfPmzSv1xGhd4YSbVRp9u0/B04EHcF/iknusXuIefLS6Ff69dVqPkTHGWOWxatUqawAICgpKLKjsIi+FQiHeeOONeAD47bffylRWYmNjow4NDb3er1+/J0ZGRppjx46ZX7582bhLly4JoaGh13PqcUvDw8ND6ebmlgEAo0aNii2q7dq1a8P79OnzRK1W0969ey03b95svXPnzppF9ckxfvz4uNOnT18PDg6OU6lUdPDgQYtz586ZSiQSjBw5MqZr1665pRdBQUFPQ0JCrrdr1y4pNTVVeujQIYvo6GiDjz/+OCokJOSWTCbT2Yjr4cOHb3755ZcPAwICkrOXwDM7cOCARUREhMLX1zf1iy++eHjt2rWrTZo0eW5VDolEgp07d95ftGjR/SZNmjy9ePGi8f79+y01Gg0NGDAgNiws7Grbtm2LXYphY2OjPnHixI033ngjXqVS0bFjx8yjo6MNPvzww+hjx47dLuw9t2rVKn3FihX3vLy8nl64cMFky5Yt1ps3b7Yuzoi3l5dX5oULF66NHDkyxsjISHPo0CGLv//+29TDwyP922+//e/o0aO3y2PTnpKi4s5AJaLhANYLIdLLN6TKpVmzZiI0NFTfYVQraalJuLZsKHySD2GMnQ1OGdWAuVqDMU7D0K/zRH2HxxgrBiIKE0I0K897XLp0Kdzb21unOzOyiunMmTM1WrVq1dDGxkb56NGjf+XyCrmKHaviLl26ZO3t7e1a0LmSjHAvA/CQiOYTUV2dRMZYKRiZmMNv/BZsrDsMf9fQfvhNkkowJ2oVZqzpA5VKJ5uKMcYYqyS++OILRwD44IMPHnOyzSqikiTcuwGYAZgA4DoR7SOinlTU4oeMlROSSDBgwI+Y4TIWVtnrdWuIsAPX8f7Klrx0IGOMVXHr168379Onj0vjxo0bHDlyxMLR0TFrypQpel17m7HCFDvhFkK8AcAdwBwATwB0AvAngPtENJWIbEobBBHJiej17NHzs0QURURZRPSIiLYSUbuX9O9PRCeIKImIUokolIjGEFGR76+0/VjF8Xb7Ufi1w+9okPlsYvZFw0wM2xeMHUeX6jEyxhhj5SksLMx4y5Yt1vfv3zcMDAxM3rt37y0zMzPd7qzCmI4Uu4b7uU5EcgB9AIwG0BKAAJAFYCuAJUKIMyW8XhCAg9kvowGEQbv9ZkMAXtnHZwkhZhTQd3F2HBkADgNQAngdgCmAHQB6CyHUuuqXH9dwVwxpGU8x/fdeOCB/NhFZKgTeoAaYMWADZDL+ipGxioRruBljVY2uarhzCSGUQoj1QojWAHwArASgAtAfwEkiCiOioURU3F2UNAC2AWgjhHAQQvQQQvQVQjQG8C4ANYDpRNQ+byciCoY2aY4G0CS739sA6gK4DuBtAGPz36y0/VjFZWRojPnD9+MTqz4wUz/bEn4HbmDwyhZ4FHlHzxEyxhhjrLoqc+mEEOISgJkAfgNA2Q8fAMsBhBPRsGJc44gQ4h0hxAvbBwohNgFYnf1yYL7T07KfpwghbufpEwPgw+yXUwsoESltP1bBDekxHb8EroRn5rNfnaMqAfJlQbh6eq8eI2OMMcZYdVWmhJKIgohoO4D7AMZAW56xCkA/AHsB2AJYRkQflzHOC9nPtfLcuxYAP2hLWbbk7yCECAHwCIA9gBZl7ccqjyZ1W2Dte6fRRVUbLkolvnwSD1skoP7+/jjz2xSoy7hbGmOMMcZYSZQ44SYicyIaT0Q3AOwH8BaASACfAaglhBguhNgkhOgJoBW0tdhlTbhzliGMynPMJ/v5ahFrg5/L17Ys/VglYmRojLnD9mK623QohSkAQEoCLSOW4tIPbXDp5kk9R8gYexVKM0+JMcZKKvtvTaF/cIqdcBORLxGtgHb0dz6AegBCAAQDcBdCfC+EiM93878B7AHgXPLQc+9rD2BI9stteU65ZT9HFNH9v3xty9KPVULNXx8A5QfHcU3ulXtsp1kURp0aicXbJ+kxMsZYeSOihKysLJ4xzRgrdxkZGQZEVOgk7ZKMcIcCGJr9v1dAO9mwgxBihxCiqGV4ngKQleA+uYhIBuB3AOYADgshduU5bZLn+oVJzX421UE/VknZOrmh3uSjOOM8AvuNjLDd1ASpUgmWpuzDyGUBiH7yQN8hMsbKgUaj+b/ExET+O84YK1dCCMTGxpqq1ep1hbUpScIdDmAStGUjI4UQV4rZ7wMApR1hWArtUn0P8OKEyZwNd0r6fWFp+z27ANGI7DW7Q2NjY0t7GfYKyeQGaDl0LjL8v4a98tmv/rQiCe/92RU7Q1boMTrGWHlQq9XLYmJiEmNiYmpmZmbKubyEMaYrQgioVCppUlKSSXh4eM2EhIR/NRrN2sLal2TkuY4oxV+r7D4vXc86PyL6GcAwaJfue10IEZ2vSUr2swkKl3MuJc+x0vbLJYRYBu1W92jWrBn/Ba9E3mz3AZp7dcGM7e/ijCIZABAlJ3x5/yf8c28PpvdfD0OFkZ6jZIzpgp+fX3hYWFivqKioETExMV2FENb6jokxVnUQURqAi0qlci+AjX5+flmFtS1Jwr2fiPYJIX58yc0nAOgqhOhUgmvnv8Z8aCdaxkKbbN8uoFl49rNLEZeqna9tWfqxKsLeujaWjTiFRds+xfqk/0OqVAI1Ef6iO7i5piUmt/wR/o1f13eYjDEd8PPzC4d2Uv9neg6FMVaNlaSkJAjPdn0sSkNoy0BKhYh+ADARQByAjkKIa4U0zVkqsBER1SikzWv52palH6tixgbPwy+tf0X9PNvC31RoMCZ0HOasH8rLBzLGGGNMJ8pjYxcDaHeOLDEimgNtnXgCtMn2pcLaCiEeADiffb/eBVyrLbTrdkcDOFPWfqxqauoZgPVD/sZboj5k2RVTGRLC/eTjuDr3dUQ/4B0qGWOMMVY2Ok24iYig3VSm0GVRiug7C8AUAInQJtvFGV3+Lvv5eyLyyHMtWwBLsl/OKWAVldL2Y1WQgYECs4ZswXf1psE5CzBVazDzSTyaZF6E8YoA/LNjIYSG/ykwxhhjrHSoqHmQRHQgz8sgaDe4KazEQwbtBjWOALYKIfoWOwiiNwDszH4ZCuBqIU1vCCHm5Ou7BNrt2DMAHAKghLakxQzAnwDeEUK8MGmztP3ya9asmQgNDX1ZM1ZJJKXGY++G8ej76C9I6Nn/N0KNWsCs7/eo59JUj9ExVnUQUZgQopm+42CMsVfhZQl33mE9gWdL6hXlXwBvCiGK2lgm/32GAPitGE1DhBDtCujfH9qt5RsDkAK4Ae0W878UNUpd2n55ccJdNd345yBM/u8j1BLazU0XWJpjk6kpBlr0wIe9vtdzdIxVfpxwM8aqk5cl3DmTHwnAAWi3cp9XSPMsAI+EEPd0GmEFxwl31ZWWmoTLaybCKOkvDHK0g5q0nzdfyzTC1M7LUc+liZ4jZKzy4oSbMVadFJlwP9eQ6ASAPflLOqo7Trirvm0HF2LJf0vxWPZsyoOpWoM+Jq/j4+D/QSKVFtGbMVYQTrgZY9VJsSdNCiECOdlm1VFwx4+w/s19CMi0zD2WIpVgZfpRDFrhj4s3T+oxOsYYY4xVdOWxLCBjVY69dW38MuI4ptkPgV2ereH/NczCiNMj8e36IVCplHqMkDHGGGMVVaElJUSUsyvXL0KIhDyvi0UI8W1Zg6sMuKSk+nmSGIVvtg7GYVkkBD2bR+yZKcEXLRegaeO2eoyOscqBS0oYY9VJUQm3BtqVSRoIIW7lef3SawIQQohqUdjKCXf1tfP4Cvx68yc8MNAm3Q0ys7D60RNccvsAzfrPhNxAoecIGau4OOFmjFUnRSXcs6FNsP8nhIjP87pYhBDTdRNixcYJd/WW8jQR324ejEO4g9+jYuCZpS0rCZc4I6PzfNRv3knPETJWMXHCzRirToq9SgkrGCfcDAAuXjwE491TUVd1O/eYGsBM+9cw7I0lcHGsp7/gGKuAOOFmjFUnPGmSMR1o2jQIblNO42zdiUgT2lKSLaYm2FEjBgP2vY2fNn8EjfqlG5cyxhhjrArihJsxHZHJDdBiwJdIHnYKJ41b4ueaFgCAJKkEK9OP4d2Vfjh5Ybeeo2SMMcbYq1bshJuIPiSiLCLqXkSbHtlthusmPMYqH3vnugiYtA8fWPeBjUqTe/y6Qo2PLk3FFPa1JgAAIABJREFUlFU9kZjyRI8RMsYYY+xVKskIdy8A8QD+r4g2/5fd5p2yBMVYVTC055fY1OsgOimdIM2eK6Eiwl5pOPpsaoe1e6vFypmMMcZYtVeShLs+gMtCCE1hDYQQagCXATQsa2CMVQU2lo6YP3wffmz4FeplPvu/W5ScMDf2D7y37DWcu3pYjxEyxhhjrLyVJOG2ARBTjHaPAdiWLhzGqqYO/u9g09BQDDJoCVP1s8+sFxQZOHlwGM6smID0pyl6jJAxxhhj5aUkCXcSgNrFaOcEILV04TBWdclkckzutwxrgzYhINMSAGCnUmFkUiJaPlyFpLlNEbb3NwhNoV8iMcYYY6wSKknCfQFACyKqU1iD7HOtAFwsa2CMVVUezl74ZcRxzKnzCUbG14BRdn23PZ7A75/xCJnXBiFhO/UcJWOMMcZ0pSQJ92oAcgB/ElHd/CeJyAPAnwCk2W0ZY0XoHjAEwRPP4VyTrxEPs9zjf5pEYdzlzzFhRSdExkboMULGGGOM6UKxd5okIgKwC0A3ACoAJwHcyD7tCSAQgAzAPiFEN92HWjHxTpNMF5ISnuD6hqlQpezGhw42ucct1Rr0rBGAccELYGCg0GOEjOkW7zTJGKtOSrS1OxEZAPgfgA+gTa7zUgFYDmCiECJTZxFWcJxwM10KCfsTS8Jm4ppC9dxxlyygv/Ng9O/8qZ4iY0y3OOFmjFUnJUq4czsR2QN4HYBL9qEIAIeFENE6jK1S4ISb6ZpGrcbSndOwNWEPYmXPV335ZCgwwv8rBPj00FN0jOkGJ9yMseqkVAk3e4YTblZeEpJiMW/7SBzALWRIKPe4VAi0VdphQrclcHXy1GOEjJUeJ9yMseqkJJMmGWOvkKW5Db55fzvWtl2D1pkWoOwPx2oiHDF4jP/WdcTZ9TORlZmh50gZY4wxVpQSJ9xE5En/396dx0dV3X0c//yysiqCyC6LglUqsojiiii4oVVUxKpVH7BSbdVWa5+21i4+tWq1tdYFtQV3qwLi0rogKoogKptSRAVlXxQQEIIkIfk9f9wbCSEZmMxMZib3+3695nWZOSd3fnNeJPlyOfccs3vMbJ6ZbQgf88zsbjP7TiqKFImyA7v04b7LpnD7AddzUHFw68TAoi0cu3Uj/Rb8lS9vOYRZrzyq9btFREQyVLw3TV4CjAIKAKumSwkw0t0fTkp1WUBTSqQulZeVMeY/N9Jj7tMcVrpyh7ZHm3ajSe8RDBnwozRVJ7L7NKVERKIknmUB+wLTCK6KTwDGAJ8RBO/OwHDgLKAMOMrd309FwZlGgVvSoaR4K7PG386Bn97LnhSxxYxT27dlXV4uvYobMKLPb+jf54x0lylSIwVuEYmSeAL308DZwIXu/q8a+nwfeBwY6+7DklZlBlPglnTauO4L5j91AzOLX+G+5ts3z8lx58iSvbh8wG306NovjRWKVE+BW0SiJJ7AvRJY7u6H7aLfu8C+7t4mCfVlPAVuyQSzPnqT+97+NdMLNuK2fbZXYblzXHkHrh58Dx1ad0ljhSI7UuAWkSiJ56bJFsCnu9FvAdC8duWISG30Pqg/D1w2lb985wZ6bC349vXiHOOVvOUMe/F0/vDIeWzYtDaNVYqIiERTPIF7PbDfbvTrEvYVkTo2qN8wHh85k+vbDGe/ku1Xujfl5jDO53HW0/0Z9/RvKC2JzGawIiIiaRdP4J4GHGZmNd6JZWanA/2AqYkWJiK1d96JP+OZ4bO5Yo/TaFu6fdpYmTmnzr+bL27uwYznR1G2bVuMs4iIiEgyxDOH+2jgTYJVSB4DHgYWAU5wVfsi4EIgF+jv7pEI3ZrDLZluy9Yi/v7MVbz4zXQu3biRi77e9G3bkpwOLOvzE/qdNIK8vPw0VilRozncIhIl8a7DfSXwV6q/Mm4EYfxn7n53csrLfArcki3WrF/J/Of/ziGLHmFPir59/Wf77M3C/EYMaTOUS065npzc3DRWKVGhwC0iURJX4AYws17AT4FjgbYEQXsFwdXvO919drKLzGQK3JJtvt6wjnnjb+bgpY+xtHAbw9ptX1Bo/2JjaMcfcN7AaxS8JaUUuEUkSuIO3LIjBW7JVhvWrmbMhCv5V958tubsuHHsgcW5DNv/Ms4+/oo0VSf1nQK3iESJAneCFLgl23265ENGvXotb+WuoqRK8O5enMfZXUYwdOBP0lSd1FcK3CISJQrcCVLglvpi7oLp3D/5V0zNX8M22/mK94X7jeT0ASOxnHgWNxKpngK3iERJjYHbzB5I4Lzu7iMT+PqsocAt9c2MeZP559QbmF6wnrJKwfvRlaspLN+fbUdfx3ePOVPBWxKiwC0iURIrcJcncF5390jccaXALfXVjHmTGTP1d7xTsI5Dt27lH6vXfNv2Sd4BbDniGg45bqhurpRaUeAWkSiJFbhHJHJidx+dyNdnCwVuqe/mfPI2iyb9ncFfTqLAyr59/cmmTXiyaQtOaz2E4YN/q+AtcVHgFpEo0RzuBClwS1SsXraQJc/9iZ5rnseslFM7tOWLvDwAOpXAqS0GM+K0P1BQUJjmSiUbKHCLSJQocCdIgVuiZs3KxUx8/lfcUTiX4iqrmrQtdQY1OZqRp99C08bN0lShZAMFbhGJkloFbjNrAhwKtASWuvu7yS4sWyhwS1R9umQO97/6S97KXb7TOt4ttpVzfP4hjBx8G61atEtThZLJFLhFJEriWmbAzJqGq5esBV4DngRGVmq/3MyWmtlhyS1TRDJNt449+culL/PUwKc4tawTTcq232e9Li+HsT6XIc+dyI1jzuarL1eksVIREZH02u3AbWaNgMnApcDXwKsE27pXNhFoDwxJUn0ikuG6dOjOrcNfYML3XuYc607zbduD96bcHFp9/R4N7+nJ9HsuZfXSBWmsVEREJD3iucJ9LdAL+BfQ2d1PrtrB3T8DFgDHJ6c8EckWrffuwO8uepLnh03hBwVH0KbUaVReznlfb6KhldBvzVhajD6c9+8YxswPJqW7XBERkToTT+A+F1gFjHD3ohj9lgCatCkSUXs2ac4vvv8AL1w0g182v5C11vnbtnwro8WWSQyf/VMuuv9Qxk66i/KyshhnExERyX7xBO79gPfcfesu+q0F9q59SSJSHxQWNGDIGb+iy/Uz+bD/aOYVHAzA6D33oNyM2Q2KuXHFA5w1uhd3j7uWLVtj/TteREQke8UTuEuB3Vlgtz2wuXbliEh9Yzk59BhwDt1//TYfnTqWVfktd2j/rNC5v2gipz1+GDc+cj6r1ixNU6UiIiKpEU/g/hToZWY1hm4zawYcAvw30cJEpP456LATGTPyPe7v+Wf6l7SgsHz7sqRrwpVNzvz3KVz7z5P48NNpaaxUREQkeeIJ3OOBVsCfYvT5I9AEGJtIUSJSvx15yCnc/cPJPD3wKU4v359mlZYU3JKTw8T8lVw07TJevON0FsyZksZKRUREErfbG9+YWWNgBtANeJsggP8NeINgPe6hwAnAPKCvuxenouBMo41vRBK3YdNa/vHv63m9aCrL84PVRntv3crDq74EYF7BwZQc+iN6HH8eueF28pLdtPGNiERJXDtNmlkHgqB9KOAE63BXnMCAOcAZ7r4syXVmLAVukeTZtq2Ux165hZdWPMPIDas4fss3O7Q/3KQNS1r1YfhJN9F+n07pKVKSQoFbRKKktlu7nwacCnQBcoFlwEvAeHcvj/W19Y0Ct0hqLJgzhY2v/42eG18nz8opAwa3b8uK/DwalDv9trVk2KHXcnSv09JdqtSCAreIREmtArdsp8Atklqrly1k0Ut3snbDC/y6VdOd2g/ems+J7Ydw4Um/JC8vPw0VSm0ocItIlNQYuM1sHDAaeNmVymukwC1SN9ZvXMPol25g8qapLCnYub1NqXNsYW9GnPJH2uy9b90XKHFR4BaRKIkVuMsJ5mevAh4GHnL3BXVYW1ZQ4BapW+VlZYx7/R5e/PwxZhduodxsh/ZG5eUMKu3AyAG30qHrIWmqUnZFgVtEoiTWsoCjgPVAW+CXwMdm9qaZXWxmjeqkOhGRKnJyczl30FU8NPI9HjzsHk4sbUfTKssKdtgyjw6PH8vcmwcw65VH2VZaksaKRUQk6mLO4TazAuAMYDgwkOAGSQeKgKeAB9090rtT6Aq3SPqt27Ca0S/ewOSid1idD68uXUGL8u0h/Eua80jHIxjc/1oO7NInjZVKBV3hFpEoiWcd7rbAxcBFwAHhy06wA+UY4FF3X52KIjOZArdI5igvK+OVN8bQZtZ4ehRNJ8eCn2+f5edxZvu25LrTq6QxJ3U6j3NPuIqc3Nw0VxxdCtwiEiW1XRbwCIKr3kOBPQiCdxnB0oBjgH+7e1kS68xYCtwimWnVkk9YPPFeuq2YwAPNc3hizx1XOGlf6hxV2IuLT7yRDq06p6nK6FLgFpEoSWhZQDNrCJwDXAIcV6lpjbu3TqiyLKHALZLZSoq38tB/fs+ktS8zv3Dn6wCF5U7f0uac8d2RnHzkBWmoMJoUuEUkSpK2DreZDQIeA1oC7u6R+L9aBW6R7PHmzGcZP+vvvJv3BVtydr5nfP9i45xmgznzlP+lcdNmaagwOhS4RSRKEr3C3QQYRnCF+0iC7d0Blrp7p0SLywYK3CLZZ836lTz40u+ZUvQOi6us6f30ilV0KM5jXotBND/mh+x/yNFYNeFcEqPALSJRUts53AOA/wHOAhoSBO1i4HmCOdwTo7JZjgK3SPYqLyvjhbcf5D+fPMiMgo10LSnhqZVf7NBndn5nXu7Uix8M+h3t9+mUnkLrIQVuEYmSeFYp6UywSsnFwL5sv5o9B3gQeMzd16eiyEymwC1SPyxaPp8Zkx/gsM8n0rF8+bevP7JHU25rsReF5U7v0mac3PUCzjz2Mq1wkiAFbhGJkl2tw92IYCWSS4BjCEK2EWyI8wQw2t3npL7MzKXALVK/eHk5H7//Kpunjab7hjc4v31zPivYcd5J+1LniIIeXDjgerp06J6mSrObAreIREmsrd1HE4TtxgQhuxx4jWDKyAR319ZtKHCL1Gfr1n3BP17+FdO2vM+igp3b89zpVdyEgZ2Gcu4JV5GXl1/3RWYpBW4RiZJYgbtim7ZFwEMEu0our7ZzhClwi9R/5WVlvDL9Cf49bzQz8tdUu8LJPtvKudqPoc/xV9Kui65674oCt4hESazA/Sgwxt3fqNuSsosCt0i0rFm/kkcn/pGpG6fyaeH27eNbbtvGxGUryQM+KjiYzQeeS/eBF2l5wRoocItIlCRtHe6oUuAWia43Zkzgudn38H7uKs7ZtImfrd+4Q/tLDffgub06cnznczlrwOWaclKJAreIRIkCd4IUuEVkU9EGPnzraRp/OJ7vbnmPvHBG3o9atWRqo4YAtC51+uZ25ZwjrqH3d45JZ7kZQYFbRKJEgTtBCtwiUtna1UtZOGk0BUvG8z/toNxspz7di/Pot9ex/GDQr2jRrHUaqkw/BW4RiRIF7gQpcItIdSputHzpoweZkfsFm3J3vtGyUXk5vUqbc1ynszh7wI/Jz69mKZR6SoFbRKJEgTtBCtwisisbN3/FExNvZeqaScwtLN7pqnfD8nKeXvoNq1qfQuujL6Jz98PTVGndUeAWkShR4E6QAreIxOOjz2fy1JQ/817pPJbnB8F78OYiblmz7ts+n+V2Zn6ngXQ+8jy671c/M6kCt4hEiQJ3ghS4RaQ2ysvK+M/Uh5n0yeN8b93nnLD1qx3ab9i7Oc81aUz3knwOb3Y03z/hl7Rq0S5N1SafAreIRIkCd4IUuEUkUaUlxXz09gS2zX6K7l9PwXO2MWDfdhRV2mCnsNzpWdqUo9udyrknXEOjBo3TWHHiFLhFJEoUuBOkwC0iybRp41dMmXQfj64by7yCUryaVU6alZXTu6wNAw84n8FHXUxObm4aKk2MAreIRIkCd4IUuEUkVeZ9NoNxU//CjOL/sriGBUzalpZzfXl/2h95MZ0P6otVs+18JlLgFpEoUeBOkAK3iNSF198fx4sfjmYmS1mbtz1UH/7NVv65+ksAluR0YGX7U2lz9Hl06tY7XaXuFgVuEYkSBe4EKXCLSF0qKSnmmTdHMXnReD7I+4qff7WeszcX7dDnjr32ZHKjPeld2J0h/a6iR7cj01RtzRS4RSRKFLgTpMAtIumyqWgD86e9QM6HEzjo67dpZMU4cHL7tqzMz/u23wHFOfRp3Iuzj76Gbh17pK/gShS4RSRKFLgTpMAtIplgy+aNzH9rHF/Mf5LftFhFaTU3W+a4c1BJPofu0Y9zj7uODq27pKHSgAK3iERJ5AO3mZ0PXA70AHKBj4EHgVHuXr6rr1fgFpFMs3LNEp5+43beXz+NeYXFlFUTvvPc6V7cgJEtz6dH/wvYs0WrOq1RgVtEoiTSgdvM7gGuALYCrwGlwAlAU2ACMNTdy2KdQ4FbRDLZ58vmMXbKHczcPIOPC7btsMzg/iUlTFixmlLPZX7DXmztejpdjx3GXi3bpLwuBW4RiZLIBm4zOxsYB6wGjnX3BeHrrYA3gAOBn7r7nbHOo8AtItnivwvfZcI7dzHrmw9ZWOhcsX4Dl2/4eoc+jzTdg5ebtKTnHodyztFX06VD95TUosAtIlES5cA9A+gDXOzuj1Rp6w9MJgjj7WJNLVHgFpFsNOujN1k780U6LXmdbts+/fb1C9q04sMGhQDkunNgST49m/RhyFFXJfWGSwVuEYmSSAZuM2sPLANKgGbu/k01fZYD7YCj3H1aTedS4BaRbLdqyScseftJWPIil7UtqnZ3yxx3vlOSR8/GvTjryKs4oHOvhN5TgVtEoiSqgft04HlgtrtXuzuEmU0AzgR+4u731HQuBW4RqU/mLpjOc9Pv5YNvPuCTgrJqw7e5c0BJHiOanEbv/iPYp13nuN9HgVtEoiRv113qpYrfDkti9Flapa+ISL13cNd+HNy1HwAffT6DZ6fdzZwtc3a44dLNWJNbwqBP7+K9rcXsM/y2dJYsIpLxohq4m4THohh9NofHplUbzOwy4DKAfffdN7mViYhkiIO6HMpBXR4C4NNFc3hm2l3MKZrF/IJSBm7ZQi7Q5ohhaa1RRCQbRDVwV/wfaa3m07j7A8ADEEwpSVZRIiKZqlvnnvyy82gAFiyZy8Lp43mv7BMOO1CzQkREdiWqgXtTeGwSo09F26YYfUREIqdrx4Pp2vHgdJchIpI1ctJdQJosDo8dY/TpUKWviIiIiEjcohq4Z4fH7mbWsIY+fav0FRERERGJWyQDt7svA2YBBcDQqu3hxjftCTa+eaduqxMRERGR+iSSgTt0c3i81cz2r3jRzPYB7g2f3hJrl0kRERERkV2J6k2TuPs4MxsFXA7MNbNJQClwArAH8CxwdxpLFBEREZF6ILKBG8DdrzCzt4EfA/2BXOBjYAwwSle3RURERCRRkQ7cAO7+BPBEuusQERERkfrJ3LVvSyLMbA2xt4iPZW9gbRLLiQKNWXw0XvHReMUnkfHq6O4tk1mMiEimUuBOIzOb4e7api0OGrP4aLzio/GKj8ZLRGT3RHmVEhERERGRlFPgFhERERFJIQXu9Hog3QVkIY1ZfDRe8dF4xUfjJSKyGzSHW0REREQkhXSFW0REREQkhRS4RURERERSSIE7SczsfDObYmYbzWyzmc0wsx+b2W6PsZnlm9kJZvYXM5tuZqvMrMTMVpjZODM7LoUfoc4lY8xinPtPZubh4+fJqDfdkj1eZtbQzH5hZu+b2QYz22Jmi8xsrJkdlez661oyx8vM2pvZXWb2iZl9Y2ZbzWyBmd1nZl1SUX9dMbMDzOxqM3vMzD42s/Lw++acBM+bsu9vEZFsozncSWBm9wBXAFuB14BS4ASgKTABGOruZbtxnoHAq+HT1cBMoAg4CPhu+Pr/uftvk/oB0iBZY1bDufsC7xD8g9KA69z99mTUnS7JHi8z6wxMBPYHvgSmA8VAJ6AncKO7/zGJH6FOJXO8zKwX8DrQDFhO8H0JcCjQDtgMnOTu05L5GeqKmf0NuLqapqHuPq6W50zZ97eISFZydz0SeABnAw6sArpWer0V8FHYdvVunut4YBxwTDVtw4Bt4fkGpPtzZ8qYVXPuQmAesILgF7sDP0/3Z86k8QIaAwvDr7sRyK/S3gLolu7PnUHjNS38mgcqjxWQD4wO2z5I9+dOYLwuBf4MnAvsB0wOP9M5mTD+euihhx714ZH2ArL9AcwIf4FcVE1b/0q/eHKS8F7/DM83Ot2fO1PHDLg1/PrTgYfqSeBO6ngBN4df83C6P1umjxfQIOzvQOtq2ttWam+U7s+epPFLNHDX2c9EPfTQQ49seWguXQLMrD3QBygBxlZtd/c3Ca60tgb6JeEtZ4fH9kk4V1qkcszM7HDgWuAJd38h8WrTL9njZWYFwA/Dp7ckr9LMkIK/X2UE/7MEwfSknU4ZHouAb+Ktt75Jw89EEZGsoMCdmF7hcZ671/TL9v0qfRPRNTyuSsK50iUlY2ZmDYCHga+ofj5qtkr2ePUhmDKyzN3nm9mR4Q2m95vZH8zsiEQLTrOkjpe7lxLMQQb4g5nlV7SFf66Y5z7a3XVDTN3/TBQRyQp56S4gy3UOj0ti9FlapW+tmFlr4JLw6fhEzpVmqRqzm4ADgPPcfW1tCstQyR6vg8PjAjN7CLi4SvtvzWw88IMYgSmTpeLv1xXAywT/M3CKmc0IX+8L7AXcCVwXZ531VZ39TBQRySYK3IlpEh6LYvTZHB6b1vZNzCwPeAzYE3gty6dLJH3MzOxI4KfAs+7+VAK1ZaJkj1fz8HgskAvcDtwHrAtfu5fgprevgeHxFpsBkv73y90/D/+OPQKcwo5TumYAb4VXwqWOfiaKiGQbTSlJTMWczlT/V/J9BEtqLQMuTPF7pVpSx8zMGgIPEgTEK5JxzgyT7L9jFd/zeQTTIK5z98/cfYO7Pw+cGb7XxVm6vnTSvyfDsP1fgiUUzwD2BloSjNVewHgzy/qlOpOkrn4miohkFQXuxGwKj01i9Klo2xSjT43M7E5gBMG63Ce4++ranCeDJHvM/gR0A65x92ye216TZI9X5T7/qNro7jMI1pnOAY7bjfNlmqSOl5k1A54luBp7srs/7+7r3H2tuz8HnExws+QNZtY11rkiIuU/E0VEspECd2IWh8eOMfp0qNJ3t5nZX4CrgDUEYXtBvOfIQIvDY7LGbAhQTnBFdnLlB0EYArg8fO2ftag33RaHx2SNV+U+i2roU/F66904X6ZZHB6TNV6DCa5mT3f3z6s2uvtC4F2C/zE4bneLrMcWh8eU/EwUEclWmsOdmIpl+rqbWcMabjLrW6XvbjGzPwPXEMytHeTuH9W+zIySijHLIVjftyZdwkez3TxfJkn2eM2q9OcWBP+Yq2rv8Li5mrZMl+zx2jc8bozRZ0N4bB6jT1Sk7GeiiEg20xXuBLj7MoIAUwAMrdpuZv0JbrBaTbDV+G4xs1sIVj1YTxC2P0hKwRkg2WPm7p3c3ap7ECwTCMHW7ubuPZP3SepGCsZrBcEVWQjuC6h6vr2A3uHTGVXbM10KvidXhsc+lZcErHS+fIKlFqHm/zGIjFT9TBQRyXYK3Im7OTzeamb7V7xoZvsQrPgAcIu7l1dqu9nMPjazm6nCzP4P+F+Cq2aD3L0+XgVK6phFQLLH66bw+Fsz61npaxoAowhWw5lJ9gaiZI7XS8AWgivdd5hZYaWvKQT+TjBFYj3wStI/SYbaxd+vuMdfRKS+05SSBLn7ODMbBVwOzDWzSUApwdXDPQhuuLq7ype1IVgzuk3lF83se8BvwqcLgSvNqtvcjo/dPWt3CUzmmEVBssfL3V8ws9uBnwPvmtm7BFOXDiPYqnwF8P1s3cglmePl7l+a2RXAaODHwBAzm0mwGkefsH8xMNzdY007yVhm1pvtQRjgoPD4JzP7ecWL7l55Z8hYf79qM/4iIvWaAncSuPsVZvY2wS/k/gTrG38MjAFGxXElp/Ic0EPDR3XeJMu35U7imEVCssfL3a8zs2nAlQQ7/jUi2JDkrwRXH6ub2501kjle7v6wmc0lWOv9GODEsGkFQRD/a5bfY7EHcHg1r9d61RV9f4uI7Miy9CKWiIiIiEhW0BxuEREREZEUUuAWEREREUkhBW4RERERkRRS4BYRERERSSEFbhERERGRFFLgFhERERFJIQVuEREREZEUUuAWySJm1sHMHjezlWa2zczczP5WTb+fhW0X1UFNOWb2YzObYWabzWyjmU0xs++n+r1FRESygXaaFMkSZmbAeKAv8BHwBsGW2e9V0/0sYBvw7xTXlAs8A3wP+BqYCBQSbOP9hJkd4e5XpbIGERGRTKedJkWyhJl1Bj4n2IJ9P3ffVkO/VsBK4HV3H5Timq4Fbif4B8Dx7v5F+HpXYArQCjjT3Z9LZR0iIiKZTFNKRLJHh/C4qKawHTqT4Ht7QiqLCa9u/yJ8enlF2AZw9wXA/4ZPr09lHSIiIplOgVskBcL50x7++ZJwfnORma02s9Fm1jJsa2BmfzCzT81sq5ktNbObzCy/0rk6hed6M3ypf8X5K96jirMAB56tdI7fh/1/b2btzewhM1tlZlvMbJaZnVOp71Fm9qKZrQvb3zCzvtW8zxHAPsByd3+rmvaxBFNe+ppZu7gGUEREpB5R4BZJITO7Fbgf+Ap4mSAIDwcmmVkT4DXgSmAe8DrQAvg1cE+l02wGHgZeCZ9/ET6veFR+vz2BAcB0d19ZTUkdgZnAMQQBfhbQC3jazM4zsyEEc8P3Bl4FlgDHAW+YWbcq5+oVHt+v7rO7+5bwcwH0rK6PiIhIFOimSZHUuhjo6e7zAcxsL+AdoEd43AB0dveNYXtPggB7qZnd5O5L3H0tcImZHQecBHzs7pfU8H6nA/nUPJ3kEuBO4Fp3Lwvf83LgXuA2oDFwgbuPDdtygCdpE4YQAAACRklEQVSAYQRTREZUOlfn8LgkxudfShC2O8foIyIiUq/pCrdIav22ImwDuPt64L7w6UHAZRVhO2yfA7wIGNC/Fu93Vnh8pob2JcAvKsJ26AFgHdAeeLkibIf1lAO3hk8HVDlXk/BYFKOezeGx6S7qFhERqbcUuEVS6+VqXlsYHpdUDuOVLAiPbeN5IzNrRHAFfK67f1ZDt9fdvaTyC2H4Xhyj3prqsYpTxFOniIhI1Chwi6TW8mpe2xyjrXJ7gzjf62SgETVf3d6d99yp3d0r2gqrNG0Kj02oWUXbphh9RERE6jUFbpEUCqdk1CRWW23sajrJ7rxnPDUtDo8dY/SpWMpwcYw+IiIi9ZoCt0g9EC4jOBj43N0/rKO3nRUeq1sysGKKy3fDp7PrpCIREZEMpMAtUj8cDzQj9tXtZHsH+BJob2bHVtM+lGDFlPfdfUUd1iUiIpJRFLhF6oeK6SQp3V2ysvBmy9vCp6PMbJ+KtnBr91vCpzfVVU0iIiKZSOtwi2S5cK3sM4BVBFed69IdwLEE638vMLPXCK5qDyS46fMud3+ujmsSERHJKLrCLZL9jgJaAc+6e50u0Rde5T6TYLfMhQTLEvYn2M3yAne/qi7rERERyURWx7+fRSTJzOwO4KfAie7+arrrERERkR3pCrdI9psP/B54I811iIiISDV0hVtEREREJIV0hVtEREREJIUUuEVEREREUkiBW0REREQkhRS4RURERERSSIFbRERERCSFFLhFRERERFJIgVtEREREJIX+Hxzynss85XxEAAAAAElFTkSuQmCC\n", | |
"text/plain": [ | |
"<Figure size 432x288 with 1 Axes>" | |
] | |
}, | |
"metadata": { | |
"needs_background": "light" | |
}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"plt.plot(num_heun[:,2]/0.25, num_heun[:,1],'-', label = 'Heun')\n", | |
"plt.plot(num_rk2[:,2]/0.25, num_rk2[:,1],'-', label = 'Runge Kutta')\n", | |
"plt.plot(d_m,vf,'--', label = 'Analytical Solution')\n", | |
"plt.ylabel('Velocity (m/s)')\n", | |
"plt.xlabel('mf/m0')\n", | |
"plt.legend(loc=9, bbox_to_anchor=(1.5,1));" | |
] | |
}, | |
{ | |
"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": 77, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"def rocket(state,dmdt=0.05, u=250,c=0.18e-3):\n", | |
" '''Computes the right-hand side of the differential equation\n", | |
" for the acceleration of a rocket, with drag, in SI units.\n", | |
" \n", | |
" Arguments\n", | |
" ---------- \n", | |
" state : array of three dependent variables [y v m]^T\n", | |
" dmdt : mass rate change of rocket in kilograms/s default set to 0.05 kg/s\n", | |
" u : speed of propellent expelled (default is 250 m/s)\n", | |
" c : drag constant for a rocket set to 0.18e-3 kg/m\n", | |
" Returns\n", | |
" -------\n", | |
" derivs: array of three derivatives [v (u/m*dmdt-g-c/mv^2) -dmdt]^T\n", | |
" '''\n", | |
" #dstate = np.zeros(np.shape(state))\n", | |
" dstate = np.array([state[1],u/state[2]*dmdt-9.81-c/state[2]*state[1]**2, -dmdt])\n", | |
" return dstate" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 78, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"y0 = 0 # initial position\n", | |
"v0 = 0 # initial velocity\n", | |
"m0 = 0.25 # initial mass\n", | |
"dmdt = 0.05\n", | |
"time = (0.25-0.05)/dmdt\n", | |
"t= np.linspace(0,time,10000)\n", | |
"dt=t[1]-t[0]\n", | |
"N=int(time/dt)\n", | |
"mf = np.linspace(0.25, 0.05, N)\n", | |
"\n", | |
"#initialize solution array\n", | |
"num_heun_drag = np.zeros([N,3])\n", | |
"num_rk2_drag = np.zeros([N,3])\n", | |
"\n", | |
"#Set intial conditions\n", | |
"num_heun_drag[0,0] = y0\n", | |
"num_heun_drag[0,1] = v0\n", | |
"num_heun_drag[0,2] = m0\n", | |
"\n", | |
"num_rk2_drag[0,0] = y0\n", | |
"num_rk2_drag[0,1] = v0\n", | |
"num_rk2_drag[0,2] = m0\n", | |
"\n", | |
"d_m = mf/m0\n", | |
"vf = -250*np.log(d_m)\n", | |
"for i in range(N-1):\n", | |
" num_heun_drag[i+1] = heun_step(num_heun_drag[i], rocket, dt)\n", | |
"for i in range(N-1):\n", | |
" num_rk2_drag[i+1] = rk2_step(num_rk2_drag[i], rocket, dt)\n", | |
" \n", | |
"#print(num_rk2[:,1])" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 102, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAtwAAAEaCAYAAAAi8IBdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd1hURxcH4N9soVfpoHRFBURAsWJBxB4TiRpLPk00GnuJJRo1tiS2mFhji9Ekdo0xlliwYC+ADXsBVLr0zpb5/lggSJOyuCDnfZ591nvvzL1nF0LOzp47wzjnIIQQQgghhFQPgaoDIIQQQggh5H1GCTchhBBCCCHViBJuQgghhBBCqhEl3IQQQgghhFQjSrgJIYQQQgipRiJVB1DbGRsbc1tbW1WHQQghtUpwcPBrzrmJquMghJB3gRLuKrK1tUVQUJCqwyCEkFqFMRah6hgIIeRdoZISQgghhBBCqhEl3IQQQgghhFSjGptwM8a+Z4zxvMe0MtoNZoxdYIylMMbSGWNBjLFxjLEyX1tl+xFCCCGEEFIRNTK5ZIy1BDADQJnrzjPG1gHYAaAFgAsATgFoBGAtgP2MMaEy+xFCCCGEEFJRNS7hZoypA9gGIBbAoTLa+QMYCyAGQDPOeW/O+UcAGgJ4AOAjAOOV1Y8QQgghhJDKqHEJN4CFAJoC+BJAShntZuU9z+ScP8nfyTmPBTAmb/PrEkpEKtuPEEIIIYSQCqtRSSVjrBWArwDs5JwfLqNdfQCeAHIB7Ct6nHMeCCASgDmA1lXtVx3SctPww7Uf8DL1ZXVehhBCCCGEqFiNSbgZYxoAtgNIBDDpLc3d857vcc6zSmlzo0jbqvRTqstRl9H3777Y+XAnvrv2HTgvs1SdEEIIIYTUYjUm4QbwHQAnABM456/f0tYu77mshRNeFGlblX5Kpa+mj4TsBADApahLOB5+vLouRQghhBBCVKxGrDTJGGsLYDKAvznne8rRRSfvOaOMNul5z7pK6PcGxtgoAKMAwNrauoxTlczZ2BmfOH2CnQ93AgCWXl+KtpZtoa+uX+FzEUIIKV1wcLChQCAYJhQKh3DOjQEwVcdECHk/MMbSOecXpFLpPgBnPT09Sy1ZUHnCzRjTBPAbgFQoZg8pV7e854rWYlS23xs455sAbAKAFi1aVOpcE9wnICAiAHFZcUjITsCKoBVY1G5RVcIihBBSSHBwsJpIJNpWr149V2Nj43QNDY0ExijfJoRUHeccUqlUmJqa6hsfH981KytrS3Bw8NLSku6aUFLyPRRzYE/lnEeXs09a3rNOGW3yj6UV2lfZfkqno6aD2a1nF2z//fRvXI66XJ2XJISQumaQnp6ea/369RM1NTVzKdkmhCgLYwxisVhmZGSU4ujomKqpqTkSQOfS2teEhPsjAHIAwxhj5wo/AHTPazMmb9+WvO3wvGebMs7boEjbqvSrFl2su8DPxq9ge+GVhciUZFb3ZQkhpE4Qi8U969WrR4k2IaRaiUQiuYmJCUQiUf/S2tSEhBtQxNGxhIdZ3nH7vO0Weds3856d80pSStKySNuq9Ks2s1rNKqjdjkyPxJqba97FZQkhpC5opq2tTaMYhJBqp6enl84Y8y7tuMoTbs65LeeclfSAYppAAJiet695Xp+XAEIAqAEo9mmCMdYRQH0oVpO8UuhalepXnYw1jTGz5cyC7R0PduBW3K13cWlCCHmvcc41hUKhXNVxEELefyKRSMo5L7VkWeUJdxX8kPe8lDHmmL+TMWYKYH3e5hLOedE/tpXtV2162/dGO6t2AAAOjnmX5yFHlvOuLk8IIe8tKichhLwLb/tbU2sTbs75fgC/QLEq5F3G2GHG2F8AnkCxNPzfANYqq191YoxhXut50BJpAQDCUsKw8fbGdxkCIYQQQgipJrU24QYAzvlYAEOgKBPpCKAbgKcAxgPw55zLlNmvOlnqWGKy5+SC7a2hW3En/s67DoMQQgghhChZjU64OefD82q3V5TRZifnvB3nXI9zrs059+Scr3tbSUhl+1WngU4D0cJMcV+ojMvwzcVvkCUtbQV6QgghhBBSG9TohLuuETABFrdfXFBaEp4ajp+Df1ZxVIQQQt5XVlZWrowxzyNHjpS6ujIAeHl5OTHGPFevXm30rmIj5H1CCXcNY6VjhZle/81asvPhTlyNvqrCiAghhBBCSFVQwl0DfeT4ETrW71iwPefiHKTmpqowIkIIIYQQUlmUcNdAjDHMbzsfBuoGAIDYzFgsvb5UxVERQgghhJDKoIS7hjLWNMbc1nMLtv959g9ORZxSYUSEEEJIcWfOnNHu3bu3vZmZWTOxWOxhaGjo5uPj43jixIlii4A8evRIjTHmaWVl5Vra+Rhjnowxz7L2b9682bB58+aNtbS03LW1td3btGnTqKTrEVJTUMJdg/nZ+qGXfa+C7W8vf4vo9GgVRkQIIYT859tvvzXz9fVtfOzYMUMTExOJr69vso2NTU5gYKB+z549nX788UdjZV9z8uTJll9++aW9WCzmnTt3TjEzM8u9evWqbp8+fRoFBARoK/t6hCiDSNUBkLLNbjUbIbEhiM6IRlpuGr6+8DV+7fYrRAL60RFCSFXYfn202ChqbRG+pFewqmPYv3+/3sKFC+ubmJhIdu/e/czHxycj/9jJkye1/f39G3799dfWXbt2TWvWrJnSlk/etm2b6blz5x54e3tnAoBMJsPQoUNtdu/ebTxv3jxLX1/fJ8q6FiHKQiPcNZyemh6WdlgKAVP8qELiQrD5zmYVR0UIIeR90qdPn0b5JRslPW7cuFGsXGPhwoWWALB27drwwsk2APj5+WVMmTIlWiqVsjVr1pgoM9YZM2ZE5ifbACAUCrFixYpIAAgODtbNyckpe41tQlSAhklrAXdTd4xxG4N1t9YBADbc2QAvCy94mtXawRlCCCE1SPv27VNNTU0lpR0PDAzUT0hIKMgZoqOjRaGhodo6Ojqyfv36lTiNVpcuXdIWLVqEoKAgpdZW+/v7pxTdZ2VlJdXT05OlpqYKY2NjhdbW1lJlXpOQqqKEu5b4wvULXI2+iuDYYMi5HF9f+Br7++yHvrq+qkMjhJBaqSaUZdQUM2fOjOndu3daace9vLycEhISChLnx48fq3HOkZ6eLhSLxWWO/iQmJio113B0dMwtab+Ojo4sNTVVmJWVRd/ekxqHEu5aQigQYon3Enx8+GOk5KQgJiMG8y/Px8pOK8EYfXtGCCHk3ZFKpQxQJLl+fn7JZbU1MjIq92izTCZ7axuhUFje0xFSY1DCXYuYa5tjQdsFmHx2MgAg4EUA9jzag08af6LiyAghhNQl9vb2uQAgEon4gQMHwsvbT11dnQNAZmZmiaPQT548UVNKgITUMPS1Sy3TxboLBjoNLNhedmMZQl+HqjAiQgghdY2dnZ2kYcOGWcnJyaIjR47olrefhYWFVCwW8+TkZFFUVFSxQb+DBw9SnSR5L1HCXQtNazENjes1BgBI5BJMPTcVydllfqNHCCGEKNW8efOiAGDEiBF2f/31l17R49nZ2WzHjh36hefGVldX5y1atEgHgOnTp1vK5fKC9idOnNBZunSp1TsInZB3jhLuWkhDpIGVHVdCV6wYVIjOiMasi7Mg5/K39CSEEEKUY+jQocnffvvtq4SEBLG/v39DW1tbFx8fH8fu3bvbN2vWrLGJiYnb0KFDHUNCQrQK91uwYEGkWCzmO3fuNGnYsKFzjx497F1dXZv07NnTadiwYXGqej2EVKcKJdyMMQPG2IeMsQWMsQ2Msd2MsV/ytvsyxgyqK1DypgZ6DbC4/eKC7YuRF7Hl7hYVRkQIIaSumT9/fuyFCxfuDxgw4LVcLsfly5f1Lly4oJ+amiry8vJK+/HHHyOGDRuWWLhP165dMw4fPvy4TZs2aTExMWrnzp3TB4C1a9eGrVq1Kko1r4SQ6sU452U3YEwIoB+AsQC8AeRPiVF4agxe6Pk8gPUADnLO3367cS3XokULHhQUpLLrrwxaid/u/QYAEDABNnbdiNYWrVUWDyGElAdjLJhz3qI6r3H79u1wNze319V5DUIIyXf79m1jNzc325KOlTlLCWNsEIAlAOpDkWC/BnAVwH0AiQBSAegBMALQFEBrAJ0AdATwkjH2Ned8t1JeBSnRRI+JuPP6TsH83DPPz8Te3nthpm2m6tAIIYQQQgjKSLgZY5egSKDjAawCsJ1zfvttJ2SMNQcwHMAgADsYYxM45+2UEy4pSiQQYXmH5eh/uD8SshOQmJ2IqeemYmv3rVAXqqs6PEIIIYSQOq+sGm4HAFMBWHPOp5Yn2QYAzvktzvlkAA0AfJV3HlKNTLRMsLzjcgiZYjGAO6/vYPHVxXhbuRAhhBBCCKl+ZSbcnPNVnPMSl1B9G855Luf8ZwD2lQuNVERL85b4qsVXBdt/P/0bOx/uVGFEhBBCCCEEKCPh5pxnKOMCnPNMZZyHvN3QJkPxgcMHBdvLbyzHtehrKoyIEEIIIYTQPNzvEcYY5rWZB1djVwCAjMswLXAaXqW9UnFkhBBCCCF1V7kTbsZYPcaYB2OsXpH9FoyxbYyxm4yxg4wxN+WHScpLXaiOnzr9BGNNYwBAck4yJp2dhEwJfdFACCGEEKIKFRnhngXgBhQ3QwIAGGNqAC4C+BSAG4C+AM4yxmhpVhUy0zbDT51+glggBgA8TnqMOZfm0EqUhBBCCCEqUJGEuzOAsCKzlQwEYAcgEEB3AOsAGAAYr7QISaU0N22OOa3nFGyfijiFtTfXqjAiQgghhJC6qSIJd30AT4vs6w3F6pIjOecnOecTAIQB6KGk+EgV9GvYD0OaDCnY3nx3Mw4+OajCiAghhBBC6p6KJNyGUKw0WVgbAI85588L7buJQmUnRLWmt5iODvU7FGwvvLIQ16OvqzAiQgghhJC6pSIJdxYUS7gDABhjDaAY9b5UpF0OAFrisIYQCoRY1mEZnAydAABSLsXkc5MRlhKm4sgIIYQQQuqGiiTcDwG0LzRLyWAoyknOF2lXH0CsEmIjSqIt1sbaLmthomkCAEjLTcO40+OQlJ2k4sgIIYQQQt5/FUm4/wCgDeA6Y2wvgIUA0gEcym/AGFMH4AHgkTKDJFVnrm2ONV3WQFOkCQB4mfYSk85OQo4sR8WREUIIIYS83yqScP8CYCcUS7V/DCAXwBec85RCbfpAkZQHKi1CojTORs5Y4r0EDAwAcDPuJmZfmA2ZXKbiyAghhLxrVlZWrowxz8IPdXV1DwsLC9eePXvaHz16VEfVMdY0oaGh6owxT5FI5Flam/Xr19cTiUQeAoHAc+HChabvMj5Sc5U74eacyznnQwE4AGgLwIpzvrdIs+cA+gPYrrwQiTL5WPvgqxZfFWyfjDiJpTeWgnOuwqgIIYSoSvv27VP79euX0K9fv4T27dunAMC///5r2Lt3b6cFCxZQwlgBS5YsMRk/frwdAPbTTz+Fz5s3L66q53Rzc2vMGPM8ffq0dknHb968qcEY87Szs3Ou6rVI9RGVdoAx9gmAo5zztML7OedhUEz9VwznPARAiFIjJEr3v6b/Q0xGDP588CcAYNfDXTDVMsVI15EqjowQQsi7NnPmzJjevXsX/L8+JyeHjRgxosGOHTtMvvvuu/pDhw5NcnBwkKgyxtpg1qxZ5kuWLLESi8V806ZNz4cPH56s6phIzVHWCPdOAHGMsSOMsRGMMZN3FRSpXowxTG85Hd1suxXsWxWyCn8//VuFURFCCKkJ1NXV+YYNG15qa2vLJRIJO3z4sJ6qY6rpxowZY7VkyRIrDQ0N+e7du59Ssk2KKivhngngFhSL2GwCEMUYO8sYm8gYs34n0ZFqI2ACfN/+e7Qyb1Wwb/7l+Tj/quikM4QQQuoaHR0dbmtrmw0AsbGx4qLHzczMmjHGPJ89e1bsGAB4eno6McY8T5w4oVPa/sDAQC0fHx9HfX395hoaGh6NGzduunr1aqOSzgcAL1++FA0ePNja1NS0mbq6uoe1tbXLpEmTLDMzM1lp1wMAuVyODRs21Gvbtm1DAwOD5mKx2MPS0tJ18ODBNk+ePFGr+LvzH5lMhsGDB1tv2LDBXEdHR3bo0KHH/fr1Sy3aLjMzk+XXyJd2LmNjYzfGmOeLFy9EALB//349xpjnnTt3tAHA19e3ceF6+9OnT2v37t3b3sPDwxkAwsPDNQofL1xi8uLFC9H8+fPN2rVr19DKyspVXV3dQ1dXt7m7u3vj5cuXG8tkdC9XdSu1pIRzvhzAcsaYBYB+eQ9vAB0B/MQYuwngAICDnPOH7yJYolxqQjX83PlnDD8+HI+SHkHGZZgWOA1b/LagmUkzVYdHCCFEhdLS0oQAYGZmpvRykn/++Ud/48aNZg4ODtkdOnRIefXqlfqtW7e0J02aZJuSkiKcO3fuG7XPz549E3t7ezeOjo5WMzIykvr4+CTn5OQINm/ebHbp0iVdqVTKSrpOTk4O69Wrl/3p06cNNDQ05M7OzpkmJiaShw8fau7atcv42LFjhkePHn3Url27rIq+BolEgn79+tkdOXKkXr169aSHDx9+3LZt2wqfpzTW1taSfv36JZw5c0Y/OTlZ1LFjxxQjIyNp/nFzc3Opt7d3Wk5ODgsICDDQ0dGR+fn5JRc6XvBz279/v8GCBQvqW1hY5FpbW+e4u7tnxMXFiW/duqU9Y8YMm3PnzukdPXr0edEYiPKUmnDn45xHA1gHYB1jzBBAXyiSb18opgBczBh7DEXy/TfnPKga4yVKpqOmg198f8Gn/36KyPRIZEmzMO70OGzrvg0OBg6qDo8QQqrX8VmWuLreolxtXfxf4+OtEW/s2/+5DUIPGJerf+ux0ej+Q9Qb+7b1dkT4Bf1y9e+6KALtJhZd8blaBAUFaURGRqqLRCLep0+fYiO2VfXLL7+Yr1q1KnzChAkJ+ftWr15tNGnSJNvly5dbfvXVV/FaWloFd/OPHDnSJjo6Wq1z584phw4deq6rqysHgPDwcLGPj0+jsLAwjZKuM2HCBKvTp08btGrVKm337t1htra2BUnookWLTOfNm9dgyJAhDk+fPg0Vid6aEr2hW7dujmfPntW3sLDIPX78+ONmzZopdZ5dLy+vrAMHDoS7ubk1Tk5OFs2dOze6S5cuGYXbODs7x/v5+aV5eHgYGBsbSw4cOBBe0rm8vb3TAwMDH3To0CGz8P6wsDCxn59fo2PHjhn++eefBkOHDqVSmGpSkWkBwTlP4pxv45x/AMAEwCcA9gGwBDAbwDXGWARj7CfGWEfGWImfOEnNYqJlgg2+G2CobggASM5JxqiTo/Ay7aWKIyOEEPIuxcfHC/fu3av38ccfO8rlcixevPhlddww2bNnz6TCyTYATJw4McHGxiYnLS1NeOnSJa38/ffu3VM/d+6cvkgk4ps2bYrIT7YBwNbWVrJ48eJXJV0jKipKtH37dlMdHR3ZwYMHnxdOtgFg7ty5ce3bt0+NiIhQ/+uvvypUpy6TyXD27Fl9AFi3bl2EspNtZWvZsmV20WQbAOzs7CQLFix4BQD79+83fPeR1R0V+zhXCOc8A8BeAHsZY2oAukIx8t0HwCQAEwF8C2CxEuIk1cxW3xbruqzDyJMjkSnNRFxWHL44+QW2dd8Gc21zVYdHCCGkmvTp06dR0X1qamp83759T/z9/ZU+ug0APXv2LHEk1d7ePjsiIkL91atXagAyACAgIEAHADw9PdMdHR2LJf+ffPJJyqhRo+QZGRlvDCIePXpUNzc3l3l7e6dZWFhIi/YDgPbt26ddvHhR7/LlyzoDBgwo92sVCARo3rx5ekhIiM6oUaNs7e3tH7m5udXopDsnJ4cdOnRI9+rVqzoxMTGi3NxcAeccKSkpQgB4/vy5uqpjfJ9VOuEujHOeC+AogKOMMQEUdd4fAajy/JPk3XE1ccXaLmsxJmAMcmQ5iEyPLEi6jTRLvY+FEEJqr+4/RBUr86iIj7dGFCszqYjhR55Wuq+StG/fPtXU1FTCOUdcXJw4KChINycnh40aNcrOycnpoYuLi9ITSTs7u9yS9uvq6soAICsrq+Ab8sjISDUAqF+/fol9BAIBLCwscp8+ffpGWUl+Ann69GkDxlipC9UAQHx8fIXyIcYYzpw588THx6dhSEiIjq+vr1NAQECNTbqDgoI0Pv74Y8eIiIhSk+r09HThu4yprlFKwl0Y51wO4Gzeg9QyLc1b4qdOP2Hi2YmQyqUITw3H6FOj8Wu3X6GvXr4yQ0IIIbVH0Xm4IyIixF27dm345MkTzcGDB9vdunXroUBQoQpUyOXyMktKK3o+QJHklnGs2OptMpmMAYCdnV22u7t7RvFe//Hy8irzeEn09fXlhZPurl27Op06darSSbdcLn97o0qQyWTIT7a7deuWNG3atFg3N7dsQ0NDmUgkwrVr1zRbt27dlBbAq16VSrgZY+ZQ1G2XeJMCAHDOL1c2KKJa3vW9sdR7Kaafnw45l+NR0iOMPT0Wm7pugra4xIWuCCGEvCdsbGwke/bsee7l5dX07t272hs2bKg3duzYxMJtxGIxB4DU1FQhgGJlHvmj0spgaWmZCwB5ZSYlio6OLnasQYMGuQDg7OycWdrNhFVV3qRbXV2dA4BEImHZ2dlMQ0Pjjew2LS1NkJycrPRBUAC4du2aZkREhLq5uXnu0aNHnwuFbw5kP3z4kEpJ3oEKfcRkjPVnjD0AEAngBoALpTxoMudazs/WDwvbLizYvhN/BxPPTES2NFuFURFCCHkX3N3dsz/99NN4AFiyZImlRPJmTm1mZpYLAKGhocUG3q5cuaIZHx9f4vzcleHr65sOAMHBwTolzfu9e/du/ZLKIfr06ZMqFAr5hQsX9BMTEys+pF5O+Um3h4dHenx8vLhr165Od+7ceSOJFQqFMDY2lnDOcffu3WIJ7oEDB/RKG2HO/3BT2tSH6urqcuC/Ef2iXr9+LQIU0zsWTbYBYOfOnVQz+g6U+xcwb6n33QCcAKQCuAPgcimPK0qPlLxzfR37Ynar2QXb12OuY8KZCZR0E0JIHbB48eJobW1t+cuXL9XXr1//RlLWsWPHNAD48ccfzZOSkgpyicePH6t9/vnndsqMw8XFJcfb2ztVIpGwUaNGWaenpxcklhEREeJvvvmmfkn9bG1tJUOGDIlPSUkRdu/evWHRJBgAYmNjhStWrDCOioqq0uhy0aTb19e3WNLdtm3bNACYN2+eZU5OTsFruHLliuasWbMalHZuCwuLUj/cAECDBg2kAoEAsbGxaoV/FvlcXV2zGWO4f/++1unTp9/4mnrZsmUmJ0+eNKjYqyWVUZFPfPmZ1yQAJpxzd865d2mPaoiVqMCgxoMw2WNywfbV6KuUdBNCSB1gaWkp/fLLL2MAYMWKFRaFR7lnzJgRa2pqKrl9+7a2k5OTi5+fn0OrVq0aubm5ORsaGkqaNWtW4Zrosvz6668RZmZmkjNnzhjY2dm59uzZ097Hx8exadOmLgYGBlIXF5dMAFBTU3ujEHrjxo2vunXrlnTjxg0dDw8PZxcXlyY9e/a079y5s6OTk1PT+vXru02fPt0mMTGxyjcMvi3pnj9/frSWlpb8+PHjho6Ojs49evSw9/DwaNyhQ4cmHTp0SC28qE1hffv2TQaAb775xtrX19dh4MCBNgMHDrS5f/++GgDo6urK27dvn5Kbm8tcXV2d+/btazdw4ECbSZMmWQKKqf8GDBjwWiKRsG7dujVu27Ztoz59+tg5ODg4f/3119bjxo2LqeprJ29XkYS7IYDLnPM1nPMSfynI+2mE6wiMbz6+YPtq9FUqLyGEkDpg7ty5sUZGRtJXr16pr127tmCBH3Nzc9mFCxce9u7dO1EikbBz587px8bGqo0fPz7m7NmzT0UikVLvwGvYsGHu9evX7w8aNOg1YwwBAQEGT5480Rg+fHhcYGDg48TExPyyiTfyEw0NDX78+PHnf/zxx7OOHTumxMbGik+dOmVw69Ytbc45Pvzww8Q//vjjmZOTk1JmFykr6XZ3d88OCAh42KlTp5Tk5GTR2bNnDTIyMgSLFy9+uWvXrlJnuhk5cmTSokWLXlpbW2dfvHhRf+/evcZ79+41jo6OLiiv2blzZ7i/v39Cbm4uO3r0qOHevXuN//nnn4J5tf/888+IH3744YWjo2PWzZs3tc+fP69vYWGR+9dffz0eMWJEQslXJsrEyntXKmMsEkAg53xw9YZUu7Ro0YIHBdWNxTU33t6ItbfWFmy3sWiD1T6roSEq9d5ZQggpEWMsmHPeojqvcfv27XA3N7d3sjIjUZ179+6pu7q6uujq6sqSkpJuVWYGFEKU4fbt28Zubm62JR2ryG/lSQAtlRIRqZVGu43GuObjCravRF/BpLOTaKSbEEJItZLJZLh48aJW0f1PnjxRGzp0qB3nHB9//HECJdukpqrITQLfArjBGFsKYDbnXFZNMZEa7Eu3LwEA626tAwBcjrqMSWcnYVXnVTTSTQghpFrk5OQwb2/vJpaWlrn29vbZ+vr6sqioKLX79+9r5eTksEaNGmUtX748UtVxElKacifcnPMXjLH2AA4B6McYOw3gFYASZ2rnnH+vnBBJTfOl25fg4Fh/az0ARdI97vQ4rPZZTfN0E0IIUTo1NTU+fvz4mMDAQN179+5ppaWlCdXU1Lijo2NWnz59kmbNmhWnp6dXPSvHEKIE5U64mWKJp/FQ3DwpBOAAoKQCcJa3nxLu99gYtzEAUJB0X4+5jlGnRuEX31+gp6anytAIIYS8Z0QiEdasWUMj2KTWqkhJydcAJgCQAjgK4CmA9OoIitQOY9zGQE2ghp9DfgagWBxnxIkR2Nh1I+pp1FNxdIQQQgghNUNFEu4RADIBtOec36qmeEgtM8J1BDRFmvjh+g8AgIeJD/HZ8c+w2W8zTLVMVRwdIYQQQojqVeR2XisA5ynZJkUNbjIYC9suhIApfp2epzzHsH+HITKdvv0jhBBCCKlIwh0JxQg3IcV81PAjLPVeChFTfGnyKv0Vhv07DGEpYSqOjBBCCCFEtSqScO8B0JExRtNQkBJ1t+uOnzr/BLFAsfhVbGYshh8fjnuv76k4MkIIIYQQ1dP6DoEAACAASURBVKlIwr0IwGMA/zDGHKopHlLLdWrQCeu6rIOmSBMAkJidiM9PfI7LUZdVHBkhhBBCiGpUJOH+B4oZSjoDeMAYe8gYC2CMnSzhcaJ6wiW1QRvLNtjUdVPB9ICZ0kyMOz0O/4b9q+LICCGEEELevYrMUuJbpF+jvEdJSpqfm9QhzU2b4/cev2P0qdGIzYyFVC7FjPMzkJidiCFNhqg6PEIIIYSQd6YiCXfXaouCvJccDBzwZ88/MfrUaDxPeQ4AWHJ9CV5nvcZE94lQrKVECCGEEPJ+q8jS7qerMxDyfjLXNsfvPX7HuNPjcDv+NgBgy90tSMhKwLw28yASVOQzHyGEEEJI7VORGm5CKkVfXR+b/TajQ/0OBfsOPj2ICWcmIEOSocLICCGEAIBMJoOFhYUrY8yzXr16bjk5OTXmK8jVq1cbMcY8/f39bd/F9aZOnWrJGPOcOnWq5bu4XlH+/v62jDHP1atXG1Wkn1QqxY8//mjcunXrRoaGhm4ikcjD0NDQzd7e3rlHjx72ixYtMo2KilLKKJeVlZUrY8zz0aNHaso4X3kwxjwZY57v6nrKRgk3eSc0RZr4ufPP6OvQt2DfxciLGPbvMMRkxKgwMkIIIQcPHtSLiYlRA4CkpCTR7t279VUdU3V49OiRGmPM08rKylXVsShTUlKSoFWrVk7Tpk2zCQoK0rW1tc3p1q1bcuvWrdPEYjE/efKk4bx58xoEBgbWyKmdK/shozYp9ZMOY+w8gK8555Wez40x1g7AD5zzDm9tTN57YoEYi9otgqmWKTbf3QwAeJT0CEOODsE633VoXK+xiiMkhJC6aevWrcYAYGpqKomLixNv27bNeNiwYcmqjksVpk+fHvfpp58mmpubS1UdS3nNmDHDMiQkRMfR0TH76NGjTxo1apRb+HhkZKRo69at9SwtLSWqirGqQkJCavWiHmWNcDcGcIExdooxNpAxpl6eEzLG1BljgxhjAQDOo/SZTEgdxBjDRI+JWNh2YcGqlHFZcRj27zCcf3VexdERQkjdExsbKzx9+rQBYwy///77c6FQiAsXLuiHh4eLVR2bKlhYWEjd3d2zLSwsak3C/c8//9QDgKVLl74smmwDgJWVlXTu3LlxHTt2rLUrhru7u2e7u7tnqzqOyior4W4IYA2AjgB2AohljB1ljM1hjPkzxjoxxjzynv0ZY3MZY8cAxAH4E4A3gFWghJuU4KOGH2G973roiHUAKObqnnBmAvY83KPiyAghpG7ZvHmzUW5uLvPy8krr1q1bert27VJkMhk2bdpU6tf7hetpN2/ebNi8efPGWlpa7tra2u5t2rRpdOLECZ2S+p05c0Z79OjR9V1cXJoYGRm5icViD1NT02bdu3e3P336dLnLHdauXWvEGPP09vZuWFqb69evazLGPE1NTZtJJBL4+/vbNm7c2BUAoqKi1PJfQ9ESk7fVcIeEhGgMGjTIxtra2kVDQ8NDT0+veaNGjZqOGjWq/uPHj9+oad62bZtB//79bR0dHZ11dXWbq6ure1hbW7t8+umn1k+fPlXaB5rExEQRAJibm1d4BFsul2PdunX1vLy8nPT09Jqrq6t7NGjQoFIxvq2228vLy4kx5nnkyBFd4L8Sn7/++ssIACZNmmRb+OdSuMSkrBru6Oho0ZgxY6zs7OycNTQ0PHR0dNzd3NwaL1myxEQiKf6WFL4vICkpSTB69Oj6VlZWrmpqah6mpqbNhgwZYh0bGyusyGt/m1ITbs55Cud8MhQj3asByAH0ALAAwF4ApwHcyHveC2A+gO4AJAB+BODEOZ/KOU9VZsDk/dHGsg3+6PEHLLUVf9PkXI7F1xZjxY0VkHO5iqMjhJC6YceOHcYAMHTo0AQAGDZsWAIA7Ny50/htfSdPnmz55Zdf2ovFYt65c+cUMzOz3KtXr+r26dOnUUBAQLEEes6cOVa//vqrmUQiYW5ubhldunRJNjAwkJ44ccKwW7dujbdu3WpYnphHjhyZWK9ePemlS5f0QkNDS/wG/ueffzYBgE8//TReLBajXbt26d26dUsCAE1NTXm/fv0S8h+9evVKKs91165da9S6deumu3fvNuaco3PnzsleXl5pnHO2efNms+PHj+sWidPh6NGjhpqamvJ27dqltmvXLjU3N1fw559/mrRs2bLpnTt3ylU98DYWFha5ea/ZVCaTlbufXC7Hhx9+aDd+/Hi7mzdvaru6umZ07do1mXPO/vzzT5MWLVo4BwYGaikjxpLo6enJ+/Xrl9CgQYMcAPDw8Egv/HNxcnLKeds5QkND1T08PJps2LDBPD09Xejj45PcsmXLtMePH2vOmjXLumPHjg2zsrJKvAk4NTVV2KpVq8Z79uwxbtq0aWb79u1Ts7OzBTt37jTx8fFppMybh996tyrn/DmAKYyx2VCMdncC0ByAKQB9AMlQjGqHADgL4ALn/K1vECEA4GjoiB29dmD86fG4l6Aoz9p+fzsiUiOwpMMSaItr5P0dhJD3gOt211o748HdYXeDlXGeS5cuaT58+FBTW1tbPmzYsCQAGDx4cPJXX30ljYiIUD9x4oROt27d0kvrv23bNtNz58498Pb2zgQUs50MHTrUZvfu3cbz5s2z9PX1fVK4/dSpU2NatWr1vEGDBm+Ua+zcuVN/2LBhDlOnTrXp379/iq6ubpmjLhoaGnzo0KHxq1evtli9erXJpk2bXhU+npiYKDh06JCRUCjkEyZMeJ137de9evVKbdy4saGhoaH0wIED4RV5rwIDA7UmT55sA4CtXLkyYtKkSa8Fgv/GLUNCQjSK9tmwYcPzgQMHvvF6JBIJpk2bZrl69WqL8ePHW58/f/5J0X4V9fnnn8cvWLCg/r59+4wvXbqk17Vr12QvL6+MVq1aZbq7u2cXjrOwZcuWmRw+fLiekZGR9Pjx449atGiRDShmPBk5cmSD7du3mw4ePNjh6dOnoZqamkpf1NDCwkJ64MCBcH9/f9uXL1+qDxs27PXEiRMTKnKOQYMG2cfExKj16NEjaf/+/WFaWlocAJ4+fSr29fV1unLlit60adMs161bF1m0b0BAgEHHjh1Tbty48VBfX18OAOHh4eI2bdo0vn//vtbWrVsNx4wZk6iM11ruWUo451mc8+Oc868559055x6ccwfOuSfnvAfn/BvOeQAl26SijDWNsbXbVnRu0Llg37lX5zD02FC8THupwsgIIeT9tnHjRhMA6NWrV2J+UqihocH79u2bCABbtmwpc5R7xowZkfnJNgAIhUKsWLEiEgCCg4N1i44Qfvzxx6lFk20AGDx4cEqPHj2SUlJShEePHtUterwkU6ZMiRcKhXzv3r3GmZmZb1znl19+Mc7MzBT4+fkl29raKuVGwUWLFlnIZDI2atSomClTprwumsR6eHhke3h4vFFjPHLkyKSiHx7EYjFWrVoVZWJiIrl06ZJeUlJSlWeMmzdvXuy0adOiNDQ05FFRUWrbt283HTdunF2LFi2cjYyM3P73v/9Zh4WFFSsPWbdunRkAzJo1KzI/2QYAkUiEDRs2vDI3N8+NiopS27ZtW7m+eXjXjh8/rhMaGqqlra0t/+233yLyk20AcHR0lCxbtuwFAGzfvt206O8IAGhpacl///338PxkGwBsbW0lI0aMiAOAM2fO6CkrVpoWkNQIWmIt/NTpJwx3Hl6w72nyUww6OgjXo6+rLjBCCHlPZWVlsfyb7UaMGPG68LEvvvjiNQAcO3bMMCUlpdRcwd/fP6XoPisrK6menp4sNzeXlVQHGx0dLVq9erXRqFGj6g8cONDG39/f1t/f3/bRo0eaAPDo0aNylVnY2tpKunXrlpySkiLcsmVLvcLHtm7dagIA48aNiyvPud5GKpXi8uXLegAwduzY129rX9idO3fUFy9ebDp8+PAG/fv3t81/vTKZjMnlcty/f7/KZSUCgQDLly+PDg8Pv7Ny5cqIjz76KMHR0TGbMYbk5GTRH3/8YeLu7u584cKFgvKQZ8+eiV+9eqUuEAgwZsyYYqPKGhoavF+/fokAEBgYWK4PQe/amTNndAHAx8cn2czMrFgtzYABA1JNTEwkGRkZgosXLxYrjXF2ds60trYu9gGwSZMm2QAQExOjtDp7WuaP1BhCgRBftfgKjgaOWHBlASRyCVJyUjDq1Ch87fU1Pmn8iapDJIS8R5RVllFb/fHHHwYpKSlCGxubHD8/vzdWIWvXrl1W48aNsx4+fKj522+/GU6ePLnEr/kdHR2LzYgBADo6OrLU1FRhVlbWG8n68uXLjefNm9cgOzu71CQ+NTW13DerTZo0KfbYsWOGmzdvNs0vRTh8+LDu8+fPNRwdHbN79epVajlMRURHR4uysrIEQqGQu7i4lOubfIlEgv/97382e/bsMea89GqM5ORkpd2cZ2ZmJpsyZcrrKVOmvAaAqKgo0a+//lpv2bJllikpKcLPPvvM7unTp/cAICIiQg0AjI2NJYVHhgtzcHDIAYDo6OgaOWNNZGSkGABsbW1L/Zk0aNAgJz4+XvzixQs1AG/8nltZWZXYL3/EOycnR2kD0zTCTWqcvo598Vv332CsqfgmU8Zl+O7ad1h0ZREk8lo7hSghhNQo27dvNwaAtLQ0oaenp1PRR0JCgggA/vjjj1LLSoTC8ueK58+f15o5c6aNVCplc+fOfXXr1q3QlJSUmzKZLJhzHjxu3LgYAOCcl/tGNT8/v4wmTZpkhoaGap0/f14LANatW2cCAPllAaqyePFis927dxsbGxtLNm3a9PzJkyd3MjMzQzjnwZzz4ObNm2cAFXu9FWVpaSmdO3du3Lp168IB4NmzZxp3795Vz7suAMV0vaUp64NCZcjlyp0QoZyvodSDpdW2Vwca4VaRbZfCIJTloEcjLRibW6s6nBrHzcQNu3rtwqSzk3A/4T4AYO/jvXie8hwrO62EoUaNLCcjhJBa4enTp+KrV6/qAYop5RITE0ucxg8AQkJCdO7cuaPerFmzKt2jtXv3bkPOOT777LO4hQsXxhY9/vz580qVVowePTpu8uTJtmvWrDG1traODAgIMNDW1paPHj26QjfflcXCwkKqoaEhz87OFty7d0/d2dn5re/F33//bQgAq1atihg0aFCx0puIiAilzFBSHh9++GHBjHExMTEiV1fXHFtb21wAiI+PF2dlZbGSbooMCwtTBwALC4tyjXaJxWIOAKmpqSVmslFRUUpdCr5+/foS4L84S/Lq1Ss1ALC2ti7x25h3hUa4VUAqk2Pt2We4fHwnDH9phrs/dMaNv9cgPbVcsxLVGeba5tjWfRt62PYo2BcUG4RPjnxSMKMJIYSQituwYYOxXC5HmzZt0vJHXEt69OjRIym/fVWvmZSUJAKABg0aFEt8oqKiRBcvXqzUDWpffPFFooGBgfTIkSOG8+fPN5fJZKxfv34JhoaGxYZT1dXVOQBIpdIKjSqLRCK0bds2FQDWr19frvciJSVFBAD5iW1hBw8e1Mt/P5ThbSPHT58+LUh0bWxsJADg4OAgqV+/fo5cLseGDRuKzbmek5PDDh48WA8AOnbsmFaeOMzMzHIBIDQ0VLPosRs3bmjExMSUmHCrqalV6ufi4+OTBgBnzpwxiI+PL/Z1y4EDB/Ti4+PFWlpa8vbt26t00R9KuFXg6vNEvE7PwQfCyxAyDtecELS8NQfCHxsh6Md+uH12H6QlTNReF2mKNLG0w1JM8pgEBsV/h1EZUfjfsf/hryd/qTg6QgipfeRyOfbs2WMEAIMGDSpzFDh/bu79+/cbSaVVW3jRyckpGwB2795tVPhGzKSkJMHQoUNt09LSKlXLrKWlxQcPHvw6OztbsH37dlMAmDRpUonlJBYWFlKxWMwTEhJEJSVoZZkzZ060UCjExo0bzQovyJLv5s2bGjdv3iyYGtDe3j4bANasWWNSeG7se/fuqU+cOFGpX227ubk1Wb58ufHr16+LvaaHDx+qjRw50javXUbhlSjHjh0bCwA//PCDZeHYpVIpxo4dWz8qKkrN0tIyd/jw4eUaEcxPzFeuXGmemJhY8DN++vSp+LPPPrMrrUTF0tIyFwAePHhQbGrFsnTv3j3dxcUlMyMjQ/D5559bF55vOywsTDx9+vQGADB8+PC40urU3xVKuFXAxUoP333oDAvNN2+o1WS5aJF2GqbnZqDtkjNYePg+QiNTlF5DVdswxjDSdSRW+6yGrlhxo3SuPBffXv4W8y/PR46MZqIkhJDyOnLkiO6rV6/UNTQ05EOHDi0zkfL39081MDCQxsfHi/ft26dfleuOGzfutbm5ee79+/e17OzsXP38/By6du3qYGdn1+zu3bta/fv3r9DsH4VNmTIlLr+e3MvLK83T07PEJcDV1dV5p06dUmQyGWvevHnTDz74wG7gwIE2Y8eOtXrbNTp37py5YsWKcECxIqK1tbVLr1697H19fR0aNWrU1MPDw/nChQsFi0fMnj07WiQS8V27dpk4ODi49O7d275du3YN3d3dnS0sLCTu7u4ZpV6sgsLDw9VnzJhhY2Fh4ebi4tKkR48e9j179rRv3rx5Y2dnZ9ebN29qm5ub5/7+++9hhfvNnDkzvnfv3onx8fHiVq1aNfX29m7Yp08fOzs7O5etW7ea6unpyXbu3PmsvHNwT58+Pc7c3Dw3NDRUy8nJycXPz8+hdevWjZo1a+aio6MjK+01+/v7JwsEAmzdutWsffv2DQcMGGAzcOBAm1OnTr11MY5du3Y9NzMzkxw5cqSejY2Na69evex9fHwcnZ2dXcLCwjTatGmTtmLFiqjyvZPVhxJuFTDQUsOQ1rZoPvssYkYE46rDJIQJbAuOH5K1Q1yGFFsvhaH3movw++k89hw5gdhXT1UXdA3QqUEn7Oq9C44GjgX7Djw5gGH/DkN0erQKIyOEkNpj69atxgDg6+ubXFLZRWHq6ur8gw8+SASA3377rUplJSYmJrKgoKAHgwYNeq2lpSU/d+6c/t27d7W7d++eFBQU9CC/HrcyHB0dJXZ2dtkA8OWXX8aX1fb3338PHzBgwGuZTMaOHTtmuHfvXuNDhw7VK6tPvsmTJydcvnz5gb+/f4JUKmWnTp0yuHHjhq5AIMDo0aNje/ToUVB64evrmxEYGPigU6dOKenp6cKAgACDmJgYtYkTJ0YHBgY+FolEShtNO3369KNvv/32Vfv27VPzpsDTO3nypEFERIS6h4dH+pw5c17dv3//XtE6fIFAgEOHDoWtXbs2rFmzZhm3bt3SPnHihKFcLmdDhgyJDw4OvtexY8dyl2KYmJjILly48PCDDz5IlEql7Ny5c/oxMTFqY8aMiTl37tyT0l5z27Zts7Zs2fLcxcUl4+bNmzr79u0z3rt3r3F5RrxdXFxybt68eX/06NGxWlpa8oCAAINr167pOjo6Zn3//fcvzp49+6Q6Fu2pKFbe0VPG2EgAOzjnWdUbUu3SokULHhQUpJRzPbt7FfEXt2FlQitcTzd949jv4h/QXhCK+xpuyGz8MZp2GQodvbp542CmJBMLrizAsbBjBfsM1A2wrMMytLFso8LICCHlxRgL5py3qM5r3L59O9zNza3So6ak9rhy5Ypm27Ztm5qYmEgiIyPviMU1chY78p67ffu2sZubm21Jxyoywr0JwCvG2I+MsYZKiYy8wcG1NVqP2YCds4Zh++de+LC5JTTFQpggCe0EoRAwDpecW/C6PQeCH51w4+dPcP/qcXAlT7NT02mJtbDEewm+9voaIqa45yQ5JxlfBnyJLXe3QM7r1vtBCCF13Zw5cywB4IsvvoijZJvURBUZ4f4HQA8AQgByAAEA1gE4wutwkbEyR7hLkpEjxaUrF2FxZT6cs29BwIq/1S+ZJV7Z9oNj11EwsbSptlhqopDYEHwV+BVeZ/03iNWhfgcsbreYpg4kpAajEW5SVTt27NA/dOiQwYMHD7RCQ0O1LC0tcx88eHBPT0+PRl2ISihlhJtz/gEAewBLALwG4AfgbwBhjLGvGWMmlQ2QMSZmjHXJGz2/yhiLZozlMsYiGWP7GWOd3tJ/MGPsAmMshTGWzhgLYoyNY4yV+foq2+9d0lYXwa9TJ7jOOofXo27iqsMkhAsavNGmAY9Cm7C1MNzYHEHL+uD43WjkSuvG3xsPMw/s7b0XHqYeBfvOvzqP/of7IyQ2RIWREUIIqU7BwcHa+/btMw4LC9Pw9vZOPXbs2GNKtklNVe4R7jc6MSYGMADAWABtAHAAuQD2A1jPOb9SwfP5AjiVtxkDIBiK5TebAnDJ27+Icz6vhL7r8uLIBnAagARAFwC6AA4C6M85lymrX1HVPcJdEi6X4/HNQCRf+g1NE05Cl/1XVv+3rC0mS8bDSFsNH7lbYUDLBmhkpvtO41MFiVyC1SGrse3etoJ9QibEePfx+NzlcwhqzmcoQghohJsQ8v4pa4S7Ugn3GydgzA3AOACDAGjl7b4FRbnJDs75W+dsY4z5QJH8ruKcXyhybCCAHVCUsvhwzs8WOuYPRZIfA6AD5/xJ3n4zAGcBNAEwmXO+qsg5K9WvJKpIuAvLykhDaMAf0Lq3C865dzA4dzYuy13eaDPfKABNbOujSdfh0DMoNnXoe+X8q/OYfXE2UnL+W9SrnVU7fN/+e9TTKNdN6ISQd4ASbkLI+6ZaE24AYIxZAZgJYHyh3RxAHIA5nPNfq3j+LQBGANjKOR9RaH8QAE8Awzjnvxfp0xHAOSiSaivO/7uTrrL9SqLqhLuwyOf3sesxw/6QKMSkKqYg1UQ2rquPgy7LQiZXR6hhFxh0GI2GzTuACd7PUd+YjBhMD5yOW/G3CvaZappiaYelaGFerf9/J4SUEyXchJD3jbJmKSmGMebLGPsLQBgUo9zZALZCMdp9DIApgE2MsYlVuQ6Am3nP9Qtduz4USXMugH1FO3DOAwFEAjAH0Lqq/WoDK/ummNa9CS597YNtn7VEL1cL9BbdKCg50WI58Eo+hkb/9MXz7zxxbe9ypKUkqjhq5TPXNsfW7lvxucvnBfvisuIw4uQIbLy9ETL5WyuFCCGEEEKUpsIJN2NMnzE2mTH2EMAJAB8CiAIwG0B9zvlIzvkeznkfAG2hqMWuasKdPw1h4dVN3POe75UxN/iNIm2r0q/WEAoYOjmZYt0QD8yePBlXnWYgTPDm7CUOsudodX8xhCsb49qqoXhy87yKoq0eYoEYUzynYH2X9TBUV8xWIudyrL21FiNOjqCFcgipI+rwJFqEkHco729NqX9wyp1wM8Y88ko7IgH8CKARgEAA/gDsOedLOedvDJdyzq8BOArAuuKhF1zXHMDwvM0DhQ7Z5T1HlNH9RZG2VelXKxkam6P1oG9gO+cWHvY6gBv63ZHN/5ujVIvloFXSYTQ81Acnf/gYO65FID1HqsKIlcu7vjf29dn3xiwmwbHB8P/HH8fDjqswMkJIdWOMJeXm5tKkzISQapedna3GGCu1hK0iI9xBAPK/o98CoBnn3IdzfvAtdc4ZAEQVuE4BxpgIwJ8A9AGc5pwfLnRYp9D5S5Oe91x4mo7K9qvVmECAxi190XLKHuRMeoCrTjMRLnjzc9DJdHt8czAUrb4LwKy/7iI0MqWUs9UuZtpm+LXbrxjrNrZgtpI0SRqmn5+Oby5+g/Tc9LecgRBSG8nl8n+Tk5Pfm7/jhJCaiXOO+Ph4XZlM9kdpbSqScIcDmA5F2choznloOft9AaCyIwwboJiq7yWAoUWOsbznin5fWNl+/52AsVF5c3YHxcfHV/Y0KqNfzwStB82GzZzbeNhzH27o+yGOG+CITFGynpErw67rL9B7zQVc+r4nru9fiYy0ZBVHXTUigQhjmo/B9u7bYaVjVbD/n2f/oP/h/rgdf1uF0RFCqoNMJtsUGxubHBsbWy8nJ0dM5SWEEGXhnEMqlQpTUlJ0wsPD6yUlJd2Ry+W/l9a+IiPPDpVZUTKvT4XvUmOMrYJiZpIYAF045zFFmqTlPeugdPnH0grtq2y/ApzzTVAsdY8WLVrU2r/gTCBAYy8/wMsPyWnpmHE7Hjuvv8DTOMWIbxvBfbTLvQSEXkL63WW4ZtID5l3Gw6aJp4ojr7zmps2xv89+/HD9B/zz7B8AwKv0Vxj27zCMdhuNL1y/gEhQqS9kCCE1jKenZ3hwcHC/6OjoUbGxsT0458aqjokQ8v5gjGUCuCWRSI4B2O3p6ZlbatsKLO1+EsBxzvnKt7SbAqAH59yvAjEXPcePAKYCiAfQiXN+v4Q2HwA4BOAm59yj6PG8Nn8B+AjABM752qr0K01NmhZQGTjnCIpIws5rL+Bzbxb6CC4Xa3NPzQ25niPg6jMIIrGaCqJUjn/D/sWiK4uQJvnvc1Vzk+b4vv33aKDXoIyehJCqehfTAhJCSE1RkYRbDmAb5/zzt7TbDOBzzrmwUgExtgyK0pUEKEa2S/yunzHWAIqbG3MBGJQ04whj7CUUUwm255xfqkq/0rxvCXdhya9j8PDEJlg82wMb+atix2NhhOe2A9Cox3gYmdUv4Qw1X1R6FGZdmIWQuP+WgdcUaeIrz68wwGkAGGNl9CaEVBYl3ISQuqQ6Vj5RA1DmYjGlYYwtgSLZTgLQtbRkGwA45y8BhORdr38J5+oIRdIcA+BKVfvVRQbG5mg9ZB6s59xFaNc/EaLdAVL+36+MGRLQJvwX6K53w/bNKxHyIqnWTcFlqWOJrd22YqL7RIiYopQkS5qFxdcWY/Sp0YjJKFrJRAghhBBSMUpNuJliONATQIVX9mKMLYJitcpkKJLtm2/pAgA/5D0vZYw5FjqXKYD1eZtLSphFpbL96iQmEMClXR94TD+MhC9u4Gr9z5AA/YLjQsiw8ZkR+q2/jA/WXsK+oJfIltSexWWEAiG+aPYF/uz1JxwNCn4dcCX6Cj469BH+fvp3rfsgQQghhJCao8ySkry67Xy+UCxwU6yepHtzGgAAIABJREFUOo8IigVqLAHs55wPLHcQ/9VVA4rpB++V0vQh53xJkb7rAYyBYpXLAAASKGY20QPwN4CPOefFsr/K9ivqfS4pKUtOdibunvwdund+Q3iODkZLpr5xvKlmMuaaX4ZNt4mwtHVSUZQVlyPLwbpb67AtdBt4oYlsOtXvhG/bfgtjTbrnihBloJISQkhd8raEu/AIL8d/U+qV5Q6AvpzzshaWKXqd4QB+K0fTQM55pxL6D4ZiaXlXAEIAD6FYYv6XskapK9uvsLqacBcWGh6D7Tdiceh2FHKlirdtpmgXxogOQ84Z7mi3hrDVKDi37wuBsFKl/e/czbibmHNxDl6kvSjYp6+ujzmt56C7bXcVRkbI+4ESbkJIXfK2hLtL/j8BnIRiKfcVpTTPBRDJOX+u1AhrOEq4/5OUkYu9QS+x58oT7M8agXrszQVlXjJLRDYaiqY9voSegZGKoiy/TEkmfg75Gbse7npjf1ebrpjdajaNdhNSBZRwE0LqkorMUnIBwNGiJR11HSXcxcmkUtw9tw+CoC1oll38vcnk6rhr3AOmXcbDrmlLFURYMVejr2LepXmIzogu2KenpoeZXjPRx74PzWRCSCVQwk0IqUvKnXCTklHCXbaXT24j8tRaOMcehi4rNgMj7qu5IrLbZnRq7gSxsDomzVGO9Nx0LLuxDAefHnxjfzvLdpjXZh4sdSxVFBkhtRMl3ISQuoQS7iqihLt8MtKSce/4Zpg8+AN28v/K+x/J66Nb7lKY62liSCtrfOJlDRNddRVGWrbLUZex8MpCRKZHFuzTFGlissdkfNL4EwhYzf3QQEhNQgk3IaQuKTXhZozNzvvnL5zzpELb5cI5/76qwdUGlHBXDJfLcf/qceRc/gXN0i7iW+lw7JD5FhwXCxkm2segazNrOHn6gAlqXgKbKcnEmptrsOPBjjdmMnE3dcf8tvNhr2+vwugIqR0o4SaE1CVlJdxyKGYmacI5f1xo+63nBMAru9JkbUMJd+XFvnqGPXfT8Hvwa7xOz8nby3FUbTacBRF4KnRAostwNOv2OTS0dFQaa0luxd3C/Mvz8SzlWcE+sUCML//f3n2HV1Xl+x9/f9MJobfQey/SERsgoGBDxTKoI7bREUenODN37v3dmTvlzujcGa96R4cpomJHAXtBRSyMgCAiNaF3QgslAULa9/fHOWBCioSck5zkfF7Pc57t2WudtddeDwkfN+usddb3ubXPrcTHxFdj70QimwK3iEST8gL3fxMI2A+7e2aR96fF3X8Zmi5GNgXuysvNL+Tdlbt4ZsEWfOsiZif+ulj5QVJIS51Au4vvo1XHHtXTyTLkFuTyzxX/5InlT5Dv+SfPd2nYhV8N/xUDmg+oxt6JRC4FbhGJJprDXUkK3KG1dvUyDr3/P/Q98D5JllesrNCN5cnDiBl2J33OvzKi1vROz0znvz7/L1btL75n08SuE/nRwB/RMKlhNfVMJDIpcItINFHgriQF7vA4tD+DNe9Opd2GF2nlu0uUb7NWbO5xB/2uuI8GdSJj6kZ+YT7Pr3mex5c9zrH8b1ZkaZTYiPsH388Vna/QEoIiQQrcIhJNFLgrSYE7vAry81n56Uzsiyfol7O4WNkT+eN5yG7hqoGtuXl4e3qk1q+mXha3K3sXD3zxAPO2zSt2fnCLwfxy+C/1pUoRFLhFJLpUZOObu4FHgavc/e0y6lwGzAamuPsTIetlBFPgrjrb1q9gxweP0Wv3G9TnKCOO/y9bPPVk+dCOjflRt0yGnDuG+ITqX1rwo60f8cAXD5BxJOPkubiYOG7rcxvf6/s9kuKSqrF3ItVLgVtEoklFAvcHQF+glbsXllEnFtgBLHP3cSHrZQRT4K56R7MP8cW8N3lwQzvSMrJOnm9BJvMTf8gBa8CGdtfS9ZIf0DS1XTX2NLCE4NSvp/Ls6mcp8IKT59uktOE/hv0H57c5vxp7J1J9FLhFJJpUJHBvA9Lcfey31PsA6O7u1Zt0qogCd/VxdxZvPsD0BZuZszKDe2Ne4Ydxs0+W53osK+qPoO75d9N98JhqXdM7PTOd3y38HV/v/brY+ZFtR/LzIT+nbb221dQzkeqhwC0i0aQigTsHmOnuN31LveeBie4eFf9ersAdGXYfzmHNq3+kz6YnacrBEuUbYjuR2etm+oy7gzp161VDD6HQC5m1bhYPf/kwWbnfPJlPiEnglj63cEffO6gTV6da+iZS1RS4RSSaVCRw7ybwhHvEt9T7GOjj7k0r373Ip8AdWXKP57D8w2epu+xJeuatLlF+iLqsaXEFbcf9mNYdu1dDD2H/sf08uvRRXl3/arHzqXVT+engn3JR+4u0monUegrcIhJNKhK43wNGAb3cfUMZdToDa4BP3X1MaXVqGwXuyLVh+efsn/c4fTPfp47lFiubnPtvxHYby3eHt2dE12bExFR9wF2+dzkPLHqAlftXFjs/LHUYvxj6C7o06lLlfRKpKgrcIhJNKhK4vwO8AKwCrnb3daeUdwFeBXoBk939uRD3NSIpcEe+Q5l7WPPOVNpueJ7WvptNhS24MPchnMCc7g5Nkrl5aCsm9m1Mg8bNqrRvhV7Ia+tf49Glj5KZk3nyfKzFMqnHJKb0n0K9hOqZAiMSTgrcIhJNKhK4DXgTuATIB+YDacHi7sD5QBzwnrtfEvquRiYF7pqjsKCAFZ/M4uM1O3l4W9diZRNi5vNg/BOsaHwRTUbdQ+d+51Rp3w7nHuavy/7KS2kvFVvNpHFSY+4dcC9XdbmK2JjI2VlTpLIUuEUkmlRo4xszSwAeBr5HIFwXlQ/8E/iJux8PWQ8jnAJ3zbR53xGeW7iFl5ds43BOPrMTfsXAmPUny9Pie5Hd7xb6XTSZhMSq+/7v2gNrefCLB1mcUXyTn66NuvKzwT9jeKvhVdYXkXBS4BaRaHJGO02aWSowGmgfPLUFmOvuGWV/qnZS4K7Zjubm8+7idPp/dBOdCzaWKN9PA9a2mUincT+gRZvOVdInd2fO5jn8ecmf2X20+Lb2I9qM4CeDf6LdKqXGU+AWkWiird0rSYG7dvDCQtIXf0j2/L/R7/DHJFhBsfJ8j2FFyrnED7+T3udcViVreh/LP8b0VdN5cuWTHMs/dvJ8rMVyXffrmHLWFBomNQx7P0TCQYFbRKKJAnclKXDXPnsztrL+3cfpvOVlmpNZrCzbk7ix/tNcc05PrhrYhpTEU2dWhd6eo3v4y1d/4fX1r+N88/NaL6Eed/W7ixt63EB8bHzY+yESSgrcIhJNKhy4zaw7cB8wEmgdPL0DmAc85u5pZXy0VlLgrr3y8nJZMfdFEpZOo09uYIfIZ/PH8Mv82wBISYzj6oGtuXloK7q0bBz2/qzZv4Y/LflTifndbeu15f5B93Nhuwu1frfUGArcIhJNKvqlyVuAqUACUNrf7LnAXe4+PSS9qwEUuKPD5jVfkjH3Mf6w5xyW57YqVvZQ/FR61DlM7sDb6Dv6BuLiE8LWD3dn3rZ5PLTkIbZmbS1W1r9Zf3486McMbDEwbNcXCRUFbhGJJhVZFnAI8DkQQ2C97SeBDQSCd0fgNuBqoAA4190Xl9FUraLAHV2ycvKYvXQHzyzYzIa9R2jMYRYk/oBEywdgN03Y2P5auo7/AU1T24atH3kFebyU/hJTv55abJt4gJFtR/LDAT/UxjkS0RS4RSSaVCRwvwxMBG5y9xfLqDMJeB54xd2vD1kvI5gCd3RydxZs2M+KD57l9ozfEmeFxcpzPZblDUaRct736T54dNi+ZHkw5yB/X/53Xkp/ifzC/JPnYyyGCZ0nMKX/FFLrpobl2iKVocAtItGkIoF7J7Dd3Yd+S71FQDt3bxmC/kU8BW7J2LaBTe89Rrcds2jCoRLlG2I7sa/HTfQZdzt164VnVZHtWdt5fNnjvL3x7WJfrEyMTeSGnjdwe5/baZDYICzXFjkTCtwiEk0qEriPE3hyfdO31HsOuNbdE0PQv4inwC0nHM85yooPniVl+dP0yFtdonyuD+GTgY9w47D2dE8Nz3btaZlpPLL0Ef6141/FztdLqMf3+n6PST0mkRRXdRv5iJRFgVtEoklFAncGsMndy93qzsw+Bzq5e1T8O7YCt5Rm/df/IvPjx+mb+QF1LBeAO3N/zPuFQwAY0qERNw5rz/g+LUiMD/3Sgot2LeLhLx9m1f5Vxc63SG7B3WfdzRVdriA+RksJSvVR4BaRaFKRwD0bmABc7e6vl1HncuB14FV3nxiyXkYwBW4pz6HMvax5728kbHifa4/8jAJiT5bFkc87Sb9kf8vzaTd2Cq079Q7ptd2d97e8z/8t/b8SK5q0rdeWu8+6m0s6XkJsTGwZLYiEjwK3iESTigTu84BPCKxC8hwwHdgEONAJuBm4CYgFRrj7v8poqlZR4JbT4e4s2Lif5xduZc6qDPILnXExX/C3hEdO1lmeNJiCgbfSd9R1IV1aMK8wj9lrZzP166nsz9lfrKxTg07c0/8exrQfQ4yFf/dMkRMUuEUkmlR0He57gf8lsDRgiWICYfzH7v5YaLoX+RS4paL2ZOXw8uJtdJv/Ey4q/LRkOY3Z2O4aOl88hWatO4bsukfzjvJC2gs8tfIpDuceLlbWo3EPftD/B1zQ5gJtniNVQoFbRKLJmew0OQD4EXAB0IpA0N5B4On3o+7+Vag7GckUuOVMFeTns+KTmdiSafQ9upgYK/6zmO8xrEg5h7hhd9D73CuIiQ3N1I+s3CyeXf0sz6x+hiN5R4qV9Wvaj3sG3MPwlsMVvCWsFLhFJJpUOHBLcQrcEgo7N6ez5f3H6bbztVKXFnwo4W7qn3cn1wxqQ6O6oZlucjDnIE+teooX017kWP6xYmWDWgzi3gH3MqjFoJBcS+RUCtwiEk0UuCtJgVtCKfd4Dis+fI46Xz9Nr9wVAOR4PMOOP84hUkiIi+HSvi35zuA2DO3YOCQb6uw7to9pK6bxcvrL5BbmFisbljqMu866iyGpQyp9HZGiFLhFJJoocFeSAreEy5a0peyaO5V1e7P5ZU7x5e/72kYeq/M3dnW6lu4X30WjZpXfZyrjSAb/XP5PZq+bTb7nFysb1GIQd591N0NTh2qqiYSEAreIRJMyA7eZ/aMS7bq731WJz9cYCtwSbsdyC3jz6508v2gLX28PTDf5Q9wT3BD3EQC5HseK+ueTOPRWep1zWaXnem/L2sbfv/47b218iwIvKFbWv1l/vn/W9zmn1TkK3lIpCtwiEk3KC9yFlWjX3T0qFvdV4JaqtHLHIWYs2sg9X19NqmWWKN9uqWzrcA1dL76LpqntKnWtbYe38cTKJ3hj/Rslnnj3bdqX75/1fc5vfb6Ct5wRBW4RiSblBe7bK9Owu0+rzOdrCgVuqQ5Hsg6y+oPp1F/9PN3z00uU53lsYIWTwbfQ+/wriY07890sd2Tv4MkVTzJ7/WzyC4sH756Ne/L9s77PqLajFLylQhS4RSSaaA53JSlwS3XbtGoRez7+Bz33vkt9jpQon5TwGMOHnc11g9uS2iDpjK+TcSSDJ1c+yay1s0p8ubJ7o+7c0e8OxrYbq50r5bQocItINFHgriQFbokUOUezWfHBM6SsfI6eeasAWFTYg+tzfwVAjMGFPZpzw4BmXNCz1RnvZrnn6B6eWvkUr6x9heMFx4uVta/fnlt738rlnS8nITZ0u2VK7aPALSLR5IwCt5mlAIOBZsBWd18U6o7VFArcEom2pC1l57x/MHNPa2YdG1isbErsa9wS/yEb2lxJuwvvoHWn3md0jX3H9jF91XRmpM8osY538zrNubn3zVzT7Rrqxtc94/uQ2kuBW0SiSUW3dq8HPATcDMQHT09399uC5XcD/w5c4+5fhLivEUmBWyLZ8fwC5qzazYuLtrJg436MQj5J+DHtYvaerLMqoR9He0+iz5jvUqduvQpfIzMnk+dWP8dLaS+RlZdVrKx+Qn0m9ZjEjT1vpFFSo0rfj9QeCtwiEk1OO3CbWTLwGTAA2AcsBS4Cni4SuDsD64A/uvu/h6XHEUaBW2qKTfuOMOezBVz79W2l7maZ7XVY3WQMDc65lW4DR1V4U53s3GxeWfsKz6x+hn3H9hUrS4pNYmK3iUzuNZmWKZVfM1xqPgVuEYkmFQncvwR+A7wI3OnuR4JLB54M3MF66cBBdx8Wjg5HGgVuqWlyj+ew6uMZ2LIX6Ht0EbFW8nfAlpi27Oo0kc6X3k+zRvUr1P7xguO8seENnlr5FNuythUri7M4Lul0Cbf3uZ1ODTtV6j6kZlPgFpFoUpHAvQJoDHR295zgudIC9/tAL3dvE4b+RhwFbqnJ9u7czPoPp9Fm00za+s5iZVsLmzE6/xFG9kjlusFtGdm9GfGxp//Uu6CwgA+2fMC0ldNIy0wrUT6yzUhu7n0zg1sM1pKCUUiBW0SiSUUC91FgjrtfVeRcaYH7BeBqdz/z9cdqEAVuqQ28sJD0xR9yeMFT9Dkwl2Q7zkN51/CXgqtP1mmaksg93Q5zYb8OtO8xsJzWTmnbnX/t/BfTVkxjye6SPyu9mvRicq/JjO0wlviY+FJakNpIgVtEoklFAvch4F/ufkmRc6UF7k8JPOFuGurORiIFbqltjmQdZNWHzzBtZ3vmbCu+Yc6MhN8yLCaNtLieHO55PT3HTKZeg8an3fayPcuYtnIaH2/7uERZy7otubHnjUzsOpGUhJTK3oZEOAVuEYkmFQnci4E2QAd3Px48Vyxwm1lDYAvwlbuPDEuPI4wCt9RmG/dm88qX25n15XbqZm9mXuL9xcpzPJ6V9S8gYdCN9D5vwmnvaLnp0CaeXf0sb2x4o8Ra3inxKUzsOpGbet1Eat3UkN2LRBYFbhGJJhUJ3L8A/gA87O73B8+dGrgfA+4G7nP3x8PT5ciiwC3RIL+gkCVLFpD46QP0zv6cBCsoUWcvjdiQegktLriVjr2GnFa7mTmZzEifwUtpL5GZk1msLNZiuajDRUzuPZneTc5srXCJXArcIhJNKhK46wJLgG7AfGAW8AgwD3gJuBYYDawChpx4Cl7bKXBLtMncs4O1H0yj+YZZdCrcXGqdZXH9+HLEdK7o35pm9RK/tc2c/Bze2vgWz6x+hk2HNpUoH5I6hO/2/C4XtLlAW8fXEgrcIhJNKrrxTVsCQXsw4IAFjwT/exkwwd23ld5C7aPALdFsw4qF7J3/NF12v0tTDp48/1z+aP4z/3ZiY4yR3Zpx9cA2jO7ZnKT48sNyoRcyf8d8nl71NIszFpcob53Smkk9JnFV16uon1Cx5Qolsihwi0g0OdOt3S8DLgE6AbHANuBdYJa7F4a0hxFOgVsE8vNyWfXZa+R/9QJ9Ds9nUu7/Y6l3K1bnT0lP0q5JXRoM+y7dB4/+1o11Vu9fzfRV05mzeQ4FXnwKS524OlzR+Qpu6HGD1vOuoRS4RSSanFHglm8ocIsUd+jAPt5Ze4TZX+1g8eYDANQnm8WJ95BoeQBst5ZsazeB9iNvpVXHHuW2tyt7Fy+lv8SsdbM4dLzkDpnntDqHG3veyHmtzyPGKrY7plQfBW4RiSZlBm4zmwlMA95zpfIyKXCLlG3L/iPMXrqDI4uf5T9z/6/UOmnxvTjU5Uq6jvoujZu3KrOtY/nHeGfjOzy35jnWH1xforxdvXZM6jGJK7tcqWUFawAFbhGJJuUF7kIC87N3AdMJrEayrgr7ViMocIt8Oy8sZM0X75O96Fl6Zs6lnh0rUSfPY1mdPIjcPt+h99jJJCeUvsSgu7M4YzHPr3mej7d/TOEps9iS45KZ0GUCk3pMomODjmG5H6k8BW4RiSblBe7HgesJbOd+otJ84EngFXc/WiU9jHAK3CIVk3M0m5XzXiR+xQx6H/uSOCsemN8vGMQP+TkX9W7Blf1bc17XpmVuKb89azsz0mcwa90ssnKzSpQPbzmc67tfz4i2I4iLOb01wqVqKHCLSDQpdw63mSUAE4DbgDEEviDpwBFgBvCUu39eBf2MWArcImcuc88O1s17lgbrX6dH3moApuTexzuFZ5+s07huAv/WZjX9eveiRxlftjyad5S3Nr7FC2teYMOhDSXKmyc355pu1zCx60SaJzcP3w3JaVPgFpFoUpF1uFsBk4Gbge7B0w6sJfDU+1l3zwhHJyOZArdIaOzclMaWT57hd/tHsnpv3snzceSzKPEemlgWO605W1pdQqvzbqZ9z0El2nB3FmUs4vnVz/PJ9k9wiv9+i7M4RrUbxfXdr2do6lDMLOz3JaVT4BaRaHKmywIOJ/DU+1qgPoHgXUBgacAngbfcveRWdLWQArdIaLk7q3cd5o1lO3nj6510z1rI0wn/U6LehtiO7G1/Oe0uuIlWHbqXKN+ZvZOZa2cya92sErtYAnRs0JHrul3HFV2u0Jre1UCBW0SiSaWWBTSzOsA1wC3AyCJFe909tVI9qyEUuEXCp7DQWb50Abmf/5UemR9RnyOl1kuP686BjpfRccSNtGjTuVhZXkEeH279kBnpM/hy95clPpsUm8T4juO5vsf12kK+Cilwi0g0Cdk63GY2FngOaAa4u0fF/ssK3CJV43jOUVZ/Ohtf/jK9sj4nyfJK1Flc2I0/tnyUy/q15JJ+LWleL6lY+boD63g5/WXe3PgmR/JKhvc+TfpwbfdrGddhHMnxyWG7F1HgFpHoUtkn3CkEVjK5BTiHwPbuAFvdvUNlO1cTKHCLVL2sQ5mkzXuBhLTX6HVsKfEWmMH267ybebpgHABmMKxjY25uf4CzB5xF4+atT37+xJcsZ6TPYO2BtSXaT45L5pJOlzCx60R6N+mtud5hoMAtItHkTOdwjwJuBa4G6hAI2seBNwjM4X4/WjbLUeAWqV6H9u8m/eMXqbP2de7Mup1dhY2KlDpzE35Ke9vNmqT+HOs2ge4jJ9GgSYtAqTtf7/2aGekzmLN5DnmFJZ+ad2/UnYndJnJpp0s11zuEFLhFJJpUZJWSjgRWKZkMtOObp9nLgKeA59z9QDg6GckUuEUix/7s47y7MoO3l+9i4ab99GAL7yb+e7E6eR7L6joDye1xJd1GfIcGjZoCkJmTyZsb3mTm2plsPry5RNuJsYlc1P4iJnabyMDmA/XUu5IUuEUkmnzbOtzJBFYiuQU4n0DINuAA8AIwzd2Xhb+bkUuBWyQy7cnKYclnc+j01YMn1/g+Va7HsqbOQI53vYyuF1xPo2YtcXeW7lnK7HWzmbN5DscLjpf4XMcGHZnYdSKXd76cxkmNw30rtZICt4hEk/J2mpxGIGzXJRCyC4G5BKaMvOruuVXVyUimwC0S+XZv38CmT56n0aa36J6fXmqdtYWt+XXbJxnfJ5WLe6fSvH4Sh3MP8/bGt5m1dhbpB0p+Li4mjgvbXsjEbhM5u+XZxFjpO2JKSQrcIhJNygvcJ/Zb3gQ8TWBXye1V1K8aQ4FbpGbZuTmdrZ89T5Mt79A1f93J84/lT+DP+dcDgS9cDmrXiBvbH2R4ny60aNuV1ftXM3PdTN7Z+A5H84+WaLdl3ZZc0fkKJnSeQNv6bavsfmoqBW4RiSblBe5ngSfdfV7VdqlmUeAWqbl2bUlny/yXaLj5Pe7PvpFV3rFY+csJv2FoTDrr4rqyv+042px7PY3bd+a9ze8xa90slu9dXmq7g1oMYkLnCVzU4SLqxtetilupcRS4RSSahGwd7milwC1SO+w+nMOcVRm8uyKDRZv208QPsijxHmKs+O/ITTEdyGhzMS2HX8fx5vV4dcOrvLXxLQ4dP1SizTpxdRjbfixXdrmSQS0GacpJEQrcIhJNFLgrSYFbpPbZn32cRV8soM0X/03PIut8n2qbtWJHi1HUPetSdqbG8ebGN5m/Yz4FXrJ+65TWTOg8gcs7X06bem3CfQsRT4FbRKKJAnclKXCL1G6HMvey9tOXiVv7Fr2OLCaxlB0u93s9xsU+waheLTm7awL7bQFvb3yDDYc2lNrm0NShTOgygTHtxkTtjpYK3CISTRS4K0mBWyR6ZB8+QPpns7A1b9AjayHJFlgy8OX8Efw8/66T9erEx3Jbm+20a7KFlU0KmJvxMVm5WSXaS45LZmz7sVze+XIGtxhMbExsld1LdVPgFpFoosBdSQrcItEp59gR0he8yfGVbzH94ADePtqjWPlf4x/hktgvKHBjRWJvPmvbk68aHOfLg8sp9MIS7TVPbs6lnS7lsk6X0a1Rt6q6jWqjwC0i0USBu5IUuEWksND5attBPli9mw9WZ7B97wGWJt5FXSu5ac6S+Ha8mdqZhXWz2Xl8d6ntdWvUjcs7Xc74juNpUbdFuLtfLRS4RSSaKHBXkgK3iJxq444M9sx9nIbbPqRb7poSK50AOPBZYhNmNWnHouRcjhSUnHJiGENbDuXyTpczpv2YWrXEoAK3iEQTBe5KUuAWkfLsy9jGxn/NJGH9e/Q8+mWJL10e9zgGHJ9KYf3tNGu5iqzYZRSUspFvUmwSo9qO4rLOlzG81XDiY+Kr6hbCQoFbRKKJAnclKXCLyOk6knWQtZ+/Tv7qt+ly6HMakcUnBf2YnPeLbyrF5NC7/rs0ariQ5UmGW8l2Gic1ZlyHcVza6VL6Nu2LWSmVIpwCt4hEEwXuSlLgFpEzUZCfz7ql8/hy60Ge3d6CtIxvppT8Ou5pbol7n92xsbyTkswbdeuxPjGu1HZap7RmfMfxjOswjm6NutWY8K3ALSLRRIG7khS4RSQUdh48xkdpe5iXtodfbLqFrra9WHl6fDxvp9Tl7ZRk9sSVHr47NejEuI7jGN9hPB0adKiCXp85BW4RiSYK3JWkwC0ioZZzNJv0he+Qs+pt2u+fTyr7TpYVAEuSEnkrpS5zk5PJii19u/iejXuefPLdMqVlFfX89Clwi0g0UeCuJAVuEQknLyxk85rFZCx5nYbb59Etdw2xwVVPcoHz4u7jcP1NxNVbjcWU3AWFJ2e9AAAQq0lEQVQTYEDzAYzrMI6LOlxE0zpNq7D3ZVPgFpFoosBdSQrcIlKVDuzdxYYFr8H6D+HwTq7J+c9AgeUSl7KGRg0WUpCykbxS5nLHWAxDU4cyvuN4RrcbTYPEBlXc+28ocItINIn6wG1mNwB3A/2AWCANeAqY6l7KdnCnUOAWkepSUOis3HGIT9fu5ZO1e1m69QB3xrzBPQkzmFs3mffqJrOwThIFpYTvuJg4hrccztj2Y7mw3YVVHr4VuEUkmkR14Dazx4EpQA4wF8gDRgP1gFeBa929oLw2FLhFJFIcOprH6oXvErtiBh0OfE5zMsmMieGDusm8WzeZpUmJeCnhO9ZiGdZyGBd3uJhRbUfRKKlR2PuqwC0i0SRqA7eZTQRmAhnABe6+Lni+BTAP6An8yN0fLa8dBW4RiUReWMjmtKVkLH2LlG0f0z1nBZlxzvt1k3k3JZmViYmlfi6GGPo36stlPSYwut1oGic1Dkv/FLhFJJpEc+BeAgwCJrv7M6eUjQA+JhDGW5c3tUSBW0RqgqPZh1j3xbvkrH6f1vs/h9g93Jw4gYx6u4mts63Uz5hDz8T2jO9xNZf1uCKkX7hU4BaRaBKVgdvM2gDbCHzJv6G7HyulznagNXCuu39eVlsK3CJSE+3YmMZne+L5bMMh5m9ey9H4ZSTVWwbJO0qtbw5dvSHnpo7k+uF30rph20pdX4FbRKJJtAbuy4E3gK/cfWAZdV4FrgR+4O6Pl9WWAreI1HSFhc7qXYdZsXwJLZf/ms1xW/goJZGvkpJKrW/udMtNoFWL7zB50LUMatOxwtdU4BaRaFL6dmW134m/HbaUU2frKXVFRGqlmBijT+sG9Gk9GsaPJudoNoO/nMvWNW+z9shiltXJLvaFSzcjPTGP9IPPcnRpPk+0+Y9qvgMRkcgWrYE7JXg8Uk6d7OCx3qkFZnYncCdAu3btQtszEZFqlpScQp/zJ9Dn/AlcAhzcl8HSha/wxfZ3WMFWViYZhcHwfX2vS6q3syIiNUC0Bu4T62Kd0Xwad/8H8A8ITCkJVadERCJRw6apXHjZvVzIvQCsTFvIS2vnsiRrJ2M6n1XNvRMRiXzRGrizgseUcuqcKMsqp46ISNTp0+Ns/rvH2dXdDRGRGiOmujtQTTYHj+3LqXPiK/iby6kjIiIiIlKuaA3cXwWPvc2sThl1hpxSV0RERESkwqIycLv7NmApkABce2p5cOObNgQ2vllQtb0TERERkdokKgN30APB4x/NrMuJk2bWHPhr8O2D5e0yKSIiIiLybaL1S5O4+0wzmwrcDawwsw+BPGA0UB94DXisGrsoIiIiIrVA1AZuAHefYmbzgXuAEUAskAY8CUzV020RERERqayoDtwA7v4C8EJ190NEREREaidz174tlWFmeyl/i/jyNAX2hbA70UBjVjEar4rReFVMZcarvbs3C2VnREQilQJ3NTKzJe4+uLr7UZNozCpG41UxGq+K0XiJiJyeaF6lREREREQk7BS4RURERETCSIG7ev2jujtQA2nMKkbjVTEar4rReImInAbN4RYRERERCSM94RYRERERCSMFbhERERGRMFLgDhEzu8HMPjOzQ2aWbWZLzOweMzvtMTazeDMbbWYPmdlCM9tlZrlmtsPMZprZyDDeQpULxZiV0/YfzMyDr5+Gor/VLdTjZWZ1zOznZrbYzA6a2VEz22Rmr5jZuaHuf1UL5XiZWRsz+4uZpZvZMTPLMbN1ZvY3M+sUjv5XFTPrbmY/NLPnzCzNzAqDPzfXVLLdsP18i4jUNJrDHQJm9jgwBcgB5gJ5wGigHvAqcK27F5xGO2OAD4JvM4AvgSNAL6BP8Pzv3P1XIb2BahCqMSuj7SHAAgL/Q2nAz9z9z6Hod3UJ9XiZWUfgfaALsAdYCBwHOgD9gd+6+3+H8BaqVCjHy8wGAB8BDYHtBH4uAQYDrYFs4GJ3/zyU91BVzOwR4IelFF3r7jPPsM2w/XyLiNRI7q5XJV7ARMCBXUDXIudbAKuDZT88zbYuBGYC55dSdj2QH2xvVHXfd6SMWSltJwKrgB0E/mJ34KfVfc+RNF5AXWB98HO/BeJPKW8CdKvu+46g8fo8+Jl/FB0rIB6YFiz7urrvuxLjdQfwP8B1QGfg4+A9XRMJ46+XXnrpVRte1d6Bmv4ClgT/Arm5lLIRRf7iiQnBtZ4Itjetuu87UscM+GPw85cDT9eSwB3S8QIeCH5menXfW6SPF5AUrO9AainlrYqUJ1f3vYdo/CobuKvsd6JeeumlV015aS5dJZhZG2AQkAu8cmq5u39C4ElrKnB2CC75VfDYJgRtVYtwjpmZDQPuB15w9zcr39vqF+rxMrME4HvBtw+GrqeRIQx/vgoI/MsSBKYnlWgyeDwCHKtof2ubavidKCJSIyhwV86A4HGVu5f1l+3iU+pWRtfgcVcI2qouYRkzM0sCpgOZlD4ftaYK9XgNIjBlZJu7rzGzc4JfMP27mf3GzIZXtsPVLKTj5e55BOYgA/zGzOJPlAX/+8Q892nuri/EVP3vRBGRGiGuujtQw3UMHreUU2frKXXPiJmlArcE386qTFvVLFxj9nugO/Add993Jh2LUKEer77B4zozexqYfEr5r8xsFvDdcgJTJAvHn68pwHsE/mVgvJktCZ4fAjQCHgV+VsF+1lZV9jtRRKQmUeCunJTg8Ug5dbKDx3pnehEziwOeAxoAc2v4dImQj5mZnQP8CHjN3WdUom+RKNTj1Th4vACIBf4M/A3YHzz3VwJfejsM3FbRzkaAkP/5cveNwT9jzwDjKT6lawnwafBJuFTR70QRkZpGU0oq58ScznD/U/LfCCyptQ24KczXCreQjpmZ1QGeIhAQp4SizQgT6j9jJ37m4whMg/iZu29w94Pu/gZwZfBak2vo+tIh/5kMhu2VBJZQnAA0BZoRGKtGwCwzq/FLdYZIVf1OFBGpURS4KycreEwpp86Jsqxy6pTJzB4FbiewLvdod884k3YiSKjH7A9AN+An7l6T57aXJdTjVbTOP08tdPclBNaZjgFGnkZ7kSak42VmDYHXCDyNHefub7j7fnff5+6vA+MIfFnyl2bWtby2okTYfyeKiNRECtyVszl4bF9Onban1D1tZvYQcB+wl0DYXlfRNiLQ5uAxVGN2FVBI4Insx0VfBMIQwN3Bc0+cQX+r2+bgMVTjVbTOpjLqnDifehrtRZrNwWOoxutSAk+zF7r7xlML3X09sIjAvxiMPN1O1mKbg8ew/E4UEampNIe7ck4s09fbzOqU8SWzIafUPS1m9j/ATwjMrR3r7qvPvJsRJRxjFkNgfd+ydAq+Gp5me5Ek1OO1tMh/NyHwP3Onaho8ZpdSFulCPV7tgsdD5dQ5GDw2LqdOtAjb70QRkZpMT7grwd23EQgwCcC1p5ab2QgCX7DKILDV+GkxswcJrHpwgEDY/jokHY4AoR4zd+/g7lbai8AygRDY2t3cvX/o7qRqhGG8dhB4IguB7wWc2l4jYGDw7ZJTyyNdGH4mdwaPg4ouCVikvXgCSy1C2f9iEDXC9TtRRKSmU+CuvAeCxz+aWZcTJ82sOYEVHwAedPfCImUPmFmamT3AKczsd8C/EXhqNtbda+NToJCOWRQI9Xj9Pnj8lZn1L/KZJGAqgdVwvqTmBqJQjte7wFECT7ofNrPEIp9JBP6PwBSJA8CckN9JhPqWP18VHn8RkdpOU0oqyd1nmtlU4G5ghZl9COQReHpYn8AXrh475WMtCawZ3bLoSTO7AvjP4Nv1wL1mpW1uR5q719hdAkM5ZtEg1OPl7m+a2Z+BnwKLzGwRgalLQwlsVb4DmFRTN3IJ5Xi5+x4zmwJMA+4BrjKzLwmsxjEoWP84cJu7lzftJGKZ2UC+CcIAvYLHP5jZT0+cdPeiO0OW9+frTMZfRKRWU+AOAXefYmbzCfyFPILA+sZpwJPA1Ao8ySk6B3Rw8FWaT6jh23KHcMyiQqjHy91/ZmafA/cS2PEvmcCGJP9L4OljaXO7a4xQjpe7TzezFQTWej8fuChYtINAEP/fGv4di/rAsFLOn/GqK/r5FhEpzmroQywRERERkRpBc7hFRERERMJIgVtEREREJIwUuEVEREREwkiBW0REREQkjBS4RURERETCSIFbRERERCSMFLhFRERERMJIgVukBjGztmb2vJntNLN8M3Mze6SUej8Olt1cBX2KMbN7zGyJmWWb2SEz+8zMJoX72iIiIjWBdpoUqSHMzIBZwBBgNTCPwJbZX5RS/WogH3grzH2KBWYDVwCHgfeBRALbeL9gZsPd/b5w9kFERCTSaadJkRrCzDoCGwlswd7Z3fPLqNcC2Al85O5jw9yn+4E/E/gfgAvdfXfwfFfgM6AFcKW7vx7OfoiIiEQyTSkRqTnaBo+bygrbQVcS+Nl+NZydCT7d/nnw7d0nwjaAu68D/i349v+Fsx8iIiKRToFbJAyC86c9+N+3BOc3HzGzDDObZmbNgmVJZvYbM1trZjlmttXMfm9m8UXa6hBs65PgqREn2j9xjVNcDTjwWpE2fh2s/2sza2NmT5vZLjM7amZLzeyaInXPNbN3zGx/sHyemQ0p5TrDgebAdnf/tJTyVwhMeRliZq0rNIAiIiK1iAK3SBiZ2R+BvwOZwHsEgvBtwIdmlgLMBe4FVgEfAU2A/wAeL9JMNjAdmBN8vzv4/sSr6PUaAKOAhe6+s5QutQe+BM4nEOCXAgOAl83sO2Z2FYG54U2BD4AtwEhgnpl1O6WtAcHj4tLu3d2PBu8LoH9pdURERKKBvjQpEl6Tgf7uvgbAzBoBC4B+weNBoKO7HwqW9ycQYO8ws9+7+xZ33wfcYmYjgYuBNHe/pYzrXQ7EU/Z0kluAR4H73b0geM27gb8CfwLqAje6+yvBshjgBeB6AlNEbi/SVsfgcUs597+VQNjuWE4dERGRWk1PuEXC61cnwjaAux8A/hZ82wu480TYDpYvA94BDBhxBte7OnicXUb5FuDnJ8J20D+A/UAb4L0TYTvYn0Lgj8G3o05pKyV4PFJOf7KDx3rf0m8REZFaS4FbJLzeK+Xc+uBxS9EwXsS64LFVRS5kZskEnoCvcPcNZVT7yN1zi54Ihu/N5fS3rP7YiSYq0k8REZFoo8AtEl7bSzmXXU5Z0fKkCl5rHJBM2U+3T+eaJcrd/URZ4ilFWcFjCmU7UZZVTh0REZFaTYFbJIyCUzLKUl7Zmfi26SSnc82K9Glz8Ni+nDonljLcXE4dERGRWk2BW6QWCC4jeCmw0d2XV9FllwaPpS0ZeGKKS5/g26+qpEciIiIRSIFbpHa4EGhI+U+3Q20BsAdoY2YXlFJ+LYEVUxa7+44q7JeIiEhEUeAWqR1OTCcJ6+6SRQW/bPmn4NupZtb8RFlwa/cHg29/X1V9EhERiURah1ukhguulT0B2EXgqXNVehi4gMD63+vMbC6Bp9pjCHzp8y/u/noV90lERCSi6Am3SM13LtACeM3dq3SJvuBT7isJ7Ja5nsCyhCMI7GZ5o7vfV5X9ERERiURWxX8/i0iImdnDwI+Ai9z9g+ruj4iIiBSnJ9wiNd8a4NfAvGruh4iIiJRCT7hFRERERMJIT7hFRERERMJIgVtEREREJIwUuEVEREREwkiBW0REREQkjBS4RURERETCSIFbRERERCSMFLhFRERERMLo/wNDU54mrtAbHQAAAABJRU5ErkJggg==\n", | |
"text/plain": [ | |
"<Figure size 432x288 with 1 Axes>" | |
] | |
}, | |
"metadata": { | |
"needs_background": "light" | |
}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"plt.plot(num_heun_drag[:,2]/0.25, num_heun_drag[:,1],'-', label = 'Heun')\n", | |
"plt.plot(num_rk2_drag[:,2]/0.25, num_rk2_drag[:,1], '--', label = 'Runge Kutta')\n", | |
"plt.plot(d_m,vf,'-', label = 'Analytical Solution')\n", | |
"plt.ylabel('Velocity (m/s)')\n", | |
"plt.xlabel('mf/m0')\n", | |
"plt.legend(loc=9, bbox_to_anchor=(1.5,1));" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 80, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"Final height reached for simplerocket integration is 597.4795789904731 meters using Heun method and 597.479575608677 meters for Runge Kutta method\n", | |
"Final height reached for rocket integration is 425.35224719755064 meters using Heun method and 425.35224735150734 meters for \n", | |
"Runge Kutta method\n" | |
] | |
} | |
], | |
"source": [ | |
"height_heun = num_heun[:,0]\n", | |
"height_rk2 = num_rk2[:,0]\n", | |
"height_heun_drag = num_heun_drag[:,0]\n", | |
"height_rk2_drag = num_rk2_drag[:,0]\n", | |
"print('Final height reached for simplerocket integration is',height_heun[-1],'meters using Heun method and',height_rk2[-1],'meters for Runge Kutta method')\n", | |
"print('Final height reached for rocket integration is',height_heun_drag[-1],'meters using Heun method and',height_rk2_drag[-1],'meters for \\nRunge Kutta method')\n" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"2. As seen by the values above the simplerocket integration reaches a height greater than the rocket integration since drag is not factored in." | |
] | |
}, | |
{ | |
"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": 81, | |
"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", | |
" y0 = 0 # initial position\n", | |
" v0 = 0 # initial velocity \n", | |
" time = (0.25-0.05)/dmdt\n", | |
" t= np.linspace(0,time,10000)\n", | |
" dt=t[1]-t[0]\n", | |
" \n", | |
" N=int(time/dt)\n", | |
" height_desired = 300\n", | |
"\n", | |
" #initialize solution array\n", | |
" height = np.zeros([N,3])\n", | |
"\n", | |
" #Set intial conditions\n", | |
" height[0,0] = y0\n", | |
" height[0,1] = v0\n", | |
" height[0,2] = m0\n", | |
"\n", | |
" for i in range(N-1):\n", | |
" height[i+1] = rk2_step(height[i],lambda state: rocket(state,dmdt=dmdt, u=250),dt)\n", | |
" height_predicted = height[:,0]\n", | |
" error = height_desired - height_predicted[-1]\n", | |
"\n", | |
" return error" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 82, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"def incsearch(func,xmin,xmax,ns=50):\n", | |
" '''incsearch: incremental search root locator\n", | |
" xb = incsearch(func,xmin,xmax,ns):\n", | |
" finds brackets of x that contain sign changes\n", | |
" of a function on an interval\n", | |
" arguments:\n", | |
" ---------\n", | |
" func = name of function\n", | |
" xmin, xmax = endpoints of interval\n", | |
" ns = number of subintervals (default = 50)\n", | |
" returns:\n", | |
" ---------\n", | |
" xb(k,1) is the lower bound of the kth sign change\n", | |
" xb(k,2) is the upper bound of the kth sign change\n", | |
" If no brackets found, xb = [].'''\n", | |
" x = np.linspace(xmin,xmax,ns)\n", | |
" f = np.empty(ns)\n", | |
" for i in range(ns):\n", | |
" f[i] = func(x[i])\n", | |
" sign_f = np.sign(f)\n", | |
" delta_sign_f = sign_f[1:]-sign_f[0:-1]\n", | |
" i_zeros = np.nonzero(delta_sign_f!=0)\n", | |
" nb = len(i_zeros[0])\n", | |
" xb = np.block([[ x[i_zeros[0]+1]],[x[i_zeros[0]] ]] )\n", | |
"\n", | |
" \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": 83, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"def mod_secant(func,dx,x0,es=0.0001,maxit=50):\n", | |
" '''mod_secant: Modified secant root location zeroes\n", | |
" root,[fx,ea,iter]=mod_secant(func,dfunc,xr,es,maxit,p1,p2,...):\n", | |
" uses modified secant method to find the root of func\n", | |
" arguments:\n", | |
" ----------\n", | |
" func = name of function\n", | |
" dx = perturbation fraction\n", | |
" xr = initial guess\n", | |
" es = desired relative error (default = 0.0001 )\n", | |
" maxit = maximum allowable iterations (default = 50)\n", | |
" p1,p2,... = additional parameters used by function\n", | |
" returns:\n", | |
" --------\n", | |
" root = real root\n", | |
" fx = func evaluated at root\n", | |
" ea = approximate relative error ( )\n", | |
" iter = number of iterations'''\n", | |
"\n", | |
" iter = 0;\n", | |
" xr=x0\n", | |
" for iter in range(0,maxit):\n", | |
" xrold = xr;\n", | |
" dfunc=(func(xr+dx)-func(xr))/dx;\n", | |
" xr = xr - func(xr)/dfunc;\n", | |
" if xr != 0:\n", | |
" ea = abs((xr - xrold)/xr) * 100;\n", | |
" else:\n", | |
" ea = abs((xr - xrold)/1) * 100;\n", | |
" if ea <= es:\n", | |
" break\n", | |
" return xr,[func(xr),ea,iter]\n", | |
"\n" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 98, | |
"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_interval = incsearch(f_m,0.05,0.4)\n", | |
"print('3a. The two closest mass change rates within that interval are',rate_interval[0],'and',rate_interval[1])\n" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 65, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"3b. The root of the function fm is 0.07910372921028956\n" | |
] | |
} | |
], | |
"source": [ | |
"fm_root = mod_secant(f_m,0.0001,0.08571429,es=0.000001)\n", | |
"print('3b. The root of the function fm is', fm_root[0])" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 88, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"y0 = 0 # initial position\n", | |
"v0 = 0 # initial velocity \n", | |
"dmdt = 0.07910372921028956\n", | |
"time = (0.25-0.05)/dmdt\n", | |
"t= np.linspace(0,time,10000)\n", | |
"dt=t[1]-t[0]\n", | |
" \n", | |
"N=int(time/dt)\n", | |
"height_desired = 300\n", | |
"t_plot = np.linspace(0,time,N)\n", | |
"\n", | |
"#initialize solution array\n", | |
"height = np.zeros([N,3])\n", | |
"\n", | |
"#Set intial conditions\n", | |
"height[0,0] = y0\n", | |
"height[0,1] = v0\n", | |
"height[0,2] = m0\n", | |
"\n", | |
"for i in range(N-1):\n", | |
" height[i+1] = rk2_step(height[i],lambda state: rocket(state,dmdt=dmdt, u=250),dt)\n", | |
"height_predicted = height[:,0]\n", | |
"\n" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 96, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"Text(0.5, 1.0, 'Height vs Time\\n')" | |
] | |
}, | |
"execution_count": 96, | |
"metadata": {}, | |
"output_type": "execute_result" | |
}, | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAb0AAAFqCAYAAACOFYTTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd5xU5dn/8c+1BVi6VKWDVFGKYFfAllgSO3bFFo09+hg1xZKfeaIx0Wis8VEsscUSe2ISjGCXIiogIKBI77K0hW3X749zZnd2mZ3dnZ3d2Z35vl+veZ0997nvMxezy157zrmLuTsiIiKZICvVAYiIiDQUJT0REckYSnoiIpIxlPRERCRjKOmJiEjGUNITEZGMoaQnTY6ZnWdmHr7GJbt+XeJJ9rkb4vxNmZk9EfW9TeT1RNS5Fodlt6buXyT1TUlPJAOY2a3hL/TFqY5FJJVyUh2AiEgdXAJcUcWxfwIHA0uAoVXUKaqPoKTxUtITqSN3fwJ4IsVhZCR33wHsiHXMzErKq/mWGpyrTxJDk0ZKtzdFRCRjKOmJhMwsz8x+ZmZTzGytmRWa2Soze9XMjonTrtqOJmaWa2b/Y2YzzWyrmW0ws/fN7JzweI2fuZlZJzP7g5ktMLPtZrbezP5pZofGqDsujOuWsKh3jM4ck6t7z/Bcx0S1OaCaur3NrDSs+5MY8d9mZtPNLN/MisxstZnNNrOnzOxsM2vwu1DxOrJEdZiZHO7vY2YvmNkKMysws7lm9kszaxHVpr2Z3Wxmc8Lv+Xoz+7uZ7VGDWHYxs1+b2adhux1mtszMnqvus5dquLteejWpF3Ae4OFrXDLqA8OAxVH1Yr0mAtnxzl/FudsCn8Q57xPAreHXi6uJfw9gWRXnKQXOrdR2XDX/Jgcm1/BzzwFWh20eqKbuL8J624H2UeV7RJ0j3qtTEn5OJlf1mVZRP/L9vzXGsScinxUwgeBZYKy4/xl+Tn2AeVXU2QQMjxPHocD6aj6f/5fq/4dN9aUrPcl4ZtYLeBfoTdDp4RKgP9AB2BO4EygBzgd+k8BbPAbsF379BLA30DHcPknwS/ScGp7rDYJfuOcCPYHOwInAUsCAB8ysY1T994E2wO3h/pJwP/p1dE3e2N2LgefD3VPNLDdO9bPC7VvuvjGq/C9AF2AtcCkwmOCz6A+MBW4iSBaN1QDgEYLkNxboBAwCHg2PHwVcCLxA8MfOBQTfp64E37PNBJ/5Q7FObmZ7EyTODsBs4GyCBNoBGEXwhxfATWZ2UTL/YRkj1VlXL71q+6Lilc/RQOtqXpdE1R8X43yvh8eWAZ2reM8Lwzo7gG5VxROj3QFR731fFed+KKrO4mr+vcuALjHq7B1V56cxjt9a1flr+dmPjnqfH1dRZ0RUnROiyttGlR/fAD8nk2vzb6ZmV3oOvEnsK/4PwuNFQD6we4w6F0WdZ3CM41+Exz4H8qqI87awztqq6uhV9UtXetLU/YPgr+d4r4eramxm/YAfhbs/c/e1VVSdCCwCmgHjaxHfueF2G/CrKurcCBTU8Hz/z93XVC5098+AL8PdfWoRX624+3TKr8TOqqLa2eF2A8H3JyI76uvlSQ6tIV3j7iUxyiNXwTnAn919UYw6fyNIWAD7Rh8In8kOC3cvcveqfiZ+B2wluMr8QW0CF3VkETmc4LZgKfCRmbWO9QJaEfwVDsHVTk0dGG4nu/umWBXcPR+YUsPz/TPOsfnhdtcanitRT4fb48ysTfQBM8sCTg93X3T3wsgxd/8e+C7cfSC8ldfUfOPuC6o4Fp3k/hWrgrtvJrhCg52/T0eE2/XA/Dg/i9mU/+FRm59FQUlPmr5D3d3ivQiexVVlULjNIrj6iHfFeFJYt3Mt4usTbufHq0TNn2OtiHNsW7htWcNzJeoZgquVPODkSscOBbqHXz/Nzq4J2+4LzAh7TD5lZj8xs771FXASxfv8o6/MVtagXl6l8sjPYkeCzi7xfhZHhXVr87MoKOmJtEugTYvqq5RpFW63VlOv2sHTAFXcVqvManKuRLn7YoLnV1B+KzMicsvzW+DDGG1fIUiMkwg6B/Um6MTzCPCNmb1nZvtWbteI1OTzr2m9yt+n+v5ZFDQji0gk2axz9/r4q3krQQeOVtXUa10P712fngYOAQ41s27uviIcoxa58nvGw14Xlbn7FGCKmbUj6OhzEPBDgmeRhwDvm9kh7j613v8VjUvkZ3G6u9fbc9lMpys9yXTfhNtOZtajHs4feYY1sJp6g6o53ti8QNCTNQs4Iyz7MUGCh9i3Nitw93x3f9vdb3L3fQnGFBYQdBb6ZdIjbvwiP4tDzKx5SiNJY0p6kun+E/V1vGd/iYrc4ju0cqePCDNrS/ALvz5FJlbOjlurhjwYe/dWuHt2pe10d6/uGWasc06h/PsxpG4RNkn/DretgFNTGUg6U9KTjObu8yjvVv8LMzswXn0z62Jmu9TiLf4ablsSjK+K5XZ27tSQbOvDbeckTvEV+beNMLNDCAZmR5dXEE4/1jHWsfB4FhDpzLK+qnpp7N8EA9IB7jKzuHcHzKyPrghrT0lPBH4KrCFIPO+a2V1mtn/kl7SZDTGzM83sOYLblbvX9MTu/hHw93D3ajN7zMyGm1kHMxthZhOByyi/tVVfZoTb5sD/M7NuFswHmmNmiV79/YNgLB7AUwS3JaNnbalsT2CpmT1rZmeY2eDwc+gWjlF7FdgrrPtcgjE1WeEz0HMJbvF2BqaZ2S1mNjL8nDqb2TAzu8DMXgcWEszuIrWgjiyS8dx9qZmNJUhOQ4Brw1fM6tR+DbbIVFT7hF9fUOn4UwSzgdxMkDSSzt2nmdlHBOMGfxG+IqaQwO1Vdy80sxcI/mjoExb/J9bg+Sh5BM8Az4hT52ngwdrGkw7cfaaZHUHwzLQ7wUw6t1ZRvYSa9yaVkK70RCi7zTmMYB7MNwjGYxUSdNZYCrwNXAX0dPcvqjpPFefOJ1jM9OcEs6YUABuBj4AL3H0C5b03N9f5H1O1Y4A/AF9RPqavrirfyozXgeUjggHYtxPMCbqY4LPYEX79N+Aodz+nhkMz0lJ4d2AgcCXBM87VBH9oFRAMBXmNYFq8ruGAf6kFq6JXsYg0IDN7DTgOeNPdf5zqeETSla70RFIsnFoqshbejHh1RaRulPRE6lk4Z2KzOFXupLxDwgsNEJJIxlLSE6l/ewLzzOyGqJ54Xc3siLAX3qVhvefc/asUximS9vRMT6Semdn+wMfVVJtCsMZcfgOEJJKxlPRE6lk4E8tZBPNLDiVYObwlwQDszwjGpD3r7qUpC1IkQyjpiYhIxtAzPRERyRhKeiIikjGU9EREJGMo6YmISMZQ0hMRkYyhpCciIhlDSU9ERDKGkp6IiGQMJT0REckYSnoiIpIxlPRERCRjKOmJiEjGUNITEZGMoaQnIiIZQ0lPREQyhpKeiIhkDCU9ERHJGEp6IiKSMZT0REQkYyjpiYhIxlDSExGRjKGkJyIiGUNJT0REMoaSnoiIZAwlPRERyRhKeiIikjGU9EREJGMo6YmISMZQ0hMRkYyRk+oApGqdOnXyPn36pDoMEZEmZcaMGevcvXOsY0p6jVifPn2YPn16qsMQEWlSzOy7qo5lxO1NM7vSzF4ws7lmtt7MisxsrZlNMrOzzczitD3TzN43s3wz22Jm083scjOL+9kl2k5EROpPplzp3QB0AWYDHwFbgd7AYcDhwClmdpK7l0Y3MrMHgMuA7cA7QFFY/37gcDMb7+4lld8s0XYiIlK/MiXpnQ7MdPet0YVmNpQgKR0PTAAejzp2MkHiWgWMcfcFYXlX4F3gROAK4N5K50yonYiI1L+MuNXm7h9UTnhh+RzggXD3yEqHfxFub4gkrrDNauDScPfGGLcrE20nIiL1TL94oTjcbo8UmFkPYBRQCLxYuYG7TwGWA7sC+9e1nYiINIyMTnpm1hf4abj7RtShkeF2jrsXVNF8WqW6dWknIiIh37QSHj8aNq9O+rkzKumZ2flm9oSZPWNmU4CvgR7A7e7+SlTVvuG2ym6vwJJKdevSTkREgJJS552/XId/9zE++Y6knz9TOrJEHETQYSWiGLgJuLtSvdbhdqfngFG2hNs2SWgnIiK/7UJ28Q6OiOzPmBi8cprDr9ck5S0y6krP3S9ydwNaAkOBe4BbgU/MrFtU1ci4Pa/lWyTarvwEZheHY/qmr127NtHTiIg0OQtO/5DXSw+kwJsBUJTVHPYaD1fPStp7ZFTSi3D3Anf/yt1/TtDbcjjBGLqIzeG29U6Ny0WObY4qS7RddGyPuPtodx/duXPMWXRERNLO9qISrnxzJZtK82hOEYXkkuNF0LwttOmatPfJyKRXSWRs3o/NLDf8enG47R2nXc9KdevSTkQko9317/nMW7WZTraJ5/wIVp36FjbqfNiS3M4smfZML5aNBM/2coAOwGpgZnhsqJnlVdETc59wOzOqLNF2IiIZ66NF63j0g28B+GnRNfzmuKH02qMP7LFf0t9LV3owhiDhbQTWAbj7UuAzoBkwvnIDMxtL0OtzFfBxpDzRdiIimSq/oIjrXvgCD3tCjBnYmXMPiHezrG7SPumZ2SFmdpaZNY9x7CDgsXD3sUrzYd4ebn9vZv2j2nQBHgx376g8X2cd2omIZJybX5vNivxgbpD2LXP5wynDiLMGQJ1lwu3N3Qme291vZp8RXGW1Ccv3COu8RTB0oYy7v2RmDxFMHTbLzCZRPnF0W+BVKnZ+qVM7EZFM8/oXK3jt8xVl+7efuBdd27ao1/fMhKQ3BbgNOAQYCBxIMLRgFfAy8LS7vxqrobtfZmYfAJcDY4FsYB4wEXioqqu1RNuJiGSKFRsL+PUr5UMRTt67B0fvtVu9v2/aJz13/xa4uQ7tnwWebah2IiLprrTUue7FL9i0PZj6uMcuedx63B7VtEqOtH+mJyIijcvED7/lo0XrATCDu08dQZsWudW0Sg4lPRERaTDzV23mzn/NL9v/6djd2bdvhwZ7fyU9ERFpEDuKS/jZ3z6nsDjo1jC0W1uuOWJgg8agpCciIg3i7n9/zdyVmwBonpPFPaeNoFlOw6YhJT0REal3HyxYx1/e+6Zs/8ajBzOga8MvNqOkJyIi9Wr9lh1c88LnZftjBnZmwgF9UhKLkp6IiNQbd+fnL33J2s07AOjUuhl3jR9OVlb9zboSj5KeiIjUmyc/Wsx/55UvAPvH8cPp3GanWSEbjJKeiIjUi7krN/G7f84r27/w4L6MG9QlhREp6YmISD0oKCzhyudmVhiecP1Rg1IclZKeiIjUg9ve+oqFa7YAkJebzZ/PGEnznOwUR6WkJyIiSfb27JU8++mSsv1bj9uD3Tu3TmFE5ZT0REQkaVZsLOCGl8tXTzh2r904dXTPFEZUkZKeiIgkRUmpc83fPie/oAiA7u3z+N1Je9XrorC1paQnIiJJ8cC7C/n02w0AZBnce/oI2uU1zOoJNaWkJyIidfbxovXcM+nrsv2rDx/I6D4Nt3pCTSnpiYhInazbsoOrn59JqQf7+/XtwBWH9U9tUFVQ0hMRkYSVhs/x1oTTjHVs1Yw/nzGS7BRNM1YdJT0REUnYg5MX8v6CdUCwCvqfThtB17YtUhxV1ZT0REQkIZ98s567/1P+HO+ycbszZmDnFEZUPSU9ERGptcrP8fbt06HBV0FPhJKeiIjUSuQ53upNwXO8DuFzvJzsxp9SGn+EIiLSqDw0ZVHZczyAu08dzq7tGu9zvGhKeiIiUmNTv93AXf+eX7Z/6bjdU75cUG0o6YmISI2s37KDK5/7rOw53ujeu/A/Rzb+53jRlPRERKRaJaXO1c+XP8fbpWUu953ZNJ7jRWta0YqISErcM+lrPlgY/RxvBLu1y0thRIlR0hMRkbjembua+/67sGz/ysP6c+jgpvMcL5qSnoiIVGnJ+m1c87fPy/YPGdCJnzWB8XhVUdITEZGYtheV8NOnZ7BpezEA3dq14N7TG++8mjWhpCciIjHd/Npsvlq5CYDcbOPBs0fRoVWzFEdVN0p6IiKyk+enLuGF6cvK9m/+8VBG9GyfwoiSQ0lPREQqmLUsn5tfn1O2f9LI7py9X68URpQ8SnoiIlJm47ZCLn1mBoXFpQAM3rUN/3viXpg13ed40ZT0REQECCaS/tnfPmfZ9wUAtGmew0NnjyKvWXaKI0seJT0REQHgT5O+ZvL8tWX7fzx1OH07tUphRMmnpCciIrw9e2WFAeg/Hbs7Pxy6awojqh9KeiIiGW7+qs1c+8IXZfuHDOjEz384KIUR1R8lPRGRDJa/rYiL/zqdbYUlAPTq0JL7zmjaA9DjUdITEclQJaXOlc/P5Lv12wBo2SybR84dRfuWTXsAejw5iTY0s87ACKAr0B74HlgDzHT3dfHaiohI6v3hX/N57+vyjit3jR/O4F3bpjCi+lerpGdmPYBLgOOBoXHqzQFeBR5x92VV1RMRkdR444sVPDxlUdn+FYf25+i9dkthRA2jRknPzHYHbgdOiGrzPTAX2ABsAtoCHYHBwJ7h60YzewX4hbt/k9zQRUQkEV+t2MT1L31Ztn/Y4C5c08RWQE9UtUnPzO4ErgKaAdOBJ4FJ7j4/TpvBwJHABGA8cLyZ/dndr09K1CIikpDvtxZy8V+nU1AUdFzp26kVfzptRNp2XKmsJh1Z/gd4Axjm7vu6+wPxEh6Au89z9/vcfTQwHHgTuLbu4YqISKKKSkq5/NnPymZcad08h/87dxTt8nJTHFnDqcntzdHuPjPRN3D3WcApZjYy0XOIiEjduDu3vj6HjxatLyu7+9Th9O/SJoVRNbxqr/TqkvDq4zwiIlJ7T338Hc98uqRs/9ojB/KDNJxxpToapycikube+3otv3mjfKmgHw/vxpWH9U9hRKmjpCciksYWrtnC5c9+RqkH+8N7tOMPpwxLm6WCaiuhpGdmu5jZr8xskpl9ZWbfVPFaVP3Z6peZ5ZrZ4WZ2l5l9YmYrzazQzJab2UtmNq6a9mea2ftmlm9mW8xsupldbmZxP7tE24mIJMv3Wwu58MlpbN5eDMCubVvwf+eOpkVu+iwVVFu1npHFzPoDU4Bdger+VPBEgkqyscB/wq9XATOArcAewMnAyWZ2m7vfXLmhmT0AXAZsB94BioDDgfuBw81svLuXJKudiEiyFJWUcukzM8qmGMvLzebRCaPp0rZFiiNLrUSmIbsL2A14H/gTsADYksygkqwUeBm4193fjz5gZqcBzwA3mdm77v5u1LGTCRLXKmCMuy8Iy7sC7wInAlcA91Y6Z0LtRESSxd25+bU5fPLNhrKyP502nD27t0thVI2DudfuYszM8oH1wGB3L6yXqBqQmT0KXAhMdPcLo8qnA6OACe7+VKU2Y4HJBImtu7uX1rVdLKNHj/bp06cn/o8TkYz0+Iff8ps3virbv+4HA7nisAEpjKhhmdmMcJz4ThJ5vuTA1HRIeKHIUIoekYJwjtFRQCHwYuUG7j4FWE5wi3f/urYTEUmWd+ev4bY3yxPe8SO6cfmhmdlTM5ZEkt7nBL+000Xkz5+VUWWRgfRz3L2ginbTKtWtSzsRkTqbsyKfK54p76k5omd7fn9y5vbUjCWRpPdH4GAzOzDZwTQ0M9sVOC/cfTnqUN9w+12c5pFRnn2jyhJtJyJSJyvzC7jgiWlsDReD7d4+j0fOHZXRPTVjqXVHFnd/08yuAd4ys/uBfwHLCDqMxKq/JFZ5qplZDvA00A54x93fiDrcOtxujXOKSOed6Dl8Em0nIpKwzduLOP/xaazetAOANi1yePz8fejSJrN7asaS6CKyM4HVwC/DV1W8Du9R3x4mGEawFDi70rHIvYDaDrlItF35CcwuBi4G6NWrV6KnEZEMEUwiPZN5qzYDkJNl/OXsUQzsqr+rY0lknN444G2CpYYg6MnZmIcs7MTM7iXosbkKONzdV1WqsjnctqZqkWObo8oSbVfG3R8BHoGg92ac84hIhnN3bnp1doXVz+84eRgH9u+Uwqgat0Suwm4jSHh3Ane4+8bkhlS/zOwugvUB1xIkvAUxqi0Ot73jnKpnpbp1aSciUmsPTVnE89OWlu1fffgAThnVI04LSaQjywhghrvf2AQT3p0E6/qtB45096+qqBoZxjDUzPKqqLNPpbp1aSciUiuvf7GCO98uX9r0pJHd+dkRmTMWL1GJJL0CgllYmhQzuwP4OfA9QcL7oqq67r4U+IzginZ8jHONJRjXtwr4uK7tRERqY9riDVz3QvmvsAP6deQODU2okUSS3vvA0GQHUp/M7DbgBmAjQcKryVXW7eH29+F8o5FzdQEeDHfviDGrSqLtRESqtWD1Zi56cjqFJcGvkP5dWvPw2aNolqO57GsikWd6NwFTzexqd2/080ea2XHAr8PdhcCVVfw1NM/d74jsuPtLZvYQcCkwy8wmUT5xdFvgVYIJpCtItJ2ISHVW5hdw7sSp5BcUAdCpdTMeP28f2rXMTXFkTUciSW808Dhwt5mdQvXj9J6KVd6AOkR9PTp8xTIFuCO6wN0vM7MPgMsJVmvIBuYBE4GHqrpaS7SdiEhV8rcVMWHiVFbmbwegVbNsHj9vX3p2aJniyJqWRCacLiUYh1ajMWnurukAEqQJp0UEYHtRCec89inTFn8PQG62MfG8fThkQOcUR9Y4xZtwOpErvadoHOvkiYikveKSUq58bmZZwgP44/jhSngJSmQasvPqIQ4REanE3bnptdn856vVZWW/PnYIx4/onsKomjZ19xERaaT+NGkBz00tH3x+8Zh+XHRIvxRG1PQp6YmINEJPf/Idf36nfEj0iSO7c+NRg1MYUXqoNumZ2cVmVqfOKGaWHU6kLCIi1fjHrJXc9Nrssv0xAztz5ynDyMrS4PO6qsmV3sPAV2Y2Ic7UWjGZWZ6ZnQfMBR5KID4RkYwyef4arn5+JpGO9cN6tOOhs/YmN1s35pKhJp/iGUALgjFmq8zsUTM7w8z6xKpsZn3N7Ewzm0gw3dZjBNNynZ6ckEVE0tO0xRv46dMzKCoJMl6/zq2YeN4+tGreWFdoa3qq/STd/W9m9hrBRM2XARcA5wOY2Q5gA7CJYLaRjpQvOWQEg9Z/B9zr7tuTHr2ISJqYvTyfCx6fxvaiYO6K7u3zePrC/ejUunmKI0svNfrzIUxYvzOz3wMnAScAY4DuQLfwFbEUeJdguq3XNfuIiEh8C9ds4dyJU9m8oxiATq2b8/RF+9Gtfa2eKEkN1Oqa2d1LgBfDF2bWCegCtCOYzHmNu69PdpAiIulq6YZtnP3op2zYWghA2xY5/PXCfenbqVWKI0tPdbpR7O7rgHVJikVEJKOs2bydcx77lFWbgqc/LZtl88QF+zJkt7Ypjix91bo7kJmda2YH1qDe/mZ2bmJhiYikt43bCjnn0aksXr8NgGbZWTxyzmj27rVLiiNLb4n0gX0CuKgG9S4kWI1BRESibNpexITHpzF/9WYAsrOM+84cycEDOqU4svRXnwM/NIpSRKSSLTuKOf/xaXyxdGNZ2R9OGcYPh+6awqgyR30mvR7Alno8v4hIk7KtsJgLnpjGjO/KV0y47YQ9OWnvHimMKrPUqCNLjGdz/eM8r8sBhhCsFD6tDrGJiKSNgsISLnpyOlO/3VBWdsuP9+Cc/XunMKrMU9Pem09QcQ29g8JXVYxgJfU/JhaWiEj62F5UwsV/nc5Hi8pHdP3qmCGcf1DfFEaVmWqa9KIXjp0ALAI+rKJuIbAceM3dv6hbeCIiTduO4hIufXoG7y8oH911/VGD+MkYLRGUCjWdkeW8yNdmNgH4wN0vqK+gRETSQVFJKVc8O5N3568tK7vmiIFcNq5/CqPKbIkMTu+LOqiIiMRVVFLK1c/PrLDq+ZWH9efqIwakMCqpddJz9+/qIxARkXRRWFzKVc/N5O05q8rKLhnbj2uPHJjCqATqMA2ZmbUARhNMNt2iqnru/lSi7yEi0tTsKC7h8mdmMmlu+RXehQf35cajBmOm4cupllDSM7NrgJsJlhOqjpKeiGSE7UVBp5XoZ3g/OaQvvzxmiBJeI1HrpGdmFwB3hbtzgXkE6+mJiGSs7UUl/OSp6RV6aV46bneu/+EgJbxGJJErvasIhi+c4+7PJjkeEZEmZ1thMRc9WXEc3lWHD+CaIwYo4TUyiSS9gcBHSngiIsFcmhc8Ma3CTCvXHjmQqw5XL83GKJGktw1YkuxARESams3bizj/8WlMj5pL8/qjBmkcXiOWSNL7CNgz2YGIiDQlG7YWct7jU/lyWX5Z2a+OGaKZVhq5RFZZ+A0wOJyZRUQk46zK385pf/m4QsK7+Ud7KOE1AdVe6ZnZmBjFdwMTzewY4C2C252lsdq7+3t1ilBEpBFZsn4bZz32CUs3FABgBv97wl6cuV+vFEcmNVGT25uTqbjCQoQBp4SvqngN30NEpNGbv2oz5zz2KWs27wAgJ8u4+7QRHDe8W4ojk5qqSUJ6j9hJT0QkY3y+dCPnPT6VjduKAGiek8VDZ+/NYYO7pjgyqY1qk567j2uAOEREGq2PF63noiensbWwBIDWzXN4dMJo9u/XMcWRSW3p1qOISBzvzF3Npc98RmFx0G1hl5a5PHnBvgzr0T7FkUkilPRERKrwwrSl/OKVWZSUBk94urZtztMX7seArm1SHJkkKpG5N2P15oylEFjn7gtr+x4iIqnk7tz/34Xc9Z+vy8p6dWjJMxftR88OLVMYmdRVIld6k6lFxxYz2wQ8Cdzk7psTeD8RkQZTUurc8vpsnv6kfOKpPXZryxMX7EOXNlWuoiZNRCJJ772w3YHh/veUj9PrDXQgSIqfAp2BPsCVwDgzO9Ddt9UxZhGRerG9qISfPf95hcVfD+rfkYfPHkWbFrkpjEySJZEZWY4Kt18Bx7h7R3cf6e6j3L0TcDQwhyDx7QUMIJi6bC+CFRpERBqd/G1FnPvY1AoJ77jh3Xj8vH2V8NJIIknv1wQJ7DB3f7vyQXf/F3AkwfycN7v7YuBMYAdwcuKhiojUj5X5BYz/y0dMXVy+UsJFB/flntNG0CwnkV+T0lgl8t08DXjX3ddUVcHdVwPvAqeG+0uBzwiWJRIRaTTmrtzESSApTbEAABrvSURBVA9+xNert5SV/eqYIfz6R3uQlaW18NJNIkmvB8FVW3V2AN2j9pcCzRN4PxGRejF5/hrGP/wxK/O3A5Cbbdx7+ghNHJ3GEunIsg4YY2Z57l4Qq4KZ5QFjgPVRxbsAGxN4PxGRpPvrJ99xy2uzCYfg0bp5Dg+dvTeHDOic2sCkXiVypfcG0BV4wcx6Vj4Ylv0N6AK8HnVoMPBNIkGKiCRLSalz25tfcdOr5Qmve/s8Xr70QCW8DJDIld4tBD00jwUWmtnHwHcEvTV7EwxlyA3LbgEws1FAL+CpJMQsIpKQbYXFXPXc50yau7qsbHiPdvzfhNEag5chap303H2tmR0IPAT8mOA2ZoUqwJvApe6+Nmwzw8xy3b2krgGLiCRi9abtXPjkNGYv31RW9sOhXbnntJHkNctOYWTSkBKae9PdVwInmFkvgqQX6bCyAng/HKZQuY0SnoikxJwV+Vz05PSyDisAF4/px41HDVYPzQxTpwmn3X0J8HSSYhERSbq3vlzJdS9+QUFR8Hd3dpZx2/F7aqXzDKVVFkQkLZWWOvdM+po//7d8zvs2zXN44Ky9GTNQHVYyVbVJL7yFCbDc3Uui9mskvBpMKTMbRDB92j7AaIJB8gaMd/eXqml7JnApMAzIBuYBjwMPuXtpstuJSN1t2VHMtX/7nH9/Vd5hpU/Hljw6YTT9u2hZoExWkyu9xQSTSe8BfB3u13SVBa/he9S3S4Gra9vIzB4ALgO2A+8ARcDhwP3A4WY2PtazykTbiUjdLVm/jZ88NZ35q8sXdTlkQCfuP2Nv2rXUHJqZriYJaQlB8iqqtN+UzAb+AEwHZgCPAWPjNTCzkwkS1ypgjLsvCMu7EkyxdiJwBXBvMtqJSN19tGgdlz3zGRu3FZWVXXRwX248ejA52ZpDU2qQ9Ny9T7z9psDdH43eN6tRb61fhNsbIokrPNdqM7uUYF3BG83svkq3KxNtJyIJcnee/Ggxt701t2yV82bZWfzupL04ZVSPFEcnjUljuPXY6JhZD2AUwervL1Y+7u5TzGw5wVCN/QmWTkq4nYgkblthMTe+PIvXv1hRVta5TXP+cs4o9u61Swojk8ZI1/uxjQy3c6qaXxSYVqluXdqJSAK+WbuFEx74sELCG9ajHW9ccbASnsSU8JWemfUHLgEOIFgh/TV3vz48tj9Br8UX3L0pTjLdN9x+F6dOpFdq36iyRNuJSC29PXsV1734BVt2FJeVnbFvT2758VBa5GqGFYktoaRnZhcCDwDNwiIHOkVV6UwwTVkRQTf9pqZ1uN0ap05k8a3o/s+JthORGiouKeWP//6ah6csKitrlpPFb4/fk1P32WkOfJEKan1708wOAv5C0B3/58B+BGPeor0NbAKOq2uAKRL599S2l2qi7cpPYHaxmU03s+lr165N9DQiaWndlh2c89jUCgmvxy55/P3SA5XwpEYSudK7nuCX+tHu/jHs3BvS3YvMbD4wpM4RpkZkgE/rOHUixzZHlSXaroy7PwI8AjB69OimNjREpN588s16rn5+Jqs3la9hPW5QZ+45bQTtWzaL01KkXCJJ7wBgaiThxbGUppv0Fofb3nHqRP6sXBxVlmg7EalCSalz/38Xcu87X5etf2cGVx8+gKsOG6AJo6VWEkl67YBlNajXLMHzNwYzw+3QOCvE71Opbl3aiUgMazZt5+rnP+fjb9aXle3SMpe7TxvBoYO6pDAyaaoSGbKwhpr1PBwELE/g/Cnn7kuBzwgS9/jKx81sLNCDYNaVj+vaTkR2NuXrtRx97/sVEt6+fTvwz6vHKOFJwhJJeh8Ce5vZ6KoqmNmRBJM6T04wrsbg9nD7+3B4BgBm1gV4MNy9I8asKom2ExGgqKSU3789jwkTp7J+ayEQ3M686vABPHvRfuzaTiucS+ISuf34J4KrmL+b2UXApOiDZjYGmAgUA/fVOcIkMLO9KU84EEyeDfA7M7suUuju+0d9/ZKZPUQwWfUsM5tE+cTRbYFXCSaQriDRdiICSzds4+rnZ/LZkvLhvZ3bNOfe00ZwYP9OcVqK1Eytk567f2pm1xNM4PxPgqEJTrCS+rEE4/UMuNbdZyUz2DpoSzC0orIB8Rq5+2Vm9gFwOcEE1ZElgiYSZ4mgRNuJZCp35+XPlnPr63MqDDY/ZEAn/nTaCDq1bp7C6CSdmHtiveLN7CjgNwTr00V3n5oF3OTur9c9vMw2evRonz59eqrDEKlX328t5FevzuIfs1aVlWVnGdf9YBCXjOmn3plSa2Y2w91jPoJLuHelu78NvG1mHQk6tmQDS919RfyWIiKB9xes5boXv6gw9q5vp1b86bQRjOjZPoWRSbqq85ACd18PrK+2oohIaHtRCXe+PZ+JH35bofzM/Xrx62OH0LJZUx3tJI1dtT9ZZnZuXd7A3Z+qS3sRSS+zl+dz7Quf8/XqLWVlHVs14/cnD+OIPbqmMDLJBDX5c+oJ6rZSupKeiFBYXMr9/13AA5MXlS30CnDY4C78/uRhdG6jzipS/2qS9N6j6qQ3FlhN0DNRRCSm2cvzue7FL5i3qnzK2Ra5Wfz62D04a79eO83fK1Jfqk167j6uqmNmVgr8090vSGZQIpIeqrq626fPLtx5ynD6dmqVwugkE+lpsYjUi6qu7q7/4WDOO7CPhiJISijpiUhSbS8q4b7/LuDhKd/o6k4aHSU9EUmaDxeu41evzGLx+m1lZbq6k8ZESU9E6mzD1kJ++9ZX/P2zigur7NunA3eeMow+urqTRkJJT0QSFpkz83/f+orvtxWVlbdpkcMvjh7C6fv01NWdNCpKeiKSkG/XbeVXr8zio0UVJ2Q6dthu3PKjPejSVksASeNTkxlZxlRTZdd4ddz9vVpHJSKN1rbCYh58dxGPvPcNhSXlC4Z0b5/Hb0/Yk0MHa4FXabxqcqU3maoHpzvww/BV1XFdTYqkAXfnn7NX8ds3v2JF/vay8iyDCw/uyzVHDtScmdLo1eQndAl1m4ZMRJq4hWs2c8vrc/hwYcVbmcN7tud/T9iTPbu3S1FkIrVTkxlZ+jRAHCLSCG3eXsSf31nA4x8upjhqzF3HVs244ajBnDKqhzqqSJOiexEispOSUuelGUv547+/Zu3m8rXusgzOPaAP1xwxkHYtc1MYoUhilPREpIIPFqzjt299VWH6MAjG3P3m+KEM2a1tiiITqTslPREBYMHqzfzuH3N5d/7aCuVd2zbnl8cM4bjh3bQagjR5SnoiGW7dlh386T9f8/y0pRXmyszLzeaSsf24eEw/9cqUtKGfZJEMtWVHMY9/8C1/ee8btuwoLis3g/GjevA/PxhEVw0wlzSjpCeSYbYXlfDMp0t48N2FrN9aWOHYQf078stjhjC0m4YgSHpS0hPJEMUlpbz82TLunbSgwuBygP5dWvPLYwZz6KAuem4naU1JTyTNlZY6b81ayd3/+Zpv122tcKx7+zyuPmIAJ43sTk52VooiFGk4Snoiaaq0NJg27L7/Lthp+EGn1s244tD+nLFfL5rnZKcoQpGGp6QnkmZKSp03v1zBff9dyMI1Wyoca9Mih5+O3Z3zDuxDq+b67y+ZRz/1ImmiuKSU1z5fwQPvLuSbSrcxWzbLZsKBffjpmN01k4pkNCU9kSZuR3EJr3y2nAcnL2LJhm0VjrVunsOEA3tz4cH96NCqWYoiFGk8lPREmqj8giKe+fQ7Hv9wcYX5MSG4jXnBQX05/6A+tG+pZCcSoaQn0sQs31jAxA++5fmpS9haWFLhWPuWuVx4UF8mHNSHti10G1OkMiU9kSbiqxWbeOS9Rbzx5coK04VBMD/m+Qf15ez9e9NaHVREqqT/HSKNWHFJKZPmrubJj77j42/W73R8YNfW/OSQfhw/ojvNcjTOTqQ6SnoijdCGrYU8N3UJz3zy3U6zpwDs368Dl4zZnXGDOmsGFZFaUNITaURmL8/niY8W8/oXKygsLq1wLDvLOGrorlw8ph/De7ZPUYQiTZuSnkiKbd1RzJtfruC5qUv5fOnGnY53aNWMM/btyVn79aZb+7wURCiSPpT0RFLA3fliWT5/m7aE1z9fsVMvTIBhPdox4YA+HDtsN1rkaqowkWRQ0hNpQPnbinhl5jKen7Z0p/kwAXKzjWP32o1zD+zDyJ7t9bxOJMmU9ETqWWFxKVO+XsurM5fzn7mrd3pWB7B751acsW8vThzZnY6tm6cgSpHMoKQnUg/cnc+WbOSVmct488uVbNxWtFOdFrlZHLtXN87Ytyejeu+iqzqRBqCkJ5JEi9Zu4bXPV/DqzOU7zYMZsWf3tpy+Ty+OG9FNs6aINDAlPZE6WrB6M2/NWsk/Z61i/uqdn9NBsFjr8SO6ccLI7gzs2qaBIxSRCCU9kVpyd+av3sw/Zq3iH7NW7rRmXUSbFjkcu9dunDCyO/v26UBWlm5fiqSakp5IDRSVlDLju+95Z+5q3pm7Zqf16iKa52QxblBnjh/RncMGd9FQA5FGRklPpAobtxUy5eu1TJq7hinz17Bpe3HMenm52Rw2uAtH77Urhw7qohXJRRox/e8UCZWWOnNWbOL9hWuZPH8tM777fqfVDCJaNcvmsCFdOWbPXRk3qAt5zXRFJ9IUKOlJRlu6YRsfLFzHBwvW8eGidTGHFkR0a9eCw4Z04fAhXTmgX0fduhRpgpT0JKOszC9g6rcbmPrtBj5YuI7v1sceVgBgBsN7tOeIIV04bHBXhuzWRmPpRJo4JT1JW+7OorVbmbZ4A9O+3cDUxRtY9n1B3DadWjfn4P4dOXhAZ8YO7EznNpodRSSdKOlJ2sgvKGLWsny+WLaRL5ZuZMZ337N+a2HcNnm52ezXrwMH9+/EwQM6MairruZE0pmSXj0yszOBS4FhQDYwD3gceMjdd56AUWqsoLCEr1Zu4sswwX25LL/KYQTR8nKzGdmrPfv06cD+/Tqyd+/2NM/RszmRTKGkV0/M7AHgMmA78A5QBBwO3A8cbmbj3X3n9WSkAndn2fcFzFu1mXkrNzFv1WbmrtrE4nVbqaJjZQXtW+YyuncH9u27C/v06cCe3duRm51V/4GLSKOkpFcPzOxkgoS3Chjj7gvC8q7Au8CJwBXAvSkLspEpKillyYZtfLN2K9+u28I3a7eycM0W5q3azJYdscfHVZaTZQzerQ3De7RneI/2jOjVnv6dW2smFBEpo6RXP34Rbm+IJDwAd19tZpcCk4Ebzey+TLrNuWVHMcu/L2D5xm0s+76ApWGS+2bdVpZs2FblmLhYsgz6dGrFiB7tGdajHcN6tmeP3dpqGIGIxKWkl2Rm1gMYBRQCL1Y+7u5TzGw50B3YH/ioYSNMvtJSZ2NBEWs372DN5u2s3byDtZt3sHrTjrIEt3xjQdwxcPG0b5nLkF3bMni3NmXbAV3aaEC4iNSakl7yjQy3c9y9qv7x0wiS3kgaQdIrLXV2FJdSUFQSvAqLKSgsZcuOYvILithUUER+QREbCwrJLygivyAo37itsCzBFdfiKq0qu7VrQb/OrejXqTX9Oreib6dWDNmtLV3aNFePShFJCiW95Osbbr+LU2dJpbpJddGT01m7ZQelpU5xqVNa6pS4U1Ja8VVYUkpBYZDoGkKz7Cy675JH9/bBq8cuefTp1KoswbVsph9HEalf+i2TfK3Dbbz+85G1aOplYbXZy/NZtWl7fZy6Sm1b5NC5TXM6t2lOlzYtyr7uFia3Hu3z6NS6uTqViEhKKeklX+S3ekL3+8zsYuBigF69eiUUQHYCiaV5ThZ5zbLJy80u27Zslk27vGa0y8ulXV4u7Vvmln3dLvy6c+sguakDiYg0BUp6yRdZOrt1nDqRYzsts+3ujwCPAIwePTqhxPnw2aMoLi0lO8vIMiMn28g2IyvLyAnLsrOM3OwsWoYJTldgIpIJlPSSb3G47R2nTs9KdZNqrx7t6uO0IiJNnqamSL6Z4XaomeVVUWefSnVFRKQBKOklmbsvBT4DmgHjKx83s7FAD4LZWj5u2OhERDKbkl79uD3c/t7M+kcKzawL8GC4e0cmzcYiItIY6JlePXD3l8zsIYIVFmaZ2STKJ5xuC7xKMPG0iIg0ICW9euLul5nZB8DlwFjKlxaaiJYWEhFJCSW9euTuzwLPpjoOEREJ6JmeiIhkDHOv+0TBUj/MbC3x5/CsTidgXZLCkcZH39/0pu9v4nq7e+dYB5T00piZTXf30amOQ+qHvr/pTd/f+qHbmyIikjGU9EREJGMo6aW3R1IdgNQrfX/Tm76/9UDP9EREJGPoSk9ERDKGkl6aMbMzzex9M8s3sy1mNt3MLjczfa+bMDMbZGZXm9nTZjbPzErNzM3slFTHJnVnZrlmdriZ3WVmn5jZSjMrNLPlZvaSmY1LdYzpQrc304iZPQBcBmwH3qF8vs82wCvAeHcvSV2Ekigzuwe4Osah8e7+UkPHI8llZkcA/wl3VwEzgK3AHsCeYflt7n5zCsJLK/rrP02Y2ckECW8VMMzdf+TuJwIDgLnAicAVKQxR6mY28AfgNKA/MCW14UiSlQIvA2Pcfbfw/+9p7r4XcDpQAtxkZoemNMo0oCu9NGFm04FRwAR3f6rSsbHAZIKE2F2TXTd9ZjaZYCJzXellADN7FLgQmOjuF6Y6nqZMV3ppwMx6ECS8QuDFysfdfQqwHNgV2L9hoxORJJgZbnukNIo0oKSXHkaG2znuXlBFnWmV6opI0zEg3K5MaRRpQEkvPfQNt/Emp15Sqa6INAFmtitwXrj7cgpDSQtKeumhdbjdGqfOlnDbpp5jEZEkMbMc4GmgHfCOu7+R4pCaPCW99GDhVr2SRNLLwwTDjpYCZ6c4lrSgpJceNofb1nHqRI5tjlNHRBoJM7uXoMfmKuBwd1+V4pDSgpJeelgcbnvHqdOzUl0RaaTM7C7gKmAtQcJbkOKQ0oaSXnqIdGceamZ5VdTZp1JdEWmEzOxO4FpgPXCku3+V4pDSipJeGnD3pcBnQDNgfOXj4eD0HgS3ST5u2OhEpKbM7A7g58D3BAnvixSHlHaU9NLH7eH292bWP1JoZl2AB8PdOzQbi0jjZGa3ATcAGwkSnu7K1ANNQ5ZGzOxB4FKCCacnUT7hdFvgVeAUTTjdNJnZ3pT/8QLBRMRtgAXAhkihu2vGnSbIzI4DXgt3pwNzqqg6z93vaJio0pOSXpoxszOBy4G9gGxgHjAReEhXeU1XuLTMu9XVc3erro40PmZ2HvB4DapOcfdx9RtNelPSExGRjKFneiIikjGU9EREJGMo6YmISMZQ0hMRkYyhpCciIhlDSU9ERDKGkp6IiGQMJT2RJs7MxpmZm9nkVMdSE2Z2QxjvUXU4x95mVmpmf0xmbJL+lPREmgAzWxwmij6pjqUuzGw34FfAe+7+dqLncffPgL8DV5nZgGTFJ+lPSU+k6ZsKDAHOTXUgNfAbgjlDf5Okc+VSPtm6SLU0DZlIE2BmiwkWCe7r7otTG01izKwjsAxYAfT3JPzyMbNpwEign7svqev5JP3pSk+kETOz88zMCRIewLfhbc7Iq09Vz/TCYx7eGs0ys2vNbI6ZFZjZMjO728xahnV3MbN7wro7zGyBmV0bJy4zs9PN7N9mti5ss8TM/i/OLdgLgBbAU7ESnpm1N7PfhTFui4pzspn9oopzPkkwsfolcT5GkTI5qQ5AROJaSPCL/RSgFfAysCXq+JZYjWJ4FvgRMDk85xjgGmCImZ0FfEJw2/EDoEN4/C4za+Huv4s+kZnlAs8DJwEFBEvhrAb2BC4CTjazH7j79EoxnBBuJ1UOLky+HxIsmbQmrLMV2C0s25/YtzEj5zqe4FmhSFy6vSnSBMS7vRm17FCFZWfCK65vw935wGHuviI81hOYCXQEZhMsQXWOu28Pjx8LvAlsBnZ1921R572DYLHT94Cz3H1Z1LErgPuARcBgdy8Oy1sSLI4K0DbyPlHtziVI7m8BJ0TahceygbHu/t8Yn4sB64FdwjhXx/4ERQK6vSmSGa6KJDwAd18KPB3u9gYujU5E7v4W8CXB1d/oSLmZdQCuIrjCHB+d8MJ29xMkrt2Bo6MODSXodPJt5YQX6hpuJ0UnvPCcJbESXnjMgbnh7ohYdUSiKemJpL8iIFbSWBhup7v7uhjHF4TbblFlhwJ5BFeVa6p4vynh9oCosi7hdn0VbaaG2xvM7Gwza19FvVgiK8d3jVtLBD3TE8kEqypfPYUizwOXxTgWfbxFVFm/cHts2MEmns5RX7cLt5tiVXT3KWZ2J3Ad8FfAzWwewTPGl939X3HeJ3LO2iRKyVBKeiLpr7SOx6Nlh9v5BJ1f4vk06uuy53lVVXb3G8zsYYJOKQcDBwE/AX5iZv8Gjq0ieUfO+X018Ygo6YlIrSwNt7Pc/bxatIvcCu0Yr5K7fwvcE74ws4OB54AfEAx5eCRGs8g5q7rdKlJGz/REmobCcJvqP1QnETwjPKKWz93mADuAvmaWV9NG7v4B8ES4O7zy8bD35uBwd2Yt4pEMpaQn0jQsD7dDUhlEOCTgAYLnZ6+b2eDKdcKB7heZWdeodgUEtztzgVEx2pxoZmPMLKtSeR5wRLj7XYyQBhMMV5gTp2ONSJlU/9UoIjXzCjAOeCZ8vhV5RnZDCmK5nqBH56nAbDP7nGA8YAugJ0FibhZuo8fNvUow6P0Igg4q0cYCVwNrzWwmsJag88uBBIPl5wF/iRFLJCG+Vud/lWQEJT2RpuF+gg4bZxHMrNI8LP9tQwfi7kXAaWb2DMFztn2BYQQD2VcSzP7yGsEA9WhPAP8LnGtmv6k0FdkTwHaCDix7Ap0IEvtCgmd6j7n75hjhTABKiJ0QRXaiGVlEpMGEvTMvAQ6vasB5Lc61F8EA+pfd/ZRkxCfpT0lPRBqMme0KfA3MdPexdTzXS8BxwFB3X1BdfRFQRxYRaUDuvorgluyYuq6cTjDh9X1KeFIbutITEZGMoSs9ERHJGEp6IiKSMZT0REQkYyjpiYhIxlDSExGRjKGkJyIiGUNJT0REMoaSnoiIZAwlPRERyRhKeiIikjGU9EREJGMo6YmISMZQ0hMRkYyhpCciIhlDSU9ERDKGkp6IiGQMJT0REckYSnoiIpIxlPRERCRjKOmJiEjGUNITEZGMoaQnIiIZQ0lPREQyhpKeiIhkDCU9ERHJGEp6IiKSMZT0REQkYyjpiYhIxlDSExGRjKGkJyIiGeP/AzfrOZ4Q+qHcAAAAAElFTkSuQmCC\n", | |
"text/plain": [ | |
"<Figure size 432x288 with 1 Axes>" | |
] | |
}, | |
"metadata": { | |
"needs_background": "light" | |
}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"plt.plot(t_plot, height_predicted,'-')\n", | |
"plt.plot(t_plot[-1],height_desired, '*')\n", | |
"plt.xlabel('time(s)\\n')\n", | |
"plt.ylabel('Height(m)\\n')\n", | |
"plt.title('Height vs Time\\n')" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"## References\n", | |
"\n", | |
"1. Math 24 _Rocket Motion_. <https://www.math24.net/rocket-motion/\\>\n", | |
"\n", | |
"2. Kasdin and Paley. _Engineering Dynamics_. [ch 6-Linear Momentum of a Multiparticle System pp234-235](https://www.jstor.org/stable/j.ctvcm4ggj.9) Princeton University Press \n", | |
"\n", | |
"3. <https://en.wikipedia.org/wiki/Specific_impulse>\n", | |
"\n", | |
"4. <https://www.apogeerockets.com/Rocket_Motors/Estes_Motors/13mm_Motors/Estes_13mm_1_4A3-3T>" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": {}, | |
"outputs": [], | |
"source": [] | |
}, | |
{ | |
"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 | |
} |