bg01_eqs_motion_central_field.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# import pyplot and set some parameters to make plots prettier\n",
    "import matplotlib.pyplot as plt\n",
    "from codes.plot_utils import plot_pretty\n",
    "plot_pretty()\n",
    "\n",
    "# setup notebook for inline figures\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Orbits in the central potentials\n",
    "\n",
    "#### Reference for reading: Binney and Tremain \"Galactic dynamics\", S 3.1, p. 143 (+ Appendix D.3)\n",
    "\n",
    "Consider motion of particle in a *central field*, where force depends only on the distance to the origin of the field, $\\vec{r}=r\\vec{e}_r$, where $\\vec{e}_r$ is a unit vector along $\\vec{r}$. The force acting on a particle in such field has only one component - along $\\vec{r}$.\n",
    "\n",
    "Cross-product of any vectors is zero, if the angle between them in the plane formed by them is $0$ or $\\pi$. Thus, in the central field,\n",
    "\n",
    "$$\\frac{d\\vec{L}}{dt}=\\frac{d}{dt}\\left(\\vec{r}\\times\\frac{d\\vec{r}}{dt}\\right)=\\frac{d\\vec{r}}{dt}\\times \\frac{d\\vec{r}}{dt} + \\vec{r}\\times\\frac{d^2\\vec{r}}{dt^2}=r\\vec{e}_r\\times g(r)\\vec{e}_r=0.$$\n",
    "The angular momentum (per unit mass) vector $\\vec{L}=\\vec{r}\\times \\vec{v}$ is thus conserved. This vector is pointing in the direction perpendicular to the plane defined by $\\vec{r}$ and $\\vec{v}$ and its conservation means that *particle orbit in the central field is confined to a fixed plane.*\n",
    "\n",
    "This means that orbit can be described using 2 coordinates instead of 3. We will use distance to the force center, $r$, and angle, $\\psi$, formed by $\\vec{r}$ relative to a reference direction, $\\psi_0$, as such coordinates."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Equation of motion in the orbital plane\n",
    "\n",
    "<a href=\"http://www.feynmanlectures.caltech.edu/II_19.html\">Principle of least action</a> in mechanics states that particles move along trajectory, for which action $\\mathcal{S}=\\int_{t_0}^{t_1}\\mathcal{L}dt$ is stationary to first order: $\\delta\\mathcal{S}=0$. Here $\\mathcal{L}(t,\\vec{q},\\dot{\\vec{q}})$ is the Lagrangian per unit mass. The principle gives rise to the <a href=\"http://farside.ph.utexas.edu/teaching/336k/Newtonhtml/node77.html\">generalized Lagrange equations</a>:\n",
    "$$\\frac{d}{dt}\\left(\\frac{\\partial\\mathcal{L}}{\\partial\\dot{q}_j}\\right)-\\frac{\\partial\\mathcal{L}}{\\partial q_j}=0,$$\n",
    "where $j$ indicates j-th component of generalized velocity, $\\dot{\\vec{q}}$, and coordinate, $\\vec{q}$.\n",
    "\n",
    "For a particle in orbital plane in the central field: $\\mathcal{L}=\\frac{1}{2}[\\dot{r}^2+(r\\dot{\\psi})^2]-\\Phi(r)$, where $\\Phi(r)$ is the potential corresponding to the central force: $g(r)=-\\partial\\Phi/\\partial r$. So that Lagrangian equations of motions are:\n",
    "<br>\n",
    "\n",
    "$$\\frac{d}{dt}\\left(\\frac{\\partial\\mathcal{L}}{\\partial \\dot{r}}\\right)-\\frac{\\partial \\mathcal{L}}{\\partial r}=0.$$\n",
    "$$\\frac{d}{dt}\\left(\\frac{\\partial\\mathcal{L}}{\\partial \\dot{\\psi}}\\right)-\\frac{\\partial \\mathcal{L}}{\\partial \\psi}=0.$$\n",
    "\n",
    "The second equation gives: \n",
    "$$\\frac{d}{dt}\\left(r^2\\dot{\\psi}\\right)=0,$$\n",
    "\n",
    "which reflects conservation of angular momentum, $L=r^2\\dot{\\psi}=\\rm const$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Effective potential\n",
    "\n",
    "The first Lagrange equation gives\n",
    "\n",
    "$$\\ddot{r}-r\\dot{\\psi}^2+\\frac{\\partial\\Phi}{\\partial r}=0.$$\n",
    "\n",
    "Given that $L=r^2\\dot{\\psi}=\\rm const$, this can be re-written as:\n",
    "\n",
    "$$\\ddot{r}=-\\frac{\\partial}{\\partial r}\\left[\\frac{L}{2r^2}+\\Phi(r)\\right]=-\\frac{\\partial\\,\\Phi_{\\rm eff}(r)}{\\partial r}.$$\n",
    "\n",
    "Here $\\Phi_{\\rm eff}$ is the *effective potential*: $\\Phi_{\\rm eff}(r)\\equiv L/(2r^2)+\\Phi(r)$. Physically, it accounts for the centrifugal force due to angular momentum of the particle motion, which counteracts the central force. Here is how it looks like for a star orbiting around the center of the Milky Way. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqUAAAG5CAYAAAC3JE1PAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAa6QAAGukBho3xbQAAIABJREFUeJzs3XlYW9edP/73kZDEjlhtvBvvxo5jwHGcpUlsyJ62SU08\nXWa62rTTdqbfLibuTGemM/MbF9ppp53OpJB070xrQ5JJ3aRJkeM6m+MY5HjfZccLZkfsSEg6vz8k\nBEICBAjulfR+PY8edI/Ovfpgg+6HswopJYiIiIiIlKRROgAiIiIiIialRERERKQ4JqVEREREpDgm\npURERESkOCalRERERKQ4JqVEREREpDgmpURERESkOCalRERERKQ4JqVEREREpDgmpURERESkOCal\nRERERKQ4JqVEREREpDgmpURERESkOCalRERERKQ4JqVEREREpLgYpQMIN0KIrVLK6nHq5AAoA9Dm\nKUoDsFtKaWZ9IiIiIn9CSql0DGFBCJEHYBeArQBSpZTWMertB7BlMCnzJG11AIqllKZork9EREQU\nCLvvxyGEqBJC1AEoBJATxClVACqHtxJKKS0Adntei/b6RERERH6YlI5DSlkspcyXUpYDsIxVVwgx\nmLjuCfByNQCjEGJrtNYnIiIiGg2T0tAq8nz1S149rYfD60RjfSIiIqKAmJSGVjCtgk9GcX0iIiKi\ngDjRaQKEEFUYY6KTEEICgJRSjHJ+OwDj4OvRVp+IiIhoNGwpJSIiIiLFMSlVgBDCyPpEREREQ7h4\nPqna4BCBINinNRAiIqLIpAdGH4o3k5iUKmC0hfdZf/L0er0eAFwul99rGo1/h4CUEiPHUwshIIT/\n72SguoGuOZPvP1rdQO/vdA6dq9VqoNGIGX3/Kf9bDUhg2NvLGDmz76/wv5XS7z9dsUb7+/NnJXze\nP9zvAeO9v8Ph8HtdKUxKQ8sKYKyu6pGvRVv9ybDr9Xq9zWYLwaWij9Ppwkce/BFsNveHzrq8BfjO\nf4TXggj1z9Tj/I7z3uOl31uKeX87T8GIiIgih8FggN1uV0VvI8eUhpYFGHcM5fD94KOtPs0wrVaD\nRTkZ3uPLl5r9/uJWu4wPZgDDGgSaX2hWLhgiIpo2TEpDa3Cfd7/tSD37wQ+vE431SQGLl2Z5n3d2\n9KGluVvBaCZOP0uPlDtTvMcdb3TA3qyKP+opgM7OTnR2diodBhGFISaloVXh+VoQ4LU8z9fhW3JG\nW/1JCTRGhoK3ZGmmz/HlS00KRTJ5GY9nQMQIpBalYtmPl0Fj4EeXWh08eBAHDx5UOgwiCpKa7rH8\nZJ+YwW7qtEAverbWrARQEuDlXQCqpZTmaK1PyhjeUgoAly6EX/d39mezcUfTHVj3p3WY+4W5iEnm\ncHgiokjDpHQcQoidQogaz+5EhZ7iS0KIOk+5z3hKKWUJAIsQokoIkeN5VAFok1IWj7x+tNWnmbd4\nSfi3lMakxECXqlM6DCIimkbcZnQKhBDG0ZY7EkLkYSiJNY3XYhht9YMlhLDFxMToBwYGQnG5qPWZ\njz6LmzfcP6pz56Xi2f/9rMIRUaTat28fAOCxxx5TOBIiCoZOp4PD4bBLKQ1Kx8I+sCkYa/1NT1IW\ndGIWbfVpZuUsyfQmpfU32tHfZ0dsnF7hqIiIiIaw+55Ub7RFiil4OcuGxpVKCVy2tCgYTWj0Wfrg\n6FbPos/kds899+Cee+5ROgwiCpKa7rHqiYSIpk3OiHGllovhN64UAOzNdlwtv4raglocXnIYLS+E\nf3IdaZKTk5GcnKx0GEQUhpiUEkWB4S2lAGC5GH4z8AHA0e6ApdSC7jr3WqvNe8Pz+yCiiSkpKUFR\nUZF3e83i4mKUl5crHRaFGMeUEkWBzKwkJCYa0N3t3q41XFtK45fHI2FdAnqO9QAA2l5tw4B1ADoj\nZ+YTRbKKCvey2KmpqSgsLERVVZXCEdF0YEspURQQQvi0ll6xtMDpVM+CyROR9eSw8bEDEq0vtioY\nDRHNBLPZjPLyclitVlgsFlRXV8NisSgdFoUYk1JSPS5bFhrD1yvt7xvAzfpRF49Qtcxi3/GxTXvC\ns9U3Up0/fx7nz59XOgyKECaTCfn5+aitrcXOnTthNBpRUFCAwsJCVFRUID8/n8npFKnpHsuklFRP\nTb8w4Sxn5HajYTquNH5ZPBLXJ3qP2/7UBnuzXcGIaLhz587h3LlzSodBEcBsNqOoqAi7du3Cjh07\nfF4zGo0oKytDQUEB8vPzYbWG5x/ZaqCmeyyTUqIokTNyu9EwHVcKAFkfHfa9OIHmqvBMsIlodLt3\n7wYAbN26ddQ6JSUlsFqtqKysnKmwaBoxKSWKEgsWpUOrHfqVD9eWUsCTlIqh48b/aVQuGCJSTFtb\nGwB3y2mwTCYTSktLUVxc7E1mA5XRzOPse6IoodfHYP7CNFzxLJwfrstCAUDsvFikfCAFHQc7AACd\nb3ei73If4hbHKRwZEYXKrl27UF1djerq6lFbS6uqqmA0GvHkk08Gfd3i4mJcvnzZZyxqoDKaeUxK\nSfWEEONXoqDkLM3yJqUtzV3o7OhDckp4JnKzPj7Lm5QCQNPvmrBw10IFIyIAiI+PVzqEqPcvf/ci\n6q+3Kx2GjznzUvGt/+9DEzonLy8PNTU1KC0tRVtbm8+4UqvVit27d8NisaCurm5CLaVWqxVGoxF5\neXljlkULNd1jmZSS6qnpFybcLfbb2akZt+YvUCiaqcn8SCYufPEC5IB7kH7j/zRiwVML+POisC1b\ntigdQtSrv96OK5cjY7ezwsJC1NXVwWQyeZeEqq2txd69e1FSUoKcnBy/cwYT1vT0dABAa2srysrK\nYLVaUVpaCsA9FnXJkiXYsWOHX9nOnTtn7htUATV9Zgo1zboiGkkIYdPr9XqbzaZ0KBHhaO37+OZX\nhxad3v7Fe/HEtgIFI5qaEx8+4bNOacGxAiTekjjGGUSR7wuf/IXqktJFizPw9C8/NeXrpKam4skn\nn/Quph9IUVERSktLUVhYCAAoLS31JqCAOwkbmfsEKosWBoMBdrvdLqU0KB0LW0qJosjipSNbSsN3\nBj7g7sIfTEpjl8TC3sSloYjmzEtVOgQ/U42ppKQEFovFO9O+ra0NGzZs8GvVtFgsMJlMnLAUpthS\nSqrGltLQ+8QTP0Fri3vv+MVLMvHfP/+kwhFNnrPPiUtfu4Ssj2Uh5c4UVXVDEdHMs1qtSE1NRXt7\n+6jjTNlS6ktNLaVcEoooygxfRP/qlVbYbQ4Fo5kabZwWy/97OYx3GZmQEhGMRiMKCwt9WkktFgtb\nTcMEk1JSvWj963W6LF0xy/vc6XSF9SL6pD6vv/46Xn/9daXDoChWU1OD1tZWlJSUoLy8HNXV1d7x\npCUlJQDcS0ANLv8UqCyaqOkey+57UjUhhC0mJkY/MDCgdCgR4523LuHbu17wHn/+bzfjQx+JvmVQ\naHrs27cPAPDYY48pHAkRBUOn08HhcLD7nohm3vJhLaUAcOFsg0KRTB/p4h/bREThhkkpUZRJy0hE\neubQsknnIyQpdTlcaH2pFaeePIVjW44pHQ4REU0Ql4QiikLLV87GoeaLAIDrV9vQ22tHfLxe4aim\n5vRfnEbLc0NrM/ac7UHCygQFIyIioolgSympHmdVh97ylbO9z6UELp4L/9bSjA9m+BzffPamQpFE\nt5ycnIC77BCROqnpHsuklFRPTb8wkWJ4UgpERhd+5kcyoU3Weo8bf9kIl92lYETRKTc3F7m5uUqH\nQURBUtM9lkkpURRaNiIpvXC2UaFIQkeboMWsjw1N4hpoGUDLi+raapGIiEbHpJQoCiUlxSJ77tBu\nJ5HQUgoA2duzfY7ZhU9EFD6YlBJFqeFd+A03O9Bh7VUwmtBIyktC4vqhlQXaa9rRd6VPwYiIiChY\nTEqJotTIcaUXzoV/Fz4worVUAg0/i4xW4HDR39+P/v5+pcMgojDEpJRUz+XiZJXpEImTnQAg66NZ\n0MQNfbTdfOYmJzzNoJqaGtTU1CgdBhEFSU33WK5TShSllizLgkYj4PLsfnT+TGQkpTqjDlnbstDw\nC/f3Y2+wo+WFFmRty1I4MiKablarFbt374bVaoXFYoHRaMQzzzwDo9E4/skKC+fYQ4UtpURRKi5e\nj/kL07zHFyJgrdJBc7801+f4xo9vKBQJEc2k7du3o6SkBBUVFaipqUFaWhq2bNmidFhBCefYQ4VJ\nKVEUG96F39bag5bmLgWjCZ2k/CQk357sPe54swNd70XG90ZEozObzTCZTN7jkpISmM3mkF2/srIS\nFoslZNcbLpjYy8vLQ/r9qA2TUqIoFqnjSgHf1tKY1Bj0XeQsfKJIV1dXhx07dniPa2trJ9T9bTKZ\nUF5eHvC18vJypKWlTduOZcHEvnPnTpSWlk5bYqw0JqWkehoNf0yny8hF9CMpKc3cmonUolSseHYF\nNl3fhKytHFM6EzZt2oRNmzYpHQZFmOLiYqSmpkIIgSVLlqCoqChgi+HIJK6srAzPPPNM0O9TUVGB\nwsJCv3Kz2YwjR45g69atfq9VV1ejuLgYxcXFyM/PR2lpKaxWa9DvOdHYKyoqUFxcPOHrj0ZN91gh\npVQ6BqJRCSFser1eb7PZlA4lItntDnzkwR/B4XDPvlxfsBD/9v3QfdgREYVKcXExqqurEWzeUlRU\nhJKSkoCJZCBWqxX5+fm4dOmS32v5+fmoqqryayWtrKxEWlqa9z2sVqt3HGhdXV1Q7zuZ2EtKSlBc\nXBwwgZ4og8EAu91ul1IapnyxKVJPekxEM06vj8HipZne4wvnGoP+wCcimklmsznoJKy4uBhlZWXY\nunWrzzjNsezduxclJSUB3xdAwG77iooKHDlyxHtsNBpRVlYGs9mM6urqoN53MrEXFxejtLR0UtdX\nMyalRFFu2YqhLvzurn7crO9QMBoiIn+DyyQVFRWNW7eoqAhFRUVoa2tDdXU1ysrKgnqPiooKnzGd\ng/bs2YNt27YFPCctLc2vq34weR2erAYr2NgLCwthNpsnNUxAzbhOKVGUW75yNl5+8Zj3+OypesyZ\nG7nr4rnsLmj0/HucKJwMthiO11Kan5/vN4s9Ly9v3OubzWakpaUFnBRlMpmwa9eugOcF2ihicBLS\nhg0bfK5RVVUFi8WCqqoqVFZWorW1FVar1dsNP9HY8/LyYDKZgh6eEA6YlBJFuVW5c3yOz5ysx+b7\nVysUzfTpu9SHG0/fQMPPG5D3Th7il8UrHVJEunLlCgBg0aJFisZBwJG1E2+py3giA4u/vdiv/Gr5\nVTT+euJbEa8/tB4xiVNPNQaTv/ESzMmO46yoqAjYdQ+4k8yJzLivqKiA0Wj0SRZrampQUVGB/Px8\n5Ofno6amBjk5OViyZAkAd7I90djT0tIibhY+k1JSPTVtgRaJ5i1IQ2JSLLq73PuVnz5Vr3BEodf6\ncitOPHoC8AyXvfFfN7DsP5YpG1SEOnHiBAAmpWrQc7JnwuckbUwKWG5vsE/qegjRx7fJZAqqxXMq\n16+oqAj4mtVqDXpZqerqaphMJp8E02q1epNPs9mMsrIyb5JbUVGBgoKCScWck5OD1tbWSZ07nJru\nsezDIopyGo3A6jVDraVXLjWjt9euYEShZ7zHCG2y1nvc8NMGDLQPKBgREQVrcDzpaF33Ux1bWV1d\nHZJZ7GazGaWlpairq/NpWTUajdixY4e3VXN4C2phYWFUbSM6HraUEhFWrZmDdw+5PzBdLolzp29i\nfcFChaMKHW2CFnO2z8G1710DADi7nah/uh4Lvxk53yPRSAlrEiZ8jmFu4FWB9LP1k7peKJq+BsdY\njjbJaffu3aiqqpr09SsqKsadDDVe0muxWLB9+3bU1dV5k8yR3f4mkwlGozFki++3tbUhPz8/JNdS\nCyalRITVa333ij994kZEJaUAMO8r83D9h9chB9x9+Nd/eB3z/t88aOO045xJFJ42nNgwfqUgLdi5\nAAt2LgjZ9SZicDxpoNZMk8k0pSTParWira1t3AlFbW1to75usVhQWlqK/fv3exNSs9mM2tpan9n8\nNTU1IWmRHf6+07W7lFLYfU9EWL5yNrTaoY+D0ycjb1ypYa4Bsz4xy3s80DSAhl9Gzg5WahETE4OY\nGLZ3UOiMNp50cCel0SYoBaOysnLU5Z4GFRYWBpxlD8C7TFVaWhp2796N0tJSlJSUYPv27X4Jo8lk\n8pmRP1Vms3nS41HVip8cpHpq2gItUsXG6rBkWZZ3m9Ezp+rhdLp8EtVIMP8b89Hw86FE9Nr3rmHO\n9jkQWqFgVJHloYceUjoEihDl5eXYs2cPLBYLjEajd2tNq9WK2tpaWK1W5OXlTam1cM+ePdi/f/+Y\ndYqKikZdqL64uBgWiwWVlZV+rwVKGEO1fNPgRgKhGI+qpnssk1IiAuAeVzqYlPb12nH1SisWL8kc\n56zwkrAqAekfSkfri+4Zq/2X+tH8fDOyirMUjoyIRtq5cyd27tw5bdcfa23S4Qa73AN1l09kGaf2\n9vaJBzmKsZawCmfqSY+JSFG5I8aVnjpxQ6FIpteCUt9xcVfLrnJrVaIoNJHE7plnnlHNtp4WiwW1\ntbURtWj+ICalRATAf7LTmQgcVwoAKZtSkHJ3ive4u64b1tcia6s+IhrfRHZDGhwmMNn97EOppKRk\nSqsNqBmTUiICAKRnJCJrdrL3+HSEtpQCgVtLiSh6TGZt0rKyMhw5ckTRXZQqKyt9Ft+PNExKSfXY\ntTpzVq8Zai1tuNmBtpZuBaOZPmkPp/msudhe044uc5eCEUWOw4cP4/Dhw0qHQTSmmpqaSXXHK50Q\n7tixI+Q7W6npHsuklFRPTb8wkW74zk5AZC4NBQBCCMz/xnz3gQaY9ZezoEvXKRtUhGhqakJTU5PS\nYRCNqaKiImJbGydKTfdYzr4nIi//caU3cNe9yxWKZnplfTQLne92Yu6X5iJh5SR2qiEiopBiUkpE\nXosWZyAuToe+Pve+8KcitKUUADQ6DZb/ODITbiKicMTueyLy0sZosDJ3qAv/0vlG2GwDCkZERETR\ngkkpqZ4Q3G1nJq0aNq7U4XDhwtlGBaOhcDN37lzMnTt3/IpEpApquscyKSXVU9MvTDQYPgMfiNxF\n9AORTonG3zaiq44z8ScrLy8v5LODiWj6qOkey6SUiHyszM2GZthe8MePXlMwmpkhXRJNVU04cssR\nnPnYGVieUm4dQiKiaMWklIh8JCQYsGz5LO/xqRPX4XA4FYxoBkjg8rcuo/d0LwCg3dQO6+vc5YmI\naCYxKSUiP2tvne99but34PyZBgWjmX5CK7Donxb5lF3+h8uqWr+PiCKLyWRCeXm50mGoCpNSIvKz\nLs93G85jUdCFn1WchfjceO9xx8EOtL3apmBE4cnlcsHlcikdBpGfyspKFBUVoaioCMXFxSguLkZl\nZSUAoLS0NOD2oSUlJViyZAmEEOPuAGUymSCEQGpqKoqKisbdjrSiomLCW51GOialpHq8wc283LVz\nR4wrjfy94YVWYPG3F/uUXfr6Jbgc/PmbiJdeegkvvfSS0mEQeVVXVyM1NRWXLl1CVVUVampqUFVV\nhaqqKuTk5CA/Px+VlZUBd3iqqKjA1q1bUVhYCLPZPOp7WK1WVFVVAYD3PcbaMcpqtcJsNqtiUqCa\n7rFMSonIT1y8HstXzvYenzlZD7vdoWBEMyPjiQwk357sPe491YuGX0T20AWiSFZaWort27dj//79\nKCsrg9Fo9Hm9sLAQBQUFePLJJ0e9hsViQVFREWpra0etU1lZifz8fO81x7N3716UlJQE+V1EDyal\nRBTQuvVDXfg2W+SPKwXcS6Ms+fclPmVXvnUFju7IT8iJIk1lZSXKy8tRVVU1Zotkfn4+ioqKRn09\nLS0NeXl5sFqtsFr9J0CaTCYUFhairq4u6JbPiooK7NixI6i60YTbjBJRQLesn489vznsPT529CrW\nrJunYEQzI+WOFGR8JAMtz7UAAOwNdlz73jUs/qfF45xJpC5r1078nCeeAL79bf/y8nLg17+e+PUO\nHQISEyd+3lRZrVaUlJSgsLBw3JbLsZJDk8mE/Px8FBQUAHC3mo5MPM1mM3bu3AmTyYStW7eOG5vZ\nbEZaWppfq63JZEJVVRUsFguqqqpQWVmJ1tZWWK1WFBcXR8X4UyalRBTQ6jVzEBOjgcMzpvL40Wv4\n+KeUjWmm5HwnB60vtkI63LPvr333GubsmAPDHIPCkREF7+TJiZ+zcWPg8oaGyV1PqeGKgxOYptpF\nXlNTg5KSEhiNRhiNRtTW1vokpeXl5di5cyesVissFgs2bNgw7jUrKioCxlVTU4OKigrk5+cjPz/f\nOy51yRJ37w2TUpoUIUQZgDwAZQBqpZRWIUQegEIAG6SUxQHOyfHUH5zumwZgt5Qy4MjqcK8/ERoN\nR5koITZOj+WrsnHas6PTmVP1sNsc0Bsi/2Mjfmk85nxxDm780P29u3pduPwPl7Hy2ZUKR6Z+g+Pq\niJRUU1MDAAFbLs1mM3bv3g2LxQKr1Qqj0Yht27Zh586dfnUtFot3wlJOTg7q6uq8r5lMJm+CajKZ\nAASXOJpMJlRUVPiUWa1Wb/JpNptRVlbmfd+KigpvS+10UNM9NvLvLsowwp2AFgI+W3hZAfh9YnsS\n1v0AtgwmcZ4kr04IUSylNEVSfQof69bP9yalA3Ynzp6uxy3rF4xzVmRY9K1FaPhFA5wd7o0DGn7W\ngHl/Mw+JtyjQFxlG5syZo3QI5LFmzcTPmTs3cPns2ZO73lTznfz8/IDjOAOpq6vzdom3tbX5dY8P\nysvL83aPl5SUoKamJqhksqCgwGeZp8FuewDeVs3R3nNQdXV1wPcyGo3YsWOH9/rDk+loaCEdxKR0\n+ljgbi00ep5XSylHW+SsCkDl8FZFKaVFCLHb81pqhNWnMHHL+gX47a/e8R4fO3otapJSXboOC/9+\nISzfcN8kEtclQrq4mD6FjxMnQnetnTvdj5k2vGVyIgoKCsZcwmn4tUdL+sxms88EqPz8fOzduxfA\nULf9oMHJTuOpqKhAWVnZqK+bTCYYjcYxl5OKZOpps408xVLKVCmlkFIuGS0hFUIUAsgBsCfAy9UA\njEKIrZFSn8LLqjXZiNFpvcfHzZG/iP5wc780F4l5iVj242XIr81H0q1JSodEREEYHLNZXV09ah2T\nyTRm8jcy0czJyYHVakV1dbXPuNLB8aRjzeAfrNfW1jbmDP1gW20jFZNS5Q3+FPtt/SCltIyoEwn1\nKYwYDDqsXJ3tPT57uh79/QMKRjSztLFa5NfmY+4X50IM20yAiNQtLy8PZWVl2L59u3e853CDk5fG\nShD37Nnjk7QOjuscmTgGO560srIS27ZtG7OOyWQKarJUpGJSqrxgWhGHr+ob7vUpzKxbP9/73OFw\n4cypegWjmXnDxoRTEG7evImbN28qHQYRdu7cif3796OiosJna9Hi4mKUlJSgqqoqYOtmeXk5lixZ\nArPZ7P0KwJvEDna/D7aObt++HQBQXFw8Zsvsnj17glqbNJhlpSKVkJJjpEJNCFEBoALuiU7pnuIc\nBJiNLoSQACClDHjnE0K0AzAOvh7u9SdKCGGLiYnRDwxET+uc2hx/7xpK/2ZodMa2T2zEp3bcrWBE\npGb79u0DADz22GMKR0KkHmazGaWlpd5VAdREp9PB4XDYpZSKr3nHltLp8wzck39KPeNJd8M9Gz16\n/wSisLRydTYMw5aBqnv3inLBqESXuQvnv3Qe/KOeiIIx2tqk5Iuz76dHDYAyKaV3HQsppVkIUQ13\nsjp6+/4ohBDG4deL9PrDORwO6HQ6n7KR66rZbDY0NDTgyJEjfuc/9NBDiInx/VE/duwYrl696lOW\nkZGBTZs2+Z3/6quvwm63+5Tl5ub6DZBva2vDW2+95Xf+li1bEB8f71N25swZXLx40acsKSkJ9957\nr9/5Bw4cQHd3t0/ZsmXLsHKl75qZ3d3dOHDggN/5d911F1JTfRdAuHTpEk6fPu1TZjAYcP/99/ud\nX1v7LmbNNeCqxb3V5sXzjXjrzSO48y7fcU8DAwN45ZVX/M6/7bbbMGvWLJ+y69ev4+jRoz5lQgg8\n+uijAd6/1q87ePbs2QHHXf3hD3/wSxTXrVuHBQt8VwxoamrC4cOHMdIDDzwAvV7vU3bixAlcuXLF\nfdANGH5jgP6PekACyRuSMfuTs711a2pq0N/f73P+qlWrsHTpUp8yq9WKN954w+/977vvPiSO2P7m\n3LlzOH/+vE9ZQkICNm/e7Hf+wYMH0dnZ6VOWk5OD3Nxcn7K+vr6A4+zuuOMOpKen+5RdvnwZJ0es\nmq7T6fDggw/6nf/OO+/A5VktfbDFdN68eVi/fr1PPafTiZdfftnv/IKCAmRnZ/uU1dfXB5yB/cgj\nj/h9DpjNZty4ccOnLCsrCxsDrAj/xz/+EQ6H7/axa9euxaJFi3zKWlpacOjQIb/zi4qKEBsb61N2\n6tQpnyWDACAlJQUf+MAH/M7fv38/ent7fcpWrFiB5cuX+5R1dnbi4MGDfuffc889SE5O9im7cOEC\nzp4961MWFxcXcKzjm2++ifb2dp+yRYsWYe2IbaBsNhv+9Kc/+Z1/++23IzMz06fs/fffx/Hjx33K\ntFotHn74Yb/z3333XTQ2NvqUzZkzJ+A6t4M/S8Pl5eVh7oj1q9R8Dxi5Nqka7wFqwJbSaSClrB42\nyWe4I3DPRldgYQ2iyVuwKMHn+MLZVoUiUZb2rBb6l90JKQBc/OpF2BvtY59ERFHthRdeiOoZ9RPB\nMaUzyNN1XwX3mqXFnjJVjfnkmFIK5Pq1Nmz/+M+8x5vvX41v/L1/60c0OPn4SbT8X4v3OLM4E7l7\nc8c4I7pwTCmRr5KSEpSWlqp27VE1jSllUhpiQggjgIJAuxh51vSsAWCRUi7xlI2ZtI1M+sK9/kQJ\nIWx6vV5vs9kmczqFiJQSn972DBob3F3DKcY4/O///TU0muibmW6rt+Hd1e96d3oCgNzncpH5ROYY\nZ0UPp9P976LVasepSURqYDAYYLfbVZGUsvs+9KoA1IzSRT+4/9jwrn0L4E1mRzN8xn6416cwJIRA\n/sbF3uONXhMTAAAgAElEQVQOax8sF5oUjEg5hjkGLP133zGi5//6PAba2JoPuJNRJqRENBlMSqeH\nFYETscG2++GvmUa85uXZP354nUioT2Eq/7ZFPse1715WJhAVmP2Z2UgtHJo4MNA4gItfvTjGGURE\nNB4mpaFXA2BLoO57DO1sVDGsbPB5QYD6g1tNDN/CM9zrU5hal7cAWu3QR0Y0Lw0lhMDyyuXQxA/9\nezT+shGtr0TnBDAiolBgUhpiUspyAGUjyz3jSQsBlA6fme95Xgkg0AJmu+CeFGWOlPoUvhISDFi1\nZo73+MzJevT0RO9Y37jFccjZ7dtBcO6z5zDQym58IqLJYFI6PUqEEDVCiB1CiELP+NIqACWepNWH\nlLIEgEUIUSWEyPE8qgC0Dc7Sj6T6E8XJeOpRMGxcqdPpwjHz1TFqR765X5qL5DuH1oq019txruRc\nVP/MHj161G8NWiJSLzV9XnH2/TTytI7mwT2GtHa8xeOFEHlwt6YCgGm8FsZwrx8MLgmlLhfPNeLL\n23/tPX74Q+vw5a/57x0dTXov9qL21lq4elzeshU/X4HsT2WPcVbk4pJQROFFTUtCcUenaeQZVxr0\nJB9PEhd0Ihfu9Sn85CzLgjE1HtZ29040dYcvQ0oJIaJvaahB8UvjsexHy3Dus+e8Zc17mzH7k7Oj\n+t+FiGii2H1PREHTaATyNizyHjc2dOLG9fbRT4gSsz89GxmPZ0DECCz+t8VYu28tE1IiogliUkpE\nEzJyaagjhwLtqBtdhBBY8cwKrD+0Hgt3LYTQMiElIpooJqWkemxxUpf82xb57OT0zluXFIxGPXTp\nOiQXJI9fMcJlZmYiM5O7WxGFCzXdYzmmlFRPTb8wBKQY47FqzRycOn4DAHDy+HV0dfYhKTlO4chI\nDW6//XalQyCiCVDTPZYtpUQ0YbffObTNpsspceSd6N3daTztB9pxcutJuAZc41cmIopi05qUCiE2\nCyF2CyFeFUIcEUK0CiGcnq8XPGV7hBCfFUKw34soTGy6y3fv90NvcovNkaRL4sq/XsGxwmNoea4F\nl/+eiTsR0VhCnpQKIZKFEN8RQrTBvRxSMYAOAHUAngHweQDfAfCcpywVwHcBtHuS1/tCHRMRhdbc\n+amYvyDNe1x3+DLsdoeCEalP57uduPKtK4CngfRa+TU0P9esaExERGoW0jGlQohvwL21ZA3cyWit\nlLIjyHMXe855TghxBMAOKeX7oYyPiELn9ruW4tr/vgsA6OsbwPGj13x2fIp2KbenYME3F+Dqvw3t\nenXmr84gblkcEm9JVDAyIiJ1CklLqad19E8AcgAsllJuk1LuDzYhBQAp5WUpZbmUMg3uFtX9QojH\nQxEfhTeXi2Px1Oj2O5f4HLML39+iby+CcbPRe+zqdeHkh07C3mJXLqhp9sorr+CVV15ROgwiCpKa\n7rGh6r5/FkCZlPILE0lERyOlrAZQAOALQojNU46OiEJuxepsGFPjvceH37oEl4vbFg+nidEgd28u\nYhfHesv6r/TjdPHpiJ34NDAwAG4LTESTMeWkVAixHsB2KeX+EMTjJaW0Sinvh3vMKRGpjFarwcY7\nhlpLW1u6cfF8o4IRqZMuXYc1L66BJmHo49b6Zysu/u1FSMkknoho0JSTUinl0VC0jo5x/eem69pE\nNDUju/DfYRd+QIlrE7Hq16t8yuqfrsf1719XKCIiIvXhOqVENGm3FiyEwTA0X5LjSkeX+XgmFv3T\nIp+yS1+/hKaqJmUCIiJSmVB03z/hWYv08RHli4UQn/OsVco1SGnSNBr+7aRWsbE65G1Y5D2+YmlB\n/Q2rcgGp3MJ/WIisj2f5lJ35yzPoeGvaOptm3Jo1a7BmzRqlwyCiIKnpHjvpSIQQKZ6lm6oBlMK9\nlFOLEGIh4J5ND/c6pOUA2kIRLBGpz6a7fRfSf+PAOYUiUT8hBFb+dCWM9w7NyJc2CctTlogZX7p4\n8WIsXsylwYho4qaSHu8CIAAUAVgC4H4ABwCYhyWmR+FOTNWzsSoRhdTtdy1FTMzQR8nrTErHpDFo\nkPt8LuJXuVcuSLknBWt+v0ZV+08TESlhKklpIYDNnvVIL0spTVLKYgAlAExCiCRPvfYpR0lEqpWU\nFOvThW+50ITr19g5MhZdqg63/PEWzP3SXKx7dR10qTqlQyIiUtxUklKLlLJzZKFnjdFtAF7jWFKi\n6HD3fSt8jt94ja2l44ldGItl/7kMGoN6xnMRESlpWj4NpZRmuBPTZwEYx6lORGFu011LEaPTeo/Z\nhT81/df7MWANzwXoW1tb0draqnQYRBSGppKU7hk28/7pkS9KKS0AnoK7m59o0tS0BRoFlpBoQMFt\ni7zHVywtuHqFiclk9F/vx3v3vofjRccx0B5+ienbb7+Nt99+W+kwiChIarrHTjop9Sxqnw73RKec\nUepY4J4AdXSy70NE4cGvC5+tpRNmu2HDsfuOof9SP7pqu3Cs8BgG2sIvMSUimowpdd9LKZ+RUhZI\nKR8Yo45FSlkwlfchIvW7/a6l0OnZhT8VXUe70He5z3vcbe7GsS3HMNDKxJSIIh9H2BNRSMTH67Fh\n41CnydUrrbhiaVYwovCT8WgGcvfmQsQMLQ/V/V433tvyHuzNdgUjIyKaftOelAohkgcfI8o/J4TY\n4xmTumi64yCi6Xf3Zt8u/Nc5C3/CMp/IxOqq1T6Jac+xHhy96yj6rvSNcSYRUXibiZbS1+Beq/QZ\nIcRmABBC7AVQASAV7oX1q4QQt85ALBSG1LQFGo1t46YcGAwx3uMDpjMRs1PRTMr8cCZyn8uF0A0l\npn3n+3D0jqPoPt6tYGTjKywsRGEh57cShQs13WNnIpI9AJZKKbdJKV8TQmwBsBVAtZTyfinlU1LK\nDXAvuk9EYSwuXo+Ndy7xHjfUd+DUiRsKRhS+Mj6YgTX/twaauKGPaftNO47efRTWg1YFIxtbXFwc\n4uLilA6DiMLQTCSlaVLKy8OOiwFIALtH1OuYgViIaJoVPpjrc7z/lVMKRRL+0h9OxzrTOsSkDrU+\nOzudOPbAMTQ/x/G6RBRZZiIpHbnfYCEASCnfG1HOPj6iCJBXsAipaQne49cPnIPNxtnjk5VyRwrW\nv7kehnkGb5m0SZwqPoUbT7MVmogix0wkpd5kUwixGO41TU0B6qXOQCxENM20MRrcV7TKe9zbY8eh\nNy8pGFH4S1idgPWH1iM+N95bponVIHkjd3ImosgxE0lphhBined5GdxJasXwCp4dofbOQCwUhjhR\nJvywCz/0YufFYv0b65FydwoAYOUvVyIpL0nhqPydOnUKp07x/5soXKjpHjvtSamU8ikA3xVCuOCe\n4FQppXxeCJEihPiJEKIN7klOnOhEAanpF4aCs3hJJnKWZXmPzUeuoK1F3bPGw4EuVYd1Neuw5sU1\nyCrOGv8EBVgsFlgsFqXDIKIgqekeOyPrAEgp74e7ez5VSvkFT1kH3C2mWwDkA3hqJmIhoplR+MBQ\na6nLJfFazRkFo4kcGoMGGR/MCPias8+Jhl81qOomQ0QUrBlbnEpK2eFJRIeXHR32uDzauUQUfu4r\nWgmNdmidzf2vnGKyNM0ufPkCzn7yLE4/eRqODofS4RARTciUk1JPF/znRu7YRETRzZiagIKNi73H\nVy634OL5RgUjimwNv2xAw08bAADN1c2oK6hD19EuhaMiIgpeKFpKKwAUALgihHhVCPFECK5J5CWE\nGL8SqVLRg2t8jl/+/XGFIol8HYd8l3ruu9gH80YzrpZfhXTOXAt1cnIykpPZRkEULtR0jxWh7E4T\nQnwEwDa4x4maAOyRUj4fsjegqCOEsOn1er3NZlM6FJoEh8OJv9pagfa2XgBAbJwOv3nu80hINIxz\nJk2UlBI3f3oTF798Ea5+l89rKR9IwapfrULswliFoiMitTIYDLDb7XYppeIfzCEdUyqlfE5K+aSU\nMh3uJZ4+L4RoFUI8PbjvPRFFj5gYLe5/ZK33uL9vAAdqTisYUeQSQmDO5+Yg7508xC333eaz4/UO\nHLnlCG7+4ibH9RKRak3bRCdPgno/3IvlmwGUexLU3UKIW6frfYlIXR569BYM7x166cVjTIymUeK6\nRBSYCzDn83N8yp2dTpz79Dkcv/84+ix9CkVHRDS6mVintENK+YyUsgDusadtAKqFEBeEEF8XQiya\n7hiISDmzslOQf9uwCU+WFpw9dVPBiCKfNkGL5U8vx9o/rIUuS+fzWrupHUfWHMHV712Fy+Ea5QpE\nRDNvxpaEAgAp5WUp5XellEsB3A8gA4BZCHGEM/iJItcjH1rnc/zy748pFEl0SX8kHRtObEDG477r\nmrr6XLj6b1fhaOOyUUSkHjOalA7nSVCfklKmAdgBzuCnUbCrN/xtuD0H6ZmJ3uPXXzuHrq5+BSOK\nHvosPdY8vwa5z+VCn633li/9wVLos/RjnDk5r732Gl577bWQX5eIpoea7rGKJaXDeRbP/7wnQS0H\n8Bee8ad7OEGK1PQLQ5OjjdHgoUdv8R7b7Q6Y/nhSwYiiT+YTmdhwegOyd2QjtSgVs/5q1rS8T09P\nD3p6eqbl2kQUemq6x6oiKR1OSrl/2Az+KgDFSsdERFP3wKNrfXZ4+v1zR+F0ckzjTNIZdVhRsQJr\nX1obcG1C2w0b6jbWofWVVlXdqIgoOqguKR1OSlktpfyC0nEQ0dRlZCbhjruXeY8bbnbgnTcvKhhR\n9NLoAn/0W3ZZ0PVuF048dALv3fee34L8RETTSdVJKRFFlseL832OX6iqUygSGqnjnQ40/npoG9iO\ngx04esdRnPjgCXSf6FYwMiKKFkxKSfXUtAUaTc2qNXOwYtVs7/Gp4zdw/myDghHRIG2cFkkbk/zK\nW/e1onZdLU49eQpd73WNe53ly5dj+fLl0xEiEU0DNd1jp5yUCiGe8CyI//iI8sWeZZ42c6knmgo1\n/cLQ1Agh8OEnC3zKXtjL1lI1SFyXiLxDech9Lhfxq+J9X5RAc1Uz6tbX4fgjx9Hx9ujd+itWrMCK\nFSumOVoiChU13WPFZAezCyFS4N7ffnh/XBuAfCnl+5466wE8A+BWKWXMFGOlKCSEsOn1er3NZlM6\nFAoRh8OJT297Fi3N7lY3rVaDn+/Zjsws/1Y6UoZ0SjT8ugFX/vEKbFcD/+6l3JOChbsWIvX+VFXd\n1IhoYgwGA+x2u11KaVA6lqm0lO4CIAAUAVgC92L4B+BeDH8h4F7qCUCdpx4REWJitPjQ1vXeY6fT\nhX3PH1UwIhpJaAWyP5WNjec3Yul/LoVhvv+9quNgB049eQrObqcCERJRJJpKUloIYLNnCafLUkqT\nlLIYQAkAkxBisNmjfcpRElFEefDRWxAbN7T95UsvvoduLqavOhqDBvO+NA8bL27Eip+tQNzyOJ/X\nsz+TjZgkdoIRUWhMJSm1SCk7RxZKKasBbAPwGseSElEgiUmxeODhtd7j3h479r3wnoIR0Vg0eg2y\nP52N207fhtV7VyPx1kRAAHO+OMevbnd3Nyw/sKBpbxNcNq5DS0TBm8qY0r1SyifHeD0HwHfgHme6\nXUqpnVyIFM2EELaYmBj9wMCA0qFQiDU3duIzH30WDoc7cUlOicMv925HbFzot76k0JJSovtoN5Ly\n/McB76veh8RPJkL0CsSkxWDWJ2Yh+zPZSFyXGOBKRKQ0nU4Hh8MR9mNK9wybef/0yBellBYAT8Hd\nzU9E5CNzVjIKH8z1Hnd29OHl3x9XMCIKlhAiYEIKADFvx0D0uqcRONocuPGjG6i9tRa1+bW49oNr\nsN3gpEUiCmzSSamU8jkA6XBPdMoZpY4F7glQnMVARH6KP3YbNJqheZDP/e4I7DaHghHRVGnPBO4U\n6zZ349JXL+HQ/EM4es9R3PjvG7A32Wc4OiJSs0l33xPNBHbfR77yf3kJB2rOeI+/+NVCPPrhWxWM\niKZi37590FzSYLllORp/0whH+xh/ZGiA1M2pSP9QOrI/nQ1tAkd5Ec20SOm+DykhxG6lYyCimbft\nExt9jqv+5zDsdraWhjPXEheW/WgZNtVvwurfrUbq/amB7zYuoN3Ujsu7LkPEcOVAominmqQUwFal\nAyB10mjU9GNKobZwcQbu/MAy73FTYxde2cexpeHq7rvvxt133w0A0MZqkbUtC+teXYc76u/Ash8v\nQ8rdKX7npD6QCo3B//e8+1g3Gn/bCHszu/mJpoua7rEz0n0vhNgM9/qlRgBpAaoYAeRwhj6NxB2d\nosPlS83460//0nucmhaPn/32c5yJH6H6r/WjuaoZTXub0HW4Cyt/uRKz/2q2X72L/+8irv/HdQBA\n4q2JSC1KRWpRKlLuSoE2jrcLolBQ045O056UCiG+AaDMc2gdpZoRgGRSSiMxKY0eZf/8B/zZdNZ7\n/Kkdd/t17VPksd20QZukRUyi/yL87655F72nev3KhV4gaUMSUu5KcT/uSIEuTedXj4jGF21J6Z8A\nlEkp949T76KUcum0BkNhh0lp9Ki/3o7tf/kzuJzuz6TERAN+tmc7kpJiFY6MlGC7acOhOYeCrp+w\nJgEpd6Ug7aE0ZHwwYxojI4osakpKZ2IggXm8hNTDNO2REJFqzZmX6rPLU3e3DdX/+66CEZGS9LP1\n2HBqA5b+cCnSH02HNnHsjrSekz2o/0k9Gv+nMeDrzl4nuNoMkbqpZtNiKeXnlY6B1Ik3kujxsU9t\ngunVUxiwOwEAL1ab8ejj65GZFXihdlKfixcvAgCWLp1ax5cQAgmrE5CwOgHz/mYeXHYXOg93ot3U\njo43OtD5Tidcff7bmCYVBP5ZufS1S2ja24TEdYlIWJuAhLUJSFybiPjc+IBDB4iihZrusTPxm3hE\nCLFZSvnaWJWEEK9KKR+YgXgozKjpF4amV0ZmEh57fD2e31MLALDZHPhF5Rv4xt8/rHBkFKwzZ9xr\nzk41KR1Jo9fAeLcRxruNAADXgAvdR7vR8UYHOt50PwZaBkZNSjuPdMLR5oD1gBXWA77TG2JzYoeS\n1NXxiF8ej7hlcYhJZrJKkU9N99iZmn2/HUAK3F30Fill54jXUzzl6dMeDE0bIUQO3JPa2jxFaQB2\nSynNU7gmF8+PMl2dffjsx36Krs5+b9n3n/4YVuXOUTAqCta+ffsAAI899tiMvq+UEn3n+2BYYPCb\nme+yufBG0huQAxO736UWpmJdzbqA7yUE11WlyKCmxfNn6s/AfADbBw/4yxx5hBB5APYD2DKYhHqS\n1DohRLGUkmOGKShJyXH4xKfvwNM/HOpcqfjPA/jB0x/jZweNSgiB+BXxAV9z9jiR/dlsdNV1oedU\nD1y9/t3+gcSkB75FXt19FTf+8wZiF8XCsNCA2EWxvo+FsVyyimgSpj0pFUL8BMCTAI4CqB2lWjqA\nJ6Y7FppWVQAqh7eKSiktnp26qgCkKhYZhZ1HPnQrXnrxGK5eaQUAnDt9E382ncV9RasUjozCkS5N\nh+VPLwcASJdEn6UPPSd6hh4ne9B7vhcYkavGLw+c5PZf7oe9wQ57gx14Z5T3zNLBMNcA/Rw9DHMM\nMMw1YME3F0CjU89C5URqMxMtpQVSykAL5vsQQlycgVhoGgghCgHkANgT4OVqAGVCiK1SyupJXn8q\n4VEY0sZosONL9+Lvv/6ct+ynPzmIjXcuQXw8F9RXs9hYdS/hJTQC8UvjEb80HpmPZ3rLnf1O9J7t\nRd/5PvSed39NviM54DX6r/QHLB9uoGkAA00D7uYYAMIgsPAfFvrV6z3Xi9MfPw19ph66LB10mTr3\n80yd9zF4rE3S8vOQQk5NP1MzkZQG221bNn4VUqkiz1fLyBc8raWDdZiUUtDyb1uMDbcvxpF3LgMA\nWpu78ZufvYUdX7pP4choLEVFReNXUiFtrBZJtyYh6dbxV3pIWJcAZ48T/Vf6Yb8Z3BaohjmGgJ9l\n/e/3o7uuO6hrCIOALkMHXaoOS36wBGmF/u09LX9oASQQY4xBTGoMYowx0KXqoInX8LOUAlLTz4Wa\nphbWKB0ATdrWIOo8CfdWs0RB2/Hl+3C07iocA0NLRG0uWo2lK2YpHBlFs6XfG1pZwNnvhO2qDf3v\n96P/ytDDdtUG200b7DfscPW7oJ8TuIXfVh/8xiDSJmG/YYf9hn3USVsXv3wxcEuuBtAmaqFN0kKb\nqEVMUoz3efLtyVj4zQCtuBd60XuuF9p4LTRxGmjiNNDGaaGJH/Y8TgOhUU9SQ+FtJpLSPUKIr0sp\nvzdOvSoAG2YgHgq9HACQUo62jawV7q1kiSZk3vw0/MUnNuI3P38bAOBySfzo3/+EHzz9cWi1HJtH\nytPGahG/PH7U8adSSjg6HHB2O0e9RmxOLAaaBsasM1JMauDb90D7KCuVuABnpxPOzuDfo+X5Flie\n8usA8yP0widxvfW1WxG70HcYh7PfifOfPw+NXgONQQOhFz5fNfrAZXEr4pC4JtHvPe1Ndjh7nBAx\nYuyHiloBaXwzkZRKAJeFEK8CqEPgyU5pAPJmIBYiCjPFH78NB0xncONaOwDgwtlGvPR/7+GDH+FH\nBqmfEAI6ow46oy7g69mfykb2p7IBuBO3geYB78PebHc/bxr2vGUADqsDugz/60mXnFDSOUibFHil\nAEeXI6jzpV3CYXe4mx8AiBj/RNDV50LjLwPvtjWWeV+dh6X/7r/mreUpCxp+3jD+BTTwJqgagwZ3\ntd3lV6X3XC9OfuSkTzKr0Wnc34cWEFrhbg3WusckC60ANED257KR/pD/SpZXv3sVfef7vPWhGXYN\nz3NofK+lS9dh3t/M87tWn6UPrftavfUh4H0MPxbCfZ3B58bNRsTO9x/f3frHVncyP6z+yAl+SpqJ\npNQMd2Iq4B5XGKjPQYxSThFECGEcozWVKCC9PgZf/loRnvrKXm/ZzyvfwG2bcjB7DhvgKXJoY7XQ\nztcGTCaCIoG1f1gLR7sDDqsDA+3uBNbR7oCzy+l+dDvh6PI9Hi0pdXZNPMEFAE2cfy+Gyza5zEdj\nCNwj4rIHeT2XO2mWdgnpCJxmOHuc6D3VO+HYUrcEXlSm9Q+t6Hi9Y0LXilsaFzAp7TnRg4tfmfg8\n8DW/XxPw5+jCFy+g/7Lv8A6povRrJpJSC9yTncYaM5oO4OkZiIXCkN1uh07n2yqg0fh+UNlsNjQ0\nNODIkSN+5z/00EOIifH9UT927BiuXr3qU5aRkYFNmzb5nf/qq6/CbvedzJCbm4ucnByfsra2Nrz1\n1lt+52/ZsgXx8b5de2fOnPFuxzgoKSkJ9957r9/5Bw4cQHe370SIZcuWYeXKlT5l3d3dOHDggN/5\nd911F1JTfT88L126hNOnT/uUGQwG3H///X7nv/3222htbfUpW7BgAdat811UfGBgAK+88orf+bfd\ndhtmzfIdA3r9+nUcPXrUp0wIgUcffdTv/NraWty8eRMrcpNx7pR7343+vgH8oOxV7P7Bk9AMG8/2\nhz/8wW93knXr1mHBggU+ZU1NTTh8+LDfez3wwAPQ633H/p04cQJXrlzxKUtLS8Odd97pd35NTQ36\n+30/8FetWuW3u5HVasUbb7zhd/59992HxETfrspz587h/PnzPmUJCQnYvHmz3/kHDx5EZ6fP3iTI\nyclBbm6uT1lfXx9MJv85qHfccQfS031bfi5fvoyTJ0/6lOl0Ojz44IN+57/zzjtobHS3hg3+js6b\nNw/r16/3qed0OvHyyy/7nV9QUIDs7Gyfsvr6etTV1fnVfeSRR/w+B8xmM27cuOFTlpWVhY0bN/qd\n/8c//hEOh29L4Nq1a7Fo0SKfspaWFhw6dMjv/KKiIr+VBk6dOgWLxbe7OyUlBR/4wAf8zt+/fz96\ne30ToRUrVmD58uU+ZZ2dnTh48KDf+ffccw+Sk31XB7houYizzrNAMtyPBUBcXBwKCwv9zn/zzTfR\n3u7ufeiQHTi/7zwWLVqEtWvXeuvM/vRsxK2Pw8m6k4AdEDYB2ADYgey0bOilHq4+l/fR29aLnvYe\nvHLgFcDza6TVavHwww9D2ieX+Ai98G7IMFx2R3aA2uNca1gL7vB7gOa8BglImPj1tL4twoP3gLiW\nOMRMML3q6evBn//8Z797gHRO7t/t3SPvwgn3HxU+9wD15J8BzURSagXwHSnllbEqCSF2zEAspCC2\nktJU3HlfFq5d6UFvj/uD9vjRa9j3wlF8iN34RFMzyrDLpFuToF+lx0C6/zjV7NuzkZmZ6VP2/vvv\n4/jx4wGvZZhvwAdsH4DL7oK0SZgPm9Fc3ww4ADEgAAeQmZKJZTnLIO0SLpsL0i4RvzI+wLouQOJD\niUhblwY54G4BlQ6Jns4eNN1sApyAcAp3t7QTmJ0x212mH+Ub1QHOhU7EiBjE6mLd1xt2XbvN7r6W\nC+6kzgUIiIDDFABAuCY+jlWKwNmidE0yixwlhElfb4ZM+zajQojFUsrLoapH6iOEaAdglFIG/DUQ\nwv3bNtrr41yb24yS1ztvXcK3d73gPTYYYvBfP/sk5s7n3gxqodQ2o0Rq4bK53AmtS7qHDzh9n8Pl\nTg69z50SIkYgbnGc37UG2gbQd6HPXVd69qn3JMbe59KTbA57npSfBH2m/4oPba+2wdnjHLqWC5j7\nsblwOKNkm9EJJJqLATApDU8WAHnjjBk1j1JOFLTb71yCwgdzYXrlFADAZnOg/F9ewvf+66PQ6bit\nIxEpT2PQACFK73RpOug2Bp4kNxlpDwTYy+jjIbv8lM3ENqN7pJTbgqhaBffYUgo/JrhXT8jBiORT\nCJEzrA7RlJV8+T4crXsfrc3ucbbnzzbgV8++ic9+4R6FIyMiUg+Hw4neHjt6e2zo7bXD2m5Dc7MD\nLS1OtLY60dbqQnu7hMupnmWzZmJM6VYhRJKUsmu0CkKIn4DrWIazCgA7ARTAv0V0cMBfoC1Ig8J1\n5mi4xKRYfP2bD+ObX92LwdFH1b89gnV5C1CwcbGywZHfRCEimjgpJWz9A+jutqG7qx/d3TZ0ddrQ\n0nL1ATAAACAASURBVGxHY4MDTU1OtLZKtLW5kJXRCCG70NtrR2+PHT09NvT22PFW3YdhsydiwBGL\nAUcCXK7RWlzVc4+diTGlLgB7pZR/EeC1RXC3kOYDkFJK9r+FKSFEBYACKWX+iPI6ABYpZfEkr2vT\n6/V6my34XU8oOvzymTfxu1+/4z1OMcbhv372SaRn+C+0TUSkBKfDha6uPnR09KGro9/9tbMPXZ39\n6Om2oavL/bW7qx/vXzXg9Lk56OqKQW+vHjZ7LAYGYj1JpQEORywk/JfI2rBmL9JTr/qVHzhcApt9\n/G1z3WMN7KoYUzpTSSkAfERK+cKw8s/B3cIm4N4TfauUklu0hDEhRJXnaannaxncE6AmvRk2k1Ia\njdPhws6/+R1On6z3lq1eMwff+eE2ji8lopAbGHCiq7MfnR296OzoQ2dHv/trZx9amm2ov+FCQyPQ\n3KJBW6sGHR0x6O7RYcARB/tAPAYGYmF3xEGv68WmW//X7/oNLcvw3pkPTTiuW1f+HrMzz/uVv1n3\nKXT3ZgRxBfUkpTPRfV8upXxKCPEdIUQH3JNiBltHTQCKpZQdQogtMxALTSMpZbEQIg/AVk/Rbikl\nJzjRtNDGaLDzHx7Blz7zK3R3u/9oOX2yHhX/+Rq+9NVJ/x1ERFFCSon+vgG0t/fC2tYDq7UX1rZe\ntLf3oKO9F41N7kSzqQlw2Frgcrb6XeO9M4+iuX0lnM7g8zmXM3DqpYvpD1g+7vUQh9S0eMTHGxCf\noEd8vB7xCQZcb3Wgra0biYlOJCW5kJwikZIMpKQIpKYJpKVpkJ6uxd9+BXAEt3nXtJv2llKfN3OP\nHd0OoAPAdinlczP25hSW2FJK43n30CX801MvYPhH2VdKH8ADj6wd/SQiilgDA060t/WgtaUbbS3d\naG3pRntbL6ztPWhv60VbWx8abjrR2KRFV3ccbLYk9NsTYLMnod+WCJs9Ef22RDicQ5sjLF/0OnLm\nv+v3XkdPfxCNrcv9yscihBP33/kDDJ8ukZBowIBzNn7/p6GRbrGxA0hMcCAx0YmUZImUFInUVCA9\nQ4P0dA2ysrSYNSsGd96lxYi9XCbEYDDAbo+ellIvKeXnPbOxvyOlfG34a0KII1LKDTMZDxGFv9s2\nLcFffuZO/OqnQ7tp/fj7Jsydl4o16/y37aPpNbj72cidsYimyulwwWrtHUo2W7s9zz0JaGs3mpt7\n0NykQZ8tGcakeoycJ9vXn4TXj2wPODZzLDZ7fMByvS647UljYx1ISXbCaHQhPR349//6S6RnGJCQ\nFIv4eD20Wg1sNuDCBSAjA0hLA/R6HYDQLQcVDkLWUuqZtBSMVACVAD6HoXVJ0wHUSim5JBT54OL5\nFAyXS+Jfv/UiDr0xtHVrUnIsvv/0xzBvfoB1+WjacPF8mgynw4W21m40NXWhuakLLU2daGrsQktT\nF1qau9Da0g1rey8cDg36bEno709Gn83z6E9Gv+d5vy0Jg3Om79v43zDofZNGp0uLmrf+34TjW73i\nKj758RNITolDSkockjxfX63JxBtvJmJ2tsC8eVrMnavBrFnArFlAVhaQmelOMA2Kt0GOTqfTweGI\nvMXzLQh+V1UBLqZORCGi0Qh87ZsP4et//VtcudwCAOjq7Mc/7Hwe33/6YzAaA7dyENH0k1Kis6MP\nzU1daG7sQnNTp/t5UxeaG93PW1u74Rpln/ebzSvw/o3N6LOlwGYPfnWNvv5kv6Q0NhYw6Ptgs/vv\nngQAiQkuzM52Yd48gfnzNcjOFpg9G1i3bgE2b17gV//u+4IOh4IQyqTUCvfEpSOTOFcA2B3CWIgo\nyiQkGPDt8ifwlZL/QXtbDwDg5g0r/nHn89j9H08iPp7dyUTTpaurH431HWi46X40Dvva1NgJm21o\nJo2UAn22ZPT2paC3PwN9/UvR22eEXt+L3KX++6w4nTpYu+ZOOKY77t2Exx4bQHp6AlLTE5CamoD4\nBD2+9jUBhwOYO9f9mDNn6HliogaYYNc+hU4ou+8vAiiUUl6Z5Pm1UsqCkARDEYPd9zRRF8414Btf\n/h1s/UM3wVvWz8c/lz8BgyG6xmcpgd33kclmG0DjzU6/hHMwCe3p9p+MOuDQo6cvDT29ab5f+4yQ\n0r9NLC7Wins2POtTFhOjwYBrGV49OPbPk9EosXChwMKF8D4efRRYsWJq33c0UFP3fSiT0i1Syv1T\nOH+9lPJoSIKhiMGklCbj8NuX8M9/938+3YEb71iCv/uXD3IN02nGpDR89fXaUX/Divob7ai/bsWN\n6+2ov96OmzesaGvtmdC1zKc+jKa2pRM6R6Nx4be/OorsuYnIzEpGVlYSjGkJqK8XKCgAFi+GT9I5\n+FiwAEhOntBb0TARmZQSTQcuCUWT9dqfTuO7//qyT9ntdy7Brm8/Br1+RhceiSpNTU0AgKysLIUj\noUD6++y4Wd+BG9fcCeeN6+3uRPR6+7iJp90eh67eDHT1ZKC7NxNdPRmYk3kGC+f6tyedPH8/rjfe\nEnRcGo07uXz7bSA7e8LfFk2BmpaEmnJSKoRYDyBfSvnsuJUnfu1FAEqklLtCfW0KD0xKaSr2PX8U\n//0fvh04t23Kwd/9yweZmFLEcjpdaGrsxLX3W3Ht/TZcv9rmafW0orWle9zzHU4dunvS0dWbie6e\nDHR7ElH7QIJf3ZwFF7D1w2bMyk7B7GGP6ucz8I//5JvjaDTu1s6VK4Hly4ElS4YeCxcCXEVMGRGV\nlAKAEOIbAArh3iq0a8oX/P/bu/Mwqao7/+Of0ys7vbCDLN0Qo2iIDbhEY6KAitGJRsBE42RVzDLO\nz0mUkExGTZwhaEzMxESBmGTiJBFBjcFJlEVRcWPpxLihQAMtiyy9sEqv5/fHvUXX1t1VRVXdW1Xv\n1/PUU12nzq37penb99tndT7zKknfkTTFWnswGZ+JzENSihO19I/r9OD9z4WUfXTiSH3/P69g8hMy\nWlNTi3bUNui92nrtcBPQ92rrtfO9BjU3x7dFj7XS6+9O14FDQ3Tkg1LFOtlnwgTp73+PLF+zRlqw\nwElAA4+xY/29NFKuyrqkVJKMMTMk/UjS/ZIWJZpIGmMulJOMWkmzrLUHkhIgMhJJKZLh8UfWa+F9\nq0PKxp08WD+46zMqKY1s/QH85EDj0eMJZ6D1873tddq756DiuYU3NffSwcOD1af3fvUs7mg/6tOn\nWMNGlOoPj12h/XWxL7lUWipNmiQ9/bQiFqlH5sjKpFSSjDElchbGv0rSUkkr5CyKH+XvqOPHjJZU\nJWmypJlyFte/3lr7WNICQ8YiKUWyPPFotR74WchGcho6vER3/OhKnTSKfTvgvSNHmlS7rU7bt+7X\ntpr92r7VeTTUx7ZrULDWtkIdPDxYBw4NVePBoTp4ZKg+ONZXkvSFz2/RdZ8/pmEjSjV8RKn69XfW\n7Pzc56SHH478rJ49pfHjpdNPl047zXmcfro0ZAjJaDbwU1Ka1EFV1tpGSbPcrURvkHSXpBJjjJWz\njml9UPUySSXu10bOGqdzrLWPJjMmAJCkT19Vpb59e+gn855SW1u7JGcd05u/9gd97wf/pDMmjfI4\nwuxQW1srSRo5MnKhcTiam1r1Xm19SPK5bet+7X0/sZFq1hodPlqmA4eGqqltlA4cGqr9df3U3h69\nC76oZ6WmXBxZfsYZzkSjM85wHhMmOAnomDFSPotWIA1SPvveGDNGTktomaRKOYloo6Q697nmRJaS\nQnZjSSgk27pXavSf//HnkHVM8/KNvvCV8zTjmjOVl0fTz4lgSagObW3t2rWjQdu27tf2Gifx3L61\nTrt2NKi9PfF7b0FhvkaMKNXAoYP02J/P0dbt/XT0aOxZ46RJ0roo29y0tzuTkZBbWBIKiBFJKVJh\n8zt7dPvcxyNmIk86a7S+9b1L2Zb0BORqUnr0aLO2btmnmk17VeM+b6vZH/eEo2DFxQUaOaZcAwYN\nVUVlicadXKqTRpZp8ND+ys/Pk7XO8kl79nT/WR/6kHTmmdLkyc7jnHMSDgtZhqQUiBFJKVKlbv9h\n3TH3cW16J/SOXj6gj+b8x6d0+kdP8iiyzJbtSam1Vvv3HtKWzfu0dfNebdm8VzWb92n3zsaEPzM/\nP08njSzTqIoBGjVmgPr0HazanYP099d6ac0ao40bnbGeV18deezMmdLSpaFl5eXSWWd1PCZPlsrK\nEg4PWc5PSSkL9QHISeUD+uju+z6rhfet1l+eeO14ed3+w/rO/3tEV193lj573dmsZ5rDWlra9N72\nOm3ZtDekFfTQwWMJfZ4x0tBhJRpVMUCjxzgJ6OiKAerZu1QvvZSvVaukhx6R3nor8tgXXoielF54\noVRbK519dkcSWlHBBCRkJlpK4Wu0lCIdnlu1UT+7e7k+ONocUn7SyDLddMtFOm3CCI8iyzxPPvmk\nJOmyyy7zOJL4HDr4gWo273Mfe1Wzea9qt9WptbU9oc/r17+nKsYOVEXlQI2uHKjRFQN00qhy9ehR\nqKNHnXU8n3lGWrVKqq52xnN2pbP1QIET5aeWUpJS+BpLQiFddu1o0H/dtkxbNu2NeG/65R/Rl288\nX3369vAgMiRTe7vV+7sPdHS9b9qnrVv2au+exPZ9MUYaNrxUFeMGqmLsICcRHTtI5QP6yERprmxu\nlgYOlA7GONF+0CDp4x+XPvEJ6ZvfpAUUyeenJaFISuFrJKVIp+amVv3Pohf0p6XVEbOjS8t66Zov\nfkyXXHa6CgpYHycTNDW1aPvWOqflMzABafO+iBbxWBUXF2h05UBVjh2kMWOd59EVA9Qzys5ghw5J\nLS3Rx3JefLG0fHn0c5SVSRdc4HTLX3ihdPLJJKJILZJSIEYkpfDCuxvf18/uXq6aKK2mw4aX6AvX\nn6ePX3By1JYweKOh/ojb7b7v+POO2vqEl14qK+8d0vJZOW6Qhg4vUX5+52smbdwoLVsm/eUvTvf8\nd78r3XFHZL1775Vuvtn5undv6fzznQR0yhSnm55lmZBOJKVAjEhK4ZXW1jb9ackG/e+vX1JTU+Sy\nPuNOHqxZ156lcz4+tstEBcnV1tquHTvqVbPJST63bnGeE9n1SHLWqD1pZJnGVA5SZVAXfCzbz7a3\nS6+8Ij3xhPSnP0nvvhv6/uTJ0tq1kcdt3Sr97ndOEnrmmVJRZEMrkDYkpUCMSErhtd27GvXbhS/o\n+Wfeifr+kGH99ZlZkzRt+nj16El2kUxHDjcdTzprNu9TzZa92l5Tl/Danz17FXW0fLpd8KPGlKu4\nuDDmzzh2zJmc9Kc/Oa2iXa0Raozz/sCBCYULpAVJKRAjY0xTYWFhUXNzYmPAgGTZ9M77+vUDz+vv\nG2qjvt+nbw9deNEpmjb9NI390OA0R+cf69ytgiZPnhzzMa2tbdq5o0G1W+u0bet+bduyXzWb9+r9\n3QcSjmPQ4L4aUzlIFeOcsZ8V4wZp8JD+Ce/Y9X//J/3mN9JTT0lHjnRdNy/PWZx++nRp9mxpwICE\nTgmkRVFRkVpaWkhKge6wJBT8pnrdNj306xe18c3dndapGDdIUy8er3M/MU6DBvdLY3Te62rx/La2\ndu3ZfUDb3eRzu/vYUVuf8NJLBYX5GjW6/HgLaMXYgRozdpD6JnmlhJtvdsaCdqZPHycJveIK57m0\nNKmnB1KGJaGAGJGUwq/een2nHn14vV5es0ld/Rodd/Jgfez8cTrznAqNrhiYcEtdpli2bJlaW9tV\n9dFztaO2Xjvfa1Dt9jrVbnMe0cbnxqpf/56qHDdIYyoHHn8+aVRZ0lZD+OADZ+H6iRMj33v5Zelj\nHwstGzJE+vSnnceFF0rFnt/SgfiRlAIxIimF3+18r0F/fqxaz654u9udfvqX9NSEM0ZqwsSROmX8\nMI0cXZ6xk6RaW9u0b+8h7drRqJ3v1WvnjgbtqG3Qpnd36tCBE7tejZGGn1SmisqBqghKQsvKeyd9\nxYP2dmn1amfi0WOPOV3ve/ZEJpjt7dLo0VKvXtKVVzotopMnM1MemS9nk1JjzCOS+ltrL07bSZHR\nSEqRKZqbW/Xqi1u0/C+vq3rd9piWIurRs1BjPzRYlWMH6aTRZTppVLlGjipT/5Jeni43Za3V0aPN\naqg/on17D2nv+we15/0DzvNu5+u6/YcTXm4p2KDBfTXK3XJz1Bhn+80Ro8rUo0fsk48S8eab0kMP\nSb//vbRjR+h7y5ZJ0Tak2rfPGR/KSmDIJrmclD4gaYG19m/u6/uttV/rpO6F1tpn0hYcfImJTshE\nDfVH9PKazXr5hU36+4bauMdL9u3XQ0OG9teAgX1VPqCPygf2UfmAPurbt4d69yl2Hr2L1aNnofLz\n81RQkKe8/Dzluw9rrVpa2tTa2q7Wlja1trSpqblVRw436eiRJh0+3OR+3azDh4/pQMNR1dcdUUP9\nEdXXO89NxxLvZo+mtKy3RlcM0MjR5Rpd4SSgI0eXq3fv9N0H9+yR/vhHJxmtru683nXXOS2nQC7I\n2YlOxpjlkq6y1h5yXy+21l7dSd2naVEFS0Ih0x053KTq9dv09w21em1DrXbuaEjp+YxRl2NcU6mo\nqEDDTyrteIwo1YiRZRpxUqn69uvpSUytrdJf/yotWuQsat/W1nndESOkz31OuuYa6aMfTV+MgJdy\ndkkoY8z1kh6QFPgbtUpSTSfVK6y17OWX40hKkW327Tmot97YpXff3q2Nb7+vLe/uOaHJP+nWr39P\nDR7aT4MG99PgIf01ZGh/jRhZpuEnlWrAwL6+msh18KA0fnxk93ywPn2kq66S/vmfnf3l87nrIMfk\nclL6iPvlCvd5tqQFUaqWSppjrS1PS2DwLZJSZLv2dqu97x9wZ6jX673tddrxXoP27zuk+v2HE14q\nKV55+Ualpb1VVt5bpWUdj7Ly3ho8tL8GD3ES0Wj7vPvZBRc4E5mC5eVJ06Y53fRXXOFs9QnkKj8l\npQVpPl+NtfY7gRfGmGnW2kXRKhpjWBYAQNbLyzMaMqxEQ4aV6MxzKkPea2+3OnjgqPbvO6z6uiM6\ncrhJR44440EPHzqm5uZWtbW1q621XW1t7Wptdb7OyzcqLMxXQUG+CgrzVViQp8LiAvXu7YxFDR6X\n2qtPkfqX9FK/fj191coZj9paqbFR+shHIt+7/vqOpPSUU6Qvf1m69lpp6NC0hgggBulOSvuHve40\n8bTW3p3iWADA1/LyjEpKe8e0D3uusVZ6/nnp5z+XHn9c+vjHI1tEJekzn5FuvNFpFT3nHGbOA36W\n7qQ05u54Y8xXrbW/SmUwyAzt7enpvgRw4p588klJ0mXR1lRKgqNHpT/8Qfrv/5Zef72j/LnnpH/8\nI7K1tEcP6f77UxIKkBX8dI9Nd1K6whjztJzJTgckVRhjviIp2t+ucySRlAJABknVPIXt26Vf/EL6\n1a+khk4WMLj/fhJQIJOlNSm11i4yxoyV9KBbVCKps2768K5+AECOee016e67pYcf7nw5p8pK6Zvf\nlL70pfTGBiC50t1SKmvtHDmtoDLGPGKtnRWtnrumKQAgx1grPfOMdNdd0vIu7gQXXyz9y79I06ez\n3SeQDdKelIaZ18V7M9MWBQDAN44dkz77WWn//sj3+vSRvvhFp2X05JPTHhqAFPI0KQ1sN9rJewfS\nGQv8K48mECBjnHHGGSf8GT17Oi2gt93WUTZ8uHTzzc4ST/36nfApALj8dI9N6+L5UQMw5qtyFtGv\nsdZebYyZImmqtXaup4HBF1g8H8heTU3SG29IEydGvldXJ40cKY0eLd16q7P9Z1FmrdsPZIRcXjw/\nhLvtaGBXp0pJstauMsbUG2O+ba39sZfxAQCS79gxZxb9j37kLPG0bVtk62d5ubR2rbPgvY8acgCk\nkNeXepW1dpK7q9Px0UNut35l54cBADLNsWPOYveVlU73/M6dzvJOP/959Prjx5OQArnE68u90ePz\nAwBSrKVFWrjQSUZvuknatSv0/Z/+1OnKB5DbvE5KS6MVGmNGK47dnwAA/rBnzx7t2bNHktTeLj3y\niNPiOXt2ZDJqjDRrlrM9aLHno9kAeM3rJaFWGGPWyR1Taoy5UFKVpLmSpngaGXzDT1ugAeja2rVr\nZa3Uo8flmjtXqq6OrGOMdPXV0r//u5OwAvCOn+6xXi8J9agxplzSQrdotiQraZa19u/eRQYASMSW\nLf3061+PD9mXPtiMGdIdd0innpreuAD4n9fd97LWLrTW5kmaJmmitTbfWvuo13EBAOJXV9dTr78+\nIKJ86lRnNv2SJSSkAKLzPCkNsNau6moxfQCA/02evEcTJuw7/nrSJGnFCucxebKHgQHwvZQnpcaY\nfsaYxcaYOvfxsDGG/TgAIIO1t0uNUdZPMUb68pff1CmnOK2ia9c6raQA0J2U7uhkjOkvab0i1xzd\nLGmStfZgyk6OrMCOToD/vPSSs87o8OHSn/8c+l5LS4skqaCgUMZ4EByAuPhpR6dUt5QukpOQLpQz\niWm2WzZAHZObAAAZYN8+6ctfls4915lVv2yZtGpVaJ3CwkIVFpKQAohfUpNSY8y8oK/HyFneqdRa\ne6O1dpH7mC2pQtIkdz3SiGMBAP7R1iY98IB08snSb34T+t7NNztd+QBwopLdUho8cugqSTOttQfC\nK1lrGyXNkjQjqLgqybEAAE7QunXS2WdLX/uasyVosNNPl+69l61AASRHsn+VVAR9XdnVbHprbbVC\nx5pWdFYXAJBejY1OInrWWdL69aHv9e3rbA1aXS1deKE38QHIPslePL/UGPOUpAOSqowxi7upH6hT\nIpJSdCKVk/EARHr8cekb35B2745873Ofk378Y2nYsOjHvvbaa5KkCRMmpDBCAMnip3tsKnZ0uijo\n6/BZ99EE6vjnuwJf8dMFA2SzPXucZPTRKNuXfPjD0i9+0X3LaG1trSSSUiBT+Okem4qkdLakBkk3\nyNnTvjMmqE65pPtTEAsAIEYtLc4i98F69pRuu82Z0FRU5E1cAHJDspPSGmvtIkkyxkyVtKKrtUiN\nMVMDW4oaY25NciwAgDiMGCHddZd0443O66lTpQULpAoGVwFIg2RPdKoJ+nqFnJbQqIwx35a0IfhY\ndnoCAG9df7105ZXSr38tLV9OQgogfZLdUnq8u95a+6gxpt4YY6219wRXMsbcIuk71tryoOL57PCE\naAyrcANJtXu3dPvt0j33SH36hL6Xlyc99ljin11eXt59JQC+4ad7bKq3Gb1B0gNyJjFVK3SW/Y2B\nrn6gM2wzCiTX4sXS178u1dc7Sz798pdeRwTAS37aZjSlSal0fKzoj8KK51hr707piZEVSEqB5Kir\nc2bWLw5bqO/pp6WLLop+DIDs56ekNBWz70NYa+8yxqyQdLVbtLirRfUBAMm1erX0+c9LO3eGlpeW\nSkePehISAERIeVIqSW4SSiIKAGnU2ir94AfSnXdK4Z1il1wiPfhg54vgA0C6pSUpBQCkV22tdM01\n0osvhpb37i395CfOLHsfzW8AAJJS+J+fdpsAMsGjj0pf/aqzf32wiROlP/5RGjcudedevny5JOki\nBqoCGcFP91iS0iQzxsyXVCVpvqT11tpGY0yVpKmSJltrZ0Y5psKtX+8WlUmaZ62t7uQcGV0/Xn66\nYAA/++ADZ+elBVH20vvWt6T/+q/U78rEpEQgs/jpHktSmnwlchLQqVLI+l+NkiaGV3YT1lWSpgSS\nODfJ22CMmWmtXZlN9QGkzhtvOONEgw0cKP3P/0jTp3sTEwDEKtk7OsFRIycJDXx9l7W21FpbE6Xu\nEkkLg1sV3Xrz3PeyrT6AFJk8WfrxjzteT5sm/eMfJKQAMgNJaWrMdJNQY62ttNbOiVbJGDNVzmYC\ni6O8vVRSiTFmRrbUB5B6N90kXXutNH++9NRT0pAhXkcEALEhKfXWNPc5ogU1qFV1WlBxptdPiJ+2\nQAP8oq5OamuLLDdGeugh6dZbnS1D0+3UU0/Vqaeemv4TA0iIn+6xJKXeiqUVcVYW1U+Iny4YwA/W\nrpUmTHD2r4/Gy0umsrJSlZWV3gUAIC5+uscy0SlF3O1Vy92XFYo+G71Ckqy1YQu3HNcoZ+JUttRP\nSHNzswoLC0PK8sKagJqamvT+++9r3bp1EcdPnz5dBQWhP+qvvfaaamtrQ8oGDBigc845J+L4p59+\nWs3NzSFl48ePV0VFRUhZfX29XgxfFFLSlClT1KtXr5Cyt99+W5s3bw4p69u3rz75yU9GHP/ss8/q\n8OHDIWXjxo3Thz/84ZCyw4cP69lnn404/rzzzlNpaWlI2ZYtW/TWW2+FlBUXF0ddxuell15SXV1d\nSNnIkSM1YcKEkLKWlhY99dRTEcefeeaZGjx4cEjZjh079Le/he6nYYzRZZddFnH8+vXrtXv37pCy\nIUOGaPLkyRF1n3zyyYiZpBMmTNDIkSNDyvbu3atXX3014viLL75YRWHT019//XVt27YtpKysrEzn\nnntuxPErVqzQsWPHQspOOeUUjR07NqSssbFRL7zwQsTxF1xwgfr06RNS9s477+jdd989/vrZZ4fr\nvvsmqKUlX3feKZ11lhT4tj333HM6ePBgyPEVFRUaP358SNkHH3yglSsj5yB+7GMfU3l5eUjZ1q1b\n9cYbb4SUFRYW6pJLLok4/pVXXtG+fftCykaMGKEzzjgjpKytrU1/+ctfIo6fNGmShg4dGlK2a9cu\nbdiwIaLupz71qYjfA9XV1doZtm3VoEGDdNZZZ0Uc/9e//lWtra0hZaeffrpGjx4dUrZ//369/PLL\nEcdPmzZNPXr0CCl78803VVMT2nHUv39/nX/++RHHr1q1SkfDttM6+eST9aEPfSik7ODBg3ruueci\njv/EJz6hfv36hZRt2rRJGzduDCnr2bOnpk6dGnH8mjVr1NDQEFI2evRonX766SFlTU1Nx5f5Cnb2\n2Wdr4MCBIWXbt2/XP/7xj5Cy/Px8XXrppRHHr127Vnv27AkpGzZsmCZOjJgPrGXLlkWUVVVVNaWV\nUAAAFmNJREFUafjw4SFl3AMSuwe0t7dHHOMVktLUWCRnNnqjdHyGemA2+lJvQwOQidrapIceOkWP\nPRaa4H7jG87e9ale6gkAUs34aX2qbOBO7KkOn2lvjFkiaaq1tjSozEqStTZq27kxpkFOS2Opu95p\nRtePVqc7xpimoqKiItY+RC47cMDZnSm8cfHDH5aeeEIKa1wDgJgVFxerubm52Vpb7HUsjClNMmvt\n0k6WflonZzb6remOCUDm2rxZOuecyIT00kulV17xX0La0NAQ0S0MALGg+17HWzEruq0YabG19q4Y\n6wYS1ciBcN2It5Ux0+uH89N4FyCdnnlGmjFDCs/xbrlFmjdPys/3Jq6urFmzRpJ0+eWXexwJgFj4\n6R5LUuqYo8Qm5IR30ZdImtTJLkaBxKwqrKyr84a/l+n1AcTod7+TvvIVKXguTnGxtGiRdN113sUF\nAKlCUqqQNTVP1BJJU40xc6K0oAYStOBz1UiqMsaUdNGaWJ1F9QF0w1rphz+UbrsttHzoUOnxx53Z\n9gCQjRhTmnyNip6IBYYHBL+3Muy949z944PrZEN9AN347W8jE9KPflRat46EFEB2IylNrhVyloKK\nlogFdjZaEFQW+HpSlPqBbv7gLTwzvT6Ablx7rXThhR2vL7lEev55KWxJRgDIOiwJlWTGmBXW2mlh\nZVPlJKwR3frGmAVyxqFODCvfIKnGWjszm+rHiyWhkIsaG6XzznNm3f/yl1LY3hG+FljkO3wTAAD+\n5KcloUhKk8zttl4gZ3xpjZwWw7lyEtKFnRyzxP1yjvs8X1JJeHKbLfXjQVKKXNXYKPXv7+2WoQCy\nH0lpDnBbR6vkjCFd392ySO6uT4G94FZG2ZI0q+rHiqQU2WztWqmkxH9rjQLIHSSlQIxISpGtVqyQ\nrrxSGjBAWrNGGjHC64gA5CI/JaVMdAKANFu6VPrUp6QjR6Tt26WLL5bq6ryOCgC8RVIK36M1H9lk\n0SJp1iyppaWj7Ngx6dAh72JKpo0bN2rjxo1ehwEgRn66x5KUwvf8dMEAJ+InP5FuuMFZID/gIx+R\nXnxRGj3as7CSatOmTdq0aZPXYQCIkZ/usSSlAJAG8+ZJ3/pWaNm550qrV0tDhngSEgD4CkkpAKSQ\ntdIdd0jf/W5o+fTp0vLlUmmpN3EBgN8UeB0AAGQra6Xvfc9pJQ12xRXS4sVSUZE3cQGAH5GUwvcM\nq4cjA1krffvbzjjSYLNmSf/7v5m1S1M82MkJyCx+useyTil8jXVKkan+7d+kn/40tOzzn5d+8xup\ngOYAAD7BOqUAkOXOPTc0+fzSl6Tf/paEFAA6Q1IKAClw1VXSY48540ZvvFH61a+k/HyvowIA/6L7\nHr5G9z0y3RtvSOPHSz4atgUAx/mp+56OJABIodNO8zoCAMgMdN/D92jNh9/de690991eR+EPq1ev\n1urVq70OA0CM/HSPpaUUvuenCwYI98AD0s03O183NUn//u/exuO1Q4cOeR0CgDj46R5LSykAJOiR\nR6Svf73j9fe/L/3xj97FAwCZjKQUABKwfLmz7mhwI8Oll0qf+Yx3MQFAJiMpBYA4vfqqk3y2tHSU\nnX++tHSpVOz5/FUAyEyMKYXv+WkLNODtt50W0SNHOsomTJD+/GepZ0/v4vKLsWPHeh0CgDj46R7L\nOqXwNdYphZ/U1jo7Ne3Y0VFWWSmtWSMNGeJdXACQKD+tU0r3PQDEYP9+6aKLQhPSIUOcsaUkpABw\n4khKAaAbH3wgXXaZ9M47HWX9+0tPPy1VVHgXFwBkE5JSAOhCe7szy/7VVzvKevSQnnxS+shHvIsL\nALINSSkAdOHAAWnbto7XeXnSkiXSeed5FpKvHT16VEePHvU6DAAZiKQUvtfe3u51CMhhpaXS889L\nl1/uvL7vPqcrH9GtWrVKq1at8joMADHy0z2WJaEAoBu9e0uPPy4tWyZdcYXX0QBAdqKlFABikJ9P\nQgoAqURSCgBBdu2SGhq8jgIAcg9JKQC4Dh1ydms65xyppsbraAAgtzCmFL6Xl8ffTki99nbpuuuk\n115zXp99trN16NlnextXpjn33HO9DgFAHPx0jyUpBQBJt98uPfFEx+vGRmfRfMSnrKzM6xAAZCj/\npMcA4JElS6Qf/jC07P77pQsu8CYeAMhFJKUActrf/iZ94QuhZTfdJH3lK97EAwC5iqQUQM7au1f6\n9KdDu+mnTJHuuce7mAAgV5GUAshJLS3SjBnSe+91lFVWSosXSwWMtk9YTU2Nali6AEAC+NUL3/PT\nFmjIHrfeKr3wQsfrPn2ciU7l5d7FlA3efPNNSVJFRYXHkQCIhZ/usbSUAsg5Dz8s3XtvaNnvfy+N\nH+9NPAAAklIAOeaNNyInMX3/+9I//ZM38QAAHCSlAHLKW29Jra0dry+5RLrtNu/iAQA4GFMK3zPG\neB0CssisWdKYMdJVV0n5+U63fX6+11Flj6KiIq9DABAHP91jjbXW6xiAThljmoqKioqampq8DgVZ\nZv9+Z0moU0/1OhIA8E5xcbGam5ubrbXFXsdCSymAnDRggPMAAPgDY0oBZLW2Nq8jAADEgqQUQNZq\nb5euvFKaOzd0chMAwH9ISgFkrZ/8RFq2TPrRj6RPfjJ09yYAgL+QlML3mIyHRLzyitNCGvDii84D\nqfXyyy/r5Zdf9joMADHy0z2WiU7wPT9dMMgM9fXS1VeHdtl/9avSZz/rXUy5Yv/+/V6HACAOfrrH\n0lIKIKtY6ySgtbUdZaedJv3sZ97FBADoHkkpgKzy4IPS4493vO7VS3rkEecZAOBfJKUAssa770r/\n+q+hZffdJ51yijfxAABix5hS+J6ftkCDf7W0SNdeKx092lE2c6b0xS96FlJOGjlypNchAIiDn+6x\nJKXwPT9dMPCv22+X1q/veD1ihPTAAxI/Puk1YcIEr0MAEAc/3WPpvgeQ8Z5/Xpo3r+O1MdLvfieV\nlXkXEwAgPiSlADJaY6N03XXOrPuAW26RLrjAu5gAAPEjKQWQ0Xr3dpLSPPe3WVWV9MMfehsTACB+\nJKUAMlphoXTnndLq1c4s+9//Xioq8jqq3NXa2qrW4F0LACBGxk8r+QPhjDFNBQUFRS0tLV6HggzQ\n3t7RYgpvLFu2TJJ0+eWXexwJgFgUFhaqtbW12Vpb7HUs/PoGkDVISAEgc/ErHEDGOXbM6wgAAMlG\nUgogo+zeLVVWSvfcI7W1eR0NACBZSEoBZAxrpdmzpV27pG9/Wzr/fGnbNq+jAgAkAzs6wffyGCgI\n10MPSe48GknS228z095vJk+e7HUIAOLgp3ssSSmAjLBjh3TTTaFlP/+5NGyYN/EguiFDhngdAoAM\n5Z/0GIiuqLm52esY4DFrpa99TTpwoKPsyiula67xLiZEV1xcrOJiz1eWARAj9x7riz4nklIAvrdk\nifTkkx2vy8ul++939rgHAGQHuu/jZIyZYa1d2k2dCknzJdW7RWWS5llrq6kPxKehIbLb/mc/kwYP\n9iYeAEBqkJTGyBhTJWmupBnGmFJrbWMX9VZJmhJIytykbYMxZqa1dmUu1wfidcst0p49Ha8vvphu\newDIRnTfd8MYs8QYs0HSVEkVMRyyRNLC4FZCa22NpHnue7leH0lSW1ur22+/XbW1tV6H0qUTiXP1\naunBBzte9+olPfBAbN32qfz+JOuzT/RzEj0+3uPiqd/e3q729va44sklmXDdehVjJlyzyfisRI5P\n1zGes9byiPEhJ8mykko6eX+q+35VlPcq3Pdm5Gr9BL/n1vkxRbg1a9ZYSXbNmjVeh9KlROM8etTa\nceOsdaY5OY977kn9edP52Sf6OYkeH+9x8dQvKCiwBQUFccWTSzLhuvUqxky4ZpPxWYkcn8pjgu6z\nnudZtJQm1zT3uSb8Deu0HgbXycX6QMzuvFPatKnj9cSJkWNLAQDZg6Q0uWbEUGdWDtcHYmKtM8Ep\nID9f+tWvpAJGwQNA1jJOyy1iYYxZIicRizrRyRjj9DdbG3XEmzGmQU7Xv8nF+okInKM7Re62PtHG\nskXbrSLQVRB2LpkogxWj1e1sB4x0nT9Qr6WlRYWFhcdfp/P80epGO78xJiLOWGO1VmppkfLyrPLy\nOmKI9/yx1o31eyVJra2tEf+mE/1exft/Za0N+d7Gev7O/k3p/LmO9TOz7fyB/7OCKH9hpfJnJZ66\n4T9X6fpeSdGvq2ScP9FrpbNY29raVFBQEFE/tt9rVm1tbRHlqfwd3NXvgNbW1kCZ54vs0e6ArNCc\nuyvs57e0tLR5HUQMTijO9nbnkch5JaXq+5Os7/2JxphoHPGeN5Xfy1yT39ra6vfvpVe/W1J53mR+\n9ol+Vn57e3u8xydyDcZyjC8WzpdISj1hjCmJ1tJK/Uh++MsNAACkHmNKAQAA4LmsbCl1x37GsqZo\nuMXW2ruSHU+4eFsNc60+AADIPVmZlEqaI6kkgeMiljaKU2M35w1/L9fqAwAARJWVSWnQGpnpViOp\nqpsxlNU5XB8AACAqxpQmV2Cf94ihA+5+8MF1crE+AABAVCSlybXAfZ4U5b0q93lxDtdPOmNMlTFm\niTFmfuCRyvMBSA5jzK3u+H8APmaMmeHeZ5cYYza499qUDM9j8fw4GGNWyNnvvbKzIQLGmAWSJllr\nJ4aVb5BUY62dmcv1k8ltjd0gaUxg+ICblFak8rwAEuPeyBbJGfozQ1I11yrgX8aYGyTVW2uXuq9L\nJK2SpPD7flLOR1LaNWPMrXL2b5+k0Ik71ZLqJc0MH08Z9Nf/HPd5vpydjaLuA59r9ZPFTYjLoiTK\nVl384QDAe+4f+Y0kpYB/uQ1MK621c4LKpkpaISf/WZrU85GUJq6rCT7GmCo5raqS8x/a5YSfXKuf\nDMaYLZKWBl8sbnmDpHnpWN4LQGJISgH/c6/TGmvt7KCyCklbJN0Vfv89UVk5+z5dulp/003KYk7M\ncqW+MWZGd39ZuT/w8+W0REtSmZwkM/zzKyTVRfmIekmV3QYNoFtJvmYBpFgyr9lOejwDE5nXnWis\n4UhKkRZuS+pcSTOMMaXdtDCvkjQlcHEExo4aY2Zaa1e6ZYGhFNE+p1GJbZ4AwJXsaxZAaqXxmp0t\np5cjqV33ErPvkWKB2XpyuvZjSRSXSFoY/NeaOzZ0nvteLFi0H0iQR9csgASl85o1xsxwz5P0SU4S\nSSlSzFo701o70R3f2eXEI3fwdIWiLyO1VFKJe0EED52IloCWKXoLKoBupOqaBZAa6bpm3RbW+ZIm\npmoiMUkp/CQwdiXihz3oAgge31IjqTzK55QoBWNdAESI95oF4K2Erlm3e3+RghLSoE1ykoakFH4S\nS4vKrKCvV6rzrgrGsQGpF+81C8BbcV+zQZOipgStCR684k7SMNEJflIhdbmqQaNCu+vnyxmYXRH0\nl9utStOSVADivmbDMf4bSK+4rlk3IV0hp6FnrjFG7vuT1LE2edKQlCJjWWtrjDFTJM03xtQHldNd\nCPiUu+lFhdxWlmjrIALwjSVyrtcbory3PtknIylFxgnetMBtEWXxbcDHwq5Zkk/A5wLXbCq2Eu0K\nY0oBAADgOZJSZJyudtIC4D9cs0Bm8eqaJSmFn3R3ETApAvAXrlkgs/j6miUphZ8EZtB3dVEwqx7w\nD65ZILP4+polKYWfBNYWjVh7NGiRXtYfBfyDaxbILL6+ZklK4ScL3OdJUd6rcp+jbY0GwBtcs0Bm\n8fU1S1KKdAp0F5RFe9NdAH+hpGhLxsyVtJRF8YG04poFMktGX7PGWuvVuZED3B2Wpsn5qyx4DEu1\npHpJM8Nn+RljlrhfBnaLmC+phEXxgdTjmgUySzZdsySl8EzwgtpR3gveV5dtQwEf4JoFMkumXbMk\npQAAAPAcY0oBAADgOZJSAAAAeI6kFAAAAJ4jKQUAAIDnSEoBAADgOZJSAAAAeI6kFAAAAJ4jKQUA\nAIDnSEoBAADgOZJSAAAAeI6kFAAAAJ4jKQUAAIDnSEoBAADgOZJSAEC3jDFVXsfQnUyIEUDnSEoB\nAF0yxtwqaZLXccSgxhizxBhT4nUgAOJHUgoA6JQxZoakadbahe7rqcaYLcaYBmOMNcZUeBzicdba\nRkkLJC3xOhYA8SMpBQBE5bY4zpc0M1BmrV0paZqklV7F1RU3vhq3dRdABiEpBQB0ZpGkBW4L5HHW\n2hpJ67wJKSZzJM2lGx/ILCSlAIAIbrf8DEkLvY4lXm4SvVLSXK9jARA7klIAQDSzJS0NbyXNIAsk\n3eB1EABiR1IKAIjmBkkrvA7iBKyXVGKMmep1IABiU+B1AACAxAVNRiqTVGOtneOu1znbrVIm6fp4\nWjzdrvsSOYldPLHMlxQ8wWiptXamG88i9zMrJFXKGRowTVKjW15trZ3TxWdXKbI7frG1dmm0+tba\nRmNMjXw8KQtAKJJSAMhsc+VM7CmTtMUYI0l11trZ7nJOS+RMSrorjs+skiRrbXU8gbgJ8a2Slioo\nEbbWVhtjZspJnivkdK3Pt9beJR1PrDe48U4MT6CNMTe4x0xzZ9cHyjcYY9RZYiqpOvBvAeB/dN8D\nQGYrcZO4wEzzqkCyJ2mynJbIzpK2zpS5x8XFbSmdaa2d2c2M/QXByaVbd5qchHVR2GdWyUlI7wpL\nSGfISTiv7iKkevffAiADkJQCQIZyE7bAuM/A2MnjXeDW2jnW2lI3IYxHpZyELp5YFkha0UWrZbCI\neNwYV0qaEbZdaCBJnRd2yEo5KwMs6OI8jXISXQAZgKQUADKUtbY6KAmcHChLwkfHs75niTFmhZyJ\nUfEmv+ECsQdvaRoYShDe8tporZ0d3HoaxRbF928B4CGSUgDIDlOVvAk98XTdL5IzVrRRJ769Z537\nPFE6PuFKSjzZLVECwxAAeIOkFAAyXNBs+WQt4VSn2MdiXu+2Vl4vqcqdlJSocvd5g3S8S19KvAu+\nXHEOQwDgHZJSAMh8gfGkyWwpjbXbu0ZSYAb8SkkLglo44xUYSxr876iWjs/Qj9DNuWgpBTIISSkA\nZL5pUtLGk0puoplAcjnTfe6uG39SeIE7uWmqpIVhE7MCE7c62zK0q4lOFTrxca4A0oSkFAAyXzLH\nkypo8lAsa3we7+Z3JyMtldON31WyODN4hr3bCrpE0kpr7ezgim4scyTdGr47k7sEVacL7stJfjN5\nVyogpxhrrdcxAABOgDGmQc7YznjXI+3qM1coyi5LQbszBZLKRjk7SU10k8bgJLBRzrqlK91jb5Uz\nKWqae3ylW69MznJSC7uIJ3hHp8A40QWdtQ67rbxbJJXGs5sVAO+QlAIAIriL08+31lZ2Wzn2zwwk\npROTONSgq3NNs9ZOS+V5ACQP3fcAgAhuq2tZeJd5BpkrJwEGkCFISgEAnbleXY/Z9CW3lXd9Nwvr\nA/AZklIAQFSBMapukpdMKduP3p00NV/S7O7qAvAXxpQCADrlJnmr5ExYSmh5pSiToyRn/dEpyZ6E\n5E7QWpDMSV8A0oOkFADQJTcxnRs+E99v3BbdRrrtgcxEUgoAAADPMaYUAAAAniMpBQAAgOdISgEA\nAOA5klIAAAB4jqQUAAAAniMpBQAAgOdISgEAAOA5klIAAAB4jqQUAAAAniMpBQAAgOdISgEAAOA5\nklIAAAB4jqQUAAAAniMpBQAAgOdISgEAAOA5klIAAAB47v8DdMSj9CqkwRMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xba5d390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rp = np.linspace(0.01,100.,1000)\n",
    "\n",
    "# rough mass of the Milky Way at r = 100 kpc\n",
    "Mdisk = 8.e11; \n",
    "# particle parameters (arbitrary but realistic)\n",
    "kpc = 3.0856e21\n",
    "rp0 = 100.; v0 = 50.; Lp = v0*rp0\n",
    "\n",
    "# compute Phi eff in (km^2/s^2); 4.30131e-6 = G*Msun/kpc in units of km^2/s^2\n",
    "Phi_eff = 0.5*Lp**2/rp**2 - Mdisk/rp * 4.30131e-6\n",
    "\n",
    "\n",
    "plt.figure(figsize=(3.5,2.5))\n",
    "plt.xlim(1,100.); plt.ylim(-1.e6,1.e6)\n",
    "plt.xscale('log')\n",
    "plt.xlabel(r'$r\\rm\\ (kpc)$'); plt.ylabel(r'$\\Phi_{\\rm eff}\\ \\mathrm{(km^2/s^2)}$')\n",
    "plt.plot(rp, Phi_eff, c='darkslateblue', label=r'$\\Phi_{\\rm eff}$')\n",
    "plt.plot(rp, 0.5*Lp**2/rp**2, '--', c='m', label=r'$L^2/(2r^2)$')\n",
    "plt.plot(rp, -Mdisk/rp * 4.30131e-6, '--', c='b', label=r'$-GM/r$')\n",
    "plt.legend(loc='upper right', frameon=False, fontsize=7)\n",
    "#plt.savefig('fig/phi_eff.png', bbox_inches='tight')\n",
    "plt.grid(linestyle='dotted')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that effective potential looks like gravitational potential at large $r$, but increases steeply at small $r$. This steep increase is kind of a \"wall\" for particle set by the centrifugal force. For a given angular momentum $L$, there is a radius where centrifugal potential $L^2/(2r^2)$ matches gravitational potential $\\Phi$. At smaller $r$ the centrifugal potential dominates and prevents motion of particles to smaller $r$. For a particle with $L>0$ there will thus be the smallest $r$ it can reach within a given gravitational potential. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Solution of the equation of motion in the central field\n",
    "\n",
    "The first Lagrange equation can be rewritten using conservation of $L=r^2\\dot{\\psi}=\\rm const$ to change variable from $dt$ to $d\\psi=L dt/r^2$. This substitution and in the first Lagrange equation gives:\n",
    "$$\\frac{L^2}{r^2}\\frac{d}{d\\psi}\\left(\\frac{1}{r^2}\\frac{dr}{d\\psi}\\right)-\\frac{L^2}{r^3}=-\\frac{\\partial\\Phi}{\\partial r}.$$\n",
    "\n",
    "Further substituting $r$ for $u\\equiv 1/r$, gives:\n",
    "\n",
    "$$\\frac{d^2 u}{d\\psi^2} + u=-\\frac{1}{L^2}\\frac{\\partial\\Phi}{\\partial u}.$$\n",
    "\n",
    "General solution of this equation for motion around a point mass, $\\Phi=-GM/r=-GMu$, has the following form:\n",
    "$$u(\\psi)=C\\cos(\\psi-\\psi_0)+\\frac{GM}{L^2},$$\n",
    "where $C>0$ and $\\psi_0$ are arbitrary constants. $C$ is related to orbit *eccentricity* $e$ and *semi-major axis* $a$:\n",
    "\n",
    "$$e=\\frac{CL^2}{GM}\\ \\ \\ \\  a\\equiv\\frac{L^2}{GM(1-e^2)}.$$\n",
    "\n",
    "The solution can then be rewritten for $r(\\psi)=1/u$ as\n",
    "\n",
    "$$r=\\frac{a(1-e^2)}{1+e\\cos(\\psi-\\psi_0)}.$$\n",
    "\n",
    "Note that if $e<1$, $r$ oscillates within a fixed range as $\\cos(\\psi-\\psi_0)$ changes between $-1$ and $1$ for varying $\\psi$. However, when $e\\geq 1$, $r$ can go to infinity. Thus, $e\\geq 1$ corresponds to *unbound* orbits, while $e<1$ orbits are *bound*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Orbit pericenter, apocenter, r as a function of orbit anomaly\n",
    "\n",
    "Bound orbits have a fixed range within which $r$ varies, *apocenter*:\n",
    "$$r_{\\rm apo}=\\frac{a(1-e^2)}{1-e}=a(1+e)$$\n",
    "and *pericenter:\n",
    "$$r_{\\rm peri}=\\frac{a(1-e^2)}{1+e}=a(1-e).$$\n",
    "\n",
    "If you think about these expressions a bit, it becomes clear that solution of the equation of motion can also be written as\n",
    "<br>\n",
    "\n",
    "$$r=a(1-e\\cos\\eta),$$\n",
    "\n",
    "where $\\eta$ is called the *eccentric anomaly* to distinguish it from the *true anomaly* $\\psi-\\psi_0$. The two are related, of course:\n",
    "$$\\sqrt{1-e}\\,\\tan\\frac{1}{2}(\\psi-\\psi_0)=\\sqrt{1+e}\\,\\tan\\frac{1}{2}\\eta.$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Orbit period\n",
    "\n",
    "Note that $t-t_0=\\int_{t_0}^tdt$ and $\\dot{\\psi}dt=d\\psi$. Assuming $t_0$ corresponds to the pericenter passage ($\\eta=0$):\n",
    "\n",
    "$$t=\\int_{\\psi_0}^{\\psi}\\frac{d\\psi}{\\dot{\\psi}}=\\int\\frac{r^2}{L}d\\psi=\\frac{a^2}{L}\\int\\limits^\\eta_0d\\eta\\frac{d\\psi}{d\\eta}(1-e\\cos\\eta)^2.$$\n",
    "\n",
    "$d\\psi/d\\eta$ can be obtained from the relation of $\\psi$ and $\\eta$ above and the integral gives: \n",
    "$$t=\\frac{a^2}{L}\\sqrt{1-e^2}(\\eta - e\\sin\\eta).$$\n",
    "$\\eta=\\pi$ is the time of the apocenter passage, so $t$ for $\\eta=\\pi$ corresponds the time between pericenter and apocenter. Orbit period will be twice that:\n",
    "\n",
    "$$T_r=2\\pi\\frac{a^2}{L}\\sqrt{1-e^2}=2\\pi\\sqrt{\\frac{a^3}{GM}},$$\n",
    "where expression for $a$ above was used to get the second equality - *the Kepler's 3rd law*.\n",
    "\n",
    "Note that from the equations above, it is clear that $\\eta$ and $\\psi-\\psi_0$ change by $2\\pi$ in one $T_r$ for an orbit around a point mass. \n",
    "\n",
    "Note that when when two bodies orbit each other, the above equations stand. The two bodies can be viewed as orbiting a point mass at their common center of mass. The mass $M$ in the equation above in this case represents the sum of the masses of the bodies, $M=M_1+M_2$. "
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}