MCRG Attempt.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MCRG Code following Swendsen Description circa 1982"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The line_profiler extension is already loaded. To reload it, use:\n",
      "  %reload_ext line_profiler\n"
     ]
    }
   ],
   "source": [
    "from __future__ import division #safeguard against evil floor division\n",
    "\n",
    "import numpy as np\n",
    "from scipy import linalg as la\n",
    "import matplotlib.pyplot as plt\n",
    "import timeit\n",
    "%load_ext line_profiler\n",
    "\n",
    "np.set_printoptions(precision = 3,suppress=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2D Ising Model, periodic B.C., initial Hamiltonian contains only nearest neighbor interactions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\mathcal{H} = -\\beta H_\\Omega = K\\sum_{\\langle i,j\\rangle} \\sigma_i \\sigma_j + h\\sum_{i = 0}^{N_s} \\sigma_i$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Block spin transform, scale factor = b "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def assignBlockSpin(total):\n",
    "    '''Rule for assigning block spin value. Random tiebreaker'''\n",
    "    if total > 0:\n",
    "        s = 1;\n",
    "    elif total < 0:\n",
    "        s = -1;\n",
    "    else:\n",
    "        s = np.random.choice([-1,1])\n",
    "    return s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RGTransform(S,b):\n",
    "    '''Take a spin config S and produce renormalized block spin config that \n",
    "    groups b*b spins together into one block. Does not modify input S'''\n",
    "    L = S.shape[0];\n",
    "    assert L//b >= 2, \"Renormalized lattice will have linear dimension <=1\"\n",
    "    newS = np.empty([L//b, L//b],dtype=int)\n",
    "    for j in np.arange(L//b):\n",
    "        for i in np.arange(L//b):\n",
    "            block = S[(b*i):(b*i+b),(b*j):(b*j+b)]\n",
    "            total = np.sum(block)\n",
    "            newS[i,j] = assignBlockSpin(total);    \n",
    "    return newS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First 7 short range even couplings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def N1Coupl(S):\n",
    "    '''for spin field config S, \n",
    "    measure sum SiSj for i,j that are nearest neighbors, dir = (1,0)'''\n",
    "    L = S.shape[0];\n",
    "    assert L >=2, \"Lattice too small for NN interaction\"\n",
    "    val = 0.0;\n",
    "    for j in np.arange(L):\n",
    "        for i in np.arange(L):\n",
    "            val += S[i,j]*(S[i,(j+1)%L] + S[(i+1)%L,j])\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def N2Coupl(S):\n",
    "    '''for spin field config S, \n",
    "    measure sum SiSj for i,j that are next-nearest neighbors, dir = (1,1)'''\n",
    "    L = S.shape[0];\n",
    "    assert L >=2, \"Lattice too small for next nearest neighbor interaction\"\n",
    "    val = 0.0;\n",
    "    for j in np.arange(L):\n",
    "        for i in np.arange(L):\n",
    "            val += S[i,j]*(S[(i+1)%L,(j+1)%L] + S[(i-1)%L,(j+1)%L])\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def N3Coupl(S):\n",
    "    '''for spin field config S, \n",
    "    measure sum SiSj for i,j that are third-nearest neighbors,dir = (2,0)'''\n",
    "    L = S.shape[0];\n",
    "    assert L >=3, \"Lattice too small for third nearest neighbor interaction\"\n",
    "    val = 0.0;\n",
    "    for j in np.arange(L):\n",
    "        for i in np.arange(L):\n",
    "            val += S[i,j]*(S[i,(j+2)%L] + S[(i+2)%L,j])\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def N4Coupl(S):\n",
    "    '''for spin field config S, \n",
    "    measure sum SiSj for i,j that are 4th-nearest neighbors, dir = (2,1)'''\n",
    "    L = S.shape[0];\n",
    "    assert L >=3, \"Lattice too small for 4th nearest neighbor interaction\"\n",
    "    val = 0.0;\n",
    "    for j in np.arange(L):\n",
    "        for i in np.arange(L):\n",
    "            val += S[i,j]*(S[(i+1)%L,(j+2)%L] + S[(i+2)%L,(j+1)%L] \n",
    "                           +S[(i-1)%L,(j+2)%L] + S[(i-2)%L,(j+1)%L])\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def N5Coupl(S):\n",
    "    '''for spin field config S, \n",
    "    measure sum SiSj for i,j that are 5th-nearest neighbors,dir = (2,2)'''\n",
    "    L = S.shape[0];\n",
    "    assert L >=3, \"Lattice too small for third nearest neighbor interaction\"\n",
    "    val = 0.0;\n",
    "    for j in np.arange(L):\n",
    "        for i in np.arange(L):\n",
    "            val += S[i,j]*(S[(i+2)%L,(j+2)%L] + S[(i-2)%L,(j+2)%L])\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def PlaqCoupl(S):\n",
    "    '''for spin field config S, \n",
    "    measure sum SiSjSlSk for i,j,k,l that form 1 plaquette,'''\n",
    "    L = S.shape[0];\n",
    "    assert L >=2, \"Lattice too small to form plaquettes\"\n",
    "    val = 0.0;\n",
    "    for j in np.arange(L):\n",
    "        for i in np.arange(L):\n",
    "            val += S[i,j]*S[i,(j+1)%L]*S[(i+1)%L,(j+1)%L]*S[(i+1)%L,j]\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SubPlaqCoupl(S):\n",
    "    '''for spin field config S, \n",
    "    measure sum SiSjSlSk for i,j,k,l that form 1 sublattice plaquette,'''\n",
    "    L = S.shape[0];\n",
    "    assert L >=3, \"Lattice too small to form sublattice plaquettes\"\n",
    "    val = 0.0;\n",
    "    for j in np.arange(L):\n",
    "        for i in np.arange(L):\n",
    "            val += S[(i+1)%L,j]*S[i,(j+1)%L]*S[(i-1)%L,j]*S[i,(j-1)%L]\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def AllEvenCoupling(S):\n",
    "    '''for spin field config S, \n",
    "    Integrate measurement of first 7 even correlation functions in one vector'''\n",
    "    L = S.shape[0];\n",
    "    assert L >=3, \"Lattice too small to fit first 7 even couplings\"\n",
    "    val = np.zeros(7,dtype = float);\n",
    "    for j in np.arange(L):\n",
    "        for i in np.arange(L):\n",
    "            val += [S[i,j]*(S[i,(j+1)%L] + S[(i+1)%L,j]), #nearest neighbor (1,0)\n",
    "                    S[i,j]*(S[(i+1)%L,(j+1)%L] + S[(i-1)%L,(j+1)%L]), #next nearest neighbor (1,1)\n",
    "                    S[i,j]*(S[i,(j+2)%L] + S[(i+2)%L,j]), #3rd nearest neighbor (2,0)\n",
    "                    S[i,j]*(S[(i+1)%L,(j+2)%L] + S[(i+2)%L,(j+1)%L] \n",
    "                            +S[(i-1)%L,(j+2)%L] + S[(i-2)%L,(j+1)%L]),#4th nearest neighbor (2,1)\n",
    "                    S[i,j]*(S[(i+2)%L,(j+2)%L] + S[(i-2)%L,(j+2)%L]),#5th nearest neighbor (2,2)\n",
    "                    S[i,j]*S[i,(j+1)%L]*S[(i+1)%L,(j+1)%L]*S[(i+1)%L,j], # plaquette\n",
    "                    S[(i+1)%L,j]*S[i,(j+1)%L]*S[(i-1)%L,j]*S[i,(j-1)%L]] # sublattice plaquette\n",
    "    return val"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First 4 short range odd couplings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Magnetization(S):\n",
    "    '''for spin field config S, \n",
    "    measure sum Si = net magnetization of entire lattice'''\n",
    "    L = S.shape[0];\n",
    "    assert L >=1, \"Lattice linear dim < 1\"\n",
    "    val = np.sum(S);\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def TripletPlaq(S):\n",
    "    '''for spin field config S, \n",
    "    measure sum SiSjSk for i,j,k that lie on 1 plaquette,'''\n",
    "    L = S.shape[0];\n",
    "    assert L >=2, \"Lattice too small to have 3 spins on a plaquette\"\n",
    "    val = 0.0;\n",
    "    for j in np.arange(L):\n",
    "        for i in np.arange(L):\n",
    "            val += S[i,j]*S[(i+1)%L,j]*S[(i+1)%L,(j+1)%L];\n",
    "            val += S[i,j]*S[i,(j+1)%L]*S[(i-1)%L,(j+1)%L];\n",
    "            val += S[i,j]*S[(i-1)%L,j]*S[(i-1)%L,(j-1)%L];\n",
    "            val += S[i,j]*S[i,(j-1)%L]*S[(i+1)%L,(j-1)%L];\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def TripletAngle(S):\n",
    "    '''for spin field config S, \n",
    "    measure sum SiSjSk for i,j,k that lie at an angle,'''\n",
    "    L = S.shape[0];\n",
    "    assert L >=3, \"Lattice too small to have 3 spins at an angle\"\n",
    "    val = 0.0;\n",
    "    for j in np.arange(L):\n",
    "        for i in np.arange(L):\n",
    "            val += S[i,j]*S[(i+1)%L,j]*S[(i+2)%L,(j+1)%L];\n",
    "            val += S[i,j]*S[i,(j+1)%L]*S[(i-1)%L,(j+2)%L];\n",
    "            val += S[i,j]*S[(i-1)%L,j]*S[(i-2)%L,(j-1)%L];\n",
    "            val += S[i,j]*S[i,(j-1)%L]*S[(i+1)%L,(j-2)%L];\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def TripletRow(S):\n",
    "    '''for spin field config S, \n",
    "    measure sum SiSjSk for i,j,k that lie in a row,'''\n",
    "    L = S.shape[0];\n",
    "    assert L >=3, \"Lattice too small have 3 spins in a row\"\n",
    "    val = 0.0;\n",
    "    for j in np.arange(L):\n",
    "        for i in np.arange(L):\n",
    "            val += S[i,j]*(S[(i+1)%L,j]*S[(i+2)%L,j] \\\n",
    "                           + S[i,(j+1)%L]*S[i,(j+2)%L])\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def AllOddCoupling(S):\n",
    "    '''for spin field config S, \n",
    "    Integrate measurement of first 4 odd correlation functions in one vector'''\n",
    "    L = S.shape[0];\n",
    "    assert L >=3, \"Lattice too small to fit first 4 odd couplings\"\n",
    "    val = np.zeros(4,dtype = float);\n",
    "    for j in np.arange(L):\n",
    "        for i in np.arange(L):\n",
    "            val += [0, #magnetization\n",
    "                    S[i,j]*S[(i+1)%L,j]*S[(i+1)%L,(j+1)%L]+\n",
    "                    S[i,j]*S[i,(j+1)%L]*S[(i-1)%L,(j+1)%L]+\n",
    "                    S[i,j]*S[(i-1)%L,j]*S[(i-1)%L,(j-1)%L]+\n",
    "                    S[i,j]*S[i,(j-1)%L]*S[(i+1)%L,(j-1)%L], #3 spin plaquette\n",
    "                    S[i,j]*S[(i+1)%L,j]*S[(i+2)%L,(j+1)%L]+\n",
    "                    S[i,j]*S[i,(j+1)%L]*S[(i-1)%L,(j+2)%L]+\n",
    "                    S[i,j]*S[(i-1)%L,j]*S[(i-2)%L,(j-1)%L]+\n",
    "                    S[i,j]*S[i,(j-1)%L]*S[(i+1)%L,(j-2)%L], # 3 spin angle\n",
    "                    S[i,j]*(S[(i+1)%L,j]*S[(i+2)%L,j] + \n",
    "                            S[i,(j+1)%L]*S[i,(j+2)%L])] #3 spin row\n",
    "    val[0] = np.sum(S);\n",
    "    return val"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Integrated MC + RG simulation function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Energy(S):\n",
    "    '''Brute force Find energy of spin configuration S for sanity check'''\n",
    "    L = S.shape[0];\n",
    "    E = 0;\n",
    "    for i in np.arange(L):\n",
    "        for j in np.arange(L):\n",
    "            E += K*S[i,j]*(S[i,(j+1)%L] + S[(i+1)%L,j])\n",
    "    E += h*np.sum(S)\n",
    "    return E"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RunMCRG(K,h):\n",
    "    '''Run MCRG simulation to find y_t exponent, keeping Nc coupling terms'''\n",
    "    print('running MCRG for linear size',L,'lattice.')\n",
    "    print('Setting K =', K, \" and h = \",h)\n",
    "    \n",
    "    #measurement accumulators for y_t\n",
    "    evenK = np.zeros(7,dtype = float)\n",
    "    evenK_1 = np.zeros(7,dtype = float)\n",
    "    mix_11 = np.zeros((7,7),dtype=float)\n",
    "    mix_01 = np.zeros((7,7),dtype=float)\n",
    "    #measurement accumulators for y_h\n",
    "    oddK = np.zeros(4,dtype = float)\n",
    "    oddK_1 = np.zeros(4,dtype = float)\n",
    "    mix_11_odd = np.zeros((4,4),dtype=float)\n",
    "    mix_01_odd = np.zeros((4,4),dtype=float)\n",
    "    \n",
    "    #Run simulation\n",
    "    k = 0\n",
    "    for n in np.arange(nmeas+nwarm):\n",
    "        # Every MC n-loop, go through every [i,j] pair and propose \n",
    "        # flipping s[i,j]. \n",
    "        # Result: A S-field config drawn with probability propto Boltzmann weight\n",
    "        for j in np.arange(L):\n",
    "            for i in np.arange(L):\n",
    "                # calculate alpha based on energy difference\n",
    "                d_hterm = -2*h*S[i,j];\n",
    "                d_tterm = -2*K*S[i,j]*(S[(i+1) % L,j]+S[(i-1) % L,j]+S[i,(j+1)%L]+S[i,(j-1)%L])\n",
    "                alpha = np.exp(d_tterm);\n",
    "                #Probability Min(1,alpha) of accepting flip \n",
    "                #If accept, update sigma field\n",
    "                #If reject, do nothing, look at next entry of S\n",
    "                if beta_array[n,i,j] < alpha:\n",
    "                    S[i,j] = -S[i,j]; \n",
    "                    \n",
    "        #sanity check\n",
    "        if n % interval == 0:\n",
    "            energy[k] = Energy(S)\n",
    "            k +=1\n",
    "        \n",
    "        # take measurements every (interval) steps if finished warmup\n",
    "        if n % interval == 0 and n >= nwarm:\n",
    "            S1 = RGTransform(S,b);\n",
    "            evenK += AllEvenCoupling(S)\n",
    "            evenK_1 += AllEvenCoupling(S1)\n",
    "            oddK += AllOddCoupling(S)\n",
    "            oddK_1 += AllOddCoupling(S1)\n",
    "            #A*B = C, B is unknown, A is symmetric\n",
    "            #Problem: C is not symmetric!?\n",
    "            mix_11 += np.outer(AllEvenCoupling(S1),AllEvenCoupling(S1))\n",
    "            mix_01 += np.outer(AllEvenCoupling(S1),AllEvenCoupling(S))\n",
    "            mix_11_odd += np.outer(AllOddCoupling(S1),AllOddCoupling(S1))\n",
    "            mix_01_odd += np.outer(AllOddCoupling(S1),AllOddCoupling(S))\n",
    "    \n",
    "    #Results\n",
    "    evenK /= ndata; evenK_1 /= ndata; mix_11 /= ndata; mix_01 /= ndata;\n",
    "    oddK /= ndata; oddK_1 /= ndata; mix_11_odd /= ndata; mix_01_odd /= ndata;\n",
    "    print('evenK = ',evenK)\n",
    "    print('evenK_1 = ', evenK_1)\n",
    "    print('mix_11 = ', mix_11)\n",
    "    print('subtract ',np.outer(evenK_1,evenK_1))\n",
    "    print('mix_01 = ', mix_01)\n",
    "    print('subtract ',np.outer(evenK_1,evenK))\n",
    "    MatA = mix_11-np.outer(evenK_1,evenK_1)\n",
    "    MatC = mix_01-np.outer(evenK_1,evenK)\n",
    "    #print('MatA (lhs) = ',MatA)\n",
    "    #print('MatC (rhs) = ',MatC)\n",
    "    LinRGMat = la.solve(MatA[0:Nc_even,0:Nc_even],MatC[0:Nc_even,0:Nc_even])\n",
    "    print('linearized RG transformation = ',LinRGMat)\n",
    "    lmbd = la.eigvals(LinRGMat);\n",
    "    print('eigenvalues are',lmbd)\n",
    "    amplitude = np.absolute(lmbd)\n",
    "    print('eigenvalue amplitudes are',amplitude)\n",
    "    #Only the eigenvalue with maximum amplitude is important.\n",
    "    #This eigenvalue should generically be real\n",
    "    imax = np.argmax(amplitude)\n",
    "    y_t = np.log(lmbd[imax])/np.log(b)\n",
    "    print('exponent y_t = ',y_t,'\\n')\n",
    "    \n",
    "    print('oddK = ',oddK)\n",
    "    print('oddK_1 = ', oddK_1)\n",
    "    print('mix_11 = ', mix_11_odd)\n",
    "    print('subtract ', np.outer(oddK_1,oddK_1))\n",
    "    print('mix_01 = ', mix_01_odd)\n",
    "    print('subtract ', np.outer(oddK_1,oddK))\n",
    "    #TODO: cancellation bad, how to avoid?\n",
    "    MatA = mix_11_odd-np.outer(oddK_1,oddK_1)\n",
    "    MatC = mix_01_odd-np.outer(oddK_1,oddK)\n",
    "    #print('MatA (lhs) = ',MatA)\n",
    "    #print('MatC (rhs) = ',MatC)\n",
    "    LinRGMat = la.solve(MatA[0:Nc_odd,0:Nc_odd],MatC[0:Nc_odd,0:Nc_odd])\n",
    "    print('linearized RG transformation = ',LinRGMat)\n",
    "    lmbd = la.eigvals(LinRGMat);\n",
    "    print('eigenvalues are',lmbd)\n",
    "    amplitude = np.absolute(lmbd)\n",
    "    print('eigenvalue amplitudes are',amplitude)\n",
    "    #Only the eigenvalue with maximum amplitude is important.\n",
    "    #This eigenvalue should generically be real\n",
    "    imax = np.argmax(amplitude)\n",
    "    y_h = np.log(lmbd[imax])/np.log(b)\n",
    "    plt.figure();\n",
    "    plt.imshow(S);plt.colorbar();\n",
    "    print('exponent y_h = ',y_h)\n",
    "    \n",
    "    return y_t,y_h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test for specific set of input parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Known: y_t = 1, y_h = 1.875 exactly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATwAAAD8CAYAAADqmhgGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnX+0XWV55z/fhh8pqA0xgDEJElZjxY4E7J1gh1n+CCjR6QLbAQt2nOiClUWXtHasLVBnsKVjV2zXEu0aRr0DUWwtAbFKRqMxBhhXlwVyEQwECsSUgWvQGALaDhrIvc/8sfdl7XvuOfu8+z3v3ufcc56P66ycvff747nnHB/e532fHzIzHMdxRoFf6LcAjuM4TeEKz3GckcEVnuM4I4MrPMdxRgZXeI7jjAyu8BzHGRlc4TmOUxuSNknaL+nBDs8l6a8l7ZG0S9LrC8/WS3osf61PIY8rPMdx6uRzwLqS528HVuWvDcCnACQtBj4CnAmsAT4i6bhehelJ4UlaJ+mRXDtf2aswjuMMF2b2beBgSZPzgc9bxl3AIklLgXOB7WZ20MyeAbZTrjiDOCK2o6QFwHXAW4FJYKekLWb2UKc+R+loW8ixALz6tOc6jv3ormNmXXdq29qurE+xbdncZZTNFypHr+OlIlSuYrvQ76VsjLJ2ZYR+Z6G/iSqffYz8KX5/ZZ93p2ePP/kCBw5OKWiCDpz7lmPt6YNTQW3v3XVoN/Dzwq1xMxuvMN0y4MnC9WR+r9P9nohWeGTLzD1mthdA0mYybd1R4S3kWM7U2QBs23Z/x4HPfeXps647tW1tV9an2LZs7jLK5guVo9fxUhEqV7Fd6PdSNkZZuzJCv7PQ30SVzz5G/hS/v7LPu9OzNec+Sa88fXCKe7adFNR2wdLHfm5mYz1M1045W8n9nujFpA3SwJI2SJqQNPECh3qYznGcJjBgOvB/CZgEVhSulwP7Su73RC8rvCANnC9vxwHGVi+0mNVap3bb9sWt1EKZ81/RwPmK/aqM0elZ6aplX8lKYl/nFUHZ51jWrxN1rFZjvt/Y30TsZxU6RsxnGvrdFp89ak8HjV2GYbxgYSZtArYAl+cW4pnAT8zsKUnbgL8oHFS8Dbiq18l6UXi1aGDHcfpPotUbkm4C3gwskTRJdvJ6JICZfRrYCrwD2AM8B7wvf3ZQ0p8DO/OhrjGzssOPIHpReDuBVZJWAj8ALgLe3atAjuP0F8OYSpQ2zswu7vLcgPd3eLYJ2JREkJxohWdmhyVdDmwDFgCbzGx3Mskcx+kb072fDwwkvazwMLOtZEvSIB7ddUzQ/kWKvbnYvafi3DF7dq1U2bPreOJXw15l6J5STLuytlX2tkL3I0Ofhdyv2jZ0/JjfSOwecq8YMOUKz3GcUcFXeI7jjAQGvDCkpR/UZE2LsdUL7Z5tK7q2S+EOEkqsmRDiItBtvlhTOJRQd4pYF4oYU7Lb3L2OEeumEzpv3b/NXse723bwUzvYU6TFaauPtK9tXRLU9qTlP7y3R8fjRvEVnuM4szGYGs4Fnis8x3Fmk0VaDCeu8BzHaUFMtQ2kmv80qvCKbimx7hpNh5MVCXVZqbIH1KtMdY+RYo8zdPxYt5SYecuowx0kxX5nyPhrzo3LBFQkO7Rwhec4zgiQ+eG5wnMcZ0SY9hVe77z6tOeCctuFZu+INSs7zVuFGFeIKnKVUbfLR4wpGUqKrYCYLCK9yJViGyL2b4tp1yu+wnMcZ2QwxNSQlrtxhec4zhzcpE1MilPDKqd1dSaTjDVfQuerw4SLIbXZ19ovNkqi7qiU1OZjbKKFMtImABXP24KexxlEfIXnOM4sMsdjN2kdxxkR/NDCcZyRwExMma/waqVf7hop9sBSZddI4XLTZERCKCmK26TI/FK3jGXE7BE2uT/bynTCFZ6kdcAnyTKjX29mG1ueXwu8Jb88BjjBzBblz6aAB/JnT5jZeb3IMjAKz3GcwSA7tEijGiQtAK4D3kpW+GunpC1m9mL9ajP7L4X2vwecURjiZ2aWTNMP57rVcZxoZg4tQl4BrAH2mNleM3se2AycX9L+YuCm3v+K9vQteUArqRMkhhJz9N8qR6qkkHVGUNRhHqXw/E+xvVC3HEVSuM6kqIHSafwUyQMAptL54S0DnixcT5LVn52DpFcBK4HbC7cXSpoADgMbzewrvQjjJq3jOLOoGGmxJFdIM4yb2Xjhup3m7JRe9CLgVrNZVcBPMrN9kk4Bbpf0gJl9P1S4VlzhOY4zh+nwU9oDXVK8TwLFug7LgX0d2l5ES41aM9uX/7tX0p1k+3vRCs/38BzHmUWWPOAXgl4B7ARWSVop6SgypbaltZGkXwGOA/6xcO84SUfn75cAZwEPtfatQt+ypaQqotLpWYw7QpW5ymgyYWkdbjUpknfWWfc2ti5t3SF5dbjVFAnplyq07IVEoWVmdljS5cA2MreUTWa2W9I1wISZzSi/i4HNNruq2KnAZyRNky3ONhZPd2Nwk9ZxnFmYkdTx2My2Altb7l3dcv2nbfp9B3hdMkEIMGklbZK0X9KDhXuLJW2X9Fj+73EphXIcp5+I6cDXfCNkhfc54H8Any/cuxLYYWYbJV2ZX19RZeJYUy91QsoyEyjU7Knyt6RO9thkvYsq5laou0bd2WNiqCPSomz8TsT8nalqWgxraFnXv8rMvg0cbLl9PnBj/v5G4J2J5XIcp48kPLQYKGL38E40s6cAzOwpSSd0aihpA7AB4KRlvmXoOIOOIU8AGkvuhDgOMLZ64YsnMFVMgzpPXOs2gcrmS3F6Gdov9vQvdK5QOVKYo3WXUUwxRtNmcUpvgKxM43AuTmL/qh9JWpqv7pYC+1MK5ThOPxneQtyxRvgWYH3+fj1wWxpxHMfpN0YWaRHymm90XeFJugl4M1nM3CTwEWAjcIukS4AngAvrFNJxnGYZ1hVeV4VnZhd3eHR21cmK2VJSuHLUHcXQz4SUZX36VZu3SnRCp7ljCy+ldp1JkSi0Cqlr59aZDNRM83L1FsJw7kw6jhNNdmjhVcscxxkJvKZFX0kdAF53jdMyUpuBsW1DI0piE1emTowZOldZ29jfTopaKSloqqZFdmgxont4juOMHvMxiiIEV3iO48zCIy0cxxkpAgv0zDv6pvBS7J3F7iml2MtJsZ+SIgtK6r2n1n5NFsipm7oLLzX5t3X63SZJAGrwwrQrPMdxRoDMpHWF5zjOiDCykRZ1kSKvfxWv/V7blclRZcwYj/5YV44m3SRiIz76Ze7WEb1Sd23lxiItGF63lOFctzqO0wNKmjxA0jpJj0jak2dIb33+Xkk/lnR//rq08Gx9XkriMUnrW/tWxU1ax3HmkKpehaQFwHXAW8lq1O6UtKVN9bGbzezylr6LyZKVjJEtPO/N+z4TK89AntLGnrCWjd9prhRRDK1yhNZwSBGwHnrKPKhJM4uE/iZiT+fLiB0/9FndURIzcyepaWHwwnSyWNo1wB4z2wsgaTNZiYiQcovnAtvN7GDedzuwDrgpVhg3aR3HmcWM43HIiyxt3EThtaFluGXAk4XryfxeK/9R0i5Jt0paUbFvMG7SOo4zhwom7QEzGyt53m4ga7n+38BNZnZI0mVkhcHWBvathK/wHMeZxcwpbeAKrxuTwIrC9XJg36z5zJ42s0P55f8Cfi20b1Xm3QovxV5OiuSdZQVy+hmREBMlUUcho5h2oa4isft0ZXPFtk1dHKrX33eKSAsgpePxTmCVpJXAD4CLgHcXG8zUx8kvzwMezt9vA/5C0nH59duAq3oRZt4pPMdx6sVMHE6k8MzssKTLyZTXAmCTme2WdA0wYWZbgN+XdB5wmKwG9nvzvgcl/TmZ0gS4ZuYAIxZXeI7jzCGl47GZbQW2tty7uvD+Kjqs3MxsE7AplSwDqfBCTcRQs7UuuULkqKMOROog+BSuLaFjpPjO6vg8yvqUfVap626EytXp70zilsLwRloMpMJzHKe/uMJzHGck8ASgjuOMFKlCywaNgcyW0kqMe0WK0LIU7g9VSBE2F9onRShV6N5WmRxl1JkZJ9V3mzqUL9Y1KXUC0MOeANRxnFFhWE3armpc0gpJd0h6WNJuSR/I7y+WtD1P27K94BzoOM48pmIs7bwiZIV3GPhDM/uupJeSpWjZTuYcuMPMNuY5rq4Erigb6NWnPce2bdkyvcqSP4X5VaQO0yZkrtb5UtRhrSPqIHVUSiwp6uN2elZHTYvQfimiXur+3do8VGYhdF3hmdlTZvbd/P2/kIV9LCNL8XJj3uxG4J11Cek4TrNMo6DXfKPSHp6kk4EzgLuBE2fi38zsKUkndOizAdgAcNIy3zJ0nEHHbHj38II1kKSXAF8C/sDMfiqFfSBmNg6MA4ytXthTahfHcZpATI3yKa2kI8mU3RfM7O/z2z+ayXIgaSmwvxdBUmcCLmsbG4KWwtUidUhX0/V9Y+ZKsecYOkYdLjYxY6TqFzNeCkZ2D0/ZUu4G4GEz+3jh0RZgpqjGeuC29OI5jtM0ifPhDRQhK7yzgPcAD0ia+c/KnwAbgVskXQI8AVxYj4iO4zSKZft4w0hXhWdm/0D7VMsAZ1eZ7NFdx7y4fK+j4EmKgi0xz2LdRlJkhWklxo0ktl2TBXJiE66mztoSW8+2TKaYLDl1Mx9PYEPwY1PHcWZho35o4TjOaDGyJm1dxC7Xm6y/EGu+pD7ZLBsjhXlUR+LNGKokGA0dI4Y6EsvWudVQbJsiASiM8Cmt4zijhVmm8EJeIUhaJ+kRSXvyMNTW5x+U9FBel3aHpFcVnk1Juj9/ben1b3OT1nGcOaRyOZG0ALgOeCtZ2cWdkraY2UOFZvcBY2b2nKTfBf4S+O382c/MLFmgsK/wHMeZg1nYK4A1wB4z22tmzwObyeLwC3PZHWY2Y4vfRVZ/thbmRQLQ1PtSqfeD6nAfSLE/FrtXGUPdWVtSuDGl+E2kcD1JsVdZJ4aYDj+lXSJponA9noeTzrAMeLJwPQmcWTLeJcDXC9cL8/EPAxvN7CuhgrXDTVrHceZQ4ZD2gJmNlTxvZxu3HV7SfwLGgDcVbp9kZvsknQLcLukBM/t+uHizcYXnOM5sLOkp7SSwonC9HNjX2kjSOcCHgTeZ2aEXRTHbl/+7V9KdZNma5ofCKyYAbSV1Tc8mzaNYM6TupJlNjd1u/BS1hOs2M3ttV9Yvxe+v6Zoqs0jnh7cTWCVpJfAD4CLg3cUGks4APgOsM7P9hfvHAc+Z2SFJS8jCXP+yF2F8hec4zhxSrfDM7LCky4FtwAJgk5ntlnQNMGFmW4C/Al4CfDFPO/eEmZ0HnAp8RtI02QHrxpbT3cq4wnMcZxYGTE+nczw2s63A1pZ7Vxfen9Oh33eA1yUTBFd4juO0YsCQRloMZGhZipqy/XQpCZU/dJ+njiI+MWPUvf/WOkaKIj5173uFJm0te5ZirzJlXVrwWFrHcUYJV3iO44wG4XGy842BibRIndc/NuNFP+sLhCYYjUlSGmti1R2tkaKWSUwERYpaJrGkzphTC77CcxxnJDCwhKe0g4QrPMdx2uAKr1ZiTutaiTEDY82jFBEfsaTw6I8dv0iK+h8pSk72y+RMNXfKpKWpEoC6Ses4zujgCs9xnJHAHY8dxxkl3PE4AaF1aVPVDO00Rmxd2tC5BpG6C/CkSIxZx35Y6sicVlJECKX+LpLgp7SO44wKGtIVXtc8zpIWSrpH0vck7Zb0Z/n9lZLulvSYpJslHVW/uI7j1I5VeM0zQlZ4h4C1Zvavko4E/kHS14EPAtea2WZJnybLRf+p0ImbrkeROnqg7mSSsWZOaOLN0Gep3XTqGD/ULO5nNEgovcqYJnmAhvbQousKzzL+Nb88Mn8ZsBa4Nb9/I/DOWiR0HKd5hnSFF1SaSNICSfcD+4HtZDnlnzWzw3mTSbLqRO36bpA0IWniBQ61a+I4zqAxHfiaZwQpPDObyovhLierM3lqu2Yd+o6b2ZiZjR3J0fGSOo7TDDN+eCGveUalU1ozezavHPQGYJGkI/JVXttKRK0Ui/iUhY/VQZ3FbeoO76rSLjSkK5QUIX+d+tQxfrf5Yki9bxebcLWsT+rQspSntJLWAZ8kq2lxvZltbHl+NPB54NeAp4HfNrPH82dXkZ0PTAG/b2bbepEl5JT2eEmL8ve/CJwDPAzcAVyQN1sP3NaLII7jDBCJ9vAkLQCuA94OvBa4WNJrW5pdAjxjZr8MXAt8LO/7WrIqZ78KrAP+Zz5eNCEm7VLgDkm7yEqubTezrwJXAB+UtAd4OXBDL4I4jjOUrAH2mNleM3se2Ayc39LmfLKDT8gOQs9WVr7sfGCzmR0ys38G9uTjRdPVpDWzXWTFb1vv7606eWikRSsxdQPqqPUQSh1mcifqqAPRa7tWYrKvtPYrux/jihKbmSXFFkLdrj4pqGDSLpE0UbgeN7PxwvUy4MnC9SRwZssYL7bJyzr+hGwRtQy4q6Vv28PRUDzSwnGc2RhVQssOmNlYyfN2A7Wq005tQvpWIuiU1nGcESOdH94ksKJw3e6A88U2ko4Afgk4GNi3EgNTprGMmOiBOiI5mkyu2alPt379SppZR3RMr+O1jlnH5x2zhVBl7lA5UpPwlHYnsErSSuAHZIcQ725ps4Xs4PMfyQ5Cbzczk7QF+DtJHwdeCawC7ulFGDdpHceZSyKFl+/JXQ5sI3NL2WRmuyVdA0yY2RayA8+/yQ9AD5IpRfJ2twAPAYeB95vZVC/yuMJzHGcuCf3wzGwrsLXl3tWF9z8HLuzQ96PAR1PJ4grPcZxZyIY3PdTAKLxQj/vQMWLbxey91LHHVkeyyhhSZ5apwy0llNjiTUXqrqcc+r13kiNNthQ8AajjOKODr/AcxxkdXOGlpW7TIDagP9TcSFEbNkUd1lgTKFSOmM+jW78YUpvxdWwLxLotpXBtSZo8wPfwHMcZKVzhOY4zKmgeJvcMwUPLHMcZGRpd4ZUlAG0ybKbuPbbQucuI3cuJqb/bz+SdKQrk1J2INLRfTLuyfjH1mZO5pbhJ6zjOSOCHFo7jjBSu8NJSxSypsx5FFVJEHcRElKQw4VJEFlT5vFMkxozJcBOboDN19EaV8UPN4rI+yd1sXOE5jjMKiOE9pXWF5zjObHwPr37qrkuQus5EbBnFmLli5UidxKCOxJv9qqcRe/IdKlfdiTzrLtPoJq3jOKODKzzHcUYFN2kdxxkdXOHVS2rP/FgXh059qjyLadfKICb5TOEm0XT0St1uKSHztlL3vmvPWDOntJIWAzcDJwOPA+8ys2da2pwOfAp4GTAFfNTMbs6ffQ54E/CTvPl7zaz0Q/JYWsdx5pKuTGMZVwI7zGwVsCO/buU54D+b2a8C64BPSFpUeP5HZnZ6/ur6X4RghSdpgaT7JH01v14p6W5Jj0m6WdJRoWM5jjPYzNS16PbqkfOBG/P3NwLvbG1gZo+a2WP5+33AfuD42AmrmLQfAB4mW1oCfAy41sw2S/o0cAnZ0rMjj+465sXlfOvyvO4A8xhSBIrX4ZZSRowpmcK8S137InaMFIlfUwT+VyHF7yU54cpsiaSJwvW4mY0H9j3RzJ4CMLOnJJ1Q1ljSGuAo4PuF2x+VdDX5CtHMDpWNEbTCk7Qc+A/A9fm1gLXArXmTttrZcZx5SKg5mynFA2Y2VnjNUnaSviXpwTav86uIJGkp8DfA+8xsZofxKuA1wL8FFgNXdBsndIX3CeCPgZfm1y8HnjWzw/n1JLCsg6AbgA0ACzkmcDrHcfqFSOeWYmbndJxH+pGkpfnqbimZudqu3cuArwH/1czuKoz9VP72kKTPAh/qJk/XFZ6k3wD2m9m9xdttmrb9iMxsfEb7H8nR3aZzHGcAaGgPbwuwPn+/HrhtjhzZ2cCXgc+b2Rdbni3N/xWZhflgtwlDVnhnAedJegewkGwP7xPAIklH5Ku85cC+bgOVJQAtUrdbSih17zfVkSQyhtgMI1XGDJmrlX7u14Y+CyXF7y/1915KM354G4FbJF0CPAFcCCBpDLjMzC4F3gW8EXi5pPfm/WbcT74g6XiyBdj9wGXdJuyq8MzsKjJbGUlvBj5kZr8j6YvABcBmOmhnx3HmKQ0oPDN7Gji7zf0J4NL8/d8Cf9uh/9qqc/bih3cF8EFJe8j29G7oYSzHcQaFQHN2PoafVYq0MLM7gTvz93uBNbET1+GJnsLMbJJQkyXWDExRszaFu0YoKdw/UmyH1F1Ttg5XpZkxvaZFOQMTWuY4zuDgCUAdxxkZ5qO5GkLfFF4dCRJjygZWMdNSmy9l0SaxAf11ll+svY5CF7k6yRHSp1u71GPWccLaWBRGmjjZgcRXeI7jzMUVnuM4o0DKSItBwxWe4zhz0PRwaryBUXipM4ykzr5SNkaoTFXGaLJoTR3uMcW2KfY4+xXt0E2OmDHriGxJWsTH9/Acxxkl3KR1HGd0cIWXlpjleitlrieh41UxlWLkiE1WWUYKk7zuhJcp6rWGJiAoo+6tkjqje6oklk0daeErPMdxRgdXeI7jjAQNVS3rB67wHMeZhfvhNUCM60KKhJRl1F3cJsVeWd01dut2G4mRq466rqnCA2OIcSWqHRtOjTcwCs9xnMHBV3iO44wG7nichmJd2jJSePDHEhOdkMrkSV13IzVVMoD0yxxLUSO4iukbut2S2izu1CdJpAXNHFpIWgzcDJwMPA68y8yeadNuCnggv3zCzM7L768kKzGxGPgu8B4ze75szl5SvDuOM6RoOuzVI1cCO8xsFXkh7Q7tfmZmp+ev8wr3PwZcm/d/Brik24Su8BzHmY2RHVqEvHrjfODG/P2NZKUWg8hLM64Fbq3Sv1GTtlimMZa6T7FCzbRGS+a1kHq+0CSiseZiTGLWKsScvqb6DFNsbcR8Bp3m6kOkxRJJE4XrcTMbD+x74kwx7bwY9wkd2i3M5zgMbDSzr5AVDns2LxMLMAks6zahH1o4jjOXcIV3wMzGOj2U9C3gFW0efbiCNCeZ2T5JpwC3S3oA+Gmbdl2ldoXnOM4sUjoem9k5HeeRfiRpab66Wwrs7zDGvvzfvZLuBM4AvgQsknREvspbDuzrJo/v4TmOMxszNB326pEtwPr8/XrgttYGko6TdHT+fglwFvCQmRlwB3BBWf9WBqaIT4q9kNT7TVXGL5La/SG2IEys20infrGuFjH7omXjhcpbB6nq2Yb0i/neU7mlNOSHtxG4RdIlwBPAhQCSxoDLzOxS4FTgM5KmyRZoG83sobz/FcBmSf8duA+4oduEbtI6jjOHJiItzOxp4Ow29yeAS/P33wFe16H/XmBNlTmDFJ6kx4F/AaaAw2Y2Fuo06DjOPMMAr2nBW8zsQOF6xmlwo6Qr8+srQgdL4Zlfh3d/3W4SneYq61eHO0hq8zy0Xx3JRkPHjzHjq4wf0y5WjtoZTn3X06FFtNOg4ziDjSzsNd8IVXgGfFPSvZI25PdmOQ0CbZ0GJW2QNCFp4sdPT/UuseM4tdPQKW3jhJq0Z+WOfycA2yX9U+gEudf1OMDY6oXz7xNynFFj1LOlFBz/9kv6MtnJSJDTYJEU2VKKxO6TDMoeWxmp3U3qzmYS6zpTJlPq7zN2T7PJescpZOyVzPF4ODVeV5NW0rGSXjrzHngb8CABToOO48xTpgNf84yQFd6JwJez5AQcAfydmX1D0k7aOA06jjP/GdYVXleFlzv3rW5zv63TYBll2VJiPO7rMKM6zRU7fhmpXTJix6j7c0wRQRFKP6M1UkfO9C2p6qjv4TmOM0rMzxPYEFzhOY4zl1E1aR3HGTGGuBC3rEFNPrZ6od2zbQVQf+hNipCr0DFj91pS71+VjVllvNThUikYFDeaOlyQOhHz+7jbdvBTO6he5n3ZS5bZmat/N6jtt77z3+4tSwA6aPgKz3GcuQynResKz3GcuWh6OG3agVF4Kep2pk6M2fqs7kiO0DFSuDiEylEmU90RH6HFhMqI6ZdiOyE2SqLuustBGPPSqTiEgVF4juMMBsJG1/HYcZwRxBVeWqqc8KXwuE9xepmC1NEDKeRPccqcwlxsOsFBp2dNRoaU9Ysx/9PVtKhf4YVkTZf0FuDawq3XABeZ2VckfQ54E/CT/Nl7zaxUsXjVMsdxZjOzh1d/8oCZrOmrgB359WxRzO4ws9PN7HRgLfAc8M1Ckz+aed5N2YErPMdx2qDp6aBXj1TNmn4B8HUzi17GusJzHKcFy0zakFdvBGVNL3ARcFPLvY9K2iXp2pn6tWU0uodXTACaosZplT2TTi4UddSv7TRvt36d/rbYiIZ+7YGVzV13FEMK95gUNXzL2pXR18I9MxhVlNkSSROF6/E8yzkAkr4FvKJNvw9XESlPMvw6YFvh9lXAD4GjyLKqXwFcUzaOn9I6jjOXcGv1QFlomZmd0+mZpCpZ098FfNnMXiiM/VT+9pCkzwIf6iasm7SO48xBZkGvHqmSNf1iWszZXEmiLDvxO8kysZfS6AqvmAC0ivmS2mxI7eYRa1rHkiLCIXT8svsxbh6DEsVQNn6VZ51IkZSijNpN32b88DbSJmu6pDHgMjO7NL8+GVgB/J+W/l+QdDxZGY77gcu6TegmreM4szGDqfpjyzplTTezCeDSwvXjwLI27dZWndMVnuM4c/FIC8dxRgZXeL1TVpc2dYGcVlLXlK07JC1FGFSROvbOYuZLlXgzZozQ8QbCNYS439+j9nTvExvgNS0cxxkNDGw480O5wnMcZzZGI4cW/aBvNS2qEBMlUTfzynO+C02bu73OVUZslpy6f1d1R2TM9EtR0+KXjjrR/t2JFwW1/cbkX3tNC8dx5jlDemgRFGkhaZGkWyX9k6SHJf26pMWStkt6LP/3uLqFdRynCRpLHtA4oSu8TwLfMLMLJB0FHAP8CVkuq42SriTLZXVF2SBlyQPKSOHp3mm8WI/41KeXVcaMkT+FmTaoyQNi63+UzR3aLnVi2Sa3CTpiwJAW8em6wpP0MuCNwA0AZva8mT1L9VxWjuPMF0Z4hXcK8GPgs5JWA/cCH6All5WktrmsJG0ANgAs5JgkQjuOUyfNhJb1g5A9vCOA1wOfMrMzgP9Hm1TMnTCzcTMbM7OxI+man89xnH4+6FuqAAAFjUlEQVRjYDYd9JpvhKzwJoFJM7s7v76VTOFVyWU1hzoyXsTOHdMuNmFk6gSjKSIXUkQ4dJMrZoxO49WRTSf2e4mJKImpxRvaL1kRnyGNtOi6wjOzHwJPSvqV/NbZwENUy2XlOM58YoT38AB+jyz31FHAXuB9ZMpyTi4rx3HmOWZDe0obpPDy8mftvKnn5LJKQR01HGKSAsTKMYjuFLHuIDFzNU0Kd5466uqGPksx/ky/JMkDYF6u3kLwSAvHcVowbGqq30LUgis8x3Fm4+mhHMcZKeahy0kIfVN4sa4cVcYMGb+KHKFjx7oWhBbPKaPukLEiMftedWePSfFbSdEv1u0q1v0mpVuKAdbACk/ShcCfAqcCa/JaFu3arSMLb10AXG9mG/P7K4HNwGLgu8B7zOz5sjm9TKPjOLOxPAFoyKs3HgR+C/h2pwaSFgDXAW8HXgtcLOm1+eOPAdea2SrgGeCSbhO6wnMcZw42NRX06mkOs4fN7JEuzdYAe8xsb7562wycn9eiXUsWCAGB8fyNJgCV9GPg/wJLgAONTdyeQZABXI5WXI7ZVJXjVWZ2fC8TSvpGPm8IC4GfF67HzWy84nx3Ah9qZ9JKugBYV6hR+x7gTDJT+C4z++X8/grg62b2b8rmanQPb+aLkDTR7yypgyCDy+FyDKIcZrYu1ViSvgW8os2jD5tZSHRWu+zNVnK/FD+ldRynNszsnB6HmASKdSGWA/vIVr2LJB1hZocL90vxPTzHcQaZncAqSSvz0NaLgC2W7cXdAVyQtwuK5++Xwqtk49fEIMgALkcrLsdsBkWO5Ej6TUmTwK8DX5O0Lb//SklbAfLV2+XANuBh4BYz250PcQXwQUl7gJeTJykunbPJQwvHcZx+4iat4zgjgys8x3FGhkYVnqR1kh6RtCevdNbUvJsk7Zf0YOFe42UmJa2QdEde6nK3pA/0QxZJCyXdI+l7uRx/lt9fKenuXI6b803i2pG0QNJ9kr7aLzkkPS7pAUn3S5rI7/XjN+IlUWukMYXXJUSkbj4HtPoWXUlWZnIVsIMKdTp64DDwh2Z2KvAG4P35Z9C0LIeAtWa2GjgdWCfpDUSE6iTiA2Qb0jP0S463mNnpBb+3fvxGZkqivgZYTfa59EOO4cTMGnmRncRsK1xfBVzV4PwnAw8Wrh8BlubvlwKPNCVLQYbbgLf2UxayGsPfJfNePwAc0e77qnH+5WT/J14LfJXMobQfcjwOLGm51+j3ArwM+Gfyw8R+yTHMryZN2mXAk4Xryfxev5hVZhJoW2ayLiSdDJwB3N0PWXIz8n6y4kvbge8Dz1rmBgDNfT+fAP4YmIlEf3mf5DDgm5LuzUuLQvPfS7Ek6n2Srpd0bB/kGFqaVHhRoSDDiKSXAF8C/sDMftoPGcxsysxOJ1thrSFL0TOnWZ0ySPoNYL+Z3Vu83bQcOWeZ2evJtlzeL+mNDczZSk8lUZ3uNKnwOoWI9Isf5eUliSkzGYukI8mU3RfM7O/7KQuAmT0L3Em2p7hI0ky4YRPfz1nAeZIeJ8uCsZZsxde0HJjZvvzf/cCXyf4j0PT30q4k6uv7IMfQ0qTCaxsi0uD8rTReZjJPaXMD8LCZfbxfskg6XtKi/P0vAueQbY5XDtXpBTO7ysyWm9nJZL+H283sd5qWQ9Kxkl468x54G1mutka/F/OSqPXT5IYh8A7gUbL9og83OO9NwFPAC2T/Fb2EbK9oB/BY/u/iBuT492Tm2S7g/vz1jqZlAU4D7svleBC4Or9/CnAPsAf4InB0g9/Rm4Gv9kOOfL7v5a/dM7/NPv1GTgcm8u/mK8Bx/ZBjWF8eWuY4zsjgkRaO44wMrvAcxxkZXOE5jjMyuMJzHGdkcIXnOM7I4ArPcZyRwRWe4zgjw/8HUyHF1IRAZ3sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Input Parameters\n",
    "L = 64; #linear dimension of square lattice\n",
    "Kc = np.arccosh(3)/4; # Critical temperature Kc assumed to be known\n",
    "K = Kc; h = 0; #Start on Critical manifold\n",
    "\n",
    "nwarm = 400; #number of warm up Monte Carlo sweeps.\n",
    "nmeas = 400; #number of measurement Monte Carlo sweeps.\n",
    "interval = 10; # Take data for correlation functions every (interval) MC-steps/site\n",
    "# RG analysis setting\n",
    "Nc_even = 1; #number of coupling constants included in calculating T\n",
    "Nc_odd = 1;\n",
    "b = 2; #scaling factor\n",
    "\n",
    "# Derived constants\n",
    "ndata = nmeas//interval\n",
    "Ns = L*L; #total number of grid points\n",
    "energy = np.zeros((nmeas+nwarm)//interval,dtype=float)\n",
    "\n",
    "#Initialize 2d spin field\n",
    "S = np.random.choice([-1,1],(L,L))\n",
    "plt.imshow(S);plt.colorbar()\n",
    "#Initialize an array of random beta in [0,1] to use in M-Hastings\n",
    "beta_array = np.random.random_sample((nwarm+nmeas,L,L))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "running MCRG for linear size 64 lattice.\n",
      "Setting K = 0.4406867935097715  and h =  0\n",
      "evenK =  [5888.5 5343.5 5027.9 9744.5 4638.8 2482.5 1996.3]\n",
      "evenK_1 =  [1463.5 1361.1 1278.8 2505.2 1205.2  602.7  524. ]\n",
      "mix_11 =  [[2149985.2 2001480.4 1882956.  3689812.  1776448.   887222.8  772544. ]\n",
      " [2001480.4 1864188.4 1754407.6 3438342.4 1655807.6  826182.8  720206.8]\n",
      " [1882956.  1754407.6 1652651.2 3238837.6 1560438.   777705.2  678089.6]\n",
      " [3689812.  3438342.4 3238837.6 6348600.8 3058908.  1524167.2 1329395.2]\n",
      " [1776448.  1655807.6 1560438.  3058908.  1474559.2  733999.6  640422. ]\n",
      " [ 887222.8  826182.8  777705.2 1524167.2  733999.6  366794.8  319370. ]\n",
      " [ 772544.   720206.8  678089.6 1329395.2  640422.   319370.   279204.8]]\n",
      "subtract  [[2141832.25]]\n",
      "mix_01 =  [[ 8634368.   7843471.6  7384754.8 14320483.2  6822541.2  3641943.6\n",
      "   2934964.4]\n",
      " [ 8033702.4  7299966.8  6874335.2 13333260.   6354066.4  3388756.4\n",
      "   2732592.8]\n",
      " [ 7553630.   6866490.8  6467980.  12547917.6  5981729.2  3186795.6\n",
      "   2571995.2]\n",
      " [14799339.6 13454086.8 12673749.2 24588677.2 11722708.8  6243774.8\n",
      "   5039979.6]\n",
      " [ 7122482.8  6476503.6  6101940.8 11840107.6  5645792.8  3005102.8\n",
      "   2426828. ]\n",
      " [ 3559280.   3234950.   3046453.6  5909340.8  2816212.   1501808.4\n",
      "   1211498.4]\n",
      " [ 3096542.8  2815773.2  2652357.6  5146578.4  2454078.   1306727.6\n",
      "   1055056. ]]\n",
      "subtract  [[8617819.75]]\n",
      "linearized RG transformation =  [[2.03]]\n",
      "eigenvalues are [2.03+0.j]\n",
      "eigenvalue amplitudes are [2.03]\n",
      "exponent y_t =  (1.0212845856618016+0j) \n",
      "\n",
      "oddK =  [2735.65 9066.   8685.8  4536.  ]\n",
      "oddK_1 =  [ 741.2 2436.  2344.4 1212.2]\n",
      "mix_11 =  [[ 554420.4 1825068.8 1756912.   907938.2]\n",
      " [1825068.8 6015683.2 5791768.  2992574.4]\n",
      " [1756912.  5791768.  5577859.2 2881453.6]\n",
      " [ 907938.2 2992574.4 2881453.6 1489042.4]]\n",
      "subtract  [[ 549377.44 1805563.2  1737669.28  898482.64]\n",
      " [1805563.2  5934096.   5710958.4  2952919.2 ]\n",
      " [1737669.28 5710958.4  5496211.36 2841881.68]\n",
      " [ 898482.64 2952919.2  2841881.68 1469428.84]]\n",
      "mix_01 =  [[ 2046538.5  6787743.6  6506208.   3396598.4]\n",
      " [ 6736866.8 22352156.8 21430579.2 11185773.6]\n",
      " [ 6485062.  21517816.  20631899.2 10768505.6]\n",
      " [ 3351501.6 11119778.4 10661840.8  5564738. ]]\n",
      "subtract  [[ 2027663.78  6719719.2   6437914.96  3362083.2 ]\n",
      " [ 6664043.4  22084776.   21158608.8  11049696.  ]\n",
      " [ 6413457.86 21254330.4  20362989.52 10634198.4 ]\n",
      " [ 3316154.93 10989805.2  10528926.76  5498539.2 ]]\n",
      "linearized RG transformation =  [[3.743]]\n",
      "eigenvalues are [3.743+0.j]\n",
      "eigenvalue amplitudes are [3.743]\n",
      "exponent y_h =  (1.904112555496899+0j)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATwAAAD8CAYAAADqmhgGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X20XXV95/H3p+EhxYeGENCYRIlrYsUuBeydaIdZPgSU6HSBMwMWbJ3ogpXFLGntWFugzsKWjl2xs0Z01lD1VqLYWgJilYyNXiMP4+qyQi6KgUCBGBmIoYYkoHaigdz7nT/2Pqxzzz3n3H3O2Xufs8/+vFh73bOff3efyze/32//HhQRmJnVwS8NOwFmZmVxwDOz2nDAM7PacMAzs9pwwDOz2nDAM7PacMAzs8JI2ixpv6T7O+yXpP8pabeknZJe27Rvg6RH0mVDHulxwDOzIn0OWN9l/9uANemyEfgkgKSlwIeB1wFrgQ9LOnHQxAwU8CStl/RQGp2vHDQxZjZeIuJbwKEuh5wPfD4S3wGWSFoOnAtsj4hDEfEUsJ3ugTOTY/o9UdIi4DrgLcBeYIekrRHxQKdzli1dFKeuOrbfWwLw8M4T+jrvFa85XOj1805Hvzqlv+j7Vl3zc2t9Vt32DUvr99xI16OPP8uBQzMa5Nrnvvl5cfDQTKZj79l5ZBfwi6ZNkxEx2cPtVgCPN63vTbd12j6QvgMeSTZzd0TsAZC0hSRadwx4p646lrunVg1wSzj3JWf0dd7U1L2FXj/vdPSrU/qLvm/VNT+31mfVbd+wtH7PjXStPffxdof35OChGe6eemmmYxctf+QXETExwO3aBefosn0ggxRpM0VgSRslTUuafvJgtn81zGx4ApjN+F8O9gLNuaCVwL4u2wcySA4vUwROs7eTABOnL67dSAVT+7LnCObkJHo4r9P9is6tVtm8HFLG593tmfb7nfXzvfd7ryyC4NkoLXOyFbg8LSG+DvhJRDwhaQr486YXFW8Frhr0ZoMEvEIisJkNX065NyTdCLwJWCZpL8mb12MBIuJTwDbg7cBu4DDw3nTfIUl/BuxIL3VNRHR7+ZHJIAFvB7BG0mrgR8BFwLsGTZCZDVcQzOQ0bFxEXLzA/gDe12HfZmBzLglJ9R3wIuKopMuBKWARsDkiduWWMjMbmtnB3w+MpEFyeETENpIs6Ugosl6jF3mkI4+6oryfRxH1V8OSV3o71ZlW7Xk0C2DGAc/M6sI5PDOrhQCeHdOpHyoX8IouKgyzWUeVi0F10e07yqMp0SgIwkVaM6uJgJnxjHcOeGY2V9LTYjw54JlZCzHTtiNV9VU64PXSTKLIJgPdrjfOTTmG1XWt325h/ara9zSo5KWFA56Z1UDSDs8Bz8xqYtY5vPJkLQb2UtTodGwvxbJ+ija9FAPHpaV+L7I+/36/9zo+00E5h2dmtRGImTGd7sYBz8zmcZG24joOfV5yMaeIlvrDUmZ6sxZN6zroaeP3fjgODnytQDwTiwa+ziiqTcAzs2yShscu0ppZTfilhZnVQoSYCefwSlNE3VDV6sfqosz6uH7rAavwt9NI49pz85k3dzbHHJ6k9cAnSEZG/0xEbGrZfy3w5nT1BOCUiFiS7psB7kv3PRYR5w2SlpEMeGY2PMlLi3xCg6RFwHXAW0gm/tohaWtEPDd/dUT8l6bjfxc4s+kSP4+I3P7lG898q5n1rfHSIsuSwVpgd0TsiYhngC3A+V2Ovxi4cfDfor1Sc3gP7zzhueJCFYoJlr88iot1bXpSppn82uGtAB5vWt9LMv/sPJJeBqwGbm/avFjSNHAU2BQRXxkkMS7SmtkcPfa0WJYGpIbJiJhsWm8XOTsNL3oRcEvEnFnAXxoR+yS9HLhd0n0R8YOsiWvlgGdm88xmf0t7ICImuuzfC6xqWl8J7Otw7EW0zFEbEfvSn3sk3UlSv9d3wHMdnpnNkQwe8EuZlgx2AGskrZZ0HElQ29p6kKRfBU4E/rFp24mSjk8/LwPOAh5oPbcXzuHVxKg0tSiiGUrWiZeq0CVt0HTk1bXs2Zy6lkXEUUmXA1MkzVI2R8QuSdcA0xHRCH4XA1si5kyXdhrwaUmzJJmzTc1vd/vhgGdmc0SQa8PjiNgGbGvZdnXL+p+0Oe/bwKtzSwgZirSSNkvaL+n+pm1LJW2X9Ej688Q8E2VmwyRmMy5VkyWH9zngfwGfb9p2JXBbRGySdGW6fkUvNx6VItY4yzrYaDf9fi9ljm6S9Rr93qvoQURHpTjdEOSbwxslC/5WEfEt4FDL5vOBG9LPNwDvyDldZjZEOb60GCn91uG9KCKeAIiIJySd0ulASRuBjQCLOaHP25lZWQJ5ANB+pY0QJwFeqKXPvYHJo2hQhWJxHsWVUe3oPopTUI7KVJLddPvORiG9yTSN4/k+s9/f6seSlqe5u+XA/jwTZWbDNL4TcfdbCN8KbEg/bwBuzSc5ZjZsQdLTIstSNQvm8CTdCLyJpM/cXuDDwCbgZkmXAI8BFxaZSDMr17jm8BYMeBFxcYddZ+ecFstBkXVARdTZ5d1kJWsah1lXlnUip2HV9UWokrm3LMazZtLM+pa8tPCsZWZWC57TonCj2MQhD72kvZ9iStbrl10czZqOvL/bbs1SRqHJR146Pbc85rRIXlrUtA7PzOqnir0osnDAM7M53NPCzGol4wQ9lTO0gDeudXaDKLLZQR7PtOxr9PM88hhEdFQMK00R8OysA56Z1UBSpHXAM7OaqG1Pizy94jWHmZoqLps+TsXkvEeTqdrv34th9dbo95qjPlLQODdLGc98q5kNQLkOHiBpvaSHJO1OR0hv3f8eSU9KujddLm3atyGdSuIRSRtaz+2Vi7RmNk9e81VIWgRcB7yFZI7aHZK2tpl97KaIuLzl3KUkg5VMkGQ870nPfarf9FQi4HUqDvQ72KOLetWT9xvi1r+Vft8CZ01X1kEBsiryu03e0ubWl3YtsDsi9gBI2kIyRUSW6RbPBbZHxKH03O3AeuDGfhPjIq2ZzdFoeJxlIRk2brpp2dhyuRXA403re9Ntrf6jpJ2SbpG0qsdzM6tEDs/MytVDkfZAREx02d/uQtGy/r+BGyPiiKTLSCYGW5fx3J44h2dmczTe0mbM4S1kL7CqaX0lsG/O/SIORsSRdPWvgF/Pem6vKpHD66eeJI/jrPqK/q7znt93VOqXc2x4vANYI2k18CPgIuBdzQc05sdJV88DHkw/TwF/LunEdP2twFWDJKYSAc/MyhMhjuYU8CLiqKTLSYLXImBzROySdA0wHRFbgd+TdB5wlGQO7Pek5x6S9GckQRPgmsYLjH454JnZPHk2PI6IbcC2lm1XN32+ig45t4jYDGzOKy2VCHjDyuZXYd7bcdapuFj172HUm1ONc0+LSgQ8MyuXA56Z1YIHADWzWsmra9moqUTAG1adTdXrivo1Kk0jijQqkysN+5rtRMBRDwBqZnUxrkXaBcO4pFWS7pD0oKRdkt6fbl8qaXs6bMv2psaBZlZhPfalrZQsObyjwB9ExHclvYBkiJbtJI0Db4uITekYV1cCVxSX1GK46cl8dXwGecwhksccwd2OK/N7iQoGsywWzOFFxBMR8d30889Iun2sIBni5Yb0sBuAdxSVSDMr1yzKtFRNT3V4kk4FzgTuAl7U6P8WEU9IOqXDORuBjQAvXeEqQ7NRFzG+dXiZI5Ck5wNfAn4/In4qZXsgETEJTAJMnL54oKFdzKwMYqbOb2klHUsS7L4QEX+Xbv5xY5QDScuB/UUlMqtuzSmyjppcR8OsK+r2nY3rd1OF36u2dXhKsnLXAw9GxMeadm0FGpNqbABuzT95Zla2nMfDGylZcnhnAe8G7pPU+Kfpj4FNwM2SLgEeAy4sJolmVqpI6vHG0YIBLyL+gfZDLQOcnW9yBtOtqDAq87z2U+wuIh3NhlnEGta9hzmHcRWaQlXxDWwWfm1qZnNE3V9amFm91LZIa8XpVqzKo+V/t+uNypvYYckrHf38bkU8g0Y6Ho6DuVyvtm9pzaxeIpKAl2XJQtJ6SQ9J2p12Q23d/wFJD6Tz0t4m6WVN+2Yk3ZsuWwf93ZzDM7N58mpyImkRcB3wFpJpF3dI2hoRDzQd9j1gIiIOS/rPwF8Av5Xu+3lE5FbccQ7PzOaJyLZksBbYHRF7IuIZYAtJP/yme8UdEXE4Xf0OyfyzhXAOrwejUv/ULGvdXB2bnpRhVOpCG+trzz3MoAIxm/0t7TJJ003rk2l30oYVwONN63uB13W53iXA15rWF6fXPwpsioivZE1YOw54ZjZPDy9pD0TERJf97crGbS8v6XeACeCNTZtfGhH7JL0cuF3SfRHxg+zJm8sBz8zmilzf0u4FVjWtrwT2tR4k6RzgQ8AbI+LIc0mJ2Jf+3CPpTpLRmhzwytZPr4i8jHMR0UZEfu3wdgBrJK0GfgRcBLyr+QBJZwKfBtZHxP6m7ScChyPiiKRlJN1c/2KQxDjgmdk8eeXwIuKopMuBKWARsDkidkm6BpiOiK3AfweeD3wxHXbusYg4DzgN+LSkWZIXrJta3u72zAHPzOYIYHY2v4bHEbEN2Nay7eqmz+d0OO/bwKtzSwgOeGbWKoAx7WkxkgGvCqNJNOtl5I0i6vSsHsqdxKe0W5VqJAOemQ2ZA56Z1UP2frJVM5IBb1SKsHkUP4sowo7i6COtRj2NVas2KZ1zeGZWCwGR41vaUeKAZ2ZtOODlqgpFiqxzTozzm9d+534Yxe+z2ainb+hcpDWz2nDAM7NacMNjM6sTNzzOWRXrULIOrjnMOr1Rbw7STRXqdWvDb2nNrC40pjm8BcdxlrRY0t2Svi9pl6Q/TbevlnSXpEck3STpuOKTa2aFix6WismSwzsCrIuIf5F0LPAPkr4GfAC4NiK2SPoUyVj0nywwrSNlVJqiFN2To9/jshZHR+U5WjON7UuLBXN4kfiXdPXYdAlgHXBLuv0G4B2FpNDMyjemObxMUxNJWiTpXmA/sJ1kTPmnI+JoeshektmJ2p27UdK0pOknD87kkWYzK9psxqViMgW8iJhJJ8NdSTLP5GntDutw7mRETETExMknLeo/pWZWjkY7vCxLxfT0ljYink5nDno9sETSMWkur+1MRFUwbk0h8k5/ty50edwrj2v0Uw9Y9e+5VeMZPBwHc7lenm9pJa0HPkEyp8VnImJTy/7jgc8Dvw4cBH4rIh5N911F8n5gBvi9iJgaJC1Z3tKeLGlJ+vmXgXOAB4E7gAvSwzYAtw6SEDMbITnV4UlaBFwHvA14FXCxpFe1HHYJ8FRE/CvgWuCj6bmvIpnl7NeA9cBfptfrW5Yi7XLgDkk7SaZc2x4RXwWuAD4gaTdwEnD9IAkxs7G0FtgdEXsi4hlgC3B+yzHnk7z4hORF6NlKpi87H9gSEUci4ofA7vR6fVuwSBsRO0kmv23dvmfQm4+Cfos2ozJaShFFs36alGSdm7foomTWdFShKqOXNDb2rT33cC737qFIu0zSdNP6ZERMNq2vAB5vWt8LvK7lGs8dk07r+BOSTNQK4Dst57Z9OZqVe1qY2VxBL13LDkTERJf97S7UGk47HZPl3J5kektrZjWTXzu8vcCqpvV2LzifO0bSMcCvAIcyntsT5/DGWL9FybyLoMMsxo6roovkOb6l3QGskbQa+BHJS4h3tRyzleTF5z+SvAi9PSJC0lbgbyV9DHgJsAa4e5DEOOCZ2Xw5Bby0Tu5yYIqkWcrmiNgl6RpgOiK2krzw/Ov0BeghkqBIetzNwAPAUeB9ETFQ7wUHPDObL8d2eBGxDdjWsu3qps+/AC7scO5HgI/klRYHPDObQzG+w0N5Ep8R0k9Tl27PMe+mIqM00Gk/qvY3NtRJkjwAqJnVhXN4ZlYfDnj5yjrnaxX0UtTLWszMer/Wc7I+xzKLS1X+bmvJdXhmVisOeGZWF6rg4J5ZuGuZmdXGyOTwiq7b6VSPlEfzmDwms+n33qNaJ1blejs3mcJFWjOrCb+0MLNaccCrtiKLJf1ee5yLTsOsGhhUFb6Hwv92HPDMrA7E+L6ldcAzs7lchze+hll8qULRqZu8576o+vMoU+HPygHPzGrDAc/M6sJFWjOrDwc8s/64bq5iopy3tJKWAjcBpwKPAu+MiKdajjkD+CTwQmAG+EhE3JTu+xzwRuAn6eHviYiuf2zuS2tm8+U3TWM3VwK3RcQa4LZ0vdVh4D9FxK8B64GPS1rStP8PI+KMdFnwX9bMAU/SIknfk/TVdH21pLskPSLpJknHZb2WmY22xrwWCy0DOh+4If18A/CO1gMi4uGIeCT9vA/YD5zc7w17KdK+H3iQJGsJ8FHg2ojYIulTwCUkWc+OHt55Qq7zIJRdVBpWh/hR6YGQh3HuXTJMjef6cBzM54LZg9kySdNN65MRMZnx3BdFxBMAEfGEpFO6HSxpLXAc8IOmzR+RdDVpDjEijnS7RqYcnqSVwL8DPpOuC1gH3JIe0jY6m1kFZS3OJkHxQERMNC1zgp2kb0q6v81yfi9JkrQc+GvgvRHRqGG8Cngl8K+BpcAVC10naw7v48AfAS9I108Cno6Io+n6XmBFh4RuBDYCLOaEjLczs2ER+TVLiYhzOt5H+rGk5WnubjlJcbXdcS8E/h74rxHxnaZrP5F+PCLps8AHF0rPgjk8Sb8J7I+Ie5o3tzm07SOKiMlG9D+W4xe6nZmNgJLq8LYCG9LPG4Bb56UjeTfwZeDzEfHFln3L058iKWHev9ANs+TwzgLOk/R2YDFJHd7HgSWSjklzeSuBfRmu9Zx+5zjNe67VXgyrvqmXiX/66e7Vy/0G5Tq7iiinHd4m4GZJlwCPARcCSJoALouIS4F3Am8ATpL0nvS8RvOTL0g6mSQDdi9w2UI3XDDgRcRVJGVlJL0J+GBE/LakLwIXAFvoEJ3NrKJKCHgRcRA4u832aeDS9PPfAH/T4fx1vd5zkHZ4VwAfkLSbpE7v+gGuZWajImNxtordz3rqaRERdwJ3pp/3AGt7Of8VrznM1FRSpCm6eULVmz8UPcJI1Z6Hddf4PteeezifC1YwmGXhrmVmNo8HADWz2qhicTWLsQ143d4CV6E4V4U0VlnZVR55//11esueS0+LfPrJjqSxDXhmNgAHPDOrgzx7WowaBzwzm0ez4xnxxirg5dF8o2p1fVmN04greSj7d67UM3YdnpnViYu0ZlYfDnj1UKmix5gZld4xZQ/I0E81Sqfj8upp4RyemdWHA56Z1UJJs5YNgwOemc3hdngFcF3ZfEU2iel3wNWilZmOMpscVf7vO8Yz4jmHZ2bzOIdnZvXghsdWtnGan2OYvVeyFpnzKFr32/RkFIu/Zby0kLQUuAk4FXgUeGdEPNXmuBngvnT1sYg4L92+mmSKiaXAd4F3R8Qz3e45yBDvZjamNJttGdCVwG0RsYZ0Iu0Ox/08Is5Il/Oatn8UuDY9/yngkoVu6IBnZnMFyUuLLMtgzgduSD/fQDLVYibp1IzrgFt6Od9F2oL10mp/FIs2/epURCz77fCgPRfKTseo6OGlxTJJ003rkxExmfHcFzUm004n4z6lw3GL03scBTZFxFdIJg57Op0mFmAvsGKhGzrgmdl82QPegYiY6LRT0jeBF7fZ9aEeUvPSiNgn6eXA7ZLuA37a5rgFU+2AZ2Zz5NnwOCLO6Xgf6ceSlqe5u+XA/g7X2Jf+3CPpTuBM4EvAEknHpLm8lcC+hdLjOjwzmysCzWZbBrQV2JB+3gDc2nqApBMlHZ9+XgacBTwQEQHcAVzQ7fxWzuEVbFTqbnqpOxtWmkflWY2qUpv3lNMObxNws6RLgMeACwEkTQCXRcSlwGnApyXNkmTQNkXEA+n5VwBbJP034HvA9Qvd0AHPzOYpo6dFRBwEzm6zfRq4NP38beDVHc7fA6zt5Z6ZAp6kR4GfATPA0YiYyNpo0MwqJgDPacGbI+JA03qj0eAmSVem61fkmroh67cJxSgWzYoYPCDvJiajMgDoqCr1eYxnvBvopUXfjQbNbLQpsi1VkzXgBfANSfdI2phum9NoEGjbaFDSRknTkqafPDgzeIrNrHAlvaUtXdYi7Vlpw79TgO2S/inrDdJW15MAE6cvrt4TMqubuo+W0tTwb7+kL5O8GcnUaLCTKtTXjGKaWlXhOfZrXOcIHnVJw+PxjHgLFmklPU/SCxqfgbcC95Oh0aCZVdRsxqVisuTwXgR8ORmcgGOAv42Ir0vaQZtGg2ZWfeOaw1sw4KWN+05vs71to0ErV79Fvebz8mh+U/QoKHkUb11EzqjudXhmVifVfAObhQOemc1X1yKtmdWMJ+LO36jWofRTz9Ot/qpbl65ReQajko6iVf337OXvbGDO4ZlZbYxnvHPAM7P5NDueZVoHvBb9zC3aS3GiCsWqTr9nv0X3rKpQ/K+FoJKNirNwwDOzOUTUt+GxmdWQA54Nq1hV9gABo/J7Fnn9KhSRhzowRAkBL8uo6ZLeDFzbtOmVwEUR8RVJnwPeCPwk3feeiOj6kDxrmZnN1ajDK37wgMao6WuA29L1uUmJuCMizoiIM4B1wGHgG02H/GFj/0LBDhzwzKwNzc5mWgbU66jpFwBfi4jD/d7QAc/MWkRSpM2yDCbTqOlNLgJubNn2EUk7JV3bmL+2G9fhVdCo1EX122Qlq+Zr5vE7F/GsivwusjaRyv3eQS/BbJmk6ab1yXSUcwAkfRN4cZvzPtRLktJBhl8NTDVtvgr4Z+A4klHVrwCu6XYdBzwzmy97afVAREx02hkR53TaJ6mXUdPfCXw5Ip5tuvYT6ccjkj4LfHChxLpIa2bzKCLTMqBeRk2/mJbibBokUTI68TtIRmLvaqxyeFmLF6NSJMyq36JN0b9n0c1I8khz1jSWWUwusvj/cBzs6/x5ymmHt4k2o6ZLmgAui4hL0/VTgVXA/2k5/wuSTiaZhuNe4LKFbjhWAc/MchABM8X3Les0anpETAOXNq0/Cqxoc9y6Xu/pgGdm87mnhZnVxpgGPEWJv9jE6Yvj7qlVpd3PypPHREBFq1rdbTednvddcRs/jUMa5Nq/cvyL49+s+J1Mx379h//jnm5vaUeNc3hm1iIgxnN8KAc8M5srKOWlxTA44BWg32YHVS5y5TEAaFb99jKo2jPtpvDnPaZ1eA54ZjbfmAa8TD0tJC2RdIukf5L0oKTfkLRU0nZJj6Q/Tyw6sWZWhtIGDyhd1hzeJ4CvR8QFko4DTgD+mGQsq02SriQZy+qKgtJZC6My8GYevQeK/F2qXjQd6sCeWQQwppP4LJjDk/RC4A3A9QAR8UxEPE3vY1mZWVWMaQ4vS5H25cCTwGclfU/SZyQ9j4xjWUnaKGla0vSTB2dyS7iZFSXtWpZlqZgsAe8Y4LXAJyPiTOD/0WYo5k4iYjIiJiJi4uSTFvWZTDMrTUDEbKalarLU4e0F9kbEXen6LSQBr5exrMZO0SOFlClrHVIRo3xUXT9NiYp4Vo1rrj2379HP55qtXnE1iwVzeBHxz8Djkn413XQ28AC9jWVlZlUypnV4Wd/S/i7J2FPHAXuA95IEy3ljWZlZxUWM7VvaTAEvnf6sXQfheWNZDVOZPRWKns+hm2H1yBjVYuswe6iM6jMZWAVzb1m4p4WZtQhiZjxbVDjgmdlcwdi+tHDAM7P5KtjkJItKBLxOdWTDrD8ZZrOUsa036lPek/3U/fkGECXk8CRdCPwJcBqwNp3Lot1x60m6ty4CPhMRm9Ltq4EtwFLgu8C7I+KZbvf0NI1mNlekA4BmWQZzP/AfgG91OkDSIuA64G3Aq4CLJb0q3f1R4NqIWAM8BVyy0A0d8MxsnpiZybQMdI+IByPioQUOWwvsjog9ae5tC3B+OhftOpKOEJCxP3+pRdp7dh45sGj57v8LLAMOZD9zd9uti5Z3Pm7+vnl6TEO2NLUqPh25qXk65v3tVPV5vGzQG/6Mp6a+Gbcsy3j4YknNRdHJiJgcNA1NVgCPN63vBV4HnAQ8HRFHm7bPm8qxVakBLyJOBpA0PeyJP0YhDU6H0zGK6YiI9XldS9I3gRe32fWhiMjSO6vdhETRZXtXlXhpYWbVFBHnDHiJvUDzVIcrgX0kud4lko5Jc3mN7V25Ds/MRtkOYI2k1WnX1ouArZHML3sHcEF6XKb+/MMKeHmW8fs1CmkAp6OV0zHXqKQjd5L+vaS9wG8Afy9pKt3+EknbANLc2+XAFPAgcHNE7EovcQXwAUm7Ser0rl/wnmVOxG1mNkwu0ppZbTjgmVltlBrwJK2X9JCk3elMZ2Xdd7Ok/ZLub9pW+jSTklZJuiOd6nKXpPcPIy2SFku6W9L303T8abp9taS70nTclFYSF07SonS+lK8OKx2SHpV0n6R7G+3KhvQ34ilRC1RawFugi0jRPge0ti26kmSayTXAbfQwT8cAjgJ/EBGnAa8H3pc+g7LTcgRYFxGnA2cA6yW9nj666uTk/SQV0g3DSsebI+KMpnZvw/gbaUyJ+krgdJLnMox0jKeIKGUheRMz1bR+FXBVifc/Fbi/af0hYHn6eTnwUFlpaUrDrcBbhpkWkjmGv0vSev0AcEy776vA+68k+Z94HfBVkgalw0jHo8Cylm2lfi/AC4Efkr5MHFY6xnkps0jbrovIgl1BCpRpmsmiSDoVOBO4axhpSYuR95JMvrQd+AF9dNXJwceBPwIaPdH76jKUgwC+IekeSRvTbWV/LwNNiWoLKzPg9dUVZBxJej7wJeD3I+Knw0hDRMxExBkkOay1JEP0zDusyDRI+k1gf0Tc07y57HSkzoqI15JUubxP0htKuGergaZEtYWVGfA6dREZlh+n00tS5jSTko4lCXZfiIi/G2ZaACLiaeBOkjrFJZIa3Q3L+H7OAs6T9CjJKBjrSHJ8ZaeDiNiX/twPfJnkH4Gyv5d2U6K+dgjpGFtlBry2XURKvH+r0qeZTIe0uR54MCI+Nqy0SDpZ0pL08y8D55BUjvfcVWcQEXFVRKyMiFNJ/h5uj4jfLjsdkp4n6QWNz8BbScZqK/V7CU+JWrwyKwyBtwMPk9QXfajE+94IPAE8S/Kv6CUkdUW3AY+kP5eWkI5/S1I82wnT9olzAAAAb0lEQVTcmy5vLzstwGuA76XpuB+4Ot3+cuBukrGSvggcX+J39Cbgq8NIR3q/76fLrsbf5pD+Rs4AptPv5ivAicNIx7gu7lpmZrXhnhZmVhsOeGZWGw54ZlYbDnhmVhsOeGZWGw54ZlYbDnhmVhv/H/b/x7jhFeqoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%lprun -f RunMCRG RunMCRG(K,h)\n",
    "#RunMCRG(K,h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x214951ff0f0>]"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD8CAYAAACRkhiPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XlYnNXZ+PHvzR4IhAQIgUASyL6vZtFoNe7axlirNa1La9u8tlq1r/212s2qb/dNbV1q1Vpbd+NurEaNW0xiQvY9QBYIBAiEJRCY7fz+mGfIADPMAMOS4f5c11wZzjwznMkM537OfZZHjDEopZTqnyJ6uwJKKaV6jwYBpZTqxzQIKKVUP6ZBQCml+jENAkop1Y9pEFBKqX5Mg4BSSvVjGgSUUqof0yCglFL9WFRvVyCQ1NRUM2rUqN6uhlJKnTLy8vKOGmPSgjm2zweBUaNGsWHDht6uhlJKnTJE5GCwx2o6SCml+jENAkop1Y9pEFBKqX5Mg4BSSvVjGgSUUqof0yCglFL9mAYBpZTqxzQIKHWK+GhvBXvL6nq7GirMaBBQ6hRgjOHW5zbxg+c3o9cFV6GkQUCpU0BVvY3qBjs7SmpZW1jV29VRYUSDgFKngMKj9c33H/90fy/WRIUbDQJKnQIKK44DcPnM4by/u6z5Z6W6KmAQEJFsEVklIrtEZIeI3GqVPy8im63bARHZbJWPEpETXo894vVas0Vkm4jki8gDIiLd99aUCh+FR+uJiYzgjosnEB0RwT9XH+jtKqkwEUxPwAHcboyZCMwHbhKRScaYrxpjZhhjZgDLgZe9nlPgecwYc6NX+cPAMmCsdbsoNG9DqfBWWFHPyJR40pPiWDIzkxfziqhusPV2tVQYCBgEjDGlxpiN1v06YBcw3PO4dTZ/FfBse68jIhlAkjFmjXFPb3gKWNKFuivVbxRWHCc3LQGAby3MpdHu4ul1h3q5VipUVu0p59IHPuGEzdnjv7tDYwIiMgqYCazzKj4TKDPG7PMqyxGRTSLykYicaZUNB4q9jinGK5gopXxzOF0cqmogN20gAOOHJXLm2FT+9dkBbA5XL9dOhcJn+UfZUVLLx/sqevx3Bx0ERGQg7rTPbcaYWq+HltKyF1AKjDDGzAT+F3hGRJIAX/l/nxOeRWSZiGwQkQ0VFT3/n6JUX1J07AR2pyE3NaG57Ntn5lJe18SKbaW9WDMVrOoGG1uLq3G6fK/xOFTVAMC7O8p6slpAkFcWE5Fo3AHgaWPMy17lUcCXgdmeMmNME9Bk3c8TkQJgHO4z/yyvl80CSnz9PmPMo8CjAHPmzNGVMapf88wE8vQEAM4am0rmoDhWbCtlyUztUPdFL2wo4sM95Ww7XENR1QkAHvr6LC6ZmtHm2EPW4+/vLsPhdBEV2XMTN4OZHSTA48AuY8yfWz18HrDbGFPsdXyaiERa93NxDwAXGmNKgToRmW+95nXAayF6H0qFrf3WGoHRaSd7AiLCoolD+TT/KI32ns8jq/ZV1dv40UtbyTt4jGnDk/nRReOJENh9pO22H8YYiqoaGJ48gOoGO+sPHOvRugYTbs4ArgUWeU37vMR67GraDgifBWwVkS3AS8CNxhjPEsfvAo8B+UAB8HZX34BSHWGM4dN9p1bDWVBRz+D4aJLjY1qUnzsxnQabk3X7dQVxX7OluBqA+746kwe/PovvnT2G4YMHcMBr0Z/HsQY7x5scfG3eCGKjInh355EerWsws4M+NcaIMWaa17TPFdZj3zDGPNLq+OXGmMnGmOnGmFnGmDe8HttgjJlijBltjLnZ6CYoKgReyivmzyv3BnXsB7vLuebxdfz9o8JurlXouGcGDWxTviA3hQHRkby/q+fzyF2xs6SWgjBf7LalqBoRmJo1qLksJ3Vgc6/Om2c8YHy6e8D/3R1lPbo/lK4YVqe8lzcW8+CqfI7Vtz9v3u508esVuwB3vtbfIF1fU3i0vsWgsEdcdCQLx6by/q7yU2ZTOafL8M0nP+empzeeMnXujC1F1YwdOpCBsSeHXXNS4jlwtL7N+/YEgREp8VwwaRiHq0+ws7SWnqJBQJ3ySqpP4HQZ3tnRfjf6ufVFFFTUs2RGJoerT/Bp/tEeqmHn1TXaqahr8tkTADh3wlAOV59gzymyxfS6/ZWU1Tax+0hdjzZ0PckYw5biGqZnJbcoz0lNoK7JQWWrk5UiKwhkD47n3IlDiZCenSWkQUCd0owxlNY0AvBWO9Ml6xrt3LdyL3NzhvC7r0xjcHw0z33etxZb7Sur4+1W76Gwwp0+yE1r2xMAWDRhKADv7yrv3sqFyOubS4iPiSQ6Uliedzjkr29zuHA4e3ftRPGxE1TV25ie3TIIjLJ6c61TQocqG0hLjGVATCQpA2OZM3II7+7UIKBUUKrqbTQ5XAyOj+azgkq/KaFHPiqgst7Gzy6dSGxUJFfMymLlzjKOHm8K+ne5XIY7lm/l9he2hKr6Lfz+nT3c/OwmjlhBDU42GL7SQQBDk+KYljXolBgXaHI4WbGtlAsnD+PcCem8tvkw9hA32L94bTtn/O4Ddpb0Xi9jc5F7UHhGdtueAPgIAlUNjBgS3/zzBZPT2VVa29xD6G4aBEKspsHO+X/+iAdX5Xcq57m3rI5dQXaTd5bUcs4fP2T1KZDW6C6eXsC180f6TQmVVJ/gsU/2s2RGJtOsLvrVc7NxuAwvbyxuc7w/f1q5h+fWF7F8YzHbD9eE5g1YnC7D2sJKnC7Dc+tP9lAKK44TIe58sT+LJgxlU1E1lR0IaK9vKeF3/93dpTp31Md7j1Lb6GDxjEyumJ1FZb2Nj/aEdjHoR3srKKtt4qt/X8OagsqQvnawthRVExMVwfhhiS3KhycPIDpSfAaB7MEDmn8+f1I6QMD0ZqhoEAixHSU17Cs/zh/e2cMPX9zaoWX9R483sfTRtXz/2U0Bjz1hc3LLc5vYf7See97YecoMcoZaSbV7kc15k9IZMSTeZ0roD+/swQA/vHB8c9mYoYnMGTmY59YXBRWsX910mAdXFbBkRiYDY6P4xye+Zxd9sLvM5wyQQLYfrqGu0cHA2Cie+7yoOaVRcLSe7CHxxEZF+n3ueRPTMQZWdaBBffjDAh77pLBT207klx9nxj3vsrawY43sa5sPMyQhhoVjUjl7fBopCTEs70AQDqSstpHSmka+c2YO6YPiuP6Jz3lra8+vqN5SXM2UzCSiWy34ioqMIHtIfItpojaHi9KaEy16AiNTEpgwLLHHUkIaBEKswPqAr50/kuUbi7nm8XUBZ62AO7d9x/JtVNbbyC8/3iIl4MuvVuwkv/w41y0YyZ6yupD+MYWCMaZH9rXx9AQyBg3g0mkZbVJCK3eW8cqmw3x7YQ5Zg1ueTV89dwSFFfUBF+dsPHSMHy3fyrycIfz+K9O5+rRs3txaSvGxlt31j/ZWcMOTG7jwvo95+MOCDuWmVxe4e3M/uWQiR2obeX+3O8dfWOF7ZpC3yZlJpCfF8sHu4BqNkuoT7Cqtxe40nZqq+eznh6husLuDa5C93eNNDt7bVcalUzOIjowgOjKCxTMyeX9Xech2Q/WkYS6aksFLNy5gWtYgbn52Iy9sKArJ6wfD4XSx7XBNm/EAj9zUhBYnCSXVJ3AZyB7S8rt5weRhOF0m5OkyXzQIhFhhxXEGREdy9+LJ3H/1DDYfqubyh1YHvAjIi3nFvLerjCtmuXfWWFPoP8WzcmcZ/1l7iO+cmcPdiyczPTuZv6zc26MLoBxOFzUNdp+P1Tc5WPqPtVx438fUNvo+JlRKak4QExlBSkIMl07NaJESOlx9gh++uIUpw5O49byxbZ57ydRhJMZG+R0gbrQ7+Xx/FcueymNYUhyPXDObmKgIbliYg0CLPf3rmxz85OVt5KYlsGj8UH73390seWh10LnpNQWVjEsfyFVzssgYFMd/1h7E5TLsP+p7jYA3EWHRhHQ+3ns0qMC7as/JQeRgU48edqeLVzcdZtCAaPIOHgt6htXKnUdotLu4bEZmc9kVs7KwOV28EaKz9c1F1URFCJMzk0iOj+E/357H/JwU7n59B2W1vk+q6hrtIZ2quq/8OI12V5vxAI9RKQkcqKzHZfXcm6eHtgoCPzhvLMu/e3qb3kR30CDQCbc8u4kHV+X7fGz/0XpyUhOIiBAumzGcZ5fNo67RweUPfeY3R1lU1cA9b+xkfu4QfnfFVJLjo1md7/vY8tpGfrx8K5MykvjhheMREe68eAKlNY08+dmBUL3FNowxvL2tlF+8tp3LH1rN5LveYca97/Knd/e0OOM93uTgG//8nPUHjnGoqoGfv7q92+oEUFLdyLBBcURYf/wjU9wpIbvTxS3PbsLpMvxt6Syf6ZT4mCgum5nJW9tKufuNHc23//fiFi6+/xMm3/UOV/19DU0OJ49fP4fBCe4Vu5nJA/jitAye+/wQNSfcQe6P7+6hpOYEv79iGo9cO5uHvj6LIzWNLP7bp2w81H5Po8nhZP2BKk4fnUpUZARL547gk31HWbu/kka7y+/MIG/nThjK8SZHUCmaVbvLGZ48gNioiA4HgVW7y6mst/G7K6aSOSiO+97bF1Qj+trmEoYnD2DWiMHNZZMzk5gwLJHleSd7scYYSqpPdKph3lJUzcSMJOKi3Z91XHQkv/nyVOxOw2/fbjv+saagkln3rgxp7n2L1RtpPT3UY1RqAo12F2V17qDkvUbAW09eb0uDQAcVH2vg9S0l/He77y9OYUV9iz/a2SOH8OpNZ5CWGMu1j6/jhfUtu6ZOl+H2F7cgwB+vnE5UZAQLclNYU1DZ5g/BGPexDTYHDyyd2dywzc9N4ZzxaTy0Kr/bLjSyak853316I8vziomOjOCa+SNZPD2Tv36Qz9cfW0dZbSPHmxx885+fs/FQNfdfPYNbzx3La5tLeGVT96WqSqtPkJkcB7j/cC6Z6k4J3fX6DvIOHuPXX57aPDXPl+sWjGLQgGheyituvn24t4K0xFhu/EIuj1wzi1U/PJux6S0H+b5zVi71NifPrDvExkPHePKzA1w7fyRzRg0B4JKpGbxz21kAvBcgt7vpUDWNdhenj04B4KunZRMZIc0L23ICpIMAFo5NJT4mkrf9fC89Gu1OPs0/yrkThzIuPZFdpR1bX/BiXjGpA2M5b2I63ztnTFC9gcrjTXyy7yiLZ2QSEXGycRMRrpiVxeaial7cUMQvX9/BWX9Yxem//YCHPyroUL2cLsPW4po2Z+CjUhP4zlk5vLLpMBsOnNxeo6y2ke8/uxG707DxUHWHfpdHXaO9zd/bluJqBg2IZqSfgXxPam+/NfW3qKqBmMgI0hPjOlWHUAhqF1F1kqfx31tW12a3vyaHk+JjDW12dcweEs/y757Ozc9s5EfLt7JufxUDYyMpr2viUFUDO0pq+eOV05tz1qePSeXt7Uc4WNnQogHbUlzDJ/uO8rNLJzJmaMsUwY8vnsDF93/CQx8W8JNLJob8fb+UV8yQhBjW3nkuMVEn3/OZY9P4+avbueT+Txg+eAA7Smq5/+oZfHFaJk6X4ZN9Ffz81R3MGTmkTd4zFEprGpmXM6T550unZvDwhwU8s+4QS+dms3h6ZjvPhnHpiXz+0/M6/HsnZw5i4ZhU/rl6Py9vLCYjKY4fXTShxTEpA2OZkJHYnKv257P8o0QIzMt1B4H0pDgumJTe3KCPDpAOAvdZ77kT03lnxxHuvWyy310o1xS6exeLJgyl0e5sXm0czJnn0eNNrNpdzg0Lc4iKjODKOVk8tCqf+97bx8IxqX5fY8X2IzhdxudncdnMTH773938v5e2EhsVwcIxqaQNjOXBD/K5cnY2aYmxAesF7jTs8SaHz1z8TeeMYXneYe56fQev37wQlzHc/MxG6pucDEuKY4+PTd2CcfMzm9hXVscb319IykB3PTcX1TAta5Df/4vmtQKV9Zw+JpVDVQ1kDRnQIjj2NO0JdJDnD7PJ4eJAZcuBwYOVDbhMy90ePQYNiOaJb5zGtfNH8vKmYl7edJi9ZXUMGhDNDy8YxxWzTgYOzxnhZ63SR8+vLyIuOoKrTstu8/oThiXx5ZlZPLn6AD99ZRurdpeHbIygusHGezvLWTw9s0UAAPjK7Cxev/kMhiTEtAgAAJERwl++OgMBbn1uU8gX8ThdhiO1jWQknzyLmpyZxOi0BManJ/KLL04O6e9rbdlZ7j3995Uf51eXT22xRYDHjOxkthbXtDt767OCSqZmJTNoQHRz2dfnjQQgISaSoUE2hJdOzaCq3sbaQv8byq3aXc6A6Ejm56YwMSOJynobFXXBTS19ddNhHC7DlbPd41axUZHctChwb+DdHUcYneae8dLa0MQ4Hr12No9dN4fNv7iAx79xGn+8cjpNDhf3vRfcflAAm/zMzQd32u8nl05kR0ktz60/xO//u5v1B47x2yumsmB0Cns7sdq6usHGp/lHKalp5LbnN+N0GRpsDvaW1fkdDwAYlhRHXHRE8wyh1msEeoP2BDrgSE0jeQePcfGUYby9/Qi7j9S2OCNv3vc91feZW3RkBPcumcIvvjSp3QGf3NQEhiXFsbrgKF+bNwJwTwl9Y0sJl0zNICku2ufz7rh4AifsDl7ZdJin1x1iQHQkF08dxm+/PK1N4w3uhTXGwD2XTW73TPCNraXYnC6+MjvL5+Nj0xN54/sLqahranO2nzU4nv+7fAq3PreZX6/Yzc8unRiys56KuiacLkPGoJNzrEWE55YtICYqggEx/qdVhsKZY1M5c2wquakJnGOt3G1tZvZg/rP2EAUVxxmX3rYRrG9ysLmomu+cldui/PTRKeSmJpA4IDro/PDZ49NIiInkrW2lLByb2uZxYwzv7yrnjDGpxEVHMjEjCYCdpbUMTWo/HWGM4cUNxUzPTm6RGrtydjYPfuC/N9Bod493LJ07wu/7OHdieoufc9MGcs38kTy15gDfOH1Um1ScL5uLqkmMi/I7k+pL0zJ4eu1BfvXWLhpsTq5bMJLLZgynpLqRVzYdpuaEvUUQDuSD3eU4XYavzRvBM+sOcf/7+zhzbCpOl/E7HgAQESGMSjk5Q6ioqoHZIwf7Pb4naE+gA/673T2L4fuLxhIVIW0G1QqsPN+o1PYje6ARfxHh9DHucQHPLIIV20o53uTgqjltewEeaYmxPPT12Wz8+fk8+c3TWDIzk5c3HuZvPgax395WylNrDvLvtQf5z9qD7dbn5Y3FjE9PZHJmkt9j4qIj/aZ7LpsxnOsXjOSJ1ftZ9u+8kM0YKqlxrxHITG7ZgKUlxnboD7qzRISnbpjL3ZdN8XvMjBHuBmGzn7zz5weqcLgMZ4xu2WhHRAiPXT+HP105Pej6eKeEfPW69pUf53D1ieatJiYOc3+ewYwLbD9cy56yujYnAjFREXz37NHkHTzGNh8L6DYeOkaj3dXm/QVyy7ljSYiN4jc+BnR92VJUzfSsZL8nGCLCLxdPptHuZHp2Mj+91J0yHT/MfcK2r4O9gXd3lDEsKY7/u2wKV8zK4q8f7OMh6+9sWvagdp/rCQI1DXZqGx293hPol0Gg8ngTNz+zsd1Vn75mJ6zYfoRx6QOZlJnE6LSB7G71x7P/aD1DE2NJ9HOm3hGnj06lqt7WvDHYCxuKGJUS3yL/7U9cdCRnjx/Kb748jctnDuehVfkt3mt1g42fv7aDKcOTOHt8Gve+uYsdJb7/LwoqjrPpUDVXzB7epRkLv1w8mV9+aRKr9pSz5MHV5Jd3fSvh0uqTawR6S6D/k5yUBJLioprTFa19ln+UmMgIn2eDuWkD24z9BHJJOymhD6y1B+dMSANgUHw0w5MHBDVD6MW8ImKiIlg8rW1e/4vTMomMEJ+bnn2WX0lkhDAvN/D31tuQhBhuPmcMH+wuD7gi/oTNye4jdUwP0PhOzEji9ZsX8tQNc5snVYy3AqGvi73402h38tHeCs6flE5EhPB/S6YwPj2RVXsqGJ48gKEBBnlz0hI4VNVA4VH330B3jJV1RL8MAu/sKOPNraUs/cda8g62nL5X02Dne0/nsehPH1Fed3JucUVdE+sPVHHxFPel4SZkJLb54rj3fQ88kyMYnnGB1flHOXC0nnX7q7hyTnaHG+K7vjSJwQkx/PDFLc1zyH/11i6ONdj43RXT+NOV0xmcEM33n9nE8SZHm+e/vLGYCIElM7p2CUMR4Rtn5PD0t+dR02BnyYOru7ysv9TTE+jFIBBIRIQwPTuZTX6miX5WUMmskckhS115p4Ra+2BXOZMykloEzYkZiQGDwI6SGl7ddJgLJw9jUHzbE5zBCTHMGTmY93zsX/Rp/lGmZw3q1InR9aePYnjyAO59cyf/XnuQe97YyTf++Tlf+8daqrwWBO4ocY+5zMgOnFaZMnxQi15i5qA4BsZGdWhc4NN9Rzlhd3LBZHcaa0BMJA9fM5vE2KigUjs5KQnYnaY5UGtPoBdsOFhFcnw0KQkxXPv4uuYzjc1F1Vz61094d0cZJdUnuOnpjc0N5zs7jmAMzdcHnTAsicPVJ5rniYO173sQMzmCkZk8gJzUBNYUVPJiXhERQvNCso5Ijo/h15dPZfeROv62Kp9P9x3lxbxi/uesXCZnDiJlYCz3Xz2TA5X1/OyVbS16QC6X4ZWNhzlrXFrAnHGw5uem8Mb3F5KWGMtPX93WpcHiw9UniI+JJGlA3x7ampmdzN6yOupbBdlj9TZ2ltZyegdTJe3xlxKqabCTd+hYcyrIY2JGEoVH631OInC6DA9/WMCSB1cTFx3J984e7ff3nj8pnd1H6lpselbbaGdrcTULx3Tu/cVFR3LHxRPYfaSOn7+6nWc/P0R5rftk7K7XdzQf55l9Fagn4IuIMC59YIdmCL278wiJcVHMy0lpLstJTWDFrWdyz2WBJyPkWCeKH+91b/OhPYFekHfwGHNHDeGF/1lA9uB4vvnken7x2naufOQzjIEXb1zAH66czvoDx7j3zZ0AvL29lNy0BMaluxv5CRnuwSrPl6eq3kZ1gz3gEv+OOH10Cuv2V/FSXjFnjx/KsEGda4jPn5TenBa6/cXN5KYmcMu5J1fQzs9N4bbzxvHq5hJ++/bu5rOsNYWVlNQ0dir4tCczeQA/vmg8hRX1vLq5pM3je47UceFfPuYLf1jFJfd/wlWPrOEHz2+mydGyoSqtbiRjUFyPLqzpjBkjknEZ2uTM1xRWYgycMSbFzzM7p3VKyOUy/HvtAZwuw6KJbYOA02XYV9YyPVd8rIGl/1jL7/67m/MmpvPObWc1DyT7csGkYQAt9rtZW1CJy7inPHfWl6Zn8tYtC1lz5yJ23H0hK249k1sWjeWNLSXNY3Sbi6qDSsP4M35YInvL6oJaoOZ0Gd7bVc6iCUPbTLbIHhLf5hKgvoxKcbcRGw5WkZIQ43NWWU8K5kLz2SKySkR2icgOEbnVKv+liBz2cd1hROROEckXkT0icqFX+UVWWb6I3NE9b6l9FXVNHKxsYM6owQxNiuO5ZfOZMCyRp9Yc5AvjhvLWLQuZOWIwi6dnsuysXP699iCPfFTA2sIqLpmS0dzgTGzOJbq70p6ZQcHM6Q7W6aNTOd7koKy2iavmdK0h9qSFymqb+M2XpzavqvS46ZwxLJmRyd8/LmTBb97nJ69s4/FP95MYF9W8q2EoXTh5GJMzk3jg/X0t9kexO1387wubqTjexPSsZDKT42h0OHll02HyWu3xU1pzgszkvpsK8vDMFmm9XuD59UWkJcY272waKidTQiXsLKnlqr+v4Y/v7mVBbkqbmSueht07JdRgc3DFw5+x01q/8tDXZzWvlvZnREo849MTWel1fdzPCiqJi45g5oiuvb/JmYPIGHRyLv2NZ49myvAkfvbqdqrqbWwuqm53WmYg49ITOdZgD2qqbN7BY1TV25qDXmekDnQ3/HanIauXewEQ3BRRB3C7MWajiCQCeSKy0nrsL8aYP3ofLCKTcF+AfjKQCbwnIuOshx8EzgeKgfUi8roxZmco3kiwPGMAs0e6B6oGJ8Tw7Hfms+lQNWeMSWlxVvmjC8ezs6S2ecn5xVNPfvDpSbEkx0c3//EEuvhHZyywxgVSEmJYNKFrDXFyfAxPfvM09h+tb16U5C0yQrjv6pncdM4YHv90Py/lFWNzuFg6d0SbgBEKIsLtF4zjhic38FJeMUvnuqfCPvJhATtKannkmllcZI2/1DbamXH3u6zbX9XirLKkppEJw/yfnfYVKQNjGTEkvsUMob1ldXy0t4IfXjAu5PvDeFJCL288zPPri0iOj+H3X5nGV2ZltZk9M3JIPPExkS2u8vXP1Qcoq23ixRsXcNqo4Ad0z5s0lEc+KqS6wUZyfAyr848yNyel3R1QOyM6MoI/XjmdL/31U259bhPFx05w3YKRnX49z5bPe8rqAqY9391xhJjICL4wPq3Tv09EyElNYNvhml4fD4DgLjRfaozZaN2vA3YB7Y0SXgY8Z4xpMsbsB/KBudYt3xhTaIyxAc9Zx/aovINVxERFMGX4ycYjITaKhWPbznGOiozgr0tnkj1kAKPTEpjk1R0WESYMO7nsvvBoPdGRwvAQnpkOSYjhyzOH871zxvic599RkzMHNS/k8mdseiK/vWIaq3+8iF9+aRI/8LHxWqicM34oM7KT+ev7+2hyONlVWssDH+zjS9MzmwMAQFJcNJMyk1i3/+RAss3h4ujxphYLxfqymSOSW/QEnvh0P3HREXxtXucbr/ZcNcd9vYRr5o9k1e1nc9WcbJ/TJyMihPHDTg4OVzfYeOSjAs6bOLRDAQDg/EnunS8/2F1OWW0j+8qPc8bo0Ka6PCYMS+KWRWP5ZJ97PC+YQWF/xqe3TO36Y4zh3Z1lnDEmpcspHM/K4RFDer8n26F3IiKjgJnAOuAM4GYRuQ7YgLu3cAx3gFjr9bRiTgaNolbl8zpV6y7YcPAY07MGBX12MjghhrduOZNGu7NNkJgwLInn1xfhchkKK44zMiXB73L9zvrzV2eE9PWClZYYyzfOyOnW3+HpDVz7+Of8e81BXrF2p7x7cdvBtbmjUnh63UGaHE5ioyIpq23EmL49M8jbjOxkXttcwpGaRqIihZc3HebK2VkMCZBm6ayFY1PZfe9FQfUJOWARAAAZuklEQVQyJmYk8eaWEowxPPJRIcebHC2uvRCsacMHMTQxlpU7y/D8qZzRhfGAQG48ezTv7DzC7tK6Fid1HZUyMJbUgTFtZgjZHC7yDh5jcEI06YlxHKlt5FBVA99tZ4A8WDnNQaD3ewJBBwERGQgsB24zxtSKyMPAvYCx/v0TcAPga5TO4LvX4XMkRkSWAcsARowYEWwVA2q0O9l+uIZvLcwNfLCXpLhon6t0J2UkccLutOb8Bt73XbW1cEwqc3OG8OsVu3AZeOSaWT4bxnm5Q3hi9X62Ftdw2qghzReTOVV6Ap6c9eaiY+w+UofN4eKGhd0bZINNM03MSOKZdYfYXFTNk5/t57LpmZ1Ks0VECOdNSufVTYeJioxgcHx0i95zqEVHRvDYdaexp6yO+JiunZmPS09kT6vB8QdX5XP/+/uafxZx386b2PUxshxrQWlvzwyCIGcHiUg07gDwtDHmZQBjTJkxxmmMcQH/wJ3uAfcZvvey1iygpJ3yNowxjxpj5hhj5qSldT731trW4hrsThOyZdqeGULbS2o4WBm66aH9iYhw+/njcBnapIG8zbVSE5/vd8948awW7s2FYh0xKTOJmMgI1hZW8e81B1k0YWhIJxF0xSTre3z7i1twOA0/OH9cgGf4d/6kdBpsTt7cWsLpo1O7fWO0YYPi+MK4rrcR49IT2VdW17xC/4TNyb/WHGDhmFT+9rWZ/OKLk1h2Vi73LJ4c9KZ27Tl/0jB+eMG4DqfcukPA8CnuHMjjwC5jzJ+9yjOMMZ4VKZcDno3jXweeEZE/4x4YHgt8jruHMFZEcoDDuAePvxaqNxKMDQfdDUiogsDYoYlEiHurYLvThHRQuD+Zl5vCqzed4XODMY/BCTGMT09kbWElN50zhhJrtXDrLSP6qtioSCZmus+4bU4X3+7mXkBHeFbNFlbUc838EYxM6fz3+PTRKSTERFJvc3J6iKe+dqfxwxJpsDkpPnaCESnxvJhXRHWDnVvPG9stDfXA2ChuXtR9420dEUxP4AzgWmBRq+mgvxeRbSKyFTgH+AGAMWYH8AKwE/gvcJPVY3AANwPv4B5cfsE6tsfkHThGblpCyPKwA2IiGZWawHu73Mvxfe0eqoIzIzs54CykeblDyDt4DIfTfV3WQQOiu5wG6Ekzs5OxOV1MzEhqnvnVFwyMjWJkSjxx0RHc0sWGKTYqsnnmTEf3C+pN3jOEnC7DY5/sZ+aIZOb08uZuPSHgX5Ax5lN85/lXtPOcXwG/8lG+or3ndSeXy5B36BgXhHjO+8RhSRRWuDtEOX52D1WhMTdnCE+tOcj2ktrmhWKnkpkjknnyM/j2wpw+t8DtB+eNw2BCsjL8e2ePYezQRL8XVumLxlp7NHmuE3KoqoE7L57Q5z6n7nDqnEZ1UeHR41Q32JkzMrRduwnDEnlrWynJ8dHdNtNDuc3N8YwLuFcyh3I6bk+4eEoGslS4dKrvcY/e1PpCSF0xZfggpgzv+BYOvSkxzr2Z3u4jdazcWcbIlHgumNz5BWGnkn6zbcQGa7Xp7FGh7d5NsGY/6Myg7jc0MY7c1ATWFVZRWnPilJkZ5BETFcHi6e4dN1XfM35YIh/sKmNzUTXfXpjTbz6n/hMEDh5jcHx0yBtrz2CmzgzqGfNyh7C2sJLqBvspMzNInRrGpSdSb3MyOD6ar8z2f92OcNNvgkDewWPMHjk45Dm+rMED+MK4tJDMHVaBzc0ZQr3NvZHcqTIzSJ0aPCd01y4Y1e1XpetL+sWYQE2Dnf1H67myi5uw+SIi/OuGuYEPVCHhvX2v9gRUKJ0zfijXLRjJt7p5pXxf0y96Ap7LGaYN7PoiD9W7MpMHkDXY3fifKltGqFPDoPho7rlsis8L54SzfhEEbNZWxaHYhE31vnk5KYhA+iAN6kp1Vb9IB3muDhYT4s3dVO+46ZzRnD469FsUK9Uf9Ysg4LloSaj3bVe9IzdtoM7GUipE+kWr2NwT0HSQUkq10C9aRZv2BJRSyqd+0Sqe7An0jxWASikVrH4RBOxO9x7hMZE6kKiUUt76SRCw0kHaE1BKqRb6RRDQKaJKKeVbv2gVdWBYKaV86xetoqcnEKtTRJVSqoV+0SrqYjGllPKtX7SKulhMKaV8C9gqiki2iKwSkV0iskNEbrXK/yAiu0Vkq4i8IiLJVvkoETnhdVH6R7xea7Z1cfp8EXlAeugCntoTUEop34JpFR3A7caYicB84CYRmQSsBKYYY6YBe4E7vZ5TYIyZYd1u9Cp/GFgGjLVuF4XiTQRis9YJREfqFFGllPIWMAgYY0qNMRut+3XALmC4MeZdY4zDOmwt0O4VW0QkA0gyxqwxxhjgKWBJl2ofJJvDRUxkRMivKqaUUqe6DuVHRGQUMBNY1+qhG4C3vX7OEZFNIvKRiJxplQ0Hir2OKbbKfP2eZSKyQUQ2VFRUdKSKPtmdLu0FKKWUD0EHAREZCCwHbjPG1HqV/xR3yuhpq6gUGGGMmQn8L/CMiCQBvlph4+t3GWMeNcbMMcbMSUtLC7aKftkcLh0UVkopH4K6noCIROMOAE8bY172Kr8e+CJwrpXiwRjTBDRZ9/NEpAAYh/vM3ztllAWUhOJNBOLuCWgQUEqp1oKZHSTA48AuY8yfvcovAn4MLDbGNHiVp4lIpHU/F/cAcKExphSoE5H51mteB7wW0nfjh/YElFLKt2B6AmcA1wLbRGSzVfYT4AEgFlhpDbiutWYCnQXcIyIOwAncaIypsp73XeBJYADuMQTvcYRuY3O6dN8gpZTyIWAQMMZ8iu98/go/xy/HnTry9dgGYEpHKhgKNoemg5RSypd+0TLanZoOUkopX/pFy2h3Gp0iqpRSPvSLIKADw0op5Vu/aBltOkVUKaV86hcto83h0msJKKWUD/2iZdTFYkop5Vu/aBltOjtIKaV86hcto13XCSillE/9omW0OY0GAaWU8qFftIw2h1MHhpVSyod+0TLqYjGllPKtXwQBHRhWSinfwr5ldLoMTpeOCSillC9h3zLanS4A7QkopZQPYd8y2jxBQHsCSinVRti3jDaH9gSUUsqfsG8ZPekgHRNQSqm2wr5ltDsMoEFAKaV8CeZC89kiskpEdonIDhG51SofIiIrRWSf9e9gq1xE5AERyReRrSIyy+u1rreO3yci13ff2zrJ5nQCmg5SSilfgmkZHcDtxpiJwHzgJhGZBNwBvG+MGQu8b/0McDEw1rotAx4Gd9AA7gLmAXOBuzyBozvZrJ5AjC4WU0qpNgIGAWNMqTFmo3W/DtgFDAcuA/5lHfYvYIl1/zLgKeO2FkgWkQzgQmClMabKGHMMWAlcFNJ344NNp4gqpZRfHWoZRWQUMBNYB6QbY0rBHSiAodZhw4Eir6cVW2X+yruVDgwrpZR/QbeMIjIQWA7cZoypbe9QH2WmnXJfv2uZiGwQkQ0VFRXBVtGn5imiGgSUUqqNoFpGEYnGHQCeNsa8bBWXWWkerH/LrfJiINvr6VlASTvlbRhjHjXGzDHGzElLSwv2vfjkSQdFazpIKaXaCGZ2kACPA7uMMX/2euh1wDPD53rgNa/y66xZQvOBGitd9A5wgYgMtgaEL7DKupVdewJKKeVXVBDHnAFcC2wTkc1W2U+A3wIviMi3gEPAldZjK4BLgHygAfgmgDGmSkTuBdZbx91jjKkKybtohw4MK6WUfwGDgDHmU3zn8wHO9XG8AW7y81pPAE90pIJdpQPDSinlX9i3jLp3kFJK+Rf2LaPN6dk2QheLKaVUa+EfBKyeQGxkZC/XRCml+p6wDwLNYwJR2hNQSqnWwj4I6GIxpZTyL+xbRrvThQhERmhPQCmlWgv7IGBzuoiOjMC95k0ppZS38A8CDhexmgpSSimfwr51tDtdum+QUkr5Efato83h0kFhpZTyI+xbR7vT6PRQpZTyI+yDgPYElFLKv7BvHT2zg5RSSrUV9q2j3ekiVgeGlVLKp7BvHW0O7QkopZQ/Yd862jUdpJRSfoV962hzuPRaAkop5UfYt442p9GegFJK+RHMheafEJFyEdnuVfa8iGy2bgc81x4WkVEicsLrsUe8njNbRLaJSL6IPCA9tJmPzeHUgWGllPIjmAvNPwn8DXjKU2CM+arnvoj8CajxOr7AGDPDx+s8DCwD1uK+GP1FwNsdr3LH2J1GryqmlFJ+BDxFNsZ8DFT5esw6m78KeLa91xCRDCDJGLPGuhD9U8CSjle343RMQCml/Otq63gmUGaM2edVliMim0TkIxE50yobDhR7HVNslXU7nR2klFL+BZMOas9SWvYCSoERxphKEZkNvCoikwFf+Rjj70VFZBnu1BEjRozoUgV1xbBSSvnX6dZRRKKALwPPe8qMMU3GmErrfh5QAIzDfeaf5fX0LKDE32sbYx41xswxxsxJS0vrbBUB63oCmg5SSimfutI6ngfsNsY0p3lEJE1EIq37ucBYoNAYUwrUich8axzhOuC1LvzuoGk6SCml/AtmiuizwBpgvIgUi8i3rIeupu2A8FnAVhHZArwE3GiM8Qwqfxd4DMjH3UPo9plBDqcLl0EHhpVSyo+AYwLGmKV+yr/ho2w5sNzP8RuAKR2sX5fYne5hB+0JKKWUb2HdOtocLkB7Akop5U9Yt442pxUEdLGYUkr5FNZBwO7UnoBSSrUnrFtHTzpIxwSUUsq3sG4dPT0BDQJKKeVbWLeOTTowrJRS7Qrr1rF5TEB7Akop5VNYt446RVQppdoX1q2jLhZTSqn2hXXraHM6Ae0JKKWUP2HdOtocnp6ALhZTSilfwjoI6MCwUkq1L6xbRx0YVkqp9oV166iLxZRSqn1h3TradO8gpZRqV1i3jrp3kFJKtS+sW0dPT0CvMayUUr6Fdetod+hiMaWUak9Yt452p4vICCEyQtcJKKWUL8FcaP4JESkXke1eZb8UkcMistm6XeL12J0iki8ie0TkQq/yi6yyfBG5I/RvpS2b06ULxZRSqh3B9ASeBC7yUf4XY8wM67YCQEQmAVcDk63nPCQikSISCTwIXAxMApZax3Yrm8OlqSCllGpHVKADjDEfi8ioIF/vMuA5Y0wTsF9E8oG51mP5xphCABF5zjp2Z4dr3AE2p0sHhZVSqh1daSFvFpGtVrposFU2HCjyOqbYKvNX7pOILBORDSKyoaKiotMVtGtPQCml2tXZFvJhYDQwAygF/mSV+0rAm3bKfTLGPGqMmWOMmZOWltbJKrp7ArpQTCml/AuYDvLFGFPmuS8i/wDetH4sBrK9Ds0CSqz7/sq7jd2pPQGllGpPp1pIEcnw+vFywDNz6HXgahGJFZEcYCzwObAeGCsiOSISg3vw+PXOVzs4NofRHUSVUqodAXsCIvIscDaQKiLFwF3A2SIyA3dK5wDwPwDGmB0i8gLuAV8HcJMxxmm9zs3AO0Ak8IQxZkfI300rNqeLaE0HKaWUX8HMDlrqo/jxdo7/FfArH+UrgBUdql0X2R0uYrUnoJRSfoV1C+nuCehiMaWU8iesg4AODCulVPvCuoW0OVw6MKyUUu0I6xZSB4aVUqp9Yd1C2nRgWCml2hXWLaSOCSilVPvCuoW0O41uG6GUUu0I6xZSt5JWSqn2hXULqesElFKqfWEbBIwxOjCslFIBhG0L6XDpReaVUiqQsG0hbQ4XgA4MK6VUO8K2hbQ73UFAewJKKeVf2LaQNqf2BJRSKpCwbSGb00HaE1BKKb/CtoW0O90Dw9oTUEop/8K2hfT0BHRMQCml/AvbFvLkwLAuFlNKKX8CBgEReUJEykVku1fZH0Rkt4hsFZFXRCTZKh8lIidEZLN1e8TrObNFZJuI5IvIAyLSra1zk04RVUqpgIJpIZ8ELmpVthKYYoyZBuwF7vR6rMAYM8O63ehV/jCwDBhr3Vq/Zkh5egI6MKyUUv4FbCGNMR8DVa3K3jXGOKwf1wJZ7b2GiGQAScaYNcYYAzwFLOlclYOji8WUUiqwULSQNwBve/2cIyKbROQjETnTKhsOFHsdU2yVdRtdLKaUUoFFdeXJIvJTwAE8bRWVAiOMMZUiMht4VUQmA77y/6ad112GO3XEiBEjOlU3uy4WU0qpgDrdQorI9cAXga9bKR6MMU3GmErrfh5QAIzDfebvnTLKAkr8vbYx5lFjzBxjzJy0tLRO1a9Jp4gqpVRAnWohReQi4MfAYmNMg1d5mohEWvdzcQ8AFxpjSoE6EZlvzQq6Dnity7Vvh2exWKz2BJRSyq+A6SAReRY4G0gVkWLgLtyzgWKBldZMz7XWTKCzgHtExAE4gRuNMZ5B5e/inmk0APcYgvc4QsjpYjGllAosYBAwxiz1Ufy4n2OXA8v9PLYBmNKh2nWBLhZTSqnAwvY0WaeIKqVUYGHbQtp0iqhSSgUUti2krhhWSqnAwraFtDlcREUIERE6JqCUUv6EbRCwO106HqCUUgGEbStpc7h0PEAppQII21bS5jQaBJRSKoCwbSVtDpeuFlZKqQDCtpW0O126UEwppQII2yBgc+jAsFJKBRK2raS7JxC2b08ppUIibFtJm04RVUqpgMK2ldQpokopFVjYtpJ2p84OUkqpQMK2lbTpmIBSSgUUtq2k3WF0iqhSSgUQtkHAPTAc2dvVUEqpPi18g4BDF4sppVQgQQUBEXlCRMpFZLtX2RARWSki+6x/B1vlIiIPiEi+iGwVkVlez7neOn6fiFwf+rdzkg4MK6VUYMG2kk8CF7UquwN43xgzFnjf+hngYmCsdVsGPAzuoIH7IvXzgLnAXZ7A0R10YFgppQILqpU0xnwMVLUqvgz4l3X/X8ASr/KnjNtaIFlEMoALgZXGmCpjzDFgJW0DS8jYHS69qphSSgXQlVYy3RhTCmD9O9QqHw4UeR1XbJX5K+8W509KZ/LwpO56eaWUCgtR3fCavkZjTTvlbV9AZBnuVBIjRozoVCXuu3pmp56nlFL9SVd6AmVWmgfr33KrvBjI9jouCyhpp7wNY8yjxpg5xpg5aWlpXaiiUkqp9nQlCLwOeGb4XA+85lV+nTVLaD5QY6WL3gEuEJHB1oDwBVaZUkqpXhJUOkhEngXOBlJFpBj3LJ/fAi+IyLeAQ8CV1uErgEuAfKAB+CaAMaZKRO4F1lvH3WOMaT3YrJRSqgeJMT7T8n3GnDlzzIYNG3q7GkopdcoQkTxjzJxgjtU5lEop1Y9pEFBKqX5Mg4BSSvVjGgSUUqof6/MDwyJSARzs5NNTgaMhrE6o9NV6Qd+tW1+tF/TduvXVekHfrVtfrRd0rG4jjTFBLbLq80GgK0RkQ7Aj5D2pr9YL+m7d+mq9oO/Wra/WC/pu3fpqvaD76qbpIKWU6sc0CCilVD8W7kHg0d6ugB99tV7Qd+vWV+sFfbdufbVe0Hfr1lfrBd1Ut7AeE1BKKdW+cO8JKKWUakdYBgERuUhE9ljXOb4j8DO6tS5BX5+5h+uVLSKrRGSXiOwQkVv7UN3iRORzEdli1e1uqzxHRNZZdXteRGJ6um5WPSJFZJOIvNnH6nVARLaJyGYR2WCV9YXPM1lEXhKR3db3bUEfqdd46//Kc6sVkdv6SN1+YH33t4vIs9bfRLd8z8IuCIhIJPAg7msdTwKWisikXqzSkwR/feae5ABuN8ZMBOYDN1n/T32hbk3AImPMdGAGcJG1LfnvgL9YdTsGfKsX6gZwK7DL6+e+Ui+Ac4wxM7ymEvaFz/N+4L/GmAnAdNz/d71eL2PMHuv/agYwG/eux6/0dt1EZDhwCzDHGDMFiASupru+Z8aYsLoBC4B3vH6+E7izl+s0Ctju9fMeIMO6nwHs6QP/b68B5/e1ugHxwEZgHu6FMlG+PucerE8W7oZhEfAm7ivm9Xq9rN99AEhtVdarnyeQBOzHGn/sK/XyUc8LgNV9oW6cvBTvENzb/b+J+xrt3fI9C7ueAD18LeNO8nd95l4hIqOAmcA6+kjdrJTLZtxXrFsJFADVxhiHdUhvfa73AT8CXNbPKX2kXuC+XOu7IpJnXaIVev/zzAUqgH9aKbTHRCShD9SrtauBZ637vVo3Y8xh4I+4r9NSCtQAeXTT9ywcg0DQ1zJWICIDgeXAbcaY2t6uj4cxxmnc3fQsYC4w0ddhPVknEfkiUG6MyfMu9nFob33fzjDGzMKdCr1JRM7qpXp4iwJmAQ8bY2YC9fROSsovK7e+GHixt+sCYI1BXAbkAJlAAu7PtLWQfM/CMQgEfS3jXuTv+sw9SkSicQeAp40xL/elunkYY6qBD3GPWySLiOdqeL3xuZ4BLBaRA8BzuFNC9/WBegFgjCmx/i3HndueS+9/nsVAsTFmnfXzS7iDQm/Xy9vFwEZjTJn1c2/X7TxgvzGmwhhjB14GTqebvmfhGATWA2OtkfQY3N2813u5Tq35uz5zjxERAR4Hdhlj/tzH6pYmIsnW/QG4/yh2AauAr/RW3Ywxdxpjsowxo3B/rz4wxny9t+sFICIJIpLouY87x72dXv48jTFHgCIRGW8VnQvs7O16tbKUk6kg6P26HQLmi0i89Xfq+T/rnu9Zbw7GdOPAyiXAXtx55J/2cl2exZ3Xs+M+K/oW7jzy+8A+698hvVCvhbi7k1uBzdbtkj5St2nAJqtu24FfWOW5wOe4r1/9IhDbi5/r2cCbfaVeVh22WLcdnu99H/k8ZwAbrM/zVWBwX6iXVbd4oBIY5FXW63UD7gZ2W9//fwOx3fU90xXDSinVj4VjOkgppVSQNAgopVQ/pkFAKaX6MQ0CSinVj2kQUEqpfkyDgFJK9WMaBJRSqh/TIKCUUv3Y/wcMtSkHhJ3jRAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(energy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}