{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Figure 1\n", "\n", "This notebook recreates Figure 1 in Rein & Tamayo 2016. The figure illustrates the use of second order variational equations in an $N$-body simulation.\n", "\n", "We start by import the REBOUND, numpy and matplotlib packages. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import rebound\n", "import numpy as np\n", "%matplotlib inline\n", "import matplotlib\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We setup a planetary system with two Jupiter mass planets. The following function takes that system, integrates it forward in time by 10 orbits and returns the inner planet's $x$ coordinate at the end of the simulation. The $x$ coordinate changes because the planet orbits the star, but also because the planet interacts with the other planet. The function takes the outer planet's initial semi-major axis, $a$, as a parameter. We setup the system using heliocentric coordinates and therefore specify the `primary` attribute when adding particles to REBOUND (by default REBOUND uses Jacobi coordinates which are not supported by variational equations)." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def run_sim(a):\n", " sim = rebound.Simulation()\n", " sim.add(m=1.)\n", " sim.add(primary=sim.particles[0],m=1e-3, a=1)\n", " sim.add(primary=sim.particles[0],m=1e-3, a=a)\n", " \n", " sim.integrate(2.*np.pi*10.)\n", " return sim.particles[1].x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now run this simulation 400 times for different initial $a$ in the range [1.4, 1.7] and store the final $x$ coordinate of the inner planet in the array `x_exact`." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [], "source": [ "N=400\n", "x_exact = np.zeros((N))\n", "a_grid = np.linspace(1.4,1.7,N)\n", "for i,a in enumerate(a_grid):\n", " x_exact[i] = run_sim(a)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we create a function that runs an $N$-body simulation including first and second order differential equations. For that we add two sets of variational particles with the `add_variation()` command (one for first order and one for second order). We then initialize the variational particles by varying the outer planet's semi-major axis. After integrating the system forward in time, the function returns the $x$ coordinate of the inner planet as well as the $x$ coordinate of the corresponding variational particles: $\\partial x/\\partial a$ and $\\partial^2 x/\\partial a^2$. Note that the variational particles' *position coordinates* are thus unit-less. " ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def run_sim_var(a):\n", " sim = rebound.Simulation()\n", " sim.add(m=1.)\n", " sim.add(primary=sim.particles[0],m=1e-3, a=1)\n", " sim.add(primary=sim.particles[0],m=1e-3, a=a)\n", " var_da = sim.add_variation()\n", " var_dda = sim.add_variation(order=2, first_order=var_da)\n", " var_da.vary(2, \"a\")\n", " var_dda.vary(2, \"a\")\n", " \n", " sim.integrate(2.*np.pi*10.)\n", " return sim.particles[1].x, var_da.particles[1].x, var_dda.particles[1].x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We run **one** simulation with variational particles at $a_0=1.56$. We then use the derivates we got from the `run_sim_var()` function to approximate the final position of the inner planet as a function of the outer planet's initial semi-major axis using a Taylor series:\n", "$$x(a) = x(a_0)+\\frac{\\partial x}{\\partial a}\\Big|_{a_0} (a-a_0)+\\frac 12\\frac{\\partial^2 x}{\\partial a^2}\\Big|_{a_0} (a-a_0)^2 + \\ldots$$\n", "We do this to both first and second order." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [], "source": [ "a_0 = 1.56\n", "x, dxda, ddxdda = run_sim_var(a_0)\n", "x_1st_order = np.zeros(N)\n", "x_2nd_order = np.zeros(N)\n", "for i,a in enumerate(a_grid):\n", " x_1st_order[i] = x + (a-a_0)*dxda\n", " x_2nd_order[i] = x + (a-a_0)*dxda + 0.5*(a-a_0)*(a-a_0)*ddxdda" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we plot the exact final position that we obtained from running a full $N$-body simulation as well as our approximation near a neighbourhood of $a_0$ which we got from the variational equations. " ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEKCAYAAADaa8itAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XmcTfX/wPHXe5iFISS7rNnJUpYKjSQSiagkS7IktCeK\nQip9S6Ks2UkqoSTKNrL+yL6Tbeyyz4yZMXPn/fvjXtOkwb0z995z753P8/E4j7n33HPO5z0t9z2f\nXVQVwzAMw3CHIKsDMAzDMAKHSSqGYRiG25ikYhiGYbiNSSqGYRiG25ikYhiGYbiNSSqGYRiG22S1\nOgBPEBEzTtowDMNFqioZfUbA1lRUNSCP999/3/IYzO9nfj/z+1l/JCcn02leJ2Ztn+WW57lLwCYV\nwzCMQDbmzzFsPLGRZmWbWR3Kv5ikYlhCFXbvtv80Asvq1fDbb1ZHEdhWR61m0IpBzH16LuEh4VaH\n8y8mqfiZiIgIq0Nwm44dITr63+cC6fdLS2b5/YIC9JvFF/79nYg+wVOzn2Jyi8mUvr201eH8h7iz\nLc1XiIgG4u9lGIbxzOxnqJSvEgMeHODW54oI6oaOepNUDMMw/MiFuAvkCstFkLi3OuiupBKglVTD\nH1y6ZO9XMQJHbCz07291FIEtT7Y8bk8o7uS7kRkBb9s2mDDB6igMd0pKgpIlrY7CsJJp/jIMwzBM\n85dhGEagS0pOYsyGMSQlJ1kditO8mlREZKKInBaRbTe5ZqSI7BeRLSJSLdX5JiKyR0T2icjb3onY\n8LTt2+HkSaujMNzlww/hr7+sjiJwvLP0HebsmYOQ4QqE13i7pjIZaHyjD0XkUaC0qpYBugNjHeeD\ngK8c91YC2opIec+Ha3jazz+bzvpAUrky5M5tdRSB4fud3/PDrh+Y9eQssgRlsTocpzndpyIibYBF\nqhotIv2BGsAQVd3kUoEixYH5qnp3Gp+NBZar6neO97uBCKAk8L6qPuo43xdQVf3kBmWYPhXDMG7q\n0qVL7Ny5k6ioKI4ePUp0dDTx8fHYbDbCw8PJkSMHBQoUoHjx4pQsWZJixYoh4p0aw44zO2gwtQG/\nP/c71QtV90qZ7upTcWWV4gGq+oOI1AUeBj4FxgC1MxpEKkWAo6neH3OcS+t8LTeWaxhGgIuNjeW3\n337jl19+Ye3atezZs8el+3Pnzk21atWoXr06devWpUGDBuTJk8ftcV6Mv0jL71ry+SOfey2huJMr\nScXm+PkYMF5VF4jIEA/ElJr/NCQa6ZKcDHPnQqtW4KU/Ag0PmTcPLl+GDh2sjuQfqsrKlSsZPXo0\nP/30E/Hx8SmfhYSEULlyZUqWLMmdd95Jnjx5CAsLIygoiCtXrhAdHc2JEyeIiopi3759nDlzhsjI\nSCIjIxk+fDgiQo0aNWjSpAmtWrWievXqbqnJqCpv3f8W7au2z/CzrOBKUjkuIuOARsAnIhKK+/tk\njgN3pnpf1HEuBCiWxvkbGjhwYMrriIgIn1izx/ivoCD46Sdo0gTCfWtdPMNFFStCqu9sS6kqc+fO\nZeDAgWzfvj3lfJ06dWjRogUNGzakatWqhISEOP3MkydPsnnzZv7880+WL1/OmjVr2LhxIxs3buTD\nDz+kZMmStGrVijZt2lCrVq10J5g82fLQ7Z5u6brXFdcSpNu5sNZ+dqAVUMbxvhDwSDrW7C8BbL/B\nZ02BBY7XdYB1jtdZgL+A4tgTzBagwk3KUMMwMqc//vhDa9asqYACWqBAAR0wYIBGRUW5tZzY2Fhd\ntGiR9ujRQwsWLJhSHqDlypXTjz/+WI8dO+bWMj3J8b2Z8X1ZnL4QPnHm3C2eMRM4ASQAUcDz2Ed5\ndUt1zVeOBLIVqJHqfBNgL7Af6HuLctz5z9owDD9w8eJF7d69e8oXe8GCBXXUqFGakJDg8bKTkpJ0\n5cqV+uqrr/4rwQQFBWmTJk101qxZGh8f7/E4MsKKpLIpjXPb3BGEuw+TVPzL0aOqCxdaHYWRESdP\nqrZta135kZGRWqRIEQU0ODhY33vvPY2JibEklsTERP3ll1+0devWGhwcnJJg8ubNq6+99pru3Lnz\nX9fHJMRoQpLnE9+teC2pAD2A7UAssC3VcQiY4Y4g3H2YpOJfduxQ/ewzq6MwMiImRnXVKu+Xa7PZ\ndOjQoRoUFKSA1q5dW3fs2OH9QG7g7NmzOnLkSK1ateq/msfuv/9+nTRpkkZHR+tTPzylH/3xkdWh\nui2p3HKeiojkAvIAHwN9U30Urarnb3qzRcw8FcMIfHFxcTz33HPMmTMHgH79+jF48GCyZnVl/JF3\nqCobN25kwoQJzJw5k2jH7nShDUK57b7bmPv4XO6vdb/X5sGkxeynchMmqRhGYDt//jyPP/44q1ev\nJleuXEyfPp3mzZtbHZZTYmJi+OGHH/hszmfsKrcLJgCXoGrVqnTt2pV27dqR24JlCby2oKSIrHL8\njBaRy46f147LGQ3AMAAWL4atW62Owkivdu1g717vlHX8+HHq1q3L6tWrKVq0KKtXr/abhAKQI0cO\nHmr5EOcePMekxybxWufXyJs3L1u3bqVXr14UKlSIDh068Mcff+CPfxybmorhE2bPhsKF4f77rY7E\nSI/t2+GuuyBbNs+Wc/LkSR588EH2799P5cqVWbhwIUWLFvVsoR7Qb0k/7sh+B2/c/wYACQkJzJs3\njwkTJrBkyZKU68qWLUuXLl3o0KEDBQoU8GhMXm/+EpEw4CWgLvbOppXAWFX1kelO/zBJxTACz+nT\np4mIiGDPnj1Uq1aNpUuXcvvtt1sdVrokazKCpNmHcvDgQSZNmsTkyZM5ceIEAFmzZqVFixZ06NCB\nRx55hLCwMLfHZEVS+R6IBmY4Tj0L5FbVNhkNwt1MUjGMwHL58mXq16/P1q1bqVy5MsuXL+eOO+6w\nOiyPSkpKYuHChUyYMIEFCxZgs9lXygoPD6dp06a0atWKpk2bctttt7mlPCuSyi5VrXirc77AJBX/\n9NFH8Nprnm9CMdxrxAj7cju9e3vm+YmJiTRv3pzffvuNsmXLsnLlSvLnz++ZwnzUiRMnmDp1KrNn\nz2bTpn8Whs+SJQt16tTh4YcfplGjRtSqVYvg4OB0lWFFUpkBfKWq6xzvawM9VdWHlo+zM0nFP40Y\nAZ06Qa5cVkdiuOLiRbh6FTzxPa+qdO/ena+//po77riDdevWUbp0afcX5EcOHz7MvHnzmDNnDmvW\nrEmpwYC9FnPPPfdQq1YtatWqRc2aNSlWrBhBQbdeptFrSUVEtmPvQwkGymFfXkWxr8O1x9RUDMPw\nlBEjRvDqq68SGhrK8uXLue+++6wOyWVXbVd5ddGrDHloCLdnc28f0KVLl1ixYgWLFy9myZIlaS7n\nny1bNsqUKUPZsmUpW7YsRYoUoUCBAuTPn58CBQqQK1cusmfPzm233ea1pFIc+9DjO4Ej13+uqv85\nZzWTVAzD/61evZqIiAiSkpKYNWsWTz/9tNUhpUuvX3tx9PJR5j49lyDx7Ga7f//9Nxs2bGD9+vWs\nX7+ejRs3cubMGafv93bz13ZVrZLRAr3BJBX/dOiQfW+V11+3OhLDWX/+Cf37w6JF7n3u6dOnqVGj\nBidOnOD1119n2LBh7i3AS6ZumcqHKz9kQ9cN5Aqzpl334sWL7N+/n3379rF//35OnTrF6dOnU47o\n6GiuXLlCbGys15PKVOx9KhsyWqinmaTin86cgSVL4NlnrY7EcFZyMpw/D+4ciGWz2WjUqBHLly+n\nXr16LF26NN2dz1badHITjWc0JrJjJJXyV7I6nFuyoqN+D3AX9iawWOy7Mqqmsde81UxSMQz/9emn\nn9KnTx/y58/Pli1bKFSokNUhuexS/CWqjq3KZ498RuuKra0OxylWJJXiaZ03fSqGkXnZbJAli/ue\nt3XrVmrWrEliYiILFiygadOm7nu4F6kqfxz5gwdLPGh1KE7z2tpf1ziSR26guePI7YsJxfBvP/0E\nP/xgdRSGM1ThzjvhwgX3PC8uLo527dqRmJhIjx49/DahgP0L2p8Sijs5nVRE5BXgGyC/45ghIh6a\n7mRkViVLQpkyVkdhOEMEDhyAPHnc87yBAweyc+dOypUrx2effeaehxpe50rz1zbgPlWNdbwPB9aa\nPhXDMDJq06ZN1KpVi+TkZNauXUvt2rWtDinT8XrzF/aOeVuq9zbHOcMwMqGYGHsTWEYlJSXRpUsX\nbDYbr7zyil8mlEvxlzh04ZDVYfgEV5LKZOD/RGSgiAwE1gETPRKVkam99Rbs2GF1FMatvPoqfPdd\nxp/zxRdfsHnzZooXL84HH3yQ8Qd6WbIm035ue8b8OcbqUHyCS/upiEgN7EvfA6xU1c0eiSqDTPOX\nf1uzBsqVg7x5rY7EuJXkZPtikul19OhRypcvz5UrV/j111959NFH3ReclwxeMZjfD/zOso7LCMkS\nYnU46eau5i+XNnNW1U3AplteeBMi0gT4AnstaaKqfnLd57mBSUBpIA7orKq7HJ+9BrwAJAPbgedV\n9WpG4jF8j9moy39kJKEA9OnThytXrtCmTRu/TCgL9i1g/MbxbOi6wa8Tijt5dedHEQkC9gENgRPA\nBuAZVd2T6pr/AdGq+oGIlANGqerDIlIYWAWUV9WrIvIdsEBVp6VRjqmpGIYHxcTYaykZ2cpj5cqV\n1K9fn7CwMPbs2UPx4mlOhfNZf53/i/sn3s+8Z+Zx/53+/5eQFR317lAL2K+qR1Q1EZgFtLjumorA\nMgBV3QuUEJF8js+yAOEikhXIjj0xGQEmLg4aNnRPJ7DhGYsXQ79+6b/fZrPx8ssvA/D222/7XUIB\n2HB8Ax80+CAgEoo7udT85QZFgKOp3h/DnmhS2wq0AlaLSC2gGFBUVTeLyDDsS+9fAX5X1SUYASdb\nNhgyxJ5U0tht1fABLVvaj/SaOHEiW7Zs4c4776RPnz7uC8yL2lZpa3UIPsmppCIi5bHXKIo4Th0H\nflbV3R6IaSgwQkQ2Ye832QzYHH0tLbDv43IJmC0iz6rqzLQeMnDgwJTXERERREREeCBUw1P8cNsM\nw0kXL17k3XffBezrfGXPnt3iiDKnyMhIIiMj3f5cZ/ZTeRtoi72p6pjjdFHgGWCWqg51ujCROsBA\nVW3ieN8X+6KUn9zknoPA3UAToLGqdnWcbw/UVtVeadxj+lQMw0OSk2HPHqiYzu35+vXrx9ChQ6lX\nrx4rVqxATHXUJ3hz58d9QCVHH0jq8yHATlV1elENEckC7MXeUX8SWA+0TV3jEZFcwBVVTRSRrsAD\nqtrJ0RQ2EagJJGCfN7NBVUelUY5JKn5uxQqYORPGjbM6EuN6p05Bu3awdKnr9548eZLSpUsTFxfH\nunXr/Gqio6oGdAL0Zkd9MlA4jfOFHJ85TVVtQC/gd2An9prObhHpLiLdHJdVAHaIyG6gMfCK4971\nwGzszWFbsc/mH+9K+Yb/qF4dBgywOgojLQULpi+hAAwZMoS4uDieeOIJv0oo8UnxNJreiL/O/2V1\nKD7PmZpKE+ArYD//dLIXw763Sm9VXejRCNPB1FQMw/ccPHiQcuXKYbPZ2L59O5Uq+f7GVWCvobzw\n8wvEJsYy68lZAVtb8drkR1VdJCJlsY/SSt1Rv8FR8zAMIxP5v/+Du++2j9Jzxfvvv09SUhIdOnTw\nm4QCMPbPsWw4sYG1L6wN2ITiThma/Cgiz6vqZDfG4xamphIYPvjAPrnulVesjsRI7amn4MsvoUAB\n5+/Zvn07VatWJWvWrOzbt48SJUp4LD53WnN0DU/MeoLVnVdTJm9g78ng9Z0fbxBElKoWy2gQ7maS\nSmC4cAHCwyHErH7h95588knmzJlDr169+PLLL60OxymJtkTKjyrPiCYjaFa2mdXheJw3R39tu9FH\nQFlVDc1oEO5mkoph+I6dO3dSuXJlQkNDOXTokF/tOX8y+iSFcvpPvBnhzQUlC2AfhXX9pqECrMlo\nAIZxM8mO8YUZXbjQcI8NG+zNXsVcaJ/46KOPAOjSpYtfJRQg0yQUd3Lmf9VfgByO9bpSH4eBSI9G\nZ2R6DRvav8gM37BmDRw86Pz1f/31F7NmzSJr1qx+uxyL4RqvrlLsLab5K3DEx0NYmNVRGOnVpUsX\nJk6cSOfOnZk40ezp58t8oqPeV5mkYhjWi4qKonTp0iQnJ7Nnzx7KlPHt0VNnr5xlddRqWpS/fuH0\nzMFfl743DJddumR1BAbA7t3wxx/OX/+///2PpKQknn76aZ9PKEnJSbT9sS1rjppu4owyScXwaYmJ\nUL48XDX7e1ru77/h8GHnrj1z5gwTJkwA4J133vFcUG7Sf1l/AD5s+KHFkfg/Z5e+F+x7mhy95cWG\n4UbBwXDihNlXxRfUr28/nDF69GgSEhJo3rw5lStX9mxgGTR712xm7ZjFn93+JGuQt7eYCjxO96mI\nyHZVreLheNzC9KkYhnXi4uIoVqwYZ8+eJTIykgcffNDqkG5o55mdREyN4LfnfqNGoRpWh2MpK/pU\nNolIzYwWaBiuSkyEo6aObKmLF2HMGOeunT59OmfPnqVGjRrUd7ZqY5GrtquMbjo60ycUd3IlqdQG\n1onIARHZJiLbbzLb3jDcZts2ePVVq6PI3OLiICHh1tclJyczfPhwAN544w2fX4CxeqHqtKnUxuow\nAoorzV/F0zqvqkfcGpEbmOYvw7DGggULaNasGUWLFuXgwYMEBwdbHZLhJCuav6KAekBHRyJR7Eu4\nGIZhADBs2DAAXnnlFZNQMilXkspo4D7s+9UDRAP/2crXMDzh77/h0CGro8icVKFvX/vqBjezefNm\nli9fTo4cOejatat3gnNR7NVYq0MIeC71qahqTyAeQFUvAGZRcsMrfv8d5s2zOorMyWaDwoUh9Bbr\nkX/11VcAvPDCC+TKlcsLkbkm9mos9028j1VRq6wOJaC50qfyf8D92Hd8rCEi+YDfVbW6JwNMD9On\nYhjedf78eYoUKUJ8fDz79u3zuRn0qkrbH9uSLTgbkx6f5PMDCKzgzaXvrxkJzAXyi8iHQGtgQEYD\nMAzD/02ZMoX4+HgeeeQRn0soAMPXDeev83+x8vmVJqF4mEsLSopIeaAh9r1Ulqrqbk8FlhGmphKY\ntm2DXLmgeJrjEA1PGTAAOnWC0qXT/jw5OZmyZcty4MABfvrpJx5//HGvxncryw8t59k5z7LuhXUU\nz23+47kRr4/+EpFPVHWPqo5S1a9UdbeIfOJqgSLSRET2iMg+EXk7jc9zi8gcEdkqIutEpGKqz3KJ\nyA8isltEdopIbVfLN/xXZCTs2WN1FJlPrVqQN++NP//99985cOAAxYoV47HHHvNeYE5QVd6LfI9v\nWn1jEoqXuNKnsklVa1x3bpuq3u10YSJBwD7stZ0TwAbgGVXdk+qa/wHRqvqBiJQDRqnqw47PpgAr\nVHWyiGQFsqvq5TTKMTUVw/CSxx9/nPnz5/PRRx/Rr18/q8P5j0RbIsFZzPDmW/FaTUVEeojIdqCc\nYyb9teMQ4OqM+lrAfsfOkYnALOD6zQsqAssAVHUvUEJE8onIbUA9VZ3s+CwprYRiGIb3HD58mF9+\n+YWQkBBeeOEFq8NJk0ko3uVM81dToBmQBWie6rhHVZ9zsbwiQOpVnI45zqW2FWgFICK1gGJAUaAk\ncFZEJovIJhEZLyLZXCzf8HPz5sHZs1ZHkXl8+in89tuNPx87diyqSps2bcifP7/3AjN8ljOjv0oD\nicBe4DL2TnoAROR2VT3v5piGAiNEZBOwHdgM2IBgoAbQU1X/FJEvgL7A+2k9ZODAgSmvIyIiiIiI\ncHOYhhV27YJKleCOO6yOJHNo2tQ+OCIt8fHxKVsE9+zZ04tRGe4QGRlJZGSk2597yz4VEXkZ6IG9\npnCCVEkFUFUt5XRhInWAgaraxPG+r+MZN+zwdzSzVQHCgbXXyhORusDbqto8jXtMn4pheNj06dPp\n0KED1apVY9OmTT4xVPdUzCm+/L8vGfLQEJ+Ix594rU9FVUeqagVgsqqWUtWSqQ6nE4rDBuAuESku\nIiHAM8DPqS9wjPAKdrzuir1jPkZVTwNHRaSs49KGwC4XyzcMw03Gjx8PwEsvveQTX+BXbVdp80Mb\nQrKE+EQ8mZWr81TyAGWAsGvnVNWFXavtQ4qBEdgT2kRVHSoi3e2P0vGO2sxUIBnYCbygqpcc91YF\nJmBvCjsIPH/ts+vKMDWVAKUK//sfvPkmZMlidTSBbfp0OHUK3nrrv5/t2bOHChUqkCNHDk6ePEmO\nHDm8H+B1ev/am8OXDvPTMz8RJGandFd5fUa9iHQBXsHeab4FqAOsBR5ypUBVXQSUu+7cuFSv113/\nearPtgJmo7BMTMS+aVdcHPjA91hAe+wxiIlJ+7NrfSnPPPOMTySUaVunsejAIjZ03WASisVc2k4Y\n+xf6OlWt5phd/5GqtvJkgOlhaiqG4TlXr16laNGi/P3336xbt47ata2dg7z11FYenv4wyzsup3L+\nypbG4s+s2E8lXlXjHYWHOiYsplmjMAwjcM2fP5+///6bypUrU6tWLavDoUTuEvz41I8mofgIV5LK\nMRHJDcwDFovIT4DP7fpoBL5z52DwYKujCGwLFtjX+0rLhAkTAOjSpYtPdIjnCstF/eL1rQ7DcHCp\noz7lJpEHgVzAIlW96vaoMsg0fwW2hASYOBFeesnqSAKXzQYXLvx3PlBUVBQlSpQgODiYEydOkPdm\ni4IZfsWKpe9TqOqKjBZsGOkVGmoSiqdlyZL2BNPJkyejqrRq1cokFCNNZpiEYRj/kpxsH2F3PZvN\nxqRJkwB705dVjlw8QnzSLfY2Nixjkorhl/btg3fftTqKwLR7NzzwwH/PL1myhKioKEqWLEmDBg28\nHxhwKf4Sj8x4hIX7F1pSvnFrLu2n4sw5w/CG/PnBLOfmGZUqwao0tnG/1kH/wgsvEBTk/b9HkzWZ\nDvM60LBkQ1pWaOn18g3neHU/FW8xHfWG4V7nz5+nUKFCJCYmEhUVRdGiRb0ew5A/hrDwr4Us77ic\nkCwhXi8/0Hmto15EegAvAaVEJPX+KTmB1RkNwDAM33L6NBQo8O9z33//PVevXqVRo0aWJJSF+xcy\n9s+xrO+63iQUH+dMHXYm9v1Tfibj+6kYhtusWAE+uNGgX7t0CerVs6+xltr06dMB6NChgwVRwYL9\nC/iu9XcUzlnYkvIN56VrnoqvM81fmcPZs/ajfHmrIwlsf/31F2XKlCE8PJzTp08THh5udUiGB3h9\nmRaxe05E3nO8L+bYmdEwLHHHHSaheMO1Wkrr1q1NQjFuyZUhHKOB+4C2jvfRwCi3R2QYhmW2bYOk\npH/eJycnM23aNADat29vUVSGP3FlRn1tVa0hIpsBVPWCY6Mtw3CKqrJr1y5WrFjBjh07OH/+PCEh\nIRQtWpT77ruPBg0auLyM+rRpsH8/fPCBh4LORJKToWdPWLQIsjq+GVavXs3hw4cpWrSo2ZLbcIor\nSSVRRLIACiAi+bBvpGUYN2Wz2Zg+fTojR45k8+bNN7wuLCyMjh078sYbb1CmTBmnnt2smdmsy12C\ngmDlyn+fu9b09dxzz5HFS/+gj146SrdfuvHzMz8TnCXYK2Ua7uPKPJV2wNPAPcAUoDXQX1V/8Fh0\n6WQ66n3H2rVr6dGjB1u3bgXg9ttv59FHH6VmzZrky5ePhIQEDhw4wNKlS1m3bh0AwcHBvP3227z7\n7ruEhYXd7PGGB8XFxVGoUCEuXbrEzp07qVixosfLjE+Kp97kejxV8SneeiCNLScNj3FXRz2q6vQB\nlAd6Oo4KrtzrzcP+axlWstls+uGHH2qWLFkU0GLFiumUKVM0Li7uhvfs3r1bO3bsqNhrw1q1alU9\nePCgU+UlJror8szr999VY2L+ef/dd98poPfee69Xyk9OTtbO8zprm+/baHJyslfKNP7h+N7M8Pev\nK6O/QoEa2Je8zwu0uTYSzDBSS0hIoG3btrz77rvYbDbefPNN9uzZQ8eOHW9a8yhfvjxTpkxh1apV\nlC5dmq1bt3Lvvfeydu3am5b34Yfw6afu/i0yn2nT/t1Jf62D3ltzU8ZvHM+64+uY1GKST+zTYqSP\nK81fi4BLwEbAdu28qg7zTGjpZ5q/rBMbG0vz5s1Zvnw5t912G7NmzeLRRx91+TkXL16kXbt2/Prr\nr4SHhzN//vwbLmIYFwdhYfb96w33OH36NEWKFEFEOHHiBPny5fNoefvO7aPupLqs6ryKsnnLerQs\nI21eb/4CdrijagQ0AfYA+4C30/g8NzAH2AqsAype93kQsAn4+SZluKEyaLgqLi5OH374YQW0UKFC\numXLlgw9LzExUdu3b6+AZsuWTdeuXeumSI1bGT58uALavHlzr5SXnJysB88719RpeAbebv4C1ohI\nlYwkMBEJAr4CGgOVgLYicv30tXeAzapaFegIjLzu81eAXRmJw3C/5ORknnvuOZYsWUKBAgWIjIyk\natWqGXpm1qxZmTJlCp06dSIuLo5mzZqxb9++NK+Nj4fLlzNUXKb27bf21Qmu8fayLCJCyTwlvVKW\n4VmuJJW6wEYR2Ssi20Rk+3ULTDqjFrBfVY+oaiIwC2hx3TUVgWUAqroXKOEYvoyIFAWaAhNcLNfw\nsMGDB/Pjjz+SK1cuFi9eTNmy7mnCCAoKYvz48TRt2pRz587RvHlzLqeRPQYPhu++c0uRmdKePf+8\n3rFjB5s2bSJ37tw0a9bMuqAMv+TKPBXXG8b/qwhwNNX7Y9gTTWpbgVbAascyMMWAosDfwHDgLeyD\nBQwfMWfOHAYNGkRQUBCzZs2iSpUMVWj/Izg4mO+//57777+fbdu20alTJ3788cd/deZ++KHpU8mI\nQYP+eX2tlvL000+bId2Gy5yuqajqEeAyUAAonupwt6FAHhHZhH3o8mbAJiKPAadVdQsgjsOw2J49\ne1KaSD755BOaNGnikXLCw8NTakJz585l5Mh/t4qahOIeNpuNGTNmAJ5t+tp2eht7z+712PMN6zhd\nUxGRLtj7M4oCW4A6wFrgIRfKO4695nFNUce5FKoaDXROVe5B4CDwDPC4iDQFsgE5RWSaqqb5X/7A\ngQNTXke3nLuAAAAgAElEQVRERJglJjzg6tWrtGvXjtjYWJ599lneeOMNj5Z31113MWXKFFq2bEnf\nvn1p3Lgx5VOtKHn0KISE/HcvEOPmvvrKvjJBiRKwbNkyTpw4QenSpbnvvvs8Ut65K+doMasF/3v4\nf5S7o5xHyjBuLTIyksjISPc/2NkefWA7EAZs0X8mQs5xZVQAkAX4C3sNJwR7cqpw3TW5gGDH667A\nlDSe8yBm9Jfl+vXrp4CWLFlSL1265LVyO3XqpIDWrFlTE1PNevzgA9XZs70WRsCYPFn11Cn762uj\n7QYOHOiRspJsSfrI9Ef0zd/e9MjzjfTDTaO/XEkIGxw/twChjtc7XS7QPqR4L7Af6Os41x3o5nhd\nx/H5bmA2kCuNZ5ikYrE//vhDRUSDgoJ01apVXi374sWLWqxYMQV08ODBXi07kEVHR2v27NkV0AMH\nDnikjL6L+2rDqQ010WaWQPA1ViSVudjnkAwE/gB+An51RxDuPkxS8awrV65oqVKlFNB33nnHkhiW\nLl2qgAYHB+vevXstiSHQTJ06VQGtW7euR54/e+dsLT68uJ6JOeOR5xsZ466k4kpHfUtVvaiqA4EB\nwETgCWfvNwLHkCFDOHjwIJUrV+b999+3JIaHHnqIzp07k5iYSO/eva/9McGff8KRI5aE5Jf694fd\nu+2vPb0sy+WEy8x+ajb5wj07O9+wltlO2HDJzp07qVatGklJSaxZs8ZjnbnO+PvvvylbtiwXL17k\nxx9/pFWrVowYAVWrghmX4ZzFi+HeeyEm5ijFixcnJCSEU6dOkTt3bqtDM7zMa9sJi0i0iFx2HNGp\n3keLiJnDnIkkJyfTvXt3kpKSePHFFy1NKAD58uXjww8/BOC1114jNjaWV14xCcUVjRpBnjzwzTff\noKo8/vjjJqEYGXLLpKKqOVX1NseRM9X7nKp6mzeCNHzDtGnTWL16NQUKFODjjz+2OhwAunfvTvXq\n1YmKiuJTs1Rxuqiq15dlMQKXK6sUhwEvYV+uRYGVwFhVjfdceOljmr/cLyYmhjJlynDq1CmmTZvm\nU/uVr1y5kvr16xMeHs6BAwfYurUA+fNDtWpWR+bb2reHt9+GhISN3HvvveTLl4/jx48THOye3RYT\nkhIIzRrqlmcZnue15q9UpmFfBPJL7ItCVgKmZzQAwz8MHTqUU6dOUatWLdq1a2d1OP9Sr149mjdv\nTmxsLB988AExMfbl8I2b69cPSpX6p4P+2WefdVtCOXD+ABVHV+RS/CW3PM/wH67UVHapasVbnfMF\npqbiXkeOHKFcuXIkJCRY3jl/Izt27KBq1aoEBQWxa9cup/e4z+wSExMpXLgwZ8+eZePGjdSoUSPD\nz4y9Gst9E++j2z3d6FWrlxuiNLzBiprKJhGpkyqA2sCfGQ3A8H19+/YlISGBZ5991icTCkDlypXp\n2LEjSUlJvPvuu1aH4/Ou/c21aNEizp49S6VKlahevbobnqt0nd+V6oWq07Nmzww/z/A/riSVe7Dv\nqXJYRA5jX/erZjqXwDf8xNatW5k1axahoaE+0zl/I4MGDSI0NJQffviBYcOOsGiR1RH5rgcfhB07\n/lmRuH379m7ZwveLdV+w99xexj421mwJnEm5klSaACWxL5HyoON1E6AZ0Nz9oRm+4NrCnC+++CLF\nihW7+cUWu/POO+nWrRsACxeOp3BhiwPyYfPnQ4ECF/j5558REbf0k529cpaR60cy56k5ZAvO5oYo\nDX9kJj8aN7Rxo31UULZs2Th48CAFCxa0OqRbOnbsGKVLlyYxMZEdO3ZQsaLPdfn5jPHjx9O9e3ce\nfvhhFi9e7JZnxifFE5bV7MHij6zoUzEymWtLsPTs2dMvEgpA0aJFeeGFF1DVlImRxr+dP2//OXXq\nVAC3Dg83CcVwqqYi9sbRoqp69JYX+wBTU8m4devWcd999xEeHs6hQ4fIl89/1muKiorirrvuIimp\nLYMHf0L//v6REL1BFapUgXHjDlG3binCw8M5deoUOXLksDo0w2Jerak4vqF/zWhhhv+4Vkt5+eWX\n/SqhABQrVoxOnTqhupLNm/9ndTg+RQS2b4fffpsMQOvWrU1CMdzKlXkqU4GvVHWDZ0PKOFNTyZhV\nq1ZRr149cubMyaFDh8ibN6/VIbns0KFDKXNV9uzZw1133WVxRL4jOTmZUqVKceTIEZYtW0aDBg3S\n9Zx1x9YRdSmKpyo95eYIDStY0adSG1grIgdEZJsZShy43nvvPcC+SKM/JhSAkiVL0r59e2w2G599\nNszqcHzGhg2wfPlKjhw5QrFixXjwwQfT9ZxTMado80Mb04di/IcrNZXiaZ1XVZ/bvcLUVNJv+fLl\nPPTQQ+TOnZtDhw759Yq1u3btolKlgWTJ8hAnTrQif/78VodkqcREaNwY7ryzG9Omfc27777LkCFD\nXH+OLZGG0xrSoEQDBjUY5IFIDSt4vaaiqkfSOjIagOE7VDWllvLGG2/4dUIBqFixIo0bKzbbG4wa\nNcrqcCwXHAy//HKFOXO+BdK/IvGbv79JztCcvB9hzQZthm9zOqmI3XMi8p7jfTERqeW50AxvW7Jk\nCatWreL222/n5Zdftjoct+jXrxdwhVGjRnHlyhWrw7Hc3LlziYmJoU6dOpQtW9bl+2dsm8GC/QuY\n0XIGQWJmJBj/5cp/FaOB+4C2jvfRgPnzL0CoKgMGDACgT58+3HZbYGyVU79+fe69tybnziUxefJk\nq8OxTGIiTJ0KU6bY56Z07NgxXc+5u8DdzHtmHnmy5XFneEYAcamjXlV7AvEAqnoBCPFIVIbXLVy4\nkP/7v/8jX7589OwZOAsBigjNmn0CTOPzzz/HZrNZHZIlzp+HDRuiWbp0CSEhITz1VPpGbN1d4G4q\n56/s5uiMQOJKUkkUkSzYN+hCRPIBya4WKCJNRGSPiOwTkbfT+Dy3iMwRka0isk5EKjrOFxWRZSKy\n0zHyLDDaZ3xA6r6Uvn37Bty8hXfeqU/Jkq9x8OBB5syZY3U4lihQAO68c3TKlsG333671SEZAcqV\npDISmAsUEJEPgVXAR64UJiJB2Df4aox9k6+2IlL+usveATaralWgo6NcgCTgdVWthL0Zrmca9xrp\n8PPPP7Nx40YKFizIiy++aHU4bhccnIU333wDgM8++4zMODJQVVOWZUlv05dhOMOV0V/fAH2wJ5IT\nwBOq+oOL5dUC9jtGjiUCs4AW111TEVjmKHMvUEJE8qnqKVXd4jgfA+wGirhYvnGd5OTklFpKv379\nyJ49u8UReUa7dp247baGrF+/njVr1lgdjlft3Qtvvx3F7t27yZcvH40bN3b63qhLUR6MzAhEroz+\nCgVqALmAvECbayPBXFAESL1+2DH+mxi2Aq0cZdYCigFFr4ulBFAN+D8XyzeuM2fOHLZt20aRIkVS\nlo0PRKrZueOO0QAMG5a5JkOGhsKOHb8B0K5dO6e3DN79927uGX8Phy8e9mB0RqDJ6sK1PwGXgI1A\ngmfCAWAoMEJENgHbgc1ASu+qiOQAZgOvOGosabq2DwhAREQEERERHgrXf9lstpQ1vt59913CwgJ3\ndnTu3LBqVU5KlAhh3rx5HDhwgNKlS1sdllcUKBDH2rX27svnn3/eqXsuJ1ym5Xct+eThTyiRu4QH\nozOsEhkZSWRkpPsfrKpOHcAOZ6+9yTPqAItSve8LvH2Lew4BORyvswKLsCeUm92jxq3NnDlTAS1W\nrJgmJCRYHY5XdOrUSQHt1auX1aF4zYwZMxTQmjVrOnW9LdmmT8x6Ql+c/6KHIzN8ieN7M0Pf8arq\nUkf9GhGpksEctgG4S0SKi0gI8Azwc+oLRCSXiAQ7XncFVug/NZJJwC5VHZHBODK9pKSklNrcgAED\nCAnJHKPD77zzAyCUSZMmceHCBavD8bivv4YPPjgGQJcuXZy6Z+iqoZyKOcWIR83/ZobrXEkqdYFN\nIrI3vQtKqqoN6AX8DuwEZqnqbhHpLiLXGvQrADtEZDf2UWKvAIjIA0A74CER2Swim0SkiSvlG/+Y\nOXMm+/bto1SpUplqNFDWrEWpX/9xrly5wrhx46wOx+PKlz/I3r1jyZ49O88888wtr49PimfpoaXM\nbjObkCyZ4w8Nw73MgpKZUGJiIhUqVODAgQNMmTIlUyUVgEWLFvHoo49SuHBhDh06FNC1tH79+jF0\n6FCef/55Jk2aZHU4hg9z14KSriSVUOBJoASpOvhVdXBGg3A3k1RubuLEiXTp0oWyZcuyc+dOsmZ1\nZbyG/1NVqlSpws6dO5k2bZpbt9P1JfHxSZQseSenTp1i1apVPPDAA1aHZPgwK/ZT+Qn7nJIkIDbV\nYfiRhIQEBg2yL1f+/vvvZ7qEAjBhgtC0qX1HyM8//zxgJ0M+8MA5Tp0qSfny5bn//vutDsfIJFz5\nRimqqqYPw8+NGzeOo0ePUqVKFafa2ANRxYrwwAMNmTo1P1u2bEnZQybQFCzYE9hAly5DEcnwH6CG\n4RRvj/4yLBQTE8OHH34IwJAhQwgKypxLlz/wAFSsGEqvXr0Ae20l0Bw/fpxFi+YSHCw3bd6LPBzJ\nW7+/5cXIjEDn6uivjRkZ/WVYa+TIkZw5c4batWvTvHlzq8Ox3IsvvkhYWBgLFixg9+7dVofjVp9+\n+gvJycm0aNHihjteHrt8jLY/tuWR0o94OTojkJnRX5nEhQsXKFmyJJcuXWLp0qUB2dzjihUrYNIk\nyJbtRcaNG0e3bt0CZojx2bM2ihbdTEJCbRYuXECTJv9ttU5ISqD+lPq0Kt+Kt+v+Z7FwIxPy+ugv\nf2KSyn+98847fPzxxzRs2JAlS5ZYHY7lYmMhJgYuXNhDhQoVCAsLIyoqinz58lkdWob98ssvNG/e\nnFKlSrF///40mzm7ze/G+bjz/NDmB9PfYgBeHP0lIqscP6NF5HKqI1pELmc0AMPzjh8/zhdffAGQ\n0qeS2YWH2/cYKV++PM2aNSM+Pp4xY8ZYHZZbfPXVVwD06NEjzYQya8csVh9dzeQWk01CMdzO1FQy\ngU6dOjF16lSefPJJZs+ebXU4PuXMGdi50z76K3/+/Bw5csSvF9b89dcjPPZYB8LC1nP8+PE0N+OK\nuRrDuSvnKJ47zRZtI5OyYp6K4Yc2btzI1KlTCQkJ4ZNPPrE6HJ9is0FEBFSpEkH16tU5c+YM33zz\njdVhZciUKfOBIjz77LM33N0xR0gOk1AMjzFJJYCpKm+8Yd/x8OWXX840S707K0sW2LED7rhDeP31\n1wH/ngwZGxvL4sUDgG/p2bOn1eEYmZRJKgFs7ty5rFixgrx58/Luu+9aHY5Putbl8NRTT1GkSBF2\n7drFb7/9Zm1Q6TRz5kwuXrxInTp1qFGjhtXhGJmUMx310x0/X/F8OIa7xMTE8OqrrwIwePBgcufO\nbXFEvmvHDti7N4TevXsD/rkzpKrSr18uoGjKpM5rVkWt4kriFWsCMzIdZ2oq94hIYaCziOQRkdtT\nH54O0EifQYMGcfToUe655x66d+9udTg+bfdu2LcPunXrRnh4OEuWLGHbNv+a17tixR+cOzeFO+5I\nonXr1innt5zaQsvvWnLowiELozMyE2eSylhgKVAe+1bCqY8/PReakV7btm1j+PDhiAhjx44lS5Ys\nVofk09q0gSefhDx58tC5c2fA/5Zu+fzzYcBCevV6kdDQUADOx52n1Xet+PLRL6mUv5K1ARqZhisz\n6seoag8Px+MWmXlIsc1mo169eqxdu5aePXumzFkwnHPw4EHKlClDlixZOHLkCIUKFbI6pFvasWMP\nVapUJCwsNGUCpy3ZRtOZTamcrzLDGvtfc57hfV4fUqyqPUSkqoj0chx3Z7Rww/0+++wz1q5dS8GC\nBRkyZIjV4fiNhATo2BGKFi1Fy5YtSUxM9JuE/MYby4DJdOjQIWVFgPeWv0eiLZFPGplh5IZ3OZ1U\nRORl4Bsgv+P4RkR6eyoww3VbtmxhwIABAEyaNMl0zrsgNBRatQJVUoYXjx07lthY394y6MyZM6xY\n8QbQKyXuZE0mPime71p/R9agzLdfjmExVXXqALYB4anehwPbnL3fm4f918pc4uLitFKlSgpojx49\nrA7H79WpU0cBHTVqlNWh3NR7772ngDZv3tzqUAw/5/jezPD3ryvzVASwpXpvc5wzfMDrr7/Ozp07\nKVOmDJ9++qnV4fi1hIR/aivDhw8nOTnZ4ojSdvnyZYYN2wYE89ZbZk8Uwze4UjeeDPyfiMx1vH8C\nmOj+kAxXTZ48mTFjxhASEsLMmTMJDw+3OiS/NW0abNwIw4a1pESJEvz111/Mnz+fFi1aWB3af4wY\nMYrY2Kd54IFY6tWrZ3U4hmHnSrUGqAG87Diqp6dqBDQB9gD7gLfT+Dw3MAfYCqwDKjp7b6rr3Fkr\n9GkrV67U0NBQBXTChAlWh+P34uJUExPtr4cPH66A1q9f39qg0hAdHa158+ZVQBcvXqxJtiSrQzL8\nHG5q/vJ2X0cQ8BdQHAgGtgDlr7vmf8AAx+tywBJn7031DDf+o/ZdO3bs0Ny5cyugL730ktXhBJzL\nly9rrly5FNA//vjD6nD+5dNPP1VA69Spo7/u+1Ufnvaw1SEZfs5dScXba3/VAvar6hFVTQRmAde3\nK1QElmH/DfcCJUQkn5P3Zhq7du2iUaNGXLx4kSeeeIKRI0daHVJAWboUYmJypix1079//2t/sFju\nypUrDBoUB9Sny1td6PRTJ96r/57VYRkG4P0FJYsAR1O9P+Y4l9pWoBWAiNQCigFFnbw3U9i4cSMP\nPvggJ0+epEGDBsycOdPMmnezrVvh+HF47bXXyJMnD3/88QdLly61OiwARowYQUzMPCpXC+XLM1/S\nv15/6hU3fSqGb3C6o15EQoEngRKp71PVwW6OaSgwQkQ2AduBzfx71JlTBg4cmPI6IiKCiIgIN4Vn\nrRkzZtC1a1fi4+Np0qQJc+bMIVu2bFaHFXAcg7+AXPTp04d+/frRv39/GjZsaOluiefOnWPo0KHA\nZfJ3fYjCBarQq1avW95nGNeLjIwkMjLS7c91ZZmWRcAl7Gt+pXzJq6rTa0CISB1goKo2cbzva3+E\n3nDar4gcAqoAlZ29NxCXafn777/p3bs33333HQBdunThq6++SlnnyfCc2NhYSpUqxZkzZ/j5559p\n3ry5ZbH06vUOo0aN4J42FbA1tLG682qyB2e3LB4jcLhrmRZXOtl3ZLQDB8jCP53tIdg72ytcd00u\nINjxuiswxdl7Uz0j/b1VGZCcnKznz5/XI0eO6M6dO3XPnj166tQpjY+PT/czT506pQMGDNCcOXMq\noNmzZ9dx48ZpcnKyGyM30pKcrNq4seqBA/+MBKtQoYJevXrVkngOHz6sWbJ0VhiqmzZt0svxly2J\nwwhMuKmj3pWaynjgS1XdnpEkJiJNgBHY+3MmqupQEenu+IXGO2ozU4FkYCfwgqpeutG9NyhDnf29\nMmrv3r3MmTOHxYsXs3XrVs6fP5/mdQULFqRUqVKULFky5WehQoUoWLAgBQoUICwsjMTERGJjYzl8\n+DAbN25k6dKlLF68GJvNXjF89NFH+eqrryhVqpRXfjcD9u6FsmXh6tUEKlWqxIEDBxg5cmTK3ive\n9Oyzz/Ltt9/Stm07Zs6c4fXyjcDmrpqKK0llF3AXcAhIwD6bXlXV5xaW9EZSWb58OUOGDGHZsmX/\nOp8jRw7y5MlDeHg4NpuNixcvcv78+ZTE4KqsWbPStGlT+vTpwwMPPOCO0I10+umnn3jiiSfIkycP\n+/fvJ2/evF4re9myZTRs2JCwsDB2795NiRIlvFa2kTm4K6m4MqP+0YwWFghOnDjBSy+9xE8//QRA\n9uzZadOmDc2aNeO+++6jcOHC/+nItdlsHD9+nEOHDnHw4EEOHjzI4cOHOXXqFKdPn+bUqVMkJiaS\nNWtWwsLCKF68OBUqVKB+/fo0adIkZeVZwxqqMH8+NG36OA0bNmTp0qUMHDiQL7/80ivlX716lY4d\nvwE6079/KZNQDJ/mdE0FQESqAtfGLq5U1a0eiSqDPFVTWbhwIR06dODs2bOEh4fTt29fevfuTa5c\nudxeluE7VKF3b+jXD86f3061atUQEdavX++VveA//vhj3pn8NQWz3M3hLd+ZwRmGR1jR/PUK9o7z\nOY5TLYHxquqdP9dc4ImkMmbMGHr16kVycjKNGjViypQpFC5c2K1lGP7h1VdfZcSIEdx9991s2LCB\nkJAQtz3bZrPx25w5rJ4yhaxXrnAJmLxrFZc7JDGp3iSef/x5t5VlGKlZMfor0y59//HHHyuggL7/\n/vtqs9nc+nzDf1y5ohodHaOlSpVSQAcOHOi2Z58+fVpfql1bl4WFabK9gqTJoAuC0Ih8OfX06dNu\nK8swroe31/7CPhExLNX7MGC7O4Jw9+HOpDJ27FgFVET066+/dttzDf/Upo3q0qWqy5cvV0CzZs2q\nW7ZsyfBzbTabvlS7tsY4ksn1RwzoS7Vrmz9oDI9xV1JxpfnrdaAjkHrp+ymq+kVGa0vu5q7mr19/\n/ZVmzZqhqowdO5bu3bu7ITrDn0VHQ86c9tc9e/Zk9OjRVKhQgfXr15MjR450P3fh7NmEtW9Pg/j4\nG16zLCyMq998Q5NWrdJdjmHciBV71H8OdAbOO47nfTGhuMuBAwdo164dqsrAgQNNQjGAfxIKwNCh\nn1CxYkV2795N165dycgfMqsmTybiJgkFoEF8PCsnTUp3GYbhDS4tKKmqG1V1pOPY7KmgrHb16lVa\nt27NxYsXadGiRcq+74Zxzfr10LNnDn788Udy5MjBrFmz+Oijj9L9vKxXrtxyG1VxXGcYvuyWSUVE\nVjl+RovI5VRHtIhc9nyI3jd48GC2bNlCqVKlmDp1KkFB3l7M2fB1NWrAe+9B+fLlmT59OiJC//79\nmTgxfZuhJmXPzq3qOeq4zjB82S2/LVW1ruNnTlW9LdWRU1Vv83yI3rVhwwaGDh2KiDB16lQzB8VI\nU9ascNdd9tfNmj3Bl1+OAqBr166MHj3a5edVa9OGRXLz/x2Xh4VRr3Nnl59tGN7k9J/gIpLWasA3\nXF3YH9lsNrp164bNZuO1116jbt26Vodk+IERIyA2tgdDhw5FVenZsydvvvkmV69eder+bdu2MWDo\nUPprCWJvcE0s8GPVqjzyxBNui9swPMLZYWLApjTOBdQ8ldGjRyugxYoV09jY2HQ9w8h84uJUo6Pt\nr7/++msNCgpSQGvUqKErV6684X3R0dE6aNAgDQ0NVcqjWbsHa+eqVXXpdfNUloaF6Uu1a5t5KoZH\n4a0hxSLSA3gJKAUcSPVRTmC1qj7n9kyXQekZUnzu3DnKlCnDhQsXmD17Nk8++aSHojMC2YEDMHz4\nYRYsaMDhw4cBqFu3Li1atKBKlSrkyJGDY8eOsXz5cn744UfOnx8Ht/ckrOdFFnVYRL1S9fht7lxW\nTZ5M1itXSMqenXqdO/PIE0+Yvj3Do7y2TIuI5ALyAB8DfVN9FK2qaa/zbrH0JJU33niDzz//nIYN\nG7J48WJLd/cz/NfRo7B5M0REXGbYsGEMG/Y9sbF/A+ccV5QGYoDTAFS4+3linl3G+40G8EKNFyyK\n2jAsWPvLn7iaVI4dO8Zdd91FQkICmzZtonr16h6MzshMPv44jqiobSQkjCMqKop9+5pRqFAiTZrE\n0LJlSz7Y+wH5wvMxttlYq0M1MjmvLX0vIqtUta6IREPKqMdrBasGwAiwwYMHk5CQwFNPPWUSiuFW\n/fplA2o7jn/b9fcuLm27xMwnZ3o9LsPwlExfU4mKiqJ06dIkJyeza9cuypUr5+HoDOMfqmqaWg2f\n4PVlWkSkjYjkdLzuLyJzRMTv/6z/4osvSEpK4umnnzYJxfA6k1CMQOPKgpLbVPVuEakLDAE+Bd5T\n1f/W6y3mbE3lwoUL3HnnncTGxpq+FMMwMjWv11SAa5usP4Z9c64FgPt2J7LAmDFjiI2NpVGjRiah\nGIZhuIErSeW4iIwDngZ+FZFQF+8HQESaiMgeEdknIm+n8fltIvKziGwRke0i0inVZ6+JyA4R2SYi\n34hIupNafHw8I0eOBKBPnz7pfYxhOG30htHM3T331hcahh9zJSk8BfwGNFbVi8DtwFuuFCYiQcBX\nQGOgEtBWRMpfd1lPYKeqVgMaAMNEJKuIFAZ6AzVU9W7sI9eecaX81GbMmMHp06epVq0aDRs2TO9j\nDMMpKw6vYNCKQVQrWM3qUAzDo1zZT+UK9hn1jUWkF5BfVX93sbxawH5VPaKqicAsoMX1RWGfrY/j\n5zlVTXK8zwKEi0hWIDtwwsXyU4wZMwawT3o0naWGJx27fIy2P7Zl2hPTKJmnpNXhGIZHuTL66xXg\nGyC/45ghIr1dLK8IcDTV+2OOc6l9BVQUkRPAVuAVAFU9AQwDooDjwEVVXeJi+QD8+eefbNq0ibx5\n89K6dev0PMIwnJKQlEDr71vTu1ZvGt/V2OpwDMPjXGn+egGorarvqep7QB2gqwdiagxsVtXCQHVg\nlIjkEJHc2Gs1xYHCQA4ReTY9BYwbNw6Ajh07EhYW5p6oDSMNfRb3oXDOwvSt2/fWFxtGALjljPpU\nhH9GgOF47Wq70XGgWKr3RR3nUnse+zpjqOoBETkElAdKAAevrTcmInOA+4E0pyMPHDgw5XVERAQR\nEREAXL58mW+//Raw731hGJ7Uo2YPCucsbJpYDZ8TGRlJZGSk25/ryjyV14GOwLXhK08AU9SFfepF\nJAuwF2gInATWA21VdXeqa0YBZ1R1kIgUAP4EqgJ3AROBmkACMBnYoKqj0ijnhvNUxowZw0svvcSD\nDz7okX+ghmEY/shra39do6qfi0gkcG3nqufVxX3qVdXm6OT/HXvT20RV3S0i3e0f63jsEyuniMg2\nx219HLWT9SIyG9gMJDp+jnexfMaPt9/SvXt3V241DMMwnJCp1v7aunUr1apV4/bbb+fEiROEhoZa\nEPGgVYoAAA2pSURBVJ1hGIbvsWLtrzARed2x5tePjomIftXLPWPGDACeeeYZk1AMj1hzdA2B+Iea\nYTjLldFf07BPWPwSx7BfYLongvIEm83GzJn2Pv3nnvO5zSqNAPDt9m9pP7c9sYk32mneMAKfK6O/\nKqtqxVTvl4vILncH5CmRkZGcOHGCUqVKUadOHavDMQLMttPbeHnRyyxpv4QcITmsDscwLONKTWWT\niKR8G4tIbewjs/zCtaav5557zgzvNNzqfNx5Wn7XkpFNRlK1YFWrwzEMS7kypHg3UA77jHawzzfZ\nCyRhH7l1t0ciTIfrO+qvXLlCwYIFiY6OZu/evZQtW9bC6IxAYku20ezbZlS4owKfN/7c6nAMI928\nPqQYaJLRwqwyf/58oqOjqVWrlkkohludiztH8VzF+eThT6wOxTB8givzVI54MhBPmjVrFgDt2rWz\nOBIj0OQPz8/YZmOtDsMwfEbAz1OJjo4mX758XL16lWPHjlG4cGGLozMMw/A9Vuz8iIiU9Le5KQsW\nLCAhIYEHHnjAJBTDMAwPc3Xnxjexr06MiNRz7Ff//+2de4xU1R3HP19BW5WHgGBVFDGNJRREEBQq\nitrUR62AiC8sEpM2jfWRVNtUrdb6aNFGTaOJbRFqGx+1FMujoBWxi4+KAXmtKOALlRrtanyAuCqw\nv/5xzshlnNmd2Z3dmXv5fZLJnjlzzj3f3/nNzu+ec+89p6aZNWsWgC9x71QEM/OHGx2nGcqa/pI0\nhbAy8RNmtkHSeDOb027qWklu+mvLli307t2bxsZGNm7cSN++fastzUk5U5+ayrambVw75tpqS3Gc\nilKV6S/gIOBz4HJJ/waGt1VAe/Lwww/T2NjIqFGjPKA4bWbhqwu5c+mdXDj0wmpLcZyapZxbigFe\nA2aZ2QOSegET2kFTxfCpL6dSbPhgA5NnT2bmxJn07eYnKI5TjHJHKn8DBsX0ocDXKiuncjQ2NrJg\nwQLAg4rTNj7Z+gkTZk7g6tFXM+aQMdWW4zg1TVkjFTPbDqyI6WXAsvYQVQkWLVrEli1bGD58OAcf\nfHDLFRynCLc8fQsDew/ksqMvq7YUx6l5yp3+Sg1z584FYNy4cVVW4qSdK0dfiWG+ZpzjlEBmH37s\n06cPDQ0N1NfXM3jw4GpLchzHqWkqdfdXZoMKQP/+/Xn11Vf9DNNxHKcFqnVLcaoYN26cBxTHcZwO\nJNNBZezYsdWW4KSQ6Sums/mzzdWW4TipJLNBpUePHhx77LHVluGkjOkrpnPbktswsjct7DgdQWaD\nymmnnUbnztm7uW3x4sXVltCuVNO+pW8t5arHr2L2ObPp9pVu7dKG+y/dZN2+SpDZoJLVW4mz/qWu\nln0NWxqYOHMid59+NwP2HdBu7bj/0k3W7asEmQ0qJ598crUlOClhW9M2zpl1DhcMuYDxA8ZXW47j\npJrszQ9FunbtWm0JTkrY3rSdMwacwcUjLq62FMdJPZl+TsVxHMcpHX/40XEcx6kpMntNxXEcx+l4\nPKg4juM4FSNVQUXSDEn/k1TfQrkRkrZKmpDIO0XSOkkvSfp5+6stjzba9rqk1ZJWSlra/mrLpyX7\nJI2R9KGkFfF1TeKzmvYdtNm+1Psvljk+2rBGUl0iP/X+i2WK2VfT/ivhu/nTqH2FpOclbZO0T/ys\nfN+ZWWpewGjgCKC+mTK7AY8D84EJibxXgH7A7sAqYEC17amEbTH/NaBHtW1oi33AGGBeEZtr2ndt\nsS9D/usOvAAcGN/vmzH/FbQvDf4r5bclUfZ7wKK2+C5VIxUzexr4oIVilwKzgIZE3lHAy2b2hplt\nBR4EaurpyDbYBiBqfNRZon2F7jyped9Bm+zL5afdf5OAh8zsrVj+vZifFf8Vsw9q3H8lfjdznAf8\nNaZb5bua7YjWIOkAYLyZ/Z6d/4EPBDYm3v835qWGZmwDMOBRScsk/bDj1VWMkXEYvkDSwJiXet8l\nKGQfZMN/hwE9JdVFOybH/Kz4r5h9kA3/IWlP4BTgoZjVKt9l7eHH3wE1OWdbAfJtSwaWY8zsbUm9\ngcckrY1nJ2liOdDPzD6RdCowh/CPnBWasy8L/usMDANOBPYGlkhaUl1JFaWgfWb2CtnwH8DpwNNm\n9mFbDpK1oDIceFBhE5V9gVMlbQPeApIb1feNeWmikG1bzWyemb0NYGbvSppNGLam6kttZh8n0o9I\nuktST7Lhu6L2mdn7WfAf4Sz2PTP7FPhU0pPAEDLiP4rb90pG/AdwLjumvqCVvkvj9JcoMjdtZofG\nV3/CtYcfm9k8YBnwdUn9JO1B6Lx5Haa4dMq2TdJekroASNobOAlY02GKy6OofZL2S6SPIjyY+z7p\n8R20wr6s+A+YC4yW1EnSXsDRwFoy4j+K2Jci/zVnG5K6E24mmZvIbpXvUjVSkfQAcDzQS9KbwHXA\nHoCZ2bS84l8sFWBm2yVdAiwkBNIZZra2Y1SXRmttA/YDZissTdMZuN/MFnaA5LIowb6Jki4CtgKN\nwDmQDt9B6+0jI/4zs3WSHgXqge3ANDN7MdZNvf+K2SepPzXuvxJ/W8YDj5pZY65ea//3fJkWx3Ec\np2KkcfrLcRzHqVE8qDiO4zgVw4OK4ziOUzE8qDiO4zgVw4OK4ziOUzE8qDiO4zgVw4PKLoykFp/6\nlTRN0oCYvqoV9Te3XmF5lKKnlcf9og86AknXSzqxjcfYQ9JjcTnzs9p4rCFxaZl2Ia6nNazCx2xX\nzU5x/DkVp2QkbTazrmXW2WRm3dpLUy0iaTcza6qyhpHADWZ2UgWONQUYbmaXllGnk5ltL7FsHXCF\nma1orcYCxyxbs1MZfKSyC5MbRShsIFUn6e+S1kq6N1GmTtIwSVOBPeOZ77159feWtEjScwqbFY1t\nod29JM1XWLG3PncmHdtZrLDa6yO5pU2ihttj/guShkt6SNJ6STfm21OgvXviWltLJL0S7Z0h6UVJ\nf0qUu0vSUoWNiq7L74OYPi9qrpd0c7JtSbdKWgmMzGv/B/G4K2MffzXmz1Fc7VbSjxL9eo/iJmyS\nblbYFGqVpN8WsK2HpNmx35+RNEhhYcN7gRHRX/3z6gyJfbEq9mP3Anb2krRBUmfgBuDs3Kgn+m+G\npGclLZd0eqwzRdJcSY8Di/La7Be/W/fFfp+Z64e8csV8sEHSr2J7qyUdFvO/pEXS7vmaC30vnHai\nLZu/+CvdL2BT/DuGsN/C/oT1gZ4BvhU/qwOGJcsXqN8J6BLTvQh7MFCoTsybAPwx8b4rYYmL/wC9\nYt7ZhGUhchqmxvRlhEXt+hCWmthI3CCpUFsx/x7ggZgeC3wEDIzvnwMOj+l94t/dYpuDkn0Q++cN\noCc7NkwbG8s0AWcWab9HIn0jcHFM9wFeImyitA7ontA7IbazLlG3W4Fj3wFcG9MnACsTPi22Kdhq\nYHRMXw/cXsDXvYDXYnoKcEei/q+BSTHdHVgP7BnLvZmzI6/NfrGPRsb3M4DLC7RbzAcbCOvdAVxE\nWCalJS13FLLfX+378pGKk2Opmb1t4b9zFXBIGXUFTJW0mnCGeoCkPs2Ufx74jqSpkkab2WbgG8Ag\nwtLhK4FfAAck6sxL1F1jZg1m9jlh172DStD4z0T9dyyuS0XYze+QmD5X0nJgJTAwvpKMAOosrCzc\nBNwPHBc/2w78o0jbgyU9qbCd6yTgmwBm1kBYh6mO8AP7UV69j4BGSdMlnUFYMyyf0YRRCWZWR9jz\no0sRHUjqRvjRz11/+kvChlI5Cbgy+mkxIbjnVrN9rIAdOd40s2dj+r6oPZ/mfDA7/l3ODp81p8Wp\nAqlaUNJpVz5LpLdT+LtRbJXT8wnL8Q81syZJG4AvTW3kMLOX4zTLd4Eb43TJHEKwOKYFfU15Wpvy\ntUq6CTgtNGW5C8DN1pd0CHAFcKSZbZJ0TxEbivVBYwzIhfgzYUSzRmGuf0zis8OB9yiw+ZGFBf2O\nAr4NnAVcEtM7FStRXylsY8eUeFH/Rc40s5d3ajhcx9lSRns7aS/BBzm/Jb+fakaLUwV8pLJrU+4P\n0Odxjj2/fnegIQaUEwhTHUXbkLQ/4Uf4AeBWwtTSeqB37sdAUmftvDtiKQjAzK4xs6GJgFKwXB7d\ngI+BzQrXcgrdObQUOE5ST0mdCFuvLm7mmDm6AO/Euf7zvxARAsbJwFDgZ5KS/YbCEuv7mNm/gMsJ\nASifp4Dvx/LHA+9aYu+WfMxsE/CBpFzwngw8EdOvE/btgRDEcmwm9E+ORwnTkDmdRxRrL4+DJR0d\n05Oi9iSl+CCfYlryNTsdhAeVXZtiZ9ZWJD0NqNeOC/m5z+4nXBReTfiBW1ukfo7BwNI4ZfFL4CYL\ne2BPBG6RtIow/TGqBZ3NaS1WpmAdM6snTPutJUzN5N+ebGb2DnAlIZCsBJ4zs/klaLyWEJCeisdH\nYX+KacCF8bhXALmbBnLH6gbMj/36JPCTAse+HjgylvkN4VpCS0wBbo39PIRwURtCgL8oTj/1TJSv\nAwYmLnrfCOweb1ZYk6jfEuuBiyW9COwD/CHml+KDYv2b1PJ8Qku+ZqeD8FuKHacF4rWQ083sjWpr\nSStxFDbfzAZXW4vTvvhIxXGaQdJCYLUHlIrgZ7C7AD5ScRzHcSqGj1Qcx3GciuFBxXEcx6kYHlQc\nx3GciuFBxXEcx6kYHlQcx3GciuFBxXEcx6kY/wdhLBqfSbPQCAAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig = plt.figure(figsize=(6,4))\n", "ax = plt.subplot(111)\n", "ax.set_xlim(a_grid[0],a_grid[-1])\n", "ax.set_ylim(np.min(x_exact),np.max(x_exact)*1.01)\n", "ax.set_xlabel(\"initial semi-major axis of outer planet\")\n", "ax.set_ylabel(\"$x$ position of inner planet after 10 orbits\")\n", "ax.plot(a_grid, x_exact, \"-\", color=\"black\", lw=2)\n", "ax.plot(a_grid, x_1st_order, \"--\", color=\"green\")\n", "ax.plot(a_grid, x_2nd_order, \":\", color=\"blue\")\n", "ax.plot(a_0, x, \"ro\",ms=10);\n", "plt.savefig('paper_test1.pdf',bbox_inches='tight'); # Save to file. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following code produces an interactive version of this graph where one can change the initial semi-major axis $a_0$ and immediately see the new plot. It uses the `ipywidgets` tool `interact`. Move the slider and see how REBOUND accurately calculates the first and second derivate using variational equations." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEKCAYAAADaa8itAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XdYleUbwPHvjYgobsMd7pk5yllqlFlamWXZUNNMrdTK\nykpbpi21ZZlmmebMrBxlOXJirsocOHFPTP25AZmH+/fHORIpyDlwznnPgedzXe8l513PjRk3zxZV\nxTAMwzDcIcDqAAzDMIzcwyQVwzAMw21MUjEMwzDcxiQVwzAMw21MUjEMwzDcxiQVwzAMw20CrQ7A\nE0TEjJM2DMNwkapKTt+Ra2sqqporj7feesvyGMz3Z74/8/3lvsNdcm1SMQwj70pJTbE6hDzLJBXD\nMHKVk3EnqTO2DnFJcVaHkieZpOJnwsPDrQ7Bo8z359984fsbGjGUe2rcQ0hQiNvf7Qvfn68Td7al\n+QoR0dz4fRmGcXU7/7eT1pNbs+uZXZQsWNLqcPyKiKCmo94wDONfryx9hVdbvmoSioVy5ZBiwzDy\nnuUHlrPjfzuY1XmW1aHkaab5yzCMXGHLiS2cvniaW6vcanUofsldzV8mqRiGYRimT8UwDMPwPV5N\nKiIyUUROiMiWq9wzWkT2iMhmEWmY7nw7EYkSkd0iMsg7ERuGYRiu8HZNZRJwZ2YXRaQ9UE1VawBP\nAV86zgcAYxzPXgc8KiK1PR+uYRiG4QqnR3+JSGdgkarGiMgbwA3Au6q60dl3qOpqEal0lVs6AlMd\n9/4pIsVEpAxQBdijqoccscx03BvlbNmGYeQ+0yKn0b5Ge64pdI3Lz54/f57t27dz+PBhjhw5QkxM\nDAkJCdhsNkJCQihcuDBlypShUqVKVKlShbCwMERy3OWQ67kypPhNVf1RRFoCtwMfAuOAZm6MpwJw\nJN3no45zGZ1v6sZyDcPwM9tObmPg4oHcU/Mep+6Pi4vjt99+49dff2XdunVERbn2O2nx4sVp2LAh\njRo1omXLltx6662UKFEiO6Hnaq4kFZvjz7uB8ao6X0Te9UBM6ZlfCwzDyNDLS17m9VavU6Jg5j/Y\nVZVVq1bxxRdf8PPPP5OQkAC0BnoTFPQa9erVo2DBTpw9ezsPP/wbwcHBxMcXIi4uEdXjHDt2jMOH\nD7N7925OnjxJREQEERERjBo1ChHhhhtuoF27dnTq1IlGjRqZmgyuJZVoEfkKaAuMFJECuL9PJhq4\nNt3nio5zQUBYBuczNXTo0LSvw8PDzZo9hpGLLN63mL1n9tL3kb4ZXldV5s6dy9ChQ9m6NQH4FPie\n5s2bc+ed91Gv3u3ce++zBAUFkZgIsbFQqpS90WXuXPjzTxgx4r/v/Oeff9i0aRN///03K1asYO3a\ntWzYsIENGzbw3nvvUaVKFTp16kTnzp1p2rSpzyeYSwnS7VxYa78Q0Amo4fhcDrgjG2v2Vwa2ZnLt\nLmC+4+vmwB+Or/MBe4FK2BPMZqDOVcpQwzBypxRbil7/xfU6e8fsDK///vvvWrdudwVRQEuXLqf9\n+4/Sw4cPZ6u8QYNUly278nxcXJwuWrRI+/btq2XLllUg7ahVq5YOHz5cjx49mq0yreD4uZnzfVmc\nvhFGOnMui3fMAI4BicBhoCf2UV5PprtnjCOBRAI3pDvfDtgF7AEGZ1GOO/+uDcPwIWsOr9HwyeGa\nmpr6n/Pnzp3Tp556yvGD/Xu95pqmOnbsWE1MTMxRebt2qZ4+ffV7UlJSdNWqVfr888//J8EEBARo\nu3btdObMmZqQkJCjODzNiqSyMYNzW9wRhLsPk1QMI3dLtiX/5/OKFRFapkwTBTR//vw6ZMgQjY2N\ndXu5UVGq3bplEVtysv7666/64IMPav78+dMSTKlSpfSFF17Q7du3uz0ud3BXUslymRYR6Qv0A6oC\n+9JdKgKsUdVuV32BBcwyLYaRN6SmpvLhhx/y6qvfozqEZs1GMHHiRK677jqPlJecDFu2wI03Onf/\n6dOnmTFjBhMnTiQyMjLt/E033UTv3r156KGHCAlx/74v2eG1tb9EpBhQAhgODE53KUZVz+Q0AE8w\nScUwcr/4+Hi6devGnDlzABg8+DXeeWcYgYHeWXxdFX76CTp2hIAshiypKhs2bGDChAnMmDGDmJgY\nAIoUKUKXLl3o3bs3N954o6Wd+2ZByaswScUwcrczZ85w662D2LKlMMWKTWLatGl06NDBqzHExcEL\nL8Cnn0KhQs4/Fxsby48//siECRNYu3Zt2vkGDRrQp08funbtSvHixT0Q8dV5s6ayWlVbikgM9rbB\n9IWqqhbNaRDuZpKKYeQuiSmJFAgsAEB0dDRt27Zl584zXHPNLUREDPFYc5enbd++nYkTJzJ16lRO\nnz4NQHBwMJ07d6Z37960atXKa7UXU1O5CpNUDCN3uXP6nQxoNoBGhRvRuvUt7N27h3r16rFw4UIq\nVqxodXicPw8DBsCXX0JwsOvPJyYm8tNPPzFhwgSWLl2adr5mzZr07t2b7t27U6ZMGTdGfCWvJxUR\nCcbeYd8Se41lFfClqibkNAh3M0nFMHKPRXsXMWDRAJZ1Wkarm5/l4MHuNGz4NsuWLaNkSd/YNlgV\nFi2Cdu0gpxWL/fv388033zBp0iSOHTsGQGBgIB07dqR79+7ccccdBGcnc2XBiqTyAxADTHec6gIU\nV9XOOQ3C3UxSMYzcISU1hYZfNuT1Fq8zsudIIiMjqVnzTtasmc4117i+iKQ/SUlJYeHChUyYMIH5\n8+djs9lXygoJCeGuu+6iU6dO3HXXXRQt6p4eCCuSyg5VrZvVOV9gkoph5A5fb/ia6VumU2BmIZYs\nXkTNmjVZtWoVpUuXtjq0TK1bB999B6NHu++dx44dY8qUKcyaNYuNG/9dGD5fvnw0b96c22+/nbZt\n29K0aVPy58+frTKsSCrTgTGq+ofjczOgv6p2z2kQ7maSimH4v9ikWGp+XpMme1sw78tBlCz5OH/9\n9QvVqlWzOrSrungR9u2D66/3zPsPHjzITz/9xJw5c1i7dm1aDQbstZgbb7yRpk2b0rRpU5o0aUJY\nWBgBWY15xrujv7Zi70PJD9TCvryKYl+HK8rUVAzD8ISYxBj6j+7PtFemERRUhYiIb2nRooXVYfmU\n8+fPs3LlSpYsWcLSpUszXM6/YMGC1KhRg5o1a1KzZk0qVKhAmTJlKF26NGXKlKFYsWIUKlSIokWL\nei2pVMK+GvG1wKHLr6tj4yxfYpKKYfi/NWvWEB4eTkpKCjNnzuThhx+2OiSXvfQSPPEE1PXSr97/\n+9//WL9+PX/99Rd//fUXGzZs4OTJk04/7+3mr62q6qEKnXuZpGIY/u3EiRPUqTOUs2eL8+KLSXz8\n8cdWh5QtK1fal3QpXNi6GM6dO8eePXvYvXs3e/bs4fjx45w4cSLtiImJ4eLFi8TFxXk9qUzB3qey\nPqeFeppJKobhv2w2G23btmXFig00btyetWunZbvz2XCeu/pUXNlkqxmwTkT2icgWEdkqIltyGoBh\nGEZ6n3zyCStWrKB06WDmzRuVKxLKrl32dcLyAldWXrvTY1EYhmEA3yz7hkEjjwCVmDTpC8qVK2d1\nSG6RnGyfdZ8XuLRMi4g0AFo5Pq5S1cir3W8V0/xlGP4nJi6G0CGhJH7blZ53leSbbz60OqQ8xevN\nXyIyAPgWKO04povIszkNwDAMA+CB9x8g8XQiNYuvZsyYYVaH4zF791odgWe50qfSC2imqkNUdQj2\nPeT7eCYswzDykt///J0lR4/DYpg6ZSqFXFlL3o8cP24fYpxuvmKu40qfigDp/yps/HcZfMMwDJel\npKTQ+ZMesHwi/R/+jWbNmlkdkseULWsfZmzhXlwe50pSmQT8KSJzHZ/vAya6PyTDMPKSj0d9zMmS\nB6kQ2p0RI66cEZ7b5OaEAq531N+Afel7sHfUb/JIVDlkOuoNwz8cOXKE2rVrczH+IgvmL6B9+/ZW\nh+QViYlw110wZw4UK2Z1NHZ+u0mXiLQDPsXenzNRVUdedr048A1QDYgHnlDVHY5rL2Dv20kFtgI9\nVTUpgzJMUjEMH2Sz2fhtzhzWTJ5M4MWL/LFzP6tO9OSuB3Ywa9YPVofnVZGRUL++79Rc/DKpiEgA\nsBtoAxwD1gOPqGpUuns+AGJU9R0RqQWMVdXbRaQ8sBqorapJIvI9MF9Vp2ZQjkkqhuFjTp48ybB7\n7+XByEjCExIQ7CvTLkSY07AB7//2m08vaZ/bWTGj3h2aAntU9ZCqJgMzgY6X3VMXWA6gqruAyiIS\n6riWDwgRkUCgEPbEZBiGj0tNTWXYvffywZ9/cqsjoYB9pM9dKJ9t3sywe+8lNTXVyjC9LjUV5s2z\n7xyZW3g7qVQAjqT7fNRxLr1IoBOAiDQFwoCKqnoM+Bj70vvRwDlVXYphGD7vtzlzeDAykpBMrocA\nD0RGsjivrGXioGpPKmfPWh2J+zg1+ktEamOvUVxKANHAPFXd6YGYRgCfichG7P0mmwCbo6+lI/Z9\nXM4Ds0Ski6rOyOglQ4cOTfs6PDyc8PBwD4RqGIYzVk+axLsJCVe959aEBN745hvaderkpaisly8f\nTJhgTdkRERFERES4/b1ZJhURGQQ8ir2p6i/H6YrAdyIyU1VHuFBeNPaaxyUVHefSqGoM8ES68vcD\n+4F2wH5VPeM4Pwe4CcgyqRiGYa3AixeznNQmjvvyKlXvdtpf/sv2sGHuWcXAmZpKL+A6Rx9IGhH5\nBNiOvWbhrPVAdcfGX/8Aj2BPWOnfWwy4qKrJItIH+F1VY0XkMNBcRIKBROyd/T6/DL9hGJBSqBDK\n1WdLq+O+vCg+Hlq3hogICMmsjdBPONOnkgqUz+B8Occ1p6mqDXgGWIw9Ic1U1Z0i8pSIPOm4rQ6w\nTUR2Yl8ZeYDj2b+AWdibwyKx//sc70r5hmFYo2XPnkQEB1/1nhXBwbR64omr3pNbFSwIM2f6f0IB\n57YTbgeMAfbwbyd7GFAdeFZVF3o0wmwwQ4oNw7ekpqbS54YbGJ1JZ30c8EqzZny+di0BAd4eP2SA\n+4YUZ9n8paqLRKQm9uHA6Tvq1ztqHoZhGFcVEBDAhRo1aL0rkveSArgzNTVtnsqK4GBmN2jAW/Pm\n5fmEkpAA69bBrbdaHUn2OTX6S1VTgT8uPy8iPVV1ktujMgwjV9m6dSuzZp0hX75CHP/8Q95YsIDA\nixdJKVSIVk88wef33ZfnEwpATAx88w2Eh/vOTHtX5WhGvYgcVtWwrO/0LtP8ZRi+pVOnB5k79356\n9drChAkjs37A8DqvLdNylX3oBaipqgVyGoS7maRiGL5j+/bt1KtXjwIFCnDgwIFcs0VwbuO1PhWg\nDPZRWJfP+RRgbU4DMAwjd3v//fcB6N27t0koTvr9d1iyBN55x+pIXOdMUvkVKKyqmy+/ICIRbo/I\nMIxcY+OOjcxYVY2AgNt45ZVXrA7Hb9StCyVLWh1F9nh96XtvMM1fhuEb6g+sz9YN8Gj5NsyYMcrq\ncIyr8Gbzl2EYhstWb1/N1sCtyCZh2NezrQ7HL124YF++xVc28nKGGcNnGIZH9Jj4DKwpzCN3P0KN\nGjWsDscvDRsGixZZHYVrTPOXYRhut2DLAu5+/RtYEM7WyHDq1atndUh+yZuLTHp1ky6xuzanhRmG\nkTe8N+c9ODqbe+5abBJKDvjjBEinkorj1/4FHo7FMIxcID4+nt1jd8NmeOmlgVaH4/dSUuCDD8Dm\nJ4tiudKnslFEmngsEsMwcoWpU6dx6tQjNGzYjNatW1sdjt8LdAyniouzNg5nOd2nIiJRQA3gIPZF\nRQV7Jaa+x6LLJtOnYhjWSE1NpU6dG9m9+36mT69O165drA7JcJIVQ4rvzGlhhmHkbgsXLmT37s1U\nrHiKhx7ab3U4hgVcaf46DLQCeqjqIeyrVpfxSFSGYfidmMQYPv74YwAGDBhA/vz5LY4od4mIgN69\nrY4ia640f43DvtPjbapaR0RKAItV1ef6WUzzl2F4174z+7jp65s4+epLBBdYx/HjkyjmTzP2/EBs\nrH0yZPmM9uF1A68OKXZopqr9gQQAVT0LBOU0AMMw/N+ry16l4vGKoJ/Ro0dtk1A8oHBhzyUUd3Kl\nTyVZRPJhb/ZCREJxcY96wzByn7VH1rLm8BpOTzgNJDJwYE+rQ8rVoqPhmmuggM9tOmLnSk1lNDAX\nKC0i7wGrgeEeicowDL+gqgxcPJDWybeQGFucO+64wyzJ4mHPPw9bt1odReacrqmo6rcisgFog304\n8X2qutNjkRmG4fN+3PEjiSmJrP7iBPAD/fufszqkXO+HH3x7pr3TNRURGamqUao6VlXHqOpOEXF5\nX1ARaSciUSKyW0QGZXC9uIjMEZFIEflDROqmu1ZMRH4UkZ0isl1EmrlavmEY7nNjuRvpdU0vjh5Z\nzrXXPsbdd99tdUi5ni8nFHCt+attBufau1KYiAQAY7DPebkOeFREal9222vAJlVtAPTA3ux2yWfA\nAlWtAzQATE3JMCxUrWQ1fpvyGwB9+z5Nvnz5LI4ob4iOhpkzrY4iY1kmFRHpKyJbgVoisiXdcQDI\nbP/6zDQF9qjqIVVNBmYCHS+7py6wHEBVdwGVRSRURIoCrVR1kuNaiqpecLF8wzDc6ODBg/zySzL5\n81enV69eVoeTZ4jAkSNWR5ExZ/pU7gLuAXYBHdKdj1HVMy6WVwFI/1dxFHuiSS8S6ASsEZGmQBhQ\nEftIs1MiMgl7LeVvYICqxrsYg2EYbvLll18CNWjbtgGlS5e2Opw8o3x5ePllq6PImDNJpRqQjD2p\nXMDeSQ+AiJTMRmLJygjgMxHZCGwFNgE2ID9wA9BfVf8WkU+BwcBbGb1k6NChaV+Hh4cTHh7u5jAN\nI29LSEhg4sSJwCneeGOt1eEYLoqIiCAiIsLt781yRr2IPAf0BaoAx0iXVLAvKFnV6cJEmgNDVbWd\n4/Ngxzsy7fB3NLNdD4QA6y6VJyItgUGq2iGDZ8yMesPwkA/XfEj3Bt1ZPHcx3bt3p2HDhmzcuBHx\n9R7kXGjAAHj8cWjUKOfv8tqMelUd7egYn6SqVVW1SrrD6YTisB6oLiKVRCQIeASYl/4Gxwiv/I6v\n+wArVTVWVU8AR0SkpuPWNsAOF8s3DCMHVh9ezZj1YyhaoChjxvwIfE2/fv1MQrFIjx5QvbrVUfyX\nS9sJO9b7qgEEXzqnqr+7VKBIO+yjuAKAiao6QkSesr9KxztqM1Ow96FsB3qp6nnHsw2ACdibwvYD\nPS9du6wMU1MxDDdL1VRaTGzBc02f48agG6lTpzHBwbfwv/99T+HCha0Oz8ghry99LyK9gQHYO803\nA82BdcBtrhSoqouAWped+yrd139cfj3dtUjA5xawNIy84IftP5CqqTx6/aMMemUQEEe3buVNQvEB\n8fFQsKDVUdi5Mk9lAPYf6IdU9VagEWCmzxpGHpCQksCry17l4zs+JiU5hcmTpwPQ2x/WYs/lFi2y\n96v4ClcWlExQ1QQRQUQKqGqUiGRYozAMI3fZcGwDN117E60rtWb27NmcOjWIChVO0rTp5TMCDG9r\n0wZuv93qKP7lSlI5KiLFgZ+AJSJyFjjkmbAMw/AlN4fdzE3X3gTAhAkTgKU899ynpoPeB/jaXmgu\nddSnPSRyC1AMWKSqSW6PKodMR71heMbhw4epXLky+fPn59ixY5QqVcrqkAxAFX7/HVq1ggBXOjXS\nsWKP+jSqujKnBRuG4X8mTZqE6vXcf/91JqH4EBEYNw5q1LB+I69s1VR8nampGIb72Ww2qlSpw5Ej\nE1mwIIX27W+1OiTDjazYTtgwjDzkbPxZ0v9ytnTpUo4c2UOVKj24885bLIzM8GUu7afizDnDMPxf\nqqbSdlpbluxfknbO3kEPvXr1IiC7DfeGRy1ZAqtWWRuDV/dTMQzDP8zYOoN8AfloW9X+v/2ZM2f4\n+eftwD306NHD2uCMq7J6QF6WHfUi0hfoB1QVkfT7pxQB1ngqMMMwrBGfHM9ry15jxgMz0oYM//DD\nDyQnB3Pdda2pWLGixREamWmb0a/+XubM6K8ZwEJgOPal5i/Jzn4qhmH4uE//+JQmFZrQMqxl2rlp\n06YBmxg8+EXrAjP8ghn9ZRhGmpNxJ6k7ti5/9P6D6iXty9/u3buXGjVqEBISwokTJwgJCbE4SuNq\n9u6F99+Hb75x7TkrFpQUoCtQVVXfFpEwoKyq/pXTIAzD8A1FCxTlh84/pCUUuFRLeY1bbgkyCcUP\nhIWBlTs7O11TEZFx2Jejv01V6ziWwV+sqj63arCpqRiGe6SmplKtWjUOHizPjz+O4MEHW1kdkuEh\nVsyob6aqN4jIJgBVPevYaMswnKKq7Nixg5UrV7Jt2zbOnDlDUFAQFStWpEWLFtx6661mGXUfs2bN\nGg4ePEjFiincf/9NVodjuCAlBQKztWZKzrhSZLKI5AMUQERCsddcDOOqbDYb06ZNY/To0WzatCnT\n+4KDg+nRowcDBw6kRo0aXozQyIy96Qu6detGvnz5LI7GcJaqfYvh+fPtzWHe5ErzV1fgYeBGYDLw\nIPCGqv7oseiyyTR/+Y5169bRt29fIiMjAShZsiTt27enSZMmhIaGkpiYyL59+1i2bBl//PEHAPnz\n52fQoEG8/vrrBAcHX+31hgfFx8dTtmwYFy78zt9/w4031rE6JMMFFy5A0aLO3++u5i9U1ekDqA30\ndxx1XHnWm4f92zKsZLPZ9L333tN8+fIpoGFhYTp58mSNj4/P9JmdO3dqjx49FHttWBs0aKD79+/3\nYtR509ydc3Xx3sVXnP/+++8V0Pr121kQleFtjp+bOf7568oyLQWAG7AveV8K6CwiQ3Kc1YxcJzEx\nkUcffZTXX38dm83GSy+9RFRUFD169LhqzaN27dpMnjyZ1atXU61aNSIjI2ncuDHr1q3zYvR5y8Xk\nizy78FlCgq4c1TV16lQAeve+y9thGW5y6hScPu3dMl1p/loEnAc2ALZL51X1Y8+Eln2m+cs6cXFx\ndOjQgRUrVlC0aFFmzpxJ+/aur+Zz7tw5unbtyoIFCwgJCeGXX37h1lvNqrju9t7v7xF5IpIfOv/w\nn/MnTpygfPkqiITwzz87CA0NtShCIydee83et9K5c9b3er35C9jmjqoR0A6IAnYDgzK4XhyYA0QC\nfwB1L7seAGwE5l2lDDdUBg1XxcfH6+23366AlitXTjdv3pyj9yUnJ+tjjz2mgBYsWFDXrVvnpkgN\nVdV/Yv7RUiNL6b4z+664NmrUKIUWWqbMHxZEZlgBbzd/AWtF5PqcJDARCQDGAHcC1wGPikjty257\nDdikqg2AHsDoy64PAHbkJA7D/VJTU+nWrRtLly6lTJkyRERE0KBBgxy9MzAwkMmTJ/P4448THx/P\nPffcw+7du90UsTE0Yig9GvSgaomqV1yzj/pax+efH/F+YIZfcyWptAQ2iMguEdkiIlsvW2DSGU2B\nPap6SFWTgZlAx8vuqQssB1DVXUBlx/BlRKQicBcwwcVyDQ97++23mT17NsWKFWPJkiXUrFnTLe8N\nCAhg/Pjx3HXXXZw+fZoOHTpw4cIFt7w7L0tMSSTqVBRvtH7jimvbtm1j48aNFC9enA4d7rEgOsOd\njhyBZcu8V54rSaU9UAO4A+gA3OP40xUVgPS/+hx1nEsvEugEICJNgTDg0rKoo4CXccyVMXzDnDlz\nGDZsGAEBAcycOZPrr89RhfYK+fPn54cffqB+/frs3r2bxx9//FIzp5FNBQILEPF4BCUKlrjimr2W\nUof27fubId25wKlTsMXVX/9zwOmkoqqHgAtAGaBSusPdRgAlRGQj9qHLmwCbiNwNnFDVzYA4DsNi\nUVFRdO/eHYCRI0fSrl07j5QTEhKSVhOaO3cuo0df3ipquIPNZmP69OlAY2rXfszqcAw3aNQIXnjB\ne+W5MvqrN/b+jIrAZqA5sE5Vb3O6MJHmwFBVbef4PBh751CmO0iKyH6gPva+lm5AClAQ+34uc1S1\newbP6FtvvZX2OTw8nPDwcGfDNJyUlJREixYt2LhxI126dGH69Olp+294yk8//cT9999PcHAwmzZt\nonbty7vkjJxYsmQJd9xxB9WqVWPPnj0e/+9pWCciIoKIiIi0z8OGDfP66K+tQDCw2fG5NvYf6q68\nIx+wF3sNJwh7cqpz2T3FgPyOr/sAkzN4zy2Y0V+We/XVVxXQKlWq6Pnz571W7uOPP66ANmnSRJOT\nk71Wbl5wabTd0KFDrQ7FcKNTp1Tfeuvq92DB6K8EVU0A+0RIVY0CarmYwGzAM8BiYDswU1V3ishT\nIvKk47Y6wDYR2Yl9lNgAV8owvGPVqlWMGDGCgIAApk2bRlFX1oPIoU8//ZSwsDDWr1/P8OHDvVau\nv7uYfJGYxJhMr8fGxjJ79mygC61b9/ReYIbHFS0KpUvb1wTzNFeav+YCPYHngduAs9hrFD433dZM\nfvSs+Ph46tWrx/79+3nttdd47733vB7D8uXLadOmDfnz52fbtm1uG22Wm7298m2iL0TzVYevMrw+\ndepUevToQfXqI1iyZBCVK3s3PsNalqz9denA3vx0LxDkjuqSuw9M85dHvfbaawpovXr1NDEx0bI4\nnnjiCQX0jjvu0NTUVMvi8AfHLhzTkiNL6oGzBzK9p02bNgro+PHjvReY4TNwU/OX2U7YcMn27dtp\n2LAhKSkprF27lhYtWlgWy//+9z9q1qzJuXPnmD17Np06dbIsFl/XZ14fShQswQdtP8jw+pEjR6hU\nqRJBQUEcP36c4sWLezlCw9NU4b77YMIEyGjVHXfVVLLsUxGRGBG54Dhi0n2OEREzCy0PSU1N5amn\nniIlJYWnn37a0oQCEBoamtb09sILLxAXF2dpPL5q64mtzNs9j9davZbpPd9++y2qSrlyC1E1CSU3\nEoHXX3dtOfzsyDKpqGoRVS3qOIqk+1xEVb3XO2tYburUqaxZs4YyZcr4TAf5U089RaNGjTh8+DAf\nfvih1eH4pJeWvMQbrd6geHDGyUJVHRMeha5dC2MqKblX06ZQoIBny3Cloz4Y6Id9uRYFVgFfqmNE\nmC8xzV90OOyQAAAgAElEQVTuFxsbS40aNTh+/DhTp07lscd8Z2LcqlWraN26NSEhIezbt48yZcpY\nHZJPWX14Nc0qNCN/vvwZXt+wYQONGzcmNDSU6Oho8ufP+D4jd1C1HwGXVSm81vyVzlTsi0B+jn1R\nyOuAaTkNwPAPI0aM4Pjx4zRt2pSuXbtaHc5/tGrVig4dOhAXF8c777xjdTg+p2VYy0wTCvy7b0qX\nLl1MQskDune3bzPsKa7UVHaoat2szvkCU1Nxr0OHDlGrVi0SExMt75zPzLZt22jQoAEBAQHs2LHD\n7HHvpOTkZMqXL8+pUwVp1Gg7GzcWsTokw8POnYNixex9LOlZUVPZ6Fhm5VIAzYC/cxqA4fsGDx5M\nYmIiXbp08cmEAlCvXj169OhBSkoKr7/+utXh+I1FixZx6tQp6tYtxuzZha0Ox/CC4sWvTCju5EpN\nZSf2GfSHHafCgF3Y1+JSVa3vkQizwdRU3CcyMpKGDRtSoEABdu/eTVhYmNUhZerIkSPUqFGDxMRE\ntmzZ4vbVknOjhx56iB9//JERI0YwaNAgq8MxvCQ+HhISoES6RaqtqKm0A6pgn/h4i+PrdmRvCXzD\nTwwdOhSAp59+2qcTCsC1117Lk0/aV/t59913LY7GOtEXonlj+ZX7pFzu7NmzzJs3Dwiic2ff6icz\nPGv4cPjhh6zvyw4z+dHI1KVRQQULFmT//v2ULVvW6pCydPToUapVq0ZycjLbtm2jbl2f6/LzuCd+\nfoIyIWUYfvvVh32PHz/eMST7ZSpV+oC5c70UoGE5Vd/oUzHymEvbB/Tv398vEgpAxYoV6dWrF6pq\nyZpkVtt8fDML9ixgcMvBWd47ZcoUAJ5/vp7Hfms1fJPlfSpi31Shoqr6xYbVpqaSc3/88QctWrQg\nJCSEAwcOEJrRug4+6vDhw1SvXh2bzcaOHTuoVculxbT9lqrSdlpbOtXpRL8m/a567969e6lRowYh\nISEcP36cwoVNJ31ec/y4fVfIevXsn71aU3H8hF6Q08IM/3GplvLcc8/5VUIBCAsL4/HHHyc1NZX3\n33/f6nC8ZuHehRy9cJQ+N/TJ8t5Lc1PuvrsnMTEmoeRFGzbAb7+5/72ujP6aAoxR1fXuD8O9TE0l\nZ1avXk2rVq0oUqQIBw4coFSpUlaH5LIDBw6kzVWJioqievXqFkfkeUNWDKFJ+SZ0qHX1cTOpqalU\nrVqVQ4cOMWTIVgID6/Hmm14K0vBZ7qqpuJJUooDqwCEgDvse8T41lPgSk1Ry5rbbbmPFihUMGTKE\nYcOGWR1OtvXs2ZPJkyfz9NNPM27cOKvD8RkrV64kPDycsLAwDhw4QMDl63UYeZIVHfV3AtWwb9DV\nATOUOFdasWIFK1asoHjx4rzwwgtWh5MjL7/8MgCTJ0/m5MmTFkfjOy510D/22GMmoeRxR47Ajz+6\n951O/4tS1UMZHe4Nx7CSqjJkyBAABg4c6Pd7atStW5e7776bhIQExo4da3U4PuHixYv86PgpcsMN\nfYiKsjggw1KpqXD0qHvf6XRSEbtuIjLE8TlMRJq6NxzDSkuXLmX16tWULFmS5557zupw3OJSbWXs\n2LFcvHjR4misN3fuXGJjY2nevDkXL1Zi3z6rIzKsVKkSuLtBwpW67xdAC+BRx+cYwPz6l0uoKm86\nemtfeeUVinp6Jx8vad26NU2aNOH06dNMmjTJ6nDcLvpCtEv3X2r66tGjB926wd13eyIqIy9zJak0\nU9X+QAKAqp4FgjwSleF1Cxcu5M8//yQ0NJT+/ftbHY7biEhabeWTTz7BZrNZHJH7bPxnI80mNCPJ\nluTU/dHR0SxdupSgoCAeeughD0dn+IsTJ+w7QrqLK0klWUTyYd+gCxEJBVJdLVBE2olIlIjsFpEr\nVrATkeIiMkdEIkXkDxGp6zhfUUSWi8h2EdkqIrmjfcYHpO9LGTx4cK6bCNepUyeqVq3K/v37mTNn\njtXhuIWqMnDxQN5s/SZB+Zz73W769OmoKvfeey8LFpRk+3YPB2n4hRIlwK07RaiqUwfQFZgHRAPv\nYV+huLOzzzveEQDsBSoB+YHNQO3L7vkAeNPxdS1gqePrskBDx9eFHeXXzqQcNZz3008/KaBly5bV\nuLg4q8PxiLFjxyqgTZs21dTUVKvDybF5UfO07ti6mmxLdur+1NRUrVOnjgL6yy+/6MyZqrt3ezhI\nw684fm46/fM8s8O1m6E20N9x1HG5MGgOLEz3eTAw6LJ7fgVuTvd5LxCawbt+AtpkUo7b/qJzO5vN\npvXr11dAP/vsM6vD8Zi4uDgtWbKkArp69Wqrw8mRpJQkrfV5LZ2/e77Tz6xfv14BDQ0N1aSkJA9G\nZ/grdyUVV0Z/FQBuAIoBpYDOl0aCuaACkH79sKOOc+lFAp0cZTbFvm9LxctiqQw0BP50sXzjMnPm\nzGHLli1UqFAhbdn43KhQoUL07dsXgI8//tjiaHJmwsYJXFvsWtpXb+/0M5MnTwaga9euZstg4wop\nKe57lysz6hcB54ENQFpvp6o6/X+oiDwA3KmqTzo+dwOaqupz6e4pAnyGPWlsxV476qOqWxzXCwMR\nwDuq+nMm5eiltasAwsPDCQ8PdzbMPMNms1G/fn127NjBF198kfZDN7f6559/qFy5MsnJyezZs4dq\n1apZHVK2nIg9QVxyHFVLVHXq/vj4eMqXL8+5c+eIjIxk9uz6PPQQXHedhwM1fFpERAQRERFpn4cN\nG+aWGfWuNF1ty2m1CHvz16J0n69o/srgmQNAYcfXgcAiYEAWz+SsHphHzJgxQwENCwvTxMREq8Px\niscff1wBfeaZZ6wOxWumT5+ugDZp0kRVVZcuVT11yuKgDJ+Dt5u/gLUiktP9WdcD1UWkkogEAY9g\n7/xPIyLFRCS/4+s+wEpVjXVc/gbYoaqf5TCOPC8lJSVtV8c333yToKC8MTr80tIz33zzDWfPnrU4\nGu+YMGECAL179wagTRvwwzVCDT/hSlJpCWwUkV0issUxrHeLK4Wpqg14BlgMbAdmqupOEXlKRC41\n6NcBtonITuzrjQ0AEJGbsY9Au01ENonIRhFp50r5xr9mzJjB7t27qVq1Kj169LA6HK+pX78+bdu2\n5eLFi3z11VdWh+Nxe/bsISIigkKFCvHII49YHY6RFzhbpcE+DPiKwx3VJXcfmOavq0pKStJq1aop\noJMnT7Y6HK9buHChAlq+fPlc3+w3ePBgBbRnz56qqtqpk+quXRYHZfgkLGj+Og5cqi30SHcYfmbq\n1Kns27ePmjVr0rVrV6vD8bo777yT6667jmPHjvH9999bHU6WVJXnFz3P8djjLj2XkpKSNuqrV69e\nALz7LlSu7OYADSMdV5LKz0BHIAX7fiqXDsOPJCYmpu2R8tZbbxEYGGhxRN4nIrz44ouAfekWdXIE\npFV+3vUzyw4sI7SQaztwLliwgOPHj1O7dm1uuukmAOrUgTzSfWZYxdkqDW4Y/eWtA9P8lanPPvtM\nAb3++uvVZrNZHY5l4uPjtXTp0grosmXLrA4nU0kpSVpjdA1dtGeRy8926NBBAf3oo49UVTUlxd3R\nGbkJfjr6y7BQbGws7733HgDvvvtunt6gKTg4mGeeeQaw11Z81Zd/f0mVElW4s/qdLj0XHR3N/Pnz\nyZ8/P4899hgA9epBtGuLGhuGy1wd/bUhJ6O/DGuNHj2akydP0qxZMzp0MJt2Pv300wQHBzN//nx2\n7txpdThXOJdwjndXvctHbT9y+dkpU6aQmppKx44dKV26NAB//gnly7s7SsP4L1eSSnugBnAHZjth\nv3P27Fk++OADAN5//31Ecj5x1t+FhoamDaf+9NNPLY7mSttPbqdHgx5cX8a1BgKbzcbXX38N/NtB\nD1C0KJj/7IanOb1Miz8REc2N31dOvPbaawwfPpw2bdqwdOlSq8PxGVFRUdSpU4fg4GAOHz5MaKhr\nneG+6Ndff6VDhw5UrVqVPXv2EBAQwIkTUKaM1ZEZvkxEUDcs05JlTUVEVjv+jBGRC+mOGBG5kNMA\nDM+Ljo5O+038Up+KYVe7dm3uueceEhISGDdunNXhuMWYMWMA6Nu3LwEBAaSmQqtWcMH832p4gamp\n5AGPP/44U6ZM4YEHHmDWrFlWh+NzVqxYwW233Ubp0qU5dOgQwcHBVoeUbXv27KFmzZoEBwcTHR1N\nyZIlAVA1TV/G1XmtpmL4tw0bNjBlyhSCgoIYOXKk1eH4pPDwcBo1asTJkyf59ttvrQ4nRy7Vtrp0\n6ZKWUMAkFMN7TFLJxVSVgQMHAvDcc8/57VLvnuZLkyEPnTuELdWW9Y0ZiIuLY9KkSQD0798/7fz6\n9e7dL8MwrsYklVxs7ty5rFy5klKlSvH6669bHY5Pe+ihh6hQoQI7duzgt99+sySGJFsSt0+7nbVH\n1mbr+RkzZnDu3DmaN2/ODTfcAIDNBoMGQXKyOyM1jMw501E/zfHnAM+HY7hLbGwszz//PABvv/02\nxYsXtzgi3xYUFMSzzz4LWLcz5Lj146hZqiatKrVy+VlVZezYsQBpkzoB8uWD5cuhYEG3hWkYV5Vl\nR72I7ABuBxYC4cB/WmdV9Yyngssu01EPL7/8Mh999BE33ngjf/75J/ny5bM6JJ939uxZrr32WuLi\n4oiMjKR+/freKzv+LLXG1GJFjxVcV9r1LRlXrlxJeHg4oaGhHDlyhAIFCnggSiM382ZH/ZfAMuzb\n+m647Pg7pwEY7rdlyxZGjRqFiPDll1+ahOKkEiVK8MQTTwDeX7rlvVXv0alOp2wlFPi3dtW/f///\nJJRZsyAhwS0hGoZTXNmjfpyq+sUm5nm5pmKz2WjVqhXr1q2jf//+aXMWDOfs37+fGjVqkC9fPg4d\nOkS5cuU8XubBcwdpPL4x2/ttp0xh12coZjaBMyUFeveGCRMgDy5GbbjI60OKVbWviDQQkWcch/fa\nBgynffTRR6xbt46yZcvy7rvvWh2O36latSr3338/ycnJXkvI1xa9luU9lmcroQCMGjUKgO7du/9n\nRYDAQJg82SQUw7tcqak8BzwJzHGcuh8Yr6qfeyi2bMurNZXNmzfTtGlTkpOTWbBgAe3bt7c6JL+0\ndu1abr75ZkqWLMnhw4cJCQmxOqRMnTx5kkqVKpGQkEBUVBS1atWyOiTDT1kx+bE30ExVh6jqEKA5\n0CenARjukZCQQLdu3UhOTqZv374moeTATTfdRPPmzTlz5gxTpkyxOpyrGjt2LAkJCXTo0OGKhPLR\nRxAfb1FgRp7lSlIRIP2sLBuXjQQzrPPiiy+yfft2atSowYcffmh1OH7v0mTIUaNGkZqaanE0Gbtw\n4QKjR48G7KP90ktJgbg48OMVZww/5UpSmQT8KSJDRWQo8Acw0SNRGS6ZNGkS48aNIygoiBkzZvh0\nc42/uP/++6lcuTJ79+7ll19+sTqcDI0ZM4Zz587RunVrWrX679yWwEB46y2zPIvhfa501H8C9ATO\nOI6equryJhQi0k5EokRkt4gMyuB6cRGZIyKRIvKHiNR19tm8aPXq1fTtax+U98UXX9C4cWOLI8od\nAgMDGTDAPt/XE8OLx28Yz6+7f83287GxsWlxvfnmm+4KyzByzh17Ejt7YE9ie4FKQH5gM1D7sns+\nAN50fF0LWOrss+ne4crWzH5r27ZtWrx4cQW0X79+VoeT61y4cEGLFSumgP7+++9ue++puFN6zQfX\n6I6TO7L9jg8//FABbd68uaampv7nWmqqat++qhcv5jRSIy/Bgj3q3aEpsEdVD6lqMjAT6HjZPXWB\n5di/w11AZREJdfLZPGPHjh20bduWc+fOcd9996W1rRvuU6RIkbSlbt544w23LTT5zu/v0LluZ+qE\n1snW8xcvXuSjj+xbDL/55ptX7OKZkgItWpilWQxreDupVACOpPt81HEuvUigE4CINAXCgIpOPpsn\nbNiwgVtuuYV//vmHW2+9lRkzZphZ8x7ywgsvUKJECX7//XeWLVuW4/ftOb2H6VumMzR8aLbf8dln\nn3HixAkaN26c4Si//PnhscdyEKRh5IDT06JEpADwAFA5/XOq+rabYxoBfCYiG4GtwCb+O+rMKUOH\nDk37Ojw8nPDwcDeFZ63p06fTp08fEhISaNeuHXPmzKGg+ZXUY4oVK8Yrr7zCq6++yhtvvEGbNm2u\nqBm4YvCywQxsMZDSIaWz9fzp06cZMWIEACNGjMhRLEbeFhERQUREhPtf7Gw7GbAI+B54BRh46XCl\nrQ373JZF6T4PBgZl8cwBoLArz5IL+1ROnjypDz/8sAIKaO/evTUhIcHqsPKE2NhYLV26tAI6b968\nbL8nLilOO//QWS8mZb+z48UXX1RA77jjjgyvp6aqduyoev58tosw8ijc1KfiSkLYluPCIB//drYH\nYe9sr3PZPcWA/I6v+wCTnX023Tvc+XfttNTUVD1z5oweOnRIt2/frlFRUXr8+PEc/fA/fvy4vvnm\nm1qkSBEFtFChQvrVV19d0TlreNaoUaMU0Dp16mhSUpIlMRw8eFCDgoIU0I0bN2Z4j82mumyZPbkY\nhivclVRcWaZlPPC5qm51rg6U6XvaAZ9h78+ZqKojROQpxzc0XkSaA1OAVGA70EtVz2f2bCZlqLPf\nV07t2rWLOXPmsGTJEiIjIzlzJuOdAMqWLUvVqlWpUqVK2p/lypWjbNmylClThuDgYJKTk4mLi+Pg\nwYNs2LCBZcuWsWTJEmw2e+tf+/btGTNmDFWrVvXK92b8KzExkeuuu459+/YxevTotL1XvKlLly58\n9913dOnSxe+3PTZ8j7uWaXElqewAqmNvjkrEPpteVdXnFpb0RlJZsWIF7777LsuXL//P+cKFC1Oi\nRAlCQkKw2WycO3eOM2fOpCUGVwUGBnLXXXfxyiuvcPPNN7sjdCObfv75Z+677z5KlCjBnj17KFWq\nlNfKXr58OW3atCE4OJidO3dSuXJlr5Vt5A1WJJVKGZ1X1UM5DcLdPJlUjh07Rr9+/fj5558BKFSo\nEJ07d+aee+6hRYsWlC9f/orOU5vNRnR0NAcOHGD//v3s37+fgwcPcvz4cU6cOMHx48dJTk4mMDCQ\n4OBgKlWqRJ06dWjdujXt2rX7z8qzhnVUlbZt27Js2TKeeeYZPv/cO2upJiUlUb9+fXbt2sW7776b\n6dbQqnDzzTB3LpTJ3oLHRh7m9aTiKLQBcGk9iFWqGpnTADzBU0ll4cKFdO/enVOnThESEsLgwYN5\n9tlnKVasmNvLMnzT1q1badiwISLCX3/9lbYXfGbOxJ8h2Zac7WXtAYYPH85rr71GzZo12bJlS6a7\nOqrCvn1QrZpZnsVwnbuSiiud7AOAbcDbjmMr8Kw7OnbcfeCBjvovvvhCAwICFNC2bdtqdHS028sw\n/MOAAQMU0Pr162tiYuJV7312wbP6/MLns13W3r17tWDBggrokiVLsv0ew8gKFoz+2gKEpPscAmxx\nRxDuPtydVIYPH542lPett95Sm83m1vcb/iU2NlarVq2qgA4dOjTT+3ad2qWlRpbSk7Ens1VOSkqK\n3nTTTQpoly5dnLg/W8UYhqpak1S2AsHpPgcDW90RhLsPdyaVL7/8UgEVEf3666/d9l7Dv61YsUIB\nDQwM1M2bN2d4z30z79ORq0dmu4z3339fAS1fvryePn06y/vr1VPdvz/bxRl5nBVJ5UXsS6gMdRyb\ngefdEYS7D3cllfnz56uIKKBffvmlW95p5B79+vVLm7sSExPzn2sRByK00qhKGp8cn613R0REaL58\n+RTQ3377zaln4uLM/BQj+7yeVOxlciPwnONo5I4APHG4I6ns3bs3bQXgqzVxGHlXTEyM1q1bVwF9\n5JFH0iakpqamarOvm+mMLTOy9d7o6GgtU6aMAjpo0CB3hmwYmXJXUnFp9Je/yOnor6SkJJo1a8bm\nzZvp2LEjc+bMISDA22tvGv4gKiqKJk2aEBsb+5/hvlGnoqhVqpbLa3OdP3+eW265hcjISG677TZ+\n++03AgOzXqLv9GkoWdKM+jKyz2t71IvIasefMSJyId0RIyIXchqAL3r77bfZvHkzVatWZcqUKSah\nGJmqXbs206ZNQ0R44403mDjRvhlq7Wtqu5xQLl68yH333UdkZCQ1a9Zk5syZTiUUgE6dYNs2l8M3\nDLczNZXLrF+/nhYtWpCamsrvv/9Oy5Yt3RydkRuNGzeOfv36ISKMGTOGfv36ufT8uXPnuOeee1iz\nZg3lypVj7dq1Ls2av/TP3dRUjOzyWk0lXYEjnTnnz2w2G08++SQ2m40XXnjBJBTDaX379mXEiBGo\nKv379+ell14iKSnJqWe3bNlCixYtWLNmDRUrVmTZsmUuL8MiYhKK4SOc7XwBNmZwLlfNU/niiy8U\n0LCwMI2Li8vWO4y87euvv06bJHvDDTfoqlWrMr03JiZGhw0bpgUKFFBA69atq4cOHXK5zF277CO/\nDCMn8FZHvYj0BfoBVYF96S4VAdaoajd3J7qcyk7z1+nTp6lRowZnz55l1qxZPPDAAx6KzsiNok5F\n8cX6LxjdfjRr166la9euHDx4EICWLVvSsWNHrr/+egoXLszRo0dZsWIFs2bN4vTp0wD07t2bTz/9\nlJCQEJfL7tULnnoKmjZ153dk5DVeW/tLRIoBJYDh2DfGuiRGVTNe591i2UkqAwcO5JNPPqFNmzYs\nWbLE7KhnuOTe7+7llkq3MPCmgQBcuHCBjz/+mFGjRhETE5Ppc82bN2f48OG5ZmdSw39ZsqCkv3A1\nqRw9epTq1auTmJjIxo0badSokQejM3KbFQdW0GteL3b230mBwP8u9nj+/HkWLlzI4sWLOXz4MDEx\nMVSoUIHrr7+eTp06Ub9+ffMLjOETvFlTWa2qLUUkBvv6V2DfSwXsbXBFcxqEu7maVJ588km+/vpr\nHnroIb7//nsPRmbkNqmaSuPxjRnccjAPXfeQ18vfsAHKl4dy5bxetJHLeG30l6q2dPxZRFWLOo4i\nlz7nNACrHT58mEmTJhEQEMDbb79tdTiGn5m+ZToFAgvQuW5nS8pfsgT27LGkaMPIkCtDijuLSBHH\n12+IyBwR8ft2ok8//ZSUlBQefvhhatWqZXU4hp+JvhDNx3d8bFkT1uDB0Lq1JUUbRoZc2flxi6rW\nF5GWwLvAh8AQVW3myQCzw9nmr7Nnz3LttdcSFxdn+lIMw8jTvD75Ebi0yfrdwHhVnQ8E5TQAK40b\nN464uDjatm1rEorhd9asMUuzGL7HlaQSLSJfAQ8DC0SkgIvPAyAi7UQkSkR2i8igDK4XFZF5IrJZ\nRLaKyOPprr0gIttEZIuIfCsi2U5qCQkJjB49GoBXXnklu68xDMscOQInT1odhWH8lyvNX4WAdtg3\n5tojIuWA61V1sdOFiQQAu4E2wDFgPfCIqkalu+dVoKiqvioi1wC7gDJAaWA1UFtVk0Tke2C+qk7N\noJwsm78mTJhAnz59aNiwIRs3bjTDOg3DyNO83vylqhexz6i/U0SeAUq7klAcmgJ7VPWQqiYDM4GO\nlxeFfbY+jj9Pq2qK43M+IEREAoFC2BNTtowbNw6wT3o0CcVwReTxSKtDMAyf5crorwHAt9hrDKWB\n6SLyrIvlVQCOpPt81HEuvTFAXRE5hn2nyQEAqnoM+Bg4DEQD51R1qYvlA/D333+zceNGSpUqxYMP\nPpidVxh51NL9S3nghwdISU3J+mYPWrkSlmbrX79heJYrfSK9gGaqOkRVhwDNgT4eiOlOYJOqlgca\nAWNFpLCIFMdeq6kElAcKi0iX7BTw1VdfAdCjRw+Cg4PdE7WR69lSbQxcPJCRt48kMMC5fU48JSDA\nrEps+CZX/s8Q/h0BhuNrV/9ZRwNh6T5XdJxLryf2dcZQ1X0icgCoDVQG9l9ab0xE5gA3ATMyKmjo\n0KFpX4eHh6etrXThwgW+++47APr08URONHKrqZFTKRJUhE51OlkdCq1aWR2B4e8iIiKIiIhw+3td\n6ah/EegBzHWcug+YrKqfOl2YSD7sHe9tgH+Av4BHVXVnunvGAidVdZiIlAH+BhoA1YGJQBMgEZgE\nrFfVsRmUk2lH/aXNlG655RaP/IUauVNcUhy1xtRi9kOzaVbR56ZmGUaOuauj3umaiqp+IiIRwKWd\nq3qq6iZXClNVm6OTfzH2preJqrpTRJ6yX9bx2CdWThaRLY7HXnHUTv4SkVnAJiDZ8ed4F8tn/Hj7\nI0899ZQrjxp53Ji/xtAyrKVPJJRVq2DvXujZ0+pIDONKeWqV4sjISBo2bEjJkiU5duwYBQoUyOBp\nw7hSTGIM8SnxlA4pbXUo7Npln59imsAMd/J6TUVEgrFv1tUS+7Df1cA4VU3IaRDeMn36dAAeeeQR\nk1AMlxQpUIQiBYpkfaMX1KplPwzDF7nSp/IDEANMd5zqAhRXVWuWZ72KjGoqNpuNsLAwjh07xtq1\na2nRooVF0RmGYfger9dUgHqqWjfd5xUisiOnAXhLREQEx44do2rVqjRv3tzqcAwjWyIj4bvvYMQI\nqyMxjIy5Mk9lo4ik/TQWkWbYR2b5hUtNX926dTMz6A2/VbEi3Hef1VEYRuZcaf7aCdTCPqMd7PNN\ndgEp2Edu1fdIhNlwefPXxYsXKVu2LDExMezatYuaNWtaGJ3hD2ypNh6Z/Qhj7xrrE53zhuFpVjR/\ntctpYVb55ZdfiImJoWnTpiahGE6ZvHkyJ+NOEloo1OpQDMOvuDJP5ZAnA/GkmTNnAtC1a1eLIzH8\nQWxSLEMihvDzIz/7VFPpP//AE0/AwoVWR2IYmcv181RiYmIIDQ0lKSmJo0ePUr58eYujM3zdWyve\nYt/ZfUzvND3rm70oKQl274Z69ayOxMiNrGj+QkSqAP/409yU+fPnk5iYSMuWLU1CMbIUfSGaMevH\nsPHJjVaHcoWgIJNQDN/n6s6NL2FfnRgRaeXYr96nzZo1C8AscW845fD5w7zR6g0qFa9kdSiG4Zdc\nai45k3AAAAv3SURBVP4SkR7YVyZeqaoHROQ+Vf3JY9Fl06Xmr7i4OEJDQ4mPj+fIkSNUrFjR6tAM\nI1uSk6F2bdixA8xiEIYneH3nR4drgSTgRRFZDjTOaQCetGDBAuLj42nRooVJKIZfCwyE1atNQjF8\nn6s7De0HZqnqDBEpBVi/scRVmKYvI7cQgXLlrI7CMLLmak3le+BSV2FVoKx7w3Gf+Ph45s+fD5ik\nYvi/5GSrIzAM57iUVFTVpqobHV+vV9V3PBNWzi1dupS4uDgaN25MWFhY1g8YeVbk8UgSUnx7QGPD\nhrB/v9VRGEbWXK2p+I2ff/4ZgI4dO1ocieHLYhJjaPdtO3b+b2fWN1to0yaoUsXqKAwja672qfiN\nX375BTBJxbi6kWtG0rZqWxqVa2R1KFcVFGR1BIbhnFw7ox6gSpUq7Nu3z6eW2jB8x9ELR2nwZQM2\nP7WZa4tda3U4mTp1CkqWhIBc265g+AKrhhT7lY4dO5qEYmTq9eWv8/SNT/t0QgF45hlYvtzqKAzD\nObm2+Qvg3nvvtToEw0dFnYpi8b7F7H5mt9WhZGnmTMiFDQpGLpVrm79KlPh/e+ceK1dVxeHvRwsR\n6Ls8wkMLxiCpQKG8qjQUMFKwQrA8hCJpTDQGARMBI0YQoWjREGL4A7GxVsMzKhRqDZZHbgUEUii3\nLYW20FAeNuCF8GjRIqVd/rHXoafDnHtn5s7cuXO6vuRk9uzZ++y19jpz1tn7nLP2aHp6ehg6tNR+\nM2gQM2P9xvXsPyJeig0CiOmvPpk2bVopHcrixYvbLUJLGSj9JLXFodSr3/r18M47rZGlFcTxGZTW\nqZT1qa+yH9Sh3/bceSfcc09rZGkFYb+gfJfyztSpU9stQhD0m8svb7cEQVAfpR2pDB8+vN0iBEEQ\n7HCU9kZ9u2UIgiDoNJpxo76UTiUIgiBoD6Wd/gqCIAgGnnAqQRAEQdPoKKciaa6kf0ta0Ue5oyVt\nljQ9l3eKpNWSXpD0o9ZLWx/91O1lScsldUta0npp66cv/SRNkfSupGd8uzL326C2HfRbv463n5c5\nwXVYKakrl9/x9vMyRfoNavvVcGxe7rI/I+lZSR9JGuW/1W87M+uYDZgMHA6s6KXMTsDDwEJgei5v\nLTAO2BlYBhzcbn2aoZvnvwSMbrcO/dEPmAIsKNB5UNuuP/qVyH4jgeeA/fz7HiWzX1X9OsF+tZxb\ncmW/BjzUH9t11EjFzB4D+nq/+BLgL0BPLu8Y4EUze8XMNgN3AYPq7ch+6AYgBvmos0b9qj15Muht\nB/3SL8vvdPvNAO42s/Ve/i3PL4v9ivSDQW6/Go/NjPOAOz3dkO0GbUc0gqR9gTPM7Dds/wfeD3gt\n9/1fntcx9KIbgAGLJD0l6TsDL13TmOTD8L9JGu95HW+7HNX0g3LY7yBgjKQu1+MCzy+L/Yr0g3LY\nD0m7AqcAd3tWQ7Yr2xv1vwYG5ZxtE6jULe9YjjOz1yXtCTwoaZVfnXQSS4FxZvZfSacC95L+yGWh\nN/3KYL+hwETgJGB34AlJT7RXpKZSVT8zW0s57AdwGvCYmb3bn52UzakcBdyltIjKHsCpkj4C1gP5\nher397xOoppum81sgZm9DmBmb0qaTxq2dtRBbWbv59L3S7pZ0hjKYbtC/czs7TLYj3QV+5aZfQB8\nIOkRYAIlsR/F+q0tif0AzmXb1Bc0aLtOnP4SBXPTZvZZ3w4k3Xv4npktAJ4CPidpnKRdSJ23YMAk\nrp26dZO0m6RhAJJ2B04GVg6YxPVRqJ+kvXPpY0gv5r5N59gOGtCvLPYD7gMmSxoiaTfgWGAVJbEf\nBfp1kP160w1JI0kPk9yXy27Idh01UpF0B3ACMFbSq8DVwC6AmdmciuIfhwowsy2SLgYeIDnSuWa2\namCkro1GdQP2BuYrhaYZCtxuZg8MgMh1UYN+Z0m6ENgMbAK+AZ1hO2hcP0piPzNbLWkRsALYAswx\ns+e9bsfbr0g/SQcyyO1X47nlDGCRmW3K6jX634swLUEQBEHT6MTpryAIgmCQEk4lCIIgaBrhVIIg\nCIKmEU4lCIIgaBrhVIIgCIKmEU4lCIIgaBrhVHZgJPX51q+kOZIO9vSPG6i/sXEJ66MWeRrc78d9\nMBBIukbSSf3cxy6SHvRw5mf3c18TPLRMS/B4WhObvM+WyhwUE++pBDUjaaOZDa+zzgYzG9EqmQYj\nknYys61tlmEScK2ZndyEfc0EjjKzS+qoM8TMttRYtgu4zMyeaVTGKvusW+agOcRIZQcmG0UoLSDV\nJenPklZJujVXpkvSREmzgV39yvfWivq7S3pI0tNKixWd3ke7u0laqBSxd0V2Je3tLFaK9np/FtrE\nZbjR85+TdJSkuyWtkTSrUp8q7c3zWFtPSFrr+s6V9Lyk3+fK3SxpidJCRVdX9oGnz3OZV0i6Pt+2\npBskdQOTKtr/tu+32/v4U55/rzzaraTv5vp1nnwRNknXKy0KtUzSr6roNlrSfO/3xyUdohTY8Fbg\naLfXgRV1JnhfLPN+HFlFz7GS1kkaClwLnJONetx+cyU9KWmppNO8zkxJ90l6GHioos1xfmzd5v3+\np6wfKsoV2WCdpJ95e8slHeT5n5BF0s6VMlc7LoIW0Z/FX2Lr7A3Y4J9TSOst7EOKD/Q48CX/rQuY\nmC9fpf4QYJinx5LWYKBaHc+bDvw29304KcTFP4GxnncOKSxEJsNsT3+fFNRuL1KoidfwBZKqteX5\n84A7PH068B4w3r8/DRzm6VH+uZO3eUi+D7x/XgHGsG3BtNO9zFbgzIL2R+fSs4CLPL0X8AJpEaXV\nwMicvNO9ndW5uiOq7Psm4CpPnwh052xatCjYcmCyp68Bbqxi67HAS56eCdyUq/9zYIanRwJrgF29\n3KuZHhVtjvM+muTf5wKXVmm3yAbrSPHuAC4khUnpS5abqukfW2u3GKkEGUvM7HVL/85lwAF11BUw\nW9Jy0hXqvpL26qX8s8BXJM2WNNnMNgKfBw4hhQ7vBn4C7JursyBXd6WZ9ZjZh6RV9z5dg4x/zdV/\nwzwuFWk1vwM8fa6kpUA3MN63PEcDXZYiC28FbgeO99+2APcUtH2opEeUlnOdAXwBwMx6SHGYukgn\n2Pcq6r0HbJL0O0lfJ8UMq2QyaVSCmXWR1vwYViAHkkaQTvrZ/ac/5nSolZOBK9xOi0nOPYtm+2AV\nPTJeNbMnPX2by15JbzaY759L2Waz3mQJ2kBHBZQMWsr/cuktVD82iqKcnk8Kx3+EmW2VtA74xNRG\nhpm96NMsXwVm+XTJvSRncVwf8m2tkHVrpaySrgOmpaYsuwHca31JBwCXAUea2QZJ8wp0KOqDTe6Q\nq/EH0ohmpdJc/5Tcb4cBb1Fl8SNLAf2OAb4MnA1c7OntitUoXy18xLYp8UL7OWea2YvbNZzu4/yn\njva2k70GG2R2yx+f6kWWoA3ESGXHpt4T0Ic+x15ZfyTQ4w7lRNJUR2EbkvYhnYTvAG4gTS2tAfbM\nTgaShmr71RFrQQBmdqWZHZFzKFXLVTACeB/YqHQvp9qTQ0uA4yWNkTSEtPTq4l72mTEMeMPn+s//\nWIjkMKYCRwA/lJTvN5RCrI8ys78Dl5IcUCWPAt/08icAb1pu7ZZKzGwD8I6kzHlfAPzD0y+T1u2B\n5MQyNpL6J2MRaRoyk/PwovYq+IykYz09w2XPU4sNKimSpVLmYIAIp7JjU3RlbQXpOcAKbbuRn/12\nO+mm8HLSCW5VQf2MQ4ElPmXxU+A6S2tgnwX8UtIy0vTHF/uQszdZi8pUrWNmK0jTfqtIUzOVjyeb\nmb0BXEFyJN3A02a2sAYZryI5pEd9/yitTzEH+Jbv9zIge2gg29cIYKH36yPAD6rs+xrgSC/zC9K9\nhL6YCdzg/TyBdFMbkoO/0KefxuTKdwHjcze9ZwE7+8MKK3P1+2INcJGk54FRwC2eX4sNivo3L8uz\nOVkqZQ4GiHikOAj6wO+FnGZmr7Rblk7FR2ELzezQdssStJYYqQRBL0h6AFgeDqUpxBXsDkCMVIIg\nCIKmESOVIAiCoGmEUwmCIAiaRjiVIAiCoGmEUwmCIAiaRjiVIAiCoGmEUwmCIAiaxv8ByM8POMbp\nNTsAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "None" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from ipywidgets import interact\n", "def generate_plot(a_0=1.56):\n", " x, dxda, ddxdda = run_sim_var(a_0)\n", " x_1st_order = np.zeros(N)\n", " x_2nd_order = np.zeros(N)\n", " for i,a in enumerate(a_grid):\n", " x_1st_order[i] = x + (a-a_0)*dxda\n", " x_2nd_order[i] = x + (a-a_0)*dxda + 0.5*(a-a_0)*(a-a_0)*ddxdda\n", " \n", " fig = plt.figure(figsize=(6,4))\n", " ax = plt.subplot(111)\n", " ax.set_xlim(a_grid[0],a_grid[-1])\n", " ax.set_ylim(np.min(x_exact),np.max(x_exact)*1.01)\n", " ax.set_xlabel(\"initial semi-major axis of outer planet\")\n", " ax.set_ylabel(\"$x$ position of inner planet after 10 orbits\")\n", " ax.plot(a_grid, x_exact, \"-\", color=\"black\", lw=2)\n", " ax.plot(a_grid, x_1st_order, \"--\", color=\"green\")\n", " ax.plot(a_grid, x_2nd_order, \":\", color=\"blue\")\n", " ax.plot(a_0, x, \"ro\",ms=10)\n", " plt.show()\n", " return \n", "interact(generate_plot,a_0=(1.4,1.7,0.01));" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.10" } }, "nbformat": 4, "nbformat_minor": 0 }