Lecture.20-Precession-of-the-Perihelion-of-Mercury.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Precession of the Perihelion of Mercury\n",
    "## Lecture 20"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Follows Giordano 4.3: Precession of the Perihelion of Mercury**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Intro Video\n",
    "\n",
    "https://www.youtube.com/watch?v=hSXNE0pNtr8\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Kepler's First Law that says planets move in elliptical orbits with the Sun at one focus is a direct result of the inverse square law of Newton's Law of Gravitation:\n",
    "\n",
    "$$ F_G = \\frac{G M_s M_e}{r^2} $$\n",
    "\n",
    "If the law was of some other form, such as\n",
    "\n",
    "$$ F_G = \\frac{G M_s M_e}{r^\\beta} $$\n",
    "\n",
    "for same value of $\\beta$ different from 2, then the orbits of the planets would not be stable ellipses.  Section 4.2 of Giordano explores this stability.  It shows that even for value of $\\beta$ close 2, (like $\\beta = 2.01$) the orbits would not be stable and the orbit itself would tend to rotate over time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solar-system itself can be thought of as an experiment to test the accuracy of the inverse-square law.  In reality, the solar system is more than the Sun and single planet -- the perturbative effects of all the other planets need to be taken into account to correctly the expected deviations from Kepler's Laws.   One such deviation is that planets with eccentric orbits will *precess* over time. That is, the axes of the elliptical orbit itself rotate. We can be specific about this by discussing the rate at which the perihelion (nearest point to the Sun) precesses around the Sun.\n",
    "\n",
    "By the middle of 19th century, the corrections to the orbits had been calculated using perturbation theory to estimate the precession rate.  If all of the planets are included in the calculation, the precession rate of Mercury's orbit was estimated as **523 arcseconds/century**. This mostly due to Jupiter. We'll try and calculate this in our next lecture."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, astronomy is a precise science can careful observations measured the actual precession of the perihelion of Mercury was **566 arcseconds/century**. The perihelion rotates around once every 230,000 years.  Those are precise measurements!\n",
    "\n",
    "The difference between theory and experiment, a mere 43 arcseconds/century, was a mystery to 19th century physics. As mentioned in the video, the solution came with Einstein's theory of general relativity (GR) in 1917."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you may know, GR describes gravity not as force but as a curvature in space-time.  Provided that objects are not too close or too massive, GR makes a similar prediction for the force of gravity that of Newton's law.  But, if the objects are close enough, GR predicts deviations from the inverse-square law.  \n",
    "\n",
    "The Sun and Mercury are close enough for those deviations to be significant and can be used to explain the discrepancy of 43 arcseconds/century.  This was one of the first major successes of the theory of general relativity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our goal in  this lecture is calculate this 43 arcseconds/century deviation.  While it can be done analytical using GR, we can also do it numerically by extending our model for the calculation of an orbiting planet.  \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The magnitude of the force law predicted by GR is\n",
    "\n",
    "$$F_G \\approx \\frac{G M_S M_M}{r^2}\\left(1+\\frac{\\alpha}{r^2}\\right) $$\n",
    "\n",
    "where $M_M$ is the mass of Mercury and $\\alpha = 1.1 \\times 10^{-8}\\;\\mathrm{AU}^2$.  This value comes from a the first  perturbation expansion in powers of $r^{-1}$.  This deviation from an inverse-square law is what leads to the precession of the perihelion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One complication is that the precession rate for Mercury is so a naive simulation would need to be computed for 230,000 years or so to be able to directly measure the precession rate.  We would never want to wait for such a model to finish running! \n",
    "\n",
    "One strategy is use much larger values of $\\alpha$ and then predict the precession rate as a function of $\\alpha$.  We can then estimate the rate of precession for $\\alpha = 1.1 \\times 10^{-8}\\;\\mathrm{AU}^2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Equations of Motion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As before, we start with Newton's Second Law:\n",
    "\\begin{align}\n",
    "\\frac{\\mathrm{d}^2x}{\\mathrm{d}t^2} &= \\frac{F_{G,x}}{M_M}  \\\\\n",
    "\\frac{\\mathrm{d}^2y}{\\mathrm{d}t^2} &= \\frac{F_{G,y}}{M_M}\n",
    "\\end{align}\n",
    "\n",
    "which can be rewritten as a system of four first order differential equations:\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{\\mathrm{d}v_x}{\\mathrm{d}t} &= -\\frac{G M_s x}{r^3}\\left(1+\\frac{\\alpha}{r^2}\\right) \\\\\n",
    "\\frac{\\mathrm{d}x}{\\mathrm{d}t} &= v_x \\\\\n",
    "\\frac{\\mathrm{d}v_y}{\\mathrm{d}t} &= -\\frac{G M_s y}{r^3}\\left(1+\\frac{\\alpha}{r^2}\\right) \\\\\n",
    "\\frac{\\mathrm{d}y}{\\mathrm{d}t} &= v_y \n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The numerical equations using the Euler Cromer method are\n",
    " \n",
    "\\begin{align}\n",
    "v_{x,i+1} &= v_{x,i} - \\frac{4\\pi^2x_i}{r^3_i}\\Delta t\\left(1+\\frac{\\alpha}{r^2_i}\\right) \\\\\n",
    "x_{i+1} &= x_i + v_{x,i+1} \\Delta t  \\\\\n",
    "v_{y,i+1} &= v_{y,i} - \\frac{4\\pi^2y_i}{r^3_i}\\Delta t\\left(1+\\frac{\\alpha}{r^2_i}\\right) \\\\\n",
    "y_{i+1} &= y_i + v_{y,i+1} \\Delta t  \\\\\n",
    "\\end{align}\n",
    "\n",
    "Recall that we are measuring distance in astronomical units (AU) and time in years (yr)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " We could go ahead and code this, but what about if we chose to attack the problem using the Runge-Kutta 4 method. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "π = np.pi\n",
    "\n",
    "def orbit(tmax = 1.0, dt = 0.001, α = 0.0, v0 = 2*π, x0 = 1):\n",
    "    \"\"\"\n",
    "    Calculate the orbit of planet\n",
    "    α is the correction due to general relativity.\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(x, y, vx, vy):\n",
    "        r = np.sqrt(x**2+y**2)\n",
    "        dxdt = vx\n",
    "        dydt = vy\n",
    "        dvxdt = - (4*π**2*x)/(r**3)*(1+α/r**2)\n",
    "        dvydt = - (4*π**2*y)/(r**3)*(1+α/r**2)\n",
    "        \n",
    "        return [dxdt, dydt, dvxdt, dvydt]\n",
    "        \n",
    "    # loop over the timesteps\n",
    "    for i in range(N-1):\n",
    "        t[i+1] = t[i] + dt\n",
    "        \n",
    "        # estimate slope at initial point\n",
    "        sx1, sy1, svx1, svy1 = f(x[i], y[i], vx[i], vy[i])\n",
    "        \n",
    "        # estimate slope at mid point\n",
    "        x1  = x[i]   + sx1  * dt/2\n",
    "        y1  = y[i]   + sy1  * dt/2\n",
    "        vx1 = vx[i]  + svx1 * dt/2\n",
    "        vy1 = vy[i]  + svy1 * dt/2\n",
    "        sx2, sy2, svx2, svy2 = f(x1, y1, vx1, vy1)\n",
    "        \n",
    "        # improved estimate slope at mid point\n",
    "        x2 = x[i] + sx2 * dt/2\n",
    "        y2 = y[i] + sy2 * dt/2\n",
    "        vx2 = vx[i] + svx2 * dt/2\n",
    "        vy2 = vy[i] + svy2 * dt/2\n",
    "        sx3, sy3, svx3, svy3 = f(x2, y2, vx2, vy2)\n",
    "        \n",
    "        # estimate slope at final point\n",
    "        x3 = x[i] + sx3 * dt\n",
    "        y3 = y[i] + sy3 * dt\n",
    "        vx3 = vx[i] + svx3 * dt\n",
    "        vy3 = vy[i] + svy3 * dt\n",
    "        sx4, sy4, svx4, svy4 = f(x3, y3, vx3, vy3)    \n",
    "        \n",
    "        # estimate final point\n",
    "        x[i+1] = x[i] + (sx1 + 2*sx2 + 2*sx3 + sx4)/6 * dt\n",
    "        y[i+1] = y[i] + (sy1 + 2*sy2 + 2*sy3 + sy4)/6 * dt\n",
    "        vx[i+1] = vx[i] + (svx1 + 2*svx2 + 2*svx3 + svx4)/6 * dt\n",
    "        vy[i+1] = vy[i] + (svy1 + 2*svy2 + 2*svy3 + svy4)/6 * dt\n",
    "        \n",
    "    return x, y, vx, vy, t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Planetary Data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = { 'Mercury' : {'m':2.4e23, 'a':0.39,  'e':0.206},\n",
    "         'Venus'   : {'m':4.9e24, 'a':0.72,  'e':0.007},   \n",
    "         'Earth'   : {'m':6.0e24, 'a':1.00,  'e':0.017},  \n",
    "         'Mars'    : {'m':6.6e23, 'a':1.52,  'e':0.093},\n",
    "         'Jupiter' : {'m':1.9e27, 'a':5.20,  'e':0.048},\n",
    "         'Saturn'  : {'m':5.7e26, 'a':9.54,  'e':0.056},\n",
    "         'Uranus'  : {'m':8.8e25, 'a':19.19, 'e':0.046},\n",
    "         'Neptune' : {'m':1.03e26,'a':30.06, 'e':0.010},\n",
    "         'Pluto'   : {'m':1.3e22, 'a':39.53, 'e':0.248},\n",
    "        }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To compute the orbit of Mercury we need to know the correct velocity to use as an initial condition "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are some of the helper functions we introduced in [Lecture 19]()."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_orbit(x, y):\n",
    "    \"Make a plot of planet showing the Sun\"\n",
    "    \n",
    "    plt.plot(x, y, 'k-', linewidth=0.5)\n",
    "\n",
    "    # Plot the Sun at the origin\n",
    "    plt.plot(0, 0, '.',color=\"yellow\", markersize=40, markeredgecolor='black')\n",
    "    # Plot the Earth at the final position\n",
    "    plt.plot(x[-1], y[-1], '.', color= 'k', markersize=15)\n",
    "\n",
    "    plt.xlabel('x (AU)')\n",
    "    plt.ylabel('y (AU)')\n",
    "    plt.axis('equal')\n",
    "    \n",
    "def calc_eccentricity(x, y, t):\n",
    "    \"Calculate the eccentricity of an orbit\"\n",
    "    \n",
    "    # major semiaxis\n",
    "    a = (max(x)-min(x))/2\n",
    "    # minor semiaxis\n",
    "    b = (max(y)-min(y))/2\n",
    "    \n",
    "    #if b > a:\n",
    "    #    a, b = b, a\n",
    "\n",
    "    # compute eccentricity\n",
    "    e = np.sqrt(1 - b**2/a**2)\n",
    "    \n",
    "    return e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e72120869f3949c1bd2e11f4a4098c13",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>interactive</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "interactive(children=(FloatSlider(value=8.2, description='v0', max=11.0, min=6.0, step=0.01), Output()), _dom_classes=('widget-interact',))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from ipywidgets import interact\n",
    "\n",
    "@interact(v0=(6, 11, 0.01))\n",
    "def investigate_orbit(v0=8.2):\n",
    "    planetname = 'Mercury'\n",
    "    planet = data[planetname]\n",
    "    \n",
    "    fig, axes = plt.subplots(1, 2, figsize=(12, 4))\n",
    "    \n",
    "    x0 = (1+planet['e'])*planet['a']\n",
    "    x, y, vx, vy, t = orbit(v0=v0, tmax = 1, x0 = x0)\n",
    "    e = calc_eccentricity(x, y, t)\n",
    "\n",
    "    plt.sca(axes[0])\n",
    "    plot_orbit(x, y)\n",
    "\n",
    "    plt.sca(axes[1])\n",
    "    plt.plot(t, x)\n",
    "    plt.plot(t, y)\n",
    "    plt.xlabel('t (yr)')\n",
    "    plt.ylabel('x or y (AU)')\n",
    "    \n",
    "    plt.suptitle('{}: v0 = {:.2f}, e = {:.3f}, e_obs = {:.3f}'.format(planetname, v0, e, planet['e']), fontsize=14)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "def eccentricity_error(v0, planetname='Mercury'):\n",
    "    planet = data[planetname]\n",
    "    \n",
    "    x0 = (1+planet['e'])*planet['a']\n",
    "    x, y, vx, vy, t = orbit(v0=v0, tmax = 1, x0 = x0)\n",
    "    e = calc_eccentricity(x, y, t)\n",
    "    \n",
    "    return abs(e-planet['e'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.163734256538973"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy import optimize\n",
    "optimize.brent(eccentricity_error, brack=(7, 11))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why the difference?  Plot error vs initial velocity:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'error')"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEKCAYAAAA4t9PUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAHWtJREFUeJzt3X+UXOV93/H3x4sXTBuwQEuL9cMSQbEjnASbQSC7JYoVsPBJEYnBCGMjEnp0cMw5bXzsE3ES7EaEYtK0uMTUsWywgRNbYDnE2wYqE2Glqc9CtQIZJLDCIggsko1AMuACkiW+/eM+Sy+j2Z3Znbkzd2Y/r3Pm7My9z73znWHQ9z7Pc5/nUURgZmY2VW/pdABmZtbdnEjMzKwpTiRmZtYUJxIzM2uKE4mZmTXFicTMzJriRGJmZk1xIjEzs6Y4kZiZWVOO6HQA7TBz5syYN29ep8MwM+sqW7ZseT4iBuqVmxaJZN68eQwPD3c6DDOzriLpnxop56YtMzNrihOJmZk1xYnEzMya4kRiZmZNcSIxM7OmOJGYmVlTnEgmMDQ0xHXXXcfQ0FCnQzEzK61pMY5kKoaGhli6dCkHDhygv7+fjRs3snjx4k6HZWbWkKGhITZt2sSSJUsK/7fLiWQcmzZt4sCBAxw6dIgDBw6wadMmJxIz6wrtvhB209Y4lixZQn9/P319ffT397NkyZJOh2Rm1pBaF8JFco1kHIsXL2bjxo1tqxqambXK2IXwWI2k6AthRUShb1AGlUolPNeWmU0nregjkbQlIir1yrlGYmbWgxYvXty2lhT3kZiZWVOcSMzMrCmFJhJJyyTtkDQiaXWN/Z+W9KikhyVtlPTO3L6Vkh5Pj5W57adJeiSd80ZJKvIzmJnZxApLJJL6gJuAc4GFwMWSFlYVewioRMSvAuuBP0vHHgd8HjgDWAR8XtKMdMyXgVXAgvRYVtRnMDOz+oqskSwCRiJiZ0QcANYBy/MFIuL7EfFKenk/MDs9/xBwb0TsjYh9wL3AMkknAsdExFBkt5vdBpxf4GcwM7M6ikwks4Bncq9H07bxXA7cU+fYWel53XNKWiVpWNLwnj17Jhm6mZk1qshEUqvvouagFUkfByrAf6pzbMPnjIi1EVGJiMrAQN216826nicZtU4pchzJKDAn93o2sKu6kKTfBP4I+PWI2J87dknVsZvS9tlV2w87p9l040lGrZOKrJFsBhZImi+pH1gBDOYLSHov8BXgvIh4LrdrA3COpBmpk/0cYENE7AZelnRmulvrUuC7BX4Gs67Q7rmVzPIKq5FExEFJV5IlhT7glojYLmkNMBwRg2RNWf8c+Ha6i/fpiDgvIvZKuoYsGQGsiYi96fkngW8AbyPrU7kHs2mu3XMrmeV5rq0mtHO+f7N6/Hu0VvNcWwVzm7S1W71E0c65lczynEimyAtfWTv5wsXKzHNtTZEXvrJ2cme6lZlrJFPkha+sndyZbmXmznazLuHOdGs3d7ab9Rh3pltZuY/EzMya4kRiZmZNcSIxM+tSZZmo030kZmZdqExji1wjMTPrQmUaW+REYmbWhco0KNpNW2ZmXahMg6KdSMzMulRZxha5acvMzJpSaCKRtEzSDkkjklbX2H+WpAclHZR0QW77b0jamnu8Jun8tO8bkp7M7Tu1yM/QjLLcmmdmVqTCmrYk9QE3AWeTrbW+WdJgRDyaK/Y0cBnwmfyxEfF94NR0nuOAEeB7uSKfjYj1RcXeCmW6Nc8MPFeXFafIPpJFwEhE7ASQtA5YDryRSCLiqbTv9QnOcwFwT0S8Ulyoref1SqxMfGFjRSqyaWsW8Ezu9WjaNlkrgG9VbbtW0sOSbpB05FQDLFKZbs0zK9OYA+s9RdZIVGPbpOasl3Qi8CvAhtzmq4AfA/3AWuAPgTU1jl0FrAKYO3fuZN62Jcp0a56Z1zOxIhWZSEaBObnXs4FdkzzHR4G7IuLnYxsiYnd6ul/S16nqX8mVW0uWaKhUKh1ZdKUst+ZZdyiyD8MXNlakIhPJZmCBpPnAs2RNVB+b5DkuJquBvEHSiRGxW5KA84FtrQjWrJPa0YfhCxsrSmF9JBFxELiSrFnqMeDOiNguaY2k8wAknS5pFLgQ+Iqk7WPHS5pHVqP5+6pT/5WkR4BHgJnAnxb1GczaxX0Y1s0KHdkeEXcDd1dt+1zu+WayJq9axz5Fjc75iPhga6M06zz3YVg38xQpZiXgPgzrZk4kZiXhPgzrVp5ry8zMmuJEYmZmTXEiMTOzpjiRmJlZU5xIzMysKU4kZmbWFCeSDvGiV2bWKzyOpAO8NoSZNaJbFiNzIukAL3plZvV00wWnm7Y6wItemVk93TSRp2skHeB5lcysnm6ayFMRHVnzqa0qlUoMDw93Ogwzs0npdB+JpC0RUalXzjUSM+v4P1hWW7dM5FloH4mkZZJ2SBqRtLrG/rMkPSjpoKQLqvYdkrQ1PQZz2+dLekDS45LukNRf5Gcw63VjnbpXX301S5cu9S3pNmmFJRJJfcBNwLnAQuBiSQurij0NXAZ8s8YpXo2IU9PjvNz264EbImIBsA+4vOXBm00j3dSpa+VUZI1kETASETsj4gCwDlieLxART0XEw8DrjZwwrdP+QWB92nQr2brtZqVX1kGovovQmlVkH8ks4Jnc61HgjEkcf5SkYeAg8IWI+BvgeOCnaT34sXMethyvWdmUeUyA7yK0ZhWZSFRj22RuEZsbEbsknQTcJ+kR4KVGzylpFbAKYO7cuZN4W7PWK/sg1G7p1LVyKrJpaxSYk3s9G9jV6MERsSv93QlsAt4LPA+8XdJYAhz3nBGxNiIqEVEZGBiYfPRmLeTmI+tlRSaSzcCCdJdVP7ACGKxzDACSZkg6Mj2fCXwAeDSyQS/fB8bu8FoJfLflkZu12Fjz0TXXXFOqZi2zVih0QKKkDwNfBPqAWyLiWklrgOGIGJR0OnAXMAN4DfhxRJwi6f3AV8g64d8CfDEibk7nPIms4/444CHg4xGxf6I4PCDRzGzyGh2Q6JHtZmZWU6OJxJM2mplZU5xISqqsYw7MzKp5rq0SKvOYAzOzaq6RlJCnrDCzbuJEUkIec2Bm3cRNWyXkKSvMrJs4kZSUp6wws27hpi0zM2uKE4mZmTXFicTMzJriRGJm1iG9MvDYne1mZh3QSwOPXSMxs7p65cq5THpp4LFrJGY2oV66ci6TsYHHY99rNw88do3ErIV68cq9l66cy6SXFjtzjcSsRXr1yr2XrpzLplcGHhdaI5G0TNIOSSOSVtfYf5akByUdlHRBbvupkoYkbZf0sKSLcvu+IelJSVvT49QiP4NZo3r1yr2XrpytGIXVSCT1ATcBZwOjwGZJgxHxaK7Y08BlwGeqDn8FuDQiHpf0DmCLpA0R8dO0/7MRsb6o2M2mopev3HvlytmKUWTT1iJgJCJ2AkhaBywH3kgkEfFU2vd6/sCI+Mfc812SngMGgJ9iQNaM4kkdy8WTbdp0VWQimQU8k3s9Cpwx2ZNIWgT0A0/kNl8r6XPARmB1ROyvcdwqYBXA3LlzJ/u2pdarbfG9wFfuNh0V2UeiGttiUieQTgRuB343IsZqLVcB7wZOB44D/rDWsRGxNiIqEVEZGBiYzNuWXq+2xZtZdyoykYwCc3KvZwO7Gj1Y0jHA3wJ/HBH3j22PiN2R2Q98nawJbVrxwldmViZFNm1tBhZImg88C6wAPtbIgZL6gbuA2yLi21X7ToyI3ZIEnA9sa23Y5ee2eDMrE0VMqrVpcieXPgx8EegDbomIayWtAYYjYlDS6WQJYwbwGvDjiDhF0sfJahvbc6e7LCK2SrqPrONdwFbgioj42URxVCqVGB4ebvnnMzPrZZK2RESlbrkiE0lZOJGYmU1eo4mkbh+JpD5Jf9CasMzMrNfUTSQRcYhs/IeZmdlhGu1s/4GkLwF3AP93bGNEPFhIVGZm1jUaTSTvT3/X5LYF8MHWhmNmZt2moUQSEb9RdCBmZtadGhqQKOlYSf9F0nB6/GdJxxYdnJmZlV+jI9tvAV4GPpoeL5GN8zAzs2mu0T6SX4yIj+Re/4mkrUUEZGZm3aXRGsmrkv7V2AtJHwBeLSYkMzPrJo3WSK4Absv1i+wDVhYTkll5eR2Y2vy91DZdvpe6iUTSW4B3RcSvpRl5iYiXCo/MrGS8Dkxt/l5qm07fSyMj218HrkzPX3IS6Q5DQ0Ncd911DA0NdTqUnuF1YGrz91LbdPpeGm3aulfSZzh8ZPveQqKypkynK6F26uU12Zvh76W26fS9NJpIfi/9/VRuWwAntTYca4VaV0JOJM3zOjC1+XupbTp9L3WnkU99JIsj4gftCan1pts08q6RmFkrtGwa+dRH8udTDGKZpB2SRiStrrH/LEkPSjoo6YKqfSslPZ4eK3PbT5P0SDrnjWmlRMsZuxK65pprnETMrHCNNm19T9JHgL+OBlfCktQH3AScTbZ++2ZJgxHxaK7Y08BlwGeqjj0O+DxQIWtC25KO3Qd8GVgF3A/cDSwD7mnwc0wbixcvdgIxs7ZoNJF8GjgaOCTpNbJlbiMijpngmEXASETsBJC0jmxdkzcSSUQ8lfa9XnXsh4B7xzrzJd0LLJO0CTgmIobS9tvI1m13IjEz65BGR7YfS1Zz+NOUPE4hq2lMZBbwTO71aNrWiPGOnZWe1z2npFVjk0zu2bOnwbc1M7PJajSR3AScCVycXr8MfKnOMbX6LhpdIH68Yxs+Z0SsjYhKRFQGBgYafFszM5usRhPJGRHxKeA1gNRX0V/nmFFgTu71bGBXg+833rGj6flUzmlmZgVoNJH8PHWeB4CkAaC6X6PaZmCBpPmS+oEVwGCD77cBOEfSDEkzgHOADRGxG3hZ0pnpbq1Lge82eE4zMytAo4nkRuAu4ARJ1wL/G/iPEx0QEQfJplbZADwG3BkR2yWtkXQegKTTJY0CFwJfkbQ9HbsXuIYsGW0G1uRG0X8S+BowAjyBO9rNzDqq7oDENwpK7waWkvVTbIyIx4oMrJWm24BEM7NWaHRAYqO3/xIRPwJ+1FRUZmbWcxpt2jIzM6vJicTMzJriRGJmZk1xIjEzs6Y4kUxTXkGxNn8vZpPX8F1b1ju8Xklt/l7MpsY1kmloOq0lPRn+Xorjml5vc41kGppOa0lPhr+XYrim1/ucSKah6bSW9GT4eylGrZper3y3Q0ND/r3gRDJteQXF2vy9tF6v1vRc0/r/nEjMrFC9WtPr5ZrWZDmRmFnherGm16s1ralwIjEzm4JerWlNhROJmdkU9WJNayoKHUciaZmkHZJGJK2usf9ISXek/Q9Impe2XyJpa+7xuqRT075N6Zxj+04o8jOYmdnECkskaWnem4BzgYXAxZIWVhW7HNgXEScDNwDXA0TEX0XEqRFxKvAJ4KmI2Jo77pKx/RHxXFGfwczM6iuyRrIIGImInRFxAFgHLK8qsxy4NT1fDyxNa7HnXQx8q8A4zcysCUUmklnAM7nXo2lbzTJpjfcXgeOrylzE4Ynk66lZ6+oaicfMzNqoyERS6x/46gXiJywj6QzglYjYltt/SUT8CvCv0+MTNd9cWiVpWNLwnj17Jhe5mZk1rMhEMgrMyb2eDewar4ykI4Bjgb25/Suoqo1ExLPp78vAN8ma0A4TEWsjohIRlYGBgSY+hpmZTaTIRLIZWCBpvqR+sqQwWFVmEFiZnl8A3BcRASDpLcCFZH0rpG1HSJqZnr8V+C1gG2Zm1jGFjSOJiIOSrgQ2AH3ALRGxXdIaYDgiBoGbgdsljZDVRFbkTnEWMBoRO3PbjgQ2pCTSB/wd8NWiPsN05YnozGwylCoAPa1SqcTw8HCnw+gKnojOzMZI2hIRlXrlvLCVvYkXdzKzyXIisTcZm4iur6+vJyei80p9Zq3nubbsTXp5Ijo325kVw4nEDtOrE9F5/QizYrhpy6aNXm+2M+sU10hs2ujlZjuzTnIisWmlV5vtup3HLnU3JxIz6yjfBNH93EdiZh3lsUvdz4nEzDqq7DdBeOxRfW7aMrOOKvNNEG52a4wTiZl1XFlvgvDYo8a4acvMbBxlb3YrC9dIzMzGUeZmtzJxIjEzm0BZm93KpNCmLUnLJO2QNCJpdY39R0q6I+1/QNK8tH2epFclbU2Pv8wdc5qkR9IxN0qqte67mZm1SWGJRFIfcBNwLrAQuFjSwqpilwP7IuJk4Abg+ty+JyLi1PS4Irf9y8AqYEF6LCvqM5iZWX1F1kgWASMRsTMiDpCtvb68qsxy4Nb0fD2wdKIahqQTgWMiYiit7X4bcH7rQ7eJ+L56M8srso9kFvBM7vUocMZ4ZdIa7y8Cx6d98yU9BLwE/HFE/EMqP1p1zlkFxG7j8H31ZlatyBpJrZpF9QLx45XZDcyNiPcCnwa+KemYBs+ZnVhaJWlY0vCePXsmEbZNxNNZmFm1IhPJKDAn93o2sGu8MpKOAI4F9kbE/oh4ASAitgBPAL+Uys+uc07ScWsjohIRlYGBgRZ8HIPy31fvZjez9iuyaWszsEDSfOBZYAXwsaoyg8BKYAi4ALgvIkLSAFlCOSTpJLJO9Z0RsVfSy5LOBB4ALgX+osDPYFXKfF+9m93MOqOwRJL6PK4ENgB9wC0RsV3SGmA4IgaBm4HbJY0Ae8mSDcBZwBpJB4FDwBURsTft+yTwDeBtwD3pYW1U1vvqPZ2FWWcUOiAxIu4G7q7a9rnc89eAC2sc9x3gO+Occxh4T2sjtV4w1uw2ViMpW7ObWa/yyHbrGWVudjPrZU4k1lPK2uxm1ss8+6+ZmTXFicTMzJriRGJmZk1xIjGzUvMg0/JzZ7uZlZYHmXYH10jMrLQ8t1t3cCIxs9Iq+9xulnHTlpmVVjsGmQ4NDXkQa5OcSMys1IocZOo+mNZw05a1nO+ysW7hPpjWcI3EWspXeNZNPNFnaziRWEt5KnfrJp7oszWcSKylfIVn3cYTfTbPicRayld4ZtNPoYlE0jLgv5KtkPi1iPhC1f4jgduA04AXgIsi4ilJZwNfAPqBA8BnI+K+dMwm4ETg1XSacyLiuSI/h01O0XfZOEmZlUthiURSH3ATcDYwCmyWNBgRj+aKXQ7si4iTJa0ArgcuAp4H/k1E7JL0HrLlemfljrskrZRo04g78s3KqcjbfxcBIxGxMyIOAOuA5VVllgO3pufrgaWSFBEPRcSutH07cFSqvdg05ls1zcqpyEQyC3gm93qUN9cq3lQmIg4CLwLHV5X5CPBQROzPbfu6pK2SrpakWm8uaZWkYUnDe/bsaeZzWEl4ugyzciqyj6TWP/AxmTKSTiFr7jont/+SiHhW0i8A3wE+QdbP8uaTRKwF1gJUKpXq97Uu5I58s3IqMpGMAnNyr2cDu8YpMyrpCOBYYC+ApNnAXcClEfHE2AER8Wz6+7Kkb5I1oR2WSKw3+VZNs/IpsmlrM7BA0nxJ/cAKYLCqzCCwMj2/ALgvIkLS24G/Ba6KiB+MFZZ0hKSZ6flbgd8CthX4GczMrI7CEknq87iS7I6rx4A7I2K7pDWSzkvFbgaOlzQCfBpYnbZfCZwMXJ36QrZKOgE4Etgg6WFgK/As8NWiPoOZmdWniN7vPqhUKjE87LuFzcwmQ9KWiKjUK+fZf83MrClOJGZm1hQnEjMza4oTiZl1NS+k1nme/dfMupbnXysH10is7XwFaa3i+dfKwTUSaytfQVoreSG1cnAisbbyUrzWSp5/rRycSKytfAVprVZv/jUvhlY8JxJrK19BWju5KbU9nEis7XwFae3iptT2cCKxUvEVpLWSm1Lbw4nESsVXkNZKbkptDycSKxVfQVqreTG04jmRWKn4CtKs+xQ6sl3SMkk7JI1IWl1j/5GS7kj7H5A0L7fvqrR9h6QPNXpO636LFy/mqquuchIx6xKFJRJJfcBNwLnAQuBiSQuril0O7IuIk4EbgOvTsQvJluY9BVgG/DdJfQ2e08zM2qjIGskiYCQidkbEAWAdsLyqzHLg1vR8PbBUktL2dRGxPyKeBEbS+Ro5p5mZtVGRiWQW8Ezu9WjaVrNMWuP9ReD4CY5t5JxmZtZGRSYS1dhWvUD8eGUmu/3wN5dWSRqWNLxnz54JAzUzs6krMpGMAnNyr2cDu8YrI+kI4Fhg7wTHNnJOACJibURUIqIyMDDQxMcwM7OJFJlINgMLJM2X1E/WeT5YVWYQWJmeXwDcFxGRtq9Id3XNBxYA/6fBc5qZWRsVNo4kIg5KuhLYAPQBt0TEdklrgOGIGARuBm6XNEJWE1mRjt0u6U7gUeAg8KmIOARQ65z1YtmyZcvzkv5pih9lJvD8FI8tmmObGsc2NY5taro5tnc2chJlFQAbj6ThiKh0Oo5aHNvUOLapcWxTMx1i81K7ZmbWFCcSMzNrihNJfWs7HcAEHNvUOLapcWxT0/OxuY/EzMya4hqJmZk1xYkEkPQuSVtzj5ck/fuqMpJ0Y5p1+GFJ7ytRbEskvZgr87l2xJbe+w8kbZe0TdK3JB1VtX/cGZ5LENtlkvbkvrd/28bY/l2Ka3v1f8+0vyO/twZja9vvTdItkp6TtC237ThJ90p6PP2dMc6xK1OZxyWtrFWmg7Edyn1/LR8LN05sF6b/pq9LGvdOLU1lhvWI8CP3IBuf8mPgnVXbPwzcQzZNy5nAAyWKbQnwPzoQzyzgSeBt6fWdwGVVZX4f+Mv0fAVwR4liuwz4Uge+t/cA24CjycZy/R2woAy/twZja9vvDTgLeB+wLbftz4DV6flq4Poaxx0H7Ex/Z6TnM8oQW9r3sw58b78MvAvYBFTGOa4PeAI4CegHfggsrPd+rpEcbinwRERUD2BcDtwWmfuBt0s6sSSxddIRwNvSFDdHc/iUNePN8FyG2Drll4H7I+KVyCYr/Xvgt6vKdOr31khsbRMR/4tssHJe/jd1K3B+jUM/BNwbEXsjYh9wL9mSFGWIrXC1YouIxyJiR51DpzTDuhPJ4VYA36qxvQwzD48XG8BiST+UdI+kU9oRTEQ8C/w58DSwG3gxIr5XVWy8GZ7LEBvAR1LT0XpJc2rsL8I24CxJx0s6mqz2Uf3enfq9NRIbdOD3lvMvImI3QPp7Qo0ynfr+GokN4Chlk8reL6kjyWYcU/renEhylM3fdR7w7Vq7a2xr2y1vdWJ7kKy569eAvwD+pk0xzSC7WpkPvAP4Z5I+Xl2sxqGFf28NxvbfgXkR8atkTTi30gYR8RjZIm73Av+TrPngYFWxjnxvDcbWkd/bJHX0/9cGzI1sRPnHgC9K+sVOB5RM6XtzInmzc4EHI+InNfY1PPNwQcaNLSJeioifped3A2+VNLMNMf0m8GRE7ImInwN/Dby/qsx4Mzx3PLaIeCEi9qeXXwVOa0NcY+99c0S8LyLOIvs+Hq8q0rHfW73YOvh7G/OTsWa+9Pe5GmU69f01EhsRsSv93UnWZ/HeNsTWiCl9b04kb3Yx4zcdDQKXprtpziRrKtndvtDGj03Svxzrd5C0iOy/6wttiOlp4ExJR6f3Xwo8VlVmvBmeOx5bVZ/DedX7iyTphPR3LvA7HP7ftmO/t3qxdfD3Nib/m1oJfLdGmQ3AOZJmpNrpOWlbx2NLMR2Zns8EPkA2QW0ZTG2G9SLvHOimB1ln7AvAsbltVwBXpOciWy/+CeARxrnroUOxXQlsJ2uGuB94fxtj+xPgR2Rt67cDRwJrgPPS/qPImuNGyJYCOKlEsV2X+96+D7y7jbH9A9k/Hj8Elpbs91Yvtrb93siS2G7g52RXy5eT9bFtJKspbQSOS2UrwNdyx/5e+t2NAL9bltjIasaPpO/vEeDyNsX22+n5fuAnwIZU9h3A3bljPwz8Y/rt/VEj7+eR7WZm1hQ3bZmZWVOcSMzMrClOJGZm1hQnEjMza4oTiZmZNcWJxKyDip6h1qwdfPuvWYdIOg4YJhtjEMAW4LTIJhk06xqukZi1gaTrJf1+7vV/AD5JwTPUmrWDE4lZe6wDLsq9/ijZ2g+dnlHarGlHdDoAs+kgIh6SdIKkdwADwD7gVbJpW95UtO3BmTXJNRKz9llPNnHlRWQ1lE7PKG3WEu5sN2uTtADUV4GZwK+TTZ63hWxJVMjW+TgtItoxzb5Zy7hpy6xNImK7pF8Ano00Jbyka8im7gZY4yRi3cg1EjMza4r7SMzMrClOJGZm1hQnEjMza4oTiZmZNcWJxMzMmuJEYmZmTXEiMTOzpjiRmJlZU/4fjQjhJsu8tBYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x12076bdd8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for v0 in np.arange(7, 11, 0.1):\n",
    "    err = eccentricity_error(v0)\n",
    "    plt.plot(v0, err, 'k.')\n",
    "plt.xlabel('v0')\n",
    "plt.ylabel('error')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For comparison, the book shows a completely analytical way of determining the initial velocity of Mercury. In any case, we can now look at the precession of Mercury's orbit over longer times if we include the effects of general relativity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1d963b4fbf15463daaa28582aa3d6642",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>interactive</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='α', max=0.02, step=0.0001), Output()), _dom_classes=('widget-interact',))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@interact(α=(0, 0.02, 0.0001))\n",
    "def investigate_precession(α = 0):\n",
    "    planet = data['Mercury']\n",
    "    x0 = (1+planet['e'])*planet['a']\n",
    "    v0 = 8.2 # AU/yr\n",
    "    \n",
    "    x, y, vx, vy, t = orbit(v0=v0, tmax = 1, x0=x0, α=α)\n",
    "\n",
    "    fig, axes = plt.subplots(1, 2, figsize=(12, 4))\n",
    "    \n",
    "    plt.sca(axes[0])\n",
    "    plot_orbit(x, y)\n",
    "\n",
    "    plt.sca(axes[1])\n",
    "    r = np.sqrt(x**2 + y**2) \n",
    "    plt.plot(t, r)\n",
    "    plt.xlabel('t (yr)')\n",
    "    plt.ylabel('r (AU)')\n",
    "    \n",
    "    plt.suptitle('Simulation of the precession of Mercury\\nα = {:.4f}'.format(α))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate the derivative of the distance from Mercury to the Sun and look for where the derivative goes from positive to negative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4c2471665a70454490151dffcb745608",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>interactive</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='α', max=0.02, step=0.001), Output()), _dom_classes=('widget-interact',))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@interact(α=(0, 0.02, 0.001))\n",
    "def investigate_precession(α = 0):\n",
    "    planet = data['Mercury']\n",
    "    x0 = (1+planet['e'])*planet['a']\n",
    "    v0 = 8.2 # AU/yr\n",
    "    \n",
    "    x, y, vx, vy, t = orbit(v0=v0, tmax = 1, x0=x0, α=α)\n",
    "\n",
    "    fig, axes = plt.subplots(1, 2, figsize=(12, 4))\n",
    "    \n",
    "    plt.sca(axes[0])\n",
    "    plot_orbit(x, y)\n",
    "\n",
    "    plt.sca(axes[1])\n",
    "    r = np.sqrt(x**2 + y**2)\n",
    "    \n",
    "    drdt = np.zeros(len(t)-2)\n",
    "    for i in range(1,len(t)-2):\n",
    "        drdt[i] = (r[i+1] - r[i-1]) / (t[i+1] - t[i-1])\n",
    "    plt.plot(t[1:-1], drdt)\n",
    "    plt.xlabel('t (yr)')\n",
    "    plt.ylabel('dr/dt (AU/yr)')\n",
    "    \n",
    "    for i in range(1, len(t)-2):\n",
    "        if drdt[i-1] < 0 and drdt[i] > 0:\n",
    "            # at the perihelion\n",
    "            plt.sca(axes[0])\n",
    "            plt.plot([0, x[i]], [0, y[i]], \"g-\" )\n",
    "            plt.sca(axes[1])\n",
    "            plt.axvline(t[i], color='green')\n",
    "            \n",
    "    plt.suptitle('Simulation of the precession of Mercury\\nα = {:.4f}'.format(α))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimate the rate of change of position of the perihelion / apihelion,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "af6dac094c714fd8b22b976b857e95f4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>interactive</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.001, description='α', max=0.01, step=0.0001), Output()), _dom_classes=('widget-interact',))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@interact(α=(0, 0.01, 0.0001))\n",
    "def investigate_precession(α = 0.001):\n",
    "    planet = data['Mercury']\n",
    "    x0 = (1+planet['e'])*planet['a']\n",
    "    v0 = 8.2 # AU/yr\n",
    "    \n",
    "    x, y, vx, vy, t = orbit(v0=v0, tmax = 2, x0=x0, α=α)\n",
    "\n",
    "    fig, axes = plt.subplots(1, 2, figsize=(12, 4))\n",
    "    \n",
    "    plt.sca(axes[0])\n",
    "    plot_orbit(x, y)\n",
    "\n",
    "    r = np.sqrt(x**2 + y**2) \n",
    "    drdt = np.zeros(len(t)-2)\n",
    "    for i in range(1,len(t)-2):\n",
    "        drdt[i] = (r[i+1] - r[i-1]) / (t[i+1] - t[i-1])\n",
    "\n",
    "    for i in range(1, len(t)-2):\n",
    "        if drdt[i-1] < 0 and drdt[i] > 0:\n",
    "            # at perihelion\n",
    "            plt.sca(axes[0])\n",
    "            plt.plot([0, x[i]], [0, y[i]], \"g-\" )\n",
    "            \n",
    "            plt.sca(axes[1])\n",
    "            # calculate degrees\n",
    "            θ = np.degrees(np.arctan2(y[i], x[i]))\n",
    "            plt.plot(t[i], θ, 'go')\n",
    "            \n",
    "    plt.sca(axes[1])\n",
    "    plt.xlabel('t (yr)')\n",
    "    plt.ylabel(r'$\\theta$ (deg)')\n",
    "            \n",
    "    plt.suptitle('Simulation of the precession of Mercury\\nα = {:.4f}'.format(α))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fit a best fit curve to estimate the rate of precession"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "fa15f3376eea4bf689f6f3073aaa2ebe",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>interactive</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.001, description='α', max=0.005, step=0.0001), Output()), _dom_classes=('widget-interact',))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@interact(α=(0, 0.005, 0.0001))\n",
    "def investigate_precession(α = 0.001):\n",
    "    planet = data['Mercury']\n",
    "    x0 = (1+planet['e'])*planet['a']\n",
    "    v0 = 8.2 # AU/yr\n",
    "    \n",
    "    x, y, vx, vy, t = orbit(v0=v0, tmax = 2, x0=x0, α=α)\n",
    "\n",
    "    fig, axes = plt.subplots(1, 2, figsize=(12, 4))\n",
    "    \n",
    "    plt.sca(axes[0])\n",
    "    plot_orbit(x, y)\n",
    "    plt.title(r'$\\alpha$ = {:.4}'.format(α))\n",
    "\n",
    "    r = np.sqrt(x**2 + y**2) \n",
    "    drdt = np.zeros(len(t)-2)\n",
    "    for i in range(1,len(t)-2):\n",
    "        drdt[i] = (r[i+1] - r[i-1]) / (t[i+1] - t[i-1])\n",
    "\n",
    "    t_vals = []\n",
    "    θ_vals = []\n",
    "    for i in range(1, len(t)-2):\n",
    "        if drdt[i-1] < 0 and drdt[i] > 0:\n",
    "            # at perihelion\n",
    "            plt.sca(axes[0])\n",
    "            plt.plot([0, x[i]], [0, y[i]], \"g-\" )\n",
    "            \n",
    "            plt.sca(axes[1])\n",
    "            # calculate degrees\n",
    "            θ = np.degrees(np.arctan2(y[i], x[i]))\n",
    "            plt.plot(t[i], θ, 'go')\n",
    "            \n",
    "            t_vals.append(t[i])\n",
    "            θ_vals.append(θ)\n",
    "    \n",
    "    p = np.polyfit(t_vals, θ_vals, 1)\n",
    "    plt.sca(axes[1])\n",
    "    plt.plot(t_vals, np.polyval(p, t_vals), 'r-')\n",
    "    plt.xlabel('t (yr)')\n",
    "    plt.ylabel(r'$\\theta$ (deg)')\n",
    "\n",
    "    plt.xlim(0, max(t))\n",
    "    plt.title('slope = {:.4}'.format(p[0]))\n",
    "            \n",
    "    plt.suptitle('Simulation of the precession of Mercury')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now that we have the algorithm worked out, let's rewrite the code without any of the graphics:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.1962447195946893"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def compute_precession(α = 0):\n",
    "    planet = data['Mercury']\n",
    "    x0 = (1+planet['e'])*planet['a']\n",
    "    v0 = 8.2 # AU/yr\n",
    "    \n",
    "    x, y, vx, vy, t = orbit(v0=v0, tmax = 2, x0=x0, α=α)\n",
    "\n",
    "    r = np.sqrt(x**2 + y**2) \n",
    "    drdt = np.zeros(len(t)-2)\n",
    "    for i in range(1,len(t)-2):\n",
    "        drdt[i] = (r[i+1] - r[i-1]) / (t[i+1] - t[i-1])\n",
    "\n",
    "    t_vals = []\n",
    "    θ_vals = []\n",
    "    for i in range(1, len(t)-3):\n",
    "        if drdt[i-1] < 0 and drdt[i] > 0:\n",
    "\n",
    "            # calculate degrees\n",
    "            θ = np.degrees(np.arctan2(y[i], x[i]))\n",
    "           \n",
    "            t_vals.append(t[i])\n",
    "            θ_vals.append(θ)\n",
    "    \n",
    "    p = np.polyfit(t_vals, θ_vals, 1)\n",
    "\n",
    "    return p[0]\n",
    "\n",
    "compute_precession(0.0001)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that from general relativity,\n",
    "\n",
    "$$\\alpha = 1.1 \\times 10^{-8}\\;\\mathrm{AU}^2$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is too small a value of $\\alpha$ to use directly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.077577246715755402"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_precession(1.1e-8) ## GIVES AN INCORRECT VALUE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explore how the rate of precession depends on $\\alpha$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgIAAAGLCAYAAAC1EvpBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xd4lGX2//H3CYISQEBlEcEkuK66otICoggsYttVUbHCYAOJvbC2r8bCLsbFFXvdoALqiL1gF5EiWEMgKOhPFElEXAXBRhAIOb8/niFLSSBDJnkmyed1XXMluZ92BsE5ucu5zd0RERGR+ikl7ABEREQkPEoERERE6jElAiIiIvWYEgEREZF6TImAiIhIPaZEQEREpB5TIiAiIlKPKREQERGpx5QIiITIzE40MzezPTZo283MxpvZMjP71cyeMrMWm1zXxsxKzaxnvM8wsx3NbISZ/Tnx76jc5x9jZi+a2RIzW2Nm35vZ82bWryaev5XY9or92fw17FhEwqJEQCRcmcAKd18IYGbtgY+AHYEIcD5wJHDvJtcdDywF3o/3GbGfbwQaVjn6LTCz7czsMWAisBq4DDgc+D+gFfCWmTWpzhgqoUvs68ehRiESou3CDkCknssE8gHMzIAJwBxggMfqf5vZXsDVZnaOu/8eu+54YKK7l8bzjJjOBB/M8xPzFip0PzAQONXdn9mgfRow1swGufvKRD7QzBYB49x9RCUv6QoscvdliYxDpDZRj4BIDYj9dnyDmS0ys1Vm9qaZ7U7wQbT+Q/oE4EDg777xJiBFQCNgt9i9dgT+ArwY7zPM7DNgNLA9sDbWLf5sNbzffsAw4JZNkoAy7v7EBuc3NLPrzGxhLPb3zWxvM8s0s2Iz2y3RMcZ0BT42sywzm29mK2PP7rLVK0XqCPUIiFSz2G/6TwD9gH8CnwBHA5OBlsCs2KlDCbr6F5rZhv82m8a+lsS+Hg2sAd7ehmecATwJzANujrV9V068DSrx1tzd11Vw7BqgGPj31m5iZg2B14EDYtd9T9CbcBOwA3Cfuy+pRDzbonMszsbAFUAT4BbgZTPbK9E9FiLJSImASPUbCpwEHOzuH8Ta3jGzg4E/AbPMrBHQl+ADaW0591gLrP8wPB54w91Xx/OMWFsB0A64Z4PzNtUHmFKJ9zWNoGdiI2bWMvZeJrj7z5W4zwXAocAh7v5e7B7dgEsBA84q76ItJCwpmyRS5SYsZvZHoAXwIdB/g6GYXwkSk57AW5WIX6RWUyIgUv3+D3ipnA/ez4B9ga+AjgRJwIUEkwU39CTBZL+SWMJwFMGHZ7zPAOhAMMyQT8VmAd229qaAXytoP4Bg2PGTStwD4DzgrfVJQMxPQHNghLv/WMF1FSUs18de65WbsBAMCwBkbzIUsz7uVluJW6ROUCIgUo1iqwD+CNxQzuG2QL67u5llxNpmuPvcDa5vDbQnmEQIQdd/Y+DVeJ8R+7kL4AQTEivy21aOr+cVtDePff1+azcws12BfYAHNjnUCFgO3L6Fy8tLWCYCrwC5G7RVlLB0Ab5y91mbtLeJfV28hWeL1BlKBESqV9vY1/9u2Bj7AOwN3BdrWv9vcdMu7NMJPnDHxX4+Hpjm7j9twzMgGBP/yt1/2ULMVRoa4H8JQLtK3GPP2Nev1zeYWQOCuQxfuXtFH+LEjuVt2GZma4Al7p5X/lUb6cr/hls2dCrwM5VbmilS6ykREKle65eldQDe2aB9BMFEuPVd9Is2OG8elH2QXw3kuvtXsTHxY4GcbXwGBMMEW1s2WNWhgXyCCYhnmtmtm8xlwMxSgQ7u/jGwfvnjThuccn4szi0NXyRCF+B3M9vO3UtisbUjGKq4zd3XVPPzRZKCEgGR6rWA4IP3htgktMXAaQTr6+F/k/hmEYzn/8vMfidY3jcS+JJgNjtAD4Ju65e28RkQjL13MbMjCX7rXbDpGHx5v2nHw93XmtkFwLPAR2Z2J7CQYMigJ8HkvysIivgUxOK4wcx+JuhF+DfwHHBsrOLflA3qJyREbChmJ4KeiHFmNjb27BGxmG6u6FqRukZ1BESqUWy2+okEH9T3A48R/Bb8AMFY/ILYeU5QR+Bb4GngDoI6AYe5e3HsdscDee6+eFueEXMDQdf9iwRd39VSZtjdXwQOIejpuBWYRDBE0RUYRTCWT2x53snAKoL3fQ1B/YHzCD6QX6FySxnjtX6i4DEEcy5ejsX1PHCkegOkPrGNJ8uKSLIys8+Bx9x906EBEZFtpkRARESkHtPQgIiISD2mREBERKQeUyIgIiJSjykREBERqcfqRB2BXXbZxTMyMsIOQ0REpEbMmjVrmbsnZD+MOpEIZGRkkJe3zfVPREREahUzK0zUvTQ0ICIiUo8pERAREanHlAiIiIjUY0oERERE6jElAiIiIvWYEgEREZF6TImAiIhIPaZEQEREpB5TIiAiIlKPKREQERGpJtFolIyMDFJSUsjIyCAajYYd0mbqRIlhERGRZBONRsnKyqK4uBiAwsJCsrKyAIhEImGGthH1CIiIiFSD7OzssiRgveLiYrKzs0OKqHxKBERERKpBUVFRXO1hUSIgIiJSDdLS0uJqD4sSARERkWqQk5NDamrqRm2pqank5OSEFFH5lAiIiIhUg8jJJ5Pbty/pgAHpaWnk5uYm1URB0KoBERGRxPvySxg4kEheHpFzz4Xbb4dNegeShRIBERGRRHrsMbjgAmjYEJ57DgYMCDuiLdLQgIiISCL88gsMHgxnnAFdukBBQdInAaBEQEREpOo+/BA6d4YJE+Cf/4R33oHddw87qkpRIiAbeeWVVzj++OPZbbfdaNSoEa1bt2bAgAFMnjy52p+9ePFiLr74Yg466CBSU1MxMxYtWpTw66v6nET55ptvOOmkk2jevDk77rgjAwYM2Or64qOOOgoz47rrrquhKEVki0pLYdQoOOQQKCmB6dPh+uuhQYOwI6s0JQICQElJCaeffjr9+/dn++23584772TSpEmMGjWKpUuXcsQRR7By5cpqjeHLL7/k6aefpmXLlvTq1avarq/qcxKhuLiYQw89lM8//5zx48fz2GOPsWDBAvr27Vvhn/OECRMoKCio4UhFpEJLlsDhh8M11wRDAAUF0LNn2FHFz91r/atr164uVTNs2DBv0KCBP/300+Uej0aj1R7DunXryr4fM2aMA/71118n/PqqPqcy0tPT/cYbb6zw+J133ukpKSm+YMGCsraFCxd6gwYN/Lbbbtvs/BUrVnjr1q39iSeecMCzs7MTGq+IxGniRPedd3ZPTXV/+GH30tIafTyQ5wn6DFWPgDB58mTGjBnD1Vdfzcknn1zuOYMGDar2OFJSqvbXsbLXx/OcgoIC+vfvT8uWLWncuDE9e/bk3Xff3dYQy0ycOJEePXqw5557lrW1b9+enj178tJLL212/lVXXUWHDh0YOHBglZ8tIlXw++9w8cXQv38wB2DWLBgyBMzCjmybKREQ/vWvf5GamspVV121Tde7OyUlJVt9rVu3LsGRV6/8/HwOPvhgli9fzpgxY3juuefYeeedOeyww5g1a1aV7j1v3jz222+/zdo7dOjA/PnzN2qbMWMGjz76KPfff3+VnikiVTR/PnTvDvfeC8OHwwcfwD77hB1VlamOQD23YsUKpkyZwsCBA2nevPk23WPatGn07dt3q+f16dOHqVOnbtMzwnDllVeSlpbGO++8Q6NGjQA48sgj2W+//Rg5ciQvvvgiECRC5SU5paWllJSUlP1sZjSITSBavnw5LVu23OyanXbaiRUrVpT9vHbtWs4991yuuOIK9t5774S+PxGpJHf4z3+CD/9mzeC11+Cvfw07qoRRIlDPzZ07l9LSUvbff/9tvkfXrl35+OOPt3pes2bNtvkZNW3VqlVMmzaNa6+9lpSUlI0+0A877DCi0WjZzxUlQiNHjmTkyJFlP2+aCFk5XYnB0N//3HLLLaxatSrpti0VqTeWL4dzzoEXXggmBj76KOy6a9hRJZQSgXru559/BqB169bbfI+mTZvSqVOnrZ5X3gdfslq+fDnr1q3b7MN8Q6WlpaSkpJSbCPXv359jjjmGrKyssrYNE6GWLVuyfPnyze65YsWKsp6CoqIicnJyeOihh1i9ejWrV68uO2/16tX89NNPNGvWrKyXQUQSbNq0oEDQ99/DrbfC3/8OVZzLlIxCSwTMbAdgOrB9LI5n3f1GM2sPPAnsBOQDp7v7mrDirOvWJwCLFy/e5nvUxaGBFi1akJKSwoUXXsgZZ5xR7jnrJx02a9aMzMzMjY41atSI3XbbbbP29Tp06MC8efM2a58/fz777rsvAAsXLuT3339n8ODBm503evRoRo8ezezZsyuVhIlIHEpKgqJAOTmwxx7w/vvQtWvYUVWbMHsEVgOHuvtvZtYQmGFmrwN/B+5w9yfN7EFgKPBAiHHWaV26dKFNmzaMHz+eK6+8ku23336j48XFxcybN49u3bpVeI+6ODTQpEkTevXqRUFBAV26dKnyioZN9e/fnyuuuIKFCxeyxx57ALBo0SJmzpzJqFGjAOjUqRNTpkzZ7Nq+ffsyePBghg4dutGqAxHZsmg0SnZ2NkVFRaSlpZGTk7P5ToCLFkEkAu+9B2eeCffcE8wLqMNCSwRi6yB/i/3YMPZy4FBg/Vq18cAIlAhUm4YNG3L//fdz0kkn0b17dy677DL22GMPfv75Z2bOnMm4ceMYPXr0FhOB8n4j3lbPPvssQNms/Ndff51WrVrRqlUr+vTpAwQ9EP369eORRx7Z7Lf1ylxf2fNuv/12evfuzZFHHsnQoUNp06YNy5YtIz8/n3Xr1pV9YG+LYcOGce+993Lcccdx0003YWZcf/317L777px77rlA0Cvxl7/8pdzr09PTKzwmIpuLRqNkZWVRXFwMQGFhYdnQXVky8NRTcO65weTAJ56A+rJcN1EFCbblBTQA5hAkBLcAuwBfbnB8d+DTCq7NAvKAvLS0tCoXZ6jv3n//fe/fv7/vvPPO3rBhQ2/Xrp3369fPb7/9dv/pp59qLA6CZHCzV58+fcrOmTJligM+duzYbbo+nvPmz5/vp556qrdq1cobNWrkbdu29WOPPdZfffXVLb6PrRUUcncvLCz0AQMGeLNmzbxp06Z+3HHHVaqwESooJBK39PT0cv/Np6enu//2m/uQIe7g3qOH+8KFYYe7VSSwoJD5JrOUw2BmLYAXgBuAse6+Z6x9d+A1d9/ilPbMzEzPy8ur/kBFRKRWSklJ2WxVDgSTmEv/9CdYsCAoFTxiRLB9cJIzs1nunpCu2KRYNeDuP5nZVKAH0MLMtnP3EqAdsCTU4EREpNZLS0ujsLBw83Z3+O03mDwZKjHpuS4KbR2EmbWK9QRgZo2Bw4DPgCnASbHTzgQ2r7cqIiISh5ycHFJTUzdqSwVyunQJNguqp0kAhFtiuA0wxczmAh8Dk9z9FeBq4O9m9iWwM/BwiDGKiEgdEIlEyM3NJf0Pf8CAdCD3zDOJ5OXBLruEHV6owlw1MBfoXE77QqB7zUckIiJ11urVRPLzifzwA3ToABMmQBUqqtYlSTFHQEREpNp88UWwFDA/H84/H267DRo3DjuqpKFEQERE6iZ3GD8eLroItt8+2C/g+OPDjirp1L2iySIiIj//DIMGwdlnQ2ZmMCFQSUC5lAiIiEjd8v770KkTPPMM3HRTsDSwXbuwo0paSgRERKRuWLcu2CioV6/g53ffhexs0A6dW6Q5AiIiUvstXgynnw5Tp8Jpp8GDD0Lz5mFHVSsoERARkdrtxRdh6FBYvRrGjg12DTQLO6paQ0MDIiJSO61aBRdcACecABkZwfLAs85SEhAnJQIiIlL7fPopdOsGDzwAl18eTBDca6+wo6qVlAiIiEjt4Q733RcsCVy2DN54A0aPhkaNwo6s1tIcARERqR2WLQvmAkycCEcdBePGQevWYUdV66lHQEREkt+UKdCxI7z+Otx+O7z6qpKABFEiICIiyWvt2qAWQL9+0KwZfPghDB8OKfr4ShQNDYiISHJauDAoE/zhh8GQwF13QZMmYUdV5ygREBGR5DNhApx7bvCb/1NPwSmnhB1RnaW+FRERSR6//hrUAhg0CPbfH+bMURJQzZQIiIhIcsjLgy5d4LHH4IYbYNq0oFCQVCslAiIiEq7S0qAWwMEHw++/BysE/vEP2E6j1zVBf8oiIhKe//432BvgrbdgwAAYMwZ22insqOoV9QiIiEg4Xn8dDjgg2C74wQfh2WeVBIRAiYCIiNSs1auJ/vWvZPztb6QsXUpGixZEmzbVZkEhUSIgAHzxxReYGa+//nrYoWxk6tSpmNlmrxYtWsR9r6OOOgoz47rrrtvs2MyZMzniiCP4wx/+wI477kiXLl145JFHEvEWNvLNN99w0kkn0bx5c3bccUcGDBhAUVHRRucsXryYiy++mIMOOojU1FTMjEWLFiU8FpFQfP450b32IuuNNygEHCj87juysrKIRqNhR1cvKREQAPLz8wHo1q1byJGU7+677+b9998ve7399ttxXT9hwgQKCgrKPTZ37lwOO+ww1q5dy5gxY3juuefo1q0bQ4cO5YEHHkhE+AAUFxdz6KGH8vnnnzN+/Hgee+wxFixYQN++fVm5cmXZeV9++SVPP/00LVu2pFevXgl7vkio3OHhh6FrV7IXL6Z4k8PFxcVkZ2eHElq95+61/tW1a1eXqrniiis8IyMj7DA2M2XKFAd80qRJ23yPFStWeOvWrf2JJ55wwLOzszc6fs0113jDhg39119/3aj9wAMP9B49emzzczd15513ekpKii9YsKCsbeHChd6gQQO/7bbbytrWrVtX9v2YMWMc8K+//jphcYjUuBUr3E8+2R3c+/VzM3OCzoCNXmYWdqS1BpDnCfoMVY+AADBr1iy6detGbm4u++67L02aNOGggw4q6ymoza666io6dOjAwIEDyz2+Zs0aGjZsSOPGjTdqb9GiBaWlpZudX1BQQP/+/WnZsiWNGzemZ8+evPvuu1uNY+LEifTo0YM999yzrK19+/b07NmTl156qawtRTXUpS6ZOTPYLOiFF2DUKHjrLdLS0so9taJ2qV76P44AMHv2bGbOnMnLL7/M6NGjGTduHN9//z3HHnvsRt3Wm3J3SkpKtvpat25dleKLRCI0aNCAnXfemUGDBm02rl6RGTNm8Oijj3L//fdXeM5ZZ50FwCWXXMKSJUv46aefGDNmDJMnT2b48OEbnZufn8/BBx/M8uXLy4YRdt55Zw477DBmzZq1xVjmzZvHfvvtt1l7hw4dmD9/fqXej0itsW4d/POf0Lt3UA9g5ky4+mpISSEnJ4fU1NSNTk9NTSUnJyekYOu5RHUthPnS0EDVfPnllw74kUce6aWlpWXtr7/+ugP+5ptvVnjt+q77rb369OmzTbHl5+f75Zdf7hMnTvSpU6f6HXfc4a1atfLddtvNv//++y1eu2bNGt933303GgqgnKEBd/ePPvrI27ZtWxZvw4YN/aGHHtrsvEMPPdT32WcfX716dVlbSUmJ77PPPn7cccdtMZ6GDRv61VdfvVl7dna2N2jQoNxrNDQgtVJhoXuvXsFQwODB7j//vNkpjz/+uKenp7uZeXp6uj/++OMhBFp7kcChARUUkrLfZHNycrANlu/sv//+ACxdurTCa7t27crHH3+81Wc0a9aswmPuvlmPwXaximKdO3emc+fOZe19+vShd+/edO/enbvvvpubbrqpwvvecsstrFq1aqsTkBYsWMCJJ55Ihw4dePDBB2ncuDEvvfQS5513HjvssAORSASAVatWMW3aNK699lpSUlIoKSkpu8dhhx1GNBrd4nsBNvrz3fD9i9QZzz0H55wDJSXw6KNw+unlnhaJRMr+bUm4lAgI+fn5/PGPf6Rr164btX/33XcAtGvXrsJrmzZtSqdOnbb6jPI+ANebNm0affv23ahtSx+OXbp0Ya+99tpiAlJUVEROTg4PPfQQq1evZvXq1WXHVq9ezU8//USzZs1o0KAB1157LQ0bNuSVV16hYcOGAPTr148ff/yRSy+9lIEDB5KSksLy5ctZt24dI0eOZOTIkeU+d8qUKfTr16/c99KyZUuWL1++2TUrVqygZcuWFb4XkVqhuBiGD4fcXOjWDZ54AjaYDyPJS4mAMGvWLHbbbbfN2p966imaN2/OQQcdVOG15X2Il6dPnz5MnTq13GOV7VXYkLtvMblYuHAhv//+O4MHD97s2OjRoxk9ejSzZ8+mU6dOfPLJJ3Ts2LEsCVive/fuPPHEE/zwww/suuuutGjRgpSUFC688ELOOOOMcp+79957V/heOnTowLx58zZrnz9/Pvvuu++W3q5IcisogIED4bPP4KqrYORIaNQo7KikkpQICPn5+eywww6UlJSUdWMvXryYBx98kMsvv5xGW/gHnYihgWbNmpGZmVnpePPy8vjiiy84ZQtbk3bq1IkpU6Zs1t63b18GDx7M0KFDy2bv77rrrsyZM4c1a9Zs9F4//PBDdthhB3aKlTxt0qQJvXr1oqCggC5dulQ4u7+i99K/f3+uuOIKFi5cyB577AHAokWLmDlzJqNGjarcmxdJJu5w771w5ZXQsiVMmgSHHRZ2VBInqwvjk5mZmZ6Xlxd2GLXSokWLaN++Pe3bt+fggw/m7LPPZvHixYwYMYK2bdvyzjvvbDERqG6RSIT27dvTpUsXWrRowezZs/nXv/5Famoq+fn57LLLLkDQM9GvXz8eeeSRCn9bh2CIIjs7e6O5Bc8++ywnn3wyRxxxBBdccAGNGzdm4sSJ3HfffQwfPpzbb7+97Nz8/Hx69+7NQQcdxNChQ2nTpg3Lli0jPz+fdevWbfEDfeXKlXTs2JHGjRtz0003YWZcf/31/Prrr8ydO5emTZtuFBPA5MmTefDBB7n//vtp1aoVrVq1ok+fPtv85ymSMEuXwpAh8MorcPTRMHYstGoVdlT1hpnNcvfK/wa1JYmadRjmS6sGtt2zzz7rgM+bN88HDBjgjRs39l133dX//ve/+2+//RZ2eH7zzTf7/vvv7zvuuKNvt9123q5dOx82bJgvWbJko/PWr14YO3bsFu9HBasGXnvtNe/Tp4/vsssu3rRpU+/YsaPfd999XlJSstm58+fP91NPPdVbtWrljRo18rZt2/qxxx7rr7766lbfT2FhoQ8YMMCbNWvmTZs29eOOO67cFQEkePWFSEK9/bZ7mzbujRq533WX+warjaRmkMBVA+oREBGRylm7Fm64AW65BfbeG558MigWJDUukT0CmiMgIiJb99VXMGgQfPQRZGXBHXfAJkWBpHZSIiAiIlsWjcL550ODBvDMM3DSSWFHJAmkEsMiIvVcNBolIyODlJQUMjIy/rcd8K+/whlnwODB0KlTsExQSUCdox4BEZF6LBqNkpWVRXFxsDFwYWEhWVlZsHAhkfHj4euvYcQIyM4O9gyQOie0HgEz293MppjZZ2Y2z8wujbWPMLNvzWxO7PW3sGIUEanrsrOzy5KA9YqLi8m+4YZgcuC0aXDjjUoC6rAw/8uWAJe7e76ZNQNmmdmk2LE73H10iLGJiNQLFe3kWQQwZ05QKEjqtNB6BNz9O3fPj33/K/AZ0DaseERE6qO0tLSK25UE1AtJMVnQzDKAzsCHsaaLzGyumT1iZvqbKCJSTXJGjCB1k27/1NRUcm6+OaSIpKaFngiYWVPgOeAyd/8FeAD4I9AJ+A64rYLrsswsz8zytrRNroiIVOCzz4jccQe5JSWkN2uGmZGenk5ubq62CK5HQq0saGYNgVeAN9399nKOZwCvuPt+W7qPKguKiMTBHR56CC69FJo0gXHjgv0CpNZIZGXBMFcNGPAw8NmGSYCZtdngtBOAT2s6NhGROmvFCjj55KA6YM+eMHeukoB6LsxVAz2B04FPzGxOrO1aYKCZdSLYZGURcG444YmI1DHvvguRCHz3Hfz733D55VDBdtpSf4SWCLj7DMDKOfRaTcciIlKnlZTATTfByJHQvj289x506xZ2VJIkVCFCRKQuKywMegFmzgzKBd97LzRrFnZUkkSUCIiI1FXPPAPDhkFpabBx0KBBYUckSUiDQyIidc3KlXDOOXDKKbD33kGFQCUBUgElAiIidcns2dC1KzzyCFxzDcyYAXvsEXZUksSUCIiI1AXucOed0KMH/PILvP023HwzNGwYdmSS5DRHQESktvvhBzj7bHjtNTj22KA3YJddwo5Kagn1CIiI1GZvvQUHHACTJwcrAl56SUmAxEWJgIhIbbRmDVx5JRx5JOy8M3z0EVx4IVh55VlEKqahARGR2mbBAhg4EGbNgvPOg9tug9TUsKOSWkqJgIhIbeEOjz4a/ObfqBE8/zyccELYUUktp6EBEZHa4OefgwqBZ50VLA8sKFASIAmhREBEJIlEo1EyMjJISUkhIyODaDQKH3wAnTvD008H+wW88w7svnvYoUodoaEBEZEkEY1GycrKori4GIDCwkKyzj4bSkqIpKXB9Olw8MEhRyl1jXoERESSRHZ2dlkSsF7x2rVkN24clAlWEiDVQImAiEiSKCoqKr991Spo0aKGo5H6QomAiEiSSGvXrvz2tLQajkTqEyUCIiLJ4NNPyXFn02oAqamp5OTkhBKS1A9KBEREwuQODzwA3boRWbOG3CuvJD09HTMjPT2d3NxcIpFI2FFKHaZVAyIiYfnxRzjnHHjxxaBU8PjxRFq3JvLvf4cdmdQj6hEQEQnD1KnQsSO8+mpQIvi116B167CjknpIiYCISE1auxauuw4OPRSaNAmKBf3975Ci/x1LODQ0ICJSU77+GgYNCj78zz4b7r4bmjYNOyqp55QIiIjUhAkTgp0C139/2mnhxiMSo74oEZHq9NtvwW//gwZBhw5BhUAlAZJElAiIiFSX/Hzo0gXGjw/mBUyfDu3bhx2VyEaUCIiIJFppabASoEcPWLUKpkwJdg3cTqOxknz0t1JEJJH++18480x46y044QR46CHYaaewoxKpkHoEREQS5fXXg9oA06cH1QKfe05JgCQ9JQIiIlW1enVQC+BvfwuKAuXlBSsEzMKOTGSrNDR6LUGzAAAgAElEQVQgIlIV/+//wcCBMHs2XHQR/Pvf0Lhx2FGJVFrciYCZbQ/sBjQGlrr70oRHJSKS7Nxh7Fi4+OLgg/+ll6B//7CjEolbpYYGzKyZmZ1vZtOBn4EvgU+B/5rZN2Y2xsy6VWegIiJJ46efgloAQ4fCgQdCQYGSAKm1tpoImNlwYBEwBJgEHAd0AvYCDgJuJOhZmGRmb5jZn6otWhGRGhSNRsnIyCAlJYWMjAyi0Si89x506hRMBLz5Zpg0Cdq2DTtUkW1WmaGBg4E+7v5pBcc/Ah4xs/OAoUAfYEGC4hMRCUU0GiUrK4vi4mIACgsLyTr7bCgpIZKRATNnBr0BIrWcuXvYMVRZZmam5+XlhR2GiNQhGRkZFBYWbtaenprKou++gx13DCEqkYCZzXL3zETcq9LLB82soZl9aGZ7J+LBIiLJrKioqPz2VauUBEidUulEwN3XAu2B2t+FICKyFWnt2pXfnpZWw5GIVK94CwqNB4ZVRyAiIklj7lxy1q0jdZPm1NRUcnJyQglJpLrEmwg0AbLMbI6ZPWxmd2/4iudGZra7mU0xs8/MbJ6ZXRpr38nMJpnZgtjXlnHGKCKybdzh3nuhe3cipaXkXn016enpmBnp6enk5uYSiUTCjlIkoeKaLGhmU7Zw2N390Dju1QZo4+75ZtYMmAUcD5wFLHf3UWb2f0BLd796S/fSZEERqbJly2DIEHj55aBU8Nix8Ic/hB2VSLkSOVkwrsqC7t43EQ+N3es74LvY97+a2WdAW4I6BX+JnTYemApsMREQEamSd96BwYPhxx/hzjvhkku0T4DUG3ENDZjZi2Z2jJkldLMiM8sAOgMfAq1jScL6ZEEpuYhUj7Vr4Zpr4LDDoHlz+PBDuPRSJQFSr8T7gb4SeApYbGY3J6KKoJk1BZ4DLnP3X+K4LsvM8swsb+lSbXcgInH66is45BAYNSooFZyXF1QMFKln4koE3D0CtAFGAocB/8/MppvZGWYW93ZbZtaQIAmIuvvzsebvY/MH1s8j+KGCWHLdPdPdM1u1ahXvo0WkPotGoXNn+OILeOYZGDMGmjQJOyqRUMTdxe/uv7j7A+7eHdifYJLffwg2IPqPmf25MvcxMwMeBj5z99s3ODQRODP2/ZnAS/HGKCJSrl9/hTPPDOYDdOwYbBZ00klhRyUSqm0e6zez3Qgm9h0DlADPArsDc83sikrcoidwOnBobDniHDP7GzAKONzMFgCHx34WEamajz+GLl3g8cfhxhthyhRQcSCR+FYNxLryjyPYifBwYDbwb2CCu/8WO+cUIBcYvaV7ufsMoKIZOf3iiUtEpEKlpXDbbXDttdCmDUydCr16hR2VSNKIKxEgWO5nwBPA/7n73HLOmQSsqGpgIiJV9t13wVDApElw4onBXICWqlEmsqF4E4HhwDPu/ntFJ7j7CoI9CUREwvPqq3DWWbByJeTmwjnnaFmgSDkqNUfAzG4zs14Es/srTAJEREK3ejVcdhkccwy0bQuzZsGwYUoCRCpQ2cmCqcAEgqV948zs+G1ZLigiUq0++wwOPBDuuiuoDvjBB/DnSi1kEqm3KpUIuPv57t4OOBr4FrgJWGZmE81siJlpIb+IhMc9GP/v2hW+/TbYL+Cuu2CHHcKOTCTpxVtQ6CN3z3b3/YCOwDSCTYIWm9kMM7vCzNpWQ5wiIuVbsQJOOQWysqBnT5g7NxgWEJFK2eY6Au7+pbvf5u69gXbAI8AhwMBEBSciskUzZgRlgV98EW65Bd58M1giKCKVFu+qgXK5+1KCROCRRNxPRGSLSkogJwf++U9o3x5mzoTu3cOOSqRWiregUEUf9A78DiwAnnb3JVUNTESkXEVFEIkEvQGnnw733QfNmoUdlUitFW+PQCugF1AKfBpr24+gyNAsYAAw0sx6ufuchEUpIgLw3HNBPYCSEnjssWDPABGpknjnCMwEXgfauXvvDeYHvAa8BaQDrwK3JTRKEanfVq4keuihZJx0Eik//UTGjjsSVV0AkYSINxG4FPinuxevb4h9nwMMd/c1wC2ANvUWkcQoKCD6pz+RNWUKhQTjkIVLlpCVlUU0Gg07OpFaL95EoClQ3pTcXWPHAH4hQZMQRaQec4e774bu3cn+/nuKNzlcXFxMdnZ2KKGJ1CXxJgIvAA+b2clmlmFm6WZ2MvAw8HzsnO7AF4kMUkTqmaVL4dhj4dJL4YgjKHIv97SioqIaDkyk7ok3ETgPeBN4HFgYez0OvAFcEDvnM2BYogIUkXpm0iQ44AB4+2245x6YOJG0tLRyT62oXUQqL97KgsXufh6wE8E8gC7ATrESxCtj58zRigERiduaNXDVVXDEEcFWwR99BBddBGbk5OSQmpq60empqank5OSEFKxI3RF3ZUEz+yvwFPAksNzdV5rZOWbWL+HRiUj98OWXQXngW2+Fc8+FvLygVyAmEomQm5tLeno6ZkZ6ejq5ublEIpEQgxapG+ItKBQBHgQeAvoBDWOHGgBXAZMTGp2I1H2PPQYXXAANGwZ1AgYMKPe0SCSiD36RahBvj8BVwDB3Hw6UbND+AVoyKCLx+OWXoCDQGWdAly5QUFBhEiAi1SfeROBPwPvltP8G7Fj1cESkXvjwQ+jcGZ58Mtgv4J13YPfdw45KpF6KNxFYAuxVTntv4KuqhyMidVppKYwaBYccAuvWwfTpcP310KBB2JGJ1FvxJgK5wN1m1jP28+5mdibwb+CBhEYmInXLkiVw+OFwzTXBEMCcOXDwwWFHJVLvxTVZ0N3/bWbNgUnADsAUYDUw2t3vq4b4RKQuePllOPtsWLUKHn44+F57BYgkhbiXD7p7NrALQQXBHkArd78+0YGJSB2walVQC6B//2AOQH4+DBmiJEAkiWzTngCxjYbyEhyLiNQl8+bBaafBp5/C8OHwr3/B9tuHHZWIbGKriYCZPVLZm7n7kKqFIyK1njv85z/Bh/+OO8Lrr8NRR4UdlYhUoDI9Aq02+bk3UAp8Evt5P4IhhukJjEtEaqPly+Gcc+CFF+DII2H8eGjdOuyoRGQLtpoIuPux6783s2uAVcDZ6/cWMLMmBLsPflL+HUSkXpg2LSgQ9P33MHp00COQEvc0JBGpYfH+K70EGLE+CQCIfT8SuDiRgYlILVFSEtQC6NsXGjeG99+Hyy9XEiBSS8T7L7UpsFs57W2A1HLaRaQuW7QIeveGm26Cs84KVgV07Rp2VCISh3hXDTwHjDWzKwn2F4BgCeEtwPOJDExEktxTTwU7BbrDhAnBCgERqXXi7RE4H3gZGEdQUvgrYDzwKnBBQiMTkeT0229E+/Qh47TTSPn5ZzKaNiW6bl3YUYnINoq3suAq4IJYj8AfAQO+3HDOgIjUYfn5RI8+mqz//pfiWFPhkiVkZWUBaJtgkVpom2bzuPtKd5/r7gVKAkTqgdJSuP126NGD7KVLy5KA9YqLi8nOzg4lNBGpmq0mAmbWvrI3s4D2EhWpS77/Ho4+OlgJcPTRFJWWlntaUVFRDQcmIolQmR6B983sYTM7qKITzKylmZ0PzAeOS1h0IhKuN9+EAw6AqVPh/vvh+edJS0sr99SK2kUkuVUmEdgHWA68amZLzewNMxtrZg+Y2ZNmNhf4ARgMXObu91ZnwCJSA1avDnoAjjoK/vAH+PhjOP98MCMnJ4fU1I1XC6emppKTkxNSsCJSFVtNBNz9J3e/EmhLsGrgc6AF0B4oIVg10Nnde7r7m9UZrIjUgC++gIMPDuYEXHghfPQR7Ldf2eFIJEJubi7p6emYGenp6eTm5mqioEgtZe4ezoODzYyOAX5w9/1ibSOAYcDS2GnXuvtrW7tXZmam5+VpM0SRKnGHcePg4ouDXQIfeQSO00ifSDIys1nunpmIe4VZA3QcUN6WZHe4e6fYa6tJgIgkwE8/wcCBMGQIdO8Oc+cqCRCpJ0JLBNx9OsHcAxEJ03vvQadO8OyzcPPNMGkStG0bdlQiUkOScVeQi8xsrpk9YmYtww5GpM5aty7YI6B372CDoBkz4JproEGDsCMTkRoUVyJgZmlmZuW0m5klYu3QAwQVCzsB3wG3bSGWLDPLM7O8pUuXVnSaiJRn8WLo1y/YNfCUU2D2bOjRI+yoRCQE8fYIfA20Kqd9p9ixKnH37919nbuXAmOA7ls4N9fdM909s1Wr8kISkXK9+CJ07Ah5ecHkwGgUmjcPOyoRCUm8iYAB5S0zaAr8XtVgzKzNBj+eAHxa1XuKSMyqVUEtgBNOgPbtg16AM8+EzTv5RKQeqdSmQ2Z2d+xbB/5lZhuWGm9A8Jv7nHgebGYTgL8Au5jZYuBG4C9m1in2nEXAufHcU0Qq8MknwaqAefPgiisgJwcaNQo7KhFJApXdfXD/2FcD/gys2eDYGiAfGB3Pg919YDnND8dzDxHZCvegNPDll0OLFkHJ4COOCDsqEUkilUoE3L0vgJmNBS5191+qNSoRqbply2DoUJg4Ef72Nxg7NigXLCKygcr2CADg7mdXVyAikkBTpsDgwUEycOedcMklmgsgIuXaaiIQKwVcKe4+pGrhiEiVrF0LN94Io0bBXnvBK69A585hRyUiSawyPQKbrs3rDZQCn8R+3o9g9cH0BMYlIvFauBAGDYIPP4Rzzgl6Apo0CTsqEUlyW00E3P3Y9d+b2TXAKuBsd18Za2tCMMnvk/LvICLVbsIEOPfcoELg00/DySeHHZGI1BLx1hG4BBixPgkAiH0/Erg4kYGJSCX8+mtQC2DQINh/fygoUBIgInGJNxFoCuxWTnsbILXq4YhIpeXlQZcu8PjjcMMNMG0apKeHHZWI1DLxJgLPAWPN7DQzy4i9TiMYGng+8eGJyGZKS+HWW+Ggg2D1aqLZ2WSMH09Ko0ZkZGQQjUbDjlBEapG4lg8C5xNsBDQOaBhrKyFIBK5IXFgiUq7vvguGAiZNggEDiB55JFnDh1NcHBT7LCwsJCsrC4BIJBJmpCJSS5h7eVsHbOWiYILgHwkqDX654ZyBMGRmZnpeXl6YIYhUv9deg7POgt9+C1YEDBtGRvv2FBYWbnZqeno6ixYtqvEQRaRmmNksd89MxL0qNTRgZreZWS8zS4FggqC7z3X3grCTAJE6b/VquOwyOPpo2HXXYG5AVhaYUVRUVO4lFbWLiGyqsnMEUoEJwPdmNs7MjjezxtUYl4gAfP459OgBd90FF18MH30E++5bdjgtLa3cyypqFxHZVKUSAXc/393bAUcD3wI3AcvMbKKZDTGzTYsOiUhVuMNDD0HXrrB4Mbz8Mtx9N+yww0an5eTkkJq68YKd1NRUcnJyajJaEanF4lo14O4fuXu2u+8HdASmAWcBi81shpldYWZtqyFOkfpjxQo45RQYNixYGVBQAMccU+6pkUiE3Nxc0tPTMTPS09PJzc3VREERqbRtmiy42U3MdgH6x14z3D2uLYmrSpMFpc6YOTMoDrRkCdx0E1x5ZVAtUERkA4mcLBjv8sFyufsy4JHYS0TiVVICOTnwz39CRkaQEHTvHnZUIlIPbPVXDTPb2cweNrP/mlmJmf1oZu+Z2a1mpv9TiVRVURH07QsjRgS9AbNnKwkQkRpTmR6Bx4E9gBxgKcEKgvtj115qZu8BQ9x9YbVFKVJXPfdcsFNgSQk8+iicfnrYEYlIPVOZwcc+wInufo+7P+nujwBrgdMI9h2YBbxnZntVY5widUtxcVAL4KST4E9/CnoBlASISAgqkwh8C+xc3gF3X+bulwOjgbsSGZhInVVQAJmZMGYMXH01zJgBe+4ZdlQiUk9VJhG4i2CjoS0NWj4D9EpMSCJ1lDvcc08w/r9iRbBfwKhR0KhR2JGJSD221TkC7n6vmbUB3jezqcALbJ5AnE4wf0BEyrN0KQwZAq+8EpQKHjsWWqkOl4iEr1LLB90928xeINhhcBTQGJhnZsuAZgQ7EZ5RbVGK1GZvvw1nnAE//hhUB7zoIjALOyoRESCOOgLungecZmbbAZ2AvYHmwDLgnVgtARFZb80auP56uPVW2HtveP116Ngx7KhERDay1UTAzCpTJOgoM8PdhyQgJpHa76uvYOBA+PjjYHXAHXfAJnsCiIgkg8r0CGw6kNkbKAU+if28H8GcgekJjEuk9nr8cTj/fNhuO3j2WTjxxLAjEhGp0FZXDbj7setfwHvAm0A7d+/t7r2B3YE3gA+rN1SRJPfLL0R79iTj9NNJ+e03Mpo0Ifr772FHJSKyRfHuZnIJMMLdV65viH0/Erg4kYGJ1CoffUR0zz3Jeu89CgEHCr/9lqysLKLRaNjRiYhUKN5EoClBNcFNtSEoPSxSv5SWBrUAevYke/lyijc5XFxcTHZ2diihiYhURryJwHMExYVOM7OM2Os04GHg+cSHJ5LEliyBI46Aa66BE06gqLS03NOKiopqODARkcqLNxE4H3gZGAd8FXuNB14FLkhoZCLJ7OWX4YAD4P33g1LBTz1FWlpauadW1C4ikgziSgTcfZW7X0Cw90BnoAuwk7tf4O6b9oqK1D2//w6XXAL9+8Puu8OsWcHugWbk5OSQuskSwdTUVHJyckIKVkRk6+LtEQCCCYLuPtfdCzacOChSp82fDwceGOwXcNll8MEHsM8+ZYcjkQi5ubmkp6djZqSnp5Obm0skEgkxaBGRLTN3DzuGKsvMzPS8vLyww5C6yh1yc2H4cGjaFMaNg7/9LeyoRKQeM7NZ7p6ZiHttU4+ASL2xfDmcdBKcdx4cckiwhbCSABGpQ5QIiFRk+vRgb4CXXw72C3jjDWjTJuyoREQSSomAyKZKSuCGG6BvX9hhB3jvPbjiCkjRPxcRqXsqvfugSL2waBFEIsGH/5lnBhMDmzULOyoRkWqjREBkvaefDnYKLC2FaBQGDQo7IhGRahdaX6eZPWJmP5jZpxu07WRmk8xsQexry7Dik3pk5UoYOhROPTVYDjhnjpIAEak3whz0HAcctUnb/wGT3f1PwOTYzyLVZ/Zs6NoVxo6Fa6+Fd9+FPfYIOyoRkRoTWiLg7tOB5Zs0H0dQspjY1+NrNCipP0pL4Y47oEcP+PVXmDwZcnKgYcOwIxMRqVHJNg26tbt/BxD7+oeKTjSzLDPLM7O8pUuX1liAUgf88AMccwz8/e9w1FFBbYC+fcOOSkQkFMmWCFSau+e6e6a7Z7Zq1SrscKS2eOutYLOgd96B++6DF1+EXXYJOyoRkdAkWyLwvZm1AYh9/SHkeKSuWLMmqAVw5JHBB//HH8MFF4BZ2JGJiIQq2RKBicCZse/PBF4KMRapK774Ag46CG67Dc4/P0gC9t8/7KhERJJCmMsHJwDvA3ub2WIzGwqMAg43swXA4bGfRbaNe7BBUJcuQaGgF16A+++Hxo3DjkxEJGmEVlDI3QdWcKhfjQYiddPPPwcbBT35JPTpA48/Du3ahR2ViEjSSbahAZGq++AD6NQJnnkGbropWBqoJEBEpFxKBKTuWLcuqAVwyCHBz+++C9nZ0KBBuHGJiCQx7TUgdcO338LgwTB1Kpx2Gjz4IDRvHnZUIiJJTz0CUvu99FJQG+Djj4NSwU88Ac2bE41GycjIICUlhYyMDKLRaNiRiogkHfUISO21ahVcfjk88ECwMmDCBNhrLwCi0ShZWVkUFxcDUFhYSFZWFgCRSCS0kEVEko16BKR2+vRT6NYtSAIuvxzef78sCQDIzs4uSwLWKy4uJjs7u6YjFRFJakoEpHZxD2oBdOsGy5bBG2/A6NHQqNFGpxUVFZV7eUXtIiL1lRIBqT1+/BFOOAEuvBD+8pdgs6Ajjyz31LS0tLjaRUTqKyUCUjtMnQodO8Jrr8Htt8Orr0Lr1hWenpOTQ2pq6kZtqamp5OTkVHOgIiK1ixIBSW5r1wa1AA49FJo0gQ8/hOHDIWXLf3UjkQi5ubmkp6djZqSnp5Obm6uJgiIimzB3DzuGKsvMzPS8vLyww5BE+/prGDQoqBQ4ZAjcdRc0bRp2VCIioTOzWe6emYh7afmgJKcJE4K9AszgqafglFPCjkhEpE7S0IAkl99+g7PPDnoC9tsP5sxREiAiUo2UCEjymDUrKAz06KNwww0wbRpkZIQdlYhInaZEQMJXWhrUAjjooKBa4JQp8I9/wHYauRIRqW5KBKRGVFj3/7//hb/+Fa68Eo49NqgN0Lt3uMGKiNQj+pVLql2Fdf8LCoiMHw+//BLsFpiVFUwOFBGRGqPlg1LtMjIyKCws3Kw9HVi0//7BCoEOHWo+MBGRWkrLB6VWqbDuP8BHH8EOO9RoPCIi8j+aIyDVrsK6/+npSgJEREKmRECqXU52NqkNGmzUprr/IiLJQYmAVK+ZM4nk5JBbWkp6ixaq+y8ikmQ0R0Cqx7p1cPPNMGIEZGQQ+eADIt27hx2ViIhsQomAJN4338DgwTB9OkQicP/9sOOOYUclIiLlUCIgifX883DOOcH2wY8+CqefHnZEIiKyBZojIIlRXBzsFnjiibDnnjB7tpIAEZFaQImAVN3cudCtG/znP3DVVTBjRpAMiIhI0lMiINvOHe65B7p3h+XL4a234JZboFGjsCMTEZFK0hwB2TZLl8KQIfDKK3D00TB2LLRqFXZUIiISJ/UISPwmT4aOHYMegLvugpdfVhIgIlJLKRGQylu7Fq65Bg4/HJo3hw8/hEsu0Y6BIiK1mIYGpHK++goGDQo2CRo2DO64A5o0CTsqERGpIiUCsnXRKJx/PjRoAM88AyedFHZEIiKSIBoakIr9+iuccUZQJbBjRygoUBIgIlLHKBGoh6LRKBkZGaSkpJCRkUE0Gt38pI8/hs6dg96AESNgyhSoYDthERGpvTQ0UM9Eo1GysrIoLi4GoLCwkKysLIBgN8DSUrj1VrjuOmjTBqZOhV69QoxYRESqk7l72DFUWWZmpufl5YUdRq2QkZFBYWHhZu3p6eksev/9YCjg7beDIYDcXGjZMoQoRURkS8xslrtnJuJe6hGoZ4qKiipuP+AAWLkSxoyBoUO1LFBEpB5IykTAzBYBvwLrgJJEZT0CaWlp5fYIpLlD27YwYQL8+c8hRCYiImFI5smCfd29k5KAxMrJySE1NXWjtlQg58gj4YMPlASIiNQzyZwISDWIRCLk5uaSvtNOGJCekkLu5ZcTeeMN2GGHsMMTEZEalpRDA4ADb5mZA/9x99xNTzCzLCALgu5uqaQVK4i88AKR5cuhXz947LFgdYCIiNRLydoj0NPduwB/BS40s96bnuDuue6e6e6ZrbThTeW8+25QGOill4Ltgt96S0mAiEg9l5SJgLsviX39AXgB6B5uRLVcSUlQFOgvf4FGjeC99+CqqyAlKf/zi4hIDUq6TwIza2JmzdZ/DxwBfBpuVLVYYWGQAPzjH0Gp4NmzoVu3sKMSEZEkkYxzBFoDL1iwhn074Al3fyPckGqpZ54JdgosLYXHH4dIJOyIREQkySRdIuDuC4GOYcdRq61cCZdeCg8/DN27wxNPwB//GHZUIiKShJJuaECqaM4c6NoVHnkErrkGZsxQEiAiIhVSIlBXuMOdd8KBB8IvvwT7Bdx8MzRsGHZkIiKSxJJuaEC2wQ8/wNlnw2uvwbHHBr0Bu+wSdlQiIlILqEegtps0KdgsaPJkuOeeoEaAkgAREakkJQK11Zo1QS2AI46AnXeGjz6Ciy7SjoEiIhIXDQ3URgsWwKBBkJcH550Ht90Gm2wkJCIiUhlKBGoTd3j0UbjwwqBC4PPPwwknhB2ViIjUYhoaqC1++SWoDHjWWcHywIICJQEiIlJlSgRqgw8+gE6d4KmnYORIeOcd2H33sKMSEZE6QIlAMlu3LqgFcMghQZng6dPhuuugQYOwIxMRkTpCcwSS1bffwumnw5QpcOqp8OCD0KJF2FGJiEgdo0QgGU2cGBQI+v33oDjQWWdpWaCIiFQLDQ2ELBqNkpGRQUpKChlpaUSPOAKOOw7S0yE/P0gIlASIiEg1UY9AiKLRKFlZWRQXFwNQ+M03ZH3zDfz1r0ReeAG23z7kCEVEpK5Tj0CIsrOzy5KA9YqB7PnzlQSIiEiNUCIQoqKiorjaRUREEk2JQFimTiUtpfw//rS0tBoORkRE6islAjWtpASuvx4OPZScVq1I3WGHjQ6npqaSk5MTUnAiIlLfKBGoSYsWQe/ecNNNcNZZRBYsIPehh0hPT8fMSE9PJzc3l0gkEnakIiJST5i7hx1DlWVmZnpeXl7YYWzZk0/CuecG3//nP3DaaeHGIyIitZaZzXL3zETcSz0C1e2332DIEBg4EDp0gDlzlASIiEjSUCJQnfLzg50Cx40L9giYPh3atw87KhERkTJKBKpDaSncfjv06AErVwa7BY4cCdupfpOIiCQXfTIl2vffw5lnwptvwvHHw0MPwc47hx2ViIhIudQjkEhvvgn/v737D7W7ruM4/nw7dW74Y23NGnM2J8OaoWupKNTADJsjXEKYcf8MhlOhQqGNUWgwKSXaH4aywFl2yZmJKAgxlrlMmhvtRxOZXvcDl0MnsTKMme7TH9+Pde6Pc87u3Tnn+73n+3zA4Z7zPZ/zOZ/vi8/dee9zvvf7vewyeP55ePBBePJJiwBJUqVZCHTC8eNw552wbBmcfz5s3w633urFgiRJledXA6fq1VeLvwLYuRNuvx3uvx+mTSt7VJIknRRXBBoMuyTw/PkMDg42b5wSbNwIS5bAoUPw1FPwwAMWAZKkScUVgWzUJYEPHWLlypUAo8/0d+xYsfS/aRNcey08+ijMndvrIUuSdMpcEcjGvCTwe++xdu3a4Q1ffBEWL4YnnoB774XNmy0CJEmTloVA1vaSwB9+WFwjYOlSOO00eOEFWLMGpkzp4SglSeosC4Gs2aV/L7zwQjh8GK67rrhq4M03FwcGXn11j0coSVLnWQhk69atY/r06cO2TZ8+nXU33QSXXw47dhSnCh4chPPOK2eQkiR1mIVANjY/7uQAAAcaSURBVDAwwIYNG/5/SeB589hwzTUMrF9fXB9g587ijIGeG0CS1EcsBBoMDAxw8OBBTuzezcFzz2Vgyxa4667iAMGFC8seniRJHeefD460cSOsWgUzZhSnDL7++rJHJElS17giMNLMmcWBgXv2WARIkvqeKwIjrVgBN97osQCSpFqo5IpARCyLiH0RMRQRq0sYQM/fUpKkMlSuEIiIKcDPgBuARcA3I2JRuaOSJKk/Va4QAK4ChlJK+1NK7wOPAStKHpMkSX2pioXAXOCNhseH87ZhImJlROyIiB1Hjx7t2eAkSeonVSwExvqCPo3akNKGlNIVKaUrZs+e3YNhSZLUf6pYCBwG5jU8vgB4s6SxSJLU16pYCGwHFkbERRFxJnAL8HTJY5IkqS9V7jwCKaUPIuIO4HfAFODhlNLLJQ9LkqS+VLlCACCl9CzwbNnjkCSp31XxqwFJktQjFgKSJNWYhYAkSTVmISBJUo1FSqPO1TPpRMRR4NAYT30ceKfHw5lMzKc9M2rNfNozo9bMp72xMvpUSqkjZ9Pri0KgmYjYkVK6ouxxVJX5tGdGrZlPe2bUmvm01+2M/GpAkqQasxCQJKnG+r0Q2FD2ACrOfNozo9bMpz0zas182utqRn19jIAkSWqt31cEJElSC5UuBCJiWUTsi4ihiFg9xvNTI2JTfn5bRMxveG5N3r4vIr7Srs98tcNtEfFa7vPMbu9fJ/Q4o0ci4kBE7Mq3xd3ev1PVpXwejoi3I2LviL5mRsTmPIc2R8THurlvndDjfO6OiL81zJ/l3dy3Tul0RhExLyKei4hXIuLliPh2Q/vaz6E2+TiHim1nRcRLEbE7Z3RPQ/uLYryfZSmlSt4orjz4OrAAOBPYDSwa0eY24KF8/xZgU76/KLefClyU+5nSqk/gceCWfP8hYFXZGVQwo0eAr5e932Xmk59bCiwB9o7o6z5gdb6/Gvhx2RlULJ+7gbvK3u+yMwLmAEtym3OAVxt+x2o/h9rk4xwq+gvg7NzmDGAbcHV+PO7PsiqvCFwFDKWU9qeU3gceA1aMaLMC+EW+/wRwXURE3v5YSul4SukAMJT7G7PP/Jov5T7IfX6ti/vWKT3LqAf70g3dyIeU0lbg72O8X2Nfk2EO9TqfyajjGaWUjqSU/gKQUnoXeAWYO0ZftZxDbfKZjLqRUUop/Su3PyPf0kQ/y6pcCMwF3mh4fJjRk+F/bVJKHwD/AGa1eG2z7bOAY7mPZu9VRb3M6CPrImJPRPw0IqZ2Yie6qBv5tPKJlNKR3NcR4PwJj7w3ep0PwB15/jw8GZa96XJGeQn4cxT/owPnULt8wDk0FyAipkTELuBtYHNKaRsT/CyrciEQY2wb+ScOzdp0anvV9TIjgDXAp4ErgZnA905umKXpRj79pNf5PAhcDCwGjgA/aTfACuhaRhFxNvBb4DsppX9OeITl6nU+zqH82pTShymlxcAFwFUR8dmTfK9RqlwIHAbmNTy+AHizWZuIOB04j2JJstlrm21/B5iR+2j2XlXUy4zIS3YppXQc2EheCq6wbuTTylsRMSf3NYeiUq+ynuaTUnor/+N1Avg51Z8/0KWMIuIMig+5wZTSkw1tnEM0z8c5NPq1KaVjwB+AZUz0s6zsAyma3YDTgf0UB0h8dIDFpSPa3M7wAywez/cvZfgBFvspDrBo2ifwG4YfYHFb2RlUMKM5+WcA64EflZ1Br/NpeN18Rh8Mdz/DD/S6r+wMKpbPnIb736X47rP0HHqdUf79+SWwfoz3q/0capOPc6jIaDYwI7eZBvwR+Gp+PO7PstJDahPgcoojRl8H1uZtPwRuzPfPyjs9BLwELGh47dr8un3ADa36zNsX5D6Gcp9Ty97/Cmb0e+CvwF7gV+SjVqt861I+v6ZYlvwPRcX+rbx9FrAFeC3/nFn2/lcsn0fz/NkDPE3DP+pVvnU6I+ALFMu1e4Bd+bbcOXRS+TiHim2XATtzDnuBHzS0H/dnmWcWlCSpxqp8jIAkSeoyCwFJkmrMQkCSpBqzEJAkqcYsBCRJqjELAUmSasxCQJKkGrMQkCSpxiwEJI1bRHw+IrZExL8jYigilkbEzRHxp7LHJml8LAQkjUtEXElxbvPnKE51+mfgHopToX6/xKFJmgBPMSxpXCJiK3AkpfSN/Hg58AywNaV0bamDkzRup7dvIkmFiPgk8EWg8QP/fYrVRVcDpEnIrwYkjcdn8s/tDdsuAfallF4oYTySTpGFgKTxmEFxidgTABFxDsWxAe+VOShJE2chIGk8dgEBrImIS4BB4AhwcUQsLHVkkibEQkDSSUspHaBYAVgF7AbeBb4M7AVeLHFokibIvxqQJKnGXBGQJKnGLAQkSaoxCwFJkmrMQkCSpBqzEJAkqcYsBCRJqjELAUmSasxCQJKkGrMQkCSpxv4LMBwyWtSTm/MAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x12075b828>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "α_vals = []\n",
    "dθdt_vals = []\n",
    "\n",
    "for α in np.arange(0.0001, 0.003, 0.0002):\n",
    "    dθdt = compute_precession(α)\n",
    "    α_vals.append(α)\n",
    "    dθdt_vals.append(dθdt)\n",
    "\n",
    "# make a linear best fit\n",
    "p = np.polyfit(α_vals, dθdt_vals, 1)\n",
    "\n",
    "fig, axes = plt.subplots(figsize=(8,6))\n",
    "plt.plot(α_vals, np.polyval(p, α_vals), 'r-')\n",
    "plt.plot(α_vals, dθdt_vals, 'ko')\n",
    "plt.xlabel(r'$\\alpha$', fontsize=14)\n",
    "plt.ylabel(r'd$\\theta$/dt (deg/yr)', fontsize=14)\n",
    "#plt.ylim(ymin=0)\n",
    "plt.title(r'$d\\theta/dt = C \\alpha + b$', fontsize=16)\n",
    "\n",
    "plt.text(0.0005, 30, '$C$ = {:.2e}'.format(p[0]), fontsize=16)\n",
    "plt.text(0.0005, 26, '$b$ = {:.2e}'.format(p[1]), fontsize=16)\n",
    "\n",
    "C = p[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, if this is going to work out, then need to calculate\n",
    "\n",
    "$$ \\frac{\\mathrm{d}\\theta}{\\mathrm{d}t} = C \\alpha $$\n",
    "\n",
    "where $$\\alpha = 1.1 \\times 10^{-8}\\;\\mathrm{AU}^2$$ is the value from general relativity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.000122397563007\n"
     ]
    }
   ],
   "source": [
    "α = 1.1e-8\n",
    "dθdt = C * α\n",
    "print (dθdt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is units of degrees / year. In arcseconds / century we have:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Precesion of rate due to general relativity is 44.1 arcseconds/century\n"
     ]
    }
   ],
   "source": [
    "print ('Precesion of rate due to general relativity is {:.1f} arcseconds/century'.format(dθdt * 3600 * 100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "as expected."
   ]
  }
 ],
 "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}