{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Stability map with MEGNO and WHFast\n", "In this tutorial, we'll create a stability map of a two planet system using the chaos indicator MEGNO (Mean Exponential Growth of Nearby Orbits) and the symplectic integrator WHFast (Rein and Tamayo 2015).\n", "\n", "We will integrate a two planet system with massive planets. We vary two orbital parameters, the semi-major axis $a$ and the eccentricity $e$. Let us first define a function that runs one simulation for a given set of initial conditions $(a, e)$." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def simulation(par):\n", " a, e = par # unpack parameters\n", " sim = rebound.Simulation()\n", " sim.integrator = \"whfast\"\n", " sim.integrator_whfast_safe_mode = 0\n", " sim.dt = 5.\n", " sim.add(m=1.) # Star\n", " sim.add(m=0.000954, a=5.204, M=0.600, omega=0.257, e=0.048)\n", " sim.add(m=0.000285, a=a, M=0.871, omega=1.616, e=e)\n", " sim.move_to_com()\n", " \n", " sim.init_megno()\n", " sim.exit_max_distance = 20.\n", " try:\n", " sim.integrate(5e2*2.*np.pi, exact_finish_time=0) # integrate for 500 years, integrating to the nearest\n", " #timestep for each output to keep the timestep constant and preserve WHFast's symplectic nature\n", " megno = sim.calculate_megno() \n", " return megno\n", " except rebound.Escape:\n", " return 10. # At least one particle got ejected, returning large MEGNO." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's try this out and run one simulation" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "2.097652943528103" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import rebound\n", "import numpy as np\n", "simulation((7,0.1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The return value is the MEGNO. It is about 2, thus the system is regular for these initial conditions. Let's run a whole array of simulations." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "Ngrid = 80\n", "par_a = np.linspace(7.,10.,Ngrid)\n", "par_e = np.linspace(0.,0.5,Ngrid)\n", "parameters = []\n", "for e in par_e:\n", " for a in par_a:\n", " parameters.append((a,e))\n", "from rebound.interruptible_pool import InterruptiblePool\n", "pool = InterruptiblePool()\n", "results = pool.map(simulation,parameters)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On my laptop (dual core CPU), this takes only 3 seconds!\n", "\n", "Let's plot it!" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAFLCAYAAACp0ulDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXd4VGXaxu8nCYGE0AKhNwFBUBEEEQUFxYJlFSso9s9d\nXWVX13XX3bXs6rrr6rpWLGBfG2ChrIJioQgKAtI70lvoJJAEQvJ+f8ww59wPyWQGJmWS53ddXJx7\nTnvnzJl5c54qzjkYhmEYRlUhobwHYBiGYRhliU18hmEYRpXCJj7DMAyjSmETn2EYhlGlsInPMAzD\nqFLYxGcYhmFUKcp84hOR/iKyTERWiMgDRazvIyJ7ROSn4L+HynqMhmEYRnhEJCH4Gz2umPUviMhK\nEZknIl18r4edA8qCpLI8mYgkABgKoB+AzQBmichY59wytelU59xlZTk2wzAMIyruAbAEQG29QkQu\nAtDWOXe8iJwO4FUAPaOYA0qVsn7i6wFgpXNunXMuH8AIAJcXsZ2U7bAMwzCMSBGR5gAuBvB6MZtc\nDuC/AOCcmwmgjog0QuRzQKlS1hNfMwAbfHpj8DVNTxGZKyKfi0inshmaYRiGESHPAvgDgOJKfxX3\nWx/pHFCqVMTgljkAWjnnuiLwSDymnMdjGIZhBBGRSwBkOufmIWCdi8RCV6GseGXq4wOwCUBLn24e\nfC2Ec26fb3mCiLwsIunOuV3+7UTEiowahlGpcc7FZMLIkGpuBw4dza6ZzrnG6rVeAC4TkYsBpACo\nJSL/dc7d5NtmE4AWPn34tz4ZJcwBZYGUZZFqEUkEsBwBx+YWAD8CuM45t9S3TSPnXGZwuQeAUc65\n1kUcy72P9hGfe9ChoaRHJA05indQNAOncxzOyF5FBjmF5RPswFVoEKshlTuRvJ8k9WfXoaP6XnoM\nOPg06dHJ95PeNKsP6eanTSGtvwmFS84NLY/pNKpSfT5A1bznSsJ/D3wx4QxatzuvgPR1J6Txzh2/\nJfnT1LNJN6/JBraG3SaTvuHOeqHlK7r0weg7n4rZxCci7sOEyH8vD3Nd4YqwYxCRPgB+r4MRg5Pi\n3c65S0SkJ4DnnHM9I5kDyoIyfeJzzhWIyBAAExEws77hnFsqIncEVrvhAK4WkV8DyAeQC2BgWY7R\nMAyjMpJwNI6twsg39f+OO+fGi8jFIrIKwH4AtwLFzwFHMbJjoqxNnXDOfQGgg3ptmG/5JQAvxfq8\nsXzC0xzNE54BtG7JetVq1voJbO0MfmI7ric/sSUOfzfs+bIPRmfdmLLxQFTbV1TOzuYn38x9WwEA\n9Z/+Fh3uPxfLm7xXHsOqkPgfbS66+Adad9ndzUmPvWJjsfsCQLc+U0kPynuOtNt8Kun7x4wPLddt\nVTeC0ZY/zrkpAKYEl4epdUX+6BY1B5Q1ZT7xGUXTEanlPYSYYu+n4tPtzOPKewgxpTJ+RrHkqJ74\nKik28VUQOlWyL629n4pP90o28VXGzyiW2MTnYROfYRhGFSChQiUUlC828VVCjts0KLS8ptmIchxJ\neHZ9rarStWFfqf6e/rnrVaRHgH18nwyZH3b/R7tdqfZnH4zevl/LGqHleM6dmVrr6ZI3OgaSk1kf\nPFiqpys71Ic+bmh4n15J+2/P30b662bPkN55b0ZoObtGdiQjjAp74vOwic8wDKMKYBOfh018hmEY\nVQCb+Dxs4quERGPevOJxjioe/dDyWA+nWLo/wuHiE9V6bV5s+ioXIdi09EbSUq8J6RGNnyKdP/zN\nqMbXqv83oeW1Ue1ZtdCmzTSV171vn7fcpjWvW722NEZ09PjvuZXTOQG9fS82jZdEQiLrSRtnFn8y\nAPm5XgWHgoOcLB8LbOLzsInPMAyjCmATn4dNfIZhGFUAm/g8bOIzDMOoAtjE51FpJr5wfgWjeMrS\np6eZ+N72sOt1uPiq1Xmk1zWuQ3pGPfbpab7/54KIxwYAa9dGtXmFQftGX3i3C+l7bpxXqucP992r\naD69cPzjw8Wk34ly/0Llprs2swbpXu/x55I020uXSKwW+1lKxBL5DmN/AxiGYRhVikrzxGcYhmEU\nj5k6PSrNxFeRTZsX5/6d9PiUh8tpJBWLhhmst4W3fGLRuPWkzzzwPumSegxs2RrhwOIcbdD6bYxN\nmzVXXkx6//Hji9kyvklvl0667+X5pCePzYrqeNKAOy6sm7+S9CFfv7+C/Cj6AUWITXwelWbiMwzD\nMIrHJj4Pm/gMwzCqADbxedjEZxiGUQWwic/DJr4ywHx6HgM/7R1aHnnltKj2/XkN66anNFZb7D7K\nUVVujjWIXadH/G8t+7a2fnIa6euvmnWMZyw79Ht79PUTQ8tzxmygdZM/i86np4/9bN5m0hOvO5P0\nGc95Rfuq16oe1bkiwSY+D5v4DMMwqgA28XnYxGcYhlEFsInPI64nPn+agJkTKyba3PNIQ+/bdyKO\njXu6N1KvLD3GIxpFoU2lL5w9gPTo6veX3WCiRN9/M7/tTbrnuWxu/+vtXrWWknoxDJx6Ken9H/5E\n+rNX2LT5m+k7SK+9U53hfN/9XAqdj60Du0dcT3yGYRhGZNgTn4dNfIZhGFUAm/g87FIYhmEYVYq4\nfuJ7bu7XoeX25TgOo3i0W+HPo9aGlscd47F+3/gk0iMwmfQFgxuQnvg++1iMo6P6Ek5XuOa1U0l/\n/Ev2dZUn+p4Z8pcZpGeXsH04Rvb5jPQ1b3fnDZSPb/st55NOvZl1zgNvhZYP7DsQxUgiI1ZPfCJS\nHQEXaDICc8jHzrlH1Tb3AxiMgLeyGoCOABo45/aIyFoAewEUAsh3zvWIzcgiJ64nPsMwDCMyEmIU\n3eKcOyAi5zjnckQkEcB0EZngnPvRt83TAJ4GABG5FMC9zrk9wdWFAPo658ot8dYmPsMwjCqAJMYu\nrNM5lxNcrI7APBIuDvU6AB/6h4JydrPF9cT36uwtoeVnynEcRuRkb9gb8bYDl91MeuQJ3ArUldAp\ntiKbNmtwT1Lk5RW9XUVk8oUjSWduK6eBHAW1m9ci7VS1n3++cXJo+cH/W0jr6nLfY+xRt/JHN2vD\nKTOlDv9KTf3qLNJJNRJDy6XSiDaG+QwikgBgDoC2AF5yzhVZrkdEUgD0B3C372UH4EsRcQCGO+de\ni9nAIiSuJz7DMAwjMmI58TnnCgF0FZHaAMaISCfn3JIiNv0FgGk+MycA9HLObRGRDABfichS51x0\n9QuPEZv4DMMwqgCRmDpn7N2HmXv3R3xM51yWiExC4KmuqIlvENjMCefcluD/20VkNIAeAGziMwzD\nMGJLJE98Z9SrhTPqeebgFzcc2R1aRBogEI25N2jKPB/Av4rYrg6APghEdx5+LRVAgnNun4jUBHAB\ngEf1vqVNXE98TerUKHkjo0IxaXTkPj7t00tM5PXjL+Zwck3XLqznxrYR+TFx+YJbSI9s/3a5jONo\naL3sLtKZ6S+X00ii59uPwwcSbt+bW+y63erWffQ1Lrr3t18uRjhO4ewbvKvOlVwzObScVKNa2GMd\nDTEMbmkC4J2gny8BwEjn3HgRuQOAc84ND243AMCXzjn/G20EYHTQv5cE4H3n3ESUMXE98RmGYRiR\nEcN0hoUATi3i9WFKvwPgHfXaGgDqT9KyxyY+wzCMKkAsg1vinbie+FavtcajVYmCAtYtVt1IunP9\nM0iv2su2zWYnDyd9+n1sqvr0/vCmqliS1bxNzI6lE6jWzuhD+rieU2J2LgDIKcgpeaM45dCBQ6Fl\nfV03/sjXddPoFWGPNWjLfaQ/2ruA9IyPV5JOrZ8aWq5eKxmxJpZ5fPGO1eo0DMMwqhRx/cRnGIZh\nRIZYe4YQNvEZhmFUAczU6RHXE9/B/QfLewhGObKz87v8wpy6JOvU5e4Ms7lYPrJmZ0Z8rgF/PZ70\nmEdXFrNlZExIfeSY9veTvPwC0q07hI8Ob5jBetuRqVphKTj77eh2qMDoqSCzldcFfRTW0rp6KfzE\nlJpRk/R9z/I94ppyibL7ngvfQ8YVeF5FVxh206PCgls84nriMwzDMCLDJj4Pm/gMwzCqAGbq9Ijr\nia9Oc3+59MjNVkZ80p4tSfjFbzqQXtuoLen7vhpF+mJ1vG592Ob35Ijiuzkcq2mzNMlXps2Sft6i\nNW1qliw9tv3Lk75ZfyQ9ufZTpD+6fmax+979OndrqN+uPh+7TTrpuj/z/flwFtfA/NMnnA6xd2NW\naDm7Vnax4zhaYpXAXhmI64nPMAzDiAwzdXrYxGcYhlEFMFOnh018hmEYVQB74vOI64mvbs3Yl/Ux\nKg66ZNSLL5xJemO/QaTP+vgt0mu2c2mtZwaxT2bT/K3HNsAKwsDvB5AeeeaYchpJ+TNwxhWkR/Yc\nTVr79I7Yf0K/0PKy7i1pXeK/x5PesYJ9wg1PYp9xi8u4e8gtt7YjnVSDf37rtKgdWvaXL4sV9sTn\nEdcTn2EYhhEZVrnFo8yvhIj0F5FlIrJCRB4Is91pIpIvIleW5fgMwzCMyk2ZPvEFGxcOBdAPwGYA\ns0RkrHNuWRHb/QvAl+GOl7k19iG/RsXh9DO5GecFfc4jffH//kv6m8x80o32sKlz19DLSaf8Z0LE\nYxm0/0nSI2oW+zdbmVOSaVObjL/+gk3G5/f/PqbjGfgFf04j+3991MfSY7/xDq7O896wPXwuZdos\n6XgPvtqJX7j4m9DiocX9aNVZPZqTnjJ9PekPZnFpoJ6L+NCbm7CpPXkqpzOkt6kXWk5Iiv0zifn4\nPMra1NkDwErn3DoAEJERAC4HsExt9xsAHwM4rWyHZxiGUTlJMB9fiLKe+JoB2ODTGxGYDEOISFMA\nA5xz54gIrTMMwzCODnvi86iIwS3PAfDbkezTMgzDOEYsqtOjrCe+TQD8McLNg6/56Q5ghIgIgAYA\nLhKRfOfcOH2wfeOXh5aXoBCdEPsQYIOpUYN1Xl7pnat+R/aJ1MnaR/rxM7nafcJHXFJKh4u3/Xgu\n6Qat64HZQsrvD7r3Rw5N71nUgCsITdZdRXpzq09Ir9lbih8ajs2np9E/1e8qn96x4tQJnO9Dv+s9\nrs025o7TSU+bxT9dtevyl6PPpbVIjxrM5dBqP8bpErmLt6Fg7W4AwNwGueEHfhTYE59HWU98swC0\nE5FWCPzKDAJwnX8D51ybw8si8haA/xU16QFAtd6tQ8udpleOnCzDMKomScfVQ9JxgT/GOp/QG4s+\nmRLbE9gTX4gynficcwUiMgTARARSKd5wzi0VkTsCq91wvUtZjs8wDKOyYk98HmXu43POfQGgg3pt\nWDHb3hbuWHVa+Lsz2BNfaaD/8hg/thfpcy+cXmrnzs85RPqC76aS/jKjIWm5bTDppNfeIH1XGzaF\np2zg9IdwfLOGTWwVydSpP6N1+3aT1vWNxlzGqbHj8VPsB1VKHOtPt95/9S/6kx5xx5LQcloTNlX2\naMRB5u/dzukyn6zMIt25F6+ftIbTrxKrJZLevcm7x3Ibx97UiURLYD9MRQxuMQzDMGKNPfGFsInP\nMAyjCmBRnR727GsYhlEVSJDo/xWBiFQXkZkiMldEForIX4vYpo+I7BGRn4L/HvKti6hsZWkS1098\nWxduK+8hVHr0rX/DrezT24zYcfWwrqQHNeaSZY2qqdv1EPsAn1j8BWld9ml417NIL+7CfuGF4BpT\n/ve+YNDNtG7kbfejojJ6Jad9DFTrc677d9kNpoIzstmzpK/Nf94TyffQupoXcomxAc25pNmvPnmF\ndOZVfflcf36b9O617Deu1SQttJxSL6XYMZc3zrkDInKOcy5HRBIBTBeRCc65H9WmU51zl/lfiLRs\nZWkT1xOfYRiGESExNHU65w5H7lRHYB4pKgK/qBNGWrayVDFTp2EYRhVAEiTqf8UeSyRBROYiEE7/\nlXNuVhGb9QyaQz8XkcPVwIsqW9ksVu8xUuL8ia9ypPmVZTWUY2VzLG2bio/v4MoqV6nvXY9dd5LO\nTKpNetpXk0g31+aidatInnjy2aQX4m3S/rvrXwu5cstxqDjon6drB2iLEzN5zN7SGwyA03fdRXpm\n+sulej4/F97EKS45L95AumnmLtIj2r9N+s2lY0PLt0+6hNZt/f2npNuewN09qteqTnpurYOk//Dy\nCaRrrGNTZ/YWz0Sd26h80hmmrt6BqWt2lLidc64QQFcRqQ1gjIh0cs4t8W0yB0CroDn0IgBjALQv\n6ljlQZxPfIZhGEYkRJLA3qddBvq08zrJ/2PSijBbA865LBGZBKA/gCW+1/f5lieIyMsiko7IylaW\nOmbqNAzDqAokSvT/ikBEGohIneByCoDzoXx0ItLIt9wDgDjndsFXtlJEkhEoW1lkScrSxJ74DMMw\nqgKxC25pAuCdYIRmAoCRzrnxqvTk1SLyawD5AHIRDC4urmxlrAYWKXE98TU60W/Pj23V9rIk1j49\nv2/qBtWxWt/6uoN1eXLmeVxSbM/73GWgZRr7SF6YP5b0kh1cIurkjJqkH8zl9JfHPx9J+vI/sedu\n7L/WhJZ7N21A68rcNhMF5Z2mHI1PT3vpf/N4K9JDH1oXdnssOZfk96d/Szrr3WdI3/HnpqS3/4Pv\nqfz2Xvf4BV3+Quv2TOT7MfUDLvV23zk89sY1G5Mec8WppM96mrtYtO3uja1xq3TEmljV6nTOLQRw\nahGvD/MtvwTgpWL2P6JsZVkT1xOfYRiGESFWuSWETXyGYRhVgQQL6ThMXE98mYu3l/cQYoI239yo\nzJPRmiP9f9dVJFNmSaTU47yOFjv5yqxN5giz52dnkj6wn8PHCxqw6fTx6tzYdlI3Pn7zf/xA2r/2\nurFc1eVplC89z/ZC52dMPVCOI4mOfteyCe+bUZxe8KIybWr0M8vpd3AloZncAOGI7Yc9wfk4vxjZ\nhPS7g58KLT/xb248++RtnL6wux2nL+Qe4say87bz/TrkmzWkW3fi1Is1i7z7eXtC7FNOrFanR1xP\nfIZhGEaEWHeGEDbxGYZhVAXsiS+EGX0NwzCMKkVcP/E1PtlvI99d7HYVjZJ8eu/G0C9X1n/jdT6J\n9QJ2jeGad70u1h/dyOX9fvyS/T2XvMGh5s/O+ZL0oUMFpHes3El6q/LxHRi3kPS0izm0veO8QtL+\nazd9zHpax97Assfv19P305LvuBTbQ99yF4rRfw1fjeNYOeeKOqQnjfb8Vd98xJ9xD1XerP1fPiQ9\n7g9nkr5oH3dFzzyO/Wqu9lOk//RyR9JP3sUpYy+dxx1BtvbpFFpe9sV8WnfWVZyucE59vr+yD/L9\nmNWOy9ytv6Me6TrNueRe7WaeTk3nY8eCWKUzVAbieuIzDMMwIiSCWp1VBZv4DMMwqgLm4wsR1xNf\nvDai1befNm3WVFaOHC5IUqHp+PR5pBf05+oUfvOmNtH96p9czGH7n7nqRg3VbWG/qtSS1pArtSQn\nJfL+915L+oOxXCJw8WvdSG8bvTK03OR53veHNq+jonBEiP/FU0mPViH+mktz/0H6s5QHSevPKWfh\nOaRrqB/UhnvYhOdGjwkt75nL+97xNXeSePdVdllkvfo5aa61cyT6Woy5kj/TJhM5nWHq9rWkayR6\nzY8/u4Lv5Q+XzyH97PdbSH/Vg89V+08/s07lxsrVlN7pM9XvSyvhQzsKzNTpEdcTn2EYhhEh9sQX\nwiY+wzCMqoA98YWI84mvfBrRNlf9gjceY8VifTvGk2lTU7BxV8kbFcNve3AVjT8OvIL03IKNpM95\nZTKfW0XVzfpuLenNl3Fll+y8Q6Sf6s6FqP/UtX9o+bPCUuzAexSck/2n0PKkWv+idZfs5trALYe/\nSHrqB6tJa9OmRt+fqSdPCrvBVaN7kPYbiet15X2H9+bqJy3f7k76o1tmhx3btQuu5/2ncuHojYWq\n2POrvySZnM2Rvv3Hec3Bk2/iiM9qauJIVabKIT8vIT3sCTaF/urPfH9nbWJzZr3WXnR3TWW2jwVW\nucUjzic+wzAMIyKsVmcIm/gMwzCqAmbqDGETn2EYRlXAnvhCxPXEV6txmk+VXReCY/XpVWaSBt/A\nL9z+U9Eb4kjf0a87c/j4qLWTSd/2FpeBqVGXuzlotAf4tol8PO0TvOtkbmo6MvWRsMc/FvTYfv6h\nD+mbH+GqNtO/Ysev9uv5GZl0N+mn1PoZSt/11+akX3l0I8JxzgBVmeWfHMZ/zYncDPaye7zqKil/\nv5PWpf+Hm9a+onx6fS7lSi2TP2O/WNo7nC4z/DYuHTSi2QukB+VyisxT33xD+vJWXirG02u+p3W/\nW5VLeph6ghq9iNOrnn+Tr8u9u/g3Knc3d6BOSPImJimNScomvhB2JQzDMIwqRVw/8RmGYRgRYj6+\nEHE98WVn7g8t/+XVTrTun3cu0ZsbpYA22XX+8G3Sf4niWAnKAPHyPE6NqJnBId75ufmk8/aw6UjT\nqxnvv3DrPtKDJ/yPtDLaEgMebkd6zN9XhT13Sbw4nc2L7b5m0+b2OX1JZ3SbHFq+9n1umDpq8EzS\n+ufu7IvTSL9cgmlTM3kMN0mVsWzavDjncdILdnkpAw0LN9C6H3/bn3T16R+QHvUomwuveK8L6QMP\nvE262bk8lrOGn0p6RMp9pN/ccRvpjmO8SjL/6cJFqT+6bjLpR9SFLVjMpvL31Zej2pT9pDNO4PSZ\n2mnJoeWaNTlVIiaYqTNEXE98hmEYRoTYxBfCJj7DMIyqgJk6Q9jEZxiGURWwJ74QcT3xie8PmBtO\nblD8hhWMk7fdQnphw7dJDzrEJac+VOHpz/33FNLDHueGmctKt88oof+GnDLwStKZzy8g3eHly0LL\nhT04fSHxgArZz2hI+vgs9uGtXRW+PFq9Vhx2P2opV/7P28vHe+AXHAq/EV5KgW4WvGgc+/S0r/O9\nT7j01qJM9u/ohqjP/YEr+Q8cy+kN31w4mfQO37L26TXM4LFs2876wJ4DpK8++Czpj5N/R/qdj9jP\ndvM13KWgYBH7ts796CPS/+7TJrR84SjuHPHYHYtJa199x7Xsh/39a++Q1p1NzlPNXjP2cfpDijr+\n/jenkC70NTfufNdkWvfc5LNIb9nHPuYeudzI+Pev8feS06+AhGo8EW3clBVa3lktCzHHJr4QcT3x\nGYZhGBESI1OniFQHMBVAMgJzyMfOuUfVNtcDeCAoswHc5ZxbEFy3FsBeAIUA8p1zXNy1DLCJzzAM\noyoQoyc+59wBETnHOZcjIokApovIBOecv7niagBnO+f2ikh/AMMB9AyuKwTQ1zm3G+VEXE98jU70\nTGEde00Ns2XFQps2NSOUaVP/nXbfaWzWDWfaPLMfN2/9/pvcYrY8OnpfxOabqz/jlIBfsrUH83p5\nzV9rpvG2Pbb8gXT6Eu6IsG0Lm7Xqt+PK+3s3cpj97nWsc1W6Q2I1blR799fLSftTMd4dzudetJTN\ntHICVxC54SrVVUB9iAN/GEB68M6tpNucypX8+1/Fptb3fCa+gds5aaTv1xNIn6Qq/S+oxoPpXJ1N\nmwNH9yb9uzU8tqfZWohFJ3EKwcz3+A/4rXWPCy0/dsdoWpey/ELSH1zwJekLrm9KWpuUr3nhZNLL\nurGp89wWbPe9fhxX/zmQzWZff3PY4Q+zibdTXf65fLozN9VN+gt3wdAk+9IVAODgfu4WciDb0/m5\n3DkkJsTQ1OmcO+yXqI7APOLUen+BoBkA/D1tBOVcPMWMvoZhGFUAEYn6X5hjJYjIXABbAXzlnJtV\n7MbA7QD8f405AF+KyCwR+WUx+5Qqcf3EZxiGYURIBE98k+esw+Q560vczjlXCKCriNQGMEZEOjnn\njqgaIiLnALgVgN+M0Ms5t0VEMgB8JSJLnXPTIn0bscAmPsMwDAMA0LdbK/Tt5lWseez18PORcy5L\nRCYB6A+AJj4R6YyAb6+/35/nnNsS/H+7iIwG0AOATXyRcsjXQTueUjO1n+L2+zls/42ntyEcyUn8\nsR0oZjsg9j49zbQJHG5++4Q5xWx5JJdO4tSHN9eyX+z6zXzsWs24Un/movDXKUF1nK7TgtMbDmTx\nlXPKtNM3+4+h5YMF7B/8w3ecQsC94o+8H3ucwSWoms3kP46f+R07anssakta+xivnXJJaPnSqV/R\nummT15K+czjHEEz7mn146ew+xKgr+Tdovq/bOwA07MDvvXpPHuuVUziFZfQkzwp233+Op3Xf1eey\nYLPW8lhydnKKy5LZHOafplICnjmrPumbxnPaSGE+pxw0asn3RN4h79vZvVEqrRvyFQ8u9e7HeLAq\nlaJWE/Z/71TpNzq9oWaGd77qtbkzfUyIkY9PRBogEI25V0RSAJwP4F9qm5YAPgFwo3PuZ9/rqQAS\nnHP7RKQmgAsAUERoWRDXE59hGIYRIbGr3NIEwDsikoBAnMhI59x4EbkDgHPODQfwMIB0AC9LwFl4\nOG2hEYDRIuIQmH/ed85NjNXAIsUmPsMwjKpA7NIZFgI4tYjXh/mWfwngiMAV59waAF3062VNXE98\n+7fvL3mjCoj+u+v1Ekyb2jT62WquRnH+sQ8pZpyqbum584rfdl4LTrX4bAZX7k99YjXpzDvZlJTe\nmm10jRqy6Wjpj9x1QJtGm57SmPTWTL6uGSleg9Y+H71F6+D4U0ljKywmjeMqHzjnO5L/+YFNmwum\n8vaLd/O16KRuguQPvCapoz9kU+ZXN7O5r/++v5G++N5hpC94gSsBvXBuR9JXjOUUg++vHEz6melj\nSH82cSXpOune55yzi03vNZLYnDho/RDSd78zjvRPqlJLrav5wl/1d77hftuQv20tunKaSK4yfW5f\n4dXESe3PZtm8bE4/OHsPd3po+xFfpz0qfSZNpZXUVHrrgszQ8oFm4RwYR4nV6gwR1xOfYRiGESFW\nsiyETXyGYRhVAZv4QpT5lRCR/iKyTERWiMgDRay/TETmi8hcEflRRHqV9RgNwzAqHQkS/b9KSlRP\nfCLSCECac+5nEWkFYL1zTrugwu2fAGAogH4ANgOYJSJjnXPLfJt97ZwbF9z+ZACjAHQ84mAADh0o\nhbI+5UC0t1fztIr7oB7Op6dZt28H6Ts7s88tfO8FYNda9vcczOFq+R17NCedobpaZ6Sy/mY3+59u\n+WpEaHkmYVglAAAgAElEQVSbKi+1ejaXU7vkHU4RqH0D+/ROyeGI7QmpfyX9p+5cuuvzlIdI63vk\n4K+8kmn5bbnDQdv6/PfsCS+9S/r517eQnjGU//7sOuJN0kt/4vf6p5bsd/vfyp2kk2pwKTi/X++U\nM1rwuV7lYyUk8djTGrAfLIXdwqhVn9f//jb24a07+wzSfUZ8QXpPJqfM+M9/1eucXqPLm3V5hcuv\nbfv1raSTHuCuF/WUT9oV8k9ncpp3PybVKIXvuD3xhYj2SlwFoJWI9EXgd+m6KPfvAWClc26dcy4f\nwAgAl/s38NWAA4A0BAqaGoZhGMdCQkL0/yop0b6z6s65bwHUdM5lI9BaIhqaAfCHq20EFy8FAIjI\nABFZCuB/AG6L8hyGYRiGxkydIaJ9nl4qIt8BWCkiSQA6A/g81oNyzo1BoP5bbwCPo5iI/fzv1oaW\nlyABnZBa1GblTpvWrFevPbbjTVjNf28UaQeOAwY07kv6p2yu+HH+TK6H8vtZXO1EpyfUa1abtE5n\n4BoewBl9WpOurcxqx9X2vh4jx3OHAm3ff6o5m/cGqDD6zu+MJP2k2v/zVDZtvjWKOwN8oZrynrvS\na1z714EdaF3D79eSfuR27khw9cucc7L3zidIZ57C17FG3RqkX/1mDensrWwu1CkDG+Z6ptUtKsR/\n+W9uIn3Ghx+TvqI9d+C4YDebE7c++DTpEcosO27aWtJv3cbdHAY+w7WV/Skum+fzZ169FndXyDyL\nj9VoODfJ1RQcLCCdv59N89X25uHQ6oCBf9n67xFzKvETXLRENfE5574QkSUABgBoCOCFKM+3CUBL\nn24efK24800TkTYiku6cO9Ll093z4XSas+eI1YZhGPFCtTbpqNYmMNGf0r4XFo6aXL4DqsRE7UF1\nzq1H9BPeYWYBaBcMjNkCYBCUn1BE2h6u7SYipwJILnLSMwzDMCJHKscTn4g0hRcbMtY5tznc9kVR\npuGBzrkCERkCYCIC/sU3nHNLVY23q0TkJgAHAeQCuLYsx2gYhlEpifOJT0ROA3AmAnEiw4MvDwhO\nhD8452YXu7OizOPinXNfAOigXvPXeHsKwFORHKt+W7/9v+KaOo/Vp6ddzA+c3oj02deyH+SbURX3\nAfmad7qHlsdtnULrumZwiSjZyPex9umF8yUVtT5nL/uXqinnfa5KWXh3sZdukaC6ABSoUlfPz+Bz\nt5rLXsB//Zr9k9pHOPNbToc4tyb7DK/bxmXJCn1ZROc/MZnH/Spvq++fj+/inJNBW+4n3XU6H+/b\ndapMWBMuE+YPwweO/BySqnvv5aWLj6N1J732AemnL+AUlBMu5+7uH3zCP1m1lV93xwpOrchRPsVz\nm59Jum4rLh3n9+vpdBjtM/6Km8MfcY+kNWKfcdZmLomXWp9zM3J2emkf++uVQjnGOPXxicgAAK0R\nCK7U7vFPgtv0FJF7AKx2zv2vpGNGfCVE5DciUq/kLQ3DMIwKhyRE/69i8L1z7jkAl4tIkf2anHMz\nnHPPA5hZ1HpNNO+sEQIJ56OC1Vcqb6yrYRhGZSNOJz7n3GFTz1MAbhCRthFsG5aITZ3OuYdE5GEE\nGgfeCmCoiIxCwE/3c/i9S4ddq3eXvFElQJvFZmzhCiMXt66ttqi4ps6EVl5Y/mXNz6N1n27g6vYF\n/bjzSepPa0lrk5pmpzLR5edy+HjXHpxCOnUSh+mf2ssLQN6ZyqHstRqzuW+fqgCiGbSHm5YO+o6z\ngJas5Qaru/ZzlZCnhyznA/r+7Fz8eCu1ir8XdVZfSjozh++fwozWpL/8ir/ONRtwmlBCIv8g1m4a\n3tx46IAXxv/QdxzEnbl0O+lNZ7J5+o2/8z3wzh9/JJ36BV/33z/fnvS2JXz8hvc9TzqjUwZpfzpD\npur+kqLSOrbk8Pq5/8fpN70//Iz0/m28fWIym7NLnQoykR0twTQ3iEg/EWnnnPuypH2KI6orESxP\ntjX47xCAegA+FpGIfHKGYRhGORHnlVuCJS/hnPsGQJ6IvCAivUSkhohcEs2xIn7iCzoObwKwA8Dr\nAP7gnMsPDmYlgD9Gc2LDMAyjDInzJz4Ar4jIFgBXAKiBQPGUmwA8h0B+eKMw+xLRRHWmA7jSObfO\n/6JzrlBELi1mH8MwDKMiEP8T37kIPHQNcs5RISYR+W00B4pm4quhJz0RedI594AeRFkRRWOICo1T\nYULXFywj/WHCCaR/dTKXqNrx+E9ldu5oabeF65hLhtdlqt7jHJn82wvakL6iPqdp+MO9i6J2M9UG\nXd0ex3doQHrLPvb5aR/Obp9PsJrq7KCPnaZ8fskrLiAt2exzn/ADd1jf+SgbTGTk+6S5kBef//s7\n+9GqltmTSI9ty76mq3L/STr97/w5JKqw/JR0Drvfu4FL5tXM4LD9jBP4Ovu7EKxVPr17Bp1E+ruN\nOaSXq47tq17l6oWJOewb3fY8fxfqtqpDes86HntSteL9bEkp/POYt5ZTI3bmcgmyx2b+QHrHRh6b\n/rmqq+6ZlHreda7TnMcdE+J/4rszaOYsimHFvF4k0VyJouplXhTNyQzDMIxyIk6jOkXkDhG5HMC3\nYTY7GGxucHskxyzxiU9Efg3gLgBtRMRfRbgWgOmRnMQwDMMwjgbn3DARaQngLhHJAfCJcy4LAESk\nNoCrEfD5jXPObQxzqBCRmDo/ADABwBMA/uR7Pbsy19DURtRts/uGlht1nxzTc4k62QeJbF7UCZMt\navL6MV8dfZWHaM8dLWtafEi6+06vl8Tbg0+kdTsOcOWU3Wx5PAJt2szalF3MlgHmq8oZjc/ndCBt\nOk8I05YlT1WB8VcnAYAe93B1/SnzJpKeNp+7EjR84hnSDRqnkZ43hKuInNfTMycWjmfz3qiZ3FVg\n/tSzSP/1fa6Wok2Vf7yITc7PTeXfkvS2bIJOTAz/ZJDna+B6IJs/46ZpfN2ebtKJdLUvuVFtP9Ux\nobfqsKHRpk1NtkpD8TeA3b02fDWop6bwdclVjYz1/VSjDude71f3kD89a29ytB3fIqCCRWlGQ7BG\n9EsikgrgGhFJQ8BiuRfASOdcVD+CJU58zrm9wYNH23TWMAzDqCCIlHHeYCkQbFQevv9TBERi6pzm\nnOstItngByEJjMPp7GnDMAyjolFBfHYVgRKvhHOud/D/Ws652r5/tWzSMwzDiBNiFNwiItVFZKaI\nzBWRhSLy12K2e0FEVorIPBHp4nu9v4gsE5EVIvJAKb3bsESTwP4OgHucc3uCuh6A/zjnbiutwZVE\nt96+Uk3D2B6vfXQ33lGX9HvDwtvvtXend1OvFPvKSAd4lGi/m34vLy/kcleqSHxMz32sXLPlT6TX\nw+t4kJTAppdEVf51Vx77c+q15s9Q+2AatK9PWpfO0mH6yco3pcvP+nXDhuxz25XM/pmhquvA1oHd\nSfdxnFv71raFpK84nX1473zKGULjhnMZsq2dvXD3jqrKf9sO7Bzd0Z0Dsseu/Ij0U+fyHTTko1Wk\ndVcB3XVgX2Z490qXM73Sbz/8ktM8er35NemHN3MXC83793IJs8HPhU/laX0ql0DL3sul4DKH/Jr0\nJeNHhJanKJ+dP90AAAoLuENHck0ua6dTYHT5NN3ZPjHZux8Tkkrh6SxGT3zOuQMico5zLkcC9tPp\nIjLBOReqJyciFwFo65w7XkROB/AqgJ7BgidDAfQDsBmB+s9jnXPLijpXaRFNHl/nw5MeADjndotI\n11IYk2EYhhFrYhjcEvS1AUB1BOYR/Sfz5QD+G9x2pojUEZFGAI4DsPJwTriIjAhuG/HEJyI1ALQL\nylXOubxw2xdFNFciwd+WSETSUQ79/AzDMIyjIIZ5fCKSICJzEajb/JVzbpbapBkCDWMPszH4WnGv\nlzx8kaRgXeiNCAS4/BfABhF5SkSqhd+biWbi+g+AH0TksJ3kGgD/iOZksWbV2uK7M2hT5bslmDZL\nonGaV7W9tE2dGv1eBnc4nfQkTEJFYdDOR0h/uovbY9Wq5pl3HpzG1Uuu6sDtHi9TlTBKCi/Xpk2N\nbh6bkcrflTSVQnDIZ8rasDCTj3WQj3WwA5sHz3mATWwPXcXVdtZlHSL90ZS1pHWT0uv2/I70oETP\nrNZ5FFd5SWvC76PB754l/fBNJ7P+jt/bjr9wFZnkIQ+T1l0uNNokvcUXtt96DDcX3rmKPzNd3aSa\nqp7y41ZOhyiJncpMu2+r6qKxjLs9TPN1/MhTZlGd9pGrqsokp7Gpc+w13Uif8xx/T7PV2Go19dJz\ntFk1JsQwuMU5VwigazCPboyIdHLOhbNTx6KN3b8RyB8/zjmXDYTy+J4O/rsn0gNF05bovyIyB8A5\nwZeuLOGNGoZhGBWFCCa+yVMWYPLUhSVudxjnXJaITALQH4B/PtgEoIVPNw++loxAQWn9eiRcCqC9\n8yVIBs//awRMpbGf+IInWQxgcTT7GIZhGBWACHx8fc/pgr7nhAIw8eg/PjxiGxFpACDfObdXRFIQ\nKGf5L7XZOAB3AxgpIj0B7HHOZYrIDgDtRKQVgC0ABiHyHHHnn/R8LxaIRBeWZ3l8hmEYVYHYmTqb\nAHgnGKGZgEDllPEicgcCc8LwoL5YRFYB2I9A8/LDk9QQABOD+74RRZODJSJyk3Puv/S2RG5AFMEx\nQGSVW3pLIK77xGDZmApD42ZeSPfA7wfQupFnjiF9rAbmZkgveaMyIqVaWskbxYhB6+/iF+qyH3pE\n7QdZN/g76Qv3c4pP7WTvOnZrPILWpVfnT+nOLpye8Lsfw5fh02H2uoSU9svtO8jV9bcv20F6ly/9\nQfsHTzmzBen7T25IOkF17/77x/y9nPvAZaRHqG4NzZpxdf5RW6aSvvSVFaHlJQfYV6R9bPpP4Svb\ntSY9eQP7vXqPfAPhOJAV3s/WsgmXktvo6zz+/Bl876bv5O+Vv5MDwGW8AOD5EYt4/zbsFz64n8eW\ns5O7PehrkdeJ007wiVd+uH47Hls15RPO2cHHbt+Sr/uwhfxzmdGGj7fB508EgFrpXqf7FOUvrEg4\n5xYCOLWI14cpPaSY/b8A0KGodSVwN4BPReQ2AHOCr3UHkIJAj76IicjU6ZxzIvI5gJNL3NgwDMOo\neMR55Rbn3CYAp4vIuQAOF/odH6ZVUbFE4+P7SUROKyJs1TAMw6joxPnEdxjn3LcI36KoRKKZ+E4H\nMFhE1iFgsz3s4+t8LAM4Flb4TF8/X35GqZ5rRO2HSvX40ZCTnxV2vd+coyvW6LSOvXPPIV2nK4dc\nV3uCndsHr+1JOoN7jqLf1pdI3zuN92+R5t1y74xh036rm/hWurAVpwQkp3EQ8UHVSFaHm6fWTyWt\nK2d8PoPNizqFoMBnCm19Cr/R9BQ2e81XpkqNrsz/jx/ZZCeqE8Tlx/PntiOPzbJ7Hr82tFzrmZG0\nrqS0j3bPTiOdehqblNvVYzPbTzW4ws4hNRZtEtx9gNcn+jpXPDlJNeD9mRu8aDNtzQz+DPdvZ/Pi\nnnWqYpMylepQiMYncwWdHzM5naFmQ+8e0g130xqxmVZ/ZmnV+ef0jCZ8Hd/+lE2b+r36TZ87XPGp\nWkdNnE98IvJImNXOOff3MOuJaCa+C6PY1jAMw6hAuFhk0pUvRdXGSwVwO4D6AEpl4rvLOUcFRUXk\nSQDlUmTUMAzDiJxAznn84pz7z+FlEamFQN7ebQBGIFBgJWKiefY9v4jXLormZIZhGEb5UOgKo/5X\n0RCRdBF5HMACBB7cTnXOPeCc2xbNcSLJ4/s1gLsAtBGRBYdfBpAG4PtidywDMk7w/C7N/s1V3mfr\njSsRG7s8FXa936KhfXra2qF9enp91nN/5hfmcRmxnje2JN3wWf7DS/twbrnS68B+5tmtaZ0Oq2+a\nxvdyrx6cQjDp29WkMxfx9rWbclh9L9Wtu2Udro7/4ThOOfCXy9qwmI+96ClO42iw8N+kD+VxSbJa\nqhzayl3sq9q7nv1Jj5zO6TnVEnisGf/2zqf9WpoOp3EKilvAPuLWqkvAwOPZD/b+11yOTXcO2LmK\nP2NdOG71Kq/0V4+Lw/cS0WW/tL9Sd7qvpVInSurYcSCby5DpBIe8PV55Ne0z1rkQTVQqxbzVfB3W\n7OF0muq1+L0VHuKJxT9Wfe/GAoeKN5FFg4j8G8CVAIYDONk5t6+EXYolElPnBwAmAHgCgL/HTLZz\nblfRuxiGYRgViYr4BBclvwdwAMBDAB70tQ2LuphKJAnsewHsFZFbEZhtWx/eT0TgnHssqqEbhmEY\nZU68P/E552IWlhpNcMsYAHsRyJjX9oJywV9pI+XJ3/PKf1Wc9INouWAwh85PfJ8riqz6OfJjlRTI\nNXAUp4GMupZNmfLIK6TzVVi+KDObrpSRosxoO3M8E+ASFYp+ent+3ztyOSx+rqrk316Z8DYqs2od\nFX4++6fNpHd35GorGn86gzbBNf33c6TX/Pn/SLd58k3S2hy4fTm/l2ZdOHVja85aPt4D75L2Nzkt\nUBVodNT6illcAzjlwhNJj1Fm3KEduM1mSZVaNGkNVaPaxy4PLe94ZCytq92MTXraRKyp24pTALK3\nZBezZQBdXaVBazZPzt/Bn0N+jpcio+/lNp0bk96jqsTs38FBh9Vr83elTguuxuM/F8DdRbJSw7+v\no6ESPPHFjGgmvubOuf6lNhLDMAyj1Ij3qE4RuRyBeeiloJ4JICO4+gHn3EfF7qyI5tHxexGxkmWG\nYRhGefBHBLo+HKY6gNMA9AVwZzQHiuaJrzeAW0VkNQKmznKv3GIYhmFERmGc+/gAJDvn/KV/pjnn\ndgLYKSI1i9upKKKZ+Cpczl6DDj6fUHKN4jeMM+pddQK/8P60ojcMooPZ73++fWj5P/esoHWytB/p\nUR3D13cd89Qa0nV+vpT0tcPnkO55Koer65QDf4i3Dj2fpHwqy5XPT3dbWKFSAHTo+j6VMqDD/jdv\nZT/KKWdwuoSf1Wu4hFTvTuwfHDxhPGndpbywIHzKQUo1Nr60zObtM5/5DenzP/K6j2zI5Khu3dEg\nowNfR3fZINIFLzxP+tCYKbz/Cby/7mKh6ag+t96vFl9WMbEapydkbeJUi4adMkj7IvkAAGed3pz0\nrOU8tr0b+HjrVdf08en8u+Fv99ZAXbc9KhVi3zb26TU9nu+/81Q3huPr8s/tIx9F1UnnmIl3UycA\nctCq7g8ZiIJoTJ3rAZwF4Gbn3DoEfm8bhd/FMAzDqAhUggT2mSLyS/1isA/gj0VsXyzRPPG9DKAQ\nwLkAHgOQDeATBGyshmEYRgUm3tMZAPwOwBgRuR7AT8HXuiHg6xtQ7F5FEFV3BufcqSIyFwCcc7tF\npFy7JTbzVd//OW9VmC3ji5FXhjdtXvksxxhtv42ryT1d55nQsjaw/W36VtYljEVXt796GNfESUhk\no0HuIQ6t181hl/hC63XHgjO6cXrCmc143+9VpRVtRvWHg0dC/bZ8vIPKHLnM1/0j/ylOl6n79Kuk\nm6kw+xpqrCWlBKxSKQd9a04n/eK5J5G+/5YFoeUb7+AQ/UYnsRl2m6pok3CIUwZ2KPNg0q19Sed/\nM5e07iqQlMzmSt3g118tJbkmd7XQIf8HsticqEP+66tqKmMuu450y+deRjSo4ilIqed16NivTJnJ\nNfnnrqG6DtnZ/BmPns/ftRu7c8qK01+uUqYCPsFFRbAs2ZmqH9/nwTZFURGNqTNfRBIR/C0VkQwg\n/v+EMAzDqAo4FET9ryIhIn8EQv34tjrnXjw86YnIP6M5VjQT3wsARgNoKCL/ADANQFQnMwzDMMqH\nSuDj80dlqSLCiCrHPGJTp3PufRGZA6AfAqkMA5xzS0vYzTAMw6gAVIKoTilmuSgdlmh8fHDOLQNQ\ntjG4YZj3/frQ8kBlLr+/jMcSDc1Ugfruq58kPaYGtzi8+28qzP53C8PqcHfAZ68vJ627WGivwz3/\nakP6nSvb8dj+x362U1UXghlT15H2+/z2ZbIPpbEqC/b5z5yusEmXQ6vHfrTc3XmIhuytnAawZVtR\nfS4DvLiK0z5S0rlb+45duhI/+64ad+YA6F2qkr/uJq8/hyW72Ac45DGvK0YD5XvKVWM5qSeH/L++\nYhzpZNU14McWPPaCfDZ5ZW8JXxQ/pwF3Tb/yHO8emqTe9ymqrNwcVbJMd0FPVb7TPVdxqeB97Xl9\nSffI9hz2yyXVKP4nsbq6P7eqEnmp6fy+W6rOEXWqs4FNb+/3C6fW53WxoBLk8blilovSYYnY1Cki\n74hIXZ+uJyJvhtvHMAzDqBg4Vxj1vwrGKSKSJSLZADqLSPbhfwCiqioWzRNfZ+dcKOM4GNXZNdwO\nhmEYRsWgAvrsosI5l1jyVpERzcSXICL1nHO7gUAn3Cj3jzn+yu7rM2Nfzby02MRNAtDlXq6c0VUV\ngXvp0Q2kB+UPJe1mTiQ9srdnytLP/526ccuqb0ady2Nr+THpb//vbNJ9XptM+uA+NhV9tjJ8i0a/\neVN3V/hq6XbSWZu56obukFCSadPfSBYAzujJJuNVyiSoK3H4uwy8Np9TJVo34LD6RYsySetQdR22\nX5jPP0LaFLp4LVdfSTuNQ+GR4Bm0dRpH3ZbcBWDLTn6f321k7Q/hB4CFO9V7URVv9HVtdRKPfdm1\n15Nu8OLw0LJOT2h5Ahfc4DpAR3Zr2KU6erz0QBfSB99TboD9fD59z21X6RN+06o2V+/JYrNr3l7e\nV1+n1Qd47F8msRNCm5CzNnm/YTn12axvACIyLtx659xlkR4rmonrPwB+EJHDFbCvgUV1GoZhxAWx\nSmAXkeYA/otA5a5CAK85515Q29wPYDACf3tXA9ARQAPn3B4RWYtAi7tCAPnOuR4RnvoMABsRaI4+\nE1EGtPiJJqrzvyIyG4HKLQBwpXNuydGe2DAMwyg7YmjqPATgPufcPBFJAzBHRCYGgx8BAM65pwE8\nDQAicimAe32uskIAfQ9bD6OgMYDzAVwH4HoAnwP40Dm3ONo3EFVwC4DNzrmhzrmhALYcTXCLiPQX\nkWUiskJEHihi/fUiMj/4b5q1QjIMwzh2HAqj/lfkcZzb6pybF1zeB2ApgGZFbhzgOgAf+rQguhzy\nw+ctcM594Zy7GUBPAKsATBaRISXsegRlGtwiIgkAhiKQC7gZwCwRGev/SwHAagBnO+f2ikh/AK8h\n8CaPwG8T792ndTRDKVW6n84+2Nkzw1dAOLiX/WStrzye9LwFK0mPSObPeeDEC0gf39Zb7vy7U2jd\noAF9SI9sQRYKnHsNV5RvWOdU0hPfepu3v55D5dfP3UI6tT77j3J8/iZd7X7/dpXeoKr8n9aUw8M/\n+Zw7T2hqqXDylcqnd0j5YHR9tlkTvPdy6+94LGuU36zBcVw2bNdG9gflK1+T9lduXcB+tdrqvd7z\n7VrSN57p+Svf/JpTSvaorhWauqdyJ/G8vewrnbyBtU47qZbKPxtrF3Bprqnn/0T6LF8ni2nL2I+b\nVo2tVdqnl6pSI6qnsd8tX/nVwt1vALBVdX/QHT8yfPdcnro/JYHHqkvD6bHrtJJV6joWHCzbyiil\nEdwiIq0BdEHA9FjU+hQEksvv9r3sAHwpIg7AcOfca1GcrzqASxCYTFvDK6wSFWUd3NIDwMpgdweI\nyAgAl8OXG+icm+HbfgbC/yVhGIZhRECs0xOCZs6PAdwTfPIril8g0DfPH5XUyzm3JVj28isRWeqc\nC1+gOHC+/wI4CcB4AI865xYd7dijDW6ZISKjgvoaAP+I8nzNAPhDFDciMBkWx+0AJkR5DsMwDEMR\nSQL7zO9WYea0n0vcTkSSEJj03nXOjQ2z6SCwmRPOuS3B/7eLyGgE5oASJz4ANwDYD+AeAL/19WY8\n3BS9dnE7aipscIuInAPgVgQ6vxdJ9RWe2WRRoiBZNYIsL0oybQ786VrSI7uNIr15Flfin3gSf55L\nt6lQ5ws4neGCwV6I+KwBHWhd5/ps5tIF4nudz+vHn3kv6ZPuPJH0wc1sutKh7rVUZY7jfaakxT9x\nXoeumtFBmbkGtOMUgslt2LxYX1W7WDmbq52cflYr0jOmsonuuO5sXGj5wq2h5YwJ/Bmt3s6fQX4B\n/6joNA9JVCa9XDaLte3O5Xx+ns3XZp9qNruzmXdPFKoWA6pXKxqo78VG1SlCd9iYMI/N1bqLhg7j\n14xdzWbbL79bG1qu3ZTv5R838/vSZlp/yhIA5KjmwtXVddWpGdrUqe+xApVWstnXUUGbm1uqNJFt\n6jro7g3Hq/335vHvwsYslY6zKQsIpvAs3RDJPBAdkTzx9ejdBj16e5V2hj45sbhN3wSwxDn3fHEb\niEgdAH0QiO48/FoqgATn3L5g1/QLADwa2fhd1H7B4oh44pPA9HoqgHTn3GMi0lJEejjnomkAuAlA\nS59uHnxNn6szgOEA+oeL/KnRz3Nm6ZvOMAwjrmhWO/APQMcuvbHs8+9jevjCGLVBEpFeCExmC4Nt\n6hyAvwBohcCT1+HEzQEAvnTO+f/6aARgdNC/lwTgfedcsbNraVHWjWhnAWgnIq0AbEHgMZgaaolI\ny+Bxb3TOlfzMbRiGYZRIrIJbnHPTAZRYRcU59w6Ad9RraxAIhilXyrQRrXOuIBh6OhGBcNY3nHNL\ng63jD/+l8DCAdAAvB58yo0lwNAzDMIogVk98lYFoJr6YNKJ1zn0BoIN6bZhv+ZcAfhnJsfy+icTk\nmJVxiznJKy4kLa2U29Kx/6jJaZNJnzGdy4bdfDUXdkpWf358+Z7nd/tTbw58uuVk9r9cx9kOGHoe\nd1+Y8yt247a6lr88OoS7UHUx36Oq62cu8camfSjdVPmqlrW5sv4VbS8mfcuGJ0jvXB0+HzbrAN+u\nV13SnvR3av8bvvwotFygwuZrqU7iylWE41QZL+1vrJbK+6ckhXdf6Oua5rvfdSV/rROq8XdjvirN\nprvFt1AdNuYrv1v12nzDOTW2j5dxCbU0n59Xlyybu5Dvx7qt2I+WnMrnataM1782m/20ieq96hJl\nK3VOuWEAACAASURBVFSne+0PpdJyat3KxdzJXrtXdFm6+aq82mlduezcGnUt/GkiidVj/3tWGF0D\ng0pNNBOfbkR7NYCHSmVUhmEYRkyxJz4Pa0RrGIZhVCniuhGt36zRQpnFSkL/7fPgq51I//PO2GVq\nHOzwJekDD3OljSsP/of0p8m/J90lg01X87/n6iuPdPkF6a8zvwstf3XtZ7Ru1FCuNrF4FFd9ufap\n70jjDk4ZaJfO9p8GNbiiybYlnN6gTaEde3iVXpb+uJHWTVEh+3m9Of0g+f0RCEcLZUraoKrIJAh/\n6mO/5tip9DZctcZfqGPilDW8TpkmazVms222Sm+or1IKdEeF4+pyGL7OzNVh/WN8JkJt5s9S1UkO\nZHP6gk5vSFIpKJt2cMpAg/bhx66P17EBv5fv13smv3ot2FTZVjX0Xa+qyGSrlIHduWwebK72X6W6\nWqzbGL6KTc2GnCLjN1/q7gm62o5OSdGmzxbquizZGL7biL8SVcGB2Fd1ife2RLGkXNsKGYZhGGWD\nmTo9bOIzDMOoAtjE52ETn2EYRhXATJ0ecT3x7Vjp+Rq6qjBo/bfNjXfUJf3eMA41nnszd40efycH\nrL7wrpdzec+N86IbqBrMp3/nbguPtnqL9HcD2TcwLo/9HMefOYX0CGGdPc/rqt7vn/1oXfOa7Je4\nfFb4FICaDdm/mKn8FF/8qhfpvk9/Q7paDQ7bL/D91ZmYzH6ygoP8xTxOhdmLulubqpSBnSp8/JQz\nuOP6etXRXftRNjZmH183X3ms7qdzF4pVm7NJZ6TzddqVw361bLV99Vrs35mqOrgfgbqHOvlKln0+\ngGpAoONb3C1s53q+Lgf389h0BwSdciAq5l/7G/XxMlUniuZtveuauZWvwzT1vhurz7RQ+Uq1z3id\nGqtOd9iWwv7KlHT+zHW5Nn+pOe2zq92A/YG56n1rX+tePTZ1ndes5/1LG3vi84jric8wDMOIDMvj\n87CJzzAMowpgpk6PSjPxNavFVRNU0QW8q0ybmgkp4XPxO/lMfucNYlPk1yN26s0J/XfWjjl9Sd/A\nFhGk3/wr0p99/h7pH4dyx4XburB5qEdNz0T4aSs2Dd3+zFzSff81mHStVVwJ4/sh55Hu9sTnpDs9\nNpV0Pvh8NVSaSbav+WaGqtSim7Fu0h0OTuGuFTn/4wr2+Xl8Ief/sIF0rSZckSShGpu5xndjU+e2\nCcWbgbVp8+fFPPZGqomuNoPpKh8abU7U1VE2+ZqkLtvFpnedxqFJVlVn6qgOGgWqg8H2ZTvCHk+b\nlFcu5+3zfSkI2qyqK9Ls36m6XihzYRPVkWO7MiGf05pNnSNUc2Nd2SVvD6dP+NMbDqo0kL0q3aZx\nM+40sUOlgdRQn7FOE9E/DPXbefefTo+JBWbq9Kg0E59hGIZRPDbxedjEZxiGUQWwic/DJj7DMIwq\ngPn4POJ64qvX2ktR+ET5c7Z+3J30jVfPjurYrVuy/scri0PLX4/hMkj676jfPtGa9OiX1pK+7X0u\nh3aVr6EuACQP+WvYsT16S2fSD3KVMMy9bXxoOffba2hd01O4w/o/Z7F/SHdMP3/kD6TTMjik+7SO\n7LPBbO54rX06W3wdrnVl/HTlv5k5hyvpzzu7I+mdq3bx2FT5qbotOYVF+6oatGOf3sPTuYRaks8f\n1Ej5xdaosmC1lb9HNQY/4jrs28r+omrq+Nr3pDsoZPn8nzmHwndE19Soy8c6cID9svrcqfW5LJj2\nu40e0Jd0+0c+IO2/Nk51udDdE3S6gyRobz3Crm9Sk314u1UJM/3Qk6pKnvnLkOlSbrp7+/oVfD/p\nrhj7lH+xgyqJt1q91+wtns6tz98jI7bE9cRnGIZhRIaZOj1s4jMMw6gC2MTnEdcTX1o9z0yxbRWn\nFGxXVdyjZd16pX0NVQflD6V1+a+/TvreHDZjzR7bl3Td79iEt0iZRJ6/g7vDPvI5dwb4eBmb+Nrd\nupD0wBlXhJbzGp5O63L2ssl3zhY2tyQq09FuVQ3lCDOZCkdv2IlTFDJUBXp/0L/+HmoTWpoyu/Z7\nnTtHNDqJO01kLuJGoftUw9XjujUlPfFq7nLR/SVO1ajtM+uuV8dq2YTDzXUIf+4uZfJVZjDdFUCb\nAGup42szW1KSZ9IbsYLTF5JUE1NtJq2lKpDUVKkW+nPIVmbZFHUP1Ejk99bqZE6v2ei7NnVbsfl5\nw0r+3mqTsE4/2KHu1yyVztAxncemr/O+TP4cc9Tn5O/IkNaQ7z9dNSapOn8mhYfYh5aizKgrlNlV\nd3846Kt4c0iZn2NBYfR9wystcT3xGYZhGJFhT3weNvEZhmFUAWzi87CJzzAMowpgE59HXE98BT6/\nSMFBtpenKb9F++N53/7/14b0V7dwl4Eljd8lPez9rqHlM87/A61rNfFJ0mN3fE/6/KFc1it3N/sV\nnvpNV9K3vMBlxeq25DJMJyvfQ98B7DdpPXd5aPn1RmzX372WfXbat9RIlUo6W4WbfzttHY9N+a7W\nDL6WdNP3RyFStC9J+8VevKY16euenRPxsQFgzZzNpL87cznpaqmcUrDbl7Kgy5stW8S+pTqqs7gO\ns9+vfITaZ1dX+baydnNKQZrafvdO7x5651v2AR9SXSd0ebQsVXrrgCr310Hdbz9u4PQd7X98Y/EM\nPr4qNefveHBQrWt7Apd226W6mu9R527Ymu/1FFU6bmceHz9HlQnTPugD2ZwK4u/AkKxyUvZpH59K\nd9BpIBrtH9ffPdpW+WVjgeXxecT1xGcYhmFEhj3xeSSUvIlhGIYR7xQ6F/W/ohCR5iLyrYgsFpGF\nIvLbIrbpIyJ7ROSn4L+HfOv6i8gyEVkhIg+U4lsulrh+4ivwhQ+37swVSe4/jUPXL72bTVPP38vN\nYLf9ezXpa/Nf4HOtGBNa/mEymyp/qM6fe73Vl5FOqcdhzVmPPUh6RLW7SSffx2H69VVI9jhVpWbw\nJDZfbm7q/T2Td15rWqcrt+xczakRa1R1/D2tuZqKDuGuoRp5NhjGnST2rGdTlb86y67VHN7dsQc3\ne136I1dSSbywFWld9WOzariqQ9drKnPimFVs8suoy5/TKp+Zrf7x3JFDf6a6EotmvzK51WnBlV62\nKRP0EWbXnXzPpfjMl00b8f2xSJ1rnzJt6jQRfa7cfDaJ6Q4KeXvZpDd+NY9dh/2n+xrR1lJm140q\nHUGnBGjz4G5lMk5X7+Xc5lzdp0ZdNmdvW6LKHCmyfI2W01WllUSVJqIb9Kar74omL4vNqrrJrt+8\nqTtoxIIYPvEdAnCfc26eiKQBmCMiE51zy9R2U51z9GMoIgkAhgLoB2AzgFkiMraIfUsVe+IzDMOo\nAhS46P8VhXNuq3NuXnB5H4ClAJoVsWlR9eZ6AFjpnFvnnMsHMALA5TF5g1FgE59hGIZxVIhIawBd\nAMwsYnVPEZkrIp+LSKfga80A+E1WG1H0pFmqxLWp0zAMw4iMWAe3BM2cHwO4J/jk52cOgFbOuRwR\nuQjAGADtYzqAYyCuJz7n+yBXzeIyYF8cx2HPQx9gnx4HPQM7VBP1ZkNfJF09xbO5L/ktd5y+sCf7\nf9Z8xn6Fr246k/SWXA4/v+tvfLzk4mwMQVJUtfyBd/MfTH++2Ov2cN3rnBrRSHVAOKC6TLc+tQnp\n7L3sl6heh30025RPsO+J7J+csId9U/7w8YQktoRon54uvfW36dwdXqN9etqfeVCV4nIl/BD4/UvZ\nJfii6qlSXLoUVoP2fI/osH7doV0biXR3eX9KzCnKN7qIdz0i1UJ3c6+uwvY3q89Ufw6axSvZT6zL\njvn9ZtVUWoYO6ddpI3q9Ts3Yp67j2SMmkS6pu0P1WlxSz19STe+bvz98GcQsVU6tukoT0T7CpAS+\nFn7/Y1YqHysWlPCzAgBY9uMmLJ+1ucTtRCQJgUnvXefcWL3ePxE65yaIyMsikg5gEwB/75vmwdfK\nlLie+AzDMIzICJM2GKL9ac0oaGzcy8W2c3sTwBLn3PNFrRSRRs65zOByDwDinNslIrMAtBORVgC2\nABgE4Loo3kZMsInPMAyjClAQI1OniPQCMBjAQhGZi0BL0r8AaAXAOeeGA7haRH4NIB9ALoCBCKws\nEJEhACYiEGPyhnNuaUwGFgVxPfH5TQP1lQlv1hY2ew1ePYj0mmYjSA9ccD3p307gahRNfZX6537K\n6QQ/VWdzyy8v5jIxF436kXTmYu4ioBuo6jDnn2eHNz1oY84jp18dWn5rwSu0Lu9Q+OoN21X6ga4w\nn72eTTBNT+JK/NNUc9jTu7EZduokz8yrzUxtOrKZVJs+t+3kz1Q3+tRsW8ah67prwHdL+HNopZrJ\nJqd649NVOrSpcnsJTUm1eTGcORA40sym00jq+irszFFNcXV1Em2mPUl9ZrNUNR5tItaVYLRZ9hK1\n/Q9qPK6pN9Z8leqgr2Na4/AdEQ7u5+11RZzzz2pNesKC8OZxff4jukqHObf+bujrUqjSQvR70ab2\n5DTPnaIr+8SCSJ74IsE5Nx1AWPu3c+4lAC8Vs+4LAB1iM5qjI64nPsMwDCMyIvHxVRVs4jMMw6gC\nxOqJrzJgE59hGEYVIFY+vspApZn4CtWfM9v3cxh+Xnf26Z2b/RfSLd8bTbpxBvvdHrxrXmi55vIb\naN0p33IA+dnN2H9zfF1OEbhH+fh0Zwlt39flrfZuYB/K4MtPIF19yF9RHLrMV14rDnXXY8lXnewL\nlN9iv6pIr31hfp+ept5x7JfVPr1quvq9Kvn0+d3cUaPfP78lrUPhE1WJKX2dV6/hEmpNfGkAO1QZ\nMF3mS/t39LF1V3M9ttrN+TPWofP6uub6fFP6WFqX1Am8RVe+PzU1M/h+1u9Fd1RYtzBT7e99l3RH\nde3T090StH+ynioLpsuKdW3E73WGuse0/1z7+Pb4yt7pdASdSqFTWnR3htqqi0W2Kh2XXJN93B1P\n9nylzVtyekwsKLTmDCEqzcRnGIZhFI/5+Dxs4jMMw6gCmI/PI64nPr8JRZsNbu/CId172TqI4x5/\nk7Tef68KPe6+8Y+h5QOXceeG6d88TbpwxDukM35m04+mmqrEnrWJzTvpylxz1SVc+eeTKWtJ+8PZ\n96pGntpkl6o6EuSoSiuaFFUtv70yw65RXQSe+dUppO8bPj+0vHUBXxdtcts4bwtp3ZzzlflsGk2p\nx+t19fxsZcZNUZ/5EWYvn9YmMp1uoLsIaNOn3l6vz1WVXrR5skClFPgrnOhmq9rUWVOZ7RuoVIts\ndeyCQ6x1WL5uHrt0B6cU6FQNbb70k6aaKuu0jQYd+FwlmXVPbsCf6aED/D3W3y19P/vvsQbKDLtf\n3T86PUH/huh5pqZKW9IpLMt917GD8PfWiC1xPfEZhmEYkWHBLR428RmGYVQBzNTpYROfYRhGFcCC\nWzzieuLb7ytZ1V2VozqtEfvFFk3tS/rVFK4Cf8dbi0m3VeWz6jzxcWh5zGQOW542ZSTp9Nbs41h7\n1a9Id3z9XdKZizi9QVOQz76FrcoXpf1HzX3dwXXH6YbtuEtAngpNr5aiunEr31Pr9uxz+Vl1RPB3\nDQCAB0dxpwo/9dtxKPpuVeqqTbempHeocPCv5nM5qtzd/Llo9Nh0J3Lto/H70XQou0430H4w7S+s\nrfxqtVWl/i0qJSBB+ScLVQcFf0qBTnVITOJjV1epFTvVZ6r9jbWUz3m38pMVOuXfVOfT5dn81zVN\n+RtL4qBKp9H3Y62mnDJwq/oea/+5TpHR6/1pKHvU/aRTWNKVz3mnur/0PaD9xHrsO5Z7Ze+0TzcW\nxLotUTwT1xOfYRiGERn2xOdhE59hGEYVwHx8HnE98fmr+3dQpqQkYfPL8PncaXbBzXeQHpK6gvRG\nlQaQu55Dxv2MVpVYjlfmnO82fEp60W3cfqrXhx+T3qG6DrRWXQP2KtOTbnq6yVdB4qSe3KQ0U5lv\ndPUJ3eohby9vv2Urh4M3aczmmjVqfW01dr/5cKfq5KBJVqalfVvZ1PnWXZwqcdPz3HRXo02bnVWo\n/ILl3GHBT0JSQrHrACA5jc2k2vy8T1X12K/WazPYHnX/NW/LZuHtPlOnTvM4oqKI6iqg0wv09rm6\ny4Bav3kLm6R1JRed2uFfn6xMtlklmKd1B4NuyqWxRKUE6DQQja4co02K1Wp495z+zKspnaXMsNrE\nq82oiYm8f576HPwm68Tk8Pfb0WBPfB5xPfEZhmEYkaHLOlZlYv9nRQmISH8RWSYiK0TkgSLWdxCR\n70UkT0TuK+vxGYZhVEYKXPT/Kitl+sQnIgkAhgLoB2AzgFkiMtY5t8y32U4AvwEwoCzHZhiGUZmx\nBz6PsjZ19gCw0jm3DgBEZASAywGEJj7n3A4AO0Tk0pIO5i9vVJvN6bhuLHdM0CkByXc/HPbY9ww6\niXTzLC8sf8pm9tEdN3kp6Rk/rCetS0Y9+P03pDeprgCi/CDzZnDHd+1PaqlKffXydRX4cNwyWlev\nNVd9174o7fPQHRQOqA4Jeco3la6q5+eoEG9/Z/L924v3mwJHVrPX1fL/+C13pq9Rp7rS7Gs6Ub2X\nxeq6a9+UvyRVoko/0CkE2h+kuwo0Vte9QP0K6b+uU9TYtyu/ry7t5Ud/Rnps2mdXvxWPbY/y0+pu\n8g2VrzRXHe+AKuWVlO5LEVDl1fRYdDqNfp8r1P66C8ZNnblU4ddr2Qe4TvtW1fn3bfPuuVqqu8Iu\n5XfVnSX0PaLZv5PHru+3Os29721q/ejSPiKhMj/BRUtZmzqbAfD/im8MvmYYhmEYZYIFtxiGYVQB\nLIHdo6wnvk0AWvp08+BrR0W1eZ6pa7wAdTp51VYaq0obSxdwlY9GJ3Fllt3K7PXeLDajHfRV6uje\nkBugzpj6Ben6KvR81MAupAd9NJ/0/PuvID1gzBTSi1V1lH1DriadNpTTIZp19Kqz1FbVIVoq883y\nJZyKUVeZvXTl/yOqm+j0B2Vmq68q0u/2pV6o4iQ4vjs//K+YxbfGGT1bkJ6hGtdq09Q+ZSqtcwKb\nwbqqKjZLtvH2flOnft8pyqyVoyqv6A4F+5X574Cq6qFTAmqq/RNUKLzfRK2rvOhKLLqDgR5b7v7w\nY6lVgmkzV72XdGXC2/izl7aizar7VbrMwRw+VvbP/Jno762umPPfBVypSDdW1pVftNvAX9lIv6+M\nlty0eft6Nn2m1ONUikLHYyupI0daVh5ygyk1i+qG75JyNJip06OsJ75ZANqJSCsAWwAMAnBdmO0l\nzDrUu6xjaLmO8v8YhmHEEyknZCAl+MfZaa1Ox+x3vylhj+iw7gweZTrxOecKRGQIgIkI+BffcM4t\nFZE7AqvdcBFpBGA2gFoACkXkHgCdnHP7ij+yYRiGEQ4Vh1alKXMfn3PuCwAd1GvDfMuZAFro/Yqi\nps908LMya+1VBY/P69WK9PiJq0jXbcVmjKzNvL8/AqvNm6+E3Vc3v7z6g59I60i1D5dzBOpmZdps\neCKbd/IbtyGtI+HyfDYNXalikzJdavNgSeSrahWZylSlizlvV9GNr1x7fGh58HNsXtamTX9DXeDI\nKh06WrafMmV+piJKv1HRttoEqM2Z/qLEOvpPV5XJUbaJHBWxqiNMNToi8JCKrk1V1zXHb1JW59bv\nSxdX1hVwDipTZ50WfD/vV/eM1nXV9joCtXZzr3pPvvr11c1c66pKQKnKfNhMXceNylxYTxdZr8Ym\n4u3LeXt9rfzfc20m3aWumy4yXVOZMn+eze4SXZRd/06smes1Xs4s5O9GLLAnPg8LbjEMw6gC6BSa\nqoxNfIZhGFWAWD3xiUhz4P/bO/cgvcr6jn+/2exms7fckw2BJIQEDBQMAYEKCi1ewHG81dYLarXV\nWi9TWjutjqPTjrZTUTtTmbZmoGrBSmVEEawXok5bhUpEkpCEZMk9Idlkc9nN3rLJbrK//nFO8j7P\n933ffd+Fl909+/4+MzvZ857znvNcztkn5/u74X4ACwAMA7jXzO6WY94N4Fxmrl4AHzWzTem+vQC6\n0+8Omdl1FWnYKPCFz3EcpwqooI3vDIBPmNlGkk0Ania5VjJw7QbwajPrJnkbgHsA3JDuGwZwi5lV\nXs8tk0wvfMeCEIMV4nK9VWwHTx+IXY8V1fOniG0gtIO8flmcAeTbkq2+99DIfjiaPeXvvxdnVznV\nHYcEKGs2/yDaVvvRxsM5G4tm7SiVXUJRm5262U4Tm8qgZMipb4q//0R7rm/qmq4FedX+07k7fk4W\nvbw12t5yJLYtqX1Tx0nPr5k0QtRtvltc/jXjjV5rvmQ/aVcbshy/WGxb2/fGfW8J7JtqN9M575aM\nI2rX0ow4kBeDRglJMbF31khMi2ZT6Qlss2pv1DHX+0dt9YMy7hrmcVhCUmZL6EbTgrgvarvtDEIv\nFkjR5fYtHdG2ZrRR9Fq9h2J7t9pSQ5v1aG3v5VCpNz4zOwzgcPp7H8ltSBKRhBm4ngy+8iTiRCXE\nOOSJDsn0wuc4juOUx0th4yO5FMAqAOtGOOyDAH4cbBuAx0gagHvM7N6KN6wEvvA5juNUAeW88R3f\nfBjH5c22GKnM+RCAO4uFm5H8HQAfAHBT8PGNZnaI5DwAPyW5zcweL+uiFSLbC1+gsOyWQp8qc80T\nye64Sn5yT9Q1xK7tjXNzMsS9Dz0b7XvFTXGoxBbJ4qHFMXtEvpmzIs4gwinxfnWzf2xPLJkMiNv+\nK6/NSYCP/+/eaF/rVXEhT83S0S+u6CrHaBFTdW3X7BZ9kv3im7/Yd/73T70hDstYI3P2fODeDQBL\nV8fJuE9LKoo/fnk8jn/1dBweodKUSoI6FmF2lEHph8p5zXLuXklIfFTmSCU2TUy+bXtcFFcl5/7g\nftdiwZqZRb+rz4ZmL1E0q5FmKNG/pxeKNHog6KtKnYO98biqFK9tpWSpmTkrlko75X7skznV82lf\nju/IFazWzC36LKh5pFuTc0vYkobfnJK/WeHxWtR2rJhzZSvmXJn7+7Hjwc0FjyM5Fcmi900ze6TI\nMVchse3dFtrzzOxQ+u9Rkg8jKV4wpgvfuOqsjuM4zthwdnj0PyPwdQBbzewrhXaSXAzguwDea2a7\ngs8b0jdFkGwE8DoAWwqd46Uk2298juM4TllUMJzhRgB3ANhMcgMSvezTAJYgzcAF4LMAZgP4Vyav\n6efCFhYAeDi1700F8C0zW1uRho0CX/gcx3GqgAp6dT4BYET3cDP7EIAPFfh8DxJnmHEl0wtfS6CJ\n14mN5IDYVBaJy/YWcT9/1VWxa3yfuFVvDuwcWvB0a1ucEd7kBlNbwHxJraUFKj92+/Jo+8sPxErA\nSbEJqt3kx7tz7uvzXha7ZKsru9qa1L28X0MExAajBX77+2ObSoOEM9y+MhfC8HvLo8x1uPvXcQUN\nHeeTYkNR29MDMseLZE5bxVbaLSmjTkjqrjCEQcdF7YEnpSqF7lcbYZjGC8i3demcqg06DLfRorjH\nth+PtpulYKpNjU+mFTjUJqi2KrVHqt23Q9oe3mN6/53YfyLaXr76gmhb06mpPbNH7seuvfH5LrlW\nzifzohUSLg6KOmvKMQ1D0ud61pzYXjgwFI/DkNxveo+EY6PpzCqBZ27JkemFz3EcxykPL0uUwxc+\nx3GcKsDf+HJkeuELs2doFg510d4kkp1mSt8ick9LXezwGsqXmllFi73mFfoUOWVI2qrS0pqf7Ym2\n514aS5vbpa3qWr/ruZwrvIY6aKFZleCaxBV9prjp94s8Uy/yYqtkINm5I5bdngvkxt3dsbRJyQCi\n2SuObI0l5dqGuN8b1z0fbb/njbGUWi/hC99/Ns4Uo7JuGIYyUqgDkB8a0SqyWPvOzmhbJWaVtlRe\nnCb3c18QdqLhMto2lVH12k0ihZ7YF8uFisqVer1hcQcMj6+TcILZUrT5lLRNs8roPBxti8M+GjR0\nSO5vzWIzb6WYHYJnUzOv6N+U0yJvH5PqDVo1Re8RHfdQMtYxrQRenSFHphc+x3Ecpzz8jS+HL3yO\n4zhVgNv4cvjC5ziOUwX4G1+OTC98oauzpv/RVERqG1B38dODUg15f6zXq91ipHMr6v6tacHUtV2P\nVxuhZvLXkIIVl+VCGLZtiu1oag+skarRee7eMq6aPk0fpv2SgV7tnXuC7Pl9Q/G5OzbHOQI1lVvz\nwtgmomnlViyKbYI/Ebf+y8R+qbZWHffQ/qn9UDuZuvSHdlYg305bqkq62oN0HsJ7X21JmparVNs1\nLETtwCcl3EYrVeRVjx8obsPWVG7azwNyv7Ysip8NTftVJ+O2cGF8/MCZ+Nk4JfeAtjXsi46r+hHo\nc6cVOtR2P1KKMiC2Z+q1K4Hb+HJ4yjLHcRynqsj0G5/jOI5THv7GlyPTC1+YqX2pFI1UiaNb5EKV\nOOo1W4VIgqFsoTKqFkhVV+S65liSU3lRpUyVOToPxefXTB3qVt0WhDCUqgzR3BrLgyrPKGFRUSBf\nptW+TFdZNpDZ1jwTS5Fa8FSzdmi4gc2JH+TVrXH1hsskFOO7T+wf8Xo1Irs1Bu7rfRISoIWKtWrA\n/MtjN3ntS2+7SMLS1rxwiRGup6EOgwOxXK3SpmaB0XtEpfuG2XHbls2Lt58VSXlKQ9zWUOJTk4He\nuyqz6nOoWWa0bwel4LSO+2L5O3FUzA41QSYuzcCkz87iK+JCyr2arUeureOqWWhaLsg9S1bBcunn\neAlOmVkyvfA5juM45eFvfDl84XMcx6kC3Kszhy98juM4VYC/8eWYNAtfj6R8GhDbgdqa+kTb7+2I\n9fdGcQGfFbjOH5Gs7GqHUFf1KTWxzaNXXP7VrqG2ALUX6fnzUm0FNki1R5ZKAaVu+eoKX9soGeel\nGoO6p5+SeQhtPPXiqv6KK+Lq8E9uiLPjD58duerFs8di++JfXhNXZ3hsV5w2TO1LGgYQVppQ13W1\nz2j1BrWjqeu6VvNWm57eU32SDqtZ0uSF6P3WKja8HZKqTUNUaqZKOITYlHdJWzQ8Z+5lsR0tr3St\nUgAAEf1JREFUDNXQkAC1jWpYh9qM1T7eImntBmReZs6Nn61DkrJM79ew8rmOS+2SuG09YvdtEVt+\npzxbag9Xe/lQYJs9I/d2JRh2G995Js3C5ziO4xRH/8NczfjC5ziOUwX4wpcj0wtfGIIwNDSyNHBW\n3J5VwlM0y8LxQALJC3UQd3KVodQFu2lBnIFEr6UVFaYPxlJVs3xf5Z4F83P7Dx+JZSmVB2fIuQZl\nHI9JBhKVfxB/PW9cR8pW0SXy9FwJ4xjoir+raDjCFqkqcOTyWE5cOTc+/jciEatrfGMQYqDVEzRT\nv2YzUUlPJTzNljJ7WVylQCVkDTkIs4SoZKvytUqRKumqDKty4im5Z1RyVmlTK36EMnBfR3w/aviM\nVk9QKVQztzRJW/vPFK8MAeRLm7PkHuoK7omLRJps74jHUeXpPs0CI20rNaehrKv3TyXwhS9Hphc+\nx3Ecpzw0LrGa8YXPcRynCvA3vhyeq9NxHMepKjL9xtf+TC6Te+tVsSu8up+f6Yo1c00RpamU1BbQ\nfyyn7/cfjW12imr/qvWrjVCvpfp+XgUFSUmlVQaOBS7defKGFHY+sjNON6XjohUStG1agb23Y+Qq\n1GG6tTaxyV2zPLZzNUpqLB13dXV/x6uXRts/3B3bZBRNA6b/I+7cm7PDaUUMtZNNF5vfgY2Hou35\nkt5KbVUzZ8Xn14rtahdmbTCRsYkvL4RFz90l46Y2uculenybpOrSkJfhrvjZ0WcptOPp/aVhRv0a\n1iH2x25JSabjrqEcedUX5Fnqk76HfTsstvlOscvOEBug2oH1HtFn0cRWGlZkzwuPqQD+xpcj0wuf\n4ziOUx6+8OXwhc9xHKcK8IUvR6YXvpZFxbNXqPu4uvEfbTsabaukpzLF8uW5/c+UkDo1q4fKWnoD\nquShbvql0BCDWYG0dFwl3CmxFDR7RSyhDcg4qXu5Spd5WWZE/tGM9eH5dIzX74qlJJXFVErqFymq\nuS7u2zce3h5tr7zuomj7Iskcc1AyksxeOuv873VSpFarffRL2IaGtNSIBKdhHyckdGPajFhC1iK5\nYTiEFmvVEAB1s1d5WqX2Nsko0ixtOSEStfalrk6qRYicGKKZVlSqbJZwBxsu/swD+fJ1ntlA7jk1\nS4Qy8TGRm/VvSimThJpbdJzV7BCG/mjoQyWo1MJH8kIA9wNYAGAYwL1mdneB4+4GcDuAfgDvN7ON\n6ee3AfgnJD4mXzOzuyrSsFGQ6YXPcRzHKY8KvvGdAfAJM9tIsgnA0yTXmlnbuQNI3g7gEjNbQfJ6\nAGsA3EByCoB/BnArgHYAT5F8JPzuWOALn+M4ThVQqYXPzA4DOJz+3kdyG4BFAMLF681I3gphZutI\nziC5AMDFAHaY2T4AIPnt9Fhf+BzHcZzK8lLY+EguBbAKwDrZtQjA88H2gfSzQp9fV/GGlSDTC19o\ny1A7maLpqOatnDfift0OS3rUSYWCZkkZpmmZ1L6Yl+5K7HBqG9DwB7WradqndrHBhKhN46y4TWv6\nK3Wj18oSedcWu4imZQr7olUl9Fxa8fp0b3FbEQAsaoptLLe/dnm0/bMn9kXbaq/UkIUzwTwcaYv7\npfaapZfH4QqHpe0acpJn55XtmWLfPLI3ntPwHlKX/SEJd1F7Y4/M4VSpvqBhIvOlLXo/1tYXD1kB\nYlvYyeOxXTa0owLAkFZZkWdH7WqdPbENetHF8fm6JPSiVOhQyMwlcco7Da1Q26WOiz73Sl7FjeDv\nSE1d5UOsK525JZU5HwJwp5n1lTq8ohd/kWR64XMcx3HKo5w3vrN7u3B2X1fJ40hORbLofdPMHilw\nyEEAoTfZhelndQAWF/h8TPGFz3EcpwooZ+GbsngmpizOJTAY+sXeYod+HcBWM/tKkf2PAvgYgAdJ\n3gDghJl1kDwGYDnJJQAOAXgngHeV24dKkemFL5QtdFLrRMbKk0JlUyUQlQSPBK7G6taskpyGI2iF\nglLFN1V+0bZrVgeVPsO+qKyqUuaAyDUqN6ocs1DCPg5ui8NC8qpeyDiHIQoarqDjqHOgaIHeL/x8\nf7R9o1QNKNU3nbdwHrQtOmc718dFczXEQN3sVSbTe+RkZywJ6vnCqhla1HaqPNZ67sH++Np6j2j4\nzaHDsTSaFz4zR6oMSKWK3vbc97Xocr9In9qWkap7APlz2qXFhSXEoF7moUbkyaOBDKzPSnNrLBmX\nMkFoX/Iyu0jbFgdS51zJvlMJKhjOcCOAOwBsJrkByVP+aQBLAJiZ3WNmPyL5BpI7kYQzfADJzrMk\nPw5gLXLhDNsq0rBRkOmFz3EcxymPCnp1PgGgpozjPl7k858AuKwijXmBeJJqx3Ecp6rwNz7HcZwq\nwFOW5cj0whfq/5rqaJro73ku1qK357kmiy0htGOcEhfpIXEnV7vEPK2uLaEStWJfHJS+XCw2lG27\nY68rdZsO01tpaIXaZ5qa4+06ufZZCUfokdRaGu6gLtNqq2qYnbPr1UpYyKnueFzFFIo5Ur1Bq57/\nllQV2HMidssPq5YDQIukV9N5C/9QzBAbW6vYMg80xHYwtfdo6rcFV0o1EZnD2a3xuHaKPXJW4Lav\n/VK7qqY/02u1iM24X9qu9ky1fQ2KDbpbbLXhOGtbNdxFw4jUDqaV69Ueqdtq0+uRUA2tJjLjwpy9\n9PJL4vtt56H4WG2b2uYbJDUcZLtlbtzWZ36VC29bNrAMlcYXvhyZXvgcx3Gc8vCFL8eY2/hI3kay\njeR2kp8scszdJHeQ3Ehy1Vi30XEcZ7JhZqP+mayM6RtfOQlKiyU3LXS+UL5UF+284q0n4/2anUIl\nEqXzaC7LSM/B3hGOzJddNXO/ZtZomNuAwe3HUHdp4n7ffzTOaLJNMtjn9U1k2rCg5TSRV1QqUklX\nqZk2csZ5lcHOyYVn9nZh6tJZeWMRZqlRafGiqxdG2+1bOqLtaSLTagaSdpHseuWe0JCC0M0eyHc/\nD13l+zd3YGogL+6UOdL/Teu4z5XQCg2lULd8zViiFT/CkJaTUqVCz6WhOufmbGhXJ2ovmZ2XFeZC\nkfh2SVFdDVnpFBlXQz/CZ03lwFJFc4+JrDpb2hb2bXDHcQzK/QiTsKZ4b570H97fW9qORfs0jEj7\nkicJy349Xufpgpe3nv9dM9pUAn/jyzHWb3zXIU1QamZDAM4lKA2JkpsCOJfcdFIzuON46YMyxNm9\npbM/ZInJ1h8AOLO7s/RBGWJo5+R6hiqNDduofyYrY73wFUtcOtIxBwsc4ziO44wCX/hyZNq5ZdXC\nFbkN9WQTqWnoVCwvnp4eSygquSgWSHQ9wyPnY60XL7ppTfG2epTWz6zHrsYjuGTeJQCAgZpYhlWJ\nTotpnqqP+9IwN9cXvXl7h2J5T70ylSmSLHeoKR5HzV4xOD2RF/c1d2PJwhV5CbhDOUjb1ihz0Hph\n7KU5e24s//QjlvjmzIkza5xsiNs6UBfLhaemx9sqVYVep3ubO7EkuN90DrUvtQ3iLdskMti0eM5q\nRXKbKhJzz6AUh52X66uOg16rTu6/oXSO9jR14uIFyzHYFEvCC+bEmWBmiPdrvciR6ompmYkYyN1n\nB1Xui8/VLF7GJ6RvNXXx/TZlau5auxuPYsmCODF5rWZXEblxqEU8WIOMTypN1knb1KuYIuvXiMSs\n3rTq1RyO2uKWyotck3khGy0cSwNmmrPtb83stnT7U0hS3NwVHLMGwH+b2YPpdhuAm82sQ87ls+g4\nzqTGzCpS1YDkXiQpxUbLPjNbWok2TCTG+o3vKZROUFowuameqFI3hOM4zmRnMi5eL4YxXfiKJSgl\n+WGUSG7qOI7jOJVgTKVOx3EcxxlvJnSSapKXktxAcn36bzfJPytwXCYC3svpD8mbSZ5Ij1lP8jPj\n1d5yIPkXJLeQ3ETyWyTrChyTifkBSvcng/NzJ8nN6U/es5Mek6X5GbE/WZgfkl8j2UFyU/DZLJJr\nST5H8jGSM4p8t2QCEKcMXkg0/3j8IFmk2wFcJJ/fDuCH6e/XA3hyvNv6IvtzM4BHx7t9ZfbhAgC7\nAdSl2w8CeF9W56fM/mRpfq4AsAnANCRlZNYCWJbh+SmnPxN+fgDcBGAVgE3BZ3cB+Ov0908C+EKB\n700BsBOJk0otgI0AXjbe/cniz4R+4xNeA2CXmT0vn2c14L1Yf4DYs3miUwOgkeRUAA1IFvOQrM1P\nqf4A2ZmflQDWmdlpMzsL4BcA3ibHZGl+yukPMMHnx8weB6AZEd4M4L709/sAvKXAV8tJAOKUQZYW\nvncA+M8Cn2c14L1YfwDghlQK/SHJy8eyUaPBzNoB/COA/UjG/YSZ/UwOy8z8lNkfICPzA2ALgFel\nMloDgDcAuEiOycz8oLz+ANmZn5D5lnqvm9lhAPMLHFNOAhCnDDKx8JGsBfAmAN8Z77ZUghL9eRrA\nEjO7Gkle0++PZdtGA8mZSP7HuQSJTNhE8t3j26oXTpn9ycz8WJID9y4APwXwIwAbAJwd8UsTmDL7\nk5n5KYF7Hb6EZGLhQ2KHeNrMjhbYdxDx//ouTD+byBTtj5n1mdnJ9PcfA6glOVuPmyC8BsBuM+tM\npafvAXilHJOl+SnZn4zND8zsG2Z2rZndAuAEgO1ySJbmp2R/sjY/AR3nJGaSrQCOFDjmIIDFwfaE\nnquJTFYWvnehuCz4KID3AeczwxQMeJ9gFO1PaF8heR2SkJOJmk14PxJZqZ5JnqpbAWyTY7I0PyX7\nk7H5Acl56b+LAbwVwANySJbmp2R/MjQ/RGyLfBTA+9Pf/xDAIwW+cz4BSOpt/M70e84omfC5OlMt\n/zUA/iT4LLMB76X6A+DtJD8CYAjAABJb4ITEzH5N8iEkktMQgPUA7snq/JTTH2RoflK+m77xDAH4\nqJn1ZHV+UkbsDzIwPyQfAHALgDkk9wP4GwBfAPAdkn8EYB+AP0iPXQjgXjN7oxVJADIefcg6HsDu\nOI7jVBVZkTodx3EcpyL4wuc4juNUFb7wOY7jOFWFL3yO4zhOVeELn+M4jlNV+MLnOI7jVBW+8DmO\n4zhVhS98juM4TlXhC58z6SD5eJbOO1Gu5zjVgmducZyXCJI0f8AcZ8Lhb3zOuEGygeR/pbXTNpH8\n/fTzO0iuI7me5FeZsITkNpLfIPkcyf8geSvJx9Pta4Pz9ha41mi+/zDJp0huJvnBQucl+Yl0/yaS\ndwbXaCN5H8nNSLLnY6TzkryW5DMk60g2ktxyrobcuesVGyfHcV4Y/sbnjBsk3wbg9Wb24XS7GUlh\nzS8CeGualPdfAPwKwC8B7ACwysy2kvwNgI1m9kGSbwLwATN7a3qeHjNrkWstGcX3Z5rZCZL1SDLi\nv9rMus6dl+Q1AL4O4HokFdvXAbgDSZmcXQB+28yeKtDfYuf9HIDp6c/zZnZX2I9C42RmvXLudwKo\nQ7LYHjGzf3tBk+I4VYC/8TnjyWYAryX5DyRvSv+Y3wpgNYCnSG4A8LsAlqXH7zGzrenvzwL4eXCe\nJWVcr9zv/znJjQCeRLKQrJDz3AjgYTM7ZWb9SOr2vSrdt6/QolfivJ8H8FoA1yBZ9JVC43Qekpci\nWRjvR1KYdXeR6zuOA1/4nHHEzHYgWeQ2A/g8yc+mu+4zs9VmdrWZrTSzz6Wfnw6+PhxsDyMusUWS\nH0mlwfVMCnuW9X2SNyNZbK83s1UANgKoH0W3+gt9WOK8cwE0AWgudC0Zp78j+Rk55D0AfpD+vhrJ\nwuo4ThF84XPGjbTW2ICZPQDgywCuRvIW9nbmCo7OYlJ0FIgLd+adLtwws6+mC+dqMztcxvfPMQNA\nl5mdJvkyADcUOOaXAN7CpGBtI5KCqL8scY2RzrsGwGcAfAvxGx+BvHH6EpLFLWQmgDaStQBaCux3\nHCdgwheidSY1VwL4EslhAIMAPmJmbekbzVqSU9LPPwagA0BokFbj9Ej7RnPMTwD8KclnATyHxL4Y\nn8RsA8l/R2KnMwD3mNkzqR1xVOcl+V4Ag2b27bS/T5C8xcz+JzhX3jjJue8H8DoAV6TnboXjOEVx\n5xbHKQOScwD8xswuHu+2OI7z4nCp03FKkEqN/4dEZnQcJ+P4G5/jOI5TVfgbn+M4jlNV+MLnOI7j\nVBW+8DmO4zhVhS98juM4TlXhC5/jOI5TVfjC5ziO41QVvvA5juM4VYUvfI7jOE5V8f8Kzg+s85OV\nBQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "results2d = np.array(results).reshape(Ngrid,Ngrid)\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "fig = plt.figure(figsize=(7,5))\n", "ax = plt.subplot(111)\n", "extent = [min(par_a),max(par_a),min(par_e),max(par_e)]\n", "ax.set_xlim(extent[0],extent[1])\n", "ax.set_xlabel(\"semi-major axis $a$\")\n", "ax.set_ylim(extent[2],extent[3])\n", "ax.set_ylabel(\"eccentricity $e$\")\n", "im = ax.imshow(results2d, interpolation=\"none\", vmin=1.9, vmax=4, cmap=\"RdYlGn_r\", origin=\"lower\", aspect='auto', extent=extent)\n", "cb = plt.colorbar(im, ax=ax)\n", "cb.set_label(\"MEGNO $\\\\langle Y \\\\rangle$\")" ] }, { "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 }