Skip to content
Permalink
master
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
Go to file
 
 
Cannot retrieve contributors at this time
{
"cells": [
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"#importing necessary libraries\n",
"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\n",
"\n",
" \n",
"#defining a simplerocket function dependent on the state of the rocet and returing its time derivitives\n",
"def simplerocket(state,dmdt=0.05, u=250):\n",
" derivs = np.array([state[1], (u/state[2])*dmdt, -dmdt])\n",
" return derivs\n",
"\n",
"#using Euler method as an explicit integration method\n",
"def eulerstep(state, rhs, dt):\n",
" next_state = state + rhs(state) * dt\n",
" return next_state\n",
"\n",
"#recalling runge-kutta method\n",
"def rk2_step(state, rhs, dt):\n",
" mid_state = state + rhs(state) * dt*0.5 \n",
" next_state = state + rhs(mid_state)*dt\n",
" return next_state\n",
"\n",
"##using Heun method as an implicit integration method\n",
"def heun_step(state,rhs,dt,etol=0.000001,maxiters = 100):\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": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0. , 0.1025, 0.205 , 0.3075, 0.41 , 0.5125, 0.615 , 0.7175,\n",
" 0.82 , 0.9225, 1.025 , 1.1275, 1.23 , 1.3325, 1.435 , 1.5375,\n",
" 1.64 , 1.7425, 1.845 , 1.9475, 2.05 , 2.1525, 2.255 , 2.3575,\n",
" 2.46 , 2.5625, 2.665 , 2.7675, 2.87 , 2.9725, 3.075 , 3.1775,\n",
" 3.28 , 3.3825, 3.485 , 3.5875, 3.69 , 3.7925, 3.895 , 3.9975,\n",
" 4.1 ])"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = np.linspace(0,4.1,41) #10 seconds, 10 time steps per second\n",
"dt = 0.1 #time step of 0.1 seconds\n",
"N = 41\n",
"t"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"#setting initial conditions\n",
"y0 = 0 # initial position\n",
"v0 = 0 # initial velocity\n",
"m0 = 0.25 #initial mass ig kg\n",
"\n",
"#initialize solution array for euler method\n",
"num_sol_euler = np.zeros([N,3])\n",
"\n",
"#Set intial conditions\n",
"num_sol_euler[0,0] = y0\n",
"num_sol_euler[0,1] = v0\n",
"num_sol_euler[0,2] = m0\n",
"\n",
"#initialize solution array for heun method\n",
"num_sol_heun = np.zeros([N,3])\n",
"\n",
"#Set intial conditions\n",
"num_sol_heun[0,0] = y0\n",
"num_sol_heun[0,1] = v0\n",
"num_sol_heun[0,2] = m0"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
"#integrating the simplerocket equation using the eulerstep and heun integration method\n",
"for i in range(N-1):\n",
" num_sol_euler[i+1] = eulerstep(num_sol_euler[i], simplerocket, dt)\n",
" num_sol_heun[i+1] = heun_step(num_sol_heun[i], simplerocket, dt)"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"#from the integrated array above, the mass reaches 0.05kg after 4.1 seconds\n",
"u = 250\n",
"v_final_euler = num_sol_euler[40,1] #grabbing the final velocity from the integrated solution\n",
"v_initial_euler = v0\n",
"v_final_heun = num_sol_heun[40,1] #grabbing the final velocity from the integrated solution\n",
"v_initial_heun = v0"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [],
"source": [
"#developing velocity and mass variables for euler method\n",
"final_mass_euler = num_sol_euler[:,2]\n",
"initial_mass_euler = m0\n",
"current_mass_euler = final_mass_euler/initial_mass_euler\n",
"current_velocity_euler = num_sol_euler[:,1]\n"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([1. , 0.98, 0.96, 0.94, 0.92, 0.9 , 0.88, 0.86, 0.84, 0.82, 0.8 ,\n",
" 0.78, 0.76, 0.74, 0.72, 0.7 , 0.68, 0.66, 0.64, 0.62, 0.6 , 0.58,\n",
" 0.56, 0.54, 0.52, 0.5 , 0.48, 0.46, 0.44, 0.42, 0.4 , 0.38, 0.36,\n",
" 0.34, 0.32, 0.3 , 0.28, 0.26, 0.24, 0.22, 0.2 ])"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"current_mass_euler"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAcsAAAE0CAYAAABDxhiLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOydd3xUxfbAv4dUAgmhBghI6BClh6IPNSBNpIgCCSASQJ5PEXn4rBRFfQ/92csTFXxSBKRXlSII2EC69Ca9BULoJIEk8/vj3t1sNrub3RQ2CfP9fO5n9945M/fceu7MnDkjSik0Go1Go9E4p5i3FdBoNBqNpqCjjaVGo9FoNNmgjaVGo9FoNNmgjaVGo9FoNNmgjaVGo9FoNNmgjaVGo9FoNNlQqI2liNQWkfdFZIuIXBSRGyJyWkS2ichcEXlWRBqJiDjIu0ZElIhM9oLqeYaIHDGPY6y3dfEWIjLZPAdrvK1LbhCRWuZxWJbjIpLtMyoi79rlG30r9NW4j4jUsLk+n3mY91MzX7qIVM9DnTaZ5f43r8rMC0TkPVOvnd7WxZZCayxFZASwC3gOaAKUAvyAikAj4FHgY2AbUNZLamqKOCLia/MSfCyPi68CtMlm/8WAvnm8X00eo5Q6BPxmrsaIiJ87+UTEF4gxV39TSh3OD/0KCyJyl83zFnUr910ojaX5UvoAwzgeBUYATYEKQDhwLzAK2OItHTWaXHLZ/O2fjVw7oLKNvKbgMtX8LQt0djNPJ6C8XX6NFyiUxhL4j/l7GGiilPpIKbVVKXVOKXVKKfWrUmqcUqoZ0Aq4Zl+AUipaKSVKqbhbqLdG4y5zzd9HRSTIhZzFmM7JZ300uWcWkGz+d7cVwnJ9k4HZea5RAUQp9bz5br7L27rYUuiMpYjUAe4wV79SSl1wJa+U+kMplZT/mmk0ecoa4BhQEnjYkYCIlAR6mKvf3Bq1NDlFKXUJWGKudhWRUFfyIhICdDNXF5v5NV6i0BlLoJzN/ys5LcSVg4+IjDXTjpjrdUXkfyJyTESSReQvEXnH9mYXkUARGWE6G10WkUsiskJE7nFXBxHpLiI/isg5EUkSkV0i8mo2NQt3jtVXRAaJyDIROWM6Qp0z9evnyAHKzBdn6R8w16uKyMcisl9ErptpoXZ5gkXkFRH5Q0QumOfrmIhMF5G/ualvIxH5XET2mOfymojsE5H5ItLX0/NhXpuFpr4pItLbiVx7EflWRI6ael8UkQ0i8pKjfYrIr8BNm03f2DnaKBFp7YmuNihguvn/cScyjwAlANv+MKeISEMRGSUia0XkrIjcNO/TP8VwEqqSTf7yIvJvEdls5rspIvEislNEppj3ko+DfA1EZKJ5Da+b5/a4iGwUkY9ExGW/rIPytpjn9ns3ZKdanmX7+1xEHhKRBaYuKSJyRUQOmefnVRGp54lebmJpSg0AemUj2wsItMuXBRG503xe9orIVfN52S2G82Pl3CgrIqHmudhkPg/J5rmcKiLN3Syjuc31v2LquFdE5ohIbxEJtJN36OAjIgnADptNGx08b+VEpKPNust3johUF8NxSonIIJcHopQqVAtQH+NFooBFuShnjVnGZAdpY820I0AH4KrNPm2XzUAIUAb43YnMDaBjdjpgNC07yq8wHJnCnJRxxJQZ6yT9DgwnJ2dlK+AHoISDvHE2Mq2ARAd5Q23kGwAnstnXu4A40dUHeA9Iz6aMh+3yTTa3r3FQZilgrZl+BWjvQKY4RhOZq30eAGrY5fs1mzwKaO3BPVnLJt9jZNzrqUBFB/IrLdce8LXJO9qBbDM3dL0EtHWi213AWTfKCLXL18/U31WebR4+uyPMfDeB8i7kgsxrroD/2KV97saxfJQP7y9fIN4sf62b76h4wNeJzBggLZtr2s5J3k2mzH+dpEcBZ1yUnQ68ns2xjnfjPLezy/eeuX2n3fYEN8oqh1EJPGquf5XNOX7dlLsKlHQpm9c3Q34vgJD5hfwlUDcH5VhuxMkO0saaaRcxDMQmMjraI4B/2+z/DWCBeVP+E6iO0YHfDThpyhxzdLPb6HDY/F0G/M3MfyeGN6/FcPyMAyODC2OJYcgPmukJGJ7D9YBQoDbwMnDdTJ/mIH+czXGeMI/jcQwvzTCgCxBoypYFTpmy14FXMF7+5YC2ZDYsLzi5Jh/ZyPyB4dFcBeNj5E7g7xiGr7tdvsk4MJYYntGWD4VzQHMn+51Pxsv3A4yXRFlz34+bx62A3UCQTb7iGMbYovNgjGZT26WYB/dkJmNpbrO80J6zkw0n4yVZi+yNZVNgHfAChgNcLfO81sMwaJvNvOeBCg7yryPjxf13oI6ZvyZwP/AqsJfMH09lyPjQ3IjRnFzd3B4JPAh8Bvzo4bNbkQwDPMyFXB+bc1LfZntHm+3fAveZ57M8hmd9b4w+4Lfy6R1muc/TgWpOZO4g49l3aLSBl2yOYy7Gc1bBPI4HMZ4hhfHBUNtBfqfGEqhEhnG6AvzLvHblMCoQG2z2/ZQT/b6ykfkF6G6e5zIYH9ZPY1Qy3DWWJYDmNmXei93zZiM71pS5jM0za1eekPH+nJTtdcuPmyG/F4yH2/6L4qh5g78E3AP4ZFPGGrI3lgrjJVLcgcw0Ml6wN4FWDmTa2ZSTpXZpo4MCljrSGeMlZJF51EG65WKPdZD2ic0N4/CDwk7HKLu0OJu0BKCKi/Np+wJwdKz+5gOjMJwVKtilt7LZ11ycfEmbsr5265OxM5YYL/G/yPhYqeekrN42++3tRKYKhrFVwD/tdbHJ/1gu72tHxnK4ub7VTtbyovzdgR5ZjKUb+/YD1pv5x9illbYp+yEPyuxh5rmBXY0ztwvGh6UC/nAh850ps8lu+8fm9g15qZMHuje1OZ8jnci8YiPT1EF6DYz3jgLed1KGv801nekg3ZWxtBi6NOBeB+nFbfJfBkrZpdu+V6bi4qPRwfPs0FiaaXfZlBvlosxqZHxM9nci84BNWVmOMYu8N26WPLrhYsioyThaTgIvAv5O8q/BPWPprAmji43MN05khIyX7KsudFDYfPnayfjZHOd3DtKP4MBYYnyFXTPTns/mXK4y5T602x5no99LLvL7YNTCFTDfhVwTm/Lsa0rfmtsv2T94btwLk7ExlkBjMpqPduPayP9myi3MZh+W5pqNdtvz21iWJ+OleJeN7E5z2z8c6OGxsTTLeNbM/7Pd9vI2ZTf2oLxeZp5EnDS95+JcPWajUy0H6eVszpv9B85n7lzz/Fww+t4UsMdJ+m6cGAwz/X0z/TiuPywfNOVSsPvox4mxxDCESWbaVBdl32tzDf5ul7bE3H4WJzU7F+Xm2liasstNuVVO0r8x0w+4o1dhdPABQCk1C+PrKgbjoA/ZiVQG/g9Y5alDiA0pGM1+jvjL5v9yJzoqG70qutjPXqXUHidl3CTDg86ps5AD7sHoswFYKyIlnS3AdlPO1SBfV84UDTCaI8HFEAal1FYyztu9dsltzd+FKhdefyJyP8Y1C8NohrpXKXXCiWxJoKW5ujqbc7TLlGvsyIklv1BKnSPj/upv6t0Uo1n6Bh4MJxCDR03HikOmI4jVOQKjxgVQ14EOJ83Vz0WksZu7tDSBlwYm5tbZxI4FZAwJczQMIwbjAyIN40PMlq3mbxcRGW5e31uNxXu5ntgNrheRZhj91eDcsaed+bsWCHRx31reK/5AQzd1a0aGY5Gr5/kXjI9SsHmezefjfkt+pdR1N/eb1/zP/G0jItVsE8TwNH7EXJ3kTmGF1lgCKKWSlVKzlVKPK6VqYrSFd8dw1kg3xVoDb+VwF+dMY+UI2+Eop12UYZEr7kJmbzZ6WG740uZFdgfbF94GjH4HZ8s/TbnyOMf+Y8QW2xtxdzZ6WYyONY+IBGP0tYDxgs0p9TGa50IwDMwDSqnzLuRrYdSKwWhGdnWOZplyvhgv/1uJ5cXaT4yIPZaxd98rpRLdKUBESgArMJq4e2L0Pzn7iCzlYJvFqaYVsFVEDpsesE+ISISjQpRSBzBqcWD0554Qka1ihG/rJSI5Po9KqWsYBhOMbhl7LAZ0pVIq3i5tCkbAEh+M654gIqtF5E0xPKL9c6qXB0wj4x1lb+wt1zedDI9oe+qYv/1wfd8etsnj6vm2xZPn2ZJum6cCEGz+z83znFsWYvTBCzDALi0G4/5Px7gfsqVQG0t7lFIXlFKLlVKxGF8NykwaIiIBOSgyLQ/lHA7PMLmaTV7b9GCnUplx9MLLjkBnCdl8HdrqlN2xWIb72OYJcZCeE4LIOIYrGC0DrsjJOQIX5ymfWITRPB2O4VzRx9zuSUSXj8joR5oEPITxsVAW41oEA8NMWV/7zEqpORh9PKswXjARGM5PE4HDYgyDctQy8SyGQ9AujGegMfAMRo34jDkEIcyD47DF8hFRS0RaWTaKSA0Mo24rY3ssN4FoYBxGzSjAXB+N8UFxRkRey0+jqZQ6hXEuAfqIEdbOEt4u1tz+k1LqpH1e812Wk3vQ3TwF5XnOFUqpGxgfJQADRDINHYozf1c4OseOKFLG0hal1CKMIRFg1OryY8xUXpFdM5Bturs3XyYDq4yIGNktER5p7Vgnd4/FNo/tf3c/BhyxGePlDEbtaYblJeQE23PUyc1zJM6adfMLpVQyGRF9PsNoYk4k4/52idkcZxmrOU4pNUgp9YNS6i+lVKJS6qpS6iquWz9QSq1WSrUjI1zbfzDOORjNbr+aTYi2eZRSaqIyorHcgWEIPsNwyPPHqEX9brYueMoqMlp1bGtnlv/XMGoXjo7lilJqFEZ3TUPgSWAGxkdJaQy/BWe1urzC8rFTAeMjCPM3zC49E0qpFDLG94714L6d66g8BxSU5zkv+Mr8rYHh9WwJbGPp0vra3YKKrLE02WXzP1cD+/OZ7Ay5pf/iglLK3Rigts2m7vYx5ZQjNv8js5G90z6PeUxnzdVc6aqU+pQMg9kL1wbzMBmtD01ys99bgKWGVMP8nWV+ObtDJIZhgqz9d7Y0cKcwpdRFpdRSpdRopVQURo0zGaOG9oqLfMeVUrOUUs9gNAO/bCbVIGszmTt62PZHxthcZ0tg+QVmc62rMpRSaodSaoJSqh9G7d3SvNtTRPIz5Np8MgzLY3a/18x0Z1iaV/Pjvj1i8z+759mSbpsnnozjyu93j0uUUjsxuqEgozY50PxNxGi1cYuibixtI5Kc8poW2VNPnEQLEWN2gq7m6u8elLmWjGbIga4E84CdGF/kYIyNdIiINMJo+gNj3KUtK83fhz3ol3WIA4M53ZHBNPv7LDWjx8WN6bAcYHFPh4z+z/zgZ4zamAVPmmBtuyAc6iiZQ6t5hFLqJ+Anc7W+K1mbPAp4h4zavVv5HGBpZisHdDSjytS1S3Mb07ja+jjkVC939nWdDIP4sOkA1d1cn5eNoV9h/nYUkUp5rNpmMnwtXD3Pf8MYjwk2z7P5EbPaXO0lIi5bLDzE1ofE3efN4ujT07zPLX3C0z344Cx8xlJEaooRcqtMNnKNyfB2OqCUOupKvgDwoZOX9Stk3JBueW2BtbZmuUniRMTpTQ/GyzKnD535cEw2Vx8RkXb2MqbR/8RcTSZrX5IlLQTDc9Lpg5BN06pFJ1uD2RvDYDoq8wPztz7Gy9spIuIjIjXt9mMZGgFGk16+YO7nHgw96yml1nuQ3baVobt9otmX81+c9OGKEebO6fNmntcIc/W8zfYa2fT7VcIY4pQpnyeYHtaWFqTHyHD2OUPGB5i9vtm15Nhe4xzp5QGWj57iGM+EpQUsu1i/n2B8qAUCU7MzSCJS11W6LWaz/wxz9TFxELJTjBB1H5qrV4CZDvQDw6lovKsPUXeeZxtsr4e7z9u3GDX1ksAEjNYD8KAJFih84yzJGGeTZJ6E/hhNAWUxvi6jMKLqXCZjPE4vB+WsIftxlkdc6BFhU360CzlX+7GkWZoDl2K8EC0RTj4idxF8SmF42iqznP9h9C2FYfTL1Mb4cvwKo2bY0y5/nOUY3bguthF8rmKMca1hbm9j6m85Xy86KeNDG5nfMT52wk1d6wODMGowbkXwMdOG2ZQ5C8eBH2bbyKzCiDRTxTx/d2D0I72NUbPLEk2FjPFc+zEcS0piOMk4Hf/m5PizjLP0IG92EXwsEZRSMBxZ6mI8L/eb953CMDoKSLXL2w7jZTMdw7monnmPhmMM+fneZt9P2eT7N0af4kcYDkURGNGjqmN4I1rGEqYCDXLxTnjZLOc6GaHkPnAh/yvGOMdRGP1Ylc17rB5GEAjLmOHj2I3TJnP0LrdDGbrQRcgIzaZs9ptt5CfgeZs8e8zno5Z5jitjjAR4CcPzd72D/O5G8LFEJ4sw75n2ZEQHUsDTTvSbaCOzBqPlwnKuLRG5fsHNCD426Zbz9TOGPShONs8bRkXD9hxv8fha5fZi3+oF4wWfYnfgzpbrOA/FtIaCYywnYzT9ODuO3MSGrYR7MUwV0M0ub5wlzc1r405s2PdwHRv2Izf0dDs2rJluazBnYmcwMZopJ7qxXwX8n4PyO7uQz3FsWA+fi+yM5Z04ju1rWaYDQ3BuLN05N5OwecmT2bA4W25iN6A9B++EqmSNJ5wl6o2NvDvPwzmghYO8eWoszTLH2e37bQ/y/hNjvG12x5MlDi3uxYaNd1FmdrFh/TBqctnp5qmxfM5FWeWc5GltJ+c0TKLT48mLi32rF4ymut4YQXrXYTiH3MBo3juD0V4+Bqjqoow15kmb7CBtrJl2xEX+CJsTH53D/WRKw6hJrcL4okvC+PJ+FRcRMMjGWJoyglFbmo3xVZaE8cFxCqOm9hKOo6DEWY7Rg2sTDIzE6FS/aO7nGEazzt/cLCMKoxZ8EOOD5wpGDXkuRo3EPhLJZFwYS1PGpcG02e+XGF/plzFe5OcxvqLfwag1OjP0D2LU0M6SETmmwBhLm3v2fxgfNDcwXoQrgb5m+hM4Npb+GAbzLYyv+cPmPZSE0cT7LY4D1JfG6DP+AuPFfNLc7xWMmt2nOIlc5eli3seW49+djWwExofBTFOPcxi120SM98kYoIyTvPlhLOuS+UUe6WH+OzBaPjaZx5CKURvcgfER2BkHkczIxliaMqHAaxj9mJcw3rFHMJqPHcZadlDGPRhjGQ+Z98wl8xmbhdGyFWAn79JYmjJx5r2YSOZA8g6NpZnH0sqW7Oz6ulrELETjBURkDUYz2BSlJ6HWaDSafENEtgGNgNlKqRhP8xc6Bx+NRqPRaDxBRJpgGErw1LHHRBtLjUaj0RR1hpu/R4Afc1KAJy67Go1Go9EUCswhKUEYHtyWsZXvKaXSnedyjjaWGo1GoylSmJGXdtht3orhnZuzMrWDT+4oV66cioiIyFHeffv2cfXqVcqWLUtOy9BoNJrCyObNmxOUUu7OhOIRNsZSYXj9/wCMUsZ0czlC1yxzSUREBJs2bfK2GhqNRlOoEJF8i6qmjJiwrmZ68hjt4KPRaDQaTTZoY6nRaDQaTTYUWGMpIuNERJnL8y7k+orILyJySUSuisgmERma3QwSOc2n0Wg0mtuPAmkYzGl2XsTonHUl9xlGTMsojIC8PwJ1MGZQmOts5oqc5tNoNBrN7UmBM5YiEoAR6zMeFxNzmlNOPY0RC7ahUqqLUqoHRqD1PUAP4Jm8yqfRaDSa25cCZywxpteKBP5BxoTCjrDMyP6SUuqAZaNSKh54ylx92UGzak7zaTQajeY2pUANHRGRlsC/gBlKqSXOJiwWkSpAM4wZDObYpyul1orISYz59lphzI2Y43wajS2XL1/m7Nmz3Lx5M3thjeY2xM/PjwoVKhASEuJtVfKMAmMszZm3p2BMuTI8G/Em5u8upVSSE5mNGEavCRlGL6f5NBrAMJTx8fGEh4dTvHhxRPJ0KJdGU+hRSpGUlMTJkycBiozBLEhNjf/BmNdtmFIqIRvZ6uavq0Gtx+xkc5NPowHg7NmzhIeHExQUpA2lRuMAESEoKIjw8HDOnj3rbXXyjAJhLEXkHowZvxcqpWa5kaWk+XvNhcxV8zc4D/JlQkT+bg412XTuXI6jJxnEx8OMGXD9eu7K0dwSbt68SfHixb2thkZT4ClevHiR6qrwurEUkeLAJIyZ6Z92N5v562lg25zmy4RSaoJSKkopFVW+fM5DG6rff4cvv4T9++HHHM0ao/ECukap0WRPUXtOvG4sgXEYYxyfU0qddjPPFfO3pAsZS9oVm205zZfnnD17lrk//cSVS6bD744dcMmV869Go9FovEVBcPDpAaQDA0RkgF1aPfP3KRHpAhxUSj2BMYEnQDUX5VY1f4/YbMtpvjxl586dLFy4kNTUVEIuXqRd/fr4dO8OpUrl1y41Go1GkwsKQs0SDD3ud7CEmek1zPUoc32r+Xun2YzriOZ2srnJl6dUqlTJ2kSxoXJlfrrjDshFc65G4wljx45FRBwu06ZN86isuLg4oqKishfMBUePHqV///7ccccdBAYGUrVqVbp3787PP//sUTmTJ09GRLh69Wr2wjbMnj2byZMnZ9keHR1Nz549PSpLU3jxes1SKRXhLE1EJgMDgBeUUu/Z5DkuIluApkAvYKpdvvuBKhhRetblNl9eU7ZsWdq3b88PP/xAmo8Pv69bR7369alatWr2mTWaPKBUqVIsW7Ysy/ZatWp5QRvnXLhwgVatWlGpUiXeeustKleuzJEjR1i8eDHr1q3jvvvuy3cdZs+eTUJCAnFxcZm2jx8/Hj8/v3zfv6Zg4HVjmQvewggs8H8i8rtS6iCAiFQAxpsybyul0vMoX57SvHlz9u7dy6FDh1BKsWDBAv7xj3/g7+8PKSng4wO+hfnyaAoyvr6+tGrVyttqWElKSnLoZTx37lzi4+P5888/qVChgnX7wIED8fbE9ZGRkV7dv+bWUlCaYT1GKTUX+ByoCOwQkSUiMh84gBEubyFGYPQ8yZfXiAjdu3cnICAAgMTERFauXAkHD8L48eBhE5NGk5esWbMGEWHnzp2ZtrvT9Hjs2DFiY2MpU6YMQUFBdOzYkX379lnTjxw5gogwffp0Hn/8cUJDQ+natavDsi5evIi/vz9lypTJkmbvbTl79mwaNGhAQEAAVatWZdSoUaSmpubqGOPi4pg3bx5r1661NlWPHTvW6bn46aefaNmyJYGBgYSFhfH0009nava17HPNmjX06tWLkiVLUqNGDcaPH4+mYFNojSWAUuppoB+wBaNPsyNwECMQ+qNKqbS8zJfXlCpVigcffNC6fmD5ci58+qnhFfvrr3Dq1K1QQ3ObkpqammXJLYmJibRu3Zp9+/bxxRdfMHv2bK5du0a7du1ISsocNOv5558nODiYOXPmMHLkSIflNW3alJSUFPr378/mzZtJT3fc4LNixQpiYmJo2rQpixYtYtiwYbz33ns880zu5kQYM2YMbdq0oUmTJqxbt45169bxxBNPOJTdvXs3nTp1oly5csybN4/XX3+dGTNmOPy4GDJkCI0aNWLBggVER0czdOhQNmzYkCtdNflLgW7nU0rFAXHZyMwAZuSg7Bzly2saNWrEnj172LdvHxdCQ9mwfz8PBAfjGxKiAxUUEiw1jcKkw/nz5x32tx0+fJiIiIgc6/Hhhx9y7do1tm3bZq0N/u1vfyMiIoKvv/6aoUOHWmVbtWrFZ5995rK8Bx54gBEjRvDRRx8xc+ZMgoODad++PU899RTt2rWzyr366qtER0czZcoUADp16gTAK6+8wujRo6lSpUqOjqdmzZqUKVOG9PT0bJut33jjDapVq8bixYvx8TFm+StTpgwxMTGsW7eOu+++2yrbp08fRo8eDRg11CVLljB//nxatGiRIz01+U+hrlkWBUSErl27EhQUBCJsveMOtty8CUOHQgFzttAUHUqVKsXGjRuzLJUrV85VuStXrqR9+/aEhIRYa6vBwcE0a9aMTZs2ZZJ96KGH3Crzgw8+YP/+/bz77rtER0ezbNkyOnTowBdffAFAWloaW7ZsoVevXpnyxcTEkJ6ezrp1+earl4kNGzbQo0cPq6EEePTRR/H19eXXX3/NJNuhQwfrfz8/P2rXrs2JEyduiZ6anKGNZQGgZMmSdOnSBYDk4sX5oXhx9ukHR5OP+Pr6EhUVlWXx9/fPVbkJCQnMmjULPz+/TMvq1as5fvx4JtmwsDAnpWSlVq1aPP/88yxevJijR4/SuHFjRo4ciVKKhIQEbt68maU8y3piYmKujsldTp8+nUUHHx8fypYtm0WH0NDQTOv+/v4kJyfnu46anFOgm2FvJyIjI2nQoAE7duwAYPHixQwdOtSocWoKNAWhGTavCQwMBODGjRuZticmJlKuXDmn+cqUKUO3bt0YM2ZMlrTg4MzhlnMaDq1cuXIMHDiQZ599lrNnz1KuXDn8/PyyBO2Oj4+36uSInB6jMypVqpRFh7S0NM6fP+9UB03hQdcsCxCdO3e2vlCuXbvGd999Z7jHX7wIv+vZwjS3Dksf3549e6zbjh8/nsmr1REPPPAAu3bt4s4778xSa61bt67HejibqODAgQMEBARQqlQpfHx8aNasGXPmZJ6idvbs2RQrVixTX6Et7h6ju7W+li1bsmDBAtLSMvwD58+fT2pqKq1bt842v6Zgo2uWBYjixYvTrVs3pk+fDhjedfvnzqXuwYPG2MvQUNBjuzR5QGpqKuvXr8+yvWrVqoSHh1OlShWaN2/OmDFjCAoKIj09nXHjxmVbQ3ruueeYNm0abdu2ZdiwYYSHhxMfH8/atWtp3bo1ffr08UjPKVOmWIeYNGrUiJs3b7Jq1SrGjx/PU089Za0dvv7663Ts2JGBAwcSGxvLjh07GDNmDEOGDHHq3OPuMdarV49FixaxcOFCqlSpQuXKlR327Y4ePZomTZrw8MMP89RTT3HixAleeuklOnbs6NRgawoPumZZwKhduzbNmjWzru9cuJDrFy4YK99/D0VoyhuN97h06RJ33313lmXSpElWmRkzZnDHHXfw2GOPMXLkSF599dVsa4flypVj/fr11KtXjxEjRtChQwdefPFFLl26RMOGDT3Ws3Pnztxzzz1MnDiR7t27ExMTw6pVq/j000/58MMPrXIdOnRg5syZbNq0ia5du/LRRx/xr3/9i//+155bs0QAACAASURBVPWQaXeO8emnn6ZDhw4MGjSI5s2bM2HCBIdl3XnnnSxdupSzZ8/yyCOPMHr0aPr06cPcuXM9Pm5NwUO8HQWjsBMVFaXsvfxyS0pKChMmTOD8+fP4pqbSZt8+Wt1zDz69ekEOXeA1ecOePXuoX7++t9XQaAoFrp4XEdmslMrfwMJ5iK5ZFkACAgLo1asXvr6+pPr68ltEBCuqV9eGUqPRaLyENpYFlIoVK9KxY0cArpcowR/btrFr1y4va6XRaDS3J9pYFmCioqIyBWtevHhxxnitGzdAj8vSaDSaW4I2lgUYEaFbt26ULl0aMPoy586dS+rJkzBhAixeDLrPWaPRaPIdbSwLOIGBgfTs2dMaQivx0CEOjRwJCQmwezds3uxlDTUajaboo41lISA8PJz27dsDkBwYyKabN0lISAA/PzCn+NJoNBpN/qGNZSGhZcuW1vFfB2vV4peLF7kUEwMNGnhZM41Goyn6aGNZSBARHn74YUqVKkV6sWL8WasWc9asyZM5CDUajUbjGqfh7kRkfx7tQymlPA8KqclC8eLF6dmzJ5MmTSI9PZ0TJ06wdOlSp7PMazQajSZvcFWzrJWHiyaPqFq1aqZJbzdv3mzME3j9OsyYAYcPe1E7TWFg7NixTmfViIuLIyqq4AVV2blzJw8//DCVKlWiePHiVK9endjYWHbu3OlROa6O3RUTJkxg4cKFWbZHRETw/PPPe1yepvCRXSD1ecDLuSj//4AeucivccDdd9/N6dOnrdN5/TxnDjUDAyktAidPwpNPQkiIl7XUaPKGgwcP0qpVK1q0aMF///tfSpcuzYEDB5gzZw7bt2/nrrvuyncdJkyYwF133cXDDz+cafuCBQsoW7Zsvu9f432yM5ZXlFJ/5bRwEbmS07wa51jGX547d44zZ86QUqwY2zdvplWDBgQAHDwITZt6W02NJk+YNGkSAQEBLF26lADT+7tt27Y8+eSTeDu2dZMmTby6f82tw1Uz7PfAtlyWvw34IZdlaBzg5+dHbGwsQUFBpAQEsLVWLbbt309qr17aUGrylGPHjhEbG0uZMmUICgqiY8eOmeZ8XLNmDSKSpUk0Ojqanj17WtctTbw//vgjDRs2pESJErRu3TrbMI4XL14kNDTUaihtsZ9A+r///S+1a9cmICCAWrVqZZqZxBGTJ09GRLh69Wqm7bbNq9HR0WzevJkpU6YgIogIkydPziJnYfbs2TRo0ICAgACqVq3KqFGjMjniWfa5Y8cO2rdvT4kSJahXrx7z5893qavGuzg1lkqprkqpT3JTuFLqY6WU9j7JJ0JDQ+nduzfFihXjYmgoK+rX57sDB7z+ta0pHKSmpmZZ7O+dxMREWrduzb59+/jiiy+YPXs2165do127diQlJXm8z2PHjvHCCy8watQovv32W86ePUvv3r1d3rNNmzbl0KFDDB8+nN27dzuVmzhxIsOGDaNbt24sWbKEXr168a9//Yu3337bYz1tGT9+PPXq1aNz586sW7eOdevW8dBDDzmUXbFiBTExMTRt2pRFixYxbNgw3nvvPZ555pkssn379qVbt24sWLCA2rVrExsby4kTJ3Klqyb/0JM/F3IiIiLo2LEjS5cuJc3Xl23btlGpUiVatmzpbdVuK9asMRZ3aNYM7B2YlyxxPxhTdLSx5Ibz58/j5+fnMM12PtUPP/yQa9eusW3bNuukyH/729+IiIjg66+/ZujQoR7tNzExkd9++43atWsDkJ6eTo8ePdi3bx/16tVzmGfAgAGsWLGCTz75hE8++YQyZcrQuXNnhg8fbnVGSk9PZ+zYscTFxfH+++8DxhyXly5d4q233uKf//yndaJoT4mMjKREiRKUL1+eVq1auZR99dVXiY6OZsqUKQB06tQJgFdeeYXRo0dnmoh6xIgRDBo0CDDOeVhYGN999x3/+Mc/cqSnJn/R4yyLAC1atKBx48bW9eXLl3P48GFIS4NVq+DSJS9qpymIlCpVio0bN2ZZunTpkklu5cqVtG/fnpCQEGvtMzg4mGbNmpGTeVwjIiKshhKwThTgqkbl6+vLrFmz+PPPP3nzzTdp1qwZs2fP5u677+b777+35j916hS9evXKlDcmJobLly9bneHyk7S0NLZs2eJQh/T0dNatW5dpe4cOHaz/y5YtS4UKFXTNsgDjds1SRKoBzYCNSqnjNtsbAJ8CjYAjwEtKqRV5rKfGBSJCly5dOHfuHCdPniQ9PZ2F06fzZJkyBJ09azj8DBwI/v7eVlVTQPD19XU4RKRs2bKcPn3aup6QkMD69euZNWtWFtkHHnjA4/2GhoZmWvc378lkN2bQadiwIQ0bNgTgyJEj3HfffYwePZqHHnrIqnNYWFimPJZ162w9+UhCQgI3b950WwdH58Kd86DxDp40wz4PPA1YAwyISDCwEihvbmoELBKRRkqpvApqoHEDX19fYmJimDBhguGskJDAzvXradq4Mb6nT8PWraCbZvON3DaNdu2atWm2IFCmTBm6devGmDFjsqQFBwcDWJs3b9y4kSk9MTExR2Ma3SEiIoJevXoxfvx4ACpVqgTA2bNnM8nFx8cDWJuQ7XGm+4ULFzzWqVy5cvj5+Xmsg6Zw4Ekz7H3AXqXUQZttj2EYytlAPeBFIAB4Ns801LhNSEgIvXv3xsfHh0uhoWytXJldu3aRdv/90KKFt9XTFEIeeOABdu3axZ133klUVFSmxRKr2NIPt2fPHmu+48ePZ/KYzQ32xsfCgQMHrLW2KlWqULlyZebMmZNJZvbs2YSEhNDASQxlR7r/8ccfXL58OZOcO7U+Hx8fmjVr5lCHYsWKcffdd7vMrynYeFKzrARstNvWEUgH/qmUOgO8JyJxQJu8UU/jKXfccYfVw+505cr8GBLCkcuX6QZItrk1msw899xzTJs2jbZt2zJs2DDCw8OJj49n7dq1tG7dmj59+lClShWaN2/OmDFjCAoKIj09nXHjxuVZTerNN9/kzz//pG/fvtSvX59r164xf/58lixZwnvvvQdAsWLFGDt2LE8++SRly5alffv2rF27ls8//5xx48Y5de5p0aIF4eHhPPvss7z55pskJibyzjvvEGIX1KNevXosX76c5cuXU7ZsWapXr+4wGMHrr79Ox44dGThwILGxsezYsYMxY8YwZMiQTM49msKHJ8YyFLBvm7gb2GEaSgu7MIyoxks0atSIixcvsnr1aq6WLMnWrVspXbo09913n7dV0xQyypUrx/r16xk1ahQjRozg4sWLVKpUidatW1v7DwFmzJjBE088wWOPPUaVKlV45513sh3j6C79+vXj6tWrvP/++5w8eZKgoCDq1KnDt99+S2xsrFVuyJAhpKSk8NFHH/Hxxx9TpUoV3n//fUaMGOG0bH9/fxYsWMDTTz9Nz549qVu3Lp9//jn9+vXLJDd69GiOHTtG7969uXz5MpMmTSIuLi5LeR06dGDmzJn8+9//Zvr06VSoUIF//etfvP7663lyLjTeQ9wdkyci54FtSqkHzPU6wF7gc6XUUBu5b4GHlFK3Rby1qKgolROvwPxGKcWiRYvYti0jrsQjjzxivOA2boSaNUH3oXjMnj17qF+/vrfV0GgKBa6eFxHZrJQqeIGIneBJzXI78DcRqa6UOgwMBhSwxk4uAjiDxquICF27duXSpUvGMBJg0YIFhG3eTNjRo1C2LAweDEFBXtZUo9FoCj6eOPhMBPyBLSKyAcM7NgH4ziIgIiWBJhhNsRov4+PjQ0xMDOXLG87KgVevsu/bb7l+/TqcPw+//OJlDTUajaZw4LaxVErNAMYBgUAUcBLopZSyjXnVC8OgrslDHTW5IDAwkH79+lGyZEmulSjBjtq12b59O8k1a0IOxslpNBrN7YhHEXyUUqOBMkBlpdQdSqmf7URWA82BSXmknyYPCA0NpW/fvvj5+XGufHnW16vHN0lJ3EhP97ZqGo1GUyhwaixF5GXTiScTSqkkO+9X27QjSqnNSqnLjtI13qNy5cr07NkTEeFiaCgnT51i9uzZpKWleVs1jUajKfC4qlmOA/aIyC4ReVNE9LxPhZy6devy4IMPWtcPHjzIggULSE9PhwMH4BbEz9RoNJrCiCtv2BjgEeBBYBQwUkSOA/OBBcCvSs8FVeho0aIF165dY+3atQDs3LmTCqdOce+FC8bcgEFBxrASjUaj0VhxNZ/lHKVUH4xwdl2ByUAQ8E8MB57TIvKliHQSEcdz/WgKJNHR0bQww99JejpXli/n8F9/QXo6rFhh/Go0Go3GSrYOPkqpm0qp75VSg4GKGKHsPgNSgCHA98A5EZkmIo+KiB64V8ARER588EEaNmyIKlaM7Q0bsj8+nr+uXoXHHoNieuY2jUajscVTb9h0pdRapdSzSqlqQEvgXYwgBH0xAqqfE5EFIvK4iIS6Kk/jPUSE7t27U6dOHVICAvizcWO+DQxky4ED3lZNo9FoChy5qkIopTYqpV5WStUD7gLGAvuB7hjDR/TsIwUYHx8fevXqRbVq1UgqXpxUX1+WLFnC7t27va2aJp8QkWyXNWvWZFvOyy+/7HFg8OTkZESEr776Kofa5x3Lli1DRDh48GD2wh5y5coVXnnlFerUqUNgYCBhYWG0adOGqVOnelTO3r17ERFWrlzpUb7ff/+df//731m25+SaaTLwJNydS5RSu4HdwJsiUh3ogRG4QFOA8fPzo0+fPkyZMoXTp0+jlGLevHkEBARQMyTE6MPs0UOHxSsirFu3zvo/KSmJtm3bWidQthAZGZltOUOHDiUmJiZfdCzsdO3alf379zN69GgiIyM5c+YMa9asYdmyZTz++OP5vn+LsRw9enSm7fqa5Y48M5a2mLFjP8iPsjV5T2BgII899hiTJk0iISGBtLQ0vpswgf4ilAkMhGnT4PHHwck0R5rCQ6tWraz/r169CkDNmjUzbXeHqlWrUrVq1TzVrSiwY8cO1q5dy+LFi+lqM5t3bGws3h48oK9Z7vC4GVZE/EWklYg8IiJ9nS35oawm/yhRogT9+/enVKlSAPhdvszOjRu5ePEiJCQYsWQ1tw3nz58nLi6OSpUqERgYSLVq1Rg61Dq5kMMmvQMHDtC1a1eCg4MJCQmhR48e1iD+zti6dSvly5fniSeeMMb7ulFOy5YtHdbQnnnmGWrXrg0Ys+688cYb1KhRg8DAQCpWrEjnzp057+I+njJlCv7+/kyaNInNmzcjIvzxxx+ZZC5evEjx4sWZMGGCwzIuXrwIQMWKFbOkiWSeUXbTpk1ER0cTFBRE2bJlGTBgAAkJCU71c9aMbXstvvjiC1544QVSUlKszeqdOnXKImchu3Nt2ecXX3zBiy++SNmyZQkLC2P48OHcvHnTqa5FEY9qliLyPDASKOWG+IwcaaTxGqVKlWLAgAFMmjSJs8CetDRS9+yhTp8+VAkP97Z6mlvIsGHD2L59O5988gkVKlTg2LFjmZpw7bE06YaEhPD1118DxhyQ0dHRbN++3foRZsuGDRvo2LEj/fr149NPP0VE3ConNjaWsWPHkpycbJ3UOT09nXnz5jF48GAAJk6cyPvvv88777xD/fr1OXfuHCtXriQpKSmLHgBffvklw4YNY8qUKfTp0weAxo0bM2nSJFq2bGmV+/bbbxERp82ZkZGRBAYG8swzz/Cf//yHe++9l4CAgCxyp0+fpk2bNjRu3JiZM2dy4cIFXnrpJXbt2sX69evx9c1Zo98jjzzC7t27+fLLL61jqUNDHftZenLNxo0bR8eOHfn222/ZvHkzo0aNombNmjz77G3klqKUcmvBcNZJN5edGIEJvnG2uFtuYV+aNWumihoJCQnqvffeU6+99pr698iR6j//+Y86evSot9UqEOzevdtxwurVSr32mrGsXp01fdmyjPTffsuavnhxRvqmTVnT587NSN++PSeqZ+HKlSsKUJMmTcqSVrNmTTVhwgSneV966SUVHh5uXf/www+Vn5+fOnbsmHXbX3/9pXx8fNQHH3yglFIqKSlJAWrixInql19+UcHBwer555/PVK475Zw4cUKJiFqwYIFV5qefflKA2rFjh1JKqcGDB6u+ffs61X/p0qUKUAcOHFAfffSR8vf3V/Pmzcsk8+mnn6pSpUqppKQk67bmzZurfv36OS1XKaUmT56sgoKCFKD8/f3V/fffr/73v/+p9PR0q8zw4cNV2bJl1dWrV63b1q5dqwA1f/58pZRSe/bsUYD68ccfs5w/W+yvxbvvvqsCAgKy6JWba9a+fftMZXXs2FHdf//9Ls+DUi6eF6UUsEkVgHe4u4snzbBDgVSgu1LqLqVUD6VUf2dLbo24xntYmoRKlizJTT8/bty4wbRp0zh27JghkJICqaneVVKTrzRu3Ji33nqLL774wi2P0Q0bNtCqVatMfWI1atSgefPm/Prrr5lkV69eTceOHRk+fDjvvvuux+WEh4fTunVrZs2aZZWZNWsWkZGR3HXXXVb9Fy5cyBtvvMGmTZusTbz2vPPOO7z88svMnz+fRx55JFNav379SElJYcGCBQDs3r2bjRs3MnDgQJfnYsCAARw5coSJEyfSq1cv9uzZw+DBgxk0aFCm4+zcuTMlSpSwbrvvvvuoWLFilvOVX3hyzTp06JBpPTIykhMnTtwSPQsKnhjLasBapdSS/FJGU3AoV64ccXFxlCxZEoAbN24wffp0jh84AFOnwuzZ2mAWYSZMmECnTp149dVXqV27NvXq1WP+/PlO5U+fPk1YWFiW7WFhYSQmJmbatmzZMooVK0b//lm/qd0tJzY2liVLlnD9+nVSU1OZN28esbGx1vSnnnqK1157jenTp9O8eXMqVqzI66+/nsVozps3j3r16tGmTZss+yxdujQ9evRg0iRjEqVJkyZRrVo12rZt6/Q8WLD0w06bNo3jx4/Tr18/Jk+ezN69ez06zvzEEx3sm3L9/f1JTk7OV/0KGp4YyzOA9vK4jShXrhwDBgywfv3evH6dnS++yOU9e2D/fsNg6vDABtHRMHassURHZ03v2DEj/Z57sqZ37ZqR3qxZ1vRHH81Ib9Agr7R2SpkyZRg/fjzx8fFs3bqVRo0a0bt3b6e1zEqVKnH27Nks2+Pj4ylTpkymbW+88QYtW7akXbt2Ga0VHpbTs2dPkpOT+f777/npp59ISEjI1I/o4+PDiy++yL59+zhy5AjDhg3j9ddfzzLWcdasWSQkJNCjRw9u3LiRZb9PPPEEq1at4vDhw0ybNo24uLgsjjrZ4e/vz/DhwwHYt2+fR8dpi5+fH8WKFcuiZ06Na050uJ3xxFguAlrrOLC3F+XLlycuLo4SJUqQXqwYZ4OD+fPPP7l8+TLUrQsevjg0hQsRoXHjxrz99tukpaWxf/9+h3ItW7Zk3bp1nDyZMbT6yJEjbNy4kdatW2eSDQgIYNGiRVStWpV27dpx5swZj8upUKECbdq0YdasWcyaNYsmTZpQp06WGQUBqFatGmPGjKFq1apZAm5ERESwcuVKtm3bRmxsbJYp69q0aUO1atV4/PHHiY+PZ8CAAS7P1+XLl0lJScmy/YAZGctSk2vZsiU//PAD169ft8r88ssvnDlzJsv5suDj40OlSpXYs2ePdVtqaiqrV6/OJOfv78/NmzedNj1b8OSaaTwzlmOBZGCyDmN3e1G+fHkGDBhAUIkSHK5enUNVqjA3OZnD+uuzyNKyZUs+/PBDVqxYwfLlyxkxYgQhISE0c1TrBf7+978TFhZGp06dmDt3LnPmzOHBBx8kPDzc6qFqS4kSJfjhhx8oWbIkHTp0sNaOPCknJiaGH374gfnz52fxTh04cCCjR49m8eLFrFmzhpEjR3L8+HGHza1169blxx9/ZM2aNQwcODDTeEgRYdCgQfz6669ER0dTvXp1l+dt+/bt1KpVi1dffZVly5axevVq3n//fZ555hlatGhhncDghRdeIDk5mQcffJAlS5YwdepUYmJiaNasWabxmfZYmoUnTJjA0qVLefTRR7MY53r16pGens7HH3/Mxo0brYbaHk+v2W2PJ95AQFlgO5AILAO+AiY4WL70tufSrVqKojesM86cOaPeeecd9dprr6nXXntNvfnmm2rfvn3eVuuW4sq7r7Dhyht2+PDh6s4771QlSpRQoaGhqm3btmrdunXWdHvPSqWU2r9/v3rooYdUiRIlVMmSJVX37t3VoUOHrOmOvDnPnTunIiMjVYsWLdTly5fdKsdCYmKi8vPzU4A6fPhwprQJEyaoVq1aqdDQUBUUFKQaNWqkpkyZYk239Ya1sH79ehUcHKyeeuqpTGXt2LFDAeqbb75xcTYNEhIS1KhRo1RUVJQqXbq0CgoKUvXq1VMjR45UFy5cyCS7YcMGdd9996nAwEBVunRp1b9/f3Xu3Dlrur03rFJKXbx4UfXt21eFhoaqihUrqrfffjvLtUhLS1PDhw9XYWFhSkRUx44dlVJ5d82cleWIouQNK8rNPicRCQRmAV2A7NrelFLKxwObXWiJiopSmzZt8rYat4yEhASmTJnClStXAChWrBiPPvood0ZGwu+/Q1QUOBhXVlTYs2cP9evX97YamlvIBx98wBtvvMGpU6cI0mEfPcLV8yIim5VSUbdYpRzjycjXNzHmtbyAEXDgIHA1P5TSFFzKlSvHoEGDmDp1KhcuXCA9PZ25c+YQGBpKzYsXYfduY5qv4sW9rapGkysOHTrE/v37eeedd3jiiSe0obzN8cRYxgAXgcZKqeP5pI+mEFC6dGkGDhzI1KlTSUhIoNSFCxxfs4bA2rUJB9i0Ce6919tqajS5YuTIkSxYsIC2bdvy6quvelsdjZfxxFiWA1ZoQ6kBCAkJYeDAgXzzzTecAfbXrg0HDnCpWjUitSedpggwc+ZMb6ugKUB44g17iOz7KjW3ESVKlCAuLo6qVatyKjycPxs1Ys7Nm6xctQp3+8I1Go2mMOCJsZwMRItIhXzSRVMICQwMpH///lSvXp0LpUujRPj1119ZvHhxxjiva9e8q2Qeoz8ENJrsKWrPiSfG8gNgBbBKRO7PJ300hRB/f3/69etH3bp1rdu2bt3KzJkzubl+PXzyCRw54j0F8xA/Pz+nM1doNJoMkpKS8PMrOjFsPDGWe4HGQCTwk4hcF5GDIrLfwbIvf9TVFFR8fX3p3bs3jRs3tm5L/O03to8bx82rV40JpIuAwaxQoQInT57k+vXrRe7LWaPJC5RSXL9+nZMnT1KhQtFpiPTEwaeW3XogUMOJrH6L3Ib4+PjQvXt3QkJC+Pnnn0kKCiIhOZktW7ZwV4cOlKhUydsq5pqQkBAATp06ddtNfqvRuIufnx9hYWHW56Uo4ImxrJ1vWmiKDCJC27ZtKVmyJEuXLmVrkybUPnCA9UlJ9E5MpFIRMZhF6SWg0Wiyx+0IPhrH3G4RfDxh9+7dzJ8/n1RzKq+AgABiYmKoUcNZg4RGo7ldKGwRfDzps9RoPCIyMpL+/fsTGBgIQEpKCtOnT2f79u1G/+XMmeBgWiSNRqMpaGhjqclXqlWrxqBBg6zNlmlpaSz/5hsOjRuH2rMHJk0CM86sRqPRFFScGksRmS8iz+SmcBEZJiLOp1fX3BZUqFCBwYMHU758eQBKJyZybP9+9uzZQ9qlS5DNvHsajUbjbVzVLB8Gmuay/KZA91yWoSkClCpVisGDB1OzZk1OVK3Kvjp1OJWYyCwfH675euJnptFoNLee7N5SQSJSORfl6zD9GiuBgYH07duXpUuXsmnTJs6VL0/qlSt89dVX9O3b11rz1Gg0moKGU29YEUknj8ZLFuW5LbU3rOcopfjjjz9Yvny5dWB/QEAAvXv3pmbp0vDnn3D//VBMd6lrNEWVwuYN66pmeQodXECTD4gIrVq1onTp0sybN48bN26QkpLCzMmTeSwlhWpBQRAfD488Av7+3lZXo9FonBtLpVSVW6mI5vajbt26DBo0iBkzZnD58mUqHz/O4SNHuBEeTs1ixSh2/jwUgSAGGo2m8KPbuTRepWLFigwZMoTKlStzrFo1jletysmTJ1kswjUdJUej0RQQtLHUeJ3g4GDi4uKIvOsu/qpZk03NmrEtLY0JEyZw+vRpb6un0Wg02lhqCgb+/v707NmTdu3aWWuUly5d4n//+58R8Sc5GX7+GdLSvKypRqO5HdED3DQFBhGhdevWhIWFMW/ePJKTk0lNTWXBvHmQkkKDgADk8GHo1QuC9KgkjUZz69A1S02Bo3bt2gwZMsQ67rLC2bMkrl/P9u3bubl/Pxw75mUNNRrN7YY2lpoCSdmyZXniiSeoW7cu8RUqcLh6dS5cuMDixEROaccfjUZzi9HGUlNgCQgIIDY2lug2bTharRpbmzRhZ4UKfP3112zevBk9vZxGo7lVaGOpKdCICNHR0cTGxpJcoQJKhNTUVJYsWcLChQu5mZICy5fD5cveVlWj0RRh3J78WUSeAKYrpZLyV6XChQ53d+s4f/48s2fPJj4+3rqt6YULtAsIIKh8ecPxJyLCewpqNBq3KWzh7jypWU4ATojI+yJSO78U0micYenHbNKkCQDFk5II3r6dzZs3c+7IEe34o9Fo8g1PjOV3QAgwAtgjIstEpKuISP6optFkxc/Pj+7du9OtWzduBgfzZ8OGJBUrxs+nT7P8+nXS9DhMjUaTD7htLJVS3YAawNtAAtABWAgcFpGXRSTH8yuJiJ+IPGDWWteLyGkRuSEiJ0VkrohEZ5O/r4j8IiKXROSqiGwSkaEi4vL4cppP432aNm3K4MGDkRo12BwVxd769Vm3fj2TJk3iwoUL3lZPo9EUMdzus8yUScQP6A08DdyNMTvJDWAuMF4ptc7D8toBP5qrZ4DNwDUgErjL3P6mUupVB3k/M/VIBlYBN4EHgGBgAdBLKZWlupHTfPboPkvvkpyczIIFC9i3b591W0BAAN26XeAEKAAAIABJREFUdePOY8cgNBRatQLdAKLRFCgKW59ljoxlpgJEGgFDgT5kTPa8DfgMwyEoxY0y2mIYro+VUr/YpcUA0wEfoK1SarVN2qMYBvoMcJ9S6oC5PQxYDdQH/qmU+tiuzBzlc4Q2lt5HKcW6detYuXIl6enpgBHI4IELF6hVqxY+d94JPXuCrw5YpdEUFAqbscx1c6NS6k/gdWASIObSBJgIHBGRwW6U8ZNSqqe9oTTTZgGTzdXH7JJfMX9fshg8M0888JS5+rKDZtWc5tMUQESEe+65h8GDB1O6dGlQikqnTnH69Gm2bNnCpWvXwKfIzj+u0WhuAbkyBiLSTkTmA4cxapfJwNcYtcwfgArABBF5Npd6bjV/rXNsikgVoBlG8+8c+wxKqbXASaAi0Cq3+TQFn/DwcJ588knuatCAHQ0bciI8nIT0dMYfP87mLVt0EAONRpNjPDaWIlJKRP4pInuB5cDDwClgJFBFKfWEUmqWUqorcA9G32NujaVlqIrtfE1NzN9dLsZ+brSTzU0+TSEgMDCQRx99lC4PP8zRyEg2NWtGCrBkyRLmzJlDUlISKAXXrnlbVY1GU4hwuxNHRJpi9CvGAsUxmlvXAJ8Ci5RS6fZ5lFJ/iMj3wKM5VVBEKgJx5uo8m6Tq5u9RF9ktA++q22zLaT5NIUFEaNq0KVWrVmXOnDmcPXsWgN27d3P8+HFia9QgfP9+6N4d6tXzsrYajaYw4InHg8WL5TrwFfCpUmqnG/muebgfKyLiC0wDSgGrlFJLbJJL2pTvjKvmb3Ae5NMUMsqXL8+QIUNYvnw5Fies9NOn2ff99yRVrkyNa9coFhsLkZFe1lSj0RR0PGmGPQK8gNHU+qSbhhJgCODnqWImX2AM5zhOVucey1gATzuicpovowCRv5tjMjedO3cup8VobgF+fn506dKFPn36EBQUhE96Ojf8/Dhx4gS/HDjAGT2DiUajcQNPjGVNpdT7SqmLnuxAGXgcVkVEPgYGYwzveEApdcZO5Ir5WxLnWNKu2GzLaT4rSqkJSqkopVSUZc5FTcGmbt26PP3004Q1a8amqCjiw8LYUKUKEydN4vfff9fOPxqNxiWeGMvlIvJcdkIiMkJEVuRCJ0TkfQynoHMYhvKAA7Ej5m81F0VVtZPNTT5NIadkyZL07duXTj16cLBhQ66VKEFaWhorVqxg6tSpXLp0CbZsgcREb6uq0WgKGJ70JbYDTrghF4nRdJojROQd4DngPNBeKbXbiahlOMmdIlLciWdrczvZ3OTTFAFEhKioKCIiIpg/fz6nTp0C4PDhw0wbN45HLl2iYtWqSKdO0LSpjvyj0WiA/JnP0h/I4hnrDiLyNka/6AUMQ/mnM1ml1HFgi7m/Xg7Kuh9jXOYZYF1u82mKFuXKlWPw4MHcd999iAiiFNW3b2ff3r3s2LyZpM2bva2iRqMpQOSpsTRnIGmGEWjd07xvAi8BFzEMpTu1urfM3/8TkVo2ZVUAxpurbzsY1pLTfJoihI+PD23btmXQoEGUKVeOPZGRXA8K4uzly3x+6hR/bt+u+zI1Gg2QTWxYu77HdhjBB5w1i/piBA+oDMxVSsW4rYRIN2CRuboJ2OVEdK9S6m27vOMxQtQlAyvJCIgegjErSk8ngdRzlM8eHRu2aHDz5k1WrVrFxt9/J+jqVa6YXrJ169alS5cuBAcHQ3IyBAZ6WVONpmhQ2GLDZmcsbWtWioxhF67YDnRXSrka9G+/nziM2LLZsVYpFe0gf1+McHsNMAKu78UIu/e5q9phTvPZoo1l0eLo0aMsXLgw0zRfxYsXp0edOtQ+cADp0kWPy9Ro8oCiZiwtjjoCrMAIb/eeE/EbwEml1KE81bCAo41l0ePGjRusXLmSDRs2AOB/4wbNN26kYkgIderUIaBnT2jWzMtaajSFm8JmLF16wyqlVln+i8hvGDW7VS6yaDSFHn9/fzp37kz9+vVZtGgRaceOkV6sGOfPn+eX7dsp3bUrTZVCtKesRnPbkOv5LG93dM2yaJOSksLKlSvZum4dNf/6i3Ply5NYpgzVqlWjW7dulC1b1tsqajSFksJWs9TGMpdoY3l7cPToURYvXsz58+et23x9fYmOjuZuHx98/P31uEyNxgOKjLEUkZHm38+VUhds1t1CKTUut8oVBrSxvH1ITU1l7dq1/Pbbb6SnG/5fxa9fp83evdSvXZvgBg2gVy8oUcLLmmo0BZ+iZCzTMTxg6yul9tusZ1smRkjY22Jqem0sbz/OnDnDokWLOH36NPX27KFifDz/3969x0dV3okf/3xzJSGBcAmEuwECIahcE7nITVEQtagVtVdtrW2xq213tbu/3d/uq93ddtu9te7a2nWr1V3XrUVFrQKi0XAxBBIhilwDIYSEAAkhgZCEQPLdP84JhDS3SWYyM5nv+/Wa1yHneeaZ7zxk5pvnnOc8B2Do9deT8rOfER0T4+cIjQl8wZYsO5rg8xOc5FjR6mdjQlpSUhKPPPII27ZtY5MIDdHRjC4pYUNUFBufeYbbbruN1NRUmwBkTB9i5yx7yEaWoa2yspK3336bkv37aYiKurx/0qRJrFixgoTiYpg0yRYzMKaVYBtZWrLsIUuWRlX57LPP2LBhA+fPX7mn+JCaGlZWVjI6NZWwFStg6lSbAGSMK9iSpSd3HTHGtEFEuO6665g4cSKZmZnk5eUhTU2M37uXwtpaTp48SXL//gy99lp/h2qM6aYuL6QuIqtFpEFEbu+gzh1unW94JzxjgkdMTAx33HEHDz/8MMOSkihKTuZCVBRn6+v5z5IS1q5dS01Njb/DNMZ0Q5cPw4rIezhrqI5sb91UEQkHSoF8VV3utSgDmB2GNW1pbGwkJyeHLe+9R7+qKqoSEgCIjo5myZIlZKSnE1ZaCmPH+jlSY/wj2A7DepIsj+Hc9eOWTuq9B0xW1ZD4FrBkaTpSXV3Nu+++y969V9+sJ62hgVtra0mYOxeWLQM3mRoTKoItWXpyzjIR2NSFeqeABd0Lx5i+ZeDAgdx3330cPnyYdevWcfr0aSIuXSIhN5f8ixcZXlbGuNhYYu+809+hGmM64MnNn6uBMV2oNwqwEzPGtDBhwgRWr17N0qVLiYyIoNJdU7a4qop/27WLrVu3cunSJT9HaYxpjyeHYTcAS4A0VT3cTp0JwD5gs6ou9VqUAcwOwxpPVVdXs3HjRkqyswlrbOTM4MEADBo0iGXLljF5/Hjk7FkYOtTPkRrjO8F2GNaTZPkA8DKwB7hHVQtalU8E1gJpwIOq+pKXYw1IlixNdxUWFrJ+/XrKy8uv2j+3vp75TU3ELV0KCxfaggamTwq2ZOnJOctXgC8DK4A9IrIV2O+WTcY5TxkBbAiVRGlMT4wfP55vf/vb5OXlkZWVRV1dHf3q64ncsYO8piZGlpQwJj6emLlz/R2qMSGvy+cs1RmC3gM84+5aDHzbfSxx9z0D3O3F+Izp08LDw7nhhht47LHHyMjIIFyVc/HxABw4d46nsrLYtm0bjY2Nfo7UmNDWreXuRCQJuBkY5+46CmSq6gkvxhYU7DCs8aZTp06xYf16arZvp75fP84OGADA4MGDWbp0KVNGj0YuXQL3PKcxwSrYDsPa2rA9ZMnSeJuqcuDAATZu3EhlZeVVZQvKy5kdE8PAW2+FBQvAbgdmglSwJUtbG9aYACMipKamMnHiRHJzc9m8eTN1dXUMOHuW8D172AUMKy5mVFISA6+/3t/hGhMSPLnOEgARmSwivxSRPSJS5T72iMjTIpLqiyCNCUURERHMnTuXxx9/nLlz50JEBNXuYdk9TU3825tvsnHjRurq6vwcqTF9n0eHYUXkIZxJPFFAW/caagC+paoveiW6IGCHYU1vqays5P333qN80yZq4uOpcw/B9uvXjwULFpAxZgyRYWEwerSfIzWmc8F2GNaT6yzTgWyc0eha4HngME7STAa+jjNbthGYr6q5vgg40FiyNL2tuLjYWdSgpOTKTlXm7d/PdQkJDF+yhLDly229WRPQgi1ZenLO8kmcRPllVf3fVmX7gfUi8gXgf4AngPu9E6IxpqWxY8fy8MMPs3fvXjIzM6msrGRoRQVRJ09y4ORJiktLGZScTEpGBmI3mzbGKzwZWR4HSlQ1o5N624GxqjrCC/EFPBtZGn9qbGxk586d5Kxfz/DPPmPYqVMcGz2awxMnMmbMGJYuXcq4ceM6b8iYXhZsI0tPkuUFYI2qfrmTei8Bq1Q12gvxBTxLliYQNDQ0kJOTwyfr1nE2PJyLkZGXyyZOnMito0YxrF8/SE+HCJsEb/wv2JKlJ5+aM8CELtQb79Y1xvSSqKgoFi5cyKxZs9iyZQu5ubmXV/0pPHiQnS+9xOi4OMZdfz3xjzwCI0LiwI8xXuPJpSPZQIaIrGyvgojcCcwBPuppYMYYz/Xv35/ly5fzJ3/yJ0yfPh0RYfiJE/Srr6eiooKcrVt5Y/PmP1rswBjTMU8Ow96Ic/PnRuAl4EXgCKA4o8mv4iy0Hg4sUtWQSJh2GNYEsvLycrIyM6nKzGRcURHF48ZRMno0YWFhzJw5kwULFjAwLg7CwsAmA5leFGyHYT29zvIx4F9pe0QqOIn0+6r6tHfCC3yWLE0wKCsrI2vjRgoOH6Yp7MrHNzw8nFsiIpjWvz8xt98O48db0jS9ok8nSwARmQF8D1gIjMRJkqU4o86nVHWXt4MMZJYsTTApLi4mMzOTo0ePAhB58SI35OQQpcqIESMY8cQTxM2Y4ecoTSjo88nSXM2SpQk2qkphYSEffvghdfn5TN27l7CmJs7HxrJrzhxmZ2Rw4403EhcX5+9QTR9myTLEWLI0wUpVOXToEB+tX09kTg5nBg2iIjERgMjISGbPns28tDTiVWHsWD9Ha/oaS5YhxpKlCXaqSkFBAVlZWRw/fvyqsrSDB5khwuiFC4lZuRKSkvwUpelrgi1ZtnudpYg824N2VVW/1YPnG2N6iYgwadIkUlJSOHjwIFlZWZSVlRFTW0vi8eOUAsd/9zsuRkYy+557SLA1Z00IandkKSJNPWhXVTW8B88PGjayNH2NqnLw4EGyN2wgZts2hp88SVVCAp9Mm0ZYWBjTpk1jwYIFDB40yHmCzZ413RBsI8uOkuXDPWlYVZ/ryfODhSVL01epKocPH2bbO+9QVlJCbf/+l8tEhBuGD2debS0D7rgDUlIsaRqP9JlkabrGkqXp61SVoqIiNm3aRFFR0eX90/LzGVRVxZAhQxj+wAMMW7XKf0GaoBNsydJWVDbGdEhESE5OJjk5maNHj7J582ZK9+xhYHU1ABWVlbyTl0fS+fMsWLCA2NjxDBsmhIfEiRgTKrqVLEUkDpgNJALFqrrdq1EZYwLSuHHj+MpXvkJpaSk5kydTl5lJWFMT9TExFBUVcejQMfbuvZVJE8fwxanHmbJqGpGxkZ03bEyA82QhdUQk3p0lWwFkAr8DvtWifLWIFItIh/e8NMYEt1GjRvH5r32NZT//OTGrVhHmLqFXVjaSyso6Tn2wnk//7tesu/VHfPqbbD9Ha0zPdTlZikgskAV8AzgLvIez1F1LG4HRwN1eis8YE8ASExO5++67efzxx8nIyCAqCiIjLnBd9W4uXbrEmZIS3n/7XbZs2UJdXZ2/wzWm2zw5DPtnwAzgf4Fvqur51peXqOphESkAbvJijMaYAJeQkMCKFStYuLCGrZuy+fT5AdTvqSPy4iWqJ0FmZiZbtmxh1qxZzJkzB7kQzYBh/fwdtjFd5sktunYDg4EJqlrv7msCXlDVr7eotxFIU9XRPog34NhsWGP+2IULF8jNySXnzY+oGXD1iLLhXD+ufb+cxFlTSHvkVkbeMAYJs8tOQk1fng07AXi3OVF2oAIY2v2QjDHBLjo6mhsX3cjcG+eye/dusrOzOXXqFAARnzTQUH2a0g+2cuTjw5Sv+i533T2cKVOEMI9mURjTezxJlheB6C7UGw3UdC8cY0xfEh4ezvTp05k2bRoFBQVs3foRJ+sOXC7PjUim4P39fJRdzNSpI7j77hHMnh1BVJQfgzamDZ4ky4PADBGJVtULbVUQkQRgGhBS97Q0xnSsef3ZSZMmUbqslKw1H1Dyxh4KNRmA2tpacnMPk59fzJ0RB0i/MYlrv7aA+JHxfo7cGIcnBz1eA4YDP+mgzt8DccCangRljOm7Ro0axZe+9xW+9cZf8Nj3hQkTjhMRcRGAyLoq4vZt4eBv1rBh+V9xKK/Az9Ea4/Bkgk9/IA+YBGzFSZ6/AD7Eud5yFXAzsAdIb2/02dfYBB9jeqa+vp7t23fx+utHiMgp5bqTewA4OyieypVDSE5O5oYbbmDSpEmXr+c0wS/YJvh4tDasiIzBSZKzAcW5zrK5AQHygZWqeszLcQYsS5bGeEdjYyOfffoZHz3/AY1bDnNmSgJNk6+cKRo0aBBjmiYzPjaOKV9IJ3pAV6ZQmEDVp5Pl5SeJ3AGsAMYD4cAxYD3wmqr25NZeQceSpTHepaqUlJSwLXsb+w/sp6mpyd0vxK49T9K5cmIHD6L/l77AnIfSGTzYzwGbbgmJZGmusGRpjO9UV1eTm5vLxx9/zOkDStoH+wFQhNdG3UPU8LHMn5/IypVJTJgQZncJCyJ9JlmKyKvAc8AGtYzaLkuWxvheQ0MDO3N2sv3ZD4jcdYwKTeSDYVcWCuvXrx/XThnCl6aWM/WB2UTF2bUnga4vJcsmnPORZcCLOCv12NS0VixZGtN7VJXCw4Wsez2fD3OEysorx2CTawpZfHoLMQkJjFh1Ewv+8k4/Rmo605eS5S+B+3GWuGuutBV4HlijqrW9EmGAs2RpjH9UVVXx/vv5rFtXztGjg1he8i6JF8oBKJuSROLnJ5Kens7UqVOJiLBb9waaPpMsAUQkClgJfB1YijOZR4HzwCvAb1U1pO+/Y8nSGP+6ePEiOz/+jLxfvUtYXhHRdRc4+vmxSJxzAjM2NpYZM2YwqDKRKUtS6D+sv58jNtDHkuVVFUVGAg8CXwUmu7sVZ2Wf54H/VtUTvggykFmyNCYwqCrHio/x0RtbOHS2kMbGxstl9dVRTHijkPjYaAbMnU7ak19m3KQYmxDkR302WV71JJG5OKPNVcAAnKTZiHP5yPPA26ra2H4LfYclS2MCz/nz59m1axd5eXlUVVXRmBPGhP2FAFRFDmR98v2kTR3K5z43knnz4oi2SzZ7XUgky8tPFokB7gUeAha3KCpX1aQeRRYkLFkaE7iampo4dOgQWc9+yIX3dhNfXcP2wRnsG5B2uU5S0mBuywhj2U2DSZo+0o/RhpaQSpZXNSRyC/ASkAioqoZ7peEAZ8nSmOBQWVnJlleyeDu7gcOliTQ1XfmKWlC+mckXiomecA3TnryHSUvTOmjJeENIJUsRicOZMfsQMA9nyTuAYlW9pqfBBQNLlsYEl8bGRj755ABvvHGEvDxoPBfOfSVrCHfPHB1dPo5Ji1KZOXMmKSkpth6tjwRbsuzWfGoRWQJ8DbgHiMFJkheAt3DOWW70VoDGGONN4eHhzJyZxsyZaVRUnOb917Zz8uWhxBefojY+Fk0SDhw4wIEDB4iPj2dK6kwSy6KY9kC6LXYQwjyZDZuMMxv2QWAsV0aR+cBvgZdU9YwvggxkNrI0JvhdunSJ/Oxd7Nr8MaWXrp7Uf2FPFJNzDxIVF0fs4rnMePJ+xo6NsJm0PdSnRpYiEosz4/UhYAFOghTgDPAy8Jyq5vs4RmOM8amIiAhmL0xn9sJ0Tp8+za5du8jPz+fcuRoGHDoLQENNDbs2l/L0/hxSUwexYkUSS5YMISbGz8GbXtHRCj7P4STK/jgJsgnIxDnMulZVG3oryEBmI0tj+qbGxkb27TvItqfe42L2fmLP17J25F1URyVcrjNwYCx3XFPO4hXjmXhLKhJmw82uCraRZWdrwwIcAV7AWa2npJfiChqWLI3p+6rOVLH591t5c4dy5Ejs5Zm0ok3cV7KG2KZ6IhOHkvbTL3L94mk2KagL+lKy/G/geVX9sHdDCi6WLI0JHarK/v1FvPVWIdnZ9QwoO83Npz4A4EJUFGUPjGBAwkCmTZvG9OnTGTJkiJ8jDlx9JlmarrFkaUxoqqurI/sPORx8aQsRe4upTB5C3fyrT2AOvDCG8eWXuP7Bm7lm/gQ7TNuCJcsQY8nSGFNSVEJ+Xj57C/dSW3vlhkzh7zUypvQYIkLl5NmMeGgly5ePJCHBDtNasgwxliyNMc0aGxspKCggPz+fzz4pZPQrR+nXdAGAjcNv4XjMKKKiIpk+fQC33TaS+fMGExUdmqPNYEuWdpM3Y4zxkvDwcFJTU0lNTaXm9ho+Gr+Jwtd3cLHwDGX9RgDQ0HCRHTtOs2N7BfeczWTUdUlc+6UlTLolzQ7TBjAbWfaQjSyNMZ0pPlrGuvVHyMqq4sQJZ9GzwQ2VfO74WwA0hocT8ecLmTFnBqmpqURF9f2VgoJtZGnJsocsWRpjuqqpqYlPPini7beLOb9uB6llewE4PWoI526JByAyMpLU1FTGj5zCxOQxxCfF+zNkn7FkGWIsWRpjuqO+rp5tr27l0GvbKRt8gaaxV0/6adoCE4qOETVlAmPuX0rG/elERfWdw7SWLEOMJUtjTE9VV1eze/duPv30U06dOkVjYxhDfn+ahAvVAGwdeiMnktLIyIhn2bKRpKcPIdjXPbBkGWIsWRpjvEVVOXnyJNuy8jnxLxuJKa/kkkTwypj7uBh25TzmkCHR3D2qiIy7ppOyODiX2bNkGWIsWRpjfKGpqYk9Wz8j+/VPeadwGGfPXkmWMZdqua9kDYIiQ4cw6ef3cv3M6+nfv78fI/aMJcsQY8nSGONrFy9eIjv7CO++e5ydOy8y/uQBMs7kAlA9eABnPjeYsLAwkpOTufbaa0mdnEpMbGDfDsWSZYixZGmM6U21tfVk/z6boteykX1HKJuWxKXrIq+q07gzkolVtQy/eQbpD93EkHGD/RRt+yxZhhhLlsYYf6mqqGLf3n3sO7SP4uJiAFSh/9pahp09BcDHQ9JpmjObhQuHsmzZGBISAuMaTkuWIcaSpTEmEFRXV7Nnzx52bNlH/K8/IlyduyyuGX0v5yPiAAgPF9LS+rFiSj3zPz+DAcMH+C1eS5ZBRkS+CKwGrgfCgf3Ab4FnVLWpo+eCJUtjTOApLShhx4ubOLytmLcb5wBXZstGNF3kgWOvEB6mhKeMY/rfr2TKtDQiIyPbb9AHLFkGERH5JfAoUA9kAheBm4F4YC2wSlUbO2rDkqUxJpAdOVLB+vVH+eijakpLw0g+f4RF5ZsAqI2L5dS9w4iKiiIlJYW0tDRSUlJ6Zbk9S5ZBQkQ+D7wKnAAWqmqBu3848CEwBfieqj7VUTuWLI0xwaKgoJycFzdTvzGbyFOnKUtN4sKcflfVaTrYjwmnLzBy6QwyvrqYgT46VGvJMkiISB4wC3hQVf+rVdkiIAsnkY7q6HCsJUtjTLBRVQ7lFnC4+DCHThymoqLiclnkOw2MKj8OwN5B11I7fz4LFgxh6dIxDBnivctRLFkGAREZDRwDGoAEVa1ro04JMAqYr6rZ7bVlydIYE8xUlfLycvbu3cuuvAP0/9XHRDU1APD2iNupiE4EICwMUlIiWZZax40rr2PY+OE9et1gS5ahej/LGe52T1uJ0pWLkyxnAO0mS2OMCWYiwrBhwxg2bBiLFy/m0MIC8l/ewoncI1REDL1cr6kJDu6/wPT3X+Wdf/sfdHQSE//fcqbfOIMBA/w3q7a3hGqyTHa3RzuoU9yqrjHG9HkTM1KYmJECwIrCSt599xjbtlVTXNzEsAuniGmsB+DiiQo+yPkA7Q+LFi3yZ8i9IlSTZZy7Pd9BnRp3+0c3kxORbwLfBBg7dqx3IzPGmAAxfvxgVq8ezOrVcOxYNdn/u53zfxiBlJygeuRACBOmTJni7zB7Ragmy+aLjrp1wlZVnwWeBeecpbeCMsaYQDVmzEDu/8Gt8INbOXn4BAf3HaSy8QyJiYn+Dq1XhGqyPOdu4zqo01x2roM6xhgTcoZPSGL4hCR/h9Grgvz2od1W5G7HdVBnTKu6xhhjQlSoJstd7naqiLR34VB6q7rGGGNCVEgmS1U9BuwEooBVrcvdRQlG4yxKsK13ozPGGBNoQjJZuv7B3f5MRCY27xSRYcCv3B9/2pXF1I0xxvRtoTrBB1V9VUSewbnjyG4ReZ8rC6kPAN4AnvZjiMYYYwJEyCZLAFV9VES2At8BFnHlFl3P08VbdBljjOn7QjpZAqjqy8DL/o7DGGNM4ArJhdS9SUTK6XjZvI4MBSo6rWVasj7zjPWXZ6y/PNOT/hqnqkGzooElSz8SkbxgWnU/EFifecb6yzPWX54Jpf4K5dmwxhhjTJdYsjTGGGM6YcnSv571dwBByPrMM9ZfnrH+8kzI9JedszTGGGM6YSNLY4wxphOWLI0xxphOWLL0EhH5oohsEZFqEakRkTwR+Y6IdLmPRSRSRG4WkX8RkRwRKRORBhEpFZFXRWSxD99Cr/NGn3XQ9k9ERN3HE96I19+83V8iEiMiPxCRXBGpEpFaETkiImtEZL634+9t3uwvERktIv8uIgdEpE5E6kWkQER+LSLjfRF/bxGRySLyXRF5SUT2i0iT+7m5t4ft+uzz7Reqao8ePoBfAgrUAW8Da4Gz7r7XgfAutrPUfY4CZW5brwC7W+z/W3+/30Dqs3baTgcuAU1ue0/4+/0GWn8ByUCB+/yTwJvA74EdQAPw//39ngOlv4AZwBn3ucdw1o1+Ayhx950D5vn7Pfegr37R4vul5ePeQOj/QHn4PYBgfwCfb5HcUlrsHw7sdcu+28W2bgJeBRa0UXa/mwAUWOLv9x0ofdZG29HAHqDU/YAGfbL0dn8B/YFDzX98AZGtyocAk/z9vgOov7Ld5zzbsq+ASOA5t+wTf78BYbd7AAAMEUlEQVTvHvTXN4B/BO4DJgBZPUmWvvx8+7Wf/B1AsD+APPc//6ttlC1q8UsT5oXX+o3b3nP+ft+B2mfAz9zn3wm80EeSpVf7C+f2dAq86O/3Fuj9BfTjykgrqY3ykS3KY/393r3Ufz1Nlr32ndibj+A8dhwgRGQ0MAvnsNWa1uWquglnhJMEzPHCS+5yt6O90JZf+LLPROQG4M+Al1X1Dz2P1v+83V8iEgU84v74U+9FGhh88PvViHNEB0DaKG++9u48ziHHkOaH78ReY8myZ2a42z2q2t4HJbdV3Z5IcbdlXmjLX3zSZyLSD3gRqAS+2/3wAo63+2sWzmHWY6q6T0TmuZOh/kNEfiQic3sasJ95tb9U9SKQ6f74IxGJbC5z//337o/PqTt0CnG9/Z3Ya0L+Fl09lOxuO7rrSHGrut0iIknAQ+6Pr/WkLT/zVZ/9GJgMPKCqfemuEd7ur+vcbYGIvAA82Kr8b0TkNeArHXzZBTJf/H49CmzAGZHfJiJ57v50YBDwFPCkh3H2Vb32ndjbLFn2TJy7Pd9BnRp3G9/dFxGRCOAlYCCQGeSHGL3eZyIyD/ge8IaqvtKD2AKRt/trsLtdiHOz838Gfg2cdvf9CmeCxlng654GGwC8/vulqoXu79h/Abdx9WmQPGCzOwI1vfSd6A92GLZnms9h+Prwy6+Bm3GmrX/Zx6/la17tMxGJAX6L8+X+qDfaDDDe/h1r/sxH4Bw6fFJVD6tqlaq+BdzlvtaDQXr9oNc/k26i/AyYCKzEuYdjIk5fDQJeE5G/8dbrBbne+k7sdZYse+acu43roE5z2bkO6rRLRJ4CHgZOADer6onutBNAvN1nPwEmAX+qqsF8Lrc93u6vlnX+s3WhquYBH+N8NyzuQnuBxqv9JSIJONdUxgPLVfUtVT2tqhWq+iawHGdiz1+LSEpHbYUIn38n+osly54pcrfjOqgzplXdLhORfwEeB8pxEmWBp20EoCJ3660+uxtn8YEHRSSr5QPniwxgtbvvN92I19+K3K23+qtlnSPt1Gnen9SF9gJNkbv1Vn/djjOKzFHVwtaFqnoI2I4zUl/c1SD7sCJ365PvRH+yc5Y903wpx1QRiWlnQkR6q7pdIiL/CPwpzrmkW1R1b/fDDCi+6LMwnOu32jPefSR0sb1A4u3+2tni30Nw/hBrbai7rWmjLNB5u7/GutvqDupUudvBHdQJFT77TvQ3G1n2gKoew/nyiQJWtS4XkUU4kwFOANu62q6I/BRndt0ZnET5iVcCDgDe7jNVvUZVpa0HzqUkAE+6+6Z77530Dh/0VynOSAic8+Ct2xsEzHR/zGtdHuh88Jk87m5ntbxspEV7kTiX40D7I/WQ4avvxIDg71URgv0B3MuVFSkmttg/DGfZtT9a2glnBZX9wD+00d7fuc85A8zy9/sLhj7r4HVeoG+s4OPt37E7ubIm7PQW+/sBv3PL8nDvdxtsD2/2l/uc8+5zngaiW5RFA8+4ZZXAQH+/dy/1XxadrODTye+Xx/0fDA+/B9AXHjjT7RXnRP8fcBYKrnb3raXVosEtvsRfaLX/c1xZOivXrdfW4y/8/Z4Dpc86eY0+kSx90V/AP7nlF4DNbhul7r4SWqzpGYwPb/YXzrWozesylwJvuW0ed/fVA3f5+z33oK9mAjktHs0Lnh9sud/D3y+P+j8YHnbO0gtU9VER2Qp8B+fcWTjOX13PA8+oalMXm2p5zmO2+2jLJoJ8qTIv9llI8HZ/qeqTIpINPIazkkoszsXi/wr8VFXbOpcZNLzZX6r6oojsxrmWdwFwq1tUirOQ+r9qcM8pGADc0Mb+bs/u7Yufb3H/CjDGGGNMO2yCjzHGGNMJS5bGGGNMJyxZGmOMMZ2wZGmMMcZ0wpKlMcYY0wlLlsYYY0wnLFkaY4wxnbBkaYwXiMgYEfkfETkuIpdEREXkF23U+75b9tVeiClMRL4jInkiUiMi1SKyRUS+0I22Frtxt3zc26rOD939P/Tam2g/nl+0jsfXr2lCm63gY0wPiYgAr+HcTWEv8CFwEdjRRvV7cJZOe9vHMYXjLDH2OZzlyzbirGV6M/CyiMxV1ce70fRJYIP77yIvhNpdO7iyUP6DfozDhAhLlsb03DU4ibIYmKaql9qqJCLDgXnAB6pa6eOYvoeTKPcCN6nqSTeGFGAL8JiIZKpzA2NP7FfVh7waaTeo6svAywAiYsnS+JwdhjWm55pvZnukvUTpugvnM7fWl8G4o8ofuD+ubk6UAOrcQPzP3R//ypdxGNOXWLI0IaXl+S0Recg9n3deRE6IyHMikuiW9RORH4nIQRGpF5FiEflxy3saisg1blub3F2LOjmHdg/OXRfeaNHG5fN8IjJaRF4QkTIRqRWRnS3PC4rIfBFZJyKn3fIPRSS9jdeZi3M7pBJV3dxG+Rqcw8TpIjLKow7sJhFZJCKVItLQeiQoIsNE5FciUuL29SG3r2NEJMvtn8W9Eacx7bHDsCYkicjPcA5VbsI5BzcP+DowW0TmA+8CU9zyQzh3TvhLIBH4pttMDc55syRgGVefz2v9egOBJTi3OjreRpVxwMdum5twbpA7H/i9iHwR51ZarwD5wHvANGAx8KGIzFTVgy3amuFuc9uKRVVrRWQPMN19lLZVz1tE5AGcWzo1ALer6nstykYCH+Ecyj6FczunaOBxnPdnTECwZGlC1YM4Nz7eByAig3Du3H69u60CklW12i2fjpN8viEiP1bVo6paATzkjnqW0fH5vDuBSNo/BPsQ8BTwZ6ra6L7mapz7Av4T0B/4kqquccvCcM7Z3Y9zWPXhFm0lu9ujHbz/YpxEmdxBnR4TkSeBn+HcCHiFqn7SqsqvcBLlemCVqp53n5cEZAJpvozPmK6yw7AmVP1Nc6IEUNUzwK/dH9OAbzYnSrc8H1gHCM4o01P3uNvX2yk/CvygOVG6ngVO44wyNzQnSjeeJpwkBM6ItaU4d3u+g3hq3G18J3F3i4iEi8gvgX/EmWQ0p3WiFJFxOJOQLgGPNidKAFU9ATzhi9iM6Q5LliZUtXW49JC7PdoykbZQ4G5HevJCIhKLM/LcraqH26n2gao2tNzhJs6iDuJtLx5pbsKTOL0oFuePgkdxLqOZr6rH2qi3ECfWbapa1LpQVdcDZ3wYpzFdZodhTagqaWNfTQdlLcv7efhay7mSQDyJp8OYVLXGucST6FZF59xtHO1rLjvXQZ3u+j7Od0s+sLz1HwEtNE8u6uxw8SAvxmZMt9jI0oQk9zBmezoq647ODsF25TU9ianI3Y7roE7z5S5FHdTprndwDh9PB/60C/U7GgF7+//CmG6xZGmMD7mXmtwOFKrqp730sjvdbVuXlTQfFr7W/XGXD14/H+c86kngH0Tkb9up1zwruKOk3lGZMb3GkqUxvnUTkEDHo0pv24ZzGcZoEVnYRvkqnJm5uarqk8tGVHU3zjnJEuCvReSf2qi2BWdUOc+d7HMVEVkGDPZFfMZ4ypKlMb7VfAjWp6v2tORODGpOTs+IyLDmMne5u5+6P/7Yx3EcxEmYR4AnRORpdx3d5vIjOIdsI4BfuiPe5jiHA//sy/iM8YQlS2N8xL0WciXONYbbevnlf45zgX8aUCAir4vIH4BPcRZR+PdurAvrMTchLgAOAN8BfuP2S7PVOJN4bgcKReT3IvImzkzfGiDHrdfeJCFjeoUlS2N8Zz4wHHhDVXv1Mg53dHkX8BjOJTHLcK4P/RhncYPu3HGku7GUuq+9G2eVpJdEJMItKwEygP/AmcyzErgOeAbnDimJbjMVvRWvMW2RXv4MGxMyROTnOEvq3dpyibdg5K5S9CGwSVUX99JrXoOT6M8Dg9qbwdy8Dq+qSlvlxniDXWdpjO/sA36Ik2T6ilQRecH999OqmteTxtxzmLNatyMiY4D/BsKB/2qdKN31cm/tyWsb4wkbWRpjOtViZNnSKlV9tYftRuDcAaUY2I+zYs8YYCbO4g+fATe2XHrQfd4vgO+23GcjS+NLliyNMX7jjiz/DlgKjMe5zOYCcBDncpunVLWm/RaM6R2WLI0xxphO2GxYY4wxphOWLI0xxphOWLI0xhhjOmHJ0hhjjOmEJUtjjDGmE/8Hvyjlmiw92DwAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"m0 = 0.25\n",
"dmdt=0.05\n",
"u = 250\n",
"t_T = np.linspace(0,4,100) #time interval\n",
"dt = t_T[2]-t_T[1]\n",
"m_T=np.zeros(len(t_T))\n",
"m_T[0]=0.25\n",
"for i in range(1,len(m_T)):\n",
" m_T[i]=m_T[i-1]-dmdt*dt\n",
"v_T = -u*np.log(m_T/m0)\n",
"\n",
"#plotting to demonstrate convergence\n",
"plt.plot(current_mass_euler, current_velocity_euler,color='k', linewidth = 3, linestyle='-', alpha=0.5, label=\"Euler Solution\")\n",
"plt.plot(current_mass_heun, current_velocity_heun, linewidth = 3, color='b', linestyle='--', alpha=0.5, label=\"Heun Solution\")\n",
"plt.plot(m_T/m0, v_T,linestyle= 'dotted', linewidth = 3, alpha = 0.5,color='r', label=\"Tsiokovsky Solution\")\n",
"plt.title(\"Simplerocket Mass vs. Velocity\")\n",
"plt.xlabel(\"mf/m0 [kg]\")\n",
"plt.ylabel(\"Velocity [m/s]\")\n",
"plt.legend(loc=0,fontsize=15);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Alternate**"
]
},
{
"cell_type": "code",
"execution_count": 33,
"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": [
"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])\n",
"\n",
"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": "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": 2
}