three-body.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "from scipy.integrate import odeint  #the differential equation integrator\n",
    "from mpl_toolkits.mplot3d import Axes3D  #Import 3-D plotting tool (might need later)\n",
    "from matplotlib import animation    # animation package (for later)\n",
    "import numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Mechanics Lab 7 – Numerical Solution to the Three-Body Problem\n",
    "\n",
    "Today we will be looking at a problem that *cannot* be solved analytically –  the three-body problem.  We already know (or will learn soon) that two  massive objects moving within their mutually described Newtonian gravitational potential will undergo orbits shaped like conic sections – i.e. circles, ellipses, parabolas, and hyperbolas.  But what are the shapes of orbits in a three-body system?  The answer here is not straightforward and requires a computer to solve precisely, which is why a system of three bodies interacting gravitationaly is one of the classic examples of *chaos*.  \n",
    "\n",
    "\n",
    "## Objectives\n",
    "\n",
    "In this lab you will \n",
    "* create an integrator to find the motion of 2 and 3-body gravitational systems;\n",
    "* test your integrator by applying it to stable situations where you can predict the orbit;\n",
    "* transform the trajectories to the center-of-mass frame;\n",
    "* create a movie of the motion;\n",
    "* and observe chaotic 3-body trajectories.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparation\n",
    "\n",
    "#### Exercise 7.1 \n",
    "\n",
    "First, we'll break the problem down into manageable pieces by defining a function. Each time we invoke the gravitational force law, we are going to need to find the distance between two objects.  Fix the code below to define a function that returns the distance between two objects with $x$, $y$, and $z$ coordinates $(x_{1}, y_{1}, z_{1})$ and $(x_{2}, y_{2}, z_{2})$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def findDist(x1,y1,z1,x2,y2,z2):\n",
    "    \n",
    "    dist = ((x2-x1)**2 + (y2-y1)**2 + (z2-z1)**2)**.5 \n",
    "    \n",
    "    \n",
    "    return dist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A warmup:  Solving the 2-body problem numerically for a circular orbit\n",
    "\n",
    "Let's ease our way into the more complex 3-body problem by first generating the orbital trajectory for a 2-body system undergoing a **circular** orbit.  We'll start with a large mass (say a star) with $M_{1} = 2000$ and a smaller mass (say a planet) with $M_2 = 1$. We'll work in a unit system where $G = 1.0$.  \n",
    "\n",
    "We will confine the motion of both bodies to the x-y plane for now.  Later if we have time we will expand our solutions into 3-dimensional space.  \n",
    "\n",
    "#### Exercise 7.2  \n",
    "In the following three cells follow the prompts to fill in the appropriate initial conditions, the deriv function, and the integration call. Choose initial conditions such that the lower mass body, $M_2$, orbits in a circle 100 distance units away from the larger body, $M_1$. Integrate the trajectory from 0 to 200 with 100000 time steps.\n",
    "\n",
    "In creating the rules for your integrator, it may help to think of the gravitational force on $M_1$ by $M_2$ as:\n",
    "$$ \\vec{F}_{G,1-2} = -\\frac{G M_{1} M_{2} \\vec{r}}{(r+\\epsilon)^{3}}$$\n",
    "\n",
    "where $\\vec{r}$ is the position of $M_1$ minus the positon of $M_2$ and $\\epsilon$ is a small number (0.1 is a good value) that will ensure that our masses never experience any extremely close encounters.  (However, in calculating your initial conditions, you should not include $\\epsilon$.) \n",
    "\n",
    "*Coding style tip: rather than doing any necessary calculations to find the initial velocity of $M_2$ elsewhere (say, a piece of scratch paper), include those calculation in the cell block. This way, you'll have a record of what you did and it will be easy to reproduce the calculations for different values, for instance a larger radius circle.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.47213595499958\n"
     ]
    }
   ],
   "source": [
    "# Defining constants\n",
    "gGrav = 1 ### FIX\n",
    "epsilon = .1 ### FIX\n",
    "\n",
    "# Masses\n",
    "m1 = 2000 ### FIX\n",
    "m2 = 1 ### FIX\n",
    "\n",
    "#Initial Conditions\n",
    "\n",
    "\n",
    "#### UNCOMMENT AND FILL IN initial conditions below\n",
    "\n",
    "\n",
    "r1_0 = array([0, 0])\n",
    "drDt1_0 = array([0, 0])\n",
    "\n",
    "r2_0_r = findDist(r1_0[0], r1_0[1], 0, 0, 100, 0)\n",
    "r2_0_v = (gGrav*m1/r2_0_r)**.5\n",
    "print(r2_0_v)\n",
    "r2_0 = array([0, r2_0_r])\n",
    "drDt2_0 = array([r2_0_v, 0])\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def deriv_2body_2D(f,t): \n",
    "    \n",
    "    x_1 = f[0]\n",
    "    y_1 = f[1]\n",
    "    \n",
    "    dxDt_1 = f[2]\n",
    "    dyDt_1 = f[3]\n",
    "    \n",
    "    x_2 = f[4]\n",
    "    y_2 = f[5]\n",
    "    \n",
    "    dxDt_2 = f[6]\n",
    "    dyDt_2 = f[7]\n",
    "    \n",
    "    #Derivative definitions\n",
    "    distance1 = findDist(x_1, y_1, 0, x_2, y_2, 0)\n",
    "    # Students fill in below:\n",
    "    \n",
    "    d2xDt2_1 = -(gGrav*m2/(distance1+epsilon)**3)*(x_1 - x_2) \n",
    "    d2yDt2_1 = -(gGrav*m2/(distance1+epsilon)**3)*(y_1 - y_2)\n",
    "\n",
    "    d2xDt2_2 = -(gGrav*m1/(distance1+epsilon)**3)*(x_2 - x_1)\n",
    "    d2yDt2_2 = -(gGrav*m1/(distance1+epsilon)**3)*(y_2 - y_1)\n",
    "\n",
    "    \n",
    "    return array([dxDt_1, dyDt_1, d2xDt2_1, d2yDt2_1, dxDt_2, dyDt_2, d2xDt2_2, d2yDt2_2], float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.00000000e+00  0.00000000e+00  0.00000000e+00 ...  1.00000000e+02\n",
      "   4.47213595e+00  0.00000000e+00]\n",
      " [ 1.08845555e-14  1.99405191e-10  8.91767126e-12 ...  9.99999996e+01\n",
      "   4.47213594e+00 -3.98806383e-04]\n",
      " [ 5.67637989e-14  7.97620746e-10  3.56706842e-11 ...  9.99999984e+01\n",
      "   4.47213588e+00 -7.97612763e-04]\n",
      " ...\n",
      " [ 4.21911941e-01  9.34776348e-02  4.16099984e-03 ... -8.69552696e+01\n",
      "  -3.84986373e+00 -2.22611559e+00]\n",
      " [ 4.21920264e-01  9.34798608e-02  4.16109847e-03 ... -8.69597216e+01\n",
      "  -3.85006099e+00 -2.22577325e+00]\n",
      " [ 4.21928586e-01  9.34820864e-02  4.16119709e-03 ... -8.69641728e+01\n",
      "  -3.85025822e+00 -2.22543089e+00]]\n"
     ]
    }
   ],
   "source": [
    "# Integrate your equation of motion\n",
    "\n",
    "#### UNCOMMENT AND FILL IN code to integrate the motion from t = 0 to 200 with 100000 timesteps\n",
    "\n",
    "\n",
    "\n",
    "# initial values\n",
    "tpts = linspace(0,200,100000)\n",
    "fInit = [r1_0[0],r1_0[1],drDt1_0[0],drDt1_0[1],r2_0[0],r2_0[1],drDt2_0[0],drDt2_0[1]]\n",
    "    \n",
    "fSoln = odeint(deriv_2body_2D, fInit, tpts)\n",
    "\n",
    "\n",
    "\n",
    "print(fSoln)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plots of Trajectories\n",
    "\n",
    "#### Exercise 7.3 \n",
    "In the cell below, create plots of the trajectories (y vs. x) of both objects on the same set of axes\n",
    "Also create plots of x vs. t and y vs. t for the planet (m2).  \n",
    "If you were unable to obtain a circular trajectory of radius 100 for the planet (m2),\n",
    "de-bug your code and try again.  If you are stuck, talk with an instructor.  \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0xeaaa198>]"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD8CAYAAAB6paOMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VNXZwPHfk31PyEISspAAYd8NmwIugOJSEVoV27r2LWrVVruptYuvra9trdpaqxYrVuteFXdFXEAR2fedELaQPSH7PnPeP2YIARIIZGbuJPN8P5/5zJ0zd+Y+3Az3ufecc88RYwxKKaV8l5/VASillLKWJgKllPJxmgiUUsrHaSJQSikfp4lAKaV8nCYCpZTycZoIlFLKx2kiUEopH6eJQCmlfFyA1QF0Rnx8vMnIyLA6DKWU6lbWrl1baoxJONV63SIRZGRksGbNGqvDUEqpbkVE9ndmPa0aUkopH6eJQCmlfJwmAqWU8nGaCJRSysdpIlBKKR/nkkQgIgtEpFhEtrQpixWRxSKy2/ncy1kuIvK4iOSIyCYRGeuKGJRSSp0ZV10R/BuYeVzZPcBnxpgs4DPna4CLgSznYx7wlItiUEopdQZcch+BMeZLEck4rngWcJ5z+XlgCXC3s/wF45gjc4WIxIhIsjGmwBWxKHW6jDHUNdmoaWyhuqGFmsYWapzPjS02mm2GZpudZpudphY7TTY7zS2GLfmVLN5WRJC/H4H+QoC/H4H+fgT5CyGB/kSEBNArLIi4iCDiwoOIDQ8mLjyIhMhg+sSEkhYbSlhQt7iVR/Vw7vwVJh45uBtjCkSkt7M8BTjYZr08Z9kxiUBE5uG4YiA9Pd2NYaqeqMVmp7i6kfyKegqrGiitbqSstonSmibKahzLR55rG1uwd2Hq7iabnSYbgM1V4RMVEsCotBhGpcYwMjWakakx9I4Mxs9PXLYNpY6w4nSkvV/yCf8NjTHzgfkA2dnZXfhvqnoiYwylNU3sLa1lb2kN+8vqyK+oJ7+igUPOg7/tuKO7n0BseBBx4cHERQQxMjWG2PAgIkMCiAgOIDIkkIiQACKDA4gICSA8KIDgQD/nGb/jrD8w4OhrAezGYDPGcdXQYqfZ7rhqaGi2U93QzOG6JspqmiivdTxKa5oorm6goLKBA2V1NNns7f77qhpa+Gp3KV/tLu1wH/SNC2NKVjyTByQwITOWXuFBrtzFyoe4MxEUHanyEZFkoNhZngektVkvFch3YxyqG7PbDQcP17G9oIpdRTXkltSwt7SW3NJaqhtaWtcL8BOSokPoExPKhMxY+sSEOh8hJEeHEh8RRExYEP4uPqP2QwgAggOAYNd9b2lNI5vzKtmYV8HGgxVsyqukrLbpmHX2l9Wxv+wAL644cMLnx2fGctnIZC4cmkRiVDAieiWhOiaOqnoXfJGjjeB9Y8xw5+uHgTJjzB9F5B4g1hjzSxG5FLgduASYADxujBl/su/Ozs42OtZQz1ffZGNbQSXbCqrZXlDFjoIqdhZWU9t0tMolJSaUzPhw+iWEkxnvePSLjyClV6jLD/LezGY35BTX8NXuEr7aXcrSXSWd+tw149P57vh0hvWJ0momHyAia40x2adczxWJQERewdEwHA8UAb8D3gZeB9KBA8CVxphycZyaPIGjl1EdcKMx5qRHeU0EPY/dbsgtrWH9gQo2HHQ8dhRWt1bnRIYEMCQ5iqHJUQxOimRIchQDEyMJDfK3OHLv19BsY/meUt7fVMAHmwpobGm/+gkcV1I3npPB9WdnkNorzINRKk/waCJwN00E3V9Ti51NeRWsyC1j5d5yNhysaK3aiQwOYGRaNKOdjaPDUqLpEx2i1RkuVlTVwBtr83h55QEOVdR3uN5FwxKZN7U/Y9Ji9Kqhm9NEoCzV9sC/IrecNfvLaWh2nJkOTopkbN9ejE6LYUxaDP0TIvSAY5GD5XW8tPIATy/d0+E6M4Ym8rMLBzI4KcqDkSlX0ESgPC7vcB1Ld5WwZGcJy3NKW+v2BydFMrFfHBP7xTE+M5ZY7d3itWobW3h7wyGe/GJPh1cNt57Xn5un9iMmTP+O3k4TgXK7ZpudVXvL+WJHMUt2lZBTXAM4GnTPG5TAlKx4JmTGabfGbqyh2cZrqw/y0EfbW6/o2hqX0Yv7Lx/GsD7RFkSnTkUTgXKLhmYbX+4q4eOthXy2vZjK+maC/P2Y0C+WcwcmcN6gBPonRGj9fg+VX1HPE1/k8PLKE7usxoQF8udvj2T6kESt6vMSmgiUy9Q2tvDZjmIWbSnki53F1DXZiAoJYPrQRC4alsSUrHgdKsEHNdvsLFx3iLvf2sTxh5GQQD+e/N5Yzh/UW08KLKSJQHVJi83OspxS3l5/iEVbi6hvthEfEcyFwxK5eHgSE/vFEeivo5grB2MMK/eWc8+bm9hXVnfMewmRwfz9mjFMyIzVpOBhmgjUaTPGsCmvkoXrD/H+pnxKa5qICgng0pF9mDW6D+MyYn3qpi11ZowxrNpbzh2vrKe4uvGY96ZkxfPIlaPoHRViUXS+RROB6rTKumbeWp/HK6sOsKuohiB/P6YN6c2s0SmcPziB4AC9iUudGbvd8PmOYm5+ce0JYz/9+tIhXDcpg6AAvbJ0F00E6qSMMazZf5hXVh7gg82Ou09HpUZz9bh0Lh2RTHRYoNUhqh6mscXG00tyeezTXceU940L44WbxtM3LtyiyHouTQSqXXVNLby57hAvLN/H7uIaIoIDuGJMH+aOS2d4inYBVJ5xoKyOO15dz8aDFceU//Xq0cwa3UfbElxEE4E6Rn5FPS98s59XVh2gsr6ZkanRfH9iXy4bmaw9fpRlbHbDy6sO8Ju3txxT/p2zUnlg1jD9bXaRJgIFwIaDFTy7bC8fbi7AGMPM4UncdE4mZ/XtpWddyqtsy6/i6vnfHDO8+IiUaJ67cRzxES4c49uHaCLwYUe68j3xeQ7LckqJDA7g6nFpXH92BmmxOsKk8m5VDc3c8+YmPtxc2FoWHRrIwh+dTb+ECAsj6340EfggYwxLd5XwxOc5rNl/mPiIYOZNzeS7E/oSEayX2Kp7abbZeWzxLp5ccuyAeIvunMqgpEiLoupeNBH4EGMMX+ws5rHFu9l8qJI+0SHccl5/rspOIyRQu36q7s1uN7y0cj+/eWfrMeWf/exc+usVwkl5RSIQkUHAa22K+gG/BWKAHwJHplX6lTHmw46+RxNBx1bvK+fPH+9g9b7DpMeGcdv5/Zk9JlX7ZqsexxjD2xsOcddrG1vLQgL9+OTOc0mP0yrP9nhFIjhmQyL+wCEc01PeCNQYY/7Smc9qIjjRtvwq/vLJTj7fUUzvyGB+PC2Lq8el6bAPqsczxtHT6L6FR3sa9UsIZ+Gt5+j9L8fpbCLwZMXxNGCPMWa/9lY5c0VVDfzp4x0sXH+IyOAA7p45mBvOztApHJXPEBG+N6Ev14xL51/Lcvm/D3eQW1LLqAc+4arsVP5wxQi9Ij5NnrwiWACsM8Y8ISL3AzcAVcAa4GfGmMMdfVavCBzDPz+7bC//+CKHFpvhpsmZ3Hpufz0DUj6v2Wbn1wu38Nqag61lv79iON+fkO7zXaS9qmpIRIKAfGCYMaZIRBKBUsAAvweSjTE3HfeZecA8gPT09LP279/v9ji9kTGGRVuLePDDbRwsr+fCoYn8+tKhWieq1HEq6pr49lPL2VNS21q2+K6pZCX6bg8jb0sEs4DbjDEXtvNeBvC+MWZ4R5/31SuCA2V1/PqdLXy5q4SBiRH89rJhTM6Ktzospbza7qJqZjz2ZevrcwbE8a/rxvlk9WlnE4GnKtKuAV458kJEktu8NxvYcsInfFizzc7TS/dw4V+Xsm7/YX73raF8+OMpmgSU6oSsxEj2/fFSHpztOLf8OqeMIb/9mLfXH7I4Mu/l9isCEQkDDgL9jDGVzrL/AKNxVA3tA242xhR09B2+dEWw8WAF97y1me0FVcwYmsgDs4aRHB1qdVhKdUs1jS1cM38Fmw9VAhDoL6y+bzoxYb4xj7ZXVQ11lS8kgsYWG48u3sUzX+aSEBnM/14+nJnDk6wOS6keYf2Bw8x+cnnr6z/OGcHc8ekWRuQZmgi6kW35Vfz09Q3sKKzmmvFp3HvJEKJCtDeQUq7UYrNz52sbeH+To/IhOMCPFfdOo1d4z7068LY2AtUOm93w5JIcZv1jGWW1TTx3wzgemjNSk4BSbhDg78cT3x3L+3dMBqCxxc6Y3y9m0dbCU3yy59NEYJGCynqumb+CP3+8kxlDE1l051TOH9zb6rCU6vGGp0Sz6w8XM835/+3m/6zl+/9aSbPNbnFk1tGqIQt8sbOYn762gaYWO3+YPZwrRqf4/I0vSlnhy10lXLdgVevrJT8/j4z4njNlplYNeaEWm50/fbyDG59bTWJUCO/dMZnZY1I1CShlkakDE9jw2xmtr8/7yxJeb3OHsq/QROAhxdUNfPeZlTy1ZA/XjE/n7dvO0Uk2lPICMWFB7H3oEn44JROAX76xif95fjU2u/fXlriKJgIP2JxXyawnvmbzoUr+evVoHpozQucJUMqLiAj3XTqU1+ZNBODT7cX0/9WHVNQ1WRyZZ2gicLN3N+bznaeX4yfCG7dO4ooxKVaHpJTqwIR+caz61bTW16MfWMyuomoLI/IMTQRuYrcbHl60gx+/sp6RqdG8c/s5DOsTbXVYSqlT6B0Vws4/zGSwczrMCx/7kk+3FVkclXtpInCDphY7d72+gX98sYe549J46X8mEh8RbHVYSqlOCg7w5+M7p3Lz1H4A/M8La3jquLmTexJNBC5W3dDMjf9exTsb8rl75mAemqOTZCjVXd17yRAevWoUAH/6eAe3vbyO7tDl/nTpEcqFiqsauPqfK1iZW86jV43i1vP6a9dQpbq5OWNT+e8tkwD4YFMB33piGfYe1qNIE4GLHCyv49tPL2dfWS3P3jCOOWNTrQ5JKeUi4zJiWfLz8wDYcqiKyX/6vEfdiayJwAX2l9Uyd/4KKuuaeeWHEzl3YILVISmlXCwjPpyVzh5F+ZUNjLh/EQ3NNoujcg1NBF2UW1LD1f9cQV1TCy//cCKj0mKsDkkp5SaJUSGs+fV0ABqa7Qz+zcc9IhloIuiCPSU1XD1/Bc02Oy//cCLDU7R7qFI9XXxE8DHDUgz/3SKaWrp3NZHbE4GI7BORzSKyQUTWOMtiRWSxiOx2Pvdydxyudqiinmv/tRK73fDqvIkMSY6yOiSllIfEhAWx8XeOKdhb7IaJD33WrdsMPHVFcL4xZnSbUfDuAT4zxmQBnzlfdxulNY1c+6+VVDe28MIPxpOVGGl1SEopD4sODWSts5qovLaJGY8u7ba9iayqGpoFPO9cfh64wqI4TltlfTPXPbuK/Mp6nrthnN4trJQPi4sIZsW9jgbkfWV1XP/cqlN8wjt5IhEY4BMRWSsi85xliUcmq3c+d4sZWZpa7Nz8nzXsLq7mn9dmk50Ra3VISimLJUWH8NUvzwfgq92l/PGjHRZHdPo8kQjOMcaMBS4GbhORqZ35kIjME5E1IrKmpKTEvRF2gjGGX7+9mRW55fz5OyO1i6hSqlVabBjv3HYOAE8v3cMba/Msjuj0uD0RGGPync/FwEJgPFAkIskAzufidj433xiTbYzJTkiw/qD79NJcXl+Tx4+nZTF7jN4sppQ61qi0GJ747hgAfv7fjaw7cNjiiDrPrYlARMJFJPLIMnAhsAV4F7jeudr1wDvujKOrPt5SyJ8+3sG3RvXhrulZVoejlPJSl43sw48vGADAnCeXU1zdYHFEnePuK4JEYJmIbARWAR8YYz4G/gjMEJHdwAzna6+0p6SGn72+gdFpMTz8nZE6dpBS6qTumjGw9cbS8Q9+1i3uMXBrIjDG5BpjRjkfw4wxDzrLy4wx04wxWc7ncnfGcabqmlq49cW1BAf689T3x+qsYkqpUxKR1pnOAH7w/GoLo+kcvbO4A8YY7lu4hd3FNfxt7miSo0OtDkkp1U2EBPqz7O6jPYm8vfFYE0EHXl9zkIXrD3HX9IFMybK+sVop1b2k9gpjwQ2Oe2h//t+N5B2usziijmkiaMeBsjoeeG8bk/rFcfv5A6wORynVTV0wOJGZw5IAmPynL2jx0mEoNBEcx2Y3/PT1Dfj5CX+5ahR+fto4rJQ6c3+7ZnTr8gPvb7Mwko5pIjjOP7/cw5r9h3lg1jBSYrRdQCnVNcEB/iy603Ef7Qvf7GdnYbXFEZ1IE0EbuSU1/HXxbi4ZkcQVo1OsDkcp1UMMSorklnP7A3DRX7/E5mWD02kicHIMIbGF4EA/7r98mN4voJRyqZ9fOLB1+W+f7bYwkhNpInB6e8Mhlu8p4+6Zg+kdGWJ1OEqpHibA3691PKLHP9tNfkW9xREdpYkAx9DSf3h/O2PSY/ju+HSrw1FK9VCj0mKYMTQRgO8/u9LiaI7SRAA8uSSH8romfj9ruPYSUkq51SNXjQIgt6SW1fu8Y1AFn08EeYfreO7rfcwZk6pzDiul3C4qJJB7Lx4MwJVPf+MVs5r5fCL4y6KdCPDziwaecl2llHKFmyZnti6/uzHfwkgcfDoRbC+o4u0N+fxgcqaOJaSU8phAfz/+ee1ZANz52gbLu5P6dCJ44oscIoIDuHlqf6tDUUr5mAudjcYAb62zdlA6n00Ee0pq+HBzAddN6kt0WKDV4SilfIyI8PT3HVcFv3hjk6VXBT6bCJ78Yg/BAX78oE1dnVJKeVLbq4JPthZaFofbEoGIpInIFyKyXUS2ishPnOX3i8ghEdngfFzirhg6UlTVwDsbDnHN+HTiIoI9vXmllALAz0947GpHd9JbX1pnXRxu/O4W4GfGmCHAROA2ERnqfO8xY8xo5+NDN8bQrpdXHsBmDNdPyvD0ppVS6hiXjujTurzxYIUlMbgtERhjCowx65zL1cB2wPKR3Jptdl5ZdYBzByaQER9udThKKR8XFODH/zirqP/3va2WxOCRNgIRyQDGAEfuqb5dRDaJyAIR6dXBZ+aJyBoRWVNSUuKyWBZtLaS4upHrJvV12XcqpVRXzDu3HwDrDlRQ1dDs8e27PRGISATwJnCnMaYKeAroD4wGCoBH2vucMWa+MSbbGJOdkOC6qSLfWJtHn+gQzh3Y22XfqZRSXdE7MoSkKMdgl+9u8PwNZm5NBCISiCMJvGSMeQvAGFNkjLEZY+zAM8B4d8bQVmlNI1/tLuXy0Sn465hCSikvcv/ljibUX7+9xePbdmevIQGeBbYbYx5tU57cZrXZgMf+1R9sKsBmN1wxps+pV1ZKKQ86f/DRWooDZZ6d6N6dVwTnANcCFxzXVfTPIrJZRDYB5wN3uTGGY7y7MZ/BSZEMTory1CaVUqpTggP8mTwgHoD3Nnm2eijAXV9sjFkGtFf/4vHuogBlNY2sO3CYn0zLsmLzSil1Sj+YksmynFIeXrST284f4LHt+sydxUt3lWAMXDBYG4mVUt5pUr+41uXSmkaPbddnEsHnO4qJjwhmeB+dc0Ap5Z1CAv1Jjnb0Hlq2u9Rj2/WJRGC3G5bllHLeoASdgUwp5dV+dJ5jNOSXVx3w2DZ9IhHklNRQUdfMxDaXXUop5Y2OHKdW7fXcNJY+kQiOzAua3bfdm5iVUsprDOgd0bqcd9gz3Uh9IhGs3XeY+Igg+saFWR2KUkqdlIi03mW8Ka/SI9v0iUSwIa+CMem9cNzjppRS3u2qcWkArMwt88j2enwiaGi2sa+0liFJkVaHopRSnTIq1dG78a11hzyyvR6fCPaU1GA3MFATgVKqmxjkPF5VN7Z4ZHs9PhHsKqoGYFCiJgKlVPeQEhPaulzjgWTQ4xPBvtI6RKBvnE5Co5TqHtq2Zx4sd3/PoR6fCAorG4iPCCYooMf/U5VSPUhilGM+9QOaCLouv7K+9ZZtpZTqLsZnOm4syztc7/Zt9fhEUFjZoIlAKdXtDEl2tGt64qayHp8IKuqbiQ0PsjoMpZQ6Lb3CHMctT0xSY1kiEJGZIrJTRHJE5B53bae2sYXwILdNu6CUUm4RGeI4bvXYNgIR8Qf+AVwMDAWuEZGhrt6OzW6oa7IREaKJQLmYSPsPpVwkIthx3CqobHD7tqy6IhgP5Bhjco0xTcCrwCxXb6SuydH/Vq8IlEud7ICvyUC5SLgzEfTk+whSgINtXuc5y1zKbhzP+n9TKdXdePKwZVUiaO/faI5ZQWSeiKwRkTUlJSVntBGdg0YppU7NqkSQB6S1eZ0K5LddwRgz3xiTbYzJTkhIOKONHLk7z27MKdZUSinfZVUiWA1kiUimiAQBc4F3Xb2RAOclQbNNE4FSqntparEDEBLo/sO0Ja2oxpgWEbkdWAT4AwuMMVtdvZ3gAD8C/YXqBs+M4Kd8hDEdNzzp1adykaqGZgCSo0NPsWbXWdadxhjzIfChO7chIkSHBrbuUKVcRg/4ys2qnCewvSOD3b6tHn9ncVRIIJX1mgiUUt1LlfO4leGBkZN7fCKIDQ+itLrR6jCUUuq07C6qASArMeIUa3Zdj08Eqb1COVTh/tH7lFLKldYdOAxAkgcGzfSBRBBGQWUDLTa71aEopVSn7S52XBEkRmki6LLUXqHY7IbCKveP16GUUq6WGa9tBF12ZCfmOLOrUkp5u/LaptblOA8Mo9/jE8Hg5CgAthdUWxyJUkp1zo7CqtZl8cBgaT0+EUSHBpISE8r2gqpTr6yUUl5gW77jeDVjaKJHttfjEwHAkOQotuZXWh2GUkp1ynsbHUOvaSJwobF9Y9hTUntMvZtSSnmrjXmOE9ex6b08sj2fSAQTMmMBWLW33OJIlFLq5ArbzEjWzwM9hsBHEsGIlBhCAv1YubfM6lCUUuqkVu1znLAmRAbj56FJVXwiEQQF+HFW314sz9FEoJTybi+u2A/Aref299g2fSIRAJw/qDc7i6o5WF5ndShKKdUuu920VmFPH+KZhmLwoURwpPX90+1FFkeilFLt25hX0bqcHhfmse36TCLoGxdOVu8ITQRKKa/15ro8AK7KTvXodt2SCETkYRHZISKbRGShiMQ4yzNEpF5ENjgfT7tj+x25cFgiK3LLKdFhqZVSXsYYw4srDgBw7cQMj27bXVcEi4HhxpiRwC7g3jbv7THGjHY+bnHT9ts1e0wKNrvhnQ2HPLlZpZQ6pSPDTgMMT4ny6LbdkgiMMZ8YY45MFLwC8Ox1TgcG9I5kVGo0b63TRKCU8i5PLdkDwE3nZHpkfKG2PNFGcBPwUZvXmSKyXkSWisgUD2z/GHPGprKtoKp1LA+llLJaQ7ONT7cXA3DD2Rke3/4ZJwIR+VREtrTzmNVmnfuAFuAlZ1EBkG6MGQP8FHhZRNq9BhKReSKyRkTWlJSUnGmYJ7h8VB+CA/x4ceV+l32nUkp1xZFaiuAAP4/2FjrijBOBMWa6MWZ4O493AETkeuAy4HvGGOP8TKMxpsy5vBbYAwzs4PvnG2OyjTHZCQkJZxrmCXqFB3HF6BTeWpdHRZ2OPaSUspYxhl8t3AzA368ZY0kM7uo1NBO4G7jcGFPXpjxBRPydy/2ALCDXHTGczA3nZNDQbOe11Qc9vWmllDrG6n1HG4kvGNzbkhjc1UbwBBAJLD6um+hUYJOIbATeAG4xxnh8JLghyVFM7BfL88v30dSicxkrpaxzy4trAfjFRYMI8Lfm1q4Ad3ypMWZAB+VvAm+6Y5un65Zz+3PDc6t5c10e14xPtzocpZQP2nKosnV4/BvPybAsDp+5s/h45w5MYFRaDE98nqNXBUopS9z28joAfjglk7Agt5yXd4rPJgIR4c7pWRyqqG+9rVsppTxle0EV+8scTai3n59laSw+mwgAzhuYwOi0GB7/bDf1TTarw1FK+ZA5Ty4HHNXU0WGBlsbi04lARLj34sEUVDbwzFce77yklPJRS3eVUN/sOPm8c7q1VwPg44kAYEK/OC4ensRTS/ZQVNVw6g8opVQX2OyG6xesAuD3VwwnJNDf4og0EQBw78VDsNkNf/54p9WhKKV6uH+1qX34rpf0WNREgGMCiJsmZ/LmujxW5Op0lkop9yiuauChj3YAsPBHZ+PvoTmJT0UTgdNPpmWRHhvGr97aTEOzNhwrpVxv7jMrADi7fxxj0ntZHM1RmgicQoP8eXD2cHJLa3ni8xyrw1FK9TCfbC0kt6QWgCe/N9biaI6liaCNKVkJzBmbwtNL97DlUKXV4SileoiKuibm/ccxlMTj14whJizI4oiOpYngOL+5dChxEUH85NX1em+BUqrLjDHM+sfXAAzoHcHlo/pYHNGJNBEcp1d4EI9cOZo9JbU8+OE2q8NRSnVzL6860HoH8Zu3nG1xNO3TRNCOyVnxzJvajxdXHGDxtiKrw1FKdVMHyuq4b+EWAJ67cZzldxB3RBNBB3524UCGJkfxizc2crC87tQfUEqpNuqbbEx9+AsAZo9J4fxB1sw10BmaCDoQHODPk98bi81uuPk/a7W9QCnVacYYrvrnN62vH/7OSAujOTVNBCeRER/O3+aOZnthFfct3Ixzxk2llDqpJ5fsYbOz5+GqX02zbMKZznJbdCJyv4gccs5QtkFELmnz3r0ikiMiO0XkInfF4AoXDE7kzmkDeWv9If69fJ/V4SilvNyy3aU8vMgxXM1r8ybSOyrE4ohOzd0zITxmjPlL2wIRGQrMBYYBfYBPRWSgMcZr617uuGAAW/Ir+f3720jrFcb0oYlWh6SU8kJ7Smr4/rMrAfjtZUOZ0C/O4og6x4rrlVnAq8aYRmPMXiAHGG9BHJ3m5yf8be5ohqdEc8cr69mUV2F1SEopL1Ne28S0R5YCcOnIZG6anGlxRJ3n7kRwu4hsEpEFInJkYI0U4GCbdfKcZccQkXkiskZE1pSUlLg5zFMLCwrgX9dnExsexE3/XqM9iZRSreqbbIz9/WIAekcG8/e5YyyO6PR0KRGIyKcisqWdxyzgKaA/MBooAB458rF2vuqEVlhjzHxjTLYxJjshIaErYbpM78gQnr9pHE0tNq5fsIqS6karQ1JKWaypxc6Q337c+vrLX56Pn5eMKtpZXUoExpjpxpjh7TzeMcYUGWNsxhg78AxHq3/ygLQ2X5MK5HclDk/zl5CeAAAUZklEQVQa0DuSBTeMo6CygWufXcnh2iarQ1JKWaTFZmfqn79ofb3p/gu9YqKZ0+XOXkPJbV7OBrY4l98F5opIsIhkAlnAKnfF4Q7ZGbH86/pscktruW7BKqoamq0OSSnlYXa74ZLHv6LQObPh6vumExXinXcOn4o72wj+LCKbRWQTcD5wF4AxZivwOrAN+Bi4zZt7DHXknAHxPP39sWwvqOKm51ZT09hidUhKKQ+xOZPArqIaAJbdfT4JkcEWR3XmpDvcJJWdnW3WrFljdRjt+nBzAXe8sp7hKdE8f+M4rxteVinlWi02OzMe+5K9pY65Bb6+5wJSYkItjqp9IrLWGJN9qvW8+3a3buCSEck89b2xbM+vYu78FdqArFQP1myzc/YfP29NAsu9OAmcDk0ELnDhsCSevSGb/WV1XP3PbyiorLc6JKWUi9U32ci67yOKnSd739x7AX16QBIATQQuMyUrgRd+MJ6S6kbmPLmcHYVVVoeklHKR8tqmY7qIrvvNDJKje0YSAE0ELjUuI5ZXb56I3RiufOobvs4ptTokpVQX7S+rbb1ZDGDbAxcRG96z2gI1EbjYsD7RLPzROfSJCeX6Bat4Y22e1SEppc7Q6n3lnPvwEgDSY8PY/eDFhAW5e4g2z9NE4AZ9YkL5762TmNgvjp//dyMPL9qBze79vbOUUkc982UuVz7tmFNg1ug+LP3FeQR6+XDSZ6rnpTYvERUSyIIbxvHbd7bwjy/2sDW/ir/NHUN0aPe84UQpX9Fis3P9c6v4OqcMgN9fMZxrJ/a1OCr36pnpzUsEBfjx0JwRPDh7OF/nlDLriWXsLKy2OiylVAcq6poYcN9HrUlg4Y/O7vFJADQRuJ2I8L0JfXl13kRqm2zMfvJr3tvYbYZWUspnrMgtY/QDRxuFV/1qGmPSe53kEz2HJgIPOatvLO/fMZnBSZHc8cp67n1rk86DrJQXsNsNv357M3PnrwBgUr84ch68uFvMLOYq2kbgQYlRIbx28yQeXbyLp5fuYfW+wzzx3TEMToqyOjSlfNLh2ibGtOka+tjVo5g9JtXCiKyhVwQeFujvx90zB/PCTeOpqGvm8ie+5j/f7KM7jPmkVE/y3sb8Y5LAV7883yeTAGgisMyUrAQ++skUJvaL4zfvbOW6Bas4VKFDUyjlbrWNLVzwlyXc8cp6AC4alkjOgxeTFhtmcWTW0dFHLWa3G15auZ+HPtqBnwi/uWwIV2WnIdK9ZjhSqjtYsrOYG55b3fr6v7dMYlxGrIURuVdnRx/VNgKL+fkJ107K4NyBvfnlmxu5+83NfLi5kIfmjOgxA1opZbXqhmZm/eNrcksco4YOSY7irVvPJjSo+80m5g5uuSIQkdeAQc6XMUCFMWa0iGQA24GdzvdWGGNuOdX39eQrgrbsdsN/Vuznjx/twE/grhkDueHsDAJ66N2MSrmbMYZXVx/k3rc2t5a9cNN4pg70jnnQ3c3SKwJjzNVtAnkEqGzz9h5jzGh3bLe78/MTrj87gwsG9+a372zhDx9s5421eTw4ewRn9fWN/sxKuUre4Tom/+nofMLnDkxg/nVnERygVwHHc2vVkDgquq8CLnDndnqatNgwFtwwjkVbC/nf97bx7aeWc834NH550WB69bBRD5VytYZmG3e8sp7F24payz748WSG9Ym2MCrv5u42gilAkTFmd5uyTBFZD1QBvzbGfOXmGLolEWHm8GSmZCXw1093seDrfXywqYAfT8vi2kl99axGqeMYY3jhm/387t2trWW/uGgQPzqvv3a+OIUzbiMQkU+BpHbeus8Y845znaeAHGPMI87XwUCEMaZMRM4C3gaGGWNOmMVFROYB8wDS09PP2r9//xnF2VPsKqrm/z7czpKdJfSNC+OemYOZOTxJf+BK4Rgu+shIoQCjUqN56YcTiQj27f4wnW0jcFv3UREJAA4BZxlj2h2UX0SWAD83xpy0JdhXGos7Y+muEh78YBu7imoYl9GLey8ZwlgfGQ9FqePlFNcw/dGlx5R9+tOpDOgdaVFE3sUbuo9OB3a0TQIikgCUG2NsItIPyAJy3RhDj3PuwATO6T+F19fk8ejincx5cjkXDO7NT2cMZHiK1oEq31BQWc+ljy+jvLaptWzBDdlcMDjRwqi6L3cmgrnAK8eVTQUeEJEWwAbcYowpd2MMPVKAvx/fnZDOrNF9+Pfyfcz/MpfL/r6MmcOSuGvGQAYl6dmQ6plKqhuZ89TXHCw/ehf+/80ewTXj9SbMrtA7i3uAqoZmnv1qL88u20ttUwuXjkjm1vP6ay8J1WPkHa5j9pPLKalubC27a/pAbr9gAP5+mgA6YnkbgStpIuicw7VNzP8qlxeW76O2yca5AxO49bz+TMiM1bMl1S3tKqrm0se/otl29Dj142lZ/GRaliaATtBE4MMq65p5ceV+nvt6L6U1TYxOi+HW8/ozY0gifvqfR3k5YwxLdpZw479XH1N+3yVD+MHkTP0NnwZNBIqGZhv/XZvHM1/mcqC8jr5xYVw7sS9XZqfp3MnK6zS12HlqyR4e+3TXMeV/uXIUc8akaAI4A5oIVKsWm52PthTy/PJ9rNl/mNBAf2aPTeH6SRnasKwsV1zVwLz/rGXDwYpjyt+8dRJn9e25I4N6giYC1a4thyp54Zt9vLMhn8YWO5P6xTF3fBoXDUsiJFDvVlaeYbMb3lyXxy/f2HRM+cjUaJ65LptEH5om0p00EaiTOlzbxGtrDvLiiv3kHa4nKiSAWaNTuCo7jeEpUdq4rNziYHkdd7yy/oSz/19cNIh5U/sRqCPtupQmAtUpdrthRW4Zr685yEdbCmlssTM4KZKrstO4bFQyvSP1zEx1TXVDM4u2FvHOhkN8tbu0tTwtNpTnbhindwG7kSYCddoq65t5b2M+/11zkI15lfgJTOofx7dG9mHm8CRiwnTkU9U5jS02luws4d0N+Xy6vYjGFjvpsWEkRAZz8fAkrpuUQVCAnv27myYC1SW7i6p5b2M+720qYG9pLQF+wtSBCXxrVDIXDE7UXkfqBA3NNr7aXcqirYV8srWQqoYW4iOCuGxkHy4f3YcxaTFa5ehhmgiUSxhj2Jpf5UgKG/PJr2wgwE+Y0C+WGUMSmT40kdRevjvpt6+rrG/m8x1FfLK1iKW7SqhrshEZEsCMIYnMGpPCOf3jdIY9C2kiUC5ntxvWH6xg8bYiPt1eRE5xDeCY/3XGkN5cMCSRESnResdnD2aMYVdRDV/uKmHJrmJW5pbTYjf0jgzmwmGJXDQsiYn94rTR10toIlBut7e0lk+3FbF4WxFr9pdjNxAdGsjZ/eOYnBXP1KwE0mL1aqG7q6xrZllOKUt3FfPlrlIKqxoAGJgYwfmDe3PRsCRGp8boDV9eSBOB8qjy2iaW5ZTy1a4SluWUUlDpOFj0jQtj8oB4JvSLY1xGL5KjQy2OVJ3K4domVu8rZ9XeclbuLWdrfiV2A1EhAUzOiufcgQlMHZigf8tuQBOBsowxhj0ltSzb7UgKK3LLqWlsASAlJpTxmbFkZ/RiXEYsAxIi9EzSQsYY8g7Xs+FgRevBf0dhNQBBAX6MSYthYr84pg6MZ1RqjNb3dzOaCJTXaLHZ2VFYzep95c6DzWFKaxzDCUeFBDAiNZrhKdGMTIlhREo0abGh2rvETQ7XNrExr4KNByudzxWUOSd3CQvy56y+vZiQGcv4zDhGpUXr3NjdnCYC5bWMMRwor2PV3nLWHahgy6FKdhRWtQ41HB0ayIiUaIb1iSIrMZKBiREM6B1BWJBvzz97OpptdvaW1rKjsJqdhVXsKKhmR2E1hyocE7qIwICECEalxTgeqdEMSY7SRt4exiOJQESuBO4HhgDj2849LCL3Aj/AMRPZj40xi5zlM4G/Af7Av4wxfzzVdjQR9HyNLTZ2Fdaw6ZAjMWzKq2R3UQ1NNjvgOHCl9gplUGIkWYmR9IsPp29cOH3jwugdGeyzVxAVdU3sLa1lX1kte0tq2VtWR05xDXuKj+47fz+hf0I4g5KiGJocxei0GEakRvv8xO6+wFNzFm8B5gD/PG7jQ3FMVTkM6AN8KiIDnW//A5gB5AGrReRdY8y2LsahurngAH9GpEYzIvXorGotNjv7yurYXVTNrqIadhVXs7uomiU7S2ixHz2BCQn0Iz02jPTYcNJjw0jtFUpiVAhJ0cEkRoXQOzKkW97FarcbymqbKKisJ7+igcLKegoqG8ivbCDvcB17S2upqGtuXV/E0QbTLyGCqVnxDEqKZHBSFP17h2sVjzqpLiUCY8x2oL2zsVnAq8aYRmCviOQA453v5Rhjcp2fe9W5riYCdYIAfz8G9HZUC1084mh5s81O3uF69pfVcqC8jv1ljseB8lqW5ZTQ0Gw/4bviwoNIjAohITKYXmGBxIQFERMWSK82z1GhgYQF+RMa6E9YkD9hQQGEBPp1+WqjxWanttFGTVMLtY0t1DQ6nmsbW6ioa6astoly58Ox3Eh5TROlNU2tZ/VHBAX4kRwdQp/oUC4ZkUxmXDgZ8eFkxoeRFhumB3x1Rtx1bZgCrGjzOs9ZBnDwuPIJbopB9VCB/n5kxoeTGR9+wnvGGA7XNVNY2UBRleNRWHVkuZHSmkZyS2uoqGumuqGlU9sLDfQnNMgfPxH8/cBfBBHB38/x8BPHsMrNNkOL3U6LzdBss9NiN7TYzAkH8462ERseRFxEEPERwQxMjCQhIpjk6BCSY0LpEx1KckwIceFBPlsNptznlIlARD4Fktp56z5jzDsdfaydMgO0d33ebiOFiMwD5gGkp6efKkylAMfVaWx4ELHhQQztE3XSdZttdirrm6moa6airomqhmbqmmzUNdmob7JR33xkuYX6Zhs2u6O6xm4MNmOw2w02A3ZjCPATAvz8CPQXAvzbLvsRGuhPeHAAEcGOZ8dyAOFBAUSFBhAXHkxokJ7JK+ucMhEYY6afwffmAWltXqcC+c7ljsqP3+58YD44GovPIAalTirQ34/4iGDiI4KtDkUpS7mrBe1dYK6IBItIJpAFrAJWA1kikikiQTgalN91UwxKKaU6oUttBCIyG/g7kAB8ICIbjDEXGWO2isjrOBqBW4DbjDE252duBxbh6D66wBiztUv/AqWUUl2iN5QppVQP1dn7CLpf52qllFIupYlAKaV8nCYCpZTycZoIlFLKx2kiUEopH9cteg2JSAmwvwtfEQ+UuigcV/LGuLwxJtC4TpfGdXp6alx9jTEJp1qpWySCrhKRNZ3pQuVp3hiXN8YEGtfp0rhOj6/HpVVDSinl4zQRKKWUj/OVRDDf6gA64I1xeWNMoHGdLo3r9Ph0XD7RRqCUUqpjvnJFoJRSqgM9KhGIyJUislVE7CKSfdx794pIjojsFJGL2pTPdJbliMg9HojxNRHZ4HzsE5ENzvIMEalv897T7o7luLjuF5FDbbZ/SZv32t13HorrYRHZISKbRGShiMQ4yy3dX84YPPrbOUkcaSLyhYhsd/7+f+Is7/Bv6sHY9onIZuf21zjLYkVksYjsdj738nBMg9rskw0iUiUid1qxv0RkgYgUi8iWNmXt7h9xeNz5e9skImNdFogxpsc8gCHAIGAJkN2mfCiwEQgGMoE9OIbB9ncu9wOCnOsM9WC8jwC/dS5nAFss3Hf3Az9vp7zdfefBuC4EApzLfwL+5CX7y9LfznGxJANjncuRwC7n363dv6mHY9sHxB9X9mfgHufyPUf+phb+HQuBvlbsL2AqMLbtb7mj/QNcAnyEYwbIicBKV8XRo64IjDHbjTE723lrFvCqMabRGLMXyAHGOx85xphcY0wT8KpzXbcTx8SzVwGveGJ7XdDRvvMIY8wnxpgjkwuvwDGrnTew7LdzPGNMgTFmnXO5GtjO0TnCvdEs4Hnn8vPAFRbGMg3YY4zpyg2rZ8wY8yVQflxxR/tnFvCCcVgBxIhIsivi6FGJ4CRSgINtXuc5yzoq94QpQJExZnebskwRWS8iS0VkiofiaOt25yXngjaX61buo+PdhOOM6Agr95c37ZdWIpIBjAFWOova+5t6kgE+EZG14piHHCDRGFMAjiQG9LYgriPmcuzJmNX7CzreP277zXW7RCAin4rIlnYeJzsbk3bKzEnKPRHjNRz7AywA0o0xY4CfAi+LyMlnX3dtXE8B/YHRzlgeOfKxdr7KpV3NOrO/ROQ+HLPdveQscvv+OlXY7ZRZ2gVPRCKAN4E7jTFVdPw39aRzjDFjgYuB20RkqgUxtEsc0+VeDvzXWeQN++tk3Pab69JUlVYwxkw/g4/lAWltXqcC+c7ljsrP2KliFJEAYA5wVpvPNAKNzuW1IrIHGAi4bGq2zu47EXkGeN/58mT7ziNxicj1wGXANOOsLPXE/joFt++X0yEigTiSwEvGmLcAjDFFbd5v+zf1GGNMvvO5WEQW4qhSKxKRZGNMgbNqo9jTcTldDKw7sp+8YX85dbR/3Pab63ZXBGfoXWCuiASLSCaQBawCVgNZIpLpPDuY61zX3aYDO4wxeUcKRCRBRPydy/2cMeZ6IJYj229b1zgbONKLoaN956m4ZgJ3A5cbY+ralFu6v7Dut3MCZ3vTs8B2Y8yjbco7+pt6Kq5wEYk8soyj4X8Ljv10vXO164F3PBlXG8dclVu9v9roaP+8C1zn7D00Eag8UoXUZZ5sIfdAC/xsHFmzESgCFrV57z4cvTx2Ahe3Kb8ERy+LPcB9Horz38Atx5V9G9iKo/fJOuBbHt53/wE2A5ucP7jkU+07D8WVg6NedIPz8bQ37C+rfjsdxDEZRxXBpjb76ZKT/U09FFc/599no/NvdZ+zPA74DNjtfI61YJ+FAWVAdJsyj+8vHImoAGh2Hrt+0NH+wVE19A/n720zbXpGdvWhdxYrpZSP85WqIaWUUh3QRKCUUj5OE4FSSvk4TQRKKeXjNBEopZSP00SglFI+ThOBUkr5OE0ESinl4/4fIS2jo3eRsMcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8ZGd95/vPT7tqUUmlKu1SS2pJvbpXudu7jRcwHmOHsARzGQzO4PBKPBNufO8NgRsg4c5cQgjcZEgg5sKAM2CWGGPjQLwvgFe1e2/1pq1brX2vKu2qZ/44R6VSW90tt1RVUun3fr3qVaWnTtV5+lT19zz1nOecR4wxKKWUSl4pia6AUkqp2NKgV0qpJKdBr5RSSU6DXimlkpwGvVJKJTkNeqWUSnIa9EopleQ06JVSKslp0CulVJJLS3QFAHw+n6msrEx0NZRSalXZt29fnzHGf6nlVkTQV1ZW0tDQkOhqKKXUqiIibYtZTrtulFIqyWnQK6VUktOgV0qpJKdBr5RSSU6DXimlktwlg15EykXkBRFpFJGjIvKndrlXRJ4RkVP2fZ5dLiLyDyJyWkQOiciuWP8jlFJKXdhiWvTTwIPGmE3AVcCfiMhm4LPAc8aYWuA5+2+A9wK19u1+4FvLXmullFKLdslx9MaYTqDTfhwQkUagFLgbuMle7AfAi8Cf2+UPG2uOwtdEJFdEiu33UUolsR+/cYaxqRk2FLqpK3Ljc2UmukqKd3jClIhUAjuB14HC2fA2xnSKSIG9WClwNupl7XbZvKAXkfuxWvxUVFRcRtWVUivNz986xxutA5G/850Z1Ba6IsG/scjNhqIcXJkr4lzNNWPRW1tEXMCjwGeMMSMicsFFFyh72wzkxpiHgIcA6uvrdYZypZLAT/7oKvqCk5zoCnCiO8BJ+/5f97UTmpyJLFfhdbCp2M2m4hw2FuWwuTiHsrxsUlIumCtqCRYV9CKSjhXyPzTG/Nwu7p7tkhGRYqDHLm8HyqNeXgZ0LFeFlVIrl4jgd2fid2dyXa0vUh4OG84NjXG8K8DxzhEau0Y43hng6WPdGLuZ58pMY4Pd6t9UnMOm4hw2FLm19b8MLrkFxWq6fxdoNMZ8PeqpJ4B7ga/Y949HlT8gIj8G9gLD2j+v1NqWkiKUex2Uex3ctrkwUj46Oc3J7iCNnSPWDqAzwBMHO/jh62ciy6zLd7C5OIetpR6usG95zoxE/DNWLTHm4r0mInId8BvgMBC2iz+H1U//U6ACOAN8yBgzYO8YvgncDowCnzTGXPSKZfX19UYvaqaUAjDGbv13Bmi0W/9Hzo1wZmA0skxpbrYV+mWeyA7AuwbDX0T2GWPqL7ncpYI+HjTolVKXMjw6xZGOYQ6fG+aIfWvtnwv/Ek8WW0s9bCvzsLMij21lHtxZ6QmscewtNui180sptSp4HOlcW+Pj2pq5vv/hsSmOdlihf/jcCEfODfP0sW4ARKCuwM2udbnsLM9jZ0Uu6/2uNXnAV1v0SqmkMjw6xYH2IfafGWT/Get+ZHwaAHdmGjsqctlZnkt9pZfd6/JwruKDvdqiV0qtSR5HOjfW+bmxzpp4KRw2NPeFrOA/O8T+M0N884XThA2kpghbSz1cVeVlb7WX3eu8eLKTr7tHW/RKqTUnODHNvrZB3mjp5/XmAQ62DzE1YxCBzcU57K3K5/paH3urvTgyVm57WA/GKqXUIo1PzfDWmUHeaBng9eYB3jozyMR0mIzUFOor87ihzs/1tT42FeWsqD5+DXqllLpM41MzvNk6wMsne/nNqT6OdwUA8LkyuXmjn/dsKeLaGh9Z6akJraf20Sul1GXKSk/l+lo/19da/fzdI+P85lQfL57o4deHu/hpQzuOjFRurLNC/5ZNBSt6KKe26JVS6h2YnA7zWnM/Tx/r4umj3fQEJshMS+HWzYX8/s5Sbqjzk54anzmdtOtGKaViLBw2HGgf4vH95/jloU4GQpN4nRm8f2cp//GqdVT6nDFdvwa9UkrF0dRMmJdP9vLoW+08fbSbGWO4qc7PfddVcV2Nj4tc8feyadArpVSCdI+M86PXz/CjN87QG5hgV0Uu//ttdcse+Br0SimVYBPTM/ysoZ1/fOE0ncPj3LKxgC/dtYVyr2NZ3n+xQR+fIwZKKbUGZaal8rGr1vHi/3kTn79jE68293PbN17i0X3tca2HBr1SSsVYZloqn7qhmmf/7EZ2lOfy4M8O8v/+qpF49aho0CulVJyU5GbzP/9wLx+7qoJ/frmZf365OS7rvWTQi8j3RKRHRI5Elf1ERA7Yt1YROWCXV4rIWNRz345l5ZVSarVJS03hy3dv5T9cUczXnjpBW38o5utcTIv++1izRUUYY/7AGLPDGLMDay7Zn0c93TT7nDHm08tXVaWUSg4iwhfet5mwMfz8rXMxX98lL4FgjHlZRCoXes6eNvDDwM3LWy2llEpuhTlZFOZkcXZw9NILL9FS++ivB7qNMaeiyqpEZL+IvCQi1y/x/ZVSKin1BMbpCUxQlpsd83Ut9aJm9wCPRP3dCVQYY/pFZDfwCxHZYowZOf+FInI/cD9ARUXFEquhlFKrhzGG/+fJRgT4vZ2lMV/fZbfoRSQN+H3gJ7NlxpgJY0y//Xgf0ATULfR6Y8xDxph6Y0y93++/3GoopdSqMhM2/NUvj/HEwQ7+9JZaqv2umK9zKS36W4HjxpjIyH8R8QMDxpgZEakGaoH4jB9SSqkVrmdknAd/dpDfnOrjk9dW8sDNNXFZ7yWDXkQeAW4CfCLSDnzRGPNd4CPM77YBuAH4axGZBmaATxtjBpa3ykoptbpMzYR5+NU2vvHMSSZnwvzNB67gD66MX5f1Ykbd3HOB8k8sUPYo1nBLpZRa8yanw/zrvna+9dJpzg6McWOdn7+6a0vML198Pp1hSimllllPYJyfNbTzP19ro3N4nO1lHv7qri28a0NBTC5XfCka9EoptQxmwoZXm/p55I0zPHW0i+mw4erqfL7ygW3cUBub69Evlga9UkpdpnDYsP/sIE8c6ODfDnfRF5zAk53OJ66p5J69FayPw4iaxdCgV0qpd2ByOkxD6wDPNvbw1NEuzg2NkZGWwi0bC3jf9hJu3lhAVnpqoqs5jwa9UkpdQn9wghdP9PL88R5ePtlLYGKajLQUrl2fz4PvruO2zYW4s9ITXc0L0qBXSqnzjE5O82brIK809fFqUz+Hzw1jDBS4M7lzezE3byzk2pp8HBmrI0JXRy2VUiqGJqfDHDg7xCtNfbxyup/9ZweZmjGkpwo7y/P4zC113LKpgM3FOaSkJO6g6uXSoFdKrTlDo5O8dWaQhtZBGtoGOdQ+xPhUGBHYWuLhvuuquGa9jysr81ZNq/1iVv+/QCmlLsIYQ2v/KA2tA+xrG2Rf2yCneoIApKUIW0py+Oiedeyp8nJ1dT4ex8rta79cGvRKqaTSG5jgUPsQB9uHOdw+xKH2YfpDkwDkZKWxa10ed+8oYfc6LzvKc8nOWFkjZGJBg14ptWoNj05x+NwwB9uHOGSHeufwOAAiUFvg4qYNBexal0v9Oi+1Ba5V2ce+VBr0SqlVoS84QWPnCMc6RjjaMcKh9iFa++dmZ6rMd1Bf6WV7mYdtZblsKcnBmakRBxr0SqkVZiZsaO0PcaxjhGOdI5Fw7wlMRJYp8WRxRZmHD9WXs63Mw7bS3KTsW18uGvRKqYQZnZzmeFcgEurHOkY40RVgbGoGsA6W1hS4uK7Wx+biHDYX57CpOIc8Z0aCa766aNArpWJuYnqGlr4QJ7oCnOwOcKIryMnuAGcHRzHGWiYnK43NJTl8ZE+5FeolOdQUuMhMS/6DpbGmQa+UWjbTM2HaBkY52RXgZLcV5ie6A7T0hZgJW4memiJU+5xcUebhA7vK2FTsZnNJDqW52Qm9wmMyW8wMU98D7gR6jDFb7bIvAZ8Ceu3FPmeM+ZX93F8Af4g1w9R/McY8FYN6K6USKBw2nBsaiwT5bLCf7g0yOR0GrFEvFV4HdYVubt9SRF2Rm7pCF1U+p7bS42wxLfrvA98EHj6v/BvGmK9FF4jIZqwpBrcAJcCzIlJnjJlZhroqpeJsYnqG1r5RTvcEaeoNRu6be0ORfnSwDo7WFbm5rtZHXaGbDYVuagpca2KM+mqwmKkEXxaRykW+393Aj40xE0CLiJwG9gCvXnYNlVIxNzw2FQnxpqhQPzMwit3jAkBpbjY1BS72VuVTU+BiQ5GL2kI3OSv4yo1qaX30D4jIx4EG4EFjzCBQCrwWtUy7XaaUSjBjDJ3D429rnZ/uCdEXnBu6mJGaQpXPyeaSHO7aXsL6Ahfr/S6q/c6kuO7LWnS5n9q3gC8Dxr7/O+A+YKEjKWaBMkTkfuB+gIqK+M2GrlSyG52cprVvlJa+EC19QZp6Q5FQH52c627JyUqjpsDFuzb4qbHDvKbARVleNmmpKQn8F6jldllBb4zpnn0sIt8BnrT/bAfKoxYtAzou8B4PAQ8B1NfXL7gzUEotbGomzNmB2TCff5u9BMCsEk8W6wtcfLi+nPUFLmr8LtYXOPG7MnWUyxpxWUEvIsXGmE77z/cDR+zHTwA/EpGvYx2MrQXeWHItlVqDwmFD18j4gmF+ZmA0MlwRwJOdTpXPydXV+VT5nFT5nVT5nFTmO/UyAGpRwysfAW4CfCLSDnwRuElEdmB1y7QCfwRgjDkqIj8FjgHTwJ/oiBulLswYw+DoFC191kiW1n4ryGcfj0+FI8tmpadQme9kU7GbO64oospnDVWs9jn1TFF1UWJM4ntN6uvrTUNDQ6KroVTMhCamaemzg7zXDnO7dT48NhVZLi1FqPA6qPRZLfLZIK/yOyl0Z63JKy+qCxORfcaY+kstp7/plFomk9NhzgyM0toXHeRBWvpCdI9MzFu2xJNFld/JnduKrTD3O6nyWQdC0/VAqFpmGvRKvQPhsKFjeMxqnUe1ylv6Qpw9b8y515lBlc/JdTV+O8jn+s31RCIVTxr0Sp3HGEN/aHL+AVC7u6W1P8TE9Fy/uSMjlSqfk62lHu7aXhIJ8yqfk1yH9purlUGDXq1ZgfEpa7x5pN88GOlyCYxPR5ZLT7X6zat8Lm6o880dBPU7KXDrEEW18mnQq6Q2MT3Dmf7RSBdLdHdLb9REFiJQ4smm2u/k93aURoYoVvuclObqCURqddOgV6veTNjQMTRmBXjvXKu8tT/EucGxef3mPpfVb35TnT8S5FU+F+vyHWSla7+5Sk4a9GpVMMbQG5iY128+2zI/0z/K5Mxcv7krM40qn5Od5Xn8/s6yuYOgPieebL34llp7NOjVijI8NrXAiJYgLb0hQlHXaclITWFdvoNqn5NbNhVEWuaVPoee2q/UeTToVdxNTodp6w/R1Bui2Q7x2VZ6f2gyslyKQFmedfJQ/TrvvBEtJbnZpOrJQ0otiga9ipmB0CTNvfY1zntDNPUEaV7gOi0F7kyqfE5u21w4dzao30m516EzESm1DDTo1ZJMz4Q5OzgWmayiuTdkB3uQwdG5U/sz0lKo9lnXablzW3Hk+uZVPidunbRCqZjSoFeLMjw2ZbfOQ/Na6W39IaZm5lrnPlcm1X4nt28tZr3faU1a4XNRmqddLUoliga9mmcgNMmp7gCneoJz9z3BeWPO01KEdfkO1vtd3La5kGrfXKB7HNo6V2ql0aBfo/qCE5zqDnKqJzDvPvpgqCvTmoHoxrq5GYjW233neuEtpVYPDfokZoyhNzjB6W6rVX7SbqGf7gkyEBXo7sw0agpd3LqpkNpCa7Ln2gIXxZ4sHaaoVBLQoE8SY5MznOwOcLxrhMZO6/5EV2DeAVF3Vhp1hW7es6WQmgIrzOsK3RTm6LhzpZLZYmaY+h5wJ9BjjNlql/0t8D5gEmgCPmmMGRKRSqAROGG//DVjzKdjUO81Kxw2tA+O0dg1wnE70I93BWjtDzE7h0x2eiobitzcvrWI2gI3dYVuagtdegEupdaoxbTovw98E3g4quwZ4C+MMdMi8jfAXwB/bj/XZIzZsay1XKNGJ6dp7AxwtGM40ko/2RWInCEqAuu8DjYW5XD3jhI2FuWwschNhdehMxEppSIuGfTGmJftlnp02dNRf74GfHB5q7X2jIxPcfTcCEc7hjlybpijHSM09QYjF+TyZKezqdjNh+rL2VjkZmNxDnWFLhwZ2vumlLq45UiJ+4CfRP1dJSL7gRHg/zbG/GahF4nI/cD9ABUVFctQjdVjIDTJkXPDHOkY5ui5EY50DNPWPxp5vigni62lOdxxRTFbSz1sKcnRA6NKqcu2pKAXkc8D08AP7aJOoMIY0y8iu4FfiMgWY8zI+a81xjwEPATW5OBLqcdKNjE9w9GOEQ6cGeLAWet2ZmAu1Mu92Wwt8fDh+nK2lOSwpcSD352ZwBorpZLNZQe9iNyLdZD2FmOsw4DGmAlgwn68T0SagDqgYRnquuIZY2jrH40E+v6zQzR2jEQuoVuYk8mO8lzu2VPB9jIPW0o8eoKRUirmLivoReR2rIOvNxpjRqPK/cCAMWZGRKqBWqB5WWq6As2EDY2dI7zRMsAbLQM0tA3QF7TGp2enp3JFmYdPXlvJjvJcdlTkUuzJTnCNlVJr0WKGVz4C3AT4RKQd+CLWKJtM4Bm733h2GOUNwF+LyDQwA3zaGDMQo7rH3cT0DIfahyPBvq9tkOCENbdoWV42N9T62V2Zx87yPOoKXTr9nFJqRVjMqJt7Fij+7gWWfRR4dKmVWimMMZzoDvDbU3385lQfb7QMMDZlDW2sK3Rx944S9lR5ubLSS0muttaVUiuTjs07z0BokhdP9PCbU3389nRf5GJe1X4nH64v45oaH3sqveQ5MxJcU6WUWhwNeqClL8Szx7p55lg3DW0DhA14nRlcW+Pj+hof19X6tMWulFq11mzQN/UGefxAB/92qIOm3hAAm4pzeOBdNdy6uZCtJR49u1QplRTWVND3BMZ54kAHjx/o4PC5YUTg6up8Pn51JbdsKqAsz5HoKiql1LJL+qA3xvB6ywD/8lobTx3pYjps2Fbm4S/v3Mz7thVTkJOV6CoqpVRMJW3Qz4QNTx7q4J9eaOJEdwBPdjr3XlPJPXsqqClwJbp6SikVN0kX9MYYfn2ki689dYLmvhC1BS6++sFt3LW9hKz01ERXTyml4i6pgr61L8TnHjvMK039bCh080//2y5u31KkB1WVUmta0gT9k4c6+OyjhxGBv757Cx/dU6FnpiqlFEkS9P/yWht/+Ysj7F6Xxz/cs5NSHfOulFIRqz7oXzndxxceP8ItGwv4p4/tIjNN++GVUiraqu7bMMbw5X9rZJ3XwX//6E4NeaWUWsCqDvpzQ2M0do5w7zWVOqWeUkpdwKoO+lF7kmyvXmBMKaUuaFUH/bp8B+7MNJ462pXoqiil1Iq1qKAXke+JSI+IHIkq84rIMyJyyr7Ps8tFRP5BRE6LyCER2RWrymempfLJayv51eEuHj9wLlarUUqpVW2xLfrvA7efV/ZZ4DljTC3wnP03wHuxphCsBe4HvrX0al7YH7+rhr1VXh786UF+9PqZWK5KKaVWpUUFvTHmZeD8KQHvBn5gP/4B8HtR5Q8by2tArogUL0dlF5KVnsp37q3nmhofn3vsMJ/58X76gxOxWp1SSq06S+mjLzTGdALY9wV2eSlwNmq5drssZnKy0vkfn7iSz9xay78d7uSWr7/E93/Xwrg97Z9SSq1lsTgYu9CFZczbFhK5X0QaRKSht7d3yStNTRE+c2sdv/ov17OxyM2XfnmMm/72RR5+tZXRyeklv79SSq1WSwn67tkuGfu+xy5vB8qjlisDOs5/sTHmIWNMvTGm3u/3L6Ea89UWunnkU1fxw/+0l7K8bL7w+FGu+m/P8eUnj9HaF1q29Sil1GqxlLOMngDuBb5i3z8eVf6AiPwY2AsMz3bxxIuIcG2Nj2vW59PQNsjDr7bxg1da+e5vW7i2Jp+7d5Ry+9YicrLS41ktpZRKCDHmbb0qb19I5BHgJsAHdANfBH4B/BSoAM4AHzLGDIiIAN/EGqUzCnzSGNNwsfevr683DQ0XXWTJekbG+fGbZ3n0rXba+kfJSEvhlo0F3LW9hBs3+PXMWqXUqiMi+4wx9ZdcbjFBH2vxCPpZxhgOtg/zi/3nePJQB33BSTLTUriuxsetmwu5ZWOBTi+olFoVNOgXYXomzOstAzzb2M0zx7ppHxwDYHuZh5s3FnJdrY/tZR69rr1SakXSoH+HjDGc6A7w7LFunmns4eDZIQDcmWlctT6f62p8XFvjY73fidU7pZRSiaVBv0QDoUleaerjd6f7+O3pPs4OWK39Yk8W16z3sacqj/pKL9U+DX6lVGJo0C+zM/2j/K7JCv1Xm/oZCE0CkO/MYPe6PPZUeamv9LKlJId07epRSsXBYoNeh5osUkW+g4r8Cu7ZU4Exhua+EA2tA7zRMkhD2wBPH+sGICs9hZ3ledRX5rG9LJdt5R4K3HpwVymVONqiXyY9I+M0tA3yZusADa2DHOscYSZsbdsSTxbby3OtW1kuV5R5cGXqPlYptTTaoo+zgpws7riimDuusK7fNjo5zdGOEQ6eHeJg+zAHzw7x6yPWdfNFoMbvsoPfw7ayXDYUuclK16kQlVLLT4M+RhwZaVxZ6eXKSm+kbCA0yaH2IQ6eHeZg+xAvHO/hX/e1A9a1etb7nWwp8bClJIfNxTlsKfHgcejZu0qppdGumwQyxtA+OMaRc8Mc6xzhaMcIRzuG6R6Zu8xyaW62FfwlOZGdQLEnS0f6KKW062Y1EBHKvQ7KvQ7ee8XcJfv7ghMc67CC39oBDPNMYzez++Q8RzqbS3LYUJjDxiI3G4rc1Ba69DIOSqkFaTKsQD5XJjfU+bmhbu6qnqOT0zR2BjjWYbX+j3WM8MgbZxizr7kvAhVeB3WFbjYWuSP3lT6nDvdUao3ToF8lHBlp7F6Xx+51eZGycNhwZmCUE90BTnQFIvfPH++JjPjJSE2h2u9kg93y31Bo3ZfmZmv3j1JrhAb9KpaSIlT6nFT6nLxnS1GkfGJ6hqaeECe7AxzvCnCyO0BD6yCPH5ibFsCVmUZtoYsav8u6L3BR43dTlpdNSoruAJRKJhr0SSgzLZXN9gHcaCPjU5zqDnCiK8iJrhFOdgd58WQvP7NH/livTaHabwV/bYG9AyhwUZnvJCNNu4CUWo006NeQnKx0dq/zsnudd1758OgUp3sDnO4JRm4Hzg7y5KGOyAHg1BRhndfB+vN2AOv9Lpx68pdSK5r+D1V4HAvvAMYmZ2jqDdLUO7cDONUT5IXjPUyH54bllniy7B2A2w5/J1V+J35Xph4HUGoFuOygF5ENwE+iiqqBLwC5wKeA2Rm/P2eM+dVl11AlTHZGKltLPWwt9cwrn5oJ09Y/yukeaydwqjvA6d7gvFFAYB0HqPI5qfY7qfJZt2qfi0qfA7dO46hU3CzLCVMikgqcw5oj9pNA0BjztcW+fq2eMJVswmFDx/AYzb0hWvqsW3NfiJa+IO2DY0R/1fzuTDv4o3YCfhcVXoceC1BqkeJ9wtQtQJMxpk1/qq9dKSlCWZ6DsjzHvHMAAManZjgzMBq1EwjS0hfi2cZu+oKTc+8hUO51RP0CcFLlc1Hld1Kck6UjgpS6DMsV9B8BHon6+wER+TjQADxojBlcpvWoVSorPZW6QutErvMNj03R2heiuS9IS+/sr4AQb7QMMDo51xWUmZZyXleQK7IzyHNmxPOfo9SqsuSuGxHJADqALcaYbhEpBPoAA3wZKDbG3LfA6+4H7geoqKjY3dbWtqR6qORjjKEnMBH5FdDcG4x0CZ0ZGJ13QDjXkT7vV0C139oJVOY7yc7Qq4Kq5BS3GaZE5G7gT4wx717guUrgSWPM1ou9h/bRq3dqaiZM++AYLX3B+ccEekN0jYzPW7bEk0VV1K+A2eMCZXnZOvG7WtXi2Ud/D1HdNiJSbIzptP98P3BkGdah1DzpqSmRFvzNG+c/F5qYprXfDv/euYPCTxzoYGR8Ouo9hAqvwwr/eSODnPjdOjRUJY8lBb2IOIDbgD+KKv6qiOzA6rppPe85pWLOmZlmX9J5/rBQYwyDo1M09wYjxwFmdwQvn+plcjo89x4ZqfavAJfdFWTtBCp9TnJ0aKhaZfR69EoxNzQ0ugvoQkNDfa7MuWGh9g5gvd9JuddBZpoeD1Dxo9ejV+odiB4aen3t24eGnh0YnfcroLkvyHPHu+lrmD80tCwvamiof+78AB0aqhJJg16pS8hKT6W20E3tRYaGzp0cZv0KaGgdIHSJoaHVfut4QK5Dh4aq2NKgV2oJPNnp1iTv5bnzyi80NPR4Z4Cnj3bPGxqa50iPDAed7Qaq8rlYl+/QCePVstCgVyoGRITCnCwKc7K4en3+vOemZsKcHRiNHA9o6rV+BfzmVG9ksnjrPaw5g88/N6DK56Q0V+cNUIunQa9UnKWnWtf8r/a73vZccGLaPkt4/glij751juDE3NDQjLQUqvLnDgjPjQxy4dWzhNV5NOiVWkFcmWkLXjHUGENvcGLeeQHNvSFO9QR47ng3UzPzzxJe77cuFz07Z8B6v4tyr4NU/RWwJmnQK7UKiAgF7iwK3FnsrZ7fFTQdOUvY2gE09QZp6gny/PFeftow1xWUYZ9ktr7ASY3fxXp7J1Dtd+LI0ChIZvrpKrXKpaWmROYOftd5z1mzhwUjE8g09QRp7Azw70e6iDoeTGluNtXn/QKoKXDhc2XoGcJJQINeqSRmzR6Wx+51efPKJ6ZnaOsfpalnbgaxpt4QP3nz7LwrhuZkpUVa/rUFLuoK3dQWuijNzdYdwCqiQa/UGpSZtvBlo40xdA6PR4V/kKaeEC+dnD8iyJlhnVtQVzgb/tbjopws3QGsQBr0SqkIEaEkN5uS3Oy3nSE8NDrJqZ4gJ7sDnOwKcLI7yPPHe+YdB3Bnpdk7EGsO4dnHepG4xNKgV0otSq4jgysrvVxZOX8S+YHQJCe7A5zqDnCi29oB/PuRLh4ZPRv12nTqCtxsKHKzsdjNpuIcNha59SBwnOhWVkotideZwVXV+VwVNRrIGENfcJJT3QFOdgc40W1NIv+wuudsAAARNklEQVSL/ecIvGadDyAClflONhZZwT8b/mV52v+/3DTolVLLTkTwuzPxuzO5psYXKTfG0D44RmPnCI2dAY53jdDYOcK/H+2KXCHUnZXGpqKcSMt/dgegl4O4fBr0Sqm4ERHKvQ7KvQ7evaUoUh6amOZEd8DeAVg7gUf3tUcuDJeaItQWuNha6uEK+4SyzcU5Ok3kImnQK6USzpmZxq6KPHZVzA0DDYet1v+xzhGOdgxzqH2YF473REb/pKYINf7Z8M/hijIPm4s9Gv4LWHLQi0grEABmgGljTL2IeIGfAJVYs0x92BgzuNR1KaXWjpQUoSLfQUW+g9u3Wq3/2eGfh88Nc+TcMIfPDfPSyR4efcsK/xSBukI3Oyty2VmRx66KXKp9rjV/AbjlmBy8Fag3xvRFlX0VGDDGfEVEPgvkGWP+/ELvoTNMKaUulzGGrpFxDrdb4X+gfZgDZwYj8wPnZKWxoyKPneW57FqXx46yXDyO5JgOMtEzTN0N3GQ//gHwInDBoFdKqcslIhR7sin2ZEf6/cNhQ3NfiLfODLL/zBD7zwzy358/FbnsQ02BiysrvVxV7WVvVT5FnqwE/gtibzla9C3AINZk4P9sjHlIRIaMMblRywwaY/LOe939wP0AFRUVu9va2pZUD6WUupjgxDSHzg6x/+wQDa0DNLQOErAv/VzhdbC3ysve6nz2VnlXzRDPxbbolyPoS4wxHSJSADwD/GfgiUsFfTTtulFKxdtM2NDYOcJrzf283jLAm60DDI1OAdZF3m6o83FDrZ9ranx4sldmV0/cum6MMR32fY+IPAbsAbpFpNgY0ykixUDPUtejlFLLKTVFItf+/0/XVxMOG072BHi9eYDfne7jlwc7eeSNs6SmCDvKc7m+1seNdX62l+WuuoO7S2rRi4gTSDHGBOzHzwB/DdwC9EcdjPUaY/6vC72PtuiVUivN1EyY/WeG+M2pXl4+2cuhc8MYA353JrdtLuQ9W4q4ujqfjLSUhNUxLl03IlINPGb/mQb8yBjzX0UkH/gpUAGcAT5kjBm40Pto0CulVrrB0CQvnezl6WNdvHiil9HJGdyZabxrYwF37yjhhjo/6anxDf249dEvBw16pdRqMj41w+9O9/HU0S6ePtbN0OgU+c4M3re9hA/sKmNraU5cDuZq0CulVBxMTod56WQvj+1v59nGHianw2wr8/Dxqyu5c1txTK/Ro0GvlFJxNjw2xRMHzvHwq22c6gmS78zgvuuquPeaSlyZy3/akga9UkoliDGGV5r6+e5vW3j+eA95jnT++KYaPnFt5bL242vQK6XUCnDg7BDfeOYkL53spa7QxVc+sG3exduWYrFBn7hxQUoptQbsKM/lB/ft4Tsfryc0McOHv/0q3/1tC/FsZGvQK6VUHNy2uZBff+Z6bt5YwJefPMbfP3cqbuvWoFdKqTjJyUrn2x/bzQd3l/H/PXuKZ451x2W9GvRKKRVHKSnCf3v/FdQVuvjqvx+PSxeOBr1SSsVZRloK9+yp4FRPkO6RiZivT4NeKaUSYPYaOVMz4ZivS4NeKaXiLBw2PPbWOUpzsynNzY75+jTolVIqzv7+uVM0tA3yn2+uicslj2M1laBSSqnzTM2E+cqvj/Pd37bwgV1l/MGV5XFZrwa9UkrFwdGOYT7/2BEOnB3iE9dU8pd3bo7bdIUa9EopFUNdw+P84wun+eHrbeQ6MvjmR3dy57aSuNZBg14ppWLg7MAo336piZ81tDNjDB+7ah0P3rYBjyP+889edtCLSDnwMFAEhIGHjDF/LyJfAj4F9NqLfs4Y86ulVlQppVa66ZkwL5zo5Uevt/HiyV7SUoQP7i7nj29aT7nXkbB6LaVFPw08aIx5S0TcwD4RecZ+7hvGmK8tvXpKKbWyGWNo7Azw5KEOHtt/js7hcQrcmTzwrho+ureCYk/sh09eymUHvTGmE+i0HwdEpBEoXa6KKaXUSna6J8iThzr45cEOmnpDpKYI19b4+OL7tnDLpoK4zx97McvSRy8ilcBO4HXgWuABEfk40IDV6h9c4DX3A/cDVFRULEc1lFIqZqZmwrzZOsDzjT08f7yH5r4QIrCn0ssnr63ivVuLyHdlJrqaC1ryxCMi4gJeAv6rMebnIlII9AEG+DJQbIy572LvoROPKKVWou6RcX57qo/nT/Tw8oleAhPTZKSmcNX6fG7e4Oe9VxRTmJOVsPotduKRJbXoRSQdeBT4oTHm5wDGmO6o578DPLmUdSilVLwMhiZ5rbmfV5r6eaWpj6beEAB+dyb/YVsx79pYwHU1PpwxmP81lpYy6kaA7wKNxpivR5UX2/33AO8HjiytikopFRsDoUneahvk9ZZ+fne6n8auEYwBZ0Yqe6q8fOTKCq5en8/m4py4XKogVpayW7oW+I/AYRE5YJd9DrhHRHZgdd20An+0pBoqpdQyMMbQ1BtkX9sgDa2D7DszSLPdYs9IS2F3RR5/dmsd19Tks60sd0UdTF2qpYy6+S2w0C5Ox8wrpRIuODHN4fZh3jozyFttVrAPjU4BkOdIZ/e6PD64u4z6dV62lXnISk9NcI1jZ3V1NCml1ALGp2Y42jHC4fYhDrUPc+jcME29QWbHmqz3O3nP5iJ2r8tjd2Ue1T5n3K4zsxJo0CulVpWpmTAnugJWoLcPcbB9mJPdAWbCVqr73ZlsL/Pwvm0lbCv3sKMslzxnRoJrnVga9EqpFWt4bIrjnSMc6xyh0b4/2R1kctqalcmTnc62Mg83b6xmW1ku28tyKczJXFOt9cXQoFdKJZwxhvbBMY51jnCsYy7Y2wfHIsvkOzPYXJLDJ66p5IpSD9vLcin3ZmuoL4IGvVIqrkIT05zqCXKia4TGzgDHOkZo7BohMD4NgAhU+5zsKM/lo3sr2FScw5biHPxubalfLg16pVRMjE/NcLonyMnuACe7Z+8D81rpjoxUNhXncPeOEjYXe9hU7GZDkRtHhkbTctKtqZRaksnpMK39IU50BTjVHeCEHext/SHs46Okpwrr/S52VuTxkSvLqSt0U1fopsLrWNUnIq0WGvRKqUWZnA7T1h+iqTc4r4Xe3Bti2k70FIFKn5ONRW7et72EDYVuNhS5WJfvTKoTkFYbDXql1Dwj41M09QQ53ROkqTfE6Z4gzb1B2gZGI0MYASq8DuoKXdy6qTDSQq/2O5P6xKPVSoNeqTXIGEPXyLgV5j1BTvcGaeqxWus9gYnIcumpQmW+kw1Fbu64opiaAhfr/S6q/c5Vd2GvtUw/KaWSWHR3y2wLvanXCvfQ5ExkOXdWGjUFLm6o80fCvKbARXleNmna5bLqadArtcoZY+gJTNDSF4rcmntDC3a3lHiyWF/g4kP15awvcLHe76SmwIXfpUMXk5kGvVKrxPDYlB3kQVp6QzTbod7aF5rXOs9IS6Eq30ldodXdsr7ASY3fTZXfiUu7W9Yk/dSVWkHGp2Zo7Q+9Lchb+kL0hyYjy6UIlOU5qPI5ubLSS7XfSZXPuhV7sknVIYsqiga9UnE2PRPm3NCYFeS9oXldLueGxuYtW+DOpMrn5N1bCu0gd1Hlc1DudZCZpqNb1OLELOhF5Hbg74FU4P83xnwlVutSaqWZ7TdvjgR50Oo77wtxdmCUqZm5fnN3VhrVfhd7qrxU5jup8jup9jmp9GlXi1oeMfkWiUgq8I/AbUA78KaIPGGMORaL9SmVKMOjUzTbIX7+bXShfvMCN+/ZUkSVzwrzKp8TrzNDD4SqmIpVc2EPcNoY0wwgIj8G7gY06NWqs1C/+ext4Lx+83Kv1W++p8prB7mLKr+T4pwsPdVfJUysgr4UOBv1dzuwN0brUiom/sfvWvjOy810DI/PKy/MsfrNrZa5w+43d1LhdZCRpmPO1coTq6BfqOli5i0gcj9wP0BFRUWMqqHU5fO7M9lbnR8ZzVKl/eZqlYrVN7YdKI/6uwzoiF7AGPMQ8BBAfX39vJ2AUivBndtKuHNbSaKrodSSxep35ptArYhUiUgG8BHgiRitSyml1EXEpEVvjJkWkQeAp7CGV37PGHM0FutSSil1cTHrbDTG/Ar4VazeXyml1OLoEAGllEpyGvRKKZXkNOiVUirJadArpVSS06BXSqkkJ8Yk/lwlEekF2i7z5T6gbxmrs1y0Xu+M1uud0Xq9M8lar3XGGP+lFloRQb8UItJgjKlPdD3Op/V6Z7Re74zW651Z6/XSrhullEpyGvRKKZXkkiHoH0p0BS5A6/XOaL3eGa3XO7Om67Xq++iVUkpdXDK06JVSSl3Eqgp6EfmQiBwVkbCI1J/33F+IyGkROSEi74kqv90uOy0in41DHX8iIgfsW6uIHLDLK0VkLOq5b8e6LufV60sici5q/XdEPbfgtotTvf5WRI6LyCEReUxEcu3yhG4vuw5x/e5cpB7lIvKCiDTa3/8/tcsv+JnGsW6tInLYXn+DXeYVkWdE5JR9nxfH+myI2h4HRGRERD6TqG0lIt8TkR4RORJVtuD2Ecs/2N+3QyKya9kqYoxZNTdgE7ABeBGojyrfDBwEMoEqoAnr8sip9uNqIMNeZnMc6/t3wBfsx5XAkQRuuy8B/8cC5QtuuzjW691Amv34b4C/WSHbK6HfnfPqUgzssh+7gZP257bgZxrnurUCvvPKvgp81n782dnPNEGfYRewLlHbCrgB2BX9Xb7Q9gHuAH6NNUPfVcDry1WPVdWiN8Y0GmNOLPDU3cCPjTETxpgW4DTWBOWRScqNMZPA7CTlMSciAnwYeCQe61uCC227uDDGPG2Mmbb/fA1rNrKVIGHfnfMZYzqNMW/ZjwNAI9a8zCvV3cAP7Mc/AH4vQfW4BWgyxlzuyZhLZox5GRg4r/hC2+du4GFjeQ3IFZHi5ajHqgr6i1hoMvLSi5THw/VAtzHmVFRZlYjsF5GXROT6ONUj2gP2T8LvRf2cTuQ2Ot99WC2aWYncXitpu0SISCWwE3jdLlroM40nAzwtIvvEmgcaoNAY0wnWTgooSEC9wJrZLrqhlehtNetC2ydm37kVF/Qi8qyIHFngdrHW1IUmI7/kJOUxrOM9zP+SdQIVxpidwJ8BPxKRnKXW5R3U61vAemCHXZe/m33ZAm+1rEOxFrO9ROTzwDTwQ7so5tvrUtVeoCyhQ9RExAU8CnzGGDPChT/TeLrWGLMLeC/wJyJyQwLq8DZiTWF6F/Azu2glbKtLidl3bsVNZ2+MufUyXnaxycgvOkn55bhUHUUkDfh9YHfUayaACfvxPhFpAuqAhqXWZ7H1iqrfd4An7T8vOZF7rOslIvcCdwK3GLuzMh7b6xJivl3eCRFJxwr5Hxpjfg5gjOmOej76M40bY0yHfd8jIo9hdXl1i0ixMabT7nroiXe9sHY8b81uo5WwraJcaPvE7Du34lr0l+kJ4CMikikiVUAt8AaJm6T8VuC4MaZ9tkBE/CKSaj+utuvYHIe6zK4/uq/v/cDsKIALbbt41et24M+Bu4wxo1HlCd1erKAJ7u3jPd8FGo0xX48qv9BnGq96OUXEPfsY68D6EaztdK+92L3A4/Gsl23eL+pEb6vzXGj7PAF83B59cxUwPNvFs2TxPgq9xCPY78fa600A3cBTUc99HmuUxAngvVHld2CNUmgCPh+nen4f+PR5ZR8AjmKN3ngLeF+ct92/AIeBQ/YXqvhS2y5O9TqN1S95wL59eyVsr0R9dy5Qj+uwfsIfitpOd1zsM41Tvartz+eg/Vl93i7PB54DTtn33jjXywH0A56osoRsK6ydTScwZWfXH15o+2B13fyj/X07TNTIwqXe9MxYpZRKcsnSdaOUUuoCNOiVUirJadArpVSS06BXSqkkp0GvlFJJToNeKaWSnAa9UkolOQ16pZRKcv8L+9uEI6hYLkwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(fSoln[:,4],fSoln[:,5])\n",
    "\n",
    "plot(fSoln[:,0],fSoln[:,1], 'or')\n",
    "figure()\n",
    "plot(fSoln[:,4], tpts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical integration of the 3-body problem\n",
    "\n",
    "Now let's extend our integrator to solve simultaneously for the orbits of three gravitationally interacting bodies.  \n",
    "\n",
    "### Two planets orbiting a central star\n",
    "#### Exercise 7.4 \n",
    "As you build your integrator, you will want to run it on a solvable problem so you can tell whether it is working correctly. Therefore, the first case you will look at is that of two planets orbiting a central star.  The star ($M_{1} = 2000$) is initially located at the origin.  The two planets ($M_{2} = M_{3} = 1$) should be set up to undergo  circular (or approximately circular) orbits with radii of 100 and 200, respectively.  Define your initial conditions in the cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Defining constants\n",
    "gGrav = 1 ### FIX\n",
    "epsilon = 1 #### FIX\n",
    "\n",
    "# Masses\n",
    "m1 = 2000 #### FIX\n",
    "m2 = 1 #### FIX\n",
    "m3 = 1 #### FIX\n",
    "\n",
    "# UNCOMMENT AND FILL IN initial conditions below\n",
    "#Initial Conditions\n",
    "\n",
    "r1_0 =  array([0,0])\n",
    "drDt1_0 = array([0,0])\n",
    "\n",
    "\n",
    "\n",
    "r2_0 = array([0,100])\n",
    "r2_0_r = findDist(r1_0[0],r1_0[1],0,r2_0[0],r2_0[1],0)\n",
    "r2_0_v = (gGrav*m1/r2_0_r)**.5\n",
    "drDt2_0 = array([r2_0_v,0])\n",
    "\n",
    "\n",
    "r3_0 = array([0,200])\n",
    "r3_0_r = findDist(r1_0[0],r1_0[1],0,r3_0[0],r3_0[1],0)\n",
    "r3_0_v = (gGrav*m1/r3_0_r)**.5\n",
    "drDt3_0 =  array([r3_0_v,0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Follow the set out steps outlined below to build your 3-body integrator. For now, we will still confine all motion to the x-y plane.  \n",
    "\n",
    "* Paste your 2-body integration from above (`deriv_2body_2D`) into the cell below. \n",
    "* Rename it as `deriv_3body_2D`\n",
    "* Update it to provide the correct position and velocity vector of a 3rd body of mass $M_{3}$, as well.   \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define deriv_3body_2D below:\n",
    "#You'll need xDoubleDot and yDoubleDot defined for each body (6 acceleration definitions).\n",
    "\n",
    "# define my bois\n",
    "def deriv_3body_2D(f,t): \n",
    "    \n",
    "    x_1 = f[0]\n",
    "    y_1 = f[1]\n",
    "    \n",
    "    dxDt_1 = f[2]\n",
    "    dyDt_1 = f[3]\n",
    "    \n",
    "    x_2 = f[4]\n",
    "    y_2 = f[5]\n",
    "    \n",
    "    dxDt_2 = f[6]\n",
    "    dyDt_2 = f[7]\n",
    "    \n",
    "    x_3 = f[8]\n",
    "    y_3 = f[9]\n",
    "    \n",
    "    dxDt_3 = f[10]\n",
    "    dyDt_3 = f[11]\n",
    "    \n",
    "    #Derivative definitions\n",
    "    distance13 = findDist(x_1, y_1, 0, x_3, y_3, 0)\n",
    "    distance12 = findDist(x_1, y_1, 0, x_2, y_2, 0)\n",
    "    distance23 = findDist(x_2, y_2, 0, x_3, y_3, 0)\n",
    "    # Students fill in below:\n",
    "    \n",
    "    d2xDt1_2 = -(gGrav*m2/(distance12+epsilon)**3)*(x_1 - x_2) \n",
    "    d2yDt1_2 = -(gGrav*m2/(distance12+epsilon)**3)*(y_1 - y_2)\n",
    "\n",
    "    d2xDt2_1 = -(gGrav*m1/(distance12+epsilon)**3)*(x_2 - x_1)\n",
    "    d2yDt2_1 = -(gGrav*m1/(distance12+epsilon)**3)*(y_2 - y_1)\n",
    "\n",
    "    d2xDt3_1 = -(gGrav*m1/(distance13+epsilon)**3)*(x_3 - x_1)\n",
    "    d2yDt3_1 = -(gGrav*m1/(distance13+epsilon)**3)*(y_3 - y_1)\n",
    "    \n",
    "    d2xDt1_3 = -(gGrav*m3/(distance13+epsilon)**3)*(x_1 - x_3)\n",
    "    d2yDt1_3 = -(gGrav*m3/(distance13+epsilon)**3)*(y_1 - y_3)\n",
    "    \n",
    "    d2xDt3_2 = -(gGrav*m2/(distance23+epsilon)**3)*(x_3 - x_2) \n",
    "    d2yDt3_2 = -(gGrav*m2/(distance23+epsilon)**3)*(y_3 - y_2)\n",
    "\n",
    "    d2xDt2_3 = -(gGrav*m3/(distance23+epsilon)**3)*(x_2 - x_3)\n",
    "    d2yDt2_3 = -(gGrav*m3/(distance23+epsilon)**3)*(y_2 - y_3)\n",
    "    \n",
    "    #addition\n",
    "    \n",
    "    d2xDt1_12_13 =  d2xDt1_2 + d2xDt1_3\n",
    "    d2yDt1_12_13 =  d2yDt1_2 + d2yDt1_3\n",
    "    \n",
    "    d2xDt2_21_23 =  d2xDt2_1 + d2xDt2_3\n",
    "    d2yDt2_21_23 =  d2yDt2_1 + d2yDt2_3\n",
    "    \n",
    "    d2xDt3_31_32 = d2xDt3_1 + d2xDt3_2\n",
    "    d2yDt3_31_32 = d2yDt3_1 + d2yDt3_2\n",
    "  \n",
    "\n",
    "    return array([dxDt_1, dyDt_1, d2xDt1_12_13, d2yDt1_12_13, \n",
    "                  dxDt_2, dyDt_2, d2xDt2_21_23, d2yDt2_21_23,\n",
    "                  dxDt_3, dyDt_3, d2xDt3_31_32, d2yDt3_31_32],  float)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Integrate your equation of motion from t = 0 to t = 1000 with 1000000 steps.\n",
    "# HINT: It may take a while to integrate, so try using a shorter time period (0 < t < 10) until you are sure everything is working.\n",
    "# BE SURE to update f_init to reflect that there are now three bodies!\n",
    "fInit = [r1_0[0],r1_0[1],drDt1_0[0],drDt1_0[1],\n",
    "         r2_0[0],r2_0[1],drDt2_0[0],drDt2_0[1],\n",
    "         r3_0[0],r3_0[1],drDt3_0[0],drDt3_0[1]]\n",
    "tpts = linspace(0,421,100000)\n",
    "fSoln = odeint(deriv_3body_2D, fInit, tpts)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0xcafce80>]"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD8CAYAAACVZ8iyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xd4VFX6wPHvmUnvvRMSQujd0EFQaSoKrIsFddVV0RUUxd+uvfd1RUSxgLqWtSC7iKigFJHegvQSEkJCeu/JJFPO748MQ5BAwCRzM8n5PM88uffMnbkvEOad04WUEkVRFKVj02kdgKIoiqI9lQwURVEUlQwURVEUlQwURVEUVDJQFEVRUMlAURRFQSUDRVEUBZUMFEVRFFQyUBRFUQAnrQO4UEFBQTImJkbrMBRFURzG7t27C6WUwRdyrcMkg5iYGBITE7UOQ1EUxWEIIdIv9FrVTKQoiqKoZKAoiqKoZKAoiqKgkoGiKIqCSgaKoigKLZAMhBCdhBDrhRBHhBCHhBBzrOUBQog1Qohk609/a7kQQiwQQqQIIfYLIQY1NwZFURSleVqiZmACHpZS9gSGAbOEEL2AR4F1Usp4YJ31HOBKIN76mAm81wIxKIqiKM3Q7HkGUsocIMd6XCGEOAJEAlOAsdbLPgV+BR6xln8m6/fb3C6E8BNChFvfR1E0U2syU2EwWR9GKg0mak0Wak0W6swW6kynHmbqzBb2nCxlfVI+zjodTnqBk16Hs06g1wvcnPR4ujrh7+FMgKcrgV4u+Hu4EODpTJCXK2G+bkT5e+Dj5oQQQus/uqK07KQzIUQMMBDYAYSe+oCXUuYIIUKsl0UCGQ1elmktOysZCCFmUl97IDo6uiVDVToAo9lCbpmB3HIDhRW1FFbWUlBZR2FlLUWVtRRW1lFcVUeFwUi5wUSdyfKH7mPgj73ufPQ6Qd9IX/pF+dIvyo++kb5E+bvj6eow80QVB9Niv1lCCC/gf8CDUsry83zbaewJ2diFUspFwCKAhISERq9ROi4pJTllBlILqkgvriKrpIas0hqySmrILq0ht9yA5Xe/NUKAv4cLQV4uBHm50jvCBx93Z7zdnPB2dcLbzXrs5oyXqxOuzjpc9DrcnHW46PW4OOlwcdLhrBfohEACZovEZLZgskiMZgtmi8RgtFBZa6SkykhxVR3F1fWJp7iqjoKKWnLLDJwsrqbGaG70z2a2SPZmlLI3oxRofBKpu7OeUfFBjOoaxJhuwcQEebbsX7DSobRIMhBCOFOfCL6QUi6zFuedav4RQoQD+dbyTKBTg5dHAdktEYfSPtWZLKTkV5KcX8HxgipSCypJLajiRGHVGR+mTjpBuJ8bEb7uDIsLJMrPnUh/d8J83Qn2ciXIy4UATxec9G1/EF25wcjRnAr2Z5ZyIKuMA5llpBZWnXFNjdHMmsN5rDmcB0CXIE+6BHvSLdSb7mHedAv1pkuwJ65Oei3+CIqDaXYyEPVVgI+AI1LKeQ2eWgHcBrxq/fldg/LZQoivgaFAmeovUE4pqarjSE45h62PIzkVpORXYDTXf8XXCYjy96BLsCfDugTSJbj+AzAm0JNQHzf0uvbR/u7j5syQ2ACGxAY0+rzFIsmrMLA7vYQtKYUUVtbhpBOkFlTxa1IBJmuVSK8TxAZ50j3Um57h3vSN8qNfpC/+ni72/OMoDkDU9+M24w2EGAVsAg6ArfH0cer7Db4BooGTwHQpZbE1ebwDTAKqgTuklE2uQJeQkCDVQnXtS63JzOHscvacrG8O2ZNRQkZxje35EG9Xeob70CvCh57hPvQI86ZzoIf6ptuEOpOFtKIqknIrOJZXQVJuBUl5FaQXVduu6RTgTj9rYugX5Ue/KF/VH9EOCSF2SykTLuja5iYDe1HJwPGVVtexPbWYHSeK2HOylMPZ5dSZ678/RPi6MSDaj/5RfvSO8KVnuDeBXq4aR9y+lBuMHMwsY39WGfszS9mfWUZmSX3y1esEfSJ8GBwTwODYAAbHBBCgag8OTyUDpU0oqzay40QR21OL2ZZaxNHccqQEVycd/aP8GBhd/xjQyZ8wXzetw+2Qiipr2Z9Vxu60EnamFbM3o9Q2qqpriBdDYgMY1TWIkXFB+Ho4axytcrFUMlA0cWoEzIakfNYnFXAwu8z24X9JZ3+GdwlkWFwg/aP8cHFq+524HVGtycz+zDJ2nihmV1oxiWklVNaa0AnoF+XHpfFBjO4WzIBOfjg7QEd8R6eSgWI3RZW1bEwuYP3RAjYmF1BabUQnYGC0P6PjgxjeJZAB0X6qnd9BmcwW9mWWsuFYIZuSC9iXUYpFgrerEyO7BjG+VyiX9whRHdJtlEoGSqvKKq3h54O5/HQwl13pxUgJQV4uXNotmMu6hzA6Pgg/D/Xh0B6VVRvZerzQ9gUgt9yAXidI6OzPhN5hTOgVSqcAD63DVKxUMlBa3InCKn46mMtPB3PYl1kGQI8wbyb2DmNcz1B6R/igayfDOpULY7FIDmSV2eY6JOVVAPW/F5P7hXNt/0iiA1Vi0JJKBkqLyK8wsGJvNsv3ZnEwqxyA/lG+TOoTzqQ+YcSqGa9KA+lFVaw5nMdPB3NJTC8BoH8nP67tH8HkfuGE+qhBAvamkoHyh1XXmVh9KI9le7LYnFyARULfSF+mDIjgyr7hRPq5ax2i4gAyS6r5YX8OK/ZmczinHCFgWGwg110SxVV9w/BwUXMa7EElA+WiSCnZk1HKVztO8uOBHKrrzET6uTN1YATTBkbSNcRb6xAVB5aSX8mKfdms2JtFWlE1Xq5OXNM/ghsGd6J/lK9atbUVqWSgXJByg5Hle7L4csdJjuZW4OmiZ3K/CP40KJLBMQGqD0BpUVJKdp4oZkliBisP5GAwWuge6s31gzvx50FRah5DK1DJQDmvg1llfLYtje/35VBjNNM7wocZQ6OZMiASL7UkgWIH5QYj3+/L5pvETPZllOLurGfqwEhuHxFD9zBVE20pKhkoZzFbJGuP5PHR5hPsPFGMh4uea/tHMGNoNP2i/LQOT+nADmeX89m2NL7dk0WtycKIuEBuGxHDuJ6h7WbhQa2oZKDYVNaaWJqYwSdb00gvqibSz53bR8Rw/eBO+LqrarnSdpRU1bEkMYPPt6WTVVpDlL87d42K5YbB0bi7qEmLf4RKBgolVXX8e8sJPtmaRrnBxMBoP+4cFcuk3mEOsZ6/0nGZzBbWHsnnw02pJKaXEOjpwh0jY7h1WIzqV7hIKhl0YIWVtSzelMp/tqVTVWdmQq9Q7h0bx6Bof61DU5SLtiutmHfXp7A+qQAvVyduHhrNnaNiCVFzFi6ISgYdUF65gfc3HOernSepM1mY3C+CWZd1VZ1xSrtwOLuc9zcc54f92Tjrddw2IoZ7x8SpZbabYPdkIIT4GJgM5Esp+1jLngXuBgqslz0upVxpfe4x4E7ADDwgpfy5qXuoZNC4smoj7288zr+3nMBolkwbGMl9Y+PoEuyldWiK0uLSi6p4a10yy/dk4e6s585Rsdw5uovq/zoHLZLBpUAl8NnvkkGllPJfv7u2F/AVMASIANYC3aSUje8MbqWSwZlq6sx8sjWN935NoaLWxJT+Ecwd312tBaN0CCn5Fby5JpkfD+Tg4+bEPWPiuHNULG7OqqO5oYtJBi0yqFxKuVEIEXOBl08BvpZS1gInhBAp1CeGbS0RS3tntkiWJmbw5tpj5JXXcnmPEP5vQnd6RfhoHZqi2E3XEG8W3jyIv2WVMW/NMV7/OYkvtqfzyJU9uLZ/hJrV/Ae09rCS2UKI/UKIj4UQp3owI4GMBtdkWsuUJuxKK+badzbz6LIDRPl78M09w/n49sEqESgdVp9IXz6+fTBf3T0Mf08X5ny9l2nvbmW3daE85cK1ZjJ4D4gDBgA5wBvW8sZSdqNtVUKImUKIRCFEYkFBQWOXdAg5ZTU88NUepr+/jeKqOhbcNJD/3jucIbEBWoemKG3C8LhAvp89itf/3I/s0hque28r93+1h7xyg9ahOYxWW3tASpl36lgIsRj4wXqaCXRqcGkUkH2O91gELIL6PoPWibTtMpotLNqYyju/pGCWkgcu78q9Y+PUio+K0gidTjA9oRNX9Q3ng42pfLDhOOuP5vPwhG78ZXiMms3chFarGQghwhucTgMOWo9XADcKIVyFELFAPLCzteJwVL+dLGHygs28/nMSl3YLYt3cMcyd0F0lAkVpgqerE3PHd2PNQ2O4pLM/z31/mCkLN7M/s1Tr0Nq0FvlkEUJ8BYwFgoQQmcAzwFghxADqm4DSgHsApJSHhBDfAIcBEzCrqZFEHUllrYnXfzrKZ9vTCfNxY/FfEhjfK1TrsBTF4UQHevDJHYNZeSCX574/xJSFW7hteAx/n9gdT7Ug41nUpLM2ZH1SPo8vO0BuuYG/DOvM/03sjrebGj+tKM1VYTDyxupjfLotjU7+Hrz+534M7RKodVit7mKGlqpFatqA6joTT3x7gDv+vQtvNyf+97cRPDelj0oEitJCvN2cefba3iyZORwh4MbF23nu+0PU1KlGiVNUXUljv50sYe6SvaQXV3P36FgentBdTZxRlFYyJDaAVXNG89qqo/x7Sxq/JhUw7/r+DFRrd6magVZMZgvzVifx5/e2YjRLvrp7GE9c3UslAkVpZR4uTjw3pQ9f3jWUOpOF6e9v4/0Nx7FYHKPJvLWoZKCBvHIDMz7cwYJfUpg2MIqfHhzNsA7QfqkobcmIrkGsnDOaCb1DeXXVUW77904KKmq1DkszKhnY2abkAq56axMHMst484b+vHF9f9U3oCga8XV3ZuGMQbw0rQ87TxRz1YJNbEkp1DosTahkYCdmi2Te6iT+8vFOAr1c+P7+kUwbGKV1WIrS4QkhuHloZ76bPRJfd2du/WgH7/16HEcZadlSVDKwgwqDkbs/S2TBLylcNyiK5bNG0jVE7TOgKG1JjzAfVsweyZV9w3ntp6PM/moP1XUmrcOyGzWaqJWlFVZx12eJnCis4oWpfbh1WGetQ1IU5Rw8XJx456aB9Inw5Z8/H+V4fiWL/5JAp4D2vzS8qhm0oq0phUxZuIXCylo+v3OISgSK4gCEEPxtbByf3DGE7NIarnlnM4lpxVqH1epUMmglSxMzuPXjnYR4u/LdrJGMiAvSOiRFUS7CmG7BrJg9igAPF2Z8uIOVB3K0DqlVqWTQwqSULFyfwt//u58RcYEsu28EnQM9tQ5LUZQ/ICbIk//9bQR9I32Z9eVvfLgpVeuQWo1KBi3IbJE8/d0hXv85iakDIvjotsFq2KiiODh/Txe+uGsok3qH8eKPR3j++8PtcoKaSgYtpM5k4f6vfuPz7encc2kX5l0/ABcn9derKO2Bm7OehTMG8deRsXy85QSPLtuPuZ0lBDWaqAXUmszM+uI31h7J58mre3LX6C5ah6QoSgvT6QRPTe6Jt5sTb61LpsZoYd71/XHWt48vfSoZNJPBaObe/+zm16QCNXRUUdo5IQQPje+Gu4ueV1cdxWA0886Mgbg6Of6aYi2S0qwb3ucLIQ42KAsQQqwRQiRbf/pby4UQYoEQIkUIsV8IMaglYtBCTZ2Zuz9LZMOxAl75U1+VCBSlg7h3TBzPXdubNYfzmPXFbxjNFq1DaraWqt98Akz6XdmjwDopZTywznoOcCX1W13GAzOB91ooBruqM1n42xe72ZxSyGvX9eOmIdFah6Qoih3dNiKGF6b0Zu2RfOZ+s8/h+xBaJBlIKTcCv5+VMQX41Hr8KTC1Qflnst52wO93+yW3eWaLZO43e/k1qYCXpvbl+oROWoekKIoGbh0ew2NX9uD7fdk8vuyAQ48yas0+g1ApZQ6AlDJHCBFiLY8EMhpcl2ktc4gZHVJKnv7uID/sz+GRST2YMVTVCNoLo9nIoYJD7Mvdx/68/RzIP8CB/APkVua2yPsPjRzKhLgJDI4YTEJEAmFeYQghWuS9Fe3cMyaOqjozC9Yl4+nqxNPX9NI6pD9Eiw7kxn77G02nQoiZ1DclER3dNj5031xzjC92nOTeMXH8bWyc1uEoF8lkMbEjcwdfHfyKb49+S3ZFtt3uvSNrBzuydpzz+fiAeOYMncOUHlOI9I5UicKBPDQungqDkY+3nCDK352/jorVOqSL1prJIE8IEW6tFYQD+dbyTKBhu0oU0Oj/SCnlImARQEJCgub1r2W/ZbLglxSmXxLFI5O6ax2O0oQ6cx0/p/zMvO3z+DXt1wt6Ta/gXgwMG0i/0H70DelLn5A+RPpEohMX16JqspgoNZRSUlNCZnkm2zK3sTZ1LevT1p/zNcnFycxeNZvZq2afUf7w8Ie555J76BrQVSWINkoIwVNX9yKn1MALPx4m0t+dib3DtA7rooiWWrNbCBED/CCl7GM9fx0oklK+KoR4FAiQUv5DCHE1MBu4ChgKLJBSDmnq/RMSEmRiYmKLxPpH7Eor5ubFO7iksz+f/nWImlDWBlXVVfHRno94av1TlNeWn/O6fqH9mNFnBlN7TCU+MP6iP+hbkkVayCrPYlf2Lj787UNWpaxq8jWXx17OK1e8QkJEgqaxK2erqTNz0+LtHM0t5+uZwxnQyU/TeIQQu6WUCRd0bUskAyHEV8BYIAjIA54BlgPfANHASWC6lLJY1H+1eYf60UfVwB1SyiY/5bVMBulFVUxduAU/Dxe+vW8Efh4umsShnElKyaaTm5j781x25+xu9JrhUcOZO3wuk7tNxs3Jzc4R/nFmi5kjhUdYemgpz298/rzXPjz8YR4Z+QjBnsF2ik45n8LKWqa9u4WaOgs/PjCKUB/tfu/sngzsQatkUFNnZtq7W8gpM/DtfSPoEuxl9xiU00wWEyuSVnDHd3c0+u1/Wo9pPDPmGfqF9mt3TSpF1UX8lPIT9/54L5V1lY1eM73XdF654hXiAlR/lpaO5VUwdeEWeob78NXdwzRrSVDJoIVIKfm/pftZtieTf98+mLHdQ5p+kdLiTtUApnw9hVJD6RnPBbgH8P7V7zOlxxRc9B2rxpZems57ie/x2pbXGn3+nkvu4aXLXyLQI9DOkSkAP+zPZvaXe7hteGeem9JHkxhUMmghS3ad5JH/HeCBK+KZO76bXe+tQH5VPvf+cC/fHv32jPI4/zg+uvYjLu18abv79v9HVdZV8vm+z7lv5X2NPv/5tM+5qc9N6HWOv2yCI3nxh8N8uPkEb97QX5M9z1UyaAGHs8uZ+u4WhsYG8MkdQ9Dr1IeOvWxK38Sln1x6Vvl3N37HNd2uUQmgCcU1xSzavYjH1j121nO39ruVNye+qWoLdmIyW5jx4Q4OZZWxcs5ou+9topJBMxmMZqa8s4Xi6jp+mjOaQC9Xu9y3I7NICx8kfnDWN9uHhj3Es2OfxcfVR6PIHNvB/IP89bu/sit71xnlYV5hbLpjE10DumoUWceRXVrDpPkbiQvxYuk9w3Gy4yqnF5MM1Li0Rsxbc4ykvAr++ed+KhG0MpPFxBtb30D/vP6MRLD8huVYnrYwb+I8lQiaoU9IH3bevZPKxyp5+fKXbeW5lbnEvx2P0/NOHMw/eJ53UJorws+dl6b1Zc/JUt7+JUXrcM5J1Qx+Z3tqETct3s5NQ6J5eVrfVr9fR2WRFhbuXMgDPz1gK3N3cidxZiK9gh1zOr8jMFvMfHPoG2Ysm3FGubPOmSOzjqhRSK1o7pK9LN+bxbL7Rtpt/oFqJvqDaurMTJi/AZ0QrHxgNJ6uaruH1rD6+Gom/mei7dzX1Zedd++kW6DqpLcXKSWrUlZx9ZdXn1E+stNIlt+4nCCPII0ia78qDEYmvLkRX3dnvr9/lF02xVHNRH/Q278kk1Fcw6t/6qcSQStIL01HPCfOSASH7jtE6aOlKhHYmRCCq+KvwvK0ha+u+8pWviVjC8GvB/PY2scwW8waRtj+eLs58/yUPhzNrWDxplStwzmLSgZWx/IqWLQxlesGRTE8To20aElGs5E5q+YQ81aMrey7G79DPiNVk5DGhBDc2OdGjE8ZefGyF23lr255FacXnEjM1m4JmPZofK9QruwTxltrk0krrNI6nDOoZEB9lfnJ5QfxcnPi8at6aB1Ou/Jbzm+4vOjCgp0LgPrRQcanjFzb/VqNI1MactI58cSlT5D7cC6jokfZygcvHsw1X11DjbFGw+jal2ev7Y2LXsez3x/SOpQzqGQA/Hggh50ninl0Ug81eqiFmC1mHlj1AJcsusRWlvpAKvMmzsNJp5rg2qpQr1A23bGJX2/71Vb2w7Ef8HjZgz05e7QLrB0J9XFjzrh4fk0qYOOxAq3DsenwyaDWZOafPyXRI8yb6WrHshaRUZaB0wtOvL3zbQBeG/calqctxPo73hrvHdWYmDFUPFbB1B5TbWWDFg3isbWPYZGOv9+v1m4d3pnOgR689OORNrNdZodPBl9sP8nJ4moeu6qnmmXcAtamriV6/umNiI7OOso/Rv5DzRp2QF4uXnx7w7esuvn0stqvbnmVyHmR510iXGmaq5Oex67sQVJeBUt2ZTT9Ajvo0MmgwmDk7V+SGdU1iEvj1VC65pBS8tQvTzH+8/EAjOsyDsMTBroHqU2AHN2krpPIfCjTtgR4bmUuvq/6crTwqMaRObaJvcNI6OzPgnXJ1Jq0H7nVoZPBZ9vSKak28o9J3dU312YwWUyM+WQML26qH40yf+J81ty6Blcn1f/SXkT6RFL+aDl/6f8XW1nPhT1Zl7pOw6gcmxCCh8Z3I7fcwDeJmVqH03GTQU2dmY83n2BMt2D6RWm7G5Ejq6qrImpeFJtObgJg4+0bmTNsjsZRKa3BWe/Mp1M/5e0r37aVjft8HF8e+FLDqBzbiLhALunsz3vrU6gzadsX0+rJQAiRJoQ4IITYK4RItJYFCCHWCCGSrT/9WzuO3/tq50mKquqYfblaqOuPKjWU4vWKF3lVeQAk35/M6M6jNY5KaW2zh8xm5YyVtvObl93MvG3zNIzIcQkhmHNFPNllBv73m7a1A3vVDC6TUg5oMC36UWCdlDIeWGc9txuT2cLiTakMiQ1gcEyAPW/dbpQZyvB/7XQOz3woU62A2YFcGX8lO+/aaTt/ePXDvLntTQ0jclyj44PoE+nDx5tPoOXyQFo1E00BPrUefwpMPc+1LW7tkTxyygzcNUoNdfwjymvL8XvtdNNawd8LiPSJ1DAiRQuDIweze+bpvafnrp7Lu7ve1TAixySE4PYRsSTnV7L1eJFmcdgjGUhgtRBitxBiprUsVEqZA2D9adf9JD/blk6knztX9Ay1523bBaPZSNcFp2sAOQ/nqEXNOrBB4YPYc8/pyWizVs5i6aGlGkbkmCb3CyfA04V/b0nTLAZ7JIORUspBwJXALCHE2VtYnYMQYqYQIlEIkVhQ0DIz9VLyK9h6vIgZQ6PVvIKLJKVk6pKpFFTX/1ukzUkjzCtM46gUrQ0IG8D629bbzq//7/XszNp5nlcov+fmrGfGkGjWHc0js6RakxhaPRlIKbOtP/OBb4EhQJ4QIhzA+jP/HK9dJKVMkFImBAcHt0g83yRm4qQT3DBYzTa+WC9sfIGVyfUdh7/N/I3Ofp01jkhpK8bGjOWbP39jOx/64VBOlp3UMCLHc8PgTkgJy/dkaXL/Vk0GQghPIYT3qWNgAnAQWAHcZr3sNuC71ozjFItFsmJvNmO7BxOk1iC6KL+m/cozvz5jOx8YPlDDaJS2aHrv6Tw5+knbefzb8RjNRg0jciydAjwYHOPPsj1ZmnQkt3bNIBTYLITYB+wEfpRS/gS8CowXQiQD463nrW5nWjG55QauHaA6Oy9GXmUel316me085f62u3Wfoq3nLnuOoZFDAagz1/HgTw9qHJFj+dOgKFILqtifWWb3e7dqMpBSpkop+1sfvaWUL1nLi6SUV0gp460/i1szjlO+25uNh4uecT3t2l/t0KSUTPpiku38s6mfqa0RlXPSCR0rbz49B+HdxHfVLOWLcFXfcFz0Or7fl233e3eYGcgWi2TN4Twu7xGCh4taQvlCfX3wa/bm7gVgYNhAbul3i8YRKW1dgHsAq29ZbTsf9/k4qo3adIo6Gl93Z4bFBfLL0Ua7UVtVh0kGh7LLKays5fIeqlZwoQqrC8/YOH31ravVGk7KBRkfN57bB9xuO39kzSPaBeNgxvUMIbWwitSCSrvet8Mkg/VJ+QgBl3ZrmVFJHcGja09PDF80eZGaT6BclPkT59uO39n1DslFyRpG4zhOfWFdd8S+tYMOlQz6RfmpUUQXKLkomY/2fGQ7v3PQnRpGozgiXzdfPpj8ge387u/v1jAaxxHl70G3UC82pRTa9b4dIhlU15nYn1nGSLXR/QWb89PplUfX37YenegQvypKC7tz4J0465wB2JC+gV1ZuzSOyDEMiQ3gt/QSTGb7rWTaIf6H788sw2yRJMTYfXFUh5RaksqqlPrdrTycPRgbM1bbgBSHpdfp+WzaZ7bzv6/5u4bROI4hsYFU1po4klNht3t2iGSwO70EgEHRKhlciH9u+afteNn1yzSMRGkPpveabjvekL6BEyUnNIzGMQyxrqa8M80uo+6BDpIMfksvIS7YEz8PF61DafOq6qr4YPfpdt4JcRM0jEZpD/Q6Pe9c+Y7tfMGOBRpG4xjCfN0I83HjUJb9Jp91iGRwNLeCXhG+WofhEJYdOV0TePWKV9VQUqVFzOh7eojy/B3zMVu03/O3resR7s2RXNVM1GKqak1kldbQLcRL61Acwts7T29p2HCcuKI0h7+7P3/q+Sfb+fq09ee5WgHoEeZDSn4FRjt1Irf7ZJCcXz9xIz7UW+NI2r4yQxm7sutHe7jqXQn1Uvs9KC1nRp/TtYNvDn1znisVgB5h3hjNkhOFVXa5X7tPBqdm8XVVNYMm/XLiF9vx85c9r2EkSns0Pm687Xjxb4s1jMQxdA70ACCj2D5LebT7ZJBTZgAg0s9d40javjWpa2zHquNYaWk+rj70CeljO08rTdMuGAcQ6V//mZVVWmOX+7X7ZJBbZsDX3Rl3F73WobR5/z38X9txw/+0itJSJsdPth2rCWjnF+TpioteR1aJSgYtIqfMQLivm9ZhtHlGs9G2naWzzhknnVrZVWl5gyMH247/IcZ1AAAgAElEQVRPrYarNE6nE4T5utlaN1r9fna5SyOEEJOEEElCiBQhxKNNv+KPKaqqVesRXYDjJcdtxw1HfShKS+oX2s92vDVzq4aROAY/D2fKDfbZLU6TZCCE0AMLgSuBXsBNQoherXGvqloTXq7qW25TGu5XOyh8kIaROAghzv1QzqmTz+m9x39N+1W7QByEt5sT5TXtOBkAQ4AU605odcDXwJTWuFFVrRkPV9Vf0JTM8kzbcbRvtIaROICmPvBVQjgnVydVS78Y3q7OVBhMdrmXVskgEshocJ5pLWtxVXWqZnAhimtOr4Hi66pmayutp+EKuFps/O5IPFz1VNfZZ7a2Vsmgsa9OZ/1WCCFmCiEShRCJBQUFf+hGZrPESdfu+8mbrdZUazv2dlUT9JTWE+YVZjs2S7UsxfmIRj8qW4dWn5KZQKcG51HAWTtASykXSSkTpJQJwcF/bIcyIcCivn00qdZ8OhmovQuU1nRqfwOoH8WmtA1a/a/fBcQLIWKFEC7AjcCK1riRTidUMrgAbk6nh98aTPYZyqZ0TJV1p/f2VUOYz8+e3U+a/EtIKU1CiNnAz4Ae+FhKeag17qUXApNFJYOmeLucbhqqMdpnkovDkvL8/0vVl4/zKqopsh07653Pc6VSZ7LgrLdPRtAsLUspVwIrW/s+Xm5OVNXapzfekfm7n974J6cyR8NIHIT6wFfsoMJgxNvNPgmz3TcO+7o7U1qt2iWbEuMXYztWO1EpraVhE1Gop1oVtykVBhPebvb5zt4hkkGZnSZtOLI4/zjb8doTazWMRGnPGi5O13AVU6VxKhm0IH8PF4qr6rQOo81rONxvZ9ZODSNR2rM9OXtsx72De2sYiWPIqzAQ7G2fiXrtPhlE+LmTU1aDWXUin5cQgtHRo23nFbX2225P6Th2ZO2wHQ+PGq5hJG1fVa2J0mojEXZafr/dJ4NOAe4YzZL8CjVcsilXxF5hO96YvlHDSJT26v3E923HDVcwVc6WU1Y/qs9ee7G0/2Tgf2q3IDVcsinjuoyzHX996GsNI1Hao5yKHNuM4xi/GDycPTSOqG079ZmlagYtJCbQEzi9/aVybsOihtmO/7P/P2rdGKVFrUw+PZJ89uDZGkbiGI7l1TfVxttpy952nwyi/N3xcnXiSE651qG0eXqdnjsG3GE735qh1ptXWs5rW16zHU/rOU3DSBxDUm4FoT6u+Hm42OV+7T4Z6HSCnuHeHFbJ4ILcPuB22/HLm1/WLhClXUkrTSO5OBmAcK9wuvh30Tiitu9obgXdw3zsdr92nwwAeoX7cCSnAosaUdSk0dGjcdHXfxNZmbySUkOpxhEp7UHDjuNXx72qYSSOwWA0k5JfSc9w+60g3CGSQd8oPyprTRzLV8MlmyKE4LVxp6vzb257U8NolPagorbijCai6b2maxiNY9iXUUqd2cLgzgF2u2eHSAZDY+v/QrcfL2riSgXg7kF3246f3/i8WrhOaZaGtYLnxz6Pu7N9Rsc4ssT0EgAu6ezfxJUtp0Mkg04BHkT6ubPjRHHTFyt4unjy8uWn+wte2PiChtEojqy8tpx/rP2H7XzOsDkaRuM4dp4opluoF/6e9uk8hg6SDACGdglgx4li1W9wgWYPOT3075XNr5Bfla9hNIqjemb9M7bjly9/GR9X+3WIOiqD0czOE8UM7xJo1/t2mGQwplswxVV17MlQHaIXwtvVm0+nfmo7n7ZEDQVULk5qSSrzd8y3nc8dPlfDaBzHtuNF1BjNXN7Tvqu6dphkMLZ7CE46werDuVqH4jBu6XeLbQjg1oytrD6+WuOIFEdhkRaGfjjUdr5yxkpcneyz4JqjW3c0Dw8Xva2v015aLRkIIZ4VQmQJIfZaH1c1eO4xIUSKECJJCDGxtWJoyNfdmeFxgaw5nGeP27ULOqFjxY2ndyOd+J+JagE75YK8tf0tCqsLAZgQN4Er46/UOCLHYLFI1h3JZ1TXINyc9Xa9d2vXDN6UUg6wPlYCCCF6Ub/ncW9gEvCuEMIuf+rxvUJJLagiKVd9oF2o3iG9eX3867bzoR8OVctUKOd1KP8Qc1efbhL6+jq1ztWFSkwvIafMwJV9w5q+uIVp0Uw0BfhaSlkrpTwBpABD7HHjq/qG46QTLPst0x63azceGvYQ/UP7A3Ck8MgZY8YVpaHy2nL6vNfHdr7m1jVnbKmqnN+3ezLxcNEzsXf7SwazhRD7hRAfCyFO/UZEAhkNrsm0lp1FCDFTCJEohEgsKChodjBBXq6M7R7Csj1ZmMyWZr9fR6HX6Vl96+n+gsfWPcam9E0aRqS0RRZpYdiHpxc7fPrSp89YCVc5P4PRzA/7c5jYOwwPF/tvT9+sZCCEWCuEONjIYwrwHhAHDABygDdOvayRt2q03UFKuUhKmSClTAgODm5OqDZ/viSSgopaNqUUtsj7dRQhniHsunuX7fzSTy7lWNExDSNS2pq7VtzFkcIjQP0KuM+OfVbbgBzMmsN5VBhMTB3Y6HfjVtesZCClHCel7NPI4zspZZ6U0iyltACLOd0UlAl0avA2UUB2c+K4GJf3CCXA04Uvtp+01y3bjYSIBJb8eYntvPs73dX8AwWAlza+xL/3/tt2vubWNQjR2Pc+5Vw+25ZGdIAHo7oGaXL/1hxNFN7gdBpw0Hq8ArhRCOEqhIgF4gG7bbrr4qTj5qHRrDuaR3pRlb1u225c3/t63pr0lu089F+htlEjSsf03q73eHL9k7bzgr8X4OVinzX424tD2WXsSivh1mGd0eu0SaKt2WfwTyHEASHEfuAy4CEAKeUh4BvgMPATMEtK6/ZHdnLLsM446QSfbE2z523bjQeGPsATo5+wnQe/HkxBVfP7dBTH887Od7hv5X228/QH0wny0OabrSP7bGs67s56rk/o1PTFraTVkoGU8lYpZV8pZT8p5bVSypwGz70kpYyTUnaXUq5qrRjOJdTHjcn9IvhmVwblBqO9b98uvHDZC2ckhJB/hZBZrkZpdSSvbHqF+1fdbzs/MecE0b7RGkbkmPLKDXy7N4tpgyLx9XDWLI4OMwP59+4aHUtVnZl/b07TOhSHJITgxctf5KXLX7KVdXqzE/vz9msYlWIPUkpuX347j//yuK0s86FMYvxitAvKgX2wIRWzRXLvpXGaxtFhk0HvCF8m9Arlw82plNWo2sEf9fjox1k0eZHtvP/7/Vl2ZJmGESmtqdZUS/zb8Xy67/S6VXn/l0ekjzYjYBxdQUUtX+5MZ+qASKIDPTSNpcMmA4A54+KpMJj4ePMJrUNxaHdfcjdrb11rO7/um+u4a8VdaqZyO3Oy7CRuL7lxvOQ4AH1D+lL9eDUhniEaR+a4Fm9Kpc5kYdZl2tYKoIMng94RvkzsHcrHm09QXFWndTgO7YouV3Bk1hHb+Ud7PkL3vI4yQ5mGUSktZfnR5XSe39l2PmfoHPbdu09tVNMMGcXVfLI1jakDI+kSrP3oqw6dDAD+b0J3qo1m3lyjJlA1V4+gHpQ/Ws6ITiNsZX6v+bEyeaWGUSnNUWuqZeJ/Jp6xhPmqm1cxf9J8NY+gmV7/OQkB/H1id61DAVQyID7Um1uGRvPFjnS1gF0L8Hb1ZvMdm3nnyndsZVd/eTXDPhyGwWTQMDLlYu3O3o3bS25nLF2eNTeLSV0naRhV+7A3o5QV+7K5e3QXwn3bRu2qwycDgAfHdcPbzZkXfjis2rlbgBCCWUNmndFstCNrB+4vubMiacV5Xqm0BQaTgclfTiZhcYKt7MnRT2J6ykSEd4SGkbUPFovkue8PEeTlwr1jte8rOEUlA8Df04UHx8WzOaWQHw/kNP0C5YL0COqB6SkTz419zlY25espiOcEeZVqX4m26OuDX+P+kjs/Jv9oKzs2+xgvXP4Cep1919dvr77YeZI9J0t5/KqeeLnaf0G6c1HJwOrWYZ3pG+nLsysOUVqtOpNbil6n5+kxT5N8fzK+rr628rA3wpi+dLpqOmojDuUfQjwnuOl/N9nK3pjwBuanzcQHxmsYWfuSV27gn6uOMrJrINM0WpDuXFQysHLS63j1ur6UVBt58ccjTb9AuShdA7pS+mjpGTun/ffwf3F/yZ0n1j2B2WLXFUkUq6zyLNxedDtjD4JhUcMo/kcxc4fPRSfUR0RLeu77Q9SaLbw4tW+b64BX/9IN9I7w5Z5Lu/Df3ZlsPKbW2mkN13S/BsMTBl647AVb2cubX8bpBSeeXv80JotJw+g6jqzyLEJeDyHqzShqzbW28oN/O8i2O7epDWlawYp92aw8kMucK+KJDfLUOpyzCEfpME1ISJCJiYmtfh+D0czVCzZRWWti1ZxLCfB0afV7dlSVdZU8+cuTvLXjrTPKbx9wO+9d/R5uTm4aRdZ+HS44TO93e59VvuOuHQyJtMuGgx1SdmkNk+ZvJC7Ei6X3DMdJb5/v4UKI3VLKhKavVDWDs7g563nrxoGUVBl55H/71eiiVuTl4sX8SfMpeaSEB4c+aCv/ZO8nuL/kTuS8SI4XH9cwwvZBSsmXB75EPCfOSgSb79iMfEaqRNCKLBbJ/y3dh8kimX/DALslgovVNqPSWJ9IX/4xqTtrDufxxQ61CU5r83Pz481Jb2J4wsCbE9+0lWdXZNP17a6I5wTzt89X/QoXqbC6kFEfj0L3vI6bl91sK/d09iRpdhLyGcnI6JEaRtgxLN6UytbjRTw9uRedA9te89ApqpnoHCwWye2f7GJHahHL7htB7wjfpl+ktAiLtLAudR0T/jOh0edX37KacV3GtbkOuLbAYDLwwoYXeHnzy2c9d9fAu3ht/GsEuAdoEFnHtD21iJs/3MHE3qEsnDHI7r+zF9NM1KxkIISYDjwL9ASGSCkTGzz3GHAnYAYekFL+bC2fBLwF6IEPpZSvXsi97J0MoH5FwWve3oyzk2DFrFH4q/4DuyuoKuDFjS+yYOeCRp9fceMKru52dYce9VJVV8XzG57nn1v/2ejz629bz5jOY1TytLO8cgNXL9iMj7sTK2aP0mROgT2TQU/AAnwA/N+pZCCE6AV8Rf2+xxHAWqCb9WXHgPHU74W8C7hJSnm4qXtpkQwA9pws4YYPtjO0SwCf3DFEsy3pFDiYf5CZ389kW+a2Rp//+4i/8+ioRzvEN9+9uXu5bflt59w/4t2r3uWOgXeoTniNGM0Wblq0ncM55Xw3ayTxod6axGG3ZNDghr9yZjJ4DEBK+Yr1/GfqaxAAz0opJzZ23flolQwAluw6ySP/O8A9Y7rw2JU9NYlBOU1Kyd7cvfxj7T9Ym7r2nNc9OfpJ5gyb0y62Ydyft5/ZK2ez6eSmc17zweQPuLnvzXi6tN126Y5ASsmj/zvAksQMFtw0kGv7a7eEx8Ukg9aqt0QC2xucZ1rLADJ+Vz60lWJoMTcMjmZ/ZhkfbEilS5AnNwxWW/tpSQjBwPCBrLl1DVDfUfr5vs+Zu3ruGde9uOlFXtz0ou08wD2AFy97kZv63oSfm59dY74YyUXJ/Gvrv1j026LzXhfrF8vS6UsZFG7/tmjl3N7fkMqSxAxmX9ZV00RwsZpMBkKItUBYI089IaX87lwva6RM0vjopXNWTYQQM4GZANHR2n4AP3ttbzJKanj824OE+rgxtrva0KOtCPII4qHhD/HQ8IcAyCjL4JO9n/D0r0+fcV1xTTH3rbzvjA3cAVz1rtybcC9Xx1/N8E7D8XJp/bXlC6oK+DH5R7448MV5azcNdfbtzKdTP2Vk9EicdG1nTRvltJUHcnjtp6Nc0z+CueO7Nf2CNkQ1E12EyloT17+/jfSiKpbcM5w+kWqEkSOoNdWyL28f87bNY8mhJc1+P2edM0EeQfi5+eHh7IGrkyvVxmoq6yopqSmhqKao2ffQCR0LJi3gpr43dYg+kPYgMa2Ymz/cQZ9IX764ayhuztov7NcW+gx6A19yugN5HRBPfY3hGHAFkEV9B/IMKeWhpu7RFpIB1I8QmLZwC0aLZOk9w4lpg9PKlQtTWVfJkYIjfLrvUxb/tpg6s/0XKOwf2p/7Bt/HtB7TCPIIUs09DupgVhk3LdpOkLcr/713OIFerlqHBNh3NNE04G0gGCgF9jb41v8E8FfABDwopVxlLb8KmE/90NKPpZQvXci92koyAEjOq+D6D7bh4eLEknuGEeWv7UbWSusxW8yU15ZTWF1IUlESx4qOkV+VT0FVARV1FVQZq6g11eLp4omXixd+rn6EeYUR5hVGz+CedA/sTqBHYIce+trepeRXcv0H23B31vPNvcOJ9Gsbm9WABjUDe2hLyQDqvwnMWLwdf08XlswcTpivGsKnKB1NRnE109/fhskiWXrv8Da3AJ1am8gO+kT68ulfh1BUWceMD7dTUFHb9IsURWk30ouquHHRdmqMZj6/c0ibSwQXSyWDZhgY7c/Htw8mp9TADR9sI7u0RuuQFEWxg+MF9U1D1XUmvrhrKD3DfbQOqdlUMmimIbEBfH7nEAoqapn+/jbSCqu0DklRlFaUlFvBDR9sx2yRfDVzWLsZVaiSQQtIiAngq5nDqDGamf7BNpJyK7QOSVGUVrDnZAk3LtqGTsDXM4fTI8zxawSnqGTQQvpE+rJk5jB0Aq7/YBs7TxRrHZKiKC1o7eE8blq8HW83Z765ZzhdQ1p/cqI9qWTQguJDvVl6zwgCPV245cMdfLc3S+uQFEVpAV/sSGfm54l0C/Xmf38b0S7nF6lk0MKiAz1Ydt8IBnTyY87Xe1m4PkXtlqYoDspikbz+81Ge+PYgY7oF8/XMYQR7t40JZS1NJYNW4Ofhwud3DWHKgAhe/zmJR/63n1qT2qVLURxJhcHIzM93s3D9cW4a0onFf0nAw6X9rgnVfv9kGnN10jP/hgF0DvBgwS8pHMur5L1bBhHu23ZmJyqK0ri0wiru/iyR1MIqnru2N38Z3rndLxWiagatSAjB3Andee/mQSTnVXDN25vZkdr8RcwURWk9G48VMGXhFgora/n8r0O4bURMu08EoJKBXVzZN5zls0bi4+bMjA938PHmE6ofQVHaGLNFMm91Erf9eyfhvm6smD2KEV0df2OkC6WSgZ3Eh3qzfPZILu8RwvM/HGbm57spqbL/KpmKopwtr9zAjMXbWfBLCtcNimLZfSPoFNCxFqBUycCOfNyc+eCWS3hqci9+Tcrnyrc2sV01GymKpjYcK+DKtzaxP7OMN6b351/T+7frjuJzUcnAznQ6wZ2jYvn2vpG4u+i5afF25q1OwmS2aB2aonQoNXVmnl1xiNs+3kmwlyvf3z+K6y6J0joszahkoJE+kb78cP8orhsUxYJfUvjTe1s5lqeWsVAUe9hzsoSrF2zik61p3D4ihuWzRra7GcUXSyUDDXm6OvGv6f1ZOGMQmSU1TF6wmYXrU1QtQVFaSZ3Jwr9+TuK697ZSa7Lw5V1Defba3ri7aL9FpdaalQyEENOFEIeEEBYhREKD8hghRI0QYq/18X6D5y4RQhwQQqQIIRaIjjBmqwlX9wtn9UOXMr5XKK//nKRqCYrSCnanl3DN25t5Z319J/GqB0d3qNFCTWluzeAg8CdgYyPPHZdSDrA+7m1Q/h4wk/o9keOBSc2MoV0I8nJl4c2DbLWEqxds4vWfj1JTp2YuK0pzlNUYeeLbA/z5/a1UGIx8dFsCr0/vj4+bs9ahtSnN6jKXUh4BLnhChhAiHPCRUm6znn8GTAVWNSeO9uTqfuEM6xLASyuPsHD9cZbvyea5a3szrleo1qEpikORUvLjgRye+/4wRZW1/HVkLHPHd8PTteONFLoQrdlnECuE2COE2CCEGG0tiwQyG1yTaS1rlBBiphAiUQiRWFBQ0Iqhti2BXq7Mu34AX88choeLnrs+S+SuTxPJKK7WOjRFcQhHc8u55aMdzP5yD2E+9RPInprcSyWC82jyb0YIsRYIa+SpJ6SU353jZTlAtJSySAhxCbBcCNEbaKwKcc6puFLKRcAigISEhA43ZXdYl0BWzhnNx5tPMH9tMuPmbeCu0bH8bWxXvNQvtaKcpbiqjnlrkvhyx0m83Zx59ppe3DKsM056NVamKU1+okgpx13sm0opa4Fa6/FuIcRxoBv1NYGGA3mjgOyLff+OxFmv454xcUzuH8HrPx1l4frjLNmVwYPjunHj4E7ql1xRqB8l9Pn2dN5ae4yqOjO3DuvMg+O64e/ponVoDqNVvl4KIYKBYimlWQjRhfqO4lQpZbEQokIIMQzYAfwFeLs1YmhvIv3cmX/jQO4YGctLPx7hyeUH+XRrGo9f1ZOx3YM7xEJaivJ7Zotkxb4s3lyTzMniakbHB/HU5F50C/XWOjSHI5qzYJoQYhr1H+bBQCmwV0o5UQhxHfA8YALMwDNSyu+tr0kAPgHcqe84vl9eQBAJCQkyMTHxD8fankgp+flQHq+uOkJaUTWXdPZn7vhujIgLVElB6RCklKw+nMcbq5M4lldJr3Af/j6pO2O7qS9GDQkhdkspE5q+spnJwJ5UMjhbncnC0t0ZvPNLCjllBobEBvDw+G4M7RKodWiK0iqklGxKLuSNNcfYl1FKlyBP5k7oxlV9wtHpVBL4PZUMOhiD0czXO0+y8NfjFFTUMqprELMu68qwLgHqW5LSLlgskjVH8li4PoX9mWVE+LoxZ1w81w2KUv1m56GSQQdlMJr5z/Z03t9wnMLKOvp38uNvY7owvlcYevWtSXFAJrOFH/bn8O6v9bsFRgd48LexcfxpUCSuTmoJiaaoZNDBGYxmlu7OZPHGVE4WV9MlyJO7L+2i/gMpDqPCYGRpYiafbE3jZHE13UK9mHVZV67uG65qAhdBJQMFqB9psepgDu9vOM7BrHKCvFyZMTSam4dGE+rjpnV4inKW9KIqPtmaxtLETCprTVzS2Z97Lu3CuJ6hqk/gD1DJQDmDlJItKUV8vOUE65Py0QvBxD5h/GVYZ4bEqn4FRVsWi2Tr8SI+2ZrGuqN5OOkEk/tFcMfIGPpF+WkdnkO7mGSgprF2AEIIRsUHMSo+iJNF1fxnRzpLdmXw4/4ceoR5c/OwzlzbPwJfd7Vwl2I/+eUGlu7OZMmuDE4WVxPo6cL9l3XllmGdCVE1V7tTNYMOqqbOzIp9WXy6NZ3DOeW4OumY2DuM6xM6MSIuUFXJlVZhtkg2HMvnq50Z/HI0H7NFMrxLIDcO6cTE3mG4Oas+rZakmomUCyal5GBWOUt3Z7B8TxblBhORfu5cNyiSPw2KIibIU+sQFQcnpeRAVhnf7c3m+33Z5FfUEuTlyp8vieKGwZ2IVb9jrUYlA+UPMRjNrDmcx9LdmWxKLkBK6BvpyzX9w7m6XwSRfu5ah6g4kPSiKpbvyea7fVmkFlThotcxtnswfxoUyRU9Q3FWo4JanUoGSrNll9bww/5sftifw/7MMgAGRfsxuV8EV/cLV6ORlEal5Ffy86FcVh/KZV9mGULA0NgApg6I5Mo+4fh6qH4pe1LJQGlRaYVV/Hggh+/3ZXM0t347zv5RvozrGcq4XqH0CPNWI5I6qFNNQD8fyuWng7kcL6gCoH8nP67qE8Y1/SOIUDVKzahkoLSalPwKfjqYy9oj+ezNKAXqV1Qd3yuUK3qGMCQ2QE1sa+fKaoxsTSlkY3IBvyYVkFNmQK8TDI0NYGLvMCb0DiXcVyWAtkAlA8Uu8ssN/HI0n7VH8tmcUoDBaMHNWcfgmABGxwcxsmsQPcN81MgkB2exSA5ll7PhWD4bjhXw28lSzBaJt6sTI7oGMr5XGFf0CFF7B7RBKhkodldTZ2br8UI2JReyJaWQ5PxKAAI9XRjRNYgRcYEMjgkgLthTNSm1cSazhcM55ew8Ucz21GJ2pRVTVmME6gcUjOkWzJjuwQzo5Kc6gds4NelMsTt3Fz1X9Azlip6hAOSWGdiSUshm6+P7ffUb2vl7OHNJ5wAGx/iTEONPn0hf1ayksapaE/szy9iTUcLOE8UkppVQWWsCICbQg0m9wxgWF8Do+GCCvFw1jlZpLc3d3OZ14BqgDjgO3CGlLLU+9xhwJ/Wb2zwgpfzZWj4JeAvQAx9KKV+9kHupmoHjklKSWlhFYlr9B01iegknCus7Gl2cdPSO8KFvpC99In3pG+lLfIiXWoyslZjMFpLyKtiXUcbejBL2ZZSRnF+BxfoxEBfsydAugQyNDWBobCBhvmrUmCOzWzOREGIC8IuU0iSEeA1ASvmIEKIX8BUwBIgA1lK/BzLAMWA89fsh7wJuklIebupeKhm0LwUVtexOLyExrZj9mWUcyi6jqs4MgKuTjp7h9Qmie5g33UK96RbqhZ+HapO+GAUVtSTlVnA0t5wjORUk5ZVzLK+SOpMFqK+l9e/kR/8oPwZ08qN/Jz8CVLt/u2K3ZiIp5eoGp9uBP1uPpwBfSylrgRNCiBTqEwNAipQy1Rro19Zrm0wGSvsS7O3KpD5hTOoTBtR3UqYWVnEwq4wD1se3e7JszRUAQV6udAv1Ij7Ei/hQb2ICPekc6EG4r1uHrUnUmSycLK4mrbCKE4VVpBZWkVZYRXJ+BYWVdbbrgr1d6RHmzW3DO9Mn0pcBnfyIDvBQ/TeKTUv2GfwVWGI9jqQ+OZySaS0DyPhd+dAWjEFxUDqdoGuIF11DvJg6sP5XRUpJdpmBY3kVJOdVkJxXybH8Sv67O9NWiwDQ6wSRfu5EB3jQKcCD6AAPIvzcCPF2I9THlTBfNzxcHLN7rKbOTHZZDdmlNeSUGk4flxk4WVxNRnG1rYkHwM/DmdggT8Z2D6FnuA89wrzpHuat2vqVJjX5P0QIsRYIa+SpJ6SU31mveQIwAV+celkj10ugsa9v52ynEkLMBGYCREdHNxWq0s4IUf8hH+nnzmXdQ2zlp5LEyaJqThZXcbK4mpPFNZwsruangzmUVBvPei9vVydCfFwJ9XEjwNMFfw8X/Dyc8XV3xs/DBT93Z/w8nPF2c8bdWY+bi7b/Z40AAAcTSURBVA53Zz3uzvpm1zqklNQYzVQaTFTWmqiqNVNRa6Sq1kxlrZGSKiNFVbUUV9VRVFlHUVWd9biWcoPpjPcSAoK9XInwc6dvpC9T+kcQE+RJrPWhmtKUP6rJZCClHHe+54UQtwGTgSvk6Q6ITKBTg8uigGzr8bnKG7v3ImAR1PcZNBWr0jE0TBLD4wLPer7CYCSv3EBeee3vftY/DmWXU1pdR1mN8Yxv1efirBe4O+txcdKjE/U1EZ0Q6HTU/xQCARgtFkxmidFswWj9aTJLjBYLTXXN6XUCfw8XgrxcCPB0oXeED4GeLoT4uBHp5064rxsRfu6E+rjh4tQxm8SU1tWsurN1ZNAjwBgpZXWDp1YAXwoh5lHfgRwP7KS+xhAvhIgFsoAbgRnNiUFRfs/brf4bftcQ7/NeZ7FIKmpNlFUbKa2po6TaSFWtiZo6MzVGMwajmWrrcU2dmVqTBSklFikxWzh9LOuPnfU6nPUCJ70OlwbHznodni56PF2d8LI+PG0/9fh7uODr7qwm5ymaam5D6juAK7DG2hG1XUp5r5TykBDiG+o7hk3ALCmlGUAIMRv4mfqhpR9LKQ81MwZF+UN0OoGve31TUTQeWoejKJpSM5AVRVHaqYsZWqoaHxVFURSVDBRFURSVDBRFURRUMlAURVFQyUBRFEVBJQNFURQFlQwURVEUHGiegRDi/9s7uxCrqiiO//5U9lCGmVqDWo5gkT2lEkYfL0WZVCYR+KRQ0FNQD0HGvAg+WdSDFEmRoCH5UpYEkRV9vKSlpqMyfoxmWA7aB5RQTEGrh70Gjpd7Zu517j37OK4fHO5m733n/Fiz7l337nMu+xfgx3H8iWnArx3S6RR1dILwaoc6OkF4tUMdnaAzXjeZ2fRWJl40xWC8SNrd6o8vqqKOThBe7VBHJwivdqijE1TvFctEQRAEQRSDIAiC4NIqBm/mFmhCHZ0gvNqhjk4QXu1QRyeo2OuSuWYQBEEQlHMpfTMIgiAISphwxUDSy5IOS+qXtE3SlMLYi5IGJR2R9GChf4n3DUpa3SWvJyQdkvSfpEWF/jmS/pa0z48NhbGFkg6413p1ePfyMicfyxarBo81kn4uxGfpWI5VUXUsxnA56bmyT9Ju75sq6VNJx/zx2i47bJR0VtLBQl9TByXWe+z6JS2o2CtrXkmaLekLSQP+GnzW+/PFy8wm1AE8AFzu7XXAOm/PB/aTNuPpBY6TNti5zNtzgUk+Z34XvG4FbgG+BBYV+ucAB0ue8y1wJ2mHuI+BhypyyhqrBsc1wPNN+ps6VphnlcdiDJ+TwLSGvpeA1d5ePfJa6KLDvcCCYj6XOQBLPacFLAZ2VeyVNa+AHmCBtycDR/3c2eI14b4ZmNkOMxvZRXwnaZ9lgGXAVjMbNrMfgEHgDj8GzeyEmf0DbPW5nfYaMLMjrc6X1ANcY2bfWMqGzcBjFTlljVWLlDlWRZ1iUcYyYJO3N9Hh/GnEzL4Gfm/RYRmw2RI7gSme81V5lVFJXpnZkJnt9fY5YACYScZ4Tbhi0MCTpGoKKdCnCmM/eV9Zf5X0Svpe0leS7vG+me6Sw6tusXrGvxpvLCx15P6/5T5/IwbskLRH0tPed72ZDUF68wFmZPAqc6hD/GqRV5LmALcDu8gYr/HugZwFSZ8BNzQZ6jOzD31OH2n/5S0jT2sy32heEC/oFqtWvJowBNxoZr9JWgh8IOm2UXyrcOp6rM472SiOwBvAWj/PWuAVUpHvSHzGQe7zN3KXmZ2WNIO0J/nhjC6tkDt+tcgrSVcD7wHPmdmfo1wW7LrXRVkMzOz+0cYlrQIeBu7zJRZIlXR2Ydos4LS3y/o76lXynGFg2Nt7JB0HbnbfWYWpF+R1IU5UEKsirTpKegv4qAXHKsh9/vMws9P+eFbSNtLSxhlJPWY25EsKZzOolTlkjZ+ZnRlp58orSVeQCsEWM3vfu7PFa8ItE0laArwAPGpmfxWGtgMrJF0pqReYR7pA+x0wT1KvpEnACp9ble90SZd5e657nfCviOckLVb6uLASKPsk32lqE6uGddHlwMgdIWWOVZE1b4pIukrS5JE26SaKg+6zyqetorr8KVLmsB1Y6XfJLAb+GFkeqYLceeWv6beBATN7tTCUL16dviKd+yBd8DkF7PNjQ2Gsj3R3wBEKd+aQrtQf9bG+LnktJ1X3YeAM8In3Pw4cIt3BsBd4pPCcRaQkPQ68hv9IsNtOuWPV4PgOcADoJ70gesZyrDDXKo3FKB5zPX/2ey71ef91wOfAMX+c2mWPd0nLnv96Xj1V5kBa9njdY3eAwt1sFXllzSvgbtIyT3/hvWppznjFL5CDIAiCibdMFARBELRPFIMgCIIgikEQBEEQxSAIgiAgikEQBEFAFIMgCIKAKAZBEAQBUQyCIAgC4H9Q65iaEHV8EAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the trajectories (x,y) for each of the 3 masses on the same set of axes\n",
    "# Also create plots of x vs. t for planets M1, M2 and M3\n",
    "\n",
    "plot(fSoln[:,4],fSoln[:,5], 'g')\n",
    "plot(fSoln[:,8], fSoln[:,9])\n",
    "plot(fSoln[:,0],fSoln[:,1], 'or')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mystery System\n",
    "#### Exercise 7.5\n",
    "\n",
    "Let's apply our new 3-body integrator to a more physically interesting problem. Try running your integrator again with the following initial conditions:\n",
    "* $M_{1} = M_{2} = 1000$\n",
    "* $M_{3} = 1$\n",
    "* $r_{1} = $[0.5, 0, 0]\n",
    "* $v_{1} = $[0, $\\sqrt{500}$, 0]\n",
    "* $r_{2} = $[-0.5, 0, 0]\n",
    "* $v_{2} = $[0, -$\\sqrt{500}$, 0]\n",
    "* $r_{3} = $[10,0,0]\n",
    "* $v_{3} = $[0,$\\sqrt{200}$,0]\n",
    "\n",
    "Integrate the equations of motion for a long enough period of time to see multiple orbits of each  \n",
    "object, and plot the trajectories of all three objects on the same set of axes.  \n",
    "\n",
    "***Write a description of the physical system these initial conditions correspond to:***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Defining constants\n",
    "gGrav = 1 ### FIX\n",
    "epsilon = 0.1 #### FIX\n",
    "\n",
    "# Masses\n",
    "m1 = 1000 #### FIX\n",
    "m2 = 1000 #### FIX\n",
    "m3 = 1 #### FIX\n",
    "\n",
    "\n",
    "\n",
    "# UNCOMMENT AND FILL IN initial conditions below\n",
    "#Initial Conditions\n",
    "\n",
    "r1_0 = array([.5,0,0])  \n",
    "drDt1_0 = array([0,500**.5,0])\n",
    "\n",
    "r2_0 = array([-.5,0,0])\n",
    "drDt2_0 = array([0,-(500**.5),0])\n",
    "\n",
    "r3_0 = array([10,0,0])\n",
    "drDt3_0 = array([0,200**.5,0])\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Integrate your equation of motion\n",
    "\n",
    "#tpts = linspace(0, 1000, 100000)\n",
    "\n",
    "fInit = [r1_0[0],r1_0[1],drDt1_0[0],drDt1_0[1],\n",
    "         r2_0[0],r2_0[1],drDt2_0[0],drDt2_0[1],\n",
    "         r3_0[0],r3_0[1],drDt3_0[0],drDt3_0[1]]\n",
    "\n",
    "fSoln = odeint(deriv_3body_2D, fInit, tpts)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x14675240>]"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd4VFX6B/DvmUmvpJEeUggloQQIgdBVIIDuothAsfxcRbHurg0RFREFdde6omJZWQuKAiIIKKGISE2QFiAhhARCCOm9TGbm/P5IUFAIycy999y59/08T54kU855n2H45s65557DOOcghBCifQbRBRBCCFEGBT4hhOgEBT4hhOgEBT4hhOgEBT4hhOgEBT4hhOgEBT4hhOgEBT4hhOgEBT4hhOiEk+gCzhcYGMijo6NFl0EIIQ4lMzOzjHMedLnHqSrwo6OjkZGRIboMQghxKIyxgo48joZ0CCFEJyjwCSFEJyjwCSFEJyjwCSFEJyjwCSFEJyjwCSFEJyjwCSFEJ1Q1D58QqXDO0Wy2oq7ZjLomc+v3ZjPqm81obLGgxWJFi5nDZLGiqcWC+mYLTBYLWiwcnHMYGAMYwMDAGOBsNMDN2QBXJyPcnA1wczLC3eX3n91cjHBzMsLbzQk+7s7wdnWCwcBEvwyEXIACnziMZrMFZ6ubcbqqESW1TSirM+F4aR2yi2uRXVyLumaz6BLbFeHnjl4hPkgI9UZUgCe6ersiyNsVXb1d4efhQn8giOwo8IlqmMxWnKxoQH5ZPfLK6rC3oAp78itQXm8SXZokCisbUVjZiPQjZ9t9XEKoD4bE+qNvuC8i/T0Q6eeBrt6u9AeB2I0CnyiurK4Z2cW1OFpci19yy7DpaInoklTl8JkaHD5Tc9H7wnzdMC4hGMnR/uje1QsxgZ5wczYqXCFxVBT4RDZNLRbknG0N9k1HSrA+q1h0SQ6vqLoJS3YUYMmOC5dOCfRyxfUDwzEk1h99wnzR1cdNUIVEzRjnXHQNv0lOTua0eJpjslg5jpfWYW9BJZbvLcSe/ErRJREANydH4sreXTEwyg9B3q6iyyEyYYxlcs6TL/c4OsInNqltakFmQSV+yCrG0t2nRJdDLuGrjFP4KuP3f58ofw/cntoNw+IC0TPEG0Y6L6ArFPikQyrrTdidX4GVe0/T0IwDO1nRgPnfH7ngtoevise43sFICPOhPwAaR0M65KJMZisyCyrxdcYprPj1tOhyiEI8XYx4ZGw8xvYORkygJxijPwCOoKNDOhT45Dcnyxuw4chZvPZjNupNFtHlEBUYGR+Iu0fGYkiMP80GUjEawyeXxTnHodM1+GxnwQXjvISc8/OxMvx8rAwA4OJkwDNX90ZanxB09aZZQI6IjvB1xmyxYnd+BT7edgLpR2j+O7HdY+N74IZBkQjxpfAXjYZ0yG8459h7sgpvbjyGrTmlosshGvTPcT1wUzKFvygU+AS5JXX4aFseTZskilowpS/+2j8Mnq40YqwUCnydqm5swbI9p/Di2iOXfzAhMuoW4IGXruuL1NgAWgdIZnTSVkfODdk8teIAcs7WiS6HEABAQXkDbv1wFwDgzmHRuH9MHC35IBgFvgOrbmzB57sK8Mr6bNGlENKuT7bn45Pt+fB1d8aiWwdiWFwAzfEXgALfARWU1+OZVVl0ApY4nOrGlt+O+h9P64k7h0XTWL+CaAzfQXDOsTOvAncv2UMXRRFN+Uv/MMy5ujeCabjHZjSGrxEWK8fq/UX4+1f7RJdCiCxW7y/C6v1FiA3yxKJbB6JXiI/okjSLAl+lzBYrvso4hadXHhJdCiGKyCutx4Q3foaHixHL7k1Fn3Bf0SVpDgW+ypjMVny+qwDPrz4suhRChGgwWXDN29sAAKseGI7+kV0EV6QdFPgqYbVyLMs4hVkrDoouhRDVmPzOLwCANQ+NoCN+CVDgC8Y5R/qREtzzPzpZTcilXPP2NkT6u2PpPUMR4echuhyHZZCiEcbYx4yxEsbYofNu82eMbWCMHWv77idFX1qy71QVYp5aS2FPSAecqmjEiJc34/aPd6OqwSS6HIckSeAD+ATAhD/cNgvARs55PICNbb8TACU1TRj/+k+4tu3jKiGk47bmlCJp3ga8tiEHZotVdDkORZLA55xvBVDxh5snA1jS9vMSANdK0ZcjazZbsGDtEaS8tJGWQFAxAwNcjAY4GxmcaA0Y1Xpr4zF0f3odth8vE12Kw5BzDD+Yc34GADjnZxhjXWXsS/V+PlaK2z7aLboMzfF0MaJniDfiu3qjW6AHAr1c4e/hAj9PF/i6O8HT1QkeLk7wdDHCySjVB9o/45yjxcLRaLKgtrkF1Y0tqGpo/apsMOF0VSPySuuQW1KH46X1stWhR7d8sAuBXq5Y+/AIWqvnMoSftGWMzQAwAwCioqIEVyO9qgYTpi7eiaPFtaJLcThhvm4YGheApMguCPN1R4ivG4J93BDg6aK61RcZY3BxYnBxMsDXwxkRHTxj1dRiQWltM4prmnC2pgnHS+qxO78cO/MqYLGq5yp4tSura0bKSxvxeFpPzBwdp7r3h1pItrQCYywawBrOeZ+237MBjGk7ug8FsIVz3rO9NrS0tALnHF9nFuKJbw6ILkXVYgI9MS4hGH3DfREd4ImoAA/4ujuLLks1LFaO4pomFJTX40RZPX7OKcOPh4tBfwva99PjY9AtwFN0GYpRfD38iwT+qwDKOecLGWOzAPhzzp9orw2tBH5ZXTOu/NcW1DSZRZeiGuFd3DE5KQxJkV3QI9gbkf4eMNJRmF1qm1pwvLQe2cU12HS0BD9knRVdkqrcOyoWT0zopYv3maKBzxhbCmAMgEAAZwE8B+BbAMsARAE4CeBGzvkfT+xeQAuBv+ZAER784lfRZQjVK8QbUwaGY2CUHxLCfODhInzkUDc45yisbMSh09XYnF2CZRmFoksSbtuTV2h+7j7teKWwumYzbv1gJ/YXVosuRXF3DY/B2ISu6B/RhZa6VSHOOU5VNGJPfgW+ySzEjrxy0SUpbsGUvpiWor1zhOdQ4Csoq6gaV7+1TXQZivD3dMGMUbEY0T0QvUK8ZZ35QuRT1WBCZkElvttfhFX7ikSXo4hQXzek/3O0Jg9KKPAVwDnH4q15WLDuqOhSZPXwld2R1icEvUN8aPaDRtU0tWDn8XJ8uecUNh0tEV2OrNb/faTmlmCmwJdZo8mCye9s0+QFVGmJwZg+tBsGR/vDzdkouhwiQFFVI7Zkl+K1DTkoq2sWXY7kXrquL24Zop0hHgp8GZ2uasTwhZtElyGpR66Kx7UDwhETqJ+pbKRjTGYrdp+owH9/OYGNGjr6HxLjj8/vHqKJYUkKfJnsPlGBm97fIboMu7kYDXh+ciLGJwQjwMtVdDnEQXDOkX22Fp/vPIlPdxaILkcS+54dhy4eLqLLsAsFvgw+3ZGPZ1ZliS7DLgum9MXV/ULh40YXNxH7ZRfXYsmOfHyx66ToUuzy8xNXINLfcaduUuBLiHOOJ5cfcNg5zU9P6o3rBoYjkI7kiUw45zh8pgb/2ZSLdYeKRZdjk2/uS0VSZBeHHOKhwJeI1cox5d3t2HeqSnQpnXL9wAg8fFV3XV1eTtTBbLHip5xSzFpxEKW1jnfC9+gLExxusgIFvgRMZitSF2xEeb1jbLbg5+GMf9/UH6PigxzyKIVoT2W9CZ/vKsC/fswRXUqnZD2f5lDz9Snw7dTUYkGvZ9aLLqND7hwWjQeu6I4gbxqyIerEOcfuExV47Jv9OFXRKLqcDjn0fBq8HCT0KfDt0Gy2oOcc9Yf9m1OTMKlvKJzpaJ44kJKaJry9KdchZvkcmDveISY4UODbyGS2osecdUJraE/3rl544+Yk9An3FV0KIXYxma1YvrcQT604KLqUdjnC8A4Fvg1aLFbEP63OsB/bOxgvXtcHwbSjD9EYzjm2ZJfi/z7ZI7qUS8qePwGuTuo9kUuB30mcc6Qu2ITimiYh/V/K9KFReHRcT/h5OvaFIYR0xMHCakxdvAP1JovoUv7k+EuTVLu2fkcDnwZ/29z1yR5Vhf1dw2OQ9Xwa5l/bl8Ke6EbfCF9kzZuANQ+NUN0J01GvbIaaDpBtQYEP4PUNOdicXSq6DADA1MGRODB3PJ79S4Lqxw0JkUufcF8cej4Nqx8cAaaSg+rTVY14+ttDosuwi+4TZf2hYry58ZjoMjA+IRgLr+8HfzqaJ+Q3fSN8cWLB1diVV46bF+8UXQ6+2HUSKdH+uHZAuOhSbKLrI/yC8nrc91mm0BqCfVyx7ckrsPj2ZAp7Qi5hSGwATiyYhNdv7i+6FPz9q33ILq4VXYZNdBv4jSYLRr+6RWgN39yXil2zx2p+v01CpMAYw3UDIpAzfyIevKK70FrS3tiKRhWeWL4cXQY+5xyD5m8Q1v+8yYnIe2kSkqP9hdVAiKNycTLgsbSeyJgzFn3Cxe1c1ftZ9V+c+Ue6DPxFW46jQcBf5wFRXbD/2fG4PTWatgokxE6BXq5Y89BILLs3VVgNr/2YLaxvW+gu8PPL6vHqD8r/Iy2fmYqV9w+Hr4f6L9MmxJGkxPgj98WJuH9MnOJ9v7UpFyfLGxTv11a6CnyzxYox/9qiaJ/TUiJx7MWJGNSNhm8IkYuT0YAnJvTCT4+PUbzvUa86zvx8XQX+k8uVXbNj7cMjsWBKP1rcjBCFdAvwRN5Lk/DouB6K9vvENwcU7c9Wukmik+UNWL5XmR2rJieFIWf+RCSEiTuhRIheGQwMD10Vj21PXqFYn19nFuJMtfqXfdZF4HPOMerVzYr0tezeVLw5dQBcnHTx0hKiWhF+Hjj+0iRMS4lSpL/UBZsU6cceukilL/eckr0PZyPD/mfHIyWGxuoJUQujgWHBlL5Yes9QRfpbd/CMIv3YSvOB32Kxyr7e9gNXxCFn/kSagUOISqXGBWD/s+Nl72fm53thtar3BK7mA3+WzCdql89MxeNpvcDUssITIeSifD2ckafAEM/r6erdv1fTgV/fbJb1RO3u2VfRdEtCHIihbYjnnVsGytbH25tyYbZYZWvfHpoOfLkWRovwc0f2/AnoSrtPEeKQru4Xio2Pjpat/e4q3TlPs4Hf1GLBz8fKJG93Wkokfn7iClVvd0YIuby4IC9Zx/V/PVkpW9u20mzgz5Fho4Jnr0nAgin9aLyeEI3w9XBG9vwJCJBhafLrFm2XvE17aTLwOef4JlPasfsPbk/GXSNiJG2TECKeq5MRe54ei5HxgZK3HT3re8nbtIfsgc8Yy2eMHWSM7WOMKbJD+eoD0s6FXXH/MIxLCJa0TUKIehgMDJ/+bQhuSo6QvO13NudK3qatlDrCv4JzntSRXdWl8PDSXyVra81DIzAwyk+y9ggh6vXKDf1x57BoSdsUsTrvpWhuSKe6oUWyttY9MhJ9wn0la48Qon5z/5qI+0ZLu9SyWoZ2lAh8DuBHxlgmY2zGH+9kjM1gjGUwxjJKS0vt7mzOKmlO1q57ZCR6h9LiZ4To0ayJvfB/w6MlbfNUhfh185UI/OGc84EAJgJ4gDE26vw7OeeLOefJnPPkoKAguztbvb/I7jZW3D+Mwp4QnXvuL4m4pl+oZO2NfEWZBRzbI3vgc86L2r6XAFgJIEWuvmqa7B/O+fD2ZBqzJ4QAAN6eNgD9I7tI1t78NYcla8sWsgY+Y8yTMeZ97mcA4wFIP0G+zTub7Dsb/sLkRIyl2TiEkDaMMaycOUyy9j7cdkKytmwh9xF+MIBtjLH9AHYD+J5zLttW7+9vzbP5udcmheG21GjpiiGEaILBwHD0hQmStXflv7dI1lZnyRr4nPM8znn/tq9EzvmLcvZnj9dvThJdAiFEpdycjcicM1aStvJK6yVpxxaamZZZWtts83OPvjCBlksghLQrwMsV3z4wXJK2Hvh8ryTtdJZmAn/J9nybnvfh7clwc6aF0MglfP45EBgIMHbhV2Bg631EV5Iiu+Cx8fZvkP69oJ2xNBP4/7Hx8mU6SUsAtIa3s/Ofg336dKC8/M+PLy9vve/c4wwG+gOgEw9eGQ8vVye729l/qkqCajpHM4Fvi/nX9hFdAhFt7Njfg91str0dzn//A3D//dLVR1QpQ4Lx/Mnv/CJBJZ2j68CfPrSb6BKISEYjsHGj9O2++25r8NMRv2a5ORuxXMLpmkrRdeATHWMMsMq8Dd306XS0r2GDuvlhTE/7VgdQeo0dTQR+rQ1X2J5YMEmGSohDUHJG1rvv0pG+hn10x2DRJXSKJgI/t6Su08+haZg65eGhfJ/TpyvfJ1GE0cCw7N5Uu9rgnEtUzeVpIvB3n6jo1OP7hNPCaLrV2Cim37HSXLRD1Cclxh9Gg+0HkDFPrZWwmvZpIvC35XZus/I1D42UqRKiaomJ4vqW4+QwUY29z4wTXUKHaCLwDxfViC6BOILDYlcqJNrl6+6MyUlhosu4LE0Efnm9SXQJhBCde+WGfjY/t6S2ScJKLk0Tgd8Z22ddKboEQogGuToZce+oWJuem/KiMkN+ugv8sC7uoksghGjU42k9RZfQLk0EvocLLX5GCBHPyWjAX/urdyxfE4E/qBttSUgIUQd7xvLlponAH08rXhJCVMLW5datVvkvwNJE4E8ZGNGhx6243/EWOyKEOJ6Nj47u9HP+a+OeHp2hicD37ODa1AOjaOiHECK/uCCvTj/nhTXyXyeiicAnhBC1UeOFWBT4hBAig3mT1bfBEgU+IYTIwNfdWXQJf0KBTwghMonyF7Acdzs0E/gj4wNFl0AIIRdYfPsg0SVcQDOBTxdfEULUxs1JXasAaCbw30g/JroEQgi5wNPfHhRdwgU0E/iEEKI2v+SWiy7hAhT4hBAigwaTWXQJf0KBTwghMlDjMDMFPiGEyGDx1jzRJfyJJgK/ox+dcktqZa6EEEKAkprOb1l43+g4GSq5kCYCf/2h4g49buxrW2WuhBBCgBvf39Hp5/xjXLwMlVxIE4G/Yu9p0SUQQggAwGLlKChv6PTzXBWYsy974DPGJjDGshljuYyxWXL0sS23TI5mCSGk0/6zKVd0CZcka+AzxowA3gEwEUACgGmMsQQ5+ySEEFE453g9PUd0GZck9xF+CoBcznke59wE4EsAk2Xus12NJovI7gkhGrZ09ynRJbRL7sAPB3D+K1DYdpswvZ9dL7J7QohGcc4xe6VtSymsfXikxNVcnNyBzy5y2wU79TLGZjDGMhhjGaWlpTKXQwgh8rBn7D4hzEfCSi5N7sAvBBB53u8RAIrOfwDnfDHnPJlznhwUFGRTJ7RSJiFEpKYWC/69Qb1j9+fIHfh7AMQzxmIYYy4ApgL4TupOUmMDOvX4RVvUexadyKhLF9EVEI2a8Wmm6BI6RNbA55ybATwI4AcARwAs45xnSd3P4Bj/Tj3+lfXZUpdAHEFlpbi+E2hymlYVlNdja47tw9E58ydKWE37nOTugHO+FsBaOftIiqAjN6JyWZIf5xAV4Jxj9Ktb7GrDxUm56181caWtr0fnNwse/GK6DJUQ1eP88o/RQp9EEYu2HBddQqdoIvBtUVrbLLoEIoqSATxzpnJ9EUUVVTXi1R/sGx4+/tIkiarpGN0GPtA69kZ0SonQnzkTWLRI/n6I4ixWjmELN9ndjtFwsZnr8tF14Ns79kYcHOfAZ59J325YWGvbFPaa9ciXv9rdxpyre0tQSedoJvCnpUTZ9Dxb1q0mGnLrra3hzDng7m5fW1dd1drOaVq9Vcu2ZJdgzYEzdrdz98hYCarpHM0E/kwbNw9IeWkjrFY6qUYANDT8Hv6ffQYEtHN9h5dX62POPZ5zIJ0mAmjdmepG3PnfPaLLsJlmAj8qwMPm517/3nYJKyGacOutQFnZhYF+/ldtbetjiG40tViQusD+cXsAyF94tSTtdJZmAt8ev56swie/nBBdBiFEpTjn6PPcD6LLsJumAn9yUpjNz527+jB+PkaLtxFC/uz2j3fDLNHQr6ije0Bjgf/MNfZdvn7bR7tx7CxtdE4I+d3L64/i52Pa2FVPU4Ef6OVqdxvjXt+KoqpGCaohhDi6z3cV4F0Jr6YVeXQPaCzwAWmWSh62cBNN1yRE55ZnFuLplYcka0/EvPs/0lzgv3PLQEnaSXlpI8rqaPkFQvRozYEiPPr1fknbFDHv/o80F/ghvm6StZU8P51CnxCdWb2/CA9+Yf+VtOfLU3jNnEvRXOAD0n50Sp6fjtM0pk+ILny6swAPLZU27AHAoPCaOZeiycC/c1i0pO0NX7gJuSU0e4cQLXsz/Rie+Va6MftzRJ+oPZ8mA9/JaMDgaGn3uR372lZkFlRI2iYhRDzOOf65bB9eT5d+T1q1DOWco8nAB4CP7hwseZvXv7sDyzMLJW+XECJGs9mCES9vxoq90i9499f+YaoZyjlHs4Hv4+aMAE8Xydt99Ov9mPPtQXDaxYgQh1ZRb0LPOetlO0f31rQBsrRrD80GPgD8+I9RsrT72c6TGPHyZjSbLbK0TwiRV87ZWgx8YYNs7atp3P58mg78AC9X9A71kaXt01WNsh4dEELk8cWukxj/+lbZ2n//tkGytW0vTQc+AHx9X6qs7Q9fuAk/ZhXL2gchxH4msxW3fLATs1celLWftMQQWdu3h+YD38vVCXePiJG1jxmfZuKRL3+FhTZSIUSViqoa0WPOOmw/Xi5rP5sfGyNr+/bSfOADwOxJ8q9hsWpfEeJmr0VhZYPsfRFCOu7L3Scl2XD8ckJ83BAT6Cl7P/bQReAbDEz2oZ1zRry8GZ/uyKdZPIQIVt3YgtGvbsasFfIO4Zyz5fExivRjD10EPgAMjvZHrEJ/fZ9ZlYWkeRtQ1WBSpD9CyIW2HStD/+d/REG5Mp+4379tENycjYr0ZQ/dBD4ArH1kpGJ9VTe2IGneBnyTWUhH+4QopLapBVMW/YLpH+1StF81n6g9n64C383ZiKX3DFW0z8e+3o8ec9bhTDVN3yRETt8fOIO+c3/E3pNVivZ7eF6aov3ZQ1eBDwCpcQEYlxCsaJ8tFo7UBZvwRnoOzeQhRGIlNU1Inp+OB77Yq3jfy2emwsPFSfF+baW7wAeA96aLuTDijfRjiJu9Frvy5J0aRogemMxWzF9zWNhmRTcMisCgbv6K92sPXQa+0cCQMWessP5vXrwTI17ehOJq2kaREFv8mFWMHnPW4cNtJ4TV8OoN/YT1bStdBj7QuuH553cPEdZ/YWUjhi7YiKdWHESDySysDkIcSW5JLaJnfY8Zn2YKrePwvDQwpq6VMDtCt4EPAMO7B+KJCT2F1rB090kkPPsD3t1yHC0Wq9BaCFGr01WNmPDGVox9Tb41cDrqp8fHONS4/flkC3zG2FzG2GnG2L62L3XtBNDm/jHdMSwuQHQZeHn9UcQ/vQ7LaRonIb8pr2vG3UsyMHzhJhwtFr/r3Ed3JKNbgLqvpm2P3Ef4r3POk9q+1srcl80++9sQqGWfgke/3o+Yp9biu/1FsNKMHqJTlfUm/HPZPgyan470I2dFlwMAmD2pF67qrewMP6npekjnHIOB4fC8CaLLuMDDS39F7Oy1+CazkKZyEt0oqW3CzM8yMeCFDbLsQmWraSmRmDEqTnQZdpM78B9kjB1gjH3MGJN2k1mJuTkbkfW8+i6geOzr/YibvRZLtufThitEs05XNeK2j3Yh5cWNWHdIXcuNp8YGYMEUx5uRczHMnvFixlg6gItdU/w0gJ0AygBwAC8ACOWc33WRNmYAmAEAUVFRgwoKCmyuRwpVDSYkzZNvJxx7/d/waDx8ZTz8ZNi+kRCl7T1ZiZmfZeJsjfLz6Duid6gP1j48QvUzchhjmZzz5Ms+TokThIyxaABrOOd92ntccnIyz8jIkL2eyymtbcbgF9NFl9GulGh/LLi+L+KCvESXQkinmC1WrDlwBn//ap/oUtrVK8Qb6x4ZqfqwB1QQ+IyxUM75mbaf/wFgCOd8anvPUUvgA60njQbIuOellP5zywCkJYbA2UinZIh6na1pwqLNuViyQ+yn+I7oG+6L7x4c7hBhD6gj8D8FkITWIZ18APee+wNwKWoKfACoaWpBv7k/ii6jw24cFIF/jOuBsC7uokshBABgsXL8lFOCR5buQ22zY1xgOLZ3V3xwe7LDhD2ggsC3hdoCHwDqm81IfO4H0WV02r9u7I9JfUMc9gIR4tgKyuvx/tY8fLHrpOhSOmXGqFhFdsiTGgW+hExmKya+uRXHS+tFl9JpvUK88ew1CRgaGwCDWi42IJpUUW/C8sxCvLj2iOhSbPLCtX1w29BuosuwCQW+xDjnePTr/aqaG9xZV/cNxT2jYtE/wtehPq4S9aprNiP98FnMXnkQDSbHnTa85K4UjO4RJLoMm1Hgy+R/O/Lx7Kos0WXYLS0xGDNGxWJApB8d+ZNOqWowYf2hYsxdnYWmFsdf/+mnx8c49HIJAAW+rDILKnD9uztElyGZ1NgAzBgdi9TYAIfYl5Mor7i6CWsPnsG8NYdFlyKpg3PHw9vNWXQZdqPAl5kjzNW3hYuTAU9O6IW0xGBE+HmILocIYjJbkVlQia/2nMS3+4pElyO5tMRgLLp1EIwa+XRLga+AFosV9/wvA1uyS0WXIpuxvYMxdXAkUmL94aOBIyFycZxznKxowE85pXgj/Rgq6k2iS5LNO7cMxNX9QkWXISkKfAWt2ncaj3yp7qsGpXJtUhiuGxiBwdF+NOXTgXHOkV/egF9yy/DRthM4UeZ4M9BssX3WlZq8ToUCX2FnqhuRumCT6DIUN7Z3V/ylfxgGRvkhws+dZv+oVLPZgiNnarH7RDmWbC/A6apG0SUpamKfELw1bYBmr0anwBfAYuWYtzrLIS4dl4u7sxF3DIvG8O4BSAzzhT8t8qY4zjkKyhuw71QVvj94BhsOq2M9eVG+uGcIhsUFii5DVhT4Ah06XY1r3t4mugzVcDEacP2gcIyMD0JimA8i/TxoKqhE6prNyC6uxeGiaqzPKsYvueWiS1KNfhG+WHrPUHi6an/okQJfsBaLFc+uOoSlu0+JLkW1ovw9MLZ3MAZEdUFskCdiA73g7kLTQv+Ic46zNc3IK6tDXmk9fsktQ/qRs2ixqOf/rtp8eHsyxiY49u5UnUGBrxIlqV5eAAAKTElEQVRHi2sw4Y2fRZfhcHqFeCMlxv+3TwQRfh4I8XWDi5P2xmA556ioN+FMdRNOVzXiZHkD9p6sxO4TFSjX8GwZOYyMD8R70wfp4qj+fBT4KmKxcrz303G8+kO26FI0pUewF3qF+KBXqDfCu7jDz8MFvu7O8HF3ho+bE3zcnYWcpOOco6nFitrmFtQ2mVFRb0J5XTNK60woKKtH9tlaHC2uRWmtOjf9cFRrHhqBPuG+ossQggJfhcrqmjH9w104WlwruhTd83AxwsfNGT7uTvB2c4aHixGuTgYwxnD+2QXGAM4Bk8WKphYLapvMbV8tqGs207CKCjye1hP3jY7TzEVUtuho4Ovrc49ggV6uWP/3UcjIr8AN72lnaQZH1GCyoMFkQXGN6EqIrYbG+uO96YPQxYNmgnUUBb4AydH+yHtpEr7OPIUnlx8UXQ4hDmfjo6Npe08bUOALYjAw3Dw4CpOTwvHahhws3ponuiRCVM/RlzEWTXtTHhyMm7MRsyf1xq/PjMPkpDDR5RCiSm9PG4ATCyZR2NuJTtqqTFldM5759hDWHSoWXQohwr1yQz/cMDCCLtS7DDpp66ACvVzx7vRBKK1txnPfHcLagxT8RH9eub4fpgwMh5NG174RhQJfpYK8XbHo1kGorDfhjfQcXa/PQ/Tj/dsGYXxCMC3CJxMa0nEQTS0WfLI9HwvXHRVdCiGScnEy4H93pWBIjD8FvY1oSEdj3JyNuG90HO4ZGYsNh4tx32d7RZdEiF1G9wjCvMmJDr+frCOhwHcwRgPDhD6hyF94NXJL6vDCmsP4KUe7O24R7XlyQi/cntpNd+vdqAG94g6se1cvLLkrBQ0mM5btOYW5q7W1wTTRjih/D7x6Qz+k0LCNUBT4GuDh4oQ7h8fgzuExOFxUg7c2HsP6LJrdQ8R7amIvTE2Jgq877YesBhT4GpMQ5oP3bhuEFosVm4+WYPbKQyiro1UZiXIm9Q3Bg1fEIyHMR3Qp5A8o8DXK2WjA+MQQjE8MQWW9Cd8fPIM53x4SXRbRqMHRfng8rReSu/nRRVIqRoGvA36eLpg+tBumD+2GinoTVu8vwnPfZYkuizi4fhG++Me4HhgeF6jJjWm0iAJfZ/w9XXDHsGjcMSwalfUmbDh8Fq9tyEFxTZPo0ogDGJ8QjHtGxWJglJ+u1593VBT4Oubn6YKbBkfipsGRMJmt2JNfgQ9/zsPmbJrmSX734BXdcU3/UPQM9qYZNg6OAp8AaL3acXj3QAzvHggAOF5ah81HS/DO5lxUNrQIro4oaVSPINyR2g2pcQHwcKGI0BL61yQXFRfkhbggL9w9MhYWK0dWUTXWHyrGoi3HRZdGJNY71Ad3DY/GsO6BCO/iLrocIiO7Ap8xdiOAuQB6A0jhnGecd99TAP4GwALgYc75D/b0RcQxGhj6RXRBv4gueGJCL5jMVmQVVWNHXjk++SUfJbQZt0MZGR+Im5IjMSTGH1193ESXQxRk7xH+IQBTALx//o2MsQQAUwEkAggDkM4Y68E5t9jZH1EBFycDBkT5YUCUH+4f0x0AUFTViMyCSqw9eIbW8leREB83TEuJwrDuAUgM86EhGp2z61+fc34EwMVO5EwG8CXnvBnACcZYLoAUALRzt0aFdXFHWBd3/KV/665dVitHfnk9sopqkH7kLFbtKxJcofZF+XvgugHhGBDV+mnM35M29yYXkuvPfTiAnef9Xth2G9EJg4EhNsgLsUFe+Ev/MLw5dQA456ioNyG3pA45JXX4KbsE6UdKRJfqcLp6u2J8YjCGxASgV4g3ogM94UwbhZAOuGzgM8bSAYRc5K6nOeerLvW0i9x20YX3GWMzAMwAgKioqMuVQxwYYwwBXq4I8HLFkNgA3Da022/3NZjMOFFWj8LKRpwoq8eO4+X4JbcMZqt69mtQUoiPG0bEByK5mx+iAz0RE+iJrt6uNC2S2OWygc85H2tDu4UAIs/7PQLART/Tc84XA1gMtG6AYkNfRAM8XJyQGOaLxDBfAMB9o+N+u49zjsqGFpytaUJJbTNKappwoqweB09X42hxLUod7KRxbJAnEkJ90DvUBxF+rUNhob5uCPZxoyN1Iiu5hnS+A/AFY+w1tJ60jQewW6a+iMYxxuDv6QJ/Txf0Dm3/sZxzNJgsqG5sueCrqcWCRpMFDSYLKhtMKKszobapBXXNZjQ0W2CyWNFiscJkbv3eYuFgDDAwBhcnA9ycDXBzMsLN2Qg3ZwO8XJ3g6+4MXw8X+Hk4w8vVCd5uTvB0dYJX25enqxO6eDjTiVKiGvZOy7wOwNsAggB8zxjbxzlP45xnMcaWATgMwAzgAZqhQ5TAGINnW9iG0ZxyQi5g7yydlQBWXuK+FwG8aE/7hBBCpEMDhoQQohMU+IQQohMU+IQQohMU+IQQohMU+IQQohMU+IQQohMU+IQQohOMc/WsZsAYKwVQIFPzgQDKZGpbC+j1uTx6jdpHr0/75Hx9unHOgy73IFUFvpwYYxmc82TRdagVvT6XR69R++j1aZ8aXh8a0iGEEJ2gwCeEEJ3QU+AvFl2AytHrc3n0GrWPXp/2CX99dDOGTwgheqenI3xCCNE1TQc+Y+xGxlgWY8zKGEv+w31PMcZyGWPZjLE0UTWqCWNsLmPsNGNsX9vXJNE1qQFjbELb+ySXMTZLdD1qxBjLZ4wdbHvfZIiuRzTG2MeMsRLG2KHzbvNnjG1gjB1r++6ndF2aDnwAhwBMAbD1/BsZYwkApgJIBDABwCLGmFH58lTpdc55UtvXWtHFiNb2vngHwEQACQCmtb1/yJ9d0fa+oamZwCdozZbzzQKwkXMeD2Bj2++K0nTgc86PcM6zL3LXZABfcs6bOecnAOQCSFG2OuIgUgDkcs7zOOcmAF+i9f1DyCVxzrcCqPjDzZMBLGn7eQmAaxUtChoP/HaEAzh13u+FbbcR4EHG2IG2j6SKf+RUIXqvdAwH8CNjLJMxNkN0MSoVzDk/AwBt37sqXYDD767MGEsHEHKRu57mnK+61NMucpsupiu193oBeBfAC2h9LV4A8G8AdylXnSrp9r3SScM550WMsa4ANjDGjrYd5RIVcfjA55yPteFphQAiz/s9AkCRNBWpW0dfL8bYBwDWyFyOI9Dte6UzOOdFbd9LGGMr0ToURoF/obOMsVDO+RnGWCiAEqUL0OuQzncApjLGXBljMQDiAewWXJNwbW/Cc65D60lvvdsDIJ4xFsMYc0Hryf7vBNekKowxT8aY97mfAYwHvXcu5jsAd7T9fAeAS41AyMbhj/Dbwxi7DsDbAIIAfM8Y28c5T+OcZzHGlgE4DMAM4AHOuUVkrSrxCmMsCa1DFvkA7hVbjnicczNj7EEAPwAwAviYc54luCy1CQawkjEGtGbKF5zz9WJLEosxthTAGACBjLFCAM8BWAhgGWPsbwBOArhR8broSltCCNEHvQ7pEEKI7lDgE0KITlDgE0KITlDgE0KITlDgE0KITlDgE0KITlDgE0KITlDgE0KITvw/gf/kbd/TU/UAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plots go here\n",
    "plot(fSoln[:,4],fSoln[:,5], 'g')\n",
    "plot(fSoln[:,8], fSoln[:,9])\n",
    "plot(fSoln[:,0],fSoln[:,1], 'or')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting in the center-of-mass frame of reference\n",
    "\n",
    "You have probably noticed an interesting precession effect.  While we start the three masses in orbits that should take place (approximately) centered on the origin, their mutual gravitational interactions slowly carry them off in one direction over time -- the exact direction will depend on exactly how you chose your initial conditions. Essentially, because of the initial conditions, the system has some non-zero momentum, which causes the center of mass to move.\n",
    "\n",
    "We can make our plots look much nicer by transforming to a coordinate system where the center of mass remains fixed at the origin.  \n",
    "\n",
    "#### Exercise 7.6\n",
    "Define two new arrays, `x_com` and `y_com` to track the location of the center of mass of the 3-body system, and then plot the trajectories of the three masses on the same set of axes.  If you have done this correctly, you should find that the precession effect has disappeared."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x148632e8>]"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHmZJREFUeJzt3Xd4FVXeB/DvLwESepHQS0BBBAR0IxZsFGlRUVddZV9d24vu6u7ququhSBWIbV1ddXdxRV17e62hCArSpASlhQ4JxSAkKig95bx/kKwBUs7cO+dOOd/P8+SBJGdmfs/N5Jtzz5w5I0opEBFR+MV5XQAREcUGA5+IyBIMfCIiSzDwiYgswcAnIrIEA5+IyBIMfCIiSzDwiYgswcAnIrJENa8LKKtx48YqOTnZ6zKIiAJl+fLl+UqppKra+Srwk5OTkZmZ6XUZRESBIiLbdNpxSIeIyBIMfCIiSzDwiYgswcAnIrIEA5+IyBIMfCIiSzDwiYgs4at5+EQmHS0sxpLs7/D8/GzM25hn/HjN6yfi93064IruzVE3sbrx4xFVRfz0TNuUlBTFG6/IiR3fH8TVzy1C/v4jXpfiitfuOBe9TmvsdRkUMCKyXCmVUlU7V3r4IjIVwOUA9iilupZ8rRGAtwAkA8gBcL1S6gc3jkf2OHi0EFc9uxAbd+/3upSY+PW/l1T6/fkP9EbrRrViVA2FjSs9fBG5GMB+AP8pE/iPAvheKZUuImkAGiqlHqxsP+zh2yt37yFckP6512UE1ut3nIsL+M7AWjHt4Sul5olI8glfHgLg0pL/vwxgLoBKA5/s8N7ynbj/nZVelxEqQ8t5Z3BVjxZ48lc9ICIeVER+ZPKibVOl1C4AUErtEpEmBo9FPrXzh4O48JE5XpdhpQ9W5OKDFbnHfe2F36Sg7xlNPaqIvOb5LB0RGQZgGAC0adPG42ooWh98/Q3ufWuF12VQBW5/+fgh08u7NcczQ8/2qBqKNZOBv1tEmpf07psD2FNeI6XUFABTgGNj+AbrIQMWbcnH0Ocrv9DoRzXi43Dbhe0wqGszdG5RD9Xjo78lpbCoGBt370fG6lw8O2eLC1Wa98mqXfhkVcZ/P7/zkvYYPugMDysik1ybllkyhv9JmYu2jwH4rsxF20ZKqQcq2wcv2vrf0cJidBw13esyTtKsXiLe/e35aNUweDNY9h48irteXY7FW7/3upSTLEzrg5YNanpdBlVB96KtW7N03sCxC7SNAewGMAbABwDeBtAGwHYA1ymlKj2jGfj+tHBzfpXTBWPhi79ciran1Pa6jJjbd7AAKRNnoaDI2zfAv0ppjUeu7eZpDVS+mAa+Wxj4/jHu4yy8uDDHk2OvGz8QNWvEe3LsICkuVrhp6hIs3PydJ8fPSU/15Lh0MgY+OZb23iq8uWxHzI43qGszPPfrszlt0GUfrvgGf3wzthfOGf7eYuCTlve/3on73orNnPiscQNQO8HziWHWUUrhzLGfYv+RQuPHOqtNA7z/u17Gj0PHY+BThfYePIoe42cZP0725MHsvfvU3a99hYzVu4we49Xbz8WFHXj3byww8OkkyWkZVTeKwvoJA5FYnWPvQaOUQrvh04weY+ukwYiL4x9/Uxj4BAD4Zu8h9DK0Rs2kq8/E0HN5s1zYbNz9E/o/Oc/Ivh+9thuuT2ltZN82Y+Bb7prnFuKr7Xtd3++ykf2QVDfB9f2SPxUVK5w6wkzvnxd63cPAt5Cpt+YcqiGA55efMfAtcuhoEc4YPcPVfa4a2x/1+JQmqoCJ8P/k9xeia8v6ru7TFgx8C7g9Pv/msPNwXvtTXNsf2eGHA0dx1gT3Zn1xnN85Bn6Ibfj2Jwz4m3sX1TiWSm55dMZ6PDfXnYXj7u3XAff26+jKvsKOgR9CW/P2o88TX7iyr6Uj+6JJ3URX9kV0Ijcv9o5KPQN3XNTelX2FFQM/RNy8UYq9eYq1m15Ygvmb8qPez8u39cQlHZNcqCh8GPghUFys0N6FXtJ9/Trij/06uFARUeTcGopcPLwvmtXnu9OyGPgB58ZdsbP/dDFOa1LXhWqI3FNQVIwOI6N/pgLv3v0ZAz+ghjy7ECt3RHfD1Oqx/VGXUyrJ59ya2slhSgZ+4Gz/7iAufiy6h31veHggEqrxBhYKnmjf0Y4Y3AnDLj7VpWqCh4EfINGe7LxTkcIi2t+FzRMHoZoLzycOGgZ+AIz9KAsvLcqJeHtOraQwcmOox7ZhHga+j0V7Qj9/cwou69zUxYqI/OfAkUJ0GTMz4u2XjuiLJvXs6BAx8H3qtpeW4fP1eyLatlXDmljwYB+XKyLytznr9+DWl5ZFvL0NvX0Gvs9E26u34aQlqkw04/srRl+GBrVquFiNvzDwfeSVxdvw0AdrItr2q4cuQ6Pa4T1RiZxgx6l8uoFv3+XsGEtOy4go7FvUT0ROeirDnqgMEUFOeipeub1nRNubfsyn3zHwDTlcUBTxyZU9eTAWDe/rckVE4XFRh6SIe+vJaRlY4MLaPkHEIR0DJk1bhynztjre7saebTD5mjMNVEQUXtn5B9D78bkRbRuWIR6O4Xskml69CNcFIYpUpL97YQh9juF7IJITrkndBOSkpzLsiaKUk56KqbdUmXknSU7LwP4jhQYq8h8Gvgv2HSqIKOzXjBuApSP7GaiIyE59OjWNqMfedcxMvPJljuv1+A0DP0ozs75F93GfOt4uJz0VdRKqGaiIiCIJ/Yc+zAr9LB4GfhSu/+eXuPOV5Y626dOpSSjGDIn8Lic9FTPuvcjxdmEOfQZ+hJLTMrA053tH22SO6oept5xjqCIiOlGnZvWQPXmw4+3CGvoM/AhEcjLkpKeicZ0EA9UQUWVKb9ZyKoyhz8B3KNKwJyJv5aSnIrG6s8gLW+gz8B1g2BMF2/oJg/DEdd0dbZOclgE/3a8UDQa+Jqdhf+cl7Rn2RD70y1+0wpfDnS0z3m74tFCEPgNfg9Owf3PYeRg+6AxD1RBRtJrXr4lNEwc52saNB657jYFfBadhv3h4X5zX/hRD1RCRW6rHxzmewRP0MX0GfiWc/nBXj+2PZvXteKQaURhEMoMnyKFvPPBFJEdEVovIChEJzMpoTn+o6ycMRN3E6oaqISKTbAn9WPXweyuleuis5uYHA56c56j9uvEDkVg93lA1RBQLNoQ+h3ROMH31LmzY/ZN2+5Vj+qNmDYY9URiEPfRjEfgKwKcislxEhsXgeBHbd7AAv33tK+32c/58KerX5DAOUZg4vZA7/uO1hipxXywCv5dS6mwAgwDcLSIXl/2miAwTkUwRyczLy4tBORXrPl5/1ct//s/ZaNe4tsFqiMgLIoINDw/Ubj91YTbyfjpisCL3GA98pVRuyb97ALwPoOcJ35+ilEpRSqUkJSWZLqdCTt6aXda5KQZ2bW6wGiLyUkK1eMz586Xa7c+ZONtcMS4yGvgiUltE6pb+H0B/AGtMHjMSTsfhnr85ENeeiSgK7RrXxv2XddRuH4TxfNM9/KYAFojISgBLAWQopWYYPqYjczbscdSeyyUQ2eP3fTs4au/30Dca+EqprUqp7iUfXZRSE00eLxK3vrhMuy3Dnsg+Tn/vs/MPGKokelZPy3Ty13j9BP2LOEQULk5Cv/fjc80VEiVrA//6f36p3fYPfTvwxioiy60ZN0C7rV+HdqwM/OJi5ejxhH9ycOGGiMKpTkI19O3URLv9Fxu9nWZeHisDv/0I/WVOOW5PRKVecPBM6t9MXWqwkshYF/iPzliv3XbJiL4GKyGiIHLSCfTb0I51gf/c3C3abZvW41LHRHSyp27ood1238ECg5U4Y1XgO/lry6EcIqrIkB4ttds6WbLFNGsCv6CoWLvt5/dfYrASIgoDJ53CKfP0RxZMsibwO4ycrt22fVIdg5UQUVh0b91Aq92kafrXDk2yIvB3/nBQuy2HcohI14d399Ju64cbsqwI/AsfmeN1CUQUUtP/eJFWOz8suRD6wN/2nf6LzN49ETl1RvN62m29nqYZ+sC/5LG5Wu2u/UUrs4UQUWhlOVh2wUuhDvz9Rwq12z5+XXeDlRBRmNVOqKbd9p7X9R+j6rZQB37XMTO12t18flvDlRBR2OmuqPvJql2GK6lYqANf1/ghXb0ugYgCzsmKulm5+wxWUrHQBr7XF0eIyD4LHuyt1S716QWGKylfaANfF2fmEJFbWjWs5XUJlQpl4H+z95DXJRCRpS49PUmrnRejEKEM/F7pn2u1Wzs+GFOpiCg4Xrq1p9clVCiUga+rVg39qVRERG4rKlYxPV7oAn/Vzr1el0BElpv/gN7F21MdPH3PDaEL/CufWajVjhdriciU1o38efE2dIFPRETlC1Xgx3o8jIioIrP/pPcgpd+/8bXhSn4WqsDXHQ/bMmmw4UqIyHanNdF7kNLHK3MNV/KzUAW+rvg48boEIqKYszLwiYhsZF3g677NIiKK1saHB2m1+/f8rYYrOSY0gf/K4m1a7T6992LDlRARHVOjml7EPpyxznAlx4Qm8B/6YI1WuziO3xORpUIT+EREVDkGPhGRJRj4RCfY0KkxIKL3QVSFRWl9tNopZf7GUasC/5Ffnul1CeRThw7v/2+In77hO/0NS7bZcFoDc8VRoLVoUFOr3efr9xiuJCSBr7ukwvUprQ1XQoEkgpo160a1i9O37ANEoIqLXSqKbHP7y5nGjxGKwP9io95fRuFbcCrjcMEh14dlJD4exTzPyKdCEfiPz9zodQkUQIk1zCxhGwfgb/edb2TfRNEwHvgiMlBENojIZhFJM3GMtbt+NLFbCjPDvfB7/7bY6P6JImE08EUkHsCzAAYB6AzgRhHpbPKYRFXJWqb5lCGlKv7QwaEd8hnTPfyeADYrpbYqpY4CeBPAEMPHJKpUl54aTzurKtSVAniBlgLGdOC3BLCjzOc7S75G5F/swVNImQ788n4jjvttEpFhIpIpIpl5eXmGyyEispfpwN8JoOzk91YAjnu8i1JqilIqRSmVkpSUZLgcoqrtPeDgxiuiADEd+MsAdBCRdiJSA8ANAD4yfEyiSv3Yqkml329Qp3GVwzUvjL6CQzoUONVM7lwpVSgi9wCYCSAewFSlVJbJYxJVpd6O3XphXUmb23UOFIO1UYicMBr4AKCUmgZAcx5cZFrUT0TuvsMmD0EhU5S9FfHt2ntdBlFMheJO25GpnNpPzsQnt8M3yaeYOwB79+RDoQj8fp0rH5MtFYvlRyk4WmbnQxUVubrPXZ1aMezJt0IR+AnV4rXaTV/zreFKKGgkLs69gFYKzdftqLodUTmeuqGH8WOEIvB1/e61r7wugfzK6bIJJQ7/kBfRdmSPwwV67yKv6NbCcCUxuGhLFDgOwjvRYBkUDjc+r7eQXlyc+Wm+VvXwiYhi7evte70u4b8Y+ERElghN4N90XlutdpypQ0S2Ck3gjx/SRatd2nurDVdCRHSMbgfzlguSzRZSIjSBr/u82rcyOW2OiGLjl/9YpNVuzBWxuXk0NIFPROQ3X2lesNXtsEaLgU9EZIlQBX7WuAFa7VIenmW4EiKy3cGjhVrtWjaoabiSn4Uq8Gsn6N1Hlr//qOFKiMh2nUfP1Gq34MHehiv5WagCn4goaGI1fg+EMPAfu7abVrs7X8k0XAkR2epIobursLoldIF/XUrrqhsBmJm123AlRGSr00fN0Gq3amx/w5UcL3SB7wTvuiUiL9VLrB7T44Uy8F+85Rytdu2GG33yIhFZaMaaXV6XUKFQBn7vTnpPwCIicttdr+o9d2PrpMGGKzlZKAPfiUnT1nldAhGFhJOLtbFY//6kY8b8iDGi+9dzyrythishIlvoXqyddPWZhispX2gD38lfz6XZ3xushIhs4GQSyNBz2xispGKhDXwA+Pz+S7TaXf+vLw1XQkRhF4RJIKEO/PZJdbTbLtqSb7ASIgozJ737nPRUg5VULtSBDwBP33iWVruhzy8xXAkRhVUQeveABYF/ZfcW2m0fm7neYCVEFEb7j+itigl427sHLAh8APj3zSla7Z6ds8VwJUQUNl3H6K2K6QdWBH6/zk212yanZRishIjC5LUl27Tbet27BywJfABYPLyvdtujhcUGKyGisBj5/hqtds3rJxquRI81gd/MwQvecdR0g5UQURg4GQ340kGH0yRrAh8Asifrr11xzXMLDVZCREGWnX9Au+07d51vsBJnrAp8EUG9RL3HIH61fS+Ki7l8MhGdrPfjc7XbnpPcyFwhDlkV+ACwaqzeg84BoP2IYMytJaLYcTKU42RUIRasC3wAWDla/ykznLVDRKXeW75Tu+2V3VvE9Hm1OqwM/Pq1nD1lJit3n6FKiCgolFK4/52V2u117/KPJSsDH3A2Jzb16QUGKyGiIHCyfILfhnJKWRv4ALDx4UHabTm0Q2QvJ7//44d08d1QTiljgS8iY0XkGxFZUfLhuz95NarF4fJuzbXbM/SJ7PO//8l01P7m85PNFOIC0z38J5VSPUo+fDnl5ZmhZztqP/ip+YYqISK/2bznJ8xau1u7vR+WT6iM1UM6pZz8kNbu+hGzHZwARBRMBUXF6PfXedrtN0/UHyL2iunAv0dEVonIVBFpaPhYUdnk4Id1x38ysefHwwarISKvdRipv8TK1FtSUC3e//3nqCoUkdkisqacjyEA/gHgVAA9AOwC8EQF+xgmIpkikpmXlxdNOVGpHh+H1+44V7t9z0mfoaCIi6wRhZHT63V9OumvyOslcfJorogPIpIM4BOlVNfK2qWkpKjMTGcXSNw2+Kn5WLvrR+322ZMH+/aKPBE55zTs/TBuLyLLlVJVPvjD5CydstNfrgagt46ox6b98SJH7dsNn+boeZZE5F9BDHsnTA46PSoiq0VkFYDeAO4zeCxXOf0hMvSJgi/sYQ8YDHyl1E1KqTOVUt2UUlcqpXaZOpYJDH0iezgNe7/eSVsV/19W9lAkoc8llYmCxWnYb5o4KLDX7Rj4VXAa+u1HTMPBo/pPsSci7zgN+3XjB6J6AKZfViS4lceQ09DvPHomtuTtN1QNEUVLKeU47FeP7Y+aNeINVRQbDHxNTkO/7xNfYMq8LYaqIaJI/XDgqKOVLwFg1dj+qJvobFl1P2LgO+A09CdNW88F14h85KOVuThrwixH26wbPxD1QhD2AAPfsUimYjH0ibzXfngG/vDG14622TJpcOCHccpi4Ecg0tDntE0ibySnZcDpBLqc9FTExwVzNk5FGPgRiiT02w2fxsclEsVQUbHzi7NAMG+q0sHAj0JOeioSqzt7CVOfXsAhHqIYeG/5Tpw6wvljOMIa9gADP2rrJwzCUzf0cLwdQ5/InOS0DEcPHC8V5rAHGPiuGNKjJbLGDXC8XXJaBuZv8m5JaKKwOVxQFFFn6qkbeoQ+7IEYLY+syw/LI0cr0p67DScbkUm90j/HN3sPOd5uw8MDkVAt2DNxPF8e2VaRBndyWga2fXfA5WqIwq/0rtlIwj4nPTXwYe8Ee/iGzNuYh5unLo1oW/b2ifT0fnwusvOdd5Sa1kvAkhH9DFTkDfbwPXZxx6SIl1BNTsvA0uzvXa6IKDxKx+ojCfulI/uGKuydYA8/BqKZkcPePtHx+Pt0MvbwfSQnPRXLR0XWo0hOy0DXMTNdrogoeN7J3BFx2I+9onNow94J9vBjLJreyXu/vQC/aNvQxWqI/O9wQRE6PTQj4u1tCHrdHj4D3wPZ+QfQ+/G5EW8fhmlkRDqi6SD97tJT8cDATi5W418M/ACI9m7b7MmDA/uoNaLKRPu7YUOvviwGfkBE+3YVsO/kpvA6Z+Js5P10JOLt3/jf83D+qae4WFEwMPADZvj/rcYbS7dHtQ8GPwXVVc8uxIode6Pah83nPwM/oNxYVM3mE5+C5YLJnyF33+Go9sGhTQZ+oBUUFaPDyOlR74e/CORXbnRsZtx7ETo1q+dCNcHHwA+Btbk/YvDT86Pez8ox/VG/ZjieyUnBpZRy/PDw8txyQTLGXtnFhYrCg4EfIlMXZGP8J2uj3s+Eq7ripvPaulARkb5d+w7h/Mmfu7IvDleWj4EfQr99dTmmr/nWlX3xF4dMu/6fX2JpjjtrQvF8rRwDP8RuemEJ5m/Kd2VfvHuX3HToaBHOGB3dNOOyGPR6GPgW+PM7K/Hu8p2u7Y8XeSlSbj+yk0HvDAPfIm9n7sAD765ybX8Na1XH16P7u7Y/CqcXF2Zj3MfRX1sq1bFpHXx63yWu7c8mDHwL7fzhIC58ZI6r+7y8W3M8M/RsV/dJwTVnwx7c+uIyV/f56LXdcH1Ka1f3aRsGvsXcmv52ou6tG+DDu3u5vl/yt7eWbceD7612fb8rR/dH/VqcLuwGBj4BcPcC74m2ThqMuDiO+YfR2RNm4fsDR43sm+Pz7mPg03F+PFyAbmM/Nbb/Z4eejdRuzY3tn8zamrcffZ74wtj+P7y7F7q3bmBs/7Zj4FOFRry/Gq8viW6htqrMf6A3WjeqZfQYFDm3p09WhL352GDgU5VMjfWX5927zkdKcqOYHItO5ubdrlXhA3pij4FPjri1YJsTmyYOQvV4PlbZbUopXP3coqiXG3ZiwYO90aoh39F5hYFPEXPjoSyRyhzVD43rJHhy7CA6UliE00d587P64i+Xou0ptT05Nh0vJoEvItcBGAvgDAA9lVKZZb43HMDtAIoA/EEpNbOq/THw/SeWwz6Vsf2awL5DBeg+ztxFd11Z4wagdkI1r8ugE+gGfrQ/uTUArgHwrxMO3hnADQC6AGgBYLaIdFRKFUV5PIoxETnuwttTszfhydkbY17HRY9WfEPZr89tg9FXdA70uHFxscIri7dhzEdZXpdynLan1MIXf+ntdRnkEleGdERkLoA/l/bwS3r3UEpNLvl8JoCxSqkvK9sPe/jBUlhUjNNiPO4frV+0bYj7L+uIc9o1Mnr9oLhYYU3uPkxdkI0PVuQaO44JHFYLnlj18CvSEsDiMp/vLPkahUi1+Ljjev8HjhSiy5gqR+48tXzbDxj67yVel+Er0/5wETq34JOjbFBl4IvIbADNyvnWSKXUhxVtVs7Xyn0rISLDAAwDgDZt2lRVDvlY7YRqx/0BKC5W6DhqOgqL/TMxgLikgc2qDHylVL8I9rsTQNnVkFoBKPd9rVJqCoApwLEhnQiORT4VFyfYPGnwcV/b/t1BXPyYuwu8UcWeGXoWLu/WwusyyCdMDel8BOB1Efkrjl207QBgqaFjUYC0OaXWSXdfBmEoKAg+u/8SnJpUx+syyMeiCnwRuRrA3wEkAcgQkRVKqQFKqSwReRvAWgCFAO7mDB2qyIlDQaUOHCnE5X9fgOz8Ax5U5V+L0vqgRYOaXpdBAcQbryiQlFKYmfUt7nr1K69Lcd39l3XE73qfhniuREqavJ6lQ2SUiGBg1+Zai3Pl7z+CN5ZsxxOzYn//QKlbeyXjtl7trL55jLzHHj4RUcDp9vC5chURkSUY+ERElmDgExFZgoFPRGQJBj4RkSUY+ERElmDgExFZgoFPRGQJX914JSJ5ALYZ2n1jAPmG9h0GfH2qxteocnx9Kmfy9WmrlEqqqpGvAt8kEcnUuRPNVnx9qsbXqHJ8fSrnh9eHQzpERJZg4BMRWcKmwJ/idQE+x9enanyNKsfXp3Kevz7WjOETEdnOph4+EZHVQh34InKdiGSJSLGIpJzwveEisllENojIAK9q9BMRGSsi34jIipKPwVVvFX4iMrDkPNksImle1+NHIpIjIqtLzhvrH2ohIlNFZI+IrCnztUYiMktENpX82zDWdYU68AGsAXANgHllvyginQHcAKALgIEAnhOR+NiX50tPKqV6lHxM87oYr5WcF88CGASgM4AbS84fOlnvkvOGUzOBl3AsW8pKA/CZUqoDgM9KPo+pUAe+UmqdUmpDOd8aAuBNpdQRpVQ2gM0Aesa2OgqIngA2K6W2KqWOAngTx84fogoppeYB+P6ELw8B8HLJ/18GcFVMi0LIA78SLQHsKPP5zpKvEXCPiKwqeUsa87ecPsRzRY8C8KmILBeRYV4X41NNlVK7AKDk3yaxLiDwDzEXkdkAmpXzrZFKqQ8r2qycr1kxXamy1wvAPwBMwLHXYgKAJwDcFrvqfMnac8WhXkqpXBFpAmCWiKwv6eWSjwQ+8JVS/SLYbCeA1mU+bwUg152K/E339RKR5wF8YricILD2XHFCKZVb8u8eEXkfx4bCGPjH2y0izZVSu0SkOYA9sS7A1iGdjwDcICIJItIOQAcASz2uyXMlJ2Gpq3HsorftlgHoICLtRKQGjl3s/8jjmnxFRGqLSN3S/wPoD5475fkIwG9K/v8bABWNQBgT+B5+ZUTkagB/B5AEIENEViilBiilskTkbQBrARQCuFspVeRlrT7xqIj0wLEhixwAd3pbjveUUoUicg+AmQDiAUxVSmV5XJbfNAXwvogAxzLldaXUDG9L8paIvAHgUgCNRWQngDEA0gG8LSK3A9gO4LqY18U7bYmI7GDrkA4RkXUY+ERElmDgExFZgoFPRGQJBj4RkSUY+ERElmDgExFZgoFPRGSJ/wcviy4Ayrnc5gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#x_com =\n",
    "#y_com =\n",
    "\n",
    "# COM Plotting\n",
    "x_com = (m1*fSoln[:,0] + m2*fSoln[:,4] + m3*fSoln[:,8])/(m1+m2+m3)\n",
    "y_com = (m1*fSoln[:,1] + m2*fSoln[:,5] + m3*fSoln[:,9])/(m1+m2+m3)\n",
    "\n",
    "\n",
    "# track\n",
    "plot(fSoln[:,4] - x_com[:],fSoln[:,5] - y_com[:], 'g')\n",
    "plot(fSoln[:,8] - x_com[:], fSoln[:,9] - y_com[:],)\n",
    "plot(fSoln[:,0] - x_com[:],fSoln[:,1] - y_com[:], 'r')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Making a movie!\n",
    "\n",
    "#### Exercise 7.7\n",
    "Now let's take a true look at the motion of the three masses.  Fill in the code below to create a movie of the orbital motion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using matplotlib backend: Qt5Agg\n",
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "# Movie\n",
    "\n",
    "%pylab\n",
    "# get rid of inline plotting -- movies won't work if they're run within the notebook\n",
    "\n",
    "# create an x-y grid for plotting.  Choose appropriate x and y limits on your axes.\n",
    "fig=figure()\n",
    "ax = axes(xlim=[-15,15], ylim=[-15,15])\n",
    "grid()\n",
    "\n",
    "# Fill in expressions for the coordinates of each object (i.e. x1 = fSoln[:,0])\n",
    "\n",
    "x1 =fSoln[:,0]\n",
    "y1 =fSoln[:,1]\n",
    "\n",
    "x2 = fSoln[:,4]\n",
    "y2 = fSoln[:,5]\n",
    "\n",
    "x3 = fSoln[:,8]\n",
    "y3 = fSoln[:,9]\n",
    "\n",
    "\n",
    "# plot orbital trajectories in the x-y plane\n",
    "\n",
    "plt.plot(x3-x_com,y3-y_com, color='red')\n",
    "plt.plot(x2-x_com,y2-y_com, color='blue')\n",
    "plt.plot(x1-x_com,y1-y_com, color='green')\n",
    "\n",
    "# initialize plotting of the positions of m1 and m2\n",
    "# also plot the position of the center of mass\n",
    "m3_position,= ax.plot([],[], 'ro', ms=6)\n",
    "m2_position,= ax.plot([],[], 'b*', ms=10)\n",
    "m1_position, = ax.plot([],[], 'g*', ms=10)\n",
    "\n",
    "\n",
    "# FuncAnimation requires an initialization function.  This just\n",
    "# lets FuncAnimation know that this line reference is the one to use\n",
    "# The function does no plotting since set_data is empty.\n",
    "def init():\n",
    "    m1_position.set_data([], [])\n",
    "    m2_position.set_data([], [])\n",
    "    m3_position.set_data([], [])\n",
    "    return m1_position, m2_position, m3_position\n",
    "\n",
    "# function used in the FuncAnimation\n",
    "def animate(i):\n",
    "    # This makes the plot, by moving data into lines's set_data method\n",
    "    \n",
    "    j=i*10  # Animate every 10 points\n",
    "    \n",
    "    m3_position.set_data(x3[j]-x_com[j],y3[j]-y_com[j]) \n",
    "    m2_position.set_data(x2[j]-x_com[j],y2[j]-y_com[j])\n",
    "    m1_position.set_data(x1[j]-x_com[j],y1[j]-y_com[j]) \n",
    "    \n",
    "    return m1_position, m2_position, m3_position\n",
    "\n",
    "num_frames = 200  # Animate 200 framces\n",
    "\n",
    "# for blit True, plot only changes from previous plot to speed up plot calls\n",
    "# blit must be False on osx (changing the backend will also work, but this is easier)\n",
    "# use non-inline plots for animation on osx\n",
    "anim= animation.FuncAnimation(fig, animate, init_func=init, frames=num_frames, blit=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# put back inline plotting\n",
    "%pylab inline "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Something a little different\n",
    "\n",
    "#### Exercise 7.8\n",
    "So far, you have just looked at stable configurations. These are actually rare in a 3-body system. Try changing the initial conditions until you get an unstable orbit. (Reducing the velocity of M_3 makes some fun plots.)\n",
    "\n",
    "You are welcome to either copy and paste your code from exercise 6.5 (and possibly 6.6) below. Or you can just alter the initial conditions in exercise 6.5 and rerun your code, provided you put it back to as it was before you submit your lab. Regardless, state what masses you used and your initial conditions ***here***.\n",
    "\n",
    "\n",
    "\n",
    "***Changing m-3 velocity to sqrt(97) causes the orbits to double back on themselves. The oribt did not loop over itself.***\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check-out\n",
    "#### Exercise 7.9\n",
    "\n",
    "Briefly summarize in the cell below the ideas in today's lab."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "We looked at orbits, 2 body and 3 body systems, seeing how the relative masses impacted the x and y center of mass coordinates. We also observed stable and unstable orbits. Finally, we learned how to animate a system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Challenge Problems\n",
    "\n",
    "If you finish early, here are some additional exercises to attempt:  \n",
    "* Play around with your integrations of two planets orbiting a single star and find  \n",
    " configurations with unstable 2-planet orbits. \n",
    "* Examine conservation of energy for the 3-body systems that we modeled today\n",
    "* 3-body 3-D integrator – Extend your 3-body integrator to solve the equations of motion in 3-D  \n",
    " and examine the motion along the z-axis. \n",
    "* Lagrange points are special locations where a small test mass will orbit the Sun with  \n",
    " exactly the same orbital period as the Earth.  Look up (google) the location of these  \n",
    " points, and generate a movie of the Earth/Sun system with a third low-mass object orbiting  \n",
    " at a Lagrange point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}