Lecture.21-Three-Body-Problem.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('default')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Three-Body Problem\n",
    "## Lecture 21"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Follows Giordano 4.4: The Three-Body Problem and the Effect of Jupiter on Earth**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Up to this point, we have only considered two-body solar systems.  The real solar system has multiple planets and this makes the an analytical theory much more difficult.\n",
    "\n",
    "Numerically, we can tackle the problem by explicitly calculating the pairwise forces of gravity in a multiple planet solar system.  For concreteness, let's start with the three-body problem of an Sun-Earth-Jupiter system.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start by reproducing the simulation of the orbit of single planet around the sun.  This code solve for the orbit of a single planet (this time ignoring the effects of general relativity).  This is the same code as given in the previous lecture with the RK4 scheme set up more compactly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "π = np.pi\n",
    "\n",
    "def orbit1(tmax = 1.0, dt = 0.001, v0 = 2*π, x0 = 1):\n",
    "    \"\"\"\n",
    "    Calculate the orbit of planet\n",
    "    \"\"\"\n",
    "\n",
    "    # number of timesteps\n",
    "    N = round(tmax/dt)\n",
    "\n",
    "    # allocate arrays\n",
    "    x = np.zeros(N)\n",
    "    y = np.zeros(N)\n",
    "    vx = np.zeros(N)\n",
    "    vy = np.zeros(N)\n",
    "    t = np.zeros(N)\n",
    "    \n",
    "    # initialise position of planet in AU\n",
    "    x[0] = x0\n",
    "    y[0] = 0\n",
    "\n",
    "    # initialise velocity of planet in AU/yr\n",
    "    vx[0] = 0\n",
    "    vy[0] = v0\n",
    "\n",
    "    def F(q):\n",
    "        x, y, vx, vy = q\n",
    "        r = np.sqrt(x**2+y**2)\n",
    "        dxdt = vx\n",
    "        dydt = vy\n",
    "        dvxdt = - (4*π**2*x)/(r**3)\n",
    "        dvydt = - (4*π**2*y)/(r**3)\n",
    "        \n",
    "        dqdt = np.array([dxdt, dydt, dvxdt, dvydt])\n",
    "        return dqdt\n",
    "        \n",
    "    # loop over the timesteps\n",
    "    for i in range(N-1):\n",
    "        t[i+1] = t[i] + dt\n",
    "        \n",
    "        # 4th order Runge-Kutta\n",
    "        q0 = np.array([x[i], y[i], vx[i], vy[i]])\n",
    "        s1 = F(q0)\n",
    "        q1 = q0 + s1*dt/2\n",
    "        s2 = F(q1)\n",
    "        q2 = q0 + s2*dt/2\n",
    "        s3 = F(q2)\n",
    "        q3 = q0 + s3*dt\n",
    "        s4 = F(q3)\n",
    "        q4 = q0 + (s1 + 2*s2 + 2*s3 + s4)/6 * dt\n",
    "        \n",
    "        x[i+1], y[i+1], vx[i+1], vy[i+1] = q4\n",
    "              \n",
    "    return x, y, vx, vy, t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our plotting tool to show the orbit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_orbit(x, y, t, t0=None):\n",
    "    \"Make a plot of planet showing the Sun\"\n",
    "    \n",
    "    # Plot the Sun at the origin\n",
    "    plt.plot(0, 0, '.',color=\"yellow\", markersize=20, markeredgecolor='black')\n",
    "  \n",
    "    plt.plot(x, y, 'k-', linewidth=0.5)\n",
    "  \n",
    "    # Plot the planet at a specific position\n",
    "    if t0 is None:\n",
    "        i = -1\n",
    "    else:\n",
    "        i = int(t0/t[-1]*len(t))\n",
    "        if i >= len(t):\n",
    "            i = len(t)-1\n",
    "    plt.plot(x[i], y[i], '.', color= 'k', markersize=15)\n",
    "\n",
    "    plt.xlabel('x (AU)')\n",
    "    plt.ylabel('y (AU)')\n",
    "    plt.axis('equal')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the planetary data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = { 'Mercury' : {'name': 'Mercury', 'm':2.4e23, 'a':0.39,  'e':0.206},\n",
    "         'Venus'   : {'name': 'Venus',   'm':4.9e24, 'a':0.72,  'e':0.007},   \n",
    "         'Earth'   : {'name': 'Earth',   'm':6.0e24, 'a':1.00,  'e':0.017},  \n",
    "         'Mars'    : {'name': 'Mars',    'm':6.6e23, 'a':1.52,  'e':0.093},\n",
    "         'Jupiter' : {'name': 'Jupiter', 'm':1.9e27, 'a':5.20,  'e':0.048},\n",
    "         'Saturn'  : {'name': 'Saturn',  'm':5.7e26, 'a':9.54,  'e':0.056},\n",
    "         'Uranus'  : {'name': 'Uranus',  'm':8.8e25, 'a':19.19, 'e':0.046},\n",
    "         'Neptune' : {'name': 'Neptune', 'm':1.03e26,'a':30.06, 'e':0.010},\n",
    "         'Pluto'   : {'name': 'Pluto',   'm':1.3e22, 'a':39.53, 'e':0.248},\n",
    "        }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Earth's orbit (two-body problem)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From conservation of energy and conservation of angular momentum, we can show that the initial velocity can be calculated as\n",
    "\n",
    "$$ v_0 = \\sqrt{G M_S \\frac{1}{a} \\frac{(1 - e)}{(1+e)} } $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEWCAYAAABBvWFzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd4VGXax/HvnUDoGDQ06c0GImpAEFwVcAVdwb7AqtgWFQuI8ArLCmJZlUVARddFRAFdhLUsqCBKsS2gBKWKQgiIoQUIRYpAkvv9Y07YIUwamZknM3N/rmuuzJzznDO/M4Hcc9rziKpijDHGlFSc6wDGGGOigxUUY4wxQWEFxRhjTFBYQTHGGBMUVlCMMcYEhRUUY4wxQWEFxUQMEWkoIioiZUK0/ttF5OsC5s8Wkd6heO9wE5HHReStULU3sckKigkaEdkoIodEZL/fY1wJ19c5mBlLQlW7quokKLz4RDIRuUxE0kO4/roi8p6I7BSRvSKyUkRuD9X7mfAJyTc9E9OuUdW5JVmBiJRR1axgBSrqe4bz/QK9f7i32aEpwHKgAXAYOBeo5TSRCQrbQzFhISJNRGS+iOzyvpm+LSKJfvM3isijIrICOCAiU4H6wIfens7/+a3uTyKyyVvP0ALe8xQRmSwiO0TkZxH5q4jEefNuF5H/isgYEckEHv/fYvKS9835RxHp5Le+z0XkbhE5G3gVaOdl25PP+98hImtE5FcRSRORe/zmXSYi6d42bwPe8Kb/WURSRSRTRGaKyOne9BMO9+Xm8duer0VklIjsFpENItLVr20jEfnCy/IZkJRP5krAbOB0v73M073ZCd7n+auIrBaRZL/lTvf2OnZ47/1Qfr8XoDXwpqoeUNUsVf1eVWf7fy55Mh3bU/UOvU3PL4dxywqKCRcBngFOB84G6vG/P+K5egJXA4mq2hPYhG+Pp7KqjvRr1wE4E+gEDPP+wAfyEnAK0Bi4FLgNuMNv/kVAGlADeDrPtCRgOPC+iJzqv1JVXQPcCyzysiUSWAbwB6Cq975jROQCv/m1gFPxfVPvIyId8X1GNwO1gZ+Bd/JZdyAXAT952UcCr4uIePP+BSz15j0JBDwXpKoHgK7AFm/bKqvqFm92Ny9PIjATGAfgFekP8e111MH3e+kvIlfmk3Mx8LKI9BCR+sXYvlwBcxj3rKCYYPuPiOzxe/wZQFVTVfUzVT2sqjuA0fj+yPt7UVV/UdVDhbzHCFU9pKrL8f0ROy9vAxGJB/4IDFHVX1V1I/A8cKtfsy2q+pL3LTn3PTOAsap6VFWn4fsDfXXxPgIfVf1YVderzxfAp8Alfk1ygOHeZ3II+BMwUVW/U9XDwBB8e0ENi/iWP6vqa6qaDUzCV5Rqen+0WwOPee/1Jb4CUFxfq+osb/1T+N/n3hqorqpPqOoRVU0DXgN65LOem4CvgMeADSKyTERaByGHccwKigm2a1U10e/xGoCI1BCRd0Rks4jsA97ixMMuvxTxPbb5PT8IVA7QJglIwPctP9fP+L5BF/R+m/X4HlN/xrdXVWwi0lVEFnuHr/YAV3H8Nu9Q1d/8Xp/un1dV9wO78mQuyLHPRVUPek8re+vd7e195PL/XIoq7+de3jsE1wDfIbJjXySAvwA1A61EVXer6mBVbe61WYbvi4gEal+MHMYxKygmXJ4BFGipqlWBW/AdBvOXt+vrknSFvRM4iu+PXa76wOZC1l8nzx+2+sCWAO0KzCYi5YD3gFFATe+w2CyO3+a869jin9c7n3Galzm3GFT0a1/UE9lbgWre+nIVdKipuJ/7L8CGPF8kqqjqVYUtqKo78X1Gp+M7/HcAv2309jSrFzOPccQKigmXKsB+YI+I1AEGFWGZ7fjOfxSbdzhkOvC0iFQRkQbAAHx7RgWpATwkImVF5CZ853tm5ZOtrogk5LOeBKAcsAPI8k6Q/76Q9/4XcIeItPIK0t+Ab1R1o3eYcDNwi4jEi8idQJNC1geAqv4MpAAjRCRBRDoA1xSwyHbgNBE5pSjrB74F9nkXGFTw8rXI7zCWiDznzS8jIlWA+4BUVd0FrMW3x3G1iJQF/orvczQRwAqKCbbcq7JyHx9400cAFwB7gY+B94uwrmeAv3qHUQaeRJYH8X3jTQO+xvcHe2Ihy3wDNMO3h/M0cKP3hy6v+cBqYJuI7Mw7U1V/BR7CV9R2A73wnUDOl6rOw3de4T18exVNOP48xJ/xFeJdQHNgYSHb4q8XvpP2mfguNphcQI4fgalAmvfZF3jIzyve1wCtgA34PrsJ+C6ICKQi8AGwB9/vpgG+E+2o6l6gr7d87p5ZyO6JMcElNsCWMcaYYLA9FGOMMUFhBcUYY0xQWEExxhgTFFZQjDHGBEVM3QyUlJSkDRs2dB3DGGMiytKlS3eqaqH3A8VUQWnYsCEpKSmuYxhjTEQRkSL1rGCHvIwxxgSFFRRjjDFBYQXFGGNMUFhBMcYYExRWUIwxxgSF04IiIhNFJENEVuUzX0TkRfENibrCf7Q7EektIuu8R8DR54wxxoSP6z2UN4EuBczviq/n12ZAH+AfAN6QrMPx9Z7aBhguItVCmtQYY0yBnBYUbyjSzAKadAcme0OoLgYSRaQ2cCXwmapmqupu4DMKLkzGGGNCrLTf2FiH44dpTfem5Tf9BCLSB9/eDfXrFzRInTEld+jQIdavX09aWhpbt25l+/btZGdnn9Aud1BIVT32PFfukBL+01UVVSUxMZE6depQp04dGjVqxOmnn37C8sa4UtoLSqD/KVrA9BMnqo4HxgMkJyfb4C8mKI4cOcKqVatYunQpmzZtOvZHvUKFCjRt2pRGjRqRnJxMjRo1KFMmOP/NVJW9e/eyefNmtmzZwkcffcSWLVuOFaAKFSpw7rnn0rp1a2rWDDicuzEhVdoLSjpQz+91XXzjbqcDl+WZ/nnYUpmYs3XrVhYsWMCaNWsAKFu2LC1atODKK6+kXr16YdlLEBESExNJTEykefPmJ8w/ePAgK1asYPr06WRkZABQtWpV2rdvT+vWrSlbtmzIM5rY5nzERhFpCHykqi0CzLsaeAC4Ct8J+BdVtY13Un4pviFlAb4DLlTVgs7HkJycrNaXlymKffv2MWfOHFasWIGqUrt2bS677DLOOeeciDrEtHfvXr766iuWLFlCVlYW1atXp1OnTrRo0SKitsO4JSJLVTW50HYuC4qITMW3p5EEbMd35VZZAFV9VXz/4sfhO+F+ELhDVVO8Ze8E/uKt6mlVfaOw97OCYgqyadMmPvzwQ9LT06latSpdunThvPPOIy7O9cWQwZORkcG8efNYuXIlCQkJdOzYkfbt2xMfH+86minFIqKghJsVFJNXRkYG06dPJz09nfr163PNNddQr169wheMAocPH2b+/PksXLgQVaVLly60b9/e9lzMCaygBGAFxYDvD+mMGTNYunQp1atX5+abb475KwCPHj3K7NmzWbhwIdWqVeOPf/wjNnaQyWUFJQArKLHtl19+YfLkyRw6dIjrrruOCy64wL6NB7Br1y6mTZtGWloaHTt2pEuXLlF12M8UnxWUAKygxB5V5YsvvmDWrFnUqVOH2267jWrVrFOFolBVFixYwCeffEL9+vW59dZbOeWUU1zHMg5YQQnACkrsUFVmzZrF3Llzufzyy/nDH/5g37JLYOPGjUyaNImyZcvSp08fkpKSXEcyYWQFJQArKNEvJyeH999/n4ULF3L11VfTsWNHO6wVRBkZGYwfPx5VpU+fPnYDZYywghKAFZTo9umnnzJr1ixuuukm2rdv7zpOVNu1axfjx48nJyeHBx54wA6FRTkrKAFYQYlOy5Yt480336Rz585cffXVtkcSRtu2bWPcuHHUqVOHu+++2+7Gj1JWUAKwghJdtm/fztixY2nWrBm33XZb0PrMMsW3evVqJk6cyGWXXcY111zjOo4JMisoAVhBiQ7Z2dlMmDCBLVu2MGDAADvcUorMnj2buXPn0r9//5i5QTQWFLWg2Fc6E1G+++473nzzTe644w7OP/9813FMHl27duXSSy9l7NixVKlShfvuu8/2HGOI7aGYiHDkyBHGjBlD1apVueeee+wS4AiwevVqXn31Ve69996AvSObyGGHvAKwghKZli9fzoQJE+jXrx9NmzZ1HccUQ1ZWFi+++CIJCQn07dvXvghEqKIWFPvtmlIrJyeHcePG8eWXXzJ27FgrJhGoTJkyDBgwgPbt29OvXz82btzoOpIJISsoplTauXMnAwYMoG3btjz44IPWvXqEO//88/n73//OG2+8wcyZM13HMSFiBcWUOl999RXPPfccTzzxBMnJhe5lmwhRvnx5RowYweHDh3nqqafIyspyHckEmRUUU2qoKq+99hpr1qxh5MiRVK1a1XUkEwI33XQTN910E/3792fr1q2u45ggclpQRKSLiPwkIqkiMjjA/DEissx7rBWRPX7zsv3m2T50hDt8+DBDhw6lWbNm9OnTx+52j3Jnnnkmzz33HKNGjeK7775zHccEibOrvEQkHlgLXAGkA0uAnqr6Qz7tHwTOV9U7vdf7VbVycd7TrvIqnXbs2MGIESMYNGgQDRo0cB3HhFFOTg6jRo3irLPOolu3bq7jmHxEwlVebYBUVU1T1SPAO0D3Atr3BKaGJZkJm/Xr1/Pkk0/y7LPPWjGJQXFxcfzf//0f+/fvZ9y4ccTSbQzRyGVBqQP84vc63Zt2AhFpADQC5vtNLi8iKSKyWESuze9NRKSP1y5lx44dwchtguTbb7/ltdde4/nnn6dy5WLtbJoo06tXL1q1asXjjz9OTk6O6zjmJLksKIEOkuf39aQH8K6qZvtNq+/tgvUCxopIk0ALqup4VU1W1eTq1auXLLEJmjlz5rBgwQKeeeYZ66HWANChQwduuOEGhgwZYleARSiXBSUd8O89ri6wJZ+2PchzuEtVt3g/04DPAevYKUJ8+OGHrFu3jkcffdROvpvjtGzZkrvvvpuBAwfy22+/uY5jisllQVkCNBORRiKSgK9onHC1loicCVQDFvlNqyYi5bznSUB7IODJfFO6vPfee2zdupUHHnjAdRRTSjVr1oxHHnmEQYMGcfDgQddxTDE4KyiqmgU8AMwB1gDTVXW1iDwhIv6Xe/QE3tHjz9adDaSIyHJgAfBsfleHmdLj3//+N/v27aNPnz6uo5hSrl69evzlL3/h0Ucf5fDhw67jmCKyziFNWHz00Uds3bqVP//5z66jmAiyadMmRo8ezd///nc71+ZQJFw2bGLE/PnzWb9+vRUTU2z169fnoYce4tFHH7UT9RHACooJqW+//ZbFixfTr18/11FMhGrcuDH33HMPw4YNs/tUSjkrKCZk0tLSeP/99xkyZIjrKCbCnXnmmfzhD39g7NixrqOYAlhBMSGxe/duxowZwxNPPGGXBpuguPjii2nQoAFTp1qHGaWVFRQTdEeOHOGxxx7jySefJCEhwXUcE0Wuv/56MjMz+fzzz11HMQFYQTFBN2LECAYOHEhiYqLrKCYK9e3bl08//ZSff/7ZdRSThxUUE1Rvvvkml19+OQ0bNnQdxUQpEWHYsGGMHDnS7qYvZaygmKBZsmQJe/fupXPnzq6jmChXvnx5Bg8ezIgRI+zKr1LECooJil27djF16lQeeugh11FMjKhXrx5XXnklb7zxhusoxmMFxZSYqvL0008zbNgwu6LLhNVll11GZmYmP/xgPS+VBlZQTIlNnDiRG264wU7CGyf69+/PK6+8Yn1+lQJWUEyJrFmzhl27dtG+fXvXUUyMKlOmDAMHDmTkyJGuo8Q8KyjmpGVlZfHKK68wYMAA11FMjGvYsCFNmzZl7ty5rqPENCso5qS9/PLL3HfffZQpU8Z1FGPo0aMHs2bN4sCBA66jxCwrKOakrF27lqysLM455xzXUYwBfPenPPLIIzz//POuo8QspwVFRLqIyE8ikioigwPMv11EdojIMu9xt9+83iKyznv0Dm/y2JaTk8NLL71klwibUqdOnTrUrl2bJUuWuI4Sk5wVFBGJB14GugLnAD1FJNDX3Wmq2sp7TPCWPRUYDlwEtAGGi0i1MEWPeVOmTOGWW26xAY9MqXTXXXcxadIkjh496jpKzHG5h9IGSFXVNFU9ArwDdC/islcCn6lqpqruBj4DuoQop/Gzd+9efvzxRy666CLXUYwJKC4ujnvvvZcJEya4jhJzXBaUOsAvfq/TvWl53SAiK0TkXRGpV8xlTZCNHTuW/v37u45hTIFatGjBli1byMzMdB0lprgsKIFuqc7bKc+HQENVbQnMBSYVY1lfQ5E+IpIiIik7duw46bAGfvjhB5KSkqhZs6brKMYUql+/frz44ouuY8QUlwUlHajn97ousMW/garuUtXc219fAy4s6rJ+6xivqsmqmly9evWgBI9VEydOpE+fPq5jGFMkSUlJJCUl8eOPP7qOEjNcFpQlQDMRaSQiCUAPYKZ/AxGp7feyG7DGez4H+L2IVPNOxv/em2ZCZNGiRVx44YV2It5ElD59+vD666+7jhEznBUUVc0CHsBXCNYA01V1tYg8ISLdvGYPichqEVkOPATc7i2bCTyJrygtAZ7wppkQUFWmTZvGH//4R9dRjCmWhIQEWrRowffff+86SkyQWBpLIDk5WVNSUlzHiDhz5swhOzubq666ynUUY4otKyuLRx991G54LAERWaqqyYW1szvlTaHmzJlD165dXccw5qSUKVOG888/3252DAMrKKZAX331FR06dLBxTkxE69mzJ9OmTXMdI+pZQTEFmjFjBtdee63rGMaUSHx8PGeeeSZr1qwpvLE5aVZQTL6WL1/OueeeS1yc/TMxke+WW27h7bffdh0jqtlfCpOv6dOn06NHD9cxjAmKChUqULlyZewG59CxgmICyszMpFKlSpQrV851FGOC5o477mDixImuY0QtKygmoMmTJ3Pbbbe5jmFMUNWsWZO9e/eSlZXlOkpUsoJiTpCTk8PWrVupW7eu6yjGBF2XLl2YM8c61ggFKyjmBF988QWXX3656xjGhMQll1zCl19+6TpGVLKCYk7w2Wef0blzZ9cxjAkJEaFmzZps27bNdZSoYwXFHGf//v2UL1+eMmXKuI5iTMj06tWLqVOnuo4RdaygmON88MEHXH/99a5jGBNStWrVIiMjw3WMqGMFxRxnzZo1tGjRwnUMY0KuUaNGbNiwwXWMqGIFxRzz66+/UrlyZdcxjAmLa6+9lv/85z+uY0QVKyjmmNmzZ1sX9SZm1KhRww57BZkVFHPMsmXLOO+881zHMCZs6tatS3p6uusYUcNpQRGRLiLyk4ikisjgAPMHiMgPIrJCROaJSAO/edkissx7zMy7rCmerKwsypYta93Um5hy5ZVX8umnn7qOETWcFRQRiQdeBroC5wA9ReScPM2+B5JVtSXwLjDSb94hVW3lPbphSmTJkiW0adPGdQxjwqpp06asX7/edYyo4XIPpQ2QqqppqnoEeAfo7t9AVReo6kHv5WLA+gIJkS+++IJLL73UdQxjwi4+Pp7s7GzXMaKCy4JSB/jF73W6Ny0/dwGz/V6XF5EUEVksIvmOACUifbx2KdZtdf72799vV3iZmJScnGzDAweJy4IS6GC9BmwocguQDPzdb3J9VU0GegFjRaRJoGVVdbyqJqtqcvXq1UuaOSodOHCASpUquY5hjBO/+93v+Prrr13HiAouC0o6UM/vdV1gS95GItIZGAp0U9XDudNVdYv3Mw34HDg/lGGjmZ0/MbEsMTGRvXv3uo4RFVwWlCVAMxFpJCIJQA/guKu1ROR84J/4ikmG3/RqIlLOe54EtAd+CFvyKPPNN99YQTExTzXgARJTDM4KiqpmAQ8Ac4A1wHRVXS0iT4hI7lVbfwcqA//Oc3nw2UCKiCwHFgDPqqoVlJO0f/9+qlSp4jqGMc7Y1V7B4bRLWVWdBczKM22Y3/OAfair6kLg3NCmiw32rcwYaNeuHQsXLqRp06auo0Q0u1M+xm3ZsoU6dQq6uM6Y6Ne0aVPS0tJcx4h4VlBi3MqVK2nZsqXrGMY4FRcXR05OjusYEc8KSoxbtWqVdVdvjAkKKygxbt++fVStWtV1DGOcO/XUU9m1a5frGBHNCkqMs84gjfFp3rw5q1evdh0johWpoHj3fTQXkcYiYkXIGBN1GjduzMaNG13HiGj5XjYsIqcA9wM9gQRgB1AeqCkii4FXVHVBWFKakNi3b5/df2KMp169emzatMl1jIhW0H0o7wKTgUtUdY//DBG5ELhVRBqr6uuhDGhCZ8OGDTRq1Mh1DGNKhYSEBI4cOeI6RkTLt6Co6hUFzFsKLA1JIhM2GzdupEGDBoU3NMaYIijokNcFeSYpsFNVfwnU3kSejIwMLrzwQtcxjDFRoqBDXs8HmHaq15FjT1VdFqJMJkwyMjKwLv2N+R+76rFkCjrkdXmg6SKSDLwI/C5UoUx4HD16lHLlyrmOYUypkTt6Y3x8vOsoEanYlwCragq+HoCNMSaqVKtWjd27d7uOEbGKXVBEpCb5jKxojDGR7LTTTrO75UugoJPyL3Fi4TgVuBjoF8pQxkS79evXM27c8/zrX2+xc+d+kpIq06vXLTzwwCM0aRJwNGsTBjZ6Y8kUtIeSgu/S4NxHCjAVaKuqH4YhmzFRafbs2bRt25IKFSawcOGvHD6sLFz4KxUqTKBt25bMnj3bdcSYZfeilJCqFuuBbxz4QcVdLp91dQF+AlKBwQHmlwOmefO/ARr6zRviTf8JuLIo73fhhReq+Z9hw4a5jhBzUlNTNSmpoi5cGPif6cKFaFJSRU1NTXWaMxYtWrRIO3XqpM2aNdMePXrookWLXEcqNYAULcLf2KL25ZUkIveJyJfA50DNkhYyEYkHXga6AucAPUXknDzN7gJ2q2pTYAzwnLfsOfjGoG+Oryi94q3PmFJt3Ljn+fOfj9KuXeD57drB3Xcf5eWXx4Q3WIwbPnw4nTp1Yv78+axbt47p06fTqVMnhg8f7jpaRMm3oIhIFRG5TUQ+Ab4FmgKNVbWJqg4Mwnu3AVJVNU1VjwDvAN3ztOkOTPKevwt0Et+F4t2Bd1T1sKpuwLen0iYImWKKXXMffv/611vcddfRAtvcffdR/vWvKWFKZBYvXsyoUaM4ePDgsSGxc3JyOHjwIKNGjWLx4sWOE0aOgvZQMvDtITwNNFHVR4BgHlysA/jfdZ/uTQvYRlWzgL3AaUVcFgAR6SMiKSKSsmPHjiBFN+bk7Ny5n8J6u6lf39fOhMcLL7zAoUOHAs777bffeOGFF8KcKHIVVFD+gq934X8AQ0Qk2JeeBPp6nPeqsvzaFGVZ30TV8aqarKrJdle4cS0pqTI//1xwm02bfO1MeKxdu/bYnkleOTk5rFu3LsyJIle+BUVVx6jqRUA3fH/A/wOcLiKPisgZQXjvdHwn+HPVBbbk10ZEygCnAJlFXNYUIr//RCZ0evW6hddfL1tgmwkTytKr161hSmTOOOMM4uIC/ymMi4vjjDOC8ecuNhR6Ut47x/G0qp4LtMb3Rz0Y1zUuAZqJSCOvf7AewMw8bWYCvb3nNwLzvSsOZgI9RKSciDQCmuE7z2NMqfbAA4/w2mtlWbQo8PxFi3wF5f77Hw5vsBjWr18/ypcvH3Be+fLleeihh8KcKHIVdFL+hMNKqrpSVf+iqk3ya1NU3jmRB4A5wBpguqquFpEnRKSb1+x14DQRSQUGAIO9ZVcD04EfgE+A+1U1+2SzGBMuTZo0YfLkd+nWrSJDhpRl/Xo4ehTWr4chQ8rSrVtFJk9+125uDKO2bdsycOBAKlaseOxClbi4OCpWrMjAgQNp27at44QRJL/rifFdHvwgUD/P9ASgI76rr24vyrXJpeVh96Ecb/jw4a4jxKzU1FR9+OH7tWbNqhofH6c1a1bVhx++3+4/cWjRokXarl07bd68ud2HkgdFvA+loO7ruwB3AlO9w0p78J2kjwc+BcaodWFvzElp0qQJo0ePY/Toca6jGE/btm3505/+xI033kjNmiW+1S4mFdR9/W/AK/huGiwLJAGHNM9wwCZyqZ2UN+Y4u3fvplq1aq5jRKyC9lCOUdWjwNYQZzFhFh8fT1ZWFmXKFOmfgTFR7+jRoyQkJLiOEbGK3X29iR41atTAbvY05n9sr71krKDEsFq1arFt2zbXMYwxUaLQgiIiD4iIHVSMQrVr12bz5s2uYxhjokRR9lBqAUtEZLqIdCnJvSemdGncuDFpaWmuYxhTKthY8iVXlDvl/4rvTvTXgduBdSLytxD07WXCLCkpiZ07d7qOYUypsG3bNmrXru06RkQr0jkU78aWbd4jC6gGvCsiI0OYzYSYiFgX9sZ41q1bZz0UlFBRzqE8JCJLgZHAf4FzVfU+4ELghhDnMyFmV7UY47N69WqaN2/uOkZEK8oNCEnA9ap6XKfbqpojIn8ITSwTLuXKleO3337Lt3M8Y2LF9u3b7Q75EirKOZRheYuJ37w1wY9kwql58+asXr3adQxjSgU7BFwydh9KjDvvvPNYvny56xjGOGWHfoPDCkqMa9iwIRs3bnQdwxinNm3aRIPCxmY2hbKCEuNEhJycHNcxjHFq0aJFNu5JEFhBMcdOzBsTq9asWcPZZ5/tOkbEc1JQRORUEflMRNZ5P0/o2kVEWonIIhFZLSIrROSPfvPeFJENIrLMe7QK7xZEl+TkZFJSUlzHMMaZnJycfMeVN0Xn6hMcDMxT1WbAPO91XgeB21S1Ob7BvsaKSKLf/EGq2sp72EBfJdCuXTsW5TfIuTFR7tChQ3bZfJC4Kijd8Q0hjPfz2rwNVHWtqq7znm8BMoDqYUsYQxITE9m7d6/rGMY48fXXX9O+fXvXMaKCq4JSU1W3Ang/axTUWETa4BvLfr3f5Ke9Q2FjRKRcAcv2EZEUEUmxsT/yV7ZsWQ4fPuw6hjFht2jRItq1a+c6RlQIWUERkbkisirAo3sx11MbmALcoaq5lyMNAc4CWgOnAo/mt7yqjlfVZFVNrl7ddnDy0759e/773/+6jmFM2B0+fJhy5fL9TmqKIWQFRVU7q2qLAI8ZwHavUOQWjIxA6xCRqsDHwF9kkpyHAAAVjElEQVRVdbHfureqz2HgDaBNqLYjVnTo0IGvvvrKdQxjwmrbtm3W3UoQuTrkNRPo7T3vDczI20BEEoAPgMmq+u8883KLkeA7/7IqpGljQPny5e3SYRNz5syZw+9//3vXMaKGq4LyLHCFiKwDrvBeIyLJIjLBa3Mz8Dvg9gCXB78tIiuBlfg6r3wqvPGjU+PGjVm/fn3hDY2JEj/++CNnnXWW6xhRoyi9DQedqu4COgWYngLc7T1/C3grn+U7hjRgjLr22muZMmUKAwYMcB3FmJDbt28fVapUcR0jqtidPOaY6tWr2wiOJmZ8/PHHXH311a5jRBUrKOY49erV4+efA45WYExUWbFiBS1btnQdI6pYQTHHufHGG/n3v/9deENjItju3bupUqWKjX8SZFZQzHGqV6/Ojh07bHwIE9WmTZtGjx49XMeIOlZQzAnat2/PwoULXccwJmQ2bNhA48aNXceIOlZQzAm6du3KrFmzXMcwJiRWrlxJixYtXMeISlZQzAnKli1LhQoV2LNnj+soxgTde++9xw033OA6RlSygmICuu2225g0aVLhDY2JIPv27UNEqFixousoUckKigmofv36bN68mezsbNdRjAmayZMn07t378IbmpNiBcXk65prruHjjz92HcOYoMjOzmbTpk00bNjQdZSoZQXF5KtDhw588cUXrmMYExQzZ86kW7durmNENSsoJl8iwiWXXMKXX37pOooxJaKqfPnllzYyY4hZQTEF6tatGzNmnDC6gDER5ZNPPqFr1652Z3yIWUExBYqLi6N169Z88803rqMYc1JUlU8//ZQrrrjCdZSoZwXFFOqmm25i2rRprmMYc1Js7yR8nBQUETlVRD4TkXXez2r5tMv2G1xrpt/0RiLyjbf8NG90RxMi8fHxXHbZZcydO9d1FGOKJScnhzlz5tjeSZi42kMZDMxT1WbAPO91IIdUtZX38L884zlgjLf8buCu0MY111xzDR999BE5OTmuoxhTZFOnTqVHjx62dxImrgpKdyD3NuxJ+MaFLxJvHPmOwLsns7w5OSLCzTffzPTp011HMaZIDh48yPLly2nbtq3rKDHDVUGpqapbAbyfNfJpV15EUkRksYjkFo3TgD2qmuW9TgfqhDauAbj44otZunQpBw4ccB3FmEK9/PLL9O3b13WMmBKygiIic0VkVYBH92Kspr6qJgO9gLEi0gQItO+a7+AdItLHK0opO3bsKOZWmLz69evHCy+84DqGMQXavHkzR44csbviwyxkBUVVO6tqiwCPGcB2EakN4P3MyGcdW7yfacDnwPnATiBRRMp4zeoCWwrIMV5Vk1U1uXr16kHbvlhVt25dKlasyNq1a11HMSZfo0ePpl+/fq5jxBxXh7xmArk9tPUGTrhzTkSqiUg573kS0B74QX1DCS4AbixoeRM6ffv25ZVXXrFRHU2pNGPGDDp16kTlypVdR4k5rgrKs8AVIrIOuMJ7jYgki8gEr83ZQIqILMdXQJ5V1R+8eY8CA0QkFd85ldfDmj7GJSQkcOONNzJlyhTXUYw5zr59+/jqq6+46qqrXEeJSWUKbxJ8qroL6BRgegpwt/d8IXBuPsunAW1CmdEUrEOHDsydO5fNmzdTp45dE2FKh+eff56BAwe6jhGz7E55c9IGDRrEqFGj7NCXKRXmzZvHGWecQa1atVxHiVlWUMxJq1SpEtdee60d+jLOZWZmMnv2bHr16uU6SkyzgmJK5NJLL2Xz5s121ZdxRlV5+umnGTp0qN0R75gVFFNiAwcO5IUXXuDw4cOuo5gY9Pbbb3P11VdTrVrALgFNGFlBMSVWtmxZBg4cyLPPPus6iokxq1atIj09nY4dO7qOYrCCYoKkUaNGtGzZknfffbfwxsYEwd69exk/fjyDBg1yHcV4rKCYoLnuuutITU1l5cqVrqOYKKeqPP744wwfPpz4+HjXcYzHCooJqkGDBjFhwgQyMzNdRzFR7IUXXqBXr16cdtpprqMYP1ZQTFDFx8czYsQIhg8fztGjR13HMVHovffeo1atWrRu3dp1FJOHFRQTdImJifTv35/HHnvMbno0QbV48WJ+/vlnevTo4TqKCcAKigmJJk2acN111/H888+7jmKiRFpaGh988AEPP/yw6ygmH1ZQTMhcdNFFNGvWjMmTJ7uOYiJcRkYGY8aM4cknn7SbF0sxKygmpLp3746q8t5777mOYiLU7t27GTFiBM888wwJCQmu45gCWEExIde7d2927tzJrFmzXEcxEebXX3/lr3/9K08//bSNbxIBrKCYsLjnnntYt24d8+bNcx3FRIiDBw8yZMgQRowYQWJious4pgisoJiw6devHytWrOCTTz5xHcWUcr/++iuDBg1i6NChJCUluY5jishJQRGRU0XkMxFZ5/08oVc3EblcRJb5PX4TkWu9eW+KyAa/ea3CvxXmZDz88MOkpaXxn//8x3UUU0plZmYyePBgRowYQe3atV3HMcXgag9lMDBPVZsB87zXx1HVBaraSlVbAR2Bg8Cnfk0G5c5X1WVhSW2Com/fvmRmZjJ16lTXUUwps337dh577DH+9re/2Z5JBHJVULoDk7znk4BrC2l/IzBbVQ+GNJUJmzvvvJP4+HhefPFFu/nRAJCamspTTz3FyJEjOeWUU1zHMSfBVUGpqapbAbyfNQpp3wPI+3X2aRFZISJjRKRcfguKSB8RSRGRlB07dpQstQmqm2++mVatWjFs2DCys7NdxzEOff3110ycOJHRo0dTqVIl13HMSZJQfTsUkblAoMGdhwKTVDXRr+1uVQ04Oo6I1AZWAKer6lG/aduABGA8sF5VnygsU3JysqakpBR7W0xorVmzhldffdUuDY1R77zzDhkZGTz44IN202IpJSJLVTW5sHYh20NR1c6q2iLAYwaw3SsKucUho4BV3Qx8kFtMvHVvVZ/DwBtAm1Bthwm9s88+m6FDhzJkyBBSU1NdxzFhkpOTw8iRIylbtiwPPfSQFZMo4OqQ10ygt/e8NzCjgLY9yXO4y68YCb7zL6tCkNGEUY0aNRg9ejRTpkyxy4pjQGZmJv3796dLly7ccMMNruOYIHFVUJ4FrhCRdcAV3mtEJFlEJuQ2EpGGQD3gizzLvy0iK4GVQBLwVBgymxArW7YsI0aMYOfOnYwaNcrOq0SppUuX8uSTT/LEE0/QsmVL13FMEIXsHEppZOdQIseqVav45z//yeDBg6lTp47rOCYIVJUJEyawf/9++vXrR1yc3VcdKZyfQzGmJFq0aMFzzz3HP/7xDz766CPXcUwJbdu2jYcffphzzz2Xhx9+2IpJlLLfqim1KlasyFNPPcWRI0cYOnQoe/fudR3JnIQZM2bw0ksv8dRTT9G2bVvXcUwIWUExpd71119P//79efLJJ+2EfQTZuXMngwf7OsGwS8Jjg51DMRHl/fff59tvv+WRRx6hevXqruOYAFSVt956i59++olHHnmEatUC3mJmIkhRz6FYQTERZ8+ePYwePZq6dety1113ER8f7zqS8axdu5Zx48Zx880306FDB9dxTJBYQQnACkp0WbZsGRMnTqRHjx5cfPHFruPEtMzMTF566SVOOeUU+vbtayMrRhkrKAFYQYk+OTk5TJs2jaVLl3L33Xdz1llnuY4UU44cOcJrr73G9u3befDBB+0wZJSyghKAFZTodeTIESZMmEB6ejr333+/3bsSYllZWUydOpXly5dz1113cfbZZ7uOZELICkoAVlCi3759+3j11VfZv38/d955Jw0bNnQdKapkZWXx9ttvs3LlSnr27MmFF17oOpIJAysoAVhBiR2//vorEydOZNu2bdx6662cc845riNFtIMHD/LWW2+xdu1aevXqxQUXXOA6kgkjKygBWEGJPb/99htTpkxh7dq1dOzYkSuvvNLu0i6GX375hTfffJMjR47wpz/9yc5RxSgrKAFYQYldqsqCBQuYM2cOtWrV4tZbb7UhZvORnZ3NvHnzmDdvHrVr16Z37952L0mMs4ISgBUUA7B582beeustdu/ezUUXXcRVV11FuXL5DvoZM9atW8f06dPZv38/nTp1omPHjrY3ZwArKAFZQTH+VJXFixcze/ZssrOz6dy5M5dccgllypRxHS1s0tLSmDFjBhkZGTRu3Jibb77ZxnM3J7CCEoAVFJOfo0ePMn/+fP773/+SlZVFy5Yt6dKlC4mJiYUvHEFycnL4/vvvmTt3Lnv27KFx48Z0796dGjVquI5mSjErKAFYQTFFoaqsXLmSTz75hH379lG+fHnatm1Lu3btqFSpkut4xaKqbNiwgXnz5rFx40bi4uJo1aoVnTp1irpiaUKnVBcUEbkJeBw4G2ijqgH/yotIF+AFIB6YoKq5Izs2At4BTgW+A25V1SOFva8VFHMyDh06xOLFi1m0aBEHDhygTJkynHHGGZx//vmceeaZpaovsX379pGSksKSJUvYv38/AA0aNKBz5852T445aaW9oJwN5AD/BAYGKigiEg+sxTdEcDqwBOipqj+IyHTgfVV9R0ReBZar6j8Ke18rKCYYsrOzWbduHd9//z0//fQT2dnZxMXFUbFiRZo0aUKzZs1o0qRJyLprP3LkCFu2bGHTpk389NNPpKeno6qICJUrVyY5OZnk5GSqVKkSkvc3saeoBcXJ2UdVXQMgIgU1awOkqmqa1/YdoLuIrAE6Ar28dpPw7e0UWlCMCYb4+HjOOuusE+7JOHDgAOvXryc1NZU5c+Zw8ODBY/P8v7iVK1eOihUrHnuICNnZ2WRnZ5OTk0N2djb79+9n3759xy2X+zwhIYHTTz+devXq0aVLF+rWrVvY/yVjwqI0X85SB/jF73U6cBFwGrBHVbP8pufbcZOI9AH6ANSvXz80SY0BKlWqRMuWLWnZsmW+bVSVo0ePcvDgQQ4cOMDBgwdRVeLj44894uLiqFy5MlWqVLHLdk1ECVlBEZG5QK0As4aq6oyirCLANC1gekCqOh4YD75DXkV4X2NCRkRISEggISHBToqbqBOygqKqnUu4inSgnt/rusAWYCeQKCJlvL2U3OnGGGMcKs3700uAZiLSSEQSgB7ATPUdSF4A3Oi16w0UZY/HGGNMCDkpKCJynYikA+2Aj0Vkjjf9dBGZBeDtfTwAzAHWANNVdbW3ikeBASKSiu+cyuvh3gZjjDHHsxsbjTHGFKiolw2X5kNexhhjIogVFGOMMUFhBcUYY0xQWEExxhgTFDF1Ul5EdgA/l3A1SfjuhYl00bAd0bANEB3bYdtQeoRiOxqoavXCGsVUQQkGEUkpytUOpV00bEc0bANEx3bYNpQeLrfDDnkZY4wJCisoxhhjgsIKSvGNdx0gSKJhO6JhGyA6tsO2ofRwth12DsUYY0xQ2B6KMcaYoLCCYowxJiisoBRCRG4SkdUikiMi+V6KJyIbRWSliCwTkVLXA2UxtqOLiPwkIqkiMjicGQsjIqeKyGciss77WS2fdtne72GZiMwMd85ACvtcRaSciEzz5n8jIg3Dn7JwRdiO20Vkh9/nf7eLnAURkYkikiEiq/KZLyLyoreNK0TkgnBnLEwRtuEyEdnr93sYFpZgqmqPAh7A2cCZwOdAcgHtNgJJrvOWZDuAeGA90BhIAJYD57jO7pdvJDDYez4YeC6fdvtdZy3u5wr0BV71nvcAprnOfZLbcTswznXWQrbjd8AFwKp85l8FzMY3Omxb4BvXmU9iGy4DPgp3LttDKYSqrlHVn1znKKkibkcbIFVV01T1CPAO0D306YqsOzDJez4JuNZhluIoyufqv23vAp1EJNBw1y6V9n8fRaKqXwKZBTTpDkxWn8X4RoitHZ50RVOEbXDCCkrwKPCpiCwVkT6uw5ykOsAvfq/TvWmlRU1V3Qrg/ayRT7vyIpIiIotFpDQUnaJ8rsfaqG9wub34Bo8rTYr67+MG71DRuyJSL8D80q60/z8oqnYislxEZotI83C8YcjGlI8kIjIXqBVg1lBVLerwwu1VdYuI1AA+E5EfvW8RYROE7Qj0jTis15UXtA3FWE1973fRGJgvIitVdX1wEp6Uonyuzj/7IihKxg+Bqap6WETuxbfX1THkyYIrEn4XhfkOX/9b+0XkKuA/QLNQv6kVFEBVOwdhHVu8nxki8gG+wwNhLShB2I50wP8bZV1gSwnXWSwFbYOIbBeR2qq61TsEkZHPOnJ/F2ki8jlwPr5j/64U5XPNbZMuImWAUyh9hzQK3Q5V3eX38jXguTDkCjbn/w9KSlX3+T2fJSKviEiSqoa080s75BUEIlJJRKrkPgd+DwS8+qKUWwI0E5FGIpKA7+RwqbhKyjMT6O097w2csNclItVEpJz3PAloD/wQtoSBFeVz9d+2G4H56p1dLUUK3Y485xq6AWvCmC9YZgK3eVd7tQX25h5qjRQiUiv3HJyItMH3t35XwUsFgeurFUr7A7gO3zeWw8B2YI43/XRglve8Mb4rXpYDq/EdYnKevbjb4b2+CliL7xt9qdoOfOcU5gHrvJ+netOTgQne84uBld7vYiVwl+vc+X2uwBNAN+95eeDfQCrwLdDYdeaT3I5nvP8Dy4EFwFmuMwfYhqnAVuCo93/iLuBe4F5vvgAve9u4kgKu7izF2/CA3+9hMXBxOHJZ1yvGGGOCwg55GWOMCQorKMYYY4LCCooxxpigsIJijDEmKKygGGOMCQorKMaEmYhUEJEvRCTeb9rDIvKbiJziN+12ERmXZ9nPc3uLFpG5+fW4bIwLVlCMCb87gfdVNdtvWk98Nw5eV4z1TMHXS7ExpYIVFGOCRERae50ilvd6T1gtIi0CNP0Tfnf5i0gToDLwV3yFpahmFrO9MSFlfXkZEySqusQb0OspoALwlqoe1wWP12VJY1Xd6De5J747n78CzhSRGqoasJ+yPO+32xuY6zQ9vg8tY5ywPRRjgusJ4Ap83cGMDDA/CdiTZ1oP4B1VzQHeB27ypufXjYX/9Ax83ecY45ztoRgTXKfiO3xVFl//XAfyzD/kTQdARFri61b8M68vvwQgDV9fUruAvCfdTwX8e4wt763TGOdsD8WY4BoPPAa8TYCu21V1NxAvIrlFpSfwuKo29B6nA3VEpAG+k/TtRaQWgHd1Vzm8wZ+83mRr4Rt+2hjnbA/FmCARkduALFX9l3dJ8EIR6aiq8/M0/RToAMzFd7ira575HwA9VPU5EekHzBKROGA/0NM7NAZwIbBYfSM8GuOc9TZsTJiJyPnAAFW9tYTreQGYqarzgpPMmJKxQ17GhJmqfg8s8L+x8SStsmJiShPbQzHGGBMUtodijDEmKKygGGOMCQorKMYYY4LCCooxxpigsIJijDEmKP4fyJNm/ddC07QAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(figsize=(6, 4))\n",
    "p = data['Earth']\n",
    "\n",
    "x0 = (1+p['e'])*p['a']\n",
    "\n",
    "## See equation (4.16) of Giordano, from conservation of energy and angular momentum\n",
    "v0 = 2*π*np.sqrt(1/p['a']*(1-p['e'])/(1+p['e']))\n",
    "\n",
    "x, y, vx, vy, t = orbit1(v0=v0, tmax = 1, x0 = x0)\n",
    "plot_orbit(x,y,t)\n",
    "plt.title('Earth orbit around the Sun')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extend to Three-Bodies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our first goal will be to rewrite the above routines to that they support multiple planets.  In terms of forces we need to add in the additional force between, say, Earth and Jupiter\n",
    "\n",
    "$$ F_{EJ} = \\frac{G M_J M_E}{r_{EJ}^2} $$\n",
    "\n",
    "Assuming the orbits are all coplanar, the components in the $xy$-plane are \n",
    "\n",
    "\\begin{align}\n",
    "F_{EJ,x} &= - \\frac{G M_J M_E}{r_{EJ}^2} \\cos \\theta_{EJ} = - \\frac{G M_J M_E( x_e - x_j)}{r_{EJ}^3} \\\\\n",
    "F_{EJ,y} &= - \\frac{G M_J M_E}{r_{EJ}^2} \\sin \\theta_{EJ} = - \\frac{G M_J M_E( y_e - y_j)}{r_{EJ}^3} \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The total force on Earth in the $x$-direction is the sum of the forces of gravity from the Sun and Jupiter. The equation of motion is of the form:\n",
    "\n",
    "$$ \\frac{d v_{x,e} }{dt} = - \\frac{G M_S x_e}{r^3} - \\frac{G M_J ( x_e - x_j)}{r_{EJ}^3}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that we are using AU and yr as are units for distance and time, respectively.  In that system,  we have\n",
    "\n",
    "$$ G M_S = 4 \\pi^2\\; \\mathrm{AU}^3/\\mathrm{yr}^2 $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we can calculate $G M_J$ as follows\n",
    "\n",
    "$$ G M_J = G M_S \\frac{M_J}{M_S} = 4\\pi^2 \\frac{M_J}{M_S} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Many more variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As soon as we start to deal with more than one planet, we will need many more variables to keep track of the entire system.  One approach is just to extend our code and use the arrays \n",
    "\n",
    "    x_E, y_E, vx_E, vy_E, x_J, y_J, vx_J, vy_J\n",
    "    \n",
    "where we track the position and velocity for both the Earth and Jupiter.  But this quickly will become cumbersome. With more than two planets, we would need even more variables.  \n",
    "\n",
    "Let's create a new data structure that holds all of the information about each planet.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Planet():\n",
    "    def __init__(self, data, N):\n",
    "        \n",
    "        # initialize with given planetary data\n",
    "        self.name = data['name']\n",
    "        self.m = data['m'] ## mass\n",
    "        self.e = data['e'] ## eccentricity\n",
    "        self.a = data['a'] ## semimajor axis length\n",
    "                \n",
    "        # arrays for storage of \n",
    "        self.x = np.zeros(N)\n",
    "        self.y = np.zeros(N)\n",
    "        self.vx = np.zeros(N)\n",
    "        self.vy = np.zeros(N)\n",
    "\n",
    "        # calculate initial position and velocity\n",
    "        self.x[0] = (1 + self.e)*self.a\n",
    "        if self.a > 0:\n",
    "            self.vy[0] = 2*π * np.sqrt(1/self.a * (1-self.e)/(1+self.e))\n",
    "        else:\n",
    "            self.vy[0] = 0\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is how our `orbit()` function changes when written with this new data structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Mass of the Sun in kg\n",
    "M_s = 2e30\n",
    "\n",
    "def orbit2(tmax = 1.0, dt = 0.001, planets=[]):\n",
    "    \"\"\"\n",
    "    Calculate the orbit of a planet\n",
    "    \"\"\"\n",
    "\n",
    "    # number of timesteps\n",
    "    N = round(tmax/dt)\n",
    "    t = np.zeros(N)\n",
    "\n",
    "    # allocate planets\n",
    "    p1 = Planet(planets[0], N)\n",
    "    p2 = Planet(planets[1], N)\n",
    "\n",
    "    def F(q):\n",
    "        x1, y1, vx1, vy1, x2, y2, vx2, vy2 = q\n",
    "        r1 = np.sqrt(x1**2+y1**2)\n",
    "        r2 = np.sqrt(x2**2+y2**2)\n",
    "        r12 = np.sqrt((x1-x2)**2+(y1-y2)**2)\n",
    "        \n",
    "        dxdt1 = vx1\n",
    "        dydt1 = vy1\n",
    "        dvxdt1 = - 4*π**2*x1/(r1**3) - 4*π**2*(p2.m/M_s)*(x1-x2)/(r12**3)\n",
    "        dvydt1 = - 4*π**2*y1/(r1**3) - 4*π**2*(p2.m/M_s)*(y1-y2)/(r12**3)\n",
    "        \n",
    "        dxdt2 = vx2\n",
    "        dydt2 = vy2\n",
    "        dvxdt2 = - 4*π**2*x2/(r2**3) - 4*π**2*(p1.m/M_s)*(x2-x1)/(r12**3)\n",
    "        dvydt2 = - 4*π**2*y2/(r2**3) - 4*π**2*(p1.m/M_s)*(y2-y1)/(r12**3)\n",
    "        \n",
    "        dqdt = np.array([dxdt1, dydt1, dvxdt1, dvydt1,\n",
    "                         dxdt2, dydt2, dvxdt2, dvydt2])\n",
    "                                          \n",
    "        return dqdt\n",
    "        \n",
    "    # loop over the timesteps\n",
    "    for i in range(N-1):\n",
    "        t[i+1] = t[i] + dt\n",
    "        \n",
    "        q0 = np.array([p1.x[i], p1.y[i], p1.vx[i], p1.vy[i],\n",
    "                       p2.x[i], p2.y[i], p2.vx[i], p2.vy[i]] )\n",
    "        \n",
    "        s1 = F(q0)\n",
    "        q1 = q0 + s1*dt/2\n",
    "        s2 = F(q1)\n",
    "        q2 = q0 + s2*dt/2\n",
    "        s3 = F(q2)\n",
    "        q3 = q0 + s3*dt\n",
    "        s4 = F(q3)\n",
    "        q4 = q0 + (s1 + 2*s2 + 2*s3 + s4)/6 * dt\n",
    "        \n",
    "        (p1.x[i+1], p1.y[i+1], p1.vx[i+1], p1.vy[i+1], \n",
    "         p2.x[i+1], p2.y[i+1], p2.vx[i+1], p2.vy[i+1]) = q4\n",
    "        \n",
    "    return p1, p2, t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the Runge-Kutta scheme remains basically unchanged.  But our 'state space' now has 8 variables."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example with Earth-Jupiter-Sun system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEWCAYAAABmE+CbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VNX5x/HPk0gMYZE1QXZBUXaEoChSK6CVCi4VrFhaERVXNg277NhCBSMiVlGpa7UoQsWiPxFU/FlBCCAQ5QdW9jXsW9gyz++PucEAWSEzZ5bn/XrlBZm5c893Msl97j3n3nNFVTHGGGNiXAcwxhgTGqwgGGOMAawgGGOM8VhBMMYYA1hBMMYY47GCYIwxBrCCYHIQkdoioiJyQQhkOSQidVznyCYi3UXkf/N5/hMRuTeYmQJFREaKyNuBWt6ELisIYURE1otIprexzP564TzX174Y830pIg8Ux7pUtbSq/uyt93URGVsc6w0UVe2gqm9AwcUjnInIr0VkcwDXX11EZojILhHZLyIrRaR7oNozp3O+J2iKrJOqfn4+KxCRC1T1ZHEFCnXn+35dHzFF2ef1FvA9UAs4BjQGqjhNFEXsCCFCiEhdEZkvIru9vat3RKRcjufXi8hAEVkBHBaRd4GawGzvSGNAjtX9QUQ2eusZeo55ztpL9rqjLvX+/7qIvCQic0XkoIh8JSK1zlxWRHoCfwAGeDlne89X9fYkM0RknYj0zvHakSLygYi8LSIHgO655LtIRN70Xr9BRJ4SkZgc2b8RkVQR2QOM/OVlMtnbc10tIu1yrO9LEXlAROoDLwHXeHn35fHzuU9EfvTe+88i8lCO534tIpu9z2s78Hfv8QdF5CcR2SMiH4lIVe/xs7r6ch6tZX8WIjJBRPZ6P68OOZa9xPv5HxSRuUClPDKXAj4BquY4Qq3qPR3n/TwPiki6iCTneF2en1UuWgKvq+phVT2pqstU9ZOcP5czMp06yvU+9+l55TAFs4IQOQT4C1AVqA/U4JcNWbauwC1AOVXtCmzEf8RRWlX/mmO564DLgXbAcG8jFwh/AMbg3wAtB945cwFVneo9/lcvZydvwz0b/55kNS9nXxH5TY6X3gZ8AJTLbb3AZOAioA5wPfAn4L4cz18N/AwkAk+f8VglYATwoYhUOCPvj8DDwLde3nLkbifQESjrtZsqIs1zPF8FqIB/T7mniLTF//neBVwMbADey2Pdubka+D8v+1+B10REvOf+AaR5z40Bch0LUdXDQAdgq/feSqvqVu/pW7085YCPgBcACvlZ5bQQmCIid4tIzSK8v2y55jCFYwUh/MwSkX05vh4EUNWfVHWuqh5T1QzgWfwbupyeV9VNqppZQBujVDVTVb/H/4fctPjfBgD/VtUFqnoMGIp/r7pGIV7XEqisqqNV9bg31vAKcHeOZb5V1Vmq6jvz/YpILPB7YLCqHlTV9cBE4I85FtuqqpO9vdTs1+8EnlPVE6r6T/wb2FvO4X2jqv9W1f+q31fAZ0CbHIv4gBHe55mJv3hOU9Wl3s9rMP6fV+1CNrlBVV9R1SzgDfxFJcnb6LYEhnltLcC/AS+q/1XVOd763+KX35nCfFY5dQG+BoYB60RkuYi0LIYcphBsDCH83J7bGIKIJALP49+olMFf7PeesdimQraxPcf/jwClzyFnYZzKo6qHvO6ZqhScsxb+bouc3TGx+DckZ607F5WAOPx72dk24N+Dze/1W/T02SA3eHmLzOuyGQHUw/9ZJQArcyySoapHc3xfFVia/Y3389rtZd5SiCZPfaaqesQ7OCiN/2ex19v7z7YB/xFmUZz5OxPvdWEV5rM6RVX3AoOAQSJSCZiAfyeo+vnkiKIxmPNiRwiR4y+AAk1UtSzQDX83Uk5nTm0byKluD+PfyAEgIrkNDNbI8Xxp/F0kW3NZ7sycm4B1qloux1cZVf1tPq/JaRdwAv/GKltNTt+w5vb6ajm6WbJfU5i8pxGRC4EZ+Dd2SV630hxO/7zOXMfWnHm9/vyKXubsjXlCjuULOxC7DSjvrS9bfl01Rf2dKcxnlXtDqrvw/4yq4v/dOPN3KhaoXMQ8Jh9WECJHGeAQsE9EqgH9C/GaHfj70AvFG9Qr7Abhe6ChiDQTkXjOHs8A+K2IXCcicfj7rhepam575mfm/A444A26lhSRWBFpVNiuBa87YTrwtIiUEf9g9hNAQefSJwK9RaSEiHTBP1YzJ4+81b33lZs44EIgAzjpHS3cVEDb/wDu836eFwJ/xv/zWu91EW4Bunk/ix5A3QLWB4CqbgCWAKNEJE5ErgM65fOSHUBFEbmoMOuniJ+ViIz3nr9ARMoAjwA/qepuYA3+Pf5bRKQE8BT+n6MpJlYQwk/2WUHZXzO9x0cBzYH9wL+BDwuxrr8AT3ljESmFWL4G8G0ByyiAqq4BRgOfA2uB3M7L/wf+bpM9QAv8/eS5eQ1o4OWc5W3QOwHNgHX49/hfxT9IXFi98O9x/uxl+wcwrYDXLAIu89p7GujsbajONB9IB7aLyK4zn1TVg0Bv/EVpL3AP/gHQPKnqPPz96jPw79XX5fR++Afx7wTsBhoC/yngveR0D/5B5z34P48388mxGngX+Nn7PPLtMjuHzyoBmAnsw//Z1MI/UIyq7gce9V6ffWQUsGsiopHYDXJMYYnIq8D7qvo/eTy/FBitqrMKsa7Xgc2q+lTxpjTGnCsbVDaFpqp5XoUsIg3xd6EsC14iY0xxsi4jc95EZDz+0yYHen3SxpgwZF1GxhhjADtCMMYY4wmrMYRKlSpp7dq1XccwxpiwkpaWtktVC7xmI6wKQu3atVmyZInrGMYYE1ZEpFBje9ZlZIwxBrCCYIwxxmMFwRhjDGAFwRhjjMcKgjHGGMAKgjHGGI8VBGOMMYAVBGNyparYtC4m2oTVhWnGFMbx48fZsWMH27dvP/W1Y8cOTpw4cdayInLahv/0G6Jx1uN5FYnsx0uUKEFiYiJJSUmnvqpUqULJkiXP920ZE3BWEEzYUVW2bNnC2rVrWbduHZs2bcLn8516/sILLzy1Ib744ou58sorSUxMJC4urxuYFZ/jx4+zc+dOduzYwc6dO1m9ejXbtm3j6NFfbo8sIlx88cXUqVOH+vXrU61atTwLkTHBZAXBhLSdO3fy/fffk56ezp49ewD/BrVatWrUq1ePdu3aUb16dWJjYx0n9YuLi6N69epUr573PeF9Ph/btm3jv//9L3PmzGHLll9u5Vy2bFkaNmxI8+bNSUxMDEZkY05xOv21iKwHDgJZwElVTc5v+eTkZLW5jCLX/v37SUtLIy0tjQMHDgCQmJhI06ZNadiwIRUrVnScMPD2799Peno6aWlpZGRkAJCQkMCVV15JcnJyVPwMTPETkbSCtq8QGgUhWVXPuu9sbqwgRJZt27axYMECVq1ahapSrlw5WrRoQYsWLShbtqzreCHj0KFDLF++nMWLF7N7t/8Wzg0aNKBNmzbUqFHDcToTDqwgmJCzZ88ePvvsM1atWgVAlSpVuP7662nYsCExMXbCW2H5fD5Wr17N119/zaZNm1BVGjVqxI033kilSpVcxzMhKFwKwjpgL6DAy6o6NZdlegI9AWrWrNliwwa7Q2O48Pl8pKWl8fnnn3Pw4EEqVKjATTfdROPGjW0QtRipKunp6Xz22Wfs2rWLhIQEbrjhBq6++mouuMCGCU34FISqqrpVRBKBuUAvVV2Q1/J2hBD6jh07xrx58/jmm28QEVq0aEG7du2sCyiIDh8+zJdffsnChQvx+Xy0atWKG2+8kfj4eNfRjCNhURByEpGRwCFVnZDXMlYQQtORI0f49NNPWbx4MXFxcbRr147WrVuHzJk/0czn87Fw4ULmzp1LZmYmzZo1o2PHjpQuXdp1NBNEIV8QRKQUEKOqB73/zwVGq+qneb3GCkLoyMrKYv78+cyfP58LL7yQDh060LJlSxsLCGGqyvfff8/s2bM5fPgw119/PTfeeKN1K0WBcCgIdYCZ3rcXAP9Q1afze40VBPdWrlzJzJkzyczMpG3btrRt29aOBMKQz+djwYIFzJ07lwsuuIBOnTrRokULG9uJUCFfEM6FFQQ3Dh06xLvvvsuaNWto0qQJd9xxh3U5RJDMzExmz57NkiVLqFu3Ll27drUxnwhjBcGct5UrVzJ9+nRiY2O5++67ueKKK1xHMgH2008/8e6773L06FHuuOMOkpML3IaYMGAFwZyTkydPMmPGDBYvXkyTJk3o0qWLTcwWhY4fP86sWbP47rvvaN68OV26dKFEiRKuY5lzZAXBFMnBgwd5/fXX2bRpE3feeSdXX32160gmRCxZsoQPPviApKQkunfvTvny5V1HMkVU2IJgpxdEuc2bN/P3v/+drKws7r33Xi655BLXkUyISU5OJjk5mc2bNzNlyhROnDhBjx49qFWrlutoppjZEUKU2rBhA6+88grly5enR48ettdnCu3gwYO89tprZGRk0KNHD+rWres6kimAdRmZXK1bt45XX32VSpUq8cADD1CmTBnXkUyYOnLkCNOmTWPLli3ce++9dtJBCLOCYE6zYcMGXn75ZZKSknjggQcoVaqU60gmQhw9epQ33niDn3/+mZ49e9oRQwiyMQQDwO7du5kyZQqlSpVi6NChVghMsYuPj+ehhx4iMzOTqVOnsnv3bh577DGSkpJcRzNFZAUhQh05coSXXnqJAwcO8Pjjj9u0yCbgSpYsSZ8+fdi7dy8vvvgisbGxPProo3aRWxixLqMI4/P5eOutt/jhhx946KGHqFOnjutIJkpt2bKFKVOmULduXbp3725TnDhkYwhR6LvvvuOtt96iW7dudh2BCRlLly7ljTfe4K677qJ169au40QlKwhRZOfOnaSmplKvXj3uvfdem3HUhBxV5b333mP58uX07t2batWquY4UVawgRIGsrCymTp3K9u3b6devH+XKlXMdyZh8HTp0iMmTJ1OqVCkeffRRm3o7SApbEGxXMkytWLGCPn36cO211zJq1CgrBiYslC5dmsGDB9OuXTv69etHWlqa60gmBztCCDNHjx4lNTWVsmXL8vDDD9tAnQlbPp+P1157jW3btvHEE0/YlOoBZNchRKCvv/6aGTNm0LdvX2rXru06jjHnJSYmhgcffJCtW7cyYsQIbr75Zm688UbXsaKadRmFgaNHjzJ69GjWr19PamqqFQMTUapWrcrEiRPZs2cPw4cP5/Dhw64jRS07Qghxy5YtY9q0afTv35+aNWu6jmNMwPz+979n69atDBkyhLvvvptrrrnGdaSoY2MIISorK4vnn3+euLg4HnnkETuV1EQNVWXatGns3LmTJ598kri4ONeRwp6dZRTGNm7cSJ8+fbjxxht57LHHrBiYqCIi3H///dx111088cQT/Pzzz64jRQ3rMgoxH3/8MYsWLWLChAnEx8e7jmOMM3Xr1uXZZ59l/Pjx1K9fn86dO7uOFPFs1zNEnDhxgtGjR5OZmcmYMWOsGBgDxMXFMWzYMOLj4xk2bBhHjx51HSmiWUEIARs2bKBfv37cc889dOnSxXUcY0JOx44d6dmzJykpKaxdu9Z1nIjlvCCISKyILBORj11ncWHevHlMnTqVCRMmcOmll7qOY0zIqlGjBqmpqbzzzjt88sknruNEJOcFAegD/Og6RLCpKi+++CLr16/n6aefti4iYwqhRIkSjBw5kt27d/Pss88STmdJhgOnBUFEqgO3AK+6zBFsmZmZDBo0iCZNmnD//fe7jmNM2OnWrRvXX389Tz75JAcPHnQdJ2K4PkJ4DhgA+PJaQER6isgSEVmSkZERvGQBsnnzZlJSUujVqxfXXXed6zjGhK0WLVowaNAgBg0axLp161zHiQjOCoKIdAR2qmq+0x2q6lRVTVbV5MqVKwcpXWAsX76c1NRUJk6cSPXq1V3HMSbsJSYmkpqayksvvcSiRYtcxwl7Lo8QWgO3ish64D2grYi87TBPQH366afMmTOHZ555xsYLjClGcXFxjBs3jq+//pqZM2e6jhPWnBUEVR2sqtVVtTZwNzBfVbu5yhNI06ZNOzVHi111bEzxExFSUlI4duwYL7zwgg02nyPbOgWQqvLnP/+Ziy++mB49eriOY0zEu/vuu2natCkjR47E58tzaNLkISQKgqp+qaodXecoTllZWQwePJibb76ZDh06uI5jTNRo06YNXbp0YeDAgZw4ccJ1nLASEgUh0hw7doz+/fvTvXt3mjdv7jqOMVGnUaNGPPLII6SkpHDkyBHXccKGFYRidvjwYVJSUujTpw9XXHGF6zjGRK06deowaNAg+vfvz759+1zHCQtWEIrRvn37GDBgAEOHDqVWrVqu4xgT9S6++GLGjh3L0KFD2bVrl+s4Ic8KQjHZt28fQ4YMYezYsVSpUsV1HGOMp3z58owbN44RI0YQCRe3BpIVhGKQXQyefvppypcv7zqOMeYMZcqUYdy4cYwcOdKKQj6sIJyn/fv3M3ToUCsGxoQ4KwoFs4JwHg4cOMCQIUMYM2aMFQNjwkCZMmUYP368FYU8WEE4R4cPH2bw4MGMHj2aChUquI5jjCmk0qVLnzpSOHDggOs4IcUKwjk4fvw4gwcPZvjw4VSsWNF1HGNMEZUpU4bRo0czePBgMjMzXccJGVYQisjn8zF48GD69etHUlKS6zjGmHNUsWJFhgwZwuDBg+2KZo8VhCJQVUaMGEH37t255JJLXMcxxpynatWq0atXL4YOHWpzH2EFoUgmTJhAhw4daNy4sesoxphiUrduXbp168bYsWNdR3HOCkIhvfnmm9SvX59rr73WdRRjTDFr0qQJrVu3ZurUqa6jOGUFoRDmz5/PkSNH6NgxoiZkNcbk0K5dO+Lj4/n4449dR3HGCkIBVq9ezZdffsnDDz/sOooxJsD+9Kc/sWrVKpYuXeo6ihNWEPKRkZHB3/72N4YPH+46ijEmSAYOHMi7777Lxo0bXUcJOisIeTh27BgjR45k7NixXHDBBa7jGGOCREQYO3Ys48ePj7prFKwg5GHMmDH079+fMmXKuI5ijAmyCy+8kCFDhjBixIiouj+zFYRcvP766/z617+mdu3arqMYYxypVq0at99+Oy+++KLrKEFjBeEMixYt4sCBA7Rv3951FGOMY9deey0lS5Zk3rx5rqMEhRWEHHbs2MH7779Pr169XEcxxoSIHj16MH/+fDZs2OA6SsBZQfBkZWUxduxYRo0ahYi4jmOMCSHDhg3jmWeeifg5j6wgeFJTU3nooYcoVaqU6yjGmBATHx9P3759GT9+vOsoAeWsIIhIvIh8JyLfi0i6iIxyleWrr76iQoUKNGrUyFUEY0yIu/TSS7nsssuYM2eO6ygB4/II4RjQVlWbAs2Am0WkVbBD7Nq1i9mzZ3PfffcFu2ljTJj5/e9/zzfffMPmzZtdRwkIZwVB/Q5535bwvoJ6wq+qMnbsWIYNG2bjBsaYQnnqqacYP348WVlZrqMUO6djCCISKyLLgZ3AXFVdlMsyPUVkiYgsKe57oL7yyit07dqViy66qFjXa4yJXCVLluTBBx9kypQprqMUO6cFQVWzVLUZUB24SkTO6sRX1amqmqyqyZUrVy62tn/66Sf27NnD1VdfXWzrNMZEhyZNmuDz+Vi5cqXrKMUqJM4yUtV9wJfAzcFoLysri+eee44nnngiGM0ZYyJQr169mDp1KsePH3cdpdi4PMuosoiU8/5fEmgPrA5G2y+88AIPPfQQcXFxwWjOGBOBYmNj6du3LxMnTnQdpdi4PEK4GPhCRFYAi/GPIQT8zhSrVq1CVe02mMaY81a3bl0SExNZuHCh6yjFwtm8zqq6ArgymG1mZWXx8ssv89xzzwWzWWNMBOvRowd9+vShefPmYd/rEBJjCMHy8ssv8+CDDxIbG+s6ijEmQogIjz/+OJMnT3Yd5bxFTUHYtGkTe/fupUmTJq6jGGMiTL169QBYs2aN4yTnJyoKgqry7LPP0q9fP9dRjDERqlevXkyePDmsb6gTFQXhww8/pEOHDiQkJLiOYoyJUHFxcfzhD3/gjTfecB3lnEXkzYIXLlzIpEmTWLNmDXXq1CEuLo533nnHdSxjTIRr1aoVH330Efv27aNcuXKu4xRZxB0hjBgxgnbt2vHPf/6TpUuXMmPGDGbOnMmIESNcRzPGRIG+ffuG7ZmMEVUQFi5cyIQJEzhy5MipfjxVJTMzkwkTJkTMucLGmNCVmJhI5cqVSU9Pdx2lyCKqIEyaNInMzMxcnzt69CiTJk0KciJjTDTq2bMnr7zyStgNMEdUQVizZk2eH4DP52Pt2rVBTmSMiUYlSpSgY8eOzJo1y3WUIilUQRCR8iLSUETqiEjIFpF69eoRE5N7vJiYmFPnChtjTKC1b9+eb775Jqzuw5znxl1ELhKRISKyElgIvAxMBzaIyPsickOwQhZWnz59iI+Pz/W5+Ph4evfuHeRExpho1r17d/7+97+7jlFo+e3tfwBsAtqo6uWqep13X4IawDjgNhG5PygpC6lVq1akpKSQkJBw6kghJiaGhIQEUlJSaNUq6HfoNMZEsUaNGrFu3ToOHTpU8MIhQMJp0CM5OVmXLFlS4HLZ1yGsXbuWyy67jD59+lgxMMY4sXHjRqZPn05KSoqzDCKSpqrJBS2X54VpItL8jIcU2KWqm843XKC1atXKCoAxJiTUrFmTI0eOsGvXLipVquQ6Tr7yu1I5t7s+VBCROKCrqi4PUCZjjIkojz32GC+++CLDhg1zHSVfeRYEVc110FhEkoHngV8FKpQxxkSSihUrEhsbS0ZGBsV5b/jiVuRTSFV1CVA6AFmMMSZiZV+sFsqKXBBEJAn/eIIxxphCqlSpEllZWezZs8d1lDzlN6g8mbM3/BWAa4E+gQxljDGRqGfPnkydOpVBgwa5jpKr/AaVzzy/U4HdwBOqujNwkYwxJjIlJSVx7Ngx9u/fz0UXXeQ6zlny7DJS1TfO+HpTVf8NXCgi/YOY0RhjIkaPHj14/fXXXcfIVWHnMqokIo+IyALgSyApoKmMMSZC1ahRg23btnHy5EnXUc6S31xGZUTkTyLyKfAdcClQR1Xrqqq7S+6MMSbMde7cmQ8++MB1jLPkd4SwE7gfeBqoq6pPAseLq2ERqSEiX4jIjyKSLiI2UG2MiQrJycksWbIk5O6XkF9BGALEA38DBotI3WJu+yTwpKrWB1oBj4lIg2JuwxhjQlKbNm345ptvXMc4TX6DyqmqejVwKyDALKCqiAwUkfO+sYCqblPVpd7/DwI/AtXOd73GGBMOOnbsyMcff+w6xmkKHFRW1Z9V9WlVbQy0BC4CPinOECJSG7gSWJTLcz1FZImILMnIyCjOZo0xxpnY2FgSExPZvn276yin5DeoLGc+pqorVXWIqtbNa5miEpHSwAygr6oeyKXNqd59GJJDeQ4QY4wpqm7duvH222+7jnFKfkcIX4hILxGpmfNBEYkTkbYi8gZw7/k0LiIl8BeDd1T1w/NZlzHGhJvExER27dpFVlaW6yhA/gXhZiALeFdEtorIDyLyM7AW6Aqkqurr59qwd3TxGvCjqj57rusxxphwdtNNN/HZZ5+5jgHkP6h8VFVfVNXWQC2gHdBcVWup6oPFcD+E1sAfgbYistz7+u15rtMYY8LKDTfcwBdffOE6BpD/XEanqOoJYFtxNqyq/4v/7CVjjIlaIkJiYmJI3CuhyNNfG2OMKV533XUX06dPdx3DCoIxxrhWs2ZNNm7c6DpGwQVBRB4XkfLBCGOMMdHqsssuY82aNU4zFOYIoQqwWESmi8jNxXHtgTHGmNOFwoR3hblS+SngMvyniHYH1orInwMwt5ExxkStcuXKceDAAacT3hVqDEH9Cbd7XyeB8sAHIvLXAGYzxpio0rx5c5YvP98z+s9dYcYQeotIGvBX4Bugsao+ArQA7gxwPmOMiRodOnRgzpw5ztovzHUIlYDfqeqGnA+qqk9EOgYmljHGRJ8yZcpw6NAhZ+0XZgxh+JnFIMdzPxZ/JGOMiV6XX345q1evdtK2XYdgjDEhpFOnTsyePdtJ21YQjDEmhFSsWJE9e/Y4adsKgjHGhJhKlSqxa9euoLdrBcEYY0LMTTfdxOeffx70dq0gGGNMiGnUqBErV64MertWEIwxJsRkzxAU7KuWrSAYY0wIatq0KStWrAhqm1YQjDEmBLVr14558+YFtU0rCMYYE4IqVqzI3r17g9qmFQRjjAlRIoLP5wtae1YQjDEmRDVs2JD09PSgtWcFwRhjQlSbNm1YsGBB0NqzgmCMMSGqatWqbNu2LWjtWUEwxpgQFsy7FjstCCIyTUR2isgqlzmMMSZUJSYmsmPHjqC05foI4XXgZscZjDEmZLVs2ZK0tLSgtFWYO6YFjKouEJHaLjOY6KKqfP/990yfPp21a9dy9OhRLrvsMjp16sT1119PTIzrfSRjTte0aVNSU1P57W9/G/C2nBaEwhCRnkBPgJo1azpOY8LJ5s2bqVGjRqGWTU1NzfXxBQsW0KZNm+KMZUyRlCxZkqNHjwalrZDfHVLVqaqarKrJlStXdh3HhLhly5YhIojIqWIwe/ZsVLVIX9lzyPzqV786tb7XXnvN5VszUSxYk9yFfEEwpjBuv/12RITmzZtTvnx5fD7fqY17x44di7y+xo0bn1Yg7rvvPh544IFTxSHYs1Ca6JaUlERGRkbA27GCYMLaXXfdhYjwr3/9i0WLFqGq7Nmzp9hP1Zs2bRqqyu7duwGIiYkJ6umAJro1aNCAH374IeDtuD7t9F3gW+ByEdksIve7zGPCR1paGiLC+++/T3p6OqrKVVddFfB2K1SogKpy+PBhwH+O+MSJEwPeroluDRo0CMoUFk4Lgqp2VdWLVbWEqlZXVeukNQUSEZKTkxkwYACqSoMGDYKeISEhAVXl888/JyUlxbqRTEAlJSUF5VqEkD/LyJhsx44dIz4+HoCsrKyQOEW0Xbt2qCoiQkxMDMuXL6dp06auY5kIE6zuSfd/UcYUQlpaGvHx8ZQqVQpVDYlikJOqMnz4cJo1a8aYMWNcxzHmnITWX5UxuZgxYwbJycn85S90TOJ2AAAN60lEQVR/4dChQ67j5GnUqFEsX76c4cOH0759e9dxTIQpW7YsBw4cCGgbVhBMSJsxYwadO3dm5syZDBo0yHWcAjVt2pTdu3czb948brrpJtdxTAS55JJLWLduXUDbsIJgQtbKlSvp3Lkzs2bN4vbbb3cdp9AqVKjArl27mDt3blgUMRMeateubQXBRKcTJ07QpEkThg8fzm233eY6TpFVrFiRH374gfHjx7N48WLXcUwEuOSSS1i/fn1A27CCYEJSXFwc4O+XD1f169dn3LhxXHXVVXZKqjlv5cuXZ9++fQFtwwqCCTnZp9hFwkZ04MCBACF3VpQJT4H+m7DfUhNS1qxZA/i7jCJF9h9xnz59HCcxJn9WEExIufzyy7nnnnu44ILIumZy2bJlPP/8865jGJMvKwgmZKSkpADwzjvvOE5S/Jo1awYE9/64JvKICD6fL2Drt4JgQsbEiRP56quvXMcImKysLICA/kGbyFahQgX27NkTsPVbQTAhYcCAAYD/hjSRKntgOTY21nESE66qVKkS0EnurCCYkPDMM8/wwQcfuI4RcMePH3cdwYSxpKQktm/fHrD1W0Ewzm3duhWAO++803GSwCtRogQAvXv3dpzEhKOKFSuyd+/egK3fCoJxrlq1aq4jBNWrr77K5MmTXccwYahs2bLs378/YOu3gmBCQvYdyKLB/ffbjQHNuQn0jKdWEIxT2X3qCQkJjpME37hx41xHMGGmTJkyVhBM5Mo+Pz/aNGzYkMGDB7uOYcJMbGysXYdgItePP/5Iz549XccIuiVLlriOYMxZrCAY5/72t7+5jhB02feGNqaoAnm1uxUE41ywZwKdP38+LVs2olQpISZGKFVKaNmyEfPnzw9qDmNCjRUEE1VGjx5Np07tuOGGdFasgOPHYcUKuOGGdDp1asfo0aODkmPhwoUAXHHFFXTt2vXU98Y4parOvoCbgf8DfgIGFbR8ixYt1EQOn8+n/l/B4Jg3b54mJKD/+U/uv2L/+Q+akIDOmzcvoDmGDx+uCQkJCiigMTExmpCQoMOHDw9ouyYyjBgxosivAZZoIbbJzo4QRCQWmAJ0ABoAXUWkgas8Jvg2bdoU1PYGDuzNY4/BNdfk/vw118Cjj8KgQYG7b8HChQuZMGECR44cOfWYz+fjyJEjTJgwwY4UjFMuu4yuAn5S1Z9V9TjwHhB+N8815yz7ZjjB8sMP6Tz0UP7LPPwwpKevCliGSZMmkZmZmetzR48eZdKkSQFr25iCuCwI1YCcu4ibvcdOIyI9RWSJiCzJyMgIWjgTeMG+CU5mJtSqlf8yNWv6lwuUNWvW5HkbRJ/Px9q1awPXuIkIf/zjHwO2bpcFIbdzp876S1HVqaqarKrJlStXDkIsEyx16tQJanslS8KGDfkvs3Gjf7lAqVevXp5nVcXExFCvXr3ANW4iQt26dQO2bpcFYTNQI8f31YGtjrIYB2rUqFHwQsWoQYOGvPxy/su89BI0bNgoYBn69OmT5zUI8fHxNguqccplQVgMXCYil4hIHHA38JHDPCbIgn07yfHjn2fKFPj229yf//ZbePFFGDcucP34rVq1IiUl5bS5m2JiYkhISCAlJYVWrVoFrG1jCuKsIKjqSeBx4H+AH4HpqpruKo+JfG3btmXgwFG0bw/9+8N//wsnTvj/7d8f2reHgQNH0bZt24DmGDVqFPPmzQOgRYsW3HXXXcybN49Ro0YFtF1jCiJ5DXCFouTkZLU5YCKLiHDy5Mmg3lZy/vz5DBrUh/T0VWRm+scMGjZsxLhxkwJeDHISkTwHmI0pTiKSpqrJBS0X3NM8jMnFY489xksvvRS09tq2bct3360MWnu5yXkdgjGhwqauME7Vq1ePlwsa6Y1AV1xxhesIxpzFCoJxauVKt3vqrmzatIlnn33WdQxjTmMFwTgVFxcHkOfVu5GsX79+riMYcxorCCYkRNMtNCdPnuw6gjG5soJgnNu+fbvrCEHVu3dv+vbt6zqGMWexgmCcS0pKAuDDDz90nCTwjh07BkBqaqrjJMaczQqCCQkDBw7kzjvvdB0j4OzWmSaUWUEwIWHcuHEAfP31146TBM7JkycByMrKcpzEmNxZQTAhIyUlhV/96leuYwRMiRIlgODfQ9qYwrLfTBMynnnmGQDuvfdex0mKX/aUKz6fz3ESY/JmBcGElNWrV/Pmm29y4sQJ11GKVcuWLRk0aFDQZ3g1pihscjsTcrI3muH0u5mfSHs/JvwUdnI7O0IwISd7wxkJe9OjR48GrKvIhAcrCCYkHT9+HPhlgxqOVq1axYgRI1i6dGlEFDcT+Wz6axOSSpQowfLly2nWrBmNGzfmjjvucB2pSDIyMmjcuDHDhw/nyiuvdB3HmEKxIwQTspo2bcrMmTP53e9+x+zZs13HKbSdO3eSmJjILbfcYndBM2HFCoIJabfffjvvv/8+t956K08//bTrOAVatmwZSUlJdOrUiY8//th1HGOKxAqCCXmdO3dm2bJlPPXUU1x00UWu4+RpwIABNG/enDFjxvDRRx+5jmNMkdkYggkLzZo1IzMzk5IlSyIiZGVlhdQVv9mDxunp6TRo0MBxGmPOTej8RRlTgPj4+FOnpMbGxpKSkuI4EXz22WenioHP57NiYMKaFQQTdlSVtLQ0Jk6ciIiwatWqoGc4dOgQIsJvfvMbnnvuOVTVTi01Yc+6jExYat68OapK586dady4MQCLFi3iqquuCmi7GRkZJCYmnvrerj42kcTJEYKIdBGRdBHxiUiBl1Mbk5cPPvgAVaVjx45cffXViAjly5cv9g11165dEZFTxUBVrRiYiOOqy2gV8DtggaP2TYSZPXv2qa6kffv2ERMTg4ggIsyaNavI61u5cuWp14sI7733Hm+99ZYVAhPRnHQZqeqPEBlz1ZjQkt2VBLB582Zq1KhxXlc5B6MbyphQEfJjCCLSE+gJULNmTcdpTDipXr36WXvzqkp6ejoffvghK1asIDY2lpo1a3LbbbfRunVr20kxUS1gBUFEPgeq5PLUUFX9V2HXo6pTgangn/66mOKZKCUiNGrUiEaNGrmOYkzICVhBUNX2gVq3McaY4mfXIRhjjAHcnXZ6h4hsBq4B/i0i/+MihzHGmF+4OstoJjDTRdvGGGNyZ11GxhhjACsIxhhjPFYQjDHGAFYQjDHGeCSc5mURkQxgQzGvthKwq5jXGU7s/Ufv+4/m9w7R9f5rqWrlghYKq4IQCCKyRFWjdsZVe//R+/6j+b2Dvf/cWJeRMcYYwAqCMcYYjxUEb+K8KGbvP3pF83sHe/9nifoxBGOMMX52hGCMMQawgmCMMcZjBSEHEUkRERWRSq6zBIuIPCMiq0VkhYjMFJFyrjMFg4jcLCL/JyI/icgg13mCSURqiMgXIvKjiKSLSB/XmVwQkVgRWSYiH7vOEiqsIHhEpAZwI7DRdZYgmws0UtUmwBpgsOM8ASciscAUoAPQAOgqIg3cpgqqk8CTqlofaAU8FmXvP1sf4EfXIUKJFYRfpAIDgKgaZVfVz1T1pPftQqC6yzxBchXwk6r+rKrHgfeA2xxnChpV3aaqS73/H8S/UazmNlVwiUh14BbgVddZQokVBEBEbgW2qOr3rrM41gP4xHWIIKgGbMrx/WaibIOYTURqA1cCi9wmCbrn8O8A+lwHCSVObpDjgoh8DlTJ5amhwBDgpuAmCp783ruq/stbZij+roR3gpnNEcnlsag6MgQQkdLADKCvqh5wnSdYRKQjsFNV00Tk167zhJKoKQiq2j63x0WkMXAJ8L2IgL/LZKmIXKWq24MYMWDyeu/ZROReoCPQTqPjwpTNQI0c31cHtjrK4oSIlMBfDN5R1Q9d5wmy1sCtIvJbIB4oKyJvq2o3x7mcswvTziAi64FkVY2KWRBF5GbgWeB6Vc1wnScYROQC/APo7YAtwGLgHlVNdxosSMS/5/MGsEdV+7rO45J3hJCiqh1dZwkFNoZgXgDKAHNFZLmIvOQ6UKB5g+iPA/+Df0B1erQUA09r4I9AW+8zX+7tLZsoZ0cIxhhjADtCMMYY47GCYIwxBrCCYIwxxmMFwRhjDGAFwRhjjMcKgjFFJCIlReQrb5K87Mf6ichREbkox2PdReSFM177pYgke///XETKBy+5MfmzgmBM0fUAPlTVrByPdcV/gdsdRVjPW8CjxRnMmPNhBcEYj4i09O4LES8ipbx7BTTKZdE/AP/K8bq6QGngKfyFobA+KuLyxgRU1MxlZExBVHWxiHwEjAVKAm+r6qqcy4hIHFBHVdfneLgr8C7wNXC5iCSq6s5CtLdXRC4UkYqqurvY3ogx58iOEIw53Wj8N0pKBv6ay/OVgH1nPHY38J6q+oAPgS7e43lNA5Dz8Z1A1XNOa0wxsiMEY05XAX/3Twn8M2EePuP5TO9xAESkCXAZ/rmgAOKAn/HfkW03cOagcQUg58SJ8d46jXHOjhCMOd1UYBj++0KMP/NJVd0LxIpIdlHoCoxU1dreV1WgmojUwj/I3FpEqgB4ZxddiHdzHm/W0SrA+sC+JWMKx44QjPGIyJ+Ak6r6D++U0v+ISFtVnX/Gop8B1wGf4+8u6nDG8zOBu1V1vHcD+zkiEgMcArp6XUsALYCFOW5haoxTNtupMUUkIlcCT6jqH89zPZOAj1R1XvEkM+b8WJeRMUWkqsuAL3JemHaOVlkxMKHEjhCMMcYAdoRgjDHGYwXBGGMMYAXBGGOMxwqCMcYYwAqCMcYYz/8D/TmGTIfHIBwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(figsize=(6, 4))\n",
    "\n",
    "p1, p2, t = orbit2(tmax=5, \n",
    "                  planets=[data['Earth'], data['Jupiter']])\n",
    "\n",
    "plot_orbit(p1.x, p1.y, t)\n",
    "plot_orbit(p2.x, p2.y, t)\n",
    "plt.axis('equal')\n",
    "\n",
    "plt.title('Earth, Jupiter orbit around the Sun')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Another option for animation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's build a little widget to explore how Earth and Jupiter move relative to each other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ipywidgets import interact"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ebd3003b15084dfaa00b826dc0db2047",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='t0', max=12.0), Output()), _dom_classes=('widget-int…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "tmax = 12.0\n",
    "\n",
    "p1, p2, t = orbit2(tmax=tmax,\n",
    "                  planets=[data['Earth'], data['Jupiter']])\n",
    "\n",
    "@interact(t0=(0, tmax))\n",
    "def investigate_orbit(t0=0):\n",
    "    plot_orbit(p1.x,p1.y, t, t0=t0)\n",
    "    plot_orbit(p2.x,p2.y, t, t0=t0)\n",
    "    plt.title('Earth, Jupiter orbit around the Sun')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Precession of Earth's Orbit due to a very large Jupiter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now observe what happen if we increase the mass of Jupiter to 300 $M_J$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7ed14a85949842d59a3446970a992577",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='t0', max=4.0), Output()), _dom_classes=('widget-inte…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "scale = 300\n",
    "tmax = 4.0\n",
    "\n",
    "jupiter = data['Jupiter']\n",
    "big_jupiter = jupiter.copy()\n",
    "big_jupiter['m'] = jupiter['m'] * scale\n",
    "\n",
    "p1, p2, t = orbit2(tmax=tmax,\n",
    "                  planets=[data['Earth'], big_jupiter])\n",
    "    \n",
    "@interact(t0=(0, tmax))\n",
    "def investigate_orbit(t0=0):\n",
    "    \n",
    "    plot_orbit(p1.x,p1.y, t, t0=t0)\n",
    "    plot_orbit(p2.x,p2.y, t, t0=t0)\n",
    "            \n",
    "    plt.suptitle(\"Effect of a very large Jupiter on Earth's orbit\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should notice that the orbit is now precessing. \n",
    "\n",
    "We can explore how the mass of Jupiter changes the rate of the precession."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "543e552a1a524a6f94611f7ed6de4c0e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=0, description='scale', max=1000, step=50), Output()), _dom_classes=('wi…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "tmax = 5.0\n",
    "\n",
    "@interact(scale=(0, 1000, 50))\n",
    "def investigate_precession(scale=0):\n",
    "    \n",
    "    # solve for the orbits of Earth and a large Jupiter\n",
    "    jupiter = data['Jupiter']\n",
    "    big_jupiter = jupiter.copy()\n",
    "    big_jupiter['m'] = jupiter['m'] * scale\n",
    "\n",
    "    p1, p2, t = orbit2(tmax=tmax,\n",
    "                  planets=[data['Earth'], big_jupiter])\n",
    "    \n",
    "\n",
    "    plot_orbit(p1.x, p1.y, t)\n",
    "    plot_orbit(p2.x, p2.y, t)\n",
    "    \n",
    "    \n",
    "    plt.suptitle(\"Effect of a large Jupiter on Earth's orbit\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that for very large Jupiters, the Earth's orbit is unstable.  At a mass of 1000 $M_J$, the Earth is actually ejected from the solar system. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that, here, we are ignoring that fact that the Sun itself should really be moving, too.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But first, let's see if we can get all of the planets into one simulation.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Mass of the Sun in kg\n",
    "M_s = 2e30\n",
    "\n",
    "def orbit3(tmax = 1.0, dt = 0.001, planets=[]):\n",
    "    \"\"\"\n",
    "    Calculate the orbit of many planets\n",
    "    \"\"\"\n",
    "\n",
    "    # number of timesteps\n",
    "    N = round(tmax/dt)\n",
    "    t = np.zeros(N)\n",
    "\n",
    "    # allocate planets\n",
    "    all_p = [Planet(data, N) for data in planets]\n",
    "    # this creates a list of all the planets\n",
    "    \n",
    "    def F(q):\n",
    "        dqdt = np.zeros(len(q))\n",
    "        \n",
    "        # for each planet\n",
    "        for i in range(len(q)//4):\n",
    "            x1, y1, vx1, vy1 = q[4*i:4*i+4]\n",
    "            r1 = np.sqrt(x1**2+y1**2)\n",
    "            dxdt = vx1\n",
    "            dydt = vy1\n",
    "            dvxdt = - 4*π**2*x1/(r1**3)\n",
    "            dvydt = - 4*π**2*y1/(r1**3)\n",
    "            \n",
    "            # for each other planet\n",
    "            for j in range(len(q)//4):\n",
    "                if i==j:\n",
    "                    continue\n",
    "                x2, y2, vx2, vy2 = q[4*j:4*j+4]\n",
    "                r12 = np.sqrt((x1-x2)**2+(y1-y2)**2)\n",
    "\n",
    "                dvxdt += - 4*π**2*(all_p[j].m/M_s)*(x1-x2)/(r12**3)\n",
    "                dvydt += - 4*π**2*(all_p[j].m/M_s)*(y1-y2)/(r12**3)\n",
    "            \n",
    "            dqdt[4*i:4*i+4] = [dxdt, dydt, dvxdt, dvydt]\n",
    "                         \n",
    "        return dqdt\n",
    "        \n",
    "    # loop over the timesteps\n",
    "    for i in range(N-1):\n",
    "        t[i+1] = t[i] + dt\n",
    "        \n",
    "        q0 = []\n",
    "        for p in all_p:\n",
    "            q0 += [p.x[i], p.y[i], p.vx[i], p.vy[i]]\n",
    "        q0 = np.array(q0)\n",
    "                    \n",
    "        s1 = F(q0)\n",
    "        q1 = q0 + s1*dt/2\n",
    "        s2 = F(q1)\n",
    "        q2 = q0 + s2*dt/2\n",
    "        s3 = F(q2)\n",
    "        q3 = q0 + s3*dt\n",
    "        s4 = F(q3)\n",
    "        q4 = q0 + (s1 + 2*s2 + 2*s3 + s4)/6 * dt\n",
    "        \n",
    "        for n, p in enumerate(all_p):\n",
    "            p.x[i+1], p.y[i+1], p.vx[i+1], p.vy[i+1] = q4[4*n:4*n+4]\n",
    "        \n",
    "    return all_p, t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmax = 6.0\n",
    "\n",
    "p, t = orbit3(tmax=tmax,\n",
    "            # planets=[data['Earth'], data['Jupiter'], data['Mars']]\n",
    "            # planets=[data['Earth'], big_jupiter, data['Mars']]\n",
    "            # planets=[p for p in data.values()]\n",
    "             planets=[p for p in data.values() if p['a'] < 8]\n",
    "             )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6f88b19814964e9bad0546ba1cd62546",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='t0', max=6.0), Output()), _dom_classes=('widget-inte…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@interact(t0=(0, tmax, 0.1))\n",
    "def investigate_orbit(t0=0):\n",
    "    for p1 in p:\n",
    "        plot_orbit(p1.x,p1.y, t, t0=t0)\n",
    "    \n",
    "    plt.title('Solar System')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Three-Body Problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the most general case, the Sun itself would move do to the gravitational forces from the other planets.  Since the Sun is so much more massive than the other planets, it is reasonable to assume that the Sun does not itself move."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, suppose we had a binary star system with two stars and a planet. Then all three bodies would be in motion. Let's rewrite our last route so that the Sun is no longer special but also interacts with all the other bodies in the solar system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Mass of the Sun in kg\n",
    "M_s = 2e30\n",
    "\n",
    "def orbit4(tmax = 1.0, dt = 0.001, bodies=[]):\n",
    "    \"\"\"\n",
    "    Calculate the orbit of many planets\n",
    "    \"\"\"\n",
    "\n",
    "    # number of timesteps\n",
    "    N = round(tmax/dt)\n",
    "    t = np.zeros(N)\n",
    "\n",
    "    # allocate planets\n",
    "    all_p = [Planet(data, N) for data in bodies]\n",
    "    # this creates a list of all the planets\n",
    "    \n",
    "    def F(q):\n",
    "        dqdt = np.zeros(len(q))\n",
    "        \n",
    "        # for each body in the system\n",
    "        for i in range(len(q)//4):\n",
    "            x1, y1, vx1, vy1 = q[4*i:4*i+4]\n",
    "            r1 = np.sqrt(x1**2+y1**2)\n",
    "            dxdt = vx1\n",
    "            dydt = vy1\n",
    "            dvxdt = 0\n",
    "            dvydt = 0\n",
    "            \n",
    "            # for each other body\n",
    "            for j in range(len(q)//4):\n",
    "                if i==j:\n",
    "                    continue\n",
    "                x2, y2, vx2, vy2 = q[4*j:4*j+4]\n",
    "                r12 = np.sqrt((x1-x2)**2+(y1-y2)**2)\n",
    "\n",
    "                dvxdt += - 4*π**2*(all_p[j].m/M_s)*(x1-x2)/(r12**3)\n",
    "                dvydt += - 4*π**2*(all_p[j].m/M_s)*(y1-y2)/(r12**3)\n",
    "            \n",
    "            dqdt[4*i:4*i+4] = [dxdt, dydt, dvxdt, dvydt]\n",
    "                         \n",
    "        return dqdt\n",
    "        \n",
    "    # loop over the timesteps\n",
    "    for i in range(N-1):\n",
    "        t[i+1] = t[i] + dt\n",
    "        \n",
    "        q0 = []\n",
    "        for p in all_p:\n",
    "            q0 += [p.x[i], p.y[i], p.vx[i], p.vy[i]]\n",
    "        q0 = np.array(q0)\n",
    "                    \n",
    "        s1 = F(q0)\n",
    "        q1 = q0 + s1*dt/2\n",
    "        s2 = F(q1)\n",
    "        q2 = q0 + s2*dt/2\n",
    "        s3 = F(q2)\n",
    "        q3 = q0 + s3*dt\n",
    "        s4 = F(q3)\n",
    "        q4 = q0 + (s1 + 2*s2 + 2*s3 + s4)/6 * dt\n",
    "        \n",
    "        for n, p in enumerate(all_p):\n",
    "            p.x[i+1], p.y[i+1], p.vx[i+1], p.vy[i+1] = q4[4*n:4*n+4]\n",
    "        \n",
    "    return all_p, t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add the Sun as an explicit member of the solar system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "data['Sun'] = { 'name' : 'Sun', 'm' : 2e30, 'a' : 0, 'e' : 0}\n",
    "data['BigJupiter'] = { 'name' : 'BigJupiter', 'm' : 1.9e+27*200, 'a' : 5.2, 'e' : 0.048}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmax = 2.0\n",
    "\n",
    "p, t = orbit4(tmax=tmax,\n",
    "              bodies=[data['Sun'], \n",
    "                      #data['Jupiter'],\n",
    "                      data['BigJupiter'],  \n",
    "                      data['Earth']\n",
    "                     ]\n",
    "             )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "908759f789ef4f7a9ddfe9ff9809ab20",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='t0', max=2.0), Output()), _dom_classes=('widget-inte…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@interact(t0=(0, tmax, 0.1))\n",
    "def investigate_orbit(t0=0):\n",
    "    for p1 in p:\n",
    "        plot_orbit(p1.x,p1.y, t, t0=t0)\n",
    "    \n",
    "    plt.title('Solar System')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.4"
  },
  "widgets": {
   "state": {
    "138d8a5ff419441d9090e1fb1d52e934": {
     "views": [
      {
       "cell_index": 32
      }
     ]
    },
    "7b81ea171a7e48779c55c15018aca96b": {
     "views": [
      {
       "cell_index": 45
      }
     ]
    },
    "7c8cf3e5d0664257b9af7a7a64504a04": {
     "views": [
      {
       "cell_index": 29
      }
     ]
    },
    "b0fc71b2ceb44e858c121a91647d3876": {
     "views": [
      {
       "cell_index": 34
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}