{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import odeint\n",
    "from scipy.integrate import solve_ivp\n",
    "from matplotlib import animation\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib.animation import PillowWriter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Newton's law of gravity $F = ma = GMm/r^2$ and hence $a = GM/r^2$. For a collection of masses, written in terms of vector calculus:\n",
    "\n",
    "$$\\frac{d^2 \\vec{r}_i}{dt^2} = \\sum_{j \\neq i} \\frac{Gm_j}{r_{ij}^3} \\vec{r}_{ij} $$\n",
    "\n",
    "where $\\vec{r}_i$ and $m_i$ is the location and mass of the ith particle, and $\\vec{r}_{ij} = \\vec{r}_j - \\vec{r}_i$ (the tail of this vector is at $\\vec{r}_i$ and points to $\\vec{r}_j$).\n",
    "\n",
    "Dimensions for this problem are gross $G = 6.67 \\times 10^{-11}$ too small for efficient computation on computer. Convert to dimensionless quantities:\n",
    "\n",
    "$$\\frac{d^2 \\vec{r}_i'}{dt'^2} = \\sum_{j \\neq i} \\frac{m_j'}{r_{ij}'^3} \\vec{r}_{ij}' $$\n",
    "\n",
    "where\n",
    "\n",
    "* $\\vec{r}_i' = \\vec{r}_i / L$ and $\\vec{r}_{ij}' = \\vec{r}_{ij} / L$ where $L$ is some characteristic length of the problem\n",
    "* $m_i' = m_i / M$ where $M$ is some characteristic length of the problem\n",
    "* $t' = t\\sqrt{GM/L^3}$ where $M$ and $L$ are the characteristic length and mass.\n",
    "\n",
    "We will drop the primes for the rest of this code, assuming everything is dimensionless. We can solve for the motion of the object, as a function of dimensionless time $t'$ and then solve for $t$ using \n",
    "\n",
    "$$ t = t' \\sqrt{L^3 / GM}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2 Body Problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we will set \n",
    "\n",
    "* $M = \\text{mass of earth} = 5.97\\times 10^{24}~$kg \n",
    "* $L = \\text{earth-sun distance} = 1.5 \\times 10^{11}~$m\n",
    "\n",
    "and simulate earths motion around the sun"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "m1 = 1\n",
    "m2 = 333000 # sun 333000x heavier than earth\n",
    "x1_0 = 1 # initial position is one earth-sun distance away from the sun\n",
    "y1_0 = 0\n",
    "x2_0 = 0\n",
    "y2_0 = 0\n",
    "vx1_0 = 0\n",
    "vy1_0 = np.sqrt(m2) #circular motion v=sqrt(a*r) with a=m2 (since G=1, r12=1)\n",
    "vx2_0 = 0\n",
    "vy2_0 = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that any second order ODE\n",
    "\n",
    "$$\\frac{d^2 y}{dt^2} = F\\left(y, \\frac{dy}{dt}, t\\right) $$\n",
    "\n",
    "can be rewritten as two coupled first order ODEs when defining $v_y = dy/dt$\n",
    "\n",
    "$$ \\frac{dy}{dt} = v_y $$\n",
    "$$ \\frac{dv_y}{dt} = F\\left(y, v_y, t\\right) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we define our system of things we need to solve for as \n",
    "\n",
    "$$\\vec{S} = (x_1, y_1, x_2, y_2, v_{x1}, v_{y1}, v_{x2}, v_{y2})$$\n",
    "\n",
    "and we need to write in a function that takes in $S$ and returns $d\\vec{S}/dt$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dSdt(S, t):\n",
    "    x1, y1, x2, y2, vx1, vy1, vx2, vy2 = S\n",
    "    r12 = np.sqrt((x2-x1)**2 + (y2-y1)**2)\n",
    "    return [ vx1,\n",
    "            vy1,\n",
    "            vx2,\n",
    "            vy2,\n",
    "            m2/r12**3 * (x2-x1),\n",
    "            m2/r12**3 * (y2-y1),\n",
    "            m1/r12**3 * (x1-x2),\n",
    "            m1/r12**3 * (y1-y2) ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get array of times "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.linspace(0,1,10000) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol = odeint(dSdt, y0=[x1_0, y1_0, x2_0, y2_0, vx1_0, vy1_0, vx2_0, vy2_0],\n",
    "             t=t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot $x$ as a function of time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAD4CAYAAAD//dEpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA+gklEQVR4nO3dd3xc5ZXw8d8ZjXqzJavLRS6y5N5NsY0bxhiCgTQIm5CErJc3Yd8k2yBvdrPJZtnNJrvJbvpCQgKBYJwEMNgGG0uWwRT3qmLLXV2WrWL1Ms/7h0ZEyJItaUZzp5zv5zOfmblzy/H11Zx5nvsUMcaglFJK3YjN6gCUUkr5Bk0YSimlBkUThlJKqUHRhKGUUmpQNGEopZQaFLvVAQzHqFGjzOTJk60O44aampqIjIy0Oowb0jjdxxdiBI3T3XwlzoMHD9YYYxKGu71PJoykpCQOHDhgdRg3lJeXx/Lly60O44Y0TvfxhRhB43Q3X4lTRC64sr1WSSmllBoUTRhKKaUGRROGUkqpQdGEoZRSalA0YSillBoUtyQMEXlGRKpF5MQAn4uI/FhETovIMRGZ1+uztSJy0vnZE+6IRymllPu5q4TxW2DtdT6/E5jifGwAfgEgIkHAz5yfTwMeFJFpbopJKaWUG7mlH4Yx5m0RmXCdVdYDz5nusdQ/EJFRIpICTABOG2POAojIRue6Bdc7Xn2b4dXDZSyYMJr00RHu+CcoZYmaxjb2nr1CRX0LrR1dJMWEkZUcw4y0GETE6vCU+ghPddxLA0p6vS91Lutv+eL+diAiG+gunRCSPJmvvXQEgPExNu6cEMyilCBsXvYH1tjYSF5entVh3JDG6T6DidEYw/GaLt4410HhFUe/68SFCSvG2rl9fDBhdvdf175wLkHj9DaeShj9XfHmOsuvXWjMU8BTAFMyp5rXvrqUPcU1bDpQwi+PNfJOTQw//NQcpiZHuy9qF/lK70+N031uFGN5XQvfePk4u09dIjkmjK+vHseyzDFMHBNFWIiNyvpWDpyvZfPRcv506hJ55cKT981g7YwUj8bpLTRO7+KphFEKjO31Ph0oB0IGWH5dNoHslBiyU2J4ZEkGW45X8C+v5/Oxn+zhu/dO59MLx7k1eKXc4Z3iS3z5hUN0OQz//LFpPLR4PCH2j95GHB8fyfj4SD4+P53DF2v51uZ8Hn3+EJ9ZPI7v3DOd4CBt2Kis46mr7zXgc87WUjcB9caYCmA/MEVEMkQkBHjAue6g2WzCPbNT2f61ZSyeGMfjfzrOD3ecRKeeVd5k0/4SPv+b/aSNCufNry7jC7dmXJMs+po7bjR/+j+38Fe3TeT3ey/yl88doLm900MRK3UtdzWrfRF4H5gqIqUi8oiIPCoijzpX2QacBU4DTwNfBjDGdAKPAduBQmCTMSZ/ODHER4XyzOcX8qkF6fw49zQ/fOuUi/8qpdzjjwdL+Yc/HeOWSfH84dGbGRc/+IYaIXYb37gzmyfvm8Hbpy7xyG8P0NrRNYLRKjUwd7WSevAGnxvgKwN8to3uhOKy4CAb//HxWQjCT3JPExsezJeWTnTHrpUalh35lfzDH4+yZPIYnv7cAsKCg4a1n4cWjyc8OIi/2XSUr208ws8emkeQzbsaeSj/53cVoiLCv90/kztnJPPktkLyTlZbHZIKUCcrr/L1l44wMy2Wpz43f9jJosf989L5p7un8WZ+Jf+146SbolRq8PwuYQAE2YT/+tRspiZF839fPMyFy01Wh6QCTENrB3/53AEiQ+089bkFRIS4p33JI0syeHDRWH6ed4bt+ZVu2adSg+WXCQMgIsTO059bAMBXNx6hs6v/9u5KjYRvvXqCsroWfvEX80mKCXPrvr99z3Rmp8fyd384Snldi1v3rdT1+G3CABgbF8GT983kSEkdP9t1xupwVIDYfKSMV4+U89VVU5g/frTb9x9qD+InD86jy2H4+z8exeHQFoHKM/w6YQB8bHYq985J5ce5xeSX11sdjvJzDW2Gb23OZ/740Xx5+aQRO864+Aj+6e5pvHv6Mr/7wKVZN5UaNL9PGADfuWcGoyOC+eYrJ+jSX2NqBP2+qI2W9i7+4+MzsY9wJ7sHFo5l6ZQx/GD7SaoaWkf0WEpBgCSM2Ihg/vGuaRwpqeP3e/XXmBoZ7xRf4oOKLr68YhKTE0d+iBoR4bvrZ9De5eC7W647XqdSbhEQCQNg/ZxUbp0cz3/uOEVdc7vV4Sg/09nl4DuvF5AUIfyfEayK6mvCmEi+snwyW45V8N7pGo8dVwWmgEkYIsK37p7O1dYO/ien2OpwlJ95Ye9FTlc38umpIYTaXetvMVR/ddtE0kaF8+S2Qr0BrkZUwCQMgKnJ0TywaBy/e/8CZy41Wh2O8hP1LR38aOcpbpkUz9xEzyYLgLDgIP5h7VTyyxt45XCZx4+vAkdAJQyAr6/OJNRu07GmlNv86p2z1DV38P/WZVs26dHHZqUyOz2W/9xxUseaUiMm4BJGQnQoX1ySwdZjFRSUN1gdjvJxlxvbeGbPOe6amcKMtFjL4rDZhH9Ym0VFfSsb9120LA7l3wIuYQB8aelEYsLsWspQLvtF3hlaOrr4+u2ZVofCLZPiWZwRx8/yzmgpQ42IgEwYseHBbFg2kZ2FVZwo0858anhqGtt4fu8F7p2bxuTEKKvDQUT42zVTuXS1jee1M58aAQGZMAA+d8sEokLtPPX2WatDUT7qmT3naOt08JUVk60O5UOLMuK4eWI8v3rnHO2dOn6acq+ATRgxYcF8ZvE4th6voORKs9XhKB9T39zBc+9fYN2MFCYlWF+66G3DbROpbGhly7Ebznas1JC4a8a9tSJyUkROi8gT/Xz+9yJyxPk4ISJdIhLn/Oy8iBx3fnbAHfEM1hdunYBN4Nd7znnysMoPPL/3Ao1tnXx5hec66Q3W8swEpiZF89TbZ3WqYuVWLicMEQkCfgbcCUwDHhSRab3XMcb8wBgzxxgzB/gGsNsYc6XXKiucny9wNZ6hSIkN557Zaby0v4TaJu39rQano8vB796/wJLJY5ieal3LqIGICF9amkFR5VXeLtbe38p93FHCWAScNsacNca0AxuB9ddZ/0HgRTcc1y02LJtIS0eXjvipBu2NE5VUNrTyhVsnWB3KgNbPSSMpJpSn9R6dciNxtcgqIp8A1hpjvuR8/1lgsTHmsX7WjQBKgck9JQwROQfUAgb4X2PMUwMcZwOwASAhIWH+pk2bXIq7tx8dbOVsXRf/tTyCkCD3dbxqbGwkKsq76rf7o3EOzXffb6Gxw/DvS8Ox9emo5y0xAmw9284fTnXwnVvCGB/z0R7o3hTn9Wic7rVixYqDLtXkGGNcegCfBH7V6/1ngZ8MsO6ngdf7LEt1PicCR4FlNzpmZmamcaf3z9SY8Y9vMb97/7xb97tr1y637m+kaJyDd/hirRn/+BbzzJ6z/X7uDTH2qGtuN9P+6Q3zf188dM1n3hTn9Wic7gUcMC5837ujSqoUGNvrfTowUPOMB+hTHWWMKXc+VwOv0F3F5VGLM+KYlhLDC3sv6k1CdV2/efccUaF2PjE/3epQbig2PJgHFo1j67EKqq/qfBnKde5IGPuBKSKSISIhdCeF1/quJCKxwG3A5l7LIkUkuuc1sAY44YaYhkREeOimcRRWNHC4pM7Th1c+oqqhla3HKvjkgnSiw4KtDmdQHlo8jk6H4Q8HSq0ORfkBlxOGMaYTeAzYDhQCm4wx+SLyqIg82mvV+4AdxpimXsuSgD0ichTYB2w1xrzpakzDsX5OGpEhQbzwgY7Do/r3/AcX6DKGz98ywepQBm1iQhQ3T4znxX0XdbZJ5TK39MMwxmwzxmQaYyYZY550LvulMeaXvdb5rTHmgT7bnTXGzHY+pvdsa4WoUDv3zk1jy7FynWBJXaOjy8GL+0pYOTWR8fGRVoczJJ9ZPI7S2hbeLr5kdSjKxwVsT+/+PLR4PG2dDv50SOcUUB+1q6iamsY2Hlg0zupQhuyO6cnER4bw+71aelau0YTRy7TUGOaOG8ULey/ozW/1ES/tLyExOpQVUxOsDmXIQuw2PrlgLLlF1VTW681vNXyaMPp4aPF4zl5q4oOzV268sgoIVQ2t7DpZzcfnp2MP8s0/mQcXjaXLYXhpf4nVoSgf5ptX/wi6e1YKMWF2XtirPb9Vtz8eLMVh4FMLxt54ZS81Pj6SpVPGsHH/RTq7dBRbNTyaMPoICw7i/nnp7Cioor65w+pwlMUcDsOmAyUszogjY4xv3ezu6zOLxlFR36o3v9WwacLox8fnpdPe6WDr8QqrQ1EW23vuChcuN/Pphb5buuixKjuJ0RHB2qhDDZsmjH7MSIthSmIUfzqknZ0C3aYDJUSH2blzRorVobgsxG7jntmpvFVQRVOHNupQQ6cJox8iwsfnp3PwQi3na5puvIHySw2tHWw7XsH6OamEhwTdeAMfcJ+z9HygstPqUJQP0oQxgHvnpCECLx/W4nugeuN4BW2dDj4x3/ero3rMTo9lYkIk75ZrwlBDpwljAMmxYSyZPIaXD5Xi0CEVAtLmI+VMiI9gdrr3TZI0XCLCx+elc6rWoVMTqyHThHEd989Lo7S2hf3ntU9GoKmsb+X9s5dZPycNEffNkeIN7p2bBsArWnpWQ6QJ4zrumJ5MZEgQL2urkoDz+tFyjPnzl6s/SRsVTnacjZcPleqIBmpINGFcR0SInTtnprD1eAUt7V1Wh6M86NUjZcxOj/X5vhcDuSXVzvnLzRy6WGd1KMqHaMK4gfvmptHY1smuk9VWh6I85HT1VfLLG1g/x/9KFz0WJNsJC7bxymFtOh4oOtzQw18Txg3cNDGeMVGhbDk20CSCyt+8ergcm8Dds32/78VAwu3Cquwk3jheqUOFBIg9p2tc3ocmjBsIsgnrZiaTU1hNY5s2RfR3xhg2Hy3j1sljSIwOszqcEfWxWSlcbmrXgTYDxGtHXP/R65aEISJrReSkiJwWkSf6+Xy5iNSLyBHn41uD3dYb3D0rlbZOBzmFVVaHokbYoYu1lFxp4V4/ro7qsXxqIpEhQVp6DgCtHV3syK90eT8uJwwRCQJ+BtwJTAMeFJFp/az6jjFmjvPxL0Pc1lILxo8mOSaM14/q2FL+7vWjFYTabdwxI9nqUEZcWHAQq6cl8WZ+pVvqt5X3yjt5iSY3NNxxRwljEXDaOd1qO7ARWO+BbT3GZhPumpXC7lPV1LfoCLb+yuEwvHGiguVTE4gKtVsdjkfcPSuVuuYOt9RvK++19XgFoyOCXd6PO/4q0oDes7KUAov7We9mETkKlAN/Z4zJH8K2iMgGYANAQkICeXl5rkc+BKmdXXR0GX78pzyWpg/uxDc2Nno8zuHQOLudqu2iqqGNjKC6YR/H586lwxBuh1/vOIxUhFod1jV87nx6ofYuw44TzdyUYueIi/tyR8Lorxts395Ah4DxxphGEVkHvApMGeS23QuNeQp4CmDq1Klm+fLlw413WG4zht+c3MXp9ij+afmiQW2Tl5eHp+McDo3Tuf/X8gmxX+Sxjy8fdgnDF8/lupqj7Cio5OYlSwm1e9cgi754Pr3Nmycqaes6yJfumMdzLu7LHVVSpUDv0dnS6S5FfMgY02CMaXS+3gYEi8iYwWzrLUSEu2el8u7pGq40tVsdjnKzD6ujMgOnOqrH3bNSuNrayTuntFrKH21zVkfdPDHe5X25I2HsB6aISIaIhAAPAK/1XkFEksU5II+ILHIe9/JgtvUmd89KodNhePOE660NlHc5XFJLVUMb62b6b9+Lgdw6eQyx4cHaWsoPtXZ0kVNYxR3Tk90yH73LezDGdAKPAduBQmCTMSZfRB4VkUedq30COOG8h/Fj4AHTrd9tXY1ppExPjWFCfARvnNDWUv5m67FKQuw2VmUnWh2Kx4XYbaydnsxbBVW0dugQOP5k96nu1lHu+iHklrK3s5ppW59lv+z1+qfATwe7rbcSEe6Ykcyv3zlHfXMHsW5odaCs11MdtWxKAtFhgfl/eufMZF46UMJ7Z2pYmZVkdTjKTbYdr2BURDA3T3K9Ogq0p/eQrZ2eTKfDkFOknfj8xeGSOirqW7lrlv/3vRjILZPGEB1q1+pWP9La0cXOgirumJZMsBuqo0ATxpDNTh9FckyY/mH5kW3HKwgJsrEqO3B/WfdUx71VUKVjS/mJt3uqo2a5776cJowhstmEtTOSu+sGdWwpn2dMdyOGpVPGEBOg1VE91s5Ipra5g306YZhf2FFQRXSYnVvcVB0FmjCG5Y7pybR1Oth96pLVoSgXFVQ0UFbXwh3TA7c6qseyzATCgm1s19Kzz+vs6h77blVWotuqo0ATxrAsnDCauMgQrZbyAzvyq7AJAdk6qq+IEDu3ZSawPb9K57H3cQcu1FLb3MEaN/8Q0oQxDPYgG7dnJ5FbVE1bpzZD9GU7CqpYMD6O+CjvGxbDCmtnJFPZ0MrR0jqrQ1Eu2JFfRYjdxrLMBLfuVxPGMK2dkUxjWyfvnb5sdShqmEquNFNY0cDt0wL3ZndfK7OSsNuEN90wFLayhjGGHQWVLJk8xu2jFmjCGKZbJscTHWrXTnw+7K2C7qbRmjD+LDY8mFsmj2H7iUqM0WopX1RYcZXS2hbWjMB1rQljmELtQazISmRnYTVdWt/rk3YUVDI1KZoJYyKtDsWrrJ2ezPnLzRRVXrU6FDUMOwoqEWFEmolrwnDB7dOSuNLUzuGLtVaHooaotqmdfeeuaOmiH7dPS0IEdhZo51Rf9FZBFfPHjSYh2v335TRhuOC2qQnYbcJbOnWrz8ktqsZhYM10TRh9JUSHMmfsKHbqde1zSmubyS9vGLHrWhOGC2LCgrlpYvyHdeHKd+woqCQ5JoyZabFWh+KVVmcncbS0nqqGVqtDUUPw5/tyI9OvSBOGi26flsTZS02cudRodShqkFrau9h96pKz6qW/ObxUT1VdTmG1xZGoodiRX0VmUhQZI3RfThOGi3o6fOVo8d1n7DldQ2uHQ6ujrmNKYhTj4iK0WsqH1Da1s+/8FdaMUOkCNGG4LH10BNkpMews0F9ivuKtgkqiw+wsznDfGDv+RkRYnZ3EntM1NLfrmGm+ILeou8XmSP4QckvCEJG1InJSRE6LyBP9fP6QiBxzPt4Tkdm9PjsvIsdF5IiIHHBHPJ52e3YiBy5c0albfUCXw7CzsJoVUxMJsevvpetZPS2R9k4H7xTr1K2+wBP35Vz+ixGRIOBnwJ3ANOBBEZnWZ7VzwG3GmFnAd4Gn+ny+whgzxxizwNV4rLB6WhIOA7uKtJTh7Q5eqOVKU7tWRw3CwglxxITZtXmtD+i5L7dm+sjel3PHT6xFwGljzFljTDuwEVjfewVjzHvGmJ7OCh8A6W44rteYkRpLUkyo1vf6gB35lYQE2bjNzWPs+KPgIBsrshI/rOpQ3uvD+3IjeP8C3DNFaxpQ0ut9KbD4Ous/ArzR670BdoiIAf7XGNO39AGAiGwANgAkJCSQl5fnSsxulx3bRW5hJTtydhES1J3hGxsbvS7O/gRKnMYYXj/UwtTRNg5+8K77AuvF385lqunkclM7z2zOZcrooJEPrA9/O58j5fkTbYTbobXkOHllI9jyzxjj0gP4JPCrXu8/C/xkgHVXAIVAfK9lqc7nROAosOxGx8zMzDTeJrewyox/fIvZVVT14bJdu3ZZF9AQBEqcp6uvmvGPbzHPvXfOLfH0x9/OZX1Lu5n0ja3me28UjmxAA/C38zkSHA6HWfivb5kvP3/whusCB4wL3/fuqJIqBcb2ep8OlPddSURmAb8C1htjPhzi1RhT7nyuBl6hu4rL59w8KZ7w4CCtlvJiuc4+BSuydO6LwYoJC2bxxDi9j+HF8ssbqL7axkoPXNfuSBj7gSkikiEiIcADwGu9VxCRccDLwGeNMad6LY8Ukeie18Aa4IQbYvK4sOAglmWOYWdBtY7y6aVyiqrISo4mfXSE1aH4lNXZSRRXN3K+psnqUFQ/cgqrEYHlU0f+vpzLCcMY0wk8Bmynu7ppkzEmX0QeFZFHnat9C4gHft6n+WwSsEdEjgL7gK3GmDddjckqt0/rnnzmRFmD1aGoPupbOth/vtYjv8L8zWrnqKdaevZOuUVVzBk7yiOTgLlldg1jzDZgW59lv+z1+kvAl/rZ7iwwu+9yX7ViagI2gbcKq5iZrmMUeZO3T12iy2F0KtZhGBsXQVZyNDsLq/jS0olWh6N6uXS1jaOl9fzdmkyPHE97LrlRfFQo88eP1vpeL5RbVM3oiGDmjB1tdSg+aXV2EvvP11LXrJ1Tvcmuk9335VZmeaZfkSYMN1uVnURBRQPldS1Wh6KcuhyGvJPVLJ+aSJBNBxscjtXTkpzn8ZLVoahecgurSYkNIzsl2iPH04ThZqucdeS52uvbaxwpqaW2uUPvX7hgVlosY6JCydHr2mu0dXbxTvElVmYlemzUZU0YbjbZOcqnjl7rPXIKqwmyCcu0d/ew2WzCyqwEdp+spqPLYXU4Cth37gpN7V0evS+nCcPNRIRV2Ym8e+YybZ3avNYb5BZVs3DCaGLDg60OxaetzEqiobWTA+d1SmJvkFNYTajdxs0Tx3jsmJowRsCqrCTaOx0UXOmyOpSAV1rbTFHlVVZ56KagP1s6ZQwhQTZyi7T0bDVjDDlFVdw6eQzhIZ4bskUTxghYlBFHVKidw9WaMKzWM4LwSm1O67LIUDs3TYrX+xhe4MylRkqutHj8vpwmjBEQYrexLHMMxy514dBRPi2VW1TN+PgIJo7QlJWBZlVWImcvNXFOe31bqmfqXE0YfmJVVhJ1bYYT5fVWhxKwmts7effMZY+2IvF3PV9Q2qjDWjlF1WSnxJA6Ktyjx9WEMUJWZCUi/PmXgPK8905fpr3Tofcv3GhsXASZSVHabNxC9c0dHLxQ+2ETfk/ShDFC4iJDmDTKRo7eILRMTlE1kSFBLMqIszoUv7IqO4l9567Q0NphdSgBaXdx9zA3VtyX04QxguYkBnGirIHK+larQwk4xhhyi6pYlpmgc3e72aqsRDodhndO6VzfVsgtrCI+MoTZ6aM8fmz9SxpBcxK6x3bU4rvn5Zc3UNXgmTkCAs3ccaMZFRGspWcLdHY5yDt1idumJlgyzI0mjBGUFiWkjw7XG4QW6GlOu3yqJgx3C7IJK6Ymknfyks717WGHS+qoa+6w7L6cJowRJCKsykpkz+kaWtq1T4Yn5RRVM3vsKBKiR36OgEC0KjuRK03tHCnRXt+elFNYjd0mLM30XO/u3jRhjLBV2Um0dTp474zW93pK9xwBdZa0IgkUS6ckYLeJtgL0sNyiKhZlxBETZs0wN25JGCKyVkROishpEXmin89FRH7s/PyYiMwb7La+bvHEOCJDgtipf1gek3eyGmM836kpkMSGB7NwQpzen/OgkivNnKpqtPS6djlhiEgQ8DPgTmAa8KCITOuz2p3AFOdjA/CLIWzr00LtQSydkkBuUZXO9e0huUXVJMWEMj01xupQ/Nqq7ESKKq9SWttsdSgBoSc5r8q2rl+RO0oYi4DTxpizxph2YCOwvs8664HnTLcPgFEikjLIbX3equxEqhrayC/Xub5HWnung3eKa1iZlaS9u0fYSp37xaNyiqqZOCaSDAuHuXHHnN5pQEmv96XA4kGskzbIbQEQkQ10l05ISEggLy/PpaA9obGxkby8PELbDAL8+o29rJ8cYnVY1+iJ09sNJs6Cy100tnWS2FlFXt5lzwTWiz+dy8FIjhD+8G4R49rOu7yv/gTa+RxIa6fhveJmVo2zW3o+3JEw+vsZ17fuZaB1BrNt90JjngKeApg6dapZvnz5EEK0Rl5eHj1x/ubMu5xpNSxfvsTaoPrRO05vNpg43369gBD7Bf7q3uVEhLjj8h4afzqXg3F3YwHPvX+BhTcvITLU/ec70M7nQLbnV9JpDvLw7fO5ZbI1LaTAPVVSpcDYXu/TgfJBrjOYbf3C6uwkjpXWU92gvb5HSs8cAbdMirckWQSildmJtHc52HNaWwGOpNzCaqJD7SyYYO0wN+5IGPuBKSKSISIhwAPAa33WeQ34nLO11E1AvTGmYpDb+gWt7x15Z2uauHC5WZvTetDCCXFEh9nJ1VaAI8bhMOSerGZp5hjLh7lx+ejGmE7gMWA7UAhsMsbki8ijIvKoc7VtwFngNPA08OXrbetqTN4oKzmatFHh2rx2BPV8aa3QhOExwUE2bstMIPdktc79MkJOlNdz6WqbV4y67JZyuzFmG91JofeyX/Z6bYCvDHZbfyQirMxK5I8HS2nt6CIs2HPTKgaKnKIqspKjSR8dYXUoAWVVdiJbjlVworyeWRYMiOfvcgqrEYHlUxOsDkV7envSquxEWjq6eP+M51vv+Lv6lg72n6/VznoWuC0zEZugpecRkltUzdyxo4iPsn6YG00YHnTTxHgiQoJ0lM8RsPtU90B4VnZqClRxkSHMGzeaXL2u3a6qoZXjZfVec11rwvCgsOAglkweQ25htfb6drPcwiriIkOYM3aU1aEEpJXZiTr3ywjY9WHvbu8oOWvC8LDV2UmU17dSUKG9vt2lZ46A5RbNEaC6r2uAXSe1WsqdcoqqSRsVztSkaKtDATRheNzyrO4bV9oM0X2sniNAwZTEKJ37xc1aO7rYU1zDyqxErxnmRhOGhyVGhzF77Ch2an8Mt7F6jgD10blfWjt07hd3+ODsZVo6uiyZu3sgmjAssCorkaMldVy62mZ1KH7B6jkCVLeV2Um0dji0FaCb5BZVEx4cxM0T460O5UOaMCzQcwNrl5YyXOYNcwSobjdNjNNWgG5ijCGnsJpbJ4/xqj5bmjAsMC0lhpTYMP3DcoOeOvPVXtLsMJB1z/2irQDd4VRVI2V1LV7TOqqHJgwL9PT6fqdY63tdlVNUzcSESCZYOEeA+rNVWd2tAAsrrlodik/r+THpbSVnTRgWWZ2dRHN7Fx+c1fre4Wps62Tv2Ss62KAX+bAVoJaeXZJbWM3MtFiSYsKsDuUjNGFY5OZJ8YQF23T0WhfsKa6hvcvBSm1O6zV6WgHm6HU9bFea2jl00TuHudGEYZHuXt8J5Gh977DlFFYRHWZnwYTRVoeielmVlciRkjpqGrUV4HDknazGYbynd3dvmjAstDo7kbK6Fk5WaX3vUDkchl0nq7ktM4HgIL2MvcnKrESM0VaAw5VTVE1CdCgzUmOtDuUa+pdmoZ4iZ472+h6yY2X11DS2a+soLzQ9NYakmFCtbh2Gji4Hb5+8xMqpidi8cJgbTRgWSowJY1Z6LDt1OIUhyy2swiZwW6b1cwSoj+puBZjE26cu0d7psDocn7L//BWutnV6Ve/u3lxKGCISJyJviUix8/maymQRGSsiu0SkUETyReSrvT77toiUicgR52OdK/H4opVa3zssOUXVzB8/mtGRIVaHovqxOjuRpvYu9p27YnUoPiW3sJqQIBtLJnvnMDeuljCeAHKMMVOAHOf7vjqBvzXGZAM3AV8RkWm9Pv+RMWaO8+H3M+/1tTo7Set7h6iivoX88gZtHeXFbpk0hlC7TUvPQ5RbVM1Nk+KJDHXLZKhu52rCWA8863z9LHBv3xWMMRXGmEPO11fpnrs7zcXj+g2t7x26XC+bI0BdKzwkiFsnjyGnqEpbAQ7S2UuNnK1pYrUXX9fiyn+miNQZY0b1el9rjBmwjaOITADeBmYYYxpE5NvA54EG4ADdJZHaAbbdAGwASEhImL9p06Zhx+0pjY2NREVF3XC9355o44OKTn6yKoJgC250DTZOq/XE+aODrZQ1OvjBsnCvGfa5h6+dy5GUe7GD5wra+bcl4aRGDe+3aSCdzzfPdbDxZDs/WBZOQsTI3F5esWLFQWPMgmHvwBhz3QewEzjRz2M9UNdn3drr7CcKOAjc32tZEhBEd0nnSeCZG8VjjCEzM9P4gl27dg1qvZ0FlWb841vM7pPVIxvQAAYbp9V27dplmto6TOY3t5l/3nzC6nD65UvncqSV1zWb8Y9vMb/MOz3sfQTS+fzkL94za//7bdeDuQ7ggBnEd+xAjxumMWPMamPMjH4em4EqEUkBcD73W68iIsHAn4AXjDEv99p3lTGmyxjjAJ4GFg0t3fmHnvperZa6sbdP1dDW6WDNNL1/4e1SYsOZlhKjzcYHoaaxjQMXrnj9de1quec14GHn64eBzX1XkO46g18DhcaYH/b5LKXX2/voLrkEnPCQ7rm+dxZqfe+N7CioJDY8mIUZcVaHogZhVXYiBy5coa653epQvFpuYXfv7jXT/TthfA+4XUSKgdud7xGRVBHpafF0K/BZYGU/zWe/LyLHReQYsAL4uovx+KyV2YmU1rZQXN1odSheq8vRPUfAqqxE7d3tI1ZmJeIwsPvUJatD8Wo7CipJG9VdIvNmLrXdMsZcBlb1s7wcWOd8vQfo986kMeazrhzfn6zKSuKbnGBnYRWZXjLhu7c5WeugvqWDNdOTrQ5FDdLs9FGMiQohp7Ca9XO0cWR/mto6ebu4hocWj/O6Rhx96c80L5EcG8aMtBhytb53QIeqOgm121imc3f7DJtNWDE1kbyT1XR0aa/v/rxT3N0jfs007/8hpAnDi6zMSuLQxVquNGl9b1/GGA5Vd7F0SgIRId7ZqUn1b1V2Ig2tnRw432+L+YC3I7+KURHBLPSBUZc1YXiR1dnd9b3a6/ta+eUNXGk1Xn9TUF1ryZQEQrTXd786uhzkFFWzKisJuw/cl/P+CAPIjNRYEqK113d/duRXIujc3b4oKtTOksljePNEpbYC7GP/uSvO+3K+cV1rwvAiNpuwKiuR3TrK5zW251eROdpGnA426JPWTk+mrK57DDD1ZzsKqggLtrFsim+MuqwJw8uszEqksa2T/ed1lM8e52uaOFl1lXlJeu/CV63KTsQmsD2/0upQvIYxhh35lSydkkB4SJDV4QyKJgwvs2TKGELsNt4q0PreHj3nYl6ib/xRqWvFR4WyKCOON09owuiRX95AeX2r1/fu7k0ThpeJCLGzbMoYtudX4nBofS90d2qalhIzYgOyKc9YOz2Z4upGzlzSzqnQfV/OJrDKh+7L6V+gF7pzRgoV9a0cLa2zOhTLdY+xU+szNwXVwHo6XGq1VLcdBVUsnBDnU/flNGF4odXZSdhtosV3IKewCmPwiU5N6vpSR4UzOz2W7fla3XrhchNFlVd9btQCTRheKDYimFsnj+ENbYbI9vwq0keHk52iw6X4gzXTkzlaUkdFfYvVoVhqhzNp+tL9C9CE4bXunJHMxSvNAd0M8WprB3uKa7hjerLXj7GjBmftjO5f1DsCvJSx9XgFM9JiGBsXYXUoQ6IJw0utmZ5MUIBXS+0srKK9y8G6mSk3Xln5hEkJUUxOjOKNExVWh2KZ0tpmjpTU+eR1rQnDS8VFhrA4I45tJyoCtlpq67FKUmLDmDt2lNWhKDdaNzOFfeeuUH211epQLPHG8e4fgXdpwlDudOfMFM5eagrIOTIaWjt4+9Ql1s1MwWbBPOdq5Nw9KwWHIWBLzz3VUePjI60OZchcShgiEicib4lIsfO53+EWReS8c6KkIyJyYKjbB6o7pichAtuOB17xfWdBd3XUXbN871eYur7MpGgyk6LYcizwrmtfro4C10sYTwA5xpgpQI7z/UBWGGPmGGMWDHP7gJMYHcbC8XEfFmEDybbjFaRqdZTfumtmKvvPX6GqIbCqpXy5OgpcTxjrgWedr58F7vXw9n5v3cxkTlZdpbjqqtWheEx9Swdvn6ph3cwUbR3lp+6alYwxgVd63nq8gumpvlkdBSCu3FAVkTpjzKhe72uNMddUK4nIOaAWMMD/GmOeGsr2zs82ABsAEhIS5m/atGnYcXtKY2MjUVFRLu2jrs3B13e1cPekYD4+ZWR6hLojTnd6t6yDp4+38483hTF51J/Hj/K2OPvjCzGCd8T5j3uaCbcL37wpfMB1vCHOwRhMnDUtDv5udwufyAzm7onW9O5esWLFwT61PENjjLnuA9gJnOjnsR6o67Nu7QD7SHU+JwJHgWXO94Pavu8jMzPT+IJdu3a5ZT+fefp9s+z7ucbhcLhlf325K053+eJv9plb/j3nmn+vt8XZH1+I0RjviPMnOafM+Me3mLLa5gHX8YY4B2MwcT61+4wZ//gWc76mceQDGgBwwAziO3agxw2rpIwxq40xM/p5bAaqRCQFwPnc78w/xphy53M18AqwyPnRoLYPdPfMTuXC5WaOl9VbHcqIq2/p4O3iS9w5Qzvr+bueG7+BUi21xcero8D1exivAQ87Xz8MbO67gohEikh0z2tgDd0llEFtr2Dt9BSCg4TXjpRbHcqI255fSUeX0dZRAWBiQhTTUmLYGgAJ43xNE0dL6rhndqrVobjE1YTxPeB2ESkGbne+R0RSRWSbc50kYI+IHAX2AVuNMW9eb3v1UbERwdyWmciWYxV+P+T5q4fLmBAfwRxtHRUQ7p6dwuGLdZRcabY6lBH16pEyROCeOQGcMIwxl40xq4wxU5zPV5zLy40x65yvzxpjZjsf040xT95oe3Wte+akUtnQyj4/nomvor6F989e5t65aVodFSA+Nqv7C/TVw2UWRzJyjDFsPlLOTRnxpMQOfIPfF2hPbx+xOjuR8OAgXjvqv9VSrx0pxxi4d06a1aEoDxkbF8GijDheOVzmt0PgHC2t51xNE/fN9f3rWhOGj4gIsXP7tCTeOF5BR5fD6nBGxCuHy5g7bhQTxvjuTUE1dPfPTeNsTRNHS/2zUcerh8sIsdtYO9O35r7ojyYMH/Kx2anUNnfwTvElq0Nxu8KKBooqr/rFrzA1NHfOTCHEbuOVQ6VWh+J2nV0OthwrZ3V2IjFhwVaH4zJNGD7ktswE4iJD+ONB//vDevVIGXab+OyQCWr4YsODuT07ideP+V/pec/pGmoa2/2mmlUThg8Jsdu4d04aOwuqqW1qtzoct3E4DJsPl3NbZgLxUaFWh6MscN/cNK40tbP7pH+Vnl89XEZseDDLpyZaHYpbaMLwMZ+Yn057l8Ovbn5/cO4ylQ2t3KvVUQHrtqndpedX/Ki1VFNbJ9vzq7hrVneVmz/wj39FAJmWGsP01Bj+cLDE6lDc5uVDZUSF2lmd7VvzGyv3CQ6y8bFZKbxVWEV9S4fV4bjF1uMVtHR0cb8f/RDShOGDPjk/nRNlDRRW+P5831dbO9h6rIKPzU4hPCToxhsov3XfvHTaOx1+M1TIS/tLmJQQyfzx/jPNjyYMH7R+ThrBQeIXN783HymnpaOLBxaOszoUZbHZ6bFkJkWxcd9Fq0NxWXHVVQ5eqOWBheP8qhOqJgwfNDoyhNXZSbx6uMznW5Vs3H+R7JQYZqXHWh2KspiI8OCicRwtreeEjw+0uXF/CcFBwv3z/Kc6CjRh+KxPLkjnclM7uUW+O8DvibJ6TpQ18OCisX71K0wN3/1z0wm129i433dLGW2dXbx8qJQ105L9rtWfJgwftWxKAonRoby033dvfr+47yJhwTbW+0kbdeW62Ihg7pqVwquHy2lq67Q6nGF5q6CK2uYOPr1wrNWhuJ0mDB9lD7LxwMKx7DpZ7ZMjfTa3d7L5SDnrZqYQG+77PWCV+3xm0Tga2zrZcsw3m45v3FdC2qhwlkweY3UobqcJw4d9ZvF4bCI8v/eC1aEM2ZZjFTS2derNbnWN+eNHMyUxit/v873Sc8mVZvacruFTC8Zis/lfNasmDB+WHBvGHdOTeGl/Ca0dXVaHMyQb911kUkIkCyf4T5ND5R4iwmcWj+NoSR0XGnzrun5x30VEuu8x+iNNGD7uszdNoK65g9d9qOf38dJ6Dl2s48FF/tXkULnPfXPTCLXbyCvxnfsYLe1d/H7fRdZMSyJ1lG/PezEQlxKGiMSJyFsiUux8vubnoohMFZEjvR4NIvI152ffFpGyXp+tcyWeQHTTxDgyk6J47v0LPjOfwK/3nCUq1M6n/PCmoHKPUREh3D0rlffKO6lv9o2e368cLqOuuYMv3pphdSgjxtUSxhNAjjFmCpDjfP8RxpiTxpg5xpg5wHygGXil1yo/6vncGLOt7/bq+kSEz948geNl9RwpqbM6nBuqrG9ly7EKPrVgrF8M96xGziNLMmjrghf2ef89OmMMz7x7jumpMSzKiLM6nBHjasJYDzzrfP0scO8N1l8FnDHGeP8V4EPum5tGVKid373v/af12ffP4zCGL9w6wepQlJeblhrD9Hgbz753nvZO7+6geqKmi9PVjTyyJMOvq1nFlWoMEakzxozq9b7WGDPgXUwReQY4ZIz5qfP9t4HPAw3AAeBvjTG1A2y7AdgAkJCQMH/Tpk3DjttTGhsbiYqK8sixni9oI6+kkx/cFs7osKH9DvBUnG2dhr/Z3UxWXBB/PTdsyNt78nwOly/ECL4T576SRn6eL/zlzBBuTfPeEun39zZS1mzjP28LJ9iLW0etWLHioDFmwbB3YIy57gPYCZzo57EeqOuzbu119hMC1ABJvZYlAUF0l3SeBJ65UTzGGDIzM40v2LVrl8eOdaGmyUz8xlbzr1vyh7ytp+J87r1zZvzjW8z+c5eHtb0nz+dw+UKMxvhOnLm5uWbND3ebO3602zgcDqvD6VdxVYMZ//gW8z87T1kdyg0BB8wgvmMHetzwp6gxZrUxZkY/j81AlYikADifrzdOxZ10ly6qeu27yhjTZYxxAE8Diwab6NRHjYuP4J7Zqbyw96JXTq7kcBieefc8s9Nj/Wr0TjWyRIRHlmZQVHmVPadrrA6nX8+8ex67DR5a7P99ily9h/Ea8LDz9cPA5uus+yDwYu8FPcnG6T66Sy5qmL68fBLN7V385t1zVodyjZyias7VNPFFP6/jVe63fk4qCdGhPP2O913X1Vdb+ePBUm5JtfvduFH9cTVhfA+4XUSKgdud7xGRVBH5sMWTiEQ4P3+5z/bfF5HjInIMWAF83cV4AtqUpGjumJ7Eb987z9VW72mKaIzhJ7nFjIuLYJ3O2a2GKNQexOdvmcDbpy6RX+5do9g+/fZZOrsc3JXhvfdX3MmlhGGMuWyMWWWMmeJ8vuJcXm6MWddrvWZjTLwxpr7P9p81xsw0xswyxtxjjPGPmVMs9NiKKTS0dvL8B94z2mfeqUscK63nKysmERykfUXV0P3FTeOJDrPzPzuLrQ7lQ5cb23j+g4usn5NGUmRgXNeB8a8MIDPTY1mWmcCv95z1iuFCjDH8OKeYtFHh3D/PP4dLUCMvNjyYR5ZksKOgymvmynj6nXO0dnbxlRWTrQ7FYzRh+KGvLJ9ETWM7z39gfb+MnMJqDl+s47GVk7V0oVzyhVsziAmz899eUMqobmjlt++d457ZqUxO9P7mye6if8F+aPHEeJZOGcNPd52mvsW6exldDsP3txcxcUwkn5yvpQvlmtjwYDYsm8jOwioOnL9iaSw/zi2ms8vwN7dnWhqHp2nC8FOPr82irrmDX+4+Y1kMrxwu41RVI3+7Zip2LV0oN3hkyUSSYkL5t22Flo2ddq6miY37Snhw0TjGx0daEoNV9K/YT81Ii+XeOak8s+ecJRMsNbV18oPtRcxKj2XdzGSPH1/5p/CQIP7m9kwOXazjjROVlsTw5NYCQu02/npV4Ny76KEJw4/9/dosbCI8ubXQ48f+6a7TVDW08c8fm679LpRbfWL+WLKSo3lyayHN7Z4d/nz3qUvsLKzmr1dNITF66MPb+DpNGH4sbVQ4j62czJv5lbx96pLHjnv2UiO/fucc989L017dyu2CbMJ3751BWV0LP8097bHjtnV28Z3X85kQHxGwg2dqwvBzX1qaQcaYSP5p8wmP/BpzOAxPvHyc0GAbT6zNGvHjqcC0cEIcH5+XztPvnOVU1VWPHPNnuac5e6mJb98znVB7kEeO6W00Yfi5UHsQ/37/TC5cbuYH20+O+PFe2HuBfeeu8I93ZZMYE3hFduU5/29dFjFhwfztpqN0dI3s8OdFlQ38PO8M989NY/nUxBE9ljfThBEAbpoYz+duHs9v3zvPB2cvj9hxztU08b03ilgyeQyfWqCz6amRFR8Vyr/eO4PjZfX8Im/kWgO2dnTxtY1HiA0P5p/unjZix/EFmjACxONrs5gQH8lXNx6mprHN7ftv6+zisd8fIthu4/ufmKU3upVH3DkzhXtmp/I/OcXsHaEfQ/++rZCiyqv85ydnMzoyZESO4Ss0YQSIyFA7P/3MXGqbO/j6S0focrivDbsxhn95vYD88gZ+8InZpI4Kd9u+lbqRJ++bwbi4CB578TDVV1vduu/NR8p49v0LPLIkgxVZgVsV1UMTRgCZnhrLv9wznXeKa/jO6/lu6/j0m3fP88Lei/zVsoncPi3JLftUarCiw4L5xV/M42prB3/53EG3Ne44fLGWv//jMRZNiONxbcABaMIIOA8sGseGZRN57v0L/NwN9b5bjpXz3a0F3DE9Sf+olGWykmP48QNzOV5ax2O/P+zyTfBTVVf50rMHSI4J45efnU+IXb8qQRNGQHpibRb3zE7lB9tP8uOc4mGXNDYfKeOrG4+wYPxo/vvTc7F58VzGyv+tmZ7Md++dQW5RNX/1u4PDHq35VNVVPvP0BwTZhN9+YSFxAX7fojdNGAHIZhN++KnZ3D8vjR++dYrf5rcP6Y/L4egesvxrL3Uni99+YRHhIYHZLl15l4cWj+fJ+2aw62Q1D/1qLxX1LUPaPreoivt//h42EV7ccBMTEwJnJNrBcClhiMgnRSRfRBwisuA6660VkZMiclpEnui1PE5E3hKRYuezdgv2EHuQjf/8xGy+vHwSu0s7uf/n73HoYu0Ntzt7qZHPPbOPH751ivvmpPHsFxcRGWr3QMRKDc5Di8fz0wfnUVjRwLr/eYfNR8pw3KCRR0NrB9/afIJHnj3AhDERvPqVW5mkyeIarv6lnwDuB/53oBVEJAj4Gd1TtJYC+0XkNWNMAfAEkGOM+Z4zkTwBPO5iTGqQbDbhH9ZmEdxQyu+L27j/5++xZloSDyway+KM+A8TQWtHF0dL6th0oJTXj5YTarfxr/fO4KHF47T5rPJKd81KISslmq9uPMxXNx7hV++c46HF41gzPfnDKiaHw1Bc3chrR8v4/d6L1Ld08PDNE3h8bZaWmAfgUsIwxhQCN/rSWAScNsacda67EVgPFDiflzvXexbIQxOGx81NtPOX65fwv7vP8MLei+woqEIEUmLCEBGqr7bS0WWIDAniUwvT+b8BOvCa8i2TEqLY/JUlvHq4jJ/tOs0TLx/niZePEx8ZQlhwEHXN7TS1d2ETWJ2dxGMrJzMrfZTVYXs1cUfTShHJA/7OGHOgn88+Aaw1xnzJ+f6zwGJjzGMiUmeMGdVr3VpjTL/VUiKyAdgAkJCQMH/Tpk0uxz3SGhsbiYry/mJt7zg7HYb8y12cr3dQ3WwQgZgQISPWxowxQYTbrStR+ML59IUYIfDiNMZwrsFB0ZUuKpsMXQ4It0NGrI2suCDiw127nesr53PFihUHjTED3j64kRuWMERkJ9DfhAbfNMZsHsQx+vuGGXKWMsY8BTwFMHXqVLN8+fKh7sLj8vLy8MU4V1sXynX5wvn0hRghMONc4Za99M9XzqerbpgwjDGufn+UAr0HFkoHyp2vq0QkxRhTISIpQLWLx1JKKTVCPNGsdj8wRUQyRCQEeAB4zfnZa8DDztcPA4MpsSillLKAq81q7xORUuBmYKuIbHcuTxWRbQDGmE7gMWA7UAhsMsbkO3fxPeB2ESmmuxXV91yJRyml1MhxtZXUK8Ar/SwvB9b1er8N2NbPepeBVa7EoJRSyjO0p7dSSqlB0YShlFJqUDRhKKWUGhRNGEoppQbFLT29PU1ErgInrY5jEMYANVYHMQgap/v4Qoygcbqbr8Q51RgTPdyNfXWY0ZOudG/3FBE5oHG6jy/E6Qsxgsbpbr4Upyvba5WUUkqpQdGEoZRSalB8NWE8ZXUAg6RxupcvxOkLMYLG6W4BEadP3vRWSinleb5awlBKKeVhmjCUUkoNik8lDBFZKyInReS0cw5wryAiY0Vkl4gUiki+iHzVufzbIlImIkecj3U32pcHYj0vIsed8RxwLosTkbdEpNj53O+shx6McWqvc3ZERBpE5GvecD5F5BkRqRaRE72WDXj+ROQbzuv1pIjcYXGcPxCRIhE5JiKviMgo5/IJItLS67z+0uI4B/x/tuJ8DhDjS73iOy8iR5zLrTyXA30Pue/6NMb4xAMIAs4AE4EQ4Cgwzeq4nLGlAPOcr6OBU8A04Nt0T11reYy9Yj0PjOmz7PvAE87XTwD/YXWcff7fK4Hx3nA+gWXAPODEjc6f8xo4CoQCGc7rN8jCONcAdufr/+gV54Te63nB+ez3/9mq89lfjH0+/y/gW15wLgf6HnLb9elLJYxFwGljzFljTDuwEVhvcUwAGGMqjDGHnK+v0j3vR5q1UQ3JeuBZ5+tngXutC+Uaq4AzxpgLVgcCYIx5G7jSZ/FA5289sNEY02aMOQecpvs6tiROY8wO0z0/DcAHdM9+aakBzudALDmf14tRRAT4FPDiSMdxI9f5HnLb9elLCSMNKOn1vhQv/FIWkQnAXGCvc9FjziqAZ6yu6nEywA4ROSgiG5zLkowxFdB90QGJlkV3rQf46B+jt51PGPj8efM1+0XgjV7vM0TksIjsFpGlVgXVS3//z954PpcCVcaY4l7LLD+Xfb6H3HZ9+lLCkH6WeVWbYBGJAv4EfM0Y0wD8ApgEzAEq6C66Wu1WY8w84E7gKyKyzOqABiLdU/reA/zBucgbz+f1eOU1KyLfBDqBF5yLKoBxxpi5wN8AvxeRGKviY+D/Z288nw/y0R80lp/Lfr6HBly1n2XXPZ++lDBKgbG93qcD5RbFcg0RCab7P+kFY8zLAMaYKmNMlzHGATyNh6ojrsd0z4aIMaaa7tkSFwFVIpIC4Hyuti7Cj7gTOGSMqQLvPJ9OA50/r7tmReRh4G7gIeOsyHZWSVx2vj5Id112plUxXuf/2avOp4jYgfuBl3qWWX0u+/sewo3Xpy8ljP3AFBHJcP7yfAB4zeKYgA/rMX8NFBpjfthreUqv1e4DTvTd1pNEJFJEonte030T9ATd5/Fh52oPA5utifAaH/n15m3ns5eBzt9rwAMiEioiGcAUYJ8F8QHdrQyBx4F7jDHNvZYniEiQ8/VEuuM8a02U1/1/9qrzCawGiowxpT0LrDyXA30P4c7r04q7+S60AlhH953/M8A3rY6nV1xL6C7KHQOOOB/rgN8Bx53LXwNSLI5zIt2tIo4C+T3nEIgHcoBi53OcF5zTCOAyENtrmeXnk+4EVgF00P0L7ZHrnT/gm87r9SRwp8Vxnqa7zrrnGv2lc92PO6+Ho8Ah4GMWxzng/7MV57O/GJ3Lfws82mddK8/lQN9Dbrs+dWgQpZRSg+JLVVJKKaUspAlDKaXUoGjCUEopNSiaMJRSSg2KJgyllFKDoglDKaXUoGjCUEopNSj/H0ADNl/DaGRKAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(sol.T[0])\n",
    "plt.xlim(0,200)\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the times in units of years"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1) Convert to seconds (SI unit) using equation above\n",
    "tt = 1/np.sqrt(6.67e-11 * 5.97e24 / (1.5e11)**3 ) \n",
    "# 2) Convert from seconds to years\n",
    "tt = tt / (60*60 * 24* 365.25) * np.diff(t)[0] # per time step (in years)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get solutions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1 = sol.T[0]\n",
    "y1 = sol.T[1]\n",
    "x2 = sol.T[2]\n",
    "y2 = sol.T[3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make animation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeIAAAHSCAYAAAAwk8gOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAX5ElEQVR4nO3de6zdZZ3v8c+DoBi2M2pbinKxCOeYA630yEaLRNgFQRzBMcYbByYQEjuKOGIQ5OIFCQ4RZvQQFQnR4Ekk1hE1csAjXrcXvMSWKSMIJXTCjMwwKhiFLTZQfM4fGwrYltLutfeXvfbrlTRZa/3WftZ3P6m8/a1bW+89AECNHaoHAIC5TIgBoJAQA0AhIQaAQkIMAIWEGAAK7VjxoPPnz++LFi2qeOjt9oc//CG77LJL9RhDzR5PP3s8/ezxzJht+7x69eq7e+8LNnesJMSLFi3KqlWrKh56u42Pj2dsbKx6jKFmj6efPZ5+9nhmzLZ9bq3925aOeWoaAAoJMQAUEmIAKCTEAFBIiAGgkBADQCEhBoBCQgwAhYQYAAoJMQAUEmIAKCTEAFBIiAGgkBADQCEhBoBCQgwAhYQYAAoJMQAUEmIAKCTEAFBIiAGgkBADQCEhBoBCQgwAhYQYAAoJMQAUEmIAKCTEAFBIiAGgkBADQCEhBoBCQgwAhYQYAAoJMQAUEmIAKCTEAFBIiAGgkBADQKGBhbi19rTW2j+31q4Z1JoAMOwGeUb8riS3DHA9ABh6Awlxa22PJK9J8ulBrAcAc0XrvU99kdauSnJhkmcleU/v/ZjN3GdFkhVJsnDhwgNXrlw55cedSRMTExkZGakeY6jZ4+lnj6efPZ4Zs22fly9fvrr3Prq5YztOdfHW2jFJft17X91aG9vS/Xrvlye5PElGR0f72NgW7/qUND4+ntk282xjj6efPZ5+9nhmDNM+D+Kp6UOSvLa1dkeSlUkOb619bgDrAsDQm3KIe+9n99736L0vSvKWJN/pvZ8w5ckAYA7wOWIAKDTl14gfq/c+nmR8kGsCwDBzRgwAhYQYAAoJMQAUEmIAKCTEAFBIiAGgkBADQCEhBoBCQgwAhYQYAAoJMQAUEmIAKCTEAFBIiAGgkBADQCEhBoBCQgwAhYQYAAoJMQAUEmIAKCTEAFBIiAGgkBADQCEhBoBCQgwAhYQYAAoJMQAUEmIAKCTEAFBIiAGgkBADQCEhBoBCQgwAhYQYAAoJMQAUEmIAKCTEAFBIiAGgkBADQCEhBoBCQgwAhYQYAAoJMQAUEmIAKCTEAFBIiAGgkBADQCEhBoBCQgwAhYQYAAoJMQAUEmIAKCTEAFBIiAGgkBADQCEhBoBCQgwAhYQYAAoJMQAUEmIAKCTEAFBIiAGgkBADQCEhBoBCQgwAhYQYAAoJMQAUEmIAKCTEAFBIiAGgkBADQCEhBoBCQgwAhYQYAAoJMQAUEmIAKCTEAFBIiAGg0JRD3Frbs7X23dbaLa21m1tr7xrEYAAwF+w4gDU2JDm9935Da+1ZSVa31r7Ze//FANYGgKE25TPi3vtdvfcbHr58X5Jbkuw+1XUBYC5ovffBLdbaoiTfT7K4937vnx1bkWRFkixcuPDAlStXDuxxZ8LExERGRkaqxxhq9nj62ePpZ49nxmzb5+XLl6/uvY9u7tjAQtxaG0nyvSQf7r1/+YnuOzo62letWjWQx50p4+PjGRsbqx5jqNnj6WePp589nhmzbZ9ba1sM8UDeNd1a2ynJl5JcubUIAwCPGsS7pluSzyS5pff+0amPBABzxyDOiA9J8jdJDm+trXn4z18NYF0AGHpT/vhS7/2HSdoAZgGAOcc3awFAISEGgEJCDACFhBgACgkxABQSYgAoJMQAUEiIAaCQEANAISEGgEJCDACFhBgACgkxABQSYgAoJMQAUEiIAaCQEANAISEGgEJCDACFhBgACgkxABQSYgAoJMQAUEiIAaCQEANAISEGgEJCDACFhBgACgkxABQSYgAoJMQAUEiIAaCQEANAISEGgEJCDACFhBgACgkxABQSYgAoJMQAUEiIAaCQEANAISEGgEJCDACFhBiYEffcc0+WLl2apUuXZrfddsvuu++epUuXZmRkJKecckr1eDn55JOz6667ZvHixVu8z+9///sce+yxOeCAA7L//vvniiuueNzxhx56KG9961tzzDHHbPKz3/jGN3LwwQen977xvkuXLs2PfvSjwf4izDpCDMyIefPmZc2aNVmzZk3e9ra35d3vfnfWrFmTiYmJXHrppdXj5aSTTsrXv/71J7zPJz/5yey333658cYbMz4+ntNPPz0PPPDAxuOXXHJJ9tprr83+7FFHHZUXvOAF+cxnPpMk+fjHP56DDjooL3/5y7d75g0bNmz3z/LUIcRAqfHx8Y1nkOedd15OPPHEHHXUUVm0aFG+/OUv58wzz8ySJUty9NFH58EHH0ySrF69OocddlgOPPDAvOpVr8pdd9015TkOPfTQPPe5z33C+7TWct9996X3nomJiTz3uc/NjjvumCS58847c+211+Y1r3nNFn/+Yx/7WC688MLcfPPN+cQnPpGPfOQjG8+UX/KSl+SNb3xjJiYmkiTnn39+DjrooCxevDgrVqzYeCY9NjaWc845J4cddlguueSSfPGLX8zixYtzwAEH5NBDD53yPjDzhBh4Slm3bl2uvfbafPWrX80JJ5yQ5cuX5+c//3me+cxn5tprr82DDz6Yd77znbnqqquyevXqnHzyyTn33HM3WefKK6/c+FT4Y/+84Q1v2O7ZTj311Nxyyy15/vOfnyVLluSSSy7JDjtM/mf0tNNOy0UXXbTx+uY873nPy2mnnZaDDz4473vf+/KnP/0pF1xwQb71rW/lhhtuyOjoaD760Y9ufKyf/exnuemmm/LHP/4x11xzzcZ1fve73+V73/teTj/99Jx//vm57rrrcuONN+bqq6/e7t+NOjtWDwDwWK9+9auz0047ZcmSJXnooYdy9NFHJ0mWLFmSO+64I2vXrs1NN92UI488Msnka63Pe97zNlnn+OOPz/HHHz/Q2a677rosXbo03/nOd7Ju3boceeSRecUrXpHvf//72XXXXXPggQfmBz/4wROu8Y53vCNnnXVWTjrppFxzzTX5xS9+kUMOOSRJ8sADD+Tggw9Oknz3u9/NRRddlPvvvz+//e1vs//+++fYY49Nkrz5zW/euN4hhxySk046KW9605vy+te/fqC/LzNDiIGnlGc84xlJkh122CE77bRTWmsbr2/YsCG99+y///758Y9//ITrXHnllbn44os3uX3ffffNVVddtV2zXXHFFTnrrLPSWsu+++6bvffeO7feemuuv/76XH311fna176We++9N+vXr88JJ5yQz33uc5usscMOO2z8nXrvOfLII/P5z3/+cfdZv359TjnllKxatSp77rlnzjvvvKxfv37j8V122WXj5csuuyw//elPc+2112bp0qVZs2ZN5s2bt12/HzU8NQ3MKi960Yvym9/8ZmOIH3zwwdx8882b3O/444/f+Oawx/7Z3ggnyV577ZVvf/vbSZJf/epXWbt2bV74whfmwgsvzJ133pk77rgjH/jAB3L44YdvNsJ/btmyZbn++utz++23J0nuv//+3HbbbRujO3/+/ExMTDzhzOvWrcvLXvaynH/++Zk/f35++ctfbvfvRw0hBmaVpz/96bnqqqvy3ve+NwcccMDAPgJ03HHH5eCDD87atWuzxx57bHx382WXXZbLLrssSfL+978/P/rRj7JkyZIcccQR+chHPpL58+dv92MuWLAgn/3sZ3PcccflxS9+cZYtW5Zbb701z372s/PWt741S5Ysyete97ocdNBBW1zjjDPOyJIlS7J48eIceuihOeCAA7Z7Hmq0R96JN5NGR0f7qlWrZvxxp2J8fDxjY2PVYww1ezz97PH0s8czY7btc2ttde99dHPHnBEDQCEhBoBCQgwAhYQYAAoJMQAUEmIAKCTEAFBIiAGgkBADQCEhBoBCQgwAhYQYAAoJMQAUEmIAKCTEAFBIiAGgkBADQCEhBoBCQgwAhYQYAAoJMQAUEmIAKCTEAFBIiAGgkBADQCEhhmE3MZF88IPJggU57PDDkwULJq9PTFRPBmRAIW6tHd1aW9tau721dtYg1gQGYGIiWbYsueii5O6703pP7r578vqyZWIMTwFTDnFr7WlJPpnk1Un2S3Jca22/qa4LDMDFFyfr1iXr1z/+9vXrJ2+/+OKauYCNBnFG/NIkt/fe/7X3/kCSlUn+egDrAlN16aWbRvgR69cnn/rUzM4DbGLHAayxe5JfPub6nUle9ud3aq2tSLIiSRYuXJjx8fEBPPTMmZiYmHUzzzb2ePAOu+eetCc43u++O9+z5wPl7/HMGKZ9HkSIN/e/877JDb1fnuTyJBkdHe1jY2MDeOiZMz4+ntk282xjj6fBvHmTrwlvQZs/354PmL/HM2OY9nkQT03fmWTPx1zfI8l/DmBdYKpOOSXZeefNH3vGM5K3v31m5wE2MYgQ/yzJf2ut7d1ae3qStyS5egDrAlN1xhnJPvtsPsa77DJ5HCg15RD33jckOTXJdUluSfJPvfebp7ouMAAjI8lPfpKceWayYEF6e8wrSb//fXLXXXWzAUkG9Dni3vvXeu//vfe+T+/9w4NYExiQkZHkQx9Kfv3rfO/b304OO2zy9oceSs45p3Y2wDdrwZzS2uSXeTziqqsmz5iBMkIMc81LX5q86U2PXj/jjKRv8kEHYIYIMcxFf//3yU47TV7+4Q+Tq72/EqoIMcxF++zz+I8uvfe9yYYNdfPAHCbEMFe9733JX/zF5MeY3vIWIYYig/hmLWA2WrAg+cIXkqVLk912q54G5iwhhrns6KOrJ4A5z1PTAFBIiIFH3XZb8pnPVE8Bc4oQA8kDDySnnprsv3/yt3+brF1bPRHMGUIMTH6m+Be/mHzn9EMPJWefXT0RzBlCDGz61Zdf+Upy/fV188AcIsTApNHR5LjjHr3uqy9hRggx8KgPf/jRr7788Y8nz4yBaSXEwKP23nvyTVuPOPvs5MEH6+aBOUCIgcc799zkL/9y8vJttyWf/nTtPDDkhBh4vHnzknPOefT6eecl991XNg4MOyEGNvXOdyZ77jl5+de/Tv7hH2rngSEmxMCmnvnM5IILJi8fdFByxBG188AQ848+AJt3/PHJc56THHPM5OeMgWkhxMDmPe1pybHHVk8BQ89T0wBQSIiBJ+ehh5Irrpj8og9gYDw1DWzdDTckJ56Y3HRT8rKXTcbY68YwEM6Iga2bN2/yyz2S5Kc/Tb70pdp5YIgIMbB1L3hB8nd/9+j1s8+e/DeMgSkTYuDJOeecyY8zJcnttyeXX147DwwJIQaenOc85/FfffmhDyX33ls3DwwJIQaevFNPTfbaa/Ly3XcnF19cOw8MASEGnrydd578N4sfccEFyQ47JAsWJB/8YDIxUTcbzFJCDGyb1752MsiP6H3y7Piii5Jly8QYtpEQA9vmH/8x+dOfNr19/fpk3TpPV8M2EmJg21x66ZY/urR+ffKpT83sPDDLCTGwbe65Z2rHgccRYmDbzJs3tePA4wgxsG1OOeXxb9Z6rJ13Tt7+9pmdB2Y5IQa2zRlnJPvss2mMd9558vYzzqiZC2YpIQa2zchI8pOfJGeeOfn54Uc+R3zmmZO3j4xUTwizin8GEdh2IyOTX3H5oQ9VTwKznjNiACgkxABQSIgBoJAQA0AhIQaAQkIMAIWEGAAKCTEAFBJiACgkxABQSIgBoJAQA0AhIQaAQkIMAIWEGAAKCTEAFBJiACgkxABQSIgBoJAQA0AhIQaAQkIMAIWEGAAKCTEAFBJiACgkxABQSIgBoJAQA0AhIQaAQkIMAIWEGAAKCTEAFBJiACgkxABQSIgBoJAQA0AhIQaAQkIMAIWEGAAKTSnErbWLW2u3ttb+pbX2ldbaswc0FwDMCVM9I/5mksW99xcnuS3J2VMfCQDmjimFuPf+jd77hoev/iTJHlMfCQDmjtZ7H8xCrf3fJF/ovX9uC8dXJFmRJAsXLjxw5cqVA3ncmTIxMZGRkZHqMYaaPZ5+9nj62eOZMdv2efny5at776ObO7bVELfWvpVkt80cOrf3/tWH73NuktEkr+9Pouyjo6N91apVWx38qWR8fDxjY2PVYww1ezz97PH0s8czY7btc2ttiyHecWs/3Ht/5VYWPzHJMUmOeDIRBgAetdUQP5HW2tFJ3pvksN77/YMZCQDmjqm+a/oTSZ6V5JuttTWttcsGMBMAzBlTOiPuve87qEEAYC7yzVoAUEiIAaCQEANAISEGgEJCDACFhBgACgkxABQSYgAoJMQAUEiIAaCQEANAISEGgEJCDACFhBgACgkxABQSYgAoJMQAUEiIAaCQEANAISEGgEJCDACFhBgACgkxABQSYgAoJMQAUEiIAaCQEANAISEGgEJCDACFhBgACgkxABQSYgAoJMQAUEiIAaCQEANAISEGgEJCDACFhBgACgkxABQSYgAoJMQAUEiIAaCQEANAISEGgEJCDACFhBgACgkxABQSYgAoJMQAUEiIAaCQEANAISEGgEJCDACFhBgACgkxABQSYgAoJMQAUEiIAaCQEANAISEGgEJCDACFhBgACgkxABQSYgAoJMQAUEiIAaCQEANAISEGgEJCDACFhBgACgkxABQSYgAoJMQAUEiIAaCQEANAISEGgEJCDACFhBgACg0kxK2197TWemtt/iDWA4C5Ysohbq3tmeTIJP8+9XEAYG4ZxBnxx5KcmaQPYC0AmFNa79vfz9baa5Mc0Xt/V2vtjiSjvfe7t3DfFUlWJMnChQsPXLly5XY/boWJiYmMjIxUjzHU7PH0s8fTzx7PjNm2z8uXL1/dex/d3LGthri19q0ku23m0LlJzklyVO/991sL8WONjo72VatWbXXwp5Lx8fGMjY1VjzHU7PH0s8fTzx7PjNm2z621LYZ4x639cO/9lVtYdEmSvZPc2FpLkj2S3NBae2nv/b+mMC8AzBlbDfGW9N5/nmTXR65vyxkxADDJ54gBoNB2nxH/ud77okGtBQBzhTNiACgkxABQSIgBoJAQA0AhIQaAQkIMAIWEGAAKCTEAFBJiACgkxABQSIgBoJAQA0AhIQaAQkIMAIWEGAAKCTEAFBJiACgkxABQSIgBoJAQA0AhIQaAQkIMAIWEGAAKCTEAFBJiACgkxABQSIgBoJAQA0AhIQaAQkIMAIWEGAAKCTEAFBJiACgkxABQSIgBoJAQA0AhIQaAQkIMAIVa733mH7S13yT5txl/4KmZn+Tu6iGGnD2efvZ4+tnjmTHb9vkFvfcFmztQEuLZqLW2qvc+Wj3HMLPH088eTz97PDOGaZ89NQ0AhYQYAAoJ8ZN3efUAc4A9nn72ePrZ45kxNPvsNWIAKOSMGAAKCfE2aq29p7XWW2vzq2cZRq21i1trt7bW/qW19pXW2rOrZxoWrbWjW2trW2u3t9bOqp5n2LTW9mytfbe1dktr7ebW2ruqZxpWrbWntdb+ubV2TfUsgyDE26C1tmeSI5P8e/UsQ+ybSRb33l+c5LYkZxfPMxRaa09L8skkr06yX5LjWmv71U41dDYkOb33/j+SLEvyDns8bd6V5JbqIQZFiLfNx5KcmcQL69Ok9/6N3vuGh6/+JMkelfMMkZcmub33/q+99weSrEzy18UzDZXe+1299xsevnxfJkOxe+1Uw6e1tkeS1yT5dPUsgyLET1Jr7bVJ/qP3fmP1LHPIyUn+X/UQQ2L3JL98zPU7IxLTprW2KMn/TPLT4lGG0f/O5AnRn4rnGJgdqwd4KmmtfSvJbps5dG6Sc5IcNbMTDacn2ufe+1cfvs+5mXyq78qZnG2Itc3c5pmdadBaG0nypSSn9d7vrZ5nmLTWjkny69776tbaWPE4AyPEj9F7f+Xmbm+tLUmyd5IbW2vJ5NOlN7TWXtp7/68ZHHEobGmfH9FaOzHJMUmO6D5fNyh3JtnzMdf3SPKfRbMMrdbaTpmM8JW99y9XzzOEDkny2tbaXyXZOclftNY+13s/oXiuKfE54u3QWrsjyWjvfTZ94fis0Fo7OslHkxzWe/9N9TzDorW2Yybf/HZEkv9I8rMk/6v3fnPpYEOkTf6/9P+T5Le999OKxxl6D58Rv6f3fkzxKFPmNWKeaj6R5FlJvtlaW9Nau6x6oGHw8BvgTk1yXSbfRPRPIjxwhyT5mySHP/x3d83DZ27whJwRA0AhZ8QAUEiIAaCQEANAISEGgEJCDACFhBgACgkxABQSYgAo9P8BLbaW2/xKl+IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def animate(i):\n",
    "    ln1.set_data([x1[i], x2[i]], [y1[i], y2[i]])\n",
    "    text.set_text('Time = {:.2f} Years'.format(i*tt))\n",
    "fig, ax = plt.subplots(1,1, figsize=(8,8))\n",
    "ax.grid()\n",
    "ln1, = plt.plot([], [], 'ro--', lw=3, markersize=8)\n",
    "text = plt.text(0.7, 0.7, '')\n",
    "ax.set_ylim(-5, 5)\n",
    "ax.set_xlim(-5,5)\n",
    "ani = animation.FuncAnimation(fig, animate, frames=200, interval=50)\n",
    "ani.save('plan.gif',writer='pillow',fps=30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3 Body Problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* https://arxiv.org/pdf/1303.0181.pdf\n",
    "* https://arxiv.org/pdf/1709.04775.pdf\n",
    "\n",
    "These two papers look at the three body problem under the following configuraton:\n",
    "\n",
    "* $m_1 = m_2 = 1$. $m_3$ can vary\n",
    "* At $t=0$ the initial conditions $x_1 = -x_2 = -1 $  and $x_3 = 0$ \n",
    "* At $t=0$ the initial conditions $y_1 = y_2 = y_3 = 0$\n",
    "* At $t=0$ the initial conditions $v_{x1} = v_{x2} = v_1$ where $v_1$ is some defined speed\n",
    "* At $t=0$ the initial conditions $v_{y1} = v_{y2} = v_2$ where $v_2$ is some defined speed\n",
    "* At $t=0$ the initial conditions $v_{x3} = -2 v_1 / m_3$ and $v_{y3} = -2v_2 / m_3$\n",
    "\n",
    "Thus the only parameters that are altered are $m_3$, $v_1$, and $v_2$. **These are states with zero angular momentum**. \n",
    "\n",
    "See Table I from paper 1 and Table I from paper 2 for interesting stable configurations (hard to find!). Here are two I found:\n",
    "\n",
    "1. From last row of Table I paper 2: $v_1 = 0.9911981217$, $v_2 = 0.7119472124$, and $m_3 = 4$\n",
    "2. From second row of Table I paper 1: $v_1=0.39295$, $v_2 = 0.09758$ and $m_3=1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PARARMS TO CHANGE\n",
    "m3 = 1\n",
    "v1 =  0.39295\n",
    "v2 = 0.09758\n",
    "\n",
    "# Everything else follows from paper\n",
    "m1 = 1\n",
    "m2 = 1 \n",
    "m3 = m3\n",
    "x1_0 = -1\n",
    "y1_0 = 0\n",
    "x2_0 = 1\n",
    "y2_0 = 0\n",
    "x3_0 = 0\n",
    "y3_0 = 0\n",
    "vx1_0 =  v1\n",
    "vy1_0 =  v2\n",
    "vx2_0 = v1\n",
    "vy2_0 = v2\n",
    "vx3_0 = -2*v1/m3\n",
    "vy3_0 = -2*v2/m3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the ODE system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dSdt(t, S):\n",
    "    x1, y1, x2, y2, x3, y3, vx1, vy1, vx2, vy2, vx3, vy3 = S\n",
    "    r12 = np.sqrt((x2-x1)**2 + (y2-y1)**2)\n",
    "    r13 = np.sqrt((x3-x1)**2 + (y3-y1)**2)\n",
    "    r23 = np.sqrt((x2-x3)**2 + (y2-y3)**2)\n",
    "    return [ vx1,\n",
    "            vy1,\n",
    "            vx2,\n",
    "            vy2,\n",
    "            vx3,\n",
    "            vy3,\n",
    "            m2/r12**3 * (x2-x1) + m3/r13**3 * (x3-x1), #mass 1\n",
    "            m2/r12**3 * (y2-y1) + m3/r13**3 * (y3-y1),\n",
    "            m1/r12**3 * (x1-x2) + m3/r23**3 * (x3-x2), #mass 2\n",
    "            m1/r12**3 * (y1-y2) + m3/r23**3 * (y3-y2),\n",
    "            m1/r13**3 * (x1-x3) + m2/r23**3 * (x2-x3), #mass 3\n",
    "            m1/r13**3 * (y1-y3) + m2/r23**3 * (y2-y3)\n",
    "           ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define times to solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.linspace(0, 20, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we use the special `DOP853` solver as recommended by the paper. We set very small values for `rtol` and `atol` to ensure a proper solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol = solve_ivp(dSdt, (0,20), y0=[x1_0, y1_0, x2_0, y2_0, x3_0, y3_0,\n",
    "                       vx1_0, vy1_0, vx2_0, vy2_0, vx3_0, vy3_0],\n",
    "                method = 'DOP853', t_eval=t, rtol=1e-10, atol=1e-13)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get solutions for the three particles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = sol.t\n",
    "x1 = sol.y[0]\n",
    "y1 = sol.y[1]\n",
    "x2 = sol.y[2]\n",
    "y2 = sol.y[3]\n",
    "x3 = sol.y[4]\n",
    "y3 = sol.y[5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot sample $x_1$ vs time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x26b0d0d5940>]"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABXFElEQVR4nO29d3Rk132g+d3KqEJODaDROTebuRkULJGiKJMc25Q0Y6+0PjZnxjNa7ZqzE8+O5kxY73j2HNt7vPbYa0uW1hrLlsPIY0vi0hxZEm1SgWLoJpudc0Js5FiofPeP916hUKgCquq9ewsE7ncODoCqV/Ve1b3v/vLvCiklBoPBYNi6+Op9AQaDwWCoL0YQGAwGwxbHCAKDwWDY4hhBYDAYDFscIwgMBoNhixOo9wXUQmdnp9y9e3e9L8NgMBjeU5w8eXJCStlV/Ph7UhDs3r2bEydO1PsyDAaD4T2FEOJWqceNa8hgMBi2OJ4IAiHEU0KIS0KIq0KIz5V4Xgghfst+/rQQ4oGC524KIc4IIU4JIYyabzAYDJpx7RoSQviB3wGeBAaBt4QQL0gpzxcc9jRwwP55BPi8/dvhcSnlhNtrMRgMBkP1eGERPAxclVJel1KmgD8Dni065lngD6XF60CrEKLXg3MbDAaDwSVeCILtwEDB/4P2Y5UeI4FvCyFOCiE+U+4kQojPCCFOCCFOjI+Pe3DZBoPBYABvBIEo8VhxJ7u1jvmAlPIBLPfRLwohPlTqJFLKL0opj0spj3d1rcp+MhgMBkONeCEIBoEdBf/3A8OVHiOldH6PAV/HcjUZDAaDQRNeCIK3gANCiD1CiBDwKeCFomNeAH7ezh56FJiVUo4IIWJCiCYAIUQM+Bhw1oNrMpRBSsl3z9/h//3+dYZmlup9OYY1mImn+IMf3uCbp4bIZHP1vhzDGpy8Nc3vvXqNi6Nz9b6UmnCdNSSlzAghngf+GvADX5ZSnhNCfNZ+/gvAS8AzwFUgDvwD++XbgK8LIZxr+RMp5bfcXpOhPP/55Sv85nevAPAb37nM7//9h3h0b0edr8pQzOxSmk/87mvcmFgE4M/eHODLf/8hGkL+Ol+ZoZgXTw/z/J+8A8Cvfusiv/p37+Gnj+9Y51UbC0/qCKSUL0kpD0op90kp/0/7sS/YQgA7W+gX7efvllKesB+/LqW81/65y3mtQQ0DU3F+6+UrPHtfHy//yw/T29rAZ796ksHpeL0vzVDEF169xq3JRf7oFx7mVz55N6/fmOTffcMYyxuNpVSWX3rhHPftaOX7/9vjvH9fJ//mL89w8tZUvS+tKkxl8Rbiq6/fwicE/+bpI+zrauRLP3+cTFby775xFrNT3cYhlcnx1ddv8fSxXn7sQBefengn/+Tx/fzF24O8fOFOvS/PUMBfnRlhYiHF554+zI72KL/zsw/Q19rAv/zauyQz2XpfXsUYQbCF+M75O3xgfyc9LREA9nTG+GcfPcArl8b57oWxOl+dweHNG1PMJzJ84v7lLOznP3KAfV0x/uOL599TC8xm59vnRultifDInnYAWhqC/KePH+PmZJz/8sOb9b24KjCCYIswMBXn+sQijx1amXr73Pt3s7+7kV/91kVyOWMVbARevTxGKODj/fuXYzehgI///Sfv4tZknD9+/XYdr87gkMnm+MHVCT5yuBs7zgnAhw528dEj2/h//uYqM/FUHa+wcowg2CK8MzADwMO25uIQ9Pv4p08c4OrYAv/97GgdrsxQzLuDsxzrayYaWpnL8aGDXTy8p50vfu+6sQo2ANfGF4mnsjy4q23Vc//qxw+ykMzwB6/d1H9hNWAEwRbh/PAcQb/gQHfTqueeubuXvV0xfvtvrhiroM7kcpLzw3Mc295S8vl/8pH9jM4l+IuTQ5qvzFDMmaFZAO7pXz1Wh3ua+eiRbfyXH95kIZnRfWlVYwTBFuHc8CwHupsIBVYPud8n+MXH9nNxdJ5Xr5j2HfXk9lSchWSGu/qaSz7/wf2d3Lujld/73jUjtOvMueFZGoJ+9nQ2lnz++Y/sZ3YpzZ+8UXILgA2FEQRbhOvjixzcVnrCAvzkvX10NYX5w/eIKbtZuTFp1Q3s7y49VkII/tEH93BrMs6rl43Qrie3J+Ps6oji95XqoAP37Wjl4T3t/NHrt8hucKFtBMEWIJ3NMTK7xI72aNljQgEf/+PDO3nl8jg37SImg34Gpqyajh1t5cfqqWM9dDeF3zP+583K7ak4O9e4pwCee99uBqaWeOXSxs7KM4JgCzA6myAnob+tYc3jfvaRnfiF4I9e3/im7GZlYCpOOOCjqylc9pig38fPPrKLVy+P5yuPDXqRUnJ7yrII1uJjd21jW3OYr/xoY99TRhBsAQanrZ5C/WtomQDdzRGevruXr50YYPE9EODajAxMLdHf1rAiHbEUn35kB0G/4I82+AKzWRmbT5LM5Na1CByh/b3L41wfX9B0ddVjBMEWwGkhsZ5FAPD337+L+USGF94tbiBr0MHgTHxdgQ3Q3RTh6WO9/PnJAZZSJpVUN6OzCQB6Wta/pz71sCW0v7qB6z+MINgCjM0nAdjWHFn32Ad2tnGgu5GvnRhY91iD90wupOhsLO8WKuTTD+9kPpHhW+dGFF+VoZjJReue6mwMrXtsd1OEjx3t4evvDJLKbMwuslteEEwvpvhHX3mLP31z40prt0wtpoiG/ESC63euFELwM8d38M7tGa7cmddwdZXz4ulhfvoLr+UDqpsNKSWTiyk6KlhcAB7Z087O9ihfe2tQ8ZVVRzqb4/k/eZtffvH8+ge/R5lYsCqGKxXaP328n+l4esP2itryguBL37/Ody+M8e+/cXbT9nyfXkzRFq1scQH4xAPbCfgEf35y4ywwyUyWf/v1s7x1c3rTBrPjqSypTI72WGVj5fMJfvrBfn50fZLbkxtHOH73/B1ePD3C7//gxqYV2pO2IKhUaP/YgS56miMb1tLe8oLgpTOWWZ3JSW5t0kk7FU9VvLiApeU8caSbv3x7kPQGEY5v3phidikNWIU8m5GpRWtxqWas/u6D/QgB/+3kxllg/uLtZQXCqb7dbEwuJGkI+le1ASmH3yf4ew/28+rl8Xx8YSOxpQXBjYlFbk7G+eQDVpfH6+ObMxVvejFFWxWLC8DPHN/BxEKKv7m4MfKfX74wRjjg46m7erg5sTkF9qQtCDqqGKu+1gZ+7EAXf35ycEMULWWyOX50bZJn7+sDYGh6c+6CV40Lz+HvPdhPTq4UlBuFLS0IXr8+CcAn7+8HrAVzMzIVT1W1uAB8+GAXnY0hXji1MbKH3rgxxUO729nTFWNsPqFl/4S3b08TT+lLo52yA5DVWARgLTAjswnevFH/zVDODs+xmMry5NFtNIUDWrZDHZtL8Ps/uKF1T42JhSQdFcYHHHZ3xnh4Tztff2dow+3/saUFwanbM7RFg9y7w2oaNf0eaRlbLdOL6apiBAABv4+nj/Xy8sU7da8pWEhmuDQ6xwO72miOBElnJYm0WpfVUirLJ3/3Nf7RV04oPU8heb9zrLoF5qNHumkI+nnxdP2F9hu2cvXwnnZaokHmbHeeSv71X5zml188z7lhffsFzy6laW0IVv26n7y3j6tjC1y+s7FqCra2IBiY4d4drTSGAwT9ghkNk/b7V8b5O7/1fW2LazKTZSGZoS1a/aT9iXt6SaRzvFxn99DpgRlyEh7Y2Upzg+WTnUuoHatMzhI0r12bVHqeQuYT1pxwPmOlREMBPnKkm2+dHa17wsPJW9Ps6YzR3RShMRzQ0nkzayvX4wtJ5edyWExmaAxXv+X7U3f14BPwVxtAaBeyZQXBQjLD5bF57tvRihCCloaQlk0k/vEfnuDc8Jy21gCLSavYqDFS/aR9aHc725rDvFjn4rK3b08DcP8OyyIAlGua9XC3O26oSgOQhfzE3b1MLqZ4o87uoYuj8/nOqboEQbut5Oh07cZTWWLh9dOxi+lqCvPo3g5ePDOyodxDW1YQnB6cQUqrQyBANOTXUqHp5PLPxNVbHwBLaeszRUPVT1qfT/DM3b28cnmcecUa+FqcHpxlb2eMlmiQZtscn1UsCOpxk8ZTWQI+UbJV+Ho8fribaMjPi6frV1y2kMxweyrOkV5bEEQCWixfJxFiSqMgWExmahLYAH/nnl6ujy9ycXTj1OlsWUHw7oCV1uYIgqBfkM6qv/kdv6JTmaiaJVvLrKSYrBQ/cU8vqUyurtlDF0bnOGJrmQ3251AdI6iPRZCtSWCDNb4fPbKNvz43WrfsoUv2wna4x9r8qDEcYF6LRaBXEEgpWazRIoBC99DGqQjfsoLg/Mgc/W0NtNqTKBTwk9LgX3XOp8uMjacci6A27eW+HW20x0K8cqk+ve/nEmkGppY4amuZAb/VjC2dUy0IlhdTXX73eCpDrAa/s8OTR7cxtZjilL0tqW4ujlrB2sP2WEVDfuJJ9Va2Y0HpEgTJTI5sTtZ8T3U0hjm+u33DpGbDFhYEF0fmONyzvAtUyC+09AFxXBu6Ju2yIKhNe/H7BI8d7OKVS2N10TQvjlhaZl4Q2JuAZBRbb4WCYFqTG28xlaWhxnEC+NCBLnyCuvW+vzgyT1MkQF+L1dPK7/OR1eBic86h+56KuRirxw51cX5kjjtzG6O4bEsKgmQmy/WJRY70Lu/fGwr4tFTROs2FJzVNWidG4GaBeexwN9PxNO8Oznh0VZVzYcTSMo/kBYE1ZVVr6YXrl6604qVUlliNWiZASzTIg7va+Nt6CYLROY70NOdbaPt9aNlO0xkrXXE3J+7hxnp7/FA3AK/WydIuZksKgqtjC2RzcoVFEPT7tFgEOc3ay5JLiwDgQwc6LU2zDqbshZE52qJBtjVbufVB2zWUUbzAFFoETn6/ahaTGVcCG+CxQ92cHZpjTLOmKaXk4sg8hwuUq4Ami8ARNjNLmi0CF4LgcE8TPc0RXrm8MdxDW1IQOO6GQz36LQLnvtBtxjbUGCwGK67xwM42/rYO2sv5kTmO9i1rmQG/bREojxEs/63LIoinsq7cDbCsab6ieT/joZkl5pOZFfeUTwiyGhIwnLHS5cJzUmLdKFdCCB471MX3L09siH5eW1MQjM4RDvjYXbDNXNDvI7kpLQJr0rrVND+wv5Nzw7PK0zYLyWRzXBqd50iB5ebECFRneBW6NHS58eKp2lMSHY70NtHZGOJ1jYVwsKxcFVrZfh96LAL7HDPxlJa0X6few41FAPChg13MJzMbojHfFhUE8xzc1pTXLkGfRaBbELjNGnJ4dG8HOQknbuorWLo5uUgyk8vHB8AS2KA+WLwiRqBxrNxomWBpmg/vaddeWOZkDBVaBH6fT0uCgXNPpbNSSwGbU6Tpdqwe3tMOsCF6RG1ZQVA4YQFCfp+W9NFlMzalJZCWDxa7cA0B3L+zlZDfp3WBOW9rmYWCIJCPEehLH9UptN0uLgCP7OlgaGZJ614AF0bn2dkeXdF2we9DqyAAq6+WavIWgUvlqrMxzL6uWL4/Uz3ZcoJgYiHJ+HwyX/TiEPL7SGd0ZDhY58hJ9dWxYAWLwwEfft/am6GvRyTo576drVon7fnhOYJ+wf7uxvxjQTtrSLlrqC6CIEPUpbsB4JG9lqapU2hb6dgr7ym/EJpcQ8t/T2mI5zjFjLUWaRbyyN4OTtycrnsL8S0nCJzqx0ItEyAYEFotAtDje/ZKywR4dE87Z4ZmtZjfYGUM7e9uWtFyIW8RKB6rFYuLhnFKZ3Oks9K15QZwsLuJ1miQN2/oEdqJdJYbE4v5QjIHn08gpfoU0pVCW33FvuNCrqUVSDGP7GlnPpnhvMbOqaXYcoLAyUsvqb1oMmMd5VzHAhNPZT1ZXADu39lGTsJZTcGtS6Pzq8YpoCl9VGq2CLxcXHw+wQM727RVGF+5s0BOwpHisbInumqroFDQTGlwDTlp5k4qsxuO77ast3cGpl2/lxu2nCC4NDpPZ2N41aYSQggtGQc5Sf7cOrSXZCbriQkLcE+/tW/DaQ2FZbPxNKNziVWxnGXXkB6LoDEc0CQIrBMG/d7ckvf0t3B1bEFL07cLRa0lHHyOIFBuESz/reOecjwHXoxVX0uEjliI04P1zRzacoLg4uj8iopiByFAh5dOSklnXhCo117S2Vxei3ZLR2OY7a0NvKth0pbKQgFrcRFCX4uJjsYQUxrSEvMWgUdjdU9/izbr7eLIPA1BPzvboyse99u1HznVFoGUNIUDhPw+bfcUWHFFtwghuKe/hTNGEOgjm5NcvrPa3QAgEGgwCMhJSWej0y1RvfaSycp8WwYvuHdHixaL4PKdlZ0sCwn6fNoqizsbw6QyORYVtyh3FpeAZxZBK4AWTfPi6BwHe5pWJST4NVkEUloKQlssqCXVN5XJEfCJvMXjlrv7W7kyNq91W9RiPJl1QoinhBCXhBBXhRCfK/G8EEL8lv38aSHEA5W+1kucvPRDPc2rnhNCTw/6XM7KNmgMB7QEi9M56Ykv0+Ge/lYGppaUu0sujs7THAnQ0xxZ9VzAL7T1GnL2ela9wGQ8dg115q23GU/erxxSSi6MzK2KD4A+QeDE3dpjYT33VDbnSSzH4V7betO51WYxrj+NEMIP/A7wNHAU+LQQ4mjRYU8DB+yfzwCfr+K1nrFc/VjKItDjGnImbVssqMX3nMnmPNMyYTlOoLoa8pJd6+G0lihER1rismvIcuOpXmCW/c7eCe27t7coH6fx+STT8XTJe0qXIMjmJD4haI8FtbQDSWVynglssMYJ4N06tQ8HbyyCh4GrUsrrUsoU8GfAs0XHPAv8obR4HWgVQvRW+FrPuDg6h9+3Mi/dwUl1U42UlhuqPRbWJAhkPnvDC5x2D5cV7q4kpeTSndVFf3kEysfKWbu6GvVYBGkPA5AOR/uauTUZV+pyuOBsRtO72sr268oasl1Duu6pVFZ6Ok7dzRG6m8J13bHMi0+zHRgo+H/QfqySYyp5LQBCiM8IIU4IIU6Mj9fWUKs9FuLH79pWMotGoD6ohX0On89yOWjJRsl5q720xUJ0N4W5dEfdpB2eTTCfyJR04cFyK2+V6LYInGJGL8fq4DZLkF65s+DZexZzsUw6NiwHi9XHCGzXUFSPlZ3O5gh76BoCa6wuK7yn1sOLT1Pqviwe+XLHVPJa60EpvyilPC6lPN7V1VXlJVr8gw/s4Xd/9sHST2rKGspJiRCCdk2CIJOVnmUNORzqacoX5qng8mh5Fx441pueOgInw0u5RZDz3jXkWFQqhfbF0Xl6miP5nfcK0Zc+6riGwswupZWnFluuIW/vqYPbmqx6jDpVGHshCAaBHQX/9wPDFR5TyWu1IDRJAimxJ22IyUX1aYmZnLdZQ2BP2rF5ZTe4YyI7Gm0xlvWm5NR5nPdvbggQ9AsNFoH3rqGd7VHCAZ9SN97F0ZV7EBTiuCQVt4WyXEN2jADUb1DjdbAY4OC2RpbSWQanlzx930rx4tO8BRwQQuwRQoSATwEvFB3zAvDzdvbQo8CslHKkwtdqwSdAapAEyxkOIS1piZms99rLoZ4mEukctxU1Nbs0OkdvS4QWe1vPYoQQysfKEXJ+IWiLhjTECLx3Dfl9ggPbGpVZBOlsjqtj8ytaTxefHzQ0CMxJhLDclqB+/4iMHZz2koO29VYv95DrWSelzADPA38NXAC+JqU8J4T4rBDis/ZhLwHXgavAl4D/Za3Xur2mWhBCvZYJhdqLprTEnPQ0awhgX1cMgBsTanzPpbrDFiLQESy2TiAKrDeVqHANAezvauT6+KKn7+lwfXyRdFaWLNAE8ouljoIyv0/QEbPceBPzautzpH0+LzlgJ7BcGVMXz1kL960OASnlS1iLfeFjXyj4WwK/WOlr64FVUKbHIhBiOT99cjHFjqKKTC9JZ3MEPZ60uzssQXBzwnuLIJ3NcW18gQ8fKh8HEkIoF9rOVPD7LEGgWstU4RoC2N0Z45vvDpNIe9dqxMGp/l7PIlDdy9FRrrqaLEEwvqBWEDjn85KmSJDOxhC3JtUI7fXYUpXFa6GvxYQ1iRwzVnV1sYpgcXssRFMkwE0Fk/bmhKVllgsUgzVWqkfL0WJ9tstBdWA/74pSILSlhMFp74X2hZF5gn7BXttCLEaba8hWrrrswP64YougsHGkl+zqiCm5pyrBCAIbq+mc+vM4kyhvESjeGD2T87agDKzvandHjBsT3k9aJ1B8aFtpLRN0uYbscwlBZyzEhAYtE7zXNHfZ27HeUGC9XRy12oSXs2Kcz6J6rBzlqrnB6jc0ofieyklKFjq6ZVdHlFuT+jYTKsQIAhtnWFW7h5ZT3fQEttJZ6blrCCyXgwrt5dLoPH6fYF93aS0TrJteV4zAJ6yCn/lEhkRaXWDfCX57vb7s6bS+RxUuh4sj8yVbSzgs31Oen3oFVmWxtTh3NYWVWwRSkUWwuyPGyGxC6TwrhxEENs4NqEPTFELQaHdLVB2E9LrFhMOejihD00v53uxecXF0nj2dMcKB8v5sK7Cvp47AJ4QWl4PMWwTevm9rNERLQ9BzoT29mGJ0LlE2dRTAyVpWneHlKFcAnU1hDTECqcwiAJRl462FEQQ2wtZfVHuHHG1CCEFHY0i59pJVkOEAsKM9Sk7C8Iy3ec+X7sytmTEEevpCOW7twiDkmMKxWhZs7w2Xw8V80d9aLjwna8jTU6+iMHjbpeGeslxR3r/vchKG/jiBEQQ2vrxFoKEvij1pu7WYsWpaMvS1NgAwPOudIFhMZhiYWuJwmUIyBx3xnOX0UZazURSPFahZYPpaGjwX2PmMoTUsAjTdU9Ju2wJocQ2psgj62+x7yuOxqgQjCGyccVWvvSz7F7uaIuoFAWoCW70tVnvokZmEZ+/pFNMcXM8i0FD8Vxi87c4LAu8+6+rzLdcteE1va4SR2YSnC/LFkXk6YqG826wU+WCxZ2ctTbbANdTVGGZqMam0rUVOkUXQHgsR8vsYmVM3z8phBIGNyE9a9furOufqbg4rdTeApS0pWFvyFsGIhxbBpXV6DDkIDd1H8zECn9V4zif0xAiUWG8tDcRTWeYS3nUhvTg6x+He0m3CHfQlYCzfv11NYXJS7T7TUnpfWQzWZ+hpiXiqXFWKEQRF6Ep1A8s1NLWY8jzgWnw+FYtLJOinPRZieNa7SXtxdJ5oyM+OtrUL7HQU/xVaBH67xbHKIKQsOJ/X9Lba1ptHQjubs9qErxUfAH3po5mCoslODYF9FQVlDj0tEUY9vKcqxQgCG12TttA11N1k3aAqc9Qt15Ca9+5tiTDioT/z0ug8B7Y1rbsFoI7iv8L0UbA0zbE59cFiFWOVd+N5tMDcnFwkkc6tH9TX5G4tLJrUUV2cU2Rlg31PzZkYQd3Ip49q8D07C123hmwUVWYsQG9Lg2eLi5T2ftLrBIpBbx2BY091K05LVPlxeltsN55HLocL9h4ER0tsRlOILtdQ4Z4bOgL7Ki2C3pYG7swmtbejNoLARlfxS6E20d1sCwKFwaGcItcQQF9rxLMMh/H5JJOLqXW1TNCzidDyHsLLmqZKi8CRBF5tiF5Id5MV4xj1yDV0fniOgN3ZdC2EpmBx4S58jmtIqZWt2CJIZXNMadhysxAjCGx0ZTisjBFYJrvqgLGqWbutOcJcIsOSB620zztaZt/aWiagZRMhpz+OU4zX3RRmYkGdppZ3DSl474DfR3dTxLN4zoWROfZ3N65Z9AeFriHFFkFB0WQsHCAa8iu2CNRZ2T22G093nMAIAhtdkzYnZf5m72wMIYQ6QSCL/Nxe0+Wh9uUIgiPruBvAXiwVS4L8/gC+ZYsgk5PMLKnZ9MT5OKoWmC5bkHnB+ZG5yscJlI9VNrdyX27VtQSqCsrAUq7ACIK6o1IOSClX7A8Q8PvoiKmrhMw3TlPkHOpssvoleSIIhufob2souxlNIT4NG9NkijaTX64uVnODqgwWA3Q0hjxpcDi1mOLOXHLd+AAsu7nUW28r99zoagwrGydQ13QOlptR6tjGthAjCGzEcrRYGRl7ZQ4VtIXubAwrK1SSqhcXeyMQLxaY8yNzFS0uYPcaUtzj3rEInGwUx42nSmirrCMAa6wmPRDYF6q13NDjGipsrLitJaI0nqOq6RxYAhtQ3oOsGCMIbPItJhRKgnR29eYj3c0Rda4h+7eqxaXT1pInXe6pEE9luDGxWFl8ALuOQLGeubxjmJ5slPxYKZLanY0hJjzYI/v8sCMIKgjqa2rkWLznRk+z95XUhaiMEURDVozDC6FdDUYQ2CxrL+rOkc6s3pe2W2E2Sr5ISZH64pixbvu/XxqdR8rKtEzQU1mcKdpDeJud4TWqKMNLufXWaO2RvZB0V118YWSObc1hOtZoLeGgLWuoaM+N3pYIS+ksc0veVVIXojJ9FNCyNWoxRhDY5CetwhUmlV29L63KbBTVJnkk6KcxHHAdIzhfYV66g7V5vVrS2RxCLO+yFQ0FaI4ElAXxdLiGwL0bryoXnv1bvWto5Z4b+YCrIqGtsqAMrJYmRhDUCQ0hgnxKYrFFkMlJpQOvctK2e7CN4/nhOZoigXz3xfUQaChSykqCvpW3h5cFdMUU7n+ggmXfc+1CO5nJcnVsoQrLTcNNxeo9N5YrqdVU6EqFwWKAzljIuIbqxbJFoO4cpVxDPXbVpwpNc1nLVDdpWxqCzLlMqXS0zEpvLj2uodyqvZ5V9oFZ3hpTydvTGrUEwUy89rG6OrZAJicrFgSOkq7cIsitjBE4FsEdhRaBqmAx6NkjuxgjCGx0lMPnXUOB5a+9z+OGYIUU98tRQUtDkFkXgiCbk1wana84UAx6eg2ls7lVe/H22e2cVaA6WOyk5c4lah8rJ1BcTVAfNAntEq4hddab2hiBF8pVtRhBYKPVNVQwafN9YFRYBPZvla6h5oaAK0Fwa3KReCpbsZYJTq8h9VpmsNgiaG5gYiGppFus6mBxcyQA4CqAemFknkjQl99Jaz103FO5nCQnIVDgxgsFfHQ2hpVaBErvqUiQxVQ2X8uiAyMIbHRoL6VcQx32ZhRe7vTlkF9cFLuGZl0sLtUGisHpNVTzKSvC0jKLYwTqXA6qg8XNtkXgRmifH5nlcE9zxVuf6qjWd2pzVgntlvB71iJosoX2vIf7R6yHEQQ2OuoISrmGfD7Btpawks0o9FgElhlbq4ZeaQOzFWjJGpIEA6tjBKAmG8WZd6oWmKDfRzTkr9nlIKXkwsh8VZabDuUqmbH6XIUCK5eynuYGhfEctTECR2gbQVAHdPROz5RIHwXLPaQkWGxbliozHFoagqSyOZI1ukvOV9jArBA9WUO5EllD6nzPqoPFYLkcarUIhmcTzC6lOVpBIZnD8tenULmy512oKJ7T0xJWmj6q0iLIu/FcxHOqxQgCm2XtRWVl8WrXEEBfS0SNawh1HS0dWly6HC5UkZfu4NOQNZTM5FZrmfnOkCrceNZvpQtMQ6DmxeXM4CwAx7a3VPwa555SqVw5VnaoSJHobWlgJp4mkXbfGbcYlb2GYNki0BkwNoLAQUM5fKkWEwC9rQ3cmUt4XlQmNWiZbgTB+HzSamBWRcYQOAVlaiVBMpMjHFy5uDRFgjSGA4osAvUbkVjZKLW5G84OzeL3iepcQxruqbxFUCS0VXbxVLkfASzHCIxFUAdUamIOjj8zHFjtckhnJRMue/YUo7q1MVjuBqhNezk7ZGmZd1ehZYLjGqr6dFWRTGdXjROo31N2o7qGTg/NcqC7kUiwcheelrhbGUGg2o2nNEaQv6dMjEA7OsrhE2lr0hbfTF5vJeigurUxWBuBACzWsDnN6cFZhIC7qhQEOraqTGZyJQVBb4uaWgLVlcVgB/Zr0DKllJwdmq1aYKPBNZQsEyNQWVSmcvtXKHANGYtAPzrM2LUsAvC+qEx1SiJALGwJtcUampmdGZplb2eMRluYVIxQ70pJZnIltd+eZjUWQU7DWEVD/pp2kxueTTC1mOKe/iott/w9pT4Tr/ie6lFuEShMHw07riFjEWhHx1aVjvZSvMD0tVoWwbDHFkE+WKxw0sZCtkVQkyCY4Z7+1qpfJ1BfWZzMlHYN9bZEGJtPeF7ssxzPUTdW0ZCfxVQN41RDoBgK7ikNMYLisWoMB2iKBJQE9lUXlPl8gkjQx1INY1XzObWdaYOjo/jFyWAonrRt0SDhgM/zdDcdweK8a6hKQTA2l+DOXLLqxQXsz6M8RpArmdLa29pATnq/vehyHYGnb7uCaChAIp0jW6WvppZAMRS0balDjABge2sDQzNqMrxUxxSjoQBxD/YCrxQjCIpQ6hpKl9ZehBD0tkQY9njS6mg6Fw3ZrqEqJ+0ZO1BcrbsB9GxVaWUNlQ4Wg/duvJwGi8Bx4y1VmVJ5poZAMei1CEoJgv62Bgan1VgEKgU2QEOwNjderRhBYLN8Ayq0CDJZAj6xomWug4oWx8uuIU/fdgXhgA+/T1RtEZwZmsUnqmst4SCE+hYT5VxD/bYbz/MFRrG7ASwtEyBexVjVHijWU6S5XEdQShBEGZxe8jxGobqgDCwFy1gEdcCnI1icLp2JAnY2iiKLQKX2IoSgIejPZ0RVypnBWfZ1NeZdS1WdE/VN56ysodUa8PY2NYIgJ9UGiqE26214NsHkYoq7a7DcHFSOleNujZQYq/62BhaSGc/TMHMS5YMVDfmJKyiGK4crQSCEaBdCfEcIccX+3VbmuKeEEJeEEFeFEJ8rePyXhBBDQohT9s8zbq7HDTqqIMtlooC1wIzOJfJFZ16QTx9VPGvDAV8+I6pSzgzN1ry4qG5DnctJUmXSR6OhAB2xkOeCQCKVuoXAcjcAVVXb1hoohuUtUlWOlePmagiVFgQAA9NxT8+pOn0UrM/zXgoWfw54WUp5AHjZ/n8FQgg/8DvA08BR4NNCiKMFh/yGlPI+++cll9dTM8stc9VqL+Usgh1tUXLS21qCvCKmWHuJVGkR3JlLMDafrMndAHZlsQZ3Q6kYATi+Z68XF7WWGyx/nmr6QjmB4ppcePZvtRaBnYlXynprjQLeW286xioaClQdy3GDW0HwLPAV+++vAB8vcczDwFUp5XUpZQr4M/t1G4rlSavuHKXaFjj0t6vRXkB9hkO1FsHpwdoDxaC+6dxyUL/cWEUVuYZUj5P1earZT6HWQDHoqc3JZ+KVENqOReB15pCOGEHI71Oy70U53AqCbVLKEQD7d3eJY7YDAwX/D9qPOTwvhDgthPhyOdcSgBDiM0KIE0KIE+Pj4y4vu+T7A+oLytayCAAGprwTBMuuIbWEq7QIlgPFG9M1VK7wz6G/rYGh6SVPe0NJ1AcJnM9TqdB2EygGPbU5iXQWIUqPVWs0SCzk99x6U910Dqzg94YSBEKI7wohzpb4qVSrL/WNOXPj88A+4D5gBPj1cm8ipfyilPK4lPJ4V1dXhaeuHD11BOUtgt6WCH6f8NQi0FFHANVbBGcGZzjQ3VTSr1sJqnsNJcsUKTn0t0VJZXPe1hJocDc4mTXJCoX2iMtAsZ62LVkiAX/JhVkIkc8c8orlViCevWVJdAuCdVM2pJQfLfecEOKOEKJXSjkihOgFxkocNgjsKPi/Hxi23/tOwXt9CXix0gv3GvUt59a2CAJ+H32tEQamPJy09m/VgiAS9FW8uEgpOTM0x2OHahfmqusIHKFWzh3Sn88ciufrCtySk1Kba6jSGMGpgRmAmqq/AS0dfZfS2TUViu229eYVjhGo3DUU8OVjVTpw6xp6AXjO/vs54JsljnkLOCCE2COECAGfsl+HLTwcPgGcdXk9NaPDNZRIl88aAss95K1FoL6RGVgLTKUWwchsgomF2gPFYNcRKLxHEmUK/xwcN563mqaGYLH9eVLZysbq3YEZQn4fR6rYjKaQ5YIytVZ2pMw4gfeB/Zwui8Dvq3mzp1pwKwh+BXhSCHEFeNL+HyFEnxDiJQApZQZ4Hvhr4ALwNSnlOfv1vyaEOCOEOA08Dvxzl9dTMzpa5pbraOmwoy3qqUWguujKIRL0VRwjcLTM+3a0ujij2q0q866hdSwCb+M56v3O+ayhCsfqnYEZjvY1V7V7XCHqSzQtiyCyhkXQ39bAXCLjWSfP5Y6+6hMwNpRraC2klJPAEyUeHwaeKfj/JWBVaqiU8ufcnN9LdFRBJjPZNS2C/rYGJhaSLKXWNncrR8+kDfh9ZCpU0U8NzBAK+KruW1OIEKqzhux9cEtUgIPlMupsDHtrESDVB/WrcA1lsjnODM7yPzy0Y91jy6GjxUTSjhGUo9+23oaml2juDbo+n664m+MasjbBUe+4NpXFNjq2qlyrshhgR7vjcvBG09RRWWy9f+V5/acGZrirr7lkS4DKz1fzSyvCyd+OriGMd7Q3MDjjbWBfx+IClRWUXRlbYCmddWW56UjAWDdG4HFLEB1bioKlhEgJGU1mvREEDnnXkDqW0lkiZYqUwFpcwLtaguUe92onrU9AtoKb3dEy3bmFrM+jcnFxerysJQj6PXbj6dD8ArYErWRx8cKFt6xc1fwW62LF3dZXrm575MbTFiNw4jma3ENGENjoMGPjqUy+f38plmsJPNJeNDSdA/CLyhbmy3fca5nguIZcvcWaOF0f19I0+9saGJ5ZqrqlczkkGsbJV7nV++7ADK3RILs6ojWfT9iri1LlKrW2a6gtGqQpEuD25KIn58tpSsAwgqBOqC6Hz+YkiXRuzcWlqylMOODz3DWk2sMohKgoi8ebQLH6gjJn85boOkI7k5OetQ7X0ePeef9KZNepgRnu7W91ZaVoaTGRWTtYLIRgd0eMm5MeW9m6BIGmFFIjCGwcbckrDa+YSvzOVgFMg3cWgaZJ66sweHtqYJr2WIid7bVrmeD0Gqqva2h3p/UZbnm2wKgPFjvujPXm+GIyw+U78x4IbA2uoXUsAoBdHVFueWQR6Coo8+eFtokRaCUvCBR98fEKtEyAne1Rbnpsxqp2OfiEqOh7s7TMFteCSXVl8VKqfNsCh90dMQDPxspyDakdqOWFee0v78zQLDnp3nLTkZKdyORoCK29jO3uiDE4veRJZ19dBWXVWG+enE/PaTY+qi2CeHJ9LRNgd2eMW5NxTzVe5ZqmT6w7YecTaa6MLXDfjrLtpCpGCLV1BPFUllgosObC3NMcIRzwcXPCO01TQ5Yg/grGynHh3etBUB/ULmbxVCbfXrscuzq8c+PpUq7yGVcma0gvygVB3t2wtkWwtzPGUjrLnTn3fWx0pbpV4ho6MziLlHDfzlYt53PDUjqzbh2Hz+et71lq2JgGrO9uPXfDqdsz7OqI0h4LuTqX6u6jTtxtvc2Ndnc61pv7sdJVUOasR8Y1pBnHJ6cuRuC4hta3CABueKBpanUNrfO9veMEimvtW1OAQK2WuZjMrjtOYGmanrmGNASLwVrA1nPjvTtoBYrdn8v6rco15Lhb18rEA/KZT57ECTTW5oBxDWlHtUWwWKlryEPfs66mc5aWufYxpwZm2NsZoyXqvrrTckWpDRav524A2NMZ4/Zk3JM5k9PkGvKtk3p7Zy7ByGzCdXwA1NcROPfUehZBV2OYaMjPzQkvLALrt2qh7by9qvWoGCMIbDaKa6ivtYFQwOeJReC4T5QXlK2zMEspOTUw48niAlZhlMobZCmdqcgi2N0ZI5XNMTLr3vcs0WMR+IVY0+/8zu0ZwDsXHqhz4zlpvrHw2mMlhGBXR8wTi0BXQVk1NR9eYASBTUBb1tDak9bvE+xqj3ojCOzfOlxDa31tw7MJxueTniwuAH6fT2npfTyVXVdgQ4H15ommqemGF2sHi98dnCHor21rymKEYvfGcgJGJWPljRtPV4zAuIbqhE+XRbCO9gKWpulFNorUNmnX/t5OOVrme8UiSFUWI3BqCTxx40nwabgbxTrB4lO3ZzjS21zT1pSrzmX/ViXjFpJOjKCSeE6MgSn3leA6EzDAuIa0ozxYXKFrCCzf862puOvUMV2Vxeu5hk4NTBMK+Djc417LBMtqyiisuIxXKAi2NUWIBL1JIc3mZH4OqmStscrmJKcHvXPh5f3ciq3s9WIEYFkEqWzOdQqpLtfQskVgBIFWHJ+cKpeD48+sNAiZyuQYdul73iiuIS86jhai2iJYSGYqWlyWU0g9EARS5q1SlazVF+rq2AKLqawnGUNgWaJBvzqhvZhygsWVWQTgPhtvucWEq7dZFx29z1acT89pNj75vF1FC8xCwgpA+iu42Xd7NGl1mrHlFpdMNseZIfcdRwvx+4XSGMFCIkNTpLLspl0d3sRzcposArFGjODkrWkAHtzlvujPIaAwnrOYrKxaH2Bft3VPXRtfcHVOXU3nHDehKmtq1fm0nOU9gOpg8UIyQ1Oksn2A9jgFMK61FydrSC1rtZi4MDJPIp3j/p1eLi7qLIJEOksqm6t4rPZ2NXJ7Ku5a683mZEVKglvWKsY7eWuajljIVcfRYoJ+4Ulrh1I4gqAS662rMUxzJOBaEOjr32VcQ3VBdbB4PpGhsYIJC7Ct2cp7vu6RRaBaEjiuoVILzMlbUwAc91DLdLKGVKTWOQHISgXB/q5G0lnJLZf97nNSakkfXav4753b09y/s83TRS7o9ykTBJU0B3QQQrCvu5GrY24Fgd4YgUkf1UxAsSCYS6QrdjcIIdjX5cGkRZMZu4Y/8+TtGXpbIvTZO0V5gV9hat18whIElQrt/d2NAK7HSpdFUK7X0NRiiusTi566hQACfkEmq841FAr4CJbZUrSY/V2NXBv3Jkag657S1IXaCAIHHRZBpVomwAFPtBfrt672xqXM2LdvTfOAgsUFqHif5GpYSDgWQWVCe59XgkCiJVhcLn30bQXxAbBiBGlFgmAukaG5intqX3cj4/NJZpdq38heX9bQyvOpxggCG9UWQTUxAoD92xoZmU0wn6h90jqfxVk4VeErE18ZmV1iaGaJBz2MD0BhYN/TtwXIf9+VjlVjOEBvS4RrLgWBFSx29RYVUS7D6+3b0wR8gnv6Wzw9XyjgUyKwwRqr5goFNlgWAbgLGGsrKDNN5+qDY4qpynCYT6RpClc+aQ90NwHuNE1HENTLNfT2rRkAju/2WstUZxHMJ6tzDYHlHrryHnENlSv+O3lrmrv6vCkkKyTgUxcsnqvSynasNzdCW18mnjplp+T59Jxm4xNQnD5ai2sIcLXA5C0CxSWr5czYk7emiQR9HPGgXUEhKvtCOTGCajTNfV2NXBtfcDV3srqCxSUKytLZHKcHZz3N7HII+NW5huYTaZobKh+nHW0NhPw+rnphEdT8DpVhXEN1QmVBWTYniaeyNFYhCHa0RwkFfK4sAuezqNY0lwNbRYLg9jT39rdWHMyrlIDCsVqwXUPVjNX+7kbiqSwjc4maz5vTZhGsdg1dHJlnKZ31PD4AKC0om1tKV6VcBfw+9nTGuDZWe8A4HyxWvHIa11CdEEJUtGlHLVQbgARr8d7X1ciVO/M1n1dXjCC/m1LBV5dIZzk3NKtkcfHbd6FKi6Ba1xC4dONJfa6h1ZableKrYqwCPnXFf/OJTFWWG1iFZW5iBPr6d5nK4rrhVzRp56oMQDoccOl7dnzoqheYUi1z3x2YIZOTyhYXUGQRJDOEA76q2mF4IQhyufrVEbx9e4aeZm9TfB2Cfh+pjKoYQXUWAVgB41uTiyQz2ZrOqS991PpdOFZziTS/+8pVBlzWrJQ8n+fv+B4m4PMpMWPzRUpVaJlgCYLB6aV8c61qWY4R6G+Ze/K2lY6owu+cjxEo8D3PVdFewqEjFqI1GuTqmAvrTZtFsLqO4OStaSUCGyxBoEJgpzI5Eulc1RbBgW1N5CQ1u4d0F5QVWm8nb03za9+6xMC0EQRKCQXUaC/zNbiGwJq0UPuk1RcjsH4Xai9v35pmb1fM9b63pVBaR1Blmi9YboID3Y1cvuMmsK9nYxqfb6XldmcuwdDMkue1Hg4BRTGCatN8HQ73WPfUpTtzNZ1Xd0FZoQw9cXOKgE942rcrfz7P3/E9TNDvI6VAy6x10h7YZrkcLtcYJ8jqEgRFriEppaVlKrAGQG1gf74GdwPA4Z5mLo3O19wSwAoW1/TSqvAVdR9944b3LUAKUVVQls/uqiJrCKy9PkJ+HxdHa7untO0D7lt5PoC3bk5z1/aWiprsVX0+z9/xPUxYsUVQTSYKwK72KEG/qDlOoE0QFGkv18YXmY6nPa8fcHCykFSM1UIVPaEKOdzbxEIyw+B0ba3DdbmGrM3rl/9/4/okjeEAd/V5m+LrEAqoqSNYjrtVJwiCfh/7uxu5OOJOEKi2CPxFrqFkJsu7AzM8pEhgG0FQQCjgI6XCjK2ykZlDwO9jb2ftmUP1qiN4/fokAI/s6VByPieQq2KsZpfStFSpZQL5TXdq1jQ1BYv9Rd1H37gxxYO72ggoMkdUtaFerveoxXpr4lKN46SrbUvxNp+nB2dJZnIc392u5HxGEBQQ8vtI1ZhNsBaz8RRAbQtMbxMXRmrzZ+qKEYiiOoI3bkyxrTnsaTvjQsIKLYKZpTSt0erH6ZDte75Y41jpDRZb4zSxkOTq2AKP7FWzuIAVI1BiEdj9gqp1DYE1VqNzCWbs+7Ia8haBpribU6T4+rVJhIBHFY2VEQQFBANCiT9zOp4mGvITDlRfvn+kt5nh2dombdYOpqrOGvIX5DxLKXn9+iSP7u1Qlmudtwg8FgRSSmbiKVoaqg9wN4YD7GyP1mwRaNuqUoh824I37fiAKssNIOjzKek+upyAUYsbr3brbbnFRNUvrQp/UUHZj65Pcrinmdao98kXYATBCkKKcp5n4mnaahzAo/akvVCDT1OXReBk8aSyOa5PLDI+n+TRveoWF1WCIJ7Kks5K2mqwCMByOVwYrTEbJadnq0ohlpsDvnljioag3/NGc4UEA0JJdletMQJYzhyqxXrT1nSuwDWUzGQ5eWua9ym8p4wgKEBV+qilZda2uDh9es7XMGmdPHvVFkHMzmJYSmV547qjZapzN6iKEczY7oZaXENgaZo3JxZZSlXvXsxKPRaB3yfyMYLXr0/y4K42z1uAFBLwqVOu/D5RdW0OQHdTmNZokEs1xN50NZ0TBa6hd27PkMzkeN8+Iwi0EAr4SSrwZ84spWmL1ba4dDWF6WwM1xQn0GURRO3NwxdTGX54bYLupnB+u00VhOyFq9bq0HJMLzqxnNqstyM9VrHSlRoKy7I5PfsROJXFEwtJLo7OK/M5OwQV7S89FU/R2hCs6TsTQljWmxsrW3XWUIFr6EfXJvEJeFihcuVKEAgh2oUQ3xFCXLF/l8xtEkJ8WQgxJoQ4W8vrdRHyC9IKtJfpeIrWGhcXgKN9zZwfrsEiyEl8Qr0Z61gEs0tpvn95nA8f7FJ6zrDdKtlrTdPZsKRm15Dje65hgclJPXUEjeEA84kMr14aB+CxQ91Kzxfwq4kRzMRTtLkoVjzc08zlO/NVd4x15lw1LUhqodA19L0r49y9vaVmr0JF53P5+s8BL0spDwAv2/+X4g+Ap1y8Xguq0kdn47Vlojgc6W3i6thC1QtfVkrlqaMAMdsi+MGVCeYSGT5yWO3iElKUNTQTd1xDtS0wO9ujNAT9NcUJdAWL22JBpuNpXrk8TldTOB+DUkXQJ0jncp7vvTu1mKpZYIOVORRPZatu1+BkQAUVN3JssPdhvjm5yKmBGZ44sk3p+dyuEs8CX7H//grw8VIHSSm/B0zV+npdqAgWSylrTkl0ONrbTCqbq7proq7NTmK2n/Ybp4YI+AQfONCp9HyONpb0eKym7cysWsfK7xMc6mmqySJIZXLKtUyw3F4TC0leuTTGYwe7lLujgn4fUnrfKdZNAgYsJ2Gcq9LS1mURNIUDNAT9/Ombt5ESPrrBBcE2KeUIgP27WlWw4tcLIT4jhDghhDgxPj5e8wWvhYpg8VwiQzYnPZm01cYJMlk9gsCpxJ1PZHjiSHfVjcCqJawoWOy4htyY4Ed6mzk/Mle1BpzK5mpKL64WR4ueT2T4+P3blZ/PKVTzOk5gWQS131OHepoI+ARnhmarep0TQ1QtCIQQdDeHmU9k2N0R5Uhvk9LzrftphBDfFUKcLfHzrNIrK0JK+UUp5XEp5fGuri4l57B6DXm8uLh0NwDs6YwRDviqjhOkslktWmZTJJhvxfzc+3YrP58611CKhqDf1XaNd29vYXYpzcBU5a0mMtkc2ZzUMlZOl9HOxrDSdEQHx4XiZVGZVe+RdhUjiAT9HNzWxNkqBYEz58J+9UL7WJ+V1vuZD+1THudbN/dKSvnRcs8JIe4IIXqllCNCiF5grMrzu329p4QCPs+DxXl3gwstM+D3cainqeoU0mQ6R0TD4gLwpZ8/zsRCkocUlcAX4vMJAj7huSCYdhnLAUsQAJwZmmVnhZXVjvIR1jBWx3e38/vPHefY9hYtWUr5vSM8DBjHU1lS2ZyrGAFYY/XX50eRUla80OpyDQH88seP8ZP39vKxoz3Kz+X207wAPGf//RzwTc2v95RQwOd5+qiTm15r+qjDXX3NnBuuzuWQyOTyGTaq2dMZ0yIEHEIBn+cxAjf1Hg4HexoJ+X2cHpqp+DXJtL7FBeCJI9vY1hzRcq6g/Zm8tAim7DRfNxYBwLH+Fmbi6aoaBeoUBO2xEE8d69WTVuzy9b8CPCmEuAI8af+PEKJPCPGSc5AQ4k+BHwGHhBCDQohfWOv19SJsB4u9zHBwm5vucE9/K7NLaW5NVp7lkExntWiZ9UBFPGdyMUVnY9jVe4QDfg73VudyWLYI9AhtnSzXfHg3Vk52l5sYASxbb+eGqxmrLH6f0BJ704mrxtZSykngiRKPDwPPFPz/6WpeXy9Cee1FEgp4M9ATC0kAulwuME4bgHcHZ9hdYbFWMpPbtIJARcvwyYUUO3e6b5R3bHsLL747XLHLwbEINuNYOWmQibR3xX9Ttru13aWVfbinCb8dMH7qWG9Fr0llcspTR+vB5pt5LnBK7b00YycWUgT9guYGd5tJHNzWRCTo492ByrWXRDqrzTWkGxU1H5MLSTpi7gQ2wD3bW5hLZLhd4d6yqay1SOpyDemkIegIAi8tAifN151FEAn6OdDdyJmhymNv6azMWzmbic33iVygopmZs7i4jfoH/T7u6mvh9OBMxa/ZzBaB1zUfS6ksi6ksnU3uuzses10OpwcrE9qJTWwROBlYSx5aBOPzlpXt1o0Hlnvo7NBsxe7gZCZHaBO68DbfzHOBimZmEwtJOhq9aR17T38LZ4dnK94D1hIEm2/Sgt0XykNB4LjwOj2wCA5uayIU8FWco57SlJteD5QIgoUkoYCvpk1pirm7v4WpxRTDs4mKjl9KZYiGNt89tflmngtU5Kd7EYB0uLe/lUQ6V/HWlcl0lkhwcw6x166hSTuo74XQDgV8HOlt5kylFoHdrdRN/cJGxZl/XsYIxueTdDW6t7Jh2Xo7U6GlvZjKGkGw2XH86V52tZyY984iuHdHKwDvDsxUdPxSOrspFxewMrySHi4uk7ZF0OGR0L57ezNnh2Yramq2aAuCWvZK3ugsxwg8FgRN3ozT0d5mQn4f79yeqej4JSMINj9Re9LGa+gnXwopJROLKdcZQw67O6I0RwK8W6GmuZCsbSP29wLhoLd1BJMLlkXQ6ZHQvm9HG/PJDFcr6A8VT1m7bW3GBUZF1pCXgiAS9HPX9mbevj1d0fHxVIZoaPPdU0YQFODciF4JgvlkhlQm55lFIITgnv7WigLGUkriqWy+M+hmoyHo93RxmVi0LQIPYgQAx+1WDidurr/ALCQtQbAZhXbEjlHVsllPOSYWvBMEAA/sbOPdwdmKXMJxYxFsfhztxatJO+FhdoPDvTtauDg6n9ciy5HMWP1rNqP2ApbQ9kpgA0zMp4iF/Pk54JZdHVE6YiFO3lpfEMST1ueIbkJBkL+nPEofzWRzTHpoZYPVfymVyVXUwsUIgi2As2h6tcBMLDgBSO8m7fHd7WTt7evWYjNrmQANocC6wrAaxheSdHqoZQoheHBXGydvleq+vpJF+3M0bMJ4jpMS65X1NrWYQko8twgA3q5EaKcym1JgG0FQwLJryJsF5s6clZK2rdlb7UUIa+PxtVi0BUFsE05agJjHFsGd2YTn/Xce3NXGzcl4PjW1HAuJDA1B/6ZrWwCWQIwEfZ4JgjHbyvZSEPS0ROhriXBynTiB0/XUTQPJjYoRBAUsm7HeTFpHEPQ2N3jyfgDNkSBHepp56+bagmA+sbktgmjIz1I6W/VWg+UYnUvQo0AQAOu6h6bjadpdNlDbyDQE/Z7dU+MKBAHAA7vaeGedcVpIZsi43Ftko2IEQQFeB4tHZxNEgj7X7SWKeXhPO+/cnlmzFYaXefEbkYZQACkh4UGqr5SS0bkEvS3eCoJj21sI+X3ruhxm4inX7a83MhEPA/t5QeChuxUs99DwbIKR2fKdSJ1mdy2bcKyMICjAyXDwShCM2Fqm15tKPLS7naV0ds0Ol07X082qaTrZUF6M1XQ8TSqT89w1FAn6ObZ9fettKp7atOMEjkXgTbB4fEGNRfDIXquF+uvXJ8se41XX042IEQQF+HzCmrRexQhmE/R4rGUCPLTHcjmstcA4FkH7Jpy0sBxY9SLDa9RuL6BirB7Z28Hpwdl88L4U0y63XdzohD22CJrCAc8LJY/0NNMaDfLa1fKCwO2e1hsZIwiK8DItUYXfGaC7KcLujihv3ijvcpheTOH3CdcbrWxUnAyvRQ+E9nJQ3/ux+uD+TjI5yZs3yi8w1v67m3OcwLmnvFGuhmeW6G31fpx8PsH79nbw2rXJsg3o8ptMbcKxMoKgiIaQ3xMtM5eT3JlL0NPiXaC4kPft6+CN65Nl4wST9uKiY3ejehD10DU0OqfOInhwVxvhgI8fltE0M9kcc4mM6922NjKN4QALSY/crbMJ+lrV3FPv39fB0MxS2fbhTvtrt5tMbUSMICjCK4tgKp4inZX0eJg6WsiHD3Yxn8yUrSeYWkxuandDvh2IBwvM6GwCIaDbY78zWHGCh3a388OrEyWfd7TMzRwjaIwEWEikPXmv4ZklepUpV50AvHattNCeXrQ+g3ENbQEaQgHiHvgzl/3OirSX/Z34fYJXLo2VfH5qcXMHIJeL/9y7HEZnE3TEwvmNibzm/fs7uDg6z9j86lbHTo3BZh6rpnBgzRhJpSTSWSYXU2xX4BoC2NcVo7spzA/KCO3RuSU6G0PK5kk92XyfyCWxkD9fjOUGlQFIsOoJHtzZxquXx0s+PzyjzoTeCDiuIS/y00fnEvS0qLHcwIoTAPzgyuoFZmDKSlfc0eZ+i8yNSmM4wELC/T01PGN9V6rmtRCCxw518b3L4yX7Dg1OL7F9k46TEQRFNEUCzHtgxub9zgoCkA4fPtTFueG5VZpmMpNleHaJne2bc9ICxGyLwAtN846ioL7Dsb4WtjWH+fa5O6uec/zROzbxWDVGAiymsmRdFv8Nz1jzXKWC8+TRHuYTmZKV+wNTcfrbNqdyZQRBEU2RYL4q1w3DM0sEfMLzfOdCHjvUBcDLF1a6hwamlpASdndu3sXFKdJzO1ZSSvsGV/dd+XyCjx3t4ZXLY6sSEQam4jSGA5syE8XBqW53m+HlWATbFQqCD+7vJBzw8d0LK4V2LicZmlnatJabEQRFWBaBe0FweyrO9rYGpf1jjvY2s7czxjfeGSo69yIAO9tjys5dbxqCfgI+wdySO+ttOp5mMZVVrpE/dayHRDrH966sdOXdnoqzoz3qedHhRsIRBG7dQ0MzSwihJs3XoSHk58cOdPLtc6Mr2peMzCVIZ6WxCLYKzZEgC8mMazN2YFq9a0YIwbP3beeNG1N5bQng0qi1Gcrezs0rCIQQNDcEmXPpxhtwXDOKb/CH97TTFg3ywqnhFY+fH55jf3ej0nPXm8aIN2684ZklupvCyvd2/sl7+xieTfB6Qe2HU8V/pLdZ6bnrhREERTRFvNFeBmxNTzXP3tcHwNcLrIKTt6bZ2xnb1LnpYI3V3JK7cdLlow/6fXzygX6+fX40nyk0MrvE6FyCB3a2Kj13vXEsAreW9shsQlnqaCE/flcPTeEA/+3EYP6xdwdmCPgEd/UZQbAlaI5Yvlo3muZCMsPUYkqLP3F3Z4wP7O/gK6/dJJG2AnJv357Od77czDRHgq4D+wPT+oK1n354B+ms5Kuv3wKWW4nfv3Nzj5WjXLnNxhuYttytqokE/Xzige38f6eHGbIt7TdvTHG4t2nT7gFuBEERXgQhHXeDrqydX3xsP2PzSb76+i1evTzG1GKKxw93azl3PWluCDDn2nJboj0W0tKue393E0/d1cMXv3edsbkEX39niJ7mCHdvb1F+7nrSGHavXKWzOQanl9jTocfd+dkP7wPgN79zmYGpOCduTfP0sV4t564Hm7NZvQuaPLAIlt0NegJL79vXwROHu/m1b12iLRZkW3OYJ49u03LuetIcCTI2t/7m8Guhy4Xn8K+fPszfXBrjk59/jcHpJf7pEwc25YY0hTgZUdPx2u+pweklsjnJbk1xr77WBn7hg3v5wqvXeO3aJCG/j0/cv13LueuBsQiKcMzY95JFIITg13/mXj5yuJumSJD//Kn7N2X1YzHNEQ+CxdNx5YHiQvZ0xvi9n3uQcMDHM3f35DXPzUyr3erEaY1eCzcnrEy43R36hPa//NhBnnvfLgJ+wa/83bs3dYGmsQiKyMcIXKQlDkzFaYoEtHb+bI2G+MLPPajtfBsBt6m+2ZxkaHqJv3O3XpP/8UPdPH5o87vuHEIBH03hQL6Ncy3cnLQFgcZMuKDfx//x7DFt56snm19trJJ2e0evKTfay2ScnZs8N3wj0BoNEk9lSda4S9ngdJxMTrJLo5a5VWmNBV1bBI3hAB2bPBOuXhhBUERTOEDQL/Ibu9TC1bEF9nVt7tzwjUB7zKrarlVoXx2z4gubPY9/I9AeDTHlIkZwczLO7k6jXKnCCIIihBB0xMJM2rne1RJPZRiaWTKLiwac/ZgnF2oTBNfGbUHQ1eTZNRlK0xYL5fv518LNyUV2a8oY2ooYQVCC9lioZi3z+rjlyzSCQD2dtiCYqFFoXx1boLMxvCk3I99otEVrv6dSGSt11AgCdRhBUIKOxhATxt2w4enwwDW0r8ssLjpoi4ZqjhHcmFgkm5Mc2GbuKVUYQVCCjliIqcXatUy/TxjtRQNuXENSSq6NLxqBrYm2aJDFGgP7l+7MA3Cox7jwVGEEQQk6GsNM1eh3vjq2wK72qPLGWAarh03I72OiBqE9sZBidiltBIEmnL5XMzUEjC+NzhHwCfZ2mrFShavVSgjRLoT4jhDiiv27ZNMUIcSXhRBjQoizRY//khBiSAhxyv55xs31eEV7LMRiKkuiht2vro4vsM8sLloQQtDRGKrJInBceCa7Sw/OVpy1uPEujS6wpzNmlCuFuP1mPwe8LKU8ALxs/1+KPwCeKvPcb0gp77N/XnJ5PZ7g5CpXm0KaSGe5MbHIoW3GhNVFR2NtQcjzI3PA5m0rvNFwBEEtQvvSnTnjFlKMW0HwLPAV+++vAB8vdZCU8nvA6r3fNigdjVYQstoU0ouj82RzkmObvInYRqLWVN9zQ7N0N4WV7iBnWMbZCnRkdmmdI1eymMwwMLVklCvFuBUE26SUIwD271rq5p8XQpy23Udl+/EKIT4jhDghhDgxPl56w3avcNISx+erW2CczSuObTdapi46GkNM1KBlnhue27S95TciPS2WIBidTaxz5Eou24Hig8YiUMq6gkAI8V0hxNkSP896cP7PA/uA+4AR4NfLHSil/KKU8riU8nhXV5cHpy6P01xquMpJe254ltZoUOmeqoaVbGuOMDafWLGt4Hok0lmuji8Yy00jkaCf9liIkbnq7qnTg5ZytdlbddebdZvOSSk/Wu45IcQdIUSvlHJECNELjJU7tsx753eIFkJ8CXixmteroqsxTNAvGJquzox1tExTBq+P7a0NpLOSsflkXutcD8eFZywCvfQ0R6q2CN4dmKGrKUxvhWNrqA23rqEXgOfsv58DvlnNi23h4fAJ4Gy5Y3Xi8wl6WxpW7AO8Hulsjosj8xzrM5qLTpwdq4Zm4hW/5tywpWXeZcZKK32tkaruKYBTgzPc299qlCvFuBUEvwI8KYS4Ajxp/48Qok8Ikc8AEkL8KfAj4JAQYlAI8Qv2U78mhDgjhDgNPA78c5fX4xl9rZH8NnWVcHFknlQ2Z9wNmum33XCDVVhv79yeoT0Wol/jPgQGK04wWoVraHYpzfXxRe7bYe4p1bjaj0BKOQk8UeLxYeCZgv8/Xeb1P+fm/CrZ3hrltWsTFR//5k0rKeqh3e2qLslQAsciqEYQnLg5xfFdbUbL1ExvSwMz8TRLqSwNofX3/j1jxwfu3dGq+MoMpkKjDNtbI9yZS5DO5io6/sTNKXa0N1TspzZ4QzQUoC0arNh6G5tPcHMybgR2HXD8/MMVppC+c3sagHu2t6q6JIONEQRl2N7WQE5Wlu4mpeStm1M8tMssLvVge1tDxYH9kzetxeX47rKZygZFOLuL3bA79K7Hj65PcqS32XSH1YARBGXY2W5NWmeLvLW4ORlnYiHFQ3uMIKgHO9qi+X2i1+Otm9NEgj4TKK4D++xeQc4+EGuRSGc5cWua9+/rUH1ZBowgKIvT8vbynfUn7Y+uTQImPlAvDnQ3cnNysaLeUK9dm+D+HW2mb00daIkG6WwM5/fsWIt3bs+QyuSMINCEuRvK0NkYpj0W4opd2bgWr1waY3trg+ltXycO9jSRk6y7wIzMLnFxdJ7HDqktSDSUZ19XrCKL4EfXJvD7BA8bK1sLRhCswYHuxnyJezmSmSw/uDrB44e7TBZKnXD60Kw3Vq9cslqTPHaolk4oBi/Y191YkSB45fI49/S30BQx8QEdGEGwBge3NXHlzgJSlm9f8NaNaeKpLI+bxaVu7O6MEfSL/AYm5Xjl0hh9LREOmp2u6sb+rkam42nG5ssnYQzPLHF6cJYnj27TeGVbGyMI1uBgTxPzycyaOerfvXCHUMDH+4wvs24E/T72djZy0W4tXYqlVJYfXJngscPdxnKrI/faxWGnbs+UPebb50YB+PG7enRckgEjCNbkwZ1WiuGJW6U7aGeyOV48PcwTh7uJhlzV5hlcct+OVt6+PVO2+dx3LtxhMZXlJ+7pLfm8QQ939bUQ9AveXkMQ/PW5O+zvbjSbBmnECII1ONTTRFMkwJs3pks+/4OrE0wspPj4/ds1X5mhmEf2tjO7lObiaGn30DffGaKnOcKje4zlVk8iQT939bXw9u3S99TwzBKv35jkmWPGGtCJEQRr4PcJHtrdzps3Jks+//V3hmiOBEwWygbgkb3WAv9GibGaXEjy6uVxnr2vD5/PuIXqzQM72zg9OFMy3fcvTg4iJfy9B3fU4cq2LkYQrMOje9u5Nr64qmBpbC7BS2dG+Pj92wkH1u+bYlDL9tYGdrQ38P0rq/tDffX122Rykp8+3l+HKzMU86GDnSTSuVW9vNLZHH/21gDv29vBzo5ona5ua2IEwTo8dZflU/6rMyMrHv/yD2+SzUl+4YN76nFZhhI8dVcP378yzkx8eceyeCrDH71+k8cPdbG/2+xytRF4374OGsMB/ur06IrHv/7OEEMzS3zmQ3vrdGVbFyMI1mFnR5Tju9r44zdukbEb0N2ejPPlH97gJ+/tY1eHKSLbKHzi/n7SWcmfvjmQf+x3//YaEwspnv/I/jpemaGQcMDPT93Xx4unh/P7TS8kM/zmdy5ztLfZuFrrgBEEFfA/fXgfA1NL/N73rrOYzPAvvnaKoE/wb54+Uu9LMxRwtK+ZDx3s4vOvXGVgKs5r1yb4ve9d49n7+njQNATcUPzDD+whk5P88ovnyWRz/PtvnGVkLsF/+sQxk95bB8RaxVIblePHj8sTJ05oO5+Ukuf/5B3+6swIDUE/iUyW3/70/fzEPX3arsFQGTcmFvmp3/4ByWyOdDbH3s4Yf/k/f8B0sNyA/MZ3LvOfX75CQ9DPUjrLv3jyIP/rEwfqfVmbGiHESSnl8VWPG0FQGalMjv/ywxvcmFjkkw/0mx4oG5irYwv8wWs3aG0I8Y9/bK8RAhsUKSX/7eQgb96Y4vHD3Tx9rMdYA4oxgsBgMBi2OOUEgYkRGAwGwxbHCAKDwWDY4hhBYDAYDFscIwgMBoNhi2MEgcFgMGxxjCAwGAyGLY4RBAaDwbDFMYLAYDAYtjjvyYIyIcQ4cKvGl3cCq3sV1x9zXdVhrqs6zHVVx0a9LnB3bbuklKu6+r0nBYEbhBAnSlXW1RtzXdVhrqs6zHVVx0a9LlBzbcY1ZDAYDFscIwgMBoNhi7MVBcEX630BZTDXVR3muqrDXFd1bNTrAgXXtuViBAaDwWBYyVa0CAwGg8FQgBEEBoPBsMXZtIJACPGUEOKSEOKqEOJzJZ4XQojfsp8/LYR4QMM17RBC/K0Q4oIQ4pwQ4p+WOOYxIcSsEOKU/fMfVF+Xfd6bQogz9jlX7fpTp+/rUMH3cEoIMSeE+GdFx2j5voQQXxZCjAkhzhY81i6E+I4Q4or9u63Ma9eciwqu6/8SQly0x+nrQojWMq9dc8wVXNcvCSGGCsbqmTKv1f19/deCa7ophDhV5rUqv6+Sa4O2OSal3HQ/gB+4BuwFQsC7wNGiY54B/jsggEeBNzRcVy/wgP13E3C5xHU9BrxYh+/sJtC5xvPav68SYzqKVRCj/fsCPgQ8AJwteOzXgM/Zf38O+NVa5qKC6/oYELD//tVS11XJmCu4rl8C/lUF46z1+yp6/teB/1CH76vk2qBrjm1Wi+Bh4KqU8rqUMgX8GfBs0THPAn8oLV4HWoUQvSovSko5IqV82/57HrgAbFd5Tg/R/n0V8QRwTUpZa0W5K6SU3wOmih5+FviK/fdXgI+XeGklc9HT65JSfltKmbH/fR3o9+p8bq6rQrR/Xw7C2jD5Z4A/9ep8lbLG2qBljm1WQbAdGCj4f5DVC24lxyhDCLEbuB94o8TT7xNCvCuE+O9CiLs0XZIEvi2EOCmE+EyJ5+v6fQGfovwNWo/vC2CblHIErBsZ6C5xTL2/t3+IZcmVYr0xV8Hztsvqy2XcHPX8vn4MuCOlvFLmeS3fV9HaoGWObVZBIEo8VpwnW8kxShBCNAJ/AfwzKeVc0dNvY7k/7gV+G/iGjmsCPiClfAB4GvhFIcSHip6v5/cVAn4K+PMST9fr+6qUen5v/xbIAH9c5pD1xtxrPg/sA+4DRrDcMMXU7fsCPs3a1oDy72udtaHsy0o8VtV3tlkFwSCwo+D/fmC4hmM8RwgRxBroP5ZS/mXx81LKOSnlgv33S0BQCNGp+rqklMP27zHg61jmZiF1+b5sngbellLeKX6iXt+XzR3HPWb/HitxTL3m2XPATwA/K21HcjEVjLmnSCnvSCmzUsoc8KUy56vX9xUAPgn813LHqP6+yqwNWubYZhUEbwEHhBB7bG3yU8ALRce8APy8nQ3zKDDrmGCqsH2Qvw9ckFL+32WO6bGPQwjxMNYYTSq+rpgQosn5GyvYeLboMO3fVwFlNbV6fF8FvAA8Z//9HPDNEsdUMhc9RQjxFPCvgZ+SUsbLHFPJmHt9XYUxpU+UOZ/278vmo8BFKeVgqSdVf19rrA165piKCPhG+MHKcrmMFU3/t/ZjnwU+a/8tgN+xnz8DHNdwTR/EMtlOA6fsn2eKrut54BxW5P914P0armuvfb537XNviO/LPm8Ua2FvKXhM+/eFJYhGgDSWBvYLQAfwMnDF/t1uH9sHvLTWXFR8XVexfMbOHPtC8XWVG3PF1/VH9tw5jbVQ9W6E78t+/A+cOVVwrM7vq9zaoGWOmRYTBoPBsMXZrK4hg8FgMFSIEQQGg8GwxTGCwGAwGLY4RhAYDAbDFscIAoPBYNjiGEFgMBgMWxwjCAwGg2GL8/8DMGbUGgu7dfkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(t, y1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the actual times (this assumes 3 suns orbiting at earth-sun distance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "tt = 1/np.sqrt(6.67e-11 * 1.99e30 / (1.5e11)**3 ) # seconds\n",
    "tt = tt / (60*60 * 24* 365.25) * np.diff(t)[0] # per time step (in years)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAHWCAYAAACBsnu3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAkKUlEQVR4nO3dfZRVdd338c+X4cEEUYkBBYHhRppusi4UFoiSHXy6ZXIJXikmpKbpLHsku328XUmZhZZlGZlOJWYSCXURlINoxsHI9EIMESJyFAgauhRQbCIZB773H+dAw8yZmTNzzjycL+/XWmfNPnv/9v79vvNj+MzZe58z5u4CAACFrVtnDwAAAOSOQAcAIAACHQCAAAh0AAACINABAAiAQAcAIICcA93MhpjZcjPbYGbrzWxmhjZmZveaWZWZrTWzU3LtFwAA/Fv3PByjTtL/dfcXzOwoSavN7El3/1O9NpMljUw/xkv6fvorAADIg5xfobv7dnd/Ib38D0kbJA1u0GyKpIc95VlJx5jZ8bn2DQAAUvJ6Dd3MSiSdLOm5BpsGS9pa7/k2NQ59AADQRvk45S5JMrM+kn4h6fPu/lbDzRl2yfiZs2ZWLqlcko444ogxQ4cOzdcQu5T9+/erW7e49yRSX2GjvsIVuTYpfn1/+ctfdrh7cVv2zUugm1kPpcJ8nrv/V4Ym2yQNqff8BEnVmY7l7hWSKiSptLTUN27cmI8hdjnJZFKJRKKzh9FuqK+wUV/hilybFL8+M9vS1n3zcZe7SfqRpA3u/q0mmi2RdHn6bvdTJe129+259g0AAFLy8Qr9dEmXSXrJzNak1/0/SUMlyd3vl1QpqUxSlaQ9kq7MQ78AACAt50B395XKfI28fhuX9Olc+wIAAJnFvbMAAIDDCIEOAEAABDoAAAEQ6EAbRH7bjER9QCEi0AEACIBABwAgAAIdAIAACHQclh566CGZmR566KHOHgoA5AWBjoJnZq16EOKH2rp1qz71qU9p/PjxOu6449SrVy8NGjRIH/zgBzV37ly98847WR/r5Zdf1l133aUzzzxTQ4YMUc+ePTVw4EBNmTJFy5cvz+oYdXV1Gj9+vMxMixYtarLdokWLZGYaP3686urqsh4jEFXe/toa0FlmzZrVaN23v/1t7d69WzNnztQxxxxzyLbRo0dr+PDhOvXUU3X88cd30Ci7rldeeUXz5s3T+PHjNXXqVPXr1087d+7U0qVLddVVV+nhhx/Wk08+qe7dW/7v4otf/KIeffRRjRo1SmVlZerXr582btyoJUuWaMmSJfrOd76jz33uc80eo3v37nrkkUd08skn65prrsk4T3//+99VXl6u3r1765FHHslqbEB47t5lH+95z3s8quXLl3f2ENpVZ9c3bNgwl+SbNm3q1HEUgr179/q+ffsara+trfVEIuGS/NFHH83qWHPnzvUXXnih0fpkMuk9evTwnj17enV1dVbHeuCBB1ySn3vuub5///5Dtk2ePNkleUVFRVbHOpx09s9ee4ten6TnvY2ZySl3HJaauoZeUlKikpIS1dTU6LrrrtOQIUP0rne9S6NHj9Yvf/lLSalTwl/72tc0cuRIHXHEERoxYoTmzJnTZF/Lli1TWVmZ+vfvr169emnEiBG64YYb9Oabb7Zfga3Qs2fPjH9fukePHpo6daqk1Kn0bHz84x/XySef3Gj9hz70ISUSCdXW1uqZZ57J6ljl5eW64IIL9MQTT+i73/3uwfX33Xefli5dqgsuuEDXXHONJGnPnj2aPXu2Ro8erd69e6tPnz6aMGGC5s+f3+i4tbW1mjNnjsrKyjRs2DD16tVL/fr109lnn62lS5dmHMuBfxdvvfWWvvCFL6ikpEQ9evTQl770JUnSP/7xD33lK1/RSSedpL59++qoo47SiBEjdMkll2j16tVZ1QvkivNUQAPvvPOOzjnnHO3atUtTpkxRbW2t5s+fr4985CN64okndN999+m5557T5MmT1atXLy1cuFCf/exnVVxcrEsuueSQY91+++2aNWuW+vXrp/PPP18DBgzQ2rVrdffdd6uyslJ/+MMf1Ldv306qtHn79u1TZWWlJOkDH/hAzsfr0aOHJLXq9PgPf/hDvf/979dNN92ks88+W927d9cNN9yggQMH6oc//KEk6c0339SZZ56pP/7xjzrllFN01VVXaf/+/Vq2bJmmT5+u9evX64477jh4zF27dmnmzJk67bTTdM4556i4uFjbt2/Xr371K5WVlekHP/iBrr766kZjqa2t1Zlnnqldu3bp3HPPVd++fTV8+HC5u8477zw988wzmjBhgq6++mp1795dW7duVTKZ1Ac/+EGNGTMmx+8e0DICHWigurpap5xyipLJpHr16iVJuuyyy3TGGWfo4osv1ogRI7Ru3bqD1+a/8IUv6L3vfa/uvPPOQwJ9+fLlmjVrliZMmKDKyspDruU/9NBDuvLKKzVr1izdc889LY5p8+bNrb6Z7+Mf/7hKSkqybr9jxw7NmTNH7q7XX39dTz75pKqqqjR9+nSdf/75req7oS1btuipp57SkUceqTPOOCPr/YqLi/Xggw/qwx/+sGbMmKEePXpoz549WrhwoYqLiyVJn//85/XHP/5Rd911l2688caD+7799tuaOnWqvva1r+miiy7S6NGjJUnHHnustmzZohNOOOGQvnbv3q3TTz9dN954o2bMmKF3vetdh2zfvn27Ro0apRUrVqh3794H17/00kt65plnNHXq1EY38e3fv1+7d+/Oul4gJ209V98RD66hF67Orq+la+hz5851ST537tyM+1VVVTXaZ/jw4S7Jn3rqqUbbEomEFxUVeV1d3cF1U6dOdUm+bt26jGMYPXq0FxcXZ1XP8uXLXVKrHq2dgw0bNhyyv5n59ddf77W1ta06TkNvv/22n3766S7Jv/71r7fpGJ/85CcPjutTn/rUwfU7duzwoqIiHzt2bMb91qxZ45L8hhtuyKqfb37zmy7JV6xYccj6A/8u1qxZ02iftWvXuiS/9NJLW1FR++nsn732Fr0+5XANnVfoQAPHHHOMRowY0Wj9oEGDtGnTpoynTwcPHqx9+/bp73//uwYPHixJ+sMf/qAePXpo4cKFWrhwYaN9amtr9frrr2vnzp1697vf3eyYEomEUj/r7ee9732v3F379u3T3/72Ny1atEi33XabVq5cqccee0z9+vVr9TH37dunyy67TL///e91ySWX6Prrr2/T2L75zW/q+9//viTp7rvvPrh+1apV2rdvn8zs4PXs+g685W7Dhg2HrF+/fr2+8Y1v6Omnn9b27dv19ttvH7L9b3/7W6NjHXHEERkvPYwaNUqjR4/W/PnztWXLFk2ZMkUTJ07U2LFj1bNnz1bXCrQVgQ40cPTRR2dcf+Dab6btB7bVf8/2zp07VVdXpy9/+cvN9ldTU9NioHekoqIiDR06VDNnztTAgQN16aWX6rbbbmv2xr9M9u3bp4997GNauHChpk2bpkceeURm1qYx1T/9XX95586dklLBvmrVqib3r6mpObj87LPP6swzz1RdXZ3OOussXXDBBerbt6+6deumNWvWaPHixdq7d2+jYwwYMCDj+IuKivTb3/5Wt99+u37+85/rpptukiQdddRRuuKKKzR79mz16dOn9UUDrUSgA+3k6KOP1v79+7Vr166cj9UR19AzmTx5siQpmUy2ar+6ujpNnz5dCxcu1PTp0/Xwww+rqKgop7FkcuCXq+uuu07f+ta3strnjjvu0L/+9S8tX7680V9dmz17thYvXpxxv+Z+GTn22GN1zz336J577lFVVZVWrFihBx54QHPmzNGbb76pn/zkJ9kVBOSAQAfayamnnqrHHntM69ev1/ve976cjrV58+YWX+k3lEgkcg70A6eeW3Nnem1traZNm6bFixfr8ssv19y5czO+LS4fxo0bp27duul3v/td1vtUVVWpX79+Gf+E6ooVK3Ie04knnqgTTzxR06dP14ABA5r8BQHIN96HDrST6667TpJ0zTXXqLq6utH2f/7zn3r22WezOtaBa+iteWT7N7+fe+457dmzp9H6mpoazZw5U5L04Q9/+JBtu3fv1p///Gdt3779kPV79+7VhRdeqMWLF+sTn/hEu4a5lDoNPmPGDD3//PP6yle+kvEjYF955RVt2rTp4POSkhLt2rVLa9euPaTdj370Iy1btqzVY9i0aZPWr1/faP0bb7yhvXv3NrpbHmgvvEIH2slZZ52lO++8U7fccotGjhypsrIyDR8+XDU1NdqyZYtWrFihiRMn6vHHH+/Ucc6ePVvJZFIf+tCHNHToUB155JHaunWrli5dqjfffFOnnXaabrnllkP2WbRoka688kpdccUVh1wKuPbaa1VZWan+/ftr8ODBuv322xv1l0gksv5lIxtz5szRyy+/rNtuu00/+clPNHHiRA0cOFDV1dXasGGDVq1apfnz52v48OGSUm9zW7ZsmSZOnKhp06bp6KOP1vPPP6+VK1fqoosu0s9//vNW9f/iiy/qwgsv1JgxY3TSSSdp0KBBev3117V48WK98847B6+pA+2NQAfa0U033aTTTz9d9957r1auXKnFixfr6KOP1uDBg1VeXq7p06d39hB1zTXXqHfv3lq1apWSyaT27NmjY489VmPGjNG0adN01VVXZX3K/cAr4R07dmQM8wPyGeh9+/bVihUrVFFRoZ/+9Kf6xS9+obffflsDBw7UyJEjdc899+icc8452P68887Tr371K91xxx169NFHVVRUpHHjxmn58uV69dVXWx3oY8eO1S233KIVK1bo8ccf1xtvvKHi4mKNGTNGn/vc5w7ehwC0N2vvt8LkorS01Ddu3NjZw2gXyWQyr/+pdTXR6wO6qug/e9HrM7PV7j62LftyDR0AgAAIdAAAAiDQAQAIgEAHACAAAh0AgAAIdAAAAiDQAQAIgEAHACAAAh1og9b+9bFCQ31A4SHQAQAIgEAHACAAAh0AgAAIdAAAAiDQAQAIgEAHACAAAh0AgAAIdAAAAiDQAQAIgEAHACAAAh0AgAAIdAAAAiDQAQAIgEAHACCAvAS6mT1oZq+Z2bomtifMbLeZrUk/bstHvwAAIKV7no7zkKQ5kh5ups3v3P38PPUHAADqycsrdHd/WtKufBwLAAC0nrl7fg5kViLp1+5+UoZtCUm/kLRNUrWk6919fRPHKZdULknFxcVjFixYkJfxdTU1NTXq06dPZw+j3VBfYaO+whW5Nil+fZMmTVrt7mPbsm9HBXpfSfvdvcbMyiR9x91HtnTM0tJS37hxY17G19Ukk0klEonOHka7ob7CRn2FK3JtUvz6zKzNgd4hd7m7+1vuXpNerpTUw8z6d0TfAAAcDjok0M3sODOz9PK4dL87O6JvAAAOB3m5y93M5ktKSOpvZtskzZLUQ5Lc/X5JF0n6pJnVSfqXpI96vs71AwCA/AS6u1/awvY5Sr2tDQAAtAM+KQ4AgAAIdAAAAiDQAQAIgEAHACAAAh0AgAAIdAAAAiDQAQAIgEAHACAAAh0AgAAIdAAAAiDQAQAIgEAHACAAAh0AgAAIdAAAAiDQAQAIgEAHACAAAh0AgAAIdAAAAiDQAQAIgEAHACAAAh0AgAAIdAAAAiDQAQAIgEAHACAAAh0AgAAIdAAAAiDQAQAIgEAHACAAAh0AgAAIdAAAAiDQAQAIgEAHACAAAh0AgAAIdAAAAiDQAQAIgEAHACAAAh0AgAAIdAAAAiDQAQAIgEAHACAAAh0AgAAIdAAAAiDQAQAIgEAHACAAAh0AgADyEuhm9qCZvWZm65rYbmZ2r5lVmdlaMzslH/0CAICUfL1Cf0jSec1snyxpZPpRLun7eeoXAAAoT4Hu7k9L2tVMkymSHvaUZyUdY2bH56NvAADQcdfQB0vaWu/5tvQ6AACQB+bu+TmQWYmkX7v7SRm2PSZptruvTD9/StKN7r46Q9typU7Lq7i4eMyCBQvyMr6upqamRn369OnsYbQb6its1Fe4Itcmxa9v0qRJq919bFv27Z7vwTRhm6Qh9Z6fIKk6U0N3r5BUIUmlpaWeSCTafXCdIZlMKmptEvUVOuorXJFrk+LXl4uOOuW+RNLl6bvdT5W02923d1DfAACEl5dX6GY2X1JCUn8z2yZplqQekuTu90uqlFQmqUrSHklX5qNfAACQkpdAd/dLW9jukj6dj74AAEBjfFIcAAABEOgAAARAoAMAEACBDgBAAAQ6AAABEOgAAARAoAMAEACBDgBAAAQ6AAABEOgAAARAoAMAEACBDgBAAAQ6AAABEOgAAARAoAMAEACBDgBAAAQ6AAABEOgAAARAoAMAEACBDgBAAAQ6AAABEOgAAARAoAMAEACBDgBAAAQ6AAABEOgAAARAoAMAEACBDgBAAAQ6AAABEOgAAARAoAMAEACBDgBAAAQ6AAABEOgAAARAoAMAEACBDgBAAAQ6AAABEOgAAARAoAMAEACBDgBAAAQ6AAABEOgAAARAoAMAEACBDgBAAAQ6AAAB5CXQzew8M9toZlVmdnOG7Qkz221ma9KP2/LRLwAASOme6wHMrEjS9ySdI2mbpFVmtsTd/9Sg6e/c/fxc+wMAAI3l4xX6OElV7v6qu9dK+pmkKXk4LgAAyFI+An2wpK31nm9Lr2togpm9aGZLzex9eegXAACkmbvndgCziyX9H3e/Ov38Mknj3P2z9dr0lbTf3WvMrEzSd9x9ZBPHK5dULknFxcVjFixYkNP4uqqamhr16dOns4fRbqivsFFf4YpcmxS/vkmTJq1297Ft2Tfna+hKvSIfUu/5CZKq6zdw97fqLVea2X1m1t/ddzQ8mLtXSKqQpNLSUk8kEnkYYteTTCYVtTaJ+god9RWuyLVJ8evLRT5Oua+SNNLMhptZT0kflbSkfgMzO87MLL08Lt3vzjz0DQAAlIdX6O5eZ2afkbRMUpGkB919vZldm95+v6SLJH3SzOok/UvSRz3Xc/0AAOCgfJxyl7tXSqpssO7+estzJM3JR18AAKAxPikOAIAACHQAAAIg0AEACIBABwAgAAIdAIAACHQAAAIg0AEACIBABwAgAAIdAIAACHQAAAIg0AEACIBABwAgAAIdAIAACHQAAAIg0AEACIBABwAgAAIdAIAACHQAAAIg0AEACIBABwAgAAIdAIAACHQAAAIg0AEACIBABwAgAAIdAIAACHQAAAIg0AEACIBAB4BCNm+eVFIidesmlZRowG9+09kjQich0AGgUM2bJ5WXS1u2SO7Sli0qvfvu1Hocdgh0AChUt94q7dlzyKqivXtT63HYIdABoFD99a+tW4/QCHQAKFRDh7ZuPUIj0AGgUH31q9KRRx6yal+vXqn1OOwQ6ABQqGbMkCoqpGHDJDNp2DBtvP761Hocdgh0AChkM2ZImzdL+/dLmzfrtbPP7uwRoZMQ6AAABECgAwAQAIEOAEAABDoAAAEQ6AAABECgAwAQAIEOAEAABDoAAAEQ6AAABECgAwAQAIEOAEAAeQl0MzvPzDaaWZWZ3Zxhu5nZventa83slHz0CwAAUnIOdDMrkvQ9SZMljZJ0qZmNatBssqSR6Ue5pO/n2i8AAPi3fLxCHyepyt1fdfdaST+TNKVBmymSHvaUZyUdY2bH56FvAAAgqXsejjFY0tZ6z7dJGp9Fm8GStjc8mJmVK/UqXsXFxUomk3kYYtdTU1MTtjaJ+god9RWuyLVJ8evLRT4C3TKs8za0Sa10r5BUIUmlpaWeSCRyGlxXlUwmFbU2ifoKHfUVrsi1SfHry0U+TrlvkzSk3vMTJFW3oQ0AAGijfAT6KkkjzWy4mfWU9FFJSxq0WSLp8vTd7qdK2u3ujU63AwCAtsn5lLu715nZZyQtk1Qk6UF3X29m16a33y+pUlKZpCpJeyRdmWu/AADg3/JxDV3uXqlUaNdfd3+9ZZf06Xz0BQAAGuOT4gAACIBAR+eYN08qKZG6dUt9nTevs0d0eOD7DoSVl1PuQKvMmyeVl0t79qSeb9mSei5JM2Z03rii4/sOhMYrdHS8W2/9d6gcsGdPaj3aD993IDQCHR3vr39t3XrkB993IDQCHR1v6NDWrUd+8H0HQiPQ0fG++lXpyCMPXXfkkan1aD9834HQCHR0vBkzpIoKadgwySz1taKCG7PaG993IDTuckfnmDGDIOkMfN+BsHiFDgBAAAQ6AAABEOgAAARAoAMAEACBDgBAAAQ6AAABEOgAAARAoAMAEACBDgBAAAQ6AAABEOgAAARAoAMAEACBDgBAAAQ6AAABEOgAAARAoAMAEACBDgBAAAQ6AAABEOgAAARAoAMAEACBDgBAAAQ6AAABEOgAAARAoAMAEACBDgBAAAQ6AAABEOgAAARAoAMAEACBDgBAAAQ6AAABEOgAAARAoAMAEACBDgBAAAQ6AAABEOgAAARAoAMAEED3XHY2s36SHpVUImmzpGnu/kaGdpsl/UPSPkl17j42l34BAMChcn2FfrOkp9x9pKSn0s+bMsndRxPmAADkX66BPkXSj9PLP5Y0NcfjAQCANsg10Ae6+3ZJSn8d0EQ7l/SEma02s/Ic+wQAAA2YuzffwOw3ko7LsOlWST9292PqtX3D3Y/NcIxB7l5tZgMkPSnps+7+dBP9lUsql6Ti4uIxCxYsyLaWglJTU6M+ffp09jDaDfUVNuorXJFrk+LXN2nSpNVtvTTdYqA3u7PZRkkJd99uZsdLSrp7aQv7fElSjbvf3dLxS0tLfePGjW0eX1eWTCaVSCQ6exjthvoKG/UVrsi1SfHrM7M2B3qup9yXSLoivXyFpMUNG5hZbzM76sCypHMlrcuxXwAAUE+ugX6npHPM7GVJ56Sfy8wGmVllus1ASSvN7EVJ/y3pMXd/PMd+AQBAPTm9D93dd0o6K8P6akll6eVXJf1HLv0AAIDm8UlxAAAEQKADABAAgQ4AQAAEOgAAARDoAAAEQKADABAAgQ4AQAAEOgAAARDoAAAEQKADABAAgQ4AQAAEOgAAARDoAAAEQKADABAAgQ4AQAAEOgAAARDoAAAEQKADABAAgQ4AQAAEOgAAARDoAAAEQKADABAAgQ4AQAAEOgAAARDoAAAEQKADABAAgQ4AQAAEOgAAARDoAAAEQKADABAAgQ4AQAAEOgAAARDoAAAEQKADABAAgQ4AQAAEOgAAARDoAAAEQKADABAAgQ4AQAAEOgAAARDoAAAEQKADABAAgQ4AQAAEOgAAARDoAAAEkFOgm9nFZrbezPab2dhm2p1nZhvNrMrMbs6lTwAA0Fiur9DXSfpPSU831cDMiiR9T9JkSaMkXWpmo3LsFwAA1NM9l53dfYMkmVlzzcZJqnL3V9NtfyZpiqQ/5dI3AAD4t464hj5Y0tZ6z7el1wEAgDxp8RW6mf1G0nEZNt3q7ouz6CPTy3dvpr9ySeWSVFxcrGQymUUXhaempiZsbRL1FTrqK1yRa5Pi15eLFgPd3c/OsY9tkobUe36CpOpm+quQVCFJpaWlnkgkcuy+a0omk4pam0R9hY76Clfk2qT49eWiI065r5I00syGm1lPSR+VtKQD+gUA4LCR69vWLjSzbZImSHrMzJal1w8ys0pJcvc6SZ+RtEzSBkkL3H19bsMGAAD15XqX+yJJizKsr5ZUVu95paTKXPoCAABN45PiAAAIgEAHACAAAh0AgAAIdAAAAiDQAQAIgEAHACAAAh0AgAAIdAAAAiDQAQAIgEAHACAAAh0AgAAIdAAAAiDQAQAIgEAHACAAAh0AgAAIdAAAAiDQAQAIgEAHACAAAh0AgAAIdAAAAiDQAQAIgEAHACAAAh0AgAAIdAAAAiDQAQAIgEAHACAAAh0AgAAIdAAAAiDQAQAIgEAHACAAAh0AgAAIdAAAAiDQAQAIgEAHACAAAh0AgAAIdAAAAiDQAQAIgEAHACAAAh0AgAAIdAAAAiDQAQAIgEAHACAAAh0AgAAIdAAAAiDQAQAIgEAHACCAnALdzC42s/Vmtt/MxjbTbrOZvWRma8zs+Vz6BAAAjXXPcf91kv5T0gNZtJ3k7jty7A8AAGSQU6C7+wZJMrP8jAYAALSJuXvuBzFLSrre3TOeTjezTZLekOSSHnD3imaOVS6pXJKKi4vHLFiwIOfxdUU1NTXq06dPZw+j3VBfYaO+whW5Nil+fZMmTVrt7k1ewm5Oi6/Qzew3ko7LsOlWd1+cZT+nu3u1mQ2Q9KSZ/dndn87UMB32FZJUWlrqiUQiyy4KSzKZVNTaJOordNRXuCLXJsWvLxctBrq7n51rJ+5enf76mpktkjROUsZABwAArdfub1szs95mdtSBZUnnKnUzHQAAyJNc37Z2oZltkzRB0mNmtiy9fpCZVaabDZS00sxelPTfkh5z98dz6RcAABwq17vcF0lalGF9taSy9PKrkv4jl34AAEDz+KQ4AAACINABAAiAQAcAIAACHQCAAAh0AAACINABAAiAQAcAIAACHQCAAAh0AAACINABAAiAQAcAIAACHQCAAAh0AAACINABAAiAQAcAIAACHQCAAAh0AAACINABAAiAQAcAIAACHQCAAAh0AAACINABAAiAQAcAIAACHQCAAAh0AAACINABAAiAQAcAIAACHQCAAAh0AAACINABAAiAQAcAIAACHQCAAAh0AAACINABAAiAQAcAIAACHQCAAAh0AAACINABAAiAQAcAIAACHQCAAAh0AAACINABAAiAQAcAIAACHQCAAAh0AAACyCnQzewbZvZnM1trZovM7Jgm2p1nZhvNrMrMbs6lTwAA0Fiur9CflHSSu39A0l8k3dKwgZkVSfqepMmSRkm61MxG5dgvAACoJ6dAd/cn3L0u/fRZSSdkaDZOUpW7v+rutZJ+JmlKLv0CAIBD5fMa+lWSlmZYP1jS1nrPt6XXAQCAPOneUgMz+42k4zJsutXdF6fb3CqpTtK8TIfIsM6b6a9cUnn66V4zW9fSGAtUf0k7OnsQ7Yj6Chv1Fa7ItUnx6ytt644tBrq7n93cdjO7QtL5ks5y90xBvU3SkHrPT5BU3Ux/FZIq0sd+3t3HtjTGQhS5Non6Ch31Fa7ItUmHR31t3TfXu9zPk3STpAvcfU8TzVZJGmlmw82sp6SPSlqSS78AAOBQuV5DnyPpKElPmtkaM7tfksxskJlVSlL6prnPSFomaYOkBe6+Psd+AQBAPS2ecm+Ou5/YxPpqSWX1nldKqmxDFxVtHFohiFybRH2FjvoKV+TaJOprkmW+7A0AAAoJH/0KAEAAXSbQo3+MrJldbGbrzWy/mTV5h6aZbTazl9L3JLT5bseO1or6CnX++pnZk2b2cvrrsU20K5j5a2kuLOXe9Pa1ZnZKZ4yzrbKoL2Fmu9NztcbMbuuMcbaFmT1oZq819bbeAHPXUn0FO3eSZGZDzGy5mW1I/785M0Ob1s+hu3eJh6RzJXVPL98l6a4MbYokvSLpf0nqKelFSaM6e+xZ1ve/lXp/YVLS2GbabZbUv7PH2x71Ffj8fV3SzenlmzP9+yyk+ctmLpS6D2apUp8lcaqk5zp73HmuLyHp15091jbWd4akUySta2J7wc5dlvUV7Nylx3+8pFPSy0cp9dHpOf/8dZlX6B78Y2TdfYO7b+zscbSXLOsr2PlTapw/Ti//WNLUzhtKXmQzF1MkPewpz0o6xsyO7+iBtlEh/1trkbs/LWlXM00Kee6yqa+guft2d38hvfwPpd4B1vATVFs9h10m0Bs4nD9G1iU9YWar05+aF0khz99Ad98upX4YJQ1ool2hzF82c1HI85Xt2CeY2YtmttTM3tcxQ+sQhTx32Qoxd2ZWIulkSc812NTqOczpbWut1dEfI9vRsqkvC6e7e7WZDVDq/f1/Tv+22unyUF/Bzl8rDtNl56+BbOaiS89XC7IZ+wuShrl7jZmVSfqlpJHtPbAOUshzl40Qc2dmfST9QtLn3f2thpsz7NLsHHZooHsHf4xsR2upviyPUZ3++pqZLVLq1GGXCIQ81Few82dm/2Nmx7v79vRpr9eaOEaXnb8GspmLLj1fLWhx7PX/A3X3SjO7z8z6u3uEzwkv5LlrUYS5M7MeSoX5PHf/rwxNWj2HXeaUu/ExsjKz3mZ21IFlpW4UjPTHaQp5/pZIuiK9fIWkRmckCmz+spmLJZIuT99te6qk3QcuOxSAFuszs+PMzNLL45T6/3Bnh4+0fRTy3LWo0OcuPfYfSdrg7t9qolnr57Cz7/ard0dflVLXC9akH/en1w+SVNngzr+/KHUH662dPe5W1HehUr9x7ZX0P5KWNaxPqTtyX0w/1kerr8Dn792SnpL0cvprv0Kfv0xzIelaSdeml03S99LbX1Iz787oio8s6vtMep5eVOpG3NM6e8ytqG2+pO2S3kn/3H0i2Ny1VF/Bzl16/BOVOn2+tl7mleU6h3xSHAAAAXSZU+4AAKDtCHQAAAIg0AEACIBABwAgAAIdAIAACHQAAAIg0AEACIBABwAggP8P2SDuwKeLkhYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def animate(i):\n",
    "    ln1.set_data([x1[i], x2[i], x3[i]], [y1[i], y2[i], y3[i]])\n",
    "    text.set_text('Time = {:.1f} Years'.format(i*tt))\n",
    "fig, ax = plt.subplots(1,1, figsize=(8,8))\n",
    "ax.grid()\n",
    "ln1, = plt.plot([], [], 'ro', lw=3, markersize=6)\n",
    "text = plt.text(0, 1.75, 'asdasd', fontsize=20, backgroundcolor='white', ha='center')\n",
    "ax.set_ylim(-2, 2)\n",
    "ax.set_xlim(-2,2)\n",
    "ani = animation.FuncAnimation(fig, animate, frames=1000, interval=50)\n",
    "ani.save('plan.gif',writer='pillow',fps=30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}