Skip to content

Instantly share code, notes, and snippets.

@hannorein
Last active January 19, 2019 23:41
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save hannorein/409e87d14791658bfee297b39a2218f6 to your computer and use it in GitHub Desktop.
Save hannorein/409e87d14791658bfee297b39a2218f6 to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 59,
"metadata": {},
"outputs": [],
"source": [
"import rebound\n",
"import numpy as np\n",
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"tmax = 1e5"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Jacobi coordinates"
]
},
{
"cell_type": "code",
"execution_count": 60,
"metadata": {},
"outputs": [],
"source": [
"sim = rebound.Simulation()\n",
"sim.add(m=1)\n",
"sim.add(m=1e-6,a=1,e=0.02,primary=sim.particles[0])\n",
"sim.add(m=1e-3,a=5,e=0.05,primary=sim.particles[0])\n",
"sim.integrator = \"whfast\"\n",
"sim.ri_whfast.coordinates = \"jacobi\"\n",
"sim.dt = np.pi*2./30.\n",
"sim.move_to_com()\n",
"sim.automateSimulationArchive('test_jac.bin',interval=tmax/500.,deletefile=True)\n",
"sim.integrate(tmax)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Jacobi coordinates, but order swapped"
]
},
{
"cell_type": "code",
"execution_count": 61,
"metadata": {},
"outputs": [],
"source": [
"sim = rebound.Simulation()\n",
"sim.add(m=1)\n",
"sim.add(m=1e-3,a=5,e=0.05,primary=sim.particles[0]) #swapped!\n",
"sim.add(m=1e-6,a=1,e=0.02,primary=sim.particles[0])\n",
"sim.integrator = \"whfast\"\n",
"sim.ri_whfast.coordinates = \"jacobi\"\n",
"sim.dt = np.pi*2./30.\n",
"sim.move_to_com()\n",
"sim.automateSimulationArchive('test_jacinv.bin',interval=tmax/500.,deletefile=True)\n",
"sim.integrate(tmax)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Democratic Heliocentric coordinates"
]
},
{
"cell_type": "code",
"execution_count": 62,
"metadata": {},
"outputs": [],
"source": [
"sim = rebound.Simulation()\n",
"sim.add(m=1)\n",
"sim.add(m=1e-6,a=1,e=0.02,primary=sim.particles[0])\n",
"sim.add(m=1e-3,a=5,e=0.05,primary=sim.particles[0])\n",
"sim.integrator = \"whfast\"\n",
"sim.ri_whfast.coordinates = \"democraticheliocentric\"\n",
"sim.dt = np.pi*2./30.\n",
"sim.move_to_com()\n",
"sim.automateSimulationArchive('test_dh.bin',interval=tmax/500.,deletefile=True)\n",
"sim.integrate(tmax)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Democratic Heliocentric coordinates, but order swapped"
]
},
{
"cell_type": "code",
"execution_count": 63,
"metadata": {},
"outputs": [],
"source": [
"sim = rebound.Simulation()\n",
"sim.add(m=1)\n",
"sim.add(m=1e-3,a=5,e=0.05,primary=sim.particles[0]) #swapped!\n",
"sim.add(m=1e-6,a=1,e=0.02,primary=sim.particles[0])\n",
"sim.integrator = \"whfast\"\n",
"sim.ri_whfast.coordinates = \"democraticheliocentric\"\n",
"sim.dt = np.pi*2./30.\n",
"sim.move_to_com()\n",
"sim.automateSimulationArchive('test_dhinv.bin',interval=tmax/500.,deletefile=True)\n",
"sim.integrate(tmax)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Comparing orbital parameters of inner planet"
]
},
{
"cell_type": "code",
"execution_count": 64,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEKCAYAAAAvlUMdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXeYXFX9/19nts22KVuzNZveGyl0MHQRpAgKiuBPilIE5WtBBQQFFFEBaYoNwQLSQSmhJPQAKSQhvW/K9jLb28z5/fG5Z++dTXYTIJsC5/U8eebOnXPv3Lvifc+nK601FovFYrHsaXz7+gIsFovF8unECozFYrFYBgUrMBaLxWIZFKzAWCwWi2VQsAJjsVgslkHBCozFYrFYBgUrMBaLxWIZFKzAWCwWi2VQsAJjsVgslkEhcV9fwL4kJydHl5WV7evLsFgslgOKhQsX1mqtc3e17jMtMGVlZSxYsGBfX4bFYrEcUCilNu/OOusis1gsFsugYAXGYrFYLIOCFRiLxWKxDApWYCwWi8UyKFiBsVgsFsugYAXGYrFYLIOCFRiLxWKxDApWYPYl3e2w4K/Q1bqvr8RisVj2OFZg9iX/Phf++z1Y98q+u4YnvgXLn9p332+xWD61WIHZU0S7ob1x1+tevw0+fEK2y9+R186mwbuugYh2w9KH4dEL9s33WyyWTzVWYPYUr9wItw6F2nXyvqsNoj3xa7SGV2+Cx/6fvFcJ8trZ3P95o93wr3Og/F15394Ayx7b9fU0VUDNmoHXNFfGv2/cIlZVe8Ouzz9YvPV7qF27777fYrHsMazAfBKaKqBqhWxveE1eFz8or7cUuEJi6PtA9zkC0zGABVO5DNY8D09dKu8f+yY8fiHUrd9xbXsDtNbJ9uMXwT0zBxaZpm3x75/8Nqx+Dja91f8x5ppisYHXfBzaG+Gl6+ChM/b8uS0Wy17HCswn4Xfj4L5DZdvn9A1trYOeTtle+Uz8+to+D3uzrrNJHq53TYc518WvqV4pr4l+eV3/qrzuTGBuLYPfT5XtxnJ5XfXf/q8/sjX+2je/Ka+xnp2vBxGXPxwBb/y2/zUflzZHHAey6CwWywGDFZhPhJaXWMy1TtrqoLXGs0S7216BifZA1CMwDRuhbh28/XuIRd11VR+639XT5e6vWysW1Es/EzeasVRMPCfJ715PfxgLJjk9fv9ALjLjvlq/G4kJ5fNdC68/OiLQ4vy9WqrlNTFl1+e2WCz7PVZgPi5N293ttlpoqXK269wHJcS7xcwDPdEP3Z7U5M5medAavAJlLJWGzfEP/tq1Evd56w5Y+xJs84wd0NoVlr4Cs+kteOce59qca47F4oWwY4BkhXrnerpa5LVqhcRN+rJ1Ifz1RHHpDcTvD4LfjJRtc98JyQMfY7FYDgiswHxcjOsKxG2kHaujrwVTv8HdNsLT0wFt9e7+jqZ4gfGKkjlXdys0ekYwNJa7v/QbNkFzhftZW70rRq218df9wMnw4k9EUMyarub4WpwBLRgniaHBuZb7DpW4SUcE5t8HNwSlvmfTG869OWK1/En4+xd3DOC3ea6v1fn7WIGxWD4VWIH5uHitlMql8hoocqwZz2deUTBWDsQH2Dub41Ocveu8AmHiKiAP7sRU2a5Z5VojIK447QTh2/oIjPe8XkvFez0DpVs3bHSuuSl+XWM5vHars2aT+73Rbnmdcx1sfM0Vnr70dLr3amJCtevEAuruiF8751o3k667HebeAp0t/V/zR2HuL2Her/bMuSyWzzhWYD4uLR4rw/wqzxouv+S91sSiB+Hxi504jVdgHBebSpCHdX8WTFstBEtk2whMuEwe7sb9Vbcu/nqqV7jnNllljeXxD+HIlnhLpXGLu93eIBbOmhflml++Ad660/3MULfO3W7YDCjZrt/oikVbXXy6tlcwTZIDiKXWa621yesjX4MPH/fEoRzevksy6QA++JcIm7m+gahd6/49+uO1X8G8X8p2LAov/tS11iwWy0fCCszHxWulmAdt1jB5rV0DSemQWSjFlMv+A9sXi2VixMJkcGUW7OgiMxZMZ7M8bMNl8r5XYIaJ9dEbZ6kXIUjPk/dGYHJGi0C1N8Adk+GXRe53RLaKSBl3lNeC6WiUAP2/vgx3TII3b4eXrpfanvZG18KIeESpsZzepIf6DR4h0SIcJjPM6z70xrFaqly3oflb1KySV6+l5BXJrjZ3rde6M5TPh3mOVVW1HO6eAf/97o7r+mPrAnjnbnj2qt0/xmKx9GIF5uPSXCkWiy/JfbiFhsprwyZIDUGo1F2/8hl52OeOkffOA311V5hoe0QelKlhCJbC3Jvhj0fDL4tlbdZweTXfkzXMsWCch3h7g1gw+RPkvYn75IwUgVr7Er0Pf0PTNhGSoPMdxupKTIW2Blc8oh4rY9MbzjF9RBIk+G9EpH5DvJA0b3eFwCvM3uNbqt3jO5vjM+laPcd43YdVy92kA5Oh11bvZq799USYd4u40Uz3hPqN9Is30SEWda01Hd35eovFMiCDKjBKqZOUUquVUuuUUtfs5PMUpdQjzufvKqXKnP3HK6UWKqWWOa/HeI6Z7uxfp5T6vVJKOfuzlFIvKaXWOq/hwbw3WqohYwj4g+7D2ePKavZl8sQHMRY3OinAb93ButYknlveTkwDkW3M2ZxF7IE65lckiED4g8SmO8WZFR+432UEJrIFUCJcOuoKTnu9CF7OKHlvLIOwY1Etf8rtGgCSxWYsGCOC5pjsEXItXjefoW49RLsgtBOBqVji1s+0VIl1le1cT+06egXO6yLzWk3NlZ76Fx3vJvSKine7o9F1TzaWi0D8erhbm2Ro3OIKXKx7x/syeBMd2uqgybm/pLT+j7FYLP0yaAKjlEoA7gE+D4wHzlVKje+z7EKgQWs9ErgdcPwZ1AKnaq0nARcAD3mOuQ+4GBjl/DvJ2X8N8IrWehTwivN+8IiUQ6YjMAZjDbRU8caabsa9tJWWt0Nw9I9ojSq6n8lm2JMb2dqZBE3bSFshtSpNTUlQv5EXVitWfutP/KfzRPjmHADeq8/gmTecGERjuXxfapa8N7+wo12SmRYqFZdXkyMOxmW3bSHkj4fjboQv3gUZ+SJWPe2uKJoHes5osT68MRmDsRKCjigZgckdG59V1xGRc2Q76ccNHqvBa9l4Baq9Ib7A0ut+6y/tuyPi3mt7g+NSc4TM635r2OQKzM6E0+BNiGipcv8GxiXYHw+eJl2xB4PuDjdRwmI5wBhMC2YWsE5rvUFr3QU8DJzWZ81pwN+d7ceAY5VSSmu9WGttnhDLgVTH2ikAAlrr+VprDTwInL6Tc/3ds3/PE9kmD/ui6b0C06jTePqaa5hfmwlA13b5NZ/VrGkYdw4LAyf2Hr4pkkascQsZTiw72pQItWvp3NKDT4Pvg41QejCVB11J5pwAo/7+Fq1RJSKSGpJ/BmOlgFhUKQHojMR/1lIpnx3xXTjofMjIcwszjQXT7Py5jQtv+2L5rORg+NJf5D6NtdDXgskd69bFJKTIw7mnwxU4EyQPlkLtavjNaDl/0zYRS1+SJDp0Nsv1Q3xMpaUfF1lnk3uvOhqfAr3lXXe7cbMrMB0Rid0YmqukLshbO2S+x4icN6W8Lx0R2DBPumJr3f+6j8vN+fCnY3a9zmLZDxlMgSkCvD+Dtzr7drpGa90DRIDsPmu+BCzSWnc66z0/e+POma+1Nj9PK4H8T3oD/fHY737AnDcK6Cw6GPzyQHxrYYjRH9SS8GYm3THIruihxwc+De/Of5s1dZm9x0ca/NQ1NZPi/DBNbfRBtJP0RkktLqlsJRaL8cJq1/e/uC6TqIZn5vbw2B330GNagRm3GIhFleJ8j/K51glApufPkZEPNY7FYUTIWAI5o+V1+yJJu75wDkw6SxII6pwHuDcG40t0XXgA4aHuL38jXiZVu2SWvLZUwTNXilAHi+Rv2OGIhbECvQLjtWbikhEiIkrGoqr3tM+pWOJuewUG4q2Y/35X6oIqPogXktZa19rqL9Xb/A0M9RukK8Hz1+yYWg2w8ll4/Tf9n6s/TBq8xXKAsV8H+ZVSExC32bc+ynGOdbPTn5NKqUuUUguUUgtqamp2tmS3KNmmeP71+bSQytpWP8PXwbZwAhkdsLY1layI5v1xktW1fcl7JG9YT0O6YmNuEr7GBCo6JHurJuAj1KCIashp0HQmQqBds3zpe3SvW9X7fY1NqZS3pzBqVYwJr2/h/XpHSIwbChyXnWMBJGdCek7vRz2pOTx94gwe/H8niwUDtEYV65ti4lZrq6U9Bi89+xxdMaSOJlDonjsjz7UejAXTViuWjXddaKi43gDSc+XcDZvk/ZjPu+uqPpQC1UCxWC0djY5YOAJjRKV4lohFZCu88guxvLJHSkypqUKuM+wkV9TtRGASksV91hFxrSNzH9Fu181Yszo+PtTZ7Fo6A7Xb8QpMRyO8+nN49z5Y/b8d1z5yHrz6C9fS6WobnKahFst+wmAKzDbA8xOaYmffTtcopRKBIFDnvC8GngTO11qv96wv7uecVY4LDefV41dx0Vrfr7WeobWekZub+7Fu7PDLb6crEcbcfC9L71rClojEUtYcf6xcZIMfH9AxeiwxBbHKbYRqG6nITqMxkIa/FRrakgDYWJJFWhdsbE8hpQc+GC2isGrh22Rur2BLTiK1mT5irQlsafb3XkP9JifwbCwOEMvEPERTMlm3egkN3RLcf+GtFYze3MrMdzbSlSIxnNffy6frohuI+ESs5m7Op/jvb/FGheOCKznEPXe652/lzY7zByXV2mAe9gD+kFyPsTrKjpRXX6IIQ/N2KHWswOYq2debKOEIzKgTxN12+wR44zew9kVJ207JdB/uxhXn7ZpQsURcb+EypxVPo3tt7Y3wwo/hpjwRObPeKyRdLW7Qv61essoql4l14hUFr3XVEXGtoL5uNSOyIKLW3SEdt1+5MX5dxRIRUq3jraCYzWSzHHgMpsC8D4xSSg1TSiUD5wB92gvzDBLEBzgLeFVrrZVSIeB/wDVa697e8Y4LrEkpdYiTPXY+8PROznWBZ/8epyAvn1eOmAxAuBmSlqYSU3DIGecB0FUj1kl6yTDqM3wk19WTE+mkKRygNRQg0AxtrRI4bh83CYDNDVKV3zZlBgBN61YxpKaZqpwg9ZnJJLUomhuT6PHB4uEBMut8dMdgTUTxSl0Rz789hPKKyl6BeX6Vj+5zvs3KF/Npiyoim1y9fWdNJU09Pso2y//8722T7+5uEjFqcTLfOosP4dEvHcmcR+4T6wgAJfU9yRny1h+EvHHuHyfg8YL6A65FBZCWDd9dBj/cAEf+Hxz7MzjsSkeEHLHoa8GM/YLEeEwtEIiV5A+4a8I7EZi2OhE+f9Dt9WbO0d4A8+8VQTOxo5rVcowvUe6xs8XTL85pq/OnY8U6+fBxeP8vIgJxFkzEdcXVrZPmpEsekYJSb1FqYzmse1m2Fz1IHP/8sghpS3V8pwXjrvMWp1os+zmDJjBOTOUK4EVgJfAfrfVypdTPlVJfdJb9BchWSq0DrsbN/LoCGAlcr5T6wPnnVBFyGfBnYB2wHnje2f8r4Hil1FrgOOf9oHH+bx7kr/93M01pkNOoqAkppkw+iM4kCFfKgzq7dCT1mcmEaxsJt2o6c/PozsomowOikUQiaZA9RgSmvU5EqWzWUbQlQ8KWTeREorTm5hIJZpDerFFNidSEFA2FhWQ1KuZtCxO9/FYKX9KUlft4/cffpk2lENUQXCRJBuEWWJc0iuTadlqd1mUbq3pYpsf03kubk5iVWS2V+P6qRDjldv77+KNMXF5L0u339D6cN3T4WbzobTeTLS07PgbjzarzBz0WVRASEsX68Qfh2OvhyKtlJo4/KPEYcAXKWDDhMrj8XbhqCZz5J2ffUDmfsYyMBdNcAWk57jWYmFRbndQDmTqlvmMTQB7gbbVyPymZrgXjd6y55kq3JuiJi+B/V4vQtHjcrB0R99w1q2W0wpOXyGyevm11Nsxz7yUOx31Wuzo+btS0Tep+bsqDVc/teP0Wy37IoMZgtNbPaa1Ha61HaK1vdvZdr7V+xtnu0FqfrbUeqbWepbXe4Oy/SWudrrWe6vlX7Xy2QGs90TnnFU68Ba11ndb6WK31KK31cVrrAVJ/Pjk5GSncdvGZVOSKMFTlppOQkEBDQKwagGGjJhIJZjByu7g6EgpLSciXB2i4IoHGoI/C4WMBSKsWURo5ZjIVWSkUb6wgQYMaUkhbVphQMyS3KJoDSVA6ggQNeqM/7pqmf1DB6pvnM3dbiHALvHSoxGeWpMwmp76NVUPD9Pigp6qaFcjcmJpMH0nNmpiG7EZ5uOU0dsOMb9L52jwA0jti9ATEsmj6Xxj/+ZdS73MevGk5oJR7EV6BSfFYMGkDlCX5g259SmpYjutpFxdXUqq7bvzpcPh35dUfcF1agWK3zic1JK5CgECB40pzhChYIskPW9+X91PPg2nnwdhTJP7SVi8Ck5whFkxXq+sO7NuuBmDzW07hqbOmI+JeU2SLmxyw6r/xLrPGcqldgvh0anAFrWZ1vCi11klnAoBVO4nveLHuNMt+wn4d5D8QqMmX+EPsi+fK+ywnTTkRhg8bQXvYfbCGh48lo0R+7Wc1Q30owNjxUwAoqFV0JMGw0jLqsgIU1YkFkl4yHJ1XQFIUCqoU7ZkpZI87CICS7Yr1eck8e9IhPPXVMwHwd0PWexKfOeiiq+nxQfealeQ2RmnOz6c2MwF/TS0JmzfSlgwbi7PIbIK67kSSo1Ad8JHWBTVVWwjUi/sos0PzYUUzFZ1JvZlv76yX2TSx1CwidRVw5Qfw7TddiwXiLZjUAQTGe0x6tueYULx4JSbD8Tc6FoyblYc/CGmOReX3CEymIzDmQZ8alrVGYI7+AZx2jxSXttWKyKRlQ0qGZLT1dLgWxs4Epn6DuM6CxSJcLdVusWlzlZs0YBINDG11rni0VEmnAYOZclq7Nt6C6WqReicQS7A/utvh51kfL1vNYtnDWIH5hJz0yz/xzmXn8eX/J/2qmoZIRtW2cBKpyUlEh7lB+MmHHseY6Uf1vm/JHUJ2OIu6DPmfoS4zkaTEBNry3ZTiISMmkFokDzkfoDOTmTLTPUfVsCJ+eMffuPyHN/DqPU+wtiCZ1C5oT4ZDDzmKynAiBevKxRrKL6I2lEqooZlQdR3bs1Noy8km1KypdLLaNhXLg3r18sVkN3WyPSwPvPWbtrO63q1o76iTh91Tf5zD+tnHUN7YDkMmeWI1iPVhLJqBBMYbpwkNdY8J9XUfeTAuOnN8mpPd7g+6CQmZBTsKXmpYHtQqwU0oSM+VfQ0bXQvG1N6Ya6haLq+BYhgyGUafBHUbRCiMcJnU6mCpjEAw6dmxHrFYUgKytm9z074tc0BcdnEC0+qmV/e1UJq2wys/l6w641p89Rf9/+0MtWuhcifCabHsIazAfEJKikr45pU/JSFBHsRq6AgAOlLkV2bZIcf3rh1VWsLUcW4zg9RJ0wHYnitB9fqAuLx8pW5MY+zEaeQNG9v7vmPIKMaUuoH0hEnTAAj4k7j82HFUlMlDsytRkZCQQH0wlbIqEYPUwmKawkFyIl3kNLZTnx2E/CISY1DRKN/dNloC9ls/XESoVbNhmFhokU3r2eLEbbZkJ5Jc182WjiTGrWojtQvm/cYZ9Zw/0f3jKOUG1vUA6bjeDLW0bLEoID62s8Mxbgo2/pArOIGCPhZMH4ExLqhAoWstmO9vqXItGJPKbETICMwlc+Hbb0DBVHGDtVSJpeUPukkHpjbJ9ESL9YiIpDpZdR1N4loz12xEJdoTPzbaG+TvanWFqKVPguST35IR1iuejq8TaquXjDdT6NrVJnN5TJr03TPgD4fv+Lf1suo5yZwzLHts11NKLRYHKzB7mMziMgAawvIr/IgjZgNQnp2IUorEBPdPfszZMu2x2+/UxIyRB3jOSPchXZSXx8jx03rfp5SOx+dz3UbTTvxS3PfPuPynAGwc6dS6hN2YSHbJcDrz8gi2aXKbYnRkZ5NRIg/zllpZV3r4cQB0LBY3UmzcJDoTgW1biEai1Gcotg0Jk1UXpbLNHW0cWu20g/H5YMaFMMyxskY4VegVOxYLNtVX0tXeJh0RDEpBniPCxirZGUZgfEkyHtq4zALFkOEIhonBGIwFA/HZbmkesQoWSQ2RGa2QliWi0FIlVo9ZmzUc0E4tkCMwxnow3RCMKIE85FOzxNoyFoxJMTdi1lJFb5C/s6mPwDS7AtO33Y3Z39HYR2Dq4JUb4M7JkqTw6i/g0W84saMmd91AHQgePlcy50wX68cv3LHXm8XSD1Zg9jCnnHsJL5w4k8Nv/TMAgdRkHrjqeipv/HPvmre+fwXPnn4MRY47rXXWEQCMP1OyrKcfcQLvl2Xyj4suRSnF8FLXVXT4KecA8L+TDmF5oZ/JE13xAZg+41DCLz7OKfc/CkB3Tl7vZ2UjJ5BU6J5LFZRQ5GSxFW3poisRjj72VLoTYMh6eWhlDRtNVSiJtJo6Ag0RaoN+2nPyyGqGxlap5VkwKkRpVTvdXU7dxim/44mOw9i0fgUUTpPU36N/xKt/uYlnzpNaoWi0h5XHHsMLZxwtLicvxnLxJgz0xTzoE+QaejO8gsVu0D1UuqPAmKB9hvt3iXPrFU0XC8YUiyaluXGYjHwRUHCtLHBdZCbWYyyYLk9vtYZNTgJDUB7u7Y3uusf+H/zhSLefmz/kDqFLTBUR7Wp1EwW8AqN1fD82b9JAa63MywERHuPCq1oBW99z1w1USGqoXhE/18di2Q120cXP8lFJS0nhe3fG1zbceum5ce8vuujyuPff/MHNbPr6VcwuFndUWUE+bX99ma/mS61JYoKP+cODNPr9fK9YXDZX/+6v9MQ0yhsEdxgy1HXDJRcUAYsAGD58JKuGu+62QNkoxk2cTj3SPaA8J4kpaanUZSYwtFrcakPHTmZtOJ28+haSe2JsKs5GZecBy2lqCNPj6yAybRYpa+ewfMHrbJ3/EpEPlnLQe+W88epLlD3yKlxfx+r351BwmwwFK1+3lI2rl5HXrglsauGNef/jyC/8lpeXbKD8+su44LrbSehsgqlfo6uzjeSUnXQzNm4tnyMwJlAeLJKCzlCpuOf8i91j/EEYMRsW/i3+QWzGHAAUHiSupt4/YLr0Y6tcJuMPDF73nXGRGfI8PV1TsyT+0tUs1lB3u2SIxbrjz1G5VNKeQVrqVK0QK8cflGSDrlY3JmPm66RkiogYMWxvlLUG7/jupgp3kNu2hfHC27A53uXoJVAk4lS5LL5rhMWyG1gLZj8gMTGBkcUFcfvGFwbi3GnT/vgCZ9zvpqf6fIrkxF3/z1c0eSYAtZk+0lJSGD5hRu9nQycdwpDsMA0mySBLBM3EgqIKxo2bQktODoX1PeREYnTk5JCSJ9eatb2b2swEhkySc6757S8Ycf9/Oeg9+aVcVF7N24/cSXdXB+t+9KPe731vzjNsePHZ3vcbFrwDMy8i+rdnOfQ/c3ns1v+DWRcz98kHWD9lOq/947c73li64z4zGVUmxuMPiVVTerC893YdSE6HYUfL9qxL3P1KwTn/khRofyBeLJIzYNLZsn3UD939ad4kg5BYJobMAvd9wRR3v0nBNpZE38SHJf8S4cwa7sZgUkNOXY4jMH2LSiOenm3tDRKfMdfvnaPTXOE2A23YKC44g+l2vWFefGKBuf+dHWOx7AZWYA4QppaEKMtJ/8jHnXDy2Sy65RY2/OwOAMaOHM3a/GT+PW0EM6cehFKK+kyxAtqz5aHZHJKHSl3AR2Z6OrHCEunyDPhKRxAsKgNgSGOUupCfkRMkbXrScnERbcxN4r1RYQoaooR/9gf+d/uPKa3s4NmJhfT4ILJkIQlVVTSlKmJAd/kGtldWUNggLpjoBxKv2TDvVQDa/uxOa1i/9A2i0R734WwsmNPuhYMu2NHdZmbSgAhJaghuiMCUr8SvG/sFSYGG+PhMSiaUHgLXVsOwI+OPGXWCTC4tnNanwDTgDikr9bTbMTEYU/PjPWbcqYCCCWfId3Y6bjR/UISxrU7cgIWOS7R6lQT2TU1MWrYIUkeTaxl5h6tFtnpiNU3xYtFSLWsfPA2e+0H8PZo4UGdLvPjYHmqW3cC6yD4DfO3MM3q301KSqPzFE5w9PKfXQlox63BGPf0qieMlHtOZlQ1U0epkwmWWjQbmApA3fhp5xW78oSkrxJiRo1mVJDU4i0aF+NITb/Lumi08+Mvvcf6CVYSefInEGBQcdhg1W58kpaaWtNZ2KsMphJK6SKmq5M0XHmNSDLoTIK9Cakb8VRL8zmyRB/KalQuJfvkSHj20jHP+9JQ8SE+4WS4kdzR88fc73rw3BXp38XahNvGXxJQd1537CKDdbgSGlIDb0qXkYHd/ajg+my4tC07/g2xPPVdayyQmw1t3ynmbtotbKtbjBu8Lp8HyJ6RDgJfiWdJNOtEvGXJJ6fGtc6pXEJ9A0CSxsViPiMjmt+WziqXw5KUw+ycSzzL1O10t8YkB7fX9u9WevhzWzIELnolvI2T5zGEtmM8glxw1gknF7gPxh7+8m8RnH+e8q28CYOQ532ZLViIbJ8rDYcKRX+hdO3nmUYwc5cZxuvLySU1OIsF5bjaWDcOflMDRE8r49p3/4I2J+eQ3yq/5Gcd+gcbMFDIjrYQjHUSC6VSHU8mqayKyTgLc74/Lp6C+h7aWRrJr5RdzuCVGZ0crC1+SVnZT3tlEdwy4cjGMPXnP/4FMPzSI7yTQF5/PTXX2CowvAWZeKNtmRAE4GWme2EdajgjLVCdGlyjZhL2p1ZEtrgVjYkaBQjjmOhEub3p3Zr4ztK1Jjk/Pju8ubepdMvIdC6bZde21ewSmZqW46h74gsRsTHFnV2sfq8czl6cvi/8h7jnvyATLZxIrMBZ8PsWoUeN7a3mO/dzxHP36B1xyl2QgTR0zkicnl7KywE9ZURE5QdcqSCqUGEdbiiQbTPj693o/G5qdTvIhR/S+HzNpFs3BDMLNnWQ1R2kPh2kKh8iKdKG2ldORBO0TpuHTsH7lIvIbumlLlv9IN6z5kPZlbsB+4TviPrv/iq/XI9zBAAAgAElEQVTw5Ekz2bxugIfZVUvgioW7/wcJlex6TV+MwJiYxYm3wDXlIg6G1HB83KU/C8CIULRLXHrJGW6mV0oAjvq+zOm53JMJlhoWoehoEqstNex2cPaH3IFy4WGScNDR6KwLyrYpCjVEtsR3HzANQw1t9VJX8/pv4htwelOeTUzohiCs6Nvn1vl8zZyd/w0snwqswFh2SkpiAkmOC00pxVfvfYKcP85BKYVSirdGS7wmq0xiHIuuvo6Hv3EBn5s1M+48J5xzEQDlBSn4fD46QyFym2IkRyGWm09Xdg6hNk3W9koqQ0kES8oAWPHOa/i7YW2xiNnGVcsIbqmgwwm5rH37Zbq7uzny5aWM3dTC3If+2P/NhMviM8B2hRELb33O7h7TGxtK2DHNOjUr3jrqr87HFIOa8xrR8n4PiEV04i0yXtsfktiO6RjgD7nZZbmuxdnbGLSpQoTMHxJh8o6iBnHleQs6d+Yie/suqa354J/ufq+V01bvpl7PvXnH+3ziEvjX2Tv2YxuIzW9LsaflgMAKjGW3KMtJ56jRrkvm6Hse5dkzj+fE078OwGXnncuN11yzw3GFxWU03vYjpvz9P7Ij122Dk1o6El++1AKN29xKXXYmeY5gtS0Ri6NhpMR76jatJdDcwYqyEO1JEF29ggXvv9l7rqQVUtT4xtz/8o/f/vST3/B3l8H5H2Hig4n1eMXBoJz/m6WG49vfmBqevnhFKDUcP420r2gderlkzHktI38gfqx2nkdgTBZaZKvTuiYkFkxLVfxcH4AaZ+BdZkF8mjSIeJgstY6IzK6JReNFqb3eTZM2DUe3L4aXbxBLZ4tjgRn33fPXwGu37fRP0svfPi/FntHugddZ9guswFg+FqNKCvnhLb8nLaWfh6SHQ0/9BkNKpWo9zel0AFA8dmpv80+Alrx8ho+WRIPsTfKrNmuGuNi6tm8lq7mHtnCIiqxkMqrrWPPuPEAadBZuFRdS+3U/ZvqfnuDZB+/8ZDfYt0hzV5j+YOk7sUqynKQIfzBePPrD28I/e5QrChAvHF68+1P6pFrnegLtvRbMVtfSadoulsfMi6RTdbETNzKdCLJGSBZZZ5Mrlu31buHnkofh5nzJQPNaQu0NbnJCV7Nknt3/OXjzdnG5mYSHqmUiOO/eB3Nv2uWfB3Abllr2a6zAWPYqUz/nJgxMOehgCke6RYmqeChlRUW0JcOI7eLemXrwkTSlKlK2biWtC6LZOTQE0wk3ttOxXn5hrx5TQlZTlPr6OgqdLtStjz+yF+8KSUee8lU4dSeZbOc9Dif+UlrYJPl3/Lwv3sSC/PHxBZkZ+Tuuhz4WTDDeksr1TD31JgakZIowmR5qmQXw5b/DcT+T96aDdNYwx4JphIwhkqnW3gB1Tl2NsXQ2vekG/1MCIkBe95fpdACS4WbcadWrXEsH4pt5Rntg8ztubMd0adj2EWJqln2GFRjLXmXCeLdOJTccZvT4g3rfh0dOJCnBR31GIj4NPT4YMWIsdYEkireIhZKUV0hLdja5kR5SampoSFd05xeSHIVXH/8rCRoa0hVjNjTQ2d6y924sMQXOuG8nA8SQfYde5r4vmOLUvewGgSLX6oD48QVe/H0sGGPR+JKkO8Gwo+GLd8WnbfsD8ccZV5xpw1O1QvqyZeQ5iQERx/2WBW0N0OAp8gQRDdPGJmd0vAUDTqq0Q81qd7u9Pv4z73jpVc/C306SKZ/e++902vAs/ofbzNOy32EFxrJXUUrx6OzpPHGYWC6FWUFWFsiv+lFTpYliRY64piqykkhMTKIxkEZ+RH7VBovLiOXLfJzirTXUBVJIduI4kcXvArD4oLH4u2HBWy/2fm9Heyv7Dd96Hb7yj4HXnPlnGSmtlNTlBIrglDv6X9+fBeNLELG54Bk46Pz4jgNp2fEWkTNUrtfKaa0WcUnOcIP+pmFopFxm5kyTMeHkT5Qkgy3viRsta5gIR0u1m6Sw/QP3u7xTRTsi8SnV3vHSJqaz4AF5NcLS2SJZbE9fLnEZy36JLbS07HWuu/eh3h5qSiliNz3I7f97grtHSzC6uWwYrGmgzhGa1nAIkIryYRNm0NzcBLxBUV0PC8cFCRTJL/zQBvklGz7yRHhjJeVLFxEpX4cvNZOhN97Fe984lQuu+XXvdcx/4l7a62qYffHP9tKdfwQmn+1u+xLg6l20yPfGYHJGSYwFdhx54E2NTsuO71pgrK/UsDTZ7GkXATIC0bRdMtISkt34TNlRIoSRrfCn2ZLllRqWc7dH5NjcMeLS8g5tMwKTlu2OLzB44zhGYHw+cZMZt1pXs2cqqMdKsuxXWAvGstfp26DzrMMncf8tP+vtrRYcL+Ocm4fLwzHmidNMmXwQxRPc6vi2rCwKhkl7/OHbWmhPhsOPP4MY0PPBQob9+gGG3ngXACVPu7Psl7/1LMGf3MWQ3z5M5eYVxGIxXj5yMu9OG8+2LZ5f0AcK3lTmtCy3NsWbogw7pjl73W8m9uPzSewHpBLfJDtEtoiLLDTUM666UKwc0wizzUwFzXTdatkjxarxDjczVkqoVNZ4M9S88RgT0+lskUahZmJoZ0v8GGrLfokVGMt+x7kXX83i73+Lr9/yJwAmnei2uklMTGD6dFdgYkVDGTl6EjElrWrqMhMpHZJLTdDHjPfiffMFDVH+e+WXaInU8sHDbsfruX++g5UfvEVRTTeBds3bLzwFyEiBl++7jobqPrGG/ZG+sZnRJ8HI4+CEX/S/Li3bHQjXFyNEo0+Kr9fxB+MbeAYKnf0B17VmpoLqmDOQLQvS86RLgMFYML0C0ySxn+RMd9Q0uBZMZ5/+aV2trgWzMz74l3SABnjsQumz1tXW/3rLoGAFxrLfkZCQwFcv+i6pfonNHDrrcJaWpDHnjM8BkJXhZllN/MJ5DAkHqAnIf8qNmX6UUmzNlV/d28MJDFu6jEW/lvqKEXNW8OJ13yZr6RpWFKUSSVN0btzEstdf6j1nZLmMN3jkzhsouvMxnv+JG6Cvq9hIS8TzANyf+PKDcMlrsp2WJdlrA6VFp2W7VktxfIEsn79NOk6PmB3vVksJQIGnoai3dsa444wFA2JxpGTK8DdDqNS1REJDxRXXVueIVE68BWNqbWI98e1puvpYMF4LqHELPHUp/NtpwfPhY9Ipuvyd/v8WlkHBCoxlv0cpxVdeWshVv7yvd9+zEwtZl5/MYTNn4fMp6pwRAw2FErSuHSlus8aAH39yIqefeBL3HjKNlhQom7eckuou6kcNozaQTHpdI20rltCVADWZPtI2ieWT8vzzAKRulQdbe1sT1bNPZt55p+y1e/9IjD8NCqfu/npjmfxwI1zwbPxnOSPh5NukGNQ7nM0flHHRY0+BU26HZM+sHtO5Ojmjz6jqgCtEyRlizYC4zUwMyBR+pufGC0yzR1RMsabyOQPZ6nf8DNy5Or5Ed04QiDUUi8nYaG8qtOHtu+DWYQNP+LR8JKzAWA5ITr7zCYJ/erl3fHRiVIr2fKNleFjeEdIEc4OTFp2eksjFN93Lm8cfQVoXJGgITJpKQzCdrEg7KTV11AQT2Z6XQbi+hUhzEyMqxKVStq2FaLSHZ+8Rd9OItX1mphyopDozbdKyBm7q6R0p7Q9ILc85/4QZ34xfN/FMeW3aJlNBDSmZrsCkZbmp0t50ahPfSc+VIH/5fMksa97ujmAw4wYCRWLBeHuleettTAp0QnIfd1uViM+j34D57o+VXuZcK6Jlkwb2GFZgLAckE4uCca1rqkPikimYIvNXvnT6WTz5k19y3s/u7l0zuTjEhNkn9L6fMftE2rKyyGnqITPSSiQjhZZggKzmbubOeZqUHlhekoa/Gyo2r6J5hWRy9figq93158cOtNkoRzvD33an6LPvOm+hZl+Gz4aDL5W5Ot4uCCkB10XmD7nWjT/oGY5WI9sZuVC7Gv56Itzl1EhNcGJwpiFnsNgJ8nsEptkjMKZwtGl7n4SBarfY08SAti6EhQ/E34fJkNsdtN65NWQBrMBYPiVM/sVf+MeXvshxx38RkGadPzn/dApC8eOWP3fcKSyZWsDKrxzBiPGz0PkFJPfA0OoOWgLpdGfnEGjX1CySPmeVU8TltHHVMtLq5IGWGIMVi9+go62Vp0+YwZzZ0w4skZn9Exm69nEYenj/n/l88PlfSZPQvgLjcyoihh3lWjDpuX3m6ARhyCT3fdEMGHm8O1PHFGAGS9wgv6njaapwjzMC09UcPxOnpcqto+luk35mfz4Gnr0q3i3mTafuy+a34bfjZNgbwBMXw+0TbW+0frACY/lUcPCEUdx8862kJCYMuC45JZVzHn6VM2+UDLX0UglMJ/dARzDU23wzbeUKYgpyDz4WgKr1qwg3ttCYJi65dUsX8cITDzC6vJWhVV0sXfTWYN3a/kV/Iwb6EtcBOgCTz4FpXxfryVgwmUP6ZKgFoPQw9/3Fr8B5j7kiZCyYUAl0t4r7KyNfhKppG2x6Cx46Q4TIDI0zM2kyC0RgjPg0lsfPq/G6xbydBPqydYFYSyudmNWyR+X96uf6P+YzjBUYy2eaglETerd1Ti4ZTtFmWXkdtZk+Rk6Sxo9tWzeTE+lmQ5E87FoqtlC7clnvscvmPAnA6y8/zZMnzqB8k6dS/dPAVUvl3+7iDfKbLLLT7o5vT5OWHR/fSQlI3U5GPpxwU/zxIK1pkjPc2FFjucR0QkMltvLQGbBe5gT1WkKmQ0DeeHGXmWSAhk3xouLtJNDuFH2+cy+sfTn+vkxBqLF4UjwuPssODKrAKKVOUkqtVkqtU0rt0MtdKZWilHrE+fxdpVSZsz9bKTVXKdWilLq7zzFfUUotVUotV0rd6tlf6hyz2Pl8EEYdWj5tjJ80o3c7ZUgRBSNEcLJaNHVBPyNKS2lJgcSNm0jvhKbhI4gqiFZX4yvfSFsytCVDz2qJz0Sv/SljN7cy94+/3un3HbCEh+68z1p/eF1kfWttzPjp5PT4FjeBAnGzfX8NHPYdd7+xmjojIi7mmPr18v6EmwAVn9E2fLa8NpZLsD9QKKnMJhmgtTY+88wkEIAkD8Ri8OKP4Z9fir92Iz5dpvWQGUO9F/veHUAMmsAopRKAe4DPA+OBc5VS4/ssuxBo0FqPBG4HjGB0ANcB3+9zzmzgNuBYrfUEYIhS6ljn42uB/2itpwHnAPfu+buyfNooHOLWZ8w6+atMmnpw71CzSFaQYGoStYEkSjfKL9SMYSNpTPeR1FBPqKqObdl+qkNJZNQ20NTcTJ7TMy20cPEO3/WZwohIsGTHDDWTOpyUJoJiCPRTs5OcLgWYIFlnxq0W6xELZuih8IO1MrnUUOZMUm3a6iYTdETczDMdhWpPMN+4ztKyRUQaN7mfVS2HRU5hrjm+q8W9Bu97SxyDacHMAtZprTdorbuAh4HT+qw5Dfi7s/0YcKxSSmmtW7XWbyJC42U4sFZrbezRlwHzE0MDxi4PAh9hTJ7ls0xlMIGtWYlMGDOG3ID7K7h9+GiUUtQF0shuliB+/shxNGQkkRZpITvSTkNWJnWhDLIa2lj6wbu9XaDz6tr7+7rPBkrBt96AS9/e8TMzabOv8AxUFNrb6TkrPm5j3GVJqSIiF70q/4wF1d7gzsfpbpMssr7xGZDiTJDPOhplhIDhvsPgme+I1WNcZJ0tkj3W3ea+749YDP71FVj3Sv9rPqUMpsAUAVs877c6+3a6RmvdA0SAfubIArAOGKOUKlNKJQKnA2aA+g3AeUqprcBzwHd2fgqLJZ6JL7zF9OfdIH13ggTyh58kleAtITfTadzkGTRlphJqaierOUZnVlbv+IDyDxcA8GFZgFCrpiVSxz1XncejJ8xgydLP4ICsgsnx4wEMh1wu4wpmXhi/P9j38eDB1NGkhkVkDN5tgOLp8s/rovMH3ZhQrAeynQFwlctcgYo47YBCJSJKtZ5xAoaWqngLxmSkgWSs9UfdOljzAjx+Uf9rPqUcUEF+rXUDcCnwCPAGsAkwSejnAg9orYuBk4GHlFI73J9S6hKl1AKl1IKaGhuYs0BuOEhW0H0Q3n/qBTw4YwxHHypulp5ct5K9qHgYreEQRXU9JMaAvCHECopJikLHImlFUjtWPMHrVy5k6luLmFjeyvKffW/v3dD+TkaujCvwxl9g4AmixuWWNSJeVFKzdr6+bxabV+hMxwFw+6qZGEywVMSjsTw+UQGkjsbEYNDxrWs6W2Q42tNXxM+6AdgurYfiOld/RhhMgdmGa10AFDv7drrGsUiCQN1AJ9VaP6u1PlhrfSiwGjDpOhcC/3HWvAP4gR1yKrXW92utZ2itZ+TmDlA0ZvnM8qsffZezbv4bqcmS8pyYL66b+gyFz+cjWuj+Z51WWErGUOkkXLh6I+3JEJwmtSJr3nuTUKsEgaetrGP7lrV78zYOHM5/Bk7+zcBrjAtr9El9Ztr0IzAJiRLjgfiCTnA7P4Nb19MrMI6brnZtfIscEBdZe4Nb0+PtHtDVIkkHix+SEdJetjvxONO14LkfwNqX+CwwmALzPjBKKTVMKZWMBN6f6bPmGeACZ/ss4FWtB24EpJTKc17DwGXAn52PyoFjnc/GIQJjTRTLRyYnI4XpQ91f10eceQEvTy0h8b6/ApA5fEzvZ3nDxlA0Rooxh1Z3UxVKpmS0tDbpfn0ePg2vHSIumUVz/7u3buHAYvjRMOvigdeccruMpC6aHp8Y0LdJpxdjEXldZBA/oqBwqghRtEtm4JhOBTWrdhxP3VIpHZ1NDMc7t6azxU0U2N4nwcOIV1udFGS+dz/886wdizM7IvDiT6Wz9KeEQRMYJ6ZyBfAisBLJ8FqulPq5UuqLzrK/ANlKqXXA1UBvKrNSahPwO+AbSqmtngy0O5VSK4C3gF9prY0F83/AxUqpJcC/gW/sSqwslt1h4pixfOfhORw8U9rQDJ/kPtQmTD2YSVNmEHW64NdlZTBi5DiiCko2izEePvkcYkD9YjcOs2GLJy3WsmvKDpeR1EZcktIl/dk7aK0vRmBS+rjIUsMw0xG0IZPcupyUTNciaq0RC+aKhXDRK5LqXLlMRhDkjJY1pj1NalgsGCMk2xfDvFvhv45b1PRDa6l25+iAuNJaquHZ70pm3dL/wDt3w+tOinu0+4BPfx7UiZZa6+eQgLt33/We7Q7g7L7HOZ+V9bP/3H72rwAG6GNhsewZpkydyUtDUlgzfSY/LCxGa82qDB+5zTHacnIpDKWzOcNHjpN5NvuYE1l6+69IW7eBOUdOJq+hm5QeeGhGCV//x5ze87a3NRHt7iYjOFCeiwWQWhnfLh5fiU4PNX8o3kXmD8Lnfy21NplDRKSat0uDzswh7rqMfOkqDWLZbFso23ljYe2LrgWTWSABfyMwHY0w7xbZPulWt/9ZW1281dNcCe/9UVKgSw91591sdjLv/nO+dAjYVVuf9XPFqsoZOfC6fcABFeS3WPYH0lNTOe7FBfzgN/cDMk5gc57za7mkjMQEH43pUkzTlKrIy8mlIiuDCWsjlNSIuADMWLCF+ipxqzTVVzL/80fy4bFHsWm5JAtUblrJi3+4ce/e3IFCSsaum3WaqZ6pIddKAREYn88tHPVaMF63mLexZ6DIncKZ5zhTTKFm5hAnMWCLjBLwUr1crCHlA7TbZBMcl5tjofR0QP1G2TYjCkz7mUZvMm4fIlvhodPhkfP6X7MPsQJjsXwMMlIS40Y/Z/zsQR790hc47pKfAtCUIQ+/ukwRmsZcN9/knWMm8PZ08ePPv/AsVr8/h3kP3sGQqi6CLTHmP/RHADacdRaldzzMknc/e/UTe4Tp35BmmRPOFJH53E9gxDHx7WnAdbOlBOKz0rx1OcYtBpDvtBcyPctCQx0LZgsUTos/95b3JLaS5aRGGxEBaK6Q2I/ZNucztTa953g3/n20G1odV9vS/8hrQhL7I1ZgLJY9wBmzRnP9zb9heIH86m3LFIsmEpRMJl0ggtLjg2/c/SizbvobnYkwbF0TG3/0I+o+WEpHEtRlKHo2bmbzhhWEW8TFtuipf++DO/oUcNgV0iwz3XE5fu5H8PUnJcPMi7FggsXxCQTezs7G/aR8bpC/YaO8zxsvQrH1fZnqedKv4GuPS6r0Zqe+yrTL8TbSbK50M9Eat7jNPLtaREQSnNTsvuMD/vVluM2ZHmpiOgPFovYhVmAslkGgvkDcL7XTpFmmf5Rkmj155DR8PsX4oYVce/qVvDYqi8KqDnI2b2dLnp+K7DTC1Q2sWuImBPidOTQrly3gtYMn8tITf93Ld/MpJzldXkOl8ftz3GxBwk7mWeFBYukon1Tx+4Nu4Wa0S0TqkEth1HGyv3y+c3yZvBpRyh0rAtPoFHg2bpYaG+Nia9gMUcfF1+YZmhaLug09uzvcYs/eMQTtrmtwP8AKjMUyCJx41a955KYH+Nb1twNw8hlf59qvXMWZ18skxQSf4okbvkXiEUeTFIWRFZ3UF+fTkJfDkLpOqtetBKAylEC4Qfz0c/92F3mRKOm3/Lb3e9YtnkdDzQA+esuuMb3RTLeAIZMka8xr6Yw4Rlxu5/xTrByTNJAadgUG4t1q2SPdYkzjYmssl2MyC6Sbs7FAmp0U6JATF/LOpPFO5az11FJ1RDwC48RyfjUU/nDER7r9wcQKjMUyCEwdmsUNZx1McqL8X6w4nMbjN36biUVuNlNyoo9px7hNv5PHjEMPKSStC7pXygNm/dA8chu7icVi5DiWTKA1Rk93F91dHXSfeymLz/zCXryzTyGdTpaWSWW+6FW4po9o+wNw6p1ulpnpQuAPuS4ziN/2FnQad1vTNhGnzAKoWkFvN2ZTQ2OsKJNQkJAcPwrAO7mzo9HNPDNCE+2URIL9pELDCozFsg+ZNt0dsDXy0KNJKZBfwPkbt9KYpugqGUpaF1RVbKKkWlrEJ2jYtmkVi96WFOeCmm56urv2+rV/ahh9krwWiBuTxORdZ6gZgUkNgy9BrBsU5I1z1wyZ7G57Rx34gyJUxgUWLJUsMu86IzhZI+IFxju502vB9O3mXLt/dI2wAmOx7EN8Ph+VOZIBNOWQ4wiXiK9/eEUnNcEUUpwBaO+89CyhNs3yUokXbFi5hBXz3BqaZfM/G61HBoWpXxOLxevq2hW9AuME10+9E66riY/jlMxyt7290VICrjsO4utXjIvMxGayhrsusm0L4wP+fQUmFnU/21mzTi8DTe3cg1iBsVj2MZP+8yRd995ISmoGpZ4Jm43hDHJHSs1F4+siIFVjJfBcvWENsfXretdudKZr/uM3P+Yv/3f+3rr0TwdK7bzr80AYwfDW1/RNFTa9zEKlbiIBOBaMp97G60rrTQZwMsqyh4t4VH4IfzoG5t/jrvUKDEhXgJ1tA/ztC/Dg6bJdsQTuOQTe/eOAt7gnGNRKfovFsmtyCkeQUyi/nkePHo/z25X27GxmTpS2NKWr5IFT8LkvwJxFtG0rJ60hQlOqItCuad5WTjQaZfqfnwJg23c3UFQyfK/fy2cGY+2U7SKgfvUqcbf5EqTnWXebiJnXgtmZwBgLxnR+fvHH7pq0HMksa2+QGExiqszY8cZnWvtknm1+030/71YpVJ1wxm7d6ifBWjAWy35EeprrSvEVFDOipISmVEVRfQ89PjjuhNNoSQFVXUWwqY3NQzLo8UG0qpJ5c57sPfbVh+7e2ekte4ojvw/fWw4Tzxx4XaDAdacZqycl6BZegrjBDOk5IkSxbum3ZsRn4+vumpKD5bWjUSyYgCNW3u7OrR4Lpnqlu93VJh0EhkzesVv0IGAFxmLZz3h/nBQGphaWkpzoozqYDEBlOJFARjp1gSRS6yNkN3fTmh2iPsNHcl0d5Uve6z1Hz2p5qCxa9A73X3Y2bW0HdtPE/Y7ktIEncO4MM8HTH5Tiz3MfhkOviG9J4w+6bjd/IL7z89HXQNEMOPk26bPWVCHNN82cGa/AeF1k3smd9euhrb7/MQd7GCswFst+xnF/eIyXTz2c075+OQB1Ifnl2+B0BWgIpJJb20ygXRPLzaMhM4X0SAvdFdJscXNOEsEqcZEs+vW1HPnqhzz600v2wZ1Y4nFSh028Z8zn4cSb46vwUwLxrWsyPM03D70cLn5FJn+mht1AvRE602xT+eJdZKZDAEj6c3t9/4Pa9jBWYCyW/YzC/CF857Y/k5nuZIzN+jwA7UERmuZQkIIGyRjyF5cQCaYTbuokqbaGSJpiW0EW+bWtaK0p3iYprrmLV+yDO7HEkT9RXnUsfn+wVCyZi1+VWI3XgvG2rvEmIqRlu5MzTTzICEzWcCnwbNgMd0yGFc+4bWeaqyQ5wFowFosF4NuXXc2j532FUT+UCv7ubNelMnTKwbSHwmQ1R8moj1CXmUTbkCKyWjRrVi9laI0MtSqq7bS1Mvuaw66U19yx8ft9PrFkiqbLe68FA/Ct1+GyPg0v07KgyREUkwhgameGHi6usP/9n1gvNSuhwKnJaXCabe5PFoxSatKuV1kslsGgMJTK9dfewGHTpBAwscTNOpp2yGx0bh5JUSitbCESTCepSGoxFvxPRvcuLcsguQfWLptPLKZ58e336Y7Gdvwiy+Ay9FD48TYYeezA64wFY6yMgikyg8ZLmmdmUI4RGEdwZnxTXtd5aqOyR0rSQN36+HMPMrtrwdyrlHpPKXWZUiq46+UWi2WwmP2Vb/Vup2cE8TudmgPtmtacbMJDpe9Vz6IFANRMkl+vaxa+zV0/vZzSb57Pn66/Yi9ftQWQ9OBdYar6+7b+9+IVmGCJtJRprgCUtKX52uNw/C9kqBrAsKMkblPvCExqeIdTDga7JTBa6yOBrwElwEKl1L+UUscP6pVZLJadMrq0kA++eixLzj8RgOyho3o/U0NHMHSsWDrFG2V6YukR0gqlcesWxrwp9RCT57xGT0+fmfCW/YM6p81L4UH9rzEC40uSQWmm+WZalsRxRh0Hh4Hzi9EAACAASURBVF8pInNtNUz9qiMwG9x1e4HdjsFordcC1wI/Ao4Gfq+UWqWU2kUiuMVi2dOce/3dnPOTOwAYNWlm7/68CQcxZvR4uhKhsL6HjiSYdeQJxIBobRX5jd00pSqym2O8NefxfXT1lgE56VYoOQSKdkdgEp1OBMat1meYmlKQ6AT4vdlq+1kMZrJS6nZgJXAMcKrWepyzffsgXp/FYtkFY0eMpCZT/q889dDZhNNTqA5Ik46qUBJ5oQBNaYqMzVvxd8N7B08gqmDzC0/ty8u29EfZ4XDhi64w7AwjMKbljLFgvDU1ffEKzH5mwdwFLAKmaK0v11ovAtBab0esGovFso9QSpH613/z6rfOp7RI4jFb8uVhUlmQhc+naExPZES5tKXPmXEYlaEEkrZIUFjvJ63dLR+BUcfDEd+DC56V90Y80rP7P8ZYOb6k+Oabg8ju9iJ7Umv9kHeHUuoqrfWdffdbLJa9z8xJk5k5yW0P35GdDdQSC0iqa1OGnzInZXn4pBlsDP6bUH0z9197KQVvvkPD6KEkd3dxzt+e3xeXb/mo+INw3A3uezOnxj9ADpaxWlLD4jrbC+yuBbOz9qzf2IPXYbFY9iDpx58LQMLsLwHQmiFdAHp8MGnydCLhIGVVXRz52DxGVnYy8/U1THlnE22tkX11yZZPgplp0zjAdFMzSiC690YqDygwSqlzlVLPAsOUUs94/s0F6vfOJVoslo/K+V/+Mg3/e4Pzvyq/DZvzpV9VW4oiLTWNzly3XfymISlszxJnxpybv8PmldLT7MVfXsaLpxzc+96yHzPqRBh/Osz+Sf9rTIPNjr33I2JXLrK3gQogB/itZ38zsHSwLspisXwylFIcNsLNKEqfcgi8vIhAu8RbEktHAQsBGPPQf9lcVUPP+V9lzBPvs2z5VTT99BeU/n0uAB/8+48M/fksFs75N1tefJpTf/0PEhLspI/9isRk+PLfB17zUQaq7SEGtGC01pu11vO01odqrV/z/Fukte7ZWxdpsVg+GUefKLM/VhfIKODRn78QgO4EGFZcxMxJk/jBFy9l3uhsijY08u4//gJAJFURXSRTFOtv+hVj/reER2757j64A8snJuA0xfTtvR8HA36TUupNrfURSqlmeluBykeA1lp/xDFwFotlXzCqtJjfXHkFM2cdDsAx44q48ce/JklHmawUGSmJPPDjS7jnV5vxr3mOsgUfUp+u+HDEEKatqaS1uYG8eullphYs2pe3Yvm4+Hxw9gM79kIbzK8c6EOt9RHOa6bWOuD5l7k74qKUOkkptVoptU4pdc1OPk9RSj3ifP6uUqrM2Z+tlJqrlGpRSt3d55ivKKWWKqWWK6Vu7fPZl5VSK5zP/rXr27dYPjt8/7LLOXqGVPn7fIobLziVa79xeu/nQ4J+Tj/rbACK6nuozk4jNqSAzA7N2y8/TbLjs8itlTG9WmtemTuHuopNe/M2LJ+ECWdA3ri99nW7W2h5iFIq0/M+Uyl18C6OSQDuAT4PjAfOVUqN77PsQqBBaz0SKdg0gtEBXAd8v885s4HbgGO11hOAIUqpY53PRgE/Bg53PrN2vMXyEZk6dSZtMt+MtiHZJA0pBGDr2/MAWFaaTn5DD12dbTz80B8pvPQqPjj7i/vmYi37Pbubpnwf4B2J1+rsG4hZwDqt9QatdRfwMHBanzWnASYy9RhwrFJKaa1btdZvIkLjZTiwVmtd47x/GfiSs30xcI/WugFAa12NxWL5SCQkJJDmdPXPO/wwAsUyUdG/Zg0AVWPGkBiDVUvnU//OqwAU1nbT2W4nZlp2ZHcFRmlPua/WOsauM9CKAG9S9lZn307XOEkDEWCAUlTWAWOUUmVKqUTgdKQBJ8BoYLRS6i2l1Hyl1Em7uD6LxbITNnx9NhtHBZn9zR9TPFL89aVbG+nxQfggieFsXPEhWZs29R6zevGbRKM9PPHlo3jiyrOIxew4AMvuV/JvUEpdiWu1XAZsGJxL6h+tdYNS6lLgESCGpFGb3LtEYBTwOaAYeF0pNUlr3eg9h1LqEuASgNLS0r105RbLgcMXfnpv7/a4cVOoBLJaNVXBBIpHTQAgsnUDI6tb2JiXzLDqLtYvnk9yMIdxS2tgaQ0rFr3GxBmz99EdWPYXdteC+TZwGLANsUQOxnlID8A2XOsC5KG/rb81jkUSBOoGOqnW+lmt9cFa60OB1cAa56OtwDNa626t9UZn/6idHH+/1nqG1npGbu4AjeEsFgvhcDb16dJWpD6QwqixU4gB3RvWE27VbBorv+9aNqzlwwXv9B63dqnNNLPs/jyYaq31OVrrPK11vtb6q7sR43gfGKWUGqaUSgbOAZ7ps+YZ4AJn+yzgVa8rbmcopfKc1zBiSf3Z+egpxHpBKZWDuMz2upVlsXzaqMpOBaC+II8hWUEa0xUF68X7nTl2Ak2pCl1TR/2a5b3HRDbI776enm46O/deaxLL/sWu6mB+qLX+tVLqLuLrYADQWl/Z37Fa6x6l1BXAi0AC8Fet9XKl1M+BBVrrZ4C/AA8ppdYhrWfO8Xz3JiAAJCulTgdO0FqvAO5USk1xlv1ca20smBeBE5RSK4Ao8AOt9YDWkMVi2TXRhATZGDGKBJ+iPjOZkZUiGqUTptCQ/jQpjU10qHI6kmSp3i6dmh8+azYFVY3MfOE1AsGBwquWTyO7isGsdF4XfJyTa62fA57rs+96z3YHcHY/x5b1s//cfvZr4Grnn8Vi2UM052XDxmbyxk4CoCac2Ssw46ceyisZfgJNbcRiMapDSWggraaeaDTK9FXyG++pe27m/J/8bl/dgmUfMaDAaK2fdepZJmmtvz/QWovF8unkjDse4vkH7uCcs6W9TH1xGayspTUFxuUX0ZKZTkltC75ojMZgGrGYJtTUxorli3sfMNHVK/s9v+XTyy5jMFrrKHD4XrgWi8WyHxIO5/DV792EzyePi8wJ0wFoSRXXWWcoRKhVkxfppj0rRGswk3BzN8vefqX3HFnbqgBo7+r+/+3deXwV1dnA8d+TnSxkIyEJCSRhCQlbgLCJIooLKiooCtSNKq61vmprxVZxaX3r0tdaq9a6V60aQK2oqKBFERQiYNjCFvawJWSFACHJPe8fM7mEmIQL5OZmeb6fz/0wc+bM3Gcyl3vumeU8LFu+lMqjdR9xU22Rq3eRZdvD9F8nIlfUvNwamVKqRRoy5kp+6hrEl1dOA8DRuQteBjocBRPdmaPhEYQeMpTbF/2zUiKIzz+Mw+HgzYfuIOiaqWTeMsGTu6CaiavPwQRg3T58bq0yA3zY5BEppVq0gT3i2fHkx0zvHQ1AWEo6YPVWOiR040hxKbCakM25HPWBAz16EbhhCfm7NxO8PgeArjk7cDgczl6RaptcPbqvGmN+WfuFdQeYUqqdERHGD+xCaAfrlrGMM4797hw4+kKCYq0HmBN2lVAY4o1fQg8ANqzMIi7fSnYVVeZgc44mMmvrXG1g/u5imVKqnemdnOScThkwkpju1pi2EQcNxR0DiOphPf2/bflSYoqrWZVsjZu7ZukiAKodjT76plqxEz0HMwLrCf4oEal9+29HrGdblFLtnIiwbsa9hEXHkQoMGDySfd7gWw0HQ0MY0m8IAH4/rcALKE7rA1uWUJa3naxVK1n9wK10uvtPXH7+eR7dD9X0TtSD8QOCsRqikFqvMqwn75VSiit+cTPnnncJAJ3DO1Jh/3Q9EhtHclwMhcFC961FAMQOG021QPW+PWS9+SxnbC6FP9/nqdCVG53oOZhvgW9F5E1jzPZmikkp1dpZw5eRdtU0Any9yQ/1J3WXdWtyesYIdgR74VtYRHSp9bXSa/cRVv30A/0HjvBUxMoNXL0G4y8iL4vIPBH5b83LrZEppVqtn+66m1nnDObMkaMBKA6zEuAe8YWEbj0oCvEjsPQACQXlrE4MBmDFuydKMaVaG1dvU54FvIQ1sGS1+8JRSrUFt9x4K9x4q3P+cKcoIJ+SIG+8vLwoCwmk254SwsoNq8/oTn7Rany36UmStsbVBqbKGKM/L5RSp6Q6NAKATb27cQ5wKCyUqA3WNRnf2Hh2R24hMt9K3WSMQUQ8FapqQq6eIvtERO4QkVgRiah5uTUypVSbMfAXv+Gfw/ozbMYrAFR3inYuC+3ancKoTsQWHuUfM+7k3UuGM+vbZSzK2eahaFVTcbUHU5OzpfatHgZIbtpwlFJt0aj0FM584328vKyeiW/nY9nTE3r1Zcua1fhlbSX9s/8SVm7g1uvYE+ZN/gf/IbpLD0+FrU6TqwnHkup5aeOilHJZTeMCEJZw7Osjre9AApKsVABh5Yb8UOtrKbakmq+ffbR5g1RNyqUGRkQCReRBEXnZnu8pIuPcG5pSqq2Kt5/2P+ILQYHB9Bh4pnNZwf0zCJz7IRti/En/ZBmL3nkagB8/eZVPb7+U4oKdHolZnTxXr8G8ARzFeqofYBfwJ7dEpJRq8/r0G8ys9CQ+/tUDAAxO7UVhsLA/xIsxF4ynW3Iq3w+yG52//4uDpYUE3P9/dF+Qy/d//b0HI1cnw9UGprsx5imgEsAYcwjno1RKKXVyOnbw49aXP+ShW64DICrEn4/S+vBev0GEB/sDkHz5PTwxahSRpdV88Ohd+DisdYO+yQZg+eIvyLziTL766E1P7IJygasX+Y+KSAesC/uISHegwm1RKaXavM4dA46bH3n3XwnyP/aVdOuoZFK8fgELF9Lrm5+oFvi2fxfOWrOL6uoqNj71R9I3FJGV+RZMmNrM0StXuNrAPAx8ASSIyL+xMlxOdVdQSqn254pB8cfNiwgjh5/BWl8IO2TYHeGDJPfAd+UudmzJobM99H9s3n5PhKtc4OpdZPOBK7AalfeADGPMN+4LSymlwM/Xl92drFNmhdEdCUmwUgOsylpM5xJrUJEu+yspKy6w6uzZSuGerZ4JVv2Mq3eRTcB6mv8zY8ynQJWIjHdvaEopBZUZ1h1nfmFBJKRYuWXyv/kCLwM/pkTgBaxe8T17CovJP+di1k641IPRqtpcvcj/sDGmtGbGGFOCddpMKaXc6pzp/8e2tAh63HQHfQdYuWVicq1xyw70HQDAntwcPvrnUwBElVRzoCTfM8Gq47jawNRXz9XrN0opdcpCI2O56MPF9D1rPJ06daasg9Bjj3WPUbezLgagfMdWHFvWOdfJ+vx9AIoLdrL8i7ebP2gFuN7ALBORZ0Sku/16BljuzsCUUqo+hR19AeshzTPOPI8jvsDe3XTct4+yDtbTE3tz1gCw7KpLCbz7fyku2O2pcNs1VxuYX2M9aJkJvA8cAX7lrqCUUqohpaFBABSG+BAe5E9+Rx8C9hcRU3SQrQlhVPhA9Z49lB8sI36v1dNZsehLT4bcbrl6F1m5MWa6MSbDGDPEGPN7Y0y5u4NTSqm6DkWEAVDSMQARobBjB8ILDxJbXEVFlxj2h3jjV1TM8h8XOdfZtWKpp8Jt11y9i2y+iITVmg8XkRP+JBCRsSKyQURyRWR6Pcv9RSTTXr5URBLt8kgRWSAiB0Xk+TrrTBKRVSKyVkSerGebV4qIEZEMV/ZNKdW6BKUPAqB67FgASsLDSNhfiY8D/JJ6UBQSQEhJObs25TjXkS2bPRJre+fqKbJO9p1jABhjioHoRuojIt7AC8BFQBowRUTS6lS7CSg2xvQA/grUNBhHgIeA39bZZiTwNDDGGNMHiBGRMbWWhwD/A+jPFaXaqKvveoweq35i0t1/BOBo5LGvopg+AykNDSG8rIIDO63nYYqChcBC6+trx/YtvP7ANA4cKPn5hlWTc7WBcYhI15oZEemGPWxMI4YCucaYLcaYo1jXbi6vU+dy4F/29GxgjIiIfUpuEVZDU1sysMkYU2DPfwVcWWv5H7EaqbrrKaXaCC8vL3z9jg0z4x3XzTndN+MsDoeFE3rQwO6dVAtsiw0lrMz6Slhz/QRGfLSYD59/vNnjbo9cbWD+ACwSkbdF5B1gIfDACdbpAtQeVzvPLqu3jjGmCigFIhvZZi6QIiKJIuIDjAcSAERkEJBgjPnMtV1SSrUFnXoNBKAgxIuozvE4omLxAjpvz6MoxIsDkRFEllVRXV1NbNFRa6W1Kz0XcDviagPzJfAg0BurJ3IWUOyuoBpin5q7Hetutu+AbUC1iHgBzwC/OdE2ROQWEVkmIssKCgpOVF0p1cINGnkhv77gWlY8+hYAQV2s4WR65R2mKDSA6sho/Kvgp6xvCKi01onZvtdT4bYrrjYwLwLDgGB7qJgDWNdXGrMLu3dhi7fL6q1j90hCgcLGNmqM+cQYM8wYMwLYAGwEQoC+wDcisg0YDsyp70K/MeZl+264jKioqBPsglKqpUuJCeHVGXfzPxdZF/9jk3s7l5XGdMI/zvoaWvPlfwDIi/AhtqiSqkqrN/PJzNfIzV2HanquNjDDjDG/wr62Yfck/E6wzo9ATxFJEhE/YDIwp06dOcAN9vRE4L/GmEav7YhItP1vOHAH8KoxptQY08kYk2iMSQSWAJcZY5a5uH9KqVYssVMQItZDlil90mstSCKiWy8AfFdbp8W29IjHtxq2b8pm8dLF9JjxFzbcPKXZY24PXG1gKu27wmrywUQBjsZWsK+p3Il1em0dMNMYs1ZEHhORy+xqrwGRIpIL3As4b2W2eyLPAFNFJK/WHWh/E5EcYDHwhDFmo4v7oJRqB5Ljj13q7dJvIN1S+gGQtM0a1t8r4ywANq1cRvZrz1jr7KnQp/3dwNXxxJ4DPgKiReRxrN7GgydayRgzF5hbp2xGrekjwFUNrJvYQPkJf2oYY0afqI5Sqm3y9hKWPzYDn0UfMeHSa6jEl83eEF5uKA4SkgaPAt6mcFMOEbvynOv9+N2XXHDFLz0XeBvkUgNjjPm3iCwHxmClSh5vjNGTlkqpFunaq6fA1dZvUX+gMMSb2JJqCsICGJw6gEJvqM7LI7awnK3RfiTlHyV/4xrPBt0GuXqKDGPMemPMC8aY57VxUUq1JsUdraRlJREdiQsLpjDEG/99+USXVrOzVxJVXlC1fRtV1Q6+ytlHeUWVhyNuG1xuYJRSqrXamWA9J14dEoyXl1DYMYDu24rxNhDQM4X8UG/89+0jc9Y7hE4Zzd/umOThiNsGbWCUUm3eZX96jW8vGET/2+4HoDQ0hCBroGVi0wayP7QDHYsPUvrDPIIrYMLiHDZv3eLBiNsGbWCUUm1ecmwEtz33bzKGWneQHYk4NmBI/6FnUR4cREh5JYF5xwYfWfj2s80eZ1ujDYxSqt0xiSnO6ZjOXagI6UhYuYPO+cWsSQ7hQIDgvT6nkS0oV2gDo5Rqd0ZebeVLPOptzVeHR+BbDYkFlZTFRLM1JojOeTqU1Oly9TmYdqOyspK8vDyOHNEBmVuDgIAA4uPj8fX19XQoqhUZ1D2Oz//vcYJCQgDwqTXkv4mJYx/e9PlhI0cPH8KvQ6Cnwmz1tIGpIy8vj5CQEBITE51DT6iWyRhDYWEheXl5JCUleToc1cpcdMkVzumgzsee/g9KSKbcyx9vs5HPZ79GuXcwU6ZM1e+DU6ANTB1HjhzRxqWVEBEiIyPRUbHV6YqIO/YDJaZnHw55BwJf0evxFwH4/bzPOGfLRrr+7QV6DzzLQ1G2PnoNph7auLQeeqxUU0jLGOWc7t13EDE9+h23/Pola0nIr2Tpv19p7tBaNW1glFLtXlJ0GJ/fcBWL+nYmNqYLKb37O5ftC/fGIVDaQYhYtZbqan3K31XawLQwJSUlvPjii6e07rPPPsuhQ4carXPxxRdTUnJ6+ci/+eYbxo0b57btK+UJ9z7wGDfP/gaAhOgwNnb252AAHHrqn0TMncXXA9PoseMQ864aReXRI+Su/Ja5D02lrEiTlzVEG5gWxt0NzNy5cwkLCzul7bvC3dtXqjn4+3hzz/BHmXThE5x/xghik/pyeNwM3hnUk8ScYlYt+pT9N95O0qylfHen5pJpiDYwLcz06dPZvHkz6enp3HfffTz99NMMGTKE/v378/DDDwNQXl7OJZdcwoABA+jbty+ZmZk899xz7N69m3POOYdzzjmnwe0nJiayf7+VF2P8+PEMHjyYPn368PLLLzvrfPHFFwwaNIgBAwYwZsyYerdTVlbGJZdcQkpKCrfddhsOh+O47W/bto3U1FRuvvlm+vTpwwUXXMDhw4eb6s+klNt9ce+5LH1wLL7e1tfk/RenUpBhZRdZ9c6bhJZbuREjN+wDwOFw8PFzf6Bgz3bPBNwC6V1kjXj0k7Xk7C5r0m2mxXXk4Uv7NLj8iSeeYM2aNWRnZzNv3jxmz55NVlYWxhguu+wyFi5cSEFBAXFxcXz22WcAlJaWEhoayjPPPMOCBQvo1KmTS7G8/vrrREREcPjwYYYMGcKVV16Jw+Hg5ptvZuHChSQlJVFUVFTvullZWeTk5NCtWzfGjh3Lhx9+yMSJE4+rs2nTJt577z1eeeUVrr76aj744AOuvfZaF/9SSnlWr84hx82HBPjy/K+n8NPbfyb9x80ALEyJZNSGQgp25fL9vE/o9eKHzFv8PddkLvBEyC2O9mBasHnz5jFv3jwGDhzIoEGDWL9+PZs2baJfv37Mnz+f+++/n++++47Q0NBT2v5zzz3HgAEDGD58ODt37mTTpk0sWbKEUaNGOZ8riYiIqHfdoUOHkpycjLe3N1OmTGHRokU/q5OUlER6upW+dvDgwWzbtu2U4lSqpfDz9WFHXAgBlVDlBZxzEQArv/8v5R/PBiBha76zR9/eaQ+mEY31NJqDMYYHHniAW2+99WfLVqxYwdy5c3nwwQcZM2YMM2bMqGcLDfvmm2/46quv+OGHHwgMDGT06NEnNXpB3duD67td2N/f3znt7e2tp8hUm3C4exJsXsnuKD/6nnEOvPQOeWtX0aXAOtsRVeYgf99OYmK7eThSz9MeTAsTEhLCgQMHALjwwgt5/fXXOXjwIAC7du0iPz+f3bt3ExgYyLXXXst9993HihUrfrbuiZSWlhIeHk5gYCDr169nyZIlAAwfPpyFCxeydetWgEZPkW3duhWHw0FmZiZnnnnmae23Uq1F1/S+ABxN6ETf1L44gKp9e4guqWJvqDW42cbV1v/JyqpqqqoqPRWqx2kD08JERkYycuRI+vbty/z58/nFL37BiBEj6NevHxMnTuTAgQOsXr2aoUOHkp6ezqOPPsqDDz4IwC233MLYsWMbvcgPVm9j7NixVFVVkZqayvTp0xk+fDgAUVFRvPzyy1xxxRUMGDCASZOsxEvLli1j2rRpzm0MGTKEO++8k9TUVJKSkpgwYYKb/iJKtSwZV97O5lHJDHn874SGdKQ0SAjbthO/atic1BmAvZtzcDgMn12Qwewrz/ZwxJ4jxhhPx+AxGRkZZtmyZceVrVu3jtTUVA9F5F7V1dVER0ezd+/eNjU4ZFs+Zqrl++zsAXQpPIp/FXx7/WWc/dYclkw4k7Bzrqb3XXdZld7+G6lDLvBsoE1IRJYbYzJOVE97MO1Inz59mDZtWptqXJTytLKQQPzth/vTzr6YI75A/h52fjXLWWfNwq89E5yH6UX+NmrYsGFUVFQcVzZr1iz69evXwBpKqVNxODQEKKHKCwYPHML3IT502F/M0cNHqfAB32o4tH0bAEs/+xd7l3zL5X983aMxNxdtYNqopUuXejoEpdqFqoRusHwnRiAoMJDCjh0ILTmI15EKdkX5E1xeie8+62FMv+lP0KsSNl21mJ79R3o4cvfTU2RKKXUaEs+2noXxrbbmy8I7ElFWSefiw5R0CiU/NICQQusWZrEveWe9/ZInQm122sAopdRpOHuMNfDrj73CAaiIiKLjYUN0mYOjcXEUh4YQWXKEXbvznNdqHDvyPBVus9JTZEopdRr8/fyo/Pg/jO1kp12OiQeyAQjs1Yd95RWEHtzHj9/NI8VeJ7jEtefVWju39mBEZKyIbBCRXBGZXs9yfxHJtJcvFZFEuzxSRBaIyEEReb7OOpNEZJWIrBWRJ2uV3ysiOfayr0VEH6NVSjWL/ikpREdaPZiw1KHO8l4jzsVEdcYLKMr+HoBdET6EllmjZmzdvpUX75xMaVlps8fcHNzWwIiIN/ACcBGQBkwRkbQ61W4Cio0xPYC/AjUNxhHgIeC3dbYZCTwNjDHG9AFiRKRmuN+fgAxjTH9gNvBU0++V+7XUfDBz5szhiSeeOKW4atu2bRt9+/atd9m0adPIyck57fdQypMyzjjfOd2v/zACYroC4L9hAwA74zsReaAah8PB/Bm/4pyvVvLhc494IlS3c2cPZiiQa4zZYow5CrwPXF6nzuXAv+zp2cAYERFjTLkxZhFWQ1NbMrDJGFOThP0r4EoAY8wCY0zNt+sSIL5pd6d5tNR8MJdddhnTp/+sE9qkXn31VdLS6v4GUap1SY0L440Jl/L+mGH4+HjTqVsPALruLKLCBw51TcKvCgp2byVql5WsLGBZlidDdht3NjBdgJ215vPssnrrGGOqgFIgspFt5gIpIpIoIj7AeCChnno3AZ+fYtwe1VLzwbz55pvceeedAEydOpW77rqLM844g+TkZGbPtkaRnTx5sjOFQE29mmW1VVVVcc0115CamsrEiROdjeLo0aOpGVkhODiYP/zhD87RnvfZt3kq1Ro88fiTPPL8GwAk9bR67J0OOCgK8cYvxvrtu2n9KpL2WgPAdtndNrPAtqqL/MaYYhG5HcgEHMD3QPfadUTkWiADqHcAIBG5BbgFoGvXro2/4efTYe/q0477ODH94KKGTzW1lnwwe/bsYdGiRaxfv57LLruMiRMnMmnSJGbOnMkll1zC0aNH+frrr/nHP/7xs3U3bNjAa6+9xsiRI7nxxht58cUX+e1vjzsbSnl5OcOHD+fxxx/nd7/7Ha+88opzzDWlWjovr2Oji/fq3p2cACHkiKE0tANhXa2vrC1ZkzTBRgAAFMZJREFU3zGkCso6CFFlDg6WFRPcMdxTIbuFO3swuzi+dxFvl9Vbx+6RhAKFjW3UGPOJMWaYMWYEsAHYWLNMRM4D/gBcZoypaGD9l40xGcaYjKioqJPcpebVkvPBjB8/Hi8vL9LS0py9i4suuogFCxZQUVHB559/zqhRo+jQocPP1k1ISGDkSOshs2uvvbbeXDJ+fn6MG2fd/qm5ZFRrFhzgx+5wPwDKEuPp2sPq0ZjVKwHYmGj9H1vz0w8YY3jhqRmsXLvKM8E2MXf2YH4EeopIElZDMhn4RZ06c4AbgB+AicB/zQlG3xSRaGNMvoiEA3cAV9vlA4F/AmONMflNsgeN9DSaQ0vOB1M710vNIQsICGD06NF8+eWXZGZmMnny5HrXdSWXjK+vr7Pc29ubqqoql2NTqqXxtb/VOg4eTq+UNPK8oMt262vqUN90WPc129csxxEQxrmvz6LqzVlUrVyJj6+fB6M+fW7rwdjXVO4EvgTWATONMWtF5DERucyu9hoQKSK5wL2A8yqyiGwDngGmikherTvQ/iYiOcBi4AljTE0P5mkgGJglItkiMsdd++ZOrSUfTEMmTZrEG2+8wXfffcfYsWPrrbNjxw5++OEHAN59913NJaPavKDf/IbNCYGcO/FGwgIDKAr2okuR9aMpbewVAJRv2Uj2/I8A8HFAdlbrT7vs1mswxpi5wNw6ZTNqTR8Brmpg3cQGyqc0UH7eKQfagtTOB3PRRRc588GAdeH7nXfeITc3l/vuuw8vLy98fX2d1zlq8sHExcWxYEHDH86afDAvvfQSqamppKSk1JsPxuFwEB0dzfz5812O/4ILLuC6667j8ssvx8/P+vW1e/dupk2bxty51kchJSWFF154gRtvvJG0tDRuv/32U/pbKdVajL70Orj0Oud8YUd/ossOU+EDI4adzdIgwWfXLrz3HLuZJXfJN2SMvNAT4TYZzQej+WBavbZ8zFTb9I9rxzJ62XY2x/ozbkE2H543CG+g0suLgKNVJOyvIHtET254pWWeiNF8MOpnNB+MUi1Erz4AVPlaKZaLIsOILjpCZNkRSiI7sjvcl8C91jWanUWH2Fvq+vXRlqRV3aasXKf5YJRqufqcdzm8O5e9o6znzCpi4ghbuQcOVbMlvROlh44QWnaYgxVVvH3bBMqDQnn8jZkejvrkaQPTRmk+GKVarlFnjGLj/G+4pYs1QKZffDKwHADvmC6UFxWTUHCQ915/jgmrdgCwftt2eie2riEW9RSZUkp5QK+Ezs4HMqN7HRufLzypOxWhYYSVG6pX/eAsX/Tvvzd7jKdLGxillPKwtPRjIzAPO+9yHOGReBmI2LGDoiChNFAI+OlHD0Z4arSBUUopD+uecOzUV0yXJHyjYgFIzitjf1gA6xI7kby1oKHVWyy9BqOUUh7m5SV88atbiU/qSSoQEpsIQFAFlEWEcLBzHOE5BRws3U9wqGtjDbYE2oNpYdw9XH9tdUdIrm/k47qCg4PrLZ8xYwZfffWVy++tlDrePb++m6vGXQJAfO+BzvKKqGi8Y61hHbesW+GR2E6VNjAtTHM2ME3pscce47zz2sRgCkp5XL/UPs5p/15pBCf0BGDzmmz++MjvWb9jL63hIXk9RdaIJ7OeZH3R+ibdZu+I3tw/9P4Gl9fOB3P++ecTHR3NzJkzqaioYMKECTz66KOUl5dz9dVXk5eXR3V1NQ899BD79u1z5oPp1KlTg0PFvPHGG/z5z38mLCyMAQMGHDdo5cKFC3nmmWfYu3cvTz31FBMnTqx3G/fccw/z5s0jJiaG999/n6ioKKZOncq4ceOYOHEiiYmJ3HDDDXzyySdUVlYya9YsevfufXp/OKXakfAgf/ba08lDzuaIt3XmIOnZN+hdBcu/n8e+ksOEvPQKgwad4blAT0B7MC3ME088Qffu3cnOzub8889n06ZNZGVlkZ2dzfLly1m4cCFffPEFcXFxrFy5kjVr1jB27Fjuuusu5xhkDTUue/bs4eGHH2bx4sUsWrToZ+mJa3K8fPrppw1mrywvLycjI4O1a9dy9tln8+ijj9Zbr1OnTqxYsYLbb7+dv/zlL6f3R1GqHRs85Ex69Eijygv87UHFB+0oJ7rMwZrXn/FscCegPZhGNNbTaA6188EAHDx4kE2bNnHWWWfxm9/8hvvvv59x48Zx1llnubS9pUuXMnr0aGry4EyaNImNG53pdOrN8VKXl5cXkyZNAqxcLldccUW99WrKBw8ezIcffujaDiulnBb+fjp7Nq3lET8/YsK8+ToxlBFbSvn+3D70XraeiLJqui5fj8PhwMurZfYVtIFpwdyZD6Y+9eV4OZH6crnU3pbmclHq1Nx6/Q3OaX8fb94a9huyYj5h6r1/YXiPaB78n5u55stFzL1xLBe//gVeXl6sW/I5IZExxPcc2MiWm0/LbPbaMXfmgxk2bBjffvsthYWFzmsjJ8vhcDjvNtNcLko1n5n3XsrAGx5mSLJ1BiL6vBs56g3dl+xky8qFrPx6Jky9l3W/mubhSI/RHkwL4858MLGxsTzyyCOMGDGCsLAw0tPTXYopPT2d7OxsAIKCgsjKyuJPf/oT0dHRZGZmNtGeK6UaEx0SwK1nd3fOXzqiP1df+Gvemvt3sma+RdgPK0gC4nccorq6Cm9vz3+9az6YdpQPpq3SY6baqzW7Stk7cSRHfbxIKqhkS5QvyQWVeP37BVIGn+u2hkbzwSilVBvXt0so+2LCSCqoBGDrcGtMs9VLFrFu9VJWpfdj9pP3eiw+bWDaqGHDhpGenn7ca/Xq1Z4OSynVxKq6Jjinh4y3MsqX5e1gyUf/JqAS+rzxOQ6HwyOxef4knXILzQejVPuQmDEIvszGITBgyNls9AZTUIBPWZmzzvYtOST16NvIVtxDGxillGrFhl55K99u20Ro734E+PlQFOyNX3Ex0ftKOOoDflWw9sdF2sAopZQ6OR0COzL2oZed8yUhfgSXlRN5oJI1SWEM2lRC4Xrr9HhxeQWO8v1ERndpltj0GoxSSrUhBzsGEVN42BrqPzmZI77gyNuOw2H4+OrR5I86j9xVi5slFm1glFKqDTnSKZKwQ9bjJwFxCRSG+OBfVMKPORsYtrkEgJ/mz2mWWLSBaWFaWj6Y3bt3Nziq8slKTExk//79Pyt/6aWXeOutt5rkPZRq7yQ23jkd3rU7xcH+hJSWk/X1x87yI+vXNUss2sC0MC0tH0xcXJxLichOx2233cb111/v1vdQqr0ITe7lnO6W0o8DHYMJP3AUx3prSKld4d6E7NzTLLHoRf5G7P3f/6ViXdPmg/FP7U3M73/f4PKWlg9m27ZtjBs3jjVr1vDmm28yZ84cDh06xObNm5kwYQJPPfUUL730Eps3b+bpp58GrJ7RsmXLeP7553/2/k899RSff/45HTp04N1336VHjx488sgjBAcH89vf/pbRo0czbNgwFixYQElJCa+99prLo0UrpaBb6rGBLvsPGML34eGEr91H6J5dFHT0Ykd0GMn7SpolFrf2YERkrIhsEJFcEflZghER8ReRTHv5UhFJtMsjRWSBiBwUkefrrDNJRFaJyFoRefJE22ptWno+mOzsbDIzM1m9ejWZmZns3LmTK6+8ko8++shZJzMzk8mTJ9e7fmhoKKtXr+bOO+/k7rvvrrdOVVUVWVlZPPvssw3mm1FK1a9PX6uBWTykG97e3lRHdsbLQM+8IvZHBHI4tCPhB6ub5eFLt/VgRMQbeAE4H8gDfhSROcaY2t9qNwHFxpgeIjIZeBKYBBwBHgL62q+abUYCTwODjTEFIvIvERljjPm6kW2dssZ6Gs2hJeaDGTNmDKGhoQCkpaWxfft2zjzzTJKTk1myZAk9e/Zk/fr1jBw5st71p0yZ4vz3nnvuqbdO7Vwy27Ztc2nflFKWwMBgohcvZGpoBAD+8T2Ab4k4aMhNDacqvBN+VVspK9pDWCf33q7szh7MUCDXGLPFGHMUeB+4vE6dy4F/2dOzgTEiIsaYcmPMIqyGprZkYJMxpsCe/wq4srFtNd3uNL+afDDZ2dlkZ2eTm5vLTTfdRK9evVixYgX9+vXjwQcf5LHHHmuS93MlH0ztOrVzvUyePJmZM2fywQcfMGHChAbzxNQu11wySrlHZGQUPj7eAMT1H+Esr46NwzsqBoC8Le6/0O/OBqYLsLPWfJ5dVm8dY0wVUApENrLNXCBFRBJFxAcYD9QMxHOy22qRWno+mIZMmDCBjz/+mPfee6/B02OAc3j/zMxMZxoCpZT7ZPQ/dk0mJuMsgmO7ApC3ZWNDqzSZVnWR3xhTLCK3A5mAA/ge6N74WscTkVuAWwC6du3a5DGerpaYD8YV4eHhpKamkpOTw9ChQ53lF198Ma+++ipxcXEAFBcX079/f/z9/Xnvvfea7P2VUvWLCe1Anj8EV8B5E25g/jdfAVCyc7vb39tt+WBEZATwiDHmQnv+AQBjzJ9r1fnSrvOD3SPZC0QZOygRmQpkGGPubOA9bgF6GGN+d6Jt1UfzwbQNesyUatzEP79LZbWDjx+8lnVbd7L0jomYSyfyyzvuO6XtuZoPxp09mB+BniKSBOwCJgO/qFNnDnAD8AMwEfhvYw0CgIhEG2PyRSQcuAO4+lS3pZRS7UHm/VOc11UTu8TxzJQXuHZ4N7e/r9saGGNMlYjcCXwJeAOvG2PWishjwDJjzBzgNeBtEckFirAaIQBEZBvQEfATkfHABfYdaH8TkQF2tceMMTUnEhvcVns0bNgwKioqjit7++236devn4ciUkp5ireXANZNNR38vHnl+hN2PpqEW6/BGGPmAnPrlM2oNX0EuKqBdRMbKJ/SQHmD22qPNB+MUsrTdKiYeuiZtdZDj5VSLZc2MHUEBARQWFioX1ytgDGGwsJCAgICPB2KUqoereo25eYQHx9PXl4eBQUFJ66sPC4gIID4+PgTV1RKNTttYOrw9fUlKSnJ02EopVSrp6fIlFJKuYU2MEoppdxCGxillFJu4bahYloDESkATnVAnk7Az/P/tm26z+2D7nP7cDr73M0YE3WiSu26gTkdIrLMlbF42hLd5/ZB97l9aI591lNkSiml3EIbGKWUUm6hDcype9nTAXiA7nP7oPvcPrh9n/UajFJKKbfQHoxSSim30AbmFIjIWBHZICK5IjLd0/GcDBFJEJEFIpIjImtF5H/s8ggRmS8im+x/w+1yEZHn7H1dJSKDam3rBrv+JhG5oVb5YBFZba/znIhI8+/pz4mIt4j8JCKf2vNJIrLUjjNTRPzscn97PtdenlhrGw/Y5RtE5MJa5S3uMyEiYSIyW0TWi8g6ERnR1o+ziNxjf67XiMh7IhLQ1o6ziLwuIvkisqZWmduPa0Pv0ShjjL5O4oWVPG0zkAz4ASuBNE/HdRLxxwKD7OkQYCOQBjwFTLfLpwNP2tMXA59jZSsaDiy1yyOALfa/4fZ0uL0sy64r9roXeXq/7bjuBd4FPrXnZwKT7emXgNvt6TuAl+zpyUCmPZ1mH29/IMn+HHi31M8E8C9gmj3tB4S15eMMdAG2Ah1qHd+pbe04A6OAQcCaWmVuP64NvUejsXr6P0FrewEjgC9rzT8APODpuE5jfz4Gzgc2ALF2WSywwZ7+JzClVv0N9vIpwD9rlf/TLosF1tcqP66eB/czHvgaOBf41P7Psx/wqXtcsbKwjrCnfex6UvdY19RriZ8JINT+spU65W32OGM1MDvtL00f+zhf2BaPM5DI8Q2M249rQ+/R2EtPkZ28mg9xjTy7rNWxTwkMBJYCnY0xe+xFe4HO9nRD+9tYeV495Z72LPA7wGHPRwIlxpgqe752nM59s5eX2vVP9m/hSUlAAfCGfVrwVREJog0fZ2PMLuAvwA5gD9ZxW07bPs41muO4NvQeDdIGpp0SkWDgA+BuY0xZ7WXG+onSZm4vFJFxQL4xZrmnY2lGPlinUf5hjBkIlGOd1nBqg8c5HLgcq3GNA4KAsR4NygOa47i6+h7awJy8XUBCrfl4u6zVEBFfrMbl38aYD+3ifSISay+PBfLt8ob2t7Hy+HrKPWkkcJmIbAPexzpN9jcgTERqciLVjtO5b/byUKCQk/9beFIekGeMWWrPz8ZqcNrycT4P2GqMKTDGVAIfYh37tnycazTHcW3oPRqkDczJ+xHoad+Z4od1cXCOh2NymX1HyGvAOmPMM7UWzQFq7iS5AevaTE359fbdKMOBUrub/CVwgYiE278cL8A6P70HKBOR4fZ7XV9rWx5hjHnAGBNvjEnEOl7/NcZcAywAJtrV6u5zzd9iol3f2OWT7buPkoCeWBdEW9xnwhizF9gpIil20RgghzZ8nLFOjQ0XkUA7ppp9brPHuZbmOK4NvUfDPHlRrrW+sO7M2Ih1R8kfPB3PScZ+JlbXdhWQbb8uxjr3/DWwCfgKiLDrC/CCva+rgYxa27oRyLVfv6xVngGssdd5njoXmj28/6M5dhdZMtYXRy4wC/C3ywPs+Vx7eXKt9f9g79cGat011RI/E0A6sMw+1v/BuluoTR9n4FFgvR3X21h3grWp4wy8h3WNqRKrp3pTcxzXht6jsZc+ya+UUsot9BSZUkopt9AGRimllFtoA6OUUsottIFRSinlFtrAKKWUcgttYJRqJmKNbnyHPR0nIrM9HZNS7qS3KSvVTOyx3z41xvT1cChKNQufE1dRSjWRJ4DuIpKN9bBaqjGmr4hMBcZjjZ3VE2vARj/gOqACuNgYUyQi3bEemosCDgE3G2PWN/9uKOUaPUWmVPOZDmw2xqQD99VZ1he4AhgCPA4cMtYglT9gDdcBVg71XxtjBgO/BV5slqiVOkXag1GqZVhgjDkAHBCRUuATu3w10N8e/foMYJYcSxzp3/xhKuU6bWCUahkqak07as07sP6femHlNUlv7sCUOlV6ikyp5nMAK031STNWzp6tInIVOHOtD2jK4JRqatrAKNVMjDGFwGIRWQM8fQqbuAa4SURWAmuxkmsp1WLpbcpKKaXcQnswSiml3EIbGKWUUm6hDYxSSim30AZGKaWUW2gDo5RSyi20gVFKKeUW2sAopZRyC21glFJKucX/A7COP0QTbBUJAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%matplotlib inline\n",
"fig, ax = plt.subplots()\n",
"files = [\"test_jac.bin\", \"test_jacinv.bin\", \"test_dh.bin\", \"test_dhinv.bin\"]\n",
"earth = [1,2,1,2]\n",
"for i,f in enumerate(files):\n",
" sa = rebound.SimulationArchive(f)\n",
" data = np.zeros((len(sa),2))\n",
" for j, sim in enumerate(sa):\n",
" data[j][1] = sim.particles[earth[i]].calculate_orbit(primary=sim.particles[0]).e\n",
" data[j][0] = sim.t\n",
" ax.plot(data[:,0],data[:,1], label=f)\n",
"\n",
"ax.set_xlabel(\"time\")\n",
"ax.set_ylabel(\"eccentricity\")\n",
"ax.legend();"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# all lines above are on top of each other except the test_jacinv.bin one"
]
}
],
"metadata": {
"anaconda-cloud": {},
"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.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 1
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment