{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] } ], "source": [ "import numpy\n", "import tqdm\n", "import wendy\n", "%pylab inline\n", "import matplotlib.animation as animation\n", "from IPython.display import HTML\n", "_SAVE_GIFS= False\n", "rcParams.update({'axes.labelsize': 17.,\n", " 'font.size': 12.,\n", " 'legend.fontsize': 17.,\n", " 'xtick.labelsize':15.,\n", " 'ytick.labelsize':15.,\n", " 'text.usetex': _SAVE_GIFS,\n", " 'figure.figsize': [5,5],\n", " 'xtick.major.size' : 4,\n", " 'ytick.major.size' : 4,\n", " 'xtick.minor.size' : 2,\n", " 'ytick.minor.size' : 2,\n", " 'legend.numpoints':1})\n", "numpy.random.seed(2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Phase-mixing and violent relaxation in one dimension\n", "\n", "We study gravitational collapse in one dimension, illustrating *phase mixing* and *violent relaxation* using the setup of [Binney (2004)](http://adsabs.harvard.edu/abs/2004MNRAS.350..939B) and [Schulz et al. (2013)](http://adsabs.harvard.edu/abs/2013MNRAS.431...49S). We start with a grid of particles on a uniform grid in $x$ with a small velocity perturbation\n", "\n", "$x \\sim [-\\pi/2,\\pi/2]\\,,$\n", "\n", "$v = -V_0\\,\\sin(x)+\\epsilon\\,V_0\\,\\mathcal{N}(0,1)\\,,$\n", "\n", "where we also add a very small $\\epsilon \\approx10^{-7}$ additional perturbation to the velocities to avoid numerical instabilities. We integrate this system with $N = 1001$ particles for $\\approx132$ time steps:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "N= 1001\n", "dx= numpy.pi\n", "V0= 0.001\n", "\n", "x= dx*(numpy.arange(N)-N//2)/N\n", "v= -V0*numpy.sin(x)+numpy.random.normal(size=N)*V0*10.**-7.\n", "m= numpy.ones(N)/float(N)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "g= wendy.nbody(x,v,m,0.05,maxcoll=100000000)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [] } ], "source": [ "nt= 2700\n", "xt= numpy.empty((N,nt+1))\n", "vt= numpy.empty((N,nt+1))\n", "Et= numpy.empty((nt+1))\n", "xt[:,0]= x\n", "vt[:,0]= v\n", "Et[0]= wendy.energy(x,v,m)\n", "for ii in tqdm.trange(nt):\n", " tx,tv= next(g)\n", " xt[:,ii+1]= tx\n", " vt[:,ii+1]= tv\n", " Et[ii+1]= wendy.energy(tx,tv,m)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, we check that the total energy is conserved during the $N$-body integration:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVsAAAFZCAYAAAA2B4zmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGvxJREFUeJzt3VuMHHda9/HvE3s3uwQyYxtxEBfsjN+ImyzrcRyEuBo8\n9vKKvUFxEoO0HC5iJ5EQcBHvJlfrq/VkE7GXcdYrJBBCspfkDiHGh8wdIJvYQuIgrT3JC+IkZHuS\nLEvYZfO8F1Vtt8vdPT2Z7n/XzHw/Ustd1XV4uqr753/9q6onMhNJ0ng9MOkCJGk7MGwlqQDDVpIK\nMGzHICLmGsPH6sfpSdWkzan5WdLmZdiOWEQsAN9qDJ/PzDPAjYg4MbHitKk0P0va3HZOuoBJi4hF\nYCkzLzXGnwBuALPAxcy8OszyMvNiRNzoGjULzAGvAivAoZEUfm+tx+qn54A9wLHMfKnPtFPA8Xrw\nALDYeW9dr60Ch4HXM/PisOuot+V1YC/wjcx8p7Hc28A0cHWdy+25LwbVO2i+HttkAZjOzDd6vT6s\nAp8lbWaZuS0fwAJwAvg2cLDx2jlgX9fw0jqX/Rd9xp8GnhjDezkBfAT8oH4/nxkw7emu5zPArc70\nVMHb/dpHwMPDrANYamyzK931NaZdXMdy++6LNeodeh8CV4BnNtNnycfme2zbboTMvJiZrwDv9Hh5\nITOvdQ2vRMTBjawvImaq1eab65jnyJCT3gamgF2Z+UhmvjughjstpaxanivAk/WoY533Wb8GVWts\n4Drq5T7W2Ga3urbZ4UYpnVbeMLUP2heD6h1qH9at2g21Hkt/lrQ5bftuhKb6y7fSGN05TL1UH/Z2\n3wkSVCH6zTUWfTwzn19nOY8DwxzaRmZ+MMR001Stylca4/fU/z7WCbuImKV6n51tMWgd+6layN1W\n6vGXgN0RsZiZL9avHcqqD3vgctfaF/3qHWK+btNUgT9yY/wsaRPaUmFb942tZKPvrb4KYLFfi69h\nuse4m1T9m3SFxMBSGus/QhVyRMRCdvUrrmHo2/si4hmq0DgAnMse/YKZeTUiHmuM3k91CExj+xwH\nvpSZ7w+xjlVgd4+y9tb/HgMuRsQh4Czw5SFrX2tf9Kw3IgbO17XeI5n5RkQ0W96t/Sxp89pSYZuZ\nr0TEqYig8yWpvxynh/xyQO/QGFodrDMR8QLwDarW6RngZkTsBk4Bw4btsM53vb83IuJ6ROzvDsqO\n7kPaiDhez/tW17gZqm6FGeCrQ67jCvcHyyx1i7EO+bNUJwcXgbeBzrIGLXfNfdGn3mHmm2JAi7aN\nn6Ve+1Obx5YKW4DMfKnzJaE6XDvd6DNbS/NwGO4eZg+z/je499D/Ih//SzdUq6bHl38VeBroezha\nt/6OZOYvNZb1DvBKHWJvd4Jv0Doy872I+FpEHMzMS/W8q9SH0F2twefr50sR8VhmXluj9jX3Ra96\nh5kPeHqtlmULP0vaxCYWtv0uk2lM0/eSoUHqL8l11v/lgOrL3uvwr9n3NnJ1EHW6DgJYqLdBJ3QT\n+Fb3NqtD5m8yszvQV7h7CN/PIvBUY/1TmfkeVCEWEavASxHxjbXWUW/zIxHxBHdPFN2oL8q/3gnV\nzHyu3jfPRsTX1ljuwH3Rr17gwhrzzQCXB2ybOzbrZ0ntUzxs65MG+4EjVJcLDXK8PsvbmXcxIi6v\ndTgVEeeoDi2Pdvrlhq0vq2sbmy3RWarLtsYqM5/rHo6IU9nnetmGLzWGpxlwhr3uj1zsbMc6EHcD\n57n/RpdpqpBfcx3d27kOtHNULcJmuJyhCvuByx20Lzo3i/Sqt56v2YLs3of7qQ7PD1H9R3YA2FV3\nGdxzNLBZP0tqn+KXfq1xmUzToEuGeqpbh1+tD1FfAg6s4xKqjgsRsa9reGZQC3yM1uxGqA+j77Se\n6u6BmU5oRMRcdN3yWW+Lt4HbETFVH3Y/RrVtm8E3Q3XC6t1B66jH3YqIh7vWcbYO8wvA0cZyD1Hd\ngLDmcum/L1b61Vs/P99vH2bmG5n5av14pV7W+R5Bu5U+S5qwtvfZDrpk6D51i+2ew72ufrcbjZND\nc1QhsEDVqjmbma/WLx8HXqwvJ3qc6mx6m52Ju7cBz3Lvf1JHqa5jfb5ubX6Le7sqEjicme9GxNX6\nZMx7VK2/Y10nzwatA6rgO1S3KLOzLev+3FNx9+4yqM7yd/bFWsvtuS/qboNB9Q61D+t1L1C1dG9l\nfR30Nv4saUwiczI/Hh4RS1SHsoP6bPdRXRe5QnXJ0LfWcSZ401tHN4Kklmv1HWR16+EsVctskbVP\n+mw1Q53EkdR+rQ7bus/s5cx8hOqkylKj/2tLy3Xc2iup3VrbZzvokiGg522vEeEfVJM0Fpm5obv5\n2tyynaX3JUMDZQt+3afX4ytf+crEa7A+62vro+31jUKrwrZxmVLfS4bKViVJGzeJmxoGXSZz5zKl\nXPuSIUnaNIqHbVa/6HQVeLHHay82hq8BWyJc5+fnJ13CQNa3Mda3MW2vbxQmdp3tOEREbqX3I6kd\n6lu5t+wJMknaMgxbSSrAsJWkAgxbSSrAsJWkAgxbSSrAsJWkAgxbSSrAsJWkAgxbSSrAsJWkAgxb\nSSrAsJWkAgxbSSrAsJWkAgxbSSrAsJWkAgxbSSrAsJWkAgxbSSrAsJWkAgxbSSrAsJWkAgxbSSrA\nsJWkAgxbSSrAsJWkAgxbSSrAsJWkAgxbSSrAsJWkAgxbSSrAsJWkAnZOasURsQgsZealNaabAZ4E\nbgORmWdK1CdJo1Q8bCNiAdgPHAGW1ph2Bng5M5+uhy9HxOXMvDb+SiVpdIqHbWZeBC5GxOEhJn8d\nONU1vJCZ74+nMkkan9b22UbEFFW4vtUZZ9BK2qwm1mc7hFlgNSIOAruAGeBq3TKWpE2l7WELcKtz\nEi0irkTEk5n57uTKkqT1a203ArAKTDdOhq0Az06oHkn62NoctitUgdscN9tjWklqtdZ2I2TmOxEx\n3Rg9TRW4fZ08efLO8/n5eebn50dem6StbXl5meXl5ZEuMzJzpAscesURS8Bi900NETEHkJlX6+FT\nwIXOSbGIuAw81a/PNiJyUu9H0tYVEWRmbGQZk7ipYQ44CiwAuyLibGa+Wr98FJgCngfIzJci4lR9\nc8Ne4JgnxyRtRhNr2Y6DLVtJ4zCKlm2bT5BJ0pZh2EpSAYatJBVg2EpSAYatJBVg2EpSAYatJBVg\n2EpSAYatJBVg2EpSAYatJBVg2EpSAYatJBVg2EpSAYatJBVg2EpSAYatJBVg2EpSAYatJBVg2EpS\nAYatJBVg2EpSAYatJBVg2EpSAYatJBVg2EpSAYatJBVg2EpSAYatJBVg2EpSAYatJBVg2EpSAYat\nJBVg2EpSAYatJBUwsbCNiMWIOLiO6Rci4sg4a5KkcSketnVongDWG5wvA7vGUJIkjd3O0ivMzIvA\nxYg4POw8EbEA3BhfVZI0Xpulz3YauD3pIiTp42p92EbEkcx8Y9J1SNJGtDpsI2IKW7SStoBWhy3w\ndGZemnQRkrRRxU+QDSsiZoDL653v5MmTd57Pz88zPz8/uqIkbQvLy8ssLy+PdJmRmSNd4NArjlgC\nFvu1XOtramc6g8BRqisSzmfmN/vMk5N6P5K2roggM2Mjy2hVyzYi5gAy82rzpFhEPM6AoJWkNise\ntnWgHgUWgF0RcTYzX61fPgpMAc835jlRTz8TEbcy882SNUvSRk2sG2Ec7EaQNA6j6EZo+9UIkrQl\nGLaSVIBhK0kFGLaSVIBhK0kFGLaSVIBhK0kFGLaSVIBhK0kFGLaSVIBhK0kFGLaSVIBhK0kFGLaS\nVIBhK0kFGLaSVIBhK0kFGLaSVIBhK0kFGLaSVIBhK0kFGLaSVIBhK0kFGLaSVIBhK0kFGLaSVIBh\nK0kFGLaSVIBhK0kFfOywjYiDXc8XIuKJ0ZQkSVvPRlq2050nmXkR2L3xciRpa1p32EbEsYi4Arwc\nEZcj4kpEXB5DbZK0ZURmrn+miCngQN2ibY2IyI/zfiRpkIggM2NDy9hIOEXEPuAQ8HZmXtpIIaNg\n2Eoah1GE7UZOkB0BngMCeDoinlnn/IvdJ9n6TDMVESfqx9mImPu49UrSJO3cyMyZ+VzneR2+a4qI\nBWA/cARYWmPylzvriIgZ4G8iYn9mvvvxKpakydjI1Qi3G8NDHb9n5sXMfAV4Z9B0dbje6JrvHWAF\neHKddUrSxG2kZbs3ImapAnD/iOrpNg0sAq80xu8Zw7okaaw+dss2M89Q99cCq5n56siqqpZ/FXis\nMXo/a3c9SFLrbLTP9gxwZkS19Fr+tc7ziDgOnM/Mt8a1Pkkal01xu25ETANHMvOXxrUOSRqnzXK7\n7iLw1BiXL0ljte5uhIg4BjwLTEXES1T9tgm8PuLaOus7ASxm5vv18Fzdn9vTyZMn7zyfn59nfn5+\nHGVJ2sKWl5dZXl4e6TJ73kEWEX9BdZXBDWAlM99svL7h23UjYokqRC91jZuDOyfHOtfurgJX6kn2\nAvsz85t9lukdZJJGbhR3kPVr2R4GZnvdPFDflDADnPs4K6wD9SiwAOyKiLNdVzIcBaaA5+vrbL/F\n3et3Oy3owx9nvZI0Sf1atucy8+mBM1aH9weA0225QsCWraRxGGfL9mbXSuaA281Wbma+EhGHgAvA\njo0UIUlbXb+rEe4keN1/+lRE3IyI17ov8crMC0Dfk1WSpEq/lu09x+J1K/ZwZj7fY9orPcZJkrr0\na9keiIifbox7u8+0zR+kkSQ19GvZPgasRMQqVZ/seaqrBHrZUKexJG0H/Vq23wD+D9C5aeFrwLN1\nv+3ZiHghIj5XT+vpf0laQ79Lv+67S6u+keEw1Z/BOQTMUgdtZrbiaoSIyK9/Pfmt36qG//AP4Y03\n4N134cd/HH7+5+G3fxt+8ifh9dfhT/8U/vmfYWoKfuVX4Dd+A2Zn4U/+pHr83d/B7t1w4AD87u/C\no4/CH/8xnDsH16/DJz8Jv/zL8MUvwuc+B3/91/D1r8O1a/Dgg/DZz8IXvgALC/Af/wF/8AfVa++9\nVy133z74hV+An/kZ+Ld/g7/6K/j7v4fvf7+a/8d+DB55BH7qpyAC/uVfqppu3ape/6Efqh4PPwwP\nPQQffgj//d/wne9U0zzwAHz609XjwQfhR34EMuF736um+6//gk98AnbuvPvvzp3w0Ufwgx/cfXzn\nO9V8Dz5YrWvHjuq9Z8L//m/1+P73q2V2puvU9+CD1b743vfurjezqq37sWNHNV4q6Xd+p/pMr2Wi\nf4OsK3xPZeYjGyliVCIiv/jF5M/+rAqAL3wBfu3XYO9e+OAD+PM/h9deqwLk85+H3/zNKhD//d/h\n7Fn4oz+qvvif/Sw8+yzMzVXBeP48/P7vV+t49FH49V+Hn/u5KqzefLOa7xOfqELq934PDh2C//kf\nuHy5Wudbb1Xh+tRT8Iu/CNPTcPMm/O3fwl/+Jdy4AXv2VP8Z/OzPwqc+Be+/XwXwP/0T/Ou/Vuv+\niZ+o1v+jP1qF23e/WwXh6moVtJ/+dDXvQw9V68uswu27363qee+9atxDD1XT/fAP3w3K7n937Lgb\ngDt2VNM/8EC1jg8+uDcwu0P6U5+qpv/ww+rx3nvVNnnggSqcP/nJapoHHqjGdz8+/LB6XSrpS1/a\nBGHbVcTp7j+PM0mdmxr+8R+rQJjr8RfL/vM/q1bmo4/e/9qtW1UreH+Pn0L/4AP4h3+Axx+vWpnd\nPvywCs79+6vQaeoEjqTNqS1hO5WZ721oISPiHWSSxqEVYdsmhq2kcZjonzKXJA3PsJWkAgxbSSrA\nsJWkAgxbSSrAsJWkAgxbSSrAsJWkAgxbSSrAsJWkAgxbSSrAsJWkAgxbSSrAsJWkAgxbSSrAsJWk\nAgxbSSrAsJWkAgxbSSrAsJWkAgxbSSrAsJWkAgxbSSrAsJWkAiYWthGxGBEHh5juREQ8EREvRMRc\nidokadR2ll5hRCwA+4EjwNIa054DvpqZ1+rhJeDzYy9SkkaseMs2My9m5ivAO0NMvtAJ2trKMK1h\nSWqb1vbZ1i3glcboVeDwBMqRpA1pbdgC0z3G3QRmSxciSRvV5rDdPekCJGlU2hy2t3qM21O8Ckka\ngeJXI6zDKr27Epr9uPc4efLknefz8/PMz8+PtChJW9/y8jLLy8sjXWZk5kgXOPSKq8u4FjPz0oBp\nbmbmnq7hc8DpfvNERE7q/UjauiKCzIyNLKNV3QgRMde4ceFCROzrGp4ZFM6S1FaTuKlhDjgKLAC7\nIuJsZr5av3wUmAKer4ePAy9GxCzwOHCsdL2SNAoT60YYB7sRJI3DlutGkKStyrCVpAIMW0kqwLCV\npAIMW0kqwLCVpAIMW0kqwLCVpAIMW0kqwLCVpAIMW0kqwLCVpAIMW0kqwLCVpAIMW0kqwLCVpAIM\nW0kqwLCVpAIMW0kqwLCVpAIMW0kqwLCVpAIMW0kqwLCVpAIMW0kqwLCVpAIMW0kqwLCVpAIMW0kq\nwLCVpAIMW0kqwLCVpAIMW0kqwLCVpAIMW0kqYOekVhwRJ4AbwCxwMTOv9pluCjgO3AamgauZebFY\noZI0AhMJ24g4B3w1M6/Vw0vA5/tMfjwzX+madzEiLmfm+wVKlaSRmFQ3wkInaGsrEXGwz7SHG8Od\n1rAkbRrFwzYiFoCVxuhV7g/Vjt0Rsdg1fKgR1JLUepPoRpjuMe4mcKDP9M8AlyLiEHAW+PK4CpOk\ncZlEN8Lu9Uxct2LPAlPAIrB3HEVJ0jhNImxv9Ri3p9/EEXEaeDkzHwHOAEsRsW9cxUnSOEyiG2GV\n3l0JzX5cImIOuJ6Z7wJk5nMRcR14Fni+18JPnjx55/n8/Dzz8/MbLljS9rK8vMzy8vJIlxmZOdIF\nDrXSiJuZuadr+BxwOjMvNaY7AmRmvtk1bgpYzMz7wjYichLvR9LWFhFkZmxkGZO69OtCoytgphO0\nETFXt2gBLgBHG/MeAl4vUKMkjcykWrZTwIvAZeBx4GzXDQ6LwFSn5VqH8q8C1+vZV5ot4K7l2rKV\nNHKjaNlOJGzHxbCVNA6buRtBkrYVw1aSCjBsJakAw1aSCjBsJakAw1aSCjBsJakAw1aSCjBsJakA\nw1aSCjBsJakAw1aSCjBsJakAw1aSCjBsJakAw1aSCjBsJakAw1aSCjBsJakAw1aSCjBsJakAw1aS\nCjBsJakAw1aSCjBsJakAw1aSCjBsJakAw1aSCjBsJakAw1aSCjBsJakAw1aSCjBsJakAw1aSCjBs\nJamAnZNacUScAG4As8DFzLw6YNoZ4EngNhCZeaZMlZI0GhMJ24g4B3w1M6/Vw0vA5/tMOwO8nJlP\n18OXI+JyZ15J2gwm1Y2w0AjLlYg42Gfa14HXBsy7KSwvL0+6hIGsb2Osb2PaXt8oFA/biFgAVhqj\nV4HDPaadogrXtzrjMvP98VY4Hm3/MFnfxljfxrS9vlGYRDfCdI9xN4EDPcbPAqt1q3cXMANczcyL\nY6xPkkZuEmG7ex3Tztb/3srMSwARcSUinszMd0demSSNSWRm2RVGHAFezMzHu8YtAjOZebQx7QKw\nlJk7usadA25k5ks9ll32zUjaNjIzNjL/JFq2q/TuSmj243bGrfYYN9tj2g1vDEkal+InyOr+1mZX\nwixwvse073B/ME/TO5glqbUmdVPDhYjY13UJ10xXn+wcQNdNDl+LiIWuk2KPAU+VLXdrqrtppjPz\nja5xQ99sonaou+GWOt+hrvF992XJ/dyrvvpKo+P14AFgsU31NV4fzfckM4s/gCngFPBE/e++rtcW\ngdca058CnmlO2/X6iXpZLwBzk3hPjfd2on6cbdbTslqvAM90DZ9r7IulCdU1U2+nZ4Bjbdp+Xfv3\nmbqGhUnVByzU6/s2cLDxWt99WWo/r1Hf6cb+vgV8pi31NaYbyfek6Ad1TBusFQHRpg/RkHUu1P8Z\ndH+Ibjbfy6AP4ZjqmgHOdQ1f7myzNmw/4ERjeBF4eJL1AUs9wqzvviy9n5v1df4zbUxzBXihfn5r\nkvU1XhvZ92Qr/BDNeu5GG6v61uIbneGs+pxXqH7XAeBQW2ql6vu+3RlYz80mYzbojsE27Ovm9ugc\nSkI76hu4L1uyn6ep/pNq2lPXd6MxfhKfw46RfU82ddi25IPTbVN8iCLiSHb1P9X63WzS88qPcRh0\nx2CL9vXuuo+v41BmXmtRfTB4X058P2fVv/lYY/R+qpPkE6+vY9Tfk00dtrRox8Dm+BDVgXa7x0vr\nudlkXO7cMRgRRyLihTrEoCXbj6qv9nh9c80J4Mv1+LbUB4P3ZRv2M91HABFxnLsnqFpR3zi+J5s9\nbFuxY7q1/UMEPJ29z7re6jFuz7iLaei+Y/CNzHwVeDkiPkNLtl+9f89SnShbBPbWL7WivtqgfdmG\n/XxHREwDRzLz/9aj2lLfyL8nmz1s27Jj7tPGD1Hdp3y5z8vrudlkXFapLrG5p98TeJYWbD+AiDhN\n9ZOfjwBngKWI2EdL6qsN2pdt2M/dFrn3Us6J1zeu78nEfjx8RCa+YwZo3YeIqktjJiIOAUF1feOu\niCAzvxkRvW42OV2wvkF3DF5g8l/COeB61r/LkZnPRcR1qv8M/nTS9XVk5sU++/K1zHyrBfsZuHOt\n6mJXv/zcgNpL1jeW78mmDtuW7Jj7tPVD1Ozsj4jHgfOZ+c16VN+bTQrV9059RNBtGlipt1+zpVh6\nX89yf3ieodrXbaivW699+daA14rtZ7jzGylvA7fr/tG9VCF3ddL1je17Mq5r10o9qPrPuq9tvDzh\neo5QXZs3VT/2U1+j16ZaqS7mvkl1uPREPa7vzSYF6zpF140CdX0/3YbtV2+fsz32d+c64KL1AXNU\nR1A/qLfTC41a+904VGQ/96uP6jrbj+rxP+h6frAN9TWmGdn3pPivfo1a/b/ii1Qb43GqL8NE/pJD\n13W2nY0a9fPDmXmpTbW2WUScotqOe+naRm3YfnX/7K8C1+tRK3n3VvOJ16f22vRhK0mbwWa/GkGS\nNgXDVpIKMGwlqQDDVpIKMGwlqQDDVpIKMGwlqQDDVpIKMGy1LUXEQkQ8POk6tH14B5m2pYj4iOrn\nHN+fdC3aHmzZatupfyrxtkGrkgxbbUcHqH4fVyrGbgRtG/XfMztM9deOb1D9nurrWf8YuDROhq22\nnYi4Bcxl5v+bdC3aPgxbbSsRMQt8OzN3TLoWbS/22Wq7WcD+Wk2AYavt5jHg/KSL0PZj2Gq7OUTd\nsq1PmElFGLbabnZx9y/kzkyyEG0vniDTthIRL9RPV/Pun6aWxs6wlaQC7EaQpAIMW0kqwLCVpAIM\nW0kqwLCVpAIMW0kqwLCVpAIMW0kqwLCVpAL+P7HjCrL0bpNWAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot(numpy.arange(len(Et))/20.,Et)\n", "xlabel(r'$t$')\n", "ylabel(r'$E_{\\mathrm{tot}}$')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we make the equivalent of Figure 1 in [Schulz et al. (2013)](http://adsabs.harvard.edu/abs/2013MNRAS.431...49S), showing the sytem's phase space at 5 different times:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAGkCAYAAAB+TFE1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8VNX9+P/XmZlsQJJJApFFWRJQ0apAcKlajUS0Ll2U\nRdvaxSqL36ptP1WJdpEuCgi/tnZTCV1c6gKEWpeqkGBU3CFBRWVNkJ1AlslCllnO74+ZhCQzmcyQ\nWe5M3s/HYx7M3DlzzjvDnHnPPffcc5XWGiGEECLaTNEOQAghhABJSEIIIQxCEpIQQghDkIQkhBDC\nECQhCSGEMARJSEIIIQzBEu0AQkkpJXPYhWForVW0Y+ggfUMYSW99I+72kLTWAd3uv//+gMueyC2c\n9cdq3QMpdiOKxffRSHXHcuxGel/8ibuEJIQQIjZJQhJCCGEIAzYh5efnx2z9sVp3uOuP5diNJJbf\nR4k98nWHsn7V15heLFFK6Xj6e0TsUkqhDTapQfqGMAJ/fWPA7iEJIYQwFklIIiyKioqYP38+mzdv\njnYoQogYIUN2EVRRUYHVamXcuHERaa+oqIisrCy01lRVVXHXXXdFpN0OV1xxBfPnz+f6668PqHyk\n35+ONhctWsTKlSu9nisqKkIpRV1dHbW1tRQWFpKenh5QvTJkJ+JVVVUVlZWVFBQUnNDr/fWNuDox\n1ujmzJnDfffdF5Ev3I4v045kUFFRwfz583n00UfD3naHKVOmBFU+ku9PRUUFzz33HODuYD0tXbqU\nefPmkZaW1rlt9uzZPhOXEOEWjR+zAJs2bUIpxeLFizt/jJWXlzNnzhxsNhtWq5WpU6eyZMkSJk2a\n1P+Gw3myVChuQDpwt+f2HDDZT1ltZCaTSTc0NESkrby8PK9tubm5EWm7w4IFC3RxcXHA5SP5/nQo\nLy/XU6dO9do+a9Ysr22FhYXaZrMFVK/ns2iY/mH0viH8y8vLC6ov9cfy5cu9Hnf97uiII9C+0JO/\nvhELx5CWaK2Xaq2XAoVAqVJqbHRDCl5JSQk5OTmkpqaGvS2bzUZFRYXXdqvVyvr168Pe/omI5PsT\niKqqKq/3sL6+vtsek0HERf8Q/lVUVDB9+vSwt2Oz2by2zZkzh9ra2s7vDu0Z+g1HXzB0QlJKjQN2\ndTzWWlcBlcDMqAUVpNLSUgoLC5k/fz5Wq5V7772XNWvWhLXNyspKrFar1/bMzEzKy8t9vqaoqIji\n4mKWLl1KcXFx5/aKigqmTp3K7Nmz2bx5M6WlpRQXFzN//nyfdaxZs4b169dTXFxMfX19n7FG4/0J\nxPLly8nLy2PZsmWAO8558+ZFOaru4qF/iL5F8sdaZWUl8+fPp6Ghodv2nJwcKisrw95+1Ifk/N2A\nyYCzx7aNwKJeyp/QLmQk5OXl6TVr1kSkrZKSEj1+/Hiv7bNmzdKFhYVe23vuok+fPl1XVVV1Pi4v\nL9e5ubndts2aNUuXlpZ2Pl6wYIEuKirqVk8wwwyRfH+66m3ITmutKysrdUZGhs7MzOz2twaCCAzZ\nBdM/jNw3hG8lJSV6wYIFOjc3V0+dOlUXFhZGZNiuoqLCa1tGRoZev3691lrr1atX65KSEl1aWqpL\nSkr00qVLdX19fcD1++sbhp7UoLWuUErl9dg8Bfd4eUwpLy8PaFbK/PnzqaurA47vGndQSnVuz8rK\n4pFHHglJbKtWrUIpxa233gq4JyOUlJR0Pgaoq6tj7NixnY87fjFNmzaN+vp6HnroIVwuV7d6p06d\nGnAMRnt/qqqqKC4uZvfu3Tz44INMnz6dRx99lDlz5pxwnaEWT/1DeCsoKKCgoICSkhLuu+8+rrvu\nuoi023NywurVq8nNzeWyyy4DIC/P/ZHr+D7Iyclh1qxZrF27tt9tGzohAWitO09kUUrNBdZprV+P\nYkhBKy0tJTc3N6Ax11DNgqutrfXaVl9fT1ZWltf21atXo7WmuLgYq9VKZWUlQ4cO7VYmJyfHZ31w\n/O87UdF4f/qyZMmSzrYWL17MjTfeSEFBAbm5uUybNi0iMQQiHvqH8C+aP9bq6+tZsmRJt2PPXX+Y\nAowbN47Kyko2b97c75l2hk9IHZRSVmCG1vpKf+UWLlzYeT8/P98Q64+tW7cu6CnQ/TF16lSfx29q\na2t9xrF27VpWrVpFUVERaWlprFu3Lug2fR2zClSk35++lJaWcsUVV3TbNmnSJFatWsW6det8JqSy\nsjLKysoiFKG3QPqHEfuG8C/aP9YKCwtZtWpVn8evrFYrGzdu9JmQgukbMZOQgMXArL4Kde10RlFS\nUtI5CcBms1FVVdXrL4muv3J609evnPT0dHJycmhoaOj2QbbZbF5fpjabjdmzZ3sNtwE0NDRQU1Pj\nNxZwD/H1NlkiEJF+fwLR8xcmuH8J+trDBO8v+F//+tcn3PYJ6rN/GLFvCP+i+WNt6dKlFBYWdtsj\nqqqqIjc31+v7IjMzk8zMTJ/1BNM3YiIhKaXuBhZrrRs8jydrrb3nNRuUzWbrPJ6yfPly7r679yH+\nUP3KWbBgAYsWLWLRokWAe7bc5Zdf7lWusrKSjIwMr21KqW7JyNcXdIdx48Yxd+5c1qxZ021VhpKS\nkoCOI0Xj/elQU1Pj9bcVFBRwww03MGPGjG7bV69ebbiZdhD7/UP0Llo/1oqKipg1a1a3ZFRaWsrU\nqVNZvny5V/mNGzeGJHEafukgpdQMoB737CGAXGCK1nqFj7LaiH/PihUr2LVrF7m5ucyePTti57Ks\nWLGCnJwc6urq/C4dtGzZMmpqajrPc5g6dSqzZs1i1qxZ5OXlsWjRIoqLi1myZAl33XUXxcXFFBYW\ndk7T7khCy5YtIycnB6vVSl1dHSUlJZSUlLBkyRK/ywdF4/2pqqriscceo6SkhIqKCubOnUteXl7n\nRI6GhgYefPBBhg4dSnp6OjabjZkzZ3qNn/cmUksHBdo/jNo3hH8TJkxg1apVTJo0iaVLl/r9sRYq\npaWlZGZmMnnyZMCdCDdu3IhSimnTprFixYpuE56KioooLy8PeETCX98wdELqcp5FR5DKc3+61trr\nDE/pdMIoIpGQgukf0jdiU6R/rHUMyXVMggD3XlXHmo4d7S9dupSMjAzq6upQSgW1TmbMJqRgSacT\nRiGLqwrhm1wPSQghhOFJQhJCCGEIkpCEEEIYgiQkIYQQhiAJSQghhCFIQhJCCGEIkpCEEEIYgiQk\nIYQQhiAJSQghhCFIQhJCCGEIkpCEEEIYgiQkIYQQhiAJSQghhCFIQhJCCGEIkpCEEEIYgiQkIYQQ\nhiAJSQghYszHTz+No7U12mGEnCQkIYSIMW/cfz+VpaXRDiPkJCEJIUSMyZ4wgd2vvhrtMEJOEpIQ\nQsSYUeeei23btmiHEXIxkZCUUouVUtOiHYcQRiT9Y+AZMWUKLXv2RDuMkDN0QlJKFSil7gZmRDsW\nIYxG+sfANfLCC9FNTdEOI+Qs0Q7AH611KVCqlJoe7ViEMBrpHwNXyrBh0NSEo60NS1JStMMJGUPv\nIQkhhPAtxWqlcefOaIcRUpKQhBAiBg3OyKBx+/ZohxFSkpCEECIGpVitNElCEkIIEW2Ds7NxHDkS\n7TBCShKSEELEoGSrlfa9e6MdRkgZepbdiVi4cGHn/fz8fPLz86MWixg4ysrKKCsri3YYfknfiC+D\nRo2iobo62mH0KZi+obTW4Y0mBJRSa4HFWuv1fZTTsfD3iPinlEJrrSLUVp/9Q/pG/KlasoTa114j\nb73fr0XD8dc3DL2HpJSaDNwAFAAZSqnntNbLohyWEIYg/WNgS8zMxGw2RzuMkDJ0QtJaVwAVQGG0\nYxHCaKR/DGwJKSmYnc5ohxFShk5IQggRTXbepoU3sbMLxVjSmI+FodEOCwBlsWCShCSEEPHL5WrC\nbv8l7aZncVgsuNQo7GjsfIKN1QzhuwzjZ9EOE7PZHHdf4PH29wghxInRLlzNc3E5N0BiM0mWhaSq\ned2KNPI/9nMPDjQjuCs6cXo4Dx5ExdkCq5KQhBCi8a/Q8mewmDENXogl6UafxVK5mglczOdcjp1j\njOZXEQ70OEtSEsknnRS19sNBTowVQgxc2gX7r4C6hyD5e5gyP8XUSzLqYCaN03iNI7xAPe9EKFBv\nzr17MZvi6ys8vv4aIYQIVHs1VOWBBka+D2n3BfzSBDI4jcfZwa9o4YvwxeiHOTWV5FNPjUrb4SIJ\nSQgx8LRsgR1TwDIGRr0GluFBVzGEMxnMBXzET8MQYN9cGzagDh2KStvhIglJCDGw1K6Bqu9C5hwY\n/TyoE19Q43QWAul8wuKQhRcoc3srSZPOiXi74SQJSQgxcNT8F74ohNRrYMT9/a7OhIXR3EQtn3CM\nyO6tmO1tWM6YGNE2w00SkhBiYGiogO1zYeS9cMrvQlbtSApQWNnMn0JWZyDMe3ZiPuWUiLYZbpKQ\nhBDx7/Dz8MEFcPpKOOnmkFc/hZ9TQyVNEdxLslhMmCeeGbH2IkESkhAivjV+DtsKYfJrkHFpWJoY\nwkiSGM4WVoal/p70gX2YUxJQiYkRaS9SJCEJIeJXUyWsnwTjH4Ss/LA2NZ5vcIhPw9pGB7V1C+Y4\nOykWJCEJIeJVy35Yfz6cWwwjrw97czlchgsL+6gIe1usexGsaeFvJ8IkIQkh4o+zFV45D8bMhZHX\nRqzZRIbzES+Hv6Hd22DaleFvJ8IkIQkh4ovW8OJUGH41nPNARJuewixsNIS/ocot8JXLw99OhElC\nEkLEl42/AFMqnP+3flfV5nRydX0zancb6jPNHQf8l89mPI3hTkg11ZBghonxdVIsSEISQsSTD34L\nW4rgypfBnNCvqvbjYHjbEd7UbSzMhJtTXfx1HxQf7f01Fiy4SGEH5f1q26+Nb0FufK1h10EuPyGE\niA9Ht8C2p+GSxyAps19VHcPFeexjSlIypYM8daXBxzVw2w6Y4eeisVbGsJPPmMCUfsXQq9WPwinj\nwlN3lElCihMazR7+x05KOIyNGhKpYzC1ZGEjkxaGM4kJ3MMZmDjxtbuEMCTtgjXfgNOug9zr+lWV\nHRcXspMzSWCtufuiqxvOgkEl8JfdcPtY369PYjBmkvoVg19ffAZ3Lwtf/VEkQ3YxytX0JC22EdTY\nk9nNULaQw15+TgsVOKkBWgEX4MCJg2aaeYF9nKbLsbbv5fTKtij/BUKEUMlPIWkoXNr/L+qfsRsH\n7azFey8k2Qx5g+D+rb2//iRyOEpdv+Pw6YvtgANOjb/jRyB7SLGjtQbevwKc5ZAKKg0SsgDSSOAi\nxrKIVM7os5p77Id4uMnENrtClcDvMuDneWGPXojw2VUCnzwHt33W76reppa3qOMZTuu1zKvnQ/Ya\nsLVDuo+FEppo5TBH+h2LT/95BM6YHJ66DUD2kIys8RA8chr8QcFTQ6G6HFQSDPsR6lQXlqGarAQb\no/hvQMkI4KHE4bRlZlM+NBHq4BcfwCWRXRNSiNBxtMML8+Di+yGlf8eNAH7JNs4lkbMY0muZrERI\ncEJdL4MMaWSQSHq/Y/Hpzedgxv8LT90GEBN7SEqpu4FdQA5QqrWOwKnQUaI1LLkSdq6DwUA6MCIL\nrloJo6aFrJnJWaBngboX3rLBQ8/APd8KWfUiggZU/+jppZ/AoOFw3m39rqqCoyThYDln9Vn29BT4\nxlr4aIb3cxoL7eE4Trvnc7CY4StfD33dBmH4PSSl1EpgndZ6jdZ6GbAk2jGFxavPQH4yXG6CD9ZB\nUhrc+gL8RsNtR0OajLrSi4A9sOCvYalehNmA6R++VG+FHa/D1/p/vhFAIZuYwmCfzyVfAmoyPPAX\n9+OTk8Hs9F3PUWxo+jfl3KdnHoRJF4e+XgMxfEICCrTWm7s8rlRKhefbOcJcLheOqy+BkQru+DYc\ns8M3fgRrNPzVBmd9LSJxVCwFDkHGyIg0J0IrbvtHn/79QxgxBUb2/wD/uxxgOGZ+i/fxmaEzoS0F\nkjUs+rN7W4IDdlb7rquBNsy9JLZ+2boBrp0f+noNJKghO6XUYmAyUI77V9n6sER1vL0CoLLH5npg\nOhDWtsOp+ZNP2HvJJQyqr8eaBGmnj4BXPoARJ0clnkkTgd1O6p2K2PiNIiB++0dAdm2Ahmr4UWlI\nqlvBZjJI9joloroOaoDd/4IdW2H6LPf2glHwcS+rNjRjZyh+TlQ6ES//GRIUnB2ey2cYxYl8+yzX\nWt8LbFJKzVBKXa+UCteys1Yf22pwj5XHnM+XLeMFs5l3zz6b5vp6Ur7zHdJaNWw+ELVk1GH75+1A\nDddf/3lU4xBBiav+EZTHb4GpN0FiSkiqU9i5k0le26/+A+RPhlHDnPxzRRUWO9TWwrptcLiXFRsO\n00gaqSGJq9Mbj8P1d4W2TgMKag9Ja13Y5b4NKAb3LzWllA7DHlP/p80YwOrZs6latYohQHZSEmc+\n9RQjZs6MdljdTJiQAmzkP/9pAiZGOxwRmLjoH0Hb8Dg4gGsWhqS6t6mkjRZyfLydnxyF678EgwdX\n0t5eDqaxbNmiOCcb2tp917efZm5gTEhiA2Dr29B4AKbPDV2dBhWSWXZa61Kl1K2Efpig1se2rBC3\nEXbbXnyR1Oxsvrd5M4NHjOjc/nBWFpm1tYzB/ZN2xDAwjQQ1AhgFjABGem4jwDkSak5K5Ig5m8Om\nk6gmm2pO4jAncZjszseHOInqxmxaDgyFPSbYZ4L9wF7c/3bcrwE4jHvUZw9wBEhBqdXAaGAckAXZ\nJjgZ920UcErX+20MGXWU7ORqTjIdIpsjZHPYE4X7fjZH3I/bqkk/5MR0ANQBoON2yBPTAXAchH21\nsBv4AtgHZOflMWfjxjD978S0uOgfQXvtYbjsRyGr7j12cSYZXtu3HYb2JPjnrfD6b07h8OGt4HIw\ncWIC1yyFtGTvuhpowYmF0xgVsvh4ZRmc/00wx8Sk6H4J9hjS3cC9wIfAKmCl1rpjaVtfwwf9Vd9L\nvT3HzTstXLiw835+fj75+fkhDypYP29p8bk9+7zzaHn1VexAmwKHCyxOMLsXWHDf7MfvKzuYnQqT\nCSw4MeHCjBMLDiw4PTc7CThINNtpMTnBrCAB9/90At3vA+5RWwtgBlpwf8c5cGcbJ6Dd/7i63Jzd\nb9qpcGmFCzMuFBoTGgsuT4TuxyacyoxWTlCe5sye5jsOXXluqsfm1JHdZ1u8f999vL9oEc1mM/c6\nHCf8/xJKZWVllJWVRbrZoPqHEftG0DY8BbWH4PI7QlblXg6T7+M8vooDYB0CyYlw+LAGzgbMZGSA\nsxlu9zF15D0qSSWERzD2fgo7yuD3e0NXZ4QF1Te01gHfgDm4z4yZC6zF/fVU47ndFUxdQbRZ0+Px\nSmBaL2W16B/4h4Znoh2GF6fTqTdecIF+CfQzoJcPGqR3v/FGtMPqleezGPL+0PMWaP+Im75x3wVa\n/+/hkFZ5m35M1+kmr+03rtR68O/c9x9/3K5hn05IqNGNx7ROmaN1c6t3XXfoYv2QLgldcL+/Qut/\n3RK6+gzAX98Idh+wEvc00+XAcgCl1DigVruPKYVDiVJqkj4+tXWcDvPsvoHs0kvP5Y036jl0qI3h\nw8O4QGSA2j/7lPqLL8ZWV48GhowZw+Uff0xSWvxdvvkEDZz+UVkBB3dDQeiOpdhoRmHH6mOadsYQ\nGDvMff/7P7QwJM3Ojm3ZjJ0L7a0wyEf3aABm+JgccUJ2vgX73ocfvxKa+mJAsJMaSpVS6UqpaR0f\neq11VXhC6zQXKFRK5QDn4t5LE2FSVvYllDrEiBHuRSOi5sGfoP/+MJYGSGyBwTNmMGH16igGZFgD\np388/Ws4cxok+jh4c4L2U00aZp/PmRLA7vmGHDYSJkwcyzs7ob4avn2Nd/lVbKUBFzmhOoy35k64\nshBMA+dUjKCPknn2hCL2C8zT3r2eh2si1e5AljF0OHVH4d3N8OUQ/dgLyN6tsOBS2FcNTaCSB8Ff\nn8H61a+H5QBlPBgw/aP1GHz6Hjz8QUirPcghkvH9y2tHE1R7VmO48uvw1Ep4705QVnjiTu/yq6ji\nfMaGJrCy30N7PVxR2HfZOBL/0zZE0GqPgBoHF94MYV8Vrb0VHi2Abe9AI+55FWdeAH96ExIS5MpN\nwu29l2HMOZA9OqTVpjOYcWT7fG7GeDjqyVVP/BmeeQ2cx8DxoXfZbdhoxcwCpvY/KHsbvPcnuPLX\n/a8rxgycfUERlCdWACeBCscAkP0YvHIJPKLgjylw5B0YnAH3vQ7/0fDIu5AQhrXAROwq+jnkhuDL\nvodKKjnMQZ/PnZYFm5vc95UCxw7QW8HsY4Tv/9jMSaE6Gbb4uzBsPJz7/dDUF0NkD0n49N0CeKMa\n/r4R1MNw/1Ww8NR+VFj3OBz4Gdhq3HtCTUByCpz3EOTdHqKoRVw61gTaDN/9ecirzmIIjb1c3fWS\n4ZCcBJtrYZKfU5Dfo44GFGs4r/8BVZXBvjdh3uY+i8Yj2UMSvVrxLdjwM8AKv96lUdvtTGjdw3aa\n/L7uMO/wGd9hK+PZqzOobzPjtP8ALDWQmAVn/xdmabj5mCQj0beX/gGWZEgeFPKqTWjc5915UwpS\nUuAX23t/fRUt3MY25jCepF4mRwTMaYeXboKL74HU4X2Xj0OyhyT8umgk6O/DfQcdLHK0cRgTV7GT\nbKrJ4ihZ1DCUWjKoxUq952bz/OvE7EongdswDXsAsuWIkDgBz/8LLg3PNYAGk0w7Db0+P2M0PN3L\nqt4Ad/AFJ5PK9wjBUvnF18CYS+CC/+t/XTFKEpIIyIMjEniQBGAIjQzjSRzspgkzSVgZxEiSGMM5\nnMmlZJF7/IUdKzIIcSIcDjh8EG7wMa0tBLLIAFp7ff6xc+DvZfCrKvjNuO7PXcMXVOLkTfozlu2x\n5R/Qsg+uXNH/umKYJCQRtFSS+H9MhVDMKBLCny0fQHIqWMOzRN8pjKeFOr9lbh4BS47Ar8aAxeRe\n3eZ651G2mh08ocYxrL8X42vcCxsXwsUPQ3poZxHGGjmGJIQwrqcfhZNDsAfSi0SSScSFA3uvZYpO\nc1/lIn2ng/OqW7A2HeUd2likT+Ii+nlcy+WEFy6EnK/D+Ov6V1cckIQkhDCuPV/AdT8IaxMZpPEh\nL/ot88mpkJoAFQ4X5+gkdqoRzDYN6X/jpd+EzNPgkr/0v644IEN2Qghjcrmgaiecc35YmxlGNof4\nDLi+1zJjk+DQOAsh/cr86H5o/Bi+2evFCwYc2UMSQhjTvj2QORxGnBLWZs7kcpqI8OUdtj4ElX+D\nayrAJLN+OkhCEkIY09oXoC3817sayUQScFHDrrC3BcC+Ytj3FFz4DCQNzIv+9kYSkhDCmPbth4sv\nD3szCsUpnMUm/hH2tji0Bj65Bcb/H5wU/r8t1khCEkIYkuOdt3AOSY9IW+cwizq24epl1YaQaNkD\n2wth4gMw9gfhayeGSUISQhiSo/EYOn96RNoawghGchofsDg8DdS/BRsnw/hfwdgfhaeNOCAJSQhh\nSPYDB1GnToxYe6dxPfV8RJufpYROSF0JbPs2jP0FjLwptHXHGUlIQgjDcTkcOJwmzNbIXZoxmzys\njOZjfhu6Sm0lUDUbxvwCTvlp6OqNU5KQhBCG49y9G4eK/GmSZ3MXjXxODe/1v7LDv4M934IJL8Pw\nef2vbwCQhCSEMB67HWdGRsSbHcQIcrmRrfwK3Z8JDnu/CQ3FMHYVpH45dAHGOUlIQgjDafv8cxym\n6Hw9jeEm0pjAFm4O/sWuZthzBtAKIx+F1PxQhxfXJCEJIQzHUV+Pc/DgqLV/Kr/GhY19/CHwF7W8\nAgdPhuQ8GPk0DA7vkkfxSNayE0IYjr2mBmdiYtTaT2Iop/Iw25iNGc0Ier9onqP9VUzNc1AuOyrj\nHzBIVu0+UTGzh6SUWqyUmhbtOIQwmnjsGy6TCaJwDKmrFMYxgeXU8QR7mYeTps7nnLTS5vwXzW25\nOBwz0SkLIPOAJKN+MvweklKqAJgCzADWRjkcIQwjnvtGe00Njra2aIfBYCZzBpv5gqvYzSSSsaI4\nQoJuIUFlkmL+PklJv4p2mHHD8AlJa10KlCqlInPKthAxIq77RkoKOopDdj2N4RXa2Ew72zCTSIq6\nEqUGxdAYU2wwfEISQgw8LqDdAHtIXSUxiSQmRTuMuCb5XQhhONpsRiclRTsMEWGSkIQQhmPOyqL1\nyJFohyEiTBKSEMJwTOnptLe0RDsMEWERP4aklJoD5AHa19Oe7Uu01rtPpP6FCxd23s/Pzyc/P/9E\nqhEiKGVlZZSVlfWrDukbx5nT02ltauq7oDC8YPqG0trXZ994lFJrgcVa6/V+yuhY+XtEfFNKobVW\nEWor7vqG7bPPeOvqq7l29+5ohyJCzF/fkCE7IYThpIweTfuxY9EOQ0SY4ad9K6UmAzcABUCGUuo5\nrfWyKIclRNTFc99IHDIEu9a4nE5MZnO0wxEREjNDdoGItWEJEb8iOWQXiFjsG8+ecgrXvP02qaNH\nRzsUEUIyZCeEiD2pqRzZtCnaUYgIkoQkhDCkxOxsDpeXRzsMEUGSkIQQhpQ5eTI127ZFOwwRQZKQ\nhBCGlHXOOVRLQhpQJCEJIQxpdH4+DYcPRzsMEUGSkIQQhpQ+dixOpWhraIh2KCJCJCEJIQzrxv/+\nl8QhQ6IdhogQOQ9JiDCQ85CE8E3OQxJCCGF4kpCEEEIYgiQkIYQQhiAJSQghhCFIQhJCCGEIkpCE\nEEIYgiQkIYQQhiAJSQghhCFIQhJCCGEIkpCEEEIYgiQkIYQQhiAJSQghhCFIQhJCCGEIkpCEEEIY\ngiQkIYQQhmCJdgB9UUqlA3M9D6cCi7XWFVEMSQhDkL4h4o3hExKwRGs9H0ApNQ7YpJSaorXeHd2w\nhIg66Rsirhh6yM7TyXZ1PNZaVwGVwMyoBSWEAUjfEPHI0AkJsAKLfWzPinQgQhiM9A0RdwydkDzj\n4Xk9Nk8B1kYhHCEMQ/qGiEdKax3tGAKmlJoLzNBaX9nL8zqW/h4Rv5RSaK1VBNuTviFigr++EQuT\nGgBQSlnx0+E6LFy4sPN+fn4++fn54Q1MCKCsrIyysrKotC19QxhZMH0j4ntISqk5uIcafDWsPNuX\n9JwppJQpQeyQAAAgAElEQVR6FLhHa93gp275FSgM4UT2kKRviIHAX9+IiSE7pdTdwKqOjqiUmuzr\nfAvpdMIoIjVkJ31DxBp/fcPQkxoAlFIzgHKgTimVrpSagvfBXCEGHOkbIt4Y+hiS51yLVRwfwugY\ntpgetaCEMADpGyIexcSQXaBkWEIYRaRn2fVF+oYwipgeshNCCDEwSEISQghhCJKQhBBCGIIkJCGE\nEIYgCUkIIYQhSEISQghhCAM2IYV73bFw1h+rdYe7/liO3Uhi+X2U2CNfdyjrl4QUg/XHat3hrj+W\nYzeSWH4fJfbI1x3K+gdsQhJCCGEskpCEEEIYQtwtHRTtGIToYLSlg6IdgxAdYvryE0IIIeKfDNkJ\nIYQwBElIIiKUUouVUtOiHYcQRiN94zhDXw9JBE4ptRhYq7Ve34867gZ2ATlAqa8rj55AnQXAFGAG\nsLa/9fXSRjow1/NwKrA4FLH3qLse97WGHtNal4aibhEZ0jdiqG9orQf8DUgH7gbmACuBgjDUfTfw\nHDA5xLEXeOreAUzrRz0rgUldHq8NcZxr+xNfH3U/2uX+OKAWGBuiuhf3qNsFpIXj7zDiTfqG9A0/\ndYe8b4T8DYjFWzi/dML5gejRTr8+1EBNz7hD2UnC1ek87+ndPbZtBO4KUf01XeP2fDYmhaLuWLhJ\n35C+4e99CXXfkGNIbnM6xnC11lWebTn9rdRzmeldHY89dVcCM/tbdyh5hg4qe2zu2A03Oiuw2Mf2\nrBDVn6c9Qz1KqRzclwnv+V7FM+kb0jd6E/K+IQnJLVxfOuH+QISK1ce2GkLwxRNu2j0entdj8xRC\nNCavtd7d5eFc4B6tdUMo6o4R0je8Sd8gPH1DEhLh+9IJ9wcihDKjHUB/aK03d9xXSs0F1mmtXw9V\n/UqpcZ6D2uOAolDVGwukb0jf8CfUfUNm2Xl4hhBm4n5jHwxVveH+QIRIrY9tRvul2iellBWYobW+\nMpT1eoaTlno+I+VKqSkDaS9J+oYX6Rseoe4bkpA8wv2lE+gHQik1B/cvR19LaCjP9iU9frn2Vz2+\nhyZi7VjJYmBWKCtUSqVrrW3g/owopeqBez23AUH6hvQNX8LRN+IyIQX7wQ3mje1HpwjoA6G1LiLC\nw0Ja61KlVM+hiRzcs4ligmfYYHHHF6VSarLu5/kWngPa6/Ae2vb1BRUTpG8ER/pGr3WGp2+Eeqph\nrN1wn6vg6rFtI/BICNu4my7TWQnx+RaeOvs7tfU5up9r8WGI4pqM+wvHCXxIiKac9mhjhuf/Md1z\nmwLcGoJ6x/WMF/cB7ctC/TcY8SZ9o/P10je86w1L3xjwi6t6hiFmaK2XddlWA8zUIRjPVkrNwL3b\nv9GzKReYorVe0d+6PfVPBm7A3bHLgee6/i1B1JMOFOLuGOd66tns/1XR12X6cMcHueOX+HTdjzPz\nu9RfgPuLw4a7M6/TWq/pb72xQPpGZz3SN3zXH/K+MeATEoTvSyfcHwghwk36hogkSUhCCCEMQc5D\nEkIIYQiSkIQQQhiCJCQhhBCGIAlJCCGEIUhCEkIIYQiSkIQQQhiCJCQhhBCGIAlJCCGEIUhCEkII\nYQiSkIQQQhiCJCQhhBCGIAlJCCGEIUhCEkIIYQiSkIQQQhhCXF3CXCkl19IQhqG1VtGOoYP0DWEk\nvfWNuNtDCvRSuffff39ILxUcyfpjte6BFLsRxeL7aKS6Yzl2I70v/sRdQhJCCBGbJCEJIYQwhAGb\nkPLz82O2/litO9z1x3LsRhLL76PEHvm6Q1m/6mtML5YopXQ8/T0idiml0Aab1CB9QxiBv74xYPeQ\nhBBCGIskJCGEEIYgCUkIIYQhSEISQghhCIZfqUEplQ7M9TycCizWWldEMSQRJhUVFVitVsaNGxfx\ntmfPns3KlSu7bSsqKiIrKwutNVVVVdx1110Rj6sv0j/iQzQ++xUVFSxatMjrcw/uz77NZkNrza5d\nu1iwYEG32IqKigDYtGkTSikWL15Menp6/4MK59m7obgBj3a5Pw6oBcb2UlaL2JWXl6eLi4sj3u6m\nTZu0yWTqtm358uW6qKio83F5ebmeN29ewHV6PouG6R/SN4wtkp/98vJyvWDBAr1gwQI9depUr+cf\neuihbo9Xr16tc3NzOx8vX7682/PLly/v9nxf/PWNqCccfzdPB7u7x7aNwF29lA/4TRHGYzKZdEND\nQ8TbLSkp8UpIeXl5XuVC1elCdQumf0jfMLZofPbLy8t9JqTx48d3S46VlZXaZDJpm82m6+vrvRKS\n1lpnZGTo0tLSgNr11zeMfgzJCiz2sT0r0oGI8CopKSEnJ4fU1NSItrtmzRoKCgo6vrQBsNlsVFR4\nj3pZrVbWr18fyfD6Iv0jDkTrs9+bdevWcf3113c+3rVrF1arlbS0NCorK5k/fz4NDQ3dXpOTk0Nl\nZWW/2zZ0QtLusfC8HpunAGujEI4Ig9LSUgoLC5k/fz5Wq5V7772XNWvWRKTtiooKpkyZArhP1utQ\nWVmJ1Wr1Kp+ZmUl5eXlEYguE9I/YFs3Pvj9jx47t9vihhx5i9erVAEyePJlNmzaRlpbWrUxlZSU5\nOTn9btvwkxq01ps77iul5gLrtNavRzEkEUIFBQUUFBRQUlLCfffdx3XXXRextquqqpg8ebLX9tra\nWjIzM722W61WampqIhFawKR/xK5ofvYDUVxcTElJCYWFhVx22WWd2ydNmtSt3OrVq8nNzWXatGn9\nbtPwCamDUsoKzNBaX+mv3MKFCzvv5+fnD5j1x2JdeXk5BQUFfZabP38+dXV1AN2G2eD4Xo7Wmqys\nLB555JFe6ykuLmbGjBn9iLi7srIyysrKQlZfsALpH9I3jCnSn/1AzZgxgxkzZrB06VLWrl3L4sXe\no8P19fUsWbLE71B2UH2jt4NLRrsBjwJpfZQJ6KCaMJaSkhI9fvz4iLVXWVmpy8vLu23rOqmhpKRE\nZ2Zmer1u+vTpeunSpQG1QYRm2WkdWP+QvhFdB95/X9tbWry2R/qz31Vvkxp86W3Swrx583RVVVVQ\n7frrGzGxh6SUuhv3+RUNnseTtZxrETfWrVvXeSwnEkpKSrDZbJSWlgLHf20uW7YMq9XKrFmzqK+v\n93pdbW1tROMMlPQP49u3YQNpY8ZgSU7utr23z76jtdWrbCRUVFRQUFBAbW1tt+05OTmsW7eu27Dc\n0qVLKSws9Drm1B+GT0hKqRlAOVDnOQkwF/eBW+lwcaKkpIT58+cD7hluVVVVXuPUHboOW/RG9zFs\nMWfOHK9thYWF3U58zcnJoaGhodvBW5vNFpJx8lCS/hEbzv2///O5vbfP/uobb+Sbjz9OcpeTTUPx\n2e9LbW0t8+bN89peWVnJ+PHjOx8XFRUxa9asbsmotLQ0oKFHv3rbdTLCDfd5Fi7A6bl13J/WS/mg\ndh2FMYwfP15XVFRorb1PyosUz+UZOhUVFenCwsLOx+Xl5Xr+/PkB10fkzkMKqH9I3wgdp8Oh35w9\nOyR1RfOzv27dOp/n23X93GvtPnE8MzNT22w2rbV7mLHrkHd9fb0uKSkJyXlIcj0kEXUrVqxg165d\n5ObmMnv2bK8ppeFUWlrKY489RnFxMTNnzmTevHmde0ErVqwgJyeHurq6oJcOkushxbdj+/czaNSo\nPsu5HA42L1vGlMJCn89H47NfVVXFY489RklJCRUVFcydO5e8vDxuvfVWABoaGnjsscc6PsNUVlay\nYMECxo4dS1VVFbm5ud1Ok9Bao5Sirq4uoPj99Q1JSEKEgSQk41l700185eGHScmK3HnDWmv2vPIK\nY66+OmJtGp0kJCEiTBKS8bTV15Pk44RnEVmSkISIMElIQvgmlzAXQog+1H/4Ic27dkU7jAFNEpIQ\nQgDOxkZcLS1+y7R8/jnOpqYIRTTwSEIawFp37qR9714ctbXsv/12tNbs/f73ox2WEL367LnnqNm2\nLSx1Z02bRuqXvuS3TFNZGfZ9+/yWcdntHHzyyVCGNmBIQhpgWrdupX3PHgDatm+n/YsvsGRmMuyu\nu1BKkf2LXwCg29vZf+ON0QxVCC/W3FySfSx8GynDbruN5NNP77OcdrkiEE38kUkNA8SRJ55g2Pe+\nh+2VVzBbrQz58pf7fI39wAESRo7EvncvjkOHSDn33AhEGh9kUkP0vfGNb3DR009jGTw42qGILvz1\nDcMvHSRCw+UZ906/6qqAX5MwcqT7tQ0NOI8cQTscaKcTU1JSWGIUIpQueuYZLIMGndBrXU1NmIYM\nCXFEoi+yhxSntl57Lae/9FJI62x6/nkce/diveOOkNYbj2QPKXa5bDYabrkFq+eidCK0ZA9pAEo6\n9dSQ1znkm98MeZ1CGI0pPV2SUZTIpIY40bBhQ7fH437/+7C256yt5cDVVyO/ukUsObRkCe3794et\n/vrZs9FtbWGrP95JQopxLs+Hv/aFFyI6s8ecmUnSGWeA3R6xNoXor7SvfhXLsGHhq/8vf0H5OMba\n8tFH1Aew17Xjr3+lfsuWcIQWE+QYUgyz19Wx9eabOev556MdiuhBjiHFmX1V8Nyj8LMlvp+/+svw\nv3d7fbn98GEc1dWknHWW32Yad+wgZeTIuJ4ZKGvZxZn26moSs7OjHUY3zupq7BUVJF95ZbRDMQRJ\nSJFX/dJLZFx8MQnBLKDa3g5fVMKEPs4tcrng6CHIHtm/IIWsZRdvtt92Gy6HI9phdOdwoI8di3YU\nYgBzNjejg+0Xpa/Cnbf0Xc5kkmQUAbKHFCPajx7lyPPPM8pzES2janv7bXRzM8lXXBHtUKJK9pCM\npfWVV0gO4hw8v2w1sGENXDMnNPUNMLKHFAfMKSmkjBsX7TD6ZD75ZMxjx+Kqq4t2KEJ0av/gg+Be\nsOIXcPSA7+cSk2H42N5fezh8s/jinSQkgzv09NPse+QRzIMHk1lQEO1w+mQZMwZlMtF0//3RDkWI\nTmm+Po8N9fDHn/t+wZXfg4yTfD+XmAR5030/19oCf1xwYkF2VLFvH3v++MeAyrbF2Q8/GbIzqJa9\ne0k66SRUQgJojTKF8beD7RCYzPDuv2D8xZA2HIYaf2/MyGTILga4XLDtY5g4KbjX/fZi+OWGvsud\nIGdLC81bt5I2eXKfZV+97jqmP/ss5hhazivmh+yUUouVUtOiHUckHXjySZp37kQpFZ5k1NIA9Qfh\nqdug8j3Y/SEMy4W2ZqjeBs/Og7YmaPd/fZi+2L797RAFLHozEPtHV+0Hehla64vJ5J2MHHY4stv/\n6040Gf34h7D10+NNvf8+jk2bvIqZU1ICSkYAX/3Pf2IqGfXF0HtISqkCYAowF5intV7fR/mY/xVY\n/fLLZF9zTXgbKfkLWEdC+zEYOwVGnuFdpv0YbHoWzIkw4RLIGH1CTTkqK7Hk5PQz4NgTiT2kYPpH\nPPSN3lTNm8fJDzxAwtCh/a/s0E54+98wo8sQ39+mw5wXISE5sDr2boP1z8D3F/ot5ti8GSwWLH1c\ng6k/6nbuJGHIEIYMHx62NoIV8+chKaXWAosHQkLa+cAD5N53H0qF4btsz0ew813IGg1nXg6WxL5f\nU/k2HN0FY86Hk047oWbbnn0G7HaSvvu9E3p9LIrkkF0g/cOofcPlcPDi177GN155JbINH9kP//0b\n3PpA32VdTtAuOPQBjLqoRz3bITMHzBZYeBEsfBucDqg5CNmnhCf2IOx4/nkGZWcz6sILox1KJ0lI\nBmdvaMBeV8egMWPC18jKe+Gyee7OdVJu8K//9w/gxiJQJvfxpiDotjbQGsxm9zGxAUASUuDaGxtJ\nTE3tfHxw/XqOvPceZ993X/gadTrhYBWcPD6w8q118EkRnHtP9+1rfw3n/RCsp7j7VpB9YyCK+WNI\n8a7p8885sm5deCpvOAp/uhEu+p57z+hEkhHAd/4F7y2Hj4uDfqlKSsJ15AjNP5TLowtvXZMRwPD8\nfM746U/D26jZ3HsycjmgenP3bckZ3snohWthWqE7GUFgycjlgprDAYWonc6AysUTSUhRZLfZ+Pi2\n28g4/3zGhOOE1w+fhzf+AXOKYNRE90Hc/rjoR3DoE9hZFvRLzaecwpAnn+5f+2JAUCYTlpSU/lfk\ncMCHbwZW9rGrj99vrYPKF3svu/NR9xDeV58BS5ATCvbshCf/EFDRAz/9Kcd8THqIZ5KQoighPZ3x\n99zTd8ETsXIhTP0GfO0eSEnts3jAvvpbGJ8Ppb8FR3vQL29/6nHsr0X4eIGIW01/+AOOqirfT7Ye\ng/J3vLcf3QP//HH3bd/99/H7g4bBBb/s/rzLDvv+6r5vGQIoSDyBfjX2VPjJ4oCKjvrTnxiUlxd8\nGzFMElIUNG7bRuVjjwEwKNSrL7Qdgz/eCKddCOGYGNEh/WTY+I+gX2a56losF30lDAGJgSjp8ssx\ne2aQOb7Z42TVIWkwr9D7RVmnwHU9jk8NynD/63LClkXer1EmSPScKDv2pr771keRu8Cfy25n/4t+\n9uhiSNxdMXbhwoWd9/Pz88nPz49aLL1JHjmSoRdfHPqKD2wHcwL88M+QFr5rvgBw1kyo+wJs+yF9\nVMAvM2VlAe49pcSb4ueYUllZGWVlZdEOw69Y6BvBSuhyOQfzU//xLlC5xT2bdHSXKygrBdZeVmFQ\nJrCe7WO7GbJn9h3Quh/A1HuhbnfvZezt7ll5XYfQ9++BtS/Azbf33UYPLoeD5t1+2ouyYPqGzLKL\noIZPP8XR2EjmBReEvvKDO+HgDvfw3MQTSHaOg+DcC7RC0iWBv+6Z2XDDs0Efn2pf/giJc28LLsYY\nIrPsTlzNhg0cq6zklO+F4DSBijcgKQXOOM/7uYMfQupoGNIjObXXwJYfwJQXoW0rNL8CmV0mWbia\noeU9GFwAu38JI28/vvfUmw9Wu09Gv/SHsGopjDkTzuty3KrlGOzcCmdNOeE/NVbE7LRvpdRk4Abg\nbqAceE5rvcxPeUN3OtuWLTibmkKfkNpb4cXfw0U3wPAgZ9HZ90DDH0DvB8t4sGSB/XlIvAYG+xju\n8OXVe+D8/wcZYwNuVldXY//lPSQ+9q/g4o0REToxNuD+YfS+0ZW9vh5HYyMpp4ThPJ5P/wcoOPMq\nePE7MPoyOMfHhCKt3XtSrhZc7etRSV9GqUz3c44j0PRfsN4KzZ/BoNPce1D+2NvcEyESg5us0b72\nNRKviK9rjMVsQgqWUTudvaGBD+fO5cJnnw1PA7+aDve/4h4GCFTLB9BaBo4dkP4zSPS+QJluvB5S\nfomy9LGMiW0ffPIMXHx3UGHrujpURkZQr4kVspadQdk8ywylB35tI6fjKZRpEiaTe0UFl/N1tGs7\n5oR54Yiwm5YHfkvyvT8P71qWESYJyQDaamtJyswMbaVN9fDPn8Edfw/udfYDUHMXpM2BQZf1Wky3\n/hNtOQv0fkwJ3/Bf56erIXc6JKcHFUr7bT8k8ZHgJ0cYnSSkGGKvhaOvQfa10LQO0q/3KtLIH0jk\nYpI4F62bgWMo5TlOq53Q8gYc2QQpubD9Kcj7C6R4kl57AySmhS/8ujoSYuiHnZwYGyXt9fUc+N//\nAEKfjOqrYe/n8J3fBf6atq3QsgmOLoThT/tNRgAq+WaUaSS49uGy93EO0RkzYNWNgcfiYflxcHtV\nQgSsqsI967QvygLmIYAG3Yqz5Ta0tncrMphbSMQ9UqDUYJQahpNq95O6Ddo+hNF3Qfb1cPGa48mo\n+RC8eUcI/6junM3NbJ87N2z1R5okpDBxtrXhbGvD0dAQngZaGuHoPsgcEfhrapZC01oYsTzglyjT\nSFTit3Dpclyuj/0UVPDdV2DTY+C0916uB9PpE3H+/Gdoe+CvEfGvfu1amgK9qN6Pb4B2H+fEff4W\nNB7t+/WWNMj+GpjTwPptTIk/4pgq5hjP08S7AJhIQ3kmJVfzJHaqqeVeABxsoD1V+Z4KPng4XP74\n8cd7P3RPLQ8R8+DBnLlqVcjqizZJSGGy9Y9/pP6jjxh9Y/B7DX65XPDyI9BUB1+ZFdhr2nbBgdth\n5N9h2L1BN6lUJuaEu3DaH0DrPn5xuhxw9PPg6p/5rQGzxp0ITOKoUSRkZwdW+BcPQ6KPhYKvvhOG\n9lilvuUw1H/mtzpl/hIpfJ0UvkoT71HFv2lkZ+fzg5iImXSG4R4qN6vLSTD/uLfqutv2CrR6fqSu\nXOC+BAzAx2XuE3kX3+5ezWGAkmNIseazd+CN5+B7v4HBARyvcbXB4d/BoPMh/dp+Na11Na32b5GS\nWNp7of0fQEstjP9qUHU7f30v5vt9nJAYo+QYUmjotjacdXVYPCe/ul57GVwuTFd9rXvBQ1/A+lXw\n7bt8V/TqfLigEJy10LQLTA2QcTEMOhlMg/3G0Eo1CaTRyGGshHAB5Opd7gthmkzuH5kXzwRzkvuE\n3j5opxPbv/+NNRRT4yNMJjVE0J7Vq0mwWhlx+eWhr9zlCu58H3s1HPoFDL0DUs7qtZjGRRXvsJdN\ntNLK2XydEUz0Wdbh/B8uvZlEi5+VmD98GDJPg9zAk5LrrTJMX8kPuLzRSUIKjbbyclpLS0m/232s\nUVcfdl9B+aQe1/dpa4V9OyC3l895q637hJuGzUA1tL6LK90CCTegTGNQ+N5Tb6WBzyhmCjf7fN7J\nYRzsJ4kpcGw9tJRA5u+g4T1ID/2lH7TTie1f/8J6yy0hrzvcJCFF0LH9+zElJZEciouF9fTU7+DU\nqXBegF/0tc+C2QrpvZc/SAXN1LGd17iQO3Dg4ACfsp/PuArvCQdaa1pdPyPRdAtmdabvStub4cgW\nGHoGJAWx3td7G+CCMKxgEQWSkKKs/RgkDgqoqHbuoM20k3ZVThMmLFyEhSFkcvwk1X18zmGqyMN9\nMusnrGI0F5LOKNrYSyPrGMRwBtHlZFdnE+xdBmMXdomrEaqeh9O+G/jf4nJBUyOkBTeD1ahkll0E\ntB45ws4nnmDQqFHhSUafbIArfxB4Mto9D1SS32S0hZVU8ToppHM1S7ByMkMZy+lMYyozeZb7aKO5\n22uUUiSbfkOr6w56/YJLHOxevr+lJsA/zmPdy+4VmoXor399C1obAyqqzBNIVleRyn1k8R1SOY3B\njO1WxsrJ7OBI5+Ph5JHIEAC28QRtpHcmI7vjTzicz7pn7nVNRuBemsjif4jQyyfl8Nj/F9xrYpQk\npBAxJSWREq7LBNcfgZoD7hUZAtF+EAadB0N8DxXYaaaStZhRXMCdjKL7isKJpDCMcXyZG9nHVuy0\ndXteqSEkmh6gVT3Yewynz4KSHwUWbwenEz7wsTqzMJzKtWujHYJ/c/8LyT32zlv2grP3PqRQJDGa\nFE7C1WWZTwdOhpDKjRxfe/FDNnCYPbhwUssRMjh+JdkEy51YzD4mM+3+JzjqIddznpNtP+zb2Pff\ncs5UuPs3fZcLI0dbG/VffBH2diQhhcCnv/89B9evZ9QVV4SngU0l7jXqRgVwdUtHA+y9FzJnQILv\n9bXK+RtNfMHpzMTkZ33dMZzNYSrZzntez1nUebg4jIteptWmZEL+Q3CovO+YOyx8CC4MYh09ETX7\n3n0X7XJFOwz/Vn0dWuuPPz60EhreBqfNZ/F9rGInz2OnmTd5gs+oYBPv8kee5F0+6lZ2POcymjMw\nYSaPBQyi+49RjWf0oPU9aPecLjF4PFi6JMnWemg42O8/80S8etVVOFoD/IEL1O7cyafPPRfGiNzk\nGFIIuJxOlMmECsflHt55yX3OUcG3Ait/4GEYdCZYfU+q2MEaBpHNKAI7VtPGMb5gC0fYz0Vc1+25\ndkpwcpAUehkPP1wBTQcg95rAYgdY9SR8+VI4eXTfZQ1MjiEdd2DdOkZOn953wf74+51w429gsPX4\nNmc7mHtMB697CizDIdW7fzTwGR+xki/xXTLIxYmTT9nPqYwgmQQ0GoXiLXZwEkm4cHF6j6E9J8c4\nxpu08z5Z3O8+YVYNgsRejrd2VfR9mPN43+VCwNHaiiU5OSJt9STHkMKkrb6eQ2+/jclsDk8yaqiF\n5CEwNcA9r6bN0LK912S0heXYsXkloybaWUE5d/Aa73OQeZR0PpfEIKwMJ4lBuOh+Qp+FqbjYj5Ne\nfuVlTYTPnwgs9g5jciA1fMusiMirfucdXOG+HPdVd3RPRgD2RjjS4+TajJvQQy5DY6eGazjKClpw\nn5eUxhl8hYU4GEw7bZgxU0cL7biPa/6J9zlII+s5QDOQSTotnuHseqppp5lWqmih3p2MAJLOPZ6M\nds4DZ/djst1c62PmaqMN6r2Pxbbc+n2cO7b3+bb0JlrJqC+yh9QPTXv2cOD11zn1+2G6rs/eHfDZ\n+3DlTYGV/+hSOHu9z5WHj3GAZg6TQjZDOH79osM08nve4zucRRaDScbCK+xmDOlMYiipuH9hvsPL\nTOAchnFyt3rb+B8umkmhl5N0v1jnXjl5bIBJtaUFnvsX/CC2L00he0hhZqtxX1Yi2c9MuoadcOR9\nyP1O5ybtagDbGFRGHRoXDg5jJhMTxy9F/i5vMJYJjMD3AqzP8hnnMoLRpPJ3XuIshnKUj0mjha9w\nO5YudTn1VjT1WNQF0LLLvdZdIGoOuFcI3/4pNDfCVQGOkMQA2UMKgwOvv46jpSV8yaj1GGwvDzwZ\n2TbAWSU+k5FGU87vUKhuyWg/NlbyCbcymbMZzihSySKFm5jIEVr4L8cPYl7INazjWbZT0a3uBC5A\n4efXlnUC6CBmzlksgGG+x0WAnO3tAR1TajtyhB2LA7uEt5e/PQhvvua+/3oxfNzLBJgqz2Wh0sZ3\nS0YAypQG1lr3fUwkMKIzGb3JEzRylESGUuWZUbeYt3idQ/yZ3Tg9x4VuYCJ/5yMSsPAtrmQ8p3I1\nc7iMn2Gm54oRGjpGFnwlo5e/6r4sRU8HdsAXW+CSa08oGbXs3s2x7cf3oNqOHMH2ySdB1xNpkpBO\nkCkxEZMljBfc1RqcQXyRb58L7Qd8PlXNO5zLgwxlUrftf+UtvsZEJuA9Tf1axpKEmYoukxau4Qfk\n8D74BUcAACAASURBVKVu5Uxk0sTfOUYv62klZ8DBDYH/HQkJsOntwMsLQ9j8l79Q5VlI2J+E9HSG\nnehJ43Puhq949rS/ORfO81GPywFfvHH88bEqaNzarUjX4XUXx9dQPJ8ZpDKUM5jAJ7iT1p1cwETS\nSSOpY5oCCsWD5AOwnm3so4E9noVW3+cRbOztrNOsJmJRF7l/lB172Tve6avcU8F7OutSOP9r3tsD\n1LZvH61dZsW1HjxI42f+l0wyBK113Nzcf074lc2fr49+/HF4G5l/idbtbYGVdbb7ffpdfbtu0Ls7\nH7dqu16iS/Qn+kC3cpc21+nfNzfrnQ6H1lrro7pFf6yP6hbtfnxMN+k/6/u86m/X27S9R13dbP6z\n1ntfD+xviROez2LU+0THLVJ9I1jtR4/2+pzr2DHtcjoDq+gfN2l9aJv39pq3tT78P613XOKuUzu1\nXX+ktda6TR/Vn/n4PHe1Xtd1e3zM6dQftB6PqVwf1k/qcl2pD2i7btdb9IbjTetXtVO3uB84j2ld\n/3DvDX2yWus97/mNJV746xtyDOkEuJxOTOY+rhDZHx+UwMm5MHJcAMG0w5vDIN/3VNadLCeTc8nk\n+EX2PucgB7BRwPGL8p3XfpRhjgQu1BZebIOZSQncNTiRn/A+MxnLxbinkB9iL220MYbjU9BdNFPH\nnWTRy3WZGr6AYwdg+JcD+OOBd8qg4n340YLAyhuQHEMKzLZvf5vcRx/FkuY9kcX+4K8xXToN80Vf\nCVl7Lupp4nc4OY8MZnd7rp4GUkgmyTPsptFcqnfTWJPNbYMSSVIwDDNLm5y8Psy9xJAdF3acDCIB\nJw52Us5pnMdHLCOLbEYwEzM+jnNp7d5zS/Ms0XV0u3tpo56XU49DsnRQiDjtdnY+8wynhXtBww9K\n4JTxMGJs32Vd7aASfC59346NLfyOKSzt3NZACx+ymzMZxXDSaNYunqCZFKeZH1jcHefBxjZOtZiY\nmeLudE9RxU24k2MtR3ifN7iKmd3aquHHpPETEvCRRBt2u5dLOecngf39Dod7bbLBQwIrb0CSkIzL\nSSNNbEaTRg17GMdVvMQzpJHDOMYypstx1m3agcUFOSb3D9CO4b4LG+o5LbWJ6SqRb+O9KrmT9s7j\nSTUsIp1bsHSUa9oMTZ9A9cdw9v/f3nmHt1ldDfx3JW/HiUf23mGThEBYhZCQQoBSyixllBYChVL4\nWhpKWGWPQtlQRsteAcImgZBAIJAQsnfi2Ikdz3hvWes93x+v7Fi2ZEu2ZEvJ/T2PHuu997733lfW\n0bnj3HMeaXNvh8x7GWacBX29261+5x3cZWWk/yV88ZdCgTZqiCbysqB4T2DKCCDnESj/2mdWBcsZ\nxMleaZ+wjjzKGYg5Ij3VKOFLaeC31kQ2lMPSIrg1JZ5FNhhTZMaYGcE+Vyfp9ENhpbyFGxWAeE6g\nkWW++xiXAsUrAnseMBXSlWcHXl6zX1Dz9tu4KypCW6ltM5R6x/+ykkIcY3BSiSKBYnZwOhcyhoPY\niXf7E1QMPys7m8TJsyUwu8bG8zYHhYUpjDFcbKUSAANhD5UUUUgN1V7GDQbVuGUdNuMGMyFuAKRM\n7JwyAhg6EhLbzrr6XHwxaddf37k6W2ErK8PtK8ZUmNEKKUCcDQ3U5uaGf3aUvRnGHhFgp6qg/Afo\n29ak2k4FtWxjAPt82VVTz8mM5wrPOaRq3JxpjeMTywASlOKyH+HJrTDqXXgm1cqrqVbqDaEPcVzA\nPoVyCEdSg/cSYTy/wI6fTdOEDNP7t9vuO79N+QR4teMNcs3+hSUlxWNl2UkWXAeuVt+xuBGQ0tb7\nRyKD6cc0+jEJOw2sZjMb2EJBC4X0Vq3Q4BTuKYilCDfPFStG2eJZaxdcqXaeyxvMnTIegAYc/MAu\nqqmmgQacNOLw+IHsx0PEq9NItDwFRj1YLZDs3/t+h5ww3e/qQajOQ25+5RX2rg3Cy0qI0AopQGpz\ncsj5woeVTKjZvQ3GBaiQcp+FkTf4yRQEl5droFXs4kuPCxQXwlzK+DXJzMuF4Z/C+yfDvZNgcjrE\nWSx81mDhH2UwlEROZ9/atgMXK/E2uY1lEL1oxwQ+/XAIZsXo+ougxve+WLgwXC4W33JLt7ap2Uev\nX/0Kq4+9pGZEYLlnoPL5k7CmlTxO+TOULIPcj6F2K1SuBGsKJJh7pY28TSObqeGT5ltS6MsojuZ4\njuIsTsJGAvexnRVUk4mNLwrdJO2J5Xh3Is8NhQ8yLQx1x3NlYwIX2WNZ6XbRKEIv4vkdRyEk4yaW\nYraQhw8/dY41YPuobfr3V0B9QZCfWPg4es4cBh97bLe3qxVSADjq6ohJTmbijQFGhewsW1bDoVMh\n1kf0S1/EZECGb5cse/mQsXjv2VRTxx885qpFuMjEwQBHHBVOOL4vvLAIFq+DN06GcW/Cn1JghwP6\nSDzTGUgOZrTYkYxmJG396hVxEW78RJSt22W+AuW/n3S7u31LTAwTztZLhRHHg3PMZVy3C3ZtNtNm\nXg0TT/Mu1/9QSDsc+k4xz/aI5/yPu+noghUrGcS1+O6WUkQWWwAooJwiFL3IYAMucNp4Lc/NtkLF\nAwvAXQyXDYS/p8ITdgv3j1Osdbl5XbLIxfQsnkgcccQyjKMYw8m4qKWMFrP9hJOg15/aPuNxz0Dy\nENg8v+0sr8mpcqOtM59eVKEVUgBUbt/Oni+/DH9DvdMgNcDQFfW7oH6nT2MGgFo2I+w7x5RHGcVU\nEuP5l39MLW8wiOxaxdoaePd42LIX3lsDs56Ex0+E0XGKq/ooFFCMnZ2eJYgYYiilDBveAjKED1H4\nORzZZ3xbNy7t8b/HYOV3HZcLIXs3bcIaH99xQU33cvp55lJeTCxcerOZFp8IVh/Le4kDIHko9D4M\n0j3e7qt/AfYlJHARsQwiAdOVTzGZ1FNGMik04mAZ27ibYziL3kwkkbvTM3jt6HiWnQIHO2H1EogH\nYhVUjYdP6qCfIwGHBSo8sjaSgfSjD5/zXwQhkzewkNjxM8Z6luBsFaYyfX8OPPNr2JsDz19n5t1+\nMdR5rxqIy4Xtmj8G+YFGLmE82Rk6lFJzgGxgNLBERNZ1cEvIEBH6TZzIgClTwt/Yf++DW/8TWFlX\nHfQ/w6dnhkp+IJ1jvMxN97CXEzkIhcKFQSF2MrBii4cbxsCsl+Cr6+Haf8GF06GvFS5fCHsmwItl\n8OyoGDLZF18miWRqqSGxhbDVsoAkTiERH+vj6UdCfU7AHwVrfoQzLwq8fAgQtxsjCuMx9aR8dAuT\n/CwdiUBdGaT0a//+9G1el1X8RDWVKNLowzDWs4vRKH7r2Vvd7II9AsfGwhslMCYRXnkfCvbA8D5Q\nvQ1mz4KNtXBWCvyaUSTgLYcncx4KxUB+QboveRA31HwDkgGxaZDssU49erb59/x/geE2le4NL5tp\nj37cphoVE0PcnLntP38UEfEzJKXUe8DXIvKhiDwKPNyd7e9ZvJifH3ww/A011ME5V0J8gE4PSxab\n6+M+iCONdLwNHdy4ORjTg3Y1BsW4iBHFhKWwYTuszYWkmVBcBE99AH99Cz7bDfOHwwfDoBcx1LSY\n/SSRTAWt93gSvXyCeSEusFcG9mwAj7wGA4d0XC6EGC4XQ6dO7dY2u0pPy0ePUpoFi9p53PLXobTt\nAM9GPikMZiCHkMYQpjKFoQzhH6zhWXYyPMbNiMIk7twKz6yH4kx4/XH4v6tg75uwNxl+KoHSElhW\nAw/8HEOjc99KxWLWNbvTaq2MbDyMQQWIA+rXg6PcDBnTGqV8zwB9YB07LqBy0UDEKyRghoisb3G9\nSyk1vTsadjudpAwbxrF33BH+xnJ2QGYQvqYK50Oa71lbHs9Q22pD9Ts2Ni/XLaKG39GHV/Ph6+Og\ndxI8dTrEV8KKb4S8EuHbTKi8DlwGjNtiGjbYkOY4L2mkU0udVxtW0qnlW9/9tZdDSRAuhCpKoToI\nBRYCNr3+uv8ouJFLj8lHj9N/HJz3KDTWwGc+nPFmXA799qUbNOKinEGcTzpHEkdv7NhIIhErFvqT\nDiRwDdlYrW6uHQkptXDdlTB/Gaz4CbJyYVstDE6EbfGwpALOHglJLXTHUNLJJYcSdpLPRgwaaSAb\ngFimo+gNlkQYPAf6z4Q+R4b1Y4omIlohKaVmAK13wquAMAdXMakvKiLzfT8+2kKOgtODWKL6xXc+\nfWAZOOnDUfRlVnNaDXXMZhYWz797PfUMIZYz+sOkFHhjPZx+JMRYhZJdBpk/wMKr4cz34bNtUOAZ\n5O3F0Wwol0QSKXjP0GIZQix+4hg5yiDtkMCfb80PsHtH4OW7iNvhYObjj4cnjEiY6Gn5CBf2i89F\ninz7ZfRJQm84YU67RRq4CRsrvCzsytnLjyyginpyKOEY+pBrQIwtgxtXWiiyw1/GQOV6qNgCb77u\nJs1i8Omv4N2fhQ0VDmb2h6n9923lVmFjCbnUUosdJ73IwEEp1Z5jEzEcjTIAt5/BlgjYAwi93tgI\n9gCPUUQRQSkkpdRDSqmvlFIPdtMoLNVHWjnmWnnY6T18ePfMjgCyt5izgkDY+w18Od5nlp1cDLy/\n0KvYziL2GRQ0AhNIYMQy+N9m+GQX/PJ6qNvthIYGqBO25cJ5o+HgNBi7HhxuSCGeBs+yXRGlrMZ7\nRueiHMGPJVD/k6Cv75DqPinfCwO6b8lu7QsvkPnpp93WXojoUfkIF/HvfIga5Dv0g1/SWzxy8UIo\n8v5fxnMdyZyChcPI4RnsVJPBAI5iJlvI5hXW0YhBH4ubZxLSefFk+MUXkJwOxz8BxSXFkG6jyALb\nc6G0FE5ONFhvg8/y4JFtdmbzPQnEcRnHMIFRFFNMKkNIYBiDaOG1v+4xqPAR4hygZDMsvbvj5/30\nHVjUdk8p2unMDOlFEZkLrFFKnaeUOlcpFa6IaulhqrdDHHV1fHDWWd3XYFE+DA8gRDlA/1NgVrbP\nLAtJ9Grhtw4gnRSmYvrMciMUeNRK9Snwt8mQdwMcMgJiLG7EsGFRNt5ZCc9/DJmVMDXZPMs3nVQS\nPV+Z8YzhOI7yaieGoVjI8N3nulwoCOIHf+KxkNHWJUu4OOq66zjoN7/puGBk0WPy0aOs/MScSfgj\ndTKkHgWlrzQnWTFDP6RyDKnMYhfmGSYbjaQRzy1MZyZ9uYURPKlqOSld8fYMGJcEhT8avDJfQa84\n6Cdc8ieDn6sVSVkJPD0MiixOXk2t4EYOIx4LvUkgjX4czXQc1FKHeb7IjiccRO+bod9XbftdsRkG\nHA6nPdrxZ3DhH+BXoTH62fP552zobEiQEBOUQhKRW0Rkvud9tYjMF5EPgaPDNGPy5UfEzy9eaInr\n1Ytz5s/vjqZMJhwJgTpszX0Lsnxb4zkooLHVKo6BmxiPFZAD4URSqHUqBv0EcxfDD0Ww/meoK7MA\niuqKRlatgbMmwpRBMDgWat1wN/lUecxb91JGQSv3QTa2YuAnImb1BkgMQsG8cH/g57G6SHVODguu\nuKJb2goxPSYfwdCQm0tNKEIf2D3ncYqyTAs0fyQMgIRB5nvnVvPVglTGMI7zqaWSWIRDmEAisXxN\nCcU08jjpxFpgZA089QzEVIIroxcUOWBzI4MGWlj8e3j2LMiIgSuGxPC/Qb05jHRUq1heDRRTiWnl\nV8CluPGzHOeywY4Xfec1UVkKxXntFnGVlFBy3XXN17bsbGxZWc3XFStXsvWuu7zuGXbmmRw+p/0l\nz+4iJGbfIrJEKXUV8E0o6mtBFb6XJfyesLyrxYc9bdo0pk2b1qmGVz32GBOv7aaopXsL4JPX4OQz\nAitfvRmG+Y7Qamcn8a2iuu5iD1M8cYyqcNGPWJSCM/vBkQ4YlQZL3oWMNAHqwNKXS8+A6RNhTB/Y\nVGSevVjSIhZSlWeNvCWCkNDCg7gX7gZIDdADBcCLCwMv20X6jBzJWW+80aU6li5dytKlS0PTocAJ\nSj5CJRvB0rh3L46KCnofYu4h1m/eTNGTTzL2pZcAqHjoIeKPPJLkWbP8V7JrK7z/LPzjWTjnJt9l\nlpwHxz4NyYPN/dV+fwDHCny5CMnhZ2qpwk1vIIYMMhhIAgYW/o9S8jMHUjkfHr8W6goV8zOtUJ+D\nGjyEbzMhNQHeK4fleZBbrjj+hHxiGchrFPMUEyiniJ2sYwKHMwwzbtNofsbGA8TxO6yMhMzzYOTT\nEDcYYhLhuKd8P1d5IaQOMN2K1VTAwGF+P6aY/v1Jv/325mtHfj6IkDjWXH1JO/poUg7x3stVSqHC\nGL0gKNnwF5fC1wuYgzkq+wq4CujdIu/vwdQVRJvlra7fA6b7KRt8cA4/ZC9YIG5PXKCwU18nkp/T\ncbkmSlaINBT7zCqSJ8Qm3nFhvpHlUi21bcquqhG5/weRXn8RiTldpN/QaoGfpN+QWrn+RZGkO0Uu\n/EDk0q0iNnHLCbKl+d4s2SObZKdXfcXyP3FIue8+Fy0WsfnucxvsjSJXTA+sbAj46NxzQ14n3RQP\nKVD5CKVstEfltm2y9IorgrrHcDh8xz1qqBd58t6OK3DaRRoqOyzmkFxxt5KD3ZIj22S7iIi4xJAt\nUiPrHS45e5sha7NFBlwvctLlLmFYpaAqJDGtTM66ReT374gszBI5bp1I0jIRl1ukQZzyuhTKu7Jc\ntsoeyZMtsk0WS4Pkyk65z9NGlhjiiV9mGCLOOpHCT9rv+Bu3ieT7iPUUpbQnG8HuIVUBo4D5wIVA\nlVKqXClVHmQ9wbBYKdUy1OkoEQn1TMyLtS++SF1JSXhjHrVkx0b4PogZQW0m2H0bQMSSTqPHFUoT\nGaQRhxlK4nPKOJHVrK+Dk7ZASQOcfiy4KqA0PhnUaCqqkumXBBW3w92nQboV4lCcxT5XPhvIol+r\nLQw7e9osWTRTvQacAfqms1jh8e6yboQZT/kZmUYH3S4f7ZF60EGc9PLLQd2jYmNRFh8/RfEJMLWt\nY9Rmdv8MjXWQ9yOsf7XdNmz8k3q+wsFur/SRjGA5Obhw04CbZVTyZKbi+1wYOgCUC/b0sWIpMmB8\nMqddmsaGeXD1CXD6GHj3IPhlGjypdqNQXMYgzmUK4xnMUA5hOIcSQwajMGd1VsagPLKIUub5PGdV\n+x/QpffBEN9GTACuN17ByNvTfh3Rgj9N5esFzADObZU2CugTTD1BttkHeBA41/N3YjtlQ6LBDbdb\nDMMISV0BkbdbZOfWwMuv/LNIle/yZfKa1Mi3XmnPyavyrfwoIiKGGFIgjc15lTaRd7eLfLBEpO9o\np1hjXZIx0i02m8hJL4i8UyjycoHIGqmXJVLTfN+nskxKpKL52hBDdspN/vu8bJY5IgyEVd+J/O/h\nwMp2kZ/uu08ayv3M6roA3TdDCkg+QiUbEcW3/xEp3xNQUaesF0PcUiPrpFDe8sorl2q/9z28UMR6\npciYE5xy5hUV8o9/iRw7S+Q/q0SmfSFSZTdnR1ulRlxifr+Xy8+yxTPrypUvpUK2Bf5M9joRhy3w\n8iLiWvGjGFVVQd3Tk7QnG0HtIYm5V9RHKTVdPKMwEdnd0X1dQUSqgSbfGB+Gsy1PexSuWsWQ7jyx\nX1YMdTUw9uDAysekgtVHFErASnwbbwm/5jSSPC5+7AhvUsTNjGTARvgoHd7aBp+eA3tvj+H/5kLN\nQTDhKsh9E+7JhhEJ5t5T0zayAxfV2OhHWnMbBg248XHiHKAuG5KH+fW714bPXodrusfcfvjMmcS3\n52E6wulu+YgopvlwUtoScQMWcK4jJm4yAMkcREKrs3LpeP//Z5HPXNI5iSRumgkHxcCqNTE8uyAN\n+wDh0F/BQ5sUg3fA1In1jCOGQYNzOJo0LiCNo5hEnOendTitHMAC2BeAUQRxF4C11Xdv8zxITIOD\nA7f4tB4bxHGKCCdos28xret6bEkg3DgbGtj52Wfd2+jEY+HEtjGN/DLuSkga6jMrmRNJwHvTsoBC\nSjE9HidgoZAG6nHxwGAY0Q/cBjz/DayvhJN+B3f+DvbEwYZCc8ni90PgXWoY5gk6pqDZY0MTRbxF\nPy723V9rAgz9XeDP948nYPCIwMt3kh/vvBOL1YqlKzF4NJGHrRD2vAFZF4M9CxrmNWdZSCCGPizH\nf5iRU+v7c82eOK5YBVtr4f1CyHHAoMGwpkgo3Gjw3hnQexY8PUjxxCDFGaTyPbtYQzZ5lPmtG4C4\nk0EdDfl/b5s3+Y9BKaP9jYj21NATVGZnMyVEURcDJnML3Hldx+Wa2PQI5PoeDFfwDCV4zy7iiaWq\nhYXwbYwhCSsL6+DnBhiUDHkVEJsAV54JN/8WBg+FpbXwaI6pfKwoJnhmXj+yhbGtLPniGEisv2Mx\nDbsC3z+qq4GXu8F3IDDlppvoP2lSxwU13YprxXKkrq7jgv6wJkHiUBj3HiSMg1Rvf3cKK8dwF9/w\nJKXks4DPAbiLNQDclBzHpmFWEuJg+gpYvwpcGZCtoLGPhQ3jrDy6Ek6Y7GCh3U2muBlLKkNJwkkC\nYxhIo48D4gZ7MSgBlQwJR8CIdsy8186H/I2d/wyiFK2QWlGRlUV9SUn3NrrxZxh/WMflmhh6OqRP\n9JnVn/sYyGNeacMYQkOL80F3k8lL7OHafvB+A5w6CnoNgY93Qp4Tqhuh4B64cQJUnwbbcDCM+GaD\nhRM5jONazcJKWUQvfJh1u22w69+Bm3zn7oAppwRWtgsULV+Oxd9muqZHMbZtQWp8LP+WFno8fFfu\nixHki7hU6Nf+dyiGBKZyKTtYw3jG8xMFHEIfSjyK5C2bQUZ/N6N/goI0OHkMzLkQ+oyBeReAXUFs\njJt1MTUUqVqWU8s9TGIao6mghCV8SCY/4MZGDZsoYg71fI5L1oI7gLhg/cdBip9zex+8Ct/upxGV\n/W0uReOLEGzc2io7Nh/tcWylIkXLfGY5JFfy5WyvNLe45VV5ufm6WhyyQ2qk0Ckyz7Ofe+9qkSVZ\nIofcITJ3vsivl4kc/JWZN1lyJMdjqpojJXK1POdVf4WslHL5yXdfHVUiu58P/Nm+mieSHYSBRyfZ\n8OyzUldUFLb66SajhkBfoZCNHueRG0T25ot88azIuq+98z6+VqRgrciWp/3e7jLWimE0tkmvlFq5\nRebLfMmUFVIiIiKGYYjbMCSvUKTf3SK3fSRyx3KR59aL/PlLkY+zRBxiyBoxDRC+kmJZ1eLIQ6PU\nySZ5W8plgxTIB/sac5eK1N5gvndUiuS8YL6vC8w4Q0REqipE6toe44gW2pONHheUUL66KnSGYcib\ns2Z1qY5OYWsQmfdyx+WaqNgqsulJn1lucUittD3X8Jl8KC5xiojIbqmTtyRXREQG7zJkTqkhB78j\nsjDHNIQrbhB5cIvIV0UidnFLubR/HitHXpACec9P5jMiBW8H9ly2epGlHZzJCAFOW3BWTJ1BK6Ru\npsmCc+ebfvMb3f+WWuMtKZd3JEfub85yiVsKpczz3pBdUi03SLZUilOyy0QG/lXknWUi/a8X2VEm\n8uqGttWXSqNUi0Pc4pZNnsFZhewQewvL1Da4GkT2LjTff3N+2/yGGpHNSzt89GijPdnQ6xUtUEpx\nyYIemAori3n2JlDEDX7O+1iIpY4PcON9NCyNdEoxlyJHkkw+DVThYOMweDBDsfW38MBmOOJjGP0p\nTE6FXw6EseSzA9OrcB12/oW3wUcDBSRxMIPx7TkCexn0OzOw5yrYDdUdbAiHgC9+8xtcje0s+WgC\npscDGjbWgr1hnwXn2Et8lys5nnjL34hXJ9CbUxjBrQgGglBGPoM8HpfqcfIdBVzFID6kjtEZ8N41\nMCoDjhsLNy+A8w+CDxtdFBpGs3FPX+J5BdM/nTSHaBlPHD5iljlLoepjsCZC/9PNtFPehx2fgrOh\nxbPVQ7Fvn5WB4KoNwGt4hKEVUiRgq4dvPg+8fJ/xULjCb3Yik3DiHSVzBCNYwL7DpkfQm+8pISNG\nMbDAyd1VLh6bCo8eA/UXwS8Hw0vuOnYzjOM8kWcTieVwvK3f6siknB99d6Q+E6xxEJMc2HN9/ioc\nE97ICavuvptZ8+YRkxBgIERNu7x/Svj3+9pl7Uew5Uv/+fUbIPtaGGDKSywjiGEAAOt4gK+5h40s\nbS7emziu4GDGEs8pnu99ehxU2ODjG+HNCyE+TlDK4D0qWEIta6n0GP6UY8HC4ZgRbmvYQLkvb2pK\n4fOnt7ECjBYKPm0gzAg8PHnR+ec3v3fX1bHz978P+N6Iwd/UKRpf7O/LEi3Z/ZmIo85nlk1WS5U8\n7pXmFqdsktXiEHtz2nzJkUZxSqPbLZ/aXLLLuc99i80wZISjSKoNM61aGuQmmdemrc3yoLg9S4Ft\nqPpZpHJFsE8WVop/+qlbXEKhl+yCwnA6pfaLL3xnXnOeSG07S18dVt72QHatbJZS2dfeMvletshm\nrzJrpEHelUpZUS9y69596bcYFXK1USoiIuulWM6W70REZKF8I9tknaf+IrFLmVTLCqmV5eKW8O9N\nO/Lywt5GKGhPNvQMKVqpyzNHVD6I53AUMV5xkSzEUE8NeS38bv5EKdU4ibdY+LTRzSW1Nj62O7jH\nVseExnJyYgfS2xMEcBOFXMkvvNop4muSGIbF1/lqwwGFr0LiyMCe557fQ37nlycCZcDUqd3nEkoT\nMOJ04sjM9J35yMvQy8fSV0dkPg6Fn7Y5kO2iGCeb6dXCWfAJnMgEDia/hQf7ccTxC5I5Ngnu7+85\nNG+4+KtK4QXVF4CxpHG+x5P+DE4gny3sZQeZfEwcGSQwGKGc+lZHMbDvDfw5LjoI3N7ezcXpRFot\nycUO9X02MZrQCilSeOweuP7Sjss1kTwctvv2UK2Iw80u3HhH3RzDIexme/P1v5jCW2SxFxsvpcbx\nq3grf6118LNhkJvYt7ncRgrJopRxnqWOJvbwBcPwc4hvz7+h1wSIHxjY84yYAEmd+NEJkOIVK1h+\nkx8v0Zoex5KYSPr//Z/vzJQgPGn88O9978f+GQa1jWlmyHriGezlsUGhqMfG+hYDthSsDCKGNc3u\n6QAAHadJREFUHVTRYBjc6axjulHG45RSjot32MRuKtlJAYtZRixxHMU0+jCYyZheJOIYTgpnkcLT\n3p3YdL05aGvNJ7e1Tbv2wTahaYwfvsf16ksBfCBRhr+pUzS+iPBliXYp3StSmB/cPVteFTF8eEkW\nEbusFJssbpP+hbwheySz+fpF2SZl4t/qzClu+ZcskdJWXpJz5EPJly/99y37n3771obGhsDKdZKq\nnTul6IcffHuUDhPoJbs22MvLpWb7dq80V0OI//eb3jf/Nvrw7ebYIUbZX8Xmvq3dKqrFJrlSIaVS\nL/lSK1fKEjmxvlJea+VjziEuMTz+63JllzjFKdWSKwXyozR6zMebMBoeEsOZLZJ9etsG178isucH\n8/22tjLbzDd+ljQ7wFFZKbWZmR2Ws9fUyKKrr+5UG8HQnmzoGVKk0Lc/DAoyZHdjObidPrOsDMfJ\nqjYufo5nFnUtfM7N5iDyqOWeFiHOm7Dh5Fo+5o9MpS+9mtPL2UgB3zPEl58ugOKXIeVw03owEOae\nAy7fzxEKxO3GcDr1IdgepiE3l8qf933P6rdvJzvUs9bDPBv7P8+F6hZLwJVfQsWPqNRbSbDc53XL\neh5mBU9TwCb2kksVNvKo5EU2MZhk/st0liWlcnlsAoKwgWIAYrHi9ASsLGQPdhqxEouDUipZ7d0v\nowCohRHvtO3zyBnQz7N8eNAM/8+2bX2bpbtAqN+1i8of/RgetSAuJYVjbvHvUqlb8KepovFFBIwC\nu8QTD4hs2Rh4+ZUPiRT4NxqoktvEKW1HRi/JHVIs+w7ibZBSuUdWia2FccIKyZN7Zanki/dI0ykN\nslPmSbVk+27Uvldk05ki7rYHEHsCt9OPwUWYQc+QuoxhGGJ/+422GXnZIu/9J7jKXLUiDt9e3V1i\nE5c4pFTyZZf4lr+tUiqrpFBc4pZ3ZZOIiDjFJc/Lu7777s/QpwlbqchW/4d4w0ljZaXYSko6Lhgm\n2pMNPWSMJEaNhWDOdRzzD8j+xG92L2Zj47k26X/kLpw0YPe4STmCvtzBFB5jDR+ygz/zJWU0cD1T\nGdIiBpKBkxLWE0MvejPad6Nlb8OEN8AS7zu/Jbs2w/nh3YhdMGMGYhhhbUMTJkTAl0+79P5wuMcb\n/7PXQdZa8/3CR/2HNrf2gljfvhbLmI+DPGrYjcK7vW/Yjh0XuVSTTSUKuMhjDBGDlau5kJ2sYDvL\nvO7L4VLqecbLsMiLmF6QfmTb9I1LoKz9MOVdpfSnnyj8+uuwttFZlKmw9g+UUhLVz7NmJaSmwRj/\nwbjakLcUhk3zm13LH0nmcSwtFAvA5zzPOCYzgWO80iuwsZoifulD4exlNdt4k2k84bux0nnQsBGG\n/RMscYH1v6EOknp1XC5IDM8yXU+dN1JKISIBxtsIPz0lG/WFhSQPHtw9ja39BCadHXiYE0Bw0sBW\nEhiPOWxRWImjnEI2sIVNNHI9Z2HHxRYK2c5extCbwaQwjAHM4xUu5kpAsDYF3vNg5zviOH5fQL7q\nNdCYBwPO8d+hdV/BkAnQf2SQDx89tCcbeoYUSRw1NThlBIAVvv6L39xEZtNIW+/Zs5hNGbn8xEde\n6ekk+lRGteRTzGr/yshVA8qAobcGpow2/gA714dFGQEULV7MtuiOBBtxrHvlFTbPm9dxwRb8/Pe/\n47bbw9SjVkz+9T5lVJUHeW33RZupvANchTTKOpz8iJVEYkmkio3YaSCLtfSlP8NIIRYrFhR27FzG\nVMYzgGQSiSWOS7mGIlZip4pc3qCU75qbiOdkUxlVvAK1iyF+MCSNMzOdNsj3cbh90mkRpYyyFi7E\nVlnZbe3pQDCRiGFAoBvwyQOhr//AfjEch4P5uFhNDFOa061YOZgTKSGbEnbTn1E+73dhJ4svyOVb\nTmttutqECFR+DCoWrAF4ZbB7XPPHhW/2MuS00xhymh+jC02nOOy3vw3aMOSUt98OU286wFEH9e24\noep1GWJJB+clpMZ925xcxVb6cQzHcTaCwUTPmD2JOFzUsoVsDmcsVRSxh00M53DiSMFKHCO4zHdb\nvc8ClQDWFEgYZKbZq6F4DQw9LlRPHBYcNTUYzvAZHLVGL9lFGmt+hj+cBxuDWEd22aFoFQw70We2\nIDS4Z5BsbevGpJw8VvExR3MOGQzzyitlOw2UY6eKkUwjDj/KJm8OGE4Y4Wf21Jql74HTDjP9CHAX\n+eK44zhzhX/XSt2BXrILE7m74O3/wtwHgrvP7TDdWC0/F8ZdAvYcGOrfwq+ANVSSy2GciwsHLlwU\nUUQ88QxmKLWUUU8lAzFnPIITFxXEtjqr1yFvz4ZT/gqDDum47H5Ce7KhFVIkIhLUOjgN5bD5NTjm\nb36LOIwPgBLiLG0DAdZSznJeYySTGcnRxBLPz7xCb/ph4OIwzm9bYRO2HRA3NLCZEcDmH8EaAwd3\nY4j4HkArpDDhdsPeIjOCZKCUb4bN/4GTn4XGYkgYCGL4PJZgYKOSBWRwXnPacuZTQwmTOIdYYkmn\nL8WsZyATcdHABm7lUP5MNV8xgOsRYxeowSgVwAqA4Q7OsTIghkHd66+TcsUVQd0XKeg9pGgjGGUE\nkJRhjgAzP/ZbJFb9BnFvxJDcNnkpZHACV9BAJV/yIIt4ACc2JnBG+8rIVQq7fwOuINygJPcxX2Fg\nwx13ULp8eVjq1kQIVmtwyggg4zBTGYGpjMBbGRmmJZypjO7H0uLMnYt6pnI24xjDAAaRTl8EodTj\n8SSGJI7iCeLoT5wnTZwfQuPXIC6o/BHK2rFoC0AZiduNZG73TuzEeaSowJ89eDS+iMKzFv5wnnua\nGI1BnuVZ9z+Rxmq/2YZRIY2OK8TwOIbsEvY9IkV3BXfP6/8UWbOo6237oD4/X5z19WL4cKTZE6DP\nIYWExq1BBmts+v/v+F7k2xcCu6f4UhFnsYiI2D1njJpYI7eLQ+olR75vTnNKvc9qXNLCA+veu0Uc\nBSL1u0Xqdnj3zR8LXhPx4fjXKC4S173te5eIJtqTDT1DilAst9+PbN4Q3E1J/cBR7zdbqTRi1EW4\nnQ93rXPOfNh7J6RdHPg9Hz0GhxwHk0MfXqKxpITN999PTFISKtjZpSZiMRobqQjWUvI/10P2Ohgx\n2TQBD4SM5yDG3PuJldGI7Ds7NJl7iSWJER7HwjYKyeQZXFSxh5twsAYXOYhUYHHn76uz/50QO9i0\nrLN7ZHL5A5D1hf9+NHld2bAS3tp3flANGIj19vv83LR/ETV7SEqph4BFIuIjwEhzGYmW5+kI2bIJ\naajHcvSxQdwk8M4s+F078WEAl/NVlDKwxgQea6UZ2xoovhmGvmEKXCAsfBFGHQljJ0NMbMflg6A2\nK4uUsWNDWmco6M49pANNNnzy7cdQmAOX+HHQ2lABa9+AE2/0Tnc3wPaL4VDzgLnL/RFCDbHW3+Om\nDHDjxkUV39C/lRWdm1pcrMFCf2LcMYh7JZa4VoY6tmLY9Roc+o+2fVr8Iky/qq1Fra0BqithYJCu\nxKKEqDZqUErNACYDVwPXHNBCFwiOelh6D/zS/yxIxI3h3oU0XkZMr58Cr7v6DbCtgAHPBrfPlb8D\n0geHxZv3qtmzmfz001gjLOBedyik/Vk2HKWlxPXrF/gNhmG+YvycZHHZYc9KGH2Sd3r+R9D7UOjt\nOf/nLgWsYE3HxrcINhKYTiO7ScL/8QovSj+Hfm29jLfhu9fhpMuC3zOOcqLaqEFElojII8Dunu5L\nT+C8/y6kpCTwG2KT4NALoNj/cp9SVizWUVji/oIYTsToIGy4uKHofHBXwoBnghOgu06HoeELLXH0\nSy9FnDLqLvZX2XA3NLDzL/4Pe/vEYvGvjABi4vcpI7cDPvS8TxwMsS3CW9iWQqN5ZCCRU0jiDCwk\nkMTB2PmBenyEfDEava9rW8heZRZU+onzdfLlnVZGuX/7G43Z4Y8f1t1EvEI60LGceBJG5raOCzah\nFGSMh00+vAp7FYvBEncJ2J+GuktBbKYpbEvclVD9byg8Ffo+A+k3BO7Bu6YMFj4Pc4I72R8oS37x\nC4xg/P5pogZrUhKHvvtuu2VqZk7vQgNxcO735vuMqZDYImZXrwsg+Uyv4m7XWwDEMpEEj4d7l+Mm\nxNhjFsg/y7Soa2J0i5hGtfnm68eXwNlKcXWBwXPnEj/ajz/JKEYrpAhHDR0OwbruSOgNU2+Et9sx\n2W6qP/FvqN5fQt3tUHcD2H+Ahreg+l4ouxZixsLgryE2wEB7ACW58Om/Ycj4kJt4N5aUIG4307//\nHkt7I2LNfk3Kp+0YB4SChizzrzMfXGbQPgu9sNIfAKvlj6A8B8mHLwbl57s4fBoMPxmS0oI+b9Qe\nsf36NRvwZN57L/ayDlY5ogStkCIcy5ixWH/VjjNGf/QeBOe8CMseCax8r3sh8WpwF4ErB+KOhgHv\nQvKv/QubL8ryTCV06mw4ogujWD9kP/UUtZmZ2pruAEclJnonFOTBzdf4LlyYbcbbqiyEFy7vuHLD\nBTkeqzbHVqyuY9oUUQW3osQJdQGuXkw6H6wegx4RuCcIC9UOGHDWWcSmpoasvp4k4o0amlBKLQIe\n2p82boOh8e//R+zlf8B6hA+X9f5wO2HPcojrDUMmha9zTTTWw3dvwKhJMH7/9sTQEd1sZbffyMae\nd9+l74knkjS0E2FJKisgzUeIibfugRmXwcBRpgl2vMerSFUeLP8PnPEArHoAJt4IsQF6HAHTQi9z\nDgy6FnofbM6AKrMgIQ0SM9q/tygHBo0MvK39iPZko9vXPJRSs4GjAF/SoTzpD4tITmfqv+uuu5rf\nT5s2jWnTpnWmmogj/v6HMXJ2Iy4XKtClKmssjDoZPr4OhrSNixRynvgt/O09iEvsuGyQLDv1VE5Y\nsABLXIBhLbqZpUuXsnTp0i7VoWUDEgYOxJqU1LmbfSkjgEvu3Pc+voXC6T0Yps423w86HqztGMc4\nSqH8SxjUwqzbmgQHPwsb/gnjr4PEAaZPybSx3grpndvg7DmQ3GIW00oZSV0dsnghlnMuaP8Zo5Cg\nZMPfidlIewGLgOkdlAn4tHA0Yrv3n+IuKOjczW9fLFLmJ8prV8nfJvKfq8JSdeNe8/R7T0V+7Sx0\no6cGLRthomqryE9XmO+d1SJlX3aunt3rRNxu77TcLJHn7m++NGprxf3+283XpTNnirumpnPtRTjt\nyYbeQ4oiEm6/C0tng51d/DZs+QRyQuzrbeETsGcDzP5PaOv1sG3uXOxlZdqAQdMG28aNGIHEWpp7\nKRTkBN9A/ECY+or5PqY3ZPgIZ1LwA1S1ML9+83IozfIuM3Ji28OvA4fBGRc2X6pevbCcv29fqe+i\nRVhSwnNUIpKJ+D0kpdQk4CJgDrAWmCcij/opK5H+PF1FnE5qzz6T3gsXBX9zXSkUrjc9DB90etc6\n8tUzMOIIOOikjst2guJPPiHt+OOJD+ZwZATRTQdjD2jZKHnkEVJ/+1vihg1rv6DbbTplDYb6Ylhx\nK5z6cvvl9nwDvYZA+gTzOphYZgcoUe2pIRj2R6Hzhdhsba2MAqUqDyr2QHUBTLqw4/L++OJRGHII\nTDyj83W0g1ZIoeVAkY1Os3OFaV034Rew5EEYPBEOnuVdZveXMPAYSPSzV9UDlK1aRc78+Ux56KGe\n7krARLWnBk1bVGIitkcexrFwQfA3pw6DAQdBcgZkL4OiLYHdZ7ghfzM8fjbUV8GZfw+5MjKcTqpX\nrwZg4K9/HbXKSNM+jro6tr7T/sHtbkEEPnvefN8rHXp5DBGmzfG9guCsNQNRBsrbj8COtV3vZztk\nTJnCxH/+0yvtpzvuoChKw7BohRSlJPz1JmJPn9VxQV8kZ8D4GaYVXu5q2Nq+M1bWfAQbF8LKd+Ev\n872thUKIq7qaskWdWIrURBdKYQl2CS1InN9/h7Nl1OCli8zAfq2J91j0DZpgzvjBDCDp65zb+Atg\n+f/apht+YhP98hIYdWhwHQeK/v1v6tesCaisUoqYVqslR82dy8DjIjs0uj/0kt1+gDQ2ojrrz60y\nzxSobx8zlyiGTYbYREjsDfcfDRc8BqmDTXPZPkF4awgCe2kpWXffzaHPPBOW+nsCvWTXs7i3boXY\nWKzjzBDjfPUZHD4p+OB+rdn8KRzWIqxF3mb46mm46oWu1dsCR2EhMenpWPZTH416D2k/xlVQQOWx\nx9IvL69rFTlssHmBOTLc9Dlc/jLkroWhR5gjxjBjy8sjsaPN6ShCK6SusfCcczj9o4869MjhbjT9\nw4XVwe66BTDkYOg/KnxtHEBohbSf4y4qwjpoUE93I2jynn2W2PR0Bl4cOjcqkYJWSF2joaSEpP79\nOyxX9Jbp+HTQJZeErzNbl8KAsZARxOzq6lnwyNvQJ828vvgseO51/4d3DyC0QtJEFDWrV5M0YQLW\n5GRQar/0S6cVUpTx5Fw47SI4aKJ3+vqlcOTJvveU6mvgn+fDowHse4qAUhgNDVg664liP0Fb2R0g\niMtF5R/+0NPd6JCatWtxlpejLJb9UhlpopA/3gLjj2ibvvE7cDp835PcOzBlBM0KrfLyyzGC9d5/\nAKFnSPsZrqwsYsaORQwDpQ/o9Rh6hqQJlooffqBmzRpG3nhjx4WjGD1DOoCIGTsW+3ffUXPbbR0X\n7kbWT5+OGEbHBTWaSMYwwBGAu6JO0GfKFAaHcy8sCtAzpP0YV2EhuN3E9KD12q7bbiPt1FPpfdxx\nB1SocT1D6j5EBFt+Pkmh+J6/9zqMGA1TT/Sd/+Mi2LYOrvpH19s6QNEzpAMU55YtODZu7PZ2Dbud\n8k8+AWDUffeRdsopB5Qy0nQvjcXFbH/Upwu/NmRPbxs0suHJJ/fN3o+fBgcd5r+CE37prYyWfQvv\nvxVEbzXtoRXSfkzizJkknXkmRkMDdR991G3tisuFs7gYQBstaMJO4qBBTH7yyYDKjvr88zZpliFD\n9lnRDR0OfYLwRHL4RDjxlHaLNGzYQMHcuYHXeQCjl+wOAIzaWurmz6f3FVeEtZ2dV13FkDlzSJow\nIaztRAN6yU7ThBgGRm0t1j59erorEYE+h6RppvLJJ0maOZP4Qw4JSX22zEzqfvqJfpdfHlw02/0c\nrZAODGxvvQUuF4m//31PdyVq0ApJ04xj505ihgzBXV1tOmYc2Dn/dNXffEOf6dNxVVZiz8khedKk\nEPc0utEK6cBAHA4QQcXH93RXogZt1KBpJm7cOCxJSTSuXk3j2uBc4zdmmZEwxTCo+fZbAGLS0rQy\n0hywqLi4Tikjl83G+ltvDUOPohutkA5QUn71K3qdYcYzyj31VIyGBsTppGkU3fLMUP7f/46rooKS\nRx/FXVeHslgYdu+9PdJvjSaUNBYU4Kqt7fZ2rQkJDD3nnIDLV+zcGcbeRA5aIWkYvmgRlqQk9t5z\nD/XLlgGw+7TTMOrrAeh71VVY09IY/vzzWHv16smuajQhpXzBAuo2bWq3jLjdlDz1VEjbVUrR95hj\nAi7/wwEyANR7SBpNGNB7SPsPIkLVBx+QdsEFfssYLheOvXtJGDKkG3sWneg9JI1Go+kkSql2lRFA\nQ1YWeS+91E092n/RMySNJgzoGZJG4xs9Q9JoNJoIYPMLL1CxbVtPdyNi0QpJo9EcMJRt2ED2Bx/0\nWPuDTzqJXkODiDx7gBHxx+qVUn2Aqz2XU4CHRGRdD3ZJo4kItGwET9LAgT0aJyz94IN7rO1oIOL3\nkJRSz4vInzzvRwFrgMkikuOjrF4n10QE3bGHpGWje6hYt46qjRsZrd0DhYSo3UPyCFl207WI7AZ2\nAef3WKc0mghAywZU7trFh90Q0C555Ej6Hndc2NvRRLhCAlKBh3ykZ3R3RzSaCOOAl4200aM5963g\nYhFtmzcv6Hbi09LoPX580PdpgieiFZJnPfyoVsmTgUU90B2NJmLQstE56ouKeroLmnaI+D2kliil\nrgbOE5HT/OTrdXJNRNDd55C0bGiihfZkI+Kt7JpQSqXSjsA1cddddzW/nzZtGtOmTQtvxzQaYOnS\npSxdurRH2tayoYlkgpGNbp8hKaVmYy41+GpYedIfbm0ppJR6HrhZRGraqVuPAjURQWdmSFo2NAcC\nUR+gTyk1B3i/SRCVUpN8nbfQQqeJFLpryU7LhibaiFqzbwCl1HnAWqBSKdVHKTWZtpu5Gs0Bh5YN\nzf5GRO8hec5avM++JYymZYuZPdYpjSYC0LKh2R+JiiW7QNHLEppIQXv71mh8E9VLdhqNRqM5MNAK\nSaPRaDQRgVZIGo1Go4kItELSaDQaTUSgFZJGo9FoIgKtkDQajUYTERywCincfsfCWX+01h3u+qO5\n75FENH+Ouu/dX3co69cKKQrrj9a6w11/NPc9kojmz1H3vfvrDmX9B6xC0mg0Gk1koRWSRqPRaCKC\n/c51UE/3QaNpItJcB/V0HzSaJqI6/IRGo9Fo9n/0kp1Go9FoIgKtkDTdglLqIaXU9J7uh0YTaWjZ\n2EdEx0PSdC+e6KPZwGhgia/Io52ocwYwGTgPWNTV+vy00Qe42nM5BXgoFH1vVXcVZqyhF0RkSSjq\n1kQPWjbarTt0siEi+rUfvICHgOlduP89YGKL60Uh7t+irvSvg7qfb/F+FFABjAzV59qqbgPo3dP/\nb/0K7n+oZSM6ZEMv2WFqeqXUHKXUbKXUe56RS6jrnqOUmqeUmhSquj31z/CM3s7rYlUzRGR9i+td\n0bCM4Imcmt10LSK7gV3A+SFqYnbT5+CpG8xR8gGBlg1Ay4Y/Qi4besnOZK6I3AKglFoMZCulUkWk\nJgR1Pywif/LUPQpYo5SaLCI5IagbMafIS5RSnQ5d7fmR2dUquWka/k0XutcdpGKOgB9plZ4RovqP\navpfKaVGY4YJb/1Z7c9o2dCy4Y+Qy4aeIZmEZRTcDSOUUJHqI62cKJgJiLkeflSr5MmEaE2+1Y/j\n1cDNIfoxjha0bLRFywbhkQ2tkEyOEpFvIOSj4KYRSmtCNUIJFek93YGu0HI5RSl1NfC1iHwbqvqV\nUqM8Sz+jgJdCVW+UoGUjiok22dAKifCNgsM9QgkhFT7SIu2HoUOUUqnAeSJyWijrFZHdIvIIcAuw\nVinVO5T1RzJaNrRstEeoZUMrJA/hGgWHe4QSIqrwvTQRbXslDwEXhLJCj2kr0LysVAXMDWUbkY6W\nDS0bvgiHbGijBg+eD/QRz9r2Ws/masj2CgIdoSilZmOOHH35dFKe9IdDtfEL5uavUqr10sRo4PlQ\ntRFuPD+YDzX9z5RSk6SL5y08G9pf03bg5usHar9Fy4aWDR91hkU29kuFFOwXVynVR0SqwRQ+pVST\npm+j7bsgFAGNUETkJXpmn2KxUmpii1HrqKa9g67gMeW9CJgBpCml5onIo12tt1Ub5wFrgUrPqG0M\n5vJPVw8A7gJubpU2CpjTxXp7DC0bnULLRlvCIhsHvHPVJk0vIpYWaauBVSJybYjamAO830LIuzxC\n8dHGIsxRUKcExfNlvQVYBRwNzGt19iIiaWGt1fRFbvrhmxmiH40ZwCSgGlOQvxaRD7tabzSgZaP5\nfi0bvusPuWxohWT+085rOTJRSpUD54diPdszQqkCVnuSxgCTReS/Xa3bU3/TKGsO5kgo5KMszYGJ\nlg1Nd3PAKyQI3yg43CMUjSbcaNnQdCdaIWk0Go0mItBm3xqNRqOJCLRC0mg0Gk1EoBWSRqPRaCIC\nrZA0Go1GExFohaTRaDSaiEArJI1Go9FEBFohaTQajSYi0ApJo9FoNBGBVkgajUajiQi0QtJoNBpN\nRLBfhp/QBI/H0WU6ZviAh4FTMWObVHnc/ms0ByRaNroPrZA0TQL3tYjUKKXGAO+JyNFKqeeBrB7u\nnkbTY2jZ6F60QtIAVLaIADoamAcgIn/quS5pNBGBlo1uRHv71nihlMrCjHcT8QHINJruRMtG+NEK\nSdOMJzJmhYhYe7ovGk0koWWje9BWdpqWzMSMrAk0B2fTaDRaNroFrZAOcJRS5ymlmkJInwpUtMju\n0wNd0mgiAi0b3Y9esjvA8Yz0zge+xhwBXoMZWjpVRB7tyb5pND2Jlo3uRyskjUaj0UQEeslOo9Fo\nNBGBVkgajUajiQi0QtJoNBpNRKAVkkaj0WgiAq2QNBqNRhMRaIWk0Wg0mohAKySNRqPRRARaIWk0\nGo0mItAKSaPRaDQRwf8DkAXLsFfcKucAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "figsize(6,6)\n", "subplot(2,2,1)\n", "scatter(xt[:,0],vt[:,0],c=numpy.fabs(x),s=1.,edgecolors='None')\n", "scatter(xt[:,360],vt[:,360],c=numpy.fabs(x),s=1.,edgecolors='None')\n", "ylabel(r'$v$')\n", "xlim(-3.49,3.49)\n", "ylim(-2.49,2.49)\n", "annotate(r'$t=0\\,\\mathrm{and}\\ t=18$',(0.95,0.95),xycoords='axes fraction',\n", " horizontalalignment='right',verticalalignment='top',size=18.)\n", "subplot(2,2,2)\n", "scatter(xt[:,500],vt[:,500],c=numpy.fabs(x),s=1.,edgecolors='None')\n", "xlim(-3.49,3.49)\n", "ylim(-2.49,2.49)\n", "annotate(r'$t=25$',(0.95,0.95),xycoords='axes fraction',\n", " horizontalalignment='right',verticalalignment='top',size=18.)\n", "subplot(2,2,3)\n", "scatter(xt[:,800],vt[:,800],c=numpy.fabs(x),s=1.,edgecolors='None')\n", "xlabel(r'$x$')\n", "ylabel(r'$v$')\n", "xlim(-3.49,3.49)\n", "ylim(-2.49,2.49)\n", "annotate(r'$t=40$',(0.95,0.95),xycoords='axes fraction',\n", " horizontalalignment='right',verticalalignment='top',size=18.)\n", "subplot(2,2,4)\n", "scatter(xt[:,2640],vt[:,2640],c=numpy.fabs(x),s=1.,edgecolors='None')\n", "xlabel(r'$x$')\n", "xlim(-3.49,3.49)\n", "ylim(-2.49,2.49)\n", "annotate(r'$t=132$',(0.95,0.95),xycoords='axes fraction',\n", " horizontalalignment='right',verticalalignment='top',size=18.)\n", "tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also make a movie of the evolution of the system in phase space:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false, "scrolled": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def init_anim_frame():\n", " line1= plot([],[])\n", " xlabel(r'$x$')\n", " ylabel(r'$v$')\n", " xlim(-3.49,3.49)\n", " ylim(-2.49,2.49)\n", " return (line1[0],)\n", "figsize(6,4)\n", "fig, ax= subplots()\n", "line= ax.scatter(x,v,c=numpy.fabs(x),s=5.,edgecolors='None')\n", "txt= ax.annotate(r'$t=%.0f$' % (0.),\n", " (0.95,0.95),xycoords='axes fraction',\n", " horizontalalignment='right',verticalalignment='top',size=18.)\n", "subsamp= 4\n", "def animate(ii):\n", " line.set_offsets(numpy.array([xt[:,ii*subsamp],vt[:,ii*subsamp]]).T)\n", " txt.set_text(r'$t=%.0f$' % (ii*subsamp/20.))\n", " return (line,)\n", "anim = animation.FuncAnimation(fig,animate,init_func=init_anim_frame,\n", " frames=nt//subsamp,interval=40,blit=True,repeat=True)\n", "if _SAVE_GIFS:\n", " anim.save('schulz_sinx.gif',writer='imagemagick',dpi=80)\n", "# The following is necessary to just get the movie, and not an additional initial frame\n", "plt.close()\n", "out= HTML(anim.to_html5_video())\n", "plt.close()\n", "out" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is clear that phase mixing happens very quickly.\n", "\n", "We can also follow the evolution of the energy distribution of the particles. Note that even though the total energy of the system is conserved, the sum of the energies of individual particles does not have to be. The following movie shows the evolution of the energy distribution and clearly shows violent relaxation in action during the phase-mixing phase, with the energy distribution becoming much broader over this period. After phase mixing is complete, the energy distribution does not evolve substantially any longer." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false, "scrolled": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "figsize(6,4)\n", "fig, ax= subplots()\n", "ii= 0\n", "Es= wendy.energy(xt[:,ii],vt[:,ii],m,individual=True)\n", "a= ax.hist(Es,bins=31,histtype='step',range=[0.,0.004],normed=True,color='k')\n", "ax.annotate(r'$t=0$',(0.95,0.95),xycoords='axes fraction',\n", " horizontalalignment='right',verticalalignment='top',size=18.)\n", "subsamp= 4\n", "def animate(ii):\n", " Es= wendy.energy(xt[:,ii*subsamp],vt[:,ii*subsamp],m,individual=True)\n", " ax.clear()\n", " a= ax.hist(Es,bins=31,histtype='step',range=[0.,0.004],normed=True,color='k')\n", " ax.set_xlim(0.,0.004)\n", " ax.set_ylim(10.,10**4.)\n", " ax.set_xlabel(r'$E$')\n", " ax.set_yscale('log')\n", " ax.annotate(r'$t=%.0f$' % (ii*subsamp/20.),\n", " (0.95,0.95),xycoords='axes fraction',\n", " horizontalalignment='right',verticalalignment='top',size=18.)\n", " return a[2]\n", "anim = animation.FuncAnimation(fig,animate,#init_func=init_anim_frame,\n", " frames=nt//subsamp,interval=40,blit=True,repeat=True)\n", "# The following is necessary to just get the movie, and not an additional initial frame\n", "plt.close()\n", "out= HTML(anim.to_html5_video())\n", "plt.close()\n", "out" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The potential evolves as follows:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false, "scrolled": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ys= numpy.linspace(-numpy.pi,numpy.pi,1001)\n", "def init_anim_frame():\n", " line1= plot([],[])\n", " xlabel(r'$x$')\n", " ylabel(r'$\\Phi(x)$')\n", " xlim(-3.49,3.49)\n", " ylim(0.,3.9)\n", " return (line1[0],)\n", "figsize(6,4)\n", "fig, ax= subplots()\n", "line,= ax.plot(ys,wendy.potential(ys,xt[:,ii],vt[:,ii],m),'k-')\n", "txt= ax.annotate(r'$t=%.0f$' % (0.),\n", " (0.95,0.95),xycoords='axes fraction',\n", " horizontalalignment='right',verticalalignment='top',size=18.)\n", "subsamp= 4\n", "def animate(ii):\n", " line.set_ydata(wendy.potential(ys,xt[:,ii*subsamp],vt[:,ii*subsamp],m))\n", " txt.set_text(r'$t=%.0f$' % (ii*subsamp/20.))\n", " return (line,)\n", "anim = animation.FuncAnimation(fig,animate,init_func=init_anim_frame,\n", " frames=nt//subsamp,interval=40,blit=True,repeat=True)\n", "# The following is necessary to just get the movie, and not an additional initial frame\n", "plt.close()\n", "out= HTML(anim.to_html5_video())\n", "plt.close()\n", "out" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As a second example we consider the initial conditions\n", "\n", "$x \\sim [-\\pi/2,\\pi/2]\\,,$\n", "\n", "$v = -V_0\\,\\sin^3(x)+\\epsilon\\,V_0\\,\\mathcal{N}(0,1)\\,,$\n", "\n", "which are such that the system first collapses into two separate clumps that later merge (see [Schulz et al. 2013](http://adsabs.harvard.edu/abs/2013MNRAS.431...49S))." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true }, "outputs": [], "source": [ "N= 1001\n", "dx= numpy.pi\n", "V0= -0.001\n", "\n", "x= dx*(numpy.arange(N)-N//2)/N\n", "v= V0*numpy.sin(x)**3.+numpy.random.normal(size=N)*V0*10.**-6.\n", "m= numpy.ones(N)/float(N)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [], "source": [ "g= wendy.nbody(x,v,m,0.05,maxcoll=100000000)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [] } ], "source": [ "nt= 2700\n", "xt= numpy.empty((N,nt+1))\n", "vt= numpy.empty((N,nt+1))\n", "Et= numpy.empty((nt+1))\n", "xt[:,0]= x\n", "vt[:,0]= v\n", "Et[0]= wendy.energy(x,v,m)\n", "for ii in tqdm.trange(nt):\n", " tx,tv= next(g)\n", " xt[:,ii+1]= tx\n", " vt[:,ii+1]= tv\n", " Et[ii+1]= wendy.energy(tx,tv,m)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We again first check that the total energy is conserved:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEhCAYAAAC6Hk0fAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHJJJREFUeJzt3duTFOd5x/HfA0hCisUukMSJ7cTaRapKUlGZ5ZCkypV4\nw4KTSlUqsVYWrlTuIpC4hxiuzJVZCe6NjP6ALIpccU5VXlg0lVzEzmJBXImclGBBiZT4IJYF2RLI\ngicX7ztL08xx35np7uX7qZramZ4+PN3T07/ptw9r7i4AAFKsKroAAED1ESYAgGSECQAgGWHSB2Y2\nlnu9Jz6OF1UTqim/LgFlRZj0mJlNSHol9/qUu5+QdNHMDhRWHColvy4BZbam6AKKZmZTkmbc/Uyu\n+wFJFyWNSpp193OdjM/dZ83sYqbTqKQxScckzUva2ZPC7651T3x6UtJGSXvc/VCTfock7Y0vt0ma\nqs9b5r1FSbskveTus51OIy7LC5I2Sfq6u1/KjfeqpGFJ5+rjzQ0/IWnY3V/NdW/6WXQwzXvms92w\nyzWAdQkoL3e/Lx+SJiQdkPSmpB25905K2px5PdPluL/VpPtxSU/1YV4OSLot6Vacn8da9Hs883xE\n0kK9f4UNbva925LWdTINSTO5ZXY2W1+u36n6eHPdz0p6ttPPos00m85nu2GrsC7x4FG2x33bzOXu\ns+5+VFKjX6MT7n4+83rezHakTM/MRsJk/RtdDDPZYa9XJQ1JWu/uT7j75RY1LP3S9fBLfF7S07HT\nnvp8+p1f6aPtphHHuzW3zBYyy2xXrpT6r/RsbRPZ2jIafhZNpnkl817T+TSz0WbDNph+W4Nel4Ay\nuu+bufLiRm0+17ne7HMmNvdkr/Q0hZB4uc2o97r7vi7L2S7p1bZ9Sebu73XQ37DCXsHRXPeN8e/W\nekjEDa7rzrJoNY0tCr/8s+Zj9zOSNpjZlLsfjO/t9HAMKV/b1WyHNp/F2QbTvBSneVWt53OsxbBn\n1CN9XJeA0llRYRLbpuf93jb34wpNOJc7GM1wg25XFNrd1WAj2LCU3PQnFTZuMrMJb3C8oImOb09g\nZs8qbES3STrpDdrl3f2cmW3Ndd6i0ESj3PLZK+kv3f16B9NYlLShQVmb4t89kmbNbKekaUlfztU+\n6e6vmll+D6bVZ3Fad8Lhrmm6+7FW89lBvaVdl4CyWlHNXLGpYVu2eSh++Y93+OWXGm9kOhanPWJm\n+81sXfx1ekLSWTO7ovCruNdOufvL7v6qh4Pir5jZukY9ZptczGxvHPa1TLeRuCEdiXV3Mo2zunfD\nOaq4LGPoTCs0k03p7o32kHJ7JBmtPouzcXzNptlqPlvWG4cv3bqUMi6g31ZUmEhS3NBtM7PJzJf/\nfLvhMvLNH1LjX8DNpv9qPKZwzN2vx/b0DbHbRnc/1kUtHf0qbbBxW5T0TMsRmw1LmnT3P8iN61Lc\nkB6U9Hp9I9ZqGu5+TdKL9WMB8ZjFomITT/wcXnD3JxQCasbMNsfxPOO5s58ymn4W7abZaj47HbZs\n61IX0wUGrrBmrmanUeb66eiU0jx3P2RmF9T9l18KG5VGzRP5tu+eixusetOWSZqIy6AeKi7plewy\nixvC77p79lfwvDK//puYkvTF3PSH4oZW7n7JzBYlHTKzr7ebRlzmk2b2lO4ciL4YL7q7UA8jd38+\nfjbPmdmLCnsJzbT8LJpNs918djFsZdclYNAGHiax2WeLpEmF0zNb2Rt/JdeHnTKzuXa/0szspMKZ\nO7vr7fGd1ufh3P5888Sowmm9feXuz2dfm9kRb3K9SM5f5l4Pq/GZUfXxHlBo978eX48pNMmc0r17\nq8MKIdZ2GtnlHEPupMLB5vzG84TCRn5M0qa4TpjCsYT1ZqbYpNb2s2gyzabzmT2W1GrYTPdKrkvA\noA28mavNaZR5bU8pzYu/7r/q7uezzRRdlnk60wwjSSOt9qD6qG0zVzztdenXb2zWGamfEWRmY5a5\nJUdcFq9LumpmQ2a2RdJWhWWbD4wRhQPtl1tNI3ZbqDeJxWlMx434aUm7c+PdqXBB5Dfc/Whsxjmq\nEDqncmczNf0sWkyz1Xy2q1eZflbSugT0VdnP5urklNIl8ZfoXc0RsZniiJldzB2UHVPYyE0o/Bqe\nzhzP2CvpYDw9drvC2UhldsLu3KZlVHeH8G6FA9X74q/vV3R3U5pL2uXul83snJntl3RNYe9xT+ag\ndatpSCGIdprZRoXTW48pPLkWl3/9anMpnCV1V5NRHPeEwgHnBb9zPU6rz6LhNFvNZ7thc/Xcj+sS\nsCzmXsw/xzKzGYUmiFbHTDYrnPc/r3A20CtdnElTeV00cwFAoUp9Nlf89dfwlNL7xFzRBQBAJ0od\nJm1OKV3xvItbrwBAkUp7zKTVKaWSGt6WxMz4h/YA0CV3T77TQpn3TEbV+JTSlrwEd89czuMrX/lK\n4TVQf/F1UH81H1Wuv1dKFSa501ibnlI62KoAAO0UcdFiq9Mol05j9Q5PKQUAFG/gYeLhCuRzCvd+\nyr93MPf6vKT7IjzGx8eLLiEJ9ReL+otV9fp7obDrTPrBzHwlzQ8A9Fu8fdGKPgAPAKgIwgQAkIww\nAQAkI0wAAMkIEwBAMsIEAJCMMAEAJCNMAADJCBMAQDLCBACQjDABACQjTAAAyQgTAEAywgQAkIww\nAQAkI0wAAMkIEwBAMsIEAJCMMAEAJCNMAADJCBMAQDLCBACQjDABACQjTAAAydYUNWEzm5I04+5n\n2vQ3IulpSVclmbufGER9AIDODTxMzGxC0hZJk5Jm2vQ7IukFd38mvp4zszl3P9//SgEAnRp4mLj7\nrKRZM9vVQe8vSTqSeT3h7tf7UxkAYLlKe8zEzIYUwuO1ejeCBADKqbBjJh0YlbRoZjskrZc0Iulc\n3LMBAJRI2cNEkhbqB+nN7KyZPe3ul4srCwCQV9pmLkmLkoZzB9vnJT1XUD0AgCbKvGcyrxAo+W6j\nDfpdcvjw4aXn4+PjGh8f73VdAFBZtVpNtVqt5+M1d+/5SDuasNmMpKlW15mY2S13X515fVzSVXc/\n1KR/L2p+AKCKzEzubqnjKVUzl5mNmdlYptOL8bqUuq0KpwsDAEqkiIsWxyTtljQhab2ZTbv7sfj2\nbklDkvZJkrsfMrMj8eLFTZL2cPAdAMqnsGaufqCZCwC6syKbuQAA1USYAACSESYAgGSECQAgGWEC\nAEhGmAAAkhEmAIBkhAkAIBlhAgBIRpgAAJIRJgCAZIQJACAZYQIASEaYAACSESYAgGSECQAgGWEC\nAEhGmAAAkhEmAIBkhAkAIBlhAgBIRpgAAJIRJgCAZIQJACAZYQIASFZYmJjZlJnt6KL/CTOb7GdN\nAIDlGXiYxFA4IKnbYHhB0vo+lAQASLRm0BN091lJs2a2q9NhzGxC0sX+VQUASFGVYybDkq4WXQQA\noLHSh4mZTbr7q0XXAQBortRhYmZDYo8EAEqv1GEi6Rl3P1N0EQCA1gZ+AL5TZjYiaa7b4Q4fPrz0\nfHx8XOPj470rCgAqrlarqVar9Xy85u49H2lHEzabkTTVbM8jXlMyUn8pabfCGV2n3P3lJsN4UfMD\nAFVkZnJ3Sx1PqfZMzGxMktz9XP6gu5ltV4sgAQAUZ+BhEgNjt6QJSevNbNrdj8W3d0sakrQvN8yB\n2P+ImS24+zcGWTMAoLXCmrn6gWYuAOhOr5q5yn42FwCgAggTAEAywgQAkIwwAQAkI0wAAMkIEwBA\nMsIEAJCMMAEAJCNMAADJCBMAQDLCBACQjDABACQjTAAAyQgTAEAywgQAkIwwAQAkI0wAAMkIEwBA\nMsIEAJCMMAEAJFt2mJjZjszzCTN7qjclAQCqJmXPZLj+xN1nJW1ILwcAUEVdh4mZ7TGzs5JeMLM5\nMztrZnN9qA0AUBHm7t0PZDYkaVvcIykNM/PlzA8A3K/MTO5uyeNJ2fia2WZJOyW97u5nUotJRZgA\nQHd6FSYpB+AnJT0vySQ9Y2bPdjn8VPYgfpN+hszsQHxMm9nYcusFAPTPmpSB3f35+vMYLm2Z2YSk\nLZImJc206f2F+jTMbETSd81si7tfXl7FAIB+SDmb62rudUftS+4+6+5HJV1q1V8Mj4uZ4S5Jmpf0\ndJd1AgD6LGXPZJOZjSps4Lf0qJ6sYUlTko7mum/sw7QAAAmWvWfi7icUj5dIWnT3Yz2rKoz/nKSt\nuc5b1L5pDAAwYKnHTE5IOtGjWhqN/3z9uZntlXTK3V/r1/QAAMtTidupmNmwpEl3/4N+TQMAsHxV\nuZ3KlKQv9nH8AIAEXTdzmdkeSc9JGjKzQwrHTVzSSz2urT69A5Km3P16fD0Wj6c0dPjw4aXn4+Pj\nGh8f70dZAFBJtVpNtVqt5+NteAW8mX1L4Syti5Lm3f0bufeTb6diZjMKIXEm021MWjr4Xr92ZVHS\n2djLJklb3P3lJuPkCngA6EKvroBvtmeyS9Joo4sD40WHI5JOLmeCMTB2S5qQtN7MpjNngu2WNCRp\nX7zO5BXduX6lvge0aznTBQD0T7M9k5Pu/kzLAUPz0zZJx8tyhhV7JgDQnX7vmVzJTGhM0tX8Xoq7\nHzWznZJOS1qdWggAoLqanc21lFLx+MUXzeyKmX0tewqwu5+W1PRgeJncuCG9/37z999/X/rxj6UP\nP2zez61b0nvvSc12fm7elD74oPn7UnjvZz/rrGYAqIpmeyZ3bQ7jXsgud9/XoN+zDboVZmRE+qM/\nkv7kT6T//E/p7/9e+u53QwisXi19/OPSr/2a9LnPSWNj0t/9nfQv/yL9139Ja9dK169LP//z0u//\nvvSFL0g/+IH0yivSm29K774rPfigtGaN9NnPSn/6p9Kv/qr0V38lzc5Kb78tPfRQCJ3HH5d+7/fC\neK5dk/7hH6R//3fpUrwj2bp10q//urR9u/Tbvx3C7J/+SZqbC9P54ANp/Xrp058O4/nsZ8P4X3tN\nOntWunpV+rmfkx55RPrlX5bGx6XRUenKlTCOf/u3MI6HHw6PkZEwno99TDp3Tvre96SFhbBM1qyR\nHn1UevJJ6Td+I3T/3vekixeln/wkzNOGDaHe3/zNEIj/8z/h8d//Hepau1b6lV+RPvUp6ZOflBYX\nw/D/939h/j/2sTCOJ58M9fz4x9Ibb4T5vnlTGhqSPvEJ6Rd/Mcz3229L//u/ob9HHgnz+vGPS489\nFgL/hz8M77/7bnjv0UfD8I88cueHw3vvSWah/nXrwnsLC2GYtWvDvD/8cKh5zZrw3g9/GLo/9FD4\nrB96SLp9+86PiPoPhfrzRq+X2497qHfVqlDDqlWh2+3boV73sG5Job+6+vNW3crw3kcfhfrNwvfM\nLCzjVauk4eHwGdT7uXUrrBcffhjmf82asEzWrbvzmfzsZ3cekvTTn0oPPBA+54ceCuP86KO7+/vp\nT8Pf+nJety4MU6+z/mj2+oEH7vxo/Ogj3aXRPLd63q9hHn1UmprSwDU7ZjIn6Wl3fyvTbcrdDzbo\n94i7H+pvmZ0xM/+P/3D9zd9I3/xm2ID+2Z9Jv/VbYUN0+7b01lthI/a3fxs2mF/4QtjIjo3d+cK+\n9Zb0j/8ofetbYSP4538ufeYz0i/8Qlipr1wJ4fHNb4YN5pe+JP3xH4eN+apV4QvwxhtSrRYC4uGH\nw/tbtoQNl1nY2L3xhvSv/yp9+9uhn9/9Xel3fkf6pV8KX4Zr16QLF8J4vvOd0P1zn5O2bQsbzvff\nD1+Oy5elf/7nsGFfty4E1NhYWKk++CA8vv/9UMvNm2FeNm8O83P7dvhSvPtuCKALF8LG/Mknw4Z7\n48YwzI9+FOr9/vfDl/pTnwrhUQ+QGzekd94JIfD226GOxx8PwbxhQ3j/Rz8Kgfrhh6HbyEgIjwcf\nDBuXt98O4bGwEMb5iU+EGm/cCKH2gx+EeV27Nnyen/xkCKH6+++8E+Z17dqwPNetC5/nzZsh3G7c\nCNOtB8nt22H5vfNOWAbr14dgvnUr1HjzZnhIYTqdbGxS+/nww7AO1TeWq1aFx40bd4JfujuMsn/L\n/F49ENzDumkWlu+NG+HzWb367kc9zFevDsviJz8Jf2/eDMvkgQfuPKTwmX/0UVgH6uPM9vPAA+GH\nx+rVYdq3boX17tatzn4IfPBBGK6+fj344L3z2enzfg6zdq30F3+hjvX1n2OZ2W2FvZNFhWMipyRt\nbbRn0ixkisABeADoTr//OdbXJT0uqX5R4ouSnovHTabNbL+ZfSb2y9YbAO5zzfZM7rnKPF6ouEvh\n3/TulDSqGCTuXoqzudgzAYDuFP4/4DPhcsTdn0gtpBcIEwDoTuH/A97dr7n7X0ta9i1VAAArw7L3\nTJZGYDbk7td6VE8S9kwAoDuFN3OVEWECAN0pvJkLAIA6wgQAkIwwAQAkI0wAAMkIEwBAMsIEAJCM\nMAEAJCNMAADJCBMAQDLCBACQjDABACQjTAAAyQgTAEAywgQAkKywMDGzKTPb0UF/B8zsqfh/58cG\nURsAoDtrBj1BM5uQtEXSpKSZNv2elPRVdz8fX89I+nzfiwQAdGXgeybuPuvuRyVd6qD3iXqQRPOd\n7M0AAAartMdM4h7MfK7zoqRdBZQDAGihtGEiabhBtyuSRgddCACgtTKHyYaiCwAAdKbMYbLQoNvG\ngVcBAGirzGGyqMZNXfnjKACAgg381OBOufusmeWbukYlHW813OHDh5eej4+Pa3x8vOe1AUBV1Wo1\n1Wq1no/X3L3nI+1owuGakSl3P5PpNiZJ7n4uvp6WdCRzncmcu29vMU4van4AoIrMTO5uqeMp4qLF\nMUm7JU1IWm9m0+5+LL69W9KQpH3x9V5JB81sVNJ2SXsGXS8AoL3C9kz6gT0TAOhOr/ZMynwAHgBQ\nEYQJACAZYQIASEaYAACSESYAgGSECQAgGWECAEhGmAAAkhEmAIBkhAkAIBlhAgBIRpgAAJIRJgCA\nZIQJACAZYQIASEaYAACSESYAgGSECQAgGWECAEhGmAAAkhEmAIBkhAkAIBlhAgBIRpgAAJIRJgCA\nZIQJACAZYQIASLamqAmb2QFJFyWNSpp193NN+huStFfSVUnDks65++zACgUAtFVImJjZSUlfdffz\n8fWMpM836X2vux/NDDtlZnPufn0ApQIAOlBUM9dEPUiieTPb0aTfXbnX9b0ZAEBJDDxMzGxC0nyu\n86LuDY26DWY2lXm9MxdEAICCFdHMNdyg2xVJ25r0/6ykM2a2U9K0pC/3qzAAwPIU0cy1oZue417I\ntKQhSVOSNvWjKADA8hWxZ7LQoNvGZj2b2XFJU+6+Lz6fMbOtzZq6Dh8+vPR8fHxc4+PjadUCwApS\nq9VUq9V6Pl5z956PtOUEwzGT4+7+RKbblCR390O5fscUDtYfy3TbL2mTu+9rMG4f9PwAQJWZmdzd\nUscz8GaueI1IvqlrVNKpBr2P6t6D9Sf6URcAYPmKOjX4tJltzrwecfczUtgbiXskknRa0u7csDsl\nvTSAGgEAHRp4M5e0dFX7QUlzkrZLms5cwDglaajejBVD50uSLsTB5+vB02C8NHMBQBd61cxVSJj0\nC2ECAN2p7DETAMDKQ5gAAJIRJgCAZIQJACAZYQIASEaYAACSESYAgGSECQAgGWECAEhGmAAAkhEm\nAIBkhAkAIBlhAgBIRpgAAJIRJgCAZIQJACAZYQIASEaYAACSESYAgGSECQAgGWECAEhGmAAAkhEm\nAIBkhAkAIBlhAgBItqaoCZvZAUkXJY1KmnX3cy36HZH0tKSrkszdTwymSgBAJwoJEzM7Kemr7n4+\nvp6R9Pkm/Y5IesHdn4mv58xsrj4sAKB4RTVzTeTCYN7MdjTp9yVJX2sx7IpQq9WKLiEJ9ReL+otV\n9fp7YeBhYmYTkuZznRcl7WrQ75BCeLxW7+bu1/tbYTGqvjJSf7Gov1hVr78XimjmGm7Q7YqkbQ26\nj0pajHst6yWNSDrn7rN9rA8A0KUiwmRDF/2Oxr8L7n5GkszsrJk97e6Xe14ZAGBZzN0HO0GzSUkH\n3X17ptuUpBF3353rd0LSjLuvznQ7Kemiux9qMO7BzgwArADubqnjKGLPZFGNm7ryx1Hq3RYbdBtt\n0G9PFggAoHsDPwAfj3fkm7pGJZ1q0O8l3Rs8w2ocPACAghR1avBpM9uceT2SOSYyZmZjmfdejM1d\ndVsVThcGAJTEwI+ZSEun/B6UNCdpu6TpzAWMU5KG3H1fpv8jClfLb8r2i2LFkB9291cz3Tq+swHu\nD/E7PVP/wZjp3nRdKdN61Kj+uA3bG19ukzRVxvqbLfvM+737Drt75R+SDkh6StJ+SWNF19Om1qFY\n7wFJ0/l6KzYvZyU9m3l9UtLmzOuZomtsUvdIXM7PStpTpeWfWX+ejTVOlLV+SROxnjcl7ci913Rd\nKct61Kb+47n1aUHSY2Wpv1Xtuf569h0ubEXr4UIr/IPrst5Sr4RdzMdEDMPsinglP6+tVuSC6h6R\ndDLzeq6+zKuw/CUdyL2ekrSuzPVLmmmwMW66rpRtPcrXX/8xkuvnrKT98flCWepvtOwz7/X0O7wS\n7hrcza1ZChXvM3ax/trDCQbzCjexlKSdVZkXhRMhrtZfdHNng4K1uj1PFdal/PKsN0dI1ai/5bpS\nkfVoWCHE8zbG+i/mupet/rqefocrHSYVWfGyVsRKaGaTnmljjZrd2aDhadxFaHV7ngqtSxtiO3jd\nTnc/X6H6pdbrSunXIw/HELbmOm9ROCO19PVL/fkOVzpMVJEPrm6FrIRDyvyayejmzgZFWbo9j5lN\nmtn+zJmClVj+CsdK9sY7QRyQ9OXYvSr1S63XlSqsR8ruAZrZXt05yF36+vv1Ha56mJT+g8ur8koY\nPeONzwxZaNBtY7+L6VL29jyvuvsxSS+Y2WOqyPKP68+0woH4KYUzHKWK1B+1WleqsB4tMbNhSZPu\n/oexUxXq78t3uOphUoUPrqEqroTxmM9ck7e7ubNBURYVToO867iCpOdUgeUvSWZ2XOH/+zwh6YSk\nmXjNViXqj1qtK1VYj7KmJH0x87rU9ffzO1zYf1rskVJ/cG1UaiWMtkgaMbOdkkzh/Pr1ZiZ3f9nM\nGt3Z4Pigi2yh1e15Tqvkyz9ezHvB401O3f15M7ugEIZ/rZLXX+fus03Wla+5+2sVWI8kLV2PMZU5\n7jbWYt7KUn/fvsOVDpMKfHANVXQlVP6AnZltl3TK3V+OnU6b2ebML/+RJrvThXD3S3GPMGtY0nxc\n/vlf8qVa/gr15MPhhMK6VIX6sxqtK6+1eK8065G0dMPa1yVdjccgNilsqM+pxPX39TtcxLnPPT6P\nelp3n1s/V3RNbeqdVDi/eyg+tiie512leVG4IOqKwi7zU7HbkKQjChfNHcnOS1kesa6JzOs5SZ+u\nwvKPy3e6wfpUv06mVPVLGlPYA78Vl/P+3Lw0XFfKsh41q1/hOpPbsfutzPMdZam/1bLP9NPT73Ah\nt1PppVa3ZimbzHUm9YVu8fkudz9TpXmpsma356nC8o/HR74k6ULsNO937mtX+vqxclU+TAAAxav6\n2VwAgBIgTAAAyQgTAEAywgQAkIwwAQAkI0wAAMkIEwBAMsIEGAAzmzCzdUXXAfQLFy0CA2BmtxXu\nWHy96FqAfmDPBOizeLffqwQJVjLCBOi/bQq3uAdWLJq5gD6J/xJ4l6SnFW4s+bqklzz+PxJgJSFM\ngD4zswVJY+7+VtG1AP1CmAB9ZGajkt5099VF1wL0E8dMgP6aEMdLcB8gTID+2irpVNFFAP1GmAD9\ntVNxzyQekAdWJMIE6K/1kubj85EiCwH6iQPwQB+Z2f74dNHdXy60GKCPCBMAQDKauQAAyQgTAEAy\nwgQAkIwwAQAkI0wAAMkIEwBAMsIEAJCMMAEAJCNMAADJ/h9w2L6P7fT3SwAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot(numpy.arange(len(Et))/20.,Et)\n", "xlabel(r'$t$')\n", "ylabel(r'$E_{\\mathrm{tot}}$')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The evolution of the system in phase space is:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false, "scrolled": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def init_anim_frame():\n", " line1= plot([],[])\n", " xlabel(r'$x$')\n", " ylabel(r'$v$')\n", " xlim(-3.49,3.49)\n", " ylim(-2.49,2.49)\n", " return (line1[0],)\n", "figsize(6,4)\n", "fig, ax= subplots()\n", "line= ax.scatter(x,v,c=numpy.fabs(x),s=5.,edgecolors='None')\n", "txt= ax.annotate(r'$t=%.0f$' % (0.),\n", " (0.95,0.95),xycoords='axes fraction',\n", " horizontalalignment='right',verticalalignment='top',size=18.)\n", "subsamp= 4\n", "def animate(ii):\n", " line.set_offsets(numpy.array([xt[:,ii*subsamp],vt[:,ii*subsamp]]).T)\n", " txt.set_text(r'$t=%.0f$' % (ii*subsamp/20.))\n", " return (line,)\n", "anim = animation.FuncAnimation(fig,animate,init_func=init_anim_frame,\n", " frames=nt//subsamp,interval=40,blit=True,repeat=True)\n", "if _SAVE_GIFS:\n", " anim.save('schulz_sin3x.gif',writer='imagemagick',dpi=80)\n", "# The following is necessary to just get the movie, and not an additional initial frame\n", "plt.close()\n", "out= HTML(anim.to_html5_video())\n", "plt.close()\n", "out" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "The energy distribution evolves to a similar final state:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false, "scrolled": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "figsize(6,4)\n", "fig, ax= subplots()\n", "ii= 0\n", "Es= wendy.energy(xt[:,ii],vt[:,ii],m,individual=True)\n", "a= ax.hist(Es,bins=31,histtype='step',range=[0.,0.004],normed=True,color='k')\n", "ax.annotate(r'$t=0$',(0.95,0.95),xycoords='axes fraction',\n", " horizontalalignment='right',verticalalignment='top',size=18.)\n", "subsamp= 4\n", "def animate(ii):\n", " Es= wendy.energy(xt[:,ii*subsamp],vt[:,ii*subsamp],m,individual=True)\n", " ax.clear()\n", " a= ax.hist(Es,bins=31,histtype='step',range=[0.,0.004],normed=True,color='k')\n", " ax.set_xlim(0.,0.004)\n", " ax.set_ylim(10.,10**4.)\n", " ax.set_xlabel(r'$E$')\n", " ax.set_yscale('log')\n", " ax.annotate(r'$t=%.0f$' % (ii*subsamp/20.),\n", " (0.95,0.95),xycoords='axes fraction',\n", " horizontalalignment='right',verticalalignment='top',size=18.)\n", " return a[2]\n", "anim = animation.FuncAnimation(fig,animate,#init_func=init_anim_frame,\n", " frames=nt//subsamp,interval=40,blit=True,repeat=True)\n", "# The following is necessary to just get the movie, and not an additional initial frame\n", "plt.close()\n", "out= HTML(anim.to_html5_video())\n", "plt.close()\n", "out" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The potential in this case evolves as:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false, "scrolled": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ys= numpy.linspace(-numpy.pi,numpy.pi,1001)\n", "def init_anim_frame():\n", " line1= plot([],[])\n", " xlabel(r'$x$')\n", " ylabel(r'$\\Phi(x)$')\n", " xlim(-3.49,3.49)\n", " ylim(0.,3.9)\n", " return (line1[0],)\n", "figsize(6,4)\n", "fig, ax= subplots()\n", "line,= ax.plot(ys,wendy.potential(ys,xt[:,ii],vt[:,ii],m),'k-')\n", "txt= ax.annotate(r'$t=%.0f$' % (0.),\n", " (0.95,0.95),xycoords='axes fraction',\n", " horizontalalignment='right',verticalalignment='top',size=18.)\n", "subsamp= 4\n", "def animate(ii):\n", " line.set_ydata(wendy.potential(ys,xt[:,ii*subsamp],vt[:,ii*subsamp],m))\n", " txt.set_text(r'$t=%.0f$' % (ii*subsamp/20.))\n", " return (line,)\n", "anim = animation.FuncAnimation(fig,animate,init_func=init_anim_frame,\n", " frames=nt//subsamp,interval=40,blit=True,repeat=True)\n", "# The following is necessary to just get the movie, and not an additional initial frame\n", "plt.close()\n", "out= HTML(anim.to_html5_video())\n", "plt.close()\n", "out" ] } ], "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 }