{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# High Order Taylor Maps I\n", "(by Dario Izzo)\n", "\n", "In this notebook we consider the system of differential equations $\\dot{\\mathbf y} = \\mathbf f(\\mathbf y)$:\n", "\n", "$$\n", "\\begin{array}{l}\n", "\\dot r = v_r \\\\\n", "\\dot v_r = - \\frac 1{r^2} + r v_\\theta^2\\\\\n", "\\dot \\theta = v_\\theta \\\\\n", "\\dot v_\\theta = -2 \\frac{v_\\theta v_r}{r}\n", "\\end{array}\n", "$$\n", "\n", "which describe, in non dimensional units, the Keplerian motion of a mass point object around some primary body. \n", "We show how we can build a high order Taylor map (HOTM, indicated with $\\mathcal M$) representing the final state of the system at the time $T$ as a function of the initial conditions. \n", "\n", "In other words, we build a polinomial representation of the relation $\\mathbf y(T) = \\mathbf f(\\mathbf y(0), T)$. Writing the initial conditions as $\\mathbf y(0) = \\overline {\\mathbf y}(0) + \\mathbf {dy}$, our HOTM will be written as:\n", "\n", "$$\n", "\\mathbf y(T) = \\mathcal M(\\mathbf {dy})\n", "$$\n", "\n", "and will be valid in a neighbourhood of $\\overline {\\mathbf y}(0)$." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# We use numpy for arrays and pyaudi for the differential algebra machinery\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline \n", "from pyaudi import gdual_double as gdual\n", "from pyaudi import sin, cos" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# This is a simple Runga Kutta fourth order numerical integrator with fixed step.\n", "# It is programmed to work both with floats and gduals. It infers the type from the initial conditions\n", "def rk4(f, t0, y0, h, N):\n", " t = t0 + np.arange(N+1)*h\n", " y = np.array([[type(y0[0])] * np.size(y0)] * (N+1))\n", " y[0] = y0\n", " for n in range(N):\n", " xi1 = y[n]\n", " f1 = f(t[n], xi1)\n", " xi2 = y[n] + (h/2.)*f1\n", " f2 = f(t[n+1], xi2)\n", " xi3 = y[n] + (h/2.)*f2\n", " f3 = f(t[n+1], xi3)\n", " xi4 = y[n] + h*f3\n", " f4 = f(t[n+1], xi4)\n", " y[n+1] = y[n] + (h/6.)*(f1 + 2*f2 + 2*f3 + f4)\n", " return y" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# The Equations of Motion of Keplerian motion in non dimensional spherical coordinates.\n", "def eom_kep_polar(t,y):\n", " return np.array([y[1], - 1 / y[0] / y[0] + y[0] * y[3]*y[3], y[3], -2*y[3]*y[1]/y[0]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### We perform the numerical integration using floats (the standard way)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# Fixed step size\n", "step = 0.01\n", "# Number of steps\n", "n_steps = 500\n", "# The initial conditions\n", "ic = [1.,0.1,0.,1.]\n", "# The intitial time (irrelevant as the system is autonomous)\n", "it = 0.\n", "y = rk4(eom_kep_polar, it, ic, step, n_steps)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0,0.5,'y')" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEWCAYAAABBvWFzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl4VOXZ+PHvnZ2E7AkhQEJANkHZ\njICodUNFq2Ld96Vaa1urfbu8au2v9fW1rbWL1dalaq1rRaVasWoRVKwoW1B2hIQ9JEBCQkISst+/\nP87Bd8AJmZCZOVnuz3XNNWfOec4595w5M/c8z1keUVWMMcaYzorwOgBjjDE9gyUUY4wxQWEJxRhj\nTFBYQjHGGBMUllCMMcYEhSUUY4wxQWEJxUMicq+IvHiY6WtE5NQwhhSQAOK+WkTeC9G6R4rI5yKy\nT0RuD8U6OkpEVESGhXgdz4rI/aFcR2eIyE9F5OlglO3I/tPevhgOvp+NiJwsIusPUzZXRGpEJDJ8\nEYaPJZQgEZEbRGSViNSJyE4ReVxEUjqzTFUdo6rz3eV7/sXxR0Ty3B/UqAPjVPUlVT0rRKv8b2C+\nqiaq6iMhWkebRGS+iNwc7vV2dar6K1UNaLv4lvVg/wkpVf1YVUceeC0iW0Rkms/0baraV1VbvIkw\ntCyhBIGI/Aj4DfATIBmYAgwG5opITBvzRPkbH8SYQrp8Dw0G1rQ1saf+8zOmO7CE0kkikgT8D/B9\nVf23qjap6hbgMpwfv2vccveKyCwReVFEqoEb3EXEicgrbhPOZyIyzmfZW0RkmohMB34KXO5Wl1e0\nEcsWEblTRFYCtSISJSIDROQfIlImIpt9m4lEZJKIFIhItYjsEpE/uONPFZFiP8uexlf9x33e68Z2\ngltbW+Azr4rIrSJSKCKVIvKoiIg7LVJEfi8i5W58tx36j9VnOR8ApwF/dtc1wm1ueFxE3hGRWuA0\nEUkWkefd97xVRH4mIhHuMm4QkU9E5CER2Ssim0Rkqjt+u4jsFpHr29i+vwRO9ln/n30mT/P3/tz5\nviki69xpc0RksL/lu2VPEpFP3di2i8gNPpNTReRtd19ZLCJH+cw3VUSWikiV+zzVZ9oN7vvc527j\nqwOJ7XCfm5+4v6xBy//VOq4XkW3uZ3uPv7IEtv887G6LahFZJiInt7X9/MQ1Q0SWu/NuFOe7hPu9\nmC0iFSJSJCLfOiS+V919aJ84Tc/5PtMniPNd3ScirwBxPtO+/O6IyAtALvCW+97+Ww6pkXUyjjtF\nZIc7bb2InBHodgkZVbVHJx7AdKAZiPIz7TngZXf4XqAJuBAnkffxGXcJEA38GNgMRLvzbAGm+cz/\nYjuxbAGWAznu8iOAZcDPgRhgKLAJONstvxC41h3uC0xxh08Fiv0s+yuxAHmA+r5/nGS5wOe1Av8C\nUnC+YGXAdHfarcBaYBCQCsw7dHmHxDEfuNnn9bNAFXCi+37jgOeBN4FEN74NwE0+sTUDNwKRwP3A\nNuBRIBY4C9gH9A1k/QG8vwuBIuBoIAr4GfBpG8vOddd9pbs/pAPjfd5nBTDJXc5LwEx3WhpQCVzr\nTrvSfZ0OJADVwEi3bDYwJpDYDve+/MTub594Cmc/HAc0AEcf4f5zjfteooAfATuBuPa+F+62qgLO\ndPeNgcAod9pHwGM4+8t4972d4bPMeuBcdx/5NbDInRYDbAX+y/2MLsH5Dt/v77uDz/fG3/vtRBwj\nge3AAJ/lHuX176HVUDovAyhX1WY/00rd6QcsVNV/qmqrqu53xy1T1Vmq2gT8AWfHmtKJeB5R1e3u\n8o8HMlX1PlVtVNVNOF/yK9yyTcAwEclQ1RpVXdSJ9bbnAVXdq6rbgA9xvjzg1OQeVtViVa0EHjiC\nZb+pqp+oaivOe7ocuFtV96lTW/w9zo/tAZtV9W/qtGO/gpOA71PVBlV9D2gEOnqQva33923g16q6\nzt1HfgWMb6OWcjUwT1VfVqemu0dVl/tMf11Vl7jLeclnHV8HClX1BVVtVtWXgS+A893prcAxItJH\nVUtV9UCTYSCxtfW+AvE/qrpfVVcAK3ASS4ep6ovutmhW1d/jJP6R7c0H3AQ8o6pz3e/cDlX9QkRy\ngJOAO1W13t3GT3PwPrJAVd9x95EXfGKfgpNI/uh+RrOApUfyvjoZRwvOdhgtItGqukVVNx5JHMFk\nCaXzyoEM8X/MItudfsB2P2W+HOf+IBYDAzoRj+86BgMD3OaTvSKyF6fpLMudfhMwAvjCbSY5rxPr\nbc9On+E6nBoROO/VN2Z/26g9vvNk8H//Ig/YivPv9IBdPsP7AVT10HF96Zi23t9g4GGf7V8ByCHx\nHJADHO5H4XDbcOshZbcCA1W1FifB3gqUuk1mozoQW1vrDERn5v2SiPzIbZarcuNM5uA/am1pa3sO\nACpUdZ/PuEP3kUNjj3O/4wOAHepWC3zmPRJHHIeqFgE/wKnF7BaRmSLSmd+NoLCE0nkLcarzF/mO\nFJEE4BzgfZ/R/m7tnOMzTwRO00+Jn3KB3hbat9x2nH/jKT6PRFU9F0BVC1X1SqAfzkkFs9y4a4F4\nn7gigcwA1nckSnHe8wE5bRU8DN8YynFqKb7/snOBHUew3PbWFYjtwLcP+Qz6qOqnbZQ9ys/49pRw\n8PsFn/esqnNU9UycPzhf4NRSOxpbqBx2e7rHS+7EqcmmqmoKTjOW32M5h2hre5YAaSKS6DMu0H2k\nFBh4yLGk3MOUP9z760wcqOrfVfUknM9ecb7DnrKE0kmqWoVzUP5PIjJdRKJFJA94Dae28UI7izhO\nRC5y//38ACc5+Wt62gXkuUknUEuAavfgXR9xDoAfIyLHA4jINSKS6daM9rrztOAcc4gTka+LSDRO\n23psG+sow2lSGdqBuHy9CtwhIgPFOc36ziNcDgBu08CrwC9FJNFtvvkhEKxTrnfRsff6BHC3iIwB\nEOeEgUvbKPsSzsH9y8Q5oSJdRAJpYnoHGCEiV7nzXQ6MBv4lIlkicoH7R6EBqMH5jDsaW6i0t/8k\n4hzzKgOiROTnQFKAy/4rcKOInCEiEe4+NkpVtwOfAr8WkTgRGYtTW38pgGUudOO53d3WF+Ecq2lL\nm/tLZ+IQ53qs00UkFuc4y37+73P1jCWUIFDVB3Gakn6HcwB0Mc6/ozNUtaGd2d/EaZI4cFD1Ivd4\nyqFec5/3iMhnAcbVgtOOPh7nYH85ThttsltkOrBGRGqAh4Er3LbcKuC7btkdODWWYvxQ1Trgl8An\nbtNJR4//PAW8B6wEPsf5cWymc1+O7+PEvAlYAPwdeKYTy/P1MHCJOGc9tXsdjKq+gfPPcaY4Z/et\nxqm5+iu7DecA7I9wmp+WE8BxB1XdA5znzrcH51qd81S1HOc7/iOcf8MVwCk4n22HYguVAPafOcC7\nOH9ytuL8eAbULKqqS3BOvngIp1bzEf9Xk7sS50B2CfAG8AtVnRvAMhtxWiNuwPnOXg68fphZfg38\nzH1vP/Yz/YjiwPmD9wDOd3onTivDTwOYL6Tk4KZAY7wlIucAT6hqm6fWGmO6JquhGE+5TXHnus0H\nA4Ff4PxTM8Z0M1ZDMZ4SkXicpohROO3AbwN3qGq1p4EZYzrMEooxxpigsCYvY4wxQdFTbyDoV0ZG\nhubl5XkdhjHGdCvLli0rV9W2rkX7Uq9KKHl5eRQUFHgdhjHGdCsiEtDdAKzJyxhjTFBYQjHGGBMU\nllCMMcYEhSUUY4wxQWEJxRhjTFBYQjHGGBMUniYUEXlGnD68V7cxXUTkEbev5ZUiMtFn2vXi9HVd\nKG30AW6MMSZ8vL4O5Vngzzh9gPtzDjDcfUwGHgcmi0gazk0E83E6llkmIrPdLmSNCTtVpb6pldrG\nZmobmqltaKG2sZmahmbq3OHmFqWltZWmFqWlVWlqbaWlRWluVaeT8QghMkIOfo6MIDpC6BMTSUJM\nFPGxznNCbCTxMVEkxETRNy6KyIhA+psyJrQ8TSiq+h+3M6q2zACed7vbXCQiKSKSDZwKzFXVCgAR\nmYvTt8fLoY3YGP8uf3IRSzZXeLJuEUjuE01afAxpCTGkJsQ4w31jSE+IoX9yHNnJfchOjqNfYixR\nkdbSbULD6xpKewZycGc6xe64tsZ/hYjcAtwCkJt7uJ46jTlyl+fncNrIfvQ9UHOIdWoRCbFOLSI+\nJpLoyAiiIp2aR1RkxEE1EYBWhebWVlpanVrLgdpLc2srdY0tX9Z06hrdGlCDUwOqrm+msraRirpG\nKmsb2V5Rx4rte6msa6Sp5eCbv0YI9EuMIzsljuzkOHLTEhiSEU9eegJDMhLITIzl4N5tjQlcV08o\n/vZsPcz4r45UfRJ4EiA/P99urWxC4uLjBnV6GZECkRGRQYjGoapU72+mtHo/pXvrKa2qZ2fVfkqq\n6tlZVc8XpfuYu3bXQUknISaSwW5yGdavL0dnJzKyfxK5afHWrGba1dUTSjGQ4/N6EE5XmcU4zV6+\n4+eHLSpjugERITk+muT4aEb1998Ne3NLKyV769m8p5Yt5bVsLq9l655a1pRU8c7qUg70bhEXHcHI\nrERG9k9kVP8kjhmYzDEDk4iP6eo/ISacuvreMBu4TURm4hyUr1LVUhGZA/xKRFLdcmcBd3sVpDHd\nVVRkBLnp8eSmx3PKiINvJru/sYXC3fv4onQfX+zcx/pd1by/bjevFhQDTvPZiKxExuekMC4nhXGD\nUhiR1deO0fRiniYUEXkZp6aRISLFOGduRQOo6hPAO8C5QBFQB9zoTqsQkf8FlrqLuu/AAXpjTHD0\niYlk7KAUxg5KOWj87n31rCquYsX2vSwvruLd1TuZudQ5pBkXHcGEnFQmD01j8pB0JuSmEBcdvGY8\n07X1qh4b8/Pz1W5fb0xwqSpb99Sxongvn2/by5LNFazbWY0qxERGMD4nhclD05gyNJ3jBqdagumG\nRGSZqua3W84SijEm2Krqmli6pYLFm/eweHMFq3dU0apODeaEoemcMiKTU0b2Y0hGgtehmgBYQvHD\nEoox3thX7ySY/2wo56MNZWwurwVgsHvs5rSR/Zg6LJ3YKKu9dEWWUPywhGJM17B1Ty0fbSjjo/Vl\nfLpxD/ubWugbG8Vpo/oxfUx/Th2ZSUJsVz9nqPewhOKHJRRjup6G5hY+3biHOat3MnftLvbUNhIT\nFcHXhmdw9pj+nDWmP8l9or0Os1ezhOKHJRRjuraWVqVgSwX/XrOT99bsYsfe/cRERXDGqH5cOGEg\np47MtGYxD1hC8cMSijHdh6qyoriKf36+g3+tLKG8ppHkPtGce2w2F44fwPF5aUTY1fthYQnFD0so\nxnRPzS2tLCgq583lJfx79U72N7UwOD2eK47P5ZLjBpGZGOt1iD2aJRQ/LKEY0/3VNTYzZ81OXl6y\nnSWbK4iKEM4ak8WVk3I58agMq7WEgCUUPyyhGNOzFO2uYeaSbfzjs2Iq65rISevDNZMHc8WkXDuQ\nH0SWUPywhGJMz1Tf1MKcNTt5afE2lmyuID4mkkuPG8SNJw4hzy6e7DRLKH5YQjGm51u9o4pnPtnM\nWytKaG5VzhiVxU0nDWHK0DTr6+UIWULxwxKKMb3H7up6Xli0lRcXbaWyrolxOSncfvowTh/VzxJL\nB1lC8cMSijG9T31TC7OWFfPERxsprtzP6Owkvn/6MM4e098O4AfIEoofllCM6b2aWlr55+c7eGz+\nRjaX1zIiqy+3nT6c847NtsTSjkATivWEY4zpFaIjI7g0P4d5PzyFh68Yjyrc/vLnfP1PC/hw/W56\n05/rUPE0oYjIdBFZLyJFInKXn+kPichy97FBRPb6TGvxmTY7vJEbY7qryAhhxviBzPnB13j4ivHU\nNjRz49+WcsWTi/hsW6XX4XVrnjV5iUgksAE4E6eP+KXAlaq6to3y3wcmqOo33dc1qtq3I+u0Ji9j\nzKEam1t5eck2/vRBIeU1jZw9Jos7p49iaGaHfl56tO7Q5DUJKFLVTaraCMwEZhym/JXAy2GJzBjT\na8RERXD91Dw++slp/Ne0ESwoLOfsP/6HX7+7jpqGZq/D61a8TCgDge0+r4vdcV8hIoOBIcAHPqPj\nRKRARBaJyIVtrUREbnHLFZSVlQUjbmNMD5QQG8Ud04Yz/yenMWP8QP7y0SZO/918/vn5Dju+EiAv\nE4q/0yra+tSuAGapaovPuFy3CnYV8EcROcrfjKr6pKrmq2p+ZmZm5yI2xvR4mYmx/O7Scbzx3alk\nJ8fxg1eWc+kTC1lTUuV1aF2elwmlGMjxeT0IKGmj7BUc0tylqiXu8yZgPjAh+CEaY3qrCbmpvPHd\nE/nNxceyubyWC/78Cb/59xfUN7W0P3Mv5WVCWQoMF5EhIhKDkzS+craWiIwEUoGFPuNSRSTWHc4A\nTgT8Hsw3xpgjFREhXH58Lh/86FQunjiQx+dv5JyHP2bxpj1eh9YleZZQVLUZuA2YA6wDXlXVNSJy\nn4hc4FP0SmCmHtyIeTRQICIrgA+BB9o6O8wYYzorOT6aBy8Zx4s3Taa5tZXLn1zEPW+sYl99k9eh\ndSl2pbwxxnRAXWMzf3hvA898spns5D784bJxTB6a7nVYIdUdThs2xphuJz4mip+dN5pZ35lKVKRw\nxVOLePDfX9DY3Op1aJ6zhGKMMUdgYm4qb99+Mpcdl8Nj8zdy8eOfsrGsxuuwPGUJxRhjjlDf2Ch+\nc8lYnrhmItsr6zjvkQW8VrC9/Rl7KEsoxhjTSdOPyWbOD77G+JwUfjJrJXe/vrJXnl5sCcUYY4Ig\nKymOF26axHdOPYqXl2zn0icWUlxZ53VYYWUJxRhjgiQqMoI7p4/iyWuPY0t5Lef9aQEfbeg9t3yy\nhGKMMUF21pj+zP7+SfRPiuPGvy3h2U82ex1SWFhCMcaYEBiSkcDr353K6aOyuPettfzizdU0t/Ts\nU4stoRhjTIjEx0Txl2uP4+aThvDcwq3c/HxBj7663hKKMcaEUGSE8LPzRnP/hcfwcWE5lz6xkF3V\n9V6HFRKWUIwxJgyumTKYZ244nm0VdVz6xEK2V/S8M8AsoRhjTJicMiKTl26eTNX+Ji554lMKd+3z\nOqSgsoRijDFhNCE3lVe+PYVWhcv+spBVxT2n4y5LKMYYE2aj+ifx2rdPICE2iqueWsSK7Xu9Diko\nLKEYY4wH8jISeO3WE0hNiOG6Z5awtqTa65A6zRKKMcZ4JDu5Dy/dPJmEmEiu+evibn9MxdOEIiLT\nRWS9iBSJyF1+pt8gImUistx93Owz7XoRKXQf14c3cmOMCY6ctHhe+tYUoiKEq55ezObyWq9DOmKe\nJRQRiQQeBc4BRgNXishoP0VfUdXx7uNpd9404BfAZGAS8AsRSQ1T6MYYE1RDMhJ46ebJtLQq1/51\nMWX7GrwO6Yh4WUOZBBSp6iZVbQRmAjMCnPdsYK6qVqhqJTAXmB6iOI0xJuSGZyXy7I3Hs6emkZue\nW0pdY7PXIXWYlwllIODbE02xO+5QF4vIShGZJSI5HZwXEblFRApEpKCsrPfc9dMY0/2MHZTCn66c\nwOodVdz+8ue0tKrXIXWIlwlF/Iw7dOu9BeSp6lhgHvBcB+Z1Rqo+qar5qpqfmZl5xMEaY0w4TBud\nxb0XjGHeut38z1trUO0+ScXLhFIM5Pi8HgSU+BZQ1T2qeqAx8SnguEDnNcaY7uq6E/K4+aQhPL9w\nKzOXdp8uhb1MKEuB4SIyRERigCuA2b4FRCTb5+UFwDp3eA5wloikugfjz3LHGWNMj3D3uUdz8vAM\nfvHmmm5z4aNnCUVVm4HbcBLBOuBVVV0jIveJyAVusdtFZI2IrABuB25w560A/hcnKS0F7nPHGWNM\njxAZITxyxQQyE2P57kufUVHb6HVI7ZLu1D7XWfn5+VpQUOB1GMYYE7CVxXu55ImFTB6SxrM3TiIy\nwt8h5NASkWWqmt9eObtS3hhjurCxg1K474IxfFxYzl/+s9HrcA7LEooxxnRxlx+fw9ePzeahuRtY\nU9J1705sCcUYY7o4EeH+C48hNT6G/3plOfVNLV6H5JclFGOM6QZSE2L47aXj2LCrht/NWe91OH5Z\nQjHGmG7ilBGZXDtlME8v2MyyrZVeh/MVllCMMaYbueucUWQnx3HPG6toamn1OpyDWEIxxphuJCE2\ninsvGMMXO/fxzILNXodzEEsoxhjTzZw9pj/Tjs7ij/MKKa3a73U4X7KEYowx3dAvzh9Niyq/7UIH\n6C2hGGNMN5STFs+NJ+bx+mc7WL2ja1ybYgnFGGO6qe+dNozU+Gh++fa6LnGbe0soxhjTTSXFRXP7\nGcNZuGkPCzft8TocSyjGGNOdXTkpl6ykWB6eV+h1KJZQjDGmO4uLjuTWU45i8eYKFnlcS7GEYowx\n3dyVk3LJTIzl0Q+LPI3D04QiItNFZL2IFInIXX6m/1BE1orIShF5X0QG+0xrEZHl7mP2ofMaY0xv\nERcdyfUnDObjwnIKd+3zLA7PEoqIRAKPAucAo4ErRWT0IcU+B/JVdSwwC3jQZ9p+VR3vPi7AGGN6\nsSsn5RIbFcHfPt3iWQxe1lAmAUWquklVG4GZwAzfAqr6oarWuS8XAYPCHKMxxnQL6X1juXD8QF7/\nrJiquiZPYvAyoQwEtvu8LnbHteUm4F2f13EiUiAii0TkwrZmEpFb3HIFZWVlnYvYGGO6sGtPGEx9\nUyuzV5Z4sn4vE4q/jpH9XpkjItcA+cBvfUbnun0cXwX8UUSO8jevqj6pqvmqmp+ZmdnZmI0xpssa\nMyCJUf0TmbWs2JP1e5lQioEcn9eDgK+kVRGZBtwDXKCqDQfGq2qJ+7wJmA9MCGWwxhjT1YkIl+bn\nsGL7Xk8OznuZUJYCw0VkiIjEAFcAB52tJSITgL/gJJPdPuNTRSTWHc4ATgTWhi1yY4zpomaMH0Bk\nhPDm8vA3e3mWUFS1GbgNmAOsA15V1TUicp+IHDhr67dAX+C1Q04PPhooEJEVwIfAA6pqCcUY0+tl\n9I1lUl4ac9bsDPu6o8K+Rh+q+g7wziHjfu4zPK2N+T4Fjg1tdMYY0z1NP6Y/v5i9ho1lNRyV2Tds\n67Ur5Y0xpoc5a0wWQNhrKZZQjDGmh8lO7sOo/oksKCwP63otoRhjTA904rAMCrZWUt/UErZ1WkIx\nxpge6MRh6TQ2t7Jsa2XY1mkJxRhjeqBJQ9KJEFiyuSJs67SEYowxPVDf2CiG9evLqjD2N28JxRhj\neqixg1JYWbw3bP3NW0IxxpgeauygZMprGimtqg/L+iyhGGNMDzUyKxGAot01YVmfJRRjjOmhhmQm\nALC5vDYs67OEYowxPVRm31j6xkZZQjHGGNM5IkJuWjxb91hCMcYY00lZSbGU1TS0XzAILKEYY0wP\nlpkYS9k+SyjGGGM6KTMxlvKaRlpbQ38tiiUUY4zpwZLiomlpVeqbQ3+TSE8TiohMF5H1IlIkInf5\nmR4rIq+40xeLSJ7PtLvd8etF5OxQxllR28ju6vBcGGSMMcEUFx0JwP7GHpxQRCQSeBQ4BxgNXCki\now8pdhNQqarDgIeA37jzjsbpg34MMB14zF1eSPz3rBXc8LeloVq8McaETFy08zNf39wa8nW1m1BE\n5DYRSQ3BuicBRaq6SVUbgZnAjEPKzACec4dnAWeIiLjjZ6pqg6puBorc5YVEhAitYboXjjHGBFNs\nlPNfuyEM/aIEUkPpDywVkVfdJioJ0roHAtt9Xhe74/yWUdVmoApID3BeAETkFhEpEJGCsrKyIwrU\nEooxprtqcQ/GR0WEvkGq3TWo6s+A4cBfgRuAQhH5lYgc1cl1+0tMh/5qt1UmkHmdkapPqmq+quZn\nZmZ2MERHZIR8+aEYY0x38mVCiQxWXaBtAaUsde59vNN9NAOpwCwRebAT6y4GcnxeDwJK2iojIlFA\nMlAR4LxBExEhWD4xxnRHTa3OsZOoiC6QUETkdhFZBjwIfAIcq6rfAY4DLu7EupcCw0VkiIjE4Bxk\nn31ImdnA9e7wJcAHbnKbDVzhngU2BKcGtaQTsRxWQkwktQ3NoVq8McaETF2Dc+ykT0zIzlv6UlQA\nZTKAi1R1q+9IVW0VkfOOdMWq2iwitwFzgEjgGVVdIyL3AQWqOhunme0FESnCqZlc4c67RkReBdbi\n1Ji+p6ohO+KUmhBDZV0jqkrwDiEZY0zoVdY1EhUh9I0N5Oe+c9pdg6r+/DDT1nVm5ar6DvBOW+tT\n1Xrg0jbm/SXwy86sP1Bp8TE0tSg1Dc0kxkWHY5XGGBMUlXVNpMRHh+XPsF0pH4C0hBgAKmubPI7E\nGGM6prK2kdT4mLCsyxJKADISYwHYaVfLG2O6meK9dQxI6ROWdVlCCcCQdKfXsy1h6qTGGGOCZXvF\nfnLSLKF0GQNS4oiOFDaHqZMaY4wJhur6Jqr2N5GTGh+W9VlCCUBUZAQ5afFsLrOEYozpPop21wCQ\nl5EQlvVZQgnQ0IwENpbVeB2GMcYEbF1pNQCjs5PCsj5LKAE6ZmAyRWU11NgFjsaYbmJdaTWJsVEM\nSrVjKF3K+JwUVGFl8V6vQzHGmICsKanm6OyksF2QbQklQONzUgD4fJslFGNM17e/sYVVxVVMHByK\n3kf8s4QSoJT4GIZkJFhCMcZ0C59vr6S5VZk8JC1s67SE0gGTh6SxeNMemlpC3/OZMcZ0xpLNFYhg\nNZSu6pQRmexraLZaijGmy/u4sJxjBiST3Cd89x+0hNIBJw7PICpCmL9+t9ehGGNMmypqG/lsWyVn\nHN0vrOu1hNIBSXHRTBycyvz1R9aVsDHGhMP89btRhTNGZYV1vZZQOmja0f1YW1rNVrsNizGmi3pv\nzS76JcYyZkB4Lmg8wBJKB503dgAi8ObykPU4bIwxR6xqfxMffLGbc4/NJiIM3f768iShiEiaiMwV\nkUL3+SunIYjIeBFZKCJrRGSliFzuM+1ZEdksIsvdx/hwxT4gpQ+T8tL45/IdOL0RG2NM1/Hv1aU0\ntrRy4YSBYV+3VzWUu4D3VXU48L77+lB1wHWqOgaYDvxRRFJ8pv9EVce7j+WhD/n/zBg/kE1ltawp\nqQ7nao0xpl3//LyEvPR4xg1KDvu6vUooM4Dn3OHngAsPLaCqG1S10B0uAXYDmWGL8DDOPbY/MVER\nvLJ0u9ehGGPMlzaV1bBw0x47a8JAAAAdBklEQVQunjgobLdb8eVVQslS1VIA9/mw57aJyCQgBtjo\nM/qXblPYQyISe5h5bxGRAhEpKCsLztlZKfExnD92AP/4rJjqeusW2BjTNby0eBtREcLlk3I8WX/I\nEoqIzBOR1X4eMzq4nGzgBeBGVT1wifrdwCjgeCANuLOt+VX1SVXNV9X8zMzgVXBumJpHXWMLswqK\ng7ZMY4w5UvsbW5i1rJjpx/SnX2KcJzFEhWrBqjqtrWkisktEslW11E0Yfq8UFJEk4G3gZ6q6yGfZ\npe5gg4j8DfhxEEMPyLGDkpmYm8LzC7dww9S8sJ9NYYwxvl7/vJiq/U1cO2WwZzF41eQ1G7jeHb4e\nePPQAiISA7wBPK+qrx0yLdt9FpzjL6tDGm0brp+ax5Y9dcxdt8uL1RtjDADNLa385aNNjMtJYVIY\nbwZ5KK8SygPAmSJSCJzpvkZE8kXkabfMZcDXgBv8nB78koisAlYBGcD94Q3f8fVjsxmcHs8j7xfa\nKcTGGM+8vaqUbRV1fPfUozw5GH9AyJq8DkdV9wBn+BlfANzsDr8IvNjG/KeHNMAARUVG8P3Th/Pj\n11Ywb91uzhwd3tscGGNMa6vy+PyNDO/XlzOP9vY3yK6U76QLxw9gcHo8f5y3wWopxpiwm72ihC92\n7uO204d5fizXEkonRUVGcNtpw1hTUs27q3d6HY4xphdpbG7l93PXM2ZAEuePHeB1OJZQguEbEwYy\nMiuRX72zjvqmFq/DMcb0En9fvJXtFfv57+mjPK+dgCWUoIiKjODn54+muHI/f12w2etwjDG9QEVt\nIw/NK2TqUel8bXiG1+EAllCC5sRhGZw5OotHPyxid3W91+EYY3q43875gtqGZu69YIynZ3b5soQS\nRPecezRNLa388p11XodijOnBlm/fy8yl27lhah4jshK9DudLllCCKC8jge+cOow3l5fwwRd2saMx\nJviaWlr56euryOgbyx3ThnsdzkEsoQTZbacNY0RWX376+mq7caQxJugen7+RtaXV3H/hMSTGRXsd\nzkEsoQRZTFQED14yjt376nng3S+8DscY04N8sbOaP31QyPnjBnD2mP5eh/MVllBCYHxOCjedNIS/\nL97Gh+v93vfSGGM6pKG5hR+9uoKkuGj+54IxXofjlyWUEPnRWSMZ1T+RH7+6ws76MsZ02m//vZ41\nJdX8+qJjSUuI8TocvyyhhEhcdCR/vmoCtY3N/Nery2lttduyGGOOzIfrd/P0gs1cd8JgzuqCTV0H\nWEIJoWH9Ern3/DF8UrSHxz/a2P4MxhhziJ1V9fz41RWM6p/IT8892utwDssSSohdfnwO543N5g9z\nN7CgsNzrcIwx3Uh9UwvffqGA+qYW/nzVBOKiI70O6bAsoYSYiPDAxWM5KjOB7/39M7buqfU6JGNM\nN6Cq3PPGalYUV/GHy8czrF/XuYCxLZ4kFBFJE5G5IlLoPqe2Ua7Fp3Ot2T7jh4jIYnf+V9zeHbus\nvrFRPH3d8YjAt54voKah2euQjDFd3LOfbuEfnxVzxxnDu+Qpwv54VUO5C3hfVYcD77uv/dmvquPd\nxwU+438DPOTOXwncFNpwOy83PZ5Hr5rIxrJafjBzOS12kN4Y04YFheXc//Y6zhydxR1ndK2r4Q/H\nq4QyA3jOHX4Op1/4gLj9yJ8OzDqS+b104rAMfn7eaOat28X/vLXGOuQyxnzF6h1VfPuFAoZl9uUP\nl43rErelD5RXCSVLVUsB3Od+bZSLE5ECEVkkIgeSRjqwV1UPtBsVAwPbWpGI3OIuo6CsrCxY8R+x\n66fmccvXhvL8wq08+mGR1+EYY7qQbXvquOFvS0nuE81z35zU5W6t0p6Q9SkvIvMAfw1/93RgMbmq\nWiIiQ4EPRGQVUO2nXJt/9VX1SeBJgPz8/C5RJbhr+ijK9jXwu/c20C8xjsuOz/E6JGOMx8prGrju\nmcU0t7Yy85YT6J8c53VIHRayhKKq09qaJiK7RCRbVUtFJBvwe38SVS1xnzeJyHxgAvAPIEVEotxa\nyiCgJOhvIIQiIoQHLxnLntpG7n5jFUl9opl+TPc46GaMCb6q/U1889ml7Kyu56Wbp3SLM7r88arJ\nazZwvTt8PfDmoQVEJFVEYt3hDOBEYK06Bx4+BC453PxdXXRkBI9fPZGxg5K57e+fMWeN9UdvTG9U\ntb+J6/66mHWl1Tx29USOG+z3pNduwauE8gBwpogUAme6rxGRfBF52i1zNFAgIitwEsgDqrrWnXYn\n8EMRKcI5pvLXsEYfJAmxUTz3zUkcMzCZ7730Ge9ZUjGmV6mub+K6Z5awtrSax68+jtNHZXkdUqdI\nbzrTKD8/XwsKCrwO4yuq65u49q9LWFtSxWNXH8eZo7v3TmWMaV93+t6LyDJVzW+vnF0p3wUkxUXz\n/DcnMTo7ie+8uIw3l+/wOiRjTAiV1zRw9VOLWVtSxaNXTezSyaQjLKF0Ecl9onnh5skcNziVO2Yu\n59lPNnsdkjEmBLZX1HHJ459SuHsff7n2uC599+COsoTShSTFOeeenzU6i3vfWssf3ltvFz8a04Os\nK63mosc/pbKuiZduntztj5kcyhJKFxMXHcljV0/k8vwcHvmgiJ++sYqmllavwzLGdNLiTXu47C8L\niRThtVtP4LjBaV6HFHQhuw7FHLmoyAgeuPhYMhNj+fOHRWzdU8djV08kJb5L3wPTGNOGV5Zu42f/\nXE1uWjzP3zSZgSl9vA4pJKyG0kWJCD8+eyR/uGwcBVsq+cZjn7KprMbrsIwxHdDSqvzvv9Zy5z9W\nMWVoOq9/58Qem0zAEkqXd9HEQfz9W5Op3t/EhY9+wseF3t+PzBjTvup65+r3vy7YzA1T8/jbDceT\nHN+97s3VUZZQuoH8vDT++b0TyU7uw3XPLOFP7xdaH/XGdGHrd+7jwkc/4ZOicn71jWO594IxREX2\n/J/bnv8Oe4ictHhe/+5ULhg3gN/P3cBNzy1lb12j12EZYw7xWsF2Zjy6gOr9zbx482SumpzrdUhh\nYwmlG0mIjeKPl4/nfy88hk+K9vD1RxawYvter8MyxgD7G1v48Wsr+MmslUzISeWdO05iytB0r8MK\nK0so3YyIcO2Uwbx26wkAXPz4pzw2v8h6gDTGQxt27WPGowv4x2fF3H7GcF68eTL9Ervf7ec7yxJK\nNzUuJ4W3bz+Js8f058F/r+fKpxZRXFnndVjG9CqtrcrTH2/ivD8tYE9NI89/cxI/PHMEkd2ol8Vg\nsoTSjaXEx/Dnqybw+0vHsbakmnMe/pg3l++wq+uNCYPtFXVc+dQi7n97HaeMyGTOf32Nk4dneh2W\np+zCxm5ORLj4uEFMGpLGD15Zzh0zl/POqlL+d8Yx9EvqfVVuY0JNVXmtoJj7/uX0pvHbS8ZyyXGD\nEOmdtRJfdvv6HqS5pZWnF2zmobkbiImK4P99fTSX5tuObkywbN1Ty8/+uZqPC8uZMjSN3106jkGp\n8V6HFXKB3r7eEkoPtLm8ljv/sZIlmys4aVgGv/rGseSm9/yd3phQaWxu5amPN/HI+4VER0bwk7NH\ncu2UwUT0kmMlXbo/FBFJE5G5IlLoPn+lz0sROU1Elvs86kXkQnfasyKy2Wfa+PC/i65rSEYCM781\nhfsvPIbl2/dy5kMf8cd5G6hvavE6NGO6nWVbKzj/Twv47Zz1nDayH/N+eArXT83rNcmkIzypoYjI\ng0CFqj4gIncBqap652HKpwFFwCBVrRORZ4F/qeqsjqy3t9RQfJVW7eeXb6/jXytLyUnrw8/PG8O0\no/tZM5gx7dhdXc+Dc9Yza1kx2clx3DfjmB7TEVZHdekaCjADeM4dfg64sJ3ylwDvqqqdF9tB2cl9\n+PNVE/n7zZOJi4rkW88XcOOzS9loN5o0xq/6phYe/bCIU383nzeX7+CWrw1l7g9P6bXJpCO8qqHs\nVdUUn9eVqvqVZi+f6R8Af1DVf7mvnwVOABqA94G7VLWhjXlvAW4ByM3NPW7r1q1Bex/dTVNLK899\nuoU/zitkf1MLV07K4Y4zRpCZGOt1aMZ4TlV5Z9VOfvXOOnbs3c9Zo7P46blHk5eR4HVonvP8oLyI\nzAP89W15D/BcoAlFRLKBlcAAVW3yGbcTiAGeBDaq6n3txdQbm7z8Ka9p4JH3C/n74m3ERkVwy9eO\n4uaTh5AQa2eRm97p043l/G7Oej7btpdR/RP5+XmjmTosw+uwugzPE8phVyqyHjhVVUvd5DBfVUe2\nUfYOYIyq3tLG9FOBH6vqee2t1xLKwTaV1fDbOet5d/VOMhNjue20YVx+fA5x0ZFeh2ZMWHy+rZLf\nvbeeT4r20D8pjjumDeey/Jxee6V7W7p6QvktsMfnoHyaqv53G2UXAXer6oc+47LdZCTAQ0C9qt7V\n3notofi3bGslv3n3C5ZsqaB/UhzfO+0oLjs+h9goSyymZ1pXWs3v39vAvHW7SE+I4TunHsU1Uwbb\nn6k2dPWEkg68CuQC24BLVbVCRPKBW1X1ZrdcHvAJkKOqrT7zfwBkAgIsd+dp9yizJZS2qSqfbtzD\nQ3M3ULC1kuzkOL57qiUW07N8tq2Sxz7cyLx1u0iMi+LbXxvKjSdac297unRC8YollPYdmlgyE2O5\n8cQ8rp48mOQ+Pbu3OdMzqSofF5bz2PwiFm2qILlPNDdMzeObJw7p8T0oBoslFD8soQTuQGJ54qON\nfFxYTkJMJFdOyuWbJw1hQA/uE9v0HM0trcxZs4snPtrIqh1VZCXF8q2Th3LlpFyrkXSQJRQ/LKEc\nmTUlVTz1n028tbIUAc4bm811U/OYkJNiF0iaLqeitpGXl2zjxUVbKa2qJy89nltPOYpvTBxozbdH\nyBKKH5ZQOqe4so5nFmzhtYLt7Gto5piBSVx3Qh4XjBtgBzON51bvqOK5T7fw5ooSGptbOWlYBtdP\nzeP0Uf3srK1OsoTihyWU4KhtaOaNz3fw/MItbNhVQ0p8NJfl53BZ/iCG9Uv0OjzTi9Q2NPP2ylJe\nKdjOsq2V9ImO5KKJA7l+ah4jsmxfDBZLKH5YQgkuVWXx5gqeX7iF99bsorlVmZibwmX5OXx9bDaJ\ncXbA0wSfqvLZtkpeWbqdf60spa6xhaGZCVw1KZdLj8uxA+0hYAnFD0sooVO2r4F/fr6DVwq2U7S7\nhj7RkZxzbH8umjCIE45KtyYH02k79u7nrRUlvFawnY1ltcTHRHLe2GwuPz6HibmpdjwvhCyh+GEJ\nJfRUleXb9/JqQTFvrSihpqGZjL6xnDc2m/PHDWBirh3IN4Er29fAO6tKmb2ihGVbKwGYmJvC5cfn\n8PWxA+hrZ2uFhSUUPyyhhFd9UwsffLGb2ctL+GD9bhqbWxmU2ofzxg7grDFZjB+UYn1KmK/YU9PA\nvHW7eGtFKZ9uLKdVYVT/RM4fN4DzxmYzON1u1hhullD8sITiner6Juau2cXsFSUsKCqnpVXJTIxl\n2tFZnDU6ixOOSrczxXqxTWU1zF27i3nrdrFsayWtCoPT47lg3ADOHzfADrB7zBKKH5ZQuoaquiY+\nXL+buWt3MX/9bmobW4iPieSUEZlMOzqLk0dk0C8xzuswTQg1tbTy+ba9fPDFbuau3cnGsloARmcn\nceboLM4cncWYAUnWPNpFWELxwxJK19PQ3MKnG/c4/07X7mL3Pqdbm1H9EzlxWAYnDc9g8pA04mOs\nrbw7U1U2ltWyoLCMjwvLWbRpD7WNLURFCFOGpnPm6Cymjc5ioN2FoUuyhOKHJZSurbVVWVtazceF\n5SwoKmPplkoam1uJiYxg4uAUTh6eyfF5aYwdlGzNY12cqrKtoo6lWypZvGkPC4rKKa2qB5ymrJOH\nZ3DSsEymDksnyU4v7/IsofhhCaV72d/YwtItFSwoKufjwnLWlVYDEBMZwTEDk8jPSyN/cCrHDU4l\nva/1Ouml5pZWvti5j6VbKijYUsnSLRVf1jaT4qK+rG2ePCyT3PR4j6M1HWUJxQ9LKN1bRW0jy7ZW\nUrClgoKtlawqrqKxxenVYGhGAuNzUhgzMJljBiQxekCSXVgZIq2tyqbyWlbvqGLVjipWFVexpqSK\n2sYWAAam9CE/L5X8vDQm5aUxvF9fO5uvm7OE4ocllJ6lvqmF1TuqWLqlkmVbK1i1o4pd1Q1fTs9L\nj3cTTDKjByQxrF9fBiTH2YHeDthX30Th7hoKd+1j/c4aVpdUsbakmpqGZgBioyIYPSCJYwcmc9xg\nJ4nYcZCeJ9CE4smRThG5FLgXOBqYpKp+f+VFZDrwMBAJPK2qD7jjhwAzgTTgM+BaVW0MQ+imC4mL\njnSavfLSgKMA2L2vnjUl1azZUcWakmpWFu/l7ZWlX84THxPJUZl9OSozgWH9+n75GJQa32uPy7S0\nKruq69m6p46te2op2l3DBjeJHDjuARAXHcHR2UlcNHEgxw5M5thByQzL7EtUZISH0ZuuxKseG48G\nWoG/4PQH/5WEIiKRwAbgTKAYWApcqaprReRV4HVVnSkiTwArVPXx9tZrNZTeaW9dI1/s3MfGshqK\ndjuPjbtrKPH5sQTolxjLoNQ+DEqNP+h5QEofMvvGktQnqlvWbuqbWthd3cCuffXsqq5nZ1U92yrq\n2Lqnju0VdRRX7v+y6RCcWsewfn0ZkZXI8Ky+DO+XyIgsJ+naLXR6py5dQ1HVdUB7X85JQJGqbnLL\nzgRmiMg64HTgKrfcczi1nXYTiumdUuJjmDI0nSlD0w8aX9vQzKayWjaW1bC9oo7tlc6P6/Lte3ln\nVSnNrQf/2YqJjCCjbwyZibFk9I398jm5TzSJcVEkxh14doaT+kTRNzaKmMgIIiOkU8mopVXZ39RC\nfVML+xud5+r6JvbWNVG133neu7+J6v1N7K1rZE9tI7uq69lV3UDV/qavLC8xLorB6fGMyk7krDH9\nyU2L//IxMLWPJQ5zRLryyf0Dge0+r4uByUA6sFdVm33GD2xrISJyC3ALQG5ubmgiNd1SQmwUxw5y\nmm4OdaAZqLhyPyV791Ne00BZTQNl+xoor2mktKqelTuq2FPTQGsAlXwRJyHFREUQGxVBTGQE0VER\nCKCAKrSqouqccqtAc6tS39RCQ1PrQTWIw60jKS6a5D7RpCbEMCQjgSlD08lKiqNfYixZSXFkJcXR\nPynO7shrQiJkCUVE5gH9/Uy6R1XfDGQRfsbpYcb7papPAk+C0+QVwHqNITJCGJDSp93ujltblZrG\nZvbVN7Ovvumg5+r6ZmobmmlsbqWppZXG5lYamp3k0OQ+q0KEOLV1ARCIcIcjI4S46Ej6xEQSFxVJ\nn5gI4qIjv3wkxUWREh9DSp9oUuKjSYyLtpqF8VTIEoqqTuvkIoqBHJ/Xg4ASoBxIEZEot5ZyYLwx\nYRcRISTFRbsX59nZTaZ368qnZywFhovIEBGJAa4AZqtzFsGHwCVuueuBQGo8xhhjQsiThCIi3xCR\nYuAE4G0RmeOOHyAi7wC4tY/bgDnAOuBVVV3jLuJO4IciUoRzTOWv4X4PxhhjDmYXNhpjjDmsQE8b\n7spNXsYYY7oRSyjGGGOCwhKKMcaYoLCEYowxJigsoRhjjAmKXnWWl4iUAbU4F0d2NRlYXIHqijGB\nxdURXTEmsLjaMlhVM9sr1KsSCoCIFARy+lu4WVyB64oxgcXVEV0xJrC4OsuavIwxxgSFJRRjjDFB\n0RsTypNeB9AGiytwXTEmsLg6oivGBBZXp/S6YyjGGGNCozfWUIwxxoSAJRRjjDFB0SMTiohcKiJr\nRKRVRNo81U5EpovIehEpEpG7fMYPEZHFIlIoIq+4/bEEI640EZnrLneuiKT6KXOaiCz3edSLyIXu\ntGdFZLPPtPHhisst1+Kz7tk+44O+vQLcVuNFZKH7Wa8Ukct9pgV1W7W1r/hMj3Xfe5G7LfJ8pt3t\njl8vImd3Jo4OxvRDEVnrbpv3RWSwzzS/n2WY4rpBRMp81n+zz7Tr3c+8UESuD3NcD/nEtEFE9vpM\nC8n2EpFnRGS3iKxuY7qIyCNuzCtFZKLPtJBtqyOmqj3uARwNjATmA/ltlIkENgJDgRhgBTDanfYq\ncIU7/ATwnSDF9SBwlzt8F/CbdsqnARVAvPv6WeCSEGyvgOICatoYH/TtFUhMwAhguDs8ACgFUoK9\nrQ63r/iU+S7whDt8BfCKOzzaLR8LDHGXExmmmE7z2Xe+cyCmw32WYYrrBuDPbezvm9znVHc4NVxx\nHVL++8AzYdheXwMmAqvbmH4u8C5O1+dTgMWh3ladefTIGoqqrlPV9e0UmwQUqeomVW0EZgIzRESA\n04FZbrnngAuDFNoMd3mBLvcS4F1VrQvS+tvS0bi+FMLt1W5MqrpBVQvd4RJgN9Du1bxHwO++cph4\nZwFnuNtmBjBTVRtUdTNQ5C4v5DGp6oc++84inO6yQy2QbdWWs4G5qlqhqpXAXGC6R3FdCbwcpHW3\nSVX/g/OnsS0zgOfVsQin+/NsQrutjliPTCgBGghs93ld7I5LB/aq02Ok7/hgyFLVUgD3uV875a/g\nqzv1L92q70MiEhvmuOJEpEBEFh1ohiN026tD20pEJuH889zoMzpY26qtfcVvGXdbVOFsm0DmDVVM\nvm7C+ad7gL/PMhgCjeti97OZJSI5HZw3lHHhNg0OAT7wGR2q7dWetuIO5bY6YlFeB3CkRGQe0N/P\npHtUNZA+5sXPOD3M+E7HFegy3OVkA8fidIF8wN3ATpwfzidxukK+L4xx5apqiYgMBT4QkVVAtZ9y\nAW2vIG+rF4DrVbXVHX3E28rfKvyMO/Q9hmR/OoyAlysi1wD5wCk+o7/yWarqRn/zhyCut4CXVbVB\nRG7FqdmdHuC8oYzrgCuAWara4jMuVNurPeHerzql2yYUVZ3WyUUUAzk+rwcBJTg3YEsRkSj3n+aB\n8Z2OS0R2iUi2qpa6P4K7D7Ooy4A3VLXJZ9ml7mCDiPwN+HE443KblVDVTSIyH5gA/IMj3F7BiElE\nkoC3gZ+5TQIHln3E28qPtvYVf2WKRSQKSMZpyghk3lDFhIhMw0nQp6hqw4HxbXyWwfiBbDcuVd3j\n8/Ip4Dc+8556yLzzgxBTQHH5uAL4nu+IEG6v9rQVdyi31RHrzU1eS4Hh4pyhFIOzE81W54jXhzjH\nLwCuBwKp8QRitru8QJb7lTZc94f1wHGLCwG/Z4aEIi4RST3QbCQiGcCJwNoQbq9AYooB3sBpY37t\nkGnB3FZ+95XDxHsJ8IG7bWYDV4hzFtgQYDiwpBOxBByTiEwA/gJcoKq7fcb7/SyDEFOgcWX7vLwA\nWOcOzwHOcuNLBc7i4Bp6SONyYxuJc5B7oc+4UG6v9swGrnPP9poCVLl/lkK5rY6c12cFhOIBfAMn\ngzcAu4A57vgBwDs+5c4FNuD807jHZ/xQnC99EfAaEBukuNKB94FC9znNHZ8PPO1TLg/YAUQcMv8H\nwCqcH8cXgb7higuY6q57hft8Uyi3V4AxXQM0Act9HuNDsa387Ss4TWgXuMNx7nsvcrfFUJ9573Hn\nWw+cE8T9vL2Y5rn7/4FtM7u9zzJMcf0aWOOu/0NglM+833S3YRFwYzjjcl/fCzxwyHwh2144fxpL\n3f24GOdY163Are50AR51Y16Fz1mrodxWR/qwW68YY4wJit7c5GWMMSaILKEYY4wJCksoxhhjgsIS\nijHGmKCwhGKMMSYoLKEYY4wJCksoxhhjgsISijEeEpHj3ZskxolIgjh9uxzjdVzGHAm7sNEYj4nI\n/ThX2vcBilX11x6HZMwRsYRijMfce0stBeqBqXrwXW6N6TasycsY76UBfYFEnJqKMd2S1VCM8Zg4\nfZTPxOnUKVtVb/M4JGOOSLftD8WYnkBErgOaVfXvIhIJfCoip6vqB+3Na0xXYzUUY4wxQWHHUIwx\nxgSFJRRjjDFBYQnFGGNMUFhCMcYYExSWUIwxxgSFJRRjjDFBYQnFGGNMUPx/AB3e1vLPbdcAAAAA\nSUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Here we transform from polar to cartesian coordinates \n", "# to then plot\n", "cx = [it[0]*np.sin(it[2]) for it in y]\n", "cy = [it[0]*np.cos(it[2]) for it in y]\n", "plt.plot(cx,cy)\n", "plt.title(\"Orbit resulting from the chosen initial conditions\")\n", "plt.xlabel(\"x\")\n", "plt.ylabel(\"y\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### We perform the numerical integration using gduals (to get a HOTM)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# Order of the Taylor Map. If we have 4 variables the number of terms in the Taylor expansion in 329 at order 7 \n", "order = 5\n", "# We now define the initial conditions as gdual (not float)\n", "ic_g = [gdual(ic[0], \"r\", order), gdual(ic[1], \"vr\", order), gdual(ic[2], \"t\", order), gdual(ic[3], \"vt\", order)]" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "--- 1.7251091003417969 seconds ---\n" ] } ], "source": [ "import time\n", "start_time = time.time()\n", "# We call the numerical integrator, this time it will compute on gduals\n", "y = rk4(eom_kep_polar, it, ic_g, step, n_steps)\n", "print(\"--- %s seconds ---\" % (time.time() - start_time))\n" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\\[ 456.247{dvt}^{3}-0.00244413{dt}^{3}+32.5178{dr}{dt}^{3}{dvt}+9.72205{dr}{dt}^{3}{dvr}-98350.9{dr}^{2}{dvr}^{3}+95097.4{dr}^{3}{dt}{dvt}+3.44559{dt}^{3}{dvt}-14318.7{dvt}^{4}+680.557{dr}{dvt}+3059.27{dr}{dt}{dvr}{dvt}-501710{dr}{dvt}^{4}-2.01229{dvt}-4.00992{dr}-58.3323{dr}{dt}{dvr}+2119.32{dr}^{2}{dvr}-1040.11{dvr}^{2}{dvt}^{2}-93373.3{dr}^{2}{dvr}{dvt}-61205.4{dr}^{3}{dvr}-185943{dr}^{4}+1.00615{dt}^{2}{dvt}+\\ldots+\\mathcal{O}\\left(6\\right) \\]\n", "xf (latex):\n" ] }, { "data": { "text/latex": [ "\\[ 456.247{dvt}^{3}-0.00244413{dt}^{3}+32.5178{dr}{dt}^{3}{dvt}+9.72205{dr}{dt}^{3}{dvr}-98350.9{dr}^{2}{dvr}^{3}+95097.4{dr}^{3}{dt}{dvt}+3.44559{dt}^{3}{dvt}-14318.7{dvt}^{4}+680.557{dr}{dvt}+3059.27{dr}{dt}{dvr}{dvt}-501710{dr}{dvt}^{4}-2.01229{dvt}-4.00992{dr}-58.3323{dr}{dt}{dvr}+2119.32{dr}^{2}{dvr}-1040.11{dvr}^{2}{dvt}^{2}-93373.3{dr}^{2}{dvr}{dvt}-61205.4{dr}^{3}{dvr}-185943{dr}^{4}+1.00615{dt}^{2}{dvt}+\\ldots+\\mathcal{O}\\left(6\\right) \\]" ], "text/plain": [ "456.247*dvt**3-0.00244413*dt**3+32.5178*dr*dt**3*dvt+9.72205*dr*dt**3*dvr-98350.9*dr**2*dvr**3+95097.4*dr**3*dt*dvt+3.44559*dt**3*dvt-14318.7*dvt**4+680.557*dr*dvt+3059.27*dr*dt*dvr*dvt-501710*dr*dvt**4-2.01229*dvt-4.00992*dr-58.3323*dr*dt*dvr+2119.32*dr**2*dvr-1040.11*dvr**2*dvt**2-93373.3*dr**2*dvr*dvt-61205.4*dr**3*dvr-185943*dr**4+1.00615*dt**2*dvt+..." ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# We extract the last point\n", "yf = y[-1]\n", "# And unpack it into some convinient names\n", "rf,vrf,tf,vtf = yf\n", "# We compute the final cartesian components\n", "xf = rf * sin(tf)\n", "yf = rf * cos(tf)\n", "# Note that you can get the latex representation of the gdual\n", "print(xf._repr_latex_())\n", "print(\"xf (latex):\")\n", "xf" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Final x from the gdual integration -0.9089833754123715\n", "Final y from the gdual integration 0.014664786455128188\n", "\n", "Final x from the float integration -0.9089833754123692\n", "Final y from the float integration 0.014664786455132996\n" ] } ], "source": [ "# We can extract the value of the polinomial when $\\mathbf {dy} = 0$\n", "print(\"Final x from the gdual integration\", xf.constant_cf)\n", "print(\"Final y from the gdual integration\", yf.constant_cf)\n", "# And check its indeed the result of the 'reference' trajectory (the lineariation point)\n", "print(\"\\nFinal x from the float integration\", cx[-1])\n", "print(\"Final y from the float integration\", cy[-1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### We visualize the HOTM" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "--- 0.34252142906188965 seconds ---\n" ] } ], "source": [ "# Let us now visualize the Taylor map by creating a grid of perturbations on the initial conditions and\n", "# evaluating the map for those values\n", "Npoints = 10 # 10000 points\n", "epsilon = 1e-3\n", "grid = np.arange(-epsilon,epsilon,2*epsilon/Npoints)\n", "nxf = [0] * len(grid)**4\n", "nyf = [0] * len(grid)**4\n", "i=0\n", "import time\n", "start_time = time.time()\n", "for dr in grid:\n", " for dt in grid:\n", " for dvr in grid:\n", " for dvt in grid:\n", " nxf[i] = xf.evaluate({\"dr\":dr, \"dt\":dt, \"dvr\":dvr,\"dvt\":dvt})\n", " nyf[i] = yf.evaluate({\"dr\":dr, \"dt\":dt, \"dvr\":dvr,\"dvt\":dvt})\n", " i = i+1\n", "print(\"--- %s seconds ---\" % (time.time() - start_time))" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0.5,1,'Stretch')" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA34AAAE/CAYAAAAZshH0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3XuclXW5///XNTOcD8NpVGA4Ch4Q\nbAwEMyPNUNRdYGohplYm6c5vm11770jL9sbDj87utqZZaFkimgZRYkqpEaIjEKOASCDKMA4i5/Np\nZq7fH/c9uBjWzKyZWbPudXg/H4/1WGvd9+e+72uVfOa+7s/J3B0RERERERHJXnlRByAiIiIiIiKt\nS4mfiIiIiIhIllPiJyIiIiIikuWU+ImIiIiIiGQ5JX4iIiIiIiJZTomfiIiIiIhIllPiJ0lnZv9t\nZr+NOg4RERERSZyZuZkNiToOaR1K/KTJzGxvzKvGzA7EfL8m6vhERBpjZtfUqctqX25mt0cdn4jk\nDjM7z8wWm9kuM9tuZi+Z2dlm9gUzW9TCc59vZhXJilUymxI/aTJ371z7AsqBT8VsezTq+EREGuPu\nj8bWZWF9NhXYDPwi4vBEJEeYWVfgT8D/AT2AvsD/AIcSPD6/9aKTbKPET1pLWzN7xMz2mNkqMxtV\nu8PM+pjZU2a2xczeNrOv1XcSM/uVmf3MzJ4Jn8a/ZGYnmdk9ZrbDzN40s7Niyk8zs7fC675hZpfH\n7PtCePz/hU/V3jSzC1vvfwIRyRRhPfITYJK7bwrrqXnh0/d1ZnZjTNl2YR1UGb7uMbN24b7zzazC\nzP7LzN43s01mNtHMLjWzf4bnuzWq3ykiaecUAHd/zN2r3f2Auz8HHAEeAD4S3v/shKP3Rfeb2Xwz\n2wdcENZJPzSzcjPbbGYPmFkHM+sEPAP0ienV0MfM8s3s1pj7pWVm1i8mpk+a2drwPus+M7MU/28i\nrUSJn7SWTwOzgW7APOBeADPLA/4IvEbwVOtCYKqZXdzAuT4LfBvoRfAE7GXgH+H3J4Efx5R9C/gY\nUEjwxOy3ZtY7Zv8YYH147HeB35tZj5b8UBHJbGbWjaAuudPdXww3PwZUAH2AK4G7Yx4U3QacA5QA\nHwJGE9RRtU4C2hPUcbcTtCB+HhhJUD/dbmaDW/EniUjm+CdQbWa/NrNLzKw7gLuvBm4CXg57JXSL\nOWYycBfQBVgEfI8ggSwBhhDWPe6+D7gEqIzp3VAJfB24GrgU6Ap8Cdgfc/5/Ac4mqN8+CzR0jyYZ\nRImftJZF7j7f3auB3xBUHhBUJEXuPt3dD7v7eoKbokkNnGuOuy9z94PAHOCguz8Snvtx4GiLn7v/\nzt0r3b3G3R8H1hLclNV6H7jH3Y+E+9cAlyXpN4tIhgmfZP8aWAl8P9zWDzgP+Ka7H3T3MuCXwLXh\nYdcA0939fXffQvCQ6dqY0x4B7nL3IwQPwHoB/+vue9x9FbAKOLP1f52IpDt3301Q3zjB/dCWsLfB\niQ0c9gd3f8ndawgeiN8I/Lu7b3f3PcDdNHxf9WXg2+6+xgOvufu2mP0z3H2nu5cDLxAklJIFCqIO\nQLLWezGf9wPtzawAGEDQ5WBnzP584O8NnGtzzOcDcb53rv1iZtcRPMkaGG7qTHDTVetdd/eY7xsI\nnuiLSG76JjAcGBlTN/QBam+gam0ARsXs31BnX2w9si18MAVBHQUN1FsiktvC1r0vAJjZacBvgXuA\nZ+s5ZGPM5yKgI7AspkemEdxb1acfQQ+p+tS9h1N9lSXU4iepthF42927xby6uPulLT2xmQ0geFp2\nC9Az7BaxkqACrNW3Tl/1/kBlS68tIpnHzM4n6LZ5pbvHPoyqBHqYWZeYbf2Bd2P2D6izT/WIiLSY\nu78J/IrggZTXVyzm81aCh0lnxNxXFYYTVtUtW2sjcHKSQpYMosRPUu1VYLeZfTMceJxvZsPN7Owk\nnLsTQQW3BcDMvkhQccY6AfiambUxs6uA04H5Sbi2iGSQcOzvbGCquy+P3efuG4HFwP9nZu3N7Ezg\nBqB21uLHgG+bWZGZ9SIYx6e1S0WkyczsNDP7hpkVh9/7EYy/e4Wgp0CxmbWt7/iwu+cvgJ+Y2Qnh\nOfrGzJ2wGehpZoUxh/0SuMPMhlrgTDPrmfxfJ+lGiZ+kVNj96VME/cXfJnhS9UuCyVhaeu43gB8R\nTP6yGRgBvFSnWCkwNLzuXQRP+rchIrnmRuBE4H/t+LX8HiC48RpI0JI3B/iuuy8Ij70TWAq8Dqwg\nmGzqzlT/ABHJCnsIJp4rDWfpfIWgt9I3gOcJxgS/Z2ZbGzjHN4F1wCtmthv4C3AqHG1BfAxYb2Y7\nzawPwaR4TwDPAbuBmUCHVvhtkmbs2OFOItnLzL4AfNndz4s6FhERERGRVFKLn4iIiIiISJZT4ici\nIiIiIpLl1NVTREREREQky6nFT0REREREJMsp8RMREREREclyBVEH0By9evXygQMHRh2GiCTRsmXL\ntrp7UdRxtITqJpHskw11E6h+EslGTa2fMjLxGzhwIEuXLo06DBFJIjPbEHUMLaW6SST7ZEPdBKqf\nRLJRU+sndfUUERERERHJckr8REREREREspwSPxERERERkSynxE9ERERERCTLKfETERERERHJckr8\nREREREREspwSPxERERERkSyXlMTPzB4ys/fNbGU9+83Mfmpm68zsdTP7cMy+681sbfi6PhnxiIiI\niIiIyAeS1eL3K2B8A/svAYaGrynA/QBm1gP4LjAGGA1818y6JykmERERERERAQqScRJ3X2hmAxso\nMgF4xN0deMXMuplZb+B8YIG7bwcwswUECeRjyYhL0lx1FezfFrwO74XD+4LXkf1QUxWUcQ/eLQ/a\ndoQ24attJ+jYEzr1gjYdovsNIiIiIiIZICmJXwL6AhtjvleE2+rbLtnq6f+AdxbBvi1Bwoe3/Jxt\nOwcJYJc+0H0AdBsQvHcfBCecBh3UiCwiIiIiuS1ViZ/F2eYNbD/+BGZTCLqJ0r9//+RFJqnVpgP0\nPBn6nwOdT4BORUHLXbuuQYte207QphPkx/6naUEL4JEDQWvgkf1By+D+bUECuW8r7H0fdlfC2wuD\n99j/jLr0gROHwQmnQ9+R0HcUFBaDxfvPT0REREQk+6Qq8asA+sV8LwYqw+3n19n+YrwTuPuDwIMA\no0aNSkIzkUTiojta/xpVh2DnRtj+Fry/Onytgrf/DtWHgjKdT4LiUTDgXBj0cThhGORpklsRERER\nyU6pSvzmAbeY2WyCiVx2ufsmM3sWuDtmQpeLgG+lKCbJVgXtoNeQ4HXKxR9srzoMm1fCu8ugYgls\nfBXe/FOwr2MvGDQWhlwIp4wPuo6KiIiIiGSJZC3n8BjwMnCqmVWY2Q1mdpOZ3RQWmQ+sB9YBvwD+\nFSCc1OUOYEn4ml470YtI0hW0hb4fhtE3wmcehH8rg39fBRN+FiR8GxbDH74KPxgCMy+Gl/4Xtr0V\nddSSADMbb2ZrwiVjpsXZP9bM/mFmVWZ2ZZ19cZeUMbORZrYiPOdPzdQ3WESaTvWTiKSLZM3qeXUj\n+x34aj37HgIeSkYcIk1WWAxnXRO83GHTa7DmGVjzNCy4PXj1+TCc+Vk44zPQ5cSoI5Y6zCwfuA8Y\nR9B9fImZzXP3N2KKlQNfAP6jzrG1S8qMIhgYuiw8dgfBsjNTgFcIHl6NB55p3V8jItlE9ZOIpBMN\nahKpZQZ9SuCCb8FNi2DqSrjozmBimT9Pgx+fBr+5HN6YB9VHoo5WPjAaWOfu6939MDCbYAmZo9z9\nHXd/Haipc+zFhEvKhDdTC4Dx4XIzXd395fDB1SPAxFb/JSKSbVQ/iUjaSNUYP5HM060fnPv/gteW\nNfD6E/DaY/DEtdClN3z4uuBVWBx1pLku3rIwY1pwbO1SMxVxtouINIXqJxFJG2rxE0lE0alw4Xfg\n316HSY/BicPhb9+He0bA7GuCiWIkKgkvC9OEY5u01IyZLTWzpVu2bEnwsiKSI1Q/iUjaUOIn0hT5\nBXDapfD5J+HfXoOPToUNL8HMcfDwpfDP54KxgpJK9S0X05JjK8LPjZ7T3R9091HuPqqoqCjhoEUk\nJ6h+EpG0ocRPpLm6D4BPfjcYCzh+Bux4B2ZdBQ+cByufgpq6wzWklSwBhprZIDNrC0wiWEImEc8C\nF5lZ93BZmYuAZ919E7DHzM4JZ8u7DvhDawQvIllN9ZOIpA0lfiIt1a4znHMzfK0MJt4fTPzy5Jfg\nwbGw9i9qAWxl7l4F3EJwk7QaeMLdV5nZdDP7NICZnW1mFcBVwM/NbFV4bENLytwM/JJgGZq30Ix5\nItJEqp9EJJ2YZ+BN6ahRo3zp0qVRhyESX0110OL3/J2wcwMMOC9oGew3OurI0pqZLXP3UVHH0RKq\nm0SyTzbUTaD6SSQbNbV+UoufSLLl5Qfr/t2yFC79IWz9ZzAG8InrYOfGxo8XEREREUkyJX4iraWg\nLYy+Ef6tDM6/NZj45b7RsPCHUHUo6uhEREREJIco8RNpbW07wfnfhFtehSEXwvN3wM/OgbULoo5M\nRERERHKEEj+RVOnWHz73W/j878Hy4NEr4akbYf/2xo8VEREREWkBJX4iqTbkQrj5Zfj4NFj1e7hv\nDKz+U9RRiYiIiEgWU+InEoWCtnDBt+DGF6DzifD4NfDkDWr9ExEREZFWocRPJEq9z4QpLwSTv7wx\nF+7/KLzzUtRRiYiIiEiWUeInErX8NsHkLzc+D206wK//BV78XrAeoIiIiIhIEijxE0kXvT8EX/kb\njLgKXrwbHpkAuzdFHZWIiIiIZAElfiLppF0XuPznMOFn8O4yeOA8df0UERERkRZT4ieSbszgrGtg\nyovQoRs88mlYMjPqqEREREQkgynxE0lXRafCl/8Kgy+Ap78Of5wKVYejjkpEREREMpASP5F01qEb\nTH4cPjoVlj0cjPvTkg8iIiIi0kRK/ETSXV4+jPsf+Mwvg3F/My+CneVRRyUiIiIiGUSJn0imOPMq\nuG4u7HsffjkO3lsZdUQiIiIikiGSkviZ2XgzW2Nm68xsWpz9PzGzsvD1TzPbGbOvOmbfvGTEI5K1\nBpwLX/wzWB48fAm8/feoIxIRERGRDNDixM/M8oH7gEuAYcDVZjYstoy7/7u7l7h7CfB/wO9jdh+o\n3efun25pPCJZ78Rh8OUF0KU3/PYzsOaZqCMSERERkTSXjBa/0cA6d1/v7oeB2cCEBspfDTyWhOuK\n5K7CYvjSn+HE4fD4tbD6T1FHJCIiIiJpLBmJX19gY8z3inDbccxsADAIeD5mc3szW2pmr5jZxCTE\nI5IbOvYIxvz1KYHfXQ9v/CHqiEREREQkTSUj8bM427yespOAJ929OmZbf3cfBUwG7jGzk+NexGxK\nmCAu3bJlS8siFskW7Qvh87+HviPhd1+ElU9FHZGIiIiIpKFkJH4VQL+Y78VAZT1lJ1Gnm6e7V4bv\n64EXgbPiHejuD7r7KHcfVVRU1NKYRbJH+67w+aeg3xh46kZ4c37UEYmIiIhImklG4rcEGGpmg8ys\nLUFyd9zsnGZ2KtAdeDlmW3czaxd+7gV8FHgjCTGJ5JZ2XeCa34XdPr+g2T5FRERE5BgtTvzcvQq4\nBXgWWA084e6rzGy6mcXO0nk1MNvdY7uBng4sNbPXgBeAGe6uxE+kOdp1hmuehB6D4LGroXJ51BGJ\n5Kyps5dTMv05ps7Wv0MRkUyXLXV6QTJO4u7zgfl1tt1e5/t/xzluMTAiGTGICMGEL9fOgZkXw2+v\ngC89C72GRh2VSE6ZOns5c8uCEQ+17/dMijuKQURE0tzEexdRVrELCOr0v6zezK2XDmPymP4RR9Z0\nSVnAXUTSSNc+wWyflhes87dXkyGJpNIfyiob/C4iIplhVmn50aSv1t5D1dw6ZwUT710UUVTNp8RP\nJBv1PBkmPw5734fHr4EjB6OOqFWZ2XgzW2Nm68xsWpz97czs8XB/qZkNDLdfY2ZlMa8aMysJ970Y\nnrN23wmp/VWSqepOa13fNNeSG1Q/iWSuh156u959ZRW7OG/GX1MYTcsp8RPJVn1HwuU/h42lMO//\ngWfn7aeZ5QP3AZcAw4CrzWxYnWI3ADvcfQjwE+B7AO7+qLuXuHsJcC3wjruXxRx3Te1+d3+/1X+M\niGQV1U8iGa6Re6eKnQcpmf4cs0rLUxRQyyjxE8lmZ0yET3wHVjwBC38QdTStZTSwzt3Xu/thYDYw\noU6ZCcCvw89PAheaWd01SK+mznIzIsmSKTcFknSqn0QyxKzScq6dWXpMff2l8wY3etzO/Ue4dc4K\nrptZ2prhJYUSP5Fs97FvwIeuhhfugtV/ijqa1tAX2BjzvSLcFrdMOBPxLqBnnTKf4/gbq4fDblTf\niXMjJnKc+hK8+15cl+JIJE2ofhLJANfNLOXWOSv4+9qt3DpnxdHZOyeP6c/Ekj4JnWPh2q2c/4MX\nWLZhR2uG2iJK/ESynRl86n+hz4dh7s2w7a2oI0q2eDc8dftmNFjGzMYA+919Zcz+a9x9BPCx8HVt\n3IubTTGzpWa2dMsWTaST6+6eH39Fok07DqQ4EkkTqp9E0tx1M0tZuHbrMdvmllUy5ZGlLNuwg3sm\nncVNYwfH/Yda1zvb9nPF/YvTdtkHJX4iuaCgHXz215CXD09cD0ey6ia0AugX870YqDuN4tEyZlYA\nFALbY/ZPos7TdHd/N3zfA8wi6LJ1HHd/0N1HufuooqKiFvwMyQZ7D1XH3V6T4jgkbah+Ekljs0rL\nj0v6aj33xmY+9/OXmVVazrRLT+fJm8/l7IHdEzrv3LLKtJz1U4mfSK7o1h8+8wvYvBKe/o+oo0mm\nJcBQMxtkZm0JbpLm1SkzD7g+/Hwl8Lx7MGLbzPKAqwjG3hBuKzCzXuHnNsC/ACsREWka1U8iaeyZ\nlZsa3F9V49w2ZwWzSssZOaA7v7vp3IS7fpZV7Eq7rp9K/ERyydBxMPY/oey3sPy3UUeTFOGYmFuA\nZ4HVwBPuvsrMppvZp8NiM4GeZrYO+DoQO6X6WKDC3dfHbGsHPGtmrwNlwLvAL1r5p4hIllH9JJI+\nlm3Ywa1zVnDbnBVHk7FLhvdu9DiHY467Z9JZ3H35CPp2a9/ose9s28+kB19Om+SvIOoARCTFzp8G\n5S/D/P+CAedCj8ZnrEp37j4fmF9n2+0xnw8SPDWPd+yLwDl1tu0DRiY9UMlq6fKHXdKL6ieR6C3b\nsIMr7l989PvsJRu5Y8JwJo/pD8BdT7/BvsPxu+rXerS0/JjjJo/pz8R7Fx23wHtdR6qdV9ZvY+SA\nxLqJtia1+Inkmrx8uPwByC+A338FqquijkgkK3zx4VejDkFEROL47AOLj/leXePcGnbhnDymP6um\nj2doUadGz1N7XG3r39xbzmPs0F6NHnfO4LoT9UZDiZ9ILioshst+DBWvwqIfRx2NSFbYfbD+hyj5\n+msrIhKJ62aWUl3POuyxSdyCb5yf8Pi9R0vLj0788sgNY7j78hH06tw2btmxQ3ulRWsfKPETyV0j\nroQRV8GLM6BiWdTRiGS1GxNYBFhERJJv8VvbGtz/aGk5Vz6wmFml5dwz6Syeuvlc2hc0niLFTvwy\neUx/ln57HCXFhceUGVrUiUduGNOi+JNJiZ9ILrv0h9Cld7C+X9WhqKMRyWgNrfE07dLTUxaHiIh8\noFO7/EbLuH/Q+gfw6I3nNHJEeBzBcbVjvOfech43jR3MwJ4duWnsYBZ84/zmht0qlPiJ5LIO3eBT\n98DWNbDoJ1FHI5LREl3fSUREWsfU2cspmf7cMQuof3N84g/eHi0t56oHFrPmvT3cffkI8hNZtR34\nztwVRz9Pu/R0XvzPC9LygZ8SP5FcN3QcDL8SFv4Q3n8z6mhEMlb7No0/VRYRkdZx3cxS5pZVsnP/\nEeaWVXLat59hxvzVTB7Tn5vGJt7dviZs/VtVuYs7Jo44OvNnQ97euq8loaeMEj8RgfEzoF1n+OPX\noKYm6mhEMtLCtVujDkFEJCct27DjuDr4YFUNDyxcz7gfvci0S0/nqZvPZUCPjgmf89HScm7/w0qG\n9ynk7stHNFh2YM/GZwRNB0r8RAQ6F8FFd8HGUlj2cNTRiGSVzm3VEigi0poaWk5n7ZZ9nPbtZ1iw\n6j3+9l8XJLRsQ63aCVxOPakLd18+goK8+H0/72wkMUwXSvxEJFAyGQZ+DJ6/A/ZvjzoakaxRzyzi\nIiKSJA0tpwPHtv4t+Mb53DR2MO0SmLkTgjr8sw8sZvKY/jz+lY8weUx/avM/A+6+fETaLNfQGCV+\nIhIwC7p8HtwVLPEgIkmx73B11CGIiGSNeBO4FHdrn9Cxa7fsY+QdzzHujJNYc+clCR9X7XDejL8y\nckB37r58BL+76Vz+8+JTefLmcxMaA5gulPiJyAdOGg4jvwhLfgnvr446GpGMEXsDIiIiraPuBC5D\nb5vP1NnLWTTtQtokmNVs23eEK+5fzHUzS1k07UImlvShIIEe+RU7Dx79PHJAd756wZCMaemrpcRP\nRI51wW3BRC9//lawsI2INGpuWWXUIYiIZLV4E7gcqXbmllVy3oy/Mvsr5ya8/AIEE3KNvOM5rv3I\nQNbddVnCrX+ZLCmJn5mNN7M1ZrbOzKbF2f8FM9tiZmXh68sx+643s7Xh6/pkxCMiLdCpJ5x/K6x/\nAdY+F3U0IiIiIlw3s7TefRU7DzLpwZf56JBeTep6uW3fEa68fzEz5q9m0bQLGTu0V71lJ5b0aVK8\n6ajFiZ+Z5QP3AZcAw4CrzWxYnKKPu3tJ+PpleGwP4LvAGGA08F0zy6w2U5FsdPYN0H0Q/PUOLe8g\n0kLZcLMgIhK1xsZLH6l2Fq7dysI17ze6/EIsBx5YuJ6ps5fzyA1jeOrmc+nRsc0xZYq7teeeSWc1\nJ+y0kowWv9HAOndf7+6HgdnAhASPvRhY4O7b3X0HsAAYn4SYRKQl8tsEXT43r4A35kQdjUhGy4ab\nBRGRVFq2YQf3vbCOZRt2NPnYip0Hmf7HVUws6dOk1r+5ZZXMKi1n5IDu/OP2i5hY0oduHdswsaQP\ni6Zd2OQ40lEyEr++wMaY7xXhtrquMLPXzexJM+vXxGMxsylmttTMlm7ZsiUJYYtIg4ZfAScMgxfu\nhuqGp0kWyWXNuTEREZH4lm3YwZX3L+YHz67hivsXc9ucFSzbsKNJvScOVtUwt6ySiu37m9T69z9/\nXHX08z2TzqLs9ouy6uFdMhK/eMMo684I8UdgoLufCfwF+HUTjg02uj/o7qPcfVRRUVGzgxWRBOXl\nwSe+DdvWwWuPRR2NSNpqaOFgERFpmqsffPmYZODR0nKuuH8xowf1bNLi6xBM4PLQS29z09jBCbX+\nHarK7uEtyUj8KoB+Md+LgWOmN3P3be5+KPz6C2BkoseKSIROvRT6joS/fR+qj0QdjUhaamzhYBER\nScyM+as5XB1/RvFb56xg9OCeTCzpk/Di6wDr3t/LAwvXs/9QVaOtf3lNmBU0EyUj8VsCDDWzQWbW\nFpgEzIstYGa9Y75+GqhdIOxZ4CIz6x5O6nJRuE1E0oEZjP0v2FUOK38fdTQiGacJ9yYiIjnvF4vW\nN7j/0dJy/lBWyawbz2ly69/csspGW/+mfGxwk86ZaVr8J8ndq4BbCBK21cAT7r7KzKab2afDYl8z\ns1Vm9hrwNeAL4bHbgTsIksclwPRwm4iki6EXBWP9Fv1EM3yKNNH0CYmPLRERyTV1J3GpTuA2w4Er\n7l/Mgm+cz01jBzer9W94n0JuGntskjd2aC+mXXp6U8LPOAXJOIm7zwfm19l2e8znbwHfqufYh4CH\nkhGHiLSCvDz46FSYMyVY1+9UTbwrkqimzCgnIpJLps5eztyyYISXAXddPoIu7fLZc6jhZRtqDZz2\nNCXFhay58xLG/ehF1m7Zl/C17/nLGl69bRz9e3bimZWbuGR475yor9UJRUQaN/wzUNgPXron6khE\nREQkw82Yv/po0gdBK96tc1bwkZPrX0A9nrKKXZxy23xmXPmhJrX+7QnHZk8e05/f3DAmJ5I+UOIn\nIonIbwPn/j8ofxk2Lok6GhEREclgDyyMP5bvuTc2Y0B+EyZZOVztXHH/Yv66ejNr7rwkobF/F59x\nUuIXyCJK/EQkMSWToW0XWPKLqCMRERGRDDWrtLzB/U4wt1xTW+HWbtnH4G89zYWnn8hNYwfHXTOu\nVjatzdcUSvxEJDHtukDJ1bBqDuzdEnU0Immvg6b0FBE5zs9eWNtomaqaIEEcWtSpSQlgjQetiU/9\no4K3Z1xGUee2x5WpO6lLLtFfJRFJ3NlfhurDsPyRqCM5hpmNN7M1ZrbOzKbF2d/OzB4P95ea2cBw\n+0AzO2BmZeHrgZhjRprZivCYn5pZlq/uI8l2IMsXApbEqH6SXDd19nJKpj/H1NnLAdh/OLHJWyBo\nxXty6cYmJ2tb9h5m4r2LWPLtcdw0djBd2ufTtX0BN40dnPUzdzYkKbN6ikiOKDoVBo2FpQ8HM33m\n5UcdEWaWD9wHjAMqgCVmNs/d34gpdgOww92HmNkk4HvA58J9b7l7SZxT3w9MAV4hmLV4PPBMK/0M\nyVC1NzIi8ah+klx33cxSFq7dCgTr6L20bivD+xYe3ZaIw9XOAwvXU9ytPWNPPaHRrqK1yip2ATDt\n0tNzOtmLpRY/EWmas2+EXRth7YKoI6k1Gljn7uvd/TAwG5hQp8wE4Nfh5yeBCxt6Qm5mvYGu7v6y\nuzvwCDAx+aFLpoudlU4kDtVPkrOWbdhxXIK3Ze9hFq7dytCiTk1afw+gYudBZr9aTklxYTLDzClK\n/ESkaU69BDr2gtdmRR1Jrb7AxpjvFeG2uGXcvQrYBfQM9w0ys+Vm9jcz+1hM+YpGziki0hjVT5Kz\nrn7w5Xr3rd2yj349OlLcrX2TzlnjQUte57b5DGlk9s544/tynRI/EWma/DYw4ipY8wzs3x51NEDc\nibs8wTKbgP7ufhbwdWCWmXVN8JzBic2mmNlSM1u6ZYsmvZEPNGU6cslaqp8kZx2ujvuf5VHr3t9L\nxc6DDOzZscmTYe09XM36rfsaTByXfHtck86ZC5T4iUjTlVwdTPKyak7UkUDwtLtfzPdioG7/u6Nl\nzKwAKAS2u/shd98G4O7LgLetTIhKAAAgAElEQVSAU8LyxY2ck/C4B919lLuPKioqSsLPkWxxx8QR\nUYcg0VP9JDljVmk5184sTXgMXq13tu3nYHUNQ4s6kdeEB2Y1HnT/LCkuPK717+7LVf/Go8ldRKTp\nTjoTThgGrz0GZ98QdTRLgKFmNgh4F5gETK5TZh5wPfAycCXwvLu7mRUR3GBVm9lgYCiw3t23m9ke\nMzsHKAWuA/4vRb9HMkRjNzdNXYNKspLqJ8kJsZO4/H3tVl59extd2uWz51BiM3i6B90/h5zQmXe2\n7qUpkyKvrNzNursvZVZpOc+s3MQlw3ur/q2HEj8RaToz+NDVsOA7sO0t6HlyZKG4e5WZ3QI8C+QD\nD7n7KjObDix193nATOA3ZrYO2E5w8wUwFphuZlVANXCTu9f2X70Z+BXQgWC2PM2YJ8e4/Q8row5B\n0pzqJ8kFM+avPm4Sl7lllRR3a59w4ldr3ft7AejWoYCdB6oSOqamJuhSOnlMfyV8jVDiJyLNM/wz\nQeL3xh/gY1+PNBR3n08wpXnstttjPh8Eropz3FPAU/WccykwPLmRSjapqml4/IoIqH6S7PfAwvVx\nt1fsPNjsc+48UMVFw07k1be3NZoAnje0V7Ovk2s0xk9EmqewGPqOhNV/jDoSERERyTLPvbGZnQeq\nGNa7S4PlHrlhTIoiynxK/ESk+U7/FFT+A3ZubLysSA4Z2sg04yIi2WDZhh1NKt+mGdMdv7FpD+3b\n5MWtV7WmX9Mo8ROR5jv908G7Wv1EjjHxw8WNFxIRyTDXzSzllNvmc/4PXmDZhh28sn5bk44/0sgS\nD/U5eKSGiR8u5p0Zl1FSXEhBnlFSXMjcW85r1vlylcb4iUjz9TwZThweJH4f+deooxFJG+cM7tl4\nIRGRDBI7c+c72/Zzxf2LKe7WHqOehSSTrLZeVbLXfGrxE5GWOfVS2PgKHGhadw+RbDZyQPeoQxAR\nSaq6M3dCMIFLqqa5Ur3ackr8RKRlhlwIXgPr/xZ1JCIiIiJSDyV+ItIyfUdBu0J4669RRyIiIiKt\nYMb81QmXLWiF7KIZc8JIHBrjJyItk18Ag8fCuufBPVjcXURERDLWdTNLeWX9Nvp068CPPlvCE8sq\nEj62qgbyDJq71GnPTm3Ytu/IMdueuOnc5p1MjpGUnNzMxpvZGjNbZ2bT4uz/upm9YWavm9lfzWxA\nzL5qMysLX/OSEY+IpNjJF8LuCtj6z6gjEUmJifcuijoEEZFWUTuJy+FqPzqJy679h5t0jtqkrzkt\ndQ9edzY3jR3MSV3bMXpgd566+VyN70uSFrf4mVk+cB8wDqgAlpjZPHd/I6bYcmCUu+83s5uB7wOf\nC/cdcPeSlsYhIhE6+RPB+/q/QdGp0cYikgJlFbuiDkFEpFXEm8SlmaswHD0u3xI/x8gB3Rk5oDvT\nLj29eReVeiWjxW80sM7d17v7YWA2MCG2gLu/4O77w6+vAFrgSCSbdOsPXfsGs3uKiIiIxKh20ECQ\n6CUj8esLbIz5XhFuq88NwDMx39ub2VIze8XMJiYhHhFJNTPoNwY2vByM8xMREZGMs2xD6y3NpLuD\n6CVjcpd4CXzc/2/N7PPAKODjMZv7u3ulmQ0GnjezFe7+VpxjpwBTAPr379/yqEUkufp/BFb9HnZt\nDFoARXJUxzaaMFtEMseyDTt4Zf02zhnck5//7bhb8FZhFv858dihvVJy/VyVjMSvAugX870YqKxb\nyMw+CdwGfNzdD9Vud/fK8H29mb0InAUc91+duz8IPAgwatQoPTQQSTf9zwney19R4ic57TdfPifq\nEEREEjLuRy+ydss+IJiJs3dh+5Rc1x0uGnYiC97YfLS1qKS4kEduGJOS6+eqZCR+S4ChZjYIeBeY\nBEyOLWBmZwE/B8a7+/sx27sD+939kJn1Aj5KMPGLiGSaE8+Atl2CxO/Mz0YdjUhkNPuciGSCifcu\nOpr0QTAT57s7D6bs+l/5+Ml85eMnH21tVN3Z+lqc+Ll7lZndAjwL5AMPufsqM5sOLHX3ecAPgM7A\n7yxY46vc3T8NnA783MxqCMYbzqgzG6iIZIq8fOj9Idj0WtSRiIiISCOinp24NtFTwpc6SVnA3d3n\nA/PrbLs95vMn6zluMTAiGTGISBrofSYsfRhqqoNEUERERETSgkagi0jynHQmVB2AbeuijkRERERi\nLNuwg/teWNeqM3dKektKi5+ICAAnhQ34m17XQu6StXTTJCKZZtmGHVx5/+KjE6lcMyZ1k7AZkJ8P\nVdUfbJtY0idl15cPqMVPRJKn6FTIbwfvaZyfZK9rZ5ZGHYKISJNc9cDiY9Zae7S0PGXXdmDdXZcx\nsaQP3Tq2YWJJH+6ZdFbKri8fUIufiCRPfhsoOgXefzPqSERazf7D1Y0XEhFJExPvXURNGiyEpmQv\nemrxE5Hk6jlUY/xERETSRNSzd2pR9vShxE9EkqvnENi5AaoORx2JiIiIREyLsqcPJX4iklw9h4DX\nwI53oo5EJOVKigujDkFEctzU2cspmf4cU2cvT+l1OxQcn1aoTkwvSvxEJLl6Dgne1d1TslBjN1Jz\nbzkvRZGIiBxv6uzlzC2rZOf+I8wtq+SM2/9MmxTd7R+oquGdGZdRUlxIQZ5RUlyoOjHNKPETkeTq\nOTh4T2HiZ2bjzWyNma0zs2lx9rczs8fD/aVmNjDcPs7MlpnZivD9EzHHvBiesyx8nZCyHyRpa25Z\nZdQhSIZR/SSpVLeO2ne4miM1KY7hlvNYd/elSvrSkGb1FJHkat8N2naG3am5QTazfOA+YBxQASwx\ns3nu/kZMsRuAHe4+xMwmAd8DPgdsBT7l7pVmNhx4Fugbc9w17r40JT9ERLKO6icRSSdq8ROR5DKD\nrn1gT8paRkYD69x9vbsfBmYDE+qUmQD8Ovz8JHChmZm7L3f32kBXAe3NrF1KohaRXKD6SUTShhI/\nEUm+Lr1T1uJH8AR8Y8z3Co59Kn5MGXevAnYBPeuUuQJY7u6HYrY9HHaj+o6ZWXLDFpEcoPpJWtWs\n0nKunVnKrNJylm3YEWksRZ3bRnp9aZy6eopI8nXtC28vTNXV4t3w1F2qtsEyZnYGQfeqi2L2X+Pu\n75pZF+Ap4FrgkeMubjYFmALQv3//pkUuItlO9ZO0mlml5dw6ZwUAf1+7lR6d2kQWS+e2+Sz59rjI\nri+JUYufiCRf1z6wZxPUVKfiahVAv5jvxUDd5sajZcysACgEtoffi4E5wHXu/lbtAe7+bvi+B5hF\n0GXrOO7+oLuPcvdRRUVFSflBIpI1VD9Jq7ktTPpqbd93JKXXv/vyEXxsaC/uvnwEK6ePT+m1pXmy\nvsVv4LSnj35+Z8ZlEUYikkM6nwBeDQd2Qqe6PZaSbgkw1MwGAe8Ck4DJdcrMA64HXgauBJ53dzez\nbsDTwLfc/aXawuHNVzd332pmbYB/Af7S2j9ERLKO6idpFTPmrz6u6TjVJo/pz+QxaknOJFnd4heb\n9MX7LiKtpH234P3gzla/VDgm5haCGe9WA0+4+yozm25mnw6LzQR6mtk64OtA7ZTqtwBDgO/UmRa9\nHfCsmb0OlBHcsP2i1X+MiGQV1U/SWn6+cH2k19d4vsyU9S1+IhKBDmHid6D1Ez8Ad58PzK+z7faY\nzweBq+IcdydwZz2nHZnMGEUkN6l+ktYQdWufxvNlJiV+IpJ8R1v8op1hTCSVJpb0iToEEclSyzbs\n4OuPl7Fp1wHOGdzqQyiO0SbfOFL9Qaqpui5zKfETkeRrXxi8p6jFTyQVJt67qMH990w6K0WRiEgu\nWbZhB1fcv/jo94Vrt2KkrtVvWO+uDOrViRf/uYXzTylSXZfBlPiJSPJ1SN0YP5FUKavYFXUIIpKD\nPvvA4uO2pbKr5+fO1iQu2SKrJ3cRkYgUtAveqw5HG4eIiEiGq454QJ+SvuyhxE9Eki8/TPyqD0Ub\nh4iIiIgASerqaWbjgf8F8oFfuvuMOvvbAY8QzEK1Dficu78T7vsWcANQDXzN3Z9NRkxTZy9PxmlE\npDnU4iciItIss0rL+fGCNew7VM3FZ5wYaSxd2uVHen1Jrha3+JlZPnAfcAkwDLjazIbVKXYDsMPd\nhwA/Ab4XHjuMYDHTM4DxwM/C87XY3LJK+ttm/q/NTznD3k7GKUUkUXn5YHlQrcRPREQkUbNKy7l1\nzgq27j3MgSPVzC2rTOn1xw07NtH81ZfGpPT60rqS0eI3Gljn7usBzGw2MAF4I6bMBOC/w89PAvea\nmYXbZ7v7IeDtcPHS0cDLSYiL7uzhU/mv8Pvqj7Eq6gVPRHJNfjt19RQREWmC2+auiPT6v7huFMs2\n7OCV9ds4Z3BPRg7oHmk8klzJGOPXF9gY870i3Ba3jLtXAbuAngke22wWvivnE4lAXj7U1EQdhUjS\nXZ73d55v+3W6sjfqUEQky3gLblqNGl5s++9cm/9ci2IYOaA7X71giJK+LJSMxM/ibKv7n219ZRI5\nNjiB2RQzW2pmS7ds2dLEEEUk5WqqguRPJMt0s70Mznsv7h8wEZGodGU/A/M205aqqEORNJWMxK8C\n6BfzvRio2yH5aBkzKwAKge0JHguAuz/o7qPcfVRRUVFCgVWFP68AtTqIpFxNNeRpqVDJPnnh80mP\nSf0KNEe2iDTTsg07uO+FdSzbsKNF5+lpuwHY6l2bdbzqseyXjLuyJcBQMxsEvEswWcvkOmXmAdcT\njN27Enje3d3M5gGzzOzHQB9gKPBqEmIC4DBtAGjLEQDemXFZsk4tIo2pqVLiJ1kj9obMwsSvOubZ\n6fQJI1Iek4hkvmUbdnDl/YuTMiypB0Hit53EEr+u7QvYffCD1kHVY9mvxXdl7l5lZrcAzxIs5/CQ\nu68ys+nAUnefB8wEfhNO3rKdIDkkLPcEwUQwVcBX3b26pTHVOhL+PDV5i6RYTQ3g6uopWeOLD3/w\nTDIv7EVSE9PipwWORaQ5rrh/cdLO1dP2ALA9wRa/MYN7csGpJ/DMyk1cMry36rEckJTH8e4+H5hf\nZ9vtMZ8PAlfVc+xdwF3JiKOuwx78vDamxE8kpWrCf3NK/CRLxD4Vr+3qWZOU0RIikqtOuW1+44Wa\noKldPW/6+MmMHNBdCV8Oyeq/WrVdPduFXT1FJEUOh7Mdtu0cbRwirSDeGD8RkaY6XJ3ceedru3ru\noEtC5TVrZ+7J2sTv7stHcIC2AHTgEHdfrn7LIilzKOhuQrvmDTAXSWfxunqKiCQiWRO5xNPTdrPb\nOxxt+BCpK2tnXpg8pj94DTV/zmPCaZ0YrmZskdQ5mvgl9tRRJJOoq6eINEcyJ3KJp4ftSXh8X0lx\nYStFIeksq/9qTT5nIHkdChneQ0u4i6SUEj/JYnmmFj8RabovPvxqqyV9AD3ZxbYGZvQsKS6kIM8o\nKS5k7i3ntWIkkq6ytsXvqPbd4MDOqKMQyS3q6ilZzHBq3ECJn4g0QewkUa2hp+2hwutf61rJnmR1\nix8AHbrBQSV+Iil1YHvw3l5dSST75OFq7RORtNPDdrPN1dNG6pcDiV93OJD8AbQi0oC9m4P3LidG\nG4dIK8ijRomfiKQZpzt7El68XXJT9id+nU+CPZujjkIkt+zZDG06aYyfZKU8HI/58zm0qFOE0YhI\nuppVWs61M0uZVVre6tfqyn7aWjXbvCvdOhw/kkv1lEAujPHr2hv2vgc1NZCX/XmuSFrYuxk6nxB1\nFCJJUfemrW6L34JvnJ/iiEQk3c0qLefWOSsA+PvarTz00tvkGyR56b6jahdv3+Zd2XmgindmXMa4\nH73IW1v3cXKvTqqnBMiFxK9Lb6ipgn1b1O1MJFX2bobO+vcm2eH2P6w85rvG+IlIY+rWG+ve39uq\n1+tOMKna9pjF25XsSV3Z3wTWtU/wvvvdaOMQySV73tODFskaVTXHP6J3JX4i0oB49UZrakM1AEdy\noE1Hmi/7E78uvYP3PZuijUMkV9TUwM5y6NY/pZc1s/FmtsbM1pnZtDj725nZ4+H+UjMbGLPvW+H2\nNWZ2caLnlNxlrboal2QT1U2SClXhLX1BmACKxJP9iV+3AcH7jg3RxiGSK/ZUQvUh6D4oZZc0s3zg\nPuASYBhwtZkNq1PsBmCHuw8BfgJ8Lzx2GDAJOAMYD/zMzPITPKfkILX2SaJUN0mqVJEPQL4SP2lA\n9rcHd+wRrCW2bV3UkYjkhu1vB+89Upf4AaOBde6+HsDMZgMTgDdiykwA/jv8/CRwr5lZuH22ux8C\n3jazdeH5SOCcIiINUd2UQybeu4iVlbsZ3if1SyrUJn4F1DCxpE/Kry+ZIftb/Myg5xDY/lbUkYjk\nhh1h4pfCFj+gL7Ax5ntFuC1uGXevAnYBPRs4NpFzSg7oUHD8n0q1+UmCVDfliIn3LqKsYhdVNU5Z\nxa6UX3/sqScBcO6gQu6ZdFbKry+ZIfsTP4AeJ8M2JX4iKbHjHcgrgMJ+qbxqvPvwuoOw6ivT1O3H\nntRsipktNbOlW7ZsaTRQyTwHqmqO+a7RfdIEkdVNoPoplaJI9mJ987IRAHzpIyn92ysZJjcSv55D\nYFcFHDkQdSQi2W/rP6H7QMhPaU/yCiD2r10xUFlfGTMrAAqB7Q0cm8g5cfcH3X2Uu48qKipq4c+Q\nTKHJXSRBkdVNoPopp+SFf3NrqqKNQ9JabiR+vYYADlvXRh2JSPbbvApOPCPVV10CDDWzQWbWlmBC\nhHl1yswDrg8/Xwk87+4ebp8Uzqw3CBgKvJrgOSUHaXIXaQLVTZIaecEYPyV+0pDsn9wF4MSg+ZvN\nK6H3mdHGIpLNDu0NJnf50NUpvay7V5nZLcCzQD7wkLuvMrPpwFJ3nwfMBH4TTpCwneBmibDcEwQT\nI1QBX3X3aoB450zpD5O0VdviN3Zor4gjkXSmuim7Lduwg1fWb+OcwT0jjWPs0F5q8ZOE5Ebi1/Nk\nKOgA762IOhKR7LblTcCjaPHD3ecD8+tsuz3m80HgqnqOvQu4K5FzisQOsXrkhjERxiGZQHVTdlq2\nYQdX3r/4aKfvru0L2H0wNUlXnsF5Q3rx6jvbGT2wR1AP7dkc7FTiJw3IjcQvLz+4EVXiJ9K6NocP\nnSNI/ERSSZ09RXLbFx9+9ZiRvqlK+gBqPM5Dp9oWv2olflK/3BjjB3DS8CDxcw3IF2k1762Atp2h\nsH/UkYi0Gv0VEZFUJnoJ0Rg/SUCLEj8z62FmC8xsbfjePU6ZEjN72cxWmdnrZva5mH2/MrO3zaws\nfJW0JJ4GnTQCDu6EneWtdgmRnFexBPqcBXm580xJcpXSPxFJIxrjJwlo6d3ZNOCv7j4U+Gv4va79\nwHXufgYwHrjHzLrF7P9Pdy8JX2UtjKd+xWcH7xVLWu0SIjnt8P5gAqV+o6OORKRVaVZPEYlSfrwq\nKL9N8K7ETxrQ0sRvAvDr8POvgYl1C7j7P919bfi5EngfSP1iMiecEXRBK38l5ZcWyQmbyoI/OLUP\nWUSymNbxE8ktM+av5py7/8JnH1jMsg07UnrtIUWdjvl+x8QRxxc6OsbvSAoikkzV0sldTnT3TQDu\nvsnMTmiosJmNBtoCb8VsvsvMbidsMXT3Qy2MKb78Aug7EjaWtsrpRXLexleDdyV+kkWmzl5+3LYa\njDwlfiI5Y8b81TywcD0A7+0+xBX3L6Zbxzbs3J+aJKtPtw586bzBPLNyE5cM783kMXHG0eflQ/tu\nsG9LSmKSzNRo4mdmfwFOirPrtqZcyMx6A78Brnf3mnDzt4D3CJLBB4FvAtPrOX4KMAWgf/9mThzR\n/xxY+AM4tAfadWneOUQkvool0H0QdNK6ZpI95pZVHretmjwKrCZOaRHJRrNePX5+iFQlfQAL127l\nkRvGxE/4YhX2g10VqQlKMlKjiZ+7f7K+fWa22cx6h619vQm6ccYr1xV4Gvi2ux/ta1nbWggcMrOH\ngf9oII4HCZJDRo0a1bxHrf3GgNcEN6gnf6JZpxCROGqq4Z1FcNplUUci0upqPBglYSj5E8kFnikz\nwhf2VeInDWrpGL95wPXh5+uBP9QtYGZtgTnAI+7+uzr7eofvRjA+cGUL42lYvzGQ1wbW/61VLyOS\nc957PZg1d/D5UUci0uqqwz+d+Ur8RHLCNWMGRHr9hG/WC4thtxI/qV9LE78ZwDgzWwuMC79jZqPM\n7Jdhmc8CY4EvxFm24VEzWwGsAHoBd7Ywnoa16xzMOLj+hVa9jEjOWf9i8D7o45GGIZIKNeGfzoln\nnhhxJCKSCtMuPT2l1yvq3PaY71PGDk7swMJiOLADDu1thagkG7Rochd33wZcGGf7UuDL4effAr+t\n5/jU97c8+QJ4/i7Yt1VjkUSSZf3foOh06KIbYcl+tS1+C1YdP/5PRLLHxHsXsbJyN8P7dGXs0F4s\nXLs1Jdfdc7CKm8YO5s+r3mP8GSclnngW9gved78LRae2XoCSsXJvleXBnwD8gxYKEWmZIweh/GV1\n85Ss1KHg+D+TtYlfdXV1qsMRkRSZeO8iyip2UVXjlFXsYuHarXRum5+Sax+sqmHapafz4n9e0LTW\nxq59g/ddG1snMMl4uZf49SkJprt9S909RZJiwyKoOhi0potkmQNVx4/jq+3qmacxfiJZq6xi13Hb\n9h6uxiBlCWCTFRYH75rgReqRe4lfXj4MuRDWPhvMRCgiLfPmfGjTEQaNjToSkZTQ5C4iucsJEsCC\nPKNbh5Yuh51kXXqD5cGud6OORNJU7iV+EEw5v2/LBwtOi0jzuMOaZ4KHKW06RB2NSErUJn6F7Szi\nSEQkKlU1zs4DVbTJT34C2KVdM1sU8wugSx+1+Em9cjPxGzIO8tvCm3+KOhKRzFa5HPZUwqlav09y\nR21Xz59OKmmkpIhkklml5Vw7s5RZpeXkJ/hc50h1yxPAuoner740plnnAcK1/DTGT+JLszbqFGnf\nNZh2/s0/wUV3gumprUizrJkPlg+nXBx1JCJJN6u0PO72GoK/GWf26ZLKcESkFc0qLefWOSsA+Pva\nrQzr3YU3Nu1J+PjaBLBDQV7cscEN2XOomqduPpdX1m/jnME9GTmge5OOP0ZhMbz7j+YfL1ktN1v8\nIOjuueMd2Lwq6khEMpM7rP4j9P8IdOwRdTQiSXfnn+L/fajt6olrnLhItqj777026evavmltJLFJ\nX7uCxBsWRg7ozlcvGNKypA/CRdwroUZjkOV4OZz4/UvQUrHyqagjEclM762ALW/C8M9EHYlIq9h/\nJP6NU42Hfzo1QZhI1qjv3/vug1UAdGzGTJ6HqhxoevLYIoX9oPoQ7E/NmoOSWXI38etcBCd/Alb8\nTk9FRJpjxROQVwBnXB51JCIpdbTFT4mfSM7Yfzj4957XjNFBtcljfQlgnOVCm692SYed8buqS27L\n3cQP4MzPBQNgy1+OOhKRzFJTDSueCiZKUjdPyVL1/YGsInzyX1OVslhEpHUlekNc482/Rm0CePpJ\nx44Pnj5hRPNPWtcJw4L3TWXJO6dkjdxO/E67FNp0gtcfjzoSkcyy4aVgNs8zr4o6EpFWM2Xs4Ljb\nj9TOi1Z9OIXRiEiyTZ29nJLpzzF19vJ6/703pluHAgqa2At0zeY93H35CD42tBd3Xz6CyWP6N+va\n8QPqD52KoGJp8s4pWSM3Z/Ws1bYTnP4pWDUXLvk+tGkfdUQimeG1x6FtZzjlkqgjEWk17+0+GHf7\n4do/nTVHUhiNiCTT1NnLmVtWCcDcsko6t8tnQI+OvLf7wNGxeYnYeSBoxevWoYA9h6qoTmD0UG2r\n4W9uaMGyDfUxg+KzlfhJXLnd4gfwoUlwaBesnhd1JCKZ4cBOWPV7GH4FtO0YdTQirab2prCuo109\nq5X4iWSqea8d++9776FqNmzfz6Eqp1vHNnRok5fwWn4QJIDVNdA5wUlgvv/sm00Jt2n6joRta+HA\njta7hmQkJX6DPg49BsPSh6OORCQzvP4EHNkPo74YdSQikfigq6cSP5FM1dBYvZ37j3DgSA2DijpT\nUlzYpK6cew9Xk2+NT9iy92Ar1h/Fo4J3recndSjxy8uDkV+A8sXw/uqooxFJb+6w7GHoXQJ9zoo6\nGsysh5ktMLO14XvcBZDM7PqwzFozuz7c1tHMnjazN81slZnNiCn/BTPbYmZl4evLqfpNkv6OeG2L\nn8b4Sf1UP6W3oUWdGi2z7v29lFXsYmDPIAHMT/Cuudqhdjm/+hoNzz25V2Ina44+Hw6urO6eUocS\nP4CSayC/rVr9RBqz8VV4/w0Y9aWoI6k1Dfiruw8F/hp+P4aZ9QC+C4wBRgPfjbkB+6G7nwacBXzU\nzGIHLT7u7iXh65et+iskoxxt8dOsntIw1U9pbME3zq83KaurNgEc1KtzQgljrNqGxbFDezGgR0fa\n5htjh/bikdYY31erfVcoOg3eVeInx8rtyV1qdeoFwybAa7Phk98NJn0RkeMtnQltuwTj+9LDBOD8\n8POvgReBb9YpczGwwN23A5jZAmC8uz8GvADg7ofN7B9AcQpilgzRsU1e3EWdPxjjpxY/aZDqpzQ0\ndfZyXvznFs4/pYgnbz6XK+5fnPCx697fC0BR57YU5Bmbdh9K+Ni/r93K2zMua3K8zVY8Et6cH/TU\nsWYsPihZSS1+tc6+MZjkZfmjUUcikp52vQsrn4KzroF2naOOptaJ7r4JIHw/IU6ZvsDGmO8V4baj\nzKwb8CmCp/K1rjCz183sSTPrl9ywJROMGhh/jcrDGuMniVH9lGZqZ/Lcuf8Ic8sqmfTgy4wd2qvJ\nyyls2Xv4aNLXu2u7hBZgb8Hyf83TdxQc2A7b16f6ypLGlPjV6j8GikfDy/cGi1OLyLFe/Tl4DZzz\nrym9rJn9xcxWxnlNSPQUcbYd/RtsZgXAY8BP3b32L+QfgYHufibwF4Kn9fXFN8XMlprZ0i1btiQY\nkmSCv6/dGne7ZvWUWqqfMsvTKzYd8/1ItbNw7VaeXLqR00/qwtkD4w7DbNCm3YeoqoFrxvRvcjfQ\nVlV8dvD+7rJo45C0omTgl4UAACAASURBVMTv/2/vvsOkqu4/jr+/21h6F+lFUFERkEUEFUUFFYli\nBxVBIYbEEqNJJCY/k9iCxliSmBjsGkVjxd4QBQtKCQIBFURQmlIEBGRh4fz+OHd1WGbZWabcKZ/X\n89xndu7cO/d7duaeueeeFqnPpbBuCcx/PuxIRNLLlg2+D+wBg6Fh25Qe2jl3nHPuoCjLBOArM2sO\nEDx+HeUtlgKRd8RbAZHjeI8DFjjnbo845hrnXHkbnruBHruJb5xzrsQ5V9K0adM9S6Skpcru0G9z\nmsBdPOVPmWX79uhn9dbtjvkrv2Xa4m+oV1xA/ZrV7wn1yAdfsHDVJm48tUvU1xvVKqz2e8Zlr85Q\nWFsDvMhOVPCLtP9JfmqH9/7q20SLiDfzISjd4G+OpJfngOHB38OBCVG2eRUYYGYNg0ETBgTrMLPr\ngfrA5ZE7lF+sBU4GNOSvfG+bJnCX2Ch/SjON6hRVuc2GLWWs/66M4oI8WjYoplur+jG/vwOufmYO\neexa0Lt7eM9qRhunvHw/+rYGeJEIKvhFysuH3hf7avEl74YdjUh6KNsKU/8JbY+AloeEHU1FY4H+\nZrYA6B88x8xKzOwegGDQhOuAacFyrXNurZm1An4LHADMrDAs+mXBEOofAZcBI1KZKElv29TUU2Kj\n/ClNPPrBFwy79wOO6Bj7FApbynawbN0WFq3eRN9OTWjZoDjmfXcAazdvo1PT2vzq+P146qd96NG2\n+s1I49aqB6yYDdu2pP7YkpbiGtUzGIb4caAdsBg4yzn3TZTttgNzgqdfOOdODta3Bx4DGgEzgWHO\nuXDbznQ7F96+Gd4aCyNeCDUUkbTw34dhw1I4+Y6wI9mFc24NcGyU9dOBURHP7wPuq7DNUiqZYsk5\n9xvgNwkNVjLKox98UelrmsBdYqH8KT08+sEXXP3MnO+fH9C8LoX5eWwuLWPBqk1V7r9hSxmTF6ym\nce1CmtYpYt13W9kW41AQC1Zt4uJ+Hfc09Pi16ulbJqycA61TXOMoaSneGr8q56gJfBcx38zJEetv\nAm4L9v8GGBlnPPErrAlH/AIWT4HPp4QdjUi4ykphyl+gdS/YZ5frF5Gsdfsbn1T62g8FP/XxE0l3\nFc/leSu+ZfbS9RzaoTGj+3agY4wDsqzZtI1VG32hr26NfIryM2CKhJYl/vHLqeHGIWkj3oLfKfww\nmtSDwOBYdzQzA44BntyT/ZOqxwiosze89Sf19ZPcNvMh2LAMjv6N5gGSnLJ6Y+WFuq0q+IlkjGjn\nssMPxjJuyiIuPKIDAw5oRpMY+v+V+7Z0O1u3O1o1KKZto1oUpmshsF5zaNoZPnkl7EgkTcRb8Itl\njhqA4mA44almVl64awysc86VBc93mbsmUkqHJC6sCUde4fv5fT45uccSSVfbtvjavjZ9oMPRYUcj\nklJ5u7mO69i0LuQXQZn6zYiku92dyzucH4zltXlfsXbTVgZ3a1Gt0TeXrtvCkrWbadu4Nt1a1Se/\nwlV1306x9ylMms6D4Iv3YFP06Wkkt1RZ8EvAHDUAbZxzJcA5wO1mtg9VzF2zywupHpL4kOFQtwVM\n/KNq/SQ3Tb8Xvl0B/VTbJ7ln0MEtKn1txYYtUFDsm0KLSFrb3bkcaYeDCbOWc1DL+tx4ahe6tqof\ncyFw4dcbmbV0PV1a1KdvpyYUF+bRt1MTHhrZK57QE2P/QX4O3k9eDjsSSQNVFvwSMEcNzrnlweMi\n4C2gO7AaaBBMTgq7zl0TrsJiOOZ3foTP/z0ddjQiqbV5rR/kaJ9joH3fsKMRSbnbh3Sv9LWNpduh\noIZq/EQywO1DujO4WwsKYmjj5oDJC1ZzzYS5tG9Sm5nXDGB03w7ULc6P6Vizlq6nUe0iPr7uxPQo\n9AE07wr128DHL4YdiaSBeJt6VjlHTTA3TY3g7ybA4cA855wDJgFn7G7/UHUdAnt3gTf+oDu7klve\nvtnP2zfg+rAjEUlPqvETyRi3D+nOwhtPirnpZdkOx7OzlrP/717m2VnLOPfQtozu24FaRVVfNr84\nZ0W84SaWmZ+n+rM3oXRj2NFIyOIt+FU5Rw3QGZgezDczCRjrnJsXvHYVcIWZLcT3+bs3zngSKy/f\nX/iu+wI++FfY0YikxuqFMO1u6D4Mmh0YdjQi6Uk1fiIZ56GRvXjqp304p1eb6HNlVLClbAcrN5Ry\n1+RF3P/eYgYcsDeDu7WgQa1CalZShejSsXtQ50GwvRQWvhF2JBKyuObxi2WOGufce0CXSvZfBBwa\nTwxJ1+Fo6DQAJt8C3c6B2mnQUVckmd74va/N6PfbsCMRCVXbRrVYsnbzLusHd2sB36jGTyQT9Wjb\nkB5tG3L6Ia14auZSxn/wReUDTEQoLdvBs7OWU7son2GHtWXMwM5c/th/eWHOcsoi5vWLtU9hSrU+\nDGo2go9fgAPTYwB9CUe8NX65YcANsG0zvH5N2JGIJNfCif6H4YjLoW6zsKMRCdWtZ3eLuv72Id2D\npp6q8RPJVD3aNuTGU7vw5E/78Kvj96Nbq/ox7bdp63bumryInte/zqwv1zHq8A7f1wIO7tZit/2D\nQ5NfAPsNhE9fgzJNQ5PLVPCLRdN9oc8lMOsRWPJe2NGIJMe27+DFK6HRPtD70rCjEQldj7YNo64f\n+9J89fETyRI92jbk4n4defaSI3jqp33o2S76eV/Rqo1bWbxmM3dNXsTe9YqZdc2A9Cz0les8CErX\nw+IpYUciIVLBL1Z9f+1HRXrxSti+LexoRBJvyl/gm89h0G1+VFsRierhqUt8H79t34UdiogkUI+2\nDXlidB9uPLULHZvWjnm/h6cuSWJUCdKhHxTW9q16JGep4Berolpw4k3w9TyY+o+woxFJrFWfwju3\nw8FnQ4ejwo5GJG0U5u86BMR327arxk8ki53Tqw1vXHk0N57ahSM7NaFVg93fDP1u2/bdvp4WCouh\n03Hw8UuwY0fY0UhIVPCrjv0Hwr4nwqQ/wZrPwo5GJDF27IAXr/A3NzR9g8hORh7efpd1xQV5GtVT\nJAec06sND4/sxTtjjuXGU7vQpE5R1O2KY5kkMB3sPwg2roRl08OOREKSId/UNDLoVsgvgmd/Bjsy\n4A6PSFWm3ePb/Pe/DursFXY0ImllzMDOuwz6cH7vdqrxE8kx5/Rqw/Tf9Y86ofv5vduFE1R1dRoA\neQUw//mwI5GQqOBXXfVawIlj4cup8MFdYUcjEp81n/nRajv2h0PODzsakbT07CVHMLpvB9o1rsXo\nvh0YM7CzavxEctSYgZ2Z84cTds0TMkHNBtC+r+/nl47zDUrSxTWPX87qOhTmTYCJ10Kn46FJx7Aj\nEqm+HdvhmdH+Avbkv4HFMp2tSG4aM7Dzzhd3ms5BJKftkidkiv0H+e4dK2dD865hRyMpphq/PWEG\ng273P/zPXKRRPiUzvfdXWPohDLwF6jUPOxqRzFJYU6N6ikjmOeg0KKgJ0+8LOxIJgQp+e6pec/jR\nHbBsBrx5XdjRiFTPspnw5g1wwCnQ5YywoxHJPEW1YMc23fgTkcxSsyF0OR1m/we2rA87GkkxFfzi\nceBg6DEC3r0DFk4MOxqR2GxZD09eAHWa+ZprNfEUqb7CWv5x66Zw4xARqa6eo2DbZpg1PuxIJMVU\n8IvX8X+Cpp19X6mNX4cdjcjuOQfPXQrrl8KZ90OtRmFHJJKZygt+au4pIpmmRXdoWeJH9dYgLzlF\nBb94FdWCM+6D0g3w1EjYXhZ2RCKVm3aPH5jo2Gug9aFhRyOSuYpq+8dtm8ONQ0RkT/QcBWsWwOeT\nw45EUkgFv0RodgCcdKs/eSb+IexoRKJbNgNevdqPRNv70rCjEclshTX9o5p6ikgmOvBUqNkIpt0d\ndiSSQir4JUr3c6Hnj+G9v8GcJ8OORmRnG1bAY+dC3b3h1LsgT6e+SFzU1FNEMllhMRwyDD5+CdYv\nCzsaSRFd/SXS8TdCm94w4RJYOSfsaES8bd/B4+fClg0w9DH16xNJhO+beqrGT0QyVI8LwO2AmQ+G\nHYmkiAp+iVRQBGc+CDUbwPhz4Nuvwo5Icp1z8NxlvpnnaeOg2YFhRySSHb5v6qk+fiKSoRq1h079\nYcYDULY17GgkBVTwS7S6zWDIo7B5NYw/W/0/JFzv3AZz/gPH/A46Dwo7mqQws0Zm9rqZLQgeG1ay\n3fBgmwVmNjxi/Vtm9omZzQqWvYL1NczscTNbaGYfmFm71KRIMkJheY2fmnpKdMqbJCP0/DFs/Ao+\nfiHsSCQFVPBLhpaH+JE+V3wET46EHdvDjkhy0UePwcQ/wkFnwJG/DDuaZBoDTHTOdQImBs93YmaN\ngN8DvYBDgd9XuAg71znXLVjK52UZCXzjnOsI3AbclMxESIYpKu/jp5t7UinlTZL+Oh4LDdr6Ub8l\n66nglyz7nQgn3ASfvgyvjNE8KZJaC9+ACRdD+74w+B/ZPkn7KUB5B4UHgcFRtjkeeN05t9Y59w3w\nOnBCNd73SeBYs+z+R0o1qKmnVE15k6S/vHwouRCWvAtfzQs7GkkyFfySqddF0PsS+HAcTP5z2NFI\nrlg2Ax4/H5p2hrMfgYIaYUeUbM2ccysAgse9omzTEvgy4vnSYF25+4OmVP8XcQH1/T7OuTJgPdA4\n0cFLhirUPH5SJeVNkhm6D4P8GjD93rAjkSSLq+AXS/t1M+sX0T59lpltMbPBwWsPmNnnEa91iyee\ntNT/Oug6FCbdAFP/GXY0ku1WL4RHzoLajeG8J6G4XtgRJYSZvWFmc6Msp8T6FlHWlVfDn+uc6wIc\nGSzDYtgnMraLzGy6mU1ftWpVjOFIxisogrwCFfxyXDrnTUF8yp+karUbw0Gn+S4iWzaEHY0kUbw1\nflW2X3fOTSpvnw4cA2wGXovY5FcR7ddnxRlP+snLg5P/Dp1/5Jt8znw47IgkW635DB4MBnA57xk/\nZ1+WcM4d55w7KMoyAfjKzJoDBI9fR3mLpUDriOetgOXBey8LHr8FHsX3s9lpHzMrAOoDa6PENs45\nV+KcK2natGkikiuZorCWmnrmuHTOm4J9lT9JbHqOgq0bYeZDYUciSRRvwS+W9uuRzgBeds7l1i9l\nfgGcfi/scyw8f5kmeJfEW/MZPDAItm+F4c9Dk45hR5RKzwHlI+ENByZE2eZVYICZNQxaJgwAXjWz\nAjNrAmBmhcAgYG6U9z0DeNM5ddaVCEV1YOu3YUch6Ut5k2SOViXQ/ih493Yo3Rh2NJIk8Rb8Ymm/\nHmkIML7CuhvMbLaZ3WZm2dsZqaAGnP1vaNMHnv4xzHo07IgkW6z9HB78EZRtgfOfg2YHhB1Rqo0F\n+pvZAqB/8BwzKzGzewCcc2uB64BpwXJtsK4G/iJrNjALWAbcHbzvvUBjM1sIXEGUFg2S42rU0QWS\n7I7yJsksx/wONq3yY1NIVrKqbhKZ2RtAtDZjvwUedM41iNj2G+dcZfPUNAdmAy2cc9si1q0EioBx\nwGfOuWsr2f8i4CKANm3a9FiyZEkVSUtTWzfDY0Nh0dvwozugx/Cq9xGpzJrP4KFTfPOM85+D5geH\nHdEeM7MZzrmSsOOIR0lJiZs+fXrYYUiqjOsHtRrBeU+FHYkkUTbkTaD8SWL0yFnw5Qdw+Wworh92\nNFKF6uZPVdb4JaD9ermzgGfKC33Be69wXilwPz+0X48WR3a0Uy+qBUMfh47H+WafH95d9T4i0ayY\nDfcdD1s3wbBnM7rQJ5KRatRVjZ+IZJd+V8OWdfD+P8KORJIg3qaesbRfLzeUCs08IwqNhu8fODfK\nftmnsBiGPAL7nQQv/RIm/Unz/En1LH4XHjgJ8ovgwlehRfYNiCuS9mrUhVL18RORLNKimx+Q8P07\nYXPUMYMkg8Vb8Kuy/XrwvB1+BKq3K+z/iJnNAeYATYDr44wncxTUgLMehG7nwdtjfe3f9rKwo5JM\n8MnL8O/ToE4zX+hrum/YEYnkJg3uIiLZ6OirfReSd+8IOxJJsIJ4dnbOrQGOjbJ+OjAq4vlidp6Q\ntHz9MfEcP+PlF8Ipf/fD7k+5BTZ+DWfc75uDikQz7R546de+Wee5T0LtJmFHJJK7NLiLiGSjZgdA\nlzP8IC+9L4Y6VY3dKJki3ho/iZcZHPt/cNJfYMFr8MBAWL8s7Kgk3Wwvg5d+BS9e6fuHDn9ehT6R\nsBXV8XfFRUSyzVFjoKwU3rkt7EgkgVTwSxc9R8GQR2H1Ari7H3w5LeyIJF18tw4ePTO483YJDB3v\n+xaJSLhq1PVzZ5aVhh2JiEhiNekI3YbCtHtVIZFFVPBLJ/udCKPegMJavubvv4+EHZGEbfUCuHcA\nfD4ZTv4bHH8D5OWHHZWIwA83YNTcU0SyUd9fg9vhuyNJVlDBL93s1Rl+/Ca06Q0Tfuab9m3bEnZU\nEoa5T8G4o/1kqsOehUPODzsiEYlUVMc/aoAXEclGDdv6+aZnPgTfLA47GkkAFfzSUa1GcN7T0OdS\nP5jHvf39RN2SG8pKfYH/yQuh2YEw+h1of2TYUYlIRTWCgp9q/EQkWx15JVg+vH1z2JFIAqjgl67y\nC2DA9TD0MVj3BfzrKPjfM2FHJcm29nM/Kfu0e3x/vhEvQv1dBsQVkXTwfY2fCn4ikqXqtfDjUHw0\n3nc/kYymgl+62+9EGD0F9tofnhgBz14MWzaEHZUkmnMw4wG46whYs8gP9HP8DX7KDxFJTzXq+UfV\n+IlINjviF1BQDG/9KexIJE4q+GWCBm3ggpd9dftHj8I/+8Cit8OOShLl25Xw6Fnw/M+h5SHw03dh\n/5PCjkpEqvJ9U0/djBORLFanKfQa7cceWDo97GgkDir4ZYr8Qjj2GrjwNcgvgodO9hN5b90UdmSy\np5zzmeg/DvOjdp54MwybAA1ahx2ZiMSiuL5/VMFPRLLdkVdA3ebwwuV+bmHJSCr4ZZrWPf1gH4f+\nBD78F9zZCz55OeyopLrWfg6PnOkHcGnYHn4yBXr9BPJ0SopkjPKC35b14cYhIpJsNerCCWNh5Rw/\nr7BkJF1lZqKiWjDwZrjgFT+4wPgh8Ni5mmAzE5Rthcm3+Fq+L96H42+Eka9D033DjkxEqquwlh/t\nTgU/EckFB5wCHY+DSTfAhuVhRyN7QAW/TNa2N/xksm8CuvANuPNQmHKr5v1LV59N8oO3vHkd7Hs8\nXPwh9L7Yj+AqIpnHzNf6qeAnIrnADAb+GXaUwStjwo5G9oAKfpmuoMgP+vKzqdC+L0z8I/y9J8x5\n0vchk/CtnAv/Ph0eHgxlW+CcJ+CshzRNg0g2KK6vkZZFJHc06gBH/hLmTYAFr4cdjVSTCn7ZolF7\nGDoezn8OataHp0b6id8XvxN2ZLlrw3I//cZdR8DSaX5exkumwb4Dwo5MRBJFNX4ikmsOvwwad4IX\nr4Rt34UdjVSDCn7ZpsNRcNHbcMqdsO5LeOAkePBHsOT9sCPLHRtWwCtXw18PgTn/8c05L5sFfS6F\nghphRyciiVRcTwU/EcktBTVg0K2wbokft0AyhjoXZaO8fOh+Hhx0Oky/H965De4/ATr0g76/grZ9\nfDttSax1X8C7d8DMh33794PPgqPHQMN2YUcmIslSXB9WLww7ChGR1GrfFw4+21/3HHy2BqnLEKrx\ny2aFNaH3z+DnH/lmhivnwAMD4Z5jYe7TmoclUb76H0y4GP7aHWY8CN2GwqUz4NS7VOgTyXbF9TWP\nn4jkpgHX+5HmX7xC40pkCNX45YKiWr6ZYclI+Gg8vH8nPHkBNGgDvUZD16FQq1HYUWaW7WXw8Qvw\n4d2w5B0oKPb/38Mvg/qtwo5ORFKluIGaeopIbqqzFxz7e1/wm/0f6Hp22BFJFVTwyyVFtaDnSOgx\nwk/6/t7f4NWr4Y0/woGDoccF0OYwNQPdnfXLYNajMON+2LDMF577Xwvdh6nwLJKLatSDrRv9zSBN\nzSIiuabHBf666NWr/eB1NRuGHZHshn6lclFePnQe5JeVc3zzxNmP+6XJvr5v2kFn+JFCBbZugvkv\nwEePwqK3AQcdjoaBt/j5+PLyQw5QREJTXN8/lm7QzR8RyT15eX6gl3FHw8RrYdBtYUcku6GCX67b\nuwucdAv0/yP87xn477/hzev90rIEupzpC4i51nyxdCMsfN3PU/Ppa7BtEzRoC0dd5ZsyNOoQdoQi\nkg7KC35b1qngJyK5qXlX33Vo6j+h6znQumfYEUklNLiLeEW1/UigF74Cl8/1zRe3l8IrV8FtB/q5\n6N68AZbNgB07wo42OdZ9CTMegMfOhT/vA0+M8PMgHnwmjHjJT8nQ7zcq9KURM2tkZq+b2YLgMWob\nEzMbHmyzwMyGB+vqmtmsiGW1md0evDbCzFZFvDYqlemSDPJ9wU8DvMjOlD9JTul3NdTdG174hQYP\nTGNx1fiZ2ZnAH4DOwKHOuemVbHcCcAeQD9zjnBsbrG8PPAY0AmYCw5xzW+OJSRKgQWs4/Od+WfUp\nfPqy7xM45RaYfDPUagLtjoD2R0K7vtCkU2b2C9y0Gr78wBfuFk6E1Z/49fVawSHD4YCToU1vNeVM\nb2OAic65sWY2Jnh+VeQGZtYI+D1QAjhghpk955z7BugWsd0M4OmIXR93zl2S7ARIhous8RPZmfIn\nyR016sIJY+GJ4TDlL3D0VVXvIykXb1PPucBpwL8q28DM8oE7gf7AUmBakKnNA24CbnPOPWZmdwEj\ngX/GGZMkUtN9/XL4z2HTGt/88bNJsHgKzHvWb1OnGbTsAS0OgZbd/WO6NXna9h18PQ9WzPa1ll9M\nhTUL/Gv5NaDd4dBjOHQ8zvdzzMSCbG46BTg6+PtB4C0qXFgBxwOvO+fWApjZ68AJwPjyDcysE7AX\nMCW54UrWKR/I4Ltvwo1D0pHyJ8ktBw6Gj8+Et2+CffpB60PDjkgqiKvg55ybD2C7v0g+FFjonFsU\nbPsYcIqZzQeOAc4JtnsQX3uogl+6qt0Yug7xi3OwdpEvAC55D5bNhE9e+mHbei19Aarp/r7g2LiT\nr0ms2wIKipITn3P+4mvtIljzGaxZCGs/g6/mwepPwW332xU38KOXdj/X1+g17waFxcmJSZKtmXNu\nBYBzboWZ7RVlm5bAlxHPlwbrIg3F30GPnIjodDPrC3wK/MI59yUiFangJ5VT/iS556S/+NZUT42C\n0e9Acb2wI5IIqRjcJVqm1gtoDKxzzpVFrK+Y2Um6MoPG+/ilxwi/bst6WD4Lls+Er+fDqk9g5kN+\nYJQfdvQ1hPVbQu2mULORv3Cq1dAPi15Qw9fAFdSA/CJfWNtRvmzzI2yWbvD9aUo3wOa1sPEr+HYF\nfLsSyrZEHCoP6reGvTr7AWr2PtgPZtOwnWr0MoiZvQHsHeWl38b6FlHWVZxpdggwLOL588B451yp\nmY3G35g6ppL4LgIuAmjTpk2MIUnWqNnAP6rgl5OUP4lUUFwfTrsb7j8RXvolnDYu7IgkQpUFv91l\nas65CTEco7JMLZbMLjIOZV7prrg+dDjKL+V27IANS30N3IZlsH6pH0Rlw1L/fOVc+G4tbNtcvWPl\nFfq7SDUbQt3m0Kqn71Rctzk0bA+NO0LDtr4AKRnNOXdcZa+Z2Vdm1jy4m94c+DrKZkv5obkVQCt8\nk6vy9+gKFDjnZkQcc03E9nfjm6VXFt84YBxASUlJpXmYZKnCmlBQ09+Ekpyj/EkkijaHQd9fw9tj\noWN/P0iepIUqC367y9RitBRoHfG8FbAcWA00MLOCoNavfH1lcSjzykR5eX6S8wZVFNa3bfGTIJeV\n+tFEy0ph+1awfMgrCJY8KKztC3wFxaq1E4DngOHA2OAx2s2oV4EbI0bUGwD8JuL1oUT0pwEov1gL\nnp4MzE9k0JJlajWC7zS4i+xC+ZPkrr6/gkWT4MUr/PQODduFHZGQmqae04BOwQiey/BNFs5xzjkz\nmwScgR/Zs7JMUXJBYbH62cmeGAv8x8xGAl8AZwKYWQkw2jk3yjm31syuw+dFANeWD6QQOAsYWOF9\nLzOzk4EyYC0wIolpkExXs6Gaeko0yp8kd+UX+Gaedx0JT1/kp8XK1/ThYbOd+wpXc2ezU4G/AU2B\ndcAs59zxZtYCP23DwGC7gcDt+Okc7nPO3RCs78AP0zn8FzjPOVda1XFLSkrc9OlRZ44QkQxlZjOc\ncyVhxxEP5U056oFBvh/yhS+HHYkkQTbkTaD8SUIy+wl4ehQcNcbPhSwJVd38Kd5RPZ8BnomyfjkR\nd6iccy8BL0XZbhF+1E8REZHMVLMBrF4QdhQiIunn4DP9VGCTb/ZTPLQ5LOyIclpe2AGIiIhktJqN\n1NRTRKQyA2/xo6w/9WM/AryERgU/ERGReJT38Yuj64SISNYqrgen3+NHc3/xyrCjyWkq+ImIiMSj\nZkM/CnF1p6UREckVrQ+Fo66COU/AR4+HHU3OUsFPREQkHjWDkfg1l5+ISOWOvBLa9Pa1fms/Dzua\nnKSCn4iISDxqNfKP6ucnIlK58ikeLM9P8bC9LOyIco4KfiIiIvEor/FTwU9EZPcatIFBt8LSD+Ht\nm8KOJueo4CciIhKPmkGN3+Y14cYhIpIJupwBXc+ByX+GTzT/aSqp4CciIhKP2k38owp+IiKxOekv\n0PxgP8XD1x+HHU3OMJeBw0+b2SpgE7A67FiSqAnZmz6lLXMlM31tnXNNk/TeKRHkTUti3DzM74qO\nnRvH1bETI+PzJqh2/pSusv03NppcTDPkZrr3JM3Vyp8ysuAHYGbTnXMlYceRLNmcPqUtc2V7+lIp\nzP+ljp0bx9WxlVdlm1z8XHMxzZCb6U5FmtXUU0REREREJMup4CciIiIiIpLlMrngNy7sAJIsm9On\ntGWubE9fKoX5v9Sxc+O4OrZkm1z8XHMxzZCb6U56mjO2j5+IiIiIiIjEJpNr/ERERERERCQGGVPw\nM7Mzzex/ZrbDzCod8cbMTjCzT8xsoZmNSWWM8TCzRmb2upktCB4bVrLddjObFSzPpTrO6qjqszCz\nGmb2ePD6B2bW+LXiXgAACJtJREFULvVR7pkY0jbCzFZFfFajwohzT5jZfWb2tZnNreR1M7O/Bmmf\nbWaHpDrGdBZvXmVm7YPzYUFwfhQF66s8X6qRj9xkZnOD5eyI9VMivrPLzezZYP3RZrY+4rVrknDs\nB8zs84hjdAvW7/b7loDjPhJ8DnOD735hCtOcis/65uD7OD/4P5qZ1Y1I1ywzW21mtwfbV5l3xXPs\nYP1bwf+8/Bh7xZLuONNcy8xeNLOPg9fGRmyfsfl1NkjAeXRJ8J1xZtYkYn2V53CYkphuszT9jU5i\nnpnW53AS013962jnXEYsQGdgP+AtoKSSbfKBz4AOQBHwEXBA2LHHmL6bgTHB32OAmyrZbmPYscaY\nnio/C+BnwF3B30OAx8OOO4FpGwH8PexY9zB9fYFDgLmVvD4QeBkw4DDgg7BjTqcl3rwK+A8wJPj7\nLuCnwd9Vni+x5CPAScDrQAFQG5gO1Iuy3VPA+cHfRwMvVJHuuI4NPACcUd3vWwKOOzB4bwPGR/y/\nU5HmZH/WfYB3g+9bPvA+cHSU7WYAfYO/R1BF3hXvsank3Kgq3fEcF6gF9Au2KQKmACfGmmYtyVsS\ncB51B9oBi4EmEftUeQ5nabrT9jc6AWmuLM9M63M4iemu9nV0xtT4OefmO+c+qWKzQ4GFzrlFzrmt\nwGPAKcmPLiFOAR4M/n4QGBxiLIkQy2cRmeYngWPL7winuUz+nlXJOTcZWLubTU4BHnLeVKCBmTVP\nTXTpL568Kvj+H4M/H2DnvCCW8yWWfOQA4G3nXJlzbhO+0HlC5AZmVjeI49kq0pHwY1fyvrv7vsV1\nXOfcS8F7O+BDoFVVCU3EsVP0WTugGF/QqQEUAl9FbmBmnYC98AWhWCXk2FW8b7R07/FxnXObnXOT\nAIJzbibV+6wleeI9h//rnFucikATLFnpTuff6GTlmekuzN+KnWRMwS9GLYEvI54vDdZlgmbOuRUA\nweNelWxXbGbTzWyqmaXzFz6Wz+L7bZxzZcB6oHFKootPrN+z04NmFk+aWevUhJYSmXyepYvK/oeN\ngXXB+RC5fqd9dnO+xJKPfAScGDR9awL0Ayp+P08FJjrnNkSs621mH5nZy2Z2YJT3TcSxbwjOmdvM\nrEbFdAcqft8SkmbzTTyHAa+kKM1J/6ydc+8Dk4AVwfKqc25+hc2G4u8SR470VlXelYhj3x80yfq/\niAuVqtKdkDSbWQPgR8DEaqRZkidR+VY0VZ3DYUpWutP5NzpZeSak9zkc5m/FTgriTEhCmdkbwN5R\nXvqtc25CLG8RZV3aDFu6u/RV423aOOeWm1kH4E0zm+Oc+ywxESZULJ9FWn9euxFL3M8D451zpWY2\nGn9H5pikR5Yamfq5JUwS86rd/W/L+0eVH7sVMMXMtpcfO4bj4px7zcx6Au8Bq/BN4coqbDYUuCfi\n+UygLb4GsDMw08wWRLyeiGP/BliJr6kZB1wFXItP921mVjPYrh3wmJltSdBxy/0DmOycK6/5Snaa\nk/5Zm1nHIPbymq3XzaxvUKtfbgi+wFvueXyT1xeB3sB8M1sc8Xoijn2uc25ZULP8VHD8hwia20b0\nf4lMd0LSbGYFQfr+6pxbFJnmLM2v00K81z8xnsMVzQTaOuc2mtlA/LncKfao4xdSukP9jU5imneX\nrtDP4ZDSXe3POq0Kfs654+J8i6XsfCekFbA8zvdMmN2lz8y+MrPmzrkVQZX815W8x/LgcZGZvYVv\n452OBb9YPovybZYGP8b12X0Tw3RRZdqcc2sint4N3JSCuFIlrc+zVEhiXrUa3yynILh7F/m/XQq0\nds4dF5wvK/H9Ar/P5KuRj9wA3BDs8yiwIOI9GuObop4asX15zd9xwTaL8f21Vifq2OV3Q4FSM7sf\n+GVEuv/tnBsf7PNJcOwVCUzz74GmwE9SmOZUfNanAlOdcxuDfV7G9/kpLwR1BQqcczMiYi3Pu44z\ns3xgrXPuoMg3jffYzrllwbG+Df4fh+ILfkuBPzjn3o+W7kSkGX9TYYFz7vYoaYbsy6/TQoKufyo9\nhyvZfkPE3y+Z2T/MrEnkOZxsYaSbkH+jk5jmSvPMdDiHw0g3e3AdnW1NPacBncyPflOEv5OZ1iNf\nRngOGB78PRzYpdbAzBpa0PwpqAY+HJiXsgirJ5bPIjLNZwBvVmhulK6qTJvt3J7+ZKBi86pM9hxw\nvnmHAesjLtolNlG/Q8H3fxL+fICd84JYzpdY8pH8oHCHmR0MHAy8FrHJmfgBEbZE7LO32fcjMh6K\n/+2I/KGN+9jl50xwnMFA+aiyVX3f4j3uKOB4YKhzbkeq0pyKzxr4AjjKzArMN2U9ip3zoqH42q/I\neGPJu/b42MHzJsGxCoFB7PxZ7y7dcaXZzK7HXxhdvgdpluRJRL61ixjP4TAlJd2k9290UvLMDDiH\nw/yt2JlLg9FuYlnwd/GWAqX4DuKvButbAC9FbDcQ+BRfC/bbsOOuRvoa4/sbLAgeGwXrS4B7gr/7\nAHPw7YDnACPDjruKNO3yWeCbbp0c/F0MPAEsxA+q0CHsmBOYtj8B/ws+q0nA/mHHXI20jcf3jdkW\nnHMjgdHA6OB1A+4M0j6HSkauzNUl3rwKP9Lnh8F58QRQI1hf5fkSYz5SjL9hNA+YCnSr8B5vASdU\nWHdJxPd5KtAn0ccG3gy+T3OBfwN1Yvm+JeC4ZcF7zwqWa1KY5mR/1vnAv/AXQfOAWyu8xyIq5E3E\nkHfFc2z8aHUzgNnBce4A8mNJd5zHbYVvAjU/4rMeFWuatSQ1z4z3PLoMn+eW4WtCyvep8hzO0nSn\n7W90AtJcWZ6Z1udwEtNd7etoC3YUERERERGRLJVtTT1FRERERESkAhX8REREREREspwKfiIiIiIi\nIllOBT8REREREZEsp4KfiIiIiIhIllPBT0REREREJMup4CciIiIiIpLlVPATERERERHJcv8PLEEC\novYbOVkAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "f, axarr = plt.subplots(1,3,figsize=(15,5))\n", "# Normal plot of the final map\n", "axarr[0].plot(nxf,nyf,'.')\n", "axarr[0].plot(cx,cy)\n", "axarr[0].set_title(\"The map\")\n", "\n", "# Zoomed plot of the final map (equal axis)\n", "axarr[1].plot(nxf,nyf,'.')\n", "axarr[1].plot(cx,cy)\n", "axarr[1].set_xlim([cx[-1] - 0.1, cx[-1] + 0.1])\n", "axarr[1].set_ylim([cy[-1] - 0.1, cy[-1] + 0.1])\n", "axarr[1].set_title(\"Zoom\")\n", "\n", "# Zoomed plot of the final map (unequal axis)\n", "axarr[2].plot(nxf,nyf,'.')\n", "axarr[2].plot(cx,cy)\n", "axarr[2].set_xlim([cx[-1] - 0.01, cx[-1] + 0.01])\n", "axarr[2].set_ylim([cy[-1] - 0.1, cy[-1] + 0.1])\n", "axarr[2].set_title(\"Stretch\")\n", "#axarr[1].set_xlim([cx[-1] - 0.1, cx[-1] + 0.1])\n", "#axarr[1].set_ylim([cy[-1] - 0.1, cy[-1] + 0.1])\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### How much faster is now to evaluate the Map rather than perform a new numerical integration?\n" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "# First we profile the method evaluate (note that you need to call the method 4 times to get the full state)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "11.8 µs ± 469 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "timeit xf.evaluate({\"dr\":epsilon, \"dt\":epsilon, \"dvr\":epsilon,\"dvt\":epsilon})" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "# Then we profile the Runge-Kutta 4 integrator" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "15.7 ms ± 1.97 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)\n" ] } ], "source": [ "timeit rk4(eom_kep_polar, 0, [it + epsilon for it in ic], step, n_steps)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0,0.5,'Error in estimating the final state (x)')" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEWCAYAAABMoxE0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl8VOXZ//HPlyUkBMKWsC8BRFwR\nFaHaumurrRVrbevWpy6t9emm3bXtU61d7GJbbfVXa6tFu7jUFa3Wfat1QRQQFBRZYyCQACEkJCHJ\n9fvjnOAwzExOMpnMhFzv12teOducc50zk3PNfe5z7ltmhnPOOddRvbIdgHPOue7NE4lzzrm0eCJx\nzjmXFk8kzjnn0uKJxDnnXFo8kTjnnEuLJ5IeQNKRkpZlOQaTtFcXb3OVpBOSzJsj6ScZ2u4xksoy\nse5skDRV0uuSaiR9LeJ7uvzzThBDzn0OksZL2iapd4Rlcy7+ZDyRZFh4MtsefnlaX9d3ZQxm9ryZ\nTY2ybHf68uaKXDhpZth3gGfMbKCZ/S5+pqRnJH0+C3F1O2a2xswGmFlztmPpTH2yHUAP8XEze6Kt\nhST1MbOmtqa1dx2drSu24XLKBOCObAfhcpeXSLJI0nmSXpD0W0mbgCuTTOsl6QeSVkvaIOk2SYPC\ndZSGv4gvlLQGeCrBdnYpZYSlpG9JWiSpWtKdkvIlFQKPAKNjSk+jJV0p6W5Jf5O0FThPUj9J10oq\nD1/XSuoXs41vS1oXzrsgLp5dfsGG+/yfmPH9JT0uaZOkCknfC6f3knSZpHclVUm6S9LQmPd9NjxG\nVZK+H+EjKA63UyPpWUkTwvXcIOnXcTE/KOnSBMf2uXBwYXi8PhMz75vh57VO0vkx0/tJukbSmnD/\nbpRUkCjA8Nj/LWa89fPuE3PsVoT7sFLSOTHLXiDpLUmbJT3aun9JtnOqpCWStoSfz77h9KeAY4Hr\nw/3bO+59PwWOjJkfW9o+QdI74fZvkKT2xhazvxeF36V1kr4ZdyyTfg9jlvu2pHvipv1e0rXh8DOS\nfqzgf69G0mOSits6PuG8VeH6F0mqlXSzpBGSHgnX9YSkIUk+v/PD41ATfo5fTPYZ5TQz81cGX8Aq\n4IQk884DmoCvEpQOC5JMuwBYDkwCBgD3An8N11EKGHAbUAgUJNjOMUBZXEyvAKOBocBbwMWJlg2n\nXQnsAE4j+PFRAFwFvAQMB0qA/wI/Dpc/CagADghj+kcY417h/GeAz8cdh/+EwwOBdcA3gfxwfFY4\n79Jwm2OBfsAfgdvDefsB24Cjwnm/CY9jsmM/B6iJWf66mBhmAuVAr3C8GKgDRiRZ1859izmGTeEx\n6gt8NHz/kHD+tcDc8NgPBB4Erk6y7iuBv8WMt37efcJjuxWYGs4bBewfDp9G8J3ZN1z2B8B/k2xj\nb6AWODGM9zvhe/MSfV4J3r/b/DDGh4DBwHhgI3BSB2Jr3d/bw/09MFzXCeH8VN/DYwi/y+GxqQUG\nh+N9gA3AoTH78G54LArC8Z9HPD6rwhhGAGPC9b4GHEzw3XoKuCL+8wvHPwZMBgQcHX5PDkn2v5ir\nr6wHsKe/wi/ZNmBLzOsL4bzzgDVxyyea9iTwpZjxqQQn9j4xX8xJKWLY5QsZxnRuzPgvgRsTLRtO\nuxJ4Lm7au8BHY8Y/AqwKh29p/ScMx/cmeiI5C3g9yX68BRwfMz4q5jj8ELgjZl4h0EjqRBK7/ACg\nGRgXs60Tw+GvAA+nOL6JEsn21pNFOG0D8IHwhFELTI6ZdziwMsm6ryR1ItkCfJK4HxAEJcsLY8Z7\nEZykJiTYxv8Bd8Ut+x5wTKLPK8H7d5sfxvihmPG7gMs6EFvr/u4T9329OcL38Bh2/d4/wvv/e6cA\nb8btww9ixr8E/Dvi8VkFnBMz/x7gDzHjXwXuj//8khzL+4FLkv0v5urLL211jdPMbHDM608x89Ym\nWD5+2mhgdcz4aoITyYg21pPK+pjhOoITaSpRYhodM29t3LyoxhGcHBKZANwXXl7YQnCybyY4Drts\n08xqgao2thW7/DZgU8w+3AqcGw6fC/y1HfsAUGW71iO1HuMSoD8wP2Y//h1Ob5dwHz8DXAysk/Qv\nSfuEsycA18VsYxNBEhuTYFW7fJZm1kJwbBIt2x7JvmPtia1V/Pcp9ruW7HsYr63PNFm8UY5PRczw\n9gTjCf+/JJ0s6SUFl3G3EJReixMtm8s8kWRfouaX46eVE/zztRpPcOkk9svaWc04J1tPlJjKw+F1\nBAkhdl6sWoKTaauRMcNrCYr6iawFTo5Lyvlm9l78NiX1B4YlWU+r2OUHEFxqat2HvwGzJR1EcAnm\n/jbWFVUlwYll/5h9GGRmyRJ5qmOFmT1qZicSlM6WAq0/UtYCX4w7VgVm9t8E29jlswzrMsYR/OqO\nor3fvfbE1ir++9T6OaX6Hsa7H5gm6QCCEsnfI8ab7vFJKKzLuQe4huCy6WDgYYKk2q14Iukebge+\nLmlieML7GXCnZebOqQpgmMLK/DZi+oGkkrBS8ocEJ18ILmOcJ2m/8IR+Rdx7FwCnS+qv4LbZC2Pm\nPQSMlHRpWJE6UNKscN6NwE/1fqV4iaTZ4by7gVMkfUhSHsG187a+3x+NWf7HwMtmthbAzMqAeQS/\nWu8xs+0p1lNBUH/VpvDX7J+A30oaHu7HGEkfSfKWBcBRCp4/GARc3jojrNA9VcFNEg0El1Bbbyu9\nEbhc0v7hsoMkfSrJNu4CPibpeEl9CeqnGgjqG6KIvP8diK3V/4Xfl/2B84E7w+mpvoe7MLN6gu/J\nP4BXzGxNxHjTPT7J5BHUoWwEmiSdDHw4zXVmhSeSrvGgdn2O5L52vv8WghPac8BKoJ7gumunM7Ol\nBP+cK8JLD8kuE/wEeBVYBLxBULn4k3AdjxBUKD9FUCkZfyfZbwnqLyoILjfs/GVoZjUElZofJ7jU\n8A7BXUMQVIjPBR6TVENQwTkrfN8S4MsEJ4l1wGagredh/kGQ5DYBhwLnxM2/laByt63LWlcCt4bH\n69NtLAvwXYLj8pKCu+CeIKj32o2ZPU5w0lwEzCdItK16EZzUysN9OJrg2j5mdh/wC+COcBuLgZOT\nbGMZwaWe3xOUmD5OcMt6Y4R9geBzOSO8A2u350wSbC9ybDGeJThmTwLXmNlj4fSk38Mkon6msfGm\ne3ySrbcG+BpBotoMnE3w/e52FFbqOOfiSDqK4NdtaViScF1MUinBj6e+nVEClzSe4BLgSDPbmu76\nXMBLJM4lEF7CuAT4syeRPYOkXsA3CO7W8yTSifzJdufihA+bvQosJLge77q5sB6pguDuq5OyHM4e\nxy9tOeecS4tf2nLOOZeWHnFpq7i42EpLS7MdhnPOdSvz58+vNLM2H5btEYmktLSUV199NdthOOdc\ntyIpUqsUfmnLOedcWjyROOecS4snEuecc2nxROKccy4tnkicc86lxROJc865tHgicc45l5Ye8RyJ\nc87tyVpajC3bd7CptnHna3Nd8HfskAJmT0+3s8vUPJE451yOaW6xnYmgclsDm2obqdrWSFVtI1Wt\n460JI0waLUmaTTz5gJGeSJxzrrszM2obm6na1kDltgYqtwWJoXJbQzAtTBCtyWJzXSOJ2tOVYHBB\nX4YW5jGssB97lQxg6MQ8hhXmMTR8Dekf/B02IBjO79s74/vXLROJpNOAjwHDgRtiektzzrkuYWbU\nNDSxsaaBypoGNm4L/laGCaJyWwMbtzVSWdNAVW0D9TsSd2tTlN+H4gH9GFqYx+SSAcwME8OwcFrr\n8LABeQwu6Euf3rlXtZ0ziUTSLcApwAYzOyBm+kkEXXn2Juhk6Odmdj9wv6QhwDWAJxLnXKeo39FM\n5bYGNtQ0sDF8tQ5Xbnt/2sZtDTQ27Z4cegmGFvajeEAeJQP7Mbm4kGED8ige0C9IGAPyKAmHhxT2\npV+fzJcYMi1nEgkwB7geuK11gqTewA0EfXiXAfMkzTWzN8NFfhDOd865lOoam6jY2sCGrfVU1AR/\nN9Y0ULG1ng0xyaJ6+46E7x9WGCSGkoH9mFRcuHO4NUEUDwySxZD+efTupS7eu+zKmURiZs+F/TPH\nmgksN7MVAJLuAGZLegv4OfCImb2WaH2SLgIuAhg/fnymwnbOZVljUwsbauqp2FrP+uogMVTU1FNR\nXU/F1gYqaurZuLWBmobdu3zP69OL4QP7MXxgUN9w+KRhwXhRkCSGD8ynZGBwialvDl5SyhU5k0iS\nGAOsjRkvA2YBXwVOAAZJ2svMbox/o5ndBNwEMGPGDO8G0rluqKZ+B+ur61m/tZ511fU7hyuqg/GK\nrfVU1Tbu9r683r0YXtSPEUX57DuyiKOmBMlhZFE+wwfmB/MG5lNU0AepZ5UeMiHXE0miT9jM7HfA\n77o6GOdc59nW0MS6Ldspr67f+Xd99XbWVb+fNLYlKEUMK8xjRFE+Iwflc9C4wYwsymdEUT9GDMoP\nh/MZ0r+vJ4gulOuJpAwYFzM+FijPUizOuYiamluoqGmgfMt2yrds573wb/mW+p3jNfW7JgkJSgb0\nY9TgAiaXFPKhvYoZNShIGKMGFTBqUFCS2BMqp/c0bSYSSb2Ag4DRwHZgiZlVZDqw0DxgiqSJwHvA\nmcDZXbRt51wSDU3NlG+pp2xzHe9t3k7Z5iA5vBf+Xb+1nua4J+SG9O/LqEEFjB3Sn5kThzJqUAGj\nB+czenCQJEYU5Xs9RDeVNJFImgx8l6Au4h1gI5AP7C2pDvgjcKuZJb45up0k3Q4cAxRLKgOuMLOb\nJX0FeJTg9t9bzGxJZ2zPOZfcjuYW1m2pZ+3mOtZuqmPt5jrKwoRRtrmOiq0Nuyzfu5cYWZTPmCEF\nzJo4lNGDCxgzpIAxg99PFv3zcv0CiOsoWaLHJ9l5Yv8D8LzFLSRpOEHJYLOZ3ZrxKNM0Y8YM8z7b\nndvVlrpGVlfVsWZT8Fq76f3hddW7lih69xKjB+czdnB/xgwpYNyQ/owdEiSLsUMKGFmUn5MPyrn0\nSJpvZjPaWi7pTwQzOyvFvA3AtR2MzTnXBcyMjTUNrKqqY1VlLas31bKqqo41VXWsrqpla1wdRfGA\nPMYN7c+hE4Ywbkh/xg0tYNzQ/owb0p9RgzxRuOSi1JH8GPiRmTWF40XAdWZ2fqaDc86lZmZU1Tay\nqrKWFZW1rKqsZVVVLSsrg+SxfUfzzmV79xJjhxQwfmh/Dho3mglDCxk/rD/jhwavwn5+6cl1TJRv\nTh/gZUnnAyOB34cv51wXqd/RzKqqWlZsrOXdDdtYESaOlRu37VKy6NNLjBvan9Jh/fnApKFMLC5k\nwrBCJgwNLkl5ZbbLhDYTiZldLulJ4GVgM3CUmS3PeGTO9UBb6hpZvmHbzte7G7fx7sZa1m6u26U1\n2FGD8plUUsjs6WMoLS5kUkkhE4cVMnZIgV+Ccl0uyqWtowgaTbwKOBC4XtIFZubPczjXQZtrG1lW\nUcM7FTW8s2Ebb1fUsHxDLZXb3r8bql+fXkwsLmTa2EF84uAxTB4+gElh0vA7oFwuifJtvAb4VGtD\niZJOB54C9slkYM7tCWobmni7ooZl62tYVlETDm/bJWEM7NeHycMHcNw+Jew1fAB7DR/AlOEDGT24\noMc1/ue6pyiJ5HAz21ljZ2b3Sno2gzE51+20tBhrN9fx1roa3lq3laXrt7J0fQ2rq+p2LlPQtzd7\njxjAsVNLmDpyIFNGDGTvEQMYWZTvzXm4bi3VA4nnAv+ITSKtzKwqfGBxlJn9J5MBOpdrGpqaeadi\nG0vKq3mzfCtvrtvKW+tqdrYLJcHEYYUcMHoQZxwylqkjB7LPyCLGDimgl5cw3B4oVYlkGPC6pPnA\nfN5/sn0v4GigErgs4xE6l0X1O5pZur6GN96rZnFZNYvLq3m7ooYdzUHNd2Feb/YdVcTph4xh31FF\n7DuqiKkjBlKQ5+1BuZ4j1QOJ10m6HjgO+CAwjaCtrbeAz5rZmq4J0bmu0dTcwtsV21hYtoVFZVtY\nVFbNsvU1NIVPeA/u35cDxwziwg9NYv/RRew/uojSYYVeynA9Xso6kvCy1uPhy7k9SvmW7by+ZgsL\n1m5mwdotvPFe9c5+tQcV9GXa2EFcdNQkDhwziAPHDmLM4AKvy3AuAb+H0PUIjU0tLCmvZv7qzby2\nZjOvrd7C+q31QNBL3gGjizh75gQOGjeI6eMGM35of08azkXkicTtkbbW7+C11Zt5ddVmXlm1iYVr\nt9DQFJQ2xg4pYNakoRw8bjAHjx/CvqOKyOvjD/E511GeSNweobpuBy+vrOKlFZt4ZVUVb5ZvpcWC\n9qUOGF3EObMmcFjpEA6dMIThRfnZDte5PUqUJ9tHAD8DRpvZyZL2I3i25OaMR+dcEjX1O5i3ahP/\nXV7Ff9+t4q31WzELngY/ZPwQvnrcFGZOHMrB4wf7U+DOZViU/7A5wF+A74fjbwN3Ap5IXJfZ0dzC\n62u28J/llfznnY0sLKumucXI69OLQ8cP4esn7M0HJg3joHGDvCtW57pYlERSbGZ3SbocwMyaJO32\nkKJznW3tpjqefXsjz769kRffrWJbQxO9BAeOHczFR0/ig5OLOWTCEPL7euJwLpuiJJJaScMAA5D0\nAaA6o1G5HmlHcwvzVm7i6WUbeHrZRpZv2AbAmMEFfPyg0Rw1pZgjJhczqH/fLEfqnIsVJZF8A5gL\nTJb0AlACfCqjUbkeo7puB08v28Djb1Xw3LKN1DQ0kde7F7MmDeWsmeM5eu8SJpcU+q24zuWwKIlk\nCUGTKFMBAcsAv1fSdVjF1noeXbKex5ZU8NKKKppajOIBeZx84EiO33cEH9qr2Hvrc64bifLf+qKZ\nHUKQUACQ9BpwSMaicnuc8i3befiNdTyyeD3zV28GYFJxIZ8/chIn7jeCg8cN9qZGnOumUrX+OxIY\nAxRIOpigNAJQBPTvgthcN7expoF/LSrnoUXreDVMHvuOKuIbJ+7NyQeMZMqIgVmO0DnXGVKVSD4C\nnAeMBX4TM70G+F4GY3LdWG1DE48uWc/9C8p5YXklzS3GPiMH8s0T9+Zj00YxqWRAtkN0znWyVK3/\n3grcKumTZnZPF8bkupmWFuOllVXcM/89Hlm8jrrGZsYOKeDioycxe/oY9vaSh3N7tDbrSMzsHkkf\nA/Yn6I+kdfpVmQzM5b711fXcPX8td766lrWbtjOwXx9OPWg0px8ylsNKh/idVs71EFGaSLmRoE7k\nWODPwBnAKxmOq62YJhE8aT/IzM7IZiw9TUuL8ew7G/n7S2t4amkFLQZHTB7Gtz48lY/sP9IfDnSu\nB4py19YRZjZN0iIz+5GkXwP3dnSDkm4BTgE2mNkBMdNPAq4DegN/NrOfJ1uHma0ALpR0d0fjcO1T\nXbeDO19dw19fWs3aTdspHpDHF4+ezJmHjWPCsMJsh+ecy6IoiWR7+LdO0migCpiYxjbnANcDt7VO\nkNQbuAE4ESgD5kmaS5BUro57/wVmtiGN7bt2WL6hhlteWMV9r73H9h3NzJw4lO98ZB8+sv9Ib3rd\nOQdESyQPSRoM/Ap4jaCplD93dINm9pyk0rjJM4HlYUkDSXcAs83saoLSS7tJugi4CGD8+PEdDbdH\nMjNeXrmJPz23gieXbqBfn16cNn0MnzuilP1GF2U7POdcjomSSH5pZg3APZIeIqhwr+/kOMYAa2PG\ny4BZyRYO2/76KXCwpMvDhLMLM7sJuAlgxowZ1rnh7pnMjCfe2sANTy9nwdotDC3M49ITpvDZD0xg\n2IB+2Q7POZejIj3ZTvgUe5hQGjLwZHui23uSnvzNrAq4uBO336O1tBgPL17H9U8tZ+n6GsYNLeDH\npx3Apw4d65Xnzrk25cqT7WXAuJjxsUB5J2/DxWlpMR57cz2/ffwdllXUsNfwAfzm0wdx6kGj6dPb\n6z+cc9FEfbL917yfSDLxZPs8YIqkicB7wJnA2Z28DRfj+Xc28ot/L2Xxe1uZVFzIdWdO55Rpo+nt\n7V0559qpy59sl3Q7cAxQLKkMuMLMbpb0FeBRgju1bjGzJSlW4zpoSXk1P39kKc+/U8mYwQX8+lMH\nMXu6l0Cccx0XpY5krKQigpLInwjqRi4zs8c6skEzOyvJ9IeBhzuyTte2ym0NXPPoMu58dS2DCvry\ng4/ty2cPn+Dd0jrn0hYlkVxgZtdJ+ggwHDifoA/3DiUS17V2NLdw639Xcd0T77B9RzMXfnAiXz1+\nCoMKvJdB51zniJJIWi+afxT4i5ktlDei1C3MX72J79+3mKXrazh2agk/OGU/Jnvru865ThYlkcyX\n9BjB0+yXSxoItGQ2LJeOrfU7uPrhpdz+yhpGD8rnps8eyof3H5ntsJxze6goieRCYDqwwszqwocB\nz89sWK6jnl62ge/d+wYVW+v5wpETufSEvb3bWudcRkVpRr6FoGmU1vEqgva2XA7Z1tDEVQ8u4a5X\ny5gyfAA3fumDHDRucLbDcs71AP5TdQ/w+prNXHLHAso21/G/x0zmkuOn+BPpzrku44mkG2tpMf7w\n7Lv85vG3GVmUz51fPJzDSodmOyznXA+TqomUlGckM9vU+eG4qDbXNnLpnQt49u2NnDJtFD/9xIF+\nS69zLitSlUjmEzScmKxBxUkZici16Y2yai7+23w21jTwk9MO4JxZ471bW+dc1qRqIiWdzqtchtz/\n+nt8555FFBfmcdfFhzPdK9Sdc1kWqY5E0hBgCkFfJEDQQVWmgnK7MzN++8Q7/O7Jd5g1cSh/OPdQ\nhhbmZTss55xrO5FI+jxwCUErwAuADxD0UXJcZkNzrep3NPOduxcxd2E5Zxw6lp994kDv5tY5lzOi\nnI0uAQ4DVpvZscDBwMaMRuV22tbQxAVz5jF3YTnfOWkqvzpjmicR51xOiXJpq97M6iUhqZ+ZLZU0\nNeOROTbXNnLenHksfq+a337mID5x8Nhsh+Scc7uJkkjKJA0G7gcel7QZ770w4zZsrefcm19mVVUd\nfzjnEG8ryzmXs6I0kfKJcPBKSU8Dg4B/ZzSqHq5yWwNn//llyrdsZ875h3HE5OJsh+Scc0lFvWur\nNzACWBlOGgmsyVRQPdmWukbO/fPLlG2u49bzZzJr0rBsh+SccylFuWvrq8AVQAXvNx9vwLQMxtUj\nba3fwf/c8gorKmu5+XMzPIk457qFKCWSS4CpYau/LkMam1q4+K/zebN8K3/87KEcOaUk2yE551wk\nURLJWqA604H0ZGbG9+97g/++W8WvP3UQx+87ItshOedcZFESyQrgGUn/AhpaJ5rZbzIWVQ/z/555\nl3/OL+Nrx0/hk4f6Lb7Oue4lSiJZE77ywpfrRA8tKudXjy5j9vTRfP2EKdkOxznn2i3K7b8/6opA\neqLlG7bxnbsXceiEIfzyjGnegq9zrltK1R/JtWZ2qaQHCe7S2oWZnZrRyPZw2xub+fLfXyO/b2+u\nP/tg+vXxHg2dc91TqhLJbeHfa7oikPaSVAg8B1xhZg9lO572+uEDi3l7Qw23nj+TUYMKsh2Oc851\nWKrW/34V/v2omT0b/+roBiXdImmDpMVx00+StEzSckmXRVjVd4G7OhpHNt09v4x/zi/jq8fuxVF7\n+22+zrnuLVWJZJSko4FTJd1BXE+JZvZaB7c5B7ie90s8rU/O3wCcCJQB8yTNBXoDV8e9/wKChyHf\nJKZ/lO5iXfV2rpy7hFkTh3LJCXtnOxznnEtbqkTyQ+Aygn5Ifs2uicToYH8kZvacpNK4yTOB5Wa2\nAiBMXLPN7GrglPh1SDoWKAT2A7ZLetjMWuKXyzXB8yKLaW4xfnnGNHr38sp151z3l6qr3buBuyX9\nn5n9OMNxjCF48LFVGTArRWzfB5B0HlCZKIlIugi4CGD8+PGdGWuHPbCgnKeWbuD/TtmPCcMKsx2O\nc851ijZ7SOqCJAJxl81aN93Wm8xsTrKKdjO7ycxmmNmMkpLs10NsrGngygeXcMj4wZx3RGm2w3HO\nuU6TK13tlQHjYsbHsof1efLjh96krqHZL2k55/Y4uZJI5gFTJE2UlAecCczNckyd5rU1m5m7sJyL\nj5nMXsMHZjsc55zrVJESiaQPSTo/HC6RNLGjG5R0O/AiMFVSmaQLzawJ+ArwKPAWcJeZLenoNnKJ\nmfGzf71FycB+fPGoSdkOxznnOl2U/kiuAGYAU4G/AH2BvwEf7MgGzeysJNMfBh7uyDpz2aNL1vPq\n6s1cffqBFPaL1I+Yc851K1FKJJ8ATgVqAcysHPDrMxE0NrXw80eWMmX4AD7lrfo65/ZQURJJo5kZ\n4V1UYdMkLoJ/vLyaVVV1fO+j+9Knd65URznnXOeKcna7S9IfgcGSvgA8Afwps2F1f/U7mrn+6Xc5\nfNIwjpma/duPnXMuU6I0I3+NpBOBrQT1JD80s8czHlk398CC96jc1sB1Z0735uGdc3u0SLW/YeLw\n5BFRS4vxp+dXst+oIo6YPCzb4TjnXEa1eWlL0umS3pFULWmrpBpJW7siuO7q2bc3snzDNr5w1EQv\njTjn9nhRSiS/BD5uZm9lOpg9xU3PrWBkUT6nTBud7VCccy7jolS2V3gSiW7xe9W8uKKK8z9YSl+/\nU8s51wOk6mr39HDwVUl3AvcDDa3zzezeDMfWLf3p+RUM6NeHs2blRovDzjmXaakubX08ZrgO+HDM\nuAGeSOJU1+3g4TfWcc6sCRTl9812OM451yVS9UfS2rbWB83shdh5kjrUPMqe7t9L1rGj2Tj9kDHZ\nDsU557pMlIv4v484rcd7YEE5E4sLOXDMoGyH4pxzXSZVHcnhwBFAiaRvxMwqIuhL3cWo2FrPiyuq\n+NpxU/yWX+dcj5KqjiQPGBAuE9tI41bgjEwG1R09uLAcMzh1ut/y65zrWVLVkTwLPCtpjpmt7sKY\nuqUHF5ZzwJgiJpcMyHYozjnXpaL02e5JpA0rK2tZWFbN7IO8kt051/P4E3OdYO6CciQ45aBR2Q7F\nOee6nCeSTjB34XvMmjiUUYMKsh2Kc851uShd7f4uweRq4FUze6DzQ+peyjbX8e7GWs6ZNSHboTjn\nXFZEKZHkA9OBd8LXNGAocKGkazMYW7fw8opNABzuzcU753qoKK3/7gUcZ2ZNAJL+ADwGnAi8kcHY\nuoWXVlQxuH9fpo7wbuydcz18sJlwAAAW5klEQVRTlBLJGCC2n/ZCYLSZNRPTiGNP9fLKTcwsHUqv\nXv4QonOuZ4raH8kCSc8AAo4CfiapkKD/9h7rvS3bWbOpjvOOKM12KM45lzVR+my/WdLDwEyCRPI9\nMysPZ387k8HlupdXVAHwgUleP+Kc67mi3v7bC9gIbAL2knRU5kLqPl5aUcWggr7sM9LrR5xzPVeU\n239/AXwGWAK0hJMNeC6DcbUV05HAOQTx72dmR2QjjpdXbmLmRK8fcc71bFHqSE4DpppZp1SsS7oF\nOAXYYGYHxEw/CbiOoGXhP5vZz5Otw8yeB56XdBowrzPiaq/yLdtZXVXH/xxemo3NO+dczoiSSFYA\nfem8O7TmANcDt7VOkNQbuIHgluIyYJ6kuQRJ5eq4919gZhvC4bOBz3dSXO3y8srW+pGh2di8c87l\njCiJpI7grq0n2bXP9q91ZINm9pyk0rjJM4HlZrYCQNIdwGwzu5qg9LIbSeOBajPbmmT+RcBFAOPH\nd37/6S+9u4lBBX3Zd2RRp6/bOee6kyiJZG74yqQxwNqY8TJgVhvvuRD4S7KZZnYTcBPAjBkzLN0A\n4720sorD/PkR55yLdPvvrV0QR6KzccqTv5ldkaFY2rSlrpHVVXWcPbPzSzrOOdfdpOpq9y4z+7Sk\nN0hwUjezaZ0YRxkwLmZ8LFCeZNmsW1FZC+CdWDnnHKlLJJeEfxPWUXSyecAUSROB94AzCSrSc9LK\njUEimVhS2MaSzjm350v6QKKZrQsHv2Rmq2NfwJc6ukFJtwMvAlMllUm6MGwQ8ivAo8BbwF1mtqSj\n28i0VVW19BKMG9I/26E451zWRalsPxH4bty0kxNMi8TMzkoy/WHg4Y6ss6utqKxl3ND+5PXxfsGc\ncy5VHcn/EpQ8JklaFDNrIPBCpgPLZasqaykd5pe1nHMOUpdI/gE8QvBA4GUx02vMbFNGo8phZsaq\nyloOK/UHEZ1zDlIkEjOrJuhS9ywAScMJekscIGmAma3pmhBzy8aaBmobm5lY7CUS55yDCK3/Svq4\npHeAlcCzwCqCkkqPtDK89bfUE4lzzgHRmpH/CfAB4G0zmwgcTw+uI2lNJJM8kTjnHBAtkewwsyqg\nl6ReZvY0MD3DceWslVW15PXuxejBBdkOxTnnckKU23+3SBpA0P/I3yVtAJoyG1buWlVZy/hh/ent\nbWw55xwQrUQyG9gOfB34N/Au8PFMBpXLVvqtv845t4sojTbWAkgqAh7MeEQ5rKXFWF1VxzFTh2c7\nFOecyxlRutr9InAVQamkhaClXgMmZTa03LNuaz0NTS1eInHOuRhR6ki+BexvZpWZDibXtTbWWFrs\nbWw551yrKHUk7xL0ktjjraxqvfXXm493zrlWUUoklwP/lfQyndDVbne2cmMtBX17M6KoX7ZDcc65\nnBElkfwReAp4g6COpMdaVVVLaXEhkt/665xzraIkkiYz+0bGI+kGVlXWss+ogdkOwznnckqUOpKn\nJV0kaZSkoa2vjEeWY5qaW1izqc4ba3TOuThRSiStXd5eHjOtx93+W1HTQFOLMdZ7RXTOuV1EeSBx\nYlcEkutqG4JWYQbmR8m9zjnXc6TqIfE4M3tK0umJ5pvZvZkLK/e0JpLCPE8kzjkXK9VZ8WiCu7US\ntatlQI9KJHWNzQD0z+ud5Uiccy63pOoh8Ypw8CozWxk7T1KPu9y1s0TSz0skzjkXK8pdW/ckmHZ3\nZweS67xE4pxziaWqI9kH2B8YFFdPUkTQd3uPUtvoJRLnnEsk1VlxKnAKMJhd60lqgC9kMqhcVNfg\nJRLnnEskVR3JA8ADkg43sxe7MKac1Foi6e93bTnn3C6i1JF8QlKRpL6SnpRUKencjEcWkjRJ0s2S\n7o6ZVijpVkl/knROV8RR19hMft9e3sWuc87FiZJIPmxmWwkuc5UBewPfjrJySbdI2iBpcdz0kyQt\nk7Rc0mWp1mFmK8zswrjJpwN3m9kXgFOjxJKu2oYmf4bEOecSiHJm7Bv+/Shwu5ltakfrt3OA64Hb\nWidI6g3cAJxIkJjmSZoL9Aaujnv/BWa2IcF6xxK0RgzQHDWYdNQ1NtO/n9ePOOdcvCiJ5EFJSwm6\n2v2SpBKgPsrKzew5SaVxk2cCy81sBYCkO4DZZnY1QaknijKCZLKAJKUqSRcBFwGMHz8+4mqT8xKJ\nc84l1ualLTO7DDgcmGFmOwh6S5ydxjbHAGtjxsvCaQlJGibpRuBgSa0NR94LfFLSH4AHk8R9k5nN\nMLMZJSUlaYQbqGts9ju2nHMugTZ/YkvqD3wZGE/wC380wa3BD3Vwm4mui1myhc2sCrg4blotcH4H\nt98htY1NDPBnSJxzbjdRKtv/AjQCR4TjZcBP0thmGTAuZnwsUJ7G+rpEXYOXSJxzLpEoiWSymf0S\n2AFgZttJXKqIah4wRdJESXnAmcDcNNbXJWobvY7EOecSiZJIGiUVEF5+kjQZaIiyckm3Ay8CUyWV\nSbrQzJqArwCPAm8Bd5nZkg5F34X8ri3nnEssyk/sK4B/A+Mk/R34IHBelJWb2VlJpj8MPBwxxpzg\nd20551xiUXpIfFzSa8AHCC5pXWJmlRmPLIc0NbfQ0NTizaM451wCkc6M4Z1T/8pwLDmrbkfwzGOh\nX9pyzrndRKkj6fHeb/nXSyTOORfPE0kE7/dF4iUS55yLF+kndtg+1ojY5c1sTaaCyjVeInHOueSi\nPNn+VYI7tyqAlnCyAdMyGFdO2Vki8QcSnXNuN1F+Yl8CTA0r3HukutZOrbyJFOec202UOpK1QHWm\nA8llteGlLS+ROOfc7qL8xF4BPCPpX8Q80W5mv8lYVDnGSyTOOZdclDPjmvCVF756HC+ROOdcclGe\nbP9RVwSSy3aWSPyuLeec203SM6Oka83sUkkPkqC/EDPrkr7Sc0FtYzN9e4u8Pv7YjXPOxUv1E/uv\n4d9ruiKQXFbX0OSlEeecSyLp2dHM5od/n+26cHJTbWOz148451wSfq0mgrrGJr9jyznnkvBEEkFt\ng5dInHMumZSJRFJvSb/qqmByVV2j15E451wyKROJmTUDh0pKp4/2bq+2odlb/nXOuSSi/Mx+HXhA\n0j+B2taJZnZvxqLKMV4icc655KKcHYcCVcBxMdMM6DGJpLbRSyTOOZdMlCfbz++KQHKZP0finHPJ\ntXnXlqSxku6TtEFShaR7JI3tiuByQUuLUbfD79pyzrlkotz++xdgLjAaGAM8GE7rEeqbmjHzln+d\ncy6ZKImkxMz+YmZN4WsOUJLhuHKGt/zrnHOpRUkklZLODZ8p6S3pXILK9x7BW/51zrnUoiSSC4BP\nA+uBdcAZ4bQuIWmSpJsl3R0zbV9JN0q6W9L/ZnL7O0skfteWc84l1OaT7cAnzexUMysxs+FmdpqZ\nrY6yckm3hJX0i+OmnyRpmaTlki5LtQ4zW2FmF8ZNe8vMLiZIcDOixNJRXiJxzrnUojzZPjuN9c8B\nToqdECanG4CTgf2AsyTtJ+lASQ/FvYYnW7GkU4H/AE+mEV+bahu9ROKcc6lE+Zn9gqTrgTvZ9cn2\n19p6o5k9J6k0bvJMYLmZrQCQdAcw28yuBk6JGDdmNheYG/Yl/4/4+ZIuAi4CGD9+fNTV7qauwUsk\nzjmXSpSz4xHh36tiphm7PuneHmOAtTHjZcCsZAtLGgb8FDhY0uVmdrWkY4DTgX7Aw4neZ2Y3ATcB\nzJgxY7ceHqPaWSLxROKccwmlPDtK6gX8wczu6sRtJmoAMumJ3syqgIvjpj0DPNOJMSW1s47EL205\n51xCbdWRtABf6eRtlgHjYsbHAuWdvI1O8/5zJF4icc65RKLc/vu4pG9JGidpaOsrjW3OA6ZImigp\nDziT4Mn5nFTX2IQE+X29DzDnnEskys/s1mdGvhwzzYBJbb1R0u3AMUCxpDLgCjO7WdJXgEeB3sAt\nZrakXVF3oaB3xD708C5ZnHMuqSit/07s6MrN7Kwk0x8mSSV5rgn6IvH6EeecSybp9RpJ34kZ/lTc\nvJ9lMqhcEvRF4vUjzjmXTKoL/2fGDF8eN+8keoigLxIvkTjnXDKpEomSDCca32PVNjb5HVvOOZdC\nqkRiSYYTje+x6hqb/RkS55xLIdVP7YMkbSUofRSEw4Tj+RmPLEfUNjQxbkj/bIfhnHM5K2kiMTP/\nGU5YIvE6EuecS8qfsmtDbUOT37XlnHMpeCJJwcy8ROKcc23wRJJCY3MLTS3mJRLnnEvBE0kKdWGD\njV4icc655DyRpFAbNiHvz5E451xynkhSqAs7tfLnSJxzLjlPJCnk9+nNx6aN8udInHMuBb9mk8L4\nYf254exDsh2Gc87lNC+ROOecS4snEuecc2nxROKccy4tnkicc86lxROJc865tHgicc45lxZPJM45\n59LiicQ551xaZLbn95oraSOwOo1VFAOVnRROd9ET9xl65n73xH2Gnrnf7d3nCWZW0tZCPSKRpEvS\nq2Y2I9txdKWeuM/QM/e7J+4z9Mz9ztQ++6Ut55xzafFE4pxzLi2eSKK5KdsBZEFP3GfomfvdE/cZ\neuZ+Z2SfvY7EOedcWrxE4pxzLi2eSJxzzqWlRycSSSdJWiZpuaTLEszvJ+nOcP7Lkkpj5l0eTl8m\n6SNdGXe6Orrfkk6UNF/SG+Hf47o69o5K57MO54+XtE3St7oq5s6Q5nd8mqQXJS0JP/P8roy9o9L4\nfveVdGu4r29JuryrY09HhP0+StJrkpoknRE373OS3glfn2v3xs2sR76A3sC7wCQgD1gI7Be3zJeA\nG8PhM4E7w+H9wuX7ARPD9fTO9j51wX4fDIwOhw8A3sv2/mR6n2Pm3wP8E/hWtveniz7rPsAi4KBw\nfFh3+I6nuc9nA3eEw/2BVUBptvepE/e7FJgG3AacETN9KLAi/DskHB7Snu335BLJTGC5ma0ws0bg\nDmB23DKzgVvD4buB4yUpnH6HmTWY2Upgebi+7qDD+21mr5tZeTh9CZAvqV+XRJ2edD5rJJ1G8M+1\npIvi7Szp7PeHgUVmthDAzKrMrLmL4k5HOvtsQKGkPkAB0Ahs7Zqw09bmfpvZKjNbBLTEvfcjwONm\ntsnMNgOPAye1Z+M9OZGMAdbGjJeF0xIuY2ZNQDXBL7Mo781V6ex3rE8Cr5tZQ4bi7Ewd3mdJhcB3\ngR91QZydLZ3Pem/AJD0aXg75ThfE2xnS2ee7gVpgHbAGuMbMNmU64E6Szjkp7fNZn/YsvIdRgmnx\n90InWybKe3NVOvsdzJT2B35B8Ku1O0hnn38E/NbMtoUFlO4knf3uA3wIOAyoA56UNN/MnuzcEDtd\nOvs8E2gGRhNc4nle0hNmtqJzQ8yIdM5JaZ/PenKJpAwYFzM+FihPtkxY3B0EbIr43lyVzn4jaSxw\nH/A/ZvZuxqPtHOns8yzgl5JWAZcC35P0lUwH3EnS/Y4/a2aVZlYHPAwckvGI05fOPp8N/NvMdpjZ\nBuAFoLu0xZXOOSnt81lPTiTzgCmSJkrKI6h0mxu3zFyg9Q6GM4CnLKidmgucGd79MRGYArzSRXGn\nq8P7LWkw8C/gcjN7ocsiTl+H99nMjjSzUjMrBa4FfmZm13dV4GlK5zv+KDBNUv/wZHs08GYXxZ2O\ndPZ5DXCcAoXAB4ClXRR3uqLsdzKPAh+WNETSEIIrDY+2a+vZvtsgmy/go8DbBHc7fD+cdhVwajic\nT3CnznKCRDEp5r3fD9+3DDg52/vSFfsN/IDgGvKCmNfwbO9Ppj/rmHVcSTe6ayvd/QbOJbjBYDHw\ny2zvS6b3GRgQTl9CkDS/ne196eT9Poyg9FELVAFLYt57QXg8lgPnt3fb3kSKc865tPTkS1vOOec6\ngScS55xzafFE4pxzLi2eSJxzzqXFE4lzzrm0eCJxKUlqlrRA0mJJ/5TUv53vv7S97wnfd2V7WtqV\nNFjSl2LGR0u6u73bbQ9JR4Yt4y6QVJAilmMkPZTGdi6W9D9tLDND0u9itndEO9/fruOdDkmlkhaH\nw2nF7XKDJxLXlu1mNt3MDiBoxO7iqG+U1JvgafD2Jp+ONN0zmKBVVwDMrNzMzkixfGc4h6A9pulm\ntj1ZLOkysxvN7LY2lnnVzL4Wjh4DHBEzr833Z0t3jdvtyhOJa4/ngb0AJJ0r6ZXw1/gfw6SBgj47\nrpL0MsFDm6OBpyU93Tq/dWWSzpA0JxyeI+k34XK/CBc5SNJTYR8JXwiXGyDpybAhwTcktbZw+nNg\nchjPr+J+9eZL+ku4/OuSjg2nnyfpXkn/Drfxy0Q7Len48H1vSLolbNHg88CngR9K+nvcW3aJJZw2\nQNLdkpZK+ru0s2XhQyU9q6B/l0cljUqw/Z2lBUnPSPpFeOzflnRkOP0YSQ8p6FvjYuDr4faPjHv/\nFyTNk7RQ0j1tlRYljZB0X7j8wtYSg6RvhKXUxZIuDaeVKujH409hSe2x1pJauJ8LJb0IfDlm/VHj\nni7pJUmLwniGtHE89o/5fi6SNCXVfro0ZftpTH/l9gvYFv7tAzwA/C+wL/Ag0Dec9/8I2t6CoLG3\nT8e8fxVQHL++cPgMYE44PAd4iLDPC4KnyBcSNOddTNA66egwjqJwmWKCJ3FF0NfC4ph17xwHvgn8\nJRzeh6ApjHzgPILm4QeF46uBcXH7nx9ue+9w/Dbg0piYz0hwzOJjOYaghdmxBD/eXiRoELEv8F+g\nJFzuM8AtCdZ3JeET9cAzwK/D4Y8CT8Rs46H45RO8f1jM9J8AX030nphl7ozZ397hsToUeAMoJHga\nfAlBXzWlQBMwPVz+LuDccHgRcHQ4/KuYzyZq3LHvvwq4to3j8XvgnHA4DyjI9v/SnvzyEolrS4Gk\nBcCrBCfgm4HjCU4m88J5xxN0qANB66n3dHBb/7Rd+7x4wMy2m1kl8DRB66wCfiZpEfAEQXPXI9pY\n74eAvwKY2VKChLF3OO9JM6s2s3qCZjEmxL13KrDSzN4Ox28FjurAvr1iZmVm1kLQtExpuO4DgMfD\n4/gDgmTTlnvDv/PD9bTHAZKel/QGwaW5/dtY/jjgDwBm1mxm1QTH8z4zqzWzbWE8R4bLrzSzBbHx\nSRoEDDazZ8Ppf21PwAneH/8ZJDoeLxI0sPldYILteunRdbKe3Iy8i2a7mU2PnRBelrnVzBJ1RVpv\nqTtAim2TJ77r1toUy7aOnwOUAIea2Q4FrfK21QVsqvbfY/tTaWb3/4nOajs+0XZE0N7R4R1cV6J4\n2zIHOM3MFko6j6BE0F7tOZ4F4fKZbItpt+NhZv8IL69+DHhU0ufN7KkMxtCjeYnEdcSTwBmShgNI\nGiop/pd8qxpgYMx4haR9JfUCPtHGdmaH9RvDCE548wgurWwIk8ixvF+CiN9OrOcIEhCS9gbGEzS2\nGcVSgl/Ve4XjnwWeTbF8W7HEWgaUSDo8jK2vgr5e0pVq+wOBdZL6Eh6TNjxJcDkTSb0lFREcz9MU\ntAxcSPA5Pp9sBWa2BaiW9KFwUrLtJow7LAVtbq3/IMJnIGkSsMLMfkfQCu60VMu79Hgice1mZm8S\nXIZ5LLzE9DiwWyVx6CbgEYWV7cBlBHUhTxH0RJfKKwTN1r8E/NiCbn7/DsyQ9CrBCWlpGFMV8EJY\n+furuPX8P6B3eDnnTuA8i9izY3jJ63zgn+H7W4Ab23hPqlhil2skqCf6haSFBJe8jki2fDs8CHyi\ntdI6bt7/AS8TfGZRmki/BDg23Pf5wP5m9hpByeaVcF1/NrPX21jP+cANYWV7sstMqeL+HPCr8Ps2\nnaCeJJXPAIvDS4b7ENRtuQzx1n+dc86lxUskzjnn0uKJxDnnXFo8kTjnnEuLJxLnnHNp8UTinHMu\nLZ5InHPOpcUTiXPOubT8f75XNeSkk3TiAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# It seems the speedup is 2-3 orders of magnitude, but did we loose precision?\n", "# We plot the error in the final result as computed by the HOTM and by the Runge-Kutta\n", "# as a function of the distance from the original initial conditions\n", "out = []\n", "pert = np.arange(0,0.1,1e-3)\n", "for epsilon in pert:\n", " res_map_xf = xf.evaluate({\"dr\":epsilon, \"dt\":epsilon, \"dvr\":epsilon,\"dvt\":epsilon})\n", " res_int = rk4(eom_kep_polar, 0, [it + epsilon for it in ic], step, n_steps)\n", " res_int_x = [it[0]*np.sin(it[2]) for it in res_int]\n", " res_int_xf = res_int_x[-1]\n", " out.append(np.abs(res_map_xf - res_int_xf))\n", "plt.semilogy(pert,out)\n", "plt.title(\"Error introduced by the use of the polynomial\")\n", "plt.xlabel(\"Perturbation of the initial conditions\")\n", "plt.ylabel(\"Error in estimating the final state (x)\")" ] } ], "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.6.4" } }, "nbformat": 4, "nbformat_minor": 2 }