{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# WHFast tutorial\n", "\n", "This tutorial is an introduction to the python interface of WHFast, a fast and unbiased symplectic Wisdom-Holman integrator. This integrator is well suited for integrations of planetary systems in which the planets stay roughly on their orbits. If close encounters and collisions occur, then WHFast is not the right integrator. The WHFast method is described in detail in Rein & Tamayo (2015).\n", "\n", "This tutorial assumes that you have already installed REBOUND." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**First WHFast integration**\n", "\n", "You can enter all the commands below into a file and execute it all at once, or open an interactive shell).\n", "\n", "First, we need to import the REBOUND module (make sure you have enabled the virtual environment if you used it to install REBOUND)." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import rebound" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we create a REBOUND simulation instance. This object encapsulated all the variables and functions that REBOUND has to offer. " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "sim = rebound.Simulation()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, we can add particles. We'll work in units in which $G=1$ (see [Units.ipynb](../Units) for using different units). The first particle we add is the central object. We place it at rest at the origin and use the convention of setting the mass of the central object $M_*$ to 1:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "scrolled": false }, "outputs": [], "source": [ "sim.add(m=1.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's look at the particle we just added:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "print(sim.particles[0])" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "The output tells us that the mass of the particle is 1 and all coordinates are zero. \n", "\n", "The next particle we're adding is a planet. We'll use Cartesian coordinates to initialize it. Any coordinate that we do not specify in the `sim.add()` command is assumed to be 0. We place our planet on a circular orbit at $a=1$ and give it a mass of $10^{-3}$ times that of the central star." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "sim.add(m=1e-3, x=1., vy=1.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Instead of initializing the particle with Cartesian coordinates, we can also use orbital elements. By default, REBOUND (as well as WHFast internally) will use Jacobi coordinates, i.e. REBOUND assumes the orbital elements describe the particle's orbit around the centre of mass of all particles added previously. Our second planet will have a mass of $10^{-3}$, a semimajoraxis of $a=2$ and an eccentricity of $e=0.1$ (note that you shouldn't change G after adding particles this way, see [Units.ipynb](../Units)):" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "sim.add(m=1e-3, a=2., e=0.1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have added two more particles, let's have a quick look at what's in this simulation by using" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "---------------------------------\n", "REBOUND version: \t3.4.0\n", "REBOUND built on: \tMay 31 2017 11:53:50\n", "Number of particles: \t3\n", "Selected integrator: \tias15\n", "Simulation time: \t0.0000000000000000e+00\n", "Current timestep: \t0.001000\n", "---------------------------------\n", "\n", "\n", "\n", "---------------------------------\n" ] } ], "source": [ "sim.status()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can see that REBOUND used the `ias15` integrator as a default. Next, let's tell REBOUND that we want to use `WHFast` instead. We'll also set the timestep. In our system of units, an orbit at $a=1$ has an orbital period of $T_{\\rm orb} =2\\pi \\sqrt{\\frac{a^3}{GM}}= 2\\pi$. So a reasonable timestep to start with would be $dt=10^{-3}$ (see Rein & Tamayo 2015 for some discussion on timestep choices)." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true }, "outputs": [], "source": [ "sim.integrator = \"whfast\"\n", "sim.dt = 1e-3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`whfast` refers to the 2nd order symplectic integrator WHFast described by Rein & Tamayo (2015). By default, no symplectic correctors are used, but they can be easily turned on (see [Advanced Settings for WHFast](../AdvWHFast)). \n", "\n", "We are now ready to start the integration. Let's integrate the simulation for one orbit, i.e. until $t=2\\pi$. Because we use a fixed timestep, rebound would have to change it to integrate exactly up to $2\\pi$. \n", "\n", "**Note: The default is for sim.integrate to simulate up to exactly the time you specify. This means that in general it has to change the timestep close to the output time to match things up. A changing timestep breaks WHFast's symplectic nature, so when using WHFast, you typically want to pass the `exact_finish_time = 0` flag, which will instead integrate up to the timestep which is nearest to the endtime that you have passed to `sim.integrate`.**" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "sim.integrate(6.28318530717959, exact_finish_time=0) # 6.28318530717959 is 2*pi" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Once again, let's look at what REBOUND's status is" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "---------------------------------\n", "REBOUND version: \t3.4.0\n", "REBOUND built on: \tMay 31 2017 11:53:50\n", "Number of particles: \t3\n", "Selected integrator: \twhfast\n", "Simulation time: \t6.2839999999992369e+00\n", "Current timestep: \t0.001000\n", "---------------------------------\n", "\n", "\n", "\n", "---------------------------------\n" ] } ], "source": [ "sim.status()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can see the time has advanced to $t=2\\pi$ and the positions and velocities of *all* particles have changed. If you want to post-process the particle data, you can access it in the following way:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.003326154866766361 0.009674635911450204 0.0005194654213133196 0.0012200269278386914\n", "1.0032694180883746 0.0366289427053581 -0.024395944012027243 0.9999782071644221\n", "-1.5284252838557046 1.496351615582015 -0.4950694773013606 -0.4364888285516785\n" ] } ], "source": [ "particles = sim.particles\n", "for p in particles:\n", " print(p.x, p.y, p.vx, p.vy)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `particles` object is an array of pointers to the particles. This means you can call `particles = sim.particles` before the integration and the contents of `particles` will be updated after the integration. If you add or remove particles, you'll need to call `sim.particles` again." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Visualization with matplotlib**\n", "\n", "Instead of just printing boring numbers at the end of the simulation, let's visualize the orbit using matplotlib (you'll need to install numpy and matplotlib to run this example, see [Installation](../Installation)).\n", "\n", "We'll use the same particles as above. As the particles are already in memory, we don't need to add them again. Let us plot the position of the inner planet at 100 steps during its orbit. First, we'll import numpy and create an array of times for which we want to have an output (here, from $T_{\\rm orb}$ to $2 T_{\\rm orb}$ (we have already advanced the simulation time to $t=2\\pi$)." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np\n", "torb = 2.*np.pi\n", "Noutputs = 100\n", "times = np.linspace(torb, 2.*torb, Noutputs)\n", "x = np.zeros(Noutputs)\n", "y = np.zeros(Noutputs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we'll step through the simulation. Rebound will integrate up to `time`. Depending on the timestep, it might overshoot slightly. If you want to have the outputs at exactly the time you specify, you can set the `exact_finish_time=1` flag in the `integrate` function (or omit it altogether, 1 is the default). However, note that changing the timestep in a symplectic integrator could have negative impacts on its properties." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "for i,time in enumerate(times):\n", " sim.integrate(time, exact_finish_time=0)\n", " x[i] = particles[1].x\n", " y[i] = particles[1].y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's plot the orbit using matplotlib." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUwAAAEzCAYAAABaGjpLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VPW9//HXJ3uAJOxJIAn7DioaFoG61KVIrVitFe2i\ndUFb7Xa7edv7895fa/fb9tdWrxaXutxWrbZUqihq1aIiSED2RUPYspAEAgkhe/L9/ZHBppCQQ2Yy\nZyZ5Px+PeeTMmW/m+8lJ8p6zfc8x5xwiItK5GL8LEBGJFgpMERGPFJgiIh4pMEVEPFJgioh4pMAU\nEfEo6MA0s2wze93MtpnZVjP7ajttzMx+Y2b5ZrbJzM4Otl8RkXCLC8F7NAHfcM6tN7MUYJ2ZveKc\n29amzWXAuMBjFnB/4KuISNQIeg3TOVfinFsfmD4KbAeGn9BsIfC4a7Ua6G9mmcH2LSISTiHdh2lm\nI4HpwJoTXhoO7G/zvJCTQ1VEJKKFYpMcADPrB/wZ+JpzriqI91kMLAbo27fvORMnTgxRhSIirdat\nW3fQOTfkdL8vJIFpZvG0huUfnHN/aadJEZDd5nlWYN5JnHNLgCUAubm5Li8vLxQlioh8yMz2duX7\nQnGU3ICHge3OuV920GwZ8PnA0fLZQKVzriTYvkVEwikUa5hzgc8Bm81sQ2Ded4EcAOfcA8ByYAGQ\nD9QAXwhBvyIiYRV0YDrn3gKskzYOuCPYvkRE/KSRPiIiHikwRUQ8UmCKiHikwBQR8UiBKSLikQJT\nRMQjBaaIiEcKTBERjxSYIiIeKTBFRDxSYIqIeKTAFBHxSIEpIuKRAlNExCMFpoiIRwpMERGPFJgi\nIh4pMEVEPFJgioh4pMAUEfFIgSki4pECU0TEIwWmiIhHCkwREY8UmCIiHoUkMM3sETMrM7MtHbx+\ngZlVmtmGwOPuUPQrIhJOcSF6n0eBe4HHT9HmTefc5SHqT0Qk7EKyhumcWwlUhOK9REQiVTj3YZ5r\nZhvN7EUzmxLGfkVEQiJUm+SdWQ+McM5Vm9kC4K/AuPYamtliYDFATk5OmMoTEelcWNYwnXNVzrnq\nwPRyIN7MBnfQdolzLtc5lztkyJBwlCci4klYAtPMMszMAtMzA/0eCkffIiKhEpJNcjN7ErgAGGxm\nhcB/AvEAzrkHgE8BXzSzJqAWWOScc6HoW0QkXEISmM656zp5/V5aTzsSEYlaGukjIuKRAlNExCMF\npoiIRwpMERGPFJgiIh4pMEVEPFJgioh4pMAUEfFIgSki4pECU0TEIwWmiIhHCkwREY8UmCIiHikw\nRUQ8UmCKiHikwBQR8UiBKSLikQJTRMQjBaaIiEcKTBERjxSYIiIeKTBFRDxSYIqIeKTAFBHxSIEp\nIuJRSALTzB4xszIz29LB62ZmvzGzfDPbZGZnh6JfEZFwCtUa5qPA/FO8fhkwLvBYDNwfon5FRMIm\nLhRv4pxbaWYjT9FkIfC4c84Bq82sv5llOudKQtG/9Gx1jc0UHq6hpLKOIzWNHKlt5MixBiprG2lq\ncbS44w+IizH6JcbRLymOlMQ4UpLiyUhLYnj/ZNJTk0iI014o6bqQBKYHw4H9bZ4XBuYpMOVDtQ3N\nbD9QxdaiSrYUVVFwsJp9FTWUVtW32z45PpaEuBhiDGLMMIOGphaONTTT3OJOam8GQ1MSGTW4LxMz\nUpmUmcLEjFTGp6eQnBDb3T+e9ADhCkzPzGwxrZvt5OTk+FyNdKfK2kZWFxxiVf5B1uyu4P3SoxzP\nuYF9Exg7tB/njRtCzsA+ZA/sw7D+yQzsG09acgJpyfEdri0656hvauFoXROVtY0cqKyjuLKW4iO1\nFB2uJb+8mj/l7aemoRloXSudlpXGrFGDmDV6ILkjBpCSFB+uxSBRJFyBWQRkt3meFZh3EufcEmAJ\nQG5u7smrCRLV8suqeXFzCa9uL2VzUSUtrnVNMXfkAC6dnM7U4WlMHZ5GZloSZtalPsyMpPhYkuJj\nGZKSyNih/U5q09Li2H+4hu0lR9lYeIR3d1fw8FsFPPCPXcTGGDNHDuSSyelcMjmd7IF9gv2xpYew\n1t2KIXij1n2Yzzvnprbz2seBO4EFwCzgN865mZ29Z25ursvLywtJfeKfgvJqlm0sZvnmEt4vrQbg\n7Jz+zBs3hLljBjE9Z0BE7FusbWhm/b7DvJ1/kFe3l35Y66TMVBaeNYyrpg9naGqSz1VKKJjZOudc\n7ml/XygC08yeBC4ABgOlwH8C8QDOuQesdVXhXlqPpNcAX3DOdZqECszoVd/UzIqtpfxxzV5WF1Rg\nBjNGDmTB1AzmT80kIy3yg2fPwWO8sq2UF7eUsH7fEWJjjAsnDOGa3Gw+OnEo8bH+h7x0ja+B2V0U\nmNGn/Gg9v397N0+t3U/FsQayByazaEYOnzoni/QoXjvbVV7Ns+sK+fO6QsqO1pOZlsRNc0exaGa2\n9ndGIQWm+KrwcA0PrizgqbX7aWhu4dLJ6Xxm1gjmjR1MTEzX9kVGoqbmFt7YWc7Db+3mnYJDpCTG\ncd2sHG6eNyqqPxB6GwWm+KK0qo5fvfI+z64rxAyump7FbeePZvSQkw+09DSbCo/w4Ju7Wb65hPhY\n48Y5o/ji+WNI66M1zkinwJSwOlbfxJKVBSxZWUBTSwufmTWCxeeNZlj/ZL9LC7t9h2r45Ss7eW5j\nMSmJcXzxgrF8Ye5IkuJ1bmekUmBKWDjnWPpeET9+cQflR+v5+BmZfOdjE8kZpFNvthVX8fMVO3h9\nZzk5A/vw/YVTuGDCUL/LknYoMKXb7a+o4btLN/PmBwc5K7s//+fyyZwzYoDfZUWcVfkH+Y/ntlBQ\nfowF0zK4+/IpUXFWQG+iwJRu09zi+P3bu/nFy+8TY3DXZRP5zKwRPepgTqjVNzXz4MoCfvtaPvGx\nMfznJybzqXOyunwyvoSWAlO6RWlVHV97agPvFBzioolD+cGVU3vlfsqu2neohm89u5E1uytYMC2D\nH145jQF9E/wuq9framBG3FhyiRyv7Sjlm89sorahmZ9/6gytIXVBzqA+/PHW2Tz4ZgG/eHkn6/Ye\n5pefPou5Ywf7XZp0gYYqyEmaWxw/Xr6dmx7NIz01ib99eR7X5GYrLLsoNsa4/fwxLP3SXPolxvG5\nh9ewZOUuInnrTtqnwJR/cbSukVseW8vvVhbw2dk5LP3SnHYvXiGnb+rwNJbdOY/5UzP40fIdfPWp\nDdQGrpgk0UGb5PKh/RU13PzYWnaVH+OeK6fy2dkj/C6px+mbGMd915/N/7yxi/9+eSf5ZdU8cuMM\nHUWPElrDFAA2F1ay8L63Ka2q54mbZiosu5GZcceFY3nkhhnsPXSMq+9fRUF5td9liQcKTGHtngqu\nf3A1yfGxLP3SHObogERYXDhxKE8tPpe6xmaueeAdNhdW+l2SdEKB2cu9+UE5n3/4XYakJPLM7ef2\nijHgkWRaVhrP3H4uSfGxLFryDmsKDvldkpyCArMX+8f75dz8aB4jBvXh6dvO1fmVPhk9pB9//uIc\nMtKSuOnRtazfd9jvkqQDCsxeat3eCm57Io+xQ/vx1OLZDElJ9LukXi0jLYk/3jqbwSmJ3PDIu2wp\n0uZ5JFJg9kLbS6r4wu/XkpGaxGM3zaR/H408iQTpqUn84ZZZpCbF87mH1/BB6VG/S5ITKDB7mcLD\nNXz+kXfpkxDHEzfP0pplhMka0Ic/3jqLuNgYvvDoWg5Wt3+LYfGHArMXqWlo4tbH11HX2MwTN8/U\n3RAj1IhBfXno87kcrK5n8eN51DXq5PZIocDsJZxzfOuZTew8UMVvr5vOuPQUv0uSUzgzuz+//PRZ\nrN93hG8/u0nDKCOEArOXuO/1fF7YXMJdl03URW2jxIJpmXzrYxNYtrGYh97c7Xc5ggKzV1i16yC/\neOV9rjxrGLd+ZLTf5chp+NIFY/jYlHR++tIONuw/4nc5vZ4Cs4c7UtPAvz29kVGD+vKjq6bpikNR\nxsz42dVnkp6axJ1/XE9lbaPfJfVqCswezDnHd5du5mB1Pb9eNJ0+CbrWSjRK6xPPb6+fzoHKOr67\ndLPf5fRqCswe7M/ri1i++QDfuHQC07LS/C5HgnB2zgC+fsl4XthUwoubS/wup9cKSWCa2Xwz22lm\n+WZ2Vzuv32hm5Wa2IfC4JRT9SscqjjVwzwvbyB0xgMXnab9lT3DbeaOZOjyV//PcViprtGnuh6AD\n08xigfuAy4DJwHVmNrmdpk87584KPB4Ktl85tR8v3051XRM/umoasbpZWY8QFxvDT68+g8M1rR+G\nEn6hWMOcCeQ75wqccw3AU8DCELyvdNGagkM8s66QW88bzXidb9mjTBmWxm3njeaZdYW8s0tXNgq3\nUATmcGB/m+eFgXknutrMNpnZs2aWHYJ+pR3NLY67n9tK1oBkvvLRcX6XI93gKxeNY1haEj9cvo2W\nFp3QHk7hOujzN2Ckc+4M4BXgsY4amtliM8szs7zy8vIwlddzLH2viJ2lR/nugkkkJ8T6XY50g6T4\nWL49fyJbiqpY+l6R3+X0KqEIzCKg7RpjVmDeh5xzh5xzx68i8BBwTkdv5pxb4pzLdc7lDhkyJATl\n9R71Tc386pX3OSMrjcumZvhdjnSjK84cxplZafx8xU7dSC2MQhGYa4FxZjbKzBKARcCytg3MLLPN\n0yuA7SHoV07wv6v3UXSklu/Mn6gT1Hu4mBjjPy6fzIGqOp5YvcfvcnqNoAPTOdcE3AmsoDUI/+Sc\n22pm3zezKwLNvmJmW81sI/AV4MZg+5V/VdfYzP1v5DNv7GDm6p48vcKMkQOZO3YQD765W1c0CpOQ\n7MN0zi13zo13zo1xzv0wMO9u59yywPS/O+emOOfOdM5d6JzbEYp+5Z+WvlfEweoG7rhwrN+lSBjd\nccFYyo/W8+y6Qr9L6RU00qcHaGlxPPRmAVOHpzJ79EC/y5EwOnfMIM7K7s8D/9hFU3OL3+X0eArM\nHuD1nWXsKj/GrR8ZrX2XvYyZ8aULxlB4uJaXt5X6XU6Pp8DsAR5dtYdhaUksmJbZeWPpcS6alM7w\n/sk8+e4+v0vp8RSYUa74SC1v5R/kmtxs4mP16+yNYmOMa2dk8+YHB9l76Jjf5fRo+g+LckvfK8I5\nuPrsLL9LER99Ojeb2BjjyXf3d95YukyBGcWcczy7rpBZowaSM0g3NOvNMtKSuHDCUJa+V6jhkt1I\ngRnFNhVWsvvgMa1dCgCfODOT0qp63tt/2O9SeiwFZhR7dXspMQaXTE73uxSJAB+dOJSE2BiWbz7g\ndyk9lgIzir2yrZTcEQMZ0DfB71IkAqQkxXPe+MG8uLlEt+XtJgrMKFV4uIYdB45y8WTdMlf+af7U\nTIor69haXOV3KT2SAjNKvb6z9dJ3H52ozXH5p4+Ma72OgC4u3D0UmFFq7e4K0lMTGTOkr9+lSARJ\nT01i9JC+rNp10O9SeiQFZpTK21NB7siBGgopJ5kzZhDv7q6gUWPLQ06BGYWKj9RSXFlH7ogBfpci\nEWjOmMEca2hmS1Gl36X0OArMKLRub+t5drkjdGUiOdmZ2f0B2KIDPyGnwIxCOw5UERdjTMjQHSHl\nZMPSkkhLjmdbsdYwQ02BGYU+KK1m5OC+JMTp1ycnMzOmDEtlm9YwQ07/cVEov6ya8en9/C5DItjk\nzFR2HDiqiwqHmAIzytQ3NbPn0DHGDtXmuHRs9JB+1De1UHq0vvPG4pkCM8oUHa6lxcFIXZ1ITmH4\ngGSg9e9FQkeBGWXKAmsM6alJPlcikSwrEJiFh2t8rqRnUWBGmeOBOTQl0edKJJIN7681zO6gwIwy\nZVV1AAxRYMopJMXHkpYcT3m19mGGkgIzylQcayAuxkhLjve7FIlwlbWNPP7OXr/L6FEUmFGmrrGF\n5PhYjSEX8UFIAtPM5pvZTjPLN7O72nk90cyeDry+xsxGhqLf3qi+qZnEeH3Oifgh6P88M4sF7gMu\nAyYD15nZ5BOa3Qwcds6NBX4F/DTYfnur+qYWEuNi/S5DpFcKxarKTCDfOVfgnGsAngIWntBmIfBY\nYPpZ4CLTNmWXNDS1EB+rRSedGzW4L+mpOjgYSqEIzOFA25shFwbmtdvGOdcEVAKDQtB3rxMXYzTr\nfi3iQdaAZIYFTi+S0Ii4nWFmttjM8swsr7y83O9yIk5ifAx1jRofLJ2raWimb0Kc32X0KKEIzCIg\nu83zrMC8dtuYWRyQBrR70xHn3BLnXK5zLnfIkCEhKK9nSYyLpa6x2e8yJAocq2+iT4L2d4dSKAJz\nLTDOzEaZWQKwCFh2QptlwA2B6U8BrzndB7RLkuJjqdcapnhwrKGJfolawwyloJemc67JzO4EVgCx\nwCPOua1m9n0gzzm3DHgYeMLM8oEKWkNVuiA1OY6G5hZqGproo80tOYXquib6KjBDKiRL0zm3HFh+\nwry720zXAdeEoq/ebmhK60U3yqrqGTlY/wzSvrrGZg7XNOooeYhF3EEfObXj/wBlus6hnEJJZes1\nBzLTdJQ8lBSYUeb4GmZp4CIcIu0pOdJ6laLM/roMYCgpMKPM8X+AQl22S06hOLCGOUxrmCGlwIwy\nqUnxDE1JJL+s2u9SJILll1UTH2s6cT3EFJhRaHx6Ch+UHfW7DIlgOw5UMXZoiu4sGmJamlFoXHo/\n8suqaWnRqazSvu0lVUzSfetDToEZhSakp1DT0Mx+3a9F2nH4WAOlVfVMzFRghpoCMwqdldMfgHV7\nD/tciUSiDfuPADB1WJrPlfQ8CswoNH5oCilJcazdo8CUk60uOERCbAzTcwb4XUqPo8CMQjExRu6I\nAazdU+F3KRKBVu+u4MzsNJJ14Y2QU2BGqdyRA8kvq+ag7goobVTXN7GlqJLZo3W52e6gwIxS541r\nvfTdGzt1zVD5p7fzD9Lc4jh3jAKzOygwo9TU4amkpyby6rZSv0uRCLJi6wHSkuOZMXKg36X0SArM\nKGVmXDwpnZUflOuCwgJAY3MLr24r5eJJ6cTH6l+7O2ipRrGLJ6dT09DMql0H/S5FIsCaggqq6pr4\n2JR0v0vpsRSYUWzOmEGkJcez9L1iv0uRCPC3jcX0SYjlvPG6tUt3UWBGscS4WBaeNYwVWw9QWdvo\ndznio+r6Jv62qZjLz8gkKV6nE3UXBWaUu+acbBqaWvjbRq1l9mbPbyympqGZRTNz/C6lR1NgRrmp\nw1OZkJ7CM3n7O28sPdZTa/czPr0f07P7+11Kj6bAjHJmxvWzcthYWMm6vRr50xttKjzChv1HuHZG\nDmbmdzk9mgKzB7gmN4u05HiWrCzwuxTxwQP/2EVKUhyfzs3yu5QeT4HZA/RJiONzs0fw8rZSdh88\n5nc5EkYF5dW8uOUAnz93BClJ8X6X0+MpMHuIG+aMJD4mRmuZvcySlQXEx8Zw45xRfpfSKygwe4gh\nKYlcOyObZ/L2s0drmb3C3kPH+PP6Qj6dm8WQFN1/PBwUmD3Ily8aS3xsDL945X2/S5Ew+PmKncTF\nxPDlj47zu5ReI6jANLOBZvaKmX0Q+NruFUvNrNnMNgQey4LpUzo2NCWJm+eN4m8bi9lSVOl3OdKN\nNuw/wvObSrj1I6NIT9W9x8Ml2DXMu4C/O+fGAX8PPG9PrXPurMDjiiD7lFNYfP5oBvSJ54cvbMc5\n3SStJ3LO8aPl2xncL4HF54/xu5xeJdjAXAg8Fph+DLgyyPeTIKUmxfNvl07gnYJDPLdBo396or+s\nL+Ld3RV8/ZLx9EuM87ucXiXYwEx3zpUEpg8AHV0mJcnM8sxstZkpVLvZZ2bmcFZ2f37w/DaO1DT4\nXY6E0KHqeu55YRvnjBjAdTM0DDLcOg1MM3vVzLa081jYtp1r3f7raBtwhHMuF7ge+H9m1uF2hJkt\nDoRrXnm5ribeFTExxo+vmsaR2kZ++tIOv8uREPrB89uorm/ix1dNIyZGo3rCrdPAdM5d7Jyb2s7j\nOaDUzDIBAl/LOniPosDXAuANYPop+lvinMt1zuUOGaLLVHXVpMxUbpk3iiff3c/K9/XB0xO8tqOU\nv24o5ovnj2F8uu457odgN8mXATcEpm8AnjuxgZkNMLPEwPRgYC6wLch+xYOvXzKe8en9+MYzGzmk\nm6VFtbKqOr71zCYmZqTwpQvH+l1OrxVsYP4EuMTMPgAuDjzHzHLN7KFAm0lAnpltBF4HfuKcU2CG\nQVJ8LL9eNJ3K2ka+/ewmHTWPUi0tjq//aQPHGpq49/rput6lj4I6xOacOwRc1M78POCWwPQqYFow\n/UjXTcpM5a75E/n+89t4dNUevjBXQ+iizf3/2MXb+Yf4yVXTGDtUm+J+0kifXuALc0dy8aSh/PCF\n7br/T5R5Y2cZv3h5Jx8/I5NrZ2T7XU6vp8DsBcyMX117FiMH9+WOP6xnf0WN3yWJB/llR/nyH99j\nQkYqP7v6DF3rMgIoMHuJlKR4Hvp8Li0Obnksj+r6Jr9LklM4fKyBmx/LIzE+hoduyKWvTlCPCArM\nXmTk4L7cd/3Z5JdXc/sT66hv0v3MI1FtQzOLn8ij5Egdv/tcLsP7J/tdkgQoMHuZeeMG87Orz+Ct\n/IN89ckNNDW3+F2StFHf1Mxt/7uOvL2H+cWnz+ScEe1ez0Z8osDsha4+J4u7L5/MS1sP8O9/2azT\njSJEU3MLX31yAyvfL+cnV03jE2cO87skOYF2jPRSN80bRWVtI7/++wfExRr3XDmNWA21801jcwvf\nfGYjL209wN2XT+ZajROPSArMXuxrF4+jxTl++1o+1fXN/PLTZxIfq42OcKtrbOaOP6zn7zvK+M78\nidw0T+fKRioFZi9mZnzj0gn0TYzjJy/uoKa+ifs+c7ZGkoRRVV0jtzyWx9o9Fdxz5VQ+O3uE3yXJ\nKWh1Qrj9/DHcc+VUXttZxmceWkP5UY07D4fiI7Us+t1q1u89zK8XTVdYRgEFpgDw2dkjuO/6s9la\nXMmV973NtuIqv0vq0fL2VHDFvW+xr6KGh27I5Qod4IkKCkz50IJpmTx7+xyaWxxX37+Kl7aUdP5N\nctqeXruP6x5cTb/EOP56xxwumDDU75LEIwWm/Iupw9NYdudcJmSkcPv/rue/lm2lrlEnuIfCsfom\nvvXMRr7z583MHj2I5+6Yp4tpRBkFppxkaGoST982m5vmjuLRVXv45P+sIr+s2u+yotrG/Uf4+G/e\n5Nn1hdx54Vh+f+MM0vrE+12WnCYFprQrMS6Wuz8xmUduzKW0qo5P/PYtHlu1h+YWneR+OhqbW7jv\n9Xyuvn8V9U0tPHnrbL75sQnE6fStqGSRPMojNzfX5eXl+V1Gr1dWVcc3n93EyvfLmZ7Tn59cdQYT\nMrQp2Zn1+w7z3b9sZseBoyyYlsGPP3mG1iojhJmtC9xn7PS+T4EpXjjn+OuGIn7w/Haqahu57fzR\n3HHhWPok6FTeE1XWNvLzFTv4w5p9pKck8X8XTuFjUzL8Lkva6Gpg6q9dPDEzPjk9i/PHD+WeF7Zx\n3+u7eCavkH+7ZDzX5GZrWCWtF8544p293Pt6PlW1jdw4ZyTfuHSC7h3eg2gNU7okb08FP1q+nfX7\njjA+vR/fmT+Rj04c2isvctvS4li2sZj/fnknhYdrmTd2MHddNpGpw9P8Lk06oE1yCTvnHC9tOcBP\nX9rBnkM1TM5M5YsXjGHBtMxescZZ39TMc+8V87uVu9hVfowpw1K567KJfGScbg8d6RSY4puGphb+\nuqGIB/6xi4LyY4wY1IdbPzKaK6cP75Gbo5U1jTy1dh+PvL2b0qp6JmemcvsFY7h8WiYxveCDoidQ\nYIrvWlocL28r5f438tlYWEnfhFiuOGsYi2bkcEZWWlRvrjvnWF1QwdNr9/HilgPUN7UwZ8wgbj9/\nDB8ZNziqf7beSIEpEcM5x3v7j/Dkmn08v6mE2sZmJmak8IkzhzF/agZjhvTzu0RPnHNsLznKiq0H\neG5DEXsO1ZCSFMeVZw3n2hnZ2kcZxRSYEpGq6hp5bkMxf1lfyHv7jgAwPr0f86dmcv74IZyRlRZR\n1+Csb2pmw74j/H1HGS9tOcC+ihrMYNaogVw7I5v5UzJJTtDl76KdAlMiXkllLSu2HODFLQdYu6eC\nFgd9E2KZOWogc8YM5uwR/ZmUmRrWczur6hrZVlzFu7srWF1wiHV7D1Pf1EJ8rDF37GDmT8ng4snp\nDO6XGLaapPv5Ephmdg3wX8AkYKZzrt10M7P5wK+BWOAh59xPvLy/ArPnqjjWwJqCQ7y96yCrdh2i\noPwYAGYwanBfpgxLY2JGCtkD+5A9IJnsgX0Y1DehS/sKW1oc5dX1FB6uofBwLXsO1rCtpJLtJUfZ\nF7hHuxlMykhl9uhBnDtmELNGDyQ1SaNyeiq/AnMS0AL8Dvhme4FpZrHA+8AlQCGwFrjOObets/dX\nYPYepVV1bCqsZGtxJVuLq9hWXEXRkdp/aZMUH8PAPgmkJseTmhxPWnI8yfGxHM9QA1ocVNc3cbSu\nkaN1TRyta6L8aD0Nbe6OaQYjB/VlcmYqk4elMikzhbNzBtC/T0IYf2Lxky8jfZxz2wOdn6rZTCDf\nOVcQaPsUsBDoNDCl90hPTeKSyUlcMjn9w3nV9U2ta4UVtew/XEPR4VoO1zRSVddIZW0j+ytqPrz0\n3PGPfQP6JcWRkhhPzsA+pCTFMyQlkeEDksnqn0zWgGSGD0jWkE7pknD81QwH9rd5XgjMCkO/EuX6\nJcYxMSOViRmpfpciAngITDN7FWjvygHfc849F+qCzGwxsBggJ0e3GhWRyNFpYDrnLg6yjyIgu83z\nrMC8jvpbAiyB1n2YQfYtIhIy4TgBbi0wzsxGmVkCsAhYFoZ+RURCKqjANLNPmlkhcC7wgpmtCMwf\nZmbLAZxzTcCdwApgO/An59zW4MoWEQm/YI+SLwWWtjO/GFjQ5vlyYHkwfYmI+C1yxqSJiEQ4BaaI\niEcKTBERjxSYIiIeKTBFRDxSYIqIeKTAFBHxSIEpIuKRAlNExCMFpoiIRwpMERGPFJgiIh4pMEVE\nPFJgioh4pMAUEfFIgSki4pECU0TEIwWmiIhHCkwREY8UmCIiHikwRUQ8UmCKiHikwBQR8UiBKSLi\nkQJTRMRiezKkAAAG/klEQVSjoALTzK4xs61m1mJmuadot8fMNpvZBjPLC6ZPERG/xAX5/VuAq4Df\neWh7oXPuYJD9iYj4JqjAdM5tBzCz0FQjIhLBwrUP0wEvm9k6M1scpj5FREKq0zVMM3sVyGjnpe85\n557z2M8851yRmQ0FXjGzHc65lR30txhYDJCTk+Px7UVEul+ngemcuzjYTpxzRYGvZWa2FJgJtBuY\nzrklwBKA3NxcF2zfIiKh0u2b5GbW18xSjk8Dl9J6sEhEJKoEe1rRJ82sEDgXeMHMVgTmDzOz5YFm\n6cBbZrYReBd4wTn3UjD9ioj4Idij5EuBpe3MLwYWBKYLgDOD6UdEJBJopI+IiEcKTBERjxSYIiIe\nKTBFRDxSYIqIeKTAFBHxSIEpIuKRAlNExCMFpoiIRwpMERGPFJgiIh4pMEVEPFJgioh4pMAUEfFI\ngSki4pECU0TEIwWmiIhHCkwREY8UmCIiHikwRUQ8UmCKiHikwBQR8UiBKSLikQJTRMQjBaaIiEdB\nBaaZ/dzMdpjZJjNbamb9O2g338x2mlm+md0VTJ8iIn4Jdg3zFWCqc+4M4H3g309sYGaxwH3AZcBk\n4DozmxxkvyIiYRdUYDrnXnbONQWergay2mk2E8h3zhU45xqAp4CFwfQrIuKHUO7DvAl4sZ35w4H9\nbZ4XBuaJiESVuM4amNmrQEY7L33POfdcoM33gCbgD8EWZGaLgcWBp/VmtiXY9wyBwcBBv4sIUC3t\nUy3tUy3tm9CVb+o0MJ1zF5/qdTO7EbgcuMg559ppUgRkt3meFZjXUX9LgCWB985zzuV2VmN3i5Q6\nQLV0RLW0T7W0z8zyuvJ9wR4lnw98G7jCOVfTQbO1wDgzG2VmCcAiYFkw/YqI+CHYfZj3AinAK2a2\nwcweADCzYWa2HCBwUOhOYAWwHfiTc25rkP2KiIRdp5vkp+KcG9vB/GJgQZvny4HlXehiSRdLC7VI\nqQNUS0dUS/tUS/u6VIu1v9tRREROpKGRIiIeRVRgRspQSzO7xsy2mlmLmXV4VM/M9pjZ5sD+2y4d\ndQthLd0+/NTMBprZK2b2QeDrgA7aNQeWyQYzC+kBvs5+TjNLNLOnA6+vMbORoez/NGu50czK2yyL\nW7qpjkfMrKyjU/Cs1W8CdW4ys7O7ow6PtVxgZpVtlsnd3VRHtpm9bmbbAv8/X22nzekvF+dcxDyA\nS4G4wPRPgZ+20yYW2AWMBhKAjcDkENcxidbztN4Ack/Rbg8wuJuXSae1hGOZBPr5GXBXYPqu9n4/\ngdequ2lZdPpzAl8CHghMLwKe9rGWG4F7u/PvI9DPecDZwJYOXl9A66ASA2YDa3ys5QLg+TAsk0zg\n7MB0Cq1Dt0/8/Zz2comoNUwXIUMtnXPbnXM7Q/meXeWxlnANP10IPBaYfgy4shv6OBUvP2fbGp8F\nLjIz86mWsHDOrQQqTtFkIfC4a7Ua6G9mmT7VEhbOuRLn3PrA9FFaz9A5cYThaS+XiArME0TDUEsH\nvGxm6wIjlPwSrmWS7pwrCUwfANI7aJdkZnlmttrMQhmqXn7OD9sEPnwrgUEhrOF0agG4OrC596yZ\nZbfzejhE0v8MwLlmttHMXjSzKd3dWWC3zHRgzQkvnfZyCeq0oq4I91DLYOrwYJ5zrsjMhtJ6LuqO\nwCesH7WExKlqafvEOefMrKNTLEYElsto4DUz2+yc2xXqWqPA34AnnXP1ZnYbrWu+H/W5Jr+tp/Xv\no9rMFgB/BcZ1V2dm1g/4M/A151xVsO8X9sB0YR5q2dU6PL5HUeBrmZktpXUz7bQDMwS1hGSZdFaL\nmZWaWaZzriSw6VLWwXscXy4FZvYGrZ/uoQhMLz/n8TaFZhYHpAGHQtD3adfinGvb70O07gP2Q8j+\nPoLVNrScc8vN7H/MbLBzLuRjzM0sntaw/INz7i/tNDnt5RJRm+QWRUMtzayvmaUcn6b1gJVfFwoJ\n1zJZBtwQmL4BOGnt18wGmFliYHowMBfYFqL+vfycbWv8FPBaBx+83V7LCfvDrqB1P5oflgGfDxwV\nng1Uttm1ElZmlnF8n7KZzaQ1g0L+gRbo42Fgu3Pulx00O/3l0t1Hq07zyFY+rfsUNgQex492DgOW\nn3B0631a11q+1w11fJLW/Rn1QCmw4sQ6aD06ujHw2NoddXitJRzLJNDHIODvwAfAq8DAwPxc4KHA\n9Bxgc2C5bAZuDnENJ/2cwPdp/ZAFSAKeCfwtvQuM7sa/185q+XHgb2Mj8DowsZvqeBIoARoDfys3\nA7cDtwdeN1ov4r0r8Dvp8MyPMNRyZ5tlshqY0011zKP1GMOmNnmyINjlopE+IiIeRdQmuYhIJFNg\nioh4pMAUEfFIgSki4pECU0TEIwWmiIhHCkwREY8UmCIiHv1/dklIGL+hxyoAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "fig = plt.figure(figsize=(5,5))\n", "ax = plt.subplot(111)\n", "ax.set_xlim([-2,2])\n", "ax.set_ylim([-2,2])\n", "plt.plot(x, y);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Hurray! It worked. The orbit looks like it should, it's an almost perfect circle. There are small perturbations though, induced by the outer planet. Let's integrate a bit longer to see them. " ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUwAAAEzCAYAAABaGjpLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8VfX9x/HXJwl77xH2UECqgBEQREFxUSvOVltt1Sq/\nVm21W+tqta0d2tpWW0Fti/0p2lpQqiiCCjgYhr1k7z3DCElI8v39kau/CPcm5+aee+56Px+PPHLu\nOd97vp9cwjtnfc8x5xwiIlK9rEQXICKSKhSYIiIeKTBFRDxSYIqIeKTAFBHxSIEpIuJRzIFpZh3N\n7D0zW2Fmy83srjBtzMz+ZGZrzWyJmQ2ItV8RkaDl+LCOUuAHzrkFZtYImG9m05xzKyq1uRToGfoa\nBPw19F1EJGXEvIXpnNvhnFsQmj4MrARyT2g2GnjeVZgDNDWzdrH2LSISJF+PYZpZF6A/MPeERbnA\nlkqvt3JyqIqIJDU/dskBMLOGwH+Au51zh2JYzxhgDECDBg3O7NWrl08ViohUmD9//l7nXKto3+dL\nYJpZLSrC8gXn3MQwTbYBHSu97hCadxLn3DhgHEBeXp7Lz8/3o0QRkc+Y2aaavM+Ps+QGPAesdM79\nPkKzycDXQ2fLBwMFzrkdsfYtIhIkP7YwhwI3AkvNbFFo3k+BTgDOuaeBKcAoYC1QCNzsQ78iIoGK\nOTCdcx8AVk0bB9wRa18iIomkkT4iIh4pMEVEPFJgioh4pMAUEfFIgSki4pECU0TEIwWmiIhHCkwR\nEY8UmCIiHikwRUQ8UmCKiHikwBQR8UiBKSLikQJTRMQjBaaIiEcKTBERjxSYIiIeKTBFRDxSYIqI\neKTAFBHxSIEpIuKRAlNExCMFpoiIRwpMERGPFJgiIh75Ephm9jcz221myyIsH25mBWa2KPT1oB/9\niogEKcen9fwDeBJ4voo27zvnLvOpPxGRwPmyhemcmwXs92NdIiLJKshjmGeb2WIze9PMTguwXxER\nX/i1S16dBUBn59wRMxsFvAr0DNfQzMYAYwA6deoUUHkiItULZAvTOXfIOXckND0FqGVmLSO0Heec\ny3PO5bVq1SqI8kREPAkkMM2srZlZaHpgqN99QfQtIuIXX3bJzWwCMBxoaWZbgYeAWgDOuaeBa4Bv\nm1kpcAy4zjnn/OhbRCQovgSmc+76apY/ScVlRyIiKUsjfUREPFJgioh4pMAUEfFIgSki4pECU0TE\nIwWmiIhHCkwREY8UmCIiHikwRUQ8UmCKiHikwBQR8UiBKSLikQJTRMSjoO64LhKz8nLH/sISdh8q\nZvfhIgqOHQcgJyuL7CzIrvQ9J8uoWyubDs3q0aphHbKyLMHVSzpQYEpSOV5Wzupdh1m8pYDl2wvY\nWVDE7sPF7Cg4xt4jJTVeb6fm9encoj4dmtWjQ7OK791bNaRX20bkZGtHS7xRYErCOOfYtK+QxVsP\nsnhLAYu2HGDB5oNx6Wvz/kI27y8Mu+zcU1oxsEszBnZtwekdmlC3VnZcapDUp8CUQJWWlTNv437e\nXLqTt5bvZM/h4kSXxKzVe5i1es9nr/t3asrZ3Vowsk8b+ndsSujpKiIKTIm/0rJy5qzfzxtLd/Dy\nx5spT/KHkyzcfJCFmw/ylxnraN+kLtfmdeSK/rl0bdkg0aVJglkyP1onLy/P5efnJ7oMqaElWw/y\nwpzNvL5kO0dLyhJdTsx6tW3EVwd14rLT29O8Qe1ElyMxMLP5zrm8qN+nwBQ/OeeYsXoPY2euY876\n/YkuJ24uPq0Nd4zowekdmia6FKmBmgamdsnFFyWl5UxevJ2n3lvLhr1HE11O3E1dvoupy3dx3imt\nuGtkTwZ0apbokiQACkyJyZHiUl6cu4mnZ65n/9GaX/aTqmau3sPM1XsY2LU5P7r4VM7q0jzRJUkc\nKTClRpxzTF68nYcmL+dg4fFEl5Nw8zbs59qnZ9OvY1N+Oqo3A7sqONORAlOitmbXYe6duJT8TQcS\nXUrSWbTlIF8eO5trzuzA/V/sTdP6OjmUThSY4tmR4lL+9M4axs1an+hSkt4r87fyyvyt/Pn6/lx2\nejtdy5kmfBkTZmZ/M7PdZrYswnIzsz+Z2VozW2JmA/zoV4LhnOP1JdsZ8ug7CssofWfCQq5/Zg7b\nDh5LdCniA78G0f4DuKSK5ZcCPUNfY4C/+tSvxFlhSSnfe3kRd764kENFpYkuJyXNWb+fob9+l/Ef\nbaQs2a/alyr5EpjOuVlAVRfdjQaedxXmAE3NrJ0ffUv8rN19hBGPzeDVRdsTXUpaeGjycq5/Zg4H\nCzPvaoJ0EdQxzFxgS6XXW0PzdgTUv0Rp8uLtfHfCwkSXQbeWDTilTSNym9WjXZO65DatR/um9WhU\nN4fsLMMwzMAMsiw0jXHwWAk7C4rYdaiInQXFbD94jEVbDrJq1+GE/jzzNuxnwCPTePOuczm1baOE\n1iLRS7qTPmY2horddjp16pTgajJPcWkZv3pjJeNnbwq872E9WzK0R0vyOjejb25sdw1q26Quvdo2\nDrvMOcf2giIWbj7A/E0H+Hf+Vo4UB3e4odzBxU/M4ukbBnBJX+1opRLfhkaaWRfgdedc3zDLxgIz\nnHMTQq9XAcOdc1VuYWpoZLB2FhQx5p/5LNlaEEh/Vw3IZXDXFpzZpRndWjZI6Jnk4tIyVmw/xPxN\nB/jHRxvZeiCYkzR3jOjODy48VTc4DliyD42cDNxpZi8Bg4CC6sJSgrVlfyGXP/kBB+J8EfpNQ7ow\nul97+iXZbdPq5GTTv1Mz+ndqxq3DunG46DjTVuxi3Kz1fLIzfrvxT723jvyNB3juprNoWCfpdvjk\nBL5sYZrZBGA40BLYBTwE1AJwzj1tFf8znqTiTHohcLNzrtpNR21hBmP9niOc//jMuK3/6gEduKJ/\ne87u1iIl725eUHicqct38vDrK+K2696hWT2m3DWMxnVrxWX98nm6W5HUyOpdh7noD7N8X2+z+rV4\n5Iq+jOzdJq3uYL73SDFjZ67jmfc3+L7u3Kb1eOvuYTRSaMadAlOitmxbAZf9+QNf19m7XWPuG9Wb\noT1aJNUut9+OFJfyz9mb+M1bn/i63raN6/LOD86jgXbP40qBKVFZuPkAV/7lI9/WN7RHC+4eeUrG\n3a2nuLSMf+dv5f5Xww5yq5FWjeow80fDqV9boRkvCkzxzM/d8DM7N+NnXzqNL3Ro4sv6UlV5ecXd\nm+5+eZEv62vRoDYf/OR86tVOn8MZyaSmgZl6R+AlJvuPlvi2Gz7+loH859tDMj4sAbKyjCv657Ly\n4Uu4cXDnmNe372gJVzz1ISWl5T5UJ35RYGaQ42XlXDdudsz/Cf/n3G6sfPgSzjullU+VpY96tbN5\n5Iq+TP/+uTGva9Wuwzz4mn+7+hI7BWYGefC1ZazedaTG7+/Soj5TvjuMe0f11q5iNXq0bsSGR0fx\niytOGscRlZc+3sKkhVt9qkpipcDMEP+cvZEJ87ZU2y6Sh77Uh3d+MJw+7cMPN5STmRk3DO7M4ocu\nokOzejVez/deXszKHYd8rExqSoGZAT5cu5cHXlte4/dP//553Dy0K9kavlcjTerV4v0fj+D+L/au\n8Tou/eP7FBzTo0ASTYGZ5vYcLuZrz86t0XsHdm3Okp9dRI/WDX2uKvOYGbcO68ak24fUeB0X/2EW\n5bqfZkIpMNPcvROX1Oh9d4/syUu3DdZQPZ/179SMBQ9cWKP37jxUxBPTV/tckURDgZnG3lq2k+kr\nd0f9vme+nsfdI0/RHXTipHmD2qz71SiG9WwZ9Xv/9O5a1iT4np6ZTIGZpgqOHedb/zs/6vdN//55\nXNinTRwqksqys4x/fnMQ942K/rjm5U9+SDIPOElnCsw09fPJ0Z/kmffTC3S8MmC3nduNx689I6r3\nHDtexssf1/yKB6k5BWYa+nDtXiYu3BbVe+b+9AJaN64bp4qkKlef2YHffzm60Lxn4lI9GygBFJhp\n5lhJWdRnxef+9ALaKCwT6qoB0Yfmzf/4OE7VSCQKzDTzv3OiexbPnHsVlski2tBcuPkgCzYfiGNF\nciIFZho5VlLGL6es9Nx+9r3n07aJwjKZXDWgA49FcUzzqr98pBNAAVJgppHnZ2/03PaDn4ygXZOa\nD9eT+LnmzA7cOaKH5/bfSYLHIWcKBWaaKDpexqNverv79wu3DqJDs/pxrkhi8cOLT/V8xcLrS3Zo\nBFBAFJhpYsK8zZ7a/ejiUxnaI/oLpiV4b941zHPbp95bG8dK5FMKzDRQdLyMn/93RbXthvVsye3D\nuwdQkfihVnYW+feP9NT28WkaMhkEBWYa+He+t4uYn77hzLR+MFk6atmwDi/eOshT2/dWRT8MVqKj\nwEwDXm7d9q6eRJiyhvRoyc1Du1Tb7ua/67rMeFNgpjgvN5b97TWn062Vhjymsoe+dJqndlsPFMa5\nksymwExxY2euq3J5o7o5XHtmh4CqkXj66J7zq23zlbFzAqgkc/kSmGZ2iZmtMrO1ZnZPmOU3mdke\nM1sU+rrVj34zXXm549VF26ts8/p3ztFxyzTRvmk9rhqQW2WbbQePBVRNZoo5MM0sG3gKuBToA1xv\nZn3CNH3ZOdcv9PVsrP0KzNmwr8rld47oQecWDQKqRoLw2DXVjwLasPdoAJVkJj+2MAcCa51z651z\nJcBLwGgf1ivV+MXrVQ+D/M4F3keLSGrIyjJeGjO4yjZffUa75fHiR2DmApWva9kamneiq81siZm9\nYmYdfeg3oxUdL2NFFSd8XrxtEHVy9CjcdDS4Wwsa1418xcOOgqIAq8ksQZ30+S/QxTl3OjANGB+p\noZmNMbN8M8vfs2dPQOWlntnrI++ON6tfiyHdNZonnc2+94Iql+87UhxQJZnFj8DcBlTeYuwQmvcZ\n59w+59yn/4LPAmdGWplzbpxzLs85l9eqVSsfyktPU5ftjLhs0u1DA6xEEqFBnZwqnwl0+wsLAqwm\nc/gRmB8DPc2sq5nVBq4DJlduYGbtKr28HPB+DzIJ66UIjyho3agOXVrqRE8mGHtjxO0O5m7YH2Al\nmSPmoR/OuVIzuxOYCmQDf3POLTezh4F859xk4LtmdjlQCuwHboq130x2vKw84rK/3XRWgJVIItWv\nnUO3lg1Yr7PigfFlrJxzbgow5YR5D1aavhe414++BJZvj3yyp29ukwArkUR75dtDGPDItLDL9h8t\noXmD2gFXlN400icFRbrZxou3ebtJg6SPqgLxmffXB1hJZlBgpqAX5oa/9+XZ3VoEXIkkgxk/HB52\n/l9nVD1sVqKnwEwTv7ryCxoCmaF0ki84CswUc7S4NOz8q8+seoyxpLcHLws3Gln8psBMMZFGcWhU\nT2a7aUiXsPN1Abu/FJgpZv2eIyfNq+p6PMkMWVnhD8es2nk44ErSmwIzxfxzzqaT5g0/VSOiBB76\n0sm75W8tjzwiTKKnwEwx76/Ze9I87Y4LwI2DO5807+UII8KkZhSYKe7P1/dPdAmSJHKyT/7vXFwa\neVSYRE+BmeIuOq1NokuQJPK9kackuoS0psBMcdodl8puO7droktIawrMFNamcZ1ElyBJpn5tPUo5\nnhSYKewBXawsEigFZgrT2HEJZ0CnpokuIW0pMFNYi4baJZeTfe9CnfiJFwWmSJrJ69w80SWkLQWm\nSJqpV1tXTsSLAjOFOOc+m76iX/sEViKSmRSYKeRIpVu7XdK3bQIrEclMCswUsv3g/9/arX3Tegms\nRCQzKTBTyI6CY59NG7q7ukjQFJgppPKNFBZsPpDASkQykwIzhWRXembPvyI8OVJE4keBmUIqXy5S\n1bPJRSQ+FJgppF2TuokuQVJAaZnugRkvvgSmmV1iZqvMbK2Z3RNmeR0zezm0fK6ZdfGj30zTurEC\nU6q3cd/RRJeQtmIOTDPLBp4CLgX6ANeb2Ym30fkmcMA51wP4A/CbWPvNRA1OGMFRWBL+kbuS2Vbu\n0IPP4sWPLcyBwFrn3HrnXAnwEjD6hDajgfGh6VeAC8xM18VE6cSPbPa6fQmqRJLZu5/sTnQJacuP\nwMwFKp+y3RqaF7aNc64UKAB0b7IYjZ21PtElSBKatHBboktIW0l30sfMxphZvpnl79mzJ9HlJLV5\nG/YnugRJMsd1wieu/AjMbUDHSq87hOaFbWNmOUATIOz+pHNunHMuzzmX16qVnrddnco35BDJ36gB\nDfHkR2B+DPQ0s65mVhu4Dph8QpvJwDdC09cA7zr9T/fFB2tPfk65ZK4J8zZ/7nWvto0SVEl6ijkw\nQ8ck7wSmAiuBfznnlpvZw2Z2eajZc0ALM1sLfB846dIjqZlfv/lJokuQJOGcY/Li7Z+bN6JX6wRV\nk558ecScc24KMOWEeQ9Wmi4CrvWjL/k8jfiRT63bc/L1l2d00PN9/JR0J30kersPF1XfSNLetBW7\nTprXr6MC008KzBRzTo+WJ8178NXlCahEks1v3jr58IyeXe8vBWaKuXVY15PmvbV8ZwIqkWSy7eCx\nsPM1PsRfCswUM7Br+CcCLt5yMOBKJJk8+e6aRJeQERSYKaZ+7fDn6a4bNyfgSiRZFJaUMmHeyfdH\n7d2ucQKqSW8KzDRx7HiZbuuVoSINhbxjRPeAK0l/CswUVCs7/HGpX05ZGXAlkmjOOe6btCzssuGn\n6hpMvykwU9Bj154Rdv7fP9yooZIZpqqRXg3r+HKZtVSiwExBVT2T/PnZmwKsRBLtl29oryJICswU\nVCcnO+KyhyYv11Zmhli+vYBPdoa/WfD3Lzwl4GoygwIzRbWv4vk+f5mxLsBKJBGcc9zw7NyIy8ec\n2y3AajKHAjNFPXr16RGX/W7qKsrLtZWZzqYs3cmBwuMRl9etFXkvRGpOgZmizu158hDJyu6duDSg\nSiRoRcfLuOPFBRGXX3xamwCrySwKzBRV3ZC3l/O3cLgo8haIpK6xM6t+NMm9l/YOqJLMo8BMYX+8\nrl+Vywf96p2AKpGgbD94jD9MX11lmy4tGwRUTeZRYKawy05vX+XywpIy3ll58i2/JHVVd6hlWDWH\naiQ2CswUlp1l1T6C4Jvj8zlarOeXp4MpS3cwc3XVDwZ84itV73VIbBSYKe7pG86sts15v5sR/0Ik\nrjbvK+T2FyKf6PlUi4a6/2U8KTBTnJfjVXuPFPPK/K0BVCPxUFJazjVPf1Rtuwm3DQ6gmsymwEwD\nv7yyb7Vtfvjvxew/WhJANeK3R99cye7DxdW2O7t7iwCqyWwKzDTwlbyO1TcCBjwyjZJS3QIulUxb\nsYu/f7ix2nY/vEhDIYOgwEwDOdlZ3Dy0i6e2l/35fY01TxFbDxRy2/P5ntreOkxDIYOgwEwTP7mk\nl6d2q3cd4WeT9dC0ZFdQeJzzH5vpqe2X8zpoKGRAFJhpom6tbO4e2dNT2/GzNzFh3uY4VyQ1dbS4\nlMuefJ8Sj3fQf+SK6o9hiz8UmGnkO+d7C0youAB69rp9caxGaqK4tIwbn5vLlv3hnwJ5ogcu61Pl\n7f7EXzEFppk1N7NpZrYm9L1ZhHZlZrYo9DU5lj4lsuws4xdRbG1c/8wcVkW4n6IEr7SsnDteWMCC\nzd6fAHqLx2PX4o9YtzDvAd5xzvUE3gm9DueYc65f6OvyGPuUKnxtUKeo2l/8xCyWbi2IUzXiVXm5\n48f/WcL0lbs9v+elMYP13PGAxRqYo4HxoenxwBUxrk9iZGa88q2zo3rPl578gHkb9sepIqlOebnj\ngdeWMXFB+Kc/htOsfi0Gd9N1l0GLNTDbOOd2hKZ3ApFuxFfXzPLNbI6ZKVTjLK9Lc87qEvboSERf\nHjubaSt0o46gHSsp45vjP+aFudGdhJt0+9A4VSRVqTYwzWy6mS0L8zW6cjtXcXFfpAv8Ojvn8oCv\nAk+YWcQHJpvZmFC45u/ZU/WNBiSyv988MOr33PZ8vs6eB2j3oSIu+eMs3lsV3e/5nSN66BZuCVJt\nYDrnRjrn+ob5eg3YZWbtAELfwx6Acc5tC31fD8wA+lfR3zjnXJ5zLq9Vq1Y1+JEEKh6x+tKY6McW\n3ztxKb+b+kkcKpLKVmw/xKBH32HTvsKo3/sDjepJmFh3yScD3whNfwN47cQGZtbMzOqEplsCQ4EV\nMfYrHgzu1oIRp0b/R+ep99ZxyROzdFu4OHn3k12M+tP71GTAVf79I3WiJ4FiDcxfAxea2RpgZOg1\nZpZnZs+G2vQG8s1sMfAe8GvnnAIzIOO+nlej932y8zCnPTSVFdsP+VxR5iovdzz7/npu+Ye34Y4n\neu4bebTU7dsSypJ5XHFeXp7Lz6/ZL5f8v7W7jzDy996G2YXzwGV9uGVoF23ZxGDL/kLufHEBi2t4\nCdcXT2/HU18d4HNVmcvM5ofOq0RFI30yQI/WDWMaPvfI6ys4//GZFBzTQ9Wi5ZzjhbmbGPbb92oc\nlgCPX3uGj1VJTSkwM8SNgztz3VnebgMXzoa9Rznj528zdflO3e3Io20Hj/G1Z+dy36RlMa3n/R+P\n0M01koQCM4P86sovcFr7xjGt43/+OZ8Rj81g7W4NqYzEOcdL8zYz9Nfv8lGM4/Xf+O45dGxe36fK\nJFYKzAySlWVMvH1IzOvZuK+Qkb+fxb0Tl1BQqN30Tznn+GDNXkY8NoN7qnm6oxcv3jqI09o38aEy\n8YsCM8PUyclmyc8u8mVdE+Zt4YyH32b8RxspK8/s3fT5m/ZzzdOzueG5uWyswbWVJ3ryq/0Z0kOP\nzE02OkueoXYWFDH40Xd8Xedvrz6dy/u1z6jjbcu2FfDY26uYEeVonar87Et9uGloV9/WJyer6Vly\nBWYGW7PrMBf+YZbv6/3uBT25YXAnWjeq6/u6k4FzjgWbDzJ25jre9nn8/beHd/d893ypOQWm1MjW\nA4Wc85v34rLuC/u04a4LetI3Nz2Ow+05XMzEBVt5euY6DsTh2O2PLj6VO0b08H29cjIFptTYgaMl\n9H9kWtzWXyvbeOCyPozs3Yb2TevFrZ94OF5WzoxVe/hX/pa43s3pj9f1Y3S/3LitXz5PgSkxOVZS\nxvmPz2BHQVFc+8ltWo+vnNWRC/u0oVfbRkk5euhQ0XHmrNvHR+v28Y+PNsa9v5fHDGaQ7m0ZKAWm\nxKy0rJw7X1zIW8t3BtbnjYM7M6JXK/q2b0Lrxok55ll0vIz5mw7w4dq9zFy9h+UBjp+f/v1z6dG6\nUWD9SQUFpvjCOccT09fwx3fWJKT/YT1b0je3Cae1b8xp7ZvQuXl9srL82QotK3dsP3iMjfuOsnHv\nUTbsLWT+5gMs3uL9GTp+aVq/Fm/ffW7C/khkOgWm+OrDtXv52rNzE13GZ/rmNia3aT2aN6hNs/q1\nP/vetH4tysodRaXlFB8v++x7cej70ZIyNu8vZPWuwzW692Q8XDUgl19c0Zf6tXMSXUrGUmCK7w4c\nLeGW8R+zMIqnGErV/n7TWYzo1TrRZWQ83a1IfNesQW0mfntIVI/ulfAGdmnO/PtHKixTnPYJpEpm\nxg2DOzOkewuuf2YOuw4VJ7qklPPba07n2jM7JOUVARIdbWGKJ91aNeSDn5zPPZdqFIpXAzo1ZdaP\nRvDlvI4KyzShLUzxrFZ2Ft86rztfzuvI42+vivrRsJnCDJ6/ZSDDeuohfulGW5gSteYNavPLK7/A\nez8czoBOTRNdTlJ54iv9WPfLUQrLNKUtTKmxri0bMPH2oeRv3M89E5eydveRRJeUMPd/sTc3nt2Z\nOjmZc6emTKTAlJjldWnOtO+dy8zVe3jugw28v2ZvoksKzA8vOoUbz+5Ck3q1El2KBECBKb4wM4af\n2prhp7Zmy/5CXpy3mb/OWJfosuJiYNfmfPu87px7SiuyfRqFJKlBF65L3BSXlvHWsp2MnbmeFTtS\n//nm3z2/B18Z2IncFLvjkpyspheuawtT4qZOTjaj++Uyul8ua3cfZvrK3UxZuoMlMTxuNmjXndWR\n83u1ZkSv1tTK1jnSTKfAlED0aN2IHq0b8a3zunOo6DgfrNnL28t38uqi7Yku7XMGdm3ORX3acN4p\nrejRuqGun5TPiWmX3MyuBX4G9AYGOufC7j+b2SXAH4Fs4Fnn3K+9rF+75OmvvNyxYsch5qzfx5pd\nR1i45QCrdwV3tv28U1rRu11jBnVtzqBuzXVDjAyRqF3yZcBVwNhIDcwsG3gKuBDYCnxsZpOdcyti\n7FvSQFaW0Te3yeceY1Fe7th28Birdx1m1a7DrN55mO0FRew9UsyGvUeJ5m987ZwsurdqSIdm9Til\nTUN6tm5Ej9YN6d6qIfVq6xIgiU5MgemcWwlUt9syEFjrnFsfavsSMBpQYEpYWVlGx+b16di8Phf0\nbhO2jXOOwpIyDheVcrjoOA5oWCeHBnVyaFA7mxwdb5Q4CGL/IxfYUun1VmBQAP1KGjOzinCsk0Pb\nJroJrwSj2sA0s+lA2zCL7nPOveZ3QWY2BhgD0KlTJ79XLyJSY9UGpnNuZIx9bAM6VnrdITQvUn/j\ngHFQcdInxr5FRHwTxIGej4GeZtbVzGoD1wGTA+hXRMRXMQWmmV1pZluBs4E3zGxqaH57M5sC4Jwr\nBe4EpgIrgX8555bHVraISPBiPUs+CZgUZv52YFSl11OAKbH0JSKSaLr2QkTEIwWmiIhHCkwREY8U\nmCIiHikwRUQ8UmCKiHikwBQR8UiBKSLikQJTRMQjBaaIiEcKTBERjxSYIiIeKTBFRDxSYIqIeKTA\nFBHxSIEpIuKRAlNExCMFpoiIRwpMERGPFJgiIh4pMEVEPFJgioh4pMAUEfFIgSki4pECU0TEo5gC\n08yuNbPlZlZuZnlVtNtoZkvNbJGZ5cfSp4hIouTE+P5lwFXAWA9tRzjn9sbYn4hIwsQUmM65lQBm\n5k81IiJJLKhjmA5428zmm9mYgPoUEfFVtVuYZjYdaBtm0X3Oudc89nOOc26bmbUGppnZJ865WRH6\nGwOMAejUqZPH1YuIxF+1gemcGxlrJ865baHvu81sEjAQCBuYzrlxwDiAvLw8F2vfIiJ+ifsuuZk1\nMLNGn04DF1FxskhEJKXEelnRlWa2FTgbeMPMpobmtzezKaFmbYAPzGwxMA94wzn3Viz9iogkQqxn\nyScBk8LQ+9SxAAAGC0lEQVTM3w6MCk2vB86IpR8RkWSgkT4iIh4pMEVEPFJgioh4pMAUEfFIgSki\n4pECU0TEIwWmiIhHCkwREY8UmCIiHikwRUQ8UmCKiHikwBQR8UiBKSLikQJTRMQjBaaIiEcKTBER\njxSYIiIeKTBFRDxSYIqIeKTAFBHxSIEpIuKRAlNExCMFpoiIRwpMERGPFJgiIh7FFJhm9jsz+8TM\nlpjZJDNrGqHdJWa2yszWmtk9sfQpIpIosW5hTgP6OudOB1YD957YwMyygaeAS4E+wPVm1ifGfkVE\nAhdTYDrn3nbOlYZezgE6hGk2EFjrnFvvnCsBXgJGx9KviEgi+HkM8xbgzTDzc4EtlV5vDc0TEUkp\nOdU1MLPpQNswi+5zzr0WanMfUAq8EGtBZjYGGBN6WWxmy2Jdpw9aAnsTXUSIaglPtYSnWsI7tSZv\nqjYwnXMjq1puZjcBlwEXOOdcmCbbgI6VXncIzYvU3zhgXGjd+c65vOpqjLdkqQNUSySqJTzVEp6Z\n5dfkfbGeJb8E+DFwuXOuMEKzj4GeZtbVzGoD1wGTY+lXRCQRYj2G+STQCJhmZovM7GkAM2tvZlMA\nQieF7gSmAiuBfznnlsfYr4hI4KrdJa+Kc65HhPnbgVGVXk8BptSgi3E1LM1vyVIHqJZIVEt4qiW8\nGtVi4Q87iojIiTQ0UkTEo6QKzGQZamlm15rZcjMrN7OIZ/XMbKOZLQ0dv63RWTcfa4n78FMza25m\n08xsTeh7swjtykKfySIz8/UEX3U/p5nVMbOXQ8vnmlkXP/uPspabzGxPpc/i1jjV8Tcz2x3pEjyr\n8KdQnUvMbEA86vBYy3AzK6j0mTwYpzo6mtl7ZrYi9P/nrjBtov9cnHNJ8wVcBOSEpn8D/CZMm2xg\nHdANqA0sBvr4XEdvKq7TmgHkVdFuI9Ayzp9JtbUE8ZmE+vktcE9o+p5w/z6hZUfi9FlU+3MCtwNP\nh6avA15OYC03AU/G8/cj1M+5wABgWYTlo6gYVGLAYGBuAmsZDrwewGfSDhgQmm5ExdDtE/99ov5c\nkmoL0yXJUEvn3Ern3Co/11lTHmsJavjpaGB8aHo8cEUc+qiKl5+zco2vABeYmSWolkA452YB+6to\nMhp43lWYAzQ1s3YJqiUQzrkdzrkFoenDVFyhc+IIw6g/l6QKzBOkwlBLB7xtZvNDI5QSJajPpI1z\nbkdoeifQJkK7umaWb2ZzzMzPUPXyc37WJvTHtwBo4WMN0dQCcHVod+8VM+sYZnkQkun/DMDZZrbY\nzN40s9Pi3VnosEx/YO4Ji6L+XGK6rKgmgh5qGUsdHpzjnNtmZq2puBb1k9Bf2ETU4ouqaqn8wjnn\nzCzSJRadQ59LN+BdM1vqnFvnd60p4L/ABOdcsZn9DxVbvucnuKZEW0DF78cRMxsFvAr0jFdnZtYQ\n+A9wt3PuUKzrCzwwXcBDLWtah8d1bAt9321mk6jYTYs6MH2oxZfPpLpazGyXmbVzzu0I7brsjrCO\nTz+X9WY2g4q/7n4Eppef89M2W80sB2gC7POh76hrcc5V7vdZKo4BJ4Jvvx+xqhxazrkpZvYXM2vp\nnPN9jLmZ1aIiLF9wzk0M0yTqzyWpdskthYZamlkDM2v06TQVJ6wSdaOQoD6TycA3QtPfAE7a+jWz\nZmZWJzTdEhgKrPCpfy8/Z+UarwHejfCHN+61nHA87HIqjqMlwmTg66GzwoOBgkqHVgJlZm0/PaZs\nZgOpyCDf/6CF+ngOWOmc+32EZtF/LvE+WxXlma21VBxTWBT6+vRsZ3tgyglnt1ZTsdVyXxzquJKK\n4xnFwC5g6ol1UHF2dHHoa3k86vBaSxCfSaiPFsA7wBpgOtA8ND8PeDY0PQRYGvpclgLf9LmGk35O\n4GEq/sgC1AX+Hfpdmgd0i+Pva3W1PBr63VgMvAf0ilMdE4AdwPHQ78o3gW8B3wotNypu4r0u9G8S\n8cqPAGq5s9JnMgcYEqc6zqHiHMOSSnkyKtbPRSN9REQ8SqpdchGRZKbAFBHxSIEpIuKRAlNExCMF\npoiIRwpMERGPFJgiIh4pMEVEPPo/KtcyO3+TPlIAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Noutputs = 1000\n", "times = np.linspace(2.*torb, 20.*torb, Noutputs)\n", "x = np.zeros(Noutputs)\n", "y = np.zeros(Noutputs)\n", "for i,time in enumerate(times):\n", " sim.integrate(time, exact_finish_time=0)\n", " x[i] = particles[1].x\n", " y[i] = particles[1].y\n", " \n", "fig = plt.figure(figsize=(5,5))\n", "ax = plt.subplot(111)\n", "ax.set_xlim([-2,2])\n", "ax.set_ylim([-2,2])\n", "plt.plot(x, y);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Oops! This doesn't look like what we expected to see (small perturbations to an almost circular orbit). What you see here is the barycenter slowly drifting. Some integration packages require that the simulation be carried out in a particular frame, but WHFast provides extra flexibility by working in any inertial frame. If you recall how we added the particles, the Sun was at the origin and at rest, and then we added the planets. This means that the center of mass, or barycenter, will have a small velocity, which results in the observed drift. There are multiple ways we can get the plot we want to.\n", "1. We can calculate only relative positions.\n", "2. We can add the particles in the barycentric frame.\n", "3. We can let REBOUND transform the particle coordinates to the barycentric frame for us.\n", "\n", "Let's use the third option (next time you run a simulation, you probably want to do that at the beginning)." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "sim.move_to_com()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So let's try this again. Let's integrate for a bit longer this time." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUwAAAEzCAYAAABaGjpLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X1U3Od5J/zvNViyXiJkwDaoaV0HPzkbDejsbg0D7ub0\ndLe7W1mABCj7bLrnWIJBQYqb2q1lCQYY7TkCLCGfuE9eHFszw2C3f/RVQDSgl6dttqcrmGFQcpKI\nQU+2jtPn1El2s232xGnXlQRz7R/M/csPjKQRMzBv3885czTAMHNr9NM198t1X7eoKoiI6P4cmW4A\nEVGuYMAkIkoSAyYRUZIYMImIksSASUSUJAZMIqIkpSVgikhQRH4kInN3+fmvishPROSbidupdLwu\nEdFGeihNz/MWgC8D+L17POa/qmpjml6PiGjDpaWHqap/BeDH6XguIqJstZFzmM+IyLdE5LKIVG3g\n6xIRpUW6huT38w0Av6iq/yAi+wCMA/j4ag8UkU4AnQCwffv2pz/xiU9sUBOJqFB8/etf/ztVfexB\nf0/StZdcRJ4EMKGq1Uk89m8A1Kjq393rcTU1NXr9+vW0tI+IyBCRr6tqzYP+3oYMyUWkQkQkcd+V\neN2/34jXJiJKl7QMyUXkDwD8KoBHReQ9AP8ZwCYAUNU3AXwKwGdFZAHABwA+rSyTREQ5Ji0BU1V/\n4z4//zKW0o6IiHIWd/oQESWJAZOIKEkMmERESWLAJCJKEgMmEVGSGDCJiJLEgElElCQGTCKiJDFg\nEhEliQGTiChJDJhEREliwCQiShIDJhFRkhgwiYiSxIBJRJQkBkwioiQxYBIRJYkBk4goSQyYRERJ\nYsAkIkoSAyYRUZIYMImIksSASUSUJAZMIqIkMWASESWJAZOIKEkMmERESWLAJCJKEgMmEVGSGDCJ\niJLEgElElCQGTCKiJDFgEhEliQGTiChJDJhEREliwCQiSlJaAqaIBEXkRyIyd5efi4h8UUTeEZFv\ni8gvpeN1KX+oKmZmZqCq1teRSASRSATxeBwzMzNYXFzE8PAwpqenEQ6Hl/0sHo8jHA5jeHgYi4uL\ny56LKG1UNeUbgF8B8EsA5u7y830ALgMQAPUAZpJ53qefflopP8XjcY1EIhqPx1VVdXp6WouLi/Xa\ntWsaCAR0enpaS0tLdefOnRoIBHTXrl3q8XgUgG7evFmLi4t1586d6vf7taysTAOBgBYXFysA3bdv\nnz7++OMaCAR0cXFRp6enNRAI6O3bt7Wnp0cXFhYy/LenTANwXdcQ69LSw1TVvwLw43s85ACA30u0\nNQLgERHZlY7Xpuynqlbvb2FhAYFAAIFAAA0NDVZPMBQK4ac//SkGBwfR2dmJUCiEoaEhFBUVwel0\nYnR0FB/72Mewfft2OBwOPP/88ygqKrJ6kU6nE5cvX8azzz6LK1euoKamBn19ffB6vdi7dy8+85nP\noLm5Ga+88grcbjf8fj8CgcCHeqO6oqdLtMxaouxqNwBP4u49zAkAn7R9/RcAau73nOxh5h7Tc1xc\nXLT+NL0/EdHDhw+riOiWLVu0uLhYw+GwhsNhLSkp0cOHD+ujjz5q9RBNz3BxcVHD4bCWlpaqx+PR\n0tJSnZ6e1kgkogsLC8seY57n8ccf156eHi0vL1efz6cej0dv3bqlhw8f1i1btigABaAej0fLysp0\nampKA4GAXrt2TYuLi9Xn82k4HLZ6wJRfsMYeZtYFTACdAK4DuP7EE0+s09tF6RSPx63ANzU1pTt3\n7lSfz2cNlSsqKtTn86nf79dr166px+PRkpIS9fv9Go/HreG4CVrl5eUaCAQ0HA7rrl27NBKJaDgc\n1rKyMitQmkAWiUTu+hgTtFc+T0lJiT733HO6Y8cOPX/+vJaVlanH41ER0UOHDikAffjhh7W0tFTD\n4fCy16P8kO0B8zyA37B9/R0Au+73nOxh5oZwOKw7d+60eoAAtLu7e1nPzfQATUAzAdYE27KyMutr\nE6Dudt9urY+x94AjkYhOTU1pcXGxvvnmm/rwww9b86MmgHs8Hp2amlr296Hcle0BswHLF32iyTwn\nA2b2MIFmYWFB/X6/Tk9PW8HGLNB4PB5rSDs1NWX16CoqKqyFHBMY7T3DuwW6TPz9zN/F9HArKiqs\naYTt27frww8/rAD02Wef1ampKfY8c1RGAyaAPwDwQwB3ALwHoAPAMQDHEj8XAK8D+C6AG8nMXyoD\nZlYxvcBDhw5ZwcOsXofD4WVDadODND06Myy3D2+zIUiuZrXeqPmQ8Pl8+sgjj+gzzzyjAHTLli3q\n8/l0enp62d+Zsl/Ge5jrcWPAzA7xeFz9fr8WFxdrSUmJHjp0SEtLS5cNT+2B0QTSSCRi/X42BscH\nZf4e165d04cffli3b99upTft2LHDWizikD37rTVgytLvZqeamhq9fv16pptRUFSX0moAwOVyYXZ2\nFqqK5uZmdHR0oKGhASICEQEAtLS0YGBgAG63GwAQjUZRW1uL2dlZuFwu63H5xLxH9v87sVgML7/8\nMj744APcvn0bHo8H+/fvR11dXV6+B7lORL6uqjUP/ItribIbdWMPc+OY9B+z+FFWVqY9PT3WUNoM\nuXt6epYt0JhepelNFirzXnzkIx/RrVu36kc+8hH2OLMY1tjDfCj9sZtyjarC6/Xi7Nmz6Orqgqqi\ns7MTgUAAR44cgcvlQl1dHQDgxIkTWFxcBACICNxuN6qrq+FyuTL5V8g4815UVVVBVTE3N4fjx4/j\nhRdewD/90z8tDedEsHv3bjgcDvY8c9VaouxG3djDXB8r5xQjkYjVe5yamtLS0lIrpcbee7TnW+b6\nfOR6Mz1O01P3+XwKQLdu3Wol7FPmgIs+lKzVErxNIFxYWLB2yJjvMUCujf2DaXp6Wnfu3Knd3d26\nc+dO9Xg8urCwkBeLYblorQGT5d0KSDwet6r53LlzB7FYDC0tLZidnYWIoKWlBadOncLw8DAGBwet\nYWNrayui0Wimm59zRMR6D+vr63H16lXs378ft2/fxpkzZ+D1etHS0sL3NocwYBaQkZERdHZ24tKl\nS3jooaXp69HRUbhcLrhcLgwMDFjB0u12Q0TgcrkwNjZW8HOUqTLBs76+Hi+++CJEBE8++SROnz6N\neDxulalb6vxQtmJaUYFQXaovGYvFsHv3bkxOTiIYDC4LjqqKaDSat+lA2SIejyMYDAIAurq6AADH\njh2zFtn6+/vhcLAvs56YVkTLmDQhk85itiKahYjS0lJrrjIQCHAebYNFIhEryd/v92tpaak+99xz\nVgUl/nusLzCtiABYvcQbN26gs7MTqoo9e/agpqYG/f39aGtrg9PpRCwWQ3t7OyorK9HX14fq6mor\ndYjWn8vlwvj4OFwuF2ZmZiAi+OhHPwoA+NKXvoTKykrr34S9/Syylii7UTf2MB+cWQH3+XxaXFys\n3d3dVk/G/Gkvd5Yv2xZzmX3PeiAQ0PPnz+u2bdu0tLS04DcErBdwlZzsqqur8dprr2FkZAQDAwNo\nb2/HwMAA+vr6AMBayLGv5FJmmH+DoqIidHR0oKioCB988AGOHj1qnVUUiUS4IJQFOCTPE2pbsJmY\nmICqwul0YmxszNr7bd+VwwCZvdrb263g2NDQgNu3b2Pz5s347Gc/ywWhDOM7n8NUf1YEIhqNLsup\nbGpqQlNTE+bn59HY2GjNk7E3mf0cDgf27NkDr9eLY8eOoaioCI2NjVbuJnuamcOAmaNUFcFgEAcO\nHEAwGEQ8Hl+WUxkKhXD27Fns3r07002lNTCLQk1NTdi0aRM++clPYvPmzXj99dcRDAYZNDOEeZg5\namZmBs3NzXC73fjyl78MALhy5QocDgdcLpfV4xwdHbUS0NmzzD1m9BCPx9HU1ITOzk74fD6EQiHr\n35r/rg9urXmY7GHmILNiNzY2hsbGRut7sVgMzc3NCAaDqK2txdjYGOrq6jgMz2H2HUITExOorKy0\n/q3NMcW0cbjok4NM79Gsfl+5csXKqxQR5lXmIbNw5/V6ce7cOezevRu3b9/Gt7/9bczNzcHtdnMx\naCOsJRdpo27Mw1yy8pwZU9DXXtzXlGFjXmX+sp94GQgEdMeOHbp9+3buDloDsLxb/lq5jc5+oJg5\n2bCnp4dVvQuEuR78fr92d3fr1q1btaSkhFtcH8BaAyaH5FnO/EP19/fj5MmTAJYKNZjJ/pmZGSsh\n3ZwhQ/nNrKCrLlXK/+IXv4jvfve7OHHiBKqqqlBfX5/pJuYtrpJnuZmZGbS0tODChQuYn59HPB7H\nqVOnMD4+DgBcCS9gmlhBV1U0Njbizp07+PznP8896ElY6yo5A2YWU9Vlq6BmocfpdFoBMp9PZ6Tk\nqC6V7guFQnjzzTfhcDgwNDRkle2jD2NaUR6KRqNobW21UkvM0Ht+fh6tra2YnZ1lT4IgInA4HBgZ\nGcGrr76Ks2fP4sSJE+jt7UU8Hs908/IKA2aWMnOXFy5csO47nU4MDAygra2NVdBpGTOvaVLLbt26\nhbNnz3IrZbqtZaVoo26FuEpuUkamp6etgr/2wr9lZWUs+UWrMidVlpeXq8/n00OHDuljjz3G1fNV\ngKvk+cGcu3P+/Hmr4G91dTVqamqsnkJtbW2GW0nZKBqNoq+vDx0dHaiursapU6fwmc98hhsZ0ogB\nM0toYsWzra0NAOB0OtHa2goAcLvdiEaj6O7uBgDs2bOHFz99iDnIrq+vD01NTRgbG7NW0E0PifPd\nqeEcZpawl2erqqoCACv3cmZmxqpzOTExwblLWpWpeTo+Pm4tBh48eBDz8/NoaGhglaM0YMDMAubT\nf3R0FADQ2NiIpqYma/8wAOts6/r6evYS6K7sNU9dLhdGR0ehqrh9+zZeeuklDA8PM2imgAEzC6xM\nHzK1LNva2jA5OcnhN62JvWDHCy+8ABHBSy+9xJ5mCjiHmUFm3tKc6FhTU4OZmRnEYjF4vV7OVVLK\nTLpRTc1SjvYXv/hFHD9+nFso14g9zAyamZlBQ0MD3nrrLXi9Xrz11ltobGxEV1cX3G43V8MpbWZn\nZxEMBvG5z30OqoobN26wl7kGDJhZYPfu3VYK0cTEBIaGhhAMBjE7O5vpplGOM4uJqoqBgQFrbry7\nuxvRaDTTzcs5HJJngBmKu1wuTE5OQlWt42/b29sBgDt5KC1cLpeVXuT1enHhwgV8/vOfx+7du5lq\ntAbsYWaA+dQ3n/C1tbVwu93o7e3FyMiItQDEC5lSZRYS6+rqrCOXvV4vbt68iZaWFi4APSD2MDOg\ntrYW/f39UFW0traiv78fwWAQHR0d1s4e9i4pnUzgVF06C6q2thaqipMnT0JV0dHRwQ/oJKSlhyki\ne0XkOyLyjoh0r/LzNhH5nyLyzcTtSDpeN1fNzs7C6/VCRDA2Noa2tja43W4MDw/j+vXrrEBE68Ze\nFtDpdOLWrVt46aWXeJhaklIOmCJSBOB1AM8CcAL4DRFxrvLQP1LVf5G4BVJ93VykibqF8XgcFy5c\nALA0x3T9+nUEg0EMDg6yZ0nrLhqNorm5GZOTkygqKsLi4iKH5UlKRw/TBeAdVX1XVW8D+EMAB9Lw\nvHlnZmYGe/fuRWNjI+bn5605pNraWoyPj7PgK20IM2ceCATwm7/5mygqKkIsFmPQTEI65jA/CuBv\nbV+/B2C1bOuDIvIrAP4bgN9R1b9d5TF5r6ioCEePHrWKbPT29gIAgyVtGJOT+corr8DpdOIrX/kK\nurq6rL3ovA7vYS014ew3AJ8CELB9/RyAL694TBmAhxP3jwL42j2erxPAdQDXn3jiiTVWu8suK49H\nLS8v10AgoAsLC9rT06Pl5eWscUkbxn49hsNhnZ6etk4jLZTrEJk6ZhfAMwCu2r72APDc4/FFAH6S\nzHPnSwHhSCSiu3bt0nA4rOFw2Lo4zXniLPBKmWC/Lk3QLJSjmtcaMNMxJJ8F8HER+RiA7wP4NID/\nZH+AiOxS1R8mvtwP4GYaXjcnmDfaVCIyB5mZr8fGxrgqThlhT2pvamqCquJ73/se+vv74XAwRXs1\nKb8rqroA4HMArmIpEP6xqsZE5LSI7E887AURiYnItwC8AKAt1dfNFfc6yOzgwYNMUKeMsSe1T0xM\n4OjRozh79izcbjcPT7sLHrO7jlSXjsnVxPYzU3lo5fcYMCmTNLFV9+mnn8aBAwdw6dIl9PT0YGBg\nIG+vTR6zm4VM79KeQgTAqoTN3iVlA7NV9+2338Y3vvEN7Nu3D1/5yleYzL4K9jDXkelhxuNxTExM\nIBgMYnx83Po5e5eUDewjoVgshq6uLiwsLODq1at5WzOTPcwsZHqQ+/fvx/nz5zE4OGgdSmV+TpQN\nYrEYWltbUVVVhXPnzuHy5csAwGT2FRgw14n51K6trbUOL6uqqoKq4s6dO7wQKWuY43nNnKWpZtTa\n2sqamSswYK4Ts193ZGQEdXV1cDgcaG1txc2bN7Fp0yb2LilrmPSiqqoq6357e7t1gBo/3H+GAXMd\nmIvMpBCZeczTp0/zYDPKOmbqyN6jNB/ojY2NXPyxYT3MdWBWx0dHRzE+Pm4lBgPgwWaUlexJ7C0t\nLRgdHUUsFkM8Hsfc3BwXKBMYMNPM9C7t5dtUFUNDQ3A6nSzfRlnJXmB4dHQUc3NzOHnyJI4ePWqd\nMvnMM89kupkZxyF5mq2WexmNRuH1euFwOPgpTVnNDM+7u7uxuLiIH/zgB3j//fcxPz+f6aZlBeZh\nptnK3Mvh4WEMDg4CWCrhxj26lO00Ueg6FApheHgYR44cQWNjI+rr6/PmA3+teZgckqeZPfdSVXHs\n2DGr1iDnLynbmW2SIoKRkRGrZmZTUxMmJibyNpE9WQyYaWIuNJfLBZfLhYmJCSsdo6GhAfPz86it\nrc1wK4nuzWyTNAuWtbW1GBkZYWpRAseHaWI/OldEUF9fD4fDgYMHD+LmzZs4deoUZmdnM91Monsy\nq+WmitHs7Cz6+vpw7NgxLliCATNtXC6Xlegbj8etw85GR0fR3t6OsbExXnCU9cypktFoFPF4HKqK\n/v5+DA8Ps6cJBsy0isViaGhowMjICBobG9HY2IhYLMYybpRTzGhpZGTE2l/e0dGBnp4eBIPBgg6a\nXCVPk5mZGTQ0NODOnTu4cuUKRARzc3Po6urC5ORkwU+WU+4w8/G1tbWYnZ21ktl//dd/HVeuXMFX\nv/rVnF+8ZLWiDDMLPVevXrV6k1VVVexVUs6xj4jMRoz+/n5cuXIFHR0dBb14yYCZRmb+Z2RkBM3N\nzRAR7hunnBWNRtHY2IimpiZrWP7GG28UdAUjphWliZn36e/vt0plZfN0B9H9uFwuhEIhxGIxAMD5\n8+exuLiY4VZlFgNmiuzzPWNjY6itrUVVVRXm5ubQ2NjIXiblJHsCu9frxcDAAC5evIj5+fmCzvbg\nkDxFpmc5Oztr1b00F9mxY8cQCoUK+gKj3GSuawBWmULmE7OHmbKV+Zezs7Oora21LrL9+/dz4Ydy\njklgd7lcqKurQ3V1NWpqlhaVuehDa2YvvmoWe0ZGRuB0OpmsTjlrtWOho9EoTp48ib6+voI9t5w9\nzBStrH/Z39+PkydPWnOX7F1SLjMr5QAwNDSEW7du4cyZM6isrMSRI0cy3LqNxx5miuz1Lw8ePIjq\n6mpMTEzg7NmzBT10ofxg8otDoRCcTid+93d/F8XFxaiqqsp00zKCATNFKw+NApaGM4U+OU75wQzN\n5+fn0draiurqarz22msFm/XBgJkik6xutpA1NjZCVTl/SXnDfgxvLBZDb29vwRbiYMBMgamuPjMz\ng5aWFquMfywWK8iLifKTyfpwOp3wer3o6OhAX19fQe744aJPCuzFVk3SOgD8zu/8DjZt2oRLly4V\n7NCFcp9ZGVdVeL1eXLhwAQMDAzh8+DAqKysLco6ePcwUmBxMAFZlF6fTic2bN+Po0aMFeUFR/rAn\nr4+NjVkbMt5++214vd6CnKNnDzMFJgfTvoe8v78fnZ2dCAaDOHDgAHuYlLPsyeumctHY2JiVwG4K\nDBdS6hzrYaYoHo9jZGQEu3fvxvz8PLq6ugAA586dg9vtLqiLifKXma8395999lkUFRXl7LQTT43M\nkNnZWStITkxMYHJyEgBYYZ3ygn0e0ySwnz17FqqKs2fPFlwmCANmilaeEGlSjIjygZnHvHDhAoaG\nhuB0OgEAmzZtgsNReEsghfc3TjP7CZFmP7k5PZIo15mFzVgshr6+PjgcDtTX12NoaKggU4vYw0yB\nfV6ntrYW/f39+MQnPoHTp09zhZzyglnYNDUxa2trMTMzU7DFZdISMEVkL4AvACgCEFDVsyt+/jCA\n3wPwNIC/B/AfVfVv0vHambSyMEFXVxc++OADbNmyBXv27MnJyXCilVwuF8bHx63jd801X4jFZVIO\nmCJSBOB1AP8OwHsAZkXkoqrO2x7WAeB/qer/JSKfBjAE4D+m+tqZZi/h39bWBgA4efIkzp07V3Cf\nvJS/7GeV19bWWtd8IY6i0jGH6QLwjqq+q6q3AfwhgAMrHnMAwNuJ+38K4NckDz6aRAQOh8MqtFFV\nVYWJiYmCreRC+ce+/dfUegWA7u7ugpu/BNITMD8K4G9tX7+X+N6qj1HVBQA/AVCWhtfOODMpPjc3\nh+bmZkxMTHDRh/LGakdVmJoJFy9eLLhCwlm36CMinQA6AeCJJ57IcGuSY1YQOzo6MDw8XPBnN1P+\nMLt9zPVstgJ3dnbi7NmzeOqpp9DR0ZHJJm6odPQwvw/gF2xf/3zie6s+RkQeArATS4s/H6KqPlWt\nUdWaxx57LA3NW1+m9NXg4CBOnz7Ns5spr5j5S5MuZ+piVlZWwufzob29PdNN3FDp6GHOAvi4iHwM\nS4Hx0wD+04rHXARwGEAYwKcAfE2zeU/mA1i5gsizmynf2OthmgDZ29uLwcFBrpI/KFVdEJHPAbiK\npbSioKrGROQ0gOuqehHAMIDfF5F3APwYS0E159nPJDd/TkxMAABTiihv2OthigjcbjcAoK+vD9XV\n1QV1rbP4RgpM4WBTqch8ApvhOPeTUz6YmZlBQ0MDAFgdAlOlKFevcRbfyAD7meQmaALASy+9hIce\neihnK7kQ2Zl6CQYT12lN7GeSm2rUu3fvxqZNmzA0NMTkdcobZvEHABPXae1M2oXZb+twOJi8Tnlh\nZdJ6MBiEqmJ+fh59fX0FeRAaA2aamOIbtbW1iMViaGhosApzEOWi1ZLWR0ZG0NfXh71796K3t7fg\n0uc4JE+Ruaj6+/vh9XoBAF1dXVhYWMhwy4hSc7ekdTNfPzg4WHDTTgyYKbJfVNXV1aipqYGqQlUL\n7mKi/GLmLYPBoJUFYv4cHx/P2RXylJj/3Nl4e/rppzXXRCIRLSsr07KyMo1EIpluDlFKIpGIVlRU\naCAQ0MXFRQ0EAlpeXq6BQEDj8Ximm7dmWMoRf+CYxB5mitSWvD47O4uamhqrlD97mJTLTJAYHR21\nMkIKOWkd4KJPyswc5sjICBoaGqxJ8fn5eczMzBTcKiLlj2g0itbWVmv/uFkxL9Rq6wDnMFNmL++m\nqnA6nXC73Xj55ZcRj8dx9epV1NfXZ7qZRA9s5fy82k6OLMSkdYA9zJTZzzw5d+4cHA4HgsEgnn/+\neWzatCnTzSNKG3PCwNmzZwsyaR1gDzMt7BWLAGB8fBw1NTWorKwsyGEL5TYzL6+qaG1tXVYrwel0\noru7G9XV1QU5cmIPMw3sZ56YifKRkRGcPHmy4BJ7KffZE9bHxsbQ3t6OgYEB9Pb2IhQKFfS8PANm\nmtgXfxobG/Hiiy/izp07mW4W0QMzc5f20VFVVRUGBgYQDAZx7ty5glsdNxgw08RsjWxra8PQ0BC2\nbNmC1157DQAK+hOZco8p2zY7O2t1AlpaWiAiGB8fh9vtLsgFH4BzmGkzOzsLr9eLqqoqVFVVYXJy\nEgDQ2tqKsbGxgv1Eptxjzy22b400VdYL+VpmDzNNzDDGpF6YuczTp08X7Ioi5SZTMNgUjzEJ64OD\ng+jr6yvoeXkGzHUyPz+PpqYmHD9+3CqLRZRLYrHYsrJuVVVVBZuwbnBIniZm0Wd0dBSTk5OoqanB\nu+++iy996Uvo6urCnj17CnooQ7mjrq4Ok5OTqK2thYigt7cX7777LoLBIMbHxwt2/hJgwEwb+8qi\niGBmZgbBYNBa+OGwnHKFWfSx71wLBAI4cuRIwV/HDJhpYi4yYHnRglgshpMnTwIAOjo6CvrTmXKD\nPXG9qakJqopjx44hGAxi//79BT1S4qmR68B+mmRXVxf+8R//EVu3bsXly5cL+mKj7LZyh48pGGxi\nRC6fErnSWk+N5KLPOliZk7lt2za8+uqrBT1ZTtkvGo2iubkZsVgMFy5cALA0n+lwOHDw4EGrbkIh\nY8BcByYnc3Z21srJNOeVZ3OPngqXmUYyVdXn5+fR0tKCYDCIeDyO0dFRfuCDAXNdmAUgYClxHQBG\nRkawb98+HoxGWcnUvqyqqsL4+Li1f/zkyZNoampi7zKBiz7ryAzNVdU6GC0Wi+XNPBDlD/upp7Oz\ns4hGo2hvb4fT6SzYM8hXw4C5DlaeJHnhwgUMDQ1hcXERJ0+eRFVVVUGWxqLsZaaRgKVTTwFgYmIC\n8/Pz8Hq9zCNOYMBcB6tVqvZ6vTh9+jQWFxc5j0lZx37NVlVVWd83dTA5f7mEc5jrwKRfOBwO60K7\ncOECVBXxeDzDrSNazl5sw+wTr6urg8vlwsDAANrb2zmFlMCAuc7sB0l5PB6ICCYmJhg4KeNU1TrY\nzF7LtaGhAcFgENFo1Mr2oCVMXF9n8XgcIyMjaGtrw+zsLC5evIihoSH4fD50dHRkunlUwEygvHDh\nAubn57F7926ICGKxGLxeL8bGxqzTBPKth7nWxHXOYa4ze51MEUF/fz+eeuopOJ1OqGreXYiUO+wl\nCbu6uqCqOHfuHNxuN/bs2ZOXgTJVDJjrzH5RNjc3Y2BgALt378bevXtx5coVPPPMM5luIhUgM29p\n5tgnJiYwNzeHvr4+68OdPoxzmOvMDGkAWKfvTUxM4P333y/4A6Uoc0zq28zMDKLRKOrq6tDR0YGx\nsTHMzc2dhhdJAAATzklEQVShubm5oAsF3w0D5gYwCz/V1dUYHx9HY2Mjtm/fjjfffJPFhSkj7Jsq\n7EWCY7EYU4nuxewhzcbb008/rfkgHo9rJBLReDxufR0Oh9Xv92tZWZmGw+EMt5AKTSQS0V27dmk4\nHNZAIKAVFRUaCAS0rKxMA4GAda3mKwDXdQ0xKaUepoiUisifichfJ/4sucvjFkXkm4nbxVReMxfZ\ny2Lpirmj27dv48aNG+xl0oZyuVwYHR21igSPjY3B6XQCAOcw7yHVIXk3gL9Q1Y8D+IvE16v5QFX/\nReK2P8XXzGnmgKlgMIju7m4sLi7i5ZdfZlEO2jDmQxsAmpqarOIa9fX1mJyc5BbIe0g1YB4A8Hbi\n/tsAmlN8voJRVVWFUCiEF198EQ899BBisRh7mbSuVBWRSATDw8NW3ctQKLRs8ZGFYe4t1YBZrqo/\nTNz/7wDK7/K4LSJyXUQiIlLQQdUcMFVfXw+Hw4GRkREcO3as4I8vpfUXjUbR2NiIrq4uuN1u9Pb2\nYn5+HiKC1tZWXn9JuG8epoj8OYCKVX7Ua/9CVVVE7tZF+kVV/b6IVAL4mojcUNXv3uX1OgF0AsAT\nTzxxv+blHPsBU/F43KrMXllZiXg8zmR2SjszBK+trcXExASApTnMyspK9PX1YXR01CrtRvexlpUi\ncwPwHQC7Evd3AfhOEr/zFoBPJfP8+bJKvppIJKJlZWXWqmRpaalu27ZNp6enM900yjORSMRaBV9c\nXLQyNuzZGhUVFRqJRDLd1A2DNa6Sp7rT5yKAwwDOJv786soHJFbO/7eq3hKRRwH8KwDnUnzdnOdy\nuZZ92r/77rs4c+YMJiYmUF9fz14mpYX5j26OngAAr9eL0dFR6xrzer3Mu0xSSsU3RKQMwB8DeALA\n/w/g/1bVH4tIDYBjqnpERH4ZwHkAcSzNmf4/qjqczPPnQ/GNZGhiMj4UCuHNN9/Eq6++CrfbzaBJ\nKVutwIbD4YDq0hG6oVDIKkNYSNdbRopvqOrfA/i1Vb5/HcCRxP1pAHtSeZ18NzMzg6amJpw5cwa3\nbt2yzjFn0KRUmR09wM8qqU9OTlo/N3PqlBxujcwiIoItW7agsbERJ06cYG4mpcxUyzJ1WE0Kkcvl\nYs7lGrBaURYwqUa1tbUQEXR1dWFxcZEHptGaaKIwMPCzHqbL5YLD4bCG6GNjYwyWa8CAmQXsqUZV\nVVW4ePEiQqEQTpw4AVVFR0cHgyYlzeRbAsDQ0NCyeqw1NTVMIUoBA2YWMVWN+vv74ff7cevWLRw/\nfhzV1dU8ZZKSYlbFQ6GQVVrQHMTX3NwMt9uNYDCI6upq9jDXgHOYWcQUG25vb8fExAS+8IUvYNOm\nTQiFQjwDiO5LVREMBtHS0oL5+Xm4XC5Eo1FrznJgYADDw8NMIUoBe5hZxL5iaf7s7OzEmTNnoKrY\nv38/5zTprqLRKHp7e7F371709i5txLOvjLvdblRXVxdcClE6sYeZpUxF7MrKSnR1deH111+3qhyl\nkjtL+cvlcmFwcBCXL19GR0cH2tralq2MAyyukSoGzCxlhlBerxdPPfUUgKXe5smTJ5luRMvYV8Xd\nbjdeeeUVDA8P46233kJdXR0cDgeLa6QJA2aWEhG43W6Mj4+jqqoKmzZtwlNPPWUdg8peZmEzQTIe\njyMYDFrHTABLQXNwcNCqgGXmxjlvmQZr2YC+Ubd8Lr7xIBYXFzUQCOidO3fU4/FoSUmJ+nw+DYfD\neX+UAH1YPB5fdqxERUWF9vT0aHl5uXW8hCmswWtkdchQ8Q3aAGa3BgD4fD7cvn0bv/3bv42tW7dy\nt0YBikaj1kFlTqfTOmpCVdHb27vsiInW1lYmqacRh+Q5wJ5uFAqF8MILL2Dz5s04c+YMvv3tbyMc\nDnOIXiBMYBwbG0NVVRUOHjyI+fl57N+/H+fPn8fg4CAAoKWlBQA4FE+3tXRLN+rGIfmHRSIRLS8v\n156eHvX5fApAi4uLC6qWYSGz17ZcWFiw/pyenla/36+Li4sfOqWUPgyZODWSNp5JHQkGgxARFBcX\n4/nnn8eNGzeY3J7HNLHIU1tba9W2fOutt+D1enH9+nU4HA6cOnUKs7Ozy04ppfTiHGaOMavn1dXV\nVrGO48eP4/3338f3vvc9DAwM8D9KHopGo2hubsbAwADa29utf/+qqiprJw+H3+uPATMHrSzWcenS\nJYRCIbzxxht48skn4XA40N7eDoeDA4hcZ3qWqj+rmm4WdcwHY2Njo1Wpn9YX/0flMFOso6ioCAcO\nHMDi4iJ+67d+C0eOHMHIyEimm0cpMIFyZmYGjY2NVvWhsbExALDyLpWLfRuKATOH2YdhdXV1eO21\n17BlyxZs27YN8Xgc4XAYkUiE/6lyiD1QNjc3Y25uDqFQCENDQ+jt7UUsFrN2gfX19UFEmFq2kday\nUrRRN66SPxiTrOzz+bS4uFi3bNmixcXFGg6HM900SoI9IX16enpZMrrZvFBaWmolo3MlfO3AxHUS\nEWseS0TgcDggIpibmwPAwgvZzp6QLiIIBoPo6OiwKg85nU7r349n8WQGA2Yeqqurw+XLlxGLxQAA\nJ06cwO3bt/GFL3yB1duzjOrPCmeYKRYAVvpQW1sbKisr0dfXh7GxMUxOTnIlPIMYMPOQ6V2eOnUK\no6Oj+OxnP4szZ87gxRdfhNPpxM2bN7mKnmGqahX3NQs6k5OTEBG0tLSgv7/fOke8ra0NAKxzeSiD\n1jKO36gb5zDXzj7Htbi4qD09PVpWVqaHDh1Sh8OhPT09nP/KoEgkort27dLp6WkNBAI6NTWl4XB4\n2e4de4GNXbt2cTdXGmGNc5gZD4r3ujFgpo9ZUHjsscf02Wef1R07dui1a9esBQXaOPF43AqU09PT\numvXLiso2oOj+dBbXFzkAk+arTVgckheIMwOIQA4fvw4fvrTn8Ln8+H3f//38e6773KH0AZQ2zC8\nqakJABAKhdDf34+2tjZUV1ejpqYGAKxdXCuPLKHMYsAsICZoOp1OzM3NQVXxJ3/yJ3jjjTcQj8dR\nWVlpnVDJ4Jk6EyBra2sxOzuLeDyOpqYmhEIhhEIha1HOzFW63W5Eo1F4vV6e6pit1tIt3agbh+Tr\nx171pqenRwGoiGhxcbH6/X4Wnk2RmQIpKyuz5iI9Ho+WlpZa729FRYX6/X4NBALL8i05/F5/4Bwm\nPYiVi0I+n0+7u7v1/PnzWlxcrNu2beMc5xqY9zUcDmtpaakWFxdb85WmLF9FRYWGw2ErkNrvc2Fn\nYzBgUkrMqm0gENDt27crAH322WdVRLiingR7oDQ9x+npaZ2entZwOGz1HM3q98q6ldy5s7EYMCkl\n9hXZ6elp9Xg8WlFRofv27dPHH39cu7u7rQK1tMT+npke4sohtgmg5qwd88HEnmRmMWBSWpk5uPLy\ncj18+LA1x+nxeJb1nAqxR2T27NvzJEtLSz803C4vL+ccZZZiwKS0Mz2ihYUF7e7u1m3btumOHTu0\nuLhYN2/erMXFxctOr8zXYaW9J2mCYVlZmZaUlGhPT49OTU1pWVmZTk1NLRtuc44yezFg0royvSoz\nXAegW7du1W3btmlZWZl6PB71+XxaWlqq09PTOR84VxtuBwIB3blzp5aUlKjH41G/329VFjJJ6Pbh\nN+cosxcDJm2YxcVF9fv96vP5tKSkRA8dOqQAdPv27bpz506rLJnH49GpqamcGL6v1ou0D7dNQCwt\nLVWPx2MNt83Z3/b5S/Yksx8DJm04+5Dd7/frtWvX1OPx6OOPP27Ne27ZskV37NihZWVl6vf7lw1b\n7c+xEcHUviVxYWHhrgHS9CLtw20T9M2Ktz0wrhx+syeZ/RgwKePsyfBTU1O6bds23blzp7XCXlpa\nqlu2bFERsb5nhvE+n89aTFotqJpgda/h7crvr9ZrLC4utlKlKioqrKIkphc5NTWlxcXFVpZAMsPt\n1V6bshsDJmXcyjm7lcnZfr9fH3nkEav4h4jotm3bdPv27VpcXKw7d+60KsWLiBWoenp6tLS01Orp\nmWBn7y2auVUzh2p/bRMUS0pKtLu7W30+n05NTanf79eSkhL1+XxWypRJODeLWRxu56eMBEwA/wFA\nDEAcQM09HrcXwHcAvAOgO9nnZ8DMfSuDqEmx8fl86vP59Pz583rt2jUriPl8Pn3kkUf00KFDVm/P\npDMFAgH1+/0qInr48GH1+XxWb7G4uNhaiPL7/VpWVqY7duxQj8ejJSUl6vf71e/3L9ttY3q9pjdp\nPwqCw+38lqmAuRvAPwPwl3cLmACKAHwXQCWAzQC+BcCZzPMzYOYfewC17y6ylzQzaTtmeO7xePTa\ntWtW6o5ZVDJnF5mFpe7ubutrs1jj8/mWzUGuzIlcuWXR3jYOt/NXRofk9wmYzwC4avvaA8CTzPMy\nYOa3u9V7DIfDWlZWZvXodu3aZfUazTymSXGybzucnp6+62PC4bDu2rXrQ71E1pwsTNkcMD8FIGD7\n+jkAX77Hc3UCuA7g+hNPPLFObxdls9XyF1cGNPsWQ3N/ZTC0P4a9RLJba8CUpd+9OxH5cwAVq/yo\nV1W/mnjMXwJ4WVWvr/L7nwKwV1WPJL5+DkCdqn7uni8MoKamRq9f/9BTEkF1qdakORDM3LfX8bQ/\nhvU9yU5Evq6qNQ/6e/ctIKyq/3ZtTbJ8H8Av2L7++cT3iNZs5TGzqxXb5VG0lG4bcQTdLICPi8jH\nRGQzgE8DuLgBr0tElFYpBUwRaRGR97C0sDMpIlcT3/85EbkEAKq6AOBzAK4CuAngj1U1llqziYg2\nXkpn+qjqGICxVb7/AwD7bF9fAnApldciIso0ngpPRJQkBkwioiQxYBIRJYkBk4goSQyYRERJYsAk\nIkoSAyYRUZIYMImIksSASUSUJAZMIqIkMWASESWJAZOIKEkMmERESWLAJCJKEgMmEVGSGDCJiJLE\ngElElCQGTCKiJDFgEhEliQGTiChJDJhEREliwCQiShIDJhFRkhgwiYiSxIBJRJQkBkwioiQxYBIR\nJYkBk4goSQyYRERJYsAkIkoSAyYRUZIYMImIksSASUSUJAZMIqIkMWASESWJAZOIKEkpBUwR+Q8i\nEhORuIjU3ONxfyMiN0TkmyJyPZXXJCLKlIdS/P05AK0Azifx2H+tqn+X4usREWVMSgFTVW8CgIik\npzVERFlso+YwFcD/KyJfF5HODXpNIqK0um8PU0T+HEDFKj/qVdWvJvk6n1TV74vI4wD+TET+P1X9\nq7u8XicAE1Rvichckq+RTR4FkKvTD7na9lxtN5C7bc/VdgPAP1vLL903YKrqv13LE694ju8n/vyR\niIwBcAFYNWCqqg+ADwBE5Lqq3nUxKVvlaruB3G17rrYbyN2252q7gaW2r+X31n1ILiLbRWSHuQ/g\n32NpsYiIKKekmlbUIiLvAXgGwKSIXE18/+dE5FLiYeUAronItwBEAUyq6pVUXpeIKBNSXSUfAzC2\nyvd/AGBf4v67AP75Gl/Ct/bWZVSuthvI3bbnaruB3G17rrYbWGPbRVXT3RAiorzErZFEREnKmoCZ\ny9ssH6Dte0XkOyLyjoh0b2Qb70ZESkXkz0TkrxN/ltzlcYuJ9/ybInJxo9tpa8c930MReVhE/ijx\n8xkReXLjW/lhSbS7TUT+p+09PpKJdq4kIkER+dHd0vtkyRcTf69vi8gvbXQb7yaJtv+qiPzE9p6f\nuu+TqmpW3ADsxlJu1F8CqLnH4/4GwKOZbu+Dth1AEYDvAqgEsBnAtwA4s6Dt5wB0J+53Axi6y+P+\nIQvaet/3EMDzAN5M3P80gD/KkXa3Afhyptu6Stt/BcAvAZi7y8/3AbgMQADUA5jJdJsfoO2/CmDi\nQZ4za3qYqnpTVb+T6XasRZJtdwF4R1XfVdXbAP4QwIH1b919HQDwduL+2wCaM9iW+0nmPbT/ff4U\nwK9J5vfuZuu//X3p0gaTH9/jIQcA/J4uiQB4RER2bUzr7i2Jtj+wrAmYDyBXt1l+FMDf2r5+L/G9\nTCtX1R8m7v93LKWBrWaLiFwXkYiIZCqoJvMeWo9R1QUAPwFQtiGtu7tk/+0PJoa1fyoiv7AxTUtZ\ntl7XyXpGRL4lIpdFpOp+D061WtED2ehtlumUprZnxL3abv9CVVVE7pY28YuJ970SwNdE5Iaqfjfd\nbS1gIQB/oKq3ROQolnrJ/ybDbcp338DSdf0PIrIPwDiAj9/rFzY0YOoGb7NMpzS0/fsA7L2Gn098\nb93dq+0i8j9EZJeq/jAxlPrRXZ7DvO/vishfAviXWJqX20jJvIfmMe+JyEMAdgL4+41p3l3dt92q\nam9jAEtzy7kgY9d1qlT1fdv9SyLyFRF5VO9RhjKnhuQ5vs1yFsDHReRjIrIZSwsSGVtttrkI4HDi\n/mEAH+oti0iJiDycuP8ogH8FYH7DWvgzybyH9r/PpwB8TRMz/Bl033avmPfbD+DmBrYvFRcBHEqs\nltcD+IltiieriUiFmd8WEReW4uG9P1wzvZJlW7FqwdL8xy0A/wPA1cT3fw7ApcT9SiytMH4LQAxL\nw+GcaLv+bEXxv2GpZ5YtbS8D8BcA/hrAnwMoTXy/BkAgcf+XAdxIvO83AHRksL0feg8BnAawP3F/\nC4A/AfAOlrbiVmb6PU6y3WcS1/S3APwXAJ/IdJsT7foDAD8EcCdxjXcAOAbgWOLnAuD1xN/rBu6R\n4ZKFbf+c7T2PAPjl+z0nd/oQESUpp4bkRESZxIBJRJQkBkwioiQxYBIRJYkBk4goSQyYRERJYsAk\nIkoSAyYRUZL+D5myyiACZE1yAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "times = np.linspace(20.*torb, 1000.*torb, Noutputs)\n", "for i,time in enumerate(times):\n", " sim.integrate(time, exact_finish_time=0)\n", " x[i] = particles[1].x\n", " y[i] = particles[1].y\n", " \n", "fig = plt.figure(figsize=(5,5))\n", "ax = plt.subplot(111)\n", "ax.set_xlim([-1.5,1.5])\n", "ax.set_ylim([-1.5,1.5])\n", "plt.scatter(x, y, marker='.', color='k', s=1.2);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That looks much more like it. Let us finally plot the orbital elements as a function of time." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA44AAAFACAYAAADziiJ5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXeYFFX297+HGRhyFBEBHRUUUREVAXNGMO8a1vBTTOua\ndvVddxV1FdecRVxFUVFUjJhAECTnnDMMQxjyMDnH8/7R1d3VPdVd1d1Vfau7z+d55pmqW7dunYp9\nzz3nnkPMDEEQBEEQBEEQBEEIRSPVAgiCIAiCIAiCIAjuRhRHQRAEQRAEQRAEISyiOAqCIAiCIAiC\nIAhhEcVREARBEARBEARBCIsojoIgCIIgCIIgCEJYRHEUBEEQBEEQBEEQwiKKoyAIgiAIgiAIghAW\nURwFQRAEQRAEQRCEsIjiKAiCIAiCIAiCIIQlXbUAKjnssMM4MzNTtRiCIAhCHFi+fPkhZu6oWo5E\nQX4jBUEQUgOrv48prThmZmZi2bJlqsUQBEEQ4gAR7VQtQyIhv5GCIAipgdXfR3FVFQRBEARBEARB\nEMIiiqMgCIIgCIIgCIIQFlEcBUEQBEEQBEEQhLCI4igIgiAIgiAIgiCERRRHQRAEQRAEQRAEISyi\nOAqCIAiCIAiCIAhhEcVREARBEARBEARBCIsojoIgCIIgCIIgCEJYRHEUBEEQBEEQBEEQwiKKoyAI\ncSHrYCl2F5SrFkMQBEFIcDbsLcaB4krVYghCyiGKoyAIceHSt2fj3NdmqhZDEARBSHCuGDEXlw+f\ng5mbD6LH05Pwx/r9qKmr923fmVeG4soahRIKgrPsK6oIeObjhSiOgiA4zsuTNqoWQRASHiIaRESb\niSiLiIYabD+fiFYQUS0R3RC0bQgRbdX+hsRPakGwj7p6xufztwMACstrcNdnS1FTx7jvy+Xo8fTv\nyBw6EZ/MzcYFb8zCTR8uVCytIDhDSWUNznplBl74bUPcjy2KY5w4WFyJbbmlqsVQTk1dPerqWbUY\nKcGktftQVVunWgwAwKg52apFEFzG375chi8X7lAtRsJARGkA3gcwGEAvALcQUa+garsA3Ang66B9\n2wMYBqA/gH4AhhFRO6dlFoRYGb96L4b9ug7nvT4DheXVGDUnG89NCN9ZfnGiZ6By0/6SeIgoCHFn\nf5HHTXvK+v3YsLc4rscWxTFO9Ht5Oi55a7ZqMZTT4+nfcdNH8R8FnLs1F2e+NA0V1e5QpJxmwbZD\neHDsCrz2+2bVoliCmfH9shzXKLqC80xZfwDP/LpetRiJRD8AWcyczczVAL4FcK2+AjPvYOY1AIL9\nly4HMJWZ85m5AMBUAIPiIbQgREpOfjkyh07E8p35+Mc3KzFm4U7k5Ffg2V/X47XJmyJq64p35zok\npSCoISe/HJe9MwcAcKC4CleMmIsVuwridnxRHONAZY10hvUs3xm/B9zLK5M2IbekKmWsvoXlnrkd\newsrFEtijcnr9uPxcWswfNpW1aKkFPX1jD0J8owI6AIgR7e+Wytzel9BiCvzsw4BAL5ZkhNQPn71\n3ojb2rCvGJlDJ+KZX9b5lFHpkwmJxsGSSizfmQ8AuOq9eQ22b9wXP6ujKI5x4NXfIxshE+JDeXUt\n/vndKuSXVasWJalgZmQOnYg3p1i3dnqDGOSVVjklVtKyKqcQ45bvjmrfD2Zl4ZxXZyA7RQZUBHOI\n6D4iWkZEy3Jzc1WLI6QguSWe34GiivDBbS44viP+NfB4S21+uWgnAOD6kQvxxI9rsLewAmt2F8Ym\nqCDEge2HytDvpem4fuRC5OSXG74XOw6VxU2e9LgdKYV49td1OK5jSww5OxNA4lh9kpkNBqMx45bv\nxk8r96BFRjpeuO7kuMqzt7AC9czo2q55XI8bT/43Mwv/uvwE1WIkPde9Px8AcMMZXcPW27ivGGt3\nF+GmM7v5yhZsywMA7CuqxLEdWzonpGAHewB006131cqs7nth0L6zjCoy8ygAowCgb9++MiFdiCvz\nsw7hralbAABTNxxosP2yXp1w4hGtcOYx7XFej46oqq3DrM25aNYkDXO3HkKrjHR8c98A7Mwrx0Nf\nrzA8xoxNBzF+9V4wAztevdLR8xGEaPB6A3Vr3xwXvTnLV37e64GR6S86oSOev/ZkdGsfv76kKI4O\n8MVCz8iWV3E0o76e8fA3KzDkrEz0P7aDg5KlJqHcU70/SgzG4uw8ZDROQ59ubeMi09mvzgAgP1qC\nOUUVNWjauBEy0tNiamewNtdHrzh6YVEPEoGlAHoQ0THwKII3A7jV4r5TALysC4gzEMCT9osoCLFx\n2yeLQ257cnBP3DbgaLTM8HddM9LTMO6BswF4vF0AgIhwcpc2mJfVDd0Pb9Ug8mRJZa0DkguCfXwy\nLxsvT9qE3/5+bsg6px3VFp/d1S+OUnkQV1UXUFJZi0lr9+OvXyxTLUpSEiogztyth3zLfxm1yGe5\nEWJDlBB7OfW/f+CWUYscadtrcQSA39fuw/Y4ursIkcHMtQAehkcJ3Ajge2ZeT0TPE9E1AEBEZxLR\nbgA3AviIiNZr++YDeAEe5XMpgOe1MkFICC49sRP+dsFxAUpjMEQEIvKtv/Ln3rjn3GPQrX2zkPvk\n5JfbKqcg2MFC7bd5VU5Dd+rb+h+Fv1/cHSNuPi3eYgEQi6Nt9Hn+D/z1vGPx0EXdG2wz60ezaY3U\npV5L3dGoEZnUFPR4lTeK02U7WFyJfi9Px8t/OgU3G1i0hNhYsSu2uTglBomwvfOIvDwwdgWIgO2v\niBXcrTDzJACTgsqe1S0vhccN1Wjf0QBGOyqgIETB98tykFdajSFnHx1QPuzqXhhyViaKK2vQtnmT\nqNuf+/jFmL0lF0NGL2mw7bzXZ2LlM5ehXYvo2xcEu6nT+nArg377M9Ib4fFBPdGmWWMFUnkQi6NN\nFJbX4I0IgoEI1jjzpWno9/K0kNtr6upjjgopFrLY2ZHnGbX9eWV0QVpS6R58u2RX3COZfmyQR9Mo\n9Ukq3QdBENRTUV2Hx8etwWuTN+HKEYHRIu865xg0akQxKY1e+h/TPuS2d6dLNG/BPRRX1mDOFk9g\nsh9XBPapNr84WKnSCIjimHBMWrsvpYLt5JVV41Bp6Kinz/66Hue8OiNs9DXpDMcXudyhKa2qxdCf\n1uLWj51xPQ1FvdwUQRBcRkFZNU58drJvXe8qP+IWe93wmjZOw5rnBuLNG0/FimcuC9iWK9G8BRfR\n+7k/DMuvOOWIOEtijCiONlNn0EOzU3F5cOwK/OmD5JmLV1/P2F9UGfX+szYfBACUVclkd7fgtjkj\nE1bvxcGS6J8xO6nXPgZ5YQZDVCDu8oIgxIuK6jrszCvDPoPf/hZN0vDjA2fhmlOPtP24rZs2xg1n\ndEX7Fk3w5OCevvLcYo/iOHtLbsiYCIKgilf+fApG3X4G3r/1dNWiAHBQcSSibkQ0k4g2ENF6InrE\noA4R0QgiyiKiNUR0um7bECLaqv0N0ZVPJqLVWpsfElGaVt6eiKZq9afqosfFlUe/W+X4MQ4UJ8/o\n2IdztmHAK9PjFpRDrI/Oc6EudLRqiitr8PdvVuKOTxvObUkljBRDs3dh+6EyvPXHZl+kQkEQBDt4\ncOxyXPDGLBwKsvTdcEZXrH9+EM44OrRbqV3oo94v2ZGP4dO2YMjoJTLlSFBKYXngoPKJnVvj5jO7\nYeBJRwQEflKJkxbHWgCPMXMvAAMAPEREvYLqDAbQQ/u7D8BIwKMEAhgGoD+AfgCG6RTBm5j5VAAn\nA+gIT/Q4ABgKYDoz9wAwXVuPOxNW7416X7c8FNFyqLQKJw+bgrW7iyzvM0+LbOqk+20qW1OceqT+\n88taDHh5ujONh6GksgaD352LjQZ5OY2o1WaY7y8ObXFclVOIpTskwGQwQ0YvwXszsgytAoIgCNEy\nc7Nn/taOPP+A8ZFtmuLNG0+NmwxNG6fhvB6H+daHT/PMcxw9fzuyQ6TwEgSnmLnpINbtKcKooHgE\nk/5xrut0A8cUR2bex8wrtOUSeMKHdwmqdi2AL9jDIgBtiagzgMsBTGXmfGYuADAVwCCtLW+PMR1A\nE/inVF0LYIy2PAbAdc6cmXPU17NPkUpE5m09hNKqWnwyr2EgDqeJRTVMXbUyer5atMtQGXPaOLVw\nWx427ivGW39ssa3N696fjxs/XGhbe+Fw1+ffj9F9q66tj78ggiAkNZU1flfQpTsKfMs9O7eOuyyf\nDjkTC5+8uIFb7MVvzY67LEJqc9fnS3HVe/PwwaxtvrKv7unvOqURiNMcRyLKBHAagODMrl0A5OjW\nd2tlocq97U0BcBBACYBxWnEnZt6nLe8H0Mke6eNHSVUt/u/TxZix6YDhdnEZC82B4krsLjCeW5eK\nly3eVlYXftsEHan4DgiC4D7G67yyJqzeiwuO74hv7xuA4Tf3ibssTdIboXObZnj6yhPR/fCWAdvq\n6xmLsvNC7CkI9qEfTPFyXZ8jcdZxHRRIY47jiiMRtQTwI4BHddbCmGDmywF0BpAB4GKD7YwQhiQi\nuo+IlhHRstzcXDvEsZ1kmsOoxyiAjV0d2j9/sADnvjbTnsaSCHKtjcsYp+R121VIlkGgksqapDkX\nQRCc5/FxawLW/335CRhwbAe0bqouxUCn1k0x7Z8XBJR9PDcbN49ahLlb3dlPFJKD1TmF6PnM5ICy\nv/TthuE3n4Y0l+Yvd1RxJKLG8CiNY5n5J4MqewDos4V31cpClftg5koAv8LjogoABzQ3V2j/DxrJ\nxMyjmLkvM/ft2LFj5CcVAXZ/cBKlfxZKzps+Cu8OuP1QGbo/NSlugXK8JMp1FfTITYsXoazJOfnl\nOOW5P/DFwp3xFUgQhISk1GDw+OQubRRIYs7ynR432oNJOpAvuANvvkYvjdMIr93QW5E01nAyqioB\n+BTARmZ+O0S18QDu0KKrDgBQpLmbTgEwkIjaaUFxBgKYQkQtdcphOoArAWzSteWNvjoEHqVSKZ/N\n32GpXqooLuv3hjc4/7xyD2rrGb+u2hO2XqRwwHKKXOwkxczfn5kDUuIkgjXs7albkDl0oqPHcOIq\neANbTN1g7FovCIKgZ/P+wD7A5EfPUySJOX9o37VN+21xlBOEAN6eugUPfLU8YGD2sl6dMOfxi9QJ\nZREnLY7nALgdwMVEtEr7u4KI7iei+7U6kwBkA8gC8DGABwGAmfMBvABgqfb3vFbWAsB4IloDYBU8\nVsUPtbZeBXAZEW0FcKm27hIi67aF6hpH0/lblVPourx6iUB9PWPm5oMJoXhEijdi3OR1+7ByV4FJ\n7dgIdf2qausM82UxGDV10Qdl+fe4NTjuqUlR7x8rB4sr8eRPa0MGljFSfEdM32rLsQcNn4PnJ2yI\nat/ke8oFQXALa3cX4ZTnpuD6kYFeRz2PiH9AnHCs++/lWPffy9GlbTNf2cdztyuUSEhWRkzfit/X\n7UdRRY2v7Mg2TdG5TbMwe7kDJ6OqzmNmYubezNxH+5vEzB8y84daHWbmh5j5OGY+hZmX6fYfzczd\ntb/PtLIDzHym1ubJzPx3Zq7VtuUx8yXM3IOZL9UUzZTnuvfn47zXQ8/9m7R2H/71w+o4ShRIQD/a\nRUra6PnbcddnSzF53X7VokRFqEv566o9uPit2Zi1+SDu/2oF/vTBAkflCHX9zn99Jk58dnKD8u+X\n7UaPp383nCxuhXHLd0e1n10MG78e3yzZhekb42+F27S/BKPnG3dyonm1xDovCIIdjJqbjZLKQDfV\n9i2aKJImNC0z0tEyIx3X9gmMsroqp1CRREKys1qXvm7wKZ0VSmKduERVTVVCWVtKq2qjyhPkhPXr\nwbErlHe2nSaa67ZLs9IeCJP/LxHx5tjccqAkLscL7ix4MQsAZWSNTAS8j5rZEycqmSAIqQAzN8hv\nnZHeCDP/daEagSzwyKU9Atave3++IkmEZCTroL//tWR7Pi49sRN2vHolBhzrziiqwYjiqIDbPl6U\nlHmCEjclQxJ344PuSeLeo9SAmaO2tgqCILiNcoNBwM0vDkabZuqiqJqRkZ6GC453NniikLq88NvG\ngPUHLzpOkSTRIYqjg4RSR/SmaTvaSzYiSclgpWYs162OPfMdo+Xv36xELwOXzFTBqrtjpFfYqhHZ\nW62g3B1pIxab5AUbOXsbej4zGfll1bYeNxq3U+/lUpXS5V8/rHY8aJAgCM4SHCX9/wYcpUiSyPj4\njr4B6zKgJ9hBfT1jdlAk1dOPaqdImugQxTGJYGbklkjoaDt54bcNuHvM0qj3n7B6r+GIqxnB0UEj\n3j/qPWNDvWoWmoXbnE/mbGbRvWfMsrDbf1npiSh8sCR6F+mi8hpkDp2In1dad0F3g1IdTLK70AtC\nslNSWYOr3psXUPb/Lj1ekTSR0Tgt8GPe85nJ+Od3q1z5rRQSh1W7A+fLPntVL0WSRI8ojgmE2fdq\n5OxtOPOlaQkVRVV/TsqUHYMD68tmbY5/AuAPZm3DcU9NQklljXnlKEjF375ih66l2/CmyQhIBxTD\n/V6zO3xgCAmiIwiCEb+sCpzb+MYNvdGhZYYiaSLDKAL2Tyv32O4NIqQWwQPYt/ZPDAu8HlEcXUi0\n89C8Cs6ewgobpYk/ds/DC1BOE6SP+82SXQCAwnJ7lR2zPIheDpVWKXXNSfa5mHY9h7UxpC6xyn1f\nLndsAEMQhORl4bZDvuWv7+2PG/t2UyhN5Gx6YRCOaN00oOysV2YokkZIdPJKq/DGlM2+9aaNG6Fp\n4zSFEkWHKI4O4ouwaFMnMRlH9heazPlyM5v2FztuxYr12Yk2J2jfF6fhL6MWxXZwGMsfiwtuLCTK\noAFgfV7hR3Oyoz5GlS7XpNmlefjrlVEfp7C8OmWsvYIgAD8u342Zmw5iZ145mjdJQ7PGaTixs7ty\nNlqhaeM0zH78woCy6jgM1gnJxwezsnDGi9MCyp6+MvHcVAEgXbUAgrvnhTlJos/HHDR8Lk7p0gYT\n/n6u7W07ZXGLpNnVDuWuevOPzeaVQhDpZUkkZTEarHoXGF2GS9+2Htk562Dk6YO89Hl+KgBgx6tX\nRt2GIAiJw2O63NDX9TkSw28+TaE0sZGRnoY7z87E5wt2+Moqa+oS0lIkqGP4tK0B689fexJuH3C0\nImliQyyODmLUWdNPrP511R7M3HTQentJ1gkOdoWM5/k9MW5N2O1Wrbtr90QXIddpzCbw232tI7Fs\nmkUWtYON+4odP0Y43PCuxiqDfndmxsJtefhu6S5fWY2MvAuCYIJZzt5EYNjVgZah60cuUCSJkIjk\nlVahujbw9/KOszLVCGMDojjGmd/X7fctP/LtKtz1+dIGnXy3KiOJxpYDJfh8/nYYqS/fLcvxLbuh\nky9Ez4Jth/BjUATOwe/ObVDPjbe5otr5wZO5W3Ntse7f8vEiPPHjWt/63Z+HjxArCELqMXre9oD1\n3t3aKJLEPogIjw86wbe+fm8xDpUmvkIsxIf1ewMHss/rcZgiSexBFMc4s8tCxNOvFu0yrROOmZsb\nWjHXpZAy+tLEDZizJRdXvDsXz03YoFqcsFTV1uHH5bslxHcM3Prx4gDXqETCGwQpGDtclb1t3P7p\nEvyspfkIiQOP34a9xfJcC0KK8fxvgb+5j112QoiaicWDF3ZH13bNfOtb9pcolEZIJOqDfgffuvFU\nRZLYgyiONmDqFmi2v22CeP59NLthwIzgUUAncaKrGMpiYhQl9OO523HH6CWoVRSEpa6e8cBXy7F8\nZ4Fp3bf/2ILHfliNaRutuyzHRJJHKw1HPHSYSBW+4B8UL1ZlNatXVlVrrR2TtzbSSzd7Sy6uGDEX\n3y7NMa8sCEJSEDz9ZN4TF6FJevJ0M4/t2NK3vC23FL+t2ZvwUewFZ9lfVIk7PwvMBX54UKTeRCN5\n3mgXYnW0PVUH5SM573u/MHaLs/LRNjuO3dFqc0uq8Pu6/Xhw7HLTugc1hThUugOnIulabXev7vp+\nvXgXxugCBBi2G3SxDXNkWjiuWVTRRH9lvNff6Xf/7FfVhI7fnusJpvPkT2tNanoC73wwK8tpkZIC\nIhpERJuJKIuIhhpszyCi77Tti4koUytvQkSfEdFaIlpNRBfGWXQhBXh9cmDgs67tmiuSxBneu/k0\njLm7H5qkN8L2Q+V4+OuVuHnUQtViCS7me920qGtOPRILn7xYoTT2IIqjDcRL8YvlOG7vaFOIZS+5\nxZUNylZZjPrJIZYN64a4yK9P3oTMoRNtDwgS6nhmilPWwVK8N31r2DpAQ4us1TQPXr5ctNO3/NTP\nazFs/HrsL2p4L0JhpqAeMLivoVidU4i9RdGP7iZSOhvrlsvQ58QMFFVYS4NhqODH6Xtz00cL8frk\nzSivtmYdTVWIKA3A+wAGA+gF4BYiCo7nfg+AAmbuDuAdAK9p5X8FAGY+BcBlAN4iIvn9F2wjO7cU\no+f7PZuuPvVIhdI4Q5vmjXHB8R1RXVvvO9d9hdZ/w4TU47c1e33Lr9/QG53bNAtTOzGQHw4XEEuH\ntqK6Dttyow+Vr5JozvvNKZuxIMuTVLjYYqc4Ej6YtQ1Aw07zZ/N3AHAukmSkLo43j1qEt6ZuiTo/\nXixKQSzpNIDAgYH+L0/HtA0HLO137fvz8eyv6yM6VvAzVllTh5cnbbTswuk0evmYrb0R+4sqfbkw\n7Rq0UuTVDaBhgCAhJP0AZDFzNjNXA/gWwLVBda4FMEZbHgfgEvKMHvUCMAMAmPkggEIAfeMitZD0\nvPXHZlz8VmB6nz7d2iqSJr7U1rN8w4SQ5ORXoHObpjivx2FJk8JFFEcbiFefy6hb+bevluOSt2aj\n3qTnl+hBKrzS/29mFm79ZHFE+74zdYvlur+t2RdR20ZEoxBHenuqaox/qMqqasNa8SJVUGNO6WDB\nVdXIDfmFic4ENRq7eBdGzcnG+zPtdY2srKlD5tCJmLTWEzXZ6jMQyfWdtuEA9hVVYMAr0yN6phse\n07nn84XfNuDKEXNTeSqtk3QBoJ80ulsrM6zDzLUAigB0ALAawDVElE5ExwA4A0A3o4MQ0X1EtIyI\nluXm5tp8CkIy8t6MwO9pk/RGuPucTDXCxIEpj54fsD5BZ1USBAB49td1eH7CBlTU1OH2s47Gl/f0\nVy2SbYji6AYsdMp25pUhO7esQfm8rZ4f9lBBNiI4RNzRi8yIXOGyWnvBttjzBlaEUNTCYcUt1Ci4\nTyxcP3IB+r88XSeDPRzUKaPxGIMYF5ReI2qCZK3VLMZ2B06KxOVWT2CuxPB152875MuJNmervR16\nvTI5Y1NDC7DVd/PTeduxfm9xRG9yIrkQJzCj4VE0lwEYDmABAMOPGjOPYua+zNy3Y8eOcRRRSBZW\nPzvQ9t82N3F8p5bmlYSUpb6e8cXCnT535lZNGyuWyF5EcbSBUCP4ds4buuCNWbjqvXkNygM+zkny\nnc4rq47r8ZxShKx0iO22BG9yKET4fV+aB/pxI8GKmVtUFLPb7uSrbGYFXp1TZFDqPJHOv01B9iDQ\nSthVKzOsQ0TpANoAyGPmWmb+f8zch5mvBdAWQPRma0EIQ7MmyeGSFwoiwk19u/rWHx+3RqE0gtsI\nNjSUVrpjaoxdiOLoAuzonpkHfdEvM/KDlLPvlu7CLIP8j04SSubPTSJ3xoNog+hES6jB2VCHsXr0\ngvLA+xyt9bVQ106iWojcKHWw1d2sbvBjYvR8RDPQbxp52I0XL/VYCqAHER1DRE0A3AxgfFCd8QCG\naMs3AJjBzExEzYmoBQAQ0WUAapnZ3UluhYSgNmje/1NX9FQkSXwZdvVJAetF5fbHXBASk+A0VBf1\nTC7PDVEcbSBUn8ofct9dva5vluTg9BemYssBv3XqiR/XNsg1EwlTNxzAtoPOBelhhi8gSKIQiwXF\nauffrNrcrYcCcmutNohEW1/PISNaxqokGu0dz9ch+Fjxsmk5Z8W2fuyQOSKjbtlZEnVAIl5ocxYf\nBjAFwEYA3zPzeiJ6noiu0ap9CqADEWUB+CcAb8qOwwGsIKKNAJ4AcHt8pReSkX1FFej+9O8BZfed\nf5wiaeJLi4z0gPXXp2xSJIngNl74zT8m9+/LT0DPI1orlMZ+0s2rCE5jRyeTPZMEQ2+HR0H474T1\nmK9ZnbJsVPT+GiLPYjgiVahv+HBBTPuHlCOafWzq48baTMhBC92Gqtr6BtG89Nfu3elb8e70rVj9\n7EC0aW7BFz9B+/dWo5bGkwZRVS0+WFYU4PV7i42PyWyxhfjgsnE1V8PMkwBMCip7VrdcCeBGg/12\nADjBafmE1OKNybFF2E4mxi7ehbGLd+GzO8/ERT0PVy2OoIjKIDfVhy7qrkgS5xCLow1E0vFRmdA9\np6AcYxbutFVhNJTFoZ7gyl3W8jaGw0w27/ZE6cyyDdlBfl3lmSaVX24wt9TidfBaSBtUt3ghI73e\nVp+xUO9FPNSm2rp609yEoU4j3nElTF1VY9jXLhkEQXAXTYPmMj5wYWpYG71M/X/nNyiTCKupzRM/\n+ue7Xtm7s0JJnEMURxdghzJpPkfJfb0yvUSEKJQHO4WJkM8szMOsrfdrdSUh8i3Gqh+4z4bmLmKd\nI2oVI7fke79Yhl7PTrF8/Gjm1UaV+iWKY0d8DJMG1+8talAmT7IgJA719YyvF+/yrU9/7AI8MSg1\n5jd66dGpFV6/vndgoXzIUpZF2Xn4dZV/4OCKk0VxFEIQqvNmGLjCZluHexzOjBm/ei8+mZutWgzb\n8N5TvQ97KA6V+i14f/7A72abV1rlb88meewgXoMLbv1dLa6saRDoAQByS6rwxcIdEbc3a7OaHHhW\nrJU3jAx2+25Yx8nH4coR/gjRPmu1Cwe3BEEwJthD5biOqZmi4qYzA9OhylcsdXnl98B5rt3aN1Mk\nibOI4hhnjJRM2+Y4hj2uGlewf3yzEi9O3Gi4LWYX3yjOx6k5hWZs1bkH3/9Vw9QWoQYUTK1QVg5u\nohREmm/L6JifzN0eURtWeOrntbhj9JKY2gh1fYzOuPdzf+DfBmHVH/p6BZ79dT225drv4q1XlszS\nhTAaKlfRvNPMwLKdBUFt6+SIvMkGRPJIib4oCInHmt3+qSPL/nOpQknUc+mJnXzLP68MzpAjpAIP\nf70iIPjgT8wsAAAgAElEQVTgGUe3Q++ubRVK5ByiONrALyE+FEYdohKH8rlE5bIW4S65JVWGycET\nHePrEH1v1uy67i00Txbv7XdbnZMZC+Ha8G4xUy5/X7c/7P7R8PXiXZizxRmrnTdYVDBGP/oFWuqa\neEb1NZ2HHOdJkJHPQW1YNtMk3Y/oj4KQONz9uT8g3mEtMxRKop5zuncIWA8OkCIkP7+t2edbJgLG\n3ttfoTTOIoqjDTzx49qw2/UdopGztoXdHil+Ny+Tihx7X/PWjxfh7s+XocbAnS86rJ+5nVaJ32KY\nvB5vdzqzo4XSZcwGEjbtLwm7vUF7uvPOdsDyFg0HiivxhkkI9GCLnpdRc7Jx7FOTDJXHeBKYx9G/\ncv9XK7TtxvLHdEyLXg92X5nXTSIwiuVREBKP4X/po1oE5bRuGhiNPCe/XJEkghvo3Lppg0j2yYQo\nji5n2oYDWJQdXdL2YGLtmGUfKrNFDtU8/PXKqPcNdwkra+psj1hr7oIc3U2duFY/OuYfUSiurEFu\nSZXRLj5W724Y2CQUTioD//x+Fd6f2XAgRs9PK/wWRKNrFSrXYayYtsoW60VxnAPF4e+f1XbijiuE\nEAQhHDV19Rg0fI5v/ZjDWiiUxh386bQuAeuXvTNHrI4pTMumyZ3p0DHFkYi6EdFMItpAROuJ6BGD\nOkREI4goi4jWENHpum1DiGir9jdEK2tORBOJaJPW5qu6+ncSUS4RrdL+7nXq3KxitVMfzop17xfL\ncPOoRRaOFbksVuSrqavHV4t2Yvuhsohc9ZjZ9MMZfNqR9hvtiihqVzt//2YlLn17Nqpq7bLIAqZX\nxYLo3vOzYi0999UZOPOlaQ3KI50HaSeZQydi6I8N5x5W1YS/znsKK/D21C0BZdHoiZat+iFYkHUI\nG/cZ51Q0I5RFMhxmOVWNA+GEb/tQaeTKqNVjiL4oCInDuj1FAR4raY3cHqLPeRo1ogbRVe34ZgqJ\nQXBQvU/uOFORJPHBSbW4FsBjzLyCiFoBWE5EU5lZH45yMIAe2l9/ACMB9Cei9gCGAegLT79iORGN\nB1AF4E1mnklETQBMJ6LBzPy71t53zPywg+fkWow6ZRe+MTPmdj+fvwMvTTIObhOOD2dn47XJ4d0I\n3UzINA5herkLsg4BaPgRcUYOTRmMYD8rik9x0Bxcf15Ltd37b5fmBKxbkabGggJvpR1v4CKGZzCk\ntp7RMsP6p/PWTxYDAHa8emXD43PQsslgSoP1JNK6QkenZqUDF4Ig+AkeP+7SNjkjR0bKTWd2Q+Zh\nLXDTRwsBAD+v2IO/X9JDsVSC06zcVYA/fRAYpfyoDs0VSRMfHLM4MvM+Zl6hLZcA2AigS1C1awF8\nwR4WAWhLRJ0BXA5gKjPnM3MBgKkABjFzOTPP1NqsBrACQFenziFWrHbqbImqalC2I8/vZx/tMQor\nGiaFt9KWN6m8XQQfsqCsOmE7zVbmrXk7yqaWZCsWR4tukXZczpW7Ck3rrNtj3dU1HGa6hNH24LJI\nn6ErRszFycOMczPGE0fVKEWDBdEM1giCEF/0OYnH3N0P7Vo0USiNu2jfwj/X8a0gbxchObnlY3OP\nwGQjLnMciSgTwGkAFgdt6gJAb0rYrZWFKte32RbA1QCm64qv11xexxFRYHKdBCMnvxxT1htHqtTj\ns4ZE0bnSJyqNBPtcRK0TPO/utBemYtJa8+tjO5ZcQyMnlBJk5xzHWBSBaC0+Rse0KzqpWU7U4O3M\nRu7R1mUZNHwusnPtnedrdnwnguMYH8fm9mysLHqjILiDGZsO4M7PlgIARtxyGi44vqNiidxF98Nb\nYczd/XzrqoOvCc6TFtQ3mvzoeYokiR+OK45E1BLAjwAeZeboJvo0bDMdwDcARjCzN7v8BACZzNwb\nHgvlmBD73kdEy4hoWW6uswm6vZ+MaDplV46Yi7992TDfn+nBQm022D51Q3SpNazkjIwVs3laALB4\ne+RBg0xlt0FRi+b6hLS2hDieV5GL5HdJxU+Yo8c00WUnWIiea5fC1NCSaaIQ+uad2nP8SDCOoBof\nQeZnHULm0ImG20K7XUvnSxDcgD4Fx4Bj2yuUxL1ccHxHtNessCNnhw/eJiQ2lTV1KKv2x/JY8cxl\n6HlEa4USxQdHFUciagyP0jiWmX8yqLIHgN4y2FUrC1XuZRSArcw83FvAzHnM7DVLfQLgDCOZmHkU\nM/dl5r4dO7p3tCx4rlkw2w+VIXPoRFRr8+ni1fGzk3D9QStKrZERbOWugoaFpoIYFEVxOZ24A6Zz\nHG3oVIfTv+LVZ3/q5/ApbaLhjSmB6R+Yo0tJ4+T0usJyvyu4WUqM/UWVMVlMwxGP+8wMfLNkV8MN\nvuBDEjxHEBKF4BQUgp/uh7cE4PkNWr4zX7E0glNM1uWvPq/HYb4Bg2THyaiqBOBTABuZ+e0Q1cYD\nuEOLrjoAQBEz7wMwBcBAImpHRO0ADNTKQEQvAmgD4NGg43XWrV4Dz5zKhCCaTtu0IMXKqbQNTmGU\nYsSOzmvwJGW7caqDHbWrqoV5YWaWbyefDCcVkmj0uQaKl6LXoqaOccuoRRizcGfYenr5pm86aKx4\nRYjbvgVe3CmVIAgA8Nz49b7l/1x5YlLnqYsVfWT160cuVCiJ4CRNG/tVqPdvOz1MzeTCyaiq5wC4\nHcBaIlqllT0F4CgAYOYPAUwCcAWALADlAO7StuUT0QsAlmr7Pa+VdQXwNIBNAFZo7nr/Y+ZPAPyD\niK6BJ5prPoA7HTw35agMMmhHZ/vmUYsS0hfcyqmbz1uL5HixX2x/BNbI23JzZz7Sd8DoXJYpHA1e\nGDR4whbu0LjluwMLorGMG+xj5PIcadNm+T9DEqP7uCAIzlJVW4fPF+wAANx/wXG497xj1Qrkcpqm\nS4r0VOD+r1b4llPJAu+Y4sjM82BiFGBPj/ahENtGAxgdVLY7VJvM/CSAJ6MS1imsRlW1Qzkw2x7l\nIYwCkLjFYmEWHMUqdp+NHR1d75nFOifTSt1ETXQQ6f03coW8/dMldonjCKHnuNp8HBse2hd+82da\nyjPIYWamGOtF0OeAdcv3RhBSlTtHL/UtS9pGc975Sx+c/eoM1WIIceKVP5+iWoS4IsMiDlLHHLfA\nDvEMIGHXodxsSQiXU8503zgGrKm3Ik+MxxDiQ/Ct1CtPth4nxu2h0D+LH8yKPCiE/p3734wsf7k8\nwIKgFL13RI2NeYqTlSPbNgsY3NtbWKFOGMER9hX57+nlJx2hUJL4I4qjgyzfWYCvFoWfwwREF900\nOD1CcN/qUNCIPxvUsXachmV29ePGr44uHUgsXD8y/BzIWKwb/nyJ9vV07Yhg65NLxRxHB9tO9pzw\nPZ+ZHPKeFVXU4PP526O6vkYupfHQzdjsI6TbVloVPjiYIAjxITh9UnWtKI5WWDNsoG/5po9knmMy\nUVJZg7Ne8VuUW2Sk1nxfURwd5scVe0zrvD8zy7ROMMF95uAOZt8Xp0XcZjAvT9qIFdFEKbXIyCis\nEnpUKA6RKGrB1NbV4/XJm1Cgi6RpdgpfL9mJD2aFfj4sWUC9qR9CSL/9UJnWllH7ps0rwa1yeVm6\nI/K5kwzrngPZuWV4bsKGqN7PW0ZZS1gc6TWO9pbYOSdYEAR7KasOHMQ5/eh2iiRJLFo1bYw7zjoa\nALC7oALr9hQplkiwi0/nbQ9Yz0gXxVGwGbN+TzQ5YhvkjbPQ+YpUzxo1Jxvzs4yin9rrrunbJ+gc\nzI6jwuAUTiSi8HWmbTyID2ZtC4i4Fqq5bE2Z+2ZJDl6fvLnBdg76HxYTi6OTvPDbBmzLLXWkbTuC\n4zjFV4tij35qhWhG//cXVzYoc0OuxFASyBxHQVDDmt2FmLhmX0DZtX26KJIm8SiuqPEtD5+2RaEk\ngp2kJbu7kwmiOLqAIt3HxSoNHlsLfatIul/hOpLeLTM3HcQPy3JiO1BYGexpx/w4DQ8UleKr7RNq\n3mFtvbVOftbB0IpWTV093pyyGeVa0lkr7qfvz8yK+kdryY48MLPPKhkNj49bE/W+4QgOjvOfX+zP\nBRlPonkO7fr5CkjfYnLM+79cHtU8p0gCeOmfNxfotIKQklzzv/l48qfE/q6q5IpT/Fnipm08KPND\nk4RmTfwWxrXPDQxTMzkRxTFJsDtESrjOmnfbXZ8vxb9tUgpqDMyuZhIHz/O0AzvyYcbazz1Y0tAi\n5OXnlXvwvwDXZvMgPmMW7sTwaVujkmvdnmKszCmMYk8/4c4nFoJvv6mVj52zXllRbga/O9f24zrx\nDkzbGH7O9eT1+zFGC80P2JdeRr9t9pZcXbkgCG7gpwfPVi1CQjHwpCNw4Qkdfes/LNsdpraQCIya\nsw0vTvSkiZ/z74vQKoXScHgRxdFh6uoZc3SdILtwosOoJ2y0Tgd6ckNGuzstgg+Dc18elAsw1LWz\nw3ISPGIZiZtztO6ItXWxCW7R0OparLxrVhTSjfuKTdqIrl070B9l0/4S0/p6F6yF2XkoKq8x/S5Y\nff7c4DYrCKlOfdCPy0MXHYfTj5L5jZEyesiZvuWnfl6LLxfuUCaLEBuVNXV4edIm3/pRHZorlEYd\nojg6TDRuqFZoMMfRtNMWWfvRzLu0m7ilMolx/+tHLsTMzQf97cXYYCT5CSOJlOrE1bRyj6y66MYD\nu3J/OkWk98iu8aNYA+G8PGmjbTLUBn18RJEUhPgTHJTtXwNPUCRJYtMoKPHlM7+uVySJECveKUKp\njiiODhMcytop7LZMhLM4qrCCxAuzvIihtq7YWYAKX969EBbH6MUKyfUjF2CZUQRPGyOkxtpxd+od\niNTqzqYp6KPHDt0mmqi2timOBtcl3H0Pfk+CO0fRyeAh+HkRtdEPEQ0ios1ElEVEQw22ZxDRd9r2\nxUSUqZU3JqIxRLSWiDYS0ZPxll1ILJbv9EdsfuSSHo57OQmC27n786W+5S/u7qdQErWI4hgjOfnl\nYbc71mkOWreS7y8SBcDKHMdIMHPTA4CPZmfHfJxY2ZZbhrzShnnuzND/yIa65VsPmLsARkppVS2e\n+rlh8AJDBcnkes7cdDB8BQOs3KNgC5IdxOPRYGZ8vywHFdXxySm4ycI7EkwjmzpzEVscg+p3advU\nfB+T43i/T8Hu2GJw9EBEaQDeBzAYQC8AtxBRr6Bq9wAoYObuAN4B8JpWfiOADGY+BcAZAP7mVSoF\nwYiZm/1TbB66qLtCSRKfFk1SK11DMlJfz1ilxXt45JIeOP/4jiZ7JC+iOMbIfyeEdzswCn1vC0Ed\nRvOIhZHZWsJbHMNTUlWLzUFKUjSBQVSF4c8pqAidGiDEBv19HjF9q2Gd92Y0zMdopMxHqgsY6WWG\n1iuT6/lSFO6GVu5QvAZPzIhUAVm8PR+Pj1uDHXnhB4cAexTZvxjkVzQNEGXDcQFgS4SDGsG3tKYu\nemtu8H1pMKdWFEcv/QBkMXM2M1cD+BbAtUF1rgUwRlseB+AS8piKGEALIkoH0AxANYDIRyqElKCg\nzJ9ruFVGOpqkS1cxFmb9+yLVIggxMnK2P+94765tFEqiHvkaxEi8XFHNsHseULjWzI7lRDAgI5xy\nnKkLMycvVOdYL8vcrYciPqaZG1C4dBjBQQwA68qkFWJ9smINrqOKsirrlkan5uHFa37fFwt3BqzX\n14dXA4O3vjt9K/JKq0PU9u1kyZOhwRxH0Ry9dAGgz3+0WyszrMPMtQCKAHSAR4ksA7APwC4AbzKz\ngY87QET3EdEyIlqWmxufb7ngLsr0XhbioRozHVtlBKxLWo7EgpnxxhR/Pu1mKW5BFsUxRpxww7NC\ndK6q1ts3m+sXD0xFsOkHbUaQi6aKb7re5XlRdl6D7X/+YH7IfZ2+V0aKqRcrio1zcxwjq29Vipq6\nerw/M0s3Z9W9ODXvqM5srq/B5li9Kw5o+wcHU3LBpygZ6AegDsCRAI4B8BgRHWtUkZlHMXNfZu7b\nsWNs7ljP/roOH+lG6oXEoFQ3aBbu+y9ER4+nf8fKXQXmFQVXsC03MLe2XVNEEhVRHGPEDQqWdSKY\n4xhGeUqkM7ZCSWWgZSlcFNBP5243zEsYawdeP5o1fFpDV9diTUajqKBGv+tGz2W0982s37C3sCLs\ndjMlJFqc+nR/tzQHb0zZjPdnhu7wPjR2Bb5fmoPF2XnYvL/Etnci+FKpetdenrSxwXuhx05LqLel\nGz5cCKChhTrZvjcxsAdAN916V63MsI7mltoGQB6AWwFMZuYaZj4IYD6Avk4LvHxnAeZuPeQazxzB\nGqW6d98ox7IQOVf27hyw/sm87YokESKlujbwHUh1i7EojjGiKtNAcJJu02igJm5iwYSd4xin35GX\nJoafc5edG9p9Mxbq60Of4yfztuOhsSsalKscfzLqlP3nl3UNyk5/YWpU7ZvNd/1TGGso4J6oqlap\n1CyNpVWhU+lMXLsPj/+4Bn8ZtQiXD5/jiBxWcGrg87P5O8Juj+aWMqy5nUo6jpAsBdCDiI4hoiYA\nbgYwPqjOeABDtOUbAMxgzwXcBeBiACCiFgAGANgEh+nQMgPzsg7hFoP5u4J7eUD3G/fc1ScplCR5\nePOGUwPW2zZLvcTxiUpFTeAgau+ubRVJ4g5EcYwRp6wpZuSWBEb+vOCNWVhilJYhSj6Zlx1yW7zm\nHH25aKd5JQcwyzuYX2Yyl8sBwnWene5Ymw0iHCiOPAqtHczaHFkUWKuXyeuGUhfB3EynboFpOg5n\nDmuKk49cbXBUVecOlVBocxYfBjAFwEYA3zPzeiJ6noiu0ap9CqADEWUB+CcAb8qO9wG0JKL18Cig\nnzHzGqdl9j6fdv42Cc5SWVMX0L+4tf9RCqVJHpo1ScPsf1/oWx+7eBeqat0/HSLV+WRuNsYu2uVb\n79K2GdqkuNKfrlqAREeV/380Cbut7lJSWRPWTS/Ze3Jm1lunLGhLd+TjzMz2htvCHdHpR9CNBp99\nhRWOnDcR4E1JGNmgkFOao5KjmmL0juw0iT67/VAZ9hU1dGsOVn6DLY5/rD+Ak7u0TvlRXgBg5kkA\nJgWVPatbroQn9UbwfqVG5U6TddA/N2j6xgO45MRO8RZBiJBvl/g7yRedkLopB5zg6A4tcOxhLZCt\nBbsrraxFRsvUDrTidl4M8nw7oo156qlkRyyOMaLK4hjN3Eqru5hViyZqaCJhFgXUKCBStC6D+pZu\n/HBhVEqp0/NsS8JEF1UV8bK8OvKRWiuyMvuT2asYE4r0erpRqQ9HZU1Da37w+xT8/j3181pc87/w\n7tCCO3npTyf7lt+f2TAdkeA+dmqB2n568Gx8dlfqJjl3ip8ePNu3HG4eueBOPvy/M1SLoBxRHGMk\nUeaNe/I42iPsYz+sbtB2MlFXH/5aORnooTyKZPNOK47/+Galo+1Hg5PPnHfuZCTeBHaJk5MfaJEz\ne2dVvXtOPXOVNXWoUTVxXLCdC084HF3aNgMArNhViHV7ihRLJIQiv6wamUMn4rP5O3BsxxY4/ah2\nqkVKSto2b+JbvvDNWeoEEUypCBqgXvzUJQ1Sq6QiojjGSKK4qka7j8p2VWGWYsVo+5YDpQY1zQk2\nVAZ/qIJZajBXSOXghap7H81hrQSIIvLfE1XeBJGg6t47dWme+nmtROBMMtIa+b9yV703T6EkQjj0\n6SES4NMnCI7zt6+W+5afHNwTnVqLmyogimPMqOrkSFJs54jnHMfglkK5YHpF+nllcPT9REsJYw9O\nnrMvOE4E97moInQE1lgwO81YcydGi1PXf3VOYcqHOk82goMdCe6kRYY/5EXXds0USiII7mDOllzf\ncp9uMsfeiyiOMaKq0x6p7lJbx5ZGEa8fuQB7CsLn5QMQkMsw2dSW2joOe1JORlWtjqKTVVfPyC+r\nxrBfG6bgSFacfO3SGkV+DG8OwnhTXaumU+7U5a+uqzedYywkFu1bNjGvJChnf5H/N/31G3orlCS1\n8KZ/EtyNBGfzI4pjjKhSHCOd27QwOw/jlu82rbd8ZwEGvzvXtF6/l6ZHLYvbUemiGG3Qow9nb8OY\nhfFPX6LqUh0scSYFCDNA8AbHUf9cq5fAGKfe+eraenFVTTJG3d5XtQiCBR79bhUAoFv7ZujcRiyO\nTjLs6l6+5ad+XqtQEiEUetftG87oimZNJPqtF1EcY0SZq2oUh524dq/9ggDYYRKGP9F4fNwazIlX\n5Nig+xjNfa1nRrvmMqpvhtVL+/Qvnh9yNygwbh2UcSp+TXVtvbiqJhlHthUlxO3szCvzLX8ukVQd\nR+/2+NOKPShQkBtaCM9aXSCvRqoSJrsUURxjRFXf0g3WEC+Xvj1btQi2c6hUTVL7cLf1vxPWG5bX\nM6OzotxCCTXX1uI7U6O5SlYpcgNNBJy67zV17AqFXXCOYb+uw45DZeYVhbhx1Qh/0KJmjcWy4jTH\nd2oVsL6n0Hx6kBBf9LqiTJ8IRBTHGFEXHCeKfeTZdx0T1+4LWA/XIf9s/g7D8vr6wMiFQvLg1lfW\nqc9edW19wqQ4EqxzXo/DfMtjFu7EX0apmRMsGKPP1SuKo/O0yEjH01ec6FsvrnQmuJpgD2aR9lMN\nURxjRN0cRyWHFRwmWldVVY9DIj2HL03amHBukG4d6bQyXzoaquvqXeVNIdjDB7edHrBeXCGJz91C\ncEoxmcsVH+497xh8e98AAEBhuSiObmNbrt8rQiKqBmKqOBLR60TUmogaE9F0Isolov+Lh3CJgKo8\njlEFUXFADsFeolYcVQ1gKDlqdFTW1GP4tK2qxYiIf3y7UrUIcUcUx+SjVdPGmPSP83zrFTV1yn47\nhUB+XR2Y4ikjXewJ8YCIcHSH5gCAB8euwCnDpiiWSPDCzPh8wQ7f+l3nZCqTxY1Y+UIMZOZiAFcB\n2AGgO4B/OylUIqEqAmc0h5X+mPuJxnZYbyGxvZCYLNmer1qEuCPPcnLS68jWAevBbvqCGvTWrktP\n7AQimfYQL/RB7fTuwoJa1u8tDliXdyIQK4qjNyvslQB+YOaicJW9EFE3IppJRBuIaD0RPWJQh4ho\nBBFlEdEaIjpdt20IEW3V/oZoZc2JaCIRbdLafFVXP4OIvtPaWkxEmVbkjJVE8nyTZ9/9RDsIr8pZ\n1a1RP4XERSyOycunQ/ypOSqqJX+dajKHTsTLkzb61kfdfoZCaVKPpo3TAt4JwR2MXexJbdapdQbm\nPn6RYmnchxXF8Tci2gTgDADTiagjgEqTfQCgFsBjzNwLwAAADxFRr6A6gwH00P7uAzASAIioPYBh\nAPoD6AdgGBG10/Z5k5l7AjgNwDlENFgrvwdAATN3B/AOgNcsyBgzqjrOEoUrOYn2eXIqPYJbjysk\nL6I3Ji+XnNjJt7xxf3GYmoLTeH9rvFGkX/rTyWgkQdbizvnHd/QtbzlQolASgZmxt7ACWw+Uov8x\n7bH4qUvRrX1z1WK5DlPFkZmHAjgbQF9mrgFQBuBaC/vtY+YV2nIJgI0AugRVuxbAF+xhEYC2RNQZ\nwOUApjJzPjMXAJgKYBAzlzPzTK3NagArAHTVtTVGWx4H4BKKg31ZZbL4SCmqkAnYbifap0nVU3jq\n838oOrKQrIjFMTUIFSVaiA/B6YZu63+0IklSm8ZpjXDvuccAAC4fPkexNKnND8t24+xXZ2DZzgIc\nc1gL1eK4lvRQG4joYmaeQUR/1pXpq/xk9SCa2+hpABYHbeoCIEe3vlsrC1Wub7MtgKsBvBvcFjPX\nElERgA4ADgXtdx881k0cddRRVk8hJImUc6xa8tK5nmj7zOIyKiQL8igLgvM8Pm6NahEEjeZaJFtm\nYFtuKY7r2FKxRKnJil0FvuWjOoilMRThLI4XaP+vNvi7yuoBiKglgB8BPKoF2YkZIkoH8A2AEcyc\nHcm+zDyKmfsyc9+OHTua72CCRIYT7CRaBVA620KyIBbH1EF+P9UxfvVe3/KEh89VKIlQpQuW8cf6\nAwolSW2aN/Hb0jLSJS1NKEJaHJl5mPb/rmgbJ6LG8CiNY5nZyEK5B0A33XpXrWwPgAuDymfp1kcB\n2MrMww3a2q0plm0A5EUru1Xkd0+wk+hdVeVBFJID+aYmN00bN0JljaejvKewQuYQuQCxrqilqsav\nOL42eRMeuPA4hdKkLvllVb7lP50WPLNO8GIlj+OXRNRGt340EU23sB8B+BTARmZ+O0S18QDu0KKr\nDgBQxMz7AEwBMJCI2mlBcQZqZSCiF+FRCh81aGuItnwDgBkcB/+9RJrjKLifaB8n6WwLyYJYHJOb\n+U9c7Fse+I7M6VJB8BSbNs0aK5JEAIDqRArPn6TU1NXjl1UeK/wxh7VA+xZNTPZIXaxEVZ0HYDER\nXUFEf4UnUM1wk30A4BwAtwO4mIhWaX9XENH9RHS/VmcSgGwAWQA+BvAgADBzPoAXACzV/p5n5nwi\n6grgaQC9AKzQ2rxXa+tTAB2IKAvAPwEMtSBjzIirjWAn0Xaapa8tJAv/+WWdahEEB+nQMsO3XFFT\nJ0HbFPDejK2+5bOP66BQEgEA7jgrMDBRrSiScWfrgVLf8rj7z1IoifsJ6arqhZk/IqL1AGbCE2jm\nNGbeb2G/eQDCRjXVLIIPhdg2GsDooLLdodpk5koAN5rJZTdicRTsJOrgOOKqKghCAvLtkl342wXi\nmhdPlu/0BwEJVlqE+NPziNYB62XVdWjTzIpdR7CLK0bM9S3rB7eEhlhxVb0dHgXuDgCfA5hERKc6\nLFfCIHqjYCc788qi2k8M34KQOBDRKaplUMmX9/TzLTfPMB2/FhxEfjvcwZi7/e/E14t3YeE2x0N0\nCEJUWBnSuB7Aucz8DTM/CeB++PMlCoJgI0N/WhvdjjKCIQiJxAdEtISIHtTHEEgVzuvREb88dA4A\noLSyVrE0qU0ipRRLZi443h/l/7XJm3DLx4sUSiMIoTFVHJn5OmY+qFtfAqBfmF0EQYgz8tMvCIkD\nM58H4DZ4IoEvJ6KviegyxWLFlVO7tkFGeiN8Mjcb578+U3LRxonaunrM3epPby3BqNzDnWdnqhYh\nJVqHU10AACAASURBVFmzu1C1CAmFFVfVpkT0EBF9QESjiWg0gA/jIJsgCBaRIE2CkFgw81YA/wHw\nBDx5k0cQ0SYi+nOofYhoEBFtJqIsImoQAI6IMojoO237YiLK1Mpv0wWpW0VE9UTUx5kzswYRoUl6\nI+SVVWNXfrkvRYfgLJv2lwSsn9cj9nzWgj0MODYwUNHGfbakPhdMGDZ+vW/5kp6HK5QkMbDiqvol\ngCMAXA5gNjw5FUvC7iEIQlwRtVEQEgci6k1E7wDYCOBiAFcz84na8jsh9kkD8D6AwfBEFr+FiHoF\nVbsHQAEzd9faeQ0AmHksM/dh5j7wRDvfzsyrHDi1iNDnq6usqVMoSWrw+9p9+GBWlm+9aeNGknbA\nRQzs1QnpjfzxH9+cslmhNMlPXT3j+QkbUFbld5d/52al42kJgRXFsTszPwOgjJnHALgSQH9nxRIE\nIRL+O2GDahEEQbDOewBWADiVmR9i5hUAwMx74bFCGtEPQBYzZzNzNYBvAVwbVOda+GMQjANwiZZT\nWc8t2r7KeUAXTXV/caVCSVKDB8auwKS1/qD4/9dfIqq6iUaNCPedf2zAuuAcS3fkY/T87diipeL4\n+8Xd0bqp5DQ1w4ri6E2yVEhEJwNoA0BsuYIgCIIQHT8z85fMXOEtIKJHAICZvwyxTxcAObr13VqZ\nYR1mrgVQBCA4Ud9fAHwTSjAiuo+IlhHRstzcXCvnEjV6nXbwu3PD1BTspEvbZtjy4mA8feWJqkUR\ngrjuNP8rvf1QdFHWBWuk6RTzTq0z8OilxyuUJnGwojiOIqJ28IyCjgewAZr7iyAIgiAIEXOHQdmd\nTh+UiPoDKGfmdaHqMPMoZu7LzH07dnR+/ttHt5/hWy6vlgirTrFDp4S0apqOJumN0NAYLajm+E6t\n0CTN0zXPOliKksoakz2EaNlT4Bu3Q2aHFgGKpBAa0wRKzPyJtjgHwLHh6gqCIAiCYAwR3QLgVgDH\nENF43aZWAPJNdt8DTxRWL121MqM6u4koHR4PIX1CuJsRxtqogqaN03zLz/yyHm/dJGminWDIZ0v8\nyxK909WkNSJAm/JbXl2HVuI+6QiPfuef5v3cNScplCSxkMy7giAIghAfFgDYB+AwAG/pyksArDHZ\ndymAHkR0DDwK4s3wKKF6xgMYAmAhgBsAzGAtzwURNQJwE4DzYjwHW8lI9zs+ZR8qVShJcqMPAHJ0\n++YKJRHM0KdIqaiWoFHxoOcRrVSLkDCI4igIgiAIcYCZdwLYCeCsKPatJaKHAUwBkAZgNDOvJ6Ln\nASxj5vEAPgXwJRFlwWPBvFnXxPkAcpg5O9bzsBN9J3nlrkKUVdWiRYZ0TewmvZFfQW/bXCKpuhl9\nas21e4qQeVgLdcKkCOK2bZ2wX2dthHIAMy+IkzyCIAiCkJQQ0TxmPpeIShCYRYcAMDO3Drc/M08C\nMCmo7FndciWAG0PsOwvAgChFd4z6oPSNewsr0KOTjP7bjT5qbbMmaWFqCqqp02mOf/9mJa4+9UiF\n0iQnBWXVqkVIWMIGx2HmenjyRgmCIAiCEAPMfK72vxUzt9b9tTJTGpOVAce2D1ivrZestHazOqcw\nYD2zg7iqupkRN5+mWoSk55aPF/mWH9TlkxXMsRJVdToRXW+QC0oQBEEQhAghogFE1Eq33kqLeJpy\npKc1wpFtmvrW9XPxBHvYkeePqPruzX3ELc/lXNm7c8A6swym2M2m/SW+5Vv6HaVQksTDiuL4NwA/\nAKgmomIiKiGiYoflEgRBEIRkZSQAfSSYMq0sJfn0zjN9yy9O3Ih//7BaoTTJx668ct/ytX2CU38K\nbmf6xoOSqsZB9AG6BHNMr5bmQtOImRunukuNIAiCINgAsc6MoE0LSdmIMCd29ncpVuUU4ofluxVK\nk3y8NXWLahGECJn8qD/48b1fLMOj364KU1uIhSaiOEaEpatFRNcQ0Zva31VOC5UoiPuAIAiCEAXZ\nRPQPImqs/T0CwFXRTuPNzw+erVqEpET6KYlJzyNaY/srV/jW/9hwQKE0yUVtXWBErsZpojhGgunV\nIqJXATwCYIP29wgRveK0YImAzOEXBEEQouB+AGfDk49xN4D+AO5TKpFigiOpfr80R5EkyUWpzBlN\nWIgIbZs39q1vPVASprZglR9XBHo0NGssUYYjwYqafQWAy5h5NDOPBjAIwJXOipUY1AbHERcEQUhC\nPrjtdNUiJBXMfJCZb2bmw5m5EzPfyswHVculkpYZ6XhiUE/f+uM/rlEoTfJwoLgKAJDWiDDiFonW\nmWi0bupXHHcXVCiUJDnYlluKZTsKfOvX9jkSjRpJsKhIsDqnoi08yYQBoI1DsiQclTWiOAqCkPzI\nz6o9ENHjzPw6Eb2HwDyOAABm/ocCsVxD/6DUHEJsFJXXYMr6/QCAif84Fz2PkPAUiUZzXc7NjMbi\nUhkrl7w1O2D98pOOUCRJ4mJFcXwFwEoimglP/+F8AEMdlSpBKKmsUS2CIDhOs8ZpqKipUy2GoBCJ\n3m8bG7X/y5RK4VIyO7RQLULSMHXDAfz1i2Vo27wxju/UUpTGBOWc7of5UkfIdNXY0M/37X9Me4y9\ntz/SZX5jxFiJqvoNgAEAfgLwI4CzmPk7pwVLBIoqUlNxbCqjXoKQUkjeN3tg5glElAbgFGYeE/yn\nWj7VtG/RBC0z/OPZ6/YUKZQmsVmUnQcAKCyvwfFB80eFxOHJwX737ds+WaxQksRn/V5/JsH8smpR\nGqMk5FUjop7a/9MBdIZnAv9uAEcS0WlEdHR8RHQvxRWpOek8vZG8bIKQSojaaB/MXAfgHNVyuJWB\nJ3XyLa/eXahQksQmPc3/1p4gimPCIsqNfVz13jzf8iOX9lAoSWITzlX1n/BEeXsrxPYORLSamW+3\nX6zEoDhFXVVlHnH8ufPsTHy+YIdqMYQURSyOtrOKiMYD+AFAmbeQmX9SJ5I7qJdw5bbwzeJdvuWj\nOjRXKIkQK78/ch4GvzsXgMdS1r5FE8USJT5nHN1OtQgJS0jFkZnv0/5fFKoOEf3hhFCJQnGKuqpK\nBKr4c/7xh4nimII0Inek/ZE33naaAsgDcLGujOGZEpLS1LngeU90mBnFlX6PqEEnSwCQREavKN74\n4QJMf+xCdcIkCW2aNTavJBhiyQZORCcT0U1EdIf3DwCYeaCz4rmb045qh3vPPUa1GHEnXRRHJYy6\n/QzD8lO7SqDjZMUtlj6XiJFMfMLMd+n/AHyqWig3oLc4vjN1q0JJEpfaoNGmjHTJU5fIZKT7u+rb\ncsvC1BSsIrkbo8dUcSSiYQDe0/4uAvA6gGsclish6H54S9x+VupN9WwUp16kdFYDGRgibPRLfzol\nzpII8cItYzTyLtrOexbLUo4zM/0uZIdKqwIiIQrWqJQo2ElFU1FyYmazFpnWi1sGZRMRKxbHGwBc\nAmC/Nip6KiSXow9KQSeutDj1ZlPvykaHyu+ffHudJRW/L3oeuPA41SLYChGdRUSPAehIRP/U/T0H\nQHqHAIacnRmwXlCemlNCYmF+Vp5vWW+tEhKT4HtYJO9ExIyY7vdeeGJQzzA1BTOsfFEqmLkeQC0R\ntQZwEEA3s52IqBsRzSSiDUS0nogeMahDRDSCiLKIaI0WwdW7bQgRbdX+hujKXyKiHCIqDWrrTiLK\nJaJV2t+9Fs4tZlKx4xyvU5YRIWs4rVyEuw1pco+cxSWXV5UC2yH5gkA0AdASnvgCrXR/xfAM0qY8\nwd/9PQUViiRJXO7/arlvefpjFyiURLADIsLYe/v71sev2atQmsSE4fdcuCcFp5jZiRXFcRkRtQXw\nMYDlAFYAWGhhv1oAjzFzL3jyQD5ERL2C6gwG0EP7uw/ASAAgovYAhgHoD6AfgGFE5PVfmaCVGfEd\nM/fR/j6xIKPgYlzSZ27Af648UbUIAajU3SRQkrO45uoqEiTZBo+YeTYz/xfAAGb+r+7vbWaWCX0a\nH+nmc//z+1VYsO2QQmkSm67tJKJqMnBO98N8y/X1jOzc0jC1BT319YxJa/cDAJqkN0ITscLHhOnV\nY+YHmbmQmT8EcBmAIZrLqtl++5h5hbZcAmAjgC5B1a4F8AV7WASgLRF1BnA5gKnMnM/MBQCmAhik\ntbWImfdFcI6OkmT9Glfhlmt7XZ8jA9ZP7dZWkSTGqLxOEihJEKIig4hGEdEfRDTD+6daKLdw+UlH\nYPOLg9C1XTNsPViKWz+WxOeC8MtDnvSvw8avx8VvzVYsTeJQUF7tW66tq1coSXJgNapqbyK6BsDp\nALoT0Z8jOQgRZQI4DUDw178LgBzd+m6tLFS5GddrLq/jiMjQnZaI7iOiZUS0LDc31+IZhCbZRsTd\nhFvmdx3VoUXAeryCA+kJdy1UXqdEclVd/9/LI6r/0p9OdkgS66h41oxQJUUSj0v8AGAlgP8A+Lfu\nT9DISE8LGKTbcUiiSQqpzald26Bxmv+jKIGjrHGo1K84uiG9VaJjJarqaACjAVwP4Grt7yqrByCi\nlgB+BPAoMxdHKacVJgDIZObe8FgoxxhVYuZRzNyXmft27Ngx5oMmSr/mxwfOsmSeH2IhSmzclGWX\nXtx4BQeyitO3I1zzJ3dJnDhZLTJCpq015Lb+6iMmu0RvVDZA5pLTd4JaZh7JzEuYebn3T7VQbqNL\n22a+5f/NzFIoSWKwOqcQi7PzzCsKCQkRoY9uMCUnX+b/WiG3pMq3fF6Pw8LUFKxgxeI4QFO0huhy\nTt1tpXEiagyP0jiWmY0SG+9BYKCdrlpZqPKQMHMeM3ufjk8AGCe9S1E6tmyK9s3NA01cf0ZX9Mts\nH7ZOi4zA4H9OKVIu0898uM3KplKaPke1xXNXB09dFpINVc9YEnt0TCCiB4moMxG19/6pFspt6KNJ\nSnTQ8OwprMC178/HX0Yt8pW9f+vpYfYQEpFHLjnet3z+GzMVSpI4bD/kmQ/65T39AuZPC9Fh5Uu8\n0CCojSnk+cX/FMBGZn47RLXxAO7QoqsOAFCkzV+cAmAgEbXTguIM1MrCHa+zbvUaeOZUOk6i9Gus\nykkg015i+xZNMP7hc/Dn0z3ew04peG5xVQ2mkQv6L2/c0BsdW2UAiO8zeFT7hoEW2rfMiJ8ALmPC\nw+di9bMDHWvfLW+Aqu9c5zZNAbgvIJUNDIHHNXUBPEHnlgNYplQiF6LPXycBLcJTUikpGlKBoztI\nsCOr7C4ox0ezt2HK+gPofnhLnNv9MDRvEpnnkdAQK1fwC3iUx/0AquDpy7DmEhqOcwDcDmAtEa3S\nyp4CcBQ8DXwIYBKAKwBkASgHcJe2LZ+IXgCwVNvveWbOBwAieh3ArQCaE9FuAJ8w83MA/qHNw6wF\nkA/gTgvnFjNuVW6MsNL5s9pB7N21LTLS07R9CID9juNuVcpVzjv74LbTMWntPtzYtxvem+F13Yqf\nPNJ58/PmjafilK4eV90ubZthT6H9bkOJYnHr060tVuUU2t7uZb064at7+uPs4zrY3rZKmFniwVug\n/zF+I2xwAm8hEH1fpP8x7dHnqLYYeFInhRIJTuAdMBbM+esXy7Fxn2eG3I1ndE2Y31O3Y0Vx/BSa\nAgjAcjgiZp4Hkx4te2b2PhRim3duZXD54wAeNyh/EsCTVuWzi2R8Dq2eknditttcN50mXorj/Rcc\nhw9nb/OsaIe84pTOuOKUzgH1VF7+1Lrzfs7p3gHXn+6P1+VUkAK3XF+zAbJfHjoHmUMn2n9cIpyb\nhHNSiKg5gH8COIqZ7yOiHgBOYObfTPYbBOBdAGnwDJq+GrQ9A57B3jMA5AH4CzPv0Lb1BvARgNbw\n/JafycyVtp6YzfTVTZtYsC0P1bX1MngVgjpd1I/C8ho8OTjprPQCAq3wgOe+uy3uglsoq6r1LR/f\nqZVCSZILK1/gXGYez8zbmXmn989xyRKERHpdrcgaiRJSr3WWnfpmuSWiZDDx+kZfqVMQTzqydch6\nTosTbpSOzD2bE5qP7+hrWH7piZ3iM3rpkosbz1fxtetPid/B1PEZgGoAZ2vrewC8GG4HIkoD8D48\n+Y97AbjFYBrJPQAKmLk7gHcAvKbtmw7gKwD3M/NJAC4EkBC+jXeenelb1ncEhUBmbDrgW66qrVMo\nieA0PY/wK0H/myFBo0KhTxfWo1NLhZIkF1YUx5VE9DUR3UJEf/b+OS5ZouCSjp0ZRNbd3jKCRrQa\ntKWdtHeA0ykFz62XNm5BZbXjdG7TFIe3ahqmnrorlUiu2tFwWS9rrl5O3QO3XN3/3959x0lRn38A\n/zzXG3dcoRx3B0fvvVcpckdRUQQbIlZijcaowWgs+NOgMdZYE4klGrsJsSFGjIIV6YhIEQVEeu8H\nz++PndubnZ3dnd2duvu8ed2L2dnZmWdnZ3e+35nv9/naGUdFcgxY3pKZ74NSeWPmg4i8m/sAWMPM\n65j5KICX4RsLWW0c6jKKvw5ghJJvoArAUmZeomxvBzN7onZxx2kd/dMvf70hzJLJ7f4PvvdPN1Fl\noxWJ5+rhrfzTC3/a5WAk7rV172GsUw3h07409MV3ER0jFcds+Po2ViGG4TgSXaIVnAmE+yd0weUn\ntYy47Aml5phiXXYcV0qzKTuO8YRG9pFxo3zs2udOXRQY310zbK6dYbj0e2+yo0SUDaVzOBG1hO88\nG46R8Y39yzBzDYA9AIoBtAHARDSbiBYSUVB3j1pmj3Vspnvf/87pEDzh4XO6Ox2CsJC6WbKckfV9\nv2V/wONG+aEvvovoRCwBq4bgUP8ZGo4jGdhRrsvNCH8H0AijBVAioGF+FqaNbhdxWaubqjpVfpw/\nbXjY59NSCRVF1l/RNXpRwtE+jpSY/XzdwqmuKw+c3c2ZDSPxLsaFcDuA9wFUENGLAP4Lnb77JkoD\nMAjAJOX/M4hohN6CZo91LOwnCVQS2wnVBdzP1253MBL3Ul/fl0y05orq1gkRLbQqEK+yo4hTmBt5\n/EWzRJN5tfany7Kmqg7VSMoMNPPpUl43CO8bVwwIs2Tsan/43HSTT/uZEJKmoB/ArmPTLVngtJ/x\nQxZWLF3yli3FzHMAjIcv+/c/AfRi5o8jvMzI+Mb+ZZR+jQXwJcnZCOATZt6uNIt9F4AM8ieEB51Q\npak8dpyx68BR54JxKXW5NE2SB5kq2jZ3svc17CjYmbEJMrgeI5WA2vX4+zha9KWM5X3nZ9k/Rk/3\nivqRF4qBWypk6iiYGXecGvWwrlErybPvYonWVcNaosjGizVeoP0umvmzN6Zz44DHbk2KZSYiOgNA\nDTO/o2RSrSGi0yO87GsArYmoORFlADgHvrGQ1WbBN0YkAEwA8JGSvXw2gM5ElKNUKE8C8K1Z78dO\nSywY9sXrlm6UfZJMqjsF/mZu3uPq5MiOeG3BRv+09Pk1V7QVR/PzrXuctogzfVxH3eXi24Z9BSmv\nZ1XNNqFZr1uk+O/suuiWI4DR6uFAEuxuMwDcWN0OC/8wMuwytiVIivD8xQPtGQ5QG0e/FuaNq3he\nn2aB20r8eiMA3M7Me2ofMPNu+JqvhqT0WbwavkrgSgCvMvMKIpqujGEM+IbPKiaiNfAN9zFNee0u\nAA/AV/lcDGAhM3vmfL78zmr/tPRzDHbO01/4p6f0bxZmSZEI8jLTkKkalmbMI5/izCc+czAi93lj\nYV3FUfr8miuqiiMz32pVIF6lLeRkWjDGlCl3HA32RYtmU8PaNgQAtG1sTbYqN5cf1bHF8/n8akiL\n0Ntw4Q6wrYmmLVuJnX3JcUI/t37GWFw+NPj4seJCjvZzb5SfhYxUc37rgu5mmrJW19PbeRGbSzDz\nu8zchplbMvPdyrzbmHmWMn2YmScycytm7sPM61Sv/Qczd2TmTsp4yJ6Rl5mGyf18FaKs9FQ8M+8H\nSdSlcvBoXYLcy4dGTmwnvE+dIAcAvvlRsquGIi2IzBXyzE9E85T/9xHRXtXfPiLaa1+I7uaW5oRO\nmNCzHMvvrEabhtaMjxNLHcWsz+P8fk2NbzOOytRQpfIdYs0AQvdxdOJOJDMHVpot2o4bK83hWBdv\n9Ct2S7/IWHk8fKMWENEDRNRS+XsAwDdOB+Vmt4z1DWj/0Xdbcdfb32LhT9I8U09+VrrTIQgb3Htm\nF6dDcC25qGStkBVHZh6k/F+PmfNVf/WYWQZEqeWRMcBJ+RdxuSg3lpdpf59CO/xuVOSssmYIt7+N\nD8dhzkFYXhh9PwCrsqq6/oJMlG+6MsasbpE2Y9d+srIyF7xql3/25rgGwFEAr8A3HuNhAFc5GpHL\nZaWnonNZgf+xFA6DVXVohNwEPSeLQGf2LMf6GWOdDsOV1HfghfkMtTUiokIi6kJEPWr/rA7MK4Kb\nWUVf6ImU0MWs06Oxwp+ByqVthVXvFCDP7FHun+7VrNDw68K9w2je/bhuTdC8JDeKV+hsL9QGSb0M\n2VKu1za3bBHne3NarMdyLK+y4uOx9CO3MPGOWzHzAWaepgx70ZuZf8/MByK/MrmNUiUF+e93Wx2M\nxD2++6WuAdiRmhNhlhSJ7p9f/eR0CI5iZsxfsx3PfrbeP29oWxlSyGwRK45EdBeApQAeBfBn5e9+\ni+PyDG0ZJ5bmg1ZlJVUzfPcqluahCTaOYyzum1DXbCSapD5GPnsjR9TD53TH3BuGGt6uHiMXBLRX\n+a26iBBU0XLZweCGPo6hnrfkDrCBdT56bl0CgjO6a8elD7Nuzd502UdtCSKaQ0T1VY8LiWi2kzF5\nwWRV4pcnPl7rYCTu8boqe2TPKC5aisRw75md/dMz3kvuxFFvLdqESX/7En+avco/TxLjmM9Im4az\nALRkZhkoRocZd8UircHOgpSbCm2xZFU1q9Ac7eeaqq4ARvFSN+3vVIMXMNQFfdvuDiVBq7TelYX4\nen1ggoNYKua+15i9wyLHMahViX+6Y5N8vLVIO8SgwS0lwy1HoETJpArAl/WUiMJ1eBaQ/nt6dhw4\nivLCbDx/cR80K/Z2ywwRvYqium4Q6alJ8dsZ0oadh4LmWZGwMtkZ2aPLAVgzUF0C8MrXlGAs1lgK\nbZFeE2s/SLvKj/XM7hMSRZk93HsszstEvcw03KokhYjGye1jK4OunzEWk/oaTwxkF7fVG6Md19DI\noVycm4kPrz8p/vU69KNk2kUbc1bjdieIyP9FI6JmcN9hLlxuz8FjeGvRJuw/UoMWDfIMX/wTiUN9\nMSXZP3/1OeiULqW4dWx7ZKUnzhBtbmGk4vhHAIuIaDYRzar9szowrzBnqAz39Bm0IpJYkxiEi6V+\nTjpO7doktoA0btMZ0N6+poj6W3p8Ug8UZKdj2Z3VGNctfLM/vebRj00ypxty1/KCoHnqkPUuSAxu\nXYJ4ubylqm2c7OP41pUD6tZp4weQmZ4UV4hvATCPiF4gon8A+ATAzQ7H5DmHkjgJxvT/fIuu0z8A\nAOw+eMzhaIRT1HfUtuw94mAkzlPXm9uX5uPSwaGHOxOxM3KGfg7AvQBmoK6P45+tDMpLgvvnRF/C\niuUi0VXDWqIwJ4pmO3H0cVxye1XEZawQrqK7+LaqgD5VVm1HT1pKhK9NHE1Vz+3TFO/8ehDGdC41\n8NrQG8pMi/4qm97aXrysX/gEPjpPpps0vl/ghsxfZTys69upfexcc2311WtDrRVi3Ce1Fz7aNa6H\nP47vjHYWjQvrMrMB3AqgHXxZVQcDkIHYDFCf9371j+QdwWTm/B/80xeo+n6K5NKwXlbA4yUbkneY\nGnVSnGS+qGQ1IyW8g8z8CDPPZeb/1f5ZHplHaAtpsY2tF77ApVd4vLG6HRbdVqWzdKxbCZab4at8\nFGSHr6BG7KNp4+0Ks7YULuTGBVmhn4x6O4Ebys9OQ8cmwXf5nJKqsyPUc6z6bN3cza1bRX2M6dw4\nYJ5ZFclo16K3/62o1Br6nGPdrPKTWZiTgXP7uK+ZtEUeB9AXQB4zvw1gH4DHnA3JG/591SD/9Cff\nb3MwEve4eGBzp0MQDinISUfXirreZBt2HXQwGmdt31+XiiXJW+1aykjnrk+J6I8AZgHw3wdn5oWW\nRSUCmDFeldHCpHq5eb8bjgNHa+LedsxNVe26sxnFsr0rzc1aV5ybYer6zEbkuwv61Cfr7N2u9lNx\nSe+v/Kw0/OuqgbZtz8ksxwFJkGzcbpLoy8w9iGgR4E+O4+4fA5fIkGQXQfIjXNwViS1b1bz/hEvO\nlXbTljOnntTSoUgSn5Ff4O4A+gG4BzIcRxBz+jgGz3tlar+YkqKE31B0sRTmZqC8MLaBy81gX5NY\nnXk2tI3MSE0JyIhmxnabFefghqo2sb1Y2bR2f/xuVLvAxSJ8MDIwtw7DTcV9C6YQML57GUZ2aBT1\namM5gtx2p69d43pOh2C1Y0SUCuWSCBE1ACCD8BkgWRKD1YswFrRIbH8cXzcc2IkkrTkeOx74vmNN\nyigii/gLzMzDdP6G2xGcF5hRwdC7pd63RbE/tbYpQ36Q/d3EXp7az/CyZlXerG4Wa2aFsk3jvLhe\n36zYV+lU92f8343DcPXw1nGtVyvSWJO2jBvo8btZRk5i6vfcOD8LD5zdLba+qjF8IHrdUgOSIFnf\nUjXAv68eiBV3Vse4Rk94BMBbABoS0d0A5sF3cVZEkK6pOM5bvd2hSJyzZuv+gMeW9CsXntG8pG4Y\nluteWexgJM45UiN9Gu0S8deGiBoR0TNE9J7yuAMRXWJ9aN4QnP0xhsqOi0rFhgqIBgumrRvm+Zdf\nffdo3DIm/B3Uu8/oFHUsZrBrEHWz/eW8HvjrBb1M7XMJRD4eA/s46jxvxc7z+EXUVg3z8OT5PdC+\nNHzil6D6cqThOHQvuJivyECT6ng/d/XLM9NSkZvAV4yZ+UUAN8GXtXwzgNOZ+TVno/KGbE16/dVb\n9zkUiXP+8cWP/umezcztPiG8LxnvOn65bqfTISQNI5epnoUvA1zt2AffA7jOqoC8RltUiiU5jh2V\nFILB4ThMDKZ2XcyM9NSUiO9zUt/AzHC2DYlh8pbiWVs0u78gOz1iU8Z4RR5H0AM1bBsY+dxG8e3y\nPAAAIABJREFUdSpFXmZ0dxC1q22cb+AigckfSWZaCsrqZ1u22cIcX6W09kJTsmDm75j5MWb+CzOv\ndDoer0hNIZzerW4opmfm/YAzHp/vYET2S1O1AhnYstjBSIQbbduffMNyXPr8Av+0NnmdMJeRimMJ\nM78Kpf8FM9cAkHvCClOakUY530oxbTPEi+KNP8W28S3NXZ/RSwde7Qpo14WO8DMSU7h9+9UtIzD3\nhqGBy5u0Ywa1ahDyuZYN8ixtAt6hST5euqwvbhkbPJ6qEHpqVHdUNu46hEU/JdcQBH+bVzcUxyWD\nZKw6Abx2eX9cMsiXXXfFz3scjsZZj0/q6XQICc1IxfEAERWjrhN/PwDJfVSqmFGcClUoi6ZeMfPC\nXhG38esRkfu+WdNfLbaVdmiSFOO5uU5A3zadIzwg2yb55ljO5ZVsq/aA+rPIz0pHdkbkO5axxDKq\nk/ErtLedol/Bi+e3Y0DLEsmWKQwbqzPObbIk5dL25SqIZjxnkbB6Vxbhxuq2yExLwbzVO7B172Gn\nQxIJysiZ+nr4huJoSUTzATwP4BpLo/IQR+6+6BjeLnKTxdO6Nom4jJlq903sw3EQOkZZeTRtKIKI\nfctCLxBXU9U4XusEq+K1KtmDVWM7RTrCYxkOx+ALgmc51EE38IKC145k4SWjO5di3T1jApJO1SRJ\nv65t+5KvGaIwJis9Fb0rizBz/g/oc89/8f7yX5wOyRbqMqY2V4Ywn5GsqgsBnARgAIBfAejIzEut\nDswrtAUkM/t8mVn0MrouI/EHJ7zUf028+yKWV8dSR42qkGtBeXjFndX4zckxDqFhItL8D4SoQGvm\nGen/Fq28rDS8eGlfVJnch/N/Nw4zdX1GGen7bFZmYacGPpa6orBTSgph/5G6cYanqvo4JbKtUnEU\nYbRVDWW0aMMuByOxj7qCnObUCTCJGMmqOhFANjOvAHA6gFeIqIflkXlI/xaJ0zndSOFveLuGBlcW\nxXaNL2o6/bpRfBH9dqTximBuZprucAheQAR0Li8ImGdWk7GBrUpQYPLA1tpxM41YcWc1PpsWfgQi\ns45f7fevaXH08Xrpbt95fd01fqTwlpcu7eufnrtqm4ORWI+ZceBIDTbvrmuC+Nh5UhQTgdQZsJOl\nEvXqgg3+abPLDCKYkeLqH5h5HxENAjACwDMAnrA2LG8pjXM4hFDlPCca3kT6mVlyWxUu6N8swlLO\niaXM3LmsIPJCtes3uNxZvSuiDwTuu2ujW6lWzaxNYJRroO+d6YHYJDczDTlWvz+N2srfxJ7lqnl6\ny5m/7b9f2Dvs82Zus9RIllghQmiZRJl4n/5kHTrePhvPfbbeP29sl+C+nkLUSk3x6BXpKO04cBSD\nW5fgyfN7oLqjZFS1mpGjqrYn9lgAf2XmdwBEHNSLiCqIaC4RfUtEK4joWp1liIgeIaI1RLRUfSeT\niKYQ0Wrlb4pq/t1EtIGI9mvWlUlEryjr+pKIKg28N1ewpZ+k8baqYRXkpAc3zw2VVdXE9zW4dQmq\nO1oz9ERhTgYqinzNLc/v57sDEm9l3mX1P0NqP9do7lhph1CJ9vWG6ez4f17WD09MsueKe8SxLS36\nEsey3ngjGRahRYH/hrJmQ2676CESX72sxB3rU+vtpZsBAF+t941X93+nS18uESxTlWTsiY/XOBiJ\nPVZu3oulG/cgIzUFozqVeqrFjVcZqThuIqKnAJwN4F0iyjT4uhoAv2XmDgD6AbiKiLTp+EYDaK38\nTYVyJ5OIigDcDqAvgD4Abiei2lFu/6PM07oEwC5mbgXgQQD3GojR1ZwZjsOdfTR/N6odnpocPnNs\nzFSBmpbaPMo374aEgHpNTPV+hNVzjGT5tEr/lsUY1LrE1m06+THpZ7jVWc6iH44LB1Ras2IhYpSd\n7tzvj93UrQ7vOr0Tzu/n3pY/wjmTVS3Cjh3nhM42PO4v8zD64U8BAPPXbnc4muRhpAJ4FoDZAKqZ\neTeAIgA3RnoRM29WEuuAmfcBWAmgTLPYOADPs88XAOoTUSmAagBzmHknM+8CMAfAKGVdXzDzZp1N\njgPwnDL9OoARZNelBw9c4DCc1dFF74Uo+gpVrPFrtxPP2Jrx7EMzK+4JRbNb/nqB7yJCol1djDcT\nsbIWU2IJufZQLQxi2G7iFmmEHbTf/7cWbXQoEuulqGqO0o9LhJKZFngx5djxxP2VXbKxbmRAGQfY\nPkayqh5k5jeZebXyeDMzfxDNRpRmo90BfKl5qgzABtXjjcq8UPPD8b+GmWvgG2syKGsNEU0logVE\ntGDbNnd0pndTZSGWSKway/KULvYMHxI4bmH8mI19pm67EBjp7mK45Syn2Vd9mhf5YrE5DCu3F7k5\nrN48nc/MqiFHIhywCVaHFx5xcvu6ptVzv3PHOd0Ki37a7Z/OT6ImuiI+y39OjmHXz5dEa7axvOcs\nEeUBeAPAdcy81+rtRcLMTzNzL2bu1aBBA1PXff/Erqauz8+UmpkJ64h2k3GUJNfPGIuRJg/FEA1t\n7JZ9tgqX1SFdI/J4mvbEEa9ohyzx4p1U70XsPUQ0iohWKX35p+k8r9vXn4gqiegQES1W/p60O3ar\nDGpV11z9hNuuxplkw86DAY8zvJqGW9hi1tUD/dPjH//MwUjs48VzpldZ+utDROnwVRpfZOY3dRbZ\nBECdfrJcmRdqfjj+1xBRGoACADtii9w8X9w8IuIykY73WE4Sg2Ps+2Xml8+JtvVxjx2p8/LJ/Zqh\nSf2skM/rrSPa3dhYyS7ZKM4MvdG6ZnirsM9r38f9E7vaWkEIdQg5dY4IdUSHCufaEa3x8Dnd8OsR\nrePetoEhNS0V7W/Dw+d0syiS5EREqQAegy83QAcA5+rkDQjX138tM3dT/i63JWgbDG1bd8cxUSuO\n6vEqASBfmqqKMLqU13c6BJHALKs4Kv0LnwGwkpkfCLHYLAAXKNlV+wHYo/RfnA2giogKlaQ4Vcq8\ncGYBqM2+OgHAR2xzzUVvc2mpsRftatdWVj8bf57YFd0qjP0Y5GWmoZNmiAmjZb6Ymqq66EJPLLFE\nKpA3qJdp6JbgfWd2iX7jiom9yvHU5J6Y1Mfe5ha/rWrrnzay67o3LYy8UJTO6lUeeSED/eoW/mGk\nSRGZ5zcj22BctzKkKxd/wl7YcNH3SE/kpqqBbyA1ScYQs1EfAGuYeR0zHwXwMnx9+9Wc6+vvkMqS\nXP/0u8t+Qc3xEw5GYw3tJ6g9vwuRbN5bppfqRNjByjuOAwFMBjBc1TxmDBFdTkS1VzvfBbAOwBoA\nfwVwJQAw804AdwH4WvmbrswDEd1HRBsB5BDRRiK6Q1nXMwCKiWgNgOsBBDXjsUq4wqBZZ+wze5aj\nODfiKCgR2TVG298v6u3JpgPx3LEsL6xrjhjtWogI1R0bByRAcEpAn08b+tD96qSWhmIJmK/zfJEJ\n349I4n3rHEejZDv7M4aNw8B854/ihGSk73+4vv7NiWgREf2PiAaH2ogb8wBEY/Oew06HYLoaVYKT\nkjzrf+dEYknEzKp3vf2tfzpSyylhLst6WDPzPEQoPyh3BK8K8dxMADN15t8E4Cad+YcBTIwpWJfS\n7rzaQuL47pHyBOm8Vvk/0hhtZiXqGda2YVDzGjsY+X3MSEvB0ZrwV6XjLZB7sdJshtO7NcHew8Y/\n96cm90TLBgYG8U68855jOjbJx4qfo+9uHukjSNJD3is2A2jKzDuIqCeAfxFRR728A8z8NICnAaBX\nr16e++Zt2HUQFUU5TodhqsPHjvunk/XcImL35sJN6N+yGE2i7GvvZuqL7FUdGjsYSfKRHtYWM/Ij\nH+1pYHTn0ojLxHy2j6mpZ/wnMrNOhZF2d6P8TFyuucNl2onY4+fzpyb3NLRcuM/71lPMTYmdn+Xr\ny5OlGa8tUQtP4S58GH3HkZbrXVlkNBzhLkb6/uv29WfmI8y8AwCY+RsAawG0sTximzw1uSfaNa4H\nAFj1yz6HozHfhCc/90/nZUpGVRHZI+d290//9rUlGDDjI2zbd8TBiMzDzNi465D/cWObc0MkO6k4\nWswtxVujBW2vl8drw79yaEsMbRucNZdgbHDIUPsrbCWZ1ctF5rZL+S2UO3+xXAiwqiL326q2uHVs\ne5zaVX9oFtccrw7FEemzmnvD0KB5qSmENo0M3OUN2laE513zYSSsrwG0JqLmRJQB4Bz4+var6fb1\nJ6IGSnIdEFELAK3h6yaSEKo7Nsb71w1BZXEOPl29Hat+2Rdwly6RROpuIkQoew4ddToEU3z5w07/\ndJfyAl8eCmEbqThazFgWTvMLXLGu0cxI0pSmBANbxZbhNRa1lbHGBVl49qI+yEqPfIibPXalm8bl\ntEK4w9Xsd56dkYpLB7dAqK6f4fZ1pD7B0WQezlSOo/E9DCTyCSOeY8PIz0Rtdt5azVWJQ2pde3Js\nGV6tuNCRgF1vLKP0WbwavkRxKwG8yswriGg6EZ2mLBaqr/8QAEuJaDF8SXMur80bkEhOatMAH323\nFdUPfYJb3lrudDimOHEi8EtSqfOdFkJrVMfGuHBAJRrl11WqEuXi3u6DdRVgGZrGftLmwQZZ6Sk4\nfCz6TG8tG/ruCoxoH/94hvH+XJTVz8aQNvrjXob6LcpKT8XcG4ai1GXNCGIpq0bzGga7505YDMyO\n/dax7fF/76w0dZ2RQjyzRzkWbdiFHQdCX2G97uQ2+HT1dkPby0pPxfI7q5GtaTJrJyPJccZ0LsXb\nS38GANx9Rifd9eRnpcdVYTNa+JBKofmY+V34ksqp592mmtbt68/Mb8A3NFZCa1ea75/+5sfEqBff\n8NoSp0MQHpSRloI7TuuI95bXZR89fiIxfpT3qfIoHEvALMpuJ1V1ixFIt7mYWn6Wfv29eUkult9Z\njXP7VATMtzJDVqhC4fxpw/HH8Z11n9P2P3vw7K7+6eYluUHPWyneOo9ZA84n+l3HaFw6uIV/ukt5\nAa4cGjqLai3t3gt11zrWSq5eM+ZI8jLTQg4x4ZZPW5219WQTLjgJ4SXqVgaJMhzMm4siDWEtRGjq\nmxZVD37iYCTxO3CkBht3HQzInDypXzMHI0pOUnE0kW51joC0lNC7uapDI/8Yb3ryMtPCXuEvzEnH\nMAOFYCvHcbxkUPOAAc7HGEje4xS9/RBN5cNInd3rlcZ4oo90N6phvSxcPKh51OvtVFaA9TPGqraj\n/B/1moD1M8bi2Yv6xPDK+IUbjsPOoyYxrjsLEag4r65ZXrjzrld9etMwp0MQHjMoii4ZbnfO019g\n0L1z8fzn6/3zzupVEXJ5YY3E+2V1G4tLaItuq8LfTSwEx3IHJys9FdePdEeCvksHt0CbRnkYG6by\nGssN277Ni3Bun6b408Quxl7g4bpjpAxlcQ9VEt/LA9fl5TbBFjO7YYI0PRVup844umrLPuw5eMzB\naOI3c94P/ulx3ZoEjBUshBF/mmCwzOIByzbtAQBs3+/rgtK/RXG4xYVFpOJoA3XZdt7vnLliaPQu\nmDlDa4RozmdDIb9ZcQ4++M1JAVee1X41pIXu/EjSUlPwx/GdUV4YenywRClXXzQw+juCteyqxnn1\nrq7Zccf6lYqnubvU1YVbabMFf7s5+vFK3WLr3sOYrhrk/KGzu8mFMhG1nIzETGXy8tR++OfUfk6H\nkZSk4mgxbdM0bcWDAfRoVmjqNuNramhaGEGs7JtZt43wz184sHnQZ2JFJcTL5/dIfYP09pdTb9fM\n7Y5o1xATesaXNTUZhPuOXTig0rY4hNAiIrRoUJd11Mv9HO99f5V/enK/ZlJpFEKldoxnYT+pOFos\nUkWGGbipuq3BtVmbyt8ssWwr1Gu0u69FiFTkhTnpmNS3KQa0lKYLVsnJiJzkKNTnqE6YZAZ/H0cT\nj+tLBjfH/RNji9OpQp1uJd5AKNFcwpHyqvCSVNUB++nqbQ5GEp83Fm70TyfqmJTCfurmz15y8GhN\nwOP87MS8k+oFUnE0QbiCFSNSdY+R5qJxaNxeSJx1zSB8+fsRQfMbF2Tj7jM6G9qXepV5M2+GEhmr\n4ttxB9YsK6ePwsI/jAQQ2zESqrlM/FlwXX7AaoRLjuN2ncp8Qx20bJAXYUkhnHPTqHb+6Uc/WuNg\nJObJNnDRTggj/jR7VeSFXOjv89cHPM7PljuOTpEquzCdlUX5vMy0gAQIZrCi7uG1Ck0kRgsuIfu3\nRrm9aKpX7UvzcflJsfVdNUusn/bTk3ti6gvf+NZh8Jj56pYR/v1s52F2Vq8K9GhaiNaN6kVc1mhc\nXq5IC3ca2SHxhqG50XCrJCGCvXXlAJzx+GcAvDuWo3a8xrwE7bvpBbLnLcbMYQuEbrvp5NWkI5dG\nMcSDy3a558RyhFi5z9+7drCFazfG6PvTfr+qOjaOelsN64XOessMZKYZqOTH8IEQkaFKY20cQojY\naAv39aQ/l4hDg3p1yQKPaipgXqE+d5bkZSDFw/2Xvc49bSQT0KS+TVGkGpBYz2+r7LmSaG8fR2Mb\ne+OKAaZt88wQSU2sLsA+dHY3/3TPZoUY3q4hpo/r6NHqtwlseuMJdkPXVM9d3AfXndwajfL1MwsL\nkUw27T7kdAhRO1rjzcK9cKdEuJCnrifGcsFVmEcqjmbSfDl/P6Z92ErU+hlj0aFJfrybMZWdBfKW\nDfQT3XhVVnoqZl7YG60aGrsr41WJ1gzXDFP6N3Nku3qfRPOSXFx3cpvwLR2sCwmAVOyFe5zx2Hw8\n/OFqHKnxToKZ1Vv3OR2CSCDa7j0fr9rqUCSx23ekLjmO3Gx0llQcTRDpGHbDMW58HEczthVifhKV\nJpPorcYl1nT5lgyhEsc6J/evxPoZY/2PHzq7Gz68/iQzwvKsBLjILTxq/rTh/umt+47gwQ+/xwuf\n/+hgRNE57S/z/dO9K80drkskn8LcDFQW1w0F94wHM6s+/ck6/7RXu1QlCqk4WsitlYe3rxmEly7t\nGzDv1K5N8PzFfSzN8KrNIqr+8tv5Q6BttuHWz8mtYtldoV4ztG1DAL7mlWZz8mM9vXsZWjWMPfuo\n0WMy1osxXsroK0S08rOC0zcc8Wjzz5en9nc6BJEAzu3T1D+988BRByOJ3t7DxwIeTx3ibDK8ZCcV\nxyShLl92KivAgFYlAc8X52ZgSJsGNkflXZHK626/Ivboud3x9jWDbNse60wX5Wb47ziepDn2It7F\nj7T/3b37dXkwZCFcKScjDS1KcpGZVlfE8cpvgnq8ujtP6xhzqwwh1KYOaYGvlKHMVvy8F+Mem4/t\n+484HJUx320ObLpdUZQTYklhB6k4WsjulPk9mta3Z0MRGH6/UewXr5z0a7k93lO7NkGnsoKYXqv7\n3lz+fhNVrLvdzPuN/7txqIlrEyJ+qSmEj24YivP71fU9fvjD1Q5GZNwj/60be9JL/TKFuxERGubX\nZeResmE3Rj30qYMRGcPMOOupz50OQ6hIxdFDagvsoVqZPas092MEF+7dWK53qnKlHTvO7XcHEwGF\nmDZjfaEkW2PM8T3KcO2I1rZvtzivLntrtE1gpcWssJJ67DevNFVdv/2Af7rGo2PuCfeqUo1z6oU7\njtrvbb8WRQ5FImpJxdFEQRUSpXTrRMXEyQJZqH5XiZQcJ5kLvLF8jkm8uwwjAhrlZ6GsfjbuPK2j\n4dfUeuCsbiiMMPxPLSPHr9FjPHG+1SLRVHXwTtp+ZgYz45e9h9GsOAeT+jbFlP6VToclEszvRrdz\nOoSoLN+0J+DxC5f0DbGksEtwD3LheW4tyAUnx3GIJMdxVG2fnUhjnIaTSBch1DLSUgIyQgohYjeo\ndWBf/uMn2LV9Bkc//Cm++8XXl6tpUQ7uPqOzwxGJRJSb4Z1i/7HjJzDhybpmqnee1hHpFiZwFMZ4\n5whKMBlp0R/80ZzugpqqurCg7caYjKqfkw7AN3ajHg+/tfiEuEul3h31czIwY3xnnNTWXcmYvPiZ\nxZxV1cR7wF7cbyI5Ha05gewM/d9sp9VWGgHgp50HHYxEJLLcTHce/3que2Wxf/rG6raYMqDSuWCE\nn1Qc7aApWM35zRAUKBWPZBK2kGtj4dOMIvO00e3QrDg3oL+AnhQCHjm3O65+aZEJW/WeB87qCiB4\nn5+jSg0eCyOHi9Rn7Ofli0Ei8X29fqcnsoffUNXG6RBEgsrR3HFcs3UfWjWs51A0oe05eAzvLN3s\nf9yhNN/BaISa3PN1QOtG9dCwXlbkBU3kxuKcW2KKJY6cjDRcMqg5UiI0e0pNIZzSpUlsgSWA8T3K\nnQ4hpDSXNlmzg5l9dNV9uGV8SOFmF8z8yukQdKmT+ADAlUNbORSJSHTaptovfvmTQ5GEd0A1LA0g\nGYbdRCqOFrL74nt6iu/j7FgmV2bCScbCbX5WGs7pXeF0GKYx47vVpbwAN1a3RYuS3PhXloQifQS1\n37OT2zfEv64aaH1AQugY6rIm8XpmvPddwONIFySFiMdpXesuZu8/XBNmSecc12QUbtPIfXdFk5VU\nHG1gVgVyYi9fwb9zuf74e9kZqXjjiv54+oJelsUQj6DkOC6IyWrhMuraWX8tLcjGjDO72LdBDyAi\nXDWsFUrqZUZeOMHYeexlpaeiW4U7xpgVyeeZKb2dDiGir37Y6XQIIoncUNXWP/3aNxtdeTFdPQxH\nVnoKWjTIczAaoSYVRxM0UgZVLcgO7Ldo9jAcIzs0wvoZY1FWPzvkMj2bFSE/Kx3dKwpN3bYVnBo/\nUfsbmQz9ssxMhhILs/bwG1cMwJT+zSIvKGyj9/VJ/G+UM4hoFBGtIqI1RDRN5/lMInpFef5LIqrU\nPN+UiPYT0Q12xew0t2ZRVVN/h0ZG6DcvRLyaFucEPN57yH13HdVNU9NSpKriJpZ9GkRUQURziehb\nIlpBRNfqLENE9IhykltKRD1Uz00hotXK3xTV/J5EtEx5zSOklPqJ6A4i2kREi5W/MVa9N61fj2iN\nB87qiuqO+mNGOXHaOrlDI3x1y4i6GFxQOQoXQ+0zZ/eqwINnd7UnIAtlpKagT2URHpvkO6SHt2vo\ncEQ+FUWhLzpYyaxqa89mhbhzXCdDx7P7rqHqq5fpS1bg/DfUfl75jNyAiFIBPAZgNIAOAM4log6a\nxS4BsIuZWwF4EMC9mucfAPCe1bG62fdb9kVeyGbq37Pz4kwcJkS0tu477HQIQdZtO+CfduMd0WRm\nZTW+BsBvmbkDgH4ArtI5yY0G0Fr5mwrgCQAgoiIAtwPoC6APgNuJqPYW2hMALlO9bpRqfQ8yczfl\n711r3lawjLQUjO9RHlSYdbquZncCnmjp7Z97J3TBGd3rEqpY8YOhXaMVHxMR4dXL+/uvHj81uSeW\n3VFlwZai8/61Q7Dg1pNNX29hjjImI7mjAnRKl1LjCzt8Tpp+ekdnAzBR/5a+cfM6lwc3TX3u4j64\nT5pLx6MPgDXMvI6ZjwJ4GcA4zTLjADynTL8OYITq4urpAH4AsMKmeF3p1n8tdzqEIAE3Rd3wAyqS\nypa9R7B80x4crTkReWGbXPPPukz0Um10F8sqjsy8mZkXKtP7AKwEUKZZbByA59nnCwD1iagUQDWA\nOcy8k5l3AZgDYJTyXD4zf8G+GsXzAE636j0IYZb01BTUy3J+CJbczDSU5Jnfp+/FS/ti+riOKMhO\n1/2Rt7ssNLRtQyl/RenuMzphTGf9VhNGjerUGEvvqELPZnVN5WsvqJ3UpgGqO8W3/iRXBmCD6vFG\nBJ9T/cswcw2APQCKiSgPwO8A3GlDnK5zx6l116y37zviqjsYzIxFP+32P05x+oqzSAovXtrXP/3c\n5+txyqPzcPc73zoXUBj/vKyf0yEIFVsaDiv9LLoD+FLzVKgTYbj5G3Xm17paafI6U3WHUhvLVCJa\nQEQLtm3bFsO7Ma72598NzUST1W2ndohq2IVk/Kg+vWkY5k8bHtc6KopycEH/ypDP21lM6xIieZRb\nuaUMO6lvMzw+qWfQ/Gi/E/lhLpAk4/fLJe6Ar0XO/kgL2nmOtMupXZugU1k+6uekY932A2h+87t4\n+St3DEOwUFVpBIDG+e5uKSQSw8BWJVhz92iM7NAIc77dAgBYsnGPw1H5qJvOntO7Al0luZqrWF5x\nVK50vgHgOmbea+GmngDQEkA3AJsB/FlvIWZ+mpl7MXOvBg3cn6Y7kREB43toL5iba1LfZlhzT2B3\n1+DkOJaG4HoVRTlhEy55ybfTq/H65QOcDkMY1K95kdMheMkmAOoxdcqVebrLEFEagAIAO+Dr9nEf\nEa0HcB2A3xPR1XobScRzZHFeJt6+ZjD2HDrmnzftzWUORuSzZe9hnPnEZwHz2jaWYQeEPdJSU9Cv\nRbH/8Y87DoRZ2j6XPLvAP33ZkBYORiL0WFpxJKJ0+CqNLzLzmzqLhDoRhptfrjMfzLyFmY8z8wkA\nf4WvP4ij5E5jZHef3tnpEAAA9bJ8CUq0zYSSeZB4M9m1F3My0pCRJhnYInFDU72V00dhQKsSp8Pw\nkq8BtCai5kSUAeAcALM0y8wCUJtMbgKAj5SuIIOZuZKZKwE8BOAeZv6LXYG7xRM6d9SdcvjYcSz8\ncVfAvES5gCe8Q32He9fBY9iy1/lEOcs21d35dMGpSmikWbVipUP+MwBWMvMDIRabBV/z0pfhuyK6\nh5k3E9FsAPeomptWAbiZmXcS0V4i6gdfs9cLADyqbK+UmTcry58BwDU94JOl6vGHUzrgLx+tjum1\nkerYZg7doTc0xTMX9sZ7yzajierE/fcLe6NVQ2vGDqoszsFtp2pzRYlkcf/ErmjTyHdsFeb6EguV\nFVpTaJzYsxyfrd2BTbsPWbL+aKi/xdkZqY7F4UXMXKPcJZwNIBXATGZeQUTTASxg5lnwnXNfIKI1\nAHbCV7kUiuqO7hnq4op/fIO5qwKbAkeV1EsIE+RlBVYD9hw65h9izm73vv8dKjVDhRQr50fhHpZV\nHAEMBDAZwDIiWqzM+z2ApgDAzE8CeBfAGABrABwEcJHy3E4iugu+K6wAMJ2Za0fIvRJDDyeLAAAa\nQklEQVTAswCy4UsrXpta/D4i6gZfd6r1AH5l1RszKlkqjLUuGdQclwxq7nQYUSMilNXPxqWDA5tE\nDLNwCI2Pbxxm2bqF+5UXZqOLknl0aJsGePL8HhjRXr9Q26phHtZsjdg1LaQ/TfQNb1M57R3/PLmI\n601KtvB3NfNuU00fBjAxwjrusCQ4D3BTKyBtpREAbhrVzoFIRDLr2CQ/4LGT35AnPl7rn75oYCVu\nGdMeaanSgshtLKs4MvM8RDgGlcyoV4V4biaAmTrzFwDopDN/cmyRWs9F5yrPsaLZoTR9MEdRgl0J\ntPNrGpB9nwijOoW+0/DmlQOwY/9R64PSqJ/jS3LTRWdoDQBoX5qPxRt2h/x9G9WpMc7uVYHfVrfx\nz3NTwV2InQeOuuZ3bHDrEqRK1whhs5K8TORlpmH/kRoAwHGXFJBqjrNUGl3KyjuOyUn1uy9lpPg9\ndl4PvPjlT2hfKgkD3OTZi3qHTOLglsO+eUku1m5zR2d/rWhOzflZ6WEzlca0fQMBVBTl4O1rBqFN\nI/3P+bmL+mDlL3uRmabf5DQzLRX3TpBxG4V7XfXiQvxzqv2p/muOB46XN/u6IZIURzhGfb3CqbEc\n1ZlUAeC6k1s7EoeITKrzZnPHxRpD1s8Yiz+c4lw/u1QDNeuKohxMG91O7lS4zNC2DVFaoN8nzy1f\ngQfP7oaZF/ZyOgxP61RWEPKuf0FOekBGPiPkWyzcZMOug45s91+Lfw54nJ0u/X2Fc1JUNccL//51\nmCWtM/aReQGPiy0Yb1qYQyqONjAzsYvZarMr2lkvm9S3KT6/ebg0Q0hwvx9T11/Hicpkvax0DG/n\nnmQYbqKXIEqIZPDU5LrMqht3HcLug/Y3Az987HjA46wMORcK56gv4u88YP/3AQC27TviyHZF9OTX\nyiITepbLXbIQ0lNT/HernCjAumEogmTQosSajLRudvHA5jg5RJIbN7nj1I4oybO/b5f8JAqnDdQM\nAePEoOfavoxZcsdROKiqY2OnQxAeIhVHYTspPNpjyW1VeOOKATG99rLBzVGQbV6/umg/cq9edLnt\n1A7425TIzWOdfnejO5diwa0jHY5CCPvlZaahV7NC/2Ntf0M7aHPg5EjFUTho+riOAY9/3OFsboA3\nr4yt3CLsIRVHs+mUCD1aBraM2Tf8MmXAd10FOenIy4wt/9UtYztgye1VMb1W73A38pE/em533DSq\nbUzb9Jpkveft5mb7InkcU1UW1+84iEl/+wKvfr3Btu2nqAoFYzuXSrcN4ah0zfF30p8+xgUzv8Le\nw8cciadpUU7khYRj5Ncqwbkl1XgksRYol91RhYV/iO7OyW9GtsGEnuUxbU9EFmul6NSuTdCnssjU\nWIyQ/n5CJJfGBXUDnN8/exXmr9mBm95Yatv2X/9mY90DuZYiXOiT77fhP0t+jrygCXZp+lVKsih3\nk4pjgnv7mkF4xkDTOadkpaWiS3kBHjqnW0yvr5eVjtwo76rVz8nA/cqg6MIeRstGLRv4+kVePLC5\ndcGE4NXmsV4iu1i4wX1ndvX3MzykSlRz4oT1F5H2H6nBlz/s9D+Wr4Rwg0fP7R407/O1Oyzf7vET\njDOf/CxgnvT5dTepOJrNZTcvmtTPxggXJ+tISSHMunoQqqVzdsKIpyBUmJuB9TPGYmyXUtPiEUII\ntYKcdKy9Zwz+fmHvgPmvLrC+uWqn22cHPD69W5nl2xQikp6qfr+13l662dJt7jpwFH/+YBXWacZb\n1iaPEu4SWwcoIeIgdx0Sm8uunQghhK6KosCxaK0eEkCb0fvTm4ahQvpzCRcoLcjCr4e3whfrduKr\n9Tsjv8AEve7+EMc1d/mdaG0koiN3HM0myXEiktEwRNKS3wIhXKOsfmClzcpz9TtLN+Pxj9cGzJNK\no3ALIsL1VW3RSNX/12raSiMA5GRIM1W3kzuOAoBkOxTJq3bYkfRUG74DctFECNfI1hRSrernfPjY\ncVz10kL/40b5mbqFZiGcdu2IVgFJcZjZku/FYVXfYrWrhrUyfVvCXFJxFEIktT9N6Io3F21Ct4r6\nToeS8KT1hXCb8T3K8ObCTQCA77fss2QbRzVjRU4b3Q5ndJfM3sJ9WjWsF/D42HFGRpr5P9x/fHel\nzrbzgi7mCPeRpqpm07mIaORu3uDWJSirnx1xOS+TMmNyiHUcR6cU5mbgkkHNJauqEEkoVfW9//di\na4YfOHLsRNjHQriV+k65mdZqEuLMnzYcb105wJJtCXPJHUeXeOGSvgCAymnv2Lrd2kxag1oXW74t\nN1ceEp2d/UrlcxahSJN44TbaDI7Hjp8IGhA9XnNXbQ143EIZdkgIt5vz7RZL1qu9TpvoN04SiVQc\nzWZCcpwMk09a4XRvWojv7hol4+YIS0l1QWHRjijJy7BmxUIkuBRNxXHvoWPIyUhDVnqKaa0Qbnp9\nacDjPs2LTFmvEHY4cYKDvifxSpEWPp4lTVVd5snze2DO9UNs3aZUGoWZ5HRgr9cv7493fz3Y6TAM\nkbKCcJtLBgWm/5/y96/Q/rb38fLX1o/pKIQbXa1JUPPiVz/hu1/24oftB0K8InoyVKN3ScXRZUZ1\nKkWz4lynw7CUFB6FME+vyiI0zLcvhboQiaSlptno8k17AQA3v7ksaNxFM5zUpoHp6xTCTDdUt0W9\nzLoGics27saohz7FsPs/NmX9J04w5q7a5n/88DndTFmvsIdUHG0g9aRAMo6jSFpJfuzLb6Hwki17\nj8T82kufW4B73l2JQ0cDhx148vye8YYlhOWmDmnhn351wUZT173KouzFwh7Sx1EIYSq31Y3mTxse\nVHgTQohafxzfGTe/uczUdX64cguwEnhF0+RVhhsQXpCbqV89WPXLPrRtXE/3uXCYGa99sxFVHRph\n9MOfxhuecJDccbSBpPkXwjll9bPRqqFLshgm+U+B/BYKNzq3T1Pd+ceOxz9sxp5Dx+JehxB2m9y/\nGe44tQPeuKJ/wPwbX18S0/rmrtqKm15fijE6lcYWJS45PwtD5I6jsI0UGZODfM5CCK8Z26UU2emp\neP2bumZ5R2piqzj+7dN1uvN/P6ZdTOsTwm7pqSm4cGDzoH6+sV7823e4BgDw857DAfPfvmYQOpUV\nxBakcITccRS2cVsTRiGEEAIAHjuvB+6f2DVg3lP/WxvTuv7vnZW68we3lsQ4wlu0FcXaRxt2Hoyq\nC8iyjXt050ul0Xuk4mgDuQMjhACQ9FdP5LcwPkQ0iohWEdEaIpqm83wmEb2iPP8lEVUq8/sQ0WLl\nbwkRnWF37F70mnL38cCRGuzYH3uinFrtS/PjXocQTlq8YTf2HDyGwffNxdQXFhh+3d/m/RDweMb4\nzvjrBb3MDk/YQCqOJhvRriHK6mcHZKQSItkleX1JiLgRUSqAxwCMBtABwLlE1EGz2CUAdjFzKwAP\nArhXmb8cQC9m7gZgFICniEi6qujQji/3w/YDqH7oE/T8vw9jXmdRbgaemNQjzsiEcIY2E/Cabb6s\nqJ+u3h62D+8Nry3Bp6u36T53evcyjOzQyLwghW2k4miy4rxMzJ82HG0aRZ91SgizSS4Sl0nyz0OO\nx7j0AbCGmdcx81EALwMYp1lmHIDnlOnXAYwgImLmg8xco8zPglzLCSlVU3O84h/fYOOuQwCA5Zv0\nm9tFsvAPIzG6c2ncsQnhBO3FlBteW+qfnrXkZ93XHKk5jte/2YjJz3wV9NzQtg2QlS7Zhb1KKo42\nkMKSj+yG5CWfvc/oTo0BABVFOQ5HIjyoDIB6bIeNyjzdZZSK4h4AxQBARH2JaAWAZQAuV1UkhUp5\nYeB3c922A/7pUx6dB8BXgXxtwQZ/1tX/rtyCymnvYNnGPbjqpYX2BSuEDbQXU37YfiDEknWG3/8/\n//Sin3YFPHf8hFy38jKpOArbyE9FcpDPObQLB1RixZ3VKKuf7XQojpDhOJzDzF8yc0cAvQHcTERZ\nessR0VQiWkBEC7Zt029mlsheuqxvwOOjOkNynPLoPNz4+lLMeO87AMC/Fvvuupz6l3l4Z+lm64MU\nwkZGEth8tnY7rn5poT8L66bdh/zPnfH4ZwHL/n5Me3MDFLayrOJIRBVENJeIviWiFUR0rc4yRESP\nKB35lxJRD9VzU4hotfI3RTW/JxEtU17zCCklESIqIqI5yvJziKjQqvcWLSksCSEA329BqIGVhYhg\nE4AK1eNyZZ7uMkofxgIAO9QLMPNKAPsBdNLbCDM/zcy9mLlXgwbJlwW0tCAb/7l6EACgYb3MoOcP\nH6vLJPndL3sjrm/Ob4aYF5wQDmiUn4VG+cHfBQDYuPMgAGDKzK/w9tLNhoawkSRR3mblHccaAL9l\n5g4A+gG4Sqcj/2gArZW/qQCeAHyVQAC3A+gLX7+O21UVwScAXKZ63Shl/jQA/2Xm1gD+qzwWLiLV\n5+Qgn7MQlvgaQGsiak5EGQDOATBLs8wsALUXWicA+IiZWXlNGgAQUTMA7QCstyds78nN9PW/6lZR\nP+i5qgc/CZoX7jevgU7lUwivubFafwzSpz4JHLP0yLHYxj4V3mFZxZGZNzPzQmV6H4CVCO6PMQ7A\n8+zzBYD6RFQKoBrAHGbeycy7AMwBMEp5Lp+Zv2Df/fDnAZyuWldtUoDnVPNdo0/zIqdDEEII4UFK\nn8SrAcyG73z6KjOvIKLpRHSastgzAIqJaA2A61F3AXUQgCVEtBjAWwCuZObt9r4D72jRIA+PndcD\nfz6ra9BzPyl3WABg/podOH6CQyYIAXwDqQvhdRN6lod87tjxEzh23NdEtev0D/DNj7tCLqt3MUZ4\niy1tppSxpLoD+FLzVKjO/uHmb9SZDwCNmLm2c8EvAHTz/BLRVPjubqJp06bRvZE4fPCbIWiSpP2a\nhHOKcjMAACPaS9rrZJeWQpjYK/TJX7gfM78L4F3NvNtU04cBTNR53QsAXrA8wAQytosvC+obV/TH\nmU98HnK5Zz9bH3Y9ORmSPVIkho5N8pFChO+37Atoktr6lvcCljvzic+0L/V768oBlsUn7GF5xZGI\n8gC8AeA6Zo7cIcAEStMc3RwdzPw0gKcBoFevXrbl8ZDhOYQTSvIy8fUtJ/srkCJ5rblnTNC85y7u\ng0NHj+ssLYQAgIrC8BmQ73r726B5eZlp2H+kBl/cPEJyHIiE8c6vBwMALn72a3z03daoX9+jaX35\nPiQASyuORJQOX6XxRWZ+U2eRUJ39NwEYqpn/sTK/XGd5ANhCRKXMvFlp0hr9US0SVrvG9XCCkzPf\np/SxEaGc1Cb5kp8IEY1ok1m1aZSHD35zkkXRCOG8lCgrfykEDGxVgvsnBjf9Ft5jWcVRyXb6DICV\nzPxAiMVmAbiaiF6GLxHOHqXiNxvAPaqEOFUAbmbmnUS0l4j6wdfs9QIAj6rWNQXADOX/f1vyxoQn\nvX+dZLYTQggRnWibmk7q28yiSIRwh2i77c6fNhylBdJVK1FY2Wt7IIDJAIYT0WLlbwwRXU5ElyvL\nvAtgHYA1AP4K4EoAYOadAO6CL4vc1wCmK/OgLPM35TVrAdQ2rp4BYCQRrQZwsvJYuEhBdjoASLPJ\nBJeZ7vtZ0Q4aLIQQXkNEAYntbj9VmxxeiORyQf/KqJaXSmNiseyOIzPPQ4TM/Epm1KtCPDcTwEyd\n+QugM/4UM+8AMCKmYIUtzuhehpoTJzC+hyToSGR3jeuEZsW5GKJqBjm4dQOc1KYBbhkrA/8KIbxl\n5oW90en22QCA7PTAO5AdSvPx7ea69A2SNVIkuoGtStCrWSEWhMie+tJlfXHeX7W5MEWikDzRwjYp\nKYSzezd1TXryejIQuyWK8zLxu1HtAu44Zmek4rmL+6B5Sa6DkQkhRPTyVOeKNM35a9bVAzG8XUN0\nLivAFzePQFepOIok8Mqv+uPG6rZ488oBGNY2sK98gzzJq5DIpOQsktYH1w/BD9sOOB2GEEIIl/v4\nhqHYtPsQDh8LzEKclpqCmRf2digqIZyRmkK4algrAMAT5/dEuz+873+uKDcDd5/RCbNXbMHjk3o4\nFaKwiFQcRdIqLci2ve19agqhQ2m+rdsUwm3GdG7sdAhCRKWyJBeVJblgZmSkpuDo8RMBdyKFSFZZ\n6alYP2Ms5ny7BWkphOK8TEzq20wSRSUo+dUTwkZrdcbSEyKZLLm9SgZFF55FRFh51yhs23fEn/BN\nCAGM7NDI6RCEDaTiKIQQwjZS2BZel5pCaFyQ5XQYQghhO3dkKRFCCCGEEEII4VpScRRCCCGEEEII\nEZZUHIUQQgghhBBChCUVRyGEEEIIIYQQYUnFUQghhBBCCCFEWFJxFEIIIYQQQggRllQchRBCCCGE\nEEKEJRVHIYQQQgghhBBhScVRCCGEEEIIIURYUnEUQgghhBBCCBEWMbPTMTiGiLYB+NHpOAwqAbDd\n6SBcRvZJMNkngWR/BEvmfdKMmRs4HYRXeOgcmczHdCiyT4LJPgkm+yRYsu4TQ+fHpK44egkRLWDm\nXk7H4SayT4LJPgkk+yOY7BORaOSYDib7JJjsk2CyT4LJPglPmqoKIYQQQgghhAhLKo5CCCGEEEII\nIcKSiqN3PO10AC4k+ySY7JNAsj+CyT4RiUaO6WCyT4LJPgkm+ySY7JMwpI+jEEIIIYQQQoiw5I6j\nEEIIIYQQQoiwpOIohBBCCCGEECIsqTjaiIhmEtFWIlqumldERHOIaLXyf6Eyn4joESJaQ0RLiaiH\n6jVTlOVXE9EU1fyeRLRMec0jRET2vsPoEVEFEc0lom+JaAURXavMT9r9QkRZRPQVES1R9smdyvzm\nRPSl8j5eIaIMZX6m8niN8nylal03K/NXEVG1av4oZd4aIppm93uMBRGlEtEiInpbeZzs+2O9clwv\nJqIFyryk/d4I75NzZCA5PwaT82Noco4MJOdIizCz/Nn0B2AIgB4Alqvm3QdgmjI9DcC9yvQYAO8B\nIAD9AHypzC8CsE75v1CZLlSe+0pZlpTXjnb6PRvYJ6UAeijT9QB8D6BDMu8XJc48ZTodwJdK/K8C\nOEeZ/ySAK5TpKwE8qUyfA+AVZboDgCUAMgE0B7AWQKrytxZACwAZyjIdnH7fBvbL9QBeAvC28jjZ\n98d6ACWaeUn7vZE/7/9BzpHa/SHnx+B9IufH0PtGzpGB+2M95Bxp+p/ccbQRM38CYKdm9jgAzynT\nzwE4XTX/efb5AkB9IioFUA1gDjPvZOZdAOYAGKU8l8/MX7DviH5etS7XYubNzLxQmd4HYCWAMiTx\nflHe237lYbryxwCGA3hdma/dJ7X76nUAI5QrX+MAvMzMR5j5BwBrAPRR/tYw8zpmPgrgZWVZ1yKi\ncgBjAfxNeUxI4v0RRtJ+b4T3yTkykJwfg8n5UZ+cIw1L2u+OWaTi6LxGzLxZmf4FQCNlugzABtVy\nG5V54eZv1JnvGUpzie7wXUFM6v2iNDlZDGArfD9UawHsZuYaZRH1+/C/d+X5PQCKEf2+crOHANwE\n4ITyuBjJvT8AX2HpAyL6hoimKvOS+nsjEpIc05Dzo5qcH3XJOTKYnCMtkOZ0AKIOMzMRJeX4KESU\nB+ANANcx8151U/Fk3C/MfBxANyKqD+AtAO0cDskxRHQKgK3M/A0RDXU6HhcZxMybiKghgDlE9J36\nyWT83ojElqzHtJwfA8n5MZCcI0OSc6QF5I6j87Yot7yh/L9Vmb8JQIVquXJlXrj55TrzXY+I0uE7\nKb7IzG8qs5N+vwAAM+8GMBdAf/iaTtRe7FG/D/97V54vALAD0e8rtxoI4DQiWg9fE5nhAB5G8u4P\nAAAzb1L+3wpf4akP5HsjEk9SH9NyfgxNzo9+co7UIedIa0jF0XmzANRmaZoC4N+q+RcomZ76Adij\n3F6fDaCKiAqVbFBVAGYrz+0lon5KW/ULVOtyLSXWZwCsZOYHVE8l7X4hogbKlVQQUTaAkfD1bZkL\nYIKymHaf1O6rCQA+UtrczwJwDvkyqDUH0Bq+ztxfA2hNvoxrGfB1jp9l/TuLDTPfzMzlzFwJX6wf\nMfMkJOn+AAAiyiWierXT8B3vy5HE3xuRsJL2mJbzYzA5PwaTc2QwOUdaiF2QoSdZ/gD8E8BmAMfg\naw99CXztyv8LYDWADwEUKcsSgMfga7u/DEAv1Xouhq/T8hoAF6nm94Lvi7EWwF8AkNPv2cA+GQRf\nO/SlABYrf2OSeb8A6AJgkbJPlgO4TZnfAr4f8TUAXgOQqczPUh6vUZ5voVrXLcr7XgVVxi9lH3+v\nPHeL0+85in0zFHUZ45J2fyjvfYnyt6I25mT+3sif9/8g50jt/pDzY/A+kfNj+P0zFHKOrH3vco60\n4I+UNy+EEEIIIYQQQuiSpqpCCCGEEEIIIcKSiqMQQgghhBBCiLCk4iiEEEIIIYQQIiypOAohhBBC\nCCGECEsqjkIIIYQQQgghwpKKoxAJgIjqE9GVynQTInrd6ZiEEEIIN5BzpBDmkOE4hEgARFQJ39hN\nnRwORQghhHAVOUcKYY40pwMQQphiBoCWRLQYvoFt2zNzJyK6EMDpAHIBtAZwP4AMAJMBHAEwhpl3\nElFL+Aa/bQDgIIDLmPk7+9+GEEIIYTo5RwphAmmqKkRimAZgLTN3A3Cj5rlOAMYD6A3gbgAHmbk7\ngM8BXKAs8zSAa5i5J4AbADxuS9RCCCGE9eQcKYQJ5I6jEIlvLjPvA7CPiPYA+I8yfxmALkSUB2AA\ngNeIqPY1mfaHKYQQQthOzpFCGCQVRyES3xHV9AnV4xPw/QakANitXIkVQgghkomcI4UwSJqqCpEY\n9gGoF8sLmXkvgB+IaCIAkE9XM4MTQgghHCTnSCFMIBVHIRIAM+8AMJ+IlgP4UwyrmATgEiJaAmAF\ngHFmxieEEEI4Rc6RQphDhuMQQgghhBBCCBGW3HEUQgghhBBCCBGWVByFEEIIIYQQQoQlFUchhBBC\nCCGEEGFJxVEIIYQQQgghRFhScRRCCCGEEEIIEZZUHIUQQgghhBBChCUVRyGEEEIIIYQQYf0/p/ou\n06UMO6sAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "times = np.linspace(1000.*torb, 9000.*torb, Noutputs)\n", "a = np.zeros(Noutputs)\n", "e = np.zeros(Noutputs)\n", "for i,time in enumerate(times):\n", " sim.integrate(time, exact_finish_time=0)\n", " a[i] = sim.particles[2].a\n", " e[i] = sim.particles[2].e\n", " \n", "fig = plt.figure(figsize=(15,5))\n", "\n", "ax = plt.subplot(121)\n", "ax.set_xlabel(\"time\")\n", "ax.set_ylabel(\"semi-major axis\")\n", "plt.plot(times, a);\n", "\n", "ax = plt.subplot(122)\n", "ax.set_xlabel(\"time\")\n", "ax.set_ylabel(\"eccentricity\")\n", "plt.plot(times, e);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The semimajor axis seems to almost stay constant, whereas the eccentricity undergoes an oscillation. Thus, one might conclude the planets interact only secularly, i.e. there are no large resonant terms." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Speeding things up and extra accuracy**\n", "\n", "There are several performance enhancements one can make to WHFast. However, each one has pitfalls that an inexperienced user can unwittingly fall into. We therefore chose safe default settings that make the integrator difficult to misuse. **This makes the default WHFast substantially slower and less accurate than it can be**, so anyone looking to use it more seriously should check out its advanced settings in [Advanced Settings for WHFast](../AdvWHFast)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Common mistakes with WHFast**\n", "\n", "If you're getting odd output, check the following:\n", "\n", "1. The Wisdom-Holman algorithm assumes that the gravitational force on a planet from the central body dominates that from all other particles. Therefore, if you have close approaches (that violate this approximation), you will get spurious results. REBOUND provides a high order integrator for close approaches. You can try it with `sim.integrator = \"ias15\"`. You can also check for close approaches following [Exceptions.ipynb](../Exceptions).\n", "\n", "2. A symplectic scheme requires a constant timestep to guarantee some of its symmetry properties. So if you call `sim.integrate(time)`, and `time` is not a multiple of `sim.dt`, your last timestep will be different (in order to reach `time` exactly). Therefore, if you need equally spaced outputs, you can make your output times be multiples of `sim.dt`, or if it doesn't matter, you can pass an optional flag like this: `sim.integrate(time, exact_finish_time=0)`, which will integrate to the nearest timestep. \n", "\n", "3. If you're somehow modifying particles or adding forces, you should make sure to read [Advanced Settings for WHFast](../AdvWHFast)." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.0" } }, "nbformat": 4, "nbformat_minor": 1 }