VMC.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:cfc2534846ee4ffdfc8587dae54ac9502fecec2ad3f8ae11e01d2a95649af36e"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Variational Monte Carlo\n",
      "=======================\n",
      "\n",
      "For details of implementation/theory see\n",
      "[Martin, Reining and Ceperley: Interacting Electrons, ch. 21+22](http://adsabs.harvard.edu/abs/2016inel.book.....M)\n",
      "\n",
      "Author: Shivesh A. Pathak"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "import scipy.integrate as integ\n",
      "\n",
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Function proportional to wave function and related quantities\n",
      "-------------------------------------------------------------"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def psiT(x):\n",
      "    \"Trial wave-function\"\n",
      "    #return np.exp(-x**2/3.) #Trial wavefunction\n",
      "    return np.exp(-x**2/2.)  #Exact ground state"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def p2L(x):\n",
      "    \"local p^2\"\n",
      "    dx=1e-7\n",
      "    numer=(2.*psiT(x)-psiT(x+dx)-psiT(x-dx))/dx**2\n",
      "    return numer/psiT(x)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def x2L(x):\n",
      "    \"local x^2\"\n",
      "    return x**2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def EL(x):\n",
      "    \"local energy\"\n",
      "    return (p2L(x)+x2L(x))/2."
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def C(O):\n",
      "    \"Autocorrelation function and time\"\n",
      "    # TODO: should be vectorized at some point\n",
      "    Ck=np.zeros(N)\n",
      "    Obar=np.mean(O)\n",
      "    denom=np.sum((O-Obar)**2)\n",
      "\n",
      "    for k in range(N):\n",
      "        numer=0\n",
      "        for n in range(N-k):\n",
      "            numer+=((O[n]-Obar)*(O[n+k]-Obar))\n",
      "        Ck[k]=numer/denom\n",
      "\n",
      "    return Ck, 1+2.*np.sum(Ck)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Metropolis sampling\n",
      "-------------------"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def T(x,sig):\n",
      "    \"Transition step: symmetric gaussian\"\n",
      "    return np.random.normal(x,sig)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def A(x,y):\n",
      "    \"Metropolis accept/reject from x->y\"\n",
      "    a=psiT(y)**2/psiT(x)**2\n",
      "    u=np.random.uniform()\n",
      "    return a > u"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def Metro(N, x0, sig):\n",
      "    \"Generate sequence of samples using Metropolis algorithm\"\n",
      "    acc=0    #Number of accepted moves\n",
      "    x=np.zeros(N)\n",
      "    x[0]=x0\n",
      "\n",
      "    #get N samples, no binning\n",
      "    for n in range(1,N):\n",
      "        xprop=T(x[n-1],sig)         #Proposed move\n",
      "        if(A(x[n-1],xprop)):\n",
      "            x[n]=xprop\n",
      "            acc+=1\n",
      "        else:\n",
      "            x[n]=x[n-1]\n",
      "\n",
      "    print \"Acceptance rate = {0}\".format((acc+0.0)/N)    #Acceptance rate, ~ 0.50 is optimal\n",
      "    return x"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Run Metropolis\n",
      "---------------"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Number of samples you want in your average\n",
      "N=100\n",
      "\n",
      "#Standard deviation of transition probability guassian.\n",
      "#Too large of a sigma will reduce your acceptance rate, too small of a sigma will make\n",
      "#the steps unecessarily small. Your acceptance rate ~ 0.5, see the printout of Metro()\n",
      "sig=1.0\n",
      "\n",
      "#The initial sample point you want to start the sampling from.\n",
      "#If you choose x0 too far off 0, then you will have to discard some steps in the \n",
      "#averaging since they will be transient (see the plots of traces below.)\n",
      "x0=1.0    "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = Metro(N,x0,sig)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Acceptance rate = 0.65\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Plot histogram of samples, compare to psiT**2/(normalization)\n",
      "z=np.sort(x)\n",
      "norm,err=integ.quad(lambda x: psiT(x)**2,-np.inf,np.inf)\n",
      "plt.plot(z,psiT(z)**2/norm,'.') #Normalized trial wavefunction\n",
      "plt.hist(x,bins=20,normed=1);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEACAYAAABWLgY0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHsdJREFUeJzt3X2UVHed5/H3pwl0YgCfEtICCYkQIMDkRLIiZoZYGt0Q\nBsWNjoadNUQm6hmXk5mNOjijc9J95sxonLMZd2SimzHxYYwyMXETZJKIMVbCjCFiyBNPoY0GgQ5N\nHtQ0DiHQfPePewuKoh+qq6u7qm5/XufUqXtv/ar625fm07/+3Xt/VxGBmZllU1OtCzAzs6HjkDcz\nyzCHvJlZhjnkzcwyzCFvZpZhDnkzswwrK+QlLZS0XdIOSSt7eP0sSfdJelzS/ZImVr9UMzMbKPV3\nnrykJmAHcAnQAWwEroiI7UVtbgPWRMS3JOWA5RFx5ZBVbWZmZSmnJz8PaI+InRFxCFgNLClpMwv4\nMUBE5Ht43czMaqCckJ8E7Cpa351uK/YYcDmApMuBsZJeW5UKzcysYtU68PopICfpEWABsAfortJn\nm5lZhU4qo80e4Kyi9cnptqMi4lngfQCSTgXeFxEvlX6QJE+UY2ZWgYhQJe8rpye/EZgmaYqkMcAV\nwJriBpJeL6lQwF8Ct/RRaN0/rrvuuprX4Dpdo+t0nYXHYPQb8hHRDawA1gFbgNURsU1Sm6TFabMc\n8JSk7cAE4G8HVZWZmVVFOcM1RMS9wIySbdcVLd8B3FHd0szMbLB8xWsPcrlcrUsoi+usnkaoEVxn\ntTVKnYPR78VQVf1iUgzn1zMzywJJxBAeeDUzswblkDczyzCHvJlZhjnkzcwyzCFvZpZhDnkzswxz\nyJuZZZhD3swswxzylhktk1uQVNGjZXJLrcs3GxK+4tUyQxK0VvjmVgY925/ZUPEVr2Zm1iOHvJlZ\nhjnkzcwyzCFvZpZhZYW8pIWStkvaIWllD6+fKel+SZskPSbpsuqXamZmA9VvyEtqAlYBlwKzgaWS\nZpY0+yzwrxExF1gK3FjtQs3MbODK6cnPA9ojYmdEHAJWA0tK2hwBxqfLrwH2VK9EMzOrVDn3eJ0E\n7Cpa300S/MXagHWSrgFeBbyzOuWZmdlglHUj7zIsBb4WEf8gaT7wLZKhnRO0trYeXc7lciPiHotm\nZgORz+fJ5/NV+ax+r3hNQ7s1Iham658GIiKuL2qzGbg0Ivak608Db4mI50s+y1e82pDxFa+WVUN9\nxetGYJqkKZLGAFcAa0ra7CQdopF0HtBcGvBmZjb8+g35iOgGVgDrgC3A6ojYJqlN0uK02SeBj0h6\nDLgVWDZUBZuZWfk8QZllRq2Ga1omt9C5p7Oi954x6Qz27t5b0Xtt5BjMcE21DryajVidezor/uXS\n2VrZLwezcnlaAzOzDHPIm5llmEPezCzDHPJmZhnmkDczyzCHvJlZhjnkzcwyzCFvZpZhDnkzswxz\nyJuZZZhD3swswxzyZmYZ5pA3M8swh7yZWYY55M3MMqyskJe0UNJ2STskrezh9RskPSppk6SnJL1Y\n/VLNzGyg+r1piKQmYBVwCdABbJR0V0RsL7SJiGuL2q8ALhiCWs3MbIDK6cnPA9ojYmdEHAJWA0v6\naL8U+E41ijMzs8EpJ+QnAbuK1nen204g6SzgbOD+QVdmZmaDVu17vF4B3N7X3bpbW1uPLudyOXK5\nXJVLMDNrbPl8nnw+X5XPUn93qJc0H2iNiIXp+qeBiIjre2i7Cfh4RGzo5bP6yn+zQZFU8Q21aYVK\nfzZr9XVt5JBERKiS95YzXLMRmCZpiqQxJL31NT0UMRN4TW8Bb2Zmw6/fkI+IbmAFsA7YAqyOiG2S\n2iQtLmr6QZKDsmZmVifKGpOPiHuBGSXbritZb6tiXWZmVgW+4tXMLMMc8mZmGeaQNzPLMIe8mVmG\nOeTNzDLMIW9mlmEOeTOzDHPIm5llmEPezCzDHPJmZhnmkDczyzCHvJlZhjnkzcwyzCFvZpZhDnkz\nswxzyJuZZVhZIS9poaTtknZIWtlLmw9I2iLpSUnfqm6ZNtxaJrcgqaJHy+SWWpc/cKOo+Ps1q2f9\n3hlKUhOwCrgE6AA2SrorIrYXtZkGrATeGhEvSTptqAq24dG5p7Pim1N3tnZWtZZh0c2gbsZtVq/K\n6cnPA9ojYmdEHCK5j+uSkjYfAf4pIl4CiIjnq1ummZlVopyQnwTsKlrfnW4rNh2YIenfJf1E0qXV\nKtDMzCpX1o28y/ycacDFwFnAg5LmFHr2xVpbW48u53I5crlclUowM8uGfD5PPp+vymeVE/J7SIK7\nYHK6rdhuYENEHAGekbQDOBd4pPTDikPezMxOVNoBbmtrq/izyhmu2QhMkzRF0hjgCmBNSZs7gbcD\npAddzwV+UXFVZmZWFf2GfER0AyuAdcAWYHVEbJPUJmlx2uYHwAuStgA/Aj4ZEb8ewrrNzKwMZY3J\nR8S9wIySbdeVrH8C+ET1SjMzs8HyFa9mZhnmkDczy7BqnUJp1rgOjgXmwEvPQMeF8PxMOG07nPMA\nNO+vdXVmg+KQt+x7qQWe+O+AYOadcOB0mLA5CfCDY+GW9cAs+OJhOHLKsfedvhWunp+02bEYpq+F\n8Xtr9V2YVcQhb9nUcT489Ak47zb47u0Qzcn2+66Hpu4kwJcvgH1z4LlZwBg4MhoomnDs+enwy4vh\n9u/C4VPgpAOw/K3QMc+Bbw3DIW/Z03E+3PRYsvzkh9KNhfBugiOj4LnzYN/spEd/+lboPC8J/+Ke\n/Gk74DdvTAIeJc9f3Zj8MjjpAFwz1UFvdc8hbxkyFja9F35cuDpQQBQ9AxyBpsNw+jaYsCUZslm+\nAD43G/58J3TMhRdmwOufgnMeTIZq7rs+CfimV+DIGI4GfvsimPIA8Kfw5AaYfq/H8K3uOOQtE7q6\nADbCmuLLOQrB/gqcug/m/RPM/i4cOO1YwEP6/DCMB8bfDdx97COa9yc99vZF8IafwS0PHRu6ef1W\nWPUU0AR3AKN/C398GZy9Yci/X7NyOeQtEzZvhmSOvMKwTMBrfg7vvRJGxfGhztMD+/Dxe+HCW5Ll\nQuCfezc8fA3JWcjp1zz0avj6T+CjF8DEJwb5HZlVh8+Tt0yYMwfg5yS997QH/4H3J73qMx+u3jBK\nIfDH74U33QwcKfqaadj/xydhV3pWjlmNOeStoXR1wX33JY9kiCYxbhzAm+E9H4LZtw5Pb/q0p2HF\nDODr6Yb0l8uz8+BrDySnZjrorcYc8tYwurrgoovgXe9KHhdddHzQw36Yeyv80YeGb7jktKeB5ckv\nlfO/DkuuhN+ckxygLZzBY1ZDDnlrCF1d8J3vwLZtx7Zt3w5bttSupuNMfAIuXw6z7kxOyWw6eOwM\nHrMa8oFXq3tdXbBgQXJw9aSToLs72T5zJsyut45y4ZTMfbOPP9h7cGxy4VXhSluzYeKQt7q3eXPS\nY+/uBgm+9CU47zyYN68wFl9nmvcnB3sLClMnPDfr2JW2DnobJmUN10haKGm7pB2SVvbw+jJJ+yRt\nSh/Lq1+qjVRz5iQ99tGjk+dly+CSS+o04HtSmDrB4/RWA/2GvKQmYBVwKTAbWCppZg9NV0fE3PRx\nS5XrtBFs3DhYvx4efDB5bphwLyhMnVAYpx91AH74d/D81FpXZiNAOcM184D2iNgJIGk1sATYXtJO\npW80q0RXVzJEM2fOsUAfNw7mz69tXRUrHqcfdQBu2gQ0wX/8BfAHta7OMq6c4ZpJwK6i9d3ptlKX\nS3pM0m2SJlelOhtxCgdZL744eT7+FMkGVhinf/xKjl0l2wQ8SEdHbUuzbKvWKZRrgLMj4gLgPuAb\nVfpcG2EKB1kPH4atW+voFMlqmVw8r42Ak7j77t4amw1eOcM1e4CzitYnp9uOiohfF61+FfhCbx/W\n2tp6dDmXy5HL5coowbKuMEQzZUpycHXrVpg1qw5PkRys6ffCa9vh1+emG15m0aJT+nyLjTz5fJ58\nPl+Vz1JE9N1AGgU8BVwCPAv8FFgaEduK2rRExN50+b8Bn4qIi3r4rOjv61l9kAStFb65FQby77xj\nR3K2zN69SajffTf86lfJ8kAOsg625mF778Gxyc1Ifns23PM9Ip6t8AvbSCGJiKjouGe/PfmI6Ja0\nAlhHMrxzc0Rsk9QGbIyItcA1kt4DHAJeBK6qpBgbeTo64Pd+D155JVnfujUJ+IY9yFqO5v0wMx2j\nuae2pVj2lXUxVETcC8wo2XZd0fJfAX9V3dJsJLj99mMBD3DGGRkcoinTpk1www1w7bUwd26tq7Gs\n8Nw1VjNdXfCVrxxbHz0afvSjBjwPvgo2bYILL4Rbb02eN22qdUWWFQ55q4nChGPt7cl6UxPccQdM\nn17bumrlhhuOX//MZzJ0+qjVlEPehl3hXPiPfzyZcGz06GRcfiSfaHXttcevr1uXzM+zY0dt6rHs\ncMjbsCuecOzwYbjxxgadrqCK5s6FRx6BhQuTv2qOHIE9e5Jffr5YygbDIW/DrnTCsQ9+cGQHfMHc\nuXDbbfCGNxzb9sor+GIpGxSHvA27hp9wbAiNGwf33w9jxiTrJ58MixbVtiZrbJ5P3mqioSccG2LT\np8Mvf5n04BctgokTa12RNTKHvFkdmjgRrr661lVYFoyokL/641dz+3dvr+i948aN45ENjzBhwoQq\nV5VtHR2wdi0sXuweqVktjKiQf3jjw/z2kt/2PFFyP7q/082+ffsc8gPQ0QFTp8LLLydjy08/7aA3\nG24j78Drq4DxA380jRp5u2qw1q5NAh6SZ58lYjb8nFw2ZBYvTnrw4LNEzGrFIW9DZuLEZIjmn//Z\nQzVmtTKixuRt+PksEbPack/ezCzDHPJmDayjA266yfPbWO/KCnlJCyVtl7RD0so+2r1P0hFJvuXB\nSNNxPtzxjeTZhkXhFNWPfSx5dtBbT/oNeUlNwCrgUmA2sFTSzB7ajQWuATaUvmYZ13E+3PQYPPmh\n5BkH/XAoPUX1e9+rbT1Wn8rpyc8D2iNiZ0QcAlYDS3po9zfA54GDVazPGsGDn00XCvcZ/vNaVTKi\nLF4Mzc3H1r/8Zd9oxE5UTshPAnYVre+m5JpRSW8CJkeEb0s80rzUAtv/MF2J9PmLtapmRJk4MZma\nuCn9X9zenszTb1Zs0AdeJQm4AfhE8ebBfq41iB2LgVM4+k++oBV4onb1jDBvf3tyY5HRo2HWrJF7\nE3TrXTnnye8Bzipan5xuKxhHMlafTwO/BbhL0nsi4oTbEbe2th5dzuVy5EbyPd+yYPpaOOkAHD4l\neX7zTbC+1kWNHIW5+bdsSQLec/NnQz6fJ5/PV+Wzygn5jcA0SVOAZ4ErgKWFFyPiJeDorF2Sfgxc\nGxGP9vRhxSFvGTB+L1wzFdoXwbl3J+s2rDw3f/aUdoDb2toq/qx+h2siohtYAawDtgCrI2KbpDZJ\ni3t6Cx6uGVnG74ULbzkW8KNAUkWPlskttf1ehpv3lQ2xsqY1iIh7gRkl267rpe07qlCXNbJuoLWy\nt3a2dlazkvrnfWVDzFe8mpllmEPezCzDHPJmZhnmkLfyeG4as4bkkLf+lc5N46A3axgOeevfQ4WL\nmdMzYzd4bhqzRuGQt/699X+nC+ncNPM9N41Zo3DIW/8mPgEfvQDO/3ryPNFz05g1Ct/j1coz8Qm4\nfHmtqzCzAXJP3noxFnbNh4Nja12ImQ2CQ95OkNx4Yj187QG4Zb2D3qyBOeQzrGVyS0UTX40f/1Zg\nFhwZA8+dB/s8SbkdU+nPlSdVqw2PyWdY557Oyia/OrgZPrcVms6D07fBBN9uyI6p+OcKT6pWCw55\nO+bgWNg3ByZsBhbAh2cnAd+8v9aVmVmFHPKWODg2GX9/bhacvhVYAGc+XOuqzGyQPCZviX1zkoAv\njMPjcXizLCgr5CUtlLRd0g5JK3t4/WOSnpD0qKQHJc2sfqk2pCZsTnrwTQeTcXg8Dm+WBf0O10hq\nAlYBlwAdwEZJd0XE9qJmt0bE/03bvxv4B+CyIajXhkrzfli+IDmTZsIW+JzH4c2yoJye/DygPSJ2\nRsQhYDWwpLhBRBQnwljgSPVKtGHTvD8Zh/eBVrPMKOfA6yRgV9H6bpLgP46kjwPXAqMB3+fVzKwO\nVO3smoi4EbhR0hXAXwNX9dSutbX16HIulyOXy1WrBDOzTMjn8+Tz+ap8Vjkhvwc4q2h9crqtN/8K\nfKW3F4tD3szMTlTaAW5ra6v4s8oZk98ITJM0RdIY4ApgTXEDSdOKVhcDOyquyMzMqqbfnnxEdEta\nAawj+aVwc0Rsk9QGbIyItcAKSe8EXgF+DSwbyqJtEIqvavUBVrPMK2tMPiLuBWaUbLuuaNn3g2sE\npVe1Ll/goDfLOF/xOpKUXtXq2SXNMs8hP5KUXtXq2SXNMs8TlI0kpVe1eqjGLPPck8+6l1rgZ1cn\nz+CrWs1GGPfkM+18+OJGODIaTjoA10yF8XtrXZSZDSP35DOqowPgp0nAIzh8CrQvqnFVZjbcHPIZ\ntXYtwBhAQEDTITj37toWZRkxFnbN9w3eG4RDPqMWLwY4QBLwr8DVb/ZQjQ1aVxfAevjaA8k1Fw76\nuucx+YyaOBFgKrx7UdKDd8BbFWzeDFByrYVvE1nX3JPPtL1w4S0OeKuaOXMA0mstXv8UvHKqe/N1\nziFvZmUbNw5gAfzxZcnhnlvv8bBNnXPIm9kA7YcxB+D5mZ4iowE45M1s4EqnyHj1Tp9xU6d84NXM\nBq54ioxX70yGbTy7aV1yT74BtExuQdKAHw1pFBV9rw37/TaywhQZvz3bs5vWMffkG0Dnnk5o7eGF\n/m4A0tN76l03lddd6ftscApDN8+d59lN61BZPXlJCyVtl7RD0soeXv9fkrZIekzSDyWdWf1S7TiF\nG4D4ohSrtcLQzYffljyDx+frSL8hL6kJWAVcCswGlkqaWdJsE3BhRFwA3AH8fbULtRK+AYjVk8LQ\nDSSdjpsfgFVb4Pmpta3LyurJzwPaI2JnRBwCVgNLihtExAMR8XK6ugGYVN0y7QS+AYjVo31zoHMW\nMAa6zoQvbz42zbXVRDkhPwnYVbS+m75D/E+AewZTlJWh9E9kn81g9WDCZhi3FwhA0N3s2U9rrKoH\nXiX9D+BC4G29tWltbT26nMvlyOVy1Swhu3o6yFr8J7JZPWjeD8vekfTgu5uT+xh49tMBy+fz5PP5\nqnxWOSG/BziraH1yuu04kt4J/CVwcTqs06PikLdypQdZfR6yNYLTnoY/OyfpwXtyvIqUdoDb2toq\n/qxyhms2AtMkTZE0BrgCWFPcQNKbgK8A74mIFyquxnrhg6zWYMZ7crx60W/IR0Q3sAJYB2wBVkfE\nNkltkhanzb4AnAp8V9Kjku4csopHmOQOTxfC69t9kNXMBqysMfmIuBeYUbLtuqLld1W5LiMJ+KlT\nAb4ELx6AD1wO5zzooRozK5unNahja9fCyy9DcpbCKfC7Fge8mQ2IQ76OLV4MJ58MED5Lwcwq4pCv\nYxMnwtNPA1wN10z1QSwzGzCHfJ1L7tXqsxSsB4OYsbNlsq9CHSk8C2Ud6epKbpQ8Z07hNmtmfRjE\njJ2drZ3VrMTqmHvydaKrCxYsgIsvTp67umpdkZllgUO+TmzeDFu2wOHDsHVrsmxmNlgO+ToxZw7M\nng2jR8OsWcmymdlgeUy+TowbB+vXJz342bM9Jm9m1eGQryPjxsH8+bWuwsyyxMM1ZmYZ5pA3M8sw\nh7yZWYY55M3MMswhb2aWYQ55M7MMKyvkJS2UtF3SDkkre3h9gaRHJB2SdHn1yzQzs0r0G/KSmoBV\nwKXAbGCppJklzXYCy4Bbq16hmVkDapncUhezhJZzMdQ8oD0idgJIWg0sAbYXGkTEr9LXomqVmZk1\nsM49nXUxS2g5wzWTgF1F67vTbWZmVueGfVqD1tbWo8u5XI5cLjfcJZiZ1b3irByMckJ+D3BW0frk\ndFtFqlW4mVmWFWdlW1tbxZ9TznDNRmCapCmSxgBXAGv6aK+KqzEzs6rqN+QjohtYAawDtgCrI2Kb\npDZJiwEk/RdJu4D3A1+R9ORQFm1mZuUpa0w+Iu4FZpRsu65o+WfAmdUtrf787ndNPPSQ78FqZo3D\nV7yWKY6M5corz/E9WM2soTjky9R96Dx+8Ytm34PVzBqKQ743L7XAz65OnoFRo7fxxjce9D1Yzayh\n+PZ/Pdn6h3DbXUATnHQArpmKmvbzzW/+kohZvgermTUM9+RLPTMfbvs+ya4RHD4F2hcBcOqpR5g/\n3wFvZo3DIV/qwb9OFwQEcATOvbuGBZmZVc4hX+riv0kX0rnWPrAExu8FIPfOXMWzyo1qHlXxe83M\nKjWixuS7u18F++ZDy2Zo3t9zo7M3wFUXwfrPwIK/TdZTL3S+UPGsckdaj1T83orfZ2Yj3ogJ+a4u\n+MX2b8Djb4SfbYXlC/oO+rPfPbwFmpkNgREzXLN5Mxw8MBWOjIHnzoN9PgfSzLJvxIT8nDnQfMrT\n0HQQTt8GE3w1k5ll34gZrhk3Dt44cxnbJgrO39L7UI2ZWYaMmJAHGDXqP2HCZmiudSVmZsNjxAzX\nmJnV1MGxsGt+8jyMRlRP3sysJg6OhVvWw3Oz4PR+zu6rsrJ68pIWStouaYeklT28PkbSakntkh6S\ndFZPn2NmNiLtm5MEfDln9x0cC8yv2nTm/Ya8pCZgFXApMBtYKmlmSbM/AV6MiHOBLwJfqE55NfLL\nWhdQJtc58jTKvmyQOvP5/PB8oQmbkx58f2f3FXr8PFC1+1aU05OfB7RHxM6IOASsBpaUtFkCfCNd\nvh24ZPCl1dAztS6gTM/UuoAyPVPrAjLkmVoXUKZnal1AeYYt5Jv3J0M0H35b30M1hR4/Y6p234py\nQn4SsKtofXe6rcc26T1hfyPpdYMvz8wsI5r3w5kP9z0WX+jxc7Bq960YqgOvdTmrVvOYZk594FRG\nnTqqz3Yvd77MyS+efNy2Ay8cGMrSzMyO9fg/N5v16zdUZVpzRUTfDaT5QGtELEzXPw1ERFxf1Oae\ntM3DkkYBz0bEhB4+q+8vZmZmPYqIijrP5fTkNwLTJE0BngWuAJaWtPk+sAx4GPgj4P5qFmlmZpXp\nN+QjolvSCmAdyRj+zRGxTVIbsDEi1gI3A/8iqR14geQXgZmZ1Vi/wzVmZta4hnRaA0lfkLRN0mOS\n7pA0vpd2fV5sNdQkvV/SZkndkub20e4ZSY9LelTST4ezxvTrl1tnrffnayWtk/SUpB9IenUv7bol\nbUr3553DVFtDXNhXRp3LJO1L998mSctrUOPNkjolPdFHm39M9+Vjki4YzvqKauizTklvk/Sbon35\n2eGuMa1jsqT7JW2R9KSka3ppN7B9GhFD9gDeCTSly58HPtdDmybg58AUYDTwGDBzKOvqoYYZwLkk\nxxLm9tHuF8Brh7O2gdZZJ/vzeuAv0uWVwOd7affSMNfV774B/hS4MV3+ILC6Bv/O5dS5DPjH4a6t\npIY/AC4Anujl9cuAf0uX3wJsqNM63wasqeW+TOtoAS5Il8cCT/Xw7z7gfTqkPfmIuC8ijqSrG4DJ\nPTQr52KrIRURT0VEO/2f+ilqOKlbmXXWfH9y/MVx3wDe20u74T4Q3ygX9pX7b1jTExki4t+BX/fR\nZAnwzbTtw8CrJZ0xHLUVK6NOqIPTviNib0Q8li7vB7Zx4jVJA96nwxlYy4F7ethezsVW9SKAH0ja\nKOkjtS6mF/WwPydERCckP7jACafTppol/VTSTyQNxy+iRrmwr9x/w8vTP9lvk9RTB6rWSr+PPdTv\n/+356bDhv0maVetiJJ1N8tfHwyUvDXifDvpiKEk/BIp/k4gkDD8TEd9P23wGOBQR3x7s16tUOXWW\n4fcj4llJpwM/lLQt7SXUW51Dro86exrP7O3o/pR0f54D3C/piYiot1lPat7D68Ua4NsRcUjSR0n+\n+mjs6URq5xGSn8X/lHQZcCcwvVbFSBpL8lfkn6U9+kEZdMhHxLv6el3SVcAi4B29NNkDFB/cmpxu\nq6r+6izzM55Nn5+T9P9I/qyuashXoc6a78/0INcZEdEpqQXY18tnFPbnLyXlgTcxtFNblbNvdgNn\nAh3phX3jI+LFIaypJ/3WGRHFww9fpT4nBdxDsi8LhuRncbCKgzQi7pF0o6TX1eDfHUknkQT8v0TE\nXT00GfA+HeqzaxYCnwLeExEHe2l29GIrSWNIzrFfM5R19aPHnpukV6W/YZF0KvBfgc3DWVhpSb1s\nr4f9uQa4Kl1eBpzwwyrpNWl9SDoNuAjYOsR1lbNvChf2QR8X9g2xfutMf3kWLGHo911vRO8/i2uA\nK+HolfO/KQzj1UCvdRaPaUuaR3Jq+bAHfOoWYGtE/J9eXh/4Ph3io8XtwE5gU/oonLXwBmBtUbuF\nJEeS24FP1+Co9ntJxrkOkFzVe09pncA5JGc5PAo8Wa911sn+fB1wX1rDOuA16fYLgZvS5bcCT6T7\n83HgqmGq7YR9A7QBi9PlZuC29PUNwNnDvf/KrPPvSDoZjwI/AqbXoMZvAx3AQeBXwIeBjwEfLWqz\niuRMocfp48y1WtYJ/M+iffkT4C01qvP3ge6inNmU/hwMap/6YigzswzzPV7NzDLMIW9mlmEOeTOz\nDHPIm5llmEPezCzDHPJmZhnmkDczyzCHvJlZhv1/1J8vqL++tt0AAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f61c80b5410>"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Postprocessing\n",
      "--------------"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Get local quantities\n",
      "p2 = p2L(x)\n",
      "x2 = x2L(x)\n",
      "E = EL(x)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Plot trace of local p^2/2., x^2/2. and E\n",
      "n=np.arange(N)\n",
      "plt.subplot(221)\n",
      "plt.plot(n,p2/2.,'b')\n",
      "plt.subplot(222)\n",
      "plt.plot(n,x2/2.,'r')\n",
      "plt.subplot(223)\n",
      "plt.plot(n,E,'k')\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEACAYAAABRQBpkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztvXmYVNW1v/+uBpoGZOpmBhkVFFAZFDWo4IwmwWii0dwk\nasw3g2OiRhOMF7n6y3XIcE00N7mKU26U6xAjMSI40A5RERlEkaEZhW7mnhhkXr8/dp2uU9VV1TVX\ndfV6n6eernPOrrP3rj61P3uttQdRVQzDMIyWTVGuC2AYhmHkHhMDwzAMw8TAMAzDMDEwDMMwMDEw\nDMMwMDEwDMMwSJMYiMh0EdkiIktipPm9iFSIyGIRGZWOfA0jk4hIPxF5U0SWisgnInJjlHT2bBvN\nnnRZBo8D50e7KCIXAENU9Wjgh8Cf0pSvYWSSg8DNqjoCOBW4TkSO8SewZ9soFNIiBqr6LlATI8lF\nwFOBtPOAziLSMx15G0amUNXNqro48H4XsAzoG5bMnm2jIMhWzKAvsMF3XEnjH5Vh5C0iMhAYBcwL\nu2TPtlEQWADZMJpARI4AngduClgIhlFwtM5SPpXAkb7jfoFzjRARWyzJyCiqKvGmFZHWOCH4i6q+\nFCFJXM+2PddGpknkuY5EOi0DCbwiMRP4LoCInALUquqWaDdS1ay/pk6d2qLybal1ToLHgM9U9cFU\nn237Pxd+3rnKNx2kxTIQkaeBiUCZiHwOTAWKAVXV/1HVV0TkQhFZBewGrk5HvoaRSURkPPBvwCci\nsghQYAowAHu2jQIjLWKgqt+KI8316cjLMLKFqv4LaBVHOnu2s0FdHfzqV3DffbkuSUFiAeQAEydO\nbFH55jLvXNa5JVIw/+dPP4Wnn85N3nHSnJ9tSZe/KV2IiOZbmYzCQUTQFANtSeZrz3Wq/O//wq23\nwubNuS5J3pGO59osA8Mwmgfr1sGBA7kuRcFiYmAYRvPAxCCjmBgYhtE8WLvWxCCDmBgYhtE8MMsg\no5gYGIaR/xw6BBs3ur8WiM8IJgaGYeQ/lZXQvTu0aWPWQYYwMTAMI/9ZuxYGDjQxyCAmBoZh5D/r\n1sGgQU4M9u/PdWkKEhMDwzDyn7Vrg2JglkFGSNceyJNEZLmIrBSR2yNcv1JEtorIwsDre+nI1zCM\nFsL69TBggIlBBkl5oToRKQIeAs4GqoD5IvKSqi4PSzpDVSNuKG4YhhGT3buhUycTgwySDstgHFCh\nqutV9QAwA7cvbDhZXw/GMIwC4cABaN3axCCDpEMMwveA3UjkPWAvEZHFIvKsiPRLQ76GYbQUDhxw\nQmBikDGyte3lTOBpVT0gIj8AnsS5lfKWbdvcsGY/n3wCf/kL1NTAWWfBFVc0/swf/+jODx2aWv6f\nfAK1tW4J923b4MQTYeRIkCTsq0OHYPZs2Ls39PxZZ0GXLjB3LrRvDyefDC++CD16wPjxqZU/GocP\nw5o1cNRRwXP797u6dujgXvHyn/8JX3wBp50G553n3ldUhKbZvx/eeAP27IFp09JTByMHHDxolkGG\nSYdlUAn09x032gNWVWsCLiSAR4GxsW541113NbzKy8sD94A5c+D0011D+9Ofplboe+6B0lI44wzX\n8Pp5/nno08c1wh6zZrnGs7jYfe5Pf2p8z5kz4a9/dY3q++/Hzn/37ujPdF0dnHIKTJkCDz/sGrOv\nfQ2OPBKuvRaqqmDLFrjoIvcX3Pn16yPf78kn4cYb3QrA3utnP4Nnnglev/hi+OAD+Ld/c3XNFHPn\nusbbP4n0Rz+Co4928cGNG+O7z/79cNddUFQEN9zg7tm/P1x+OXz728HX5ZeX89xzd7F+vXuejGaK\nWQaZJw17b7YCVuG2AiwGFgPHhqXp5Xt/MfBejPvpPfdoI55/XrV/f9X/+z/VF15QHTu2cZpEOO88\n1aeeUn30UdVu3VSvuUb1hz90f3v0UO3eXfXTT4Ppf/c71RtvdO9XrlQdPLjxPb/7XdU//Un12992\n9w6nslL1/fddnt27q37nO5HLNnOm6tlnh547fFh1+XLV229X7dXL5d++vepbb7nrAweqvvtu43sd\nPKg6bJjq3Lmh5++4Q/Wuu9z7Cy9UPeEE1datVY8/XvXaayOXKx38+7+rgurSpcFzxx2n+tFHqnfe\nqXrJJfHdp6LC1VlVdd8+1RkzVFevbvpz7pGP+9meDmwBlkS5PgGoBRYGXr+Mca/4KmZE5vTT3cN+\n8smq772X69LkHYk819FeKbuJVPWQiFwPzMFZGtNVdZmITAPmq+rLwI0iMhk4AFQDV8W65x//CL/4\nhev1eTz3HNx5J1x2mRty7O+1J8OSJTBxouttn3oqvP128NpPfgI//jHs2BE8t2MHlJW59337utnx\nqqFum7ffhp//HD77LPjZV16B8nJX5jffdD3grl3h2Wfh//0/55a5+OLQss2d66wQPyIwbBjcey98\n5StuDs6LLzrLQNX9jdRheukl5wqaMCH0fPfusGqVe79tGzz4ICxb5tw04ZZBdbWzhtLBO+84C+CN\nN2D4cNfDr6hw70eMgOOOc1ZQu3bBzxQVQc+ezhU0aZI7t2YNDBni3hcXwze/mZ7yhfE48AfgqRhp\n3lbVyRnJ3QjiuYmKi80yyBDp2gP5VWBY2LmpvvdTcBuJx0VVFXz0EYwb547373c+7wcfdMfdu6cm\nBlu3Ov95v0AYe/hw9/JTVgbbtwePt293DRU4/3qHDu6cF1f4/HPn+jnmGPfZ6mp3/qGHoHdv15BN\nn+5Gx3k8+aQ7X1LixO+WW9z5uXOdIEbjtNPc6733nAjs2uX85ZEmZj7+uHMRhccauncPurK2b3ff\nxYQJ8M9/Ov+9h6qb67N+vROVVNi/Hz78EB54wLn8brgBVqxwqwx4jf+cOc5d5efgQbe51U03ue/3\nb3+D1ath8ODUytMUqvquiAxoIpmNkssG5ibKONkKICfEwIGu1+uJwdtvu0agZ0933KGDa6R27w4G\nHKdMgb//HcaMcT7xWHzyCRx/fOxgbLdu0S0DcI3nxo1BMXjnHRd/EHG96GXL3Plt22DqVBdHCOdL\nX4JNm5yQnHGGs3rat3e93hNPjF0HcN/H1q3BuEEkMVizJihifvyCum2bqy84q6WmJpiuvt691q6F\n0aObLlMsFixw8Z6vfQ3uuMMFtpcscf8Lj0GD3CsSN90Exx7rvlu/ZZBjThGRRbg5Nj9T1c9yXaCC\nxALIGScvxeD733ejdqYEbIm//Q0m+wxxEdd4bd8eFIP333c965/+tGm3xscfhzZAkSgrCxWD7duD\nDSYExcBrIN9+2zXo4PL2LIPwz4XTsaNzj1x3nespjxvnRKJNm9jlAzfqZ/Hi6GKg6oSmf//Gn/XE\nYO9e9znPYgkXA886SocYvPOOGwDQuzf06gWLFjUWg1gUF7vA+ocfOsvA6yzkkAXAAFXdIyIXAH8H\noo4j8wewJ06c2Kw3T886fsvA1iaivLy8YXBNushLMbjgAjcyp08fd9y+Pbz7bmgarzEbEDDid+92\nQy9PP9355r/xjej3X7LEuVliUVYW6oqKZhl4fPIJfOc7wc96YhBpiGokbr8drr7alf3HP246PTjL\nYMsWZx1A4w5Tba3zt3fu3Piz3vfniZVnJUUTg3Xr4itTJJYtc0NVd+1yI7XA/X/++EdnGcVbX3AC\nMH++swwy7SZqClXd5Xs/S0T+KCKlqlodKb2NZkoBcxOFEN6ZmJaGcdN5KQY9e8K8ebHThMcNdu92\nonHuuc7v3JQYXHtt7PuXlcFy34IaTYnB9u2upw7OMtixw/nxDxxwvf+m6NDBBZUTwRODaJZBNKsA\ngpbV1q2hlosnBl5w3LOO1q5NrGx+Kiudq+of/whaILfc4lxGe/bAf/93/Pc66SRnNa5enTU3kRAl\nLiAiPVV1S+D9OECiCYGRIuYmyjh5KQax3Cr+NH4x2LPHNajnnusCzf6RPvv2Ob+9N+lq2TLnmmnq\n/l6vWNW9DxcDv5Xmdwd5biLPKkhmolg8pCIGbds68Vy1KtRyadvW/ea873P7dmdZpCIGe/Y4EfAH\nzzt3dm7AO+8MWnfxMGqUs8I6dEg9oN0UIvI0MBEoE5HPgam44dOqqv8DfENEfowbJfcFkJkxTYZZ\nBlkgL8Wgbdum03TvHjraxwsmDxrkGsXVq4OzXBcscENTb7jBHT/2WNMzXf0xgz17nLulffvgdb9l\ncPAg7NwZbJy8zzYVL0gVfwBZJDExAPcdfvZZYzdWly7OOvDEYOzY1MTAH+j38+MfO7dPImLZvr0L\nIscTU0kVVf1WE9cfBh7OfEkME4PMk5diEA+R3EQdOriG5Zxz4PXXQ8XgnHPc/IF48YtBpEbdLwbV\n1c694s2L6NTJlWfTpvjiBclyxBFueYe1a52LKvw3Eo8YLFsWHKXl0bWrizf06+fqftJJ8Ic/NJ5X\nES979oQKqUdxsZvjkSgnneRGOBktCHMTZZxmu7mNXwwOH3b+eW+s+oQJoZPIPvrI9W4TwS8G4fEC\nCIqB50Lyi0VRketdr1yZWTEQcSLwySeuPMlaBuFC5w8ib9/urK127ZKf2+G5nNLFt74VOyZkFCBm\nGWScZisG/pjB3r3OtdSqlTs+4wx4663g+jcLFsQ3bt9PaalrEA8fjmwZdOzoOiq1tZGvl5W5CVWZ\nFANwvfrKSjeTOhkxiCRY4WLQrZub+5Gsq8gL7qeLM8+ESy9N3/2MZoBZBhmn2YqBP2YQ7pP2Rpms\nWeOurVnjhp0mQuvWrsGvrY1sGYDrjW/YEFkMSkudGGQyZgBBF0+yYrB/f3xiMGhQ8mIQzU1kGHFj\nlkHGadZi4FkG4WIgErQOFi92I4eKixPPo6lA8IABrsGNJQbZsAzatHH5+H8jBw+6wLI3VyMSXtli\nuYk8IRw0yMVh3nor8TJGCyAbRtxkSwxqa0OX1G1BFIQYROp5TpjgGq5kXEQenhhEswy83nI0N1Gm\nA8jgxKBHD+cm81sGVVVBoYiGV7Z4LIMLLnDuqPPOC05yixezDIyUOHzY/S0qyvxCdV/9auMZri2E\ntIiBiEwSkeUislJEbo9wvVhEZohIhYi8LyIxnBfx0bWrG8554EDknuc557j1+m+6KbhMRKJ4i9VF\nsww8P3o0ywAyLwY9erhGv7g4VAwqK2NbBf6yRRpaWlvrfoPV1e57OPNMt5ppnz7u+04EswyMlPCs\nAsi8ZVBX59araYGkLAYiUgQ8BJwPjACuEJFjwpJdA1Sr6tHAfwH3p5pvUVFwpm+kxuaoo1zjqNp4\nR7J48Rari8cyCL+eLTHo08et9RMuBnv3Nt0b98oWvo6TZxnU1bnv1W9ddOjgevqJYJaBkRLhYpDJ\ntYn27nVD7Fog6bAMxgEVqrpe3W5mM4CLwtJchNvqEuB50rTlpecqylTP0x8ziCYG69ZFdxNB5gPI\nkyfDo4827jDt29d0nKRXL/cdtg6bbeKJQaR6tW+fuGWQ7qGlRgvDG0kEmbcMvvjCxCAF+gIbfMcb\nA+ciplHVQ0CtiKS8XYo3vDRTjY0/ZpCMm8hbzjqTFBe7Rj3cMti/v2kx6NfPzcEIpykxSNQySPfQ\nUqOFkU03UQu2DHI1AznmPNZ4l/otK3ONVqYam/793XLaqpH976Wlzq++alVkMSgrC859yDSRxCCe\nZT0iDT1Ntxjk0k2UiaV+jSyTbcvg4MH4lxsuINIhBpWAv0npFzjnZyNwJFAlIq2ATrFWd4x3qV8v\nZnDgQGYsg+98x23E0rp15MZMxLmKlixp3Gj26OF67NkiXAzicRNFwxODSBZRhw7NK4CciaV+jSyT\nbctgzBhnHYTvFVvgpMNNNB84SkQGiEgxcDkwMyzNP4ArA+8vBd5MQ74Nq4NmqrERcesMxerVDhzo\nxMK/Iie41TVnz05/maIR/huJx00Uja5dnVXwl780jpU0N8vAKACyJQYHD7q/J5zQIl1FKVsGqnpI\nRK4H5uDEZbqqLhORacB8VX0ZmA78RUQqgB04wUgZTwy8PYlzwaBBoZvDeIg0PbQznSTrJopEu3Zu\n7+aaGvjyl0OvJRszsACykTTZchN98YXbkHz4cBODZFHVV4FhYeem+t7vAy5LR15+ysqcv141d+49\nTwxyTTrdRCLwgx9EvmZDS42sky3LYO9e1xMaPRpuuw2eeqpxmr594emnnelfYDTbJawhaBkUFzt3\nTS4YMqTxEtC5IJ2WQSwSHVqq6jpcJgZG0hw4kF3LYOJE17B4M5/9vPKK20Fr9GhnQfzXf2WmLDmg\n2S5HAZmPGcTDBRfA//5vbvL2Ez5LPxXLIBaJuom++MKJUlEzfNJEZLqIbBGRJTHS/D4ws36xiBRe\ndzEfOHgwu5YBuFUqO3du/LriCvjwQ7jqKrdjVgHRDH+iQTwxyOWkplatsjtqKBrhEzNTCSDHIlEx\naOYuosdxM+sjIiIXAEMCM+t/CPwpWwVrUWTLTeRZBk0xaBBMmhTcR7dAaPZi4C1H0YwbnLSQLTdR\nokNLm3PwWFXfBWpiJLkIeCqQdh7QWUTywGlYYPgDyOEPejrZuzc+MQCX7osvMlOOHNGsxaCsLPdu\nonwhnQHkWLQwy6ApwmffV9J49n3yzJzptrFr6WTTMvDcRE1RUuLEo4CWu27WAeR27YLbTpoYpG+e\nQSwSFQMT6iDxzqxv4Jln4PTT4bjjMlquvCebo4nitQyKioJlycQPrQkyMbO+WYuBt/bPxo3W4ESK\nGWRqNJFZBg1U4mbWe0Safd9AvDPrG6irKzi/dFJka56BP4AcD56rKAdikImZ9c3aTQRODOrrTQyy\n5SZKNGZQACuWCtHX0poJfBdARE4BalV1S9pyrq8vOL90UuRbANnDcxUVCM3aMoDgqqAF3PuMi2RW\nLU2GZNxEzfV/IyJPAxOBMhH5HJgKFAOqqv+jqq+IyIUisgrYDVyd1gKYZeDI1jyDRC2Ddu0K6v/T\n7MXAWzunmfc+UyZSzMDcRKmhqt+KI831GStAXZ1ZBpC9eQbJWAYF9P8pCDeRSGKCXoh4MQNvcEMm\nRxO1lKGlOcfEwJHt5SjipcDcRCmJgYh0FZE5IrJCRGaLSOco6Q6JyEIRWSQif08lz3BKS10DFb5Q\nXEujVSv3HRw65I4z5SZKdG2i5mwZ5JTDh90m3wXU2CRNtheqi5d27QpKrFO1DH4OvK6qw3DLUv8i\nSrrdqjpGVUer6tdSzDMETwyM0LhBvriJzDJIkl27ggs7tXTycWgpmGUQhn9v4yeBaA19xvrtpaXW\n2Hj44waZchMVF7uOmrf0e1OYZZAkdXXur4lBdheqMzdR0vTwhtKp6magR5R0bUXkQxF5T0QuSjHP\nEEwMgvjnGmTKTSSSmKuoAIaW5gZPDAqosUmabC5U14LdRE2OJhKR1wD/eisCKPDLCMmjzc0eoKqb\nRGQQ8KaILFHVtdHyTGSmZlmZNTYe2XATQdBVFL67WyRyPbS02e6BXF/v/hZQY5M0kdxEqukPFLbw\nAHKTYqCq50a7Fljet6eqbhGRXsDWKPfYFPi7VkTKgdFAXGLQFOPHw8MPx528oPGLQabcRJBY3CDX\nbqJmuwdyXZ1bMtnEINRNVFTkRkv4rYV00cInnaXqJpoJXBV4fyXwUngCEekS2BsZEekGfAlI255y\nbdvCiSem627NG3/MINOWQbzDSy2AnCR1dW7XpAJqbJImvOHPlKsomUlnBSTWqU46uw94VkS+B6wn\nsLWliIwFfqiqPwCOBf4sIodw4vOfqro8xXyNCPhjBpm0DLyYwfz5sGlT4+uHDsHnn0NlJSxZYgHk\npPDEYPPmXJck9xw4ENpjz5QYtHDLICUxUNVq4JwI5xcAPwi8fx84PpV8jPgIjxlk0k20ezd8/evO\nTRe+i5kIHHkk9OsH118PJ5+cmXIUNPX1btektVG9qS0H/zwDyB/LwMTAyFeyGUCuqHD3f/nlzOTR\n4jE3URB/ABnyxzIoMDdRs1+OwgiSjXkG4MTgww9hxIjM3N8gKAYF1NgkjT+ADJm1DFqwm8jEoIDw\nr08U3plKJx06ODEYOTIz9zdwYtCrV0E1NkmTrQByMpPOCkisTQwKCM9N5AlBuC8/XbRvD0uXmhhk\nlPr64JK8mZpk1RTTpsE77+Qmbz/ZchMlM+msgMTaxKCA8MQgky4icGKgamKQUbx5BrnsfZaXw2dp\nGwWePNl0E7XgALKJQQHhxQwyOZIIgkNFLWaQQerq3BTvXPY+q6qgtjY3efvJppvI9jMwCgEvZpDJ\nkUTgYgZHHhnfchTNHRGZJCLLRWSliNwe4fqVIrI1sET7wsCcm9Spr3eWQa5GrKi6iSL5IAbptAx2\n74Y334R33w1u/uFhO50ZhUI23UQtwUUkIkXAQ8DZQBUwX0ReijBpcoaq3pjWzHPtJtq50zWc+SAG\nkSwD/x6vifD003DXXa5eH38MRx3lzh865AQmkR9OgbmJTAwKCE8MMu0mOvtsOProzN0/jxgHVKjq\negARmYFbtj1cDNK7Yppq7t1EVVXubyQxePBBN6Pw61/PTlnSGUCurYVvfcvFQ6qrg+f37XPmdCKL\n35mbyMhX/DGDTLqJhg+HSZMyd/88oi+wwXe8MXAunEtEZLGIPCsi/VLOde9e1yiVlGTOTXTggHNF\nRWtUY4nB4sVuOFm2CHcThW/4nQj19dCxI3TpElq3RIeVQsG5iUwMCohsuYmMEGYCA1V1FPA6wc2e\nkseLF0Dmep/nnQc9esDFF0e+XlkJAwdGFoP6eqipSX+ZopHOAPLOnc7i6tIluGcEJD6sFMxN5EdE\nvgHchVuM7iRVXRgl3STgv3DiM11V70slXyMy2QogtyAqgf6+436Bcw2oqr9VfBS4P9rN4t6nw3MR\nQeZ6n1VV8Nhj8OtfR78+YgSsXh25fNmMJaTTTZROyyCHbqJM7NORaszgE+Bi4M/REiQQhDNSxCyD\ntDMfOEpEBgCbgMuBK/wJRKRXYJc/cPGEqAPz496no6YGunZ17zPlJqqthQEDQnvHfqqqnD9wwYLG\n1+rrsy8GiY4mev5597195zuh5+vrndB27hxah2Qsgxy6iTKxT0dKbiJVXaGqFcQOoDUE4VT1AOAF\n4Yw0k615Bi0FVT0EXA/MAZbiRg0tE5FpIvKVQLIbReRTEVkUSHtVyhnX1gbFIFO9z7o6JwbRGnVP\nDCJdr6vLjpvok0+cWycZN9GHH8K8eY3P+91E4WLQjCyDTJCN0USRgnDjspBvi6O42O0zYG6i9KGq\nrwLDws5N9b2fAkxJa6Y1Na6xgsz0PvfudSOWevZ0DXukLSQrK2HIEDfkMrzXXFeXeMOZDLfeCtdc\nk5xlUFMTeTs+v5vIvxlHohPOoOXFDGLsgXyHqv4jUwUzEseLGZibKE8JX5c/Gn7LIJKbaMECmDnT\npbnppsT3AvbmMLRp4xq03bvhiCNC01RVQd++wUCrv6Gsr89Ob2PHDjf8M9wyaNu26UY4lhhECyAn\nKnBt2zpROnw4cwuBZZGU9kCOkyaDcOHEHWgzQsjWPIPmRCYCbUlTVxdcfC4WfssgUu/zN79xPfZX\nX4XLL3ermyZaDm+0kuc794uBqus19+kTdKf0DPQHDxxw4pSNmMGOHe4VHkDu0QO2RtxuPUg0MYjk\nJpowwc1IvuKKxuljIRIUpgLYzi+dbqJo3ZMmg3DhxB1oM0LI1jyD5kQmAm1JU1MTnxg0ZRnMmwf/\n+Ads2QKffpq4GNTWBsXGaxT7+aZHbN/uxKGkpLFvvb7enauvd4LUqlVieSdCdXVQDPwWVe/eTc9z\nqKmJ7M/33ET+APKnn8KGDYl/jxB04xWAGKRk24jI10RkA3AK8LKIzAqc7y0iL0P0IFxqxTYiYaOJ\n8px4g67+xjpcDLZtcw3kMce4NUE+/TTxcvgtg3B3CTiroHfv4HW/GNTVuXOdOkUfiZQOvElxkdxE\nvXpF3nzbT01N4+9bFXbtCh1aum+fO9erV3KungKKG6S6B/Lfgb9HOL8J+IrvuFEQzkg//nkGJgZ5\nSLyulXA30bZtwWvz5sFJJ7mGa+RI+OijxMsRyU3kZ8uWoFsokhh07uxcJDU1UFqaeP7x4C0VEclN\n1KsXbN4c+XMenmXgD47v2eN+GK1bB0Vw0yZX12R9/gU0oqj5Rz2MBvwxA3MT5SGJWAbR3ETz5sHJ\nJ7v3yVoGtbWhlkG4GGzdGl0M/AHYTMYNwsUg3E0USwwOHXKxAQjttXtlh2D5/VZQMhTQkhQmBgWE\nFzMwN1GeEq8YhFsG0cRgxAjnOz98OLFyeK4eiC4GPXpEvu5ZBl27ZlYMduxweURyE/Xs6cQgfAlq\nDy8uUFoa+p17wWNwMZEvvnCxglTEIN1uoqlTXcwmB5gYFBA2mijPSdYy8Bqbw4fdZCpPDLp2dQ3z\n558nVo6mYgZNiYFnGWRy4tmOHW5p3EiWQYcOThyixSy8GdzhZfREApzrqFMnWL7cjZpKlnS6iVTh\nd7+DJUvSc78EMTEoIDp2dL9bcxPlKfH2pKMFkJcvd71dr6GG5FxFfjdRpJhBLDHwFtHLhpvoqKNc\nY75/f6hlALHjBp4YhFsvfssAXB2WLcsfN9GmTa6MlTFH3mcME4MCYtgw116YmyhPiacnrRoqBv6e\n5/vvw6mnhqY/7jj45S/hkkvc629/azqPptxEW7bE5ybKtGXQq5cbshlJDGLFDfxiEM0ygPSIQTTL\nYMYMNw8kEZYFBllu3Jh8eVLAxKCA6N7dWb8bNpgY5CXxNJ67d7t/nvcP9Pc8I4nBrbfCv/87fPvb\nrif9zDNN55GKmyibAeSysuC8jPD5DPFYBpHcROGWwYoVqbuJwi2DXbvghhvg740GWsZm+XLnDjPL\nwEgVEbe22OLF5ibKS+IRA3/wGELdRJHEoGfPoFXwzW9CRUXTecTjJvKPJvKLRTYtg9JSJwatWzde\nciPWXIN43USdO7uGPFU30ZYtwdFLAA895M5HWv47FsuXwymn5EwMbNvLAuPYY+Ff/zLLIC+prXVz\nBh56KDgC6Ljj4LLLQtN4wWMIuiFqa2H9ejj++Oj3P/poWLUq8sJzfmK5iVQbWwb+Rt8LIIvEbxnU\n1cHs2e7w6DG2AAAgAElEQVShnDw5vjH9O3Y4ISgtjbyeUzxuorZtm3YTefdKlmHD4O67Ydo0973t\n3w+//a3bK+KGGxK71/Llbk/Zf/4z+fKkgFkGBcbw4e73bJZBHlJTA889B3PmuIZRFa67LnSIZCTL\nYO9eN6R07NjGvnM/nTq5kTbelpXRiOUm2rXLNfQdOrjjfv3caKVDh9xxMgHkn/zEjZL5wQ8i748Q\nierqoGUQqc7pchMVFTn/arL88pfOMigudmK9dKmzqs4/3/0f9u2L/16eGJibyEgHxx7r/pplkIfU\n1MDcuXDttXDnnXDPPa5x+sy3H064ZdCunWv0pk5t7CKKxNChTbuKYrmJ/FYBuPL17g0rV7pjzzKI\n1020erVbR2nWLPjqV+MXA88ySEUM4hlN1KtXetZXOuEE+PhjWLQIRo92ZT7ySFi3Lr7P79zpBPDk\nk908g2R3cksBE4MCY/hw99fEIA+pq3NicOaZwXMTJoB/VVX/SCJwwc1HH3W96htvbDqPo4+OLQaq\noXssh/fww8UAXOO2eLF7H69lsHOnE4DbboPrr3fpx45NXAxiuYniiRk05SZKxUXkJ1wMwAX0440b\nrFjhhLy4GLp1a3q5jQxgMYMCo18/N7nS3ER5SLt2rnHzrxA6cSK8/LJzF0HolpfgXDaJLK3clBjs\n3u0eDq+3XVLiBMLbwMYfPPYYNco1cldcEXQxeWsmLVzoGj1/jxuccH38sbt2003u3NixTtjioSk3\nUd++zi0TaYal9x22aRN9BjKkXwyee8414t/4hjs3ZIiL4YD7jh9+uLELqE8f97388pcueOzVrbLS\nWRZZJNVVS78R2PLvkIiMiZFunYh8LCKLROTDVPLMFLla8z7d+Yq4eUieyzebecdL3uwvEAciMklE\nlovIShG5PcL1YhGZISIVIvK+iPSPdB/ANVBnnRV6bsIEeOutYNwg3DJIlKOPDrp0AHbtovyFF1wv\neseOUBeRq4A79uIGkSyDUaOCloHnJiorc8thfO97rsG/5x544gmX94YNMHMm5dOmwbPPBsXt+OOD\nE2GisWqVy+PwYTfHIJoYdOvmGuBXXml8raaG8jVrgm6irVtdDz3cMvjKV1y508Hxx8OiRZQvXOi+\nL3DfiycG8+bBAw+4787/WrwYrrwSJk2CP/zBpe3XLydxg1TdRJ8AFwNvNZHuMDBRVUeral5ueVlI\nDePLL8P48bnJOx6aixiISBHwEHA+MAK4QkSOCUt2DVCtqkcD/wXcH/WGkcRgwABnMXz8sTsODyAn\nit8yWLoUBg6k/MorYcwYl9cTTzS+v9/l459w5uG5iTwXU6dOrjf+1lvu/Ny5Llj6+utw2mluq8rv\nfY/ycJdQu3aufJ98Erns//ynG1116aXOKhCJ7iYCuOoqV59wamoo/+yzoJvommtcvKK6OtQyKCtz\ngpIOhg6FqirKvT0gIFQMHn7YWUu/+EXoa/p09/+67bZgPT3LIMukJAaqukJVK4i+sY2HpJqXET9l\nZQWxC18+MA6oUNX1qnoAmAFcFJbmIuDJwPvngbOj3u2Pf3RDK8O55RY491yYMgXeey/UTZQoRx0F\na9a4EUtf/aobxXPrrc4yeOUV+I//CLUMIFQMIlkGniuloiLUxeQxYoSr2//+Lzz/vGvIfvKTyOWL\nFDdYutQ1lN/7noszeBPOILplAM4dU14eusQ3OAHwGuWNG91aP927O1dXuDsrXbRu7Uxy/wY5nhhs\n3eqC6FdfHd+9+vZ15VaF3/8eTj8dXnjBucX8rx070luFtN4tOgrMFhEF/kdVH8lSvoaRCn2BDb7j\njTiBiJhGVQ+JSK2IlKpqdaO7RTPXbrgBzjkHnnzSxRD8AeZE6dDBTUC7914XuP3Od8DbOfCMM9xM\n5fDA65lnwte/7nra77wT9F17iDjXx8SJ7h6xOOOM2LuQjR3r3CX//KcTvQ0bnOvo+993ezMceaTb\n3/mDD1z644+P7srp1MmV+5hjQvdo9ja379jRlf2uu1z5x47NnBh4dfMsAYBBg1z9Roxwc0ni3fuh\nXz83k3zZMjca6fbb3dyF8GUqamvd8MHf/CY95VfVmC/gNWCJ7/VJ4O9XfWnmAmNi3KN34G93YDFw\nWoy0ai97ZfLV1DPvexa/juu8eMffBn4fluYToI/veBVQas+1vbL9ive5jvZq0jJQ1XObShPHPTYF\n/m4TkRdxvat3o6RtyuVkGNmiEvAHhPsFzvnZCBwJVIlIK6BTJKvAnmsj30mnZzniwy4i7UXkiMD7\nDsB5QBLbMxlG1pkPHCUiA0SkGLgcmBmW5h/AlYH3lwJvZrF8hpE2Uh1a+jUR2QCcArwsIrMC53uL\nyMuBZD2Bd0VkEfAB8A9VnZNKvoaRDVT1EHA9MAdYCsxQ1WUiMk1EvD2+pwPdRKQC+Anw89yU1jBS\nQzTa1nGGYRhGiyFvBiA2Nbknjfn0E5E3RWSpiHwiIjcGzncVkTkiskJEZotI56bulUIZikRkoYjM\nDBwPFJEPAnV/RkTSPspLRDqLyHMisixQ95OzVWcR+WlgcuISEflrYKJW2ussItNFZIuILPGdi1pH\nEfl9YLLYYhEZlWr+McrVIp7tXDzXgXxy8mxn67kO5JXxZzsvxEDim9yTLg4CN6vqCOBU4LpAXj8H\nXlfVYTi/7y8ylD/ATYBvdTLuA36jqkOBWtxEpnTzIPCKqh4LnAAsJwt1FpE+wA240WbH44YzX0Fm\n6vw47hnyE7GOInIBMCQwWeyHwJ/SkH8jWtiznYvnGnLwbGf5uYZsPNupDkdKxwsXc5jlO/45cHuW\n8v47cA7uAeoZONcLWJ6h/PrhhutOBGYGzm0DinzfxatpzrMTsDrC+YzXGegDrAe64n4wM4Fzga2Z\nqDMwAFgSo47LAu//BHzTl26Zly7N9W8Rz3YunuvAfXPybGf7uQ7cL6PPdl5YBkSe3NM305mKyEBg\nFC6w3VNVtwCo6magR/RPpsTvgJ/hxgYjImVAjaoGdjthI+5BSyeDgO0i8njAjP8fEWlPFuqsqlXA\nb4DPccMy64CFQG2G6+zRI6yO3ips4c9cJZl55lrKs52L5xpy9GznwXMNaX6280UMso644a7PAzep\n6i4CD7GPtEfWReTLwBZVXUzoUNxMj0FvDYwBHlbVMcBuXA81G3XugluyYQDuh9EBmJTufBKg4EdM\nZPvZzuFzDTl6tvPwuYYU65gvYhDP5J60EQjqPA/8RVVfCpzeIiI9A9d74cy9dDMemCwia4BngLNw\n/s7OAd8yZKbuG4ENqvpR4PgF3A8oG3U+B1ijqtXqhmq+iPseumS4zh7R6liJmyzmkakytIRnO1fP\nNeTu2c71cw1pfrbzRQzimdyTTh4DPlPVB33nZgJXBd5fCbwU/qFUUdUpqtpfVQfj6vimqn4bt5zH\npZnKO2BKbhCRoYFTZ+PGzWe8zjgz+hQRKRER8eWdqToLoT1Sfx2v8uUzE/gugIicgjPvt6SpDH4K\n/tnO1XMdyDtXz3a2n2vI9LOd7oBOCsGRScAKoAL4eQbzGQ8cwq2RtAjn55sElAKvB8owB+iS4fpO\nIBhoGwTMA1YC/we0yUB+J+AapsXA34DO2aozMBUXxFqCW+GzTSbqDDwNVAH7cD/Wq3EBvoh1xI3y\nWQV8TIy1tezZzt/nOpfPdrae62w92zbpzDAMw8gbN5FhGIaRQ0wMDMMwDBMDwzAMw8TAMAzDwMTA\nMAzDwMTAMAzDgKa3vcw2ImJjXY2MorYFpWE0Ii8tg0xNiIn1mjp1aovKt6XW2TCMyOSlGBiGYRjZ\nxcTAMAzDMDHwmDhxYovKN5d557LOhmFEJu/WJhIRzbcyGYWDiKAWQDaMRuTdaCKj+fD+++9TWdl4\nmfTWrVvz5S9/mTZt2uSgVIZhJENCloGITAL+C+demq6q94VdvxJ4ALfhBMBDqvqYiJwA/DfQEbfE\n7q9U9dkoeZhl0EwYMGAAI0aMoH379iHn586dy8svv8ypp56ao5JFxywDw4hM3JZBYPeeh3CbOFQB\n80XkJVVdHpZ0hqreGHZuN/AdVV0tIr2BBSLyqqrWp1J4I7fs37+f6dOn07t375DzZ511Fnv27MlR\nqQzDSIZEAsjjgApVXa+qB4AZuD1Aw2nU61LVVaq6OvB+E257tu5JlNfIIw4cOBDRFdS+fXu++OKL\nHJTIMIxkSUQM+gIbfMcbA+fCuUREFovIsyLSL/yiiIzD7f6zOrGiGvnGwYMHad26sXHZrl07swwM\no5mR7qGlM4GBqjoKtx3bk/6LARfRUwT37TSaMWYZGEbhkMhookqgv++4X+BcA6pa4zt8FLjfOxCR\njsDLwC9UdX6sjO66666G9xMnTrRx6XlKLMsgX8SgvLyc8vLyXBfDMPKeuEcTiUgr3MbLZwObgA+B\nK1R1mS9NL1XdHHh/MfAzVf2SiLQBXgVeUtXfN5GPjSZqJhQVFXHw4EGKikINzJtvvpl+/fpx8803\n56hk0bHRRIYRmbgtA1U9JCLXA3MIDi1dJiLTgPmq+jJwo4hMBg4A1QTdQZcBpwFdReRqQIGrVHVJ\n+qpiZJNDhw4hIo2EAPLLMjAMIz4SmnSmqq8Cw8LOTfW9nwJMifC5vwJ/TbKMRh4SzUUELmZgAWTD\naF7Y2kRGUkQLHoNZBobRHDExMJIilmVgQ0sNo/lhYmAkRSzLwIaWGkbzw8TASIqm3ERmGRhG88LE\nwEiKpgLIZhkYRvPCxMBICrMMDKOwMDEwksIsA8MoLEwMjKQwy8AwCgsTAyMpbDSRYRQWJgZGUtg8\nA8MoLEwMjKQwy8AwCgsTAyMpzDIwjMLCxMBIiliWQUlJCfv27ePw4cNZLpVhGMliYmAkRSwxKCoq\nom3btuzduzfLpTIMI1kSEgMRmSQiy0VkpYjcHuH6lSKyVUQWBl7f812bJSI1IjIzHQU3ckssNxFY\n3MAwmhtxi4GIFAEPAecDI4ArROSYCElnqOqYwOsx3/n7gW+nVFojb4hlGYDFDQyjuZGIZTAOqFDV\n9ap6AJgBXBQhXcQtBVV1LrAr8SIa+YhZBoZRWCQiBn2BDb7jjYFz4VwiIotF5FkR6ZdS6Yy8JVXL\n4Pe//z0fffRRJopmGEYSJLTtZRzMBJ5W1QMi8gPgSeDsRG9y1113NbyfOHEiEydOTFf5jDTRlBg0\nZRnMnj2bkpISTjzxxEwUr4Hy8nLKy8szmodhFAKJiEEl0N933C9wrgFVrfEdPoqLEySMXwyM/KQp\nN1FTlkFtbS3V1dWZKFoI4Z2JadOmZTxPw2iOJOImmg8cJSIDRKQYuBxnCTQgIr18hxcBn4XdQ4gS\nUzCaF6laBnV1dVkRA8Mw4iNuy0BVD4nI9cAcnIhMV9VlIjINmK+qLwM3ishk4ABQDVzlfV5E3gaG\nAUeIyOfANar6WvqqYmSTVC0DEwPDyC8Sihmo6qu4Bt1/bqrv/RRgSpTPnpFMAY38JB2WwY4dOzJR\nNMMwksBmIBtx8f7773PFFVdwzTXXsH///pRGEx06dIidO3dGtAzuvvtu2rVrF/Lq27evLW1hGBnG\nxMCIi9dff519+/bxwgsvsH379pTmGdTX1wNEFINVq1bx4IMPUl1d3fCqqamxCWyGkWFMDIy4qKmp\n4Utf+hJdunRh3759KVkGdXV1tGnTJqIY1NXV0b179xDLoFOnTuzcuTNtdTEMozHNWgwqKys5dOhQ\nrovRIqipqaFr164NC9AlYhls376d2267jVtvvZXKykrq6uoYMGBARDGor6+nU6dOIec6duzYYE0Y\nhpEZmrUYXHLJJcydOzfXxSgo/vSnP/GPf/yj0fna2lq6du3asDx1IpbB888/z3vvvUd5eTlvvfUW\ntbW19OrlRiGHWw91dXV07tw55FzHjh1TsgxqamqYP39+0p83jJZAsxUDVWXFihVs3Lgx10UpGNat\nW8dPf/pT3njjjUbXampq6NKlS4NlkMhoolmzZnHddddx+umns3nzZurq6ujSpQulpaWNrINIlkGq\nbqJXX32VX/3qV0l/3jBaAs1WDLZv305dXR2bNm3KdVFSYvfu3dTU1OTFaJlbb72VY489lsrKykbX\nwi2DeOcZ7Nu3j/Lycs477zx69erVIAadO3eOKAbRLINU3ERvvfUWEyZMSPrzhtESaLZisGrVKoAG\nMTh8+DB79+5tVjGEnTt3UlZWRp8+ffjtb3+b07JUVVUxd+5c7r///ohiEB4ziNcyePvttxk5ciRl\nZWX07NmzSTGIFjNIxTJ4++23OeMMm+ZiGLFotmJQUVFBSUlJgxhceumldOzYkZNPPjlrZfjoo49S\nGvK4detW+vTpw4MPPsjy5cvTWLLEWbBgAePGjWPIkCFRxaBLly4JWQaff/45jz/+OBdccAEAvXr1\nYsuWLdTW1tK5c2fKyspCxMDbKrOkpCTkXqm4ibZu3UpVVRUnnHBCUp83jJZCsxaDU089tUEMlixZ\nwkcffcTSpUuztt3ixRdfzN13383hw4e58cYbWb16dcM1VWX79u0xP79jxw7KysoYNGgQa9euzXRx\nY7JgwQLGjBlDnz592Lx5c4jb6uDBg+zZs4dOnTrFbRkce+yxdOnShR07dvDNb34TIMRN5MUM/LOQ\nPatAJHT5qlTcRO+++y7jx4+nVatWSX3eMFoKeSkGBw4caDLNqlWrOP3006mqquLQoUN8/vnnDBs2\njKOOOorPPgtfHy/9bN26lbq6Oh555BF+9KMf8Yc//IE333yz4fq8efOa9FN7YjBw4MC8EIOxY8fS\ntm1bOnfuzNatWxuu1dXV0alTJ4qKiuIeTTRkyBBmz57N7NmzOfroowGajBlEihdA/JaBqnLw4MGQ\nV3l5ubmIDCMO8lIMamtrm0xTUVHBGWecwaZNm9iwYQM9evSgpKSEE044gY8//jih/F5++WXOPvvs\nhteMGTOa/MzChQs56aSTuO2225g5cyY333wzS5cubbg+f/58Vq1aFTOGsWPHDkpLS+nfvz+VlZUc\nPHgwoXKnE88yAOjbt2+Iq8hzEQFxzzOIRPfu3amurmbHjh0RxSBSvADiixmce+65lJSU0LZtW0pK\nShpef/7znzn33HMTKqdhtETyUgyaWs1SVVm1ahWjRo2iTZs2LFy4kMGDBwMkJQZz5sxh+PDhTJky\nhTFjxjBnzpwmP+M1nj/72c9Yvnw55557Lp9++mnD9UWLFrF///6I/ncPzzJo27YtPXr0iJnWT11d\nHStXrgw5t337dlS1UdpbbrmF//zP/4x4zWPTpk3s27ePAQMGAI3FwBtJBMRtGUSiVatWlJWVUVFR\nkZBlEO4meuGFF3jggQcajrdu3cr8+fOpqanh0KFDIZbBvn37GkTOMIzoJCQGIjJJRJaLyEoRuT3C\n9StFZKuILAy8vhd2baWIrBCR78bKpykx2L59O0VFRZSVldG7d2/effddhgwZAiQnBlVVVZxxxhmc\nffbZnHbaaU36+sFZBmPHjkVE6NKlCyNHjgwRg4ULF9KxY8eQOEI4nhgACcUNbr/9diZNmhTiTrvw\nwguZN29eo7SzZ8/mkUceYeLEiVxxxRUhr8WLFzeUdcyYMQ2++n79+jWyDDwxSMUyAOcqqqioiDjP\nIJplEO4mWrJkCU899VTD8bx58xg3bhzt27dPuDyGYTjiFgMRKQIeAs4HRgBXiMgxEZLOUNUxgddj\ngc92Bf4dOAk4GZgqIo27gAFiicFDDz3ED3/4Q4466iiAqGIQqyccTmVlJX37uu2cu3XrxrZt25r8\njOdj9+jbty979+5l27Zt7Nu3jxUrVnDBBRekXQw2btzIs88+S8+ePXnyyScbzm/ZsoX169eHpFVV\n1qxZwwcffMD111/P5MmTG1579+5l1qxZUevin8zndxOlYhmAE4O9e/c2jCbyB5BjWQZ+MaipqeHT\nTz9tEKx58+ZxyimnJFwWwzCCJGIZjAMqVHW9qh4AZuB2Mwsn0k5m5wNzVLVOVWtxG+RMipZRNDH4\n4osvmDJlCueddx4PPvgg4MRg0aJFDW6inj17UlxcnNDM5KqqKvr06QM4v3Y8o4BqamoaBAhARBg5\nciRLly5l6dKlDBkyhJEjR7JmzZqY9/HEIN4g8v33388111zDr3/9a+655x72798PuO/Mq/N9993H\n2rVr2bJlCx06dKBHjx5ceumlIVbB+eef35Df0qVLGTlyZEMesdxE8Y4mioa3DEWiMQO/m6i6upq2\nbdvy2mtub6QPPvggq0OKDaMQSUQM+gIbfMcbA+fCuUREFovIsyLiXQ//bGWUzwLRxeCtt95i1KhR\n/OhHP+LUU08FnBgcPHgwpGE+4YQTWLJkSczK7Nmzh71793L48GE2b95M7969gdiWwYoVK7jyyiuZ\nPHkyo0ePpqgo9OvzXEWLFi1izJgxDBkyJK2Wwauvvsqzzz7LLbfcwvjx4zniiCP49NNP2b9/P7t2\n7WoQg+nTp/P666+zdu3aBpEMZ9CgQQ1CVVFRwdChQxuuRQogh8cMUnETgRODnj17UlVV1XAt3tFE\nNTU1XHjhhcyePZvDhw8zf/58xo0bl3BZDMMIku4A8kxgoKqOAl4HnmoifUSiicE///lPLrzwwpBz\nXo/eLwa9evVqsnf/gx/8gAceeIBt27Y1jJ8H6NKlC7t37240vHXu3LmMHz+e4cOHc8cddzB9+vRG\n9xw5ciQLFy7krbfeYvTo0QwePDghMVi3bl3UtJ999hnf/e53eeGFFxoa1N69e7N9+3ZqamqA4Cqu\n69evZ8mSJaxZs4ZBgwZFvN/gwYNZu3YtqkpFRUXD8E+IbzRRKpZBq1at6NChAz179uSLL75o6PV7\nQ1jDCXcTVVdXc/nll/Paa6/xxhtv0K1bN7p3755wWQzDCJJI164S6O877hc414Cq1vgOHwXu8312\nYthnoy43+sorrzQEMydOnMjEiRNRVV555RVefPHFkLS9e/emU6dOlJaWNpw74ogj2LVrV9SKbN26\nlWeffZbJkydTVVXVEC8AKCoqorS0lO3btzdYC+Cskh//+MfcfnujuHkDY8aM4aabbmLUqFHceeed\ndOnSJW430aBBg6ioqGD+/PkMGzasUaP45z//meuuu47x48c3nOvevTvbtm1rEM+NGzdSVVXF/v37\nWbJkCT169IhqGQwYMKAhfdu2bRt6/uACyOvWreOOO+7ga1/7GrW1tQ0jjVK1DHr27Ennzp0b/r+D\nBw9mzZo1jBo1ivr6+ojiFe4mqqmp4bjjjuPqq6/msssuY/LkyVHzKy8vp7y8POFyGkZLI5Ff83zg\nKBEZAGwCLgeu8CcQkV6qujlweBGwLPB+NvD/BYLGRcC5wM+jZVRfX8/UqVNDZqKuXLmS/fv3c9xx\nx4WkPfLIIxk6dGhI2qbE4LHHHmPUqFEsXbqUysrKBuvCw4sb+MVgxYoVfPnLX456T4BTTz2VPXv2\nNFgZ3iQoz82ydetWVq9eTVlZGUOHDg0Rg759+3Lcccdx6aWXcskllzRaq2jWrFn83//9X8g5z6VV\nXV1N7969qaysZO3atQwdOpQlS5YwaNAgTjvttIhlLS4upmfPnrz55pshVgE4F86vf/1r5s2bx69+\n9SuKi4vTGjPwu4I8V9qoUaOatAxUFRGhpqaG0tJSHnjgAe65556Yi/x5nQmPadOmJVxmw2gJxC0G\nqnpIRK7HBX+LgOmqukxEpgHzVfVl4EYRmQwcAKqBqwKfrRGRu4GPAAWmBQLJEdmzZw/Lli1j+PDh\nDeeeeOIJvv71rzdaquC0005j5syZIefCxeDQoUM89NBDDctUPPzwwzzzzDOcc845rFu3LsQygMhx\ng+XLl3PLLbfE/I5EpEEIvGPPVbRnzx4uu+wyBg4cyLJly9i8eTP79++nY8eOgBuD/9prr/H6669z\n9913h9x39erV7Ny5k1GjRoWc90SrurqakSNHUl5eTkVFBePGjeP111/n7bff5rvfjT6Kd/DgwSEz\nhP3lvvbaa7nssssYMmQIo0ePTttoopEjR/Jv//ZvIWXwXGn19fURYwbFxcW0bt2avXv3UlJS0mio\nq2EYqZNQzEBVX1XVYap6tKreGzg3NSAEqOoUVR2pqqNV9WxVXen77BOBzw1V1ZixhEsuuYQXXnih\n4XjXrl088sgj3HjjjY0rUFQU0oOHxmKwYcMGpk6d2rCn7m233cZpp53GwIEDKS8vj2oZeBw+fJiV\nK1cybNiweL6mEEaOHMkZZ5zBJZdcwl/+8hc++OADysrK+OCDDygtLW0kbqNHj2bx4sUhvd1Zs2Zx\nwQUXNErruYlqamro1asXXbt2Zd68eQwePJgTTjiBtWvXRo0ZgHNNzZkzp5EYeHTr1o3hw4fzr3/9\nK23zDLp37x4idkOGDGlwpUWzDCDoKtq9ezdt2rShuLg44bwNw4hOXs5AvuSSS/jb3/7WcPzEE08w\nYcKEqP7vcMLFoL6+nv79+3Pfffdx3333ccMNNwAwYsQI3njjjYiWgV8MNm7cSJcuXRp68YnwxBNP\nsHnzZiorKxuWRRgzZgyvvfZaSJzDo6ysjM6dO4eMLHrllVcaVv4ML6fnJiotLaVfv3688847DBo0\niOOPP55WrVpx5JFHRi3b4MGD2bZtW1QxADeZ7eDBg2mZgRwJ/4iraJYBBEcUeXU1DCO95KUYnHba\nadTX19OhQwc6dOjAbbfdxs9+9rO4Px9JDCL1OIcPH05tbW1Ey8DvJlq+fDnHHBNpfl3TtG7dOmS0\nEsDYsWN57bXXGuIF4YwePZqFCxcCLlj63nvvRVxfxx9A9sRg+fLlDZbBgAEDYvbePavBP6w0HG/0\nVjrWJoqE303UlGWwc+fOEBeRYRjpIy/FoFWrVixfvpytW7eydetWqqurE5phGq8YjBgxAqBJy2D5\n8uVJuYiiMWbMGBYsWBBTDBYtWgTAX//6Vy644IKGxtiPP2bQtWvXhnoMHjyYs88+m5/+9Kcxy+FZ\nWqwHU5QAAArrSURBVLEsg9GjR3Pttdc29MbTbRkMHDiQyspKDhw4ENMy8NxEJgaGkRnyUgwA2rRp\n02AZhG920hSJikFTlsGKFSuStgwiMWbMGFS1STFQVR555BG+//3vR0wXyU1UXFxMnz596NWrF9df\nf33McgwbNozx48fHdH8VFRXx8MMPN1gBqY4mCqe4uJjevXuzfv36qP8nMDeRYWSavBWDVIgkBpEa\nvKFDh3LKKac0mrDktwz27dvHsmXL0ioG3bt358gjj4wqBp7l8Oijj7Jz507OPPPMiOnKysqora1l\n27ZtDWIwcODARjOjo1FWVsa7776bUNlTnWcQicGDB1NeXk5JSUnUe5qbyDAyS4sRg0g9zrZt2/L+\n++83ajw9y+DBBx+kQ4cOfPTRR43mN6TK2LFjo4pBv379uPDCC3nxxRe5//77ozburVq1okuXLqxa\ntYrS0lJGjRrVaIZ2ukm3ZQBw5pln8h//8R8NS4xEolOnTtTX15tlYBgZIj1duzwjXjGIRrdu3aiq\nquLee+9l0aJFaRcCgDvuuCNqD1dEeOKJJ+K6T7du3Vi5ciWlpaUcffTR/O53v0tjKRuTCcvgzjvv\n5M4774yZxiwDw8gsBWsZ+NeySUYMtm/fzpgxYzIiBAAnnnhiyHpKydK9e3dUNWu95UxYBvHQuXNn\nqqurTQwMI0MUpBh06NCBXbt2NexpkKgYlJSU0Llz55jrEOULXrwj0mijTFBSUpLWoaXxcvLJJ/P2\n22+bm8gwMkRBuomKi4tp1aoV+/bto6SkhJ07dyYkBuB20+rfv3/TCXNMt27d6NKlC61atcpKfm3b\ntmXPnj20atWq0YzoTDJhwgQuvfRS9u7da5aBYWSAgrQMwPmYvbhBopYB0CyEAJxlkM2ecnFxMaqa\nVRcROItkwoQJLFiwwMTAMDJAwYqBP4icjBg0F7ItBt5ifNl0EXl4I6XMTWQY6acg3UTQcsSgW7du\nWW8cS0pK4p7LkE689ZnMMjCM9GNi0Mw5//zzs+7SytWy0QMHDuS5554zy8AwMkBC3TsRmSQiy0Vk\npYhEHWojIl8XkcMiMiZw3EZEHhORJSKySEQmpFrwpmgpYtC9e3fOOOOMrOYZa6ZwpvnGN76R1cC1\nYbQU4v5Fi0gR8BBwNlAFzBeRl1R1eVi6I4AbgQ98p/8foKp6vIh0B2YBJ6Za+Fh4YnD48GF27drF\nEUcckcnsWhRt27a1BtkwCoxELINxQIWqrlfVA8AM3NaW4dwN3Avs850bDrwJoKrbgFoRyYoY7N69\nm3bt2mVt6GVLIJeWgWEYmSERMegLbPAdbwyca0BERgP9VHVW2Gc/BiaLSCsRGQSMBaLvupIGPDEo\nZBdRrmjbtm3Wh5YahpFZ0ta9E+c3+C1wpf904O9jwLHAfGA98C/gULR73XXXXQ3vwzc0jxcTg8zh\nrU/UHCgvL6e8vDzXxTCMvCcRMagE/MNW+gXOeXQERgDlAWHoBbwkIpNVdSFws5dQRP4FrCQKfjFI\nFk8Mkpl9bMQmV/MMkiG8MzFt2rTcFcYw8phE3ETzgaNEZICIFAOXAzO9i6par6o9VHWwqg7CBZC/\nqqoLRaSdiLQHEJFzgQPhged0Y5ZB5igpKTE3kWEUGHF371T1kIhcD8zBich0VV0mItOA+ar6cvhH\nCLqJegCzReQQzpr4TupFj423cqmJQfqxALJhFB4J/aJV9VVgWNi5qVHSnuV7vx5I31ZhcWCWQeaw\nALJhFB4FvzaRiUH6MTeRYRQeJgZGwjSnALJhGPHRIsSgY8eOuS5OQWGWgWEUHgUtBp9//jn/+te/\nzDJIM2YZGEbhUbC/6CFDhnDllVeyf/9+zjnnnFwXp6Awy8AwCo+CFYN27dpx77335roYBYmNJjKM\nwqNg3URG5rB5BoZReJgYGAnTvn17SkpKcl0MwzDSiKhqrssQgohovpXJCMVbGrxnz565LkrCiAiq\napsxGEYYJgZGi8LEwDAiY24iwzAMw8TAMAzDMDEwDMMwSFAMRGSSiCwXkZUicnuMdF8XkcMiMiZw\n3FpEnhCRJSKyVER+nmrBDcMwjPQRtxiISBHwEHA+bkezK0Sk0bLUInIEcCNucxuPS4FiVT0eOBH4\noYj0D/9sLsnV1oi53JKxJdbZMIzIJGIZjAMqVHW9qh4AZgAXRUh3N3Av4N8kV4EOItIKaB+4Vp9c\nkTNDS2wYW2KdDcOITCJi0BfY4DveGDjXgIiMBvqp6qywzz4P7AE2AeuAX6tqbcKlNQzDMDJC2tYU\nEBEBfgtcGeHyOOAg0AsoA94RkddVdV268jcMwzCSJ+5JZyJyCnCXqk4KHP8cUFW9L3DcCVgF7MLt\nfdwL2AFMBr4HvK+qfw2knQ7MUtXnI+RjM86MjGKTzgyjMYlYBvOBo0RkAM7dczlwhXdRVetxG98D\nICJzgZtVdZGInAOcBfxVRDoApwC/i5SJ/VANwzCyT9wxA1U9BFwPzAGWAjNUdZmITBORr0T6CM5C\nAHgY6CginwLzgOmq+mlqRTcMwzDSRd6tTWQYhmFkn7yZgRzvhLY05NNPRN4MTH77RERuDJzvKiJz\nRGSFiMwWkc4ZLEORiCwUkZmB44Ei8kGg7s+ISNo3CxCRziLynIgsC9T95GzVWUR+KiKfBiYd/lVE\nijNRZxGZLiJbRGSJ71zUOorI70WkQkQWi8ioVPM3jOZMXohBvBPa0sRBXCxjBHAqcF0gr58Dr6vq\nMOBN4BcZyh/gJuAz3/F9wG9UdShQC1yTgTwfBF5R1WOBE4DlZKHOItIHuAEYE5h02BoXa8pEnR/H\nPUN+ItZRRC4Ahqjq0cAPgT+lIX/DaLbkhRgQ/4S2lFHVzaq6OPB+F7AM6BfI78lAsieBr2UifxHp\nB1wIPOo7fRbwgi/vi9OcZyfgdFV9HEBVD6pqHVmqM9AKN+mwNdAOqALOJM11VtV3gZqw0+F1vMh3\n/qnA5+YBnUWk+W3QYBhpIl/EoMkJbZlARAYCo3BLZ/RU1S3gBAPfyKg08zvgZ7gAOyJSBtSo6uHA\n9Y1AnzTnOQjYLiKPB9xT/yMi7clCnVW1CvgN8DlQCdQBC4HaDNfZo0dYHb0GP/yZqyQLz5xh5Cv5\nIgZZJ7CG0vPATQELITySnvbIuoh8GdgSsEz8Q2gzPZy2NTAGeFhVxwC7ce6TbNS5C64XPgDX4HcA\nJqU7nwSwEROGEYF8EYNKwL9wXb/AuYwQcFc8D/xFVV8KnN7iuQlEpBewNQNZjwcmi8ga4Bmce+hB\nnIvC+19kou4bgQ2q+lHg+AWcOGSjzucAa1S1OjA8+UXc99Alw3X2iFbHSuBIX7qMPnOGke/kixg0\nTGgTkWLchLaZGczvMeAzVX3Qd24mcFXg/ZXAS+EfShVVnaKq/VV1MK6Ob6rqt4G5uJVdM5J3wE2y\nQUSGBk6djZsrkvE649xDp4hISWDJEi/vTNVZCLW0/HW8ypfPTOC70DC7vtZzJxlGSyRv5hmIyCRc\nL7kINynt3gzlMx54G/gE5zJQYArwIfAsrre4Hrgsk4vpicgE4BZVnSwig3BB867AIuDbgUB6OvM7\nARe0bgOsAa7GBXYzXmcRmYoTvwO4+n0f1xNPa51F5GlgIm79qy3AVODvwHNEqKOIPIRzWe0GrlbV\nhankbxjNmbwRA8MwDCN35IubyDAMw8ghJgaGYRiGiYFhGIZhYmAYhmFgYmAYhmFgYmAYhmFgYmAY\nhmFgYmAYhmEA/z+uGx70JPLdZQAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f61c80b5150>"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Caclulate and plot autocorrelation function and autocorrelation time\n",
      "ap2, kp2 = C(p2/2.)\n",
      "ax2, kx2 = C(x2/2.)\n",
      "aE, kE = C(E)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.plot(ap2,'b')\n",
      "plt.plot(ax2,'r')\n",
      "plt.plot(aE,'m')\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEACAYAAAC+gnFaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd0E1faBvDnSnKTewU3wHSHXkxIKIEkBFJJI71u2mY3\nvSxkQxLSl2RTN213v/TGhlSSQEIoLrjgbnDvvVu2XGSrzfv9YWNc5CrZlu33d44P1szVnSshz6O5\nd+aOICIwxhibvGRj3QDGGGNji4OAMcYmOQ4Cxhib5DgIGGNskuMgYIyxSY6DgDHGJjmLBIEQ4kMh\nRJUQ4kQ/Zd4WQuQIIZKFEEstsV3GGGPms9QRwccANve1UghxIYBZRDQHwD0APrDQdhljjJnJIkFA\nRMcA1PdTZCuAzzrKHgfgKoSYYoltM8YYM89ojRH4Ayjp8risYxljjLExxoPFjDE2ySlGaTtlAAK7\nPA7oWNaNEIInPmKMsWEgIjHc51ryiEB0/JiyD8AtACCEWA2ggYiqTBVsKGgAEU36n2eeeWbM22At\nP/xe8HvB70X/P+ayyBGBEOIrABsAeAohigE8A8AWABHRf4hovxDiIiFELoAWALf3VVdTgRquM1wt\n0SzGGGODYJEgIKIbBlHmvsHUpSluNL9BjDHGBs3qBovbypvGuglWYcOGDWPdBKvB78Vp/F6cxu+F\n5QhL9C9ZihCCYh88iJA3N411UxhjbNwQQoCsZLDYIvQ1mrFuAmOMTSpWFwQGFQcBY4yNJqsLAlK3\njnUTGGNsUrG+IGjSjnUTGGNsUrG6IECLbqxbwBhjk4rVBYHQcBAwxthosrogkGkNY90ExhibVKwv\nCPTGsW4CY4xNKlYXBHK9NNZNYIyxScX6goAPCBhjbFRZXRDI+ICAMcZGlfUFAVldkxhjbEKzur2u\nHDLAiibCY4yxic7qgsAIe+ib2sa6GYwxNmlYXRDooERzOd+chjHGRovVBYERSrRUcBAwxthosb4g\nEPZoreIgYIyx0WJ1QSDBDm3VHASMMTZarC4ICLbQ1XIQMMbYaLG6IABsoK/jIGCMsdFidUFAJIe+\nvmmsm8EYY5OGRYJACLFFCJEphMgWQmw3sT5QCHFECJEohEgWQlzYd4uMoLpmSzSLMcbYIJgdBEII\nGYB3AGwGsADA9UKI+T2K7QTwPyJaDuB6AO/1WaGcAHWLuc1ijDE2SJY4IlgFIIeIiohID2APgK09\nykgAXDp+dwNQ1mdtCoJo5BvYM8bYaFFYoA5/ACVdHpeiPRy6ehbAQSHEAwCUAM7vqzJhKyBr4SBg\njLHRYokgGIzrAXxMRG8IIVYD+ALt3Ui9fGL4HIqaNPyxaxc2bNiADRs2jFITGWNsfAgNDUVoaKjF\n6hNk5kyfHTv2XUS0pePxDgBERLu7lEkFsJmIyjoe5wE4k4hqe9RFYTOOQFn3KVY2fmJWuxhjbLIQ\nQoCIxHCfb4kxgjgAs4UQ04UQtgCuA7CvR5kidHQHCSGCAdj1DIHOBjkqYKPTW6BZjDHGBsPsICAi\nI4D7ABwEkAZgDxFlCCGeFUJc0lHsMQB3CSGSAXwJ4Na+6pM720BhNJjbLMYYY4NkkTECIvoNwLwe\ny57p8nsGgLWDqcvGzRa2Rr5xMWOMjRaru7LY1sMONsQ3LmaMsdFidUFg52ELAQWIxwkYY2xUWF0Q\n2DjL0QpXtFbzfEOMMTYarC4IZEoZ2uDMdyljjLFRYnVBIFfKoZM5Q1PJQcAYY6PB6oJAppTBIHNE\nK9+ljDHGRoXVBYFcKYdBpoSWg4AxxkaF1QWBTCmDJHPgu5QxxtgosbogkCvlIJkDDCoOAsYYGw3W\nFwSOckDYQWrgIGCMsdFgdUEgU8pAwgak5iBgjLHRYHVBIFfKAdgATRwEjDE2GqwuCGRKGQQpIGvm\nIGCMsdFgdUEgV8ohJDkULRwEjDE2GqwuCGRKGWAUsGnlIGCMsdFgfUHgIAP0AjZtPOkcY4yNBusL\nAoUMkAP2Os1YN4UxxiYFqwsCABD2MtgZ2sa6GYwxNilYZRDIlHLYG3Vj3QzGGJsUrDIIFE5y2MMI\nSHzLSsYYG2lWGwQtcIWxsWWsm8IYYxOeVQaBXClHMzz4LmWMMTYKrDIIZEoZNDJ3DgLGGBsFFgkC\nIcQWIUSmECJbCLG9jzLXCCHShBAnhRBf9FefXCmHVu6GNr45DWOMjTiFuRUIIWQA3gFwHoByAHFC\niJ+IKLNLmdkAtgM4i4gahRBe/dUpU8qgU7igrUptbvMYY4wNwBJHBKsA5BBRERHpAewBsLVHmbsA\nvEtEjQBARLX9VSh3lEOvcIW2usECzWOMMdYfSwSBP4CSLo9LO5Z1NRfAPCHEMSFElBBic7+NUspg\nsHWFobLOAs1jjDHWH7O7hoawndkA1gOYBiBcCLHw1BFCV7t27YIqSYV6bQHOSvfBylFqIGOMjReh\noaEIDQ21WH2CiMyrQIjVAHYR0ZaOxzsAEBHt7lLmfQAxRPRpx+NDALYTUUKPuoiIULCrALmfxcBl\nynGcGf2mWe1jjLGJTggBIhLDfb4luobiAMwWQkwXQtgCuA7Avh5lfgSwEQA6BornAMjvq0K5Ug5h\np4SsgbuGGGNspJkdBERkBHAfgIMA0gDsIaIMIcSzQohLOsr8DqBOCJEG4DCAx4iovs9GKWWQ2Slh\n28hBwBhjI80iYwRE9BuAeT2WPdPj8aMAHh1MfXKlHHJbe9i3cBAwxthIs9orixUKWyjbOAgYY2yk\nWWUQyJVyyBUKuOg5CBhjbKRZZRDIlDLIIYej1AQYDGPdHMYYm9CsMgjkSjlkegkNcIOxts8xZcYY\nYxZglUEgU8ogtRrRIPNEYwF3DzHG2EiyyiCQK+WQNBIabT3RVMhBwBhjI2m0ppgYEplSBqPGCL29\nJ1qKOQgYY2wkWfURQZuTJ7TlHASMMTaSrDMIHOXtRwQuntDzDKSMMTairDIIhK0A6QmSqyeoloOA\nMcZGknUGgRDtE8+5e0GoOAgYY2wkWWUQAB0Xlbl5Qq7mIGCMsZFktUEgV8ph6+EBuyYOAsYYG0lW\nGwQypQz2Xq5QtvZ7e2PGGGNmssrrCID2IwIbDzc4afmIgDHGRpJVHxE4ejjB1VAHmHk7TcYYY32z\n2iCQK+VwcpBDB1uguXmsm8MYYxOW1QaBjY8NRKMOKniitZS7hxhjbKRYbRDY+dtBW6qF2sYT6nwO\nAsYYGynWGwQBdtCV6dBsxxPPMcbYSLLeIOg4Imh18ERbGQcBY4yNFOsNggA7aMu00Dl7QlvBQcAY\nYyPFuoOgVAuDmyeMVRwEjDE2UiwSBEKILUKITCFEthBiez/lrhJCSEKI5QPVaeNtA4PaALh5AXUc\nBIwxNlLMDgIhhAzAOwA2A1gA4HohxHwT5ZwAPAAgZlD1ygRsfW0hd/SGrIGDgDHGRooljghWAcgh\noiIi0gPYA2CriXLPA/gHAO1gK7bzt4ONvQdsGjkIGGNspFgiCPwBlHR5XNqxrJMQYhmAACI6MJSK\n7QLsYGPrDPsWDgLGGBspIz7pnBBCAHgdwK1dF/dVfteuXZ2/z6E5mCnOgGMbBwFjjJ0SGhqK0NBQ\ni9UnyMwJ3YQQqwHsIqItHY93ACAi2t3x2AVALoBmtAfAVAB1AC4josQedVHX9pS8XoLKODXmf7MG\nzka1We1kjLGJSggBIurzC/ZALNE1FAdgthBiuhDCFsB1APadWklEjUTkQ0QziSgI7YPFl/YMAVPs\nAuwgbwbsJQ2g01mgqYwxxnoyOwiIyAjgPgAHAaQB2ENEGUKIZ4UQl5h6CvrpGurKzt8OVKODCh4w\n1qjMbSpjjDETLDJGQES/AZjXY9kzfZQ9d7D1ts83pEWDzBP2RXVw9Z9qZksZY4z1ZLVXFgOAra8t\ndFU6qG280FjIA8aMMTYSrDoIZLYyKDwU0NgFQsMzkDLG2Iiw6iAAOrqH7P2hLecgYIyxkTAugkCy\nnQI9TzzHGGMjwvqDwN8OsPEE1XAQMMbYSLD+IAiwg5C7QtTVjnVTGGNsQrL+IPC3g5C7w60qc6yb\nwhhjE5L1B0GAHeQyJ/jVngD0+rFuDmOMTTjjIgjsW40oxjRQatpYN4cxxiYc6w8CfzsYK7VIVqyC\n6rfYsW4OY4xNOFYfBHJHOWT2MqiDVkN9OG6sm8MYYxOO1QcBANj620K5cCnsT3AQMMaYpY2LILAL\nsMOUM2bAoy4b0GjGujmMMTahjJsgmOEmIR0LYIxPGuvmMMbYhDI+gsDfDooGLTKcQlCzn7uHGGPM\nksZHEATYoa2oDY3zV6EljIOAMcYsaVwEgds5blD9poLj+pVwzuBTSBljzJLGRRAo5yphP80eM6ZM\nhbKpCqivH+smMcbYhDEuggAAptw4Be5JNUjCMuii4se6OYwxNmGMmyDwvtYb6v11yHdfjepfuHuI\nMcYsZdwEgd1UOzivcoZiytnQRfGAMWOMWcq4CQKgvXvIU5oK95xYQJLGujmMMTYhjKsg8LrCC/al\nOuTSIuCNN8a6OYwxNiFYJAiEEFuEEJlCiGwhxHYT6x8WQqQJIZKFEH8IIQKHsx2FswJeF3niC5tX\noX1+N5CcbH7jGWNskjM7CIQQMgDvANgMYAGA64UQ83sUSwSwgoiWAvgOwKvD3d7Um31w5QwDHpLe\ngO7qG3juIcYYM5MljghWAcghoiIi0gPYA2Br1wJEFEZEbR0PYwD4D3djHps9ICtqwdKd2/BH7TLo\nH/7bsBvOGGPMMkHgD6Cky+NS9L+jvwPAgeFuTGYjg+vZrrhyjhq/XvQuGr74BbR/2NUxxtikpxjN\njQkhbgKwAsA5fZXZtWtX5+8bNmzAhg0bepVxXe8KdYQar380Gw8t/giv3ngbnItSARcXyzeaMcas\nTGhoKEJDQy1WnyAi8yoQYjWAXUS0pePxDgBERLt7lDsfwFsA1hNRXR910WDao45UI+fBHKyMX4nc\nXCBq4d24+DIZPL/5wKzXwhhj45EQAkQkhvt8S3QNxQGYLYSYLoSwBXAdgH1dCwghlgH4AMBlfYXA\nUDivdIYmUwNDowGzZwPy11+F/odfoTsYam7VjDE26ZgdBERkBHAfgIMA0gDsIaIMIcSzQohLOoq9\nAsARwF4hRJIQ4kdztimzk8F5pTPUUWoAwA33uuKTM9+H+po7+SwixhgbIrO7hixpsF1DAFDwdAHI\nSJj54kwA7ROShgbciMWXTMes/700ks1kjDGrYg1dQ2PCdb0r1OHqzsfu7oDrzgdgOHDQrHoNzQZI\nWp6+gjE2eYzbIHBZ7YKmxCYYW42dy+ZcOh8BzZlmzUOU+0AuSt8stUQTGWNsXBi3QaBwUsBxoSOa\nYps6lwUscIUarlCnlvTzzP41hDagMabREk1kjLFxYdwGAQC4rXdDQ3hD52MhgBKnYJQfzhhWfW2l\nbdBV6NAYy0HAGJs8xnUQuK5rv7Csq0a/YDQeH14QqCPU8LjQA5JWgrZMa4kmMsaY1RvfQbDWFY0x\njZD0p8cEpHnBoIxhBkG4Gq7rXeES4oLGOD4qYIxNDuM6CGw8bGA/wx7Nic2dyxxXBsO5pD0ICp8r\nhF6lH3R9DRENcFvnBudVzmiKaxr4CYwxNgGM6yAAAJezXLr16U899wz4qdOhrdCi8JlCnLz0JIwa\nYz81tNPX6aEt0cJxiSOcQzgIGGOTx7gPAqdlTmhOPn1EMGOVD0giqP4og9sGNzjMckDaNWnduo9M\nUR9Tw+UsF8gUMriEuKAprgnWdLEdY4yNlPEfBEu7B4HCRqDIIRgVPxXDaakT5n04D5CArLuy+t2x\nn+oWAgDbKbaQu8jRmtvarYw6Sg1drW5kXghjjI2R8R8Ei5ygydB0+8ZfNyUYLclqOC11gsxGhgV7\nF0CTpkHlJ5V91nNqoPgUl1Uu3a5RaCtpw4ktJ3B89nGkXZMG1UEVSOIjBsbY+Dfug0DuKIfdNDto\nsk5PNmeYHQypQsBxiWNnmcDHA1H9v2qTdRiaDWhJb4FziHPnMucQ525nDpW+UQrfu3yxunA13Da6\nIX97PpLWJ6GtpM1UlYwxNm6M+yAAencP2S8MBrU6wDHYsXOZx4UeaIxqhL6h91lEjdGNcF7uDLm9\nvHOZc4hz5xGBvl6Pyk8qEfBwAGzcbOB/rz9WJKyA16VeSAhJQN2vZs+s3V1rK1rf/5hnUmX9enbp\nD3jlzmw08XkNzEwTMgjcp8+EHSogszv98hTOCrid4wbVr6pez1dHqOG6zrXbMucVzmg+0QxJL6H8\nvXJ4bfWCfYB953ohE5i2fRoWfLsA2fdmo/CFQsu8GIMBuSuuRclfXka19xkofft7gAetWQ95J1rw\nyIlb8ecv1uDXqXfgxzcL+WPChm1iBMGS7kHgZusMJ+TAoOp+UZjXFV6o+aGm1/Mbwhvgtt6t2zKF\niwL20+zRFN+E0n+VIvDxQJPbdlvrhuXRy1HySgnI2P9fIhFQ2fcwBUCEtLPvQlmhHjiZisM3foym\nR57GSf/NaCq08FHHJJLy3E+IXr99rJthUWlv/oES31VwqcjG2m2+2PDYCuw5++2xbhYbpyZGEHQc\nEZw6K0iXqUGrvBllhzO7lfO8zBP1f9R3m7FUk6VBS2oLXM7ufb9j51XOyH0wFy6rXOB4hmOv9afY\n+dvBzt8OzSeb+ywDAL++W4jfpt0NXR8nHiVs2gHdiQzMSvoWcxfa4vr/bIRfZRIayA3J93/Yb92s\nb5pPv8HyiDdRl9ZfCo8vtvt/hPHirYC7OwI+eQE2JxJxXtzLOPjE0bFuGhuHJkQQ2E61hZCLzvmB\nmlOa0eKmQG1E96kmbL1s4bzcGfUH6zuX5T6Si+lPTIfCWdGr3lPXEwT+zfTRQLeya1zQGNn3tBQ6\nHaD9+y7cpv8v4j/vPQVG3ENfwiX8Z3jF/IqAeadDx9XLBg6P/gVTDn81YBuYCUSYWXQUyW4bkHbf\n+2PdGotQVRsQUv0LZj+6tXOZ4xnT0fzuZ1j8yo3IDS8fw9ax8WhCBIEQottRQXNKM9pmeUKX0nuH\n63Xl6e6huv11aM1phf/9/ibrdd/kjqm3TYXbWjeT67tyPdu189aZpnz7fAbOa9uPE8tvheq/3/Za\n7/DJe6h+7FUELvXstW7p/evgpK1F0f60bsu1lVqkXp06YJfUZFYemg09KeD+2dsIDv8AbQ3j/yyv\nhLcjUe8yHQ7zpnVbPvOeTSi9+M9QX3QdWpsMY9Q6Nh5NiCAATncPtRW1Qe4oh+2KObDLNxEEl3uh\n7pc6GDVG5D6Si9lvzIbM1vTboJyrxPyP5w9q+65rXKGONB0ETU2A86tPQ3Pvo3B+8A7MTPy228Be\nZWQefBpzseKJC0w+X2EnR+aS61C8++tuy/O356P2u9o+t2tJJBG0leNvRtbiz44iN3Aj5l46DwUe\nKxH38Jdj3aRejEagoGDw5XXf/Aj1OVtNrlvxw07InZWIWPd3C7WOTQYTLghaUlrgtMQJnmuD4VnT\nOwjsA+3hEOSA9OvT4RDkAI+LPCyyfYe5DjA2G01OX/3lY0lYK4uE30v3YcaNa+Ap1SDz5+zO9QXP\nfY6U+dfB3tmmz/q9H7wBM6K/6ryITR2tRv2hegQ+Hoia73sPgFtazfc1iJ0f2+16jfFAFhYK49oN\nAAD5ow9hytdvWt2FgN89m4rweXehrnbgdum0hIW5P2LGQ5ebXC/kMsyK+gIL0v6HiIe/t3RT2QQ1\n4YKgOaUZTkucMOP82fDRlkBq6N1v73WlF1T7VZj1+iwIMez7PXcjhDDZPVRVBcz8eCek7X8HHB0h\n5DLkLLwS5W9/116ACIGhn8P1/lv6rX/hzcuggy2yPj8OMhJy7s/BrFdmwWHrFNT+UDvi8yI1HG2A\nY7AjTl52Evr6wc/oOqaIEFQcioCbNwIAlv/tfAgA8bsPj227ujAaAeUbL+AGw6f46c6fByyf8MlJ\nKGwFPDcs6rOMc5AXmj/5FsFv3YP8A1mWbK5FSQYJ1d9WI+mcJMQtjUP69enIfzIHpf85OdZNm3Qm\nTBAo5yqhq9BBHaGG4xJHuPnYIsrlQmQ9u6dXWd87fBH8VXC3C84swVT30IGnorDCPg2eO+7qXOZ0\n29Xwj9oLACj6OgoayQ4r7lreb91CJlB41g2oeesrVHxUAbmDHGFHUlC8dhNILtCc1P8ZS+ZqCG3A\nnHfmwGOLB9KvS4dkGP59oUdLxZEMaMgBczbNAND+Hlbf8BCk198Y24Z18ccHeVjbdgjN736GNT/v\nQFFe/337qo9+RNnKre234+vHvBtDkH7DizBecRU01SP72RiOsnfLcHzWcZS9VQb/+/0x7//mQT5N\ng4ZXPkLhPXmIv/RDq7l+RtJLVncUaWkTJgiEXMBxoSPqD9fDaakTAKDtlrth+8m/e5W19bGFzzYf\ni7fBZY0LGqN6XLtwdC9KL7oHsLPrXLbgz+vg0VaGimN5qH7tc+SceTPkioGPTKY/cQPmJP+CgqcK\nUG6Tjo2f3oIA+xpofVpQ+0OtxV/PKboaHbRlWjgtdcKs12YBBOQ9ljdi27OUks+OIi9wI2RdPuUr\n37gRQaoEZH6fPnYN60AEtDz7T1RtvQfuf74WtoE+OHzzJ/2Wn570I7zvNN0t1NO6z+5Cuf8qnFx9\nV++dKhGyd/+AhAW3oDq2cPgvYhhaMlpQ+HwhFny3AMsilsHnKm/kf/g1PF+9EGWXzYPd32dD84sX\n4s+4H8aWsRncr/qyChm3ZCBuaRwinCJw8pKTE3s2YiIy+wfAFgCZALIBbDex3hbAHgA5AKIBTOuj\nHjJH5t2ZFOYQRpJBIiKiJrWRimTTqfLXeLPqHSyDxkBhyjAytBg6l0Upz6P8d/f3Knto9j10/JJn\nSSXzoJRfiga9jeM2j1GY/btUbBNEVeGZFHv1bjro/yQdX3C8WzlJksioMw7/xXRR/W01pVyU0vlY\np9JRdFA01f1WZ5H69U16ilsRRw1RDb3W1dcT3X5ZLd2/NJxeCvoPfei3k9KOVpGxzUh1f9RRc2pz\nn/XGz7iKfr/ps17Lwza/QEdn3m6Rtpsj8vtKqpe5k6GskoiImo4cpzKZP52IaelVtq2N6KNt+6lO\n7kWSTj/obTRVayhBuYaK7OdS4hXPkSYtn2oOxFGmzzpKUyyi/Yv+RlWyKRT51AGLva6B5DyaQ7nb\nc4mIqKWkjhKDrqQ02yWU+FVGZ5n0J3PpiHwvZTqvppP3vktNx5KJDIa+quxOrycyDv+zX/BcAcXM\ni6HyD8tJHacmvVpPccvjqPyj8n6fZ9QZqDK2iKpPVFBdroqa69qG3Yah6th3Dn8fbs6T27cPGYBc\nANMB2ABIBjC/R5l7AbzX8fu1APb0UZdZb0bpe6UUH9J9p79v9QsUu+xus+odioTVCVQfWk9E7X+8\nFZhC2rySXuWOPvkHaYQDRTtsJEkiakxspNr9tQPWf3jKLxRn+yeqPlFBRER1MdlUBj86NjWSWrLa\ndyCSJFHm3ZmUtCGJJEky+zVl359NRbu7h1XNvhqKmRdDRq35YZN9fzZFz4qmxPWJvdr78z8zqVnu\nTA3Bq6niotspeeojdMjhMwp3CafYRbEUuzCWJKOJ12g0Up3MkzIOFlNbaRvpVLrOVeqCOlIJdyqO\nLjW77ebYM+sJSjvn3m7LspdeTR/PeZH0Xfb1cWEt9IXHfVTtEEi1e48MeTs6rUSh/4imX/wep0zc\nQMnYTgf8fqLEi05Q+i3plPJWOJXL/emnZU9T+P/KKSHWQAUFRDrdgFUPmVFrpGM+x6glu4Wy3j9M\n5YoA2j//YWqobO1WTpIkSrs9g/7wO0C/u75AibiHssXtlDD/MVJFZ5qsWzIYKe7+T6ncJpBynJZQ\nZXjWkNomSRLlP51Px4OPU1tF9514U3ITHfM+Rm1lpnfuLdXNFO+9maqFN1ULH2oQrtQGW8qzmUsR\nAdfRoS2vUMwz+6kwtIAMFvqC1pU1BMFqAAe6PN7R86gAwG8Azuz4XQ6gpo+6zHoztNXazp3wKdmh\nZVQv3Ehb22hW3YOV80gOFb5YSEREqUerSS1zJTKxM1bX6akGnvTj5R8REVHu47kUFRDV77d4o9ZI\n4c7h1FrT/S80T7mQws4Ko6J/tO+sC18spLhlcRQzN4bqfjf/W3vs4lhSx6i7LZMkiVIuTKHifxab\nVXd9RD1F+kaStkpLMfNieh1lfL3iVTq57l4yNBso/eZ0ipkfTTniSop77yRJkkTxIfFU9U1Vr3qz\nn4mj43iFIgOiKMIjgqKDoklbre1cH7HiQTq0/HGz2m6O1Cg11QlPakvP67Zcm5pNDQoP+o+4m95U\nPkGv+75CWfL5VHD2DSTVqYa8HckoUf7OfIoOiqbIqZEUf3kqpfy9gMo/Lqfq76vpxCUnKPexXGrM\nrqCTQZdSnY0P6aCgStlU2rPkJbNeoyqlmGoiMrp9O6/6XyVFLwin6NUPUbnMj47u+K3P5xu1Rsp+\nMJsybsug9D9nUdhlKXRY8TudxAN0fOb1lPL6ITrxbjglv3GEoh/5hjKVSylZuZqO7T5GP1/yAdXK\nvOjkzj2De58kifKezKPYhbGkrdKaLJO/M59ObD3R68tKQ34dnXReTaEz7qTGjKbO5YY2PeXvO0lR\nf/6Uwpc9QAme51O53J+aoaR0h+UUPudPFHrV25T0xlHKP5BJqpxaMur7DwlJkqghuoFyHsmhuoOn\n/1bMDQJBZvZ7CSGuArCZiO7ueHwTgFVE9ECXMic7ypR3PM7pCAZVj7rI3PaYcsz7CthdfiFC/nu3\nxevuqeb7GlR8WIHFvy7GoSePIuDfT2F+7TGTZT9+LA2b/joXAUE2SL0iFaqDKsz9YC6m3jzVZPmG\n8AbkPpKLlfEruy2PueBpNOY6wtPnfPjf54+CnQVYHrUcDeENKH2jFMtjlg/77Ci9So+YGTFYU7cG\nMpvuQ0qaLA0S1yQiJDUEdlPt+qihb8ZWI+KXxmPmyzPhfaU3qvdWo/gfxVgRvwJCCEgSEGV/LgJ3\nPoLqvVPgtNwJc9+fi4Q/vQX9H6E4u3YfVAfqkPe3PIScCIGQtb/GxuONSNpwHHrHOJwd/RAcZjug\nYGcB1OHB5bTaAAAgAElEQVRqLDm0BDI7Gapii2CzejmQlw+PINcBWmp5/wv5JxbqErAg5eveKyMj\nYUw+iZaiWrSW1kG5aQ2cb796yNsgIyHzjky0FbZhzjtz4LjAsdfnQFerQ/ySeAR/EQz3je7tC/V6\nNMVmQL9uI1BaBg8/exO19y/vyxi43LIVLXCEB9WhyHM59A4uoJL1kMlSUT3PDQu+/DsClnkPqV5d\ntQ4nb0pHY3Q5nAw/QCJnGGkKjHABbp6CNf++AjJ5+2s8/kESfO7bBtXCczDzv0/APWS2yTqNrUZk\n35GOhuP1kK3KhJQcA8+iRDjq6yEnAxSkh05mD5X9NNS1PQwE5EO+xhluZwXDabonDNuuQ9G8S+Fm\nczVaMlrge4cvZjwzAzYepk8HbypvQuEvqagPTQGSk+FamgrXtiq4GWrhRE2oFx5oUExFk50/Wu39\nobcNgGTjAwFviHov2HrawfcmX0y9bSocZjoAaD9rkYiGfQrkWAVBbkeZXkHwzDPPdD7esGEDNmzY\nYFb7ACDiyd/g/fZOzG+KN7uugWgrtYg7Iw5ratfgly3vwLc+HSvjBp7aIG5RHKbcNAVVX1ZhZcpK\nkzvugl0FkFolzNo9q9vy0p+ToLv8OpS7fghhI7D0yFI4LnAESYT4JfEIejkIXpd4Dev11P5Ui7L3\nyrDk9yUm1+f9LQ/6Gv2gL7zr9twdeWjLa8OCvQsAtF+0lrAyAdOfnA7vq7yREqGG6/rbUOb1KIJe\nmgnfO30hhIChRYtK9/moePlTrHxkHRJXJyLw0UD4XOMDQ7MBCcsSIBr2QnXhcmz57IbOutOuToPC\nVYF5H82DEALRM+6GxmYd1hzaBvvpQ9/ZDVddlQEa31lwPvgd3M5fOfAThkEySMi8LRO6Ch0W/bwI\ncqW8z7J1v9Uh++5srExZCRv30zuvE1M3oeqSO7Hp/64d0razP4yAx11XIfWxT3DO7otQEF+H4u/j\noS9ohO3+KVhTsQYKx77bMxAiQtWXVaj8qBL2M+3hGOyIlrQWkJEQ/Glwt7IVmWokbnsZZ6Z9iDK/\nENjcfy8U9gq0pORCys6FIq8cmqqL4EjFkMt/QWXgYkjLQ+B94Uo4B3lBbqeA3N4GOnUr1OllaAmv\nhOEnF8g8YuHafABTNXlIPecReOjOg9MyJ8x4egYKdxWi5tsaBDwSAKfFTrDxsYGtty3aStrQnNCM\npoQmtBW3tf+NywAQYGwywtBggEFtgKHB0P4t3R6ArQGwbwPJmyGMNXBsjkRQ0yEcIz32rd6MwM1n\nAACeffZZs4LAUl1Dv3V5bKpr6AC6dw1V91FXv4dFw6VrNVCxfDpl7xmdQeOYOTFUH1ZPB6bfTYl3\nvDNgeckoUZhDGOkb9RS7OJZqD5geK0hcm2i6q0eSqNRmOkXfdJzqw7t3jVV/X01xy+JM96MPQs7D\np7u6TNE36inSL5JyH8+lvCfyKOfRHKr4pGLAejW5GjrmdYy0ld0Pw2sP1NLx+ceptbiVfvc7SBE2\nX1PTiaZez4++/0s66biKJKNEtftr6fgZx0kySJR5VyalX59CauFCWdHd3ytDs4HilsZR9n3ZdPLy\nkxTmFEqR+Jyyn8gd5LvRm0pF9LdHDXTwoMkeQJN+vPlbyvY+e9jbHIhklCjtujRK3pTc7cSF/mT9\nNYvSrk/rtizp8S8p2nXzkLad8c4hqhHeFPb0oV7rCp4roKx7h9ZvP1j6Rj1F+kaSOk5tcn1tiYYO\nbPuQ4pVrKdr5fDoQdC/9vOx9Oqw8TMduSKTaqsEPwDfGN1Lc0jhK3pRMzanNlLQhiTLvzOz2N9ac\n2kyZd2ZS8uZkilsWR5H+kRS/Mp4y78mksv+Wkeqwqv3nUPuP+riaWjJbqK2ijQwaQ/9je5JEVUll\nVJd9ej8BKxgjkOP0YLEt2geLg3uU+QtODxZfhxEaLO7P4TN3UOi6nSNWf1dVe6oodlEsxduuocLP\nwgYs31rSSpFTI4mIqOLzCko6N6lXGX2TnsIcw/r8w4468yH6/exdvZZLkkRxy+Oo+tvqIb6KdnHL\n46jhWO+zebpqiGyg/GfyqfDFQip6tYiO+Rzr92weIqLS90sp/ZZ0k+1NXJdIYY5h9LvTW3TyT2+Y\nfL5Rb6Q0h+UU89g37WMFZ8ZT+q3pFD0zmlIf+pSOuFxm8nmtxa2Uui2Vyv5TRvomPR13vYeOePY+\nq2swiouJnvT7iLRyezquPIeenP45fftBI2lb+hnnMRLF2K+n7OcH13c9HI1JjRQdFE0GzSDPsiEi\nQ4uBYubFUM2PNZ3L9OoWUgl3ygsf3KB68hN7qFp4U8RLvT/zklGi6BnR1JgwcmN1Zf8to8S1vU84\nMKW1qJUifSOp5qeaAcuaYtQZqegfRRTmEEbpN6cP+4uWpYx5ELS3AVsAZKH99NAdHcueBXBJx+92\nAL7pWB8DYEYf9YzYGxX70h90wmXkvoV1JUkSJW5KpizcSvrKgc8EUh1RUeK6RCJq/4BFBUZRY3z3\nP5jaA7WUuD6xzzqy/i+c0myXmPxWWvtrLR1fcHzIH1ZdvY7CncKHfGZQ8evFdOLSE/2WSbsujco/\nNH06niZfQ1VRjVQq/Kk5se9vkJE791O2w6L2o4IDtXRUdpQaIhsoNeAC+uH6we1osz6KoTDsp/qc\n3kcd/UlN0tGHTg9QndcckpJTSP/Vt5Tuv5PC8Qv9bLeHEo6aDs+IfyVRpY0/SdoROCWnQ+VXlZR6\ndeqQn1d7oJZi5sZ0O2EhauFddHDjAIPGkkTRl71EJbJASvgo2WSR0vdLKW5Z3JDbNBSSQaLYxbFU\ntbf3yQNdGVoMFLc8jopeHfxp233R1mg7T1cfS1YRBJb6GckgaKpqoSY4kqZqdM4eSvpvDoVjH2ny\nNAOWLft3GWXcfvoc6uLXiin1mu5/yLmP5VLBswV91iHpDVQt86HM/Xm910kSxS6KJdWRoZ11UvNz\nDSWd1/voZCDGNiNFTY/q1U3VtT2RvpGkye37vTn6VgqV2c/st7/FoJcow3YRpb9+gCRJIk2BhnQl\nldQgXKkwvfe5+H0Jt/uIDi37aNDlT4TXU7jNuVS2eAtRfX37DnReDKVsSaHGX7PomNM3lCruob3n\nvU9N9d27HA4G3E5xV7446G0NR/5T+ZS/M3/Iz5MkiZLPT6bSd08fAWR+HEV5irlk7GNnV/VjOYV7\nv0mpdssoP6L3kYMkSVSwq4Cig6KpOaP/o0RLUB1SUXRQNBnbTH95kSSJ0m5Io7Qb0yxyarW1MDcI\nJsyVxQNx8lEi22Ulsj40fQaPpdWmZcHWKR7Zf8k+FXJ9as1thcMch87Hvnf5oim2CVVfV3Uuqz9c\nD/fz3PusQyjkyA26AOWf955HRwgB3zt9UfF/FUN6DarfVHA7Z+ApuHuS2ckQ9HwQ8v+Wb/K1a7Jb\n0dIs4a2t36OtjwtHG77ej4rlF/U7lYJcIVB4zd9g/McrEELAYYYDsl/4BtGel2J6sHLQ7XW5eQE8\nkwrQUDLwzX81LYSqi27H1DWz4LHvO6TeWYqc+3Iw+7XZWLR/EZwvmouVmZeh3uN6zExMRLnPEoQ9\nGwoioCi+BiFlP+CMN+6GockAY4txwO0NhyZTA+UQXv8pQgjMfGUmCp8rhKFjGuu5t6yGkAmc+Hc0\nAEDfZkTE9p9xdNEDSFOGoPDyX2GsOQM+cUcQtLb7dO6STkLWn7JQ90sdlkcvh+N8y07pYor7ee5w\nXOSIlE0pKN5dDHW0GoYmA1rzWlEfWo+8x/OgydRg3n/nWWyesYlg0gQBAKiWnIumfUdGZVu6xFTQ\nohboynSo/qq637KtOd2DQOGswMJ9C5H7YC4ajzdCX6dHa24rnEOc+63HuHIVREKcyXVTbpyCul/r\nek0YV/VVFRLPToShsfscNw1hDaj9rha+f/Ltd5t9mXLDFBg1RtT+2H3qi/x84MNzI+DVHIo/5zyM\nr7eZniHTL/kAPG64cMDtnP3mtXCryUXVL+2v2/a7r6C7+oYhtXXO4wvRIF+FhJvfGrDsvkv+jZmK\nEjhc8nckrEyC4wJHhKSGwPNiz84di52/HRb+uAStilvR+peXMOfFW3HE70ak3/wSshZcBeU0L+T8\nNQfRAdHI35kPXVUft6wbJk3G8IIAAJyXOcPjAg+UvFICoH1+pqKNt6P5+a9x5OLXUOE0G57vvwib\nmYFou/110ILF8LnBH22Heodo1l1Z0NfqsTR0KWyn2Jr1moYi+MtgBDwYAG25Fjl/zUGUTxRSNqWg\n8OlC6Kv0WPjjQsgdhn/W0oRkzuGEpX8wgl1DRESxb0RQhuPyYT1XXauj7y94j5pqB3fZ+GG/Gynx\n/o+oMaGRjnkd6/ew+PiC49SY1LvLqmZfDUX6RVLJmyWUcmGKiWd2l/FJNKXbL+1zfdp1aVTyr9NX\nORuaDRTpH0kpW1Io+fzkzrEAbaWWIv0izb4Y7VSXyamLZA7s09G/7B+jZJsXqOSRUGr4PYaqZT4U\n8VX3i9LyEuqpSTiR1Dy47p3v179BKfO3kSY1j6qFN1UUD63/XZIkivANp1IspZqsvl9z+HsnqVbm\nRQXPJVHs4lhqyey/fSVvl1DiOYmkq2+m2PN2UCvsqejnFJIMEkV4RpDqqIqy7s2iCLcIKnypcEht\n7otRb6Qw+zAyNA9+oLin1qJWivCIIHWcmvKfyacIr3AKw36K8Pk7ZX3ePpVJ58V8e6uo4VgDxcyN\n6dbV0pjQSJG+kWa1gw0eeIxg8FrqtaSGMzUVD62vvCipjmIcN5IRgmJfOTqo56TaLKGSH9oHx8r/\nr5yOzz9O+sbep6h1njraZPr0taJXiugojg7qCl6tupVa4EBN1ab73lWHVBS7OLbzD7bguQJKvTaV\njHojnbjsBKXflE5GvZGSzk2i/KeG3sfc06kzgCo+bz+d9KuAx6li4fkU5RfRuRPNuu0lirZbT6qa\n9h1GbS3Rx2f9m1L8twx6O3kpTVQjvCh73e30o/9fhtXWrHuzKMznNTrmczmF/yOSmtTd+5irizSU\noVhAGds/bh/o/r/+550hah8riXCLoLbyji8Pbe3/NkQ3UOzC2M5yzenNFOkXOax299SS00JR06PM\nridvRx6FOYZR1r1Z1JzaTOrjajrmffoLTe2BjlN2jVL7GNTC7mNQyRd0H2tgI4uDYIhi3DdT4tM/\nDLp8wleZlCufQwkbH6WIpffRoc27B3yOuk5PGtiTsfH0UUDmnZmUenVqrwGq1uL209j6IkkSFb5c\nSK3FrX2W6SpDuYyS3ze9I5CMEkUHRZM6Tk1tFW0U4RlBmvz20DC0GChhdQLFLomlpI1JFjsTovZA\n+xlLVeUGKhX+VP/TCTo25djp98FgoCz/jfS/4Gfok63fU4zibKp2DhrynDr/m/MkEUD7th8bVjtr\nfqqh+DVxFHX+U5TnuJAqxFT61edW+sX7Nvrd5So6IV9KScHXk6Q3UoRnBLWWDO7/I+2GNCp9v/sO\nMf+pfMr92+lrFySjRGGOYaRvGPy57H2+jn01lLJl4KPHgRj1xl7f5sv+W0bHg9u/0CSclUCVX1V2\nrit9p5RSt7Wf4KA6rKLomdEWmYeKDY65QTCpxggAQL38XGh+Hdw4QeYvuQi4cT1a/rody4/8E4q1\nq+FwMnbA5+UfzEGtnT9kzqcHx2b/azZaC1pR+nppt7Ktua1wmO3Qs4pOQghM3zEd9oGDu/K1dkYI\nVL+bHicQMoGpf5qKyg8rUfh0IXxv94VDUPu25Uo5Fv68EI4LHRH8VTCE3DIDaR6bPSCzlSH+/igY\nnd3QWucFtw1upwfq5HIEhn6OC3LfxTkx/8Csfz0M7/oceF69cUjb8Xv5fnwuuwVrHz9rWO102+gG\nTUorQn54GjObT8IpIRyzbjwLc+9Yh3lPXwfvj3ZjScJHaIxvgp2fHewDBvf/4X2lN2q/6z5OUre/\nDp4Xnb43tZAJKOcpock0/+5vmgwNlPOHNz7QlUwhg7zH1b9+d/rB5WwXJG9Ihr5OD59rTk/lPuWm\nKaj/ox7aCi3yn8hH0AtBfd4CllkfxVg3YLR5bjsXbg/fNqiyZU99gLo1t2PNW3cAAPwuXwXZ+wPf\nC7YuNBVGn0UI7LJMbi/Hgr0LEL8sHn73+nVe8t9zoNhsq0JgExrW5+qpt01F3II4yOxlWJW5qts6\nWy9bnPHFGZZrC9qDbNqOaWi4JQINl1wDm7CGXmciOcz2h0NDEdwcHAa84Upf1lw5BdMLP4W758Bl\nTVE4K+Ac4oyG0AZ4XeIFp2VzMG/ZnF7lVAcq4LFl8Lc39djigcw/ZUKv0sPGwwbaSi3a8trgcrZL\nt3LKYCVaMlrgcqZLHzUNjiZTA5fV5tXRnznvzMGJTScwbfu0bl8WFK4KeG/zRvo16SAdweday9/v\ng42cSRfZC25aBo/WMjTlVvVbrq1Rh0Upn2PGc3/qXBa4cTYcpSZUn+z/uYbkVGjnLuy13CHIAc4r\nnaE6cHqKJUsHwZSLQ+BXZvqIAADsA+zhscUDQc8FdZtX5pSoKMDQ/02yhszuPE84aHXw3rQV6jC1\n6VNSlcphhwDQ/tTAwIHL9cdjswdU+1X9llH9poLHhYMPArmjHO7nuaN2X23n893Pd+81gZ9jsCM0\nGRY6IhjmGUODIbeXt99M5preO3q/e/ygPqZG0MtBnRMAsvFh0gWBvaMc6V7rkfOfo/2WS3j+V1S6\nzIP/xrmdy4RMIN99JQr39r2jBQDHgpNwWNk7CADAZ5sPqveePp10oK6hoQq6ZAF89KVoKFL3WeaM\nPWfA7x6/XsuLc7SQrVmN55b/gOr+z3gdkrh/RcPF7ijKPtDD2GIc0R2VOXyu80H1/6phaDadhLoa\nHTSZGriuGdpspd5XeaP2+44g2K+Cx0W9g0QZrDQ7CIio/RoCC3QNDYfzCmcsjVgKj82DD0pmHSZd\nEABAc8i5wM/7+r0nqs2nH6Lpmjt6LW+cvwqa0L7HCXRawoyaeAReutTkeq8rvaD6TQWjpv1iIk2O\nBso5lvvDVdgrkO+8BAXfJvRZpq8LaZJeOoAAt2Y8UvQQvp29AzHHLHNooPtyL7SbpkFXrYPreler\nvZDHfro93Na7oeoL00d89Qfr4bbRbch93x4Xe6AhtAF6lR71f9Sb7FqyRBDoqnSAHLD1Hr1z9nty\nW+tmtf+/rG+TMgimP34NFDmZyJ++EZq4tF7rKxPKMLc2Ektf6D0HvMP6EDhn9B0ESR8ng2xt4XX2\nXJPrbb1tO7uHSCK05bdZ9IgAAOpnhUB9uP+jFlOUP32Nplvvh1tuPLYFxUN37mb89KF590JubZGw\nNO9bzHziasx5Zw78/tz7SMSa+N/vj7J3yk6dxdaN6jfVkMYHTrFxs4HrWlcUPFUA+5n2sPPtfe8G\nh9kOaCtpg7Ft+FcbazI0cAwe+at32cQzKYNg/kZfeBfGIcJ3G1pXb0DmxY+CWk5/G8t68jOkzN0G\nR5/ef1TTrl6FoLo4kGT6aKL+ox9QFnJFv/3dp7qHtGVaKNwUvc7OMJfirBDYnxhaEBScaMKZDb9h\nzhNXA97e8E78HcHXLkHQny9A3OHGYbcl4V9RaFV6wuPs+fC+whse51t3t4Hbxvbxi4ajDd2Wk0RQ\n/a6C54XDG432utIL5e+XdztbqCuZjQwOQQ5ozWkdVv0AxrRbiI1vkzIIAMA3QI5bj/8VeT+lIS+q\nEmVTlqPm9wSQRJhx5CO4PdK7WwgAfJf7olUoURGZ32sdETA9+Uf43H15v9s+1T3UcqLFsmcMnWrj\nZSEIrBxaEKS++BOKAtdBMaVjRyWXw/uz1+BywVnQX3QZirOGt4PSfPINas7ZNqznjgUhBPzv80fZ\nv8q6LW9KaIKNl82wb2DjtdULEDA5PnCKcr553UMjPVDMJq5JGwSnrLrEB5uqvkTUBbsgLrwQSUtv\ng1bYYfGdq0yWFwIo8glB6fe9u4ey9ufB21iFoOtX97vNU91DpW+XjkgQzDh/NhyNatSmD37E1/3A\nV5Df3GOOHiEw4+d/wWOBL/JCrkOjauAxA31BKQpCtiHfKwRVCj+clf0Jgp64bqgvYUxNuWkKGsIb\n0FbUPiMeEaH2+9phdQudYutti+XRy/s9PdTccYLhTjbH2KQPAgCwtQWu+f46lP8cD01BFQoufbDf\n0980C1ZBF9n7G3fJOz8i+4ytEIqBu3p8tvmg/mC9xccHgPaLgfLcV6Lw28HdmjMvpgaLmyMx7/HL\nTFQmw7zoT+HtpkPK3G0oiyzst66sm55HcoErUu56B82HY+GkVWHK2t7n41szhZMCU2+ZirJ3ylDz\nYw2Szk5C9d5qTP2T6XtJD5bLKpd+P1fKYPMuKrPUxWRs8uEg6GLxxdOwpvE3XLD3rn7LOZ+3Cm7Z\nvY8IvI/9CMcb+u8WOsXrSi9ABoueMdSVem4IWkIH1z2U9eK3yJl9EeSuTibXCztbzE/9DmLRQjis\nW4HU9X+BVFLWq5w2twQB0XsRtOcfuOLlMzHrnAAIm/F5zaLfX/1Q8noJil8sRsCjATgz60w4LTT9\n/ljKqYvKhsPQZIBepR/Vey+ziYODoAchBr6uaea2FZihTgbpT3eVVJ6oRlDLSZxx37mD2o6tty2m\nbZ8G5zP7n1p6uJTnr4FHwh+DKutz5Gso77i+3zIKFyXWHn0eNRFZSM51QtPMxWj6Papbmdy7X8Gh\n6Xdg6flew263tVDOVmJ10Wosj10On6t9LDblRr/bnK9Ea3YryNj//StM0WRpoJyr5Au52LBwEAyD\n92xXVCoCUfLb6VNPs/75M9IDN8PGefDfyGa+NHPQc9YM1bLtF2BqUzZKjub2Wy77UDFmtaZh3oNb\nBlXvvDVeuL7kFXyx6TPotm6DVFoOADCUVMA/7EtMe+tRs9tuLewD7Ef1nHiFkwI2XjadYxODZWwz\nQh2h5m4hNmwcBMNU6heC+lf/C6mqBgDg8PsPoMsG1y00GuycbZG2+AbkPfNZv+UKX92LrDOugMx+\n8BchyeXA3T9djO+n3Ivy1VcCWi2y73kNh3xvxqrLzOtHn+wGO2BMRMj+SzZiZsUg0j0SFf+pgPe1\n3qPQQjYRjc8OXCsgf/IJlDz+PGb4zkGx72rMr46C9NiXY92sbnx33ArlTVdAMuyCTGE686dG7IX8\nhWeHXLeNDXBp9N+RODsJOP9W+EX9gbovT5jb5Env1DiB58X9X69Qs7cG6kg1Fv+2GPZB9n3+/zI2\nGMLUFZRjRQhB1tSewShIbUH2az9Drlbh/O//MtbN6Y4I2cqlaHr+Tax4rPe0znlHi+F+3jK4tVZC\nZtd7ArrBiD3cBKdNq5HpvR5XVL5vzrxxDED5v8vRGNuI+R/O77OMscWI2OBYBH8RDLf1Q7+nNJt4\nhBAgomH/9fERgZmCFjoi6GMrPU9eCFRdeBuk9z4FTARB/j+/g9O8rThrmCEAAKvOc8avX0Zj+nQ7\nDgELUAYrUflpZb9lil4ugutaVw4BZjF8RDDBqTKrIQueByougXtg99MfU5zOBnY+hSU7Br5JPBsd\nuhodYufGYo1qjcmB6ta8ViScmYCQlBDY+fees4hNTuYeEZjVsSiEcBdCHBRCZAkhfhdC9JqfVwix\nRAgRJYQ4KYRIFkJcY8422dB4zPdBru86JO38rtvyoshSTNNkYsED541Ry5gptt62EAqBptgmk+tz\nH85F4GOBHALMoswdYdoB4BARzQNwBMATJsq0ALiZiBYBuBDAm0KIkbuFEutFfsdtcP3+426zbue9\n8h0yZl8GhXLspixmps15fw5OXnISVV+fng5bW6lF+k3paM1rReDDZt6Bh7EezOoaEkJkAjiHiKqE\nEFMBhBJR36Nc7c9JBnAVEeWZWMddQyPA2KpDqcdinJi/DRfFPQe5QiDZZR2kx3dg+VMXj3XzmAnN\nKc1IvSIV3ld5w266HYqeLcLUO6ZixlMzLD5bLRv/zO0aMjcIVETk0ddjE+VXAfiYiBb0sZ6DYIQ0\nF9SgdOnFKPNcgtmfPgWXc5bCqakSNo58RGCt9HV6ZNyaAalNwpy358DxDL7XADNtxM8aEkL8AWBK\n10UACMBOE8X73IsLIXwBfAbg5v62t2vXrs7fN2zYgA0bNgzURDYITkHemJF3BJoFV8Nmw9lIn3kp\n1nAIWDUbTxss/mXxWDeDWaHQ0FCEhoZarD5zjwgyAGzo0jV0lIiCTZRzBhAK4AUi+qGf+viIYIQZ\nNDpErfsbPO+9FgvuPGusm8MYs4Cx7hraDUBFRLuFENsBuBPRjh5lbAD8BuAnInp7gPo4CBhjbIjG\nOgg8AHwDIBBAEYBriKhBCLECwD1EdLcQ4kYAHwFIw+lupduIqNd8BBwEjDE2dGMaBJbGQcAYY0M3\npheUMcYYG/84CBhjbJLjIGCMsUmOg4AxxiY5DgLGGJvkOAgYY2yS4yBgjLFJjoOAMcYmOQ4Cxhib\n5DgIGGNskuMgYIyxSY6DgDHGJjkOAsYYm+Q4CBhjbJLjIGCMsUmOg4AxxiY5DgLGGJvkOAgYY2yS\n4yBgjLFJjoOAMcYmOQ4Cxhib5DgIGGNskjMrCIQQ7kKIg0KILCHE70II137KOgshSoQQb5uzTcYY\nY5Zl7hHBDgCHiGgegCMAnuin7PMAwszc3qQRGho61k2wGvxenMbvxWn8XliOuUGwFcCnHb9/CuBy\nU4WEECsA+AA4aOb2Jg3+kJ/G78Vp/F6cxu+F5ZgbBD5EVAUARFSJ9p19N0IIAeCfAB4DIMzcHmOM\nMQtTDFRACPEHgCldFwEgADtNFCcTy/4C4FciKm/PBA4DxhizJoLI1L57kE8WIgPABiKqEkJMBXCU\niIJ7lPkCwFoAEgBnADYA3iOiv5uob/iNYYyxSYyIhv0l29wg2A1ARUS7hRDbAbgT0Y5+yt8KYAUR\nPTDsjTLGGLMoc8cIdgPYJITIAnAegH8A7YPDQoj/mNs4xhhjI8+sIwLGGGPjn9VcWSyE2CKEyBRC\nZH25m8YAAAPjSURBVHd0M00aQogAIcQRIUSaEOKkEOKBjuWDvmBvohFCyIQQiUKIfR2PZwghYjo+\nH18LIQY80WEiEEK4CiH2CiEyOj4fZ07Wz4UQ4mEhRKoQ4oQQ4kshhO1k+VwIIT4UQlQJIU50Wdbn\n50AI8bYQIkcIkSyEWDpQ/VYRBEIIGYB3AGwGsADA9UKI+WPbqlFlAPAIES0AcBaAv3a8/qFcsDfR\nPAggvcvj3QBeI6K5ABoA3DEmrRp9bwHY33ESxhIAmZiEnwshhB+A+wEsJ6LFaD/j8XpMns/Fx2jf\nP3Zl8nMghLgQwCwimgPgHgAfDFS5VQQBgFUAcoioiIj0APag/WK1SYGIKokoueP3ZgAZAAIwyAv2\nJhohRACAiwD8X5fF5wL4ruP3TwFcMdrtGm1CCBcA64joYwAgIgMRqTFJPxcA5AAcO771OwAoB7AR\nk+BzQUTHANT3WNzzc7C1y/LPOp53HICrEGIK+mEtQeAPoKTL49KOZZOOEGIGgKUAYgBMGeiCvQnq\nDQCPo+O6FCGEJ4B6IpI61pcC8Bujto2mIAC1QoiPO7rJ/iOEUGISfi6IqBzAawCKAZQBUANIBNAw\nCT8Xp/S8oPfUzr7n/rQMA+xPrSUIGAAhhBOAb/+/nbtnjSKKozD+HFERbCK2mqCIrZAqkEJBi4CQ\nTrEImoC9hV2atLYBLbURFHxBsx9AwTKKFqLlgsQUayE2qUSOxZ3FNWgEi51h7/lVO5fZ5e5w2P99\nmVngRjMz2L2TP/E7+5IuAoNmhjR6X3SNDyLuB2aBO7ZngR3KckCNuZiijHRnKD/2h4GFVjvVPf+d\ng64Ugm1geuT4WNNWjWa6+wS4b3ujaR4Mp3TNA3tf2urfGM0Di5L6wEPKktA6ZXo7zGst+fgMbNl+\n0xw/pRSGGnNxAejb/mr7B/CMkpWpCnMx9LccbAPHR87753XpSiF4DZySNCPpIHAF6LXcp3G7B3y0\nvT7S1gOWm9fXgI3db5o0tldtT9s+ScnBC9tLwEvgUnNaLddiAGxJOt00nQc+UGEuKEtCc5IONf9f\nNrwWNeVC/D4zHs3BMr++ew+4CiBpjrJ8Ntjzg7vyHIGkBcrIbx9w1/atlrs0NpLmgVfAe8r0zsAq\nsAk8olT3T8Bl29/a6ue4SToL3LS9KOkE5SaCI8A7YKm5sWCiSTpD2TQ/APSBFcqmaXW5kLRGGRx8\np2TgOmW0O/G5kPQAOAccBQbAGvAceMwfciDpNmXpbAdYsf12z8/vSiGIiIh2dGVpKCIiWpJCEBFR\nuRSCiIjKpRBERFQuhSAionIpBBERlUshiIioXApBRETlfgKdNM3WuCnXyAAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f61c809dbd0>"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Get and print correlation time\n",
      "print r\"k_{p^2/2.}\"+\" = {0}\".format(kp2)\n",
      "print r\"k_{x^2/2.}\"+\" = {0}\".format(kx2)\n",
      "print r\"k_{E_v}\"+\" = {0}\".format(kE)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "k_{p^2/2.} = 2.0\n",
        "k_{x^2/2.} = 2.0\n",
        "k_{E_v} = 2.0\n"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Print out averages and errors with autocorrelation included\n",
      "print \"p^2 = {0} +/- {1}\".format(np.mean(p2)/2.,np.sqrt(np.mean((p2/2.-np.mean(p2/2.))**2)*kp2/N))\n",
      "print \"x^2 = {0} +/- {1}\".format(np.mean(x2)/2.,np.sqrt(np.mean((x2/2.-np.mean(x2/2.))**2)*kx2/N))\n",
      "print \"E_v = {0} +/- {1}\".format(np.mean(E),np.sqrt(np.mean((E-np.mean(E))**2)*kE/N))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "p^2 = 0.240048056082 +/- 0.0550919868639\n",
        "x^2 = 0.261419866594 +/- 0.0553305863695\n",
        "E_v = 0.501467922676 +/- 0.000918372564392\n"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 18
    }
   ],
   "metadata": {}
  }
 ]
}