{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Additional forces\n",
    "REBOUND is a gravitational N-body integrator. But you can also use it to integrate systems with additional, non-gravitational forces.\n",
    "\n",
    "This tutorial gives you a very quick overview of how that works. Implementing additional forces in python as below will typically be a factor of a few slower than a C implementation.  For a library that has C implementations for several commonly used additional effects (with everything callable from Python), see [REBOUNDx](https://github.com/dtamayo/reboundx).\n",
    "\n",
    "**Stark problem**\n",
    "\n",
    "We'll start by adding two particles, the Sun and an Earth-like planet, to REBOUND."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "import rebound\n",
    "sim = rebound.Simulation()\n",
    "sim.integrator = \"whfast\"\n",
    "sim.add(m=1.)\n",
    "sim.add(m=1e-6,a=1.)\n",
    "sim.move_to_com()  # Moves to the center of momentum frame"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could integrate this system, and the planet would go around the star at a fixed orbit with $a=1$ forever. Let's add an additional constant force that acts on the planet and points in one direction $F_x = m\\cdot c$, where $m$ is the planet's mass and $c$ a constant. This is called the Stark problem. In python, we can describe this with the following function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "ps = sim.particles\n",
    "c = 0.01\n",
    "def starkForce(reb_sim):\n",
    "    ps[1].ax += c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we need to tell REBOUND about this function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "sim.additional_forces = starkForce"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can just integrate as usual. Let's keep track of the eccentricity as we integrate as it will change due to the additional force."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "Nout = 1000\n",
    "es = np.zeros(Nout)\n",
    "times = np.linspace(0.,100.*2.*np.pi,Nout)\n",
    "for i, time in enumerate(times):\n",
    "    sim.integrate(time, exact_finish_time=0) # integrate to the nearest timestep so WHFast's timestep stays constant\n",
    "    es[i] = sim.particles[1].e "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let's plot the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA24AAAE4CAYAAAAjGaCkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XmgndPVx/HvMrU1lVJD0ZqJMeYiSMSQkBAqbWOsMVQ6\n0vltaatUX5RWDTXP8yzGRG5EEPMsZkq0qKl9VQnW+8fe1Su5ueM5Zz3D7/NXb9z33t9LzjnP2nvt\ntc3dERERERERkeKaLTqAiIiIiIiIdE6Fm4iIiIiISMGpcBMRERERESk4FW4iIiIiIiIFp8JNRERE\nRESk4FS4iYiIiIiIFFyXhZuZnW5mr5jZw518zx/M7Ckze9DM1mpsRBERERERkXrrzo7bGcCQWf1D\nM9sGWN7dVwD2A05sUDYRERERERGhG4Wbu08C3uzkW7YDzsrfOwVYwMwWbUw8ERERERERacQZtyWA\nF9t9/RKwZAN+roiIiIiIiNC44SQ2w9feoJ8rIiIiIiJSe3M04GdMA5Zq9/WS+c8+wcxUzImIiIiI\nSK25+4ybXt3SiMLtamAMcKGZfRl4y91f6egbextSisWMzwEHAN8C7gOOh7k2dH//5z34GUsCuwOj\ngZeB84A7gPvctWNbJGZ2qLsfGp1DmqM7/33z6/UgYA/gUuCPwCPdfa2aYcAGwG7AV4HHgUnApe7c\n3/v00hm9dqutmf99zdgM+B3wGeAE4FpgWg9e8wsAI4CvAxsBjwATgOPcebUZmatEr91q68tmVneu\nA7gAuB1YycxeNLO9zGy0mY0GcPfrgGfN7GngZOCbvQ0jxWXGHGYcZMYDwF+AFYAt3NnGnetg+oc9\n+XnuvOTO4cCywG+BNYELgSlmbN7o/CLSM2YsbMbBZkwAHgY+AlZ3Zz93Hu7JAos77s6d7hxIOhf9\nu/yPrjHjAjMWa/z/ByLSE2YsYsbXzGgDzgSOA/q7c1L+zO7Ja/4td850ZwipE+snwILAVDO+Y6Z7\nhEV6o8sdN3cf1Y3vGdOYOFJEZiwCXAR8ABwI3O3O+4342e58CFwFXJXfyHcEzjbjHODn7nzQiN8j\nIt2Td8f2Ao4grbIfDdzizr8a8fPze8e1wLVmHA78DLjHjFHuTGrE7xCR7jNjXtLrfVdS58vJpN3w\n6Y34+e78A5gITDTjGOAcYFsz9nDnr434HSJ1oRUP6ZQZ6wP3AJOBIe5MnkXR1tbX3+XOR+5cCqwF\nrA3cYMZCff250mdt0QGkqdr+8z/MmJ20yv59YLA7e7lzbaOKthm58447PwX2AS4148c5gzRGW3QA\naaq2vv4AM5YjHXmYH1g2d9Fc0KiibUbuPA1sQurkus+MYXmxSD6pLTqAFJO5t+Y4kZm5zriVixn7\nAr8B9nPnyhb/7tlJK4A7ATu7c2crf79I3ZgxD3A+MC/wFXfeavHv/yJwNjAXqXi82p13W5lBpE7M\nWBO4DjjMnRMDfv8mwOnAv0mtmX9oVsEoUiR9qYlUuMlM8i7bT4EVgR3ceSIwy0jg96RzdUe7c1lU\nFpGqysNHrgYeBEY3qhW6FzlmI7Vr7QYsD+zkzr0RWUSqzIzBpIWab7lzcWCO2YENgUNI52hH5tZK\nkcpS4SYNYcbnSateq5HOtZzarBapnjBjTmAIcAxwvjuHBEcSqQQzlgV+AIwkDQz536JMdTVjJ9I0\nuxHu3B6dR6QKZpgK/XX3YrTkmTEHcDywDqlNW8WbVFZfaiKdcRMAzFiUNJ77YWB5d/5QhKINwJ3p\n7lwDbAyMNOMn0ZlEys6MjYE7gVeBddz5XVGKNoB83nV34EozVo/OI1J2ZnwHeAZYGdi4KEUbQB5E\ndgBwN3C1GZ8JjiRSSNpxE8yYC7gFmOjOz6LzdMaMxUmDUg515+zoPCJlZMYGwDXAru7cFJ2nM2aM\nAo4ENnRnWnQekTLKC567AsPdeTY6z6zkdulzgfmAHXXmTapIrZLSJ2acBCxGepP8KDpPV8xYhTRx\naaQ7E4PjiJRK3l2/B/hm3skuPDN+TLrId1O1UIn0jBl7kq7d2KQM4/fz8YgrgLeB3crwXCLSEyrc\npNfM+CbpbrYNy/RA1O5g9aaRw1NEyiSfI7kZmOTOL6LzdFceF34SsAJp4uWbwZFESsGMNYDxpM/K\nx6PzdFdulbwemA78jztTgiOJNIzOuEmvmDEc+DmwXZmKNgB3xgM/AW4yY3PdAyPSLf8LvAf8MjpI\nT+SzdwcCDwBTzfh+bqkSkVkw47PApcB3y1S0AeSrQLYi5b/CjDHBkUQKQTtuNZTbEPYCDgOGlXkl\ny4wRpIfRR0nndf4vOJJIIZnxfdJF1wPceSM6T2+ZsRpwMvCwO/tH5xEporywcTkwzZ0Do/P0hRlf\nIl3Yvbc7N0TnEekrtUpKt7W7r+kfpPtbHg6O1Gd5uMopwMKkg9fqhxfJ8kCfI4FNSe1SfwmO1Gdm\nzAvcBfxWQ4pEPikXbUcBXwYGRt3L2Ej5su5LgfXdeSE6j0hfqFVSusWMeUg945cBg6pQtAHkD6W9\ngQVJd9OICB8P8rkbeBlYvQpFG0DeWf8acLQZK0TnESkKM9YFppDuQxtWhaINwJ1JwLHA6WqTljrT\njluN5OmRcwN7FOm+pkYxY3nSvVSblK2fX6TRzJgfuA84zJ0zg+M0hRnfBXYgLURpp11qrd01HwcD\n51Ttcz4PV5oMnO3On6LziPSWWiWlS2YMAC4EVinbIJKeMONAYGdS8aYHOaktM44D5nNnr+gszWLG\n7KQHuTPcOTk6j0gUMxYAHiQdgbg6Ok+zmLEycBuwgTvPROcR6Q0VbtKpPIzkfuCX7lwSnaeZcgvF\nbaTVxhOj84hEMGMdYCywqjuvR+dppjysZALQXxd0S12Z8UdgLndGR2dptjxoaQTp/J4WaKV0VLhJ\np/LltQOBoVVrneiIGauSLujWg5zUTt6FuhP4U1VbJGdkxqHARqT3uA+D44i0lBn9gRtJHTWVXqiB\nj9/jbgaeBMa480FwJJEe0XASmSUzliP1ux9Qh6INwJ1HgRNIh5hnj84j0mIHAP8CzooO0kKHkT7P\nzs5DmERqId9hejzw8zoUbQB5cWYEsDQw2YxlYxOJtI4KtwrLB3lPAY5057noPC32a8BIF3cuFB1G\npBXMWAo4BBhdl4UagLziPhyYDtyfF6xE6mB34FPAadFBWimf1d8GuAS4yYwFgyOJtIRaJSvKjLlJ\nl9QuRmofql0rQb7f7XDgK6Spc8/HJhJpnny+82rgLnd+FZ0nihljgDHAuvnaAJFKMmNR4CHSZ/x9\n0Xmi5InZH7nzzegsIt2hM27yCWYsBtwEPExaea/1w4sZ3wP2IF3cWYk7bUTaywOI/gCsDgx2573g\nSKHMOAf4qzs/jM4i0gx5oeZi4Gl3fhydJ1LebXsc2KbOBayUhwo3+Vh+gJtMumj70Dq1S81KPgNw\nDXCbO7+NziPSSPmhZSzwOrC7O28GRwpnxhdIC1fr1rBNXCrOjE+RzrAuAWzlzrvBkcKZsTewJ+kq\noNo/90ixaTiJtPc94A1UtH0s/3v4PnCwzrtJleRFifNIF21vr6Itcedl4FjQQo1UixnzkiZIzgFs\nqaLtY2cCCwDbBucQaSrtuFVIvoDzKWCAO09E5ykaM04E3nHn4OgsIo1gxq6kRYn163iOtTP5nO9T\nwHbu3BudR6QRzDgDmBPYQ1dffJIZ25MGk/XX/W5SZGqVFODju4yWducbwVEKyYzFgUeB1XW/m5Rd\nHr7zJLCrO7dF5ykiMw4ARrizdXQWkb4yYzNSi+Sq7rwTnadocgfC7cDx7pwXnUdkVlS4CWZ8jvQQ\nt747z0bnKSozjgI+7c6Y6CwifWHGPsDX3NkyOktR5TO/jwP7ujMhOo9Ib+WiZApwjDsXRucpKjMG\nAqeTLiP/d3AckQ6pcBPMOBxY2J39orMUmRmLkB7k+rvzYnQekd7IBcmTwG7abeucGTsD3wY21Llf\nKSsztgaOIXWMqA2wE2ZcArygYxFSVBpOUnN5/P9o4LDoLEXnzquk++1+HZ1FpA92B55R0dYtF5I+\n6/aPDiLSG3m37efAb1S0dcuBwHZm/CQ6iEijzREdQBriEOAMd/4SHaQkjgAeN2MLd8ZFhxHpiTwK\n/GekuwmlC+58ZMZuwEQz3nLnguhMIj00EFgEuCg4Rym482pumbzFjNnc+U10JpFGUeFWcmb0A3YC\nVo7OUhbu/DNP47vIjG3duSc6k0gP7A887s6k6CBl4c4TZmwJ3GCGuXN+dCaRHjgUOExTJLvPnZfN\nGATcYcY97twYnUmkEXTGrcTMmB+YBJzozknRecomjw4+BRjuzpToPCJdyVd+PAEMdueR6DxlY8Zq\nwAR0ZYqURN45OgXopys/es6MocBxpH9/KnylEHTGrYbyFMlxpMLt5OA4peTOVcC+pJ23+aPziHTG\njNmAE4ErVLT1Tv73diRwVHQWka7kuwh/D/xaRVuv3QC8RupMEik97biVUH6AGwc8BHxPk9L6xoyz\ngBfd+Z/oLCIdycMJjgfWBLbQmOvey2cEnwW2cefB6DwiHcmv+QuA6cDu+pzvPTO2BQ4nTZPWv0cJ\npx23+tkN+DRwkN6EGuLnwIFmLBgdRGQWfg2sSyo2VLT1gTvvkVqnfhCdRaQTPwOWId1BqM/5vrkO\ncGCb6CAifaUdt5LJq3APAd935+boPFVhxrnAA+5qoZJiMWNt0oPH6u68Fp2nCsz4LGnXbS1N45Wi\nMWME8EdgfXf+Gp2nCvJ9jvu6Myg6i4h23OplK9LKkcbYN9ZxwBgzTVqVwjkM+KWKtsZx523gTOBb\nwVFEPiEPIDoJGKmiraEuAZY1Y93oICJ9ocKtfA4CjlHrRGO5czcwDdg+OovIf5jRH1gDOD06SwX9\nEdjTjHmjg4i08wNgrDt3RgepEnemA8cCB0dnEekLtUqWiBmrAzcCy+RzGtJAZnwVGOPOptFZRADM\nOBV41p3Do7NUkRmXAePdOSE6i4gZ8wHPkVokn43OUzX53+/zwLruPBccR2pMrZL18X3geBVtTXM5\nsLQZ60UHEclXfnwFODU6S4UdC3xPLdJSELsDbSramsOdf5LuxPtZdBaR3lLhVhJmLA6MAF203Sz5\nnpzfkaZMikTbC7jGnVejg1TYbaQW6d2jg0i95cFjBwB/is5ScUcAQ800YVLKSYVbeYwBznfnjegg\nFXcqsK4Za0UHkfoyY3bgQNLdbdIk+azwD4AjzFg6OI7U2wBgdqAtOEel5cFEI4Gz8hlikVJR4VYC\n+fD8fsDvo7NUXb4j60jgaLVPSaDtgL+5c1d0kKrLg4mOBC7Ol3OLRPgmcJIGjzWfO7cD3wYuyVeD\niJSGhpOUgBk/B/q5s3N0ljrIBdsNwH3u/DA6j9SPGW3Aie5cFJ2lDnKb2kXAdGBXPTxLK5mxFPAg\nafDY29F56sKMk4DZ3NkvOovUS19qIhVuBWfGasAEYAMdWG4dMxYG7gZ+rIdnaSUzvgxcDCyXR1hL\nC5jxaeA+4BB3LonOI/VhxpHAXO58LzpLnZixIPAU6fnqmeg8Uh8q3CrKjHlIxcOR7pwVnaduzFgT\nGA+s6c606DxSfXm3dzJwgl7zrWfGIODPwMrufBidR6ovH4V4HlhPI+pbz4wjSEXzQdFZpD5UuFWQ\nGbORxtO/4c5e0XnqyozDgUXd2Ts6i1SfGb8jXbi9jTsfReepm9wyeRtwrHbdpBVy4bCUO7tGZ6kj\nM5YFppD+G/w7Oo/Ug+5xq6bdgC8A+0cHqbmjgR1z66RI05ixI/BVYBcVbTHy2bbjSGPZRZrKjCHA\nnsDB0VnqKh9BuQ/YKTqLSHeocCsgM+YEfgV81533o/PUmTuvA1eBdj2lefJwgpOAkfnvnMS5EljF\njJWjg0h1mbEEcDawozt/i85TcycDo6NDiHRHl4WbmQ0xs6lm9pSZ/aiDf76wmd1gZg+Y2SNm9o2m\nJK2XHYAX8shaiXcCsH++W0ukGX4BnJpH00ugvFh2OmjSnDRHbsk9HThen/OFcA2wghkrRQcR6Uqn\nhZuZzU66AHYIsAowysz6zfBtY4D73b0/MBA42sx0/1Xf7I8u3i2MfJfWG8CW0VmkesxYjNSmc1R0\nFvnYKcDuZnwmOohU0jBgUeCI6CACeXrv2aizRkqgqx239YGn3f15d58OXAhsP8P3/BWYP//v+YHX\n3f2DxsasDzO+BKwJXB2dRT7hNNCAEmmKfYCL3XkjOogkebrfXcDI6CxSSfsAx+m6j0I5nbRYM2d0\nEJHOdFW4LQG82O7rl/KftXcKsKqZvUy6QPI7jYtXS7uRHuI03ahYLgC21JASaaTcMrUnaQS9FMtJ\naDiUNJgZXwA2BU0tLRJ3pgLPAMOjs4h0pqvCrTt3BfwUeMDdvwD0B/5kZvP1OVkN5Ye43UH3NxWN\nO2+R+uA1slkaaSPgPdJUMymW64Cl8n2OIo2yB3CJO/8XHURmcizwg/wsJlJIXZ1FmwYs1e7rpUi7\nbu1tBPwGwN2fMbPngJWAe2b8YWZ2aLsv29y9rYd5q24A8BHpThEpntOAP5pxXB4bLtJXuwPn6O9T\n8bjzgRmnkKbNfTM6j5RfLgj2BnaJziIdugI4lHQty0WxUaRKzGwgaQ5I339WZxdw5yEjTwCDgZdJ\nPf+j3P3xdt9zDPC2u//SzBYF7gXWcPc3ZvhZuoC7C2acBTzkztHRWWRm+VL0p4BReWCJSK+Z8WnS\n4tha7vwlOo/MLI9sfxhYzp03o/NIuZkxEPgjsIYWa4rJjA1IVwCt7s5r0Xmkmpp2AXceMjIGuBF4\nDLjI3R83s9Fm9p87Lw4H1jWzB4FxwA9nLNqka2YsQBr8cnZ0FulYvhT5ROAH0VmkEkYAD6hoKy53\npgHnoSm/0hj7kK79UNFWUO5MAc4BTlLLpBRRpztuDf1F2nHrlBkHApu589XoLDJrZswNPAns4c74\n6DxSXmbcApzsrpacIstXAjwK7OPOLdF5pJzMWBB4jrR7+3p0Hpm13A0xBTjKnXOi80j1NG3HTVrD\njE8BB5NaKKTA3PkX6XD5OXmXVKTHzFgRWJV0pkIKzJ13SdOSTzfj89F5pLR2Bm5Q0VZ8ear3fsDv\nzJgnOo9IeyrciuFnwKPuTIoOIl3LO21jgV9EZ5HS2h84y533o4NI19y5hnTu5bDoLFI+ueVuH9KA\nKymB3DJ5J2mAlEhhqFUyWD4IewWwtjt/i84j3WPGIsD9pEElt0bnkfIwYyVgMtDffaYpvVJQZnyO\n1Ca9rjvPB8eREjFjGHAkaeDFR9F5pHvMGEQ637qaziVKI6lVsqTylMJjgJ+raCsXd14ltbf+VgeY\npbvy35WTgMNUtJWLO28AF6IVeOkBM2YHfgv8WEVb6bQBnwbWCc4h8jEVbrG+Tbrk/MzgHNI7FwOf\nBzaMDiKlsRvwWTSlsKzOAL6RF91EumMP4E3g2ugg0jN5l+1stFgjBaJWySB5OuFzwCB3HovOI71j\nxk+BJd11Qa90zoyFSNMJh7lzT3Qe6bm8Y/oQcKBapKUreSLpk8BX3bkjOo/0nBnLks66LakzydIo\napUsp52BKSraSu984KtmzBUdRArvV8DFKtrKK6/An0naRRHpym+B21S0lZc7zwJPAEOis4iACrdI\no1CLZOnlIQWPozd16UReed+Z9CAn5XYRMEKLNdIZM9YDRoK6MSrgbLRYIwWhwi2AGYsDawPXR2eR\nhjgX2CU6hBTaMOBed16ODiJ9k4fKPA5sEZ1FCu03wCHuvBkdRPrsYmCLPFlWJJQKtxg7Adfki12l\n/C4DhuiiTunEzsB50SGkYS4GvhodQorJjBWBNUk7NVJy7rxNWmj/WnQWERVuMb5OGistFeDO34E7\ngG2js0jxmLEgsDlweXQWaZjLgO3M+FR0ECmkvYCz3HkvOog0jKZLSiGocGsxM74IrASMi84iDaUV\neJmVnYCb86qtVIA704AHSS2wIh/Lk0e/ClwQnUUa6iZgaTNWig4i9abCrfW+ClyhsbKVcyWwpRnz\nRQeRwtkFtUlW0enA3tEhpHDWIt3P+kB0EGkcdz4gTZHeLzqL1JsKt9ZTm2QFufMGMJE0RUwEADOW\nAlYHrovOIg13GbC2GatFB5FCGQlckq+OkGo5CtjNjDWig0h9qXBrITNWAJYE2oKjSHP8Gdg/t8qI\nQFqouVxnXarHnX8BRwNH6jUv8HGb5Ejgkugs0nju/BU4HPhFdBapLxVurbUbcJE7H0YHkaa4HpgX\nDSmR/1KbZLUdByxOmhoqsjZgwH3RQaRpTgEGmLFRdBCpJxVuLWLGnMA+pF0ZqaBckB8MHJX/e0uN\nmbEqsDBwa3QWaY58VvlHwE+06yakcfEXqk2yutx5BzgAOMOMOaLzSP2ocGudfYBH3Xk0Oog01fXA\nC8Do6CAS7gDgHHc+ig4iTTUOeBfYMzqIxDFjNnLhFp1FmsudK4C/kjoqRFrK3FuzMGRm7u61XJE0\nY37gSWCoO/dH55HmMmN10sPcyu68GZ1HWs+ML5HapVZ257XoPNJceVjBeGBZd/4ZnUdaL7fOnQKs\nph236jNjKHCoOxtEZ5Hy6UtNpB231vgRcIOKtnpw52HS9QD/E51FwvwCOFFFWz248xAwgXTxstTT\nN4ALVLTVxk3AkmasEh1E6kU7bk2Wx4E/AKzpzkvReaQ1zFgUmEracXklOo+0jhkrApOBFdx5KzqP\ntIYZmwHHudM/Oou0lhmLAY8BK2mxpj7MOBp4x11TJqVntONWbIcBJ6hoq5dcrF2GVuDr6FDgWBVt\ntXMbsGi+9kXq5TjgFBVttXMZsGN0CKkX7bg1kRn9SO0zK+jcQ/2YsQFwljsrR2eR1sjnG28Glnfn\n/6LzSGuZ8SfgJXeOiM4irWFGf+Ba0uf8u9F5pHXyQJoXgc3deSI6j5SHdtyKazRwqoq22robmDcX\n8FIPvwJ+p6Ktti4BdooOIS31LdJ5VhVtNZMnBl8BfCU6i9SHCrcmyfd77AqcFp1FYrR7U98hOos0\nX54suD5wYnQWCTMJWMKMZaODSPOZsSEwFDgpOouEUbuktJQKt+bZBHjeneeig0ioK9Cbel3sQ9ph\n18p7TbnzIXA1sH10FmmJI4EfufN6dBAJMwn4khlLRweRelDh1jwjSCPhpd5uBZbO93pJReUd9lHA\nmcFRJN5VpPd/qTAz1gOWBC6IziJx3PmA9JrXAq20hAq3JjDDUOEmfPymfi1aga+6AcAL2mEX0kXc\n/c34fHQQaaqRwLn5PV7q7TJ0zk1aRIVbc/QH3gcejQ4ihaBzbtW3HalFTmrOnX+TJosOi84izZEX\nZ3dAi7OS3AKsYsbi0UGk+lS4NccI4Cp3WnPXghTdTcA6ZiwUHUQaLz/EqXCT9q5A0yWrrB/wKeD+\n6CASz533SEX8HtFZpPpUuDXHDqQPbhHysIpxwPDoLNIU/YA5gQejg0hhXAlspBX4yhoBXKnFWWnn\nGOA7ZswTHUSqTYVbg5mxHLAIcGd0FikUrcBX1/bANXqIk/9w5x3SnW4HRGeRptAZdvkEdx4mtUwe\nEp1Fqs3cW/Os0ZdbwsvEjIOAldzZLzqLFIcZ8wJ/AdZw56XoPNI4ZtxNGgl+S3QWKY58l9vdQD93\nXo3OI41hxpKk3fXF3JkenUeKw4zFgKnAl9x5OzqPFFdfaiLtuDWe2iRlJu78H3Ae8J3oLNI4ZnwR\nWIZ07YPIx9x5Fjgf+Gl0Fmmo7YCxKtpkRu78jXQs4uvRWaS6tOPWQGasRHqA+2I+rCrysXze5WFg\nPY2NrwYzvg2s5c6e0VmkeMxYlLQCv4I7f4/OI31nxk3ASe5cHp1FiseMTYFzgFXzgq3ITLTjVhxj\ngNNVtElH3PkrcBpwYHQWaZgdQQ9w0jF3XiFdzrt3dBbpOzMWAL4M3BidRYrJnVuBycC3o7NINWnH\nrUHMWBmYBKyWP6xFZmLGMqRzL8u488/oPNJ7+QzTXcBSeXKoyEzMWIN0JcjK7rwVnUd6z4w9gB3d\n2T46ixSXGWsCY0mf82qplZlox60Yvgccp6JNOpNbJG8AvhudRfrsu8CpKtqkM+48RHqI0wp8+e1C\nOqssMkvuPAhMAzaLziLVox23Bmg3MXA1d16OziPFZsbypOsiVnTnjeg80nPtdthXzwfSRWbJjP6k\nlsll3fkwOo/0nBkrArejHXbpBjN+Cizuzreis0jxaMct3hjgJhVt0h3uPE26A0hn3crrYOD3Ktqk\nO9x5AHgd2CQ6i/TaL4FjVLRJN10FbG9GJTcsJI523PrIjM8CTwObuDM1Oo+UgxnrAheSps3p4uYS\nMePTwMuk3bZp0XmkHMw4FJjHnR9EZ5GeyWeWbiC9X2tSoHQpF2xPATvlhRuRj2nHLdZ3SXe6qGiT\nnrgX+ABYPzqI9Ni+wBQVbdJD1wLDokNIrxwEHK2iTborL8heBRpkI42lHbc+MGM24AVgqDuPROeR\ncjHjCOBDd/4nOot0jxnzA88AA915NDqPlEf+vJhG6s54OjqPdE8+w/4S6Uzyq9F5pDzMGAj8rzvr\nRWeRYtGOW5xNgDdUtEkvaQW+fPYFxqtok55y5yPSdMlto7NIj+wDTFTRJr0wGVjejMWig0h1qHDr\nmxHApdEhpLTuBJY0Y8noINJtewJ/jA4hpaXFmhLJZ9h/Bvw4OouUT77DbRwwJDqLVIcKt74ZClwX\nHULKKY8FvwGtwJeCGV8CFiEV3CK9MQ7YwIz5ooNIt+wG3OLO49FBpLSuA7aJDiHVocKtl8xYBlgQ\nuD86i5SaVuDLY1vget3DJb2Vh1vcAWwZnUW6ZS/g5OgQUmrXA1uaMWd0EKmGLgs3MxtiZlPN7Ckz\n+9Esvmegmd1vZo+YWVvDUxbTEODGfG5BpLduBDYz4zPRQaRLw0hnlET6Qos1JZBb2L8ITIzOIuWV\n7/qcihZrpEE6LdzMbHbgeFKRsgowysz6zfA9CwB/Aoa7+2rATk3KWjRDSSspIr3mzpvAfcCg6Cwy\na2bMDQwAborOIqU3Ftg2T5mU4toGuEE77NIA55DabkX6rKsPjvWBp939eXefTroweMY7KXYGLnP3\nlwDc/e/JlkSOAAAgAElEQVSNj1ksZnwKGIge4qQxxqIV+KLbHLjXnbeig0i5ufMs8BqwcXQW6dQw\n0u6oSF9dCGytQWTSCF0VbksAL7b7+qX8Z+2tAHzOzCaY2T1mVodVhU2Ax9x5PTqIVMK1wDAzKnXP\nYcWoTVIa6UzSmHkpoNy6PpDUyi7SJ+68AZwBHBydRcpvji7+eXdu554TWBsYDMwN3GFmd7r7UzN+\no5kd2u7LNndv62bOolGbpDTSVGA60B8NuymcXFBvC2wVnUUq4yxgqhlLuX9icVSKYSDwQG5lF2mE\no4FHzDjCnVeiw0hrmdlA0vtKn3VVuE0Dlmr39VKkXbf2XgT+7u7vAu+a2a3AmsBMhZu7H9r7qMWQ\nH+JGUJ+zfNJk7rgZZwP7A6Oj88hMVgfeJxXYIn3mzmtmnAgcgnbeikhX/UhDufOyGZeSXu+/ic4j\nrZU3qtr+87WZHdLbn9VVq+Q9wApmtrSZzQV8Dbh6hu+5ChhgZrOb2dzABsBjvQ1UAmuRdiIfiA4i\nlXISMNKML0QHkZlsC4x171YHgkh3/QHYyYx5ooPITNRVI81wBrCbjkVIX3RauLn7B8AYUp/3Y8BF\n7v64mY02s9H5e6aSLhF+CJgCnOLuVS7c9gIu0EOcNFJunTgF+GV0FpnJ9mhIgTSYO68Ck4GvRGeR\n/zJjeWAe0jONSCPdSVr419UA0mvm3pr6w8zc3Uu9ymDGQsDTwKruvBydR6rFjEWBJ4DF3Pl3dB4B\nM5YhLUgt4c706DxSLWZsAxwB9NdiYDGY8S1gLXf2is4i1WPGV4CfAuvqNV9ffamJdI9Mz/wEuFhF\nmzRD3nV7iDToR4phN+AyFW3SJP9pxxsamkLa2wa1SUrzXAEsAKwbHUTKSTtu3ZQv4P0r0E+FmzSL\nGaOB7d3ZJjpL3eWzR88Bm7nzeHQeqSYzRgEHuLNpdJa6M2MR4ElgKXf+GZ1HqsmMnwOLujMmOovE\n0I5ba2wJ3K2iTZrsDKCfmS7nLYDRwEQVbdJklwArmLFidBBhV+BqFW3SZJcD22pIifSGCrfu+xpp\ngqZI07jzPvC/6KLOUGbMARwEHB6dRarNnQ+Ai4Cdo7PUWb50+yDguOgsUnmPka7jWiE6iJSPCrdu\nMGNVYAvg7OgsUgtnApuZsXh0kBobAPzNXReiS0tcQppeKnF2A+5x597oIFJteSjJDcCQ6CxSPirc\nuueHwDHuvB0dRKrPnX+RLn/dITpLje1IamcRaYUpwJe0WBNqT+DP0SGkNlS4Sa+ocOuCGZ8HhqM3\ndGmty0jFg7RYfs3vDFwQnUXqIbdLjgO2is5SR2b0A5Ym3Vkr0grjgQG5RVek21S4dW0bYLw7b0QH\nkVoZB2yQp5lKax0EXOTOs9FBpFZuQNcCRNkDOCcX0CJN585bwIPAJtFZpFxUuHVtCOkDVaRl8lSz\nB9CbekuZ8SlgLzSgQFrvRmBLM2aPDlInZsxHapM8PTqL1I7aJaXHVLh1In+AbonaJyTGzaS/f9I6\nWwOPu/NkdBCpF3emAdOA9aKz1MwYYJw7U6ODSO1ol116TIVb59YhTZZ7KTqI1JIKt9YbAlwbHUJq\nSyvwLZTv0dob+H10Fqml+4GFzFg6OoiUhwq3zm2Ndtskzt3AF81YNDpIHeSHuKGoNVriqHBrrQ2A\n6aArAKT13PmI9Iy5dXQWKQ8Vbp3bGj3ESZB8UL6NdIegNN8KwFzAI9FBpLYmA/3MWCg6SE1sA1yZ\n79USiTAWXf0jPaDCbRbMWABYA5gUnUVqTe2SrTMEuEEPcRLFnfeAieg13ypbkN5jRaJcBaxnxlLR\nQaQcVLjN2mBgsjv/jg4itXYzadKcRQepgSGoNVriXQ9sGx2i6sz4LLAacHt0Fqkvd94FzgO+E51F\nykGF26zpGgApgqeBd4B1o4NUWb4GYADp/jyRSJcA25qxWHSQihsI3KHFWSmA3wJ7mrFIdBApPhVu\nHci7GxpMIuFy295FwNeis1TcRqRrAN6IDiL15s7fgfOBb0dnqbgt0UKNFIA7L5OeN0dEZ5HiU+HW\nsVWBj4AnooOIkNoodjXj09FBKmwr4KboECLZ0cB+ZswfHaTCtkCFmxTH5cBXokNI8alw69jXgEs1\npECKIF8Mex+wS3SWCtsSFW5SEO48B9xBmnooDZYHQSwEPBidRSS7HljfjMWjg0ixqXCbgRmzATsD\nF0RnEWnnaOCg/PdTGsiMz5OuArgzOotIO9ehwq1ZtgDG53u0RMK58w5wKfCN4ChScHoInNmOwN9J\nOxwiRXEL8D4wKDpIBQ0nPcRNjw4i0s51wBAz5ooOUkFD0DUAUjwnAGN0LEI6o8JtZgcDh6tNUook\n/328hFRkSGPtDpwTHUKkPXdeILXyjYrOUiW5EN4auDY6i0h77twPPADsGp1FikuFWztmrAksgd7Q\npZiuI40J151uDWLGGsBKpH+3IkVzFHCwXvMNNQiY6s4r0UFEOnAq6biOSIdUuH3SzsDZ7nwYHUSk\nAw8AjtolG+lnwFHuvBcdRKQDN5Fe81tFB6mQvdEOuxTX9cBausdRZkWF2ycNA66KDiHSkdwu+Svg\nF9FZqiCPWh8KnBadRaQj+TV/ErBbdJYqMGNJ0gTZc6OziHQkXwg/Hi3WyCyocMvMWBFYGLgnOotI\nJy4CVjJjleggFTAcmOjOW9FBRDpxJbCNhpQ0xI+BU9x5OzqISCduJi0wiMxEhdt//QQ4SeOBpcjy\n5MMz0MjgRhgJXBwdQqQz7rwMPAUMiM5SZmbMSboL89joLCJduAnYUmdbpSMq3IDcSzwC+H10FpFu\nuBxNl+yT3CY5CLg6OotIN9wMbB4douQ2Ap7JhbBIYbnzHPB/wOrRWaR4VLglewCXq2VKSuI+YAEz\nlosOUmJ7AOPUMiUlcQsq3PpqOGnwg0gZqF1SOqTCLRmJDitLSeR23pvQm3qv5LNCPwYOi84i0k13\nAGvknWLpoXZtkudFZxHpJn3GS4dqX7iZsRCwIjA5OotID4wHBkeHKKkhwLP5slORwnPnXeAuYJPo\nLCW1PfC0O1Ojg4h00wRgIzM+HR1EiqX2hRup/WSSO+9HBxHpgVuAQWZ6DffCrugeJykftUv2Qh7w\n8CPg6OgsIt2Vj+48CmwcnUWKRQ99sAUwLjqESE+48xLwOrBGdJYyMWN2UvuJ7muUslHh1jtrA59D\ng4ikfNQuKTNR4ZZeFDdHhxDpBbVL9txawDR3XokOItJDdwPL5vZ+6b6RwMW66kdKSANKZCa1LtzM\nWBaYm7QdLVI2WoHvuc1JZwdESiXf4Xgb6RoL6b4RwKXRIUR6YQqwnBmLRgeR4qh14UZuk3THo4OI\n9MIEYECemCbdszmp4BUpoxuBYdEhysKMJYCFQYOIpHzyYs1lwD7RWaQ4VLipTVJKyp3XgWeB9aKz\nlEG+BmAjYGJ0FpFeugzYzoxPRQcpiUHARLVJSokdB3wzf36J1Ldwy9P4NiedExIpq/GoB7671gOe\ncueN6CAiveHONOABYMfoLCUxCLVGS4m58xAwlXRWU6S+hRvQH3gtT+cTKatLgVF55LV0TufbpAqO\nAX6o13y3qHCTKjgB+EZ0CCmGOhduugZAqmAK6XW8YXSQEtD5NqmCscB8qEW6U2YsDcwDPBYcRaSv\nxgFf1mXcAvUu3LZF59uk5PJgnT8CP4zOUmRmzE160J0UnUWkL/Jr/ky0At+VQUCbho9J2bnzNvAI\nuoxbqGnhlq8B6AfcEJ1FpAFOJU2XXCo6SIENBu5255/RQUQa4BJguNolO6U2SamSq4FR0SEkXi0L\nN2B/4Fx33o8OItJX7rwLXIfGhHdmGHBtdAiRBnkScGCl6CBFlAtaFW5SJacBXzFjoeggEqt2hZsZ\nCwJ7A8dGZxFpoKuB7aJDFFFuk/wKcHl0FpFGyO1/44Ah0VkKanXgI1KBK1J67rxK+pzfOzqLxKpd\n4QZsT+p7/0t0EJEGuhHY2Iz5ooMU0C7A7e48Fx1EpIHOAfZTu2SH9gTO1vk2qZg/AgfoNV9vdSzc\nhpNWLUQqI5/duh3YKjpLkeQPuG8Df4jOItJgbcCHpAnJkuWLinchDXARqZJ7genA2tFBJE6tCrc8\nSnUw6TyQSNVcSlpplv/aGJiDdFG5SGXk3aTjSAsT8l/DgMfceSY6iEgj5df8lcAO0VkkTpeFm5kN\nMbOpZvaUmf2ok+9bz8w+MLMdGxuxoQYCj7jzWnQQkSY4F1jbjDWigxTIUOBytUxJRZ0HbGLGItFB\nCmQv4IzoECJNMhZ11tRap4Wbmc0OHE86AL0KMMrM+s3i+44kjdcvcu/tcOCa6BAizeDOv4FTSA8u\nkmyB7muUisoTZW8Bto7OUgRmLE7aZb80OotIk0wBVjFj/uggEqOrHbf1gafd/Xl3nw5cSBruMaNv\nkd4oC7uTZcZnga8BF0VnEWmic4CdzZgjOki0vAuxMnBHdBaRJroO2CY6REGMAK5x553oICLNkBdo\n7wIGRGeRGF0VbksAL7b7+qX8Zx8zsyVIxdyJ+Y+K2pK0L3CDO89HBxFpFneeBv4GrBOdpQB2Asa6\n8150EJEmuhnYXJPmANicdE2CSJW1kY7+SA11tSrfnSLsWODH7u5mZnTSKmlmh7b7ss3d27rx8xtl\nN3SIW+phAuny2SnRQYLtDvwmOoRIM7nzghnvAP2Ax6LzRDFjNtLD7PeDo4g0WxtwVHQI6T4zG0iD\niu2uCrdpwFLtvl6KtOvW3jrAhalmY2FgqJlNd/eZRu67+6G9j9p7ZvQDPgdMivj9Ii02ATgA+G10\nkChmDCC9H2mCrNRBG2mxpraFG7AJ8Ir7J7qERKroLvI5N3f+ER1GupY3qtr+87WZHdLbn9VVq+Q9\nwApmtrSZzUU6I/aJgszdl3X3Zdx9GdI5twM6KtqCbUFqk/woOohIC0wENsz3GdXV/sAf3PkwOohI\nC0xArVMHAidFhxBpNp1zq7dOCzd3/wAYA9xIWsm7yN0fN7PRZja6FQEbZCDtKl2RKnPnTeBpYL3o\nLBHyfY3bAhdHZxFpkTZgYG4XrB0zliAt0J4dnUWkRdrQYk0tmXtrZomYmbt7yw9P5w+yV4H+7jO1\neYpUkhlHA2+6c1h0llYzYyjwE3c2jc4i0ipmPA2McOeR6CytZsYhwCLuHBidRaQVzNgUONq9ngu0\nZdeXmqgOq3OrAG+paJOaaYPaFi6D0GQ5qZ//DCWqo52Ac6NDiLTQFGA5M74QHURaqw6F20DUJin1\ncxvwZTPmjA4SYCB6zUv9tFHD1ikzlgUWIZ35EamFfM3NNcDI6CzSWnUo3DZDD3FSM/mc27PU7D43\nM+Yj7bLrIU7qZgLpnFtX06KrZghwvQYRSQ2dB+ylOxzrpdKFW/7LvBlpyp5I3dxMeqipk42Be/LU\nLZHacOdl4Clgy+gsLabFWamrm0nP8VtEB5HWqXThRlp5/6fudZGaGgsMiw7RYgPRQ5zU15nA3tEh\nWiUvzm6KFmelhtxx0hUYe0RnkdapeuE2EL2hS31NBpbOZ0DqQjvsUmfnA4PMWDo6SIusAHwAPB+c\nQyTKZcAwMz4THURao+qF2zDg+ugQIhHcmQ6cBZTpzsVeM2NeYHXgzugsIhHc+Qdp123f4CitMhCY\nkHceRGrHnb8BD1PfKdK1U9nCLQ8pGEC6PFykrk4C9syXUlfdRsB97rwbHUQk0JXU52yrzreJpMFE\nm0WHkNaobOEG7AjcmlcgRWrJnaeA+6nHyOBBqE1S5E7S/U6LRgdppny+bSB6zYu0UcOrQOqqyoXb\nAcDJ0SFECuBs0kJGZeWHuB2Aq6OziETKLdK3AFtFZ2my5YGPSNeeiNTZHcBKZnwxOog0XyULNzPW\nAr5AmqonUncTgM3Mqvl6z1YB5gbuiQ4iUgA3AFtHh2iygUCbzrdJ3eXjAecC+0Vnkear6oPcaOAU\nXcgp8vH9Tq8Ba0ZnaaL9gXP1ECcCpLPdW1f8Mu6tSTuLIgInAvuYMVd0EGmuyhVuuWVqW+Ci6Cwi\nBXITsE10iGbIg4h2Bf4UnUWkCNx5AXgc2Ck6SzOY8VnSReNXRmcRKQJ3pgKPko4MSIVVrnADlgHm\nAJ6KDiJSIFdQ3Tf0HYFJ7kyLDiJSIEcD34kO0SRfB8a782Z0EJECuQAYER1CmquKhdtmwES1TIl8\nwq3AUmasHB2kCUYB50WHECmYscDSZqwUHaSR8lnd7wF/iM4iUjDXA1tVvEW69qpYuA0nHcwWkcyd\nD0g98AdFZ2mk3M+/MXrNi3xCfs2fD+wcnaXBvgx8gK4BEPmE3HXyIrBedBZpnkoVbvmsy2Dgqugs\nIgV0IjDSjM9EB2mgdYCn3Xk7OohIAY2letMlhwNXqatGpENt6DLuSqtU4UZqmWpT37vIzNx5hTQu\nv0pDSjYjtYGKyMxuB1Y1Y8HoIA20LXBtdAiRgmpDl3FXWmUKtzxN8nvA76OziBTYZcD20SEaaCf0\nECfSIXf+Tbqcd9PoLI1gxiLAF4G7o7OIFNQkYEOdc6uuyhRuwIrAPKjvXaQz44DBeaGj1MzoB3wB\n3eUk0plJwIDoEA2yKXBbPr8nIjNw53XgBWDt6CzSHFUq3AaTxgOr711k1p4GPoRKTJr7JnCGOx9G\nBxEpsNuoTuG2JTAhOoRIwbWhdsnKqlLhtjlaeRfpVF7YGE96vZSWGfMAuwAnRGcRKbgpwBpmzB0d\npC/ya34n4OLoLCIFNxENKKmsShRu+V6XQahwE+mOW0g71GW2NXCPLt0W6Zw7/wIepvwjwr9OapN8\nMTqISMHdCmysc27VVInCDVgTeE0PcSLdMh4YmBc8ymo74OroECIlUYV2yQOAk6JDiBSdO68BLwH9\no7NI45X5wa09tUmKdJM7LwOvUtI39TxYZSvguugsIiUxiRK3TpmxCrAIcGN0FpGSmED17nAUqlO4\nDSbtIohI95S5XXIZwIHnooOIlMQE0ojweaOD9NJQYKw7H0UHESmJs4G9St5ZIx0o/X9QM+YktYC0\nBUcRKZMyDygZQDrrogmyIt3gzj9IQ0q2iM7SS0OAG6JDiJTIPcDblPc1L7NQ+sINWB94Jt9dISLd\n00Y6vDxXdJBeGEA6syMi3Xc5acBHqeT3qA3R4qxIt+WFzZOB0dFZpLGqULhtjtokRXrEnTeAp0gL\nH2Wjwk2k584Hhpjx+eggPbQW8LQ7b0cHESmZ84GtStwiLR2oQuE2GA0mEemN8aQLbUvDjIWBJUjj\nzUWkm9x5C7gK2D06Sw9poUakF9z5J3AvsEl0FmmcUhduZswPrEOamCUiPXMBsGc+J1oWA4Ap7nwQ\nHUSkhP4M7Jsns5aFCjeR3ivzeXbpQKkLN2B7YEJeVRCRHnDnfuB5YHhwlJ7QBFmR3rsdmAdYPjpI\nd+QCc2NUuIn0VpknSEsHyl647QJcGB1CpMTOA0ZGh+iBwcC46BAiZZQHFtwKbBqdpZtWAN5156Xo\nICIldRewnBkLRQeRxiht4WbGSsDawBXRWURK7EpgqBmfjg7SFTOWBRYGHojOIlJikyjPmZdN0W6b\nSK+5M530GhoUnUUao7SFG7APcJo770YHESkrd14hTZdcNzpLN+wGXOjOh9FBREpsArBFSc65DQFu\njA4hUnJXU8KrQKRjZS7ctgMujQ4hUgGTSQMACsuM2YFvAGcHRxEpuyeB94DVo4N0Jg9N2gIVbiJ9\ndQEw2IzFooNI35WycDNjRWBe4L7oLCIVcBsFL9yAbYFX3bknOohImeVzbteRFj+LbCjwaO4KEJFe\ncucfpI2OPaOzSN+VsnAj9eqOyx9AItI3E4EBZnwmOkgnRgGnRocQqYgzgf3MmCM6SCf2A06JDiFS\nEScD+5SkRVo6UdbCbWN0d5tIQ7jzGumSzqHRWTqSP2gGATdHZxGpAnfuBV4AhkVn6YgZSwIbAZdE\nZxGpiHuB2YGVo4NI35S1cNOFnCKNdSGwe3SIWegH/Mud56ODiFTIOaSd7CLaE7jAnXeig4hUQe5Q\nuxnYMjqL9E3pCjczVgXmAp6IziJSIeeT2iWXiQ7Sgc1Jl4iKSONcBgwxY+7oIB0YDlwcHUKkYlS4\nVUDpCjfSCuGFOt8m0jh5ZfsCYOfoLB3YnDTCXEQaxJ3XgUdILYmFYcZnSbvsd0ZnEamYScDGZqV8\n9pesVP/x8kHqPUgtHiLSWJcDO0SHaC9fAzAQFW4izTAeGBwdYgYDgLvceS86iEiVuPNX4C10zq3U\nSlW4AdsDz7vzYHQQkQqaBCxjxheig7QzGHjOnZejg4hU0C0Ur3AbBLRFhxCpqNtIA/6kpMpWuI0C\nTo8OIVJF7nxAKt42jc7Szj7oGgCRZrkD6GfGAtFB2hmEdthFmmUyxb+3VTpRmsLNjE+RDlWOjc4i\nUmG3UpDCzYzPk17z50dnEami3I54B7BZdBYAMz4HrAjcFZ1FpKImox23UutW4WZmQ8xsqpk9ZWY/\n6uCf72JmD5rZQ2Y22czWaHxUBgCPufNqE362iCSFKdyA3YCr3Hk7OohIhRWpXXIUcJ0770cHEamo\nx4CFzFgsOoj0TpeFm5nNDhwPDAFWAUaZWb8Zvu1ZYFN3XwP4NfDnRgcFNgEmNuHnish/PQB80YyF\no4OQBqVcGB1CpOLGkya3hjLDgH2BU6KziFSVOx+RzpDqWoCS6s6O2/rA0+7+vLtPJz1Ibd/+G9z9\nDnf/z6r4FGDJxsYE0tbu5Cb8XBHJ8jm32wnugTdjQWBNtFgj0mz3AUsUYAV+XWB+dGejSLNdzQzP\n8VIe3SnclgBebPf1S/nPZmVv4Lq+hJpRvgZgA9IDpYg0VxHaJQcDt7nzbnAOkUpz50PSAsmg4Ch7\nAqflHQERaZ5rgS3MmD86iPRcdwq3bl90bWaDgL2Amc7B9dHqwLR8YaiINFcRCrcBOYeINF8R7nMb\nTHqgFJEmcuc1YBzpHLmUzBzd+J5pwFLtvl6KtOv2CXkgySnAEHd/s6MfZGaHtvuyzd3buplzI9Qm\nKdIqdwMrmzG/O/8IyrARcFDQ7xapm/HA96N+eW7TXAR4JCqDSM2cCBwN/Ck6SB2Y2UBgYEN+lnvn\nG2pmNgfwBGk17GXSmN5R7v54u+/5IqkvfVd3v3MWP8fd3XoV0jgfuNmdM3rzfy8iPWPGBOAo99Zf\nv2HG3MBrwMJqlRRpvjwY5Bng6+6tH8Vvxkhgd3eGt/p3i9SRGbOTjkFt7s7U6Dx105eaqMtWSXf/\nABgD3EgaI3qRuz9uZqPNbHT+tl8ACwInmtn9ZtawN34zZiNNvNKQApHWuQzYJeh3rw88rKJNpDXc\nceAE4FtBETZFrdEiLZPPtl4C7BidRXqmyx23hv2iXlaXZvQHLnFnhSbEEpEO5KmOzwHLtfpsqRk/\nAz7nrlZJkVYxYxHgSWCRVt+jZsaDwH7uTGnl7xWpMzOGA99219UArdbUHbcCGEra7RORFnHnTWAs\nMbtuOtMq0mLuvAo8TosHE+VFomVJ1xKISOtMAr5sxlzRQaT7Cl245b77PYALorOI1NBppBHdLWPG\nvKQ7G9UaLdJ6Y0mLpa20MTDFnekt/r0itebOW6Rd9vWis0j3FbpwI72hO7q/TSTCROBLZizewt+5\nHTBZV3+IhJhE2vFuJZ1vE4nTRoOmHUprFL1w2wc4NR+cFpEWyoeXb4GW9r+PQjvsIlHuBtYw4zMt\n/J0q3ETiTESFW6kUtnDLHxw7AOdEZxGpsZuArVrxi8z4HOkh7qpW/D4R+ST///buPViPujzg+Pfh\nEi6BgqggTSKgXASVIdAGSsgNCMMdqnJTnAwwWrVepp3pVOxMtX+1/aPTy7SiIySgZUBADSBEgyEB\nEjDcEkDuUKOAcpkpQRCxIE//2A2cnCTnuvvuvnu+n5nMed99d3/7zDzZc95n93dJXqWYPbonXaci\n2AP4ALDZZYQk1e424PAIdm46EI1Maws34HDg4XLAtKRm3ATML5flqNtHgaWZvNyDc0navKX0bpzb\n6cCSTF7r0fkkDVBORLYcOLPpWDQybS7c7D4hNSyTnwMvAR/uwensJik17zqKsaa98Engqh6dS9Lm\nXQxc0HQQGpk2F25zcGY5qQ1q7y4ZwbuBQ4Eb6zyPpGHdBewcwfQ6T1Ku0boXcH2d55E0rCXA3hEc\n1HQgGl4rC7dyTYkZuJaT1Aa9GOc2m2I2SbtMSQ3K5E3gIuALNZ/qXGBRJm/UfB5JQyivwcuBs5uO\nRcNrZeEGHAY8Ua4xIalZKygW6Zxc4znsGi21x0LgIxFsX0fj5RqtpwPfr6N9SaO2DJjVdBAaXlsL\nt3n4JU5qhUx+Q9Ft+awaT2PhJrVEJs8Ba6nvSfsBwCTgvpralzQ6dwB/WvZ4U4u1rnAr78SdC1zT\ndCyS3vJfwGfraDiCXYF9gXvqaF/SmCwGTq6p7ZnALa7RKrVDJi8Bj1P0eFOLta5wo1gGYFtgZdOB\nSHrLTcD7IphaQ9szgdWZ/F8NbUsamxUUk4TV4Ujg9praljQ2K4Gjmg5CQ2tj4XY+sNA7cVJ7lIOX\nfwScWEPz8ygWAZXUHg8Au0fwnhranknRNUtSe9yG49xar1WFWwTbAWcAlzUdi6RNLAHmV9lgBNtQ\nrN9m12ipRTL5A8Ud+Eq/yEUwBXg3RWEoqT1WAjMj2lUbaGNtS84Mitkkf9V0IJI2seGXelTY5vHA\nU5k8WGGbkqpxC9V3l5wP/KQsDCW1RPnd+3+BDzUdi7asbYXbHJxZTmqrX5Q/96qwzQuAiytsT1J1\nbqWY8bVKJ1GsDSmpfZZTDF9QS7WxcLul6SAkbaocd7qKYnzKuEXwToo/EN+toj1JlVsD7B3BblU0\nVl7z84HvVdGepMqtAI5pOghtWWsKtwi2BY7ASQqkNrudigo3irEzP83k5Yrak1ShTF6nmESkqnFu\nHwduyGR9Re1JqtYSYE4Ef9R0INq81hRuFGtHPJnJi00HImmLVlFM5V0Fn7BL7VdJd8lybOwFwMJx\nR6hBwoUAAArcSURBVCSpFuV38NuAU5uORZvXpsLNL3FS+60F3h/BLhW0NRvHtEptt4JqZpM9ENiN\nYgyNpPb6LnBW00Fo89pUuPklTmq5cpHslcBx42mnLPwOAO6uIi5JtbkD2CWC6eNsZzawLJM3K4hJ\nUn2upegu+Y6mA9GmWlG4RbA1xbgZCzep/a5j/N0ojgTuzOT3FcQjqSZlobUIOH+cTR1FcdNHUotl\n8htgGXB607FoU60o3IBDgGcyeaHpQCQN61rgpAgmj6MNl/6Q+sci4JwIth9HGxZuUv+wu2RLtaVw\nOwX4cdNBSBpeuUjnrcAnxtHMbBzTKvWFTH4B3AecMJbjI5gGTAYeqzIuSbW5HpgVwU5NB6KNtaVw\n+yiu6yL1k28BC8ZyYAQ7AgcDqyuNSFKdlgDHjvHYmcCqci1ISS2XyW+BB4A/aToWbazxwi2CQ4Bd\nKQZAS+oPNwEfiGDqGI79M+C+TF6tOCZJ9bkZOHqMx9pNUuo/d1D8vVaLNF64AZ8DvulMU1L/KGeX\nvIGim/NonU5x915S/1gL7B7BlDEcOw+7Rkv9xsKthRot3CLYCvhz4NtNxiFpTFYAs0ZzQDmD7BnA\nlXUEJKke5c3V5YzyqVsEewJ7AvfWEZek2iwFvtp0ENpY00/cPgSsz+SXDcchafRWUnSBGo0TgXWZ\nPFFDPJLqNZbukvOBFZn8oYZ4JNUkk/WZrGk6Dm2s6cLtaIo/BJL6z+PAdhG8dxTHfBa4qKZ4JNVr\nGXBMBDGKY84Frq4pHkmaUCzcJI1JOUPcSkbYXbKcTXI2ziAr9avHKL43vH8kO0ewF3AYsLjOoCRp\nomiscItgG4ovccubikHSuN3GyLtLzgTWZvJKjfFIqkl5s+Zm4JgRHrIAuCKT39UXlSRNHE0+cTsU\n+GUmzzcYg6TxWQnMHeG+8/BGjdTvrgfOGm6nsjvlAmBR7RFJ0gTRZOFmN0mp/90L7BDB9BHsa+Em\n9b9rKdZw/OAw++0HbIezSUpSZSzcJI1ZOUX4pcB5Q+0Xwc7AhynWhZHUp8o1HP+b4Z+6HQMsK7tX\nSpIq0EjhFsEOFIv63drE+SVV6lLg4xFsN8Q+s4C7HOsidcK1wGnD7DOPYhZKSVJFmnri9hFgVSbr\nGzq/pIpksg64Dzh5iN3sJil1x0+BKRFMHWKfI4FVPYpHkiaEpgq384GFDZ1bUvW+D5w0xOd2jZY6\nolxMexXFTLGbiGAaMAn4n17GJUld1/PCLYL3AQdTdLWQ1A0/Ao7f3MK8EbwT2B+4s+dRSarLFgs3\niqEQdzi+TZKq1cQTt3Mp1nX5fQPnllSDTJ4EfgebnWnuHOCH5aQGkrphqMLtWGBF70KRpImhicJt\nPvDDBs4rqV63svnFuC8ALulxLJLqdQ/FsgA7DdxYPnU/AVjSSFSS1GE9Ldwi2JFi4e3be3leST2x\nyR34cn23d+D4NqlTMnkNWAvMGPTRdOAN4NGeByVJHdfrJ25HAPdn8kqPzyupfivZ9InbAuDScr03\nSd2yik2v+c8ACx3fJknV63XhNhf7vUtd9SiwcwRTBmyza7TUXUsZMJtsBJOBM4GLG4tIkjqs14Xb\nHCzcpE4q77C/1V0ygt2BKRTdqSR1zy3AvgPWczsFWJ3JrxuMSZI6q9eF22E4vk3qspUUN2gAZgOr\nMnmjwXgk1SST14HFwCfKTWcAVzYXkSR127CFW0QcHxGPRMTjEfG3W9jnP8rP74uI6UM097NMXh5z\ntJLabjHwsQi2xa7R0kTwDeAz5TV/NMWajpKkGgxZuEXE1sB/AscDBwHnRMSBg/Y5Edg3M/cDPg1c\nNESTt4wvXLVVRMxtOgbVYzS5zeRxirFuJ2HX6L7gtdtdvchtJncB64G/Bp6xm2TveO12l7nVlgz3\nxG0G8ERmrsvM1ym6QJw2aJ9TgcsAMnM1sGtE7LGF9izcumtu0wGoNnNHuf8i4CvAe4E1lUejqs1t\nOgDVZm6PznM58A/AtT06nwpzmw5AtZnbdABqp+EKtynAUwPeP11uG26fqWyeazlJ3Xc1cCDwVce3\nSRPCt4F/Ab7WcByS1GnbDPP5SNdhiZEcVy7YKanDMnklggPALlPSRJDJ88DfNR2HJHVdZG65NouI\nI4CvZebx5fsLgTcz858H7PMNYEVmXlm+fwSYk5nPDWrLxTglSZIkTWiZOfih14gM98TtbmC/iNgb\n+BVwFnDOoH2uAz4PXFkWeusHF23jCVCSJEmSJrohC7fMfCMiPg/8GNgauCQzH46Ivyg//2Zm3hgR\nJ0bEE8BvgfNqj1qSJEmSJpAhu0pKkiRJkpo37ALc4zWSBbzVbhGxMCKei4gHBmzbLSJuiojHImJp\nROw64LMLy3w/EhHHNRO1RiIipkXE8oh4MCJ+FhFfLLeb3w6IiO0jYnVErI2IhyLiH8vt5rcjImLr\niFgTEdeX781tR0TEuoi4v8zvneU289sREbFrRFwTEQ+Xv58PN7/9LyIOKK/ZDf9eiogvVpXbWgu3\nkSzgrb6wiCKHA30ZuCkz9weWle+JiIMoxkIeVB7z9Yio/QaBxux14K8y84PAEcBflteo+e2AzHwN\nmJeZhwAHA/Mi4ijMb5d8CXiIt2dzNrfdkcDczJyemTPKbea3O/4duDEzD6T4/fwI5rfvZeaj5TU7\nHTgMeBX4ARXltu6kj2QBb7VcZt4GvDho81sLr5c/Ty9fnwZckZmvZ+Y64AmK/wdqocx8NjPXlq9f\nAR6mWJvR/HZEZr5avpxEMVb5RcxvJ0TEVOBE4GLeXpbH3HbL4IndzG8HRMQuwKzMXAjFnBKZ+RLm\nt2uOpaiDnqKi3NZduI1kAW/1pz0GzB76HLBH+fqPKfK8gTnvE+XssdOB1ZjfzoiIrSJiLUUel2fm\ng5jfrvhX4G+ANwdsM7fdkcBPIuLuiPhUuc38dsM+wAsRsSgi7o2Ib0XEZMxv15wNXFG+riS3dRdu\nznwyAWQxw81Qufb/QctFxE7A94AvZebLAz8zv/0tM98su0pOBWZHxLxBn5vfPhQRJwPPZ+YaNn0q\nA5jbDphZdrc6gaIb+6yBH5rfvrYNcCjw9cw8lGJW9i8P3MH89reImAScAlw9+LPx5Lbuwu0ZYNqA\n99PYuKpU/3ouIt4DEBF7As+X2wfnfGq5TS0VEdtSFG3fyczF5Wbz2zFlN5wbKPrcm9/+dyRwakT8\nnOKO7tER8R3MbWdk5q/Lny9QjJGZgfntiqeBpzPzrvL9NRSF3LPmtzNOAO4pr1+o6Nqtu3B7awHv\nsvI8i2LBbvW/64AF5esFwOIB28+OiEkRsQ+wH3BnA/FpBCIigEuAhzLz3wZ8ZH47ICLetWHmqojY\nAZgPrMH89r3M/EpmTsvMfSi649ycmZ/E3HZCROwYETuXrycDxwEPYH47ITOfBZ6KiP3LTccCDwLX\nY3674hze7iYJFV27Qy7APV5bWsC7znOqehFxBTAHeFdEPAX8PfBPwFURcQGwDjgTIDMfioirKGY5\newP4XLpYYJvNBM4F7o+INeW2CzG/XbEncFk5Q9VWFE9Vl5W5Nr/dsiFPXrvdsAfwg+LeGtsAl2fm\n0oi4G/PbFV8ALi8fbDwJnEfxXdn89rnyZsuxwKcGbK7kd7MLcEuSJElSy7kGhCRJkiS1nIWbJEmS\nJLWchZskSZIktZyFmyRJkiS1nIWbJEmSJLWchZskSZIktZyFmyRJkiS1nIWbJEmSJLXc/wPtHhHH\nT+GJogAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x106c2da90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "fig = plt.figure(figsize=(15,5))\n",
    "ax = plt.subplot(111)\n",
    "plt.plot(times, es);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that the eccentricity is oscillating between 0 and almost 1. \n",
    "\n",
    "Note that the function `starkForce(reb_sim)` above receives the argument `reb_sim` when it is called. This is a pointer to the simulation structure. Instead of using the global `ps` variable to access particle data, one could also use `reb_sim.contents.particles`. This could be useful when one is running multiple simulations in parallel or when the particles get added and removed (in those cases `particles` might change). The `contents` has the same meaning as a `->` in C, i.e. follow the memory address. To find out more about pointers, check out the [ctypes documentation](https://docs.python.org/3/library/ctypes.html#pointers).\n",
    "\n",
    "\n",
    "**Non-conservative forces**\n",
    "\n",
    "The previous example assumed a conservative force, i.e. we could describe it as a potential as it is velocity independent. Now, let's assume we have a velocity dependent force. This could be a migration force in a protoplanetary disk or PR drag. We'll start from scratch and add the same two particles as before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "sim = rebound.Simulation()\n",
    "sim.integrator = \"ias15\"\n",
    "sim.add(m=1.)\n",
    "sim.add(m=1e-6,a=1.)\n",
    "sim.move_to_com()  # Moves to the center of momentum frame"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But we change the additional force to be"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "ps = sim.particles\n",
    "tau = 1000.\n",
    "def migrationForce(reb_sim):\n",
    "    ps[1].ax -= ps[1].vx/tau\n",
    "    ps[1].ay -= ps[1].vy/tau\n",
    "    ps[1].az -= ps[1].vz/tau"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to let REBOUND know that our force is velocity dependent. Otherwise, REBOUND will not update the velocities of the particles. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "sim.additional_forces = migrationForce\n",
    "sim.force_is_velocity_dependent = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we integrate as before. But this time we keep track of the semi-major axis instead of the eccentricity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA30AAAFCCAYAAABFObToAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYXVWZ7/HvmwFQBAKiyBBllCHMEESZggwGGcQJpMX5\n2nbbDte+2ordrXjVFm1taRttFRG9iNAygyAySBwQgTAGCJFRmZUpIAgk8N4/1i6qKqlA5aR27XP2\n+X6eZz91hl1nv3lWpl+ttdcbmYkkSZIkqZ0mNF2AJEmSJKk+hj5JkiRJajFDnyRJkiS1mKFPkiRJ\nklrM0CdJkiRJLWbokyRJkqQWm9R0AaMREfaVkCRJktTXMjM6+b6eCH3Q+S9Q3S8iDs/Mw5uuQ2PP\nsW03x7e9HNt2c3zby7Ftt2WZCHN5pyRJkiS1mKFPkiRJklrM0KduMKvpAlSbWU0XoFrNaroA1WZW\n0wWoVrOaLkC1mdV0AepOkdn9e6RERHpPnyRJkqR+tSyZyJk+SZIkSWoxQ58kSZIktZihT5IkSZJa\nzNAnSZIkSS1m6JMkSZKkFjP0SZIkSVKLGfokSZIkqcUMfZIkSZLUYoY+SZIkSWoxQ58kSZIktZih\nT5IkSZJazNAnSZIkSS1m6JMkSZKkFjP0SZIkSVKLGfokSZIkqcUMfZIkSZLUYoY+SZIkSWoxQ58k\nSZIktZihT5IkSZJazNAnSZIkSS1Wa+iLiO9HxH0RMec5zvlGRNwUEddExDZ11iNJkiRJ/abumb5j\ngZlLejMiXg9smJkbAX8L/HfN9UiSJElSX6k19GXmr4GHnuOUA4AfVudeCkyJiDXqrEmSJEmS+knT\n9/StDdwx5PmdwDoN1SJJkiRJrdN06AOIRZ7niCcFK41DLZIkSZLUKpMavv5dwNQhz9epXhvBp++J\nmHsJnHcZPH5+Zs6qvzxJkiRJGn8RMQOYMRaf1XToOxP4EHBiROwIPJyZ94186r/tABwOvBf4cwSX\nZvLXcapTkiRJksZNNck1a+B5RHy208+KzBFXU46JiDgB2A1YHbgP+CwwGSAzv1OdcxRlh8/HgPdk\n5pUjfE5mZpTHbAV8Dtge+BLwvUyerO0XIUmSJEkNG5qJlvp76wx9Y2WkX2AE2wH/F9gc+ALwg0wW\nNFGfJEmSJNVpWUJfN2zk0pFMrshkX+BtwFuBGyN4V0TjS1YlSZIkqWv07Ezf4uewK2Xmb03K8s//\nyeTp8ahPkiRJkurUl8s7Rz6PAF4LfB5YhXIP4amZPFNziZIkSZJUG0PfYucTwOso4W8yJfydmTly\nD0BJkiRJ6maGviV+HwHsT1n2uZCy7POnhj9JkiRJvcTQ97zfzwTgQOAzwNOUEOjMnyRJkqSeYOgb\n9ecwATiAEv6CEv7O8J4/SZIkSd3M0LfUn/fsss/PApMo4e80w58kSZKkbmTo6/hzCWBfSvhbgRL+\nTjH8SZIkSeomhr5l/nwC2IcS/lak7Pp5sn3+JEmSJHUDQ9+YXefZVg+fpfT5+zzwE8OfJEmSpCYZ\n+sb8egSwFyX8rQZ8ATjR8CdJkiSpCYa+2q5LAHtQwt9LKeHvhEwWjnctkiRJkvqXoa/26xPA7pTw\ntybwb8DxmSxoqiZJkiRJ/cPQN44imAH8K7AB8BXg+5k80WhRkiRJklptWTLRhLEupu0ymZXJHsDb\nKDt+3hrBxyN4UcOlSZIkSdJiDH0dyuR3mexPCX7TKeHvMxGs2nBpkiRJkvQsQ98yyuSaTA4GdgHW\nBW6O4IgI1mi2MkmSJEky9I2ZTOZl8l5gW2AlYG4E34hgasOlSZIkSepjhr4xlskfMvkHYBrwJHBN\nBN+LYMOGS5MkSZLUhwx9Ncnknkw+AWwE3AVcEsGPI9i84dIkSZIk9RFDX80yeSCTz1JaPFwDXBDB\n6RFMb7g0SZIkSX3A0DdOMnkkky8D6wEXAqdEcF4EM6rm75IkSZI05mzO3pAIlgMOBf4JmA8cAZyR\nyTONFiZJkiSp6yxLJjL0NSyCCcAbgE8BqwBfAY7P5MlGC5MkSZLUNQx9LVAt8ZwBfBLYHDgS+E4m\njzZZlyRJkqTmLUsm8p6+LpFJZnJRJjOB/YHtgNsi+KKN3iVJkiR1ytDXhTK5KpNDgFcBq1IavX8r\ngvUbLk2SJElSjzH0dbFMbsnkg8CmwEPAZRGcEMHWDZcmSZIkqUcY+npAJvdl8s/A+sAVwNkRnBvB\n7rZ7kCRJkvRc3MilB0WwPIPtHh4BvgqcksnCRguTJEmSVAt37+xTVbuH/YCPAy8Hvg583x0/JUmS\npHZx984+lckzmZyZya7AwcBOlB0/vxTBWg2XJ0mSJKkLGPpaIpNLMzkI2AFYEbgugmMj2Lzh0iRJ\nkiQ1yNDXMpncmslHgA2Bm4HzI/hZBHu46YskSZLUf7ynr+UiWAF4O/B/gCcpm778JJMFjRYmSZIk\nadTcyEXPq9r0ZR/Kpi8bAP8JHJ3JI40WJkmSJOl5uZGLnle16cvZmewOvAnYnrLpy79HMLXh8iRJ\nkiTVxNDXhzKZnckhwHbAJOCaCI6LYOuGS5MkSZI0xgx9fSyT2zP5GGW553XATyO4KIIDIpjYcHmS\nJEmSxoD39OlZEUwG3gp8DFiVct/fD2z2LkmSJDXLe/o0JjJZkMmPKb3+3gXsCtwewVcjWLfR4iRJ\nkiR1xNCnxWSSmVycyVuBbauXr4jgpAh2st+fJEmS1Dtc3qlRiWAl4N3AR4EHgSOBk+z3J0mSJNXP\nPn0aN9UGL/sC/xt4JfBN4LuZPNBoYZIkSVKLeU+fxk0mT2dyZiavBfajBL+bI/h2BJs2XJ4kSZKk\nRRj61LFMrs7kPcCmwL3ARRH8LIJ9Ivy9JUmSJHUDl3dqzESwAnAI8GHgRZSlnz/IZH6jhUmSJEk9\nznv61FWq3T1fQwl/ewMnAEdlMrfRwiRJkqQe5T196ipDWj68DdiCstvnRRGcF8H+1WYwkiRJksaB\nM30aFxEsDxwEfARYjbL08/uZPNxoYZIkSVIPcKZPXS+TJzM5DtgBeDuwHXBbtevntGarkyRJktrL\n0KdxVS39/F0mbwc2A+4BLojgwggOdOmnJEmSNLZqDX0RMTMiboyImyLikyO8v3pEnBsRV0fEdRHx\n7jrrUXfJ5J5MPge8AjgG+CRwSwT/FMFqzVYnSZIktUNt9/RFxERgHrAncBdwOXBIZs4dcs7hwPKZ\neVhErF6dv0ZmLlzks7ynr09EMJ2y6+f+wCnANzO5qtmqJEmSpGZ16z19OwA3Z+btmbkAOBF4wyLn\n3AOsXD1eGXhg0cCn/pLJ5Zm8E9gYuA04I4JLInhH1QdQkiRJ0lKoM/StDdwx5Pmd1WtDHQ1Mi4i7\ngWuAj9ZYj3pIJn/K5IvA+sCXKJu//DGCL0ewfrPVSZIkSb2jztA3mnWjnwauzsy1gK2Bb0bESjXW\npB6TycJMzsxkJqXh+0TgsgjOjmBfN36RJEmSntukGj/7LmDqkOdTKbN9Q70G+CJAZt4SEbdRlvXN\nXvTDqvv/BszKzFljWay6XyY3Ax+P4F+Bg4HPAkdF8G1Kz78/N1qgJEmSNEYiYgYwY0w+q8aNXCZR\nNmbZA7gbuIzFN3L5D2B+Zn4uItYArgC2zMwHF/ksN3LRiKqNX/4eeCNwNvAt4JLMUc00S5IkST1h\nWTJRbaEPICL2AY6kLMk7JjO/FBEfAMjM71Q7dh4LvJyy1PRLmfnjET7H0KfnVLV4eBfwQeAxSvj7\ncSZ/abQwSZIkaQx0begbK4Y+jVYEEyizyx8EdgWOB76TyfWNFiZJkiQtg25t2SCNu0yeyeT8TN5I\n2RxoPnB+BL+p2j68oOESJUmSpHHlTJ9aL4LJwL7AB4DpwI+A72ZyQ6OFSZIkSaPkTJ/0HDJZkMnp\nmewDbA/8Bbgggl9HcKizf5IkSWozZ/rUl6rZv/0os3/bA8dRZv/mPuc3SpIkSQ2odaYvIv49IlaO\niMkRcWFE3B8R7+jkYlK3qGb/Tquavk8HHgd+EcGvInh7BCs0XKIkSZI0JkazvHPvzHyEMityO7AB\n8Ik6i5LGUya3ZfLPlNYhRwLvBO6I4D8i2KTZ6iRJkqRlM5rQN6n6uh9wcmbOBxtfq32q2b9TM3kd\n8CrgCeCiCH4Zwd84+ydJkqReNJrQd1ZE3AhsB1wYES+l/GdYaq1Mbs3k08ArgP8C3k2Z/fuas3+S\nJEnqJaPayCUiXgw8nJlPR8SKwEqZeW/t1Q1e341c1LgINgDeD7wLuBX4HnBSJn9ptDBJkiS13rJk\noiWGvojYIzMvjIg3M7icc+AimZmndnLBThj61E2qnT9fD7wP2AU4GTgGuDTTpc+SJEkae8uSiSY9\nx3u7AhcC+zPyPXzjFvqkbpLJAuAM4IwI1qTM/B0HPBXB94DjMrm/yRolSZKkAfbpk8ZABEH5Qcn7\ngAOA8ymzf+dn8nSTtUmSJKn31d2n70cRMWXI83Uj4hedXExqq0wyk19m8k7K5i8XAp8HbovgcxGs\n22iBkiRJ6luj2b3z18ClEbFvRPwtcB7w9XrLknpXJvMz+XYm0ynLo1cFZkdwfgRvs/WDJEmSxtNo\nd+/cBfgFcD+wbWbeU3dhi1zf5Z3qaVXQO5Cy/HMb4HjgmEyubbQwSZIk9YS6l3e+A/g+8E7gB8A5\nEbF1JxeT+lUmT2RyYiZ7AdOB+cDZEVwewd9FsGrDJUqSJKmlnnemLyJOB/42M/9UPd8B+G5mjlvw\nc6ZPbRTBRGAv4L3A3sC5lB+suPmLJEmShqmlT9/zXHC5zHyqkwt2wtCntotgNeBg4N3AVEoLiB9m\nckOTdUmSJKk71Br6IuIFlPuQpsGzG1BkZr63kwt2wtCnfhLBZpTef+8A7qTM/p2YyYNN1iVJkqTm\n1HpPH2XGYQ3gdcAsYB3gL51cTNLzy+SGTD4JvBz4DKX/320R/CSC10cwqdkKJUmS1EtGM9N3dWZu\nHRHXZuaWETEZ+E1mvmp8SnSmT6o2ejmYMgO4LvAjyvLP65qsS5IkSeOj7pm+gXv35kfEFsAU4CWd\nXExSZzJ5qOr992pgd2Ah8PNq988PRfDihkuUJElSlxpN6Ds6IlYD/gU4E7gB+EqtVUlaokxuzOQw\nyvLPfwF2Am6N4OQI9o9gcrMVSpIkqZt0tHvneHN5p/TcIpgCHETZ/XN9SvP3H9r8XZIkqR3GvWXD\neDP0SaMXwcaUe/8OBR6mbMZ0QiZ3NlqYJEmSOmbok7SYCCZQdv48FHgTcCVlA5hTMnm0ydokSZK0\ndGoLfRExAdgxM3/baXFjwdAnLZsIXgDsR+n9txtwDmUG8PxMFjRZmyRJkp5f3c3Zr87MrTuqbIwY\n+qSxE8HqlPYPh1Lu/zuRMgM4O5Pun/qXJEnqQ3W3bLggIt4SEYYuqQUyuT+Tb1btH3YGHgJOAOZG\n8C8RrNdshZIkSRpLo5np+wvwQuBp4Inq5czMlWuubWgNzvRJNYoggB0pyz8PAuZSln+elMlDTdYm\nSZIkN3KRNIYiWA7Yh7L8c2/gAsryz3MyebLJ2iRJkvpV7aEvIt5A2QUwgV9m5lmdXKxThj6pGVX/\nv7dQZgC3AE6jLAW9KJOnm6xNkiSpn9S9kcsRwHRKs+cA3gbMzszDOrlgJwx9UvMiWIeyAczfAGsB\nPwF+DFzmBjCSJEn1qjv0zQG2zsynq+cTgaszc4tOLtgJQ5/UXaoG8IdUxyTK7N8JmVzfaGGSJEkt\nVffunQlMGfJ8SvWapD6VybxMDgc2Ad4KrAD8PIJrIvhUBOs2WZ8kSZIGjWam7xDgCGBW9dJuwKcy\n88R6SxtWgzN9UpeLYAKwC2X27y3APMryz5My+VOTtUmSJPW68djIZS3KfX0JXJaZ93ZysU4Z+qTe\nUu0Auhfl/r99gd9RloCelskjTdYmSZLUi2oJfRGxaWbOjYjtKGFv4AJZHQ9m5h86uehSF2nok3pW\nBCsC+1NmAGcA51NmAM/JfLb3pyRJkp5DXaHv6Mx8f0TMYuR7+F4MXJuZh3Zy4aVh6JPaIYLVgDdR\nZgC3Ac4G/gc4zx6AkiRJS9ZYc/aIOC8z9+74A0Z/HUOf1DIRvAx4M3AQpQfgWZQAeEEmTzVZmyRJ\nUrcZj3v6tgA2pezQB0Bm/r9OLtgJQ5/UbhGsRdn85SDK3zVnUALgLzJZ0GRtkiRJ3aDuPn2HU3bs\nnEZZirUP8JvMfEsnF+yEoU/qHxFMZTAAbgicRmkEPyuThU3WJkmS1JS6Q991wFbAlZm5VUSsARyf\nmXt2csFOGPqk/hTBKyh9AA8CXgGcSgmAv8rk6SZrkyRJGk91N2f/a2Y+DSyMiFWAPwFTO7mYJC2N\nTP6QyVcz2QF4NXA78FXgrgiOimDXqj+gJEmSlmA0/1m6PCJWBY4GZgNXAb+ttSpJWkQmt2by5Uy2\nozSBvxv4BnBHBP8ZwU4GQEmSpMUt1e6dEbEesFJmXltfSSNe1+WdkkYUwSaUJaAHA6sAJwGnAJdk\n8kyTtUmSJI2V8di9cytgXWAipUl7ZuapnVywE4Y+SaMRwTRKAHwzpZfoqZQA+Gs3gZEkSb2s7o1c\njqX00LoeBn9qnpnv6eSCnTD0SVpaEWxMCX9vptyHfAZwMnCRfQAlSVKvqTv03QBMy2Xp4r6MDH2S\nlkUE6wFvogTAjYGfUmYAz8vkiSZrkyRJGo26d++8HNiskw+XpG6QyW2ZfC2T11Ba0MwG/hG4N4IT\nI3hLBCs2W6UkSVI9RjPTNwM4E7gXeLJ6OTNzy3pLG1aDM32SxlwELwUOpDSDfxVwIWUG8KeZzG+y\nNkmSpKHqXt55C/Ax4DqG39N3eycX7IShT1LdIlgNOICyBHQ34FeUAHhmJg80WZskSVLdoe+SzHx1\nR5WNEUOfpPEUwcrAfpQAuCdwGYMB8O4ma5MkSf2p7tD3LWAKcBY8u+PdqFo2RMRM4EhKq4fvZeaX\nRzhnBvB1YDJwf2bOGOEcQ5+kRlT3+s2kBMB9gBuB04HTM5nXZG2SJKl/1B36fgAsdtLztWyIiInA\nPMpPye+ibAhzSGbOHXLOFOBi4HWZeWdErJ6Z94/wWYY+SY2LYDnK0s83Am8A5lMFQGC2zeAlSVJd\nam/O3tEHR7wa+GxmzqyefwogM48Ycs4HgZdl5mee57MMfZK6SgQTgO0pAfBAYCVKL8DTgVmZLGiw\nPEmS1DJ1t2wYeqErl+L0tYE7hjy/s3ptqI2A1SLiooiYHRHvWJp6JKkpmTyTyWWZHJbJppRVDX8E\nPg/cF8GPqlYQL2q2UkmS1O+WKvQBS5MsRzOFOBnYFng98DrgXyNio6WsSZIal8mNmXw5kx2BzYHf\nAP8LuDuCsyJ4X9UiQpIkaVxNWsrzz1mKc+8Cpg55PpUy2zfUHZTNW/4K/DUifkVpnHzToh8WEYcP\neTorM2ctRS2SNG6qHT6/DXw7gimUDWAOBL4WwRzgNMpGMLc2WKYkSepi1YaXM8bks2q8p28SZSOX\nPYC7KVueL7qRyybAUZRZvuWBS4GDM/OGRT7Le/ok9bwIlqf8nXggpSfg/cCZ1XGZG8FIkqQlqWUj\nl4i4ODN3ioi/sPhSzczMlUdR2D4Mtmw4JjO/FBEfqD7gO9U5HwfeQ2n8fnRmfmOEzzH0SWqVaiOY\nHYD9KQHwpcBPKe1xzs/ksQbLkyRJXaYrd+8cS4Y+SW0XwXoMBsAdgF9RAuBZNoSXJEm1h76IWJVy\nT96z9wBm5tLs5LlMDH2S+kl1H+BMSgjcB7iFsgT0LOCazFFtlCVJklqk7ubsnwfeDdwKg/ebZObu\nnVywE4Y+Sf0qgsnAzpQZwAMoux4PBMBZmTzZYHmSJGmc1B36fg9snplPdXKBsWDokySIIIBNKeFv\nf2AacAElBJ6Tyf0NlidJkmpUd+g7Dfi7zLyvkwuMBUOfJC2u6vu3LyUEvha4Hji7OlwGKklSi9Qd\n+qYDZwDXwbPLiDIzD+jkgp0w9EnSc6vaQexGCYH7UdrgnEMJgBe4G6gkSb2t7tA3F/hvSugbuKcv\nM/OXnVywE4Y+SRq9ahnoxpQAuC8wHfgtpSXE2TaFlySp99Qd+i7PzOkdVTZGDH2S1LkIVgH2ogTA\nfYCHGFwG+ptMFjRYniRJGoW6Q99/UJZ1nsng8k5bNkhSD6qawm/H4CzghpTNYM4GfpZJY/dvS5Kk\nJas79M2CxTcDsGWDJPW+CF5Gmf3bF9gT+D3VMlDgqszBVj2SJKk5tTdnb5qhT5LqF8FylJ6AA7OA\nU4Bzq+O8TB5ssDxJkvpa3TN9LwO+CKydmTMjYjPg1Zl5TCcX7IShT5LGXwQbAK+jzATuRmkJcS7w\nM+CKTJ5usDxJkvpK3aHvXOBY4J8zc8uImAxclZmbd3LBThj6JKlZVUuInSkBcCbwMuA8SgA8z3sB\nJUmqV92hb3Zmbh8RV2XmNtVrV2fm1p1csBOGPknqLhFMpYS/mcAewM0MLgX9XSYLGyxPkqTWWZZM\nNGEU5/wlIl485GI7AvM7uZgkqR0yuSOTozN5M/AS4B8p/6b8F/DnCE6K4H0RrN1ooZIkaVQzfdtR\n/hGfRrmf4yXAWzLzmvrLe7YGZ/okqUdEsCawN2Up6F7AXZRloOcCF2fyVIPlSZLUk2rfvbO6j2/j\n6um8zBzXRr6GPknqTRFMBKYzeC/gJsBFDO4IemuD5UmS1DPqvqfvIODczHwkIv4V2Ab4gs3ZJUlL\nK4LVKbN/A7OAjwHnUzaFuSiThxssT5KkrlV36JuTmVtExM7AF4CvAp/JzB06uWAnDH2S1D4RBLA5\nZSnoXsBOwBxKADwPuMwNYSRJKuoOfVdn5tYRcQQwJzOPH7qT53gw9ElS+0WwAqUtxEAIXBeYRRUC\nM7mlseIkSWpY3aHvbMpN+HtRlnY+AVyamVt1csFOGPokqf9EsAawJ4Mh8AkGZwF/4VJQSVI/qTv0\nrUi5+f7azLwpItYEtsjM8zq5YCcMfZLU36qloNMoAXBvylLQ6xi+FHRcNxmTJGk81b57Z9MMfZKk\noaqloDsxGALXA34JXABcCMzNpPv/gZMkaZQMfZKkvhbBSylLQfeojuWAX1AC4IWZ/LHB8iRJWmaG\nPkmSKtVS0PUZDICvBR6iCoCU1hAPNFehJElLz9AnSdISRDAB2JLBELgzcDODIfDXmTzWXIWSJD0/\nQ58kSaMUwXLADgyGwG2BKxm8H9BNYSRJXcfQJ0lShyJYEdiFwRC4AfAbBmcC52TyTHMVSpJk6JMk\nacxEsDqwO4MhcBXgIgZD4K3uDCpJGm+GPkmSahLByxm+KcxCYBYlCM7K5LbmqpMk9QtDnyRJ46Da\nGXQjykzgjOrrEwwPgX9oqj5JUnsZ+iRJakAVAjdmMADOAB6jhMBZlPYQdzRTnSSpTQx9kiR1gSoE\nbkoJfwPHowyfCbyzmeokSb3M0CdJUheqQuBmDM4C7gY8zPAQeHdT9UmSeoehT5KkHlA1ip/G4HLQ\n3YAHqAIgJQTe01R9kqTuZeiTJKkHVSFwc4bPBN4P/GrI8QdbREiSDH2SJLXAkJnAXYccCxkeAm80\nBEpS/zH0SZLUQtU9gRsyPASuCPyawRB4bSZPN1akJGlcGPokSeoTVbP4XapjV2At4GIGQ+AVmTzV\nXIWSpDoY+iRJ6lMRvBTYmcGZwI2AyxgMgZdm8nhzFUqSxoKhT5IkARDBFOA1DIbALYFrKUtCfwP8\nNpMHmqtQktQJQ58kSRpRBC8EXg3sRJkRfBVwF2VJ6MWUIHiLm8NIUncz9EmSpFGJYBKwBYMhcCdg\nMsND4FWZLGisSEnSYgx9kiSpY9XmMAMBcGdgfWA2gyHwkkzmN1ehJMnQJ0mSxkwEqzB8Sej2wK0M\nhsCLgT+6JFSSxo+hT5Ik1SaCycDWDJ8NXMDwJaH2C5SkGhn6JEnSuKmaxq/P8BC4NnApJQReQmkV\n4ZJQSRojhj5JktSoCF7M4JLQVwPbAX+gBMCBY14mzzRWpCT1MEOfJEnqKtWS0C0oAXDgWI0yGzgQ\nAp0NlKRRMvRJkqSuF8EawI4MhsDtgNtxNlCSnpehT5Ik9ZxqNnBLhs8GTmHx2cBHGitSkrqEoU+S\nJLVCBC9j+GzgtsBtDIbA3+FsoKQ+ZOiTJEmtVM0GbsXis4GXU2YEL6XMBv65sSIlaRwY+iRJUt+o\n7g3cAXhVdUwHHgQuYzAIXpXJXxsrUpLGWNeGvoiYCRwJTAS+l5lfXsJ50ylLNg7KzFNHeN/QJ0mS\nRhTBBOCVlAA4EAY3BW5kyGwg8HuXhUrqVV0Z+iJiIjAP2BO4i7IM45DMnDvCeecDjwPHZuYpI3yW\noU+SJI1aBC8AtmH4jOBqDF8Welkm9zVWpCQthWXJRJPGupghdgBuzszbASLiROANwNxFzvswcDJl\naYYkSdIyq5Z2/rY6AIjgJQyGwA8BO0QwnyoAVl+vzOTx8a9YkupTZ+hbG7hjyPM7KX/JPisi1qYE\nwddSQl/332AoSZJ6UrXZy9nVQQQBbMTgstCDgWkR3ATMpswKzgbmZPJUI0VL0hioM/SNJsAdCXwq\nMzMiAnAJpyRJGheZJPD76jgOIILlKb0DpzM4I7h+BNczGAIvB+Zm8nQTdUvS0qoz9N0FTB3yfCpl\ntm+o7YATS95jdWCfiFiQmWcu+mERcfiQp7Myc9aYVitJkvpeJk9SQt3lA69FsCLl/sDtKXsVHAas\nGcHVDA+Ct7hRjKSxEhEzgBlj8lk1buQyibKRyx7A3ZS18ott5DLk/GOBs9y9U5IkdbsIplB+eL09\nZVZwe2AV4AqGLw39YzWjKEnLpCs3csnMhRHxIeDnlJYNx2Tm3Ij4QPX+d+q6tiRJUp0yeRi4sDoA\niOCllCA4HXgXcBQwMWJYCLw8k3vHv2JJ/czm7JIkSTWoNopZixICB2YDt6e0qboCuHLgayb3NFWn\npN7QlX36xpKhT5IktUEVBNcHtq2O7aqvCxgMggPHHS4NlTTA0CdJktSjqiA4lcEAOBAGJzI8BF4J\n3GoQlPqToU+SJKllIliLwRA4cKwEXMVgCLwCuMldQ6X2M/RJkiT1gWqzmG0YPiP4EuBqhtwjCNyY\nycKm6pQphyFNAAAK90lEQVQ09gx9kiRJfSqCVRkMggNLRNcB5jAYAq8Crq/6EErqQYY+SZIkPSuC\nlYGtGAyCWwMbAjdTZgUHjmsyeaCpOiWNnqFPkiRJzymCFYDNKAFw4NgKeIThQfBq4DbvE5S6i6FP\nkiRJSy2CCcC6DA+CWwNTgGsYHgSvz+SJZiqVZOiTJEnSmIngxZRZwKFBcCNGXh56f1N1Sv3E0CdJ\nkqRaDVkeumgYHLo8dGB28FaXh0pjy9AnSZKkcVc1ll+X4SFwG2BV4Hrg2uqYA8zJ5MFmKpV6n6FP\nkiRJXSOCKcAWwJaLfH2YwRA4EAjnZbKgoVKlnmHokyRJUlerNo15BSUADg2DrwB+z+Jh8J5Muv8/\nqtI4MfRJkiSpJ0XwAsq9gkPD4FbABAYD4EAgvD6TxxoqVWqUoU+SJEmtUd0ruAaDs4EDYXAT4E6G\nzwjOwY1j1AcMfZIkSWq9CCZTWkcMDYNbAi+mbBxzPXBddVwP3O0SUbWFoU+SJEl9q9o4Zlp1bF4d\n04DlGAyAz37N5M8NlSp1zNAnSZIkLSKCl7B4ENwceIrFZwWvz+ThhkqVnpehT5IkSRqF6n7BNVk8\nCG4GzGd4ELwOuMHNY9QNDH2SJEnSMqhaSryc4UFwc2Bj4F4WXyY6L5O/NlOt+pGhT5IkSapBBBOB\nDRgeBDcDNgTuBm6ojrkDXzN5tJlq1WaGPkmSJGkcRTCJEgY3pYTAzarHmwAPsngYvCGTB5upVm1g\n6JMkSZK6QLVM9BUsHgY3A55ghDAI3GdrCT0fQ58kSZLUxaoNZNZieAgcOCYychi8wzCoAYY+SZIk\nqUdVrSUWnRXcDFgJmAfcOOSYB9yUyRPNVKumGPokSZKklqmazm9K2UF0Y8r9gpsA6wF3MTwQDjz+\nk7OD7WTokyRJkvpEBJMpwW8gBA4NhBMYHgIHHt+SyVONFKwxYeiTJEmSRASrM3IYnAr8gRGWi2by\nQDPVamkY+iRJkiQtUQTLU1pMjBQIn2Lk2cHbMlnYSMFajKFPkiRJ0lKrdhVdg5HD4JrAbcDvhxzz\nqq+2mRhnhj5JkiRJYyqCF1BmB19ZHRsPebw8w0Pgs0cmjzZScMsZ+iRJkiSNmwhWYzAADg2EGwE/\nz+SNDZbXSoY+SZIkSY2LYAKwciYPN11L2xj6JEmSJKnFliUTTRjrYiRJkiRJ3cPQJ0mSJEktZuiT\nJEmSpBYz9EmSJElSixn6JEmSJKnFDH2SJEmS1GKGPkmSJElqMUOfJEmSJLWYoU+SJEmSWszQJ0mS\nJEktZuiTJEmSpBYz9EmSJElSixn6JEmSJKnFDH2SJEmS1GKGPkmSJElqMUOfJEmSJLWYoU+SJEmS\nWszQJ0mSJEktZuiTJEmSpBarPfRFxMyIuDEiboqIT47w/tsj4pqIuDYiLo6ILeuuSZIkSZL6Ra2h\nLyImAkcBM4HNgEMiYtNFTrsV2DUztwQ+D3y3zprUfSJiRtM1qB6Obbs5vu3l2Lab49tejq2WpO6Z\nvh2AmzPz9sxcAJwIvGHoCZl5SWbOr55eCqxTc03qPjOaLkC1mdF0AarVjKYLUG1mNF2AajWj6QJU\nmxlNF6DuVHfoWxu4Y8jzO6vXluR9wDm1ViRJkiRJfWRSzZ+foz0xInYH3gvsVF85kiRJktRfInPU\nuWzpPzxiR+DwzJxZPT8MeCYzv7zIeVsCpwIzM/PmET6nviIlSZIkqQdkZnTyfXXP9M0GNoqIdYG7\ngYOBQ4aeEBEvpwS+Q0cKfND5L06SJEmS+l2toS8zF0bEh4CfAxOBYzJzbkR8oHr/O8BngFWB/44I\ngAWZuUOddUmSJElSv6h1eackSZIkqVm1N2dfFs/X2F3dLyK+HxH3RcScIa+tFhHnR8TvI+K8iJgy\n5L3DqvG+MSL2bqZqjUZETI2IiyLi+oi4LiI+Ur3u+LZARKwQEZdGxNURcUNEfKl63fFtiYiYGBFX\nRcRZ1XPHtiUi4vaIuLYa38uq1xzfloiIKRFxckTMrf5+fpXj2/siYuPqz+zAMT8iPjJWY9u1oW+U\njd3V/Y6ljOFQnwLOz8xXAhdWz4mIzSj3fW5Wfc+3IqJrf4+KBcDHMnMasCPwD9WfUce3BTLzCWD3\nzNwa2BLYPSJ2xvFtk48CNzC407Zj2x4JzMjMbYbcMuP4tsd/Audk5qaUv59vxPHteZk5r/ozuw2w\nHfA4cBpjNLbdPOjP29hd3S8zfw08tMjLBwA/rB7/EDiwevwG4ITMXJCZtwM3U34fqAtl5r2ZeXX1\n+C/AXEofTse3JTLz8erhcpT7sh/C8W2FiFgHeD3wPWBgszTHtl0W3QTP8W2BiFgF2CUzvw9l/4zM\nnI/j2zZ7UnLQHYzR2HZz6Fvaxu7qHWtk5n3V4/uANarHa1HGeYBj3iOqHXq3AS7F8W2NiJgQEVdT\nxvGizLwex7ctvg58AnhmyGuObXskcEFEzI6I91evOb7tsB7w54g4NiKujIijI2JFHN+2eRtwQvV4\nTMa2m0OfO8z0gSw7CT3XWPv7oMtFxIuAU4CPZuajQ99zfHtbZj5TLe9cB9g1InZf5H3HtwdFxH7A\nnzLzKhafDQIc2xbYqVoitg9l6f0uQ990fHvaJGBb4FuZuS3wGNVyvwGOb2+LiOWA/YGTFn1vWca2\nm0PfXcDUIc+nMjzNqnfdFxEvA4iINYE/Va8vOubrVK+pS0XEZErgOy4zT69ednxbplo6dDblHgPH\nt/e9BjggIm6j/CT5tRFxHI5ta2TmPdXXP1PuCdoBx7ct7gTuzMzLq+cnU0LgvY5va+wDXFH9+YUx\n+rPbzaHv2cbuVeI9GDiz4Zo0Ns4E3lU9fhdw+pDX3xYRy0XEesBGwGUN1KdRiIgAjgFuyMwjh7zl\n+LZARKw+sENYRLwA2Au4Cse352XmpzNzamauR1lC9IvMfAeObStExAsjYqXq8YrA3sAcHN9WyMx7\ngTsi4pXVS3sC1wNn4fi2xSEMLu2EMfqzW2tz9mWxpMbuDZelpRQRJwC7AatHxB3AZ4AjgJ9ExPuA\n24GDADLzhoj4CWU3uYXAB9NGkt1sJ+BQ4NqIuKp67TAc37ZYE/hhtRPYBMps7oXVWDu+7TIwTv7Z\nbYc1gNPKz+WYBByfmedFxGwc37b4MHB8NSlyC/Aeyv+VHd8eV/2gZk/g/UNeHpO/m23OLkmSJEkt\n1s3LOyVJkiRJy8jQJ0mSJEktZuiTJEmSpBYz9EmSJElSixn6JEmSJKnFDH2SJEmS1GKGPklS34qI\nVSLi76vHa0bESU3XJEnSWLNPnySpb0XEusBZmblFw6VIklSbSU0XIElSg44ANoiIq4CbgE0zc4uI\neDdwIPBCYCPga8AKwN8ATwKvz8yHImID4CjgJcDjwPszc974/zIkSVoyl3dKkvrZJ4FbMnMb4BOL\nvDcNeCMwHfgi8EhmbgtcAryzOue7wIczc/vq+781LlVLkrQUnOmTJPWzWMJjgIsy8zHgsYh4GDir\nen0OsGVErAi8Bjgp4tlvXa7OYiVJ6oShT5KkkT055PEzQ54/Q/n3cwLwUDVLKElS13J5pySpnz0K\nrLSU3xMAmfkocFtEvAUgii3HuD5JkpaZoU+S1Lcy8wHg4oiYA3wFGNjSOoc8ZoTHA8/fDrwvIq4G\nrgMOqLdiSZKWni0bJEmSJKnFnOmTJEmSpBYz9EmSJElSixn6JEmSJKnFDH2SJEmS1GKGPkmSJElq\nMUOfJEmSJLWYoU+SJEmSWszQJ0mSJEkt9v8Bh4vTbEqwOfcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x108c955c0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Nout = 1000\n",
    "a_s = np.zeros(Nout)\n",
    "times = np.linspace(0.,100.*2.*np.pi,Nout)\n",
    "for i, time in enumerate(times):\n",
    "    sim.integrate(time)\n",
    "    a_s[i] = sim.particles[1].a \n",
    "fig = plt.figure(figsize=(15,5))\n",
    "ax = plt.subplot(111)\n",
    "ax.set_xlabel(\"time\")\n",
    "ax.set_ylabel(\"semi-major axis\")\n",
    "plt.plot(times, a_s);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The semi-major axis decays exponentially on a timescale `tau`.\n",
    "\n",
    "In the above example, REBOUND is calling a python function at every timestep. This can be slow. Note that you can also set `rebound.additional_forces` to a c function pointer. This let's you speed up the simulation significantly. However, you need to write your own c function/library that knows how to calculate the forces. Or, you use Dan Tamayo's new migration library (in preparation)."
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}