{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Langevin Integrator Check\n", "\n", "The `toy_dynamics` subpackage provides an integrator called `LangevinBAOABIntegrator`, which is based on [a paper by Leimkuhler and Matthews](http://dx.doi.org/10.1093/amrx/abs010). This notebook uses the `toy_dynamics` package to check that the integrator gives the correct position and velocity distribution for a harmonic oscillator.\n", "\n", "Note that this particular test does not make use of the trajectory storage tools. It is mainly to show how to use the `toy_dynamics` subpackage, and has little connection to the main package. The trajectory generated here is extremely long, so in this case we choose not to store it. For an example using `toy_dynamics` with the storage tools, see ???(to be added later)???." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Imports" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import openpathsampling.engines.toy as toys\n", "import openpathsampling as paths\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Set up the simulation\n", "\n", "This the potential energy surface is $V(x,y) = \\frac{A[0]}{2}m[0] \\omega[0]^2 (x-x_0[0])^2 + \\frac{A[1]}{2}m[1] \\omega[1]^2 (y-x_0[1])^2$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "my_pes = toys.HarmonicOscillator(A=[1.0, 1.0], omega=[2.0, 1.0], x0=[0.0, 0.0])\n", "\n", "topology=toys.Topology(n_spatial=2, masses=[1.0,2.0], pes=my_pes)\n", "\n", "template = toys.Snapshot(coordinates=np.array([[0.0, 0.0]]), \n", " velocities=np.array([[0.1, 0.0]]), \n", " topology=topology)\n", "\n", "my_integ = toys.LangevinBAOABIntegrator(dt=0.02, temperature=0.5, gamma=1.0)\n", "\n", "sim = toys.Engine(options={'integ' : my_integ, 'nsteps_per_frame' : 50}, template=template)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "nframes = 250000" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set the initial conditions for the system, and initialize the sample storage." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "sim.current_snapshot = template" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "x1 = []\n", "x2 = []\n", "v1 = []\n", "v2 = []" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Run the simulation\n", "\n", "This might take a while..." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "for i in range(nframes):\n", " # generate the next frame (which is sim.nsteps_per_iteration timesteps)\n", " snap = sim.generate_next_frame()\n", " # sample the information desired to check distributions\n", " pos = snap.coordinates[0]\n", " vel = snap.velocities[0]\n", " x1.append(pos[0])\n", " x2.append(pos[1])\n", " v1.append(vel[0])\n", " v2.append(vel[1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Run analysis calculation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Build the 1D histograms we'll use:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": true }, "outputs": [], "source": [ "nbins = 50\n", "rrange = (-2.5, 2.5)\n", "rrangex1 = ((min(x1)), (max(x1)))\n", "rrangev1 = ((min(v1)), (max(v1)))\n", "rrangex2 = (min(x2), max(x2))\n", "rrangev2 = (min(v2), max(v2))\n", "dens = True\n", "(x1hist, binsx1) = np.histogram(x1, bins=nbins, range=rrange, density=dens)\n", "(x2hist, binsx2) = np.histogram(x2, bins=nbins, range=rrange, density=dens)\n", "(v1hist, binsv1) = np.histogram(v1, bins=nbins, range=rrange, density=dens)\n", "(v2hist, binsv2) = np.histogram(v2, bins=nbins, range=rrange, density=dens)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Build the 2D histograms:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true }, "outputs": [], "source": [ "(hist1, xb1, yb1) = np.histogram2d(x1, v1, [nbins/2, nbins/2], [rrangex1, rrangev1])\n", "(hist2, xb2, yb2) = np.histogram2d(x2, v2, [nbins/2, nbins/2], [rrangex2, rrangev2])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Run the analysis of the kinetic energy:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": true }, "outputs": [], "source": [ "instantaneous_ke = []\n", "cumulative_ke_1 = []\n", "cumulative_ke_2 = []\n", "tot_ke_1 = 0.0\n", "tot_ke_2 = 0.0\n", "for v in zip(v1, v2):\n", " local_ke_1 = 0.5*sim.mass[0]*v[0]*v[0]\n", " local_ke_2 = 0.5*sim.mass[1]*v[1]*v[1]\n", " instantaneous_ke.append(local_ke_1+local_ke_2)\n", " tot_ke_1 += local_ke_1\n", " tot_ke_2 += local_ke_2\n", " cumulative_ke_1.append(tot_ke_1 / (len(cumulative_ke_1)+1))\n", " cumulative_ke_2.append(tot_ke_2 / (len(cumulative_ke_2)+1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot our results" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Imports for the plots we'll use, as well as some parameter adjustment." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true }, "outputs": [], "source": [ "%matplotlib inline\n", "import matplotlib\n", "import matplotlib.pyplot as plt\n", "import matplotlib.pylab as pylab\n", "from matplotlib.legend_handler import HandlerLine2D\n", "import numpy as np\n", "\n", "pylab.rcParams['figure.figsize'] = 12, 4\n", "matplotlib.rcParams.update({'font.size' : 18})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we plot the distributions of the positions and velocities. These should match the exact Gaussians they're paired with." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtIAAAEjCAYAAAAfTOaLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl4VdXZ/vHvczISCFOYZUZkFAUhggjG2VptRRFbtfoq\nztYBrf211irqi7bOxalKFccqOCMOdQyigICRGWSeZQpDAiFkOOv3R0LfGBlOQnLWGe7PdZ3rhH32\nXuvekmwfdtZey5xziIiIiIhI1QR8BxARERERiUYqpEVEREREqkGFtIiIiIhINaiQFhERERGpBhXS\nIiIiIiLVoEJaRERERKQaVEiLiIiIiFRDSIW0mf3ZzMab2TIzC5rZ8qp0YmatytvINrP1ZrbTzOaZ\n2QNm1rh60UVERERE/LFQFmQxsyCQC+QAfYEdzrmOIXdidjXwGPAB8DWQD2QClwHrgUzn3KYqpxcR\nERER8STUQrq9c25l+ddzgbpVLKS7AbmVi2UzGw6MAR5yzv2xKsFFRERERHwKaWjH3iK6upxzC/dz\nx3lc+XvPQ2lfRERERCTcfD9s2Kb8fYPXFCIiIiIiVeS7kL4bcMCLnnOIiIiIiFSJt0LazG4FhgLP\nOOcm+cohIiIiIlIdXgppM7sCeAB4H7jBRwYRERERkUORGO4Ozexy4BngY2Coc670IPsffFoREZEI\n5Zwz3xnCSddsEYlmVb1mh/WOtJldBjwLfAIMcc4Vh3Kccy4mXnfddZf3DDoPnUukv2LlPJyL33rS\n9393fS/G7rnEynnoXCLzVR01XkibWRsz62JmCZW2/w9lc0Z/DpzjnCuq6b5FRERERMIlpKEdZnYx\n0A4woCmQZGZ/Kf94lXPulQq7vwwMBtoDq8uP/xXwL2AH8AYw1Ownd853Oufeq/5piIiIiIiEV6hj\npIdTVhxXdE/5+ySgYiHtgGClfXtTVoQ3pGx8dGWrgJgvpLOysnxHqBGxch6gc4lEsXIeEv1i6Xsx\nVs4lVs4DdC6xIqQlwn0yMxfpGUVE9sXMcHH4sKGu2SISjapzzfa9IIuIiIiISFRSIS0iIiIiUg0q\npEVEREREqkGFtIiIiIhINaiQFtmzB159FYKVJ5sRERER2T8V0hLfVq6E44+HSy+FGTN8pxEREZEo\nokJa4lJeXh5vX345u448kt3nnAO33Qbvvus7loiIiEQRFdISd9avXs0rbdsyeNw4Hh44kK7PPsum\ngQPhvZhfE0hERERqkBZkkbgS3LOH75s3p2HTpnSaMgWaNuW+++7jw4kTmbxqFfbFF9Cli++YEiO0\nIIuISPTQgiwiB/HN7beTXFJChwULoGlTAP70pz9R6hxLu3XTXWkREREJmQppiRvBYBD3z39i119P\nICnpv9sDgQD33Xcf9y1YgHvnHY8JRUREJJpoaIfEjUkvvshRw4fTIC8PS0v72efH9e3LpB9+IGnJ\nEmjRwkNCiTUa2iEiEj00tEPkADb/7W+sHjx4n0U0wGVXX820hg3h/ffDnExERESike5IS1zYtmkT\nhS1akD5tGvUyM/e5T35+Pjc1b84/Bw4k+dNPw5xQYpHuSIuIRA/dkRbZj9mjRrGtUaP9FtEA6enp\nFJ18Mnz9NeTnhzGdiIiIRCMV0hIX6r/+OrlDhx50vzMuuIB56enwn/+EIZWIiIhEMxXSEvMKFy6k\n3ebN9LjzzoPue9ZZZ/FiXh4lb74ZhmQiIiISzVRIS8xbf++9fNqiBY0PO+yg+zZs2JBVvXrhPvgA\niovDkE5ERESilQppiW3FxTSZMIEtQ4aEfEifs87ix3r1YNKkWgwmIiIi0U6FtMS2999nKdD7wgtD\nPuS0007jHefg3XdrL5eIiIhEPU1/JzGt+KSTuG7qVJ7KyyOpwmqGB1JSUsKAhg2Z1qABCWvXgsXV\n7GVSgzT9nYhI9ND0dyIVbdyImz6d9ccdF3IRDZCYmEjGwIEUlpTADz/UYkARERGJZiqkJXbNncva\nRo049oQTqnzooMGDWZGeDvPn10IwERERiQUhFdJm9mczG29my8wsaGbLq9OZmV1iZjlmVmBmG8xs\njJk1qU5bIge1cCGzi4oYNGhQlQ89/vjjmVFQAPPm1UIwERERiQWh3pEeBZwILAW2VacjMxsBvFB+\n/I3AP4HfAF+aWZ3qtClyICXz5/PNtm0ce+yxVT42MzOTSVu2UDJnTi0kExERkViQGOJ+HZ1zKwHM\nbC5QtyqdmFkGcC/wLXDK3idRzGwmMAG4CfhbVdoUOZiCmTPZ1aYNaWlpVT62Tp067O7UieKcnJB/\nSERERCS+hHRHem8RfQiGAHWAxys+zu2cmwgsBy4+xPZFfiZhyRLS+/Wr9vEZxx1H0rp1sGdPDaYS\nERGRWBGuhw37lr9P28dn04CuZlb124Yi+7N9OwkFBXQYPLjaTRx97LFsSkuDxYtrMJiIiIjEinAV\n0q3K39ft47N1gFXYR+TQLVrEsqQk+h7CHem+ffsyzznN3CEiIiL7FK5COg3AObev35EXVtxHpCYU\nzZ7N7KIievXqVe02evbsyYyCAopycmowmYiIiMSKcBXSBQBmlrKPz1Ir7iNSE7ZMnszW5s1JSdnX\nt1xokpOTyW/blvxvv63BZCIiIhIrwjUhwfry98Moe7iwosMAV2Gfnxk5cuR/v87KyiIrK6tm00nM\n2TN7Nok9ehxyO3X69SPh889rIJHEg+zsbLKzs33H8E7XbBGJBjVxzbYKk2iEdkD59HfOuY5VOGY4\nMAb4nXPu1UqfLQGKnHP7rHrMzFU1o8iG+vX56tZbGXbXXYfUzvPPPMPF119Pcn4+1NF051I1ZoZz\nznznCCdds0UkWlXnml3jQzvMrI2ZdTGzhAqb3wN2A783M6uw79lAJ+CVms4hcaywkEY7d3L4GWcc\nclPH9O/PisREWLiwBoKJiIhILAlpaIeZXQy0o2x2jaZAkpn9pfzjVc65ioXwy8BgoD2wGsA5t8XM\n/go8CHxuZq8BrYFbgAXAPw79VETKFM6dy0qgZ+/eh9xW9+7deaekhPY5OaT06XPI7YmIiEjsCHWM\n9HDKiuOK7il/n8RP7yg7IFi5AefcI2a2BRhBWeGcB7wO/Nk5pwcNpcas//xz1qan0zU5+ZDbSkpK\nYlPTpmydPJmWV1xRA+lEREQkVoRUSDvnTgy1wQPt65x7CXgp1LZEqiNv+nR2tW5dY+2VdOlCyaxZ\nNdaeiIiIxIZwTX8nEjZu4UKse/caa69Ov37UXbmyxtoTERGR2KBCWmJO/XXraDhgQI2112bwYNJ2\n7oT8/BprU0RERKKfCmmJLaWltNq5k7annVZjTfY86igWBwKwYEGNtSkiIiLRT4W0xJQdc+eSC7St\nwaEdbdq0Yb4ZO7XCoYiIiFSgQlpiytpPP2VdejqBQM19a5sZW1u1YvvXX9dYmyIiIhL9wrVEuEhY\n7Jg2jYI2bWq83WC3bgTnzq3xdkVERCR66Y60xBS3cCGBGhzWsVf9AQOov3p1jbcrIiIi0UuFtMSU\nmp6xY692gwaRXFgI27bVeNsiIiISnVRIS8xwwSCt8/Npc+qpNd52zyOPZL4Zbt68Gm9bREREopMK\naYkZm8unp2vao0eNt92kSRMWJyWxY8qUGm9bREREopMKaYkZ6z//nDX16oFZrbS/tVUr8qdNq5W2\nRUREJPpo1g6JGXnTp2MtW9Za+6Vdu2JalEVERETK6Y60xAy3YAGlRxxRa+2n9etHgzVraq19ERER\niS4qpCVm1FuzhrrHHFNr7bfJzITiYti0qdb6EBERkeihQlpiRovt22l2wgm11n7Xbt34IRCARYtq\nrQ8RERGJHiqkJSYUbNpEo9JSWg8cWGt9tG3blkWlpRRqCjwRERFBhbTEiDWff87alBQSkpNrrY+E\nhAR2NGnC9pkza60PERERiR6atUNiwubp00nKyKj1fko6dGCPZu4QERERdEdaYsTuuXMpbtu21vtJ\n7dGDpFWrar0fERERiXwqpCUm2IoVJHfrVuv9NOnfn4ZbtoBztd6XiIiIRDYV0hIT6m3aRKM+fWq9\nn059+1LkHGzeXOt9iYiISGRTIS1RLxgM0nzXLloNGlTrfR1xxBEsCQYp/eGHWu9LREREIpsKaYl6\na1esoJVz1A3D0I60tDTWpqaydfr0Wu9LREREIpsKaYl6a6dMYWtKCtTi1HcV5TVvTl5OTlj6EhER\nkcgVciFtZUaY2UIz221mq83sITNLC/H4umZ2u5nNMbM8M9tsZt+Y2aXVjy8CW2fOZHvjxmHrr7RD\nBw3tEBERkSrdkX4MeBiYB/weGA/cCEw42IFmZsDHwN3At8AtwL3l/Y81s/urFlvk/+xZsICi1q3D\n1l9qz56krF0btv5EREQkMoVUSJtZd8qK5zedc+c7555zzv2BsoL4JDP7zUGaOBYYCPzDOXelc+5f\nzrnRwCBgBXB19U9B4p2tWEHCEUeErb/GmZk03ro1bP2JiIhIZAr1jvSF5e+PVdo+BigALj7I8fXL\n33+suNE5VwJsAXaFmEPkZ+pu2kT6UUeFrb92ffrgSktBxbSIiEhcC7WQ7gsEgRkVNzrn9gCzgH4H\nOX46sB34o5kNNbM2ZnZE+ZCOPsBdVYstUqa0tJTmO3fS7LjjwtZnx06dWOIcJYsWha1PERERiTyh\nFtKtgC3OueJ9fLYOaGJmifs72Dm3HTgb2EbZ2OpVwCLgWuA859zzVUotUm7tmjV0BOp07x62PlNS\nUliflsaWqVPD1qeIiIhEnlAL6TRgz34+K6ywz4HsouxBxQeBIcBwYCnwmpmdHGIOkZ9Y+d13WEIC\nNGoU1n7zmzUj7/vvw9qniIiIRJb93kWupABoup/PUivss09mdiQwBbjJOTemwvbXKSuux5hZJ+ec\n29fxI0eO/O/XWVlZZGVlhRhbYt2W6dPJbdiQ9DD3qynwZF+ys7PJzs72HcM7XbNFJBrUxDXb9lO7\n/nQns4+Bk4G0ysM7zOxroLNzrvkBjn8euBRo4pzbVumz0cD1wOHOuRX7OHZ/9bUIL595Jv3WrqXr\nnDlh7feNm27imHHj6LhhQ1j7lehiZjjnzHeOcNI1W0SiVXWu2aEO7ZhRvm9mpQ5TgKOp9BDiPrQq\nf0/Yx2eJld5FQrd8Oda5c9i71RR4IiIiEmohPa78/eZK268C6gCv7t1gZh3NrEul/RYABvxPxY1m\n1hA4h7KHEJeFmEXkv9I2bKDekUeGvd92mZkklZTA9u1h71tEREQiQ0h3gZ1z88zsSeB6M3sL+BDo\nDtwAZDvnXquw+xdAG3569/kx4BLgb2bWC/gGyACuAJoD1znngod6MhJfSktLaZqXR0Zm5sF3rmHt\nO3RgPtBt4UKSBwwIe/8iIiLiX1WGU9xE2SqEVwFnUraQyj/4+RzQrvz1fxucW21m/YA7KRtrfQGw\nm7I5qEc4596rVnqJa6tXr6ZTIEBqGKe+2ysxMZEf69YlY8oUWquQFhERiUshF9LlT488Wv460H4d\n9rN9BXBZldKJHMDS+fM5wTlo3dpL//nNm5OvKfBERETiVqhjpEUizuYZM9hWrx4k+nlOVVPgiYiI\nxDcV0hK18mfPpqBFC2/9p/boQcratd76FxEREb9USEvUckuX4jp29NZ/xrHHago8ERGROKZCWqJW\nyrp1pPXo4a3/tv37k1ZcDDt3essgIiIi/qiQlqhUUlJCk7w8Gvft6y1Dm3btWAbsnjfPWwYRERHx\nR4W0RKXVq1fTOSGB5G7dvGUIBAJsqFePTd984y2DiIiI+KNCWqLSksWLaVdaCh7HSAPkN2umKfBE\nRETilAppiUrrcnIoTk6G9HSvOTQFnoiISPxSIS1RacesWexs1sx3DFJ69CBVU+CJiIjEJRXSEpVK\nFy8m2L697xg0zszUFHgiIiJxSoW0RKXkNWtI8fig4V5tBgwgvagICgp8RxEREZEwUyEtUaekpITG\n27fTsE8f31Fo3a4dK83YPX++7ygiIiISZiqkJeqsXLmSLklJJHXt6jtK2RR4deuyUVPgiYiIxB0V\n0hJ1lixZQkfnvE99t9fOpk3Jmz3bdwwREREJMxXSEnVWLlhAemkptGrlOwoApe3aUawp8EREROKO\nCmmJOttycsjPyIBAZHz7pnTtSqKmwBMREYk7kVGJiFRB8cKFlLRt6zvGfzXs04cGubm+Y4iIiEiY\nqZCWqJO8ejVJ3bv7jvFfhx1/PM0LCsA531FEREQkjFRIS1QpLi6m6bZtpPft6zvKfx3WpQv5QMHy\n5b6jiIiISBipkJaosnLlSnokJ5MUAYux7BUIBPixTh3Wf/217ygiIiISRiqkJaosWbKEw52Dzp19\nR/mJHY0bs+2773zHEBERkTBSIS1RZcW8eTQoLYU2bXxH+YmiNm0oXLjQdwwREREJIxXSElXycnLI\na9IkYqa+2yuxc2cCq1b5jiEiIiJhFHI1YmVGmNlCM9ttZqvN7CEzS6tCG43Kj1lS3sYmM/vCzAZW\nL77Em5KFCylp3953jJ+pf9RR1Nu40XcMERERCaPEKuz7GHAD8BbwENANuBE4GjjlYAebWVtgEpAG\nPAcsBhoAvYDDqpRa4lbqmjUkn3uu7xg/03zAAGznTt8xREREJIxCKqTNrDvwe+BN59ywCttXAqPN\n7DfOudcP0syrlN0BP9I5t6maeSWOFRUV0WzHDur36+c7ys+07NeP4mCQXbm51M3I8B1HREREwiDU\noR0Xlr8/Vmn7GKAAuPhAB5vZYGAg8Hfn3CYzSzSzOlVKKnFv+fLl9EhOJrFrV99RfiaQlMTG5GTW\naAo8ERGRuBFqId0XCAIzKm50zu0BZgEHu0X4C8ABa83sfWA3sMvMfjCzi6oWWeLVkiVL6BQMRtzU\nd3ttbdCALTNmHHxHERERiQmhFtKtgC3OueJ9fLYOaGJmBxom0gUwyu5gNwR+B1wO7AFeNrNLQ48s\n8Wr13LmkOQctW/qOsk+FLVuya+5c3zFEREQkTEJ92DCNsqJ3Xwor7JO3n33Sy9/zgBOdcyUAZvYu\nsBy4D3gxxCwSp/Jzcshr1oymZr6j7FunTrBkie8UIiIiEiahFtIFQNP9fJZaYZ/92U3Z0I7X9hbR\nAM657WY2AfidmXVxzv2wr4NHjhz536+zsrLIysoKMbbEkuAPP1DSoYPvGPtV98gj2TN5su8Y4lF2\ndjbZ2dm+Y3ina7aIRIOauGabc+7gO5l9DJwMpFUe3mFmXwOdnXPND3D8U8DVwA3OuacqfXY/8Edg\noHNu2j6OdaFklNj3SMOG/M+FF9L4qacOvrMHGyZOZMuQIfQs3tcIKIlHZoZzLkJ/hVI7dM0WkWhV\nnWt2qGOkZ5Tvm1mpwxTK5pE+2BNW0ykbI916H5/tXetZU+LJfhUWFtIyP58Gffv6jrJfzfr3p01J\nCbs0n7SIiEhcCLWQHlf+fnOl7VcBdSibIxoAM+toZl0q7fcukA9cXHElRDNrCfwaWOycW16V4BJf\nli9fTvekJBIicOq7vQIZGVggwIqcHN9RREREJAxCKqSdc/OAJ4FzzewtMxtuZg8DDwPZzrnXKuz+\nBbCg0vHbgT9QtoLht+VLjf8JmAokUbbYi8h+LVm8mE6lpXDEEb6j7J8Zm9PT2Th1qu8kIiIiEgZV\nWSL8JmAFZXehzwS2AP8A7qq0nyt//XSjc2PMbDNl46HvoWxe6inAb/Y1NlqkorWzZxNISIAIXzWw\noFkz8mfP9h1DREREwiDkQrr86ZFHy18H2m+/0yo4596lbJiHSJXs/P578lu0IC1Sp74rF2zfnpLF\ni33HEBGJCdtWrGDDrbeyeds2FiYmsq5RI3a1bEnnbt3o27cvvXv3JiEhwXdMiWNVuSMt4o1bvJhg\nx46+YxxUavfuJI8bd/AdRURkn0pLS5kwYQLf//WvXLNgAatatya5TRtO27qVjFmzqLNjBz+mpzPF\nOZ4HWv3mN1z65z/Tpk2bg7YtUtNCfdhQxKs669aR2quX7xgH1fiYY2i4davvGCIiUemTTz7hxCOO\noNHw4fxh2zYyPv2UM1av5qRvvqHD7NnU37yZpK1bafvhh/zmjjv4W+/ejBg7lqL27Znasyd5H3/s\n+xQkzqiQloi3e/duDtu5M6Knvtsro18/WhcXs1NT4ImIhGz7tm3c/ItfMOO3v+XTLVvIuuEG6i9f\nTsrJJ/9853r14Nhj4ZZbqP/559TdtYsGn33G2gYNyD/rLDYMGACLFoX/JCQuqZCWiLds2TK6JyUR\n6FJ5VsXIE+jQgdbAsh/2uUiniIjsNXMm3H8/2wcNItikCX+dMoU/nXYaKd98A3ffDSkpobUTCNDk\nxBM5/5tvWPbBB4z94Qd29ulD8Mor4ccfa/ccJO6pkJaIt2TxYtqXlEDnzr6jHFxKCnmpqaz79lvf\nSUREItekSXDmmSydOpVb5sxh2jPPkLFjBwmvvQY9e1a72cGnn87VS5dySWYmb3/2GcEePeBf/6rB\n4CI/pUJaIt76nBxKk5OhYUPfUUKSl5HBdi3KIiKyb9u2wSWX8OGwYZzw3Xdc/8UXnHnFFTXWfOPG\njRn/2Wd8fsYZnH/YYZT+7//Cgw/WWPsiFamQloi3a9YsdrZs6TtGyIrbtGGPhnaIiPycc3DVVczp\n1InrP/iAyZMnc8wxx9R4N4mJiTz11FN0+sUvOD01lZIxY+COO8r6F6lBKqQl8i1Zgjv8cN8pQpbc\npQsJq1b5jiEiEnmee46t06Zx3pIlfP7553SsxWlNzYy///3vZJ57LmfUqUPpxIlw000QDNZanxJ/\nVEhLxKu7fj11jjrKd4yQ1T/6aNI3b/YdQ0QksixaRNGttzKksJAParmI3svMGDVqFIcPGMA56ekE\nZ86E4cOhpKTW+5b4oEJaIlp+fj6tCwpoUAu/+qstjY45hsOKisjLy/MdRUQkMuzZQ8E553B7aSl/\nmzCBI444ImxdmxlPPvkkSU2bcl2nTrh16+DPfw5b/xLbVEhLRFu8eDHdk5OjYuq7vQKHH87hgQA/\naJy0iAgABSNGMHn1avqPHcuAAQPC3n9CQgIvvfQSU2bP5rkTToCxY2H58rDnkNijQloi2qIFC2hb\nXAxRNEaa5s2p6xxLZ83ynURExLvgs8+SN3YsM666iqHnn+8tR7169Xj33Xf5y+jRLP/1r3VXWmqE\nCmmJaD/OmMGetLSylayihRk7Gjdm8/TpvpOIiPj16qvsvPVW/l/v3vz54Yd9p6Fjx468+uqrnDxx\nIiWTJ8PUqb4jSZRTIS0RbdesWRS2bu07RpWVtGnD7nnzfMcQEfHnnXcovOEGzq1XjwfffZeEhATf\niQA45ZRTuHrECB5IT8fdequmxJNDokJaIlry0qUEevTwHaPKEnv1InXZMt8xRET8+OgjSq+6il8C\n9779Ns2aNfOd6Cduu+02Pm/Vio0rVsBbb/mOI1FMhbRErNLSUppt2kT9gQN9R6myhoMGcVhuLiWa\nYklE4k12Nu7SS/lDp04MuvFGLw8XHkxCQgIvvvwy1xcWUnjzzVBU5DuSRCkV0hKxVq5cSZ+EBJKj\naOq7vZL79uWohARWrFjhO4qISPgsXQrDhjHhoov4urSUv/zlL74T7Vfr1q25+PnnmbJ1K4WPPOI7\njkQpFdISsRYtWECX0lI48kjfUaqua1falpayeM4c30lERMLDObj+ejZfdhnDX36ZV155haSkJN+p\nDmjIkCF8+YtfUHz33bBtm+84EoVUSEvE+nHqVPbUqQMNG/qOUnXJyWxt1IgtX33lO4mISHi88QZu\n/XrOnTSJkSNH0iVK5v//w/PP825CAmuHD/cdRaKQCmmJWIUzZ5LXvr3vGNW2s1MnSr7/3ncMEZHa\nl5cHt9zCW6ecgiUnc9111/lOFLIGDRrQ7NlnSZwwgd2ffOI7jkQZFdISsVJ++AGLxmEd5RKPPpq0\npUt9xxARqX133snO44/nmpdfZsyYMQQC0VVenH7hhfx74EB2XnABFBT4jiNRJLq+0yWuNN2wISpn\n7Nir0Qkn0GLzZpzmKBWRWJaTg3vtNa7IzWXEiBFRM6SjskvffpuvCgv58fLLfUeRKKJCWiLSli1b\n6FpSQoPjj/cdpdoaDBpE99JSNm/e7DuKiEjtKC2Fa69l+jnnsGDjRm677TbfiaotIyOD4GOPkfD2\n25R8+aXvOBIlQiqkrcwIM1toZrvNbLWZPWRmaVXt0MzqmNkKMwua2eiqR5Z4sHjOHNo7h3Xr5jtK\ntVnr1tQJBFimJWhFJFaNGUNJIMC5Eybw7LPPkpyc7DvRIRl61VX8o0sXdg4bBrt2+Y4jUSDUO9KP\nAQ8D84DfA+OBG4EJ1ejzXqAxoN93y35tmjSJLQ0aQDRflM34sUkTcrOzfScREal5ublw55083r07\np55+Ov379/ed6JCZGRePH88n+fnsuvFG33EkChy0kDaz7pQVz2865853zj3nnPsDcAtwkpn9JtTO\nzKwPcBNwF2DVzCxxoHDGDHa0a+c7xiHb1bGjZu4Qkdg0Zgx5gwcz6r33uO+++3ynqTHdunVj0TXX\nsOe112DSJN9xJMKFckf6wvL3xyptHwMUABeH0pGZBcqP+RB4J9SAEp+Sf/ghOhdiqSShd2/N3CEi\nsae0FJ5+mrtzc7n11ltp1aqV70Q16tZRo7itbl0KLroIdu/2HUciWCiFdF8gCMyouNE5tweYBfQL\nsa9bgCMou7stckBNNmyg/nHH+Y5xyBqdcAIt9bChiMSaiRPZXq8e76xaxYgRI3ynqXF169bl7DFj\n+GrHDkruvtt3HIlgoRTSrYAtzrnifXy2DmhiZokHasDMOgAjgbudc2uqnFLiyu7duzm8oIAWp57q\nO8oha3XqqXQqKmJXXp7vKCIiNSb4+OP8LS+Phx56iNTUVN9xasWvf/1rXs3MpOiJJ2DePN9xJEKF\nUkinAXv281lhhX0O5GlgGfBoiLkkji2ZNo30QICkTp18RzlkiY0akZuczHKtliUisWLRIgpnzCCn\nY0eGDBniO02tMTPueOop7gSKLrsMgkHfkSQCHfBOcrkCoOl+PkutsM8+mdnFwCnAIOdcadXilRk5\ncuR/v87rW9yVAAAgAElEQVTKyiIrK6s6zUiU2PDpp6Q0bkwXi43nUTc2a0b+l1/C0KG+o0gty87O\nJluztOiaHeMKH3mEf5aU8ODo0ViMXKf3p0uXLrgrr2T1v//N4c88A9de6zuS1KCauGbbwVZdM7OP\ngZOBtMrDO8zsa6Czc675fo5NBtYA3wIVB1G1Br4EXgbuoWzoyI79tOG0Mlx8ee+UU2iTl0ef6dN9\nR6kRX590Ert37+ZUzScdd8wM51xsVxqV6Jod4/Lz2dW0KfcOHcrfXnnFd5qw2L59O78+/HA+Ky0l\naf58iLEHK+X/VOeaHcrQjhnl+2VW6iwFOJpKDyFWUoeyu9m/BJZUeH1J2TzSvwMWA8OrElpiW8ri\nxST07u07Ro1Jzcyk7vLlvmOIiByyHx94gC+c49ZH42ekZsOGDblw1CheS0/HaW5pqSSUQnpc+fvN\nlbZfRVmh/OreDWbW0cy6VNhnFzAUOL/8fe/rWsrmkf6o/M/VWdhFYlSzjRtpcuKJvmPUmBannsph\nubm+Y4iIHBIXDFL06KMUXnEFTZvub8RnbLriiit4okEDdn7zDUxQySL/56BDOwDKl/K+HniXsnmg\nuwM3AJOdcydX2G8l0MY5l3CQ9toBK4AnnHMH/Oedfk0YX/K2b4dGjai3ZQuBjAzfcWpEsLiY3cnJ\nFK9eTcM2bXzHkTDS0A6JJdPuv5/GI0fSPi+P5JQU33HCLjs7m+cuuIAXGzYksHAhBEJdHFqiRW0N\n7YCy1Qj/QFkB/QQwDPgHcHal/RyhL/1dlX0lTiz9/HN2JyXFTBENEEhKYnVaGqs//NB3FBGRaikq\nKmLHqFEUXXFFXBbRUPbg7O7jj2fjzp3w0Ue+40iECKmQdmUedc51c87Vcc61cc7d5pwrqLRfB+fc\nQWcCcc6tcs4lOOduqm5wiU2bv/iCjc2a+Y5R43IPO4ztX33lO4aISLW8eM89DCgqouff/+47ilcP\nPvQQI3fsoPD++31HkQih30tIRCnJyaGwc2ffMWpcSffuMHeu7xgiIlW2adMm0h96qGwu5Xr1fMfx\nqkOHDjS7/noKvv8eZs3yHUcigAppiShpy5eT2i/UVeejR/pxx9Fg9WrfMUREquzZa6/lzIQEmjz4\noO8oEeGPd9zB0wkJbLn9dt9RJAKokJaI4ZyjVW4uzWNgafDKWp95Ju3y8nBaGUtEosis778n6/33\nSbj3Xqhf33eciJCenk77++4j+dNPcWvX+o4jnoU0a4dPegI8fvy4ciUNO3QgdfduLDX14AdEEecc\nmxMSCH77LS1i8I677Jtm7ZBo5pzj9iOP5LatW2m8ejUkhrIYcnwIBoOMb9GCXscdR/d33/UdR2pI\nbc7aIVLrlr7/Phvq1o25IhrKfjjXNG7M6okTfUcREQnJ2+PGcdXSpTR45hkV0ZUEAgE6jR5N8/ff\nZ/eWLb7jiEcqpCVibJ00ie0xPM/yro4dyybzFxGJcLt372b29dfToEcPEs46y3eciNTvN79hacuW\nfHnZZb6jiEcqpCViJM2eDTG0NHhlSZmZpC5a5DuGiMhBPXnffdyyaxeNn38eLK5GJ1VJ20cfpeuH\nH7JOD5PHLRXSEjFarFlDo9NO8x2j1jQ//XRabtzoO4aIyAGtW7eOpIceIuHss+Goo3zHiWgthw4l\nuXlz3rj0Ut9RxBM9bCgRIS83l0CTJtTZto2Ehg19x6kVpXv2UJiaqqXC44geNpRodMu55/K/H31E\n2tKlcNhhvuNEvN1jxzL/6qsJTp5M5rHH+o4jh0APG0rUWvree2xOTY3ZIhogISWFVXXrsuK993xH\nERHZp2lTp/LrDz8k4Y9/VBEdojqXXEL7Vq34+OKL0T8i448KaYkI2z/9lA1t2/qOUeu2tm3L9uxs\n3zFERH4mGAzy5UUX0aNVK1L++lffcaJHQgKNX3+d61au5O3nn/edRsJMhbREhEBODqUx/KDhf/Xu\nTWDOHN8pRER+5t3HHuOaNWtoPGGCprurokD//hT/4hfk33wzBQUFvuNIGKmQlojQfM0aGsXgioaV\nNT7pJJpoJSwRiTA78/Jodvvt5F15JYGePX3HiUotx47lnOJiXr31Vt9RJIz0sKF4V7h9O8FGjSA3\nl7TGjX3HqVWFubkEmzTBduygjpbbjXl62FCixYRf/ILuM2Zw+IYNuht9CLb8/e+svOMOWi5fzmF6\nqDzq6GFDiUrL3nqLVXXqxHwRDZCakcGmlBQWT5jgO4qICACrsrM57j//od6bb6qIPkRNbruNps2b\n8/GwYb6jSJiokBbvNn/0EZvbtfMdI2y2tGnDjx995DuGiAgEg+wYOpTZZ5xBi6ws32miXyBAk3Hj\n+NX06eR8+qnvNBIGKqTFu4ScHALxNPfm0UcT/O473ylERJgzYgQlu3Zx/Jtv+o4SM+oOHMjmk05i\njabDiwsqpMW7VmvX0uLss33HCJtmp51Go1WrfMcQkTi3JzeX5k8+ScHf/kZKWprvODGl6/jxDNy6\nlY8feMB3FKllethQvNq6ciXJHTpQp7CQhJQU33HCIvjjj+xo1YrSTZto0rSp7zhSi/SwoUSyr08+\nmeIffuBEzSRUK5aOGMG6f/6Tfrm5pOkfKlFBDxtK1Fk6fjwr69ePmyIaINCyJaXJycz94APfUUQk\nTq2fOZPuX35Jx3//23eUmHX4Aw/QOSmJ8cOH+44itUiFtHi147PP2N65s+8YYZfbujUbPv7YdwwR\niVOLLriA+f37027wYN9RYldSEqmPPUbmG2+waP5832mklqiQFq9S580jZeBA3zHCr3dvPXAoIl7M\nGDuWI1eu5Bg9YFjrGl92GQ3ateP9oUP14GGMUiEt3jjnaLtxI23OPdd3lLDb+8BhMBj0HUVE4khR\nURG7b7yRdZdeSlqrVr7jxD4zmr/8MpcuW8b4sWN9p5FaEFIhbWVGmNlCM9ttZqvN7CEzO+joeTPr\nbGb3mNlUM9tkZnlm9r2Z3R7K8RK7lk6dSiPnaDFokO8oYdfopJM4yjkWLVrkO4qIxJE3rr6aTsEg\nRz39tO8ocSPxuOMIZGWx+qab2LFjh+84UsNCvSP9GPAwMA/4PTAeuBEIZXm2y4GbgKXA3cAfgEXA\n/wLfmFn8PGUmP7H09ddZ07w5BOLwFyMdO9LYjG+1MIuIhMmSH37gyJdfJvmhh7A4esA7EjR55hmu\nLS7mgREjfEeRGnbQCsbMulNWPL/pnDvfOfecc+4PwC3ASWb2m4M08QbQ2jn3O+fck865Z51zvwVG\nAb0APc4ap3ZNmkTJ0Uf7juFHIMCO9u1Z9+GHvpOISBxwzjHx17+mSevWNL3mGt9x4k+HDiRceSXd\nXn+dGTNm+E4jNSiUW4EXlr8/Vmn7GKAAuPhABzvncpxz+fv4aBxgQM8QMkiMcc7RYPFimp55pu8o\n3tTp35/SmTP1AIqI1Lr//PWvXLpsGc0mTACLq6nNI0ade+7hvORkHvztbykqKvIdR2pIKIV0XyAI\n/OSfUM65PcAsoF81+25T/r6xmsdLFFu1ciW9iopoGUcrGlZWf/BgehQXs3z5ct9RRCSGbc3Joff9\n95P7yCMk9urlO078atSI1Acf5J7cXO4fNcp3GqkhoRTSrYAtzrnifXy2DmhiZolV6dTMAsCdQDGg\n2eDj0Mz33iM1KQlr1853FG+sTx8yk5P56quvfEcRkViVn8+uk0/mm0GD6HzDDb7TxD0bPpyO7duT\n+8gjzJ0713ccqQGhFNJpwJ79fFZYYZ+q+AeQCfzVObekisdKDNg4cSLbDj88vn/F2KMHLQsKmPLF\nF76TiEgsKi3lx5NOYnIwyOkTJ/pOIwCBAMljxvC3QICbLrmEkpIS34nkEIVSSBcA+3u8N7XCPiEx\ns3uB64FnnHMPhHqcxJY6M2eSeuKJvmP4lZJCaYcObFQhLSK1oODGG1k2Zw7tJ0ygbr16vuPIXn37\nUue3v+XGLVt49NFHfaeRQxTKkIz1QDczS9rH8I7DKBv2EdI/qcxsJPAX4Dnn3HWhhhw5cuR/v87K\nyiIrKyvUQyUCrVi+nJPy82l2+eW+o3iX0r8/nd99l2XLltGpUyffceQQZWdnk52d7TuGd7pm++de\neokdL77IZ1ddxcgTTvAdRyqx++7j7Dff5JRRozj77LPp2rWr70hxqSau2XawGQPK7yDfDgx2zn1T\nYXsKkAtkO+fOOmhHZncBdwEvOOdCrqDMzGlWg9gy7s47OemRR2ianx/fQzsA/vUvvh01iu9uu43r\nrgv535YSJcwM51xcfZPrmh0B1q6lsHt3ftu0Ka/Nn09qaurBj5Hw+9e/2HD//ZzVoAFTpk0jOTnZ\nd6K4V51rdihDO8aVv99caftVQB3g1QoBOppZl30Eu5OyIvrFqhTREpuK33yTrQMHqogG+NWv6L1p\nE59/8IHvJCISC5yj8NJLGV1ayh3jx6uIjmSXX07zJk24qLSUu+++23caqaaDFtLOuXnAk8C5ZvaW\nmQ03s4cpW+kw2zn3WoXdvwAWVDzezK4HRgKrgC/M7KJKr1Nq6mQk8hUVFdFl8WKaX3ml7yiRoVkz\n7OijSfzyS/bs2d8zvSIioQm+9BLrpk9nzy23cMwxx/iOIwcSCGBPPsmNP/7IxDFjmDx5su9EUg0H\nHdoBYGZG2R3pq4D2wBbgdeAu51xBhf1WAG2cc4kVto0FLjlA85OccycdoG/9mjCGTHnrLXpecAH1\nd++GpCTfcSLD6NFMvPde0saN46ST9vujIFFIQzskrH78kV2dO3PD4Yfz7MyZJCZWaWZa8eWee9j6\n+uv0Lyhg+qxZNGzY0HeiuFWda3ZIhbRPuijHlrfOOIPOa9fSa94831Eix5o1FHTpwj3XXsvfHn7Y\ndxqpQSqkJWycIzcri5dnzuTchQtp27at70QSqmAQzjmHSStXMrpzZ958801MQx+9qK0x0iI1pumU\nKaQMG+Y7RmRp04Zg+/Zseftt30lEJErtfP55cqdO5fCXXlIRHW0CAXj5ZQbt3k3PnBxGjx7tO5FU\ngQppCZuV8+fTe+dOOml2ip9Ju/hijtu4kWXLlvmOIiJRJrhxI8XXX8/EIUM467zzfMeR6mjQgMA7\n73BnXh7vjRzJtGnTfCeSEKmQlrCZ9cADrGnVisQmTXxHiTiBoUM514w3x4/3HUVEoolz/JCVxUcZ\nGfz+5Zd9p5FD0bMnCU8/zfvJyVw9dChbtmzxnUhCoEJawibhww8JnHOO7xiR6YgjSGzWjEUvvug7\niYhEkXlXXEHRsmWcNGWK5iGOBcOGUfeSS3gjOZlh551HUVGR70RyECqkJSxWLlvGgNxcDh8xwneU\niJV20UX0Xb2a5cuX+44iIlFg1Wuv0XzsWNy4cbRo1853HKkp999P544duXHtWm644Qb08G5kUyEt\nYTHtkUcoaNyYRC2DvV+BoUMZlpjIW2++6TuKiES4LYsWkXTJJcz6/e85esgQ33GkJiUmYm+8wa/M\naD5xIk8++aTvRHIAKqQlLEreeYeiM87wHSOyHXUUdevWZdZLL/lOIiIRrGDnTpb278/Svn05VTM8\nxKZGjQh88AF37tnDZ3feyYcffug7keyHCmmpdYsWLmTApk20u/FG31Eimxmpv/0tx6xaxYIFCw6+\nv4jEndLSUt459lgaJSczaNIk33GkNnXpQuJrrzHejLt+9zvN5BGhVEhLrXvvgQfISEsjqV8/31Ei\nXmDoUC6sU4cxY8b4jiIiEcY5xxPnnccZS5fS4dtvMT1cGPtOPZXku+/mi3r1uPhXv2LhwoW+E0kl\nKqSlVu3Zs4fg+PHY2WeDVmo6uP79aQJMefFFCgsLfacRkQjhnOORyy7j4g8+IOWNN0ju0MF3JAmX\n668n/cwzyW7Vil+efjorV670nUgqUCEtteqjF17g2uJiGvzlL76jRIdAgMShQ7mmcWPeeecd32lE\nJEI8duONXPzvf5Py9NPU+9WvfMeRcDKD0aNp3bo1X9aty68GD1YxHUFUSEutqnP33Ww47TTo3t13\nlOhxzTX8JjeXl55+2ncSEfHMOccjt93Gec88Q9qoUdS74grfkcSHpCR4913anXoq2UVFXDpwICtW\nrPCdSlAhLbVo8Qsv0GvjRjq88ILvKNGlVy+STz2VQbNnk5OT4zuNiHgSDAa547rrOPvxx2l0662k\n33ab70jiU2IijB5N47vu4sP8fP40YACLFy/2nSruWaRP9G1mLtIzyj6UlrKiSRMW/PKX/PKVV3yn\niT6LFlHQrx9XnXgir0yY4DuNVJOZ4ZyLq4cDdM2uGcXFxVx7ySXc/MEHHH7JJaQ+/rieM5H/89ln\n7D73XG4344L//If+/fv7ThQTqnPN1h1pqRUbRo5k/a5dZP3zn76jRKeuXUk6+2yO/OILTYUnEmd2\n7NjBJWecwY2ffELXIUNURMvPnXIKdWbM4N569Zh94olMfPdd34nilgppqXmbNpH24IPMu/pq6tar\n5ztN1Eq6915ucI7RI0f6jiIiYbJkyRKu69WLp6ZPp8fVV5P4/PMqomXfunSh3rx5XNC7Nw0uuIAn\nRo7UcuIeqJCWGrf1qqv4dyDABffe6ztKdOvUicShQ+n6wQfMmjXLdxoRqWX/+fhjXunThzHbt9No\n/HgS7rsPEhJ8x5JI1qgRDSdPptfw4Qy5/37+ePrp5Ofn+04VV1RIS41yU6ZQ+tFHJI0aRcOGDX3H\niXrJ997LNWb89ZprdKdBJEYVFRVx54034oYM4dY2bUibOxd+8QvfsSRaJCTQ4KmnaDpmDHdOnszI\nbt2YM2eO71RxQ4W01Jxt28g/7zwebtWK/9Fy4DWjbVuSL72UYStW8O9//9t3GhGpYUuWLOHyo4/m\n2uee44SLL6b+rFnQtq3vWBKFki+5hPSpUxm5Zw+LMjP55+23U1pa6jtWzNOsHVIzSkooOvVUXpox\ng27/+Q8DBw70nSh2rF9PSdeu9K1Th0/nzaNp06a+E0mINGuH7E9xcTEPP/QQW0aN4h4z6jz3HDZs\nmO9YEgvy89l+++3YP//JxBYtyHznHTr37es7VVTQrB3ijfvTn1g4dy5Lrr5aRXRNa9WKxCuvZGxG\nBv9zySUEg0HfiUTkEEydOpVBRx/NsY8/zqi2bUnLyVERLTUnPZ2Gjz9OvWXLOLpVKxpkZjLh1FPZ\ntW2b72QxSYW0HLpXXmHHCy9wW7t23Hv//b7TxKa77+aohg35n++/5+GHHvKdRkSqYcmSJZx//vnc\nOWQI/9m2jaxf/pKU776Dzp19R5MYlNC2LT2+/Rb75BPazZ/PriZNmH3WWZQsXeo7WkxRIS2HZuZM\n9lx/PecAz7z5JsnJyb4TxaZ69Qh89BG/at6cenffzTtvv+07kYiEaPny5VxzzTVclJnJHStW8Elx\nMQ3+/ndszBioU8d3PIlxTU85haPWr2fNSy+xeNYsdnbpwvo+fSidMAE0hvqQhVxIW5kRZrbQzHab\n2Woze8jM0sJxvESgDRsoPPNMrgYe+vhjOnTo4DtRbGvQgJQvv+TS1q1Z/rvf8dlnn/lOJCIHkJOT\nw4UXXsjv+vTh2q+/ZlogwFFnnYUtWQK/+53veBJnjrnoIoauWcPMd97h+R07mDdsGDszMii6+WbQ\nLB/VFvLDhmb2D+AG4C3gY6AbcCPwlXPulNo6Xg+uRKCiInj5ZXb+5S88kZ/PgA8/5IQTTvCdKn5s\n3syuzExGb9rE4S+8wPnnn+87keyHHjaMP/n5+bz++uu8+tRT9FyzhpuaN6fT5s0EbrkFrrsO6tf3\nHVEE5xxTp05l/MiRtJs8mUsTE0lt0YI6V1yBXXghtGnjO6IX1blmh1RIm1l3YC7wlnNuWIXtvwdG\nAxc6516vjePj/aIcUXbtgjFjcA8/zLKUFP6cl8efP/qIPscc4ztZ/NmwgcL+/Xlh61byrriCWx54\ngMTERN+ppBIV0vFh165dfPDBB3zxwgvU++ILLq5fn547d5KQlYUNGQIXXgh16/qOKbJPq1ev5sWx\nY5n/zDMMLSzkzMJCgt27U/fKK7Hzz4fGjX1HDJvaLKT/F/gzMMg5N6XC9hQgF8h2zp1VG8fH40U5\n4uTmwpNPwhNPsK1XL36/bh2bWrfmpZdeomXLlr7Txa9169j9xz9SMn4836Wlcdhf/0rnm26CpCTf\nyaScCunYFAwGmTNnDp988gkL336b9t9/z29TUmgbDGK//CUpw4bB6adDvXq+o4qELBgMMmXKFCa8\n8QY7xo3j7Px8Ti4pYftRR1Fv+HDSL7gAYnyhtdospD8GTgbSnHPFlT77GujsnGteG8fH0kU5Ozub\nrKws3zFCt2oVPPII7uWX2TJwIPcUFvLG3Llcdtll3HfffZhFf30QdX8n++Dy85l6663Mfv55LkhM\npPD882l5551YlM4EEAt/J3upkI5u2dnZnHDCCWzcuJGcnBxmT5rE+smTKZgzh8zkZH4NNAgESDjv\nPJIvuAAGD4YI/c1QrPxcxcp5QGSfi3OO+fPn85833sC9/TY9Fi5kkHOsbtWK7YMHU++iizjipJNI\nTU0FIvtcqqI61+xQf+JbAVsqF8Hl1gEDzCzROVdSS8fHhKj4Risqgu+/Jzh6NMGJE5neqxd/adKE\nlfPmMWLECJa/9x4PPPBATBTRECV/Jwdh6ekc9+yzfNC0KR+asfvxxzn3tdfY0a4dweHD6TBiBAlR\nNDNALPydSHQJBoNsXr+eHxcsYMuCBeyYM4eChQt5Z948cvfsoa1zHGfGyWYUHnYYKSeeSGrv3nDW\nWdC3LwQifwKsWPm5ipXzgMg+FzOjZ8+e9OzZE+6+m5KSEuZMmcL6sWNp9NVXtB03jhXBIOvq1WNH\ny5a8XlrK9nPPJSMzk9a9e9OqdWtSUlJ8n0ZYhFpIpwF79vNZYYV98mrpeKkpwSDk5eG2b2fPxo1s\nX7aM3TNnYrNnk7Z4MQ03bGBdUhJjg0E+6dCBY/v04X+HDeO4446LmeI5ViUlJXHxyJG4e+9l2qRJ\nLH7gAbredx/177iDr5s1w9q1I71TJ9Latye9Uycade5MRrt2pGZkYHXrgv5+JdqUlsLOnRRv3Urh\n5s3syc2lKDeX4m3bKMnLo2THDkrz8ynesYPibdtw27ZBXh4J+fkkFRSQUlhI3T17qF9aSmMgOTGR\njNRUdmVk4Nq25bsjjuC0K68kvWdPOPxwaN6cFP2cSBxKTEykz+DB9Bk8uGxDcTF1vv+e1C+/ZGdO\nDk0nTeKosWNp/OijpJWUsAXYnJDAjtRUCurVo6hRI4IZGQQzMqBZMwJNm5KakUFao0bUbdKEuhkZ\npDRoQGr9+qSmp1MnLY3k5OSoqDtCLaQLgP2tS5xaYZ/aOj4ssvv1I231av77S8lKv578yZ8O9Nk+\nPgdYU1DAlCee2P+x+zjmYO06wJwjAAScI8G5sncgJRj8v1dpKcnBIHWdYxewHcgzY1dyMuvr12dD\ny5bkZ2aSPmAA3fr25aYjj+SejIyf5ZHIZ2YMyMpiQPmdju3Tp3PUiy+yffFiir7/nqQvviCwcyeu\nsJD80lIckAIUmlFoRkkgQEkgQKkZpRX/XOG9NBDA1eQFrkJba3fu5Nunn/75LjXQTUHr1mR9910N\ntCSRYNKgQfSeOpV8YFcgwO5AgILERIoSEylKSqI4OZmS5GRcWhpWvz6BRo1I7NCB5KZNSWnWjLRW\nrUju1Imkbt1IatqURmY0qtB+w5EjSb/mGl+nJxK5kpJIycykQ2YmAM1GjqTDyJFlnxUX03zDBuos\nWULe4sXsWr6cwjVrKP3xRxJXriRl9mxSdu0ioaiIxOJikkpKSC4tJdk5Ep0jGSih7C5rMZCyYwcp\nETzbTVSMkT5oQBGRCBWPY6R9ZxARqa7aGiM9AzgVyAS+2buxfNaNo4Hs2jo+3v4nJCISzXTNFpF4\nEuoTEuPK32+utP0qoA7w6t4NZtbRzLpU93gRERERkWhQlZUNRwPXA+8CHwLdKVupcLJz7uQK+60E\n2jjnEqpzvIiIiIhINKhKIW2U3VG+CmgPbAFeB+5yzhVU2G8FZYV0YnWOFxERERGJBiFPfunKPOqc\n6+acq+Oca+Ocu61yEeyc61C5iK7K8ftjZk3NbKyZzTazXDPbbWZLzOxfZtYp1POIBGbWysz+bGbZ\nZrbezHaa2Twze8DMom4tTjO72sxeMbOFZlZiZqW+Mx2IlRlRnne3ma02s4fMLM13tqoo/x4ab2bL\nzCxoZst9Z6oOM+tsZveY2VQz22RmeWb2vZndHoV/J0eU/ywsMLPtZrar/PvsYTNr4TtfuMXKdVvX\nbL90zY4sumZXaiNaVqAysyOA54CpwCpgN9AZGE7Z7F3HOucW+UsYOjO7GngM+AD4Gsin7EHMy4D1\nQKZzbpO/hFVT/luIxsD3QEfgsMpDeyKJmf2DsmFFbwEfA92AG4GvnHOn+MxWFWYWBHKBHKAvsMM5\n19Fvqqozs/uB64AJwDTKZjw6EbgAmA30d87tbx76iGJmJwG3U3YeaymbxelI4HJgB3C0c26Lv4Th\nFSvXbV2z/dI1O7Loml2Jcy6qX5R9MwaBJ3xnqULmbkCzfWwfXn4uD/jOWMXzaVvh6/eBUt+ZDpC1\nO1AKjK+0/ffl/+1/4ztjFc6lfYWv5wLLfWeq5nn0AdL3sf3e8r+r63xnrIFzHFr+/fUH31ki4RVt\n121ds71m1TU7wl66Zv/0Ffnrmh7c6vL3RgfcK4I45xa6fd+92Du7Sc9w5jlUzrnVB98rYlxY/v5Y\npe1jKFsU6OLwxqk+59xK3xlqgnMuxzmXv4+PxlG2DktU/TzsR9Rdp2pZVP330DXbK12zI4yu2T8V\n6jzSEcPMEoH/3979h95V13Ecf773AzbLlrQtMmRW9Gtapv2giFK3atEPcxUhtX5ZW4T/FJISkTOE\nkOgHWlYEo0aaBeEaUcrIclitmMR0xIIQHTjdcMkIGyua7/743C+73H1v7nv87nzOuff5gC/bPedy\neZ7MTIAAAAZaSURBVH22e9/f9/2czzlnGbCYcojwesoN/n5VMdZ8OWfw58GqKSbbzEzY7uGNmfnv\niNgDvL5KKs1m5vNwqGqKBqJcI//ZlDu3ngfcSKlTv66Zq5YJrtvW7NPPmt0fU1mze9dIA+soh6Jm\nHASuzsyfVMozn75C+Y/bWjvIBDsbOJwjd9gcOAC8KSIWZeZ/W86lIRGxALiOsvauj5/tTwPfHnr8\nELAhM/8w5vmTblLrtjX79LNm98A01+zWG+mIWAZ8nlJ8TsVNmXlk6PEu4G2UG7mspixuPysiFmZm\nq2cez8NYhl/rasqanO9n5s55injK5nMsHXcGMO4kiGNDz/lnO3E0xk2Uk7m+mJl/rx2mgW3APsoM\nx4XAZcCKqomegUmp29Zsa7ZOm+mt2RUWcK+iHKY5foo/L36a13sB8Djwvb6OhfJN6DiwHVjY9jjm\neSxdP3HlAeCxMft+Nhjboto5G4yrtyeuzDKWGwbvxe/WzjKPY3oV5Zf+tbWzNMw/EXXbmj3r61iz\n64zLmt3hn7nU7N5c/u7/iYjbgQ8Az8rZD/90VkRcSTlp4i7g8r7lHxURvwTelR29lFJE3AWsBc4Y\n/beOiN8DL83M51cJ9wxExF7K+793l1IaFhHXUw4PbsnMjZXjzKuI2AWcnZmramfpgr7WbWt2u6zZ\n3WbNnsMNWTpuKbAQeE7tIHMREZ8EfgDsANb3vSD3xG7K+/4NwxsHJxq8hpETWtSeiNhMKcg/mrSC\nPLCUcu1eFb2r29bsKqzZHWXNLnrTSEfEyjHbV1O+rT6Ymf9oN1VzEfEJyqzG3ZRZjf/UTTQ1Zi5X\n9bmR7ZsoH5rb2o0jgIi4DtgMbM3MK2vnaSoiZp0Zi4hLKZeE2tVuoromqW5bs6uxZneQNXvouX1Z\n2hER3wLeTrlc0sOcuFbhRyknTb43M39TLeAcRMRlwB2Uu+ZcS7nb17AnM3N768Eaioj3ABcMHm4A\nXkb5lgpwJDNvqRJsjIi4GbgK+AXl0jarKXfNujcz19bMNhcRsYGyTjIoNydYDHxzsHt/Zt5aK9tc\nRMRVlLOl91PeN0+NPOVQjz7bd1DW//6WMp4lwGuBK4AngUsyc2+9hO2alLptza7Lmt0t1uyR1+hR\nI70G+CxlgCsphwQPAPcA38jMffXSzc3Q4ZBx9vdp3VRE/BD42JjdnRtLRARldmMTcC5wGPgpsDkz\nj1aMNicR8TvgrWN278zMNW3maepp3j/Qr7F8EPg48GrKGd9JKc47gK9n5iMV47VuUuq2Nbsua3a3\nWLNHXqMvjbQkSZLUJb1ZIy1JkiR1iY20JEmS1ICNtCRJktSAjbQkSZLUgI20JEmS1ICNtCRJktSA\njbQkSZLUgI20JEmS1ICNtCRJktSAjbQkSZLUgI20JEmS1ICNtCRJktSAjbQkSeqUiHh/RHw5Iu6M\niLMG214XEYcjYnXtfNIMG2lJktQZEfFC4JzMvAE4F7h4sOtRIIFXVoomnWRR7QBS2yJiE7AceDnw\nY2AVsBI4H7gmMw9UjCdJ0+6dwNaIOB94CfAngMx8NCK2APtrhpOGOSOtqRIRG4H7M/OrwHeAnwOH\ngfuAKyjNtCSpkszckplHgI3Ajsw8OLT7X8Bf6iSTTuaMtKbN8zLzz4O/rwKOZ+b2iFgKXJKZ91bM\nJkk64UPAl0a2LcjMp2qEkWYTmVk7g1RFRNxMWYe3vnYWSdIJEbECOARclJl7BtvWAccyc2dErKUs\nyXsj8EfgPODBzNxaK7Omk0s7NM0uBe6pHUKSdJIjwFFgIUBEPBd486CJXg4syczbgWWUXmYfZZme\n1CpnpDU1ImIBsAa4G1gBHAQuzMz7B/uvycyvVYwoSRqIiPXAh4E9wGLgxsw8FhFLMvPY4Dl7gYsz\n84mKUTXFXCOtafIZ4BbgFcA7KLMdjwBExPuAv9aLJkkalpnbgG2zbJ9popdTJgSfiIgAlmbm0ZZj\naso5I62pEREXAF8A/gY8AJxJmaF+GHgoM2+tl06SdCoGEx9nUpZ0rMvMj0TEu4H7MvNQ3XSaNjbS\nkiSpNyLiU8CLgL3AW4DdwGOZuaNqME0lG2lJkiSpAa/aIUmSJDVgIy1JkiQ1YCMtSZIkNWAjLUmS\nJDVgIy1JkiQ1YCMtSZIkNWAjLUmSJDVgIy1JkiQ1YCMtSZIkNfA/RVYveKOTsxgAAAAASUVORK5C\nYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Boltzmann info as a in exp(-ax^2)\n", "boltzmann_vel1 = 0.5*sim.integ.beta*sim.mass[0]\n", "boltzmann_pos1 = 0.5*sim.integ.beta*sim.mass[0]*sim.pes.omega[0]**2\n", "plotbinsx1 = [0.5*(binsx1[i]+binsx1[i+1]) for i in range(len(binsx1)-1)]\n", "plotbinsx2 = [0.5*(binsx2[i]+binsx2[i+1]) for i in range(len(binsx2)-1)]\n", "plotbinsv1 = [0.5*(binsv1[i]+binsv1[i+1]) for i in range(len(binsv1)-1)]\n", "plotbinsv2 = [0.5*(binsv2[i]+binsv2[i+1]) for i in range(len(binsv2)-1)]\n", "lx1 = np.linspace(min(plotbinsx1), max(plotbinsx1), 5*len(plotbinsx1))\n", "lx2 = np.linspace(min(plotbinsx2), max(plotbinsx2), 5*len(plotbinsx2))\n", "lv1 = np.linspace(min(plotbinsv1), max(plotbinsv1), 5*len(plotbinsv1))\n", "lv2 = np.linspace(min(plotbinsv2), max(plotbinsv2), 5*len(plotbinsv2))\n", "f, (ax1, av1) = plt.subplots(1,2, sharey=True)\n", "px1 = ax1.plot(lx1, np.sqrt(boltzmann_pos1/np.pi)*np.exp(-boltzmann_pos1*lx1**2), 'k-', plotbinsx1, x1hist, 'r-')\n", "px1 = ax1.set_xlabel('$x$')\n", "pv1 = av1.plot(lv1, np.sqrt(boltzmann_vel1/np.pi)*np.exp(-boltzmann_vel1*lv1**2), 'k-', plotbinsv1, v1hist, 'r-')\n", "pv1 = av1.set_xlabel('$v_x$')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the above, you should see that the exact answer (black line) matches up reasonably well with the calculated results (red line). You might notice that the left graph, for position, doesn't match quite as well as the right graph, for velocities. This is as expected: the integrator should impose the correct velocity distribution, but sampling space correctly requires more time to converge.\n", "\n", "The plots above check the $x$ degree of freedom; the plots below do the same for $y$." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtIAAAEjCAYAAAAfTOaLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl4VdW9//H3NwmZmIcIhFkRZJAZAgQwDCpiHYoDamm1\ntbV1qNhe/XW8lVo76xW9na22tmqVKooDAiKEIUAIIAphEASZZAqKEBIyrt8fJ3BjGpKTkGSdc/J5\nPU+eE/bZe6/PuRd3v6yzBnPOISIiIiIiNRPlO4CIiIiISDhSIS0iIiIiUgsqpEVEREREakGFtIiI\niIhILaiQFhERERGpBRXSIiIiIiK1oEJaRERERKQWgiqkLeA7ZrbFzPLNbI+ZPWJmiUFef56Z/ans\nugIz221ms8ys5bnFFxERERHxw4LZkMXMHge+DbwMzAf6APcCy5xzk6q5NgnIAjoAfwKygf7At4BN\nQKpz7tQ5fAYRERERkQYXU90JZtYXuAd4yTl3Y7njHwFPmNlNzrkXqrjFj4AuwM3Oudnlrl8FPA98\nF/hF7eKLiIiIiPhRbY+0mT0M/AAY65xbWe54HHAUSHfOfaGK6zcAPZ1zzSocN+AksN85d2HtP4KI\niIiISMMLZoz0MKCUwPCMM5xzBcAGYHg118cB/zF0wwUq+HzgfDNrE1RaEREREZEQEUwhnQzkOOeK\nKnlvP9DOzKoaIpINtDazAeUPmtkgoHXZH7sGE1ZEREREJFQEU0gnAgVnee9UuXPOZhbggH+b2RVm\n1sXMrgBeAAqDuF5EREREJORUO9kQyAOSzvJefLlzKuWcW2Fm04AngDcAA4qBvwLtgWuB42e73syq\nX1ZERCREOefMd4aGpGe2iISzmj6zg+mR/pjA8I0mlbzXicCwj+JqQr0MdAYGA2OBZOfcXWXHioEd\n1Vwf9j8PPvig9wz6LPos4fATKZ/DucZbT/r+v7v+Lkbm59BnCd2fSPkstRFMIZ1Vdt6I8gfLVu0Y\nRIVJiGfjAt53zmU453LMrAOBwjrdaR1pEREREQkzwRTSL5a93lfh+B1AAvDc6QNmdr6Z9a7uhmVL\n3z1R1v7Pg4sqIiIiIhI6qh0j7ZzbZGa/B+42s5eBeUBfAjsdpjvn/lXu9MUENl+JPn3AzJoCa4BX\ngF1AS+BmYAjwQ+fcsjr6LCEtLS3Nd4Q6o88SmiLls0TK55DwFyl/FyPlc4A+S6iKpM9SU8FuEW4E\neqTvALoDOQRW3XjQOZdX7rxdQBfnXEy5Y02AZ4AUoCOBiYlZwKPOuUVBtO1qO25FRMQnM8M1wsmG\nemaLSDiqzTM7qELaJz2URSRcqZAWEQkftXlmBzNGWkREREREKlAhLSIiIiJSCyqkRURERERqQYW0\niIiIiEgtqJAWEREREakFFdIiIiIiIrWgQlpEREREpBZUSIuIiIiI1IIKaRERERGRWlAhLSIiIiJS\nCyqkRURERERqIahC2gK+Y2ZbzCzfzPaY2SNmlhjk9U3N7Idm9r6ZHTezI2aWYWa3nlt8ERERERE/\nzDlX/UlmjwPfBl4G5gN9gHuBZc65SdVca8AyYCTwdyATSARuBlKAXzvnflDF9S6YjCIiocbMcM6Z\n7xwNSc9sEQlXtXlmV1tIm1lfYCPwsnPuxnLH7wGeAG5xzr1QxfUjgZXA/zjn7i93PAbYBrR2zrWp\n4no9lEUkLKmQFhEJH7V5ZgcztOOWstdZFY4/CeQB06u5vkXZ64HyB51zxUAOcDKIDCIiIiIiISUm\niHOGAaVAVvmDzrkCM9sADK/m+jXAMeD/mdluAkM7EoCvAkOAb9Y0tIiIiIiIb8EU0slAjnOuqJL3\n9gOjzCymrIf5PzjnjpnZVcBTwOxybx0HrnPOvVbT0CIiIiIivgUztCMRKDjLe6fKnVOVk8Am4LfA\nF4HbgR3Av8xsYhAZREREpLH79FP49a+hXz947z3faUSC6pHOA5LO8l58uXMqZWYXE5hsOMM592S5\n4y8QKK6fNLMLqpqdMnPmzDO/p6WlkZaWFkRskcodOHCARx99lOXLlxMfH89VV13F3XffTUJCgu9o\nEubS09NJT0/3HcM7PbOlLuXn5/P8zJm0+ec/mXD4MNt79+aisWNp9sMfwptv+o4nYawuntnBrNox\nH5gIJFYc3mFmK4ALnXPtq7j+aeBWoJ1z7tMK7z0B3A30dM7tOsv1mgEudSb9+efZdvvtbP/Wt7j2\nuus4efIkf/nLX9i0aRNvvfUW559/vu+IEkG0aofIudm5bRtbhw0jtaiInKlTybnxRmYvX84LzzzD\njpgYEl56CcaM8R1TIkRtntnB9EhnAZcCI4CMco3FAYOA9GquTy57ja6i/WByiJyTJUuWsPnrX+fO\n4mKievQ48/C9/PLL+cMf/sAll1xCZmYmycnJ1dxJRETq2/79+/n3iBHc0rkzLdeto2ViIhcAKdde\ny9SpU/l/V1zBz+66i1bvvQfWqP69KiEkmDHSL5a93lfh+B0EVt947vQBMzvfzHpXOG8zYMBt5Q+a\nWSvgWuBT4MPgI4vU3N69e7lj2jTuiIkh6q234KGHYOfOM+/fdddd3HnnnUydOpWiosrm1YqISEMp\nKiri/iuu4NvFxXRZsAASPz8VKzU1lS++/DKHN2/myLPPekopEvzOhqeHYLwKzAP6EtjpcLlzbmK5\n8z4Cujjnossd6wqsB1oBzxPo1W4LfB3oBtzlnPtzFW3ra0I5J845pkyZwndLS7m0bVt4/nl45BGY\nNw/eeedMT4ZzjsmTJzNu3Dh+9KMfeU4tkUBDO0Rq5xcPPcR1jz5Kr1//GvvWt8563uybbmLIW29x\nwSefYNGVffEtErx62dmw7MZGoEf6DqA7gY1UXgAedM7llTtvF4FCOqbC9T2AnxAYa90eyAc2AI85\n5+ZW07YeynJOXnrpJX4xcybrPvsMmzsXhgyBkhIYPRpuvx3uuOPMuXv27GHw4MGsX7+ebt26eUwt\nkUCFtEjNffTRR/y9b1++P3gw8StWVDlso6iwkA9ateLEnXcy8tFHGzClRKJ6K6R90kNZzkVhYSF9\n+vRh7vXX03/t2kAP9GnZ2ZCWBuvXQ5cuZw4/+OCD7Nixg+eee+4/byhSAyqkRWru/i98gZ8uWULT\n7Gzo3r3a89f98pe0evBBuhw7RmxidavxipxdfW0RLhK2/vWvf9Gje3f6z58PDzzw+Tf79YMZMwI9\n0uX+h/+BBx7g7bffZsuWLQ2cVkSkcdu+dSvTFiwg+uGHgyqiAYZ+//ucaNaM1XffXb/hRCqhQloi\nlnOORx55hF9PmgSlpXD55f950ve+BwcPwj/+ceZQs2bNmDFjBr/85S8bMK2IiLx3++207dSJ+Bkz\ngr/IDPfww1zw7LMU5+bWXziRSqiQloi1cOFCzIwhS5bA/fdXPs6uSRN4+mn47ncDPdaHDgFw9913\n8/rrr3Oo7M8iIlK/juzaxaUrV9Ji9myIqll5Mviuu9jTtCkbf/zjekonUjkV0hKxHnnkEX5+ww1Y\ndjbcfPPZTxw8OLDVbEEB9OkD991Hq7w8rr/+ep5++umGCywi0ogt+u//5kDHjrQbMaJW15+67TZi\ny327KNIQVEhLRNq8eTPZ2dlM2boV7r0XYmOrvqBzZ3jiicAExKgo6N+fnx87xvLf/Y6SkpKGCS0i\n0kgVFRVxcs4cWk+bVut7jHzoIdoeO8a+hQvrMJlI1VRIS0T6xz/+wT3XXEP0W2/BN78Z/IUdO8L/\n/A9s3cp5ffrwVE4OJ3v1glmz4PDh+gssItKILVy4kEudo/1XvlLreyS0aMHGYcP4SMM7pAGpkJaI\nU1JSwrPPPsvXS0pg+nRo1armNznvPHjoId783//lN+3bw7vvQq9ecMMNUFhY96FFRBqxeX/8I0lN\nmsDAged0n/Mffpg+a9dSeOJEHSUTqZoKaYk4S5YsoVNSEue99hrcddc53eumW27hf7OzOfo//wN7\n98IHH0BmZh0lFRGRY8eOEb14MdGTJ9d4kmFFF1x2GTtbtCD74YfrKJ1I1VRIS8R55plnmDlgQGCd\n6IsuOqd7NWvWjMsuu4w5c+ZA8+aBJfSWLKmjpCIiMnv2bG5p04a4L3yhTu6Xc+21mnQoDUaFtESU\n3NxcXn/9dSZ+8ME590afNm3aNGbPnh34w/jxsHhxndxXRETguWeeYchnn8Gll9bJ/QY++CDnHTpE\nwebNdXI/kaqokJaIMn/+fG7s04fYPXvg6qvr5J5TpkwhKyuLw4cPw5gxsHYt5OfXyb1FRBqzAwcO\nEL9xIzE9egQme9eB5B49WJyczJ6f/KRO7idSlaAKaQv4jpltMbN8M9tjZo+YWbWb2pvZg2ZWWsVP\nwbl/DJGAV199lXubNAls+92kSZ3cMzExkSuuuOL/hncMGACrVtXJvUVEGrO5c+dyR/fuRFW28+w5\nKP3a10h6800oKqrT+4pUFGyP9CzgUWATcA8wG7gXeC2Ia18Gplfy89uy94O5h0i1CgsLWf7mm/Td\nuBG+8Y06vfe0adN48cUXA38YP17jpEVE6sCrr75KWkFBYP5JHZpw991sLSqi4OWX6/S+IhWZc67q\nE8z6AhuBl51zN5Y7fg/wBHCLc+6FGjds9mfg68CVzrn5VZznqssoAvD222+z4Y47eGDYMPj3v+v0\n3vn5+bRv356PPvqINuvXw4MPQkZGnbYhkcfMcM5Vsjd95NIzW4L12Wef0b9zZ/aYYYcPQ3x8nd7/\nF7178/UWLTgvK6tO7yuRqzbP7GB6pG8pe51V4fiTQB6B3uUaMbMEYBqwH1hQ0+tFKvPqK69wW35+\nnU0yLC8hIYFLLrmEBQsWwOjRgS3Fc3PrvB0RkcZi3rx53NmrF5aaWudFNEDc9Ok03bgR9uyp83uL\nnBZMIT0MKAU+908651wBsAEYXot2pwEtgKfVdSF1obS0lAOzZ9OiaVNIS6uXNq688krefPNNSEyE\noUPVIy0icg5eeeUVvtisGVx2Wb3cf/LUqSyMicFpy3CpR8EU0slAjnOushH7+4F2ZhZTw3ZvJ1Cc\n/62G14lUasOGDdxeWEjcffeB1c836VOmTGH+/PmUlJRoGTwRkXNQVFTEwgUL6LVzZ52Pjz6tb9++\nbIyN5dNFi+rl/iIQXCGdCJxtZY1T5c4Jipn1AlKBd5xzu4O9TqQqK156ifGFhfCVr9RbG127dqVT\np05kZmZqwqGIyDlYvXo1aZ06Ee0c9OlTL22YGc0uuYRTK1fWy/1FAILpSc4Dks7yXny5c4L1dcAB\nfw32gpkzZ575PS0tjbR6+upewlfunDkcGzGCxJYt67WdK6+8kjfeeIPRDz4IW7bAZ59BPbcp4SM9\nPZ309HTfMbzTM1uqM3/+fL7WuTN06VJv3yIC9J8+ndavvQbFxRBT0y/PJdLVxTM7mFU75gMTgcSK\nwzvMbAVwoXOufVCNmUUD+wgU8MlnGS5S8RoNo5YqHT9+nKfbteOuH/yA2J/+tF7bysjI4K677uK9\n996DSZPgvvugjra1lcijVTtEKjds2DAWxMbSdsYMmDat3trJz89nb7NmtF+2jJapqfXWjkSG+lq1\nI6vsvBEVGosDBlFhEmI1rgbaA/8IpogWCcaSJUsY1aIFsUOG1HtbKSkp7N69m0OHDmmctIhILRw5\ncoSPPviANtnZgQ6JepSQkMD+885jxws1XqVXJCjBFNJlu1BwX4XjdwAJwHOnD5jZ+WbWu4p73U5g\nWMfTNQkpUpUFCxZwUWkp9OtX723FxMRwySWXsHjxYo2TFhGphbfffpt7+/XD+veHtm3rvT0bOpQT\nS5fWezvSOFVbSDvnNgG/B6aa2ctmdruZPUpgp8N059y/yp2+GNhc2X3MLBm4HMh0zmWfe3SRgBVv\nvUXz/Hzo0aNB2ps0aRKLFi2C4cPhww/h6NEGaVdEJBIsWLCAG8zg+usbpL3O11xDi+3bG6QtaXyC\n3SJ8BnA/0Bf4HXAj8DhwVYXzXNlPZW4tay/oSYYi1dm5cyddTpzA+vaF6OgGafN0Ie1iYgKbs6in\nQ0QkKM45lixcSK+tW+G66xqkzfOvu47ep06x96OPGqQ9aVyCKqRdwGPOuT7OuQTnXBfn3APOubwK\n5/VwzlU6LdY590vnXLRzTsM6pM4sXryYa3r2xBpgWMdpF110EcXFxXz44YcwYYKGd4iIBGnbtm2M\nKS0lqmdP6Nq1QdqMatOGE82asfb55xukPWlcgu2RFglJS5YsYVTz5tC/f4O1aWb/N7xD46RFRIK2\nZMkSvt6qFdZAwzpOy+vdm4NvvtmgbUrjoEJawpZzjvT0dC44dapBJhqWd6aQHjwY9u2DI0catH0R\nkXC0dPFiRh082GDDOk5rNXEiURs2oKUZpa6pkJawtX37dqKjo0nYubPBC+kJEyaQnp5OaVQUDBgA\nmzY1aPsiIuHGOUf+okVEd+4MF1zQoG23mTSJQSUlbNmypUHblcinQlrC1pIlS5gyejR2/HiDjbU7\nrVOnTrRs2TLwUL7oosAuhyIiclabN2/mi6WlxN58c8M3PmQIA0pLWa7J4VLHVEhL2FqyZAlXnX8+\n9O0LUQ3/V3ns2LEsX74c+vSBrVsbvH0RkXCy5J13uLq4uMGWvfuctm0pbtGC7W+91fBtS0RTIS1h\n6fT46JRmzRp0omF548aNCxTS6pEWEanW/jlziGrdOvDM9GHoUPJXrvTTtkQsFdISlrZt20ZcXBzt\nDh5s8PHRp40dO5Zly5bhLrpIPdIiIlVwztElM7PBV+sor9nYsfTNz2f37t3eMkjkUSEtYSkjI4Ox\nY8cGJvl5KqR79uxJcXExu52DTz6BEye85BARCXUfbNvG1UVFtPza17xlsGHDGNe0KcuWLfOWQSKP\nCmkJSytWrGDMmDGQne1taIeZBXqlV6yAXr1g2zYvOUREQt3W554jOjERLr7YX4jBg7kwN5flKqSl\nDqmQlrC0YsUKxvXpAwUFkJzsLceZCYcaJy0iclYxc+dyYNQoMPMXon17opo3Z+c77/jLIBFHhbSE\nnUOHDpGTk8NFJSWBYR0eH8xnJhxq5Q4Rkco5R7+tW2nucVjHaU1SUuh06BCHDx/2HUUiRFCFtAV8\nx8y2mFm+me0xs0fMLDHYhsysddk128vucdjMFptZau3jS2OUkZHB6NGjidqyxduwjtP69+/PoUOH\nONahg3qkRUQqkbNqFTHFxVzgcaLhaTZ0KFd26BDoABGpA8H2SM8CHgU2AfcAs4F7gdeCudjMugLr\ngS8D/wbuBH4O7AI61SyyNHYZGRmkpqYGxkd7mmh4WnR0NKmpqWTl5qpHWkSkErtmz2ZX+/ZERUf7\njgJDhzI8KkqFtNSZmOpOMLO+BIrnl5xzN5Y7/hHwhJnd5Jx7oZrbPEegaL/YOafvU+ScrFixgt/+\n9rfw4x/7Wdi/grFjxzL/ww+5dNcuKCqCJk18RxIRCRl5y5dTMnCg7xgBQ4bQ+cgRlmmHQ6kjwfRI\n31L2OqvC8SeBPGB6VReb2TggFfi1c+6wmcWYWUKNk4oAeXl5bNq0ieHDhoVEjzQExkkvWbUKOnWC\nnTt9xxERCSkttm+n7WWX+Y4RkJxMTEIC+du2cfz4cd9pJAIEU0gPA0qBrPIHnXMFwAZgeDXXXwE4\nYJ+ZvQ7kAyfNbJuZfanmkaUxW7NmDQMHDiThs88CkwzPO893JIYMGcK2bdsovvBCjZMWESkn78QJ\nep44wYXTpvmOcoYNGcLUbt1Ys2aN7ygSAYIppJOBHOdcUSXv7QfamVlVQ0R6A0agB7sVgXHSXwMK\ngH+a2a01iyyN2YoVK/5vfHT//n6XUioTFxfHgAED+LhlS42TFhEpZ9OcORyLiyOxUwhNh0pJ4dLm\nzVm1apXvJBIBgimkEwkUvZU5Ve6cs2le9nocGO+ce8E593dgHHAM+EUQGUSAwETDMxuxhMCwjtNG\njRrFpqIi9UiLiJRz4PXXOdqtm+8Yn5eaysXHj7N69WrfSSQCVDvZkMA46KSzvBdf7pyzyScwtONf\nzrni0wedc8fM7DXgy2bW2zl31m3hZs6ceeb3tLQ00tLSgogtkaakpIRVq1bxj3/8A954A4YM8R3p\njJEjR7I0K4sp+/b5jiIepaenk56e7juGd3pmy2lu7VpiLrnEd4zPS0mh9Z49vHv4MM45LAS+2RQ/\n6uKZbc65qk8wmw9MBBIrDu8wsxXAhc659lVc/wfgm8C3nXN/qPDeL4H/B6Q65yr9p6GZueoySuOw\nceNGrr/+erZt2wajR8OvfgXjxvmOBcDevXuZNGQIWwsLsWPHQmLIifhnZjjnGtVfBj2z5bTS0lLW\nxMbS68UXaXPddb7jfN6QIUz9+GN+tWwZvXr18p1GQkRtntnBDO3IKjtvRIXG4oBBVJiEWIk1BMZI\nd67kvS5lr1oST6qVmZlJSkoKOBdyQzs6d+5MbmwsJU2awIEDvuOIiHj3wZYt9C8tpc3Eib6j/KfU\nVKa2b6/hHXLOgimkXyx7va/C8TuABAJrRANgZuebWe8K570KnACml98J0cw6AtcAHzjntGaYVOtM\nIb1vHyQmQtu2viOdYWaMGjWKo0lJmnAoIgJsnTuXE02bQqtWvqP8p9RURpWWasKhnLNqC2nn3Cbg\n98BUM3vZzG43s0cJ7HSY7pz7V7nTFwObK1x/DLifwA6GmWVbjX8fWAU0IbDZi0i11qxZw4gRI0Ku\nN/q0kSNHsj06WhMORUSAz955h8969vQdo3KjR9N1/35Wq5CWcxTsFuEzCBTDfYHfATcCjwNXVTjP\nlf18/qBzTwLXEeiZfgj4AbAFSHPOvVOr5NKo5ObmsmPHDgYOHAhZWTB4sO9I/2HUqFGsOX5cPdIi\nIkDspk3EjR7tO0blunYlJjGRkm3bOHnypO80EsaCKqRdwGPOuT7OuQTnXBfn3APOubwK5/VwzlW6\nEohz7lXn3GjnXHPnXEvn3BVnm2AoUtG6desYMGAAsbGxsGwZhNoscAIbsyw7fJiS7GzfUUREvMrP\nz6d7Tg7JV1XsbwsdNmYMNyQnk5VV3VQvkbMLtkdaxKsz46OLimD1akhN9R3pPyQkJFDSqxfFmzb5\njiIi4tW7WVkMBOJGjfId5exSU5mUmKgJh3JOVEhLWDhTSK9bBxdcAK1b+45UqfMvuYSoY8fg+HHf\nUUREvNnxxhucaN4cWrb0HeXsRo+m77FjmnAo50SFtISFM4X00qUhOazjtJRRo9ibmAjbzrq/kIhI\nxMtdupS8iy7yHaNqAwfS/NNP2bpyJVr7XGpLhbSEvP3791NQUECPHj0C46NDZBOWyowaNYr3Cgpw\nmzdXf7KISIRK3LKFpiHc6QFATAyWkkJKaSm7du3ynUbClAppCXmne6OttBQyMkK6kO7evTvboqM5\nrskrItJIHT58mL55ebS7/HLfUaplo0dzTbt2GicttaZCWkLemWEd770HycmQlOQ70lmZGVF9+nBi\nzRrfUUREvFizciUXA1HDhvmOUr3UVEYUF6uQllpTIS0h70whHeLDOk5rk5pKkx07fMcQEfFi17x5\nnGzZElq08B2leqNGkXzgAGsyMnwnkTClQlpCWklJCevWrWP48OEhu350Rb2uvJJWx44FluoTEWlk\nTmVkcKp/f98xgtOyJXbBBcRmZ5Ofn+87jYQhFdIS0jZv3kxycjKtW7YMmx7pIaNGsRco1IRDEWlk\nSktLabljBy0nTvQdJWhRqalck5TE+vXrfUeRMKRCWkLamWEdW7ZAq1bQqZPvSNVq1qwZO1q0YN+c\nOb6jiIg0qA8++IBhZjQPg28Pz0hNZXxcHJmZmb6TSBgKqpC2gO+Y2RYzyzezPWb2iJklBnl96Vl+\ntGuFVCncxkef9km/fuQvXOg7hohIg1qzciV9i4th8GDfUYKXmspFR4+SqQmHUgvB9kjPAh4FNgH3\nALOBe4HXatDWMmB6hZ/ba3C9NEKf24gljArphEsvpbWGdohII7N3wQLy2rYNj4mGp/XoQWxMDPs1\n4VBqodpC2sz6EiieX3LO3eCce8o5dz/wXWCCmd0UZFs7nXPPV/j59zlklwiXm5vLzp07GXDxxWEz\n0fC0XlOn0iw3Fw4f9h1FRKTBlK5eTXE49UYDmBE9bhz9jh3j4MGDvtNImAmmR/qWstdZFY4/CeQR\n6FkOipk1MbOmwZ4vjdvatWsZMGAAsXv3QnQ0dO/uO1LQLurXj9VRUZyYP993FBGRBpGfn0+n/ftp\nddllvqPUmI0dy9Vt2mictNRYMIX0MKAU+NxWbc65AmADMDzItq4nUHifMLNDZvaEmYXRdz/S0P5j\nWIeZ70hBi46OZk/XruS88orvKCIiDWL9+vWkxsYSO2aM7yg1N3Ysw0+dUiEtNRZMIZ0M5DjnKlsU\ndz/QzsxiqrlHJvAgcB3wFeAdAsNFlgU7YVEan89NNAyjYR2nlY4eTRPtcCgijcS6ZcvoUVQEAwf6\njlJzAwfSOi+PbcuX+04iYSaYQjoRKDjLe6fKnXNWzrlRzrnHnHOvOeeedc7dAvwIGADMCDqtNCqZ\nmZmMGDEi7FbsOK3j1VfT9tAhyM31HUVEpN7lvP02x7t0gbg431FqLiYGN3IkCevWUVJS4juNhJFg\nCuk84Gz/VcSXO6emfgsUAlfW4lqJcPv27aOwsJAe0dGQlwe9e/uOVGPDx47lPcBpSSURaQSabNhA\n9KhRvmPUWuykSUyMjWXr1q2+o0gYqW5IBsDHQB8za1LJ8I5OBIZ9FNe0YedcsZl9DLSr7tyZM2ee\n+T0tLY20tLSaNidh5vSwDlu+POzGR5/WoUMH5jRtyoVz59J20iTfcaQBpKenk56e7juGd3pmNz6H\nDh2iX24uLS+91HeU2hs3jrTf/pYlmZn069fPdxppAHXxzDbnXNUnmP0M+CEwzjmXUe54HHAUSHfO\nfaHGDQeuPwGscs6ddQCsmbnqMkrk+d73vkezZs347z17YMAA+Pa3fUeqlV+OHcttn3xCx+xs31HE\nAzPDORd+/wo8B3pmN06vv/46w2+8kQ4bNoTlN4gAFBRQ2KIFD9xyC4//7W++04gHtXlmBzO048Wy\n1/sqHL8HK7TLAAAgAElEQVQDSACeKxfgfDP73H9BZtbmLPd9GIimZpu6SCPxuRU7wnCi4WnNL7uM\n1tu3Q1Flc3VFRCLDpsWLaekcXHih7yi1FxfHqf79KV661HcSCSPVFtLOuU3A74GpZvaymd1uZo8S\n2Okw3Tn3r3KnLwYqbuf2YzNbaWY/N7Nvmtl/mdk7wH8Bq4Hf1c1HkUhRUlLCunXrSOnaFY4ehf79\nfUeqtUHjx7M3JgY2bPAdRUSk3uSmp5Pbpw9EBbthcmhKnDyZHnv3cvLkSd9RJEwE+zd+BnA/0JdA\n4Xsj8DhwVYXzXNlPeenAZwSWvXsMmAm0Bn4AjC9bj1rkjOzsbDp16kTLDRtg7NiwfjAPGTKEJYWF\nFC1e7DuKiEi9KC0tpcWWLSRGwFj4mAkTuCw+nrVr1/qOImEiqArFBTzmnOvjnEtwznVxzj3gnMur\ncF4P51xMhWOvOeeuKLsm0TnX3Dk3xDn3a+dcYV1+GIkMkTKsAyAxMZGPunTh+Ftv+Y4iIlIvtm3b\nRkpUFE0joJBm5Eh6nzrFOq0nLUEK364+iViRVEhDYOvZhHXrQBOwRCQCZa5ezdCSEhgxwneUc9e0\nKSe6deOzhQt9J5EwoUJaQk5mZiapF14IBw6E5w5ZFfSaOJGTpaXwwQe+o4iI1Lmdb78NiYnQsaPv\nKHUiesIEmmteiwRJhbSElBMnTrBz50765uRAaipER/uOdM5SUlLIMAN9VSgiEah45UoKBw3yHaPO\ntLrqKoafOsX+/ft9R5EwoEJaQsratWsZOHAgTTIyImJYB0CvXr1YWlpK/qJFvqOIiNSpvLw8kvft\no0UEbTplY8YwvLSUNStW+I4iYUCFtISUSBsfDRAVFUXuoEGUam1SEYkw69evZ2x8PE1SU31HqTut\nW3M8KYmP33jDdxIJAyqkJaRkZmYytk8f2L0bhgzxHafOdJwwAfvsM/j4Y99RRETqTNbKlVxUUABD\nh/qOUqcKUlKIzsio/kRp9FRIS0hZs2YNqaWlMHo0xMRUf0GYGDFyJO81bQr6qlBEIsihRYvIa98e\nmjf3HaVOtfviF+m+Zw/FxcW+o0iIUyEtIWPfvn0UFRVx3pYtETOs47SUlBTm5ebili3zHUVEpM7E\nrFtHVEqK7xh1runkyYx2juz33/cdRUKcCmkJGatXr2bkyJHYsmURV0gnJSWR3bo1BZpwKCIR4uDB\ng/TNzaX5xIm+o9S99u3JbdaMHa+84juJhDgV0hIyMjMzuWTAANixA4YN8x2nzjUdO5bo3bvh2DHf\nUUREzllmZiapTZoQNWqU7yj14tOLL6bg7bd9x5AQp0JaQsbq1auZFB8PKSkQG+s7Tp0bNno0O9u0\ngZUrfUcRETln7y5dSsfCQujf33eUepF4+eUkbd7sO4aEuKAKaQv4jpltMbN8M9tjZo+YWWJNGzSz\nBDPbZWalZvZEzSNLJCoqKuLdd9/lokOHIm5Yx2kpKSmkl5RoYxYRiQjHFy/mZM+e0KSJ7yj1ovOX\nvsSgEyc4/tlnvqNICAu2R3oW8CiwCbgHmA3cC7xWizZ/BrQBXC2ulQi1ceNGunfvTtzq1RFbSA8a\nNIjXjx2jROtJi0iYKykpoe2WLcRH0EYsFTU5/3wK4+LY8vLLvqNICKu2kDazvgSK55ecczc4555y\nzt0PfBeYYGY3BduYmQ0BZgAPAlbLzBKBVq9ezSWDB8OWLTBihO849SI+Pp4T/frBhg1w6pTvOCIi\ntZadnc2E6GgSLr/cd5R6tb9nT47Nnes7hoSwYHqkbyl7nVXh+JNAHjA9mIbMLKrsmnmApsHK52Rm\nZvKF1q0Dkwzj433HqTcDRo/mcLt2sGaN7ygiIrWWlZHBoKIiiKQdDSszbhxN1671nUJCWDCF9DCg\nFMgqf9A5VwBsAIYH2dZ3gV4EerdFPmf16tUMPXkyYod1nJaSksLa+HiNkxaRsHZo3jxyO3SAVq18\nR6lXnW65hV4HD+JKS31HkRAVTCGdDOQ454oqeW8/0M7MqtyCzsx6ADOBnzrn9tY4pUS0Tz75hAMH\nDpC0bRuMHes7Tr1KSUlh7iefqJAWkbAWv2YNpWPG+I5R75JHj6YQ+Dg93XcUCVHBFNKJQMFZ3jtV\n7pyq/BH4EHgsyFzSiKxZs4aRgwZhGzYElr6LYD179mRpSQmlK1dCSYnvOCIiNXb8+HH6Hj1K22uv\n9R2l3llUFNuTkznwwgu+o0iIqrInuUwekHSW9+LLnVMpM5sOTALGOudqVTnMnDnzzO9paWmkpaXV\n5jYSojIzM5narRvk5UHz5r7j1Cszo+fIkeRu2kSL996DIUN8R5I6lJ6eTrp6rvTMjnBZq1eTCsRM\nmOA7SoPIGzaMlvoWMSLVxTPbnKt6FTozmw9MBBIrDu8wsxXAhc659me5NhbYC2QC3yn3VmdgCfBP\n4CECQ0cqXajRzFx1GSW8XXHFFTySlES/Nm1gVsU5rZFn5syZXDF3Lim33QYzZviOI/XIzHDONaoV\nivTMjnxP3X03X3jhBdofPeo7SoNY+cwz9LzjDs47dQqsUf3n3OjU5pkdzNCOrLLzPrcmmZnFAYOo\nMAmxggQCvdlXAtvL/SwhsI70l4EPgNtrEloih3OONWvWcMGBA9AIxtsBjBw5kncKCjROWkTCklu6\nlNyhQ33HaDD9r72W4qIiirZt8x1FQlAwhfSLZa/3VTh+B4FC+bnTB8zsfDPrXe6ck8D1wA1lr6d/\n7iSwjvRbZX+uzcYuEgG2b99O86ZNiV+/PvKXUSozYsQIntuzB7d8OajnTkTCiHOOjtu30+qqq3xH\naTAtWrbk3ebN+fhf//IdRUJQtWOknXObzOz3wN1m9jKBdaD7At8G0p1z5f9mLQa6ANFl1xYDcyre\n08y6lf36oXNOa0o3YqtXr2Zq376wfTt07Og7ToNo06YNRR07UnTiBLE7dsCFF/qOJCISlF07dzKq\nuJjW11zjO0qDOtq/P6fmz4ef/tR3FAkxwW4RPgO4n0AB/TvgRuBxoOI/SR3Bb/1dk3MlQmVmZjKl\nRYtGM6zjtJSRI9nTrZuGd4hIWNkyZw6F8fFY166+ozSouMsvp212tu8YEoKCKqRdwGPOuT7OuQTn\nXBfn3APOubwK5/VwzgXTy73bORftnNNMq0Zu9erVDDhxovEV0ikprI6JUSEtImHl5Lx5HOnTx3eM\nBnfR1VfjTp2Cjz7yHUVCTLA90iJ1Li8vjy1bttBu69ZGWUjPOXJEhbSIhJVWGzfSZNIk3zEaXL/+\n/Vluxsl583xHkRCjQlq8ycrKYkLv3kSdOAEXXeQ7ToMaOHAgC/ftw33yCRw44DuOiEi18vPy6P/J\nJ3SbPt13lAYXExPD7u7dOTZ3ru8oEmJUSIs3K1eu5Ibk5MBqHY1sbc7Y2FguHjiQo336wNKlvuOI\niFRr02uvERMTQ0K/fr6jeFE6ZgwJWVWt+CuNkQpp8SYjI4PRzjW6YR2npaSksKFtW1i0yHcUEZFq\nHX7pJfb16NHoOj5O6z5lCtG5ubBvn+8oEkJUSIsXpaWlrFy5ku779jXqQvrV/HxYuFDrSYtIyIvL\nzMSNG+c7hjcpo0ax3Ax3jltKS2RRIS1ebN26lU4tWtBk504YMsR3HC9GjRrFS++/jzMD7ZglIiGs\ntLSUCz7+mC633OI7ijedO3cmKyGB46+/7juKhBAV0uLFypUrmd6zZ6CIjovzHceLbt26EdOkCSdG\njgz0SouIhKidy5bRAjgvLc13FK8KR44M7EorUkaFtHiRkZHBxLi4RjusA8DMGDNmDOvbtVMhLSIh\nbe/zz7MzObnRjo8+rcuVVxJ79CgcPuw7ioQIFdLiRUZGBr1zchp1IQ0wZswYXjl+HJYtg4IC33FE\nRCoVtXQpBSNH+o7hXeq4cWQ1aQIrVviOIiFChbQ0uCNHjvDJoUM027wZRo3yHcerMWPGsHDtWujT\nB1at8h1HRKRSXXftIumGG3zH8K5///6kFxeTp28RpYwKaWlwK1eu5Ja+fbEePaB1a99xvLr44ov5\n+OOPOTlmjIZ3iEhIOrppE62Kiuj5xS/6juJddHQ0JwYNokDLlkqZoAppC/iOmW0xs3wz22Nmj5hZ\nYhDX9jKzZ81ss5kdM7OTZfd51Mw6nPtHkHCTkZHBla1aNfphHRB4KI8aNYr1bduqkBaRkLTr739n\nS1IS0U2a+I4SEtpdcQWJe/bA8eO+o0gICLZHehbwKLAJuAeYDdwLvBbEtZ2BDsAc4PvADGAhcAew\n1sza1TCzhLmMjAyGHD8OY8f6jhISxowZwxs5ObB9Oxw54juOiMjnFC9axIlGukxpZUalpZEdH6/h\neAIEUUibWV8CxfNLzrkbnHNPOefuB74LTDCzm6q63jm32Dk3yTn3Y+fcn5xzf3XOzQC+CiQDt537\nx5BwUVBQwMZ336VddjZMmOA7TkgYM2YMy1atgrQ0eOcd33FERD6nwwcf0Oraa33HCBnDhw/n7fx8\nihYv9h1FQkAwPdKnV1+fVeH4k0AeML2Wbe8pe23cg2QbmczMTKZ264Z16wbt2/uOExJGjBjB+++/\nT2FamoZ3iEhIyd+1i5b5+fS7+WbfUUJGYmIiB3r25OT8+b6jSAgIppAeBpQCWeUPOucKgA3A8GAa\nMrM4M2trZp3M7DLgT4AD5tUssoSzpUuXMq1dO5g40XeUkJGYmMjFF1/Mu0lJ2i5cRELKzr/9jU2t\nWtGsZUvfUUJKs8suI3HrVi1bKkEV0slAjnOuqJL39gPtzCwmiPt8HTgC7AXmAy2B6c65jGDDSvhb\nunQpwz/7DCZN8h0lpIwZM4a3d+2CJk1gyxbfcUREAMh76y0+GzzYd4yQM3zCBD6Kj4esrOpPlogW\nTCGdCJztn1ynyp1TnVeAScC1wE+BY0BSENdJhCgsLGTj6tW03bkTxo3zHSekjBkzhuUrVsBll2l4\nh4iEjKTNm2mtZe/+w+jRo3k7P5+SpUt9RxHPgulJzuPsBW98uXOq5Jz7GPi47I+vmdkcIMvMEpxz\nv67q2pkzZ575PS0tjbS0tOqakxCUlZXF9R06YMnJ0KyZ7zghZezYsdx6660Uf+1rxDzzDNx3n+9I\nUgvp6emkp6f7juGdntmR4dSePbTMy6PN9NpOhYpcSUlJfNC+PblvvknLH/3Idxyppbp4ZpurZjym\nmc0HJgKJFYd3mNkK4ELnXK1mjZnZKiDZOdetinNcdRklPPziF79g1Jw5jL/6avjJT3zHCTkDBw7k\nyd/8hhE33BBYBi8uznckOUdmhnPOfOdoSHpmR47NP/0pRx97jLHHjvmOEpK+f/vtzHz+eeJzcyE6\n2nccqQO1eWYHM7Qjq+y8ERUaiwMGUWESYg0lAG3O4XoJI0uXLmXwp59qfPRZjB8/nkXr1kG/fpCh\nqQMi4tfJefM4NnCg7xgha/iUKRyOjob33/cdRTwKppB+sey14nfNdxAohJ87fcDMzjez3uVPMrNK\ne6vNbDzQH9CK5o1AUVER2zIyaHHkCAwPaqGXRmfChAksWbIkME56wQLfcUSkkWuXna31o6twySWX\n8E5hISUaztWoVTu0A8DMngDuBl4lsFxdX+DbwHLn3MRy530EdHHORZc7NgfoCCwGdhMYVz0UuAnI\nBdKccxuraFtfE0aA1atXM+emm/jNxRfD66/7jhOSjh07RpcuXTg6dy6x998P69f7jiTnSEM7JFwV\n7t9PXufOcOQIrdppA+Kz+VHXrnz3ggtou2SJ7yhSB+praAcEtvW+n0AB/TvgRuBx4KoK57myn/Ke\nB3IIbNwyC/glgbWn/wgMrKqIlsixdOlSrm3WTMM6qtCqVSt69+7NGjP48ENtFy4i3uz8+9/Z2KKF\niuhqxE6aRHxWltb/b8SCKqRdwGPOuT7OuQTnXBfn3APOubwK5/VwzsVUOPaSc+4q51w351yic66p\nc66vc+4+59y+uvwwErqWLl3KwCNHVEhXY/z48SxevlzbhYuIV7lvvMGnAwb4jhHyBl1zDSeLi2H7\ndt9RxJNge6RFaq24uJj9K1aQANC3r+84IW38+PGBcdKXXqr1pEXEmzabNtHymmt8xwh5l1xyCUtK\nSihSx0ejpUJa6t3atWuZ2qIFUZdeCtaohovW2NixY8nKyqJg3Dh4+219XSgiDa5g/37a5uYy4NZb\nfUcJea1ateLD5GQ+mTvXdxTxRIW01LuFCxdyjcZHB6V58+ZcfPHFrMzJCfyjY9s235FEpJH54K9/\nJbtlS1onafPhYERddhlNMzKgtNR3FPFAhbTUu0ULF9L3wAGYOLH6kyUwvCM9XduFi4gXJ+fO5TMt\nUxq0AV/8Ip+WlMB77/mOIh6okJZ6dfz4cUrWryemQwfo0sV3nLAwfvx4Fi9eHBgn/fbbvuOISGPi\nHF2zsznvK1/xnSRsjB07lteLiih67TXfUcQDFdJSr9LT07m1U6fA+GgJypgxY3j//fc5Pnw4LFsG\nRUW+I4lII/HJypUUFxUxYNo031HCRvPmzdnZqxcnXnrJdxTxQIW01KuFCxcyubQ0MExBgpKQkMDo\n0aNZtGEDXHghrF7tO5KINBK7//Qnsrt2pUlsrO8oYSXphhtI3LYNjh3zHUUamAppqVcb580jOScH\nLr/cd5SwMnnyZObPn69l8ESkQcW98w4l+gaxxiZddRVZsbFa/78RUiEt9Wb37t1MOHyY6Ouvh7g4\n33HCyulC2k2apHHSItIg3IkTdD14kF7f+pbvKGFn8ODBzAdOzJ7tO4o0MBXSUm/efvttbo2NxaZP\n9x0l7PTu3ZuoqCi2tm0LmzfDp5/6jiQiEW7/P//Je7GxXDhkiO8oYScqKorCiROxhQu1/n8jo0Ja\n6k32Sy+R5ByMG+c7StgxMyZPnsxbixdDaiosXuw7kohEuE+ef579AwZg2jirVgbdcAO5hYWwaZPv\nKNKAgiqkLeA7ZrbFzPLNbI+ZPWJmiUFce6GZPWRmq8zssJkdN7N3zeyHwVwv4amkpITOy5ZROm0a\nREf7jhOWPjdOWsM7RKQ+OUf79etpfuONvpOErcsuv5zXi4spfv1131GkAQXbIz0LeBTYBNwDzAbu\nBYJZNPFrwAxgB/BT4H5gK/AwkGFmGjwbgdZkZjKtuJjmd9zhO0rYmjBhAqtWrSIvNVWFtIjUq1Mb\nNlB46hQjbrvNd5SwlZSUxJauXcn99799R5EGVG0hbWZ9CRTPLznnbnDOPeWcux/4LjDBzG6q5hb/\nBjo7577snPu9c+4vzrmbgZ8DA4Dbz/EzSAh6789/JrZFCxg82HeUsNWiRQuGDh1Kek4O5OXBhx/6\njiQiEWrXH/7A+vbtaduune8oYa311KkkZGfDiRO+o0gDCaZH+pay11kVjj8J5AFVziRzzq13zlX2\nN+pFwID+QWSQMNPyzTc59cUvgsbanZPJkyczf8ECDe8Qkfo1f76WvasDE66+mg1xcVoGrxEJppAe\nBpQCWeUPOucKgA3A8Fq2fXq/6EO1vF5C1O4dO5j4ySd0euAB31HC3uTJk3nzzTdxl18Ob73lO46I\nRCB34gRd9u+n77e/7TtK2EtJSeGNkhJOvvyy7yjSQIIppJOBHOdcZfsU7wfamVlMTRo1syjgJ0AR\n8HxNrpXQ9/5jj3G8bVuie/XyHSXsDRw4kOLiYrZ27w5LlkB+vu9IIhJhPvrb33g/Lo7ew4b5jhL2\nYmJiKJg4kdJ587QMXiMRTCGdCBSc5b1T5c6piceBEcB/O+e21/BaCXFNX32VE1dd5TtGRDAzrrnm\nGuakp8OgQZCe7juSiESYT557jsNDh2rZuzoy/MtfJj8vD7Zs8R1FGkAwPcl5QNJZ3osvd05QzOxn\nwN3An5xzvwnmmpkzZ575PS0tjbS0tGCbkwaWe/gwQz7+mKgf/MB3lIhx7bXX8r3vfY8fXX89vPkm\nXHGF70hyFunp6aTrHzt6ZocT5+j43nsUPPGE7yQRY/IVV/DvkhJuefVV4vv29R1HqlAXz2xz1Xz1\nYGbzgYlAYsXhHWa2ArjQOdc+qMbMZhIY0vGUc+4bQV7jqssooWP1ffcR8+yzDMvJ8R0lYhQVFdGh\nQweyZ8+mw+23w65dmsQZJswM51yj+n+Wntnh5dDSpRSPH09Sfj6xcVqNtq78bOhQ7igooL02Zwkr\ntXlmBzO0I6vsvBEVGosDBlFhEmIV4R4kUET/PdgiWsJP7OzZ5F59te8YEaVJkyZMmTKFOVu3Bgro\n7GzfkUQkQmyfNYsPLrhARXQdS/7qV2n2wQdw8KDvKFLPgimkXyx7va/C8TuABOC50wfM7Hwz613x\nBmb2E+BB4Bnn3NdqmVVCXMG2bXQ7eJB+Dz7oO0rEueaaa3h17ly48kp44w3fcUQkEjhH54ULaaJN\nWOrclOuu43Uzip/XegqRrtqhHQBm9gSBcc2vAvOAvsC3geXOuYnlzvsI6OKciy537G7gf4HdBHqk\nSyvc/pBzblEVbetrwjCx/brr2LBmDTfs3es7SsTJzc2lU6dO7PnLX2j5u9/B8uW+I0kQNLRDQtmx\nuXM5MHUq3U+cICGxpmsGSHXu79ePH5WU0HrrVt9RJEj1NbQDAlt830+ggP4dcCOBlTcqLs3gyn7K\nG1Z2rCvwd+AfFX5+WJPAEqJOnaL9vHmc+upXfSeJSM2aNePyyy/npZwceP99OHrUdyQRCXNHH36Y\n1QMHqoiuJz2+8Q2iP/oIdu70HUXqUVCFtAt4zDnXxzmX4Jzr4px7wDmXV+G8Hs65mArHvuqci67i\nZ0JdfiDxo/j558kqKWHCN7/pO0rEmjZtGv965RVIS4P5833HEZFwdvAg523YQKsZM3wniVhTp03j\n385R9M9/+o4i9SjYHmmRKp341a9Y2LMnnTp18h0lYk2ZMoW1a9dyfNy4wDJ4IiK1lPv448wx47Lr\nr/cdJWJ17NiRDX36kP/0076jSD1SIS3nbu1aivbv5wJtL1uvEhISuPLKK3mloCDQI11c7DuSiISj\nkhJK//xntk+aRNOmTX2niWh9v/ENCnNyYONG31GknqiQlnNWOGsWfygp4fpp03xHiXjTpk3jqfnz\noUcPWLnSdxwRCUfz5rG3sJBUdX7Uu+tuuIHniospfOYZ31GknqiQlnNz9CjMmcOHaWm0adPGd5qI\nN3nyZLZs2cKno0dreIeI1MrJRx/lz1FRXHrppb6jRLzzzjuP7cOGBQpprWYTkVRIy7l5+mmWtmrF\nVV/T8uANITY2lptvvpmXTp3SetIiUnO7dmFZWTSZPp2YmJjqz5dzNu7eezmalwerV/uOIvVAhbTU\nXkkJxb/7Hb/JzeWqqyquhCj15dZbb+VXixbhcnIC24WLiATJ/elPvNCkCTfeeqvvKI3G1ddcw/PO\nceLPf/YdReqBCmmpvfnzOVxcTM9bbiEhIcF3mkZjyJAhJDZrxqEhQ9QrLSLBKyig6MkneTkpiREj\nRvhO02jEx8dz6tprsZdf1iTxCKRCWmrN/f73PFZQwDfuuMN3lEbFzLjtttt4trgY/vIXjbsTkeC8\n/DIfxMVx2T33YNaoNtz07sr77mNnURFu8WLfUaSOBbVFuE/abjZE7d5N4YABXNKjB6s2bPCdptE5\ndOgQvXv14mjXrkT/6ldw5ZW+I0kltEW4hJKiUaO4/f33mbV3ryaHNzDnHL/t2JEvDxpER22oFbLq\nc4twkc976ikWnXceX9FOhl60b9+eK6ZMYeGQIfCrX/mOIyKh7v33KdiyhdKrrlIR7YGZkTRjBs0X\nL4bDh33HkTqkQlpqrriYkief5GcHD3LzzTf7TtNo3XXXXfzXqlW4jz+GFSt8xxGREOb+8Af+HhvL\n7er88Oa6u+/mxagoTjz0kO8oUoeCKqQt4DtmtsXM8s1sj5k9YmaJQV7/AzObbWYfmlmpme08t9ji\n1Ztvsj8mhkHTp9OqVSvfaRqtMWPGEBMfz9YvfAF+/WvfcUQkVB0/TvFzzzGnTRvS0tJ8p2m0WrRo\nwUc33ED000/Dp5/6jiN1JNge6VnAo8Am4B5gNnAv8FqQ1/8cGA/sAPS3J8yV/OlP/Pazz5gxY4bv\nKI2amXHXXXfx4K5dsHYtbNrkO5KIhKJ//pPVzZrxpQce0CRDz6Z973u8DhTPmuU7itSRaicbmllf\nYCPwsnPuxnLH7wGeAG5xzr1QzT26O+c+Kvt9I9DUOXd+UAE1cSW07NnDqb59uTk1lVcWLPCdptHL\nzc2lR48ebL71VpIOH4Z//MN3JClHkw3FO+co6NWLm3JyeP7jj7VUaQi4bdQo/pSdTfz+/dC8ue84\nUk59TTa8pey14j+fngTygOnV3eB0ES3hz/31r8yJj+eu++/3HUWAZs2aceedd/JwTk5gy/Ddu31H\nEpFQsmwZnxw9ysX33KMiOkR86aGHWGxG6R//6DuK1IFgeqTnAxOBROdcUYX3VgAXOufaB92geqTD\nV3Ex+R06ML1NG17atk1fEYaII0eO0KtXL/Z+6Us0i4qCJ57wHUnKqEdafMu/+mpmLlrEd3bupEOH\nDr7jCIGl8G7q25dnDh0K9ErrHzgho756pJOBnIpFdJn9QDszi6lJoxKe3Lx57Dh1ipt/8QsV0SEk\nKSmJ6dOn8zjAs8/CkSO+I4lIKDh4ELdwIXzlKyqiQ4iZMXXmTLIAnnrKdxw5R8H0SO8AYpxz3St5\n7xkCQztaO+eOB9WgeqTD1uGUFB7ft4+f7d1LVJRWTgwlu3fvZujQoey98koSunSBhx/2HUlQj7T4\nlfv97/PSY49x2a5dJCcn+44j5RQXF3PT+efz7KlTxO/bB7GxviMJtXtmB9OTnAckneW9+HLn1JuZ\nM2ee+T0tLU3L93jg9uwhbv16Bv31ryqiQ1C3bt246aabePTUKX78xz/C3XdDx46+YzU66enppKen\n+47hnZ7ZIaC4mOI//IED116rIjoExcTEcP1vfsOGb36TlGeewb7xDd+RGqW6eGZrjLQEZdvNN7Nu\nwYOahFAAABeYSURBVAKmHTlCdHS07zhSiQMHDtC/f38+uvFGmhcXw5NP+o7U6KlHWnzJ+etf+fDO\nO0neuZMuXbr4jvP/27v36Kjqc//j72dyGUggIYCgKMhFQBCL1xzxQk+l1B7By7G1pRYqggKKCliV\nn9oWKVQ9iheQWpeKRQsV7x6wXBQsFVrwgJeKi4sUIlQBJRJygSRMZr6/P/ZgQxogGSbZM5nPa61Z\nyeyZPXx2sufhyZ7v/m6pRSQS4foePXi8rIzsggKNlU4ADTVGek30efk1/rEgcEb0cWnCQsXFtHr5\nZTpNnaomOoGdcMIJ3HDDDfxi/36YPx/WrfM7koj44cAByiZOpGDwYDXRCSwQCPCD6dNZWV5ORMPx\nklZdjkj3Af4OvOacu7ra8lvwpsQb6px7IbqsK5DhnNt0hNfTEekks2rwYEJr1nDRrl06yTDBFRUV\n0atXL9b+7GectG4dLFrkd6SUpiPS4oddI0bw8R//SL+vvqJlTo7fceQInHNcfcEF/GHdOpq/9x70\n7u13pJTWIEeknXOfAL8FrjKzV81spJk9jHelw+UHm+iod4D1tQQbamb3mNkv8MZb50bv32NmR52H\nWvyzZ/16eixcSLtnn1UTnQTy8vKYOnUq17z7Lu4f/4C33vI7kog0IrdiBRlz57JryhQ10UnAzLj3\nqae41zlCI0ZAJOJ3JKmnox6RBjCvgxoPjAI6A4XAPGCSc25/tecVAB2dc+k11v8z0P8wL/8X59zF\nR/i3dXTDR3/u2ROys/nOBx/4HUXqKBKJcN555/FAfj4Xr1gBH3wAGpLjCx2RlkZVUkJJt25MbtWK\nBzdu1FC8JHLbuHHc/MILdL3/fhg50u84KSuWml2nRtpPKsr+WfX003QfM4ZgQQEtO3XyO47Uw9q1\naxk8aBDbu3Qhc/RouO46vyOlJDXS0pj2/+QnvPrGG3xr9Wr69u3rdxyph+LiYq7u3p03QyEyN22C\ndu38jpSSGupkQ0lBZWVlHBg3ji9HjlQTnYTOOeccRowcyT2Zmbhf/hL27fM7kog0IPfqq5QsWMA/\nbrpJTXQSys3NZfzs2cyORDhwyy1+x5F60BFpqdVDAwfys/feo31hoSaKT1KVlZXk5+fzRvPmdBk4\nEKZM8TtSytERaWkUO3ey/9RTGdW2LU9/8gnNNY1a0ho7fDiTX32Vtq+9BgMH+h0n5eiItMTF3Oef\n5/J33yX3ySfVRCexYDDI888/z5WffsqBZ5+F6dP9jiQi8VZczL7LLuN3oRD3vPmmmugkd/+MGUxs\n0YLSoUOhqMjvOFIHaqTlEOvWrWPt2LGceNppNPvJT/yOI8eob9++THjkES4JBgnPmAH33ed3JBGJ\nl88/J3z++czfupW86dPp1auX34nkGOXk5HDr4sXMKSmh7Hvfg8pKvyPJUaiRlm/s2rWL4Zdeym/S\n02nx5JOg6e6ahOHDh3PqJZdwXdeuuLlz4Z57QB+9iyS3detw55/P76uqWDlkCNddf73fiSRO+vbt\nS85TT7Fy/XoqrrlGU+IlODXSAngnF/735ZfzckYGWTfeCPn5R19Jksb06dPZZcaEM87ALVoEEyao\nmRZJVsuW4QYM4Pc9evBKly5MnzFD8/w3MT8dNoz3brmFT5csoXLiRL/jyBHoZEOhrKyMSy+9lDtL\nShjUpg22ZAmkpx99RUkqpaWlDBgwgO+fdx6T167FTjsNnnxSc0w3IJ1sKHE3Zw7uttv4bf/+PLdt\nG0uXLiU3N9fvVNIAnHPcfu21THjlFY576CGCY8f6HanJ0zzSUm/FxcVcccUVXJmRwbiNG7H339f8\nlU1YYWEhAwcO5OL8fKZt3owdfzw89xxkZPgdrUlSIy1x4xxMm4Z7/HEe/Pa3eWn9epYuXUpeXp7f\nyaQBhcNhJl51FXctWkSzF14g+wc/8DtSk6ZZO6Retm/fzoUXXshFJ53EuI8/xubNUxPdxLVt25bl\ny5ezZsMGhrdtS7ioCK6+Gioq/I4mIocTicBttxGZPZsxp5/Om599xttvv60mOgWkpaXx4Ouv89zl\nl1M5ZAi7Z8/2O5LUoEY6Rb3zzjv069ePG4YN49cbNmB33QUXXOB3LGkEubm5LF68mFB6Ov127qQs\nFILLLtNFW0QSUWUl/PSnlK9cycBgkLJWrVi6dCmtW7f2O5k0kkAgwISXX2bh6NFUjRzJltGjdY5L\nAqlzI22eCWa2wczKzWy7mU0zs6zGWF/io6Kigrvvvpthw4Yx+5lnuHXzZqxbNxg3zu9o0oiysrKY\nO3cu1wwfTrf33mNzeTnukkuguNjvaCJyUEkJbtAgtm/eTI+CAi6/9lrmzJlDMBj0O5k0MjNj6MyZ\nbJk7l33PPssHffpQsXev37GE+h2Rfgx4GPgEuBl4CbgVmN9I68sxWrJkCaeffjr//PhjNl1/PQPH\njIGtW2HWLE11l4LMjPHjx/PWsmX8uKyM+du3U3HuufDWWzraIeK3FSuoOOccFmzcyJWhEPPffptx\n48Zpdo4Ud+GQIRy/eTPlJSVsOuEE3pkzB52T4K86nWxoZr2BdcCrzrkfVVt+MzADuMY5N68h1teJ\nK8fGOcfKlSuZNGkSGVu28FSfPpz8t795H+WPGwdnn+13REkAoVCI386cySeTJjElECC3e3eyHn4Y\n+vf3O1pS08mGUm9bt1J2001UrFzJ3Wb0mTqVm8aOJV0zKUl1zrFh2DDyXnyR17p2pe+0aZw/eLD+\n0DpGDTZrh5lNBe4CLnLO/a3a8iDwNbDcOTe4IdZXUY5N8d69LJ05kw1PP03vPXv4bvPmtAwEsNGj\nYcwYOOEEvyNKAioqKmLGI49QOH0694TDWPfu5D32GJn9+0NAp1TUlxppqasDu3fzzzFjOG7BAh5P\nT8dNmMDNd9xBq1at/I4mCaxq5UoKbr+ddmvWsDIvjwMjRvDtiRNp3aaN39GSUkPO2nEOEAHWVF/o\nnKsEPgLObeD1k97y5csb7sWdw23ZwpdPPMFHgwaxtl07qlq35oL772f4Kadw5UMPkbN8ObZjB0ye\nfMxNdINuSyPTthwqLy+PSVOm8MDOnSx+9FHmlpezfcAASrKz2X7++ZQ++qg3HKgBNaXfiSS3ht4X\nC7dvZ9XEiaw45RRK27dn/YoVLLjvPiYUFvKL3/wmbk10U3pPaVsOlX7hhXRfvZoWO3bQ5Yor+I/f\n/Y4d7dox79RTeePOO9n66afHHrQOmtLvpb7q2kh3AAqdc6FaHvsCaGtmR/rc6VjXT3rHvJOVlnpj\nV6dNI/Lzn1N21VXsPussvuzQgbJgkJ09e/Lhbbexeds29l97LcGNGzl+3z5OWraMwJgx0Lt33I4o\nNqU3jLaldtnZ2QwfNYrbNm0is6CAV+68k9f37WPRnXeyu2dPdrdsyfp+/dh4113sXL2aSBwvYduU\nfieS3OK1L0YiEXZ98QV/e+klXhk/nucuuohlOTlknnwyzZ55BnfmmYRWr+ayr77ip7ffTlZWfM/B\nb0rvKW1L7dLat6f3rFl0KC6m8/z59OnWjTOfeYbWp57Kouxsnjr7bF644Qb++sQTFKxYwYE4n1je\nlH4v9VXX5jULqDzMYxXVnlPSQOsnDOdcnW6RSOSQ+xUVFRQVFX1zv7KykoqKim++VpSXEyopoaqo\niKqiIsq/+opAQQFtNmzgpIIC2u3dy4asLN6PRNiyfz8VublkdelCmwED6JCfzzn/9V9c0q2bxkdJ\n3HXq1IkRkyfD5MmEQiE++vBDChYuJLxsGe2efpo2Dz7I9kiEj1q2pLJ1a1xuLpaXR3qbNmQedxzZ\nubkEg0GCwSCZwSDBZs1Ib9WKjLZtyWjblkDr1lh2Nmnp6aSlpVFZWcm+fftIS0sjLS2NQCBAIBDQ\nvi0xqbVGRyLeraoKFw57t4PfRyK4QABnRnlpKXu+/BIXiVBVVkZlcTGh0lIOlJZSVVJCeO9ewsXF\nREpKCO3ZQ+jrr4kUFUFxMYGSEigpIVhaSruKCk4EemRm0j4vj6pOnci84w6yRo/mTM3dL/EUCNBi\n0CD6DBoEgNuxg75/+AM933yTyJIlBObNo3l5OZFwmHIzCps1Y0/LlhS3aUNF+/aEO3YkvX17MnJy\nyGjZkoycHDJzcgjm5hJs2ZJgTg6BZs1Iy8oikJ5OWiBAWiDAgYoKyktKvLqdmUkgLQ1LkeGAdW2k\n9wPHHeaxZtWe01DrN4rtGRm0raqKy2vV3H3SgOD//M8397PgkJkyMpwjHAhQnpZGeUYGoWCQ0lat\n2HnKKfzfsGFU9e3LiV27MrBTJ67t0IEMXYlOfJCRkcG5+fmcm58P997rLXSOlh99RP6f/kRZQQGh\n3buJ7NkDmzYRWLuW8IEDhCMRIuEwleEwFeEwmVVVtKiqIjMcpmUkQgZw8OOqAGAPPEAEbzxYPH3Y\nuTMXFBTE+VXFL3/t0oUzP/vsqM+z6C0Nb/9Kg2/2r3D0q+Nf+1sgessAWj7yCAAVZqSbcSAQID0Q\nIJSWRkV6OlUZGRzIzCSclQXZ2QRatcI6diS9TRtanHgirbt1o81ZZxHs0YNgs2a0jetPQOTIrEMH\nOkycCBMnHrK8KhTi802b+PL996n89FPYupUWn39O89WrySgtJS0UIj0UIr2qioxw2LtFImQ4R4Zz\n3zRu4L1v0oGMaI9z8P128P3lqt2o5euRFKan0ylU22CGxFHXkw0XAwOArJrDM8xsJdDdOde+IdY3\nM521IiJJKxVPNvQ7g4hIrOpbs+t6RHoNMBDIB/56cGF01o0zgOUNtX6q/SckIpLMVLNFJJXUdQDL\ni9Gv42ssHwU0B+YeXGBmXc2sZ6zri4iIiIgkgzoN7QAwsxnAWOANYCHQG7gFWOGcG1DteZ8BHZ1z\nabGsLyIiIiKSDOrTSBveEeVRQGegEJgHTHLO7a/2vAK8Rjo9lvVFRERERJJBnecmcZ5HnXO9nHPN\nnXMdnXN31GyCnXNdajbR9Vn/cMzsODP7vZn93cy+NrNyM9tsZs+YWbe6bkciMLMOZnaXmS03sx1m\nVmZmn5jZg2bW2u989WVmo81sjpltMLMqMwv7nelwzDMhmrXczLab2TQzi+/ErY0gug+9ZGZbzCxi\nZg17pZQGYmbdzezXZrbKzL4ysxIz+9DM7k6234uZ9Yi+F9ab2V4z2xfd1x42s+P9zteYVLMTVzLV\nbGg6dVs1O/HEo2bX+Yi038ysBzALWAVsA8qB7sBIIAj8h3Nuo38J687MRgOPAX8CVgKleCdiXgfs\nAPKdc1/5l7B+op9CtAY+BLoCJ9Yc2pMozGw63pCiV4HFQC/gVuBd59x3/cxWX2YWAb4GPsC7emix\nc66rv6nqz8zuB24C5gOr8WbC+w7wY+DvwHnRq6AmPDO7GLgbbzs+B6qA04ERQDFwhnOu0L+EjUc1\nO3ElU82GplO3VbMTT1xqdl0vMJKoN/51+fGZfmepR+ZeQLtalo+MbsuDfmes5/Z0qvb9AiDsd6bD\n5OyNN63lSzWW3xz9uQ/xO2M9t6dzte/XAVv9zhTjdpwFtKxl+ZTo7+smvzPGYRt/GN3Hbvc7i983\n1Wz/b8lSs6P5mkzdVs1Onlt9anZTuOzM9ujXPF9T1INzboOr/ejFwdlN+jRmnmPlnNt+9GclhGui\nXx+rsfxpvAsCDW3cOMfGOfeZ3xniwTn3gXOutJaHXsSb1z+p3g+HkXR1qgEl3c9CNdtXTaZuq2Yn\nlTrXqbrOI50wzCwdyMW76FR34F68C+T8ycdY8dIx+nWXrymaroNHwtZUX+icqzSzj4BzfUklh3Pw\n/fClryliYN4c+S3wrtx6GvAAXp1a6GcuP6hmyzFS3U4eKVmzk66RBi7B+yjqoF3Az51zf/QpTzxN\nxvvFPed3kCaqA1DoalxdM+oLoJ+ZpTvn4nOdeImZmQWAX+GNvUvG9/b1wOPV7hcAQ51zfz3M85sy\n1Ww5FqrbSSCVa3ajN9JmlgtMoG6XWQeY7pzbW+3+KuC7eBdy6Y03uD3PzNKcc4165nEctqX6a/0c\nb0zOk865v8QpYp3Fc1sSWBZwuBMgKqo9p6Rx4sgRTMc7mesu59xmv8PE4HVgA94RjjOBy4HjfE0U\nI9Vs1WyfqW4nh9St2T4M4D4Z72OacB1vXY/yeicAu4HfJeu24P0lFAb+F0hr7O2I87Yk7IkrwMfA\nzsM89mJ0u9L9zhnjtiXtiSu1bMuU6L74hN9Z4rhNp+P9pz/R7ywxZFfN/vfXUc1uvO1sknVbNTux\nb/Wp2Ukz/d2RmNkLwA+AbFf7xz8Jy8xG4J00sRi4Mtny12RmC4BLXQJOpWRmi4EBQFbNn7OZrQS6\nO+fa+xLuGJnZOrz9P+mmUqrOzO7F+3hwlnPuBp/jxJWZrQI6OOdO9juL31SzE0ci12xounVbNTvx\n1bVmN4VZO8D7yDANyPE7SH2Y2XXAU8BbwH8ne0FOAmvw9vn86gujJxmcQY2TWaRxmdkkvII8u6kV\n5KjmeHP3imq21J3qdoJSzfYkTSNtZu0Os7w33l+rW5xzXzduqtiZ2XC8oxrL8I5qHPA3UUo4OFXV\n+BrLR+G9YeY2bhw5yMx+BUwCnnPOjfA7T6zMrNYjY2b2HbwpoVY1biL/qGZLnKhuJyDV7GrPTZah\nHWb2KDAQb8qkz/jXXIXD8E6avMw5t9S3gPVgZpcDr+FdNWci3hW/qitzzv1voweLkZkNBvpG7w4F\neuD9lQqw1zn3W1+C1cLMZgBjgTfwprXpjXfFrBXOuQF+ZqsvMxuKN07S8C5OkAE8En14m3Nujl/Z\n6sPMxuKdLb0Nb7+J1HjKl0n03n4NbwzwO3jb0ww4GxgClAH/6Zxb51/CxqOanbiSqWZD06nbqtmJ\nJx41O5ka6YuBG/E2sB3ex4JfAMuBh51zG/xLVz/VPg45nG3JNG7KzH4P/OwwDyfUtpiZ4R3ZGAV0\nBgqBecAk59x+H6PVm5n9Geh/mIf/4py7uDHzxOoo+w8k17b8ELgW+BbeGd8Orzi/BUxzzn3uY7xG\npZqduJKpZkPTqduq2YknHjU7aRppEREREZFEkjRjpEVEREREEokaaRERERGRGKiRFhERERGJgRpp\nEREREZEYqJEWEREREYmBGmkRERERkRiokRYRERERiYEaaRERERGRGKiRFhERERGJgRppEREREZEY\nqJEWERGRhGJmV5nZL81skZnlRZedY2aFZtbb73wiB6mRFhERkYRhZicCHZ1zU4DOwLejD+0AHNDL\np2gi/0aNtKQUM/u+mf3CzJaYWetqy68xs9f9zCYiIgB8H3jOzPoA3YDVAM65HcAsYJuP2UQOoUZa\nUka0cT7NOTcV6AT0r/bwj4D9vgQTEZFvOOdmOef2AjcAbznndlV7eB/wgT/JRP6dGmlJJd8DXjCz\nbwGnAP9X7bELgXd9SSUiIrX5EfBajWUB51zEjzAitVEjLSnDOTcv+tHgdcA70e+JNtZ5qJEWEUkI\nZnYc0J5qR5/N7BJguZk1N7MbzWyemaVHH5tlZj19iispTI20pKIfAq9Uu98f+No5t8GnPCIicqi9\neMPt0gDMrBVwgXPuL8Bg4FmgL5AZff730Nhp8UG63wFEGlN0GqUTOXRYR39gpT+JRESkJudcyMyG\nAf/PzD4CMoD7og8vxJu54x/Ouf1m1hXY4Zyr8CmupDA10pJqDkRvDiD6UeD3gV/5GUpERA7lnHsd\n+LfZlJxz+8xsELAguqgfsKoxs4kcpEZaUkq0AI/iX0c5TgGygRX+JhMRkXo4HlgT/X4AsNjHLJLC\nzDnndwYR35jZvcCNwPFObwYRkaRgZucBQ/BORpyJN7XpP/1NJalIR6QlpZjZFGCVc26hmRleIZ6p\nJlpEJHk451YDq82sM7BLTbT4RbN2SMows7bAnUDb6KLb8S45+4BvoUREpF7M7HwzWxq9ezvwaz/z\nSGrT0A5JKWZ2KxDEm590PzDFORfyN5WIiNRVdJaOHwN78PqYJ32OJClMjbSIiIiISAw0tENERERE\nJAZqpEVEREREYqBGWkREREQkBmqkRURERERioEZaRERERCQGaqRFRERERGKgRlpEREREJAZqpEVE\nREREYvD/Ac44/Jk+vsmuAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "boltzmann_vel2 = 0.5*sim.integ.beta*sim.mass[1]\n", "boltzmann_pos2 = 0.5*sim.integ.beta*sim.mass[1]*sim.pes.omega[1]**2\n", "f, (ax2, av2) = plt.subplots(1,2, sharey=True)\n", "px2 = ax2.plot(lx2, np.sqrt(boltzmann_pos2/np.pi)*np.exp(-boltzmann_pos2*lx2**2), 'k-', plotbinsx2, x2hist, 'r-')\n", "px2 = ax2.set_xlabel('$y$')\n", "pv2 = av2.plot(lv2, np.sqrt(boltzmann_vel2/np.pi)*np.exp(-boltzmann_vel2*lv2**2), 'k-', plotbinsv2, v2hist, 'r-')\n", "pv2 = av2.set_xlabel('$v_y$')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we plot the 2D histograms for each degree of freedom. These should be reasonably circular." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj0AAAEdCAYAAADn1daqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xm0JXV16PHvpiFMytDiyNQiiCKKAxI1PmyNinHWEHUZ\nnoKovIjj0/gEspjUJI+YODxNHJ+oqA+c5yGKjYLgEFDBYCLibaQJGCAICCjQ+/1Rdem61fdUndvn\n9K1bt76ftWqdqvOrX9XvDHfXvnXq7BOZiSRJ0nK3RdcDkCRJWgwmPZIkaRBMeiRJ0iCY9EiSpEEw\n6ZEkSYNg0iNJkgbBpEeSJA2CSY8kLXERsU9EnBwR50bEryPi+oi4ICKOjYjtuh6f1BdhcUJJWtoi\n4m+AlwGfB84DbgUeCzwX+DHwiMz8XXcjlPphkElPRAzvQXcoM6PrMUh9FhEPBX6emTfU7n8jcCzw\nisz8x3n6Ges0WPMde7bsYiBLwwkLXH8NsHqC/XXZv8t9nzTBfiUBZOb5I5pOB44D9h/de6GxbtYa\nJosb7sf9dLmf+Y89XtMjSf21e3l7VaejkHrCpEeSeigitgCOp7i+52MdD0fqhQF/vLVQq3rcv8t9\nS9pM3g4cBByTmT+f/uZXTX+T7sf9dLwfz/SMbVWP+3e5b0nTVl7AfDTwnsw8ZfPsZdXm2az7cT8d\n7sczPZLUIxFxIsXFyx/IzJe191hTmV+F/8RoeZopp2a9Tnoi4r4Un2k/FLgXsBVwGfBl4O8y88oO\nhydJUxURJ1DEvFMz8yXj9Vq9GUckLRWrmJvQnzXvWr1OeoDdgHsAnwYuB24DHgi8FHhuRDw4M6/u\ncHySNBURcTzF988/lJkv6no8Uh/1OunJzDOBM+v3R8R3gDOAw4G3LPKwJGmqIuJo4ERgLXBmRPx5\nbZWrMvMbiz4wqWd6nfQ0uKy83bnTUUjSdBwIJLAHcOo87WcBJj1Si2WR9ETE1sCdgG2ABwB/SxEg\nvtzluCRpGjLzCOCIrsch9d1y+cr6i4H/BH4FfBXYETgsM8/pdFSSJGnJWBZneoDPABdTnO15CPB0\n4K6djkiSJC0pyyLpycwrgCvKxc9HxKeBH0TEtpn5vzscmiRJWiKWRdJTl5kXRsQFwMuAEUnPmsr8\nKizYNS0zjFMgSpKkxbYsk57StsDK0c2rF2scA7OKcQpESZK02Hp9IXNE3H3E/Y8F9gfOXdwRSZI0\nylYtkza3vp/p+aeIuCdFgcK1FF9ZfxjwPOA3wGs7HJskSVpC+p70fAx4IXAYxbe1kiL5+SfgLZl5\neYdjkyRJS0ivk57M/CTwya7HIUmSlr5eX9MjSZI0LpMeSZI0CCY9kiRpEEx6JEnSIJj0SJKkQej1\nt7ckSZqupiKBbYfM2zbjvm+dcNsCz/RIkqSBMOmRJEmDYNIjSZIGwaRHkiQNgkmPJEkaBJMeSZI0\nCCY9kiRpEEx6JEnSIFicsDNNRagmYQErSX3XFh+bDl1tfbfdjNu+uaW9LT43je2GCbfdVDhxOMcN\nz/RIkqRBMOmRJEmDYNIjSZIGoddJT0TsExEnR8S5EfHriLg+Ii6IiGMjYruuxydJkpaOXic9wIuA\nVwGXACcBrwN+BrwJOCcitu5wbJIkaQnp+7e3PgH8dWZWL2t/b0RcAhwLHAn8YycjkyRJS0qvz/Rk\n5vm1hGfW6UAA+y/ykCRJ0hLV66Snwe7l7VWdjkKSJC0Zff94ayMRsQVwPEW1pY8t7t4XUnCw6alv\nK55VLzK1qYWlhlOQStJSMmkBwTtP0PfuLe1N2rY9afHC6xvadpigb9u+28bVVNgQ+nQsWXZJD/B2\n4CDgmMz8edeDkSRJS8Oy+ngrIt4IHA28JzNP6Xo8kiRp6Vg2Z3oi4kTgOOADmfmy9h5rKvOrykmT\nmyknSZKWlmWR9ETECRTX8ZyamS8Zr9fqzTiiIVvF3ATyrG6GIUlSTe8/3oqI44ETgA9l5ou6Ho8k\nSVqaen2mJyKOBk4E1gJnRsSf11a5KjO/segDk6Qpi4hjgIcADwPuDcxk5l7djkrql14nPcCBQAJ7\nAKfO034WYNIjaTl4M3ANcD6wU8djkXqp10lPZh4BHNH1OCRpEeyVmTMAEXEhsH23w5H6p9dJT/fq\nhajqT+dWDW31IlfVwlNN25lPtShV/Vc56kWjmopQ9afAlDQ0swlPfzTFrbYifytb2puKE+7W0vdu\nLe1NxQu3a+k7qaYCg5e39F3X0t70AwXz/ZpTVVvhwyZL67jS+wuZJUmSxmHSI0mSBsGkR5IkDYJJ\njyRJGgQvZJakZW1NZX4V/uSOlqcZxvkJJJMeSVrWVnc9AGkRrGKcn0Dy4y1JkjQInumRpB6IiMOA\nPYEA7gpsFRHHlc1rM/O0zgYn9YRJz4JtasHBXWttezesW99OvSBWvVBUdblegKpesKraXi9U2FbY\nUFKHjgQOrt13cnl7FrDISU9b0dSmAoI7NLQBPLClvakAYVNxwZausHFoHrcNYJuW9lta2i9qeF5m\n9mvue+U+LRv/SUNbW2HDNv0pXmjSI0k9kJmP7XoMUt95TY8kSRoEkx5JkjQIJj2SJGkQTHokSdIg\nmPRIkqRBMOmRJEmD0OukJyKOiYgzIuIXEbE+Ii7tekySJGlp6nudnjcD1wDnAzttnl3UC3BtO2Ie\nNi7G9aAR88CDa6vuX5mvF8CqF7y6rlbA6pLK8nm1yluX1wtWnV+Zn6HZzSPmJS1/bcUH6/FvIf3r\nxVrr2ioINhQgrMfWusNa2lePbtrtYT9v7HoXrmls//G6lsFd0lDd8IfNXflsy+t19sMaGle2bLzt\nvXB5Q9uvJ9z2dIsX9j3p2SszZwAi4kJg+26HI0mSlqpef7w1m/BIkiS16XXSI0mSNC6THkmSNAgm\nPZIkaRBMeiRJ0iD0/dtbE1hTmV9VTprcDO1fhZckafENOOlZ3fUAlqlVzE0gz+pmGJIk1Qw46RlX\nvQBXdflutbaH1pbvvWH2mbWm2vI+L/zxHfMH8i+NI/oVu89Z/tfb97tj/tq31Yp+fbZWyPDs6hjb\nij7dVpm3OKG0/DQVhms7PLS1NxW8W9Xc9U4NxQeh+X/WtzR3PWDf8xrbH883Rrbdh0sa+27N7xvb\nb991RWP7Bbs+ZGTbxx7x/Ma+v9n7Ho3t7NLQ9tl7NzQC3NDSfn1DW9uxo63d4oR3iIjDgD2BAO4K\nbBURx5XNazPztM4GJ0mSlpReJz3AkcDBtftOLm/PAkx6JEkS0POkJzMf2/UYJElSP/iVdUmSNAgm\nPZIkaRBMeiRJ0iCY9EiSpEEw6ZEkSYPQ629vbR71Yl31p6i6vKrWVivwdGhl/sVzm/7kkE/PWX4Z\n77pj/qk/P3Puyr+s7WafuYv/fO9H3zF/8muPn9N29p2eMHflGyvFCn/0oNqG6wWmqsv1Io0WK5T6\nr+kQUP+bX2j7rqObttlndBvAo5ubmwoQPnHfzzV2/eOG4oMAr//1O0c3fr6xK2zT0r57c3M+cHSV\nld1X/qqx7zue8YrG9iuv3mt04w8bu8Lle7essDmLEzYV0Fx44ULP9EiSpEEw6ZEkSYNg0iNJkgbB\npEeSJA2CSY8kSRoEkx5JkjQIJj2SJGkQTHokSdIgWJxwI/WnpF4YqVLcj9vmNt2jtmqlwNZBh5w1\np+kkTpiz/PBTLrpj/ppj527mwtvnLtfLRD3hg2ffMX/F4e+f03b24x8/d+VLYsP8j+5W21L9sVaX\nLUYo9VNTcbfbGtratPVdNbppp5auuzU377nvz0a2HcLXGvs+l9ObN/6U0U1XtRTx23JFc3ubu7x6\ndNu+b/m3xr5P5wuN7e99xKtGN96vsStc3vZaN73H2rRte+EFCJt4pkeSJA2CSY8kSRoEkx5JkjQI\nvU96ovCaiLg4Im6OiMsi4i0RsV3XY5OkaTHWSZPrfdIDvA34e+Ai4OXAGcAraf89XEnqE2OdNKFe\nf3srIvaj+OP/ZGY+p3L/DPCOiHheZv6/rsYnSdNgrJOmY6IzPRHxwoh4Sjm/XUTcfTrDGtvzy9u3\n1e5/H3ATcNjiDkfSchIR+0bENl2PA2OdNBWTnum5D/D7cv5m4LERcWtmfn3C7Y7rQGA98IPqnZn5\nu4j4EfDwRRqHpOXp9cATI+IK4Fzgu8C5mfmrRR6HsU6agkmTnp8CW0TEysy8FvhSRLxwCuMa172A\nqzNzvupF64BHRsSWmTlBBa76prcdvWr92axUEdyduTHywGsvmrNcLUi4rlaMcA3Nyyf+esP8o/ju\nnLY97zO3oNXaOVWogrnqj23LhjaLFWr5y8wjASJiL+CRwJOA90fElcDJmXnaIg1lgli3qWG+7W+8\nIRa2qRdyrWspXrgdN41suw+/aOy760XXNrZ/qqEA4YWNPWHl7c3tbc/Yn797dNuzX/WVxr6f2f1Z\nzRvf6ZbRbXdqO5m5Q0t703ul7VG3vT+nW5xw0qTnUcCOwLEREcCPgf8CPjTpwMa0HfC7EW23VNa5\nfnGGI2k5iYgtgT/IzEuBS4GPRsQa4NvA6yLiTpnZcKiaGmOdNAWTJj3nZebHASJiF+DJjP7D3Bxu\nAu46om2byjqStClOAw6IiEuALwAzwEMz88PAyyPiqEUah7FOmoJJk55fR8RTMvNLmXl1RKwHHght\nP24yNVcA94+IreY57bsrxengER9tranMr6Lxd2K0ADPlJC0LZwEvBO5P8atMDwQ+ChAR3wC+uUjj\nmCDW/XNlfi+KSzGl5WaGcY49Yyc95Vcm7w+sBc7PzPWZ+c2I2CUitsnMW4AbgF9u0ng3zQ+AJwAH\nAedUxro18GA2vvylYvVmHdhwrWJuAnnW/KtJS8x8MQ44FXg68OXMfHOty/HAlYs0vAli3RM278ik\nJWEV4xx7FnKm59tAAh8A/iQitgDWZOYdW87Mzy10mBM6HTgWeDWVQAC8lOLqqY8u8ngk9deoGPeJ\n+VbOzO/Od/9mYqyTpmAhSc+TgCsy84rZOyLikIg4Gfjr8kzPosrMiyLiXcDREfEp4MvAfsArKILV\nxxd7TJJ6a8nFuFnGOmk6xk56MnOjL/Jl5tci4nvAG4ATpziuhXgVxUdqL6W4kPpq4O3ACR2NR1IP\nLeEYN8tYJ01oIdf07Aw8huIbW3d8jp2Z10XEdL9IvwCZmcBby2kKFlLS585zF+v/B169YfaXtQul\nr1i5cs7yrnffUDviC1fMaeL+tc0+bfvaHTtumP0Vu89pWkGtcMSdaFB/7Lc1tEnLy1KNcZVxbGKs\na/rbbaqh0lanZ6uW9pWjm37W0vXBzc0Xr9tvZNu/7jq6DeCQezfXzn10Q1tzhZ92bRVrZn47um3l\n7juObgRWtMXo6zZnYfGmR9b2p9P2rDQ9roX/WS7k462PUVwldJ+I+DrwaeB8YAWtb1FJWvKMcdIy\nt5Df3vpuZt4feCjwC+DNFAHhLIpgIUl9ZoyTlrmFJD0/jYhXAZdn5qsy857A3YGdM/Mzm2d4krRo\njHHSMjd20pOZnwY+Cfxp5b7/HPFbMJLUK8Y4aflbUEXmzFxHUcNCkpYdY5y0vC3k4y1JkqTeMumR\nJEmDYNIjSZIGYdJfWR+ganmqG+Y2XVdb9YsbZs9/9B/NafrwfV4wZ/mYH7ztjvnDT6ttp1askHvP\nXfzt4Rty1w9yxJy2Sz/1gLkrr6ku1K/PrL8dbh0xL6k/mv52mw4BbUXj2ooXfn900y2rm7vWY2nd\nV0cX2vvgkYc3dt17+0sa2//snV8c2faSyxq7tv++cr3S7ALa/4HnN3b9Dgc3b7vpYc80d4WbWtqv\nb2i7oaFtnPbpHns80yNJkgbBpEeSJA2CSY8kSRoEkx5JkjQIJj2SJGkQTHokSdIgmPRIkqRBMOmR\nJEmD0OvihBFxFPDfgIcB+wCRmSumu5fbasvVglyX11ZdOXf5h/tsmD815jQde+hb5yyvOWD1HfOr\nX79mTtvetapSv2L3Ocvf4PF3zH/lnGfPHcOpcxfnFic8p9Y4g6Qhqce3qrbihG2aihdmc9dLorl9\nzeimn+92QGPXdx9yVGP7FUffa2TbM/lMY982O9/eXHXxrSteM7Lt3TSP+8qv7dW88882tF3U3LX9\n2ND0PmorLtiWhky3OGGvkx7gDcBK4AJge2DXbocjSZKWqr5/vPWYzNwxM1cDP+56MJIkaenqddKT\nmW2/hCJJkgT0POmRJEkal0mPJEkahM4vZI6IHYHX0Ho5/x3enpnNl8BLkiTVdJ70ADsBxzN+0vMR\nwKRHkiQtSOdJT2aupZOP2dZU5leVkyY3g/V+JElLUedJT3dWj7h/IYWQrq0tzy0iyEyluNc7d5vb\n9rO5i1+/3zM2zD/4GXMbt6nt5uracrWw1Nm1tvNqy/xLZf6qWlu9mNjNDW2jrGJuAnnWmP0kLS1t\nf/NNBekArmloO7+560UPa25vOnI1FeEDztzmqY3tP3r0Q0a2/XDFgY19b24p6PivK/ZrbL947QNH\nN351q8a+fKO5mS82tN32y5bO61raL29om/R9NF1eyCxJkgah12d6IuKpwGzN8b3L+44rl6/LzHd1\nMjBJkrTk9DrpAf4UeEHtvpPL27WASY8kSQJ6/vFWZh6RmStGTC2/viZJkoak10mPJEnSuEx6JEnS\nIJj0SJKkQTDpkSRJg9D3b28tgnqxwupTVi+6VC/2V+l7Xa3tk3evrVspXninWlN9+cam5XqRqHpR\nqZnKfL24Yr1I1OIWjZK02BZSjHWh6vGl6sKWvs1F/vjRPqPbLm8p4vfD5uZrH7HryLbTdnlJc+ed\nmpvrRWkXZKalvV6Ytu6W+vGp6pKGNmgvTtj0WrcdRzbne3BjnumRJEmDYNIjSZIGwaRHknogIo6K\niNMi4uKIuC0ibu96TFLfeE2PJPXDG4CVwAXA9sDoi08kzcszPZLUD4/JzB0zczXw464HI/WRSY8k\n9UBmXtb1GKS+M+mRJEmDYNIjSZIGwQuZF6xaaKlenLBehKladOn6Wlu9mNPMhtkbV85tqhcj3Gi/\nlb4bjaG+n+qYFjJ+SZOKiB2B1wA5Zpe3Z+Z1m3FIDSb9+6/Hl6qWAoKcv+nbvrqhcCHA1S2b/mFD\nYcRtWsa9Tcu2W1/JprdFUwFAaC8g2NQ+09K36bXsF5MeSVo8OwHHM37S8xHGOFRKGo9JjyQtksxc\ny6JfVrCmMr+qnKTlZob2M1YmPZK0zK3uegDSIljF3IT+rHnX8kJmSZI0CL090xMR9wJeCBwC3BfY\ngeLc1peBv83Mtqu+JKk3IuKpwAHl4t7lfceVy9dl5rs6GZjUI71NeoCnUVwQ+CXgs8ANwEHAq4Hn\nRMRBmfnrDscnSdP0p8ALavedXN6uBUx6pBZ9Tnq+DexZS2w+EBHfB94HvA54fScjk6Qpy8wjgCO6\nHofUZ729piczLx5xJuf08nb/xRyPJEla2vp8pmeU3cvbKzfP5m8dMQ8bF/erFnSqP9X1S46ail61\nvUxNhaPaChBWWYxQ0rja4sUk8WSSYnhtVzXcuaW9IRbf0jKuW9pidVtRxnoR26qm2A1wTUv7JFd7\ntL0e/Tl29PZMT4OTKAp/fajrgUiSpKWj8zM90yzLHhGvBQ4F3p2Z839JX5IkDVLnSQ9TKsseES8G\nTgG+ALxiaqOTJEnLQudJzzTKskfEi4D3AF8FDs3M29t7ranMr8LS7NMywzilwCVJWmydJz2Tiogj\ngPcCXweelZljXlG1evMNatBWMU4pcEmSFluvL2SOiMMpavJ8E3hmZv6+2xFJkqSlqrdneiLi6cD7\ngd8AnwAOjYjqKjdm5ue6GJskSVp6epv0AA8BguJC6PfM074WMOmRJElAj5OezDyJoibPEtJ0OVG9\nrV7sqa1o1abqT9EoSUPRVuyurRDfuoa2tiJ827a0b87DYlucb4rXkz5nTduepG+/9PqaHkmSpHGZ\n9EiSpEEw6ZEkSYNg0iNJkgbBpEeSJA2CSY8kSRoEkx5JkjQIJj2SJGkQeluccPlZPsWfJGkybfFw\nkkJ7bUX+luphcTgFBDcnz/RIkqRBMOmRJEmDYNIjSZIGwaRHkiQNgkmPJEkaBJMeSZI0CCY9kiRp\nEEx6JEnSICzVKkytIuKuwCnAQ4HdgO2Ay4GzgL/JzF90ODxJUie6LNJngcClrrdJD7AzsDfwNWAt\nRZnNfYAjgT+LiD/MzJ91OD5JkrSE9Dbpycx/B/5b/f6I+BTwfeDl5SRJkrQsr+m5rLzdudNRSJKk\nJaW3Z3pmRcSWwI7AVhQfb50IJPClDoclSZKWmN4nPcAhwBcqy1cCr83Mj3U0HkmStAR1nvRExI7A\nayjOzozj7Zl5XWX5XODxwLbAfsBzgZ0jYkVm3j7VwUqSpN7qPOkBdgKOZ/yk5yPAHUlPZl4LnFku\nfikiTgN+AtwN+IspjlOSJPVY50lPZq5lihdUZ+Z/RMQ3gCMj4pWZOaJwwprK/Kpy0uRmykmSpKWl\n86RnM9kWWAHsAFwz/yqrF280g7KKuQnkWd0MQ5I2iQUGl7PefmU9Iu424v79gD8GfpGZIxIeSeqP\niLhXRBwTEWsi4oqIuDEiLoqIUyJiZdfjk/qiz2d6jomIJ1B8NX0GCGB/4L9TPK6XdTc0SZqqp1Fc\n+/gl4LPADcBBwKuB50TEQZn56w7HJ/VCn5OeL1D85tafUVy0vAJYB5wO/H1mXtzh2CRpmr4N7FlL\nbD4QEd8H3ge8Dnh9JyOTeqS3SU9mnsmGb21J0rLV8E/c6RRJz/6LOBypt3p7TY8kid3L2ys7HYXU\nEyY9ktRfJ1HUOPtQ1wOR+qC3H29JUt9MoQJ9dVuvBQ4F3p2Z1oaQxmDSI0mLZ6IK9LMi4sXAKRRf\n6HjF1EYnLXMmPZK0SKZRgT4iXgS8B/gqcGj7bwyuqcyvwurzWp5mGOfXAEx6JKknIuII4L3A14Fn\njf6ZnarVm3dQ0pKwinF+DcALmSWpByLicIqvp38TeGZm/r7bEUn945keSVriIuLpwPuB3wCfAA6N\niOoqN2bm57oYm9QnJj2StPQ9hOKndnaiuJ6nbi1g0iO18OOtsc30uH+X+5Y0qcw8KTNXNEx7TX+v\nM9PfpPtxPx3vx6RnbDM97t/lviX104z7cT/Lbj8mPZIkaRBMeiRJ0iBE5riFQZePiBjeg+5QZkb7\nWpKmzVinIZvv2DPIpEeSJA2PH29JkqRBMOmRJEmDYNIjSZIGwaRnHhFxVEScFhEXR8RtEdHyK8bz\nbuMFEXFJRNweEesj4paI+GxErCrbIyJeU+7j5oi4LCLeEhHble0zZb/6dHtErCz7XlW5//qIeNts\n/5axrY+IrE3rI+L6TXiM50fETRFxZUS8LyJ2WehzJalbU4p5axpi1kOntZ9yO0+OiHMi4saIuCYi\nzpiNrZV1JopP88Tgarx85zixttxOY6xf6HoN+5nvud8orkfEMeXz9Yuy/dJxn5Padhqf32nsp+04\nuEnj9kLmjUXEL4GVwAXAXsCumbliAf1fA/w9Rdn4tcBFwBMoksyrgAOBY4BXAJ8CvgrcH3gl8O3M\nfHw5hpuAN5XbqXok8PJy/ofANeX2AzgzMx/fMr4EbgUuLh/fTcBrgFsz8xMLfIzfAj4O7Aa8lqLC\n1EGZefM425HUvUljXrmNbwH7Aa9m45j15cy8bkr7eTbF749dQPF7ZDtSxK/bgAMz88ppxKdKDL4M\nOAT4PvAT4F7AEylj9RjbeTsNsX6h6zXsZz3wbeC9taY5cb1c7xrgfIpj0W8WWtF7nOd3SvtpOg5+\nIjNvXcj2AMhMp9oE7FGZ/wJw+wL63gW4EfgdcCmwbXn/U4H1wO3A6eXtGbW+Ly/XeR7wS4oEpr79\n/cq+66v9K33XA89rGeN64P+W8xcCly7w+Zl9jOdSJs61x/iGrl9DJyen8adJYl6l37faYsmk+6H4\nvch11dha3n8ARdLz7mnFpzIGn9cWq1u2sd84/cddr2Vfd8T1lvVWVeY3W/yfdD+V12Cj4+Akkx9v\nzSMzL5ug+7OAbYGtgPdn+R9FZn6R4g/1ZuDp5bpvq/V9H0VWe9jsHRGxIiLuXFnn+RQZb9b6z/Zd\nX+3fJCK2YuPseRyzj/H/ZPnOhDmPcaz9S1oaJox5c5Qf09x5vrYp7OcxwD2pxNZyuz8G1gDPBZ7N\n9OLTPcrb+tmTjWL1CM8vb9ti/bjrtYqIrSJi+1HtmTkz7rZGGCv+T2E/d5jnOLjJTHqm70A2JCXn\n1drOA7YDtqFITn5QbczM3wE/Ah5e3vWHFG/430TEf0XEqcAflW1z+lf6rq/0b3Joue0HAHtExDsi\nYofxHiIHVh5P3XnA/cb9HFrSsrIrxVmA35TX23wqIvad4vYfzvyxlfK+HYA/rizPt85C4tNuFMfJ\nr8/G4Ii45zyxepQDGS/Wj7tem9m4fkMU13wuJK6Pa7Hj/0bHwYi456ZubMvpjUule1Xm19XaqsvX\n5/yfR66juGbnKxRZ/s8oXqfVwEso/uBvB66ep/+6ct1dImLLzLxtxBi/B5wB/AJ4F7AzxWnUgyPi\nUZl5U+Mj3PAY649v9r4o17mkZTuSlo9LgbMprnu5neJg9QrgcRHx6Mz86RT20RZ7APZoWWfc+HQR\nxT+o2wBHsSEGPy4iDiq39ciWWHsv5o/Vs2N5ZERsOe56DfuBuXF9B+DJLCyuj2sx4/9FFB+j1Y+D\nj4uIgzLzyoVucNkmPRExe3HbuFdqvz0zr5un/32LxThhzP7VDPd3tXVuqcyPevPOrvP8zKxedX9G\nRHwH+BjFY6pvu7797YB5v42VmY+cnY+IN1Nc1Pw+4M3Aq4C/GTG26rZn/wsZNQbP9EiLaNKYN4X9\nzH50tV95exPFRa4vBP6B4mLgSffzRxQH1f8REfWx/2d5O5X4lJlPi4hLgN9l5hlsiMEfBU6qbWvU\nN1+3Y/5YXR/LuOuN/IZtNa6XTouICxk/ro9r0eJ/Zj6tdlf9NThqodtctkkPsBNwPOMHgI8A1T+i\n2f6z17wcP2b/aja9dW2dbSrzo5772XU2ysoz8+PlR1x/MM+269tfaFb/d8AJwFNo/+O4CSAitp7n\njT9y/JIiX3+HAAAGgklEQVQ2q0lj3ubaz/eAx46IFwvdT5TTK+dZ/8Tydprx6SbgrrMLZQx+M0Wc\nPGeMbc3p3zCWcddbqIXE9XF1Gv9rr8GCLdtrejJzbWZukZkrxpwuna8/8EVg/QL6X1HZzK61YVWX\ndygvJK7bleI056gzQTeWt7vM039XijNITf3nVa5/BTBOHYvZx1h/fLP3JXOfB0mb2aQxb3Pth+Kj\niRUUH6NPtB/gWIprX54wz35m/xGcPeM0jfh0BRvH2hmKONkWq0f1r45ltv+46y3IAuP6uJZC/J9h\nEx/Tsk16OjR7IVpQXJtTNXtB1i0Uz/1B1caI2Bp4MLWL2WpmzxDN6V/pu0VL/3mV/XejqCPU5gfM\n//igeIz/NsXPjyX1230p/hm7dgrbaoo9j6T4+OcbDessND79gI1j9d4UcbItVo/qP1+sH3e9BVlg\nXB/XUoj/s6/Bgpn0TCgi7hER+0bEtuVdn6P4WvqtwItjQ4XlpwH3ofiq3xfKdV8dEbuX/VcALy3b\nP1vZ/l6z336IiKOBO1Nk0ltQFAGb9VKKz1G3oPi8c07fyvZGVbF8E8V/Y5+vrV8d36zZx/jyiIjK\nurOP8bQR+5DUc/PEPCJih4jY6HgSEU8BHgV8PTN/P+l+gLOA/6ASW8t1D6D4OvsZFPFz7Pg0X4yL\niNmzUqeXt68u7z+aIomYoYjVH6302Sje1vtXzMb6jy5wvanE9XEtVvxveQ3q686+Bpv0mKZW8Gc5\nTRRFlo4rp4spvokwu3x0bd1TKU63Hly573+yoYDgWuBLFBep3Upx2u+ewDvKdf6zXO+9wO+Bb1Jc\ndPYT4BSKipbrgc+Ut/9WrrueIvn5F+Brlf19szKOGWrFvyguKPx3ipoWn6P4z+jWylgPr62/pmzb\no3b/7GM8k+Jq+pOAGyiutt+u69fQyclp/GkKMe8ZFN8aehvF9TYvAz5EcYbnSmDvaeynvP/Qcrvn\nA38BvKHcxzrgnuU6Y8en+WJcLQZ/q2xfV95eOxura9vZKN6W98/G+k8BR1JUMp6v/7jrjYrr36W4\naPkoiurI3yzHew6wdWXdw8rn+6/K5+2aymtwWNtzM+7zO+l+aq/ByyjeV9Xj4F026b3e9R/bUpyA\nD5Yv6HzTpfOse9s8f5gvoAgCs8nILeULdu+yPSi+mfBbNnwG+ncUZ2seRfHfygwbkpuflm/oHSp9\nr2JDFebrgbdS+YOmqGZ5W21cTy/ffFmbZqtFn1lb/1vl49tjnufpBRSl4G8q39TvA3bp+vVzcnJa\n2DRpzAPuR3G24udlLLq5nH8HZSIyjf1U2p5McZC/sYxnp8/G1so6Y8Wn+WJcLQb/liL5mP3H9fLZ\nWF3bzkbxtrx/Nl5fXD4vvxrRf9z1RsX1r5R9bqJIQM4H/hfwB/M83lGvwdTi/6T7mec1uInKcXBT\n3+v+9pYkSRoEr+mRJEmDYNIjSZIGwaRHkiQNgkmPJEkaBJMeSZI0CCY9kiRpEEx6JEnSIJj0SJKk\nQTDpkSRJg2DSI0mSBsGkR5IkDYJJjyRJGgSTHkmSpiAinhQRfxURX4uIlZX7nx8Rn+lybCqY9EiS\nNKEyyXlAZr4J2AM4uNL8HOCmTgamObbsegBqFhEvBXYB9gU+AuwJ3A3YH3h9Zq7rcHiSpMITgY9H\nxIOAvYHvV9oeDRzXyag0R2Rm12PQCBHxEuAnmfm9iHg48M/ACyn+Y/gq8OTM/FqXY5QkbRARbwX2\ny8xDyuUHARcA+2fmxZ0OTp7pWeLukpnfK+f3BG7PzM9FxLbA6sz8TodjkyRt7FDg5MrywcA1JjxL\ng0nPEpaZf1tZPBj4dnn/zYAJjyQtIRGxM7Arcz/aOhg4u5sRqc4LmfvjscCargchSRrp9+WUABGx\nL/Akyn9Y1T2TniUqIraIiMdH4W7AA6gkPRHx+s4GJ0naSGb+Fngp8IYyRr8W2B7PzC8ZJj1L11HA\n14F92PB1x8sBIuIZwE+7G5okaT6Z+eHMfH5mngJcAVwNnN/xsFTy21tLVEQcAPwl8DPgJ8CdgccB\nM8AvM/O07kYnSaqLiDcC52bmlyMigIuBj2bmGzsemkomPZIkTSgidgHWAS/JzA9HxF8CfwIckpm3\ndjs6zTLpkSRpCiLilcDWwN0pLkl4ownP0mLSI0mSBsELmSVJ0iCY9EiSpEEw6ZEkSYNg0iNJkgbB\npEeSJA2CSY8kSRoEkx5JkjQIJj2SJGkQ/j8v2IOfSg9tCAAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "f, (ah1, ah2) = plt.subplots(1,2)\n", "ah1.set_xlabel('$x$')\n", "ah1.set_ylabel('$v_x$')\n", "ah2.set_xlabel('$y$')\n", "ah2.set_ylabel('$v_y$')\n", "hist1plt = ah1.imshow(hist1.T, extent=[xb1[0],xb1[-1],yb1[0],yb1[-1]], interpolation='nearest')\n", "hist2plt = ah2.imshow(hist2.T, extent=[xb2[0],xb2[-1],yb2[0],yb2[-1]], interpolation='nearest')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The two plots above should look reasonably similar to each other, although the axes will depend on your choice of $m$ and $\\omega$.\n", "\n", "The final plot is of the kinetic energy information:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwUAAAEjCAYAAABjD617AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8VPW9//HXNyirhCZAUTEsxaKCtXUBK6KdUK+Ie0RE\nELCKVm/VXLEVdwXFpVxta2zVnywVRBSspmoLKGAiUGvRihvlWhFZFFAQCmENkO/vjzMTh2G2Mzkz\nc2byfj4e80jmbPM5+/dzvud7jrHWIiIiIiIijVdBtgMQEREREZHsUlIgIiIiItLIKSkQEREREWnk\nlBSIiIiIiDRySgpERERERBo5JQUiIiIiIo2ckgIRERERkUbOd0mBMaaFMeZzY0ydMaYiyXGqg8NH\nfvYZY05Id8wiIiIiIrnsoGwHEMV9QDHg5q1qFtgA3AiYiH4rPIpLRERERCQv+SopCF7V/x/gZuA3\nLkffbq19zvuoRERERETym29uHzLGFAATgFlAZYrTMMaY1p4GJiIiIiKS53yTFAA3Ad2B61McvyOw\nDdhijNlmjHnRGHOUZ9GJiIiIiOQpX9w+ZIzpCowBxlhr1xhjOrucxApgEfAhsA84GbgB6GeM6Wut\nXeplvCIiIiIi+cRY66Y9b5qCMGYOcBhwgrV2XzAp+Bz4vbW2PMVp9gWqgfnW2v6eBSsiIiIikmey\nXlNgjBkGnAGcZq3d59V0rbWLjDELgFJjTDNr7W6vpi0iIiIikk+ymhQYY5oCj+A0Lv7aGNMt2OuI\n4N82wW4brbVbUviJlcBPgCJgfZTfz341iYiIiIg0CtbayEfn+0ZWbx8yxrQBNuO8ZyByIYW6WeBm\na63bR5RijFkE9AJaW2tro/S3frh9SvxlzJgxjBkzJtthiM9ou5BotF1INNouJBpjjK+TgmzfPrQd\nuDhK9/bAE8BsYCLwEYAx5lCgDbDaWrsz2K0Q2GatrQufgDHmHKAP8NdoCYGIiIiIiDiymhRYa/cC\nL0V2D3v60GfW2vB3FjwEjAACwIJgt1LgN8aYV3GeQrQX5+lDlwFfA6PSEryIiIiISJ7Idk1BPDb4\niexWF9HtE+Bd4BygA3Aw8AXwOPCgtXZdmuOUPBMIBLIdgviQtguJRtuFRKPtQnKRLx5Jmi1qUyAi\nIiIimaA2BSIiIiI+1aVLF1atWpXtMCSPdO7cmZUrV2Y7DNdUU9CI519ERKSxC169zXYYkkdibVN+\nrykoyHYAIiIiIiKSXUoKREREREQaOSUFIiIiIiKNnJICEREREZ869thjWbBgQeIBRRpISYGISJ6o\nrKyktLSUysrKxAOLSE74+OOPOf300xs0jbFjxzJixAiPIoIrrriCu+++27PpZdubb75JSUlJ/fc9\ne/Zw0UUXcdppp7Ft2zbGjh1L06ZNKSwspLCwkNatW1NcXJzFiNNDjyQVEckTFRUVVFdXA1BWVpbd\nYEREcogxzkOBamtrueiii9i1axdz586lefPmAFx66aVMnTo1myGmnWoKRETyRHl5OYFAgPLy8myH\nIiIe6dq1K2+88QbgXPEfPHgwl19+OYWFhfzgBz/gvffeqx/217/+NUcccQSFhYUcc8wxVFVV8dpr\nr/HAAw8wY8YMWrduzfHHHw/A008/TY8ePSgsLOTII4/kqaeeqp9O6Mr5b37zGzp06EDHjh15+umn\nAZgwYQLPPvss48ePp7CwkAsuuKD+t4888kgKCws59thj+fOf/1w/vSlTpnDaaadx8803U1xcTLdu\n3ZgzZ059/61bt3LVVVdx+OGHU1JSwl133VX/SE9rLePGjaNLly4ceuih/OxnP6Ompma/OGMtr3fe\neYdevXrRpk0bDjvsMH71q1/FXdY7d+7k3HPPpa6ujr/+9a/1CUFjoaRARCRPlJWVUVVVpVoCkQzK\n9G17r776KkOHDmXLli2cd955XHfddQD8+9//5g9/+AP//Oc/2bp1K6+99hpdunShf//+3H777Qwe\nPJiamhqWLFkCQIcOHZg1axZbt27lj3/8I6NGjeL999+v/53169dTU1PD2rVrmThxItdddx1btmzh\n6quv5rLLLmP06NFs3bqVl19+GYAjjzySv/3tb2zdupV77rmHYcOG8dVXX9VPb/HixRxzzDF88803\n3HzzzYwcObK+34gRI2jatCkrVqxgyZIlzJ07l4kTJwLwxz/+kalTp/Lmm2+yYsUKampq6ucZvr3C\nH83//M//cOONN7JlyxY+++wzLrnkkpjD7tq1iwEDBtCyZUtefvllmjVr5ma15AUlBSIiIiIpCt22\nV1FRkZHf69u3L/3798cYw/Dhw/nwww8BaNKkCbW1tXz88cfs3buXTp060bVr15jTGTBgAF26dAHg\ntNNO48wzz2ThwoX1/Zs2bcpdd91FkyZNGDBgAIcccgiffPJJzOkNHDiQDh06ADBo0CC+//3vs3jx\n4vr+nTt35sorr8QYw+WXX866dev4+uuv+frrr5kzZw6//e1vad68Oe3atePGG2/k+eefB2D69Onc\ndNNNdO7cmZYtW/Lggw8yY8YM6urqEi6rpk2bsnz5cr755htatmxJ7969Yw5bU1PD22+/zeWXX87B\nBx98QP8ZM2ZQXFxc//npT3+a8Pdzje+SAmNMC2PM58aYOmNM0nuYMeZsY8zfjDHbjDHfGGNmGmO6\npC9SERERaewyfdveoYceWv9/y5Yt2bVrF3V1dXTr1o3f/e53jBkzhg4dOjB06FDWr18fczqzZ8/m\nlFNOoW3bthQVFTF79mw2btxY379t27YUFHxbTGzZsiXbtm2LOb2pU6dy/PHHU1RURFFREUuXLt1v\neuFxt2jRAoBt27axatUq9uzZw2GHHUZxcTFFRUVce+21bNiwAYC1a9fSuXPn+nE7d+7Mnj179quF\niGXSpEl88sknHH300Zx88sn89a9/jTls+/btef755xkxYgSvv/76Af0HDx7Mpk2b6j/z589P+Pu5\nxndJAXAfUAwk/c5xY8xFwKtAM+BXwHjgdGCRMebQeOOKiIiIpMpPt+1deumlLFy4kFWrVgFwyy23\nAAfeYlNbW8vFF1/M6NGj2bBhA5s3b2bAgAH19/EnEjm91atX8/Of/5zHH3+czZs3s3nzZnr27JnU\n9EpKSmjevDnffPMNmzZtYvPmzfznP/+prwE5/PDD6+cHYNWqVRx88MF06NCBVq1asWPHjvp++/bt\nq08mALp168b06dPZsGEDo0eP5uKLL2bnzp0xY7nwwguZMGECgwYNqn9oQ2Piq6TAGHMC8D/APUDs\nm8T2H+cg4DFgFXCatfZJa+2vgf7AocCY9EQrIiIikl2hgve///1vqqqqqK2tpWnTprRo0YImTZoA\nTvuBlStX1g9bW1tLbW0t7dq1o6CggNmzZ0e9Oh5Lhw4dWLFiRf337du3U1BQQLt27airq+OPf/wj\nH3/8cVLTOvTQQznzzDMZNWoUNTU1WGtZsWJF/bsZhgwZwm9/+1tWrlzJtm3buOOOO7j00kspKCig\ne/fu7Nq1i9mzZ7N3717GjRtHbW1t/bSfffbZ+tqKNm3aYIypXyaxXHrppTz22GNccMEFvPXWW0kv\nk3zgm6TAGFMATABmAW5a6/wEOAyYaK2tT/+stR8A1cBgY0z8LUBERETEh+I1pA3vv3v3bm699Vba\nt2/P4YcfzoYNG3jggQcA5x5/ay1t27blpJNO4pBDDuHRRx9l0KBBFBcX8/zzz9c/RSiZOEaOHMnS\npUspLi7moosu4phjjuGmm27ixz/+MYceeihLly6lb9++SU9v6tSp1NbW0qNHD4qLixk0aFD9rU9X\nXnklw4cP5/TTT6dbt260bNmyvv1GYWEhjz/+OCNHjuSII46gdevWHHHEEfXTnTNnDj179qSwsJBR\no0YxY8YMmjZtGjcucBo+P/LII5x77rm8++67gNOmIPw9BYWFhfvdHpUPTLJVRelmjPkVTg1BD5xk\n5XPg99bauDfpGWNuBe4H/sta+0ZEv3HAbcCx1tplUca1fpl/ERERyTxjTNK3zYgkI9Y2Feye1J0w\n2eCLmgJjTFec23zGWmvXuBz98ODfL6P0C3XrmGJoIiIiIiJ5zxdJAfAE8Bnw2xTGbRn8uztKv10R\nw4iIiIiISISDsh2AMWYYcAZOI+F9KUwi1Ow82lsmmkcMc4AxY8bU/x8IBAgEAimEICIiIiLyrerq\n6px6ilFW2xQYY5oCa4B/AKPCeh0BVAHPAPcCG621W2JMQ20KREREJCVqUyBeU5uC1LQA2gPnAJ+G\nfapw3lMwHPg3MDLWBIB3cB5fekqUfqcAW4PTEBERERGRKLJdU3AQcH6UXu1x2hnMBiYCH1lrlwdf\nRNYGWB16/GhwGquAWqCntXZHsPsPgX8Ck6y118T4fdUUiIiINGKqKRCv5WpNgW8eSRrOGNOZKI8k\nNcY8DYwAAtbaBWHdLwaeBz7EeddBG+BGYB9wkrV2XYzfUVIgIiLSiCkpEK/lalKQ7duH4rHBT2S3\nugMGtPZPODUOu4D/BW4G3gT6xkoIRERERETE4cuagkxRTYGIiEjjppqCxKZMmcLEiRNZuHBhSuOf\nffbZDBkyhOHDh3scmT+ppkBERERE8pIxyZVlx44dy4gRI/brNmvWrLQkBAUFBaxYsaL++8MPP0zH\njh1ZtmwZb775Jk2aNKGwsJDCwkJat25NYWEh//jHPzyPI19k/T0FIiIiIiJuhScq48aN46mnnmLB\nggV069aNr7/+mo4dO7J69eosRphbVFMgIiKSoyorKyktLaWysjLboUiafPHFFwwcOJDvfve7tG/f\nnvJy5/krY8eO3e/q+6pVqygoKKCuzml6WVpayl133cWpp55K69atueCCC9i0aRPDhg2jTZs2nHzy\nyfUF5shxQ+NPnjw5akw33ngjnTp1ok2bNvTq1YtFixYB8Nprr/HAAw8wY8YMWrduzfHHH7/ftGpr\naykqKuJf//pX/bQ2btxIy5Yt2bhxIwB/+ctfOP744ykqKqJv37589NFHMZdN6BadO++8k8mTJ7Nw\n4UK6devmbgFHMXPmTLp168a2bdsAmD17NocddhjffPNNg6ftZ0oKREREclRFRQXV1dVUVFRkO5T8\nZYy3Hxfq6uo499xz6dq1K6tXr+bLL7/k0ksvDQtt/+lFfp8xYwbPPvssa9euZfny5fTp04eRI0ey\nefNmjj76aMaOHRtz3Hh69+7Nhx9+yObNmxk6dCiDBg2itraW/v37c/vttzN48GBqampYsmTJfuM1\nbdqUgQMH8txzz9V3mzlzJoFAgHbt2vHee+8xcuRIJkyYwKZNm7jmmms4//zz2bNnT8xYbrnlFl54\n4QUWLlxI586dk56HeC655BL69OlDeXk5mzZt4qqrrmLy5Mm0bdvWk+n7lZICERGRHFVeXk4gEKi/\neiz5ZfHixaxbt47x48fTvHlzmjZtSp8+fZIe/4orrqBLly60bt2aAQMG0K1bN0pLSykoKGDQoEEH\nFNqTNXToUL7zne9QUFDAqFGj2L17N5988klS4w4ZMoTp06fXf58+fTqXXXYZABMnTuTaa6/lpJNO\nwhjD8OHDadasGW+//XbM6c2dO5ezzjqLjh07HtDvyy+/pLi4mOLiYoqKiiguLmbnzp1Jxfn73/+e\n+fPnEwgEuOCCCxgwYEBS4+UyJQUiIiI5qqysjKqqKsrKyrIdiqTBmjVr6Ny5MwUFqRXXOnToUP9/\nixYtDvgeuj3GrUceeYQePXpQVFREUVERW7durb/9J5F+/fqxa9cu3nnnHVavXs0HH3zAhRdeCDi3\nMT3yyCP7FeS/+OIL1q5dG3N6zz//PC+88AJjxow5oF/Hjh3ZtGkTmzZtYvPmzWzatIkWLVokFWeb\nNm0YNGgQS5cu5aabbkpqnFynhsYiIiIisWTxcaUlJSWsXr2aurq6AxKDVq1asWPHjvrv69al/lqm\nVq1aAbBjxw4OOeQQANavXx912IULFzJ+/Hiqqqro0aMHAMXFxfX39ye6DckYwyWXXML06dPp0KED\n5557bv3vl5SUcMcdd3DbbbclHXv37t2ZN28epaWltGjRgltuuSXpceN5//33mTx5MkOGDOGGG25g\n9uzZnkzXz1RTICIiIuJDvXv35rDDDuPWW29lx44d7N69m7feeguAH/3oRyxYsIA1a9awZcsWHnro\noZR/p127dnTs2JFp06ZRV1fH5MmT+eyzz6IOu23bNg4++GDatm1LbW0t9957LzU1NfX9O3TowMqV\nK+O++2HIkCHMmDGD6dOnM3To0PruV199NU8++SSLFy8GYPv27cyaNYvt27fHjb9Hjx7MnTuXhx9+\nmEcffbS+e7wYxo4dS79+/aL227VrF8OHD+ehhx5i8uTJrF27lieeeCJuDPlASYGIiIiIDxUUFPDq\nq6/y6aef0qlTJ0pKSpg5cyYAZ5xxBoMHD+a4446jV69enHfeefuN66bhMMCECRMYP3487dq1Y9my\nZZx66qlRh+vfvz9nnXUW3bt3p2vXrrRs2ZKSkpL6/oMGDcJaS9u2bTnppJOixtK7d29atWrFunXr\n9rtX/8QTT2TChAlcf/31FBcX0717d6ZMmRIz5vDpHnfcccyZM4d7772Xp556CnBqTyLfUxB6Utea\nNWtizuPtt99Op06d+PnPf07Tpk155plnuOuuu2ImSvlCbzRuxPMvIiLS2OmNxo3TCSecwPz58ykq\nKvJ82rn6RmMlBY14/kVERBo7JQXitVxNCrJ++5AxprsxZpox5l/GmP8YY7YbY5YZYx4xxhya5DSq\njTF1UT77jDEnpHseRERERERymR+ePnQEcCjwEvAFsBf4AfBzYLAx5kfW2kTPubLABuBGIDIDW+Ft\nuCIiIiIi+cW3tw8ZYy4GZgKjrbUPJxi2Cuhsrf2ey9/Q7UMiIiKNmG4fEq/p9iHvrQ7+TboFiHG0\nTlM8IiIiIiJ5yTdJgTGmmTGmrTGmozHmTOBJnNuCZiU5iY7ANmCLMWabMeZFY8xR6YpXRERERCRf\n+KFNQchVwGNh3z8Hhllr/5bEuCuARcCHwD7gZOAGoJ8xpq+1dqnXwYqIiIiI5Iuk2xQYY0qstWvS\nFogxhwNHA4cAxwPnA1OstRUpTq8vUA3Mt9b2jzGM2hSIiIg0YmpTIF7L1TYFbmoKVhpjXgcmAi9b\na/d6GYi1di2wNvj1FWPMS8A7xpgW1tpfpzC9RcaYBUCpMaaZtXZ3tOHGjBlT/38gECAQCLiOXURE\nRCRfTZkyhYkTJ7Jw4cKUxj/77LMZMmQIw4cP9zgyf6uurqa6ujrbYSTNTU3B48ClQBvgG2AqMMla\nuyxtwRnzd+Bwa23nFMefDFwOdLTWro/SXzUFIiIijZhqChKbMmUKkyZNYsGCBQmHHTt2LJ999hlT\np05Ne1wFBQUsX76c733Pefjkww8/zG9/+1vmzZvH119/Tb9+/WjVqhUA1lqMMcydO5eTTz45rXHl\nak1B0g2NrbW/AA4DRgAfA6OAj40xbxljrjDGtExDfC2A4gaM3x3nvQebvAlHRERERPzAmG/L1+PG\njaOiooIFCxZwzDHHANCxY0e2bt3K1q1bqampYevWrWlPCHKZq6cPWWt3W2uftdb2A44EHsR5+dhE\nYL0x5iljTG830zTGdIjRvRQ4Fvh7WLdDjTFHGWNahHUrNMYcMB/GmHOAPsDr1tpaNzGJiIiI+MEX\nX3zBwIED+e53v0v79u0pLy8HnCvy4bfjrFq1ioKCAurq6gAoLS3lrrvu4tRTT6V169ZccMEFbNq0\niWHDhtGmTRtOPvlkVq9eHXXc0PiTJ0+OGtONN95Ip06daNOmDb169WLRokUAvPbaazzwwAPMmDGD\n1q1bc/zxx+83rdraWoqKivjXv/5VP62NGzfSsmVLNm503lP7l7/8heOPP56ioiL69u3LRx99FHPZ\nhK7G33nnnUyePJmFCxfSrVs3dws4inPPPZc//OEP+3X74Q9/yCuvvNLgaftZyo8ktdZ+bq29E6dx\n8LM4DYSvAv5ujFlijBmU5KSeMMb83RhzvzHm58aYcmPMFGAOsAX4VdiwDwHLgF5h3UqBT40xvwuO\n+4vg+C8DX+PUaIiIiIi4Zoy3Hzfq6uo499xz6dq1K6tXr+bLL7/k0ksvDYtt/wlGfp8xYwbPPvss\na9euZfny5fTp04eRI0eyefNmjj76aMaOHRtz3Hh69+7Nhx9+yObNmxk6dCiDBg2itraW/v37c/vt\ntzN48GBqampYsmTJfuM1bdqUgQMH8txzz9V3mzlzJoFAgHbt2vHee+8xcuRIJkyYwKZNm7jmmms4\n//zz2bNnT8xYbrnlFl544QUWLlxI584p3W1+gMsvv5xnnnmm/vsHH3zA2rVrOfvssz2Zvl+lnBQY\nY44zxjyK85KxYcAq4G7gNqAQeN4Yc3cSk5oObAxO43c4tQ+9gCeAH1prPwwb1gJ1EeN/ArwLnAOM\nAx7BqSF4HDjeWrs8pRnMMZWVlZSWllJZWZntUERERMQDixcvZt26dYwfP57mzZvTtGlT+vTpk/T4\nV1xxBV26dKF169YMGDCAbt26UVpaSkFBAYMGDTqg0J6soUOH8p3vfIeCggJGjRrF7t27+eSTT5Ia\nd8iQIUyfPr3++/Tp07nssssAmDhxItdeey0nnXQSxhiGDx9Os2bNePvtt2NOb+7cuZx11ll07Njx\ngH5ffvklxcXFFBcXU1RURHFxMTt37kwY4wUXXMDy5cv57LPPAJg2bRqDBw/moIP89CR/77maO2NM\nITAUGAmcgPNOgFeBCcBroVa7xphHcAr71wH3xpumtfZPwJ+S+X1r7RXAFRHd/g8Y7GY+8lFFRUV9\nC/eysrLsBiMiIiINtmbNGjp37kxBQWrXcDt0+PYO7RYtWhzwfdu2bSlN95FHHmHSpEmsW7cOgJqa\nmvrbfxLp168fu3bt4p133qFDhw588MEHXHjhhYBzG9PUqVN57DHntVXWWvbs2cPatWtjTu/555/n\nyiuvpKioaL8nSoLTpiB0i5QbTZs25ZJLLmHatGncfffdPPfcc7z44ouup5Nrkk4KjDFTgYE4jX8/\nB+4EJltrv4oc1lq7zxjzMpDsLUTSQKF7DEN/RUREpOGy+WCikpISVq9eTV1d3QGJQatWrdixY0f9\n91ABPRWhJ/Ts2LGDQw45BID16w94aCMACxcuZPz48VRVVdGjRw8AiouL6+/vT3QbkjGGSy65hOnT\np9OhQwfOPffc+t8vKSnhjjvu4Lbbbks69u7duzNv3jxKS0tp0aIFt9xyS9LjxjNixAiGDx/Oqaee\nSqtWrRpFA2U3qedgYBbQ31rbzVr7YLSEIMxbRFzVl/QpKyujqqpKtQQiIiJ5onfv3hx22GHceuut\n7Nixg927d/PWW28B8KMf/YgFCxawZs0atmzZwkMPPZTy77Rr146OHTsybdo06urqmDx5cv2tM5G2\nbdvGwQcfTNu2bamtreXee++lpqamvn+HDh1YuXJl3Me8DhkyhBkzZjB9+nSGDh1a3/3qq6/mySef\nZPHixQBs376dWbNmsX379rjx9+jRg7lz5/Lwww/z6KOP1nePF8PYsWPp169fzP4//vGPKSgo4Je/\n/GWjeb+Cm6SgxFo7yFo7N5mBrbUrrbVTUoxLREREpFErKCjg1Vdf5dNPP6VTp06UlJQwc+ZMAM44\n4wwGDx7McccdR69evTjvvPP2G9dNw2GACRMmMH78eNq1a8eyZcs49dRTow7Xv39/zjrrLLp3707X\nrl1p2bIlJSUl9f0HDRqEtZa2bdty0kknRY2ld+/etGrVinXr1jFgwID67ieeeCITJkzg+uuvp7i4\nmO7duzNlSuyiZPh0jzvuOObMmcO9997LU089BTi1J4WFhRQWFtK6dWsKCwvr216uWbMm5jyGjBgx\ngo8//phhw4bFHS5fJP3ysnykl5eJiIg0bnp5WeN0wgknMH/+fIqKimIO88wzzzBhwoSkXtoWLldf\nXuamTUH0h9V+ywI7cZ5GNNdam1qTdhERERGRNHrvvffi9t+xYwePP/44119/fYYiyr6kawqMMXU4\nBX+AyCwnsrsFngdGWGv3NTTIdFFNgYiISOOmmgKJ9Prrr3PRRRdx5pln8qc//cn1059ytabATVLQ\nFpgNrMB5F8AnOIX/Y4BfAl2AS4B2wGjgYuBOa+2DnkftESUFIiIijZuSAvFaY0gK/gi0t9aeG6P/\nX4Gvg+8SwBizAGhrre3pVbBeU1IgIiLSuCkpEK/lalLgpj7kPJxHksby1+AwIa8AXVMJSkRERERE\nMsdNUtAcODxO/yOCw4RsB/amEpSIiIiIiGRO0k8fwnkZ2Q3GmL9Ya98O72GMOQW4PjhMyA+ANQ0P\nUURERCQ9Onfu7PqZ/iLxdO7cOdshpMRNm4LjgIXAIcBinIbGAEcBvYFtwGnW2g+NMc2B94E/W2tv\nTTDd7sDdwAk4NREH4zzWdBbwv9ba6O/ZPnA6ZwN3AD8EdgPzgdHW2pVxxlGbAhERERFJO7+3KXD1\n8jJjzPeB+4EBQKtg5+04TyW601r7b9cBGNMPuB14G/gC55ajHwBXAluAH1lrNyaYxkXAC8ASYCLQ\nBhgVnNZJsRILJQUiIiIikgl5lRTUj2RMAdAe570EX1tr6zwPzJiLgZk4V/sfjjPcQcAqnNqBntba\nncHuPwT+CUy01l4bY1wlBY1QZWUlFRUVlJeXU1ZWlu1wREREpBHwe1KQVENjY8whxpg3jDEjAay1\nddbar6y169OREAStDv6N/f5px0+Aw3AK/ztDHa21HwDVwGBjTJO0RCg5qaKigurqaioqKrIdioiI\niIgvJNXQ2Fq7zRjTC3g2XYEYY5rhtFdoDvQEHsJ5OVq8x6AC9AoO93aUfm8DpUB3YJlnwUpOKy8v\n3++viIiISGPn5ulD7+O8vThdrgIeC/v+OTDMWvu3BOOFHpP6ZZR+oW4dUVIgQWVlZbptSERERCSM\nm6TgHqDSGPNXa21VGmKpxCm4HwIcD5yP024hkZbBv7uj9NsVMYyIiIiIiERwkxQMw7nPf54x5gPg\n38COiGG2djtqAAAgAElEQVSstXZkKoFYa9cCa4NfXzHGvAS8Y4xpYa39dZxRQzE0i9KvecQwBxgz\nZkz9/4FAgEAgkGzIIiIiIiJRVVdXU11dne0wkubmPQXJNCi21lrPGvUaY/4OHG6tjfkWCGPMrTiP\nSf0va+0bEf3GAbcBx1prD7h9SE8fEhEREZFMyIunDwFYawuS+Hj9lJ8WQHGCYd7BeTTqKVH6nQJs\nxanVEBERERGRKJJOCtLFGNMhRvdS4Fjg72HdDjXGHGWMaRE26JvAOuAqY0zLsGF/iPO40pnW2n1p\nCV5yXmVlJaWlpVRWVmY7FBEREZGscf3yMmNMK5wr8B2AedbarxoUgNN24DDgDZyXkDUHTgQuBbYB\npdbaD4PDPg2MAALW2gVh07gYeB74EJiA80bjG4F9OG80Xhfjt3X7UCNXWlpKdXU1gUCAqqp0tJ8X\nERERyaPbhwCMMf+N85jP14GpOO8TwBjT3hizyxjz8xRimA5sxGnI/DvgQZx3DzwB/DCUEARZ4IC2\nDdbaP+E8rWgX8L/AzTg1CH1jJQSSHrl25b28vJxAIKB3FoiIiEij5qah8UDgBeBl4FVgInBGqHGv\nMebPwMHW2nPSFKvnVFPgPV15FxERETmQ32sK3DyS9GagylpbZoxpi5MUhHsXuNqzyCQn6W3BIiIi\nIrnHTU3BduAWa+3vg0nBBvavKRgJ/N5a2yLedPxENQUiIiIikgl+rylw06ZgX4LhDwe2NywcERER\nERHJNDdJwQdA/2g9jDEFwCCcdwaIiIiIiEgOcZMU/B4YYIy5j29fKFZgjDkKpwFyT6DC4/gkx+Xa\n04gkOq1HERGR/ObqPQXGmHHA7TiPBS0I/jXBzz3W2vvSEWS6qE1B+ulpRPlB61FERKRh/N6mwM3T\nh7DW3hl82dhlwNE4ycCnwDPW2nfTEJ/kOD2NKD9oPYqIiOQ31280zieqKRARERGRTPB7TYGrNxqL\niIiIiEj+cZUUGGM6GWPuN8bMNMbMN8a8EfGZn65ARfKFGu2KiIiI3yTdpsAYMwCoBJoCNcAmLwIw\nxnwfGA78F9ANaA58hvNEo99Za3ckMY1q4PQovSzQy1r7nhexinihoqKC6upqAMrKyrIbjIiIiAju\nGho/CGwELvS4UfGVwC+AV4BpwB6gFBgHDDLG/NhauzvBNCzOG5ZvxGn8HG6Fh7GKNJga7YqIiIjf\nJN3Q2BizC7jTWvuwpwEYcwLwqbW2JqL7fTiPP73BWvt4gmlUAZ2ttd9z+dtqaCwiIiIiaZdPDY03\nALVeB2CtfS8yIQiagXPV/9hkp2UcrT0LTkRERESkEXCTFDwDDExXIFGUBP9+leTwHYFtwBZjzDZj\nzIvBty2LiIiIiEgcbm4f6g5MAb4GHgU+B/ZFDmetXd3goIwpAP4GnAAca639NMHwk4C1wIfBmE4G\nbgB2A32ttUtjjKfbh0REREQk7fx++5CbpKAOp0GvCf6NylrbpMFBGfMYTuPj26y141OcRl+gGphv\nre0fYxglBSIiIiKSdn5PCtw8fehe4iQDXgk2ML4OeDLVhADAWrvIGLMAKDXGNEviCUYiIiIiIo1S\n0kmBtXZMGuMAwBgzBrgDmGSt/YUHk1wJ/AQoAtZHG2DMmDH1/wcCAQKBgAc/KyIiIiKNWXV1df17\niXJB0rcPpZsx5h7gHuBpa+2VHk1zEdALaG2tPeDJSbp9SEREREQaqrKykoqKCsrLy2O+mNTvtw+5\nefoQxpjWxpi7jTGLjDGfGmNOCXZvF+x+dCpBGGPuxkkIpsRLCIwxhxpjjjLGtAjrVhhsmBw57DlA\nH+D1aAmBiIiIiJ9VVlZSWlpKZWVltkORBCoqKqiurqaioiLboaQs6duHjDHtgUXA94Dlwb8tAKy1\nG40xlwPfAW5yE4Ax5jpgDLAKeMMYc1nEIF9Za+cF/38IGAEEgAXBbqXAb4wxr+K8vXgvztOHLsN5\nUtIoN/GIiIiI+EGooAnEvPos/lBeXr7f31zkpqHxOOBQnAL3apwCd7iXgZ+mEMNJOA2YOwFPR+n/\nJhBKCixQF9H/E+Bd4BygA3Aw8AXwOPCgtXZdCjGJiIiIZFU+FDQbi7KyspxP3Nw8kvRLYKq19jZj\nTFucNxyfYa19I9j/euBea21x2qL1mNoUiGRXMvdgioiI5AO/tylwU1PQDue2oVjqgOYNC0dEGhNV\njYuIiPiDm6RgPdAtTv/jcW4rEhFJiqrGRURE/MHN7UNPABcBPwJqCbt9yBhzMk7D399Za29JV7Be\n0+1DIiIiIpIJfr99yM0jScfiPNlnCfAgTqPfy40xz+EkBGuBX3seoYiknR57JyHaFkTyU67s27kS\nZz5y9fIyY0wJ8HucJ/2EEgoLzAL+21r7hecRppFqCkQcpaWlVFdXEwgEqKqqynY4kkXaFkTyU67s\n27kSZyr8XlPgpk0B1to1wAXGmELgKMAAy621m9IRnIhkhu7tlxBtCyL5KVf27VyJMx+5qinIN6op\nEBEREZFM8HtNgZs2BSIiIiIikoeUFEhOUQMkERER/9P5Ove4alMgkm162ZWIiIj/6Xyde5QUSE5R\nAyQRERH/0/k692S9obEx5vvAcOC/cN6Y3Bz4DHgB52VoO5KcztnAHcAPgd3AfGC0tXZlnHHU0FhE\nclJlZSUVFRWUl5frKpyISA7we0NjPyQFDwK/AF4B3gb2AKXAYOAD4MfW2t0JpnERThKxBJgItAFG\n4bxs7SRr7foY4ykpEJGclM/P8hYRyUd+TwqSvn3IGDMWGGitPTZG/w+BmdbacS5jeAF4wFpbE9bt\nKWPMcuB2YCTweJy4DgIeA1YBp1lrdwa7zwH+CYwBrnUZk4iIr6lqXkREvJR0TUGw0D/fWjsqRv9H\ngJ9aa3/kSWDGHAt8CDxprf1FnOF+CswF7rTWPhDRbx5wItDOWrsvyriqKRARERGRtPN7TYGbR5J2\nBf4vTv9PgsN4pST496sEw/UCLM6tR5HeBgqB7h7GJSIikhI9plFE/Mrt04e+E6dfEdCkAbHUM8YU\nAHfjtC+YnmDww4N/v4zSL9StI7DMi9hERERSpcc0iohfuakpWApcEK2HMcYA5xO/JsGNR4HewF3W\n2k8TDNsy+DdaY+RdEcMcQFdsRETES/FqA8rLywkEAmoLIvtRDZL4gZuagknA/zPGPA3cbK3dAGCM\naQ+MB34MXN/QgIwx9wHX4bQlGJ/EKKFHljaL0q95xDAH0BUbERHxUrzagLKyMp1v5ACqQRI/SDop\nsNZOMMb8BBgBDDfGrMO5l/9wwAAzrLVPNCQYY8wYnHcNTIrXuDjC2uDfjjjtGsJ1DP6NdmsRAJ07\nd+a73/0uY8aMIRAIEAgEXEQsIiKyPz0ZStzSNpOfqqur65O9XOD6PQXGmEuAy4AjcZKBT4BnrbV/\nalAgxtwD3AM8ba290sV4oacP3WWtvT+i33zgBPT0IRERERHJIr8/fSjrLy8DMMbcjfM+gSnW2ivi\nDHcozovJVoe9j+AgnHcU1AI9Q29ANsb8EOc9BZOstdfEmJ6SAhERERFJOyUFiQIw5jq+ffnY3UBd\nxCBfWWvnBYd9Guf2pYC1dkHYNC4Gnsd5r8EEnMThRmAfzhuN18X4bSUFIiIiIpJ2fk8KYrYpMMaM\nCP77jLXWhn2Py1o71WUMJ+G0TegEPB2l/5vAvNDkOTBpwFr7J2PM+cCdwP/iPIloHnBrrIRARERE\nREQcMWsKjDF1OIXwFtba2rDv8TIca6315F0FmaCaAhERERHJBL/XFMR7T0Ep0M9aWxv+Pfg31qdf\n+kIV0bOcRUREJDkqM7gTMymw1r5prX0z8nuiT2bClsYq9CzniooKV+PpwCAiIrku2XOZznmOVMsM\njVXS7ykwxkwG/p+19h8x+vcGrnXzOFERt1J9lrNeDCMiIrku2XOZznkOvf/BHTdvNP4ZTuPdqEkB\n0BW4HFBSIGmT6ttAdWAQkXxWWVlJRUUF5eXljboQmO+SPZfpnOfQG8TdSfqRpMGGxsOstdNj9L8S\neNxa29zD+NJKDY1FRCQflJaWUl1dTSAQoKqqKtvhiEgUfm9oHLemwBjTCegS1uloY8zpUQYtBv4b\nWO5daCIiIpIMXRluXFQzJOkQt6bAGHMPcA/Oo0jjTgfn/QFXWGuf8S689FJNgfiJXw7yfolDRESi\nU81QbsrpmgLgz8BKnEL/ZOAp4O8Rw1hgG/COtXaN1wGK5JtYhW6/NAzzSxwiXlKyK/lENUOSDnGT\nAmvtB8AHAMaYzsCL1tqPMxGYSL6KVej2y0E+nXGoYCbZomRX8oka0Eo6JN3QOB/p9iHJhmQLxvlY\ngFaVt2RLPu5PIpJb/H77ULw3Gh/AGFNijJlsjPnCGFNrjOkX7N4+2L1XesL0F70URBqirKyMqqqq\nhAWTfHrpSmif6d27N4FAIOu1IZmk44U/JLvfiYg0Vm5eXtYVeBtoHvx7WKiftXaDMeYk4CrgHTcB\nGGNuA44HTsR518FKa+33XE5jJdApSi8LtLfWbnIzvURUDS2Z4JfbibwQvs+ks4bAj1eDdbwQEZFc\n4OblZffjPGHoWGAn8HVE/1nAeSnEcD/wDfAe8J0Uxgen8L8MGIfTKDpcTYrTjCkbhTU/FnYkvfLp\nntFM7TN+LIDnU3InIuIHKhOlh5uXl30NPGatvc8Y0xbYAJxhrX0j2P+/gYestW1cBWBMF2vtyuD/\nHwGtUqgp+Bz43Frbz+V4OdOmIN/vxdYOLl7QdiQikv9ytUyUT20KCoF1cfo3xV3NAwChhMALxpgm\nxpjWXk3PT8rLy+vvxc7He5Rz6f75fFz+uS60TgDdNy4ikufCy0TiHTeF+DVAzzj9f0x232h8MrAD\nONgYswV4GbjNWhsvkckZ4beShDLkUPd8kEu3WPjxFpXGTusk/VQLIyJ+kU+31/qJm6TgJeBaY8wk\nvq0xsADGmIHAIJy3H2fDxzgvVfs/nHkKAFcD/Ywxva2167MUV1rkUgE6Wbm0g+fj8s91WifJS7Vw\nr8RLRCS/uWlTUIhT8O4CLADOBObh3FbUG3gfONVauyvlYFJsUxBjWkOAZ4EJ1tprYgyTM20KRES8\nkOq9uG6TCdUsiIjsz+9tCpKuKbDWbjXGnALcBwzFecrPfwH/AR4H7mhIQuA1a+1zxpj7gXPiDTdm\nzJj6/wOBAIFAIL2BifiEF4U2FfxyT6q1Km5r81SzICKNXXV1df1xMBek/EZjY0x7nMRgg1eX272s\nKQhO7w2gj7W2eYz+qimQRsuLpzfk6hMg3MqX5CeT85Evy0xEcp9fjkd+rylw9UbjcNbaDdbar31e\nqj4S+CrbQUh26WlB0Xnx9IbG8gSIXHo6VjyZnA+9QVgyQcd3SUa+HMPTzfUjRI0x3XEK22058EVh\nWGunehBXrN8uAVoCy621+4Ldiqy1m6MMex1wBPCHdMUjuUG3MUTnRePuXGog3hD50pA5X+bDr/xy\nNdJP0r1MdHyXZOjYlxw3DY07AFNw2hFAlIQAsNbaJq4CMGYY0Dk4veuBg4HfBHuvstZOCxu2Gjgd\n6GKtXR3s9j/ASGAOsBIn0SkFLgA+xbl96JsYv+3zig7xgk7Ukg+0HftfY7mdzo10LxPtF+I38bZJ\nv98+hLU2qQ/wArAP+D1wEfCTaJ9kpxc23argdKN93ogy7F6gU1i3PsCfcRKC7TjvKlgK3A8UJvht\nK43XSy+9ZAOBgH3ppZeyHYpIQoFAwAI2EAhkOxSJwU/HFL/E4nUcfpmvfKPlmrrIZRfvWB0sd7oq\nJ2fy46bw/h/gD9kO2NOZV1LQqOV7IUsH+fySifXZ0N/QNpc5iZZ1vh7fkpkvbYfu5ev2kgmRyy7e\n9pdPScFW4OfZDtjTmVdS0Kil68QRa7qZPlHpIN+4pbK9NXSb0TaXOYmWdb4VjEPzM3r06ITzpe3Q\nvXzbXtLBi3N7PiUFf1VNgTvayRqnWCekTJ+o3JxEJTm5tE+nsr2ppiB3NLZl7WZ7bmzLprHNb7Z4\ncQ7Pp6TgKGAtMDDbQXs28ykkBW52Pl2tSE2uH+D8UlMQou3QO7m0LHN9P/ILvyxHv8SRLY19/uPJ\npeNSLvNiG8ynpOAN4N/BBsBrgDeD3cI/87M9Q65mPoWkQFcr0k8HOG9pO/ROJpel1ps/+OV45Jc4\nxH90rMgd+ZQUrAQ+T/TJ9gy5mvk01xRIarSMRb4tBPbs2TMv94dc2c/9Eqdf4hAJp+3SHb8nBUm/\npyAf5cp7CvQcZpHGJ7Tfb9iwgaVLl+bds+/1TP8D6VgvuUb7sTt+f0+B6zcaS+bdddddLF26lA0b\nNuhEIZJlmSq4hd4WHf57+URvGD2Q3s4ruUb7cX4pyHYAflBZWUlpaSmVlZXZDkUaET9ud36MKZps\nxhkquFVUVGTk98rKyqiqqsq7QmK+zldDlJeXEwgEcrqAlSvHEPFGOvdjbUtZEOu+IoINh4GDwr4n\n+uRkQ2O/N+DSPXv5yY/bnR9jiiYbj9v0ejoi+ShXjiHif/m4LeHzNgXxagq+B3QFTMT3eJ/vNThL\nyQI/XZ2Jlhnrilrq/HylwU/bXYgfY4omPM5k13GsK/xutxHtj9kTvq78vG83ZrlyDBH/y/VtKSeP\nUdnOSrL5wYdvNM7HzDibGsPydHPlOh+vcie7jmPNey5vI/m4PuPNU/i6SuWNvvm4vBqqsR8/8pnW\nV+ZEW9bRjlH4vKbA60J2YYrj3QbMBD4D6oAVKU5nBPAesANYD0wA2sUZPuYKzpZ82Yn9Mh9+iSOe\nhsboplCbywXgWBrzW3jzcX3Gm6fwdZVovUWbTj4uLzeSLbjEko7ll8v7n9819u09k6It62jbdt4k\nBcCjCfq3Bt5KKQgnEdgAvAZ8k0pSAIwKTmc+cBUwBqgBPgJaxBgn7kqW1OlglLyGLitd6Wu88nF9\nprPtRy4uLy9jTrbgkolY4sUk3vDr9p5qXH6en549e9qePXsmjC2fkoI64OYY/VoCi4AdKQUBXcL+\n/8htUgC0BbYBfwfn3QvB7ucG4741xnhxV14u8stO45c4coGWlaSDtqv84GWh2Y/bhB9jSrfGOM/h\nUt2m/ZpAuokrn5KCO4B9wGUR3VsA1cAuYECDA0otKbgqGNvQKP2WAx/HGC/hCsw1ft1pxP8a+4kq\n30QeC7R+c1O211u2fz8fNfbzdD7WFCQbV94kBc688Hiw8H9G8HtzYB6wGzjPk4BSSwqeDCYF34vS\nbxqwF2gZpV/CFZhr/LrTiP819hNVvok8Fmj9Siq03Xgv0Xla5/HU+X3Z5VtSYIBKYAvQJ9gGoBYY\n6FlAqSUFrwSTgmZR+v062O/IKP2SX5MpyuQG6vedQfwtl7cfv8fuh/j8EIPkHm03madELHV+X3Z5\nlRQ480PzYPuBvcGE4BJPA0otKZgH7IvRb2wwKTguSr8kV2PqMrmBZuq3dJIQv/HjiSB8P/FjfBKb\njnGSTdr+Uuf3Zef3pOAgYjDGnB6rH/AI8DTwPLA+fFhr7YI446XLDgBjTDNr7e6Ifs3Dh4k0ZsyY\n+v8DgQCBQMDTwEIv3fDq5RuVlZVUVFRQXl5+wMuTvP6tWEIvgdqwYUPMWMS9eOtW4svUtu9GaD8B\nf8YnsYWvu7KyMu2bklFlZWXazlLkt2VXXV1dfyzJCbGyBZyn9uyL86mLGKaOGFfr3XxQm4K4/HDF\nMZSJ9+zZM+ux5BM/rNvGJp1XldJ9xcrvV8Rymdpj5D7tH9/K9WWR6/GHw+c1BfEK2pen8mlwQKkl\nBSODScllUfp9CiyNMZ7L1Zl9fto5vI7FT/OWDY19/rMhlwt7uRx7NH7e/v0cm0SXb/tHQ+T6ssj1\n+MPlbFKQtYASJAVACXAU0CSsWztgOwe+p+C8YLJwW4xpJbkaJVkNOXn6bcdP50uUsskv8WSrEX6s\n/3NNLsceTb4/j18yS9vAt3J9WeR6/OGUFCSXCAzDeQ/CncB6nLca3xH8DIsYtjpY0O8U0f2m4C1E\nbwBX4zQwrgE+jnbrUHCcvNnQ/KIhJ3a/7fheFVL8luz4JZ5sNcL3y/zL/rzc/0PruGfPnr46pohI\n46akILmkoCpO24U3ogy7NzIpCPYbASzBaVS8HpgAtIvzu3ELB34rpPpJrGWTT8sscl7y5YUrmV53\nfthWotUOjB49OqPrxavtKVdlY337vd1TNo8pjW37k8SibRPaTrylpMDHn0Q1BbqiGFvksknHgcPt\nNNN98Mr37SFd89fQ6aZrvbqNq6FxRP5eOreneLFm+ra4bBbQs93YO1H/VLcBL7adfD+eNRbpqGEL\n3yZU6+YtJQU+/iRqU5CpDDkXM/FMPJ3D7TTTfZLLxfXkRqZrCpKVrvXqNi6vk5uGLJeGFDYzfVtc\nrEJFPuxPiZZBov65VFOQ7QQr0fD5sD2lwssLdPFqCvxe65YrlBT4+OOXhsbZuGLj9QG0MdQU+E1j\nmV8v5zOdBfFMakhh06tbBNzWFEQOlw9XqpOtKcj0rWrp0JD1lcy2Epp+mzZtklpODa1589P+3BCZ\nuEAX7XckNUoKfPzxS1KQjZ0tkydkHUzSIx8KVZnm1TLLxDadiVuAQrJRoAqfRj4VnqPJh321Ies8\nmfl/6aWXbJs2bZJeTsnWFORzUhqNzrcOvy4HJQU+/vglKUhFQzf4TO4w+XrwzbZMH/T8epB1w6t5\nyMQ27efE3evYwq8Sp3IrhN9vJfFbPJmWbO1UOpZTrG3VzW/5pbYyl387HeLNj98uAIWmo6TAx590\nJQWZ2PFyqaCdbwciP9I2F126lku2awqyLV23H0bWFCS7zaVyK4mfl6+fpGs5paN2Kl2JhpfHvmwe\nR3PxGB5PvPnx2wWg0HSUFPj4k66kIBM7nk5o/pSpq1+R8mmbS/Y+7VTuO5bsSHXbGT16tG3Tpo0d\nPXq0q+m7uYdd20Z8frlHPZk4/BJrOqbl9S17+SCXLs6opiAHPpmuKUjHPcL5tpNHyrUrvdFOSn4p\nsKfr6m6mryC6WZ75vn+kgxe34kSOk+o+kM59x4v5agyyfdyIVYuUiVjdxJfuc7eS2PygpMDHn0y3\nKYi3U/vxpJkOfjkQZvKKktt7V3v27Gl79uzp+UnG63lO5zKMtwyykQClyi9xuJHKrTiJppEvFz0S\nLQu/xZtt0ZaH2+3J7+e4hp67k33aUviy1HaWHpmoEVJS4ONPppMC1RSkdh9pz5497RFHHJFSQTne\ndP243MLvO3R7JTzTBZaGTC/dCUymChKZTsQyIR01BYm654JkEna/JQ1eJdCpxh1teaRaU+DXbaYh\n5243T1sKF1qufn2hWEMvjmWLl8frWNNSUuDjTy4/fShXpXJgcFtQzmRs6YjBTU1B+IHHD/Eny23h\nya8FCb8VAt3KdHx+S5LczH8ysfstSUzm97waJhovEsx81tCLgX59oVi07SXeNuSX7UI1BUoKItej\nZ/yykeeDhtxS44YXJ/1My0Q86ajFcjvNTBemko3Pb9uDW35drpn6DTfz70Xsja2mIJLfksJc59c7\nDNzWFKRzu/CyptMLSgqSK5wbYBSwDNgJrAYeBlomOX5djM/WBOMlvSLdauhGnuuFjVwUWubxGrRl\n8qTml20g1Xn2clk19iva6eKXbcwNL6/GezX/ubgcU9WQefVbAa2xyubxLdO3Frmd14bc1pUMJQXJ\nFeofDRbiXwBGBhOCWmBekuPXAdXA0IjPoATjxV15mT74hfOqulbc80s1p18SS79ejUqnbBRe/La8\noj0O1A8xpnLLVrrjzsckMtYyy/TxMR+XbbZlcz9OpZCeyeNqKL5kG4C7jUFJQeICfQ9gHzAzovv1\nwcL+pUlMow6YnMJvx12RyRz8knlMWipS3RFy6QDqhwJGNH6JK1uJpaSmocs7NL5fGg+GrpaFTo5+\nuYc5lf0i3fuCX44ZXoq1zGLNa7qusPpt2b70UmZuZ02F35ZVNMm+eyQkW7c4RpbtGrJsw+dBSUHi\ngvm4YFLQJ6J7M2Ab8JckplEHTAYOBlq5+G1nje3cae3nnx+wIiM3gvCNJfyEme2TZLhcOChYm/4q\nOrex5MIyS1a6E9Z8l62akciCt5dXqlIRfvL2W8Lidll7vY/HOjfEm366jzPpnsdE/dJxhdWPQvPp\nh3NXJD9dCEqlpsnNdNItMs6GLFvVFLhLCuYAe4CDo/RbBHyVxDTqgK3B6dQBXwEVQGGC8azdvt3a\nkhJnUTz2WNwVG37QC/31a8HL7wVdr04gXsynnw6kXsjl+Yl3FTJXbtuKJ9nCo18S5hAvl38+7LOp\nFBiSGcarq5HJ8rJG2u/nHK9EqynI9Lz74TiZiNuapmxLlOh7FbeSgsQF+g+BdTH6zQjWIhyUYBp/\nx2mofD4wDJgeTA7eJ05jZcDayZOdxRD6xJGoWslNguDXe1wzdfL36ne8KBz49SCVqlybn/B4Y61P\nN+vZqyv26Vh+8eYj/HfjHVP8vn4TxefnfTbZW0MaWlOQ7FVUN/OZyjLxw3kil3h19buhEtXc+WH9\n+Ok4HG9aoRrRkpKSjKxDJQWJk4LlwMoY/aYEk4K4V/xjjHtbMDG4Lc4w1j7++P5Jwc6dSazW/UXW\nICSzUaX7IOLlFaBUZeJA6YeDnzRM+HbixRWwbF9FjideoTNa3Ml285NE8aWyz6Zy4SWV3wzFnu7l\nm+xVVL+eJxqrdF39djt+aPhYbXxC3Xv27Nmg38/m9pGpskio3NaiRYu01AxEUlKQuPDe4JqCGOMe\nBOwCFsUZxt5z9tn2HrD3gK0Ca995J4nVuj83J6yGntzSzW+3Cbj9Hb+c5NJ9hS/XpXOd+X15uilY\nJNvNTxoSX6IrsQ19A2yi8SKTtnTWSCQz3UTD+X1byDfpWt5e19i4TQq8qKH1Wqx5c9tQOd60wqd3\n4YUXpiUhr6qqsvfcc0/9R0lB4sJ7g9sUxJn2CuD/4vS39vDD7X41BSNGJF7LDZCNncyLA1mmTz5u\nTmmdQekAACAASURBVJrhBYVUlm865s1NHJnaJvxUgEjXPPtpHmPJhRizJdZtEemuKfDLbSFu+T2+\nWBrrPpCpK/JuLyZks6bA7b4Z3qbTS5H7kmoKspcU3BesDTg1onvSTx+KMd1mOO86eDPOMPsnBEm0\nK2ioTO1k4Ve74p043Farp+PkEy2GZH8vNFyLFi1sz549Uyo0eH2VJtSvZ8+e9ogjjkjYIC3VbcLt\neH4qQKTrCmj4PHpxEszFmjO3/BBXZKE/U48+Df1urN/zw7KJJdYxxg8SLbd0FMDScWz1ukY/m8dg\nPx3/wyWKK7J/KjUFycjUvq6kIHHh/dhgUvBCRPcbgt2HhHX7HnBUxHDFMab7v8Hxfxnnt6MmBUNP\nPtnVhuG3E0doJ0qmcNTQQpSX8YYfFJItNHpRiEhXATRyPcSa11S5nZbfttN4vEjU4k0jVr/I7snE\n4YfEuiHc3mqQDpm6Shfrd/3yqFU3QrG7vaUqE8s20bYeGYMX+0a0aaQ63VSXbSLZPAb79fifrgtE\n6YqnocMrKUguMagIFuBfxHmj8SPBq/zzI4ZbCeyL6PYb4C3gfuAa4JfAfJxGxn8DmsX53ahJwQMu\nDwLhBxA/7HCRNQWJhk2mEJXM+KkOG9491QJWNg54iZZdtPXQGK48eyFdVw5j9Yt1ZTCZOJLdb/y6\nvhqSFHh1IaGhyybZ8dNZE9SQuBoybbdXszORBKa7YJXsNLyqKUilxsCv+7uf4vJTLCFuy0Fuh1dS\nkFxSYHAeKboM2AmsCV7pbxkx3OfA3ohu5wOzg+PsAGqA94BbgKYJfjdqUmDBvvTii0mtYGsPvK89\nl8QqkEf+7+aJKbEkM2yuFrD8EE9DT4B+WZbJ8jJut1c20xVHQ2Rj/adSG5OOmJL9rUzX1mT695JZ\nbrGSAjcXk7ySrX3H6+0r2rLza82gV3F5se6ysYwaWjPR0AsLSgp8/ImXFEzq0iUjVwUasmN58Zux\ndsrw7qH/ow0Xq+V+qvH6pYDlF26Wh9sDbGjaqdx65Yf15OUJJdH8+PUEHy6TyyOZ4dKVbKQr3nTI\n9PE9meUWa7rxjvPpkq39yuvtK9qyS8e2lq4alVSG82LdxboQ2RDpjjt8/FRiVlLg488BScFZZ+33\n/ay+feOsWm80ZAN1M26sRCBW1hteXZpMTUH4EwH8UKBvaDbvF6muYzfTTuV+aj8UkjO5TlP5rVwq\ngEbK9vr1+7rNhFTWQUOTELc1BdF+z00M2a4p8KohcTpqWaItm0zul8nUjni57ryat3THncxF1XiU\nFPj4s19ScOKJ1m7YsH+SANbOmxdn9TZcpq4kJZvdRm7ksYaNPKiGagwydcByezUglQJwOmt/4p1E\nEl098fKKSqonsmyczNO5PrzUkBoYv0j3MvNTQTzbCVAs6VxGiaad7G9HW3Z+XZ7RJHu+S5aXyUa0\n5ZiLybJXNRNe/54XUpk3JQU+/gB2AiPtJ3zf2iVLnLX34IMHJgZgbfv21j79tLUvv2zt+edbe8cd\n1u7aFWMTiM7txurlxp3qjhnrAO/2gNXQA0OyccUaPpWCWqont2TGCw2TyrL18qQbb1p+KrhZm/pJ\nMtOFlIbUwDQWfio4er2dJ3tRIBMXg2JJtPyTXT/JXLTwMpn3+iJJvPNKKtMNje/FBTK/HX8TiRWv\nn/b1TItW7vB7UnAQjdzVTATg+knw2GPA6NEwZw68+eb+A27YAD/72bffX3kF7r8fTj4ZmjWDJUuc\n/7t1g+9/H77zHdi8GZo2haOOgu9/n98/+ijVwemWlZVRWVlJRUUF5eXllJWVHRBbRUUF1dXV9cOn\nIvw3qqqqEg5fVla232+Vl5fv9zde98hxwyU7L7GGi+weK65YsYS+hy+PSJWVldx1110A3HfffTF/\nI9F6SxQbdXWU33ADGzZsiDpctPHD5793794sWbKE3r171/e3FvbscT4Ar7zyCk888QTXXvvfnH/+\n+fsNF+6qq25i796WXHXVtWzZArt2wd69znQeeuhPLF78BTt2vEj37mXs2AH79jnTOPhgZ7Nv2hRe\nf30u06Y9zZVXjuCcc/rTvLnTv6AAjHH+AtTWfvv7oWw7xJhvP6Hvu3dDXd230/jZz37Fnj2tGTny\n52zZ4nT7zW+eYtGif7Bv35P071+G3VeHrdmGrd1DnTVYC9decg3NdzRnxKCfsfHTzdgmB2F312Jr\ntlG3czfWFGBNAdULFzHjhT8xuOwi+vbqzaK/vcWfX3mV8y+8kL6nnYats1hTsP+MBT/GQJMCS5MC\ny9UDrqSwpjmXDR/KGWeeCU2asOHfm7GmgIKDm1DQxGAOakLBQc73v8z6C48/8STXXPMLzjvvPGfB\nhFbCnj1g66JvR/XbU9jCDC3A2lrnb1iMFBgwBRHdvv0UHFRAnTXs2gU7djiTCG1P4esl9Jk3by7P\nPDOVESOGc+aZZ9av07lz5zJ16jMMHz4Ca2HatGlcdtkwfvrTM+rDPP/8O9m2rRvnnTeU998/cHsI\nF/7b8f5GjhO+ikL/G+NsU3PmzGPq1Ge57LJh9OtXxu9+V4a1zmE82lWheDGFPvv2OavtvvveYMmS\nvaxc+QpjxvyNCy8s49RTT60f11q4556/89FHzfjmm7/TsmVZ1OkCLFq0iD//uZKysjJOO60vAGPH\nLuKDDwybNy+kTZuyqMth0aIFvPjiSwwceBGnn3468G18+/ZBnz4PsHHjq5xyynm8+uqB83fKKQf2\nDy2Lujrnb0EBNGlSxi9/WUZBAbz2mtOtsLCMe+5x4lq4EO69t4r339/Lf/5TTdu2ZfXT2rfv25j2\n7nWmW1fndKurgzFjPmLZssNZv34Z1dUlvP7669TU1PDll0ewfv0yduxwfnfcuA/517/asnjxS9xw\nw2GceuqP99u0Q78X+t6kybf/L1tWzJYtPais/IJx467h2GMv5MQTu3H++UN44IFpvPvu12zf/gpH\nHVWW1DZ45JEDePfdLfTqdTJbt27lkksu57PPkt+Gw//26lXGtGnO706aNJtJkyYxcuSVnH322fsP\nby3G1n37N2w64dvvAf/X2ajdQ0LLLrJH5P4WGv+3Dz3LPxavYO/26Zxw7DlgDKZJAUOHjGbXrvYc\nc8yPOOWUQYwceSUDBgw44Hgf+VOzZs1i0qTJjBw5kgEDBsQ8PsTqHs3s2XP+f3tnHmZHVSb833tv\n315uku4QdtIxAZEIAUJkCSGBgMuAASQow8B8IyoKCYR1/EQENQou87EMH4oQgl+COA4fW4CsiCyR\nAEICIiYMDiKgkLCGpJekl9v3nvmjqjrV1VX3VtWt27fbfn/PU8/tPvXWqbfOec/ynjrnFIsWLeLs\ns88G6P37hBNO6PMs3jSL+/+PfrSM557bzD777M8nPjGD448/k9Wrjw2vcBUQEyVF/84QEWMN1lqs\nXg0zZ9r/LFkCX/iC73UvcQCb2Iu7+Cfu5Ezq6OKTPMYsVnIMT7APr9FDDXdzOm2Moos6HuYfeF32\nZoRsZURtgb0yW2huW89Wunm/bgId9ROply4aG6FxlzrGjdlGvuUVNrz5GpMPnMi+21rofPlltnR2\n0DbmIDY2TGGLjGFj1y48/+5Y6tM59sy2steIFhprO+kytWztrOev73VQyLUwJrudQw8eS7quhrps\nClNTS6EAKZMnXQNIiu35OjrytTTUG+qlizGZNhrSXbTmshTqs7R01LKtq4aGrDCqUUjVpNjaXkN2\nZIqRo6yOTken0NFhVSjZrNU5FIFXXvkTTz+9gUMOOYSDD96X9nZ4/90C5HtIFfKkCjm6t+VoeedN\nPtz4X3x83z2ZMG432rpr2d7Sw+aNb9Py9jvslG0k1S30FNIUSNFNLR2mHqFATw4KPXlEhEKqhrpU\njpGpbYys6SSTNmTSeWoKOavDSR35dIaeQpqufJrthXo2b4O2QgPd1COSpqGuBhGDtTmWU9KFXFcX\nxhTISJ7G2jw1hQ5S+Q5G1OSpTxfozGfopIFMOo8YQ1c+TWc+Qxd1dFFPV6GGVhrZwk4USJHCqsh7\nqKGAVRNnyPWG10ieOrrImC5qJUfKdNNJA+000i315EmTM5mKlBFFURRFUZJCMMb4DGUMDtQpcDkF\nBxzwAbvt9o9cdNFF7LffqUyZYsjlhAa2UyDFWdzBbZxbRY0VRVEURVGUoYk6BYMWr1OgKEOZGnJk\nyCEem47yfx1d1NJNhhxp8qQoYBAy5BjBNtLkEQw5MnRTSxd1CIY8aXpsDTqpp4caDEKBFAVSGIRa\nuhFM7/3cfzuyzt+OfJp87zkn3B2n+1d6JQwpCjv+F8AUXOdAxJAS65z7OkfeG+7gvL2xggymYCgU\nCpBKU5Aa8qTJm5SVHibde23K1rpgLM0KvSGp3ueuoadf3jg51BfT/5xbxJheCfsdl33a+a8/7jSs\np5Ms2623U+T6yAQdfnnqTb+w/3v1CvvrvtbPXpywNPk+b+KCn8r009F7T+cAK49r6OktIzX0UEPP\nDpuxr3Fs0x2X3zOVCot6jaNbmnw/vUrhpK+Tkk6Yk6550n3+zpPuY1MO3rKUIddHLyd+99/u9HLr\nY5B+93Lr4fzvTmvveb9r3GF+6TmUfouVsWJ/u+3am+47/rak/epcdx4F6VasDAfZjH/9GBweRzZu\nvVXsf6cuSJPnWaYNaqdg2K8pSIJ99rmFOXPO48EH4emnd4SPGWPNX5wxA4wxfObg95gweisL//3n\nvPhWgZrdpjN2j2M5cPTvWPvfb3HoJ/bliL3q2fJmO6++O4r27lrefP1vdHT2UGgYybjJh0B9PY2Z\nDj46ejN7Zj7gud/cy87vr2D8njtz2GkXsmn7aNq7MtTTSUNdAcn3UEcX9y97mvc2w+iGLBlpYPem\nRppGZilk6ujJWxNtC905GnMfkjcpPsyl+Ut7hsxOYxnbkCPduY1G08JI08b2XIb2nHXd6MKHdPWk\nac1nkUKeBjpooIMCKbYxorfCzpNmJO2kKLCdLDX00CybkNoMpq6eQm096RH15OpH0dqwO9t66ugo\n1DKKNkbU9lCzxy682d7OH157hckzp3LwkYeRokBtxtCQ6QERUtl60ruM5qnVT7Dqwfv52ISJ/OnP\nbzPjyE/z8X0/Tq7L0GPSFIxQn86RNj2k01DXkCKbhWx9gRENBbJ1VqOe7zGsW7eOhx96iJa2Nv66\n8S0m7X8A3/3hDyGTIdcjdHSl6OoWunIpawp4j1A3ooa6dA+5zjzGQH1THc9teJ7lS+/hjNNO4jOf\nO54RozPsPDpPupDDdHVT6O4hk4G0FDCm7/z7FSsf5vbFv+KMM7/EjKOOo7sjT0Omh8b6burTOeZ+\n7Su8+PtnmXroFBbedJM1Gb+rC3I5rrj8cl7asIGGbJbW7ds58OCDueb66y3DdI5UCmpqehcJPPKb\n3/CrO+7g4AMP5KX16znji1+05sb31PWdCJzv4pJ5Z7H+D39g/PjxNI0ezamnncYxM2fy2zVruOfe\ne/lw61b+8vrrHDplCjffdNOOicPuCb6ODn5/O5O13ROinYUGnonfD69bxzV33cW8Sy7pt9bDWQdy\nxBFHsHbt2sD1IElRat2Jw3HHHcfq1asZ3dhIe2srxxxzDI+uXu0/UT7kPd3P6KxFOfbYY3v/D9Sp\nd9J4AxSaevNq2QMPcPMtt3Deedb6lKVLl3LLggWcd/75fO6UUzj++OP57Zo11NfXM2Hvvfne97/P\n7FNP7Tth2G+Cvncic4xzy5Yt4/zzzqOltZWjZ8xAgDVPPsnR06ezfPlyMIaTTz6ZJ596iqOnT2fp\ngw+yYtky/s+//RsCfPMb32DWZz8LxnD6aafxzDPPMHrUKK695hruuP12nn32WQ6fOpU7774bUilW\n/vrXXHzppWxta2Pa9OksdU+4d+zSeWb7d/ny5dx6663MmTOHk046qX96AyuWL2fhwoWce845nHji\niX3Ty5nE7/3fu7ChpgbSaZY/9BALbruNQw47jLXPP8+8Cy4A4Gc/+xnz5s3jlFNOCZ6IHSH8hBNO\nYM2aNRxz9NGsWrUqdDyzZs2y8mjGDFauXNlPftmyZSxYsIC5c+cCcOstC5g7Zw4nzZq1Y9GBO/3c\nNuJaTLBi5UoW3nYb586ZgwFuXbiQOXPnctLJJ/dPO6/OPnZ30okn8uRTTyHAjOnTWbZ0qX85Dci7\nVStW8PVLL6WtrY0jp03jnnvv7bfgZfnKlfzoxz/GAFd8+9sgws1OWZs9u+/iGOdXpGR943c+bB3l\nF9d3vvMdWlpaaGpq4uqrr+bUU0/trcuamppYvHhx3/vceKN1n9mze9P0gfvv56c//SkXXnABs085\npX+6+9UBPmHLli5l3rx5tLa2MmP6dJYvW1by2hUrVnDJxRfT1tbGUdOmseS++3zlVq1axeJFizj7\nK1+x1h2E0O/hhx/m/95wAwCXXnIJv7zjDp5du5aphx/OHb/4BXJA6KSuDtVe6VzNA+jbqyhy1Na+\nZyZO/MDsuuuj5rOffdW0tRXbR7jOHH30Z0wQUXYVKLXbQiW2+4q9i1B3tzEffGDMxo3mn6ZONfuB\nmTxypKkHc9j++5sTp0830ydONFkwOzU2RvpqdBSS3AEiqZ1kou7AEGW3oTD5FWeLvDC7/QTt6uRc\n29zcbJqamsxll10W+r5xibLrU6V3woi6e8tAbGEY99m913n/X7Jkxxfd49ppXNxlPag+DNphxi+t\n3HVGUD0b9Tmi2mVS24VG1SEqcXcZKiXn1rUcvZN+5iVLdmzhHLe8lirvQTqHrSei2EZS9YH73uW0\nuXHrB28dUCp+9zbqcesjty146xpvXeh9LqvbXf3+b9BRdQWq+vC2U/Cpmd1myRI/Z2B25M5gOY1H\nWPkwnbWoON8ZcApKUMEOe29vAXQ6jZMmTarI9wzc95s0aZJpbm42zc3NZX9Mpljhr0TeOkTpJJS6\nVxJbAgZ1/t2y7nR30n4g9+oP85xJpU0pPQbiQ0be+xWz93J0CpNmYeKvRMfUXW+FpZiuQelcju5R\n865YJydOxyupvA96jlJ1Q1SitKFhBkSSLNsOxewhjP6lOv9B8qXyPoqdxk27sOfitnnl9LnCxJ/k\nYKE3Hie8oaEhsLypUxCucy7ApcDLQAfwN+A6IFvJ6x2nYPr0N0uOLsUhbsEvZax+15fb4Hq/SByl\nQghz76DKIk6lX0wH93Mk1QHxe74wFVicTmqUa8PqWyp/wtwrzPO6K8q4Dk3SlLp/JTqqScfp7gj6\npX9QA1VKpzgd6iD9wpThpL4e66YS+edHuc5kqTrf3ZEo9iHIYg5DEO4Bmai6utM3KK0r3Zkrdn3c\nej7KvYu1QX427U4nJ+2bm5sTsZ+k3hSUauOT7DiHbQeSdi6D4k+iHgpytMM8pzoF4Tr1NwIF4B7g\nq3aHvht4pJLXO05BOr2t3+vwcgqwQ1KjBOXcK+y1zc3NJpvNFu0gDEQHNsx5vzBvgU+yA1KsUXBX\nYF65MB2WSnRqotpvsfSMUgacijJoxDrO6FMcp7HUs7mvi2on5Th67vPuBiXsMzQ1Nfk2mEENVCmd\n3A50GKLWZ2HPl0OlnM6knclidjhp0iSTzWZDOdNxnjfIKQiTn26ZIPnLLrvMNDQ0mHHjxhXtIMWp\nH4sRx0Hye8ZyZP3S1v2cznknf/3iiFMvlkuxZyonXd0EOZdRrgsjU6k0iqpTWNQpKN2hPwDIA3d7\nwi+wO/pnVOp6xynIZruKVlblVl5ROwHOSNHs2bPLHkkJc33Y5wtbWZQ7ClPqfNJOWjly7nz1pmMS\nHci4slHi8Yu3Eh25qI2r3zVROxZhO7LllpUoeePEUWoE1hu3MzWu1DSQsLpEfVMQ99kHotFOmrh2\nFURURzFJojpzUZ/Na89B8cctc1GfK8w1UQYDit2n1FuYMPfzc8K8026TdqijOCJJOC1J1wHuNCvH\njsup9+M+kzoFpZ2CH9id+qM84XVAO7C8Utc7TsGDDxb3PMs1aKdghy3czgheOp1OZCSl1PVhC/2S\nJeFeKybRoaxER6KYXnFGNZJu5MrRP+l44jSepSjW8QkacfPKJ5XmYZ2LsHVBlLxxnst5oxI2jZO0\ny3I6VAPVafXKVKrT7KdD0g591GsqkdZRnzHu4Ejc+5WKLwmcshG0XqMaZcL9vE49WGy6YCUo1v9x\n0szRxV1fRW2Tit0n6rXeRd/eejSMbuX0keK2x+oUlHYKHgJyQMbn3JPAu5W63nEKDjvsy30KY9Ij\nF+5XiGEqzdmzZw/om4JixOmERblvUOczqQ5wWL3c97vqqqsiNdyVmLIUVf9KxRN3FCZs59mdhkEO\ngLtMVsJZCaNvU1OTueqqq/rIu/WNMz8/6uhT3A6l332SHPkqhzADDW6ZSoyaGhOvvnGumTx5cq+e\npdKvmgMOpe7tLVvljlLHdQocPcMM4PhRrCz62ZsTX9Kj8o8//nho2WId70ri9+xBtud1WBz9ojhw\n7rijljmvfLH/w+pVTrmKe606BaWdgj8Cbwecu8t+C1BTiesdpwAmlNwVJ8xom7dgl6oMSxl5GB2i\nEGekqNKvtP0qmjC6Rb1PlI7U/PnzQ+W3g7vjGDf/ojpbA9WBC7IZbwMapiPqyAVNufI6AcU6jJVw\nGoOe39Fh/Pjx/e7v1d1vZC8or4qNPkUZnIjjoCU58lUOxZ7Xa2vuHT3CPlNcRyoMji3vuuuu/Ww5\nzH2i1AdR8iaKvbnPeTdpCBqljtq5D9u+ueMvtjapVDyl1soUa7eTbO/mz58feE9v+EA5AV786q2o\nukaps8ppw7zOnvf6sDsoVht1Cko7Ba8CbwSc+4XdqW+sxPWOUzBzZmkj9TOyUh51mMovSie9WEcr\nSiXtFJwktzeLykCNjPg9Q6mKeP78+f1GHYrla9iKvVh6Fjvnd/9ybSEuzn29W64FlQU/2/ZL2yAn\nICieStmLX/xO2Omnn95Hzm+dkJ+NRO3QRd35o1T85XSMo9hT3I51qTq32I5Lfuns7pxUunNQzJZL\npUUU3Uq9iQoqX0EyQbq4O1Vx6rFi94vqvBZ7jmKd96C0CrrOL+2idnD9zo0fP77ks1TaPkulVRLt\nhbfOSmpXMy+l0iqKDXupdLvpRp2C0k5B1d8UhCFKg1kpAwuKN0yH0tvp8m5BGuV+lXyWpOPwk3FX\nHn7p5n1TUKxzkpS+UTtHcWwhCYI6raU660GdJ7/4k+iQxKWYbbhH/oIot2NdCfm4nbiga8u1vXKc\nlGLnwjiVSePt/DmESYsouoXtEBUrX2HbiVJUKk29+oWpE6PUA841pRxod9sYdoFw3GdJMi2Lld84\ndWZcmyhnV7Ny2oZyyn9SZSMMg90pEEvH6iEiDwGfwvqmQM5z7kngY8aY3StxvYhU9+EVRVEURVGU\nYYMxRqqtQxA11VYAWAd8BjgCeMoJFJE64BBgdaWuH8wZoyiKoiiKoigDRaraCmBN8QG4xBN+LtAA\n/MoJEJF9RGRi3OsVRVEURVEURelP1acPAYjIT4B5wAPASqwPkl0IrDHGfMol9wYwzhiTjnO9oiiK\noiiKoij9GSxOgWCN9J8LTAA+AP4/MN8Ys90l9zqWU1AT53pFURRFURRFUfozGKYPOfv/3GCM2d8Y\n02CMGWeM+Ya3Q2+M2dvrEES5XiwuFZGXRaRDRP4mIteJSLbSz6iUj4gUAo5WH9n9ROQBEflQRNpF\n5AkROS4g3kYR+amIvGXbxQYRmRsgG8mGRGSWiDxl67BZRO4WkQnlpMNwRUS+ZaffX+x8f62E/FQR\neUREWkWkRURWicjkANk9ReQOEXlPRLaLyDoROS1AtlZErhKR10SkU0ReFZErRcR3jZaInCUiv7fj\nfUdEbhORXcrVWbGIYhcisjigDsmLyOd95AdFXkexTwVE5GN2vv3OTrNWEXlBRK7wq6uHYnsRRWfF\nIopdiMj8InXFv/rEPSjyOop9+lLt7Y8G8gBuBArAPcBXgeuAbuCRauumR6j8K2AtHP9nz/GPHrl9\ngM3A28BlwFzgeTuvP+mRzQBrgS7gWtsu7rXv9d1ybAj4PNaWuM/ZOnwTeAd4C9ij2uk51A473d8H\nfm3n72tFZI8EOoA/AxcBF9t/twKTPLI7Aa/Z5+YDXwMes+/3JZ+4H7DzdSFwNnCbLbvIR/ZS+9yj\ndrzfA9qA9UBDXJ31iG0Xi+28O9OnHmkejHkd1T71MAA/BlqAX2JNLT4XuNNOsxeAOpfskGsvouis\nR2y7mG/nx4U+dcXEwZjXUe3TN42qnUkDaAwH2Jlwtyf8AjvBzqi2jnqUzEPfxthH7m4gBxzkChsB\nvAG87JE93473fE/4vUAn1nS1yDaEtbPXRqzGvMEVPhnoARZUOz2H2gFMcP29nuKdv7XAVncFC+xl\nNwgPeWSvsfN1lissBTyL1dnMusJn2Xl9jSeO6+w4jnSF7Qy0A7/Dnqpph59kx3F5XJ31iG0Xi4F8\nyHgHRV5HsU89etPnE8Aon/Cr7bQ83xU25NqLKDrrEdsuHKfgIyHiHRR5HcU+A5+l2pk0gMbwAzvT\njvKE19mV+fJq66hHyTwsAIuwvOERATJZrBG4h33Ofdu2gcNcYU9ijebVemRn2Pf733FsCOvbGQXg\nCh89HgG2AOlqp+lQPSjS+QM+aqf9Qp9zP7cr3t1cYW8Cr/jI/oud36e5wv7DDhvrkW2273mTK+xr\ntuw/+8T9KrAhrs56RLcL+3yvUwCMwtWB95EdFHkdxT71KGkfB9ppf7P9/5BrL6LqrEd0u7DDHKdg\nvF1XBLbXgyWvo9hn0DEo1hQMEIdhJco6d6Axpgv4A3B4NZRSInMasB1oE5F3ReQnItLoOn8wVkF8\nxufaZwDBzmsREWAK8IIxptsjuxbLXtx2EcWGDsf6kmSQHo3AfsGPqZSBkw/FbOBQABHZAxhbQtZr\nAxuNMRvdgsaYt4BNPrLF9Pi4a75paJ2V8hGRFqxR+Q4ReVhEjvARq3pex7BPpTjj7N937N+hFgse\nrQAACRRJREFU2F6E1lkJjWMX73rCBfgjVl3Raa8BOMHn+qrndQz79GU4OQV7AR8Yz1ePbTYCuwQt\nHlMGDc9iee9fAM7Cmrt7AfCEq8Hdy/7d2P/y3rCx9u9OWN+y6CdrF6rNLlkn7rA2FEUPJVn2wqp0\ng9Je2JH2UfNprwBZR94rWyxucclE0VmJz9vADVjzcmcDP8TqgK8RkU96ZAdDXms9khAikgK+izUV\n4047eCi2F2oTCeKxi/90ndoK3IrVx/gccDnwEWCFiJzliWYw5HVU+/RlOHWCs1iLL/zodMn028lG\nGRwYY6Z5gv5DRNZjNewXYy0icpwDv7x25zMlZB15984BUWwoih5KsiRlA375VMoGvLLOaFHSeigx\nMMZc4QlaKiJ3Yo3m3QK4P445GPJa7SI5bgSOAL5ljPmzHTYU2wu1iWTxswuMMTd65JaLyCLgJeAG\nEbnX7NjhcjDkdVT79GU4vSnYjvUaxo96l4wytLgWaxX+ifb/Th765bU3n4vJOvJum4hiQ1H0UJIl\nKRvwy6dSNuCVRUQqoYeSEMaYV7EW8+0rIvu6Tg2GvFa7SAARuRprt5lbjTHXuE4NxfZCbSIhitiF\nL8aYLcACYDRwlOvUYMjrqPbpy3ByCjZhvcLJ+Jwbi/Xqp2eAdVLKxM6zTYCzF/gm+9fvNZkT5rxe\n24K1iKefrIjUYu0o4n4VF8WGouihJMsmgqfbjKXv1I2o+bQpQNaR98oWi9u4ZKLorCTPG/av+5sC\ngyGvtR4pExH5HnAl8P+MMed7Tg/F9kJtIgFK2EUx3rB/vXVFtfM6qn36MpycgnVYz9tnQZk9snMI\nngUiytDAzr9mdiwQWo/1+sw71Qg7zGDtDYyxluX/HpjiU5inYtmL2y5K2dBzHlkpokcr8Erxp1Ni\n4uRZMRt4HsAY8w5WRXlkgCz0z9exItKn4hWRZqz5n157CbKBqcB/u14/h9ZZqQjOwj73QsOq53UM\n+1RciMh8rPnitxtjzvERGYrtRWidFX9C2EUxguqKquZ1DPv0p9pbQQ3UgbXlVB64xxN+oR1+ZrV1\n1KNo/o0JCL/Wzr+vu8L89vUdCfyV4H1953nC77ML4/g4NsSOvYhfp+8+985exLdWO02H8kF53yn4\ntUfW2Qf+RFdYyo5jM67tb9mxd/21njiut+OY5grbBdhG/73rT7bj+FZcnfWIbhdY82nrfMKnYM23\nXe8JHxR5HcU+9eiTbt+1035xCbkh115E0VmP6HYBpIFGn/Bxdpl7l74fOhsUeR3FPgOfvdoZNMDG\n8BM7g+7D+tLb9Vjz0R+ttm56lMy7fweexlpUPAf4OtbuQwXgKU8B/SjwAda2c98EzsP6WmE38GlP\nvBks77kL66NEXwWW2HbyvXJsCGv71B4s7/08rN0LnJG/PaudpkPtwNqX/Uqs/ZnfsSvnK+3jXzyy\n07Bepb6KtQj9EvvvVuBAj+wYu4JuwfoK7TnA43Y+f9lHj6X2uduwvnL7c9sOb/eR/Vdb9jE73u9j\n7SO9wd0gRNVZj+h2gdXobgJuxvr68Ln23x12nkzzibvqeR3VPvUwYM0TL9jp9kXgf3mOT7tkh1x7\nEUVnPaLbBdAEfIj1XaRvYH2H5DqsKTrdwOcHY15HtU/fNKp2Jg2wQQhWY/AyVoX8JtZIs34RcpAf\nWFuCrbLzbDtWY/t7u5DU+shPBO63C3Y78FvguIC4G+0C/ZZtFxuA85KwIazRxqdtHTYDdwF7Vzs9\nh+LBjo6Q3/GYj/xU4DdYHa0WYCUwOSDuPYFfAO/Z9vUcAR+FAmqBq7C+SOl07K4g4OM2WNvnvmDH\n+w5WB3OXANnQOusRzS6A3e08/i+sUfourPnBi4D9BnNeR7FPPQzYH6kLW18wBNuLKDrrEc0u7HK/\nEHjRzocurA77XcChgzmvo9in3yF2JIqiKIqiKIqiDFOG00JjRVEURVEURVF8UKdAURRFURRFUYY5\n6hQoiqIoiqIoyjBHnQJFURRFURRFGeaoU6AoiqIoiqIowxx1ChRFURRFURRlmKNOgaIoiqIoiqIM\nc9QpUBRFURRFUZRhjjoFiqIowxgRmSkiBRE5q9q6KIqiKNVDnQJFUZRhgIhMFpH5IvIRn9P6aXtF\nUZRhjhijbYGiKMrfOyLyJWAxcKwx5gnPuVogZ7RBUBRFGbbUVFsBRVEUZUAQAt4IGGO6B1gXRVEU\nZZCh04cURVH+zhGR+cAi+9/V9hqCgogs8ltT4A4TkfNF5E8i0iEifxSRWbbMQSKySkRaROQDEblR\nRNI+995XRH4pIptEpEtEXheRa0QkOzBPryiKooRB3xQoiqL8/XMfsCdwDvAD4E92+F+AeoLXFFwA\njAZ+DnQCFwH3i8jpwG3AfwL3A/8AXAi8C/zIuVhEDgUeBbYAC4CNwGQ7nqNEZKYxJp/YUyqKoiix\n0TUFiqIowwB7TcEi4Dj3mgIRmQk8DnzZGHOHJ2wjsL8xpt0OPwh4ESgAXzDGPOiK5zlgT2PMWFfY\ni1iDT4cbY7a7wk/BciZ676koiqJUF50+pCiKogSx2HEIAIwx64FWYJPbIbB5EtjDmRYkIgcCBwF3\nAg0isrNzAE8D27DeMCiKoiiDAHUKFEVRlCBe9wnbUiQcYGf7d3/79/vA+57jXSAL7J6YpoqiKEpZ\n6JoCRVEUJYig+f7F1gGI5/d64KEA2S0B4YqiKMoAo06BoijK8GCgF5D92f7NG2MeG+B7K4qiKBHR\n6UOKoijDg3as0fsxA3EzY8wLwAZgrojs7T0vImkR2WkgdFEURVFKo28KFEVRhgfrsHYNulJExmAt\n9PVbG5AkX8TakvSPIrIIeAlrLcG+wOeBywHdfUhRFGUQoG8KFEVRhgHGmDeBrwANwM1Y3xiY65z2\nu6RYdCHv+SIwBfglcDLwE+BKYCrW9qiPholHURRFqTz6nQJFURRFURRFGebomwJFURRFURRFGeao\nU6AoiqIoiqIowxx1ChRFURRFURRlmKNOgaIoiqIoiqIMc9QpUBRFURRFUZRhjjoFiqIoiqIoijLM\nUadAURRFURRFUYY56hQoiqIoiqIoyjBHnQJFURRFURRFGeaoU6AoiqIoiqIow5z/AaplA5w5gGmt\nAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "timeseries = [sim.integ.dt*sim.nsteps_per_frame*i for i in range(nframes)]\n", "inst_KE, = plt.plot(timeseries[::nframes/1000], instantaneous_ke[::nframes/1000], 'ko', label='instantaneous KE',markersize=2)\n", "\n", "ke_1 = plt.plot(timeseries, cumulative_ke_1, 'r-', label='cumulative KE, x', linewidth=3)\n", "ke_2 = plt.plot(timeseries, cumulative_ke_2, 'b-', label='cumulative KE, y', linewidth=3)\n", "leg = plt.legend(prop={'size' : 12}, handler_map={inst_KE: HandlerLine2D(numpoints=1)})\n", "plt.xlabel('time');\n", "plt.ylabel('kinetic energy');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The cumulative kinetic energy for each degree of freedom should converge to the same value, which should be half the temperature (since we have $k_\\text{B}=1$). The instantaneous values of the kinetic energy should suggest a longer tail." ] }, { "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.11" } }, "nbformat": 4, "nbformat_minor": 0 }