MC_v3_r2_CL.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "a681aef5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import random as rd\n",
    "import numpy as np\n",
    "from time import time\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "bafbce71",
   "metadata": {},
   "outputs": [],
   "source": [
    "class site:\n",
    "    def __init__(self,lattice,x,y,z):\n",
    "        self.x = x\n",
    "        self.y = y\n",
    "        self.z = z\n",
    "        self.lattice = lattice\n",
    "\n",
    "    def setspin(self,s):\n",
    "        self.spin = s\n",
    "    \n",
    "    def setrandomspin(self): # set random spin for this site\n",
    "        self.spin = rd.choice([-1,1]) # set spin = -1 or 1\n",
    "     \n",
    "    def getspin(self):\n",
    "        return self.spin\n",
    "\n",
    "    def flipspin(self):\n",
    "        self.spin *= -1\n",
    "\n",
    "    def getposition(self):\n",
    "        return np.array([self.x,self.y,self.z])\n",
    "\n",
    "    def getzeemanenergy(self,B):\n",
    "        return -g*muB*self.getspin()*B\n",
    "\n",
    "    def getexchangeenergy(self, J):\n",
    "        '''\n",
    "        This is a very similar method to the ones you wrote below, but just calculates the exchange energy of all exchange interactions involving the site, not the change.\n",
    "        Accordingly we have not multiplied by 2. We can use this to calculate the changes in energy at a later point, in the MC algorithm itself. We can also use this to \n",
    "        calculate the total lattice energy \n",
    "        '''\n",
    "        E = 0\n",
    "        NN = self.getnearestneighbourscubic()\n",
    "        for nn in NN:\n",
    "            E += -J*self.spin*nn.spin\n",
    "        return E\n",
    "        \n",
    "\n",
    "    def getdeltaexchangeenergy(self,J):\n",
    "        # This function is fine, but will only work for Ising-like spins (i.e. just 2 possible states). The reason is because you need no prior information to know what the state \n",
    "        # before was, and therefore working out the change in energy is easy (every positive contribution to energy turns into a negative and vice versa). This will not be the \n",
    "        # case for orbitals, which will be able to take 3 states, so you will need to explicitly use the pre-flipped state in calculating your answer. It will also not be the case \n",
    "        # more spin models, and certainly not for heisenberg spins (freely rotating in 3 dimensions).    \n",
    "        E = 0\n",
    "        NN = self.getnearestneighbourscubic()\n",
    "        for nn in NN:\n",
    "            E += -J*self.spin*nn.spin\n",
    "        E *= 2 # You forgot to multiply by 2 here\n",
    "        return E\n",
    "    \n",
    "    def getdeltaenergy(self,J,g,muB,B):\n",
    "        # \n",
    "        E = 0\n",
    "        NN = self.getnearestneighbourscubic()\n",
    "        for nn in NN:\n",
    "            E += -J*self.spin*nn.spin\n",
    "        E += g*muB*B*self.spin\n",
    "        E *= 2\n",
    "        return E\n",
    "    \n",
    "    def getsiteenergy(self,J,B):\n",
    "        E = 0 \n",
    "        E += self.getexchangeenergy(J)\n",
    "        E += g*muB*B*self.spin\n",
    "        return E\n",
    "    \n",
    "    \n",
    "    def getlistposition(self,x,y,z):\n",
    "        return x*self.lattice.N2*self.lattice.N3 + y*self.lattice.N3 + z\n",
    "\n",
    "    def getnearestneighbourscubic(self):\n",
    "        #i-1 i+1, j-1 j+1, k-1 k+1 with periodic boundary conditions\n",
    "        NN1 = self.lattice.sites[self.getlistposition((self.x-1) % self.lattice.N1, self.y, self.z)]\n",
    "        NN2 = self.lattice.sites[self.getlistposition((self.x+1) % self.lattice.N1, self.y, self.z)]\n",
    "        NN3 = self.lattice.sites[self.getlistposition(self.x, (self.y-1) % self.lattice.N1, self.z)]\n",
    "        NN4 = self.lattice.sites[self.getlistposition(self.x, (self.y+1) % self.lattice.N1, self.z)]\n",
    "        NN5 = self.lattice.sites[self.getlistposition(self.x, self.y, (self.z-1) % self.lattice.N1)]\n",
    "        NN6 = self.lattice.sites[self.getlistposition(self.x, self.y, (self.z+1) % self.lattice.N1)]\n",
    "        return [NN1,NN2,NN3,NN4,NN5,NN6]\n",
    "    \n",
    "    def getnumberofalignedspins(self):\n",
    "        NN = self.getnearestneighbourscubic()\n",
    "        a = 0\n",
    "        for nn in NN:\n",
    "            a += self.spin*nn.spin\n",
    "        a = (a+6)/2\n",
    "        return a \n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "b82f2b72",
   "metadata": {},
   "outputs": [],
   "source": [
    "class lattice:\n",
    "    def __init__(self,N1,N2,N3):\n",
    "        self.N1 = N1 # Dimension of the lattice\n",
    "        self.N2 = N2\n",
    "        self.N3 = N3\n",
    "        self.sites = [] # Creating a list to take in the information of sites\n",
    "\n",
    "    def getnumberofsites(self):\n",
    "        return self.N1*self.N2*self.N3\n",
    "    \n",
    "    def populatesites(self,lat):\n",
    "        for i in range(self.N1):\n",
    "            for j in range(self.N2):\n",
    "                for k in range(self.N3):\n",
    "                    self.sites.append(site(lat,i,j,k)) # Assign each site a coordianate\n",
    "\n",
    "    def initialisespins(self): # Set spin for each site\n",
    "        for i in range(self.N1*self.N2*self.N3):\n",
    "            self.sites[i].setrandomspin()\n",
    "\n",
    "    def getlistposition(self,x,y,z):\n",
    "        return x*self.N2*self.N3 + y*self.N3 + z\n",
    "            \n",
    "    def gettotalenergy(self,J, B):\n",
    "        '''\n",
    "        This now calculates the zeeman energy too, and has been reformulated in a simpler way. \n",
    "        '''\n",
    "        E = 0 #meV\n",
    "        for i in range(self.N1):\n",
    "            for j in range(self.N2):\n",
    "                for k in range(self.N3):\n",
    "                    # Why perform a new energy evaluation calculation when you can call your site methods to calculate it, as below\n",
    "                    site0 = self.sites[self.getlistposition(i, j, k)]\n",
    "                    E += site0.getexchangeenergy(J)/2 + site0.getzeemanenergy(B) # divide exchange energies by 2 to account for double counting of spin pairs\n",
    "                    # There is no need to do 3 separate sets of nested for loops, you can perform the same operation within 1, as I have done now\n",
    "                    # Note that the .getexchangeenergy() method calculates exchange energies with neighbours along all 3 axes in 1 go.\n",
    "        return E\n",
    "    \n",
    "    def pickrandomsite(self):\n",
    "        site_f = rd.randrange(self.N1*self.N2*self.N3)\n",
    "        return site_f\n",
    "    \n",
    "    def filpspin(self,site_f):\n",
    "        self.sites[site_f].flipspin()\n",
    "        \n",
    "    def getmagnetization(self):\n",
    "        M = 0\n",
    "        for i in range(self.N1):\n",
    "            for j in range(self.N2):\n",
    "                for k in range(self.N3):\n",
    "                    M += -self.sites[self.getlistposition(i, j, k)].getspin()\n",
    "        return M/(self.N1*self.N2*self.N3)\n",
    "\n",
    "    def tensorrepresentation(self):\n",
    "        '''\n",
    "        A method returning a tensor representation of the lattice i.e. a 3D array\n",
    "        '''\n",
    "        tensor = np.zeros((self.N1, self.N2, self.N3)) # initialize empty array\n",
    "        for i in range(self.N1):\n",
    "            for j in range(self.N2):\n",
    "                for k in range(self.N3):\n",
    "                    # Why perform a new energy evaluation calculation when you can call your site methods to calculate it, as below\n",
    "                    tensor[i,j,k] = self.sites[self.getlistposition(i, j, k)].getspin()\n",
    "        return tensor\n",
    "\n",
    "    def latticeslice(self, axis = 0, index = 0):\n",
    "        '''\n",
    "        Slices the tensor representation of the lattice at position 'index' along a particular axis\n",
    "        '''\n",
    "        latt = self.tensorrepresentation()\n",
    "        lattslice = latt.take(indices = index, axis = axis)\n",
    "        return lattslice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "01f054a8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1789.5815999999882\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-14-df1f52430166>:56: RuntimeWarning: divide by zero encountered in log\n",
      "  CL_1 = -1/np.log(CF)\n",
      "<ipython-input-14-df1f52430166>:56: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  CL_1 = -1/np.log(CF)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Runtime:  8.332524\n",
      "-23.65982328888739\n"
     ]
    }
   ],
   "source": [
    "starttime = time()\n",
    "\n",
    "MCS = 1000 # MC steps per temperature\n",
    "J = 10 # meV\n",
    "T = 500 # Kelvin\n",
    "T_p = 0.97\n",
    "kB = 8.617e-2 # meV/K\n",
    "g = 1\n",
    "muB = 5.788e-2 # meV/T\n",
    "B = 2 #T\n",
    "\n",
    "MCA_plot_1 = []\n",
    "T_plot_1 = []\n",
    "E_plot_1 = []\n",
    "Mag_plot_1 = []\n",
    "CL_plot_1 = []\n",
    "\n",
    "# set up the sites\n",
    "lat1 = lattice(30,30,30) # set up the size of lattice\n",
    "lat1.populatesites(lat1)     # Assign each site a coordinate\n",
    "lat1.initialisespins()   # Assign each site a random spin\n",
    "slc0 = lat1.latticeslice() # print lattice slice to check random initial config\n",
    "E_i = lat1.gettotalenergy(J,B) # get initial energy per site of the system\n",
    "print(E_i)\n",
    "NoS = lat1.getnumberofsites()\n",
    "E_1 = E_i/NoS\n",
    "\n",
    "while (T>0.05):\n",
    "    \n",
    "    MCA = 0    # Initialise MC accepts\n",
    "    Mag_1 = 0  # Initialise Magnetization\n",
    "    AS_1 = 0 # Initialise number of aligned spins\n",
    "    #CF_w = []\n",
    "    #CF_avg_1 = []\n",
    "    #CF = 0\n",
    "    cs = lat1.pickrandomsite() # pick a random site to calculate correlation length\n",
    "    \n",
    "    for i in range(MCS):\n",
    "        fs = lat1.pickrandomsite() # pick a random site to flip its spin\n",
    "        E_b = lat1.sites[fs].getsiteenergy(J,B) # Preflip energy of the site\n",
    "        lat1.filpspin(fs) # flip spin\n",
    "        E_a = lat1.sites[fs].getsiteenergy(J,B) # Postflip energy of the site\n",
    "        dE = E_a - E_b     \n",
    "        \n",
    "        if dE <= 0 or np.exp(-dE/(kB*T)) > rd.random():\n",
    "            E_1 += dE/NoS\n",
    "            MCA += 1 # MC Accepts\n",
    "            \n",
    "            # correlation length\n",
    "            AS_1 += lat1.sites[cs].getnumberofalignedspins() # calculate the number of aligned sites\n",
    "            \n",
    "        else:\n",
    "            lat1.filpspin(fs) # flip the spin back\n",
    "    \n",
    "    CF = AS_1/(6*MCA) # calculate the correlation function\n",
    "    CL_1 = -1/np.log(CF)\n",
    "    Mag_1 = lat1.getmagnetization() # magnetization of the system\n",
    "    \n",
    "    MCA_plot_1.append(MCA)\n",
    "    T_plot_1.append(T)\n",
    "    E_plot_1.append(E_1)\n",
    "    Mag_plot_1.append(Mag_1)\n",
    "    CL_plot_1.append(CL_1)\n",
    "    \n",
    "    T = T_p * T\n",
    "    \n",
    "print (\"Runtime:  %.6f\"  %(time()-starttime))\n",
    "slc1 = lat1.latticeslice()\n",
    "print(E_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "a6878592",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAY2UlEQVR4nO3de5BsZ1nv8d9vZu8t6SSWZGZIpYTpRk1ZIsUJZIh4olY8CRg4FIgXBMdtUHR0b/SQg7fgnKPwx5R3xFLr6HjIMcW0iJdEIiCeGBMoC4swG0IuFTARZrbIPtlbvACOBZr9nD969aanp3t69eXty1rfT9Wq6X57rXc97+qeZ9a871pvOyIEACiPuUkHAAAYLxI/AJQMiR8ASobEDwAlQ+IHgJI5MukA8lhcXIxarTbpMABgppw6deofImKpvXwmEn+tVtP29vakwwCAmWJ7t1M5XT0AUDIkfgAoGRI/AJQMiR8ASobEDwAlU9jEX6/XVavVNDc3p1qtpnq9PumQAGAqzMTlnP2q1+taW1vT3t6eJGl3d1dra2uSpNXV1UmGBgATV8gz/vX19QtJv2lvb0/r6+sTiggApkchE//p06f7KgeAMilk4l9eXu6rHADKpJCJf2NjQ5VKZV9ZpVLRxsbGhCICgOlRyMS/urqqzc1NVatV2Va1WtXm5iYDuwAgybPwnbsrKyvBJG0A0B/bpyJipb28kGf8AIDuSPwAUDIkfgAoGRI/AJQMiR8ASobEDwAlQ+IHgJIh8QNAyZD4AaBkSPwAUDIkfgAomWSJ3/bTbN9j+xHbD9t+bVZ+me27bD+a/XxyqhgAAAelPOP/D0k/FhFfI+l5kl5j+xmSbpF0d0RcKenu7DkAYEySJf6IOBMRH8oef1bSI5K+XNJLJd2WrXabpG9NFQMA4KCx9PHbrkl6tqQPSLo8Is5IjT8Okp4yjhgAAA3JE7/tSyT9saSbI+IzfWy3Znvb9va5c+fSBQgAJZM08ds+qkbSr0fE7Vnx47avyF6/QtLZTttGxGZErETEytLSUsowAaBUUl7VY0lvkfRIRLyp5aU7Jd2UPb5J0jtSxQAAOOhIwrqvlXRc0oO278/KflrSz0v6A9uvlnRa0ncmjAEA0CZZ4o+Iv5LkLi9fn2q/AIDDcecuAJQMiR8ASobEDwAlQ+IHgJIh8QNAyZD4AaBkSPw4VL1eV61W09zcnGq1mur1+qRDAjCklDdwYcbV63Wtra1pb29PkrS7u6u1tTVJ0urq6iRDAzAEzvjR1fr6+oWk37S3t6f19fUJRQRgFEj86Or06dN9lQOYDSR+dLW8vNxXOYDZQOJHVxsbG6pUKvvKKpWKNjY2JhQRgFEg8aOr1dVVbW5uqlqtyraq1ao2NzcZ2AVmnCNi0jH0tLKyEtvb25MOAwBmiu1TEbHSXs4ZPwCUDIkfAEqGxA8AJUPiB4CSIfEDQMmQ+AGgZEj8AFAyJH4AKBkSPwCUDIkfAEqGxA8AJUPiB4CSIfEDQMmQ+AGgZEj8AFAyJH4AKBkSPwCUDIkfAEqGxA8AJUPiB4CSIfEDQMmQ+AGgZJIlftu32j5r+6GWsjfY/nvb92fLi1LtHwDQWcoz/t+VdGOH8l+NiKuy5d0J9w8A6CBZ4o+I90n6x1T1AwAGM4k+/h+x/UDWFfTkbivZXrO9bXv73Llz44wPAApt3In/f0n6SklXSToj6Ve6rRgRmxGxEhErS0tLYwoPAIpvrIk/Ih6PiCci4ryk35F0zTj3DwAYc+K3fUXL05dJeqjbugCANI6kqtj22yRdJ2nR9icl/ayk62xfJSkk7Uj6oVT7BwB0lizxR8QrOxS/JdX+AAD5cOcuAJQMiR8ASobEDwAlQ+IHgJLJPbhre17S5a3bRMTpFEEBANLJlfht/6gal2M+Lul8VhySnpUoLgBAInnP+F8r6asj4tMpgwEApJe3j//vJP1LykAAAOOR94z/45Lutf0uSZ9vFkbEm5JEBQBIJm/iP50tx7IFADCjciX+iHijJNm+tPE0Ppc0KgBAMrn6+G0/0/aH1ZhN82Hbp2x/bdrQAAAp5B3c3ZT0uoioRkRV0o+pMZ8+AGDG5E38F0fEPc0nEXGvpIuTRAQASCr3VT22/6ekt2bPv0fSJ9KEBABIKe8Z//dLWpJ0u6Q7ssfflyooAEA6ea/q+SdJ/y1xLACAMTg08dt+c0TcbPtP1ZibZ5+IeEmyyAAASfQ642/26f9y6kAAAONxaOKPiFPZw6si4tdaX7P9WknvTRUYACCNvIO7N3Uoe9UI4wAAjEmvPv5XSvpuSU+3fWfLS5dKYopmAJhBvfr43y/pjKRFSb/SUv5ZSQ+kCgoAkE6vPv5dSbuSvn484QAAUss7SdvzbH/Q9udsf8H2E7Y/kzo4AMDo5R3c/Q1Jr5T0qKSLJP2ApF9PFRQAIJ28c/UoIh6zPR8RT0j6P7bfnzAuAEAieRP/nu1jku63/YtqDPgyOycAzKC8XT3HJc1L+hFJ/yrpaZK+PVVQAIB08k7Stps9/DdJb0wXDgAgtV43cD2oDpOzNUXEs0YeEQAgqV5n/C8eSxQAgLHJcwOXJMl2VdKVEfEXti/qtS0AYDrlvYHrByX9kaTfzoqeKulPEsUEAEgo71U9r5F0raTPSFJEPCrpKamCAgCkkzfxfz4ivtB8YvuIDhn0BQBMr7yJ/722f1rSRbafL+kPJf3pYRvYvtX2WdsPtZRdZvsu249mP588eOgAgEHkTfw/JemcpAcl/ZCkd0v6Hz22+V1JN7aV3SLp7oi4UtLd2XMAwBj1vDLH9pykByLimZJ+J2/FEfE+27W24pdKui57fJuke9X4owIAGJOeZ/wRcV7SR2wvj2B/l0fEmazeMzpkgNj2mu1t29vnzp0bwa4BAFL+a/GvkPSw7fvUmKtHkhQRL0kSVaPuTUmbkrSyssJAMgCMSN7EP6r5eR63fUVEnLF9haSzI6oXAJBT3j7+38z6+Id1p6SbJP189vMdI6gTANCHZH38tt8m6a8lfbXtT9p+tRoJ//m2H5X0/Ow5AGCMkvXxR8Qru7x0ff7wAACjNu4+fgDAhOX9Ipb32r5c0nOzovsigoFZAJhBeWfnfLmk+yR9p6SXS/qA7e9IGRgAII28XT3rkp7bPMu3vSTpL9SYqhkAMEPyztUz19a18+k+tgUATJG8Z/zvsf3nkt6WPf8uNSZqAwDMmF5ftv5Vasyv8xO2v03SN0iyGtfn18cQHwBgxHp117xZ0mclKSJuj4jXRcR/V+Ns/81pQwMApNAr8dci4oH2wojYllRLEhEAIKleif9Jh7x20SgDAQCMR6/E/0HbP9hemM27cypNSACAlHpd1XOzpDtsr+qLiX5F0jFJL0sYFwAgkUMTf0Q8Luk/2/5mSc1pmd8VEX+ZPDIAQBJ55+q5R9I9iWMBAIwBd98CQMmQ+AGgZEj8AFAyJH4AKBkSPwCUDIkfAEqGxA8AJUPiB4CSIfEDQMmQ+AGgZEj8AFAyJH4AKBkSPwCUDIkfAEqGxA8AJUPiB4CSIfEDQMkUNvHbHvly5MgR3XDDDVpcXDx0vcXFRdXr9X3x1Ov1fds112kv71VPe521Wk1zc3Oq1WoX1m2WN2O2ve/1Xtv3em1UBt1HP3GfPHny0Ocp2pUnTmCiImLql6uvvjr6IWniy9GjR2NraysiIra2tuLYsWMH1pmfn4/5+fnc9bTa2tqKSqWyb91KpRInTpw4UN76emtMnbbf2to69LVRGXQf/cbdaxl1u4ZtHzBKkrajQ06deFLPs8xi4pcU1Wo1IiKq1epI6mnVrc5ef0h6xVStVg99bVQG3ccgcQ9yfCfVPmCUuiV+N16bbisrK7G9vZ17fdsJo8nPts6fP6+5uTkNc5yb9bQatM5eMTWPXbfX2uMY1GH7P2wfg8Tdyyjb1TRo+4BRsn0qIlbayyfSx297x/aDtu+3nT+jz5jl5eV9P4etp1eZJM3Pzw8V0/Ly8qGvjcqg+xgk7kFjGcY4jiEwsE7/BqReJO1IWsy7/ix29dDHfzj6+IH0NE19/KkTf9bgkS/z8/Nx/fXXx8LCwqHrLSwsHPgF39ra2rddc5328l71tNdZrVbDdlSr1X1JvdnH3PzD0vp6r+17vTYqg+6jn7hPnDhx6POUiXgcxxA4TLfEP5E+ftufkPRPWYL77YjY7LDOmqQ1SVpeXr56d3d3vEECwIybqj5+SddGxHMkvVDSa2x/U/sKEbEZESsRsbK0tDT+CAGgoCaS+CPiU9nPs5LukHTNJOIAgDIae+K3fbHtS5uPJb1A0kPjjgMAyurIBPZ5uaQ7suuuj0j6vYh4zwTiAIBSGnvij4iPS/pP494vAKChsJO0AQA6I/EDQMmQ+AGgZEj8AFAyJH4AKBkSPwCUDIkfAEqGxA8AJUPiB4CSIfEDY1Sv11Wr1TQ3N6daraZ6vT7pkFBCk5irByiler2utbU17e3tSZJ2d3e1trYmSVpdXZ1kaCgZzviBMVlfX7+Q9Jv29va0vr4+oYhQViR+YExOnz7dVzmQCokfGJPl5eW+yoFUCpv4bede5ubmZFu1Wk0nT548MPjWaUCuXq9rcXHxQh2XXHLJhefz8/MHypvbNutvX29xcfHAQN+oBwKHra/bcWiWLS4u7mtrs/48+21fp9P7kMI497uxsaFKpbKvrFKpaGNjY6h6ux3fft/vafu8javOWdr/yHT6BvZpW66++up+v1l+ZMuxY8fi6NGjB8rm5+dHuh9JcfTo0dja2oqIiK2trahUKvter1QqF17v17D1ddq+07Fpr//EiRM999up7k51Ddr2ftqUer9bW1tRrVbDdlSr1aHr7va+5jnueeqZ1OdtXHXO0v4HIWk7OuTUiSf1PMskE/+4l2q1GhER1Wr10Nf7NWx93bbvtXT7A9m637x1D9r2Yds06v2OUrc25DnueeqZ1OdtXHXO0v4H0S3xu/HadFtZWYnt7e3c62df6ziTbOv8+fOam5tTp/em+Xq/hq2v2/aDat1v3roHbXs3k9rvKPX7vnRry7R93sZV5yztfxC2T0XESnt5Yfv4Z1VzoG/UA4HD1jfofufn53vWlzqGYeubpsHX9j7myy67rON6eY77MOW9pBjInvTg+KT3P1Kd/g2YtoU+/sZCH//s9/EPI+97QB9/GpPe/yBUpj7+rMG5F9sX+upOnDhxYPCt04Dc1tZWLCwsXKjj4osvvvB8bm7uQHlz22b97estLCwc+AD1OxDYa/1hBxa7HYdm2cLCwr62tv4R67Xf9nU6vQ8pTGq/g+jWx7ywsNAx5lF/fvo16vpS1TlL++9Xt8RfyD7+MmqfDkBqXCq4ubnJdAAFMYt9zJgs+vjHbNzX+zIdQPEVqo8ZE1XYxH/DDTf0faNWU6ebelpv1up0s1XrdrZ1/Phx7e7uKiK0u7ur48eP6+TJk7nj73ZjUbcbvw6bDiDvTUqj/GM1rroKc0NNm07tGtUNYEU9ZlLvto2j7TNxfDv1/0zb0m8f//XXX9/3wGpzkCbPgJ/aBmIj8g0U2s7VJ5g3htZYWscbWpeFhYVcA5j9Dgb2G3+KumZxsC2PXm0edpymiMcsonfbxtH2aTu+KtPgbr9Jv7lUq9W+blRKdRPSIDdLdUrwlUql6x+E9qXfG34GiX/Udc3iDTV5pGxXUY9ZRO+2jaPt03Z8uyX+Qg7uDnoDV3O7vMck1U1Ig9wsZVtvfetbtb6+rtOnT2t5eVkbGxs6fvx433X1G2+7UQ5CHlaX1Pm9mvXBzpSDuEUeIO7VtnG0fdqOL4O7OSwvL/c1UJbqJqRBBuuWl5e1urqqnZ0dnT9/Xjs7O1pdXc1dV783/PSKZRx1FXWwM2W7inrMpN5tG0fbZ+b4dvo3YNqWcfTx244TJ04k7ePP29c3SB9/PzfmdIqLPv7hDNP33uleglTtGuSYHda2Qdud6hr/Yfv4izaGojL18W9tbfWd+NsTSvsvYmtfeaebrZr7HdXNQM26mn+UusXcLZZB4hrlL+O46pqGG2qG+WXvtm3KG8n6OWYp/vCmTI7D3MQ4qrim4TPZVKrEP8jgaHOZtkGuaRsswkHDvEfT/v6mGFyf1jZPa1zD6Jb4C9nHP8xX2U3b1+DxdX2T1+u67GHeo2l/fw+Lb9DYp7XN44prGq7zL2Ti7zZjYR7TNggzM4NFBdWcCqP1Zry1tbV9v6zDvEfT/v6mGFyf1jaPI648n6ex6PRvwLQt/Xb15Ll23XbHWQ2n7UaWaRssKps8//6n6OOflvd31vr4hzGOuMbdnaQy9fEfNhgqNWbMHMVdkOMyK3EWUbfPku19643yqp5pe39n5aqeUUgdV97P06h0S/yFvIGrVqtpd3f3QHm1WtXOzs4II0PR8VnCKI378zRVN3DZvtH2x2w/ZvuWUdc/qsmsMHrTMLDVDz5LGKV+Pk9Jf1c6/RuQcpE0L+lvJX2FpGOSPiLpGYdtM8gXsUzrv5JlNq19u73wWcIo5fk8jep3RdPS1WP76yW9ISK+JXv++uwP0M9124YvYikGuk2AfEb1uzJNXT1fLunvWp5/Mivbx/aa7W3b2+fOnRtbcEhnWq/fBqZN6t+VSST+TlNnHvi3IyI2I2IlIlaWlpbGEBZSm9brt4Fpk/p3ZRKJ/5OSntby/KmSPjWBODBmDJQC+aT+XZlE4v+gpCttP932MUmvkHTnBOLAmK2urmpzc1PValW2Va1W+TJ4oIPUvysTuY7f9oskvVmNK3xujYhD/4wxuAsA/es2uHtkEsFExLslvXsS+waAsivkJG0AgO5I/ABQMiR+ACgZEj8AlMxMzM5p+5ykg/cv57Mo6R9GGM4soM3lQJvLYZg2VyPiwB2wM5H4h2F7u9PlTEVGm8uBNpdDijbT1QMAJUPiB4CSKUPi35x0ABNAm8uBNpfDyNtc+D5+AMB+ZTjjBwC0IPEDQMkUOvGn/lL3SbF9q+2zth9qKbvM9l22H81+Prnltddnx+Bjtr9lMlEPzvbTbN9j+xHbD9t+bVZe5DY/yfZ9tj+StfmNWXlh29xke972h22/M3te6Dbb3rH9oO37bW9nZWnb3OmLeIuwaIAvdZ+VRdI3SXqOpIdayn5R0i3Z41sk/UL2+BlZ279E0tOzYzI/6Tb02d4rJD0ne3yppL/J2lXkNlvSJdnjo5I+IOl5RW5zS9tfJ+n3JL0ze17oNkvakbTYVpa0zUU+479G0mMR8fGI+IKk35f00gnHNBIR8T5J/9hW/FJJt2WPb5P0rS3lvx8Rn4+IT0h6TI1jMzMi4kxEfCh7/FlJj6jxPc1FbnNExOeyp0ezJVTgNkuS7adK+q+S/ndLcaHb3EXSNhc58ef6UvcCuTwizkiNRCnpKVl5oY6D7ZqkZ6txBlzoNmddHvdLOivprogofJvV+IKmn5R0vqWs6G0OSf/X9inba1lZ0jZP5ItYxiTXl7qXQGGOg+1LJP2xpJsj4jN2p6Y1Vu1QNnNtjognJF1l+8sk3WH7mYesPvNttv1iSWcj4pTt6/Js0qFsptqcuTYiPmX7KZLusv3RQ9YdSZuLfMZfti91f9z2FZKU/TyblRfiONg+qkbSr0fE7VlxodvcFBH/LOleSTeq2G2+VtJLbO+o0TX7X2xvqdhtVkR8Kvt5VtIdanTdJG1zkRN/2b7U/U5JN2WPb5L0jpbyV9j+EttPl3SlpPsmEN/A3Di1f4ukRyLiTS0vFbnNS9mZvmxfJOkGSR9VgdscEa+PiKdGRE2N39e/jIjvUYHbbPti25c2H0t6gaSHlLrNkx7RTjxa/iI1rgD5W0nrk45nhO16m6Qzkv5djTOAV0takHS3pEezn5e1rL+eHYOPSXrhpOMfoL3foMa/sw9Iuj9bXlTwNj9L0oezNj8k6Wey8sK2ua391+mLV/UUts1qXHX4kWx5uJmnUreZKRsAoGSK3NUDAOiAxA8AJUPiB4CSIfEDQMmQ+AGgZEj8KDzbC9nMh/fb/n+2/77l+bG2dV9ve7Wt7FW2fyN7PGf7NjdmSO1467Dt65ozSx4S0w/b/t5h2wYMoshTNgCSpIj4tKSrJMn2GyR9LiJ+ucvqL5D08k4vZIn+t9SYMO37YohroSPitwbdFhgWZ/xAxvaXSjoWEee6rPJratxY870RcT676/JW2x/M5o/fN/tr9t/BTvMO3KzsMduX236D7R/Pyu61/QvZ/Pt/Y/sbEzURkETiB1rdoMZdkp18t6SrJb0iIv4jK1tXY1qB50r6Zkm/lN12L0mKiPNq3Gr/Mkmy/XWSdiLi8Q71H4mIayTdLOlnR9AWoCsSP/BFN0r6sy6vfUhSVfvnPn+BpFuyqZPvlfQkSctt271d0ndlj1+RPe+kOfHcKUm1PmIG+kbiB77oGnWf8OqjavT9v93212ZllvTtEXFVtixHxCNt2/21pK+yvaTGl2ncrs4+n/18Qoy9ITESPyApS+YfjcYc+B1FxPsl/bCkd9lelvTnkn60eXWP7Wd32CbUmGr3TWrMLvrpFPED/eDMAmh4oaT39FopIt6Znb2/R9LzJf2MpAey5L8j6cUdNnu7GtOEv2pUwQLDYHZOQJLtu9S4WufMpGMBUiPxA0DJ0McPACVD4geAkiHxA0DJkPgBoGRI/ABQMiR+ACiZ/w8FthJarOmnLQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.close('all')\n",
    "plt.figure()\n",
    "plt.plot(T_plot_1,CL_plot_1,'ok')\n",
    "#plt.gca().invert_xaxis()\n",
    "plt.xlabel('T / Kelvin')\n",
    "plt.ylabel('Correlation')\n",
    "#plt.xlim([1,0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "5252c135",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABSNklEQVR4nO3dd3gU1frA8e+kNxISaiDAhiI7dAEBAZWiUoKgKIiV5sWC4LWgwa6/q+ZauNgRQUGlCwoaQAXEhkDobRakLCmEJJDes9n5/bFLTCSkQDab8n6eJ8/unGnv4Xrz5sycoui6jhBCCFHTuDg7ACGEEKI0kqCEEELUSJKghBBC1EiSoIQQQtRIkqCEEELUSJKghBBC1Ehujry4ZlQfBx4AdOAgMBnwAVYABsAMjFdNWor9+NnAVKAQmKmatB8cGZ8QQoiay2EtKM2otgRmAr1Vk9YFcAUmAOHAZtWkdQA227fRjGon+/7OwHDgI82oujoqPiGEEDWbox/xuQHemlF1w9ZyOgOMARbb9y8GbrV/HwMsV01anmrSTgHHgT4Ojk8IIUQN5bBHfKpJi9OM6ttANJAD/KiatB81o9pMNWnx9mPiNaPa1H5KS2B7sUvE2stKUBRlGjDNvtnLx8fHUVUQQohaKzs7W9d1vVb3M3BYgtKMaiC2VlEokAqs0ozqvWWcopRSdtE8TLquzwfmA/j6+upZWVlXHqwQQtQxiqLkODuGK+XI7HojcEo1aUmqSSsA1gD9gQTNqAYD2D8T7cfHAq2KnR+C7ZGgEEKIesiRvfiigX6aUfXB9ohvKLALyAImAhH2z7X249cBSzWjOgdoAXQAdjowPiGEEDWY4sjZzDWj+gpwJ2AB9mLrcu4HrARaY0ti41STlmw//jlgiv34f6smbUNZ15dHfEIIUTpFUbJ1Xfct6xhDeORw4F1svawXmCPCIqoluApyaIJyNElQQghRuvISlCE80hU4BtyE7RVLFHCXOSLsSDWFWK5a3cNDCCHEZesDHDdHhJ00R4TlA8uxdWyrMSRBCSFE3eSmKMquYj/T/rG/JRBTbLvUoT3O5NCpjoQQQjiNRdf13mXsr9DQHmeSFpQQQtRQBYVWR16+xg/tkQTlAIqicN999xVtWywWmjRpwqhRo4rKNmzYQO/evVFVFaPRyFNPPXXJ640ZM4Zrr73WoTEXt2/fPtavX19t9xOivilISOD0ffeT9t33pe7Pzrew4LeTDIjYwr6YVEeFEQV0MIRHhhrCIz2wzYW6zlE3uxySoBzA19eXQ4cOkZNjG8j9008/0bLl3492Dx06xKOPPspXX32FpmkcOnSItm3blnqt1NRU9uzZQ2pqKqdOnaqW+CVBCeE4uq5z9uVXyI6K4sysWSR98CEXelNn5Bbw4c/HGfjfn/lPpEb7pn64uZT2JO7KmSPCLMCjwA+ABqw0R4QddsjNLpN0M3cAPz8/Zs6cSc+ePbnjjju4//776dy5M7/99hvff/89999/P4MGDWLKlCnlXmvhwoXs3r2bZs2a4eHhwezZswE4fvw4Dz30EElJSbi6urJq1SratWvHm2++yZdffomLiwsjRowgIiKCEydOMH36dJKSkvDx8eHTTz/FaDQyadIkvLy8OHz4MAkJCcyZM4ebb76Z9u3bk5OTQ8uWLZk9ezbNmzfnscceA2ytw19//ZUGDRo49N9QiLoq7ftIzjz1FE2eeIL8EydIW7sWRt/G+iH3sXh7NOm5FgZ1bMKMIe3p1Sbosu9TkXFQNV2d7iRx9vXXydNMVXpNT9VI82efLfe4CRMm8OqrrzJq1CgOHDjAlClT+O233wBbC+rJJ5+s0P2WLVvGSy+9RLNmzbjjjjuKEtQ999xDeHg4t912G7m5uVitVjZs2MC3337Ljh078PHxITk5GYBp06Yxb948OnTowI4dO3jkkUfYsmULAGazmV9++YUTJ04wePBgjh8/zquvvsquXbv44IMPALjlllv48MMPGTBgAJmZmXh5eVX6300IAZbkZBJeew2v7t1oNHUK57ILWB7QnZXnPcjZepKbrwpixs2d6BoS4OxQa4Q6naCcqVu3bpjNZpYtW8bIkSMv6xoJCQkcP36cgQMHoigKbm5uHDp0iDZt2hAXF8dtt90GUJQwNm3axOTJk7kww3tQUBCZmZls27aNcePGFV03Ly+v6Pv48eNxcXGhQ4cOtG3bFpPp4oQ+YMAAnnjiCe655x7Gjh1LSEjIZdVHiPou4T+vUZiZidfsl3h1vYllO6PJs/gzrAXcuvpd2u9RaDV4HiAJCup4gqpIS8eRRo8ezVNPPcXWrVs5f/58UXnnzp3ZvXs33bt3L/P8FStWkJKSQmhoKADp6eksX76cp59+utTjdV1HUUo+r7ZarTRs2JB9+/aVes4/j//nNkB4eDhhYWGsX7+efv36sWnTJoxGY5mxCyFKytiyhWNbt/Pdnc/y/arTFOo6t13dkocHtaNdEz+yB7cgdvp0zBPuIuT99/DtI8vhSScJB5oyZQovvvgiXbt2LVE+a9YsXn/9dY4dOwbYksicOXMuOn/ZsmVs3LgRs9mM2Wxm9+7dLF++HH9/f0JCQvj2228BW4soOzubm2++mc8++4zs7GwAkpOT8ff3JzQ0lFWrVgG2JLZ///6ie6xatQqr1cqJEyc4efIkHTt2pEGDBmRkZBQdc+LECbp27cozzzxD7969S21lCVFf5B49xslbbyNj06YKn3P81Fme+HInU28K5/scf8b1DmHrU4N4e1x32jXxA8Cn59UYVq7ArVEjoqc+QOo33zqoBrWHJCgHCgkJKepcUFy3bt2YO3cud911F6qq0qVLF+Lj40scYzabiY6Opl+/fkVloaGh+Pv7s2PHDr788kvee+89unXrRv/+/Tl79izDhw9n9OjR9O7dmx49evD2228DsGTJEhYuXEj37t3p3Lkza9euLbpmx44dueGGGxgxYgTz5s3Dy8uLwYMHc+TIEXr06MGKFSuYO3cuXbp0oXv37nh7ezNixAgH/YsJUbPpuk7C66+TZzIRO2MmyV9+VebxprPpzFi2l5s+2cUvjY3cowbwy9ODee22rrQKunixVY9WrTAsX4ZPr17Ez55N4ty56FaHjoWq0aQXXz02adIkRo0axR133OHsUISoFTI2byZ2+qM0nfUU2Xv2krl5M0GTJtH06VkoLn//vb8/JpUPfz7Oj0cS8HGFMG0LU65thTrr3xW6j15QQPwrr5D29WqC33iDhrfdWulYpRefEELUE3p+PolvvoVHu3YETZxI0KRJJLwRQfKiRRScOUOTN97gh+OpLPrjFHuiU/H3cmPGdW0Y/O4zBLhYCZ3xWoXvpbi7E/x//4ffwIE0uOkmB9aqZpMWlBBCVEDy4sUkvBFBq/mf4Hf99YDtkd+JhV/w5booIq+6nnNuvrRp5MPEaw2M6x1C1jtvkfLVV7T56kt8evWq1nilBSWEEPWAJSWFpA8/wnfAAHyvuw6AI2fSWbTtFN+am5LfaQRXn/uLf6dp3P7gbLxDDWTv2UPKV18RePfd1Z6c6gppQQkhRDnO/uc1UpYupfWaNfxmCeDzP06x41Qy3u6ujO3Zkkn9DbSMO07sI4+AotBy7lzOvvwyel4ebb9bh4tv9Tdk6kILShKUEEKUIe/kSfaPvZNfR03lG7+riEvNoWVDbyb2b8OdvVsT4ONedGy+2Uz0tAcpiI4GoNWCBfgNHOCUuCVBOZkkKCFEZaVv2EBCxH9pNns2/sOHlXnsXwkZfBDxBRvdWpDn6kHf0CAmDwjlRrUpbq6lj9KxJCdzZtbTeF51Fc2eKX1QfXWQBOVkkqCEEJWRazJhnnAXWK3o+fk0fuQRGj86vUQXcatV5+ejiSzaZua3v87hXljAyMACHpw0jE4t/J0YfeXUhQTlsE4SmlHtCKwoVtQWeBH4wl5uAMzAeNWkpdjPmQ1MBQqBmapJ+8FR8Qkh6hdLSgqxj87ANSAAw7KlJL3/Aec++oi8v/6iRcQbZLt58vWuGBZtM2M+n01zf0+mJu5k5Jk9XL1uNS4eHs6uQr1TLS0ozai6AnFAX2A6kKyatAjNqIYDgapJe0Yzqp2AZUAfoAWwCbhKNWmFl7qutKCEEBWhWyzETHuQ7F27aPPVl3h364au6yQvXsyBDxaw/uowNgT3IDPfytWtGzJ5QCj9TL9z/qWXaDl3brmPAmuiutCCqq6pjoYCJ1STdhoYAyy2ly8GbrV/HwMsV01anmrSTgHHsSUrIYS4Iklz55K1bRvNX3qpKDntPJXMs65dmHrjbNYEduaamAMsGxzEN48MIKydP6nvvYd37140GHazs8Ovt6prHNQEbK0jgGaqSYsHUE1avGZUm9rLWwLbi50Tay8rQVGUacA0AA9pcgshypG+fj3nFywk8O678LxlDCujYlj8p5nDZ9Jp6OPOw4PbcWcLFwpmfUz+U0tISXuegphYCs+fp9m8eaXO8C+qh8MTlGZUPYDRwOxyDi3tv4KLnj/quj4fmA+2R3xXHKAQos7KPXqUM889T1LfQazqfhtfv76J9FwLVzXz442xXbm1R0u8PVwBKFyxnLgnn+Lsiy+BiwsBY8bg3bWLk2tQv1VHC2oEsEc1aQn27QTNqAbbW0/BQKK9PBZoVey8EOBMNcQnhKiD8pJTWPb8XNb1ncrehgbcd8YwvEsw9/VrwzWGwItaRq4BAbSa9zGJ78wh/YeNNHn8384JXBSpjgR1F38/3gNYB0wEIuyfa4uVL9WM6hxsnSQ6ADurIT4hRB2SkJ7Lsu2n+eqng5xrN4pgH1dmXdee8b1b0aSBZ5nnKm5uNHvmadvs5PJoz+kcmqA0o+oD3AQ8WKw4AlipGdWpQDQwDkA1aYc1o7oSOAJYgOll9eATQogLdF1n24nzfLX9ND8eScBqtdLr3Gme7x/KLVNvw9WlcslGklPNIAN1hRC1Vlp2AV/viWXJjtOcTMoiwNXKsNg9DDv4E51G30jwyy87O0SnqQvdzGU2cyFEraLrOvtj01i2I5q1++PILbDSxT2HWUd+YODx7QT06EajN1/Gb8gQZ4cqrpAkKCFErZCSlc83e+NYERXD0YQMvN0UhhWe5abfVtAuOYYGN95I0Ctf4HP11c4OVVQRSVBCiBrLatX548Q5VkTF8OPhBPILrXRt5sMs6zH6fvsFfq46AbfdSqNJn+BhMDg7XFHFJEEJIWqcM6k5rNoVy8pdMcSl5hDg7c7d17RkWOxuAj/7AAoKCJo8iaBJE3Fr1MjZ4QoHkQQlhKgR8i1WNmkJrIiK4de/ktB1GNi+Mc+MMDLg/DFS//s0+adP4ztkCM3Cn8GjdWtnhywcTBKUEMKp/krIYEVUDGv2xpGclU9zfy9mDG7PuN6taJaRRMLrr5L4yy94GAy0+nQ+fvYl10XdJwlKCFHtYpKzWX8wnvUH49kfm4abi8KNajPu7NOK6zs0wdVF4fznizg5Zw6KuztNZ80i6L57UWT+zXpFEpQQolpEn88m0p6UDsalAdC1ZQDPjjQytmcIjf3+nuUhZdkyEv/7X/xuHErzF1/EvWnTS11W1GEyUFcI4TBZeRa+3h3Lqt0xHIpLB6B7SAAjuwYzsmswrYJ8LjonfcMG4p54Er9Bgwh5/z0UN/k7+nLUhYG6kqCEEFUuLjWHxdvMLNsZTUauhW4hAdzSrQUjujYnJPDipHRB5u9/EPPww3h370brBQtw8fKqxqjrlrqQoORPEyFEldB1nT3RqXz2+yk2Hj4LwIguzZkyMJSerQPLPT9n/35iZ87Es107Wn30kSQnIQlKCHH5dF3neGImm7RENhyK50BsGv5ebjxwXSj3X2ugZUPvCl0n78QJYqY9iFujRrT+dD6u/v4OjlzUBpKghBCVUlBoZeepZDZpCWzWEolOzgagcwt/Xh3Tmdt7huDrWfFfLQVnzhA99QFwd6f1wgW4NWniqNBFLSMJSghRLl3XORCbxtId0aw/FE9GrgUPNxcGtGvEtOvbMlRtSnBAxVpLxVlSUoie+gDWrCzafPmFDL4VJUiCEkJcUkZuAWv3nWHpjmiOxKfj4+HKyK7B3NypGQM7NMbH4/J/heiFhcTNmEnBmTO0XrgAL6OxCiMXV8oQHvky8C8gyV70rDkibH11xiAJSghRwoXlLFZERbN23xmy8wtRg/35v1u7cGuPFjTwcq+S+yQvWkz2rl0Ev/EGPr17V8k1RZX7nzki7G1n3VwSlBACgPi0HL7ZG8eaPXEcT8zE292VW7oHc3ffNnQPCajSVWZzjx0jae5cGtx0IwG3jqmy64q6RcZBCVGPZedb+OHwWVbvjuOPE+fQdbjGEMjtPUMY2S0Y/ypqLRWn5+dz6s4JWBISaPv9d7gFBVX5PQQoipIPHCxWNF/X9fkVPd/+iG8SkA7sAp40R4SlVGWM5ZEWlBD1UHpuAV9sM7Pg91OkZhfQKsibmUM6MLZnS9o0uvTYzvQffyRlyVJ8BwygwY1D8WzbttL3TvroI/I0jZAPP5Dk5FgWXdfLfHZqCI/cBDQvZddzwMfA/wG6/fMdYEpVB1kWh7agNKPaEFgAdMFWySnAUWAFYADMwHjVpKXYj58NTAUKgZmqSfuhrOtLC0qIyknLLuDzbaf47PdTpOdaGGJsyrTr29I3NKjcR3i6xcKJYcOxpKSgZ9u6lnsYDPgNHUKDoUPx7t4dxdW1zGvk7NuH+e57CBgzhhZvvF5l9RIXq8qZJAzhkQbge3NEWJequF5FuTj4+u8CG1WTZgS6AxoQDmxWTVoHYLN9G82odgImAJ2B4cBHmlEt+792IUSFJGbk8s6PRxn43y3M3fQX/do24rtHB/LZpGvo17ZRhd4vpW/8gYK4OFq+9Sbtt/5MsxdfwL1lS5K/+JLTd9/DX9ffwLlPP8Wal1fq+dacHM48E45b82Y0e3Z2VVdRVDFDeGRwsc3bgEPVHYPDHvFpRtUfuB7bM0xUk5YP5GtGdQwwyH7YYmAr8AwwBliumrQ84JRmVI8DfYA/HRWjEHVZXGoOGw+d5YdDZ4k6nQzAyC7BPDqkPWpw5WZq0HWd8wsW4NGuHX6DB6O4uBB0990E3X03hRkZZP76K2lr15L0zhxSli2j6eNP4B82EsXl77+BE99+h/zTp2m9aBGuDRpUaV2FQ7xpCI/sge3plxl4sLoDcOQ7qLbY+s9/rhnV7sBu4DGgmWrS4gFUkxavGdUL8+i3BLYXOz/WXlaCoijTgGkAHrI2jBAlJGXksXpPLBvs6ywBGJs34LGhHbilewvaNfG7rOtm/f47eSYTwa+/XiLpALg2aEBAWBgBYWFkbd9OwptvcmbWLJK/+IJm4c/g06sXWdu2kbJkCYH334dvv75XXE/heOaIsPucHYMjE5Qb0BOYoZq0HZpRfRf747xLKO0Zw0UvyOy9UOaD7R1UVQQqRG13KC6Nz/8w893+M+QXWukWEsDTwzsyokswoY2v/DXE+U8X4Na8OQGjwso8zrdfP0K//pq0tetImjuX0/fcS4ObbiLn4EE82ral6RNPXHEsov5wZIKKBWJVk7bDvv01tgSVoBnVYHvrKRhILHZ8q2LnhwBnHBifELVaoVXnpyNn+ewPMztPJePj4cqEPq2Y2N9w2S2l0uTs30/2zp00DX+mQivaKi4uNLztVvyHD+P8559zfsFC9Lw8DMuXyQzlolIclqBUk3ZWM6oxmlHtqJq0o8BQ4Ij9ZyIQYf9caz9lHbBUM6pzgBZAB2Cno+ITorbKsxSyZk8c8345wenz2bRs6M1zI1XGX9OKAO+qH7d0fsECXAICCBw3rlLnuXh70+SRRwgcPx5LUhJeqlrlsYm6zdHjoGYASzSj6gGcBCZj6zm4UjOqU4FoYByAatIOa0Z1JbYEZgGmqyat0MHxCVFr5OQXsjwqmvm/niQ+LZduIQF8fE9PburUDDdXx3TIzTt5koxNm2n88EO4+F7eo0K3xo1xa9y4iiMT9YHMJCFEDZeWXcDSndEs/P0k5zLz6WMIYvqQ9lzfoXGVTj9UmjPPPUd65Hrab9ksg2prGVlRVwjhMEfPZrBom5lv9saSW2Dlug6NeXRwe/q2bVQt9y84e5a0dd8ROG6cJCfhFJKghKhBLIVWNmmJLN5m5s+T5/F0c+HWHi2Z2N9ApxbVu8ps8uIvwGolaMrkar2vEBdIghKiBohJzmZFVAyrdseQkJ5HiwAvnhluZMI1rQj0ddx4P2t+PlgsuPj4lCgvTEsjdcUK/EeMwCMkxGH3F6IskqCEcBJLoZUfjySwbGc0vx8/hwLccFUTXh3TmqHGplXe8cGalUXu0aPkHj5C7pEj5GoaecePFyUotyZNcG3SGLfGTbCmp2HNzqbRvx6o0hiEqAxJUEJUs5z8QlZERfPpb6eIS82hRYAXjw3twPjerWjRsPLLppcn99gxkt6ZQ+avv4K9U5Rro0Z4deqE3w034OLrS+H5c1iSkrAknSPv6FEsSUn4jxqFV8eOVR6PEBUlCUqIapKanc8Xf55m0TYzyVn59G4TyMujOzPE2BRXl6rvjVeQkEDSe++R9s23uPj50WjaNLx7dMerU2fcmjZxeA9AIa6UJCghHMx0Np0v/jzNN3viyCkoZKixKQ8Nasc1Bsf0jCvMzOL8wgUkf74IvbCQoPvuo9FDD+IWGOiQ+wnhKJKghHCAgkIrPxw+yxd/nmbnqWQ83VwY3b0FD1zXlo7NHTeTd8bWrcQ/9zyF58/jP3IkTZ54XDo5iFpLEpQQVSgtp4BlO6NZ9IeZs+m5tAry5tmRRsb1cmxvPL2wkKQPPuD8x/PwVFVazfsY765dHXY/IaqDJCghqkBMcjaf/XGKlVExZOUXMqB9I167rQuDOjrm/VJxluRkzjz1FFnb/iTgjttp/vzzMimrqBMkQQlxmQoKrWw9msSKqBi2mBJwURRGd2/B1OtC6dwioFpiyN67l7jHn6AwOZng//wfDe+4o1ruK0R1kAQlRCWdTMpkxa4Y1uyJIykjjyYNPHnwhnbcf20bggOqvpv4paSuXk38Sy/jHhyMYfkyvDp1qrZ7C1EdJEEJUQGWQiubTYl8+edpfj9+DlcXhSHGpozv3YrBHZs4bDbxS8k9doz4l1/B55rehMydi2tA9bTYhKhOMpu5EGVIycpn6c5olmw/zZm0XIIDvLinb2vGX9OKpg2c855HLyzEfNfdFMTE0Dbye5nIVZRKZjMXoo6KS81hwW8nWb4zhpyCQga2b8xLozs7ZAqiykr+8ktyDxygxdtvS3ISdZokKCGKOXo2g09+OcG6/WcAGN2jBQ9e385hY5d0qxXFpeIJLz8mhqS57+I3aBD+YSMdEpMQNYUkKCGAKHMy87aeYLMpER8PV+6/1sDU60Jp6YC58S7I2raNuMefIDjiDRoMHlzu8bquE//CiyhubjR/+SWZqkjUeZKgRL1ltepsNiUy75cT7D6dQpCvB4/feBX3X9vGoYNqwZZsEue+S2FaGnFPPEmbxYvw7tatzHNSv/6a7O3baf7KK7g3b+7Q+ISoCaSThKh38i1W1u6LY/6vJ/krMZOQQG/+dV1bxvduhbeHa7XEkPnb78T86180njmDtG++xZqZiWH5MjzatCn1+IKEBE6GjcKrUydaL/q8Uo8FRf1UFzpJODRBaUbVDGQAhYBFNWm9NaMaBKwADIAZGK+atBT78bOBqfbjZ6om7Yeyri8JSlRGTn4hy6Oimf/rSeLTcjE2b8DDg9oR1jW4Wjs+6LrO6bvvoeDsWdr9sJGCuDhO33U3Lv7+GJYvu6jjgzUri7gnnyJr+3barv32kklMiOLqQoKqjkd8g1WTdq7YdjiwWTVpEZpRDbdvP6MZ1U7ABKAz0ALYpBnVq1STVlgNMYo6LDPPwlfbT7Pgt5Ocy8ynjyGI18d2ZdBVzllyInvHDnL27qXZC8/j4uGBZ2goreZ9zOlJk4l56GHaLPocxdubnL17SV29mvQNG9Gzs2k2O1ySk6hXnPEOagwwyP59MbAVeMZevlw1aXnAKc2oHgf6AH86IUZRB+i6zspdMURsMJGSXcB1HRrz6OD29G3byKlxnfvwI9yaNCkxLZF3jx60nPMOsY/OIHryFArT0sg3m1F8fPAfMZyGt9+OT8+eToxaiOrn6ASlAz9qRlUHPlFN2nygmWrS4gFUkxavGdWm9mNbAtuLnRtrLytBUZRpwDQADw/HvsgWtdfJpExmrznIjlPJ9DEEMXukkatbO389pOyoKLKjomj27GxcPD1L7GswZAjNX3yBs6+8infPngT/61/4Dx+Gi2+tfkojxGVzdIIaoJq0M/Yk9JNmVE1lHFvas5aLXpDpuj4fmA+2d1BVE6aoK7LyLHz620k+2noCLzcXIsZ2ZXzvVrg4eEbxijr38ce4Nm5Mw3HjSt0fOGECAaNH4+LjU82RCVHzODRBqSbtjP0zUTOq32B7ZJegGdVge+spGEi0Hx4LtCp2eghwxpHxibojt6CQpTui+fDn45zPyiesWzAv3dLJadMRlSZ7716ytv1J01mzcPG+9PgqSU5C2DgsQWlG1RdwUU1ahv37zcCrwDpgIhBh/1xrP2UdsFQzqnOwdZLoAOx0VHyibrAUWlm9J5Z3N/3FmbRc+rdrxFPDOtKzGh/n6VYr1vR0XBs2vPQxum5rPQUGEjjhzmqLTYjarNwEZU8uOapJs2pG9SrACGxQTVpBOac2A77RjOqF+yxVTdpGzahGASs1ozoViAbGAagm7bBmVFcCRwALMF168IlLsVp1Ig/G87+fjnHyXBY9WjXkrXHdGdC+cbXHkvz55yS+9TaeqkqDIUNoMHQInqoKFgvZu3eT+fPPZGzdSsHpaJo88YS8UxKigsodB6UZ1d3AdUAgtk4Mu4Bs1aTd4/jwyibjoOqnP46f47VIjSPx6XRs1oCnhnXkRrWpU7qM61YrJ266GcXTE9fAQHL27AFdxy04GGtGBtbMTBQPD3z69aXBkKE0vON2FDeZwEU4Xn0ZB6WoJi3b3uJ5XzVpb2pGda+jAxPin44nZvD6ehNbTImEBHoz984ejO7ewqkdILL+/JOCuDhavPM2AWFhWM6fJ3PrVjJ/+RWXBn40GDwY3/795b2SEJehQglKM6rXAvdgm+WhoucJccWsVp3fjp9j2Y5oftIS8HF3JXyEkUn9DXi5V8+0RGVJ/fprXAMCaHDTTQC4NWpEw9tvp+Httzs5MiFqv4okmseA2cA39vdEbYGfHRuWqO8shVaW7LBNSxSXmkOQrwdTB4by4PVtaeTnWf4FqoElOZmMTZsJuvsuXGRMnhBVriIJqplq0kZf2FBN2knNqP7mwJhEPbfLnMwLaw+jxafTxxBE+AgjN3duhqeb81tMxaWtXQcFBSVmhBBCVJ2KJKjZwKoKlAlxRZIy8nhjg8aaPXG0CPDi43t6MrxL8xq57pGu66SuWoV3jx54dujg7HCEqHKG8MhxwMuACvQxR4TtKravxMTe5oiwMif2vlyXTFCaUR0BjARaakb1vWK7/LF1AxeiSlgKrXy5/TRzfjxGrqWQ6YPbMX1we3w8au6rzpy9e8k/eZLg1/7j7FCEcJRDwFjgk+KFhvDIiyb2NoRHXmWOCKvyYUFl/QY4g61L+Whgd7HyDODxqg5E1E/7Y1IJX3MQLT6d6zo05pXRnWnbxM/ZYZUrddXXuPj64j98uLNDEcIhzBFhGoAhPPKfu8YAy80RYXnAKUN4pMMm9r5kglJN2n5gv2ZUl2KbJ8+IbW68o6pJy6/qQET9kpVn4e0fj7J4m5kmDTxr9OO8fyrMyCB9wwbbnHky6FbUXG6Kouwqtj3fPpfplarQxN6G8MjexR8LXo6KPEO5CVsT7wS2RBWqGdUHVZO24UpuLOqvn02JPP/tIc6k5XBv3zbMGt4Rfy93Z4dVYemRkei5uTQcJ50jRI1m0XW9d1kHGMIjNwHNS9n1nDkibG0p5VDBib2BTw3hkX7AMmwtriNlRluKiiSoOdgWHTwOoBnVdkAkIAlKVEpSRh6vfHeY7w/Ec1UzP75+6Fp6tQkq/8QaxJqTQ8qSJXgajXh16eLscIS4IuaIsBsv47QKTextjgi72hAe2RHb+6qvDeGR+fydrE5X5EYVSVCJF5KT3Un+noFciHKl5xaw4NeTLPz9FAWFOk/edBUP3tAOD7fqW2a9Kljz84md+Rh5x08Q8sH7teJxpBAOsA5YagiPLHdib3NE2FHgFeAVQ3hkd2zJaoshPPKsOSJsQHk3qkiCOqwZ1fXASmzNuHFAlGZUxwKoJm1NBa4h6qFCq85X20/zv03HSM0uIKxbME/d3JHQxrXvvY1usXDmySfJ+u03gl/7Dw2GDnV2SEI4lCE88jbgfaAJEGkIj9xnjggbZo4IO2wIjywxsXd5PfgM4ZEuQFNsk4j7AkkViaEik8V+XsZuXTVpUypyI0eQyWJrrmMJGTyz+gB7o1MZ2L4x4SOMdGkZ4OywLotutRI/ezZpa9fR7NlnCbr/PmeHJES5asJksYbwyOuAu4BbsXVbXw6sNkeEpVXk/HITVE0mCarmycyz8MGW4yz8/SR+nm68eEsnbu3RskY/DtN1/ZLx6brO2VdfJXXZcpr8+980fujBao5OiMvj7ARlCI+Mwbak0nJgpTkiLKGy16jIelBXAR9jm/Koi2ZUuwGjVZMmIxRFEV3X+XZfHG+sN5GYkcfYni15dqRK4xoyb96lnPv0U1KWLKXFf/+Lb98+Jfbp+fnEv/AiaWvX0uhf/5LkJETlDCzeGcIQHulrjgirVIuiIm+pP8U2tVEBgGrSDmB70SUEAAdj07hj3p88vmI/zQO8WPNIf+aM71Hjk1PaunUkvTOHwrQ0oqdOJWXV37N3FWZkEPPQQ6StXUvjmTNo8oSMTReiMi4kJ0N45LWG8MgjwIWBv90N4ZEfVeQaFUlQPqpJ+2cPDZnqSHAuM4/w1QcY/eHvnD6fxZt3dOPbRwZU63Lrlytrx07OPPc8Pn370n7zJnz79ePsCy+S8EYEBWfOcPre+8jaGUXwG2/Q5JFHavQjSiFquLnAMOA8gDkibD9wfUVOrEgvvnP2sU86gGZU7wDiLytMUScUWnW++NPMnJ+OkZNfyNQBocy8sUOtGWybd+IEsTNm4NG6NSHvvYtrQACt5n1MQsR/SV68mOSlS3Hx8KDVJ/PwG1BuT1ghRDnMEWEx/5gyqULz9lUkQU0H5gNGzajGAaeAeysdoagTjiVk8PTXB9gXk8p1HRrz0i2daN+0gbPDqrD82Fhipj2I4u5Oq08+wTXA1rNQcXOj+fPP4dm+HalrviH41VfwMhqdHK0QdUKMITyyP6AbwiM9gJnYH/eVp8K9+DSj6gu4qCYtozKRaUbVFduks3GqSRulGdUgYAVgAMzAeNWkpdiPLTGFu2rSypzCXXrxVZ+CQivztp7g/S3H8fV05eXRnRndvUWtefRlzc3l/MKFnJ//KYqrK60XL8K7a1dnhyWEwzi7F98FhvDIxsC7wI3Ypkn6EXjMHBF2vrxzK9KL73XgTdWkpdq3A4EnVZP2fAXjewxbtvS3b4cDm1WTFqEZ1XD79jOaUb1oCnfNqF6lmrQqn8JdVM6huDRmfX0ALT6dUd2CeXl05xrfAeICXdfJ3LLF9m4pNhb/kSNoOmsW7sHBzg5NiHrBHBF2Drjncs6tyCO+EapJe/bChmrSUjSjOhIoN0FpRjUECANeA56wF48BBtm/Lwa2As/Yy5erJi0POKUZVYdN4S4qJregkPc2/8Unv54kyNeDT+7rxbDOpc0rWfMUnD1L2rffkrrmGwqio/Fo347Wixbh26+vs0MTol4xhEeGAjOwPTUryjnmiLDRlzrngookKFfNqHraEweaUfUGKvrn81zgaaD4S4pmqkmLB1BNWrxmVJvayys0hbuoHrtPJ/P01wc4kZTFuF4hPB/WiQCfmt8JwpqdTfzzz5O+YSPoOj59+tBk+iP4jxyJ4l7z4xeiDvoWWAh8B1grc2JFEtRXwGb7lEc6MAVby6dMmlEdhW2i2d2aUR1UgftUaAp3RVGmAdMAPDw8KnBZURnZ+Rbe+uEoi7aZaRHgzeIpfbjhqibODqtCbGOXHiZn714aPfAADcePw6NVq/JPFEI4Uq45Iuy98g+7WLkJSjVpb2pG9QB/v+D6v/I6L9gNAEbbHwd6Af6aUf0KSNCMarC99RTM3zOjV2gKd/uCW/PB1kmiAnGICtp2/BzPrDlATHIO91/bhqeHG/HzrLnLrhdnSUkh5oF/kXv0KC3nvCMr3QpRc7xrCI98CVvniLwLheaIsD3lnViRThKhwFbVpG20b3trRtWgmjRzWeepJm02thkosLegnlJN2r2aUX0LmAhE2D8vLIq1DliqGdVyp3AXVUvXdeb8dIz3txzH0MiHFdP60bdtI2eHVSF6fj45Bw8S/9JLFMTE0urDD/C74QZnhyWE+FtX4D5gCH8/4tPt22WqyJ/Hq4D+xbYL7WXXVC7GIhHASs2oTsU2keA4ANWkHdaMaokp3KUHn+PlWQoJX32Qb/bGMa5XCK+O6YK3h6uzw7oka14euQcOkBUVRfbOKHL27UPPzcXF15dW8+dfNJ+eEMLpbgPamiPC8it7YkUSlJtq0oourJq0fM2oVurlj2rStmLrrYdq0s4DpS6mo5q017D1+BPVwHwui5nL93IgNo0nb7qKR4e0r7HjmvJjY4l//gVy9uxBz88HRcHTaKTh+HH4XHMNvtdcg2vDhs4OUwhxsf1AQy5joduKJKgkzaiOVk3aOgDNqI4BzlX2RqJm+XZvHM99cxBXF4V59/ZkeJeaPS4o48efyN6+ncD778O3Xz98evUqmgVCCFGjNQNMhvDIKEq+g6qSbuYPAUs0o/oBtk4SMdieJ4paKCvPwgtrD7FmTxzXGAKZO+FqWjb0dnZY5co9fBi34GCaP/ts+QcLIWqSly73xIr04jsB9NOMqh+gqCYtQzOq1wAnLvemwjlSsvKZ9PlODsal8djQDswY0h4314pMaO98uUeO4NWpk7PDEEJUkjki7JfLPbcyfYhbAxM0ozoBSAd6X+5NRfWLS81h0mc7OZ2czSf39eamTs2cHVKFFWZmkW8243/LKGeHIoSoRmUmKM2otsG2nvxd2HrWtQF6l9fFXNQsu0+n8OCXu8krKGTx5D5c2652dCG/IO+oCXRdWlBC1DOXTFCaUd0GBGBbT/4O1aT9pRnVU5Kcag9d1/n8DzMRG0wEN/Ri+bS+tWppjAtyDx8BkAQlRC1iCI+cD2wANpkjwiq1CsYFZbWgkrDN5tAMaAL8RSlTD4maKbegkOlL9rDZlMhQY1PeHtedQN/aOTVU7uHDuDZpjHvTpuUfLISoKT4DhgNPGMIj87HNJLHRvqJuhZS5HpRmVAOA27E94muPrS/7sFKWgHcKWQ+qdLkFhfzri138fvwcL47qxKT+hho7vqkiTt4yGrcWwbT+5BNnhyJErVFT1oMCMIRHNgJuBkZgm1liL7ZktbKs8yqzYGFT4E5syaqVatKcPgunJKiLJWfl8/BXu9lpTubN27sxrrfT/2e6ItacHI726k2jB6fR9LHHnB2OELVGTUpQ/2QIj+wFDDdHhJU5MUOFe/GpJi0ReB943955QtQwp85lcf9nO0hIz2PunT0Y06P2r1aSd+wYWK14d+7s7FCEEFXEHBG2G9hd3nGXNVW1atJOX855wnHOpOZw74Id5BQUsmJaP65uHejskKpEzuHDgHSQEKI+qh1rKYgymc6m89CXu0nPKWDZtH50aVl3pgDKPXIE14YNcZMl2oWod2rHNALiktbsiWX0B3+QlV/IoinX1MrklL13L6lff130k3/67wb6hRkkanMnDyHqI0N45NPFvo/7x77XK3KNssZBvQmcVE3avH+UPw40V03aM5ULV1S1RX+c4uXvjtC/XSPeu+tqGvt5OjukStF1neSFC0l8+50S5S4BAbT54gs8Qg3k/XUcv0kTnRShEOIKTADetH+fjW2ZpguGA+VOrFlWC2oU9pVr/+FdIKyCAQoHsBRaeeW7w7z83RGGdW7G55OvqVHJSc/PJy0ykpyDh2xLY5R2TEEB8S+8QOLb7+A/cgTtNm2i/c9bCF2zGhcvL6IfmErm5s1QUICXdJAQojZSLvG9tO1SlfUOSldNmvWfhapJs2pGVZ63OEmepZBHl+7lpyMJTBkQyrMjjTVuwtf0H37gzCxb615xd8enTx+aPv00Xh2vAiD32DHOvvIqObt30+jhh2gyYwaKi60O7sHBtP5sIafvuZe4p2YB0kFCiFpKv8T30rZLVVaCytaMagfVpP1VvFAzqh2AnIrFJ6pSTn4h077cxW9/neOV0Z2Z2N/g7JBKlXPgIIqPDy1ef42cAwdJW7OGU2PHEnTvvehWKylLl+Li50eLt94k4JZbLjrfs107Wn36KdGTJoGLC+6tavdYLiHqqe6G8Mh0bK0lb/t37NteFbnAJQfqakZ1BLZxT//h7/7qvbE9S/y3atLWX0HgVaI+DdTNzLPwwOIodpxK5r9juzH+mpr7S9t8193g6oLhq68AsKSkkPS/uaSuWgWKQuCEO2k8YwZugWV3hc85fBhLYiINBg+ujrCFqFNq8kDdiipvqqMuwCygi73oMPCWatIOVkNs5aovCSotp4DJn+9kf2wac8Z3r1EDcPXCQqw5Obj6+dm2LRaO9r6GwDvvpNns8BLH5p04AS4ueIaGOiNUIeqVupCgyhwHpZq0Q4B0oXKilKx87vtsB0fPZvDh3VfXuKXZz308j5Rly2j/8xZcPDzIO3ESPTcXry5dLjrWs107J0QohKityupmvq6sE1WTVuZ68ppR9QJ+BTzt9/laNWkvaUY1CFgBGAAzMF41aSn2c2YDU4FCYKZq0n6ocE3qoOOJmTyyZDfm89nMv683g401bzbv9PXrKTx/nuyoKPwGDCD30CEAvLpIzzshxJUpqwV1LRADLAN2UMFugcXkAUNUk5apGVV34HfNqG4AxgKbVZMWoRnVcCAceEYzqp2w9ZvvDLQANmlG9SrVpBVW8r51wi/Hknjwy114u7uyaNI19G/f2NkhXSTfbCb/5EkAsn791ZagDh/Cxc8PjzYyXaMQtZl9cO3LgAr0MUeE7bKXGwANOGo/dLs5IuwhR8RQVoJqDtyEbfbyu4FIYJlq0g5X5MKqSdOBTPumu/1HB8YAg+zli4GtwDP28uWqScsDTmlG9TjQB/iz4tWpG46cSWf6kj2ENvZj0eRraOZfoQ4v1S7j560AeF51FZlbf6HZ7NnkHDpsm/nBpWZ1fRdCVNohbA2K0ta5OWGOCOvh6AAu+VtENWmFqknbqJq0iUA/4DiwVTOqMyp6cc2oumpGdR+QCPykmrQdQDPVpMXb7xEPXHhu1RJbi+2CWHtZCYqiTFMUZZeiKLssFktFQ6k19sekcveC7TTwcuPzSdWfnAozs0ic8z8KM8pfADNzyxY8O3ak4fjx5J8+Td7x4+SZTKW+fxJC1C7miDDNHBF2tPwjHafMP3M1o+qpGdWxwFfAdOA9YE1FL25Pcj2wrczbx94r8FJKe4R4URdDXdfn67reW9f13m5udWuu2+0nz3P3p7bktPLBa2keUP0tp4yNGzg/fz4pS5aUeVxhairZe/bgN3gQfoNuAOD8pwvQ8/PxlvdPQtQEbhf+mLf/TKvCa4cawiP3GsIjfzGER15XhdctoaxOEouxdS/fALxi79F3WVSTlqoZ1a3Y5l9K0IxqsGrS4jWjGoytdQW2FlPxwT0hwJnLvWdtsyc6hcmfRxES6M2XU/s6JTkBZG3bBkDKkqU0mjIFxaP0ZeIzf/sNCgtpMGQIHiEheLRrR9r33wNIC0qImsGi63rvsg4whEduwvY655+eM0eErb3EafFAa3NE2Hn7woPfGsIjO5sjwtIvcfxlK6sJch+QBVwFzNSM6oVyBds0SP5lXVgzqk2AAnty8gZuBP4LrMPWdT3C/nnhH2EdsFQzqnOwdZLoANSIpeUdbX9MKlMWRdHU35Ol/+pHkwbOmVdPt1rJ+nM77q1aURATQ/rGjQSMLr2zZsaWLbg2aVyUjPxuuIHkEydw8feXmR+EqCXMEWE3XsY5edg6wWGOCNttCI88gS1P7Kri8C6doFSTdqVvuYOBxZpRdcX2KHGlatK+14zqn8BKzahOBaKBcfb7HdaM6krgCGABpteHHny//3WOaV/uopGfB19O6VvtyUnXdSgoQPHwIFfTKExJoVn4M5z7ZD7Jixbjf8stFy11UZCQQNavv+E/ckRRZwi/G24g+bPP8O7SWZbGEKIOM4RHNgGSzRFhhYbwyLbYGhMnHXEvh73EUU3aAeDqUsrPA0Mvcc5rQJlr1NclW0wJPPTlHto28eWLKX1o6oTeeqkrV5E4Zw5tv1lD1h+2x3u+/ftjzcnh7MuvkLNnDz69ehUdb83KIubhhwEImvj3GG6fnlfj1rw5PtdeW70VEEI4hCE88jZs0901ASIN4ZH7zBFhw4DrgVcN4ZEWbGNWHzJHhCU7IoYypzqq6WrzVEdbjyYy7YvddGzegK+m9iXAx90pccQ89DCZW7fS4KabKMzIoDA5mbbr1mLNyeGvQYPx7dOHkPffA2zTGsU+OoPMX36h1byP8bv++hLXsublobi7SxdzIWqAujDVkfwmcQLbY73dtG/qx5dT+zgtOemFhWTv3o1LQAAZP/1E9s6d+PbvD4CLtzeB48eTsXkz+bGxAJyfP5/Mn3+m2XPPXpScAFw8PSU5CSGqjPw2qWZ7olN44Iso2jb2ZckDfWnoU3ovueqQd/Qo1owMmj09C4+2bcFqxXdA/6L9gffcDS4upHy1hLyTJzn30cf4jxxB0D33OC1mIUT9UbcGEtVwMcnZ/GvxLpr5e/HVA30J9K2+5JR75AiJ775L8Kuv4t6sGQDZu2ydbnz798ezXTvOL1yIzzXXFJ3j3rw5/sOGkfr11+Ts24fi40OzZ8tdpVkIIaqEtKCqyenzWUxZFEVBoZWFE6t3iXbdaiX+hRfJ+uVXEt6IKCrPjorCPSQE9+BgvHv0IOT993HxKtlRI2ji/VgzM8nZt49mT8/CrXHNmxNQCFE3SYKqBtHnsxnz4R+cTc9l3n29aN/Uz6H3S168mPgXXuBCB5jU1avJPXwY7969yNi4kczf/0C3WsmO2lWixVQa727d8Bs0CL+hQwkYO9ahcQshRHHyiM/BsvIsTPtyF7oO6x4dSGhjx3aq0XWd5MVfUHDmDL4DBuDdvTtJc/6Hd69etF64kFOjx3D2pZcImjyZwtRUfHqXOdAcgJCPPgRFkfFNQohqJd3MHSglK59Hluxhx6nzLJrch+uvauLwe+adPMnJkWEo7u64NmyI4u5OYXo6bZYswavjVWTv2UPcvx/HkmibYardTz/iITM/CFHn1IVu5tKCcpDsfAuTFkWhxafz1h3dqyU5AWT99hsALd78L3FPPIlrQACtFy/Cq+NVAPj07En7zZtI/+FHLElJuIeEVEtcQghRWZKgHMBSaOXRpXs5GJvKvHt7cXPn0uZidIzMX3/Do21b/EeMwMWvAR5tWuPRunWJYxR3dwJGhVVbTEIIcTkkQVUxXdd59puDbDEl8tptXao1OVlzcsiOiiLwrrsA8LtuYLXdWwghqpr04qtin/x6kpW7Ypk5pD339K3eZc+z/tyOnp+P73UOW55FCCGqjSSoKrT7dApv/XCUsK7BPH7TVdV+/5Rly3Bt3BifPmV3HRdCiNpAElQVScsuYOayvQQHePHG7V0d1iU75+AhClNT0XWdrD//pDDdtkZY3okTZP32G4F334XLJRYZFEKI2kTeQVWRF9cdIiE9l1UPXYu/l2Mmf807fhzz+PF4GAz4XnstKUuX4t6iBU1nh5OxYSOKhweBEyY45N5CCFHdJEFVgR8Pn2XtvjP8+8YOXN060GH3OTfvExQvLywJCaQsXYr/LbeQs3cvcTNmAtBw/HjcgoIcdn8hhKhOkqCukPlcFk+t2k+nYH8eGdTeYffJPXaM9PXraTRlMv4jR5KraQSMHYuenU3OwYOAbVoiIYSoK2QmiStgKbQy/pM/OZ6YSeTM62gV5OOQ+2Tv2UvsI48A0Dbye9waNXLIfYQQdUddmElCOklcgU9+Pcme6FT+79YuDktOAAn/+Q8uvr4YViyX5CSEqDcc9ohPM6qtgC+A5oAVmK+atHc1oxoErAAMgBkYr5q0FPs5s4Gp2Na5n6matB8cFd+V2nHyPP/76RhhXYMZ3b1Fpc7VrVYyf/4Zv8GDy12B1pqbS+7RozR64AE82lTvuCohhHAmR7agLMCTqklTgX7AdM2odgLCgc2qSesAbLZvY983AegMDAc+0oyqqwPju2yp2fk8vfoAIYHel9WlPHv7dmKnP0r2zqiS1129mtyjxwCwJCdzbv6n5Bw4AIWFeHfrWmXxCyFEbeCwFpRq0uKBePv3DM2oakBLYAwwyH7YYmAr8Iy9fLlq0vKAU5pRPQ70Af50VIyXI99iZeLnUcSn5fLllD6X1aU8z2wGoOBsPLquk/nzz6DrxD/3PP4jR9ByzhwS33yLtG+/xbtnTwC8ukqCEkLUL9XSi08zqgbgamAH0MyevFBNWrxmVJvaD2sJbC92Wqy9rARFUaYB0wA8nDAgdeHvp9gfk8pH9/Skb9vLex9UEBsHgCUpiZx9+4h9ZHrRvuy9+8g1mUhbuxaAnD17cGveHPemTUu9lhBC1FUO7yShGVU/YDXwb9WkpZdxaGnPyS7qYqjr+nxd13vrut7bza16e8nHJGfz7uZjDOvcjJFdgy/7OgUxMYA9Qe3eDYB7y5b43nA9lvh4zr7yKi4NGtBwwp0AeEvrSQhRDzk0QWlG1R1bclqimrQ19uIEzagG2/cHA4n28lig+Mp5IcAZR8ZXGbqu8+LaQ7goCi/d0rlS5xbExWE5d65oOz82FgBLYhLZu/fgERpK+82baPLoDABy9u6l4bg7CLp/IigK3ldfXXUVEUKIWsKRvfgUYCGgqSZtTrFd64CJQIT9c22x8qWaUZ0DtAA6ADsdFV9lrd4Tx89Hk3g+TKVFQ+9KnXt86I0o7u4YDx5A1/W/W1CJieSfOIHfTTcC4GXsiOLtjZ6bS+Bdd+EREkLomtV4tG1b5fURQoiazpHPyAYA9wEHNaO6z172LLbEtFIzqlOBaGAcgGrSDmtGdSVwBFsPwOmqSSt0YHwVlpqdT8QGE73bBDJ5QGilztUtFttnQQEA1rQ0rJmZgO39EoBPz16AbSHBBkOHori54WFf6dZLVaukDkIIUds4shff75T+Xglg6CXOeQ14zVExXa5XvztCanY+L0++BleXynUpL7A/zrsgbd13ALi3aU3B6WgAfPr0Kdrf8u23rjBaIYSoG2QmiXLsPp3Cmr1x/Ov6tnRpGVDp8/NOnSr6nrF5Mwmvv47v9dfR8LbbAHBr2hSPkIs6KwohRL0nCaoMuQWFPL5iHy0bevPIoHaXdY38k38nqNjpj+LRrh0hH3yAq33KIs+OHaskViGEqGtkNvMyzPvlBNHJ2Xw1tS8NLnONp3zzqRLbTWbOxMXDA/dg2/RIAbeMuuI4hRCiLpIW1CWcTMrkw5+PM7p7CwZ2aFyhc3RdJ2X5iqJVbsHWpdzFz69o2/fafgD4DRxA2/WRBIweXbWBCyFEHSEJ6hLmbvoLNxcXXhjVqcLn5B4+wtmXX+bM7GeLyixn4kuMY3L19y/67indx4UQ4pIkQZVi69FE1u0/w9SBoTRp4Fnh8wribFMYZf7yCynLl5P47rvknz6NZ4cOADQcN84h8QohRF0k76D+IT23gPDVB+nQ1I9Hh5S/Qm7ykiV4tm2L77XXkh992lZosXD25VeKjnFv0QLj4UNQyVnPhRCiPpME9Q9f/nmas+m5zLtvAF7u5a/2kfB//wFANWkUREeXeox7i2AU1xq5cogQQtRYkqCKybdY+eJPM9d1aEyPVg3LPd6an1/0Pe2770ld9TVenTrhc801WPPzSF22HLC1oIQQQlSOJKhi1h+MJyE9j4jbu13ymPT16/Fo1w6vjh0pPH++qDx58WIAAu+7j4a33QrYZivP3LRZEpQQotYxhEe+BdwC5AMngMnmiLBU+74Sq5+bI8Icsvq5dJIoZumOaNo29uWGDk1KlBemp6NbrQDEPfEkp8bcClBihvLcQ4do9txzRckJoOXbb2NYtapEzz0hhKglfgK6mCPCugHHgNkAhvDIi1Y/N4RHOuQdhiQou9iUbHaak7m9VwguxebbKzhzhmN9+pKyZGmJR3oAlqS/E5RnJ5XAuyaU2O/i5YV31y6ODVwIIRzAHBH2ozkizGLf3I5tCSSwr35ujgjLM0eEnQIurH5e5eQRn913++MBGN295OO4zF9+ASB7926sOTlF5fGvvIJXsWmKPEPbolTzAopCCFEGN0VRdhXbnq/r+vzLvNYUYIX9e4VWP68K8hvV7rv9Z+jRqiGtgnxKlOfbZxzP2LiRjI0bi8pTly3HNTCwaNutccVmmxBCiGpi0XW9d1kHGMIjNwHNS9n1nDkibK39mOewLYG0xL6vQqufVwVJUIDpbDpH4tNLnTUi//TpS55XmJJS9N3F19chsQkhhKOYI8JuLGu/ITxyIjAKGGqOCLuQhKpt9XN5BwWs2ROHu6vCbVf/3UrVrVaydu7EkpBw0fGG1V/j1rzkHx2Kt5fD4xRCiOpiCI8cDjwDjDZHhGUX27UOmGAIj/Q0hEeG4sDVz+t9gtJ1nfUH4xnQvjFBvh5F5Ulz5hB9/0Ryjxy56BzPDh2g0LbYr+JjeyTo2qBB9QQshBDV4wOgAfCTITxynyE8ch6AOSLsMHBh9fONwHRzRJhDVj+v94/4DsWlE5uSw8whHTi/YAF+Q4bgERpK1vYdlzzHxcOjqMNEi//8H1lRUQTcems1RSyEEI5njgi75Fxv5oiwaln9vN4nqPWH4nF1URjaviGJ971D4tvv4Nq4MR4hIaUe79GmDUBRgvLq1g3/kSOrLV4hhKgvHJagNKP6GbaXa4mqSetiLwvC1lXRAJiB8apJS7HvKzEyWTVpDhmZXJyu63x/4Az92zUi0MOFRHt54blzYE9EAL4DB+Lbry9BkyaBbntP2ODmm8jYsBH35qV1gBFCCHGlHPkOahG2UcbFhQObVZPWAdhs30YzqheNTNaMqsNnVz18Jp2Y5Bxu6d4C3WIpsS9n9+6i782enU2jBx5AcXNDcbetrNsiIoL2P2+RsU9CCOEgDktQqkn7FUj+R/EYYLH9+2Lg1mLly1WTlqeaNIeOTC5u5ylbeNd1aHxRgrqgwc034xEaelG5i6cn7sHBDo1PCCHqs+ruxddMNWnxAPbPpvbylkBMseMcNjK5uChzMiGB3vhs/ZEzs56+aL/fjUMJee9dFFnHSQghql1NeT5V4ZHJiqJMA6YBeHh4lHZIhei6TpQ5hes6NObMM9NLPabl229f9vWFEEJcmepuQSVoRjUYwP55oV9ChUcm67o+X9f13rqu93a7gvc/sSk5nMvMo2dzn0se4+Ilg2+FEMJZqjtBrQMm2r9PBNYWK5+gGVVPzag6dGTyBXtjUgFo+sHrF+1r+swzePfq5cjbCyGEKIcju5kvAwYBjTWjGgu8BEQAKzWjOhWIBsYBqCbtsGZUL4xMtgDTVZPmkJHJFxyOS8PD1YUWh6JKlDccP55GkyfRaPIkR95eCCFEORyWoFSTdtcldg29xPHVMjL5gtiUHFoGeuOmW0uU+/R1eOdBIYQQFVBv5+KLTc2huZ5zUXlAWJgTohFCCPFP9TdBJWfTXM91dhhCCCEuoaZ0M69W+RYr57PyCTxvLlHe8cB+5wQkhBDiIvWyBZWWUwCA99FDRWX+I0fgcgXjqoQQQlStepqg8gHwy/97Da4Wb73lrHCEEEKUol4mqNRsWwuqgT1BtfthI4qrw+emFUIIUQn1MkFdeMTnV2Drxefi5+fMcIQQQpSiXiaofItt7JOH1ZaocKmX/wxCCFGj1cvfzAVW2zy0blbbZBWuvr7ODEcIIUQp6meCsregXK22T0V67wkhRI1TLxPUuWXLAHDTC8G+Qq6oH1xdXenRo0fRT0REhLNDqjb9+/cHwGw2s3TpUidHI0T56uVA3cxDR6BHJ9yshbgFBjo7HFGNvL292bdvX5Ve02KxcCVLv1TX/bdt2wb8naDuvvtuR4cmxBWplwnKoti6lLvoVgwrVzg5mvrple8Oc+RMepVes1MLf166pfNlnWswGJg4cSLfffcdBQUFrFq1CqPRSFZWFjNmzODgwYNYLBZefvllxowZw6JFi4iMjCQ3N5esrCy+//57Jk2ahMlkQlVVzGYzH374Ifv37+fQoUP873//A+DTTz9F0zTmzJlT4v5+fn48+OCD/PzzzwQGBrJ8+XKaNGnCiRMnmD59OklJSfj4+PDpp59iNBqZNGkSQUFB7N27l549e/LOO+8UXevw4cNMnjyZ/Px8rFYrq1evpkOHDvj5+ZGZmUl4eDiaptGjRw8mTpzIzJkzCQ8PZ+vWreTl5TF9+nQefPDBy/8fQogqUi8f8RXae+256YW4N2/u5GhEdcrJySnxiG/Fir//QGncuDF79uzh4Ycf5m37asqvvfYaQ4YMISoqip9//plZs2aRlZUFwJ9//snixYvZsmULH330EYGBgRw4cIAXXniB3bt3AzBhwgTWrVtHQYGtx+jnn3/O5MmTL4orKyuLnj17smfPHm644QZeeeUVAKZNm8b777/P7t27efvtt3nkkUeKzjl27BibNm0qkZwA5s2bx2OPPca+ffvYtWsXISEhJfZHRERw3XXXsW/fPh5//HEWLlxIQEAAUVFRREVF8emnn3Lq1Kkr/acW4orV6xaUm9VazpHCUS63pXOlynrEN3bsWAB69erFmjVrAPjxxx9Zt25dUcLKzc0lOjoagJtuuomgoCAAfv/9dx577DEAunTpQrdu3QDw9fVlyJAhfP/996iqSkFBAV27dr3o3i4uLtx5550A3HvvvYwdO5bMzEy2bdvGuHHjio7Ly8sr+j5u3DhcSxlgfu211/Laa68RGxvL2LFj6dChQ5n/Jj/++CMHDhzg66+/BiAtLY2//vqL0NDQMs8TwtHqZ4Jysf2f2tXq0DURRS3j6ekJ2DpSWCwWAHRdZ/Xq1XTs2LHEsTt27MC32PAEXdcved0HHniA119/HaPRWGrrqTSKomC1WmnYsOElE6rvJYZH3H333fTt25fIyEiGDRvGggULGDJkyCXvpes677//PsOGDatQbEJUl3r6iM/egnKvl/lZVMKwYcN4//33ixLQ3r17Sz1u4MCBrFy5EoAjR45w8ODBon19+/YlJiaGpUuXctddpa/jabVai1owS5cuZeDAgfj7+xMaGsqqVasAWyLZv7/8GfdPnjxJ27ZtmTlzJqNHj+bAgQMl9jdo0ICMjIwSdfz444+LHkMeO3as6DGmEM5UL39Du4S0wkW30vbrlc4ORVSzC++gLhg+fHiZXc1feOEF/v3vf9OtWzd0XcdgMPD9999fdNwjjzzCxIkT6datG1dffTXdunUjICCgaP/48ePZt28fgZfoNerr68vhw4fp1asXAQEBRe/GlixZwsMPP8x//vMfCgoKmDBhAt27dy+zjitWrOCrr77C3d2d5s2b8+KLL5bY361bN9zc3OjevTuTJk3isccew2w207NnT3Rdp0mTJnz77bdl3kOI6qCU9WiipvP19dUv5y+9QqtOQaEVL3eZIFZUjcLCQgoKCvDy8uLEiRMMHTqUY8eO4WEfBD5q1Cgef/xxhg4dWur5F3rYCVFVFEXJ1nW9Vk+TUy9bUK4uCq4ukpxE1cnOzmbw4MEUFBSg6zoff/wxHh4epKam0qdPH7p3737J5CSEKF2Na0FpRnU48C7gCixQTdoln79cbgtKCCHqurrQgqpRnSQ0o+oKfAiMADoBd2lGtZNzoxJCCOEMNSpBAX2A46pJO6matHxgOTDGyTEJIYRwgpqWoFoCMcW2Y+1lRRRFmaYoyi5FUXZdGKsihBCi7qlpCUoppazESzJd1+frut5b1/XezpygUwghhGPVtAQVC7Qqth0CnHFSLEIIIZyopjVBooAOmlENBeKACYCsCSCEEPVQjWpBqSbNAjwK/ABowErVpB12blRCCCGcocaNg6oMRVGsQM5lnu4G1KdeFlLfuqs+1RWkvhXlret6jWqEVFatTlBXQlGUXbqu93Z2HNVF6lt31ae6gtS3PqnV2VUIIUTdJQlKCCFEjVSfE9R8ZwdQzaS+dVd9qitIfeuNevsOSgghRM1Wn1tQQgghajBJUEIIIWqkepmgFEUZrijKUUVRjiuKEu7seK6UoiifKYqSqCjKoWJlQYqi/KQoyl/2z8Bi+2bb635UUZRhzon68imK0kpRlJ8VRdEURTmsKMpj9vI6V2dFUbwURdmpKMp+e11fsZfXuboWpyiKq6IoexVF+d6+XWfrqyiKWVGUg4qi7FMUZZe9rM7Wt1J0Xa9XP9gWQjwBtAU8gP1AJ2fHdYV1uh7oCRwqVvYmEG7/Hg781/69k73OnkCo/d/C1dl1qGR9g4Ge9u8NgGP2etW5OmObQNnP/t0d2AH0q4t1/Ue9nwCWAt/bt+tsfQEz0PgfZXW2vpX5qY8tqD7AcV3XT+q6XifWnNJ1/Vcg+R/FY4DF9u+LgVuLlS/XdT1P1/VTwHFs/ya1hq7r8bqu77F/z8A2LVZL6mCddZtM+6a7/UenDtb1AkVRQoAwYEGx4jpb30uob/UtVX1MUOWuOVVHNNN1PR5sv9CBpvbyOlV/RVEMwNXYWhZ1ss72x137gETgJ13X62xd7eYCTwPWYmV1ub468KOiKLsVRZlmL6vL9a2wmjabeXUod82pOq7O1F9RFD9gNfBvXdfTFaW0qtkOLaWs1tRZ1/VCoIeiKA2BbxRF6VLG4bW6roqijAISdV3frSjKoIqcUkpZramv3QBd188oitIU+ElRFFMZx9aF+lZYfWxB1Zc1pxIURQkGsH8m2svrRP0VRXHHlpyW6Lq+xl5cp+us63oqsBUYTt2t6wBgtKIoZmyP34coivIVdbe+6Lp+xv6ZCHyD7ZFdna1vZdTHBBUFdFAUJVRRFA9sa06tc3JMjrAOmGj/PhFYW6x8gqIonoqihAIdgJ1OiO+yKbam0kJA03V9TrFdda7OiqI0sbecUBTFG7gRMFEH6wqg6/psXddDdF03YPv/5hZd1++ljtZXURRfRVEaXPgO3Awcoo7Wt9Kc3UvDGT/ASGw9v04Azzk7niqozzIgHijA9hfWVKARsBn4y/4ZVOz45+x1PwqMcHb8l1HfgdgeaxwA9tl/RtbFOgPdgL32uh4CXrSX17m6llL3Qfzdi69O1hdbb+L99p/DF34f1dX6VvZHpjoSQghRI9XHR3xCCCFqAUlQQgghaiRJUEIIIWokSVBCCCFqJElQQgghaiRJUKJeURSlkX3W6H2KopxVFCWu2LbHP46drSjKPf8om6Qoygf27y6KoixWbLPJlzqNhaIogy7MyF1GTA8pinL/ldZNiLqmPk51JOoxXdfPAz0AFEV5GcjUdf3tSxx+MzC+tB32hDQP2+Stk/UrGK+h6/q8yz1XiLpMWlBClEJRFH/AQ9f1pEsc8i62wZT367putc8I8JmiKFH2dYxKzJBvb22ZL8wKYS87rihKM0VRXlYU5Sl72VZFUf5rXwPqmKIo1zmoikLUeJKghCjdjdhG8JfmbqAXMEHXdYu97Dls0/JcAwwG3rJPXQOArutWbNPV3AagKEpfwKzrekIp13fTdb0P8G/gpSqoixC1kiQoIUo3HNhwiX17gDaUXIfnZiDcvizGVsALaP2P81YAd9q/T7Bvl+bC5Le7AUMlYhaiTpEEJUTp+nDpSThN2N5NrVAUpbO9TAFu13W9h/2nta7r2j/O+xNoryhKE2wL0K2hdHn2z0LkPbGoxyRBCfEP9qRj0m3rMJVK1/VtwENApKIorYEfgBkXevMpinJ1Kefo2JZTmINtJvbzjohfiLpC/joT4mIjgI3lHaTr+vf21tBG4CbgReCAPUmZgVGlnLYC25Ivk6oqWCHqKpnNXIh/UBTlJ2y98+KdHYsQ9ZkkKCGEEDWSvIMSQghRI0mCEkIIUSNJghJCCFEjSYISQghRI0mCEkIIUSNJghJCCFEj/T8CeH4fNNDAnAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax1 = plt.subplots()\n",
    "\n",
    "color = 'tab:red'\n",
    "ax1.set_xlabel('T / Kelvin')\n",
    "ax1.set_ylabel('MC Accepts', color=color)\n",
    "ax1.plot(T_plot_1,MCA_plot_1, color=color, label=\"MC Accepts\")\n",
    "ax1.tick_params(axis='y', labelcolor=color)\n",
    "ax1.legend(loc=\"best\",frameon=False)\n",
    "ax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis\n",
    "\n",
    "color = 'tab:blue'\n",
    "ax2.set_ylabel('E / meV', color=color)  # we already handled the x-label with ax1\n",
    "ax2.plot(T_plot_1,E_plot_1, color=color, label=\"Energy per site\") \n",
    "ax2.tick_params(axis='y', labelcolor=color)\n",
    "ax2.legend(loc=\"lower center\",frameon=False)\n",
    "\n",
    "fig.tight_layout()  # otherwise the right y-label is slightly clipped\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "31f70b51",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAYT0lEQVR4nO3df5AcZ33n8fdnV16RFTi5rJbAWfauIL5yyTkw9sYxgSQYDLFFKsYFCTJjsFEqKkkxGDjqkEspApXoDxJC4SRItgqU6NAeJnUxF8XYOMQXIBUSohWxjWVbIIy2LOyLpeXOxsGxvNI3f0yvMhrNzPbMTvfszvN5VU3t9NPdM99nbPVn+sc8rYjAzMzSNdDrAszMrLccBGZmiXMQmJklzkFgZpY4B4GZWeKW9bqAdq1cuTLGx8d7XYaZ2ZKyf//+YxEx2mjekguC8fFxpqamel2GmdmSImm62TwfGjIzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS1wSQTA5Ocn4+DgDAwOMj48zOTnZ65LMzBaNJXf5aLsmJydZv349x48fB2B6epr169cDUKlUelmamdmioKU2DPXExES08zuClStXMjMzc0b7yMgIx44d62ZpZmaLlqT9ETHRaF7fHxpqFAKt2s3MUtP3QWBmZq05CMzMEucgMDNLnIPAzCxxDgIzs8Q5CMzMEucgMDNLnIPAzCxxDgIzs8Q5CMzMEucgMDNLXKFBIOlKSQclHZK0pcVyPyvphKS3FVmPmZmdqbAgkDQIfAq4ClgDXCtpTZPlPgbcU1QtZmbWXJF7BJcChyLi0Yg4DtwOXN1gufcAfwE8WWAtZmbWRJFBcA7wWM30kaztFEnnANcAt7Z6IUkbJE1Jmjp69GjXCzUzS1mRQaAGbfV3wfkk8KGIONHqhSJiZ0RMRMTE6Ohot+ozMzOKvVXlEeDcmulVwON1y0wAt0sCWAmslTQbEf+7wLrMzKxGkUGwDzhf0mrg+8A64B21C0TE6rnnkv4MuNMhYGZWrsKCICJmJd1I9WqgQWBXRByQtDGb3/K8gJmZlaPIPQIi4i7grrq2hgEQETcUWYuZmTXmXxabmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmlrikg2BycrLXJZiZ9VzSQbB169Zel2Bm1nN9HwQjIyNN501PT5dYiZnZ4tT3QXDLLbf0ugQzs0Wt74OgUqn0ugQzs0Wt74PAzMxacxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSWu0CCQdKWkg5IOSdrSYP7Vkh6QdJ+kKUmvLbIeMzM707KiXljSIPAp4I3AEWCfpL0R8VDNYvcCeyMiJL0C+HPggqJqMjOzMxW5R3ApcCgiHo2I48DtwNW1C0TEMxER2eQKIDAzs1IVGQTnAI/VTB/J2k4j6RpJjwBfBNY3eiFJG7JDR1NHjx4tpFgzs1QVGQRq0HbGN/6I+EJEXAC8BfjdRi8UETsjYiIiJkZHR7tbpZlZ4ooMgiPAuTXTq4DHmy0cEV8DXi5pZYE1mZlZnSKDYB9wvqTVkoaAdcDe2gUk/bQkZc8vBoaAmQJrMjOzOoVdNRQRs5JuBO4BBoFdEXFA0sZs/q3AW4F3SXoeeBZ4e83JYzMzK4GW2nZ3YmIipqam2lon2+loaKn138ysE5L2R8REo3m59ggkvQb4CDCWrSMgIuJl3SrSzMx6I++hoc8A7wf2AyeKK8fMzMqWNwieioi7C63EzMx6Im8Q/K2kPwDuAJ6ba4yIbxZSlZmZlSZvEPxc9rf2REMAr+9uOWZmVrZcQRARlxddiJmZ9UauH5RJ+nFJn5gb70fSH0r68aKLMzOz4uX9ZfEu4IfAr2ePp4E/LaooMzMrT95zBC+PiLfWTH9U0n0F1GNmZiXLu0fwbO3dw7IfmD1bTElmZlamvHsEm4Dd2XkBAT8AbiiqKDMzK0/eq4buA14p6exs+ukiizIzs/K0DAJJ10XEHkkfqGsHICI+UWBtZmZWgvn2CFZkf1/UYJ6H7TQz6wMtgyAibsue/k1E/H3tvOyEsZmZLXF5rxr645xtZma2xMx3juDVwM8Do3XnCc6metcxMzNb4uY7RzAEvDBbrvY8wdPA24oqyszMyjPfOYKvAl+V9GcRMV1STWZmVqK8Pyj7UXY/gguBF8w1RoSHoTYzW+LyniyeBB4BVgMfBQ4D+wqqyczMSpQ3CEYi4jPA8xHx1YhYD1xWYF1mZlaSvIeGns/+PiHpzcDjwKpiSjIzszLlDYLfywac+29Ufz9wNvD+wqoyM7PS5A2C+yPiKeAp4HIASS8prCozMytN3nME35P0OUnDNW13FVGQmZmVK28QfAv4O+DvJL08a1MxJZmZWZnyHhqKiNgu6X7gryR9CI8+ambWF/IGgQAi4u8lvQH4PHBBYVWZmVlp8gbB2rknEfGEpNdTHYzOzMyWuFx3KAOunbsrWZ2vFVKVmZmVxncoMzNLnO9QZmaWON+hzMwscb5DmZlZ4nyHMjOzxLV1hzJJKyLiX0uqzczMSpD3HMF/lvQQ8DCApFdK2j7fSpKulHRQ0iFJWxrMr0h6IHt8XdIr2yvfzMwWKm8QfBL4ZWAGICLuB36x1QqSBoFPAVcBa6j+FmFN3WLfA34pIl4B/C6wM3flZmbWFXmDgIh4rK7pxDyrXAociohHI+I4cDtwdd1rfj0i/l82+Y/4ZjdmZqXLGwSPSfp5ICQNSfog2WGiFs4BasPjSNbWzG8AdzeaIWmDpClJU0ePHs1ZspmZ5ZE3CDYCv0V1Q34EuCibbqXRmBQNf40s6XKqQfChRvMjYmdETETExOjoaM6Szcwsj1yDzkXEMaDS5msfAc6tmV5F9V7Hp5H0CuDTwFURMdPme5iZ2QLlCgJJo8BvAuO160TE+har7QPOl7Qa+D6wDnhH3eueB9wBvDMivt1W5WZm1hV5h6H+S6p3KPsb5j9JDEBEzEq6EbiH6q+Qd0XEAUkbs/m3Ah8GRoDt2eimsxEx0V4XzMxsIRQx/yCiku6LiIuKL2d+ExMTMTU11dY6TYbQBiBP/83MljpJ+5t90c57svhOSWvnX8zMzJaavEFwE9UweFbS05J+KOnpIgszM7Ny5L1qqNGNaczMrA/kvWro4gbNTwHTETHb3ZLMzKxMea8a2g5cDHwrm/6vwP3AiKSNEfHXRRRnZmbFy3uO4DDwqoi4JCIuofrL4geBK4DfL6Y0MzMrQ94guCAiDsxNRMRDVIPh0WLKMjOzsuQ9NHRQ0g6qI4gCvB34tqTlwPOFVGZmZqXIu0dwA3AIeB/wfuDRrO154PIC6jIzs5LkCoKIeDYi/jAiromIt0TExyPiRxFxMiKeKbrIIk1OTva6BDOznsoVBJLOl/S/JD0k6dG5R9HFlWHr1q29LsHMrKfyHhr6U2AHMEv1UND/AD5bVFHdNjg42HTe9PR0iZWYmS0+eYPgxyLiXqqD1E1HxEeA1xdXVndt2LCh6byBgdx36zQz60t5t4L/JmkA+I6kGyVdA7y4wLq6avv27U3nnTx5ssRKzMwWn7xB8D5gGHgvcAnwTuD6gmoyM7MS5R10bl/29Bng3cWVY2ZmZWsZBJL2tpofEb/a3XLMzKxs8+0RvBp4DPgc8A2g+a2+zMxsSZovCF4CvBG4luqN578IfK523CEzM1vaWp4sjogTEfGliLgeuIzqMBNfkfSeUqozM7PCzXuyOBtY7s1U9wrGgT8C7ii2LDMzK8t8J4t3Az8D3A18NCIeLKUqMzMrzXx7BO8E/hX4L8B7pVPnigVERJxdYG1mZlaClkEQER5/wcysz3lDb2aWOAeBmVniHARmZolzEJiZJc5BgG9XaWZpcxAAN910U69LMDPrGQcBMDMz0+sSzMx6JpkgWLFiRa9LMDNblJIJgttuu63XJZiZLUrJBEGlUul1CWZmi1IyQWBmZo05CMzMEucgMDNLXKFBIOlKSQclHZK0pcH8CyT9g6TnJH2wyFrm4x+VmVmqCgsCSYPAp4CrgDXAtZLW1C32A+C9wMeLqiMv/6jMzFJV5B7BpcChiHg0Io4DtwNX1y4QEU9GxD7g+QLrOKXmxjpn8I/KzCxVRQbBOcBjNdNHsra2SdogaUrS1NGjRzsuaOPGjR2va2bWr4oMgkZfv6OTF4qInRExERETo6OjHRe0ffv2jtc1M+tXRQbBEeDcmulVwOMFvp+ZmXWgyCDYB5wvabWkIWAdsLfA9zMzsw4UFgQRMQvcCNwDPAz8eUQckLRR0kYASS+RdAT4APDbko5IOruomubjS0jNLEWK6Oiwfc9MTEzE1NRUx+u3unJoZGSEY8eOdfzaZmaLlaT9ETHRaF5yvyxuNRy1LyE1sxQlFwQejtrM7HTJBYGHozYzO11yQTCfzZs397oEM7NSOQjq7Nixo9clmJmVKskgaHXlEPgyUjNLS5JBMN+YQzfccEM5hZiZLQJJBsF8Yw7Nzs56r8DMkpFkEAAMDQ21nO+9AjNLRbJBsGvXrpbzZ2dnfQWRmSUhuSEmai1btowTJ040nS+JkydPduW9zMx6yUNMNLF79+6W8yPC5wrMrO8lvUcAsHz5co4fP950/uDgILOzs117PzOzXvAeQQvznSs4ceIEF154YUnVmJmVL/kgqFQqLUckBXjooYe44oorSqrIzKxcyQcB5BuR9N5770WSryQys77jIKC6V7BmzZpcy+7YscOBYGZ9xUGQOXDgwLxjENXasWMHy5cv91VFZrbkOQhqfPazn21r+ePHj3PdddchCUkMDAx4T8HMlhwHQY1KpcKmTZs6Xj8i2LFjB8uWLTsVDitXrvReg5ktag6COtu3b19QGACn/Vp5ZmbmtL0GB4OZLTYOgga6EQbNOBjMbLFxEDSxfft29uzZw9jYWKHvMxcMK1euZPPmzYyPjzMwMMD4+LgDwsxKkfwQE3lNTk5y/fXXtxykrigjIyPccsstVCqV0t/bzPqDh5jogkqlwu7duxkZGSn9vesPJ9U+BgYGkMTg4KAPN5lZRxwEbahUKhw7doyIOO3xhje8oWc1ze3R1Q6XXXu4aS4QJicnfdjJzBryoaEu2bx5M7fddtupDfLQ0BCSeO6553pcWXOSiAjGxsZYu3Ytd911F9PT0wwODnLixIlTez8zMzMt28bGxti2bZsPXZktYq0ODZ3x7XaxPy655JJYivbs2RMjIyMB9OVDUgAxNjYWe/bsiT179sTY2FhIOtW2kM+uW69llipgKppsV3u+YW/3sVSDoJFNmzb1fANedlCMjIzEyMhISMr9/KyzzmoZOotNo+DKG2bthF4779Nuu/UfB8Ei1mhPYcWKFbF8+fKeb7yXwmNoaKhhgIyNjcWmTZva2sjl2bDWv2aj6eHh4TNqrA+z4eHhXOs2Wm6urnbep9nrNmqfL7gcKkuTg2CJmvuHBcTAwEDLb9rN5vtx5kau/vNttufRaMM632Puv0cnyzZbt759eHi4rcOMg4ODbbWPjIy0DIh2Q8UWBxwE6ej3cxELfQwODjbd8PvR+jE2Nnbqi0mjz7XZOnP/X3ayt9DJnok1hoMgPbV7E3P/SGsPZcD8exlzh1pqX6NRWzvfgv0o/9HuHkGzh6S2/1tLaroXkef8Rbt7Ju0ERWpBgoPAmunGP4ZmeyFzG42FBsXceYBuvFYZj/oaGx1iynsYqNlyjQ7fdOscQbM9yk73CJqtM/f/W7MNeqv1Ws2b+39yvvBJLUhwEFgZ8pxEbPeqoVaHA4o+vNPJOYJWJ3c7OTHc7gndblw11Mk38VZ1tgq9Vhv0Vuu1mhcR8wZFGUHS6nPvBQeB9a1mIVN7GExSrsMgtVcgtdqAt5pu5x96EZeUdks3rxrqdGO/kD2C+YKi6CCZ+zwWU1j0LAiAK4GDwCFgS4P5Av4om/8AcPF8r+kgsE40u9yy2Z6HdU+nh38Wco5goRvyhQZJnveY77Op/wwXGhY9CQJgEPgu8DJgCLgfWFO3zFrg7iwQLgO+Md/rOgisU4tpNz01nZwQbrVennkLObTTjT2CssNiPr0KglcD99RM3wzcXLfMbcC1NdMHgZe2el0HgVl/KSqgF3KytxvnCMoMizx6FQRvAz5dM/1O4E/qlrkTeG3N9L3ARIPX2gBMAVPnnXdeW503M+vEQq8aKjMs8uhVEPxagyD447plvtggCC5p9breIzCzpaKssMijVRAUeT+CI8C5NdOrgMc7WMbMbEmqVCocPnyYkydPcvjw4TOGaq9UKuzcuZOxsTEkMTY2xs6dO09bbtu2bQwPD5+23vDwMNu2betanUUGwT7gfEmrJQ0B64C9dcvsBd6lqsuApyLiiQJrMjNbVLoRFgu1rGuvVCciZiXdCNxD9QqiXRFxQNLGbP6twF1Urxw6BPwIeHdR9ZiZLVWVSqXQGz8VFgQAEXEX1Y19bdutNc8D+K0iazAzs9Z8z2Izs8Q5CMzMEucgMDNLnIPAzCxxqp6vXTokHQWmO1x9JXCsi+UsBe5zGtznNCykz2MRMdpoxpILgoWQNBURE72uo0zucxrc5zQU1WcfGjIzS5yDwMwscakFwc5eF9AD7nMa3Oc0FNLnpM4RmJnZmVLbIzAzszoOAjOzxCUTBJKulHRQ0iFJW3pdT7dI2iXpSUkP1rT9pKQvS/pO9vc/1cy7OfsMDkr65d5U3TlJ50r6W0kPSzog6aasvZ/7/AJJ/yTp/qzPH83a+7bPcyQNSvpnSXdm033dZ0mHJX1L0n2SprK24vvc7I41/fSgOgz2d4GXAUPA/cCaXtfVpb79InAx8GBN2+8DW7LnW4CPZc/XZH1fDqzOPpPBXvehzf6+FLg4e/4i4NtZv/q5zwJemD0/C/gGcFk/97mm7x8A/idwZzbd130GDgMr69oK73MqewSXAoci4tGIOA7cDlzd45q6IiK+BvygrvlqYHf2fDfwlpr22yPiuYj4HtX7QFxaRp3dEhFPRMQ3s+c/BB4GzqG/+xwR8Uw2eVb2CPq4zwCSVgFvBj5d09zXfW6i8D6nEgTnAI/VTB/J2vrVT0V2p7fs74uz9r76HCSNA6+i+g25r/ucHSK5D3gS+HJE9H2fgU8C/x04WdPW730O4K8l7Ze0IWsrvM+F3phmEVGDthSvm+2bz0HSC4G/AN4XEU9LjbpWXbRB25Lrc0ScAC6S9BPAFyT9TIvFl3yfJf0K8GRE7Jf0ujyrNGhbUn3OvCYiHpf0YuDLkh5psWzX+pzKHsER4Nya6VXA4z2qpQz/IumlANnfJ7P2vvgcJJ1FNQQmI+KOrLmv+zwnIv4/8BXgSvq7z68BflXSYaqHcl8vaQ/93Wci4vHs75PAF6ge6im8z6kEwT7gfEmrJQ0B64C9Pa6pSHuB67Pn1wN/WdO+TtJySauB84F/6kF9HVP1q/9ngIcj4hM1s/q5z6PZngCSfgy4AniEPu5zRNwcEasiYpzqv9f/ExHX0cd9lrRC0ovmngNvAh6kjD73+ix5iWfj11K9wuS7wNZe19PFfn0OeAJ4nuo3hN8ARoB7ge9kf3+yZvmt2WdwELiq1/V30N/XUt39fQC4L3us7fM+vwL456zPDwIfztr7ts91/X8d/3HVUN/2mepVjfdnjwNz26ky+uwhJszMEpfKoSEzM2vCQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWDJkTSSje54n6T/K+n7NdNDdcveLKlS13aDpD/Jng9I2q3qKLANf94s6XVzo2e2qGmjpHcttG9mnUhliAmzUyJiBrgIQNJHgGci4uNNFn8T8OuNZmQb/lupDgL37ljAtdgRcWun65otlPcIzJqQdDYwFBFHmyxyC9Uf+7wrIk5mvwzdJWlfNob+aSPcZnsPh+d+JZy1HZL0U5I+IumDWdtXJH0suwfBtyX9QkFdNAMcBGatXEH1l5yNvAO4BFgXEbNZ21aqQyH8LHA58AfZUAEARMRJqsMDXAMg6eeAwxHxLw1ef1lEXAq8D/idLvTFrCkHgVlzVwJ3N5n3TWCM08d/fxOwJRsu+ivAC4Dz6tb7PPD27Pm6bLqRucH09gPjbdRs1jYHgVlzl9J8EK9HqJ47+LykC7M2AW+NiIuyx3kR8XDdev8A/LSkUao3GLmDxp7L/p7A5/KsYA4CswayjfsjUb0PQEMR8XVgI/BFSecB9wDvmbt6SNKrGqwTVIcX/gTVEVRniqjfrB3+pmHW2FXAl+ZbKCLuzL7dfwl4I/Bh4IEsDA4Dv9Jgtc9THRr9hm4Va7YQHn3UrAFJX6Z6NdATva7FrGgOAjOzxPkcgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJ+3dJXq4PUWw3vAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.close('all')\n",
    "plt.figure()\n",
    "plt.plot(T_plot_1,Mag_plot_1,'ok')\n",
    "#plt.gca().invert_xaxis()\n",
    "plt.xlabel('T / Kelvin')\n",
    "plt.ylabel('Magnetization')\n",
    "#plt.xlim([1,0])\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  },
  "vscode": {
   "interpreter": {
    "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}