Monte Carlo v2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b3efdd2a",
   "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": "markdown",
   "id": "5aec40fb",
   "metadata": {},
   "source": [
    "## Define Sites"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "864f29b4",
   "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",
    "        g = 2.0\n",
    "        muB = 5e-2 #mev/T\n",
    "        return -g*muB*self.getspin()*B\n",
    "\n",
    "    def getdeltaexchangeenergy(self,J):\n",
    "        E = 0\n",
    "        NN = self.getnearestneighbourscubic()\n",
    "        for nn in NN:\n",
    "            E += -J*self.spin*nn.spin\n",
    "        return E\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]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2df6ff72",
   "metadata": {},
   "source": [
    "## Define Lattice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2e7363fb",
   "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 populatesites(self):\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 gettotalenergy1(self,J):\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",
    "                    E += -J*self.sites[self.getlistposition(i, j, k)].getspin() * self.sites[self.getlistposition(i, j, (k+1)% self.N3)].getspin()\n",
    "        \n",
    "        for i in range(self.N1):\n",
    "            for k in range(self.N3):\n",
    "                for j in range(self.N2):\n",
    "                    E += -J*self.sites[self.getlistposition(i, j, k)].getspin() * self.sites[self.getlistposition(i, (j+1)% self.N2, k)].getspin()\n",
    "                \n",
    "        for j in range(self.N2):\n",
    "            for k in range(self.N3):\n",
    "                for i in range(self.N1):\n",
    "                    E += -J*self.sites[self.getlistposition(i, j, k)].getspin() * self.sites[self.getlistposition((i+1)% self.N1, j, k)].getspin()\n",
    "                \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 plotDvM(self,MCS,data,dataname):\n",
    "        MCS_plot = np.arange(0,MCS)\n",
    "        plt.close('all')\n",
    "        plt.figure()\n",
    "        plt.plot(MCS_plot,data,'ok')\n",
    "        plt.xlabel('MC Accepts')\n",
    "        plt.ylabel(dataname)\n",
    "        plt.show()\n",
    "    \n",
    "    def plotEvT(self,E,T):\n",
    "        plt.close('all')\n",
    "        plt.figure()\n",
    "        plt.plot(T,E,'ok')\n",
    "        plt.xlabel('Temperature')\n",
    "        plt.ylabel('Energy')\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "53cb8508",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "200.00  64\n",
      "184.00  59\n",
      "169.28  51\n",
      "155.74  47\n",
      "143.28  48\n",
      "131.82  46\n",
      "121.27  39\n",
      "111.57  30\n",
      "102.64  28\n",
      "94.43  37\n",
      "86.88  28\n",
      "79.93  33\n",
      "73.53  23\n",
      "67.65  30\n",
      "62.24  26\n",
      "57.26  17\n",
      "52.68  17\n",
      "48.46  25\n",
      "44.59  30\n",
      "41.02  25\n",
      "37.74  28\n",
      "34.72  13\n",
      "31.94  21\n",
      "29.39  18\n",
      "27.04  16\n",
      "24.87  18\n",
      "22.88  18\n",
      "21.05  18\n",
      "19.37  12\n",
      "17.82  16\n",
      "16.39  17\n",
      "15.08  18\n",
      "13.88  15\n",
      "12.77  11\n",
      "11.74  14\n",
      "10.80  19\n",
      "9.94  18\n",
      "9.14  16\n",
      "8.41  14\n",
      "7.74  16\n",
      "7.12  9\n",
      "6.55  19\n",
      "6.03  10\n",
      "5.55  13\n",
      "5.10  14\n",
      "4.69  12\n",
      "4.32  12\n",
      "3.97  14\n",
      "3.65  12\n",
      "3.36  12\n",
      "3.09  9\n",
      "2.85  12\n",
      "2.62  6\n",
      "2.41  8\n",
      "2.22  15\n",
      "2.04  20\n",
      "1.88  15\n",
      "1.73  15\n",
      "1.59  17\n",
      "1.46  16\n",
      "1.34  13\n",
      "1.24  13\n",
      "1.14  13\n",
      "1.05  17\n",
      "0.96  10\n",
      "0.89  10\n",
      "0.81  10\n",
      "0.75  7\n",
      "0.69  10\n",
      "0.63  14\n",
      "0.58  12\n",
      "0.54  8\n",
      "0.49  5\n",
      "0.45  10\n",
      "0.42  12\n",
      "0.38  4\n",
      "0.35  11\n",
      "0.33  4\n",
      "0.30  5\n",
      "0.28  10\n",
      "0.25  14\n",
      "0.23  12\n",
      "0.21  5\n",
      "0.20  10\n",
      "0.18  10\n",
      "0.17  6\n",
      "0.15  13\n",
      "0.14  4\n",
      "0.13  6\n",
      "0.12  19\n",
      "0.11  7\n",
      "0.10  10\n",
      "0.09  5\n",
      "0.09  11\n",
      "0.08  4\n",
      "0.07  9\n",
      "0.07  5\n",
      "0.06  3\n",
      "0.06  3\n",
      "0.05  2\n",
      "Runtime:  0.062050\n"
     ]
    }
   ],
   "source": [
    "starttime = time()\n",
    "\n",
    "MCS = 100 # MC steps per temperature\n",
    "J = 100 # meV\n",
    "T = 200 # Kelvin\n",
    "T_p = 0.92\n",
    "kB = 8.617e-2 # meV/K\n",
    "\n",
    "MCA_plot = []\n",
    "T_plot = []\n",
    "E_plot = []\n",
    "\n",
    "# set up the sites\n",
    "lat = lattice(10,10,10) # set up the size of lattice\n",
    "lat.populatesites()     # Assign each site a coordinate\n",
    "lat.initialisespins()   # Assign each site a random spin\n",
    "E_0 = lat.gettotalenergy1(J) # get initial total energy of the system\n",
    "E_t = E_0\n",
    "\n",
    "while (T>0.05):\n",
    "    \n",
    "    MCA = 0    # Initialise MC accepts\n",
    "    \n",
    "    for i in range(MCS):\n",
    "        fs = lat.pickrandomsite() # pick a random site to flip its spin\n",
    "        lat.filpspin(fs) # flip spin\n",
    "        dE = lat.sites[fs].getdeltaexchangeenergy(J)*2\n",
    "\n",
    "        if dE <= 0 or np.exp(-dE/(kB*T)) > rd.random():\n",
    "            E_t += dE\n",
    "            MCA += 1 # MC Accepts\n",
    "            \n",
    "        else:\n",
    "            lat.filpspin(fs) # flip the spin back\n",
    "    \n",
    "    print('%.2f  %d' %(T, MCA))\n",
    "    MCA_plot.append(MCA)\n",
    "    T_plot.append(T)\n",
    "    E_plot.append(E_t)\n",
    "    T = T_p * T\n",
    "    \n",
    "print (\"Runtime:  %.6f\"  %(time()-starttime))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "d3c72fd9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEHCAYAAACp9y31AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcvElEQVR4nO3df5Akd3nf8ffnVncSI4Gl+yHl+LGzCA6IisKAFgeCDQEkkGQsyTi4UK2pK6yqrRyEEvlhc8okFE68KWFsB0NyUGtb4so7YDBGJRVOCKoLgtgkwJ6QhM5CPiR2V4Lj7iRBkFmwhO7JH9Nzmtubnp2ZnZ7umf68qqZm5jvdPc/1zD3T+/S3v19FBGZmVh6b8g7AzMyGy4nfzKxknPjNzErGid/MrGSc+M3MSuaMvAPoxvbt22NqairvMMzMRsrBgwcfjogda9tHIvFPTU2xuLiYdxhmZiNF0nK7dpd6zMxKxonfzKxknPjNzErGid/MrGSc+M3MSmZsE3+9XmdqaopNmzYxNTVFvV7POyQzs0IYie6cvarX68zOzrK6ugrA8vIys7OzAMzMzOQZmplZ7sbyiL9Wq51M+k2rq6vUarWcIjIzK46xTPwrKys9tZuZlclYJv7Jycme2s3MymQsE//c3ByVSuWUtkqlwtzcXE4RmZkVx1gm/pmZGebn56lWq0iiWq0yPz/vE7tmZoBGYc7d6enp8CBtZma9kXQwIqbXto/lEb+ZmaVz4jczKxknfjOzknHiNzMrGSd+M7OSyTTxSzpX0mckfUvSvZJeJWmrpNskHU7uz8syBjMzO1XWR/x/BHw+Il4E/DxwL7AXOBARu4ADyXMzMxuSzBK/pGcArwH+FCAiHo+IHwJXAfuTxfYDV2cVg5mZnS7LI/4LgePATZK+IelPJJ0NXBARRwCS+/PbrSxpVtKipMXjx49nGGb/POa/mY2iLBP/GcDLgY9GxMuAH9NDWSci5iNiOiKmd+zYkVWMfWuO+b+8vExEnBzz38nfzIouy8T/EPBQRHw1ef4ZGj8ERyXtBEjuj2UYQ2Y85r+ZjarMEn9EfB94UNILk6Y3AH8L3ArsTtp2A7dkFUOWPOa/mY2qrKdefDdQl7QFeAB4B40fm09LuhZYAd6acQyZmJycZHl5uW27mVmRZZr4I+JO4LSR4Wgc/Y+0ubm5U+b1BY/5b2ajwVfu9slj/pvZqPJ4/GZmY8rj8ZuZGeDEb2ZWOk78ZmYl48RvZlYyTvxmZiXjxG9mVjJO/GZmJePEb2ZWMk78ZmYl48RvZlYyTvxmZiXjxG9mVjJO/GZmJePEb2ZWMk78ZmYl48RvZlYyTvxmZiXjxG9mVjJO/AVUr9eZmppi06ZNTE1NUa/X8w7JzMbIGXkHYKeq1+vMzs6yuroKwPLyMrOzswCeyN3MBsJH/AVTq9VOJv2m1dVVarVaThGZ2bjJ9Ihf0hLwGPAk8LOImJa0FfgUMAUsAb8eET/IMo5RsrKy0lO7mVmvhnHE/7qIeGlETCfP9wIHImIXcCB5bonJycme2s3MepVHqecqYH/yeD9wdQ4xFNbc3ByVSuWUtkqlwtzcXE4Rmdm4yTrxB/AFSQclzSZtF0TEEYDk/vx2K0qalbQoafH48eMZh1kcMzMzzM/PU61WkUS1WmV+ft4nds1sYBQR2W1cemZEfE/S+cBtwLuBWyPi3JZlfhAR53XazvT0dCwuLmYWp5nZOJJ0sKXMflKmR/wR8b3k/hhwM/ALwFFJO5OgdgLHsozBzMxOlVnil3S2pKc3HwNvBO4BbgV2J4vtBm7JKgbLhi8wMxttWXbnvAC4WVLzfT4REZ+X9HXg05KuBVaAt2YYgw2YLzAzG32Z1vgHxTX+4piammJ5efm09mq1ytLS0vADMrNUudT4bfz4AjOz0efEXxKDqsv7AjOz0efEXwLNuvzy8jIRcbIu30/y9wVmZqPPib8EBjnwmy8wMxt9TvwlMOi6/MzMDEtLS5w4cYKlpaW+kr67hJrlx4m/BIpWlx9k6cnMeufEXwJFq8t7zgGzfDnxl0DR6vLuEmqWL1/AZUPni8DMhsMXcFlhFK30ZFY2Tvw2dEUrPZmVjUs9ZmZjyqUeMzMDnPjNzErHid/MrGSc+M3MSsaJ38ysZJz4zcxKxonfzKxknPjNzErGid/MrGTWTfySzpa0KXn8AklXStqcfWhmZpaFbo74vwycJelZwAHgHcDHswzKzMyy003iV0SsAm8BPhIRvwpc1O0bSJqQ9A1Jn0ueb5V0m6TDyf15/YVuZmb96CrxS3oVMAP8VdJ2Rg/vcR1wb8vzvcCBiNhF4y+IvT1sy8zMNqibxH8dcD1wc0QcknQh8MVuNi7p2cAvA3/S0nwVsD95vB+4uutozcxsw7pJ/BdExJUR8QGAiHgA+N9dbv9DwG8DJ9Zs70iyrSPA+e1WlDQraVHS4vHjx7t8OyuLer3O1NQUmzZtYmpqyhO1m/Wgm8R/fZdtp5D0ZuBYRBzsOSogIuYjYjoipnfs2NHPJmxM1et1ZmdnWV5eJiJYXl5mdnbWyd+sS6m1ekmXA1cAz5L04ZaXngH8rIttvxq4UtIVwFnAMyQtAEcl7YyII5J2Asf6D9/KqFarsbq6ekrb6uoqtVrNs3iZdaHTEf/3gEXgp8DBltutwJvW23BEXB8Rz46IKeBtwP+KiN9I1t+dLLYbuKXv6K2UVlZWemo3s1OlHvFHxF3AXZI+AQh4ERDAfRHx+Abe8wbg05KuBVaAt25gW1ZCk5OTLC8vt203s/V1U+O/FLgf+DDwX4FvJ2WgrkXE7RHx5uTxIxHxhojYldw/2nPUVmpzc3NUKpVT2iqVCnNzczlFZDZaukn8fwi8LiL+WUS8Fngd8F+yDcss3czMDPPz81SrVSRRrVaZn58vXH3fPY+sqBQRnReQvhwRr2l5LuBLrW1Zm56ejsXFxWG9ndmGNXsetZ6ErlQqhfyBsvEl6WBETJ/W3kXi/yhQBT5No8b/VuA+4G8AIuKzA492DSd+GzVTU1Ntz0NUq1WWlpaGH5CVUlri72bohbOAo8Brk+fHga3Ar9D4Icg88ZuNGvc8siJbN/FHxDuGEYjZOHHPIyuybsbjf4GkA5LuSZ6/RNK/zz40s9HlnkdWZN306vljGkM0PAEQEXfTuCDLzFKMSs8jK6duavyViPhaozPPSd0M2WBWajMzM070VkjdHPE/LOl5NE7kIumfA0cyjcrMzDLTzRH/u4B54EWSvgt8B/iNTKMyM7PMdNOr5wHgEklnA5si4rHswzIzs6x006vnP0s6NyJ+HBGPSTpP0u8OIzgzMxu8bmr8l0fED5tPIuIHNMbpNzOzEdRN4p+QdGbziaSnAWd2WN7MMuTB32yjujm5uwAckHQTjZ49v8lTk6Wb2RCtHfytOe0k4K6j1rV1B2kDkHQZcAmNCVm+EBH/M+vAWnmQNrMGD/5mveh7kDZJzwVuj4jPJ8+fJmkqIpYGH6aZdeLB32wQuqnx/wVwouX5k0mbmQ1Z2iBvHvzNetFN4j+jdY7d5PGW7EIyszQe/M0GoZvEf1zSlc0nkq4CHs4uJDNL48HfbBC6mYHreUAdeCaNk7sPAm+PiPuzD6/BJ3fNzHrX98ndJMG/UtI5NH4oHpP0CmBoid/MzAanm1JP0yTwW5L+DvjoegtLOkvS1yTdJemQpN9J2rdKuk3S4eT+vD5jNzOzPnRM/JKqkvZKugv4M+CdwBvb/enQxj8Ar4+InwdeClwm6ZXAXuBAROwCDiTPzYbGV74+pUj7okixjL2IaHsDvgIcAv4DsCtp+07a8p1uQAW4A/gnwH3AzqR9J3DfeutffPHFYTYICwsLUalUgsZV6AFEpVKJhYWFvEMbuiLtiyLFMk6AxWiTU1NP7kq6BXgZcCvwiYj4iqQHIuLCbn9UJE0AB4HnA/8tIt4r6YcRcW7LMj+IiI7lHp/ctUHxla9PKdK+KFIs4yTt5G7HXj2Sfg74NeAaGsn7XOBNEfG1Ht/8XOBm4N3AX3eT+CXNArMAk5OTF7f7Upj1atOmTbT7zkvixIkTbdYYX0XaF0WKZZykJf6ONf6I+H8RcWNEXEqjTPM+4EOSHuzlzaMxrPPtwGXAUUk7k6B2AsdS1pmPiOmImN6xY0cvb2eWyle+PqVI+6JIsZRB1716IuJYRHwkIv4p8IvrLS9pR3Kk3xzK+RLgWzRKR7uTxXYDt/QatFm/fOXrU4q0L4oUSym0K/wP4ga8BPgGcDdwD/C+pH0bjd48h5P7rettyyd3bZAWFhaiWq2GpKhWq12dQOxnnaxjGuX33UgsRYq56Eg5uZtZ4h/kzYnf8pR1jxP3aOme91Vv0hJ/V+Px5829eixPWfc4cY+W7nlf9abnXj2Sfg94ICI+tqb9XwH/KCLem0mkbTjxW56y7nHiHi3d877qTT+9et4MzLdp/yPglwcVmFnRZd3jxD1auud9NRidEn9ExGk/oUmbsgvJrFiy7nHiHi3d874akHaF/+RPqa+TDNWwpn0XKScMsrr55K7lbVx79ayniHEVMaaiSsvVnWr8lwMfAX6XxrALANPA9cB7IuK/Z/uT9BTX+M2Gr16vMzs7y+rq6sm2SqXiiV9GSL9DNrwY+C3gxUnTIeCDEfHNTKJM4cRvNnzuQTP6+pqIJSLu4amrbM2sRFZWVnpqt9GRenJX0q2dbsMM0syGL6seNL2Mu+8x+rPR6Yj/VTTm1/0k8FXck8esVObm5trW+DfSg2bteYPl5WVmZ2cBTjtv0Muy1ptOJ3cngEtpDMn8EuCvgE9GxKHhhdfgGr9ZPur1OrVajZWVFSYnJ5mbm9tQ0u3lvIHPMWxcXyd3W1Y+k8YPwAeB/xgRHxl8iOmc+M3GQy9X3voq3Y3razx+SWdKeguwALwL+DDw2WxCNLNx18t5A1+lm51OJ3f305h39+XA70TEKyLiP0XEd4cWnZmNlV6uvPVVutnpdMT/duAFwHXAVyT9KLk9JulHwwnPzPpVxB4xMzMzzM/PU61WkUS1Wk29IKyXZa03HpbZbAz5qluDPmv8ZjaaarXaKUkfYHV1lVqtllNEViRO/GZjyFfdWidO/GZjyD1irBMnfrMx5B4x1okTv9kYKnuPmCL2aCoS9+oxs7HiHk1Pca8eMysF92haX2aJX9JzJH1R0r2SDkm6LmnfKuk2SYeT+/OyisHMysc9mtaX5RH/z4B/ExH/GHgl8C5JFwF7gQMRsQs4kDw3MxuIUe7RNKxzE5kl/og4EhF3JI8fA+4FngVcBexPFtsPXJ1VDGZWPqPao6l5bmJ5eZmIODn/QBbJfyg1fklTwMtoTOhyQUQcgcaPA3D+MGIws3IY1R5Nwzw3kXnil3QO8JfAeyKi68HdJM1KWpS0ePz48ewCNLOxMzMzw9LSEidOnGBpaWnoSb9TySbttbRzEMvLywM/6u842fpGSdpMI+nXI6I5jv9RSTsj4oikncCxdutGxDwwD43unFnGaWY2KJ2mjARSX5ucnGw741hzHRjclJOZ9eOXJBo1/Ecj4j0t7R8EHomIGyTtBbZGxG932pb78ZvZqOg0ZSSQ+lq7OY7XLtPrlJN59ON/NY0x/V8v6c7kdgVwA3CppMM05vS9IcMYzMy6Lr1s376d7du3n/Y4rYdNu+12Ktl06mraPDeRZqDdUSOi8LeLL744zMz6sbCwEJVKJYCTt0qlEgsLC21fS7s111lvu9u2bWu7vqTU16rV6sntVqvVdZfpFrAYbXJq7km9m5sTv5n1q1MiTXst7dZNgt62bVtISn0t7UeoqdMPVa/SEr+HbDCzvvTTcyUPncorvZZPmsvX6/XUE7GPPvpo46g65bVOXU3r9TrXXXfdKXX+bdu2Db47artfg6LdfMRvViy9lk/6PWIdhEEf8a9XHuq03U7lmoWFhdi8efNp62zZsqXvfYdLPWY2KP0k035q1IMw6Bp/px+Ljfz4ddpuv/vOid/MBiathi2p42t5aSbs5gnWZh2+Wq3Gnj172r6WtlynH4e1tfrm8hMTE6fcN/9yaJW23zay75z4zWxgRumIv1W/Zahu/jJo9+/rtN7a903r8eMjfjMrhFGq8bfq90dpvSP9tH/feuu1vm+nxO8av5kVQmv5ZG3potNreeq3DNWpDNPp39dpvbXv22nZfjnxm1npDfqIv5n0037kBnHEv5ESWVridz9+MyuNfsfqT1vviiuu6DiGfrv12r1vvV7nkUceabvc85///K7+bT1p92tQtJuP+M1sUPotQ7Vbr5u/INr17unlL4OJiYm+/62kHPFnNjrnIHl0TjMrok2bNpGWQyUxOTl58qi+VquxsrJyWlvaFcCt9uzZw759+3qOL210Tid+M7M+pQ3B3Grz5s1I4vHHH+/Ytp5+kn8ewzKbmY21TjX8pieeeOK0BN+ubT2dhmzulRO/mY28vAaFWzu/b5aefPLJgW3Lid/MRlpzqsO0njVZa53ft5rMspWFQf6wOPGb2Uir1WqnTVe4urpKrVYbeiztSj+bN29my5YtG9722WefveFtNDnxm9lI6zTe/iD0UkZaW/qpVqvcdNNN3HjjjRsuB/34xz/ue9213KvHzEZap8nNe52cfK1mGan1L4pKpbKhiVG66QnUzllnncVPfvKTntZxrx4zG0v9Xo3bjSzKSHNzc30d+f/0pz/t+z3XcuI3s5HWrrwyqKkKeykjdVsSmpmZSb3oa1hc6jEzS9FtGanXklC/5Z5e87VLPWZmPeq2jNRrSWhubo7NmzcPNtgeZJb4Jd0o6Zike1ratkq6TdLh5P68rN7fzGyjui0j9dqzaGZmhptuuolt27YNPOZuZHnE/3HgsjVte4EDEbELOJA8NzMrrNYLtJaWltqWbiYnJ9uu2669eS7g7W9/O+eccw4LCwunTpKSYmJiov9/xBqZJf6I+DLw6Jrmq4D9yeP9wNVZvb+Z2bB0WxJa7yrjSy65JPU9Nm0aXLoedo3/gog4ApDcnz/k9zczG7huS0LrnQs4cOBA6ns88cQTA4s30149kqaAz0XEi5PnP4yIc1te/0FEtK3zS5oFZgEmJycv7ucMuJlZkXTqv1+tVtft6TOqvXqOStqZBLQTOJa2YETMR8R0REzv2LFjaAGamWWhXq+nJn5JfXXv7NewE/+twO7k8W7gliG/v5lZLmq1WuoRezdH8iMxOqekTwL/B3ihpIckXQvcAFwq6TBwafLczGzsDWrQuEE4I6sNR8Q1KS+9Iav3NDMrqsnJybblnImJia4mWdm6devAYvGVu2ZmQ5DW5XN2djbz2bvWcuI3MxuCtC6f+/bt66rG/+ijay+L6p8HaTMzy1k3g7Zt27aNhx9+uKftFqU7p5mZrdGuDJQlJ34zs5zNzMywe/fujssMstTjxG9mlrN6vc7+/fs7LuPJ1s3Mxki7MXzWGuRk6078ZmY56+birkF2xHHiNzPLWdp4/mulzePbKyd+M7OcddurJ20qx1458ZuZ5az14q5OBjXejxO/mVkBNKd47KTbktB6nPjNzApk8+bNqa+tncqxX078ZmYF0mmKxXYTvffDid/MrGSc+M3MRsQ73/nOgWzHid/MrEAuuuii1Nfm5+cH8h5O/GZmBXLo0KHU17qZqasbTvxmZiXjxG9mVjBpUzEOaopGJ34zs4JJG5BtUAO1OfGbmRXMxMRET+29cuI3MyuY2dnZntp7lUvil3SZpPskfVvS3jxiMDMrqn379rFnz56TR/gTExPs2bOHffv2DWT7GuTg/l29oTQB/B1wKfAQ8HXgmoj427R1pqenY3FxcUgRmpmNB0kHI2J6bXseR/y/AHw7Ih6IiMeBPweuyiEOM7NSyiPxPwt4sOX5Q0nbKSTNSlqUtHj8+PGhBWdmNu7ySPztOqKeVm+KiPmImI6I6R07dgwhLDOzcsgj8T8EPKfl+bOB7+UQh5lZKeWR+L8O7JL0XElbgLcBt+YQh5lZKQ29Vw+ApCuADwETwI0R0XFaGUnHgeU+32478HCf62bJcfXGcfXGcfWuqLFtJK5qRJxWK88l8Q+TpMV23Zny5rh647h647h6V9TYsojLV+6amZWME7+ZWcmUIfEPZsqawXNcvXFcvXFcvStqbAOPa+xr/GZmdqoyHPGbmVkLJ34zs5IZ6cQv6TmSvijpXkmHJF2XtG+VdJukw8n9eS3rXJ8MB32fpDdlGNuNko5Juqel7VOS7kxuS5LuTNqnJP2k5bWPDTmu90v6bsv7X9HyWp7764OSviXpbkk3Szo3ac97f+X+/VoT4wtb9sWdkn4k6T2dPtdhSb7n30zefzFpS91/Q4opLW/kvr/WxJnd8PURMbI3YCfw8uTx02kM93wR8HvA3qR9L/CB5PFFwF3AmcBzgfuBiYxiew3wcuCelNf/AHhf8ngqbblhxAW8H/i3bZbNdX8BbwTOSB5/oOVzzHt/5f796hDvBPB9oJr2uQ45niVg+5q2tvtviDGl5Y3c99eaz/F+4EJgS/K9umhQ2x/pI/6IOBIRdySPHwPupTHS51XA/mSx/cDVyeOrgD+PiH+IiO8A36YxTHQWsX0ZeLTda2rMmPzrwCezeO9OOsXVRq77KyK+EBE/S57+XxrjOg1Vyv7K/fvVwRuA+yOi3yvdhyFt/w1Fh7xRJJkOXz/Sib+VpCngZcBXgQsi4gg0PmTg/GSxroaEHoJfAo5GxOGWtudK+oakL0n6pRxi+pdJSeXGlj+9i7K/AH4T+B8tz/PcX0X+fr2NUw8o2n2uwxTAFyQdlNScNzBt/w3dmrwB+e+vpky/S2OR+CWdA/wl8J6I+FGnRdu05dGf9RpO/c95BJiMiJcB/xr4hKRnDDGejwLPA16axPIHSXsh9pekGvAzoJ405b2/0uS6v9QY9PBK4C+SprTPdZheHREvBy4H3iXpNTnE0FabvFGE/dWU6Xdp5BO/pM00Prx6RHw2aT4qaWfy+k7gWNKe+5DQks4A3gJ8qtmWlAYeSR4fpFHbe8GwYoqIoxHxZEScAP6Yp8oTRdhfu4E3AzORFD/z3l8U9/t1OXBHRByFjp/r0ETE95L7Y8DNSQxp+29o2uWNIuyvFpl+l0Y68Se18j8F7o2IP2x56VZgd/J4N3BLS/vbJJ0p6bnALuBrw4o3cQnwrYh4qNkgaYcacxEj6cIkrgeGFVDzP2HiV4FmD5Zc95eky4D3AldGxGpLe677i+J+v075S7LD5zoUks6W9PTmYxon6+8hff8NK662eSPv/bVGtsPX5332eoNnvn+Rxp8/dwN3JrcrgG3AAeBwcr+1ZZ0ajSPE+4DLM4ztkzT+XHyCxq/3tUn7x4F/sWbZXwMO0ThzfwfwK8OMC/gz4JvJfrwV2FmE/UXj5OiDLZ/txwqyv3L/frWJswI8AvxcS1vq5zqkmC5MPqO7ks+rlrSn7r8hxZWWN3LdX23ivIJGj6P7m/tuUDcP2WBmVjIjXeoxM7PeOfGbmZWME7+ZWck48ZuZlYwTv5lZyTjxm5mVjBO/lZ6kbS1D8X5/zdC8W3rYzu2SptdZ5t9tPGKzjXE/frMWkt4P/H1E/H4f695OY1jfxQ7L/H1EnNN/hGYb5yN+sx6oMQnMtyTtT0Zx/IykSpvlrkkmILlH0geSthuApyV/SdRP27jZkDjxm/XuhcB8RLwE+BHwztYXJT2TxsQxr6cx0uMrJF0dEXuBn0TESyNiZsgxm53kxG/Wuwcj4m+Sxws0xn5p9Qrg9og4Ho2JZOo0ZvIyKwQnfrPerT0xtvZ5u7HUzQrDid+sd5OSXpU8vgb46zWvfxV4raTtyfDR1wBfSl57IhkL3iw3TvxmvbsX2C3pbmArjZmbTorGdILXA18kGTo6Ippjzs8Dd/vkruXJ3TnNepDM0fq5iHhx3rGY9ctH/GZmJeMjfrM2JDVniVrrDZHM92s2qpz4zcxKxqUeM7OSceI3MysZJ34zs5Jx4jczK5n/D9uaUKq5VZ+IAAAAAElFTkSuQmCC\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,MCA_plot,'ok')\n",
    "plt.gca().invert_xaxis()\n",
    "plt.xlabel('T / Kelvin')\n",
    "plt.ylabel('MC Accepts')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "040ddf81",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAEGCAYAAAAjc0GqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbvklEQVR4nO3df5BdZZng8e8DAWYzAwqBHSkg3WFltio4I0gWnfVHucKGQLmCju5mtgPxR1UEnF0sylqgunZwtbK7jOOwuG6SiiM1gWQFxxkXylExjjpW7aoYEIGISNBuyMAg6aACmUKSPPvHfS/c3Nx707fT9557u7+fqlM5/Z73Pf326Vv95D3vc94TmYkkSf12RNUdkCTNTwYgSVIlDECSpEoYgCRJlTAASZIqsaDqDgyTE088MUdHR6vuhiQNlXvuuWdXZp7UXG4A6sLo6Cjbtm2ruhuSNFQiYrJVubfgJEmVMABJkiphAJIkVcIAJEmqhAFIklQJA1CPbdmyhdHRUY444ghGR0fZsmVL1V2SpIFgGnYPbdmyhTVr1rBnzx4AJicnWbNmDQBjY2NVdk2SKucIqIfGx8dfCj51e/bsYXx8vKIeSdLgMAD10GOPPdZVuSTNJwagHlq8eHFX5ZI0nxiAemjt2rUsXLjwgLKFCxeydu3ainokSYPDANRDY2NjbNy4kZGRESKCkZERNm7caAKCJAGRmVX3YWgsW7YsXYxUkroTEfdk5rLmckdAkqRKGIAkSZUwAEmSKmEAkiRVwgAkSaqEAUiSVAkDkCSpEgYgSVIlDECSpEoYgIacL7yTNKx8Id0Q84V3koaZI6Ah5gvvJA0zA9AQ84V3koaZAWiI+cI7ScPMADTEfOGdpGFmABpivvBO0jDzhXRd8IV0ktQ9X0gnSRooBiBJUiUMQJKkShiAJEmVMABJkiphAJIkVcIApLZcaVtSL7katlpypW1JveYISC250rakXhu4ABQRH42Iv4+I+8p2UcOx6yJiR0Q8HBEXNJSfExEPlGOfiogo5cdExO2l/HsRMdrQZnVEPFK21X39IYeAK21L6rWBC0DFjZl5Vtm+DBARS4GVwJnACmBdRBxZ6q8H1gBnlG1FKf8A8Exmvhq4EbihnOsE4Hrg9cC5wPURcXxffrIh4UrbknptUANQKxcDt2XmC5n5M2AHcG5EnAwcl5nfydrCdrcAlzS02VT2vwCcV0ZHFwBbM3N3Zj4DbOXloCVcaVtS7w1qAPqjiLg/Im5uGJmcAjzeUGdnKTul7DeXH9AmM/cCvwQWdTjXQSJiTURsi4htTz/99OH9VEPElbYl9VolWXAR8XXgVS0OjVO7nfZxIMu/nwTeD0SL+tmhnBm2ObAwcyOwEWqrYbeqM1eNjY0ZcCT1TCUBKDPPn069iPgM8KXy5U7gtIbDpwJPlPJTW5Q3ttkZEQuAVwC7S/lbm9p8q5ufQZJ0eAbuFlyZ06l7J/Bg2b8TWFky25ZQSza4OzOfBJ6NiDeU+Z3LgDsa2tQz3N4NfKPME90FLI+I48stvuWlTJLUJ4P4IOqfRMRZ1G6JTQAfBMjM7RHxeeBHwF7gQ5m5r7S5AvgL4J8AXykbwGeBWyNiB7WRz8pyrt0R8XHg+6XexzJzd29/LElSI9+I2gXfiCpJ3fONqJKkgWIAUt+4uKmkRoM4B6Q5yMVNJTVzBKS+cHFTSc0MQOoLFzeV1MwApL5wcVNJzQxA6gsXN5XUzACkvnBxU0nNfBC1Cz6IKknd80FUSdJAMQBJkiphAJIkVcIAJEmqhAFIQ8l15aTh51pwGjquKyfNDY6ANHRcV06aGwxAGjquKyfNDQYgDR3XlZPmBgOQho7ryklzgwFIQ2c215Uzm06qjmvBdcG14OaW5mw6qI2kXCRVml2uBSc1MZtOqpYBSPOW2XRStQxAmrfMppOqZQDSvGU2nVQtA5DmLd/SKlXLLLgumAUnSd0zC07qAZ8jkmbO1bClGXJVbunwOAKSZsjniKTDYwCSZsjniKTDYwCSZsjniKTDYwCSZuhwnyMygUHznQFImqHDeY6onsAwOTlJZr6UwGAQ0nzic0Bd8DkgzZbR0VEmJycPKh8ZGWFiYqJtuy1btjA+Ps5jjz3G4sWLWbt2rRl3GnjtngMyDVuqwEwSGEz71lzjLTipAjNJYDDtW3ONAUiqwEwSGEz71lxTSQCKiPdExPaI2B8Ry5qOXRcROyLi4Yi4oKH8nIh4oBz7VEREKT8mIm4v5d+LiNGGNqsj4pGyrW4oX1LqPlLaHt2HH1t6yUwSGEz71lxT1QjoQeBdwLcbCyNiKbASOBNYAayLiCPL4fXAGuCMsq0o5R8AnsnMVwM3AjeUc50AXA+8HjgXuD4iji9tbgBuzMwzgGfKOaS+GhsbY2Jigv379zMxMXHIeRxfH6G5ppIAlJkPZebDLQ5dDNyWmS9k5s+AHcC5EXEycFxmfidraXu3AJc0tNlU9r8AnFdGRxcAWzNzd2Y+A2wFVpRjbyt1KW3r55IG1kzTvn3eSINq0LLgTgG+2/D1zlL2YtlvLq+3eRwgM/dGxC+BRY3lTW0WAb/IzL0tznWQiFhDbeTlrQ5VbmxsrKuMNzPnNMjajoCa52a6FRFfj4gHW2wXd2rWoiw7lM+kTadzHXwgc2NmLsvMZSeddFK7atJAMnNOg6zTLbjPlEn6j5W5ma5k5vmZ+ZoW2x0dmu0ETmv4+lTgiVJ+aovyA9pExALgFcDuDufaBbyy1G0+lzSndJs55+069VPbAJSZZwNvB/YBX4iI+yLimogY6WF/7gRWlsy2JdSSDe7OzCeBZyPiDWUO5zLgjoY29Qy3dwPfKPNEdwHLI+L4knywHLirHPtmqUtp2ykoSkOrm8w5lwdS32XmtDbgtcB/Ax4F/u9027U51zupjVBeAJ6iFhjqx8bL93gYuLChfBm17LlHgU/z8jJCvwH8JbWEhbuB0xvavL+U7wDe11B+eqm7o7Q9Zjr9Puecc1IaJps3b86FCxcmtdvMCeTChQtz8+bNB9UdGRk5oF59GxkZ6X/HNacA27LF39RprQUXEUcA5wF/CFwEfDczL+ku1A0/14LTMJru+nFHHHEErf4eRAT79+/vR1c1R7VbC65jAIqIN1MLOpdQG33cBvxVZv6yR/0caAYgzWUzXSBVOpR2AahTFtzjwH8HHgLOzszlmXnzfA0+0lzng67qt05ZcG/KzDdm5v/MzKci4jf71itJfdftg65mzOlwHXIOKCJ+H/gs8FuZuTgiXgt8MDOv7EcHB4m34KSa5gdcoTZamu4L+TS/dH0LrsH/oLaszRRAZv4QeMus9k7SUPEBV82Gaa0Fl5mPNxXt60FfJA2Jbh5w9Vad2plOAHo8Iv4lkBFxdER8hFpigqR5aroPuPpwqzqZTgC6HPgQtQU7dwJnla8lzVPTzZjzVp06OeRq2Jm5C3BWUdJL6okGh3rA1be4qpO2WXAR8alODTPzP/akRwPMLDipOz7cKphZFtzlwJuorRS9DbinaZOkjny4VZ10CkAnAxuppWBfChwF3JmZmzJzU4d2kgRM/+FWM+Xmp+kuRnoKtTXhrgauycxbe92xQeQtOGn2+VDr3DfjB1Ej4nXAh4FVwFfw9pukWWSm3PzVNgsuIv4LtRfSPURtFezrMnNvvzomaX4wU27+6jQC+s/UXm9dfxHdvRFxf0Q8EBH396V3kua86TzU6hzR3NTpOaAlfeuFpHlr7dq1LeeA6plyzXNE9dUUAOeIhlzbEVBmTnba+tlJSXPXoTLlDjVH5OhoeE0rC041ZsFJ/dfpVeG33nqrGXRD4HBexyBJlek0R+ToaLh1eiX3xoh4Z0Qc288OSVKjTqspdMqgcyXuwddpBHQztQy4L0fE30bENeVtqJLUN53miA5ndKTqTXclhEXAcuBC4HeBHwBfzczP97Z7g8U5IGmwdFpF4dJLL207d7R///5+dnPeO6w5oMycyszPZeZlmXk28L+AM2a7k5LUjZmOjsD5oUFgFlwXHAFJw6PT6Agwe66PzIKTNK90Gh05PzQYHAF1wRGQNDd0erbI+aHZ1/UIKCL+U8P+e5qO/dfZ7Z4k9U+n+SHnhvqn0y24lQ371zUdW9GDvkhSX7R7tuiiiy7y2aE+6hSAos1+q68laWi0mx/68pe/7NxQH7WdA4qIezPzdc37rb6eL5wDkuY254Z6YyZZcK+NiF9FxLPA75X9+te/27OeSlJF2s0NnXDCCc4L9UCn1zEcmZnHZeaxmbmg7Ne/PqqfnZSkfmg1N3TUUUfx7LPPOi/UAz4HJElFq7mh4447jl//+tcH1HNeaHb4HFAXnAOS5h/nhQ6fKyFI0gw4L9Q7BiBJ6qCbeaErr7zSoNSFBVV3QJIGWX1x0vHxcR577DEWL17Mc889x9TU1AH19uzZw4YNG166XVcPSo3n0IEqGQFFxHsiYntE7I+IZQ3loxHxjxFxX9k2NBw7JyIeiIgdEfGpiIhSfkxE3F7KvxcRow1tVkfEI2Vb3VC+pNR9pLQ9uk8/uqQhNDY2xsTEBPv372diYoLdu3e3rNc8V2SyQmdV3YJ7EHgX8O0Wxx7NzLPKdnlD+XpgDbX3EJ3By8sBfQB4JjNfDdwI3AAQEScA1wOvB84Fro+I40ubG4AbM/MM4JlyDkmalnbzQq1MTk56S66NSgJQZj6UmQ9Pt35EnAwcl5nfydp/MW4BLimHLwY2lf0vAOeV0dEFwNbM3J2ZzwBbgRXl2NtKXUrb+rkk6ZBazQuVmzIt+fxQa4OYhLAkIn4QEX8XEW8uZacAOxvq7Cxl9WOPA2TmXuCXwKLG8qY2i4BflLrN5zpIRKyJiG0Rse3pp58+vJ9M0pzQ6nmhyy+//KCg1GzPnj2sXr3aEVHRsySEiPg68KoWh8Yz8442zZ4EFmfmVEScA/yfiDiT1ouf1m+2tjvWbXlLmbkR2Ai154Da1ZM0v4yNjR2UXPDGN77xpWSFds9Y7tu3DzBJAXo4AsrM8zPzNS22dsGHzHwhM6fK/j3Ao8DvUBulnNpQ9VTgibK/EzgNICIWAK8AdjeWN7XZBbyy1G0+lyTNWGOywsjIyCHr79mzh1WrVs3b0dBA3YKLiJMi4siyfzq1ZIOfZuaTwLMR8YYyh3MZUA9kdwL1DLd3A98o80R3Acsj4viSfLAcuKsc+2apS2nbNihK0ky0midqZ3Jykve9732ceOKJ8+r2XFVp2O+MiJ3A7wN/ExF3lUNvAe6PiB9SSxK4PDPr+Y5XAH8O7KA2MvpKKf8ssCgidgBXA9cClHYfB75fto81nOsa4OrSZlE5hyTNmuZ5oiOPPLJj/RdffJGpqal5lbDgWnBdcC04STO1ZcsW1qxZc9AL7w5lZGSEtWvXDvU8kWvBSVKFGkdE3ZicnGTVqlWceOKJc25EZACSpD6pJyls3rx52vNDdVNTU3PutpwBSJL6rHl+aNGiRRx99KFXBJtrWXMGIEmqQGPK9q5du7j55punfXturtyWMwBJ0gCYye25Yb8tZwCSpAFSvz23aNGiadUf5hW3DUCSNGDGxsbYtWsXmzdvntZtucnJSRYsWEBEDNX8kAFIkgZUN7flGteYu/TSS7nyyiv70cXDYgCSpAHX7W25zGT9+vUDn6RgAJKkIdDtbTkY/CQFA5AkDZH6bbnpBqH6O4gGMQgZgCRpCHWz2va+ffsGciRkAJKkIdTt2nJ79uzhqquu6nGvumMAkqQhVb8dl5ls3rz5kEkKU1NTAzUKMgBJ0hzQmKTQ6d1DgzQKMgBJ0hwyNjbGpk2b2h6fmpoamPRsA5AkzTFjY2Mdb8dNTU0NxMOqBiBJmoNuuummjsfrD6see+yxlY2GDECSNAcdahRU99xzz7Fq1apKRkMGIEmao2666aZpPyu0fv36vo+EDECSNEd1u4bcBz/4wR736EAGIEmaw+rp2VdcccUh6z7//PN9vRVnAJKkeWDdunXTCkIbNmzo2604A5AkzRPr1q1j8+bNRETbOpnJpZde2pcgtKDn30GSNDDGxsYAWLVqVds6mcl73/veA+r3giMgSZpnxsbGDnk7bu/evYyPj/e0HwYgSZqH1q1bd8g6k5OTPe2DAUiS5qnpvMrhzDPP7Nn3NwBJ0jy1du1ajjiicxj40Y9+1LOEBAOQJM1TY2Nj3HLLLR2z4oCezQWZBSdJ89h0suJ6NRfkCEiS5rmxsTGWLl3a9ninF9wdDgOQJInt27e3PbZv376efE8DkCQJ6DzS6UUiggFIkgR0Hun0YqVsA5AkCej8XNDzzz8/69/PACRJAmrPBfWTAUiSBBx64dHZngeqJABFxCci4scRcX9EfDEiXtlw7LqI2BERD0fEBQ3l50TEA+XYp6I8ORURx0TE7aX8exEx2tBmdUQ8UrbVDeVLSt1HStuj+/OTS9Lwuuqqq2b1fFWNgLYCr8nM3wN+AlwHEBFLgZXAmcAKYF1E1NMy1gNrgDPKtqKUfwB4JjNfDdwI3FDOdQJwPfB64Fzg+og4vrS5AbgxM88AninnkKR5r9Pru6empmb1e1USgDLza5m5t3z5XeDUsn8xcFtmvpCZPwN2AOdGxMnAcZn5ncxM4BbgkoY2m8r+F4DzyujoAmBrZu7OzGeoBb0V5djbSl1K2/q5JGleu+mmm/r2vQZhDuj9wFfK/inA4w3HdpayU8p+c/kBbUpQ+yWwqMO5FgG/aAiAjec6SESsiYhtEbHt6aef7vqHk6Rh0ssX0DXr2VpwEfF14FUtDo1n5h2lzjiwF6jPbLVaES87lM+kTadzHXwgcyOwEWDZsmVt60mSutOzAJSZ53c6XpIC3g6cV26rQW00clpDtVOBJ0r5qS3KG9vsjIgFwCuA3aX8rU1tvgXsAl4ZEQvKKKjxXJKkPqkqC24FcA3wjszc03DoTmBlyWxbQi3Z4O7MfBJ4NiLeUOZwLgPuaGhTz3B7N/CNEtDuApZHxPEl+WA5cFc59s1Sl9K2fi5JUp9UNQf0aeBYYGtE3BcRGwAyczvweeBHwFeBD2VmfW2IK4A/p5aY8Cgvzxt9FlgUETuAq4Fry7l2Ax8Hvl+2j5UyqAW/q0ubReUckqRDmM1ngeLlu186lGXLluW2bduq7oYk9VSnF9SNjIwwMTHR7fnuycxlzeWDkAUnSRogndaEm82X0xmAJEkH6NeacAYgSdIB+vUskAFIklQJA5AkqRIGIElSJQxAkqRKGIAkSZUwAEmSKmEAkiRVwgAkSaqEAUiSVAkDkCSpEgYgSVIlDECSpEoYgCRJlTAASZIqYQCSJFXCACRJqoQBSJJ0kKVLl3ZVPhMGIEnSQbZv335QsFm6dCnbt2+fte+xYNbOJEmaU2Yz2LTiCEiSVAkDkCSpEgYgSVIlDECSpEoYgCRJlYjMrLoPQyMingYmZ9j8RGDXLHZnttiv7tiv7tiv7szVfo1k5knNhQagPomIbZm5rOp+NLNf3bFf3bFf3Zlv/fIWnCSpEgYgSVIlDED9s7HqDrRhv7pjv7pjv7ozr/rlHJAkqRKOgCRJlTAASZIqYQCaBRFxWkR8MyIeiojtEXFVKT8hIrZGxCPl3+Mb2lwXETsi4uGIuKCHfbs5In4eEQ82lN0eEfeVbSIi7ivloxHxjw3HNvS5Xx+NiL9v+P4XNRyr8np9IiJ+HBH3R8QXI+KVpbzq61X556upj/+84VrcFxG/iogPd/q99kv5nD9Qvv+2Utb2+vWpT+3+blR+vZr6uaJ8jnZExLWzevLMdDvMDTgZeF3ZPxb4CbAU+BPg2lJ+LXBD2V8K/BA4BlgCPAoc2aO+vQV4HfBgm+OfBP647I+2q9ePfgEfBT7Som6l1wtYDiwo+zc0/B6rvl6Vf7469PdI4B+AkXa/1z73ZwI4sams5fXrY5/a/d2o/Ho1/R4fBU4Hji6fq6WzdX5HQLMgM5/MzHvL/rPAQ8ApwMXAplJtE3BJ2b8YuC0zX8jMnwE7gHN71LdvA7tbHYuIAP4t8LlefO9OOvWrhUqvV2Z+LTP3li+/C5zai+/dbb8YgM9XB+cBj2bmTFcO6Yd2168vOvzdGCTnAjsy86eZ+WvgNmrXbVYYgGZZRIwCZwPfA347M5+E2ocN+Kel2inA4w3NdlLNB+/NwFOZ+UhD2ZKI+EFE/F1EvLmCPv1RudV1c8MtkUG5XgDvB77S8HWV12uQP18rOfA/Nq1+r/2UwNci4p6IWFPK2l2/vmv6uwHVX6+6nn6WDECzKCJ+C/gr4MOZ+atOVVuUVZEP/4cc+EfiSWBxZp4NXA3874g4ro/9WQ/8M+Cs0pdPlvKBuF4RMQ7sBbaUoqqvVzuVXq+IOBp4B/CXpajd77Wf3piZrwMuBD4UEW+poA8ttfi7MQjXq66nnyUD0CyJiKOofYi2ZOZfl+KnIuLkcvxk4OelfCdwWkPzU4En+tXX0p8FwLuA2+tl5ZbNVNm/h9q939/pV58y86nM3JeZ+4HP8PJto0G4XquBtwNjWW6OV329GNzP14XAvZn5FHT8vfZNZj5R/v058MXSh3bXr29a/d0YhOvVoKefJQPQLChzKZ8FHsrMP2s4dCewuuyvBu5oKF8ZEcdExBLgDODufvW3OB/4cWburBdExEkRcWTZP73066f96lD9j0HxTqCe8VXp9YqIFcA1wDsyc09DeaXXi8H9fB0wsu7we+2LiPjNiDi2vk8tqeRB2l+/fvWr5d+Nqq9Xk+8DZ0TEkjKyXUntus2OqrMs5sIGvInasPR+4L6yXQQsAv4WeKT8e0JDm3Fq/2N+GLiwh337HLVh/IvU/jfzgVL+F8DlTXX/ANhOLdPlXuDf9LNfwK3AA+U63gmcPAjXi9ok/uMNv9sNA3K9Kv98tejnQmAKeEVDWdvfa5/6dHr5Hf2w/L7GS3nb69enfrX7u1Hp9WrRz4uoZeg9Wr92s7W5FI8kqRLegpMkVcIAJEmqhAFIklQJA5AkqRIGIElSJQxAUp9ExKKGFY7/oWnF46Ob6l4XEWNNZe+NiE+X/SMiYlNZqqXV0+pExFsj4kuH6NPlEXHZ4f5s0kwsqLoD0nyRtVUTzoLakvvAc5n5p22qL6e2UOxBSsDZABwFvC8P41mKzOzZKySkQ3EEJA2Ysp7c0Zn5dJsqN1F7iPKyzNxfnvS/OSK+XxZGPWC14jJamojyHqNStiMifru8e+YjpexbEXFDRNwdET+paDFazSMGIGnwnE/tyfxW/j1wDrAyX35FxDjwjcz8F8C/Aj5RlpwBIGtrit1BbVkXIuL1wESWtdqaLMjMc4EPA9fPws8itWUAkgbPCg585UOje6m95K1xgcrlwLVRe7Ptt4DfABY3tbsd+HdlfyUNi9A2qS+kew+1F+5JPWMAkgbPubRfPPTH1OaGbo+IM0tZAH+QmWeVbXFmPtTU7jvAqyPiJGovXvtrWnuh/LsP54jVYwYgaYCUoPLjzNzXrk5m/j/gcuBvImIxcBfwH+rZcBFxdos2Se01BH9GbfXlqV70X+qG/8ORBsuFwFcPVSkzv1RGM18F/jXwx8D9JQhNUHt3UbPbqS2v/97Z6qx0OFwNWxogEbGVWnbbk1X3Reo1A5AkqRLOAUmSKmEAkiRVwgAkSaqEAUiSVAkDkCSpEgYgSVIl/j8MCMTqN1APjQAAAABJRU5ErkJggg==\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,E_plot,'ok')\n",
    "plt.gca().invert_xaxis()\n",
    "plt.xlabel('T / Kelvin')\n",
    "plt.ylabel('E / MeV')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "08aa6925",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "200.00  132\n",
      "184.00  141\n",
      "169.28  125\n",
      "155.74  143\n",
      "143.28  130\n",
      "131.82  137\n",
      "121.27  143\n",
      "111.57  128\n",
      "102.64  134\n",
      "94.43  140\n",
      "86.88  135\n",
      "79.93  119\n",
      "73.53  144\n",
      "67.65  134\n",
      "62.24  124\n",
      "57.26  127\n",
      "52.68  133\n",
      "48.46  127\n",
      "44.59  141\n",
      "41.02  136\n",
      "37.74  137\n",
      "34.72  133\n",
      "31.94  129\n",
      "29.39  128\n",
      "27.04  146\n",
      "24.87  120\n",
      "22.88  124\n",
      "21.05  139\n",
      "19.37  125\n",
      "17.82  123\n",
      "16.39  129\n",
      "15.08  127\n",
      "13.88  136\n",
      "12.77  129\n",
      "11.74  131\n",
      "10.80  119\n",
      "9.94  137\n",
      "9.14  128\n",
      "8.41  116\n",
      "7.74  135\n",
      "7.12  131\n",
      "6.55  133\n",
      "6.03  135\n",
      "5.55  126\n",
      "5.10  121\n",
      "4.69  129\n",
      "4.32  132\n",
      "3.97  124\n",
      "3.65  123\n",
      "3.36  144\n",
      "3.09  133\n",
      "2.85  128\n",
      "2.62  140\n",
      "2.41  119\n",
      "2.22  132\n",
      "2.04  125\n",
      "1.88  130\n",
      "1.73  136\n",
      "1.59  121\n",
      "1.46  130\n",
      "1.34  120\n",
      "1.24  129\n",
      "1.14  137\n",
      "1.05  126\n",
      "0.96  138\n",
      "0.89  125\n",
      "0.81  132\n",
      "0.75  133\n",
      "0.69  134\n",
      "0.63  126\n",
      "0.58  126\n",
      "0.54  134\n",
      "0.49  131\n",
      "0.45  123\n",
      "0.42  124\n",
      "0.38  144\n",
      "0.35  123\n",
      "0.33  133\n",
      "0.30  140\n",
      "0.28  123\n",
      "0.25  132\n",
      "0.23  127\n",
      "0.21  128\n",
      "0.20  126\n",
      "0.18  125\n",
      "0.17  132\n",
      "0.15  137\n",
      "0.14  135\n",
      "0.13  125\n",
      "0.12  131\n",
      "0.11  131\n",
      "0.10  128\n",
      "0.09  130\n",
      "0.09  120\n",
      "0.08  126\n",
      "0.07  127\n",
      "0.07  133\n",
      "0.06  135\n",
      "0.06  113\n",
      "0.05  122\n",
      "Runtime:  4.041039\n"
     ]
    }
   ],
   "source": [
    "starttime = time()\n",
    "\n",
    "MCS = 200 # MC steps per temperature\n",
    "J = 100 # meV\n",
    "T = 200 # Kelvin\n",
    "T_p = 0.92\n",
    "kB = 8.617e-2 # meV/K\n",
    "\n",
    "MCA_plot = []\n",
    "T_plot = []\n",
    "E_plot = []\n",
    "\n",
    "# set up the sites\n",
    "lat = lattice(100,100,100) # set up the size of lattice\n",
    "lat.populatesites()     # Assign each site a coordinate\n",
    "lat.initialisespins()   # Assign each site a random spin\n",
    "E_0 = lat.gettotalenergy1(J) # get initial total energy of the system\n",
    "E_t = E_0\n",
    "\n",
    "while (T>0.05):\n",
    "    \n",
    "    MCA = 0    # Initialise MC accepts\n",
    "    \n",
    "    for i in range(MCS):\n",
    "        fs = lat.pickrandomsite() # pick a random site to flip its spin\n",
    "        lat.filpspin(fs) # flip spin\n",
    "        dE = lat.sites[fs].getdeltaexchangeenergy(J)*2\n",
    "\n",
    "        if dE <= 0 or np.exp(-dE/(kB*T)) > rd.random():\n",
    "            E_t += dE\n",
    "            MCA += 1 # MC Accepts\n",
    "            \n",
    "        else:\n",
    "            lat.filpspin(fs) # flip the spin back\n",
    "    \n",
    "    print('%.2f  %d' %(T, MCA))\n",
    "    MCA_plot.append(MCA)\n",
    "    T_plot.append(T)\n",
    "    E_plot.append(E_t)\n",
    "    T = T_p * T\n",
    "    \n",
    "print (\"Runtime:  %.6f\"  %(time()-starttime))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "d9bb8d54",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAERCAYAAABl3+CQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAdPUlEQVR4nO3df5Ac9Xnn8fdHK4NPCneOVwLza0c4ln2lcEbAnhLHMRfOGAuVg2K7nMMZBTl23ZbKwSdXyjmgtspOLrVX9pGkTinDcUqOOnKSDc5ddFBG5pfzg5ALgRVBWDIIhKwFRQSJxTYOikOkfe6P6cWjVXfvzM5M9/z4vKqmdubbPd2Pekfz7NPfb39bEYGZmVmWRWUHYGZm3c2JwszMcjlRmJlZLicKMzPL5URhZma5nCjMzCxX3yYKSbdJOiJpT4Pr/6Kkb0vaK+krnY7PzKxXqF+vo5B0GfD3wB9GxIXzrLsS+BrwbyPiu5LOjIgjRcRpZtbt+raiiIiHgFfq2yT9hKR7Je2S9BeS/mWy6N8DN0fEd5P3OkmYmSX6NlFk2Ap8JiIuBT4H3JK0vxN4p6S/lPSIpLWlRWhm1mUWlx1AUST9GPAzwB9Jmm0+Pfm5GFgJ/BxwHvAXki6MiO8VHKaZWdcZmERBrXr6XkSsTll2CHgkIv4J+I6kfdQSx2MFxmdm1pUG5tRTRLxKLQl8DEA1FyWL/y9wedK+jNqpqANlxGlm1m36NlFI+irwV8C7JB2S9CmgCnxK0m5gL7A+Wf0+YFrSt4E/BX49IqbLiNvMrNv07fBYMzNrj76tKMzMrD36sjN72bJlsWLFirLDMDPrGbt27Xo5IpanLevLRLFixQomJyfLDsPMrGdImspa5lNPZmaWy4nCzMxyOVGYmVkuJwozM8vlRGFmZrmcKBLbt29nxYoVLFq0iBUrVrB9+/ayQzIz6wp9OTy2Wdu3b2dsbIxjx44BMDU1xdjYGADVarXM0MzMSueKAhgfH38jScw6duwY4+PjbduHKxYz61WuKIDnn3++qfZmuWIxs17migIYGRlpqr1Zna5YXK2YWSc5UQATExMsWbLkpLYlS5YwMTHRlu13smKZrVampqaIiDeqFScLM2sXJwpqp3+2bt1KpVJBEpVKha1bt7bttFAnK5ZOVCuuUMysXl/ej2J0dDS6aVLAuX0UUKtY2pGMFi1aRNrvUBIzMzNdFauZdS9JuyJiNG1ZqRWFpLWS9knaL+mGlOWS9HvJ8iclXVJGnK3qZMXS7mql3RWKqxOzPhARpTyAIeA54O3AacBuYNWcddYB3wAE/DTw141s+9JLL41BsW3btliyZEkAbzyWLFkS27ZtW9D2JJ20rdmHpNJjM7POASYj4zu1zIpiDbA/Ig5ExOvAHfzoHtaz1gN/mPw7HgHeIunsogPtZu2uVtpZobSzOnFlYlaeMhPFucALda8PJW3NrgOApDFJk5Imjx492tZAu121WuXgwYPMzMxw8ODBlk5ptXMEWLtGe3lkl1m5ykwUSmmb2yvbyDq1xoitETEaEaPLl6fezc8a0M4KpV3VSRFXzptZtjITxSHg/LrX5wGHF7COtVm7KpR2VSftqEx86sps4cpMFI8BKyVdIOk04Brg7jnr3A1cm4x++mng+xHxYtGB2sK0qzpptTLxqSuzFmX1chfxoDaq6Rlqo5/Gk7ZNwKbkuYCbk+XfAkYb2e4gjXoaBK2OnqpUKqkjuSqVSkP7rlQqISkqlYpHbFnfImfUU6mJolMPJ4r+08oX9kKH/Hp4rw2SvEThKTysJ7TSb7LQU1etdKK7T8T6iROF9b2FdqovtBPdfSLWb5worO8ttFO9yErEFYh1M08KaJZhoRMkNjtRoyditG7QtZMCmnWzoioRVyDW7VxRmLVZsxWCKxDrBq4ozArUbCXS6QrE1Ye1yhWFWck6WYG4+rBGuaIw62KdrEA8oaK1gxOFWRdo5oLCZq4LafRaEJ+esjxOFGY9ppkKpJHqwxcI2nycKMx6UKMVSCPVR6Onp1x1DC4nCrM+1kj10cjpKVcdg82JwqzPzVd9NHJ6qpGqwxVH/3KiMBtwjZyemq/qcMXR35wozAZcI6en5qs6XHH0N19wZ2bzmu/CvfkuAvSFf93PF9yZWUvmqzparThcbXS5rFvfdfIBvBV4AHg2+fnjGesdpHav7CfIuU3f3IdvhWpWrPluG5t3O1rfcrY75H3HllVR3AB8MyJWAt9MXme5PCJWR0ZJZGbla6XiyKs2XGl0iawM0skHsA84O3l+NrAvY72DwLJmt++Kwqy75FUNWdXG7Dpp79m2bVtUKpWQFJVKxdVHG5BTUZSVKL435/V3M9b7DvA4sAsYm2ebY8AkMDkyMtLmQ2hmrcr6cq9UKqlJYmhoKLV9eHjYp6o6oJREATwI7El5rG8iUZyT/DwT2A1c1si+XVGY9Y6saiOrysh6zCYfVxoLk5coOtZHERFXRMSFKY+7gJcknQ2Q/DySsY3Dyc8jwA5gTafiNbNyZPVvVCqVprYze5GfL/prv7I6s+8GNibPNwJ3zV1B0lJJZ8w+B66kVpGYWZ9Jm2Yk64rx4eHh1G0MDQ2ldopv3rzZHeItKitRfBH4gKRngQ8kr5F0jqSdyTpnAQ9L2g08CtwTEfeWEq2ZFS6r0tiyZUtqAjlx4kTqdqanp11ltMhXZptZz9m+fTvj4+M8//zzjIyMMDExwfj4OFNTUw29f2hoiJmZmTfe66vD86/MXlx0MGZmrapWq6lf7nOnCckyW33MVhiz27R0nsLDzPpC2qmqrP6Meu7HmJ9PPZlZ30qbjLARgzhhoScFNLOBNLfKGBoaauh9abeCHWROFGbW1+qH3t5+++2njJjKMjU15VNRCScKMxsYzfZjeEhtjROFmQ2UuRf3pV2XMdexY8fYsGHDwFYXThRmNtDmVhl5BrW6cKIws4FXX2XMN8fUIFYXThRmZnXS5phKM0jVhROFmVmd+lNR8xmU6sKJwsxsjtlTUdu2bXN1gROFmVkmVxc1ThRmZjlcXThRmJk1ZJCrCycKM7MGDWp14URhZtakZquLjRs39nSyKCVRSPqYpL2SZiSlTmubrLdW0j5J+yXdUGSMZmZ5mqkuTpw40dOVRVkVxR7gI8BDWStIGgJuBq4CVgEfl7SqmPDMzBrTaHXRy5VFKYkiIp6KiH3zrLYG2B8RByLideAOYH3nozMza06j1UWvVhbd3EdxLvBC3etDSZuZWVearS7ybpDUi5XF4k5tWNKDwNtSFo1HxF2NbCKlLfO+rZLGgDGAkZGRhmI0M2u32dun5t2CdbayqF+/m3UsUUTEFS1u4hBwft3r84DDOfvbCmyF2j2zW9y3mdmCzX75b9y4kRMnTqSuM1tZ1K/frbr51NNjwEpJF0g6DbgGuLvkmMzMGlKtVue99Wqv9FmUNTz2w5IOAe8B7pF0X9J+jqSdABFxHLgOuA94CvhaROwtI14zs4Xolz6Ljp16yhMRO4AdKe2HgXV1r3cCOwsMzcysrfqhz6KbTz2ZmfWFRiuLzZs3FxhV45wozMwK0EifxfT0dFeegirl1JOZ2SBqZDRUN46EcqIwMyvQbALYsGFD6vJu7K9QRP9dcjA6OhqTk5Nlh2FmlmnZsmVMT09nLh8eHubll18uLB5JuyIidZJW91GYmZVgy5YtPdNf4URhZlaCRkZCdcsoKCcKM7OSzI6EytItVYX7KMzMSpbXX1FUX4X7KMzMutiWLVsyl01PT/PpT3+6wGhO5URhZlayarXK8PBw5vJbb7211FNQThRmZl0gr6qIiFI7tp0ozMy6wHxVRZkd204UZmZdYsuWLUhpN/esKauqcKIwM+sS1WqVTZs2ZS4vq6rITBSSUodJmZlZ59xyyy25p6DKqCryKorfl/SspP8kaVVhEZmZDbj5hssWXVVkJoqIuBj4EHAC+N+SnpB0vaRKYdGZmQ2g+Tq2i64qcvsoImJfRPxmRKwCNgJvAf5E0l+2slNJH5O0V9JM3ikuSQclfStJUr7U2swGRjdVFQ11ZktaBJwJnAUsBY62uN89wEeAhxpY9/KIWJ11abmZWT+ar6oYHx8vLJbcRCHpfZJuAQ4Bvw48DLwrIn6hlZ1GxFMRsa+VbZiZ9bu8qmJqaqqwOPJGPb0AfBF4Crg4Iq6MiNsi4vuFRQcB3C9pl6SxAvdrZla6arXKokXpX9N505O3W96tUH82It5IWZKWRsRrjW5Y0oPA21IWjUfEXQ1u5r0RcVjSmcADkp6OiNTTVUkiGQMYGRlpNEwzs642MzOT2p51z+1OyBv1NAUg6T2Svk2tskDSRcnpqFwRcUVEXJjyaDRJEBGHk59HgB3Ampx1t0bEaESMLl++vNFdmJl1tUolfaCppMI6tBvpzP6vwAeBaYCI2A1c1sGYgFoFI+mM2efAldQ6wc3MBsbExETqtB5FThTY0KiniHhhTlNLNY+kD0s6BLwHuEfSfUn7OZJ2JqudBTwsaTfwKHBPRNzbyn7NzHpNtVol6wZzRQ2TbSRRvCDpZ4CQdJqkz5GchlqoiNgREedFxOkRcVZEfDBpPxwR65LnByLiouTxkxEx0co+zcx6VdbpJyhmmGwjiWIT8KvAudSGya5OXpuZWQEmJrL/Ti5imOy8iSIiXo6IavKX/5kRsSEi0m/uamZmbZc3TBbo+OknZZ37kvR7eW+MiP/QkYjaYHR0NCYnPeOHmfWPvPtUVCoVDh482Or2d2XNgJF3HcUmaqOMvgYcBrKjNDOzjqpUKpmnmTp9+inv1NPZwFZqQ2N/GXgTcHdE3B4Rt3c0KjMzO0leP0VetdEOeRfcTUfErRFxOfAJajPH7pX0yx2NyMzMTlGtVjOXRURH+ynm7cyWdAnwWWAD8A1gV8eiMTOzTGUNk82bFPA3Je0Cfg34c2A0Ij4VEd/uWDRmZpaprGGyeaOeZoADwD8kTbMrCoiIeHfHomqRRz2ZWb8aGhpKnShwaGiI48ePL3i7Cx31dMGC92hmZh1RxmyymYmifopxMzPrDsPDw0xPn3rN89KlSzu2z4YmBTQzs+722muvdWzkkxOFmVkPeeWVVzKXdWrkU96op63JdOBndGTPZmbWtLw7eHZq5FNeRXEbcBGwU9I3JV0v6aKORGFmZg3JGyLbqftoZw6PPWklaZjaHeauAv4V8DfAvRHxtY5E1SIPjzWzfpY3ZUcj3+kZ21zQ8Nj6HU8DX00eSLoUWLugaMzMrCVZI5+Gh4c7sr+GEsVcEbELT+VhZlaKH/7wh021t6qUUU+SbpL0tKQnJe2Q9JaM9dZK2idpv6QbCg7TzKwrvfbaa021t6qs4bEPABcm04A8A9w4dwVJQ8DN1PpFVgEfl7Sq0CjNzCx3eOx/rHv+sTnL/nMrO42I+yNidlKSR4DzUlZbA+yPiAMR8TpwB7C+lf2amfWDvL6ITlx0l1dRXFP3fO5f/O3syP4ktenL5zoXeKHu9aGkLZWkMUmTkiaPHj3axvDMzLrLli1bMpdt3ry57fvLSxTKeJ72+tQ3Sw9K2pPyWF+3zjhwHEhLgWn7yBz3FRFbI2I0IkaXL18+X3hmZj0r7yZGaaOhWpU36ikynqe9PvXNEVfkLZe0EfgQ8P5IH/h7CDi/7vV51O7dbWZmBcpLFBdJepXaX/b/LHlO8vrNrexU0lrgeuDfRMSxjNUeA1ZKugD4W2qnwn6plf2amfWLpUuXpo5y6sQssnnTjHfmWvCaLwOnAw8kVxg+EhGbJJ0D/EFErIuI45KuA+4DhoDbImJvB2MyM7MUC7rgrlUR8Y6M9sPAurrXO4GdRcVlZtYriryWwtOMm5lZLicKM7MetGhRcV/fThRmZj0o697Z0P6L7pwozMx6UKVSyVzW7jvdOVGYmfWgvBsYtftOd04UZmY9qFqtZvZTtPtOd04UZmY9Kquf4sSJE23djxOFmVmPyppFtt13unOiMDPrUUXd6c6JwsysRxV1dbYThZmZ5XKiMDPrUVmjntp91bYThZlZj8oa9ZR31fZCOFGYmfWorOslfB2FmZkB2ddL+DoKMzMDXFGYmdk8XFGYmVmuoiqKUm6FKukm4OeB14HngF+JiO+lrHcQ+AFwAjgeEaMFhmlm1tX6vaJ4ALgwIt4NPAPcmLPu5RGx2knCzOxkfd1HERH3R8Tx5OUjwHllxGFm1sv6vaKo90ngGxnLArhf0i5JY3kbkTQmaVLS5NGjR9sepJlZt+n5PgpJDwJvS1k0HhF3JeuMA8eBrBu8vjciDks6E3hA0tMR8VDaihGxFdgKMDo6Gi3/A8zMulxRFUXHEkVEXJG3XNJG4EPA+yMi9Ys9Ig4nP49I2gGsAVIThZnZoBkaGkpNCpLaup9STj1JWgtcD1wdEccy1lkq6YzZ58CVwJ7iojQz625ZlUNEsH171oma5pXVR/Fl4Axqp5OekHQrgKRzJO1M1jkLeFjSbuBR4J6IuLeccM3Muk+lUslcNj4+3rb9lHIdRUS8I6P9MLAueX4AuKjIuMzMesnExAQbNmxIXTY1NdW2/XTDqCczM1uAarWaee+Jdo58cqIwM+thWfeeaOfIJycKM7MeVsS1FE4UZmY9rIhrKZwozMx6mCsKMzPL5YrCzMxyuaIwM7NcrijMzCyXKwozM8vlisLMzHK5ojAzs1yuKMzMLNfw8HBT7QvhRGFmZrmcKMzMetj09HRT7QvhRGFm1sPcmW1mZrncmW1mZrn6tqKQ9FuSnkzul32/pHMy1lsraZ+k/ZJuKDpOM7Nu188VxU0R8e6IWA18Hfj83BUkDQE3A1cBq4CPS1pVaJRmZl1OUlPtC1FKooiIV+teLgUiZbU1wP6IOBARrwN3AOuLiM/MrFdEpH19ZrcvxOK2balJkiaAa4HvA5enrHIu8ELd60PAT+VsbwwYAxgZGWlfoGZmA65jFYWkByXtSXmsB4iI8Yg4H9gOXJe2iZS2zBQZEVsjYjQiRpcvX96ef4SZWZdbtCj9azyrfSE6VlFExBUNrvoV4B7gC3PaDwHn170+DzjchtDMzPrGzMxMU+0LUdaop5V1L68Gnk5Z7TFgpaQLJJ0GXAPcXUR8Zma9ooi5nsrqo/iipHcBM8AUsAkgGSb7BxGxLiKOS7oOuA8YAm6LiL0lxWtmNrDUzp7xbjE6OhqTk5Nlh2Fm1nF5w2Cb+X6XtCsiRtOW+cpsM7Me1rfXUZiZWXsUcR2FE4WZmeVyojAz62FFXEfhRGFm1sP69joKMzNrj76dZtzMzNqjn6cZNzOzNnBFYWZmuVxRmJlZLlcUZmaWyxWFmZnl8hQeZmaWy1N4mJlZ6ZwozMx6mKfwMDOzXJ7Cw8zMcnl4rJmZ5SpieGwp98yW9FvAemr3zD4CfCIiDqesdxD4AXACOJ51mz4zs0E1NDSUmhT6oaK4KSLeHRGrga8Dn89Z9/KIWO0kYWZ2qr694C4iXq17uRRo34BfM7MBUkQfRSmnngAkTQDXAt8HLs9YLYD7JQXw3yNia872xoAxgJGRkTZHa2bWnYqoKNTOq/dO2rD0IPC2lEXjEXFX3Xo3Am+OiC+kbOOciDgs6UzgAeAzEfHQfPseHR2NycnJFqI3M+sNixcvzuyjOH78eMPbkbQr6xR/xyqKiLiiwVW/AtwDnJIoZju4I+KIpB3AGmDeRGFmNij6to9C0sq6l1cDT6ess1TSGbPPgSuBPcVEaGbWGyqVSlPtC1HWqKcvStoj6UlqCWAz1E41SdqZrHMW8LCk3cCjwD0RcW854ZqZdaeJiQmWLFlyUtuSJUuYmJho2z5K6cyOiI9mtB8G1iXPDwAXFRmXmVmvqVarAIyPj/P8888zMjLCxMTEG+3t0LHO7DK5M9vMrDl5ndmewsPMzHI5UZiZWS4nCjMzy+VEYWZmuZwozMwsV1+OepJ0FJha4NuXAS+3MZx2cVzNcVzNcVzN6ce4KhGxPG1BXyaKVkia7MYpzR1XcxxXcxxXcwYtLp96MjOzXE4UZmaWy4niVJn3vCiZ42qO42qO42rOQMXlPgozM8vlisLMzHI5UZiZWa6BShSSzpf0p5KekrRX0ux9MN4q6QFJzyY/f7zuPTdK2i9pn6QPdiiu2yQdkbSnru1OSU8kj4OSnkjaV0j6h7plt3Yippy4fkPS39btf13dso4fq3liu0nS05KelLRD0luS9rKPWamfr5QY31V3LJ6Q9Kqkz+b9bouSfNa/lex/MmnLPH4FxZT1vVH68ZoT59rkc7Rf0g1t3XhEDMwDOBu4JHl+BvAMsAr4L8ANSfsNwJeS56uA3cDpwAXAc8BQB+K6DLgE2JOx/HeAzyfPV2StV0RcwG8An0tZt5BjNU9sVwKLk+dfqvs9ln3MSv18zRPvEPB3QCXrd1twPAeBZXPaUo9fgTFlfW+Ufrzm/B6fA94OnJZ8rla1a/sDVVFExIsR8Xjy/AfAU8C5wHrg9mS124FfSJ6vB+6IiH+MiO8A+6ndt7vdcT0EvJK2TJKAXwS+2u79zicvrhSFHKu82CLi/oiYvZv8I8B5ndp/M3FR8udrHu8HnouIhc5kUISs41eInO+NbrIG2B8RByLideAOasetLQYqUdSTtAK4GPhr4KyIeBFqHwrgzGS1c4EX6t52iOI/IO8DXoqIZ+vaLpD0N5L+XNL7Co4H4Lrk9M5tdacBuuFY1fsk8I2612Ues27+fF3DyX+EpP1uixTA/ZJ2SRpL2rKOX+HmfG9A+cdrVkc/SwOZKCT9GPB/gM9GxKt5q6a0FT2e+OOc/B/5RWAkIi4Gfg34iqR/XmA8/w34CWB1EsvvJO3dcKwAkDQOHAe2J01lH7MspR4zSacBVwN/lDRl/W6L9N6IuAS4CvhVSZeVEEOqlO+Nbjheszr6WRq4RCHpTdR+2dsj4o+T5pcknZ0sPxs4krQfAs6ve/t5wOECY10MfAS4c7YtOU0xnTzfRe285DuLiikiXoqIExExA/w+PzpVUuqxmiVpI/AhoBrJyduyjxld+vmi9mX8eES8BLm/28JExOHk5xFgRxJD1vErTNr3Rjccrzod/SwNVKJIzvf/D+CpiPjdukV3AxuT5xuBu+rar5F0uqQLgJXAo0XFC1wBPB0Rh2YbJC2XNJQ8f3sS04GiApr9D5v4MDA7uqfsY4WktcD1wNURcayuvdRjRvd+vk6qVnN+t4WQtFTSGbPPqQ1O2EP28SsqrtTvjbKP1xyPASslXZBUitdQO27tUXZvfZEP4GeplWNPAk8kj3XAMPBN4Nnk51vr3jNO7S/QfcBVHYrrq9RK13+i9pfBp5L2/wlsmrPuR4G91EY1PA78fAeP1ylxAf8L+FZyDO8Gzi7yWM0T235q52lnf7e3dskxK/XzlRHnEmAa+Bd1bZm/24JienvyO9qd/L7Gk/bM41dQXFnfG6Uer5Q411EbkfXc7LFr18NTeJiZWa6BOvVkZmbNc6IwM7NcThRmZpbLicLMzHI5UZiZWS4nCrM5JA3XzQj6d3NmCD1tzro3SqrOafuEpC8nzxdJuj2Z4iHt6lkk/Zykr88T0yZJ17b6bzNbiMVlB2DWbaJ2FfdqqE0lDfx9RPx2xupXUpu08RRJYrgVeBPwK9HCWPSI6NjU6GbzcUVhtkDJfFGnRcTRjFW2ULtY7NqImEmuPL5N0mPJBIUnze6ZVB8HldxHI2nbL+ms5N4Hn0va/kzSlyQ9KumZkiaGtAHiRGG2cFdQu1I4zS8BlwLXxI+mPh8H/iQi/jVwOXBTMlUFAFGbM+guatNBIOmngIORzMU0x+KIWAN8FvhCG/4tZpmcKMwWbi0nT2Ve73FqNwOqnyjuSuAG1e5W+GfAm4GROe+7E/h3yfNrqJsQco7ZCS13Ubsxk1nHOFGYLdwasifxe5pa38Wdkn4yaRPw0YhYnTxGIuKpOe/7K+AdkpZTu0HPH5PuH5OfJ3Bfo3WYE4XZAiRf/k9HxImsdSLi/wGbgHskjQD3AZ+ZHf0k6eKU9wS16bV/l9pspdOdiN+sGf5LxGxhrgLunW+liPh6Uh3cC3wA+DzwZJIsDlK7d8Zcd1KbNvoT7QrWrBWePdZsASQ9QG0004tlx2LWaU4UZmaWy30UZmaWy4nCzMxyOVGYmVkuJwozM8vlRGFmZrmcKMzMLNf/Bz9ZjBQjAPYYAAAAAElFTkSuQmCC\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,E_plot,'ok')\n",
    "plt.gca().invert_xaxis()\n",
    "plt.xlabel('T / Kelvin')\n",
    "plt.ylabel('E / MeV')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "598afa8d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAgPUlEQVR4nO3dfZRkdX3n8fenn9CelYhNa+RhugcX1owkK8yEbJaFmGMC4ipEMbuYPusYwT40moR1SZTtPWrU2T1ITDa6h4XeDUfcbtCTo24IatTjrrBHk8UemMFhBZ3BGRhEaGA3xNMuCzPf/aNuFzU1t6qrqu9DPXxe59wz1b+6Vfd7H6q+U/X71u+niMDMzAxgqOwAzMysezgpmJlZlZOCmZlVOSmYmVmVk4KZmVWNlB3ARpx44okxPT1ddhhmZj1l165dT0bEZNp9PZ0UpqenWV5eLjsMM7OeIulgo/v89ZGZmVU5KZiZWZWTgpmZVTkpmJlZlZOCmZlVOSmYmSWWlpaYnp5maGiI6elplpaWyg6pcD1dkmpmlpWlpSVmZ2dZXV0F4ODBg8zOzgIwMzNTZmiF8icFMzNgfn6+mhDWrK6uMj8/X1JE5XBSMDMDHn744bba+5WTgpkZsHnz5rba+5WTgpkZsHPnTsbHx49qGx8fZ+fOnSVFVI7ckoKkmyU9IWlvyn3XSApJJyZ/T0v6qaTdyXJjXnGZmaWZmZlhYWGBqakpJDE1NcXCwsJAdTIDKK85miWdD/wE+ExEnFnTfirwn4FXA9si4klJ08Adteu1Yvv27eEB8czM2iNpV0RsT7svt08KEXEX8HTKXX8C/AGQTzYyM7OOFdqnIOli4NGI2JNy9xZJ90q6U9J5TZ5jVtKypOWVlZX8gjUzG0CF/XhN0jgwD1yQcvdjwOaIeErSNuC/SnpNRDxTv2JELAALUPn6KM+YzcwGTZGfFF4FbAH2SDoAnALcI+lnI+LZiHgKICJ2AfuBMwqMzczMKDApRMR3I+LlETEdEdPAIeDsiPixpElJwwCSTgNOBx4qKjazQeLxfXpH/bm66qqr8j93EZHLAtxG5Wuh56gkgMvr7j8AnJjcvhS4H9gD3AO8uZVtbNu2LcysdYuLizE+Ph5UCj0CiPHx8VhcXCw7NKuTdq7ql07PHbAcDd5XcytJLYJLUs3aMz09zcGDx07POzU1xYEDB4oPyBpqdK7qdXLuSilJNbPu4/F9eker5yTrc+ekYDZAPL5P72j1nGR97pwUzAaIx/fpHWnnql4e585JYcD1WyVKv+1P1jy+T+9IO1dzc3O5nzt3NA+w+pmmoPI/j159k+i3/THLS7OOZieFAdZvlSj9tj9meXH1kaXqt0qUftsfszI4KQywfqtE6bf9MSuDk8IA67dKlH7bH7MyOCkMsH6rROm3/TErgzuazcwGjDuazcysJU4KZmZW5aRgZmZVTgpmZlblpJADj79jlh2/noo1UnYA/aZ+/J2DBw8yOzsL4NJIszb59VQ8l6RmzOPvmGXHr6d8uCS1QB5/xyw7fj0Vz0khYx5/xyw7fj0Vz0khYx5/xyw7fj0VL7ekIOlmSU9I2pty3zWSQtKJNW3XSton6UFJF+YVV948/o6tx9U0rfPr6djr5aqrrsr3+omIXBbgfOBsYG9d+6nAV4GDwIlJ21ZgD3AcsAXYDwyvt41t27aFWS9ZXFyM8fHxAKrL+Ph4LC4ulh2adaG066V+6eT6AZajwftqbp8UIuIu4OmUu/4E+INkh9ZcAnw2Ip6NiB8C+4Bz8orNrCzz8/NHTRcKsLq6yvz8fEkRWTdLu17qZX39FNqnIOli4NGI2FN318nAIzV/H0ra0p5jVtKypOWVlZWcIjXLh6tprB2tXhdZXj+FJQVJ48A88MG0u1PaUn9AERELEbE9IrZPTk5mGaJZ7lxNY+1o9brI8vop8pPCq6j0F+yRdAA4BbhH0s9S+WRwas26pwA/KjA2s0K4msbakXa91Mv6+iksKUTEdyPi5RExHRHTVBLB2RHxY+B24DJJx0naApwO3F1UbGZFcTVN7yuyeiztepmbm8v1+sltmAtJtwGvA04EHgc+FBF/VnP/AWB7RDyZ/D0PvAt4Hrg6Ir6y3ja6cZgLM+tf9WMxQeV/6r2W2JsNc+Gxj8zMWtQvYzF57CMzswwMQvWYk4KZWYsGoXrMScHMrEWDUD3mpGBm1qKyq8eKqHxyR7OZWQ/IsvLJHc1mZj2uqHGznBTMzHpAUZVPTgpmZj2gqMonJwUzsx5QVOWTk4J1Fc9KZhvRK9dPp3G++MUvbvp3JhrNvtMLi2de6y+elcw2oleun07ibDYDW9Yzr7kk1bpGv4wrY+XoleunkzgbPaaVx6bxgHjWE4aGhki7HiVx5MiREiKyXtIr108ncTZ6TCuPbbC+f6dg3W8QxpWx/PTK9dNJnOvtQ6/OvGbW1CCMK2P56ZXrp5M4m83A1rMzr5mtp+xxZay31FfwAF13/aRVGXV6nddWGg0NVd66c9nHRj3QvbC4+shsMPVCpVFWMeaxr7j6yMz6SS9UGmUVYx776uojM+srvVBplFWMeeyrq4/MrK/0QqVRVjEWva9OCmbWc3qh0iirGAvf10adDRtdgJuBJ4C9NW0fBe4DdgNfA05K2qeBnybtu4EbW9mGO5ptECwuLsbU1FRIiqmpqa7qTC1TLxyXVmKsXWdiYiImJiaOWb92nU2bNsXQ0FAAMTw8HHNzc23HRZOO5jyTwvnA2XVJ4fia27+79uafJIW97W7DScH6XS9U2Vjnmo1plHau5+bmUtdrNzE0Swq5djRLmgbuiIgzU+67FtgcEXPN1mvGHc3W73qhysY6t96YRnD0uR4ZGeHw4cPHrDM8PMzzzz/f8na7qqNZ0k5JjwAzwAdr7toi6V5Jd0o6r8njZyUtS1peWVnJPV6zMhU125aVo5XzWLtOWkJo1t6JwpNCRMxHxKnAEvDepPkxKp8azgLeB9wq6fgGj1+IiO0RsX1ycrKYoM1K0gtVNta5Vs5j7TrDw8Op6zRq70SZ1Ue3ApcCRMSzEfFUcnsXsB84o8TYzLpCL1TZWOeajWkEx57r2dnZ1PUatXei0KQg6fSaPy8GHkjaJyUNJ7dPA04HHioyNrNO5TnbV1njQfXKDGa9rv78TkxMMDEx0fK5lsTc3Bw33HBDdkE16oHe6ALcRuVroeeAQ8DlwOeBvVTKUv8SODlZ91LgfmAPcA/w5la24eojK1s/Vgf14z71g6wqjyJKrD7Km6uPrGz9WB3Uj/vUD7KqPIINVh9J2iRpKLl9hqSLJY22FYFZn+rH6qB+3Kd+UETlEbTWp3AX8CJJJwPfAH4b+HSmUZj1qH6sDurHfeoHRVQeQWtJQRGxCrwV+FREvAXYmmkUZj2qH6uD+nGf+kERlUfQYlKQ9MtUfmz2paRtJNMozHpUP84Wl/c+FVnZ1EtVVGmx1rZ9+ctf5qSTTjrqMSMjI5x77rnZBtKoBzpeqCI6H7gdeH/y92nAJ9d7XBGLq4/MekuRlU29VEWVFuvo6GiMjY01HBdpbRkZGWl7n9hI9ZGk34yIP1+vrQyuPjLrLUVWNvVSFVUrYyA10+4+bWjmNUn3RMTZ67WVwUnBrLcUOWNaL8zOtqZRrK1qd5+aJYWGfQOSLgLeCJws6ZM1dx0PtFcUa2ZGpYIp7X/EeVQ2FbmtjWoUazuPz0qzjuYfAcvA/wV21Sy3AxdmFoGZDYwiK5t6qYoqLdbR0VHGxsbWfezIyEi2+9Sos2FtAUaBMeAXgJ8HxtZ7TFFLpx3NvTBjk1m/KvL110uv9bQZ2EhmV6v9t355/etf3/a22MjMa1S+QnoE+CZwJ/AwcNF6jyti6SQp9FJFgpkNnvVmY0tbCp15TdIDwJsiYl/y96uAL0XEqzf2GWXjOulo7qWKBDMbPJ1UIhU989oTawkh8RDwRMtb7zIe18XMulkn70VFz7x2v6QvS3qnpB1Uhrz+jqS3SnprZpEUxOO6mFk36+S9qOiZ114EPA78CvA6YAV4GfBm4E2ZRVKQXqpIMLPBs95sbGkyHf+oUWdDLyyuPjKzXtHsfaf+vrm5udS/Selk3rp1a9uxsMHqozOoDJm9N/n7F4B/s97jilg89pGZ9YJmVY+tVkQuLi7G0NBQV1Qf3Qn8PnBTRJyVtO2NiDOz+7zSGQ9zYWa9oFnVI9BSRWSzqqQsq49aGQJ7PCLullTb5mEuzMxa1EnVY/19zdYtuvroyeS3CQEg6W3AY5lFYGbW55pVPbZaEdmsKqno6qP3ADcBr5b0KHA1MJdZBGZmfa5Z1WOrFZE7d+5kaCj9LbuU6iNgE/CSNta/mcqP3PbWtH0UuA/YDXwNOKnmvmuBfcCDwIWtbMMdzd3FVV1mjbVTfZTWybw2FlL9UsbYR/8WeGnN3ycAH2vhcecDZ9clheNrbv8ucGNyeyuwBzgO2ALsB4bX24aTQvfwmFJm+VhcXIzR0dGG4x518jprlhRa+froooj4PzWfLP43lUHy1vsEchfwdF3bMzV/bkp2CuAS4LMR8WxE/JDKJ4ZzWojNusT8/Dyrq6tHta2urjI/P19SRGb9YX5+nueee67h/Vm/zlqpPhqWdFxEPAsg6cVU/kffEUk7gXcAfwv8atJ8MvA3NasdStrSHj8LzIKHpugmHlPKLB+tvIayfJ218klhEfiGpMslvQv4OnBLpxuMiPmIOBVYAt6bNCtt1QaPX4iI7RGxfXJystMwLGMeU8osH628hoqaeQ2AiPg48DHg54DXAB9N2jbqVuDS5PYh4NSa+06hMvOb9QiPKWWWj507dzI6Otrw/qxfZ+smBUlbgG9GxDUR8a+AuyRNd7IxSafX/Hkx8EBy+3bgMknHJds7Hbi7k21YOWZmZlhYWGBqagpJTE1NsbCwwMzMTNmhmVUtLS0xPT3N0NAQ09PTLC0tFfL87W63dv35+XmuuOIKNm3adMx6ubzOGvVAxwtVQsvUTMFJZWrO77TwuNuo/MjtOSqfBC4HPg/spVKW+pfAyTXrz1OpOnqQFmd2c/WRmbUq7wq5Rs8/NzfX1nbTnmd0dDTGxsYyi50Njn20OyJeW9e2JyL+YWdpKDse+8jMWpX3rIuNnn94eDh1GIpG221n5rVOY9/ozGsrki6uebJLgCfbjsLMrER5V8g1ep5G4xJlEU8e1X2tJIUrgX8t6WFJjwDvJykJNTPrFXlXyDV6nkbjEmURTx7Vfa1UH+2PiH9E5VfHWyPiH1OZec3MrGfkXSHX6PlnZ2fb2m7a84yOjjI2Ntbyc2xIo86G+oVKUvgI8H2adFIUubij2bLksZv633oznG30nDe6htLaWx0LaWJiojru0fDw8FGdzcPDw21PsBOxgbGPgCngA1TGJdpFpS9hutljilycFCwrHrtp8JR5ztuZba1+vbSlkJnXJH0b+Bngs1TGJfqBpB9GxJbMPqZskKuPLCt5V6ZY9ynznLe67VYrkbKcea1Zn8IK8BLgFcDaeBLN61fNepTHbho8ZZ7zVrfdaiyFzLwWEZcAPw/cA/yhpB8CJ0jy6KXWdzx20+Ap85xnMdtarcJmXouIv42ImyPi14FfAj4I/PukNNWsb3jspsFT5jlvZ7a1+vXSlDLzWu0CTHXyuKwXdzRbllx9lJ1OjmUZxz+vbbbyvK1uu369k0466ahO5q1bt7YdHxuZea2bFycFs+7TSVVPP1V/5bkvc3Nz5VUf9QJXH5l1n06qevqp+ivPfRkZGUntVM6y+shJwcwyNTQ0RNr7iiSOHDmS2WO6VZ77IqXNR1bRznt5RyWpkj4u6cqU9n8p6bqWt25mA6WTqp5+qv7Kc18aVRkVVX30JmAhpf1PgX+aWQRm1lc6qerpp+qvPPelUZVRIdVHwP2d3Ffk4o5ms+7UK9VHednIvjQbP2nTpk3HdDJnXX3UbJiL7wC/FRE/qGs/HbgtGnwfVST3KZhZP1laWmJ2dpbV1dVq2/j4ODt27OCmm25q2CcxNzfHDTfc0PJ2OupolnQR8CngY1QGwwPYDlwLXB0RX245gpw4KZhZP2l39rba+7OqPhpp9KCI+Iqk3wB+H/idpPl+4NKI+G7LWzczs5a0O3tbq/e3o2FSAIiIvcCOzLZmZmYNbd68ueNPCllpVpJ6e7MlswjMzAxoPnvb0FDjYtEsq4+alaT+MnAK8D+APwI+Ubc0JelmSU9I2lvTdr2kByTdJ+mLkl6atE9L+qmk3cly4wb2qW8tLS0xPT3N0NAQ09PTLC0tlR2S2UDL+jU5MzPDwsICU1NTSGJqaoqFhQXOPfdcTjjhhNTHbN26ta1O5nU1KksChoE3ALcA91LpcH5No/VTHn8+cDawt6btAmAkuX0dcF1ye7p2vVaXQSpJ7aexYcz6QVGvyVZmXyt87CNJxwFvB64HPhIRn2ol4UiaBu6IiDNT7nsL8LaImGm2XjODVH3UT2PDmPWDol6Trcy+VtTMa0g6TtJbgUXgPcAngS+0vOXm3gV8pebvLZLulXSnpPOaxDQraVnS8srKSkahdD/PDGbWXYp6TbbyfIXMvCbpFuDbVL4C+sOI+MWI+GhEPLrRjUqaB54H1r6AewzYHBFnAe8DbpV0fNpjI2IhIrZHxPbJycm0VfpSP40NY9YPinpNtvJ8RY199C+AM4DfA74t6Zlk+TtJz3S6QUk7qIyrNJN8t0VEPBsRTyW3dwH7k21bop/GhjHrB0W9JluZfa30mddaXajrQKbScf2/gMm69SaB4eT2acCjwMvWe/5B6miO6K+xYay3Ddq12Gw8oryOQ+1zT0xMxMTExDEdzJLa7mSOaN7RnGdCuI3K10LPAYeAy4F9wCPA7mS5MVn3Uiq/lt4D3AO8uZVtDFpSMOsGg1YJV8b+pm1zdHQ0xsbGMomjWVLwJDtm1pZBq4QrY39bqTjaSByeec3MMtNPs6S1ooz9bbTNNJ3E0XFJqplZvUGrhCtjf9t57qzjcFIws7YMWiVcGfubts3R0VHGxsZyj8NJwcza0mh8npmZmbJDy0Ur+5vHGEg7duyo/v5geHiYK664gvPOO/p3vXl8feU+BTOzDWg0W9pGEmXacw4NDaUmgZGRET796U+3tS13NJuZ5SSP6qR2qo862ZY7ms3McpLHGEjtPjbL8ZacFMzMNiCP6qR2H5tlBZKTgpnZBuRRnZT2nI1mXhsZGcm0AslJwcxsA/Koxkp7zs985jPMzc0dtZ4k3v3ud2da+eWOZjOzHpFVpZM7ms3M+sD8/PxRCQFgdXWV+fn5zLbhpGBm1iOKmO3NScHMrEcUMQ6Tk4KZWY8oYhwmJwUzsy7SaBylpaWlap/C2phIeYw7NZLZM5mZ2YbUVxcdPHiQ2dlZvvWtb3HLLbdU2w8fPlz9hJD1QIQuSTUz6xKNxjwaHh7m8OHDx7R3Or6SS1LNzHpAoyqitITQbP2NcFIwM+sSjaqI1voQWl1/I5wUzMy6RKPqotnZ2cJmf8stKUi6WdITkvbWtF0v6QFJ90n6oqSX1tx3raR9kh6UdGFecZmZtSLr2dRa0WgcpRtuuOGomdjghV8yZx5XROSyAOcDZwN7a9ouAEaS29cB1yW3twJ7gOOALcB+YHi9bWzbti3MzLK2uLgY4+PjAVSX8fHxWFxc7Jp4NhIXsBwN3ldzrT6SNA3cERFnptz3FuBtETEj6VqAiPh3yX1fBT4cEX/d7PldfWRmechjNrU84lnTLzOvvQv4SnL7ZOCRmvsOJW3HkDQraVnS8srKSs4hmtkgKmKMoXast92eH/tI0jzwPLD2ZZhSVkv9CBMRCxGxPSK2T05O5hWimQ2wIsYYasd62+3psY8k7QDeBMzEC99dHQJOrVntFOBHRcdmZgbFjDG00XjW9PTYR5LeALwfuDgiagcFvx24TNJxkrYApwN3FxmbWa8ro1qmX+Uxm9qads9T2phHeY59lGf10W3AY8BzVD4JXA7so9J3sDtZbqxZf55K1dGDwEWtbMPVR2YV3VYtY+naPU95nVfKqj7Km6uPzCq6rVrG0rV7nvI6r82qj5wUzPrA0NAQaa9lSRw5cqSEiCxNu+cpr/ParSWpZpaRbquWsXTtnqcyzquTglkf6LZqGUvX7nkq5bw26mzohcUdzWYvWFxcjKmpqZAUU1NT7mTuUu2epzzOK+5oNjOzNe5TMDOzljgpmJlZlZOCmZlVOSmYmVmVk4KZmVU5KZiZWZWTgpmZVTkpmJlZlZOCmZlVOSmYmVmVk4KZWUm6cba8kbIDMDMbREtLS8zOzrK6WpmZ+ODBg8zOzgJkO71mm/xJwcysBGvzLtdaXV1lfn6+pIgqnBTMzErw8MMPt9VeFCcFM7MSdOtseU4KZmYl6NbZ8nJLCpJulvSEpL01bb8p6X5JRyRtr2mflvRTSbuT5ca84jIz6wYzMzMsLCwwNTWFJKamplhYWCi1kxnIb+Y1SecDPwE+ExFnJm0/BxwBbgKuiYjlpH0auGNtvVZ55jUzs/Y1m3ktt5LUiLgrebOvbfteElBemzUzsw3opj6FLZLulXSnpPMarSRpVtKypOWVlZUi4zMz63vdkhQeAzZHxFnA+4BbJR2ftmJELETE9ojYPjk5WWiQZmb9riuSQkQ8GxFPJbd3AfuBM8qNysxs8HRFUpA0KWk4uX0acDrwULlRmZkNntw6miXdBrwOOFHSIeBDwNPAp4BJ4EuSdkfEhcD5wEckPQ8cBq6MiKfzis3MzNLlWX309gZ3fTFl3c8Dn88rFjMza01XfH1kZmbdwUnBzMyqnBTMzKzKScHMrIfkPVubZ14zM+sRRczW5k8KZmY9oojZ2pwUzMx6RBGztTkpmJn1iCJma3NSMDPrEUXM1uakYGbWI4qYrS23mdeK4JnXzMza12zmNX9SMDOzKicFMzOrclIwM7MqJwUzM6tyUjAzs6qerj6StAIc3MBTnAg8mVE4WXJc7XFc7XFc7enHuKYiYjLtjp5OChslablRWVaZHFd7HFd7HFd7Bi0uf31kZmZVTgpmZlY16ElhoewAGnBc7XFc7XFc7RmouAa6T8HMzI426J8UzMyshpOCmZlV9W1SkHSqpP8u6XuS7pf0e0n7yyR9XdIPkn9PqHnMtZL2SXpQ0oU5xXWzpCck7a1p+5yk3clyQNLupH1a0k9r7rsxj5iaxPVhSY/WbP+NNfflfqyaxHW9pAck3Sfpi5JemrSXfbxKvbZSYvwHNcdit6RnJF3d7LwWJbnOv5tsfzlpa3j8Coqp0XtG6cerLs43JNfRPkkfyHwDEdGXC/BK4Ozk9kuA7wNbgY8DH0jaPwBcl9zeCuwBjgO2APuB4RziOh84G9jb4P5PAB9Mbk83Wq+IuIAPA9ekrFvIsWoS1wXASHL7uppzWPbxKvXaWifeYeDHwFSj81pwPAeAE+vaUo9fgTE1es8o/XjVncf9wGnAWHJdbc1yG337SSEiHouIe5Lbfwd8DzgZuAS4JVntFuA3ktuXAJ+NiGcj4ofAPuCcHOK6C3g67T5JAv4ZcFvW211Ps7hSFHKsGsUVEV+LiOeTP/8GOCWPbbcbFyVfW+t4PbA/IjYyAkDeGh2/QjR5z+gm5wD7IuKhiPh/wGepHLfM9G1SqCVpGjgL+J/AKyLiMahcBMDLk9VOBh6pedghir8gzgMej4gf1LRtkXSvpDslnVdwPADvTb6mubnm43w3HKs17wK+UvN3mcerm6+tyzj6Pxtp57VIAXxN0i5Js0lbo+NXuLr3DCj/eK3J/Vrq+6Qg6e8Bnweujohnmq2a0lZ0ve7bOfqF+xiwOSLOAt4H3Crp+ALj+Y/Aq4DXJrF8ImnvhmOFpHngeWApaSr7eDVS6vGSNAZcDPx50tTovBbp3Ig4G7gIeI+k80uIIVXKe0Y3HK81uV9LfZ0UJI1SOblLEfGFpPlxSa9M7n8l8ETSfgg4tebhpwA/KjDWEeCtwOfW2pKvG55Kbu+i8l3iGUXFFBGPR8ThiDgC/Cde+Mqj1GMFIGkH8CZgJpIvW8s+XnTptUXljfeeiHgcmp7XwkTEj5J/nwC+mMTQ6PgVJu09oxuOV43cr6W+TQrJ9/N/BnwvIv645q7bgR3J7R3AX9S0XybpOElbgNOBu4uKF/g14IGIOLTWIGlS0nBy+7QkpoeKCmjtBZp4C7BWaVPqsZL0BuD9wMURsVrTXurxonuvraM+gTY5r4WQtEnSS9ZuUykc2Evj41dUXKnvGWUfrzrfAU6XtCX5BHgZleOWnbJ70/NagH9C5WPVfcDuZHkjMAF8A/hB8u/Lah4zT+V/lw8CF+UU121UPoI+RyXrX560fxq4sm7dS4H7qVQY3AO8OcfjdUxcwH8Bvpscw9uBVxZ5rJrEtY/K96pr5/XGLjlepV5bDeIcB54CfqamreF5LSim05JztCc5X/NJe8PjV1Bcjd4zSj1eKXG+kUpl1P61Y5fl4mEuzMysqm+/PjIzs/Y5KZiZWZWTgpmZVTkpmJlZlZOCmZlVOSnYQJM0UTP65Y/rRsMcq1v3WkkzdW3vlPQfkttDkm5JhkJI++Upkl4n6Y51YrpS0js2um9mnRgpOwCzMkXlF9CvhcoQycBPIuKPGqx+AZUBC4+RJIEbgVHgt2MDtd4RkduQ32br8ScFsxYkYyiNRcRKg1X+lMqPr94REUeSX+3eLOk7yQB9R41kmXyqOKBkLoikbZ+kVyTj91+TtH1T0nWS7pb0/ZIGRbQB4qRg1ppfo/Ir2zS/BWwDLosXhvSeB/5bRPwi8KvA9cmQDgBEZRydv6AybAKSfgk4EMn4RHVGIuIc4GrgQxnsi1lDTgpmrXkDRw/RXeseKpPX1A6UdgHwAVVm0fsm8CJgc93jPgf88+T2ZdQMhlhnbTDHXVQmEjLLjZOCWWvOofEgdg9Q6Wv4nKTXJG0CLo2I1ybL5oj4Xt3j/hr4+5ImqUwo8wXSPZv8exj3A1rOnBTM1pG80T8QEYcbrRMR3wauBL4kaTPwVeB31qqQJJ2V8pigMmz0H1MZmfOpPOI3a4f/12G2vouAv1pvpYi4I/lf/18Bvw58ELgvSQwHqMz/UO9zVIZDfmdWwZpthEdJNVuHpK9TqSp6rOxYzPLmpGBmZlXuUzAzsyonBTMzq3JSMDOzKicFMzOrclIwM7MqJwUzM6v6/++AxQtVOetFAAAAAElFTkSuQmCC\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,MCA_plot,'ok')\n",
    "plt.gca().invert_xaxis()\n",
    "plt.xlabel('T / Kelvin')\n",
    "plt.ylabel('MC Accepts')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "2534c150",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "200.00  191\n",
      "184.00  203\n",
      "169.28  190\n",
      "155.74  197\n",
      "143.28  199\n",
      "131.82  193\n",
      "121.27  202\n",
      "111.57  189\n",
      "102.64  206\n",
      "94.43  195\n",
      "86.88  199\n",
      "79.93  180\n",
      "73.53  185\n",
      "67.65  190\n",
      "62.24  194\n",
      "57.26  181\n",
      "52.68  193\n",
      "48.46  194\n",
      "44.59  195\n",
      "41.02  191\n",
      "37.74  205\n",
      "34.72  197\n",
      "31.94  201\n",
      "29.39  177\n",
      "27.04  190\n",
      "24.87  183\n",
      "22.88  188\n",
      "21.05  187\n",
      "19.37  175\n",
      "17.82  189\n",
      "16.39  186\n",
      "15.08  187\n",
      "13.88  181\n",
      "12.77  182\n",
      "11.74  171\n",
      "10.80  179\n",
      "9.94  196\n",
      "9.14  176\n",
      "8.41  189\n",
      "7.74  173\n",
      "7.12  172\n",
      "6.55  183\n",
      "6.03  186\n",
      "5.55  162\n",
      "5.10  170\n",
      "4.69  185\n",
      "4.32  189\n",
      "3.97  189\n",
      "3.65  176\n",
      "3.36  191\n",
      "3.09  175\n",
      "2.85  179\n",
      "2.62  182\n",
      "2.41  163\n",
      "2.22  164\n",
      "2.04  163\n",
      "1.88  167\n",
      "1.73  188\n",
      "1.59  170\n",
      "1.46  177\n",
      "1.34  160\n",
      "1.24  172\n",
      "1.14  179\n",
      "1.05  162\n",
      "0.96  157\n",
      "0.89  174\n",
      "0.81  179\n",
      "0.75  157\n",
      "0.69  167\n",
      "0.63  165\n",
      "0.58  192\n",
      "0.54  177\n",
      "0.49  160\n",
      "0.45  172\n",
      "0.42  173\n",
      "0.38  181\n",
      "0.35  172\n",
      "0.33  168\n",
      "0.30  144\n",
      "0.28  175\n",
      "0.25  169\n",
      "0.23  169\n",
      "0.21  158\n",
      "0.20  182\n",
      "0.18  153\n",
      "0.17  168\n",
      "0.15  153\n",
      "0.14  169\n",
      "0.13  165\n",
      "0.12  162\n",
      "0.11  153\n",
      "0.10  170\n",
      "0.09  161\n",
      "0.09  159\n",
      "0.08  154\n",
      "0.07  162\n",
      "0.07  155\n",
      "0.06  170\n",
      "0.06  154\n",
      "0.05  159\n",
      "Runtime:  0.847407\n"
     ]
    }
   ],
   "source": [
    "starttime = time()\n",
    "\n",
    "MCS = 300 # MC steps per temperature\n",
    "J = 100 # meV\n",
    "T = 200 # Kelvin\n",
    "T_p = 0.92\n",
    "kB = 8.617e-2 # meV/K\n",
    "\n",
    "MCA_plot = []\n",
    "T_plot = []\n",
    "E_plot = []\n",
    "\n",
    "# set up the sites\n",
    "lat = lattice(50,50,50) # set up the size of lattice\n",
    "lat.populatesites()     # Assign each site a coordinate\n",
    "lat.initialisespins()   # Assign each site a random spin\n",
    "E_0 = lat.gettotalenergy1(J) # get initial total energy of the system\n",
    "E_t = E_0\n",
    "\n",
    "while (T>0.05):\n",
    "    \n",
    "    MCA = 0    # Initialise MC accepts\n",
    "    \n",
    "    for i in range(MCS):\n",
    "        fs = lat.pickrandomsite() # pick a random site to flip its spin\n",
    "        lat.filpspin(fs) # flip spin\n",
    "        dE = lat.sites[fs].getdeltaexchangeenergy(J)*2\n",
    "\n",
    "        if dE <= 0 or np.exp(-dE/(kB*T)) > rd.random():\n",
    "            E_t += dE\n",
    "            MCA += 1 # MC Accepts\n",
    "            \n",
    "        else:\n",
    "            lat.filpspin(fs) # flip the spin back\n",
    "    \n",
    "    print('%.2f  %d' %(T, MCA))\n",
    "    MCA_plot.append(MCA)\n",
    "    T_plot.append(T)\n",
    "    E_plot.append(E_t)\n",
    "    T = T_p * T\n",
    "    \n",
    "print (\"Runtime:  %.6f\"  %(time()-starttime))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "281b4ce8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAERCAYAAABy/XBZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYbklEQVR4nO3de5BkZXnH8d9vdwUyoDHsLmiE7UEFUwQFZEJKUUsiIlgKXkpFe8OaUNmivGQpQwQyVcZUalIaosmmoqEmFSokOwIaNRjvl3irRCKzCMhVEWeWFZTdQSXuemP3yR99Rnp7+/T07Vy6z/dT1bU9p7vnPJxu+pnnfd7zHkeEAADVs6roAAAAxSABAEBFkQAAoKJIAABQUSQAAKgoEgAAVNTIJQDbV9t+yPbtXT7/tbbvtH2H7Q9kHR8AjAqP2nkAtl8g6SeS/jUiTlrhucdL+qCk34uIH9o+KiIeyiNOACi7kasAIuIrkh5u3mb7abY/bXu77a/a/q3koT+S9L6I+GHyWr78ASAxcgkgxaykt0bEaZIulfT+ZPsJkk6w/d+2b7R9TmERAkDJrCk6gEHZPkLScyV9yPby5kOTf9dIOl7SCyUdI+mrtk+KiB/lHCYAlM7IJwA1qpgfRcQpbR7bKenGiPilpO/avkeNhHBTjvEBQCmN/BBQRDyixpf7ayTJDScnD/+HpDOT7evUGBK6r4g4AaBsRi4B2L5W0tckPcP2TtsXSapLusj2rZLukHR+8vTPSFqyfaekL0r604hYKiJuACibkZsGCgAYjpGrAAAAwzFSTeB169bF5ORk0WEAwEjZvn377ohY37p9pBLA5OSk5ufniw4DAEaK7cV22wsdArJ9ju17bN9r+/IiYwGAqiksAdheLel9ks6VdKKk19s+sah4AKBqiqwATpd0b0TcFxG/kHSdHpu+CQDIWJEJ4CmS7m/6eWey7QC2N9uetz2/a9eu3IIDgHFXZAJwm20HnZQQEbMRMRURU+vXH9TEBgD0qcgEsFPSsU0/HyPpgWHvZG5uTpOTk1q1apUmJyc1Nzc37F0AwEgqchroTZKOt32cpO9JukDSG4a5g7m5OW3evFl79+6VJC0uLmrz5s2SpHq9PsxdAcDIKawCiIhHJb1FjfV67pL0wYi4Y5j7mJ6e/tWX/7K9e/dqenp6mLsBgJFU6HkAEfHJiDghIp4WETPD/v07duzoaXu/GGYCMIrGei2gDRs29LS9H8vDTIuLi4qIXw0zkQQAlN1YJ4CZmRlNTEwcsG1iYkIzM8MrNrIcZqKyAJClsU4A9Xpds7OzqtVqsq1arabZ2dmhNoCzGmaisgCQtZG6HsDU1FSUbTG4yclJLS4evM5SrVbTwsJC6X4vgOqxvT0iplq3j3UFkIeshpmyrCwYVgIgkQAGltUwUxYNbIaVADQjAQxBvV7XwsKC9u/fr4WFhaH0GLKoLIbdsKaaAEYbCaCksqgshjmsRDUBjD6awBUyzMYyTWpgdNAExlCHlYZVTTCMBBSHBFAhwxxWGkaTmmEkoFgkgIoZVsN6GNXEMJrSVBBA/0gA6MswqolBh5GoIIDBkADQt0GriUGHkQatIKgeUHUkABRm0GGkQSoIqgeABIACDTqMNEgFwcWCABIACjbIMNIgFUS/1QPDRhgnJACMrEEqiH6qB4aNMG5IABhp/VYQ/VQP/Q4bUTWgrEgAqKR+qod+ho2oGlBmrAUEdKmf9Y9YMwllwFpAwID6GTbqtWpguAh5IgEAXepn2KiXZnOvw0UkCwwsIkbmdtpppwUwSrZt2xYTExMh6Ve3iYmJ2LZt20HPrdVqBzxv+Var1Qb6vYCk+WjznUoPAMjY3NycpqentWPHDm3YsEEzMzNtq4ZVq1ap3f+PtrV///4DttFbQC/oAQAF6Xaqai/DRb30FhgqQhoSAFASvTSZu00WTENFJyQAoCR6aTJ3myy6PXmNKqGi2jUGynqjCQw8Ztu2bVGr1cJ21Gq1tg1g220by7YP+D00lMebaAID1dNNs5iG8vijCQxUUDdDRd00lBkiGk8kAGCMddNXWKmhTCN5fJEAgDG30jTUlaqElRrJVAejiwQAVNxKVUKnISKqg9FGExhAR52axJJoII+AUjWBbb/G9h2299s+KCgA5dFpiGilBjLDQ+VW1BDQ7ZJeJekrBe0fQJc6DRF1aiAzPFR+hSSAiLgrIu4pYt8AepfWSO5UHdA8Lr/SN4Ftb7Y9b3t+165dRYcDoEmn6oDmcfll1gS2/XlJT2rz0HRE3JA850uSLo2Irjq7NIGB0UHzuDxybwJHxFkRcVKb2w1Z7RNAefTbPGZoKD+lHwICMJr6aR4feeSRDA3lqKhpoK+0vVPScyR9wvZniogDQLZ6bR5Lats43rJlC1VBBoqaBfTRiDgmIg6NiKMj4iVFxAGgGGnVwcMPP9z2+UtLS1QFGeBMYAClkdY4boeGcfdKdSYwALTTbmgozeLiIsNCAyIBACiNdkNDa9eubftc2wwLDYgEAKBUWhvHW7duPagqsK3W4Wuaxb0jAQAotXZVQVrvkmZxb0gAAEqvtSpYPpt4JXv37tWmTZuoCFKQAACMnF6axfv27aMiSEECADByemkWN6MiOBAJAMBI6qZZ3A4VwWNIAADGQmtVsHr16hVfs3fvXm3cuLGy1QAJAMDYaK4Krrnmmp5OKqtiNUACADCWeq0IqtgfIAEAGFu9VgRV6w+QAABUQnNF0I0q9AdIAAAqY7ki2LZtG/0BkQAAVFA//YHp6emcossPCQBAJfXaH1hefnqcKgESAIDK67Y/sLi4qI0bN2rdunVjkQhIAACg3voDS0tLY9EXIAEAQJNuq4Hl8wZGOQmQAACgxXI1sFIS2Ldv30hXAiQAAEjRzbLTo3y+AAkAAFIsDwd1s9T0KJ4vQAIAgA7q9bp2796tbdu2db2e0KgkARIAAHShXq93vZ7QqFQCJAAA6FIvM4S2bNmSU1T9IwEAQA+6PV9gaWmp9FXAmqIDAIBRVK/XJUmbNm3Svn372j5n06ZNBzy3bEgAANCn5S/2jRs3tn18uR/Q/NwyYQgIAAZQr9c7ThMt80qiJAAAGNDWrVs79gMWFxdL2Q9gCAgABtRNP6CMQ0FUAAAwBCudJ1DGqaEkAAAYkuXzBNKUbWooCQAAhqher3c8UaxMVUAhCcD2lbbvtn2b7Y/afmIRcQBAFmZmZlIfK1MVUFQF8DlJJ0XEsyR9S9IVBcUBAEO30tTQsiwYl5oAbE9ltdOI+GxEPJr8eKOkY7LaFwAUYevWramPlWXBOEdE+wfsb0g6QtK1kq6LiDszCcD+T0nXR8S2lMc3S9osSRs2bDhtcXExizAAYOjWrVunpaWl1MfXrl2r3bt3Zx6H7e0RcdAf9akVQEScKullkvZJ+nfbt9i+zHbnZfAe2+Hnbd/e5nZ+03OmJT0qKTUNRsRsRExFxNT69eu72TUAlMJKJ4gV3Q9IrQAOeqJ9sqQLJL1W0vcj4oyBdmxvknSxpBdFxN5uXjM1NRXz8/OD7BYAcjU3N9fxBLE8qoCeK4CWF6+SdJSkoyUdLmnXgMGcI+kySed1++UPAKNo+QSxNEVWAR0rANvPl/R6Sa+QdLuk6yR9OCJ+PNBO7XslHSppeXDsxoi4eKXXUQEAGFWd+gFZVwE9VwC275f0Lkl3STo1Is6OiKsH/fKXpIh4ekQcGxGnJLcVv/wBYJR1mhVUVBXQaRZQLSIWm34+PCL25BZZG1QAAEZZUVVAP7OAFpMXPsf2nWpUArJ9su33ZxIlAIyxslUB3TSB/07SS5SM10fErZJekGFMADCWVjpDOO8Lx3Q1Cygi7m/Z1H4+EwCgo05VwI4dO3KMpLsEcL/t50oK24fYvlTJcBAAoDf1el2HH35428eOPPLIXGPpJgFcLOnNkp4iaaekU5KfAQB9OOyww9pu/9nPfpZrHCsmgIjYHRH1iDg6Io6KiI0Rkb64BQCgo4cffrjt9j179uhNb3pTbnGkXhPY9t93emFE/PHwwwGA8bdhwwalLWx51VVX6Ywzzsjl2sGdKoCLJT1P0gOS5iVtb7kBAPrQ6YIxEZHbbKBOJ4KtlfQaSa9TY8XO69VYBuKHuUTWBieCARgXnU4Ks639+/cPbV/9nAi2FBFXRcSZkt4o6YmS7rD9+0OLCgAqqtN00LxmA63YBLb9bEmXSNoo6VNi+AcABlav13XEEUe0fSyv2UCdFoP7C9vbJb1N0pclTUXERVldGQwAqmbPnvbLq+3ZsyeXZSE69QD2S7pP0k+TTctPtKRILuieK3oAAMbJ5ORk6mygYS4O188FYY6T9CI1Lgv5MkkvT27L9wEAA+g0GyiPxeG6viRkGVABABg3nWYD1Wo1LSwsDLyPgS4JCQDIRqfZQGnDQ8NCAgCAAtXrda1a1f6rePXq1Znuu9MsoFnbr7T9+EwjAICKSzvpa9++bFfe71QBXC3pZEmftP0F25fZPjnTaACggtL+0s+6AkhdDC4ibpR0o6R3JstCnC3pT2w/U9I3JH06Ij6YaXQAUAFpf+kXWQH8SrIsxLURcWFEnCrpfZKOzzQyAKiIWq3WdrvtTKeC9tUEjojtEZE+gRUA0LWZmRnZPmh71iuDMgsIAApWr9eVdk5WllNBSQAAUAKdGr5ZDQN1mgb69qb7r2l57K8yiQYAKqpTwzerYaBOFcAFTfevaHnsnAxiAYDKSmsES9kNA3VKAE653+5nAMAAOi0Ml9X5AJ0SQKTcb/czAGAAnS4Cn9X5AKkngkk62fYjavy1/2vJfSU/H5ZJNABQYatXr277Zd9uiugwdLom8OqIeEJEPD4i1iT3l39+XCbRAECFpf2lHxGZzARiGigAlESnRnAWM4FIAABQEp0awVnMBCIBAEBJ5H1tABIAAJRIntcGKCQB2P5L27fZvsX2Z23/ZhFxAEDZ5HltgKIqgCsj4lkRcYqkj0t6R0FxAECp5HltgEISQEQ80vTj4eLEMgCQlG8F0OlEsEzZnpF0oaQfSzqzqDgAoEzGogKw/Xnbt7e5nS9JETEdEcdKmpP0lg6/Z7Ptedvzu3btyipcACiFPJeFdtpFCPJiuybpExFx0krPnZqaivn5+RyiAoBidFr2oVaraWFhoZ/fuT0iplq3FzULqPl6wudJuruIOACgbPJcFrqoWUDvSoaDbpN0tqQtBcUBAKWS57LQhQ8B9YIhIABV0GkYqJ/v7FINAQEA0uU1FZQEAAAlk9dUUBIAAJRM2hDQsC8MQwIAgJJJG+cfds+WBAAAFUUCAICSSbsmQNr2vvcz1N8GABhY2jUB0rb3iwQAACVDExgAKoomMAAgUyQAACiZtWvXpj42zCWhSQAAUDJbt25NfWzLluGtnUkCAICSqdfrqY8tLS0NbT8kAACoKBIAAJRQHieDkQAAoITyOBmMBAAAJZTHNQFIAABQQnlcE4AEAAAlRAUAABVFBQAAFZXHgnAkAAAooTwWhCMBAEBFkQAAoIQ4EQwAKooTwQCgopgGCgAVxTRQAKgopoECQEUxDRQAkBkSAACUENNAAaCimAYKABVFExgAKoomMAAgMyQAACihsW8C277UdtheV2QcAFA2Y90Etn2spBdL2lFUDABQVuO+FtDfSnq7pOF1NABgTIztWkC2z5P0vYi4tYvnbrY9b3t+165dOUQHAMXLowJYM7Tf1ML25yU9qc1D05L+TNLZ3fyeiJiVNCtJU1NTVAsAKiGPCiCzBBARZ7XbbvuZko6TdGtyQsMxkm62fXpEfD+reABglNhuO+d/mCeCZZYA0kTENyUdtfyz7QVJUxGxO+9YAKCsOBEMAJCZ3CuAVhExWXQMAFA2q1atajvnf2xOBAMAtDfWJ4IBANKxGigAVBRNYABAZkgAAFBCY78aKACgPZrAAFBRNIEBoKJoAgMAMkMCAIASogkMABVFExgAKqrThV/m5uaGsg8SAACUUKcLv0xPTw9lHyQAACihWq2W+tji4uJQ9kECAIASmpmZSX1sWNcFJgEAQAnV6/XUx4Z1XWASAACUVNZnA5MAAKCksj4bmAQAABVFAgCAksr6bGASAACUVNZnA5MAAKCk0s4F6HSOQC9IAABQUjMzM5qYmDhg28TERMdzBHpBAgCAkqrX65qdnVWtVpNt1Wo1zc7OdjxHoBce5sUFsjY1NRXz8/NFhwEAI8X29oiYat1OBQAAFUUCAICKIgEAQEWRAACgokgAAFBRIzULyPYuSf1eCWGdpN1DDGdYiKs3xNWbssYllTe2cYyrFhHrWzeOVAIYhO35dtOgikZcvSGu3pQ1Lqm8sVUpLoaAAKCiSAAAUFFVSgCzRQeQgrh6Q1y9KWtcUnljq0xclekBAAAOVKUKAADQhAQAABU1FgnA9rG2v2j7Ltt32N6SbD/S9udsfzv59zeaXnOF7Xtt32P7JRnGdrXth2zf3rTtetu3JLcF27ck2ydt/7Tpsatyjuudtr/XtP+XNj1W5PG60vbdtm+z/VHbT0y2F328Cv98tcT4jKZjcYvtR2xf0ul9zUvyOf9msv/5ZFvq8cspprTvjcKPV0uc5ySfo3ttXz7UXx4RI3+T9GRJz07uP17StySdKOmvJV2ebL9c0ruT+ydKulXSoZKOk/QdSasziu0Fkp4t6faUx98j6R3J/cm05+URl6R3Srq0zXMLPV6Szpa0Jrn/7qb3sejjVfjnq0O8qyV9X1It7X3NOZ4FSetatrU9fjnGlPa9UfjxankfvyPpqZIOST5XJw7r949FBRARD0bEzcn9/5N0l6SnSDpf0jXJ066R9Irk/vmSrouIn0fEdyXdK+n0jGL7iqSH2z1m25JeK+naLPbdSae42ij0eEXEZyPi0eTHGyUdk8W+e41LJfh8dfAiSd+JiH7PnM9D2vHLRYfvjTI5XdK9EXFfRPxC0nVqHLehGIsE0Mz2pKRTJf2vpKMj4kGp8WZLOip52lMk3d/0sp0q5o1/vqQfRMS3m7YdZ/sbtr9s+/kFxPSWZKjl6qaSvCzHS5L+UNKnmn4u8niV+fN1gQ78w6Ld+5qnkPRZ29ttb062pR2/3LV8b0jFH69lmX6WxioB2D5C0oclXRIRj3R6apttRcyHfb0O/J/0QUkbIuJUSW+T9AHbT8gxnn+U9DRJpySxvCfZXorjZXta0qOS5pJNRR+vNIUeL9uHSDpP0oeSTWnva57OiIhnSzpX0pttv6CAGNpq871RhuO1LNPP0tgkANuPU+NNnIuIjySbf2D7ycnjT5b0ULJ9p6Rjm15+jKQH8oo1iWeNpFdJun55WzJksJTc367G2N8JecUUET+IiH0RsV/SP+mxYYsyHK9Nkl4mqR7J4GjRx0vl/XydK+nmiPiB1PF9zU1EPJD8+5CkjyYxpB2/3LT73ijD8WqS6WdpLBJAMpb+z5Luioj3Nj30MUmbkvubJN3QtP0C24faPk7S8ZK+nle8ibMk3R0RO5c32F5ve3Vy/6lJXPflFdDy/4yJV0panvFS6PGyfY6kyySdFxF7m7YXerxU3s/XAZVlh/c1F7YPt/345ftqNPVvV/rxyyuutt8bRR+vFjdJOt72cUlld4Eax204iu5yD6lT/jw1yqLbJN2S3F4qaa2kL0j6dvLvkU2vmVbjL8Z7JJ2bYWzXqlFG/lKNbH5Rsv1fJF3c8txXS7pDjU7/zZJenmdckv5N0jeT4/gxSU8uw/FSo4l6f9N7e1VJjlfhn682cU5IWpL0603bUt/XnGJ6avIe3Zq8X9PJ9tTjl1Ncad8bhR6vNnG+VI0ZSt9ZPnbDurEUBABU1FgMAQEAekcCAICKIgEAQEWRAACgokgAAFBRJABUhu21TSs8fr9lxcdDWp57he16y7Y32v6H5P4q29ckSwW0O1tTtl9o++MrxHSx7QsH/W8D+rGm6ACAvETjrOFTpMaSv5J+EhF/k/L0s9VYqO8gyRf+VZIeJ+kPYoC51BGR2RLWwEqoAIAWyXpCh0TErpSnbFXjJKYLI2J/cqbr1bZvShamO2C1xqRaWHByHYNk2722j07Wnr802fYl2++2/XXb3ypoMUBUCAkAONhZapyZ2s4bJJ0m6YJ4bInqaUn/FRG/I+lMSVcmSx5IkqKxpswNaiwrINu/K2khkrV6WqyJiNMlXSLpz4fw3wKkIgEABztHBy453exmNS6y0rxA2NmSLnfjym5fknSYpA0tr7te0uuS+xeoaRHAFssLGW5X44I3QGZIAMDBTlf64m13q9EbuN72byfbLOnVEXFKctsQEXe1vO5rkp5ue70aFz75iNr7efLvPtGjQ8ZIAECT5Ev97ojYl/aciPgfSRdL+oTtDZI+I+mty7OBbJ/a5jWhxjLI71Vj9cmlLOIHesFfGMCBzpX06ZWeFBEfT/6a/7SkF0t6h6TbkiSwoMa1C1pdr8byvm8cVrDAIFgNFGhi+3NqzO55sOhYgKyRAACgougBAEBFkQAAoKJIAABQUSQAAKgoEgAAVBQJAAAq6v8BmQZ2C6++8TQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAgnklEQVR4nO3dfZRkdX3n8fene6YndAIJ0zNyEOhqQDA7JAZhJJqoMUFFXCMRXTNuSxA82zoQVzYnMc72HjVHe48Pa3Z9yOhpdZTsVHjICW5IgkZxJazHB5whKDM8yCA9wwiBgcmu6CAD09/9o27V1PTUvfV066nr8zqnzlT/6t663751p759f4+KCMzMzABGeh2AmZn1DycFMzOrcFIwM7MKJwUzM6twUjAzs4oVvQ6gHWvWrImpqaleh2FmNlC2b9/+WESsrfXaQCeFqakptm3b1uswzMwGiqTdaa+5+sjMzCqcFMzMrMJJwczMKpwUzMyswknBzMwqnBTM2lQsFpmammJkZISpqSmKxWKvQzJr2UB3STXrtWKxyMzMDAcOHABg9+7dzMzMADA9Pd3L0Mxa4jsFszbMzs5WEkLZgQMHmJ2d7VFEZu1xUjBrw549e5oqt85yVV77nBTM2jA5OdlUuXVOuSpv9+7dRESlKs+JoTlOCmZtmJubY3x8/Iiy8fFx5ubmehTR8HJVXj6cFMzaMD09zfz8PIVCAUkUCgXm5+fdyNwDrsrLhwZ5jeb169eHJ8QzMyhNkLl799HzvBUKBRYWFrofUB+TtD0i1td6zXcKZrYsuCovH04KZrYsuCovH04KHeBucWa9MT09zcLCAouLiywsLDghtMAjmnPmEa5mNsh8p5Azd4szs0HmpJAzd4vrLVfdmbXHSSFnHuHaOx7RatY+J4WcuVtc77jqzqx9Tgo5c7e43nHVndXj6sX6PKLZlg2PaLUsS3sGQukufhj/aPOIZhsKrrqzLK5ebEzHkoKkUyR9XdLdknZKemdSvlrSVyXdl/x7fNU+myTtknSvpAs6FZstT666syyuXmxMx6qPJJ0InBgRt0s6FtgO/B7wFmB/RHxQ0ruB4yPiTyWtA64BzgOeDdwMnBkRh9KO4eojM2uUqxcP60n1UUQ8HBG3J8+fAO4GTgIuAq5ONruaUqIgKb82Ip6KiAeAXZQShJlZ21y92JiutClImgKeD3wHOCEiHoZS4gCelWx2EvBg1W57k7Kl7zUjaZukbfv27eto3Ga2fLh6sTEdTwqSfgH4G+CqiPhx1qY1yo6q24qI+YhYHxHr165dm1eYlgN39+sMn9f8eMK8+jo6IZ6klZQSQjEibkiKH5F0YkQ8nLQ7PJqU7wVOqdr9ZOChTsZn+fFEgJ3h82rd1smGZlFqM9gfEVdVlX8EeLyqoXl1RLxL0lnAX3G4oflrwBluaB4MbsTrDJ9X64SshuZO3in8JnAJcKekO5Ky/wx8ELhe0luBPcC/A4iInZKuB+4CngGuzEoI1l/c3a8zfF6t2zqWFCLiG9RuJwA4P2WfOcBdAQbQ5ORkzb9oPRFge3xerds8otly4e5+neHzat3mpDDk8urZ4u5+neHz2h733GqeJ8QbYp4gzJYzX9/pshqanRSGmHu22HLm6zudZ0m1mtyzxZazfr6++7lay0lhiHnpUFvO+vX67vdlY50Uhph7tthy1q/Xd7+v6+CkMMTcs6U3+rnqoB+1er769fru52otACJiYB/nnntumA2SrVu3xvj4eFCa7DGAGB8fj61bt/Y6tL60HM9XoVA44vcpPwqFQtdiALZFyveq7xTMuqjfqw76zXI8X/1arVXmpGDWRX1fddBnluP56tdqrTInBbMu6tceMc3qRLtIrfdcLudrqb5e1yGtXmkQHm5TsEGzHOrIO/E7pL3nxo0bB/589SMy2hR6/sXezsNJwQbR1q1bo1AohKQoFAoD9wXXiYbSrPcc9PPVj5wUbGD4C6D/Sar5BS6pof1rfcbtvqc1JyspuE3B+ka/j/S0knbq+dM+49WrV7f8npYvJwXrG8ux++Fy1E6XyrTPuPwerbyn5ctJwfrGcux+uBy106Uy7bPcv39/X3fT7BfdGA3vqbOtb3iq4+XPn3Hr8lwfwlNn20Do95Ge1j5/xq3rVvXqUCYFT0jWn/p9pKe1L8/PeBj+H1f/jrXusKAD1atp3ZIG4dFKl9TlMHjIbNgNw//jWr9jrUcr40PI6JI6dG0KrtM0G3zD8P847XesNlBtCpK2SHpU0o6qsl+T9C1Jd0r6O0nHVb22SdIuSfdKuqBTcQ1yD5dhuF22/tDv11qn/h/30++d9bt0tHo17Rai3QfwUuAcYEdV2XeB30qeXw68P3m+DvgesAo4FbgfGK13jFaqj/phLvNWDMPtsvWHQbjWOvH/uN9+705+V9GraS6AqSVJ4ccc7gZ7CnBX8nwTsKlqu38EXlTv/YepTWFQk5kNnkG41jrx/7jffu9Oflf1U1L4JnBR8vyPgCeS558E3ly13eeAN6S85wywDdg2OTnZ0gkZxPl1PDdM5w3iddEJg3Kt5f15tfN7txpL2n5bt26NiYmJI+KYmJjI7Zrsp6Twy8BXgO3Ae4HHk/K/qJEUXl/v/YdpQrx++ytmuRnUO8hOGNZrrdXfu9VrJ2u68JUrVx4Vx9jY2PJLCkteOxO4LXneteqjQeUvrc4a1i/CWob1Wmv192712knbb3R0tGZ5ntdj3yQF4FnJvyPAXwKXJz+fxZENzT+kQw3Ng8zVG50zKFUm3TKs11orv3er107aflmPvK7HniQF4BrgYeBpYC/wVuCdwA+SxwdJGp2T7Wcp9Tq6F7iwkWMMW1KwzvGdgmXJShbdvFMov2+7SbpndwqdfjgpWF6GtcrE6qt3bXSrTSHPa9NJwawBw1plYtkauRPoRu+jPO9is5LC0E1zYWbWjJGREdK+JwuFAnNzc0eNKi4Wi8zOzrJnzx4mJydrbtPu8SWxuLjY0nt66mwzsxZlLQlaa8nYvJeVbWf501Y4KZiZZai1BkS1pWsa5L3uQbfXoHBSMDPLUL0GRJry5HXFYrGtdQ9qTcjX7XVG3KZgZtagrCm75+bmjlouc+k2WdN657ncZj1uUzAzy0FWVU6taqOl22Tp1nKb9TgpmA2ZflozYNBkVeVkVQ818td+2v67d+/u6mfk6iOzIdLNKoph0+5qcFkrreX9Gbn6yMyA/qmiWI7a7SWU1cupm5+Rk4LZEBnk5Wj7Xbu9hMr7p+nWZ+TqI7MhMgwL3g+6bnxGrj4yM6D7A6Gseb3+jJwUzIZItwdCWfN6/Rk5KZgNmenpaRYWFlhcXGRhYcEJIUW3u+5WH292dpa5ubmefEYr6m0g6eeBJyNiUdKZlNZZ/lJEPN3x6MzMemBp193ypHZAR76gu328LHUbmiVtB14CHA98G9gGHIiInv954YZmM+uEbjfId/t47TY0KyIOABcDn4iI1wHr8gzQzKyf5N11t15VVNbxlu57xRVXdLZaK231nfID+GfgRZTuEs5Kyu6st183Hl55zcw6Ic81uxtZrjPteBMTE0ftu/TRytKcZKy81sidwjuBTcAXI2KnpNOAr7efjszM+lOe3UIbGUWedrzytlnyHu3cSFI4ISJeGxEfAoiIHwL/J7cIzGzZGfRJ9/LqFtro+gppx9u/f39Dx8l1tHPaLUT5AdzeSFkvHq4+Mus/jVSXDINa54Emq6LSqpVaea9qtFJ9JOlCSZ8ATpL08arHF4Bn8ktLZraceNK9knbXV4D6S4ECSMp1tHNW9dFDlLqf/gzYXvW4Ebig3htL2iLpUUk7qsrOlvRtSXdI2ibpvKrXNknaJeleSXXf38yO1C9VNp50r6Td9RWgsaVAIyLfsQxptxDlB7ASGAOeB/wqMFZvn2S/lwLnADuqyr4CXJg8fzVwS/J8HfA9YBVwKnA/MFrvGK4+MivppyqbPHvuDLK8z0Oe70ebvY9ekXxJfxz4JLBL0oUNJJtbgaWtJAEclzz/RUp3IwAXAddGxFMR8QCwCzgPM2tIP1XZ9HpCt36R93no2nlNyxZx+K/7e4DnVP18OnBPvf2Sbac48k7h3wB7gAeBHwGFpPyTwJurtvsc8IaU95yhVK21bXJysukMabYcSar5V6SknsSzdevWKBQKISkKhcLQNTKX1ToPzZ6b8vZLP9uRkZHYuHFjS3GRcafQyBf7rUt+1tKyjH2XJoWPA69Pnr8RuDl5/hc1ksLr672/q4/MSlxlMxiarear14NpbGyspYSblRQaqT7aKekmSW+RdCnwd8B3JV0s6eIG9q92KXBD8vyvOVxFtBc4pWq7kzlctWRmdbjKZjA0W82X1YMJ4ODBg7lXETaSFH4OeAT4LeBlwD5gNfC7wGuaPN5DyfsA/A5wX/L8RmCDpFWSTgXOAG5r8r3Nhlav5+C3xjTbM6uRHlt59+qqO3V2RFzWyhtLuoZSElkjaS/wXuA/AB+TtIJSV9eZ5Bg7JV0P3EVpDMSVEXGoleOaDavp6WkngT43OTlZc4Tz5ORkU9s3sm+r6t4pSDpT0tfK4w0kPU/Sf6m3X0S8KSJOjIiVEXFyRHwuIr4REedGxK9FxK9HxPaq7eci4vSIeG5EfKm9X8vMrP80W81Xb/Da2NhY7lWEjVQffYbShHhPA0TE94ENuUZhZjYEmq3mWzp4TVLltYmJCbZs2ZL73WEjSWE8IpbW73uaCzMbeL0YBd7scqjl7bdu3crq1as7Hl/dNgXgMUmnU+oChaQ3AA93NCozsw7rpyUw6ykWi1x22WU8/fThVZAff/xxLr/8ciDfeBtZjvM0YB74DeBfgQcojSlYyC2KFnk5TjNrVbeXwGxHWqzQWrxtLccZET+MiJcDa4FfjogX90NCMDNrR68n7qtXdVX9elYPpLzjbaT30X+V9EsR8dOIeELS8ZI+kGsUZmZdltUNtNPKVVe7d+8mIipVV+XEsPT1LF3vkkppVtP/W/4hIv6V0gynZmYDq5ejwOuNbK43krmsV11SRyWtKv8g6RhKU1ybmQ2sPEaBt9p7qV7VVaNVQscee2xjgTYjbVKkODw53buAbwBvBS5Pnr+r3n7deHhCPDPrlXbWsKg3gWHWMpxjY2Ntr5tBOxPiRcSHgQ9Qmvb6LOD9SZmZ2dBqZw2LelVXc3NzRwxUq3bw4MGWjtmoRrqkngo8HBE/S34+Bjgh+qAHkrukmlmvjIyM1GwElsTi4mLd/YvFIrOzs+zZs4fJyUnm5uaOqLpKSwq1NHrMqu1b75JKaYrr6qMdSsrMzIZWu72X6o1szlqXudVjNqKRpLAiIir3K8nzsdwiMDMbQJ3uvVRvMrwySbn2QGokKeyT9NqqAC4CHsstAjOzAZTXGhZpPZiWvn+aiOj6NBenA0Xg2ZSW4nwQuCQi7s8tiha5TcHMBtnS+ZegdLdRK7msWbOGxx9//Kj3mJiY4LHHmvs7vd1pLu6PiBcC64B1EfEblFZeMzOzNrTTg6lTGpkltWyS0pKZG4AfAzWzjJmZNaaZ+Zf2799fc9u08lZlJgVJBeBNyeMZoACs74fuqGZmgy5tuc1a6yasXr26ZvVR1+Y+kvRN4CZgJfCGiDgXeMIJwcwsH3Nzc6xcufKo8ieeeOKIKTOKxWLNhADwnOc8J9eYstoU9gHHAidQmjYbkoV2zMysfdPT0xx33HFHlR88ePCIdoWsNoZbbrkl15hSk0JEXAT8KnA78GeSHgCOl3RerhGYmQ2BtK6naW0C1e0KWRPkHTp0KNc4M9sUIuL/AVuALZKeBfw+8D8knRIRp+QaiZnZMpW19Gdau0J1W0HaNgCjo6O5xtrI4DUAIuLRiPhE0iX1xblGYWa2jGV1PW1kZPTc3BwjI7W/rsvJJS8NJ4VqEZG+NlxC0hZJj0raUVV2naQ7kseCpDuqXtskaZekeyVd0EpcZmb9KKvr6dKRyxMTExxzzDFccskllWqm6elp3va2tx21fzOT5jUsbU7tdh/AS4FzgB0pr38UeE/yfB3wPUqL95wK3A+M1juG11Mws0FQb/2EsrQ1GjZu3HhUefVj48aNTcVDO+sptJFsbgVqtqColN7eCFyTFF0EXBsRT0XEA8AuwA3aZrYsNDp5Xlo10/z8fObynPPz87nFmjVO4cOS3l6j/D9J+lCbx30J8EhE3Jf8fBKlOZXK9iZlteKakbRN0rZ9+/a1GYaZWec1OnleWjVTvR5GefZAyrpTeA1QK/18DPi3bR73TRy+S4DSRHtL1RwTERHzEbE+ItavXbu21iZmZn2n3voJkD46Oe8eRlmykkJExFFL+SRlLbduSFoBXAxcV1W8F6ju4noy8FCrxzAzG0Rp1UwzMzMNra2Qh6ykcEDSGUsLk7In2zjmy4F7ImJvVdmNlCbbW5Us/3kGcFsbxzAzGzhp1UybN2/Otd0gS9bgtfcAX5L0AWB7UrYe2ARcVe+NJV0DvAxYI2kv8N6I+BywgSOrjoiInZKuB+6iNPHelRGR7zA9M7MBMD09neuiOc3KXGRH0q8AfwL8SlK0E/hIRNzZhdjq8iI7ZjYspqamUkc1r1q1ip/97GcNv1fWIjv1prnYAVza8JHMzKwjsuY/OnjwYG7HSU0Kkm7M2jEiXpv1upmZ5aNYLDIyMpLa9TTPNRWy7hReRGnswDXAd2ijx5GZmbWmPJleWkJYsWLFUYPg2pHapiBpFHgFpTEFzwP+AbgmInbmdvQ2uU3BzJa7rLYEgImJCR577LGm3jOrTSFrPYVDEfHliLgUeCGlqSdukfSOpo5uZmYty2pLgPzXaM6c+ygZN3AxsBW4Evg4cEOuEZiZWap67QXdXKP5auCblGY6/bOIeEFEvD8ifpRrBGZmlqrWKOeyWpPqtSurofkS4KfAmcB/rJq3W5SmwDh6YVEzM8tVeSDb7OzsUW0LaQvvtCM1KUREx6bVNjOzxpUTw2WXXcbTTz9dKf/JT37C5ZdffsQ27fIXv5nZAJidnT0iIZQdPHiQ2dnZ3I7jpGBm1ueKxWJmt9R6PZSa4aRgZtbHyoPXsuTZA8lJwcysj9VaorPa2NhYrj2QnBTMzPpYvdHMW7ZsyXWq7cxZUs3MrHeKxWLqa61Mb9EI3ymYmfWpPHsVNcpJwcysT2X1Ksp7zqMyJwUzsz5SLBaZmppiZGQkc8RyRLBixQquuOKKXI/vNgUzsz5R7n5a7m2UtoZC2aFDh/jUpz4FwObNm3OJIXON5n7n9RTMbDmpt3ZCmtHRUZ555pmGt29pPQUzM+uuVkcm17ujaIaTgplZn0gbmTw6Opq5X56zpTopmJn1iVprJ4yPjzMzM5O6pgLAMccck1sMHUsKkrZIelTSjiXl75B0r6Sdkj5cVb5J0q7ktQs6FZeZWb+anp5mfn6eQqGAJAqFAvPz82zevJlLL700db+saTCa1cneR18APgn8ZblA0m8DFwHPi4inJD0rKV8HbADOAp4N3CzpzIjIr6LMzGwATE9PHzVtRbFY5LOf/WzqPgMxIV5E3AosHV2xEfhgRDyVbPNoUn4RcG1EPBURDwC7gPM6FZuZ2SBJW0sBQNJAT4h3JvASSd+R9E+SXpCUnwQ8WLXd3qTsKJJmJG2TtG3fvn0dDtfMrPeyeiVFRK4T4nU7KawAjgdeCPwJcL1Kiz+rxrY1B1BExHxErI+I9WvXru1cpGZmfaJe9VCeo5q7nRT2AjdEyW3AIrAmKT+laruTgYe6HJuZWV+am5tj5cqVqa/Pz8/ndqxuJ4X/BfwOgKQzgTHgMeBGYIOkVZJOBc4AbutybGZmfWl6eprPf/7zqa8PxOA1SdcA3wKeK2mvpLcCW4DTkm6q1wKXJncNO4HrgbuALwNXuueRmdlh09PTqYPU8hy81rEuqRHxppSX3pyy/RyQXxO6mdkyc8wxx/DTn/60ZnlePKLZzGxA1EoIkO/gNScFM7MBkNXDaCAGr5mZWX6yehgN8uA1MzNrQVYPo0EevGZmZi1Imz673rTazXJSMDMbADMzM02Vt8prNJuZDYDyGszz8/McOnSI0dFRZmZmclubucx3CmZmA2Lz5s1cffXVFAoFFhcXuemmmygWi7kew3cKZmYDolgsMjMzUxmXsHv37kr1UV6Nzb5TMDMbELOzs0cNVDtw4ACzs7O5HcNJwcxsQKStq5C13kKznBTMzAbE+Ph4U+WtcFIwMxsQTz75ZFPlrXBSMDMbEIuLi02Vt8JJwczMKpwUzMyswknBzMwqnBTMzKzCScHMbECcf/75TZW3wknBzGxA3HzzzUclgPPPP5+bb745t2N47iMzswGSZwKoxXcKZmZW4aRgZmYVHUsKkrZIelTSjqqy90n6kaQ7kserq17bJGmXpHslXdCpuMzMLF0n7xS+ALyqRvl/j4izk8dNAJLWARuAs5J9NkvKd+FRMzOrq2NJISJuBfY3uPlFwLUR8VREPADsAs7rVGxmZlZbL9oU/lDS95PqpeOTspOAB6u22ZuUHUXSjKRtkrbt27ev07GamQ2VbieFTwGnA2cDDwMfTcpVY9uo9QYRMR8R6yNi/dq1azsSpJnZsOpqUoiIRyLiUEQsAp/hcBXRXuCUqk1PBh7qZmxmZtblpCDpxKofXweUeybdCGyQtErSqcAZwG3djM3MzDo4olnSNcDLgDWS9gLvBV4m6WxKVUMLwNsAImKnpOuBu4BngCsj4lCnYjMzs9oUUbPqfiCsX78+tm3b1uswzMwGiqTtEbG+1mse0WxmZhVOCmZmVuGkYGZmFU4KZmZW4aRgZmYVTgpmZlbhpGBmZhVOCmZmVuGkYGZmFU4KZmZW4aRgZjZAisUiU1NTjIyMMDU1RbFYzPX9OzYhnpmZ5atYLDIzM8OBAwcA2L17NzMzMwBMT0/ncgzfKZiZDYjZ2dlKQig7cOAAs7OzuR3DScHMbEDs2bOnqfJWOCmYmQ2IycnJpspb4aRgZjYg5ubmGB8fP6JsfHycubm53I7hpGBmNiCmp6eZn5+nUCggiUKhwPz8fG6NzOCV18zMho5XXjMzs4Y4KZiZWYWTgpmZVTgpmJlZhZOCmZlVDHTvI0n7gN1tvMUa4LGcwsmT42qO42qO42rOcoyrEBFra70w0EmhXZK2pXXL6iXH1RzH1RzH1Zxhi8vVR2ZmVuGkYGZmFcOeFOZ7HUAKx9Ucx9Ucx9WcoYprqNsUzMzsSMN+p2BmZlWcFMzMrGLZJgVJp0j6uqS7Je2U9M6kfLWkr0q6L/n3+Kp9NknaJeleSRd0KK4tkh6VtKOq7DpJdySPBUl3JOVTkp6seu3TnYgpI673SfpR1fFfXfVax89VRlwfkXSPpO9L+qKkX0rKe32+enpt1YjxuVXn4g5JP5Z0Vdbn2i3JdX5ncvxtSVnq+etSTGnfGT0/X0vifFVyHe2S9O7cDxARy/IBnAickzw/FvgBsA74MPDupPzdwIeS5+uA7wGrgFOB+4HRDsT1UuAcYEfK6x8F3pM8n0rbrhtxAe8D/rjGtl05VxlxvRJYkTz/UNVn2Ovz1dNrq068o8C/AIW0z7XL8SwAa5aU1Tx/XYwp7Tuj5+dryed4P3AaMJZcV+vyPMayvVOIiIcj4vbk+RPA3cBJwEXA1clmVwO/lzy/CLg2Ip6KiAeAXcB5HYjrVmB/rdckCXgjcE3ex60nK64aunKu0uKKiK9ExDPJj98GTu7EsZuNix5fW3WcD9wfEe3MANBpaeevKzK+M/rJecCuiPhhRBwErqV03nKzbJNCNUlTwPOB7wAnRMTDULoIgGclm50EPFi12166f0G8BHgkIu6rKjtV0j9L+idJL+lyPAB/mFTTbKm6ne+Hc1V2OfClqp97eb76+drawJF/bNT6XLspgK9I2i5pJilLO39dt+Q7A3p/vso6fi0t+6Qg6ReAvwGuiogfZ21ao6zb/XXfxJH/cR8GJiPi+cAfAX8l6bguxvMp4HTg7CSWjybl/XCukDQLPAMUk6Jen680PT1fksaA1wJ/nRSlfa7d9JsRcQ5wIXClpJf2IIaaanxn9MP5Kuv4tbSsk4KklZQ+3GJE3JAUPyLpxOT1E4FHk/K9wClVu58MPNTFWFcAFwPXlcuS6obHk+fbKdUlntmtmCLikYg4FBGLwGc4XOXR03MFIOlS4DXAdCSVrb0+X/TptUXpi/f2iHgEMj/XromIh5J/HwW+mMSQdv66ptZ3Rj+cryodv5aWbVJI6uc/B9wdEX9e9dKNwKXJ80uBv60q3yBplaRTgTOA27oVL/By4J6I2FsukLRW0mjy/LQkph92K6Dyf9DE64ByT5uenitJrwL+FHhtRByoKu/p+aJ/r60j7kAzPteukPTzko4tP6fUcWAH6eevW3HV/M7o9fla4rvAGZJOTe4AN1A6b/npdWt6px7AiyndVn0fuCN5vBqYAL4G3Jf8u7pqn1lKf13eC1zYobiuoXQL+jSlrP/WpPwLwNuXbPt6YCelHga3A7/bwfN1VFzA/wTuTM7hjcCJ3TxXGXHtolSvWv5cP90n56un11ZKnOPA48AvVpWlfq5dium05DP6XvJ5zSblqeevS3GlfWf09HzViPPVlHpG3V8+d3k+PM2FmZlVLNvqIzMza56TgpmZVTgpmJlZhZOCmZlVOCmYmVmFk4INNUkTVbNf/suS2TDHlmy7SdL0krK3SPpk8nxE0tXJVAi1Rp4i6WWS/r5OTG+X9Aft/m5mrVjR6wDMeilKI6DPhtIUycBPIuK/pWz+SkoTFh4lSQKfBlYCl0Ubfb0jomNTfpvV4zsFswYkcyiNRcS+lE0+Rmnw1R9ExGIyaneLpO8mE/QdMZNlclexoGQtiKRsl6QTkvn7/zgpu0XShyTdJukHPZoU0YaIk4JZY15OaZRtLf8eOBfYEIen9J4F/ndEvAD4beAjyZQOAERpHp2/pTRtApJ+HViIZH6iJVZExHnAVcB7c/hdzFI5KZg15lUcOUV3tdspLV5TPVHaK4F3q7SK3i3AzwGTS/a7Dvj95PkGqiZDXKI8meN2SgsJmXWMk4JZY84jfRK7eyi1NVwn6aykTMDrI+Ls5DEZEXcv2e9bwHMkraW0oMwN1PZU8u8h3A5oHeakYFZH8kV/T0QcStsmIr4JvB34B0mTwD8C7yj3QpL0/Br7BKVpo/+c0sycj3cifrNm+K8Os/ouBL5cb6OI+Pvkr/4vA68A3gN8P0kMC5TWf1jqOkrTIb8lr2DN2uFZUs3qkPRVSr2KHu51LGad5qRgZmYVblMwM7MKJwUzM6twUjAzswonBTMzq3BSMDOzCicFMzOr+P9xUw0Yrya+nQAAAABJRU5ErkJggg==\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,E_plot,'ok')\n",
    "plt.gca().invert_xaxis()\n",
    "plt.xlabel('T / Kelvin')\n",
    "plt.ylabel('E / MeV')\n",
    "plt.show()\n",
    "plt.close('all')\n",
    "plt.figure()\n",
    "plt.plot(T_plot,MCA_plot,'ok')\n",
    "plt.gca().invert_xaxis()\n",
    "plt.xlabel('T / Kelvin')\n",
    "plt.ylabel('MC Accepts')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "8c9f0831",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "200.00  206\n",
      "194.00  194\n",
      "188.18  203\n",
      "182.53  188\n",
      "177.06  191\n",
      "171.75  201\n",
      "166.59  184\n",
      "161.60  210\n",
      "156.75  196\n",
      "152.05  191\n",
      "147.48  179\n",
      "143.06  204\n",
      "138.77  194\n",
      "134.61  186\n",
      "130.57  201\n",
      "126.65  193\n",
      "122.85  195\n",
      "119.17  181\n",
      "115.59  193\n",
      "112.12  189\n",
      "108.76  185\n",
      "105.50  186\n",
      "102.33  193\n",
      "99.26  194\n",
      "96.28  175\n",
      "93.39  173\n",
      "90.59  180\n",
      "87.88  190\n",
      "85.24  195\n",
      "82.68  168\n",
      "80.20  190\n",
      "77.80  198\n",
      "75.46  172\n",
      "73.20  187\n",
      "71.00  178\n",
      "68.87  196\n",
      "66.81  173\n",
      "64.80  180\n",
      "62.86  169\n",
      "60.97  186\n",
      "59.14  188\n",
      "57.37  175\n",
      "55.65  173\n",
      "53.98  180\n",
      "52.36  189\n",
      "50.79  180\n",
      "49.26  171\n",
      "47.79  170\n",
      "46.35  176\n",
      "44.96  169\n",
      "43.61  174\n",
      "42.30  179\n",
      "41.04  165\n",
      "39.80  171\n",
      "38.61  182\n",
      "37.45  159\n",
      "36.33  169\n",
      "35.24  171\n",
      "34.18  163\n",
      "33.16  176\n",
      "32.16  175\n",
      "31.20  182\n",
      "30.26  173\n",
      "29.35  162\n",
      "28.47  171\n",
      "27.62  160\n",
      "26.79  166\n",
      "25.99  177\n",
      "25.21  171\n",
      "24.45  157\n",
      "23.72  164\n",
      "23.01  174\n",
      "22.31  157\n",
      "21.65  168\n",
      "21.00  163\n",
      "20.37  168\n",
      "19.76  159\n",
      "19.16  160\n",
      "18.59  173\n",
      "18.03  174\n",
      "17.49  182\n",
      "16.96  177\n",
      "16.46  167\n",
      "15.96  172\n",
      "15.48  168\n",
      "15.02  155\n",
      "14.57  179\n",
      "14.13  165\n",
      "13.71  177\n",
      "13.30  165\n",
      "12.90  159\n",
      "12.51  164\n",
      "12.13  158\n",
      "11.77  163\n",
      "11.42  163\n",
      "11.08  165\n",
      "10.74  165\n",
      "10.42  152\n",
      "10.11  163\n",
      "9.80  159\n",
      "9.51  163\n",
      "9.23  169\n",
      "8.95  152\n",
      "8.68  157\n",
      "8.42  152\n",
      "8.17  159\n",
      "7.92  151\n",
      "7.68  164\n",
      "7.45  150\n",
      "7.23  167\n",
      "7.01  163\n",
      "6.80  143\n",
      "6.60  139\n",
      "6.40  163\n",
      "6.21  151\n",
      "6.02  154\n",
      "5.84  162\n",
      "5.67  160\n",
      "5.50  164\n",
      "5.33  160\n",
      "5.17  151\n",
      "5.02  147\n",
      "4.87  162\n",
      "4.72  176\n",
      "4.58  156\n",
      "4.44  159\n",
      "4.31  143\n",
      "4.18  139\n",
      "4.05  157\n",
      "3.93  137\n",
      "3.81  145\n",
      "3.70  150\n",
      "3.59  147\n",
      "3.48  142\n",
      "3.38  150\n",
      "3.28  139\n",
      "3.18  157\n",
      "3.08  153\n",
      "2.99  148\n",
      "2.90  157\n",
      "2.81  151\n",
      "2.73  134\n",
      "2.65  147\n",
      "2.57  144\n",
      "2.49  130\n",
      "2.42  147\n",
      "2.34  143\n",
      "2.27  132\n",
      "2.20  154\n",
      "2.14  149\n",
      "2.07  160\n",
      "2.01  144\n",
      "1.95  132\n",
      "1.89  136\n",
      "1.84  143\n",
      "1.78  148\n",
      "1.73  136\n",
      "1.68  153\n",
      "1.63  140\n",
      "1.58  143\n",
      "1.53  123\n",
      "1.48  154\n",
      "1.44  126\n",
      "1.40  143\n",
      "1.35  140\n",
      "1.31  138\n",
      "1.27  139\n",
      "1.24  148\n",
      "1.20  139\n",
      "1.16  136\n",
      "1.13  141\n",
      "1.09  127\n",
      "1.06  162\n",
      "1.03  140\n",
      "1.00  134\n",
      "0.97  126\n",
      "0.94  156\n",
      "0.91  160\n",
      "0.88  139\n",
      "0.86  142\n",
      "0.83  128\n",
      "0.81  129\n",
      "0.78  134\n",
      "0.76  136\n",
      "0.74  137\n",
      "0.71  135\n",
      "0.69  128\n",
      "0.67  132\n",
      "0.65  132\n",
      "0.63  114\n",
      "0.61  153\n",
      "0.59  144\n",
      "0.58  139\n",
      "0.56  135\n",
      "0.54  127\n",
      "0.53  144\n",
      "0.51  130\n",
      "0.50  136\n",
      "0.48  143\n",
      "0.47  138\n",
      "0.45  130\n",
      "0.44  125\n",
      "0.43  127\n",
      "0.41  142\n",
      "0.40  126\n",
      "0.39  133\n",
      "0.38  125\n",
      "0.37  133\n",
      "0.35  137\n",
      "0.34  120\n",
      "0.33  124\n",
      "0.32  122\n",
      "0.31  138\n",
      "0.30  116\n",
      "0.30  124\n",
      "0.29  132\n",
      "0.28  128\n",
      "0.27  137\n",
      "0.26  131\n",
      "0.25  143\n",
      "0.25  121\n",
      "0.24  130\n",
      "0.23  128\n",
      "0.22  140\n",
      "0.22  129\n",
      "0.21  117\n",
      "0.20  124\n",
      "0.20  131\n",
      "0.19  147\n",
      "0.19  121\n",
      "0.18  125\n",
      "0.18  123\n",
      "0.17  140\n",
      "0.17  110\n",
      "0.16  122\n",
      "0.16  131\n",
      "0.15  119\n",
      "0.15  137\n",
      "0.14  129\n",
      "0.14  117\n",
      "0.13  126\n",
      "0.13  120\n",
      "0.13  140\n",
      "0.12  132\n",
      "0.12  121\n",
      "0.11  116\n",
      "0.11  116\n",
      "0.11  124\n",
      "0.10  118\n",
      "0.10  125\n",
      "0.10  133\n",
      "0.10  126\n",
      "0.09  133\n",
      "0.09  109\n",
      "0.09  119\n",
      "0.08  120\n",
      "0.08  120\n",
      "0.08  111\n",
      "0.08  115\n",
      "0.07  118\n",
      "0.07  125\n",
      "0.07  112\n",
      "0.07  129\n",
      "0.07  130\n",
      "0.06  115\n",
      "0.06  119\n",
      "0.06  128\n",
      "0.06  122\n",
      "0.06  127\n",
      "0.06  124\n",
      "0.05  109\n",
      "0.05  132\n",
      "0.05  113\n",
      "Runtime:  0.998150\n"
     ]
    }
   ],
   "source": [
    "starttime = time()\n",
    "\n",
    "MCS = 300 # MC steps per temperature\n",
    "J = 100 # meV\n",
    "T = 200 # Kelvin\n",
    "T_p = 0.97\n",
    "kB = 8.617e-2 # meV/K\n",
    "\n",
    "MCA_plot = []\n",
    "T_plot = []\n",
    "E_plot = []\n",
    "\n",
    "# set up the sites\n",
    "lat = lattice(50,50,50) # set up the size of lattice\n",
    "lat.populatesites()     # Assign each site a coordinate\n",
    "lat.initialisespins()   # Assign each site a random spin\n",
    "E_0 = lat.gettotalenergy1(J) # get initial total energy of the system\n",
    "E_t = E_0\n",
    "\n",
    "while (T>0.05):\n",
    "    \n",
    "    MCA = 0    # Initialise MC accepts\n",
    "    \n",
    "    for i in range(MCS):\n",
    "        fs = lat.pickrandomsite() # pick a random site to flip its spin\n",
    "        lat.filpspin(fs) # flip spin\n",
    "        dE = lat.sites[fs].getdeltaexchangeenergy(J)*2\n",
    "\n",
    "        if dE <= 0 or np.exp(-dE/(kB*T)) > rd.random():\n",
    "            E_t += dE\n",
    "            MCA += 1 # MC Accepts\n",
    "            \n",
    "        else:\n",
    "            lat.filpspin(fs) # flip the spin back\n",
    "    \n",
    "    print('%.2f  %d' %(T, MCA))\n",
    "    MCA_plot.append(MCA)\n",
    "    T_plot.append(T)\n",
    "    E_plot.append(E_t)\n",
    "    T = T_p * T\n",
    "    \n",
    "print (\"Runtime:  %.6f\"  %(time()-starttime))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "539bfe72",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAERCAYAAABl3+CQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcEUlEQVR4nO3da5Bc9X3m8eeRhAwjiTWMBCuwNcJZLVnWFW6zxPfaBOGAllh4U/aKHWHZuEoWsh2xiTfGpaqs88Jb9ibxRtnNaGqMtSjWLJAQu1A5im0kx/FeYscjjDFYYGEsCQUZRsPaIIkCW/rtiz4DrVb36cv0Oacv309V13Sfc7r7x5mhH/0v59+OCAEAUMucogsAAHQ2ggIAkIqgAACkIigAAKkICgBAKoICAJCqZ4PC9jbbz9p+pIFj/6vth5LbD23/NIcSAaAruFevo7D9DknHJP15RLyxied9VNKVEXFrZsUBQBfp2RZFRHxT0nPl22z/ku2v2N5r+3/Z/uUqT71Z0t25FAkAXWBe0QXkbFzShojYb/tXJY1K+vWZnbaHJF0i6esF1QcAHadvgsL2QklvkfSXtmc2v6bisDWS7ouIk3nWBgCdrG+CQqVutp9GxBUpx6yR9OF8ygGA7tCzYxSVIuJ5ST+2/R5JcsnlM/ttXyrpPEl/X1CJANCRejYobN+t0of+pbYP2/6gpBFJH7T9PUmPSlpd9pSbJd0TvToNDABa1LPTYwEA7dGzLQoAQHv05GD24sWLY/ny5UWXAQBdY+/evUcjYkm1fT0ZFMuXL9fk5GTRZQBA17B9sNY+up4AAKkICgBAqkKDwvb1th+3/YTtO6rst+0/TfY/bPuqIuoEgH5WWFDYnivpzyTdIOkySTfbvqzisBskrUhu6yVtzbVIAEChLYprJD0REU9GxMuS7tHpF8ApefznUfItSa+1vTTvQgGgnxUZFBdLeqrs8eFkW7PHSJJsr7c9aXtyamqq6WImJia0fPlyzZkzR8uXL9fExETTrwEAvajIoHCVbZWXiTdyTGljxHhEDEfE8JIlVacC1zQxMaH169fr4MGDiggdPHhQ69evPyMsCBMA/ajIoDgs6fVlj18n6ekWjpm1zZs368SJE6dtO3HihDZv3vzK40bDBAB6TZFB8R1JK2xfYnu+Skt876w4Zqek9yWzn94k6WcRcaTdhRw6dKju9kbCRKLVAaD3FBYUEfELSR+R9FVJ+yT9RUQ8anuD7Q3JYbskPSnpCUmfk7Qxi1qWLVtWd3sjYdJIq4MgAdB1IqLnbldffXU0Y8eOHTEwMBAqjX+EpBgYGIgdO3a8cszQ0NBp+2duQ0NDDR/TyPvMHDc0NBS2Y2ho6Iz9ANBukiajxmdq4R/qWdyaDYqI+h/OjXzI264aFLYjorGwaeR9CBIA7UZQtEm9D+h6QVAvSBp5DYIEQBYIipzU+xBvpEUx21ZJvRoIEQDVEBQ5SvsgbsdYyGyCpNExEgD9h6DoILMdC5lNkDQ6RpIWdLRGgN5EUHSZ2bRK0sKgXmsk7bVpjQC9jaDoMa0GSb0WRdr+RsZGaG0A3Yug6DO1PrTrtQrSWhxp++q1RAgQoPMRFHhF2gd3qy2KWvsGBweZgQV0CYICDWl1jKJWa6PWrd4MLAIEyB9BgYa1MuupVoui1i1tBlZaK4QAAbKTFhQu7e8tw8PDMTk5WXQZfWNmMcTy1XUHBgZ0zjnnaHp6+ozjh4aGdOjQITXztzc4OKgXX3zxjPcYHx+XVFrd99ChQ1q2bJk+9alPaWRkZBb/RUD/sb03Ioar7qyVIN18o0WRv2r/2m9lBlazN1ogQHuIricUpdkZWIODgwQIUIC0oKDrCYWZmJg4o8tIUlPdWM2q1YW1bt067dq1i+4r9K20rieCAh2niACxfdqYycz4B2GBfpEWFEV+FSpQ1cjIiA4cOKBTp07pwIEDGhkZ0cjIiMbHxzU0NCTbGhoa0vj4uLZs2aKBgYHTnj8wMKDBwcGm3rPyH0wzX3PLNxICtCjQA7JsgQwMDNBNhb5A1xP6UqMBUtntNGPu3Lk6efLkGdurdVMRHuh2BAVQpjJAVq1ape3bt5/Rcih/XE9leJx11lk699xz9dxzzxEc6AqMUQBlKsdARkdHq45/DA0NNfyalf/g+vnPf67p6WlFhA4ePKj169dr48aNjHegK9GiAGqodsV5rW6qRtBlhU5GiwJoQbWZVhs2bDhjlpXthl6v2syqsbExHTx48JWWxwc+8AEtXryYVgc6CkEBpGikm6paeDSqXpcVwYFOQNcT0AblA+Tnn3++XnjhBb388suv7J9Nl1U5BsmRFbqegIyVtzyOHj2qbdu2ta3LqhwtDhSBoAAykHWX1QyCA3mg6wkoUL0uq9miqwqNousJ6FBpXVaDg4OaP3/+rF6fFgfagaAAOkjewXHLLbfINqGBVAQF0MGyDo6ZrmdaG0hDUABdJMvgoLWBWggKoItlGRzlrQ1Co78VEhS2z7f9gO39yc/zqhzzett/a3uf7UdtbyqiVqCbZBUcdFH1t6JaFHdI2hMRKyTtSR5X+oWk342IfyHpTZI+bPuyHGsEul6t4JBau+BPoouqHxUVFKslbU/ub5d0U+UBEXEkIh5M7r8gaZ+ki/MqEOhFM8EREfrCF77Q9tYGodGbigqKCyPiiFQKBEkXpB1se7mkKyV9O+WY9bYnbU9OTU21s1agJ2XR2igPjbVr12rx4sUERg/ILChs77b9SJXb6iZfZ6Gkv5J0e0Q8X+u4iBiPiOGIGF6yZMlsywf6TrXWhtR6aEjS9PQ0rYwekFlQRMTKiHhjldv9kp6xvVSSkp/PVnsN22epFBITEfHFrGoFcLp2dlFV65pavHgxA+JdpKiup52S1iX310m6v/IAl/4Z83lJ+yLisznWBqBMO7uoZkJjenr6tAFxuqk6W1FB8WlJ19neL+m65LFsX2R7V3LMWyXdIunXbT+U3FYVUy6AGVl0UUl0U3UyVo8F0BYTExPatGmTpqen2/J6AwMDGh8fZ7XbnLB6LIDMjYyM6OjRo9qxY0dbWhknTpzQ2rVraWF0AIICQFtl0TXFOEaxCAoAmUmbPbVgwYKmX49xjGIQFAByUTl76tixYy11U3FRX/4ICgCFqWxxzJ07t+nXmJ6eJjAyRlAA6AgjIyPavn27BgYGWno+gZEdggJAxxgZGdH4+PisBsAJjPYjKAB0lPLuqFOnTrU83ZbAaB+CAkBHqxzHGBwcbOr5MzOlNm7cmFGFvY+gANA1yi/qayYwIkJbt26lddEiggJA12k1MGa6o+bOncu1GE0gKAB0rVYD49SpU5JeXfqcbql0BAWArtdqYEh0SzWCoADQM2YC47bbbmt6ai2D3rURFAB6zujoaEszpGZaFzPfwkcLo4SgANCTZtMdJXEdRjmCAkBPq/Y9Gc2gS4qgANAnyi/ca/U6jH4NC4ICQN9pddC7X8OCoADQt1oZ9O7HqbQEBYC+1sqgd7+NWxAUAKBXAyMidNttt9U9fmbcYtGiRT3fuiAoAKDC6OhoQ2EhSceOHdPatWt7unUxr+gCAKATjY6OSpLGxsZe+Z7uNFu3bj3teb2EFgUA1NDsYHevzooiKAAgRflU2kb0YlgQFADQgGbGLXotLAgKAGjQ6OioduzYoQULFtQ9tpdmRBEUANCEkZERHTt2rKHWRa/MiCIoAKAF/dQVRVAAQIv6JSwICgCYhWbDohvHLAgKAJilZsLiQx/6UMbVtF8hQWH7fNsP2N6f/Dwv5di5tr9r+8t51ggAzWh0RtTx48e7rguqqBbFHZL2RMQKSXuSx7VskrQvl6oAYBYanRG1detWrVy5MqeqZq+ooFgtaXtyf7ukm6odZPt1kv6NpDvzKQsAZm90dFQLFy5MPWbPnj1dExZFBcWFEXFEkpKfF9Q47k8k/Z6kU/Ve0PZ625O2J6emptpWKAC0YmxsrO4xe/bs6YpuqJpBYXt4Ni9se7ftR6rcVjf4/BslPRsRexs5PiLGI2I4IoaXLFkym9IBYNZGRkYaGuDuhplQacuMf872Qkl3S7onIn7QzAtHRM02le1nbC+NiCO2l0p6tsphb5X0LturJJ0t6VzbOyJibTN1AEBRZpYcn1mCvJZbbrlFUilcOlHNFkVEXCnpRkknJd1n+yHbH7c91Ib33SlpXXJ/naT7q7z/JyLidRGxXNIaSV8nJAB0m9HRUV177bWpx0SE1q1b17Eti9Qxioh4PCL+ICIuU+kD/bWSvm77/8zyfT8t6Trb+yVdlzyW7Yts75rlawNAR9m9e3fdsDh58mTHXmPR0GC27TkqDThfKGmBpFmNFkfEdERcGxErkp/PJdufjohVVY7/RkTcOJv3BIAi7d69u+5MqE69xiI1KGy/3faopMOS/qOk/y3p0oi4KYfaAKCnNDITqhMHt9NmPT2lUpfQPklXRsQ7I2JbRPwst+oAoIc0OhOq07qg0loUb4uIt0bEf4uIZ2zX/6YOAECqRtaFOn78eEe1KtJmPR2UJNtvtv0DJcto2L486Y4CALSgkSu3O6lV0chg9p9I+g1J05IUEd+T9I4MawKAnldvvKKTBrYbmvUUEU9VbDqZQS0A0DcaGa/olIHtRoLiKdtvkRS259v+mFjNFQBmrVu6oBoJig2SPizpYpWmyV6RPAYAzFIjXVBFtyocEYUWkIXh4eGYnJwsugwAaMiiRYt07Nixmvvnz5+vl156KdMabO+NiKqLwdZcFND2n6a9aET89mwLAwCUWhVr19Zeyu7ll1/WxMREYYsGpnU9bZD0NklPS5qUtLfiBgBog0YGtoscq0gLiqWSxlWaGnuLpLMk7YyI7RGxPeV5AIAmjY6OynbN/UWOVaRdcDcdEWMR8WuS3q/SyrGP2r4lp9oAoK9s2LAhdf+tt96aUyWnqzvryfZVkm6XtFbS34huJwDIxOjoqM4+++ya+2fGKvKWtijgH9jeK+l3JP2dpOGI+GCz33QHAGjcnXfembp/06ZNOVXyqprTY22fkvSkpBeTTTMHWlJExK9kX15rmB4LoJvNmTNHaZcuZHFZQ9r02LSup0skXavS16HeKOk3k9vMfQBABuqNVeTd/cQFdwDQgdJmQC1YsCD1Ar0W36+lFgUAoCBz5tT+eM57qixBAQAdqN4FdnkOaqfNehq3/W7bi3KrBgAgqTRVNs309HROlaS3KLZJulzSLtt7bH/c9uU51QUAfW9wcDB1f17dT2lXZn8rIj4ZEW+X9F5JhyT9ru3v2t5m+725VAgAfWrLli2p+/Pqfmpp1pPtqyVdHxGfan9Js8esJwC9Im32k9S+aypaWmY8TUSwgiwA5GBwcDDX8YhqmPUEAB2sXvdTHuMUXHAHAB0urftpcHBQR48ebcd7NH/Bne3fK7v/nop9/3nWVQEAZi2Pbqm0rqc1Zfc/UbHv+gxqAQBUUW+abNbSgsI17ld7DADISL1xio0bN2b6/mlBETXuV3sMAMjIyMhI6v6xsbFM3z9teuzltp9XqfVwTnJfyePaX8EEAGi7hQsX1lwxNutJSWlXZs+NiHMjYlFEzEvuzzw+K9OqAACnybrVkKaQ6yhsn2/7Adv7k5/n1Tjutbbvs/2Y7X2235x3rQDQCep1P2U5TlHUBXd3SNoTESsk7UkeV7NF0lci4pdVWqBwX071AUBXybLFUVRQrJa0Pbm/XdJNlQfYPlfSOyR9XpIi4uWI+GlO9QFAx1m4cGHNfVmOUxQVFBdGxBFJSn5eUOWYN0iakvQ/khVr77S9oNYL2l5ve9L25NTUVDZVA0CBihqnyCwobO+2/UiV2+oGX2KepKskbY2IKyUdV+0uKkXEeEQMR8TwkiVL2vBfAACdpd44RVYyC4qIWBkRb6xyu1/SM7aXSlLy89kqL3FY0uGI+Hby+D6VggMAUEVWCwQW1fW0U9K65P46SfdXHhARP5H0lO1Lk03XSvpBPuUBQPep9z3brSoqKD4t6Trb+yVdlzyW7Yts7yo77qOSJmw/LOkKSSxGCKCvzZlT+2P7+PHj2bxnJq9aR0RMR8S1EbEi+flcsv3piFhVdtxDybjDr0TETRHx/4qoFwA6RVathjR8cREAdJHR0dHc35OgAACkIigAoIdkMfOJoACAHrJp06a2vyZBAQA9JIuvRiUoAKDLpK35lAWCAgC6TN5rPhEUANBl8l7ziaAAAKQiKAAAqQgKAEAqggIAkIqgAIAe0+6rswkKAOgx7b46m6AAgB7T7quzCQoA6EKDg4O5vRdBAQBdaMuWLbm9F0EBAF0oz6uzCQoAQCqCAgCQiqAAAKQiKAAAqQgKAEAqggIAkIqgAACkIigAAKkICgBAKoICAJCKoAAApCIoAKAHbdy4sW2vRVAAQA/aunVr216LoACALrVw4cJc3oegAIAuNTY2lsv7FBIUts+3/YDt/cnP82oc9x9sP2r7Edt32z4771oBoFPl9Z0URbUo7pC0JyJWSNqTPD6N7Ysl/bak4Yh4o6S5ktbkWiUAoLCgWC1pe3J/u6Sbahw3T9I5tudJGpD0dPalAQDKFRUUF0bEEUlKfl5QeUBE/KOkP5J0SNIRST+LiK/VekHb621P2p6cmprKqGwA6D+ZBYXt3cnYQuVtdYPPP0+llsclki6StMD22lrHR8R4RAxHxPCSJUva8x8BANC8rF44IlbW2mf7GdtLI+KI7aWSnq1y2EpJP46IqeQ5X5T0Fkk7MikYAFBVUV1POyWtS+6vk3R/lWMOSXqT7QHblnStpH051QcASBQVFJ+WdJ3t/ZKuSx7L9kW2d0lSRHxb0n2SHpT0/aTW8WLKBYD+5Ygouoa2Gx4ejsnJyaLLAIDMlTpcqmvm89323ogYrraPK7MBAKkICgBAKoICAJCKoAAApCIoAACpCAoAQCqCAgCQiqAAAKQiKAAAqQgKAEAqggIAkIqgAACkIigAAKkICgBAKoICAJCKoAAApCIoAACpCAoAQCqCAgCQiqAAAKQiKAAAqQgKAEAqggIAkIqgAACkIigAAKkICgBAKoICAJCKoAAApCIoAKCLLVy4sKntrSAoAKCLjY2Nad68eadtmzdvnsbGxtr2HgQFAHSxkZER3XXXXRoaGpJtDQ0N6a677tLIyEjb3sMR0bYX6xTDw8MxOTlZdBkA0DVs742I4Wr7aFEAAFIVEhS232P7UdunbFdNsOS4620/bvsJ23fkWSMAoKSoFsUjkv6tpG/WOsD2XEl/JukGSZdJutn2ZfmUBwCYMa/+Ie0XEfskyXbaYddIeiIinkyOvUfSakk/yLxAAMArOnmM4mJJT5U9Ppxsq8r2etuTtienpqYyLw4A+kVmLQrbuyX90yq7NkfE/Y28RJVtNadoRcS4pPHkvadsH2yo0DMtlnS0xedmibqaQ13Noa7m9GJdQ7V2ZBYUEbFyli9xWNLryx6/TtLTDb73klbf1PZkrSliRaKu5lBXc6irOf1WVyd3PX1H0grbl9ieL2mNpJ0F1wQAfaeo6bHvtn1Y0psl/bXtrybbL7K9S5Ii4heSPiLpq5L2SfqLiHi0iHoBoJ8VNevpS5K+VGX705JWlT3eJWlXjqVJyThHB6Ku5lBXc6irOX1VV08u4QEAaJ9OHqMAAHQAggIAkKqvgsL2623/re19yVpTm5Lt59t+wPb+5Od5Zc/5RLLW1OO2fyOjurbZftb2I2Xb7rX9UHI7YPuhZPty2y+W7WvfovON1fVJ2/9Y9v6ryvZlfq7q1PaHth+z/bDtL9l+bbK96HNW6N9XlRovLTsXD9l+3vbtab/bvCR/699P3n8y2Vbz/OVUU63PjcLPV0Wd2a2NFxF9c5O0VNJVyf1Fkn6o0jpS/0XSHcn2OyR9Jrl/maTvSXqNpEsk/UjS3AzqeoekqyQ9UmP/H0v6/eT+8lrH5VGXpE9K+liVY3M5V3Vqe6ekecn9z5T9Hos+Z4X+fdWpd66kn6h0sVXV323O9RyQtLhiW9Xzl2NNtT43Cj9fFb/HH0l6g6T5yd/VZe16/b5qUUTEkYh4MLn/gkrTbi9WaQ2p7clh2yXdlNxfLemeiHgpIn4s6QmV1qBqd13flPRctX0uLYj1Xkl3t/t960mrq4pczlVabRHxtShNq5akb6l0kWauapyzQv++6rhW0o8iotWVDPJQ6/zlIuVzo5O8sjZeRLwsaWZtvLboq6AoZ3u5pCslfVvShRFxRCr9UUi6IDmsqfWmMvJ2Sc9ExP6ybZfY/q7tv7P99pzrkaSPJN0728q6ATrhXJW7VdLflD0u8px18t/XGp3+j5Bqv9s8haSv2d5re32yrdb5y13F54ZU/PmakenfUl8Ghe2Fkv5K0u0R8XzaoVW25T2f+Gad/j/yEUnLIuJKSb8j6X/aPjfHerZK+iVJVyS1/HGyvRPOlSTJ9mZJv5A0kWwq+pzVUug5c2nFg3dJ+stkU63fbZ7eGhFXqfT1Ah+2/Y4CaqiqyudGJ5yvGZn+LfVdUNg+S6Vf9kREfDHZ/Iztpcn+pZKeTba3vN5Um2qdp9L3dtw7sy3ppphO7u9VqV/yn+dVU0Q8ExEnI+KUpM/p1a6SQs/VDNvrJN0oaSSSztuiz5k69O9LpQ/jByPiGSn1d5ubKF10q4h4VqWLcq9R7fOXm2qfG51wvspk+rfUV0GR9Pd/XtK+iPhs2a6dktYl99dJur9s+xrbr7F9iaQVkv4hr3olrZT0WEQcntlge4lLX+ok229Ianoyr4Jm/odNvFulL6GSij9Xsn29pI9LeldEnCjbXug5U+f+fZ3WWk353ebC9gLbi2buqzQ54RHVPn951VX1c6Po81Uh27Xxih6tz/Mm6W0qNccelvRQclslaVDSHkn7k5/nlz1ns0r/An1c0g0Z1XW3Sk3Xn6v0L4MPJtvvkrSh4tjfkvSoSrMaHpT0mxmerzPqkvQFSd9PzuFOSUvzPFd1antCpX7amd/tWIecs0L/vmrUOSBpWtI/KdtW83ebU01vSH5H30t+X5uT7TXPX0511frcKPR8ValzlUozsn40c+7adWMJDwBAqr7qegIANI+gAACkIigAAKkICgBAKoICAJCKoAAq2B4sWxH0JxUrhM6vOPYTtkcqtr3f9n9P7s+xvT1Z4qHa1bOy/a9tf7lOTRtsv2+2/21AKwr5KlSgk0XpKu4rpNJS0pKORcQf1Tj8nSot2niGJBjGJJ0l6QMxi7noEZHZ0uhAPbQogBYl60XNj4ipGodsUelisfdFxKnkyuNttr+TLFB42uqeSevjgJPv0Ui2PWH7wuS7Dz6WbPuG7c/Y/gfbPyxoYUj0EYICaN1Kla4UrubfS7pa0pp4denzzZK+HhH/StKvSfrDZKkKSVKU1gy6X6XlIGT7VyUdiGQtpgrzIuIaSbdL+k9t+G8BaiIogNZdr9OXMi/3oEpfBlS+UNw7Jd3h0rcVfkPS2ZKWVTzvXkn/Lrm/RmULQlaYWdByr0pfzARkhqAAWneNai/i95hKYxf32v6XyTZL+q2IuCK5LYuIfRXP+3tJ/8z2EpW+oOeLqu6l5OdJMdaIjBEUQAuSD//HIuJkrWMi4v9K2iDpr20vk/RVSR+dmf1k+8oqzwmVltf+rEqrlU5nUT/QDP4lArTmBklfqXdQRHw5aR18RdJ1kn5f0sNJWBxQ6bszKt2r0rLR729XscBssHos0ALbD6g0m+lI0bUAWSMoAACpGKMAAKQiKAAAqQgKAEAqggIAkIqgAACkIigAAKn+P5zmWrFOEkd1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAiB0lEQVR4nO3df5AkdZnn8fdTTffc9SAn1IweB0zVwKJxg3oK7exy6xqoKOh54Im7O25zjqsRHY6eoXfnKnMd4WrE9d2qq3e4Hhq94Ry43fJjY/UkXBHRWzRWRexBQBDQQbuHURQYI1aO3hCcfu6PyiqqayqzsrLyV1V9XhEZU51VmflkVk49+c3vjzR3R0REBKBSdAAiIlIeSgoiItKipCAiIi1KCiIi0qKkICIiLScUHcAgtm3b5vV6vegwRESGysGDBx9z9+3d3hvqpFCv11lZWSk6DBGRoWJma2Hv6faRiIi0KCmIiEiLkoKIiLQoKYiISIuSgoiItCgpjLDl5WXq9TqVSoV6vc7y8nLRIYlIyQ11k1QJt7y8zNzcHOvr6wCsra0xNzcHwOzsbJGhiUiJqaQwoubn51sJoWl9fZ35+fmCIhKRYaCkMKIOHz7c13wREVBSGFk7duzoa76ICCgpjKyFhQWmp6c3zZuenmZhYaGgiERkGCgpjKjZ2VkWFxep1WqYGbVajcXFRVUyi0gkG+ZnNM/MzLgGxBMR6Y+ZHXT3mW7vjWVJQe33RUS6G7t+Cmq/LyISbuxKCmq/LyISbuySgtrvi4iEG7ukoPb72VFdjcjwG7ukoPb72WjW1aytreHurboaJQaR4TJ2SSHN9vu6Mn6a6mpERoP6KSTU2YoJGiWOce0gVqlU6HYumRkbGxsFRCQiYdRPIQO6Mt5MdTUio0FJISG1YtpMdTUio0FJIaG8royHpd4i7bGWhmW/RUaOuw/tdN5553lRlpaWfHp62oHWND097UtLS0O1jTIa1/0WyQuw4iG/q5n9YANnAH8H3AfcC7wrmH8KcAvwo+Dfk9uW2Q8cAh4ALuq1jSKTgnvjx6tWq7mZea1WS/1Hq1arbfphbE61Wi3V7ZTNuO63SF6ikkJmrY/M7FTgVHe/w8yeARwEXge8Gfilu/+ZmV0RJIX3mdku4FpgN/AvgK8Cz3H3Y2HbGPVRUse1Rc+47rdIXgppfeTuD7v7HcHrx2mUGE4DLgWuCT52DY1EQTD/Onf/tbv/hEaJYXdW8Q2DcW3RM677LVIGuVQ0m1kdeBHwHeDZ7v4wNBIH8KzgY6cBD7UtdiSY17muOTNbMbOVRx99NNO4izauLXrGdb9FyiDzpGBmJwJ/A7zb3X8V9dEu8467h+Dui+4+4+4z27dvTyvMUhrXp6eN636nTS24JIlMezSb2STwReBmd/9YMO8B4AJ3fziod7jV3Z9rZvsB3P2/B5+7GfiAu387bP2jXqcgkpR63EuUQuoUzMyATwP3NRNC4EZgb/B6L/CFtvl7zGyLme0EzgZuzyo+kVGmHveSVJZPXvtd4N8D3zezO4N5/wX4M+AGM3srcBj4fQB3v9fMbgB+APwGeEdUyyMRCace95JUlq2P/t7dzd1f4O4vDKYvuftRd3+Fu58d/PvLtmUW3P0sd3+uu9+UVWyi+82jLmkLLp0XUniv5EGmojuvDSv1GB59Sb5jnRfjgyJ6NOcxKSkkM4o9hrPuXV4G/e5jv58fxfNCulNSkE3MrOt/fjMrOrRExuEKN499HLXzQsJFJQWNkhpilO+tFtFjOMvjWcaWNmnvbx77qJ7kAqik0M2oX3nmvX9Zb69sV7hZ7G8e+zjq5708Dd0+6s843FvN8x58msezW9z9rj/tfe9cX7VaTf38yeucHIe6GVFS6FvZrjyHXVrHM+xKdt++fbGvcNO+Gu62vrBpkPNHV/GSJiWFPo1DSSFPaR3PqPXEvcJN+7sNW18W54+u4iUtSgp9yuuqrMz/ydOMLa3jmUaJI+1SYNj6Oidd1UuZKCkkkPUPdplvB2QRWxrHM42r/LxKCtVqtbQJX0RJoYTKfIuqrLGlkazyqFNIK7ln3VlNxpeSQgmVuTK7zLGl8cOXdeujtBJCP8mmzCVPKR8lhRIq69V42WPLStmusvv9DobpOyvbsR5HSgolVOYruzLHloUy7m+/pbUyl+7alfFYjyMlhYL0uiIq8xVTmWNLWxmvsoexpBDnnClDnKKkUAhdEQ2PMl5lD1udQtztl/FYjyMlhQLoimh4lPW7yrL1UdolwbjHsKzHetwoKRRAV0TDo+ir7LwVOWDfuB3rslJSKICuiOIpS91F3DjKEu8g0jw3owYkDFvnKBzDYaekUABdEfU2bMdo2OINk+UAhcN+bMaFkkJBdEUUreylqTyGxC5CWse9VwlB53t5RSUFa7w/nGZmZnxlZaXoMCShSqVCt/PPzNjY2CggoqctLy8zNzd33NPOuilDvP3otm/T09MsLi4yOzsbez1l/v4kmpkddPeZbu/pcZxDbNgfGVrmxz92e/xlmDLE24/Z2VkWFxep1WqYGbVare+EAOX+/mQAYUWIYZjKfvsoS6Nwf7vM+6AhsXsr8/cn0VCdwugp+/34uMpa7xJ2fDUk9mZl/f4kWlRSUJ3CkNL93Gyldd9dpIxUpzCCdD83W2ndd5d89FO/lkVd3LDX720SVoQYhmmcbx/pfq5IQz//F7J6quCw/V9EdQqjSfdze9MxGn391K9lURc3jPV7UUlBdQoyslQvMB76qV/Loi5uGOv3CqlTMLMDZvaImd3TNu+FZnabmd1pZitmtrvtvf1mdsjMHjCzi7KKS8ZHt74G6+vrzM/PFxRRtDLfl04jtm7rSGO9/dSvZVEXN3L1e2FFiEEn4KXAucA9bfO+Arw6eP0a4Nbg9S7gLmALsBN4EJjotY1xv30k0YZppNoy35dOI7Zu65icnPSpqamB91l1Cv2jqDoFoN6RFG4G/jB4/Ubgs8Hr/cD+js+d32v9Sgr5G6Z79MN0r7fMsaYRW9Q4Sf2ut9s5WOSzJNKIKW9lSgr/EjgMPAT8FKgF8z8BXN72uU8DbwhZ5xywAqzs2LEju6Mmxxm2K6JhirfMpZo0YovbQ7zXeoflOy17nGVKCh8HLgte/wHw1eD1/+qSFC7rtX6VFPJV5qvZMGW+WmtX5LHtdYzKVFIYlnOw7HGWKSn8A7RaPBnwq+C1bh8NgTJfzQ67oq4s42y3THUKw3IOlj3OMiWF+4ALgtevAA4Gr89hc0Xzj1FFc+mU/epn2BVRqon7naYRWxr33fM4B5Pua/tyExMTXeOsVqupxTmIQpICcC3wMPAUcAR4K/AS4GCQAL4DnNf2+XkarY4eIGih1GtSUshX2e+TSv/KfkXbKetzMOn6ez2FrjlNTU2V4v9LYSWFrCclhfwNyz16iaespb+o8yyNK/mw5ZIejzRbVw2yj3EpKYhIV2Us/RXZlyBpySmt1lVZ7X8nJQWRgpW5hJVXu/24sii9hK2z8x5/GUoKeZTelBREClTGq/EsDbq/WdRzRF3Jp9HSKs0e23nU8ygpiBSorPftszLo/kZd1add+ugWV5JWUktLS16tVjfFmrRXs0oKSgoy4oathc+gBt3fLMZJWlpayuwef9olQdUpKCnIiMviyjcPSesFwvZ3YmIicUuh9qvwpFfPYetoricstl5X7mnE1mv/1fpISUFGSJYjhGZlkKvVqDb7SfcxjdJWr74EYbFFbXuQEkiRBkoKwFagErx+DnAJMNlruTwmJQUZFllc+Q6y/ax7Di8tLYX26k2yj2ndZ28eh6gSQz/b7tXqqIwlQPfBk8JBYBo4jcbopp8Hlnstl8ekpCDDKs96hiRX/VmOjJpkH9O+z95PbFHbjtM/oUwlwKZBk8Idwb/vBN4bvP5er+XymJQUZFjl2SIpybbSiK+fdcQpyfTbKijqvV6lhc51V6tVr1arx62rV0kh6xJgUoMmhe8B5wO3AecE877fa7k8JiUFGVZ59l1IcsWe1cio3daRZt+A5nJJWgy1T3HrfOKOeVS2+oVBk8JLgRuB9wV/nwl8vNdyeUxKCjLM8urlnPSqP6uRUePG16u1UpJ7/e373Kt+IUkpJ816lCwNmhR+P868IiYlBZHesmhHn2Yyi7ovPz097fv27eu6vagSUD9X7EWNW1Tk0CeDJoU74swrYlJSEIknrR+gLG579bpa7/zRbm4vbLnmvf+4V+z9lBbiXPHHrR8pcuiTREkBeDXwF8AvaDxGszldDdwetlyek5KCjLMirjSzqCCPe1++c3thP6xhzX2b/QribD+sTqFZagFat4r6qeDO8jj2I2lS+FfAXmAt+Lc5vR44OWy5PCclBRlXRV1pZtWUNqpfQ9T2uv0AR90Oitp+r5ZN+/bt69kpL+73UvTQJ4mSQusDMAlMAS8Ang9M9Vomr0lJQcZVUVeaWW632w9qP7eBso6x122mWh+dEqPWlUepb9Ck8BoandZuBb4OHCbm4zKznpQUZFwVdaWZdQklztV5r+1lFWM/FdK9vpekQ26kZdCkcD/wW21/nwXc32u5PCYlBRlXRd6TzrsuI8n2ksaYtMMbEHn7K6opaz/LpGXQpPCNjr+tc15Rk5KCjKuiW6+MokE6vPWqKM96SJF+DZoUPgl8CXgzjYrmLwIfpVHh/Ppey2c5KSnIOCuynXsZpL3/cUpf7Vf3na2Poq76m5+tVqu+devW1vzm8BllKilY4/1wZva/I952d39L5AoyNDMz4ysrK0VtXkQKsry8zNzcHOvr661509PTLC4uMjs7m2idlUqFbr+HZsbGxkaimOKYmJhgYmKCJ598sjVv0H3pxcwOuvtM1zfDssUwTCopiIynLOpUojrDxS2RJBk6Aza3stq6dWvXwffSxIAlhefQuIX0bHd/npm9ALjE3f9r0iyVFpUURMbToFf13XS70p+cnMTM+r6KD4sviSxKDVElhUqM5f8S2A88BeDudwN7UotORKRPO3bsiD1/eXmZer1OpVKhXq+zvLzcddnZ2VkWFxep1WqYGbVajZNOOmlTQgBYX19n7969kesLiy+J9fV15ufnU1tfT2FFiOYEfDf493tt8+7stVwek24fiYynrIflbkr6EJ0kQ3dETWm3RGLA1kc30eib0HzYzhuAm3otl8ekpCAyvtIYY6jXOuLWDyTph9DPVK1WU21tNWhSOBP4KrAO/BT4e6Dea7k8JiUFEYkS1QcgTikijYfoLC0tDZwUJiYmYj30J66opNCzornJzLYCFXd/PNYCOVBFs4hEqdfrrK2tHTe/VqsBhL63urra+nt5eZn5+XkOHz5MpVLh2LFjPZfptG3bNo4ePdr/DvTQa7thBqpoNrP/ZmbPdPcn3P1xMzvZzApveSQi0svCwgLT09Ob5k1PT7OwsMDhw4e7LtM5f3Z2ltXVVTY2NrjmmmtC1xflyiuvZHJyMsEeRAvbh4GEFSGaE20VzG3zej5kBzgAPALc0zH/ncADwL3Ah9vm7wcOBe9d1Gv9rttHIhJD2L34pCOVtq+v2SM5bh+GsN7LW7duTTTgXrVaTXRMGLBO4W5gS9vf/xS4N8ZyLwXObU8KwMto1E9sCf5+VvDvLuAuYAuwE3gQmOi1DSUFEUlq0JFK0x5/at++fX0nhampqUTbGzQpvJdG5fJbgbcEr9/ba7lg2XpHUrgBuLDL5/YD+9v+vhk4v9f6lRREpF2/LXQGGak07phFvR7gEzX+UZwpSS/ugZJCY3kuBv6cxkB4sW7tePekcCfwQeA7NJ7N8OJg/ieAy9s+92ngDSHrnANWgJUdO3b0fTBEZDQNcuXe70ilUS2K2peJ+6jPQaYkfRiikkKciuadwK3u/h53/8/AN8ys3mu5ECcAJwO/A/wJcIOZGY3huDt5txW4+6K7z7j7zPbt2xOGISKjZn5+/rjB6OL2Bu6nh3RzW3HW1S2mp5566rhe0oNIs/c0xBvm4q+B9sFEjgXzkjgCfC5IVrcH690WzD+j7XOnAz9LuA0RGUNxWxN1E9VKqd91ti+TSeugNnFaPvUrTlI4wd1baS14PZVwe/8HeDm0BtqbAh4DbgT2mNmWoGRyNnB7wm2IyBjq92q/XbdxjxYXFwFa4yZt27aNbdu2UalUqFS6/3RWq9VNA9elfRUPtLbdjDH14bXD7iv50/fwb6ExKmrz70uBr8VY7lrgYRoD6R2hUVE9BSwB9wB3AC9v+/w8jVZHDxDzGdCqaBaRprRbA/U7flFaYyBNTExEPtozaYujdgzY+ugs4DbgMPAQ8C3grF7L5TEpKYhIuzTHB4ozbtHExETPbS0tLUX+yHdO1WrVTzzxxFitjooe5uJEwLzRq/nF7v7dWAtmSMNciEhW4jwTIeunsvWS9FkLgz5PoWkH8Cdm9kMaD90RERlZceoDwj7T+QwH4Lg6i2q12nXZiYmJ2DFm8ayFyKRgZjUzu8LM7gL+Cng78KqwDCMiMiq6tUhqF9byp1kqWFtbw91ZW1tjbm4OoDWG0urqKldeeWXXFk/dBtyLknYLp9CkYGbfAr4ETNLoSHYe8Li7r6YagYhICXW2SKpWq1Sr1U2tk7rdtonbX6Jbi6e9e/f2VVKA9Fs4hdYpmNkXgBfRaC76WXf/lpn92N3PTDWCAahOQUTKJunzo5PUO0xNTXHgwIF86hTc/VLg+TSajn7QzH4CnGxmu/vauojIGEnaX6JbCaOXycnJ1PspRNYpuPs/uPsBd38l8NvA+4H/aWYPpRqFiMiI6Ld3dFOSuoEnnniC5eXlvpeLErv1kbs/4u5/4e7/GnhJqlGIiIyIsN7RzSv6ZsskM6NSqWBmrSmJtFsfxe6nUEaqUxCRYZJFf4W4fSU6lkmln4KIiAwgSb1BL0WMkioiIilIu09BrqOkmtmHzextXeb/RzP7UKpRiIiMgTSv6rMaJTWqpPBaYLHL/CuBf5NqFCIiY6BXL+m4zIzV1dX0h80mOim4ux9XexHMS1ZNLiIyxtpbJg3ilFNOSSmi40UlhXUzO7tzZjDvHzOLSERkhM3OzrK6uoq7s7S0lGgdjz/+eOr9E5qiksL7gZvM7M1m9vxg+mPgb4P3RERkALOzs6GjpUZ58sknede73pVBRNHDXNwEvA54GXB1ML0MuMzdv5RJNCIiY6bbaKmTk5M9O7MdPXo0k9LCCVFvuvs9wN7UtyoiIgCtyuL5+XkOHz7cqi84evRoz2X37t27aR1piBol9caoBd39ktSiSEg9mkVklCTp8Zzk6WtRPZqjksKjNJ7JfC3wHTpaHLn712NHkBElBREZJfV6nbW1tb6Xq9VqrK6uxv58VFKIun30z4FXAm8E/ohGBfO17n5v/FBFRCSupD2e0+wpHVXRfMzdv+zue4HfAQ4Bt5rZO1PbuoiItCTt8ZxmT+lez2jeYmavB5aAdwAfBz6X2tZFRKQlSY/ntMc/Cr19ZGbXAM8DbgI+GLREEhGRjDQriy+//PLYy5x//vm5tT7aAJ4I/mz/kNEYAuOk1KJISBXNIjIqlpeXW81S+3lGQm7PU3D3irs/I5hOapueUYaEICIyKppNUdfW1nD3vn7k3T3VTmx6noKISMEGffhOmo/kVFIQESnYoE1Kc2mSKiIi+Ri0SWluTVJFRCR7gz58J80mqZklBTM7YGaPmNlxTVnN7D1m5ma2rW3efjM7ZGYPmNlFWcUlIlI27Q/fMTOq1SqVSvyf5zSbpGZZUrgauLhzppmdQWP4jMNt83YBe4BzgmWuMrOJDGMTESmV5sN3NjY2eOyxx/jMZz4Tu/Rw4YUXphZHZknB3b8B/LLLW/8DeC+b+z5cClzn7r9295/QGFJjd1axiYiUzfLyMvV6nUqlQr1e55vf/GbPZyo0fe1rX0stjsjnKaTNzC4Bfurud3Xs7GnAbW1/HwnmdVvHHDAH6VauiIgUpXPI7LW1NT75yU8WEktuFc1mNg3M0/1Rnt3SYdeu1u6+6O4z7j6zffv2NEMUESnEoP0U0pRnSeEsYCfQLCWcDtxhZrtplAzOaPvs6cDPcoxNRKQwg/Yz2LVrV0qR5FhScPfvu/uz3L3u7nUaieBcd/85cCOwJxiVdSdwNnB7XrGJiBRp0FvhTzzxRO8PxZRlk9RrgW8DzzWzI2b21rDPBg/uuQH4AfBl4B3ufiyr2EREymTQfgpJntYWJrPbR+7+xh7v1zv+XgDS64EhIjIkmv0MmqOkVioVjh0r5rpYPZpFREqgvZ9Cv0Nhp0lJQUSkZPqtY5iYSK+vr5KCiEjJ9FvHMDc3l9q2lRREREqmfSykOK666qrUtp1rj2YREYmnWfn8pje9Kdc6BpUURERKan5+PvdKZyUFEZGSitvT+e1vf3tq21RSEBEpqbitkBYXF1PbppKCiEhJLSwsxHrYTpod3ZQURERKKs0nqsWlpCAiUmKqaBYRkZZ+ntWcyvZy3ZqIiMS2vLyc+zaVFERESkr9FEREpGXQJ7IloaQgIlJSgz6RLQklBRGRkhr0iWxJKCmIiJRUv6OlpkGjpIqIlFizA9vll1+ey/ZUUhARKbn5+fnctqWkICJScmtra7ltS0lBRKTk0nwGcy9KCiIiJZfmKKi9KCmIiJRcnq2PlBREREpuYWEht20pKYiIlFyez1VQUhARKbk8R0tVUhARKTn1UxARkZY8R0tVUhARKbk8R0vNLCmY2QEze8TM7mmb9xEzu9/M7jazz5vZM9ve229mh8zsATO7KKu4RESGzcLCAmaWy7ayLClcDVzcMe8W4Hnu/gLgh8B+ADPbBewBzgmWucrM8uvCJyJSYrOzs7z85S/PZVuZJQV3/wbwy455X3H33wR/3gacHry+FLjO3X/t7j8BDgG7s4pNRGTYHDp0KJftFFmn8BbgpuD1acBDbe8dCeYdx8zmzGzFzFYeffTRjEMUESmHvCqbC0kKZjYP/AZoNr7tdrPMuy3r7ovuPuPuM9u3b88qRBGRUsmrsjn3pGBme4HXArPu3vzhPwKc0fax04Gf5R2biEhZ5TXURa5JwcwuBt4HXOLu621v3QjsMbMtZrYTOBu4Pc/YRETKLK+hLjJ7HKeZXQtcAGwzsyPAn9JobbQFuCVoXnWbu7/N3e81sxuAH9C4rfQOd89vrFgREQHAnr6DM3xmZmZ8ZWWl6DBERHIR1Vehn99yMzvo7jPd3lOPZhERaVFSEBGRFiUFERFpUVIQEZEWJQURkSGQ14N2lBRERIZA1IN2JibSGz9USUFEZAhEjX10wQUXpLYdJQURkSEQNfZRmiOoKimIiAyBqLGP1tbWUtuOkoKIyBCYnZ2lUun+k606BRGRMbSxsdF1/rFj6Q0Vp6QgIjIkwkoEKimIiIyhsBKBSgoiImNIJQUREWlRSUFERHKlpCAiIi1KCiIiQ0J1CiIi0jI3N9fX/CROSG1NIiKSqauuugqAxcVFjh07xsTEBHNzc635abB+HvZcNjMzM76yslJ0GCIiQ8XMDrr7TLf3dPtIRERalBRERKRFSUFERFqUFEREpEVJQUREWoa69ZGZPQoM8sihbcBjKYWTJsXVH8XVH8XVn1GMq+bu27u9MdRJYVBmthLWLKtIiqs/iqs/iqs/4xaXbh+JiEiLkoKIiLSMe1JYLDqAEIqrP4qrP4qrP2MV11jXKYiIyGbjXlIQEZE2SgoiItIysknBzM4ws78zs/vM7F4ze1cw/xQzu8XMfhT8e3LbMvvN7JCZPWBmF2UU1wEze8TM7mmbd72Z3RlMq2Z2ZzC/bmb/2Pbep7KIKSKuD5jZT9u2/5q29zI/VhFxfcTM7jezu83s82b2zGB+0cer0HOrS4zPbTsWd5rZr8zs3VHfa16C8/z7wfZXgnmhxy+nmMJ+Mwo/Xh1xXhycR4fM7IrUN+DuIzkBpwLnBq+fAfwQ2AV8GLgimH8F8KHg9S7gLmALsBN4EJjIIK6XAucC94S8/1Hg/cHretjn8ogL+ADwni6fzeVYRcT1KuCE4PWH2r7Doo9XoedWj3gngJ8DtbDvNed4VoFtHfO6Hr8cYwr7zSj8eHV8jw8CZwJTwXm1K81tjGxJwd0fdvc7gtePA/cBpwGXAtcEH7sGeF3w+lLgOnf/tbv/BDgE7M4grm8Av+z2npkZ8AfAtWlvt5eouLrI5ViFxeXuX3H33wR/3gacnsW2+42Lgs+tHl4BPOjug4wAkLWw45eLiN+MMtkNHHL3H7v7k8B1NI5bakY2KbQzszrwIuA7wLPd/WFonATAs4KPnQY81LbYEfI/IX4P+IW7/6ht3k4z+56Zfd3Mfi/neAD+Q3Cb5kBbcb4Mx6rpLcBNbX8XebzKfG7tYfPFRrfvNU8OfMXMDppZ81mSYccvdx2/GVD88WrK/Fwa+aRgZicCfwO8291/FfXRLvPybq/7Rjb/x30Y2OHuLwL+E/BZMzspx3g+CZwFvDCI5aPB/DIcK8xsHvgNsBzMKvp4hSn0eJnZFHAJ8NfBrLDvNU+/6+7nAq8G3mFmLy0ghq66/GaU4Xg1ZX4ujXRSMLNJGl/usrt/Lpj9CzM7NXj/VOCRYP4R4Iy2xU8HfpZjrCcArweub84LbjccDV4fpHEv8Tl5xeTuv3D3Y+6+AfwlT9/yKPRYAZjZXuC1wKwHN1uLPl6U9Nyi8cN7h7v/AiK/19y4+8+Cfx8BPh/EEHb8ctPtN6MMx6tN5ufSyCaF4P78p4H73P1jbW/dCOwNXu8FvtA2f4+ZbTGzncDZwO15xQtcCNzv7keaM8xsu5lNBK/PDGL6cV4BNf+DBv4d0GxpU+ixMrOLgfcBl7j7etv8Qo8X5T23NpVAI77XXJjZVjN7RvM1jYYD9xB+/PKKq+tvRtHHq8N3gbPNbGdQAtxD47ilp+ja9Kwm4CU0ilV3A3cG02uAKvA14EfBv6e0LTNP4+ryAeDVGcV1LY0i6FM0sv5bg/lXA2/r+OxlwL00WhjcAfzbDI/XcXEBfwV8PziGNwKn5nmsIuI6ROO+avN7/VRJjleh51ZInNPAUeCftc0L/V5ziunM4Du6K/i+5oP5occvp7jCfjMKPV5d4nwNjZZRDzaPXZqThrkQEZGWkb19JCIi/VNSEBGRFiUFERFpUVIQEZEWJQUREWlRUpCxZmbVttEvf94xGuZUx2f3m9lsx7w3m9kngtcVM7smGAqhW89TzOwCM/tij5jeZmZvGnTfRJI4oegARIrkjR7QL4TGEMnA/3P3Pw/5+KtoDFh4nCAJfAqYBP7YB2jr7e6ZDfkt0otKCiIxBGMoTbn7oyEfuZJG56s3uftG0Gv3gJl9Nxigb9NIlkGpYtWCZ0EE8w6Z2bOD8fvfE8y71cw+ZGa3m9kPCxoUUcaIkoJIPBfS6GXbzR8B5wF7/OkhveeB/+vuLwZeBnwkGNIBAG+Mo/MFGsMmYGa/Dax6MD5RhxPcfTfwbuBPU9gXkVBKCiLxXMzmIbrb3UHj4TXtA6W9CrjCGk/RuxX4J8COjuWuB/4weL2HtsEQOzQHczxI40FCIplRUhCJZzfhg9jdT6Ou4XozOyeYZ8Bl7v7CYNrh7vd1LPdt4LfMbDuNB8p8ju5+Hfx7DNUDSsaUFER6CH7o73f3Y2GfcfdvAW8D/tbMdgA3A+9stkIysxd1WcZpDBv9MRojcx7NIn6RfuiqQ6S3VwNf7vUhd/9icNX/ZeCVwPuBu4PEsErj+Q+drqcxHPKb0wpWZBAaJVWkBzO7hUarooeLjkUka0oKIiLSojoFERFpUVIQEZEWJQUREWlRUhARkRYlBRERaVFSEBGRlv8Pv4oV4DJMbIkAAAAASUVORK5CYII=\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,E_plot,'ok')\n",
    "plt.gca().invert_xaxis()\n",
    "plt.xlabel('T / Kelvin')\n",
    "plt.ylabel('E / MeV')\n",
    "plt.show()\n",
    "plt.close('all')\n",
    "plt.figure()\n",
    "plt.plot(T_plot,MCA_plot,'ok')\n",
    "plt.gca().invert_xaxis()\n",
    "plt.xlabel('T / Kelvin')\n",
    "plt.ylabel('MC Accepts')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df4e59a0",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}