runge_kutta.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import packages\n",
    "import numpy as np\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import astropy.units as u\n",
    "import astropy.constants as const"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Two-body dynamics: setting up Runge-Kutte method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(rs, **kwargs):\n",
    "    \"\"\"\n",
    "    Evaluate the second derivative of the position vector in 2D.\n",
    "    \n",
    "    Inputs:  array with x-, y- components of the position vector\n",
    "    Returns: array of x-,y- components of acceleration \n",
    "             evaluated at these points\n",
    "    \"\"\"\n",
    "        \n",
    "    # position of mass1, mass2    \n",
    "    x1, y1 = rs[0], rs[1]  \n",
    "    x2, y2 = rs[2], rs[3] \n",
    "    r = np.sqrt((x2 - x1)**2 + (y2 - y1)**2)\n",
    "    \n",
    "    # x, y components of Acceleration a(t) = x''(t), which is also v'(t)\n",
    "    ax1 = G*m2*(x2 - x1)/r**3    \n",
    "    ay1 = G*m2*(y2 - y1)/r**3    \n",
    "    ax2 = G*m1*(x1 - x2)/r**3    \n",
    "    ay2 = G*m1*(y1 - y2)/r**3    \n",
    "    \n",
    "    return np.array([ax1, ay1, ax2, ay2], float)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_f(r, f):\n",
    "    \"\"\"\n",
    "    Calculate updated x, y, vx, vy for each particle using a function \n",
    "        f(x,y,t) = v'(t)\n",
    "    \n",
    "    Inputs:\n",
    "        r, 1-D array [x1,y1,x2,y2,vx1, vy1, vx2, vy2]\n",
    "        \n",
    "    Returns:\n",
    "        fxy, 1-D array v(t), f(x,y,t)\n",
    "        \n",
    "    \"\"\"\n",
    "    \n",
    "    x1, y1 = r[0], r[1]\n",
    "    vx1, vy1 = r[4], r[5]\n",
    "    x2, y2 = r[2], r[3]\n",
    "    vx2, vy2 = r[6], r[7]\n",
    "    \n",
    "    fxy = np.array([vx1, vy1, vx2, vy2])\n",
    "    fv = f(np.array([x1, y1, x2, y2]))\n",
    "    \n",
    "    return np.array(np.concatenate([fxy, fv]))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_energy(r, vi):\n",
    "    \"\"\"\n",
    "    Calculate vis-viva energy for two masses\n",
    "    \n",
    "    Inputs:\n",
    "        r, 1-D array of (x1, y1, x2, y2)\n",
    "        vi, 1-D array of (vx1, vy1, vx2, vy2)\n",
    "    Returns:\n",
    "        energy, v**2/2 - mu/r\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    x1, y1, x2, y2 = r\n",
    "    vx1, vy1, vx2, vy2 = vi\n",
    "    \n",
    "    r12 = np.sqrt((x2 - x1)**2 + (y2 - y1)**2)\n",
    "    v12 = np.sqrt((vx2 - vx1)**2 + (vy2 - vy1)**2)\n",
    "    \n",
    "    vve = 0.5*v12**2 - mu/r12  \n",
    "    E = 0.5*mu*v12**2 - G*m1*m2/r12\n",
    "    \n",
    "    return E\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_err(truth, computed):\n",
    "    \"\"\"Calculate relative error\"\"\"\n",
    "    if truth == 0.:\n",
    "        print(\"using alt err\")\n",
    "        err = truth - computed\n",
    "    else: \n",
    "        err = (truth - computed)/truth\n",
    "    return err"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rk4(t0, tmax, r0, phi0, dphi0, dr0, f, h, err_tol):\n",
    "    \"\"\"\n",
    "    Calculate orbit of a body about center of mass using 4th-order Runge-Kutte\n",
    "    method.\n",
    "    \n",
    "    Inputs:\n",
    "      t0, tmax: start, end time of simulation, Myr\n",
    "      r0, phi0: initial position, angle of reduced mass, mu\n",
    "      dphi0, dr0: initial angular, radial velocity of mu\n",
    "      f: function for getting vx', vy'\n",
    "      h: time interval\n",
    "      err_tol: amount of allowed deviation from r_ref\n",
    "    \n",
    "    Returns:\n",
    "      R: 2-D array of x1,y1,x2,y2 coordinates at each step\n",
    "      err: relative error of total Energy of system\n",
    "      \n",
    "    \"\"\"\n",
    "    \n",
    "    # Initialize arrays for time, position vector, & velocity vector\n",
    "    ti = t0\n",
    "    t = np.array(ti)\n",
    "    \n",
    "    x0 = r0 * np.cos(phi0)\n",
    "    y0 = r0 * np.sin(phi0)\n",
    "    vx0 = -r0 * np.sin(phi0)*dphi0 + np.cos(phi0)*dr0\n",
    "    vy0 = r0 * np.cos(phi0)*dphi0 + np.sin(phi0)*dr0\n",
    "    \n",
    "    # convert from one-body reduced-mass to two-body problem\n",
    "    x1, y1 = (mu/m1)*x0, (mu/m1)*y0 \n",
    "    x2, y2 = -(mu / m2)*x0, -(mu / m2)*y0 \n",
    "    vx1, vy1 = (mu/m1)*vx0, (mu/m1)*vy0 \n",
    "    vx2, vy2 = -(mu / m2)*vx0, -(mu / m2)*vy0 \n",
    "     \n",
    "    # initialize arrays of positions at each step\n",
    "    ri = np.array([x1, y1, x2, y2], float)\n",
    "    \n",
    "    # initialize arrays of velocities at each step\n",
    "    vinit = np.array([vx1, vy1, vx2, vy2], float)\n",
    "    vi = vinit\n",
    "    v_ = np.array([vi]) # Initial vx, vy\n",
    "    \n",
    "    # initialize position-velocity mega-vector\n",
    "    r = np.concatenate([ri,vi])\n",
    "    R = np.array([r])\n",
    "    \n",
    "    debug = False\n",
    "    \n",
    "#     # Initialize steps for error check for circular orbit\n",
    "#     r_check = np.sqrt(r[0]**2 + r[1]**2)  # current radius, r\n",
    "#     rel_err = (r_ref - r_check)/r_ref\n",
    "\n",
    "    # track relative error between initial Energy, calculated Energy \n",
    "    E_now = get_energy(ri, vi)\n",
    "    # total energy of system, E = KE + PE\n",
    "    E_ref = 0.5*mu*(vx0**2 + vy0**2) - G*m1*m2/r0\n",
    "    \n",
    "    rel_err = get_err(E_ref, E_now)\n",
    "    #print(f\"t=0: E_ref: {E_ref}, rel_err: {rel_err}\")\n",
    "    err = np.array([rel_err])\n",
    "    global E\n",
    "    E = np.array([E_now])\n",
    "    \n",
    "    # solve eqns of motion using R-K method\n",
    "    while np.abs(rel_err) < err_tol:\n",
    "        # break if t = tmax\n",
    "        if ti > tmax:\n",
    "            if debug: print(f\"Time: {ti} > {tmax}\")\n",
    "            break\n",
    "        else:\n",
    "            # calculate r(t), r'(t)\n",
    "            \n",
    "            # compute k1, k2, k3, k4 for r\n",
    "            k1  = h*get_f(r, f)\n",
    "            k2 = h*get_f(r + k1*0.5, f) \n",
    "            k3 = h*get_f(r + k2*0.5, f) \n",
    "            k4 = h*get_f(r + k3, f)  \n",
    "\n",
    "            # update ti, Ri\n",
    "            ti += h\n",
    "            r += (1/6)*(k1 + 2*k2 + 2*k3 + k4)\n",
    "\n",
    "            # calculate new rel error\n",
    "            rel_err = get_err(E_ref, get_energy(r[:4], r[4:]))\n",
    "\n",
    "            # save new ri to r\n",
    "            t = np.append(t, ti)\n",
    "            R = np.append(R, np.array([r]), axis=0)\n",
    "            err = np.append(err, rel_err) # relative error\n",
    "            E = np.append(E, get_energy(r[:4], r[4:]))  # total energy\n",
    "            \n",
    "    # components of cartesian coords\n",
    "    # R: [x1, x2, y1, y2]\n",
    "    X = R[:,[0,2]]  # x1, x2\n",
    "    Y = R[:,[1,3]]  # y1, y2\n",
    "    # show error\n",
    "    print(f\"E_ref = {E_ref}, |mean(E)| = {np.abs(np.mean(E))}\")\n",
    "    \n",
    "    return t, X, Y, err\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "parabolic orbit\n",
      "initial r0, phi0: 50.0 kpc, 1.5708 rad\n"
     ]
    }
   ],
   "source": [
    "# runtime initial conditions\n",
    "t0, tmax = 0, 800 # 5000\n",
    "h = 1 # 1e-3 Myr\n",
    "\n",
    "# UNITS: kpc, Msun, Myr\n",
    "# constants\n",
    "G = (const.G).to(u.kpc**3 * u.Msun**-1 * u.Myr**-2).value # kpc^3/Msun/Myr\n",
    "\n",
    "# orbit parameters\n",
    "m1 = 1e12 # Msun\n",
    "m2 = 1e12 # Msun\n",
    "M = m1 + m2\n",
    "mu = m1*m2 / M\n",
    "\n",
    "# can get a = (G*M*T**2/(4*np.pi**2))**(1/3)\n",
    "rmin = 25 # kpc\n",
    "a = rmin   # semi-major axis (kpc, Kepler's Law)\n",
    "ecc = 1. # 0.5\n",
    "\n",
    "# define values for p, error calculation\n",
    "if ecc == 1.:      # parabolic orbit\n",
    "    print(\"parabolic orbit\")\n",
    "    p = 2 * rmin  \n",
    "    err_tol = 1e10  # arbitrary value because E = 0\n",
    "elif (ecc > 0.) & (ecc < 1.): # closed orbit, elliptical\n",
    "    print(\"elliptical orbit\")\n",
    "    p = a * (1 - ecc**2)\n",
    "    err_tol = 1e-5\n",
    "else: # closed orbit, circular\n",
    "    print(\"circular orbit\")\n",
    "    p = a\n",
    "    err_tol = 1e-5\n",
    "\n",
    "\n",
    "# choose initial position, motion of reduced-mass galaxy pair relative to \n",
    "#    semi-major axis or axis of closest approach\n",
    "phi0 = 0.5*np.pi # np.pi*(2*np.random.rand() - 1) \n",
    "r0 = p / (1 + ecc * np.cos(phi0))\n",
    "dphi0 = np.sqrt(G*M*p) / r0**2\n",
    "dr0 = ecc * r0**2 * np.sin(phi0) * dphi0 / p\n",
    "\n",
    "print(f\"initial r0, phi0: {r0:.5} kpc, {phi0:.5} rad\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E_ref = 1.52587890625e-05, |mean(E)| = 6.37734529205094\n",
      "runtime: 801 kpc\n"
     ]
    }
   ],
   "source": [
    "t, X, Y, err = rk4(t0, tmax, r0, phi0, dphi0, dr0, f, h, err_tol)\n",
    "print(f\"runtime: {t[-1]} kpc\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAX8ElEQVR4nO3dfZAc913n8fd3nvb5Sfugx5VXihQHxVFssXZsAglJlEdyNiFAKVWBUHcgSGGKXI6HGFcBqeOufCEkVxSQIMDcwZmYEOIQnAQ/hCQQU3a8tiVbliVLsmVZj7taSSvtrrQPs9/7Y2atkbSrfZie7enuz6tqamZ6Zro/Wq0+89NvunvM3RERkehKhR1ARETKoyIXEYk4FbmISMSpyEVEIk5FLiIScZkwNtrR0eE9PT1hbFpEJLKeeuqpU+7eeeXyUIq8p6eHvr6+MDYtIhJZZvbKTMs1tSIiEnEqchGRiFORi4hEnIpcRCTiVOQiIhEXSJGb2fvMbJ+ZHTCzTwWxThERmZ+yi9zM0sCfAu8HNgEfMbNN5a5XRETmJ4gR+S3AAXd/yd3HgfuBOwJY71W+/cJJ/uy7ByqxahGRyAqiyFcDr5bcP1Jcdhkz225mfWbWNzAwsKgN/fv+U3zxuwcXl1JEJKaCKHKbYdlV31bh7jvcvdfdezs7rzrCdF6aajMMj02iL8MQEbkkiCI/AnSX3F8DHAtgvVdpqs0w5TAynq/E6kVEIimIIn8S2Ghm68wsB2wDvh7Aeq/SVJsF4PzFiUqsXkQkksoucnefBO4EHgJeAL7s7s+Xu96ZNNUWzvF1/uJkJVYvIhJJgZz90N2/CXwziHVdi0bkIiJXi9SRndMj8nMXNCIXEZkWqSJvni5yjchFRF4TqSK/NLWiEbmIyLRIFXmzilxE5CqRKvLabIpMyvRhp4hIiUgVuZnRVJvRiFxEpESkihwK8+QakYuIXBLBIteIXESklIpcRCTiIljkWe1HLiJSIoJFrhG5iEipyBV5s0bkIiKXiWCRF75cYmpKXy4hIgIRLPKm2izuMDKu6RUREYhkkeuc5CIipSJY5DrfiohIqQgW+fSIXB94iohApItcI3IREYhkkRemVrQLoohIQeSK/NK3BGlELiICUSzyuuKI/IJG5CIiUGaRm9nPmNnzZjZlZr1BhbqW2myamkxKRS4iUlTuiHw38FPAvwWQZd5a6rIMqchFRADIlPNid38BCt/cs5RU5CIilyzZHLmZbTezPjPrGxgYKGtdLXVZzo6qyEVEYB5FbmaPmtnuGS53LGRD7r7D3Xvdvbezs3PxidGIXESk1JxTK+6+dSmCLERLfZa9J86HHUNEpCpEbvdDKIzItdeKiEhBubsffsjMjgC3Ad8ws4eCiXVtLXVZzo9NMpmfWorNiYhUtXL3WnkAeCCgLPPWMn1Q0MVJljXklnrzIiJVJbJTK4A+8BQRQUUuIhJ5kSzy1noVuYjItEgWuUbkIiKXRLLIp8+AODQ6HnISEZHwRbLINSIXEbkkkkVek0lTm02pyEVEiGiRA7TW5VTkIiJEuMh14iwRkYJIF7lOZSsiEuEib9aIXEQEiHCR6wyIIiIFkS3y1vosZ1XkIiLRLfK2+iyj43nGJvNhRxERCVV0i7x4+lp94CkiSRfdIq8vFPnpER2mLyLJFvkiP6MiF5GEi2yRT38z0BlNrYhIwkW2yNsaCifOOq0zIIpIwkW2yFvrNLUiIgIRLvJcJkVTTUYfdopI4pVV5Gb2h2a218yeNbMHzKw1oFzz0taQ46ymVkQk4codkT8C3ODum4EXgbvKjzR/bQ05TuvDThFJuLKK3N0fdvfJ4t3HgTXlR5q/tvqs5shFJPGCnCP/z8C3ZnvQzLabWZ+Z9Q0MDASywWX1Oc2Ri0jiZeZ6gpk9CqyY4aG73f2fis+5G5gE7pttPe6+A9gB0Nvb64tKewXNkYuIzKPI3X3rtR43s48BHwTe5e6BFPR8tdVnGRnPc3EiT202vZSbFhGpGuXutfI+4LeB2919NJhI86cTZ4mIlD9H/idAE/CIme00sy8GkGnelunEWSIic0+tXIu7bwgqyGJcGpGryEUkuSJ7ZCeUnMpWRS4iCRbtIi+eOEv7kotIkkW7yIsj8kEVuYgkWKSLPJtO0VKXZXBYRS4iyRXpIgfoaMxxangs7BgiIqGJQZHXqMhFJNGiX+RNNZpaEZFEi36RN+QY0IhcRBIs+kXeWMP5i5NcnMiHHUVEJBTRL/KmGkC7IIpIckW/yBsLRX7qvKZXRCSZIl/k7Y3TBwWpyEUkmSJf5J2vjcg1tSIiyRT5Ip+eWtGeKyKSVJEv8rpcmoZcWgcFiUhiRb7IAdobdVCQiCRXLIpc51sRkSSLSZHrfCsiklzxKPKmGk5pakVEEioeRd5Yw5nRcSbyU2FHERFZcmUVuZn9dzN71sx2mtnDZrYqqGALsby5Bnc0vSIiiVTuiPwP3X2zu98IPAj8bvmRFm55Uy0AJ8+pyEUkecoqcnc/V3K3AfDy4izOipZCkZ8YuhjG5kVEQpUpdwVm9j+AnweGgHeUnWgRupoLR3f2n1eRi0jyzDkiN7NHzWz3DJc7ANz9bnfvBu4D7rzGerabWZ+Z9Q0MDAT3JwA6GmpIp4yT51TkIpI8c47I3X3rPNf1d8A3gN+bZT07gB0Avb29gU7BpFJGV1MNJ4Y0Ry4iyVPuXisbS+7eDuwtL87idTXXampFRBKp3Dnye8zsemAKeAX4lfIjLc6K5hpePjUS1uZFREJTVpG7+4eDClKu5c21PP7S6bBjiIgsuVgc2QmFIh+6MKEvYRaRxIlVkQPac0VEEidGRV7Yl1xHd4pI0sSmyFcUR+QnNCIXkYSJTZF3FYu8X0UuIgkTmyJvrs1Qn0tz7KyKXESSJTZFbmasaq3j2NkLYUcREVlSsSlygNWtdRwbUpGLSLLEqshXtdZx9IyKXESSJVZFvqatjsGRcR0UJCKJEqsiX9Va2HPlqObJRSRBYlXkq1vrATS9IiKJEqsinx6Ra88VEUmSWBX5iuZaUqapFRFJllgVeSadYkVzrYpcRBIlVkUOsLpNuyCKSLLErshXtdZpRC4iiRK7Il/dWseJoYvkpwL9fmcRkaoVuyLvXlbP5JRzXIfqi0hCxK7Ir1tW2Jf8lcHRkJOIiCyN+BV5RwOgIheR5AikyM3sN8zMzawjiPWVY0VzLbl0ilcGR8KOIiKyJMoucjPrBt4NHC4/TvnSKaN7WZ1G5CKSGEGMyD8P/BZQNbuJ9LQ3cEgjchFJiLKK3MxuB466+655PHe7mfWZWd/AwEA5m53T2vZ6Dp8exb1q3ltERComM9cTzOxRYMUMD90N/A7wnvlsyN13ADsAent7K9qwPe0NjI7nGRgeo6uptpKbEhEJ3ZxF7u5bZ1puZm8C1gG7zAxgDfC0md3i7icCTblAa9sLuyAeHhxVkYtI7C16asXdn3P3Lnfvcfce4AiwJewSh8KIHOCQPvAUkQSI3X7kUDhMP2VoF0QRSYQ5p1bmqzgqrwq5TIruZfW8NKAiF5H4i+WIHGBjVyP7+8+HHUNEpOJiW+Sv62rk5VMjTOanwo4iIlJRsS3yDZ2NTOSdw6f1gaeIxFtsi3zj8iYADvQPh5xERKSyYlvkr+ss7IJ4YEBFLiLxFtsib6rNsqK5lgMnVeQiEm+xLXKADV2NGpGLSOzFvsgP9g/r5FkiEmuxLvKNyxsZGc9z5Iy+v1NE4ivWRf5DK5sBeOH4uZCTiIhUTqyL/A0rmjCDPSpyEYmxWBd5fS7Duo4G9hxTkYtIfMW6yAE2rWzWiFxEYi3+Rb6qmSNnLjB0YSLsKCIiFRH/Ii9+4LlXo3IRian4F/mqQpFrekVE4ir2Rd7VVEtHYw3PHR0KO4qISEXEvsgBbuxuYderZ8OOISJSEYko8pvWtnFwYIShUX3gKSLxk4giv7G7FYCdR86GmkNEpBISUeSb17RgBjsPnw07iohI4MoqcjP7fTM7amY7i5cPBBUsSE21WTZ2NfLMq2fCjiIiErhMAOv4vLt/NoD1VNRN3W08tOcE7o6ZhR1HRCQwiZhaAbhpbStnRyc4ODASdhQRkUAFUeR3mtmzZnavmbXN9iQz225mfWbWNzAwEMBmF+Yt69sBeOLlwSXftohIJc1Z5Gb2qJntnuFyB/AF4HXAjcBx4I9mW4+773D3Xnfv7ezsDCr/vPW017O8uYbHXzq95NsWEamkOefI3X3rfFZkZn8BPFh2ogoxM25b3873DwxqnlxEYqXcvVZWltz9ELC7vDiVdev6dk4Nj2meXERipdy9Vj5jZjcCDhwCfrncQJV02+sK8+SPvzTIhq7GkNOIiASjrCJ3958LKshSWLusnlUttTx24BQfvfW6sOOIiAQiMbsfQmGe/O3Xd/L9/aeYyE+FHUdEJBCJKnKAH7++i/Njk/Qd0lGeIhIPiSvyH93QQTZtfGdff9hRREQCkbgib6jJ8JZ17Xxnr4pcROIhcUUO8I43dLG/f5jDg6NhRxERKVsii/w9m5YD8M3dx0NOIiJSvkQWefeyet7c3co3nlWRi0j0JbLIAT74ppU8d3SIVwZ1lKeIRFtii/wDmwtnF3hQo3IRibjEFvnq1jp++Lo2HnjmKO4edhwRkUVLbJED/GzvGg70D/O0vstTRCIs0UX+E5tXUZ9L8/dPHg47iojIoiW6yBtrMnxw80oefPY4w2OTYccREVmURBc5wEduWcvoeJ6v9L0adhQRkUVJfJHftLaNLWtbufexQ+Sn9KGniERP4osc4Bd/bD2HT4/yyJ6TYUcREVkwFTnw3jeuoHtZHV/43kHtiigikaMiB9Ip4853bGDXq2d59AWdFVFEokVFXvThLWtY19HAHz28jynNlYtIhKjIizLpFJ/YupG9J87zz88eCzuOiMi8qchL/KfNq9i0spl7vrVX+5WLSGSUXeRm9mtmts/MnjezzwQRKiyplPEHH7qBE+cu8vlHXgw7jojIvJRV5Gb2DuAOYLO7vxH4bCCpQrRlbRsfuWUtf/3Yy+w+OhR2HBGROZU7Iv84cI+7jwG4eyx2+fjt976BZQ01fPLLO7kwng87jojINZVb5K8HfszMnjCz75nZzbM90cy2m1mfmfUNDAyUudnKaqnP8rmffTMvnhzmD76xJ+w4IiLXNGeRm9mjZrZ7hssdQAZoA24FfhP4spnZTOtx9x3u3uvuvZ2dnYH+ISrhba/v5Jfftp77njjM13dpLxYRqV6ZuZ7g7ltne8zMPg581QuHQ/7AzKaADqC6h9zz9N/ecz3PHD7Lb/zDLla11NLbsyzsSCIiVyl3auVrwDsBzOz1QA44VeY6q0Yuk+LPf+6HWd1axy/9TR8HB4bDjiQicpVyi/xeYL2Z7QbuBz7mMTtZSVtDjnt/4WZSZmzb8TgH+lXmIlJdyipydx9394+6+w3uvsXd/zWoYNVkXUcD92+/FXfYtuNx7ZYoIlVFR3bO08blTdy//S3k0sZPf/E/+NZzx8OOJCICqMgXZENXE1+786380MpmPn7f09zzrb2MT06FHUtEEk5FvkBdTbV86ZduZdvN3Xzxewf5yT99jL0nzoUdS0QSTEW+CLXZNPd8eDN/8fO9nDx3kZ/44+/z6X9+nqHRibCjiUgCqcjL8O5Ny3nkk29n283d/J//OMTbP/sd/vjb+xm6oEIXkaVjYewt2Nvb6319fUu+3Up6/tgQn3v4Rb69t5/Gmgzbbu5m2y3dbOhqCjuaiMSEmT3l7r1XLVeRB2vPsXP82XcP8C+7TzA55WxZ28pPbVnDuzctZ3lzbdjxRCTCVORL7NTwGA88fZT7nzzMwYERAN7c3crWN3TxIxvaedPqVnIZzWyJyPypyEPi7uzvH+aRPSd5eM9Jdr16FoDabIqbutu4cW0rb1zVzKaVzfS0N5BKzXjOMRERFXm1GBwe48lDZ3ji5UF+8PJp9p04z2Txy57rc2k2djVyXXsDPe31heuOela31tPRmCOT1gheJMlU5FVqbDLP/pPD7Dl+jj3HznFwYJhDgyMcPXOBqZK/GjNob6hheXMNy5trWd5cQ3tDDa31WZrrsrTWZWmpy9Jan6OlLktzXYa6bJpZziosIhE0W5HPeRpbqayaTJobVrdww+qWy5aPT05x5MwohwZHOHb2Iv3nx+g/V7g+ee4izx4Z4vTI2GVlfyUzqMumqc+lqc9lqM+lqculachlqMulqcumyWVShUv60nW2eDubNmoypfcLy1JmpFMlFzNSKSOTKlynr3h8+vmlj6dSFK7NMAPDsBQYYGYYvPbY9J8lVVxuZqQMvUmJFKnIq1Quk2J9ZyPrOxtnfY67Mzw2ydnRCYYuFC7Tt89dnGB0PM/o2CSjE3kujOcZHZ8sLBvPMzgyzoXxSSbyztjkFOOTeSbyznh+ivy13h2q0OUlf+mNYPoNInXlsiveAErv2mXLbZbll716luWzv8Zmec1s254ta5DrXaxA1xbw+3KQqwvy5/Y/P/QmblkX7HcbqMgjzMxoqs3SVJulO8D15qecifwUY5NTTOSnGL/s2plyJz/lTE5duj1VvJ/3wu389GX6cXcmX3strz3P3ZlycApvTO7gTF/DVHEZUHg+vPacqeKd6WVTJbf9tdtefOzSekuVTi36ZctLbpc8cvnymZ9/5aOXvaaM9TpXhp/tNXP/mYIQ5OqCnuINdG0B/9waatLBrhAVucygMCWSpjYb/C+ciARPu0GIiEScilxEJOJU5CIiEaciFxGJOBW5iEjEqchFRCJORS4iEnEqchGRiAvlpFlmNgC8ssiXdwCnAowTFOVaGOVaGOVamLjmus7dO69cGEqRl8PM+mY6+1fYlGthlGthlGthkpZLUysiIhGnIhcRibgoFvmOsAPMQrkWRrkWRrkWJlG5IjdHLiIil4viiFxEREqoyEVEIi5SRW5m7zOzfWZ2wMw+tcTbvtfM+s1sd8myZWb2iJntL163lTx2VzHnPjN7b4UydZvZd8zsBTN73sx+vUpy1ZrZD8xsVzHXp6shV8m20mb2jJk9WC25zOyQmT1nZjvNrK+KcrWa2VfMbG/x9+y2sHOZ2fXFn9P05ZyZfSLsXMXt/Nfi7/xuM/tS8d9C5XMVvgqr+i9AGjgIrAdywC5g0xJu/23AFmB3ybLPAJ8q3v4U8L+KtzcV89UA64q50xXItBLYUrzdBLxY3HbYuQxoLN7OAk8At4adqyTfJ4G/Ax6shr/H4rYOAR1XLKuGXP8X+MXi7RzQWg25SvKlgRPAdWHnAlYDLwN1xftfBn5hKXJV7AdcgR/SbcBDJffvAu5a4gw9XF7k+4CVxdsrgX0zZQMeAm5bgnz/BLy7mnIB9cDTwFuqIRewBvg28E4uFXk15DrE1UUeai6guVhMVk25rsjyHuCxashFochfBZZR+BrNB4v5Kp4rSlMr0z+kaUeKy8K03N2PAxSvu4rLlzyrmfUAN1EY/Yaeqzh9sRPoBx5x96rIBfxv4LeAqZJl1ZDLgYfN7Ckz214ludYDA8BfF6ei/tLMGqogV6ltwJeKt0PN5e5Hgc8Ch4HjwJC7P7wUuaJU5DbDsmrdd3JJs5pZI/CPwCfc/dy1njrDsorkcve8u99IYQR8i5ndEHYuM/sg0O/uT833JTMsq9Tf41vdfQvwfuBXzext13juUuXKUJhO/IK73wSMUJgaCDtXYWNmOeB24B/meuoMyyrx+9UG3EFhmmQV0GBmH12KXFEq8iNAd8n9NcCxkLJMO2lmKwGK1/3F5UuW1cyyFEr8Pnf/arXkmubuZ4HvAu+rglxvBW43s0PA/cA7zez/VUEu3P1Y8bofeAC4pQpyHQGOFP83BfAVCsUedq5p7weedveTxfth59oKvOzuA+4+AXwV+JGlyBWlIn8S2Ghm64rvxNuAr4ec6evAx4q3P0Zhjnp6+TYzqzGzdcBG4AdBb9zMDPgr4AV3/1wV5eo0s9bi7ToKv+B7w87l7ne5+xp376Hw+/Ov7v7RsHOZWYOZNU3fpjCvujvsXO5+AnjVzK4vLnoXsCfsXCU+wqVplenth5nrMHCrmdUX/22+C3hhSXJV8oOICnyY8AEKe2YcBO5e4m1/icK81wSFd9L/ArRT+OBsf/F6Wcnz7y7m3Ae8v0KZfpTCf8WeBXYWLx+oglybgWeKuXYDv1tcHmquKzL+OJc+7Az757Wewt4Lu4Dnp3+3w85V3M6NQF/x7/JrQFuV5KoHBoGWkmXVkOvTFAYtu4G/pbBHSsVz6RB9EZGIi9LUioiIzEBFLiIScSpyEZGIU5GLiEScilxEJOJU5CIiEaciFxGJuP8PeTami/K2OxsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(E);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "error at 801 Myr: 4.34779e+05 kpc\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtYAAAE9CAYAAAAveyUaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAABNR0lEQVR4nO3dd3yV9fn/8deVQdiyh0QIS5QtBgRxoChKtSK2KFoVJ47aOmqrfP1VrV9na+torV+pA5wsBzhbRHALBMXBKiABIitsjEgY1++P+w4GCCEJ5+Q+Sd7Px+M8zjmfc9/3ue7k5M6VT67P52PujoiIiIiIHJykqAMQEREREakMlFiLiIiIiMSAEmsRERERkRhQYi0iIiIiEgNKrEVEREREYkCJtYiIiIhIDKREHUAsNGrUyDMyMqIOQ0SkTGbNmrXW3RtHHUd50nVbRCqq4q7ZlSKxzsjIICsrK+owRETKxMyWRh1DedN1W0QqquKu2SoFERERERGJASXWIiIiIiIxoMRaRERERCQGlFiLiIiIiMSAEmsRERERkRhQYi0iIiIiEgNKrEVEREREYkCJtYiIiIhIDCixFhERERGJgUhXXjSzG4ErAAe+Bi4FagJjgQwgGzjX3TfE+r1nLV1P7pZtnHJkU1KS9feFiIhIolmfl8//TVvEZ9+uY+XGH3F8v9smmVGvZiob8rbvd7uSbKNjVY1jJZmR3qAmfxnSjbaNa+/3GKVl7vsPKJ7MrAXwEdDR3bea2TjgLaAjsN7d7zezW4H67n5LccfKzMz00i6Ne9PY2bzyxXc0q1ud83u15Pxeh9GkbvUyno2ISNmZ2Sx3z4w6jvJUluu2HNjnSzdw3Yufk79j5+62WCY0Jd0uVsfaun0Xefk7i9hTJDZO6tCYZy7tVap9irtmR9pjHb5/DTPbTtBTvQIYAfQLXx8NTAOKTazL4i9DujGwS3Oe/2wpD737X/7+3kJO69SMC3u3onebBphZrN9SRESkxNbn5fPw5AX8e85qdu7aVeQ2eyemG7buYOeufRPUNd/nl+g9Y7ldLI9Vq1oSNVNTEibh17Eq/rEKeqz/35kd97t/WUSWWLv7d2b2ILAM2Ar8x93/Y2ZN3X1luM1KM2sSj/dPTjJO7diUUzs2JXttHi9MX8q4rBze/Hol7ZrU5qLerRjcowV1q6fG4+1FRET2mzwnmYHBmi2lT2ANaFjrp99dFbnHOl7/rheJl8gSazOrDwwCWgMbgfFmdmEp9h8ODAdo2bLlQcWS0agWt53Rkd8N6MDrX67g+c+WcsekOTzwznzOPqoFFx7Tio6H1j2o9xAREfl86QauH/MFNVKT2JC3nR+27zxgqUOqwSE1i+7k2TsxTUtN5u/n96BHq/rxCF9EDiDKUpBTgCXungtgZq8AxwKrzax52FvdHFhT1M7uPhIYCUGtXiwCqp6azJDMwxiSeRhfLt/I858t5eVZObw4fRmZrepzUZ9WnN65GWkpybF4OxERqeT27pHeX6nG3slzkhmNalejTo1q3HtOF/XWilQQUSbWy4DeZlaToBSkP5AF5AHDgPvD+4lRBNftsHp0O6wet51xJBNm5fD8Z0u5fsxsGtaqxgXHtOTC3q1oqsGOIiJShIJBhJu2bt+nRzoJaNekFhvytmOGSh1EKpEoa6ynm9kE4HNgB/AFQQ90bWCcmV1OkHwPiSpGgHo1q3HF8W24rG9rPlq0lmc/zeYfUxfx+LTFnNG1OZf2bU33w+pFGaKIiCSAwr3Ta/Py9+iZLuiRVqmGSOUW6awg7n4HcMdezdsIeq8TSlKSccLhjTnh8MZkr81j9KfZjM/KYeLsFRzVsh6X9m3NwM7NSNWc2CIiVUpBQj3xyxVs2rpjd7sBjWtXU4+0SBUS9XR7FVJGo1rc8fNO3HTq4UyYlcPoT7L57Utf0KxudS7q04rze7WkQa1qUYcpIiJx9vnSDVw6asYeCXWttCTq1aimnmmRKkiJ9UGoUz2VS/u2ZlifDKYuWMMzH2fzl38v4NEpCzm7ewsu6ZvBkc01m4iISGVT0Ev94ozl7AhLPtKSjfN6tuSGUw9X54pIFaXEOgaSkoz+Rzal/5FNWbh6C898ks0rn+cwNms5vds04Mrj23BShyYkJWnRGRGRim5x7vdc8vQMlm/YCgSDEbu3rKdyDxFRYh1r7ZvW4d7BXfjDaR0YM3M5oz/J5vLRWbRrUpsrj2/NoO4tqJ6q6fpERCqiz5du4FdPfsbW7cFiLtVTknjxyt4q+RARIPhDW+KgXs1qXH1iWz74w0k8MrQ71ZKTuOXlrznugan8472FbMgr2VKvIiKSGD5fuoFzn/h0d1J9WP0avHn98UqqRWQ39VjHWWpyEoO6t+Csbofy6eJ1PPHBtzz4n//y2NTFnJuZzuXHtaFlw5pRhykiIsVYnPs9v3ryM3bscpINfnVMK9VSi8g+lFiXEzPj2HaNOLZdIxas2sK/PvyWF2cs47nPlnJ652ZceXwbjmqpXg8RkURTuPyjRmoSL1yh0g8RKZoS6wh0aFaHB4d04/endWDUJ9m88NlS3vp6FT0z6nP1iW05+YgmmGmgo4hI1Ap6qguS6jd+e7wGKIrIfqnGOkJN61bnltOP4NMR/bnj5x1ZsfFHLh+dxcBHPmTi7O/YsXNX1CGKSAViZslm9oWZvRE+b2Bmk81sYXhfv9C2I8xskZktMLPTCrUfbWZfh689auFf+WaWZmZjw/bpZpZRaJ9h4XssNLNh5XjKcbU493t+8c9P9uipVlItIsVRYp0AaqWlcGnf1kz7fT/+dm43duxyrh8zm/5/e5+XZixj246dUYcoIhXD9cC8Qs9vBaa4e3tgSvgcM+sIDAU6AacD/zSzgumKHgeGA+3D2+lh++XABndvBzwEPBAeqwHBCrrHAL2AOwon8BXV+rx8Lh81k41btwNwUe9WKv8QkQNSYp1AUpOTOKdHOv+54QSeuOho6tVIZcQrX3PCn6fy5Iffkrdtx4EPIiJVkpmlA2cATxZqHgSMDh+PBs4u1D7G3be5+xJgEdDLzJoDdd39U3d34Nm99ik41gSgf9ibfRow2d3Xu/sGYDI/JeMV1v+9v5jsdT8AcEzr+lzdr13EEYlIRaDEOgElJRmndWrGa7/uy/OXH0PbxrW5+8159H3gPR5+979s/EFT9YnIPh4G/gAUriFr6u4rAcL7JmF7C2B5oe1ywrYW4eO92/fYx913AJuAhsUcq8JanPs947KWAdCqQU0evzBTs3+ISIkosU5gZsZx7Rvx4pW9eeXaY8ls1YCH311I3/vf49635pG7ZVvUIYpIAjCzM4E17j6rpLsU0ebFtJd1nz3f1Gy4mWWZWVZubm6JAo3CHRPnsPGHHdSvmcrTl/ZUUi0iJabEuoLo0bI+Tw7L5N83nMCpHZvy5Iffcvyf3+PuN+ayZsuPUYcnItHqC5xlZtnAGOBkM3seWB2WdxDerwm3zwEOK7R/OrAibE8von2PfcwsBTgEWF/Msfbh7iPdPdPdMxs3bly2M42zxbnfs3RdHgBDjk7XYEURKRUl1hVMh2Z1eHjoUUz5XT9+1qU5T3+8hBP+PFUJtkgV5u4j3D3d3TMIBiW+5+4XApOAglk6hgETw8eTgKHhTB+tCQYpzgjLRbaYWe+wfvrivfYpONYvw/dw4N/AADOrHw5aHBC2VUh3TJzD8g1byWhYU3XVIlJqmse6gmrdqBZ/O7c7vzm5Pf94bxFPf7yE5z5byoW9W3HViW1oUqd61CGKSPTuB8aZ2eXAMmAIgLvPMbNxwFxgB/Brdy+YfugaYBRQA3g7vAE8BTxnZosIeqqHhsdab2b/C8wMt7vL3dfH+8TioXBv9YCOTVUCIiKlZkGHQ8WWmZnpWVlZUYcRqSVr8/jHe4t4bfZ3pCSZEmyRCsTMZrl7ZtRxlKdEvG5f+OR0Plq0loyGNXnl2r5KrEWkSMVds1UKUkm0blSLv57bjSk3ncjPux3KqE+yOf6Bqdz1+lzWfq9BjiIixVFvtYjEghLrSiajUS0eHFI4wQ5qsP/6nwVs/nF71OGJiCQk1VaLSCwosa6kChLsyTedyElHNOHv7y3i+Aem8n/vL2ZrvlZyFBEpsD4vn/xwhVv1VovIwVBiXcm1bVybxy7owRu/OY6jWtbj/rfnc8JfpvLcp9nk79h14AOIiFRyoz/JZkb2Bo5r10i91SJyUCJNrM2snplNMLP5ZjbPzPqYWQMzm2xmC8P7+lHGWFl0bnEIoy7txbir+pDRsCZ/nDiH/n+bxiuf57BzV8UfwCoiUhbr8/L57Nu1AHRsXle91SJyUKLusX4EeMfdjwC6AfOAW4Ep7t4emBI+lxjp1boB467qwzOX9qRu9VRuGvclAx/5gCnzVlMZZogRESmN0Z9kM33JBgBqVIv6V6KIVHSRXUXMrC5wAsHcqLh7vrtvBAYBo8PNRgNnRxFfZWZmnNShCa9fdxyPXdCD7Tudy0dnccG/pvPNd5uiDk9EpFysz8snKzuYcvu4dg0ZdmzriCMSkYouyj/P2wC5wDNm9oWZPWlmtYCm4epfhPdNIoyxUktKMs7o2pz/3HgCdw3qxILVWzjz7x9xw5gvyNnwQ9ThiYjE1ehPsvl48TqOa9eIR8/voTIQETloUSbWKUAP4HF3PwrIoxRlH2Y23MyyzCwrNzc3XjFWCanJSVzcJ4Npv+/Htf3a8vY3qzj5r+9z39vz2LRVU/SJSOVUMEOSaqtFJFaiTKxzgBx3nx4+n0CQaK82s+YA4f2aonZ295HununumY0bNy6XgCu7utVT+cPpRzD15n6c2bU5Iz/4ln5/mcozHy/RDCIiUqmsz8tnzoqg9E211SISK5FdTdx9FbDczDqETf2BucAkYFjYNgyYGEF4Vdqh9Wrwt3O78/p1x3Fk87r86fW5nPbwB7w3f3XUoYmIxEThMhDVVotIrET9Z/pvgBfM7CugO3AvcD9wqpktBE4Nn0sEOrc4hBeuOIZnLumJGVw2KotLnpnB4tzvow5NROQgBbMgHd2qnspARCRmUqJ8c3efDWQW8VL/cg5F9sPMOOmIJvRt14hnP83mkXcXctpDH3Bp3wx+2789daqnRh2iiEipndW9BV/lbOKs7i2iDkVEKpGoe6ylgqiWksQVx7fhvZv78Yse6Tz50RJOevB9xmUtZ5cWmBGRCmR9Xj53vzGXqQtyeXeuStxEJHaUWEupNK6TxgO/7MrEX/elZYMa/GHCVwz+58d8vmxD1KGJiJTI6E+ymbogl+PaNWJI5mFRhyMilYgSaymTrun1ePmaY3n4vO6s2vwj5/zzE0a88hUbf8iPOjQRkWJpmj0RiRcl1lJmZsbZR7Xgvd/1Y/gJbRiXlUP/v77Py7NytDy6iCQkTbMnIvGkq4octFppKfzPz47kjd8cR6uGNfnd+C8ZOvIzFq3ZEnVoIiJ7GJ+1nI8Xr+OkDo01zZ6IxJwSa4mZI5vXZcLVx3LfOV2Yv2oLAx/5kL/8e/7uf7uKiETtlI5NOalDY/7fmR1VBiIiMafEWmIqKck4v1dLpvzuRH7e7VAem7qYAQ+/z4cLtey8iERv0uwVTF2Qy6TZK6IORUQqISXWEheNaqfxt3O789KVvUlNTuKip2bw+/FfsumH7VGHJiJVmu91LyISO0qsJa76tG3IW789nmv7teWVL77j1Ife5z9zVkUdlohUQevzglmLru/fXvXVIhIXSqwl7qqnJvOH049g4q/70rB2GsOfm8V1L37Ouu+3RR2aiFQhoz/J5pEpiwBUXy0icaHEWspN5xaHMOm6vvzu1MP5z5zVnPK395k4+ztNzSci5URlICISX0qspVylJifxm/7tefO3x9GqYS2uHzOba57/fPe/aEVE4uWs7i04qUNjzureIupQRKSSUmItkWjftA4vX3Mstw48gvfmr2HAQx/w3vzVUYclIpXU+rx87n5jLlMX5PLuXF1rRCQ+lFhLZJKTjKtPbMvE6/rSqHY1LhuVxYhXviZv246oQxORSmZ81nKmLsjlpA6NGZJ5WNThiEglpcRaIndk87pMvK4vV53QhjEzl/GzRz9k1tINUYclIpWIFoYRkfKgxFoSQlpKMiN+diRjruzNjp3OkP/7hL/9ZwE7du6KOjQRqQTenbtaZSAiEndKrCWhHNOmIe/ccDyDj0rn0fcWccG/prNy09aowxKRCq6gx/qUjk2jDkVEKjEl1pJw6lRP5a/nduOh87oxZ8UmBj7yoXqZRKTMNHBRRMqLEmtJWIOPSuf13xxHi3o1uOLZLO56fS7bduyMOiwRqWA0cFFEyosSa0lobRrX5pVrj+WSYzN4+uMl/PLxT1m27oeowxKRCkQDF0WkvCixloSXlpLMnWd1YuRFR7N0XR5n/v1Dpi5YE3VYIlJBaOCiiJQXJdZSYQzo1Iw3fnM8LerX5LJRM3nk3YXs2qWliUWkeBq4KCLlJfLE2sySzewLM3sjfN7AzCab2cLwvn7UMUriaNmwJq9ccyyDu7fgoXf/y5XPZrFp6/aowxKRBKYeaxEpL5En1sD1wLxCz28Fprh7e2BK+FxktxrVkvnrud3401mdeP+/uZz1j4+Yv2pz1GGJSAJan5fPD/k7uL5/ew1cFJG4izSxNrN04AzgyULNg4DR4ePRwNnlHJZUAGbGsGMzGDO8N1vzd3LOPz9Rb5SI7GN81nIembKImtWSNXBRROIu6h7rh4E/AIWX12vq7isBwvsmRe1oZsPNLMvMsnJzc+MeqCSmzIwGTLruONo2rs2Vz2Xx5Iff4q66axEJqL5aRMpTZIm1mZ0JrHH3WWXZ391Hunumu2c2btw4xtFJRdLskOqMu6oPp3dqxt1vzmPEK1+Tv0NLoYuI6qtFpHxF2WPdFzjLzLKBMcDJZvY8sNrMmgOE95pXTQ6oRrVkHrugB9ed1I4xM5dz8dPT2fhDftRhiUjE1GMtIuUpssTa3Ue4e7q7ZwBDgffc/UJgEjAs3GwYMDGiEKWCSUoybj6tAw+d143Pl27kF49/wncbt0YdlohESD3WIlKeoq6xLsr9wKlmthA4NXwuUmKDj0rnuct7sWbLNs7558csWLUl6pBEJCJDMg9jxMAjNCOIiJSLhEis3X2au58ZPl7n7v3dvX14vz7q+KTiOaZNQ8Zf3QeAIf/3CTOW6GMkUtWsz8tnfNZyhmQephlBRKRcJERiLRIPRzSry8vXHEujOmlc9NR0/j1nVdQhiUg5Gp+1nPvens/4rOVRhyIiVYQSa6nU0uvXZMLVx3Jk87pc8/wsXv0iJ+qQRKScaOCiiJQ3JdZS6TWoVY0XrzyG3m0actO4Lxmn3iuRKkEDF0WkvCmxliqhZrUUnr6kJ8e3b8wfJnzFC9OXRh2SiMSRljIXkSgosZYqo3pqMiMvOpr+RzThtle/YdTHS6IOSSQmzKy6mc0wsy/NbI6Z/Slsb2Bmk81sYXhfv9A+I8xskZktMLPTCrUfbWZfh689amYWtqeZ2diwfbqZZRTaZ1j4HgvNbBgJQEuZi0gUlFhLlVI9NZnHLzya0zs1487X5/Lsp9lRhyQSC9uAk929G9AdON3MegO3AlPcvT0wJXyOmXUkWD+gE3A68E8zSw6P9TgwHGgf3k4P2y8HNrh7O+Ah4IHwWA2AO4BjgF7AHYUT+Khomj0RiYISa6lyqqUk8fcLjuLUjk25feIcJszSgEap2Dzwffg0Nbw5MAgYHbaPBs4OHw8Cxrj7NndfAiwCeoWr3dZ190/d3YFn99qn4FgTgP5hb/ZpwGR3X+/uG4DJ/JSMi4hUKUqspUpKTU7i7+cfxXHtGvGHCV/y9tcrow5J5KCYWbKZzQbWECS604Gm7r4SILxvEm7eAig8ijcnbGsRPt67fY993H0HsAloWMyxIqWp9kQkCkqspcqqnprMyIuP5qiW9fntmC+YtmBN1CGJlJm773T37kA6Qe9z52I2t6IOUUx7WffZ803NhptZlpll5ebmFhPewdNUeyISBSXWUqUVzBbSvkkdrn5+FrOXb4w6JJGD4u4bgWkE5Rirw/IOwvuCvx5zgMLFx+nAirA9vYj2PfYxsxTgEGB9MccqKraR7p7p7pmNGzcu2wmWkKbaE5EoKLGWKu+QGqmMvqwXjeukccXomSxf/0PUIYmUipk1NrN64eMawCnAfGASUDBLxzBgYvh4EjA0nOmjNcEgxRlhucgWM+sd1k9fvNc+Bcf6JfBeWIf9b2CAmdUPBy0OCNsipcGLIhIFJdYiQOM6aTxzSS/yd+zi0lEz2fTD9qhDEimN5sBUM/sKmElQY/0GcD9wqpktBE4Nn+Puc4BxwFzgHeDX7r4zPNY1wJMEAxoXA2+H7U8BDc1sEXAT4Qwj7r4e+N/wfWcCd4VtkVmfl8/4rOUMyTxMU+2JSLmyoMOhYsvMzPSsrKyow5BK4LNv13HRU9Pp0bI+z17ei7SU5APvJHKQzGyWu2dGHUd5iud1+4n3F3Pf2/MZMfAIrjqxbVzeQ0SqruKu2eqxFimkd5uG/OWX3Zi+ZD13TpobdTgiUgYqAxGRqKREHYBIojn7qBYsWL2Fx6ctplv6IQzt1TLqkERERKQCUGItUoSbB3Tgm+82cfvEOWzP387mLVvp0aYxfTpo6i6RRFcwhzWgUhARKVcqBREpQnKS8ejQozikRgp3vTmfV2Yu5b5Xv+DTBZq6SyTRqRRERKKixFpkP+rXqsbAwxuyC8jfBdu27+Tzb+O7qIWIiIhUXEqsRYoxsFsLjqybRM1kSEtNpkeb+C5qISIHT8uZi0hUVGMtUow+HZpy2zk9+PzbXNVYi1QQBSUgKgURkfKmxFrkAPp0aKqEWqSC0OIwIhKlyEpBzOwwM5tqZvPMbI6ZXR+2NzCzyWa2MLyvH1WMIiJSsagMRESiFGWP9Q7gd+7+uZnVAWaZ2WTgEmCKu99vZrcSLJt7S4RxiohIBaEyEBGJUmQ91u6+0t0/Dx9vAeYBLYBBwOhws9HA2ZEEKCIiIiJSCgkxK4iZZQBHAdOBpu6+EoLkG2gSYWgiIuXCzJLM7Juo46joVAoiIlGKfPCimdUGXgZucPfNZlbS/YYDwwFattSS0yJSsbn7LjP70sxauvuyqOOpqFQKIiJRijSxNrNUgqT6BXd/JWxebWbN3X2lmTUH1hS1r7uPBEYCZGZmerkELCISX82BOWY2A8graHT3s6ILSURESiqyxNqCrumngHnu/rdCL00ChgH3h/cTIwhPRCQKf4o6gIquoBQE4KoT20YcjYhUNVH2WPcFLgK+NrPZYdv/ECTU48zscmAZMCSa8EREype7v29mTYGeYdMMdy/yv3ZSNJWCiEiUIkus3f0jYH8F1f3LMxYRkURgZucCfwGmEVwf/25mv3f3CZEGVkFocRgRiVrkgxdFRGS324CeBb3UZtYYeBdQYl0CKgMRkagpsRYRSRxJe5V+rCNBpkWtCFQGIiJRU2ItIpI43jGzfwMvhc/PA96KMJ4KpUGtauqpFpFIqSdERCQBhDMlPQo8AXQFugEj3f2WSAOrQNbn5fPE+4tZn5cfdSgiUkWpx1pEJAG4u5vZa+5+NPDKAXeQfajGWkSipsRaRCRxfGZmPd19ZtSBVESqsRaRqCmxFhFJHCcBV5nZUoKVF42gM7trtGFVDKqxFpGoKbEWEUkAYY311cDSqGOpqDSPtYhETYm1iEgCCGusHwprrKUMVGMtIlHbb2JtZo+WYP/N7v7/YhiPiEhVphrrg6AaaxGJWnE91oOA2w+w/62AEmsRkdg4CbjazLJRjXWpqAxERBJBcYn1Q+4+uridzax+jOMREanKBkYdQEWlMhARSQT7Tazd/eED7VySbUREpGTcfamZHQe0d/dnzKwxUDvquCoClYGISCI44MqLZjbazOoVel7fzJ6Oa1QiIlWQmd0B3AKMCJtSgeeji6jiKJhqT2UgIhKlkixp3tXdNxY8cfcNwFFxi0hEpOoaDJxFUF+Nu68A6kQaUQWh5cxFJBGUJLFOKlxLbWYN0DR9IiLxkO/uDjiAmdWKOJ4Ko6DGenzW8qhDEZEqrCQJ8l+BT8xsQvh8CHBP/EISEamyxpnZE0A9M7sSuAz4V8QxVQiqsRaRRHDAxNrdnzWzLODksOkcd58b37BERKoed3/QzE4FNgMdgNvdfXLEYVUIWs5cRBJBSUs6qgPJwK7wsYiIxEGYSCuZLiXNYy0iiaAks4LcDowG6gONgGfMTIvCiIhIwlCNtYgkgpL0WJ8PHOXuPwKY2f3A58Dd8QxMRESkpFRjLSKJoCSzgmSzZ/lHGrA4LtEUYmanm9kCM1tkZrfG+/1EihJM0LD/5yKxZmY1zKxD1HFUNJrHWkQSQUkS623AHDMbZWbPAN8A35vZo2b2aDyCMrNk4DGC5X07AuebWcd4vJfI/mRnZ7N48eLdybS7s3jxYrKzs6MNTCotM/s5MBt4J3ze3cwmRRpUBaA5rEUkUZSkFOTV8FZgWnxC2UMvYJG7fwtgZmOAQYBmI5Fy4e7s2LGDnJwcANq2bcvixYvJyckhPT0dd8fMIo5SKqE7Ca5/0wDcfbaZZUQYT4VQUF8NaGYQEYlUSRLrb9x9VuEGM/u5u78ep5gAWgCFR6DkAMfE8f1E9mBmtG0b/ILOycnZnWCnp6fTtm1bJdUSLzvcfZM+X6Wj+moRSRQlKQX5l5l1KXhiZucD8Z4VpKjfKnsUt5rZcDPLMrOs3NzcOIcjVVHh5LqAkmqJs2/M7AIg2czam9nfgU+iDirRqb5aRBJFSRLrXwKjzezIcCWwa4EB8Q2LHKBw10M6sKLwBu4+0t0z3T2zcePGcQ5HqqKCmurCCtdci8TBb4BOBGNbXgQ2ATdEGZCIiJRcSVZe/NbMhgKvEZRnDHD3rXGOaybQ3sxaA98BQ4EL4vyeIrsVJNUFNdWFa6xBPdcSNx3c/TbgtqgDqUi0OIyIJIr9JtZm9jV7ll80IFh9cbqZ4e5d4xWUu+8ws+uAf4fv+bS7z4nX+4nszcxISUnZo6a6oCwkJSVFSbXEy9/MrDkwHhij617JaPCiiCSK4nqszyy3KIrg7m8Bb0UZg1RtGRkZe8z+UZBcK6mWeHH3k8ysGXAuMNLM6gJj3V0LchVDgxdFJFHsN7F296XlGYhIIto7iVZSLfHm7quAR81sKvAH4Ha00m2xCgYviohEbb+DF83s8wPtXJJtRESkZMJB4nea2TfAPwhmBEkvwX6HmdlUM5tnZnPM7PqwvYGZTTazheF9/UL7jAhXtl1gZqcVaj/azL4OX3vUwr8mzSzNzMaG7dMLz69tZsPC91hoZsNi9xUpGS0QIyKJorhSkCPN7KtiXjfgkBjHIyJSlT0DvEQwSHzFgTYuZAfwO3f/3MzqALPMbDJwCTDF3e83s1uBW4FbwpVshxLMQHIo8K6ZHe7uO4HHgeHAZwTleKcDbwOXAxvcvV04oP0B4DwzawDcAWQSjMuZZWaT3H3DwX0pSk411iKSKIpLrI8owf47YxWIiEhV5+69y7jfSmBl+HiLmc0jWGhrENAv3Gw0wYqOt4TtY9x9G7DEzBYBvcwsG6jr7p8CmNmzwNkEifUggpUhASYA/wh7s08DJrv7+nCfyQTJ+EtlOZeyUI21iCQK1ViLiESvDRQ5G5MBXppZmMISjaOA6UDTMOnG3VeaWZNwsxYEPdIFcsK27eHjvdsL9lkeHmuHmW0CGlL0SrktKEeqsRaRRFGSJc1FRCS+ChLTg5qNycxqAy8DN7j75mIG2+5vddviVr0tyz57xzecoMyEli1b7i82EZEKqyQrL4qISHxtD++vdfelhW8Eq90ekJmlEiTVL7j7K2Hz6nBebML7NWH7/la3zWHPwZKFV73dvY+ZpRCMsVlfzLH2EY8VczVwUUQSyQETazO7rvBIchERiZtTi2gbeKCdwlrnp4B57v63Qi9NAgpm6RgGTCzUPjSc6aM10B6YEZaNbDGz3uExL95rn4Jj/RJ4z92dYCGvAWZWP/xdMSBsKxcFAxfHZy0/8MYiInFWklKQZsDMcGq9p4F/hxdTERGJjcZhfXWbvWZjqgN8XIL9+wIXAV+b2eyw7X+A+4FxZnY5sAwYAuDuc8xsHDCXYEaRX4czggBcA4wCahAMWnw7bH8KeC4c6LieYFYR3H29mf0vMDPc7q6CgYzlQQMXRSSRWEly5LDnYgBwKcGUSuOAp9x9cXzDK5nMzEzPysqKOgwRkTIJk+HBwH0EU+IV2FKeSWp50nVbRCoqM5vl7plFvVaiGuuwh3pVeNsB1AcmmNmfYxaliEjVtdPds939/LCueivBAMDaZqZRfiIiFURJaqx/a2azgD8T/Euyi7tfAxwN/CLO8YmIVBlm9nMzWwgsAd4HsvmpFEOKoMGLIpJISlJj3Qg4Z+95rd19l5kd1NRQIiKyh7uB3sC77n6UmZ0EnB9xTAlNqy6KSCI5YGLt7rcX89q82IYjIlKlbXf3dWaWZGZJ7j7VzB6IOqhEpsGLIpJItECMiEji2Bgu8vIB8IKZrSEY1yL7oVUXRSSRaIEYEZHEMYhg4OKNwDvAYuDnkUYkIiIlph5rEZEE4e55hZ6OjiyQCmR9Xj7js5YzJPMwGtSqFnU4IlLFKbEWEYneUWa2udBzI5huzwhmPK0bTViJT4MXRSSRKLEWEYneF/tbbECKp8GLIpJIVGMtIpJAzOw4M7s0fNzIzFpHHVMiKxi8qDIQEUkESqxFRBKEmd0B3AKMCJuqAc9HF5GIiJRGJIm1mf3FzOab2Vdm9qqZ1Sv02ggzW2RmC8zstCjiExGJyGDgLCAPwN1XAHUijSiBadVFEUk0UfVYTwY6u3tX4L+EvTNm1hEYCnQCTgf+aWbJEcUoIlLe8t3dCQYuYma1Io4noRUMXByftTzqUEREgIgGL7r7fwo9/Qz4Zfh4EDDG3bcBS8xsEdAL+LScQxQRicI4M3sCqGdmVwKXAf+KOKaEpYGLIpJoEmFWkMuAseHjFgSJdoGcsE1EpFIzMyO4Fh4BbAY6ALe7++RIA0tgWnVRRBJN3BJrM3sXaFbES7e5+8Rwm9sIlut9oWC3Irb3/Rx/ODAcoGXLlgcdr4hIlNzdzew1dz+aoFxOREQqmLgl1u5+SnGvm9kw4Eygf1hTCEEPdeH/6aUDK/Zz/JHASIDMzMwik28RkQrmMzPr6e4zow5ERERKL6pZQU4nmFLqLHf/odBLk4ChZpYWzt3aHpgRRYwiIhE4CfjUzBaHsyZ9bWZfRR1UotKsICKSaKKqsf4HkAZMDsoK+czdr3b3OWY2DphLUCLya3ffGVGMIiLlbWDUAVQkWs5cRBJNVLOCtCvmtXuAe8oxHBGRhODuS6OOoSLRrCAikmgSYVYQERGRUtOsICKSaJRYixzInJHww1fQfACknxV1NCIiIpKgolp5UaRi+PR/YfzN8PkT8PH5kDMp6ohEJKTBiyKSaJRYixTHVkD1nbCqBvzwI6z8z4H3EZFyoSXNRSTRqBREpDjpAyH9WVjisLIWnHhc1BGJSEiDF0Uk0ajHWqQ46WdBv5eg1+mwPQWmT4Rdu6KOSkT4afBig1rVog5FRARQYi1yYOlnwWnj4fT7YMFbMOVPUUckIiIiCUilICIl1Ws4rP0vfPww1G8FmZdFHZGIiIgkEPVYi5SUGZz+ALQ/Dd68GRZOjjoikSpLM4KISCJSYi1SGskp8MunoWknGDcMVsyOOiKRKkkzgohIIlJiLVJaabXhgnFQswE8/wvI/W/UEYlUOUMyD2PEwCM0I4iIJBQl1iJlUbc5XDwRLAmeHQQblkYdkUiVohlBRCQRKbEWKauGbeGiV2H7D/DsWbB5ZdQRiYiISISUWIscjGad4cKXIW8tPHd2cC8iIiJVkhJrkYOVngnnj4EN2TDqDNiyOuqIRCo9zQoiIolIibVILLQ+PhjQuHFZkFxvXhF1RCKVmmYFEZFEpMRaJFbanBiUhWxZCc/8DDbqF75IvGhWEBFJREqsRWKp1bFw0Wvww3oY9TNYvyTqiEQqJc0KIiKJSIm1SKwd1hOGTYRtW+Dp02DlV1FHJCIiIuVAibVIPBx6FFz6DiSlBGUh374fdUQiIiISZ0qsReKlyRFw+WQ4JB1e+CV880rUEYlUGpoVREQSkRJrkXg6pAVc9ja0OBomXAbTn4g6IpFKQbOCiEgiijSxNrObzczNrFGhthFmtsjMFpjZaVHGJxITNeoHKzR2+Bm8/Qd46w+wc0fUUYlUaJoVREQSUWSJtZkdBpwKLCvU1hEYCnQCTgf+aWbJ0UQoEkOpNeC856DPdTDjCXjpPPhxU9RRiVRYmhVERBJRlD3WDwF/ALxQ2yBgjLtvc/clwCKgVxTBicRcUjKcdg/8/BH4dho8NSBYrVFEREQqhUgSazM7C/jO3b/c66UWQOGCuZywrahjDDezLDPLys3NjVOkInFw9CVw4SuwZRX862RY+mnUEYmIiEgMxC2xNrN3zeybIm6DgNuA24varYg2L6INdx/p7pnuntm4ceNYhi4Sf21OhCumQPV6MPrnMONf4EV+1EVERKSCiFti7e6nuHvnvW/At0Br4EszywbSgc/NrBlBD3XhkSjpwIp4xSgSqUbt4Mr3oO3J8NbN8No1sH1r1FGJJDxNtSciiarcS0Hc/Wt3b+LuGe6eQZBM93D3VcAkYKiZpZlZa6A9MKO8YxQpNzXqwfljoN//wJdj4KlTVXctpWZmT5vZGjP7plBbAzObbGYLw/v6hV4rcvYlMzvazL4OX3vUzCxsTzOzsWH7dDPLKLTPsPA9FprZsPI4X021JyKJKqHmsXb3OcA4YC7wDvBrd98ZbVQicZaUBP1ugQvGwcZl8MSJsPDdqKOSimUUwUxKhd0KTHH39sCU8PmBZl96HBhO0KnRvtAxLwc2uHs7goHnD4THagDcARxDMND8jsIJfLxoqj0RSVSRJ9Zhz/XaQs/vcfe27t7B3d+OMjaRcnX4ABg+7aeVGqfcpfmupUTc/QNg/V7Ng4DR4ePRwNmF2veZfcnMmgN13f1Td3fg2b32KTjWBKB/2Jt9GjDZ3de7+wZgMvsm+DGnqfZEJFFFnliLSCEN2gTLoB91IXz4Vxh1BmzUv7ulTJq6+0qA8L5J2L6/2ZdahI/3bt9jH3ffAWwCGhZzLBGRKkmJtUiiqVYTBv0DfvEUrJ4D/3cczH8z6qik8tjf7EvFzcpUln32fWNNkyoilZwSa5FE1eWXcNX7UD8DxlwQLIW+Y1vUUUnFsTos7yC8XxO272/2pZzw8d7te+xjZinAIQSlJ6WayUnTpIpIZafEWiSRNWwLl/8Hel8bLIX+ZH9YMy/qqKRimAQUzNIxDJhYqH2f2ZfCcpEtZtY7rJ++eK99Co71S+C9sA7738AAM6sfDlocELbFlabbE5FEpcRaJNGlpMHp98H5Y2HzymDWkE8fg127oo5MEoSZvQR8CnQwsxwzuxy4HzjVzBYCp4bPDzT70jXAkwQDGhcDBQPInwIamtki4CbCGUbcfT3wv8DM8HZX2BZXmm5PRBKVeSVY7S0zM9OzsrKiDkMk/r7Phdd/Cwvegozj4ezHoZ6mHKvozGyWu2dGHUd5Opjr9vq8fMZnLWdI5mGaGUREyl1x12z1WItUJLUbw9AX4ax/wIov4PFjg4VlKsEfyCIlpen2RCRRKbEWqWjMoMdFcPVH0LQTvHoVjB8GeWsPvK+IiIjEjRJrkYqqQWu45E045U6Y/xY81gu+nqDeaxERkYgosRapyJKS4bgb4eoPg2n5Xr4cXjofNu93xjMRERGJEyXWIpVBkyODFRsH3APfToPHesOs0eq9lkpHU+2JSCJTYi1SWSQlw7HXwTUfQ7Muwewhzw6CDdlRRyYSM5pqT0QSmRJrkcqmYVsY9jqc8Tf47nP4Zx/4+BHYuT3qyEQO2pDMwxgx8AiGZGqaSRFJPEqsRSqjpCToeTn8+jNofSJMvh2eOAGWfRZ1ZCIHRVPtiUgiU2ItUpkdkg4XjAnmvv5xMzx9Gkz6DfwQ98XxREREqhwl1iJVwRFnwK+nw7G/hS9egH9kBvca3CgiIhIzSqxFqoq02jDgf+GqD6BhO5h4LTzzM1g9N+rIREREKoWUqAMQqSq2b99OTk4OP/74Y8SRJEOfhyEzD7ZuhG9z4LuNUP0QMP2tHU/Vq1cnPT2d1NTUqEOpsNbn5TM+azlDMg9TnbWIJBwl1iLlJCcnhzp16pCRkYGZRR1OYOcO2LISflgLthPqNoWaDYNl0yWm3J1169aRk5ND69atow6nwiqYbg/gqhPbRhyNiMielFiLlJMff/wxsZJqgOQUqHcY1GoIm76DTcshb20w6DGtdtTRVSpmRsOGDcnNzY06lAqtYJo9TbcnIolIibVIOUqopLqw1JpB3fWPG4MEe91CqF4f6h4KKfp3e6wk7Pe/AimYbk9EJBFFVlBpZr8xswVmNsfM/lyofYSZLQpfOy2q+ESqHDOoUZ/5640+Z19JWvPDefCuW4NSkV07D+rQ8+fPp0+fPqSlpfHggw/GJNwPPviAHj16kJKSwoQJE3a3z549mz59+tCpUye6du3K2LFji9z/kksuoWbNmmzZsmV32/XXX4+ZsXbt2pjEKCIiVUskibWZnQQMArq6eyfgwbC9IzAU6AScDvzTzJKjiFGkqmrQqDGPPvY4N990I6RUhy2rYM28oESkjNPzNWjQgEcffZSbb7651PtmZ2fTr1+/fdpbtmzJqFGjuOCCC/Zor1mzJs8++yxz5szhnXfe4YYbbmDjxo1FHrtdu3ZMnDgRgF27djF16lRatGhR6hhFREQguh7ra4D73X0bgLuvCdsHAWPcfZu7LwEWAb0iilGkUsnOzuaII47giiuuoHPnzvzqV7/i3XffpW/fvrRv354ZM2YA0KRJE3r27ElqWg2o2QAatofkakH9de78YKGZUibYu49ZxGwYzz//PL169aJ79+5cddVV7NxZst7xjIwMunbtSlLSnpexww8/nPbt2wNw6KGH0qRJk/3WNZ9//vm7e7SnTZtG3759SUkJKuTy8vI444wz6NatG507d969XUZGxu4e7aysrN1J/5133slll11Gv379aNOmDY8++miJzkNERCqPqGqsDweON7N7gB+Bm919JtACKLzmck7YJlKp/On1OcxdsTmmx+x4aF3u+HmnYrdZtGgR48ePZ+TIkfTs2ZMXX3yRjz76iEmTJnHvvffy2muv7btTWm1o1B5+3ASbv4P1i6Fabc67ZgQLFi7aZ/ObbrqJiy++uEQxz5s3j7Fjx/Lxxx+TmprKtddeywsvvFDi/Q9kxowZ5Ofn07Zt0TW57du3Z+LEiWzYsIGXXnqJCy+8kLfffhuAd955h0MPPZQ333wTgE2bNh3w/ebPn8/UqVPZsmULHTp04JprrtHUejGm6fZEJJHFLbE2s3eBZkW8dFv4vvWB3kBPYJyZtQGKGtlTZNeYmQ0HhkPwL2ERObDWrVvTpUsXADp16kT//v0xM7p06UJ2dvb+dzSDGvWget2gJGTLKsY+chvUaAB1mpd5gOOUKVOYNWsWPXv2BGDr1q00adIEgMGDB7NkyRLy8/NZtmwZ3bt3B4I66EsvvfSAx165ciUXXXQRo0eP3qdXu7BzzjmHMWPGMH36dJ544ond7V26dOHmm2/mlltu4cwzz+T4448/4HueccYZpKWlkZaWRpMmTVi9ejXp6ekH3E9KTtPtiUgii1ti7e6n7O81M7sGeMXdHZhhZruARgQ91IXnUEoHVuzn+COBkQCZmZlal1kqlAP1LMdLWlra7sdJSUm7nyclJbFjx44DH8CSoHYTqNmA8355DgsWLAgPlhJM3YeVqsfa3Rk2bBj33XffPq+9+uqrQFDCcskllzBt2rQSHRNg8+bNnHHGGdx999307t272G2HDh1Kjx49GDZs2B4J+OGHH86sWbN46623GDFiBAMGDOD2228nJSWFXbt2Aeyz2E/hr29ycnLJvqZSKppuT0QSWVSlIK8BJwPTzOxwoBqwFpgEvGhmfwMOBdoDMyKKUUT2JymFsa9Mgh3bgsGNW9f/lHTXalLiw/Tv359BgwZx44030qRJE9avX8+WLVto1apVmUPLz89n8ODBXHzxxQwZMuSA27ds2ZJ77rmHU07Zsy9gxYoVNGjQgAsvvJDatWszatQoIKixnjVrFgMHDuTll18uc5xSNppuT0QSWVSJ9dPA02b2DZAPDAt7r+eY2ThgLrAD+LW7H9w8XyJSKqtWrSIzM5PNmzeTlJTEww8/zNy5c6lbt+6+G6ekQf1WQUK9ZWWQZOethdpNoWYjCHuA93fMjh07cvfddzNgwAB27dpFamoqjz32WIkS65kzZzJ48GA2bNjA66+/zh133MGcOXMYN24cH3zwAevWrdudDI8aNWp3KUlRrrrqqn3avv76a37/+9+TlJREamoqjz/+OAB33HEHl19+Offeey/HHHPMgb+gIiJSZZiXcfqsRJKZmelZWVlRhyFSrHnz5nHkkUdGHUZ85efB5pWQvwWSUqFOs2BmEYtsyvyEU9TnwMxmuXtmRCFFQtdtEamoirtm67ediMROtVrQqF2wimNyajBF35r58MP6Ms+BLSIiUlEosRaR2EurA40OhwZtghlFNi4N5sDeukEJtoiIVFpKrEUkPsyg+iHQ+Aio1wpw2JCtBFvKbH1ePk+8v5j1eflRhyIiUqSoBi+KSFVhFtRZ16gfJNTfrwoS7JTqQQ129XrBNiIHoDmsRSTRKbEWkfKxT4K9Wgm2lIrmsBaRRKfEWkTKlxJsKSPNYS0iiU411iKyhxdeeIGuXbvStWtXjj32WL788suDPub8+fPp06cPaWlpPPjgg0FjQYK9uwabIMFeMw9+WAe+q9hjfvDBB/To0YOUlBQmTJiwu3327Nn06dOHTp060bVrV8aOHVvk/pdccgk1a9Zky5Ytu9uuv/56zIy1a9ce1PmKiEjVpMRaRPbQunVr3n//fb766iv++Mc/Mnz48IM+ZoMGDXj00Ue5+eab932xcIJdPyOY83rjsiDB/j6X7G8X069fv312a9myJaNGjeKCCy7Yo71mzZo8++yzzJkzh3feeYcbbriBjRs3FhlXu3btmDhxIgC7du1i6tSptGjR4iDPVkREqiol1iJVRHZ2NkcccQRXXHEFnTt35le/+hXvvvsuffv2pX379syYMQOAY489lvr16wPQu3dvcnJyDvq9mzRpQs+ePUlNTd3nteeff55evXrR/aijuOqGW9nZoF0wTV9SKmzOgXULYWc+7Nqxx34ZGRl07dqVpKQ9L2OHH3447du3B+DQQw+lSZMm5ObmFhnX+eefv7tHe9q0afTt25eUlKBCLi8vjzPOOINu3brRuXPn3dtlZGTs7tHOysranfTfeeedXHbZZfTr1482bdrw6KOPlvGrJSIiFZVqrEWi8PatsOrr2B6zWRcYeH+xmyxatIjx48czcuRIevbsyYsvvshHH33EpEmTuPfee3nttdf22P6pp55i4MCBRR7rvPPOY8GCBfu033TTTVx88cUlCnnevHmMHTuWjz/+mNTUVK699lpeePHFYP/qh8C272H1piCxXj0XajWCWo2DxWdKYMaMGeTn59O2bdF1ue3bt2fixIls2LCBl156iQsvvJC3334bgHfeeYdDDz2UN998E4BNmzYd8P3mz5/P1KlT2bJlCx06dOCaa64p8o8JERGpnJRYi1QhrVu3pkuXLgB06tSJ/v37Y2Z06dKF7OzsPbadOnUqTz31FB999FGRx9pf7XJpTJkyhVmzZtGzZ08Atm7dSpMmTQAYPHgwS5YsIT8/n2XLltF9wFDYtZPrr7iASy+7DGo1gZS0/R575cqVXHTRRYwePXqfXu3CzjnnHMaMGcP06dN54okndrd36dKFm2++mVtuuYUzzzyT448//oDnc8YZZ5CWlkZaWhpNmjRh9erVpKenl/TLIQewPi+f8VnLGZJ5GA1qVYs6HBGRfSixFonCAXqW4yUt7adENCkpaffzpKQkduz4qdTiq6++4oorruDtt9+mYcOGRR4rFj3W7s6wYcO477779nnt1VdfBYISlksuuYRp06bB9h8hbzXkrYO8tcHMInuViABs3ryZM844g7vvvpvevXsXG8PQoUPp0aMHw4YN2yMBP/zww5k1axZvvfUWI0aMYMCAAdx+++2kpKSwa1cwsPLHH3/c41iFv77Jycl7fE3l4GkeaxFJdEqsRWQPy5Yt45xzzuG5557j8MMP3+92seix7t+/P4MGDeLGG2+kSZMmrF+/ni1bttCqVauid0itHswgUrs55K0JZg/5cTNsXhncp9Uhf/t2Bg8ezMUXX8yQIUMOGEPLli255557OOWUU/ZoX7FiBQ0aNODCCy+kdu3ajBo1CghqrGfNmsXAgQN5+eWXD/ZLIKWgeaxFJNEpsRaRPdx1112sW7eOa6+9FoCUlBSysrIO6pirVq0iMzOTzZs3k5SUxMMPP8zcuXPp2LEjd999NwMGDGDXrl2kpqby2GOP7T+xLpBSjZn/XcngwYPZsGE9r7/7AXc88DBzPnyDca9P44MPPmDdunW7k+FRo0bRvXv3/R7uqquu2qft66+/5ve//z1JSUmkpqby+OOPA3DHHXdw+eWXc++993LMMceU8SsiZaF5rEUk0Zm7Rx3DQcvMzPSD/cUvEm/z5s3jyCOPjDqMysl3hYvNrIEdPwYzitRqFNySEqv/oKjPgZnNcvfMiEKKhK7bIlJRFXfNTqzfOCIiZWFJULMh1GgA27YECfaWlcGqjjUbBjOJFDPQUUREJBaUWItI5WEG1esGt+0/wPe5wSDHvNxgqfRajaFaLS2ZLiIicaHEWkQqp9SaUL8V1GkOP+QGM4n8uBFSawQJdo36QU+3iIhIjOi3iohUbinVoG4LaNoJDjkM3IMl01fPCWYT2bk96gilhNbn5fPE+4tZn5cfdSgiIkVSj7WIVA1JycFgxpoNgzrsvFz4flVQh12j3k9lIpKwNI+1iCQ6JdYiUrUUrsPe8WNQg/3DumBWkdRaULtxsJy6ykQSjuaxFpFEF8lvDjPrbmafmdlsM8sys16FXhthZovMbIGZnRZFfCKV1apVqxg6dCht27alY8eO/OxnP+O///1vkdtmZ2djZvzxj3/c3bZ27VpSU1O57rrryivk+EqpDoekQ9POQbnIru2wIRtWz4Utq1QmkmAK5rHWcuYikqii6pL5M/And+8O3B4+x8w6AkOBTsDpwD/NLDmiGEUitfcc8wc757y7M3jwYPr168fixYuZO3cu9957L6tXr97vPm3atOGNN97Y/Xz8+PF06tTpoOJISEnJULsJNOkIDdoECfeWlUEd9vps2PZ9UJstIiJSjKgSawfqho8PAVaEjwcBY9x9m7svARYBvYrYX6RSy87OZvHixbuTaXdn8eLFZGdnl/mYU6dOJTU1lauvvnp3W/fu3TnuuOP4/e9/T+fOnenSpcseS5XXqFGDI488cvfKi2PHjuXcc88tcwwJzywoA2nUDhofGdRkb9sM6xZC7vygLnvXzqijFBGRBBVVjfUNwL/N7EGC5P7YsL0F8Fmh7XLCNpEqw93ZsWMHOTk5ALRt25bFixeTk5NDeno67o6VYR7mb775hqOPPnqf9ldeeYXZs2fz5ZdfsnbtWnr27MkJJ5yw+/WhQ4cyZswYmjVrRnJyMoceeigrVqzY5ziVTmpYJlKneVB/nbcWNuXA5hXBQjS1GgVT9wkAZnY68AiQDDzp7vdHHJKISLmLW2JtZu8CzYp46TagP3Cju79sZucCTwGnAEVlC0X+/9XMhgPDAVq2bBmTmEUSgZnRtm0w40FOTs7uBDs9PZ22bduWKakuzkcffcT5559PcnIyTZs25cQTT2TmzJl07doVgNNPP50//vGPNG3alPPOOy+m710hFJ5NZPsPPw12/GFtMItIzUbBrCJVeLBjWLL3GHAqQYfITDOb5O5zo41MRKR8xe03gbuf4u6di7hNBIYBr4Sbjuenco8coPBw73R+KhPZ+/gj3T3T3TMbN24cr9MQiUTh5LrAwSbVnTp1YtasWfu0H6h2u1q1ahx99NH89a9/5Re/+EWZ37/CMwsS6fqtwsGOhwaDGzcuDefEXgE7quz8yr2ARe7+rbvnA2MISvtiRnNYi0hFEFUXywrgxPDxycDC8PEkYKiZpZlZa6A9MCOC+EQiVVBTXVjhmuuyOPnkk9m2bRv/+te/drfNnDmT+vXrM3bsWHbu3Elubi4ffPABvXrtObThd7/7HQ888AANGzYs8/tXKskpULtpONixbbDK4/erYUsVKJEpWgtgeaHnRZbxmdnwcCaorNzc3FK9QcEc1uOzlh94YxGRiERVY30l8IiZpQA/EpZ0uPscMxsHzAV2AL92d40UkiqlIKkuqKkuXGMNZe+5NjNeffVVbrjhBu6//36qV69ORkYGDz/8MN9//z3dunXDzPjzn/9Ms2bN9hgo2alTp8o5G8jB2mNO7G1RRxOlEpXxuftIYCRAZmZmqf5K1BzWIlIR2MFO4ZUIMjMzvWDWApFENW/ePI488sgSbZudnc2OHTt2J9EFyXZKSgoZGRnxDVTiqqjPgZnNcvfMiEI6aGbWB7jT3U8Ln48AcPf79rePrtsiUlEVd83WyosiCSgjI2OP2T8Kaq5jPXBRJEZmAu3DEr7vCNYjuCDakEREyp8Sa5EEtXcSraRaEpW77zCz64B/E0y397S7z4k4LBGRcqfEWkREDpq7vwW8FXUcIiJRqroTr4pEoDKMaZCy0/dfRKRyU2ItUk6qV6/OunXrlFxVUe7OunXrqF69etShiIhInKgURKScpKenk5OTQ2nn75XKo3r16qSnp0cdhoiIxIkSa5FykpqaSuvWraMOQ0REROJEpSAiIiIiIjGgxFpEREREJAaUWIuIiIiIxEClWNLczHKBpVHHUYRGwNqog4izqnCOUDXOU+cYnVbu3jjqIMpTGa/bifT9S5RYEiUOUCxFSZQ4IHFiSZQ4oOyx7PeaXSkS60RlZln7W0u+sqgK5whV4zx1jpLoEun7lyixJEocoFgSOQ5InFgSJQ6ITywqBRERERERiQEl1iIiIiIiMaDEOr5GRh1AOagK5whV4zx1jpLoEun7lyixJEocoFiKkihxQOLEkihxQBxiUY21iIiIiEgMqMdaRERERCQGlFjHiJkNMbM5ZrbLzDL3em2EmS0yswVmdlqh9qPN7OvwtUfNzMo/8rIxszvN7Dszmx3eflbotSLPtyIys9PD81hkZrdGHU+smFl2+NmbbWZZYVsDM5tsZgvD+/pRx1kaZva0ma0xs28Kte33nCrT57QqKO+fxUT5PJnZYWY21czmhb9jro8iFjOrbmYzzOzLMI4/RRFHoWMnm9kXZvZGxHGU6loa51jqmdkEM5sffl76RPA56VAoL5htZpvN7IYIvyY3hp/Xb8zspfBzHN9Y3F23GNyAI4EOwDQgs1B7R+BLIA1oDSwGksPXZgB9AAPeBgZGfR6lON87gZuLaN/v+Va0G5Acxt8GqBaeV8eo44rRuWUDjfZq+zNwa/j4VuCBqOMs5TmdAPQAvjnQOVWmz2lVuEXxs5gonyegOdAjfFwH+G/4fuUaS/h7qnb4OBWYDvSO6mcMuAl4EXgjqu9NePwSX0vLIZbRwBXh42pAvSivgeHP7SqgVUQ/Oy2AJUCN8Pk44JJ4x6Ie6xhx93nuvqCIlwYBY9x9m7svARYBvcysOVDX3T/14Dv6LHB2+UUcN0Web8QxlVUvYJG7f+vu+cAYgvOrrAYRXJgJ78+OLpTSc/cPgPV7Ne/vnCrT57QqKPefxUT5PLn7Snf/PHy8BZhHkDCUaywe+D58mhrevLzjADCzdOAM4MlCzYn0sx7F16QuwR+DTwG4e767b4wilkL6A4vdfWmEcaQANcwsBagJrIh3LEqs468FsLzQ85ywrUX4eO/2iuQ6M/sq/Jdpwb9S9ne+FVFlOpe9OfAfM5tlZsPDtqbuvhKCX+ZAk8iii539nVNl/t5WRony/Yr082RmGcBRBL3F5R5LWH4xG1gDTHb3SOIAHgb+AOwq1BbV96Y019J4xtIGyAWeCUtknjSzWhHFUmAo8FL4uNzjcPfvgAeBZcBKYJO7/yfesSixLgUzezes09n7VlzPSVF1015Me8I4wPk+DrQFuhN8YP9asFsRh0qo8yqFynQue+vr7j2AgcCvzeyEqAMqZ5X5e1sZJfr3K+7xmVlt4GXgBnffHEUs7r7T3bsD6QT/ee1c3nGY2ZnAGnefVdJd4hFHIaW5lsYzlhSC0qXH3f0oII+gzCGKWDCzasBZwPgDbRqvOMIOv0EEZR2HArXM7MJ4x5JS2h2qMnc/pQy75QCHFXqeTvCviJzw8d7tCaOk52tm/wLeCJ/u73wrosp0Lntw9xXh/Roze5Xg312rzay5u68MS5XWRBpkbOzvnCrt97aSSpTvVySfJzNLJUiqX3D3V6KMBcDdN5rZNOD0COLoC5xlwYD56kBdM3s+gjiAUl9L4xlLDpAT/hcBYAJBYh3V52Qg8Lm7rw6fRxHHKcASd88FMLNXgGPjHYt6rONvEjDUzNLMrDXQHpgR/vthi5n1NjMDLgYmRhloaYQfxgKDgYKR80Web3nHFyMzgfZm1jr863sowflVaGZWy8zqFDwGBhB8/yYBw8LNhlGBPo/F2N85VabPaVWQKD+L5f55Cn8/PAXMc/e/RRWLmTU2s3rh4xoEScv88o7D3Ue4e7q7ZxB8Dt5z9wvLOw4o07U0brG4+ypguZl1CJv6A3OjiCV0Pj+VgRS8X3nHsQzobWY1w5+j/gRjFOIbS2lHO+q239Gngwn+2tkGrAb+Xei12whGly6g0MwfQCbBD+Fi4B+EC/ZUhBvwHPA18FX4YWx+oPOtiDfgZwSj8BcDt0UdT4zOqQ3ByOcvgTkF5wU0BKYAC8P7BlHHWsrzeomgLGl7+LN4eXHnVJk+p1XhVt4/i4nyeQKOI/h39FfA7PD2s/KOBegKfBHG8Q1we9ge2c8Y0I+fZgWJ4ntT6mtpPL8mBKWZWeH36DWgfkRfl5rAOuCQQm1RfU3+RPAH4DcEeUtavGPRyosiIiIiIjGgUhARERERkRhQYi0iIiIiEgNKrEVEREREYkCJtYiIiIhIDCixFhERERGJASXWUuWZWYaZbQ2X6i14/s0BdivJcW80s2Vm9o+DDlJEpBIws3pmdm2h54ea2YQ4vM+dZvadmd0VPr/EzNzM+hfaZnDY9stSHvsFM1tf2v2kalBiLRJY7MFSvTHj7g8Bt8fymCIiFVw9YHdi7e4r3D1eCepD7l74Gvw1wcIlBYYSzEFdYmaW7O6/ohIsFibxocRaKjUz62lmX5lZ9XCVrDlm1rkU+7cxsy/C41xiZhPN7B0zW2BmdxTa7uLwfb40s+ficzYiIhXe/UBbM5ttZn8p/B/C8Br7mpm9bmZLzOw6M7spvAZ/ZmYNwu3ahtfhWWb2oZkdUcL3/hDoZWapZlYbaEew4A5m1j9ckpzw+anhEtiY2fdmdpeZTQf6xOwrIZVSStQBiMSTu880s0nA3UAN4Hl3L1GZR7g07BjgUnefbWadgF5AZ+AHYKaZvQlsJVitqa+7ry24+IuIyD5uBToX/IfQzDL2er0zcBRQHVgE3OLuR5nZQ8DFwMPASOBqd19oZscA/wROLsF7O/AucBpwCEGvc+vwtfeAx8yssbvnApcCz4Sv1QK+2av3W6RISqylKrgLmAn8CPy2hPs0BiYCv3D3OYXaJ7v7OoCwN+M4YCcwwd3XArj7+lgFLiJSxUx19y3AFjPbBLwetn8NdA17mo8FxptZwT5ppTj+GILfA4cAvwP+B8DdPfxv44Vm9gxBz/TF4T47gZfLfkpSlSixlqqgAVAbSCXoBckrwT6bgOVAX6BwYu17beeAFdEuIiKlt63Q412Fnu8iyFmSgI1lHRPj7jPCcsCt7v7fQsk5BD3UrxN0wox39x1h+4/uvrMs7ydVj2qspSoYCfwReAF4oIT75ANnAxeb2QWF2k81swZmViN8/WNgCnCumTUEUCmIiMh+bQHqlHVnd98MLDGzIQAW6FbKw4wg7Kne69grgBXA/wNGlTVGqdrUYy2VmpldDOxw9xfNLBn4xMxOdvf3DrSvu+eZ2ZnAZDMr6OX+CHiOYNDLi+6eFb7PPcD7ZrYT+AK4JA6nIyJSobn7OjP7OByw+DbwWBkO8yvgcTP7fwT/iRxDKWb3cPe3i3n5BaCxu88tQ1wimLv+gy1VWzh45g13L3a2EDO7BMh09+tKcexS7yMiIgfHzO4Evnf3B0u53z+AL9z9qQNsN4rg90bM5+CWik2lICLBwJRDChaIiRUzu5HgX46bY3lcERE5oO+B4QULxJSEmc0CugLPH2C7F4ATCWqxRfagHmsRERERkRhQj7WIiIiISAwosRYRERERiQEl1iIiIiIiMaDEWkREREQkBpRYi4iIiIjEgBJrEREREZEY+P8z7vM/cMN0KAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(12,5))\n",
    "ax1 = fig.add_subplot(121, aspect='auto')\n",
    "# ax1.scatter(X[:,0], Y[:,0], s=1, label=\"m1\")\n",
    "# ax1.scatter(X[:,1], Y[:,1], s=1, label=\"m2\")\n",
    "plt.scatter(0,0, marker='x', color='grey', alpha=0.5, label=r'CoM')\n",
    "plt.scatter(X[0,0], Y[0,0], s=10, c='steelblue') # initial position\n",
    "plt.scatter(X[0,1], Y[0,1], s=10, c='orange') # initial position\n",
    "plt.plot(X[:,0], Y[:,0], label=f\"m1 = {m1:.3} Msun\")\n",
    "plt.plot(X[:,1], Y[:,1], label=f\"m2 = {m2:.3} Msun\")\n",
    "ax1.set_xlabel(\"x [kpc]\"), ax1.set_ylabel(\"y [kpc]\")\n",
    "plt.legend(loc='lower right')\n",
    "\n",
    "ax2 = fig.add_subplot(122,aspect='auto')\n",
    "ax2.scatter(t, err, s=1)\n",
    "#ax2.set_ylim(err_tol*np.array([-1,1]))\n",
    "ax2.set_xlabel(\"time [Myr]\"), ax2.set_ylabel(\"relative error\")\n",
    "print(f\"error at {t[-1]} Myr: {err[-1]:.6} kpc\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "UFRC Python-3.8",
   "language": "python",
   "name": "python3-3.8-ufrc"
  },
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}