HyperbolicOrbits.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Loading Hyperbolic Orbits into REBOUND\n",
    "\n",
    "Imagine we have a table of orbital elements for comets (kindly provided by Toni Engelhardt).  Let's first load it and display it.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table><tr><td>tracklet_id</td><td>e</td><td>q[AU]</td><td>i[deg]</td><td>Omega[deg]</td><td>argperi[deg]</td><td>t_peri[MJD, days]</td><td>epoch_of_elements[MJD, days]</td><td>epoch_of_observation[MJD, days]</td></tr><tr><td>168026</td><td>12.181214</td><td>15.346358</td><td>136.78247</td><td>37.581438</td><td>268.412314</td><td>54776.806093</td><td>53371.0</td><td>55516.41727</td></tr><tr><td>21170</td><td>2.662235</td><td>2.013923</td><td>140.646538</td><td>23.02949</td><td>46.292039</td><td>54336.126288</td><td>53371.0</td><td>53673.44043</td></tr><tr><td>189298</td><td>15.503013</td><td>11.550314</td><td>20.042232</td><td>203.240743</td><td>150.855761</td><td>55761.641176</td><td>53371.0</td><td>55718.447145</td></tr><tr><td>72278</td><td>34.638392</td><td>24.742323</td><td>157.984412</td><td>126.43154</td><td>178.612758</td><td>54382.158401</td><td>53371.0</td><td>54347.240445</td></tr><tr><td>109766</td><td>8.832472</td><td>9.900228</td><td>144.857801</td><td>243.102255</td><td>271.345342</td><td>55627.501618</td><td>53371.0</td><td>54748.37722</td></tr></table>"
      ],
      "text/plain": [
       "[['tracklet_id',\n",
       "  'e',\n",
       "  'q[AU]',\n",
       "  'i[deg]',\n",
       "  'Omega[deg]',\n",
       "  'argperi[deg]',\n",
       "  't_peri[MJD, days]',\n",
       "  'epoch_of_elements[MJD, days]',\n",
       "  'epoch_of_observation[MJD, days]'],\n",
       " (168026, 12.181214, 15.346358, 136.78247, 37.581438, 268.412314, 54776.806093, 53371.0, 55516.41727),\n",
       " (21170, 2.662235, 2.013923, 140.646538, 23.02949, 46.292039, 54336.126288, 53371.0, 53673.44043),\n",
       " (189298, 15.503013, 11.550314, 20.042232, 203.240743, 150.855761, 55761.641176, 53371.0, 55718.447145),\n",
       " (72278, 34.638392, 24.742323, 157.984412, 126.43154, 178.612758, 54382.158401, 53371.0, 54347.240445),\n",
       " (109766, 8.832472, 9.900228, 144.857801, 243.102255, 271.345342, 55627.501618, 53371.0, 54748.37722)]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "class ListTable(list):     # this is just to nicely display the table (credit: Caleb Madrigal)\n",
    "    def _repr_html_(self):\n",
    "        html = [\"<table>\"]\n",
    "        for row in self:\n",
    "            html.append(\"<tr>\")\n",
    "            \n",
    "            for col in row:\n",
    "                html.append(\"<td>{0}</td>\".format(col))\n",
    "            \n",
    "            html.append(\"</tr>\")\n",
    "        html.append(\"</table>\")\n",
    "        return ''.join(html)\n",
    "\n",
    "table = ListTable()\n",
    "# this is what each column represents\n",
    "table.append(['tracklet_id', 'e', 'q[AU]', 'i[deg]', 'Omega[deg]', 'argperi[deg]', 't_peri[MJD, days]', 'epoch_of_elements[MJD, days]', 'epoch_of_observation[MJD, days]'])\n",
    "\n",
    "import numpy as np\n",
    "comets = np.genfromtxt(\"comets.txt\", skip_header=1, dtype=None) # load the table\n",
    "for comet in comets:\n",
    "    table.append(comet)\n",
    "table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now want to add these to simulation(s).  The first thing to do is set the units, which have to be consistent throughout.  Here we have a table in AU and days, so we'll use the gaussian gravitational constant (AU, days, solar masses).  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import rebound\n",
    "sim = rebound.Simulation()\n",
    "k = 0.01720209895 # Gaussian constant\n",
    "sim.G = k**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We would then add the giant planets to the simulation.  You could for example query JPL HORIZONS for the states of the planets at each comet's corresponding epoch of observation (see [Horizons.ipynb](Horizons.ipynb)).  Here we set up toy masses and orbits for Jupiter & Saturn:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sim.add(m=1.) # Sun\n",
    "sim.add(m=1.e-3, a=5.) # Jupiter\n",
    "sim.add(m=3.e-4, a=10.) # Saturn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Let's write a function that takes a comet from the table and returns an orbit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def getOrbit(comet_elem):\n",
    "    o = rebound.Orbit()\n",
    "    o.e = comet_elem[1]\n",
    "    o.a = comet_elem[2]/(1.-o.e) # q = a(1-e)\n",
    "    o.inc = comet_elem[3]*np.pi/180. # have to convert to radians\n",
    "    o.Omega = comet_elem[4]*np.pi/180.\n",
    "    o.omega = comet_elem[5]*np.pi/180.\n",
    "    o.T = comet_elem[6] # time of pericenter passage\n",
    "    return o"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now set up the first comet.  We begin by setting the simulation time to the epoch of observation (when our orbital elements are valid):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "comet_elem = comets[0]\n",
    "sim.t = comet_elem[-1] # last column is the time of observation\n",
    "o = getOrbit(comet_elem)\n",
    "sim.add(a=o.a, e=o.e, inc=o.inc, Omega=o.Omega, omega=o.omega, T=o.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default, REBOUND adds and outputs particles in jacobi orbital elements.  Typically orbital elements for comets are heliocentric.  Mixing the two will give you relative errors in elements, positions etc. of order the mass ratio of Jupiter to the Sun ($\\sim 0.001$).  If that matters to you, you should additionally pass `primary=sim.particles[0]` in the `add` call above, and you'd output with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<rebound.Orbit instance, a=-1.3693531710375346 e=12.201677899305905 inc=2.386856395628858 Omega=0.6563769038524557 omega=-1.5984922294899382 f=0.662947294342688>\n"
     ]
    }
   ],
   "source": [
    "comet = sim.particles[3]\n",
    "print(comet.calculate_orbit(primary=sim.particles[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You'll notice that, as expected, these elements differ from the ones in the table above by $\\sim 0.001$ because we added in Jacobi, and outputted heliocentric.  If you're consistent, they'll be the same.  We can quickly plot what the system looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU0AAAFCCAYAAAB1po8RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4lNXd//H3yb5MVgJhC2uAKAoimwLaiCLiUrRad5G2\n0lqXqi0WsVrA5bHaxx0ViwpaHwRUhEaLgEhEFAIKYQk7gYQQErJMJpN1tvP7g8AvQIBMSHLfM/m+\nrmsukslM8mGSfHLOvZxbaa0RQgjROAFGBxBCCF8ipSmEEF6Q0hRCCC9IaQohhBekNIUQwgtSmkII\n4YUgowOcC6WUHC8lhGgRWmvV0P0+P9LUWrf6bdq0aYZ8XV/KZLY8Zsx0tjyVlZVkZWWRnZ2Ny+Uy\nRSYzvk4tcTsTny9NIfxRQUEBX3/9NR06dKBnz54EBgYaHUnU8enpuRD+xu12k5OTQ01NDVdffTUW\ni8XoSOIkUppNkJqaanSEU5gtk9nygPkynZynsrKS/fv3Ex0dTUpKCgEBrT8RNNtrBObLpM42fzcz\npZT25fxCHLN9+3asVivnn38+cXFxRsdp85RS6NPsCJKRphAGcrlcHDhwgJqaGgYOHCjTcR8gI00h\nDGK329m/fz/t2rWjc+fOKNXgwEYYQEaaQpiIx+Nh79691NTU0KNHD6Kjo42OJLwghxwJ0YqcTic7\nd+4kLy+PlJQUKUwfJNNzIVqJzWYjJyeH9u3b07FjR5mOm5hMz4UwkMfj4dChQ5SVldGrVy/Z2ePj\npDSFaEG1tbWsXbuWyMhIBg0aRFCQ/Mr5OpmeC9FCrFYrubm5dOjQgU6dOhkdR3hBpudCtCKPx0Ne\nXh7l5eUkJycTGRlpdCTRjGTvuRDNyG638+233+J0OjnvvPOkMP2QTM+FaCalpaXk5uZisVhITk42\nOo44BzI9F6IFeTweDh48SEVFBf369SM8PNzoSKIFyUhTiHNQWlrKwYMHiY6Opnv37oasTCSa35lG\nmvIdFqKJiouL2bhxIzExMfTs2VMKs42QkaYQXnK73eTm5lJdXU2vXr0ICwszOpJoZmcaaUppCnEG\nbrebpUuXsmnTJgYNGsQvfvGL4zt7kpKSZHTpp2RHkBBN4Ha7GTt2LBkZGVRWVhIeHk7v3r1ZuXIl\n7du3NzqeMIj8mRTiNJYuXUpGRgYVFRVoramqqiI7O5uMjAyjowkDSWkKcRqbNm2isrLyhPuqqqrI\nzMw0KJEwA8NLUyn1vlKqUCm1pd5905RSeUqpjXW3a4zMKNoerTUrV6485RrYkZGRXHTRRQalEmZg\neGkCc4CxDdz/itb64rrb160dSrRdNpuN22+/ndzcXC666CIsFgtKKSwWC8OHD2fcuHFGRxQGMnxH\nkNZ6jVKqewMfkhVaRavbuHEjt99+O2PGjGH79u0EBwezdOlSMjMzueiiixg3bhyBgYFGxxQGMsUh\nR3Wlmaa1HlD3/jRgImADfgL+orW2NfA8OeRINAuPx8Nzzz3HW2+9xZtvvsmtt95qdCRhINMfp9lA\nabYHirXWWin1HNBJa/27Bp4npSnOWXl5OXfffTerV69m/fr19O3b1+hIwmA+d5ym1rqo3ruzgbTT\nPXb69OnH305NTSU1NbXFcgn/k5mZyZ133smVV17J4cOHZbGNNio9PZ309PRGPdYsI80eHB1pXlj3\nfketdUHd248BQ7XWdzbwPBlpiibRWjN79mxmzJjBq6++KtNxcQJTjzSVUvOAVKCdUioXmAZcoZS6\nCPAAB4A/GBZQ+J2KigqGDx+O2+0mPT2dPn36GB1J+BBTjDSbSkaawlvbtm3j7rvvpn///rz55pvE\nx8cbHUmYkCwNJ9o8rTVz587l2muvZerUqfzf//2fFKZoEsOn50K0tCNHjnDvvfdSVFTEN998I3vH\nxTmR0hR+befOndx2220EBASwZs0audCZOGcyPRd+a/78+Vx77bU8/vjjbNq0SQpTNAsZaQq/U15e\nzlNPPcUPP/xAWloa/fv3NzqS8CMy0hR+JTs7mz59+rBv3z6+++47KUzR7KQ0hd9IS0tj3LhxPPnk\nk6SlpWGxWIyOJPyQTM+Fz3M4HEyfPp2lS5eyYMECWe9StCgZaQqf9vPPPzNkyBBycnJYtWqVFKZo\ncVKawmd9++233HXXXVxwwQV8/PHHxMbGGh1JtAEyPRc+x+12889//pP58+czZ84cLr30UqMjiTZE\nSlP4lAMHDjB58mQ8Hg8rVqyQS+mKVifTc+EzNmzYwODBg+nUqROffvqpFKYwhIw0hekdW/vyrbfe\n4uOPP5YLmwlDSWkKU7Pb7fz5z38mLy+PL7/8kqSkJKMjiTZOpufCtJYtW8agQYNo164dixcvlsIU\npiClKUzp888/Z/Lkydx666384x//IDQ01OhIQgAyPRcmU1tby/Tp01m3bh2fffYZ/fr1MzqSECeQ\nkaYwja1bt3LzzTdjtVpJS0uTwhSmJKUpTGH16tVcccUVDB06lHfeeUcW2xCmJdNzYSiPx8PMmTOZ\nP38+aWlpcnaPMD0pTWGYsrIyHnvsMaqrq1m8eDEdOnQwOpIQZyXTc2GIr776iuHDh9O3b18+/vhj\nKUzhM+S656LVffrpp7z44ovcfPPNTJ061eg4QpziTNc9l+m5aDW1tbU888wzbNq0iYULF9KrVy+j\nIwnhNZmei1axbds27rjjDioqKvj000+lMIXPktIULe7HH3/kqquuYujQobz22mtyKV3h02R6LlqM\nx+Nh9uzZfPLJJyxevJhLLrnE6EhCnDPDS1Mp9T5wPVCotR5Qd18csADoDhwAbtVa2wwLKbxmt9uZ\nOnUqVquVBQsWkJiYaHQkIZqFGabnc4CxJ933BPCN1rof8C0gu1h9yKpVq7jsssvo2LEjc+fOlcIU\nfsUUhxwppboDafVGmjuBX2itC5VSHYF0rXVKA8+TQ45MZtmyZcyYMYOrr76a6dOnGx1HiCbxxUOO\nOmitCwG01gVKKTny2eTcbjevvfYaK1eu5IMPPiAl5ZS/ca2itrYWp9NJQEAAbrcbt9uNzWYjMjKS\noKCjP+5aa2pra4mJiSEoKIigoCCUavD3Q4hTmLU0T3ba4WT90UxqaiqpqamtEEfUl5eXx7Rp0wgJ\nCeGTTz4hJiamRb+e0+mksrISu91OdnY2HTp0wOFw4HA4sFqtKKXo1KkTAQEBBAQEUFxcTLt27QgO\nDgaOFuuRI0dITEzE5XLhdrspKysjLCyMTp06ERISQmhoKG63m/j4+ONlK/xXeno66enpjXqsWafn\nO4DUetPzVVrr8xp4nkzPDZaVlcW4ceO49957mTFjBgEBzbuZ3OPxYLfbKS8vZ/fu3WitiY6OJiIi\ngvDwcKqrq+natSshISGEhIQ0ueBqa2txOBy43W5qa2upra0lOzub+Ph4QkNDiYiIIDg4mNDQUOLi\n4pr1/yjMxxem56rudsx/gInAi8C9wBIDMomzWLx4MW+++SazZ89m7NiT9+U1XXV1Nfn5+ezcuZO4\nuDgiIiKIjo4mJSWF2NjYFjnOMzQ09JTV4bt164bWmpqaGqqrqykuLqasrIzo6GgsFgsWi4Xw8HBZ\nVb6NMXykqZSaB6QC7YBCYBqwGPgUSAJyOHrIUVkDz5WRpgEcDgf/+7//y08//cQrr7xCjx49zvlz\nFhQUUFZWRkVFBS6X6/i0uFu3bqabHjscDioqKqioqGDXrl306tWL+Ph4oqKiTJdVNM2ZRpqGl+a5\nkNJsfTt37uThhx9mwIABPPPMM+c06nM6nRQWFlJQUEBhYSGdOnUiOTmZqKgon9kx4/F4qKysxGaz\nUVFRgdaauLg4EhMTfeb/IE4lpSmaxebNm/nLX/5Cr169mDVrVpO3X9psNr7//nuioqJo3749nTp1\nIjY21udLxuPxcOTIEex2O0opYmNjiY2NldGnD5LSFOdEa80XX3zBv/71L/7+978zYsSIJn2OkpIS\ncnNzqampISoqir59+/rt9sDa2lrKysooKipCa+3X/1d/JKUpmqympoaXX36ZrKwsXnrpJbp27erV\n87XWHDx4kDVr1pCcnEz37t3p0KGDz48qG8vlclFUVER1dTXh4eHH98YLc5PSFE1y5MgRrr76agYP\nHszMmTMJDw/36vklJSXs3bsXrTWJiYn07NmzhZKan8fjoby8nOLiYqqrq0lJSTl+3KgwHylN4bUt\nW7bw1FNPcc0113D//fd7tf2yuLiYDRs2EB0dTe/evWWnSD1ut5vCwkJqa2uJiooiLi6OwMBAo2OJ\nk0hpikbTWrN48WLee+89nn76aa+Wc3O73ezbt4/c3FxiYmIYPHhwsx/s7i9cLhelpaVUVVURFhZG\nx44djY4k6pHSFI1SVVXFlClTKC0t5R//+AdJSUmNfm5+fj579+4lJiaGlJQU2W7XSJWVlWRnZ9O+\nfXsSEhJkT7tJSGmKsyopKeHPf/4zBw8eJC0trdHHX7pcLjIzM9m0aRM33XQTCQkJLZzU/2itsdls\n2Gw24uLiiI6ONjpSmyelKc5o+/btPP3004wfP56777670VPqsrIytmzZQlxcnOzYaAYOh4PDhw9T\nUVFBv379ZNRpIClNcVpLly5l1qxZTJkypdHHX3o8HtatW0d1dTX9+/eX7XHNSGtNUVERDoeDhIQE\nwsLCjI7UJklpilO4XC4mTpxISUkJs2bNonv37o1+XmZmJvv27ePaa68lKiqqhZO2TdXV1ZSUlGCx\nWIiNjTU6TpsjpSlOYLPZmDFjBpWVlcyYMaPRI8XKykp+/vlnYmNj6d+/vxwq08JcLhf79+8nKCiI\nHj16yGFbrUhKUxy3b98+nn76aa644gp+97vfNXr75a5du9i5cycDBw5sllWNRON4PB5KSkrweDwk\nJCTIH6pWIqUpAFiyZAlz587l4YcfZvTo0Y1+3uHDh1mzZg2DBg0iOTm5BROK07HZbFRVVdGuXTtC\nQkKMjuP3pDTbOK0177zzDq+99hpLlizhvPNOWQT/tHJzc9m1axfDhg1r8ctYiDMrLy/nwIED9OvX\nT46DbWG+sHK7aCFVVVW8+OKL2O12MjIyvLpUQ0ZGBna7nUsvvRSLxdKCKUVjHDst1Wq1Eh8fLyNO\ng8g5bn4sPz+fRx55hLi4OF566SWvCnPXrl1s27aNYcOGSWGaSGRkJLGxsZSWluJwOIyO0ybJ9NxP\nzZ8/nzlz5vDQQw9x/fXXe7Xndf/+/ezbt4+RI0d6vbKRaB01NTXk5eXRsWNH+aPWAs40PZeRpp/R\nWrNkyRI+/PBDJk2axA033OB1Ye7Zs4dLL71UCtPEwsLCaN++PXa7HZfLZXScNkW2afoRp9PJ22+/\nzc6dO5kzZ47XZ+rs37+fxYsXM2nSpBa54qNoXjExMQQHB1NaWkq7du3kcKRWItNzP2G1Wnn++eeJ\niYlh8uTJXo8Sy8vLWbNmDYMHDyYxMbGFUoqWcOzKmLJKUvOR6bmf27NnD6NGjaJPnz489dRTXhdm\nTU0N69at48ILL5TC9EEWiwWPx0NxcbHRUdoEGWn6uLVr1/LWW29xxx13cN1113n9fIfDwX/+8x/6\n9+/v1fGbwly01litVkJCQmTHUDOQ4zT9kNaahQsXsmzZMqZPn97kM3U2b96M2+0mJSWlmROK1qSU\nIiYmBqvVSnBwsBz83oKkNH1QeXk5L7zwAk6nk5dfftmr4y/rO3z4MCUlJdx0002yGIQfCAwMJDo6\nmoKCAhITE6U4W4iUpo8pKSnhL3/5C7W1tXz44YdNPiukurqan3/+mUsuuUTOLPEjISEhREVFUVFR\nIaXZQmRHkA/ZvXs3U6ZMYfz48cybN6/JZefxeJg3bx5JSUlyeQo/dOwKl1VVVUZH8UumHmkqpQ4A\nNsADOLXWw4xNZJzvv/+e999/nwcffJChQ4ee0+fKzs4mLi6OCy+8sJnSCbOJioo6vmNIDkNqXmZ/\nNT1AqtbaanQQo2it+fOf/0xRURHPP/98o1dYP52qqiqysrK46qqr5PK6fiwwMJCIiAgOHz7s1VVF\nxdmZvTQVbXgTQk1NDTNnzqSyspJ//OMfdO3a9Zw/58aNG+nTp49c8bANiIiIoLq6mpqaGrnWUDMy\neyFpYIVSaoNSapLRYVpTcXExTz/9NJGRkbz11lvNUpjr169n165dcnhRGxITE0N1dTUej8foKH7D\n7CPNkVrrw0qp9hwtzx1a6zVGh2ppWVlZvPHGG1x//fVer1B0Om63m4KCAkaPHu2303KXy4XD4aCo\nqIiwsDBsNhsAAQEBVFRU0K5dO4qKikhMTMTpdBIWFkZ4eDiBgYF+e0B4UFAQISEhVFVV+e3/sbWZ\nujS11ofr/i1SSn0BDANOKM3p06cffzs1NZXU1NRWTNj8li1bxl//+ldmz57NsGHNt99rz549xMfH\n07Nnz2b7nEY7dupgeXk5DoeD/Px8evbsSVVVFSEhIURHRxMcHExwcDAJCQmEhIQQHh5OSEgIFRUV\neDwerFYrdrudyMhIKisrSUpKIjg42K8WLDm2bTMoKEim6aeRnp5Oenp6ox5r2tMolVIRQIDWukIp\nFQksB2ZorZfXe4zfnEaptebTTz/lhx9+4MEHH6Rv377N9rlra2tZunQpo0eP9vltmR6Ph4qKCvbu\n3YvT6cRisRAeHk5iYiJhYWFNXunH5XJRXl6O1pq8vDxiYmIIDw8nKiqKiIiIZv5ftL7Kyko8Ho9c\ncrmRfPU0ykTgC6WU5mjO/6tfmP6ktraWd955h7KyMp599tlmL7a0tDRiY2N9ujCdTidWq5XNmzfT\nsWNHEhMTSUhIaLYDuIOCgoiPjwegXbt2OBwOjhw5QkFBAaGhoVgsFp++RlJkZCQ2mw2XyyWHIJ0j\n0440G8MfRprZ2dm8+uqrnH/++dx3330EBwc36+evqanh008/5YYbbiA2NrZZP3dr8Hg8HDlyhJ9/\n/pm+ffuSmJjYquWvtaa0tJSKigqCgoJ8ujxra2txOBwy2mwEuRqlSe3fv5+pU6fSv39/nnrqqRY5\n/3vz5s04nU6GDBnS7J+7pZWUlLBu3Tr69OlD586dDd+RYbPZsFqtREZGHt8s4GuKi4uJjIz0yeyt\nSUrThDZs2MDcuXP57W9/y+DBg1vkazgcDtLS0hg7dqzhheON2tpaMjIyiImJISkp6fi02Qy01pSU\nlFBeXk6HDh2IjIz0qcVOqqqqcLlcPr2ppjVIaZqI1pq0tDS+/fZbHn30UXr06NFiX+ubb77B4/Fw\n9dVXt9jXaG4HDx6kqKiI8PBw+vTpY9rtb263m/z8fAA6d+7sU5easNvtRERE+FTm1uarO4L8jsvl\n4v777ycmJoZp06Y1eUm3xvB4POTk5DBmzJgW+xrNSWvNoUOHyMrKYvDgwbRv397oSGcUGBhIUlIS\n5eXl5OfnExcX5zOj+eDgYBwOh0zRm0hKs5XY7XZmzpxJx44deeKJJ1r8Fyw/P58uXbrQrVu3Fv06\nzcHj8bBmzRpCQkIYPXq0Ty1VFx0dTUBAAJWVlQQEBPjE4UkhISGUlZURGhrqtyc6tCQpzVaQn5/P\n66+/zvDhw1ttwd+9e/fSp0+fFv8658rhcJCVlUV8fDwpKSmmnY6ficViISwsjIMHD9K+fXvTjzgD\nAgIIDAyktrZWRptNIH9mWtimTZt48cUXGT9+PL/61a9apTBLSkrIy8s75xWRWprH4yEzM5PKykr6\n9+/vk4V5TFBQEElJSdTW1h4/fdPMIiMj8bX9AWYhpdmCvvzySx588EEeeOABRowY0Wpfd8uWLQQE\nBJi6hFwuF+np6URFRTFy5Eif2gN9OkFBQURHR1NVVUVtba3Rcc4oKCgIrbUs5NEEUpotwOPxsGDB\nAtavX8/ChQvp169fq379qqoqRo8e3apf01s5OTmEhobSp08fvyjMY4KDg+nQoQOlpaWmXzk9KCgI\nh8NhdAyfY96hiI+qqalh9uzZOJ1OnnjiiVbfMVBRUUFVVRUdOnRo1a/rjW3btpGfn8/o0aNNPRpu\nqsDAQGJjY4+vY2nWnS0BAQGUl5fLIh5eMud300dlZ2fz1FNPERsbyyOPPGLIntScnBySkpJM+4tq\nt9ux2WwMHTrULwvzmGOrKRUVFRkd5bSCg4OJiIiQKbqXzPmb5YNycnJ4/vnn6dy5M3fffbdhBw5/\n/fXXpjqDpj6Px8Py5cvp3Llzix6jahbHTlesrq42OsppBQQE4Ha7jY7hU/z3T30ryszMZN68efzp\nT39i4MCBhuVwOBzExsaads3MAwcO0Lt3b9PmawkWi4XCwkKCgoKafTGW5hAYGIjL5TJlNrOSkeY5\n0FqzYsUKFi5cyMMPP2xoYQIUFhbSrVs3U/4C2Gw2fvzxx2ZdJ9QXBAQEEB8fb9odLsd2BskUvfFk\npNlEbrebyZMnExYWxpQpU0yx7Nrhw4fp2LGj0TEalJOTw8iRI33ijJnmFhoait1uJyAgwJQHkwcF\nBeHxeEy7Hdxs5FVqgurqat555x3i4+OZOnWqKQoTjl5bKDEx0egYpzh8+DBZWVk+cUpnS4mJiTHt\naM6XTls1AxlpeqmkpIRZs2aRkpLCAw88YJq/ztXV1eTk5NCuXTujo5yiqKiIESNGtOlVdYKDg6mo\nqDDlaFMpJSu6e8Ecv/E+YteuXbz66qtcdtll3HzzzaYpTDi6zXDIkCGmK6YjR46wbdu2ZrkEsa8z\n6znpAQEBckqlF8zzW29y69at4/777+fXv/41l19+udFxTlFcXExCQoLRMU5RVlbGoEGDTFfmRggO\nDqa2than02l0lBMopXA6nbhcLqOj+AQpzbPQWrNy5UqWL1/OBx98YPge8tMpLS013fGZDoeDDRs2\ntOltmScLDw835aguPDzcr05nbUlSmmfgdrtZuHAhmZmZPProo6Y+vvCHH35otiszNher1UqXLl38\n6hri5yo4OJiamhqjYzTIjGVuRrLl9zRKS0t5++236dSpEw899JDpCulkUVFRpttzvnPnTrp06WJ0\nDFMJCAggODjYdDte6i7vYHQMnyAjzQaUlpby8ssv43K5+M1vfmP6wqyqqiI2NtZ0Cy+Ulpaa/rIV\nRggKCjLdqYsBAQGyTbORzPOnziRyc3P56KOPuO6661p1DcxzYbfbTXct65KSEiwWi2mOYTUTpZTp\njtmU7ZmNJ6VZz9atW1m0aBG33nor5513ntFxGq2kpMR0o+FjV2oUpzo20tRam6qszLS5wMzkVeLo\nBvBZs2axe/duJk+e7HPb4fbt20d5ebnRMU7R1s4z94bH48HtdpuuqMxW5GbU5rdpejwelixZQlFR\nEY899pjPFSZAYmIi/fv3NzrGCTIzM42OYGohISGmOjlCNJ65/sy1straWubNm4dSiilTpphuittY\nDofDdGebdOzY0dSrxxvt2PV5zFScxzLJiQhnZp7vWAOUUtcopXYqpXYrpaY05+cuLi7m3XffJS4u\njgkTJvhsYcLR8jfTogs1NTXk5eWZ7hxrM1FKmaow4egedLNlMiPTvkJKqQBgJjAW6A/coZRKaY7P\nnZ+fz+9//3t69+7N+PHjff4HxeFwmKr0nU6n6c5OstvtDB06lIqKCqOj4Ha7+e9//8tzzz3Hl19+\naYrDj9xuN1999RXPPvusaTKZltbalDfgEmBpvfefAKac9BjtrV27dukXXnhBb9y40evnmtXu3bv1\nG2+8YXSM48rKyvTq1auNjnGC+fPna0AvWLDA0Bwul0tfeeWV2mKxaKWUtlgs+sorr9Qul0symUhd\ntzTYTWcdYimlHlZKGXFBly7AwXrv59Xd12Tr169n8eLF3HXXXQwaNOicwplJnz59TDU9DwwMJDs7\n2+gYAEyaNIno6GjuueceAO6++26io6P5/e9/b0iepUuXkpGRQUVFBVprKioqyMjIYOnSpYbkMWsm\nM2vMjqBEYINSaiPwAbCsrolNYfr06cffTk1NJTU19ZTHeDwe/v3vf1NWVsZ9991nuqmjv4mMjDTN\nQe1PPvkkGRkZ7N27F6fTSVBQED169GDq1KmG5Nm0aROVlZUn3FdZWUlmZibXX3+9ZDJIeno66enp\njXvw6Yag+sRpsOLotsX5wF7gf4DejXluU28cnZ5/Xe/9Jk3PnU6n/vjjj/UDDzygKyoqmjBQ9w2z\nZs0yOsJxJSUleubMmUbHOG7RokU6KChIWywWHRQUpBctWmRYlrS0NG2xWDRw/GaxWHRaWppkMhHO\nZXp+rJmAgrqbC4gDPlNKvdS4am6SDUCyUqq7UioEuB34jzefoKqqio8++ojw8HDeeOMNv11tZ/ny\n5abYwXFMUFCQqc6o+uyzz4iKiuL555/HYrHw+eefG5Zl3LhxDB8+/PiRBRaLheHDhzNu3DjDM4WE\nhKCUMkUmMzvr9Fwp9QgwASgG3gMe11o76/Zu7wH+2hLBtNZupdRDwHKO7uV/X2u9o7HPLykpYd68\neVxwwQWkpqb69VkObrfbdAvb7t+/3+gIxz3xxBO8/vrrJCQkcOedd1JQUGBYlsDAQJYtW8asWbNI\nS0vjoYceYty4cYYeG3ks08SJE3G5XNx1112GZzKzxmzTjAd+pbXOqX+n1tqjlGrRDR5a66+Bft4+\n76effuLDDz/kt7/9rV/t8DmdkJAQbrzxRqNjHBcREUFkZCRut9sUv3gXXnjh8bcTEhIMX+E+MDCQ\nbt26cdVVV5lmm2FgYCChoaFcfvnlpslkVmednmutp51cmPU+1uiRX2vZvn073377Lb/+9a/bRGHC\n0UvE1tbWGh3juKCgIAIDA7HZbEZHMa2qqirDy/tkTqfTNDvwzMyvTqNcu3YtP//8MxMnTmxTp/B5\nPB5KS0uNjnGCwMBArFarHKlwGkeOHDHdH3WtNeeff77RMUzPt0+FqeNyufjqq6/Iysri3nvvbVOF\nCUfPdsnLyzM6xgmSk5Ox2+1GxzCt3Nxc0+2YLC0tNeUloM3G50vT4XDwwgsvsGXLFiZMmGC6xXhb\nQ//+/U03rXK73ezatcvoGKZUXl6OUspUF+lzu91YrVYpzUbw+en5vHnzGDRoENdcc43p1iZsLVFR\nUaZbT/OtB7MHAAAWNUlEQVT8888nNzfXNDuDzOTHH3803YIdOTk5OJ1OgoODjY5ieub5rjVRSkoK\n119/fZstTIDo6Ojjp8CZRWhoKEVFRaY69MgsbDYbY8aMMTrGCaxW6wlHGYjT8/nSvOSSS4yOYLjg\n4GAOHjxISUmJ0VFOMGTIEFMddG8GbrebH3/8keTkZKOjnCAnJ4devXoZHcMn+HxpiqOSkpJMd4hP\n165dycjIMN1FxIy0b98+2rdvT/fu3Y2OcoIDBw6YLpNZSWn6ifPOOw+r1Wp0jBPEx8eTkJBg6Bk4\nZpOWlsaFF15oujPUMjMz29xRJ00lpekn2rdvT3FxsdExThAQEEBKSgrr1q0zOoopVFZWUlBQwLBh\nw4yOcgKtNQcOHPDJ62MZQUrTT7Rv355Dhw4ZHeMUPXv2xGazyTGbwKJFi0hJSaFTp05GRznBoUOH\niI+PN93F+cxKStNPxMTEsH79etPteImIiGDgwIF8//33RkcxVE1NDdu2bTPdKBNg27ZtXHDBBUbH\n8BlSmn4iODiYa665xpTbD5OTk8nLyzPdsaSt6fPPPyclJcWUh/Vs3brVlLnMSkrTjyQlJZnudEo4\nehzpiBEjWLlypdFRDGGz2cjMzGTEiBFGR2nQf//7XxITE42O4TOkNP1IUlISBw8ePPsDDdC3b1+K\niorIzc01OkqrmzNnDqmpqfTr5/Uqhy2usLCQkJAQLr30UqOj+AwpTT/Svn17du7cSVFRkdFRThES\nEsL48eP55ptvTLfdtSVt3ryZwsJCLr/8cqOjNGjt2rWMGDHCVJeANjspTT8SEBBA7969OXDggNFR\nGpSYmEh8fDzr1q0z1SmfLeXw4cN88skn/OEPfzDtQjLHSlM0npSmn7nsssvIyWlwzWhTuOGGGygu\nLmb16tVGR2lRbrebzz//nL59+9KjRw+j4zTIbreTkZHB8OHDjY7iU6Q0/UyPHj0oLCykurra6CgN\nCgwMZNSoURw4cIDDhw8bHadFaK2ZOXMmABMnTjQ2zBmsWrWKuLg4oqOjjY7iU6Q0/UxwcDBdunRh\n8+bNRkc5ra5du3LVVVfx8ccfm/IQqXP13XffUV1dzZ133mmq5d9OlpGRwX333Wd0DJ9j3u+oaLKI\niIjGX/jeIF26dGHs2LGsXr3ar4pzxYoVfPnll/zxj3809aU+bDYbmzZtIjU11egoPkdK0w9deeWV\nREdHm24Bj5MNGDCArl27smDBAsrKyoyOc86WLFnCpk2beOyxx4iJiTE6zhktX76c4cOHm+6SG75A\nStMPBQYG0r9/f7Zs2WJ0lLMaMWIE1113HQsWLGDHDtNd3LRRPB4PK1asICMjg9tuu830C19orfnX\nv/5Fnz59jI7ik6Q0/dSAAQPIysrC6XQaHeWskpOTGTx4MOvXr/eJoq+vtraWV155hczMTJ588kmf\nWJMyKyuLhIQEbrnlFqOj+CQpTT8VFxdHSUmJz5y6OGTIEMaPH8+yZctYtmwZDofD6EhnlZWVxeuv\nv05SUhIPPvggFovF6EiN8tlnn3HzzTcTEhJidBSfJKXpx2666Sby8/N95kDy2NhY/vSnPxEWFsbr\nr79Odna20ZEa5HA4WL58OZ999hkpKSncdtttREREGB2rUQoKCti4cSPXXnut0VF8lvKVX6iGKKW0\nL+dvaVprPvzwQ0aNGmW6a9KczebNm8nKysLlcnH99debZk/07t27+fzzzwkICOD+++83/Q6fk/3x\nj38kNjaWF154wegopqaUQmvd4PL6pixNpdQ0YBJwpO6uJ7XWXzfwOCnNs9i5cyc//fQTd911l+ku\nsXA2NTU1fPXVV9jtdqKjoxk9erRh13c/cOAAS5YsobS0lDvuuIO+ffua+hjMhpSUlDBhwgRee+01\n2Ql0Fr5amnat9StneZyU5ll4PB5mzJjBuHHjfPbKnXa7nS+++AKn04nFYuGSSy6hW7duLf5HwOFw\nsGPHDlasWEFVVRXjxo1jwIABPru4xZtvvklgYCAPPPCA0VFM70ylaeaLhfvWsMikAgICuPLKK9my\nZQvDhw/3udEmQFRUFBMmTKC6upqvvvqKVatWUVlZSUpKCoMGDWrWqbvb7eann36iqKiIdevWkZiY\nyJgxYzjvvPN8esdJYWEhK1euZO7cuUZH8XlmHmlOBGzAT8BftNanXJ9WRpqNo7Xm3//+N4MHD/aL\n68Bordm8eTNHjhwhKysLt9vNBRdcgFKKgQMHEhcX16jRoNYaq9XKoUOHsFqtZGRkEBgYSHV1NcOG\nDePiiy+mXbt2rfA/anmTJk2iT58+/PWvfzU6ik8w5fRcKbUCqL9ctAI08DdgHVCstdZKqeeATlrr\n3zXwOaQ0GykvL4+0tDQmTZpEUJCZJxje8Xg8HDlyhNzcXLZu3Up4eDhbt24lKSkJl8tFZGQkFosF\nm81Gly5d2LVrF506dSInJwePx0NQUBBKKS6++GLcbjcjR470u7Nkdu3axeTJk3n33Xfp3Lmz0XF8\ngilLs7GUUt2BNK31gAY+pqdNm3b8/dTUVDmX9gzee+89wsLCuPvuu42O0qJcLhfV1dUcPnwYp9OJ\n2+2mtLSUjh07cvjwYXr27IlSivj4eNOuc9lcPB4PjzzyCNdddx3XXHON0XFMKz09/YT1GmbMmOFb\npamU6qi1Lqh7+zFgqNb6zgYeJyNNL+Tk5PDRRx/x4IMPmuYQHtGyli5dytKlS3nttdd8bm+/kXxu\npKmU+gi4CPAAB4A/aK0LG3iclKaXMjIy2LNnj08egiS8k5eXx4QJE3j77bdJSUkxOo5POVNpmvJP\nj9Z6gtZ6gNb6Iq31jQ0VpmiaoUOH4nQ6Tb3epjh3Wmvee+89Ro0aJYXZzExZmqLlBAQEMG7cOGbP\nnu23K6cL+Pbbbzl06BBPPfWU0VH8jpRmG9SxY0fGjBnDsmXL8Hg8RscRzayoqIh3332XKVOm+PSx\npWYlpdlGjR8/nqCgIL7//nujo4hm5HK5uO+++xgxYoTPrTfgK6Q02yilFDfeeCObNm1iz549RscR\nzeTDDz+kZ8+e/PGPfzQ6it+S0mzDLBYLV1xxBS+88AI22yknXAkfs2bNGlavXs2MGTN89vx4XyCl\n2cYNHDiQCRMmsGDBAp9Y+Fc0bNOmTbzxxhv87W9/87nl6nyNlKbgF7/4BR07dmTx4sU+s2Cx+P/K\nysr4+9//zqhRo+jbt6/RcfyeKQ9ubyw5uL35uFwu5s6di1KK3/3ulNP8hUlVV1czZcoUhg4dyj33\n3GN0HL/hcwe3i9YXFBTE+PHj2blzJ2vWrDE6jmgEl8vF//zP/9C7d2+/X0/ATGSkKU5gs9mYM2cO\nI0eOZOjQoUbHEafhdru57777iI2N5Z///KdfrVxlBjLSFI0WExPDhAkT+P7771m3bp3RcUQDXC4X\nr776KtHR0UybNk0Ks5VJaYpTxMfHc+edd/L222+zdu1ao+OIelwuFy+//DI2m42XXnrJsGsmtWUy\nPRendeTIET7++GOGDRvGqFGjjI7T5lVXV/Poo4/SsWNH/va3v8kpki3I55aGaywpzZZns9n46KOP\n6Nu3L2PGjJE1GQ1it9t57rnnKC4u5o033vC71eXNRkpTnJPKykqeffZZEhISePTRR2UbWisrKCjg\nmWeeYciQIUycOFH+cLUCKU1xzux2O2lpadhsNu666y6io6ONjtQmLF++nA8++IB77rmH6667zug4\nbYaUpmgWWmtWr17NunXruOWWW+jdu7fRkfyW1pqvvvqKOXPmcPPNN3Pnnadc7UW0IClN0ay2b9/O\ns88+y3333cfo0aPlshnNrLKykpkzZ1JUVMTjjz9OYmLi2Z8kmpWUpmh2BQUFLF68mJCQEG655RaZ\nrjeTtWvX8s9//pOxY8fym9/8RvaQG0RKU7QIt9tNeno669atY8iQIYwdO9boSD7L4XCwcOFCvvnm\nGy677DI5/99gUpqiRW3ZsoW33nqLUaNGccMNN8gB117avHkz77//PklJSUyaNElePxOQ0hQtzuVy\nsWrVKtauXcvo0aO55JJL5NCks7DZbMydO5clS5bwzDPPyAkEJiKlKVpNUVERH330ETt27ODxxx+n\nb9++sqPoJC6Xi6+//ppFixZx+eWXc+ONN8ro0mSkNEWr8ng8rF27lh9++IGYmBiuvfZakpKSjI5l\nOI/Hw8qVK3nzzTcZNmwYEydOpGvXrkbHEg2Q0hSGcLvdbNiwga+//hq73c4jjzxCt27djI7V6jwe\nDz/++COLFi0iPDycoUOH8stf/lLO7DExKU1hKLvdzpdffkl2djZxcXFcccUVpKSk+P20vaamhsWL\nF/PFF1/Qt29fbrnlFgYMGOD3/29/IKUpTMHtdrN582ZWrVpFVlYWt912G6NGjfK7xSfy8/NZsWIF\na9asoXfv3qSkpDB+/HgpSx8ipSlMxePx8N1333Ho0CGysrLo168fAwcO5IILLiA4ONjoeE1SUVHB\n+vXrmT9/PlarlZtvvpmrrrqKhIQEo6OJJjBlaSqlbgGmA+cBQ7XWG+t9bCrwW8AFPKK1Xn6azyGl\n6eOqq6vZuHEjCxYsoKamhtTUVC688EL69etn+rNhrFYrixYtoqSkhN27dzNgwAAGDhzIkCFD/G70\n3NaYtTT7AR7gXWDysdJUSp0HzAOGAl2Bb4A+DbWjlKZ/KS4uJisri23btpGdnU1tbS2/+tWvSE5O\npkuXLgQGBhqar7q6mqysLLKzs9m6dSuFhYV4PB5uv/12Ro4cSXh4uKH5RPMxZWkeD6DUKuAv9Urz\nCUBrrV+se38pMF1rndHAc6U0/VRpaSlr1qyhurqaffv2YbVaCQ4OJjk5mYsuuojOnTvTvn37FitS\nh8PB1q1bqaysZP/+/WRnZ7Njxw4iIiK4/vrrueCCC0hOTpYD+P3UmUrTjN/xLkD9C9McqrtPtCHx\n8fH88pe/PP5+RUUF3333HeXl5fz000/k5+djtVopKCggNTWVuLg44uLiCAgIICwsjC5duhAaGkpo\naChaa0JCQtBa4/F4cDgclJaWEhwcjN1ux263s2PHDiorK/F4PBQWFmKz2SgsLOSqq66iX79+XHrp\npXTv3p3Q0FADXxVhBi1amkqpFUD9da0UoIG/aa3TmuNrTJ8+/fjbqamppKamNsenFSZjsVhOWYS3\nsrKSbdu2ERoaitVqxWq1sn37dsrKykhISKC2tpaamhqysrLo3LkzFouFgIAAnE4nubm5DBkyBIvF\nQlRUFE6nk8TERC666CISExOJj483fHOAaD3p6emkp6c36rG+MD3/Gpgm03MhRGvxheue1w/3H+B2\npVSIUqonkAysNyaWEEKcyLDSVErdqJQ6CFwCfFm3wwet9XZgIbAd+C/wgAwnhRBmYfj0/FzI9FwI\n0RJ8YXouhBA+QUpTCCG8IKUphBBekNIUQggvSGkKIYQXpDSFEMILUppCCOEFKU0hhPCClKYQQnhB\nSlMIIbwgpSmEEF6Q0hRCCC9IaQohhBekNIUQwgtSmkII4QUpTSGE8IKUphBCeEFKUwghvCClKYQQ\nXpDSFEIIL0hpCiGEF6Q0hRDCC1KaQgjhBSlNIYTwgpSmEEJ4QUpTCCG8IKUphBBeMKw0lVK3KKW2\nKaXcSqmL693fXSlVpZTaWHd726iMQghxsiADv/ZW4Cbg3QY+tldrfXED9wshhKEMK02t9S4ApZRq\n4MMN3SeEEIYz6zbNHnVT81VKqVFGhxFCiGNadKSplFoBJNa/C9DA37TWaad5Wj7QTWttrdvWuVgp\ndb7WuqIlswohRGO0aGlqrcc04TlOwFr39kal1D6gL7CxocdPnz79+NupqamkpqY2JaoQog1LT08n\nPT29UY9VWuuWTXO2AEqtAiZrrX+uez8BKNVae5RSvYDvgAu11mUNPFcbnV8I4X+UUmitG9y3YuQh\nRzcqpQ4ClwBfKqWW1n3ocmCLUmojsBD4Q0OFKYQQRjB8pHkuZKQphGgJphxpCiGEL5LSFEIIL0hp\nCiGEF6Q0hRDCC1KaQgjhBSlNIYTwgpSmEEJ4QUpTCCG8IKUphBBekNIUQggvSGkKIYQXpDSFEMIL\nUppCCOEFKU0hhPCClKYQQnhBSlMIIbwgpSmEEF6Q0hRCCC9IaQohhBekNIUQwgtSmkII4QUpTSGE\n8IKUZhOkp6cbHeEUZstktjxgvkxmywOSqTGkNJvAbN9EMF8ms+UB82UyWx6QTI0hpSmEEF6Q0hRC\nCC8orbXRGZpMKeW74YUQpqa1Vg3d79OlKYQQrU2m50II4QUpTSGE8IKUpheUUrcopbYppdxKqYvr\n3d9dKVWllNpYd3vbyDx1H5uqlNqjlNqhlLq6NfI0kG+aUiqv3utyjUE5rlFK7VRK7VZKTTEiw8mU\nUgeUUpuVUpuUUusNyvC+UqpQKbWl3n1xSqnlSqldSqllSqkYg/OY4meoPilN72wFbgK+a+Bje7XW\nF9fdHjAyj1LqPOBW4DxgHPC2UqrBjdqt4JV6r8vXrf3FlVIBwExgLNAfuEMpldLaORrgAVK11oO0\n1sMMyjCHo69LfU8A32it+wHfAlMNzgMG/wydTErTC1rrXVrrPUBDBdTqpXSGPOOB+Vprl9b6ALAH\nMOoX06iyPmYYsEdrnaO1dgLzOfr6GE1h8O+f1noNYD3p7vHAh3VvfwjcaHAeMP5n6ARSms2nR930\nYZVSapTBWboAB+u9f6juPiM8pJTKVEq915pTvXpOfi3yMO61qE8DK5RSG5RSk4wOU08HrXUhgNa6\nAOhgcB4w/mfoBEFGBzAbpdQKILH+XRz9Af+b1jrtNE/LB7ppra112xYXK6XO11pXGJSn1ZwpH/A2\n8IzWWiulngNeAX7X+ilNaaTW+rBSqj1Hy3NH3UjLbIw+JtF0P0NSmifRWo9pwnOc1E0rtNYblVL7\ngL7ARiPycHRkmVTv/a519zU7L/LNBowo+UNAt3rvt9hr4Q2t9eG6f4uUUl9wdDOCGUqzUCmVqLUu\nVEp1BI4YGUZrXVTvXaN+hk4g0/OmO76dRSmVULfDAaVULyAZyDYqD/Af4HalVIhSqmddnlbfQ1v3\nS3fMr4BtrZ0B2AAk1x3hEALcztHXxzBKqQillKXu7Ujgaox5beDoz83JPzsT696+F1hiZB6T/Ayd\nSGstt0beOLpR/CBQDRwGltbdf+ybuRH4CbjWyDx1H5sK7AV2AFcb9Hp9BGwBMoHFQKJBOa4BdnF0\nh9gTJvg56ln3mmzi6BEQhmQC5nF001ItkAv8BogDvql7vZYDsQbnMcXPUP2bnEYphBBekOm5EEJ4\nQUpTCCG8IKUphBBekNIUQggvSGkKIYQXpDSFEMILUppCCOEFKU0hhPCClKbwa0qpIXWL/YYopSLr\nFm0+3+hcwnfJGUHC7ymlngHC624HtdYvGhxJ+DApTeH3lFLBHF24oxoYoeWHXpwDmZ6LtiABsABR\nQJjBWYSPk5Gm8HtKqSXAJxxdXaiz1vphgyMJHyaLEAu/ppS6B3BorefXrXn6g1IqVWudbnA04aNk\npCmEEF6QbZpCCOEFKU0hhPCClKYQQnhBSlMIIbwgpSmEEF6Q0hRCCC9IaQohhBekNIUQwgv/D6nl\nHp2ofPPZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1080a5c88>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "fig = rebound.OrbitPlot(sim, trails=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we just integrate until whatever final time we’re interested in:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU0AAAFCCAYAAAB1po8RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl0Tef+BvDnmwSRRGJIBRESU+tqTUVQNDVVWlSxXKUt\n1Utpq8O6pSjLcH9quGpVLx0UlaJmJe0SUpq41JQbGYoiQgwRQuaT+eS8vz+kVpQMO3Kyz/B81srq\nsc/O3t+cxuN93/3ud4tSCkREVDEOehdARGRNGJpERBowNImINGBoEhFpwNAkItKAoUlEpIGT3gU8\nChHhfCkiMgullDxsu9W3NJVSNvU1d+5c3Wuw1S9+tvxsK/pVFqsPTSKi6sTQJCLSgKFpYQICAvQu\nwWbxszUfe/pspbz+uyUTEWXN9RORZRIRKFu9EEREVJ0YmkREGjA0iYg0YGgSEWnA0CQi0oChSUSk\nAUOTiEgDhiYRkQYMTSKiElJSUsp836qXhiMiqiomkwnnzp1DQkJCmfvxNkoisnvZ2dm4fPky3Nzc\n4O3tjZo1a5Z6GyVbmkRkt5RSuHnzJpKTk9GsWTPUq1ev3O9haBKRXSooKMChQ4fQuHFjtG3bFjVr\n1qzQ97F7TkR2JzU1FdeuXYO7uzv8/PweeL+sVY7Y0iQiu1FUVISrV68iJycHrVu3houLi+ZjsKVJ\nRHbh1q1buHjxInx8fNC0aVM4OJQ+45ItTSKyW0opJCUl4fr16/Dx8UGzZs0e6XhsaRKRzcrPz8fl\ny5fh5OSE5s2bo0aNGhX6PrY0icjuJCQkICMjA40bN0bDhg2r7Li8jZKIbIrRaERcXBzOnTuHVq1a\nVWlgAuyeE5ENyczMREJCAurXrw9vb2+IPLSHXS52z4nIpplMJiQmJiItLQ2+vr5wd3c327kYmkRk\n1XJzc3Hs2DHUqVMHnTp1gpOTeWON3XMislrJyclISkqq8os9ZXXPGZpEZHUKCwtx5coVGI1G+Pn5\noVatWlV6/LJCk1fPiciqJCcnIzQ0FC4uLnj88cerPDDLw5YmEVkFk8mE69evIyMjA15eXlU+lagk\nds+JyKrl5OTg8uXLcHFxQbNmzeDo6GjW83HKERFZrdOnTyMvLw8tWrRA/fr19S6HoUlElqmwsBCX\nL1+GwWDAk08+CTc3N71LAsDuORFZoLS0NFy9ehUNGzZEo0aNKn1nT2Wxe05EVsFoNCI+Ph6FhYVo\n1aoVXF1d9S7pAQxNIrII2dnZOH36NJRS6Nq1q9kv9lQWu+dEpKvKPBHS3DjliIgsRlFREUJCQhAV\nFYV27drhiSeegJOTE3x9fSv8REhz45gmEVmEoqIiPP/88zhx4gSys7NRq1YtPP300zh06JDFdsf/\nirdRElG1CQkJwfHjx2EwGKCUQl5eHmJiYhASEqJ3aRXG0CSiarN7925kZ2ffty07OxvR0dE6VaQd\nu+dEZHaFhYUYPnw4fv31Vzg7OyMvL+/ee66urujYsaOO1WnDliYRmVVcXBz69OmDzMxMHDlyBM88\n8wzc3NwgInBzc4O/vz8CAwP1LrPCePWciMxCKYX169fj448/xpw5c/DOO+/AwcHh3tXz6OhodOzY\nEYGBgRZ3EYhTjoioWl2+fBkjR45EUVERNm7ciCeffFLvkjSx6EWIRSRBRGJEJEpEThZvqycioSJy\nXkT2i4iH3nUSUcWEhYUhICAALi4uOHr0qNUFZnl0D00AJgABSqlOSqluxdtmADiglHocwK8AZupW\nHRFVSH5+PmbOnInx48djzZo1OHz4MFxcXPQuq8pZQmgKHqzjJQBBxa+DAAyr1oqISJN9+/ahZ8+e\nuHDhAiIjIzFgwAC9SzIbSwhNBeAXEYkQkX8Ub/NSSt0CAKXUTQDmW9eeiCpNKYU1a9Zg5MiReP75\n57Fjxw54enrqXZZZWcI8zWeUUkki8hiAUBE5j7tBWlKpV3vmzZt373VAQAACAgLMUSMR/cWdO3cw\nefJk3LhxA5GRkXj88cf1LqnSwsPDER4eXqF9LerquYjMBWAA8A/cHee8JSKNAIQppdo+ZH9ePSfS\nwd69ezF16lSMHTsWs2fPtpiFNqqKxS7YISIuAByUUgYRcQUwEMB8AMEAxgNYAmAcgD26FUlE9+Tl\n5WHcuHHYt28ffvrpJ/Tp00fvkqqd3t1zLwA/iogqrmWTUipURP4HYJuITABwBcAoPYskIuDMmTOY\nMGECWrdujdjYWDRv3lzvknRhUd1zrdg9JzI/k8mEr7/+GkuXLsXChQsxZsyYan9mT3Wz2O45EVm2\nxMRE9O7dG97e3ggLC4Ofn5/eJenOEqYcEZEFCgkJQd++ffHss89i//79DMxibGkS0X1yc3PxySef\n4ODBgwgKCkL37t31LsmisKVJRPfs3LkT/v7+MBgMOHz4MAPzIRiaRASTyYRVq1Zh6tSpGDp0KFav\nXg13d3e9y7JI7J4T2bmkpCS8/fbbKCgowPHjx9GsWTO9S7JobGkS2bENGzbgueeeQ69evbBnzx4G\nZgWwpUlkh7Kzs/HBBx9g+/btOHDgALp06aJ3SVaDLU0iOxMdHY0BAwbA1dUVV65cYWBqxJYmkZ0w\nGo348ssvsXr1aixatAhDhgzRuySrxNAksgN/3tnj5+eH0NBQNGnSRO+SrBa750Q2Ljg4GIGBgXj9\n9dcREhLCwHxEbGkS2SiDwYBPPvkEkZGR2LBhAzp06KB3STaBLU0iG7Rt2zZ07doVNWvWxP79+xmY\nVYgtTSIbYjQasWrVKnzzzTcYM2YM5syZo3dJNoehSWQjrl+/jvfffx81a9bEr7/+ikaNGuldkk1i\n95zIBvznP//B0KFDMXDgQGzatImBaUZsaRJZsaysLMyaNQs7duzATz/9xInq1YAtTSIrderUKQwe\nPBgeHh64dOkSA7OasKVJZGUKCgrw5ZdfYtOmTVi4cCEGDhyod0l2haFJZEUSExMxcOBANG7cGD//\n/DO8vLz0LsnusHtOZCWCg4MxYsQITJ48Gfv27WNg6oQtTSILl5WVhfnz5+P06dNYt24d/va3v+ld\nkl1jS5PIgu3evRvdunWDq6srfvzxRwamBRCllN41VJqIKGuun6g0RUVF+Oqrr7Bx40aMHDkSH330\nkd4l2RURgVJKHvYeu+dEFiYxMRHTpk2Ds7MzgoOD0bBhQ71LohLYPSeyIJ9//jlGjRqFAQMGYM2a\nNQxMC8SWJpEFMBgM+Ne//oXt27dj8+bN8Pf317skKgVbmkQ6i4mJwciRI+Hq6oqzZ88yMC0cW5pE\nOjEajfj222+xZcsWzJ07F3379tW7JKoAhiaRDm7evIkhQ4agXr162L59O8curQi750TVLDQ0FK+8\n8gomTJiAvXv3MjCtDFuaRNUkJycHS5YswalTp7Bq1SpOVLdSbGkSVYPQ0FD07NkTIoKtW7cyMK0Y\n7wgiMiOTyYQNGzZg/fr1CAwMxPTp0/UuiSqAdwQR6SA5ORlz586F0WjE5s2b+QgKG8HuOZEZbNiw\nAa+++iqefvppfPPNNwxMG8KWJlEVys/Px/Lly/Htt99i48aN6Nmzp94lURVjS5OoisTHx2P8+PHI\nzs5GbGwsA9NGsaVJ9IhMJhN27NiBtWvX4r333sMLL7wAkYdeQyAbwNAkegTp6ekYMWIEatSogTVr\n1sDHx0fvksjM2D0nqqTIyEiMGzcOAwYMwK5duxiYdoItTSKNjEYj1qxZg71792L27Nno1q2b3iVR\nNWJLk0iDU6dOYfjw4YiLi0NQUBAD0w4xNIkq6MCBA/joo4/Qpk0bLFu2DPXq1dO7JNIBu+dE5cjN\nzcWKFSsQGxuLr7/+Gm3atNG7JNIRW5pEZTh9+jQmTpwIk8mEtWvXMjCJLU2ih1FK4YcffsDMmTOx\nevVqDBo0SO+SyEIwNIn+IjMzE0uWLEFKSgrCw8PRokULvUsiC8LuOVEJsbGxmDhxIpo2bYqVK1cy\nMOkBbGkS4e6tkO+99x6uXLmC2bNn84mQVCqGJtm91NRULFq0CElJSVi+fDlat26td0lkwdg9J7sW\nGRmJKVOmoEOHDti+fTsDk8rFlibZpfz8fAQFBSEsLAyzZs1Chw4d9C6JrARDk+zO7du3MW7cOLi4\nuGDt2rXw8PDQuySyInywGtmViIgILF++HIMHD8Yrr7wCBweOUNGD+GA1sntFRUX4/vvvcfDgQXzy\nySd48skn9S6JrBT/mSWbFxcXh7Fjx+LixYtYuXIlA5MeCbvnZNNiYmKwYMEC+Pr6YtmyZXwMBVVI\nWd1zhibZJKUUtm7dir179+Kjjz5C+/bt9S6JrAjHNMmuZGZm4rPPPkNeXh4+//xz1K9fX++SyIZw\nTJNsSmxsLHr37o2mTZvi008/ZWBSlWP3nGyCUgqhoaHYuHEjRowYgWHDhuldElkxds/JphUUFOCb\nb77B+fPnsXjxYnh7e+tdEtkwtjTJql24cAFfffUVmjRpgqlTp8LZ2RlKKcTGxuL48eNISEhAeno6\nCgoKICJwcXGBp6cnnnjiCfTv35/dd3ooXj0nmxQbG4s33ngDkyZNwqRJk3DhwgV8+eWX+P3332Ey\nmdC0aVP4+PigUaNGcHNzg8lkQkZGBm7duoWEhAQkJyejcePGGDhwIEaNGgU3Nze9fySyEAxNsilK\nKfz000/Ys2cPPvjgA7i5uWHOnDm4dOkSOnTogJEjRyIgIACOjo5lHsdgMCA4OBh79+7FnTt3MHTo\nUEyYMAHOzs7V9JOQpWJoks3Iz8/HV199hevXr2P69OlYt24dgoOD0aVLF8yYMQNNmjSp1HEjIiKw\nYsUK5OfnY/bs2Vz1yM4xNMkm/PHHH/jss8/g4OCAxx57DEeOHIGTkxPmzJmDgICARz6+yWTC999/\nj23btmHs2LEYO3bsoxdNVomhSVbv3Llz+Pjjj3H69GncunUL2dnZcHR0RM+ePREWFlZuV1yLU6dO\nYdGiRQgICMA777xTZccl61FWaHJyO1m8gwcPYtmyZejfvz+Sk5ORnZ0N4O7KRVFRUQgJCanS83Xu\n3BkLFy5EeHg4tm7dWqXHJuvH0CSLVVRUhO+++w7BwcFYsGAB0tPT7wXmn7KzsxEdHV3l527Tpg3m\nzJmDTZs2ISYmpsqPT9aLoUkW6fbt23jttdeQmJiITz/9FE2aNEH79u3h4uJy336urq7o2LGjWWpo\n3749xo4di2XLliE3N9cs5yDrwzFNsjg3b97E4sWL4eTkhEWLFqFGjRoA7i7EMWzYMERERCA7Oxuu\nrq7w9/fH/v37HxjTNJlMuHbtGlJSUpCXlwcAqFGjBtzc3ODh4YF69eqhdu3a5dailMK7776L9u3b\n46233qr6H5YsEi8EkdU4e/YsvvjiC4waNQp9+/a9tz0/Px95eXlwc3NDSEgIoqOj0bFjRwQGBj4Q\nmImJiYiNjYWLiwt8fHzQoEGDewGZm5uL3Nxc5OXlwcXFBR4eHvjll18QFRWFTp06PfR4Z8+exbx5\n87B69WrUrVvX/B8C6Y6hSVYhLCwMW7duxdSpU9GuXbv73ktPT0ft2rVRq1atMo8RGxuL+Ph4dO7c\nGc2bNy91P6UU0tLSMGzYMERFRZXbcv3nP/+JDh064PXXX6/8D0hWg1fPyaKZTCZ8+OGHCAoKwty5\ncx8ITKPRiPz8/HID848//sD58+cREBBQZmACd/9SHD16FFFRUTAYDFBKwWAw4MSJEw+9Gv/iiy/i\n6NGj2n84sjkMTdJVfn4+VqxYAQcHByxevBiNGzd+YJ+cnJx745qlycjIQEREBHr16oV69epV6Nx/\ntjBLKu1qfJ8+fZCamoqkpKQKHZtsl0WHpogMEpFzInJBRD7Wux6qWmlpaVi4cCFcXV2xZMkSNGrU\n6KH75ebmwsmp7FUMT548CW9v74eGbmk6deoEV1fX+7aVdjXeyckJPj4+iIiIqPDxyTZZbGiKiAOA\nlQCeB9AOwCsi8oS+VVFVuXLlCubPn48uXbpg4sSJZYZiQUFBmV3zwsJCJCQklHm/eFZWFrp27QqD\nwXBvW2BgIPz9/e+tblSjRg34+/sjMDDwocdo2bIl4uLiyvvRyMZZbGgC6AYgTil1RSlVCGALgJd0\nromqQFhYGMaMGYO///3vGDp0aLlPiMzNzS2ze56UlARHR0d4enqWus/evXvxv//9D3v37r23zdHR\nEfv378fmzZsxceJEDBgw4KEXgf7k7e2N9PT0cn46snXlhqaITBWRig0SVS1vANdK/Pl68TayYuHh\n4fjxxx+xbNky9OjRo9z9i4qKYDKZ4OBQ+q9qampqqXMuJ06cCHd3d7z22msAgFdffRXu7u6YNGkS\ngLvBOXjwYHzxxRdo1qwZCgoKSj2Ph4cHCgsLy62ZbFtFHnfhBSBCRE4BWAdgvyXN85k3b9691wEB\nAVWy2g1VPaUUduzYgcjISMyePRsNGzas8PcWFRWV+X5BQQFK+5WcNWsWTpw4gYsXL6KwsBBOTk7w\n9fXFzJkz79vvz+5/YWFhqQHs5OTE56bbqPDwcISHh1do33JDUyk1W0TmABgI4A0AK0VkG4C1Sqn4\nRym0HIkAmpX4c9PibfcpGZpkmXJzc7Fu3TpkZGRg9uzZmlZIFxEYjcYy96ldu/Z9Y5Ul+fn5Yf78\n+fdWZs/Ly8P8+fPh5+d3336JiYkQkTJry8/P5wLFNuqvDa758+eXum+FxjSLW5Y3i7+MAOoB2CEi\nSx+l0HJEAGglIs1FpCaA0QCCzXg+MgODwYBx48bh6tWrmDZtmuZHSvzZLS+r2+zj44P09PRSu847\nduxAnTp1sHDhQri5uWHnzp0P7PPbb7/B2dm53GEADw8PTfWT7Sm3pSki7wN4HcAdAGsATFNKFRZf\n3Y4DMN0chSmlikTkXQChuBvua5VSf5jjXGQet2/fxueff45hw4Zh9OjRZQZSWZycnJCbm4uaNWs+\n9P26devC3d0dJ06cQK9evR54f8aMGVixYgU8PT0xZswY3Lx584F9YmJi0KpVqzLrSEpKqvTK8GQ7\nKjKmWR/AcKXUlZIblVImERlsnrLunWMfgMfNeQ4yj8uXL2PlypUYMmTII48zu7m5IT09vcxWXqdO\nnXD48GF069btgXB96qmn7r329PR84Cr7b7/9hszMTIwcObLMOhITE9GnT59K/ARkS8r9p18pNfev\ngVniPbb86AGHDx/GF198gXHjxlXJhbk6deogNTW1zH26desGV1dXbN68GSaTqcLHTk5Oxo4dO+Dv\n71/mxambN28iMzPzgVs8yf5Y8jxNskL79+/HzJkz8fbbb6N9+/ZVcsy6desiJyfngVseSxIRvPLK\nK0hKSsL69evLHAP9U3x8PP7973/Dy8ur3OcBhYeHo23btuXezkm2j6FJVUIphV27duHo0aPYtWsX\nWrduXWXHdnBwgJeXF65ceWiH5566devirbfeQmZmJpYuXYqwsLCHhmdKSgpWr16NlStXws/PD9Om\nTStzvLWgoAARERH3LVVH9otLw9EjMxqNCAoKQmpqKt5+++0H7ueuCvn5+Thy5Ai6d+9e7vGLiooQ\nFhaG48ePIzs7G3Xq1IGrqyuMRiMyMjJgMBhQr149DB48GJ06dSr33Hv27MHvv/+O2bNnV9WPQxaO\n62mS2aSkpGDVqlXw8vLChAkTzNp9vXjxIu7cuYPu3btX+HsSEhJw7tw5pKWl3Xv0b4cOHdCgQYMK\nff/t27exePFivPXWW2jTpk1lSycrU1ZoVuTqOdFDpaWlYe7cufDw8MDEiRMrPaWoolq0aIHExERc\nvHix3OlBf/L19YWvr2+lzmc0GrFu3Tp07tyZgUn3MDSpUm7cuIFVq1bh5ZdfRr9+/arlnA4ODujc\nuTP++9//wtXVVdMycJXx/fffQ0QwatQos56HrAsvBJFm58+fx4oVKzB8+PBqC8w/1alTB127dkVk\nZCRu3LhhtvP88MMPSEpKwuTJk3nFnO7DMU3SZPPmzYiIiMCUKVOq9Aq5VklJSTh58iTatGmDtm3b\nVtlxc3NzsXHjRty5cwdTpkzhg9TsFC8E0SNTSiEsLAxr167F9OnTy1zwt7qkp6fj8OHDcHFxuW8x\n4co6c+YMdu7ciaZNm+LVV18t9bZNsn0MTXokSin8+OOPOH/+PCZPnlzhZ/BUB6PRiKioKJw/fx4t\nWrRAu3btNC+qceHCBRw8eBB37typ8DQksm0MTaq0wsJCbNy4EZmZmZg0aVKpa03qzWAw4NSpU4iL\ni0ODBg3QtGlTNG7cGI899tgDLcasrCxcv34d8fHxOHv2LIC7t2H27NmTrUsCwNCkSkpLS8OMGTPQ\nrVs3vP7661ZxQaSwsBDx8fG4evXqvfvFHR0d4eDggKKiIuTm5sJkMsHDwwONGjXCU089BV9fXy4u\nTPdhaJJm6enp+PrrrwEA06dPN/scTHMxmUzIzs5Gfn4+HBwcULt2bYttLZPlYGiSJklJSfj222/R\nu3dvBAQEsBVGdod3BFGFxcfHIygoCMOGDUPnzp31LofI4jA06Z6ff/4ZmzZtwoIFC3Sdg0lkydg9\nJwB3Vy8PDQ3FsGHDOOWG7B7HNKlUSimEhITg9OnTmDhxokXNwSTSC0OTHqqwsBA7duxASkoK3njj\nDbOsg0lkjXghiB6Qn5+PadOmoUGDBpg1a5ZVzMEksgRsadohg8GA9evXo379+hg1ahScnPhvJ1FJ\nbGnSPSkpKfjuu+/QsWNH9OvXj3MwiTRiaNqR8+fPY+fOnejXrx/8/f31LofIKlnnvXGk2YULFzB9\n+nQGJtEj4pimHYiOjkZISAhGjx4NPz8/vcshsngc07RTSikcOXIEJ06cwJtvvomGDRvqXRKR1WNL\n00aZTCYEBQXBYDBg3LhxcHd317skIqtRVkuTY5o2yGg0YuPGjTh58iQmTJjAwCSqQmxp2pj8/Hxs\n2bIFzs7OGDFiBOdgElUCxzTtRFZWFjZt2gQfHx8EBgZa7cLBRJaMoWkj4uLi8N1332HQoEHo3bs3\nJ60TmQlD0wYkJiZi8+bN6NChA/r06aN3OUQ2jWOaVu7ixYvYs2cPhg4dyoWDiaoIxzRtVExMDA4e\nPIjRo0fD29tb73KI7AJD00qtXbsWCQkJeP/99+Hp6al3OUR2g6FpZZRSOHDgADIyMjBlyhQGJlE1\n45imFSkqKkJwcDCys7MxcuRIODs7610SkU3i4y5sQF5eHnbu3InatWtj6NChnLROZEa8jdLKGQwG\nfPjhh3BycsLLL7/MwCTSEVuaFi41NRXbt29Hy5Yt0b9/f73LIbILnHJkpZKSkrBr1y707t0b7du3\n17scIgJD02KdOXMGv/76K1544QW0bNlS73KIqBjHNC1QTEwMli5diqFDhzIwiSwMxzQtTEREBKKj\no/Hyyy9zDiaRTjjlyAoopXDo0CFcuXIFw4cPR506dfQuichuMTQtXEFBAbZs2YJatWrhpZde4qR1\nIp1xnqYFKygowNatW5GQkIARI0YwMIksHFuaOsrJycGePXvg5eWF5557jgsHE1kIztO0QOnp6di9\nezfatm0Lf39/vcshogpiaOrg7Nmz2L17N4YMGYKnnnpK73KISAOOaVaza9euITQ0FD179mRgElkh\njmlWo7i4OBw+fBiDBg1CkyZN9C6HiErBMU0LEBMTg+joaAwbNgz169fXuxwiqiSGZjUICgpCWloa\n3nzzTU5aJ7JyDE0zMplMCA8Ph6OjI8aPH8/AJLIBHNM0E6PRiNDQUIgIBgwYwIWDiawIxzSrWXZ2\nNvbv34/69eujT58+cHDgJAUiW8G/zVUsKysLCxYsgJubGwICAhiYRDaG3fMqlJqain379qFt27bo\n1KmT3uUQUSWxe14NkpKScODAATzzzDNo0aKF3uUQkZkwNKtAbGwsYmJi0K9fP05aJ7JxHHB7RGfO\nnMGGDRvQv39/BiaRHeCY5iOIiorCpUuXMHDgQM7BJLIhHNOsYiaTCceOHUNKSgpefPFFLhxMZEfY\n0tSooKAA27dvh6enJ/r27YsaNWpU6/mJyPzY0qwiBQUF2LdvH3JyctC/f384OjrqXRIRVTO2NCso\nJycHBw4cQJMmTfD000/z0RRENowtzUeUkZGBAwcO4IknnkC7du30LoeIdMQpR+WIi4vDmjVr0LFj\nRwYmETE0y5KYmIjIyEj07dsXLVu21LscIrIAHNMsxaVLlxAdHY1nn30WDRo0MMs5iMgycUxTozNn\nzuDixYvo378/3N3d9S6HiCwIQ7MEpRSCg4ORlJSEcePGoXbt2nqXREQWht3zYiaTCSdOnEBWVhae\neeYZuLq6Vslxicj6lNU9Z2gCKCwsxG+//YYaNWqge/funLROZOc4plmGjIwMHDt2DJ6enpy0TkTl\nsuspRwaDAUFBQXB2dkaXLl0YmERULt1CU0Tmish1ETlV/DWoxHszRSRORP4QkYHmOH96ejoOHTqE\nwYMHIyAgwBynICIbpHf3fLlSannJDSLSFsAoAG0BNAVwQERaV+WEzOTkZJw8eRKdO3fmwsFEpIne\nofmw/vBLALYopYwAEkQkDkA3ACeq4oRRUVG4du0aevbsCU9Pz6o4JBHZEb3HNN8VkWgRWSMiHsXb\nvAFcK7FPYvG2RxYfH49jx46hS5cuDEwiqhSztjRF5BcAXiU3AVAAPgHwJYAFSiklIv8H4DMA/9B6\njnnz5t17HRAQUOr45NmzZ3Hjxg2MHz8eLi4uWk9DRDYsPDwc4eHhFdrXIuZpikhzAD8ppdqLyAwA\nSim1pPi9fQDmKqUe6J5XZJ6myWRCdHQ0srKy0KNHD9SsWdMcPwIR2ZCy5mnqefW8UYk/Dgdwuvh1\nMIDRIlJTRPwAtAJwsjLnKCgowLZt25CTk4NevXoxMInokel5IWipiHQEYAKQAOAtAFBKnRWRbQDO\nAigE8HZlrpwXFBQgIiIC7u7u6NGjB+/yIaIqYRHd88oqrXuem5uLEydOwMvLC23bttWhMiKyZnZ1\nG2VWVhZOnjyJFi1awM/PT+9yiMjG6D3lqErduHEDQUFBePzxxxmYRGQWNtM9T05Oxu+//46WLVvC\n19dX38LdYglGAAAEW0lEQVSIyKrZfPf82rVriIuLQ5cuXeDh4VH+NxARVZLVh+a5c+dw8+ZN+Pv7\nc+FgIjI7qx/TPH78OHr06MHAJKJqYfVjmvn5+Zy0TkRVio+7ICLSwCJvoyQiskYMTSIiDRiaREQa\nMDSJiDRgaBIRacDQJCLSgKFJRKQBQ5OISAOGJhGRBgxNIiINGJpERBowNImINGBoEhFpwNAkItKA\noUlEpAFDk4hIA4YmEZEGDE0iIg0YmkREGjA0iYg0YGgSEWnA0CQi0oChSUSkAUOTiEgDhiYRkQYM\nTSIiDRiaREQaMDSJiDRgaBIRacDQJCLSgKFJRKQBQ5OISAOGJhGRBgxNIiINGJpERBowNImINGBo\nWpjw8HC9S7BZ/GzNx54+W4amhbGnX77qxs/WfOzps2VoEhFpwNAkItJAlFJ611BpImK9xRORRVNK\nycO2W3VoEhFVN3bPiYg0YGgSEWnA0LQAIjJXRK6LyKnir0El3pspInEi8oeIDNSzTmslIoNE5JyI\nXBCRj/Wux9qJSIKIxIhIlIicLN5WT0RCReS8iOwXEQ+96zQXjmlaABGZCyBLKbX8L9vbAvgBQFcA\nTQEcANBa8X9ahYmIA4ALAPoBuAEgAsBopdQ5XQuzYiJyCcDTSqm0EtuWAEhRSi0t/oepnlJqhm5F\nmhFbmpbjYVfqXgKwRSllVEolAIgD0K1aq7J+3QDEKaWuKKUKAWzB3c+VKk/wYHa8BCCo+HUQgGHV\nWlE1YmhajndFJFpE1pTo2ngDuFZin8TibVRxf/0Mr4Of4aNSAH4RkQgR+UfxNi+l1C0AUErdBNBQ\nt+rMzEnvAuyFiPwCwKvkJtz95fsEwJcAFiillIj8H4DPAPzjwaMQWYRnlFJJIvIYgFAROY+7v8sl\n2ewQEkOzmiilBlRw128B/FT8OhGAT4n3mhZvo4pLBNCsxJ/5GT4ipVRS8X9vi8hu3B0CuSUiXkqp\nWyLSCECyrkWaEbvnFqD4l+xPwwGcLn4dDGC0iNQUET8ArQCcrO76rFwEgFYi0lxEagIYjbufK1WC\niLiIiFvxa1cAAwH8jruf6fji3cYB2KNLgdWALU3LsFREOgIwAUgA8BYAKKXOisg2AGcBFAJ4m1fO\ntVFKFYnIuwBCcbeRsFYp9YfOZVkzLwA/Ft/C7ARgk1IqVET+B2CbiEwAcAXAKD2LNCdOOSIi0oDd\ncyIiDRiaREQaMDSJiDRgaBIRacDQJCLSgKFJRKQBQ5OISAOGJhGRBgxNsmki0qV4wdyaIuIqIqdF\n5G9610XWi3cEkc0TkQUAahd/XVNKLdG5JLJiDE2yeSJSA3cX7sgF0JP379OjYPec7IEnADcAdQA4\n61wLWTm2NMnmicgeAJsB+AFoopSaqnNJZMW4NBzZNBF5DUCBUmpL8UPWfhORAKVUuM6lkZViS5OI\nSAOOaRIRacDQJCLSgKFJRKQBQ5OISAOGJhGRBgxNIiINGJpERBowNImINPh/k7IfAY74qWUAAAAA\nSUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1088e8470>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "tfinal = 60000 # Final time in days since J2000\n",
    "sim.integrate(tfinal)\n",
    "fig = rebound.OrbitPlot(sim, trails=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<rebound.Orbit instance, a=-1.3725067517418916 e=12.181305454717414 inc=2.387303295568815 Omega=0.6559203102943116 omega=-1.598504146646605 f=1.445539842442101>\n"
     ]
    }
   ],
   "source": [
    "print(comet.orbit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The osculating elements change slightly because of the presence of the giant planets.  For fun, let's add them all to a simulation (wouldn't actually make sense because the elements are for different observation epochs!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU0AAAFCCAYAAAB1po8RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XlUHNeZP/xv0QvQNPtOI5AAsWtBkoUky1qi2JIc2Uoc\nx46cTCbjxE7s5HXseJLYTmbsmUnsOJMzOc7Ec5yMfk6cjGRb3lckWQuyZAu0gYRYhQABTe9N73vV\nff9AECSxNFDdVd3czzk6BzVF1dMN/fStuzyXIYSAoiiKCk6M0AFQFEVFEpo0KYqiZoAmTYqiqBmg\nSZOiKGoGaNKkKIqaAZo0KYqiZkAqdABzwTAMnS9FUVRIEEKYiR6P+JYmISQs/55++umwXUuIf9H+\n/ObDcwzF83O73ejo6EBPTw9YlhXsubW0tKCvry9sv8OpRHRLk6Ko0CCEwGq1wmq1Ij8/HwkJCYLF\n4vV64ff7sWDBAsFiGI8mTYqirhEIBNDR0YGUlBSoVCpIpcKmicHBQWRmZiImRhw3xuKIIgJs2rRJ\n6BBCKtqfHxD9z5GP5+d0OqFWq5GZmYm8vDzBEybLsjh58iRSUlIAiON3yEx3/y5mDMOQSI6fosQi\nEAjAYDAgEAggMzMTcXFxQocEANBoNDCZTKiurg7rdRmGAZlkIIjenlPUPOfz+TAwMACWZVFSUiKa\n22AAUKvVKCkpETqMa4jn1aEoKuxsNhs0Gg2ysrJQWloqqoR55coVGI3GsVtzsaAtTYqah1iWxaVL\nlxAXFweVSgWZTCZ0SDcYGhqCSqUSOowb0KRJUfOM2+2GwWBASkoKsrKyRNW6HOV2u+F2u7F69Wqh\nQ7kBTZoUNU9wHAe9Xg+fz4fMzEzEx8cLHdKkent7kZeXB4lEInQoNxDfRwxFUbzz+/0YGBiA2WyG\nSqUSdcJkWRZ1dXXIysoSOpQJ0ZYmRUU5h8MBk8mElJQUFBYWCh3OtNRqNdauXYu0tDShQ5kQTZoU\nFaU4jsPly5cRExOD/Px8xMbGCh1SUHp6elBRUSF0GJOit+cUFYW8Xi/UajUUCgUKCwsjKmGOToES\nK9rSpKgoo9fr4Xa7kZ6eDqVSKXQ4M9Ld3Y2CggIwzISLcUSBtjQpKkqwLIuhoSFotVrk5uZGXMJ0\nu91wOByinGY0Hm1pUlQUGJ17qVQqsWTJElG31CbT2dmJwsJCUU60H4+2NCkqghFCcOXKFVy5cgWZ\nmZlIS0uLyITJsixeffVVUa4Auh5taVJUhPL7/TAYDJBKpSgqKoJcLhc6pFnr6+vDhg0bkJOTI3Qo\n06JJk6IikMlkgsPhQEpKCpKTk4UOZ87a29uxdOlSocMICr09p6gIwnEcdDod+vv7kZWVFZKEGe4a\ntV1dXRgYGBDNdhbTETRpMgyTzzDMEYZhWhmGaWEY5pGrj6cyDHOQYZhOhmEOMAwT+R+lFDVHXq8X\nQ0NDkEgkWLZsWciWQg4NDcFisYTk3BM5f/686KcZjSd0SzMA4MeEkCoAawH8gGGYcgBPADhECCkD\ncATAkwLGSFGC02g0uHz5MlJSUpCRkRHSykQsy4at8pHL5YLNZsMXv/jFsFyPD4ImTUKIlhDSfPVr\nB4B2APkAdgJ45ephrwD4sjARUpSwWJaFTqeD3+/HwoULwzL30mAwgOO4kF8HAFpbW1FcXBwxK5YA\n4VuaYxiGWQhgOYAGANmEEB0wklgBiHdNFUWFiNVqhVqthkwmw4IFC6BQKMJyXYlEEpa5kizL4k9/\n+hOKi4tDfi0+iSJpMgyjBPAmgB9dbXFe3xNNd0+j5g1CCMxmM7q7u5GWlhb2uZcxMTFhSZrd3d3Y\nsGFDxAwAjRJ8yhHDMFKMJMy/EULeu/qwjmGYbEKIjmGYHAD6yX7+mWeeGft606ZNotjik6JG+Xw+\nyGSyoJPe6K6QDMNg2bJlgmyhq9Vqw7L747lz57BmzZqQXycY9fX1qK+vD+pYwbfwZRjmrwCMhJAf\nj3vseQBmQsjzDMP8DEAqIeSJCX6WbuFLiZrH4wl6O1yTyQStVguVSiXYZmKjeweVl5eH9DodHR14\n77338JOf/ESU221MtYWv0FOObgbwDQBfYBimiWGYcwzDbAPwPIBbGYbpBLAFwK+FjJOiZoNlWXi9\n3mmP4zgORqMRNpsN+fn5gu6+6PF4oNdPemPHm88//xwqlUqUCXM6gt6eE0I+AzDZJiCRMweBoibA\ncdy0t+VutxtmsxkymQyFhYWiSCIZGRkhPb/L5YLZbMbDDz8c0uuEivC/IYqKUj6fb8qpO3a7Ha2t\nrYiPjxfNrpA+ny/kcTQ1NaGsrCxsswH4JvxviaKiFMuyE7Y0R3eFtNlsWLJkiaj2wvF6vfD5fCE7\nP8uyeOGFF1BVVRWya4QaTZoUFSJer/eGLWjtdju6uroQExOD3Nxc0U3qlslkyMzMDNn5Ozo6cPPN\nN6OoqChk1wg1mjQpKkRcLtc1LU2r1Qq9Xo/09PSQL4WcLaPRCKvVGrLznzx5Elu2bAnZ+cNBfL81\niooSHMdBLpfD7/dDp9PB6XSioKAgpC25uVIqlSHrLmhvb0d/f7+od5oMBk2aFBUCHMfBYDDA7/fj\n4sWLAIDc3FzRb+Vgs9lCVhqurq4OWVlZN3RZRBqaNCkqBLxeL+RyOYxGI0pLS5GdnR0Rpc84jgtJ\nt4HT6YTZbMa3v/1t3s8dbjRpUhTPPB4PLl68CL1ej7y8PCQkJAgdUtBiY2NDEm9jYyOqq6sjbofM\nidCkSVE8crlcGBoaQlxcHKqqqgRZOz4X/f39vJeFCwQC+N3vfodVq1bxel6h0KRJUTwYrUxkMpmQ\nl5cHlmXhcDiEDmvG8vLygl4rH6yWlhYsXboUJSUlvJ5XKDRpUtQcjQ72OByOsaSTmJgY8uWIfOM4\nDlqtlveBmvr6etx11128nlNINGlS1Bw4nU4MDQ1BpVKhoKBgLOF0d3cLHNnMsSyLhIQEXges2tra\noFarsXz5ct7OKbTI6nChKJEIBALQ6XTgOA45OTk3rOxJSEiIuEGPQCDA+3rwPXv2ICkpKeKnGY1H\nW5oUNUM+nw9qtRoOhwMqleqGhOnxeKDT6US3RHI6brcbNpuNt/M5HA5YLBY89thjvJ1TDGhLk6Jm\nwOFwwGQyITU1FUlJSRMe4/F4kJqaGubI5k4qlUKlUvF2vmPHjmHFihVITEzk7ZxiQFuaFBUEjuPQ\n3d0NnU6H3NzcSRMmMFKBPdKmGgEja+P5Wnfu9/vxu9/9Dhs2bODlfGISeb9ZigoxlmVRV1eHpqYm\n1NTUYMuWLTCZTIiPj0d2dva0CdHtdkMul4cpWv7IZDIkJyfzcq5Tp05h0aJFWLx4MS/nExOaNClq\nHJZlsXXrVjQ2NsLpdEKhUGDZsmX46KOPgt6Gwm63i6pGZrC8Xi9va+MPHTqE+++/n5dziQ29Paeo\ncerq6tDY2AiHwwFCCJxOJ86fP48TJ04EfQ69Xi/Ksm/TcblcYFl2zudpbW2FRqPB6tWreYhKfCLv\nN0tRIdTU1ASn03nNYy6XC83NzUGfIyYmBllZWXyHFnIJCQlT9tUG64UXXkBqampUTTMajyZNirqK\nEAKdTnfDmz0hISHoydmEELS3tyM+Pj4UIYbU4ODgnLe6sNlscDqdeOKJG3bcjho0aVIURgZvXnrp\nJUgkEqxYsQJKpRIMw0CpVKK2thbbt28P6jxOpxOlpaURORCUnp4+5wn5Bw8exJo1a3gbUBIjOhBE\nzXu9vb146aWXUFpait/85jeQSqWoq6tDc3Mzli9fju3btwd9qzk4OIi+vr7QBhwiVqt1Tn2xXq8X\nv//97/Hyyy/zGJX40KRJzVuEENTV1WHfvn347ne/i/Xr1499b8eOHdixY8eMzymTySK2BJrP55vT\n6Pmnn36KjIyMqKlmNBmaNKl5yePx4JVXXsHAwAB+8Ytf8PZGP3/+fEQOgLAsi/T09Dm1NOvq6vDj\nH/+Yx6jEiSZNat45c+YM9u3bh4qKCvzLv/wLr2vEk5KSkJuby9v5wsXj8UCtVs9607P29nYYDAas\nXbuW58jEhw4EUfMGIQQnTpzA//7v/2LFihX4p3/6J96LarS0tETkEkqGYbBo0aJZ//zPfvaza0rj\nRbPI++1S1Cx4PB7s2bMHfX19ePrpp5GXlxeS60ilUlFv0TsZp9M560rzw8PD8Pv9ePTRR3mOSpxo\nS5OKet3d3XjsscfAcRx+/vOfhyxhejweGAyGiKxwRAiZ9TYX7733HrZs2RKRHxazQVuaVFRrbGzE\na6+9ho0bN+Lee+8N6Ta6w8PDvFc+DxdCyKx2oXS5XPjjH/+IPXv2hCAqcaJJk4pKLpcLb7zxBrq7\nu/H4448jPz8/5NfU6XQR29oaHh6e1cj5Bx98gNjYWBQVFYUgKnGiSZOKOnq9Hs899xzS0tLw1FNP\nhW1Jo0aj4X3723BJSkqa8e05IQQff/wx/v3f/z1EUYkTTZpUVGlqasLevXuxY8cObN68OazVhrxe\nL8rKysJ2PT7p9Xrk5OTM6GcuXrwIq9V6zaKA+YAmTSoqeDwe/OlPf4JOp8MPfvADLFy4MOwxXLx4\nEbW1tWG/Lh9iY2Nn3CJ/+OGHsXXr1ogsgzcXNGlSEW94eBi7d++GyWTCk08+KVixiOTkZCxZskSQ\na89VIBCY0RJKnU4Hj8eDBx98MIRRidP8+oigok57ezt+85vfYOnSpXj22WcFS5gcx+Hzzz+PyJJw\nAGA2m2c06v/aa6/hrrvuisi6oXNFW5pUROI4Di+99BK6urrw8MMPo7S0VNB4nE4nKioqIrYkWmZm\nZtADQVarFX/961/x9ttvhzgqcaItTSriOBwOvPTSS3A4HHj88ccFT5jAyBYPRqNR6DBmxePxQKPR\nBN3SfPnll0EIQWFhYYgjEyfa0qQiSnt7O1577TWsWLECO3bsENVa58rKSqFDmLXs7OygjiOE4ODB\ng/j9738f4ojEiyZNKiIQQlBfX4+9e/fi+9//PlauXCl0SNeI1JJwwMhUqWBbmQ0NDeA4DjfffHOI\noxIvmjQp0fN4PNi3bx+0Wi1++ctfBt0qCie/3z+nKkFC8vl8CAQCQR37ox/9CN/4xjcicqkoX2if\nJiVqra2t+PWvfw2ZTIZHH31UlAkTAAYGBkJWCCTUgq3MdOXKFRiNRnzzm98MQ1TiRZMmJVrNzc34\n85//jPLycuzatUvUm5Xp9XqkpaUJHcas6PV6WK3WaY97+eWX8d3vfhfp6elhiEq86O05JTosy+LD\nDz9ES0sLfvSjH2HBggVChzQll8sFQkhYioKEglKpnLZwssFgwOuvv44jR46EKSrxoi1NSlSGh4fx\n5JNPQqPR4LHHHhN9wgRGdrP0er0ROxBkt9unPeY//uM/IJVKI7YLgk+0pUmJRk9PD/bs2YNly5bh\n3nvvjZhtI5xOJ2pqaoQOY9b8fv+UCZ9lWZw4cQK7d+8OY1TiFRl/lVRU4zgO9fX1OH78OL7+9a9H\nXKWg48ePQ6PRCB3GrCkUiikLEB88eBCpqakRW4yEb4InTYZh/h+AHQB0hJClVx9LBfA6gEIAfQDu\nIYRM31NNRRyv14sXX3wRBoMBP/nJTyJyMCUlJUUUq5Jmq6enZ8odNH/2s5/h8ccfn9fTjMYTQ5/m\nnwFsve6xJwAcIoSUATgC4MmwR0WFnMFgwP/8z/9ApVLhF7/4RUQmTGCkpen3+4UOY9ZUKtWk686b\nmpqgVqtx9913hzkq8RK8pUkIOcEwzPWLWHcC2Hj161cA1GMkkVJR4siRI/jss89w2223YfXq1RHd\nipHJZKiurhY6jFnx+/1TbnXx0ksv4YknnpjV/kHRSvCkOYksQogOAAghWoZh5l/9qSjFcRwOHTqE\njz/+GA8++GBEr9cGRpKOWq2OiFH+ibAsO+ne73q9Hu+99x7a2trCHJW4iTVpXo8IHQA1d06nE/v2\n7UMgEMAvf/lLKJVKXs8fCATAsiy8Xi/sdjsUCgU4joPT6QTDMIiPjwfDMLBarUhLSxsbnR/93kwr\nkLMsi71796Kvrw+HDx/G9u3bI27aUSAQQEpKyoTf+9a3voXS0tKI7TYJFbEmTR3DMNmEEB3DMDkA\n9JMd+Mwzz4x9vWnTJmzatCn00VEzduHCBbzxxhtYu3Yttm3bNqctEgghcLlcMJvNkMlk8Hg8GBoa\ngkKhQHJyMhiGgdPphEwmQ0xMDCQSCWJiYq65ps/ng8fjgcPhgNfrhUKhgEQigcPhQGZmJhQKBeRy\n+aTVzFmWxdatW/H555/D7XZj165dqK2txYEDByIqcdrtdpjN5mseY1kW77//Pk6cOIGnn34aLMtG\n1HOajfr6etTX1wd1LEOI8I04hmEWAviAELLk6v+fB2AmhDzPMMzPAKQSQm7o02QYhoghfmpqzc3N\neP/997F8+XLceeedszqHz+fD4OAgWJYFy7IghCAQCEClUiE2NhYymWxOyyxHz2ez2cAwDPx+PywW\nC3w+H/Ly8hAXF3dNEv3www+xa9cuOByOsXMolUq8+uqr2LFjx6zjCLfh4WG43e6xSeujHwYnTpyA\n1+uFUqmMyA+DuWIYBoSQCTvaBW9pMgyzF8AmAOkMw/QDeBrArwG8wTDM/QCuALhHuAip2WJZFvv3\n70dXVxe+973vzbjYhsvlwsDAACQSCQKBABiGQWpq6jW31nxhGAYymeyaddXZ2dnwer3wer1wOp3o\n7OxEXl4elEolzp07B6fTec05nE4nmpubIy5pjm941NXVobGxEV6vF8BIwefGxkbU1dVF1PMKJcGT\nJiHkvkm+9cWwBkLxyuFw4LnnnkNBQQG+//3vB713js/nw/DwMBwOBzweDziOw6JFi5CQkCDICHts\nbCxiY2ORlJSE9PR0eL1euFwuZGZmIj4+Hi6Xa+zYhIQELF++POwxzgXDMNcMBDU1NUXFh0EoCZ40\nqeijVqvx+uuvY/Xq1fjSl74UVKswEAjAaDSir68PcrkcxcXFSExMFNX2sBKJBAqFAgqFAvfffz9e\nf/11nDp1Ch6PBwqFAqtXr8b27duFDnPGFArF2Nc1NTVR8WEQSqLo05wt2qcpPo2NjTh27BjuuOMO\nVFRUTHu83+9Ha2srAoEAcnJykJmZOekUGLFhWRa/+c1v4HA4UFJSgttuuw0pKSlQKBQRM+/01KlT\nWLhw4diuki6XC5s3b0ZbWxucTicSEhJon+b134vkpEOTpniwLIs33ngDZ86cwU9/+tNpt3b1+Xww\nGo0wmUxQKBRTrkoRq0AggJ6enrEllD6fDzabDXq9HoWFhRExIfzy5cvXvPZnz55FQkICuru70dzc\njOXLl0fkVKq5EvVAEBX5HA4H3nzzTcTGxuKZZ56Zdv6lyWTCqVOnsHjxYpSVlYm6uPBURvtdR8nl\ncmRkZCAhIQFOpxNutxsJCQmi3gt9cHAQKpUKAMbmuFZUVKC8vJz2YU6CJk1qTpqbm3H48GGsXLkS\nGzZsmLIP0mazQavVguM4rF+/HomJiWGMlH8cx034HOLj4xEfHw+n04n+/n5kZ2cjKSlJVP2zwMg0\nq4SEhLEPLaPRiNzc3Gv6OKkbieu3SEWUCxcuYN++fVi1ahU2bdo0ZVLQ6XQ4ceIEFAoFysrKIj5h\nAiP9f1P1vyYkJKCkpAQMw8BkMsHtdocxuun5fD5IpVLExMSA4zhotVrR7sEkJrSlSc0Yx3E4fPgw\nurq68Oijj07Zf+l2uzE4OIhAIIDNmzeL+lZ1plwuF5KTk6c8RiKRIDk5GR6PB/39/cjJyUFSUlKY\nIpyaz+cbGyVXq9XQ6/VYtmyZwFGJH21pUjNiNpvxX//1XzAYDLj//vunTJgOhwNHjx4FwzAoKyuL\nqoQJjCTEYFvMcXFxKC4uBsuyGB4eDnrL3FBiGGZsJZDT6cSqVasEjigy0KRJBU2v12Pv3r1IS0vD\nPffcM2USNBgMuHTpEtasWYOSkhLR9efxYWhoaEYrkyQSCVJTUxETE4P+/n7BE6fFYoHNZoPD4Ziy\ncAd1LXp7TgWlo6MD+/fvx6233oqqqqpJjyOE4Ny5c/B4PFi5cmXETSMKFsuyyMvLm9XIf3JyMuRy\nOaxWK5RKpWDzUqVSKVJSUtDR0QGVShUxc0uFRpMmNSVCCP7v//4P/f39uP/++6fcFoHjOPT29kIq\nlWLVqlURM0l9Nmw2GywWy6x/Pj4+HlKpFCaTCXFxcYK08rxeLwghGB4epit+ZiD67pko3vh8Przz\nzjvw+Xx44IEHpkyYgUAA586dQyAQQHV1dVQnTACIiYmZdgL/dGQyGdLS0sbK1IWbw+GAwWBAZWVl\nxOz8KQb0laImZDQa8d5770GlUmHnzp1TrgjhOA7t7e3weDxYsWJFVPZfXm94eHjSWpszIZfLkZ6e\nDpvNBo7jwjpHMj4+Hlardc7Jf76J/r9uasZ6enrw3HPPobKyElu3bp02YV6+fBlyuRzr1q2bFwkT\nGHnefPXXSiQSJCUlwWg0XlMoI9SOHz8On8/HS/KfT2hLkxpDCEFzczNOnjyJ733ve0FtS3vy5EnI\n5XKsXLly3iRMYKTQCJ+tQolEgry8PDgcDkil0pAvLSWEIC4uLuL3aBICTZoUgJH+y7q6OpjNZtx3\n331BDUyo1WrExsZi6dKl8yphAiPr5xctWsTrOaVSKZRKJSwWC5RKZUhnHphMJkgkEmRmZobsGtGK\nJk0KLpcLf/vb3+BwOPDDH/4wqEEcg8EArVaLJUuWRGzBjdny+/3IyckJyfOWSqVQKBSw2+2Qy+Uh\n+zC6fPmyYIWdIx1NmvOc0WjEu+++i6qqqqD7JJ1OJ44dO4YvfvGLUT9KPhGLxTKn6UbTUSgUiImJ\ngdPphFKp5D2xmUwmqNVqumRylmg9zXns+PHjOH/+PLZs2RJUwWBgZADkwoULyMrKGluCN98MDw/D\n4/FMOQWLD1arFYQQ3udwdnd3j23bcdNNN/F67mgxVT3N+dURRQEYGQQ4ffo0jh8/ji984QtBJ0wA\n6OzsRFxc3LxNmABu2PI2VBISEsCyLK/LLT0eD1wuF7KysuhUo1miSXOeCQQCOHDgALq7u/HII4/M\naPTUYDCgra0NRUVFIYxQ/DiOC0vxEalUiuTkZLhcLnAcx8s529rakJSUBLfbTfszZ4kmzXnEaDTi\nv//7v8FxHO6+++5pK6yPx7Is+vr6sHnz5nk38HM9vqcbTWW03qXVap3zufx+P6xWK3JycgRZgRQt\naNKcJ/R6Pd577z0UFhZi27ZtM57Q3NXVhaSkJKSlpYUowsgQCAQwPDwc1g+O0Y3a5nqbbjAYUFxc\njLi4OMHWu0cDOno+D3R3d+Po0aPYvHkzSkpKZvzzNpsNFy5cwM6dO0MQXWgEAgEEAgE4nU5IpVJ4\nvV7ExsaC4zjIZDLIZDLI5fIZ36L6/X7k5+eHKOqJxcTEQKlUwuv1znqNOMdxMBgMKC8vBwD09/fP\n6E6D+juaNKMYIQTvvfceuru78a1vfWvWHf9XrlzB+vXrRV3mzePxwGw2w+12w+12w+fzITExET6f\nDykpKXC73VAoFGOVfQghcDqdSE1NBcdxyMjIQGxs7LRJyWw2w+FwhOlZ/Z1UKoXb7R7bVnemOjs7\n4XQ6x6aIpaen06Q5SzRpRqlAIICjR4/C6XTOKWEaDIawTK+ZDY/HA7vdjt7eXiQmJoJhGGRnZ4/d\nfk7XiuQ4Dl6vF3a7HVarFXa7HUlJSUhISBg73/WkUqlgq2hGN2sjhMyohUwIQSAQuKYOqsVioWvO\nZ4kmzShkt9tx8OBBKJVK3HPPPXN6c3z22WdYtmyZqJZJut1uDAwMwGazITk5GaWlpbPa7TEmJmZs\n50hgJIlaLBa4XC7odLqxnRnHn9disQi2+ZhUKkVsbCz8fv+M+lQtFgvkcvnYfkY+nw8cx9HR81mi\nSTPKGI1G/OEPf8CmTZuwcePGOb0xTCYTVCoVFi5cyF+Ac8BxHDo6OsCyLORyOZYvX85rHciYmJix\nga7k5GTYbDaYTCZkZmaOjZZ7PB5BZw/I5XLY7faxUfVg9PT0oLi4eOz/LMsiJycnVCFGPZo0o8iV\nK1dw9OhRfO1rX5tyS4pgDQwMoKCgQBQtEoPBgJ6eHqSmpqKoqCjk/aujLVCv1wuNRoPExEQkJCRA\noVAIui94TEwMJBIJ/H5/UEtYh4eHodfrsXLlyrHHLBbL2Eg6NXM0aUYBjuNw+vRpdHZ2Yvv27bzc\nPg4MDODKlStYunQpDxHOHiEEra2tYFkWixYtCvsqltjYWBQWFsJqtaK7uxsWiwWLFy8OawzXUygU\nQc+zNJlMNyyVjIuLC/sMgGhCk2aE4zgOBw8eRFNTE37wgx/wtqe2zWZDZWWloH2ZXq8XbW1tkMvl\nKC8vF6w4CMMwSElJgd/vH9srXOjWZkxMDHw+35RdBW63Gy6X64YVXFarVVR91JGGvnIRzOPxoK6u\nDgDw4x//mLeE6ff7YTabBe3LdLlc6OrqgkwmQ0VFhSiqKdlsNuTk5MBut4e0ylEwGIaB1+ud8piL\nFy8iOTn5hgTp8/mmrMZPTY0mzQg1ODiIN998E2lpabjtttt4TSodHR1QKBSCTUnx+Xxoa2uDRCJB\ndXW1aFpFgUAACQkJyMzMHNvFUSixsbGQyWSTrkn3+Xxwu90TThWLjY3l7QN2PhLHXyM1I1qtFu+8\n8w5yc3Oxdu1a3pOKXq8XbC4ix3FobGxESkqKqLZiGJ0bqVQqERMTMzZ9R4iJ7qNiYmImTZo6nQ5F\nRUUT3r739PSAZdlQhxe1aNKMMJcuXcKRI0dw1113YcuWLbyff3RKjVADBWfPnkVmZuaslnuGksvl\ngtFoHLutjYmJQVJSElwuF6+l22ZCIpFMuBGb3++HyWSadEAwJydnVquKqBF0IChCcByH999/Hx6P\nB7fffnvIii2o1WpkZGQIcks8ODgIhmFElzABjI3ejyeRSJCWlgaTyYSMjIyw9xNKJJKxW/Txv6/2\n9nb4/f4Ju1cCgQBdDTRHtKUZAQKBAOrr62EymbBt27aQVqdpaWkRZPK22+3GxYsXUVpayuuEdb4Y\njcYJP0ggdwPsAAAgAElEQVSkUimSkpLgdDoFiGrk+uNv0TmOA8dxk3Zt0EGguaNJU+Tcbjf2798P\nqVSKf/zHfwxpwgwEAkhOTkZhYWHIrjGZM2fOoKqqSrQDFF6vd9IJ9fHx8eA4Dj6fL8xRjXQTjO8e\nMBqNYxPxJ0IIoRXb54gmTRFTq9X405/+BJVKhVtuuSXkLTCr1YqkpKSwt/T0ej1iY2OhUqnCet1g\njZaTS0xMnPQYpVIJu93OW4X1YMXExIwN6nAch97e3imXSBqNRphMpnCFF5Vo0hSpoaEhHD9+HLW1\ntaipqQnLUsaurq6wJ0xCCC5cuIC8vDzRTC26nsVigdlsnjK+0a13/X5/GCMbma8pl8vBcRw0Gs20\npeMUCoVoP5wihfg6jyh0dnbiwoUL2Lx5c1gr6vj9fqSnp4ftesDI4E9iYqKol/URQlBQUDDtcTKZ\nDE6nc1bFjefi6s6JsFgsWLVq1ZTHms1m0XaBRAqaNEWE4zh89NFH0Gq1uPfee8P6x00IAcuyYe/v\n0uv1gpVaC5bFYglq3qpUKoVMJoPP5wvrCiZCCEwm09g0qKn4/X5RDrRFEnHeD81DgUAAJ06cgEQi\nCXvCBEYmaY/e6oWL1WqFRqMJqhUnFI7joNfrg66qJJPJQAgJcVTXkkqlaGlpQUZGxrTHxsXFTdk3\nS02PJk0RsNvt+OSTTyCVSrFt2zZBbp80Gg0vOx7OxPDwMIqLi0XblwmMzF5QqVRBf5jIZDL4/f6w\nDgiN1tcMpmvl0qVLon69IwF99QRmsViwe/dupKSkhGRJZLDkcjlKS0vDes3u7m7Rr0zR6/Uzngg+\nuoFbuOh0OpSXl0/bj8qyLBYuXDjvt2CeK9q5ISCNRoOTJ0/ijjvuEHwVzOjGY+EyWqFH7CO5Tqdz\nxoNjo9OAwtF36HQ64fF4UFBQMO3eQS6XC263O+QxRTva0hRIS0sLGhoasH79esETJjBSdDicbDYb\npFKpqFeneL1eJCQkzLi7JJzP6cyZM0FXpJrp3kLUxGjSDDNCCM6ePYsjR45g8+bNolmdwXFcWLfo\n1ev1op/6Mjg4OKvuktHWXqgHhNxuN5KSkoJewcVxXNinlEUjUSdNhmG2MQzTwTBMF8MwPxM6nrkK\nBAJoaGiA1WrFgw8+GNIlkTOlUCjCug82IWSsvJoYjW57O7rR2kzFxMSEPGlqtVrk5eVBKpWO7eU+\n3fE2my2kMc0Hok2aDMPEAPgDgK0AqgDsYhimXNioZs9qteKjjz4CAGzcuHFs21ixCAQCYZ2/p1ar\nw3o9u92Om266Kej6ly6XC/Hx8bOansOyLD7++GP88pe/xIcffhiS2pV2u32sulKw10tOThb9nNhI\nIOaBoNUALhFCrgAAwzCvAdgJoEPQqGbBbrfj0KFDiI+Px5o1a0Sxu+P1JislFirp6elhvT3/+OOP\ncebMGXz88ce45557pj3+8uXLs+o6YVkWW7duRWNj49iSxtraWhw4cIDXvs7Lly+P9U8Gez2DwSCq\nu5tINW1Lk2GY/49hmNRwBHMdFYDxoxODVx+LKEajEceOHcPKlStx++23izJhAsCtt96K/fv3h+16\nMpkMZ86cCfl1HnjgASQlJeEf/uEfAADf/OY3kZSUhAcffHDSnxndsCyYyeLXq6urQ2NjIxwOBwgh\ncDgcaGxsHNvLiQ9erxcMw6CqqmpG16OrgfgRzCuYDeA0wzDnALwM4AAJ95KHKTzzzDNjX2/atAmb\nNm0SLJbrnTlzBn19fbjllltEf1sU7vmhBQUFYam289RTT6GxsRHd3d1jSWPhwoV48sknJ/0Zk8mE\nrKysWSWYpqamG2prOp1ONDc3Y8eOHTM+30S0Wi0yMzMhkUiCvh7HcUhMTBRdt5BY1NfXo76+PriD\nRzuQp/oHgMFI3+JrALoBPAugOJifne0/AGsA7B/3/ycA/Oy6Y4hYdXR0kNdff5309/cLHcq0OI4j\nfX195MMPPwzbNU+dOkUOHDgQlmu9/fbbRCqVEqVSSaRSKXn77bcnPdbpdJITJ06QQCAwq2t98MEH\nRKlUEgBj/5RKJfnggw9mG/417HY7OXz4MPH7/TO6ntVqJYcPH+Ylhvngam6ZMDcF1by4ehLt1X8B\nAKkA3mQY5jfBpeZZOQ2ghGGYQoZh5AC+DuD9EF6PFxzH4dy5c1Cr1bjjjjuwYMECoUOaFsMw6Ojo\nCOuaaYfDgSVLloTlWm+++SYSExPxq1/9CkqlEm+99dakx5rNZpSWls66/3H79u2ora1FbGzs2EZs\ntbW12L59+2zDv4ZWq0V2dvZYK3j0evHx8VNej2GYG/Y/p2aHme6NwjDMjwB8C4ARwG4A7xJC/FdH\nty8RQopDFhzDbAPwAkb6Xv8fIeTX132fhPONPh2fz4c333wTKpUKN998c0T1Hx0+fBjr1q0L2+3b\nZ599hgULFoSlWEdLSwtyc3ORkZEBo9EIrVaL6urqG46z2Wwwm80oLCycU98zy7J45JFHEAgEcMcd\nd2D79u28DAL5/X60traiurr6mr8tlmXxr//6r2hubsZDDz004fUGBgbg9XpFsZAiElwttzfhH0Ew\n7+o0AHeRq6PYowghHMMw/HTSTIIQsh9AWSivwRe3243Gxkbk5+dj3bp1EZUwgb8XmghX0pRKpdDr\n9WFJmuNbtBkZGRMO8LAsi4sXL6KiomLOg3USiQTl5eX40pe+xGvrTqvVIiMj44a/LYlEgptuugll\nZWWT9psGAgG6Gogn096eE0Kevj5hjvteO/8hRR6r1Yrjx48jPz8fGzZsiLiECYysgw/nNghZWVmi\nmklgMBhQVFSE1NS5TxTx+XywWq3Iy8vjIbIRbrcbp0+fnnQaFCFkyq4giUQi6sUEkUS0k9sjRX9/\nP/bu3YvKysqIvvUpLCwM6wh6QkICLl26FLbrTcVoNMLlcvG2pHVwcBAZGRm8Lks1Go1Yvnz5pK3F\n4eHhKee99vb2hn3/omhFk+Yc9Pf3o7OzE3feeaeot2sIRkpKyljloXDIyMhAamoq7HZ72K45EZ/P\nh8uXLyMtLY23D42WlhZe7zb8fj/MZvOULcnu7u4JX0uWZfHhhx/irbfewqeffhqS1UnzTeTdR4pE\nW1sbtFot1q5dG9Y126EikUig0+nCVlMzJiYGDMNArVajvFyY1bEsy2JoaAjl5eW83rpqtVps2LCB\nt/P19fUhPT19yhVbEokEK1asuOax0dVJDQ0NcLlc+POf/xyS1UnzDW1pzhDHcThy5AhaWlpw8803\nR0XCBP6eNMOpvLwcV65M2F0ecoQQtLa2Qi6X85owbTYbLBYLb3vH+3w+tLW1TVmdyOFwjE1eH290\ntZDT6QzZ6qT5iCbNGfD5fDh16hRSUlLwla98JaybZ4Vafn4+4uLiwtrvlZOTA4/HA5fLFbZrjhrd\na32qPcJn49ChQ1CpVLwVdNZqtVi5cuWUsxra29vH5mmON9VqIWr2aNIMktlsxmeffYbk5GTU1NRE\n3fQNuVyO2NjYsCYwuVyO1NTUsKxBH0+r1cLj8aCkpITXwa/R2/3Vq1fzcj6fz4fh4eFpE7vBYJhw\nEKumpuaG7UQSEhKwfPlyXuKbr2jSDILVasXBgweRnJzMyzw+sXK5XGEf0V62bBlsNlvYtmEwmUzo\n6elBXl4e7/16p0+fhs/n461f+OLFi0hKSpp2UMlgMEy43/noaiGZTBaS1UnzFU2a0zAYDDh79iw2\nb958Q0d7tCkqKgr76GpycjKkUik6OkJf8c9ms8Fms2HVqlW8l8ELBAI4deoUdu7cycv5vF4vbDYb\ncnNzpzyOEIK+vr4J+1AlEgkOHDiAJUuW4Fvf+hZeffVVOgjEA5o0p3D+/Hm0tLRgxYoVoq9SxIcF\nCxZgeHg47NfdsGED2traQnptg8GAjo6OGW3HOxOnT59GXFwcFi1axMv5NBoNysrKpp3rOTAwAKvV\nOumkfIlEgvj4eDz11FPYsWMHTZg8oElzEt3d3ejq6sLy5cvnTeHWpKSksdHfcFIoFKisrMSFCxdC\nMhBlt9ths9lQWVkZkoRpNptRX1+PO+64g5c+Uo/HA6vVGtRk+/7+fqxevXrS69rtdhBCInrhhdjQ\npHkdjuPQ0tICs9mMnTt3znqPmEjEMAxkMhn6+/vDfu2lS5diaGgIly9f5vW8Q0ND6OvrQ35+fkim\nh7Esi7feegvr16+f9lY6WKdOnYJCoQiqVdje3o7KyspJv3/mzBn4/f6w10uNZvSVHMfj8WD//v1g\nWRYrV66MuhHyYKxatUqQzbckEgl27NiBkydPQqPR8HLO0T7MwsLCkE0P+/TTT8FxHNauXcvL+Vwu\nF5RKZVC3+TabDZ2dnVi8ePGkxxBCRFWYOxrQpHmVx+NBU1MTkpKSsHTp0nnb95OdnQ2z2Ry20ezx\nEhMTsXbtWjQ2Ns65i2BgYAA6nQ5FRUUh24vo/PnzOH36NO677z7elk0ODQ1hwYIFQZ2vvb0deXl5\nU87h/PTTT+ft33Ko0KSJkX6fs2fPQqVSYf369fP6VkYikcDv9ws2AXrx4sVIS0vDZ599Nut16Waz\nGRaLBXl5eSG7Wzh37hwOHz6Mb3/727PasXIiJpMJDocj6L2JLly4gC984QtTHmO1WsNSfm8+mb/Z\n4ar+/n4cOHAAixcvpn9cV9XW1gqySmfUhg0boFQqsX///hnH0dvbC6vVioqKihsmdvOloaEBDQ0N\n2LZtG2+VkQCgo6MDCQkJQc0D1mg06O/vn7I/ExhZAVRbW8tXiBTmedLUarXo7+/HmjVreP3jj3Qq\nlQp+v1+Q6UejbrnlFpSVleG1114L6ladEAKdTgebzQaVShWSmqaEEBw5cgTnzp3Dxo0bp01YMzE6\nbSjYUe7jx48jNTV12pb0uXPneFvSSY2Yt0mzt7cXV65cQU1NTcSXdeMbwzDIzMzE559/LlgMMTEx\nWLp0KTZu3Ih3330Xra2tk+5hxHHc2G6T1dXVIbkld7lc+OCDD9Da2oq77roLVVVVvJ2bEAK1Wo28\nvLygWpmBQACXLl3CfffdN+VxDocDmZmZKCuLiM0PIsa8Kw3HcRw+//xzOJ1ObN68eV6OkAejpKQE\nnZ2d8Hg8vBbTnani4mLExsbi3Llz6OzsxO23335NPCzLQq1Ww+v1oqioKCSDHi0tLWhsbIRcLsf3\nv/993lcT9fX1wel0Bj0fuKmpCdnZ2dOuST9z5szYEkqKP/MqaQYCAXR0dCAxMRE33XQTTZhTSE5O\nRmVlJbq6urB06VJBY8nPz0dmZibOnz+PPXv2oKamBkuWLAEhBN3d3UhPT0dVVRXvyUGj0aChoQH9\n/f3YsmXLhJuxzRUhBHq9PuiVRBzH4Z133glq/bjZbKb99CEwb5Kmy+VCR0cHlEolFi9eTD99g1Ba\nWop3330XxcXFIRtUCVZsbCxWr14NlUqF9vZ27N69GwsWLEBhYSHvS1x1Oh0++eQTuFwu5OTk4Ac/\n+EHI9n0yGAxIT08Puk+9tbUVOTk5WLdu3bTH6vV61NTUzDVE6jrzImm6XC588sknKCkpCVtl8mig\nUCigVCpx4cIF3iZvz5VKpUJCQgLi4uKgVqvR1NSE7u5uLFmyZKwm6GwYDAZcunQJbW1t4DgOKpUK\nX/nKV0L6YREIBKDRaIL+m+Q4DocOHcLtt98eVDfEsWPH6Mh5CER90hxdNVFbW8t7wdn54NZbb8XB\ngwfhdDoFb20CI1NoNBoNiouLsX79eng8Hpw6dQpdXV3Yv38/VCoVlEol0tPTkZGRgfj4eCiVSkil\nUhBC4PV64XA44HK5cOXKlbEuG6VSifj4eGzcuBGLFi0Ky46ira2t8Pl8QW+b/P7778PlcgWdZNPS\n0nD77bfPJURqAlGdNPv7+6HValFaWjpvim7wLT4+HsXFxTh58iS++MUvChqLWq3G8PAwiouLxxJN\nXFzc2H48Pp8PV65cgclkGpsdMTg4iOTkZBiNRiQlJcHr9UKhUEChUEAul2PhwoUoKSnBwoULw7qo\nwe/3gxASdD+px+NBW1sbdu7cGfQIe0tLC2/r4am/i9qkOTQ0hObmZtxyyy00Yc5RSUkJPv30U1RW\nVvK6l/dMmM1m2Gw2FBQUTNoyk8vlWLx48ZRrscViaGgImZmZQbcyDx8+jPLy8qCTrNFoRGpqKm+r\nlai/i7p5mqNFWQ0GA7Zu3TppnUEqeHFxcfjyl7+MpqYmBAKBsF6bEILOzk4YDAYsXrw4ZOvIw8li\nsUCr1QbdXdTS0oITJ05g69atQV+jt7eXzs8MkahKmoFAAGfOnIHT6URVVVVUbXwmtLy8PKSmpqKh\noSFs1+Q4DhqNBlKpNGz9jOHQ09ODpKSkoAZzAoEAjh49ii1btsyoT7mhoQFWq3UuYVKTiJqkOdqh\n7/P5UF5eHjVvMDGpqqpCc3MztFptyK/FsizOnTuHQCCAoqKiqJlTa7PZIJfLg97r/fDhw8jIyMCW\nLVtmdJ2EhATceeedswmRmkZUJE2v14v29nYkJSVh3bp1tBRWiCQnJ+Ouu+5CY2MjvF5vyK7j8XjQ\n398PlUqFgoKCqJlTSwjB4OAgVCpVUM+psbERx44dC3rwZ7xPP/1U0JVc0Szik6bJZMLnn3+OrKys\nqHqDiVVeXh4KCgqwf//+kGxN4XK5cObMGSQmJkbdyG93d/eU+/mMZ7VaUV9fP+u5oizL0i0uQiTi\nk2Z/fz8WLVo0LzY+E4uqqiqYzWbea26aTCZoNBosW7Ys6JqSkYLjONhstqDmWHq9Xrz66qtYu3Yt\nbrrpphlfy2q1ghCCBQsWzCZUahoRnzRH59hR4SOXy7Fr1y4MDg6ir6+Pl3OaTCZ0dXUhJycnKqfJ\naLVapKWlIT09fcrjCCHYvXs33G43brnlllldq62tDR6PZ14X0w6liH9Vo/ENFglGJ5V/8skn6Ozs\nnPV5OI5Db28vHA4HVq1aJYpVR3xzuVzQarVQqVTTHnvkyBFIpVI88MADs+5q8vl8uO2222b1s9T0\nIj5pUsJJSUnBmjVr0NTUBLPZPOOf5zgOarUadrsdKpWK95JrYtHR0QGGYaadAfDRRx/h4sWL+OY3\nvzmnnTOPHDkS9vm08wlNmtScLFmyBDfffDOOHj0KvV4f9M95vV60tLRALpdjyZIlUTtFzOl0QiqV\nTlte7/Tp02hoaMCXv/zlObe2JRIJr0WSqWvRpEnN2YIFC1BWVoY9e/bAZDJNe7zX68Xg4CBSUlKQ\nnZ0d1TMeBgcHkZeXN+U0uNOnT+PEiRN47LHHUFhYOOdrajQauhoohKLz450Ku+rqaiQkJODIkSPY\nsGHDpLMZLBYLOjo6UFlZGRVLIqfS19cHi8Uy6Yg5IQRvv/02Ll68iEceeYSXJb8ejweXL1+me16F\nEE2aFG8WLVqEuLg4fPDBB6ipqcHKlSuv+f7o1rrl5eVRnzA5joNer0dhYeGELWmO41BXVwedTof7\n77+ftxoJer0ea9eujZoVVGJEkybFq9zcXGzcuBGnT58e68tjGAaDg4NQq9WoqamZF29onU6HjIyM\nCYty2Gw2vPLKK0hNTcV3vvMdXmskHD9+HA6Hg7fzUTeiSZPi3eLFi6FSqXDs2DH09PSM7UG+YsWK\nqB0hH8/tdkOn0024xa9Wq8Ubb7wBqVSKe++9l/fXQyqVYsWKFbyek7oWHQiiQkKhUGDLli0YHBzE\n/v37IZPJ5kXCBIALFy4gJibmmhY1x3E4cOAA9uzZg82bN+Ohhx4KyevR3d0d9CZt1OzQliYVEg6H\nA319fdi1axecTieOHTuG3NxcrFu3LmqnFwEjz1uhUFzTyrTb7fjwww/R1dWF7373u0FNcp+trq4u\nWt0oxKL3r5cSjMvlwuDgIFJTU5GRkYGMjAykp6djz5490Ov12LBhQ1SO7hJCMDAwAJVKBYlEMrYR\nWmtrK1auXIl77rknpBW4/H4/0tPTJ+wWoPhDkybFK41Gg6GhIVRUVEChUIw9rlQq8eCDD6K3txdH\njx6FQqHAhg0bkJycLGC0/Orp6YHH40FaWhqGhoZw8OBBqNVqfOMb3whLfYSWlhbo9XpaGjHEGEKI\n0DHMGsMwJJLjjzZ6vR52ux3Z2dlTLgP0eDx455134PP5sGzZMlRXV0f8LTvLsmhoaEB6ejrq6+vh\n8/mwceNGLFmyJGyFMxobG3Hu3Dk89NBDYbleNGMYBoSQCVdd0KRJzRnHceju7obL5cKSJUuCbunY\nbDacPXsWFy9exPLly1FbWxux05FaWlpw9uxZDA8PIz09HXffffc1Le1wePnll5GUlIS77747rNeN\nRjRpUiHj8/kwMDCA+Ph4ZGVlzarFeOnSJfT29o5NBi8vL0dmZmYIouXXaCX2gwcPoqGhAbt27cLa\ntWuD3mGSb//2b/+G2tpabNu2TZDrR5OpkqZg90QMw9wN4BkAFQBuIoScG/e9JwHcDyAA4EeEkIOC\nBElNyefzoaurC/Hx8XPa2nd021273Y76+nrs2bMHZWVlKCsrQ0FBgehu3UdnA/T09CA2NhYsy+Lh\nhx9GTU2NoHHJ5fIbVmFR/BOspckwTBkADsAfAfzzaNJkGKYCwF4ANwHIB3AIwOKJmpS0pSkco9GI\noaEhLFy4kPclkYFAAL29vWhtbUVLSwtuvvlmLFq0CPn5+YLN9bRarWOrbfR6PTIzM7F48WIUFRVB\no9GgoqJC0KK/DocD9913H9555x06EMQDUbY0CSGdAMDcuDB3J4DXCCEBAH0Mw1wCsBpAY5hDpCZh\ntVphMBiQl5cXkjXkUql0rPW5YcMGDA0Noa2tDfv27UN+fj7KysqQk5ODjIyMkG0e5nA40NXVhb6+\nPlitVjidTgDAxo0bcffdd0MqlYLjOLS2tmLRokWCV0m3WCzYtGkTTZhhIK77nhEqACfH/V999TFK\nBHp6emCz2VBZWRmWQZu0tDSkpaWhuroaGzZsgF6vh8FgwLFjx9Dd3Y3S0lKkp6eD4zhkZ2cjOzsb\nCQkJiI+PnzI+Qgi8Xi8sFguGh4cRCAQwPDyMpqamsdZsYmIiFAoFtm3bhuzs7BsSY2trK2Qy2ZwK\nBvPlwoULdPfJMAlp0mQY5hMA42uEMQAIgJ8TQj4I5bUpfgUCAWg0GsTExIQtYV4vMTERiYmJKC4u\nxpo1a+D3+2Gz2WAymdDS0oLe3l709PRgeHgYfX19yMvLg0wmg91uh1QqRWJiIliWxZUrV5CRkQG5\nXI5AIACGYVBaWorU1FSsWrUKpaWlyMjImLLOp8fjgdlsntXGZ6FgNBqRlpYmdBjzQkiTJiHk1ln8\nmBrA+G308q8+NqFnnnlm7OtNmzZh06ZNs7gkNZVAIID29nbIZDKUlZWJpmiwTCZDeno60tPTb6hZ\nyXEcOI6Dz+eDzWYDwzCIj49HTEwM3G43UlNT5zTANDAwgPLy8rBPK5oMy7K0UMcc1NfXo76+Pqhj\nBZ9yxDDMUYwMBJ29+v9KAHsA1GLktvwT0IEgwdhsNgwNDSEzM3PanRTnC4PBAIPBgIqKCtF8gHzj\nG9/A888/j/z8fKFDiQpTDQQJ1nvNMMyXGYYZALAGwIcMw9QBACGkDcA+AG0APgbwMM2MwnA6nejv\n70dSUhJNmFdxHIeGhoZpb9/DyeVyYenSpXOa9kUFT8jR83cBvDvJ954D8Fx4I6LGGxgYwPDwMBYv\nXizYZG0x0mg0qKqqElXBkfb2dvj9fsFH8OcL+ipT1+A4DlqtFh6PB8XFxTRhjuPxeGA0GrFgwYLp\nDw4jm80W0nJz1LXEOOWIEgjLsujq6gLHcSgvL6dz/q5z6tQpFBcXi66Y8vDwMN2yN4xoS5MCMNIv\n1tvbi9TUVFRUVNCEeR2z2QypVIrc3FyhQ7nB6dOnhQ5hXqFJk4Lb7UZ3dzfkcjlycnJo39h1WJbF\n4OAgli5dKrrXhhCC3NxcLFmyROhQ5g16ez7P6XQ66HQ6FBUViWJlixj19PQgOTlZlK+PWq2GWq2m\nfc9hJK6PTSqsTCYTLBYLCgoKRJkQxMDhcKC1tXXCrXjFwG633zCxnwot2tKchwghuHTpEtxuN6qq\nqkRXek1MBgcHsX79el73JufTaKUpij/nz5+f8vv03TLP+Hw+DA0NQaFQoKioiCbMKRiNRgBARkaG\nwJFMrqurC4sXLxY6jKhACMH777+Pt956a8rj6O35PDJaNDgQCEClUtGEOQWPx4OTJ0+Kbk7m9ZKT\nk7F8+XKhw4h4FosFjz32GOrr6/Hb3/52ymPpu2aeMJlMGBoaQkFBQVTtABkqWq0WpaWlSEhIEDqU\nSXk8HjQ2NuKee+4ROpSIptPp8Itf/AIOhwMvv/zytINqtKU5D9hsNuj1euTk5NCEGQSHwwG73Y6S\nkhKhQ5mSzWZDTU0NvWOYg+bmZjz22GP40pe+hL179wY1C4G+2lGut7cXw8PDqK6ujtidHsOJ4zhc\nvHgRZWVlop/gf+XKFaSkpAgdRkTiOA5//etfcfDgQfz0pz+dURcHTZpRanzRYJowgzc4OAi/34/U\n1FShQ5mWRqMR7ai+mHk8Hjz99NM4ceIE3n77bWRnZ0//Q+PQ2/MoFAgEcOnSJdhsNhQUFNCEGSSv\n14vh4WGsXr1a6FCCIpPJ6EqgGdLpdHj00UeRmpqKTz75ZMYJE6Atzahjs9nQ39+P3NxcWgNzhvr7\n+5GdnR0xrbf6+nrU1tYKHUbEOHv2LH7961/j61//Ou66665Z10OlSTOKOJ1ODA0NjW0BQQVPrVaj\nu7sb27ZtEzqUoIwWHqb7Ak2PEILf/e532LNnD3bv3j3n/elp0owSarV6rEAuXYc8MyzLwmw2Y926\ndaKpxj6d/v7+iIlVSB6PB7/97W/R19eHv/3tb6isrJzzOWmfZoTjOA5DQ0Pw+XyorKykCXMWhoaG\nkJiYGFHTsSwWS0QMVglJo9Hghz/8IRiGwYsvvshLwgRo0oxoLMuip6cHBoMBBQUFotkZMZIYjUZc\nvrEYMcMAAA/3SURBVHw54jYk83q9tFDHFF577TV87Wtfw/bt2/HUU0/x2k9Nk2aEcjqdaG9vR1JS\nEpYuXSr6OYViRAiBWq3GwoULI26C+NmzZ+nt+QQIIdi7dy9eeeUV/PCHP8RXv/pV3l+nyPpLoQCM\nFA1Wq9VISkoS1QZfkcZgMCA+Pj7iqgSxLIvCwkIUFhYKHYqouN1uPPvss9BoNPjLX/4yq+lEwaAt\nzQhjMBhw4cIFqFQqFBQUCB1OxBrdnjgSX0ONRgOr1UrvLsY5ffo07rvvPsTHx+PFF18MWcIEaNKM\nKBqNBna7HZWVlaIuJBEJOjo6IJVKERcXJ3QoM+Z2uyOuDzaUPv/8czz66KOorq7Gk08+GfJ5tvT2\nPAJwHIeBgQEYDAYsX7484vrfxMZisSAuLg4VFRVChzIrZrOZJk2MvC/+8pe/4M0338SLL74YthJ5\n9N0ncl6vF93d3UhLS8OKFStEt7FXpGFZFgMDA1i4cGHEvpa9vb3zvoamTqfDc889B5vNhr/85S9h\n7duPzL+aecLr9WJgYADx8fHIzc2N2De5mDQ1NSEQCCAxMVHoUGYtLS0t4gav+NTb24udO3dieHgY\nL730UtgHQ+m7UKSGh4dx9uxZZGVloaioSOhwooLL5YJEIuFtkrMQLBYLenp6IrIvlg9HjhzBd77z\nHfz0pz/FK6+8IkgxGnp7LkI6nQ5OpxPV1dVISkoSOpyoQAjBlStXUFBQEDEFOSbi9Xojti92Lvx+\nP/74xz/io48+wgsvvCBodSeaNEWEEAKNRoOBgQHU1NTQkm486uzshN/vj/hCJlqtdt59kNpsNtxx\nxx1ISkrCnj17BC9SQm/PRcLn86GtrQ0cx2HVqlU0YfLI5/PBYDBExWRwk8kU0S3lmerq6sLXv/51\nbN68GW+88YbgCROgLU1R8Pl8UKvVkMvlyMvLowM+POvv70dJSUlEFeSYTFxcXEROyJ+Njz76CM8/\n/zyeeOIJ3H777UKHM4YmTYE5HA6cP38e5eXlEX/rKEYajQYejwfFxcVChzJnPp8P58+fj/rCw36/\nH7t27YLRaMTu3btFV5iENmkEZLPZoNPpUFVVRRNmCLAsi6amJmRlZUVFcQuPx4NVq1ZF9fJJo9GI\n+++/f2zjM7ElTIAmTcFwHAebzYb8/Hy6o2CIqNVqVFdXR01RE41GI3QIIdXc3Iy7774bN910E954\n4w3RdkPQ23OBxMTE0KVwIeRwOGCxWFBVVSV0KLzxeDxROXLOcRyeeuop1NfX49lnn8UXvvAFoUOa\nEm1pUlGHEIJjx44hKysrqm5lA4FASKv3CMHj8eDJJ5/EgQMH8MILL4g+YQI0aVJRSKvVoqioCDk5\nOUKHwhtCCLq7u6NqupFarca9994Lt9uNEydORMwAF02aVFTxeDzQ6/VRMVo+ntfrRXl5edTsAbVv\n3z589atfxZ133okXXnghokod0j5NKqqcOnUKixcvjrrFATqdDm63W+gw5owQgj/84Q/4z//8T+ze\nvRu33Xab0CHNGE2aVNQwGAyQSqVR1+8HjCSbSJ8FYLfb8fjjj0On06GhoQF5eXlChzQr9Pacigp+\nvx9DQ0OoqamJyhVVdrs9olc0tba24o477kBGRgb27dsXsQkToEmTihKdnZ3IyMiImj6/60XymvPX\nX38dt956Kx544AE8++yzEfs8RtHbcyrimc1m9Pf3R23JtEAggMzMTCiVSqFDmZFAIIBf/epX+Oij\nj/Dmm29i3bp1QofEC5o0qYjGsizUajVuueWWqJqTOZ7ZbIbdbhc6jBnR6/V46KGHEBsbi/3794ui\nOhFfaNKkIpparUZiYmJEb18xnZiYGCxYsEDoMCbFsizq6urQ1NSEmpoacByHRx55BA899BD++Z//\nOeo+zGjSpCKW2WzGuXPnsGPHDqFDCSmz2SzalhrLsti6dSsaGxvhdDohl8shkUjw3HPP4ZFHHhE6\nvJCgA0FURCKEQKvV4uabb476LY19Pp9oB0/q6urQ2NgIh8MBQgi8Xi8ARPW+VoIlTYZhfsMwTDvD\nMM0Mw7zFMEzSuO89yTDMpavfj7zZr1TIabVayOVyZGZmCh1KSLEsi7i4ONEOAjU1NcHpdF7zmNvt\nRnNzs0ARhZ6QLc2DAKoIIcsBXALwJAAwDFMJ4B4AFQC2A/gfJhqKIVK8sVgsaGlpEW3pMD65XC7Y\nbDbR1gMtKyu7YWfMhISEqN6XXbCkSQg5RAjhrv63AcBonbQ7AbxGCAkQQvowklBXCxAiJVJDQ0Mo\nKSmJuqWSkxHjh8No90hFRQVqa2uhVCrBMAyUSiVqa2uxfft2oUMMGbF0Bt0P4NWrX6sAnBz3PfXV\nxygKRqMRMTExWLRokdChhMXw8DAUCoXQYVzD4/Ggp6cHsbGxqKiowKFDh1BXV4fm5mYsX74c27dv\nj7oR8/FCmjQZhvkEwPiFwAwAAuDnhJAPrh7zcwB+QsirE5yCosa4XC50dXWhpqZGtLerfJNIJKJa\n5eR0OtHW1oa4uDgUFRWN/R527NgR9bMYRoU0aRJCbp3q+wzDfBvA7QDGVx5VAxg/KS3/6mMTeuaZ\nZ8a+3rRpEzZt2jTzQKmIcPnyZSiVSlElkVALBAKieb5GoxFarRalpaURvQ5+IvX19aj//9u7v9C4\n0jKO498nMalhatjQ0jRp869Jqv2zdClhpVrKUlC7XuyuIMt6IS56V1FBEK17IXilghfLinerLKK7\nrAquguDuslRR+mcp9J/ZdoZCMklDpu2maZu2mZnMPF7MuAzdtJNpk3nnnPl9IDB5zyTzPHOSZ973\nnPO+5+jRFT3X3H1to7nfC5sdAn4JHHD3DyvadwK/Bz5LaVj+DjDqywRqZss1SwzNz89z+fJlduzY\nEcsFOZaTzWaZmppiZGQkaBxLS0tMTU2Ry+UYGBj42ImfODIz3H3Z4UzIY5qvAO3AO+Uu/nF3P+zu\n42b2JjAO5IHDqozNLZfLMTExwcjISNMUTCidbAl9SVU2m+XUqVN0dHSwZ8+epnr/7ydYT3M1qKfZ\nHE6fPk0ul+PJJ5vrIoqrV6/S2toabDbQzZs3SafTbNq0KfJredaqUXuaIlXdunWLlpYW9u7dGzqU\numtvbw8yFC4Wi2QyGa5fv87Q0FCkbkVRDyqa0rCKxSKTk5P09/fHfqrkcrLZbN0XIllaWuLEiRMk\nEgl2797dlO97NXpHpGGdOXOGQqHAY489FjqUuisWS/M+6nkMcWFhgXQ6TV9fH1u3btXxy/tQ0ZSG\ndPfuXdydXbt2hQ4liEKhUNf55levXiWTydDf309nZ2f1H2hiKprScNydiYkJ+vv7G+YaxXrL5XIf\n9TbXUrFY5OTJkyQSCbZv396wqyk1EhVNaTgXLlwgm82ycePG0KEE09bWtubD48XFRSYmJujq6mJ4\neFjHL1dI75I0lGw2y9zcHI8//njoUILK5/Nr2suen59nenqanp4eNmzYsGavE0cqmtJQ0uk0w8PD\nTX9craWlZU16msVikfPnz9Pa2sq2bdsabjGQKNDpMWkYk5OTZLNZuru7qz85xtydtra2Vf+9uVyO\nS5cukcvlGB0dVcF8SCqa0hDy+TzJZJLNmzc3zQpGD7Lavcxbt26RSqXo7OxkbGysadYiXQuaRikN\n4f/rM27ZoqVTV9ulS5dYXFxkYGCgYW+b0WgeNI1SPU0J7tq1a9y5c4eenp7QocTK0tISqVSK6elp\nRkdHVTBXiYqmBFUoFDh27Bjd3d2agbKKbt++TTKZJJFIcODAAQ3HV5GG5xJUOp2mWCwyODgYOpTY\nmJmZYW5ujr6+vtgtFlwvGp5LQ1pYWGB+fp6+vr7qT5aqisUiyWSS8fFxRkdHVTDXiIqmBFEoFDh6\n9Cg9PT2xvglXvSwuLpJMJuno6ODgwYOaDrmGdHG7BDE7O8vw8HDwlcnjYG5ujpmZGXp7e4MtWNxM\nVDSl7u7cucO1a9fYuXNn6FAizd25ePEiMzMz7Nu3r2kXN6k3Dc+lrorFIidOnKC3t3dNZr00i1wu\nRyqVYt26dezfv18Fs45UNFdopbf3jKp65XflyhUSiUSQYXlc9uHNmzdJJpN0dXUxNDT00eVEccnv\nQRohRxXNFWqEnbWW6pFfNpslk8mwZ8+eNX+t5UR9H7o74+PjnD17lsHBwY998EQ9v5VohBx1TFPq\n5ty5c/T19enM7kPI5/Ok02laW1sZGxtrinuPNyr1NKUuZmdnuXHjhs6WP4SFhQWSySTr169n+/bt\nKpiBRX5GUOgYRCSe7jcjKNJFU0Sk3jQ8FxGpgYqmiEgNVDSrMLNfmNkHZnbazP5sZp0V246YWaq8\n/Ysh43xYZvZVMztvZgUz23vPtsjnB2Bmh8zsgpklzeyHoeNZDWb2qpllzOxsRVuXmb1tZhfN7B9m\nFtkVO8xsq5m9Z2b/NbNzZvbdcnvwHFU0q3sb2OXuTwAp4AiAme0Engd2AE8Dv7Zo3qfhHPAV4J+V\njWa2gxjkZ2YtwK+ALwG7gK+Z2WfCRrUqfkspp0o/At51908D71H+W42oJeD77r4L2Ad8u7zfgueo\nolmFu7/r7sXyt8eBreXHzwBvuPuSu09QKqhPBgjxkbj7RXdPAfcWxGeJQX6UYk65+6S754E3KOUW\nae7+b+D6Pc3PAq+VH78GPFfXoFaRu8+6++ny4wXgA0r/e8FzVNGszTeBv5cfbwGmKrZdLrfFRVzy\nuzePaaKZx0pscvcMlIoOsClwPKvCzAaBJyh1WrpD56gZQYCZvQNU3jfWAAdecve/lZ/zEpB399cD\nhPhIVpKfxFLkryc0s/XAn4DvufvCMtdm1z1HFU3A3b/woO1m9iLwZeBgRfNloHLJ8a3ltoZTLb/7\niEx+VVwG+iu+j2oeK5Exs253z5jZZuBK6IAehZl9glLB/J27v1VuDp6jhudVmNkh4AfAM+6erdj0\nV+AFM2s3syFgBDgZIsZVVHlcMy75vQ+MmNmAmbUDL1DKLQ6Mj++zF8uPvwG8de8PRMxvgHF3f7mi\nLXiOmhFUhZmlgHbgw3LTcXc/XN52BPgWkKc0fHg7TJQPz8yeA14BNgLzwGl3f7q8LfL5wUcffC9T\n6iS86u4/CxzSIzOzPwBPARuADPAT4C/AHymNECaB5919PlSMj8LMPg/8i9LVHV7++jGlD+43CZij\niqaISA00PBcRqYGKpohIDVQ0RURqoKIpIlIDFU0RkRqoaIqI1EBFU0SkBiqaIiI1UNGUWDOzMTM7\nU54OmigvuLwzdFwSXZoRJLFnZj8FOspfU+7+88AhSYSpaErsmVkbpYU77gKfc/3RyyPQ8FyawUZg\nPfAp4JOBY5GIU09TYs/M3gJeB4aAXnf/TuCQJMK0CLHEmpl9Hci5+xvlm6z9x8yecvejgUOTiFJP\nU0SkBjqmKSJSAxVNEZEaqGiKiNRARVNEpAYqmiIiNVDRFBGpgYqmiEgNVDRFRGrwPxk2rGb5nC2e\nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10806fba8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import rebound\n",
    "sim = rebound.Simulation()\n",
    "k = 0.01720209895 # Gaussian constant\n",
    "sim.G = k**2\n",
    "sim.add(m=1.) # Sun\n",
    "sim.add(m=1.e-3, a=5.) # Jupiter\n",
    "sim.add(m=3.e-4, a=10.) # Saturn\n",
    "for comet_elem in comets:\n",
    "    sim.t = comet_elem[-1] # last column is the time of observation\n",
    "    o = getOrbit(comet_elem)\n",
    "    sim.add(a=o.a, e=o.e, inc=o.inc, Omega=o.Omega, omega=o.omega, T=o.T)\n",
    "fig = rebound.OrbitPlot(sim, trails=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}