{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Wang Landau algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import *\n",
    "from numpy import random\n",
    "from numba import jit\n",
    "\n",
    "@jit(nopython=True)\n",
    "def CEnergy(latt):\n",
    "    \"Energy of configuration for the 2D Ising model\"\n",
    "    N = shape(latt)[0]\n",
    "    Ene = 0\n",
    "    for i in range(len(latt)):\n",
    "        for j in range(len(latt)):\n",
    "            S = latt[i,j]   # Spin, can be either +1 or -1\n",
    "            WF = latt[(i+1)%N,j]+latt[i,(j+1)%N]+latt[(i-1)%N,j]+latt[i,(j-1)%N]\n",
    "            Ene += -S * WF\n",
    "    return Ene/2.\n",
    "\n",
    "def RandomL(N):\n",
    "    \"Random lattice corresponding to infinite temperature\"\n",
    "    return array(sign(2*random.random((N,N))-1),dtype=int) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def PrepareEnergies(N):\n",
    "    Energies = (array(4*arange(-int(N*N/2),int(N*N/2)+1),dtype=int)).tolist() # -2 N^2...2N^2 in steps of 4\n",
    "    Energies.pop(1)  # take out E[1]\n",
    "    Energies.pop(-2) # take out E[-2]\n",
    "    Energies = array(Energies) # make array of energies again\n",
    "    Emin, Emax = Energies[0],Energies[-1]\n",
    "    #index array to energies\n",
    "    indE = -ones(Emax+1-Emin, dtype=int) # index table to get index to particular energy g(E)~g[indE[E]]\n",
    "    for i,E in enumerate(Energies):\n",
    "        indE[E-Emin]=i\n",
    "    # Ising lattice at infinite T\n",
    "    return (Energies, indE, Emin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Emin= -200\n",
      "[-200 -192 -188 -184 -180 -176 -172 -168 -164 -160 -156 -152 -148 -144\n",
      " -140 -136 -132 -128 -124 -120 -116 -112 -108 -104 -100  -96  -92  -88\n",
      "  -84  -80  -76  -72  -68  -64  -60  -56  -52  -48  -44  -40  -36  -32\n",
      "  -28  -24  -20  -16  -12   -8   -4    0    4    8   12   16   20   24\n",
      "   28   32   36   40   44   48   52   56   60   64   68   72   76   80\n",
      "   84   88   92   96  100  104  108  112  116  120  124  128  132  136\n",
      "  140  144  148  152  156  160  164  168  172  176  180  184  188  192\n",
      "  200]\n",
      "[ 0 -1 -1 -1 -1 -1 -1 -1  1 -1 -1 -1  2 -1 -1 -1  3 -1 -1 -1  4 -1 -1 -1\n",
      "  5 -1 -1 -1  6 -1 -1 -1  7 -1 -1 -1  8 -1 -1 -1  9 -1 -1 -1 10 -1 -1 -1\n",
      " 11 -1 -1 -1 12 -1 -1 -1 13 -1 -1 -1 14 -1 -1 -1 15 -1 -1 -1 16 -1 -1 -1\n",
      " 17 -1 -1 -1 18 -1 -1 -1 19 -1 -1 -1 20 -1 -1 -1 21 -1 -1 -1 22 -1 -1 -1\n",
      " 23 -1 -1 -1 24 -1 -1 -1 25 -1 -1 -1 26 -1 -1 -1 27 -1 -1 -1 28 -1 -1 -1\n",
      " 29 -1 -1 -1 30 -1 -1 -1 31 -1 -1 -1 32 -1 -1 -1 33 -1 -1 -1 34 -1 -1 -1\n",
      " 35 -1 -1 -1 36 -1 -1 -1 37 -1 -1 -1 38 -1 -1 -1 39 -1 -1 -1 40 -1 -1 -1\n",
      " 41 -1 -1 -1 42 -1 -1 -1 43 -1 -1 -1 44 -1 -1 -1 45 -1 -1 -1 46 -1 -1 -1\n",
      " 47 -1 -1 -1 48 -1 -1 -1 49 -1 -1 -1 50 -1 -1 -1 51 -1 -1 -1 52 -1 -1 -1\n",
      " 53 -1 -1 -1 54 -1 -1 -1 55 -1 -1 -1 56 -1 -1 -1 57 -1 -1 -1 58 -1 -1 -1\n",
      " 59 -1 -1 -1 60 -1 -1 -1 61 -1 -1 -1 62 -1 -1 -1 63 -1 -1 -1 64 -1 -1 -1\n",
      " 65 -1 -1 -1 66 -1 -1 -1 67 -1 -1 -1 68 -1 -1 -1 69 -1 -1 -1 70 -1 -1 -1\n",
      " 71 -1 -1 -1 72 -1 -1 -1 73 -1 -1 -1 74 -1 -1 -1 75 -1 -1 -1 76 -1 -1 -1\n",
      " 77 -1 -1 -1 78 -1 -1 -1 79 -1 -1 -1 80 -1 -1 -1 81 -1 -1 -1 82 -1 -1 -1\n",
      " 83 -1 -1 -1 84 -1 -1 -1 85 -1 -1 -1 86 -1 -1 -1 87 -1 -1 -1 88 -1 -1 -1\n",
      " 89 -1 -1 -1 90 -1 -1 -1 91 -1 -1 -1 92 -1 -1 -1 93 -1 -1 -1 94 -1 -1 -1\n",
      " 95 -1 -1 -1 96 -1 -1 -1 97 -1 -1 -1 -1 -1 -1 -1 98]\n"
     ]
    }
   ],
   "source": [
    "Ene,indE,Emin = PrepareEnergies(10)\n",
    "print('Emin=',Emin)\n",
    "print(Ene)\n",
    "print(indE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import random\n",
    "def WangLandau(Nitt, N, flatness):\n",
    "    \"Wang Landau in Python\"\n",
    "    (Energies, indE, Emin) = PrepareEnergies(N)\n",
    "    latt = RandomL(N)\n",
    "    return RunWangLandau(Nitt,Energies,latt,indE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "@jit(nopython=True)\n",
    "def RunWangLandau(Nitt,Energies,latt,indE):\n",
    "    N   = len(latt)\n",
    "    Ene = int(CEnergy(latt))\n",
    "    #min,maximum energy\n",
    "    Emin, Emax = Energies[0],Energies[-1]\n",
    "    # Logarithm of the density of states\n",
    "    lngE = zeros(len(Energies))\n",
    "    # Histogram\n",
    "    Hist = zeros(len(Energies))\n",
    "    # modification factor\n",
    "    lnf = 1.0   # f = exp(lnf)=e\n",
    "    N2 = N*N\n",
    "    for itt in range(Nitt):\n",
    "        t = int(random.rand()*N2)\n",
    "        (i, j) = (int(t/N), t%N)\n",
    "        S = latt[i,j]\n",
    "        WF = latt[(i+1)%N,j]+latt[i,(j+1)%N]+latt[(i-1)%N,j]+latt[i,(j-1)%N]\n",
    "        Enew = Ene + int(2*S*WF) # the energy change if we flip the spin\n",
    "        # P = g(E)/g(Enew) = exp(log(g(E))-log(g(Enew)))\n",
    "        # P = exp(lngE(E)-lngE(Enew))\n",
    "        lgnew = lngE[indE[Enew-Emin]]\n",
    "        lgold = lngE[indE[Ene-Emin]]\n",
    "        P = 1.0\n",
    "        if lgold-lgnew < 0 : P=exp(lgold-lgnew) # P = g_old/g_new = exp(log(g_old)-log(g_new))\n",
    "        if P > random.rand():\n",
    "            # accept the step\n",
    "            latt[i,j] = -S\n",
    "            Ene = Enew\n",
    "        Hist[indE[Ene-Emin]] += 1\n",
    "        lngE[indE[Ene-Emin]] += lnf\n",
    "        \n",
    "        if (itt+1) % 1000 == 0: # checking for flatness of the histogram\n",
    "            aH = sum(Hist)/N2 # average\n",
    "            mH = min(Hist)\n",
    "            if mH > aH*flatness:  # histogram is flat\n",
    "                Hist[:] = 0\n",
    "                lnf /= 2.\n",
    "                print(itt, 'histogram is flat', mH, aH, 'f=', exp(lnf))\n",
    "    return (lngE, Hist, Energies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10649999 histogram is flat 24269.0 26625.0 f= 1.6487212707001282\n",
      "13240999 histogram is flat 5840.0 6477.5 f= 1.2840254166877414\n",
      "15756999 histogram is flat 5673.0 6290.0 f= 1.1331484530668263\n",
      "18197999 histogram is flat 5495.0 6102.5 f= 1.0644944589178593\n",
      "20460999 histogram is flat 5104.0 5657.5 f= 1.0317434074991028\n",
      "22623999 histogram is flat 4874.0 5407.5 f= 1.0157477085866857\n",
      "25381999 histogram is flat 6222.0 6895.0 f= 1.007843097206448\n",
      "28732999 histogram is flat 7571.0 8377.5 f= 1.0039138893383475\n",
      "33813999 histogram is flat 11457.0 12702.5 f= 1.0019550335910028\n",
      "41300999 histogram is flat 16862.0 18717.5 f= 1.0009770394924165\n",
      "49078999 histogram is flat 17505.0 19445.0 f= 1.0004884004786945\n",
      "60349999 histogram is flat 25418.0 28177.5 f= 1.0002441704297478\n",
      "70789999 histogram is flat 23491.0 26100.0 f= 1.0001220777633837\n",
      "92226999 histogram is flat 48270.0 53592.5 f= 1.0000610370189331\n"
     ]
    }
   ],
   "source": [
    "from numpy import random\n",
    "\n",
    "flatness = 0.9\n",
    "N = 20\n",
    "Nitt = int(100e6)\n",
    "\n",
    "#N=10\n",
    "#Nitt = int(10e6)\n",
    "\n",
    "(lngE, Hist, Energies) = WangLandau(Nitt, N, flatness)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Proper normalization of the Density of states\n",
    "# g *= 4/(g[0]+g[-1])\n",
    "# log(g) += log(4)-log(g[-1])-log(1+exp(log(g[0])-log(g[-1])))\n",
    "# \n",
    "if lngE[-1]>lngE[0]:\n",
    "    lgC = log(4)-lngE[-1]-log(1+exp(lngE[0]-lngE[-1]))\n",
    "else:\n",
    "    lgC = log(4)-lngE[0]-log(1+exp(lngE[-1]-lngE[0]))\n",
    "    \n",
    "lngE += lgC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.0"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "exp(lngE[0])+exp(lngE[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fda096c9fd0>]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAArIklEQVR4nO3deViVdf7/8eebXREQFBAVBRVF3JVMsyy31Ky0mWmyaXGmGmvGapqtbJ1qWqa+bTPNtI7N2KJOm2XpmGlmi5XiDiKKioIiqCAgO5zP7w9Pv6FCReCczzmH9+O6uDjn5hzul/cFL2/u5fMRYwxKKaV8i5/tAEoppVqflrtSSvkgLXellPJBWu5KKeWDtNyVUsoHBdgOANC5c2eTkJBgO4ZSSnmVDRs2HDHGRDf2NY8o94SEBNLS0mzHUEopryIi+072NT0so5RSPkjLXSmlfJCWu1JK+SAtd6WU8kFa7kop5YO03JVSygdpuSullA/yiOvclbKtqraefUcrKCitorCsmuLyGmodDowBYwztggIICw6gc1gQXcLbERcRQsf2gYiI7ehKNUrLXbU5dfUOMg6W8vWeo2zaf4ydBWXkHC3HcYZTG3QKDWJAtwgGdQtnRM9IRiZ2okOw/kopz6A/iapNKKuq5ZMdhSxPP8Tnu45wvLoOgIRO7ekfF84lQ7rSO6YDcREhRHcIJqpDEEH+fny7Y15ZU09ZVR1HjldzqKSKA8cq2VlQxrYDpby4Zg91DkOAnzA0viPjkmO4eHAcPTuFWvwXq7ZOPGEmptTUVKPDD6jWVlfvYNWOQt5cn8vnu45QU+8gOiyYif1jOad3J87uFUVMWEiL11NVW8+GfcV8mX2EL7KPsDWvBIDB3SO4ZHBXpg/r2irrUer7RGSDMSa10a9puStfU1haxaL1uSxct5/8kiq6hIcwbXAcUwd2YXiPSPz8XHucPK+4gmXb8vlgSz7bDpQQ5O/HJUO6ct25CQzoGuHSdau2RctdtQk5R8p5/tPdvLspj9p6w3lJnblmVE/GJ8cQ4G/nwrDdh4/z6toc3tqQR0VNPWcnRnHL+CTG9OmkJ2NVi2m5K5+WX1LJUyt28s7GPAL8/bgiNZ7rzk0ksbPnHPMuqazlzfW5vPLlXvJLqjg7MYo/TO7HWQlRtqMpL6blrnxSWVUtL6zZzbwv9uJwwLWjezL7/F4efXy7uq6eRety+fvqbA6XVTO2bzT3TOtP39gw29GUF9JyVz7FGMO7Gw/wyLJMjpbXMH1oV/5wYT/io9rbjtZklTX1vPZ1Dv9YvZvj1XVcO7onv53Ul/CQQNvRlBfRclc+Y9/Rcu5avI0vs48yvEdH7r90AIO7d7Qdq9mKymt4YkUWC9ftJyYsmIdmDGJSSqztWMpLaLkrr1fvMLz8+R6e/ngnQf5+3D41matG9nD5lS/usiX3GHe8s5Udh8qYNjiOBy4dQOcOwbZjKQ93qnLXm5iUxztwrJLfLtrMupwiJg+I5YFLB9IlwnOPqzfHkPiOLLn5XF5cs5tnP8nmmz1FPHH5YC7oF2M7mvJSOnCY8mhLt+Yz9ZnP2J5fytNXDOHFa1J9rti/FRTgxy0Tklhyyxg6hQbx83+t58EPtlNVW287mvJCWu7KI1XV1jP3na3MWbCRXtEdWHrruVw2rLvtWG6R3CWc928ew6zRPXnly7385IW15BZV2I6lvIyWu/I4+SWVXPHiVyxan8uvL+jNWzeNbnPjtIQE+vPA9IH889pU9h+tYNrfPmfl9gLbsZQX0XJXHmV9ThGXPPsl2YXHefGaEdw+JZlAS3eXeoKJKbF8eMt59OjUnhteTePJFVk4znT4StUmnfa3RkTiRWS1iGSKSIaI/Ma5/H4ROSAim50fFzV4z50iki0iWSIy2ZX/AOU7Xv96H1e+9DVhIQG8N2cMkwd0sR3JI/To1J63bzqHn6Z259lPsrll4SY9Dq9OqylXy9QBvzfGbBSRMGCDiHzs/NrTxpgnGr5YRFKAmcAAoCuwUkT6GmP0p1E1qt5heOCDDF79ah8X9IvmrzOHEdFOb+ZpKCTQn8d+PJg+MR149L87yDtWycvXjvDou3GVXafdczfG5BtjNjoflwGZQLdTvGU6sMgYU22M2QtkAyNbI6zyPZU19dz0+gZe/Wofs8f2Yt6ss7TYT0JEmD22Ny9cPYKdh8qY8fcvyTpUZjuW8lBndDBTRBKAYcA3zkU3i8hWEXlFRCKdy7oBuQ3elkcj/xmIyGwRSRORtMOHD595cuX1jhyvZubLX7Mys4AHpw/grov64+8jNyW50uQBXXjrptHUG8PlL6xlw75i25GUB2pyuYtIB+Ad4DZjTCnwPNAbGArkA09++9JG3v6DM0DGmJeMManGmNTo6Ogzza283N4j5fzoubXsyC/lhatHcO3oBNuRvMrAbhG8fdM5RIUGcfU/v2HNTt1BUt/VpHIXkUBOFPsbxph3AYwxBcaYemOMA3iZ/x16yQPiG7y9O3Cw9SIrb5d+oIQfP7+W49V1LJw9Sk+cNlN8VHveuukcEjuHcsP89XywRX/N1P805WoZAeYBmcaYpxosj2vwssuAdOfjJcBMEQkWkUQgCVjXepGVN9uwr5grX/6adoH+vPOrcxjeI/L0b1InFR0WzMLZoxga35FbF23ijW/22Y6kPERTrpYZA1wDbBORzc5ldwFXishQThxyyQFuBDDGZIjIm8B2TlxpM0evlFEA3+w5yi/+vZ6YsGDe+OUounVsZzuST4hoF8ir153NnAUbuXtxOrV1Dn4+JtF2LGWZjgqp3GJ9ThGzXllH147tWHDD2cSE6yV8ra2mzsGcBRv5ePuJE9R6HsP3nWpUyLZ7659ym437i/n5K+voEh6ixe5CQQF+/ONnw5nYP5b73s/gta9ybEdSFmm5K5faknuMWfPW0TksmAW/HKXF7mJBAX48d9VwJvaP4d73M3j9az0G31ZpuSuXST9QwjXzvqFjaCALfznKZ4fq9TRBAX7846rhTEiO4Z730lm4br/tSMoCLXflErsPH+eaed8QFhLIghtG0VVPnrpVcIA/z109nAv6RXPX4m0s25ZvO5JyMy131eoKS6uY9co6/P2EBb8826smrvYlwQH+PH/VCEb0iOS2RZv5MvuI7UjKjbTcVasqq6rl5/9aT1F5Da/8/Kw2Nw67p2kX5M+8WWfRKzqU2a+msTXvmO1Iyk203FWrqalz8KvXN5JVUMY/rhrO4O4dbUdSQET7QOZfN5JI59R92YXHbUdSbqDlrlqFMYY73tnKF9lHePRHgxinEzt7lNjwEF67/mz8BGa9so78kkrbkZSLabmrVvH4R1ks3nSA30/qy09T40//BuV2iZ1D+fcvRlJSWcu189ZRXF5jO5JyIS131WLz1+bw/Ke7+dnZPbh5fB/bcdQpDOwWwcvXprKvqILr5q+noqbOdiTlIlruqkWWp+dz/wcZTOwfy4OXDuDEOHPKk43u3YlnrxzGltxj3LJgE/U6J6tP0nJXzbY+p4hbF21maHxHnr1yGAFteCJrbzN5QBfuv3QAq3YU8vjyHbbjKBdoyqiQSv1AXnEFN762gW4d2zFv1lm0C/K3HUmdoWtHJ7Cr4DgvfraH3jEd9FyJj9FdLXXGKmrqmP3qBmrrHfxzVipRoUG2I6lmuu+SFM7t05m7F29j3d4i23FUK9JyV2fEGMMf395K5qFS/nblMHpHd7AdSbVAoP+JkSTjI9tz42tp7D9aYTuSaiVa7uqMPPfpbpZuzeeOKcl6LbuPiGgfyLyfn4XDwPXz13O8Wq+g8QVa7qrJVm4v4IkVWUwf2pUbx/ayHUe1osTOoTx31XD2HCnnD29uwRMm8VEto+WumiS7sIzb/rOZAV3DeezHg/WSRx80pk9n7pyazPKMQzz36W7bcVQLabmr0yqpqOWXr24gJNCPl65JJSRQr4zxVdefm8glQ7ryxIosPtt52HYc1QJa7uqU6h2GWxZtIq+4gheuHqHjsvs4EeGxHw+ib0wYtyzcRG6RnmD1Vlru6pS+3YN7cPpAUhOibMdRbtA+KIAXrxmBMYYbX9tAZU297UiqGbTc1Umtyizg+U93c+XIHlw5softOMqNEjqH8teZw8g8VMpdi7fpCVYvpOWuGpVXXMHv3txCSlw4f7okxXYcZcG45Bh+O7EvizcdYP7aHNtx1BnSclc/UFPnYM6CTTgchuevHq4nUNuwm8f1YWL/WB5amklajt7B6k203NUPPLIsky25x/i/ywfrNHltnJ+f8NQVQ+jasR23LtzEsQodA95baLmr71i2LZ9/r83hujGJTBkYZzuO8gDhIYE8e+UwDh+v5va3t+rxdy+h5a7+v31Hy7n97a0M69GRuVOTbcdRHmRIfEfumJLMiu0FvPrVPttxVBOcttxFJF5EVotIpohkiMhvnMujRORjEdnl/BzZ4D13iki2iGSJyGRX/gNU66ipc3DLwk34Cfz9Z8MJCtD/99V3XX9uIuOTY3h4aSbpB0psx1Gn0ZTf4Drg98aY/sAoYI6IpABzgVXGmCRglfM5zq/NBAYAU4DnRETPyHm4//toB1vzSnj8J0PopjcqqUaICE9cPoTI0EBuWbhJBxjzcKctd2NMvjFmo/NxGZAJdAOmA/OdL5sPzHA+ng4sMsZUG2P2AtnAyFbOrVrR6h2FvPz5Xq4d3ZMpA7vYjqM8WFRoEH+dOYx9R8u5771023HUKZzR394ikgAMA74BYo0x+XDiPwDg2/FfuwG5Dd6W51z2/e81W0TSRCTt8GEdw8KWgtIqfv/WFpK7hHHXRf1tx1FeYFSvTtw6IYl3Nx3g7Q15tuOok2hyuYtIB+Ad4DZjTOmpXtrIsh+cXjfGvGSMSTXGpEZHRzc1hmpF9Q7DbYs2U1lTz99/NkyvZ1dNdsv4JEb1iuLe99LZffi47TiqEU0qdxEJ5ESxv2GMede5uEBE4pxfjwMKncvzgIaTMXYHDrZOXNWanludzVd7jvLA9AH0iQmzHUd5EX8/4a8zh9EuyJ85b2ykqlbHn/E0TblaRoB5QKYx5qkGX1oCzHI+ngW832D5TBEJFpFEIAlY13qRVWtIyyni6ZU7mT60K5eP6G47jvJCseEhPHn5EHYcKuPhpZm246jvacqe+xjgGmC8iGx2flwE/AWYJCK7gEnO5xhjMoA3ge3AcmCOMUb/W/cgx6vr+O2bm+kW2Y6HZgzUiTdUs41LjuGX5yXy2tf7WJ6ebzuOaiDgdC8wxnxB48fRASac5D0PAw+3IJdyoYeXZpJXXMmbN44mLCTQdhzl5f44OZl1e4uY++42hvWIJDY8xHYkhd6h2uas3lHIwnX7mT22F2fp+OyqFQQF+PHMzGFU1zr4w1tbcDh0eAJPoOXehhSX13D7O1vpFxvG7yb1tR1H+ZDEzqHcc3F/Pt91hPlf5diOo9BybzOMMdzzXjrHKmp4+oqhBAfoZY+qdf1sZA8mJMfwl//uYFdBme04bZ6WexuxZMtBlm7L57eT+pLSNdx2HOWDRIS//HgwHYID+M2izdTUOWxHatO03NuA/JJK7n0vnRE9I7lxbG/bcZQPiw4L5rEfD2Z7filPr9xpO06bpuXu44wx3P72VmrrDU9ePgR/P73sUbnWxJRYrhwZzwtrdvPNnqO247RZWu4+7vWv9/H5riPcPa0/CZ11ViXlHvdMS6FnVHt+9+YWyqpqbcdpk7TcfdjeI+U8vCyT8/tGc9XZPWzHUW1IaHAAT/50KPkllTyybIftOG2SlruPqqt38Ls3NxMc4M/jPxmsd6EqtxvRM5Jfju3FwnX7WbNTR351Ny13H/XiZ3vYtP8YD80YqHcMKmt+O7EvSTEduOPtrZRU6uEZd9Jy90FZh8p4ZuVOLh4cxyVDutqOo9qwkEB/nrh8CIePV/PQh9ttx2lTtNx9TF29g9vf3kJ4SCAPTh9oO45SDInvyK/O781bG/L4ZEeB7Ththpa7j5n3xV625JXw4PSBRIUG2Y6jFAC3TOhDcpcw5r6zjWMVNbbjtAla7j5k9+HjPPnxTiYPiOWiQToXqvIcwQEnDs8UldfwwAd6eMYdtNx9hMNhuOPtrbQL9OfP03WMduV5BnaLYM64PizedICPMg7ZjuPztNx9xKtf5ZC2r5j7Lk4hRq+OUR7q5vF9SIkL5+7F2ygq18MzrqTl7gNyiyp4bHkWF/SL5kfDu9mOo9RJBfr78eRPh1BSWct976fbjuPTtNy9nDGGO97Zir+f8Mhlg/RwjPJ4/ePCuXV8Eh9uzWdVpl494ypa7l5u0fpc1u4+yl0X9adrx3a24yjVJDee35u+sR249710yqvrbMfxSVruXiy/pJKHl2ZyTu9OXDky3nYcpZosKMCPR380mPzSKp5coUMDu4KWu5cyxnDP4nTqHYa//EjHjlHeZ0TPSK4+uyf/XruXLbnHbMfxOVruXmrptnxW7Sjk9xf2pUen9rbjKNUsf5zSj+iwYOa+u43aep25qTVpuXuhkopa7l+yncHdI/jFmETbcZRqtvCQQB64dCCZ+aW88sVe23F8ipa7F3r0v5kUV9Tw6I8G6cxKyutNGdiFC1NieXrlTvYfrbAdx2douXuZr/ccZdH6XG44L5EBXSNsx1GqVTwwfQABfn7c/d42jDG24/gELXcvUlVbz12LtxEf1Y7bJvS1HUepVhMX0Y4/Tu7H57uO8P7mg7bj+AQtdy/y3Ops9hwu5+EZg2gX5G87jlKt6upRPRka35EHP9xOsQ5N0GKnLXcReUVECkUkvcGy+0XkgIhsdn5c1OBrd4pItohkichkVwVva3YWlPH8mt1cNqwbY/tG246jVKvz9xMe/dEgSitreXhZpu04Xq8pe+7/BqY0svxpY8xQ58cyABFJAWYCA5zveU5EdBezhRwOw53vbqNDcAD3TOtvO45SLtM/LpzZY3vx9oY81mYfsR3Hq5223I0xnwFFTfx+04FFxphqY8xeIBsY2YJ8Cliwbj8b9hVzz7QUOnUIth1HKZe6dUISPTu1567F26iqrbcdx2u15Jj7zSKy1XnYJtK5rBuQ2+A1ec5lPyAis0UkTUTSDh/WmdFP5lBJFY/9dwdj+nTSER9VmxAS6M8jlw0i52gFz36yy3Ycr9Xccn8e6A0MBfKBJ53LG7voutHrmowxLxljUo0xqdHRegz5ZO5fkkFNvYOHZ+iIj6rtGNOnMz8e3p0X1+xhx6FS23G8UrPK3RhTYIypN8Y4gJf536GXPKDhCFbdAb2uqZk+yjjE8oxD/GZiEgmdQ23HUcqt7p7Wn/B2gdz57jYcDr32/Uw1q9xFJK7B08uAb6+kWQLMFJFgEUkEkoB1LYvYNpVV1fKn9zNI7hLGL8/rZTuOUm4XFRrEvRf3Z9P+Y7zxzT7bcbxOwOleICILgQuAziKSB/wJuEBEhnLikEsOcCOAMSZDRN4EtgN1wBxjjJ4RaYZnVu6ioKyK568eTqC/3o6g2qYZQ7vx7sYDPLY8i0kpXegSoVNINpV4wq2+qampJi0tzXYMj7H9YCmX/P0LrjgrnkcuG2Q7jlJW7T9awYXPrOH8vtG8eE2q7TgeRUQ2GGMa3Si6S+hhHA7Dve+nE9EukNsn97MdRynrenRqz20T+/JRRgHL0w/ZjuM1tNw9zNsb8tiwr5i5U5Pp2D7IdhylPML15ybSPy6cPy1J57hOy9ckWu4epLi8hkf/m0lqz0h+Mry77ThKeYxAfz8evmwgBaXV/HWlTsvXFFruHuTxj7IorarjzzMG4qfjtCv1HcN7RDLzrHhe+TKHrENltuN4PC13D7FxfzGL1u/nF+ck0D8u3HYcpTzS7VOSCQsJ4N7303Xc99PQcvcA9Q7Dve+lExMWzG2TdJx2pU4mKjSIO6Yks25vEYs3HbAdx6NpuXuA17/eR8bBUu69OIUOwae99UCpNu2K1HiGxnfkkWWZlFTW2o7jsbTcLSssq+KJj7I4L6kz0wbFnf4NSrVxfn7CQzMGUlRew1MrsmzH8Vha7pY9sjST6joHD1w6QAcGU6qJBnaL4OpRPXnt632kHyixHccjablbtHb3Ed7bfJAbz+9Fr+gOtuMo5VV+f2E/okKDuOe9dB1YrBFa7pbU1Dm47/0M4qPaMWdcH9txlPI6Ee0Cueui/mzOPcZ/0nJP/4Y2Rsvdknlf7CW78Dj3XzKAkECdiVCp5rhsWDdGJkbx2PIdFOmk2t+h5W5BXnEFf1u1iwtTYpnQP9Z2HKW8lojw5+kDKauq4/HlO2zH8Sha7hY8+MF2DIb7LkmxHUUpr9evSxjXjUlg0fpcNu4vth3HY2i5u9knOwpYsb2AWyck0T2yve04SvmE30zsS2x4MPe+l069nlwFtNzdqqq2nj8tyaB3dCg3nKuzKynVWjoEB3DvxSlkHCzl9a911ibQcner51Znk1tUyZ9nDCQoQDe9Uq1p2qA4zu3TmSdWZFFYVmU7jnXaMG6y90g5L6zZw/ShXTmnd2fbcZTyOSLCA9MHUFVbz6PL9OSqlrsbGGO47/10ggP8uPui/rbjKOWzekd3YPbYXizedICv9xy1HccqLXc3WLbtEJ/vOsLvL+xLTLhO8KuUK908LoluHdtx3/vp1NY7bMexRsvdxY5X1/HghxkM6BrO1aN62o6jlM9rF+TP/ZcOYGfBcf715V7bcazRcnexZz7eSUFpNX+eMZAAf93cSrnDpJRYJiTH8MzKXeSXVNqOY4W2jQvtOFTKv9bmcOXIeIb3iLQdR6k25f5LB1DvMDz0YabtKFZoubuIw2G4Z3E64SEB3D452XYcpdqc+Kj2zBnXh6Xb8vls52HbcdxOy91F3tmYR9q+Yu6c2p/I0CDbcZRqk2aP7UVCp/b8aUkG1XX1tuO4lZa7CxyrqOHR/+5geI+O/GREd9txlGqzQgL9eWD6QPYeKeelNXtsx3ErLXcX+L+PsjhWUcNDMwbh56ezKyll0/l9o5k6sAt/X51NblGF7Thuc9pyF5FXRKRQRNIbLIsSkY9FZJfzc2SDr90pItkikiUik10V3FOlHyhhwbr9XDs6gZSu4bbjKKWAey9Owd9PeOCDDNtR3KYpe+7/BqZ8b9lcYJUxJglY5XyOiKQAM4EBzvc8JyJtZiYKYwz3L8kgsn0Qv53U13YcpZRT147tuHVCEiszC1m5vcB2HLc4bbkbYz4Dir63eDow3/l4PjCjwfJFxphqY8xeIBsY2TpRPd+SLQdJ21fM7ZP7EdEu0HYcpVQD141JpE9MB+7/IIPKGt8/udrcY+6xxph8AOfnGOfybkDDyQzznMt+QERmi0iaiKQdPuz9lymVV9fxyLJMBnWL4PLUeNtxlFLfExTgx5+nDySvuJLn1+y2HcflWvuEamNnDxsdOd8Y85IxJtUYkxodHd3KMdzvH6uzKSit5v5LTxzbU0p5ntG9O3HJkK68uGY3ecW+fXK1ueVeICJxAM7Phc7leUDD3dbuwMHmx/MOOUfK+efne/nRsG6M6BllO45S6hTmTk1GBJ8fFri55b4EmOV8PAt4v8HymSISLCKJQBKwrmURPd9DS7cT6C/Mnap3oirl6bp1bMdN5/dm6bZ8nx4WuCmXQi4EvgL6iUieiFwP/AWYJCK7gEnO5xhjMoA3ge3AcmCOMcanz1yszipkZWYht05I0uF8lfISN47tTdeIEB74YLvPzrkqxtj/h6Wmppq0tDTbMc5YTZ2DKc98BsDy28bq1HlKeZEPtx7k5gWbeOSyQfzs7B624zSLiGwwxqQ29jVtoxb499q97DlSzr2XpGixK+Vlpg2KY2RCFE+syKKkstZ2nFanjdRMhaVV/HXlLiYkxzCuX8zp36CU8igiwn2XpFBcUcPfVu2yHafVabk302PLs6itN9x7cYrtKEqpZhrYLYKZZ8Uzf20O2YVltuO0Ki33Zti4v5h3NuZx/XmJJHQOtR1HKdUCf7iwH+2C/Hnww0w84Rxka9FyP0MOx4nxY2LDg7l5XB/bcZRSLdSpQzC/mZDEZzsPszqr8PRv8BJa7mfo7Q15bM0r4c6p/QkNDrAdRynVCq4dnUCv6FD+/GEmNXUO23FahZb7GSiprOWx5TsY0TOS6UO72o6jlGolQQF+3HdxCnuPlDN/bY7tOK1Cy/0M/G3VLooqanjg0gGI6PgxSvmSC/rFMD45hr+t2sXhsmrbcVpMy72J9h4p59WvcrgiNZ6B3SJsx1FKucA90/pTWVvPEx9l2Y7SYlruTfSX/2YS6O/H7y7USTiU8lW9ojvwizEJvLkhl+0HS23HaREt9yb4Zs9RPsoo4Ffn9yYmTMePUcqX3TwuiYh2gTyyzLsvjdRyPw2Hw/DQ0kziIkK44bxetuMopVwson0gt45P4ovsI3y603snEtJyP433txxg24ES/jj5xI0OSinfd/WoniR0as8jSzOpq/fOSyO13E+hsqaex5dnMahbBDOGNjpboFLKBwUF+DF3ajK7Co/z1oY823GaRcv9FOZ9sYf8kirumdYfP506T6k2ZfKALqT2jOTJFTs5Xl1nO84Z03I/icKyKp77dDeTB8Rydq9OtuMopdxMRLh7Wn+OHK/mJS+cUFvL/SSe/ngnNXUO5k7tbzuKUsqSYT0iuXhwHC99vodDJVW245wRLfdG7Cwo4z/rc7lmdE8SddRHpdq0O6Yk43DAEyu868YmLfdGPL58B6FBAdw6Psl2FKWUZfFR7fn5mATe2ZhHxsES23GaTMv9e9btLWJlZiE3XdCbyNAg23GUUh5gzgV9vO7GJi33BowxPPrfTGLDg7luTKLtOEopDxHRPpDfTEjiy+yjXnNjk5Z7Ax9lFLBp/zF+O7Gv3rCklPqOq84+cWPTo8syqXd4/t67lrtTXb2Dxz/aQe/oUH4yorvtOEopDxMU4MfvL+zHzoLjLNlywHac09Jyd3prQx57Dpdz+5RkAvx1syilfmjaoDhS4sJ5ynmptCfTFgMqaup4+uOdjOgZyYUpsbbjKKU8lJ+f8MfJ/cgtquQ/abm245ySljvwry9zKCyrZu7UZJ1hSSl1Shf0i+ashEieXbWLypp623FOqs2Xe1F5DS98upuJ/WM5KyHKdhyllIcTEf44OZnCsmrmf5VjO85JtajcRSRHRLaJyGYRSXMuixKRj0Vkl/NzZOtEdY1/rM6mvKaOO6b0sx1FKeUlRiZGcUG/aJ7/dDcllbW24zSqNfbcxxljhhpjUp3P5wKrjDFJwCrnc4+UW1TBa1/t4/IR8STFhtmOo5TyIn+4sB8llbX88/M9tqM0yhWHZaYD852P5wMzXLCOVvHUxzsRgdsm6TADSqkzM7BbBNMGxzHvi70cLqu2HecHWlruBlghIhtEZLZzWawxJh/A+TmmsTeKyGwRSRORtMOH3X/HV9ahMt7bfIBfjEkkLqKd29evlPJ+v5vUl+o6B899mm07yg+0tNzHGGOGA1OBOSIytqlvNMa8ZIxJNcakRkdHtzDGmXtyRRYdggK46XydF1Up1Ty9ozvwk+HdeePr/Rw8Vmk7zne0qNyNMQednwuBxcBIoEBE4gCcnwtbGrK1bc07xortBdxwXi86ttfBwZRSzXfLhD44jOFFD5vQo9nlLiKhIhL27WPgQiAdWALMcr5sFvB+S0O2tidW7CSyfSDXnZtgO4pSyst1j2zPj4d3Z+H6XApLPWdCj5bsuccCX4jIFmAdsNQYsxz4CzBJRHYBk5zPPca6vUV8tvMwv7qgN2EhgbbjKKV8wK/H9abeYXjxM8+5ciaguW80xuwBhjSy/CgwoSWhXMUYwxMrsogJC+aaUQm24yilfETPTqFMH9KVN77Zx68u6E3nDsG2I7WtO1S/2nOUdXuLmDOujw7pq5RqVb8e14fqOgf//Hyv7ShAGyv3Z1dlExMWzBVnxduOopTyMX1iOjBtUByvfZVDcXmN7Thtp9zTcor4as9RZo/tRUig7rUrpVrfzeP7UF5Tz7++tL/33mbK/W+fZNMpNIirzu5pO4pSykcldwln8oBY/rU2h9Iqu2POtIly35x7jM92HuaG83rpsXallEvdMj6Jsqo65n+ZYzVHmyj3v3+yi47tA7lmtO61K6Vca2C3CMYnxzDvy70cr66zlsPny33HoVJWZhbyi3MS6RDc7Cs/lVKqyW4Z34djFbW8/vU+axl8vtxfWrOH9kH+zDpH99qVUu4xrEck5yV15p+f77E2W5NPl/uBY5Us2XKQmWf10DFklFJudeuEJI4cr2HBuv1W1u/T5f6y81bg689LtJxEKdXWnJUQxaheUby4ZjdVte7fe/fZcj94rJIF3+znJyO6062jjteulHK/W8cnUVhWzZtpuW5ft8+W+4Jv9lPncDBnXB/bUZRSbdTo3p0Y0TOSlz/fQ73DuHXdPlnudfUO3t6Qx9i+0cRHtbcdRynVRokIvxiTQG5RJWt2undqC58s93c25nGotErvRlVKWTd5QBfiIkJ4bvVujHHf3rvPlbvDYXj2k2yGxndkYv9Gp29VSim3CfT349fj+pC2r5iv9hx123p9rtw37i8mr7iSa0f3RERsx1FKKS4f0Z2w4ADe3pDntnX6XLm/mZZLcIAfk1JibUdRSikAQgL9uWhQHMvTD1FS4Z4BxXyq3HccKuXtDXlcObKHTqGnlPIos85JoLK2nuc+zXbL+nyq3N/deAB/P+G2iUm2oyil1HekdA3n4sFdWbBuP7X1Dpevz6fKffWOQkYmRulQA0opjzRtUBfKqupYn1Pk8nX5TLlvzj3GrsLjXNBXr5BRSnmm85KiCfL34600159Y9Zly/+NbW+gaEcLlqd1tR1FKqUaFBgdww3mJLN50gM92Hnbpunyi3A8cq2RX4XFuOK+XHpJRSnm02yb2JTjAjzVa7qeX5jx+NTIxynISpZQ6taAAP4Z070javmKXrsfry7223sGybfm0D/InuUuY7ThKKXVaIxIiyThQwu7Dx122Dq8v90Xrc/koo4BrRycQ4O/1/xylVBtw6ZCutA/y5/p/r3fZOry+DTfvP0ZMWDBzpybbjqKUUk3SPy6cm8f3IedoBcXlNS5Zh8vKXUSmiEiWiGSLyFxXrSeroJR+ejhGKeVl+nUJB2DHoTKXfH+XlLuI+AP/AKYCKcCVIpLS2uupdxh2FRynX6yWu1LKu3x7jjDrUKlLvr+r9txHAtnGmD3GmBpgETC9tVey72g51XUO3XNXSnmdmLBgOrYP9K49d6Ab0HDSwDznsv9PRGaLSJqIpB0+3LzrPR0Gpg7swuDuHZsdVCmlbBARfjayB4O6R7jk+we45LtCYwOpf2cKEmPMS8BLAKmpqc2anqRPTAeev3pEc96qlFLW3T7FdReCuGrPPQ+Ib/C8O3DQRetSSin1Pa4q9/VAkogkikgQMBNY4qJ1KaWU+h6XHJYxxtSJyM3AR4A/8IoxJsMV61JKKfVDrjrmjjFmGbDMVd9fKaXUyXn9HapKKaV+SMtdKaV8kJa7Ukr5IC13pZTyQWJMs+4fat0QIoeBfS34Fp2BI60UpzVprjOjuc6cp2bTXGemubl6GmOiG/uCR5R7S4lImjEm1XaO79NcZ0ZznTlPzaa5zowrculhGaWU8kFa7kop5YN8pdxfsh3gJDTXmdFcZ85Ts2muM9PquXzimLtSSqnv8pU9d6WUUg1ouSullA/y6nJ31yTcTcySIyLbRGSziKQ5l0WJyMcissv5OdJNWV4RkUIRSW+w7KRZRORO5zbMEpHJbs51v4gccG63zSJykYVc8SKyWkQyRSRDRH7jXG51m50il9VtJiIhIrJORLY4cz3gXG57e50sl/WfMee6/EVkk4h86Hzu2u1ljPHKD04MJbwb6AUEAVuAFIt5coDO31v2ODDX+Xgu8JibsowFhgPpp8vCiQnMtwDBQKJzm/q7Mdf9wB8aea07c8UBw52Pw4CdzvVb3WanyGV1m3FiprUOzseBwDfAKA/YXifLZf1nzLm+3wELgA+dz126vbx5z90tk3C30HRgvvPxfGCGO1ZqjPkMKGpilunAImNMtTFmL5DNiW3rrlwn485c+caYjc7HZUAmJ+b8tbrNTpHrZNyVyxhjjjufBjo/DPa318lynYzbfsZEpDswDfjn99bvsu3lzeV+2km43cwAK0Rkg4jMdi6LNcbkw4lfVCDGWrqTZ/GE7XiziGx1Hrb59k9TK7lEJAEYxom9Po/ZZt/LBZa3mfMQw2agEPjYGOMR2+skucD+z9gzwO2Ao8Eyl24vby73007C7WZjjDHDganAHBEZazHLmbC9HZ8HegNDgXzgSedyt+cSkQ7AO8BtxpjSU720kWUuy9ZILuvbzBhTb4wZyon5kUeKyMBTvNx2LqvbS0QuBgqNMRua+pZGlp1xLm8ud4+ahNsYc9D5uRBYzIk/owpEJA7A+bnQVr5TZLG6HY0xBc5fSAfwMv/789OtuUQkkBMF+oYx5l3nYuvbrLFcnrLNnFmOAZ8CU/CA7dVYLg/YXmOAS0UkhxOHj8eLyOu4eHt5c7l7zCTcIhIqImHfPgYuBNKdeWY5XzYLeN9GPqeTZVkCzBSRYBFJBJKAde4K9e0Pt9NlnNhubs0lIgLMAzKNMU81+JLVbXayXLa3mYhEi0hH5+N2wERgB/a3V6O5bG8vY8ydxpjuxpgETvTUJ8aYq3H19nLVmWF3fAAXceIKgt3A3RZz9OLE2e0tQMa3WYBOwCpgl/NzlJvyLOTEn5+1nNgLuP5UWYC7ndswC5jq5lyvAduArc4f6jgLuc7lxJ+9W4HNzo+LbG+zU+Syus2AwcAm5/rTgftO9/NuOZf1n7EG67uA/10t49LtpcMPKKWUD/LmwzJKKaVOQstdKaV8kJa7Ukr5IC13pZTyQVruSinlg7TclVLKB2m5K6WUD/p/uYCF+NxScPgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from pylab import *\n",
    "%matplotlib inline\n",
    "\n",
    "plot(lngE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Thermod(T, lngE, Energies, N):\n",
    "    Z = 0.\n",
    "    Ev = 0.  # <E>\n",
    "    E2v = 0. # <E^2>\n",
    "    for i,E in enumerate(Energies):\n",
    "        w = exp(lngE[i]-lngE[0]-(E-Energies[0])/T) # g(E)/g0 Exp(-(E-E0)/T)\n",
    "        Z += w\n",
    "        Ev += w*E\n",
    "        E2v += w*E**2\n",
    "    Ev *= 1./Z\n",
    "    E2v *= 1./Z\n",
    "    cv = (E2v-Ev**2)/T**2\n",
    "    return (Ev/(N**2), cv/(N**2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "Te = linspace(0.5,4.,300)\n",
    "\n",
    "Thm=[]\n",
    "for T in Te:\n",
    "    Thm.append(Thermod(T, lngE, Energies, N))\n",
    "Thm = array(Thm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEGCAYAAABmXi5tAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsx0lEQVR4nO3deXxV1b3//9cnA5kISQiBMIUwySigRgYVilOdsbW02FoV2yuda39+v99brfd7W++9/dbba++trW0tdaqz1jqglTrWoVXQgKDMQpgCBELmkDln/f7YBwjxhAznJOecnPfz8TjuffZeZ+9PtuRzVtZeey1zziEiIv1fXLgDEBGRvqGELyISI5TwRURihBK+iEiMUMIXEYkRCeEO4GSGDBni8vPzwx2GiEjUWLNmzWHnXE6gfRGd8PPz8yksLAx3GCIiUcPMdne0T006IiIxQglfRCRGKOGLiMSIiG7DFxHpqubmZoqLi2loaAh3KH0iOTmZUaNGkZiY2OXPKOGLSL9QXFxMeno6+fn5mFm4w+lVzjnKysooLi5m7NixXf6cmnREpF9oaGggOzu73yd7ADMjOzu723/NKOGLSL8RC8n+qJ78rEr4Ij3hHOxZDe/eDQc3hjsakS5RG75IdzkHb/w7vPOL49su+TnM+Ub4YhLpAtXwRbrro6e8ZH/6dXDTRzD5clj5z7DtlXBHJmEWHx/PrFmzjr3uuOOOY/sWL15MUVERc+bMYdasWeTl5ZGTk3Os7K5du7jggguoqKjotfhUwxfpjvpKeOU2GFkAl98FcXGw+H64Zz689L8gfzUMSA13lBImKSkprFu37lPbN27cSGtrK+PGjWP16tUAPPjggxQWFnL33XcfK3fttdfy29/+lttuu61X4lPCF+mOD/4AR0rhmj95yR4gIQkuuxP+eAUU3gdnfS+8MQq3v7CRTfurQ3rMqSMG8eMrpvXos48++ihXXnllp+UWLVrE/Pnzey3hq0lHpKuaG2D1cphwAYw47cR9YxdA/nx47zfQ0hie+CTs6uvrT2jSefLJJwH4xz/+wRlnnNHp57OysmhsbKSsrKxX4lMNX6SrNq+AI4dg3ncC7z/nB/DIF2DjszDz6j4NTU7U05p4sDpq0jlw4AA5OQFHLP6UoUOHsn//frKzs0McnWr4Il237lHIHANjFwbeP/58yBoLHz7Sl1FJFEhJSenyQ1INDQ2kpKT0ShxK+CJdUbUPit6CmV8+3nbfnhnMugZ2vQMVu/o0PIlsU6ZMYfv27Z2Wc85RUlJCb038pIQv0hWbngcczPjSycvN+jJgsO7xvohKIkz7NvxbbrkFgMsuu4w333yz08+vWbOGuXPnkpDQO63tasMX6YptKyFnCmSPP3m5jFEwbiGseww+88OO/xqQfqm1tTXg9sWLF3Puuedy++23Ex8fD8DSpUtZunTpCeUefvhhvv3tb/dafPrXKNKZ+grY9Q+YdEnXys+6Bqr2eE07Inht+Lfffjv79u07abnp06dz/vnn91ocIUn4Zna/mR0ysw0d7F9oZlVmts7/+tdQnFekT2x/HVxr1xP+lMthwEDY8HTvxiVR5aKLLiIvL++kZW688cZejSFUNfwHgYs7KfOOc26W//VvITqvSO/buhLScmBk5/2oAUhMgcmXwaYV0NLUu7GJdENIEr5z7m2gPBTHEokorc3wyasw8SKIi+/656YvhoZK2PFGr4Um0l192YY/z8zWm9lKM+vwqQgzW2ZmhWZWWFpa2ofhiQSw+11orOp6c85R4xZCShZs+HOvhCXSE32V8NcCY5xzM4FfA891VNA5t9w5V+CcK+jqk2kivWbbXyE+Ccaf273PJQyAKYtg60vQVNc7sYl0U58kfOdctXOu1r/+EpBoZkP64twiPeYcbPkLjPsMDEjr/uenfwGaauETDZssgf3gBz/g7bff5vOf/zyzZs1iwoQJZGRkHOvH/+6773L11VfzySefhOR8fdIP38xygYPOOWdms/G+aHpndCCRUCndApW7vTFyeiL/HBg4zOutM+1zoYxM+oHy8nJWrVrFL3/5SxYsWADAm2++yZ133smLL754rFxzczM///nP+cMf/hD0OUOS8M3scWAhMMTMioEfA4kAzrl7gMXAt8ysBagHrnbOuVCcW6TXbF3pLU/prANaB+LiYdrnofABaKiG5EGhi01ObuUtUPJxaI+ZeypcckenxR566CHuvPNOzIxx48bx4YcfUlRURFxcHHV1dUyaNImioiKefvppLr64839b8+fPZ+nSpbS0tAT9BG6oeul82Tk33DmX6Jwb5Zy7zzl3jz/Z45y72zk3zTk30zk31zn3bijOK9Krtq6E4bNg0IieH2P6F6C10WvLl35v48aN/PSnP+WNN95g/fr13HfffcycOZO33noLgBdeeIGLLrqIxMTELg+ZHBcXx4QJE1i/fn3Q8WloBZFAakuh+ANYeGtwxxl1JmTkeb11NGRy3+lCTbw3vPHGGyxevJghQ7xblIMHD2bJkiU8+eSTnHvuuTzxxBPHhk7oyZDJXfmCOBkNrSASyCcvAw4m9bA55ygzmP55rz9+nR5V6e+cc5jZCdsWLVrEypUrKS8vZ82aNZx33nlAeIZMVsIXCWTrShg0EnJnBH+s6V8AX4t/xE3pz84//3yeeuqpYzNWlZeXM3DgQGbPns1NN93E5ZdffmzwtK4OmQywbds2pk0LflIXJXyR9pobvBr5KRd7NfRg5c6A7Il6CCsGTJs2jdtuu43PfOYzzJw5k5tvvhmAJUuW8Mgjj7BkyZJjZbs6ZPLBgwdJSUlh+PDhQcenNnyR9na+Dc11MOnS0BzPzKvlv/WfUFMC6bmhOa5EpOuvv57rr7/+hG2LFy+mfcfE+fPnc+utt1JZWUlmZiYACxcuZOHChSeUe+yxx/jGN74RkthUwxdpb8uL3miX+eeE7pjTrwIcbHwudMeUqPeLX/yCPXv2nLRMZmbmp75AekoJX6QtX6v3dO3Ez0JicuiOmzMJhp2qZp1eFm2P98yZM4cZM05+n+iGG24I2P++Jz+rEr5IW3tXQ91hmHJF6I89/Soofh8qdof+2EJycjJlZWVRl/R7wjlHWVkZycndq5SoDV+krc0veIOlTbww9MeefhW8fjt8/CdY8L9Df/wYN2rUKIqLi4mVUXaTk5MZNWpUtz6jhC9ylHOw+UVvZMyk9NAfPysf8ufD2ofgnJs1322IJSYmMnbs2HCHEdH0L07kqAPrvbloJ1/ee+c4Y6k3IFuRJkaRvqeEL3LUpufA4rs/2Ul3TLkCUrO9AdVE+pgSvgiAzwcfPw3jz4O0XpyqISEJZl3jPclbfaD3ziMSgBK+CMCed6FqL8xY0nnZYJ2xFFwrrP1j759LpA0lfBGAj56ExDSYHKKna08me7zXz//9P0Bzfe+fT8RPCV+kuQE2Pg9TF/VsKsOeOPsmr7//usf65nwiKOGLeEMhN1bBjC/13TnHnA0jTof37vae7hXpAyFJ+GZ2v5kdMrMNHew3M/uVmW03s4/M7PRQnFckJAofgPQRMPYzfXdOMzj7+1Be5I3dI9IHQlXDfxA42UwRlwAT/a9lwO9CdF6R4BzeDkV/g4KveXPQ9qUpi2DwOHjr514vIZFeFqo5bd8GTjadz5XAQ86zCsg0s+AHdxYJVuF9EJcIp1/X9+eOi4eFP4KDG2DjM31/fok5fdWGPxLY2+Z9sX/bp5jZMjMrNLPCWBkTQ8Kk6Qh8+Kh3szZ9WHhimP4FGDoN/vZTaG0OTwwSM/oq4QeaNijgkHbOueXOuQLnXEFXJ/gV6ZGPnvJu1p55Y/hiiIuD8/7Fa8v/8JHwxSExoa8SfjEwus37UcD+Pjq3yKe1tsA/7oLhsyBvbnhjmXQJjJ4Lb/wH1FeENxbp1/oq4a8ArvP31pkLVDnn9Fy5hM+GP0PFTljwf0Izb20wzODSn0N9Ofzt/4U3FunXQjI8spk9DiwEhphZMfBjIBHAOXcP8BJwKbAdqANuCMV5RXrE54N37oShU0M3b22whs+Egq/DB/fCaV/13ouEWEgSvnPuy53sd8B3QnEukaBtfh4Ob4PF90fWmPTn3QYbn4UV34d/eg3iE8MdkfQzEfSvXaQPtDTB6/8OQybB1M+FO5oTpWTB5f8NB9bBO78IdzTSDynhS2wpvB/Kd8Bn/73vH7TqiqlXwqlf8h7G2rc23NFIP6OEL7GjvgLeugPGLfRGq4xUl/4cBg6DP38d6ivDHY30I0r4EjvevMNLoJ/9j/D3zDmZlCz44gNQuQee/aaGXZCQUcKX2LBvDby/3BszJ/fUcEfTuby5cNH/g20rvR5FIiEQkl46IhGttRlW3ARpQ+GCH4c7mq6bvcxrx//bTyEzD2ZeHe6IJMop4Uv/995v4ODH8KWHITkj3NF0nRks+jXU7IfnvwNpOTDh/HBHJVFMTTrSv5Vs8GrIky+HKVeEO5ruSxgASx6BnMnwxDVQ9Fa4I5IopoQv/VdzAzxzo3cT9IpfRfaN2pNJzoBrn4PBY+GxL8GON8IdkUQpJXzpv177MRzaBFf+BtKywx1NcAbmwPUvQPYEePRL8NGfwh2RRCElfOmfPn4aVt8Dc74FEy8MdzShkTYElr4Io+fAM/8Eb/2XumxKtyjhS/9zaLM3Hs3ouXDhv4U7mtBKyYJrn/Gexv3bf8DjS6DuZJPNiRynhC/9S0M1PPlVGJAGX3zQu+nZ3yQkwVXL4dI7oehNuGc+7H0/3FFJFFDCl/7D1wrPfgPKd3rJflA/njbZDGbfCF972Rvx8/6L4OXboKku3JFJBFPCl/7BOVj5z7D1Jbj4Dsg/O9wR9Y2Rp8M3/+5Nwv7e3fC7ebDjb+GOSiKUEr70D/+4y5s85Kzvw5xl4Y6mbyVnwBV3wdK/gMXBw5+Dx66G0m3hjkwijBK+RL+P/uR1wZz+Bbjg9nBHEz7558C33oXzfwy7/g6/nQsv3gxVxeGOTCJESBK+mV1sZlvNbLuZ3RJg/0IzqzKzdf7Xv4bivCJseMZrtx9zDnzud5E1g1U4JKbA/Jvh+x/CGUth7R/hrlnwwk1QsSvMwUm4Bf3bYWbxwG+AS4CpwJfNbGqAou8452b5X/2sr5yExcdPw5//yeuX/pUnvN4r4hmY482e9b21cPq1sO4x+NXp8Oy34MD6cEcnYRKK6tBsYLtzrsg51wQ8AVwZguOKdOyjp7xhE/LmwTV/gqT0cEcUmbLGwOX/Azet90bf3PQc/H4B3HcRbPizN5KoxIxQJPyRwN4274v929qbZ2brzWylmU0LwXklFjkHf/+ll+zHnA3XPAVJA8MdVeQbNAIuuQNu3uSNs19bAk9/Df5nujcxTOXezo8hUS8UCT/QiFSu3fu1wBjn3Ezg18BzHR7MbJmZFZpZYWlpaQjCk36jpdFriz56g/aap70HrKTrUrJg3nfgex/CV56C3Onw5s/gl6fCHxfB+ifVl78fM+fa5+ZuHsBsHvAT59xF/ve3AjjnfnaSz+wCCpxzh0927IKCAldYWBhUfNJPVO6Bp66H/WvhnJvhvP+rG7ShUrEL1j/htfNX7oYB6TDtSu9LNX8BxGvajGhiZmuccwUB94Ug4ScA24DzgX3AB8BXnHMb25TJBQ4655yZzQaexqvxn/TkSviCc157/cr/461/7rfROa59NPD5YM97XuLf9Bw01ULKYJhyOUz9HIxdAPGJ4Y5SOnGyhB/0V7dzrsXMvgu8DMQD9zvnNprZN/377wEWA98ysxagHri6s2QvQvlOePlH3tOzo+d43S6zx4c7qv4rLs57Qjn/bLjsTtj+upf4NzwDax/ymoMmXAinXATjz4PUweGOWLop6Bp+b1INP0Y11sA7/+0NFRCXAAtv9dqd4+LDHVlsaq73kv/mFfDJq1Bf7j3RO2q2N/T0xM/CsOlqYosQvdqk05uU8GNMQxV8cB+s+i0cKYUZS+CCn3g9TCQy+Fq9idU/ecV7HVjnbU/N9npN5c/3nvjNmawvgDBRwpfIVr4T1jwAhQ9AY7XXXLDwRzD6zHBHJp2pOQjbX/OGctj1DlT5u3emDPaahkbPgZEFMHwmDEgNb6wxolfb8EV6pOkIbPur1zZc9KbXRDBlEZzz/8GIWeGOTroqfRicdo33AqjYDbv/cfwLYPML3naLh2HTYOQZMKrAWw45Rc10fUw1fOk79ZVekt/8glcrbGmAjDzv0f9Z10BGoOf1JKrVHoJ9a6C40FvuWwuNVd6+xFQYOtV7FmDYdMidAcOm6qnpIKlJR8Kn9hBs+YuX5He+Bb4WSB/hda2ccgWMOUu1vFji80H5Du8LoOQjKPnYezVUHi+TNfb4l8CQU7z7AdnjNVZSF6lJR/pW5V7Y8qKX5He/Czjvl3jed7xmmxGn64ZerIqLgyETvRdf9rY5B9X7/Ml/g/dFcHADbH6RYw/tWzwMHgtDJkGO/0tgyCneS0NrdJkSvoRG5V7Y+Kz32r/W2zZ0Knzmh15Nftg0b1o+kfbMIGOU95p0yfHtTXVQth1Kt8LhrVC6xZvU5ZOXvb8Ujxo0yvsyGDzO+0tg8DjvlZWvoTfaUcKXnqs+cPzBnGL/JNrDZ3kTcExZBEMmhDM6iXYDUmH4DO/VVmszlBcd/yI4/In3fsuLUFd2Ytn04f4vAP8XQtZYyMyDjNGQlhNzf2kq4Uv31JbC5udhw7Nebwyc19Z63v+F6Vd5v1QivSk+EXImea/26iuhYqf3BVBe5HX5LS/yHhirPdjuOAOO/2WRkectM0f734+G9Nx+9xeCEr50rq7cqz1teAZ2vg2u1Ws7XXgLTLvKa1MViQQpmZByGow47dP7Gmu9geKq9nrTPlbu8ZZVe2HH61BTwqcG+k0aBAOHecn/ZMvkjKhoslTCl8Aaqr0xbDY8AzveAF+z9+fwOT/wkrza5CXaJA30ev/kTg+8v6XJu3l89Auh5oD3YFltibfcV+gtW+o//dm4RO9p49Rsb4yhY+vZgbenZHl/PfTx75ASvhzXUO09Lr/xWe9P4NZG70/bud/0hsodPktJXvqvhAH+tv6xHZdxznsavO0XQW0JHDns3T+oK/eWBzd6y/oKPj09iF9cAiRnen8dpGSeuD4wFxb+MPQ/YsiPKNGlvgK2roRNK7w/a1ubvH9sBV/z2uRHFsTcjS2RDpl5STk5o2tNmb5W775Cvf+LoK7M+3Kor/DGjmqo9PYfXa/Y5a0PSFPClxAp3+kl9y1/8drkfS1e17Yzb4Spi7xREJXkRYIXFw9p2d6LiV3/XC89EKuEHwsaa70eNdtf84a5Ld/hbT/6MNTUK72HodRcIxIZeul3UQm/P6raB3tXwZ7V3rJkg9ezJjHVG7p2zjdg/PneQypK8iIxQwk/mvl8ULnr+HgkJR/DgY+gZr+3PzHVG5Vw/s3eWOV58yAxOawhi0j4KOFHOue8mzllRd5j5ie8dkDzEa+cxXsPooyd7/VBHj0Hck/VHKQickxIEr6ZXQzchTen7b3OuTva7Tf//kuBOmCpc25tKM4d1Vqbj9+1P3IIqvd7zTFVe/39gfd5y6ba45+xOMgcA9kTvFr70CleYh86VbV3ETmpoBO+mcUDvwEuBIqBD8xshXNuU5til+Ddop4IzAF+519GD+e83iwtjV7XxZZGr596S9PxZUu9d4O0scbrq9tY0+5V7e+ve9ibwq+hKvC50oZ6j3fnnOLN/pQx0j8w1ATvRmvCgL792UWkXwhFDX82sN05VwRgZk8AVwJtE/6VwEPOG3x/lZllmtlw59yBEJz/0+45x5t42TlwPsC/dLR5335fgLJH148m+o4eoDgZi/MmdEga5C1Ts70aeVoOpA7xd9nyrw8a4b007reI9IJQJPyRwN4274v5dO09UJmRwKcSvpktA5YB5OXl9SyiodO8WrjFeb1QLA6w4+8x//YA+469b7MeFwcJyRCf5NWuT1gmeYMwJfjXkwYdT+5J6ZCYop4wIhIRQpHwA2Wz9lXhrpTxNjq3HFgO3oxXPYroqt/36GMiIv1ZKB6nLAZGt3k/CtjfgzIiItKLQpHwPwAmmtlYMxsAXA2saFdmBXCdeeYCVb3Wfi8iIgEF3aTjnGsxs+8CL+N1y7zfObfRzL7p338P8BJel8zteN0ybwj2vCIi0j0h6YfvnHsJL6m33XZPm3UHfCcU5xIRkZ7RkIgiIjFCCV9EJEYo4YuIxAglfBGRGKGELyISI5TwRURihBK+iEiMUMIXEYkRSvgiIjFCCV9EJEYo4YuIxAglfBGRGKGELyISI5TwRURihBK+iEiMUMIXEYkRSvgiIjEiqBmvzGww8CSQD+wCvuScqwhQbhdQA7QCLc65gmDOKyIi3RdsDf8W4HXn3ETgdf/7jpzrnJulZC8iEh7BJvwrgT/61/8IfC7I44mISC8JNuEPc84dAPAvh3ZQzgGvmNkaM1t2sgOa2TIzKzSzwtLS0iDDExGRozptwzez14DcALtu68Z5znbO7TezocCrZrbFOfd2oILOueXAcoCCggLXjXOIiMhJdJrwnXMXdLTPzA6a2XDn3AEzGw4c6uAY+/3LQ2b2LDAbCJjwRURiiXOOirpmSqoaOFjdQEl1Ay2tPq6dlx/ycwXVSwdYAVwP3OFfPt++gJmlAXHOuRr/+meBfwvyvCIiEa+l1cehmkZKqhsoqWo4ltQPVDUc31bdQFOL74TPZaYmRmTCvwN4ysy+DuwBvghgZiOAe51zlwLDgGfN7Oj5HnPO/TXI84qIhJVzjrIjTeyrqGdfZT3FFXXsq6jnQJukfri2EV+7hukB8XHkZiSTOyiZWaMzj63nZnivYYOSyRmY1CsxB5XwnXNlwPkBtu8HLvWvFwEzgzmPiEhfa/U5DlY3sK+ynn0V/oReWU+xP8Hvr6ynofnEmvnApARGZHpJe1Juuj+Rp5CbkcSwQckMz0ghKzURfwW4zwVbwxcRiVp1TS3sLqvzv46wu7yOPWV17C4/woHKBlraVc+z0wYwMiuFScPSOW/SUEZmpTAyM4WRWSmMykxlUEpC2JJ5Vyjhi0i/VlXfzM7DR7yEXlbHrrIj/qReR2lN4wlls1ITyctO47TRWVwxI+VYQh+VlcKIzBRSB0R3yozu6EVEAJ/Psa+ynh2ltewoPeItD3nrh2tPTOq5g5IZk53KuZNyGJOdxpjsVMYMTiMvO5WMlMQw/QR9QwlfRKJGS6uPnYePsKWkhk8O1R5L7DsPH6GxTU+XjJREJgwdyHmTcxifM5CxQ9IYOySN0YNTSU6MD+NPEF5K+CIScZxzHKxuZEtJNVtLathaUsPmkhp2HKqlqdVL7HEGowenMj5nIPMnDmFczkDG5wxkfE4ag9MGRHRbergo4YtIWDW1+Nh2sIaP91Wx+UA1W/wJvqq++ViZXH+vlwUThzApN51ThqUzYejAmK6t94QSvoj0mYbmVraWeMl94/4qPt5XxdaSGppbvd4wA5MSmJSbzmUzhjM5N51Jw9KZlJtOZuqAMEfePyjhi0ivaGn1sfVgDR/uqeSj4ko+3lfNJwdrjnV1zEhJ5NSRGXz9nHFMHzmIU0dmkDc4VU0xvUgJX0RCovxIEx/uqWDtngrW7q5kfXEldU2tgNfdcfrIDM6dNI5TR2YwfWQGo7JSlNz7mBK+iHSbc46iw0dYXVTOmt0VfLingqLDRwCIjzOmDE9n8RmjOD0vi9Pzshg9WMk9Eijhi0injib4VUVlrCoqZ1VR2bGHlrLTBnBaXhaLC7wEP2NURtQ/oNRf6f+KiAS0p6yOd7aXfirBD01P4qzx2cwdl82csYMZOyRNtfcooYQvIoA3rsyqojLe3naYt7aVstPfRNM2wc8dl01+tm6sRislfJEYtqO0ltc3H+StbaV8sLOCplYfyYlxzB2XzXXzxjB/Yg7jc1SD7y+U8EViiM/nWF9cySubDvLKxhJ2lHq1+FOGDeT6s8aw4JQczswfrAea+iklfJF+rqnFx3tFZbyysYRXNx3kUE0jCXHGnHGDuW5ePhdMHcbIzJRwhyl9QAlfpB/y+Ryrd5azYv1+Vm44QGVdM6kD4lk4KYfPTs3l3ElDyUjt3yNDyqcFlfDN7IvAT4ApwGznXGEH5S4G7gLi8aY+vCOY84rIpznn+HhfFSvW7eeFj/ZzsLqR1AHxfHbqMC6fMYJzJg5RU02MC7aGvwG4Cvh9RwXMLB74DXAhUAx8YGYrnHObgjy3iABltY08s3YfTxbuZfuhWhLjjc+cMpR/uWwE508Zqj7xckywc9puBjq7gz8b2O6f2xYzewK4ElDCF+mhVp/jnU9KeapwL69uOkhzq+P0vEx+dtWpXDp9uJprJKC++OofCext874YmNNRYTNbBiwDyMvL693IRKJMTUMzTxUW8+C7O9lbXk9WaiLXz8tnyZmjmTgsPdzhSYTrNOGb2WtAboBdtznnnu/COQJV/12Abd4O55YDywEKCgo6LCcSS3aXHeHBd3fxp8JiahtbODM/ix9ePJkLpw4jKUHt8tI1nSZ859wFQZ6jGBjd5v0oYH+QxxSJCRv3V/Hr17fz8qYS4s24YuYIbjg7nxmjMsMdmkShvmjS+QCYaGZjgX3A1cBX+uC8IlFrw74q7nr9E17ddJD05AS+vXA8183LZ9ig5HCHJlEs2G6Znwd+DeQAfzGzdc65i8xsBF73y0udcy1m9l3gZbxumfc75zYGHblIP7Tz8BHuWLmZlzceZFByAj+4YCI3nD2WjBTdhJXgmXOR20xeUFDgCgsDdu0X6Veq6pq56/VPeOi9XSQlxLFswXhuOCefQclK9NI9ZrbGOVcQaJ866IqEUavP8ejq3fz3q9uorm9myZmjufnCSeSkJ4U7NOmHlPBFwmRrSQ0//PNHrNtbydkTsvmXy6YyZfigcIcl/ZgSvkgfa2718es3tvO7N7eTnpzIXVfPYtHMERqCWHqdEr5IH9pbXsf3Hv+QdXsr+dysEfzrFdMYnDYg3GFJjFDCF+kjK9bv57ZnPgaDu79yGpfPGBHukCTGKOGL9LJWn+NnL23m3r/v5IwxWfxyySxGD04Nd1gSg5TwRXpRVX0z33/8Q97aVsrSs/K57bIpJMbHhTssiVFK+CK9ZG95Hdc/8D57y+v42VWn8uXZGgxQwksJX6QXbDtYw7X3raa+qZVHvj6HOeOywx2SiBK+SKh9uKeCGx78gAHxcTz1zXlMzlXfeokMSvgiIfT+znKWPvA+OelJPPy1OeRl6+asRA4lfJEQWbO7ghseeJ/hGck8fuNchmpkS4kw6i4gEgIfFVey9H6vZv+Ykr1EKCV8kSBt3F/FV+9dTWZaIo/dOFdj1kvEUsIXCcK2gzV89d7VDExK4LF/msuIzJRwhyTSISV8kR7adfgI19y7msT4OB67ca6enpWIp5u2Ij2wr7Kea+5dTavP8eSyueQPSQt3SCKdUg1fpJsOVTdwzR9WUd3QzENfm83EYenhDkmkS4JK+Gb2RTPbaGY+Mws4pZa/3C4z+9jM1pmZ5iyUqFVxpImv3reaQzWNPHjDbKaPzAh3SCJdFmyTzgbgKuD3XSh7rnPucJDnEwmbqvpmrrv/fXaV1fHg0jM5Y0xWuEMS6ZagEr5zbjOgmXqk3ys/0sS1961m28Eafn/tGZw1YUi4QxLptr5qw3fAK2a2xsyWnaygmS0zs0IzKywtLe2j8EQ6dqi6gSW/f4/th2pZfm0B500eFu6QRHqk0xq+mb0G5AbYdZtz7vkunuds59x+MxsKvGpmW5xzbwcq6JxbDiwHKCgocF08vkiv2Ftex7X+NvsHbjiTs8arZi/Rq9OE75y7INiTOOf2+5eHzOxZYDYQMOGLRIoPdpXzjYfX0NLq4+Gvz1GbvUS9Xm/SMbM0M0s/ug58Fu9mr0jEemZtMdf8YTUZKYk8952zleylXwi2W+bnzawYmAf8xcxe9m8fYWYv+YsNA/5uZuuB94G/OOf+Gsx5RXpLQ3MrP1mxkZufWk9BfhbPfvssxuUMDHdYIiERbC+dZ4FnA2zfD1zqXy8CZgZzHpG+sP1QDd97fB2bD1Rzw9n5/OhSzT8r/YuGVpCY1+pzPPzeLv7zr1tJGRDPfdcXcP4U9cSR/kcJX2Laxv1V/OjZDazfW8mCU3L4r8UzNLyx9FtK+BKTyo808avXP+HhVbvJSk3krqtnsWjmCD1EKP2aEr7ElPqmVh54dye/+9sOjjS1cPXsPP75oklkpg4Id2givU4JX2JCTUMzj6zaw31/L+JwbRMXTBnKLZdMZsJQjXQpsUMJX/q1stpG/vjuLh58dxfVDS0sOCWH7547gdljB4c7NJE+p4Qv/Y5zjjW7K3hk1W5e+riEplYfF0/L5dvnjmfGqMxwhycSNkr40m8crm3kxfX7eeKDvWwpqSE9KYEvzx7NtfPGqOlGBCV8iXK1jS28srGE59ft5+/bD9Pqc0wbMYifXXUqi2aOIC1J/8RFjtJvg0SdfZX1vL75IK9tPsSqHWU0tfoYmZnCNxaM48pZI5mUq9q8SCBK+BLx6ptaKdxdzns7ynhjyyG2lNQAMHZIGtfNG8NF03M5Iy+LuDj1oRc5GSV8iTgNza2s3V3Be0VlvLejjPXFlTS3OhLijNPHZPGjSydz/pRhjNegZiLdooQvYeXzOYoO17JubxXr91aybm8lmw9U0+JzxBmcOiqTr58zjrnjBnNm/mC1yYsEQb890mdafY6dh4+wpaSaTfurWV9cyUd7q6hpbAFgYFICM0ZlcOOCcczOH0xBfhbpyYlhjlqk/1DCl15RVdfM5pJqNh+oZsuBGjaXVLO1pIbGFh8A8XHG5Nx0Fs0awczRmZw2OpNxOQOJVzu8SK9Rwpcea271sbe8jqLSIxQdrvWW/vXDtU3Hyg1OG8CU4el8de4YpgwfxOTcdCYOG0hSQnwYoxeJPUr40iGfz1Fa20hxRR3FFfXsq6z3lhX17C2vY095HS2+4/PMD04bwLghaZw3eSjjcgYyOTedqcMHkZOepFEoRSJAUAnfzP4LuAJoAnYANzjnKgOUuxi4C4gH7nXO3RHMeaXnGltaqa5vobqhmYojTRyubaS0xv+qbfIvGzns39bU6jvh81mpiYzMSuGUYelcPD2XcTkDGZeTxrghaRpxUiTCBVvDfxW41TnXYmb/CdwK/LBtATOLB34DXAgUAx+Y2Qrn3KYgz90v+HyOZp+PllbnvXw+WnyO5lbfsffN/n3NPh+NzT4amlupb26lvslbNvhf3jYf9c2tHGlsoaq+meqGZqrrm6luaKG6vvlYG3p7ZpCdNoAhA5PISU9i/JA0cgYlMSorlVGZKYzMSmFkZop6yYhEsWDntH2lzdtVwOIAxWYD2/1z22JmTwBXAr2W8C//9TvUN7VyrLHBcWzdOYcDnDu6y+Hc8fdHy/g/5u3zf9qdcBz/p48dJ8Dn2rynXblmn6Ol1UebFpGgxccZqYnxJCXGk56cwKDkBAalJDIiI4VBKQkMSk5kUErise0ZKYnkpHsJfnDqABI0f6tIvxbK6trXgCcDbB8J7G3zvhiY09FBzGwZsAwgLy+vR4FMHJpOU4sPDI62HJtZm3Vv+9F2ZfP/52gJO+Fz3vajTdBmxz7R5jj+fW3LtTl+2+McXU+INxLj4rxlfBwJcUZ8nH+9zb6E+DgS47xlQryRnBBPyoB4UhLjSU6M85b+95pwW0ROptOEb2avAbkBdt3mnHveX+Y2oAV4NNAhAmzrsF7rnFsOLAcoKCjoUf33f5bM6snHRET6tU4TvnPugpPtN7PrgcuB851zgRJ0MTC6zftRwP7uBCkiIsELqg3A3/vmh8Ai51xdB8U+ACaa2VgzGwBcDawI5rwiItJ9wTb63g2kA6+a2TozuwfAzEaY2UsAzrkW4LvAy8Bm4Cnn3MYgzysiIt0UbC+dCR1s3w9c2ub9S8BLwZxLRESCo24dIiIxQglfRCRGKOGLiMQIJXwRkRhhgbvORwYzKwV29/DjQ4DDIQynN0VTrBBd8UZTrBBd8UZTrBBd8QYT6xjnXE6gHRGd8INhZoXOuYJwx9EV0RQrRFe80RQrRFe80RQrRFe8vRWrmnRERGKEEr6ISIzozwl/ebgD6IZoihWiK95oihWiK95oihWiK95eibXftuGLiMiJ+nMNX0RE2lDCFxGJEVGd8M3sYjPbambbzeyWAPsXmlmVfyTPdWb2r+GI0x/L/WZ2yMw2dLDfzOxX/p/lIzM7va9jbBdPZ/FG0rUdbWZ/M7PNZrbRzG4KUCYirm8XY42ka5tsZu+b2Xp/vLcHKBMp17YrsUbMtW0TU7yZfWhmLwbYF9pr65yLyhcQD+wAxgEDgPXA1HZlFgIvhjtWfywLgNOBDR3svxRYiTdD2FxgdYTHG0nXdjhwun89HdgW4N9CRFzfLsYaSdfWgIH+9URgNTA3Qq9tV2KNmGvbJqabgccCxRXqaxvNNfxjk6M755qAo5OjRyTn3NtA+UmKXAk85DyrgEwzG9430X1aF+KNGM65A865tf71Grx5F0a2KxYR17eLsUYM//Wq9b9N9L/a9/SIlGvblVgjipmNAi4D7u2gSEivbTQn/ECTowf6xZnn/xNvpZlN65vQeqSrP08kibhra2b5wGl4tbu2Iu76niRWiKBr629yWAccAl51zkXste1CrBBB1xb4JfDPgK+D/SG9ttGc8LsyOfpavHElZgK/Bp7r7aCC0K3J3iNAxF1bMxsI/Bn4gXOuuv3uAB8J2/XtJNaIurbOuVbn3Cy8+ahnm9n0dkUi5tp2IdaIubZmdjlwyDm35mTFAmzr8bWN5oTf6eTozrnqo3/iOW/WrUQzG9J3IXZLVE32HmnX1swS8RLoo865ZwIUiZjr21mskXZtj3LOVQJvAhe32xUx1/aojmKNsGt7NrDIzHbhNUmfZ2aPtCsT0msbzQm/08nRzSzXzMy/Phvv5y3r80i7ZgVwnf+u/Fygyjl3INxBdSSSrq0/jvuAzc65/+6gWERc367EGmHXNsfMMv3rKcAFwJZ2xSLl2nYaayRdW+fcrc65Uc65fLz89YZz7qvtioX02gY1p204OedazOzo5OjxwP3OuY1m9k3//nuAxcC3zKwFqAeudv5b333NzB7H6yEwxMyKgR/j3VQ6GutLeHfktwN1wA3hiPOoLsQbMdcWr6Z0LfCxv/0W4EdAHkTc9e1KrJF0bYcDfzSzeLzk+JRz7sV2v2eRcm27EmskXduAevPaamgFEZEYEc1NOiIi0g1K+CIiMUIJX0QkRijhi4jECCV8EZEYEbXdMkX6mpllA6/73+YCrUCp//1s/5hOIhFL3TJFesDMfgLUOufuDHcsIl2lJh0RkRihhC8iEiOU8EVEYoQSvohIjFDCFxGJEUr4IiIxQt0yRURihGr4IiIxQglfRCRGKOGLiMQIJXwRkRihhC8iEiOU8EVEYoQSvohIjPj/AdKxx9ss9FqLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from pylab import *\n",
    "%matplotlib inline\n",
    "\n",
    "plot(Te, Thm[:,0], label='E(T)')\n",
    "plot(Te, Thm[:,1], label='cv(T)')\n",
    "xlabel('T')\n",
    "legend(loc='best')\n",
    "show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Homework\n",
    "\n",
    "* Estimate the error of the density of states $log(g(E))==S(E)$ for 32x32 Ising model. \n",
    "\n",
    "To do that, you can use five independent runs of WL agorithm, which will produce $S_i(E)$ for i=1..5. \n",
    "In each WL run use $10^9$ Monte Carlo steps. From the series of $S_i(E)$ you can calculate \n",
    "\\begin{eqnarray}\n",
    "&& <S(E)> = \\frac{1}{n} \\sum_{i=1}^n S_i(E) \\\\\n",
    "&& \\sigma_E = \\sqrt{\\frac{<S^2(E)>-<S(E)>^2}{n}} = \\sqrt{\\frac{1}{n^2} \\sum_{i=1}^n (S_i(E)-<S(E)>)^2}\n",
    "\\end{eqnarray}\n",
    "\n",
    "Plot $S_i(E)-<S(E)>$ and $\\sigma_E$. \n",
    "\n",
    "* According to the publication \"Exact Distribution of Energies in the Two-Dimensional Ising Model\", PRL 76, 78 (1995), the first 10 values for the density of states in 32x32 Ising model is:\n",
    "$g_{Exact}=[2,2048,4096,1057792,4218880,371621888,2191790080,100903637504,768629792768,22748079183872]$\n",
    "\n",
    "Plot $<S(E)>$ and $log(g_{exact})$ for the first 10 coefficients, as well as their difference.\n",
    "\n",
    "Verify if your estimated $<S(E)>$ is approximately $\\sigma_E$ from exact result. To do that, plot $<S(E)>-S_{exact}$ and $\\sigma_E$ for the first 10 energies.\n",
    "\n",
    "* Calculate thermodynamics from your best estimation of $g(E)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}