{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sympy as smp\n",
    "from scipy.integrate import odeint\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import animation\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib.animation import PillowWriter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Position confined to some 1D path: $x=x(\\theta)$ and $y=y(\\theta)$ A few examples\n",
    "\n",
    "* Parabola: $x = \\theta$ and $y=\\theta^2$\n",
    "* Simple Pendulum: $x = \\cos(\\theta)$ and $y=\\sin(\\theta)$\n",
    "* Tautochrone: $x = \\sin(2\\theta) + 2\\theta$ and $y=1-\\cos(2 \\theta)$\n",
    "* Many more ...\n",
    "    \n",
    "\n",
    "Kinetic Energy: $ T = \\frac{1}{2}m (\\dot{x}^2 + \\dot{y}^2) $\n",
    "\n",
    "Potential Energy: $ V = mgy $\n",
    "\n",
    "The Lagrangian: $L = T- V $\n",
    "\n",
    "Lagranges Equation: $\\frac{dL}{d\\theta} - \\frac{d}{dt} \\frac{dL}{d\\dot{\\theta}} = 0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "t, m, g = smp.symbols('t m g')\n",
    "the = smp.symbols(r'\\theta', cls=smp.Function)\n",
    "the = the(t)\n",
    "the_d = smp.diff(the, t)\n",
    "the_dd = smp.diff(the_d, t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d^{2}}{d t^{2}} \\theta{\\left(t \\right)}$"
      ],
      "text/plain": [
       "Derivative(\\theta(t), (t, 2))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "the_dd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make sure $x$ and $y$ are functions of $\\theta$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = smp.symbols('x y', cls=smp.Function)\n",
    "x = x(the)\n",
    "y = y(the)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define $x$ and $y$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "path='taut'\n",
    "if path=='taut':\n",
    "    x = smp.sin(2*the) + 2*the\n",
    "    y = 1 - smp.cos(2*the)\n",
    "    x_f = smp.lambdify(the, x)\n",
    "    y_f = smp.lambdify(the, y)\n",
    "if path=='parab':\n",
    "    x = the\n",
    "    y = the**2\n",
    "    x_f = smp.lambdify(the, x)\n",
    "    y_f = smp.lambdify(the, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define $T$, $V$, and $L$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 1/2 * m * (smp.diff(x,t)**2 + smp.diff(y,t)**2)\n",
    "V = m*g*y\n",
    "L = T-V"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - g m \\left(1 - \\cos{\\left(2 \\theta{\\left(t \\right)} \\right)}\\right) + 0.5 m \\left(\\left(2 \\cos{\\left(2 \\theta{\\left(t \\right)} \\right)} \\frac{d}{d t} \\theta{\\left(t \\right)} + 2 \\frac{d}{d t} \\theta{\\left(t \\right)}\\right)^{2} + 4 \\sin^{2}{\\left(2 \\theta{\\left(t \\right)} \\right)} \\left(\\frac{d}{d t} \\theta{\\left(t \\right)}\\right)^{2}\\right)$"
      ],
      "text/plain": [
       "-g*m*(1 - cos(2*\\theta(t))) + 0.5*m*((2*cos(2*\\theta(t))*Derivative(\\theta(t), t) + 2*Derivative(\\theta(t), t))**2 + 4*sin(2*\\theta(t))**2*Derivative(\\theta(t), t)**2)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute $\\frac{dL}{d\\theta} - \\frac{d}{dt} \\frac{dL}{d\\dot{\\theta}}$ (which is equal to zero:Lagranges equation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "LE = smp.diff(L, the) - smp.diff(smp.diff(L, the_d), t)\n",
    "LE = LE.simplify()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle m \\left(- 2.0 g \\sin{\\left(2 \\theta{\\left(t \\right)} \\right)} + 8.0 \\sin{\\left(2 \\theta{\\left(t \\right)} \\right)} \\left(\\frac{d}{d t} \\theta{\\left(t \\right)}\\right)^{2} - 8.0 \\cos{\\left(2 \\theta{\\left(t \\right)} \\right)} \\frac{d^{2}}{d t^{2}} \\theta{\\left(t \\right)} - 8.0 \\frac{d^{2}}{d t^{2}} \\theta{\\left(t \\right)}\\right)$"
      ],
      "text/plain": [
       "m*(-2.0*g*sin(2*\\theta(t)) + 8.0*sin(2*\\theta(t))*Derivative(\\theta(t), t)**2 - 8.0*cos(2*\\theta(t))*Derivative(\\theta(t), (t, 2)) - 8.0*Derivative(\\theta(t), (t, 2)))"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve for $d^2 \\theta / dt^2 $ so we can solve it with an ODE solver. Our system of equations will then be\n",
    "\n",
    "$$ \\frac{d \\theta}{dt} = \\omega$$\n",
    "$$\\frac{d \\omega}{dt} = \\frac{d^2 \\theta}{dt^2} =  \\text{whatever is returned} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "deriv_2 =smp.solve(LE, the_dd)[0]\n",
    "deriv_1 = the_d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0.25 \\left(- g + 4.0 \\left(\\frac{d}{d t} \\theta{\\left(t \\right)}\\right)^{2}\\right) \\tan{\\left(\\theta{\\left(t \\right)} \\right)}$"
      ],
      "text/plain": [
       "0.25*(-g + 4.0*Derivative(\\theta(t), t)**2)*tan(\\theta(t))"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "deriv_2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convert this into a numpy expression:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "deriv2_f = smp.lambdify((g, the, the_d), deriv_2)\n",
    "deriv1_f = smp.lambdify(the_d, the_d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$S = (\\theta, \\omega) $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create ODE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dSdt(S, t):\n",
    "    return [\n",
    "        deriv1_f(S[1]), #dtheta/dt\n",
    "        deriv2_f(g, S[0], S[1]) #domega/dt\n",
    "    ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.linspace(0, 20, 1000)\n",
    "g = 9.81\n",
    "ans1 = odeint(dSdt, y0=[np.pi/4, 0], t=t)\n",
    "ans2 = odeint(dSdt, y0=[np.pi/5, 0], t=t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1fdef50f520>]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABnlklEQVR4nO39d3hkyX3fC3+qczdCNxpo5IzJm5ezicslKe1S3CUpLiVLMuVEW5b47pXoK1v2tXgtWb7v49fXtH1tS7Ik8+WVKNGSrkhJlESaWjEtc1juzqbZMAEYZKABdAMd0DnV/aPOwWAxwCB0n+5zMP19HjzoUH3q1Kmqb/2q6vv7lZBS0kQTTTTRxPGHrdE30EQTTTTRRH3QJPwmmmiiiVsETcJvookmmrhF0CT8JppooolbBE3Cb6KJJpq4ReBo9A3cDF1dXXJ0dLTRt9FEE000YRm88MILUSllaLfvTE34o6OjXLhwodG30UQTTTRhGQgh5vb6rrmk00QTTTRxi6BJ+E000UQTtwiahN9EE000cYugJoQvhHhcCHFFCDElhPjoLt/7hRD/UwjxihDidSHEP6pFvk000UQTTRwcVRO+EMIO/DbwBHAO+GkhxLkdyX4BeENKeRfwTuA/CyFc1ebdRBNNNNHEwVELC/9+YEpKOS2lLACfBp7ckUYCbUIIAbQCG0CpBnk30UQTTTRxQNSC8AeAhW3vF7XPtuO3gLPAMvAq8ItSyspuFxNCfFgIcUEIcSESidTg9ppoookmmoDaEL7Y5bOdMZffDbwM9AN3A78lhGjf7WJSyk9IKc9LKc+HQrv6DuyL33xmkm9evXUGi2yhzB98d4YvvrbS6FupG/KlMn96YYEvXFzmVgnxXa5IPv3cPJ97eemWKvOfXVjgM8/PU67cGmXOFcsUSrvaw1WjFo5Xi8DQtveDKEt+O/4R8DGpWumUEGIGOAM8V4P834RCqcLvfnuaZK7Ef/yJO/mp80P7/8jCqFQkP/c/LvCdqSgAv/a+c/zM28YafFfGQkrJL/7Jy3zxdTXATb8rzf/66MkG35Xx+OhnL/JnLywC8NpSgl95786tsuOH/99fv8Hvf3cWgJfm43zsb93Z2BuqA/7rV6/ytUtrfO4jD+Nz1dY3thYW/vPASSHEmLYR+0Hg8zvSzAOPAggheoDTwHQN8r4BLoeN53/1MR4cD/Jvv/AGiWzRiGxMg/95cZnvTEX5Nz96jkfPdPOfvnSFtc1co2/LUHzjSoQvvr7CL73rFD96Vz//7WuTzK9nGn1bhuK5mQ3+7IVFnnrHBD99/zC/+50ZLq8kG31bhuL15QS//91Z/t6Dw3z47eN8+vkFXpjbaPRtGYrZaJrf+/YM9w531JzsoQaEL6UsAR8BvgRcAv5USvm6EOIpIcRTWrJ/C7xVCPEq8Azwy1LKaLV57wW3w86vvOccm7kSn9UsouOK//H9OSZCLXzooVF+9X3nyBbLfPq5hf1/aGF88rszDAS8PPWOCX71vWeREv7w2dlG35ah+KNn5/B7nfzTx07yy4+fxu2w8fvfmW30bRmKP/z+HD6Xnf/t3Wf4Z4+dIuBz8knN2j+u+JPn55HAL/3IKUOuXxMdvpTyaSnlKSnlhJTy32mffVxK+XHt9bKU8keklHdIKW+XUv5RLfK9Ge4Y9HPnoJ+/fGnJ6KwahtlomhfmYvzk+SFsNsFYVwtvnejkz19YPLZrvJHNPN+divJj9wzgctjoaffw2Nke/vKlJSrHdI03kS3ypddXePLufjxOOwGfi/fd2c/Tr4Uplo1Z6200csUyf/1qmMdv68XvdeJ12fmxewb48usrJHPHc9ZeLFf47AtL/NDpbnraPYbkcaw9bd99Wy+vLiVYTR7PJY4vXFxGCPjA3ddFUe+9s4/5jQzXIqkG3plxePrVMBUJT97dv/XZ47f3Ek0VuLiUaOCdGYdnLq2SL1X48XsHtz579229bOZK/GD6eC5xfOtqhM1ciSfvud62n7i9j2JZ8p1JwxYHGooX52JEU3n+1r07RY61w7Em/MfO9gDwzKW1Bt+JMfj2ZJTb+tvp9V+3Bt55uhtQ69zHEd+6GmGsq4WTPW1bn73jVAibgK9dWm3gnRmH70xGCba4uHPAv/XZIye78DrtfPmN46nM+vZkFJ/LzkPjnVuf3TscoN3j4GuXj29/ttsED5/sMiyPY034p3paGQh4+fbk8SO/TKHEi/MxHp54c+MYCHg52d16LAm/VK7w3MwGD010vunzjhYX9w53HEsprpSS70xFeetEJzbbdQW0x2nnwfEg37u23sC7Mw7fuxbl/rEgLsd1inLYbbz9VIhvXo0cyyXLb01GuGcoQLvHaVgex5rwhRDcN9rBhbnYsWsgz8/GKJYlD5+40Rp4+EQXL8zFKB2z9d3XlpNs5ktvsvp03D8W5LXlJJnC8XLgnlpLsbaZ55FdrL77xoJMraXYSBcacGfGYSWR41okfYMxA/DQRCeRzTzzG8dLlZXMFXl1KbFrf64ljjXhA5wfDR7LBvLiXAybgLeMdNzw3T3DAbLFMpdXNhtwZ8bhuRllzT64C+GfH+2gXJG8snC81vFfmo8Dqh3vxH3aZxdmj9c6/gVNevnA+I1l1tv7C3Oxut6T0Xh1MYGUu/fnWuIWIHz1AC/MHq8GcnExzonuVlrcN2p17x1WZX5p/riVOcFAwEuozX3Dd3qZXzxmZX5lMU6b28FYZ8sN390x4Mdlt/H8MSP8VxcTuOw2zvTe6Ix/sruNNrfj2NXzywtxAO4c9N88YZU49oR/qruNVrdj64EeB0gpeXUpwZ2DgV2/H+xQpPiiZh0eF7y2lOCOgd07RMDn4kR3Ky8eM8vv4mKCOwb9b1q/1+Fx2jnb386rx0yd9OpSgrN9bW9av9dhtwnuHg7wwly8/jdmIF5ZiDPa6SPgMzaI8LEnfJtNcLavjUvh4+OVuJzIEU0V9rQGhBDcNRg4VkSQyBaZXc9wx00soNv6249VPedLZS6vJPcc2EGV+Y3l5LHZo6pUlDFz+x4DO8DtA34mVzfJl8p1vDNjcXExwV1DAcPzOfaED3CuTxHBcXHMeXUxDrCntQtwrq+N6UiKXPF4dIrXtcHrZkRwrq+d5USOeOZ4bGJeCm9SLEvuuskgd66vnWSuxGIsW8c7Mw5zGxk2c6WbLm2c62unVJFMrR0PX5PVZI6VZI67bjKw1wq3BuH3t5MulI/Nxu3FxQQOm+Bs364BRwE409dORcLk6vHoFPps5WaDnP483jgmVv7WwL7PrAaOT5kvbhkzgT3T6PV8KXw8RAmv6W3b4PV7uFUIv089yOMy3b8UTnKiuxWP075nmq1OcUwCbL22nGQg4CXYsvca53Ejgiurm7S5HQwEvHumOdPbjk3AG8vHo57fCCdx2W2c7GndM81YVwsep+3Y9OermlF2apszoVG4JQj/ZE8rdps4NlbQ1dXUmzxNd8NI0IfXaT82nWJydZPTvTcvc6jNTajNfWzI7+pqilO9baiD4naH12VnPNTK68ekzFOrKca6WnDa96Ymu01wuuf47MtNrm7S0+7G7zXO4UrHLUH4HqedE6HWY0EE6XyJpXiWU917W0CgNqtP97Zx+RhYu6VyhelI+qZWn45zfe3HYmCXUnJ1dZNTByjz2b72YxMqeXItxYkDlvlS+HhsVk+upepi3cMtQvigrPypYxBQbFLbqNrPwgetU6xYv1PMbWQolCuc7N6/zGd627gWSVn+dKRIKk88UzxQmU+EWlmKZ8kWrL1BnyuWWYhlOLmPMQOqnmOZIpHNfB3uzDhUtM3ng9RzLXDLEP5EqJWFjYzlVStXV5XFfhDL72R3K/FM0fKu9/rG80GIYDzUQqFUYcniqhW9zPstYwGc6G5FSiwfIfVaJIWUHIj8JrS2cC2SNvq2DMViLEu2WD5Qf64Fjifh/9FPwJd/FfLXO8BEdysVCbPr1m4gk6ubuBw2RrZ7Xi6/DH/+M/A3vwzZ645H4yGVZjpq7TJPralB7sR2wl+7DJ/9OfjCL0HqetC0iZBGBFFrk58+sL9pGWv9GvzVz8Pn/wkkr58iemKL/Kxd5qmt2eu2MicW4XMfUeXeuH5Inl7P08exnjdmYPKrhuRXkzO0hBCPA78B2IHflVJ+bJc07wR+HXACUSnlO2qR9w0oZEAI+N5vKSL8B58Dm1rDB9WodnPZtgom11JMhNQmNKCI7w/eC8IGhTQsPAc/8yVwuBjv0jpFJLUVd8WKuLqaYiDgvR5GIjYHn/wRqJShlIe578LPfR1cPsZ1wl9L8UNaqGgr4upqioDPSahVCyOxuQqffBwKKVXu6W/CU98Gj5/RLh82ocpsZUyuprDbBKO6MZONwSefgPQaIGDqGXjqO9Aaorfdg9dp59qatY2Zq1vGjDaryafg//nbkFmHX3wF3LW1/Ku28IUQduC3gSeAc8BPCyHO7UgTAH4HeL+U8jbgJ6vNd0+4fPB3/wze/99g9tvw0h8CytoVAss3kKm11Jst3b/+52B3wS/8AH7y92H5RfjBxwEY6PDictiYtvi0d3It9WYL6Ev/CioVRXh/5zMQuQzf+S8ABFtcdPiclp/VTEfUwL6l0PnKr0EuAT/7DPzDLyjL9+v/J6CO9BwO+iy/RzW5tslIp+96SIWv/TvYXIZ/+Nfws1+B7IZ6DrB1wpvVLfyp1RS97Z7rCp3v/xZEr8BP/F7NyR5qs6RzPzAlpZyWUhaATwNP7kjzd4C/kFLOA0gpjT/B4J6/BwPn4Vv/GSplPE47gx1eS3eKfKnMcjzLWKdPfbDwPMx9B97xL6G9H849CRM/DN/9DSgVNGvJZ+l1TiklM9HU1myFtUtw+Qvw0C9AcBxOPArnPgA/+P8rQgTGQ62Wt3Zn19PXLd34Arz6Z3D/z0HPORi6H+7+O3Dhk1vLWSe6Wy3veTq5lrq+T5Nagxf+AO75+zB4HnrvgPs/DBc/A7FZQBlxVjdmpiLbjJlsXK1MnP1RGH+nIfnVgvAHgO2nZi9qn23HKaBDCPENIcQLQoh/UIN8bw4h4K0fgcQ8zHwTUOt+ViaChY0MFQmjXRoR/ODj4G5Xg5uOB38eMlG48jQA412tlraC1jbz5IoVxrq0Qe65T4DDCw/8f64neus/gXwS3vgcABOhFktb+JlCidVkfmsPhuf/b/X/gaeuJ3roF6BcgFf/FFB7VDPRtGXPQCiWK8ytZ67PXl/5E6gUVXvW8eDPg6zAS38MaEKMmHWFGMqY2Tawv/4XUNiEt/2SYXnWgvB38wrZqYlzAG8B3gu8G/jXQohdj2UXQnxYCHFBCHEhEqnyBKNTT4DHDy//CaAayHQ0ZVmZ4kxUhYYY7WpR6/VXnoY7fhLc21QNEz8M7QPwkjonfjzUwvx6xrKHXc9qxD3S2QLloiL1M+8B37Y9iYG3QOcJeOXTgLLwI5t5yx52PavXc2cLSAmv/QWceAwCQ9cTdZ+F/nu3yG+8q4ViWbIct+b5zYuxLOWKvF7ml/4Ihh6E0Daa8A/A+DvUICcl46EWpIS5dWuGTIlnimzmSozoM/aX/wRCZ6H/HsPyrAXhLwLbWiKDwPIuab4opUxLKaPAt4C7druYlPITUsrzUsrzoVCoujtzeuC2H4PLfw2lPKOdPnLFCmsW1e7q5Dfe1QJTX4ViRi3jbIfNDrf/OEx/A3JJxkOtlCrSsnGE9M482tmiZmqZdbj9b705kRBw5wfV5m0yfF3BYdHp/oxWz6NdPlh6ARILcNsHbkx41wdh7XWITm2ptuY2rFnmuXW9zC0QnYToVbhzl62+O35KLeksvbitnq05g9UVg6OdLZBYgsXn4M6fUu3ZINSC8J8HTgohxoQQLuCDwOd3pPkc8IgQwiGE8AEPAJdqkPf+OPUEFNMw9z2GtU4xa9Hp/sx6moDPqWJmX/oC+Dph5OEbE556Qk2Hp79+XZppVfJbT+OwCfoDHrj6JXD6YOLRGxOefkL9n/rqVpmtunz3JiK4/AWwOeD0e25MeOrd6v/kl7aWBWYtau3qA/tI0AdXv6g+PPX4jQlPP6EUaZNfZkxb2rSqHHXLmOnywbVn1Id6nRqEqglfSlkCPgJ8CUXifyqlfF0I8ZQQ4iktzSXgi8BF4DmUdPO1avM+EMYeAbsbJr/CqDZ1mrOstZu+PuWd/oZavrHvoqwdekAtZV35IhNd1raC5tbTDAV9OOw2VeaRt6qZ20703AZt/TD1FYY6lExxzqI+FzPRND3tbiVDnf4mDN4H3sCNCTtGoes0XP0S3W1uPE4bcxY1ZubWM3iddnWa2dUvQc/t4B+8MaEvqJbwpr5Ki9tBb7vHsvs1s+tphIDBDp+SnLb1Qfe5/X9YBWrieCWlfFpKeUpKOSGl/HfaZx+XUn58W5r/JKU8J6W8XUr567XI90BwtcDow3DtGQYCXhw2YVkimI1mlFUTuaK0yWNv3z2h3aF2+We+hd/nJOBzWnZJZzaaUQN1YklN8/dSLwihFDvXvoHLJunzey1cZm1gz8Yh/DKM3cRl5eS7YP772Mp5hoM+SxszI50+RCkHCz9QxsxeOPGYWurKbDDc6WNxw5pe1XPrGfr9Xjw23YB71NDlHDiunrY7MfIwRC7jyMcZ6PBacpMnVyyznMgqIpj9tvpw9JG9fzDyMCQXIT7PcNBnSfKTUmpE0LKltLqpXG3s7ZBPwOrrjHRas8ygLL+xrhaY+55Spew1sIOa8ZQLsPwiI50tljVm5jYyDAd9sHhBLUfutlSpY+ztgIT5Zy3btkHV80inTw3quTicuMkgVyPcGoQ//JD6v/ADRjpbLNlA5jcySKmt9818E/zDakq/F/Qyz32foaCPBQuWOZoqkC6UlYU/8y21Z9F9294/GHpA/V/4gUYE1rP8krki0VRBbV7OfFNJUAfP7/2DoQfV/7nvMaoNclY72a2iiQpGu1pg/vuAgOEH9v5B/z1gc8KCIvyVZM6S0sy59YwyZuafVR+MvM3wPG8Nwh+4VzWQ+e8zEvRZctNWV26MdfqU5Tf6tptP/3puUxr9ue8yHPSxFM9aLoKkbq2OdLWokBFDD4LtJk02MKzW8efVIBdN5ckUSnW629pgfkuV5FPkN3QfONx7/6ClU63jzz/LSGeLJVVoK8kchVJFWfhz31Nt19ux9w+cXui/G+bVwA5Y7ojHZE4FNRzp9Cl1TmAY2noMz/fWIHynV1kF888y0ukjmStZ7txT3UIftUeVNPFmVh8oeebQAzD/fYaDPoplyUrSWhptXXEy3lKAjWsw+Jab/0BoluE2IliwmJWv1/Nwuw1WX1cblPth5CE1e+1QA4PVAgRuqVU63LD4/PXZ6c0w9AAsv8SwX536ZrUZ7NbAHvQqY2bw/rrke2sQPsDwg7D0ImMBpWqxmnxtMZal1e2gbUMTNw3cu/+Phh+A6FXGWpWVO2+xMs+tp7HbBP3pN9QHA/sMcqBmAclFJlyxrWtYCQsxjfCLM1ApKeeq/TD8EOSTnBRLgPXqeV7zHTjBggoON3ST5Rwdww9COc9YYVK7hrXKrA/KE+4EbIZVuIw64NYh/IG3QKXISeYB6xHBYizDYIcXsfyiCpZ2s7VsHZrH3nhxCrCeFTS7nmEg4MUZfhEQB/NAHLwPgOHcFcB6RLAYy9LucdAavag+OMjArj2XUOoSTrtgxmJte3Y9g8MmCKUuqw/6797/R5pFHNh4Ba/Tbrl61mc1wxnNgGsSfo3Rpxx7+zIaEVjMClrYyCq97vJLSqPs2Psw7y30KSLoTL6B3SYs1ynmNzIMBb1KuRE6A54DhLXuuQ2EnZaN12nzOCw3yC1sZBgK+lTU05ZuFSZjP3SeAGcL9pVX6A94LbeePb+uymxfvQiuVghO7P+jth5o7UWsvGpJpc7ceppQmxv32ivKT6jn9rrke+sQfscoePw41y4SanNbqlNIKVmMZRjqcKsY/weNtdHSCf4h7CuvMBCwni59KZZlKOBV5HeQtWxQTlndZxHhVyxJBIuxLIMdXlh6UVn3B9Fl2+wqmmT4FQY7vCzGrFXm2fW02nNZfhl677z5xvx29N0F4VcsqUJbjGUZ6vDCyqsqLpLd+APM4VYifCFUA1l+mcEOL0tx6xB+LFMkXShzmzuioukdZJqvo+8uCL/MUNBahJ8rlomm8pzypdUmde8dB/+xVubhDmuVWQ3sWSbapXIyO8j6vY7+u2HlVYb81jJmQJHfcIdLkV/friG2dkffXRC9wnjApsmWraNCW4xlGQhohH+Ytl0lbh3CB+i7G9beYNjvsJQVpN/r6bLaoDocEdwDG9Oc8pctZQXpA/Jp25z6oPcQU96+uyAd4VxbmoVY1jK69GiqQLZY5g77HCAPObDfDcUMd7hXiWzmLaNLT+VLJLJFbnNFoJQ92Pq9jr67QFa407lIplBm3SJnN5crknAiy5nWlGbM3Fm3vG8xwr8LygXudq+wFLcOEejSwv7cNbXe17VrZOndoVlMdzkXWU8XSOWtoUvXrdTh4oz6oOcAm9Q6tDLfbpulUKqwumkNOao+sI9XZtUHh1nX1YjyNOrc12WLzGD1w+ZPV5Sw4NAWPnCyfA2wzgb92maOYllyDt2YaVr4xkAbSU+JeYplaRkHFZ0I/KkpFR98t4Bpe0ELxjShDhuzjJWvl7kzdRX8Qzd3xNkJjSjHior8rLJBv6CRX09uBjwBaOs9+I87T4DdRX9hFrCOI5Jez/25q8qruPPkwX/sHwRvB73Zq4B16lkf5EZK2qHshzFmqsStRfjBcbC7GCrOAlhmWWchliHgc+KIXj58NL32fvD46c/PqmtZhPCXYlkcNoF349LhO4S7FQIjdGZVh7LKfo1eN23JSVXPhwmkZXdC1ymCKWUpW4Xw9brp2JzUNi8PYcwIAd230ZZQS51W6c963YTSV6Fj7GDqsxrh1iJ8uwO6TtGZsRYRLMaynPaXIbmkOsVhIAR0n6N9U3UKK5V5xG9DRCePJlnrPkuLRgRWWd5YjGXp9DmxRy6rs2sPi+6zuGNXcNiEZchvKZbF5bDh3Lh6+LYN0H0GW/QKnT4nSxY57Uvvgy2xS3VdzoFbjfABus/ijaspoFWsoIWNDOd92rnvR4mX3X0WR/QSHqfYmk6aHUvxLPe3RECWD7dhq6P7LLb1SXpbhIUGuQx3+9Mq4ueRyO8sIrHICb+0TNtejGc54y8hUqsQOn34C3SfhXySu/xpS9XzgE9i25ium/5exy1J+LbkIiMtJUtYQbpU7zbHovrgSERwDpFPcmd7huWEdTrF3S4VKuBoFv45qJR4S1vMOuQXy3Kvb0W9OdLArn7zQMuqJdo2qDJvGTOhI7Rt7Td3e1ZYslCZ72+PAhK6z9Q171uQ8LVO0bpmCSKIpPLkSxXG5bzyQvQP7f+jndDKfN67YgkLP18qs5rMMyEWVRiJjrHDXySkOtI97rAlLL9KRbIUy3LWVs3Arn5zh2vZEm0b1JLO7U7tCOyjkJ9W5rP2RZbjOUto8ZdiWe50r6o3XUeY1VSBmhC+EOJxIcQVIcSUEOKjN0l3nxCiLIT4iVrkeyRoDeQud9gS5Le4XbnRffZoJ+JoZb7dsWiJdc6wdo/9pQVNfXKIjTwdXadA2DhtW2Q5njU9Eaxu5iiUK4yWZ1WI58OoknT4h8HVyoRcYM0CWnzduW6cxaMbM74gtPYwUp4jWywTyxRrf6M1hJSSpXiWU7YlEHYlJKkjqiZ8IYQd+G3gCeAc8NNCiBvmo1q6/4A6+7Zx8A+Ds4VTYoFFC2jxlXJD0p484qYWqE7R1sdYZY5oyvxEoA9yHZlZ6DqETG87nB4ITjBUmiNXrLBhcqecLeVGdvro9WyzQegMA5rvgtk3q/X76y/MagP0EY/3C51RBhGY3ojTZ+yD5QXonDhYTKwaohYW/v3AlJRyWkpZAD4NPLlLun8CfBZYq0GeR4fNBt1nGCwqp5xoytxa/MVYli6S2HOx6g44Dp2mp6C0+OGEua38pXgGF0U8qfnDOZntRPcZQtlr2jXNTQRLsSw2KrQkrx2d8AG6z9ChSTPNXmZ9kAukqy3zOdo2ryGomL7M+oDUlZ2prm0fEbUg/AFgYdv7Re2zLQghBoAfAz7OPhBCfFgIcUEIcSESidTg9nZB1yk6ssrLbcHkFkE4keUeb1i9CVWxwdN1ivb0LCBNbwUtxrKM21cRslLdGmfoDL70Ai6Kprd2w4kcQ2INWzlfHfl1ncKZW6edtOnX8ZfiWfykcGUj1bXt7jPYSlkGRcT0hL8Yy+KkhC81fzRVUpWoBeHvNg/buU7y68AvSyn3XUuQUn5CSnleSnk+FArV4PZ2QecJPNlVWsiaXs2wkshxlzeq3lRjEXSexFFMESLOUtzcZV6MZXmLTxvsj7qkA9B5EiErDItV05NfOJHldrc2+T2Mt+lOaL89aV8xfdteimU5bdeUWNUQvlbmM4410xszS/EsI2IFIcvVlfmIqAXhLwLbd1sGgeUdac4DnxZCzAI/AfyOEOIDNcj7aNCIc0yYX8GxHM9xyrECTp/ymj0quk4AcMK2bPqN26VYlju2VAxVkJ9W5rPOVUvU810ebZDrPHH0C2nP6x5fxPTHOy7Fs5z3avVcjTxRK/PdvojpZ3KLsQx3ezTprUWXdJ4HTgohxoQQLuCDwOe3J5BSjkkpR6WUo8CfAz8vpfyrGuR9NGgN5A6P+aWZ4USWEbmsNniOuqkFW43rbm/U9FbQYizDSXtYqTZcLUe/kGb53eWNmL7MK8ksJx0rSp3T0nn0C3WMgs3Bba5VS5DfOfeqZswMHv1CLSFw+zljgYF9KZblLs8qIKxJ+FLKEvARlPrmEvCnUsrXhRBPCSGeqvb6hiA4DsLGHW5zWwQ5TWbWW1yobpoPSurn9HGbe8XUZS6WK6wkcwyWFqqz7kHFKGnt4bRjxfQOZ+F4jmG5XH09253QMca4LWz+zflYljGxovrjQQ892Q1CQNcJRlg2ddsGLUyKfRkCQ+Dy1T3/Iwicb4SU8mng6R2f7bpBK6X8h7XIsyo43BAY5kQhzKdMvLwRTuRwUaQ9F65umg+qQ3WeYHzT3MtYK4kcUlbozM1B1w9Vf8HOkwxHlk1t4eeKKpZ7j30But5d/QW7TtI/f5mVZI5yRWK3VTEzNAj6wN7nX4LOA57gdjN0nqQ38jXW0wWyhTJel736a9YYugZ/qGUReuu/YQu3oqetjs6TDFWWTG35hRNZhsUqgkr11i4oIigtEk6Y1/9gOZ6ljw0c5WyNynyCnsICsUyRTMGcZwGsJnO0kKW1EFVLd9Wi8wTB3AKyUiZi0hDgK4kcNlkikF+u3pgB6DpBa2ENHznTGjRxrQ12FpZqU89HwK1L+F0n6covkMqZ91CQcDzHuNAkmTXpFKfw58PYynkiJvU/CCdyTNi0Pf9ayNY6T+IpJeggadrp/nI8x9hWPddikDuFXRYZFBHTGjSLMSWjtMlSbdq29tzGhHmXLMOJHJ0kcZUzaq+lAbilCd9ZydHHBmHTNpBsbQm/8wQCyahYMe1mdTiRY0JohF+LTS1tljAuwqYt80pyWz3XaCYHMC7Mu6YdTmjr91AjY+Z6mc1q4YcTWUaEpko6SnyoGuDWJXzNIhi3hVk26eZWOJHjrGsVWntqc0iC1ikmTE4EZ5wr4PYr9UW16Lpez2YlguV4jnFbGImoDRF06vUc3opLZDaEEzWevQbHkQgmbOaNkRVO5BgWmq9FsEn49cU2i8C8Fn6OE7bV2nQI2LrOuIn9D8KJHBOOKARHq5Oh6giMIO0uTpiaCLKcdqwiAsMqBlC1aOlEeoOcdoRNu6QTTmQ57VxVMlRfsPoLOr0I/xDnXGumNmZGbWtqYA+MNOQebl3Cb+1ButqYMLGFvxzPMiyXakf4rhZoH+CU07wa7XAiyxCrtZvy2uyI4Lipna9WEjlO2ldqs5yjQXSe4JRj1bQW/koixylHDY0ZgK4TnLCFWTRpPYcTOU67ooj2/toM7EfArUv4QiA6xznpMK8WP5OI0FZJ1LZTBMeZsJvXEWktnqarvFrbTa3OE4zZVkyrS1+OZRmsLNVmw1ZHcJxBVgmb1MJfjucYlmEI1lCt0nmS/soSyyYNKRGO5xi1rzVswxZuZcIHCI4zajNnp8gWynTmtJh0NbT86BhlUK6YclaTL5VxpldwyFJt1ziDY/SWw6zE07W7Zg1RSizhkbnaSvWC43SVI6zHE7W7Zg0RTyQIliM1N2Y8lSylzVVTyo5Xkjn6KysN27CFJuHTU15lNWY+IlAqhhpuaukIjuOvxEjEN2p3zRphNZFn2KZtatXSCuoYwymLyOSK6YggVyzTmdcG9hrXM4AnvUi+ZK7zD3LFMu1Zvcy1HeQA+iurRNPmkh1LKYnFY/jLG2p/qkG4tQm/Yww7ZUgumO5EpHAix4htFSlstd3g0TpFILdItmAuIlhOZK+rGGppBWllHpBh1k12EIphyg2tzKNildWEuchvJZEzyJhRz29EmG/vIp4p0l1urCQTbnXC1zpFXzlM3GRHoy3HlWa33Npf21NxtDKPiFVWkubqFCuJHCNiFWlzgL+KYFo7sb3MJlvKCmv1XLE5qwsgthNb5Ge+OELL22evtTziLzCMFDZtmdZc9by8XYPfIEkmNAkfUERgtk6hyG8N0VnjMy+1xjYqVk0nR9UtfOkfBlsNY6H4B6nYnKasZ2Xhr1JuGzra2b17wRek7A4wItZMt0elLPwVSi294G6t3YUdbiptAwyLVVZMVuZwXNUz0LTwG4a2Xip2jymngMuJHKO2Ney1tgbcbZR8IUaE+VQr4XiOMfsatlqX2Wan4h9mRKyYz8LXLD9bpwEkEBxlVKywbLK2HU7kGLatIQywdG2d44wJ81n44aQy4Cru9qMdUF8j3NqELwSVwKiydk1mEcRjUTpIGjL9EyZVJ+nWrhGyNXvnBGM2E1r4ceWMY6/1TA6tzHbzWfi6A1LNjRlABMcYta2Zj/DjWUZtq2qQq4VD4RFxaxM+YO+aYNS2YrpToGyxWfXCgOmfvXPClJ1iMxahTaYNGuTGGBFrrJhsGSsZW6ONjDHrusFx+omwGtus/bWrQDSWpIcNY/TowXECJEnEorW/dhVYSeQYtUcQDVzOgSbhIzrHGRZrptNou1Pz6oVBRNDLOtGYuTTazuSsemFEpwiO00KW9MZq7a9dBWyxGfXCoDLbqVDemK/9tatAJTanXhjRtrXnaE/M1v7aVSAcT9EnIw3dsIUaEb4Q4nEhxBUhxJQQ4qO7fP93hRAXtb/vCSHuqkW+NUFwHDdFchtLjb6TLaTzJUJFTcVgCBFo14zP1v7aR0SuWMaf06JkGmT5ATiTM7W/dhXwpY0d2AHcm3O1v3YV8OjGjIH13JqeN5XPhYwv4qTUUC9bqAHhCyHswG8DTwDngJ8WQpzbkWwGeIeU8k7g3wKfqDbfmkEjVKeJLIKwJk/MuzpqEyVzJzRy8ZqICFaTuW2hY0drn4FGBC3pBdMQQbZQpqtg4CCnte1Qcck0h7/kimU68kaWWV1zUJrH+UpKeX3GfgyWdO4HpqSU01LKAvBp4MntCaSU35NSxrS3zwI1FBxXCY0I2jLmsQj0k66K7QZF1NPK3FVYMo3z1bImWyt4Omsr1dMRGKaCjSFWTON8tZzIMmJbI+vpBqe39hm0dlOyexkVq6ZR6qxojmYlu7c24a93wt1K3hMylfIulimqkApwLJZ0BoCFbe8Xtc/2wj8G/mavL4UQHxZCXBBCXIhEIjW4vX3gH6QsHAyyStQkp0CFtzT4BjUObwd5p59RsWIa56uVpNLgVwKjxmTgcJHz9WlyVHNs3K5oqiTDBnYhyLePKvIzSZmXNWMm3zZsmFql5B9hxETOV8txrW0LJ7TfjBqNRy0If7da29VUFkL8EIrwf3mvi0kpPyGlPC+lPB8KGWAB7ITNTq51SPNINEcDWd3YpF9EcYeMO/eyoBOBSVQry3GlzXYYIE/UUQ6MaY5I5qhn3ZvaEA2+BhEc17T45qhn3cI3ci3b3jVhqkFOH9gL7UO1dSg8AmpB+IvA0Lb3g8DyzkRCiDuB3wWelFKu1yDfmqESGDOV52l+fQa7kIaSn+gc1zqFOcgvEtukX6zj6DKuzI6ucWXhm6SeIxtxekQcd3cN48nsgLvnBIMiQtgkAQLDmrXr7jbOmHGHJugTG0RjccPyOAyUc90aosEbtlAbwn8eOCmEGBNCuIAPAp/fnkAIMQz8BfD3pZRXa5BnTeEMKYtgySxxtDdm1X8D1/s83RMMiCirsaRheRwGxY1ZbEhDN7Xc3ScJihQbG3VYKjwACtFpAJxdxpGfvXMctyiRjZpjg35zfRGvKBhszKjnWYjOGpbHYaDHS3IZOGM/KKomfCllCfgI8CXgEvCnUsrXhRBPCSGe0pL9GtAJ/I4Q4mUhxIVq860l3N0TtIocifVwo28F2CajM5D8HF0T2IUkZ5JOYdclogZaQfrSSTl6zbA8DgOha/CN3MjT25DuyNdgyPVZ9cJIa1crsyNuDglucmOVNpE1JJTEYVGTaE1SyqeBp3d89vFtr38W+Nla5GUEhKZaqaxPA29r7M0A/twiReHG2dZrXCY6EWyYo1P40tq+fx3IzxE3hyOSd7MOUj3tebqS5iizM1m/MntT5igzRjrXHRK3vKctcF2Ln2z8tDeVL9FXXmHTN2hszA2tU7g3G98pcsUyXcVlijY3tPYYl5FmVXrTjS8zQHtugay9rTaHeO+ZyQBl4aA9t2iKMx9aMwtUEBAY2j/xUeELkrO3EcwvmkJq7do00LnukGgSPmgabUFLemH/tAZDbWppsjUj0dpDwea5fvJQA6ErNzItQ8YOcu5WUs4gwcJSw4kgnS/RXw6z6TXYJcVmZ9M7QH8lTDLXWOerXLFMqBwm7e4Gh9vQvFItQwyZQGotpaQts6je1PIgoyOiSfgATg9pV4hgYZlSudLQW9FVDKKWB0PsBiHY9A7QW15puPOVHge/7B81PK9MyzBDJvDC1E+6KvqNJ4F827Ap4uLrZc62GmzMoLT4wyZQoW2kCwzKFTKuELh8Db0XaBL+FjKtwwyLVdY2G0sEsbV5vKJgqAZfR75VdYpGO1+taIOc3UA9uo6Sf4RhW+NPvgrHNhkQ0bqs64rgmCK/BstRdXmiNMq5bhvsneMMiijheMrwvG4GPfZ/zugZ+wHRJHwNZf+oKayg/JpSkLT2nzQ+s+Aow2KNcIPlqInoEj6Rx9tjnB5dhz04Rh8bhNcbK0dNrszgFOW6lNkVOkG7yLIRWTE8r5thLRqjW8RxhgyevQLe7gmcoszm6qzhed0MelwsM6zfQ5Pwt+AMjdMt4qytbzT2RjTVjJHabB2u0AReUSC2tmh4XjdDMaIGOZeBTlc6fL0nsAlJerWx0sx8RGnw2/pOGZ5XW58aVPJrU4bndTNktTK39ho/yLVoeRSjjVWhRTZi9IqYKTT40CT8LbT0Kos6vdJYInBuzlPGBn4DVQwa2rVZRK7BRGCLGxgffQf0ei5pTk8Ng6aLd3YZX2b9NC3ZYAmu/sxddTBmdM27aLD/QWZNldnX0yR8U8GruXqX1htLBO2ZRWKObnC4DM/LFVJWUKOJwJua16R6xq9z2jTyszX4LADX5jxFHNDeb3xmmhzV1eBw2I6E7lA4anxm7QOUcOBNNVhqrfUtm9EijAOiSfgadIvAEW9sA+kqLpP01iminn+IMjZcDfY/8OcWSTpDhkv1AGgJkROehhNBW3aRDWdPfYJpOb3E7Z3X5YENgi+zSEb4jPU70GGzs+Hqoy3b2IONtnx7mmv4JoMvSFq0Xj+BqAHYzBXppw4afB0OFxuObtoaqMXPFcv0VlZJ+YxfwgJACDZc/ddP12oQugp10OBvQ8I7RFdxuaHOV8H8InHPQN0O8U55B+kphRvqc9GaWSBn84Gvs2H3sB1Nwt+GDfcAHfnGWQQr0XVCIgmB+lkDSc8AnYXGkZ+uYijVQY+uI+Uboqe83DAi2MwVGWCVfB306DrybcMMskoy2xjnq2yhTF9llWxLnQZ2oNA+zJBonPOVlJLOQpi4x2Cv+UOgSfjbkGkZorccptgg56v40iQArlD9CD/XNky/XG2Y89VqdJ2QSNQ1sFTRP8IQa0RTjZHgrq2tEBBpZEf9BjnZMUaviBFukAptJZFhSEQo13FgF8Ex/CLD2lpj5Kjr6QJDrJBrrd8gtx+ahL8NZf8IgyLKarwxscNzmgZfl9HVA5XAKF0iyWo0Wrc8t2MzrBRC9dCj67AFx/GIIpHlxizfxZdVmT0GxoTfCV37nlierFue27G+PItbFLHXQXqrw6OdM5BcbkxE9pVYmkERQZogDr6OJuFvg71rHKcos7HcGNVKWQsd2zFgvDZbh64Pji9eqVue21HQNPj+ejiaadAHl1S4MUTQiIFdl+BmVxsjwdX9HlrqKE/U21Qh0hjl3cbKLG5RqotPzUHRJPxt8GnOGqmVxlhBjuQcKXw4W+u3wdPa11giEJo8sh6hJHT4+1U964NNvVHRBvbgQP0GuY6BMwDI9cYYMyXtDIJ6GjN6PW+FJ64z9D6l9zEzoEn42xAYOA1cb5z1RmtmkTVHb103eDqHVJkrG42xgjyb86REC3g76pZnoG+CkrRdd/iqM5zJOeK04WypX5ntLUFlTDRIgmuLz1HGhrtrtG55Cncr66LjenjiOkOfsbcfN8IXQjwuhLgihJgSQnx0l++FEOI3te8vCiHurUW+tUZbaJiitDeMCDoKyyTc9T3V3tPaQYw2nInGlLktu8S6s7+ug5xwuFizhRp2FoAvu8Sao6++mQrBmqOP1kxjJLie9AIRWwjszrrmG3X00d4g/wNncpYSdmxGxv4/JKomfCGEHfht4AngHPDTQohzO5I9AZzU/j4M/Pdq8zUENjsr9t6GnJQjK2V6y6tkG7Cjv2ZvHBF0l5brqkfXEXX24881hgiChWUSnvoO7ABxzyDBBklwO3JLxFx18CregaRvkK5iY8rckl4gau8Ge00OFqwJamHh3w9MSSmnpZQF4NPAkzvSPAn8D6nwLBAQQtTZxDkYNlz9BHL11+JvRpdwiyKVOsrWdMQ8A3Q0gAiyuQJ9co18e/1Dx6Z8Q4QaQQSVMj3lNfKt9R/ksq3D9FZWkeX6a/G7y2FSLfUvc75thG65TrlQ/3DYwcISMU/9y3wz1ILwB4Dt5uGi9tlh0wAghPiwEOKCEOJCJBKpwe0dDsopJwx19kiMLSmVTD1Cx+5EtnWYUCUC5WJd840sTeMSZWwNcDvPtw8TYJNKJl7XfDcjczhFmUodDnvZCdkxikuUSa7Wd/kuuxknSJJSewNOfOoYwSbklhS2XqhUJH3llbo6mh0EtSD83RZfd7LlQdKoD6X8hJTyvJTyfCgUqvrmDoti+witZChs1leXnlnRZWv1k+rpqARGcVAhF5mta75JTRbZiNCxNu3gkXiddemxRVXmekTJ3AldHhhbqq8cdV2T/DZiYHdrAQLjddbixzbW8Is05Toc9nIY1ILwF4Htw9ggsHOufJA0poBd64h6x6wXSuszVKQgOFB/8nNsEUF9tfg5LXSsv45SPR0uzemp3k45GU2P3lpHDb6OFi3P7Gp9B7lUA5zrdFw/C6C+yrvYgupLjjo6mh0EtSD854GTQogxIYQL+CDw+R1pPg/8A02t8yCQkFKGa5B3zeHVvPM2w/XtFPbEHGE66Q601TVfuO5/kF2psxZ/Y4aitNPVX3/Lz9+vBpl6E0F5fYaStNHZX38i6OobpyDtlOusxS9oh5AE+utP+N19w2Skm0qdy5zWNfg95pFkAlS9fSylLAkhPgJ8CbADn5RSvi6EeEr7/uPA08B7gCkgA/yjavM1CgHNGaZYZy2+N73Iiq2XAXv9XSNCvSPkpLPuZwG4N+cIixDD7jqERd6BnlAXUdled6cce2KOJboY8LfUNV+AkN/HPCGcidm65itiMySkj+6e+us0OlpcXKW77lr8snbYS2DwmBE+gJTyaRSpb//s49teS+AXapGX0ejtCrIiOxB1JoJAbolp91vqmqeO3oCXedmNo85E0JpdJOLspxHHOwdbXLxMDx11JgJfeoEVWy8jDRjY7TbBir2PoXR9Jbju1DxLopdzzjrE/t8BIQQRRx/jdZYd2xOzRKSfzo46xP4/BJqetjvQ6nawJHrw1JMIChk6KhtkGrSj73HaCdv68NWZCLqKyyQbJFsTQhB19NFWZyII5JeJ1dm5bjvi7kGChaW6qtDas0usN0CDryPhGaSzUF/lnS+9wIq9D5vNHGGRdTQJfxesO/vx11GLLzXP3lID9Og61FkAy/XrFNkYbTJNrl6HveyCpHeQQCkCpUJ9Mswlaa8kyDRAj64j0zqET2YgG6tPhpUynaVVUvU6xW0X5FqHcZOH1Grd8uzILRFzN26Q2wtNwt8Fm74hOspRKNYnXrp+cLqts3HHoGVah3DLHKTW6pKfHjGSOsaE34l82wh2KpCoj5UvtQO16xkTficqmkxQ1it2UnIJJyXy7aP1yW8XVLTwxOV67VGV8gQrUTL1OsXtEGgS/i4o6g4isfo4qKQ0dUwjNPg6tk6cqtPeRUKTgDZCg69D1pkIMqsqH3sDB3a7doh7uk6KrLwmvRUNjAmvOzOm6lTmysYsNiTFBjjX7Ycm4e8G7YT5eil1itFp0tJNsLtxU0B7pyLeeoUMzmpE0NbASIJuTYufqVM47JQm1fM1cGDX9f/pOsmOk1o+nu7G6dHbe8YpS1G3etYl3WbT4EOT8HeFV2ucqTp1ChGfZV520x/w1SW/3dDSM0ZFiroRgdyYISLb6elq3OHOgdAQWemq21kAxegMSekj1N1Tl/x2Q3ewgxXZUTcJbn7tGkVpp72ncbOanmA7YTrrpsXX+1AjHM32Q5Pwd0Ew1Mem9JKr00k5ns0FFukh1FZ/PbqO3g4/y3RSitanzI7EHPOyh552T13y2w26HFVfWzca+sDe529cmfv9HuZkT90kuDI2y5Lsoq+jtS757YZ+v5f5SjeOOp0FUFxXM/auHnMFToMm4e+KvoCPedldHy2+lLTnllh39mFvoISrz+9hvtKNrU5E0JpZYNXei6cB2mwd/QEv87KnbmcBeFMLLNBNd1vjCL+r1c2i7MaXqs9GtWtzjnnZTW8DB7mAz8mi6KWlTrJjW0wN7L0Bb13yOwyahL8L+jQryF0PiyC1hkvmSfkaaw30amX21oMISgX8xTXiDQ4d2+Fzsih6aM0uGi9HrVRozy2z7uxv6MBuswk23AO0FSN1UaG1aae4NXJgF0KQ8AzQWopBPmV4ft7UPIv00tniMjyvw6JJ+Lugxe1g1d5LW24ZKmVjM9OWEwoN1KODcr6KOPvxFTcgv2lsZvF5bEhyDQ4dq4hgEFclZ7xGezOMQxZJ+RqnR9ex5eBn9FJWNo6vnCRhgpjwdStzpYI/rxzNRB1PcTsomoS/Bza9gzhkETaNjfEmtWUjW3DU0HwOgpSvTp1CK3OlgVI9HTn9hLENg5fvtGdabERM+B3YkgsaXc+aQ2G+rfF69K0wxUYv026Gccri9b5kMjQJfw/k27SOaTAR6CGCvQ04+GQntrT4BpdZl36aQbYmA5p6xGDy2xrYTTDIOTQtvuHOV3o7CjROoaPDoYU9r9RpYG+k1/zN0CT8PSD1wxoMtghykWnCMkh30G9oPgeB6KxPmTOrk2Sli0Co8VaQKzRKWQoqBpNfbm2ashT4uhtPfoHOHpLSe93b2SAUtbDIzlDjyxzs7CEuW8itGltmvR0JExgzu6FJ+HvA1zVMUdqNlynGZlmQIfpMsKMfDIaIyVbDy1yJzjAne0yhYujuUBpto+Pi5yPThOmkp6P+5x3sRG/Ax7zsMbyes2tTbMhWurrqf3LdTvQHPMzLbsP9D7Kr1yhJG77QqKH5HBVNwt8DvR1tLMku8mvGOuW4NudZaLA2W4dSJ3Ub7m1rT8wyJ3tMUeZ+v4e5Sg8Vox2RYrPMV7rpN8Eg1x9Q9WwzeCZX2dDkie2NL3Nvu5Lg2uPGKu+KkWssy056O9oNzeeoaBL+Hujze4x3yinm8OXWWKSbUGvjnK509PlVpxBGdopKhZb0AnMNdrrS0avVs9FOOe7NeeVNbYJBrs/vZUH24EkvGqpCcyTmWJDd9AcaX2Z9kPNmlqBcMiwfEZthVvaawpjZDVURvhAiKIT4ihBiUvvfsUuaISHE14UQl4QQrwshfrGaPOsFXYvvMpIIEgsIJEnPoCniZm/5H6SXoFw0JpPNMA5ZIOrsb6g2W0e/Nsi58wbKUQsZvIV1lkQPXSYY2DtbXCzSg12WIGlQGPByCV9m2TQDu9/rJGzrxWZkmQGP5mh2LAkf+CjwjJTyJPCM9n4nSsA/l1KeBR4EfkEIca7KfA1Hn9+rCL+YNC52+JYGv/Gbl3Dd2rXJsnEhg7VNrUyrOVQMAZ+TZVuvemPUbE677qZnwBQDu80mSLUYLEdNLmKjzLpJBnYhxPVwxUYtZWU2cJeSLIo+giZ0uoLqCf9J4FPa608BH9iZQEoZllK+qL3eBC4Bjfc+2Qdel/36KT1GdQqN/GTQHDv6HqedmFtzkjGqzFpnK/kbr9wARQRZffAxuJ5zDYwJvxPFNoPDYWvPcsvPwQS4HgJ81pgMtGeZ9A6Z0ukKqif8HillGBSxA903SyyEGAXuAX5wkzQfFkJcEEJciEQiVd5edcjpRGBQp5Dr10hJD61B85yMU2g3tsxsTFPEjqfLTEQwql4YWGYAguYY5ACcnUMUcRg+qykHGu9opsMV1Mps2MCurlswgXPdXtiX8IUQXxVCvLbL35OHyUgI0Qp8FvinUsrkXumklJ+QUp6XUp4PhRos59KdZAzqFIXINeZkD/0djVcx6HAH+ingNKxTlKLTLFRCdAcaFz1xJwLBThK0GlZmuTHNhmwlELypPVRX9AZaWZRdSMNmcrNqYA+aZ2DvDbSwIEPGOV/pGnwTeM3vBcd+CaSUj+31nRBiVQjRJ6UMCyH6gF3PxxNCOFFk/8dSyr848t3WGcGODqLzfrqMaiDr15iVPfT7zUP4vYEWFuhmwqBBrhy9xrzJytzn9zBb6eHOjRmMmIgXI1PMyV5TqFV09Gly1MHoNE4Drl9eVwN7T6DFgKsfDX0BFRF2aH3GEHmi3JhmVQbpCt6gXTENqi3354EPaa8/BHxuZwKhFrN+D7gkpfwvVeZXV/QHvMxWeigbcXBCuYRzc4E52UufiYig1+9hptxtjC5dSuyJWWZlT0PD5e6E2qDvpmzgVN9sA7vuc2GLzxoSKbS0PmMaSaYOvcwiPmNImYsRZcCZVaED1RP+x4B3CSEmgXdp7xFC9AshntbSPAz8feCHhRAva3/vqTLfuqC3XdPiG+F2n1jAJkvMyF5TEYHySOxRy1i17hSZdRzFlCkt/DnZgz25WHs5aimPM7XEbMVcA7vuc2EvGKBCkxJ7fIYZ2WsKpysdqszdOAqbhijvRGyGuUoPfSZq2zux75LOzSClXAce3eXzZeA92uvvgCEzZcPRF/DwnOzGnvoulPLgqKGGWhtEwrY+Aj4jJtVHQ2+7l5dlD7ZSBtIRaK3hurNmQc/KHnr8jdej6+jze/mi7EboctRaqqZicwik+Sx8LdQAoOrFF6zdxdMRHMUUs7KXHzLRIKf7XABqg76WZc6ncGYjzMkebjvGFv6xRr/fy1ylB4GE+HxtL64Rfr5txFQSLt0jEaj9JqZW5qR3ELej8dpsHeq0L40Ial5mFaYibDfXwN7Z4mJZ9Kk3tVYnrasyz5rE6UpHu9fBql0rc63rOXbdmDnOSzrHGvopUIAh5JfDjavDPJJMgJ52z5utoFoiNkMFQdlvDqcrHQGfkxWH7nxlzCBXaB811cAuhLguwTVoYI95hk3hdKVDCEFJb3u1FiVoZV62mdfpCpqEf1N4nHaSXs0RyQAiWBC99Ad8tb1ulfA47WS8A1QQhhBBVHTSFWh8KOjtEEJgb++nKAyQo25MkxKttAQaHzFyJ4IBPxu2oAFt+xplbEi/eSSZOjoCATZsHQaUWdfgD5tqYN+JJuHvA1d7DznhqTkRyPVppkrdpgiLvBOdgTZi9i5DOoVZp7y9AR8rtl5DLL95eugz2cAOSoW2QE/tB7n1a6zYekzla6Gjz+9lUXZDrMYxsjamSdj8tPk7a3vdGqNJ+Pugv8NLWNSYCCpliM8qbbYJya/P72VR9NZ+kNuY5lqp25Qqhq3luxoTvl5mM9Zzr9/DtVKo9hFhN64xU+k1ZT33+T1MlUK1V95tTDNPrynCX98MTcLfB31+LzOVUG2t3eQSolxQ1q4JG0ifRgQ1Jb9cEpGJmiYO/k70+71MFjXP01rJUUsFiM8zLc1JBP2az4XYXIZitjYXlRK5Ps1kqdtUvhY6+jQhBsllpbyrEeTGDNdKIVOWeTuahL8Pev0epkpaXPxKpTYX1VQMc7LHlJZfX8DDVDEE6TXIp2pzUW3AnDOZ05WOXr+HuUo3ophWctRaILGAkBWlzTYh4fdqEWGB2i1xpFYRxbTyLzGRJFOHfs5FTZV3xRwkl5gpm7M/b0eT8PeBfjSaKOUgtVKbi2rTyZlKrymJoG+7OqlWVn50EoAZ2ceACcus5Kg1VmRp9Txr2qU7Dwu6BLdWM9htZR4w4b5FnxH1HJvZ8rXoNeEy1nY0CX8f6N55QE2JoChcZD0hWt1V+b4ZgjeVuVZEsH4NiWCOXnNa+O3byzxbm4tq5Ddn0qW7/oC39uSnzV5Na+G3e7cNcrO1uahmzEzLflMaM9vRJPx9YIi1uzHDmqOfPhMFltoOPeYIUEMimCTm7CHQ3o7Tbr5m1x9Q1q5E1HSQy9l8FDxBUw7sHT4naYefvM1XQwv/GmXhYIUuek3kdKWj3esg7QpSsHlrWM/67LWXARNFvt0N5ut5JkOv38OS7KKCvaYNZA7znnvZ0+4hSSs5R3vtyhydZME2YMrNS1BH4NmcbpKu7poOcivOIfr95lvaAOV/0Of3subsr6mFv+7sI9TegsOEA7sQgl6/l1XnwJZlXjWiU2w6uxDudvxe83hT7wbz1YjJ4HbYCbT6iNeKCMpF2JjmcqnPlNN8UM5XnS0uorUiAilhfYrJcq9pp7xCCPr9XlbtvTUd5KZlv6mCpu1Er9/DIj01XcNfEH2mtnT7/V5mGIDo1dpccH2SJbt5jZntaBL+AdDr9xC21YgIYrNQKfFa3tw7+gMd3toRweYKFFK8ng+ZulP06ZFCazHIFdKQWOBS0Zx6dB39fi/XSpojUqVc3cUqFdiY5lq5x9z17PdwpdSrVDrVylGlVAN7xdyDnI4m4R8ASovfC+tT1Wu0Navimuw3NREMdniZKocgvgDlUnUXW58CYLJs7k7Rp2nxlRx1s7qLaWV+Nd/DgIkt/L6Ah0v5TqgUlTa9GiQXoZjhlVyPaWdyoAj/Yq4bkFv1dGRk1iEX5/WCucuso0n4B0Cf38NrhR7IJSC166FeB8e2HX0zW0EDAS9v5LpAliFepUZb29SarvSZmvwGO7y8ktM2q6td392q5z5TD+y9fi8zFX2Dvkrv04gyZq6Wzd22Bzt8TFYG1Jtql3W0en690G3qMutoEv4B0Of38lpBC6savVLdxaKTZN0hNvGZUramYyDg5VJRL3O1nWKKst1DmKAptdk6BgLe60QQqb6epbBpkkzz1nO/38O1ihaxtep6Vs9sSvabeiY30OFlRvYqRVa1A/v6NkmmicusoyrCF0IEhRBfEUJMav/3PMxRCGEXQrwkhPhCNXk2Av0BD1N6p6iaCK4S9agQrWbUo+sY6PAxJXXyu1zdxdYniXmHkdhMPcgNdviYkz1UhKMGA/tVUt4B8rgYNPEg1+f3skoHRUdr9W07coW8q4MY7aZe3hjs8JLHRdo3UJOBvWJzsihDpi6zjmot/I8Cz0gpTwLPaO/3wi8Cl6rMryHobfewQpCSo6U6K0hKiF5l0TZIT7vbVIeA7MRgh5dNfOQ83TXpFGHHIO0eB20e88rWBju8lHCw2TKytTxxZEQnibhHsAlzD+xKGiyItYxXP7BHr7LuHQUw9fJGn9+LELDmHqmBhX+NpG+YCrZbgvCfBD6lvf4U8IHdEgkhBoH3Ar9bZX4NgWq8gkTLWHXkl45CLs5kpY/BDvNafcDW9DTqHauOCEp5iM8xK/sYMHmZe/0ebAJWXSPVlblSgfVJ5kU/ve0eXA7zrpwGfE48ThvLjqGaWPhLjiH8XqcpHc10uBw2eto8LNgG1JJMNTGy1idZdQ7htAu628xzbOdeqLYl9kgpwwDa/70OQP114F8C+z5ZIcSHhRAXhBAXIpEaBbGqEvoxbavu0eo6hTY7MLuKAaDd46TN42DBMays3aOqk6KTICu8Vuw39YYtgNNuo7fdw6xtUMlRjxpNMbEApRxXTK5KguvOV9NiSKmTMhtHu1A6CtkNpirmDy8AajZ3pdQHpRwkjhhErVSAjWnmxAB9fi82m3kPPtGxL+ELIb4qhHhtl78nD5KBEOJ9wJqU8oWDpJdSfkJKeV5KeT4UMscpQS6HjVCbm1kxoAKo5RJHu5BG+M9vdjJociKAbZuYxTQkFo92Ec1Sfj5j/kEO1Dr+pVIfyMrRJXvaMsEr2W7Tz+RAkd9rBe2Ix6MaNNrvXiuYMxT0Tgx0eHm5WkXW+hRUSrxeGjT13tR27Ev4UsrHpJS37/L3OWBVCHUSsvZ/N83iw8D7hRCzwKeBHxZC/FENy1AXDHV4eUNXrRx1fTd6lYrDy0IlaBEi8PFKvkoiWLuEtDl4NddtemsXFPm9lNFiJx2Z/LRBLtVliYF9sMPHcymN/I66lKVtcj+fClmkzF6e3+xSb466L7f2BgAvZHtMrT7bjmqXdD4PfEh7/SHgczsTSCn/dynloJRyFPgg8DUp5d+rMt+6Y7DDxwsZbcZxVAXH2htk/CeR2CzTKZ7TO8VRiWDtEgX/GEUclrH8nkt1apK9oxNB2RciUmmzTD2/nmlHOn1VDHJXkU4fk3m/JazdgYCPSKWNsidYnTEj7FxIdVrCmIHqCf9jwLuEEJPAu7T3CCH6hRBPV3tzZsJQ0MsLST/S7jo6+a2+TtQ3AWAJIhgIeFnI+6j4uo5e5sgl4q0ntq5ndgx2eMlUnJTbh6uq5832UwCWsPyGgj4kNvKBE0c3ZqJXyPknkNgsUWa9/2Xax6uayZUC4+Sk0/T7UzqqInwp5bqU8lEp5Unt/4b2+bKU8j27pP+GlPJ91eTZKAx1+ChUbBQDE0drIKk1SEeYd44B5pat6dCtlqz/xNGs3UIGNmYIu0bV9axQZo2sNtuPWM+VMkQus2ahgV2/x5ivChXa6uvEWlSZrWDh62Ve851QSzNHESWsvUGy/SRgjYEdmp62B4a+5p5oPwUrrx3+AquvA3C5MkSozY3HaV4Nvg69U2z4xmDt8uE7RfQqIJm2DeGy2+hqNb9sTS/zqmdCbeaVCoe7wMYMlHLM2kcRAlN72eoY0tr2smsEkkuHFyWkIpBaZcmjCN8Kyxu6wTXnGIN88vDhQzRjZtWtDDgrlBmahH9gDAVVhS57TsDmMqTXD3cBbYPnpXy/Jaw+uG6RzzsnIJ84/Bmga8rP7vXiAP0BjyVkazpBX7OPqYBih13WWVMD++vlIXraPKZ2rtPR1erC47QxKZUHuG6cHBirrwIwyaga2FvMP7B7nHZCbW4uVUbUB4cts2bMzNjVMzPr2RY70ST8A6LP78UmYFKoEV1v5AfG6hvQ0s0bSbclFDoAwRZFBJcZVR+sXDzcBdbeALuLl1JBy1hAboednnY3r5U1Ilg5Qj0LGy9mraFKAqXFH+zw8VJRI/zwIetZm/G+VBxkMGgNPToog+aFXC8gDj9r14yZN0qDdLVaY8YOTcI/MFwO5ZTzSnFIfXBoIngN2X2O5XjWMha+TgQv5/tB2A5f5pWL0H2OuVie4aA1BjlQRPBKugOcLUeqZ4LjzCYqlqln0LT4SS+0hA4/sK++Bm19XEo4LVXPgx1ephNAcOzwBtzKq+BQ7WQ4aJ16bhL+ITAY9HF50wVtfYezCLSNvHTHaYplaSkiGAh4mU1K6Dp1OMtPSghfpNhzJ+vpAsNBc57fuxsGO3wsJgrQc+7whL/2BpXQOcLxnKXqeajDx0IsC713Hs3C77md+fWMxQjfx1I8i+y5/fAWfvgV6L2dmQ1rGTNNwj8Ehjp8LGxkoef2wxHB+hSUcqz51I6+VZZ0QDv5KpaB3jsOZ/klFiG7QbT1NIClOsVAh5fleJZKzx2qng+6WZ1LwMY0qY6zlCrSUvU82OElmSuRC92u9i0OulldyitJZuc5krnS1gawFTDQ4aVYlqQCZ1UojXzqYD+sVGDlIuWeOwknspZq203CPwQGO7ysbuYodd+u9MoHjbWy9CIA005dwmUdy2+ww0ssUyQful0pOA66WR1+BYA5lyqzlTrFYIeXUkWyGThzuM1qrczLvrOAtep5SKufSMtJbbP6gIFtI1egUtoyZoYsVs8Aq5qE9sAbt7EZyCeJ+c9SkdYqc5PwD4GhoA8pYaPtNFRKW8qbfbH8IjhbuFxUYQqsNNUf0ZZiVjyqQx/Yyg+/AsLOpYra87AS4etEvehWDmMHLrM2sE/aNUczC9Wz3iZnnVqZD7qsE34ZgGnHOGCteh7U6vmaQ2vbyy8e7Ie6MeNWznVWKnOT8A+BrU7hVssULB0oHpwigv67mY/nLbWjDzDSqRrzpE1TJ2mNfV+EX4HQaWYSFdo9Dvw+88bB34mRTjXIXWEEbI4tIt8Xyy9BYISptBshrDWw60sxVwohcLUevJ4XnwdPgEsFFYtnyEIbmPqS29VMG7T1w+KFg/0w/DLYXVwuqwOChjubhH8soU/dpgqdSs2w8Pz+PyoV1Dpw/z3MrWcYtVDjgOuEP5V2Q3BcdfD9ICUsXYD+e5lbz1iqQ4Cy8O02wUy8rPYuDlJmUBZi/z3Mrafp93stocHXEfA5aXHZWYznof+eg5d58QUYPM98LEewxWXqA252wutSEtzZ9QwMvkW12YNg6UXouY25eHErtr5V0CT8Q6C33YPDJliIZ2Hw/oN1irU3oJyHAUV+uvVoFbR5nARbXMytp2HoAVj4wf6bmOtTkFmH4QdY2LCWcgOUBHcg4NWI4H41kyuXbv6j9Lpa6x+4l7mNzNZAaRUIIRgK+tQG/dD9ykgppG/+o/ymat+D97GwkbHUWraOkc4W1bYH74PYrIrrfzOUi2omMPQg8+sZhjqs43cATcI/FOw2QX/Ay8JGBobug41r+29iLjwHQL7nXlaSOctZ+KCs/Ll1jQjSEbVpdTPMPwtAefABFmNZixKBTxvk7odiZsuDdk/Mf1/9H7zPkgM7qCWO+Y0MDD0Isrz/kuXSi4CEgfPMW3BgBxjrbFED+8B59cF+ZQ5fhFIWhh+0ZJmbhH9IbJHf4H3qg/2s/Nlvg3+YuYoKM2y15Q2AkaBO+A+oD7RBbE8sPAveDladgxTKFct1CoDRzhZmo+nr9bxfmee+Cw4Pyc472UgXLGfhA4xqbbsyoJV5/gc3/4HW9kt997AUz1rKAUnHSJePaCpPqvN2tV+zsE+ZF5QxI4ceaBL+rYCxLkUEsv8esDlh7jt7J65UFBGMvk2RB4pIrIaRzhaWE1nywVPgbr9uze6F+R/A0APMx3KAtVQMOkY6fSRzJeKuXmjt3b/Ms9+GwfuYT5QBLDmTG+1qIV+qsFL0QOjsFrntidlvQ/dthAteyhVpyXrW++PcJsrKn/7GzX8w/33oGCVm7ySVL1lu9tok/ENitLOFzXyJ9YIDhh+Ea1/fO3HkslrLHn2bspCxKuErOepCLA8jb1WdYq91/GRYHQw9/JBaHsCahK/X0+xGFsbfqeq5Ut49cTamPDW31bOVPIt1jHVpZY6mVT3PP7u3A1Yxp74ff8dWPVuN/OC6KGFuPQPj71BKq2xs98SVCsx9/01t22pLd03CPyTe1ClOPKriiGyu7J54VrP+R9/G3EaagM9pKXmiDr1Rz2+k4cRjanNrr/Nerz2j/p94lIWNDDZhjdj/OzHapROBVubsxpbm/AbMfR+QW/UMWHNJR2vbM+ta2y6kYHGPpazF59QB4GNv3zbIWa/MetueiabVwC4r1/vtToRfhkwUJn7YssZMVYQvhAgKIb4ihJjU/nfskS4ghPhzIcRlIcQlIcRD1eTbSGx1imgaJh5VH+5l5U9+CTrGoGNEbeRZrHHo0MlrNpqBk+9SH05+ZffEU1+F1h7ouZ3paJrBDh9Ou/XsisEOH0JoZZ74IUDA1Nd2T3z1i+BqUxu20QxdrW5a3I663m8t0Nfuwe2wKWNm9BG1pj311d0TT35FLWmOPMx0JIXbYaPfb72BvdXtoKvVrQb2gfMqYN5e/Xnqq4CAiUeZ05ZoreR3ANVb+B8FnpFSngSe0d7vht8AviilPAPcBRzQb9t8GOxQGu3Z9bSKqdPSrYh9J3IJmP4mnFUHfFlVuQHQ2eKi1e1QVk3HqAqkNvnlGxMWc4oUTzwGQnBtLcVEyJpl9jjt9Pu9ighauqD/7t3LXCnDlafVQOhwM7eRtuT6PYDNJhjp9DETzYCnHYYfgstP37h8JyVc+p9qCcTTznQ0zVhXi6Xkidsx1uVTSh2HS1n5V55Wyzc7ceVvYOBeaOnkWiRFv9+Dz2Wtgb1awn8S+JT2+lPAB3YmEEK0A28Hfg9ASlmQUsarzLdhcNptDAd9yvKz2eDc+1VDyCXfnPDql1RMkjM/SqFUYTFmPacrHUIIVeZ1TZd96nG1YZeKvDnh5JdV7Jnbf5xKRTK7nmY81Fr/G64RRjq3lfnsj6pljI0dktT5Z5VU9cx7ASzpaLYdo50t18t824+pmFE7AwWuvqakuWffD8B0JMW4RQd22KbFB7j9x2EzfOMmfXRKOdad+wAA01Frtu1qCb9HShkG0P5375JmHIgAvy+EeEkI8btCCOu2DpQCY1qb0nHnB9Va5hufe3OiF/4AAiPKKSWWoSKtt8GzHaNdvi2lEXf/XRVL6OKn35zo1T9VHshj72Q5kSVXrBwDIlBrtdz5twEBFz/z5kQv/IFSLp16nHS+RDiRY8KCRKBjrKuF+fUM5YpUhG9z7FLmT4HdBWfeS6FUYSGWZbzLumUe7fSxmsyTKZSUMePwwit/8uZEFz+jzoS44yeRUlp29rov4QshviqEeG2XvycPmIcDuBf471LKe4A0ey/9IIT4sBDighDiQiQS2StZQzHapSwCKSUMnofQGfj+b19Xcay8puSY9/1jsNmYWlNhV090W7dTnAi1Mr+RIV8qQ/cZpU+/8MnrHqixWbj813DXT4PdwXREDQ5WJr/RTh/r6QLJXBH8g2ot//nfu+6BmlqDN/4K7voguFu3ldl6RKBjtKuFQrnCcjwLvqCa2bzwKchsqATZOLz8/8DtPwEtXcxvpClXpOUHdkAtWbpbVX1e/IxSnIE6v/aF31d7du19rG3mSRfKTFiwP+9L+FLKx6SUt+/y9zlgVQjRB6D9X9vlEovAopRS92j4c9QAsFd+n5BSnpdSng+FQocvUR0w1tVCplBmbTMPQsA7flmFk335jxXpf+lfgdsP9/x9gC3Ct2ID0THR3UpFapvVAG/7Z7AxDc99Qq3pfulXlNX34M8DcC2iymxlItCn7Hr98c7/HdJr8N3fUO+/8m+UquOBp1S6yCZg7YF9S46qL3G8/V9CYRO+9Z/U+2/8e+V5/OD/AsA1bZCz4vKGDl15d21NK/PDv6hmsF//d+r9d/6rWrZ75Je0dFrbtuCsptodh88DHwI+pv3/3M4EUsoVIcSCEOK0lPIK8ChwwLjC5sRWA4mk6Gn3qHW9kd+Fv/7nygIMvwzv+6/KQkI1kN52D60WVG7oONndBsDkaoozve1w+j1w8t3w5V+F1/5cuaQ/9v+F9j4ApiNp2twOQq3mP9B6L5zqUR16cnWTe4c7VJiFO34Kvvkf4drXlKfpI/8COlU89Wtraew2YUkNvo6JbnXvU2spHjkZUqd+nf/H8OzvqBgyi8/B/R+GvjsBtmY1Vh7YJ0KtCAGTa5tAnzry8K3/K3z315UvzeLzaul25K0AXNOMHv1ZWQnVruF/DHiXEGISeJf2HiFEvxDi6W3p/gnwx0KIi8DdwP9ZZb4Nxame6+QHqM3bv/1HavpbzMLjH4PzP7OV/lokZWmrD1SHtgmY1K1dIeBv/a6a/uY34Yf/tbKMNExHU4x3tyKENZUboEIGe5w2rur1DPCjvwHn/5Fa2njkX8AP/autr6bWUowEfbgc1pOh6gi1ugn4nG8u8+P/Xs3csjF46CPw7n+/9dV0JEVXq5t2C0XJ3Amvy85Qh+962wbVnh/556qe7/tZeP9vbn11bS2Fz2Wnt906UTJ1VGVySinXURb7zs+Xgfdse/8ycL6avMyE7jY37R4HV1c3r3/oC8JPfPKGtFJKrkXS/MRbBut4h7WHx2lnOOjbms6qD9vhA7+za/rpSJqHxjvrdHfGwGYTnOhufTMRuHxq9rYLrkVSll62A6XIOtXdxuT2tu1wK9J//N/fkH46mrb0noWOk92tTG0f5OwOePTX1N8OKIVOiyWNGeuaIg2EEILTvW3XLfybYCWZI5UvHYtOcaK7TZv23hy6WsXK03wdN5DfHiiVK8yupy29Sa3jZI8a5OQ+YbCVMZOy9Pq9jpM9bUxHU5TKu+jvd0ApdKxZ5ibhHxEne9q4srq5b6fY2rC1aAPZjpM9rcxE0xT36RSXVxRB6ktfVsaJnlbCiZxS6twEs+sZimVp+aU7UPWWyBaJbN78zOaVZI54psjZPuvX88nuVoplqRywboJUvsRSPMsJi/bnJuEfEae6Ww/UKS6HFfmd7rV+pzgRUp1ibp9OcSmsnNDO9rXX47YMxanuHfs1e+B6ma1fzye1zeqrByzzmV7r17Ne5ql9ZrCXLd62m4R/RJzSCHy/TvH6coLedg+dFlar6Di5TbVyM1wKJ2nzOCx1puteuL5Bf/MyvxFO4rSLLTWTlaGX+eq+9ay+P3MMBrkT3Uqpc2XlYIPcuf4m4d9SOK11issryZumeyOc5DaLNo6dONndht0meCN88zJfCic529tuyU2tnRjs8OJ12vcd2N9YTnKiu83SCh0dXa1uulpdB6rnwQ6vpRU6OnwuB2NdLby2nLhpujfCm/i9Tvr81lPoQJPwj4zOVjd9fg8XF/duILlimWuRtGWtgZ3wuuyc7G69aZkrFcnllc1jsbQBSqlzpq/tAESQ5JxFp/m74Y4BP6/epJ5BEf5xWM7RcceAn9eWDlbPVjVmmoRfBe4c9HNxMb7n91dWNilX5LGx8OF6p9hrs3p+I0OmULbsGuduuGswwKuLiT0VHGubOSKb+WMzsAPcMRhgcm1TxZfZBblimZlomnPHZGAH1bbDidye+3LliuTKStLSbbtJ+FXgzsEAs+sZEpndFRz6lPi2fn89b8tQ3DnoZz1dYDmR2/V73RI+TuR391CAbLHMVGT3ZR19Lfs4Wfh3DvipSHh9efdlndeXk1QknDtGbfv2AVWWvaz8mWiKXLFi6dlrk/CrwF2DAQAuLsV3/f7iYpz2Y7J5qUPvFK/uMbN5cS6Ox2mztBW0E3cOqjK/shDf9fuX5+MIAbcNHL8y77V899K8Ogbw3pFAvW7JcOgz8Vf3IPwX5+IA3DMcqNMd1R5Nwq8Cd+zTKZ6b2eD8aNCy63274WxfOw6b2LPML8zHuHMgYMlTrvbCaGcL7R4Hr+xR5gtzG5zpbT8Wm5c6uts99LZ79lyyfGEuxmCHl+42a25e7oY2j5PxrpY9y3xhboOAz2nJoGk6jk+vbAD8XifjoRZenLvx0OP1VJ5rkTT3jQYbcGfGweO0c9uAn+dmNm74Llcs88ZygntHdj3p0rKw2QR3DQV2redSucKLczHOH7Myg7JkL8zGbtivkVLy4nxMBZQ7ZrhvNMhzMxvqPIAdeGEuxluGOyx7shc0Cb9qPDTeybPT6zd4n17QyOG+0ePXKd460cnLC3HS+Tdv6L00H6dYlrzlGJLfg+OdXF7ZJJp684be5ZVN0oUy549jPZ/oYimevcHRbnY9w2oyf0zL3EkyV+KNHXsXugH3FouXuUn4VeKRk12kC2Vemo+/6fNvXY3gc9m3ln2OEx6e6KJUkTw3+2Yr/+tX1nDaBQ9NWDto2m545GQXAN+dir7p829eVYf0WD1Q3G54WKvH7+wo8zeuqGMv3nlqtwPurA29Hr8/vbPMqp4fOWHOMzoOiibhV4mHJrqwCUXwOioVyVcvrfL2kyHcDnsD784YnB/twO2w8fXLbz7v5uuX13hgrNPScf/3wm39fjp8Tr555c2nsH3ljVXuGgrQbcFQufthrKuFfr9ni+x0fONKhPGuFkuf3bsXuts9nO5p46uX3ty2n7m8Sneb2/IS6ybhVwm/18lDE5381ctLVLR1v4tLCVaTeX7ktp4G350x8DjtPHauh7++GN5aypqOpJhcS/HO09a2gPaC3SZ49GwPX35jdUubvpbM8fJCnHedPX6WLqiosE/c0cc3r64RSxcASGSKfH96nXeePp5lBnjvnX08P7tBOJEF1N7Ut65G+eEz3ZZev4cm4dcEP/mWIRZjWZ6dWQfgM8/P43bYePTM8SR8gA/cPcB6urBl8X7m+QXsNsH77+pv8J0Zh586P0QqX+LpV1cAVWaA99zR18jbMhQ/fu8AxbLkf15cBuCvXl6iUKrw4/cONPjOjMP77+pHSvjLl5YA+MLFMKl8ifffbf223ST8GuDdt/XS4XPy61+d5OrqJp99cYkfv3cAv+/4yPR24h2nQgwEvPz6M1eZW0/zh8/O8fjtvcdyaUPHfaMdnOhu5Xe+PsVSPMsnvzvD20+FjkU8+L1wrq+duwb9fPwb11iOZ/mdb0xx73Bgyx/jOGK0q4VHTnbxe9+eYSme5be+NsmpntZjsU9TFeELIYJCiK8IISa1/7tuYQsh/pkQ4nUhxGtCiD8RQhwrVvC67Hz0iTM8N7PBj/zXb+F12vlnj51q9G0ZCpfDxq+89yyvLSV55//1DQTwL999utG3ZSiEEPzr951jOprm4Y99jXS+zK++92yjb8tQ6GVe28zz8H/4GhvpAr/y3nONvi3D8dEnzpDKl3j4Y19jbiPDr73vtmPhTyP2O8Djpj8W4j8CG1LKjwkhPgp0SCl/eUeaAeA7wDkpZVYI8afA01LKP9jv+ufPn5cXLlw48v3VG599YZEfzKzzc4+Mc/IYHP5xEHz+lWW+eSXCh946wp2a5/FxxzOXVnn61RX+zgPDx1KCuhu+dTXCX728xI/dM6AON78F8PzsBn/x4iLvOtfDD1toeVYI8YKUctcjZasl/CvAO6WUYSFEH/ANKeXpHWkGgGeBu4Ak8FfAb0opv7zf9a1G+E000UQTjcbNCL/aNfweKWUYQPt/w9a9lHIJ+L+AeSAMJG5G9kKIDwshLgghLkQikb2SNdFEE000cUjsS/hCiK9qa+87/548SAbauv6TwBjQD7QIIf7eXumllJ+QUp6XUp4PhW6NqWMTTTTRRD2wr4eMlPKxvb4TQqwKIfq2Lems7ZLsMWBGShnRfvMXwFuBPzriPTfRRBNNNHEEVLuk83ngQ9rrDwGf2yXNPPCgEMIn1Db3o8ClKvNtookmmmjikKiW8D8GvEsIMQm8S3uPEKJfCPE0gJTyB8CfAy8Cr2p5fqLKfJtoookmmjgkqlLpGI2mSqeJJppo4nAwUqXTRBNNNNGERdAk/CaaaKKJWwSmXtIRQkSAuSP+vAuI7puq/mje1+HQvK/DoXlfh8NxvK8RKeWumnZTE341EEJc2Gsdq5Fo3tfh0Lyvw6F5X4fDrXZfzSWdJppooolbBE3Cb6KJJpq4RXCcCd+sWv/mfR0Ozfs6HJr3dTjcUvd1bNfwm2iiiSaaeDOOs4XfRBNNNNHENjQJv4kmmmjiFoGlCV8I8bgQ4ooQYko7cWvn90II8Zva9xeFEPfW6b6GhBBfF0Jc0o52/MVd0rxTCJEQQrys/f1ane5tVgjxqpbnDXErGvHMhBCntz2Hl4UQSSHEP92Rpi7PSwjxSSHEmhDitW2fHfQoz5u2RwPu6z8JIS5r9fSXQojAHr+9aZ0bcF//hxBiaVtdvWeP39b7eX1m2z3NCiFe3uO3Rj6vXbmhbm1MSmnJP8AOXAPGARfwCuoYxe1p3gP8DSCAB4Ef1One+oB7tddtwNVd7u2dwBca8Nxmga6bfN+QZ7ajXldQziN1f17A24F7gde2ffYfgY9qrz8K/IejtEcD7utHAIf2+j/sdl8HqXMD7uv/AP7FAeq5rs9rx/f/Gfi1BjyvXbmhXm3Myhb+/cCUlHJaSlkAPo06aGU7ngT+h1R4FghocfsNhZQyLKV8UXu9iQoHPWB0vjVCQ57ZNjwKXJNSHtXDuipIKb8FbOz4+EngU9rrTwEf2OWnB2mPNb0vKeWXpZQl7e2zwGCt8qvmvg6Iuj8vHVqY9p8C/qRW+R0UN+GGurQxKxP+ALCw7f0iN5LqQdIYCiHEKHAP8INdvn5ICPGKEOJvhBC31emWJPBlIcQLQogP7/J9o5/ZB9m7IzbiecEBjvKk8c/tZ1Azs92wX50bgY9oS02f3GN5opHP6xFgVUo5ucf3dXleO7ihLm3MyoQvdvlsp8b0IGkMgxCiFfgs8E+llMkdX7+IWra4C/hvqMPd64GHpZT3Ak8AvyCEePuO7xv2zIQQLuD9wJ/t8nWjntdB0cjn9itACfjjPZLsV+e1xn8HJoC7UedY/+dd0jSyb/40N7fuDX9e+3DDnj/b5bNDPTMrE/4iMLTt/SCwfIQ0hkAI4URV6B9LKf9i5/dSyqSUMqW9fhpwCiG6jL4vKeWy9n8N+EvUNHE7GvbMUB3sRSnl6s4vGvW8NKzqy1pi76M8G/LchBAfAt4H/F2pLfTuxAHqvKaQUq5KKctSygrwf++RX6OelwP4ceAze6Ux+nntwQ11aWNWJvzngZNCiDHNMvwg6sjF7fg88A805cmDQEKfNhkJbY3w94BLUsr/skeaXi0dQoj7UXWxbvB9tQgh2vTXqE2/13Yka8gz07Cn5dWI57UNBznK8yDtsaYQQjwO/DLwfillZo80B6nzWt/X9j2fH9sjv7o/Lw2PAZellIu7fWn087oJN9SnjRmxE12vP5Si5Cpq5/pXtM+eAp7SXgvgt7XvXwXO1+m+3oaaal0EXtb+3rPj3j4CvI7aaX8WeGsd7mtcy+8VLW8zPTMfisD92z6r+/NCDThhoIiyqP4x0Ak8A0xq/4Na2n7g6Zu1R4Pvawq1pqu3sY/vvK+96tzg+/pDre1cRBFSnxmel/b5H+htalvaej6vvbihLm2sGVqhiSaaaOIWgZWXdJpoookmmjgEmoTfRBNNNHGLoEn4TTTRRBO3CJqE30QTTTRxi6BJ+E000UQTtwiahN9EE000cYugSfhNNNFEE7cI/l856Yyt0jkevAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(t,ans1.T[0])\n",
    "plt.plot(t,ans2.T[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_xy(theta):\n",
    "    return x_f(theta), y_f(theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1, y1 = get_xy(ans1.T[0])\n",
    "x2, y2 = get_xy(ans2.T[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAM7UlEQVR4nO3cX4hc93nG8efRn6C46+ILb2lq/TPUK2LU1KmEm+KLalW3KIlRaGjBRjGFtuxNDS6kJHEFDaEYCoE0F3UpIjYpRI1ZSEKDkuC4eFUTaBJLju1KlW1MaBw1KSaEkCwCB1dPL3YMsrSrmZ357Z45b74fGHZn9sxvXlbar47OnD1OIgBAHVu6HgAA0BZhB4BiCDsAFEPYAaAYwg4AxRB2ACimWdhtb7X9HdunWq0JAFi/lnvsD0q60HA9AMAYmoTd9k5J75f0mRbrAQDGt63ROp+W9BFJN661ge0FSQuStGPHjgO7d+9u9NIb5/Lly9qyZfrfhmDOdvowo8ScrfVlzpdffvlHSWaHbphkopukeyT94+DzQ5JODXvO3Nxc+mBpaanrEUbCnO30YcaEOVvry5ySzmSELrf4J+ouSUdt/7ekxyUdtv25BusCAMYwcdiTPJRkZ5K9ku6V9FSSD008GQBgLNN/UAkAsC6t3jyVJCU5Lel0yzUBAOvDHjsAFEPYAaAYwg4AxRB2ACiGsANAMYQdAIoh7ABQDGEHgGIIOwAUQ9gBoBjCDgDFEHYAKIawA0AxhB0AiiHsAFAMYQeAYgg7ABRD2AGgGMIOAMUQdgAohrADQDGEHQCKIewAUAxhB4BiCDsAFEPYAaAYwg4AxRB2ACiGsANAMYQdAIoh7ABQDGEHgGIIOwAUQ9gBoBjCDgDFTBx22ztsf9v287bP2/5Ei8EAAOPZ1mCN1yUdTrJse7ukb9j+WpJvNlgbALBOE4c9SSQtD+5uH9wy6boAgPF4pcsTLmJvlXRW0q9LeiTJR1fZZkHSgiTNzs4eWFxcnPh1N9ry8rJmZma6HmMo5mynDzNKzNlaX+acn58/m+Tg0A2TNLtJuknSkqT919tubm4ufbC0tNT1CCNhznb6MGPCnK31ZU5JZzJCi5ueFZPkJ5JOSzrScl0AwOhanBUza/umwedvl3S3pBcnXRcAMJ4WZ8W8Q9I/D46zb5G0mORUg3UBAGNocVbMC5Le3WAWAEAD/OYpABRD2AGgGMIOAMUQdgAohrADQDGEHQCKIewAUAxhB4BiCDsAFEPYAaAYwg4AxRB2ACiGsANAMYQdAIoh7ABQDGEHgGIIOwAUQ9gBoBjCDgDFEHYAKIawA0AxhB0AiiHsAFAMYQeAYgg7ABRD2AGgGMIOAMUQdgAohrADQDGEHQCKIewAUAxhB4BiCDsAFEPYAaCYicNue5ftJdsXbJ+3/WCLwQAA49nWYI03JH04ybO2b5R01vaTSf6rwdoAgHWaeI89yQ+TPDv4/GeSLki6ZdJ1AQDjcZJ2i9l7JT0taX+Sn171tQVJC5I0Ozt7YHFxsdnrbpTl5WXNzMx0PcZQzNlOH2aUmLO1vsw5Pz9/NsnBoRsmaXKTNCPprKQPDtt2bm4ufbC0tNT1CCNhznb6MGPCnK31ZU5JZzJCj5ucFWN7u6QvSDqZ5Ist1gQAjKfFWTGW9KikC0k+NflIAIBJtNhjv0vS/ZIO235ucHtfg3UBAGOY+HTHJN+Q5AazAAAa4DdPAaAYwg4AxRB2ACiGsANAMYQdAIoh7ABQDGEHgGIIOwAUQ9gBoBjCDgDFEHYAKIawA0AxhB0AiiHsAFAMYQeAYgg7ABRD2AGgGMIOAMUQdgAohrADQDGEHQCKIewAUAxhB4BiCDsAFEPYAaAYwg4AxRB2ACiGsANAMYQdAIoh7ABQDGEHgGIIOwAUQ9gBoBjCDgDFEHYAKKZJ2G0/Zvs12+darAcA13XypLR3r7Rly8rHkye7nmiqtNpj/6ykI43WAoC1nTwpLSxI3/uelKx8XFgg7ldoEvYkT0v6cYu1AOC6jh+XLl1662OXLq08DkmSk7RZyN4r6VSS/Wt8fUHSgiTNzs4eWFxcbPK6G2l5eVkzMzNdjzEUc7bThxmlX+w5f/fwYXmVbsXWvz/11Fhr9uX7OT8/fzbJwaEbJmlyk7RX0rlRtp2bm0sfLC0tdT3CSJiznT7MmPyCz7lnT7JyEOattz17xl6yL99PSWcyQmM5KwZAvzz8sHTDDW997IYbVh6HJE53BNA3x45JJ05Ie/ZI9srHEydWHockaVuLRWx/XtIhSTfbvijp40kebbE2AFzj2DFCfh1Nwp7kvhbrAAAmx6EYACiGsANAMYQdAIoh7ABQDGEHgGIIO4DNx9UZN1ST0x0BYGRvXp3xzQt5vXl1Rolz0xthjx3A5uLqjBuOsAPYXK++ur7HsW6EHcDm2r17fY9j3Qg7gM3F1Rk3HGEHsLm4OuOG46wYAJuPqzNuKPbYAaAYwg4AxRB2ACiGsANAMYQdAIoh7ABQDGEHgGIIOwAUQ9gBoBjCDgDFEHYAKIawA0AxhB0AiiHsAFAMYQeAYgg7ABRD2AGgGMIOAMUQdgAohrADQDGEHQCKaRJ220dsv2T7Fdsfa7EmAGA8E4fd9lZJj0h6r6TbJd1n+/ZJ1wUAjKfFHvudkl5J8t0kP5f0uKQPNFgXADAGJ5lsAfuPJB1J8ueD+/dL+u0kD1y13YKkBUmanZ09sLi4ONHrbobl5WXNzMx0PcZQzNlOH2aUmLO1vsw5Pz9/NsnBYdtta/BaXuWxa/61SHJC0glJ2rdvXw4dOtTgpTfW6dOnxZzt9GHOPswoMWdrfZlzVC0OxVyUtOuK+zsl/aDBugCAMbQI+zOSbrN9q+23SbpX0pcbrAsAGMPEh2KSvGH7AUlPSNoq6bEk5yeeDAAwlhbH2JXkq5K+2mItAMBk+M1TACiGsANAMYQdAIoh7ABQDGEHgGIIOwAUQ9gBoBjCDgDFEHYAKIawA0AxhB0AiiHsAFAMYQeAYgg7ABRD2AGgGMIOAMUQdgAohrADQDGEHQCKIewAUAxhB4BiCDsAFEPYAaAYwg4AxRB2ACiGsANAMYQdAIoh7ABQDGEHgGIIOwAUQ9gBoBjCDgDFEHYAKIawA0AxhB0Aipko7Lb/2PZ525dtH2w1FABgfJPusZ+T9EFJTzeYBQDQwLZJnpzkgiTZbjMNAGBiE4V9PWwvSFoY3H3d9rnNeu0J3CzpR10PMQLmbKcPM0rM2Vpf5tw3ykZDw2773yT96ipfOp7kX0edJskJSScGa55JMvXH5JmzrT7M2YcZJeZsrU9zjrLd0LAnuXvycQAAm4XTHQGgmElPd/xD2xcl/Y6kr9h+YsSnnpjkdTcRc7bVhzn7MKPEnK2VmtNJNnoQAMAm4lAMABRD2AGgmM7DbvuvbMf2zV3Pshrbf2v7BdvP2f667V/reqar2f6k7RcHc37J9k1dz7Saab8Ehe0jtl+y/Yrtj3U9z2psP2b7tWn/PRDbu2wv2b4w+DN/sOuZrmZ7h+1v235+MOMnup7pemxvtf0d26eGbdtp2G3vkvT7kl7tco4hPpnkXUnukHRK0t90PM9qnpS0P8m7JL0s6aGO51nL1F6CwvZWSY9Ieq+k2yXdZ/v2bqda1WclHel6iBG8IenDSd4p6T2S/mIKv5+vSzqc5Dcl3SHpiO33dDvSdT0o6cIoG3a9x/73kj4iaWrfwU3y0yvu/pKmcNYkX0/yxuDuNyXt7HKetSS5kOSlrudYw52SXkny3SQ/l/S4pA90PNM1kjwt6cddzzFMkh8meXbw+c+0EqRbup3qrbJieXB3++A2dT/fkmR7p6T3S/rMKNt3FnbbRyX9T5Lnu5phVLYftv19Scc0nXvsV/pTSV/reogeukXS96+4f1FTFqK+sr1X0rslfavjUa4xOLzxnKTXJD2ZZOpmHPi0VnaCL4+y8YZeK+Z6lyOQ9NeS/mAjX39Uwy6bkOS4pOO2H5L0gKSPb+qAGu3SDraPa+W/wCc3c7YrtboERQdWu5LdVO699YntGUlfkPSXV/3vdyok+T9Jdwzel/qS7f1Jpur9C9v3SHotyVnbh0Z5zoaGfa3LEdj+DUm3Snp+cGXInZKetX1nkv/dyJlWs47LJvyLpK+og7APm9H2n0i6R9LvpcNfTujxJSguStp1xf2dkn7Q0Swl2N6ulaifTPLFrue5niQ/sX1aK+9fTFXYJd0l6ajt90naIemXbX8uyYfWekInh2KS/GeSX0myN8lerfxQ/VYXUR/G9m1X3D0q6cWuZlmL7SOSPirpaJJLXc/TU89Ius32rbbfJuleSV/ueKbe8soe26OSLiT5VNfzrMb27JtnkNl+u6S7NYU/30keSrJz0Mp7JT11vahL3b952gd/Z/uc7Re0cuho6k7bkvQPkm6U9OTgtMx/6nqg1UxwCYoNN3jz+QFJT2jljb7FJOe7nepatj8v6T8k7bN90fafdT3TGu6SdL+kw4O/k88N9jinyTskLQ1+tp/RyjH2oacS9gGXFACAYthjB4BiCDsAFEPYAaAYwg4AxRB2ACiGsANAMYQdAIr5f44gfzNwumB7AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def animate(i):\n",
    "    ln1.set_data([x1[i]], [y1[i]])\n",
    "    ln2.set_data([x2[i]], [y2[i]])\n",
    "    \n",
    "fig, ax = plt.subplots(1,1)\n",
    "ax.grid()\n",
    "ln1, = plt.plot([], [], 'ro')\n",
    "ln2, = plt.plot([], [], 'ro')\n",
    "ax.set_ylim(-1, 4)\n",
    "ax.set_xlim(-4,4)\n",
    "ani = animation.FuncAnimation(fig, animate, frames=1000, interval=50)\n",
    "ani.save('pen.gif',writer='pillow',fps=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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
}