Skip to content

Instantly share code, notes, and snippets.

@cjtu

cjtu/AST501_hw1.ipynb

Created February 5, 2020 18:14
Show Gist options
  • Save cjtu/682d4938456908ca9e9307423ca87fc2 to your computer and use it in GitHub Desktop.
Save cjtu/682d4938456908ca9e9307423ca87fc2 to your computer and use it in GitHub Desktop.
Download ZIP
AST501 Homework 1
Display the source blob
Display the rendered blob
Raw
AST501_hw1.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Homework 1: Shapes and Orbits\n",
"\n",
"Christian Tai Udovicic (cjt347)\n",
"\n",
"\n",
"Code, plots, derivations and solutions to the problems are presented below."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# Imports\n",
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib import cm\n",
"\n",
"# Constants\n",
"G = 6.674e-11 # [m^3/(kg s^2)]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Q1) Maclaurin vs Jacobi\n",
"\n",
"Input parameters of each body"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>name</th>\n",
" <th>a</th>\n",
" <th>b</th>\n",
" <th>c</th>\n",
" <th>rho</th>\n",
" <th>omega</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>Mercury</td>\n",
" <td>2439700.0</td>\n",
" <td>2439700.0</td>\n",
" <td>2439690.0</td>\n",
" <td>5427.0</td>\n",
" <td>1.240010e-06</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>Venus</td>\n",
" <td>6051800.0</td>\n",
" <td>6051800.0</td>\n",
" <td>6051790.0</td>\n",
" <td>5423.0</td>\n",
" <td>2.992400e-07</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>Earth</td>\n",
" <td>6378100.0</td>\n",
" <td>6378100.0</td>\n",
" <td>6356790.0</td>\n",
" <td>5514.0</td>\n",
" <td>7.292120e-07</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>Mars</td>\n",
" <td>3396200.0</td>\n",
" <td>3396200.0</td>\n",
" <td>3376190.0</td>\n",
" <td>3933.5</td>\n",
" <td>7.088218e-05</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>Jupiter</td>\n",
" <td>71492000.0</td>\n",
" <td>71492000.0</td>\n",
" <td>66853900.0</td>\n",
" <td>1326.0</td>\n",
" <td>1.773408e-04</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
" <td>Saturn</td>\n",
" <td>60268000.0</td>\n",
" <td>60268000.0</td>\n",
" <td>54363900.0</td>\n",
" <td>687.0</td>\n",
" <td>1.636246e-04</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
" <td>Uranus</td>\n",
" <td>25559000.0</td>\n",
" <td>25559000.0</td>\n",
" <td>24972900.0</td>\n",
" <td>1270.0</td>\n",
" <td>1.041366e-04</td>\n",
" </tr>\n",
" <tr>\n",
" <th>7</th>\n",
" <td>Neptune</td>\n",
" <td>24764000.0</td>\n",
" <td>24764000.0</td>\n",
" <td>24340900.0</td>\n",
" <td>1638.0</td>\n",
" <td>1.083383e-04</td>\n",
" </tr>\n",
" <tr>\n",
" <th>8</th>\n",
" <td>Pluto</td>\n",
" <td>1188300.0</td>\n",
" <td>1188300.0</td>\n",
" <td>1188290.0</td>\n",
" <td>1854.0</td>\n",
" <td>1.139000e-05</td>\n",
" </tr>\n",
" <tr>\n",
" <th>9</th>\n",
" <td>Ceres</td>\n",
" <td>487300.0</td>\n",
" <td>487300.0</td>\n",
" <td>454690.0</td>\n",
" <td>2162.0</td>\n",
" <td>1.923000e-04</td>\n",
" </tr>\n",
" <tr>\n",
" <th>10</th>\n",
" <td>Eris</td>\n",
" <td>1163000.0</td>\n",
" <td>1163000.0</td>\n",
" <td>1162900.0</td>\n",
" <td>2520.0</td>\n",
" <td>6.734000e-05</td>\n",
" </tr>\n",
" <tr>\n",
" <th>11</th>\n",
" <td>Makemake</td>\n",
" <td>717000.0</td>\n",
" <td>717000.0</td>\n",
" <td>710900.0</td>\n",
" <td>1700.0</td>\n",
" <td>7.646000e-05</td>\n",
" </tr>\n",
" <tr>\n",
" <th>12</th>\n",
" <td>Haumea</td>\n",
" <td>980000.0</td>\n",
" <td>759000.0</td>\n",
" <td>498000.0</td>\n",
" <td>2600.0</td>\n",
" <td>4.457600e-04</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" name a b c rho omega\n",
"0 Mercury 2439700.0 2439700.0 2439690.0 5427.0 1.240010e-06\n",
"1 Venus 6051800.0 6051800.0 6051790.0 5423.0 2.992400e-07\n",
"2 Earth 6378100.0 6378100.0 6356790.0 5514.0 7.292120e-07\n",
"3 Mars 3396200.0 3396200.0 3376190.0 3933.5 7.088218e-05\n",
"4 Jupiter 71492000.0 71492000.0 66853900.0 1326.0 1.773408e-04\n",
"5 Saturn 60268000.0 60268000.0 54363900.0 687.0 1.636246e-04\n",
"6 Uranus 25559000.0 25559000.0 24972900.0 1270.0 1.041366e-04\n",
"7 Neptune 24764000.0 24764000.0 24340900.0 1638.0 1.083383e-04\n",
"8 Pluto 1188300.0 1188300.0 1188290.0 1854.0 1.139000e-05\n",
"9 Ceres 487300.0 487300.0 454690.0 2162.0 1.923000e-04\n",
"10 Eris 1163000.0 1163000.0 1162900.0 2520.0 6.734000e-05\n",
"11 Makemake 717000.0 717000.0 710900.0 1700.0 7.646000e-05\n",
"12 Haumea 980000.0 759000.0 498000.0 2600.0 4.457600e-04"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"d = {\n",
" 'name': ['Mercury', 'Venus', 'Earth', 'Mars', 'Jupiter', 'Saturn', 'Uranus', \n",
" 'Neptune', 'Pluto', 'Ceres', 'Eris', 'Makemake', 'Haumea'],\n",
" 'a': [2439.7e3, 6051.8e3, 6378.1e3, 3396.2e3, 71492e3, 60268e3, 25559e3, \n",
" 24764e3, 1188.3e3, 487.3e3, 1163e3, 717e3, 980e3],\n",
" 'b': [2439.7e3, 6051.8e3, 6378.1e3, 3396.2e3, 71492e3, 60268e3, 25559e3, \n",
" 24764e3, 1188.3e3, 487.3e3, 1163e3, 717e3, 759e3],\n",
" 'c': [2439.69e3, 6051.79e3, 6356.79e3, 3376.19e3, 66853.9e3, 54363.9e3, \n",
" 24972.9e3, 24340.9e3, 1188.29e3, 454.69e3, 1162.9e3, 710.9e3, 498e3],\n",
" 'rho': [5427, 5423, 5514, 3933.5, 1326, 687, 1270, 1638, 1854, 2162, 2520, \n",
" 1700, 2600],\n",
" 'omega': [1.24001e-6, 2.9924e-7, 7.29212e-7, 7.088218e-5, 1.773408e-4, \n",
" 1.636246e-4, 1.041366e-4, 1.083383e-4, 1.139e-5, 1.923e-4, \n",
" 6.734e-5, 7.646e-5, 4.4576e-4]\n",
"}\n",
"\n",
"df = pd.DataFrame(d)\n",
"df.head(13)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Maclaurin spheroids"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x7f2951a77ed0>"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAFNCAYAAABIc7ibAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd1xX1f/A8deRLcOFe+EWEUQFFLUcufe3NDUzR2mmllpWZkPzV33N+mb5rb5li4arNMu0Ya40UxmKiHvhwoUKggyBz/n9wQiUaR+4H/i8n48Hj7j3nnvO+35SeXPOuecorTVCCCGEEKJ0VTA6ACGEEEIIayRJmBBCCCGEASQJE0IIIYQwgCRhQgghhBAGkCRMCCGEEMIAkoQJIYQQQhhAkjAhhCGUUkFKqddKqO4GSqkEpZRNSdRfQLseSimtlLItpfYOKKW65XOtm1LqXGnEIYS4O5KECSHypJSKUkrdUkq533Y+PDPR8DAmssJprc9orV201ul3c79Sao5S6lRmIndOKbXS3DGag9baS2u91eg4hBB3R5IwIURBTgGjsg6UUt6Ak3HhFO6f9kIppcYCY4CeWmsXwA/YZI7YihlHqfSmCSGMI0mYEKIgXwOP5DgeC3yVs4BSaoBSaq9S6oZS6qxSat5t17sopf5SSsVmXh93eyNKqSpKqXVKqStKqeuZ39fLcT1KKdUzx/E8pdQ3md9nDQE+qpQ6A2y+fVhQKbVVKfV/SqkdSql4pdSG23v4cvAHftNanwDQWl/UWi/J0fZWpdS/lVLBSqk4pdSPSqmqt9UxWil1RikVo5R6Mce9FZRSs5VSJ5RSV5VS32bdm9dzZJ4fnDnsGJvZtmden4tSyilziPe6Uupg5nPk/IyfV0qdz3z+I0qp+/J5fiFEKZEkTAhRkF2Am1LKM3N+1Qjgm9vK3CQjUasMDACeUEoNhYy5WcAvwH+B6oAvEJ5HOxWAL4CGQAMgCXi/mLF2BTyBPvlcfwgYD9QA7IFZ+ZTbBTyilHpWKeWXz7yyR4AJQB0gDVh82/UuQAvgPuCVHInTU8DQzFjrANeBD/J7DqVUc2A5MIOMz+9n4CellH0eMc0FmmR+9SEjYQZAKdUCmAb4a61dM69H5fP8QohSIkmYEKIwWb1hvYDDwPmcF7XWW7XW+7XWJq11BBlJQ9fMy6OBjVrr5VrrVK31Va31HUlY5vnVWutErXU88HqOOopqntb6ptY6KZ/rX2itj2Ze/5aMhPAOWutvgCfJSFT+AC4rpWbfVuxrrXWk1vom8DLw4G3J2qta6ySt9T5gH9Am8/zjwIta63Na6xRgHjDstqHHnM8xAlivtf5da50KvE3GcHCnPEJ/EHhda31Na32W3IlhOuAAtFJK2Wmto7J6+oQQxpE5B0KIwnwNbAMacdtQJIBSqgOwAGhNRg+TA/Bd5uX6QKE/7JVSFYFFQF+gSuZpV6WUTTEm158t5PrFHN8nAi75FdRaLwWWKqXsyOi5WqqU2qu1/i2Ptk4DdkDO4c382moIrFFKmXJcTwdq5vMcdTLrz4rLpJQ6C9TNI+w6ecSVdd9xpdQMMpI+L6XUb8DTWuvoPOoRQpQS6QkTQhRIa32ajAn6/YHv8yiyDFgL1NdaVwI+AlTmtbNkDI8V5hkyhu86aK3dgHszz2fVcxOomKN8rbxCLUI7xZLZe/cdEEFGkpmlfo7vGwCpQEwRqjwL9NNaV87x5ai1ztm7mPM5oslI3ABQSqnMtnP1Rma6kEdcOZ9lmda6S2Z9GnizCPEKIUqQJGFCiKJ4FOiROfx2O1fgmtY6WSkVQMbcqyxLgZ5KqQeVUrZKqWpKqbyGAV3JmAcWmzlRfe5t18OBkUopO6WUHzDsHz9RPpRS4zJfNnDNnEjfD/ACduco9rBSqlVmD958YFURe+w+Al5XSjXMbKu6UmpIAeW/BQYope7L7JV7BkgB/sqn7AuZLznUI2NINeuZWiileiilHIBkMj7ru1q+QwhhPpKECSEKpbU+obUOzefyFGC+UioeeIWMZCDrvjNk9KA9A1wjI5lqk0cd75Ix1ymGjInxv952/WUyetSuA6+S0ftWUm4Ac4AzQCywEHhCa/1njjJfA0FkDDs6kjHhvijeI6PXcEPm57UL6JBfYa31EeBhMl5siAEGAYO01rfyKP4qGUOQp4ANmTFmcSBjyDgmM+Yamc8ohDCQ0trsPfhCCFFuKaW2At9orT81OhYhRNkmPWFCCCGEEAaQJEwIIYQQwgAyHCmEEEIIYQDpCRNCCCGEMIAkYUIIIYQQBihzK+a7u7trDw8Po8MQQgghhChUWFhYjNa6el7XylwS5uHhQWhofssVCSGEEEJYDqXU6fyuyXCkEEIIIYQBJAkTQgghhDCAJGFCCCGEEAYoc3PC8pKamsq5c+dITk42OpRyy9HRkXr16mFnZ2d0KEIIIUS5UC6SsHPnzuHq6oqHhwdKKaPDKXe01ly9epVz587RqFEjo8MRQgghyoVyMRyZnJxMtWrVJAErIUopqlWrJj2NQgghhBmViyQMkASshMnnK4QQQphXuUnC7sai34+arS6lFGPGjMk+TktLo3r16gwcONBsbQghhBCi/LDqJOy9TcfMVpezszORkZEkJSUB8Pvvv1O3bt1i1ZGWlvaPYvin9wshhBCi9Fh1EmZu/fr1Y/369QAsX76cUaNGZV+7efMmEyZMwN/fn7Zt2/Ljjz8CEBQUxPDhwxk0aBC9e/cGYOHChXh7e9OmTRtmz54NQLdu3bJ3CoiJiSFr66bb7x8zZkx23QCjR49m7dq1Jf7sQgghhCiecvF2pKUYOXIk8+fPZ+DAgURERDBhwgS2b98OwOuvv06PHj34/PPPiY2NJSAggJ49ewKwc+dOIiIiqFq1Kr/88gs//PADu3fvpmLFily7dq3QdnPe/8cff7Bo0SKGDBlCXFwcf/31F19++WWJPrcQQgiRn7R0E5fiU7gYl0xqugmAas721KzkiJujdS97VO6SsFd/OsDB6BtFLj/i452FlmlVx425g7wKLefj40NUVBTLly+nf//+ua5t2LCBtWvX8vbbbwMZb3SeOXMGgF69elG1alUANm7cyPjx46lYsSJA9vmC5Ly/a9euTJ06lcuXL/P999/zwAMPYGtb7v43CyGEsFCX45PZevgKoaevsedMLCevJGDSeZf1qFYRf4+q9Pepzb3NqmNTwbpeArOan87nridyPvbOJRZ2n8rd01S3siP1qlS863YGDx7MrFmz2Lp1K1evXs0+r7Vm9erVtGjRInf7u3fj7Oycq1xebyLa2tpiMmX8BnH7UhE57wcYM2YMS5cuZcWKFXz++ed3/SxCCCFEUcQlprJm7zl+irjAnjPX0RoqV7SjXYMq9PWqRb0qTtSs5IiDbQW0hqs3b3HueiJ7z8Sy4eAlvgs7Ry03Rybd25iHOjTA0c7G6EcqFeUuCStKj1UWj9nriVowwKztT5gwgUqVKuHt7c3WrVuzz/fp04f//ve//Pe//0Upxd69e2nbtu0d9/fu3Zv58+fz0EMPZQ9HVq1aFQ8PD8LCwggICGDVqlUFxjBu3DgCAgKoVasWXl5F/zyEEEKI4og4F0vQjijW779ASpqJVrXdmNmzOb1a1aRlLdciLW+UkpbOlsOXCforivnrDvLVzigWDmtDQKPCR4LKunKXhBmtXr16TJ8+/Y7zL7/8MjNmzMDHxwetNR4eHqxbt+6Ocn379iU8PBw/Pz/s7e3p378/b7zxBrNmzeLBBx/k66+/pkePHgXGULNmTTw9PRk6dKjZnksIIYTIsvvkVd7fcpztx2JwcbBluF89Rvo3oHXdSsWuy8HWhr6ta9O3dW22Hb3Ciz/sZ8SSnczq3YIp3ZqU63Uqldb5DNRaKD8/P531lmCWQ4cO4enpWey6SqInzBIkJibi7e3Nnj17qFSp+H8h8nO3n7MQQojyYc+Z6yz4+TDBUddwd7Hn0S6NebhjA1zNOMH+Zkoac9bs58fwaB5oV4+Fw3zK9FwxpVSY1tovr2vSE1bObNy4kQkTJvD000+bNQETQghhvc5eS+TNXw+zLuIC1V0dmDeoFSMDSmbulrODLe+O8KWxuwuLNh7FpDVvD29TphOx/Fh1Ejb9vmZGh2B2PXv2zH7rUgghhPgnklPTeX/zcZZsO4lNBcX0+5ox6d7GODuUbPqglGJ6z2bYVIC3Nxylppsjs/u1LNE2jWDVSdjMXs2NDkEIIYSwSNuOXuGlHyI5cy2R+9vW5bm+LalVybFUY5jWoxkXbyTz0R8naFHLhX+1rVeq7Zc0q07ChBBCiPLu8scRANR43KdI5WMSUpj/00HW7oumsbszyyZ2oFMT95IMsUBzB3lx9FICL62JpH2DqjSodvfLSFka2bZICCGEEAD8GnmRPou28WvkRWb0bMYvM+4xNAEDsLOpwKIRvlSooHj623DS81v5tQySJEwIIYSwcnFJqTy9MpzJ34RRu7Ij657qwoyezXGwtYxFU+tWduLVwV6Enr7OipDyM+/ZqpOwQ2cTzFJPt27d+O2333Kde/fdd5kyZYpZ6hdCCCFKyo7jMfR9dxs/7ovmqfuasWZKZ5rXdDU6rDv8q21dAhpV5e3fjhCXmGp0OGZh1UnY4fM3zVLPqFGjWLFiRa5zK1asYNSoUWapXwghhDC3tHQTb/12mIc/242TvQ2rn+jE072aY2djmamBUoq5g1oRl5TKe5uOGR2OWVjmJ13GDBs2jHXr1pGSkgJAVFQU0dHRdOnShbfeegt/f398fHyYO3du9nVPT08mTpyIl5cXvXv3JikpCcjoVctajDYmJgYPDw8ADhw4QEBAAL6+vvj4+HDsWPn4AyiEEKL0RccmMXLJLj7YcoIRfvVZ/+Q9+NavbHRYhfKqU4n729Vj6e7TXL5x537QZY0kYWZQrVo1AgIC+PXXX4GMXrARI0bw+++/c+zYMYKDgwkPDycsLIxt27YBcOzYMaZOncqBAweoXLkyq1evLrCNjz76iOnTpxMeHk5oaCj16pWv13SFEEKUjo0HL9F/8XYOXbjBeyN9WfCAD072ljH3qyimdW9KarqJJdtOGh3KP1bulqiIiIon7mbRx4q3H7hWaJlKznb4eBQ8Pp41JDlkyBBWrFjB559/zrJly9iwYUP2Rt0JCQkcO3aMBg0a0KhRI3x9fQFo3749UVFRBdYfGBjI66+/zrlz57j//vtp1qz8LTQrhBCi5NxKM/Hmr4f57M9TeNVx4/2H2tHI3dnosIrNw92Zob51+Wb3aaZ0b0pVZ3ujQ7prVtMTdjMlnZj41FxfwB3nbqak31X9Q4cOZdOmTezZs4ekpCTatWuH1poXXniB8PBwwsPDOX78OI8++igADg4O2ffa2NiQlpYGgK2tLSaTCYDk5L+7Wh966CHWrl2Lk5MTffr0YfPmzXcVpxBCCOsTk5bOQ5/s4rM/TzGukwffT+lUJhOwLJO7NSE51cTKkLNGh/KPlLuesMJ6rHJas+sS/+pY0yzturi40K1bNyZMmJA9Ib9Pnz68/PLLjB49GhcXF86fP4+dXcGbnHp4eBAWFkZAQACrVq3KPn/y5EkaN27MU089xcmTJ4mIiKBHjx5miV0IIUTZlLUQa0FCT13lZZJIVPB69Wr0upBO3GcH8ixb1AVdjda8pisdG1dl6e7TTLq3cZndV9JqesJKw6hRo9i3bx8jR44EoHfv3jz00EMEBgbi7e3NsGHDiI+PL7COWbNm8b///Y9OnToRExOTfX7lypW0bt0aX19fDh8+zCOPPFKizyKEEKJs01qz6kY8T5KIE4rP69Sgl0v5WW1+TEcPzl1P4o+jl40O5a4prcvWyrN+fn466+3BLIcOHcLT07PYdZmzJ8wa3O3nLIQQonQlp6bz0g+RrAo7R2cnR+ZXr0aTKb5Gh2VWqekmOi3YTJt6lfl0rJ/R4eRLKRWmtc4zwHI3HCmEEEJYs/OxSUz+Ooz95+OYfl8zRp1IooIqm8N1BbGzqcCQNnX4cmcUcYmpVKpY8HQfS2TVw5Et65bdSYlCCCHE7UKjrjH4v38SFXOTTx/xY2av5uUyAcsy2LcOqemaXyIvGB3KXSnRJEwp1VcpdUQpdVwpNbuAcsOUUlopVar9iZ71XUqzOSGEEKLErAo7x0Of7MbNyY4fpnWmZ6vyP93Gu24lGrk7s3ZftNGh3JUSS8KUUjbAB0A/oBUwSinVKo9yrsBTwO6SikUIIYQor0wmzYJfDjPru334eVRhzZRONKluHZ0MSikGt6nDzpNXuRxf9lbQL8mesADguNb6pNb6FrACGJJHuf8DFgJl79MTQgghDHQzJY3Hvwnjoz9OMLpDA76cEEDlimV38dK70cerFlrD1sNXjA6l2EoyCasL5FxF7VzmuWxKqbZAfa31uoIqUkpNUkqFKqVCr1wpex+yEEIIYW7nrifywP/+YtOhS7w62IvXhra22M23S5JnbVdqV3Jk8+Gyt1RFSf7fymsmYPZ6GEqpCsAi4JnCKtJaL9Fa+2mt/apXr26+CLf822xV2djY4Ovrm/21YMGCYt3/ww8/cPDgwezjnBt5CyGEEDmFnb7O0A92cD42iaDxAYzt5IEqxxPwC6KUonvLGmw/doWUtLvb9cYoJZmEnQPq5ziuB+ScOecKtAa2KqWigI7A2lKdnP9H8RKlgjg5OWVvTxQeHs7s2fm+h3CHtLS0O5IwIYQQIi8/7D3PqCW7cHGwZc2Uztzb3IydE2XUfS1rcPNWOiGnrhsdSrGUZBIWAjRTSjVSStkDI4G1WRe11nFaa3ettYfW2gPYBQzWWper7p/58+fj7+9P69atmTRpElmL43br1o05c+bQtWtX3nzzTdauXcuzzz6Lr68vJ06cAOC7774jICCA5s2bs337diMfQwghhMG01nyw5TgzVobTrmFlfpjamaY1Cp+AX+NxnzKzHdHd6tTEHQfbCmVuSLLEkjCtdRowDfgNOAR8q7U+oJSar5QaXFLtGiUpKSnXcOTKlSsBmDZtGiEhIURGRpKUlMS6dX9Pf4uNjeWPP/7gxRdfZPDgwbz11luEh4fTpEkTIKOHLDg4mHfffZdXX33VkOcSQghhvLR0Ey/9EMlbvx1hqG8dvprQweom4BfEyd6G9g2rsOvkVaNDKZYSXTFfa/0z8PNt517Jp2w3szT6y2y4uL/o5b8YUHiZWt7Qr+Chy6zhyNtt2bKFhQsXkpiYyLVr1/Dy8mLQoEEAjBgxosA677//fgDat29PVFRU4XEKIYQodxJvpfHksr1sOnyZJ7o14dneLahQRjesLkkdGlXj3U1Hy9Tq+dazbVHsaYg7e+f503/mPq5UHyo3NEuTycnJTJkyhdDQUOrXr8+8efNITv57JQ5n54JX7HdwcAAyJv2npaWZJSYhhBBlR0xCCo8GhbD/fBz/N7Q1Yzqa5+dTedShcVX0RgiJulZmFqotf0lYIT1WucyrBPPiSiyUrITL3d2dhIQEVq1axbBhw/Is6+rqSnx8fInFIoQQomw5FXOTsZ8Hczk+mY/H+NGrjCQWRvGtXxl72wrsOnlVkjBrkzUnLEvfvn1ZsGABEydOxNvbGw8PD/z9/fO9f+TIkUycOJHFixezatWq0ghZCCGEhQo7fZ3HvgyhglIsn9iRtg2qGB2SxXO0s8G3fmV2n7pmdChFprLe1isr/Pz89O3rZx06dAhPT8/iV1bCPWHlzV1/zkIIIYrstwMXeWr5XmpXciRofAAe7gVPXRF/e2fDEd7fcpx9c3vj6mgZ88KUUmFa6zyX37K+pXWFEEIIC/XlX1FM/iaMVnXcWP1EJ0nAiqltwyqYNByIvmF0KEVi3UlY16IvqCqEEEKUFJNJ8++fDzF37QF6etZk2WMdqebiYHRYZY533UoA7D9XNka5rHtOWPcXjI5ACCGElUtJS2fWdxH8tC+aRwIbMneQFzayBMVdcXdxoG5lJyLOSxImhBBCiALEJaYy6etQdp+6xgv9WjLp3sZWuwekufjUq0TEuVijwygS6x6OFEIIIQxyPjaJYR/9xZ4z13lvpC+Pd20iCZgZeNerxOmricQlphodSqEkCRNCCCFK2YHoOP71wQ4u3kjmqwkdGOJb1+iQyg2fupUB2F8GhiStOwk7sdFsVSmlGDNmTPZxWloa1atXZ+DAgWZrQwghRNm37egVHvxoJ7YVFKuf6ERgk2pGh1SuZE3O31cGhiStOwk7tclsVTk7O2dv0g3w+++/U7du8X6zka2JhBCifPsu9CwTgkKoX7Ui30/pTPOarkaHVO5UqmhHvSpOHLlo+bvQWHcSZmb9+vVj/fr1ACxfvpxRo0ZlXwsODqZTp060bduWTp06ceTIEQCCgoIYPnw4gwYNonfv3ly4cIF7770XX19fWrduzfbt2w15FiGEEOajtWbxpmM8uyqCjo2r8d3kQGpVcjQ6rHKrWQ0Xjl6SJMyqjBw5khUrVpCcnExERAQdOnTIvtayZUu2bdvG3r17mT9/PnPmzMm+tnPnTr788ks2b97MsmXL6NOnD+Hh4ezbty/XVkhCCCHKntR0Ey98v593fj/K/e3q8vk4f4tZzb28al7TlZMxN0lLNxkdSoHK3xIVR36C+AtFLx+6pPAyrrWhxaBCi/n4+BAVFcXy5cvp379/rmtxcXGMHTuWY8eOoZQiNfXvtzZ69epF1apVAfD392fChAmkpqYydOhQScKEEKIMu5mSxtRle9h65ApP9WjKzF7N5Q3IUtC0hgu30kycuZZI4+ouRoeTL+vpCUu6DrGncn/BneeSrv+jZgYPHsysWbNyDUUCvPzyy3Tv3p3IyEh++uknkpOTs685O/+9LcW9997Ltm3bqFu3LmPGjOGrr776R/EIIYQwxuX4ZEYs2cn2YzH8+35vnu7dQhKwUpI11+7opQSDIylY+esJK0KPVbaNL0DPf5u1+QkTJlCpUiW8vb3ZunVr9vm4uLjsifpBQUH53n/69Gnq1q3LxIkTuXnzJnv27OGRRx4xa4xCCCFK1vHLCYz7IpirCbf49BE/uresYXRIVqVpjYzer+OX44FaxgZTAOvpCSsl9erVY/r06Xecf+6553jhhRfo3Lkz6enp+d6/detWfH19adu2LatXr86zLiGEEJYrJOoaD/zvL5JT01n5eEdJwAzg7GBL3cpOFt8TprTWRsdQLH5+fjo0NDTXuUOHDuHp6Vn8ykqgJ6w8u+vPWQghrMTP+y8wY2U49So7ETQ+gAbVKhodktUa90Uwl26k8Mv0ewyNQykVprX2y+ua9IQJIYQQZvDZn6eYumwP3nUrsfqJTpKAGax5TVdOXEkg3WS5nU3lb05YcTS6z+gIhBBClHEmk+a19Yf4fMcp+rWuxaIRvjja2RgdltVr7O7MrTQT0bFJ1K9qmQmxdSdhTXoaHYEQQogyLDk1nae/Defn/RcZ39mDlwa0wqaCvAFpCbJ6Ik9fTZQkTAghhChPYhNvMfGrUEKirvPSAE8eu6ex0SGJHDyqZSz/dPraTbrgbnA0eZMkTAghhCims9cSGfdFMGevJ/HBQ+0Y4FPb6JDEbWq5OWJvW4HTVxONDiVfkoQJIYQQxbD/XBzjg0JITTfxzaMdCGhU1eiQRB4qVFA0qFqR01dvGh1Kvqz67cgPwz80a30uLne3NcJHH32UvTJ+UFAQ0dHR5gxLCCGEmWw5cpkRS3biYFuB1U8ESgJm4RpWrWjRPWFWnYT9b9//jA4BgMmTJ2evin83SVhaWlpJhCWEECKHlSFneOzLUBq5O7NmSiea1nA1OiRRiIbVnDl9NRFLXRPVqpOwkrB161YGDhyYfTxt2rTsbYo8PDx4/vnnCQgIICAggOPHjwMwb9483n77bVatWkVoaCijR4/G19eXpKQkwsLC6Nq1K+3bt6dPnz5cuJCxOXm3bt2YM2cOXbt25b333iv15xRCCGuhteadDUd4fvV+Ojd1Z+XjgdRwczQ6LFEEDao6kZSazrWbt4wOJU+ShJUyNzc3goODmTZtGjNmzMh1bdiwYfj5+bF06VLCw8OxtbXlySefZNWqVYSFhTFhwgRefPHF7PKxsbH88ccfPPPMM6X9GEIIYRVupZl45tt9LN58nBF+9flsrB8uDjKduqyoU9kJgPOxSQZHkrdy9yfpzeA3OXztcJHLj/91fKFlWlZtyfMBz/+TsLKNGjUq+78zZ84ssOyRI0eIjIykV69eAKSnp1O79t9v4IwYMcIsMQkhhLjTjeRUnvgmjB3Hr/JMr+ZM69EUpWQNsLIkKwmLjk3Cp15lg6O5U7lLwvJzPuE8F25euON86KXc+1DWdq5NXZe6d92Ora0tJpMp+zg5OTnX9Zx/gQv7y6y1xsvLi507d+Z53dnZ+a7jFEIIkb/o2CTGfxHCiSsJ/Gd4Gx5oX8/okMRdqJvdE5ZcSEljlLskrDg9Vt5ferN/7H6ztt+wYUMOHjxISkoKycnJbNq0iS5dumRfX7lyJbNnz2blypUEBgbecb+rqyvx8fEAtGjRgitXrrBz504CAwNJTU3l6NGjeHl5mTVmIYQQfzsYfYMJQSHcTEnjywkBdG5qmQt9isJVrmiHk50N0TIcWb6lpaXh4OBA/fr1efDBB/Hx8aFZs2a0bds2V7mUlBQ6dOiAyWRi+fLld9Qzbtw4Jk+ejJOTEzt37mTVqlU89dRTxMXFkZaWxowZMyQJE0KIErL92BWe+GYPLg62fPdEIC1ruRkdkvgHlFLUqexosUmYstTXNvPj5+enQ0NzDyEeOnQIT0/PYtdlzp6wffv2MXHiRIKDg/Mt4+HhQWhoKO7uZfO3qrv9nIUQoiz4LvQsL3y/n6Y1XPhivD+1KzkZHZIwgzGf7eZGUio/TutSeOESoJQK01r75XVN3o40g48++ohRo0bx2muvGR2KEEKIYtJa8+7Gozy7KoLAJtX4bnKgJGDlSL0qTjInzBI90eYJs9QzefJkJk+eXGi5qKgos7QnhBDCPFLTTcz5fj/fhZ1jWPt6/Pt+b+xspH+iPKldyYmYhBSSU9NxtLMxOpxcrDoJm+I7xegQhBBCGCQ+OZUpS/ew/VgM0+9rxoyezWQJihdxgiMAACAASURBVHKoppsDAFfiU6hftaLB0eRm1UmYEEII63QxLpnxQSEcuxTPwmE+POhX3+iQRAmp4Zqxu8FlScKEEEIIYx2+eIPxX4RwIymVz8f5c2/z6kaHJEpQddesnjDLmxcmSZgQQgirseN4DJO/DqOigw3fTg7Eq04lo0MSJaxG5nDk5fgUgyO5k1XPPgy9vsOs9b3++ut4eXnh4+ODr68vu3fvzrdsUFAQ0dHRZm1fCCFE/laGnGHs58HUruzImimdJQGzEtWcHaig4PINK0vClFJ9lVJHlFLHlVKz87g+WSm1XykVrpT6UynVqiTjuV1Y3F9mq2vnzp2sW7eOPXv2EBERwcaNG6lfP/85BneThKWlpf3TMIUQwuqYTJo3fz3M86v3E9ikGque6JS9p6Ao/2wqKNxdHLhsgcORJZaEKaVsgA+AfkArYFQeSdYyrbW31toXWAi8U1LxlLQLFy7g7u6Og0NGt6e7uzt16tRh/vz5+Pv707p1ayZNmoTWmlWrVhEaGsro0aPx9fUlKSkJDw8PYmJiAAgNDaVbt24AzJs3j0mTJtG7d28eeeQRgoKCuP/+++nbty/NmjXjueeeM+qRhRDC4iWnpvPk8r38b+sJHurQgM/H+ePmaGd0WKKU1XBzsLrhyADguNb6pNb6FrACGJKzgNb6Ro5DZ6BsLd+fQ+/evTl79izNmzdnypQp/PHHHwBMmzaNkJAQIiMjSUpKYt26dQwbNgw/Pz+WLl1KeHg4Tk4F/0YWFhbGjz/+yLJlywAIDw9n5cqV7N+/n5UrV3L27NkSfz4hhChrrsSnMHLJLn6OvMCL/T15fWhrWQPMStV0dbTI4ciSnJhfF8iZHZwDOtxeSCk1FXgasAd6/NNGd1zdzNVbl4tcfu2FFYWWqWZfg87VCg7NxcWFsLAwtm/fzpYtWxgxYgQLFizA1dWVhQsXkpiYyLVr1/Dy8mLQoEFFjg9g8ODBuRK1++67j0qVMuYytGrVitOnTxc49CmEENbm6KV4xn8RwtWbKXz0cHv6eNUyOiRhoBpuDuw7F2d0GHcoySQsrxXv7ujp0lp/AHyglHoIeAkYe0dFSk0CJgE0aNDgroKJT4sjIf3GHecvpOTuRXKxccPV9u4ma9rY2NCtWze6deuGt7c3H3/8MREREYSGhlK/fn3mzZtHcnLeY9K2traYTCaAO8o4OzvnOs4a8sxqU+aKCSHE37Yfu8KUb/bgaG/Dt48H4lOvstEhCYNVd3Xk6s0U0tJN2FpQb2hJJmHngJzdM/WAgmairwD+l9cFrfUSYAlkbOBdUKOF9Vjl9HHUWzzu8WyRyxfkyJEjVKhQgWbNmgEZQ4YtWrQgIiICd3d3EhISWLVqFcOGDQPA1dWV+Pj47Ps9PDwICwujX79+rF692iwxCSGEtVkefIaXfoikWQ0XPhvnT12ZgC+Aas72aA1xSalUc3Eo/IZSUpJJWAjQTCnVCDgPjAQeyllAKdVMa30s83AAcIwyKiEhgSeffJLY2FhsbW1p2rQpS5YsoXLlynh7e+Ph4YG/v392+XHjxjF58mScnJzYuXMnc+fO5dFHH+WNN96gQ4c7Rm2FEEIUIOsNyI+3naRr8+q8/1BbXGUCvshUxdkegGs3b1lUEqa0Lrm58Eqp/sC7gA3wudb6daXUfCBUa71WKfUe0BNIBa4D07TWBwqq08/PT4eGhuY6d+jQITw9PYsdnzl7wqzB3X7OQghRkpJupTNzZTi/HrjImI4NmTuolUUNOQnj7Tgew+hPd7NyUkc6NK5Wqm0rpcK01n55XSvRFfO11j8DP9927pUc308vyfaFEEKUb5fjk5n4ZSgR5+N4ZWArxnf2kE24xR2qVPy7J8ySWPW2Re0rdTI6BCGEEHfpYPQNJn4VyvXEW3wyxo+erWoaHZKwUNVcMpOwREnCLIZflc5GhyCEEOIu/HbgIjNXhuPmaMe3jwfSuq5sQSTyV7lixvzA69ITJoQQQtwdrTUfbj3B2xuO4FOvMp+MaU8NN0ejwxIWzsHWBhcHW65KEiaEEEIUX3JqOi98v581e88zuE0dFg7zwdHOxuiwRBlR1dleesKEEEKI4rocn8zjX4ex90wss3o3Z2r3pjIBXxRLFWd7i+sJs+p3eK/8932z1RUVFUXr1q1znZs3bx5vv/222doQQghrFHk+jqHv7+DwhXg+ergd03o0kwRMFFs1Z3uuW9jEfKtOwmI++KDU25QthoQQouh+jbzA8I92ArDqiUD6tq5tcESirKpS0Z5rCZKEWZ1u3boxZ84cunbtynvvvcdPP/1Ehw4daNu2LT179uTSpUtARs/ZhAkT6NatG40bN2bx4sXAnb1sb7/9NvPmzQNg8eLFtGrVCh8fH0aOHFnqzyaEECVBa837m48x+Zs9tKjlyg/TOuNVR96AFHevqrOdxQ1HypywUhIbG8sff/wBwPXr19m1axdKKT799FMWLlzIf/7zHwAOHz7Mli1biI+Pp0WLFjzxxBMF1rtgwQJOnTqFg4MDsbGxJf4cQghR0pJT03luVQRr90Uz1LcOCx6QCfjin6tc0Z6UNBPJqekW8+ep3CVhF994g5RDh4tc/vSYRwot4+DZklpz5hRYJr/5CVnnR4wYkX3u3LlzjBgxggsXLnDr1i0aNWqUfW3AgAE4ODjg4OBAjRo1snvJ8uPj48Po0aMZOnQoQ4cOLfRZhBDCkl2+kczEr8OIOBfLc31b8ETXJjL/S5iFm1PGWmE3klMlCSttt86fJy06+o7ziSEhuY5t69TBvm7dYtdfrVo1rl+/nuvctWvXshMsZ2fn7PNPPvkkTz/9NIMHD2br1q3ZQ4sADg5/byxqY2NDWloatra2mEym7PPJycnZ369fv55t27axdu1a/u///o8DBw5ga2s1/1uFEOXInjPXmfx1GAkpaXz0cHv6eNUyOiRRjrg5ZvxsvJGURg1Xg4PJVO5+WhfWY5XToZaeeB4+ZJZ2XVxcqF27Nps2beK+++7j2rVr/Prrr0yfPp0vvvgiV9m4uDjqZiZ6X375ZaF116xZk8uXL3P16lVcXFxYt24dffv2xWQycfbsWbp3706XLl1YtmwZCQkJVK5c2SzPJIQQpeXbkLO89EMktSo58tWjAbSs5WZ0SKKcyeoJi09ONTiSv5W7JMxIX331FVOnTuWZZ54BYO7cuTRp0uSOcvPmzWP48OHUrVuXjh07curUqQLrtbOz45VXXqFDhw40atSIli1bApCens7DDz9MXFwcWmtmzpwpCZgQokxJTTfx2rqDfLnzNF2auvP+Q22pnLnZshDmlN0Tlmw5qxQorbXRMRSLn5+fDg0NzXXu0KFDeHp6Frsuc/aEWYO7/ZyFECIvVxNSmLJ0D7tPXWPiPY14vm9LbG3kpX1RMo5diqfXom38d1RbBrWpU2rtKqXCtNZ+eV2TnjAhhBClLvJ8HI9/HUZMQgqLRrThX23rGR2SKOdyTsy3FFadhLlPnWp0CEIIYXV+DD/P86sjqFLRnlWTO+FdT9b/EiXPzTEzCUuynOFIq07Cqj85zegQhBDCaqSbNAt/PczH207i71GFD0e3p7qrQ+E3CmEGjnYVsLNR0hMmhBDCusQlpjJt+R62H4vh4Y4NeGWgF/a2Mv9LlB6lFG6OdtxIkiRMCCGElTh88QaPfx1GdGwS/77fm1EBDYwOSVgpNyc7i3o7UpIwIYQQJebH8PPMXr0fF0dblk/siJ9HVaNDElbMzdHWonrCrLsvODSk8DJFpJTKXh8Mcm+yXVxRUVEsW7bMTJEJIUTpu5VmYt7aA0xfEU7rum6sf7KLJGDCcK6OdhY1J8y6k7Cw0MLLFJGDgwPff/89MTEx/7guScKEEGXZpRvJPPTJLoL+iuLRLo1YNrEjNdwcjQ5LCNycpCesXLK1tWXSpEksWrTojmtXrlzhgQcewN/fH39/f3bs2AFkrJw/ZswYevToQbNmzfjkk08AmD17Ntu3b8fX15dFixYRFBTEtGl/v8k5cOBAtm7dCmRsl/Tiiy/Spk0bOnbsmL3hd35tCiFESdp98ioDFv/JwQs3+O+otrw8sBV2sgCrsBBujpY1J0z+ZpjR1KlTWbp0KXFxcbnOT58+nZkzZxISEsLq1at57LHHsq9FRESwfv16du7cyfz584mOjmbBggXcc889hIeHM3PmzALbvHnzJh07dmTfvn3ce++92YlcQW0KIYS5aa35dPtJHvp0N26OtvwwtXOprkouRFG4OcnbkSVrx59w9WrRy6/9sfAy1apB5y6FFnNzc+ORRx5h8eLFODk5ZZ/fuHEjBw8ezD6+ceMG8fHxAAwZMgQnJyecnJzo3r07wcHBxdr/0d7enoEDBwLQvn17fv/99wLbdHW1kK3jhRDlRkJKGs+vimD9/gv08arJ28Pb4Jq5MKYQlsTZ3paUNBNp6SaL2CKr/CVh+Ym/AQkJd56/EJ372MUFXN3uupkZM2bQrl07xo8fn33OZDKxc+fOXIlZFqVUgceQMdRpMpmyj5OTk7O/t7Ozy77HxsaGtLS0QtsUQghzOX45gcnfhHHySgKz+7Xk8Xsb5/nvmBCWwNnBBoDE1HTcJAkrAUXoscr28f/g8SfM2nzVqlV58MEH+eyzz5gwYQIAvXv35v333+fZZ58FIDw8HF9fXwB+/PFHXnjhBW7evMnWrVtZsGABFy5cyO4pA/Dw8ODDDz/EZDJx/vx5goODC42joDaFEMIcftl/gVnf7cPRzoZvHu1Ap6buRockRIGc7DOTsJT07G2MjGR8GlgOPfPMM7nekly8eDGhoaH4+PjQqlUrPvroo+xrAQEBDBgwgI4dO/Lyyy9Tp04dfHx8sLW1pU2bNixatIjOnTvTqFEjvL29mTVrFu3atSs0hoLaFEKIf+JWmon5Px3kiaV7aFbTlZ+e7CIJmCgTnO0z+p4Sb1nG5Pzy1xNmkIQcQ501a9YkMTEx+9jd3Z2VK1fmeV/z5s1ZsmRJrnN2dnZs2rQp17mlS5cW2u6wYcMYNmxYoW0KIcTdOnc9kWnL9hJ+NpaxgQ2ZM8ATB1sbo8MSokiye8JupRscSQZJwoQQQhTJpkOXePrbfZhMmg9Ht6O/d22jQxKiWP7uCZMkzHjt/Qxt/m5X1BdCiNKUmm7i7d+O8PG2k7Sq7caHo9vh4e5sdFhCFNvfPWEyHGk8P3+jIxBCCIt2IS6Jacv2Enb6Og93bMBLA1rhaCfDj6JsqijDkUIIIcqCrUcuM3NlOLfSTCwe1ZbBsviqKONkOFIIIYRFS0s3sWjjUT7YcoKWtVz5YHQ7mlR3MTosIf4xSxuOLNISFUqpnnmcG2v+cIQQQhjp0o1kHvp0Nx9sOcFI//r8MLWzJGCi3MherNVCesKKuk7YK0qp/ymlnJVSNZVSPwGDSjKw0hD800mz1WVjY4Ovry+tW7dm+PDh2UtUuLgU/o/Xu+++m2tJCyGEMMLmw5fo99529p+LY9GINix4wEfmf4lyxdE2a7HWMtQTBnQFTgDhwJ/AMq31sBKLqpSErI8yW11OTk6Eh4cTGRmJvb19sRZHlSRMCGGklLR0Xv3pABOCQqnp5shPT3bhX23rGR2WEGZXoYKior1NmesJqwJ0ICMRSwEaKtkcLF/33HMPx48fz3Vu69at2RttA0ybNo2goCAWL15MdHQ03bt3p3v37gAsX74cb29vWrduzfPPP1+qsQshrMuJKwn864O/+GJHFOM7e7BmSiea1pDhR1F+VbS3ITG1bCVhu4BftNZ9AX+gDrCjxKIqw9LS0vjll1/w9vYuUvmnnnqKOnXqsGXLFrZs2UJ0dDTPP/88mzdvJjw8nJCQEH744YcSjloIYW201nwbepaBi//kQlwSn431Y+4gLxl+FOWek72NxQxHFvXtyJ5a6zMAWusk4Cml1L0lF9bd2/7tUWLOJhReMNOa/+wptIx7fRfuebB5gWWSkpKyN8i+5557ePTRR4scQ04hISF069aN6tWrAzB69Gi2bdvG0KFD76o+IYS43Y3kVF5aE8nafdEENq7GuyN9qenmaHRYQpQKZ3tbixmOLGoS1l4pFae1jgNQSlUGqpZcWOZ342oSCddS7jgffSw217FLVQfcqjkVu/6sOWH5sbW1xWQyZR8nJyfnWU5rXey2hRCiqPaeuc5TK/YSHZvMs31aMLlrE2wqyOwSYT2cLGhOWFGTsLla6zVZB1rrWKXUXMDixskK67HK6YPJm5n6UY8SjOZvDRs25ODBg6SkpJCcnMymTZvo0qULAK6ursTHx+Pu7k6HDh2YPn06MTExVKlSheXLl/Pkk0+WSoxCiPLLZNJ8vO0k/9lwhJpujnz7eCDtG1YxOiwhSl1GT1jZGo7Ma+5YofcqpfoC7wE2wKda6wW3XX8aeAxIA64AE7TWp4sYU5lSv359HnzwQXx8fGjWrBlt27bNvjZp0iT69etH7dq12bJlC//+97/p3r07Wmv69+/PkCFDDIxcCFHWXYxLZtZ3+/jzeAwDfGrzxr+8qeRkZ3RYQhjCyd6GmIQ7R8aMoIoy/KWU+hyIBT4ANPAkUEVrPa6Ae2yAo0Av4BwQAozSWh/MUaY7sFtrnaiUegLoprUeUVAsfn5+OjQ0NNe5Q4cO4enpWehz3K40e8LKg7v9nIUQxlkXEc2LayK5lWZi7qBWjPCvj7zcLqzZjBV72XMmlm3PdS+V9pRSYVprv7yuFfXtyCeBW8AK4FsgCZhSyD0BwHGt9Umtdda9ubp0tNZbtNZZC2TtAmRhGiGEMIMbyak8vTKcacv20sjdmZ+n38PIgAaSgAmr51QGJ+Z7Ai0yy9uSsVr+QMCngHvqAmdzHJ8jY62x/DwK/FLEeMzCf4BHaTYnhBClIvjUNWauDOfijWRm9GzGtO5NsbUp6u/cQpRvzvY2ZW5O2FJgFhAJmAopmyWvX7fyHPtUSj0M+JGxMn9e1ycBkwAaNGhQxOYLFzCosdnqEkIIo91KM/HO70f5eNsJGlatyKrJgbRtIJPvhcjJyd6GpNR0tNaG9wwXNQm7orX+qZh1nwPq5ziuB0TfXihzc/AXga5a6zxnymmtlwBLIGNOWDHjEEKIcu/YpXimrwjn4IUbjAqoz0sDWuHsUNR/4oWwHg62FdAaUtM19rZlIwmbq5T6FNhExrZFAGitvy/gnhCgmVKqEXAeGAk8lLOAUqot8DHQV2t9uTiBCyGEyFh64qudUfz7l8O4ONjyySN+9GpV0+iwhLBYDpmbeKekpWNva+wwfVGTsPFAS8COv4cjNZBvEqa1TlNKTQN+I2OJis+11geUUvOBUK31WuAtwAX4LrNL8IzWevBdPYkQQliZi3HJPLtqH9uPxdCjZQ3efMCH6q4ORoclhEVzsMtIvFLSTLgaHEtRk7A2WuuibYaYg9b6Z+Dn2869kuP7nsWt05zOHL5Mg5Y1zFbfxYsXmTFjBiEhITg4OODh4cG7775L8+ZFX0BWCCEKo7Vmzd7zzFt7gFvpJl4b2prRHeTNRyGKwjG7J6yoU9xLTlGTsF1KqVY51/gqD84diTFbEqa15l//+hdjx45lxYoVAISHh3Pp0qVCkzCtNVprKlSQt5eEEAW7HJ/Mi2si+f3gJfwaVuHt4W3wcHc2OiwhyozsnrBU45epKOpP/S5AuFLqiFIqQim1XykVUZKBlTVbtmzBzs6OyZMnZ5/z9fXlnnvu4a233sLf3x8fHx/mzp0LQFRUFJ6enkyZMoV27dpx9uxZNmzYQGBgIO3atWP48OEkJGRsRD579mxatWqFj48Ps2bNMuT5hBDG+2lfNL0XbeOPo1d4aYAnKx8PlARMiGJyyJwHlpxadnrC+pZoFOVAZGQk7du3v+P8hg0bOHbsGMHBwWitGTx4MNu2baNBgwYcOXKEL774gg8//JCYmBhee+01Nm7ciLOzM2+++SbvvPMO06ZNY82aNRw+fBilFLGxsXm0LoQoz64mpPDKjwdYv/8CbepX5j/D29C0hovRYQlRJuWcmG+0IiVhZWk/x1P7L3IzLrnI5SP/jCq0jHMlRxp517qreDZs2MCGDRuy94pMSEjg2LFjNGjQgIYNG9KxY0cAdu3axcGDB+ncuTMAt27dIjAwEDc3NxwdHXnssccYMGAAAwcOvKs4hBBl06+RF3lxzX7ik9N4rm8LJt3TWBZeFeIfyOoJK0tzwsq85MRb3Eq6c4XcG1cTcx3bO9niWNG+2PV7eXmxatWqO85rrXnhhRd4/PHHc52PiorC2dk5V7levXqxfPnyO+oIDg5m06ZNrFixgvfff5/NmzcXOz4hRNkSm3iLeWsP8EN4NF513Fg20ZcWtYx+l0uIss/BruxNzC8zitNj9dePB+k0pJVZ2u3Rowdz5szhk08+YeLEiQCEhITg5ubG559/zujRo3FxceH8+fPY2dndcX/Hjh2ZOnUqx48fp2nTpiQmJnLu3Dnq1KlDYmIi/fv3p2PHjjRt2tQs8QohLNfvBy/x4pr9XLt5i5k9mzOlexPspPdLCLPI7gmzgIn55S4JM4pSijVr1jBjxgwWLFiAo6Nj9hIVlStXJjAwEAAXFxe++eYbbGxsct1fvXp1goKCGDVqFCkpGevhvvbaa7i6ujJkyBCSk5PRWrNo0aJSfzYhROmISUhh3toDrIu4QMtarnw+zp/WdSsZHZYQ5Ypj5tuRydITVr7UqVOHb7/99o7z06dPZ/r06Xecj4yMzHXco0cPQkJC7igXHBxsviCFEBYna92v+esOkpiSzqzezXm8q/R+CVESsifmS0+YEEJYt/OxScz5fj9/HL1CuwaVWTjMh6Y1ZO6XECVFJuZbiHot3I0OQQhhpUwmzTe7T/PmL4fRwLxBrRgT6IFNBVn1XoiSJBPzLYQ5tywSQoiiOnElgdmrIwiJus49zdx541/e1K9a0eiwhLAKfy/WKsORQghhNVLTTXyy/STvbjyGk50Nbw9vwwPt6sqej0KUoqy5lmnp2uBIJAkTQohSEXb6Oi+u2c/hi/H0a12LV4d4UcPV0eiwhLA6NhUUFVTGL0VGkyRMCCFKUFxiKm/+dphlu89Qp5IjS8a0p7fX3e3AIYQwDzubChaRhFn1+89/fbfUbHXZ2Njg6+ub/bVgwYI8y73yyits3LjRbO0KISyT1pofw89z3ztbWRF8hse6NOL3p7tKAiaEBbC3qUCqDEcaa+eq5XQaPtosdTk5OREeHl5gmfT0dObPn2+W9oQQluv01Zu89EMk24/F0KZeJYLGB8iiq0JYEDtb6QmzCh4eHsyfP58uXbrw3XffMW7cuOw9JmfPnk2rVq3w8fFh1qxZBkcqhPinbqWZeH/zMXov2kb4mVjmD/Hi+ymdJQETwsLYVlAWkYRZdU+YOSUlJeHr65t9/MILLzBixAgAHB0d+fPPPwH49ddfAbh27Rpr1qzh8OHDKKWIjY0t/aCFEGYTfOoac9bs5/jlBAZ41+aVQa2o6SYT74WwRHYyHFkytgQt4fLpk0Uuv/LV2YWWqdGwMd3HTSqwTEHDkVnJWE5ubm44Ojry2GOPMWDAAAYOHFi0gIUQFuVyfDILfj7M93vPU6+KE1+M86e7rEEohEWzt5DhyHKXhOUn7vIl4mOu3HH+3MHc+ze6ulenUo2aZm3b2dn5jnO2trYEBwezadMmVqxYwfvvv8/mzZvN2q4QouSkppv4audp3v39KClpJqZ1b8qU7k2oaG81/6wKUWbJcGQJKazHKqf/jBjIMyvXlWA0+UtISCAxMZH+/fvTsWNHmjZtakgcQoji23XyKnN/PMCRS/F0a1GduYO8aOR+5y9bQgjLZClLVJS7JMwot88J69u3b77LVADEx8czZMgQkpOT0VqzaNGi0ghTCPEPXLqRzOvrD7F2XzT1qjixZEx7erWqKSveC1HGZLwdKXPCyo309Lz3oIqKisp1HBQUlP19cHBwCUYkhDCXW2kmvthxisWbjpFq0jx1XzOmdGuCY+ZGwEKIssVOhiOFEMLy/XkshrlrIzlx5SY9PWvw8sBWNKwmQ49ClGUyHGkBAoeNMjoEIYSFOhVzkzd+PsTvBy/RoGpFPh/nR4+W5n1pRwhhDDvbCiQl5T2CVZqsOgkz12r5Qojy40ZyKu9vPs4XO05hb1OB5/q2YELnRjL0KEQ5Ym8jw5FCCGEx0k2aFSFneGfDUa4l3mJ4+3rM6t2CGrLgqhDljm0FGY4UQgiL8NfxGOavO8jhi/EEeFTly0GtZKshIcoxO9sKpMnbkUIIYZyozHlfGw5eol4VJz4c3Y5+rWvJkhNClHN2NopbFtATZtUbeJ/ZY74lIpRSjBkzJvs4LS2N6tWrF7odUVBQENOmTTNbHIXp1q0boaGhpdaeEJYoLjGVN34+RK9Ff7DjeAzP9mnBxqe70t+7tiRgQlgBOxmONN7ZvSE0aBdglrqcnZ2JjIwkKSkJJycnfv/9d+rWrWuWuoUQ5pGcms7XO0/z/pbj3EhOZVi7ejzbR+Z9CWFt7GyVRSzWatU9YebWr18/1q9fD8Dy5csZNervJTCCg4Pp1KkTbdu2pVOnThw5cuSO+9evX09gYCAxMTFcuXKFBx54AH9/f/z9/dmxYwcA8+bNY+zYsfTu3RsPDw++//57nnvuOby9venbty+pqakAzJ8/H39/f1q3bs2kSZPQOvcfNpPJxNixY3nppZcA2LBhA4GBgbRr147hw4eTkJBQIp+REEYwmTRr9p7jvv/8wes/H6Jtg8r8/NQ9vDW8jSRgQlghS1knTJIwMxo5ciQrVqwgOTmZiIgIOnTokH2tZcuWbNu2jb179zJ//nzmzJmT6941a9awYMECfv75Z9zd3Zk+fTozZ84kJCSE1atX89hjj2WXPXHiBOvXr+fHH3/k4Ycfpnv37uzfvx8nJ6fsJHDatGmEhIRk986tW/f3HplpaWmMHj2a5s2bxrDMPgAAG9xJREFU89prrxETE8Nrr73Gxo0b2bNnD35+frzzzjsl/GkJUTr+PBbDoPf/ZObKfVRxtmPpYx0IGh+AZ203o0MTQhjEUpKwcjcceXLXdm5ejSly+f3r1xRaxrmaO4073lNoOR8fH6Kioli+fDn9+/fPdS0uLo6xY8dy7NgxlFLZPVYAW7ZsITQ0lA0bNuDmlvGDYePGjRw8eDC7zI0bN4iPjwcyetzs7Ozw9vYmPT2dvn37AuDt7Z29TdKWLVtYuHAhiYmJXLt2DS8vLwYNGgTA448/zoMPPsiLL74IwK5duzh48CCdO3cG4NatWwQGBhb6vEJYsoPRN1jw62G2Hb1C3cpOvDfSl0E+dahQQeZ8CWHtKiiFyfjRyPKXhOUnOf4Gt27eOcR242J0rmN7ZxccXe/+N+TBgwcza9Ystm7dytWrV7PPv/zyy3Tv3p01a9YQFRVFt27dsq81btyYkydPcvToUf6/vTsPj7K89z/+vifrEEIgG0s2QHYIEgiL4FrB9Qh6hKp1KYp6vKz91aNtr3q1PaXt1erxp+f86nJ+yrEW9Ce1hR4Fl1YUwQ0XgiCyEyFAIIQsJJBlkszM/ftjQgj7IJl5JpnP67pyzTPzPHme73AzyTf3fT/fu7CwEAgMF3766ae43e4TrpGQkACAy+UiLi6ubSKxy+XC6/Xi8Xi4//77KSoqIicnh7lz5+LxeNq+f/LkyaxYsYKHH36YxMRErLVMmzaNP//5z9/6fYt0hKVlrwIwve/N3/oce2sa+Y9l2/iftaX0SIzjF9cO5/YL8kiIVbFVEQlwGU6YpuOELpeEBdNjdcQnf3yWKXN+0KHXv+uuu0hJSSE/P5+VK1e2vV5bW9s2Ub/9It4AeXl5PPHEE9xwww0sWrSIkSNHcsUVV/DMM8/wk5/8BIB169YxZsyYoGI4knClp6dTV1fH4sWLmTlzZtv+OXPm8OGHHzJr1ixee+01Jk2axA9+8AOKi4sZNGgQDQ0NlJaWMmTIkHP4lxAJr4rDTTy7opiFn+8GA/deNJD7Lx1ESrc4p0MTkQgTKT1hmhPWwbKzs/nRj350wus//elPeeSRR5gyZQo+34nrVQ0dOpRXXnmFWbNm8c033/DUU09RVFTE6NGjGTFiBM8991zQMfTs2ZN77rmH/Px8rr/+esaPH3/CMQ899BBjx47l9ttvJy0tjfnz53PLLbcwevRoJk2axJYtW87ujYs4pLahhcf/sYWLH1/By5/t4oaCLFb8+FIeuWa4EjAROSmXAX8E9ISZSOiOOxuFhYX2+DpXmzdvZvjw4Wd9rlD0hHVl3/bfWSRYZzMcWd/k5U+f7OT5D3dw2OPluvP78a9TBzMwo3uowxSRTu4/393GH5ZvZ+ej14S8NqAxZo21tvBk+7rccKSIdG2eFh+vfL6b/1pRTFV9M1OHZ/LQtKGM6Ke7HUUkOK7WxMtacLI+c1QnYTkFJw7TiUhkavH5WbymlKeWb6es1sPk89L48ZVDGZvby+nQRKSTOXKTtN9aXDiXhUV1EtZR1fJFJHRafH5e+3Ivz6woZnd1AwW5PXly1vlMHpTudGgi0kkdKVXj9OT8kCZhxpirgD8AMcAL1trHjtt/MfB/gNHAzdbaxaGMR0Q6j2avn//5spRnVxazp7qR/KwUXrijkMuHZ2p9RxE5J6ZdT5iTQpaEGWNigGeBaUApsNoYs9Rau6ndYbuB2cCPQxWHiHQuXh8s/Hw3z64oZm9NI+dnp/Dr6SO5bKiSLxHpGO3nhDkplD1hE4Bia+0OAGPMq8AMoC0Js9aWtO5zfu0AEQmpI3c+norXB8vWH2bVl+kcqvuavD4+7pvRxPC8w9SZUt7Yf+zx51LQVUSim6ur94QBWcCeds9LgYmnONYRte/uImVaXoecq3v37scsej1//nyKiop45plnOuT8Il1Vixc+2xTHe6vjOViXTFbvBm6d6mVYrs/Ru5ZEpOs60hPWlZOwk/34/Fbv1hhzL3AvQG5u7rnEdIzDy3d3WBImIqd3fM9VXZOXhZ/v4oWPdnLgcBOFeb2YMHUvQ3N9zOinXi4RCR1jImNifigr5pcCOe2eZwP7TnHsaVlr51lrC621hRkZGR0SXDi98cYbTJw4kYKCAqZOnUp5eTkAc+fO5Yknnmg7btSoUZSUlFBSUsKwYcO4++67GTVqFLfeeivvvfceU6ZMYfDgwXzxxRcA1NfXc9dddzF+/HgKCgpYsmQJACUlJVx00UWMHTuWsWPHsmrVqvC/aZFTqK5v5j+WbWXKY+/z+7e3MLh3d165eyKL7ruAYXnq/RKR0GsbjnQ4CwtlT9hqYLAxZgCwF7gZ+F4Ir+eoxsbGY9Z2rK6uZvr06QBceOGFfPbZZxhjeOGFF3j88cd58sknT3u+4uJiFi1axLx58xg/fjwLFy7k448/ZunSpfz+97/n9ddf53e/+x3f+c53ePHFF6mpqWHChAlMnTqVzMxM3n33XRITE9m+fTu33HILx68yIBJue2sa+e8Pd/Dq6t14WvxcObI39186iPNzejodmohEmS4/HGmt9RpjHgDeIVCi4kVr7UZjzG+AImvtUmPMeOA1oBdwnTHm19bakedy3Zo3vqF5X33Qxx94fv0Zj4nvl0TP68477TFut5t169a1PT8yJwygtLSUm266ibKyMpqbmxkwYMAZrzlgwADy8/MBGDlyJJdffjnGGPLz8ykpKQFg2bJlLF26tK03zePxsHv3bvr168cDDzzAunXriImJYdu2bWe8nkioFB+o47kPvuH1tXsBuL4gi/suGcigzGSHIxORaHV0Yr6zcYS0Tpi19m3g7eNe+7d226sJDFOGnPegB39N0wmvN++sPea5q2cCsb0SO/TaP/zhD3nooYeYPn06K1euZO7cuQDExsbi9x+9MdTj8bRtJyQkHI3J5Wp77nK58Hq9AFhr+dvf/sbQoUOPud7cuXPp3bs3X331FX6/n8TEjn0/IsH4cvdB5n2wg3c27Sch1sVtk/K45+KBZPV0Ox2aiEQ501aioov2hDnlTD1W7ZX+7COyH7sohNEE1NbWkpWVBcCCBQvaXu/fvz9vvvkmAF9++SU7d+48q/NeeeWVPP300zz99NMYY1i7di0FBQXU1taSnZ2Ny+ViwYIF+Hy+jnsz0qVs+LgEgFEX9u+Q8/n8lmUb9/PCxztZs+sgPRJjeeCyQcye3J+07glnPoGISBjEREPFfAmYO3cus2bNIisri0mTJrUlWzfeeCMvvfQSY8aMYfz48QwZMuSszvvLX/6SBx98kNGjR2OtbUvq7r//fm688UYWLVrEZZddRlJSUijelnRSd/7jTgD+dNWfOuyc9U1e/lq0hxc/2cme6kZyU7sx97oRzCrMISlBP2ZEJLJEQ52wqNK+RhjA7NmzmT17NgAzZsxgxowZJ3yP2+1m2bJlJz3fhg0b2rbnz5/ftt2/f/+2fW63m+eff/6E7x08eDDr1x+d6/boo48G/T5EzkZZbSMLVu1i4ee7OOTxMi6vFz+/ZjjTRvRp+0tTRCTSmK4+MV9Euq4Ne2v548c7eeOrffit5epRfZlz0QDG5vZyOjQRkTOKhmWLIl7y5R1X+FWkq2vx+Xln435eWrWLL0qqSYqP4Y4L+nPnlP7kpHbrkGtoKSIRCQcNR0YAVcsXObMDhz28+sUeXvl8F+WHmshJdfPza4bz3fE5pLjjnA5PROSsuSKkYn6XScKstW1jvNLxnL6NV8LLWsvaPTW8tKqEt74uo8VnuXhIBo/+cx6XDMnUfC8R6dSMesI6TmJiIlVVVaSlpSkRCwFrLVVVVao3FgU8LT7eXF/GS5+WsL60luSEWG6blMftk/IYmNHd6fBERDqES3XCOk52djalpaVUVFQ4HUqXlZiYSHZ2WOrqigN2VNTx6uo9LF5TSnV9M4Mzu/Pb60dxQ0EW3VViQkS6GA1HdqC4uLiglgISiSZHCrEe79rmmwBY+8EOlhdX8n75YTa9uZEYAxN7J3Ntfl9GpyVhvJaS1aVt39dRBV1FRJymifki0uGOFGKFo8nW8UoO17GmIpb1lZup90JGQixX5xrGZxp6xDcADWw9GDh2WOqwMEQtIhJebXXC/Gc4MMSUhIl0UW/F/6Vt2+93UVHRj31lA6mpyQR8ZKSX8r3MXAan+Pl7wiI+Osk5Zl7YcVX1RUQihXrCRKTDtV+KyFrLxn2HWLymlCXr9nKwoYWcVDepA76gT58SFs54rm3I8rtKtkQkiqhYq4iERGVdE6+v3cviNaVs2X+Y+BgX00b05qbxOVw4KJ05y152OkQREUe5XIFH9YSJyDlr9vpZsfUAi9eUsmLLAbx+y/nZKfx2xkiuO78fPbvFOx2iiEjE0NqRInJOrLVsKjvE39bs5fV1e6mubyYjOYE5Fw7gxnHZDOmd7HSIIiIRSSUqRORb2VPdwJJ1e1mybh/bD9QRH+Ni6ohMZo7L5uLBGcTGuJwOUUQkoh2ZmK9irSJyRpV1Tby1vowl6/by5e4aAMb378Vvrx/FP+X3pVeShhtFRIKlnjAROa26Ji/LNu5nybp9fFxcic9vGdYnmZ9eNZTp5/cju1c3p0MUEemUtHakiJygodnLii0VvP11Gcu3lONp8ZPV082/XDyQ6WP6MaxPD6dDFBHp9FyamC8iAPVNXpZvOcDfvy5jxdYDeFr8pCXFM3NcNtePyWJsbi9crtAsTK+liEQkGqlOmEgUO+xp4f0tB3hrfRkfbKugyesnIzmBWeNyuCa/LxMGpBITosSrfUFXEZFopIr5IlGmur6Z97cc4B8byvhwWyXNPj+9eyRwy4Rcrsnvy7i8XiFLvERE5CijifkiXd/Oynre21TOu5vKKdpVjd9C35REbpuUx7Wj+1CQE7qhRhEROTn1hImE2K7b7wAg7+WXwnZNv9+ydk8N724q573N5RQfqANgeN8ePHDZIKaN6MOorB5tf4WJiEj4HZ0TpiRMpFOra/KyqriS5ZsPsHxLOZV1zcS6DBMHpnLbxFwuH96bnFSVkxARiRRtd0f6nY1DSZjIWbLWsrX8MCu3VvDB1gqKdlXT4rMkJ8Ry6bBMpg7P5NKhmaS445wOVURETkJ1wkSCsXRJ4HH6DEfDOORp4ZPtlYHEa1sF+w95ABjWJ5k5Fw7kkiEZjMvrRXyslgwSEYl0qpgvEsG8Pj9fldayqriSj7ZXsmb3QXx+S3JiLBcNTufSIZlcPCSDPimJTocqIiJnydX697LmhIlEAL8/MMS46psqVhVX8vnOauqavBgDI/r24L5LBnLp0EwKcnpqgWwRkU5OPWEiDrLWsqe6kU++qeST4ko+/aaKqvpmAAakJzFjTD+mDErngoFpWhxbRKSLUYkK6bS+fus1APKvvcHhSILn91uKK+r4Ymc1q0uq+WJnNWW1gXldmckJXDwkg8nnpTFlUDr9erodjlZERELJaO1IkdBp8fnZFJfK+oQMti8oomhXNTUNLUAg6Ro/IJWJA1KZfF4a52V0V90uEZEoorUjJfIUzQs8Ft7rbBxBOFKI9YgaVwKb4lPZHJfGxvg0Nsan48mcBkDW18Vc0FTB6OZKzm86QL+99ZgtR793N+Et6CoiIs7ScKSEzp+uDTze+ZazcYRIs9fP5rhUNsWnsSk+jc3xqeyNTQbAZf0MbKnl6oadDNv0GaOqdpA9epjDEYuISCTRxHyRIPgs7GiEDWtL+br0EOv2HGTDvkM0t/Zy9e6RQEFOL27P7UlBTk/ys1PoFh/4b73r9jsga5h6uURE5Bgq1tqJfbSxmk1lh6j8Rwkj+vbghofHAvCXX/8MgJt+9ZiT4XUeRwqxtmr2w7YG2Fhv2FBv2FBn2FzvwmNdsO4rEl2W/CSYnWkZ091SkGzpm+AF6qG2FC5xtqCriIh0Dlo7UqKWtZbKuma21hi2NsDWBsPGesO2BmixgQ9G9xjLiCT4XnIdo+KbGZXdk4FuiNX8eREROUcajpSocMjTwvbyw2zdX8fW/YfYWn6YbeV1VNc3AzEApCbFM7JfD+b0S2FUVg9G9UshN7UbLpdp7S3r5viyRSIi0nVoYr50GT6/ZV9NI99U1LGzsp4dFfXsqKxjR0V9Wy0ugKT4GIb0SeaKEb0Z0juZYX2SGdInmfTuCQ5GLyIi0caoJ0wi0dfFXih/7YTXfRaqmqDcYygpP8T+llgqNr3B3kYoazRtw4gASTGWrG6WoW7L7VeNYGjvZIb0TiarpzvQuyUiIuKgI7+KNCdMIoLPb6lqdLGlPp7yRhflHihvMpR7Al8VTeBrS7TSicHSxw3Zbkthqp9+bku225LltqTEHb3zJP/SQY69JxERkZNpmxPmcFdYSJMwY8xVwB8ITP55wVr72HH7E4CXgHFAFXCTtbYklDFFI0+Lj4pV8ylriGG/J4b9DTHsbzzy5WJ/YwwHPDF4bd9jvi89wUd2kpfxGV5yuvnISfKSk+SjobKJzLQMxvxT51m2SERE5IguPzHfGBMDPAtMA0qB1caYpdbaTe0OmwMctNYOMsbcDPw7cFOoYuoK7IvXUueP45A/gVpfAjX+BKq83aj0uanyuqnyuak8PIEq24OqXy6kyuemzh8P9DnmPN1cXvokeOjT3TAps5k+iT76Uk7fxGZyM9PITvLRLfbk/zu/PuwFjSqKiEgnZVyBx648MX8CUGyt3QFgjHkVmAG0T8JmAHNbtxcDzxhjjHV6kPYcWAtenx+v39Ls8+P1WVp8fpq9ra95/TQ0e2ls9tHQ7KOxxde67aWhbTvwvLaxhUONrY+elsBjw334cZ302i78pMZ4SLeVpJlDjEksJy22kbSYRjJiG+gTW0/f2Dp6x9aT7GoODBm2r6pfNA9wQ+Fdp3+TJ5kzFolUpFVERE4mGtaOzAL2tHteCkw81THWWq8xphZIAypDGNdpHSnE6l2575jXm0r/X9v2kRTI3eJnZxk8eG8mb/Uch89Mxo+Ln/3879/6+sZAt7gYuiXEkuKOo0diLOnd4xmYkdT6PC7w6G7d744jo3sCad0T6OmOC0x87+LLFomIiJyLaChRcbIBq+PfbTDHYIy5F7gXIDc399wj62DdfY2MbNyNq6WGGOsnJjYZF35c1mL8PhJjDSm94onBTwyWeOPjyu/fTbf4GLrFx+COj6FbfCzd4mNIiHW13TorIiIiHS/W5eKOC/IY0a+Hs3GE8NylQE6759nAvlMcU2qMiQVSgOrjT2StnQfMAygsLAxp2nrRyFQuGpkKU/sft2dy29bxyxY9dMKyRb8KZYjRRUVaRUSkg8XHuvjNjFFOh3GKyUUdYzUw2BgzwBgTD9wMLD3umKXA91u3ZwLvd+b5YCIiIiLBCllPWOscrweAdwiUqHjRWrvRGPMboMhauxT4I/CyMaaYQA/YzaGKR0RERCSShLROmLX2beDt4177t3bbHmBWKGMQERERiUSqmC9nLf9aFWkVERE5V6GcEyYiIiIip6CeMDmq8F6nIxAREYkaSsK6IhVpFRERiXgajhQRERFxgOlsZbkKCwttUVGR02GIiIiInJExZo21tvBk+9QTJiIiIuIAJWEiIiIiDlASJiIiIuIAJWEiIiIiDlASJiIiIuIAJWEiIiIiDlASJiIiIuIAJWEiIiIiDlASJiIiIuIAJWEiIiIiDlASJiIiIuIAJWEiIiIiDlASJiIiIuIAJWEiIiIiDlASJiIiIuIAJWEiIiIiDlASJiIiIuIAJWEiIiIiDlASJiIiIuIAJWEiIiIiDjDWWqdjOCvGmApgl8NhpAOVDscgR6k9IovaI7KoPSKL2iOyhKM98qy1GSfb0emSsEhgjCmy1hY6HYcEqD0ii9ojsqg9IovaI7I43R4ajhQRERFxgJIwEREREQcoCft25jkdgBxD7RFZ1B6RRe0RWdQekcXR9tCcMBEREREHqCdMRERExAFKwk7DGHOVMWarMabYGPOzk+xPMMb8pXX/58aY/uGPMnoE0R4PGWM2GWPWG2OWG2PynIgzWpypPdodN9MYY40xuiMshIJpD2PMd1s/IxuNMQvDHWM0CeLnVa4xZoUxZm3rz6xrnIgzGhhjXjTGHDDGbDjFfmOMeaq1rdYbY8aGKzYlYadgjIkBngWuBkYAtxhjRhx32BzgoLV2EPCfwL+HN8roEWR7rAUKrbWjgcXA4+GNMnoE2R4YY5KB/wV8Ht4Io0sw7WGMGQw8Akyx1o4EHgx7oFEiyM/HL4C/WmsLgJuB/wpvlFFlPnDVafZfDQxu/boX+L9hiAlQEnY6E4Bia+0Oa20z8Cow47hjZgALWrcXA5cbY0wYY4wmZ2wPa+0Ka21D69PPgOwwxxhNgvl8APyWQDLsCWdwUSiY9rgHeNZaexDAWnsgzDFGk2DawwI9WrdTgH1hjC+qWGs/BKpPc8gM4CUb8BnQ0xjTNxyxKQk7tSxgT7vnpa2vnfQYa60XqAXSwhJd9AmmPdqbA/w9pBFFtzO2hzGmAMix1r4ZzsCiVDCfjyHAEGPMJ8aYz4wxp+sZkHMTTHvMBW4zxpQCbwM/DE9ochJn+/ulw8SG4yKd1Ml6tI6/lTSYY6RjBP1vbYy5DSgELglpRNHttO1hjHERGKKfHa6Aolwwn49YAsMtlxLoJf7IGDPKWlsT4tiiUTDtcQsw31r7pDHmAuDl1vbwhz48OY5jv8vVE3ZqpUBOu+fZnNhd3HaMMSaWQJfy6bo85dsLpj0wxkwFfg5Mt9Y2hSm2aHSm9kgGRgErjTElwCRgqSbnh0ywP6+WWGtbrLU7ga0EkjLpeMG0xxzgrwDW2k+BRALrGEr4BfX7JRSUhJ3aamCwMWaAMSaewMTJpccdsxT4fuv2TOB9q8JroXLG9mgd/nqeQAKm+S6hddr2sNbWWmvTrbX9rbX9CczRm26tLXIm3C4vmJ9XrwOXARhj0gkMT+4Ia5TRI5j22A1cDmCMGU4gCasIa5RyxFLgjta7JCcBtdbasnBcWMORp2Ct9RpjHgDeAWKAF621G40xvwGKrLVLgT8S6EIuJtADdrNzEXdtQbbH/wa6A4ta74/Yba2d7ljQXViQ7SFhEmR7vANcYYzZBPiAn1hrq5yLuusKsj0eBv7bGPOvBIa+ZuuP+NAwxvyZwDB8euscvF8BcQDW2ucIzMm7BigGGoA7wxab2lxEREQk/DQcKSIiIuIAJWEiIiIiDlASJiIiIuIAJWEiIiIiDlASJiIiIuIAJWEiIiIiDlASJiIiIuIAJWEiEvWMMbcZY74wxqwzxjxvjIlxOiYR6fqUhIlIVGtdMuYmYIq1dgyBavK3OhuViEQDLVskItHucmAcsLp1uSs3oLVHRSTklISJSLQzwAJr7SNOByIi0UXDkSIS7ZYDM40xmQDGmFRjTJ7DMYlIFFASJiJRzVq7CfgFsMwYsx54F+jrbFQiEg2MtdbpGERERESijnrCRERERBygJExERETEAUrCRERERBygJExERETEAUrCRERERBygJExERETEAUrCRERERBygJExERETEAf8fARnslGC8LR4AAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 720x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"def maclaurin(e):\n",
" t1 = 2 * np.sqrt(1 - e**2) / e**3\n",
" t2 = (3 - 2 * e**2)\n",
" t3 = np.arcsin(e)\n",
" t4 = (6 / e**2) * (1 - e**2)\n",
" return t1 * t2 * t3 - t4\n",
"\n",
"def maclaurin_e(a, c):\n",
" return np.sqrt(1 - (c**2/a**2))\n",
"\n",
"def maclaurin_y(omega, rho):\n",
" return omega**2 / (np.pi * G * rho)\n",
"\n",
"# Compute e and the Maclaurin x, and y values of each object\n",
"df['e'] = maclaurin_e(df['a'], df['c'])\n",
"df['mcx'] = maclaurin(df['e'])\n",
"df['mcy'] = maclaurin_y(df['omega'], df['rho'])\n",
"\n",
"# Color setup\n",
"cmap = cm.get_cmap('tab20')\n",
"df['n'] = np.arange(len(df))\n",
"colors = [cmap(n) for n in df['n'].values]\n",
"\n",
"# Maclaurin Plot\n",
"f, ax = plt.subplots(figsize=(10, 5))\n",
"ax.set_title('Maclaurin Spheroids')\n",
"ax.set_ylabel(r'$\\frac{\\Omega^2}{\\pi G \\rho}$')\n",
"ax.set_xlabel('e')\n",
"\n",
"# Plot objects\n",
"df.plot('e', 'mcx', 'scatter', c=colors, marker='+', s=400, ax=ax, \n",
" colorbar=False, legend=True)\n",
"\n",
"# PLot curve\n",
"e = np.linspace(0.01, 1, 1000)\n",
"ax.plot(e, maclaurin(e))\n",
"\n",
"# legend\n",
"from matplotlib.lines import Line2D\n",
"legend_entries = [Line2D([0], [0], marker='+', color=cmap(i), label=label,\n",
" markerfacecolor=cmap(i), markersize=10) \n",
" for i, label in df[['n', 'name']].values]\n",
"plt.legend(legend_entries, df['name'])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Jacobi ellipsoid"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x7f29519115d0>"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAmoAAAG5CAYAAAA+kBhjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeXhU1f3H8ffJHgIIJMgOERdkC4EECSIluIAoIlUQKeKCSK17larVqtSftu7ibi1StEVQURRBLaKgqERMkCLIjgmrQEIIgayTOb8/kowJ2SZhJjNJPq/n4TH33nPP/d7LOPlyzrnnGGstIiIiIuJ/AnwdgIiIiIhUTomaiIiIiJ9SoiYiIiLip5SoiYiIiPgpJWoiIiIifkqJmoiIiIifUqImIg2OMcYaY06r4tgkY8zSOtY7xxjzSMnPicaY3WWObTDGJNYpYPeuPdQYs9md2ESk6VCiJiIeZ4xZYYyZ6otrW2vnWmtHVHbMGBNdkuQdPe7PBDfq7W2tXeHxgH+tf6W1toe36heRhinI1wGIiPhAK2utw9dBiIjURC1qIuI1xpjWxpjFxpiDxpjMkp87lznexhjzL2PM3pLjH5Q5doMxZpsx5pAxZpExpuNx1V9kjNlhjEk3xjxpjAkoOe9aY8zXXriXVGPM+SU/zzDGLDDGvG2MyTbGrDHG9CtT9h5jzJ6SY5uNMeeV7A81xswsud+9JT+Hlhw7vqu1f0m92caYt4GwMseiSp7l4ZLns7L0/kWkcdH/2CLiTQHAv4BuQFcgF3ixzPF/A82A3sDJwLMAxphzgb8DVwAdgDRg/nF1/xaIBwYAlwJTvHUTVbgUeBdoA7wFfGCMCTbG9ABuAQZaa1sAI4HUknPuBxKAWKAfcBbwl+MrNsaEAB9Q/HzalFzn8jJF7gJ2A22BdsB9gNYDFGmElKiJiNdYazOste9Za3OstdnAo8AwAGNMB2AUcKO1NtNaW2it/bLk1EnAbGvtGmttPvBnYLAxJrpM9Y9baw9Za3cCM4GJtQgtvaQ1qvRPzzrcXoq1doG1thB4huIWrwSgCAgFehljgq21qdba7WXu62Fr7QFr7UHgr8DkSupOAIKBmSXPZQHwfZnjhRQnsN1Kjq+0WrhZpFFSoiYiXmOMaWaM+YcxJs0YcwT4CmhljAkEugCHrLWZlZzakeJWNACstUeBDKBTmTK7yvycVnKOu6Ksta3K/NlYi3MrXN9a66S4haujtXYbcAcwAzhgjJlfptu23H1VE3dHYM9xyVfZ854EtgFLS7p/761D/CLSAChRExFvugvoAQyy1rYEflOy31Cc6LQxxrSq5Ly9FHeXFhc2JgKIBPaUKdOlzM9dS86pT67rl4wP61wag7X2LWvtORTfgwUeLyla7r6oOu59QCdjjDmuLCX1Z1tr77LWdgcuAe4sHQcnIo2LEjUR8aYWFI9LO2yMaQM8VHrAWrsP+AR4ueSlg2BjTGki9xZwnTEmtmSw/d+A76y1qWXq/lPJeV2A24G36+F+yoozxlxmjAmiuAUtH0gyxvQwxpxbEncexfdfVHLOPOAvxpi2xpgo4EHgP5XUvQpwALcZY4KMMZdRPJ4NAGPMaGPMaSWJ3JGS+osqqUdEGjglaiLiLZbisWPhQDqQBHx6XJnJFI+32gQcoDjhwVr7OfAA8B7FrUunAlced+6HQAqwFlgCvF6L2A4fN4/anbU4t+z1JwCZJfdxWcl4tVDgMYrv+ReKX5K4r+ScR4BkYB3wI7CmZF851toC4DLg2pL6JwDvlylyOrAMOEpxUveyN+d4ExHfMRp/KiKeZoxZQ/Gg+Q9qLNwAGWNmAKdZa6/ydSwi0ripRU1EPMoY0xvoCfzg61hERBo6JWoi4jHGmMeBpcA91tq0msqLiEj11PUpIiIi4qfUoiYiIiLipxrkouxRUVE2Ojra12GIiIiI1CglJSXdWtu2Luc2yEQtOjqa5ORkX4chIiIiUiNjTJ3H7KrrU0RERMRPKVETERER8VNK1ERERET8VIMco1aZwsJCdu/eTV5enq9DERFp9MLCwujcuTPBwcG+DkWkUWs0idru3btp0aIF0dHRFK9TLCIi3mCtJSMjg927d3PKKaf4OhyPS0nLJGlHBgndI4nr1trX4UgT12gStby8PCVpIiL1wBhDZGQkBw8e9HUoHpeSlsmkWUkUOJyEBAUwd2qCkjXxqUY1Rk1JmohI/Wis37dJOzIocDhxWih0OEnakeHrkKSJa1SJWk3at2+PMabGP+3bt/d1qI3Kxl1HfR1Cg5ac+Y2vQ2gydm464OsQxMcSukcSEhRAoIHgoAASukf6OiRp4ppUorZ//36PlhP3bNpzzCP1JCYm8t///rfcvpkzZ3LTTTd5pH5/lZL1rUfre/TRR+nduzcxMTHExsby3XffVVl2zpw57N2716PX92e7N6d7rK5ffvmFK6+8klNPPZVevXpx0UUXsWXLFo/VL94R1601c6cmcOeIHur2FL/QaMaoSeM3ceJE5s+fz8iRI1375s+fz5NPPunDqBqWVatWsXjxYtasWUNoaCjp6ekUFBRUWX7OnDn06dOHjh07un0Nh8NBUFDT/mqx1vLb3/6Wa665hvnz5wOwdu1a9u/fzxlnnFHjudZaAgKa1L+j/Upct9ZK0MRv6JvAw8aOHUtcXBy9e/fmtdde83U4jcq4ceNYvHgx+fn5AKSmprJ3717OOeccnnzySQYOHEhMTAwPPfSQ63jPnj254YYb6N27NyNGjCA3Nxcobp0rXYYsPT2d0rVjN2zYwFlnnUVsbCwxMTFs3bq1/m/Ui/bt20dUVBShoaEAREVF0bFjRx5++GEGDhxInz59mDZtGtZaFixYQHJyMpMmTSI2Npbc3Fyio6NJTy9udUpOTiYxMRGAGTNmMG3aNEaMGMHVV1/NnDlzuOyyy7jwwgs5/fTTufvuu311yz6xfPlygoODufHGG137YmNjGTp0aLWf1ZtuuokBAwawa9culi5dyuDBgxkwYADjx4/n6NHiIQT33nsvvXr1IiYmhunTp/vk/kSk/jTKf/becccdrF279oTqKP0FVCo2NpaZM2fWeN7s2bNp06YNubm5DBw4kMsvv5zIyMY3xmFdajZZxwrdLr9yw6Eay5wUEUxMdIsqj0dGRnLWWWfx6aefcumllzJ//nwmTJjAZ599xtatW1m9ejXWWsaMGcNXX31F165d2bp1K/PmzeOf//wnV1xxBe+99x5XXXVVldd49dVXuf3225k0aRIFBQUUFRW5fY+18U3GF2QUuD8eatG++TWWiQw5mSGR51ZbZsSIETz88MOcccYZnH/++UyYMIFhw4Zxyy238OCDDwIwefJkFi9ezLhx43jxxRd56qmniI+Pr/H6KSkpfP3114SHhzNnzhzWrl3LDz/8QGhoKD169ODWW2+lS5cu7t2wB/384y8cy3J/fsX1X6fWWCbipDBO6Vv1WNb169cTFxdXYf/SpUur/Kxu3ryZf/3rX7z88sukp6fzyCOPsGzZMiIiInj88cd55plnuOWWW1i4cCGbNm3CGMPhw4fdvi8RaZgaZaLmS88//zwLFy4EYNeuXWzdurVRJmpVOZZfRG6Bs8L+9OzySV14SAARoYG1rr+0+7M0UZs9ezZvvfUWS5cupX///gAcPXqUrVu30rVrV0455RRiY2MBiIuLIzU1tdr6Bw8ezKOPPsru3bu57LLLOP3002sd44nIdmRxtOhIhf378neV224e2JIWQSfVuv7mzZuTkpLCypUrWb58ORMmTOCxxx6jRYsWPPHEE+Tk5HDo0CF69+7NJZdcUqu6x4wZQ3h4uGv7vPPO46STimPs1asXaWlpPknUqpKXU0BBrqPC/iMZOeW2Q8KDCGsW4pFrLl26tMrPardu3UhISAAgKSmJn376iSFDhgBQUFDA4MGDadmyJWFhYUydOpWLL76Y0aNHeyQuEfFfjTJRq6rlqzavk69YsaLW112xYgXLli1j1apVNGvWjMTExEa7UkJ1LV/HW5i0n98mtPPIdceOHcudd97JmjVryM3NZcCAAcydO5c///nP/P73vy9XNjU11dXFBxAYGOjq+gwKCsLpLE4oy/4d/e53v2PQoEEsWbKEkSNHMmvWLM49t/pWqrqoqeWrrH+kPsnvo//ksWsHBgaSmJhIYmIiffv25R//+Afr1q0jOTmZLl26MGPGjCo/t1U9N4CIiIhy28c/e4ejYlJUH6pr+Tretx/+xNmX9jrha/bu3ZsFCxZU2G+trfKzWvb5WWu54IILmDdvXoU6Vq9ezeeff878+fN58cUX+eKLL044XhHxXxqj5kFZWVm0bt2aZs2asWnTJpKSknwdUqPTvHlzEhMTmTJlChMnTgRg5MiRzJ492zWGZ8+ePRw4UH23YnR0NCkpKQDlfqHu2LGD7t27c9tttzFmzBjWrVvnpTvxjc2bN5cbd7d27Vp69OgBFI9XO3r0aLnn0aJFC7Kzs13bZZ/be++9V09RNzznnnsu+fn5/POf/3Tt+/7772nZsqVbn9WEhAS++eYbtm3bBkBOTg5btmzh6NGjZGVlcdFFFzFz5swTHuIhIv6vUbao+cqFF17Iq6++SkxMDD169HB1Y4hnTZw4kcsuu8z1Nt2IESPYuHEjgwcPBoqTuf/85z8EBlbdtTp9+nSuuOIK/v3vf5drMXv77bf5z3/+Q3BwMO3bt3eN22osjh49yq233srhw4cJCgritNNO47XXXqNVq1b07duX6OhoBg4c6Cp/7bXXcuONNxIeHs6qVat46KGHuP766/nb3/7GoEGDfHgn/s0Yw8KFC7njjjt47LHHCAsLIzo6mpkzZ9KqVasaP6tt27Zlzpw5TJw40fXyzCOPPEKLFi249NJLycvLw1rLs88+W+/3JiL1y1hrfR1DrcXHx9vSN/ZKbdy4kZ49e1Z7Xm26Phvic/FXnuz6bIo83fUpVfNU12dT4c73roiAMSbFWlvzW1mVaFJdn+3auZcsuFtO3HNmp4iaC0mV4k4629chNBmde0T5OgQRkXKaVNfnL7/84usQmqSeXZr7OoQGLb71EF+H0GR0PfNkX4cgIlJOk2pRExEREWlIlKiJiEijlpKWyUvLt5GSlunrUERqrcknajNmzPB1CCIi4iUpaZlMmpXE00s3M2lWkpI1aXCafKL217/+1dchNH7L/+7rCBq0gy+86OsQmoxv353r6xDEw5J2ZFDgcOK0UOhwkrQjw9chidRKk0/UPCk1NZU+ffr4Ogz/8+VjHqsqMDCQ2NhY15/HHqtd3R988AE//fSTa7vs4uz+Kv2llzxWV2Wf0RkzZvDUU0957BoN2aoFFVcCqCt3P6sPPvggy5Yt89h1pbyE7pGEBAUQaCA4KICE7k1nST9pHJrUW5/S8IWHh9d5NnaHw8EHH3zA6NGj6dVLc2VVx+FwEBSkr4cT4c5ntaioiIcffrieImqa4rq1Zu7UBJJ2ZJDQPZK4bq19HZJIrahFzcMcDgfXXHMNMTExjBs3jpycnJpPkhP28MMPM3DgQPr06cO0adNcExYnJiZy3333MWzYMB5//HEWLVrEn/70J2JjY9m+fTsA7777LmeddRZnnHEGK1eu9OVt+FTZZ/Xcc8/x0UcfMWjQIPr378/555/P/v37geIWuClTppCYmEj37t15/vnngYqtdU899ZRrDOjzzz9Pr169iImJ4corr6z3e/Mn0dHRPPzww5xzzjm8++67XHvtta5lu+69917Xc5o+fbqPI2084rq15ubhpylJkwapUf6T+Y477qhVq0tiYmKNZWJjY6tc7L2szZs38/rrrzNkyBCmTJnCyy+/3Di/cD+5F3750f3y/7q45jLt+8Ko6rsyc3NziY2NdW3/+c9/ZsKECdxyyy2u5Z4mT57M4sWLueSSSwA4fPgwX375JQBbt25l9OjRjBs3zlWHw+Fg9erVfPzxx/z1r3+tl26oX/72N/I3bnK7fNrkq2ssE9rzTNrfd9+JhFXuWWVmZpKUlIQxhlmzZvHEE0/w9NNPA7Bp0yaWL19OdnY2PXr04A9/+EO19T722GP8/PPPhIaGcvjw4ROKsbaWz3mNA2k73C7/9l/vrbHMyd26M/zaadWWqeqzChAWFsbXX38NwKeffgrAoUOHWLhwIZs2bcIYU+/PSUT8U6NM1NzhdDrZtWsXAGlpaXTp0oWAgBNvYOzSpQtDhhRPUHrVVVfx/PPPN85ErSqH0yBrV8X9aV+X3z6pC7TqVuvqq+pOWr58OU888QQ5OTkcOnSI3r17uxK10l+OVbnssssAiIuLIzU1tdYxeVLBnj049u6tsD/n++/LbQd17EhIp061rr+qZdRK95d9Vrt372bChAns27ePgoICTjnlFNexiy++mNDQUEJDQzn55JNdrW1ViYmJYdKkSYwdO5axY8fWOm5vyDqwn+z0gxX27/5pfbntFlFtOenk2q9WUl3XZ2WfyZYtWxIWFsbUqVO5+OKLGT16dK2vKSKNT6NM1Gpq+fr444+57bbbiI2NZc6cOcycOZMff/yRF154gVGjRp3QtY//RVib9UUblBpavsqZcRLMyPJaKHl5edx0000kJyfTpUsXZsyYQV5enut4RET1S1iFhoYCxYO/HQ6H1+IsqzYtXxvP7EnPTRs9ct3IyEgyM8tPT3Do0CFXElb2Wd16663ceeedjBkzhhUrVpSbyqb0mcGvzy0oKAin0+naX/bvYMmSJXz11VcsWrSI//u//2PDhg31Ngauppavsp6eMJq73l7sxWiKVfaZDAoKYvXq1Xz++efMnz+fF198kS+++MLrsYiIf2tSY9RSU1MZO3Yst99+Oy+88AILFy5k2LBhLFy4kBdeeIHbbruNsWPHnlCrys6dO1m1ahUA8+bN45xzzvFQ9FKV0oQgKiqKo0ePusb7VKZFixZkZ2fXV2h+p3nz5nTo0IHPP/8cKE7SPv3000o/p1lZWXQqabV74403aqy7Xbt2HDhwgIyMDPLz81m8uDjhKW29Hj58OE888QSHDx/m6NGjHryrxuHo0aNkZWVx0UUXMXPmzDq/NCMijUuTStTeeecdBg4cyI8//lih5WzUqFH8+OOPxMfH884779T5Gj179uSNN94gJiaGQ4cO1Th2R2qndNxP6Z97772XVq1accMNN9C3b1/Gjh3LwIEDqzz/yiuv5Mknn6R///6ulwmamjfffJNHHnmE2NhYzj33XB566CFOPfXUCuVmzJjB+PHjGTp0KFFRNS9WHhwczIMPPsigQYMYPXo0Z555JlD8ZuNVV11F37596d+/P3/84x9p1aqVx+/L31T2Wa1OdnY2o0ePJiYmhmHDhvHss8/WU6Qi4s9M6dtxDUl8fLw9fu6rjRs30rNnTx9FJNXyctdnY+fJrk+pXn11fTYW+t4VcY8xJsVaG1+Xc5tUi5r4yLCa36KTqkXdfLOvQ2gyBo+b6OsQRETKUaIm3jf8z76OoEFre+stvg6hyTh7/CRfhyAiUo4SNRERERE/pURNRERExE8pURMRERHxU00yUdu+fTs33XQTLVu2JCAggJYtW3LTTTc12ekavG6795dkatSSv6+5jHjEzjWrfR2CiEg5TS5R++STT4iJiWHWrFlkZ2djrSU7O5tZs2YRExPDJ5984usQG5+fP/dYVcYYJk+e7Np2OBy0bdu2cS+3k5Jccxk3GWO46667XNtlF06vrdTUVN566y0PReYfdv3guaS4rp/VOXPmcMst9fcCSWJiIsdPdyQi/qNJJWrbt29n3Lhx5OTkUFhYWO5YYWEhOTk5jBs3Ti1rfiwiIoL169eTm5sLwGeffeaaPd9d9bVMlD8KDQ3l/fffJz09/YTraoyJmid54rMqItKkErWnn366QoJ2vMLCwhOaEfzNN98kJiaGfv36lfvXtHjOqFGjWLJkCVC8TNfEib/OfbV69WrOPvts+vfvz9lnn83mzZuB4laK8ePHc8kllzBixAj27dvHb37zG2JjY+nTpw8rV670yb3Ut6CgIKZNm1bpZ/zgwYNcfvnlDBw4kIEDB/LNN98AxSsUTJ48mXPPPZfTTz+df/7znwDce++9rFy5ktjYWJ599tkKLUGjR49mxYoVQPHSVffffz/9+vUjISHBtYh7VddsLOryWS1ryZIlDB48mPT09Gr/fq655hpGjBhBdHQ077//PnfffTd9+/blwgsvdH3nPfzwwwwcOJA+ffowbdo0jp/s3Ol0cs011/CXv/wFgKVLlzJ48GAGDBjA+PHjteyXiI94dVVkY8xsYDRwwFrbp4oyicBMIBhIt9YOO9Hr3nHHHZWuk/f1119TVFRU7bmFhYW8+uqrrF+/vtz+2NjYGhd737BhA48++ijffPMNUVFRHDp0qPbBNxSbP4Lsfe6XT36t5jItOkCPS2osduWVV/Lwww8zevRo1q1bx5QpU1yJ1plnnslXX31FUFAQy5Yt47777uO9994DYNWqVaxbt442bdrw9NNPM3LkSO6//36KiorIyclx/1484ZuvISPD/fKLPqy5TGQkDKl5bdmbb76ZmJgY7r777nL7b7/9dv74xz9yzjnnsHPnTkaOHMnGjcUrIqxbt46kpCSOHTtG//79ufjii3nsscd46qmnXGt6zpkzp8prHjt2jISEBB599FHuvvtu/vnPf/KXv/yl2mt6yo6klRzLcL8F8cclC2ssExEZRfeEoTWWq+tnFWDhwoU888wzfPzxx7Ru3Zrf/e53VT6r7du3s3z5cn766ScGDx7Me++9xxNPPMFvf/tblixZwtixY7nlllt48MEHAZg8eTKLFy/mkkuK/39zOBxMmjSJPn36cP/995Oens4jjzzCsmXLiIiI4PHHH+eZZ55xnS8i9ceriRowB3gReLOyg8aYVsDLwIXW2p3GmJO9GUxNSVptyx3viy++YNy4ca51Edu0aVOnehq03EzIP1xx/+Gfy2+HtoLw1nW6RExMDKmpqcybN4+LLrqo3LGsrCyuueYatm7dijGmXAvqBRdc4Po7GThwIFOmTKGwsJCxY8cSGxtbp1g8LvsIVNZysW9v+e3mzaFFyzpdomXLllx99dU8//zzhIeHu/YvW7aMn376ybV95MgR1wL2l156KeHh4YSHhzN8+HBWr15dq/U6Q0JCXGOz4uLi+Oyzz6q9ZosWLep0b7WRl32EgmMVn/WRX8o/65CI5oTV8VnX9bO6fPlykpOTWbp0KS1bFl+7ur+fUaNGERwcTN++fSkqKuLCCy8EoG/fvqSmprrqfOKJJ8jJyeHQoUP07t3blaj9/ve/54orruD+++8HICkpiZ9++okhQ4YAUFBQwODBg+v0DLwhJS2TpB0ZJHSPJK5b3b5HRBoKryZq1tqvjDHR1RT5HfC+tXZnSfkDnrhuVS1fLVu2dH2xVadly5auLpvasNZijKn1eQ2SGy1fLsv+DOf/3aOXHzNmDNOnT2fFihVklGmZeuCBBxg+fDgLFy4kNTWVxMRE17GIiAjXz7/5zW/46quvWLJkCZMnT+ZPf/oTV199tUdjrJYbLV8u/3gFfv8Hj17+jjvuYMCAAVx33XWufU6nk1WrVpVL3kod/7mu7HMeFBSE0+l0befl5bl+Dg4Odp0TGBjoGidY3TU9xZ2Wr1LfvP4SQ6737JJddfmsdu/enR07drBlyxbi44uXB6zuWYWGhgIQEBBQ7lkHBATgcDjIy8vjpptuIjk5mS5dujBjxoxyfz9nn302y5cv56677iIsLAxrLRdccAHz5s3z6LPwhJS0TCbNSqLA4SQkKIC5UxOUrEmj5usxamcArY0xK4wxKcaYKn9TGmOmGWOSjTHJBw8erNPFrrrqKoKDg6stExwcXOexZeeddx7vvPOO68u4UXd9+tiUKVN48MEH6du3b7n9WVlZrgHb1XXFpaWlcfLJJ3PDDTdw/fXXs2bNGm+G63fatGnDFVdcweuvv+7aN2LECF588UXXdtnhAx9++CF5eXlkZGSwYsUKBg4cSIsWLcr9wyc6Opq1a9fidDrZtWsXq1fXPNVFdddsLOryWe3WrRvvv/8+V199NRs2bABO7FmVJmVRUVEcPXqUBQsWlDt+/fXXc9FFFzF+/HgcDgcJCQl88803bNu2DYCcnBy2bNni9vW8KWlHBgUOJ04LhQ4nSTtqMYRApAHydaIWBMQBFwMjgQeMMWdUVtBa+5q1Nt5aG9+2bds6Xeyuu+5yK1H74x//WKf6e/fuzf3338+wYcPo168fd955Z53qkZp17tyZ22+/vcL+u+++mz//+c8MGTKk2i7sFStWEBsbS//+/Xnvvfcqrauxu+uuu8q9/fn888+TnJxMTEwMvXr14tVXX3UdO+uss7j44otJSEjggQceoGPHjsTExBAUFES/fv149tlnGTJkCKeccgp9+/Zl+vTpDBgwoMYYqrtmY1HXz2qPHj2YO3cu48ePZ/v27Sf0rFq1asUNN9xA3759GTt2LAMHDqxQ5s4772TAgAFMnjyZyMhI5syZw8SJE4mJiSEhIYFNmzbV7sa9JKF7JCFBAQQaCA4KIKF7pK9DEvEqc/ybPx6/QHHX5+LKXiYwxtwLhFlrZ5Rsvw58aq19t7o64+Pj7fHz/mzcuJGePXvWGM8nn3zCuHHjKCwsLDcmJDg4mODgYBYsWMCoUaNqrEdqwQtdn02KF7o+a2PGjBk0b96c6dOn+yyG+uKNrs/GzN3vXU/TGDVpaIwxKdba+Lqc6+sWtQ+BocaYIGNMM2AQ4NlXvo4zatQo1q1bx7Rp08qtTDBt2jTWrVunJM0bTjnP1xE0bHF1+n9b6qBL/4otTeJ/4rq15ubhpylJkybBqy1qxph5QCIQBewHHqJ4Gg6sta+WlPkTcB3gBGZZa6ufA4MTa1ETERHP0PeuiHtOpEXN2299TnSjzJPAk96MQ0RERKQh8nXXp4iIiIhUQYmaiIiIiJ9q8onas5/5x9xAjdnLa1/2dQgN2uqPdvg6hCYj67M0X4cgIlJOk0/Unvt8q8fqMsaUmyzX4XDQtm1b19I5TdUr/3vFY3U1b968Tue9+uqrvPlm8Upmc+bMYe/evTWc4T++X5LqsWWJdQAAACAASURBVLoCAwNdC9GPHz/etcapO8915syZ9b8maj3L/nynx+o6/pkev2i9iIg7mnyi5kkRERGsX7+e3NxcAD777DPXzOPuKl1ap65O9PzG6sYbb3QtEVWXRK2xPNfw8HDWrl3L+vXrCQkJqdWkqU0hURMR8TdK1Dxs1KhRLFmyBIB58+YxceKvL74eO3aMKVOmMHDgQPr378+HH34IFCcO48eP55JLLmHEiBEAPPHEE/Tt25d+/fpx7733ApCYmEjptCTp6elER0dXev7kyZNddQNMmjSJRYsWef3e68uKFSvKtVLecsstriV4oqOjueeeezjrrLM466yzXEvgzJgxg6eeeooFCxaQnJzMpEmTiI2NJTc3l5SUFIYNG0ZcXBwjR45k3759QPHzvu+++xg2bBjPPfdcvd+ntw0dOtT1fEpV9Wyff/559u7dy/Dhwxk+fDhQ/Pnu27cvffr04Z577qnX2Bu6jz76iEGDBtG/f3/OP/989u/fD/z6OS3Vp08fUlNTSU1N5cwzz2Tq1Kn06dOHSZMmsWzZMoYMGcLpp5/uWq6rqu+Y1NRUhg4dyoABAxgwYADffvtt/d+0iNSJV6fn8JW/frSBn/Yecbv8hH+sqrFMr44teeiS3jWWu/LKK3n44YcZPXo069atY8qUKaxcuRKARx99lHPPPZfZs2dz+PBhzjrrLM4//3wAVq1axbp162jTpg2ffPIJH3zwAd999x3NmjVza83Qsud/+eWXPPvss1x66aVkZWXx7bff8sYbb9RYR208vvpxNh1yf0mZ6z69rsYyZ7Y5k3vOOvFf+C1btmT16tW8+eab3HHHHSxevNh1bNy4cbz44os89dRTxMfHU1hYyK233sqHH35I27Ztefvtt7n//vuZPXs2AIcPH+bLL7884ZiOt/KdLaTvOup2+YVP17wWaVSX5gy9otIV2CpwOBx88sknXHjhhW6Vv+2223jmmWdYvnw5UVFR7N27l3vuuYeUlBRat27NiBEj+OCDDxg7dqxb9dWnwx9tp2DvMbfLH/jHuhrLhHSMoNUlp1ZbJjc3l9jYWNf2oUOHGDNmDADnnHMOSUlJGGOYNWsWTzzxBE8//XS19W3bto13332X1157jYEDB/LWW2/x9ddfs2jRIv72t7/xwQcfVPkdc/LJJ/PZZ58RFhbG1q1bmThxIsfPRSki/qlRJmpV2Z2Zw57DeRX2f/dz+USoU6swOrduVqdrxMTEkJqayrx587jooovKHVu6dCmLFi1y/Ys5Ly+PnTuLx8RccMEFtGnTBoBly5Zx3XXX0axZcQyl+6tT9vxhw4Zx8803c+DAAd5//30uv/xygoLq5696z9E97Du2r8L+5P3lfyl0iOhAp+a16xZ2V2kr5sSJE2tct3Xz5s2sX7+eCy64AICioiI6dOjgOj5hwgSvxFiVIxm5HD2UX2H/3q2Hy203bxNKy8jwWtdfNnkYOnQo119/fZ3i/P7770lMTKR03d1Jkybx1Vdf+WWiVhVHZh7OwxWfdcHPWeW2A1qFEtQ6rNb1l3Yzl5ozZ44rOdq9ezcTJkxg3759FBQUcMopp9RYX+k6qlC8rvB5552HMYa+ffuSmpoKVP0d07FjR2655RbWrl1LYGCg3yywLiI1a5SJmjstX6Wi711C6mMXe/T6Y8aMYfr06axYsYKMjAzXfmst7733Hj169ChX/rvvviMiIqJcOWNMhXqDgoJwOp1A8RdwWWXPB5g8eTJz585l/vz5rtYhT6pNy1ffN/ry4zU/euzaZZ8DVHwWZZ9dZc+xLGstvXv3ZtWqyltVj3+unuJuyxfASzd+wc2vnuuR6x6fPByvpmdbyttrBHtSTS1fZe2+dyWdHxvqxWiK3Xrrrdx5552MGTOGFStWMGPGDKD65x8aGur6OSAgwLUdEBDgGkNZ1XfMjBkzaNeuHf/73/9wOp2EhdU+8RQR39AYNS+YMmUKDz74oOtfv6VGjhzJCy+84Pol98MPP1R6/ogRI5g9e7Zr4HZp12d0dDQpKSkALFiwoNoYrr32WmbOLF6Nq3dv9xPXhqBbt2789NNP5Ofnk5WVxeeff17u+Ntvv+367+DBgyuc36JFC7KzswHo0aMHBw8edCVqhYWFbNiwwct34L+qe7Zln9ugQYP48ssvSU9Pp6ioiHnz5jFs2DBfhd3gZGVluV40KjssITo6mjVriru516xZw88//1yreqv6jsnKyqJDhw4EBATw73//m6KiIk/chojUAyVqXtC5c2duv/32CvsfeOABCgsLiYmJoU+fPjzwwAOVnn/hhRcyZswY4uPjiY2NdXVjTJ8+nVdeeYWzzz6b9PT0amNo164dPXv25Lrrah4b1lA4HA5CQ0Pp0qULV1xxBTExMUyaNIn+/fuXK5efn8+gQYN47rnnePbZZyvUc+2113LjjTcSGxtLUVERCxYs4J577qFfv37ExsY26YHW1T3badOmMWrUKIYPH06HDh34+9//zvDhw+nXrx8DBgzg0ksv9WHkDcuMGTMYP348Q4cOJSoqyrX/8ssv59ChQ8TGxvLKK69wxhnut7xC1d8xN910E2+88QYJCQls2bLFay3FIuJ5Xl2U3Vs8uSi7N7o+/UFOTg59+/ZlzZo1nHTSST6NxVNdn//73/+44YYbXG+4VSY6Oprk5ORyv/waOk92fUr16qvrs7HQouwi7jmRRdmbfIva7eed7usQPG7ZsmWceeaZ3HrrrT5P0gD+0O8PJ1zHq6++ysSJE3nkkUc8EFHDMvDiaF+H0GS0OK+rr0No8lLSMnlp+TZS0jJ9HYqIX2jyLWoiIlI3nv7eTUnLZNKsJAocTkKCApg7NYG4bq09Vr+Ir6hFrURDTDpFRBoib3zfJu3IoMDhxGmh0OEkaUdGzSeJNHKNJlELCwsjIyNDyZqIiJdZa8nIyPD4NB8J3SMJCQog0EBwUAAJ3SM9Wr9IQ9Ro5lHr3Lkzu3fv5uDBg74ORUSk0QsLC6Nz584erTOuW2vmTk0gaUcGCd0j1e0pQiNK1IKDg92a3VtERPxXXLfWStBEymg0XZ8iIiIijY0SNRERERE/pURNRERExE8pURMRERHxU0rURERERPyUEjURERERP6VETURERMRPKVETERER8VNK1ERERET8lBI1ERERET+lRE1ERETETylRExERr0pJy+Sl5dtIScv0dSgiDU6jWZRdRET8T0paJpNmJVHgcBISFMDcqQladF2kFtSiJiIiXpO0I4MChxOnhUKHk6QdGb4OSaRBUaImIiJek9A9kpCgAAINBAcFkNA90tchiTQo6voUERGvievWmrlTE0jakUFC90h1e4rUkhI1ERHxqrhurZWgidSRuj5FRERE/JQSNRERERE/pa5PERERES9IScskaUcGJiQ8oq51KFETERER8bCycwgGt+54Rl3rUdeniIiIiIeVnUMQMHWtR4maiIiIiIeVnUMQsHWtR4maiIiIiIeVziF454geFGbu3VLXejRGTURERMQLSucQvKUg91hd61CLmoiIiIifUqImIiIi4qeUqImIiIj4KSVqIiIiIn5KiZqIiIiIn1KiJiIiIk1CSlomLy3fRkpapq9DcZum5xAREZFGr+ySTiFBAcydmkBct9a+DqtGalETERGRRq/skk6FDidJOzJ8HZJblKiJiIhIo1d2SafgoAASukf6OiS3eLXr0xgzGxgNHLDW9qmm3EAgCZhgrV3gzZhERMTzUtIySdqRQUL3yAbRnSRNT+mSTg3tc+rtMWpzgBeBN6sqYIwJBB4H/uvlWERExAsa6tgfaXpKl3RqSLza9Wmt/Qo4VEOxW4H3gAPejEVERLyjoY79EWkIfDpGzRjTCfgt8KobZacZY5KNMckHDx70fnAiIuKWhjr2R6Qh8PX0HDOBe6y1RcaYagtaa18DXgOIj4+39RCbiIi4oaGO/RFpCHydqMUD80uStCjgImOMw1r7gW/DEhGR2miIY39EGgKfJmrW2lNKfzbGzAEWK0kTERERKebt6TnmAYlAlDFmN/AQEAxgra1xXJqIiIhIU+bVRM1aO7EWZa/1YigiIiIiDY6vx6iJiIif2vPQt9j8ohrLmdBAOv317HqISBoLTZDsPiVqIiJSKXeStNqUEwFNkFxbWutTRERE6o0mSK4dJWoiIiJSbzRBcu2o61NERETqjSZIrh0laiIiIlKvNEGy+9T1KSIiIuKnlKiJiIiI+CklaiIiIiJ+SomaiIhIE5OSlslLy7eRkpbp61CkBnqZQEREKmVCA12T2a7HwQ8U0Z9A+hz3q8OEBvoiPKkjTTjbsChRExGRSpUuC5WSlskd+sXeaFQ24az+Pv2Xuj5FRKRamkm+cdGEsw2LWtRERKRapb/YCx1O/WJvBDThbMNirLW+jqHW4uPjbXJysq/DEBFpMlLSMvWLXaSOjDEp1tr4upyrFjUREamRZpIX8Q2NURMRERHxU0rURERERPyUEjURERERP6VETURERMRPKVETERER8VNK1ERERET8lBI1ERERET+lRE1ERMSHUtIyeWn5NlLSMn0divghTXgrIiLiIylpmUzSgvdSDbWoiYiI+IgWvJeaKFETERHxkdIF7wMNWvBeKqWuTxERER+J69aauVMTtOC9VEmJmoiIiA9pwXupjro+RURERPyUEjURERERP6VETURERMRPKVETERER8VNK1EREpJymPFN+U7538U9661NERFya8kz5TfnexX+pRU1ERFya8kz5TfnexX8pURMREZemPFN+U7538V/GWuvrGGotPj7eJicn+zoMEZFGKSUts8nOlN+U7128xxiTYq2Nr8u5GqMmIiLlNOWZ8pvyvYt/UteniIiIiJ9SoiYiIiLip5SoiYiIiPgpJWoiIiIifkqJmoiI+DWtFiBNmd76FBERv6XVAqSpU4uaiIj4La0WIE2dEjUREfFbWi1Amjp1fYqIiN+K69aauVMTtFqANFlK1ERExK9ptQBpytT1KSIiIuKnlKiJiIiI+CklaiIiIiJ+yquJmjFmtjHmgDFmfRXHJxlj1pX8+dYY08+b8YiIiOdpQloR7/H2ywRzgBeBN6s4/jMwzFqbaYwZBbwGDPJyTCIi4iGakFbEu7zaomat/Qo4VM3xb621pf8ESwI6ezMeERGpWl1axjQhrYh3+dP0HNcDn1R10BgzDZgG0LVr1/qKSUSkSahry1jphLSFDqcmpBXxAr9I1IwxwylO1M6pqoy19jWKu0aJj4+39RSaiEiT8P6a3eQXOrH82jLmTqKmCWlFvMvniZoxJgaYBYyy1qrNXESknqWkZfJu8i5K/wUcGFi7ljFNSCviPT6dnsMY0xV4H5hsrd3iy1hERJqqpB0ZOJzFaZoBxsV1VuIl4ie82qJmjJkHJAJRxpjdwENAMIC19lXgQSASeNkYA+Cw1sZ7MyYRESnv+HFmlw/Qe10i/sJY2/CGe8XHx9vk5GRfhyEi4rdS0jJrNW6stuVFxH3GmJS6NkT5fIyaiIh4Vl3e4NQ4MxH/pCWkREQaGc1tJtJ4KFETEWlkSsecBRo0t5lIA6euTxERP+KJsWKa20yk8VCiJiLiJzy5bqbGnIk0Dur6FBHxExpbJiLHcytRM8ZEGGMCymwHGGOaeS8sEZGmR2PLROR47nZ9fg6cDxwt2W4GLAXO9kZQIiJNkcaWicjx3E3Uwqy1pUka1tqjalETEfE8jS0TkbLcHaN2zBgzoHTDGBMH5HonJBEREREB91vU7gDeNcbsLdnuAEzwTkgiIiIiAm4matba740xZwI9AANsstYWlh43xlxgrf3MSzGKiPgNrYkpIvXJ7XnUShKz9VUcfhxQoiYijZon5zkTEXGHp+ZRMx6qR0TEb2meMxGpb55K1KyH6hER8Vua50xE6puWkBIRcVNTmudMY/FE/IOnErVUD9UjIuLXmsI8ZxqLJ+I/3F1CKtkYc7MxptL/U621l3k2LBER8RWNxRPxH+6OUbsS6Ah8b4yZb4wZaYzRCwQi4pdS0jJ5afk2UtIyfR1Kg6SxeCL+w1jr/nsAJQuzjwZeAZzAbOA5a+0h74RXufj4eJucnFyflxSRBkLddp6hMWoinmOMSbHWxtflXLfHqBljYoDrgIuA94C5wDnAF0BsXS4uIuJplXXbKdGovaYwFk+kIXArUTPGpACHgdeBe621+SWHvjPGDPFWcCIitVXabVfocKrbTkQaPHdb1MZba3dUdkAvEoiIP2lKU2iISOPn7ssEU40xrUo3jDGtjTGPeCkmEZETEtetNTcPP01Jmog0eO4maqOstYdLN6y1mRSPVRMRERERL3E3UQs0xoSWbhhjwoHQasqLiGiaDBGRE+TuGLX/AJ8bY/5Vsj0FeNM7IYlIY6BpMkRETpxbLWrW2ieAR4AzgZ7AX621j3szMBFp2DS7feNz6J3NZK/c7eswRJqUalvUjDFfW2vPMcZkAxYoXY3g98YYJ3AIeNJa+7KX4xSRBkbTZDQ++TuyfB2CSJNTbaJmrT2n5L8tKjtujIkEvgWUqIk0Ap6cjV7TZIiInDi3VyaojLU2wxiT6KFYRMSHvDGmTLPbi4icGHff+qyStXafJwIREd/SmDIREf9zwomaiDQOpWPKAg0aUyYi4idOqOtTRLzD4SzEYgkOCKm3a2pMmYiI/1GiJuKHvkhfQlZhJuM7XVev19WYMqnOj45C1v1yiOFpmfqciNQTdX2KiEiNUtIyufXoYV7ee4hJs5K02oRIPVGiJuLHtAST+IukHekUAk70solIfVLXp/i1Y0fySN99hA6ntiEk9MQ/rvuO7mPr4a0M6jCI0ED/Xq42dV8QD3yoJZjE9woP5nDm+sMEA4VAkDEM6tzK12GJNAlqURO/lpudz56t6TjyizxS39d7v+bmz28mK9//Z1jfsSdE02WIT1mHkyOf72T/c2s4M6OQWQmncXPXtjznDKfjO9vIWXcQa62vwxRp1NSiJuJBdZ3Zf3bacxTagnL7otqFExDQDOs0mADLL+GfMjutiCndbvd02CIV5O84TObCbTgO5hIeE0Wr0afSqWUIQ4GC3dlkvr+VQ29tIqxHa1pdehpBbcJ8HbJIo6RETcRDTmRm/+OTNIDO7XOZeEkaO/c2o2vHHDq3z6VQjRfiZUXHCsn6+GdyUvYT2DqUyOt6E96jTbkyIZ1bcPLN/Tn67V6OfJbK/mdTaHl+N5qf0xETqI4aEU9SoibiIZXN7H+iY8o6t8+lc/tcD0UoUr381Cwy/v0TzlwHLYZ1psV5XQkICay0rAk0tBjaifC+kRz+cDtZn/xMzg8HiJrah8Dm9Tf/n0hjp0StqTpyBAICoHnzcrttYSFpV19DqwlX0Grs2HoP67slmyhyOCvsX7t8e7ntwKAABl18Zn2F5ZbSmf0LHU7N7C8NUlDbZoR0aUnLkdGEdIhw75xWYURe3Yu8DRnkbsggICLYy1GKNC1K1DzNWsACBow54eqy9u0hKDSMiDYe/qX/6cfQqjWMGFl+v7Xk/vADzYcN8+z13FRZknYi5RLmJnDMcezX83K64shJJPGN6whsttO1PyIogqRJSbUL9jia2V8ausCIYKKu7V3r84wxhPeJIrxPlBeiEmnaGvVggr0r32L/8+ezP6Pi23LbfzhA6rp071z48/vh58/L7dq1YR0/Ll9a66q2fLmMvRv+56nImpzjk7ScnVMpOHgBOTunUpTTtdJyJyKuW2tuHn6akjQREfGIRp2oZe3fSbtD37Mr/WiFYz8s3cm6FbvrLZbNq1by9bw36+16UpEjpzvYQCAQbGDxtoiIiB9r1Ima1C9vzaK/JTuPD/ccZkt23gnVE9RsB5gioAhMUfG2jxXZIv6XtbpW5xQ6K74hKiIijVOjHqPmtJBP5W8fWWvrdaLG+r5ebWyOi8d5rHzX38GZMzk4c2a5fQEREfRISa60jhOZmqI6W7Lz+NtPv+BwWoICDPf1as8ZLeo2X1Ngs5006zoLR053gprtKDdGzRf25+9lZfpSMgoP1uq8d/b8i3Miz6dbs1O9FJmIiPiLRpeoFTqK+HJjBpn79tAppAtrYx4ne1MKT205RlzXNjTfk8uPy3eTm10IwJx7v6HfeZ3peXZHwjz8tlJhQT6bv13J94sWcGjPHsDy0vUTiTl/JLEjRtMish4H3s6eRUq2ISkvjISwPOLCCiAzE/7xCj2mXe8qtjrL8vHKDcSkb6dnZlq5Ko5P5sryxtQUABuP5OFwWpyAw2nZeCSvzokaFCdrvk7QALIdR/hw31s0C4xgRNtLWXrwQ7fPDQ4I4dMD73NJuwl0DO9a8wkiItJgNapE7cDM81lz2p/o+/PLtM/8FmdASMlbmBCdN4DPPrmNIhuAsb/2+B47nM/qj37m+8WpXHRTDJ171CG5WD4DivLL7cpIWcI7f32BQoelsMwbinlHs0lZ/CFrPl7E8Gt/T8x55d+6THrzNYoKC8vf15aNHNiysdy+7YURZPc4x+23C1OyDZN+OZkCawgxLZnb/kBxsla2TF4IV2eeTEHPzgQ7i/j7N/+okKxVxRNTU1TW4tizZRhBAcbVotazZeOY/bxFUEvOa3sJXcKjCQkIJTg9pNJJb48XbEK4vOPVbD+2iQ5hXeohUhER8aVGk6gdySlkc9drOH/N1QQ68zFYAkvG8uzO78PnmTfjJIjKJsxwFBQnUh8++wOBwQHc+EJi7S5+XJJ2JLuAeYu2k19Q+RQSRY7iROyz117gi9mvcsfchb8eOy5Jq8zW3FD+ticSR9pmt7sZk/LCKLAGJ4ZCW7x9fKLmKhMQQCGwLupUtxO1E52aIi+ngB3rfqmw/4wWYdzXqz0bj+TRs2VYuda07Wv30q1XO4IqmZAzz5HHa+teq1UM/039LyO6jcB4YFoVd5wa0cP1c22XhTqjee2nUBDxprounyYi1Ws0iVrKqm8Yun46Qc7yA86thc+zbsNBqFv1FBW6Nz+Xy8GNFXYtX7WPAjfrKXIUUliQT3BIKNbp3jk/5YbhsAYn7nczJoTlEWJaUmgh2FgSwioOzE8IyyOElhQ4nQQ7i4hJ315JTVWL69a6Tl/Qe7dlsHPTAag0jS5O1irr7ty/8zCHfsmme0wHIju2dO3/ds+3/F/S/7H7aO3e6p3+5XSGdhrK/Qn306l5p1qdK9KUeWuMqoh4+a1PY8xsY8wBY8z6Ko4bY8zzxphtxph1xpgBdblOXkERp29+iQBnfoVjewt6k2+bV3JW1Q7uzHavYFEBrJ9fbldOroOfd2dTm/cGNn+7EoAdSSvL7d+aG8qHh05ia275JLNXeB5BxhJoKNfNWN1bl3FhBcxtf4A7W2dV2u0JcMa+LTy2+nWu3vjfWnV7nqi8nAJatW1O//NOJTDIvY9kYFAA/YZ1JzQ8mMKCItd+ay0v/+9lggKCeH3E67WK456B95CyP4WPd3xcq/NEmrrKxqiKiGd4u0VtDvAiUNUEYqOA00v+DAJeKflvrezatYfu6SsIoGJ2tCVvKIXWvda0UtvWHGCndVTZjO9q4m++n7jjzt2xK5sAYygqiWVfaDv2hHWkU95eOuTvr/R6G1Ys44xBZ7N/86+tc8Xdm+1wWEOQsdzXaT+nhxcnoqeH53Nfp/0c653oiq/af9EWFCdlcWEFlSZoAM7CQgzQY98WetjNtXpeJ+qUPu0xAcWtabVdFqrvb04pt22M4elhT9M6rDUhgSFEBEW4NZltRFAEV/W6ivO7nU9kmJZ+EqkNLZ8m4j1eTdSstV8ZY6KrKXIp8KYtHkWeZIxpZYzpYK3dV5vrOPeupSgg1DUmrawcZytq23D44y9HeGbWlkqTnnIJkXEyd0Aoca1+vW5enoOiki7MfaHt+KD9JRSZQAJtEWN/+ajSZC03+wj7N/9Urufv1+5Ng8MWb5cmalCcrPXr25LmUcVxVfvW5RfLarzngOBgmnXuRNTgQaR/e2JLKdVWaZJWp3MrGU/WLqKd6+faLgvVPqJ9nWMRaaq0fJqI9/h6wttOwK4y27tL9lVgjJlmjEk2xiQfPFh+3qkw46CqvsZQc4w9gUUkhRayJ7Co0jLH21ZYUGUzfrmEyGlIyjwJgJTDLXjp585sd0YSUJI87AnrSJEJxJoAikwge8I6VrjWvtB2fBvSg1UbdmKLfo2vtHszAEuQsfQKrzim7Mgvv+azpf+iPb47lAP7Ic296SgCgoNp0z8WE1hxcH6pQ/+ZWy5OERHQ8mki3uLrlwkqa0qpNOOy1r4GvAYQHx9frkxkh04UD62vKCdsN++E5lOEIRC44mgInYqKE5E9gUXsCnLSxRHg2gcwPKY9iw4cqrQZv1wTv7EktM4i5XALJq3pQ4EzgGDThbEhGZyc+wud8vYSaIsoAgJtEZ3y9paLrbTFzWmDWLEe7usUWqF786fcMHqF55VrTStVkJfr+rnCv2i7toIN6+Hbb6p6pFVqcdqpHNm8pdJj+x95hCOLF9Pt329igj0775yIiIiU5+tEbTdQdjKozsDeKspWzlqK9v2IIyCM4KKcCod3BlqKHGANFFnYFeSkU1EgewKLeKd5QXESRfkELrbzSVU248d1a83c6weR9L/1JDi+Je6kbF76uTMFzoCSbkpDZpuunLznFzrk72fsLx9VOUbN1eKGwWFtpd2blSVopfat/x9hLVrS7oyeGGPKv3W5dy98vRJCQyG/6jqOFxAcTEibNlUeb/fQgzgPH1aSJiIiUg983fW5CLi65O3PBCCrNuPTihwO0l+6kJOW3UVWWFccARWncBgcuIFgijC2OCHr4ii+5V1BToooSeBKtku9/9Qajn5zgKkJ0ZU248cFbOLmiE+I6xgGAUEktM4iJMBJIE6CAyyXn1lIUGBxY2GH/P3EZ/1Q6di04hY3JwFYAo2lV/PareHYLDKStO+/xZGfX/GNz44d4ZJLIap22lFqpQAAGJtJREFUqx8YY6ocMxbQvDnB7doT9Yc/1KpOERERqRuvtqgZY+YBiUCUMWY38BAQDGCtfRX4GLgI2AbkANfVpn5rAtnb8ixWNzuftudMpOU754B1ElRmhve4gK3MC/0/3ssdR2BOX9oXFSchXRwBBFLcylY2gQPod24XNif9QvzF0YSGV/KIOsZBYAh06A/J/yDOuYu5A9aTlHkSCa2ziGvloLvpwscrdlFUZKucqqO4xW0RIYPHwA9L6XFmDE6HmzdvDDGjLycvO4t1+3Mrf+OzY0f4uU1x65qb84U4CwoozDpS+UFrsQXut86JiIjIiTH+ulB4deLj421ycvHi4EdyCwkPCSQ4MIAFK7Zyzv9uITJrLYHOAgLKjFsrNGH8kn8anx++lWPOk4HKx6gFhwYy7blh5Oc6Kk/SjldwDL5/BfIOgy0/yD4jM4/V6w6yeXsW1loqm882JCycW994l/ycHHIOHWTTsk9xFlWXrRmCQkOJ/e0EQiOK54d7afk2nl66GaeFQAN3jujBzcNPKy6emQnvvQtuvgDgLCxk62uv4yyo2LoXEBFBl3+8SrP4eLfqEhERETDGpFhr6/TL09dj1E5Yy/Bfx0qNSzwdEv8Le9dS9M0LFP68EorysWFtCD7rOrr0v4prm7Uh71gh2YfyMAZatAkjtFnF8VZuJWkAIREw6FbYsQz2rAYMWCcYQ2RkEKMuOZ3zOg/jcEE4RY5CmreJpEWbit2Roc2aEdqsGz3OHcHm5UsBcDrKJ2wBwcEEh4TS+6KxriQNapjDqHVriIyCgwcqtKql5IWUW6Td6XCQtWlzpUkagAkJIbx/f/eeiwcUOZwcy8qjyOEkKCSQ5ieF1Xkqj6MFR1m0fRFf7f6KHEcOUeFRjD1tLEM6DiEwoOq3XEXEc7TMlEjtNfgWNb/idEDGVijIhoAgaNkFItrWuhpHQT4Htm5i34Z1FOTmYEwAzVq3oVPf/rTpGo0JqDi0sNovwGPHilvV8vJcyVpKXkiZRdot/2m7jz65+0md9w7WUbFFz4SGEnXjjUT94Ub3r1tH+bmF7NmazoGdh13zpFkgwBg6dG9Dh1PbEBTsXnJlreWltS8xZ8McAgggt+jXN2WbBTX7//buPcjuurzj+Oc5Z3ezm2wumwsQ2FwQiSYQAs0iAYdbRQqtBS1YwXFKrMo4I8UZL61tHUfsdLTFsWOnaadCrRaxoLRSrKlgC61VAdmgAQlSYyAkBMxusglJdje7Z8/TP/Ys2cvZPb9z9vzO+Z7f7/2ayWQvv3PyJN/89jzn+32+z1etTa367MWf1UWnXlSV2AEUxzFTSLNUz6gFJdMkLVs766dpapmjU8/aoFPP2hD5MTOeszlvnnTdO6XvPST19kj5/JRD2r//8qDav/vN4klaS4taTl+txe/dPOHrcfzgPXZ4UD/74QsayeULXUVOvJHIS9r7i17t33NI6y9erZbWmXeeurtue/Q2fef57+j4yNTauv5cv/pz/frwwx/W5y/9vC5dcemsYgcwvRmbcgOYVr13faJW5s2T3v6O0YTt7PXatHKhWjKjveWaMtLqH39PapqUt2cysrY2tW3YoFVf+5oyrRN31Vb7fL+hwdxokjacn7b1m+ddxweG9bMf7FY+P/Ns8MN7HtbW57dqMDe1WfB4gyOD+vj3P67Dxw9XGjqAEqZtyg1gRsyopU1Hh3TRm7XxIunui04sW577x5fq0Lfu16F779VIX5+suVlzuzZq8eb3qm392UWfqtrn+72864DyIxGW4l0aGhzWwVeOaOmpC6a97I6n7tBAbmDa7094Snfdv/N+3XTWTVHDBVAGjpkCKkONGmalWjVq+bzr7m9s188O9mvtglatmT+1J95k8xa2asNlryv6vZeOvqRr77+26JLndJbPW66Hrn8o8vUAAERBjRrqZsbauDI89tx+/dlT+5TLu5oypj9Zd0rJZO3Y4UG5e9GD2fcd3afmTHNZiVrvQG/ZcQMAECcSNQTh8RcOKpd35SXl8q5nXx2MNKvmLhXJ05S18ltuZKz+JZv7j7+spw93q3d4v9xdC5s7tH7BRp3WuqpoQgoASDYSNQThwjOW6u/+9/nXZtTWLiidpFnGlJmmr9qqBas0NFLekVwrF6ws6/pqOpY7qu/u/xcdGu7TiOfkhd0Uh3N92je4R62ZNl198nVa3FLekWAAgMZW/ykEQNKmNct0W9dKvXNlR6RlT5lm3EiwpG2JNi3fJFO0Wai5TXP13rPKOsGsagZG+vWvL9+lg0O9yvnwa0namJwP6+jIq7r/5bvVNzS7nbUAgMZCooZgXHF+p96xsiPSkmcmY1p+xuIZr7l5w82ak50T6c9ubWrVlauvjHRttX3/wEMaHOlXXkXOGBtn2If0vZ5/q1FUAIAQkKghGEtPW6j2jraSx0RlsqaTVi5S+6K2Ga/bsGyDPtb1sRmTtYwyam9u1x1X3hE5qaumgZFj2tO/q2SSNuZI7rB6jr8Sc1QAgFCQqCEYljGt3bRSC5fOUyZbJFmz0WtOXtWh09efEuk53/XGd+n2S25XZ3un2pralCn8l2/JtKgl06LzTzlf977tXq3pWFPNv0pkvzz2XPHdENMY8RHtOLI9xogAACFhMwGCks1mtO7ClTrSN6B9O3t1qOeY8iOubFNGS09boOWvW6K29paynvPylZfrshWXaXvPdnX/qlv9w/1a3LpYV6y6QqfMi5bwxeVo7lWN+NRju6bjch3JHYoxIgBASEjUEKT5HW16w/krqvZ8T754SI/tatem110XVEf0rJV/C1byGABAY+InPhIvjsPjq2XZnJPVbM0a9uFI12etSafMOS3mqICwVOsEFKARkagh8YodHh/KD/uVbWcoY1kpYqImd71x/jnxBgUEJOQ3WkAtsJkAiTd2eHzWVJXD46spYxn92sIL1RRhObPJmvT69nVqy86tQWRAGIq90QLShBk1JN7GVR26+/2bgl06Wb9go/qGe7Xz2M+Vm2ZmrcmadFLLcl285K01jg6or7E3WsO5fHBvtIBaMHcvfVVgurq6vLu7u95hAFXj7vr50af15KFHNZgfGD2dwEdn3DKW0TkLunTuwguCOI8UqDVq1NDozGybu3dV8lhm1IAAmJnWzj9Hb2xfr1eOv6S+4QNyz2t+0yJ1tq0iQUOqbVzVQYKG1CJRA6qgWu/4zUzLWzu1vLWzitEBABoViRowS+xKA8LBMimShkQNmKWQ238AacKbJiQRhS/ALIXc/gNIE1p5IImYUQNmKfT2H0Ba0MoDSUR7DgBAYlCjhhDRngNIIF5wgPLRygNJQ6IGBIiiaACAxGYCIEgURQO1s213n7Y8slPbdvfVOxRgCmbUgABRFA3UBrPXCB2JGhAgdpICtUEfRISORA0IFEXRQPyYvUboSNQAAKnF7DVCR6IGpBCtP4ATmL1GyEjUgJSheBoAGgftOYCUofUHMHu09ECtMKMGpAzF08DsMCuNWiJRA1KG4mlgdmjpgVoiUQNSiOJpoHLMSqOWSNQAACgDs9KoJRI1AADKxKw0aoVdnwAAAIEiUQMAAAgUiRoAAECgSNQAzAqNPwEgPmwmAFAxGn8C0XC+LipFogagYjT+BErjDQ1mg6VPABUba/yZNdH4E5gG5+tiNphRA1AxGn8CpXGSAWbD3L3eMZStq6vLu7u76x0GAACRUKOWbma2zd27Knls7DNqZnaVpC9Kykq6090/N+n7KyV9VdKiwjWfcPetcccFAECtcJIBKhVrjZqZZSVtkXS1pHWSbjSzdZMu+6Skb7j7eZJukPS3ccYEIGy0+wBGcS9Ain9G7U2Sdrr7Lkkys3skXStpx7hrXNKCwscLJe2LOSYAgWJ3HDCKewFj4t71eZqkPeM+31v42niflvQeM9sraaukPyj2RGZ2s5l1m1l3T09PHLECqDN2xwGjuBcwJu5EzYp8bfLuhRslfcXdOyX9pqS7zGxKXO7+JXfvcveuZcuWxRAqgHqj3QcwinsBY+Je+twracW4zzs1dWnzfZKukiR3f9TMWiUtlbQ/5tgABKaSdh/spkMS0foGY+JO1J6QdKaZnS7pJY1uFnj3pGtelPQWSV8xs7WSWiWxtgmkVDm746jjQZKxUxRSzEuf7p6TdIukByU9q9Hdnc+Y2WfM7JrCZR+V9AEz2y7pnyVt9kZs7gag5qjjAZB0sfdRK/RE2zrpa58a9/EOSW+OOw4AyUPHd2B6lAUkA0dIAWhY1PEAxVEWkBwkagAaGnU8wFTFygK4TxpT3O05AKBh0AkeSUF7j+RgRg0AxFIRkoWygOQgUQMAsVSE5KlGWQAbEuqPRA0AxA5SYDJmmcNAogYAYqkImIxZ5jCQqAFAATtIgROYZQ4DiRoAAJiCWeYwkKgBQAwowkYSMMtcfyRqAFBlFGEjrXiDUn0kagBQZRRhI414gxIPTiYAgCqjKzzSqNgbFMweM2oAUGUUYSON2CUaD3P3esdQtq6uLu/u7q53GAAAYJxq1Kglsc7NzLa5e1clj2VGDQAAVMVsd4lS5zYVNWoAACAI1LlNRaIGAACCwEacqVj6BIBAJbFWB5hJpRtxknyvkKgBQICo1UFalVvnlvR7hUQNAAIUUtPcoX1Hdezxl5XrHZAyppaVC9R+wSnKLphTl3iA8UK6V+JAogYAAQqhJ1XuwIAOfO1Z5XoH5Lm8VOjmdPz5wzryP3vUtm6JOq5fo0xLtuaxAWNCuFfiRB81AAhUPetucgcG9Ku/+al8MPdagjZFk6n55Hk66YMbZM3sTUP9hF6jRh81AEig2fakqpS7q/efdsycpElSzjX8q2M6/ODzWvS2M2oWHzDZbO6V0JM8EjUASIFyXoyG9x7VyMHBmZO0MTnXsR+/ooW/sVrWzBIoGksjbEQgUQOAhCv3xejoj/aN1qRFZhrYcVBzNyybfbBADTXCRgSKCgAg4crt9j7c0x9tNq3Ah0eU6xucZZRA7TVCg11m1AAg4crdFWcZK+8PMMmszMcUhF4fhGSrtMHumFr8/yVRA4CEK/fFaMeCrB614zrPszo7wsuENWXVdPLcsuOqdn0QSR8qUelGhFrVt5GoAUAKRH0x2ra7Tx/8+V4NeV7Nkr6ouSWTNWsyta4p/wWqmvVBjVAUjmSpVX0bNWoAgNc8tuuAhkbyyksalvQTjcx4vTVnNP/SzvKXS1Xd+qBy6/CA2apVfRszagCA14yvZ2ty6bxsk6bL1aw5o9a1i9V+cWdFf9Zs64OmizvUonAkSzX//86EkwkAABOM1Xpd0LlIpz/Zq4Gne2Vm8uHRlh02JyuZNP+STs2/bMWE2bR61olV48+mzg1xmM3JBCRqAIAZ5fuH1f90r3IHB2VZU/PydrWtWyzLTqyeafQ6sUaPH+HiCCkAQGwyc5vVfsHyktc1QvPQmTR6/EgmNhMAAKqiEZqHzqSa8W/b3actj+zUtt19VYwQacTSJwCgahq9xqtadW4soWI8lj4BAEGotHloKKoRP0uoqCaWPgEAQWn0ZcNGXwJGWJhRAwAEIwnLhrXqr4V0IFEDAAQjKcuGjb4EjHCw9AkACAbLhsBEzKgBAILBsiEwEYkaACAoLBtO1OgtTzA7JGoAAAQqrs0VJH+Ng0QNAIBAxbG5Igk7a9OEzQQAAAQqjs0VxZI/hIsZNQAAAhXH5oqx5G84lw92Zy1Lsydw1icAINXSmBRU++9czedL4tIsZ30CAFCBJCYFUVRzZ221/w2T0vS4WqhRAwCkFvVas1ftf0OaHk/EjBoAILUaoV4rdNX+N6Tp8USx16iZ2VWSvigpK+lOd/9ckWt+V9KnJbmk7e7+7pmekxo1AEC1pLFGrdr4N5zZbGrUYk3UzCwr6f8kvVXSXklPSLrR3XeMu+ZMSd+Q9Ovu3mdmJ7n7/pmel0QNAABEEUISGfJmgjdJ2unuuyTJzO6RdK2kHeOu+YCkLe7eJ0mlkjQAAIAokrBZJO7NBKdJ2jPu872Fr423RtIaM/uhmT1WWCqdwsxuNrNuM+vu6emJKVwAAJAUSdgsEneiZkW+NnmttUnSmZIuk3SjpDvNbNGUB7l/yd273L1r2bJlVQ8UAAAkSxw7SLft7tOWR3Zq2+6+KkRYWtxLn3slrRj3eaekfUWueczdhyU9b2bPaTRxeyLm2AAAQIJVewdpPZZS455Re0LSmWZ2upm1SLpB0gOTrrlf0uWSZGZLNboUuivmuAAACEqtZ2rSYuOqDn3o8tdXJaGqx1JqrDNq7p4zs1skPajR9hxfdvdnzOwzkrrd/YHC9640sx2SRiR93N0bbxEZAIAKJaHoPQ3q0Xcv9oa37r5V0tZJX/vUuI9d0kcKvwAASB2OTWoM9WjGy8kEAADUGSckNI5qnpMaBYkaAAB1xrFJmA6JGgAAAaj1TA0aQ9y7PgEAAFAhEjUAAIBAkagBAAAEikQNAAAgUCRqAAAAgSJRAwAgpTi2Kny05wAAIIU4tqoxMKMGAEAK1eOAcZSPRA0AgBQaO7Yqa+LYqoCx9AkAQApxbFVjIFEDACClOLYqfCx9AgAABIpEDQAAIFAkagAAAIEiUQMAAAgUiRoAAECgSNQAAAACRaIGAAAQKBI1AABQFxwKXxoNbwEAQM1xKHw0zKgBAICa41D4aEjUAABAzXEofDTm7vWOoWxm1iNpd8TLl0rqjTEcVB9j1pgYt8bDmDWmxIybtbTNy7TMnZ8f6j/iQwPH6h1PjN7g7vMreWBD1qi5+7Ko15pZt7t3xRkPqosxa0yMW+NhzBoT49Z4zKy70sey9AkAABAoEjUAAIBApSFR+1K9A0DZGLPGxLg1HsasMTFujafiMWvIzQQAAABpkIYZNQAAgIZEogYAABCoxCRqZnaVmT1nZjvN7BNFvj/HzO4tfP9xM1td+ygxXoQxu8TMnjSznJldX48YMVWEcfuIme0ws6fM7L/MbFU94sQJEcbsg2b2tJn91Mx+YGbr6hEnTig1ZuOuu97M3Mxo1xGACPfaZjPrKdxrPzWz95d6zkQkamaWlbRF0tWS1km6scgPmvdJ6nP310v6K0l/UdsoMV7EMXtR0mZJX69tdJhOxHH7iaQudz9H0n2S/rK2UWK8iGP2dXdf7+7nanS8vlDjMDFOxDGTmc2XdKukx2sbIYqJOm6S7nX3cwu/7iz1vIlI1CS9SdJOd9/l7kOS7pF07aRrrpX01cLH90l6i5lZDWPERCXHzN1fcPenJOXrESCKijJuj7h7f+HTxyR11jhGTBRlzF4d9+k8Sewyq68or2mS9GcaTawHaxkcphV13MqSlETtNEl7xn2+t/C1ote4e07SYUkcLFY/UcYM4Sl33N4n6T9ijQilRBozM/uQmf1Soy/8t9YoNhRXcszM7DxJK9z932sZGGYU9efjdYXSkPvMbEWpJ01KolZsZmzyO8Io16B2GI/GFHnczOw9krok3R5rRCgl0pi5+xZ3P0PSH0n6ZOxRYSYzjpmZZTRawvPRmkWEKKLca9+WtLpQGvKfOrHSN62kJGp7JY3PSjsl7ZvuGjNrkrRQ0sGaRIdioowZwhNp3MzsCkl/Kukadz9eo9hQXLn32j2S3h5rRCil1JjNl3S2pP82sxckbZL0ABsK6q7kvebuB8b9TLxD0sZST5qURO0JSWea2elm1iLpBkkPTLrmAUk3FT6+XtLDTrffeooyZghPyXErLMn8vUaTtP11iBETRRmzM8d9+luSflHD+DDVjGPm7ofdfam7r3b31RqtBb3G3Ss++BtVEeVeWz7u02skPVvqSZuqGmKduHvOzG6R9KCkrKQvu/szZvYZSd3u/oCkf5B0l5nt1OhM2g31ixhRxszMzpf0LUkdkn7bzG5z97PqGHbqRbzXbpfULumbhf06L7r7NXULOuUijtkthVnQYUl9OvGmFnUQccwQmIjjdquZXSMpp9FcZHOp5+UIKQAAgEAlZekTAAAgcUjUAAAAAkWiBgAAECgSNQAAgECRqAEAAASKRA1AKpnZj+odAwCUQnsOAACAQDGjBiCVzOxo4fc/NLOnzWy7mX3OzJrM7Akzu6zw/c+a2Z/XNVgAqZWIkwkAoBJmdrVGz7W8wN37zWxxobv4Zkn3mdmtkq6SdEE94wSQXiRqANLsCkn/6O79kuTuBwu/P2Nmd0n6tqQL3X2ojjECSDGWPgGkmUmarlB3vaRDkk6uXTgAMBGJGoA0e0jS75vZXEkys8WF339H0hJJl0j6azNbVL8QAaQZuz4BpJKZHXH3+Wb2CUm/J2lI0lZJX5D0I0lvcfc9hTq1je5+Ux3DBZBSJGoAUsfMlkh60t1X1TsWAJgJS58AUsXMTpX0qKTP1zsWACiFGTUAAIBAMaMGAAAQKBI1AACAQJGoAQAABIpEDQAAIFAkagAAAIH6fykf2Jt1x6/oAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 720x504 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# Read curves taken with Data Theif\n",
"curves = pd.read_csv('Jacobi_Ellipsoid_Data.txt')\n",
"\n",
"def jacobi_x(a, b, c, omega, rho):\n",
" t1 = (np.sqrt(3)/10)\n",
" t2 = (a**2 + b**2) / (a*b*c)**(2/3)\n",
" t3 = np.sqrt(omega**2/(np.pi*G*rho))\n",
" return t1*t2*t3\n",
"\n",
"def jacobi_y(a, b, c, ax='a'):\n",
" d = {'a': a, 'b': b, 'c': c}\n",
" num = d[ax]\n",
" denom = (a*b*c)**(1/3)\n",
" return num/denom\n",
"\n",
"# Compute the jacobi x and a, b, c values\n",
"df['jcx'] = jacobi_x(df['a'], df['b'], df['c'], df['omega'], df['rho'])\n",
"df['jcy_a'] = jacobi_y(df['a'], df['b'], df['c'], 'a')\n",
"df['jcy_b'] = jacobi_y(df['a'], df['b'], df['c'], 'b')\n",
"df['jcy_c'] = jacobi_y(df['a'], df['b'], df['c'], 'c')\n",
"\n",
"# Jacobi Plot\n",
"f, ax = plt.subplots(figsize=(10, 7))\n",
"ax.set_xlim(-0.01, 0.5)\n",
"ax.set_ylim(0.5, 1.75)\n",
"ax.set_title('Jacobi Ellipsoids')\n",
"ax.set_ylabel('Semi-Principal Axis Lengths')\n",
"ax.set_xlabel('Normalized Angular Momentum')\n",
"\n",
"# Plot Jacobi a, b, c for each object\n",
"df.plot('jcx', 'jcy_a', 'scatter', c=colors, marker='s', s=100, \n",
" facecolors='None', ax=ax, \n",
" colorbar=False, legend=True)\n",
"df.plot('jcx', 'jcy_b', 'scatter', c=colors, marker='2', s=400, ax=ax, \n",
" colorbar=False)\n",
"df.plot('jcx', 'jcy_c', 'scatter', c=colors, marker='o', s=150, ax=ax,\n",
" colorbar=False)\n",
"\n",
"# Plot curves\n",
"ax.scatter(curves['x'], curves['y'], marker='.')\n",
"\n",
"# Legend\n",
"new_legend = [Line2D([0], [0], marker='s', color='k', label='a',\n",
" markerfacecolor='k', markersize=10),\n",
" Line2D([0], [0], marker='2', color='k', label='b',\n",
" markerfacecolor='k', markersize=12),\n",
" Line2D([0], [0], marker='o', color='k', label='c',\n",
" markerfacecolor='k', markersize=10)]\n",
"new_legend.extend(legend_entries)\n",
"plt.legend(new_legend, np.hstack(['a', 'b', 'c', df['name'].values]), ncol=4)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The 8 bodies all plot very close to the predicted curves in both the Maclaurin and Jacobi cases. In the Jacobi plot, Jupiter and Saturn fall off the curve, but this may be because I used DataTheif to determine the standard curves. Haumea is interesting because it is in the Jacobi regime and its shape is well predicted by the Jacobi solution, but also appears to fit on the Maclaurin plot. This may be because Haumea is near the point of divergence of the two solutions, so Maclaurin is not yet a poor assumption."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Q2) Shape vs J2\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"def J2toW(J2, a, M):\n",
" return np.sqrt(2*G*M*J2/a**3)\n",
"\n",
"def J2toW_flat(J2, a, M, f):\n",
" t1 = f - 1.5*J2\n",
" t2 = 2*G*M/a**3\n",
" return np.sqrt(t2*t1)\n",
"\n",
"orbit_conversion = (60*60*24) / (2*np.pi) # rad/s -> rotations/day"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### a) Moon"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The rotation period of the moon should be 3.724 days\n"
]
}
],
"source": [
"moon_J2 = 2.037e-4\n",
"moon_r = 1.7371e6 # m\n",
"moon_M = 7.35e22 # kg\n",
"moon_w = J2toW(moon_J2, moon_r, moon_M)\n",
"\n",
"print('The rotation period of the moon should be {:.3f} days'.format(1/(moon_w*orbit_conversion)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The rotation period is smaller than the observed 28 day rotation period of the current lunar day. This implies that the Moon was spinning more rapidly in the past, potentially while still molten, freezing in the shape of a more rapidly rotating body as it cooled, and maintaining that shape as it slowed in rotation over tha last several billion years."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### b) Mars"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The rotation period of mars is 1.280 days\n"
]
}
],
"source": [
"df = df.set_index('name')\n",
"mars_f = 3.352811e-3\n",
"mars_J2 = 1.247e-3\n",
"mars_r = df.loc['Mars', 'a'] # m\n",
"mars_M = 6.39e23 # kg\n",
"\n",
"mars_w = J2toW_flat(mars_J2, mars_r, mars_M, mars_f)\n",
"\n",
"print('The rotation period of mars is {:.3f} days'.format(1/(mars_w*orbit_conversion)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The roation period is larger than the current 1.03 day rotation period of Mars. This implies that Mars was spinning more slowly in the past and has sped up since it cooled, potentially due to impacts imparting more angular momentum to the planet."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Q3) Mercury\n",
"\n",
"### a) Density and gravitational acceleration"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The bulk density of Mercury is: 5423.203 kg/m^3\n",
"The acceleration due to gravity on Mercury is: 3.699 m/s^2\n"
]
}
],
"source": [
"def volume_sphere(r):\n",
" return (4/3) * np.pi * r**3\n",
"\n",
"def g(M, R):\n",
" return G * M / R**2\n",
"\n",
"R = 2.44e6 # m\n",
"M = 3.3e23 # kg\n",
"rho = M / volume_sphere(R)\n",
"g = g(M, R)\n",
"\n",
"print('The bulk density of Mercury is: {:.3f} kg/m^3'.format(rho))\n",
"print('The acceleration due to gravity on Mercury is: {:.3f} m/s^2'.format(g))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### b) Core thickness\n",
"\n",
"Knowing $R = 2.44 \\times 10^6$, $\\rho = 5423$, $\\rho_c = 7250$ and $\\rho_m = 3500$, compute the mantle thickness.\n",
"\n",
"Assuming the crust is negligible, $M = M_c + M_M$.\n",
"\n",
"Rewriting in terms of density: $\\rho V = \\rho_c V_c + \\rho_m (V - V_c)$\n",
"\n",
"Let $\\eta = \\frac{4 \\pi}{3}$. Then isolating $R_c$, we can write:\n",
"\n",
"$\\rho \\eta R^3 = \\rho_c \\eta R_c^3 + \\rho_m \\eta (R^3 - R_c^3)$\n",
"\n",
"Isolating $R_c$ and cancelling $\\eta$:\n",
"\n",
"$R_c^3 (\\rho_m - \\rho_c) = R^3 (\\rho_m - \\rho$)\n",
"\n",
"And since $R = R_m + R_c$, we can write:\n",
"\n",
"$R - R_m = R \\left(\\frac{(\\rho_m - \\rho)}{(\\rho_m - \\rho_c)}\\right)^{1/3}$\n",
"\n",
"Finally,\n",
"\n",
"$R_m = R \\left(1 - \\left(\\frac{(\\rho_m - \\rho)}{(\\rho_m - \\rho_c)}\\right)^{1/3}\\right)$"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The radius of the mantle is about 4.87e+02 km\n"
]
}
],
"source": [
"rho_c = 7250 # kg/m^3\n",
"rho_m = 3500 # kg/m^3\n",
"\n",
"R_m = R * (1 - ((rho_m - rho)/(rho_m - rho_c))**(1/3))\n",
"print('The radius of the mantle is about {:.3} km'.format(R_m/1000))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### c) Normalized moment of inertia"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The normalized moment of inertia of Mercury is about 0.349\n"
]
}
],
"source": [
"R_c = R - R_m\n",
"\n",
"def get_C(rho_c, R_c, rho_m, R):\n",
" t1 = 8 * np.pi / 15\n",
" t2 = rho_c * R_c**5\n",
" t3 = rho_m * (R**5 - R_c**5)\n",
" return t1 * (t2 + t3)\n",
"\n",
"def get_M(rho_c, R_c, rho_m, R):\n",
" t1 = 4 * np.pi / 3\n",
" t2 = rho_c * R_c**3\n",
" t3 = rho_m * (R**3 - R_c**3)\n",
" return t1 * (t2 + t3)\n",
"\n",
"def norm_MOI(C, M, R):\n",
" return C / (M * R**2)\n",
"\n",
"norm_moi = norm_MOI(get_C(rho_c, R_c, rho_m, R),\n",
" get_M(rho_c, R_c, rho_m, R), R)\n",
"\n",
"print('The normalized moment of inertia of Mercury is about {:.3}'.format(norm_moi))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### d) Comparison\n",
"\n",
"The normalized moment of inertia we calculated for Mercury is close, but outside the confidence levels derived from spacecraft data (0.349 is about twice the confidence level outside 0.353 +/- 0.017). A possible cause of this is that the mantle is not homogeneous in density, causing us to slightly overestimate the core radius which would cause the normalized moment of inertia calculation to report a lower value."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 4) Backtracking to Gaspra\n",
"\n",
"Using the M, e, and a 50 days after the fly-by, and knowing the mean motion we will backtrack to the mean anomaly during the fly-by and use the 3rd order approximation of the eccentricity anomaly to solve for heliocentric distance."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The heliocentric distance during the Gaspra fly-by was 2.222 au\n"
]
}
],
"source": [
"def get_r(a, e, E):\n",
" \"\"\"Compute heliocentric distance.\"\"\"\n",
" return a * (1 - e * np.cos(E))\n",
"\n",
"\n",
"def E3(M, e):\n",
" \"\"\"\n",
" The 3rd order approximation of the eccentric anomaly from Murray and \n",
" Durmott).\n",
" \"\"\"\n",
" t1 = M\n",
" t2 = (e - (e**3 / 8)) * np.sin(M)\n",
" t3 = (e**2 / 2) * np.sin(2 * M)\n",
" t4 = (3 / 8) * e**3 * np.sin(3 * M)\n",
" return t1 + t2 + t3 + t4\n",
" \n",
"\n",
"M0 = 293.23 # degrees\n",
"e = 0.174\n",
"a = 2.21 # au\n",
"M_motion = 0.302 # degree/day\n",
"\n",
"# Back out mean motion, compute E and use it to compute r\n",
"M = M0 - 50 * M_motion\n",
"E = E3(M*(np.pi/180), e)\n",
"r = get_r(a, e, E)\n",
"\n",
"print('The heliocentric distance during the Gaspra fly-by was {:.3f} au'.format(r))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 5) Tisserand's parameter"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>a</th>\n",
" <th>e</th>\n",
" <th>i_deg</th>\n",
" <th>i</th>\n",
" <th>tiss</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>mystery</th>\n",
" <td>4.210</td>\n",
" <td>0.580</td>\n",
" <td>5.50</td>\n",
" <td>0.095993</td>\n",
" <td>2.694795</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Object 1</th>\n",
" <td>10.381</td>\n",
" <td>0.649</td>\n",
" <td>2.38</td>\n",
" <td>0.041539</td>\n",
" <td>2.648444</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Object 2</th>\n",
" <td>5.325</td>\n",
" <td>0.205</td>\n",
" <td>11.13</td>\n",
" <td>0.194255</td>\n",
" <td>2.920185</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Object 3</th>\n",
" <td>8.650</td>\n",
" <td>0.783</td>\n",
" <td>19.49</td>\n",
" <td>0.340165</td>\n",
" <td>2.113595</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Object 4</th>\n",
" <td>9.238</td>\n",
" <td>0.495</td>\n",
" <td>7.89</td>\n",
" <td>0.137706</td>\n",
" <td>2.856712</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" a e i_deg i tiss\n",
"mystery 4.210 0.580 5.50 0.095993 2.694795\n",
"Object 1 10.381 0.649 2.38 0.041539 2.648444\n",
"Object 2 5.325 0.205 11.13 0.194255 2.920185\n",
"Object 3 8.650 0.783 19.49 0.340165 2.113595\n",
"Object 4 9.238 0.495 7.89 0.137706 2.856712"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def tisserand(a, e, i, ap):\n",
" t1 = (ap / a)\n",
" t2 = 2 * np.cos(i)\n",
" t3 = np.sqrt(a*(1 - e**2)/ap)\n",
" return t1 + (t2 * t3)\n",
" \n",
"d = {\n",
" 'mystery': (4.21, 0.58, 5.5),\n",
" 'Object 1': (10.381, 0.649, 2.38),\n",
" 'Object 2': (5.325, 0.205, 11.13),\n",
" 'Object 3': (8.65, 0.783, 19.49),\n",
" 'Object 4': (9.238, 0.495, 7.89)\n",
"}\n",
"\n",
"df = pd.DataFrame.from_dict(d, orient='index', columns=('a', 'e', 'i_deg'))\n",
"df['i'] = df['i_deg'] * (np.pi/180)\n",
"df['tiss'] = tisserand(df['a'], df['e'], df['i'], 5.2044)\n",
"df.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"According to the Tisserand's parameter computed above, since it is between 2 and 3, it is likely a Jupiter family comet. Although it has a similar Tisserand's parameter as Objects 1 and 4, the 2-3 degree difference in inclination and 5-6 au difference in semi-major axis in both cases seem too large to be the same object. Conclusions about whether they are the same would depend on the precision with which the orbital elements were determined in both cases."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 6) Three-body problem\n",
"\n",
"### a) Solve equations of motion\n",
"\n",
"Compute some partials to help later:\n",
"\n",
"$\\frac{\\partial}{\\partial x} \\frac{1}{r_1} = \\frac{\\partial}{\\partial x} \\frac{1}{((x + \\mu)^2 + y^2)^{1/2}} = -\\frac{\\mu + x}{((\\mu + x)^2 + y^2)^{3/2}}$\n",
"\n",
"$\\frac{\\partial}{\\partial x} \\frac{1}{r_2} = \\frac{\\partial}{\\partial x} \\frac{1}{((x + \\mu - 1)^2 + y^2)^{1/2}} = -\\frac{\\mu + x - 1}{((\\mu + x - 1)^2 + y^2)^{3/2}}$\n",
"\n",
"$\\frac{\\partial}{\\partial y} \\frac{1}{r_1} = \\frac{\\partial}{\\partial y} \\frac{1}{((x + \\mu)^2 + y^2)^{1/2}} = -\\frac{y}{((\\mu + x)^2 + y^2)^{3/2}}$\n",
"\n",
"$\\frac{\\partial}{\\partial y} \\frac{1}{r_2} = \\frac{\\partial}{\\partial y} \\frac{1}{((x + \\mu - 1)^2 + y^2)^{1/2}} = -\\frac{y}{((\\mu + x - 1)^2 + y^2)^{3/2}}$\n",
"\n",
"#### Compute $U_x$:\n",
"\n",
"$\\frac{\\partial U}{\\partial x} = \\frac{\\partial}{\\partial x} \\left(\\frac{1}{2}(x^2 + y^2) + \\frac{1-\\mu}{r_1} + \\frac{\\mu}{r_2}\\right) = \\frac{\\partial}{\\partial x} \\frac{x^2}{2} + \\frac{\\partial}{\\partial x}\\frac{1-\\mu}{r_1} + \\frac{\\partial}{\\partial x}\\frac{\\mu}{r_2}$\n",
"\n",
"By chain rule with $u = \\frac{1}{r_1}$ and $v = \\frac{1}{r_2}$:\n",
"\n",
"$U_x = x + \\frac{\\partial}{\\partial u}(1-\\mu)u \\cdot \\frac{\\partial}{\\partial x}\\frac{1}{r_1} + \\frac{\\partial}{\\partial v}(\\mu v) \\cdot \\frac{\\partial}{\\partial x}\\frac{1}{r_2}$\n",
"\n",
"$U_x = x - \\frac{(1-\\mu)(\\mu + x)}{((\\mu + x)^2 + y^2)^{3/2}} - \\frac{\\mu(\\mu + x - 1)}{((\\mu + x - 1)^2 + y^2)^{3/2}}$\n",
"\n",
"$U_x = x + \\frac{\\mu^2 - \\mu + \\mu x - x}{((\\mu + x)^2 + y^2)^{3/2}} - \\frac{\\mu^2 + \\mu x - \\mu}{((\\mu + x - 1)^2 + y^2)^{3/2}}$\n",
"\n",
"Plug in $x_0=\\frac{1}{2}-\\mu$ and $y_0=\\sqrt{3}/2$\n",
"\n",
"$ U_x = \\frac{1}{2} + \\frac{\\mu^2 - \\mu + \\mu/2 - \\mu^2 - 1/2 + \\mu}{(1/4 + 3/4)^{3/2}} - \\frac{\\mu^2 + \\mu/2 - \\mu^2 - \\mu}{(1/4 + 3/4)^{3/2}}$\n",
"\n",
"$ U_x = \\frac{1}{2} + \\frac{\\mu}{2} - \\frac{1}{2} - \\frac{\\mu}{2} = 0$\n",
"\n",
"#### Compute $U_y$:\n",
"\n",
"$\\frac{\\partial U}{\\partial y} = \\frac{\\partial}{\\partial y} \\left(\\frac{1}{2}(x^2 + y^2) + \\frac{1-\\mu}{r_1} + \\frac{\\mu}{r_2}\\right) = \\frac{\\partial}{\\partial y} \\frac{y^2}{2} + \\frac{\\partial}{\\partial y}\\frac{1-\\mu}{r_1} + \\frac{\\partial}{\\partial y}\\frac{\\mu}{r_2}$\n",
"\n",
"By chain rule with $u = \\frac{1}{r_1}$ and $v = \\frac{1}{r_2}$:\n",
"\n",
"$U_y = y + \\frac{\\partial}{\\partial u}(1-\\mu)u \\cdot \\frac{\\partial}{\\partial y}\\frac{1}{r_1} + \\frac{\\partial}{\\partial v}(\\mu v) \\cdot \\frac{\\partial}{\\partial y}\\frac{1}{r_2}$\n",
"\n",
"$U_y = y - \\frac{(1-\\mu)(y)}{((\\mu + x)^2 + y^2)^{3/2}} - \\frac{\\mu y}{((\\mu + x - 1)^2 + y^2)^{3/2}}$\n",
"\n",
"$U_y = y - \\frac{y - \\mu y}{((\\mu + x)^2 + y^2)^{3/2}} - \\frac{\\mu y}{((\\mu + x - 1)^2 + y^2)^{3/2}}$\n",
"\n",
"Plug in $x_0=\\frac{1}{2}-\\mu$ and $y_0=\\sqrt{3}/2$\n",
"\n",
"$ U_y = \\sqrt{3}/2 - \\frac{\\sqrt{3}/2 - \\mu \\sqrt{3}/2}{(1/4 + 3/4)^{3/2}} - \\frac{\\mu \\sqrt{3}/2}{(1/4 + 3/4)^{3/2}}$\n",
"\n",
"$ U_y = \\sqrt{3}/2 - \\sqrt{3}/2 + \\sqrt{3}/2 \\mu - \\sqrt{3}/2 \\mu = 0$\n",
"\n",
"\n",
"Since $U_x = U_y = 0$, the equations of motion have equillibrium solutions at $x_0=\\frac{1}{2}-\\mu$ and $y_0=\\sqrt{3}/2$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### b) Lagrange points\n",
"\n",
"The solutions to the equations of motion in the restricted 3-body problem are called Lagrange points and are shown in the sketch below. \n",
"\n",
"The L1 and L2, are near the smaller mass, M2, where the gravitational forces of M1 and M2 cancel. These are unstable equillibrium points, meaning that if some M3 placed at L1 or L2 begins to move away it will continue to accelerate away from these points. \n",
"\n",
"The L3 point is on the opposite side of M1 than M2. Any M3 at this point is also in an unstable equillibrium. \n",
"\n",
"The L4 and L5 points are equidistant from M1 and M2, so the x-components of the gravity from both masses cancel, leaving only an effective acceleration towards the barycenter. This force keeps masses in orbit in much the same way that M1 and M2 orbit the barycenter. L4 and L5 are technically unstable gravitational equillibrium points (i.e. a \"hill\" in gravitational potential rather than a \"well\"), but the rotational motion of the system applies a Coriolis force which acts to accelerate M3 when it lags behind L4 or L5, and to deceletate it when it gets ahead. This makes L4 and L5 the most stable of the Lagrange points, allowing them to capture objects of non-negligible mass, e.g., the Jupiter Trojans."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"![Lagrange points](lagrange_points.png)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### c) Tadpole and horseshoe orbits\n",
"\n",
"The tadpole and horseshoe orbits are named for their general shapes the are described and drawn below:\n",
"\n",
"Tadpole: These two orbits surround the L4 and L5 points. They are wider towards the M2 body and taper off as they extend along the orbit, giving them the distinct tadpole shape. Objects in orbit here appear to orbit the L4 or L5 points in a rotating frame of reference that holds M1 and M2 fixed.\n",
"\n",
"Horseshoe: The horseshoe orbit is a long stable orbit following the orbital path of M2, extending from near M2, past L4, L3, and L5, and ending symmetrically back near M2, giving it a horseshoe shape. Objects on this orbit oscillate between prograde and retrograde motion with respect to M2. The object \"catches-up\" to M2 in its orbit, is deceletated by the gravity \"hill\" (transferring angular momenentum to M2), then returns on its orbit along the horseshoe until M2 catches-up, pushes it with its gravity hill (imparting angular momentum), and sending the object back along the horseshoe again. This repeating pattern makes the horseshoe orbit stable."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"![Tadpole and horseshoe orbits](tadpole_horseshoe.png)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### d) Pan\n",
"\n",
"Pan's gravity clears a gap in Saturn's rings called the Encke gap. The presence of a ringlet in the same orbit of Pan probably indicates that particles are trapped in horseshoe orbits around Saturn due to the solution to the 3-body problem discussed in c). This is drawn in the figure below (very not to scale)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"![Pan](pan.png)"
]
}
],
"metadata": {
"gist": {
"data": {
"description": "AST501_hw1",
"public": true
},
"id": ""
},
"kernelspec": {
"display_name": "Python [conda env:ml]",
"language": "python",
"name": "conda-env-ml-py"
},
"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.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment