ipython_examples/WHFast.ipynb

{
 "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": [
      "<rebound.Particle object, m=1.0 x=0.0 y=0.0 z=0.0 vx=0.0 vy=0.0 vz=0.0>\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 center 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",
      "<rebound.Particle object, m=1.0 x=0.0 y=0.0 z=0.0 vx=0.0 vy=0.0 vz=0.0>\n",
      "<rebound.Particle object, m=0.001 x=1.0 y=0.0 z=0.0 vx=0.0 vy=1.0 vz=0.0>\n",
      "<rebound.Particle object, m=0.001 x=1.800999000999001 y=0.0 z=0.0 vx=0.0 vy=0.7835163064519433 vz=0.0>\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",
      "<rebound.Particle object, m=1.0 x=0.003326154866766361 y=0.009674635911450204 z=0.0 vx=0.0005194654213133196 vy=0.0012200269278386914 vz=0.0>\n",
      "<rebound.Particle object, m=0.001 x=1.0032694180883746 y=0.0366289427053581 z=0.0 vx=-0.024395944012027243 vy=0.9999782071644221 vz=0.0>\n",
      "<rebound.Particle object, m=0.001 x=-1.5284252838557046 y=1.496351615582015 z=0.0 vx=-0.4950694773013606 vy=-0.4364888285516785 vz=0.0>\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": "iVBORw0KGgoAAAANSUhEUgAAAUwAAAEzCAYAAABaGjpLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VPW9//HXJ3uAJOxJIAn7DioaFoG61KVIrVitFe2idUFb7Xa7edv7895fa/fb9tdWrxaXutxWrbZUqihq1aIiSED2RUPYspAEAgkhe/L9/ZHBppCQQ2YyZyZ5Px+PeeTMmW/m+8lJ8p6zfc8x5xwiItK5GL8LEBGJFgpMERGPFJgiIh4pMEVEPFJgioh4pMAUEfEo6MA0s2wze93MtpnZVjP7ajttzMx+Y2b5ZrbJzM4Otl8RkXCLC8F7NAHfcM6tN7MUYJ2ZveKc29amzWXAuMBjFnB/4KuISNQIeg3TOVfinFsfmD4KbAeGn9BsIfC4a7Ua6G9mmcH2LSISTiHdh2lmI4HpwJoTXhoO7G/zvJCTQ1VEJKKFYpMcADPrB/wZ+JpzriqI91kMLAbo27fvORMnTgxRhSIirdatW3fQOTfkdL8vJIFpZvG0huUfnHN/aadJEZDd5nlWYN5JnHNLgCUAubm5Li8vLxQlioh8yMz2duX7QnGU3ICHge3OuV920GwZ8PnA0fLZQKVzriTYvkVEwikUa5hzgc8Bm81sQ2Ded4EcAOfcA8ByYAGQD9QAXwhBvyIiYRV0YDrn3gKskzYOuCPYvkRE/KSRPiIiHikwRUQ8UmCKiHikwBQR8UiBKSLikQJTRMQjBaaIiEcKTBERjxSYIiIeKTBFRDxSYIqIeKTAFBHxSIEpIuKRAlNExCMFpoiIRwpMERGPFJgiIh4pMEVEPFJgioh4pMAUEfFIgSki4pECU0TEIwWmiIhHCkwREY8UmCIiHoUkMM3sETMrM7MtHbx+gZlVmtmGwOPuUPQrIhJOcSF6n0eBe4HHT9HmTefc5SHqT0Qk7EKyhumcWwlUhOK9REQiVTj3YZ5rZhvN7EUzmxLGfkVEQiJUm+SdWQ+McM5Vm9kC4K/AuPYamtliYDFATk5OmMoTEelcWNYwnXNVzrnqwPRyIN7MBnfQdolzLtc5lztkyJBwlCci4klYAtPMMszMAtMzA/0eCkffIiKhEpJNcjN7ErgAGGxmhcB/AvEAzrkHgE8BXzSzJqAWWOScc6HoW0QkXEISmM656zp5/V5aTzsSEYlaGukjIuKRAlNExCMFpoiIRwpMERGPFJgiIh4pMEVEPFJgioh4pMAUEfFIgSki4pECU0TEIwWmiIhHCkwREY8UmCIiHikwRUQ8UmCKiHikwBQR8UiBKSLikQJTRMQjBaaIiEcKTBERjxSYIiIeKTBFRDxSYIqIeKTAFBHxSIEpIuJRSALTzB4xszIz29LB62ZmvzGzfDPbZGZnh6JfEZFwCtUa5qPA/FO8fhkwLvBYDNwfon5FRMImLhRv4pxbaWYjT9FkIfC4c84Bq82sv5llOudKQtG/9Gx1jc0UHq6hpLKOIzWNHKlt5MixBiprG2lqcbS44w+IizH6JcbRLymOlMQ4UpLiyUhLYnj/ZNJTk0iI014o6bqQBKYHw4H9bZ4XBuYpMOVDtQ3NbD9QxdaiSrYUVVFwsJp9FTWUVtW32z45PpaEuBhiDGLMMIOGphaONTTT3OJOam8GQ1MSGTW4LxMzUpmUmcLEjFTGp6eQnBDb3T+e9ADhCkzPzGwxrZvt5OTk+FyNdKfK2kZWFxxiVf5B1uyu4P3SoxzPuYF9Exg7tB/njRtCzsA+ZA/sw7D+yQzsG09acgJpyfEdri0656hvauFoXROVtY0cqKyjuLKW4iO1FB2uJb+8mj/l7aemoRloXSudlpXGrFGDmDV6ILkjBpCSFB+uxSBRJFyBWQRkt3meFZh3EufcEmAJQG5u7smrCRLV8suqeXFzCa9uL2VzUSUtrnVNMXfkAC6dnM7U4WlMHZ5GZloSZtalPsyMpPhYkuJjGZKSyNih/U5q09Li2H+4hu0lR9lYeIR3d1fw8FsFPPCPXcTGGDNHDuSSyelcMjmd7IF9gv2xpYew1t2KIXij1n2Yzzvnprbz2seBO4EFwCzgN865mZ29Z25ursvLywtJfeKfgvJqlm0sZvnmEt4vrQbg7Jz+zBs3hLljBjE9Z0BE7FusbWhm/b7DvJ1/kFe3l35Y66TMVBaeNYyrpg9naGqSz1VKKJjZOudc7ml/XygC08yeBC4ABgOlwH8C8QDOuQesdVXhXlqPpNcAX3DOdZqECszoVd/UzIqtpfxxzV5WF1RgBjNGDmTB1AzmT80kIy3yg2fPwWO8sq2UF7eUsH7fEWJjjAsnDOGa3Gw+OnEo8bH+h7x0ja+B2V0UmNGn/Gg9v397N0+t3U/FsQayByazaEYOnzoni/QoXjvbVV7Ns+sK+fO6QsqO1pOZlsRNc0exaGa29ndGIQWm+KrwcA0PrizgqbX7aWhu4dLJ6Xxm1gjmjR1MTEzX9kVGoqbmFt7YWc7Db+3mnYJDpCTGcd2sHG6eNyqqPxB6GwWm+KK0qo5fvfI+z64rxAyump7FbeePZvSQkw+09DSbCo/w4Ju7Wb65hPhY48Y5o/ji+WNI66M1zkinwJSwOlbfxJKVBSxZWUBTSwufmTWCxeeNZlj/ZL9LC7t9h2r45Ss7eW5jMSmJcXzxgrF8Ye5IkuJ1bmekUmBKWDjnWPpeET9+cQflR+v5+BmZfOdjE8kZpFNvthVX8fMVO3h9Zzk5A/vw/YVTuGDCUL/LknYoMKXb7a+o4btLN/PmBwc5K7s//+fyyZwzYoDfZUWcVfkH+Y/ntlBQfowF0zK4+/IpUXFWQG+iwJRu09zi+P3bu/nFy+8TY3DXZRP5zKwRPepgTqjVNzXz4MoCfvtaPvGxMfznJybzqXOyunwyvoSWAlO6RWlVHV97agPvFBzioolD+cGVU3vlfsqu2neohm89u5E1uytYMC2DH145jQF9E/wuq9framBG3FhyiRyv7Sjlm89sorahmZ9/6gytIXVBzqA+/PHW2Tz4ZgG/eHkn6/Ye5pefPou5Ywf7XZp0gYYqyEmaWxw/Xr6dmx7NIz01ib99eR7X5GYrLLsoNsa4/fwxLP3SXPolxvG5h9ewZOUuInnrTtqnwJR/cbSukVseW8vvVhbw2dk5LP3SnHYvXiGnb+rwNJbdOY/5UzP40fIdfPWpDdQGrpgk0UGb5PKh/RU13PzYWnaVH+OeK6fy2dkj/C6px+mbGMd915/N/7yxi/9+eSf5ZdU8cuMMHUWPElrDFAA2F1ay8L63Ka2q54mbZiosu5GZcceFY3nkhhnsPXSMq+9fRUF5td9liQcKTGHtngquf3A1yfGxLP3SHObogERYXDhxKE8tPpe6xmaueeAdNhdW+l2SdEKB2cu9+UE5n3/4XYakJPLM7ef2ijHgkWRaVhrP3H4uSfGxLFryDmsKDvldkpyCArMX+8f75dz8aB4jBvXh6dvO1fmVPhk9pB9//uIcMtKSuOnRtazfd9jvkqQDCsxeat3eCm57Io+xQ/vx1OLZDElJ9LukXi0jLYk/3jqbwSmJ3PDIu2wp0uZ5JFJg9kLbS6r4wu/XkpGaxGM3zaR/H408iQTpqUn84ZZZpCbF87mH1/BB6VG/S5ITKDB7mcLDNXz+kXfpkxDHEzfP0pplhMka0Ic/3jqLuNgYvvDoWg5Wt3+LYfGHArMXqWlo4tbH11HX2MwTN8/U3RAj1IhBfXno87kcrK5n8eN51DXq5PZIocDsJZxzfOuZTew8UMVvr5vOuPQUv0uSUzgzuz+//PRZrN93hG8/u0nDKCOEArOXuO/1fF7YXMJdl03URW2jxIJpmXzrYxNYtrGYh97c7Xc5ggKzV1i16yC/eOV9rjxrGLd+ZLTf5chp+NIFY/jYlHR++tIONuw/4nc5vZ4Cs4c7UtPAvz29kVGD+vKjq6bpikNRxsz42dVnkp6axJ1/XE9lbaPfJfVqCswezDnHd5du5mB1Pb9eNJ0+CbrWSjRK6xPPb6+fzoHKOr67dLPf5fRqCswe7M/ri1i++QDfuHQC07LS/C5HgnB2zgC+fsl4XthUwoubS/wup9cKSWCa2Xwz22lm+WZ2Vzuv32hm5Wa2IfC4JRT9SscqjjVwzwvbyB0xgMXnab9lT3DbeaOZOjyV//PcViprtGnuh6AD08xigfuAy4DJwHVmNrmdpk87584KPB4Ktl85tR8v3051XRM/umoasbpZWY8QFxvDT68+g8M1rR+GEn6hWMOcCeQ75wqccw3AU8DCELyvdNGagkM8s66QW88bzXidb9mjTBmWxm3njeaZdYW8s0tXNgq3UATmcGB/m+eFgXknutrMNpnZs2aWHYJ+pR3NLY67n9tK1oBkvvLRcX6XI93gKxeNY1haEj9cvo2WFp3QHk7hOujzN2Ckc+4M4BXgsY4amtliM8szs7zy8vIwlddzLH2viJ2lR/nugkkkJ8T6XY50g6T4WL49fyJbiqpY+l6R3+X0KqEIzCKg7RpjVmDeh5xzh5xzx68i8BBwTkdv5pxb4pzLdc7lDhkyJATl9R71Tc386pX3OSMrjcumZvhdjnSjK84cxplZafx8xU7dSC2MQhGYa4FxZjbKzBKARcCytg3MLLPN0yuA7SHoV07wv6v3UXSklu/Mn6gT1Hu4mBjjPy6fzIGqOp5YvcfvcnqNoAPTOdcE3AmsoDUI/+Sc22pm3zezKwLNvmJmW81sI/AV4MZg+5V/VdfYzP1v5DNv7GDm6p48vcKMkQOZO3YQD765W1c0CpOQ7MN0zi13zo13zo1xzv0wMO9u59yywPS/O+emOOfOdM5d6JzbEYp+5Z+WvlfEweoG7rhwrN+lSBjdccFYyo/W8+y6Qr9L6RU00qcHaGlxPPRmAVOHpzJ79EC/y5EwOnfMIM7K7s8D/9hFU3OL3+X0eArMHuD1nWXsKj/GrR8ZrX2XvYyZ8aULxlB4uJaXt5X6XU6Pp8DsAR5dtYdhaUksmJbZeWPpcS6alM7w/sk8+e4+v0vp8RSYUa74SC1v5R/kmtxs4mP16+yNYmOMa2dk8+YHB9l76Jjf5fRo+g+LckvfK8I5uPrsLL9LER99Ojeb2BjjyXf3d95YukyBGcWcczy7rpBZowaSM0g3NOvNMtKSuHDCUJa+V6jhkt1IgRnFNhVWsvvgMa1dCgCfODOT0qp63tt/2O9SeiwFZhR7dXspMQaXTE73uxSJAB+dOJSE2BiWbz7gdyk9lgIzir2yrZTcEQMZ0DfB71IkAqQkxXPe+MG8uLlEt+XtJgrMKFV4uIYdB45y8WTdMlf+af7UTIor69haXOV3KT2SAjNKvb6z9dJ3H52ozXH5p4+Ma72OgC4u3D0UmFFq7e4K0lMTGTOkr9+lSARJT01i9JC+rNp10O9SeiQFZpTK21NB7siBGgopJ5kzZhDv7q6gUWPLQ06BGYWKj9RSXFlH7ogBfpciEWjOmMEca2hmS1Gl36X0OArMKLRub+t5drkjdGUiOdmZ2f0B2KIDPyGnwIxCOw5UERdjTMjQHSHlZMPSkkhLjmdbsdYwQ02BGYU+KK1m5OC+JMTp1ycnMzOmDEtlm9YwQ07/cVEov6ya8en9/C5DItjkzFR2HDiqiwqHmAIzytQ3NbPn0DHGDtXmuHRs9JB+1De1UHq0vvPG4pkCM8oUHa6lxcFIXZ1ITmH4gGSg9e9FQkeBGWXKAmsM6alJPlcikSwrEJiFh2t8rqRnUWBGmeOBOTQl0edKJJIN7681zO6gwIwyZVV1AAxRYMopJMXHkpYcT3m19mGGkgIzylQcayAuxkhLjve7FIlwlbWNPP7OXr/L6FEUmFGmrrGF5PhYjSEX8UFIAtPM5pvZTjPLN7O72nk90cyeDry+xsxGhqLf3qi+qZnEeH3Oifgh6P88M4sF7gMuAyYD15nZ5BOa3Qwcds6NBX4F/DTYfnur+qYWEuNi/S5DpFcKxarKTCDfOVfgnGsAngIWntBmIfBYYPpZ4CLTNmWXNDS1EB+rRSedGzW4L+mpOjgYSqEIzOFA25shFwbmtdvGOdcEVAKDQtB3rxMXYzTrfi3iQdaAZIYFTi+S0Ii4nWFmttjM8swsr7y83O9yIk5ifAx1jRofLJ2raWimb0Kc32X0KKEIzCIgu83zrMC8dtuYWRyQBrR70xHn3BLnXK5zLnfIkCEhKK9nSYyLpa6x2e8yJAocq2+iT4L2d4dSKAJzLTDOzEaZWQKwCFh2QptlwA2B6U8BrzndB7RLkuJjqdcapnhwrKGJfolawwyloJemc67JzO4EVgCxwCPOua1m9n0gzzm3DHgYeMLM8oEKWkNVuiA1OY6G5hZqGproo80tOYXquib6KjBDKiRL0zm3HFh+wry720zXAdeEoq/ebmhK60U3yqrqGTlY/wzSvrrGZg7XNOooeYhF3EEfObXj/wBlus6hnEJJZes1BzLTdJQ8lBSYUeb4GmZp4CIcIu0pOdJ6laLM/roMYCgpMKPM8X+AQl22S06hOLCGOUxrmCGlwIwyqUnxDE1JJL+s2u9SJILll1UTH2s6cT3EFJhRaHx6Ch+UHfW7DIlgOw5UMXZoiu4sGmJamlFoXHo/8suqaWnRqazSvu0lVUzSfetDToEZhSakp1DT0Mx+3a9F2nH4WAOlVfVMzFRghpoCMwqdldMfgHV7D/tciUSiDfuPADB1WJrPlfQ8CswoNH5oCilJcazdo8CUk60uOERCbAzTcwb4XUqPo8CMQjExRu6IAazdU+F3KRKBVu+u4MzsNJJ14Y2QU2BGqdyRA8kvq+ag7goobVTXN7GlqJLZo3W52e6gwIxS541rvfTdGzt1zVD5p7fzD9Lc4jh3jAKzOygwo9TU4amkpyby6rZSv0uRCLJi6wHSkuOZMXKg36X0SArMKGVmXDwpnZUflOuCwgJAY3MLr24r5eJJ6cTH6l+7O2ipRrGLJ6dT09DMql0H/S5FIsCaggqq6pr42JR0v0vpsRSYUWzOmEGkJcez9L1iv0uRCPC3jcX0SYjlvPG6tUt3UWBGscS4WBaeNYwVWw9QWdvodznio+r6Jv62qZjLz8gkKV6nE3UXBWaUu+acbBqaWvjbRq1l9mbPbyympqGZRTNz/C6lR1NgRrmpw1OZkJ7CM3n7O28sPdZTa/czPr0f07P7+11Kj6bAjHJmxvWzcthYWMm6vRr50xttKjzChv1HuHZGDmbmdzk9mgKzB7gmN4u05HiWrCzwuxTxwQP/2EVKUhyfzs3yu5QeT4HZA/RJiONzs0fw8rZSdh885nc5EkYF5dW8uOUAnz93BClJ8X6X0+MpMHuIG+aMJD4mRmuZvcySlQXEx8Zw45xRfpfSKygwe4ghKYlcOyObZ/L2s0drmb3C3kPH+PP6Qj6dm8WQFN1/PBwUmD3Ily8aS3xsDL945X2/S5Ew+PmKncTFxPDlj47zu5ReI6jANLOBZvaKmX0Q+NruFUvNrNnMNgQey4LpUzo2NCWJm+eN4m8bi9lSVOl3OdKNNuw/wvObSrj1I6NIT9W9x8Ml2DXMu4C/O+fGAX8PPG9PrXPurMDjiiD7lFNYfP5oBvSJ54cvbMc53SStJ3LO8aPl2xncL4HF54/xu5xeJdjAXAg8Fph+DLgyyPeTIKUmxfNvl07gnYJDPLdBo396or+sL+Ld3RV8/ZLx9EuM87ucXiXYwEx3zpUEpg8AHV0mJcnM8sxstZkpVLvZZ2bmcFZ2f37w/DaO1DT4XY6E0KHqeu55YRvnjBjAdTM0DDLcOg1MM3vVzLa081jYtp1r3f7raBtwhHMuF7ge+H9m1uF2hJktDoRrXnm5ribeFTExxo+vmsaR2kZ++tIOv8uREPrB89uorm/ix1dNIyZGo3rCrdPAdM5d7Jyb2s7jOaDUzDIBAl/LOniPosDXAuANYPop+lvinMt1zuUOGaLLVHXVpMxUbpk3iiff3c/K9/XB0xO8tqOUv24o5ovnj2F8uu457odgN8mXATcEpm8AnjuxgZkNMLPEwPRgYC6wLch+xYOvXzKe8en9+MYzGzmkm6VFtbKqOr71zCYmZqTwpQvH+l1OrxVsYP4EuMTMPgAuDjzHzHLN7KFAm0lAnpltBF4HfuKcU2CGQVJ8LL9eNJ3K2ka+/ewmHTWPUi0tjq//aQPHGpq49/rput6lj4I6xOacOwRc1M78POCWwPQqYFow/UjXTcpM5a75E/n+89t4dNUevjBXQ+iizf3/2MXb+Yf4yVXTGDtUm+J+0kifXuALc0dy8aSh/PCF7br/T5R5Y2cZv3h5Jx8/I5NrZ2T7XU6vp8DsBcyMX117FiMH9+WOP6xnf0WN3yWJB/llR/nyH99jQkYqP7v6DF3rMgIoMHuJlKR4Hvp8Li0Obnksj+r6Jr9LklM4fKyBmx/LIzE+hoduyKWvTlCPCArMXmTk4L7cd/3Z5JdXc/sT66hv0v3MI1FtQzOLn8ij5Egdv/tcLsP7J/tdkgQoMHuZeeMG87Orz+Ct/IN89ckNNDW3+F2StFHf1Mxt/7uOvL2H+cWnz+ScEe1ez0Z8osDsha4+J4u7L5/MS1sP8O9/2azTjSJEU3MLX31yAyvfL+cnV03jE2cO87skOYF2jPRSN80bRWVtI7/++wfExRr3XDmNWA21801jcwvffGYjL209wN2XT+ZajROPSArMXuxrF4+jxTl++1o+1fXN/PLTZxIfq42OcKtrbOaOP6zn7zvK+M78idw0T+fKRioFZi9mZnzj0gn0TYzjJy/uoKa+ifs+c7ZGkoRRVV0jtzyWx9o9Fdxz5VQ+O3uE3yXJKWh1Qrj9/DHcc+VUXttZxmceWkP5UY07D4fiI7Us+t1q1u89zK8XTVdYRgEFpgDw2dkjuO/6s9laXMmV973NtuIqv0vq0fL2VHDFvW+xr6KGh27I5Qod4IkKCkz50IJpmTx7+xyaWxxX37+Kl7aUdP5NctqeXruP6x5cTb/EOP56xxwumDDU75LEIwWm/Iupw9NYdudcJmSkcPv/rue/lm2lrlEnuIfCsfomvvXMRr7z583MHj2I5+6Yp4tpRBkFppxkaGoST982m5vmjuLRVXv45P+sIr+s2u+yotrG/Uf4+G/e5Nn1hdx54Vh+f+MM0vrE+12WnCYFprQrMS6Wuz8xmUduzKW0qo5P/PYtHlu1h+YWneR+OhqbW7jv9Xyuvn8V9U0tPHnrbL75sQnE6fStqGSRPMojNzfX5eXl+V1Gr1dWVcc3n93EyvfLmZ7Tn59cdQYTMrQp2Zn1+w7z3b9sZseBoyyYlsGPP3mG1iojhJmtC9xn7PS+T4EpXjjn+OuGIn7w/Haqahu57fzR3HHhWPok6FTeE1XWNvLzFTv4w5p9pKck8X8XTuFjUzL8Lkva6Gpg6q9dPDEzPjk9i/PHD+WeF7Zx3+u7eCavkH+7ZDzX5GZrWCWtF8544p293Pt6PlW1jdw4ZyTfuHSC7h3eg2gNU7okb08FP1q+nfX7jjA+vR/fmT+Rj04c2isvctvS4li2sZj/fnknhYdrmTd2MHddNpGpw9P8Lk06oE1yCTvnHC9tOcBPX9rBnkM1TM5M5YsXjGHBtMxescZZ39TMc+8V87uVu9hVfowpw1K567KJfGScbg8d6RSY4puGphb+uqGIB/6xi4LyY4wY1IdbPzKaK6cP75Gbo5U1jTy1dh+PvL2b0qp6JmemcvsFY7h8WiYxveCDoidQYIrvWlocL28r5f438tlYWEnfhFiuOGsYi2bkcEZWWlRvrjvnWF1QwdNr9/HilgPUN7UwZ8wgbj9/DB8ZNziqf7beSIEpEcM5x3v7j/Dkmn08v6mE2sZmJmak8IkzhzF/agZjhvTzu0RPnHNsLznKiq0HeG5DEXsO1ZCSFMeVZw3n2hnZ2kcZxRSYEpGq6hp5bkMxf1lfyHv7jgAwPr0f86dmcv74IZyRlRZR1+Csb2pmw74j/H1HGS9tOcC+ihrMYNaogVw7I5v5UzJJTtDl76KdAlMiXkllLSu2HODFLQdYu6eCFgd9E2KZOWogc8YM5uwR/ZmUmRrWczur6hrZVlzFu7srWF1wiHV7D1Pf1EJ8rDF37GDmT8ng4snpDO6XGLaapPv5Ephmdg3wX8AkYKZzrt10M7P5wK+BWOAh59xPvLy/ArPnqjjWwJqCQ7y96yCrdh2ioPwYAGYwanBfpgxLY2JGCtkD+5A9IJnsgX0Y1DehS/sKW1oc5dX1FB6uofBwLXsO1rCtpJLtJUfZF7hHuxlMykhl9uhBnDtmELNGDyQ1SaNyeiq/AnMS0AL8Dvhme4FpZrHA+8AlQCGwFrjOObets/dXYPYepVV1bCqsZGtxJVuLq9hWXEXRkdp/aZMUH8PAPgmkJseTmhxPWnI8yfGxHM9QA1ocVNc3cbSukaN1TRyta6L8aD0Nbe6OaQYjB/VlcmYqk4elMikzhbNzBtC/T0IYf2Lxky8jfZxz2wOdn6rZTCDfOVcQaPsUsBDoNDCl90hPTeKSyUlcMjn9w3nV9U2ta4UVtew/XEPR4VoO1zRSVddIZW0j+ytqPrz03PGPfQP6JcWRkhhPzsA+pCTFMyQlkeEDksnqn0zWgGSGD0jWkE7pknD81QwH9rd5XgjMCkO/EuX6JcYxMSOViRmpfpciAngITDN7FWjvygHfc849F+qCzGwxsBggJ0e3GhWRyNFpYDrnLg6yjyIgu83zrMC8jvpbAiyB1n2YQfYtIhIy4TgBbi0wzsxGmVkCsAhYFoZ+RURCKqjANLNPmlkhcC7wgpmtCMwfZmbLAZxzTcCdwApgO/An59zW4MoWEQm/YI+SLwWWtjO/GFjQ5vlyYHkwfYmI+C1yxqSJiEQ4BaaIiEcKTBERjxSYIiIeKTBFRDxSYIqIeKTAFBHxSIEpIuKRAlNExCMFpoiIRwpMERGPFJgiIh4pMEVEPFJgioh4pMAUEfFIgSki4pECU0TEIwWmiIhHCkwREY8UmCIiHikwRUQ8UmCKiHikwBQR8UiBKSLikQJTRMRiezKkAAAG/klEQVSjoALTzK4xs61m1mJmuadot8fMNpvZBjPLC6ZPERG/xAX5/VuAq4DfeWh7oXPuYJD9iYj4JqjAdM5tBzCz0FQjIhLBwrUP0wEvm9k6M1scpj5FREKq0zVMM3sVyGjnpe85557z2M8851yRmQ0FXjGzHc65lR30txhYDJCTk+Px7UVEul+ngemcuzjYTpxzRYGvZWa2FJgJtBuYzrklwBKA3NxcF2zfIiKh0u2b5GbW18xSjk8Dl9J6sEhEJKoEe1rRJ82sEDgXeMHMVgTmDzOz5YFm6cBbZrYReBd4wTn3UjD9ioj4Idij5EuBpe3MLwYWBKYLgDOD6UdEJBJopI+IiEcKTBERjxSYIiIeKTBFRDxSYIqIeKTAFBHxSIEpIuKRAlNExCMFpoiIRwpMERGPFJgiIh4pMEVEPFJgioh4pMAUEfFIgSki4pECU0TEIwWmiIhHCkwREY8UmCIiHikwRUQ8UmCKiHikwBQR8UiBKSLikQJTRMQjBaaIiEdBBaaZ/dzMdpjZJjNbamb9O2g338x2mlm+md0VTJ8iIn4Jdg3zFWCqc+4M4H3g309sYGaxwH3AZcBk4DozmxxkvyIiYRdUYDrnXnbONQWergay2mk2E8h3zhU45xqAp4CFwfQrIuKHUO7DvAl4sZ35w4H9bZ4XBuaJiESVuM4amNmrQEY7L33POfdcoM33gCbgD8EWZGaLgcWBp/VmtiXY9wyBwcBBv4sIUC3tUy3tUy3tm9CVb+o0MJ1zF5/qdTO7EbgcuMg559ppUgRkt3meFZjXUX9LgCWB985zzuV2VmN3i5Q6QLV0RLW0T7W0z8zyuvJ9wR4lnw98G7jCOVfTQbO1wDgzG2VmCcAiYFkw/YqI+CHYfZj3AinAK2a2wcweADCzYWa2HCBwUOhOYAWwHfiTc25rkP2KiIRdp5vkp+KcG9vB/GJgQZvny4HlXehiSRdLC7VIqQNUS0dUS/tUS/u6VIu1v9tRREROpKGRIiIeRVRgRspQSzO7xsy2mlmLmXV4VM/M9pjZ5sD+2y4ddQthLd0+/NTMBprZK2b2QeDrgA7aNQeWyQYzC+kBvs5+TjNLNLOnA6+vMbORoez/NGu50czK2yyLW7qpjkfMrKyjU/Cs1W8CdW4ys7O7ow6PtVxgZpVtlsnd3VRHtpm9bmbbAv8/X22nzekvF+dcxDyAS4G4wPRPgZ+20yYW2AWMBhKAjcDkENcxidbztN4Ack/Rbg8wuJuXSae1hGOZBPr5GXBXYPqu9n4/gdequ2lZdPpzAl8CHghMLwKe9rGWG4F7u/PvI9DPecDZwJYOXl9A66ASA2YDa3ys5QLg+TAsk0zg7MB0Cq1Dt0/8/Zz2comoNUwXIUMtnXPbnXM7Q/meXeWxlnANP10IPBaYfgy4shv6OBUvP2fbGp8FLjIz86mWsHDOrQQqTtFkIfC4a7Ua6G9mmT7VEhbOuRLn3PrA9FFaz9A5cYThaS+XiArME0TDUEsHvGxm6wIjlPwSrmWS7pwrCUwfANI7aJdkZnlmttrMQhmqXn7OD9sEPnwrgUEhrOF0agG4OrC596yZZbfzejhE0v8MwLlmttHMXjSzKd3dWWC3zHRgzQkvnfZyCeq0oq4I91DLYOrwYJ5zrsjMhtJ6LuqOwCesH7WExKlqafvEOefMrKNTLEYElsto4DUz2+yc2xXqWqPA34AnnXP1ZnYbrWu+H/W5Jr+tp/Xvo9rMFgB/BcZ1V2dm1g/4M/A151xVsO8X9sB0YR5q2dU6PL5HUeBrmZktpXUz7bQDMwS1hGSZdFaLmZWaWaZzriSw6VLWwXscXy4FZvYGrZ/uoQhMLz/n8TaFZhYHpAGHQtD3adfinGvb70O07gP2Q8j+PoLVNrScc8vN7H/MbLBzLuRjzM0sntaw/INz7i/tNDnt5RJRm+QWRUMtzayvmaUcn6b1gJVfFwoJ1zJZBtwQmL4BOGnt18wGmFliYHowMBfYFqL+vfycbWv8FPBaBx+83V7LCfvDrqB1P5oflgGfDxwVng1Uttm1ElZmlnF8n7KZzaQ1g0L+gRbo42Fgu3Pulx00O/3l0t1Hq07zyFY+rfsUNgQex492DgOWn3B0631a11q+1w11fJLW/Rn1QCmw4sQ6aD06ujHw2NoddXitJRzLJNDHIODvwAfAq8DAwPxc4KHA9Bxgc2C5bAZuDnENJ/2cwPdp/ZAFSAKeCfwtvQuM7sa/185q+XHgb2Mj8DowsZvqeBIoARoDfys3A7cDtwdeN1ov4r0r8Dvp8MyPMNRyZ5tlshqY0011zKP1GMOmNnmyINjlopE+IiIeRdQmuYhIJFNgioh4pMAUEfFIgSki4pECU0TEIwWmiIhHCkwREY8UmCIiHv1/dklIGL+hxyoAAAAASUVORK5CYII=",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1089cde10>"
      ]
     },
     "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": "iVBORw0KGgoAAAANSUhEUgAAAUwAAAEzCAYAAABaGjpLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8VfX9x/HXJwl77xH2UECqgBEQREFxUSvOVltt1Sq/Vm21W+tqta0d2tpWW0Fti/0p2lpQqiiCCjgYhr1k7z3DCElI8v39kau/CPcm5+aee+56Px+PPHLuOd97vp9cwjtnfc8x5xwiIlK9rEQXICKSKhSYIiIeKTBFRDxSYIqIeKTAFBHxSIEpIuJRzIFpZh3N7D0zW2Fmy83srjBtzMz+ZGZrzWyJmQ2ItV8RkaDl+LCOUuAHzrkFZtYImG9m05xzKyq1uRToGfoaBPw19F1EJGXEvIXpnNvhnFsQmj4MrARyT2g2GnjeVZgDNDWzdrH2LSISJF+PYZpZF6A/MPeERbnAlkqvt3JyqIqIJDU/dskBMLOGwH+Au51zh2JYzxhgDECDBg3O7NWrl08ViohUmD9//l7nXKto3+dLYJpZLSrC8gXn3MQwTbYBHSu97hCadxLn3DhgHEBeXp7Lz8/3o0QRkc+Y2aaavM+Ps+QGPAesdM79PkKzycDXQ2fLBwMFzrkdsfYtIhIkP7YwhwI3AkvNbFFo3k+BTgDOuaeBKcAoYC1QCNzsQ78iIoGKOTCdcx8AVk0bB9wRa18iIomkkT4iIh4pMEVEPFJgioh4pMAUEfFIgSki4pECU0TEIwWmiIhHCkwREY8UmCIiHikwRUQ8UmCKiHikwBQR8UiBKSLikQJTRMQjBaaIiEcKTBERjxSYIiIeKTBFRDxSYIqIeKTAFBHxSIEpIuKRAlNExCMFpoiIRwpMERGPFJgiIh75Ephm9jcz221myyIsH25mBWa2KPT1oB/9iogEKcen9fwDeBJ4voo27zvnLvOpPxGRwPmyhemcmwXs92NdIiLJKshjmGeb2WIze9PMTguwXxERX/i1S16dBUBn59wRMxsFvAr0DNfQzMYAYwA6deoUUHkiItULZAvTOXfIOXckND0FqGVmLSO0Heecy3PO5bVq1SqI8kREPAkkMM2srZlZaHpgqN99QfQtIuIXX3bJzWwCMBxoaWZbgYeAWgDOuaeBa4Bvm1kpcAy4zjnn/OhbRCQovgSmc+76apY/ScVlRyIiKUsjfUREPFJgioh4pMAUEfFIgSki4pECU0TEIwWmiIhHCkwREY8UmCIiHikwRUQ8UmCKiHikwBQR8UiBKSLikQJTRMSjoO64LhKz8nLH/sISdh8qZvfhIgqOHQcgJyuL7CzIrvQ9J8uoWyubDs3q0aphHbKyLMHVSzpQYEpSOV5Wzupdh1m8pYDl2wvYWVDE7sPF7Cg4xt4jJTVeb6fm9encoj4dmtWjQ7OK791bNaRX20bkZGtHS7xRYErCOOfYtK+QxVsPsnhLAYu2HGDB5oNx6Wvz/kI27y8Mu+zcU1oxsEszBnZtwekdmlC3VnZcapDUp8CUQJWWlTNv437eXLqTt5bvZM/h4kSXxKzVe5i1es9nr/t3asrZ3Vowsk8b+ndsSujpKiIKTIm/0rJy5qzfzxtLd/Dyx5spT/KHkyzcfJCFmw/ylxnraN+kLtfmdeSK/rl0bdkg0aVJglkyP1onLy/P5efnJ7oMqaElWw/ywpzNvL5kO0dLyhJdTsx6tW3EVwd14rLT29O8Qe1ElyMxMLP5zrm8qN+nwBQ/OeeYsXoPY2euY876/YkuJ24uPq0Nd4zowekdmia6FKmBmgamdsnFFyWl5UxevJ2n3lvLhr1HE11O3E1dvoupy3dx3imtuGtkTwZ0apbokiQACkyJyZHiUl6cu4mnZ65n/9GaX/aTqmau3sPM1XsY2LU5P7r4VM7q0jzRJUkcKTClRpxzTF68nYcmL+dg4fFEl5Nw8zbs59qnZ9OvY1N+Oqo3A7sqONORAlOitmbXYe6duJT8TQcSXUrSWbTlIF8eO5trzuzA/V/sTdP6OjmUThSY4tmR4lL+9M4axs1an+hSkt4r87fyyvyt/Pn6/lx2ejtdy5kmfBkTZmZ/M7PdZrYswnIzsz+Z2VozW2JmA/zoV4LhnOP1JdsZ8ug7CssofWfCQq5/Zg7bDh5LdCniA78G0f4DuKSK5ZcCPUNfY4C/+tSvxFlhSSnfe3kRd764kENFpYkuJyXNWb+fob9+l/EfbaQs2a/alyr5EpjOuVlAVRfdjQaedxXmAE3NrJ0ffUv8rN19hBGPzeDVRdsTXUpaeGjycq5/Zg4HCzPvaoJ0EdQxzFxgS6XXW0PzdgTUv0Rp8uLtfHfCwkSXQbeWDTilTSNym9WjXZO65DatR/um9WhUN4fsLMMwzMAMsiw0jXHwWAk7C4rYdaiInQXFbD94jEVbDrJq1+GE/jzzNuxnwCPTePOuczm1baOE1iLRS7qTPmY2horddjp16pTgajJPcWkZv3pjJeNnbwq872E9WzK0R0vyOjejb25sdw1q26Quvdo2DrvMOcf2giIWbj7A/E0H+Hf+Vo4UB3e4odzBxU/M4ukbBnBJX+1opRLfhkaaWRfgdedc3zDLxgIznHMTQq9XAcOdc1VuYWpoZLB2FhQx5p/5LNlaEEh/Vw3IZXDXFpzZpRndWjZI6Jnk4tIyVmw/xPxNB/jHRxvZeiCYkzR3jOjODy48VTc4DliyD42cDNxpZi8Bg4CC6sJSgrVlfyGXP/kBB+J8EfpNQ7owul97+iXZbdPq5GTTv1Mz+ndqxq3DunG46DjTVuxi3Kz1fLIzfrvxT723jvyNB3juprNoWCfpdvjkBL5sYZrZBGA40BLYBTwE1AJwzj1tFf8znqTiTHohcLNzrtpNR21hBmP9niOc//jMuK3/6gEduKJ/e87u1iIl725eUHicqct38vDrK+K2696hWT2m3DWMxnVrxWX98nm6W5HUyOpdh7noD7N8X2+z+rV45Iq+jOzdJq3uYL73SDFjZ67jmfc3+L7u3Kb1eOvuYTRSaMadAlOitmxbAZf9+QNf19m7XWPuG9WboT1aJNUut9+OFJfyz9mb+M1bn/i63raN6/LOD86jgXbP40qBKVFZuPkAV/7lI9/WN7RHC+4eeUrG3a2nuLSMf+dv5f5Xww5yq5FWjeow80fDqV9boRkvCkzxzM/d8DM7N+NnXzqNL3Ro4sv6UlV5ecXdm+5+eZEv62vRoDYf/OR86tVOn8MZyaSmgZl6R+AlJvuPlvi2Gz7+loH859tDMj4sAbKyjCv657Ly4Uu4cXDnmNe372gJVzz1ISWl5T5UJ35RYGaQ42XlXDdudsz/Cf/n3G6sfPgSzjullU+VpY96tbN55Iq+TP/+uTGva9Wuwzz4mn+7+hI7BWYGefC1ZazedaTG7+/Soj5TvjuMe0f11q5iNXq0bsSGR0fxiytOGscRlZc+3sKkhVt9qkpipcDMEP+cvZEJ87ZU2y6Sh77Uh3d+MJw+7cMPN5STmRk3DO7M4ocuokOzejVez/deXszKHYd8rExqSoGZAT5cu5cHXlte4/dP//553Dy0K9kavlcjTerV4v0fj+D+L/au8Tou/eP7FBzTo0ASTYGZ5vYcLuZrz86t0XsHdm3Okp9dRI/WDX2uKvOYGbcO68ak24fUeB0X/2EW5bqfZkIpMNPcvROX1Oh9d4/syUu3DdZQPZ/179SMBQ9cWKP37jxUxBPTV/tckURDgZnG3lq2k+krd0f9vme+nsfdI0/RHXTipHmD2qz71SiG9WwZ9Xv/9O5a1iT4np6ZTIGZpgqOHedb/zs/6vdN//55XNinTRwqksqys4x/fnMQ942K/rjm5U9+SDIPOElnCsw09fPJ0Z/kmffTC3S8MmC3nduNx689I6r3HDtexssf1/yKB6k5BWYa+nDtXiYu3BbVe+b+9AJaN64bp4qkKlef2YHffzm60Lxn4lI9GygBFJhp5lhJWdRnxef+9ALaKCwT6qoB0Yfmzf/4OE7VSCQKzDTzv3OiexbPnHsVlski2tBcuPkgCzYfiGNFciIFZho5VlLGL6es9Nx+9r3n07aJwjKZXDWgA49FcUzzqr98pBNAAVJgppHnZ2/03PaDn4ygXZOaD9eT+LnmzA7cOaKH5/bfSYLHIWcKBWaaKDpexqNverv79wu3DqJDs/pxrkhi8cOLT/V8xcLrS3ZoBFBAFJhpYsK8zZ7a/ejiUxnaI/oLpiV4b941zHPbp95bG8dK5FMKzDRQdLyMn/93RbXthvVsye3DuwdQkfihVnYW+feP9NT28WkaMhkEBWYa+He+t4uYn77hzLR+MFk6atmwDi/eOshT2/dWRT8MVqKjwEwDXm7d9q6eRJiyhvRoyc1Du1Tb7ua/67rMeFNgpjgvN5b97TWn062Vhjymsoe+dJqndlsPFMa5ksymwExxY2euq3J5o7o5XHtmh4CqkXj66J7zq23zlbFzAqgkc/kSmGZ2iZmtMrO1ZnZPmOU3mdkeM1sU+rrVj34zXXm549VF26ts8/p3ztFxyzTRvmk9rhqQW2WbbQePBVRNZoo5MM0sG3gKuBToA1xvZn3CNH3ZOdcv9PVsrP0KzNmwr8rld47oQecWDQKqRoLw2DXVjwLasPdoAJVkJj+2MAcCa51z651zJcBLwGgf1ivV+MXrVQ+D/M4F3keLSGrIyjJeGjO4yjZffUa75fHiR2DmApWva9kamneiq81siZm9YmYdfeg3oxUdL2NFFSd8XrxtEHVy9CjcdDS4Wwsa1418xcOOgqIAq8ksQZ30+S/QxTl3OjANGB+poZmNMbN8M8vfs2dPQOWlntnrI++ON6tfiyHdNZonnc2+94Iql+87UhxQJZnFj8DcBlTeYuwQmvcZ59w+59yn/4LPAmdGWplzbpxzLs85l9eqVSsfyktPU5ftjLhs0u1DA6xEEqFBnZwqnwl0+wsLAqwmc/gRmB8DPc2sq5nVBq4DJlduYGbtKr28HPB+DzIJ66UIjyho3agOXVrqRE8mGHtjxO0O5m7YH2AlmSPmoR/OuVIzuxOYCmQDf3POLTezh4F859xk4LtmdjlQCuwHboq130x2vKw84rK/3XRWgJVIItWvnUO3lg1Yr7PigfFlrJxzbgow5YR5D1aavhe414++BJZvj3yyp29ukwArkUR75dtDGPDItLDL9h8toXmD2gFXlN400icFRbrZxou3ebtJg6SPqgLxmffXB1hJZlBgpqAX5oa/9+XZ3VoEXIkkgxk/HB52/l9nVD1sVqKnwEwTv7ryCxoCmaF0ki84CswUc7S4NOz8q8+seoyxpLcHLws3Gln8psBMMZFGcWhUT2a7aUiXsPN1Abu/FJgpZv2eIyfNq+p6PMkMWVnhD8es2nk44ErSmwIzxfxzzqaT5g0/VSOiBB760sm75W8tjzwiTKKnwEwx76/Ze9I87Y4LwI2DO5807+UII8KkZhSYKe7P1/dPdAmSJHKyT/7vXFwaeVSYRE+BmeIuOq1NokuQJPK9kackuoS0psBMcdodl8puO7droktIawrMFNamcZ1ElyBJpn5tPUo5nhSYKewBXawsEigFZgrT2HEJZ0CnpokuIW0pMFNYi4baJZeTfe9CnfiJFwWmSJrJ69w80SWkLQWmSJqpV1tXTsSLAjOFOOc+m76iX/sEViKSmRSYKeRIpVu7XdK3bQIrEclMCswUsv3g/9/arX3TegmsRCQzKTBTyI6CY59NG7q7ukjQFJgppPKNFBZsPpDASkQykwIzhWRXembPvyI8OVJE4keBmUIqXy5S1bPJRSQ+FJgppF2TuokuQVJAaZnugRkvvgSmmV1iZqvMbK2Z3RNmeR0zezm0fK6ZdfGj30zTurECU6q3cd/RRJeQtmIOTDPLBp4CLgX6ANeb2Ym30fkmcMA51wP4A/CbWPvNRA1OGMFRWBL+kbuS2Vbu0IPP4sWPLcyBwFrn3HrnXAnwEjD6hDajgfGh6VeAC8xM18VE6cSPbPa6fQmqRJLZu5/sTnQJacuPwMwFKp+y3RqaF7aNc64UKAB0b7IYjZ21PtElSBKatHBboktIW0l30sfMxphZvpnl79mzJ9HlJLV5G/YnugRJMsd1wieu/AjMbUDHSq87hOaFbWNmOUATIOz+pHNunHMuzzmX16qVnrddnco35BDJ36gBDfHkR2B+DPQ0s65mVhu4Dph8QpvJwDdC09cA7zr9T/fFB2tPfk65ZK4J8zZ/7nWvto0SVEl6ijkwQ8ck7wSmAiuBfznnlpvZw2Z2eajZc0ALM1sLfB846dIjqZlfv/lJokuQJOGcY/Li7Z+bN6JX6wRVk558ecScc24KMOWEeQ9Wmi4CrvWjL/k8jfiRT63bc/L1l2d00PN9/JR0J30kersPF1XfSNLetBW7TprXr6MC008KzBRzTo+WJ8178NXlCahEks1v3jr58IyeXe8vBWaKuXVY15PmvbV8ZwIqkWSy7eCxsPM1PsRfCswUM7Br+CcCLt5yMOBKJJk8+e6aRJeQERSYKaZ+7fDn6a4bNyfgSiRZFJaUMmHeyfdH7d2ucQKqSW8KzDRx7HiZbuuVoSINhbxjRPeAK0l/CswUVCs7/HGpX05ZGXAlkmjOOe6btCzssuGn6hpMvykwU9Bj154Rdv7fP9yooZIZpqqRXg3r+HKZtVSiwExBVT2T/PnZmwKsRBLtl29oryJICswUVCcnO+KyhyYv11Zmhli+vYBPdoa/WfD3Lzwl4GoygwIzRbWv4vk+f5mxLsBKJBGcc9zw7NyIy8ec2y3AajKHAjNFPXr16RGX/W7qKsrLtZWZzqYs3cmBwuMRl9etFXkvRGpOgZmizu158hDJyu6duDSgSiRoRcfLuOPFBRGXX3xamwCrySwKzBRV3ZC3l/O3cLgo8haIpK6xM6t+NMm9l/YOqJLMo8BMYX+8rl+Vywf96p2AKpGgbD94jD9MX11lmy4tGwRUTeZRYKawy05vX+XywpIy3ll58i2/JHVVd6hlWDWHaiQ2CswUlp1l1T6C4Jvj8zlarOeXp4MpS3cwc3XVDwZ84itV73VIbBSYKe7pG86sts15v5sR/0IkrjbvK+T2FyKf6PlUi4a6/2U8KTBTnJfjVXuPFPPK/K0BVCPxUFJazjVPf1Rtuwm3DQ6gmsymwEwDv7yyb7Vtfvjvxew/WhJANeK3R99cye7DxdW2O7t7iwCqyWwKzDTwlbyO1TcCBjwyjZJS3QIulUxbsYu/f7ix2nY/vEhDIYOgwEwDOdlZ3Dy0i6e2l/35fY01TxFbDxRy2/P5ntreOkxDIYOgwEwTP7mkl6d2q3cd4WeT9dC0ZFdQeJzzH5vpqe2X8zpoKGRAFJhpom6tbO4e2dNT2/GzNzFh3uY4VyQ1dbS4lMuefJ8Sj3fQf+SK6o9hiz8UmGnkO+d7C0youAB69rp9caxGaqK4tIwbn5vLlv3hnwJ5ogcu61Pl7f7EXzEFppk1N7NpZrYm9L1ZhHZlZrYo9DU5lj4lsuws4xdRbG1c/8wcVkW4n6IEr7SsnDteWMCCzd6fAHqLx2PX4o9YtzDvAd5xzvUE3gm9DueYc65f6OvyGPuUKnxtUKeo2l/8xCyWbi2IUzXiVXm548f/WcL0lbs9v+elMYP13PGAxRqYo4HxoenxwBUxrk9iZGa88q2zo3rPl578gHkb9sepIqlOebnjgdeWMXFB+Kc/htOsfi0Gd9N1l0GLNTDbOOd2hKZ3ApFuxFfXzPLNbI6ZKVTjLK9Lc87qEvboSERfHjubaSt0o46gHSsp45vjP+aFudGdhJt0+9A4VSRVqTYwzWy6mS0L8zW6cjtXcXFfpAv8Ojvn8oCvAk+YWcQHJpvZmFC45u/ZU/WNBiSyv988MOr33PZ8vs6eB2j3oSIu+eMs3lsV3e/5nSN66BZuCVJtYDrnRjrn+ob5eg3YZWbtAELfwx6Acc5tC31fD8wA+lfR3zjnXJ5zLq9Vq1Y1+JEEKh6x+tKY6McW3ztxKb+b+kkcKpLKVmw/xKBH32HTvsKo3/sDjepJmFh3yScD3whNfwN47cQGZtbMzOqEplsCQ4EVMfYrHgzu1oIRp0b/R+ep99ZxyROzdFu4OHn3k12M+tP71GTAVf79I3WiJ4FiDcxfAxea2RpgZOg1ZpZnZs+G2vQG8s1sMfAe8GvnnAIzIOO+nlej932y8zCnPTSVFdsP+VxR5iovdzz7/npu+Ye34Y4neu4bebTU7dsSypJ5XHFeXp7Lz6/ZL5f8v7W7jzDy996G2YXzwGV9uGVoF23ZxGDL/kLufHEBi2t4CdcXT2/HU18d4HNVmcvM5ofOq0RFI30yQI/WDWMaPvfI6ys4//GZFBzTQ9Wi5ZzjhbmbGPbb92oclgCPX3uGj1VJTSkwM8SNgztz3VnebgMXzoa9Rznj528zdflO3e3Io20Hj/G1Z+dy36RlMa3n/R+P0M01koQCM4P86sovcFr7xjGt43/+OZ8Rj81g7W4NqYzEOcdL8zYz9Nfv8lGM4/Xf+O45dGxe36fKJFYKzAySlWVMvH1IzOvZuK+Qkb+fxb0Tl1BQqN30Tznn+GDNXkY8NoN7qnm6oxcv3jqI09o38aEy8YsCM8PUyclmyc8u8mVdE+Zt4YyH32b8RxspK8/s3fT5m/ZzzdOzueG5uWyswbWVJ3ryq/0Z0kOPzE02OkueoXYWFDH40Xd8Xedvrz6dy/u1z6jjbcu2FfDY26uYEeVonar87Et9uGloV9/WJyer6VlyBWYGW7PrMBf+YZbv6/3uBT25YXAnWjeq6/u6k4FzjgWbDzJ25jre9nn8/beHd/d893ypOQWm1MjWA4Wc85v34rLuC/u04a4LetI3Nz2Ow+05XMzEBVt5euY6DsTh2O2PLj6VO0b08H29cjIFptTYgaMl9H9kWtzWXyvbeOCyPozs3Yb2TevFrZ94OF5WzoxVe/hX/pa43s3pj9f1Y3S/3LitXz5PgSkxOVZSxvmPz2BHQVFc+8ltWo+vnNWRC/u0oVfbRkk5euhQ0XHmrNvHR+v28Y+PNsa9v5fHDGaQ7m0ZKAWmxKy0rJw7X1zIW8t3BtbnjYM7M6JXK/q2b0Lrxok55ll0vIz5mw7w4dq9zFy9h+UBjp+f/v1z6dG6UWD9SQUFpvjCOccT09fwx3fWJKT/YT1b0je3Cae1b8xp7ZvQuXl9srL82QotK3dsP3iMjfuOsnHvUTbsLWT+5gMs3uL9GTp+aVq/Fm/ffW7C/khkOgWm+OrDtXv52rNzE13GZ/rmNia3aT2aN6hNs/q1P/vetH4tysodRaXlFB8v++x7cej70ZIyNu8vZPWuwzW692Q8XDUgl19c0Zf6tXMSXUrGUmCK7w4cLeGW8R+zMIqnGErV/n7TWYzo1TrRZWQ83a1IfNesQW0mfntIVI/ulfAGdmnO/PtHKixTnPYJpEpmxg2DOzOkewuuf2YOuw4VJ7qklPPba07n2jM7JOUVARIdbWGKJ91aNeSDn5zPPZdqFIpXAzo1ZdaPRvDlvI4KyzShLUzxrFZ2Ft86rztfzuvI42+vivrRsJnCDJ6/ZSDDeuohfulGW5gSteYNavPLK7/Aez8czoBOTRNdTlJ54iv9WPfLUQrLNKUtTKmxri0bMPH2oeRv3M89E5eydveRRJeUMPd/sTc3nt2ZOjmZc6emTKTAlJjldWnOtO+dy8zVe3jugw28v2ZvoksKzA8vOoUbz+5Ck3q1El2KBECBKb4wM4af2prhp7Zmy/5CXpy3mb/OWJfosuJiYNfmfPu87px7SiuyfRqFJKlBF65L3BSXlvHWsp2MnbmeFTtS//nm3z2/B18Z2IncFLvjkpyspheuawtT4qZOTjaj++Uyul8ua3cfZvrK3UxZuoMlMTxuNmjXndWR83u1ZkSv1tTK1jnSTKfAlED0aN2IHq0b8a3zunOo6DgfrNnL28t38uqi7Yku7XMGdm3ORX3acN4prejRuqGun5TPiWmX3MyuBX4G9AYGOufC7j+b2SXAH4Fs4Fnn3K+9rF+75OmvvNyxYsch5qzfx5pdR1i45QCrdwV3tv28U1rRu11jBnVtzqBuzXVDjAyRqF3yZcBVwNhIDcwsG3gKuBDYCnxsZpOdcyti7FvSQFaW0Te3yeceY1Fe7th28Birdx1m1a7DrN55mO0FRew9UsyGvUeJ5m987ZwsurdqSIdm9TilTUN6tm5Ej9YN6d6qIfVq6xIgiU5MgemcWwlUt9syEFjrnFsfavsSMBpQYEpYWVlGx+b16di8Phf0bhO2jXOOwpIyDheVcrjoOA5oWCeHBnVyaFA7mxwdb5Q4CGL/IxfYUun1VmBQAP1KGjOzinCsk0PbJroJrwSj2sA0s+lA2zCL7nPOveZ3QWY2BhgD0KlTJ79XLyJSY9UGpnNuZIx9bAM6VnrdITQvUn/jgHFQcdInxr5FRHwTxIGej4GeZtbVzGoD1wGTA+hXRMRXMQWmmV1pZluBs4E3zGxqaH57M5sC4JwrBe4EpgIrgX8555bHVraISPBiPUs+CZgUZv52YFSl11OAKbH0JSKSaLr2QkTEIwWmiIhHCkwREY8UmCIiHikwRUQ8UmCKiHikwBQR8UiBKSLikQJTRMQjBaaIiEcKTBERjxSYIiIeKTBFRDxSYIqIeKTAFBHxSIEpIuKRAlNExCMFpoiIRwpMERGPFJgiIh4pMEVEPFJgioh4pMAUEfFIgSki4pECU0TEo5gC08yuNbPlZlZuZnlVtNtoZkvNbJGZ5cfSp4hIouTE+P5lwFXAWA9tRzjn9sbYn4hIwsQUmM65lQBm5k81IiJJLKhjmA5428zmm9mYgPoUEfFVtVuYZjYdaBtm0X3Oudc89nOOc26bmbUGppnZJ865WRH6GwOMAejUqZPH1YuIxF+1gemcGxlrJ865baHvu81sEjAQCBuYzrlxwDiAvLw8F2vfIiJ+ifsuuZk1MLNGn04DF1FxskhEJKXEelnRlWa2FTgbeMPMpobmtzezKaFmbYAPzGwxMA94wzn3Viz9iogkQqxnyScBk8LQ+9SxAAAGC0lEQVTM3w6MCk2vB86IpR8RkWSgkT4iIh4pMEVEPFJgioh4pMAUEfFIgSki4pECU0TEIwWmiIhHCkwREY8UmCIiHikwRUQ8UmCKiHikwBQR8UiBKSLikQJTRMQjBaaIiEcKTBERjxSYIiIeKTBFRDxSYIqIeKTAFBHxSIEpIuKRAlNExCMFpoiIRwpMERGPFJgiIh7FFJhm9jsz+8TMlpjZJDNrGqHdJWa2yszWmtk9sfQpIpIosW5hTgP6OudOB1YD957YwMyygaeAS4E+wPVm1ifGfkVEAhdTYDrn3nbOlYZezgE6hGk2EFjrnFvvnCsBXgJGx9KviEgi+HkM8xbgzTDzc4EtlV5vDc0TEUkpOdU1MLPpQNswi+5zzr0WanMfUAq8EGtBZjYGGBN6WWxmy2Jdpw9aAnsTXUSIaglPtYSnWsI7tSZvqjYwnXMjq1puZjcBlwEXOOdcmCbbgI6VXncIzYvU3zhgXGjd+c65vOpqjLdkqQNUSySqJTzVEp6Z5dfkfbGeJb8E+DFwuXOuMEKzj4GeZtbVzGoD1wGTY+lXRCQRYj2G+STQCJhmZovM7GkAM2tvZlMAQieF7gSmAiuBfznnlsfYr4hI4KrdJa+Kc65HhPnbgVGVXk8BptSgi3E1LM1vyVIHqJZIVEt4qiW8GtVi4Q87iojIiTQ0UkTEo6QKzGQZamlm15rZcjMrN7OIZ/XMbKOZLQ0dv63RWTcfa4n78FMza25m08xsTeh7swjtykKfySIz8/UEX3U/p5nVMbOXQ8vnmlkXP/uPspabzGxPpc/i1jjV8Tcz2x3pEjyr8KdQnUvMbEA86vBYy3AzK6j0mTwYpzo6mtl7ZrYi9P/nrjBtov9cnHNJ8wVcBOSEpn8D/CZMm2xgHdANqA0sBvr4XEdvKq7TmgHkVdFuI9Ayzp9JtbUE8ZmE+vktcE9o+p5w/z6hZUfi9FlU+3MCtwNPh6avA15OYC03AU/G8/cj1M+5wABgWYTlo6gYVGLAYGBuAmsZDrwewGfSDhgQmm5ExdDtE/99ov5ckmoL0yXJUEvn3Ern3Co/11lTHmsJavjpaGB8aHo8cEUc+qiKl5+zco2vABeYmSWolkA452YB+6toMhp43lWYAzQ1s3YJqiUQzrkdzrkFoenDVFyhc+IIw6g/l6QKzBOkwlBLB7xtZvNDI5QSJajPpI1zbkdoeifQJkK7umaWb2ZzzMzPUPXyc37WJvTHtwBo4WMN0dQCcHVod+8VM+sYZnkQkun/DMDZZrbYzN40s9Pi3VnosEx/YO4Ji6L+XGK6rKgmgh5qGUsdHpzjnNtmZq2puBb1k9Bf2ETU4ouqaqn8wjnnzCzSJRadQ59LN+BdM1vqnFvnd60p4L/ABOdcsZn9DxVbvucnuKZEW0DF78cRMxsFvAr0jFdnZtYQ+A9wt3PuUKzrCzwwXcBDLWtah8d1bAt9321mk6jYTYs6MH2oxZfPpLpazGyXmbVzzu0I7brsjrCOTz+X9WY2g4q/7n4Eppef89M2W80sB2gC7POh76hrcc5V7vdZKo4BJ4Jvvx+xqhxazrkpZvYXM2vpnPN9jLmZ1aIiLF9wzk0M0yTqzyWpdskthYZamlkDM2v06TQVJ6wSdaOQoD6TycA3QtPfAE7a+jWzZmZWJzTdEhgKrPCpfy8/Z+UarwHejfCHN+61nHA87HIqjqMlwmTg66GzwoOBgkqHVgJlZm0/PaZsZgOpyCDf/6CF+ngOWOmc+32EZtF/LvE+WxXlma21VBxTWBT6+vRsZ3tgyglnt1ZTsdVyXxzquJKK4xnFwC5g6ol1UHF2dHHoa3k86vBaSxCfSaiPFsA7wBpgOtA8ND8PeDY0PQRYGvpclgLf9LmGk35O4GEq/sgC1AX+Hfpdmgd0i+Pva3W1PBr63VgMvAf0ilMdE4AdwPHQ78o3gW8B3wotNypu4r0u9G8S8cqPAGq5s9JnMgcYEqc6zqHiHMOSSnkyKtbPRSN9REQ8SqpdchGRZKbAFBHxSIEpIuKRAlNExCMFpoiIRwpMERGPFJgiIh4pMEVEPPo/KtcyO3+TPlIAAAAASUVORK5CYII=",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x108966940>"
      ]
     },
     "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": "iVBORw0KGgoAAAANSUhEUgAAAUwAAAEzCAYAAABaGjpLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X1U3Od5J/zvNViyXiJkwDaoaV0HPzkbDejsbg0D7ub0dLe7W1mABCj7bLrnWIJBQYqb2q1lCQYY7TkCLCGfuE9eHFszw2C3f/RVQDSgl6dttqcrmGFQcpKIQU+2jtPn1El2s232xGnXlQRz7R/M/csPjKQRMzBv3885czTAMHNr9NM198t1X7eoKoiI6P4cmW4AEVGuYMAkIkoSAyYRUZIYMImIksSASUSUJAZMIqIkpSVgikhQRH4kInN3+fmvishPROSbidupdLwuEdFGeihNz/MWgC8D+L17POa/qmpjml6PiGjDpaWHqap/BeDH6XguIqJstZFzmM+IyLdE5LKIVG3g6xIRpUW6huT38w0Av6iq/yAi+wCMA/j4ag8UkU4AnQCwffv2pz/xiU9sUBOJqFB8/etf/ztVfexBf0/StZdcRJ4EMKGq1Uk89m8A1Kjq393rcTU1NXr9+vW0tI+IyBCRr6tqzYP+3oYMyUWkQkQkcd+VeN2/34jXJiJKl7QMyUXkDwD8KoBHReQ9AP8ZwCYAUNU3AXwKwGdFZAHABwA+rSyTREQ5Ji0BU1V/4z4//zKW0o6IiHIWd/oQESWJAZOIKEkMmERESWLAJCJKEgMmEVGSGDCJiJLEgElElCQGTCKiJDFgEhEliQGTiChJDJhEREliwCQiShIDJhFRkhgwiYiSxIBJRJQkBkwioiQxYBIRJYkBk4goSQyYRERJYsAkIkoSAyYRUZIYMImIksSASUSUJAZMIqIkMWASESWJAZOIKEkMmERESWLAJCJKEgMmEVGSGDCJiJLEgElElCQGTCKiJDFgEhEliQGTiChJDJhEREliwCQiSlJaAqaIBEXkRyIyd5efi4h8UUTeEZFvi8gvpeN1KX+oKmZmZqCq1teRSASRSATxeBwzMzNYXFzE8PAwpqenEQ6Hl/0sHo8jHA5jeHgYi4uLy56LKG1UNeUbgF8B8EsA5u7y830ALgMQAPUAZpJ53qefflopP8XjcY1EIhqPx1VVdXp6WouLi/XatWsaCAR0enpaS0tLdefOnRoIBHTXrl3q8XgUgG7evFmLi4t1586d6vf7taysTAOBgBYXFysA3bdvnz7++OMaCAR0cXFRp6enNRAI6O3bt7Wnp0cXFhYy/LenTANwXdcQ69LSw1TVvwLw43s85ACA30u0NQLgERHZlY7Xpuynqlbvb2FhAYFAAIFAAA0NDVZPMBQK4ac//SkGBwfR2dmJUCiEoaEhFBUVwel0YnR0FB/72Mewfft2OBwOPP/88ygqKrJ6kU6nE5cvX8azzz6LK1euoKamBn19ffB6vdi7dy8+85nPoLm5Ga+88grcbjf8fj8CgcCHeqO6oqdLtMxaouxqNwBP4u49zAkAn7R9/RcAau73nOxh5h7Tc1xcXLT+NL0/EdHDhw+riOiWLVu0uLhYw+GwhsNhLSkp0cOHD+ujjz5q9RBNz3BxcVHD4bCWlpaqx+PR0tJSnZ6e1kgkogsLC8seY57n8ccf156eHi0vL1efz6cej0dv3bqlhw8f1i1btigABaAej0fLysp0ampKA4GAXrt2TYuLi9Xn82k4HLZ6wJRfsMYeZtYFTACdAK4DuP7EE0+s09tF6RSPx63ANzU1pTt37lSfz2cNlSsqKtTn86nf79dr166px+PRkpIS9fv9Go/HreG4CVrl5eUaCAQ0HA7rrl27NBKJaDgc1rKyMitQmkAWiUTu+hgTtFc+T0lJiT733HO6Y8cOPX/+vJaVlanH41ER0UOHDikAffjhh7W0tFTD4fCy16P8kO0B8zyA37B9/R0Au+73nOxh5oZwOKw7d+60eoAAtLu7e1nPzfQATUAzAdYE27KyMutrE6Dudt9urY+x94AjkYhOTU1pcXGxvvnmm/rwww9b86MmgHs8Hp2amlr296Hcle0BswHLF32iyTwnA2b2MIFmYWFB/X6/Tk9PW8HGLNB4PB5rSDs1NWX16CoqKqyFHBMY7T3DuwW6TPz9zN/F9HArKiqsaYTt27frww8/rAD02Wef1ampKfY8c1RGAyaAPwDwQwB3ALwHoAPAMQDHEj8XAK8D+C6AG8nMXyoDZlYxvcBDhw5ZwcOsXofD4WVDadODND06Myy3D2+zIUiuZrXeqPmQ8Pl8+sgjj+gzzzyjAHTLli3q8/l0enp62d+Zsl/Ge5jrcWPAzA7xeFz9fr8WFxdrSUmJHjp0SEtLS5cNT+2B0QTSSCRi/X42BscHZf4e165d04cffli3b99upTft2LHDWizikD37rTVgytLvZqeamhq9fv16pptRUFSX0moAwOVyYXZ2FqqK5uZmdHR0oKGhASICEQEAtLS0YGBgAG63GwAQjUZRW1uL2dlZuFwu63H5xLxH9v87sVgML7/8Mj744APcvn0bHo8H+/fvR11dXV6+B7lORL6uqjUP/ItribIbdWMPc+OY9B+z+FFWVqY9PT3WUNoMuXt6epYt0JhepelNFirzXnzkIx/RrVu36kc+8hH2OLMY1tjDfCj9sZtyjarC6/Xi7Nmz6Orqgqqis7MTgUAAR44cgcvlQl1dHQDgxIkTWFxcBACICNxuN6qrq+FyuTL5V8g4815UVVVBVTE3N4fjx4/jhRdewD/90z8tDedEsHv3bjgcDvY8c9VaouxG3djDXB8r5xQjkYjVe5yamtLS0lIrpcbee7TnW+b6fOR6Mz1O01P3+XwKQLdu3Wol7FPmgIs+lKzVErxNIFxYWLB2yJjvMUCujf2DaXp6Wnfu3Knd3d26c+dO9Xg8urCwkBeLYblorQGT5d0KSDwet6r53LlzB7FYDC0tLZidnYWIoKWlBadOncLw8DAGBwetYWNrayui0Wimm59zRMR6D+vr63H16lXs378ft2/fxpkzZ+D1etHS0sL3NocwYBaQkZERdHZ24tKlS3jooaXp69HRUbhcLrhcLgwMDFjB0u12Q0TgcrkwNjZW8HOUqTLBs76+Hi+++CJEBE8++SROnz6NeDxulalb6vxQtmJaUYFQXaovGYvFsHv3bkxOTiIYDC4LjqqKaDSat+lA2SIejyMYDAIAurq6AADHjh2zFtn6+/vhcLAvs56YVkTLmDQhk85itiKahYjS0lJrrjIQCHAebYNFIhEryd/v92tpaak+99xzVgUl/nusLzCtiABYvcQbN26gs7MTqoo9e/agpqYG/f39aGtrg9PpRCwWQ3t7OyorK9HX14fq6mordYjWn8vlwvj4OFwuF2ZmZiAi+OhHPwoA+NKXvoTKykrr34S9/Syylii7UTf2MB+cWQH3+XxaXFys3d3dVk/G/Gkvd5Yv2xZzmX3PeiAQ0PPnz+u2bdu0tLS04DcErBdwlZzsqqur8dprr2FkZAQDAwNob2/HwMAA+vr6AMBayLGv5FJmmH+DoqIidHR0oKioCB988AGOHj1qnVUUiUS4IJQFOCTPE2pbsJmYmICqwul0YmxszNr7bd+VwwCZvdrb263g2NDQgNu3b2Pz5s347Gc/ywWhDOM7n8NUf1YEIhqNLsupbGpqQlNTE+bn59HY2GjNk7E3mf0cDgf27NkDr9eLY8eOoaioCI2NjVbuJnuamcOAmaNUFcFgEAcOHEAwGEQ8Hl+WUxkKhXD27Fns3r07002lNTCLQk1NTdi0aRM++clPYvPmzXj99dcRDAYZNDOEeZg5amZmBs3NzXC73fjyl78MALhy5QocDgdcLpfV4xwdHbUS0NmzzD1m9BCPx9HU1ITOzk74fD6EQiHr35r/rg9urXmY7GHmILNiNzY2hsbGRut7sVgMzc3NCAaDqK2txdjYGOrq6jgMz2H2HUITExOorKy0/q3NMcW0cbjok4NM79Gsfl+5csXKqxQR5lXmIbNw5/V6ce7cOezevRu3b9/Gt7/9bczNzcHtdnMxaCOsJRdpo27Mw1yy8pwZU9DXXtzXlGFjXmX+sp94GQgEdMeOHbp9+3buDloDsLxb/lq5jc5+oJg52bCnp4dVvQuEuR78fr92d3fr1q1btaSkhFtcH8BaAyaH5FnO/EP19/fj5MmTAJYKNZjJ/pmZGSsh3ZwhQ/nNrKCrLlXK/+IXv4jvfve7OHHiBKqqqlBfX5/pJuYtrpJnuZmZGbS0tODChQuYn59HPB7HqVOnMD4+DgBcCS9gmlhBV1U0Njbizp07+PznP8896ElY6yo5A2YWU9Vlq6BmocfpdFoBMp9PZ6TkqC6V7guFQnjzzTfhcDgwNDRkle2jD2NaUR6KRqNobW21UkvM0Ht+fh6tra2YnZ1lT4IgInA4HBgZGcGrr76Ks2fP4sSJE+jt7UU8Hs908/IKA2aWMnOXFy5csO47nU4MDAygra2NVdBpGTOvaVLLbt26hbNnz3IrZbqtZaVoo26FuEpuUkamp6etgr/2wr9lZWUs+UWrMidVlpeXq8/n00OHDuljjz3G1fNVgKvk+cGcu3P+/Hmr4G91dTVqamqsnkJtbW2GW0nZKBqNoq+vDx0dHaiursapU6fwmc98hhsZ0ogBM0toYsWzra0NAOB0OtHa2goAcLvdiEaj6O7uBgDs2bOHFz99iDnIrq+vD01NTRgbG7NW0E0PifPdqeEcZpawl2erqqoCACv3cmZmxqpzOTExwblLWpWpeTo+Pm4tBh48eBDz8/NoaGhglaM0YMDMAubTf3R0FADQ2NiIpqYma/8wAOts6/r6evYS6K7sNU9dLhdGR0ehqrh9+zZeeuklDA8PM2imgAEzC6xMHzK1LNva2jA5OcnhN62JvWDHCy+8ABHBSy+9xJ5mCjiHmUFm3tKc6FhTU4OZmRnEYjF4vV7OVVLKTLpRTc1SjvYXv/hFHD9+nFso14g9zAyamZlBQ0MD3nrrLXi9Xrz11ltobGxEV1cX3G43V8MpbWZnZxEMBvG5z30OqoobN26wl7kGDJhZYPfu3VYK0cTEBIaGhhAMBjE7O5vpplGOM4uJqoqBgQFrbry7uxvRaDTTzcs5HJJngBmKu1wuTE5OQlWt42/b29sBgDt5KC1cLpeVXuT1enHhwgV8/vOfx+7du5lqtAbsYWaA+dQ3n/C1tbVwu93o7e3FyMiItQDEC5lSZRYS6+rqrCOXvV4vbt68iZaWFi4APSD2MDOgtrYW/f39UFW0traiv78fwWAQHR0d1s4e9i4pnUzgVF06C6q2thaqipMnT0JV0dHRwQ/oJKSlhykie0XkOyLyjoh0r/LzNhH5nyLyzcTtSDpeN1fNzs7C6/VCRDA2Noa2tja43W4MDw/j+vXrrEBE68ZeFtDpdOLWrVt46aWXeJhaklIOmCJSBOB1AM8CcAL4DRFxrvLQP1LVf5G4BVJ93VykibqF8XgcFy5cALA0x3T9+nUEg0EMDg6yZ0nrLhqNorm5GZOTkygqKsLi4iKH5UlKRw/TBeAdVX1XVW8D+EMAB9LwvHlnZmYGe/fuRWNjI+bn5605pNraWoyPj7PgK20IM2ceCATwm7/5mygqKkIsFmPQTEI65jA/CuBvbV+/B2C1bOuDIvIrAP4bgN9R1b9d5TF5r6ioCEePHrWKbPT29gIAgyVtGJOT+corr8DpdOIrX/kKurq6rL3ovA7vYS014ew3AJ8CELB9/RyAL694TBmAhxP3jwL42j2erxPAdQDXn3jiiTVWu8suK49HLS8v10AgoAsLC9rT06Pl5eWscUkbxn49hsNhnZ6etk4jLZTrEJk6ZhfAMwCu2r72APDc4/FFAH6SzHPnSwHhSCSiu3bt0nA4rOFw2Lo4zXniLPBKmWC/Lk3QLJSjmtcaMNMxJJ8F8HER+RiA7wP4NID/ZH+AiOxS1R8mvtwP4GYaXjcnmDfaVCIyB5mZr8fGxrgqThlhT2pvamqCquJ73/se+vv74XAwRXs1Kb8rqroA4HMArmIpEP6xqsZE5LSI7E887AURiYnItwC8AKAt1dfNFfc6yOzgwYNMUKeMsSe1T0xM4OjRozh79izcbjcPT7sLHrO7jlSXjsnVxPYzU3lo5fcYMCmTNLFV9+mnn8aBAwdw6dIl9PT0YGBgIG+vTR6zm4VM79KeQgTAqoTN3iVlA7NV9+2338Y3vvEN7Nu3D1/5yleYzL4K9jDXkelhxuNxTExMIBgMYnx83Po5e5eUDewjoVgshq6uLiwsLODq1at5WzOTPcwsZHqQ+/fvx/nz5zE4OGgdSmV+TpQNYrEYWltbUVVVhXPnzuHy5csAwGT2FRgw14n51K6trbUOL6uqqoKq4s6dO7wQKWuY43nNnKWpZtTa2sqamSswYK4Ts193ZGQEdXV1cDgcaG1txc2bN7Fp0yb2LilrmPSiqqoq6357e7t1gBo/3H+GAXMdmIvMpBCZeczTp0/zYDPKOmbqyN6jNB/ojY2NXPyxYT3MdWBWx0dHRzE+Pm4lBgPgwWaUlexJ7C0tLRgdHUUsFkM8Hsfc3BwXKBMYMNPM9C7t5dtUFUNDQ3A6nSzfRlnJXmB4dHQUc3NzOHnyJI4ePWqdMvnMM89kupkZxyF5mq2WexmNRuH1euFwOPgpTVnNDM+7u7uxuLiIH/zgB3j//fcxPz+f6aZlBeZhptnK3Mvh4WEMDg4CWCrhxj26lO00Ueg6FApheHgYR44cQWNjI+rr6/PmA3+teZgckqeZPfdSVXHs2DGr1iDnLynbmW2SIoKRkRGrZmZTUxMmJibyNpE9WQyYaWIuNJfLBZfLhYmJCSsdo6GhAfPz86itrc1wK4nuzWyTNAuWtbW1GBkZYWpRAseHaWI/OldEUF9fD4fDgYMHD+LmzZs4deoUZmdnM91Monsyq+WmitHs7Cz6+vpw7NgxLliCATNtXC6Xlegbj8etw85GR0fR3t6OsbExXnCU9cypktFoFPF4HKqK/v5+DA8Ps6cJBsy0isViaGhowMjICBobG9HY2IhYLMYybpRTzGhpZGTE2l/e0dGBnp4eBIPBgg6aXCVPk5mZGTQ0NODOnTu4cuUKRARzc3Po6urC5ORkwU+WU+4w8/G1tbWYnZ21ktl//dd/HVeuXMFXv/rVnF+8ZLWiDDMLPVevXrV6k1VVVexVUs6xj4jMRoz+/n5cuXIFHR0dBb14yYCZRmb+Z2RkBM3NzRAR7hunnBWNRtHY2IimpiZrWP7GG28UdAUjphWliZn36e/vt0plZfN0B9H9uFwuhEIhxGIxAMD58+exuLiY4VZlFgNmiuzzPWNjY6itrUVVVRXm5ubQ2NjIXiblJHsCu9frxcDAAC5evIj5+fmCzvbgkDxFpmc5Oztr1b00F9mxY8cQCoUK+gKj3GSuawBWmULmE7OHmbKV+Zezs7Oora21LrL9+/dz4Ydyjklgd7lcqKurQ3V1NWpqlhaVuehDa2YvvmoWe0ZGRuB0OpmsTjlrtWOho9EoTp48ib6+voI9t5w9zBStrH/Z39+PkydPWnOX7F1SLjMr5QAwNDSEW7du4cyZM6isrMSRI0cy3LqNxx5miuz1Lw8ePIjq6mpMTEzg7NmzBT10ofxg8otDoRCcTid+93d/F8XFxaiqqsp00zKCATNFKw+NApaGM4U+OU75wQzN5+fn0draiurqarz22msFm/XBgJkik6xutpA1NjZCVTl/SXnDfgxvLBZDb29vwRbiYMBMgamuPjMzg5aWFquMfywWK8iLifKTyfpwOp3wer3o6OhAX19fQe744aJPCuzFVk3SOgD8zu/8DjZt2oRLly4V7NCFcp9ZGVdVeL1eXLhwAQMDAzh8+DAqKysLco6ePcwUmBxMAFZlF6fTic2bN+Po0aMFeUFR/rAnr4+NjVkbMt5++214vd6CnKNnDzMFJgfTvoe8v78fnZ2dCAaDOHDgAHuYlLPsyeumctHY2JiVwG4KDBdS6hzrYaYoHo9jZGQEu3fvxvz8PLq6ugAA586dg9vtLqiLifKXma8395999lkUFRXl7LQTT43MkNnZWStITkxMYHJyEgBYYZ3ygn0e0ySwnz17FqqKs2fPFlwmCANmilaeEGlSjIjygZnHvHDhAoaGhuB0OgEAmzZtgsNReEsghfc3TjP7CZFmP7k5PZIo15mFzVgshr6+PjgcDtTX12NoaKggU4vYw0yBfV6ntrYW/f39+MQnPoHTp09zhZzyglnYNDUxa2trMTMzU7DFZdISMEVkL4AvACgCEFDVsyt+/jCA3wPwNIC/B/AfVfVv0vHambSyMEFXVxc++OADbNmyBXv27MnJyXCilVwuF8bHx63jd801X4jFZVIOmCJSBOB1AP8OwHsAZkXkoqrO2x7WAeB/qer/JSKfBjAE4D+m+tqZZi/h39bWBgA4efIkzp07V3CfvJS/7GeV19bWWtd8IY6i0jGH6QLwjqq+q6q3AfwhgAMrHnMAwNuJ+38K4NckDz6aRAQOh8MqtFFVVYWJiYmCreRC+ce+/dfUegWA7u7ugpu/BNITMD8K4G9tX7+X+N6qj1HVBQA/AVCWhtfOODMpPjc3h+bmZkxMTHDRh/LGakdVmJoJFy9eLLhCwlm36CMinQA6AeCJJ57IcGuSY1YQOzo6MDw8XPBnN1P+MLt9zPVstgJ3dnbi7NmzeOqpp9DR0ZHJJm6odPQwvw/gF2xf/3zie6s+RkQeArATS4s/H6KqPlWtUdWaxx57LA3NW1+m9NXg4CBOnz7Ns5spr5j5S5MuZ+piVlZWwufzob29PdNN3FDp6GHOAvi4iHwMS4Hx0wD+04rHXARwGEAYwKcAfE2zeU/mA1i5gsizmynf2OthmgDZ29uLwcFBrpI/KFVdEJHPAbiKpbSioKrGROQ0gOuqehHAMIDfF5F3APwYS0E159nPJDd/TkxMAABTiihv2OthigjcbjcAoK+vD9XV1QV1rbP4RgpM4WBTqch8ApvhOPeTUz6YmZlBQ0MDAFgdAlOlKFevcRbfyAD7meQmaALASy+9hIceeihnK7kQ2Zl6CQYT12lN7GeSm2rUu3fvxqZNmzA0NMTkdcobZvEHABPXae1M2oXZb+twOJi8TnlhZdJ6MBiEqmJ+fh59fX0FeRAaA2aamOIbtbW1iMViaGhosApzEOWi1ZLWR0ZG0NfXh71796K3t7fg0uc4JE+Ruaj6+/vh9XoBAF1dXVhYWMhwy4hSc7ekdTNfPzg4WHDTTgyYKbJfVNXV1aipqYGqQlUL7mKi/GLmLYPBoJUFYv4cHx/P2RXylJj/3Nl4e/rppzXXRCIRLSsr07KyMo1EIpluDlFKIpGIVlRUaCAQ0MXFRQ0EAlpeXq6BQEDj8Ximm7dmWMoRf+CYxB5mitSWvD47O4uamhqrlD97mJTLTJAYHR21MkIKOWkd4KJPyswc5sjICBoaGqxJ8fn5eczMzBTcKiLlj2g0itbWVmv/uFkxL9Rq6wDnMFNmL++mqnA6nXC73Xj55ZcRj8dx9epV1NfXZ7qZRA9s5fy82k6OLMSkdYA9zJTZzzw5d+4cHA4HgsEgnn/+eWzatCnTzSNKG3PCwNmzZwsyaR1gDzMt7BWLAGB8fBw1NTWorKwsyGEL5TYzL6+qaG1tXVYrwel0oru7G9XV1QU5cmIPMw3sZ56YifKRkRGcPHmy4BJ7KffZE9bHxsbQ3t6OgYEB9Pb2IhQKFfS8PANmmtgXfxobG/Hiiy/izp07mW4W0QMzc5f20VFVVRUGBgYQDAZx7ty5glsdNxgw08RsjWxra8PQ0BC2bNmC1157DQAK+hOZco8p2zY7O2t1AlpaWiAiGB8fh9vtLsgFH4BzmGkzOzsLr9eLqqoqVFVVYXJyEgDQ2tqKsbGxgv1Eptxjzy22b400VdYL+VpmDzNNzDDGpF6YuczTp08X7Ioi5SZTMNgUjzEJ64ODg+jr6yvoeXkGzHUyPz+PpqYmHD9+3CqLRZRLYrHYsrJuVVVVBZuwbnBIniZm0Wd0dBSTk5OoqanBu+++iy996Uvo6urCnj17CnooQ7mjrq4Ok5OTqK2thYigt7cX7777LoLBIMbHxwt2/hJgwEwb+8qiiGBmZgbBYNBa+OGwnHKFWfSx71wLBAI4cuRIwV/HDJhpYi4yYHnRglgshpMnTwIAOjo6CvrTmXKDPXG9qakJqopjx44hGAxi//79BT1S4qmR68B+mmRXVxf+8R//EVu3bsXly5cL+mKj7LZyh48pGGxiRC6fErnSWk+N5KLPOliZk7lt2za8+uqrBT1ZTtkvGo2iubkZsVgMFy5cALA0n+lwOHDw4EGrbkIhY8BcByYnc3Z21srJNOeVZ3OPngqXmUYyVdXn5+fR0tKCYDCIeDyO0dFRfuCDAXNdmAUgYClxHQBGRkawb98+HoxGWcnUvqyqqsL4+Li1f/zkyZNoampi7zKBiz7ryAzNVdU6GC0Wi+XNPBDlD/upp7Ozs4hGo2hvb4fT6SzYM8hXw4C5DlaeJHnhwgUMDQ1hcXERJ0+eRFVVVUGWxqLsZaaRgKVTTwFgYmIC8/Pz8Hq9zCNOYMBcB6tVqvZ6vTh9+jQWFxc5j0lZx37NVlVVWd83dTA5f7mEc5jrwKRfOBwO60K7cOECVBXxeDzDrSNazl5sw+wTr6urg8vlwsDAANrb2zmFlMCAuc7sB0l5PB6ICCYmJhg4KeNU1TrYzF7LtaGhAcFgENFo1Mr2oCVMXF9n8XgcIyMjaGtrw+zsLC5evIihoSH4fD50dHRkunlUwEygvHDhAubn57F7926ICGKxGLxeL8bGxqzTBPKth7nWxHXOYa4ze51MEUF/fz+eeuopOJ1OqGreXYiUO+wlCbu6uqCqOHfuHNxuN/bs2ZOXgTJVDJjrzH5RNjc3Y2BgALt378bevXtx5coVPPPMM5luIhUgM29p5tgnJiYwNzeHvr4+68OdPoxzmOvMDGkAWKfvTUxM4P333y/4A6Uoc0zq28zMDKLRKOrq6tDR0YGxsTHMzc2dhhdJAAATzklEQVShubm5oAsF3w0D5gYwCz/V1dUYHx9HY2Mjtm/fjjfffJPFhSkj7Jsq7EWCY7EYU4nuxewhzcbb008/rfkgHo9rJBLReDxufR0Oh9Xv92tZWZmGw+EMt5AKTSQS0V27dmk4HNZAIKAVFRUaCAS0rKxMA4GAda3mKwDXdQ0xKaUepoiUisifichfJ/4sucvjFkXkm4nbxVReMxfZy2Lpirmj27dv48aNG+xl0oZyuVwYHR21igSPjY3B6XQCAOcw7yHVIXk3gL9Q1Y8D+IvE16v5QFX/ReK2P8XXzGnmgKlgMIju7m4sLi7i5ZdfZlEO2jDmQxsAmpqarOIa9fX1mJyc5BbIe0g1YB4A8Hbi/tsAmlN8voJRVVWFUCiEF198EQ899BBisRh7mbSuVBWRSATDw8NW3ctQKLRs8ZGFYe4t1YBZrqo/TNz/7wDK7/K4LSJyXUQiIlLQQdUcMFVfXw+Hw4GRkREcO3as4I8vpfUXjUbR2NiIrq4uuN1u9Pb2Yn5+HiKC1tZWXn9JuG8epoj8OYCKVX7Ua/9CVVVE7tZF+kVV/b6IVAL4mojcUNXv3uX1OgF0AsATTzxxv+blHPsBU/F43KrMXllZiXg8zmR2SjszBK+trcXExASApTnMyspK9PX1YXR01CrtRvexlpUicwPwHQC7Evd3AfhOEr/zFoBPJfP8+bJKvppIJKJlZWXWqmRpaalu27ZNp6enM900yjORSMRaBV9cXLQyNuzZGhUVFRqJRDLd1A2DNa6Sp7rT5yKAwwDOJv786soHJFbO/7eq3hKRRwH8KwDnUnzdnOdyuZZ92r/77rs4c+YMJiYmUF9fz14mpYX5j26OngAAr9eL0dFR6xrzer3Mu0xSSsU3RKQMwB8DeALA/w/g/1bVH4tIDYBjqnpERH4ZwHkAcSzNmf4/qjqczPPnQ/GNZGhiMj4UCuHNN9/Eq6++CrfbzaBJKVutwIbD4YDq0hG6oVDIKkNYSNdbRopvqOrfA/i1Vb5/HcCRxP1pAHtSeZ18NzMzg6amJpw5cwa3bt2yzjFn0KRUmR09wM8qqU9OTlo/N3PqlBxujcwiIoItW7agsbERJ06cYG4mpcxUyzJ1WE0KkcvlYs7lGrBaURYwqUa1tbUQEXR1dWFxcZEHptGaaKIwMPCzHqbL5YLD4bCG6GNjYwyWa8CAmQXsqUZVVVW4ePEiQqEQTpw4AVVFR0cHgyYlzeRbAsDQ0NCyeqw1NTVMIUoBA2YWMVWN+vv74ff7cevWLRw/fhzV1dU8ZZKSYlbFQ6GQVVrQHMTX3NwMt9uNYDCI6upq9jDXgHOYWcQUG25vb8fExAS+8IUvYNOmTQiFQjwDiO5LVREMBtHS0oL5+Xm4XC5Eo1FrznJgYADDw8NMIUoBe5hZxL5iaf7s7OzEmTNnoKrYv38/5zTprqLRKHp7e7F371709i5txLOvjLvdblRXVxdcClE6sYeZpUxF7MrKSnR1deH111+3qhylkjtL+cvlcmFwcBCXL19GR0cH2tralq2MAyyukSoGzCxlhlBerxdPPfUUgKXe5smTJ5luRMvYV8XdbjdeeeUVDA8P46233kJdXR0cDgeLa6QJA2aWEhG43W6Mj4+jqqoKmzZtwlNPPWUdg8peZmEzQTIejyMYDFrHTABLQXNwcNCqgGXmxjlvmQZr2YC+Ubd8Lr7xIBYXFzUQCOidO3fU4/FoSUmJ+nw+DYfDeX+UAH1YPB5fdqxERUWF9vT0aHl5uXW8hCmswWtkdchQ8Q3aAGa3BgD4fD7cvn0bv/3bv42tW7dyt0YBikaj1kFlTqfTOmpCVdHb27vsiInW1lYmqacRh+Q5wJ5uFAqF8MILL2Dz5s04c+YMvv3tbyMcDnOIXiBMYBwbG0NVVRUOHjyI+fl57N+/H+fPn8fg4CAAoKWlBQA4FE+3tXRLN+rGIfmHRSIRLS8v156eHvX5fApAi4uLC6qWYSGz17ZcWFiw/pyenla/36+Li4sfOqWUPgyZODWSNp5JHQkGgxARFBcX4/nnn8eNGzeY3J7HNLHIU1tba9W2fOutt+D1enH9+nU4HA6cOnUKs7Ozy04ppfTiHGaOMavn1dXVVrGO48eP4/3338f3vvc9DAwM8D9KHopGo2hubsbAwADa29utf/+qqiprJw+H3+uPATMHrSzWcenSJYRCIbzxxht48skn4XA40N7eDoeDA4hcZ3qWqj+rmm4WdcwHY2Njo1Wpn9YX/0flMFOso6ioCAcOHMDi4iJ+67d+C0eOHMHIyEimm0cpMIFyZmYGjY2NVvWhsbExALDyLpWLfRuKATOH2YdhdXV1eO2117BlyxZs27YN8Xgc4XAYkUiE/6lyiD1QNjc3Y25uDqFQCENDQ+jt7UUsFrN2gfX19UFEmFq2kdayUrRRN66SPxiTrOzz+bS4uFi3bNmixcXFGg6HM900SoI9IX16enpZMrrZvFBaWmolo3MlfO3AxHUSEWseS0TgcDggIpibmwPAwgvZzp6QLiIIBoPo6OiwKg85nU7r349n8WQGA2Yeqqurw+XLlxGLxQAAJ06cwO3bt/GFL3yB1duzjOrPCmeYKRYAVvpQW1sbKisr0dfXh7GxMUxOTnIlPIMYMPOQ6V2eOnUKo6Oj+OxnP4szZ87gxRdfhNPpxM2bN7mKnmGqahX3NQs6k5OTEBG0tLSgv7/fOke8ra0NAKxzeSiD1jKO36gb5zDXzj7Htbi4qD09PVpWVqaHDh1Sh8OhPT09nP/KoEgkort27dLp6WkNBAI6NTWl4XB42e4de4GNXbt2cTdXGmGNc5gZD4r3ujFgpo9ZUHjsscf02Wef1R07dui1a9esBQXaOPF43AqU09PTumvXLiso2oOj+dBbXFzkAk+arTVgckheIMwOIQA4fvw4fvrTn8Ln8+H3f//38e6773KH0AZQ2zC8qakJABAKhdDf34+2tjZUV1ejpqYGAKxdXCuPLKHMYsAsICZoOp1OzM3NQVXxJ3/yJ3jjjTcQj8dRWVlpnVDJ4Jk6EyBra2sxOzuLeDyOpqYmhEIhhEIha1HOzFW63W5Eo1F4vV6e6pit1tIt3agbh+Trx171pqenRwGoiGhxcbH6/X4Wnk2RmQIpKyuz5iI9Ho+WlpZa729FRYX6/X4NBALL8i05/F5/4BwmPYiVi0I+n0+7u7v1/PnzWlxcrNu2beMc5xqY9zUcDmtpaakWFxdb85WmLF9FRYWGw2ErkNrvc2FnYzBgUkrMqm0gENDt27crAH322WdVRLiingR7oDQ9x+npaZ2entZwOGz1HM3q98q6ldy5s7EYMCkl9hXZ6elp9Xg8WlFRofv27dPHH39cu7u7rQK1tMT+npke4sohtgmg5qwd88HEnmRmMWBSWpk5uPLycj18+LA1x+nxeJb1nAqxR2T27NvzJEtLSz803C4vL+ccZZZiwKS0Mz2ihYUF7e7u1m3btumOHTu0uLhYN2/erMXFxctOr8zXYaW9J2mCYVlZmZaUlGhPT49OTU1pWVmZTk1NLRtuc44yezFg0royvSozXAegW7du1W3btmlZWZl6PB71+XxaWlqq09PTOR84VxtuBwIB3blzp5aUlKjH41G/329VFjJJ6PbhN+cosxcDJm2YxcVF9fv96vP5tKSkRA8dOqQAdPv27bpz506rLJnH49GpqamcGL6v1ou0D7dNQCwtLVWPx2MNt83Z3/b5S/Yksx8DJm04+5Dd7/frtWvX1OPx6OOPP27Ne27ZskV37NihZWVl6vf7lw1b7c+xEcHUviVxYWHhrgHS9CLtw20T9M2Ktz0wrhx+syeZ/RgwKePsyfBTU1O6bds23blzp7XCXlpaqlu2bFERsb5nhvE+n89aTFotqJpgda/h7crvr9ZrLC4utlKlKioqrKIkphc5NTWlxcXFVpZAMsPt1V6bshsDJmXcyjm7lcnZfr9fH3nkEav4h4jotm3bdPv27VpcXKw7d+60KsWLiBWoenp6tLS01OrpmWBn7y2auVUzh2p/bRMUS0pKtLu7W30+n05NTanf79eSkhL1+XxWypRJODeLWRxu56eMBEwA/wFADEAcQM09HrcXwHcAvAOgO9nnZ8DMfSuDqEmx8fl86vP59Pz583rt2jUriPl8Pn3kkUf00KFDVm/PpDMFAgH1+/0qInr48GH1+XxWb7G4uNhaiPL7/VpWVqY7duxQj8ejJSUl6vf71e/3L9ttY3q9pjdpPwqCw+38lqmAuRvAPwPwl3cLmACKAHwXQCWAzQC+BcCZzPMzYOYfewC17y6ylzQzaTtmeO7xePTatWtW6o5ZVDJnF5mFpe7ubutrs1jj8/mWzUGuzIlcuWXR3jYOt/NXRofk9wmYzwC4avvaA8CTzPMyYOa3u9V7DIfDWlZWZvXodu3aZfUazTymSXGybzucnp6+62PC4bDu2rXrQ71E1pwsTNkcMD8FIGD7+jkAX77Hc3UCuA7g+hNPPLFObxdls9XyF1cGNPsWQ3N/ZTC0P4a9RLJba8CUpd+9OxH5cwAVq/yoV1W/mnjMXwJ4WVWvr/L7nwKwV1WPJL5+DkCdqn7uni8MoKamRq9f/9BTEkF1qdakORDM3LfX8bQ/hvU9yU5Evq6qNQ/6e/ctIKyq/3ZtTbJ8H8Av2L7++cT3iNZs5TGzqxXb5VG0lG4bcQTdLICPi8jHRGQzgE8DuLgBr0tElFYpBUwRaRGR97C0sDMpIlcT3/85EbkEAKq6AOBzAK4CuAngj1U1llqziYg2Xkpn+qjqGICxVb7/AwD7bF9fAnApldciIso0ngpPRJQkBkwioiQxYBIRJYkBk4goSQyYRERJYsAkIkoSAyYRUZIYMImIksSASUSUJAZMIqIkMWASESWJAZOIKEkMmERESWLAJCJKEgMmEVGSGDCJiJLEgElElCQGTCKiJDFgEhEliQGTiChJDJhEREliwCQiShIDJhFRkhgwiYiSxIBJRJQkBkwioiQxYBIRJYkBk4goSQyYRERJYsAkIkoSAyYRUZIYMImIksSASUSUJAZMIqIkMWASESWJAZOIKEkpBUwR+Q8iEhORuIjU3ONxfyMiN0TkmyJyPZXXJCLKlIdS/P05AK0Azifx2H+tqn+X4usREWVMSgFTVW8CgIikpzVERFlso+YwFcD/KyJfF5HODXpNIqK0um8PU0T+HEDFKj/qVdWvJvk6n1TV74vI4wD+TET+P1X9q7u8XicAE1Rvichckq+RTR4FkKvTD7na9lxtN5C7bc/VdgPAP1vLL903YKrqv13LE694ju8n/vyRiIwBcAFYNWCqqg+ADwBE5Lqq3nUxKVvlaruB3G17rrYbyN2252q7gaW2r+X31n1ILiLbRWSHuQ/g32NpsYiIKKekmlbUIiLvAXgGwKSIXE18/+dE5FLiYeUAronItwBEAUyq6pVUXpeIKBNSXSUfAzC2yvd/AGBf4v67AP75Gl/Ct/bWZVSuthvI3bbnaruB3G17rrYbWGPbRVXT3RAiorzErZFEREnKmoCZy9ssH6Dte0XkOyLyjoh0b2Qb70ZESkXkz0TkrxN/ltzlcYuJ9/ybInJxo9tpa8c930MReVhE/ijx8xkReXLjW/lhSbS7TUT+p+09PpKJdq4kIkER+dHd0vtkyRcTf69vi8gvbXQb7yaJtv+qiPzE9p6fuu+TqmpW3ADsxlJu1F8CqLnH4/4GwKOZbu+Dth1AEYDvAqgEsBnAtwA4s6Dt5wB0J+53Axi6y+P+IQvaet/3EMDzAN5M3P80gD/KkXa3Afhyptu6Stt/BcAvAZi7y8/3AbgMQADUA5jJdJsfoO2/CmDiQZ4za3qYqnpTVb+T6XasRZJtdwF4R1XfVdXbAP4QwIH1b919HQDwduL+2wCaM9iW+0nmPbT/ff4UwK9J5vfuZuu//X3p0gaTH9/jIQcA/J4uiQB4RER2bUzr7i2Jtj+wrAmYDyBXt1l+FMDf2r5+L/G9TCtX1R8m7v93LKWBrWaLiFwXkYiIZCqoJvMeWo9R1QUAPwFQtiGtu7tk/+0PJoa1fyoiv7AxTUtZtl7XyXpGRL4lIpdFpOp+D061WtED2ehtlumUprZnxL3abv9CVVVE7pY28YuJ970SwNdE5IaqfjfdbS1gIQB/oKq3ROQolnrJ/ybDbcp338DSdf0PIrIPwDiAj9/rFzY0YOoGb7NMpzS0/fsA7L2Gn098b93dq+0i8j9EZJeq/jAxlPrRXZ7DvO/vishfAviXWJqX20jJvIfmMe+JyEMAdgL4+41p3l3dt92qam9jAEtzy7kgY9d1qlT1fdv9SyLyFRF5VO9RhjKnhuQ5vs1yFsDHReRjIrIZSwsSGVtttrkI4HDi/mEAH+oti0iJiDycuP8ogH8FYH7DWvgzybyH9r/PpwB8TRMz/Bl033avmPfbD+DmBrYvFRcBHEqsltcD+IltiieriUiFmd8WEReW4uG9P1wzvZJlW7FqwdL8xy0A/wPA1cT3fw7ApcT9SiytMH4LQAxLw+GcaLv+bEXxv2GpZ5YtbS8D8BcA/hrAnwMoTXy/BkAgcf+XAdxIvO83AHRksL0feg8BnAawP3F/C4A/AfAOlrbiVmb6PU6y3WcS1/S3APwXAJ/IdJsT7foDAD8EcCdxjXcAOAbgWOLnAuD1xN/rBu6R4ZKFbf+c7T2PAPjl+z0nd/oQESUpp4bkRESZxIBJRJQkBkwioiQxYBIRJYkBk4goSQyYRERJYsAkIkoSAyYRUZL+D5myyiACZE1yAAAAAElFTkSuQmCC",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x108299048>"
      ]
     },
     "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": "iVBORw0KGgoAAAANSUhEUgAAA44AAAFACAYAAADziiJ5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXeYFFX297+HGRhyFBEBHRUUUREVAXNGMO8a1vBTTOuadvVddxV1FdecRVxFUVFUjJhAECTnnDMMQxjyMDnH8/7R1d3VPdVd1d1Vfau7z+d55pmqW7dunYp9zz3nnkPMDEEQBEEQBEEQBEEIRSPVAgiCIAiCIAiCIAjuRhRHQRAEQRAEQRAEISyiOAqCIAiCIAiCIAhhEcVREARBEARBEARBCIsojoIgCIIgCIIgCEJYRHEUBEEQBEEQBEEQwiKKoyAIgiAIgiAIghAWURwFQRAEQRAEQRCEsIjiKAiCIAiCIAiCIIQlXbUAKjnssMM4MzNTtRiCIAhCHFi+fPkhZu6oWo5EQX4jBUEQUgOrv48prThmZmZi2bJlqsUQBEEQ4gAR7VQtQyIhv5GCIAipgdXfR3FVFQRBEARBEARBEMIiiqMgCIIgCIIgCIIQFlEcBUEQBEEQBEEQhLCI4igIgiAIgiAIgiCERRRHQRAEQRAEQRAEISyiOAqCIAiCIAiCIAhhEcVREARBEARBEARBCIsojoIgCIIgCIIgCEJYRHEUBEEQBEEQBEEQwiKKoyAIcSHrYCl2F5SrFkMQBEFIcDbsLcaB4krVYghCyiGKoyAIceHSt2fj3NdmqhZDEARBSHCuGDEXlw+fg5mbD6LH05Pwx/r9qKmr923fmVeG4soahRIKgrPsK6oIeObjhSiOgiA4zsuTNqoWQRASHiIaRESbiSiLiIYabD+fiFYQUS0R3RC0bQgRbdX+hsRPakGwj7p6xufztwMACstrcNdnS1FTx7jvy+Xo8fTvyBw6EZ/MzcYFb8zCTR8uVCytIDhDSWUNznplBl74bUPcjy2KY5w4WFyJbbmlqsVQTk1dPerqWbUYKcGktftQVVunWgwAwKg52apFEFzG375chi8X7lAtRsJARGkA3gcwGEAvALcQUa+garsA3Ang66B92wMYBqA/gH4AhhFRO6dlFoRYGb96L4b9ug7nvT4DheXVGDUnG89NCN9ZfnGiZ6By0/6SeIgoCHFnf5HHTXvK+v3YsLc4rscWxTFO9Ht5Oi55a7ZqMZTT4+nfcdNH8R8FnLs1F2e+NA0V1e5QpJxmwbZDeHDsCrz2+2bVoliCmfH9shzXKLqC80xZfwDP/LpetRiJRD8AWcyczczVAL4FcK2+AjPvYOY1AIL9ly4HMJWZ85m5AMBUAIPiIbQgREpOfjkyh07E8p35+Mc3KzFm4U7k5Ffg2V/X47XJmyJq64p35zokpSCoISe/HJe9MwcAcKC4CleMmIsVuwridnxRHONAZY10hvUs3xm/B9zLK5M2IbekKmWsvoXlnrkdewsrFEtijcnr9uPxcWswfNpW1aKkFPX1jD0J8owI6AIgR7e+Wytzel9BiCvzsw4BAL5ZkhNQPn713ojb2rCvGJlDJ+KZX9b5lFHpkwmJxsGSSizfmQ8AuOq9eQ22b9wXP6ujKI5x4NXfIxshE+JDeXUt/vndKuSXVasWJalgZmQOnYg3p1i3dnqDGOSVVjklVtKyKqcQ45bvjmrfD2Zl4ZxXZyA7RQZUBHOI6D4iWkZEy3Jzc1WLI6QguSWe34GiivDBbS44viP+NfB4S21+uWgnAOD6kQvxxI9rsLewAmt2F8YmqCDEge2HytDvpem4fuRC5OSXG74XOw6VxU2e9LgdKYV49td1OK5jSww5OxNA4lh9kpkNBqMx45bvxk8r96BFRjpeuO7kuMqzt7AC9czo2q55XI8bT/43Mwv/uvwE1WIkPde9Px8AcMMZXcPW27ivGGt3F+GmM7v5yhZsywMA7CuqxLEdWzonpGAHewB006131cqs7nth0L6zjCoy8ygAowCgb9++MiFdiCvzsw7hralbAABTNxxosP2yXp1w4hGtcOYx7XFej46oqq3DrM25aNYkDXO3HkKrjHR8c98A7Mwrx0NfrzA8xoxNBzF+9V4wAztevdLR8xGEaPB6A3Vr3xwXvTnLV37e64GR6S86oSOev/ZkdGsfv76kKI4O8MVCz8iWV3E0o76e8fA3KzDkrEz0P7aDg5KlJqHcU70/SgzG4uw8ZDROQ59ubeMi09mvzgAgP1qCOUUVNWjauBEy0tNiamewNtdHrzh6YVEPEoGlAHoQ0THwKII3A7jV4r5TALysC4gzEMCT9osoCLFx2yeLQ257cnBP3DbgaLTM8HddM9LTMO6BswF4vF0AgIhwcpc2mJfVDd0Pb9Ug8mRJZa0DkguCfXwyLxsvT9qE3/5+bsg6px3VFp/d1S+OUnkQV1UXUFJZi0lr9+OvXyxTLUpSEiogztyth3zLfxm1yGe5EWJDlBB7OfW/f+CWUYscadtrcQSA39fuw/Y4ursIkcHMtQAehkcJ3Ajge2ZeT0TPE9E1AEBEZxLRbgA3AviIiNZr++YDeAEe5XMpgOe1MkFICC49sRP+dsFxAUpjMEQEIvKtv/Ln3rjn3GPQrX2zkPvk5JfbKqcg2MFC7bd5VU5Dd+rb+h+Fv1/cHSNuPi3eYgEQi6Nt9Hn+D/z1vGPx0EXdG2wz60ezaY3UpV5L3dGoEZnUFPR4lTeK02U7WFyJfi9Px8t/OgU3G1i0hNhYsSu2uTglBomwvfOIvDwwdgWIgO2viBXcrTDzJACTgsqe1S0vhccN1Wjf0QBGOyqgIETB98tykFdajSFnHx1QPuzqXhhyViaKK2vQtnmTqNuf+/jFmL0lF0NGL2mw7bzXZ2LlM5ehXYvo2xcEu6nT+nArg377M9Ib4fFBPdGmWWMFUnkQi6NNFJbX4I0IgoEI1jjzpWno9/K0kNtr6upjjgopFrLY2ZHnGbX9eWV0QVpS6R58u2RX3COZfmyQR9Mo9Ukq3QdBENRTUV2Hx8etwWuTN+HKEYHRIu865xg0akQxKY1e+h/TPuS2d6dLNG/BPRRX1mDOFk9gsh9XBPapNr84WKnSCIjimHBMWrsvpYLt5JVV41Bp6Kinz/66Hue8OiNs9DXpDMcXudyhKa2qxdCf1uLWj51xPQ1FvdwUQRBcRkFZNU58drJvXe8qP+IWe93wmjZOw5rnBuLNG0/FimcuC9iWK9G8BRfR+7k/DMuvOOWIOEtijCiONlNn0EOzU3F5cOwK/OmD5JmLV1/P2F9UGfX+szYfBACUVclkd7fgtjkjE1bvxcGS6J8xO6nXPgZ5YQZDVCDu8oIgxIuK6jrszCvDPoPf/hZN0vDjA2fhmlOPtP24rZs2xg1ndEX7Fk3w5OCevvLcYo/iOHtLbsiYCIKgilf+fApG3X4G3r/1dNWiAHBQcSSibkQ0k4g2ENF6InrEoA4R0QgiyiKiNUR0um7bECLaqv0N0ZVPJqLVWpsfElGaVt6eiKZq9afqosfFlUe/W+X4MQ4UJ8/o2IdztmHAK9PjFpRDrI/Oc6EudLRqiitr8PdvVuKOTxvObUkljBRDs3dh+6EyvPXHZl+kQkEQBDt4cOxyXPDGLBwKsvTdcEZXrH9+EM44OrRbqV3oo94v2ZGP4dO2YMjoJTLlSFBKYXngoPKJnVvj5jO7YeBJRwQEflKJkxbHWgCPMXMvAAMAPEREvYLqDAbQQ/u7D8BIwKMEAhgGoD+AfgCG6RTBm5j5VAAnA+gIT/Q4ABgKYDoz9wAwXVuPOxNW7416X7c8FNFyqLQKJw+bgrW7iyzvM0+LbOqk+20qW1OceqT+88taDHh5ujONh6GksgaD352LjQZ5OY2o1WaY7y8ObXFclVOIpTskwGQwQ0YvwXszsgytAoIgCNEyc7Nn/taOPP+A8ZFtmuLNG0+NmwxNG6fhvB6H+daHT/PMcxw9fzuyQ6TwEgSnmLnpINbtKcKooHgEk/5xrut0A8cUR2bex8wrtOUSeMKHdwmqdi2AL9jDIgBtiagzgMsBTGXmfGYuADAVwCCtLW+PMR1AE/inVF0LYIy2PAbAdc6cmXPU17NPkUpE5m09hNKqWnwyr2EgDqeJRTVMXbUyer5atMtQGXPaOLVwWx427ivGW39ssa3N696fjxs/XGhbe+Fw1+ffj9F9q66tj78ggiAkNZU1flfQpTsKfMs9O7eOuyyfDjkTC5+8uIFb7MVvzY67LEJqc9fnS3HVe/PwwaxtvrKv7unvOqURiNMcRyLKBHAagODMrl0A5OjWd2tlocq97U0BcBBACYBxWnEnZt6nLe8H0Mke6eNHSVUt/u/TxZix6YDhdnEZC82B4krsLjCeW5eKly3eVlYXftsEHan4DgiC4D7G67yyJqzeiwuO74hv7xuA4Tf3ibssTdIboXObZnj6yhPR/fCWAdvq6xmLsvNC7CkI9qEfTPFyXZ8jcdZxHRRIY47jiiMRtQTwI4BHddbCmGDmywF0BpAB4GKD7YwQhiQiuo+IlhHRstzcXDvEsZ1kmsOoxyiAjV0d2j9/sADnvjbTnsaSCHKtjcsYp+R121VIlkGgksqapDkXQRCc5/FxawLW/335CRhwbAe0bqouxUCn1k0x7Z8XBJR9PDcbN49ahLlb3dlPFJKD1TmF6PnM5ICyv/TthuE3n4Y0l+Yvd1RxJKLG8CiNY5n5J4MqewDos4V31cpClftg5koAv8LjogoABzQ3V2j/DxrJxMyjmLkvM/ft2LFj5CcVAXZ/cBKlfxZKzps+Cu8OuP1QGbo/NSlugXK8JMp1FfTITYsXoazJOfnlOOW5P/DFwp3xFUgQhISk1GDw+OQubRRIYs7ynR432oNJOpAvuANvvkYvjdMIr93QW5E01nAyqioB+BTARmZ+O0S18QDu0KKrDgBQpLmbTgEwkIjaaUFxBgKYQkQtdcphOoArAWzSteWNvjoEHqVSKZ/N32GpXqooLuv3hjc4/7xyD2rrGb+u2hO2XqRwwHKKXOwkxczfn5kDUuIkgjXs7albkDl0oqPHcOIqeANbTN1g7FovCIKgZ/P+wD7A5EfPUySJOX9o37VN+21xlBOEAN6eugUPfLU8YGD2sl6dMOfxi9QJZREnLY7nALgdwMVEtEr7u4KI7iei+7U6kwBkA8gC8DGABwGAmfMBvABgqfb3vFbWAsB4IloDYBU8VsUPtbZeBXAZEW0FcKm27hIi67aF6hpH0/lblVPourx6iUB9PWPm5oMJoXhEijdi3OR1+7ByV4FJ7dgIdf2qausM82UxGDV10Qdl+fe4NTjuqUlR7x8rB4sr8eRPa0MGljFSfEdM32rLsQcNn4PnJ2yIat/ke8oFQXALa3cX4ZTnpuD6kYFeRz2PiH9AnHCs++/lWPffy9GlbTNf2cdztyuUSEhWRkzfit/X7UdRRY2v7Mg2TdG5TbMwe7kDJ6OqzmNmYubezNxH+5vEzB8y84daHWbmh5j5OGY+hZmX6fYfzczdtb/PtLIDzHym1ubJzPx3Zq7VtuUx8yXM3IOZL9UUzZTnuvfn47zXQ8/9m7R2H/71w+o4ShRIQD/aRUra6PnbcddnSzF53X7VokRFqEv566o9uPit2Zi1+SDu/2oF/vTBAkflCHX9zn99Jk58dnKD8u+X7UaPp383nCxuhXHLd0e1n10MG78e3yzZhekb42+F27S/BKPnG3dyonm1xDovCIIdjJqbjZLKQDfV9i2aKJImNC0z0tEyIx3X9gmMsroqp1CRREKys1qXvm7wKZ0VSmKduERVTVVCWVtKq2qjyhPkhPXrwbErlHe2nSaa67ZLs9IeCJP/LxHx5tjccqAkLscL7ix4MQsAZWSNTAS8j5rZEycqmSAIqQAzN8hvnZHeCDP/daEagSzwyKU9Atave3++IkmEZCTroL//tWR7Pi49sRN2vHolBhzrziiqwYjiqIDbPl6UlHmCEjclQxJ344PuSeLeo9SAmaO2tgqCILiNcoNBwM0vDkabZuqiqJqRkZ6GC453NniikLq88NvGgPUHLzpOkSTRIYqjg4RSR/SmaTvaSzYiSclgpWYs162OPfMdo+Xv36xELwOXzFTBqrtjpFfYqhHZW62g3B1pIxab5AUbOXsbej4zGfll1bYeNxq3U+/lUpXS5V8/rHY8aJAgCM4SHCX9/wYcpUiSyPj4jr4B6zKgJ9hBfT1jdlAk1dOPaqdImugQxTGJYGbklkjoaDt54bcNuHvM0qj3n7B6r+GIqxnB0UEj3j/qPWNDvWoWmoXbnE/mbGbRvWfMsrDbf1npiSh8sCR6F+mi8hpkDp2In1dad0F3g1IdTLK70AtCslNSWYOr3psXUPb/Lj1ekTSR0Tgt8GPe85nJ+Od3q1z5rRQSh1W7A+fLPntVL0WSRI8ojgmE2fdq5OxtOPOlaQkVRVV/TsqUHYMD68tmbY5/AuAPZm3DcU9NQklljXnlKEjF375ih66l2/CmyQhIBxTD/V6zO3xgCAmiIwiCEb+sCpzb+MYNvdGhZYYiaSLDKAL2Tyv32O4NIqQWwQPYt/ZPDAu8HlEcXUi089C8Cs6ewgobpYk/ds/DC1BOE6SP+82SXQCAwnJ7lR2zPIheDpVWKXXNSfa5mHY9h7UxpC6xyn1fLndsAEMQhORl4bZDvuWv7+2PG/t2UyhN5Gx6YRCOaN00oOysV2YokkZIdPJKq/DGlM2+9aaNG6Fp4zSFEkWHKI4O4ouwaFMnMRlH9heazPlyM5v2FztuxYr12Yk2J2jfF6fhL6MWxXZwGMsfiwtuLCTKoAFgfV7hR3Oyoz5GlS7XpNmlefjrlVEfp7C8OmWsvYIgAD8u342Zmw5iZ145mjdJQ7PGaTixs7tyNlqhaeM0zH78woCy6jgM1gnJxwezsnDGi9MCyp6+MvHcVAEgXbUAgrvnhTlJos/HHDR8Lk7p0gYT/n6u7W07ZXGLpNnVDuWuevOPzeaVQhDpZUkkZTEarHoXGF2GS9+2Htk562Dk6YO89Hl+KgBgx6tXRt2GIAiJw2O63NDX9TkSw28+TaE0sZGRnoY7z87E5wt2+Moqa+oS0lIkqGP4tK0B689fexJuH3C0ImliQyyODmLUWdNPrP511R7M3HTQentJ1gkOdoWM5/k9MW5N2O1Wrbtr90QXIddpzCbw232tI7FsmkUWtYON+4odP0Y43PCuxiqDfndmxsJtefhu6S5fWY2MvAuCYIJZzt5EYNjVgZah60cuUCSJkIjklVahujbw9/KOszLVCGMDojjGmd/X7fctP/LtKtz1+dIGnXy3KiOJxpYDJfh8/nYYqS/fLcvxLbuhky9Ez4Jth/BjUATOwe/ObVDPjbe5otr5wZO5W3Ntse7f8vEiPPHjWt/63Z+HjxArCELqMXre9oD13t3aKJLEPogIjw86wbe+fm8xDpUmvkIsxIf1ewMHss/rcZgiSexBFMc4s8tCxNOvFu0yrROOmZsbWjHXpZAy+tLEDZizJRdXvDsXz03YoFqcsFTV1uHH5bslxHcM3Prx4gDXqETCGwQpGDtclb1t3P7pEvyspfkIiQOP34a9xfJcC0KK8fxvgb+5j112QoiaicWDF3ZH13bNfOtb9pcolEZIJOqDfgffuvFURZLYgyiONmDqFmi2v22CeP59NLthwIzgUUAncaKrGMpiYhQl9OO523HH6CWoVRSEpa6e8cBXy7F8Z4Fp3bf/2ILHfliNaRutuyzHRJJHKw1HPHSYSBW+4B8UL1ZlNatXVlVrrR2TtzbSSzd7Sy6uGDEX3y7NMa8sCEJSEDz9ZN4TF6FJevJ0M4/t2NK3vC23FL+t2ZvwUewFZ9lfVIk7PwvMBX54UKTeRCN53mgXYnW0PVUH5SM573u/MHaLs/LRNjuO3dFqc0uq8Pu6/Xhw7HLTugc1hThUugOnIulabXev7vp+vXgXxugCBBi2G3SxDXNkWjiuWVTRRH9lvNff6Xf/7FfVhI7fnusJpvPkT2tNanoC73wwK8tpkZICIhpERJuJKIuIhhpszyCi77Tti4koUytvQkSfEdFaIlpNRBfGWXQhBXh9cmDgs67tmiuSxBneu/k0jLm7H5qkN8L2Q+V4+OuVuHnUQtViCS7me920qGtOPRILn7xYoTT2IIqjDcRL8YvlOG7vaFOIZS+5xZUNylZZjPrJIZYN64a4yK9P3oTMoRNtDwgS6nhmilPWwVK8N31r2DpAQ4us1TQPXr5ctNO3/NTPazFs/HrsL2p4L0JhpqAeMLivoVidU4i9RdGP7iZSOhvrlsvQ58QMFFVYS4NhqODH6Xtz00cL8frkzSivtmYdTVWIKA3A+wAGA+gF4BYiCo7nfg+AAmbuDuAdAK9p5X8FAGY+BcBlAN4iIvn9F2wjO7cUo+f7PZuuPvVIhdI4Q5vmjXHB8R1RXVvvO9d9hdZ/w4TU47c1e33Lr9/QG53bNAtTOzGQHw4XEEuHtqK6Dttyow+Vr5JozvvNKZuxIMuTVLjYYqc4Ej6YtQ1Aw07zZ/N3AHAukmSkLo43j1qEt6ZuiTo/XixKQSzpNIDAgYH+L0/HtA0HLO137fvz8eyv6yM6VvAzVllTh5cnbbTswuk0evmYrb0R+4sqfbkw7Rq0UuTVDaBhgCAhJP0AZDFzNjNXA/gWwLVBda4FMEZbHgfgEvKMHvUCMAMAmPkggEIAfeMitZD0vPXHZlz8VmB6nz7d2iqSJr7U1rN8w4SQ5ORXoHObpjivx2FJk8JFFEcbiFefy6hb+bevluOSt2aj3qTnl+hBKrzS/29mFm79ZHFE+74zdYvlur+t2RdR20ZEoxBHenuqaox/qMqqasNa8SJVUGNO6WDBVdXIDfmFic4ENRq7eBdGzcnG+zPtdY2srKlD5tCJmLTWEzXZ6jMQyfWdtuEA9hVVYMAr0yN6phse07nn84XfNuDKEXNTeSqtk3QBoJ80ulsrM6zDzLUAigB0ALAawDVElE5ExwA4A0A3o4MQ0X1EtIyIluXm5tp8CkIy8t6MwO9pk/RGuPucTDXCxIEpj54fsD5BZ1USBAB49td1eH7CBlTU1OH2s47Gl/f0Vy2SbYji6AYsdMp25pUhO7esQfm8rZ4f9lBBNiI4RNzRi8yIXOGyWnvBttjzBlaEUNTCYcUt1Ci4TyxcP3IB+r88XSeDPRzUKaPxGIMYF5ReI2qCZK3VLMZ2B06KxOVWT2CuxPB152875MuJNmervR16vTI5Y1NDC7DVd/PTeduxfm9xRG9yIrkQJzCj4VE0lwEYDmABAMOPGjOPYua+zNy3Y8eOcRRRSBZWPzvQ9t82N3F8p5bmlYSUpb6e8cXCnT535lZNGyuWyF5EcbSBUCP4ds4buuCNWbjqvXkNygM+zknync4rq47r8ZxShKx0iO22BG9yKET4fV+aB/pxI8GKmVtUFLPb7uSrbGYFXp1TZFDqPJHOv01B9iDQSthVKzOsQ0TpANoAyGPmWmb+f8zch5mvBdAWQPRma0EIQ7MmyeGSFwoiwk19u/rWHx+3RqE0gtsINjSUVrpjaoxdiOLoAuzonpkHfdEvM/KDlLPvlu7CLIP8j04SSubPTSJ3xoNog+hES6jB2VCHsXr0gvLA+xyt9bVQ106iWojcKHWw1d2sbvBjYvR8RDPQbxp52I0XL/VYCqAHER1DRE0A3AxgfFCd8QCGaMs3AJjBzExEzYmoBQAQ0WUAapnZ3UluhYSgNmje/1NX9FQkSXwZdvVJAetF5fbHXBASk+A0VBf1TC7PDVEcbSBUn8ofct9dva5vluTg9BemYssBv3XqiR/XNsg1EwlTNxzAtoPOBelhhi8gSKIQiwXFauffrNrcrYcCcmutNohEW1/PISNaxqokGu0dz9ch+Fjxsmk5Z8W2fuyQOSKjbtlZEnVAIl5ocxYfBjAFwEYA3zPzeiJ6noiu0ap9CqADEWUB+CcAb8qOwwGsIKKNAJ4AcHt8pReSkX1FFej+9O8BZfedf5wiaeJLi4z0gPXXp2xSJIngNl74zT8m9+/LT0DPI1orlMZ+0s2rCE5jRyeTPZMEQ2+HR0H474T1mK9ZnbJsVPT+GiLPYjgiVahv+HBBTPuHlCOafWzq48baTMhBC92Gqtr6BtG89Nfu3elb8e70rVj97EC0aW7BFz9B+/dWo5bGkwZRVS0+WFYU4PV7i42PyWyxhfjgsnE1V8PMkwBMCip7VrdcCeBGg/12ADjBafmE1OKNybFF2E4mxi7ehbGLd+GzO8/ERT0PVy2OoIjKIDfVhy7qrkgS5xCLow1E0vFRmdA9p6AcYxbutFVhNJTFoZ7gyl3W8jaGw0w27/ZE6cyyDdlBfl3lmSaVX24wt9TidfBaSBtUt3ghI73eVp+xUO9FPNSm2rp609yEoU4j3nElTF1VY9jXLhkEQXAXTYPmMj5wYWpYG71M/X/nNyiTCKupzRM/+ue7Xtm7s0JJnEMURxdghzJpPkfJfb0yvUSEKJQHO4WJkM8szMOsrfdrdSUh8i3Gqh+4z4bmLmKdI2oVI7fke79Yhl7PTrF8/Gjm1UaV+iWKY0d8DJMG1+8talAmT7IgJA719YyvF+/yrU9/7AI8MSg15jd66dGpFV6/vndgoXzIUpZF2Xn4dZV/4OCKk0VxFEIQqvNmGLjCZluHexzOjBm/ei8+mZutWgzb8N5TvQ97KA6V+i14f/7A72abV1rlb88meewgXoMLbv1dLa6saRDoAQByS6rwxcIdEbc3a7OaHHhWrJU3jAx2+25Yx8nH4coR/gjRPmu1Cwe3BEEwJthD5biOqZmi4qYzA9OhylcsdXnl98B5rt3aN1MkibOI4hhnjJRM2+Y4hj2uGlewf3yzEi9O3Gi4LWYX3yjOx6k5hWZs1bkH3/9Vw9QWoQYUTK1QVg5uohREmm/L6JifzN0eURtWeOrntbhj9JKY2gh1fYzOuPdzf+DfBmHVH/p6BZ79dT225drv4q1XlszShTAaKlfRvNPMwLKdBUFt6+SIvMkGRPJIib4oCInHmt3+qSPL/nOpQknUc+mJnXzLP68MzpAjpAIPf70iIPjgT8wsAAAgAElEQVTgGUe3Q++ubRVK5ByiONrALyE+FEYdohKH8rlE5bIW4S65JVWGycETHePrEH1v1uy67i00Txbv7XdbnZMZC+Ha8G4xUy5/X7c/7P7R8PXiXZizxRmrnTdYVDBGP/oFWuqaeEb1NZ2HHOdJkJHPQW1YNtMk3Y/oj4KQONz9uT8g3mEtMxRKop5zuncIWA8OkCIkP7+t2edbJgLG3ttfoTTOIoqjDTzx49qw2/UdopGztoXdHil+Ny+Tihx7X/PWjxfh7s+XocbAnS86rJ+5nVaJ32KYvB5vdzqzo4XSZcwGEjbtLwm7vUF7uvPOdsDyFg0HiivxhkkI9GCLnpdRc7Jx7FOTDJXHeBKYx9G/cv9XK7TtxvLHdEyLXg92X5nXTSIwiuVREBKP4X/po1oE5bRuGhiNPCe/XJEkghvo3Lppg0j2yYQoji5n2oYDWJQdXdL2YGLtmGUfKrNFDtU8/PXKqPcNdwkra+psj1hr7oIc3U2duFY/OuYfUSiurEFuSZXRLj5W724Y2CQUTioD//x+Fd6f2XAgRs9PK/wWRKNrFSrXYayYtsoW60VxnAPF4e+f1XbijiuEEAQhHDV19Rg0fI5v/ZjDWiiUxh386bQuAeuXvTNHrI4pTMumyZ3p0DHFkYi6EdFMItpAROuJ6BGDOkREI4goi4jWENHpum1DiGir9jdEK2tORBOJaJPW5qu6+ncSUS4RrdL+7nXq3KxitVMfzop17xfLcPOoRRaOFbksVuSrqavHV4t2Yvuhsohc9ZjZ9MMZfNqR9hvtiihqVzt//2YlLn17Nqpq7bLIAqZXxYLo3vOzYi0999UZOPOlaQ3KI50HaSeZQydi6I8N5x5W1YS/znsKK/D21C0BZdHoiZat+iFYkHUIG/cZ51Q0I5RFMhxmOVWNA+GEb/tQaeTKqNVjiL4oCInDuj1FAR4raY3cHqLPeRo1ogbRVe34ZgqJQXBQvU/uOFORJPHBSbW4FsBjzLyCiFoBWE5EU5lZH45yMIAe2l9/ACMB9Cei9gCGAegLT79iORGNB1AF4E1mnklETQBMJ6LBzPy71t53zPywg+fkWow6ZRe+MTPmdj+fvwMvTTIObhOOD2dn47XJ4d0I3UzINA5herkLsg4BaPgRcUYOTRmMYD8rik9x0Bxcf15Ltd37b5fmBKxbkabGggJvpR1v4CKGZzCktp7RMsP6p/PWTxYDAHa8emXD43PQsslgSoP1JNK6QkenZqUDF4Ig+AkeP+7SNjkjR0bKTWd2Q+ZhLXDTRwsBAD+v2IO/X9JDsVSC06zcVYA/fRAYpfyoDs0VSRMfHLM4MvM+Zl6hLZcA2AigS1C1awF8wR4WAWhLRJ0BXA5gKjPnM3MBgKkABjFzOTPP1NqsBrACQFenziFWrHbqbImqalC2I8/vZx/tMQorGiaFt9KWN6m8XQQfsqCsOmE7zVbmrXk7yqaWZCsWR4tukXZczpW7Ck3rrNtj3dU1HGa6hNH24LJIn6ErRszFycOMczPGE0fVKEWDBdEM1giCEF/0OYnH3N0P7Vo0USiNu2jfwj/X8a0gbxchObnlY3OPwGQjLnMciSgTwGkAFgdt6gJAb0rYrZWFKte32RbA1QCm64qv11xexxFRYHKdBCMnvxxT1htHqtTjs4ZE0bnSJyqNBPtcRK0TPO/utBemYtJa8+tjO5ZcQyMnlBJk5xzHWBSBaC0+Rse0KzqpWU7U4O3MRu7R1mUZNHwusnPtnedrdnwnguMYH8fm9mysLHqjILiDGZsO4M7PlgIARtxyGi44vqNiidxF98NbYczd/XzrqoOvCc6TFtQ3mvzoeYokiR+OK45E1BLAjwAeZeboJvo0bDMdwDcARjCzN7v8BACZzNwbHgvlmBD73kdEy4hoWW6uswm6vZ+MaDplV46Yi7992TDfn+nBQm022D51Q3SpNazkjIwVs3laALB4e+RBg0xlt0FRi+b6hLS2hDieV5GL5HdJxU+Yo8c00WUnWIiea5fC1NCSaaIQ+uad2nP8SDCOoBofQeZnHULm0ImG20K7XUvnSxDcgD4Fx4Bj2yuUxL1ccHxHtNessCNnhw/eJiQ2lTV1KKv2x/JY8cxl6HlEa4USxQdHFUciagyP0jiWmX8yqLIHgN4y2FUrC1XuZRSArcw83FvAzHnM7DVLfQLgDCOZmHkUM/dl5r4dO7p3tCx4rlkw2w+VIXPoRFRr8+ni1fGzk3D9QStKrZERbOWugoaFpoIYFEVxOZ24A6ZzHG3oVIfTv+LVZ3/q5/ApbaLhjSmB6R+Yo0tJ4+T0usJyvyu4WUqM/UWVMVlMwxGP+8wMfLNkV8MNvuBDEjxHEBKF4BQUgp/uh7cE4PkNWr4zX7E0glNM1uWvPq/HYb4Bg2THyaiqBOBTABuZ+e0Q1cYDuEOLrjoAQBEz7wMwBcBAImpHRO0ADNTKQEQvAmgD4NGg43XWrV4Dz5zKhCCaTtu0IMXKqbQNTmGUYsSOzmvwJGW7caqDHbWrqoV5YWaWbyefDCcVkmj0uQaKl6LXoqaOccuoRRizcGfYenr5pm86aKx4RYjbvgVe3CmVIAgA8Nz49b7l/1x5YlLnqYsVfWT160cuVCiJ4CRNG/tVqPdvOz1MzeTCyaiq5wC4HcBaIlqllT0F4CgAYOYPAUwCcAWALADlAO7StuUT0QsAlmr7Pa+VdQXwNIBNAFZo7nr/Y+ZPAPyDiK6BJ5prPoA7HTw35agMMmhHZ/vmUYsS0hfcyqmbz1uL5HixX2x/BNbI23JzZz7Sd8DoXJYpHA1eGDR4whbu0LjluwMLorGMG+xj5PIcadNm+T9DEqP7uCAIzlJVW4fPF+wAANx/wXG497xj1QrkcpqmS4r0VOD+r1b4llPJAu+Y4sjM82BiFGBPj/ahENtGAxgdVLY7VJvM/CSAJ6MS1imsRlW1Qzkw2x7lIYwCkLjFYmEWHMUqdp+NHR1d75nFOifTSt1ETXQQ6f03coW8/dMldonjCKHnuNp8HBse2hd+82dayjPIYWamGOtF0OeAdcv3RhBSlTtHL/UtS9pGc975Sx+c/eoM1WIIceKVP5+iWoS4IsMiDlLHHLfADvEMIGHXodxsSQiXU8503zgGrKm3Ik+MxxDiQ/Ct1CtPth4nxu2h0D+LH8yKPCiE/p3734wsf7k8wIKgFL13RI2NeYqTlSPbNgsY3NtbWKFOGMER9hX57+nlJx2hUJL4I4qjgyzfWYCvFoWfwwREF900OD1CcN/qUNCIPxvUsXachmV29ePGr44uHUgsXD8y/BzIWKwb/nyJ9vV07Yhg65NLxRxHB9tO9pzwPZ+ZHPKeFVXU4PP526O6vkYupfHQzdjsI6TbVloVPjiYIAjxITh9UnWtKI5WWDNsoG/5po9knmMyUVJZg7Ne8VuUW2Sk1nxfURwd5scVe0zrvD8zy7ROMMF95uAOZt8Xp0XcZjAvT9qIFdFEKbXIyCisEnpUKA6RKGrB1NbV4/XJm1Cgi6RpdgpfL9mJD2aFfj4sWUC9qR9CSL/9UJnWllH7ps0rwa1yeVm6I/K5kwzrngPZuWV4bsKGqN7PW0ZZS1gc6TWO9pbYOSdYEAR7KasOHMQ5/eh2iiRJLFo1bYw7zjoaALC7oALr9hQplkiwi0/nbQ9Yz0gXxVGwGbN+TzQ5YhvkjbPQ+YpUzxo1Jxvzs4yin9rrrunbJ+gczI6jwuAUTiSi8HWmbTyID2ZtC4i4Fqq5bE2Z+2ZJDl6fvLnBdg76HxYTi6OTvPDbBmzLLXWkbTuC4zjFV4tij35qhWhG//cXVzYoc0OuxFASyBxHQVDDmt2FmLhmX0DZtX26KJIm8SiuqPEtD5+2RaEkgp2kJbu7kwmiOLqAIt3HxSoNHlsLfatIul/hOpLeLTM3HcQPy3JiO1BYGexpx/w4DQ8UleKr7RNq3mFtvbVOftbB0IpWTV093pyyGeVa0lkr7qfvz8yK+kdryY48MLPPKhkNj49bE/W+4QgOjvOfX+zPBRlPonkO7fr5CkjfYnLM+79cHtU8p0gCeOmfNxfotIKQklzzv/l48qfE/q6q5IpT/Fnipm08KPNDk4RmTfwWxrXPDQxTMzkRxTFJsDtESrjOmnfbXZ8vxb9tUgpqDMyuZhIHz/O0AzvyYcbazz1Y0tAi5OXnlXvwvwDXZvMgPmMW7sTwaVujkmvdnmKszCmMYk8/4c4nFoJvv6mVj52zXllRbga/O9f24zrxDkzbGH7O9eT1+zFGC80P2JdeRr9t9pZcXbkgCG7gpwfPVi1CQjHwpCNw4Qkdfes/LNsdpraQCIyasw0vTvSkiZ/z74vQKoXScHgRxdFh6uoZc3SdILtwosOoJ2y0Tgd6ckNGuzstgg+Dc18elAsw1LWzw3ISPGIZiZtztO6ItXWxCW7R0OparLxrVhTSjfuKTdqIrl070B9l0/4S0/p6F6yF2XkoKq8x/S5Yff7c4DYrCKlOfdCPy0MXHYfTj5L5jZEyesiZvuWnfl6LLxfuUCaLEBuVNXV4edIm3/pRHZorlEYdojg6TDRuqFZoMMfRtNMWWfvRzLu0m7ilMolx/+tHLsTMzQf97cXYYCT5CSOJlOrE1bRyj6y66MYDu3J/OkWk98iu8aNYA+G8PGmjbTLUBn18RJEUhPgTHJTtXwNPUCRJYtMoKPHlM7+uVySJECveKUKpjiiODhMcytop7LZMhLM4qrCCxAuzvIihtq7YWYAKX969EBbH6MUKyfUjF2CZUQRPGyOkxtpxd+odiNTqzqYp6KPHDt0mmqi2timOBtcl3H0Pfk+CO0fRyeAh+HkRtdEPEQ0ios1ElEVEQw22ZxDRd9r2xUSUqZU3JqIxRLSWiDYS0ZPxll1ILJbv9EdsfuSSHo57OQmC27n786W+5S/u7qdQErWI4hgjOfnlYbc71mkOWreS7y8SBcDKHMdIMHPTA4CPZmfHfJxY2ZZbhrzShnnuzND/yIa65VsPmLsARkppVS2e+rlh8AJDBcnkes7cdDB8BQOs3KNgC5IdxOPRYGZ8vywHFdXxySm4ycI7EkwjmzpzEVscg+p3advUfB+T43i/T8Hu2GJw9EBEaQDeBzAYQC8AtxBRr6Bq9wAoYObuAN4B8JpWfiOADGY+BcAZAP7mVSoFwYiZm/1TbB66qLtCSRKfFk1SK11DMlJfz1ilxXt45JIeOP/4jiZ7JC+iOMbIfyeEdzswCn1vC0EdRvOIhZHZWsJbHMNTUlWLzUFKUjSBQVSF4c8pqAidGiDEBv19HjF9q2Gd92Y0zMdopMxHqgsY6WWG1iuT6/lSFO6GVu5QvAZPzIhUAVm8PR+Pj1uDHXnhB4cAexTZvxjkVzQNEGXDcQFgS4SDGsG3tKYuemtu8H1pMKdWFEcv/QBkMXM2M1cD+BbAtUF1rgUwRlseB+AS8piKGEALIkoH0AxANYDIRyqElKCgzJ9ruFVGOpqkS1cxFmb9+yLVIggxMnK2P+94765tFEqiHvkaxEi8XFHNsHseULjWzI7lRDAgI5xynKkLMycvVOdYL8vcrYciPqaZG1C4dBjBQQwA68qkFWJ9smINrqOKsirrlkan5uHFa37fFwt3BqzX14dXA4O3vjt9K/JKq0PU9u1kyZOhwRxH0Ry9dAGgz3+0WyszrMPMtQCKAHSAR4ksA7APwC4AbzKzgY87QET3EdEyIlqWmxufb7ngLsr0XhbioRozHVtlBKxLWo7EgpnxxhR/Pu1mKW5BFsUxRpxww7NCdK6q1ts3m+sXD0xFsOkHbUaQi6aKb7re5XlRdl6D7X/+YH7IfZ2+V0aKqRcrio1zcxwjq29Vipq6erw/M0s3Z9W9ODXvqM5srq/B5li9Kw5o+wcHU3LBpygZ6AegDsCRAI4B8BgRHWtUkZlHMXNfZu7bsWNs7ljP/roOH+lG6oXEoFQ3aBbu+y9ER4+nf8fKXQXmFQVXsC03MLe2XVNEEhVRHGPEDQqWdSKY4xhGeUqkM7ZCSWWgZSlcFNBP5243zEsYawdeP5o1fFpDV9diTUajqKBGv+tGz2W0982s37C3sCLsdjMlJFqc+nR/tzQHb0zZjPdnhu7wPjR2Bb5fmoPF2XnYvL/Etnci+FKpetdenrSxwXuhx05LqLelGz5cCKChhTrZvjcxsAdAN916V63MsI7mltoGQB6AWwFMZuYaZj4IYD6Avk4LvHxnAeZuPeQazxzBGqW6d98ox7IQOVf27hyw/sm87YokESKlujbwHUh1i7EojjGiKtNAcJJu02igJm5iwYSd4xin35GXJoafc5edG9p9Mxbq60Of4yfztuOhsSsalKscfzLqlP3nl3UNyk5/YWpU7ZvNd/1TGGso4J6oqlap1CyNpVWhU+lMXLsPj/+4Bn8ZtQiXD5/jiBxWcGrg87P5O8Juj+aWMqy5nUo6jpAsBdCDiI4hoiYAbgYwPqjOeABDtOUbAMxgzwXcBeBiACCiFgAGANgEh+nQMgPzsg7hFoP5u4J7eUD3G/fc1ScplCR5ePOGUwPW2zZLvcTxiUpFTeAgau+ubRVJ4g5EcYwRp6wpZuSWBEb+vOCNWVhilJYhSj6Zlx1yW7zmHH25aKd5JQcwyzuYX2Yyl8sBwnWene5Ymw0iHCiOPAqtHczaHFkUWKuXyeuGUhfB3EynboFpOg5nDmuKk49cbXBUVecOlVBocxYfBjAFwEYA3zPzeiJ6noiu0ap9CqADEWUB+CcAb8qO9wG0JKL18CignzHzGqdl9j6fdv42Cc5SWVMX0L+4tf9RCqVJHpo1ScPsf1/oWx+7eBeqat0/HSLV+WRuNsYu2uVb79K2GdqkuNKfrlqAREeV/380Cbut7lJSWRPWTS/Ze3Jm1lunLGhLd+TjzMz2htvCHdHpR9CNBp99hRWOnDcR4E1JGNmgkFOao5KjmmL0juw0iT67/VAZ9hU1dGsOVn6DLY5/rD+Ak7u0TvlRXgBg5kkAJgWVPatbroQn9UbwfqVG5U6TddA/N2j6xgO45MRO8RZBiJBvl/g7yRedkLopB5zg6A4tcOxhLZCtBbsrraxFRsvUDrTidl4M8nw7oo156qlkRyyOMaLK4hjN3Eqru5hViyZqaCJhFgXUKCBStC6D+pZu/HBhVEqp0/NsS8JEF1UV8bK8OvKRWiuyMvuT2asYE4r0erpRqQ9HZU1Da37w+xT8/j3181pc87/w7tCCO3npTyf7lt+f2TAdkeA+dmqB2n568Gx8dlfqJjl3ip8ePNu3HG4eueBOPvy/M1SLoBxRHGMkUeaNe/I42iPsYz+sbtB2MlFXH/5aORnooTyKZPNOK47/+Galo+1Hg5PPnHfuZCTeBHaJk5MfaJEze2dVvXtOPXOVNXWoUTVxXLCdC084HF3aNgMArNhViHV7ihRLJIQiv6wamUMn4rP5O3BsxxY4/ah2qkVKSto2b+JbvvDNWeoEEUypCBqgXvzUJQ1Sq6QiojjGSKK4qka7j8p2VWGWYsVo+5YDpQY1zQk2VAZ/qIJZajBXSOXghap7H81hrQSIIvLfE1XeBJGg6t47dWme+nmtROBMMtIa+b9yV703T6EkQjj06SES4NMnCI7zt6+W+5afHNwTnVqLmyogimPMqOrkSFJs54jnHMfglkK5YHpF+nllcPT9REsJYw9OnrMvOE4E97moInQE1lgwO81YcydGi1PXf3VOYcqHOk82goMdCe6kRYY/5EXXds0USiII7mDOllzfcp9uMsfeiyiOMaKq0x6p7lJbx5ZGEa8fuQB7CsLn5QMQkMsw2dSW2joOe1JORlWtjqKTVVfPyC+rxrBfG6bgSFacfO3SGkV+DG8OwnhTXaumU+7U5a+uqzedYywkFu1bNjGvJChnf5H/N/31G3orlCS18KZ/EtyNBGfzI4pjjKhSHCOd27QwOw/jlu82rbd8ZwEGvzvXtF6/l6ZHLYvbUemiGG3Qow9nb8OYhfFPX6LqUh0scSYFCDNA8AbHUf9cq5fAGKfe+eraenFVTTJG3d5XtQiCBR79bhUAoFv7ZujcRiyOTjLs6l6+5ad+XqtQEiEUetftG87oimZNJPqtF1EcY0SZq2oUh524dq/9ggDYYRKGP9F4fNwazIlX5Nig+xjNfa1nRrvmMqpvhtVL+/Qvnh9yNygwbh2UcSp+TXVtvbiqJhlHthUlxO3szCvzLX8ukVQdR+/2+NOKPShQkBtaCM9aXSCvRqoSJrsUURxjRFXf0g3WEC+Xvj1btQi2c6hUTVL7cLf1vxPWG5bXM6OzotxCCTXX1uI7U6O5SlYpcgNNBJy67zV17AqFXXCOYb+uw45DZeYVhbhx1Qh/0KJmjcWy4jTHd2oVsL6n0Hx6kBBf9LqiTJ8IRBTHGFEXHCeKfeTZdx0T1+4LWA/XIf9s/g7D8vr6wMiFQvLg1lfWqc9edW19wqQ4EqxzXo/DfMtjFu7EX0apmRMsGKPP1SuKo/O0yEjH01ec6FsvrnQmuJpgD2aR9lMNURxjRN0cRyWHFRwmWldVVY9DIj2HL03amHBukG4d6bQyXzoaquvqXeVNIdjDB7edHrBeXCGJz91CcEoxmcsVH+497xh8e98AAEBhuSiObmNbrt8rQiKqBmKqOBLR60TUmogaE9F0Isolov+Lh3CJgKo8jlEFUXFADsFeolYcVQ1gKDlqdFTW1GP4tK2qxYiIf3y7UrUIcUcUx+SjVdPGmPSP83zrFTV1yn47hUB+XR2Y4ikjXewJ8YCIcHSH5gCAB8euwCnDpiiWSPDCzPh8wQ7f+l3nZCqTxY1Y+UIMZOZiAFcB2AGgO4B/OylUIqEqAmc0h5X+mPuJxnZYbyGxvZCYLNmer1qEuCPPcnLS68jWAevBbvqCGvTWrktP7AQimfYQL/RB7fTuwoJa1u8tDliXdyIQK4qjNyvslQB+YOaicJW9EFE3IppJRBuIaD0RPWJQh4hoBBFlEdEaIjpdt20IEW3V/oZoZc2JaCIRbdLafFVXP4OIvtPaWkxEmVbkjJVE8nyTZ9/9RDsIr8pZ1a1RP4XERSyOycunQ/ypOSqqJX+dajKHTsTLkzb61kfdfoZCaVKPpo3TAt4JwR2MXexJbdapdQbmPn6RYmnchxXF8Tci2gTgDADTiagjgEqTfQCgFsBjzNwLwAAADxFRr6A6gwH00P7uAzASAIioPYBhAPoD6AdgGBG10/Z5k5l7AjgNwDlENFgrvwdAATN3B/AOgNcsyBgzqjrOEoUrOYn2eXIqPYJbjyskL6I3Ji+XnNjJt7xxf3GYmoLTeH9rvFGkX/rTyWgkQdbizvnHd/QtbzlQolASgZmxt7ACWw+Uov8x7bH4qUvRrX1z1WK5DlPFkZmHAjgbQF9mrgFQBuBaC/vtY+YV2nIJgI0AugRVuxbAF+xhEYC2RNQZwOUApjJzPjMXAJgKYBAzlzPzTK3NagArAHTVtTVGWx4H4BKKg31ZZbL4SCmqkAnYbifap0nVU3jq838oOrKQrIjFMTUIFSVaiA/B6YZu63+0IklSm8ZpjXDvuccAAC4fPkexNKnND8t24+xXZ2DZzgIcc1gL1eK4lvRQG4joYmaeQUR/1pXpq/xk9SCa2+hpABYHbeoCIEe3vlsrC1Wub7MtgKsBvBvcFjPXElERgA4ADgXtdx881k0cddRRVk8hJImUc6xa8tK5nmj7zOIyKiQL8igLgvM8Pm6NahEEjeZaJFtmYFtuKY7r2FKxRKnJil0FvuWjOoilMRThLI4XaP+vNvi7yuoBiKglgB8BPKoF2YkZIkoH8A2AEcycHcm+zDyKmfsyc9+OHTua72CCRIYT7CRaBVA620KyIBbH1EF+P9UxfvVe3/KEh89VKIlQpQuW8cf6AwolSW2aN/Hb0jLSJS1NKEJaHJl5mPb/rmgbJ6LG8CiNY5nZyEK5B0A33XpXrWwPgAuDymfp1kcB2MrMww3a2q0plm0A5EUru1Xkd0+wk+hdVeVBFJID+aYmN00bN0JljaejvKewQuYQuQCxrqilqsavOL42eRMeuPA4hdKkLvllVb7lP50WPLNO8GIlj+OXRNRGt340EU23sB8B+BTARmZ+O0S18QDu0KKrDgBQxMz7AEwBMJCI2mlBcQZqZSCiF+FRCh81aGuItnwDgBkcB/+9RJrjKLifaB8n6WwLyYJYHJOb+U9c7Fse+I7M6VJB8BSbNs0aK5JEAIDqRArPn6TU1NXjl1UeK/wxh7VA+xZNTPZIXaxEVZ0HYDERXUFEf4UnUM1wk30A4BwAtwO4mIhWaX9XENH9RHS/VmcSgGwAWQA+BvAgADBzPoAXACzV/p5n5nwi6grgaQC9AKzQ2rxXa+tTAB2IKAvAPwEMtSBjzIirjWAn0Xaapa8tJAv/+WWdahEEB+nQMsO3XFFTJ0HbFPDejK2+5bOP66BQEgEA7jgrMDBRrSiScWfrgVLf8rj7z1IoifsJ6arqhZk/IqL1AGbCE2jmNGbeb2G/eQDCRjXVLIIPhdg2GsDooLLdodpk5koAN5rJZTdicRTsJOrgOOKqKghCAvLtkl342wXimhdPlu/0BwEJVlqE+NPziNYB62XVdWjTzIpdR7CLK0bM9S3rB7eEhlhxVb0dHgXuDgCfA5hERKc6LFfCIHqjYCc788qi2k8M34KQOBDRKaplUMmX9/TzLTfPMB2/FhxEfjvcwZi7/e/E14t3YeE2x0N0CEJUWBnSuB7Aucz8DTM/CeB++PMlCoJgI0N/WhvdjjKCIQiJxAdEtISIHtTHEEgVzuvREb88dA4AoLSyVrE0qU0ipRRLZi443h/l/7XJm3DLx4sUSiMIoTFVHJn5OmY+qFtfAqBfmF0EQYgz8tMvCIkDM58H4DZ4IoEvJ6KviegyxWLFlVO7tkFGeiN8Mjcb578+U3LRxonaunrM3epPby3BqNzDnWdnqhYhJVqHU10AACAASURBVFmzu1C1CAmFFVfVpkT0EBF9QESjiWg0gA/jIJsgCBaRIE2CkFgw81YA/wHwBDx5k0cQ0SYi+nOofYhoEBFtJqIsImoQAI6IMojoO237YiLK1Mpv0wWpW0VE9UTUx5kzswYRoUl6I+SVVWNXfrkvRYfgLJv2lwSsn9cj9nzWgj0MODYwUNHGfbakPhdMGDZ+vW/5kp6HK5QkMbDiqvolgCMAXA5gNjw5FUvC7iEIQlwRtVEQEgci6k1E7wDYCOBiAFcz84na8jsh9kkD8D6AwfBEFr+FiHoFVbsHQAEzd9faeQ0AmHksM/dh5j7wRDvfzsyrHDi1iNDnq6usqVMoSWrw+9p9+GBWlm+9aeNGknbARQzs1QnpjfzxH9+cslmhNMlPXT3j+QkbUFbld5d/52al42kJgRXFsTszPwOgjJnHALgSQH9nxRIEIRL+O2GDahEEQbDOewBWADiVmR9i5hUAwMx74bFCGtEPQBYzZzNzNYBvAVwbVOda+GMQjANwiZZTWc8t2r7KeUAXTXV/caVCSVKDB8auwKS1/qD4/9dfIqq6iUaNCPedf2zAuuAcS3fkY/T87diipeL4+8Xd0bqp5DQ1w4ri6E2yVEhEJwNoA0BsuYIgCIIQHT8z85fMXOEtIKJHAICZvwyxTxcAObr13VqZYR1mrgVQBCA4Ud9fAHwTSjAiuo+IlhHRstzcXCvnEjV6nXbwu3PD1BTspEvbZtjy4mA8feWJqkURgrjuNP8rvf1QdFHWBWuk6RTzTq0z8OilxyuUJnGwojiOIqJ28IyCjgewAZr7iyAIgiAIEXOHQdmdTh+UiPoDKGfmdaHqMPMoZu7LzH07dnR+/ttHt5/hWy6vlgirTrFDp4S0apqOJumN0NAYLajm+E6t0CTN0zXPOliKksoakz2EaNlT4Bu3Q2aHFgGKpBAa0wRKzPyJtjgHwLHh6gqCIAiCYAwR3QLgVgDHENF43aZWAPJNdt8DTxRWL121MqM6u4koHR4PIX1CuJsRxtqogqaN03zLz/yyHm/dJGminWDIZ0v8yxK909WkNSJAm/JbXl2HVuI+6QiPfuef5v3cNScplCSxkMy7giAIghAfFgDYB+AwAG/pyksArDHZdymAHkR0DDwK4s3wKKF6xgMYAmAhgBsAzGAtzwURNQJwE4DzYjwHW8lI9zs+ZR8qVShJcqMPAHJ0++YKJRHM0KdIqaiWoFHxoOcRrVSLkDCI4igIgiAIcYCZdwLYCeCsKPatJaKHAUwBkAZgNDOvJ6LnASxj5vEAPgXwJRFlwWPBvFnXxPkAcpg5O9bzsBN9J3nlrkKUVdWiRYZ0TewmvZFfQW/bXCKpuhl9as21e4qQeVgLdcKkCOK2bZ2wX2dthHIAMy+IkzyCIAiCkJQQ0TxmPpeIShCYRYcAMDO3Drc/M08CMCmo7FndciWAG0PsOwvAgChFd4z6oPSNewsr0KOTjP7bjT5qbbMmaWFqCqqp02mOf/9mJa4+9UiF0iQnBWXVqkVIWMIGx2HmenjyRgmCIAiCEAPMfK72vxUzt9b9tTJTGpOVAce2D1ivrZestHazOqcwYD2zg7iqupkRN5+mWoSk55aPF/mWH9TlkxXMsRJVdToRXW+QC0oQBEEQhAghogFE1Eq33kqLeJpypKc1wpFtmvrW9XPxBHvYkeePqPruzX3ELc/lXNm7c8A6swym2M2m/SW+5Vv6HaVQksTDiuL4NwA/AKgmomIiKiGiYoflEgRBEIRkZSQAfSSYMq0sJfn0zjN9yy9O3Ih//7BaoTTJx668ct/ytX2CU38Kbmf6xoOSqsZB9AG6BHNMr5bmQtOImRunukuNIAiCINgAsc6MoE0LSdmIMCd29ncpVuUU4ofluxVKk3y8NXWLahGECJn8qD/48b1fLMOj364KU1uIhSaiOEaEpatFRNcQ0Zva31VOC5UoiPuAIAiCEAXZRPQPImqs/T0CwFXRTuPNzw+erVqEpET6KYlJzyNaY/srV/jW/9hwQKE0yUVtXWBErsZpojhGgunVIqJXATwCYIP29wgRveK0YImAzOEXBEEQouB+AGfDk49xN4D+AO5TKpFigiOpfr80R5EkyUWpzBlNWIgIbZs39q1vPVASprZglR9XBHo0NGssUYYjwYqafQWAy5h5NDOPBjAIwJXOipUY1AbHERcEQUhCPrjtdNUiJBXMfJCZb2bmw5m5EzPfyswHVculkpYZ6XhiUE/f+uM/rlEoTfJwoLgKAJDWiDDiFonWmWi0bupXHHcXVCiUJDnYlluKZTsKfOvX9jkSjRpJsKhIsDqnoi08yYQBoI1DsiQclTWiOAqCkPzIz6o9ENHjzPw6Eb2HwDyOAABm/ocCsVxD/6DUHEJsFJXXYMr6/QCAif84Fz2PkPAUiUZzXc7NjMbiUhkrl7w1O2D98pOOUCRJ4mJFcXwFwEoimglP/+F8AEMdlSpBKKmsUS2CIDhOs8ZpqKipUy2GoBCJ3m8bG7X/y5RK4VIyO7RQLULSMHXDAfz1i2Vo27wxju/UUpTGBOWc7of5UkfIdNXY0M/37X9Me4y9tz/SZX5jxFiJqvoNgAEAfgLwI4CzmPk7pwVLBIoqUlNxbCqjXoKQUkjeN3tg5glElAbgFGYeE/ynWj7VtG/RBC0z/OPZ6/YUKZQmsVmUnQcAKCyvwfFB80eFxOHJwX737ds+WaxQksRn/V5/JsH8smpRGqMk5FUjop7a/9MBdIZnAv9uAEcS0WlEdHR8RHQvxRWpOek8vZG8bIKQSojaaB/MXAfgHNVyuJWBJ3XyLa/eXahQksQmPc3/1p4gimPCIsqNfVz13jzf8iOX9lAoSWITzlX1n/BEeXsrxPYORLSamW+3X6zEoDhFXVVlHnH8ufPsTHy+YIdqMYQURSyOtrOKiMYD+AFAmbeQmX9SJ5I7qJdw5bbwzeJdvuWjOjRXKIkQK78/ch4GvzsXgMdS1r5FE8USJT5nHN1OtQgJS0jFkZnv0/5fFKoOEf3hhFCJQnGKuqpKBKr4c/7xh4nimII0Inek/ZE33naaAsgDcLGujOGZEpLS1LngeU90mBnFlX6PqEEnSwCQREavKN744QJMf+xCdcIkCW2aNTavJBhiyQZORCcT0U1EdIf3DwCYeaCz4rmb045qh3vPPUa1GHEnXRRHJYy6/QzD8lO7SqDjZMUtlj6XiJFMfMLMd+n/AHyqWig3oLc4vjN1q0JJEpfaoNGmjHTJU5fIZKT7u+rbcsvC1BSsIrkbo8dUcSSiYQDe0/4uAvA6gGsclish6H54S9x+VupN9WwUp16kdFYDGRgibPRLfzolzpII8cItYzTyLtrOexbLUo4zM/0uZIdKqwIiIQrWqJQo2ElFU1FyYmazFpnWi1sGZRMRKxbHGwBcAmC/Nip6KiSXow9KQSeutDj1ZlPvykaHyu+ffHudJRW/L3oeuPA41SLYChGdRUSPAehIRP/U/T0HQHqHAIacnRmwXlCemlNCYmF+Vp5vWW+tEhKT4HtYJO9ExIyY7vdeeGJQzzA1BTOsfFEqmLkeQC0RtQZwEEA3s52IqBsRzSSiDUS0nogeMahDRDSCiLKIaI0WwdW7bQgRbdX+hujKXyKiHCIqDWrrTiLKJaJV2t+9Fs4tZlKx4xyvU5YRIWs4rVyEuw1pco+cxSWXV5UC2yH5gkA0AdASnvgCrXR/xfAM0qY8wd/9PQUViiRJXO7/arlvefpjFyiURLADIsLYe/v71sev2atQmsSE4fdcuCcFp5jZiRXFcRkRtQXwMYDlAFYAWGhhv1oAjzFzL3jyQD5ERL2C6gwG0EP7uw/ASAAgovYAhgHoD6AfgGFE5PVfmaCVGfEdM/fR/j6xIKPgYlzSZ27Af648UbUIAajU3SRQkrO45uoqEiTZBo+YeTYz/xfAAGb+r+7vbWaWCX0aH+nmc//z+1VYsO2QQmkSm67tJKJqMnBO98N8y/X1jOzc0jC1BT319YxJa/cDAJqkN0ITscLHhOnVY+YHmbmQmT8EcBmAIZrLqtl++5h5hbZcAmAjgC5B1a4F8AV7WASgLRF1BnA5gKnMnM/MBQCmAhiktbWImfdFcI6OkmT9Glfhlmt7XZ8jA9ZP7dZWkSTGqLxOEihJEKIig4hGEdEfRDTD+6daKLdw+UlHYPOLg9C1XTNsPViKWz+WxOeC8MtDnvSvw8avx8VvzVYsTeJQUF7tW66tq1coSXJgNapqbyK6BsDpALoT0Z8jOQgRZQI4DUDw178LgBzd+m6tLFS5GddrLq/jiMjQnZaI7iOiZUS0LDc31+IZhCbZRsTdhFvmdx3VoUXAeryCA+kJdy1UXqdEclVd/9/LI6r/0p9OdkgS66h41oxQJUUSj0v8AGAlgP8A+LfuT9DISE8LGKTbcUiiSQqpzald26Bxmv+jKIGjrHGo1K84uiG9VaJjJarqaACjAVwP4Grt7yqrByCilgB+BPAoMxdHKacVJgDIZObe8FgoxxhVYuZRzNyXmft27Ngx5oMmSr/mxwfOsmSeH2IhSmzclGWXXtx4BQeyitO3I1zzJ3dJnDhZLTJCpq015Lb+6iMmu0RvVDZA5pLTd4JaZh7JzEuYebn3T7VQbqNL22a+5f/NzFIoSWKwOqcQi7PzzCsKCQkRoY9uMCUnX+b/WiG3pMq3fF6Pw8LUFKxgxeI4QFO0huhyTt1tpXEiagyP0jiWmY0SG+9BYKCdrlpZqPKQMHMeM3ufjk8AGCe9S1E6tmyK9s3NA01cf0ZX9MtsH7ZOi4zA4H9OKVIu0898uM3KplKaPke1xXNXB09dFpINVc9YEnt0TCCiB4moMxG19/6pFspt6KNJSnTQ8OwprMC178/HX0Yt8pW9f+vpYfYQEpFHLjnet3z+GzMVSpI4bD/kmQ/65T39AuZPC9Fh5Uu80CCojSnk+cX/FMBGZn47RLXxAO7QoqsOAFCkzV+cAmAgEbXTguIM1MrCHa+zbvUaeOZUOk6i9Gusykkg015i+xZNMP7hc/Dn0z3ew04peG5xVQ2mkQv6L2/c0BsdW2UAiO8zeFT7hoEW2rfMiJ8ALmPCw+di9bMDHWvfLW+Aqu9c5zZNAbgvIJUNDIHHNXUBPEHnlgNYplQiF6LPXycBLcJTUikpGlKBoztIsCOr7C4ox0ezt2HK+gPofnhLnNv9MDRvEpnnkdAQK1fwC3iUx/0AquDpy7DmEhqOcwDcDmAtEa3Syp4CcBQ8DXwIYBKAKwBkASgHcJe2LZ+IXgCwVNvveWbOBwAieh3ArQCaE9FuAJ8w83MA/qHNw6wFkA/gTgvnFjNuVW6MsNL5s9pB7N21LTLS07R9CID9juNuVcpVzjv74LbTMWntPtzYtxvem+F13YqfPNJ58/PmjafilK4eV90ubZthT6H9bkOJYnHr060tVuUU2t7uZb064at7+uPs4zrY3rZKmFniwVug/zF+I2xwAm8hEH1fpP8x7dHnqLYYeFInhRIJTuAdMBbM+esXy7Fxn2eG3I1ndE2Y31O3Y0Vx/BSaAgjAcjgiZp4Hkx4te2b2PhRim3duZXD54wAeNyh/EsCTVuWzi2R8Dq2eknditttcN50mXorj/Rcchw9nb/OsaIe84pTOuOKUzgH1VF7+1Lrzfs7p3gHXn+6P1+VUkAK3XF+zAbJfHjoHmUMn2n9cIpybhHNSiKg5gH8COIqZ7yOiHgBOYObfTPYbBOBdAGnwDJq+GrQ9A57B3jMA5AH4CzPv0Lb1BvARgNbw/JafycyVtp6YzfTVTZtYsC0P1bX1MngVgjpd1I/C8ho8OTjprPQCAq3wgOe+uy3uglsoq6r1LR/fqZVCSZILK1/gXGYez8zbmXmn989xyRKERHpdrcgaiRJSr3WWnfpmuSWiZDDx+kZfqVMQTzqydch6TosTbpSOzD2bE5qP7+hrWH7piZ3iM3rpkosbz1fxtetPid/B1PEZgGoAZ2vrewC8GG4HIkoD8D48+Y97AbjFYBrJPQAKmLk7gHcAvKbtmw7gKwD3M/NJAC4EkBC+jXeenelb1ncEhUBmbDrgW66qrVMoieA0PY/wK0H/myFBo0KhTxfWo1NLhZIkF1YUx5VE9DUR3UJEf/b+OS5ZouCSjp0ZRNbd3jKCRrQatKWdtHeA0ykFz62XNm5BZbXjdG7TFIe3ahqmnrorlUiu2tFwWS9rrl5O3QO3XN3/3959x0lRn38A/zzXG3dcoRx3B0fvvVcpckdRUQQbIlZijcaowWgs+NOgMdZYE4klGrsJsSFGjIIV6YhIEQVEeu8Hz++PndubnZ3dnd2duvu8ed2L2dnZmWdnZ3e+35nv9/naGUdFcgxY3pKZ74NSeWPmg4i8m/sAWMPM65j5KICX4RsLWW0c6jKKvw5ghJJvoArAUmZeomxvBzN7onZxx2kd/dMvf70hzJLJ7f4PvvdPN1FloxWJ5+rhrfzTC3/a5WAk7rV172GsUw3h07409MV3ER0jFcds+Po2ViGG4TgSXaIVnAmE+yd0weUntYy47Aml5phiXXYcV0qzKTuO8YRG9pFxo3zs2udOXRQY310zbK6dYbj0e2+yo0SUDaVzOBG1hO88G46R8Y39yzBzDYA9AIoBtAHARDSbiBYSUVB3j1pmj3Vspnvf/87pEDzh4XO6Ox2CsJC6WbKckfV9v2V/wONG+aEvvovoRCwBq4bgUP8ZGo4jGdhRrsvNCH8H0AijBVAioGF+FqaNbhdxWaubqjpVfpw/bXjY59NSCRVF1l/RNXpRwtE+jpSY/XzdwqmuKw+c3c2ZDSPxLsaFcDuA9wFUENGLAP4Lnb77JkoDMAjAJOX/M4hohN6CZo91LOwnCVQS2wnVBdzP1253MBL3Ul/fl0y05orq1gkRLbQqEK+yo4hTmBt5/EWzRJN5tfany7Kmqg7VSMoMNPPpUl43CO8bVwwIs2Tsan/43HSTT/uZEJKmoB/ArmPTLVngtJ/xQxZWLF3yli3FzHMAjIcv+/c/AfRi5o8jvMzI+Mb+ZZR+jQXwJcnZCOATZt6uNIt9F4AM8ieEB51Qpak8dpyx68BR54JxKXW5NE2SB5kq2jZ3svc17CjYmbEJMrgeI5WA2vX4+zha9KWM5X3nZ9k/Rk/3ivqRF4qBWypk6iiYGXecGvWwrlErybPvYonWVcNaosjGizVeoP0umvmzN6Zz44DHbk2KZSYiOgNADTO/o2RSrSGi0yO87GsArYmoORFlADgHvrGQ1WbBN0YkAEwA8JGSvXw2gM5ElKNUKE8C8K1Z78dOSywY9sXrlm6UfZJMqjsF/mZu3uPq5MiOeG3BRv+09Pk1V7QVR/PzrXuctogzfVxH3eXi24Z9BSmvZ1XNNqFZr1uk+O/suuiWI4DR6uFAEuxuMwDcWN0OC/8wMuwytiVIivD8xQPtGQ5QG0e/FuaNq3hen2aB20r8eiMA3M7Me2ofMPNu+JqvhqT0WbwavkrgSgCvMvMKIpqujGEM+IbPKiaiNfAN9zFNee0uAA/AV/lcDGAhM3vmfL78zmr/tPRzDHbO01/4p6f0bxZmSZEI8jLTkKkalmbMI5/izCc+czAi93ljYV3FUfr8miuqiiMz32pVIF6lLeRkWjDGlCl3HA32RYtmU8PaNgQAtG1sTbYqN5cf1bHF8/n8akiL0Ntw4Q6wrYmmLVuJnX3JcUI/t37GWFw+NPj4seJCjvZzb5SfhYxUc37rgu5mmrJW19PbeRGbSzDzu8zchplbMvPdyrzbmHmWMn2YmScycytm7sPM61Sv/Qczd2TmTsp4yJ6Rl5mGyf18FaKs9FQ8M+8HSdSlcvBoXYLcy4dGTmwnvE+dIAcAvvlRsquGIi2IzBXyzE9E85T/9xHRXtXfPiLaa1+I7uaW5oROmNCzHMvvrEabhtaMjxNLHcWsz+P8fk2NbzOOytRQpfIdYs0AQvdxdOJOJDMHVpot2o4bK83hWBdv9Ct2S7/IWHk8fKMWENEDRNRS+XsAwDdOB+Vmt4z1DWj/0Xdbcdfb32LhT9I8U09+VrrTIQgb3HtmF6dDcC25qGStkBVHZh6k/F+PmfNVf/WYWQZEqeWRMcBJ+RdxuSg3lpdpf59CO/xuVOSssmYIt7+ND8dhzkFYXhh9PwCrsqq6/oJMlG+6MsasbpE2Y9d+srIyF7xql3/25rgGwFEAr8A3HuNhAFc5GpHLZaWnonNZgf+xFA6DVXVohNwEPSeLQGf2LMf6GWOdDsOV1HfghfkMtTUiokIi6kJEPWr/rA7MK4KbWUVf6ImU0MWs06Oxwp+ByqVthVXvFCDP7FHun+7VrNDw68K9w2je/bhuTdC8JDeKV+hsL9QGSb0M2VKu1za3bBHne3NarMdyLK+y4uOx9CO3MPGOWzHzAWaepgx70ZuZf8/MByK/MrmNUiUF+e93Wx2MxD2++6WuAdiRmhNhlhSJ7p9f/eR0CI5iZsxfsx3PfrbeP29oWxlSyGwRK45EdBeApQAeBfBn5e9+i+PyDG0ZJ5bmg1ZlJVUzfPcqluahCTaOYyzum1DXbCSapD5GPnsjR9TD53TH3BuGGt6uHiMXBLRX+a26iBBU0XLZweCGPo6hnrfkDrCBdT56bl0CgjO6a8elD7Nuzd502UdtCSKaQ0T1VY8LiWi2kzF5wWRV4pcnPl7rYCTu8boqe2TPKC5aisRw75md/dMz3kvuxFFvLdqESX/7En+avco/TxLjmM9Im4azALRkZhkoRocZd8UircHOgpSbCm2xZFU1q9Ac7eeaqq4ARvFSN+3vVIMXMNQFfdvuDiVBq7TelYX4en1ggoNYKua+15i9wyLHMahViX+6Y5N8vLVIO8SgwS0lwy1HoETJpArAl/WUiMJ1eBaQ/nt6dhw4ivLCbDx/cR80K/Z2ywwRvYqium4Q6alJ8dsZ0oadh4LmWZGwMtkZ2aPLAVgzUF0C8MrXlGAs1lgKbZFeE2s/SLvKj/XM7hMSRZk93HsszstEvcw03KokhYjGye1jK4OunzEWk/oaTwxkF7fVG6Md19DIoVycm4kPrz8p/vU69KNk2kUbc1bjdieIyP9FI6JmcN9hLlxuz8FjeGvRJuw/UoMWDfIMX/wTiUN9MSXZP3/1OeiULqW4dWx7ZKUnzhBtbmGk4vhHAIuIaDYRzar9szowrzBnqAz39Bm0IpJYkxiEi6V+TjpO7doktoA0btMZ0N6+poj6W3p8Ug8UZKdj2Z3VGNctfLM/vebRj00ypxty1/KCoHnqkPUuSAxuXYJ4ubylqm2c7OP41pUD6tZp4weQmZ4UV4hvATCPiF4gon8A+ATAzQ7H5DmHkjgJxvT/fIuu0z8AAOw+eMzhaIRT1HfUtuw94mAkzlPXm9uX5uPSwaGHOxOxM3KGfg7AvQBmoK6P45+tDMpLgvvnRF/CiuUi0VXDWqIwJ4pmO3H0cVxye1XEZawQrqK7+LaqgD5VVm1HT1pKhK9NHE1Vz+3TFO/8ehDGdC418NrQG8pMi/4qm97aXrysX/gEPjpPpps0vl/ghsxfZTys69upfexcc2311WtDrRVi3Ce1Fz7aNa6HP47vjHYWjQvrMrMB3AqgHXxZVQcDkIHYDFCf9371j+QdwWTm/B/80xeo+n6K5NKwXlbA4yUbkneYGnVSnGS+qGQ1IyW8g8z8CDPPZeb/1f5ZHplHaAtpsY2tF77ApVd4vLG6HRbdVqWzdKxbCZab4at8FGSHr6BG7KNp4+0Ks7YULuTGBVmhn4x6O4Ebys9OQ8cmwXf5nJKqsyPUc6z6bN3cza1bRX2M6dw4YJ5ZFclo16K3/62o1Br6nGPdrPKTWZiTgXP7uK+ZtEUeB9AXQB4zvw1gH4DHnA3JG/591SD/9Cffb3MwEve4eGBzp0MQDinISUfXirreZBt2HXQwGmdt31+XiiXJW+1aykjnrk+J6I8AZgHw3wdn5oWWRSUCmDFeldHCpHq5eb8bjgNHa+LedsxNVe26sxnFsr0rzc1aV5ybYer6zEbkuwv61Cfr7N2u9lNxSe+v/Kw0/OuqgbZtz8ksxwFJkGzcbpLoy8w9iGgR4E+O4+4fA5fIkGQXQfIjXNwViS1b1bz/hEvOlXbTljOnntTSoUgSn5Ff4O4A+gG4BzIcRxBz+jgGz3tlar+YkqKE31B0sRTmZqC8MLaBy81gX5NYnXk2tI3MSE0JyIhmxnabFefghqo2sb1Y2bR2f/xuVLvAxSJ8MDIwtw7DTcV9C6YQML57GUZ2aBT1amM5gtx2p69d43pOh2C1Y0SUCuWSCBE1ACCD8BkgWRKD1YswFrRIbH8cXzcc2IkkrTkeOx74vmNNyigii/gLzMzDdP6G2xGcF5hRwdC7pd63RbE/tbYpQ36Q/d3EXp7az/CyZlXerG4Wa2aFsk3jvLhe36zYV+lU92f8343DcPXw1nGtVyvSWJO2jBvo8btZRk5i6vfcOD8LD5zdLba+qjF8IHrdUgOSIFnfUjXAv68eiBV3Vse4Rk94BMBbABoS0d0A5sF3cVZEkK6pOM5bvd2hSJyzZuv+gMeW9CsXntG8pG4YluteWexgJM45UiN9Gu0S8deGiBoR0TNE9J7yuAMRXWJ9aN4QnP0xhsqOi0rFhgqIBgumrRvm+Zdfffdo3DIm/B3Uu8/oFHUsZrBrEHWz/eW8HvjrBb1M7XMJRD4eA/s46jxvxc7z+EXUVg3z8OT5PdC+NHzil6D6cqThOHQvuJivyECT6ng/d/XLM9NSkZvAV4yZ+UUAN8GXtXwzgNOZ+TVno/KGbE16/dVb9zkUiXP+8cWP/umezcztPiG8LxnvOn65bqfTISQNI5epnoUvA1zt2AffA7jOqoC8RltUiiU5jh2VFILB4ThMDKZ2XcyM9NSUiO9zUt/AzHC2DYlh8pbiWVs0u78gOz1iU8Z4RR5H0AM1bBsY+dxG8e3yPAAAIABJREFUdSpFXmZ0dxC1q22cb+AigckfSWZaCsrqZ1u22cIcX6W09kJTsmDm75j5MWb+CzOvdDoer0hNIZzerW4opmfm/YAzHp/vYET2S1O1AhnYstjBSIQbbduffMNyXPr8Av+0NnmdMJeRimMJM78Kpf8FM9cAkHvCClOakUY530oxbTPEi+KNP8W28S3NXZ/RSwde7Qpo14WO8DMSU7h9+9UtIzD3hqGBy5u0Ywa1ahDyuZYN8ixtAt6hST5euqwvbhkbPJ6qEHpqVHdUNu46hEU/JdcQBH+bVzcUxyWDZKw6Abx2eX9cMsiXXXfFz3scjsZZj0/q6XQICc1IxfEAERWjrhN/PwDJfVSqmFGcClUoi6ZeMfPCXhG38esRkfu+WdNfLbaVdmiSFOO5uU5A3zadIzwg2yb55ljO5ZVsq/aA+rPIz0pHdkbkO5axxDKqk/ErtLedol/Bi+e3Y0DLEsmWKQwbqzPObbIk5dL25SqIZjxnkbB6Vxbhxuq2yExLwbzVO7B172GnQxIJysiZ+nr4huJoSUTzATwP4BpLo/IQR+6+6BjeLnKTxdO6Nom4jJlq903sw3EQOkZZeTRtKIKIfctCLxBXU9U4XusEq+K1KtmDVWM7RTrCYxkOx+ALgmc51EE38IKC145k4SWjO5di3T1jApJO1SRJv65t+5KvGaIwJis9Fb0rizBz/g/oc89/8f7yX5wOyRbqMqY2V4Ywn5GsqgsBnARgAIBfAejIzEutDswrtAUkM/t8mVn0MrouI/EHJ7zUf028+yKWV8dSR42qkGtBeXjFndX4zckxDqFhItL8D4SoQGvmGen/Fq28rDS8eGlfVJnch/N/Nw4zdX1GGen7bFZmYacGPpa6orBTSgph/5G6cYanqvo4JbKtUnEUYbRVDWW0aMMuByOxj7qCnObUCTCJGMmqOhFANjOvAHA6gFeIqIflkXlI/xaJ0zndSOFveLuGBlcWxXaNL2o6/bpRfBH9dqTximBuZprucAheQAR0Li8ImGdWk7GBrUpQYPLA1tpxM41YcWc1PpsWfgQis45f7fevaXH08Xrpbt95fd01fqTwlpcu7eufnrtqm4ORWI+ZceBIDTbvrmuC+Nh5UhQTgdQZsJOlEvXqgg3+abPLDCKYkeLqH5h5HxENAjACwDMAnrA2LG8pjXM4hFDlPCca3kT6mVlyWxUu6N8swlLOiaXM3LmsIPJCtes3uNxZvSuiDwTuu2ujW6lWzaxNYJRroO+d6YHYJDczDTlWvz+N2srfxJ7lqnl6y5m/7b9f2Dvs82Zus9RIllghQmiZRJl4n/5kHTrePhvPfbbeP29sl+C+nkLUSk3x6BXpKO04cBSDW5fgyfN7oLqjZFS1mpGjqrYn9lgAf2XmdwBEHNSLiCqIaC4RfUtEK4joWp1liIgeIaI1RLRUfSeTiKYQ0Wrlb4pq/t1EtIGI9mvWlUlEryjr+pKIKg28N1ewpZ+k8baqYRXkpAc3zw2VVdXE9zW4dQmqO1oz9ERhTgYqinzNLc/v57sDEm9l3mX1P0NqP9do7lhph1CJ9vWG6ez4f17WD09MsueKe8SxLS36Esey3ngjGRahRYH/hrJmQ2676CESX72sxB3rU+vtpZsBAF+t941X93+nS18uESxTlWTsiY/XOBiJPVZu3oulG/cgIzUFozqVeqrFjVcZqThuIqKnAJwN4F0iyjT4uhoAv2XmDgD6AbiKiLTp+EYDaK38TYVyJ5OIigDcDqAvgD4Abiei2lFu/6PM07oEwC5mbgXgQQD3GojR1ZwZjsOdfTR/N6odnpocPnNszFSBmpbaPMo374aEgHpNTPV+hNVzjGT5tEr/lsUY1LrE1m06+THpZ7jVWc6iH44LB1Ras2IhYpSd7tzvj93UrQ7vOr0Tzu/n3pY/wjmTVS3Cjh3nhM42PO4v8zD64U8BAPPXbnc4muRhpAJ4FoDZAKqZeTeAIgA3RnoRM29WEuuAmfcBWAmgTLPYOADPs88XAOoTUSmAagBzmHknM+8CMAfAKGVdXzDzZp1NjgPwnDL9OoARZNelBw9c4DCc1dFF74Uo+gpVrPFrtxPP2Jrx7EMzK+4JRbNb/nqB7yJCol1djDcTsbIWU2IJufZQLQxi2G7iFmmEHbTf/7cWbXQoEuulqGqO0o9LhJKZFngx5djxxP2VXbKxbmRAGQfYPkayqh5k5jeZebXyeDMzfxDNRpRmo90BfKl5qgzABtXjjcq8UPPD8b+GmWvgG2syKGsNEU0logVEtGDbNnd0pndTZSGWSKway/KULvYMHxI4bmH8mI19pm67EBjp7mK45Syn2Vd9mhf5YrE5DCu3F7k5rN48nc/MqiFHIhywCVaHFx5xcvu6ptVzv3PHOd0Ki37a7Z/OT6ImuiI+y39OjmHXz5dEa7axvOcsEeUBeAPAdcy81+rtRcLMTzNzL2bu1aBBA1PXff/Erqauz8+UmpkJ64h2k3GUJNfPGIuRJg/FEA1t7JZ9tgqX1SFdI/J4mvbEEa9ohyzx4p1U70XsPUQ0iohWKX35p+k8r9vXn4gqiegQES1W/p60O3arDGpV11z9hNuuxplkw86DAY8zvJqGW9hi1tUD/dPjH//MwUjs48VzpldZ+utDROnwVRpfZOY3dRbZBECdfrJcmRdqfjj+1xBRGoACADtii9w8X9w8IuIykY73WE4Sg2Ps+2Xml8+JtvVxjx2p8/LJ/ZqhSf2skM/rrSPa3dhYyS7ZKM4MvdG6ZnirsM9r38f9E7vaWkEIdQg5dY4IdUSHCufaEa3x8Dnd8OsRrePetoEhNS0V7W/Dw+d0syiS5EREqQAegy83QAcA5+rkDQjX138tM3dT/i63JWgbDG1bd8cxUSuO6vEqASBfmqqKMLqU13c6BJHALKs4Kv0LnwGwkpkfCLHYLAAXKNlV+wHYo/RfnA2giogKlaQ4Vcq8cGYBqM2+OgHAR2xzzUVvc2mpsRftatdWVj8bf57YFd0qjP0Y5GWmoZNmiAmjZb6Ymqq66EJPLLFEKpA3qJdp6JbgfWd2iX7jiom9yvHU5J6Y1Mfe5ha/rWrrnzay67o3LYy8UJTO6lUeeSED/eoW/mGkSRGZ5zcj22BctzKkKxd/wl7YcNH3SE/kpqqBbyA1ScYQs1EfAGuYeR0zHwXwMnx9+9Wc6+vvkMqSXP/0u8t+Qc3xEw5GYw3tJ6g9vwuRbN5bppfqRNjByjuOAwFMBjBc1TxmDBFdTkS1VzvfBbAOwBoAfwVwJQAw804AdwH4WvmbrswDEd1HRBsB5BDRRiK6Q1nXMwCKiWgNgOsBBDXjsUq4wqBZZ+wze5ajODfiKCgR2TVG298v6u3JpgPx3LEsL6xrjhjtWogI1R0bByRAcEpAn08b+tD96qSWhmIJmK/zfJEJ349I4n3rHEejZDv7M4aNw8B854/ihGSk73+4vv7NiWgREf2PiAaH2ogb8wBEY/Oew06HYLoaVYKTkjzrf+dEYknEzKp3vf2tfzpSyylhLst6WDPzPEQoPyh3BK8K8dxMADN15t8E4Cad+YcBTIwpWJfS7rzaQuL47pHyBOm8Vvk/0hhtZiXqGda2YVDzGjsY+X3MSEvB0ZrwV6XjLZB7sdJshtO7NcHew8Y/96cm90TLBgYG8U68855jOjbJx4qfo+9uHukjSNJD3is2A2jKzDuIqCeAfxFRR728A8z8NICnAaBXr16e++Zt2HUQFUU5TodhqsPHjvunk/XcImL35sJN6N+yGE2i7GvvZuqL7FUdGjsYSfKRHtYWM/IjH+1pYHTn0ojLxHy2j6mpZ/wnMrNOhZF2d6P8TFyuucNl2onY4+fzpyb3NLRcuM/71lPMTYmdn+Xry5OlGa8tUQtP4S58GH3HkZbrXVlkNBzhLkb6/uv29WfmI8y8AwCY+RsAawG0sTximzw1uSfaNa4HAFj1yz6HozHfhCc/90/nZUpGVRHZI+d290//9rUlGDDjI2zbd8TBiMzDzNi465D/cWObc0MkO6k4WswtxVujBW2vl8drw79yaEsMbRucNZdgbHDIUPsrbCWZ1ctF5rZL+S2UO3+xXAiwqiL326q2uHVse5zaVX9oFtccrw7FEemzmnvD0KB5qSmENo0M3OUN2laE513zYSSsrwG0JqLmRJQB4Bz4+var6fb1J6IGSnIdEFELAK3h6yaSEKo7Nsb71w1BZXEOPl29Hat+2Rdwly6RROpuIkQoew4ddToEU3z5w07/dJfyAl8eCmEbqThazFgWTvMLXLGu0cxI0pSmBANbxZbhNRa1lbHGBVl49qI+yEqPfIibPXalm8bltEK4w9Xsd56dkYpLB7dAqK6f4fZ1pD7B0WQezlSOo/E9DCTyCSOeY8PIz0Rtdt5azVWJQ2pde3JsGV6tuNCRgF1vLKP0WbwavkRxKwG8yswriGg6EZ2mLBaqr/8QAEuJaDF8SXMur80bkEhOatMAH323FdUPfYJb3lrudDimOHEi8EtSqfOdFkJrVMfGuHBAJRrl11WqEuXi3u6DdRVgGZrGftLmwQZZ6Sk4fCz6TG8tG/ruCoxoH/94hvH+XJTVz8aQNvrjXob6LcpKT8XcG4ai1GXNCGIpq0bzGga7505YDMyO/dax7fF/76w0dZ2RQjyzRzkWbdiFHQdCX2G97uQ2+HT1dkPby0pPxfI7q5GtaTJrJyPJccZ0LsXbS38GANx9Rifd9eRnpcdVYTNa+JBKofmY+V34ksqp592mmtbt68/Mb8A3NFZCa1ea75/+5sfEqBff8NoSp0MQHpSRloI7TuuI95bXZR89fiIxfpT3qfIoHEvALMpuJ1V1ixFIt7mYWn6Wfv29eUkult9ZjXP7VATMtzJDVqhC4fxpw/HH8Z11n9P2P3vw7K7+6eYluUHPWyneOo9ZA84n+l3HaFw6uIV/ukt5Aa4cGjqLai3t3gt11zrWSq5eM+ZI8jLTQg4x4ZZPW5219WQTLjgJ4SXqVgaJMhzMm4siDWEtRGjqmxZVD37iYCTxO3CkBht3HQzInDypXzMHI0pOUnE0kW51joC0lNC7uapDI/8Yb3ryMtPCXuEvzEnHMAOFYCvHcbxkUPOAAc7HGEje4xS9/RBN5cNInd3rlcZ4oo90N6phvSxcPKh51OvtVFaA9TPGqraj/B/1moD1M8bi2Yv6xPDK+IUbjsPOoyYxrjsLEag4r65ZXrjzrld9etMwp0MQHjMoii4ZbnfO019g0L1z8fzn6/3zzupVEXJ5YY3E+2V1G4tLaItuq8LfTSwEx3IHJys9FdePdEeCvksHt0CbRnkYG6byGssN277Ni3Bun6b408Quxl7g4bpjpAxlcQ9VEt/LA9fl5TbBFjO7YYI0PRVup844umrLPuw5eMzBaOI3c94P/ulx3ZoEjBUshBF/mmCwzOIByzbtAQBs3+/rgtK/RXG4xYVFpOJoA3XZdt7vnLliaPQumDlDa4RozmdDIb9ZcQ4++M1JAVee1X41pIXu/EjSUlPwx/GdUV4YenywRClXXzQw+juCteyqxnn1rq7Zccf6lYqnubvU1YVbabMFf7s5+vFK3WLr3sOYrhrk/KGzu8mFMhG1nIzETGXy8tR++OfUfk6HkZSk4mgxbdM0bcWDAfRoVmjqNuNramhaGEGs7JtZt43wz184sHnQZ2JFJcTL5/dIfYP09pdTb9fM7Y5o1xATesaXNTUZhPuOXTig0rY4hNAiIrRoUJd11Mv9HO99f5V/enK/ZlJpFEKldoxnYT+pOFosUkWGGbipuq3BtVmbyt8ssWwr1Gu0u69FiFTkhTnpmNS3KQa0lKYLVsnJiJzkKNTnqE6YZAZ/H0cTj+tLBjfH/RNji9OpQp1uJd5AKNFcwpHyqvCSVNUB++nqbQ5GEp83Fm70TyfqmJTCfurmz15y8GhNwOP87MS8k+oFUnE0QbiCFSNSdY+R5qJxaNxeSJx1zSB8+fsRQfMbF2Tj7jM6G9qXepV5M2+GEhmr4ttxB9YsK6ePwsI/jAQQ2zESqrlM/FlwXX7AaoRLjuN2ncp8Qx20bJAXYUkhnHPTqHb+6Uc/WuNgJObJNnDRTggj/jR7VeSFXOjv89cHPM7PljuOTpEquzCdlUX5vMy0gAQIZrCi7uG1Ck0kRgsuIfu3Rrm9aKpX7UvzcflJsfVdNUusn/bTk3ti6gvf+NZh8Jj56pYR/v1s52F2Vq8K9GhaiNaN6kVc1mhcXq5IC3ca2SHxhqG50XCrJCGCvXXlAJzx+GcAvDuWo3a8xrwE7bvpBbLnLcbMYQuEbrvp5NWkI5dGMcSDy3a558RyhFi5z9+7drCFazfG6PvTfr+qOjaOelsN64XOessMZKYZqOTH8IEQkaFKY20cQojYaAv39aQ/l4hDg3p1yQKPaipgXqE+d5bkZSDFw/2Xvc49bSQT0KS+TVGkGpBYz2+r7LmSaG8fR2Mbe+OKAaZt88wQSU2sLsA+dHY3/3TPZoUY3q4hpo/r6NHqtwlseuMJdkPXVM9d3AfXndwajfL1MwsLkUw27T7kdAhRO1rjzcK9cKdEuJCnrifGcsFVmEcqjmbSfDl/P6Z92ErU+hlj0aFJfrybMZWdBfKWDfQT3XhVVnoqZl7YG60aGrsr41WJ1gzXDFP6N3Nku3qfRPOSXFx3cpvwLR2sCwmAVOyFe5zx2Hw8/OFqHKnxToKZ1Vv3OR2CSCDa7j0fr9rqUCSx23ekLjmO3Gx0llQcTRDpGHbDMW58HEczthVifhKVJpPorcYl1nT5lgyhEsc6J/evxPoZY/2PHzq7Gz68/iQzwvKsBLjILTxq/rTh/umt+47gwQ+/xwuf/+hgRNE57S/z/dO9K80drkskn8LcDFQW1w0F94wHM6s+/ck6/7RXu1QlCqk4WsitlYe3rxmEly7tGzDv1K5N8PzFfSzN8KrNIqr+8tv5Q6BttuHWz8mtYtldoV4ztG1DAL7mlWZz8mM9vXsZWjWMPfuo0WMy1osxXsroK0S08rOC0zcc8Wjzz5en9nc6BJEAzu3T1D+988BRByOJ3t7DxwIeTx3ibDK8ZCcVxyShLl92KivAgFYlAc8X52ZgSJsGNkflXZHK626/Ivboud3x9jWDbNse60wX5Wb47ziepDn2It7Fj7T/3b37dXkwZCFcKScjDS1KcpGZVlfE8cpvgnq8ujtP6xhzqwwh1KYOaYGvlKHMVvy8F+Mem4/t+484HJUx320ObLpdUZQTYklhB6k4WsjulPk9mta3Z0MRGH6/UewXr5z0a7k93lO7NkGnsoKYXqv73lz+fhNVrLvdzPuN/7txqIlrEyJ+qSmEj24YivP71fU9fvjD1Q5GZNwj/60be9JL/TKFuxERGubXZeResmE3Rj30qYMRGcPMOOupz50OQ6hIxdFDagvsoVqZPas092MEF+7dWK53qnKlHTvO7XcHEwGFmDZjfaEkW2PM8T3KcO2I1rZvtzivLntrtE1gpcWssJJ67DevNFVdv/2Af7rGo2PuCfeqUo1z6oU7jtrvbb8WRQ5FImpJxdFEQRUSpXTrRMXEyQJZqH5XiZQcJ5kLvLF8jkm8uwwjAhrlZ6GsfjbuPK2j4dfUeuCsbiiMMPxPLSPHr9FjPHG+1SLRVHXwTtp+ZgYz45e9h9GsOAeT+jbFlP6VToclEszvRrdzOoSoLN+0J+DxC5f0DbGksEtwD3LheW4tyAUnx3GIJMdxVG2fnUhjnIaTSBch1DLSUgIyQgohYjeodWBf/uMn2LV9Bkc//Cm++8XXl6tpUQ7uPqOzwxGJRJSb4Z1i/7HjJzDhybpmqnee1hHpFiZwFMZ45whKMBlp0R/80ZzugpqqurCg7caYjKqfkw7AN3ajHg+/tfiEuEul3h31czIwY3xnnNTWXcmYvPiZxZxV1cR7wF7cbyI5Ha05gewM/d9sp9VWGgHgp50HHYxEJLLcTHce/3que2Wxf/rG6raYMqDSuWCEn1Qc7aApWM35zRAUKBWPZBK2kGtj4dOMIvO00e3QrDg3oL+AnhQCHjm3O65+aZEJW/WeB87qCiB4n5+jSg0eCyOHi9Rn7Ofli0Ei8X29fqcnsoffUNXG6RBEgsrR3HFcs3UfWjWs51A0oe05eAzvLN3sf9yhNN/BaISa3PN1QOtG9dCwXlbkBU3kxuKcW2KKJY6cjDRcMqg5UiI0e0pNIZzSpUlsgSWA8T3KnQ4hpDSXNlmzg5l9dNV9uGV8SOFmF8z8yukQdKmT+ADAlUNbORSJSHTaptovfvmTQ5GEd0A1LA0gGYbdRCqOFrL74nt6iu/j7FgmV2bCScbCbX5WGs7pXeF0GKYx47vVpbwAN1a3RYuS3PhXloQifQS137OT2zfEv64aaH1AQugY6rIm8XpmvPddwONIFySFiMdpXesuZu8/XBNmSecc12QUbtPIfXdFk5VUHG1gVgVyYi9fwb9zuf74e9kZqXjjiv54+oJelsUQj6DkOC6IyWrhMuraWX8tLcjGjDO72LdBDyAiXDWsFUrqZUZeOMHYeexlpaeiW4U7xpgVyeeZKb2dDiGir37Y6XQIIoncUNXWP/3aNxtdeTFdPQxHVnoKWjTIczAaoSYVRxM0UgZVLcgO7Ldo9jAcIzs0wvoZY1FWPzvkMj2bFSE/Kx3dKwpN3bYVnBo/UfsbmQz9ssxMhhILs/bwG1cMwJT+zSIvKGyj9/VJ/G+UM4hoFBGtIqI1RDRN5/lMInpFef5LIqrUPN+UiPYT0Q12xew0t2ZRVVN/h0ZG6DcvRLyaFucEPN57yH13HdVNU9NSpKriJpZ9GkRUQURziehbIlpBRNfqLENE9IhykltKRD1Uz00hotXK3xTV/J5EtEx5zSOklPqJ6A4i2kREi5W/MVa9N61fj2iNB87qiuqO+mNGOXHaOrlDI3x1y4i6GFxQOQoXQ+0zZ/eqwINnd7UnIAtlpKagT2URHpvkO6SHt2vocEQ+FUWhLzpYyaxqa89mhbhzXCdDx7P7rqHqq5fpS1bg/DfUfl75jNyAiFIBPAZgNIAOAM4log6axS4BsIuZWwF4EMC9mucfAPCe1bG62fdb9kVeyGbq37Pz4kwcJkS0tu477HQIQdZtO+CfduMd0WRmZTW+BsBvmbkDgH4ArtI5yY0G0Fr5mwrgCQAgoiIAtwPoC6APgNuJqPYW2hMALlO9bpRqfQ8yczfl711r3lawjLQUjO9RHlSYdbquZncCnmjp7Z97J3TBGd3rEqpY8YOhXaMVHxMR4dXL+/uvHj81uSeW3VFlwZai8/61Q7Dg1pNNX29hjjImI7mjAnRKl1LjCzt8Tpp+ekdnAzBR/5a+cfM6lwc3TX3u4j64T5pLx6MPgDXMvI6ZjwJ4GcA4zTLjADynTL8OYITq4urpAH4AsMKmeF3p1n8tdzqEIAE3Rd3wAyqSypa9R7B80x4crTkReWGbXPPPukz0Um10F8sqjsy8mZkXKtP7AKwEUKZZbByA59nnCwD1iagUQDWAOcy8k5l3AZgDYJTyXD4zf8G+GsXzAE636j0IYZb01BTUy3J+CJbczDSU5Jnfp+/FS/ti+riOKMhO1/2Rt7ssNLRtQyl/RenuMzphTGf9VhNGjerUGEvvqELPZnVN5WsvqJ3UpgGqO8W3/iRXBmCD6vFGBJ9T/cswcw2APQCKiSgPwO8A3GlDnK5zx6l116y37zviqjsYzIxFP+32P05x+oqzSAovXtrXP/3c5+txyqPzcPc73zoXUBj/vKyf0yEIFVsaDiv9LLoD+FLzVKgTYbj5G3Xm17paafI6U3WHUhvLVCJaQEQLtm3bFsO7Ma72598NzUST1W2ndohq2IVk/Kg+vWkY5k8bHtc6KopycEH/ypDP21lM6xIieZRbuaUMO6lvMzw+qWfQ/Gi/E/lhLpAk4/fLJe6Ar0XO/kgL2nmOtMupXZugU1k+6uekY932A2h+87t4+St3DEOwUFVpBIDG+e5uKSQSw8BWJVhz92iM7NAIc77dAgBYsnGPw1H5qJvOntO7Al0luZqrWF5xVK50vgHgOmbea+GmngDQEkA3AJsB/FlvIWZ+mpl7MXOvBg3cn6Y7kREB43toL5iba1LfZlhzT2B31+DkOJaG4HoVRTlhEy55ybfTq/H65QOcDkMY1K95kdMheMkmAOoxdcqVebrLEFEagAIAO+Dr9nEfEa0HcB2A3xPR1XobScRzZHFeJt6+ZjD2HDrmnzftzWUORuSzZe9hnPnEZwHz2jaWYQeEPdJSU9CvRbH/8Y87DoRZ2j6XPLvAP33ZkBYORiL0WFpxJKJ0+CqNLzLzmzqLhDoRhptfrjMfzLyFmY8z8wkAf4WvP4ij5E5jZHef3tnpEAAA9bJ8CUq0zYSSeZB4M9m1F3My0pCRJhnYInFDU72V00dhQKsSp8Pwkq8BtCai5kSUAeAcALM0y8wCUJtMbgKAj5SuIIOZuZKZKwE8BOAeZv6LXYG7xRM6d9SdcvjYcSz8cVfAvES5gCe8Q32He9fBY9iy1/lEOcs21d35dMGpSmikWbVipUP+MwBWMvMDIRabBV/z0pfhuyK6h5k3E9FsAPeomptWAbiZmXcS0V4i6gdfs9cLADyqbK+UmTcry58BwDU94JOl6vGHUzrgLx+tjum1kerYZg7doTc0xTMX9sZ7yzajierE/fcLe6NVQ2vGDqoszsFtp2pzRYlkcf/ErmjTyHdsFeb6EguVFVpTaJzYsxyfrd2BTbsPWbL+aKi/xdkZqY7F4UXMXKPcJZwNIBXATGZeQUTTASxg5lnwnXNfIKI1AHbCV7kUiuqO7hnq4op/fIO5qwKbAkeV1EsIE+RlBVYD9hw65h9izm73vv8dKjVDhRQr50fhHpZVHAEMBDAZwDIiWqzM+z2ApgDAzE8CeBfAGABrABwEcJHy3E4iugu+K6wAMJ2Za0fIvRJDDyeLAAAaQklEQVTAswCy4UsrXpta/D4i6gZfd6r1AH5l1RszKlkqjLUuGdQclwxq7nQYUSMilNXPxqWDA5tEDLNwCI2Pbxxm2bqF+5UXZqOLknl0aJsGePL8HhjRXr9Q26phHtZsjdg1LaQ/TfQNb1M57R3/PLmI601KtvB3NfNuU00fBjAxwjrusCQ4D3BTKyBtpREAbhrVzoFIRDLr2CQ/4LGT35AnPl7rn75oYCVuGdMeaanSgshtLKs4MvM8RDgGlcyoV4V4biaAmTrzFwDopDN/cmyRWs9F5yrPsaLZoTR9MEdRgl0JtPNrGpB9nwijOoW+0/DmlQOwY/9R64PSqJ/jS3LTRWdoDQBoX5qPxRt2h/x9G9WpMc7uVYHfVrfxz3NTwV2InQeOuuZ3bHDrEqRK1whhs5K8TORlpmH/kRoAwHGXFJBqjrNUGl3KyjuOyUn1uy9lpPg9dl4PvPjlT2hfKgkD3OTZi3qHTOLglsO+eUku1m5zR2d/rWhOzflZ6WEzlca0fQMBVBTl4O1rBqFNI/3P+bmL+mDlL3uRmabf5DQzLRX3TpBxG4V7XfXiQvxzqv2p/muOB46XN/u6IZIURzhGfb3CqbEc1ZlUAeC6k1s7EoeITKrzZnPHxRpD1s8Yiz+c4lw/u1QDNeuKohxMG91O7lS4zNC2DVFaoN8nzy1fgQfP7oaZF/ZyOgxP61RWEPKuf0FOekBGPiPkWyzcZMOug45s91+Lfw54nJ0u/X2Fc1JUNccL//51mCWtM/aReQGPiy0Yb1qYQyqONjAzsYvZarMr2lkvm9S3KT6/ebg0Q0hwvx9T11/Hicpkvax0DG/nnmQYbqKXIEqIZPDU5LrMqht3HcLug/Y3Az987HjA46wMORcK56gv4u88YP/3AQC27TviyHZF9OTXyiITepbLXbIQ0lNT/HernCjAumEogmTQosSajLRudvHA5jg5RJIbN7nj1I4oybO/b5f8JAqnDdQMAePEoOfavoxZcsdROKiqY2OnQxAeIhVHYTspPNpjyW1VeOOKATG99rLBzVGQbV6/umg/cq9edLnt1A7425TIzWOdfnejO5diwa0jHY5CCPvlZaahV7NC/2Ntf0M7aHPg5EjFUTho+riOAY9/3OFsboA3r4yt3CLsIRVHs+mUCD1aBraM2Tf8MmXAd10FOenIy4wt/9UtYztgye1VMb1W73A38pE/em533DSqbUzb9Jpkveft5mb7InkcU1UW1+84iEl/+wKvfr3Btu2nqAoFYzuXSrcN4ah0zfF30p8+xgUzv8Lew8cciadpUU7khYRj5Ncqwbkl1XgksRYol91RhYV/iO7OyW9GtsGEnuUxbU9EFmul6NSuTdCnssjUWIyQ/n5CJJfGBXUDnN8/exXmr9mBm95Yatv2X/9mY90DuZYiXOiT77fhP0t+jrygCXZp+lVKsih3k4pjgnv7mkF4xkDTOadkpaWiS3kBHjqnW0yvr5eVjtwo76rVz8nA/cqg6MIeRstGLRv4+kVePLC5dcGE4NXmsV4iu1i4wX1ndvX3MzykSlRz4oT1F5H2H6nBlz/s9D+Wr4Rwg0fP7R407/O1Oyzf7vETjDOf/CxgnvT5dTepOJrNZTcvmtTPxggXJ+tISSHMunoQqqVzdsKIpyBUmJuB9TPGYmyXUtPiEUIItYKcdKy9Zwz+fmHvgPmvLrC+uWqn22cHPD69W5nl2xQikp6qfr+13l662dJt7jpwFH/+YBXWacZb1iaPEu4SWwcoIeIgdx0Sm8uunQghhK6KosCxaK0eEkCb0fvTm4ahQvpzCRcoLcjCr4e3whfrduKr9Tsjv8AEve7+EMc1d/mdaG0koiN3HM0myXEiktEwRNKS3wIhXKOsfmClzcpz9TtLN+Pxj9cGzJNKo3ALIsL1VW3RSNX/12raSiMA5GRIM1W3kzuOAoBkOxTJq3bYkfRUG74DctFECNfI1hRSrernfPjYcVz10kL/40b5mbqFZiGcdu2IVgFJcZjZku/FYVXfYrWrhrUyfVvCXFJxFEIktT9N6Io3F21Ct4r6ToeS8KT1hXCb8T3K8ObCTQCA77fss2QbRzVjRU4b3Q5ndJfM3sJ9WjWsF/D42HFGRpr5P9x/fHelzrbzgi7mCPeRpqpm07mIaORu3uDWJSirnx1xOS+TMmNyiHUcR6cU5mbgkkHNJauqEEkoVfW9//dia4YfOHLsRNjHQriV+k65mdZqEuLMnzYcb105wJJtCXPJHUeXeOGSvgCAymnv2Lrd2kxag1oXW74tN1ceEp2d/UrlcxahSJN44TbaDI7Hjp8IGhA9XnNXbQ143EIZdkgIt5vz7RZL1qu9TpvoN04SiVQczWZCcpwMk09a4XRvWojv7hol4+YIS0l1QWHRjijJy7BmxUIkuBRNxXHvoWPIyUhDVnqKaa0Qbnp9acDjPs2LTFmvEHY4cYKDvifxSpEWPp4lTVVd5snze2DO9UNs3aZUGoWZ5HRgr9cv7493fz3Y6TAMkbKCcJtLBgWm/5/y96/Q/rb38fLX1o/pKIQbXa1JUPPiVz/hu1/24oftB0K8InoyVKN3ScXRZUZ1KkWz4lynw7CUFB6FME+vyiI0zLcvhboQiaSlptno8k17AQA3v7ksaNxFM5zUpoHp6xTCTDdUt0W9zLoGics27saohz7FsPs/NmX9J04w5q7a5n/88DndTFmvsIdUHG0g9aRAMo6jSFpJfuzLb6Hwki17j8T82kufW4B73l2JQ0cDhx148vye8YYlhOWmDmnhn351wUZT173KouzFwh7Sx1EIYSq31Y3mTxseVHgTQohafxzfGTe/uczUdX64cguwEnhF0+RVhhsQXpCbqV89WPXLPrRtXE/3uXCYGa99sxFVHRph9MOfxhuecJDccbSBpPkXwjll9bPRqqFLshgm+U+B/BYKNzq3T1Pd+ceOxz9sxp5Dx+JehxB2m9y/Ge44tQPeuKJ/wPwbX18S0/rmrtqKm15fijE6lcYWJS45PwtD5I6jsI0UGZODfM5CCK8Z26UU2empeP2bumZ5R2piqzj+7dN1uvN/P6ZdTOsTwm7pqSm4cGDzoH6+sV7823e4BgDw857DAfPfvmYQOpUVxBakcITccRS2cVsTRiGEEAIAHjuvB+6f2DVg3lP/WxvTuv7vnZW68we3lsQ4wlu0FcXaRxt2HoyqC8iyjXt050ul0Xuk4mgDuQMjhACQ9FdP5LcwPkQ0iohWEdEaIpqm83wmEb2iPP8lEVUq8/sQ0WLlbwkRnWF37F70mnL38cCRGuzYH3uinFrtS/PjXocQTlq8YTf2HDyGwffNxdQXFhh+3d/m/RDweMb4zvjrBb3MDk/YQCqOJhvRriHK6mcHZKQSItkleX1JiLgRUSqAxwCMBtABwLlE1EGz2CUAdjFzKwAPArhXmb8cQC9m7gZgFICniEi6qujQji/3w/YDqH7oE/T8vw9jXmdRbgaemNQjzsiEcIY2E/Cabb6sqJ+u3h62D+8Nry3Bp6u36T53evcyjOzQyLwghW2k4miy4rxMzJ82HG0aRZ91SgizSS4Sl0nyz0OOx7j0AbCGmdcx81EALwMYp1lmHIDnlOnXAYwgImLmg8xco8zPglzLCSlVU3O84h/fYOOuQwCA5Zv0m9tFsvAPIzG6c2ncsQnhBO3FlBteW+qfnrXkZ93XHKk5jte/2YjJz3wV9NzQtg2QlS7Zhb1KKo42kMKSj+yG5CWfvc/oTo0BABVFOQ5HIjyoDIB6bIeNyjzdZZSK4h4AxQBARH2JaAWAZQAuV1UkhUp5YeB3c922A/7pUx6dB8BXgXxtwQZ/1tX/rtyCymnvYNnGPbjqpYX2BSuEDbQXU37YfiDEknWG3/8///Sin3YFPHf8hFy38jKpOArbyE9FcpDPObQLB1RixZ3VKKuf7XQojpDhOJzDzF8yc0cAvQHcTERZessR0VQiWkBEC7Zt029mlsheuqxvwOOjOkNynPLoPNz4+lLMeO87AMC/Fvvuupz6l3l4Z+lm64MUwkZGEth8tnY7rn5poT8L66bdh/zPnfH4ZwHL/n5Me3MDFLayrOJIRBVENJeIviWiFUR0rc4yRESPKB35lxJRD9VzU4hotfI3RTW/JxEtU17zCCklESIqIqI5yvJziKjQqvcWLSksCSEA329BqIGVhYhgE4AK1eNyZZ7uMkofxgIAO9QLMPNKAPsBdNLbCDM/zcy9mLlXgwbJlwW0tCAb/7l6EACgYb3MoOcPH6vLJPndL3sjrm/Ob4aYF5wQDmiUn4VG+cHfBQDYuPMgAGDKzK/w9tLNhoawkSRR3mblHccaAL9l5g4A+gG4Sqcj/2gArZW/qQCeAHyVQAC3A+gLX7+O21UVwScAXKZ63Shl/jQA/2Xm1gD+qzwWLiLV5+Qgn7MQlvgaQGsiak5EGQDOATBLs8wsALUXWicA+IiZWXlNGgAQUTMA7QCstyds78nN9PW/6lZRP+i5qgc/CZoX7jevgU7lUwivubFafwzSpz4JHLP0yLHYxj4V3mFZxZGZNzPzQmV6H4CVCO6PMQ7A8+zzBYD6RFQKoBrAHGbeycy7AMwBMEp5Lp+Zv2Df/fDnAZyuWldtUoDnVPNdo0/zIqdDEEII4UFKn8SrAcyG73z6KjOvIKLpRHSastgzAIqJaA2A61F3AXUQgCVEtBjAWwCuZObt9r4D72jRIA+PndcDfz6ra9BzPyl3WABg/podOH6CQyYIAXwDqQvhdRN6lod87tjxEzh23NdEtev0D/DNj7tCLqt3MUZ4iy1tppSxpLoD+FLzVKjO/uHmb9SZDwCNmLm2c8EvAHTz/BLRVPjubqJp06bRvZE4fPCbIWiSpP2ahHOKcjMAACPaS9rrZJeWQpjYK/TJX7gfM78L4F3NvNtU04cBTNR53QsAXrA8wAQytosvC+obV/THmU98HnK5Zz9bH3Y9ORmSPVIkho5N8pFChO+37Atoktr6lvcCljvzic+0L/V768oBlsUn7GF5xZGI8gC8AeA6Zo7cIcAEStMc3RwdzPw0gKcBoFevXrbl8ZDhOYQTSvIy8fUtJ/srkCJ5rblnTNC85y7ug0NHj+ssLYQAgIrC8BmQ73r726B5eZlp2H+kBl/cPEJyHIiE8c6vBwMALn72a3z03daoX9+jaX35PiQASyuORJQOX6XxRWZ+U2eRUJ39NwEYqpn/sTK/XGd5ANhCRKXMvFlp0hr9US0SVrvG9XCCkzPfp/SxEaGc1Cb5kp8IEY1ok1m1aZSHD35zkkXRCOG8lCgrfykEDGxVgvsnBjf9Ft5jWcVRyXb6DICVzPxAiMVmAbiaiF6GLxHOHqXiNxvAPaqEOFUAbmbmnUS0l4j6wdfs9QIAj6rWNQXADOX/f1vyxoQnvX+dZLYTQggRnWibmk7q28yiSIRwh2i77c6fNhylBdJVK1FY2Wt7IIDJAIYT0WLlbwwRXU5ElyvLvAtgHYA1AP4K4EoAYOadAO6CL4vc1wCmK/OgLPM35TVrAdQ2rp4BYCQRrQZwsvJYuEhBdjoASLPJBJeZ7vtZ0Q4aLIQQXkNEAYntbj9VmxxeiORyQf/KqJaXSmNiseyOIzPPQ4TM/Epm1KtCPDcTwEyd+QugM/4UM+8AMCKmYIUtzuhehpoTJzC+hyToSGR3jeuEZsW5GKJqBjm4dQOc1KYBbhkrA/8KIbxl5oW90en22QCA7PTAO5AdSvPx7ea69A2SNVIkuoGtStCrWSEWhMie+tJlfXHeX7W5MEWikDzRwjYpKYSzezd1TXryejIQuyWK8zLxu1HtAu44Zmek4rmL+6B5Sa6DkQkhRPTyVOeKNM35a9bVAzG8XUN0LivAFzePQFepOIok8Mqv+uPG6rZ488oBGNY2sK98gzzJq5DIpOQsktYH1w/BD9sOOB2GEEIIl/v4hqHYtPsQDh8LzEKclpqCmRf2digqIZyRmkK4algrAMAT5/dEuz+873+uKDcDd5/RCbNXbMHjk3o4FaKwiFQcRdIqLci2ve19agqhQ2m+rdsUwm3GdG7sdAhCRKWyJBeVJblgZmSkpuDo8RMBdyKFSFZZ6alYP2Ms5ny7BWkphOK8TEzq20wSRSUo+dUTwkZrdcbSEyKZLLm9SgZFF55FRFh51yhs23fEn/BNCAGM7NDI6RCEDaTiKIQQwjZS2BZel5pCaFyQ5XQYQghhO3dkKRFCCCGEEEII4VpScRRCCCGEEEIIEZZUHIUQQgghhBBChCUVRyGEEEIIIYQQYUnFUQghhBBCCCFEWFJxFEIIIYQQQggRllQchRBCCCGEEEKEJRVHIYQQQgghhBBhScVRCCGEEEIIIURYUnEUQgghhBBCCBEWMbPTMTiGiLYB+NHpOAwqAbDd6SBcRvZJMNkngWR/BEvmfdKMmRs4HYRXeOgcmczHdCiyT4LJPgkm+yRYsu4TQ+fHpK44egkRLWDmXk7H4SayT4LJPgkk+yOY7BORaOSYDib7JJjsk2CyT4LJPglPmqoKIYQQQgghhAhLKo5CCCGEEEIIIcKSiqN3PO10AC4k+ySY7JNAsj+CyT4RiUaO6WCyT4LJPgkm+ySY7JMwpI+jEEIIIYQQQoiw5I6jEEIIIYQQQoiwpOIohBBCCCGEECIsqTjaiIhmEtFWIlqumldERHOIaLXyf6Eyn4joESJaQ0RLiaiH6jVTlOVXE9EU1fyeRLRMec0jRET2vsPoEVEFEc0lom+JaAURXavMT9r9QkRZRPQVES1R9smdyvzmRPSl8j5eIaIMZX6m8niN8nylal03K/NXEVG1av4oZd4aIppm93uMBRGlEtEiInpbeZzs+2O9clwvJqIFyryk/d4I75NzZCA5PwaT82Noco4MJOdIizCz/Nn0B2AIgB4Alqvm3QdgmjI9DcC9yvQYAO8BIAD9AHypzC8CsE75v1CZLlSe+0pZlpTXjnb6PRvYJ6UAeijT9QB8D6BDMu8XJc48ZTodwJdK/K8COEeZ/ySAK5TpKwE8qUyfA+AVZboDgCUAMgE0B7AWQKrytxZACwAZyjIdnH7fBvbL9QBeAvC28jjZ98d6ACWaeUn7vZE/7/9BzpHa/SHnx+B9IufH0PtGzpGB+2M95Bxp+p/ccbQRM38CYKdm9jgAzynTzwE4XTX/efb5AkB9IioFUA1gDjPvZOZdAOYAGKU8l8/MX7DviH5etS7XYubNzLxQmd4HYCWAMiTxflHe237lYbryxwCGA3hdma/dJ7X76nUAI5QrX+MAvMzMR5j5BwBrAPRR/tYw8zpmPgrgZWVZ1yKicgBjAfxNeUxI4v0RRtJ+b4T3yTkykJwfg8n5UZ+cIw1L2u+OWaTi6LxGzLxZmf4FQCNlugzABtVyG5V54eZv1JnvGUpzie7wXUFM6v2iNDlZDGArfD9UawHsZuYaZRH1+/C/d+X5PQCKEf2+crOHANwE4ITyuBjJvT8AX2HpAyL6hoimKvOS+nsjEpIc05Dzo5qcH3XJOTKYnCMtkOZ0AKIOMzMRJeX4KESUB+ANANcx8151U/Fk3C/MfBxANyKqD+AtAO0cDskxRHQKgK3M/A0RDXU6HhcZxMybiKghgDlE9J36yWT83ojElqzHtJwfA8n5MZCcI0OSc6QF5I6j87Yot7yh/L9Vmb8JQIVquXJlXrj55TrzXY+I0uE7Kb7IzG8qs5N+vwAAM+8GMBdAf/iaTtRe7FG/D/97V54vALAD0e8rtxoI4DQiWg9fE5nhAB5G8u4PAAAzb1L+3wpf4akP5HsjEk9SH9NyfgxNzo9+co7UIedIa0jF0XmzANRmaZoC4N+q+RcomZ76Adij3F6fDaCKiAqVbFBVAGYrz+0lon5KW/ULVOtyLSXWZwCsZOYHVE8l7X4hogbKlVQQUTaAkfD1bZkLYIKymHaf1O6rCQA+UtrczwJwDvkyqDUH0Bq+ztxfA2hNvoxrGfB1jp9l/TuLDTPfzMzlzFwJX6wfMfMkJOn+AAAiyiWierXT8B3vy5HE3xuRsJL2mJbzYzA5PwaTc2QwOUdaiF2QoSdZ/gD8E8BmAMfgaw99CXztyv8LYDWADwEUKcsSgMfga7u/DEAv1Xouhq/T8hoAF6nm94Lvi7EWwF8AkNPv2cA+GQRfO/SlABYrf2OSeb8A6AJgkbJPlgO4TZnfAr4f8TUAXgOQqczPUh6vUZ5voVrXLcr7XgVVxi9lH3+vPHeL0+85in0zFHUZ45J2fyjvfYnyt6I25mT+3sif9/8g50jt/pDzY/A+kfNj+P0zFHKOrH3vco604I+UNy+EEEIIIYQQQuiSpqpCCCGEEEIIIcKSiqMQQgghhBBCiLCk4iiEEEIIIYQQIiypOAohhBBCCCGECEsqjkIIIYQQQgghwpKKoxAJgIjqE9GVynQTInrd6ZiEEEIIN5BzpBDmkOE4hEgARFQJ39hNnRwORQghhHAVOUcKYY40pwMQQphiBoCWRLQYvoFt2zNzJyK6EMDpAHIBtAZwP4AMAJMBHAEwhpl3ElFL+Aa/bQDgIIDLmPk7+9+GEEIIYTo5RwphAmmqKkRimAZgLTN3A3Cj5rlOAMYD6A3gbgAHmbk7gM8BXKAs8zSAa5i5J4AbADxuS9RCCCGE9eQcKYQJ5I6jEIlvLjPvA7CPiPYA+I8yfxmALkSUB2AAgNeIqPY1mfaHKYQQQthOzpFCGCQVRyES3xHV9AnV4xPw/QakANitXIkVQgghkomcI4UwSJqqCpEY9gGoF8sLmXkvgB+IaCIAkE9XM4MTQgghHCTnSCFMIBVHIRIAM+8AMJ+IlgP4UwyrmATgEiJaAmAFgHFmxieEEEI4Rc6RQphDhuMQQgghhBBCCBGW3HEUQgghhBBCCBGWVByFEEIIIYQQQoQlFUchhBBCCCGEEGFJxVEIIYQQQgghRFhScRRCCCGEEEIIEZZUHIUQQgghhBBChCUVRyGEEEIIIYQQYf0/p/ou06UMO6sAAAAASUVORK5CYII=",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10896bb70>"
      ]
     },
     "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 [Close Encounters](../CloseEncounters).\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)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}