1D-Ising-Model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1D Ising Model\n",
    "\n",
    "This notebook will guide you through some aspects of the one-dimensional Ising model. The Ising Model originated from magnetism, but has many other applications as well. It models the interactions between neighbooring sites as\n",
    "\n",
    "\n",
    "$$\n",
    "H(\\{\\sigma\\}) = \\sum_{i,j} J \\cdot \\sigma_i\\sigma_j - \\mu \\sigma_i \\quad \\sigma_i\\in\\{-1,+1\\}\n",
    "$$\n",
    "\n",
    "Each site can be in one of two states (e.g., spin up $\\uparrow$ and spin down $\\downarrow$). If the exchange energy $J$ is positive, anti-aligned spin pairs ($\\uparrow\\downarrow$) are favoured, otherwise aligned spin pairs ($\\uparrow\\uparrow$) are energetically advantageous. The point term represents interactions with an external field $\\mu$. The spin variables $\\sigma_i$ can take on two values: $+1$ and $-1$, respectively. These choices to designate the spin states are not unique, but they make some of the math simpler."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Helping functions\n",
    "\n",
    "Before turning to solution schemes for the Ising model, we load needed libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import time\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and define some helpful functions.\n",
    "\n",
    "### Lattice state iterator\n",
    "\n",
    "For some of the solution approaches, we need to iterate over all possible states of the lattice. This is best done with an iterator. Itertors are a handy feature of python that allow to generate sequences of values that are returned by `yield` each time the iterator is accessed.\n",
    "\n",
    "This iterator returns an array containing spin variables $[\\pm 1,\\pm 1,\\dots,\\pm 1]$. \n",
    "\n",
    "> The implementation utilised the bit representation of integers, e.g., 3 -> 00000010, which allows to iterate through all possible states by simply incrementing a number.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 i =  (array([-1.]), 0)\n",
      "1 i =  (array([1.]), 1)\n"
     ]
    }
   ],
   "source": [
    "## Create a 1D lattice of spin variables\n",
    "def lattice_iterator(N):\n",
    "    '''\n",
    "    Create an iterator that iterates over all possible states {-1,+1} of a set of spin variables\n",
    "    \n",
    "    This simple implementation expoits binary representations of intergers with some padding\n",
    "    '''\n",
    "    \n",
    "    n=-1\n",
    "    M=2**N\n",
    "    d = np.zeros(N) - 1;\n",
    "    \n",
    "    while n<(M-1):\n",
    "        n=n+1\n",
    "        bit = [1 if digit == '1' else -1 for digit in bin(n)[2:]]\n",
    "        d[N-len(bit):] = bit\n",
    "        yield d, n\n",
    "        \n",
    "for n,i in enumerate(lattice_iterator(1)):\n",
    "    print(n, 'i = ', i)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Operators\n",
    "\n",
    "The matrix transform formalism explained below utilises **operators** acting on a function space. In computational studies, these operators take the form of matrices (similar to the formalism that you have already seen in the FEM method). Generally, an operator matrix contains the value of the operator for two states that are combined with each other through the operator in each of its matrix cells.\n",
    "\n",
    "### Magnetisation\n",
    "\n",
    "The Magnetisation operator only depends the state at hand. Drawing another analogy to the FEM formalism, it is similar to a linear form (while the Hamiltonian that we look at next is similar to a bilinear form). The matrix is of size $2^M\\times 2^M$, because this is the numer of possible states for a chain of length $M$. The magnetisation operator holds the total spin on its main diagonal for each possible spin configuration, which is constructed using the `lattice_iterator` that has been defined earlier. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-2.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  2.]]\n"
     ]
    }
   ],
   "source": [
    "def magnetisation(M):\n",
    "    m=np.zeros((2**M,2**M))\n",
    "    for s,i in lattice_iterator(M):\n",
    "        m[i,i] = np.sum(s)\n",
    "    return np.matrix(m)\n",
    "\n",
    "print(magnetisation(2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hamiltonian\n",
    "\n",
    "The Hamiltonian is a sum over point and pair terms. The `point` function sums all spin values over an array of spins. The `pair` function sums products of spin variables $\\sigma_a\\sigma_b$ over pairs of neighbouring lattices. The benefit of this approach will become apparent later. It can be thought of as \"summing over the rugs of a ladder\".\n",
    "\n",
    "The `hamiltonian` itself computes the energy of all possible arrangements of two spin chains, and returns the resulting energies as a matrix. For, example the Hamiltonian matrix with just one spin per chain looks like this:\n",
    "\n",
    "$$\n",
    "\\begin{array}{r|cc}\n",
    "           & \\uparrow             & \\downarrow\\\\\n",
    "  \\hline\n",
    "  \\uparrow & J(+1)(+1) - \\mu(+1) & J(+1)(-1) - \\mu(-1)\\\\\n",
    "  \\downarrow & J(+1)(-1) - \\mu(+1) & J(-1)(-1) - \\mu(-1)\\\\\n",
    "\\end{array}\n",
    "\\quad\\rightarrow\\quad{\\bf H}_1 = \\left[\n",
    "\\begin{array}{cc}\n",
    "J-\\mu & -J+\\mu \\\\\n",
    "-J-\\mu & J+\\mu\\\\\n",
    "\\end{array}\n",
    "\\right]\n",
    "$$\n",
    "\n",
    "In this representation, spin states are described by a state vector $|\\sigma\\rangle$ (using Dirac notation, which we will learn more about later). For now, you can think of it as a fancy way of writing a vector:\n",
    "\n",
    "$$\n",
    "\\begin{array}{rcl}\n",
    "|\\uparrow\\rangle & = & [1,0]\\\\\n",
    "|\\downarrow\\rangle & = & [0,1]\\\\\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "The elements of the Hamiltonian, are then given by $H_{ij} = \\langle\\sigma_i|H|\\sigma_j\\rangle$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def point(lat,na):\n",
    "    '''\n",
    "    Sum the spins over the lattice\n",
    "    '''\n",
    "    res = 0\n",
    "    for i in range(na):\n",
    "        res += lat[i]\n",
    "\n",
    "    return res\n",
    "\n",
    "## Define functions for summing over lattice states for \n",
    "## intra-layer and trans-layer interactions\n",
    "def pair(lat_a,lat_b,n):\n",
    "    '''\n",
    "    Sum intra-layer interactions over a lattice state\n",
    "    '''\n",
    "    res = 0\n",
    "    for i in range(n):\n",
    "        res += lat_a[i]*lat_b[i]\n",
    "    return res\n",
    "\n",
    "def hamiltonian(J,mu,M):\n",
    "    H=np.zeros((2**M,2**M))\n",
    "    for s,i in lattice_iterator(M):\n",
    "        for t,j in lattice_iterator(M):\n",
    "            H[i,j] = J*pair(s,t,M) - mu*point(s,M)\n",
    "    return H"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solutions\n",
    "\n",
    "### Brute force\n",
    "\n",
    "The aim is to compute the thermoodynamic properties of the spin chain, namely the free energy and average spin as a function of temperature and the model parameters. Statistical Mechanics provides the tools to establish this link.\n",
    "\n",
    "The central quantity in Statistical Mechanics is the **Partition Function**:\n",
    "\n",
    "$$\n",
    "Z = \\sum_{\\{\\sigma\\}}\\exp(-\\beta\\cdot H(\\{\\sigma\\}))\n",
    "$$\n",
    "\n",
    "The partition function sums over *all possible* states $\\{\\sigma\\}$ of the spin chain, weighting each state by the exponential of the states energy $H(\\{\\sigma\\})$ divided by $kT$. It is common to use $\\beta = 1/kT$ to save some writing. This quickly becomes computationally *very* expensive. The 1D spin-chain with $N$ sites has $2^N$ possible configurations. To make this explicit, a chain with 4 sites has 16 states. A chain with 10 sites has 1024 sites, and a chain with 20 sites has 1,048,576 states.\n",
    "\n",
    "The free energy $g$ is directly computed from the partition function using the thermodynamic identity\n",
    "\n",
    "$$\n",
    "g = -kT\\cdot\\ln(Z)/N\n",
    "$$\n",
    "\n",
    "The internal energy and magnitisation is computed using explicit summation:\n",
    "\n",
    "$$\n",
    "u = \\langle H\\rangle = \\sum_{\\{\\sigma\\}}H_i\\cdot\\exp\\{-\\beta\\cdot H_i\\}\\\\\n",
    "m = \\langle \\sigma\\rangle = \\sum_{\\{\\sigma\\}}\\left(\\sum\\sigma\\right)\\cdot\\exp\\{-\\beta\\cdot H_i\\}\n",
    "$$\n",
    "\n",
    "where $\\left(\\sum\\sigma\\right)$ is the sum of spins for a given configuration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def brute_force(M,J,mu,kTs):\n",
    "    '''\n",
    "    Solve the one-dimensional Ising model by explicitly summing over all states to compute\n",
    "    partition function, free energy, internal energy, and magnetisation.\n",
    "    \n",
    "    The function returns the results as a Pandas dataframe.\n",
    "    \n",
    "    Parameters:\n",
    "        \n",
    "        M   - length of spin chain\n",
    "        J   - exchange interaction\n",
    "        mu  - external field (same units as J)\n",
    "        kTs - array with temperatures (in units of J)\n",
    "    \n",
    "    Returns:\n",
    "    \n",
    "        A Pandas dataframe with free energy g, internal energy u, magnetisation m\n",
    "        \n",
    "    Usage: \n",
    "        \n",
    "        brute_force(10,-1,0.1,np.linspace(0.1,4,30))\n",
    "    '''\n",
    "    \n",
    "    # Create a new data-frame to take results\n",
    "    df = pd.DataFrame(columns=['kT','g','u','m','mu','J'],dtype=object)\n",
    "    \n",
    "    # Iterate over temperatures\n",
    "    for i,kT in enumerate(kTs):\n",
    "        beta =1/kT\n",
    "        Z = m = g = u = 0\n",
    "        \n",
    "        # Iterate over all possible spin states\n",
    "        # of a spin chain with length M\n",
    "        for l,j in lattice_iterator(M):\n",
    "            H = J*pair(l,np.roll(l,1),M)-mu*point(l,M)\n",
    "            Z += np.exp(-beta*H)\n",
    "            m += np.sum(l)*np.exp(-beta*H)\n",
    "            u += H*np.exp(-beta*H)\n",
    "\n",
    "        # Normalise sums and compute free energy\n",
    "        m = m/Z/M\n",
    "        u = u/Z/M\n",
    "        g = -1/beta*np.log(Z)/M\n",
    "        \n",
    "        # Append result to dataframe\n",
    "        s = pd.DataFrame({'g': [g], 'u': [u], 'kT': [kT], 'm': [m], 'mu': [mu], 'J': [J], 'M':[M]})\n",
    "        df = pd.concat([df,s])\n",
    " \n",
    "    return df\n",
    "\n",
    "#brute_force(10,-1,0.1,np.linspace(0.1,4,30))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Execise 1: AFM solution of the 1D Ising model\n",
    "\n",
    "Let's investigate a simple 1D Ising chain with direct summation of the partition function. We will look at the anti-ferro-magnetic solution by setting $J<0$ and vary the "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnQAAAJRCAYAAADMCCpKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA6pklEQVR4nO3deZzkdX3v+9enezZmpocBZkBgGCBKVBIRcAQ9auJGAkaDZjmRaB7c5OTwINEb4zlZyM2enCUmJ8m5OSExmGOOua5ZJHBycEGjEmNQYBiQTUV0mGFQupgZqJqZru7q/tw/qnpoht6qp2uqfr/f6/l4tF31W6o+vYz95rtGZiJJkqTiGup3AZIkSTo6BjpJkqSCM9BJkiQVnIFOkiSp4Ax0kiRJBWegkyRJKrgV/S6gnzZt2pRnnXVWv8uQJEla0B133FHLzM2znat0oDvrrLO4/fbb+12GJEnSgiJi51zn7HKVJEkqOAOdJElSwRnoJEmSCs5AJ0mSVHAGOkmSpIIz0EmSJBWcgU6SJKngDHSSJEkFZ6CTJEkqOAOdJElSwRnoJEmSCs5AJ0mSVHAGOkmSpIIz0EmSJBWcgU6SJKngDHSSJEkFZ6CTJEkqOAOdJElSwRnoJEmSCm6gAl1EvDciHouIe+Y4HxHxJxHxYETcHREXzjh3aUR8pXPummNXtSRJUn8NVKAD/hdw6TznLwPO6XxcBfw5QEQMA9d2zp8LXBER5/a0UkmSpAExUIEuM28B9s5zyeXAX2fbrcDGiDgVuAh4MDMfysxx4MOdayVJkkpvRb8L6NLpwK4Zz3d3js12/OJjWJckqeTGW1N87qujfOHrNRpjLQ6OT3JgvMWBZvvxVEJmkglTmUxl9rtkHUO/9vpzedVzT+7b+xct0MUsx3Ke4898gYiraHfXsnXr1uWrTJJUOlNTyZe+uZcbdjzCTV/+Fk8cmuC4lcOcsHYla1evYN2qYdauWsEpG1YyPBQEMBTB0BAEMftfJ5XShjUr+/r+RQt0u4EzZjzfAuwBVs1x/Bky8zrgOoBt27b5n0+SpFn969cf550f2cG3nhxj7aphvv+7nsUPnn8aL3/OJlYOD9SIJalwge5G4O0R8WHaXapPZOajETEKnBMRZwOPAG8GfryPdUqSCu5dH3+A4aHgT664gNc+/2TWriran0xVyUD9dkbEh4BXApsiYjfwm8BKgMx8N3AT8DrgQeAg8JOdc62IeDvwCWAYeG9m3nvMvwBJUinc88gT7Ni1n998w7n84AtP63c50oIGKtBl5hULnE/gbXOcu4l24JMk6ai8/9adHLdymB+6cEu/S5EWxUEAkiTN8OTYBDfs2MPl55/G8cf1d6C7tFgGOkmSZvjoHbs5NDHJW19yZr9LkRbNQCdJUkdm8v4vPswLz9jId59+fL/LkRbNQCdJUsetD+3lwccavPVi1ylVsRjoJEnqeP8Xd3L8cSt5gzNbVTAGOkmSgMeeHOMT93yLH33RFtasHO53OVJXDHSSJAEfuW0XrankLU6GUAEZ6CRJlTc5lXzoSw/zinM2cfamdf0uR+qagU6SVHn/9MBj7HlijLdcbOucislAJ0mqvI99+VFOWreK1z7/5H6XIi2JgU6SVHnfenKMM09ay4ph/yyqmPzNlSRVXq3RZNP61f0uQ1oyA50kqfJqjXE2jxjoVFwGOklSpbUmp9h3cNwWOhWagU6SVGl7D4yTCZtsoVOBGegkSZX2WL0JwOb1q/pcibR0BjpJUqXVGu1AZ5erisxAJ0mqtFpjHDDQqdgMdJKkSptuoXOWq4rMQCdJqrRavclxK4dZt3pFv0uRlsxAJ0mqtFqjyaYRJ0So2Ax0kqRKG3WXCJWAgU6SVGm1uosKq/gMdJKkSnMfV5WBgU6SVFmtySn2HnQfVxWfgU6SVFl7D7a3/XKXCBWdgU6SVFm1uosKqxwMdJKkyhqd3vbLLlcVnIFOklRZtbr7uKocDHSSpMqa3vZrk2PoVHAGOklSZdUaTdasHGK9236p4Ax0kqTKqjXaiwpHRL9LkY6KgU6SVFmjdRcVVjkY6CRJleUuESoLA50kqbJqjSabR5wQoeIz0EmSKmlyKtl7YNwWOpWCgU6SVEl7D4wzlbiPq0rBQCdJqqSn1qAz0Kn4DHSSpEoadZcIlYiBTpJUSe4SoTIx0EmSKulwoHMMnUrAQCdJqqRaY5xVK4YYcdsvlYCBTpJUSbV6k81u+6WSMNBJkipptNG0u1WlYaCTJFXSaL3JZidEqCQMdJKkSqo13CVC5WGgkyRVTnvbr6aBTqVhoJMkVc6+g277pXIx0EmSKsdtv1Q2BjpJUuXU6uOAu0SoPAx0kqTKGW2MAe4SofIw0EmSKuepFjoDncrBQCdJqpxao8mq4SE2rHHbL5XDQAW6iLg0Ir4SEQ9GxDWznP/FiNjR+bgnIiYj4sTOuW9GxJc7524/9tVLkopitNFk84jbfqk8BuY/TSJiGLgWuATYDdwWETdm5n3T12TmHwB/0Ln+DcA7M3PvjJd5VWbWjmHZkqQCai8q7IQIlccgtdBdBDyYmQ9l5jjwYeDyea6/AvjQMalMklQqtbqLCqtcBinQnQ7smvF8d+fYM0TEWuBS4O9nHE7gkxFxR0Rc1bMqJUmFN9ow0KlcBqbLFZhtIEPOce0bgH85orv1ZZm5JyJOBm6OiAcy85ZnvEk77F0FsHXr1qOtWZJUMFNTyd4D42wasctV5TFILXS7gTNmPN8C7Jnj2jdzRHdrZu7pfH4MuJ52F+4zZOZ1mbktM7dt3rz5qIuWJBXLvoPjTE6lLXQqlUEKdLcB50TE2RGxinZou/HIiyLieOB7gRtmHFsXESPTj4HvA+45JlVLkgql1mivQec+riqTgelyzcxWRLwd+AQwDLw3M++NiKs759/dufRNwCcz88CM208Bru9MP18BfDAzP37sqpckFYX7uKqMBibQAWTmTcBNRxx79xHP/xfwv4449hDwwh6XJ0kqAQOdymiQulwlSeq50Xo70G020KlEDHSSpEoZnd7267iB6qSSjoqBTpJUKbX6OCetX+W2XyoVA50kqVJqnX1cpTIx0EmSKmX/wXE2rnVRYZWLgU6SVCn1ZouRNY6fU7kY6CRJlVIfa7HBQKeSMdBJkiqlMdZi/WoDncrFQCdJqoyJySkOTUwysmZlv0uRlpWBTpJUGQeaLQBb6FQ6BjpJUmXUx9qBzkkRKhsDnSSpMgx0KisDnSSpMupjEwCOoVPpGOgkSZXRcAydSspAJ0mqDLtcVVYGOklSZdSnW+gMdCoZA50kqTKmx9BtcAydSsZAJ0mqjMZYixVDweoV/vlTufgbLUmqjPpYi5E1K4iIfpciLSsDnSSpMhrNluPnVEoGOklSZdTHJhhZ7fg5lY+BTpJUGfUxW+hUTgY6SVJl1MdabDDQqYQMdJKkymg0W+4SoVIy0EmSKqM+NuE+riolA50kqRIy01muKi0DnSSpEpqtKSYm031cVUoGOklSJdTH2vu4jjiGTiVkoJMkVcL0Pq6OoVMZGegkSZXQaLZb6JzlqjIy0EmSKuFwl6tj6FRCBjpJUiVMBzpnuaqMDHSSpEo4PIbOvVxVQgY6SVIlTI+hs8tVZWSgkyRVgl2uKjMDnSSpEhrNFmtWDrFy2D99Kh9/qyVJlVAfm2C94+dUUgY6SVIl1MdabLC7VSVloJMkVUJ9rOX4OZWWgU6SVAmNZssZriotA50kqRLaY+gMdConA50kqRIaYy1G1jgpQuVkoJMkVUJ9rGULnUrLQCdJKr2pqaQx7ixXlZeBTpJUegfGW2S6S4TKy0AnSSq9p/ZxdQydyslAJ0kqvcP7uDqGTiVloJMkld50oHMdOpWVgU6SVHpPdbka6FROBjpJUunVxyYAx9CpvAx0kqTSaziGTiVnoJMklZ5j6FR2AxXoIuLSiPhKRDwYEdfMcv6VEfFEROzofPzGYu+VJFVXvdkiAtatMtCpnAbmNzsihoFrgUuA3cBtEXFjZt53xKX/nJmvX+K9kqQKqo9NsH7VCoaGot+lSD0xSC10FwEPZuZDmTkOfBi4/BjcK0kqucZYy10iVGqDFOhOB3bNeL67c+xIL42IuyLiYxHxXV3eK0mqoPpYy/FzKrVB+u2erR08j3i+HTgzMxsR8TrgH4BzFnlv+00irgKuAti6deuSi5UkFUej2XKGq0ptkFrodgNnzHi+Bdgz84LMfDIzG53HNwErI2LTYu6d8RrXZea2zNy2efPm5axfkjSg6mMTrkGnUhukQHcbcE5EnB0Rq4A3AzfOvCAinhUR0Xl8Ee36H1/MvZKk6qo3HUOnchuY3+7MbEXE24FPAMPAezPz3oi4unP+3cCPAD8TES3gEPDmzExg1nv78oVIkgZOfazFBgOdSmygfrs73ag3HXHs3TMe/ynwp4u9V5Ik6MxydQydSmyQulwlSVp2E5NTHJqYdAydSs1AJ0kqtQNN93FV+RnoJEml5j6uqgIDnSSp1Ax0qgIDnSSp1OpjEwCOoVOpGegkSaXWcAydKsBAJ0kqNbtcVQUGOklSqdWnW+gMdCoxA50kqdSmx9BtcAydSsxAJ0kqtcZYixVDweoV/slTefnbLUkqtfpYi5E1K4iIfpci9YyBTpJUao1my/FzKj0DnSSp1OpjE4ysdvycys1AJ0kqtfqYLXQqPwOdJKnU6mMtNhjoVHIGOklSqTWaLXeJUOkZ6CRJpVYfm3AfV5WegU6SVFqZ6SxXVYKBTpJUWs3WFBOT6T6uKj0DnSSptOpj7X1cRxxDp5Iz0EmSSmt6H1fH0KnsDHSSpNJqNNstdM5yVdkZ6CRJpXW4y9UxdCo5A50kqbSmA52zXFV2BjpJUmlNj6Hb4Bg6lZyBTpJUWo6hU1UY6CRJpWWXq6rCQCdJKq1Gs8WalUOsHPbPncrN33BJUmm5j6uqwkAnSSqt+ljLXSJUCQY6SVJp1cdarkGnSjDQSZJKq9FsOSFClWCgkySVVn1sgpHVjqFT+RnoJEml9eQhW+hUDQY6SVIpZSZ7D4xz0vpV/S5F6jkDnSSplJ481GJ8corN61f3uxSp5wx0kqRSGm00AdhkoFMFGOgkSaVU6wS6zSMGOpWfgU6SVEo1W+hUIQY6SVIp1erTgc5JESo/A50kqZRqjXGGh4IT1hroVH4GOklSKY3Wm5y4bhVDQ9HvUqSeM9BJkkqp1mg6fk6VYaCTJJVSO9DZ3apqMNBJkkqp1hh3UWFVhoFOklQ6mcloo8km16BTRRjoJEmlU2+2GG+57Zeqw0AnSSqdw2vQjTiGTtWwYqELIuLERbzOVGbuP/pyJEk6erXGOOAuEaqOBQMdsKfzMd9CPsPA1mWpSJKkozRad9svVctiAt39mXnBfBdExJ3LVI8kSUfNfVxVNYsZQ/fSZbpGkqRjotZoMhRw4jrH0KkaFgx0mTl25LGIeNZC10iS1C+1Rnvbr2G3/VJFLHWW603LWkVHRFwaEV+JiAcj4ppZzr8lIu7ufHwhIl4449w3I+LLEbEjIm7vRX2SpGIYrY/b3apKWcwYutks+3/yRMQwcC1wCbAbuC0ibszM+2Zc9g3gezNzX0RcBlwHXDzj/Ksys7bctUmSisV9XFU1S22he8+yVtF2EfBgZj6UmePAh4HLZ16QmV/IzH2dp7cCW3pQhySp4GqNJpvdJUIVsqRAl5l/ttyFAKcDu2Y83905Npd/B3xsZlnAJyPijoi4qgf1SZIKIDM7LXROiFB1LGZh4e2ZeeHRXrMIs3Xj5hzv9yrage7lMw6/LDP3RMTJwM0R8UBm3jLLvVcBVwFs3erSeZJUNgfGJxmbmLLLVZWymDF0z4+Iu+c5H8Dxy1DLbuCMGc+30F7Q+OlvFnEe8JfAZZn5+PTxzNzT+fxYRFxPuwv3GYEuM6+jPfaObdu2zRoYJUnF5aLCqqLFBLrnLeKayaMtBLgNOCcizgYeAd4M/PjMCyJiK/BR4Ccy86szjq8DhjKz3nn8fcDvLENNkqSCObyosGPoVCELBrrM3HksCsnMVkS8HfgE7a3E3puZ90bE1Z3z7wZ+AzgJ+LOIAGhl5jbgFOD6zrEVwAcz8+PHom5J0mCpHW6hcwydqmOpy5b0RGbexBFr3HWC3PTjnwZ+epb7HgJeeORxSVL1TLfQbbbLVRWy1GVLJEkaSKONccJtv1QxXQW6iFjZq0IkSVoOtUaTE9auYsWwbRaqjkV3uUbEXwI/FBEHaM8+vRu4OzP/R6+KkySpW7V60+5WVU43Y+heAZySmRMRcTrtMWvn9aYsSZKWptZosmnE7lZVSzeB7lbgBOCxzHyE9tIiN81/iyRJx9Zoo8mFW0/odxnSMdXNAIPrgM9FxC9ExCsiYjkWE5YkaVnV6uMuKqzK6SbQvR/4G9qtej8LfCEivt6TqiRJWoIDzRaHJiYNdKqcbrpcd2fmb848EBH+i5EkDYzDu0S4qLAqppsWuh0R8Y6ZBzKzucz1SJK0ZG77parqpoXuFOC1EfHLwHbgLmBHZv5tTyqTJKlLo/VxwF0iVD2LDnSZ+W/hcDfrdwEvAC4CDHSSpIHwVJergU7V0vVerp1u1u2dD0mSBsZ0oDvJMXSqGPdFkSSVxmi9yQlrV7LSbb9UMf7GS5JKo9Zo2t2qSupmL9fVwA8DZ828LzN/Z/nLkiSpe7WGiwqrmrppobsBuBxoAQdmfEiSNBDa+7ga6FQ93UyK2JKZl/asEkmSjlKt3nRRYVVSNy10X4iIF/SsEkmSjsKh8UkOjLvtl6qpmxa6lwP/V0R8A2gCAWRmnteTyiRJ6sL0kiUuKqwq6ibQXdazKiRJOkqjh7f9sstV1bPoLtfM3AlsBN7Q+djYOSZJUt+N1qdb6Nb0uRLp2Ft0oIuIdwAfAE7ufLw/Iv7vXhUmSVI3arbQqcK66XL9d8DFmXkAICLeBfwr8D96UZgkSd2o1ccBOGmdY+hUPd3Mcg1gcsbzyc4xSZL6rtZocvxxK1m1wk2QVD3dtND9FfDFiLi+8/yNwP9c9ookSVqC9rZfdreqmhYd6DLzjyLic8DLaLfM/WRm3tmzyiRJ6oL7uKrKummhIzPvAO7oUS2SJC1ZrTHOuadt6HcZUl8sGOgi4vOZ+fKIqAM58xTthYX91yNJ6rtavemiwqqsBQNdZr6883mk9+VIktS9vQfGqTdbnHq8a9CpmrpZh+5dizkmSdKxdufD+wA4/4yN/S1E6pNu5nZfMssxtwOTJPXdnQ/vZ3goeMGW4/tditQXixlD9zPAzwLfERF3zzg1AnyhV4VJkrRY2x/ex/NPHWHtqq7m+kmlsZjf/A8CHwP+K3DNjOP1zNzbk6okSVqkyankrl37+aELt/S7FKlvFjMp4gngCeCKiDgBOAdYAxARZOYtvS1RkqS5ffXbdQ6MT3LhmRv7XYrUN4tum46InwbeAWwBdgAvob2X66t7UpkkSYtw58P7AbjgjBP6W4jUR91MingH8GJgZ2a+CrgAGO1JVZIkLdL2h/dx4rpVnHnS2n6XIvVNN4FuLDPHACJidWY+ADy3N2VJkrQ4dz68jwvO2EhE9LsUqW+6CXS7I2Ij8A/AzRFxA7CnF0VJkrQY+w+O8/XRA1x4pt2tqrZFj6HLzDd1Hv5WRHwG2AB8oidVSZK0CDt27QfgAhcUVsV1s1PEj0bE9PZf3wP8JPBdPalKkqRFuPPh/QwFnGegU8V10+X665lZj4iX09414n3Au3tTliRJC9v+8D6+85QR1q92QWFVWzeBbrLz+QeAd2fmDcCq5S9JkqSFTU0lO3btd/ycRHeB7pGI+Avg3wI3RcTqLu+XJGnZfH20QX2s5fg5ie4C2b+lPQni0szcD5wI/GIvipIkaSHTCwrbQid1N8v1IPDRGc8fBR7tRVGSJC1k+8P7OP64lZx90rp+lyL13YKBLiI+n5kvj4g6kDNPAZmZG3pWnSRJc7jz4f2cf8ZGhoZcUFhaMNBl5ss7n0cWulaSpGPhybEJvvpYnde94NR+lyINhG7WoXvXYo5JktRrd+96gky4YOvGfpciDYRuJkVcMsuxy5arEEmSFmv7w/uIgPMNdBKwuDF0PwP8LPAdEXH3jFMjwBd6VZgkSXO58+F9PGfzejasWdnvUqSBsJhZrh8EPgb8V+CaGcfrmbm3J1VJkjSHzOTOXfv5/nOf1e9SpIGxYJdrZj6Rmd8E3gK8ArgyM3cC6yPioh7XJ0nSYZnJn/7Tg+w/OMGLzz6x3+VIA6ObMXTXAi8Frug8r3eOSZLUc1NTye/843384c1f5U0XnM7l55/W75KkgdFNoLs4M98GjAFk5j6WeS/XiLg0Ir4SEQ9GxDWznI+I+JPO+bsj4sLF3itJKq6JySn+49/exV/9yzf5qZedzR/+6AtZOezuk9K0Re8UAUxExDCdxYUjYjMwtVyFdF77WtqzaXcDt0XEjZl534zLLgPO6XxcDPw5cPEi75UkFdCh8Une9sHt/NMDj/GL3/9cfvaVzybCxYSlmboJdH8CXA+cEhH/GfgR4NeWsZaLgAcz8yGAiPgwcDkwM5RdDvx1ZiZwa0RsjIhTgbMWca8kqQDGW1M8sv8Qu/Ye5OG9B/n77bvZsWs//+VNL+DHL97a7/KkgdTNXq4fiIg7gNd0Dr0xM+9fxlpOB3bNeL6bdivcQtecvsh7Jc2iNTnFo0+MsWvvQXbta/8B3bN/jPHWFBOTU0xOJRNTyeTUFFOdNvkkyZz/daWZsvM/M393pjKZnEpaU0lrMmlNTXFwfJJvPznG1Izfr/WrV3Dtj1/orhDSPBYd6CJiNXAhcHznvh+NCDLzd5apltnaz4/8kzHXNYu5t/0CEVcBVwFs3ep/6ana/r9bd/Kf/vE+mq2nRk+sGApO2bCG41YNs2IoWDEcDA8NsWIoGAqI6X9uMfs/PGk2AcQQBENE0P6g/fu1YihYMTTE8HCwZsUwp59wHFtPXHv44+SR1e7XKi2gmy7XG4AngDuAZg9q2Q2cMeP5FmDPIq9ZtYh7AcjM64DrALZt22YbgyopM/l/P/01/vunvsYrztnED7zgVLaeuJYzTlzLqcevYYWDzSWpULoJdFsy89KeVQK3AedExNnAI8CbgR8/4pobgbd3xshdDDyRmY9GxOgi7pVEe+mH3/7f9/K+f93JD1+4hXf98AsMcJJUcN0Eui9ExAsy88u9KCQzWxHxduATwDDw3sy8NyKu7px/N3AT8DrgQeAg8JPz3duLOqUiG29N8Qt/exc33rWHf/+Ks/mVy55vV5YklUDkIkc2R8R9wHOAb9Ducg0gM/O83pXXW9u2bcvbb7+932VIx8TB8RZXv387t3x1lGsuex5Xf++z+12SJKkLEXFHZm6b7Vw3LXSXLVM9kvrgr/7lm9zy1VHe9cMv4Mde7IQgSSqTbpYt2dnLQiT11pe+sZfnnjJimJOkEupm2ZL/MMvhJ4A7MnPHslUkadlNTSV3PryPHzjPdbwkqYy6mdq2DbiapxbyvQp4JfCeiPil5S9N0nJ5qHaAJ8daXHDGCf0uRZLUA92MoTsJuDAzGwAR8ZvA3wHfQ3ttut9f/vIkLYc7H94HwIVnbuxvIZKknuimhW4rMD7j+QRwZmYeojcLDUtaJtsf3s+GNSv4jk3r+12KJKkHummh+yBwa0Tc0Hn+BuBDEbEOuG/ZK5O0bO58eB/nbz3BNeckqaS6meX6uxFxE/By2mvQXZ2Z04u4vaUXxUk6eo1mi69+u873f9ez+l2KJKlHummhA3iI9k4Ma4C1EfE9mXnL8pclabncvWs/UwkXnumECEkqq26WLflp4B20N77fAbwE+Ffg1T2pTNKy2N6ZEHH+lo39LUSS1DPdTIp4B/BiYGdmvgq4ABjtSVWSls2dD+/nOSev5/i1K/tdiiSpR7oJdGOZOQYQEasz8wHgub0pS9JyyEzu3LWfC87Y2O9SJEk91M0Yut0RsRH4B+DmiNgH7OlFUZKWx87HD7L3wLjj5ySp5LqZ5fqmzsPfiojPAMcDH+9JVZKWxZ272uPnLti6sb+FSJJ6qttZrgBk5ueWuxBJy2/7zv2sX72Cc04e6XcpkqQe6maW6zbgV4EzZ96Xmef1oC5Jy+DOXft44RnHM+yCwpJUat200H0A+EXgy8BUb8qRtFwOjre4/9E6P/O9z+53KZKkHusm0I1m5o09q0TSsvry7ieYnEouPHNjv0uRJPVYN4HuNyPiL4FPA83pg5n50WWvStJRu3PXfgDOP8MZrpJUdt0Eup8Enges5Kku1wQMdNIA2r5zH2dvWseJ61b1uxRJUo91E+hemJkv6FklkpbN9ILCr3jOpn6XIkk6BrrZKeLWiDi3Z5VIWja79x1itN7kAhcUlqRK6KaF7uXAlRHxDdpj6AJIly2RBs/0+Dm3/JKkaugm0F3asyokLau7du1nzcohnvcsFxSWpCroZuuvnb0sRNLy2bP/EFtOWMuK4W5GVUiSisr/t5dKqNZosmm9s1slqSoMdFIJ1RrjbFq/ut9lSJKOkUUHumh7a0T8Ruf51oi4qHelSVqqWr1poJOkCummhe7PgJcCV3Se14Frl70iSUdlbGKSerPF5hEDnSRVRTezXC/OzAsj4k6AzNwXEQ7SkQZMrdHemc8xdJJUHd200E1ExDDt7b6IiM08tQWYpAFRa4wD2OUqSRXSTaD7E+B64JSI+M/A54H/0pOqJC1ZrT7dQmegk6Sq6GYdug9ExB3AazqH3piZ9/emLElLNTrd5eoYOkmqjEUHuoj4D0ccuiwi/g1wR2buWNaqJC3ZUy10jqGTpKropst1G3A1cHrn4yrglcB7IuKXlr80SUtRazTZsGYFq1cM97sUSdIx0s0s15OACzOzARARvwn8HfA9wB3A7y9/eZK6VWuM290qSRXTTQvdVmB8xvMJ4MzMPAQ0l7UqSUs22nBRYUmqmm5a6D4I3BoRN3SevwH4UESsA+5b9sokLUmt0eT5z9rQ7zIkScdQN7NcfzcibgJeDgRwdWbe3jn9ll4UJ6l7tXqTTc9xQoQkVUk3LXQADwHDwBpgbUR8T2besvxlSVqKZmuSJ8dadrlKUsV0s2zJTwPvALYAO4CXAP8KvLonlUnq2uPTu0Q4KUKSKqWbSRHvAF4M7MzMVwEXAKM9qUrSkoy6S4QkVVI3gW4sM8cAImJ1Zj4APLc3ZUlailrDRYUlqYq6GUO3OyI2Av8A3BwR+4A9vShK0tJMB7rNdrlKUqV0M8v1TZ2HvxURnwGOBz7ek6okLUltegydXa6SVCmLCnQREcCWzNwFkJmf62lVkpZktN5kZPUK1qx02y9JqpJFjaHLzKTd1SppgNUaTWe4SlIFdTMp4taIeHHPKpF01GqNphMiJKmCupkU8Srg6oj4JnCA9m4RmZnn9aIwSd2rNcY55+T1/S5DknSMdRPoLutZFZKWRa3R5KXfcVK/y5AkHWPddLk+DLwCuDIzdwIJnNKTqiR1bbw1xf6DE85wlaQK6ibQ/RnwUuCKzvM6cO2yVyRpSR4/0FlUeMQxdJJUNd10uV6cmRdGxJ0AmbkvIvzLIQ2IWt016CSpqrppoZuIiGHaXa1ExGZgajmKiIgTI+LmiPha5/MJs1xzRkR8JiLuj4h7I+IdM879VkQ8EhE7Oh+vW466pCJxlwhJqq5uAt2fANcDJ0fEfwY+D/yXZarjGuDTmXkO8OnO8yO1gP+Ymc8HXgK8LSLOnXH+jzPz/M7HTctUl1QYo9OBzhY6Saqcbrb++kBE3AG8hvaSJW/MzPuXqY7LgVd2Hr8P+Czwy0e8/6PAo53H9Yi4HzgduG+ZapAKbbqFzi5XSaqeRbfQRcQ7gUZmXpuZf7qMYQ7glE5gmw5uJy9Qy1nABcAXZxx+e0TcHRHvna3LViq7Wn2cdauGOW6V235JUtV00+W6AfhERPxzRLwtIrpasiQiPhUR98zycXmXr7Me+Hvg5zPzyc7hPweeDZxPuxXvD+e5/6qIuD0ibh8dHe3mraWB5rZfklRd3XS5/jbw2xFxHvBjwOciYndmvnaR9895XUR8OyJOzcxHI+JU4LE5rltJO8x9IDM/OuO1vz3jmvcA/zhPHdcB1wFs27YtF1O7VATtbb8MdJJURd200E17DPgW8DgLdI124Ubgys7jK4EbjrwgIgL4n8D9mflHR5w7dcbTNwH3LFNdUmGM1t3HVZKqqpsxdD8TEZ+lPQt1E/Dvl3Ef198DLomIrwGXdJ4TEadFxPSM1ZcBPwG8epblSX4/Ir4cEXfT3nP2nctUl1QYttBJUnV1s7DwmbTHre0AiIiXRcS1mfm2oy0iMx+nPXv2yON7gNd1Hn+e9uza2e7/iaOtQSqyickp9rntlyRV1qJb6DLzGoCIeFdEfBP4T8ADPapLUhf2HmjvEuGiwpJUTQu20EXEdwJvpr2H6+PAR4DIzFf1uDZJizRadw06SaqyxXS5PgD8M/CGzHwQDq9JJ2lAPLXtl5MiJKmKFtPl+sO0Z7V+JiLeExHTO0VIGhC1RrvL1RY6SaqmBQNdZl6fmT8GPI/2llzvBE6JiD+PiO/rcX2SFsFtvySp2rqZFHEgMz+Qma8HtgA7gGt6VZikxavVmxy3cph1q7uZuC5JKoulLCxMZu7NzL/IzFcvd0GSujfaaLLJ8XOSVFlLCnSSBouLCktStRnopBKo1ccNdJJUYQY6qQRsoZOkajPQSQXXmpxi78Fxd4mQpAoz0EkFt/fgOJmweb2TIiSpqgx0UsHV6i4qLElVZ6CTCu7wosJ2uUpSZRnopIJzlwhJkoFOKrjR+nSgcwydJFWVgU4quFqjyeoVQ6x32y9JqiwDnVRwtUZ7UeGI6HcpkqQ+MdBJBff4gXG7WyWp4gx0UsHVxybYcNzKfpchSeojA51UcPWxFiNrHD8nSVVmoJMKrjHWckKEJFWcgU4quPrYBCNr7HKVpCoz0EkFNjmVHBiftIVOkirOQCcVWKPZAnAMnSRVnIFOKjADnSQJDHRSodXHJgAcQydJFWegkwqsMdZuoXMMnSRVm4FOKrD6mF2ukiQDnVRodcfQSZIw0EmF5hg6SRIY6KRCcwydJAkMdFKh1cdaDAWsXTXc71IkSX1koJMKrD42wfrVK4iIfpciSeojA51UYPVmy/FzkiQDnVRk9bGWM1wlSQY6qcgaBjpJEgY6qdDqzQlnuEqSDHRSkbVb6BxDJ0lVZ6CTCqw+1mK9Xa6SVHkGOqnA2rNcDXSSVHUGOqmgmq1JxltTjDiGTpIqz0AnFdT0tl+OoZMkGeikgqq7j6skqcNAJxVUozndQmegk6SqM9BJBfXk2ASAs1wlSQY6qaimx9BtcAydJFWegU4qKMfQSZKmGeikgnIMnSRpmoFOKqi6Y+gkSR0GOqmg6s0Wq1YMsXrFcL9LkST1mYFOKqj6WMtdIiRJgIFOKqzGmPu4SpLaBiLQRcSJEXFzRHyt8/mEOa77ZkR8OSJ2RMTt3d4vlUl9bMLxc5IkYEACHXAN8OnMPAf4dOf5XF6Vmedn5rYl3i+VQqPZYmS1a9BJkgYn0F0OvK/z+H3AG4/x/VLh1MdattBJkoDBCXSnZOajAJ3PJ89xXQKfjIg7IuKqJdwvlUbdMXSSpI5j9tcgIj4FPGuWU7/axcu8LDP3RMTJwM0R8UBm3tJlHVcBVwFs3bq1m1ulgVIfm3CWqyQJOIaBLjNfO9e5iPh2RJyamY9GxKnAY3O8xp7O58ci4nrgIuAWYFH3d+69DrgOYNu2bbn0r0jqn8xsj6FzH1dJEoPT5XojcGXn8ZXADUdeEBHrImJk+jHwfcA9i71fKpOD45NMpbtESJLaBiXQ/R5wSUR8Dbik85yIOC0ibupccwrw+Yi4C/gS8H8y8+Pz3S+Vlfu4SpJmGoi/Bpn5OPCaWY7vAV7XefwQ8MJu7pfKanofV7tcJUkwOC10krpQH+u00DkpQpKEgU4qpMOBzi5XSRIGOqmQpsfQOSlCkgQGOqmQHEMnSZrJQCcV0HSX63rH0EmSMNBJhWSgkyTNZKCTCqjRbLFu1TDDQ9HvUiRJA8BAJxVQfWzC8XOSpMMMdFIBNZotZ7hKkg4z0EkFVB9ruQadJOkwA51UQPWxlhMiJEmHGeikAqqPTbDBMXSSpA4DnVRAjaYtdJKkpxjopAJyDJ0kaSYDnVQwk1PJwfFJZ7lKkg4z0EkF0+jsEuE6dJKkaQY6qWDqzQkARhxDJ0nqMNBJBVM/3EJnoJMktRnopIJpNNuBzjF0kqRpBjqpYOpjnS5Xx9BJkjoMdFLBTHe5ug6dJGmagU4qmOlAt8EuV0lSh4FOKhjH0EmSjmSgkwqmPjbB8FBw3MrhfpciSRoQBjqpYBpj7X1cI6LfpUiSBoSBTioY93GVJB3JQCcVTL3ZcoarJOlpDHRSwdTHJtjgGnSSpBkMdFLBNJotZ7hKkp7GQCcVjGPoJElHMtBJBTM9y1WSpGkGOqlg2i10jqGTJD3FQCcVSLM1yfjklF2ukqSnMdBJBTK9j6uBTpI0k4FOKpBGJ9A5hk6SNJOBTiqQp1roHEMnSXqKgU4qkHpzArCFTpL0dAY6qUAcQydJmo2BTiqQhoFOkjQLA51UIPWxdperY+gkSTMZ6KQCaTSd5SpJeiYDnVQg9bEWq1cMsWqF/3QlSU/xr4JUIPVmy/FzkqRnMNBJBeI+rpKk2RjopAJpjE04fk6S9AwGOqlA2i10BjpJ0tMZ6KQCaTRbttBJkp7BQCcVyN4D45ywdlW/y5AkDRgDnVQQU1PJ4wfG2TRioJMkPZ2BTiqIfQfHmZxKNq1f3e9SJEkDxkAnFUStMQ5goJMkPYOBTiqIWqMJwOYRA50k6ekGItBFxIkRcXNEfK3z+YRZrnluROyY8fFkRPx859xvRcQjM8697ph/EVKPTQc6W+gkSUcaiEAHXAN8OjPPAT7def40mfmVzDw/M88HXgQcBK6fcckfT5/PzJuORdHSsTRa77TQGegkSUcYlEB3OfC+zuP3AW9c4PrXAF/PzJ29LEoaJKONJquGh9hwnOvQSZKeblAC3SmZ+ShA5/PJC1z/ZuBDRxx7e0TcHRHvna3LViq6Wn2ck9avIiL6XYokacAcs0AXEZ+KiHtm+bi8y9dZBfwg8LczDv858GzgfOBR4A/nuf+qiLg9Im4fHR3t/guR+qTWaDp+TpI0q2PWd5OZr53rXER8OyJOzcxHI+JU4LF5XuoyYHtmfnvGax9+HBHvAf5xnjquA64D2LZtW3bxJUh9VWs0OWXDmn6XIUkaQIPS5XojcGXn8ZXADfNcewVHdLd2QuC0NwH3LGt10gBot9C5S4Qk6ZkGJdD9HnBJRHwNuKTznIg4LSIOz1iNiLWd8x894v7fj4gvR8TdwKuAdx6bsqVjY2oqebwxbperJGlWAzFdLjMfpz1z9cjje4DXzXh+EDhplut+oqcFSn22/9AELbf9kiTNYVBa6CTN4/Ciwu4SIUmahYFOKoBafXqXCMfQSZKeyUAnFcBop4XuZFvoJEmzMNBJBVBrjAPu4ypJmp2BTiqAWqPJyuHg+ONW9rsUSdIAMtBJBTBab3LSutVu+yVJmpWBTiqAWqPJphEnREiSZmegkwrAfVwlSfMx0EkFUKuPs9lAJ0mag4FOGnBTU8njB5ouKixJmpOBThpwTxyaYGLSbb8kSXMz0EkD7vC2X+4SIUmag4FOGnDTu0Q4hk6SNBcDnTTgDu8S4Rg6SdIcDHTSgKvVbaGTJM3PQCcNuFqjyYoht/2SJM3NQCcNuNF6k5PWr2JoyG2/JEmzM9BJA85dIiRJCzHQSQOu1hg30EmS5mWgkwacLXSSpIUY6KQBlpk83hhns0uWSJLmYaCTBtiTh1qMT065S4QkaV4GOmmAjTbGAGyhkyTNy0AnDbDRemeXCMfQSZLmYaCTBlits4+rgU6SNB8DnTTApgOdXa6SpPkY6KQBVms0GR4KNrrtlyRpHgY6aYDV6uOctM5tvyRJ8zPQSQNs1EWFJUmLYKCTBlit0WST4+ckSQsw0EkDrFZvuqiwJGlBBjppQGUmNbf9kiQtgoFOGlBPjrW3/drsGDpJ0gIMdNKAclFhSdJiGeikATVaN9BJkhbHQCcNqMMtdCNOipAkzc9AJw2omi10kqRFMtBJA6rWGGd4KDhhrS10kqT5GeikAVVrNDlx3SqG3fZLkrQAA500oGpu+yVJWiQDnTSgRt0lQpK0SAY6aUDVGuMuKixJWhQDnTSAnhyb4LH6GKccv6bfpUiSCsBAJw2g67c/wsRkctl3P6vfpUiSCsBAJw2YzOT9t+7khVuO57wtG/tdjiSpAAx00oD54jf28rXHGrzlJWf2uxRJUkEY6KQB8/5bd7JhzQrecN5p/S5FklQQBjppgDxWH+MT936LH3nRGRy3arjf5UiSCsJAJw2Qv7ltFxOTyVtesrXfpUiSCsRAJw2IyankQ1/axcuecxLP3ry+3+VIkgrEQCcNiM888BiP7D/EWy92MoQkqTsGOmlAvP+LOzl5ZDWvPfeUfpciSSoYA500AB5+/CCf++ooV1y0lZXD/rOUJHVnIP5yRMSPRsS9ETEVEdvmue7SiPhKRDwYEdfMOH5iRNwcEV/rfD7h2FQuLY8PfGknQxFccZGTISRJ3RuIQAfcA/wQcMtcF0TEMHAtcBlwLnBFRJzbOX0N8OnMPAf4dOe5NPCeODTB39y2i4/ctovXPv9knuXerZKkJVjR7wIAMvN+gIiY77KLgAcz86HOtR8GLgfu63x+Zee69wGfBX65N9VKR2dsYpJ/euAxbtjxCJ95YJTxySnOOmktP//a7+x3aZKkghqIQLdIpwO7ZjzfDVzceXxKZj4KkJmPRsTJx7q42Xzmgcf4u+27+12GeiFnO5RMTXU+J2TC5NQUB8cnOTg+yYFmiwPjLfYfnKDZmmLzyGre+pIzufz80zhvy/EL/QeNJElzOmaBLiI+BTxrllO/mpk3LOYlZjk2y5/VBeu4CrgKYOvW3o5XevzAOA88+mRP30P9M1sAGwoY6hwfimBoCNauWsHmkdWcedJa1q1awYbjVvC933kyL332SQwPGeIkSUfvmAW6zHztUb7EbuCMGc+3AHs6j78dEad2WudOBR6bp47rgOsAtm3b1nUg7MaPvGgLP/KiLb18C0mSpIGZFLEYtwHnRMTZEbEKeDNwY+fcjcCVncdXAotp8ZMkSSqFgQh0EfGmiNgNvBT4PxHxic7x0yLiJoDMbAFvBz4B3A/8TWbe23mJ3wMuiYivAZd0nkuSJFVCZPa013Ggbdu2LW+//fZ+lyFJkrSgiLgjM2ddr3cgWugkSZK0dAY6SZKkgjPQSZIkFZyBTpIkqeAMdJIkSQVnoJMkSSo4A50kSVLBGegkSZIKzkAnSZJUcAY6SZKkgjPQSZIkFZyBTpIkqeAMdJIkSQVnoJMkSSo4A50kSVLBGegkSZIKzkAnSZJUcAY6SZKkgjPQSZIkFZyBTpIkqeAiM/tdQ99ExCiws991FMgmoNbvIvQ0/kwGkz+XwePPZDD5c+nOmZm5ebYTlQ506k5E3J6Z2/pdh57iz2Qw+XMZPP5MBpM/l+Vjl6skSVLBGegkSZIKzkCnblzX7wL0DP5MBpM/l8Hjz2Qw+XNZJo6hkyRJKjhb6CRJkgrOQKcliYhfiIiMiE39rqXqIuIPIuKBiLg7Iq6PiI39rqmqIuLSiPhKRDwYEdf0ux5BRJwREZ+JiPsj4t6IeEe/a1JbRAxHxJ0R8Y/9rqUMDHTqWkScAVwCPNzvWgTAzcB3Z+Z5wFeBX+lzPZUUEcPAtcBlwLnAFRFxbn+rEtAC/mNmPh94CfA2fy4D4x3A/f0uoiwMdFqKPwZ+CXAA5gDIzE9mZqvz9FZgSz/rqbCLgAcz86HMHAc+DFze55oqLzMfzcztncd12gHi9P5WpYjYAvwA8Jf9rqUsDHTqSkT8IPBIZt7V71o0q58CPtbvIirqdGDXjOe7MTgMlIg4C7gA+GKfSxH8d9oNA1N9rqM0VvS7AA2eiPgU8KxZTv0q8P8A33dsK9J8P5PMvKFzza/S7l76wLGsTYfFLMdsxR4QEbEe+Hvg5zPzyX7XU2UR8Xrgscy8IyJe2edySsNAp2fIzNfOdjwiXgCcDdwVEdDu2tseERdl5reOYYmVM9fPZFpEXAm8HnhNuhZRv+wGzpjxfAuwp0+1aIaIWEk7zH0gMz/a73rEy4AfjIjXAWuADRHx/sx8a5/rKjTXodOSRcQ3gW2Z6cbKfRQRlwJ/BHxvZo72u56qiogVtCelvAZ4BLgN+PHMvLevhVVctP/r833A3sz8+T6XoyN0Wuh+ITNf3+dSCs8xdFLx/SkwAtwcETsi4t39LqiKOhNT3g58gvbA+78xzA2ElwE/Aby68+9jR6dlSCoVW+gkSZIKzhY6SZKkgjPQSZIkFZyBTpIkqeAMdJIkSQVnoJMkSSo4A50kSVLBGegkSZIKzkAnaUkiYnLGQq07IuKaea7dGBE/e4zqasxx/Oci4v6I+EDn+ReO4rV+KyJ+YRH3H8uv+2nvFRFnRcQ9y/TaZ0XEoYjYscB1n42I7z/i2M9HxJ9FxHGd35PxiNi0HHVJeoqBTtJSHcrM82d8/N48124Eugo20bac/x/1s8DrMvMtAJn5b5bxteeykTm+7h58fXO+1zL5emaev8A1HwLefMSxNwMfysxDnfvd31bqAQOdpGUTES+OiLsjYk1ErIuIeyPiu4HfA57daaH5g861b42IL3WO/UVEDHdagu6PiD8DtgOv6Dx/T+e1PhkRx3Xu/4eIuKNz/KoF6no38B3AjRHxzs6xxozzz6hlltf41Yj4SkR8CnjuLOfXRcT/iYi7IuKeiPixI7/uWb6+Mxb4Pjzj6+68169HxAMRcXNEfKjTWviM7zEwPNdrHFH7ZyPiuZ3HJy2mZW+O79nfAa+PiNWda84CTgM+v9DrSTo6BjpJSzXdhTb98WOZeRtwI/CfgN8H3p+Z9wDX0GnhycxfjIjnAz8GvKzTajMJvKXzus8F/jozLwB2AucA12bmdwH7gR/uXPdTmfkiYBvwcxFx0lyFZubVtFuGXpWZfzzz3AK1TF/zItotTRcAPwS8eJa3uRTYk5kvzMzvBj5+5Nc9y9e3dp73nvXrjohtncfTtWzrXD/be831vTvSc4CvdR6fB3x5juumvx+zfs8y83HgS53vBbS/Zx9J95iUem5FvwuQVFiH5uiC+x3gNmAM+Lk57n0N8CLgtogAOA54DLgF2JmZt8649huZuaPz+A7grM7jn4uIN3Uen0E7vDy+hK9jrlpmegVwfWYeBIiIG2d5nS8D/y0i3gX8Y2b+c0ScMMt1M7+++b4Pc33dLwduyMxDnVr+9zxf21yvcVhEnAk8kplTnUPnAXfP85rz1Q1Pdbve0Pn8Uwu8lqRlYKCTtNxOBNYDK4E1wIFZrgngfZn5K0872O6iO/L65ozHk7RbBl8JvBZ4aWYejIjPdt5rKWatZRbztjJl5lc7LXmvA/5rRHwS+OtZLp359c33fXjG1z3jnsWa6zVmOp+nB7gXAR9Z4HXn+579A/BHEXEhcFxmbl90tZKWzC5XScvtOuDXgQ8A7+ocqwMjM675NPAjEXEyQESc2GkpWqzjgX2dMPc84CVHUe9iarkFeFO0Z2qOAG848kUi4jTgYGa+H/hvwIU88+teynsf6fPAGzrjFNcDP9A5vtB7zeWFdMJwRJwDXM4CXa7z1Z2ZDeCzwHtpt9ZJOgZsoZO0VMfF05ex+DhwH9DKzA92Bsl/ISJenZn/FBH/0hls/7HOOLpfAz4Z7ZmeE8DbgG8t8r0/DlwdEXcDXwFuXeD6OWXmfXPUsnPGNdsj4iPAjs7xf57lpV4A/EFETHVe42cy8/GZXzdw7SLfe87vQ2be1unyvatTy+3AEwu91zzOBw5FxF20W+ruB64EfneeGhb6nn0I+CjPnPEqqUfCsaqSVCwRsT4zGxGxlnbr4VVL7dqMiAeBCzKzPs81Z9EeF/jdSyr46a/1TWBbZtaO9rUkPcUuV0kqnus6raPbgb8/ijA3AkzNF+Y6JoHjY4GFhRd4r+kW3ZXA1AKXS+qSLXSSJEkFZwudJElSwRnoJEmSCs5AJ0mSVHAGOkmSpIIz0EmSJBWcgU6SJKngDHSSJEkFZ6CTJEkquP8faBWEoblmv20AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Parameters\n",
    "M = 10      # Size of the spin chain \n",
    "J = 1      # Exchange interaction\n",
    "kT = 0.1   # Temperature\n",
    "\n",
    "# Prepare plots\n",
    "df = pd.DataFrame()\n",
    "\n",
    "# Brute-Force (green)\n",
    "for mu in np.linspace(-5,5,100):\n",
    "    bf = brute_force(M,J,mu,[kT])\n",
    "    df = pd.concat([df,bf])\n",
    "\n",
    "plt.figure(figsize=(10,10));\n",
    "plt.plot(df.mu,df.m);\n",
    "plt.xlabel('External field strength $\\mu$ [eV]');\n",
    "plt.ylabel('Average magnetisation $m$ [-]');\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transfer Matrix\n",
    "\n",
    "The Transfer Matrix approach exploits a *duality* between the 1D Ising Model and a 0D quantum mechanical problem. The problem is recast into an eigenvalue problem of a $2\\times 2$ matrix, which is very efficient to solve (even analytically as we will see below). \n",
    "\n",
    "Lets start with the usual partition function\n",
    "\n",
    "$$\n",
    "Z= \\sum_{\\{\\sigma\\}} \\exp(-\\beta H_i) = \\sum_{\\{\\sigma\\}} \\exp\\left(-\\beta \\sum_j(J\\sigma_j\\sigma_{j+1} +\\mu\\sigma_j) \\right)=\\sum_{\\{\\sigma\\}}\\prod_j\\exp(-\\beta(\\underbrace{J\\sigma_j\\sigma_{j+1} +\\mu\\sigma_j}_{=H_{j,j+1}}))\n",
    "$$\n",
    "\n",
    "Now we switch to the matrix representation that is described earlier by letting\n",
    "\n",
    "$$\n",
    "T_{j,j+1} = \\exp(-\\beta\\cdot H_{j,j+1}) = \\rightarrow \\langle j|T| j+1\\rangle = \\langle j|\\exp(-\\beta H)| j+1\\rangle\n",
    "$$\n",
    "\n",
    "With this definition, the partition function becomes\n",
    "\n",
    "$$\n",
    "Z = \\sum_{\\{\\sigma\\}}\\langle 1|T| 2\\rangle\\langle 2|T| 3\\rangle\\dots\\langle N|T|1\\rangle\n",
    "$$\n",
    "\n",
    "where the last term exploits periodic boundary conditions. Now we recall that $|i\\rangle$ really only is a vector and $\\langle i|$ is its (conjugate) transpose. So $|i\\rangle\\langle i|$ is the *outer* product (or tensor product) of a unit vector. In quantum mechanics, it is called the projection operator. An important property of the projection operator is that \n",
    "\n",
    "$$\n",
    "I = \\sum_n |n\\rangle\\langle n|\n",
    "$$\n",
    "\n",
    "the sum over all projections is the identity. This can be understood as taking a vector and projecting out all its components and then summing up the components, which restores the vector. \n",
    "\n",
    "Because the sum of the partition function is over all spin configurations, it contains sums over all projection operators apart from $|1\\rangle\\langle 1|$. This greatly simplifies the sum over products, as can be demonstrated with a simple two site chain:\n",
    "\n",
    "$$\n",
    "\\begin{array}{rcll}\n",
    "Z & = & \\langle1_+|T|2_+\\rangle\\langle2_+|T|1_+\\rangle & (+,+)\\\\\n",
    "  & + & \\langle1_+|T|2_-\\rangle\\langle2_-|T|1_+\\rangle & (+,-)\\\\\n",
    "  & + & \\langle1_-|T|2_+\\rangle\\langle2_+|T|1_-\\rangle & (-,+)\\\\\n",
    "  & + & \\langle1_-|T|2_-\\rangle\\langle2_-|T|1_-\\rangle & (-,-)\\\\\n",
    "  & = & \\langle1_+|T|\\left[\\underbrace{|2_+\\rangle\\langle2_+|+|2_-\\rangle\\langle2_-|}_{=I}\\right]|T|1_+\\rangle\\\\\n",
    "  & + & \\langle1_-|T|\\left[\\underbrace{|2_+\\rangle\\langle2_+|+|2_-\\rangle\\langle2_-|}_{=I}\\right]|T|1_-\\rangle\\\\\n",
    "  & = & \\langle1_+|T^2|1_+\\rangle + \\langle1_-|T^2|1_-\\rangle\\\\\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "In general you have\n",
    "\n",
    "$$\n",
    "\\begin{array}{rcl}\\displaystyle\n",
    "Z & = & \\sum_{\\{\\sigma\\}}\\langle 1|T| 2\\rangle\\langle 2|T| 3\\rangle\\dots\\langle N|T|1\\rangle\\\\\n",
    "  & = & \\sum_{\\sigma = \\pm 1}\\langle1|T^N|1\\rangle\\\\\n",
    "  & = & {\\rm Tr}({\\bf T^N})\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "The partition function can, therefore, be computed from the trace of the N-th power of the transfer matrix through this interesting duality with quantum mechanics.\n",
    "\n",
    "Taking the $N$-th power of a matrix, however, is still an expensive operation if $N$ is large. This can be avoided by rotating the basis states so that ${\\bf T}$ becomes diagonal. In this basis, the eigenvalues $\\lambda_i$ of ${\\bf T}$ occupy the main diagonal and the trace becomes:\n",
    "\n",
    "$$\n",
    "{\\rm Tr}({\\bf T}^N) = \\sum \\lambda_i^N\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def transfer_matrix(M,J,mu,kTs):\n",
    "    '''\n",
    "    Solve the one-dimensional Ising model by computing the partition function through the \n",
    "    transfer matrix formalism.\n",
    "    \n",
    "    The function returns the results as a Pandas dataframe containing the free energy g,\n",
    "    the internal energy u, and the magnetisation m\n",
    "    \n",
    "    Parameters:\n",
    "        \n",
    "        M   - length of spin chain\n",
    "        J   - exchange interaction\n",
    "        mu  - external field (same units as J)\n",
    "        kTs - array with temperatures (in units of J)\n",
    "    \n",
    "    Returns:\n",
    "    \n",
    "        A Pandas dataframe with free energy g, internal energy u, magnetisation m\n",
    "        \n",
    "    Usage: \n",
    "        \n",
    "        transfer_matrix(10,-1,0.1,np.linspace(0.1,4,30))\n",
    "    '''\n",
    "\n",
    "    # Create a new data-frame to take results\n",
    "    df = pd.DataFrame(columns=['kT','g','u','m','mu'],dtype=object)\n",
    "    \n",
    "    # Iterate over all tempertures\n",
    "    for i,kT in enumerate(kTs):\n",
    "        # Initialise variables\n",
    "        Z = m = u = 0\n",
    "        beta = 1/kT\n",
    "\n",
    "        # Compute the Hamiltonian and Transfer Matrix\n",
    "        H = hamiltonian(J,mu,1)\n",
    "        T = np.exp(-beta*H)\n",
    "        \n",
    "        # Compute eigenvalues of T and compute its eigenvalues\n",
    "        v,l = np.linalg.eig(T)\n",
    "        L = np.matrix(np.diag(v)**M)\n",
    "        Z = np.sum(v**M)\n",
    "\n",
    "        # Rotate the the Pauli matrix into the eigenbasis of T\n",
    "        U  = np.matrix(l)\n",
    "        Ui = np.matrix(np.linalg.inv(U))\n",
    "        Mg = Ui*magnetisation(1)*U\n",
    "                \n",
    "        # Compute thermodynamic quantities\n",
    "        m = np.trace(L*Mg)/Z\n",
    "        \n",
    "        g = -1/beta*np.log(Z)/M\n",
    "        \n",
    "         # Append result to dataframe\n",
    "        s = pd.DataFrame({'g': [g], 'kT': [kT], 'm': [m], 'mu': [mu], 'J':[J], 'M':[m]})\n",
    "        df = pd.concat([df,s])\n",
    "    \n",
    "    # The internal energy is not easily obtained in the transfer matrix formalism,\n",
    "    # we therefore use U = -dZ/dbeta to calculate the internal energy numerically\n",
    "    if len(kTs) > 1:\n",
    "        lZ = (1/kTs)*df.g\n",
    "        df['u'] = np.gradient(lZ,1/kTs)\n",
    "    \n",
    "    return df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2: Magnetisation of odd AFM chains\n",
    "\n",
    "The Transfer Matrix method is computationally more efficient, which allows to study longer chains. The below cell calculates the magnetisation for a longer chain with an odd number of spins."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['kT', 'g', 'u', 'm', 'mu', 'J', 'M']\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnQAAAJRCAYAAADMCCpKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA9A0lEQVR4nO3de5jkZX3n/fe3e7rn2DMjMAcERtCgBhFRB5CARlRcMLro5qDGdYnG8JjoxmQ32bCPu0n2lESzmzxXEoyLhmjWYzYRYbOAp2iUEJQBEYaDkXCQYaanhwGmama6+vh9/qjqoenpQ1VPV1f9qt6v6yq7fqeqb3UP9qfv+/7dd2QmkiRJKq6eVhcgSZKkY2OgkyRJKjgDnSRJUsEZ6CRJkgrOQCdJklRwBjpJkqSCW9HqAlrphBNOyFNPPbXVZUiSJC3o9ttvfzwzN812rKsD3amnnsqOHTtaXYYkSdKCIuKRuY7Z5SpJklRwBjpJkqSCM9BJkiQVnIFOkiSp4Ax0kiRJBWegkyRJKjgDnSRJUsEZ6CRJkgrOQCdJklRwBjpJkqSCM9BJkiQVnIFOkiSp4Ax0kiRJBWegkyRJKjgDnSRJUsEZ6CRJkgrOQCdJklRwBjpJkqSCM9BJkiQVXFsFuoi4JiKGImLnHMcjIv4oIh6IiLsi4mXTjl0SEd+vHbty+aqWJElqrbYKdMAngEvmOX4pcHrtcQXwpwAR0QtcVTt+BvD2iDijqZVKkiS1ibYKdJn5TeCJeU65DPiLrLoV2BgRJwLnAg9k5oOZOQp8rnauJElSx1vR6gIadBLw6LTtXbV9s+0/bxnrkiR1uf0HR/jyvXu5b0+JQyMTDI+NV7+OTnB4bJyJSchMMmEyk6S6rc7wH954Bhe9YHPL3r9ogS5m2Zfz7D/6BSKuoNpdy7Zt25auMklS19l/cISb7hnkhrv3cOuDTzAxmaxftYKBVX2s6e+tPVaweWAVvT1BABHQE0EEBDH7bzAVzvpVfS19/6IFul3AKdO2TwZ2A/1z7D9KZl4NXA2wfft2/zSSJDVsdHySf/3ZO/jKvXuZTDjthLW898efyxtefCJnnLieCFOallfRAt31wPsj4nNUu1QPZOaeiNgHnB4RpwGPAW8DfraFdUqSOtjf3LWbL92zl3ddcCo/s/0UXrh1wBCnlmqrQBcRnwVeDZwQEbuA3wL6ADLzo8ANwBuAB4DDwLtqx8Yj4v3Al4Be4JrMvGfZP4AkqeNlJn9280P8yOZ1/OYbzzDIqS20VaDLzLcvcDyB981x7AaqgU+SpKa57eEnuWd3id95y4sNc2obbTVtiSRJ7e6amx9i45o+3vLSk1pdinSEgU6SpDo9+sRhvnzvID977jZW9/e2uhzpCAOdJEl1+sQtD9MTwb86/9RWlyI9g4FOkqQ6lCtjfP62R/mJs05k64ZVrS5HegYDnSRJdfir23dxcGScd11wWqtLkY5ioJMkaQETk8knbnmYlz/nWZx9ysZWlyMdxUAnSdIC/vb+IR7Zf5h32zqnNmWgkyRpAdfc/BDP3rCKf/aiLa0uRZqVgU6SpHncu7vEPzy4n8t/7FRW9PprU+3Jf5mSJM3jb+7aTV9v8LZztrW6FGlOBjpJkuYxeKDC5oFVbFjT1+pSpDkZ6CRJmsdQeYTN61e2ugxpXgY6SZLmMVSusHnAQKf2ZqCTJGke+8ojbDLQqc0Z6CRJmsPo+CRPHh5j84BLfam9GegkSZrDvoMjAHa5qu0Z6CRJmsNQqQJgl6vanoFOkqQ57CtPtdDZ5ar2ZqCTJGkOQ1OBzmlL1OYMdJIkzWGoPEIEHL+2v9WlSPMy0EmSNId95RGOX9vvGq5qe/4LlSRpDvvKFTY5fk4FYKCTJGkOQ04qrIIw0EmSNId95RHnoFMhGOgkSZrF5GQa6FQYBjpJkmbx5OFRxifTLlcVgoFOkqRZPL3slzdFqP0Z6CRJmsVQyUmFVRwGOkmSZjG1SsSmdQY6tT8DnSRJs9jnsl8qEAOdJEmzGCpXWLdyBWv6V7S6FGlBBjpJkmbhpMIqEgOdJEmz2GegU4EY6CRJmoWTCqtIDHSSJM1iqFSxhU6FYaCTJGmGQyPjHBqdcFJhFYaBTpKkGY5MWWILnQrCQCdJ0gxHJhU20KkgDHSSJM3gpMIqGgOdJEkzDJUrAI6hU2EY6CRJmmGoPMKKnmDj6r5WlyLVxUAnSdIMU5MK9/REq0uR6mKgkyRpBpf9UtEY6CRJmmGoVHHKEhWKgU6SpBkePzjCJm+IUIEY6CRJmmZ8YpL9h0btclWhGOgkSZrm8YOjZLpKhIrFQCdJ0jQu+6UiMtBJkjTN1KTCdrmqSAx0kiRNM3Rk2S9vilBxGOgkSZpmqsv1hHX9La5Eqp+BTpKkaYbKFTau6WPlit5WlyLVzUAnSdI0Q6URb4hQ4bRVoIuISyLi+xHxQERcOcvxX4+IO2uPnRExERHH1Y49HBF3147tWP7qJUmdYN/BETY7qbAKpm0CXUT0AlcBlwJnAG+PiDOmn5OZv5+ZZ2fm2cC/B/4uM5+YdspFtePbl6tuSVJnGSq5jquKp20CHXAu8EBmPpiZo8DngMvmOf/twGeXpTJJUlfITPaV7XJV8bRToDsJeHTa9q7avqNExBrgEuCvp+1O4MsRcXtEXNG0KiVJHas0PM7oxKQtdCqcFa0uYJqYZV/Oce6bgL+f0d16QWbujojNwFci4v7M/OZRb1INe1cAbNu27VhrliR1ECcVVlG1UwvdLuCUadsnA7vnOPdtzOhuzczdta9DwLVUu3CPkplXZ+b2zNy+adOmYy5aktQ5jkwq7E0RKph2CnS3AadHxGkR0U81tF0/86SI2AD8OHDdtH1rI2Jg6jnwemDnslQtSeoYR9ZxXW8LnYqlbbpcM3M8It4PfAnoBa7JzHsi4r214x+tnfoW4MuZeWja5VuAayMCqp/pM5l50/JVL0nqBHa5qqjaJtABZOYNwA0z9n10xvYngE/M2Pcg8JImlydJ6nBDpRFW9fUwsLKtfj1KC2qnLldJklpqalLhWo+PVBgGOkmSapxUWEVloJMkqWaoXHFSYRWSgU6SpBpXiVBRGegkSQLGJyYpVcZ51tr+VpciNcxAJ0kSUK6MA7B+VV+LK5EaZ6CTJIlpgW61gU7FY6CTJAkoVcYAGFjlHHQqHgOdJEk8HejsclURGegkSQJKw9UuV1voVEQGOkmSgHKthW6DY+hUQAY6SZKAUsUWOhWXgU6SJJ5uoVu30kCn4jHQSZJEddqStf29rOj1V6OKx3+1kiQBpeEx56BTYRnoJEmi2kLn+DkVlYFOkiSq89ANOAedCspAJ0kS1Ra69bbQqaAMdJIkYQudis1AJ0kStRa61bbQqZgMdJKkrpeZlIZtoVNxGegkSV2vMjbJ+GSy3kCngjLQSZK6Xqm2SoTTlqioDHSSpK43teyXEwurqAx0kqSuV6qMA7bQqbgMdJKkrlcarrXQOYZOBWWgkyR1vXKthc6JhVVUBjpJUtd7+qYIW+hUTAY6SVLXO9JC58TCKigDnSSp65WGx+jtCVb39ba6FGlRDHSSpK5XroyzftUKIqLVpUiLYqCTJHW9UsVlv1RsBjpJUtcrV8YdP6dCM9BJkrpeaXiMgZW20Km4DHSSpK5nC52KzkAnSep6ZcfQqeAMdJKkrleqjLvslwrNQCdJ6moTk8nBkXEGXPZLBWagkyR1tYO1VSIMdCoyA50kqatNreO6frVdriouA50kqasdCXS20KnADHSSpK5WrnW5elOEisxAJ0nqaqXhagud05aoyAx0kqSudqSFzomFVWAGOklSV5saQ2cLnYrMQCdJ6mplpy1RBzDQSZK6Wrkyxuq+Xvp6/ZWo4vJfrySpq5WGXSVCxWegkyR1tfLImJMKq/AMdJKkrmYLnTqBgU6S1NXKlTEnFVbhGegkSV2tVLGFTsXXVoEuIi6JiO9HxAMRceUsx18dEQci4s7a4zfrvVaSpNmUK46hU/G1zZ8kEdELXAVcDOwCbouI6zPz3hmnfisz37jIayVJegbH0KkTtFML3bnAA5n5YGaOAp8DLluGayVJXaoyNsHoxKRj6FR47RToTgIenba9q7ZvpvMj4nsRcWNEvKjBayVJOmJq2a/1ttCp4NrpX3DMsi9nbN8BPCczD0bEG4AvAqfXeW31TSKuAK4A2LZt26KLlSQV39SyX46hU9G1UwvdLuCUadsnA7unn5CZpcw8WHt+A9AXESfUc+2017g6M7dn5vZNmzYtZf2SpIJxHVd1inYKdLcBp0fEaRHRD7wNuH76CRGxNSKi9vxcqvXvr+daSZJmKg1Xu1wHHEOngmubP0kyczwi3g98CegFrsnMeyLivbXjHwV+CvjFiBgHhoG3ZWYCs17bkg8iSSqMI12uBjoVXNsEOjjSjXrDjH0fnfb8T4A/qfdaSZLmM3VThF2uKrp26nKVJGlZlafucvWmCBWcgU6S1LVKw+P0BKzt7211KdIxMdBJkrpWuTLGwKo+avfbSYVloJMkda1SxWW/1BkMdJKkrlWujHmHqzqCgU6S1LVKw7bQqTMY6CRJXatUGfMOV3UEA50kqWuVHUOnDmGgkyR1rZJj6NQhDHSSpK40OZkcHBlnvS106gAGOklSVzo4Ok4mDNhCpw5goJMkdaVyZRyA9attoVPxGegkSV2pNFxdx9UWOnUCA50kqSsdaaEz0KkDGOgkSV3p6RY6u1xVfAY6SVJXKo9UA50TC6sTGOgkSV2pNFztcrWFTp3AQCdJ6krlil2u6hwGOklSVypXxlm5ooeVK3pbXYp0zAx0kqSuVKqMOWWJOoaBTpLUlUqVcScVVscw0EmSulJp2BY6dQ4DnSSpK5Ur46z3hgh1CAOdJKkrlSpjrhKhjmGgkyR1pbJj6NRBDHSSpK7kGDp1EgOdJKnrjIxPMDI+6Rg6dQwDnSSp65QrU8t+2UKnzmCgkyR1nalA5xg6dQoDnSSp6xxZx3WlLXTqDAY6SVLXKQ1PdbnaQqfOYKCTJHWdqRa69attoVNnMNBJkrpOaarL1RY6dQgDnSSp6zx1uBroNthCpw5hoJMkdZ29pRHW9PeybqUtdOoMBjpJUtfZW6qwdf0qIqLVpUhLwkAnSeo6e0sVNq9f2eoypCVjoJMkdZ3BWgud1CkMdJKkrpKZDJVG2LLBQKfOYaCTJHWVJw+PMToxyZYBA506h4FOktRVBg9UANhqC506iIFOktRV9paqgW6LY+jUQQx0kqSuMliyhU6dx0AnSeoqUy10m9Y5bYk6h4FOktRV9pYqnLCun/4V/gpU5/BfsySpqwweqDh+Th3HQCdJ6iqDpREnFVbHMdBJkrrKUKnCZgOdOsyKhU6IiOPqeJ3JzHzq2MuRJKl5RsYn2H9o1BY6dZwFAx2wu/aIec7pBbYtSUWSJDXJUGkEgK0bvMNVnaWeQHdfZr50vhMi4rtLVI8kSU0zVK5OWWKXqzpNPWPozl+icyRJaqnBA7UWOgOdOsyCgS4zKzP3RcTWhc6RJKndHFklwkCnDrPYu1xvWNIqaiLikoj4fkQ8EBFXznL8HRFxV+1xS0S8ZNqxhyPi7oi4MyJ2NKM+SVKx7S1V6F/Rw8Y1fa0uRVpS9Yyhm818N0gs7gUjeoGrgIuBXcBtEXF9Zt477bSHgB/PzCcj4lLgauC8accvyszHl7o2SVJn2FuqsGX9SiKW/NeY1FKLbaH72JJWUXUu8EBmPpiZo8DngMumn5CZt2Tmk7XNW4GTm1CHJKlDDR6o2N2qjrSoQJeZH1nqQoCTgEenbe+q7ZvLzwM3Ti8L+HJE3B4RVzShPklSwVVb6Ax06jwLBrqIuGMpzqnDbO3fOcf7XUQ10P3GtN0XZObLgEuB90XEq+a49oqI2BERO/bt23esNUuSCiIzGTTQqUPVM4buRyPirnmOB7BhCWrZBZwybftkqhMaP/PNIs4CPg5cmpn7p/Zn5u7a16GIuJZqF+43Z16fmVdTHXvH9u3bZw2MkqTOU6qMUxmbtMtVHameQPfCOs6ZONZCgNuA0yPiNOAx4G3Az04/ISK2AV8A3pmZ/zht/1qgJzPLteevB/7zEtQkSeoQe2tTlmzZYKBT51kw0GXmI8tRSGaOR8T7gS9RXUrsmsy8JyLeWzv+UeA3geOBj9TuUBrPzO3AFuDa2r4VwGcy86blqFuSVAyDB5yDTp1rsdOWNEVm3sCMOe5qQW7q+XuA98xy3YPAS2bulyRpypEWuvWu46rOs9hpSyRJKpSnA50tdOo8DQW6iHBqbUlSIQ2WKmxc08eqvt5WlyItubq7XCPi48C/iIhDVO8+vQu4KzP/uFnFSZK0VAYPjDh+Th2rkTF0rwS2ZOZYRJxEdczaWc0pS5KkpTVUrrDZQKcO1UiguxV4FjCUmY9RnVrkhvkvkSSpPQweqPDCrQOtLkNqikbG0F0N/F1E/FpEvDIilmIyYUmSmm58YpLHD9rlqs7VSKD7FPCXVFv1fgm4JSL+qSlVSZK0hB4/OMpkYperOlYjXa67MvO3pu+ICCfzkSS1vcGSkwqrszXSQndnRHxg+o7MHFnieiRJWnJHVolw2S91qEZa6LYAr4uI3wDuAL4H3JmZ/7splUmStEScVFidru5Al5k/A0e6WV8EvBg4FzDQSZLa2t5ShRU9wfFr+1tditQUDa/lWutmvaP2kCSp7Q2WKmweWElPT7S6FKkpXMtVktTx9pYqbHH8nDqYgU6S1PEGD1S8w1UdrZG1XFcCPwmcOv26zPzPS1+WJElLZ6g0witP39TqMqSmaWQM3XXAAeB2wOlKJEmFcGhknPLIuHe4qqM1EuhOzsxLmlaJJElNcGRS4Q3Oha/O1cgYulsi4sVNq0SSpCY4MgfdgC106lyNtNBdCPxcRDxEtcs1gMzMs5pSmSRJS+BIoPMuV3WwRgLdpU2rQpKkJhk8UB327V2u6mR1d7lm5iPARuBNtcfG2j5JktrW3lKFgZUrWLuy4bn0pcKoO9BFxAeATwOba49PRcS/blZhkiQthb2lCpvXe0OEOlsjf678PHBeZh4CiIgPAf8A/HEzCpMkaSkMlipsdfycOlwjd7kGMDFte6K2T5KktpSZPPrEMFvXr251KVJTNdJC9+fAtyPi2tr2m4E/W/KKJElaInftOsDjB0d4xXOPa3UpUlPVHegy8w8i4u+AC6i2zL0rM7/btMokSTpGN+4cZEVPcPEZW1pditRUDd3yk5m3U136S5KktpaZ3LRzD+c/73g2rulvdTlSUy04hi4ibq59LUdEadqjHBGl5pcoSVLj7h8s8/D+w1xy5tZWlyI13YItdJl5Ye3rQPPLkSRpady4c5AIeP0ZBjp1vkbmoftQPfskSWoHN+3cw7mnHsemAeegU+drZNqSi2fZ53JgkqS288DQQf5x70EutbtVXWLBLteI+EXgl4DnRsRd0w4NALc0qzBJkhbrpp17ALjkzBNbXIm0POq5y/UzwI3A7wJXTttfzswnmlKVJEnH4Madg7x020ZXiFDXWLDLNTMPZObDmfl2oARsAZ4DnBkRr2p2gZIkNeKH+w9zz+6S3a3qKnXPQxcR7wE+AJwM3Am8guparq9pSmWSJC3CTfdUu1svtbtVXaSRmyI+AJwDPJKZFwEvBfY1pSpJkhbpxp2DnHnSek45bk2rS5GWTSOBrpKZFYCIWJmZ9wMvaE5ZkiQ1bs+BYb77w6dsnVPXaWTpr10RsRH4IvCViHgS2N2MoiRJWoybdg4CuDqEuk7dgS4z31J7+tsR8XVgPfClplQlSdIi3LhzkOdvWcfzNq1rdSnSsmpkpYifjoip5b9eBbwLeFFTqpIkqUH7yiPc9vATzj2nrtTIGLr/mJnliLiQ6qoRnwQ+2pyyJEmqT2Zy/fd28+ar/p4A3niWgU7dp5FAN1H7+hPARzPzOqB/6UuSJKk+tz/yBG/5yC388me/y/rVfXzq58/j+VsGFr5Q6jCN3BTxWET8T+B1wIciYiWNBUJJkhYtMymPjPPkoVGGyiP8+d8/xA13D7J5YCUf/qmz+MmXnUxvT7S6TKklGgl0PwNcAvz3zHwqIk4Efr05ZUlS56uMTXD/YJl7d5f4x71lDo2MMzoxyeh47TExycRkkglJAtSeF8wsBeccn2IyYWIymczqY2ISxicmeWp4jCcPjTI++fR1q/t6+dXXPZ9feNVprOlv5NeZ1Hkaucv1MPCFadt7gD3NKEqSOtWOh5/gU7c+wj27S/zTvoNM5ZO1/b2sX91H/4oe+nt76OvtoX9FDyt6gggIAmqNTwFEkRqipuqfsXO2zxABPRH09gS9EUQEfb3BxjV9PGtNP8et7T/y9cyTNrBpYOWyfASp3S0Y6CLi5sy8MCLKPPPvrAAyM9c3rTpJ6iDfeegJ/tU132ZN/wpeespGLj1zK2c8ez0vevYGTn7WaqJQKU1SO1kw0GXmhbWvjjKVpEW6a9dTvPsTt3HSxtV8/v85nxPW2bIkaek0Mg/dh+rZJ0l6pu8PlvlX13yHjWv6+NR7zjPMSVpyjdylevEs+y5dqkIkqRM99Pgh3vHxb7NyRQ+fec8rOHHD6laXJKkD1TOG7heBXwKeGxF3TTs0ANzSrMIkqegee2qYd3zsViYz+dx7XsG249e0uiRJHaqeu1w/A9wI/C5w5bT95cx8oilVSVLBVcYmeOfHv015ZJzP/sIr+JHNDkOW1DwLdrlm5oHMfBh4B/BK4PLMfARYFxHnNrk+SSqkv75jFw8+fog/fvtLOfOkDa0uR1KHa2QM3VXA+cDba9vl2j5J0jSTk8mffeshXnzSBn78+ZtaXY6kLtBIoDsvM98HVAAy80mWeC3XiLgkIr4fEQ9ExJWzHI+I+KPa8bsi4mX1XitJy+Vv7x/iwccP8Z5XnubccpKWRSOBbiwieqlNLhwRm4DJpSqk9tpXUb1z9gzg7RFxxozTLgVOrz2uAP60gWslaVl87FsP8uwNq3jDi09sdSmSukQjge6PgGuBLRHx34Cbgd9ZwlrOBR7IzAczcxT4HHDZjHMuA/4iq24FNtbWlK3nWklqurt2PcW3H3qCd194Gn29jfxfrCQtXiNruX46Im4HXlvb9ebMvG8JazkJeHTa9i7gvDrOOanOayWp6T72rYcYWLmCt55zSqtLkdRFGlkpYiXwMmADcDzw0xHxm0tYy2wDTbLOc+q5tvoCEVdExI6I2LFv374GS5SkuT321DA33L2Ht517CgOr+lpdjqQu0kh/wHVUuzHHgUPTHktlFzD9T9qTgd11nlPPtQBk5tWZuT0zt2/a5N1nkpbOn9/8EAA/d8FpLa5EUrepu8sVODkzL2laJXAbcHpEnAY8BrwN+NkZ51wPvD8iPke1S/VAZu6JiH11XCtJTVOqjPG52x7ljWedyEkbXd5L0vJqJNDdEhEvzsy7m1FIZo5HxPuBLwG9wDWZeU9EvLd2/KPADcAbgAeAw8C75ru2GXVK0mw+950fcnBknF945XNbXYqkLhSZsw41O/rEiHuBHwEeAkaojlvLzDyreeU11/bt23PHjh2tLkNSwY1NTPKqD3+d5xy/hs9dcX6ry5HUoSLi9szcPtuxRlroLl2ieiSpo1x7x2PsOVDhv775zFaXIqlLNTJtySPNLESSiujmHzzOf/jiTl62bSMXvWBzq8uR1KXqDnQR8W9m2X0AuD0z71yyiiSpIHY8/AS/8Bc7eO6mtVzzc+fQ0+MyX5Jao5FpS7YD7+XpiXyvAF4NfCwi/t3SlyZJ7WvnYwd415/fxtYNq/hfP38eG9cs6dLWktSQRsbQHQ+8LDMPAkTEbwF/BbwKuB348NKXJ0nt5wd7y7zzz77N+tV9fOo957FpYGWrS5LU5RppodsGjE7bHgOek5nDVO96laSO98j+Q7zj499mRW8Pn37Pec45J6ktNNJC9xng1oi4rrb9JuCzEbEWuHfJK5OkFpuYTB7ef4h7dpe4d3eJe3Yf4M5Hn6K3J/j8Fedz6glrW12iJAGN3eX6XyLiBuBCqnPQvTczpyZxe0czipO0dHY/Ncy9u0s8cXiUJw+NHvl6YHiMiclkYjKZTJjM6vPZpqjM2ZdIbktzTbGZR/6n+nkyq5vjk8no+CSj4xOMTkwyOj7JgeExKmOTAPT1BqdvHuCSF23l3Reexgu2DizHx5CkujTSQgfwINWVGFYBayLiVZn5zaUvS9JSeerwKH/8tw/wF//wMGMTT6ec/t4ejlvbz/rVK1jR00NvT9DTE/QE9EbQE3PcsVmgGzlnKzWA6IGoPiMCIqC3p4f+3h5Wruihf0UPfb3BwKo+Xrh1gDOevZ7TNw/Qv6KRUSqStHwambbkPcAHqC58fyfwCuAfgNc0pTJJx2R0fJL/desj/NHXfkC5MsbPbD+Ft55zCiesW8lxa/tZ099LzBXaJEmF0kgL3QeAc4BbM/OiiHgh8J+aU5akY3HTzj387o3388j+w7zy9BP44E/8KC/cur7VZUmSmqSRQFfJzEpEEBErM/P+iHhB0yqTtCi3PPA47/3UHTx/yzo+8a5zeLWrF0hSx2sk0O2KiI3AF4GvRMSTwO5mFCVp8a67czfrVq7g+vdfyKq+3laXI0laBo3c5fqW2tPfjoivAxuAm5pSlaRFGZ+Y5Mv3DvKaF242zElSF2n0LlcAMvPvlroQScfuOw89wZOHx7j0zK2tLkWStIwauct1O/BB4DnTr8vMs5pQl6RFuHHnIKv7eh03J0ldppEWuk8Dvw7cDUw2pxxJizU5mXzpnkFe/YJNrO63u1WSukkjgW5fZl7ftEokHZM7fvgkQ+URLrG7VZK6TiOB7rci4uPA14CRqZ2Z+YUlr0pSw27cOUh/bw+veaHdrZLUbRoJdO8CXgj08XSXawIGOqnFMpObdg7yytNPYGBVX6vLkSQts0YC3Usy88VNq0TSot216wCPPTXMr7zu9FaXIklqgUZWmr41Is5oWiWSFu3GnYOs6AkuPmNLq0uRJLVAIy10FwKXR8RDVMfQBZBOWyK1VrW7dQ/nP+94Nq7pb3U5kqQWaCTQXdK0KiQt2v2DZR7ef5grXvW8VpciSWqRRpb+eqSZhUhanBt3DtIT8PoX2d0qSd2qkTF0ktrQTTv3cM6px3HCupWtLkWS1CIGOqnA/mnfQf5x70HXbpWkLld3oIuqfxkRv1nb3hYR5zavNEkLuWnnIACXnHliiyuRJLVSIy10HwHOB95e2y4DVy15RZLqduuD+3nRs9ezdcOqVpciSWqhRgLdeZn5PqACkJlPAs6RILXQngMVth23ptVlSJJarJFANxYRvVSX+yIiNvH0EmCSWmDvgQpb1ts6J0ndrpFA90fAtcCWiPhvwM3A7zSlKkkLOjQyTnlk3EAnSWpoHrpPR8TtwGtru96cmfc1pyxJC9lbqgCwdYPTlUhSt6s70EXEv5mx69KI+DHg9sy8c0mrkrSgwVqgs4VOktRIl+t24L3ASbXHFcCrgY9FxL9b+tIkzWevgU6SVNPIWq7HAy/LzIMAEfFbwF8BrwJuBz689OVJmsve0ggAWw10ktT1Gmmh2waMTtseA56TmcPAyJJWJWlBgwcqDKxcwdqVjfxdJknqRI38JvgMcGtEXFfbfhPw2YhYC9y75JVJmtfeUoUtTigsSaKxu1z/S0TcAFwIBPDezNxRO/yOZhQnaW57SxW2rPcOV0lSY12uAA8C/wDcAayJiFctfUmS6rG3NOINEZIkoLFpS94DfAA4GbgTeAXVcPeaplQmaU6Tk8neUsUbIiRJQGMtdB8AzgEeycyLgJcC+5pSlaR57T80yvhk2kInSQIaC3SVzKwARMTKzLwfeEFzypI0H+egkyRN18hdrrsiYiPwReArEfEksLsZRUma39PLfhnoJEmN3eX6ltrT346IrwMbgJuaUpWkeU0t++UYOkkS1BnoIiKAkzPzUYDM/LumViVpXnsPVOgJOGFdf6tLkSS1gbrG0GVmUu1qldQG9pZGOGHdSlb0NjrzkCSpEzXy2+DWiDinaZVIqttgqeL4OUnSEY3cFHER8N6IeBg4RHW1iMzMs5pRmKS57S1VOPlZa1pdhiSpTTQS6C5tWhWSGrK3VGH7qc9qdRmSpDbRSJfrD4FXApdn5iNAAluaUpWkOVXGJnjy8Jh3uEqSjmgk0H0EOB94e227DFy15BVJmtdQaQSAzQY6SVJNI12u52XmyyLiuwCZ+WREOGeCtMycg06SNFMjLXRjEdFLtauViNgETC5FERFxXER8JSJ+UPt61OCgiDglIr4eEfdFxD0R8YFpx347Ih6LiDtrjzcsRV1SO3KVCEnSTI0Euj8CrgU2R8R/A24GfmeJ6rgS+Fpmng58rbY90zjwbzPzR4FXAO+LiDOmHf/DzDy79rhhieqS2o7ruEqSZmpk6a9PR8TtwGupTlny5sy8b4nquAx4de35J4FvAL8x4/33AHtqz8sRcR9wEnDvEtUgFcLggQqr+npYv6qREROSpE5WdwtdRPwqcDAzr8rMP1nCMAewpRbYpoLb5gVqORV4KfDtabvfHxF3RcQ1s3XZSp1isFRh6/pVVFfkkySpsS7X9cCXIuJbEfG+iGhoypKI+GpE7JzlcVmDr7MO+GvgVzKzVNv9p8DzgLOptuL9j3muvyIidkTEjn379jXy1lJbGCqN2N0qSXqGRrpc/xPwnyLiLOCtwN9FxK7MfF2d1895XkTsjYgTM3NPRJwIDM1xXh/VMPfpzPzCtNfeO+2cjwF/M08dVwNXA2zfvj3rqV1qJ4OlCmefsrHVZUiS2shiVvYeAgaB/SzQNdqA64HLa88vB66beUJU+5f+DLgvM/9gxrETp22+Bdi5RHVJbSUzXcdVknSURsbQ/WJEfIPqXagnAL+whOu4/h5wcUT8ALi4tk1EPDsipu5YvQB4J/CaWaYn+XBE3B0Rd1Fdc/ZXl6guqa0cGB5jdHzSLldJ0jM0cpvcc6iOW7sTICIuiIirMvN9x1pEZu6nevfszP27gTfUnt9M9e7a2a5/57HWIBXB4JEpS1a2uBJJUjupu4UuM68EiIgPRcTDwH8F7m9SXZJmMXjAVSIkSUdbsIUuIp4PvI3qGq77gc8DkZkXNbk2STM4qbAkaTb1dLneD3wLeFNmPgBH5qSTtMz2lkYA2GyXqyRpmnq6XH+S6l2tX4+Ij0XE1EoRkpbZYKnCcWv7Wbmit9WlSJLayIKBLjOvzcy3Ai+kuiTXrwJbIuJPI+L1Ta5P0jR7D1TsbpUkHaWRmyIOZeanM/ONwMnAncCVzSpM0tH2litstbtVkjTDYiYWJjOfyMz/mZmvWeqCJM1t8IDLfkmSjraoQCdp+Y1NTLL/kIFOknQ0A51UEEPlETJx2S9J0lEMdFJBTM1B56TCkqSZDHRSQeytrRLhHHSSpJkMdFJBDNpCJ0mag4FOKojBUoW+3uC4tf2tLkWS1GYMdFJBDJVG2DywiggXapEkPZOBTiqIwQMV73CVJM3KQCcVxN5SxfFzkqRZGeikgthbqniHqyRpVgY6qQDKlTEOjU7YQidJmpWBTiqAvaURAJf9kiTNykAnFcCB4TEANqzpa3ElkqR2ZKCTCqBcqQa69asMdJKkoxnopAIoVcYBWL9qRYsrkSS1IwOdVABHWuhW20InSTqagU4qgNJwtYVuwBY6SdIsDHRSAZQrY6zoCVb39ba6FElSGzLQSQVQrowzsGqF67hKkmZloJMKoFQZY8A7XCVJczDQSQVQroyzfrXj5yRJszPQSQVQGh5jYKUtdJKk2RnopAKwhU6SNB8DnVQAjqGTJM3HQCcVQLky7rJfkqQ5GeikNjcxmRwcGXdSYUnSnAx0Ups7OLWOq8t+SZLmYKCT2lypto6rLXSSpLkY6KQ2NxXoHEMnSZqLgU5qc6XhWperLXSSpDkY6KQ2V55qoXMMnSRpDgY6qc2VajdFOIZOkjQXA53U5sqOoZMkLcBAJ7W5cq2Fbp0tdJKkORjopDZXGh5jdV8vfb3+5ypJmp2/IaQ2V66Ms361rXOSpLkZ6KQ2V6qMMeD4OUnSPAx0UpsrV8adg06SNC8DndTmbKGTJC3EQCe1ueoYOgOdJGluBjqpzZWGx5xUWJI0LwOd1MYyszaGzhY6SdLcDHRSGxsZn2R0YtIWOknSvAx0UhsrTS375Rg6SdI8DHRSGysNV5f9ctoSSdJ8DHRSGytPtdA5hk6SNI+2CHQRcVxEfCUiflD7+qw5zns4Iu6OiDsjYkej10tFU65UW+gcQydJmk9bBDrgSuBrmXk68LXa9lwuysyzM3P7Iq+XCsMxdJKkerRLoLsM+GTt+SeBNy/z9VJbsoVOklSPdgl0WzJzD0Dt6+Y5zkvgyxFxe0RcsYjrpUIpDVdb6Fz6S5I0n2X7sz8ivgpsneXQBxt4mQsyc3dEbAa+EhH3Z+Y3G6zjCuAKgG3btjVyqbTsypVxegLW9ve2uhRJUhtbtkCXma+b61hE7I2IEzNzT0ScCAzN8Rq7a1+HIuJa4Fzgm0Bd19euvRq4GmD79u25+E8kNV+pMsbAqj4iotWlSJLaWLt0uV4PXF57fjlw3cwTImJtRAxMPQdeD+ys93qpiMqVcdavdvycJGl+7RLofg+4OCJ+AFxc2yYinh0RN9TO2QLcHBHfA74D/N/MvGm+66WiKw2PMbDS8XOSpPm1xZ/+mbkfeO0s+3cDb6g9fxB4SSPXS0VnC50kqR7t0kInaRZTY+gkSZqPgU5qY+XKuMt+SZIWZKCT2lhpeMxJhSVJCzLQSW1qcjI5ODrusl+SpAUZ6KQ2dXB0nExYbwudJGkBBjqpTU0t++UYOknSQgx0UpsqV8YBHEMnSVqQgU5qU1MtdE5bIklaiIFOalNTLXROLCxJWoiBTmpTpYotdJKk+hjopDZ1pIXOMXSSpAUY6KQ25Rg6SVK9DHRSmyqPjLOqr4f+Ff5nKkman78ppDZVXfbL1jlJ0sIMdFKbKlfGHT8nSaqLgU5qU6WKLXSSpPoY6KQ2VaqMs361gU6StDADndSmysNjLvslSaqLgU5qU6XKOOvtcpUk1cFAJ7WpcmXMmyIkSXUx0EltaGR8gpHxScfQSZLqYqCT2tDUsl+OoZMk1cNAJ7Whp5f9MtBJkhZmoJPa0FQLnTdFSJLqYaCT2lCpMtVCZ6CTJC3MQCe1oSMtdKvtcpUkLcxAJ7Whp8fQ2UInSVqYgU5qQ0+PobOFTpK0MAOd1IZKlTEiYG2/gU6StDADndSGypVxBlauoKcnWl2KJKkADHRSGyoNjzl+TpJUNwOd1IZKlXGX/ZIk1c1AJ7WhcmXMVSIkSXUz0EltqFQZd5UISVLdDHRSGypXxpyyRJJUNwOd1IZKw2OOoZMk1c1AJ7WZycnk4Mi4Y+gkSXUz0Elt5tDoOJOJgU6SVDcDndRmnl72yy5XSVJ9DHRSmylVxgCcWFiSVDcDndRmjrTQrbbLVZJUHwOd1GZKw7bQSZIaY6CT2szTY+hsoZMk1cdAJ7UZx9BJkhploJPazFQLndOWSJLqZaCT2kxpeIz+FT2s6uttdSmSpIIw0EltplQZd/ycJKkhBjqpzTxxaITj1va3ugxJUoEY6KQ2M1QeYdPAylaXIUkqEAOd1GaGSiNsHljV6jIkSQVioJPaSGay7+AIm22hkyQ1wEAntZHS8Dij45N2uUqSGtIWgS4ijouIr0TED2pfnzXLOS+IiDunPUoR8Su1Y78dEY9NO/aGZf8Q0hIYKlcADHSSpIa0RaADrgS+lpmnA1+rbT9DZn4/M8/OzLOBlwOHgWunnfKHU8cz84blKFpaavvKIwCOoZMkNaRdAt1lwCdrzz8JvHmB818L/FNmPtLMoqTlNlQLdLbQSZIa0S6Bbktm7gGofd28wPlvAz47Y9/7I+KuiLhmti5bqQimulw3rzfQSZLqt2yBLiK+GhE7Z3lc1uDr9AP/HPjf03b/KfA84GxgD/A/5rn+iojYERE79u3b1/gHkZpoX3mEVX09DKx0pQhJUv2W7bdGZr5urmMRsTciTszMPRFxIjA0z0tdCtyRmXunvfaR5xHxMeBv5qnjauBqgO3bt2cDH0FquqlJhSOi1aVIkgqkXbpcrwcurz2/HLhunnPfzozu1loInPIWYOeSVictEycVliQtRrsEut8DLo6IHwAX17aJiGdHxJE7ViNiTe34F2Zc/+GIuDsi7gIuAn51ecqWlpaTCkuSFqMtBupk5n6qd67O3L8beMO07cPA8bOc986mFigtk6FShR973lH/xCVJmle7tNBJXa8yNkGpMm4LnSSpYQY6qU04qbAkabEMdFKbcFJhSdJiGeikNrHPdVwlSYtkoJPaxJEuV1eJkCQ1yEAntYmh8gg9AcevNdBJkhpjoJPaxFBphOPXraS3x1UiJEmNMdBJbcJJhSVJi2Wgk9rEULniDRGSpEUx0EltorqOq4FOktQ4A53UBiYmk/2HRp1UWJK0KAY6qQ08cWiUicm0y1WStCgGOqkNPL3sl4FOktQ4A53UBoZqq0Q4qbAkaTEMdFIbOLKO6zrH0EmSGmegk9qAy35Jko6FgU5qA/vKIwysWsGqvt5WlyJJKiADndQGnFRYknQsDHRSG9hXdlJhSdLiGeikNjBUHnFSYUnSohnopBbLTIZKI3a5SpIWzUAntdih0QmGxybscpUkLZqBTmqxoZKTCkuSjo2BTmoxJxWWJB0rA53UYk4qLEk6VgY6qcWmWugcQydJWiwDndRiQ+UK/b09bFjd1+pSJEkFZaCTWmxfuTplSUS0uhRJUkEZ6KQWmwp0kiQtloFOajEnFZYkHSsDndRi+w66jqsk6dgY6KQWGh2f5IlDo67jKkk6JgY6qYUeP1ibVNgWOknSMTDQSS20zznoJElLwEAntdCQq0RIkpaAgU5qoaFyBbDLVZJ0bAx0UgvtK48QASesM9BJkhbPQCe10FB5hOPW9NPX63+KkqTF87eI1EJOKixJWgoGOqmF9h000EmSjp2BTmqhfaWKkwpLko6ZgU5qkcysLvvllCWSpGNkoJNa5KHHDzE2kWxdbwudJOnYGOikFvlftz5CX29wyZlbW12KJKngDHRSC5QrY/zvHbt441nPZostdJKkY2Sgk1rgL3fs4uDIOO++4LRWlyJJ6gAGOmmZTUwmn7jlIc459Vm8+OQNrS5HktQBDHTSMvvqfXt59IlhW+ckSUvGQCcts2tufoiTNq7m4jO2tLoUSVKHMNBJy2jnYwf49kNP8HM/diorXL9VkrRE/I0iLaNr/v4h1vT38jPnnNLqUiRJHcRAJy2ToXKF//O93fz0y09mw+q+VpcjSeogBjppmXzq1h8yNpH8nDdDSJKWWFsEuoj46Yi4JyImI2L7POddEhHfj4gHIuLKafuPi4ivRMQPal+ftTyVS/WpjE3w6Vsf4bUv3MxpJ6xtdTmSpA7TFoEO2An8C+Cbc50QEb3AVcClwBnA2yPijNrhK4GvZebpwNdq21LLZSZ3/PBJ/t9r72b/oVHefaGtc5Kkpbei1QUAZOZ9ABEx32nnAg9k5oO1cz8HXAbcW/v66tp5nwS+AfxGc6qV5jc5mXz30Se54e5Bbrx7D7sPVOjv7eFntp/Mjz3v+FaXJ0nqQG0R6Op0EvDotO1dwHm151sycw9AZu6JiM3LXdxsvn7/EH91x65Wl6FG5Wy7kkzIhMlMkmrr2/hkcnh0guHRCQ6PjnN4dIJyZZyDI+P09/bwquefwK/9sxfwujO2sH6VN0JIkppj2QJdRHwV2DrLoQ9m5nX1vMQs+2b51btgHVcAVwBs27at0csbsv/QKPfvKTX1PdQcs7UWB9ATQUT1eAC9PcGa/l5OWNfPmpVrWNPXy9qVK3jJKRt47Y8a4iRJy2PZAl1mvu4YX2IXMH3yrpOB3bXneyPixFrr3InA0Dx1XA1cDbB9+/aGA2EjfurlJ/NTLz+5mW8hSZLUNjdF1OM24PSIOC0i+oG3AdfXjl0PXF57fjlQT4ufJElSR2iLQBcRb4mIXcD5wP+NiC/V9j87Im4AyMxx4P3Al4D7gL/MzHtqL/F7wMUR8QPg4tq2JElSV4jMpvY6trXt27fnjh07Wl2GJEnSgiLi9sycdb7etmihkyRJ0uIZ6CRJkgrOQCdJklRwBjpJkqSCM9BJkiQVnIFOkiSp4Ax0kiRJBWegkyRJKjgDnSRJUsEZ6CRJkgrOQCdJklRwBjpJkqSCM9BJkiQVnIFOkiSp4Ax0kiRJBWegkyRJKjgDnSRJUsEZ6CRJkgrOQCdJklRwBjpJkqSCi8xsdQ0tExH7gEdaXUeLnAA83uoiVDd/XsXiz6tY/HkVSzf/vJ6TmZtmO9DVga6bRcSOzNze6jpUH39exeLPq1j8eRWLP6/Z2eUqSZJUcAY6SZKkgjPQda+rW12AGuLPq1j8eRWLP69i8ec1C8fQSZIkFZwtdJIkSQVnoBMR8WsRkRFxQqtr0dwi4vcj4v6IuCsiro2Ija2uSc8UEZdExPcj4oGIuLLV9Wh+EXFKRHw9Iu6LiHsi4gOtrknzi4jeiPhuRPxNq2tpNwa6LhcRpwAXAz9sdS1a0FeAMzPzLOAfgX/f4no0TUT0AlcBlwJnAG+PiDNaW5UWMA7828z8UeAVwPv8mbW9DwD3tbqIdmSg0x8C/w5wMGWby8wvZ+Z4bfNW4ORW1qOjnAs8kJkPZuYo8DngshbXpHlk5p7MvKP2vEw1KJzU2qo0l4g4GfgJ4OOtrqUdGei6WET8c+CxzPxeq2tRw94N3NjqIvQMJwGPTtveheGgMCLiVOClwLdbXIrm9v9RbYCYbHEdbWlFqwtQc0XEV4Gtsxz6IPD/Aq9f3oo0n/l+Xpl5Xe2cD1LtKvr0ctamBcUs+2z5LoCIWAf8NfArmVlqdT06WkS8ERjKzNsj4tUtLqctGeg6XGa+brb9EfFi4DTgexEB1e67OyLi3MwcXMYSNc1cP68pEXE58EbgtemcQ+1mF3DKtO2Tgd0tqkV1iog+qmHu05n5hVbXozldAPzziHgDsApYHxGfysx/2eK62obz0AmAiHgY2J6Z3brgcduLiEuAPwB+PDP3tboePVNErKB6s8prgceA24Cfzcx7WlqY5hTVv2Y/CTyRmb/S4nJUp1oL3a9l5htbXEpbcQydVBx/AgwAX4mIOyPio60uSE+r3bDyfuBLVAfX/6Vhru1dALwTeE3tv6k7ay1AUuHYQidJklRwttBJkiQVnIFOkiSp4Ax0kiRJBWegkyRJKjgDnSRJUsEZ6CRJkgrOQCdJklRwBjpJixIRE9MmY70zIq6c59yNEfFLy1TXwTn2/3JE3BcRn65t33IMr/XbEfFrdVy/nJ/7Ge8VEadGxM4leu1TI2I4Iu5c4LxvRMQ/m7HvVyLiIxGxuvbvZDQiTliKuiQ9zUAnabGGM/PsaY/fm+fcjUBDwSaqlvL/o34JeENmvgMgM39sCV97LhuZ43M34fPN+V5L5J8y8+wFzvks8LYZ+94GfDYzh2vXu76t1AQGOklLJiLOiYi7ImJVRKyNiHsi4kzg94Dn1Vpofr927r+MiO/U9v3PiOittQTdFxEfAe4AXlnb/ljttb4cEatr138xIm6v7b9igbo+CjwXuD4ifrW27+C040fVMstrfDAivh8RXwVeMMvxtRHxfyPiexGxMyLeOvNzz/L5Tlng+3DU566913+MiPsj4isR8dlaa+FR32Ogd67XmFH7NyLiBbXnx9fTsjfH9+yvgDdGxMraOacCzwZuXuj1JB0bA52kxZrqQpt6vDUzbwOuB/4r8GHgU5m5E7iSWgtPZv56RPwo8FbgglqrzQTwjtrrvgD4i8x8KfAIcDpwVWa+CHgK+Mnaee/OzJcD24Ffjojj5yo0M99LtWXoosz8w+nHFqhl6pyXU21peinwL4BzZnmbS4DdmfmSzDwTuGnm557l862Z571n/dwRsb32fKqW7bXzZ3uvub53M/0I8IPa87OAu+c4b+r7Mev3LDP3A9+pfS+g+j37fLrGpNR0K1pdgKTCGp6jC+4/A7cBFeCX57j2tcDLgdsiAmA1MAR8E3gkM2+ddu5DmXln7fntwKm1578cEW+pPT+FanjZv4jPMVct070SuDYzDwNExPWzvM7dwH+PiA8Bf5OZ34qIZ81y3vTPN9/3Ya7PfSFwXWYO12r5P/N8trle44iIeA7wWGZO1nadBdw1z2vOVzc83e16Xe3ruxd4LUlLwEAnaakdB6wD+oBVwKFZzgngk5n575+xs9pFN/P8kWnPJ6i2DL4aeB1wfmYejohv1N5rMWatZRbztjJl5j/WWvLeAPxuRHwZ+ItZTp3++eb7Phz1uaddU6+5XmO6s3lmgHs58PkFXne+79kXgT+IiJcBqzPzjrqrlbRodrlKWmpXA/8R+DTwodq+MjAw7ZyvAT8VEZsBIuK4WktRvTYAT9bC3AuBVxxDvfXU8k3gLVG9U3MAeNPMF4mIZwOHM/NTwH8HXsbRn3sx7z3TzcCbauMU1wE/Udu/0HvN5SXUwnBEnA5cxgJdrvPVnZkHgW8A11BtrZO0DGyhk7RYq+OZ01jcBNwLjGfmZ2qD5G+JiNdk5t9GxN/XBtvfWBtH9x+AL0f1Ts8x4H3AYJ3vfRPw3oi4C/g+cOsC588pM++do5ZHpp1zR0R8Hriztv9bs7zUi4Hfj4jJ2mv8Ymbun/65gavqfO85vw+ZeVuty/d7tVp2AAcWeq95nA0MR8T3qLbU3QdcDvyXeWpY6Hv2WeALHH3Hq6QmCceqSlKxRMS6zDwYEWuoth5esdiuzYh4AHhpZpbnOedUquMCz1xUwc98rYeB7Zn5+LG+lqSn2eUqScVzda119A7gr48hzA0Ak/OFuZoJYEMsMLHwAu811aLbB0wucLqkBtlCJ0mSVHC20EmSJBWcgU6SJKngDHSSJEkFZ6CTJEkqOAOdJElSwRnoJEmSCs5AJ0mSVHAGOkmSpIL7/wF2+RZGaTghDQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Parameters\n",
    "M = 19      # Size of the spin chain \n",
    "J = 1      # Exchange interaction\n",
    "kT = 0.1   # Temperature\n",
    "\n",
    "# Prepare plots\n",
    "df = pd.DataFrame()\n",
    "\n",
    "# Brute-Force (green)\n",
    "for mu in np.linspace(-5,5,100):\n",
    "    bf = transfer_matrix(M,J,mu,[kT])\n",
    "    df = pd.concat([df,bf])\n",
    "\n",
    "plt.figure(figsize=(10,10))\n",
    "print(list(df))\n",
    "plt.plot(df['mu'],df['m'])\n",
    "plt.xlabel('External field strength $\\mu$ [eV]')\n",
    "plt.ylabel('Average magnetisation $m$ [-]')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analytical solution\n",
    "\n",
    "The transfer matrix can be solved analytically. The eigenvalues of the transfer matrix\n",
    "\n",
    "$$\n",
    "{\\bf T} = \\left[\n",
    "\\begin{array}{cc}\n",
    "\\exp(J-\\mu) & \\exp(-J+\\mu) \\\\\n",
    "\\exp(-J-\\mu) & \\exp(J+\\mu) \\\\\n",
    "\\end{array}\n",
    "\\right]\n",
    "$$\n",
    "\n",
    "are\n",
    "\n",
    "$$\n",
    "\\lambda_\\pm = \\exp(-\\beta J)\\cdot\\cosh(\\beta\\cdot\\mu) \\pm \\sqrt{\n",
    "    \\exp(-2\\beta J)\\cdot\\sinh^2(\\beta\\cdot\\mu) + \\exp(2\\beta/ J)\n",
    "}\n",
    "$$\n",
    "\n",
    "Once the eigenvalues are known, the partition function becomes \n",
    "\n",
    "$$\n",
    "Z = \\lambda^N_+ + \\lambda^N_- \\quad {N\\rightarrow\\infty} \\quad Z \\approx \\lambda^N_+\n",
    "$$\n",
    "\n",
    "In the thermodynamic limit ($N\\rightarrow\\infty$) the largest eigenvalue dominates and the thermodynamic properties can be computed from\n",
    "\n",
    "$$\n",
    "g = -kT\\cdot\\ln(\\lambda_+)/N\\quad\n",
    "u = -\\frac{\\partial \\ln(\\lambda_+)}{\\partial \\beta}/N\\quad\n",
    "m = \\frac{\\partial g}{\\partial \\mu}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def analytical_solution(J,mu,kTs):\n",
    "    '''\n",
    "    Solve the one-dimensional Ising model analytically through the \n",
    "    transfer matrix formalism.\n",
    "    \n",
    "    The function returns the results as a Pandas dataframe containing the free energy g,\n",
    "    the internal energy u, and the magnetisation m\n",
    "    \n",
    "    Parameters:\n",
    "        \n",
    "        J   - exchange interaction\n",
    "        mu  - external field (same units as J)\n",
    "        kTs - array with temperatures (in units of J)\n",
    "    \n",
    "    Returns:\n",
    "    \n",
    "        A Pandas dataframe with free energy g, internal energy u, magnetisation m\n",
    "        \n",
    "    Usage: \n",
    "        \n",
    "        analytical_solution(-1,0.1,np.linspace(0.1,4,30))\n",
    "    '''\n",
    "    \n",
    "    betas = 1/kTs\n",
    "\n",
    "    # Create a new data-frame to take results\n",
    "    df = pd.DataFrame(columns=['kT','g','u','m'],dtype=object)\n",
    "    \n",
    "    ## Eigenvalue\n",
    "    lmbda = np.exp(-betas*J)*np.cosh(betas*mu) \\\n",
    "          + np.sqrt(np.exp(-2*betas*J)*np.sinh(betas*mu)**2 + np.exp(2*betas*J))\n",
    "    \n",
    "    # Free Energy\n",
    "    g = -1/betas*np.log(lmbda)\n",
    "\n",
    "    # Internal energy (fairly long expression)\n",
    "    ua = np.exp(J*betas)*np.cosh(betas*mu)\n",
    "    ub = np.exp(J*betas)*mu*np.sinh(betas*mu)\n",
    "    uc = mu*np.cosh(betas*mu)*np.sinh(betas*mu) + J*(np.exp(4*J*betas)-np.sinh(betas*mu)**2)\n",
    "    ud = np.sqrt(np.exp(-2*betas*J)*(np.exp(4*J*betas)+np.sinh(betas*mu)**2))\n",
    "    ue = np.exp(-J*betas)*np.cosh(betas*mu)\n",
    "    u = -np.exp(-2*J*betas)*(-J*ua+ub+uc/ud)/(ue+ud)\n",
    "\n",
    "    # Magnetisation\n",
    "    m = np.exp(-betas*J)*np.sinh(betas*mu)/np.sqrt(np.exp(-2*betas*J)*np.sinh(betas*mu)**2 + np.exp(2*betas*J))\n",
    "\n",
    "    df['u'] = u\n",
    "    df['g'] = g\n",
    "    df['m'] = m\n",
    "    df['kT'] = kTs\n",
    "    # df['mu'] = np.ones(len(mus))*mus\n",
    "    return df\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Excercise 3: Scaling analysis\n",
    "\n",
    "Let's first check that all three methods yield similar results and then see how the \"Brute Force\" method scales with chain length compared to the \"Transfer Matrix\" method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Brute force solution took 0.49805712699890137 sec\n",
      "Avg. error Brute Force free energy       : 0.0006439824293777207\n",
      "Avg. error Brute Force magnetisation     : 0.01060166266750696\n",
      "Transfer matrix solution took 0.02721095085144043 sec\n",
      "Avg. error Transfer Matrix free energy   : 0.0006439824293775987\n",
      "Avg. error Transfer Matrix magnetisation : 0.010601662667506898\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4cAAAFzCAYAAACNcXywAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAACAYUlEQVR4nOzdd5gUVdrG4d87Mww5SM5BkmQQRBBFJQiigigoKIJgzmHXVdfV1dVdc86oIApiAAkqggQxB4JIFMk5KRkkDPN+f3T7OeIAA8zMmel+7uuqq7urT3c95bI1/VadOsfcHREREREREYlvCaEDiIiIiIiISHgqDkVERERERETFoYiIiIiIiKg4FBEREREREVQcioiIiIiICCoORUREREREBEgKHSA7lSxZ0qtWrRo6hohkomnTpv3i7qVC5zgaOjaJxCYdn0QkJzrYsSmuisOqVasyderU0DFEJBOZ2bLQGY6Wjk0isUnHJxHJiQ52bFK3UhEREREREVFxKCIiIiIiIioORUREREREBBWHIiIiIiIigopDERERERERQcWhiIiIiIiIoOJQREREREREUHEoIiIiIiIiqDgUERERERERVByKiIiIiIgIKg5FREREREQESAodICcaPmg8c+ZuOGS7BPvjeWICJNqfH5MSjQJ5EyiYN5ECeRMplC+JQvnzULJYfiqWL0FCkcJQsCAUKgR584LZgTcmIiIiucKu7ZsZ897/QscQETlsKg7T8c/HN/PzrIuydiMJe0nI9yt58q4jb/JcCub5hfL5V9K42GpaV9pDx/olKF23BtSsCQ0bQpL+pxIREcmp9u3dw2cfPMuQb15mmK9k69a6oSOJiBw2VRzpePnBSixaOPqgbdz/eJ7q4O7sSyW6OPsc9uxN5bc9zm+7U9m1B37b7fy2B7Zsh183J7J5azJbthVg545CbNxZlzVrOjFtXz5emwK8n0pS0cWUKjSL0/K/zh3H76Thee2hY0coWjRr/wOIiIjIIaXuS+WDoR8wYMQXfLq2ENs21cfWfYBvrEnkzh31CBKR3EXFYTpOO6sFpwXY7p69+/h0+hImfreOqdN2sGBOHtYuacLQVV0ZujCVwpO/pk2B+/h75SWc3K0N9OsX6ZYqIiIi2Wb6l7O4+4EpTJzSgt0buwBdgFTKlf6F5q2K0vj4BBo0gG7dQicVETk8Kg5zkOQ8iXQ4sRodTqz2/+tSU51RXyzg2ddX8s0nZRm19AlGLYWC877jxmev4IErm5Jw7TVQoEC44CIiIjFu47qN/Puu0QydUIlfl58G3oBC5b7j3K6vc/nl7Wh5akUKFiwdOqaIyFHRaKU5XEKC0fXUmkwaeDq/rarDxKnL6HTtZPaklubBBW9R+OFTubdJX1Kfehp27QodV0REJGak7ktl8ItjqNfkbUpUSeS51y5l88aatGo1mLHDvmDb6hN5+/1LadepojryiEhMUHGYy7RpWoWPnj+Nrasrcdl9X7IntRT3/fwORf7Xgv82voTU1wb8+YZIEREROSxLf1rOBV1fo1C5WVxybSfmzu5MteoTeezed9n1a3m+/KIPHc4/JXRMEZFMp+Iwl8qXnMSr95zMllXl6f2vz9m1tzz/mv8eJf9VhW8694GNG0NHFBERyTVS96Xy7IPDObbuCKo1KMV7Iy8DnB7nvcbSOb+weM55/O3fF5CUR3fkiEjsUnGYyxXIl4dB97dm8+oydP/7Z2za2JyTxj/DDa3+AZMnh44nIiKSo61cvIoLu75GgTLzuPGf57NkcRvq1R/F4BfHsHN9Y4YOv4wqtSqHjikiki1UHMaIQvmTeffRU5n49UYKlVvKcz+9SrXeK1n293/B3r2h44mIiOQo7w34hHqN36HScUV5d+RlJCbupk+PAaxbuJvZP/Tg4qs7hY4oIpLtVBzGmDZNq7Bhfn3a9J7I0pU9qN7/Cp495TJYtix0NBERkaB27dzFLVcNomj5b7jgsjOYO6czNWqN4/VnP2THuuN5fWg/SlfMnSOOmtkAM1tvZrMP8L6Z2TNmttDMZprZ8dmdUURyPhWHMShfchITB7Xl1VE/kZC8jxu/e50zznid1Bk/ho4mIiKS7RbNXsyZZwygUNn1PNW/Dzt3lqZzp9dYOGMtC2afT5/rzw4dMTO8DnQ8yPtnAjWjy5XAi9mQSURyGRWHMeyyc+qz/OfSVDnxc8b//G/qn/c9uyZODh1LREQkW4weOom6jd+hRpOyjB3fjyJFl/DPm99kx9pKjProMqrXq3boL8kl3P1z4GCj0XUB3vCIb4FiZlYue9KJSG6hIbdiXNnihVj89am06jaOb0dcQbU+o5n12HBK9jg/dDQREZFM9/uoo/8bUJT1S86ApN+oVfcD/nfrMZzfp33oeCFVAFakeb0yum5NZnx58cqfsXl9Q7BUMAcczDFSSc6/jiL5V1Iqz0oqJq6kOitptncpfXbMJiFPHsiTB5KSIkuxYlC/PjRsGFkaNICSJTMjoohkgIrDOJCQYHzzfgd63PwJ7zx9NtVu+o6pq1+i9q1Xh44mIiKSKfbs2sOdtwzl5RH12LGuO1ZgPaefPoAXn2xN7UYXhI6XE1g669KdGNnMriTS9ZTKlTM2Umur+gtZvn4V7oZj0UdITU1g885ibNlamQ2/tmDObyX+/zM3VpjCAye8z03F1kJKSmTZsAE++AAGDPjjy8uWhb594f77ITExwzssIodPxWEcefupM6hS+Qseub0ZDe4vwSerH+S0R+8AS+/vhYiISM63af0mrrn6fYZPOpWULX3IU2wBF3V7jWef60rxMv1Cx8tJVgKV0ryuCKxOr6G79wf6AzRr1izdAnJ/H4y5LEMh1m/awdSf1jHkg5W8/UJNbh75IPc3mMrTjxbk4g51/mi4bh3MnAmzZsEXX8CDD0aev/UWFC6coW2JyOHTPYdx5uFbT+Glt34iZXcJ2rzYjxHXPRA6koiIyGFbt2I9nTu9Rsljd/HOiMvIm38Tt133BtvXVGHIe5dRvEzx0BFzmtFA7+iopS2ALe6eKV1KD0fpYwrSqeWxDPlfazYsL8ZZ101m46Jj6dWxDpVafsP475dGGpYpA+3bw623wogR8Pzz8PHHcPLJsHx5dscWiRsqDuPQVd2b8NGEzVhiCt0GXcG4O58IHUlERCRDli9YTvu2AyhXO4EPPr6MosV/5vH73mPryqY88lxvkvMlh44YhJkNBb4BapvZSjO7zMyuNrPf7yEZAywGFgKvANcGivr/ihfJz4fPncbSxYmccslkVk5vwBmtizPq84V/bXzttfDRR7B0KTRvDt9/n+15ReKBisM4deZJ1fngox04CZz1XHe+fKh/6EgiIiIHtGjOEk455XWqNCjChEn9KFX2B15+YhQbl5/Krfd0JyExvn/SuHtPdy/n7nncvaK7v+buL7n7S9H33d2vc/fq7t7A3aeGzvy7ymWK8vkbp/HllM0kJP9Gt/PysGBFOgOvdugA33wDBQrAqafCe+9lf1iRGBffR9I41+mUWrw7YgOpKQU5/cG2TH1xSOhIIiIif7J84QpOP20gNZoW58svL6V8pa8Y8vLHrFvcnitv6RI6nmSiVg0r8vLgDaRsKUPzDsvYuWvvXxvVrQvffQdNm8IFF8CoUdkfVCSGqTiMc93OqMdrg5eSsqsUre46njlDdJAVEZHw1q1YT8czBlC1QUEmf9aXcpW+Ytjr41m14CwuuvLM0PEki1zeuT5X/Hsqm+c14cRuX6ffqFQpmDgRataEhx/O3oAiMU7FodD3/ON55uU57NlWhWY3VWHRBxNDRxIRkTi1af0mupz1GuWOg3Hj+1Gy3FQGvziG1Qs6xfs8hXGj/79Oplm3ycz+6FQuuvPz9BvlzQvXXRfpZjptWvYGFIlhKg4FgBsubcn/nvyeXZvr0Ojywqz/fkboSCIiEkf27NrDFb0GUKrGdkaPuYxjSs7l5SdGsX7xGVx8dafQ8SSbffXWKZRqPIWhj7TkqbdnpN/o0kuhYEF47rnsjCYS04IUh2bW3czmmFmqmTU7QJtKZvapmc2Ltr0pzXv3mtkqM5sRXfRXIxPcef1p3HHfZHasb0aTC5aRsuGX0JFERCTGpe5L5f7b36Jopfm8OqQf+Qqs4/H73uPXZafpnsI4lpwnkalja5Fcejm3XlaJyT+kM31F0aLQuzcMHQobNmR/SJEYFOrK4WzgPOAAfQUASAH+5u51gBbAdWZWN837T7p74+gyJguzxpUH7+pAl14jWb2sCye3fRf2pnMzuIiISCZ48/mPKFH1S+555CJS9hbg5isHsXXV8dx6T/fQ0SQHqFymKGM+SAJP4LxeBzhhfd11sHs3vPZa9oYTiVFBikN3n+fu8w/RZo27T48+3wbMAypkR7549/6grtRtOZrvZl1L305PgHvoSCIiEkOmTJ5BtToj6X39WWz5tQ7nd3mNTcvK8OTLfeJ+Sgr5s7bNqtCp349smns8gz6a+9cG9epBmzbwwguQkpL9AUViTK44AptZVaAJ8F2a1deb2UwzG2Bmxxzks1ea2VQzm7pBXQ4yJCHBmDLpTIof+wWvT7qVJ655OnQkERGJARvX/kq7NgNofkZNli7oyEmtBrF89h6GjbyMQkULhY4nOdRL/z4e8m3hrvu3pN/g+uthxQoYPTp7g4nEoCwrDs1sgpnNTmc5rBsIzKwQMBy42d23Rle/CFQHGgNrgMcP9Hl37+/uzdy9WalSpY5sZ+JQgXx5mP5FffIWX8jfB/Xm42feCR1JRERyqdR9qfzt6tcpXXsbEz/tR4Wqk/nsw1l89WUfKh6rTkFycBVLFeGk835g1Xcn8vG3i//a4JxzoHJlDUwjkgmyrDh093buXj+dJcMT6ZlZHiKF4RB3fz/Nd69z933ungq8AjTP/D2QKuWP4ZPRiZCQQud7mrJg4rehI4mISC4zbMA4SlT9iidevpQ8yVt5/P5hrPz5LFp3PCF0NMlFXr6/HiTt5pZ7Vv71zaQkuOYa+PRTmD07+8OJxJAc263UzAx4DZjn7k/s9165NC+7EhngRrJA65a1ePHZn0jZXpHmvXeza71GMBURkUNbvXQ1TZsPpvvlbdnyax16dhvApuW1ufVf3UJHk1yo/rGlqHfGFOZPasHUn9b8tcHll0fmPnz++ewPJxJDQk1l0dXMVgItgY/MbFx0fXkz+33k0VbAJUCbdKaseMTMZpnZTOB04Jbs3od4clW/1lx1zUdsXn0qrTu8D6mpoSOJiEgOlbovlduueZ3KDVKZPuUi6jQczoKp23nrvX7ky583dDzJxZ57oCqkJnLN3emMaViyJPTsCW+8AZs3Z3c0kZgRarTSEe5e0d3zunsZd+8QXb/a3TtFn3/p7ubuDfefssLdL3H3BtH3Ort7OqeQJDO99Oz51DtxBFNmXMltvZ4NHUdERHKgz8d8S9kaE3nspUvJk/9X+j/5IXNnXEj1ulVDR5MYcFqTylQ+6Tumjm7KktWb/9rghhtg5054/fXsjiYSM3Jst1LJeb6ecCYFy03lseF9GfPyyNBxREQkh9i5bScd2w3g1C6N2LC6BV3Ofo1NS2tzxc2dQ0eTGPPIvSVgT2GuunfGX988/nho2TLStVS9nESOiIpDybAihfIxcVRRLHE3Xe+sxdq5C0NHEhGRwN55ZSwlq//MuIn9KF/1U74bv4iRH1xGvgL5QkeTGHRhu9qUaDCVCUPrsXHrb39tcPXVsHAhTJuW/eFEYoCKQzksJ55Qk4f/M4U9m2vTvPMsUnfvCR1JREQC2LZxK61PGUiPq9uye2dZ/n79G6xa0InmrRuHjiYx7l//TMK3l+K6/07565tt2kQev/oqe0OJxAgVh3LYbvt7Jzqc/Q4rFnWl21n9Q8cREZFsNujp0ZSqtYIvvuzLsXU+YP53v/Hos71Dx5I4ceMFjShYbTbDXqvCrj0pf36zYsXInIdffx0mnEgup+JQjsiYkT0pW/MTRnx6FS/8a1DoOCIikg22/rKZFi0GcektZ5Gypwj/vn0oi2afR4161UJHkziSkGBcd8t2Un6twl3Pp3P1sFWryJVD9+wPJ5LLqTiUI5KQYHw3oQlJRZZz4zOns2jKnNCRREQkCw0fOI4ydVbw3Xd9qN1gOItnJHDvQz1Dx5I4df/VJ2CFNjBseDoFYKtWsHo1LFuW/cFEcjkVh3LEKlcuxWtPrWDfjrKccuESUvfsDR1JREQyWcqevXTp9CrdrjiNPb+V4t+3D+WnHy+g8rEVQkeTOJacJ5GaLX5i+dR6bP9tv/EPWrWKPOq+Q5HDpuJQjkrvPqdxdud3WbPkbHqdp/sPRURiyfcTplG6+teM/vhyylb9lNlfbdfVQskxup2XDLuL8ty7s/78RoMGULiwikORI6DiUI7aiGEXUaLqpwwd148RL40OHUdERI6WOzdfPpAWZx/LpnVN6X3RQFbNP4M6jWqETiby/265qAHk2clbw7b/+Y3ERGjRQsWhyBFQcShHLSkxgcmjqpKQdws9/3UsvyxdHTqSiIgcoV9XbaB2g2E8/VpfCh6zgAnDZzNoSF8SEvWTQXKWkkULULbRLOZ+WYPU1P3uPWzVCmbNgi1bwoQTyaV0pJdMUb9hNR64/Vt2/1qf1udM1ghhIiK50Jihk6jQaB0/z+nOSSe/wfpFDWl7TovQsUQO6Myz9rBvcwXemTD/z2+0ahX5LfLtt2GCieRSKg4l09x597k0P/kt5s2+iNsu0/2HIiK5hjvX9nmVs/o0Z8/OMjxw1zt89UVv8hfIFzqZyEH9rc9xQCovv7X2z2+ceCIkJKhrqchhUnEomWri2PMoWHYaj791Ad98pAOyiEhOt3H1Bo5r+C4vvnE5RUrP5LtP1nLXAxeGjiWSIfWqlaJwjdl8P6nMn98oXBgaNVJxKHKYVBxKpipUMB+j3zDckzjr6r2k7NodOpKIiBzAhGGfUaHhWubPvpCTTnmDdQuacsLJDULHEjksp5yxkd9W1OGb2av+/EarVvDdd5CSEiaYSC6k4lAyXZv2x9P7wpFsWnka3c8dGDqOiIik4/6/vUn7Xo3YvbMc//nnO3z1eW/y5c8bOpbIYbv+ksoAPDFo4Z/fOOkk2LEDfvwxQCqR3EnFoWSJga/3okz1Txg54VJGDxwbOo6IiESlpuzjzHavcs8TF5OvyHI+G72Mu/+rbqSSe53Z4liSSy9m0seF/vxGq1aRx6+/zv5QIrmUikPJEgkJxoRhx2LJ2+hxe2m2b9RQ0iIioW1Yvo5qdT9k7MTLqVT7I5bPrMAp7ZqGjiVy1BqftpyN8xqyfF2a3xuVK0PFirrvUOQwqDiULFO/cQ1uu2YCv204njPPHh46johIXPvyo2+o3HQtyxecQ8czX2fpnLMoVbZE6FgimaJfjxKQmofH3pjz5zdatVJxKHIYVBxKlnr48Z5Urz+CL7/tw6sPvxc6jsgRM7PiZjbezBZEH485SNtEM/vBzD7MzowiB/LMfW/TuntNdm2rwn/uGsbHYy7VpPYSU/qeXRcruIFRo/Z7o1UrWLkSli8Pkkskt9FfBslykz5sSWLhFVz7UGPWL10dOo7IkboDmOjuNYGJ0dcHchMwL1tSiRzCtZe8yk3/OY88+X9h3Hs/cfcDF4SOJJLpkvMkUqPFTyyfVo/tv+35443f7zvU1UORDFFxKFmucpWyPH7Xj+zdXJ02XSaHjiNypLoAg6LPBwHnptfIzCoCZwGvZk8skfSlpuyj3amv8uLgyylW/jt++rYgZ5zTInQskSzT/fxk2FWUF4bN/mNlw4ZQsKCKQ5EMUnEo2eKmf3Shacu3mTPzIh65483QcUSORBl3XwMQfSx9gHZPAf8AUrMpl8hf7Ny0leMavc/Ezy+nWr2RLJ9zPNVqVgodSyRL3XJRA8izk8HvbvtjZVIStGih4lAkg1QcSrb5ZPRZJB8zn3++0JrVi1aEjiPyF2Y2wcxmp7N0yeDnzwbWu/u0DLS90symmtnUDRs2HHV2kd+tmLeESo2msWBud045/U0W/tiZwkUKho4lkuVKFi1A2YazmPtlDVJT/Y83WrWCmTNh27YDf1hEABWHko2KlyzKs/ctZN+2SrQ/78vQcUT+wt3buXv9dJZRwDozKwcQfVyfzle0Ajqb2VLgbaCNmQ0+wLb6u3szd29WqlSpLNojiTffjvmWmqfsYOOqk+l36Zt8PukSDTwjcaVDpz3s21yBUV8s/GNlq1aQmgrffhsumEguob8Ykq2uvOEsmrZ4m7kze/LoXUNCxxE5HKOBPtHnfYD9x8TD3e9094ruXhXoAUxy917ZF1Hi2cdvfsLJF1Zg9/ZKPPyfD3ht4CWhI4lku95dKgPw3tg1f6xs0QISEtS1VCQDVBxKtvtkdCfyHPMzdz7XitWLV4aOI5JRDwHtzWwB0D76GjMrb2ZjgiaTuDf0mVGcdVUTUlOTefvV6fzjrvNCRxIJ4rQmlUkovI5vvk78Y2WRItCggYpDkQxQcSjZrnipYjx993z2ba1Mh/M+Dx1HJEPc/Vd3b+vuNaOPG6PrV7t7p3TaT3b3s7M/qcSbF+4bykX/OI2EPDsZ++4yLux1euhIIsEkJBjl6i1mxawqf36jZUuYMgXc0/+giAAqDiWQa245h+NPfJvZP17E43cPDR1HRCRX+u8tA7nugS7kKbiWLz/axhlnNQ8dSSS45i32sG9TRb6bk2Zu5YYNYcsWWKkeSyIHo+JQghk3+kzyFFvIHc+0YO2S1Yf+gIiI/L+/9e3Pv565mPzHLOCHycm0OLl+6EgiOcJ5HSIzDQ0Zs/SPlQ0aRB5nzcr+QCK5iIpDCaZk6WN44q45pGytwhnnTQ4dR0Qk17is24s8MegyCpedzrzvSlOvQbXQkURyjG5takLyNiZ/vvePlfWjJ09UHIoclIpDCer6v3eh0QnvMGtGD158aFjoOCIiOV6/815kwPBrKFHlcxb9UJsq1cqFjiSSo+RLTqJ4rZ9ZMKPsHyuLFYNKlVQcihyCikMJbuz77UgsvIKbH63P1l82h44jIpJjXdn9RQaOuIoSVT9l4Q8nUKr0MaEjieRIjU7Yxq6VNVm+bssfKxs0UHEocghBikMz625mc8ws1cyaHaTdUjObZWYzzGxqmvXFzWy8mS2IPuqvYy5WtmIp7r3hW/ZsPI6zu44IHUdEJEe6tseLvDLsKopX+YKfpzWjWLFCoSOJ5FhntS0KJPDGRwv+WNmgAcybB3v3HvBzIvEu1JXD2cB5QEbmMTjd3Ru7e9oi8g5gorvXBCZGX0su9q//Xki1uiP54uteDB84LnQcEZEc5aaLX+DFd6+kWOWvmD+tCcWLFw4dSXIYM+toZvPNbKGZ/eV3kZkVNbMPzOzH6An6viFyZpdLOtWChBTGTtr+x8oGDSKF4c8/hwsmksMFKQ7dfZ67zz+Kr+gCDIo+HwSce9ShJLiP32mI5d1MnzuPYc9vu0PHERHJEf7W+0WeGXolRSp8x/wpDShZokjoSJLDmFki8DxwJlAX6Glmdfdrdh0w190bAacBj5tZcrYGzUaljylIwcrzmT212B8rNWKpyCHl9HsOHfjEzKaZ2ZVp1pdx9zUA0cfSB/oCM7vSzKaa2dQNGzZkcVw5GrXrH8s1F3/MjnXN6dH9jdBxRESCu/3SF3li8BUULj+Vn6YcR+nSxUJHkpypObDQ3Re7+x7gbSIn0tNyoLCZGVAI2AikZG/M7FW7yQa2LKrN1h3RE87HHQdJSTBzZthgIjlYlhWHZjbBzGans+x/sDqYVu5+PJEzYdeZWevDzeHu/d29mbs3K1Wq1OF+XLLZsy/1olTVCYz4pCdfjvk2dBwRkWAevulVHnnzMgqVm86872tSrmzx0JEk56oArEjzemV0XVrPAXWA1cAs4CZ3T03vy2LlxHrb0/JBSn7eHh/tRpqcDLVr68qhyEFkWXHo7u3cvX46y6jD+I7V0cf1wAgiZ8YA1plZOYDo4/rMzi9hJCQmMGpASQDOvWYHqfvS/bslIhLT3nxoKHe82IN8JeYz++uqVChfInQkydksnXW+3+sOwAygPNAYeM7M0u2jHCsn1nufXR2AUeN//WOlRiwVOagc263UzAqaWeHfnwNnEBnIBmA00Cf6vA+Q4YJTcr6Wpzfm/LPe49flbbmu38DQcUREstX4QR/T5z/tSCzwC1+NzU+VKge8c0LkdyuBSmleVyRyhTCtvsD7HrEQWAIcl035gqh/bCnylFrC9O/y/7GyQQNYtgy2bg0XTCQHCzWVRVczWwm0BD4ys3HR9eXNbEy0WRngSzP7Efge+Mjdx0bfewhob2YLgPbR1xJD3n73EgqV/Z6X3+vMglkLQ8cREckWP477hk43Rn6vf/j2Ro4/vkbgRJJLTAFqmlm16CAzPYicSE9rOdAWwMzKALWBxdmaMoCqDVeybl4NUn7vifT7oDSzZx/4QyJxLNRopSPcvaK753X3Mu7eIbp+tbt3ij5f7O6Noks9d/9vms//6u5t3b1m9HFjiP2QrJOUJ4mBD2/Ddx3D2Rf9EDqOiEiWW/HDPFpcXJCU30ox4MU5dOx4fOhIkku4ewpwPTAOmAe86+5zzOxqM7s62ux+4CQzm0VkGrDb3f2XMImzT6tW4DtLMPbbJZEVGrFU5KBybLdSkW6929K8xdv8PLs7zz7wXug4IiJZZsuyVTQ8ewO7NtXhgfs+5dI+p4WOJLmMu49x91ruXv33E+ru/pK7vxR9vtrdz3D3BtExIAaHTZw9enSqCMC7H6+KrKhSBQoXVnEocgAqDiVH+2BYR5KKLOHvTzZg669bQscREcl0e7dso97p09m8ujXXXD+Su+48J3QkkZjR/oSqWKENfPVV9CevmQalETkIFYeSo5UuX5K7r/2OPRuPo2u34aHjiIhkrtRUTjntHVYtOYfO3YfywtPdQycSiSkJCUa5uotYPivNeD2/F4e+/4CuIqLiUHK8ex7sQaXaHzDpi4sYN+yz0HFERDLNlec+w3czLqfuicMZ+U6P0HFEYtIJLXaR8msVpv+8NrKiQQPYtAlW7z+gq4ioOJRcYfSgGpC0i563oLkPRSQmvHz7q7wy5hqKVf2CKZ+eg1l6U9WJyNHqekZkOpjBH2pQGpFDUXEouULjE+vQ85z32bTyVK7tq7kPRSR3+/qtj7nmuc7kKbKCKZNqUiB/cuhIIjGre9uakGcHkz7bE1mh4lDkgFQcSq7xxlu9KVT2e/oP68yCmQtCxxEROSLrZsyj7Y0V8H3JjBi6jRrVyoaOJBLTCuTLwzE15/PzD5EriBxzDFSooOJQJB0qDiXXSMqTxOuPROY+PKeX5j4UkdwnZfMWmpy7iF0b6/HAf77mrA5NQkcSiQt1m2zlt1U1+WXLzsgKjVgqki4Vh5KrnH9JW5o2f5v5sy6g/2MjQscREcm4ffs49fS3WLPsbM69aDh3/aNT6EQicaN1y/yQmsTwSYsiKxo0gLlzISUlbDCRHEbFoeQ6H753BgmFVnLTw8eya8dvoeOIiGTIrRc8xdczrqF281EMf0NTVohkpwvOqArA2M82RlY0aAB79sDPP4cLJZIDqTiUXKdspdLc2udTdv3SiJ4XDgkdR0TkkMY/+y5PfnglhSpMYerkM0lI0MikItmpcc0yJBRdww/TkyIrNCiNSLpUHEqu9PDTF1Oq6gRGjr+QqZ/PCB1HROSANs75mc73Hocl7mb86BIU0sikIkGUqbWM1T+Vj7yoUwcSE1UciuxHxaHkSgmJCbz9XCFITeK8y1eEjiMikr7duzmp6/fs2tiQe++dQovjjw2dSCRuNWiyi70bqrFk9WbImxdq1VJxKLIfFYeSa7U5qwVtWr/NigXn8OA/3wodR0TkL67q/BTzF/TipI7DuOcfZ4aOIxLX2p5cBIB3xqcZlEbFocifqDiUXG3Ee+eTp9gC/v3iCWzesCl0HBGR/zfq/gH0n3QDRat8zcRR54aOIxL3LmxfHYCJX26LrGjQAJYsgW3bAqYSyVlUHEquVqR4Ef5z/XT2bq7Jed01tYWI5Azrvp3OBU+2JCHvVj4dU5l8yUmhI4nEvSpli5Kn9BJmz8gfWfH7oDRz5oQLJZLDqDiUXO+O+y+kUu3RfPrlRUwa/XXoOCIS51K3bqPFRQvZs7kWjzw6jyZ1K4aOJCJRFeusYv3PVSIvNGKpyF+oOJSY8P4rVSFxLz1v3gbuoeOISBy7vOtzLF1yAW3OHcnfrjk9dBwRSaPJ8ftI3VqWqT+tgapVoVAhFYciaag4lJjQ7JSGnNX2PdYv6cDdtw4OHUdE4tRnTw9m4JfXU6zqV4x999zQcURkPx1aFwdg+ITlkJAA9evD7NmBU4nkHCoOJWa8+04P8hafy0MDW7Fx3cbQcUQkzuxdsYrOD1UBS+WjERXIk5QYOpKI7KdbmxqQsJfPvv4tsqJ+/ciVQ/U6EgFUHEoMKVC4AP+7aQ4pW47l/AtHhY4jIvHEnfPOG8LWtadw+fWfc1LjqqETiUg6ihfJT/4KC5k/s3BkRf368MsvsH592GAiOYSKQ4kpt97Tncq1RzH5y55M/vCb0HFEJE6MuP1ZPvzhJsrVm8DLj5wdOo6IHETVehvYuKgGKftSNSiNyH5UHErMGf5KNUjcS48bt4aOIiJxYPsPM7l4QGss32YmfNCQhAQLHUlEDqL5CQa7ijJx6rLIlUPQfYciUSoOJeY0O6Uhndq8x7olHbjnbxqcRkSy0J49dOg1md9+bcyd98yibrXSoROJyCGcfXoZAEZNWg2lS0cWXTkUAVQcSox6550LSS4+jwdfa8nmDZtCxxGRGNX/ykf4et511DxxDP/9R7vQcUQkA85udSzk2cHX36ZEVmjEUpH/p+JQYlKhIgX5342zSNlSna7dR4aOI9nEzI7PwNIgdE6JDevHTub6kT1ILLKCiaNbhY4jIhmULzmJIlUXsmhOZFoL6teHOXMgNTVsMJEcICl0AJGs8rd/X8DTb33A5K968OW4KZzc4YTQkSTrfQZMAQ5201c1oGq2pJHYtWsXZ9z4M3u3tObJgd9RqXTV0IlE5DDUbLCJaaNPZOeuvRRo0AB27IClS+HYY0NHEwlKVw4lpr33UiWwVC68XkNUx4kp7t7G3U8/0AIsDh1Scr9B1z7Gjwsup8FpY7j50pah44jIYTqpRTKk5GfU54s0KI1IGioOJaadeHpj2rV+l9ULz+KRu98OHUeymLu3yYw2Igeza/oMrhl9FgmF1/LRO6eEjiMiR+DcNhUA+Oiz9VCvXmSlBqURUbdSiX3vvXMepaov4u4XGnH9bTsoUKRg6EiSjczsXne/N3QOiRGpqfS47AN++/VubnlgApVKaxAayRxmdmsGmu1w95ezPEwcOK1JZazARqZ8b1C4MFStqiuHIujKocSBYiWKcke/b9mzsQ49e+rqYRzqHDqAxI7v/vMco+bcQqnjPuexO9uGjiOx5TagEFD4IMvfgqWLMQkJRvHqi1k+Lzr9TP36unIogq4cSpy477GLeHnEBEZP6M7Mb+fQsEW90JEk+2hGcskcK1fS9bUq4AkMG1pFk91LZnvT3f9zsAZmpq4vmahu4+18MbgJ6zftoHSDBjB2LOzZA8nJoaOJBBPkyqGZdTezOWaWambNDtCmtpnNSLNsNbObo+/da2ar0rzXKVt3QHKdhATjzcfzQUo+zrv8p9BxJHs1DR1AYoA7/7rkadas7EKHiz6ldeMqoRNJjHH3f2RGG8m4U1sWAE9k2KTooDQpKfDzz6FjiQQVqlvpbOA84PMDNXD3+e7e2N0bE/lxtxMYkabJk7+/7+5jsjStxIQO553MiS3eZtGc8+n/+MjQcSQLmNmV+69z99RDtRE5lHWD3+OhaTeQt/Q8hr/cPnQciRNm9mHoDLGs+xlVAfjk843QIDoFrrqWSpw7YLdSMzsvA5/fdSSFmbvPi24jox9pCyxy92WHuy2RtEYO7UCFequ45ZHKXHrdHpLzqetIjLnDzH45yPsG3AT0z6Y8Egs2b+bs+9ayb1tlnh0wnYI6bkj2qRA6QCxrWL00icVWMX1qMjxaG5KSNCiNxL2D3XP4CjCKg9+v0xrIjqt2PYCh+6273sx6A1OBv7n7pmzIIblc2cpluPbC13nutUu5vM8A3ninX+hIkrk+A845RJvx2RFEYse71/yPqYv/R+1Tx3FNtw6h40h8+SF0gFhXtvZy1syvELnPsFYtXTmUuHew4vBjdz/oL2czG3yQ9yYAZdN56y53H5XBfJhZMpHRBu9Ms/pF4H7Ao4+PA+lmjXYhuxKgcuXKGd2sxLCnX+7NGx9/zeAPz+b+n5dTpZb+XcQKd+8bOoPEltQfZnD5pHOxAr/y0ZATQseROHOo32Fy9Boev5uPv6vC/OW/UrtBA/j++9CRRII64D2H7t7rUB8+WBt3b+fu9dNZMlwYRp0JTHf3dWm+e52774veS/QK0PwgOfq7ezN3b1aqVKnD3LTEooTEBF66dyu+syTn9Trgba8iEu/c+fu1b7Jt/Ulcct1MqlcoHjqRxDAzO2R394y0kcNzRutiALzzyeLIoDRLlsC2bWFDiQR0WAPSBDoo9WS/LqVmVi7Ny65EBrgRybCeV3TkuIbDmD7tQj58e3LoOCKSA219ZzjPzr2OfGXm8soDbULHkdh3rpk9c5DlWeD00CFjTY8zaoDt49OvdvwxKM3cuWFDiQR0uKOVpjvtxOEys65mthJoCXxkZuOi68ub2Zg07QoA7YH39/uKR8xslpnNJHKgvCUzckl8ef/1xpBnB5fesY/UfamHbC8icWT3bnrdN5OUrcdyz/+2k5wnMXQiiX23AdMOskwF7gqWLkaVLV6IvOUWM+eHgpErh6BBaSSuHeyew/Ssz4yNuvsI/jwtxe/rVwOd0rzeCZRIp90lmZFD4ludJrU4t/1rjPzwMv5922Duf+KQPaklhzOzWw/y9m5gEfDJ/tNbiOxvwX+f5sMlt1Kyzlfc2a9V6DgSB9x9UOgM8apynbUs/KYuqVWKkVCggAalkbh2WFcO3b1jVgURCWHIWz1JLv4TDw1owdaNW0PHkaNX+CBLWSKDU70dLJ3kDuvWccHg/PiewvR/Ufeqi8S6E5o7vrMEn89aBfXq6cqhxLUDFodmdu+hPpyRNiI5WYHCBbjnqh9I2VKDC3sOCx1Hjt4C4Dl3vy+d5W53Pw84LnRIydnG3/Q/Ziy7muNOn0TXU2uFjiMiWezs00oD8P6ElZGupbpyKHHsYN1KLzezg11KMSLzD96bqYlEstld/+vJs++MY+zk7sz4ZjaNW9YPHUmOXGXgPTPLA0wEPga+d3f/vYG7NwwVTnKBH37g0i9Pg6RdvNNf/1RE4kGX1tUh6Te++mYPnNIABg6E9euhdOnQ0USy3cG6lb7CwbtoFYq2Ecn1hjxRGPYl0+2K+aGjyFFw94fcvQ2Re5d/JDL/6XQze8vMeptZmbAJJUdz59nrnmP1qq606fkNDavrn4tkPzNrZmYjzGy6mc1MMwCfZJEC+fJQuOoCFswq/segNHPmhA0lEsgBrxy6+33ZGUQkpLZdTuLE5m/w3Te9ee3p0Vx2U+fQkeQouPs2IoNejQAws7pE5kx9A+gQMJrkYCnvDePOBVeTUGQlbz+lQWgkmCFERi6dBWjwrGxSs8FGpn9wIjtrFqMARLqWnq6ZQyT+HO5UFiIx6/232pNQcC03P1RWU1vkcmbWOu0ClASmuLsKQ0nfnj3cfO8kdvxyApf9bRGlihUMnUji1wZ3H+3uS9x92e/LoT5kZh3NbL6ZLTSzOw7Q5jQzm2Fmc8zss8yPnnu1apkMKfkZ9fMOKFFCg9JI3FJxKBJVvmo5encew/a1zbnhCo0onsvdlma5G/iAo7w/2syKm9l4M1sQfTzmAO2KmdkwM/vJzOaZWcuj2a5kj+0vv8pLK+8gf4XZPP/Pk0PHkfj2bzN71cx6mtl5vy8H+4CZJQLPE+khURfoGe0xkbZNMeAFoLO71wO6Z0383Om8dhUB+HDyBg1KI3FNxaFIGq8M6k3+0j/w8rB2rF+ZKdN6SgDufk6apT1QH1h3lF97BzDR3WsSGewm3TPzwNPAWHc/DmgEzDvK7UpW27WLa55fwL5tVbjzvh3kSdKE9xJUX6Ax0BE4J7qcfYjPNAcWuvtid99DZMqeLvu1uQh4392XA7i7/sil0bpRJazAr0z53qBBg8iVwz/GMhOJGwe859DMngUO9P+K3yeTHhK9t0ckJiTlSeKRW5Zww51NuKDXQCZP7hs6kmSOlUQKxKPRBTgt+nwQMBm4PW0DMysCtAYuBYj+SNtzlNuVLLbjuZcYuvpmClaZwV19m4eOI9LI3Rsc5mcqACvSvF4JnLhfm1pAHjObTGRgwafd/Y0jThljEhKMkjUXs3xeWeheH7Zvh2XLoGrV0NFEstXBprKYeojP1QPeB9pnaiKRwK6/4zweGfgBn311Id9NnM6JbY8PHUkO034ntxKInIX/8Si/toy7rwFw9zVmlt4Y58cCG4CBZtYImAbc5O47jnLbklV27ODaVxazb1sV/vHEOhISLHQikW/NrK67zz2Mz6T3D3f/E/xJQFOgLZAf+MbMvnX3n//yZWZXAlcCVK5c+TBi5G51G+/gs0FNWV15BeUh0rVUxaHEmYN1K60AzHT3Qeksr7n7rcDebMopkq2GPl8eMC68bsUh20qONJVIYTYN+Aa43d17HepDZjbBzGans+zfPetAkoDjgRfdvQmwgwN0PzWzK81sqplN3bBhQwa/XjLbjqee4601t1Kw8o/8q98JoeOIAJwMzIgOLpPRqSxWApXSvK4IrE6nzVh33+HuvwCfE+n6/hfu3t/dm7l7s1KlSh3hbuQ+bU4uCCTw7roCYAYzZoSOJJLtDnblcDFwU/Ts949EJpP+xN03/d7A3TtlcT6RIFq1a8rJLV/nyy8u5eXHRnLV388NHUkOg7sf0YhC7t7uQO+Z2TozKxe9algOSO9+nZXASnf/Lvp6GAcoDt29P9AfoFmzZrqxJYStW7n+9eWkbKvKXY9/p6uGklN0PILPTAFqmlk1YBXQg8g9hmmNAp4zsyQgmUi30yePJmisufCMY/k38Mn3u7i5Vi2YerBOdCKx6YBXDt39bXe/NHr2+2ki3aXeN7PPzeweM9ONGRLThr/ViYRCK7n1sUqk7E0JHUfCGw30iT7vQ+SH1p+4+1pghZnVjq5qCxxO1zDJRr89/hSD1/+dApVncs9l+pMmOUPa6SsyOpWFu6cA1wPjiAyC9a67zzGzq83s6mibecBYYCbwPfCqu2u+hjRqVy5BUollzJyeF5o2hWnTQkcSyXYZGq3U3X9w9wfd/XQiI2bNAS7P0mQigZWuWJorzx/PznVNuaaf7tkXHgLam9kCIvdaPwRgZuXNbEyadjcAQ6LdwBoD/8vuoJIBmzZx3ZDVpGytxt/+uVNXDSXXc/cx7l7L3au7+3+j615y95fStHnU3eu6e313fypY2Bys/HGrWDu/cqQ4XLUK1h3tQNciucshi8O0c+xE59lpB2wG/pXV4URCe/aVSyhY9nsGjOjI2uX6A5FbmNnRjkz6F+7+q7u3dfea0ceN0fWr03axd/cZ0Xt1Grr7uWm74kvO8dsjj/Pmhn9QoNIs7r1i/0EdRSReNWm6h32bKzCj3HGRFbp6KHEmI1cOLwNeBS6OLq8AtwJfmdklWZhNJLikPEk89Y81pO4ozwW9Pg4dRzLuJTP73syujU78LPKHDRu44Z11pGw9llvu3KGrhiLy/zq0Lg7Ae+sKRVaoOJQ4k5HiMBWo4+7nu/v5QF0i8xyeyH5zfInEostv6UKV2qP44psL+WaC/kjkBu5+MpGTWZWAqWb2lplp2h0BYNdDj/DGxtvJX2k2/7lKVw0lZ4n21FpgZlvMbKuZbTOzraFzxYsL29eEhBQmT0mFWrVUHErcyUhxWNXd0/anWw/Uinap0lQWEheGPlcRgJ43aGqL3MLdFxDp/n47cCrwjJn9FO0eL/Fq/Xpuen8je7fU4OY7tuuqoeREjwCd3b2ouxdx98LuXiR0qHhRvEh+8ldYwE8/FtagNBKXMlIcfmFmH5pZHzP7fYS+z82sIJF7D0ViXst2TTm55Tss++lcXn3yL4NUSg5jZg3N7Ekio/a1Ac5x9zrR5xq6PY7tfuY5Bm66g/wV5/LA1bpqKDnSuujIohJI1Xob2LioBimNm8DKlbA+vZmLRGJTRorD64CBREbdawK8AVwXnUT19CzMJpKjvDf4TBIKrubmR8qRmrIvdBw5uOeAH4BG7n6du0+HyOAxaDCt+LVjB/e8vYi9W2py7d8266qh5FRTzewdM+u534CAkk1ObG6wqyjji9WMrNDVQ4kjhywOPWK4u9/i7je7+zB314TNEnfKVi5Dn3PHsmNtc264UlNb5GTu3trd33D339J5780QmSQHGDCAl7ZdQ1LxpfzvOl01lByrCLATOAM4J7qcHTRRnOncpiwA76+I9uZVcShxJCl0AJHcpP/A3rz9yQxeHtaW+/77CyXLlQwdSdJhZrOA/U9ibQGmAg+4+6/Zn0qCSknhzefHsHX9DXS56VOS81QNnUgkXe7eN3SGeHfOydWxfJv54rskDUojcScj3UpFJCopTxIP37KIfdsqc+HFH4SOIwf2MfARf0zB8wHwObAWeD1cLAlm2DDu2XEJ5N3Cc3cdHzqNyAGZWUUzG2Fm681snZkNN7OKoXPFk6TEBEod9zNLZlbQoDQSdw5ZHJrZ2WamIlIk6oY7z6d8jY+Y9OUFTP9iZug4kr5W7n6nu8+KLncBp7n7w0DVwNkku7kz7dHXWLq6O006TadiqaKhE4kczEBgNFAeqEDk5NbAoIni0PEtdrJnXXXm1WoCK1ZoUBqJGxkp+noAC8zsETOrk9WBRHKDN58sDqlJXHDNgtBRJH2FzOz/byozs+ZAdEZjUsJEkmA+/ZQbf20DnsCT9xwbOo3IoZRy94HunhJdXgdKhQ4Vbzq3Kw7A6zuqRFbo6qHEiYwMSNOLyCili4CBZvaNmV1pZoWzPJ1IDtXm7JaccMK7LJrTlSEvjQkdR/7qMuBVM1tiZkuAV4HLo1PwPBg2mmS3Xx56km/WXUX5E77l1MZVQscROZRfzKyXmSVGl16A7pPOZhd3rA1Ju5i4oERkhYpDiRMZ6i7q7luB4cDbQDmgKzDdzG7IwmwiOdp7g07F8v/K9f8tTOq+1NBxJMrMEoFT3L0B0Sl43L2hu0+JTsHzbtiEkq1mzuSmnyriu4rzz38UDJ1GJCP6ARcQuUd6DdAtuk6yUZGCeSly7Hx++rEs1Kyp4lDiRkbuOTzHzEYAk4A8QHN3PxNoBPw9i/OJ5FhValXmgo6j2bzyFP71tyGh40iUu+8DukSfb3H3zWETSUj7Hn2MYdtupkCV2VzTtVHoOCKH5O7L3b2zu5dy99Lufq67LwudKx7Vb7aJHctrs7ZRcxWHEjcycuWwO/Bk9Mz7o+6+HsDdd6IzWRLnXh98McnF5/HYoBZs37I9dBz5w1dm9pyZnWJmx/++hA4l2WzFCv43aRN7Ntemz9WbNOm95Ghm9o/o47Nm9sz+S+h88ejMtoUgNYk3izWNDEqzYUPoSCJZLiP3HPZ2988P8N7EzI8kknvkK5CPf142nb2ba9L7kndCx5E/nATUA/4DPB5dHguaSLLfU0/xdOqNJBRdxaM3a9J7yfHmRR+nAtPSWSSb9T6rJpDKR+urRlbo6qHEgaRDNTCzbRx4Mum/ufvirAgmklvc/WBPnnnnU0aO78qiOUuoXq9a6Ehxz91PD51BAtu8mVFvfcGvax+n/VUTKZivQuhEIgfl7r9PnrvT3d9L+56ZdQ8QKe5VLlOU/JV+YuaCypEV06ZBx45hQ4lksYx0K30CuI3IXDsVidxn+AqRwWkGZF00kdwhITGB/v9JwXcXpfulX4eOI4CZlTGz18zs4+jrumZ2Wehcko3eeIPb81wFeXbw/L8ah04jcjjuzOA6yQY1m6xj08Lj2Fmjtq4cSlzISHHY0d1fdvdt7r7V3fsDndz9HeCYI9momT1qZj+Z2UwzG2FmxQ7QrqOZzTezhWZ2R5r1xc1svJktiD4eUQ6RzHJ+n/Yc12A4P0y/gHHDvwgdR+B1YByRSaQBfgZuDhVGspk7P738FvNXX0yddlOoWbFE6EQih2RmZ5rZs0CF/e43fB3NzxpMm9Z5YG9B3qndXsWhxIWMFIepZnaBmSVElwvSvLd/d9OMGg/Ud/eGRH60/eWMWHQ4+ueBM4G6QE8zqxt9+w5gorvXBCZGX4sE9c4rDSBpF31v18A0OUDJ6JQVqQDungLsCxtJss0XX3Db9lNgXz4evLP8oduL5Ayridyys4s/32s4GugQMFdc69O5OgAjvREsXw6//BI4kUjWykhxeDFwCbAeWBd93svM8gPXH8lG3f2T6I81gG+JdFfdX3Ngobsvdvc9RLqxdom+1wUYFH0+CDj3SHKIZKaGzevQvvUw1iw6k6ceGBY6TrzbYWYliJ7AMrMWRO6VljiQ+uKLfLL1MgrXmEGXU2qFjiOSIe7+o7sPAmq4+6A0y/vuvil0vnjVuGYZkkouY8rK6LFEVw8lxh20OIxevbvG3c9x95LROXfOcfeF7v6bu3+ZCRn6AR+ns74CsCLN65XRdQBl3H0NQPSxdCbkEDlqbw/uQmLh5fzz2eqk7FUvoIBuJXK2vbqZfQW8AdwQNpJki/XreeXT1ezZfBznX7w5dBqRI1HVzIaZ2VwzW/z7EjpUPKvSYAVrF9WNdEWZOjV0HJEsddDiMDqZdNMj+WIzm2Bms9NZuqRpcxeRfvTpzSCe3oRUh92N1cyuNLOpZjZ1g+ankSxWvExxrjh/Ir+tb8L1l78ZOk7ccvfpwKlEprS4Cqjn7jPDppJs8frrPJ6nL+TdyoM3HtGfL5HQBgIvEvl9dDqRk1v6gxLQySc7vqMkY+ueriuHEvMy0q30BzMbbWaXmNl5vy+H+pC7t3P3+uksowDMrA9wNnCxu6dX9K0EKqV5XZFIf3yAdWZWLvo95Yh0eT1Qjv7u3szdm5UqVSoDuytydJ7tfwkFykzj1RHt2Lj219Bx4llzoBFwPJF7lnsHziNZLTWVVa8MYcG67tQ6dTplixcOnUjkSOSPziNt7r7M3e8F2gTOFNd6nhW5++ntkm1VHErMy0hxWBz4lciB6ZzocvbRbNTMOgK3A53dfecBmk0BappZNTNLBnoQ6SZG9LFP9HkfYNTR5BHJTEl5knjopqXs21aJHr30TzMEM3uTyKT3JwMnRJdmQUNJ1ps4kX/uagF7C3Lb9RrEWnKtXWaWACwws+vNrCu6fSao9idUxQpt4Kst0UFp1h/wmoRIrpd0qAbu3jcLtvsckBcYb2YA37r71WZWHnjV3Tu5e4qZXU9kOPpEYIC7z4l+/iHg3ei8ZcsBTQ4rOcoNd57Pg699zPgvujP7+3nUb14ndKR40wyoe4BeCRKrXnqJ4bvvJG/5+fQ7q2HoNCJH6magAHAjcD+Rk/N9DvYByVoJCUa5uotYvrBBZMUXX8D554cNJZJFDnnl0MxqmdlEM5sdfd3QzP51NBt19xruXsndG0eXq6PrV7t7pzTtxrh7LXev7u7/TbP+V3dv6+41o48bjyaPSFYY+HAhSMnHhVf+GDpKPJoNlA0dQrLR6tWMmLyYHRua0fGClSQkpHfbukjO5+5T3H27u68ELgMudfdvQ+eKdye02EXKxipMLXEsTJ4cOo5IlslIt9JXiMxDuBcgOqhDj6wMJRILOpx/Cg2bvMfcmd0Z/dak0HHiTUlgrpmNi94zPdrMRh/yU5J7vfYa9xfsC0m7ePhWXTWU3MvM3jKzImZWEJgLzDez20LninfndYz07H2z9nnw2WeB04hknYwUhwXc/fv91mmMfpEMePe15pC8jcvv0vzr2exeIvOf/g94PM0isSglha2vvM6MDb2o2HwqtStp8DHJ1eq6+1Yix7AxQGUic0xLQBe0rQXJ2/k0pTnMmgW//BI6kkiWyEhx+IuZVeePyaS7AWuyNJVIjKjdqAZntRnOhqXteeSed0LHiRvu/hmwFMgTfT4FmB40lGSdjz/m3/ua47uKc82VeUKnETlaecwsD5HicJS77+UIpvKSzJUvOYniNeezYEX9yIovvggbSCSLZKQ4vA54GTjOzFYRuVH6mqwMJRJL3hrcnaSii/n3i3XYs2tP6DhxwcyuAIYROXYBVABGBgskWeull3jTrySx+DJuu1iD0kqu9zKRk1sFgc/NrAqwNWgiAaBR823sWlObJUXK6r5DiVmHLA7dfbG7twNKAce5+8nuvjTLk4nEiCLFi3D9hV+w65eGXHv54NBx4sV1QCuiP6jcfQEaCj42LVvGV5/P49c1p3NylwXkSUoMnUjkqLj7M+5eITpyu7v7MuD00LkEup9VHEjgufo9dd+hxKxDTmVhZnmB84GqQFJ06gnc/T9Zmkwkhjz+wiW8Mvp7Bo7qwENrfqFkuZKhI8W63e6+5/fjlZkloW5Zsen117nrmH6wYx//vbVm6DQiR8zMern7YDO79QBNnsjWQPIXfc+py7X5tvDhb615fOZTsHEjFC8eOpZIpspIt9JRQBcig9DsSLOISAYlJCbwyC0rSd1egR69PggdJx58Zmb/BPKbWXvgPUD/4WONO3sGD+XLTX0p0XAKrepXCZ1I5GgUjD4WTmcpFCqU/CFfchIVGs9l4aLmpLrrvkOJSYe8cghUdPeOWZ5EJMZd+4/z+O8rHzHxy27M/n4e9ZvXCR0plt1BZH6wWcBVREb8ezVoIsl8U6bw+LZj2be9Apf0WRw6jchRcfff75Ge4O5fpX3PzFoFiCTpaH9GCq9/W54RZZpw/uTJ0KVL6EgimSojVw6/NrMGWZ5EJA68/kgxSMnHhVf+GDpKTHP3VHd/xd27u3u36HN1K401Q4bwWvLFkH8j916pgWgkZjybwXUSwHU9awAwoMIFuu9QYlJGrhyeDFxqZkuA3YAB7u6aZVjkMLXv2oqGTd5i5vQLGD10Ep17tgkdSSR3Sklh8zsjWbTpAWqcOpWiBTVeh+RuZtYSOAkotd99h0UAjbSUQzQ7rhz5yv/MN5taw9J/wqZNcMwxoWOJZJqMXDk8E6gJnAGcA5wdfRSRI/Dua80heQdX3JUSOopI7jVhAg8lNIc9henXq0DoNCKZIZnIvYVJ/Pl+w61At4C5ZD8NWq1m0/JmrE0qAF9+GTqOSKbKyFQWy4BKQJvo850Z+ZyIpK92oxp0PHU465ecweP3vRs6jkjuNGQIbyVehBVax809moROI3LU3P0zd78PaBF9fMzd73P3J6LT8UgO0ePcorAvmWfLttd8hxJzDlnkmdm/gduBO6Or8gCarE3kKAx541wSCy/j7udrkbJXVxAzi5l9YGajD7SEzieZZMcO1o0ez4o1Z1Ln1NnkT04OnUgkM5U3s7nAPAAza2RmLwTOJGlc2bUeJG9nVKFzVRxKzMnIPYddgSbAdAB3X21mhbM0lUiMK16mOFecP4qXXu/LDVcM5MXX+4aOFCseCx1AssGoUfy3UHvYmo+r+hQNnUYksz0FdABGA7j7j2bWOmgi+ZNC+ZMp2+AH5i86jdQtfUnYvBmKFQsdSyRTZKR76J7oKH8OYGYFD9FeRDLg2f6XkL/0D7zyfls2rtsYOk5MiHbLOuASOp9kkiFDeM8uJqHYCq7u2jh0GpFM5+4r9lu1L0gQOaDT2u4mZXMVxhetofsOJaZkpDh818xeBoqZ2RXABOCVrI0lEvuS8iTxwHUL2betMhf3Hhk6Tkwxs5pmNszM5prZ4t+X0LkkE6xfz9JJ37N2TVsatf2J5KSMdIARyVVWmNlJgJtZspn9nWgXU8k5rulZDYD+x3TSlBYSUzIyIM1jwDBgOFAbuMfdNd+OSCa49Z7ulK42nrGfnc/8HxeGjhNLBgIvAinA6cAbwJtBE0nmePddHihxLqTm4bq+JUOnEckKVwPXARWAlUDj6GvJQVo3rkRy6cV8vu9s3XcoMSVDo466+3h3v83d/+7u47M6lEg86X9/EuwpRI/LvwsdJZbkd/eJgLn7Mne/F9CkkrFg8GBG+cUkllzMpWc2Cp1GJNO5+y/ufrG7l3H30u7ey91/DZ1L/qpuyxX8sroVG3+cA1u3ho4jkik0JYVIYF0uPp3aDYYzY3p3Jn3wTeg4sWKXmSUAC8zsejPrCpQOHUqO0sKFzJuxjF9Wt+aEDgtJTNCfMIk9ZlbKzP5pZv3NbMDvSwY+19HM5pvZQjO74yDtTjCzfWamuROP0vmdC0JKfl4seYruO5SYob+sIjnA2y/VhcS99PnbL6GjxIqbgQLAjUBT4BKgT8hAkgneeov7S3YHErj5snKh04hklVFAUSJjPHyUZjkgM0sEngfOBOoCPc2s7gHaPQyMy+TMcenabvUg6TeG5z1T9x1KzMhQcWhm+c2sdlaHEYlXjVvWp3XL91i54BxefXJU6Di5nrtPcfft7r7S3fu6+3nu/m3oXHIU3GHwYMbs60Vy+XlccFr90IlEskoBd7/d3d919+G/L4f4THNgobsvdvc9wNtAl3Ta3UBkDIn1mZw5LhUvkp+SdWczZ2tHmDQpdByRTHHI4tDMzgFmAGOjrxtrMmmRzDf0jY4kFFzL3x4rReq+1NBxcjUzq2Vmr5jZJ2Y26fcldC45ClOnMm3FHrasbU7Ljssxs9CJRLLKh2bW6TA/UwFIO/3Fyui6/2dmFYjMXf3S0cWTtE5ps5M9G4/jy/kbYNWq0HFEjlpGrhzeS+SM1GYAd58BVM2qQCLxqnyVslx45hi2rj6Ju24ZEjpObvceMB34F3BbmkVyq6FD+W+pHgD846oqgcOIZKmbiBSIv5nZVjPbZmaHGu0kvbMlvt/rp4Db3f2Qcyaa2ZVmNtXMpm7YsCFjqePUFRdWAuCFYzrCKPX8kdwvI8VhirtvyfIkIsKAQReRXPwnHh98Iju37QwdJzdLcfcX3f17d5/2+xI6lBwhdxgxgvG7e5Kvyo+ceYLucpDY5e6F3T3B3fO7e5Ho6yKH+NhKoFKa1xWB1fu1aQa8bWZLgW7AC2Z27gEy9Hf3Zu7erFSpUke2I3GiQ/NqJBVfwaecDSNGhI4jctQyUhzONrOLgMToxNLPAl9ncS6RuJSvQD7+3nsKezfVol+foaHj5GYfmNm1ZlbOzIr/voQOJUdo5kw+35DM9vWNOPWstepSKjHNzI5PZ6luZkkH+dgUoKaZVTOzZKAH8KdbgNy9mrtXdfeqROavvtbdR2bVfsSLhASjzkmLWbu6Dau/mgIbN4aOJHJUMlIc3gDUA3YDbwFbiIwEKCJZ4P7HLqZoha94d9zZrF66JnSc3KoPkW6kXwPTosvUoInkyI0YwUPHXACkcsdV1UOnEclqLwDfAq9El2+JDDDzs5mdkd4H3D0FuJ7IKKTzgHfdfY6ZXW1mV2dP7Ph1zaXHQEoBHihxFnz4Yeg4IkflkMWhu+9097uA09z9BHf/l7vvyoZsInEpITGBJ/+xEd9Zhp59NNr44YrOb3hH9Cx52uXY0NnkCI0cyRe7u5G/2kxOa1gjdBqRrLYUaOLuTd29KdAYmA20Ax450IfcfYy713L36u7+3+i6l9z9LwPQuPul7j4sS9LHoSvOrU9isVUMt4vVtVRyvYyMVnqSmc0lciYKM2tkZi9keTKRONb3xnOoVOsDPv+mOzO+mR06Tq7i7qnAdaFzSCZZsoQff9rI9g2NaNF2beg0ItnhOHef8/sLd59LpFhcHDCTHERSYgLHt1/A+lXtWfDp97BjR+hIIkcsI91KnwQ6AL8CuPuPQOusDCUi8PpjJWFfHnpePefQjWV/483s72ZWSfcc5nIjR/JEic4AXN2rfOAwItlivpm9aGanRpcXiHQpzQvsDR1O0nfL5WUhNQ/3FT0LxqnXj+ReGSkOcfcV+6065DDIInJ02pzTkkZNhvHTrG58+Pbk0HFym35Erh5+ju45zN1GjmQs55NYaiHdWmvie4kLlwILiYzvcAuwOLpuL3B6qFBycBe2q01y6cV8tPcidS2VXC0jxeEKMzsJcDNLNrO/E+1iKiJZa+grzSDPDq64S7f5Ho507jfUPYe50YYNrPt6JuvXnky9kxeRYBk6nymSq7n7b+7+uLt3dfdz3f2x6PgPqe6+PXQ+SV9CgnFSp+VsXtOaqROmwJ49oSOJHJGM/KW9msgZ+ApE5tFpjO7nEckWdZrUov0pw1i7uCPPPTQ8dJxcw8wKmNm/zKx/9HVNMzs7dC45TB98wJPFO0JqHi7uXih0GpFsET1eDTOzuWa2+PcldC45tDuurgIk8EDeM2Hy5NBxRI7IQYtDM0sEnnL3i929jLuXdvde7v5rNuUTiXtvvdmZhEKruOPpSqTuSw0dJ7cYCOwBToq+Xgk8EC6OHJGRIxme93wouIFrux4fOo1IdhkIvAikEOlG+gbwZtBEkiEdTqxG/kpzmbCzp7qWSq510OLQ3fcBpaITqmYaM3vUzH4ys5lmNsLMiqXTppKZfWpm88xsjpndlOa9e81slZnNiC6dMjOfSE5SslxJep8zjh1rm3Pb9fp9kEHV3f0RooM3uPtvgGZOz022b2f3uEksWt+eyif8SKF8+UMnEsku+d19ImDuvszd7wXaBM4kGdSm83p2bGjOxHHTIVUndCX3yUi30qXAV2Z2t5nd+vtylNsdD9R394bAz8Cd6bRJAf7m7nWAFsB1ZlY3zftPunvj6DLmKPOI5GgvD+hF3hJzePbtk9i+RbecZMAeM8sPOICZVQd2h40kh2XcOF4tdBK+uyhduoQOI5KtdkXna11gZtebWVegdOhQkjF3X1MLgAdT28N33wVOI3L4MlIcrgY+jLYtnGY5Yu7+ibunRF9+C1RMp80ad58efb6NyCA4FY5muyK5VXK+ZO7sN4O9m2vSt887oePkBvcCY4FKZjYEmAjcfjRfGJ0OY7yZLYg+HnOAdrdEezvMNrOhZpbvaLYbt0aO5PVC50HSTm69uGHoNCLZ6WagAHAj0BS4BOgTMpBk3In1ylOk5gy+3NqT1OEaK0BynwMWh2b2e/+1ze5+3/5LJmboB3x8sAZmVhVoAqQ9BXN9tFvqgAP9SIt+9kozm2pmUzds2JApgUVCuPvBnhSr+AXDx5/DysWrQsfJ0dz9E+A8IsO/DwWaufunR/m1dwAT3b0mkWLzjv0bmFkFIj/omrl7fSAR6HGU240/e/fiH3zIj5vOoXj9qVQtpYsmEj/cfYq7b3f3le7e193Pc/dvQ+eSjOvUdSu7N9Xj/TE/gnvoOCKH5WBXDpuaWRWgn5kdk3Yi6YxMJm1mE6JnzvdfuqRpcxeR7qNDDvI9hYDhwM3uvjW6+kWgOpGRU9cAjx/o8+7e392buXuzUqVKHSq2SI6VkJjAU7dvwXeW5qJLx4eOk6OZ2UR3/9XdP3L3D939FzObeJRf2wUYFH0+CDj3AO2SgPxmlkTk7P/qo9xu/PnsM8akHMvebRVoc6a6UUt8MbNm0fEYpkdPgs80s5mhc0nG3XNNXbAUHt/ZBmbPDh1H5LAkHeS9l4h0yzqWyATSaQdz8Oj6A3L3dgd738z6AGcDbd3TP61iZnmIFIZD3P39NN+9Lk2bV4h0exWJeX2uP5u7n/2AL77tzoxvZtO4pSYFTyvahbMAUDLao+D341YRoPxRfn0Zd18DkW7vZvaXy1nuvsrMHgOWA78Bn0SvYqaX9UrgSoDKlSsfZbQYM3Ikzx/TFXbu45beNUOnEcluQ4DbgFmARjTJhepULUmJut8zdUUPUt8bSEKDBqEjiWTYAa8cuvsz0cFgBrj7sZk5mbSZdSRy/09nd995gDYGvAbMc/cn9nuvXJqXXQGdlpG4MfCRErAvmYuumRM6Sk50FZGTWcdFH39fRgHPH+rDGenxcIjPH0PkCmM1IsVoQTPrlV5b9Wo4gNRUGDmSL3efS75jp9Oydo3QiUSy2wZ3H+3uS6KjlS5z92WhQ8nh6XrBHlK2VmPAu1Nh377QcUQy7JAD0rj7NVmw3eeIDGozPjoVxUsAZlbezH4febQVkZuw26QzZcUjZjYr2s3idOCWLMgokiO17XISDZu8x7yZ5/PxsM9Dx8lR3P1pd68G/H2/k1qN3P25DHy+nbvXT2cZBaz7/cRU9HF9Ol/RDlji7hvcfS/wPn/MtSgZMW0aszYksW1DfU5ss5bIeUKRuPJvM3vVzHqa2Xm/L6FDyeH59zUNIWkXT/12NozRoPqSexysW2mWcfd0TwW7+2qgU/T5lxxgXjJ3vyTr0onkfG+/0oy6LXZy2R3bWd0tdJqcx92fNbOTgKqkOc65+xtH8bWjiYwY+FD0cVQ6bZYDLcysAJFupW2BqUexzfgzYgRPlOgCa+Cqi8uGTiMSQl8ivR/y8Ee3UidysklyiYqlilD95C+Y81Vvlr1wCVXOOSd0JJEMychUFiKSw9RpUou2Jw9nzaJOvPzYyNBxcpzoaMuPAScDJ0SXZkf5tQ8B7c1sAdA++vpPPR7c/TtgGDCdyP1CCUD/o9xufPnoIz5O6EZi6Z/ofkqT0GlEQmgU7XLeJzpaaV937xc6lBy+++8sBXsLceO8arBkSeg4Ihmi4lAkl3pr0FkkFFzNbU+UIXWfxizYTzOglbtf6+43RJcbj+YLo6OftnX3mtHHjdH1q929U5p2/3b346LdUS9x991HuS/xY/Vq1s9azrrVLal7ygKSEoJ0bhEJ7Vszqxs6hBy9nmccR5EaPzBm0w3seVnnCSV3UHEokkuVrlianp0+ZtualvzrbwecDSZezQbUJzG3+eQTni7RCTyJi7oVCp1GJJSTgRlmNj86jcXvYyxILnTldbtI2VqNu99fAXv2hI4jckgqDkVysVdfv5jk4j/x+JsnsGvnrtBxcpKSwFwzG2dmo39fQoeSQxg7lmF5u0KhNVx77gmh04iE0hGoCZwBnENk2i/dsJZL3X/tCSQWW8HLWy+HESNCxxE5JBWHIrlYvgL5uPXiKezZeBxX9n0rdJyc5F4ik9T/D3g8zSI51b59pIybyMJf2lHh+JkUyacrhxKf0k5foakscr98yUl06LmILetO453n0hvHTCRnUXEoksv998mLKVzuW4Z81JFf1vwSOk6O4O6fpbeEziUHMXUqw1OOI3V3Mc7omBI6jYhIpnnuX40hzw7uXt4O5s0LHUfkoFQciuRyCYkJPHzLGlJ3lOeiS+K756SZfRl93GZmW9Ms28xsa+h8chBjxzKw8JmQsJdrL6wVOo2ISKapVr4YdU/7jgWrLmLeU6+EjiNyUCoORWLANbd1peyxYxn/5fnM/3Fh6DjBuPvJ0cfC7l4kzVLY3YuEzicHMW4c3+w9m7zVptO0WrpT4YqI5FoP/7sK7MvHDZOKwM6doeOIHJCKQ5EY8fIDybCnED0v/y50FJHDs3Ejy75fytZfGtL45DWYWehEIiKZ6uxW1Slx3Ld8uuZKtg8eGjqOyAGpOBSJEZ17tqF2/eH8ML07n4/5PnQckYybMIHni3UAoGfXYmGziIhkkRtuTSB1R3lu7z8rdBSRA1JxKBJDBj9fGxJS6PO31aGjiGTcuHGMTj4bCq+ib4emodOIiGSJu/o1I0+JhQxafhGp3+kkruRMKg5FYkizUxrRsvl7LP3pXIa+MjZ0HJFDcyfl4/Es2NiOck1mUCRf4dCJRESyRFJiAl16LWPHhua8dOdroeOIpEvFoUiMGTrwdCzfr9zwQHLoKCKHNmcOozdXIXV3UdqesTd0GhGRLPXsPSeSUHAdd869mNQvvgwdR+QvVByKxJgqtSrTue1Ifl3ehkfueSd0HJGDGzuW14pEprC4+gKNUioisa1s8UJccPUctq5rzb9vHwruoSOJ/ImKQ5EY9MYb3UksspR7X65Fyl5NKC452NixfL33bJKrTuWkGvVCpxERyXKvPXAKeUos5NGfrmXX2E9CxxH5ExWHIjGoSPEiXH7up/y2vgl/u25w6Dgi6duxg+VfLGDzxoY0aLVKU1iISFwokC8PN925lt2b6nH1XRN19VByFBWHIjHqmZcvJl/Jmbzw7sls37I9dByRv5o8mZcKtgXgwi5FAocREck+D9/SikKVfuDNn2/il7eHh44j8v9UHIrEqOR8ydzZbzYpW2rQ71Ldeyg50LhxjMzTCQqv4rIzTwidRkQk2yQkGA8+Cqk7KnDxf+fDvn2hI4kAKg5FYtq//teDohW+ZNgnZ7N66ZrQcUT+ZO/HE/h5S3tKN55O8QLHhI4jIpKtrr+wCWXqfs4nC69n/ktvho4jAqg4FIlpCYkJPH7bRnxnGS7uq3kPJQdZvJiPV5Rk3+6inN5uV+g0IiJB9H+xHOwpRM+nd8BeTecj4ak4FIlxl93UmQo1PmLy192YM/Wn0HFEIsaN49XoFBZXda8eOo2ISBCdW9ekVotP+GHRFXz+YP/QcURUHIrEg1ceLAx783PR1dNDRxGJGDuWL1POIqnKFE49rnHoNCIiwQx5rREkpNBnYHHYrgHkJCwVhyJx4MxurTmuwXBm/tCNSR98EzqOxLuUFFZOnMumTQ2pf9IKEkx/ikQkfjWrU56TzpzA0qU9eaLvI6HjSJzTX2SRODHk+bqQkEK/29aHjiLxbsoU+iedAkD3cwsGDiMiEt6oN9uRp/hCbv+kL6s/mhA6jsQxFYciceL4kxtw0onvsWx+F4a8NCZ0HIlnEyfyfp5OUHgll3VsHjqNiEhwJYsW4NnnN5OyrTJtb14BO3eGjiRxSsWhSBwZ+nobLN+v3PRgvtBRJI6ljP+Un7a2o0TDaZQpVDp0HBGRHOGqHs044YyP+GlhX/576eOh40icUnEoEkcq16jEWaeP5NflbXj03++EjiPxaOdOJkzZx749xWjddkfoNCIiOconw84gb8mfuOfjviwb82noOBKHVByKxJk33zifxMLLufflGqTuSw0dR+LNV1/xSoH2YClcdn6V0GlERHKUYoXy0f+lHaTuKEvbG9fA7t2hI0mcUXEoEmeKlSzGpZ0nsHNdU/5+/Zuh40i8mTiRz1M6klj5e86op/sNRUT21/v8ppzc6QMWLbqIuy55KnQciTMqDkXi0Auv9iJviTk89/ZJ7Nymm94l+6wdO51ftjSlVovF5EnMEzqOiEiO9PF7Z5Gv9Cwe+qgPC8Z9HjqOxBEVhyJxKDlfMn+75Af2bq7J5X2Hho4j8WLTJt5YVBaALmflDRxGRCTnKpQ/mUGv7CV1V3FOv3ojqZu3hI4kcULFoUicuv+xiyhc7hveGduJ9Ss196Fkg8mTeTe5IxRcx2WdmoZOIxJTzKyjmc03s4Vmdkc6719sZjOjy9dm1ihETsm4Czofz1ndP2DV0nNp3+51SNU4AZL1ghSHZvaomf0UPUCNMLNiB2i31MxmmdkMM5uaZn1xMxtvZguij8dkW3iRGJGQmMDDt6wldUc5el36Yeg4Egf2jZ/EzB1nULjuV9QocWzoOCIxw8wSgeeBM4G6QE8zq7tfsyXAqe7eELgf6J+9KeVIjH7rPI5t8hGTpt3Erd01vYVkvVBXDscD9aMHqJ+BOw/S9nR3b+zuzdKsuwOY6O41gYnR1yJymK65rStljh3H+C/PZ8GsRaHjSIz7+uO17N1dkhanbQ4dRSTWNAcWuvtid98DvA10SdvA3b92903Rl98CFbM5oxyBhARj2mdtKFRhCk+Ovo7B9w0MHUliXJDi0N0/cfeU6MsjOUB1AQZFnw8Czs2kaCJx58X7kmB3YXpe/k3oKBLLVq3itc11gVR6n1cudBqRWFMBWJHm9crougO5DPj4QG+a2ZVmNtXMpm7YsCGTIsqRKlY4P99MLE9igV/o81h7ZozW/IeSdXLCPYf9OPAByoFPzGyamV2ZZn0Zd18DEH0sfaAv1wFO5OC69mpLjXojmDatO99MmBY6jsSqSZMY7x2h4hS6NDk5dBqRWGPprPN0G5qdTqQ4vP1AX+bu/d29mbs3K1WqVCZFlKNRv3YF3h60ltTdxTj5ikJsXLgsdCSJUVlWHJrZBDObnc7SJU2bu4AUYMgBvqaVux9PpA/9dWbW+nBz6AAncmhvPFUVcHrfvDRwEolVv370Lau3nEiVpvMonLdw6DgisWYlUCnN64rA6v0bmVlD4FWgi7v/mk3ZJJN0O7c5/7xzEjvWN6VJx2mk/rYrdCSJQVlWHLp7O3evn84yCsDM+gBnAxe7e7pnt9x9dfRxPTCCSJ96gHVmVi76PeUADbUochRatmtK06bDWDj3XEYNnhQ6jsQad977LBVIoNOZ6V3gEJGjNAWoaWbVzCwZ6AGMTtvAzCoD7wOXuPvPATJKJvjvfZ3peO47LF90Hied/Aapu/eEjiQxJtRopR2JdGfo7O7pzsBtZgXNrPDvz4EzgNnRt0cDfaLP+wCjsjaxSOwb8koLSN7O1f/WHxrJZAsWMHh3C8j/K33PahA6jUjMiY7jcD0wDpgHvOvuc8zsajO7OtrsHqAE8ML+o8BL7vLR8B7Ubf4+302/khYtVSBK5gp1z+FzQGFgfPQA9RKAmZU3szHRNmWAL83sR+B74CN3Hxt97yGgvZktANpHX4vIUajdsAZtW73P2sUdefmxkaHjSAxJHT+RaTs7kLfmZJpWaBw6jkhMcvcx7l7L3au7+3+j615y95eizy9392OiI8DvPwq85CIJCcaPX59L/RPfZ8oPl9PsxCEqECXTJIXYqLvXOMD61UCn6PPFQLoTtEb7ybfNsoAiceqtQWdR7ri1/OPJUlxxSyoJiTlhzCrJ7X4YsYhdu8vS8tQNJJj+TYmIHK2kxAR+/LorTU95jx++7kvjE95gxvcXkpAvb+hoksvpr7SI/L/SFUvT/YyP2Lq6FfffOTR0HIkFqakM/LEIABedWyJwGBGR2JGQYEz7ohtNT3mXWbN60+CE90jRIDVylFQcisifvPp6D5KKLuKhgQ1I2Zty6A+IHMyMGXy05zQoN41uzU8JnUZEJKYkJBjfT+7Oiae+zdzZvajfbAR7tqc7nIdIhqg4FJE/KVSkIFd3+4JdvzTkhiveDB1HcrktH37B0m0nUbrBVMoWKhs6johIzElIML6edCEnt3mL+XN7Uqn+FyyftTB0LMmlVByKyF88+WIv8peawasjTmP7lu2h40gu9sHIzeBJtO2QGjqKiEjMSkgwvph4ERf2HML6FW2oeWoK44aMDx1LciEVhyLyF0l5krjr8nmkbK3GpX3eCR1Hcqs9e3hjWUXIu5m+5xwXOo2ISMx7+62LeeqRcezdVYKOl53AA7cMDB1JchkVhyKSrjvvv5BiFb/g/fHnsHrpmtBxJBfyr77mm51nkFh9Eq2rnRQ6johIXLjpb2fz+Zg15CuynLuf7s1Z7V4mNWVf6FiSS6g4FJF0JSQm8NjfN+E7S9Or79hDf0BkPzPenMX2XZWo13IpeZM0vLqISHY5+bSGLJtViYq1xjJm4lUcW+8Dls9bEjqW5AIqDkXkgC67qTPlqo/h06/PZ94PP4eOI7nMaxMTgVR6dCsaOoqISNwpXeYYls05k/Yd3mDZgnM49sQEHr1zcOhYksOpOBSRg+r/QAHYW5BeV30fOorkJhs2MOLXllD5Ky5q0S50GhGRuJSQmMAnY3vz6vPjSUjazT8e6kXDJkNYv3xt6GiSQ6k4FJGDOrvHadSs9z7Tp3fnmwnTQseRXGLZ0K9ZvaMJ5Rp+RZViVULHERGJa5dd05HV80vTuPlQZs24mIqNtvHyQ++FjiU5kIpDETmkQU9VB5zeNy8NHUVyiaFDNwHQtVty4CQiIgJQslQxfviuJ4//bySpqXm4+p/n0bTZmyyfvyx0NMlBVByKyCG1bHs8TZsOY+Hccxk1eFLoOJLTufPmT9Wg5FwubXdK6DQiIpLGrXeey9KZ+Tmu0ftMn3YxVY/PT58LXmHPb7tDR5McQMWhiGTI0FdbQvJ2rrlXfzzk4DZ98xPztrSiYJ2xNC3fNHQcERHZT8UqZZj3Q3fe7D+BQsUW8cZ7V1CsyjyevV9zG8c7FYcikiE1G1Snbav3WbPoTPo/MTJ0HMnBRj2/CPck2rXfRILpz4yISE7V64oz2Lz8RG68+k327CrOjfdcSMWaHzLx/S9CR5NA9FdbRDJs8MBOWIF1/OOJUqTuSw0dR3KoAV/khUKr6df9xNBRRETkEBISE3j6xUtYu6AQp7cZxKqlbWjX7SSq132fj9/5NHQ8yWYqDkUkw8pWLsP5Z3zIllWt+N9db4eOIznQ7i27+HZ1CxJrfUT76m1DxxERkQwqWaY4kyb24cevl9P0xLdZvLADnXqcSpXaoxj15vjQ8SSbqDgUkcMy8PULSSq6mP8NqEfK3pTQcSSHmfj8PPbuK0yzZnPJnyd/6DgiInKYGp5wHFO/uZi5U9bQotUQli9tw7m921Ox5oe8+PBw9RyKcSoOReSwFCpaiCvO+4zfNjTi1msHh44jOczA9zZD8jYu7Vk7dBQRETkKdRrV4JsvL+Hnab9wSutBrFpxMtfecT6Fys3k0h6v8cvqDaEjShZQcSgih+2Zly8hX8mZvPjeyezctjN0HMkhUlNh3E91oMZYujbrGjqOiIhkgpr1q/H5Z31Yt2gvF18wgNTUJAa9cxmlqydyQos3mfzhN6EjSiZScSgihy0pTxK39ZlFypYaXN53aOg4kkNMHfcr23aVpfpxX1KmUJnQcUREJBOVrlCKwe/0Y+e6ujz1v+GUr/Q1U7/vyenntKRYxS+5tOdrLF+4InRMOUoqDkXkiNz7cE8Kl/+Gd8aeyS9rfgkdR3KAwS8sBkuhR6e8oaOIiEgWSUhM4KY7z2flz2fz/eTZtGs7kJ07SzPo7cuoUqckVY4bxX/+8Ra7du4KHVWOgIpDETkiCYkJPHTjGlJ3lOeSS0eHjiM5wLCvSkCVz+lxVq/QUUREJBuc0Lox4yf0ZdeGGrz69Ac0aDiSFctO5t+PXkSBktup1eA97v37YLb+uiV0VMkgFYcicsSuvf08Slcbz7jPurJk3rLQcbKUmXU3szlmlmpmzQ7SrqOZzTezhWZ2R3ZmDGnhAmfNpmMpXvNj6pVpEDqOiIhko4TEBC678RxmTuvJ1rV5ueOGN6lS9UsWLDiD+x7vRdHyRpXjRnHLVa+zctHK0HHlIFQcishRef4e8N1Fuejyz0JHyWqzgfOAzw/UwMwSgeeBM4G6QE8zq5s98cJ674XlAHRuvhozC5xGRERCKVS0EA8+cwlL5p7LtnWJ3H/HUI47biwrVrTkqf6XUqlWGYpV/JIO7Qcw9JWxpKbsCx1Z0lBxKCJHpdul7al63Gi+/b4bM76ZHTpOlnH3ee4+/xDNmgML3X2xu+8B3ga6ZH268AaP3AVlZtCrc+fQUUREJIcoVLQQ/3qwJ/N+vIBdvxbj2Qffp2XLoezZU5hPJvTjois7kqfYeqrXe58rLhnAV59MCR057qk4FJGjNvCxspCah17Xxm5xmEEVgLRDta2MrvsLM7vSzKaa2dQNG3L3XFHr18O8ZTVIrjmaU5pqCgsREfmr5HzJXH/HeXz9ZW92rm/ED1/P4fJeA6lY8XuWLGnNq4P7cXKHE8hTdAk16g/n8l4D+HzMd6TuSw0dPa4khQ4gIrnfaWe1oEHjocya0Y2Jo76mbZeTQkc6ImY2ASibzlt3ufuojHxFOus8vYbu3h/oD9CsWbN02+QWb7y2G/e8tKnzLcmJyaHjiIhILtC4ZT1eaVkPgJS9KYwaOoGh7y/jmzllWbz4VBbNKclrQyCh4FpKl5lB4xqr6XRaUXpeeioly5UMnD52qTgUkUwx+IVGNDplD5fdvoGlubQjpbu3O8qvWAlUSvO6IrD6KL8zR3OHZ5/dBpWmcslpjULHERGRXCgpTxLn927H+b0jr38vFoeNXsb384qzYnV9xn7SkbGfwI13pZC/5I9UKvszx9fcyhmnlqJrz1MoVuqYsDsRI1QcikimaHhiXU5qPoivv+rDO6+O48LLO4SOFMIUoKaZVQNWAT2Ai8JGylpffQXL15Qk31m3cW7Xp0PHERGRGLB/sQgw/8eFDHnjWyZ+t4/5KyqxYOHp/DyrJG+/D/1u2Uu+4j9SvvQi6lbZRKumBTn7nEbUb14n3E7kUioORSTTvDXgdKo12sgN/03gwstDp8lcZtYVeBYoBXxkZjPcvYOZlQdedfdO7p5iZtcD44BEYIC7zwkYO8s9/8x2yLuPXrXnUiB/kdBxREQkRtVuVIP/PF6D/0Rfp+5L5ZtJUxk1ai5fzXB+XlmRpctPZPHcCnz4Mdz5AFiBdRQrMY+KpVZTr+pvnNikKO071qdOk1okJGrolfSYe66+1eWwNGvWzKdOnRo6hkhMO7PDAMZ+0o9nHxzO9Xecn+XbM7Np7n7AeQdzg9x6bNq8GUqX3s3eRgOY9t8CHH9Gn9CRRHIUHZ9Est/8Hxfy0agf+HLqVuYsLcKq9dXYsak27Cn8/20s368UOmYBZYqv4thyW6hfI5GWJ5bntA5N4uJ+xoMdm3TlUEQy1Zuvd6ZMzdX885nyXHtbqs7MxbAhQ5y9e/NSu/ogjm//Teg4IiIi1G5Ug9qNanBrmnWp+1KZ9uUMJk74ie9/3MFPywux+pcKLF7WgoVzKvDJBOClSNuEQqsoVGQZJYutpXLpLdSuDI3qF+OEFsfSuEU9kvLEdvkUZO/M7FHgHGAPsAjo6+6b92tTG3gnzapjgXvc/Skzuxe4Avh9/Pd/uvuYrM4tIodWslxJLuz4GkOHX8Z/bh/CvY9dHDqSZAF3eOa5bVB2ITfWOgY08b2IiORQCYkJnHBqY044tfFf3lu7fB2ffTKDb6esZdbCVJavL8z6TWVZtqIZi+dWZHLamf8S9pBUeBGFCq+iZNENVCi5jarlUql9bAHq1S9Ns5Z1KF+1XPbtWBYI0q3UzM4AJkXvz3kYwN1vP0j7RCKDO5zo7suixeF2d3/scLarrhEi2WPntp0UrbSapDw72ba6bpaeZVO3rTCmTYNmzSBPh2tZ/3AfijU6MXQkkRxHxyeR3G3rxq18NWkGU6as5KdFv7F4dR7WbCzKxq1l2bm1Eqk7yv/1Q3k3k7fgagoWWscxhTdSpth2KpRKoVrFPNQ8tih16legXpMaFCtZLNv353c5rlupu3+S5uW3QLdDfKQtsMjdl2VdKhHJLAUKF+Dqbl/y3GuXctNVA3l+QN/QkSSTvfDSXkjay/mlPqJYoxdCxxEREcl0RYoX4cxurTnzAJXK9i3bmf7VXGbOXMG8n7eyZLWzckM+ftlSjK3bS7F0Uw0WzS0H/teSy/JtJE+BdRQosIEiBTdRvPB2ShXdTflSUKlcMlWrFKF6jdLUqleFspXKZNttOjmh02w//tx9ND09gKH7rbvezHoDU4G/ufumrAgnIkfmyRd78eqoH+n/fmsefnI7hYoWCh1JMsmOHfDWW6lQ712urdspdBwREZEgChUtROtOzWndqfkB2+zZtYeffpzLzOlLWbhoE0tX7WblemP95nxs3FqErTtKsGZdXZYvLQ27i6b/JYm7SMj/C3nz/0r+fJsoVGALRQvupHjh3ZQsto8yxRPo0L4KnXu2Oep9yrLi0MwmAGXTeesudx8VbXMXkAIMOcj3JAOdgTvTrH4RuB/w6OPjRIrM9D5/JXAlQOXKlQ97P0TkyCTlSeKOvrO599GLuaJf5B5EiQ3vvgu7dualUvVXOfmiAx6+RURE4l5yvmQanliXhifWPWTbjes2Mn/WEn7+aQ1Llm5h5drdrPnV+WVLHjZtK8CWHUXY+Vsx1qyrwIpdJfBdxf//s7MWvk7nnkefN8uKQ3dvd7D3zawPcDbQ1g9+4+OZwHR3X5fmu///uZm9Anx4kBz9gf4Q6TefsfQikhnufrAnTwz5mnfHdeLplespXbF06EiSCZ55YSeUXMb1+VZjVaqEjiMiIhITipcpTssyxWl50CrqD7t27mLJ/OUsWbCGsuWbZkqGIGPMm1lH4Hags7vvPETznuzXpdTM0g4D1BWYnbkJRSQzJCQm8NDN60jdUY7e/Q54DkdykTlzYMbUAiQ0fpU+J10dOo6IiEjcylcgH3Wa1KLTBady/MkNMuU7Q01A9hxQGBhvZjPM7CUAMytvZv8/JYWZFQDaA+/v9/lHzGyWmc0ETgduyabcInKYrrmtK6WrjeeTz7uyZJ7GlMrtXu6/DxL2cGbBNylzYbq9+UVERCSXCjVaaY0DrF8NdErzeidQIp12l2RdOhHJbM/+K5ULLzuGiy8fzddf9QkdR47Q7t0wcFAK1BnJdQWrQsmSoSOJiIhIJgp15VBE4sgF/TpQpfYovpnSjZnfzQ0dR47QwIGwfUteStZ+lTM63Rg6joiIiGQyFYciki1ee7gU7Eum17U/ho4iR2DDBrjjzn1Q9VOu2TqZxC7nho4kIiIimUzFoYhki7ZdTqJOg/eZNeN8Ph/zfeg4cpjuuAO2bnPydLiOy8t1gkKat1JERCTWqDgUkWzzxrPHQcI++t62KnQUOQxffw0DBoCf+AR3z5tH5e6Xh44kIvsxs45mNt/MFprZHem8b2b2TPT9mWZ2fIicIpKzqTgUkWzT7JRGnNBsGIvndmHEmxNCx5EMSEmBq6/ZR2Kx1dRu9D/+sb0hnHlm6FgikoaZJQLPE5kbui7Q08z2n3H7TKBmdLkSeDFbQ4pIrqDiUESy1ZD+J0HebVx7b0roKJIBL74Is2Ymsq/9jfQfs4W8z70ESUEGuhaRA2sOLHT3xe6+B3gb6LJfmy7AGx7xLVBsv3mjRURUHIpI9qrZoDptW73P2sUd6f/EyNBx5CDWroV/3rUPqo+n3+73ad3uMmjZMnQsEfmrCsCKNK9XRtcdbhsRiXMqDkUk2w0e2AkrsI7bnyxB6r7U0HHkAP7+91R2/LaPY9rdzKNTisJDD4WOJCLps3TW+RG0iTQ0u9LMpprZ1A0bNhx1OBHJPVQciki2K1u5DOe1/5DNK0/h4XveDR1H0vHZZzBkSALe8hGe/Xouxe99WJPei+RcK4FKaV5XBFYfQRsA3L2/uzdz92alSpXK1KAikrOpOBSRIAYM6E5ikaU88GotXT3MYfbuhauu2YsVW8bpxz7KRflPgMs1QqlIDjYFqGlm1cwsGegBjN6vzWigd3TU0hbAFndfk91BRSRnU3EoIkEUKV6ES8+ZxM71x3P7DYNDx5GolBS4/nqYPy8PSR1upv+H27AXXoQE/bkQyancPQW4HhgHzAPedfc5Zna1mV0dbTYGWAwsBF4Brg0SVkRyNA05JyLBvPBqL978eB7PvNOc+x/bRb4C+UJHimtbt0KPHvDxx0Crh/n36lHU6HEtNG0aOpqIHIK7jyFSAKZd91Ka5w5cl925RCR30algEQkmOV8yN144hT0bj+P6K98KHSeurVgBrVqlMnbcPjj7Sk4//n/ctqAkPPBA6GgiIiKSTVQcikhQDz/biwJlpvH66LZs3bg1dJy4NG0aNGuewtyFO/GLO3J73Wl88vxWkh9+DIoVCx1PREREsomKQxEJKiExgbuvWMC+bVXo2/e90HHizqhR0OqUFDbsWk3BK05n9PI5PPTIdJKuuwEuuSR0PBEREclGKg5FJLh/3HsBxSp+zogJZ7N2+brQcWJeaip8/TVcfU0q53ZNZXfxadS/9ExmDJ7KOSvyw6RJ8MwzYOlNiyYiIiKxSsWhiASXkJjAo7duwneW4ZJ+Yw79ATls7vDDD3D9LTsoXXEHrVrBy6/ugcYD6HfSOXz/3FyO7X0TzJwJp58eOq6IiIgEoNFKRSRHuPyWLtz93FgmfNmVRXOWUL1etdCRcqWdO2H5in38tHgb85fsYPGyPSxfkcr3XxRh48pSkJAM1cdRqMtwOh7zKZd8s4zOP9aCT9+Hk08OHV9EREQCUnEoIjnGc/ck0u3SIlx0xWi++zq+i8N3B4zjsruKA5GrfmCRxQ3H8NQ8pO5LJnVfMr4vL6mpefCUvPjewkAiUCy6APk2YmV/oMZJb3NB6vt0Xb6R478rQ0KdutCzL/zjH5A/f4C9FBERkZxExaGI5Bjn92nPcU+8Q6F8+0JHCS5PnkTyJm8D8+gax/BIjYiTmLiHxIS9JCbuISlhL0mJe0hM2EPRfFsoV+BXKubfxLEFt1KjwHYq5U3k2CJVKVKnJdTpB8cdB8ccE3DvREREJCdScSgiOcq8Hy8MHSFH6HpJO7pqsFARERHJRhqQRkRERERERFQcioiIiIiIiIpDERERERERQcWhiIiIiIiIoOJQREREREREUHEoIiIiIiIiqDgUERERERERVByKiIiIiIgIKg5FREREREQEFYciIiIiIiKCikMRERERERFBxaGIiIiIiIig4lBEREREREQAc/fQGbKNmW0Alu23uiTwS4A42S0e9jMe9hHiYz8PZx+ruHuprAyT1cxsGzA/0OZD/nuKx23H4z7H87Zru3vhQNvOFPrtFPP7qX2MHZny2ykp8/LkfOn9RzCzqe7eLESe7BQP+xkP+wjxsZ/xsI/7mR9qf0P+t47HbcfjPsf7tkNsNzPpt1Ns76f2MXZk1n6qW6mIiIiIiIioOBQREREREREVhwD9QwfIJvGwn/GwjxAf+xkP+5hWyP3VtuNju9p2/G07K8Xqfu0vHvZT+xg7MmU/42pAGhEREREREUmfrhyKiIiIiIhIfBeHZtbRzOab2UIzuyN0nqxgZgPMbL2ZzQ6dJauYWSUz+9TM5pnZHDO7KXSmzGZm+czsezP7MbqP94XOlFXMLNHMfjCzD0NnyWyHOuZYxDPR92ea2fHZuO3TzGyLmc2ILvdk0nYPegzK4n0+1Lazap8PeUzKqv3O4Lazar8PeZzKwv3OyLazZL+j333A41ZW/hsPQb+dYoN+O8WWTP3t5O5xuQCJwCLgWCAZ+BGoGzpXFuxna+B4YHboLFm4j+WA46PPCwM/x9r/loABhaLP8wDfAS1C58qifb0VeAv4MHSWTN6vQx5zgE7Ax9H/vVsA32Xjtk/Liv/mhzoGZdU+Z3DbWbXPhzwmZeH/1hnZdlbt9yGPU1m43xnZdpbsd/S7D3jcysp/49m9ZORYEgvLoY4dsbBk5FiR25eMHBdiZTnYMehwl3i+ctgcWOjui919D/A20CVwpkzn7p8DG0PnyEruvsbdp0efbwPmARXCpspcHrE9+jJPdIm5G4bNrCJwFvBq6CxZICPHnC7AG9H/vb8FiplZuWzadpbIwDEoq/Y52PEvg8ekLNnvkMfDDB6nsmq/gx0jM3DcyrJ/4wHot1OM0G+n2JHZv53iuTisAKxI83olMfZ/inhkZlWBJkTODsWUaJeBGcB6YLy7x9w+Ak8B/wBSA+fIChk55mTVcSmj39sy2v3mYzOrlwnbzYjQx+Is3eeDHJOyfL8PcTzMkv3OwHEqy/Y7g8fIrNjvpzj4cSv0v/HMFEv7IlH67ZTrPUUm/naK5+LQ0lkXc2cT4omZFQKGAze7+9bQeTKbu+9z98ZARaC5mdUPHClTmdnZwHp3nxY6SxbJyDEnq45LGfne6UAVd28EPAuMzITtZkTIY3GW7vMhjklZut+H2HaW7XcGjlNZtt8Z2Ham73cGj1ux9HsjlvZF0G+n3C4rfjvFc3G4EqiU5nVFYHWgLHKUzCwPkYPbEHd/P3SerOTum4HJQMewSTJdK6CzmS0l0lWpjZkNDhspU2XkmJNVx6VDfq+7b/29+427jwHymFnJTNj2UWfLKlm5zxk4JmXZfh9q29nxv/VBjlNZ/r/3gbadRfudkeNWLP3eiKV9iXv67RQTMv23UzwXh1OAmmZWzcySgR7A6MCZ5AiYmQGvAfPc/YnQebKCmZUys2LR5/mBdsBPQUNlMne/090runtVIv9/nOTuvQLHykwZOeaMBnpHRzdsAWxx9zXZsW0zKxv9/xJm1pzI34dfM2Hbh5JV+3xIWbXPGTwmZcl+Z2TbWbjfGTlOZdV+H3LbWbHfGTxuBfs3ngX02ylG6LdTbMiK305JmZIsF3L3FDO7HhhHZPStAe4+J3CsTGdmQ4mM0FbSzFYC/3b318KmynStgEuAWdF+5QD/jJ4ZjhXlgEFmlkjkB8277h5zUz3EsgMdc8zs6uj7LwFjiIxsuBDYCfTNxm13A64xsxTgN6CHux91d7H0jkFEBgXI0n3O4LazZJ85wDEJqJxm21m13xnZdlbtd7rHqez4N57BbWfVfv9FNu1zttNvp5ii306SLsui46KIiIiIiIjkIvHcrVRERERERESiVByKiIiIiIiIikMRERERERFRcSgiIiIiIiKoOBQRERERERFUHEo6zKyYmV0bOkdGmNnNZlYgi777XjP7+37rnjezGWY218x+iz6fYWbdou+3NLNXzOw0M/vLcMlmdqmZbTCzVw+wzUfNbO3+2xURHZvSfLeOTSI5jI5P///dOj7lcnE7z6EcVDHgWuCFwDl+n6TV3D31AE1uBgYTmTsqo9+Z5O4pR5LH3a+LfkdV4EN3b7xfk47A2EN8zTvufv0Bvv82M9txJNlE4kAxdGxKl45NIsEVQ8endOn4lLvoyqGk5yGgevSszqMAZnabmU0xs5lmdl90XVUz+8nMXjWz2WY2xMzamdlXZrbAzJpH291rZm+a2aTo+it+39BBvneemb0ATAcqmdmLZjbVzOakaXcjUB741Mw+ja7bnua7u5nZ69Hnr5vZE9F2D5tZdTMba2bTzOwLMzvuYP9BzOwKM/vYzPIf4r9dW2BChv9Li8jh0LFpPzo2ieQYOj7tR8en3ElXDiU9dwD1fz+zY2ZnADWB5oD9X3t3ExpXFYZx/P9UxYLRjaJoRRREqUis+ImIWNFAXWhd1IUfYF1UEUEQcaEuCgpuxIUgoiBuC10phqKQnaLoIta2K8VglYjgV2slSGlfF/fWTJK5zYy1OGP/P5jFTc5959w7Mw+cc8/cAd5LchuwH7gc2AJsAz4HHgBuBe4BngM2tzUngZuBs4DZJNPA1cepeyWwtaqeaPvwfFX9kuQ0YCbJZFW9luRpYGNV/TTAcV0B3FlVR5LMAI9X1VdJbqKZ6buj305JngSmgM1V9WdX8STnAYer6kCSAbojaUhmUw+zSRop5lMP82l8OTjUIKbax2y7PUETTPuBuaraA5BkHzBTVZVkD3BpT413q2oBWGhnoG6kCcKuut9W1ac9+9+fZBvNe/ZC4CrgyyGPY2cbbhPALcDOniA6s2Ofh4HvacLt8Cr1p4APh+yTpH/ObDKbpFFlPplPY8nBoQYR4OWqenPJH5u1472zQUd7to+y9P1Vy2rWKnX/6Nm+DHgGuKGqfm2XO6zt6Gvv8yxvc6zmGuC3Pmve+9kLbAAuBuZWabsJeHWAmn9L8g5wLTBfVXcPs68kswmzSRpV5pP5NJb8zqH6+R04u2f7A+DRdtaIJOuSnD9kzXuTrE1yLnA7zTKKQeueQxNOB5JcQBMkXX39Mcn6JGuA+/p1pKoOAnNJtrTPmyTXdPR7FniMZtnGRV0Hl2YabRL4oqtNR1+2VtUGw00aiNm0yGySRov5tMh8GmNeOdQKVfVzmi9G7wV2tXeBWg980i4lOAQ8BBwZouxnwDRwCfBiVc0D84PUrardSWaBfcA3wMc9/34L2JXkh6raSLPm/33gO5qZq4mO/jwIvJHkBeAMYAewu+N8fJTm9sjTSe7qWKN/HTBbVcdm305n6cygpBNkNq04H2aTNCLMpxXnw3waU1l8PaSTI8l24FBVvfJf9+VkaYPy66ra0W4/BayrqmeXtXsEuL7rdsxtm+38z8+XNApOhc+a2SSNp1Ph82Y+jSaXlUr/gqp6qSfc3qa589jrfZouAJtynB9ypZkB9Pd6JJ0ws0nSqDKfRpNXDiVJkiRJXjmUJEmSJDk4lCRJkiTh4FCSJEmShINDSZIkSRIODiVJkiRJODiUJEmSJAF/AUyu4uzc8TdpAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Parameters\n",
    "kTs =np.linspace(0.1,4,40)\n",
    "M = 10\n",
    "J = -1\n",
    "mu = 0.1\n",
    "\n",
    "# Prepare plots\n",
    "fig = plt.figure(figsize=(15,6))\n",
    "ax1 = plt.subplot(1,3,1)\n",
    "ax2 = plt.subplot(1,3,2)\n",
    "ax3 = plt.subplot(1,3,3)\n",
    "\n",
    "# Analytical solution (red)\n",
    "an = analytical_solution(J,mu,kTs)\n",
    "ax1.plot(np.abs(an.kT/J),an.g/np.abs(J),'r')\n",
    "ax2.plot(np.abs(an.kT/J),an.u/np.abs(J),'r')\n",
    "ax3.plot(np.abs(an.kT/J),an.m,'r')\n",
    "\n",
    "# Brute-Force (green)\n",
    "start = time.time()\n",
    "bf = brute_force(M,J,mu,kTs)\n",
    "end = time.time()    \n",
    "print(\"Brute force solution took {} sec\".format(end-start))\n",
    "ax1.plot(np.abs(bf.kT/J),bf.g/np.abs(J),'g-')\n",
    "ax2.plot(np.abs(bf.kT/J),bf.u/np.abs(J),'g-')\n",
    "ax3.plot(np.abs(bf.kT/J),bf.m,'g-')\n",
    "bfg_arr = bf.g.to_numpy()\n",
    "ang_arr = an.g.to_numpy()\n",
    "\n",
    "bfm_arr = bf.m.to_numpy()\n",
    "anm_arr = an.m.to_numpy()\n",
    "\n",
    "print(\"Avg. error Brute Force free energy       : {}\".format(np.average(np.abs(bfg_arr - ang_arr))))\n",
    "print(\"Avg. error Brute Force magnetisation     : {}\".format(np.average(np.abs(bfm_arr - anm_arr))))\n",
    "\n",
    "    \n",
    "# Transfer Matrix (blue)\n",
    "start = time.time()\n",
    "tm = transfer_matrix(M,J,mu,kTs)\n",
    "tmg_arr = tm.g.to_numpy()\n",
    "tmm_arr = tm.m.to_numpy()\n",
    "end = time.time()    \n",
    "print(\"Transfer matrix solution took {} sec\".format(end-start))\n",
    "ax1.plot(np.abs(tm.kT/J),tm.g/np.abs(J),'b-')\n",
    "ax2.plot(np.abs(tm.kT/J),tm.u/np.abs(J),'b-')\n",
    "ax3.plot(np.abs(tm.kT/J),tm.m,'b-')\n",
    "print(\"Avg. error Transfer Matrix free energy   : {}\".format(np.average(np.abs(tmg_arr-ang_arr))))\n",
    "print(\"Avg. error Transfer Matrix magnetisation : {}\".format(np.average(np.abs(tmm_arr-anm_arr))))\n",
    "\n",
    "ax1.set_xlabel('temperature kT/J [-]')\n",
    "ax1.set_ylabel('free energy g/J [-]')\n",
    "ax2.set_xlabel('temperature kT/J [-]')\n",
    "ax2.set_ylabel('internal energy u/J [-]')\n",
    "ax3.set_xlabel('temperature kT/J [-]')\n",
    "ax3.set_ylabel('magnetisation m [-]')\n",
    "\n",
    "ax2.set_xlim([0,4])\n",
    "ax3.set_ylim([-0.1,1.1])\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}