CTHYB.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Horrible CT-HYB\n",
    "# Authors: Markus Wallerberger, Patrik Gunacker\n",
    "# Permission to use for anything but serious science."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from __future__ import division\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as pl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "u = 2.00\n",
    "mu = 1.00\n",
    "beta = 10.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def hyb_pos(cspin, aspin, tau):\n",
    "    global beta\n",
    "    cspin, aspin, tau = np.broadcast_arrays(cspin, aspin, tau)\n",
    "    return np.where(cspin == aspin, \n",
    "                    np.exp(-tau) + np.exp(-beta+tau), \n",
    "                    0\n",
    "                   )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x110e005d0>]"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgkAAAFkCAYAAACq4KjhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xl4lOXZ/vHvxS4ooKKgFQUFBTckcResVsUKri0uUbTV\nbhRbNX1bl9rWqlWrvq5Vq1b7U4umUnvUqqHibqMIasLiiktBEAtF4A27LLl/f1yZJsQJZCaTuWc5\nP8cxR8iTZ2YuHsLMOfdqIQREREREmmoXuwARERHJTQoJIiIikpRCgoiIiCSlkCAiIiJJKSSIiIhI\nUgoJIiIikpRCgoiIiCSlkCAiIiJJKSSIiIhIUgoJIiIiklTKIcHMhpvZE2Y238zqzOzEzZx/ipk9\nY2b/MbNaM5tsZiPSL1lERESyIZ2WhG7AdGAc0JKNHw4HngGOA0qAF4EnzWxIGs8tIiIiWWKt2eDJ\nzOqAk0MIT6R4v7eBP4cQfpP2k4uIiEibyvqYBDMzYCtgSbafW0RERFquQ4Tn/BneZTGhuRPMbFvg\nWGAOsCY7ZYmIiBSELkA/YFIIYXFrHiirIcHMzgR+CZwYQvh8E6ceCzycnapEREQK0lnAI615gKyF\nBDM7A7gXGB1CeHEzp88BGD9+PIMHD27r0qReeXk5t9xyS+wyioquefbpmmefrnl2vffee4wZMwbq\n30tbIyshwczKgPuA00MIT7fgLmsABg8eTElJSZvWJg169Oih651luubZp2uefbrm0bS6uz7lkGBm\n3YABgNUf2rV+OuOSEMI8M7sO2DGE8K36888EHgAuAN4ws97191sdQljW2r+AiIiItI10ZjfsD0wD\nqvF1Em4CaoAr63/eB+jb6PzvAe2BO4HPGt1uTa9kERERyYaUWxJCCC+ziXARQji3yfdHplGXiIiI\nRKa9G+S/ysrKYpdQdHTNs0/XPPt0zfNXq1ZcbCtmVgJUV1dXa7CLiIhICmpqaigtLQUoDSHUtOax\n1JIgIiIiSSkkiIiISFIKCSIiIpKUQoKIiIgkpZAgIiIiSSkkiIiISFIKCSIiIpKUQoKIiIgkpZAg\nIiIiSSkkiIiISFIKCSIiIpKUQoKIiIgkpZAgIiIiSSkkiIiISFIKCSIiIpKUQoKIiIgkpZAgIiIi\nSSkkiIiISFIKCSIiIpKUQoKIiIgkpZAgIiIiSSkkiIiISFIKCSIiIpKUQoKIiIgkpZAgIiIiSSkk\niIiISFIKCSIiIpKUQoKIiIgkpZAgIiIiSSkkiIiISFIKCSIiIpKUQoKIiIgkpZAgIiIiSSkkiIiI\nSFIKCSIiIpKUQoKIiIgkpZAgIiIiSSkkiIiISFIKCSIiIpJUyiHBzIab2RNmNt/M6szsxBbc5wgz\nqzazNWb2gZl9K71yRUREJFvSaUnoBkwHxgFhcyebWT/gKeB5YAhwG3CfmR2TxnOLiIhIlnRI9Q4h\nhKeBpwHMzFpwlx8C/wohXFz//SwzGwaUA8+m+vwiIiKSHdkYk3Aw8FyTY5OAQ7Lw3CIiIpKmbISE\nPsDCJscWAt3NrHMWnl9ERCTvvfMOjBsHixdn7zlT7m7IpvLycnr06LHRsbKyMsrKyiJVJCIiEkdV\nFfzhD3DzzQ3HKioqqKio2Oi82trajD1nNkLCAqB3k2O9gWUhhC82dcdbbrmFkpKSNitMREQkX9TU\nwF57QZcuDceSfXCuqamhtLQ0I8+Zje6G14CjmhwbUX9cREREWqCmBrL9uTmddRK6mdkQM9uv/tCu\n9d/3rf/5dWb2YKO73F1/zvVmtoeZjQNGAzezGcuXp1qdiIhI4Vm7Ft56Kw9CArA/MA2oxtdJuAmo\nAa6s/3kfoG/i5BDCHGAUcDS+vkI58J0QQtMZD1/ywQdpVCciIlJg3n3Xg0K2Q0I66yS8zCbCRQjh\n3CTH/gmk3EHy/vup3kNERKTw1NSAGQwZkt3nzem9G957L3YFIiIi8dXUwB57QLdu2X3enA4JakkQ\nERGJM2gRcjwkzJkDq1bFrkJERCSeDRtgxgyFhC8JAWbOjF2FiIhIPB984B+YFRKa6NDBm1hERESK\nVeJ9cOjQ7D93ToeE3XZTSBARkeJWUwO77go9e2b/uXM6JAwapJAgIiLFLdagRciDkPD22/DFJnd4\nEBERKUx1dQoJzRo0CNat8+0xRUREis3s2bBsmUJCUgMHQrt26nIQEZHiFHPQIuR4SNhiC41LEBGR\n4lVTAzvtBNtvH+f5czokgDexKCSIiEgxijkeAfIkJMyYAevXx65EREQke0JQSNiskhJYs0b7OIiI\nSHH59FP4/HOFhE3abz//qi4HEREpJon3PYWETejRAwYMUEgQEZHiUlPjAxZ33DFeDTkfEsBT1LRp\nsasQERHJnmnTfOqjWbwa8iok1NXFrkRERCQ7Yg9ahDwKCcuXw8cfx65ERESk7S1cCPPnKyS0SGKl\nKY1LEBGRYpDoYldIaIFevaBvX4UEEREpDjU1PnC/f/+4deRFSACtvCgiIsUjMR4h5qBFyMOQEELs\nSkRERNpWLgxahDwLCUuWwNy5sSsRERFpO0uX+hbRCgkpSFwsdTmIiEghy5VBi5BHIWGHHaB3b4UE\nEREpbDU10K0bDBwYu5I8CglmGrwoIiKFr6bG9y1q3z52JXkUEkAhQURECl+uDFqEPAwJCxbAv/8d\nuxIREZHMW7ECPvhAISEtiYumzZ5ERKQQzZjhU/0VEtKwyy6w9dbqchARkcJUUwOdO8PgwbErcXkV\nEjR4UURECllNDeyzD3TsGLsSl1chARQSRESkcOXSoEXIw5AwdCh88gksXhy7EhERkcxZswbeeUch\noVU0eFFERArRW2/Bhg0KCa0ycCBsuaW6HEREpLDU1PgCSvvsE7uSBnkXEtq185WoFBJERKSQ1NTA\nXntBly6xK2mQdyEBNHhRREQKT64NWoQ8DgkffgjLlsWuREREpPXWrYOZMxUSMiJxEadPj1uHiIhI\nJrz7LqxdWyAhwczON7PZZrbazKaY2QGbOf8sM5tuZivN7DMzu9/MtkmvZF+JqksXdTmIiEhhqKnx\nBQOHDIldycZSDglmdjpwE3AFMBSYAUwys17NnH8Y8CDwB2BPYDRwIHBvmjXToQPsuy9UV6f7CCIi\nIrmjuhp2391n7+WSdFoSyoF7QggPhRDeB8YCq4Dzmjn/YGB2COHOEMInIYTJwD14UEjbwQfDq6+2\n5hFERERyw6uvwiGHxK7iy1IKCWbWESgFnk8cCyEE4Dmgub/ea0BfMzuu/jF6A6cClekUnDB8OMye\nDfPnt+ZRRERE4qqt9d0fhw+PXcmXpdqS0AtoDyxscnwh0CfZHepbDsYAj5rZWuDfwFLgRyk+90YS\nF7OqqjWPIiIiEtfkyb49dCGEhJSZ2Z7AbcCvgRLgWKA/3uWQtt69vf9GIUFERPJZVZW/pw0YELuS\nL+uQ4vmfAxuA3k2O9wYWNHOfS4FXQwg313//tpmNA6rM7PIQQtNWif8qLy+nR48eGx0rKyujrKwM\n8NSlkCAiIvmsqgoOP9xnN6SqoqKCioqKjY7V1tZmqLIUQ0IIYZ2ZVQNHAU8AmJnVf397M3frCqxt\ncqwOCMAmL8ktt9xCySYmjQ4fDvffD0uWwDZpT6gUERGJY80aeP11+N//Te/+jT84J9TU1FBaWpqB\n6tLrbrgZ+J6ZnWNmg4C78SDwAICZXWdmDzY6/0ngm2Y21sz610+JvA2YGkJorvWhRRL9N5rlICIi\n+ej1130RpVwcjwBphIQQwgTgp8BVwDRgX+DYEMKi+lP6AH0bnf8g8BPgfOAt4FHgPeCbraoc6N8f\ndtxRXQ4iIpKfqqqge/fc2vmxsVTHJAAQQrgLuKuZn52b5NidwJ3pPNemmHk/jkKCiIjko6oqOOww\n3yI6F+Xl3g2NDR8Ob74JK1fGrkRERKTl1q/37vLDD49dSfMKIiSsXw9Tp8auREREpOVmzIAVK3J3\nPAIUQEjYay/Yemt1OYiISH6pqoLOnWH//WNX0ry8Dwnt2sGwYQoJIiKSX6qq4KCDPCjkqrwPCeBN\nNa+9BuvWxa5ERERk80JoWEQplxVMSFi1yvfjFhERyXWzZsGiRbk9HgEKJCSUlMAWW6jLQURE8kNV\nlXeX5+L20I0VREjo1MkvtEKCiIjkg6oqGDoUttoqdiWbVhAhARo2e6qri12JiIjIpv3zn7nf1QAF\nFhKWLoV3341diYiISPPmzYNPPsn9QYtQQCHh4IOhQwd1OYiISG5LvE8NGxa3jpYomJDQrRuUliok\niIhIbquqgkGDYLvtYleyeQUTEsC7HP75T59/KiIikovyZTwCFGBImD8f5syJXYmIiMiXLV7sY+cU\nEiI47DD/qi4HERHJRa+84l/zYdAiFFhI2HZb2HtvhQQREclNVVXQty/sskvsSlqmoEICNKyXICIi\nkmuqqvKnqwEKNCTMmgULF8auREREpMGKFVBdrZAQVeLiJ/p9REREcsGUKbBhg0JCVDvtBP36qctB\nRERyS1WVj50bPDh2JS1XcCEBfNSoQoKIiOSSqipfZbFdHr3z5lGpLTd8OEyfDsuWxa5EREQE1q6F\n117Lr64GKOCQUFcHkyfHrkRERMQHLK5Zo5CQE3bfHbbfXl0OIiKSG6qqfI+hoUNjV5KaggwJZlov\nQUREckdVFRxyCHTsGLuS1BRkSAAPCVOnevOOiIhILHV1Pi0/37oaoMBDwtq18MYbsSsREZFi9vbb\n8H//p5CQU4YMga22UpeDiIjEVVXl3QwHHRS7ktQVbEho3953hVRIEBGRmKqqYP/9oWvX2JWkrmBD\nAnjTzquvwrp1sSsREZFiFAK8/HJ+djVAgYeEESNg+XIPCiIiItk2fTosWODvR/mooENCSQn07g2V\nlbErERGRYlRZ6ePj1JKQg9q1g5EjFRJERCSOyko45hjo1Cl2Jekp6JAAMGoUvPcezJ4duxIRESkm\nixb5ej2jRsWuJH0FHxKOOcannqg1QUREsunpp33g4siRsStJX8GHhO7dfetohQQREcmmykqf+tin\nT+xK0lfwIQG8qefFF2HlytiViIhIMVi3DiZNyu+uBiiikPDFF/D887ErERGRYjB5si/FrJCQB3bf\nHQYMUJeDiIhkR2WlT8EvLY1dSesURUgAT3MTJ/ogEhERkbZUWQnHHedT8fNZnpffcqNGwaefwsyZ\nsSsREZFCNmcOvPtu/nc1QJohwczON7PZZrbazKaY2QGbOb+TmV1jZnPMbI2Z/cvMvp1WxWk6/HDo\n1k1dDiIi0rYqK6FDB5+Cn+9SDglmdjpwE3AFMBSYAUwys16buNtfgCOBc4HdgTJgVsrVtkLnzv4P\nppAgIiJtqbLSl2Hu0SN2Ja2XTktCOXBPCOGhEML7wFhgFXBespPN7OvAcGBkCOHFEMLcEMLUEMJr\naVedpuOPhylTYPHibD+ziIgUg1WrfMr98cfHriQzUgoJZtYRKAX+O5kwhBCA54BDmrnbCcCbwCVm\n9qmZzTKzG82sS5o1p23kSKir81WwREREMu2FF2DNmsIYjwCptyT0AtoDC5scXwg0t6bUrnhLwl7A\nycCFwGjgzhSfu9V22MF3hlSXg4iItIXKSthtN596Xwg6ZOE52gF1wJkhhBUAZvYT4C9mNi6E8EVz\ndywvL6dHk06dsrIyysrK0i5m1Ci44w5Yv94HloiIiGRCCB4STjkFzLLznBUVFVRUVGx0rLa2NmOP\nbyGFhQPquxtWAd8MITzR6PgDQI8QwilJ7vMAcGgIYfdGxwYB7wC7hxA+TnKfEqC6urqakpKSlv9t\nWmDqVDj4YKiqgmHDMvrQIiJSxN56C/bd15djHjEiXh01NTWU+ipOpSGEmtY8VkrdDSGEdUA1cFTi\nmJlZ/feTm7nbq8COZta10bE98NaFT1OqNgMOOAC22w6eeirbzywiIoXsqad8qv1Xvxq7ksxJZ3bD\nzcD3zOyc+haBu4GuwAMAZnadmT3Y6PxHgMXA/zOzwWZ2OHADcP+muhraSrt2vgqWxiWIiEgmVVbC\n0Uf7lPtCkXJICCFMAH4KXAVMA/YFjg0hLKo/pQ/Qt9H5K4FjgJ7AG8CfgL/jAxijGDUK3n4b5s6N\nVYGIiBSSxYvhtdcKZ1ZDQlpD90IIdwF3NfOzc5Mc+wA4Np3nagsjRkD79p76fvjD2NWIiEi+mzTJ\np9iPHBm7kswqmr0bGuvZ01fDUpeDiIhkQmUlDB0KX/lK7EoyqyhDAniT0AsvwOrVsSsREZF8tmGD\nL9JXaF0NUOQhYfVqXz5TREQkXVOmwJIlCgkFZdAg6N9fXQ4iItI6lZXQq5dPsS80RRsSzDz1VVb6\nKlkiIiLpqKz0qfXt28euJPOKNiSAh4RPPoF3341diYiI5KN582DmzMLsaoAiDwlHHAFdu6rLQURE\n0jNxorcgHJszk/wzq6hDQpcucNRRWqJZRETS89RTcNhhPrW+EBV1SABvIpo8GZYujV2JiIjkk9Wr\n4fnnC7erARQSGDXK57g++WTsSkREJJ8884wHhRNOiF1J2yn6kLDTTr5j18MPx65ERETyycMPw377\nweDBsStpO0UfEgDGjIHnnoN//zt2JSIikg9qa+GJJ/z9o5ApJACjR0OHDlBREbsSERHJB489BmvX\nQllZ7EralkICPir1hBNg/PjYlYiISD4YP95nx+24Y+xK2pZCQr0xY2DaNHjnndiViIhILps7F156\nqfC7GkAh4b+OOw623loDGEVEZNMqKmCLLeCUU2JX0vYUEup17gynneYhoa4udjUiIpKLQoA//QlO\nOgm6d49dTdtTSGhkzBhvRnrlldiViIhILpo507uli6GrARQSNnLoodCvnwYwiohIcuPHw3bbwYgR\nsSvJDoWERtq1g7POggkTYM2a2NWIiEgu2bABHnkEzjgDOnaMXU12KCQ0cdZZvkjGxImxKxERkVzy\n0kvw2WfF09UACglfMngwlJaqy0FERDY2fjwMHAgHHBC7kuxRSEhizBjf/nPJktiViIhILli1yldZ\nHDMGzGJXkz0KCUmccYb3Pf3lL7ErERGRXPDEE7BihXdJFxOFhCT69PGRq+pyEBER8PeDQw6B3XaL\nXUl2KSQ0Y8wYXy9h9uzYlYiISEyLFsHTTxfXgMUEhYRmnHwydOvm011ERKR4Pfqoj0M47bTYlWSf\nQkIzunXzdbnHj/dlOEVEpDiNH+/7+/TqFbuS7FNI2IQxY+D996GmJnYlIiISw4cfwtSpcPbZsSuJ\nQyFhE446Cnr31gBGEZFi9fDDvpHT8cfHriQOhYRN6NABysp8W9D162NXIyIi2RSCf0gcPdq3hi5G\nCgmbMWYMLFwIzz8fuxIREcmmqVPh44+Lc1ZDgkLCZpSUwKBBvn+4iIgUjz/9CXbaCb761diVxKOQ\nsBlmniL/9jdfbUtERArf2rU+9fHMM32H4GJVxH/1ljvzTF+3+29/i12JiIhkw6RJsHhxcXc1gEJC\ni/Tv781N994buxIREcmGe+6B/faDffaJXUlcCgktdMEFvkzzm2/GrkRERNrSrFlQWQkXXhi7kvgU\nElropJO8ReGWW2JXIiIibem223yNnLKy2JXEp5DQQu3be2vChAkwf37sakREpC0sWQIPPgjjxkHn\nzrGriU8hIQXnnecLatxxR+xKRESkLdx7L2zYAGPHxq4kNygkpKB7d/jud31Ay8qVsasREZFMWrfO\nPwSOGQPbbx+7mtyQVkgws/PNbLaZrTazKWZ2QAvvd5iZrTOzvN0y6YILoLYWHnoodiUiIpJJf/mL\ndydfdFHsSnJHyiHBzE4HbgKuAIYCM4BJZrbJTTTNrAfwIPBcGnXmjH79fAvpW2+FurrY1YiISCaE\n4APTjzkG9t47djW5I52WhHLgnhDCQyGE94GxwCrgvM3c727gYWBKGs+ZU8rL4YMP4B//iF2JiIhk\nwquv+hT38vLYleSWlEKCmXUESoH/bncUQgh468Ahm7jfuUB/4Mr0yswthx4KBxyg6ZAiIoXillt8\nn55jj41dSW5JtSWhF9AeWNjk+EKgT7I7mNlA4FrgrBBCQTTQm3nafP55mDkzdjUiItIas2fD44/7\nWIRi3qchmQ5t+eBm1g7vYrgihPBx4nBL719eXk6PHj02OlZWVkZZDqxwMXo0XHyxj0344x9jVyMi\nIum6/Xbo2RPOPjt2JamrqKigoqJio2O1tbUZe3zz3oIWnuzdDauAb4YQnmh0/AGgRwjhlCbn9wCW\nAutpCAft6v+8HhgRQngpyfOUANXV1dWUlJSk8vfJquuvh1/9CubO9dW5REQkvyxb5ttB//jHcM01\nsavJjJqaGkpLSwFKQwitmk2YUsNKCGEdUA0clThmZlb//eQkd1kG7A3sBwypv90NvF//56lpVZ0j\nvv996NABfv/72JWIiEg67r8f1qyB88+PXUluSqf35Wbge2Z2jpkNwt/0uwIPAJjZdWb2IPigxhDC\nu41vwH+ANSGE90IIqzPz14hj663h29+Gu+7yXzIREckfGzZ4V8Ppp8OOO8auJjelHBJCCBOAnwJX\nAdOAfYFjQwiL6k/pA/TNWIU57sILYdEieOSR2JWIiEgqHn8c5szRtMdNSWlMQrbky5iEhBNOgE8+\ngRkzfOaDiIjkvmHDfPO+l1+OXUlmRRuTIMmVl8Nbb/mUSBERyX1vvOELKKkVYdMUEjLgyCNh3319\nOqSIiOS+W2+FXXf1lmBpnkJCBiQWV6qshFmzYlcjIiKbMn8+TJjgY8rat49dTW5TSMiQsjJfK+G2\n22JXIiIim3LHHdC1K5x7buxKcp9CQoZ07gw/+pGvvjhvXuxqREQkmc8/hzvvhB/8ALbaKnY1uU8h\nIYMuvBC6d/dVGEVEJPf85jfeRXzxxbEryQ8KCRm01VYeEB580Gc7iIhI7pg92xe/u+QS6NUrdjX5\nQSEhw77/fR8xe9llsSsREZHGfvELDwcXXRS7kvyhkJBhnTrBtdf6TIdCW6BDRCRf1dT4yrhXXumD\nFqVlFBLawOjRsP/+3ueVgwtaiogUnUsvhUGDNKMhVQoJbaBdO7jhBnj9dfjrX2NXIyJS3J591m/X\nXec790rLKSS0kSOPhK9/HX7+c1i3LnY1IiLFqa7OByoeeiicdFLsavKPQkIbuv56+OgjuO++2JWI\niBSnP/8Zpk3z1l1twJc6hYQ2tO++cPbZPlBmxYrY1YiIFJcvvoDLL/cWhMMOi11NflJIaGNXXQVL\nl8LNN8euRESkuNx9N8yd62MRJD0KCW1sl13gxz+GG2+E//wndjUiIsWhthauvhrOOw8GD45dTf5S\nSMiCyy7zncauvjp2JSIixeHGG2HVKvj1r2NXkt8UErJg2209KNx9tw9kFBGRtvPZZ97Fe9FF8JWv\nxK4mvykkZMkFF/hW0pdfHrsSEZHC9utfwxZb+NRHaR2FhCzZYgsfxDhhArzxRuxqREQK0/vvw/33\n+z4NPXrErib/KSRk0be+BXvt5elWyzWLiGTeZZdB374wblzsSgqDQkIWtW8Pv/0tvPgiPPZY7GpE\nRArLM8/A44/Db34DnTvHrqYwKCRk2ahR8I1vwPnnw6JFsasRESkMy5bBd78LRx8NZ50Vu5rCoZCQ\nZWZw112wYYOvnyAiIq33s5/5wnV/+IOWX84khYQIeveGO+6ARx/VLpEiIq313HNw772+NkK/frGr\nKSwKCZGccQacfLIPrvn889jViIjkp+XL4Tvfga99Db7//djVFB6FhEjM4Pe/922kL7ggdjUiIvnp\n4oth8WLfbbed3tEyTpc0oj594He/g4oKH5ErIiIt98ILvpLtDTdA//6xqylMCgmRnXkmnHgijB0L\nS5bErkZEJD+sWOHdDEcc4a+f0jYUEiIz8yT8xRdw4YWxqxERyQ+XXuo7695/v7oZ2pIubQ7YYQe4\n/XYYPx6eeCJ2NSIiue2ll+DOO+H662HXXWNXU9gUEnLEmDFw/PHwgx+o20FEpDkrV8J558Hhh2vp\n5WxQSMgRiW6H1auhvDx2NSIiuemyy2DBAnUzZIsucQ75ylfg1lvhoYfgqadiVyMiklteftlnhF13\nHQwYELua4qCQkGO+9S047jjvdli6NHY1IiK5YeVKn80wbJiWtM8mhYQcY+bLiyb63erqYlckIhJX\nCD7+4LPP4I9/VDdDNulS56CddoI//Qn+/ne46qrY1YiIxJXohr3vPhg4MHY1xUUhIUedcILviX7l\nldoESkSK1zPPwE9/6ssvn3lm7GqKT4fYBUjzLrsMZszwcQoDB8K++8auSEQkez76CE4/HY49Fq69\nNnY1xUktCTnMzPvfBg6Ek07SbpEiUjyWLfMl67ffHh55BNq3j11RcUorJJjZ+WY228xWm9kUMztg\nE+eeYmbPmNl/zKzWzCab2Yj0Sy4u3br55k8rV8Jpp/mukSIihayuzheYmz/fV6Ht2TN2RcUr5ZBg\nZqcDNwFXAEOBGcAkM+vVzF0OB54BjgNKgBeBJ81sSFoVF6FddoHHHoOqKvjJT2JXIyLStn71K18r\npqIC9tgjdjXFLZ2WhHLgnhDCQyGE94GxwCrgvGQnhxDKQwj/G0KoDiF8HEK4HPgQOCHtqovQ4Yf7\nIiJ33OEjfEVECtGECXDNNb5g0siRsauRlEKCmXUESoHnE8dCCAF4DjikhY9hwFaAdihI0dixfhs3\nDl59NXY1IiKZNX06fPvbUFbmsxkkvlRbEnoB7YGFTY4vBPq08DF+BnQDJqT43ALcdhscfDB84xsw\nb17sakREMuM///EB2oMHe2upWeyKBLI8u8HMzgR+CZwaQtBY/TR06uTjE7p0gZNPhlWrYlckItI6\na9fC6NGwZo0P1O7aNXZFkpDqOgmfAxuA3k2O9wYWbOqOZnYGcC8wOoTwYkuerLy8nB49emx0rKys\njLKyshYXXIi2395XYzz0UF9DoaICOmjFCxHJQ3V18MMfwpQp8OKL0Ldv7IryS0VFBRUVFRsdq62t\nzdjjmw8pSOEOZlOAqSGEC+u/N2AucHsI4cZm7lMG3AecHkLY7P6GZlYCVFdXV1NSUpJSfcXk8cc9\nfZ9+OjxaSZ+LAAAQgElEQVT4oIKCiOSXujofY3Xvvf4advbZsSsqDDU1NZSWlgKUhhBqWvNY6byt\n3Aw8YGbVwOv4bIeuwAMAZnYdsGMI4Vv1359Z/7MLgDfMLNEKsTqEsKw1xRe7k0+GP/8ZzjjDv1dQ\nEJF80Tgg3H+/AkKuSvktJYQwoX5NhKvwbobpwLEhhEX1p/QBGjcYfQ8f7Hhn/S3hQZqZNiktN3q0\ngoKI5JemAeHcc2NXJM1J6+0khHAXcFczPzu3yfdHpvMc0nIKCiKSLxQQ8ov2bigQiaDw6KM+mHH9\n+tgViYhsTAEh/ygkFBAFBRHJVQoI+UkhocAoKIhIrlFAyF8KCQVIQUFEcoUCQn5TSChQjYPCOef4\nimYiItm0YYMvlKSAkL8UEgpYIig89hgcfbSvjS4ikg1Ll8KoUb4PgwJC/lJIKHCjR/tSp7NmwQEH\n+C5rIiJt6f334aCD4PXXYdIkBYR8ppBQBA47DN58E3r18v0eJmj/TRFpI5WVHhA6dvSQcPTRsSuS\n1lBIKBJ9+0JVlS/lfPrp8Itf+IAiEZFMCAF++1s44QQ44gh47TUYMCB2VdJaWpeviHTtCg8/DEOG\nwGWXwcyZMH48dO8euzIRyWerVsF3v+s70v7iF3DlldBOH0ELgv4Zi4wZXHIJPPkkvPwyHHIIfPRR\n7KpEJF/NmwfDh/v29RMmwNVXKyAUEv1TFqlRo2DqVF9D4cAD4dlnY1ckIvnm1Vdh//3h88/9z6ee\nGrsiyTSFhCI2aJAHhYMOgq9/Ha64QuspiMjmrV8PN9wARx7pryNvvgn77Re7KmkLCglFrmdPeOop\n+OUv4dprfZrktGmxqxKRXPXuuz5j6tJL4cILvRVyu+1iVyVtRSFBaN8efv1rn65k5kHhV79Sq4KI\nNFi/3mcvDB0KtbXevXDjjdCpU+zKpC0pJMh/DR3qQeGXv4TrrvO+xurq2FWJSGxvv+2DnC+/HC66\nyFsbDzkkdlWSDQoJspFOnXxswptvegvDQQf5lKYvvohdmYhk2/r13g1ZWgorV8LkyXD99bDFFrEr\nk2xRSJCkhgzxVoUrrvABSqWlHhxEpDi89RYcfLC3LP7kJ1BT4x8apLgoJEizOnb0F4g334TOnf0F\n43/+B5YsiV2ZiLSVZcu89bC0FNasgSlTvPuxS5fYlUkMCgmyWfvu6y8UV18N99wDu+3mrQurV8eu\nTEQy5Ysv4Pbb/f/3TTf57IXqah/ILMVLIUFapGNHX8r544/hrLN8ANPAgfDHP/qe8SKSn+rq4JFH\nYPBgKC+Hk06CDz+Eq67yFkQpbgoJkpLeveGOO+C993yu9He+4y0NTz7pG7yISP549lmfxXTWWf7/\n+K234L77YKedYlcmuUIhQdIyYAA8+qgPbuzdG048EQ4/3Ec/i0huq66GY46BESN8pkJVFTz+OOy5\nZ+zKJNcoJEirHHAAPP88PP00LF/urQsnnQSvvKKWBZFc88YbvlX8/vvDp5/C3/7m/1eHDYtdmeQq\nhQRpNTM49lifIvXQQzBrlu8Kd+CB3te5bl3sCkWK1/r18NhjHuAPPNCDwr33etfCySf7/1+R5igk\nSMa0awdnn+1ru0+cCFtv7X2d/fr5FKrFi2NXKFI8amt9lsKAAb47Y8eO3nLw4Yfwve9Bhw6xK5R8\noJAgGdeuHRx3HDzzjH9aGTkSrrwS+vaFsWPh/fdjVyhSuD7+GC64wAcfXnaZjxWqroaXXvKWg/bt\nY1co+UQhQdrU3nvDH/4A8+bBz38Of/+7T7U67jhvAl2zJnaFIvlv7Vr/v3XyyT41uaLC91j45BPv\nAiwpiV2h5CuFBMmK7bbzVdzmzPEXrSVLvAm0d2+fRvnCC1pvQSQVdXU+K2HsWOjTxwPCJ5/4eIO5\nc33xsx12iF2l5DuFBMmqzp193MLUqT7A8aKL4OWX4aijYJdd4Gc/g+nTNTNCpDnvvuutcrvu6l0J\n//iHB4W33/bdGb/7XW3AJJmjkCDR7L67j1X48EN47TU45RR44AHfsnqffXyw40cfxa5SJL45c3wQ\n4tChsNdecPfdPqPon/+E2bN9p8a99opdpRQihQSJzsw3j/rd7+Czz6Cy0ld/u/pq71/dfXdvcXjm\nGY1hkOKwdi28+CL89Ke+wFH//r4U+oABPkPh3//2fVSGD/eBwiJtRZNgJKd07OizIUaOhBUrfKzC\nxInw17/CbbdB165w9NH+8+OOg513jl2xSGZ89pl3HUyc6MslL1/uYw1GjvTAfMwx0L177Cql2Cgk\nSM7acktf7vnEE32Mwttv+wvoxIlw/vk+0HHvvb3ZdfhwXyymV6/YVYu0zJIlvoz5K6/ApEk+Fqdd\nO29Vu+QSDwdDhqilQOJSSJC8YObjFPbZx19A/+///NPWxIm+h8RNN/l5gwb5ErPDhnlo2G03rSgn\n8YXgMw9eeaXh9s47/rMddoCvfQ0uvtj3Uth227i1ijSmkCB5qWdPn0J56qn+/dy5G78A33+/vzD3\n7t0QGEpL/ZNZjx5xa5fCt3w5zJzpS5W/+qr/Ts6f7z/bc0//nbzkEv/ar5+CrOQuhQQpCDvvDGee\n6TeApUt9xsQrr/iL9M9/3jDosX9/2G8/vw0d6l932kkv1JK6EHwQ4fTpPv1w+nS/JWbldOrUsBXz\nsGFw6KFqKZD8opAgBWnrrRsGQIJvMjVrVsOL+LRpPhByyRL/+TbbeFjYc0/YYw+fUbHHHr6UtPqE\npa7OWwI++MB/jz74wNcrmD4dFi3yc3r29N+h449vCKGDB3tQEMlXCglSFDp29EGOe+8NY8b4sRB8\nu9zGweH5531qWWLnyi5dGqZhJsLDgAG+8NMOO2gd/EJSVwcLFvjYgY8/bggDs2b5Wh6rVvl5HTr4\nWJdBg2DcuIZAsMsuao2SwqOQIEXLzFsK+vaFE05oOL5+vb9RNP7UOGuWLyf96acN53Xo4N0UO+/s\nbxCNv+68s09f69lTbxy5IATfFXHBAt9HZO5c/zdu/HXevI23Nd9hBw+GBx0E55zTEBT79fPQKVIM\nFBLkvyoqKigrK4tdRnSJT4q77eZrMTS2YoWvfpd4Y0m8yfzrX77L3vz5/ok0oVMnHzzZp0/D18Sf\ne/eGGTMqOPXUMrbZxrs8unZVqGiJEGD1au8uStwWL4aFCz0ILFjQ8OfE1y++SNy7ArMydtyxIdAd\neODGIa9fP61JkEl6bclfaYUEMzsf+CnQB5gB/DiE8MYmzj8CuAnYC5gLXBNCeDCd55a2o//Im7fl\nlg3dFsmsW+eL4syb9+U3qwULvFsj8b1/aq3gmmsarnmnTj6eIhEattnGZ2NstZXfttyy+T936eJr\n9nfp0nDr3Dk3QkcIvorgmjX+5r5mTcOfV6702QDLl3sIS/y58fe1tQ1hYOlS/9rwpt+gQ4eNw9je\ne/viW43D2VVXVTBxYpnGCmSRXlvyV8ohwcxOx9/wvw+8DpQDk8xs9xDC50nO7wc8BdwFnAkcDdxn\nZp+FEJ5Nv3SR3NOxo38a3WWXTZ8Xgr/ZnXqq71HR9A2w8afjTz758ptmon+8JRqHho4d/Y10U7eW\nDNSsq/NumU3d1q3zN/JEKGipLbb4chDq3t1npZSWenBqGqQSt6233nz9XbtqMKFIS6XTklAO3BNC\neAjAzMYCo4DzgBuSnP9D4F8hhIvrv59lZsPqH0chQYqSmb+pdevmTd2p2rDBA0MiNKxY0fDpvPGt\n8af2NWs2fgNv7o29JTtwmrUscDRt2UjcEsc7d964RaRbN7+fiOSGlP47mllHoBS4NnEshBDM7Dng\nkGbudjDwXJNjk4BbUnluEWnQvr13Q2hhKBFpS6lm9l5Ae2Bhk+MLgT2auU+fZs7vbmadQwhJehbp\nAvDee++lWJ60Rm1tLTU1NbHLKCq65tmna559uubZ1ei9s0trHytXG/b6AYxJTGiXrCktLY1dQtHR\nNc8+XfPs0zWPoh8wuTUPkGpI+BzYAPRucrw3sKCZ+yxo5vxlzbQigHdHnAXMAVIY8iQiIlL0uuAB\nYVJrHyilkBBCWGdm1cBRwBMAZmb139/ezN1eA5rMNmdE/fHmnmcx8EgqtYmIiMh/taoFISGdVelv\nBr5nZueY2SDgbqAr8ACAmV1nZo3XQLgb2NXMrjezPcxsHDC6/nFEREQkR6U8JiGEMMHMegFX4d0G\n04FjQwj125zQB+jb6Pw5ZjYKn81wAfAp8J0QQtMZDyIiIpJDLLRkUrSIiIgUHW2CKyIiIkkpJIiI\niEhSORcSzOx8M5ttZqvNbIqZHRC7pkJhZsPN7Akzm29mdWZ2YpJzrjKzz8xslZk9a2YDYtRaKMzs\nMjN73cyWmdlCM/ubme2e5Dxd9wwxs7FmNsPMautvk83s603O0fVuI2Z2af3ry81NjuuaZ5CZXVF/\nnRvf3m1yTquveU6FhEabR10BDMV3mJxUP1BSWq8bPtB0HPClwShmdgnwI3zzrgOBlfj113Y46RsO\n/A44CN/crCPwjJltkThB1z3j5gGXACX4MvIvAH83s8Gg692W6j/UfR9/7W58XNe8bbyNTyDoU38b\nlvhBxq55CCFnbsAU4LZG3xs+G+Li2LUV2g2oA05scuwzoLzR992B1cBpsestlBu+tHkdMEzXPavX\nfTFwrq53m17jLYFZwNeAF4GbG/1M1zzz1/sKoGYTP8/INc+ZloRGm0c9nzgW/G+2qc2jJEPMrD+e\nRBtf/2XAVHT9M6kn3oqzBHTd25qZtTOzM/C1XCbrerepO4EnQwgvND6oa96mBtZ3H39sZuPNrC9k\n9prn0t4N6WweJZnTB3/zSnb9+2S/nMJTvzrprcArIYRE36Guexsws73xVV27AMuBU0IIs8zsEHS9\nM64+iO0H7J/kx/odbxtTgG/jrTc7AL8G/ln/u5+xa55LIUGk0N0F7AkcFruQIvA+MAToga/w+pCZ\nHR63pMJkZjvh4ffoEMK62PUUixBC430Z3jaz14FPgNPw3/+MyJnuBtLbPEoyZwE+BkTXvw2Y2R3A\nSOCIEMK/G/1I170NhBDWhxD+FUKYFkK4HB9IdyG63m2hFNgOqDGzdWa2DvgqcKGZrcU/veqat7EQ\nQi3wATCADP6e50xIqE+gic2jgI02j8rIRhXSvBDCbPyXp/H1746Pytf1b4X6gHAScGQIYW7jn+m6\nZ007oLOud5t4DtgH724YUn97ExgPDAkh/Atd8zZnZlviAeGzTP6e51p3w83AA/U7Tb4OlNNo8yhp\nHTPrhv8SWf2hXc1sCLAkhDAPbzL8hZl9hG/TfTU+u+TvEcotCGZ2F1AGnAisNLNEsq8NISS2Qdd1\nzyAzuxb4BzAX2Arfdv6r+O6zoOudUSGElUDT+fkrgcUhhPfqD+maZ5iZ3Qg8iXcxfAW4ElgH/Ln+\nlIxc85wKCWHzm0dJ6+yPT00K9beb6o8/CJwXQrjBzLoC9+Cj8KuA40IIa2MUWyDG4tf6pSbHzwUe\nAtB1z7jt8d/pHYBaYCYwIjHqXtc7KzZah0XXvE3sBDwCbAssAl4BDg4hLIbMXXNt8CQiIiJJ5cyY\nBBEREcktCgkiIiKSlEKCiIiIJKWQICIiIkkpJIiIiEhSCgkiIiKSlEKCiIiIJKWQICIiIkkpJIiI\niEhSCgkiIiKSlEKCiIiIJPX/AQi0yEXL0TKiAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110b93b10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pl.plot(hyb_pos(0, 0, np.linspace(0, beta)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "c_up = np.array([\n",
    "        [0, 1, 0, 0],\n",
    "        [0, 0, 0, 0],\n",
    "        [0, 0, 0, 1],\n",
    "        [0, 0, 0, 0],\n",
    "      ])\n",
    "c_down = np.array([\n",
    "        [0, 0, 1, 0],\n",
    "        [0, 0, 0,-1],\n",
    "        [0, 0, 0, 0],\n",
    "        [0, 0, 0, 0],\n",
    "      ])\n",
    "cdagger_up = c_up.T\n",
    "cdagger_down = c_down.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Operator:\n",
    "    def __init__(self, spin, typ, tau):\n",
    "        self.spin = spin\n",
    "        self.typ = typ\n",
    "        self.tau = tau\n",
    "        self.matrix = {\n",
    "                ('cdagger', 0): cdagger_up,\n",
    "                ('cdagger', 1): cdagger_down,\n",
    "                ('c', 0): c_up,\n",
    "                ('c', 1): c_down\n",
    "            }[(typ, spin)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "h = - mu * cdagger_up.dot(c_up) - mu * cdagger_down.dot(c_down) \\\n",
    "    + u * cdagger_up.dot(c_up).dot(cdagger_down).dot(c_down)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def evolve(h, tau):\n",
    "    return np.diag(np.exp(-h.diagonal() * tau))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_local_weight(opers):\n",
    "    global h\n",
    "    last_tau = 0\n",
    "    s = np.eye(4)\n",
    "    for oper in opers:\n",
    "        s = evolve(h, oper.tau - last_tau).dot(s)\n",
    "        s = oper.matrix.dot(s)\n",
    "        last_tau = oper.tau\n",
    "    s = evolve(h, beta - last_tau).dot(s)\n",
    "    return s.trace()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_hyb_antiperiodic(cspin, aspin, tau):\n",
    "    global beta\n",
    "    cspin, aspin, tau = np.broadcast_arrays(cspin, aspin, tau)\n",
    "    where_pos = tau >= 0\n",
    "    vals = np.empty_like(tau)\n",
    "    vals[where_pos] = hyb_pos(cspin[where_pos], aspin[where_pos],\n",
    "                              tau[where_pos])\n",
    "    vals[~where_pos] = -hyb_pos(cspin[~where_pos], aspin[~where_pos],\n",
    "                                tau[~where_pos] + beta)\n",
    "    return vals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_hyb_matrix(opers):\n",
    "    opers = sorted(opers, key=lambda op: op.typ)\n",
    "    spin = np.asarray([op.spin for op in opers])\n",
    "    tau = np.asarray([op.tau for op in opers])\n",
    "    nopershalf = spin.size//2\n",
    "    hyb_matrix = get_hyb_antiperiodic(\n",
    "        spin[nopershalf:, np.newaxis],\n",
    "        spin[np.newaxis, :nopershalf],\n",
    "        tau[nopershalf:, np.newaxis] - tau[np.newaxis, :nopershalf]\n",
    "        )\n",
    "    return hyb_matrix, tau, spin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_bath_weight(opers):\n",
    "    try:\n",
    "        hyb_matrix, tau, spin = get_hyb_matrix(opers)\n",
    "        return np.linalg.det(hyb_matrix)\n",
    "    except np.linalg.LinAlgError:\n",
    "        if np.size(opers): raise\n",
    "        return 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def add_pair(spin, tau1, tau2):\n",
    "    global opers, prop_opers, prop_weight\n",
    "    prop_opers = list(opers) + [\n",
    "                    Operator(spin, 'c', tau1),\n",
    "                    Operator(spin, 'cdagger', tau2)\n",
    "                    ]\n",
    "    prop_opers = sorted(prop_opers, key=lambda op: op.tau)\n",
    "    prop_weight = get_local_weight(prop_opers)\n",
    "    if prop_weight != 0:\n",
    "        prop_weight *= get_bath_weight(prop_opers)\n",
    "    # Ignore sign problem\n",
    "    prop_weight = np.abs(prop_weight)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def remove_pair(index1, index2):\n",
    "    # Convert bath to local index\n",
    "    global opers, prop_opers, prop_weight\n",
    "    lookup = np.argsort([op.typ for op in opers], kind='mergesort')\n",
    "    index1 = lookup[index1]\n",
    "    index2 = lookup[len(opers)//2 + index2]\n",
    "    prop_opers = np.delete(opers, [index1, index2])\n",
    "    prop_weight = get_local_weight(prop_opers)\n",
    "    if prop_weight != 0:\n",
    "        prop_weight *= get_bath_weight(prop_opers)\n",
    "    # Ignore sign problem\n",
    "    prop_weight = np.abs(prop_weight)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def qmc_step(rng):\n",
    "    global beta, opers, prop_weight, weight, prop_opers\n",
    "    move_typ = rng.randint(0, 2)\n",
    "    k = len(opers)//2\n",
    "    if move_typ == 0:\n",
    "        # INSERT\n",
    "        tau1, tau2 = rng.rand(2) * beta\n",
    "        spin = rng.randint(0, 2)\n",
    "        add_pair(spin, tau1, tau2)\n",
    "        ratio = (2 * beta**2)/(k + 1)**2\n",
    "    else:\n",
    "        # REMOVE\n",
    "        if not k:\n",
    "            return        \n",
    "        index1, index2 = rng.randint(0, k, 2)\n",
    "        remove_pair(index1, index2)\n",
    "        ratio = k**2/(2 * beta**2)\n",
    "    ratio *= np.abs(prop_weight/weight)\n",
    "    \n",
    "    if rng.rand() < ratio:\n",
    "        weight = prop_weight\n",
    "        opers = list(prop_opers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def meas_gtau():\n",
    "    global opers, gtau\n",
    "    hyb_matrix, tau, spin = get_hyb_matrix(opers)\n",
    "    nopershalf = tau.size//2\n",
    "    inv_matrix = np.linalg.inv(hyb_matrix)\n",
    "    indices = spin[nopershalf:, np.newaxis] * 2 * gtau.shape[-1] + \\\n",
    "              spin[np.newaxis, :nopershalf] * gtau.shape[-1] + \\\n",
    "              (tau[nopershalf:, np.newaxis] - \n",
    "               tau[np.newaxis, :nopershalf] + beta) * gtau.shape[-1]/(2*beta)\n",
    "    gtau += np.bincount(np.asarray(indices,int).ravel(), weights=inv_matrix.T.ravel(),\n",
    "                        minlength=gtau.size).reshape(gtau.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "opers = []\n",
    "weight = get_local_weight([])\n",
    "prop_opers = []\n",
    "prop_weight = 1.0\n",
    "\n",
    "gtau = np.zeros((2, 2, 100))\n",
    "histo = np.zeros(1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nmeas = 10000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rng = np.random.RandomState()\n",
    "for _ in xrange(nmeas):\n",
    "    for _ in xrange(10):\n",
    "        qmc_step(rng)\n",
    "    histo[len(opers)//2] += 1\n",
    "    meas_gtau()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ntau = gtau.shape[-1]//2\n",
    "gtau = gtau / (nmeas * beta * (beta/ntau))\n",
    "gtau = -gtau[0,0,:ntau] + gtau[0,0,ntau:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "histo /= nmeas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x110f72e50>]"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhYAAAFkCAYAAAB8RXKEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XuUHOV55/HvowuSJVkDWHgGDFjCFwnWIKMJCdhJwCGA\nbY4d78Y2HjZrHzDE2Ox6M44PCZsLa9gcTHKQAg4yJF4bE+xZs84mxiRBBgxmbW4nEogslsCAEFdd\nZgRSQEjo8u4f1W1a4+nRVE91V/fM93NOn1ZXv1X1qLjMb956630jpYQkSVIRppRdgCRJmjgMFpIk\nqTAGC0mSVBiDhSRJKozBQpIkFcZgIUmSCmOwkCRJhTFYSJKkwhgsJElSYQwWkiSpMA0Fi4i4MCLW\nRcSrEXFfRJwwStueiPhWRDwaEXsiYmmddl0RcU1EPB8ROyJibUS8v5H6JElSOXIHi4g4C7gSuAQ4\nHlgNrIiIeXV2mQFsAi4DHqpzzOnA7cCRwH8A3gmcDzyXtz5JklSeyLsIWUTcB9yfUvqvlc8BPANc\nnVL68/3seyfwYErpC8O2XwD8PrAopbQnV0GSJKlt5OqxqPQs9AJ3VLelLJncDpw0jjo+BNwLLI+I\nDRHxrxFxcUQ4BkSSpA4yLWf7ecBUYOOw7RuBheOo4yjgN4AbgQ8Abwe+WqnvspF2iIg3AWcATwE7\nxnFuSZImm5nAfGBFSmmoyAPnDRbNMoUsnPxupQfkwYg4HPgidYIFWaj4VovqkyRpIvqPwLeLPGDe\nYDEI7AG6h23vBjaMo44XgNfSvgM+1gA9ETEtpbR7hH2eArjxxhs5+uijx3Fq5dHf38+yZcvKLmNS\n8Zq3nte89bzmrbVmzRp+53d+Byo/S4uUK1iklHZFxErgVOBm+PngzVOBq8dRx0+AvmHbFgIv1AkV\nULn9cfTRR7NkyZJxnFp5dHV1eb1bzGveel7z1vOal6bwoQSNDI5cCpwfEZ+MiEXAtcAs4HqAiLg8\nIr5Zu0NELI6IdwNzgEMqn2u7Gb4KHBwRV0fEOyLiTOBi4K8aqE+SJJUk9xiLlNJNlTkrLiW7BfIQ\ncEZKaXOlSQ9wxLDdHgSqtzmWAGcD68kGbZJSejYizgCWkc2L8Vzlz6M+vipJktpLQ4M3U0rLgeV1\nvjtnhG377RlJKd0PvKeReiRJUntwngjl0tc3fCiMms1r3npe89bzmk8cuWfebBcRsQRYuXLlSgf8\nSJKUw6pVq+jt7QXoTSmtKvLY9lhIkqTCGCwkSVJhDBaSJKkwBgtJklQYg4UkSSqMwUKSJBXGYCFJ\nkgpjsJAkSYUxWEiSpMIYLCRJUmEMFpIkqTAGC0mSVBiDhSRJKozBQpIkFcZgIUmSCmOwkCRJhTFY\nSJKkwhgsJElSYQwWkiSpMAYLSZJUGIOFJEkqjMFCkiQVxmAhSZIKY7CQJEmFMVhIkqTCGCwkSVJh\nGgoWEXFhRKyLiFcj4r6IOGGUtj0R8a2IeDQi9kTE0v0c+xMRsTci/k8jtbW7556Dxx4ruwpJkpoj\nd7CIiLOAK4FLgOOB1cCKiJhXZ5cZwCbgMuCh/Rx7PvAXwN156+oU550HS5bA/feXXYkkScVrpMei\nH7gupXRDSmktcAGwHTh3pMYppfUppf6U0o3AtnoHjYgpwI3AnwLrGqir7W3fDnfeCRHwwQ/CI4+U\nXZEkScXKFSwiYjrQC9xR3ZZSSsDtwEnjrOUSYGNK6RvjPE7b+tGPYOdO+MEP4PDD4fTT4amnyq5K\nkqTi5O2xmAdMBTYO274R6Gm0iIj4VeAc4LxGj9EJbr0VjjwSTjwx+/PMmVm42LSp7MokSSrGtLIL\niIg5wA3A+SmlF/Pu39/fT1dX1z7b+vr66OvrK6jC4tx6K7z//dmtkEMPhdtug/e+N9t2110wd27Z\nFUqSJpqBgQEGBgb22bZ169amnS+yOxljbJzdCtkO/HZK6eaa7dcDXSmlf7+f/e8EHkwpfaFm22Jg\nFbAHiMrmak/KHmBhSukXxlxExBJg5cqVK1myZMmY/w5lefJJeNvb4O//Hj7ykde3P/xwNpjzmmvg\nM58prz5J0uSxatUqent7AXpTSquKPHauWyEppV3ASuDU6raIiMrnexqsYQ1wLPBuYHHldTPww8qf\nn2nwuG3l1lth2jT4jd/Yd/txx0FPT/YYqiRJna6RWyFLgesjYiXwANlTIrOA6wEi4nLgsJTSp6o7\nVHolApgDHFL5/FpKaU1K6TXgp7UniIiXyMaFrmmgvrZ0663ZbY+Rbnf09MCGDa2vSZKkouUOFiml\nmypzVlwKdJPNTXFGSmlzpUkPcMSw3R4EqvdclgBnA+uBoxoputPs3Ak//CH8yZ+M/H13N2wcPhxW\nkqQO1NDgzZTScmB5ne/OGWFb3lsuv3CMTvbjH8Mrr2SDNEfS0wP/7/+1tiZJkprBtUJa4NZbs/Bw\n3HEjf2+PhSRpojBYtEDtY6YjqY6xyPGAjiRJbclg0WTPPJPd5vjAB+q36e7OxmFsqzvhuSRJncFg\n0WQrVsCUKfCbv1m/TU9lzlKfDJEkdTqDRZPdeiv8yq/AwQfXb9Pdnb07zkKS1OkMFk20a1c2bXe9\np0Gq7LGQJE0UBosm+pd/ycZN7C9YdHXBjBn2WEiSOp/BoonWr8/ejz569HYR2e0QeywkSZ3OYNFE\nQ0MwfTrMmbP/tk7rLUmaCAwWTTQ0BG96U/35K2o5SZYkaSIwWDRRNViMhT0WkqSJwGDRRHmChT0W\nkqSJwGDRRHl7LDZudFpvSVJnM1g0Ud4ei1274MUXm1uTJEnNZLBoorw9FuA4C0lSZzNYNNHQ0OhT\neddyWm9J0kRgsGiS3bth61Z7LCRJk4vBokm2bMnexxos5syBWbPssZAkdTaDRZPkDRZO6y1JmggM\nFk0yNJS9jzVYwOuPnEqS1KkMFk3SSLCwx0KS1OkMFk1SDRZjfSoE7LGQJHU+g0WTDA3B3LnZ6qZj\nZY+FJKnTGSyaJM/kWFU9PbBpE+zd25yaJElqNoNFkzQSLLq7Yc+e12+jSJLUaQwWTZJn1s0qJ8mS\nJHU6g0WTNNpjAQ7glCR1LoNFk4wnWNhjIUnqVA0Fi4i4MCLWRcSrEXFfRJwwStueiPhWRDwaEXsi\nYukIbc6LiLsjYkvlddtox+wEjQSL2bPhjW+0x0KS1LlyB4uIOAu4ErgEOB5YDayIiHl1dpkBbAIu\nAx6q0+Zk4NvAKcCJwDPADyLi0Lz1tYOUsim98wYL8JFTSVJna6THoh+4LqV0Q0ppLXABsB04d6TG\nKaX1KaX+lNKNwLY6bf5TSunalNLDKaXHgPMqtZ3aQH2l274ddu5sLFg4SZYkqZPlChYRMR3oBe6o\nbkspJeB24KQC65oNTAe2FHjMlmlkOu8qeywkSZ0sb4/FPGAqMPx36o1ATyEVZa4AniMLLB1nPMHC\nHgtJUiebVnYBw0XEHwIfB05OKb22v/b9/f10dXXts62vr4++vr4mVbh/9lhIktrFwMAAAwMD+2zb\nunVr086XN1gMAnuA7mHbu4Fx/ziMiC8CFwGnppQeGcs+y5YtY8mSJeM9daHG22MxOJjNwDl1arF1\nSZImn5F+2V61ahW9vb1NOV+uWyEppV3ASmoGVUZEVD7fM55CIuIi4I+AM1JKD47nWGUbGoJp07JH\nR/Pq7s7WCtm8ufi6JElqtkZuhSwFro+IlcADZE+JzAKuB4iIy4HDUkqfqu4QEYuBAOYAh1Q+v5ZS\nWlP5/g+ALwF9wNMRUe0ReTml9Eojf7EyVafzjsi/b3Va740bX/+zJEmdInewSCndVJmz4lKyWyAP\nkfUyVH/H7gGOGLbbg0Cq/HkJcDawHjiqsu0CsqdAvjtsvy9VztNRGpkcq6p29s3Fi4urSZKkVmho\n8GZKaTmwvM5354ywbdRbLimlBY3U0a6KCBY+GSJJ6kSuFdIE4wkWM2fCgQf6ZIgkqTMZLJpgPMEC\nfORUktS5DBZN0Og6IVVOkiVJ6lQGiyawx0KSNFkZLAq2Zw+89JI9FpKkyclgUbAXX8yWTbfHQpI0\nGRksCjae6byrenqy4+zaVUxNkiS1isGiYNVgcfDBjR+jOpeF03pLkjqNwaJgRfRYvPnN2bvjLCRJ\nncZgUbAigkW1x2LTpvHXI0lSKxksCjY0BHPmwAEHNH4MeywkSZ3KYFGw8c5hAdm03nPn2mMhSeo8\nBouCFREsIOu1sMdCktRpDBYFKypYdHfbYyFJ6jwGi4KNd52QKnssJEmdyGBRMHssJEmTmcGiYI6x\nkCRNZgaLgg0NjW/Wzapqj0VK4z+WJEmtYrAo0PbtsGNHcT0Wu3dni5pJktQpDBYFKmLWzSpn35Qk\ndSKDRYGKDBbOvilJ6kQGiwLZYyFJmuwMFgUqMlh0dWXrjdhjIUnqJAaLAg0NwdSpWSgYr4jsdog9\nFpKkTmKwKFD1UdOIYo7nXBaSpE5jsChQUdN5Vzn7piSp0xgsClTUrJtV9lhIkjpNQ8EiIi6MiHUR\n8WpE3BcRJ4zSticivhURj0bEnohYWqfdxyJiTeWYqyPiA43UVqaig4U9FpKkTpM7WETEWcCVwCXA\n8cBqYEVEzKuzywxgE3AZ8FCdY74H+DbwN8C7ge8B/xARx+Str0xFTedd1d1tj4UkqbM00mPRD1yX\nUrohpbQWuADYDpw7UuOU0vqUUn9K6UZgW51jfh7455TS0pTSoymlPwVWAf+5gfpK04xbIS+/nE0V\nLklSJ8gVLCJiOtAL3FHdllJKwO3ASeOo46TKMWqtGOcxW64Zt0LA2yGSpM6Rt8diHjAVGN5BvxHo\nGUcdPU04Zkvt3ZstGFZ0jwUYLCRJncOnQgry0kvZEudFj7EAx1lIkjrHtJztB4E9QPew7d3AhnHU\nsaHRY/b399M1bKrLvr4++vr6xlFOflu2ZO9F9ljMm5dNtmWPhSSpUQMDAwwMDOyzbevWrU07X65g\nkVLaFRErgVOBmwEiIiqfrx5HHfeOcIzTKttHtWzZMpYsWTKOUxejuk5IkT0W06ZlQcUeC0lSo0b6\nZXvVqlX09vY25Xx5eywAlgLXVwLGA2RPicwCrgeIiMuBw1JKn6ruEBGLgQDmAIdUPr+WUlpTaXIV\ncFdEfAH4R6CPbJDo+Y38pcrQjB4LcL0QSVJnyR0sUko3VeasuJTsdsVDwBkppc2VJj3AEcN2exBI\nlT8vAc4G1gNHVY55b0ScDfxZ5fUz4LdSSj/NW19ZmtFjAc5lIUnqLI30WJBSWg4sr/PdOSNs2+8g\n0ZTS3wF/10g97WDLFpg5E2bNKva49lhIkjqJT4UUZMuW4nsrwB4LSVJnMVgUpOjJsarssZAkdRKD\nRUGa2WMxOAi7dxd/bEmSimawKEgzeyxSen1wqCRJ7cxgUZBm9liA4ywkSZ3BYFGQZvZYgOMsJEmd\nwWBRkGb1WFSDhT0WkqROYLAowO7d2SJkzQgWs2dnL3ssJEmdwGBRgJdeyt6bcSsEnMtCktQ5DBYF\nqK4T0oweC3AuC0lS5zBYFKD6KKg9FpKkyc5gUQB7LCRJyhgsCtCslU2r7LGQJHUKg0UBtmzJVjWd\nObM5x6/2WKS0/7aSJJXJYFGAZk2OVdXdDTt3wrZtzTuHJElFMFgUoFmTY1U5+6YkqVMYLArQih4L\ncJyFJKn9GSwKYI+FJEkZg0UBtmxpbo/FQQfBtGn2WEiS2p/BogBDQ83tsZgyBQ45xB4LSVL7M1gU\noNm3QsC5LCRJncFgMU67dmWPgTbzVghk4ywMFpKkdmewGKcXX8zem91jccQR8MwzzT2HJEnjZbAY\np2YvQFa1YAGsW9fcc0iSNF4Gi3Fq9gJkVfPnw+AgvPxyc88jSdJ4GCzGqZU9FmCvhSSpvRksxqna\nY3HQQc09TzVYPPVUc88jSdJ4GCzGacsWeOMb4YADmnue7m6YMcMeC0lSe2soWETEhRGxLiJejYj7\nIuKE/bQ/JSJWRsSOiHgsIj41Qpvfi4i1EbE9Ip6OiKURMaOR+lqp2ZNjVU2Zko2zMFhIktpZ7mAR\nEWcBVwKXAMcDq4EVETGvTvv5wC3AHcBi4CrgaxFxWk2bs4HLK8dcBJwLfBz4s7z1tVqzp/OutWCB\nt0IkSe2tkR6LfuC6lNINKaW1wAXAdrIwMJLPAk+mlC5KKT2aUroG+G7lOFUnAT9OKX0npfR0Sul2\n4H8Bv9xAfS3Vqh4LsMdCktT+cgWLiJgO9JL1PgCQUkrA7WThYCQnVr6vtWJY+3uA3uotlYg4Cvgg\n8I956itDK6bzrqrOZZFSa84nSVJe03K2nwdMBYZPLr0RWFhnn5467edGxIyU0s6U0kDlVsqPIyIq\n57g2pXRFzvpabmgI3vnO1pxrwYJs+vCXXmr+UyiSJDUib7Boiog4BfhvZLdVHgDeDlwdES+klP7H\naPv29/fT1dW1z7a+vj76+vqaVO2+WtljMX9+9r5uncFCkjQ2AwMDDAwM7LNt69atTTtf3mAxCOwB\nuodt7wY21NlnQ53221JKOyufLwX+NqX0jcrnRyJiDnAdMGqwWLZsGUuWLBlj+cUbGmrt4E3IgkWJ\nf2VJUgcZ6ZftVatW0dvb25Tz5RpjkVLaBawETq1uq9y6OJVsnMRI7q1tX3F6ZXvVLGD3sDZ7a47f\nlnbuhFdeaV2PxZveBHPmOIBTktS+GrkVshS4PiJWkt226CcLBtcDRMTlwGEppepcFdcCF0bEFcDX\nyULGR8kGZ1Z9H+iPiNXA/cA7yHoxbq4MDm1L1ZVNW9VjEZHdDvGRU0lSu8odLFJKN1UGWl5Kdkvj\nIeCMlNLmSpMe4Iia9k9FxJnAMuDzwLPApyuPlFZdRtZDcRnwFmAzcDPwx7n/Ri1UXSekVT0W4Cqn\nkqT21tDgzZTScmB5ne/OGWHb3WSPqdY7XjVUXNZIPWWprhPSqh4LyILFbbe17nySJOXhWiHjUEaP\nRfVWSPveIJIkTWYGi3Fo1cqmtRYsgFdfhU2bWndOSZLGymAxDkND0NUF01o4G0jtI6eSJLUbg8U4\ntHIBsqrqJFk+GSJJakcGi3Fo5QJkVV1d2a0XeywkSe3IYDEOrZzOu5aPnEqS2pXBYhzKuBUCLp8u\nSWpfBotxKONWCGQ9Fo6xkCS1I4PFOJTVY7FgAaxfD3v2tP7ckiSNxmAxDmX1WMyfD7t2wfPPt/7c\nkiSNxmDRoFdfzV5l9ViAt0MkSe3HYNGg6qybZfVYgAM4JUntx2DRoDIWIKuaNQve/GaDhSSp/Rgs\nGlTGAmS1fDJEktSODBYNKrPHApwkS5LUngwWDRoagohsiu0yOEmWJKkdGSwatGVLtmbH1KnlnH/B\nAnj22eyxU0mS2oXBokFlrRNStWAB7N0LzzxTXg2SJA1nsGhQWZNjVfnIqSSpHRksGjQ4CIccUt75\njzwye7fHQpLUTgwWDRochHnzyjv/jBkwdy5s3lxeDZIkDWewaFDZwQKyHhODhSSpnRgsGjQ4WN4c\nFlUGC0lSuzFYNGDPHnjxRXssJEkazmDRgJdeyh71NFhIkrQvg0UDBgezd4OFJEn7Mlg0wGAhSdLI\nDBYNaKdg8fLLsGNHuXVIklTVULCIiAsjYl1EvBoR90XECftpf0pErIyIHRHxWER8aoQ2XRFxTUQ8\nX2m3NiLe30h9zTY4mC1AdtBB5dZRnaDLXgtJUrvIHSwi4izgSuAS4HhgNbAiIkb8/T0i5gO3AHcA\ni4GrgK9FxGk1baYDtwNHAv8BeCdwPvBc3vpaYXAwCxXTppVbh8FCktRuGvnR2A9cl1K6ASAiLgDO\nBM4F/nyE9p8FnkwpXVT5/GhE/GrlOLdVtn0aOBA4MaW0p7Lt6QZqa4l2mMMCDBaSpPaTq8ei0rPQ\nS9b7AEBKKZH1NpxUZ7cTK9/XWjGs/YeAe4HlEbEhIv41Ii6OiLYcA9IOs27C6zVUx3xIklS2vD+4\n5wFTgY3Dtm8Eeurs01On/dyImFH5fBTwsUo9HwAuBX4f+KOc9bXE0FB7BIvZs+ENb7DHQpLUPkoe\nJfBzU8jCxu9WekAejIjDgS8Cl422Y39/P11dXfts6+vro6+vr1m1MjgIixY17fC5+MipJGk0AwMD\nDAwM7LNt69atTTtf3mAxCOwBuodt7wY21NlnQ53221JKOyufXwBeq4SKqjVAT0RMSyntrlfQsmXL\nWLJkyVjrL0S73AoBg4UkaXQj/bK9atUqent7m3K+XLdCUkq7gJXAqdVtERGVz/fU2e3e2vYVp1e2\nV/0EePuwNguBF0YLFWUxWEiSNLJGBkcuBc6PiE9GxCLgWmAWcD1ARFweEd+saX8tcFREXBERCyPi\nc8BHK8ep+ipwcERcHRHviIgzgYuBv2qgvqbavbs9FiCrMlhIktpJ7jEWKaWbKnNWXEp2S+Mh4IyU\nUvXHWw9wRE37pypBYRnweeBZ4NMppdtr2jwbEWdU2qwmm79iGSM/vlqqLVuy93YKFvffX3YVkiRl\nGhq8mVJaDiyv8905I2y7m+wx1dGOeT/wnkbqaaV2mc67yh4LSVI7act5ItpZNVi0wwRZkAWLF1+E\nXbvKrkSSJINFbu3YYwHZ3BqSJJXNYJHT0BBMmQIHHlh2JRmn9ZYktRODRU6Dg3DwwTB1atmVZAwW\nkqR2YrDIqZ3msACDhSSpvRgscmq3YDF3LkyfbrCQJLUHg0VO7RYsInzkVJLUPgwWObVbsACDhSSp\nfRgscjJYSJJUn8Eip8HB9pkcq8pgIUlqFwaLHHbtgm3b7LGQJKkeg0UO1dktDRaSJI3MYJFDu03n\nXXXIIVno2bu37EokSZOdwSKHdg4We/e+vqS7JEllMVjk0M7BArwdIkkqn8Eih8HBbI2Qrq6yK9mX\nwUKS1C4MFjlU57CIKLuSfRksJEntwmCRQztOjgVw0EFZT4rBQpJUNoNFDu04ORbAlClZXQYLSVLZ\nDBY5DA21Z48FOJeFJKk9GCxyaNdbIWCwkCS1B4NFDgYLSZJGZ7DIwWAhSdLoDBZjtGMHvPxy+waL\nefMMFpKk8hksxqhdFyCrOuSQrEclpbIrkSRNZgaLMWrX6byrDjnk9WXdJUkqi8FijDohWIC3QyRJ\n5TJYjFH1Vkg7TpAFBgtJUntoKFhExIURsS4iXo2I+yLihP20PyUiVkbEjoh4LCI+NUrbT0TE3oj4\nP43U1iyDgzB9OrzxjWVXMjKDhSSpHeQOFhFxFnAlcAlwPLAaWBERI94kiIj5wC3AHcBi4CrgaxFx\nWp22fwHcnbeuZmvXBciqqj0pBgtJUpka6bHoB65LKd2QUloLXABsB86t0/6zwJMppYtSSo+mlK4B\nvls5zs9FxBTgRuBPgXUN1NVU7TyHBWS9KQcdZLCQJJUrV7CIiOlAL1nvAwAppQTcDpxUZ7cTK9/X\nWjFC+0uAjSmlb+SpqVXaPViAk2RJkso3LWf7ecBUYOOw7RuBhXX26anTfm5EzEgp7YyIXwXOIbtV\n0pYMFpIk7V/eYFG4iJgD3ACcn1J6Me/+/f39dHV17bOtr6+Pvr6+girMDA7C299e6CELZ7CQJA03\nMDDAwMDAPtu2bt3atPPlDRaDwB6ge9j2bmBDnX021Gm/rdJbsQh4K/D9iJ8PjZwCEBGvAQtTSnXH\nXCxbtowlS5bk+1s0oFN6LFatKrsKSVI7GemX7VWrVtHb29uU8+UaY5FS2gWsBE6tbquEgVOBe+rs\ndm9t+4rTK9sB1gLHAu8muxWyGLgZ+GHlz8/kqbFZBgfbdw6LKnssJElla+RWyFLg+ohYCTxA9nTH\nLOB6gIi4HDgspVSdq+Ja4MKIuAL4OlnI+CjwQYCU0k7gp7UniIiXsq/SmgbqK9z27fDqq53RY7F5\nc7ZeSLs+FitJmthyB4uU0k2VOSsuJbul8RBwRkqp+rtyD3BETfunIuJMYBnweeBZ4NMppeFPirSt\ndl+ArOqtb80C0IYNcOihZVcjSZqMGhq8mVJaDiyv8905I2y7m+wx1bEe/xeOUaZ2Xyek6rjjsveH\nHzZYSJLK4VohY/Dcc9l7T0+5dezPggUwZw6sXl12JZKkycpgMQZPPAEzZ7Z/L8CUKXDssVmPhSRJ\nZTBYjMHjj8Pb3pb94G53xx1nsJAklacDflSW74knsmDRCY47DtasgddeK7sSSdJkZLAYg04LFrt3\nw9q1ZVciSZqMDBb7sWcPrFvX/tN5Vx17bPbuAE5JUhkMFvvxzDOwa1fn9Fh0dcH8+Y6zkCSVw2Cx\nH088kb13SrAAB3BKkspjsNiPxx+HqVOzWS07hcFCklQWg8V+PPFEFiqmTy+7krE77rhsWu9Nm8qu\nRJI02Rgs9qOTngipWrw4e//Xfy23DknS5GOw2I/q5Fid5G1vgze8wSdDJEmtZ7AYRUpZj0WnPGpa\nNXUqvOtdjrOQJLWewWIUmzbBK690Xo8FOIBTklQOg8UoHn88e+/UYPHII9ksnJIktYrBYhTVOSyO\nOqrcOhqxeHG2Xshjj5VdiSRpMjFYjOKJJ7Kl0mfPLruS/JzaW5JUBoPFKDrxiZCqgw+Gww93nIUk\nqbUMFqPoxCdCajmAU5LUagaLUXTi5Fi1DBaSpFYzWNSxdSsMDnZ2sFi8GJ59FrZsKbsSSdJkYbCo\noxNXNR3uuOOyd6f2liS1isGijmqw6OQxFu98JxxwgE+GSJJax2BRx+OPw4EHZk9XdKpp0+Df/TvH\nWUiSWsdgUUenD9ysOuYYWLu27CokSZOFwaKOTn/UtGrhQoOFJKl1DBZ1dPLkWLUWLoShoewlSVKz\nNRQsIuLCiFgXEa9GxH0RccJ+2p8SESsjYkdEPBYRnxr2/XkRcXdEbKm8btvfMZtpxw547rmJEywA\nHn203DokSZND7mAREWcBVwKXAMcDq4EVETGvTvv5wC3AHcBi4CrgaxFxWk2zk4FvA6cAJwLPAD+I\niEPz1leEdesgpYlxK+Qd78jeDRaSpFZopMeiH7gupXRDSmktcAGwHTi3TvvPAk+mlC5KKT2aUroG\n+G7lOACklP5TSunalNLDKaXHgPMqtZ3aQH3j1snLpQ83axYceaTBQpLUGrmCRURMB3rJeh8ASCkl\n4HbgpDqL+pX6AAAMz0lEQVS7nVj5vtaKUdoDzAamA6XMGfnEEzBzZray6USwaJHBQpLUGnl7LOYB\nU4GNw7ZvBHrq7NNTp/3ciJhRZ58rgOf4xUDSEtVHTadMkKGtPhkiSWqVtvvRGRF/CHwc+EhK6bUy\napgoT4RULVyYhaXdu8uuRJI00U3L2X4Q2AN0D9veDWyos8+GOu23pZR21m6MiC8CFwGnppQeGUtB\n/f39dHV17bOtr6+Pvr6+sew+oscfhw99qOHd287ChbBrVzYotTqYU5I0OQwMDDAwMLDPtq1btzbt\nfLmCRUppV0SsJBtUeTNARETl89V1drsX+MCwbadXtv9cRFwEXAycnlJ6cKw1LVu2jCVLloy1+X5t\n3pwFi+OPL+yQpat95NRgIUmTy0i/bK9atYre3t6mnK+RWyFLgfMj4pMRsQi4FpgFXA8QEZdHxDdr\n2l8LHBURV0TEwoj4HPDRynGo7PMHwKVkT5Y8HRHdldfshv5W43DXXdn7+97X6jM3z1vekj0d4gBO\nSVKz5b0VQkrppsqcFZeS3dJ4CDgjpbS50qQHOKKm/VMRcSawDPg88Czw6ZRS7cDMC8ieAvnusNN9\nqXKelrnrrmz+isMPb+VZm2vKlGylU4OFJKnZcgcLgJTScmB5ne/OGWHb3WSPqdY73oJG6miGO++c\nWL0VVQsXGiwkSc3Xdk+FlGnDBlizZmIGi0WLfORUktR8BosaP/pR9n7KKaWW0RQLF8KmTfDSS2VX\nIkmayAwWNe68M/sBPFFm3KzlYmSSpFYwWNSYqOMrIBu8CQYLSVJzGSwqnn8eHntsYt4GAZgzJ3vs\n1GAhSWomg0VFdf6KiRoswCdDJEnNZ7CouPNOOOYY6B4++fgEYrCQJDWbwaLirrsm7viKqkWL4Gc/\ngz17yq5EkjRRGSyAZ5/N1geZyLdBIOux2LkT1q8vuxJJ0kRlsCC7DQKTI1iAt0MkSc1jsCC7DXLs\nsTBvXtmVNNeRR8LMmQYLSVLzGCyY2PNX1JoyJVs23WAhSWqWSR8s1q+Hdesm/m2QKp8MkSQ106QP\nFnfdBRFw8sllV9IaBgtJUjNN+mDxT/8Exx8PBx9cdiWtsXBhNsvotm1lVyJJmogmdbB45RW45Rb4\n2MfKrqR1Fi3K3h97rNw6JEkT06QOFv/4j7B9O3z842VX0jrHHANz58LXv152JZKkiWhSB4vvfAd+\n6ZfgqKPKrqR1Zs+GL30JrrsOVq8uuxpJ0kQzaYPFv/1bNr7irLPKrqT1LrwwG2vx+c9DSmVXI0ma\nSCZtsPj+92HHjsk1vqJq+nS46iq4+2646aayq5EkTSSTNljcdBOceCK89a1lV1KO006Dj3wEvvjF\nbBCrJElFmJTBYutW+Od/nlyDNkdy5ZWweTN8+ctlVyJJmigmZbC4+WZ47bXJeRuk1lFHZT0Wf/EX\n8OSTZVcjSZoIJmWw+M534L3vhcMPL7uS8l18cbb42n/5L7BrV9nVSJI63aQLFi++CD/4weR8GmQk\ns2fD8uXZNTn5ZHj66bIrkiR1skkXLP7hH2D3bvjt3y67kvbx4Q/D//2/8Nxz2fTmt9xSdkWSpE41\n6YLFTTfBr/86HHZY2ZW0lxNPhAcfzG4RfehDcNFF3hqRJOU3qYLF+vVw++0+DVLPwQfD976XPS2y\nbBm8733wwgv7thkYGCinuEnMa956XvPW85pPHA0Fi4i4MCLWRcSrEXFfRJywn/anRMTKiNgREY9F\nxKdGaPOxiFhTOebqiPhAI7XV8+STcMop8Ja3wCc+UeSRJ5YI+MIX4Ec/yq5Zby/cc8/r3/sff+t5\nzVvPa956XvOJI3ewiIizgCuBS4DjgdXAioiYV6f9fOAW4A5gMXAV8LWIOK2mzXuAbwN/A7wb+B7w\nDxFxTN76RrJ2Lfzar8EBB2RjCSbLEunj8Z73wMqV2SOpp5wCX/2q039LkvavkR6LfuC6lNINKaW1\nwAXAduDcOu0/CzyZUroopfRoSuka4LuV41R9HvjnlNLSSps/BVYB/7mB+vbx8MPZmIqDD85+Cz/i\niPEecfI49FD44Q/hM5+Bz30Ozj3XcReSpNHlChYRMR3oJet9ACCllIDbgZPq7HZi5ftaK4a1P2kM\nbXLZtg1uuy37bfuII+Cuu6Cnp9GjTV4HHABf+Qp885vZ/B+33Qa/93vZeBVJkoablrP9PGAqsHHY\n9o3Awjr79NRpPzciZqSUdo7SZrQoMBPg6qvXMGNGNj/F5s3ZYMPnn89WLwU49lhYujT7QegPw8a9\n613ZwM7zztvKN76xiq98BU4/Hc48Mwsfap7167fy13+9quwyJhWveet5zZtr3jyYP//1z2vWrKn+\ncWbhJ0spjfkFHArsBX5l2PYrgHvr7PMo8AfDtn0A2APMqHzeCZw1rM1ngRdGqeVsIPny5cuXL1++\nGn6dnScHjOWVt8dikCwQdA/b3g1sqLPPhjrtt1V6K0ZrU++YkN0q+Y/AU8COUauWJEm1ZgLzyX6W\nFipXsEgp7YqIlcCpwM0AERGVz1fX2e1esh6KWqdXtte2GX6M04a1GV7LENmTJJIkKb979t8kv0ae\nClkKnB8Rn4yIRcC1wCzgeoCIuDwivlnT/lrgqIi4IiIWRsTngI9WjlN1FfD+iPhCpc1/Jxsk+lcN\n1CdJkkqS91YIKaWbKnNWXEp2u+Ih4IyU0uZKkx7giJr2T0XEmcAyssdKnwU+nVK6vabNvRFxNvBn\nldfPgN9KKf20sb+WJEkqQyRnPZIkSQWZVGuFSJKk5jJYSJKkwnRksMi7CJrGLiIujogHImJbRGyM\niL+PiHeO0O7SiHg+IrZHxG0R8fYy6p1oIuIPI2JvRCwdtt3rXbCIOCwi/jYiBivXdXVELBnWxute\nkIiYEhGXRcSTlev5eET88QjtvOYNiohfi4ibI+K5yv9HPjxCm1Gvb0TMiIhrKv9d/FtEfDci3pyn\njo4LFnkXQVNuvwZ8BfgV4DeB6cAPIuIN1QYR8Qdk67j8LvDLwCtk/wycg3McKgH5d8n+na7d7vUu\nWEQcCPyEbHK+M4Cjgd8HXqxp43Uv1h8CnwE+BywCLgIuioifrwnlNR+32WQPVHyObPKrfYzx+v4l\ncCbw28CvA4cBf5eriqJn3Gr2C7gPuKrmc5A9aXJR2bVNxBfZNO57gV+t2fY80F/zeS7wKvDxsuvt\n1Bcwh2yW2t8A7gSWer2ber2/DPxoP2287sVe8+8DfzNs23eBG7zmTbnee4EPD9s26vWtfN4J/Pua\nNgsrx/rlsZ67o3osGlwETeNzIFny3QIQEQvIHimu/WewDbgf/xmMxzXA91NKP6zd6PVumg8B/xIR\nN1Vu+a2KiPOqX3rdm+Ie4NSIeAdARCwG3gv8U+Wz17yJxnh9f4lsGoraNo8CT5Pjn0HueSxK1sgi\naGpQZVbVvwR+nF6fU6SHLGjkXTROdUTEJ4B3k/1HPZzXuzmOIluP6EqyuXN+Gbg6InamlP4Wr3sz\nfJnsN+K1EbGH7Fb8H6WU/lfle695c43l+nYDr1UCR702+9VpwUKttRw4huy3CjVBRBxOFt5+M6W0\nq+x6JpEpwAMppT+pfF4dEe8CLgD+tryyJrSzyBaP/ATwU7IwfVVEPF8Jc5ogOupWCI0tgqYGRMRf\nAR8ETkkpvVDz1QaycS3+MyhGL3AIsCoidkXELuBk4L9GxGtkvyl4vYv3ArBm2LY1wJGVP/vvefH+\nHPhySul/p5QeSSl9i2xG5osr33vNm2ss13cDcEBEzB2lzX51VLCo/EZXXQQN2GcRtKYspjIZVULF\nbwHvSyk9XftdSmkd2b9gtf8M5pI9ReI/g/xuB44l++1tceX1L8CNwOKU0pN4vZvhJ/zi7dOFwHrw\n3/MmmUX2i2GtvVR+DnnNm2uM13clsHtYm4VkgbvuoqDDdeKtkKXA9ZVVVh8A+qlZBE3jExHLgT7g\nw8ArEVFNt1tTStXl6f8S+OOIeJxs2frLyJ7M+V6Ly+14KaVXyLqFfy4iXgGGUkrV36i93sVbBvwk\nIi4GbiL7n+t5wPk1bbzuxfo+2fV8FngEWEL2/++v1bTxmo9DRMwG3k7WMwHZAqCLgS0ppWfYz/VN\nKW2LiP8JLI2IF4F/I1t1/CcppQfGXEjZj8Q0+BjN5yoX5VWyFPVLZdc0UV5kv0HsGeH1yWHt/jvZ\no0vbgRXA28uufaK8gB9S87ip17tp1/mDwMOVa/oIcO4IbbzuxV3v2WS/GK4jmz/hZ8CXgGle88Ku\n8cl1/h/+9bFeX2AG2VxGg5Vg8b+BN+epw0XIJElSYTpqjIUkSWpvBgtJklQYg4UkSSqMwUKSJBXG\nYCFJkgpjsJAkSYUxWEiSpMIYLCRJUmEMFpIkqTAGC0mSVBiDhSRJKsz/B05cp7qYp4U6AAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110b938d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pl.plot(histo[:100])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1111eded0>]"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhIAAAFkCAYAAAB1rtL+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XucnPP5//HXlchBhFAhkUqJarBIsllpE4cg6tiWUodu\naJWvoChdp9YXVVryaysJRSpFhWLbqNahqnEuIaESu4icRAgiZ9mQgxz2+v3xmflmsvYwc8/s3DOz\n7+fjMY/J3HPfn/vaO7M71/05mrsjIiIiEkW7uAMQERGR4qVEQkRERCJTIiEiIiKRKZEQERGRyJRI\niIiISGRKJERERCQyJRIiIiISmRIJERERiUyJhIiIiESmREJEREQii5RImNn5ZjbPzNaY2RQzG9TM\nvgebWX2Dx0Yz27HBfieZ2YxEmbVmdnSU2ERERCR/Mk4kzOwUYBRwDVAO1AITzax7M4c58DWgZ+Kx\nk7svTilzf+AB4A5gAPAI8LCZlWUan4iIiOSPZbpol5lNAV5x94sSrw34APi9u/+2kf0PBp4FtnP3\nlU2U+Regi7sfm7JtMvC6u5+XUYAiIiKSNxnVSJhZB6ACeCa5zUMm8jQwpLlDgRozW2BmTyZqIFIN\nSZSRamILZYqIiEjMtshw/+5Ae2BRg+2LgD2aOOZj4BzgNaATMAJ43sy+7u41iX16NlFmz6YCMbPt\ngSOB94C16f8IIiIibV5nYFdgorsvy6agTBOJjLn7bGB2yqYpZvZVoAo4PYuijwTuzyY2ERGRNu5U\nQh/FyDJNJJYCG4EeDbb3ABZmUM6rwAEprxdGKPM9gPvuu4+99torg1NLNqqqqhgzZkzcYbQpuub5\np2uef7rm+TVjxgxOO+00SHyXZiOjRMLd15vZVOAw4FH4v86WhwG/z6CoAYQmj6TJjZRxeGJ7U9YC\n7LXXXgwcODCDU0s2unXrpuudZ7rm+adrnn+65rHJumtAlKaN0cD4RELxKqGJogswHsDMRgK93P30\nxOuLgHnAdEKbzAjgUEKikHQzod/ExcDjQCWhU+eICPGJiIhInmScSLj7hMScEdcRmh9qgCPdfUli\nl55A75RDOhLmnegFrAbeAA5z9xdSypxsZsOB6xOPOcBx7v525j+SiIiI5EukzpbuPhYY28R7ZzR4\n/Tvgd2mU+RDwUJR4REREJB5aa0MyUllZGXcIbY6uef7pmuefrnnxynhmy0JhZgOBqVOnTlUHHRER\nkQxMmzaNiooKgAp3n5ZNWaqREBERkciUSIiIiEhkSiREREQkMiUSIiIiEpkSCREREYlMiYSIiIhE\npkRCREREIlMiISIiIpEpkRAREZHIlEiIiIhIZEokREREJDIlEiIiIhKZEgkRERGJTImEiIiIRKZE\nQkRERCJTIiEiIiKRKZEQERGRyJRIiIiISGRKJERERCQyJRIiIiIl4pNPYOFCcM/fOZVIiIiIlIh7\n7oE+fYogkTCz881snpmtMbMpZjYozeMOMLP1ZjatwfbTzazezDYmnuvNbHWU2ERERNqq2bPha1+D\ndnmsJsj4VGZ2CjAKuAYoB2qBiWbWvYXjugH3AE83sUsd0DPlsUumsYmIiLRlc+aERCKfouQsVcA4\nd7/X3WcC5wKrgTNbOO524H5gShPvu7svcffFiceSCLGJiIi0WQWfSJhZB6ACeCa5zd2dUMswpJnj\nzgD6ANc2U3xXM3vPzOab2cNmVpZJbCIiIm3Z2rUwfz707Zvf82ZaI9EdaA8sarB9EaE54gvM7GvA\nDcCp7l7fRLmzCDUaxwKnJuJ62cx6ZRifiIhImzR3buhkme8aiS1as3Aza0dozrjG3ecmNzfcz92n\nkNLkYWaTgRnAOYS+GCIiItKMOXPCc6EnEkuBjUCPBtt7AAsb2X9rYD9ggJndltjWDjAzWwcc4e7P\nNzzI3TeY2evA7i0FVFVVRbdu3TbbVllZSWVlZUuHioiIlIzZs2HrraFHg2/o6upqqqurN9tWV1eX\ns/OaZzjY1MymAK+4+0WJ1wbMB37v7r9rsK8BezUo4nzgUOB7wHvuvqaRc7QDpgOPu/ulTcQxEJg6\ndepUBg4cmNHPICIiUmpGjIBp02Dq1Jb3nTZtGhUVFQAV7j6tpf2bE6VpYzQw3symAq8SRnF0AcYD\nmNlIoJe7n57oiPl26sFmthhY6+4zUrZdTWjaeAfYFrgc+ApwZ4T4RERE2pzkHBL5lnEi4e4TEnNG\nXEdo0qgBjkwZrtkT6J1hsdsBf0wc+wkwFRiSGF4qIiIiLZgzBw4+OP/njdTZ0t3HAmObeO+MFo69\nlgbDQN39YuDiKLGIiIi0dZ99Bh9/HE+NhNbaEBERKXJxjdgAJRIiIiJFL5lI5HsyKlAiISIiUvTm\nzIEvfSk88k2JhIiISJGbPTue2ghQIiEiIlL04lisK0mJhIiISJFTIiEiIiKRfPIJLF2qpg0RERGJ\nIM6hn6BEQkREpKgpkRAREZHIZs+Gnj3Dyp9xKPpEYs0X1g4VERFpO+LsaAklkEjMnRt3BCIiIvFR\nIpGl2bPjjkBERCQe7vFORgVKJERERIrWkiWwcqVqJLIya1bcEYiIiMQj7hEbUAKJxOzZUF8fdxQi\nIiL5l6yV3333+GIo+kRi7Vp49924oxAREcm/OXOgd2/Ycsv4Yij6RAKgpibuCERERPJv9ux4mzWg\nBBKJ7t2VSIiISNs0Z068IzagBBKJvn2VSIiISNtTXw/vvKMaiaztsQfU1sYdhYiISH4tWACrVyuR\nyFrfvvDhh2EJVRERkbYiOfRTTRtZSl5A1UqIiEhbMmcOtGsHffrEG0fRJxK9e0OXLuonISIibcvs\n2SGJ6Ngx3jgiJRJmdr6ZzTOzNWY2xcwGpXncAWa23symNfLeSWY2I1FmrZkdnU6Z7dtDv36qkRAR\nkbYl7sW6kjJOJMzsFGAUcA1QDtQCE82sewvHdQPuAZ5u5L39gQeAO4ABwCPAw2ZWlk5M/furRkJE\nRNqWok0kgCpgnLvf6+4zgXOB1cCZLRx3O3A/MKWR9y4EnnD30e4+y91/AUwDLkgnoAEDYMYM+Pzz\ntH8GERGRorVxI8ydG39HS8gwkTCzDkAF8Exym7s7oZZhSDPHnQH0Aa5tYpchfLGmYmJzZaYaMAA2\nbIC3305nbxERkeI2fz6sW1ecNRLdgfbAogbbFwE9GzvAzL4G3ACc6u5NLa/VM5MyG9p3XzBT84aI\niLQNhbDqZ9IWrVm4mbUjNGdc4+5zk5tzeY6qqiq6detGly7wq1/BP/4BlZWVVFZW5vI0IiIiBWP2\nbOjQAXbZpeV9q6urqa6u3mxbXV1dzmLJNJFYCmwEejTY3gNY2Mj+WwP7AQPM7LbEtnaAmdk64Ah3\nfz5xbLplbmbMmDEMHDiQU06BhQvh0UfT/llERESK0pw58NWvhpGLLWns5nratGlUVFTkJJaMmjbc\nfT0wFTgsuc3MLPH65UYOWQnsQxiJ0T/xuB2Ymfj3K4n9JqeWmXB4YntaBgwIQ0Dd0z1CRESkOBXK\niA2I1rQxGhhvZlOBVwmjOLoA4wHMbCTQy91PT3TE3KwLpJktBta6+4yUzTcDz5vZxcDjQCWhU+eI\ndIPq3x/q6uD992HXXSP8VCIiIkVi9mz47nfjjiLIePinu08ALgWuA14H+gFHuvuSxC49gd4ZljkZ\nGA6cDdQAJwDHuXva4zAGDAjP6nApIiKlbP16eO+9wqmRiDSzpbuPdfdd3X1Ldx/i7q+lvHeGuw9r\n5thr3X1gI9sfcvc9E2X2c/eJmcS0006www6a4VJERErbvHlhHomiTiQKkVmolVCNhIiIlLLZs8Nz\nIUxGBSWUSICmyhYRkdI3Zw5suSX06hV3JEFJJRIDBoR2oxUr4o5ERESkdcyeDbvvHpYQLwQFEkZu\nJDtcvvFGvHGIiIi0ljlzCqdZA0oskdhjD+jUSc0bIiJSugppDgkosURiiy1gn32USIiISGlasyYs\n2KVEohUlZ7gUEREpNXMTq1apaaMVDRgAb70VJuwQEREpJYW06mdSySUS/fuHNdpnzow7EhERkdya\nPRu22QZ23DHuSDYpuUSiX7/wrH4SIiJSapIdLc3ijmSTkkskunWD3XZTPwkRESk9hTZiA0owkQDN\ncCkiIqVp9uzC6mgJJZpIJNfccI87EhERkdz49FNYuFA1EnkxYAAsWwYffRR3JCIiIrnxzjvhWYlE\nHiSnylY/CRERKRXJVT+VSORB796w7bbqJyEiIqVjzhzYfnv40pfijmRzJZlImG3qJyEiIlIKCnHE\nBpRoIgGaKltERErLjBlhccpCU7KJRP/+oWPKp5/GHYmIiEh2NmyAN98M322FpmQTiQEDwvDPN9+M\nOxIREZHszJkDa9duGkxQSEo2kSgrgw4d1E9CRESKX7KpXjUSedSxI+y5J7zxRtyRiIiIZKemJoxI\nLLQRG1DCiQTA3nvD22/HHYWIiEh2amsLszYCIiYSZna+mc0zszVmNsXMBjWz7wFmNsnMlprZajOb\nYWY/bbDP6WZWb2YbE8/1ZrY6Smyp9t4bpk/XVNkiIlLcamoKs38EwBaZHmBmpwCjgLOBV4EqYKKZ\n9XX3pY0csgq4BXgj8e8DgT+a2WfufmfKfnVAXyC5OGrWX/9lZbB8OSxZUlhrt4uIiKRr0aKwxkYp\n1UhUAePc/V53nwmcC6wGzmxsZ3evcfe/uvsMd5/v7g8AE4GDvrirL3H3xYnHkgixbaasLDyreUNE\nRIpVsqNlodZIZJRImFkHoAJ4JrnN3R14GhiSZhnliX2fb/BWVzN7z8zmm9nDZlaWSWyN+epXw8gN\nJRIiIlKsamqga1fYbbe4I2lcpjUS3YH2wKIG2xcBPZs70Mw+MLO1hOaQ29z97pS3ZxFqNI4FTk3E\n9bKZ9cowvs106BDWbZ8+PZtSRERE4lNbC/36QbsCHR6RcR+JLBwIdAUGA78xs3fc/a8A7j4FmJLc\n0cwmAzOAc4BrsjlpWZlqJEREpHjV1MDBB8cdRdMyTSSWAhuBHg229wAWNnegu7+f+Od0M+sJ/BL4\naxP7bjCz14HdWwqoqqqKbt26bbatsrKSyspKIIzcGDu2pVJEREQKz5o1MGsWXHRR9DKqq6uprq7e\nbFtdXV2WkW2SUSLh7uvNbCpwGPAogJlZ4vXvMyiqPdCpqTfNrB2wL/B4SwWNGTOGgQMHNvl+WRks\nXgxLl0L37hlEKCIiErPp02HjxuxGbKTeXCdNmzaNioqKLKMLojRtjAbGJxKK5PDPLsB4ADMbCfRy\n99MTr88D5gMzE8cfDFwC3JQs0MyuJjRtvANsC1wOfAVIHR4aSXLkxowZcFDDcSIiIiIFrLY29I3Y\nd9+4I2laxomEu08ws+7AdYQmjRrgyJThmj2B3imHtANGArsCG4C5wGXu/seUfbYD/pg49hNgKjAk\nMbw0K1/7GrRvH7I6JRIiIlJMamrCoIEuXeKOpGmROlu6+1ig0Z4H7n5Gg9e3Are2UN7FwMVRYmlJ\nx44hmVCHSxERKTaFPDV2UoEOJsktjdwQEZFi4x4SiUKdiCqpTSQSWrxLRESKzXvvwcqVqpEoCGVl\n8PHH8MkncUciIiKSnpqa8KwaiQKQOnJDRESkGNTWwg47QM9m542OX5tIJPr2DcNnNFW2iIgUi+TS\n4WYt7xunNpFIdO4cFvBSPwkRESkWNTWF3z8C2kgiARq5ISIixWPFCnj//cLvHwFtKJHQyA0RESkW\ntbXhWTUSBaSsDD78MAylERERKWS1tdCpE+yxR9yRtKxNJRKgWgkRESl8NTWwzz7QoUPckbSszSQS\ne+wRer4qkRARkUJXDFNjJ7WZRKJLF+jTR4mEiIgUtvXr4a23iqOjJbShRALU4VJERArfrFmwbp1q\nJAqShoCKiEihS06NrUSiAJWVhXG5n30WdyQiIiKNq62FXXeFbt3ijiQ9bS6RAK25ISIihSs5NXax\naFOJxJ57hmc1b4iISCFyL64RG9DGEomuXWGXXZRIiIhIYfr4Y1iyRDUSBU0jN0REpFAVW0dLaIOJ\nhEZuiIhIoaqthW22CZ0ti0WbTCTmzYPVq+OOREREZHPJjpZmcUeSvjaZSLjDzJlxRyIiIrK5Yuto\nCW0wkdhrr/Cs5g0RESkkq1bB7NnF1dES2mAisc02sPPOSiRERKSwvPVWqDFvEzUSZna+mc0zszVm\nNsXMBjWz7wFmNsnMlprZajObYWY/bWS/kxLvrTGzWjM7Okps6dDIDRERKTQ1NdC+ffiOKiYZJxJm\ndgowCrgGKAdqgYlm1r2JQ1YBtwAHAXsCvwJ+bWZnpZS5P/AAcAcwAHgEeNjMyjKNLx1lZTB9emuU\nLCIiEk1tbZg4sXPnuCPJTJQaiSpgnLvf6+4zgXOB1cCZje3s7jXu/ld3n+Hu8939AWAiIbFIuhB4\nwt1Hu/ssd/8FMA24IEJ8LSorg3ffhTVrWqN0ERGRzBXb1NhJGSUSZtYBqACeSW5zdweeBoakWUZ5\nYt/nUzYPSZSRamK6ZWaqrAzq60OnFhERkbjV18MbbxRf/wjIvEaiO9AeWNRg+yKgZ3MHmtkHZrYW\neBW4zd3vTnm7Z5Qyo9LIDRERKSRz54ZRG8VYI7FFHs91INAVGAz8xszecfe/ZltoVVUV3RqstVpZ\nWUllZWWTx2y3Hey0kxIJEREpDLW14bk1aiSqq6uprq7ebFtdXV3Oys80kVgKbAR6NNjeA1jY3IHu\n/n7in9PNrCfwSyCZSCyMUibAmDFjGDhwYEu7fcHee6vDpYiIFIaamnCDu+OOuS+7sZvradOmUVFR\nkZPyM2racPf1wFTgsOQ2M7PE65czKKo90Cnl9eTUMhMOT2xvFVpzQ0RECkWxdrSEaKM2RgMjzOyH\nZrYncDvQBRgPYGYjzeye5M5mdp6ZfdvMdk88/ge4BPhzSpk3A0eZ2cVmtoeZ/ZLQqfPWSD9VGsrK\n4J134PPPW+sMIiIi6SnGqbGTMu4j4e4TEnNGXEdofqgBjnT3JYldegK9Uw5pB4wEdgU2AHOBy9z9\njyllTjaz4cD1iccc4Dh3b7U6g7Iy2LgR5syBffZprbOIiIg0b/Fi+PDD4q2RiNTZ0t3HAmObeO+M\nBq9vJY2aBXd/CHgoSjxRlCWmunr7bSUSIiISn5deCs8HHBBvHFG1ubU2krbfHnr0UD8JERGJ16RJ\nsMsuYR2oYtRmEwnQVNkiIhK/SZPgwAPjjiK6Np9IqEZCRETismoVTJumRKJolZWFabLXr487EhER\naYtefRU2bFAiUbTKysJ/4DvvxB2JiIi0RZMmwbbbbhoAUIzafCIBat4QEZF4TJoURmu0K+Jv4yIO\nPXs77gjdu6vDpYiI5N+GDfDyy8U77DOpTScSoA6XIiISjzffhM8+K+7+EaBEgrKy8J8pIiKSTy+9\nBB07wqBBcUeSnTafSOy/f6iRWLo07khERKQtmTQJ9tsPOneOO5LstPlE4tBDw/Pzz8cahoiItCHu\n8OKLxd+sAUok2Hln6NsXnn027khERKSteP99WLBAiUTJOPRQeO65uKMQEZG2YtKk8Lz//vHGkQtK\nJIBhw2DmzJAdioiItLZJk0Jn/+23jzuS7CmRAA45JDyrVkJERPKh2BfqSqVEgjAx1b77qp+EiIi0\nvuXLw0SIxT4RVZISiYRhw5RIiIhI65s8OTyrRqLEDBsG770H8+bFHYmIiJSySZNgp52gT5+4I8kN\nJRIJQ4eGRVNUKyEiIq0p2T/CLO5IckOJRMK228LAgUokRESk9axdC6++WjrNGqBEYjPDhoWRG+5x\nRyIiIqVo6lRYt06JRMkaNgw+/hhmzYo7EhERKUWTJkHXrtCvX9yR5I4SiRQHHghbbKHmDRERaR2T\nJsGQIeG7plRESiTM7Hwzm2dma8xsipk1uQiqmR1vZk+a2WIzqzOzl83siAb7nG5m9Wa2MfFcb2ar\no8SWja22gsGDlUiIiEju1deHpcNLqVkDIiQSZnYKMAq4BigHaoGJZta9iUOGAk8CRwMDgeeAx8ys\nf4P96oCeKY9dMo0tF5L9JOrr4zi7iIiUqhkz4JNPSmciqqQoNRJVwDh3v9fdZwLnAquBMxvb2d2r\n3P1Gd5/q7nPd/UpgDvCdL+7qS9x9ceKxJEJsWRs2LMw69sYbcZxdRERK1UsvQfv28I1vxB1JbmWU\nSJhZB6ACeCa5zd0deBoYkmYZBmwNLG/wVlcze8/M5pvZw2ZWlklsuTJ4MHTurOYNERHJrUmToLw8\ndLYsJZnWSHQH2gOLGmxfRGiOSMdlwFbAhJRtswg1GscCpybietnMemUYX9Y6dQrVTkokREQkl0pp\noa5UeR21YWbDgauBk9x9aXK7u09x9/vc/Q13fxE4AVgCnJPP+JKGDYMXXoANG+I4u4iIlJqPPgpL\nMJRiIpHpAJSlwEagR4PtPYCFzR1oZt8H/gic6O7NLtjt7hvM7HVg95YCqqqqolu3bpttq6yspLKy\nsqVDmzRsGFx5ZZg4pNTaskREJP9eeik8x9HRsrq6murq6s221dXV5az8jBIJd19vZlOBw4BH4f/6\nPBwG/L6p48ysErgTOMXd/93SecysHbAv8HhL+44ZM4aBAwem9wOkab/9YOutQ/OGEgkREcnWpEmw\n++7QM91OADnU2M31tGnTqKioyEn5UZo2RgMjzOyHZrYncDvQBRgPYGYjzeye5M6J5ox7gEuA/5pZ\nj8Rjm5R9rjazw82sj5mVA/cDXyEkH3m3xRZhES/1kxARkVwo1f4RECGRcPcJwKXAdcDrQD/gyJTh\nmj2B3imHjCB00LwNWJDyuClln+0IzR5vE2ohugJDEsNLYzFsWPiP//zzuCIQEZFSsHIl1NaWbiIR\naZJOdx8LjG3ivTMavD40jfIuBi6OEktrOfTQsErblClw8MFxRyMiIsVqypQwyWGpTUSVpLU2mtC/\nP2y3nZo3REQkOy+9BNtvD3vsEXckrUOJRBPatQu1EkokREQkG8n+EWZxR9I6lEg0Y9gweOUVWLUq\n7khERKQYrV8fmjZKtX8EKJFo1rBh4UOQHP8rIiKSiddfh9WrlUi0WXvuGcb8qnlDREQy5Q7XXgu9\nekGOpzsqKJFGbbQVZqFWIpeJxIQJYZ6KE07IXZkiIlJ47r8f/vUveOQR6Ngx7mhaj2okWnDooWGq\n7BUrsiunvj5Mu33KKXDFFbmJTURECtPChXDhhTB8OBx7bNzRtC4lEi0YNiwkAS+8EL2Mzz+H006D\nG24I5c2ZA59+mrsYRUSksFxwQah9vvnmuCNpfUokWtCnD+yyS/TmjWXL4PDD4e9/D80ao0eHdrPa\n2tzGKSIiheFvf4OHHoLbboPu3eOOpvWpj0QLsuknMXcuHHMMLF8ejt9//zAKpGPH0JO3lHvxioi0\nRUuXwvnnh35wJ54YdzT5oRqJNAwbBm++CVdfDTPTXP1jyhQYPDjUPkyeHJIIgA4dYJ99QiIhIiKl\n5ac/DTeMt91WuhNQNaREIg3f+x6cfTbccgvstVcYxvO738EHHzS+/9//Hjpp7rFHSCJ2333z98vL\nlUiIiJSaf/4zjNS46aZ4lguPixKJNGy5JYwbF3rh/v3vITH4xS/gK18JC3qNGxf6QrjDmDGhOuu4\n4+Dpp8P86g2Vl8P06bBuXf5/FhERyb0VK+Ccc+Doo+EHP4g7mvxSIpGBzp3h+ONDp8lFi+Cee6BL\nl9Ae1rMnDBoEF18MP/sZPPBA2L8x5eWh6mv69PzGLyIireOyy8JovHHj2k6TRpISiYi22QZ++EN4\n4glYsCBUZfXoAXfdBSNHhkW/mtKvX/igqXlDRKT4PfUU3Hkn3Hgj9O4ddzT5p1EbObDjjqFW4vzz\n09u/a1fo21eJhIhIsfvsMxgxInTKHzEi7mjioUQiJupwKSJS/K64ApYsCUP821qTRpKaNmJSXg41\nNWHWTBERKT4vvgi33hpmLd5tt7ijiY8SiZiUl8OqVWG6bBERKS6rV8P//E+YI+iCC+KOJl5KJGJS\nXh6e1bwhIlJ8fvELmD8f/vQnaN8+7mjipUQiJt27h969SiRERIrLlClhzqDrrgsTD7Z1SiRipA6X\nIiLFZe1aOOMMqKgI8waJEolYJRMJ97gjERGRdFx3Hbz7Ltx9d1gmXJRIxKq8PKwU99FHcUciIiIt\nmToVfvvb0D9i773jjqZwREokzOx8M5tnZmvMbIqZDWpm3+PN7EkzW2xmdWb2spkd0ch+J5nZjESZ\ntWZ2dJTYiok6XIqIFId160KTRr9+cPnlcUdTWDJOJMzsFGAUcA1QDtQCE82sexOHDAWeBI4GBgLP\nAY+ZWf+UMvcHHgDuAAYAjwAPm1lZpvEVk9694UtfUiIhIlLorr8eZswITRodOsQdTWGJUiNRBYxz\n93vdfSZwLrAaOLOxnd29yt1vdPep7j7X3a8E5gDfSdntQuAJdx/t7rPc/RfANKCkR+eaqcOliEih\nq6kJk0797/9C//4t79/WZJRImFkHoAJ4JrnN3R14GhiSZhkGbA0sT9k8JFFGqonpllnMlEiIiBSu\n9evhzDNhzz3hyivjjqYwZVoj0R1oDyxqsH0R0DPNMi4DtgImpGzrmWWZRau8HN5/H5YtizsSERFp\n6Le/hTfeCE0aHTvGHU1hyuvgFTMbDlwNHOvuS3NRZlVVFd26ddtsW2VlJZWVlbkovtUlO1zW1MBh\nh8Ubi4iIbDJ9ehjuedllsN9+cUcTXXV1NdXV1Zttq6ury1n5mSYSS4GNQI8G23sAC5s70My+D/wR\nONHdn2vw9sIoZQKMGTOGgQMHtrRbwerbF7p0Cc0bSiRERArDhg1hlMZXvwrXXBN3NNlp7OZ62rRp\nVFRU5KT8jJo23H09MBX4v6+8RJ+Hw4CXmzrOzCqBu4Dvu/u/G9llcmqZCYcntpe09u1D5x31kxAR\nKRyjRsFrr4W1NDp3jjuawhalaWM0MN7MpgKvEkZxdAHGA5jZSKCXu5+eeD088d6FwH/NLFnzsMbd\nVyb+fTPwvJldDDwOVBI6dY6IEF/RKS+H5xrW0YiISN5Nnw5XXQUPPwyXXgqDB8cdUeHLePinu08A\nLgWuA15yPF/XAAAgAElEQVQH+gFHuvuSxC49gd4ph4wgdNC8DViQ8rgppczJwHDgbKAGOAE4zt3f\nzjS+YlReDrNmhWVpRUQk/957D04/HfbdN/RZu/de+H//L+6oikOkzpbuPhYY28R7ZzR4fWiaZT4E\nPBQlnmJXXg719aFnsLJfEZH8WbQoTDZ1++1hgsBbboERIzRCIxNaa6MA7LNPWPxF/SRERPKjrg6u\nvjp0prz3XvjlL2HuXDj/fCURmdLaZQWgUycoK1MiISLS2urr4aabQi3E6tVw4YXws5+F2giJRolE\ngdAMlyIire+uu+CSS+Ccc0KNxJe/HHdExU9NGwWivBzefDNMxyoiIrn3+efw61/DKaeEPhFKInJD\nNRIForw8fMhnzAjL1IqISG796U/wwQfw78ZmM5LIVCNRIAYMCM9q3hARyb21a0O/iOHDYa+94o6m\ntCiRKBDbbBN6DyuREBHJvTvugI8/hl/8Iu5ISo8SiQKiDpciIrm3Zg3ccAP84AdhfSPJLSUSBWTg\nwDCjWn193JGIiJSO22+HJUvCKA3JPSUSBaS8HFauhHnz4o5ERKQ0rFoVprr+0Y9C87HknhKJAlJe\nHp7VvCEikht/+AMsXx4W4pLWoUSigPToATvtlNtEYskSWLEid+WJiBSLzz6D3/wGzjwTdt017mhK\nlxKJApPLDpfTp4ept08+OTfliYgUk1tvDWtqXHll3JGUNiUSBSZXicTs2XDYYaHj5lNPwcyZ2Zcp\nIlIsVq6E3/0urOT5la/EHU1pUyJRYMrLYeHCMN45qnffhWHDwiI0tbXQvXtoJxQRaSt+//vQ0fKK\nK+KOpPQpkSgw2Xa4/OCDUBOx5ZbwzDOw885w1lkwfnz4pRIRKSSrV4eFtFavzl2ZK1bAqFFhYa6d\nd85dudI4JRIFpk8f6NYtWiLx8cehJgLg2WdDx02Ac8+FTz+FBx7IXZwiItmaMwcGDw43O9dem7ty\nb745TIn985/nrkxpmhKJAmMW1t3INJFYvDjURKxdG5KI3r03vbfLLvDtb8PYseCe23hFRKJ4+GHY\nb7+wWOFZZ8GYMTBrVvblfvIJjB4N55236WZKWpcSiQKUaYfL5cvh8MPDL9Azz4RajYbOOy/Mmjll\nSu7iFBHJ1IYN8LOfwfHHh79b//0v3HJLuPm56KLsb3ZGjw7nuPzy3MQrLVMiUYAOOCB0mNxtNzj7\nbHjwwZAsNKauDo44AhYsgKefbnoe+SOOCOWNHdt6cYu0htdeg/Xr445CcmHhwpA8jBoFN94Y/rZt\nsw107hxqJCZOhEcfjV7+smVw001wwQVhXh7JDyUSBeh734PHHw/NES++GOaB6N4dBg2C//1feO65\nUB346adw9NEwd24Y4rn33k2X2a4d/PjHMGFCmKRKpBh8+CF8/esadVQKXnoprCc0c2Zofr3kktCU\nm/Sd78BRR0FVVVhkK4orrgg1GpddlpuYJT1KJAqQGRxzTBi+NGMGzJ8fejX37Ruehw2D7baDffeF\nt96CJ58M/SpacsYZIaG4667W/xlEcuGpp8IXw4QJcUciUbmHWoJDDglrXUybBkOHfnE/s9BJ8sMP\nw/wPmRo7NiwVftNN4cZL8keJRBHo3TskAfffH0Zm1NTAddfB/vuHqsBBg9IrZ/vt4fvfDyvhbdzY\nujGL5MKTT0L79uFu9qOP4o5GMrVqVfibU1UFF164+WiyxvTtCxdfDCNHwvvvp3+eZ54J5V90Uei4\nKfkVKZEws/PNbJ6ZrTGzKWbW5FeZmfU0s/vNbJaZbTSz0Y3sc7qZ1Sfer088cjiquHS0awf9+8Ol\nl4bhnEOGZHb8eeeFX9AnnsjsuI8/Dr+gCxZkdpxIVMlZWc85Bzp0gL/9Le6Iit+GDaHD9bp1+Tnf\nb34Djz0W+kKMGhX+H1ty1VVhMr1LLknvHHPmwEknhVFrN96YXbwSTcaJhJmdAowCrgHKgVpgopk1\nVZnUCVgM/AqoaaboOqBnymOXTGOTlg0aFIZcZdLpcvVqOPbY0CQyZkzrxSaS6vXXQ+e5k08OnYUf\nfDDuiIrfXXeFm48ddoDTToN//CN6f4SWuIda1OHD4cQT0z+ua9eQEDz0UOhA3pwVK0Lfih12gL/+\nFbbYIruYJZooNRJVwDh3v9fdZwLnAquBMxvb2d3fd/cqd78PWNlMue7uS9x9ceKhLoGt5Lzz4N//\nDp00W1JfDz/4Abz99qZkIpcz0DXlvvtgr71g6tTWP5cUpiefhK22Cl98J52k5o1cmDAhTABVVRWm\nzz/hhNCf4KST4C9/CR24c+XVV8Pos+HDMz/2+98P/Sh+8pOma082boTKSli0KNR6bLttdvFKdBkl\nEmbWAagAnkluc3cHngYyrGT/gq5m9p6ZzTezh82sLMvypAmnnBJ+6caNa3nfK64Idy3V1aET04oV\nrTtDZn09XH11SF4WLw53Gx980Hrnk8L11FNw6KHQsSMcd1yoFn/oobijKl5LlsB//hP6W/3yl/Dm\nm2ECqKuugnnzwpfyDjuEG4Y//zn7flT33w+9esHBB2d+rFmYW2L27PDcmMsvD5+RCROaHvYu+ZFp\njUR3oD2wqMH2RYTmiKhmEWo0jgVOTcT1spn1yqJMaUKXLuGPyV13NV+teccd8Nvfhglejj02THT1\n7W+HX+zWmCFz9epwJ3L99aFtdfr08OXxne/k9k5JCt+qVTBpUmjSgJD4qnkjO488En5vv/vdTdv6\n9g03C6+9FpKJG24IzUk//GHTX+Dp2LAhNDV8//uhs2wU/fqF2tNrr/3iIoZ/+lP4u3TTTWFeComZ\nu6f9AHYC6oFvNNj+G2ByGsc/B4xOY78tgDnAtc3sMxDwqVOnumRuzhx3cB8/vvH3n3rKfYst3H/8\nY/f6+k3bn3wyHPfCC7mNZ8EC90GD3Lt0cf/73zdtf/NN9623dv/Wt9w3bMjtOaVwPf54+JzNnLlp\n2/jx7mbuH30UX1zF7Mgj3Q85JL19Kyvdv/Y1940bo51r4sTw//faa9GOT1q+3L17d/cf/nDTthdf\ndO/Qwf3sszf/2ySZmTp1qgMODPQM8oDGHuYZ3FommjZWA99z90dTto8Hurn78S0c/xzwurtfnMa5\nJgDr3f3UJt4fCEwdOnQo3bp12+y9yspKKisrWzpFm3fUUWFa7Vde2Xz7jBmhXXrwYPjnPzfvwFRf\nD2Vl4W4hV2P7a2tDrcPGjaGtc+DAzd+fOBG+9S04//wwzlxK309/GprU3ntv06RFK1bAjjuG3v8/\n+Ums4eXdiy+G0VrbbBPt+E8+CdfuppvC71FLJk2Cgw4Kv3vJWqFM/OhHMHlymHwqddKpKO68E0aM\nCH1kevUKE5SVlYU+NB07Zld2W1FdXU11dfVm2+rq6njhhRcAKtx9WlYnyDTzAKYAN6e8NuAD4LI0\njk23RqIdMAO4sZl9VCORpUceCXcN//3vpm2LFrn36eO+zz7udXWNH3fLLe7t27t/8EH2MTz6qPtW\nW7kPHOj+4YdN7zd2bIj1llsyK3/qVPeKCvdevdyvvdZ94cLs4pX82Gsv97PO+uL2Y45xP+ig/MfT\n0DnnNF2bl2t//3v47J99dvQyMq3Nqa9333df9+OOy/xcq1eHWsRf/jLzYxuzcaP7fvu5Dxjg3q9f\n+Pu0ZEluym7LclkjESWROJlQK/FDYE9gHLAM2CHx/kjgngbH9AcGAP8F/px4vVfK+1cDhwN9CENK\nq4FVwJ7NxKFEIksbNrh/5SvuZ54ZXq9Z4z5kiPuOO7q/917Tx9XVuXft6n7VVdHPXV/vfuON4Y/b\n8ce7f/ZZy8dUVbm3axeqvVuyerX75ZeHhKdfP/cRI9y33NK9Y0f3008PCYYUpvnzw1+mCRO++N7d\nd4fPzIIFeQ/r/3z4YYhv553d161r3XNNnx5+1770pZBwN5Xct+Tb33Y/4IDMjrn99vD79v77mR03\nYUK4PrNnZ3Zcc6ZMCWV27er+1lu5K7ctizWR8PAlfh7wHrAGmAzsl/Le3cCzDfavBzY2eLyb8v5o\nYF6ivAXAY0C/FmJQIpED11/v3rmz+9Kl7t//fvj3lCktH3fBBSHhWLs283OuWxe+2MH9Zz9Lvx12\nwwb373wn/DGprW16v+efD+27nTqFny/5x37ZMvff/jYkT+B+4IHhj9769Zn/DG3V2rUh8XvqqdY7\nx113hWRh2bIvvrd8eWgfz7RmKpd+//uQoIL7X/7Seuf55JPwOd5779BXpH1791tvzbycurqQQI8Z\nk9lxn34aahauvDKz47773dDfKdfuuCP3fbPastgTiUJ4KJHIjYULwx/mfv3Cp+HBB9M7bsaMsP+f\n/5zZ+err3U89NZzz7rszDtc//TRUcfbu/cW70hUrQpVzMklI7aiXav1694cecj/44E13liNHhmSq\nVNXXhxqnbF1+ebhmp5+efVlNOeUU969/ven3jz463uaNoUNDDMOGuQ8e3Drn2LgxdDDedtvQMdrd\n/YQT3MvKMu9geP/94f8s05oF98xvGJYvj5a0SP4pkVAikVOVleGTMHJkZscdfnjzf/Abc/vt4VwP\nPJDZcak++CD0edhvP/dVq8K2Rx91//KXQ23FbbelX8vx+uuhaadTpzBKZfBg95//3P2JJ9xXrowe\nY6E599yQMGXTJPDii6GmYKedQjt1a9iwIVTjN3cX/Kc/xde8sWBBOPddd4XPHLhPnpz781x1VTjP\nE09s2vb00+F8//lPZmUdf3zmv6dJ06dn9vt6553xNz1JepRIKJHIqQ8+CH+cM73TSf4hfeWV9Paf\nOjXcrfz4x5nH2NC0aWGo6HHHhTtYCB3xotx1ubsvXhw6dJ58snuPHqG89u1DFe2ll7r/85+hxqMY\nvfJK+Hk6dgw1NVHa9VeuDMnDAQdsagOfPz/3sf73vy1/WcbZvHHbbSHhXLYsJKu77x6aBHPpoYfC\nNbjhhs2319e79+0bPu/p+vTT0Fz5299Gj+fQQ8PnJh3Dhrkfdlj0c0n+KJFQIlEQNmxw33VX9x/8\noOV9V6xw3223MDojF1Xs7mHUiZn79tuH6ttcjSmvrw/NIuPGuQ8fHmo6IHQ8q6jYfJ6LQrdxY7gb\n7d8/fDlvsYX7hRdmXs7ZZ4fOfu+8E5IuCNc8166/PtQqff558/sdfXRoYsi3Qw91P+KITa9zOYLJ\nfVPnyhNPbPzzPGZMSKI+/ji98pJJ39y50WN68MFQRnP9ktzDiJBkbY0UPiUSSiQKxu9+F+50Fy1q\nep/6+tC+261bdn/QGvPqq60/FKy+PnyB3nWX+1FHhd+a4cOLo0/F3Xf7ZhOI3XJLeH3ffemXkZwc\naty4Tdv22iv0R8m1gw8OHWpbEkfzxqJFIZn84x83bfv00/C5/vnPsy8/2blyn31CuY1ZvjyMPvr1\nr9Mr8+ST3cvLs4tr3brQnNXS//eoUaGJsFhr7toaJRJKJArGsmXhD9uvftX0PjffHD5pxXQn35T6\n+vAlvN127j17hlqRQrViRWimqazctK2+PtQgbbmle01Ny2UsXRp+zqOP3vwO+ZxzQjKRSytXpt9k\nsWxZqF2JMoohqnHjQu3D4sWbb7/kkvB5SPbXiWLDhtA0t+22IWltzplnhpFHLc30unp1qEW6/vro\ncSVdc00oq7kkoaIi3DBIcVAioUSioJx1Vqj+b6ztfcqU8OVQVZX/uFrTRx+FsfkQpu9dvjzuiL7o\n4ovDH/+GE32tWhWaOnbbrfm46+vdTzopdH5seOefHAnQXE1Uph57LJQ5a1Z6+x91VKjByJdvfrPx\n9v9580JNxR/+EL3sK6/8YufKpiT7kTz6aPP7/eMfmV3P5nz0UUiimkryZs4M5/rb37I/l+SHEgkl\nEgWlpsYbnUBo2bJw5zR4cMtt3sWovj7MGNitWxhFks5EWfny9tvhjr1hh72kuXPDXfQxxzQ9wuWB\nB8L/61//+sX3kpNG5fKL48IL3XfZJf2+Lsn5JtLtL5CNJUvCF2lTycL3vue+557R1qZIdq7MZNTU\noEEhkWrOaaeFZpJcOfHEUAvV2P/PNde4b7NNqAWR4qBEQolEwTnooM07v23cGO7Yv/Sl6CMpisUH\nH4QFkSBUO8fdRlxfHzoEfvWrzY//f+KJ8EV8zTVffO/DD0M1e3MjEvr0idZxsyl77hkmKktXsnnj\ntttyF0NTksMam5pi/cUXw/9/OjUKqWpqQufKk07KrLNwso9IU80ga9eGL/ZcTVPt7v7ss+FnfPbZ\nzbfX14e+HT/6Ue7OJa0vtkW7Ckly0a6pU6cysOEqT5J3Dz4IJ58cFuDq1y8sP/6zn8Hjj8Mxx8Qd\nXetzD8uyX3wxdOsWFi2CsBDZxo1hsbOG/+7QAb78Zdh5502PL38ZOnXKLpZHHglLRT/2WFj2vTm/\n/jVcffXm+7qHBd3efBPeegu+9KXGj/3Rj8L/9+uvZxcvwPz5sMsu4XN04onpH3fUUfD55/Dcc9nH\n0Jyjj4Y1a+D55xt/3x0GDYLu3eHf/06vzHffhQMOCAtR/ec/0LVr+vGsXh0+L2edFX7XGnr88fD/\n+dZbsPfe6ZfbHPdQ1t57b76c+3//GxbSevJJLeldTKZNm0ZFRQXEsWhXoTxQjURBWbcu9JM466ww\nQqB9e/crrog7qvx7//1QE7PzzqFZp0+fUDPQt2+oFt577zCL6IAB4XW3buEuL/Wx445hmOyxx4Y7\n/nffTf/8q1eHczbsHNmUjRvDKIlu3TbNoJhcIK2lu+vkXfEnn6QfX1PuvDP0M2hsWuyWjmvt5o3l\ny0PNR0udQP/853Ddpk9vucyFC8PnYvfdo/czqaoKQ58bG059xhnue+yR+2W2k8NdUxf/+ulPQ6fe\nljp/SmFR04YSiYL0q1+F0QC9eoVmDq1hkZ6VK0OfhiefDF/O110X5m045pgwYmKrrULbfDpfCr/6\nVejcmkkHuxUrQtX0vvuGqvYuXcJMmC15553wF+Sf/0z/XE05+WT3b3wj8+OWLm395o3kQmEtrZz5\n+efh/6ulYZIrVoREcqedMksSG5o1K1z/e+/dfPu6daH/S6ZrZKRjxYrweUw2h23YEH7miy7K/bmk\ndSmRUCJRkBYtCnNK7Lhj+ssVS/Pq6jYtcHb44c3PJjl/fkjkLr888/O8+Wb4gujUKdwpNzWPQar6\n+pA0RjlfquS02FdfHe34I490P+SQ7GJozre+lf7Mjslkuqk5RtasCbFuu637G29kH9s3vxlW7E01\ncWL4vLz+evblN+acc0IStG7dpmm7053dVgpHLhOJdlm1i4ik2HFH+MtfYOLE0O4r2dtmG/jjH+GJ\nJ+Dtt2GffeDuu0MjSEOXXhr6Z1x1VebnSZbbuTPcc0967fVmMHQovPBC5udLNW0aLF8ORxwR7fiT\nTgoxLFqUXRyNqasLbf/p9ts455zQB+aOO7743saNMHw4TJkS+qTsu2/28f34xzB5MtTUbNr20EOw\n227Qv3/25TfmvPPg449DX5wHHoCvfjX0D5G2S4mE5NTxx8OAAXFHUXqOOip0nPvud+HMM+HYY8Mf\n86Tnn4cJE0LHu623jnaOk06CZctCB8B0DR0Kr70Gq1ZFOyeEL+qtt4ZvfCPa8d/9bkhqmkqwsvHY\nY7B+PZxwQnr777ADnHYa3HprOC7JPXzpP/po+H868MDcxHfssSFp/8MfwusNG+Af/wiJj1luztFQ\nv34h/jFjQtIyfHjrnUuKgxIJkSKx7bahtuCRR0JP+b33DneEGzbAhRfCkCHhSywb7dtntv/QoeH8\nU6ZEP+eTT8Khh4ZRLFFsv3344rziCth1V7jsspDc5CKpePBBGDwYevdO/5iLLoKPPgpfsklXXx1q\nKe66C77znezjStpiCzj7bLj//lB78uKLsGRJZiNfojjvPHj55XDO4cNb91xS+JRIiBSZY4+F6dPh\nyCPh1FOhvDzUVtxyS/7vDPfaK3yR/+c/0Y7/9NPwhRS1WSPp/vvDENBjjgnJ1qBBsPvuIbmoqYmW\nVKxcGZrpTjops+P23RcOOwxuuim8vvlmuP56+N3v4PTTM4+jJSNGwNq18Oc/h+TlK1+B/fbL/XlS\nnXBCaMocOBD23LN1zyWFT4mESBHafnuorg53zAsXhmrzMCQ8v9q1g4MOit5P4vnnQ41GtolE+/Zw\nyCGhin/BglDLMWxY6F9SXh6+7H7xC5g9O/0yH388zFHxve9lHs9PfwqvvAKXXBL+fdlloQ9La+jV\nKzTv/OEP8Pe/h3hbO6Hs1An++le4/fbWPY8UB01IJVLk1q8PX6TtYrotGDMm3PnX1WU+mdZPfgL/\n/GeYnKk1vvzWr4ennw5fev/4R5hU6rbbwl18S773PfjgA3j11czPW18Pe+wB77wTJu76059a98v9\n2WdDLQjASy/B/vu33rmkNORyQirVSIgUuQ4d4ksiIPST+Pzz0G8jU08+GWojWutLtkOHMCvl+PFh\nVMdZZ4U+BT/5yeadIRv67DP417+i9zVo1y40bVx4Yegb0do1BIceGhKXXr1Cnw6RfNoi7gBEpLj1\n7x9GXbzwQmajEd5/PzQ13HBD68WWqnNnGDs29GG48MIwnHbChNBM1NC//hX6HWTTafFb3wqPfDAL\ntR6ffhpvUiltkz5yIpKVLbYIQ0Yz7Scxfjx07Bj6MuTTj38cmjtqa8MaEdOnf3Gfv/0tdCTcbbf8\nxpaN/fcPHXBF8k2JhIhk7eCDQ9v8hg3p7b9wYRjF8JOfwHbbtW5sjTn44NAUs9VWoSngscc2vbd6\ndeho2dpDKEVKhRIJEcna0KGhX0HqDIvNueaaUBtx5ZWtG1dz+vQJQ08PPxyOOw5GjgzDRP/975BM\nRBmtIdIWqY+EiGRtv/1CH4QXXmh5DoPp0+HOO2HUqHhqI1J17RqaMa69Fv73f+GNN0LfiH79oG/f\neGMTKRaRaiTM7Hwzm2dma8xsipk1OdO6mfU0s/vNbJaZbTSz0U3sd5KZzUiUWWtmR0eJTUTyr2PH\nMLNmOv0kLr881Aacd17rx5WOdu1CIvHgg2EK64cfznwSKpG2LONEwsxOAUYB1wDlQC0w0cy6N3FI\nJ2Ax8Cug0YpPM9sfeAC4AxgAPAI8bGZlmcYnIvEYOjRM0Vxf3/Q+zzwTRkSMHBmSj0Jy4omhn8dx\nx7XODJQipSpKjUQVMM7d73X3mcC5wGrgzMZ2dvf33b3K3e8DVjZR5oXAE+4+2t1nufsvgGnABRHi\nE5EYDB0aVvFsbBQEhATj0ktD58ZC7cg4YECokchkbQ2Rti6jRMLMOgAVwDPJbR6mxnwaGJJFHEMS\nZaSamGWZIpJHgweHoaBNNW/cd1/ojDlqlFaLFCklmdZIdAfaA4sabF8E9Mwijp6tUKaI5FGXLmGx\nrMYSiTVrwgiNE0/U9M0ipaboR21UVVXRrVu3zbZVVlZSWVkZU0QibdfQoWH1TffNax1uuilMUT1y\nZHyxibRV1dXVVFdXb7atrq4uZ+VnmkgsBTYCPRps7wEszCKOhVHLHDNmjBbtEikQQ4fCb34TFqv6\n2tfCtsWLQwJx3nlhaW8Rya/Gbq5TFu3KWkZNG+6+HpgKHJbcZmaWeP1yFnFMTi0z4fDEdhEpEgcc\nEIZTpjZvXHtt2Hb11fHFJSKtJ8qojdHACDP7oZntCdwOdAHGA5jZSDO7J/UAM+tvZgOArsAOidd7\npexyM3CUmV1sZnuY2S8JnTpvjRCfiMSkW7cw8iGZSMycCePGwVVXNb44logUv4z7SLj7hMScEdcR\nmh9qgCPdfUlil55Aw8FTrwOe+PdAYDjwPrBboszJZjYcuD7xmAMc5+5vZxqfiMRr6NAwhBLg5z8P\nQykv0EBukZIVqbOlu48Fxjbx3hmNbGux5sPdHwIeihKPiBSOoUND58r77oNHHoEHHgjTZ4tIadKi\nXSKSUwceGJ5HjAjDQU85Jd54RKR1KZEQkZzaYQcoKwuLX914Y+hoKSKlq+jnkRCRwnPmmTB3bmjm\nEJHSpkRCRHLukkvijkBE8kWVjiIiIhKZEgkRERGJTImEiIiIRKZEQkRERCJTIiEiIiKRKZEQERGR\nyJRIiIiISGRKJERERCQyJRIiIiISmRIJERERiUyJhIiIiESmREJEREQiUyIhIiIikSmREBERkciU\nSIiIiEhkSiREREQkMiUSIiIiEpkSCREREYlMiYRkpLq6Ou4Q2hxd8/zTNc8/XfPiFSmRMLPzzWye\nma0xsylmNqiF/Q8xs6lmttbMZpvZ6Q3eP93M6s1sY+K53sxWR4lNWpd+2fNP1zz/dM3zT9e8eGWc\nSJjZKcAo4BqgHKgFJppZ9yb23xX4J/AM0B+4GbjTzA5vsGsd0DPlsUumsYmIiEh+RamRqALGufu9\n7j4TOBdYDZzZxP4/Bt5198vdfZa73wb8LVFOKnf3Je6+OPFYEiE2ERERyaOMEgkz6wBUEGoXgPDt\nDzwNDGnisMGJ91NNbGT/rmb2npnNN7OHzawsk9hEREQk/7bIcP/uQHtgUYPti4A9mjimZxP7b2Nm\nndz9c2AWoUbjDaAbcBnwspmVufuCJsrtDDBjxowMfwTJRl1dHdOmTYs7jDZF1zz/dM3zT9c8v1K+\nOztnW1amiUSrcPcpwJTkazObDMwAziH0xWjMrgCnnXZaa4cnDVRUVMQdQpuja55/uub5p2sei12B\nl7MpINNEYimwEejRYHsPYGETxyxsYv+VidqIL3D3DWb2OrB7M7FMBE4F3gPWNh+2iIiIpOhMSCIm\nZltQRomEu683s6nAYcCjAGZmide/b+KwycDRDbYdkdjeKDNrB+wLPN5MLMuAB9IOXkRERFJlVROR\nFGXUxmhghJn90Mz2BG4HugDjAcxspJndk7L/7cBuZvYbM9vDzM4DTkyUQ+KYq83scDPrY2blwP3A\nV4A7I/1UIiIikhcZ95Fw9wmJOSOuIzRR1ABHpgzX7An0Ttn/PTP7FjAGuBD4EPgfd08dybEd8MfE\nsQYNLQkAAAQZSURBVJ8AU4EhieGlIiIiUqAsjN4UERERyZzW2hAREZHIlEiIiIhIZEWZSGS6aJik\nz8wOMrNHzeyjxOJpxzayz3VmtsDMVpvZU2bW3DBdaYGZXWFmr5rZSjNbZGb/MLO+jeyn654jZnau\nmdWaWV3i8bKZHdVgH13vVmJmP0/8fRndYLuueQ6Z2TUpC2EmH2832Cfra150iUSmi4ZJxrYidKA9\nD/hCBxoz+xlwAXA28HVgFeH6d8xnkCXmIOAW4BvAN4EOwJNmtmVyB133nPsA+BkwkDDt/7PAI2a2\nF+h6t6bEjd/ZhL/dqdt1zVvHW4SBEckFMQ9MvpGza+7uRfUgzIB5c8prI4wEuTzu2ErtAdQDxzbY\ntgCoSnm9DbAGODnueEvlQZiKvh44UNc9r9d9GXCGrnerXuOuhCURhgHPAaNT3tM1z/31vgaY1sz7\nObnmRVUjEXHRMMkRM+tDyGhTr/9K4BV0/XNpW0Jt0HLQdW9tZtbOzL5PmA/nZV3vVnUb8Ji7P5u6\nUde8VX0t0VQ918zuM7PekNtrXhBrbWQgyqJhkjs9CV9wjV3/nvkPp/QkZoq9CZjk7sm2TF33VmBm\n+xBm2O0MfAoc7+6zzGwIut45l0jWBgD7NfK2PuOtYwrwI0It0E7AL4EXEp/9nF3zYkskRErdWKAM\nOCDuQNqAmUB/worDJwL3mtnQeEMqTWa2MyFB/qa7r487nrbC3VPX0XjLzF4F3gdOJnz+c6KomjaI\ntmiY5M5CQp8UXf9WYGa3AscAh7j7xylv6bq3Anff4O7vuvvr7n4lofPfReh6t4YKYAdgmpmtN7P1\nwMHARWa2jnAXrGveyty9DphNWBAzZ5/zokokEplsctEwYLNFw3Ky+Ig0zd3nET5gqdd/G8JoA13/\nLCSSiOOAQ919fup7uu550w7opOvdKp4mLMQ4gFAL1B94DbgP6O/u76Jr3urMrCshiViQy895MTZt\njAbGJ1YhfRWoImXRMMmOmW1F+KBZYtNuZtYfWO7uHxCqJ68ys3cIS7j/ijBq5pEYwi0JZjYWqASO\nBVaZWfIOoc7d1yb+reueQ2Z2A/AEMB/YGjiVcId8RGIXXe8ccvdVQMP5C1YBy9x9RmKTrnmOmdnv\ngMcIzRlfBq4F1gN/SeySk2tedImEt7xomGRnP8KwLE88RiW23wOc6e6/NbMuwDjC6IIXgaPdfV0c\nwZaIcwnX+vkG288A7gXQdc+5HQmf6Z2AOuAN4IjkaAJd77zYbJ4aXfNWsTPwALA9sASYBAx292WQ\nu2uuRbtEREQksqLqIyEiIiKFRYmEiIiIRKZEQkRERCJTIiEiIiKRKZEQERGRyJRIiIiISGRKJERE\nRCQyJRIiIiISmRIJERERiUyJhIiIiESmREJEREQi+/8PJS8s9PxV2QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110fa8310>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pl.plot(gtau)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}