exam/2020/Exam.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Python Exam (Master, 13/10/2020): Ising Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> + **Allocated time:** 2h\n",
    "> + **Send your final notebook** at the end of the exam at [romain.madar@cern.ch](mailto:romain.madar@cern.ch)\n",
    "> + **Allowed documents:** documentation of [python](https://docs.python.org/3/), [numpy](https://numpy.org/devdocs/user/index.html), [matplotlib](https://matplotlib.org/contents.html), [pandas](https://pandas.pydata.org/pandas-docs/stable/) and [scipy](https://docs.scipy.org/doc/scipy/reference/) as well as the documentation from the notebook interface using `Shift+tab` keyboard shortcut, `help(module)` or `help(function)`.\n",
    "> \n",
    "> The final mark will be a number between 0 (very bad) and 20 (very good). The evaluation of this exam takes into account mainly the correctness of the answers, but also the clarity of the explanations and the quality of the code. *There is no required knowledge about Ising model to answer the questions asked in this exam.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## General informations\n",
    "\n",
    "### A bit of context ...\n",
    "\n",
    "The Ising model is famous piece of statistical physics, which is originally meant to describe magnetic material made of interacting magnetic dipole moments at given temperature. This problem is *a priori* extremely complex and the Ising model offers a simplified approach to model macroscopic magnetic phenomena, in particular phase transition. In practice, this model describes a system of coupled discrete system, such as spin 1/2, forming a n-dimensional lattice. Interactions are assumed to take place only between closest neighbours, accouting for the short-range nature of magnetic interactions.\n",
    "\n",
    "This exam proposes to implement a simulation of such a system and explore some its dynamical properties. In a first part, a 1 dimension lattice is considered (*i.e.* spins are disposed on a line) and the Metropolis algorithm will be implemented for the thermal evolution. In a second part, a 2 dimension lattice will be investigated using a provided python package."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import packages and plot cosmetics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "mpl.rcParams['legend.frameon'] = False\n",
    "mpl.rcParams['legend.fontsize'] = 'xx-large'\n",
    "mpl.rcParams['xtick.labelsize'] = 16\n",
    "mpl.rcParams['ytick.labelsize'] = 16\n",
    "mpl.rcParams['axes.titlesize'] = 18\n",
    "mpl.rcParams['axes.labelsize'] = 18\n",
    "mpl.rcParams['lines.linewidth'] = 2.5\n",
    "mpl.rcParams['figure.figsize'] = (10, 7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Ising model in one dimension (12pts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**1.1. [1pts]** We want to numerically create a system of N spins in one dimension (spatial positions are determined by one coordinate). Spin can be 'up', *i.e.* with a value $+1$ or 'down', *i.e.* with a value $-1$. Using the numpy function `np.random.randint()`, creates a numpy array corresponding to `N=100` spins in a random state. Create a function `initState(N)` which return a spin system of N spins and print the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "scrolled": false
   },
   "source": [
    "**1.2 [1pts]** Write a function `magnet(s)` which compute the total normalized magnetization of the system `s`, defined as \n",
    "$$\n",
    "M = \\frac{1}{N} \\sum_{i} s_i\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**1.3 [2pts]** Using an explicit `for` loop, create 10000 times a new spin system, compute $M$ for each of them, and plot and histogram of 10000 $M$ values, using a 100 bins between -1 and 1. Repeat this for N=10, 100, 1000. We will keep N=100 after this question, unless specifed otherwise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**1.4 a) [2pts]** We now consider a case of 100 spins. The goal of this question is to compute the energy of the system. Close-by aligned spins leads to a lower energy, while anti-aligned close-by spins has a energy cost. At the end, we can define the normalized energy of the system as:\n",
    "$$\n",
    "E = -\\frac{1}{2} \\frac{1}{N}\\sum^{N-1}_{i=1} ( s_{i}s_{i-1} + s_{i}s_{i+1})\n",
    "$$\n",
    "\n",
    "Write a function `energy(s)` returning the total energy of the system `s`. Since the first (last) spin has not left (right) neighbour, we will simply exclude them from sum in the energy calculation (but we keep them in the normalization factor $N$).\n",
    "\n",
    "**HINT:** the correct approach doens't use explicit for loop. One can start by defining a 'bulk' spin system, exluding the first and the last element. In a next step, one can try to get all left neighbours in one array and all the right ones in another array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**1.4 b) [1pts]** Which value of magnetization and energy do you expect for a system with all spins up? Check the result explicitely by generating a state with all spins at 1 and by calucating energy and magnetization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**1.4 c) [1pts]** Repeat the question for a system of spins where they are all anti-aligned (*i.e*, one up, one down, one up, etc ...).\n",
    "\n",
    "**HINT:** It can be useful to remember how to select every two elements in an array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**1.5** The goal of this question is to implement a thermal evolution of the spin system, using a Metropolis algorithm, defined by the following steps:\n",
    "  1. flip randomely $n$ spins\n",
    "  2. compute the energy difference $\\Delta E$ due to the flipping step\n",
    "  3. accept the new configuration with a probability 1.0 if $\\Delta E<0$, and with a probability $e^{\\frac{-\\Delta E} {T}}$ if $\\Delta E>0$, where $T$ is a new parameter of the simulation, obviously called temperature. If the new configuration is not accecped, the lattice state is unchanged but the procedure still count as one thermal evolution step."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**1.5 a) [3pts]** Write a function `thermalEvolution(spins, T, n)` which operate these three steps on a system `spins` and return the evolved lattice state.\n",
    "\n",
    "   + **HINT1:** flipping a spin can be seen by a multiplication by $-1$.\n",
    "   + **HINT2:** to generate $n$ random integers, between $0$ and $N$, one can use `np.random.randint(low=0, high=N, size=n)` command.\n",
    "   + **HINT3:** an acceptation rate following $e^{\\frac{-\\Delta E} {T}}$ can be obtained by accepting if $r<e^{\\frac{-\\Delta E} {T}}$, where $r$ is a random number between $0$ and $1$ (which can be obtained with the function `np.random.rand()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**1.5 b) [1ps]** Plot the total magnetization as function of time for 5000 thermal evolutions (we will use an explicit `for` loop), where for each thermal evolution, 30 spins are flipped. Compare the evolution of the magnetization with time for three temperature values, 0.02, 0.03 and 0.1 (the spin system has to be initialized for each temperature before starting the \"thermalization\")."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Ising model in two dimensions (8pts)\n",
    "\n",
    "Since the 2D simulation is a bit more tricky, a python package is provided for the exam. The main object of this package is a spin lattice, which can be created and vizualized as follow:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAGOCAYAAACABOR8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAd+0lEQVR4nO3dPa4cV3MG4KHhBTD4HNuAtAMSZEBBS3DoRSr0EgQpoEDugAwcO2HkwDAwjiRfN+7tnneq6vTww/NkxHSfPjPTc4t9fqpeXa/XCwDc6h/O7gAA3xeBA4CIwAFAROAAICJwABAROACI/GNy8N/+9rfrP//Lv9x1oa+fv/y/f//w5se72um49p5tvzrPTV4/um5nP1N71zpqe/K7SN5X5+c1+T0nfUt/Y8nx1d9vcs886m+uIu3n1uRn8tKx/3X5r8t/X//71XPHvUr2cbx5+/b6+8c/bj7+qX/7x3/9f//+5X/+/a52Oq69Z9uvznOT14+u29nP1N61jtqe/C6S99X5eU1+z0nf0t9Ycnz195vcM4/6m6tI+7k1+Zm8dOyvl98u367fng0chqoAiERDVVuV/4V0/u/xyF4ErvxvZ3t+5/9YOv+3eKTS71Xf0y2vV56E9to+6sfk03PlPaf93vtdVL+bvWPT1/e+myPJd1V54u1+eqk8We65t5+eOACICBwARAQOACLRqqrXr15ff7789Ne/q/MDT03Oj1RWdDzKaocjne+x83uszGWtPPdR5uvOXH3YKVm9NDlG33lvH7W958wVW/eyqgqANgIHAJHSBsBVQxrdj+uV5birlvZuTQ4DVobnOr+blY/zRyrv8XsZbqoMJ3UOIz7KUGll8+BW58bObXuTmy+f+vD+3eXzp0+GqgCoEzgAiAgcAERKKUcqy9Eq45yTY6SpzmWGnckDJxMmdo6FJ/3aWjkv1qlzju3e6x7pTv/SOedRSRtS8aj3YyW1zF5bXy8vZ+T1xAFAROAAICJwABBZlnLke0lR3NnWZIGalValzJicIzpSWTd/VoGpbdvVfkzOvTzKfqfO+3GyENbUfGfStn0cALQROACICBwAREq5qrZWjTF3joV3jheuTF++9Shjt0fOGu9fuacmaXtlHqyz9jdtr91ZymBlWYSts/YKdZYM2GtLWnUA2ggcAEQEDgAiY/s4OsduJ3PbHx2beKS12smxK8eYz6q58ajlSs/8fB51X1HXGP0917r1ukeqc3+T382t78M+DgDaCBwAREpp1bdWDRVUh0emHjnT6+4d35kqefLzepTyuNW2JlN3r7zXu9ruvmdWLXtNr7v3nivvcTJdzuQS5Fv77YkDgIjAAUBE4AAgMpZyZOU8Q2fKkcpSy85UFEc6xzGnUsXc8/qeyX5W/L0sc7332FuO70zPsSpF+951u9s6Ky39Xj8sxwWgjcABQETgACBS2scxmcZhr63OLfbVccupcc2j99j5nivpTVam8khMpWF4rq2VKUgqey8S1e911bzPynmHlWnUk/ec/n47fmeeOACICBwARAQOACKt+zgSZ67JT8boJ/PTbE3tL0mPfZRcX2lbq9bkP9J3sXfspMo9M3m/paZKG6wsBTHVltKxALQROACICBwARKJ9HF8/f7l5THByTXR1DDSZlzhr/LWznsRk/qOjfmwlNUiSc4/6Uf2ek2Mn65kc9aVL9+93aq7rrLofz+mcL0nPX80TBwARgQOAiMABQCTax/H61evrz5efWi78vdYGWJW36ei6icp69Oq1j/py1nj15HtKrOxH8rup3hOT9bgrv7lHrcOy9Qi1Z+zjAKCNwAFAJFqO+8ObHy+/fLztMbEzrUD34+jK5al71105XNd1bNpWJaXz5D1TSZFxZHLIcjJd+WT68s7l4cmwV2fKoGTZdfWe2ZpaXn/v/eSJA4CIwAFAROAAIFJKqz6Vknj7evfStUfbvv+nR1nOV732nkcpQ9uZrnwyRftRW1Np1qtzGlOlYlOdpWOPjk/6sfL+3Lp1CbzluAC0ETgAiAgcAERaU45MpSxYOfaYtjWZMjuxakx0e/7KdNCT+2Am00fs6RzrPrIyNcrkfpPKuavux+7f/hlzgeY4AGgjcAAQETgAiES5qhKdpTqPrBxnn1z/35mva3KvSmeK7Keq8xBT+yW6c1VN5VpKTabuPisn2VHbyfnpZz25F2gqLf3W09c+vH/34nGeOACICBwARAQOACKlOY5VZUDTfhw5a114Z52QR/o8nzozB8/3qPo9duZimpx72Zsz6tyfk6rUrfle93F09M0TBwARgQOAyFjp2COV5WmT6SI6Sz52liutLFXtTrfRmQ6mM33J1JLFzqWn2/ZWfn5HbSU676lqyvGKs1K4TLZ11Hbld/InTxwARAQOACICBwCRaI7j6+cv0dzDnqMx0c7lZ1OlTo+OTdtKxh5XpivvbKtz+WNiMt12JTXFyn7t9SN57Z5rdX7vnfMlybmd5Rw652nTOaGOvx2eOACICBwARAQOACJR6dg3b99ef//4x1//XlWWsXvscWUJzVXX/R5L2KZtHamkzEiu3bkP5sjkGP3KdC+r5iVSq/aXdJbe7Zxr2bvuh/fvLp8/fVI6FoA6gQOAiMABQCSa43j96vX158tPd12oc7x6spxmRbUfZ829HJmcP1k13l/RXZ50VX6uzhT3RyZTjneWjl1ZljbRmXtu697v+dfLb5dv12/mOACoEzgAiAgcAERK9Ti2OkvHTuboqajk90ly8Kwsidl5ncn8Pp3j6JM5otJrJdftrIFT2cPQmZ9r79hbrl3pV+LMHG6dtY5uPffD+3cvHueJA4CIwAFAZNly3K3JMqrptROr0jScuUwzuVbn0NTKZdcr0790Dk2l5yfOWo47mXLkrCXc3Val4n/aluW4ALQROACICBwARKLluIkzxx4nr905zl4Zm6yU3t1KS3U+PX7lstaKzrQgnXNG6TLNzuWlnfNNlWtVrCrtcItVc4ErS1S/xBMHABGBA4CIwAFApHWOYyr9QXW7fmeK561kvPredm9pu3O+JDm+c0w+fQ9nvcdUkiKj0lbnXEz3nMRUKvmqybmsxMq2O+aBPHEAEBE4AIgIHABESnMclfXBjzqm15kv6cw5j+Tc9NoVlc8vafvMNfhHpsboO+eEUunegq62O8uodn6v1bm/yRT4t15XWnUA2ggcAEQEDgAiy+pxTJZdTM/vbKuyDnzlOvG9625N5hGbqh2QmqztUfmeK59ndax7qq2tlXVVVtYzmexH5bu592/Wh/fvLp8/fVKPA4A6gQOAyLK06p0pMyaHhCpprtPlkJVUKJVH21TyXRydm3x+K0u0Vt7jyiHIzvQ5lSXdR211Oqu0bOV+nEyrUl1q3vFdeeIAICJwABAROACIlJbjdo49TqVhOGp7ZRnQyjxO6qz0EpX5pkcq+/m93o9d/Tg6N7UqHdGjvMejtleWQbj3N2c5LgBtBA4AIgIHAJHSPo7KuvmV8yGd+xC2x3eut+7c03DvdW9pe2rvwGQa+kq/jo5N2967zqOk2+7eEzL1PVd1zlFO7j3bk/6t7fguPHEAEBE4AIgIHABEWvdxPNU5VnvkrHw0fy8eNWV2airFffcc21SupckyqtXfc2cepz0rc+JV+nF0/Fbn39pb2ccBQBuBA4CIwAFA5LR6HMn53eurnx4/WcOg0tbKsdrKOvDUqjmiyRokneee2c/O655VLjft11lWzp2uyNfliQOAiMABQETgACASzXH88ObHyy8f789T9FRnneZ0fHCqDkZnDq2qqToXadtn1ZdIr7uybsNUbepkbq/aryOdeeySvFede1fSfk3q7Kea4wAsJ3AAEImGqr5+/nLzY86jLou7XNaldO9MLb9yOXOns0pkTqbenxxuq5xfXXa9d2zne+7s51ZaBmFqOX362+9Meb9iCM0TBwARgQOAiMABQKS0HHdr1djvZIrnzpKOqcoyzcn5kqNr77XdmaK9U+d4f+U9dqffvvfY1ORnsLXXdnUpfrKceWty7mtP5z1yL08cAEQEDgAiAgcAkdI+jsl14Z1pLSbLgu5dZ6tzP8WZewkq6UyS+abJFC6TKUY6U3us2ndQPbeSJiS99p7vZe9UxWSZg6dtf718efE8TxwARAQOACICBwCRUunYynjh5H6JlfMWXdc5uu6jzGl0t520ldwjnXMv1XP35m46v+fqb6iSoj15vbPsbKqzpMKqsghH153cC/QSTxwARAQOACICBwCR0hzHqnw/3WuiO2tVJDn4k351nruyFsrKWgGT8xJ7x26dOWZ/lsl5su5r76nMJey11Vked3t8tQzyrXODH96/e/E4TxwARAQOACKtadUTZy1l274+mZb5qF8rh0eSfiXXWtnPI8ky1zOHzCrLiDvT0FfSvU+Wf+1MGZQMN1VTG+3p3I5QGU5P+/USTxwARAQOACICBwCR1rTqU+P9R85MV16Zp5gsm9qZlj5JhZ7qnJfYm3tZmcZi5VLzyncxtXy52tbKVDN7OtO/VHWmbJdyBIDlBA4AIgIHAJFSypE91dQTlTXmR6ZK3lbXslf2RyT9XLnvIHFmKo7KeH/lfqyuyZ/6zKqpujv3qnTud1pVcmFyvqmzZMC9798TBwARgQOAiMABQKR1jqMydlYpT3pmuc1KivbO/D570v0kyRh05/xSZz8rKbGP2kqvlex3Wrn/aVLnPEXl+LPSuad7o5JrVecob53//Hr58uJxnjgAiAgcAEQEDgAipTmOzvHWs3IJda6jX1kD4sjkmHLn+OvU57dyHX1lbqtzH8LkXovq5zm136SyHyxte7LWx1alpO2KvzOeOACICBwAREqlYyuPU1tnpZvovG7nssIzh22O+pKYSulQNTnMkJxf/Z1Uhr2S63TqHAbrTNFS/fz2+jXZ9hlLsj1xABAROACICBwARB6mdOyezpKiR20n5x6ZTPdeaffM8eu911emy6+oppbZ07mk/cheP6tLeyu/74pVadOP2l75m+u8H2/tpycOACICBwARgQOAyKvr9XrzwW/evr3+/vGPv/49OXfQuT49Wfu+clxzcn9EJS1IJS345HtOTKYcOStV93M65xIe5X1NpmTv3L8zubeqMn9cvfaffr38dvl2/fbqudc8cQAQETgAiAgcAERKadU7125Pln9NcgWlbSfHpvlrbr3uUV/S8qQr18JXdKYJn5wreJRcXyvLle6dP7lvY++6z/2789qVEgGVcrhT9+OH9+9efM0TBwARgQOAiMABQKQ1V9Wezrz5nePVVUl+/0epodG5j2MyZ1HnHpvOfR1n7v1JdO4VmNyv0/k7SefvOlXm2Cr5zFbm7vuTJw4AIgIHABGBA4BIqeZ44sz9Ep25bibzXiX1ELZWjpV3jscm5ybHT+bjSvpxpDpHVHnPlX5Vzl+5J2QyZ1Zn/feVv7Gkny/xxAFAROAAIDK2HLc7FXpictimc8ndZIqHrnOfO3/V8tLOJdyVYZyjflTa7lxamX72K9NcTN4HlbZWLRs+ait9PTl24u+hJw4AIgIHABGBA4DIactxO5cOptdOHI3VVuY89lQ/r8rSwM7ymiuXpib9mGzrrGWxk2lV0n51vufEZHmGyvGdf/86073c+7fREwcAEYEDgIjAAUDk1fV6vfng169eX3++/PTXv1ft45hMhXxkciy8UnoyufZkmvnqfpxV3+3Kz2AynUkllUdnKorJND8r09Q/6v131j3z1K+X3y7frt9ePfeaJw4AIgIHABGBA4DI2BzH1soSoxWTbVf6sbVyT0Pnvo7J73kyd1DlOp35kfZMzotN3jPpuZ2lDFaZnAucKjNhjgOANgIHABGBA4BIlKvqyNTeizPH5CtWrqPfmix9utfPyRKYnfU30r7cet2qznxTR8eu/F107iua2u+0d52jY4/6lfZj73eV5gW79z1+eP/uxfM8cQAQETgAiETLcd+8fXv9/eMfg935P8njaGV55OTSyslU6EdWLVl81KW9nUMt1eWOZy1fPmqrc1nwqqXQ22udOaQ7uUz4EVKhWI4LQBuBA4CIwAFAJFqO+/Xzl7vHa9MxvMl0Jp1LLZMx0kdJ0X7Ur8n30Tkv0bk8t7M87qOUBU3PnZxPSc7tfI9pids9nWmTqqV5K20/bWuvXctxAWgjcAAQETgAiLSmHNmzcnz1SDL+f9Zegu69K5W2HiWtdWfKlsl+HB2/9xlNrt9PxvtXplWp3H+PUgbhcqml4umcP9nr19a9f1c8cQAQETgAiAgcAERKpWMTnXMF1dTJj7Je/VHLqG5NjhtPzZdM5pPqHI+u5kdKTKXfvqVfle95z2RutDNTtO/p/D3uvWe5qgBoI3AAEBE4AIhE+zh+ePPj5ZePt43NpeuYp0o4PqdzvHVyr8bea5VSk2nbW5Nj4YlVNQuqJse+k7Yr+c0m905VPoPOub/qvTtV3vWo7VTHb9ATBwARgQOASCnlyGQqgMqwzco0BMkj+JFVQ2idywwndabLr6ZGr7TVdd3UZCqZyuudfwuqw7Cdw8NTS47Tfh25tZ/SqgPQRuAAICJwABApzXGsKlfanao7SX+c9KW6bHhVGdrU5HjtXrvpd3PW57fyntk7f3K+pHPJduffhup7rKTiXzlvUfmbtXXrXPTXy5cX2/DEAUBE4AAgInAAEInmOL5+/jJWKvWsNdHVtNbJ+Z3v46xU8dvXq6kVptqqmty7Ukkf0ZmKv3O/yco9DlN/g47aSlO43Hvs0fGde+T2fnP2cQDQRuAAICJwABBp3cfR2dZkWdWKyt6Vrc4ylp25lpLjO1PeV3OQVfIhdeqcWzhqu2LlfonOtODJnobJ/TyV39zKeduJ34InDgAiAgcAEYEDgEhr6dipdfST9SOqbU/tN6nWLFi1pvxI57xE9fV7VeeXOsfCE+kc0eQ+mErblT0NnTnxkn4dHZt6tP1OnjgAiAgcAERaS8fu6Vyq2llS9MjkcNNZy1zT15PH+862jkymXalcZ2XbncubO00uh55axt45nHTUryNnpTKSVh2AEQIHABGBA4DIq+v1evPBb96+vf7+8Y8XX1+V1rpzHDg991HSbVfSlR+1tTW5nHRqCePKVN1nlo7d07ksuHrtTsm8xMql0JPzElPpTPaO/fD+3eXzp0+vnjvPEwcAEYEDgIjAAUAkmuN4/er19efLT3/9e9XYbfe+jc4Uz4nJ8f6z5oSOTKbQWHmPJeduJfNPK9tKrnNkcs5t8u9M59zB5L625Dpd8zq/Xn67fLt+M8cBQJ3AAUBE4AAgMpZWfWUOo6l1zNMq48CVeZw09XRlHHjVuG9VZz6zSn6kFWU/77lWNe/V1PxT5x6byTTzqUffe+aJA4CIwAFAROAAIFLax1HROUafXitp+1HXW1es/PyStjv7lerct5Fca+Ucx+S+g7Nybj3K53PU1pEz51pfIlcVAG0EDgAiAgcAkWX7OFLJOPBeP6om65evmk9J19wn+zyq301lvqlzrfteP7trdVf261T2KB3p3O/UOS+x6nvv7Ff3nMVUTq17++WJA4CIwAFAJBqqSkyWXUzbepR0x53nJ8Mn1bT0lVQMlWNXLk3d0z20UhnqW/U9p7+D5J7pHLLtTKdT/V0kw3Gd3/uRidQpnjgAiAgcAEQEDgAirXMclXHjqXHz50yVWZ1MX36kM5X8WUs+j9ruHBtP0oRXU8dMjqt3qpQ2WNl20lbn0vLJOaIzl57fwxMHABGBA4CIwAFAJJrj+Pr5y90ptc8qLVk9t7Oc5KPsS9jqXEefvt65zyNpd2Va/840F3sm+9m9ryM5t9NkiYXK/pwjq9L+Pz336+XLi8d54gAgInAAEBE4AIiMlY5duea+Mp49uS4+6cflUhtXr6wpX7mvo3NebOrcybbSa5312VfPnUzRfu9102tPfp6VtibngJ5e99fLb5dv129KxwJQJ3AAEBE4AIiM5arqLDXZOd5fbfvoWklbnWOiicmcPJV9MGk/KmVBJ2tmbFXmnzrnCo+cNXeVfs9751b+7nTnl9o7Nm3r6eudc757bX14/+7FNjxxABAROACIlIaqKumNO7fQd6amriwNXLm0d+VjcnLuVmdbU+lJjnSm/T56feXw5coloEma+so9tTKtT2Vo9Ejnvb43nLd3rJQjALQROACICBwAREopRyZTFO9ZuWR2cjlf0lbnXELSj+m2O6/Tef91lndd9flVdaYBmjr3ufMTk0uOV7XdudVhz4f37y6fP32ScgSAOoEDgIjAAUAk2sfxw5sfL798vG2sbTKtxVY6H1BZI51YmeaiM0V2pZ+dJucSOsenJ62cS5ico6yk+d+q3H+VVCidfzcqe33Sa+0df+/8mycOACICBwARgQOASDTH8fXzl1NyNa1MJV3ZI9K5v+TI5L6OlaU699JFp21VdL7ns+aIquP9lXxSj1JG9aw5ySOdv5Mz93j9yRMHABGBA4CIwAFApLUex6r8PpVSk0dtbXWufa+MC1f6OTn3Mjm+Ork/p9KPIytrgUzWbUj6caTye+66znPXqpicT6lct2veRz0OANoIHABEorTqb96+vf7+8Y+//j2ZDqGSoqDyqDaZ2uOorc5lmpMpC/auMzkkmfTlzM9v5dLU73GIsnN4+MiZqdBvvc49197T1S9p1QFoI3AAEBE4AIi0znE81Tmm1zlmXLUq/fZkioIjlbTWZ47Jd6SLfk51DmNVapnOkqLdSz6n5k9Wzi8dmbr/Ul1/h369/Hb5dv1mjgOAOoEDgIjAAUCklFa9M43A5DxG5x6RybQMt17nuWtNprtfNV5b3ROyZzLVRDXVzN6xW5O/k0rq7q1Vc1lHbXfutVg5j7PXdppyaWJO2BMHABGBA4CIwAFAJJrj+OHNj5dfPt421taZw2hlSdajtibHq299bbXKuPBk3rDq8ava7SwxmhxfmS/pLqM6tU+rut+pcu5UGYRpt177w/t3L77miQOAiMABQETgACBSylXVaXJ9ddJeZ52BzrHbzvXpk+Vyj0zOcexdpzpGP2VyrmCllXMJlZxaR9dOzq20tbJWyr3XlqsKgDYCBwCRsaGqlUMBK1N7r0qNMpmW4VHTSaxKM3PP63vHVvu2Z7KEQJKKp/r5TS2x/Xsdnjs6fsrTfigdC0AbgQOAiMABQCRKOXLkUcYej1TGC1eVh03HiDtTtiR9qaYcqaS4r1g5d5D04yhFdqWtybTgK8s7nzVHOVkCeHKpfnKtp699vXx58ThPHABEBA4AIgIHAJHSHEeldOfK8ezJMdHKuPCqtdlH/ZhMX/69vMfOMeak7c77K51v2jor5X1l/u7oOivvv7NSClXsfZbSqgPQRuAAICJwABAp5apatXfgzPXpiWoa5lXr06t7bDr3XnSmqa+k2z5rLuZIZz8rn19nidbqHNv3mPfqzH1D937eclUB0EbgACAicAAQKe3jWDX+Ws2HtFVp+9Z2u9vuLIFZrb8xNZcwWU6zs61HKvdamZeo/OaO2tqaykmW5EKr6pxb6dzHkb7HjvvREwcAEYEDgIjAAUAkmuP4+vnL7ljlVD6azrHbW87fM1mbunPt++T8SSVXUOe4eqVfadtJW2flR5rs12Qupc59WKnKfN2eyb0qj7AHyRMHABGBA4BIa1r1VctLq4/NyXDb5CN6JU3IyqGDZHhp5ZBF5drJtY6WeHYvF79Xd4ngvdcq6WCS6x61VR1e6hwq3btuZ/mG6hB4R5oVTxwARAQOACICBwCRUlr1ramlld1p1c9KjTK57HCy7Odk+u1OnffM3rlHzkqXn5r8XWxVfieTZRE6l9iueo+pe+dLfr38dvl2/SatOgB1AgcAEYEDgEhpH0eiMpY7WSq202RakLPSQx+ZLCma6iwDuqqfyXU7jr+33UcqJTuVdr0zldHk53UkTe90D08cAEQEDgAiAgcAkVJa9YrJ3Cwrx/ArOY0m04LvXWdS5bvqzId0dG5l38GZecMqks+kOk5emdPo3GMzlb+s22QpiHv3m3x4/+7F1zxxABAROACICBwAREr7ODrrSWydlTdncox+qzOHUXJ+Z92BI5NzWVuVGi5blTX5e209116Xzv0Q3fM4k/V1Km09Sp2LzjxYRzrmajxxABAROACIjKUcWZlC/MwytInOVAmVto90pouYHOZa1Vb30F7nEu7OodLOpauJyb8Fye9mZaqiyne1Kv3L18uXF4/zxAFAROAAICJwABCJSse+fvX6+vPlpxdf71xaOVXScft6Z1qQo35tdaaSr5gs97qyJOZUPx5pvq7izJQaU3MxZ5W/7W67cq3O+bun5yodC0AbgQOAiMABQCSa43j16tV/Xi6X/5jrDgAP4p+v1+s/PfdCFDgAwFAVABGBA4CIwAFAROAAICJwABAROACICBwARAQOACICBwCR/wVNdMROb87HmQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import ising\n",
    "lat = ising.spinLattice(100)\n",
    "lat.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a quick description of the main functions availaible in the the `ising` package, which can be called on a `spinLattice` object:\n",
    "\n",
    " + `lat.spins` is directly the 2D numpy array of spin values\n",
    " + `lat.plot()` plots the current state of the lattice\n",
    " + `lat.energy()` computes and returns the total (normalized) energy of the lattice\n",
    " + `lat.magnetization()` computes and returns the total (normalized) magnetization of the lattice \n",
    " + `lat.randomize()` randomizes the state of the lattice (*i.e.* randomely flip all spins)\n",
    " + `lat.align()` aligns all spin of the lattice at +1\n",
    " + `lat.antialign()` anti-aligns all spins of the lattice\n",
    " + `lat.thermalEvolution(T)` performs one single thermal evolution under the temperature T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**2.1 [2pts]** Using an explicit `for` loop, vizualize the state of a $100\\times 100$ spins lattice after 0, 100 and 1000 thermal evolution steps, for a temperature of 0.5. One can use the instruction `spins.thermalEvolution(T)` to perform one single thermal evolution at a temperature T. Compare the result of the two following situations:\n",
    " + the lattice is randomized before thermal evolutions\n",
    " + the lattice is aligned before thermal evolutions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**2.2 a) [1pts]** Write a function `plotEvolution(T, mode)` which produce the three plots before for a give temperature `T` and a given initialization mode (*i.e.* randomized or aligned). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**2.2 b) [1pts]** Compare the evolution of the lattice, when it is aligned at first, for four temperatures, namely 1.0, 1.5, 2.0 and 2.5. What do you notice?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**2.3** The goal of this question is to investigate the behaviour of the magnetization with the temperature. \n",
    "\n",
    "**2.3 a) [2pts]** Write a function `magThermalEq(lattice, T)` which aligns the lattice `lattice`, performs 1000 thermal evolutions at a temperature `T` and return the final magnetization at thermal equilibrium."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**2.3 b) [1pts]** Plot the magnetization after thermalization, as a function of the temperature. One would consider 30 points between 1.5 and 3.0, and an explicit `for` loop over temperatures is allowed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**2.3 c) [1pts]** Does this behaviour depends on the lattice size? One could try a couple of smaller sizes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**BONUS QUESTION [3pts]** \n",
    "\n",
    "The goal of this question is to use some already-simulated data to investigate the correlation distance between pair of spins in the lattice. This can be done measuring the correlation between pairs of spins as function of the distance separing the two spins. It's not needed to fully understand this concept to answer the following question. At the end, all you need to know is that you have a function of r for different temperatures: $\\mathcal{C}_T(r)$.\n",
    "\n",
    "In order to make the use of these data easier, a dataframe `df` containing values of $r$ and correlations for several temperatures is loaded below. This dataset corresponds to a $30 \\times 30$ spin lattice. Also, the column name of correlations values for each temperature is given via a dictionnary `corrTs`. An example of how to access the information is given below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simulated temperatures \n",
    "Ts = [0.5, 1., 1.5, 1.75, 2, 2.1, 2.15, 2.175, 2.2, 2.225, 2.25, 2.269, 2.30, 2.35, 2.40, 2.5, 2.6, 2.8, 3, 3.5, 4]\n",
    "\n",
    "# Column name of the correlations data for each temperatures\n",
    "corrTs = {T: 'MeanT'+str(T) for T in Ts}\n",
    "\n",
    "# Full data frame\n",
    "df = pd.read_csv('TwoPointCorr_30x30.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmYAAAGgCAYAAAAAdau+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeXxcV33//9dHm2XLuyTHS7zEibOSxYkJgQCBEAhbEpYCLRD2tdAWeLTfkm8phQAtv2+/UPrtBqEk7CUUEhKWQmg2IE0CjrMQZ3XiLV5iWd5lWev5/TEzzkiWrStppBnJr+fjocfYd+6ZOff6euatc849J1JKSJIkqfyqyl0BSZIk5RjMJEmSKoTBTJIkqUIYzCRJkiqEwUySJKlC1JS7AqXS1NSUlixZUu5qSJIkDeqee+7ZnlJq7r99wgSzJUuWsHLlynJXQ5IkaVARsX6g7XZlSpIkVQiDmSRJUoUwmEmSJFUIg5kkSVKFMJhJkiRVCIOZJElShTCYSZIkVQiDmSRJUoUwmEmSJFUIg5kkSVKFMJhJkiRVCIOZJElShTCYSZIkVQiDmSRJUoWoKXcFxoM12/bxr7euoWVfB1e84hROnT+93FWSJEkTkC1mGezv7Oa6ezfx68e3s3Hn/nJXR5IkTVAGswyapk46+Oft+zrKWBNJkjSRGcwyaJxad/DP2/d2lrEmkiRpIssUzCJiYUT8ICJ2R8SeiLguIhZlLPu3EXFTRLRGRIqIdwywz4kR8Y8R8UBE7IuILRFxY0ScOcTjGRWTaqqZXp8bjtfaZouZJEkaHYMGs4iYAtwCnAy8HbgcWAbcGhENGd7jT4DJwE+OsM/LgBcD3wAuAf4YaAbujohzMrzHqCt0Z9qVKUmSRkuWuzLfCywFTkoprQGIiAeAx4H3A18cpPyMlFJvRJwAvO0w+3wP+JeUUipsiIhbgHXAnx2h3JhpmjqJJ7e32ZUpSZJGTZauzEuBuwqhDCCltBa4A7hssMIppd4M+2wvDmX5bbuBx4AFGeo46pqm5caZ2WImSZJGS5Zgdhrw4ADbVwOnlrY6z4iI2cCzgIdH6z2GotCV2WIwkyRJoyRLMJsN7Bxg+w5gVmmr08c/AQF86XA7RMT7ImJlRKxsaWkZxapAY0MumO090E1Hd8+ovpckSTo6ZZ0uIw2wLUpZkT4vHHEF8Gbgw8VdqIdUKqWrUkorUkormpubR6s6wDNdmQCt+xxnJkmSSi9LMNtJrtWsv1kM3JI2IhHxAeBvgU+klK4u9esPl5PMSpKk0ZYlmK0mN86sv1OBh0pZmYi4HPhX4Asppc+V8rVHymAmSZJGW5ZgdiNwXkQsLWyIiCXA+fnnSiIiXgtcA/x7SunPS/W6pdLcJ5jZlSlJkkovyzxmXwU+DNwQEZ8gN97sM8BG4CuFnSJiMfAEcGVK6cqi7ReQmyx2bn7TiojYB5BS+kF+nxcC/wE8AHw9Is4rev+OlNK9wzu80umzLJMtZpIkaRQMGsxSSm0RcSHwD8C3yA36vxn4SEppX9GuAVRzaCvcp4ELiv7+ofxPoQzAhcAkYDm5+dGKrQeWDFbP0dYwqYbJtdW0d/U4yawkSRoVWVrMSCltAF4/yD7rGOBOzZTSizK8/qeAT2WpSzk1Tatj4452W8wkSdKoyDpdhnC9TEmSNLoMZkNQCGbOYyZJkkaDwWwImqa6XqYkSRo9BrMhKLSY7djfSXfPoGuzS5IkDYnBbAgKwSylXDiTJEkqJYPZEBTP/u84M0mSVGoGsyFwkllJkjSaDGZD4HqZkiRpNBnMhqDPepnO/i9JkkrMYDYE0yfXUFedO2W2mEmSpFIzmA1BRBwcZ7bdwf+SJKnEDGZD1Ogks5IkaZQYzIbI9TIlSdJoMZgNkcFMkiSNFoPZEBUvZJ5SKnNtJEnSRGIwG6LCQubdvYnd7V1lro0kSZpIDGZD5CSzkiRptBjMhqg4mLU4yawkSSohg9kQNU1zvUxJkjQ6DGZDVNxi1mowkyRJJWQwG6JZU+qoityfnf1fkiSVksFsiKqrgtkNzv4vSZJKz2A2DE4yK0mSRoPBbBgKwazFrkxJklRCBrNhKEwy6+B/SZJUSgazYWgs6sp0WSZJklQqBrNhKHRlHujqpa2zp8y1kSRJE4XBbBgKXZkA2/fanSlJkkrDYDYMTdOKJpltM5hJkqTSMJgNQ1OD62VKkqTSM5gNg+tlSpKk0WAwG4bGohYzg5kkSSoVg9kw1NVUMWNyLWAwkyRJpWMwG6ZnJpl1jJkkSSoNg9kwNbpepiRJKjGD2TA1HwxmtphJkqTSMJgNU6Er0wlmJUlSqRjMhqmwLNPejm4OdLkskyRJGjmD2TAVxpgBtLbZnSlJkkbOYDZMrpcpSZJKzWA2TMXrZXpnpiRJKgWD2TA1TzWYSZKk0soUzCJiYUT8ICJ2R8SeiLguIhZlLPu3EXFTRLRGRIqIdxxmv6qIuCIi1kXEgYi4PyJeP4RjGVONxV2ZTpkhSZJKYNBgFhFTgFuAk4G3A5cDy4BbI6Ihw3v8CTAZ+Mkg+30G+BTwz8ArgLuA/4yIV2Z4jzE3pa6GKXXVgC1mkiSpNGoy7PNeYClwUkppDUBEPAA8Drwf+OIg5WeklHoj4gTgbQPtEBFzgD8HPp9S+r/5zbfmy3we+FmGeo65pqmT2LBjvy1mkiSpJLJ0ZV4K3FUIZQAppbXAHcBlgxVOKfVmeI+LgTrg2/22fxs4PSKOy/AaY85JZiVJUillCWanAQ8OsH01cGqJ6nEa0AGs6bd9df6xVO9TUoVJZlvbDGaSJGnksgSz2cDOAbbvAGaVqB6zgV0ppTTAexSeP0REvC8iVkbEypaWlhJVJbtG18uUJEkllHW6jP6BCSBKWI8YznuklK5KKa1IKa1obm4uYXWyac53Ze7c30l3T5YeW0mSpMPLEsx2MnCL1SwGbkkbjh3ArIjoH8RmFT1fcQqTzKYEO1yWSZIkjVCWYLaa3Biw/k4FHipRPVYDk4DjB3gPSvg+JdVUNMlsi1NmSJKkEcoSzG4EzouIpYUNEbEEOD//XCn8HOgE3tJv+1uBB/N3gVac4mDW6jgzSZI0QlnmMfsq8GHghoj4BLmxYJ8BNgJfKewUEYuBJ4ArU0pXFm2/AGgG5uY3rYiIfQAppR/kH7dFxD8AV0TEXmAV8CbgQjJMyVEufWf/t8VMkiSNzKDBLKXUFhEXAv8AfIvcgPybgY+klPYV7RpANYe2wn0auKDo7x/K/xTKFPwVsA/4M3Ih7lHgjSmlH2c+mjHW5HqZkiSphLK0mJFS2gAccd3KlNI6BriLMqX0oozv0QN8Nv8zLkyvr6GuuorOnl6nzJAkSSOWdboMDSAinpn93xYzSZI0QgazESpMmWGLmSRJGimD2Qg1NrhepiRJKg2D2Qg1HVyWyWAmSZJGxmA2QoWuzNa2Tnp7B1pVSpIkKRuD2QgVWsx6ehO727vKXBtJkjSeGcxGqMlJZiVJUokYzEbI9TIlSVKpGMxGqO/s/06ZIUmShs9gNkJ9ujKdMkOSJI2AwWyEZk2poyq/EFVrm8FMkiQNn8FshKqqgtkN+bnM9tqVKUmShs9gVgKulylJkkrBYFYCzdOc/V+SJI2cwawEnlmWya5MSZI0fAazEji4kPm+DlJyWSZJkjQ8BrMSKKyX2dHdy76O7jLXRpIkjVcGsxJwkllJklQKBrMScL1MSZJUCgazEihuMWs1mEmSpGEymJVA34XM7cqUJEnDYzArgUbXy5QkSSVgMCuB2uoqZk6pBRxjJkmShs9gViKF7sxWuzIlSdIwGcxKxPUyJUnSSBnMSqRxqutlSpKkkTGYlUiz62VKkqQRMpiVSKErc19HNwe6espcG0mSNB4ZzEqk77JMdmdKkqShM5iVSKPrZUqSpBEymJVIk5PMSpKkETKYlYhdmZIkaaQMZiXSPK1oIfM2uzIlSdLQGcxKpL62mqmTagBosStTkiQNg8GshBqd/V+SJI2AwayEmpz9X5IkjYDBrISeWS/TMWaSJGnoDGYlVGgxa7XFTJIkDYPBrIQKk8zu3N9FV09vmWsjSZLGG4NZCTUXTTK7wykzJEnSEBnMSqh4klmnzJAkSUNlMCuhJieZlSRJI5ApmEXEwoj4QUTsjog9EXFdRCzKWLY+Iv4+IrZERHtE3BkRLxxgv8aI+MeIeDK/39qI+OeIaB7qQZVLn2WZbDGTJElDVDPYDhExBbgF6ADeDiTgs8CtEXFGSqltkJf4GvAq4C+AJ4EPAb+IiOemlO7Lv0cANwInAp8EHgZOBT4DnBMRz0sppWEc35hqLF7I3DszJUnSEA0azID3AkuBk1JKawAi4gHgceD9wBcPVzAizgTeDLwrpXRNftvtwGrgSuDS/K7LgOcB708pXZXfdltE9AL/Ri6wPTq0Qxt70ybVUFdTRWd3r8FMkiQNWZauzEuBuwqhDCCltBa4A7gsQ9ku4Nqist3A94CLI6LQ91doatrTr/yuIdSz7CKC5oOz/zvGTJIkDU2WwHMa8OAA21eT624crOzalNL+AcrWAScU/f1XwF9HxIqImBoR55Lr1vyvlNLDGepZEZpcL1OSJA1TlmA2G9g5wPYdwKwRlC08T3782CvJdVf+DtgL3E1uTNrrD/fiEfG+iFgZEStbWloGqcrYaLTFTJIkDVPWLsKBBt5HhnIxhLJfBc4DPgBckH9cAfwgIgasZ0rpqpTSipTSiubmyrh50xYzSZI0XFkG/+8k37LVzywGbg0rtgMYaFqNWUXPExGvAv4IuCildHP+uV9FxJPATcAlwA0Z6lp2hSkzdrR10tubqKrKkl8lSZKytZitJjdWrL9TgYcylD0uP+VG/7KdQOGGgtPzj7/rt99v84+nZKhnRSgEs57exM79dmdKkqTssgSzG4HzImJpYUNELAHOzz83WNla4A1FZWuANwE3pZQK/X1b84/n9iv/nPzjpgz1rAjO/i9JkoYrSzD7KrAOuCEiLouIS8l1K24EvlLYKSIWR0R3RHyysC0/gey1wJci4j0R8RJyU2UcB/xN0XtcB2wGvhkRH4yIF0fEB4Fv5t/n+pEc5FhqaiiaZNbZ/yVJ0hAMGszyM/tfCDwGfAv4DrAWuDCltK9o1wCqB3jNdwLXkFst4KfAQuDlKaVVRe+xh9zA//8C/lfR44+B5/Z7n4pW3GLW4g0AkiRpCLIM/ieltIEjTFuR32cdA9xtmVJqBz6W/zlS+Y3Au7PUp5L1WS/TKTMkSdIQjIsZ9ceTmZNrqc7fidlqi5kkSRoCg1mJVVUFjQ3OZSZJkobOYDYKnP1fkiQNh8FsFDj7vyRJGg6D2ShoLrSYOV2GJEkaAoPZKChMmbG9rZPc+uySJEmDM5iNgsLg/87uXvZ2dJe5NpIkabwwmI2CPnOZ2Z0pSZIyMpiNguLZ/70zU5IkZWUwGwWFuzLBSWYlSVJ2BrNR0NxnWSaDmSRJysZgNgpmNTzTYtZiV6YkScrIYDYKaqurmDWlFrDFTJIkZWcwGyVNTjIrSZKGyGA2SgrBrLXNrkxJkpSNwWyUNLpepiRJGiKD2SixK1OSJA2VwWyUNOcnmW3r7KG9s6fMtZEkSeOBwWyUFE8ya3emJEnKwmA2SpqcZFaSJA2RwWyUNE51vUxJkjQ0BrNRYlemJEkaKoPZKOnTlemdmZIkKQOD2Sipr61m2qQawElmJUlSNgazUdSUnzKjxa5MSZKUgcFsFDU25Gf/tytTkiRlYDAbRQdn/7fFTJIkZWAwG0VN03ItZo4xkyRJWRjMRlGhxWzX/i66enrLXBtJklTpDGajqHiS2VYnmZUkSYMwmI2iZieZlSRJQ2AwG0XFk8w6ZYYkSRqMwWwUNdmVKUmShsBgNooKE8yCXZmSJGlwBrNR1FBXzaSa3Cl2kllJkjQYg9koiggnmZUkSZkZzEZZoTtzu2PMJEnSIAxmo6wwZYYtZpIkaTAGs1HW2GCLmSRJysZgNsoK62XuaOugpzeVuTaSJKmSGcxGWWHwf2+CnfttNZMkSYdnMBtlTjIrSZKyyhTMImJhRPwgInZHxJ6IuC4iFmUsWx8Rfx8RWyKiPSLujIgXHmbfBRFxdURsjYiOiFgbEX83lAOqNMXBzBsAJEnSkdQMtkNETAFuATqAtwMJ+Cxwa0SckVJqG+Qlvga8CvgL4EngQ8AvIuK5KaX7it5nCXAHsBb4U+BpYAlwwpCOqMI0uZC5JEnKaNBgBrwXWAqclFJaAxARDwCPA+8Hvni4ghFxJvBm4F0ppWvy224HVgNXApcW7f5lYBPw4pRSV37b7UM6mgrUZyFzZ/+XJElHkKUr81LgrkIoA0gprSXXunVZhrJdwLVFZbuB7wEXR8QkgIg4HrgY+KeiUDYhzJhcS01VAE6ZIUmSjixLMDsNeHCA7auBUzOUXZtS2j9A2Tqe6aY8P//YHhG/zI8v2xkR34yIxgx1rFhVVUFjvjuz1a5MSZJ0BFmC2Wxg5wDbdwCzRlC28DzA/Pzj1cBjwCuAvyQ3Nu0XETFgPSPifRGxMiJWtrS0DFKV8nlmklmDmSRJOrys02UMNDNqZCgXGcsW6nFbSulDKaVbUkpXAX8MnEOum/PQSqV0VUppRUppRXNzc4bqlIfrZUqSpCyyBLOdPNOyVWwWA7eGFdtxhLKF5wFa84+/7LffTfnH5YO8T0Vrcr1MSZKUQZZgtprcWLH+TgUeylD2uPyUG/3LdgJrivaDgVvXAHoz1LNiNefvzGzd10lKLsskSZIGliWY3QicFxFLCxvyc46dn39usLK1wBuKytYAbwJuSikVmpDuArYCL+9XvvD332WoZ8UqTJnR2dPLngPdZa6NJEmqVFmC2VeBdcANEXFZRFwK3ABsBL5S2CkiFkdEd0R8srAtP4HstcCXIuI9EfESclNlHAf8TdF+3cDHgVdFxJcj4mUR8cfAvwK3kZvgdtxqdJJZSZKUwaDBLD+z/4Xk7pb8FvAdcrPzX5hS2le0awDVA7zmO4FryK0W8FNgIfDylNKqfu/zDeBtwPOBHwN/DXwbuCSN8/6/PssyOcmsJEk6jCwz/5NS2gC8fpB91jHAnZoppXbgY/mfwd7nW+TC34TSd71M78yUJEkDyzpdhkagadozXZmtbbaYSZKkgRnMxsDsKXVEvi3RrkxJknQ4BrMxUFNdxawpuVazFrsyJUnSYRjMxoiTzEqSpMEYzMZI08FJZg1mkiRpYAazMVIIZt6VKUmSDsdgNkYa7cqUJEmDMJiNkUKL2f7OHvZ3uiyTJEk6lMFsjDT3mf3f7kxJknQog9kYKZ5kdruTzEqSpAEYzMaI62VKkqTBGMzGSKPrZUqSpEEYzMZIY0NRV6Z3ZkqSpAEYzMZIfW010+prACeZlSRJAzOYjaFmJ5mVJElHYDAbQ4VJZltsMZMkSQMwmI2hZ5ZlMphJkqRDGczG0MFg5nQZkiRpAAazMVQIZnsOdNPZ3Vvm2kiSpEpjMBtDxbP/tzr7vyRJ6sdgNoYaG1wvU5IkHZ7BbAw1T3OSWUmSdHgGszFUvF6mU2ZIkqT+DGZjqDiYtTrJrCRJ6sdgNoYaJtUwubYasCtTkiQdymA2xgqz/xvMJElSfwazMebs/5Ik6XAMZmOsEMwcYyZJkvozmI2xwpQZtphJkqT+DGZjrDDJ7I62Tnp6U5lrI0mSKonBbIw15Qf/96ZcOJMkSSowmI2xpmlFyzLZnSlJkooYzMaYk8xKkqTDMZiNseJgZouZJEkqZjAbY4UxZmAwkyRJfRnMxtiMybXUVgfgQuaSJKkvg9kYi4iDU2Y4xkySJBUzmJVBk5PMSpKkARjMyqDQYmYwkyRJxQxmZXBwIfO9dmVKkqRnGMzKoNCV2drWQUouyyRJknIMZmXQnG8x6+pJ7GnvLnNtJElSpcgUzCJiYUT8ICJ2R8SeiLguIhZlLFsfEX8fEVsioj0i7oyIFw5S5o8iIkXEU1neY7wpnmTWKTMkSVLBoMEsIqYAtwAnA28HLgeWAbdGREOG9/ga8F7gk8CrgS3ALyLirMO830zgH4CtWQ5gPGp0kllJkjSAmgz7vBdYCpyUUloDEBEPAI8D7we+eLiCEXEm8GbgXSmla/LbbgdWA1cClw5Q7P8A95MLcBdlPpJxxGWZJEnSQLJ0ZV4K3FUIZQAppbXAHcBlGcp2AdcWle0GvgdcHBGTineOiPOBtwIfylT7ccqFzCVJ0kCyBLPTgAcH2L4aODVD2bUppf0DlK0DTihsiIha4Crg74tD4EQ0u6GOyK3KZIuZJEk6KEswmw3sHGD7DmDWCMoWni/4S2AS8HcZ6gRARLwvIlZGxMqWlpasxcquuiqYPcXZ/yVJUl9Zp8sYaLKtyFAuspSNiBOAvwI+nFI6kLFOpJSuSimtSCmtaG5uzlqsIhS6M1ucZFaSJOVlCWY76duyVTCLgVvDiu04QtnC8wD/j9ydn3dFxMz8nZl1QOT/PjlDPccV18uUJEn9ZbkrczW5sWL9nQo8lKHsayNiSr9xZqcCncCaor8vZuCgtxP4R+AjGeo6bhRazFrbDGaSJCknS4vZjcB5EbG0sCEilgDn558brGwt8IaisjXAm4CbUkqFVPKHwIv7/fwC2J7/8z9nqOe44nqZkiSpvywtZl8FPgzcEBGfIDdm7DPARuArhZ0iYjHwBHBlSulKgJTSfRFxLfCl/F2Xa4EPAscBbymUTSnd1f9NI+IdQEdK6bZhHVmFK0wy297VQ1tHNw2TsvxTSJKkiWzQFrOUUhtwIfAY8C3gO+QC1oUppX1FuwZQPcBrvhO4Bvgs8FNgIfDylNKqEdd+HHOSWUmS1F+mZpqU0gbg9YPss44B7tRMKbUDH8v/ZJZSesdQ9h9vmvsEs04WN2ZZ3UqSJE1kWafLUInZYiZJkvozmJVJYboMgG17DWaSJMlgVjbNUycxNT/g/9ePjZ9VCyRJ0ugxmJVJTXUVr3jWXABufXQbu/Y7bYYkSUc7g1kZvfbsBQB09SR+8sCWMtdGkiSVm8GsjM47rpF5M+oBuP7eTWWujSRJKjeDWRlVVQWXnZVrNbtn/U42tO4fpIQkSZrIDGZl9rp8dybYaiZJ0tHOYFZmJx4zjVPnTQfg+nufIqVU5hpJkqRyMZhVgEKr2brW/dy3cVeZayNJksrFYFYBLj1zPlX5xazszpQk6ehlMKsAc6bXc/4JTQD8+P7NdHb3lrlGkiSpHAxmFaLQnblzfxe/ciUASZKOSgazCnHxaXOZUlcN2J0pSdLRymBWIabU1XDxabklmn758NPsbu8qc40kSdJYM5hVkNcuz3Vndnb38vMHXaJJkqSjjcGsgpx/QhPN0yYBcN0quzMlSTraGMwqSHVVcNmZ8wG4e+0OntrpEk2SJB1NDGYV5rVFSzTdcN/mMtZEkiSNNYNZhTl13nROOmYakLs70yWaJEk6ehjMKkxE8Jr8TQBrtu1j9eY9Za6RJEkaKwazCvSa5fOJ/BJN3gQgSdLRw2BWgebNmMxzlzYCcOP9m+nucYkmSZKOBgazClXozty+r4PfrNle5tpIkqSxYDCrUK941lwm1eT+eVyiSZKko4PBrEJNq6/lZfklmn6xeiv7OrrLXCNJkjTaDGYV7LXLc5PNHujq5RcPbi1zbSRJ0mgzmFWwFyxrprGhDrA7U5Kko4HBrILVVldxSX6Jpjue2M7W3QfKXCNJkjSaDGYV7rX5uzNTghvvt9VMkqSJzGBW4c44dgZLmxsAJ5uVJGmiM5hVuIjgdflWs0e27uXhLS7RJEnSRGUwGwcuO2vBwT//yJsAJEmasAxm48DC2VM4d8lsAH503yZ6elOZayRJkkaDwWyceO3ZuVazp/d0cOcTrWWujSRJGg0Gs3Hilc+aR121SzRJkjSRGczGiRlTannJKXMA+PmDW2jv7ClzjSRJUqkZzMaRwpxmbZ093PSQSzRJkjTRGMzGkRedNIeZU2oBuzMlSZqIDGbjSF1NFa8+Yx4Av358Oy17O8pcI0mSVEoGs3HmtcuPBaCnN3Hj/ZvLXBtJklRKBrNx5uxFM1ncOAVwsllJkiaaTMEsIhZGxA8iYndE7ImI6yJiUcay9RHx9xGxJSLaI+LOiHhhv31OjIh/jIgHImJfft8bI+LM4RzURBYRvCa/EsDvN+1mzba9Za6RJEkqlUGDWURMAW4BTgbeDlwOLANujYiGDO/xNeC9wCeBVwNbgF9ExFlF+7wMeDHwDeAS4I+BZuDuiDgn89EcJQp3Z4ILm0uSNJFkaTF7L7AUeE1K6UcppRuAS4HFwPuPVDDf4vVm4KMppa+mlG4G3ghsAK4s2vV7wJkppS+klG5NKV0PvBxoB/5sqAc10S1pamD5opkA3HDfZnpdokmSpAkhSzC7FLgrpbSmsCGltBa4A7gsQ9ku4Nqist3kgtjFETEpv217SqlPukgp7QYeAxagQ7wu32q2aVc7v123o8y1kSRJpZAlmJ0GPDjA9tXAqRnKrk0p7R+gbB1wwuEKRsRs4FnAwxnqeNR59RnzqakKAK63O1OSpAkhSzCbDewcYPsOYNYIyhaeP5x/AgL40uF2iIj3RcTKiFjZ0tIySFUmllkNdbzopNwSTT/7/RYOdLlEkyRJ413W6TIGGsQUGcrFcMpGxBXkxqZ9uLgL9ZBKpXRVSmlFSmlFc3NzhupMLK87O9edubejm5sf3lbm2kiSpJHKEsx2MnDL1iwGbg0rtuMIZQvP9xERHwD+FvhESunqDPU7al148hym1dcAcP29T5W5NpIkaaSyBLPV5MaK9Xcq8FCGssflp9zoX7YT6NMaFhGXA/8KfCGl9LkMdTuq1ddW86rTc0s03fZoCzvaOstcI0mSNBJZgtmNwHkRsbSwISKWAOfnnxusbC3whqKyNcCbgJtSSh1F218LXAP8e0rpzzPW/6hXmNOsuzfxkwdcokmSpPEsSzD7KrAOuCEiLouIS4EbgI3AVwo7RcTiiOiOiE8WtqWU7iM3VcaXIuI9EfESclNlHAf8TVHZFwL/ATwAfD0iziv6WT7io5zAnr1kNgtmTgacbFaSpPFu0GCWUmoDLiQ3p9i3gO8Aa4ELU0r7ihaDyCkAAB+aSURBVHYNoHqA13wnuZawzwI/BRYCL08prSra50JgErCc3Pxodxb9XD/kozqKVFUFr1k+H4D7Nu5i7fa2MtdIkiQNV/Sb13XcWrFiRVq5cmW5q1EWa7bt46Iv3g7An75kGR976YllrpEkSTqSiLgnpbSi//as02Wogp0wZypnHDsDgB/du4mJErYlSTraGMwmiNeclbsJYMOO/azaMNgsJpIkqRIZzCaIS8+aT3V+iSZvApAkaXwymE0QTVMn8cJlTQD85IEtdHS7RJMkSeONwWwCeU1+TrPd7V3c4hJNkiSNOwazCeRlp85l6qTcEk2f+NGDPNmyb5ASkiSpkhjMJpDJddVc8cqTAWht6+Tyr/2Wp/ccKHOtJElSVgazCeYtz1nMn71kGQCbdrXztq/9lt37u8pcK0mSlIXBbAL6yEXLuPy8xQA8+vRe3vWN39He6c0AkiRVOoPZBBQRfOrS03jVGfMAuGf9Tv74O/fQ1dNb5ppJkqQjMZhNUNVVwRffeCbPPyE3hcatj7bwv37wAL29rgogSVKlMphNYJNqqvny5edwZn65puvv3cTnfvawSzZJklShDGYT3NRJNVz9jmeztLkBgK/9Zi3/dvsTZa6VJEkaiMHsKNA4dRLfevdzmDu9HoD/8/NH+d5vN5S5VpIkqT+D2VFiwczJfOvd5zJjci0A//v63/PzB7eWuVaSJKmYwewosuyYaVz9jmczubaa3gR/+r17ufOJ1nJXS5Ik5RnMjjLnLJ7Fv731bGqqgs7uXt77zZU8uGl3uaslSZIwmB2VXnTSHL7wxjMB2NfRzTuu+S3rtreVuVaSJMlgdpS67KwF/M0lpwKwfV8nl199N9tcV1OSpLIymB3F3nn+cfzJhScAsHFHO2+7+rfsbnddTUmSysVgdpT72EtP5M3PWQTAI1v38h7X1ZQkqWwMZke5iOAzlz2LVzxrLgC/W7eTD393letqSpJUBgYzUV0VfOkPz+J5xzcCcPMj2/j4D3/vupqSJI0xg5mA3LqaX7n8HJ61YDoAP1z1FJ//+SNlrpUkSUcXg5kOmlZfy9ffeS7HNeXW1bzqV0/yZdfVlCRpzBjM1EfT1El8813ncsz0SQB8/r8e4fu/21jmWkmSdHQwmOkQC2dP4Zvveg7T62sA+Ph1D/DLh54uc60kSZr4DGYa0Elzp3HNO59NfW0VvQk+9N1VXPu7DRzocioNSZJGi8FMh3XO4tn821vOoTq/ruZf/vD3nP/5W/jiTY/ytKsESJJUcgYzHdGLT57Dv7x5Oc3TcmPOWts6+X+3rOH8z9/Cn/7Hvdy7YWeZayhJ0sQRKU2MuapWrFiRVq5cWe5qTFid3b3814NbuOaOddy3cVef585cOJN3Pm8Jrzx9HnU1Zn1JkgYTEfeklFYcst1gpqG6d8NOvv4/6/jpA1voLpqEtnnaJN7ynEW85TmLD7awSZKkQxnMVHJP7znAd+7ewHfvXs/2fZ0Ht9dVV/HqM+bxjvOXcMaxM8tYQ0mSKpPBTKOmo7uHn9y/hWv+Zy0PbtrT57lzFs/iHc9bwsufNZfaars5JUkCg5nGQEqJe9bv5Jr/WcfPH9xKT1E35zHTJ3H5eYv5o3MX0TjVbk5J0tHNYKYxtWV3O9++az3fvXsDO/d3HdxeV1PFZWfO5x3nL+G0+TPKWENJksrHYKayONDVw433bebqO9byyNa9fZ4797jZvOU5i3jZqXOZXFddphpKkjT2DGYqq5QSd6/dwdfvWMdND22lqJeThrpqXnH6PF63fAHnLW2kqirKV1FJksaAwUwV46md+/nWXeu59ncb2VXUzQkwf0Y9ly1fwOuWL2DZMdPKVENJkkaXwUwVp7O7l9se3cZ1qzZxyyPb6Ozp7fP86Qtm8LqzF3DJmfNp8oYBSdIEYjBTRdu1v5OfPLCF61Y9xaoNfVcWqKkKLjixmdedfSwvOWUO9bWOR5MkjW8GM40ba7e3cf29m7j+3qfYuKO9z3PT6mt41enzeN3Zx7Ji8SzHo0nSCO3e38V9T+1i1fqd3LtxFw9u2s2caZP4wAXHc8mZ86n2c3ZUjCiYRcRC4B+AlwIB/DfwkZTShgxl64HPAG8FZgL3AX+ZUvpVv/2qgL8E3g/MBR4Frkwp/XDQCmIwm4hSSqxcv5PrVm3iJw9sZu+B7j7PL5w9mdeetYDXnn0sxzU1lKmWlaGto5ufP7iVLbvbOWHONE6dN51jZ002uEpltq+jm/WtbWxo3c/m3QdomlrHcU0NLGlqYHp97ZjXp6c38fi2vdy74ZkgtmbbvsPuf3xzAx+56ERedfo8P09KbNjBLCKmAPcDHcAngAR8FpgCnJFSahuk/HeAVwF/ATwJfAh4BfDclNJ9Rft9Dvhz4K+Ae4A/BN4LvDql9LPBDtBgNrEd6Orhlke2cd2qp7jt0ZY+a3QCLF80k9edfSyXnDGPmVPqylTLsVWY0Pf7Kzfykwe2sL+zp8/zUyfVcPLcaZw8bxqnzJvOKfOmc/LcaUypqylTjaWJJ6VEy74ONrTuZ33rftbv2M+G1rb8435a2zoPW7YQ0nI/UzmuqYGlzQ0smj2lZEM2drZ1ct/GXazasJN7N+zivo272NfRfdj9j29u4PQFM/jNmla27+s4uP2kY6bx0Zcu42WnzjWglchIgtmfAV8ETkoprclvOw54HPhfKaUvHqHsmeRayN6VUromv60GWA08mlK6NL9tDrAR+HxK6W+Kyt8MNKeUzhjsAA1mR4/WfR38+P7NXH/vJu5/anef52qrg+ULZ3Hq/Om5n3nTWXbMVCbVTJxxadv2HuC6VZv4/sqNPNlyxN+LDhEBSxobOGXeNE6Zmwtrp8yfzvwZ9UT4YSsNpKunl8272vsGr9b9bNiR++n/S9FIRcCCmZNzQa0Q3JqnsrSpgfkzJx+2a7G7p5fHnt53MITdu2EnT24//GfEtEk1nLVoJssXzeLsRTM5a+HMg7/Ytnf28O271vPl25/oEy5PnTedj770RC46Zc6E+8zo6U388qGtfPn2J/n860/n5LnTR/X9RhLMbgbqU0rn99t+O0BK6YIjlP1r4K+BmSml/UXbPw18HJieUuqIiMuBbwInppQeL9rvncDVwNKU0toj1dNgdnRas20f19/7FNev2sTm3QcG3KemKjhhztSDQa3wOJ5a1rp6ern1kW18f+VGbn20pc9yVzVVwYUnz+GNKxby7CWzeXzbXh7eupeHt+zh4S17eHTr3kG/OKbX1xxsVTsl38J24jHTvNFCR53d7V38+vEW7n5yB+vyAWzTrvY+/+cGM3d6PYsap7B49hQWN05hUWMDi2dPYf7MybS2dbC2pY0nt7exNv+zbnvbEVvWitVVV7GoccrB0LaocQqbd7Wzav0u7n9q1xH/ry+bM5WzF83i7MW5MHZC89RBW7/aOrr55p3r+cqvnugzvdEZx87goy89kRed2DzuA1pHdw/Xr9rEVb968mCQfd3yBXzxTWeN6vuOJJhtBW5IKb2/3/Z/Bd6QUmo+QtnvActTSif12/5G4FrgWSml1RHxeeAjwORUVKGIOBe4m1x35k+PVE+D2dGttzc3ge1Pf7+Z3z+1m0e27qWju/eIZRbMnMwp86YVhbUZLJw9uaI+ZNZs28v3Vz7FdaueYvu+vh/cJ8yZyptWLOQ1yxfQPO3w04n09ibW79h/MKjlfvayaVf7YcsAVAUsbZ7KafOn89yljTx/WRPHzppSkuMqlU272vnVYy3c/mgL927cSW11FdPra5kxuZbpk2uYXl/L9Mm1+ceaouf6Pt9QVz3kf/fe3sS+zm72Huhm74GuPo97DnSzp71rwOf2HugmkaipqqK2poraqqC2uoqa6qAu/1hTXZX7c1X026eK2up+++f3mV5fy4JZk1kwczLNUyfZ3ZRRSoknWvZxyyPbuPnhbaxcv3PQEFZXXcWxsyfng1eu63FRPoQtHGY35O79XaxtbWPt9n2HBLfhtMhNr6/Jt4TNYvmimZy5cCYzJg9/TNveA118/Y51fPXXT7KnaLzv8kUz+dhLT+T5JzRV1GdnFnsOdPHduzdw9W/Wsm3vM922MybX8u7nH8efvmTZqL7/SIJZJ/DFlNLH+23/LPDxlNJhB6xExE3kWsXO67f9IuCXwAtTSr+OiKuAS1NKc/vtdwK5LtO3pZS+NcDrvw94H8CiRYvOWb9+/RGPRUeP7p5e1m5v46Ete3ho856Dj4P9VjptUq7lqLh1bay7Qvd1dPOT+zfz/ZUbD5k6ZOqkGi45cx5vWLGQ5QtnjuiDcPf+Lh7emgtqj2zZy8Nbc61rRwq0xzU18PwTmnj+siaee3zjmA9ePtDVw2/X7uD2x1q4/bGWIw5aHorqqmBafU0utBWFuGn1NXR09x4SrPa0d7Gvs5tKvam9rrqKeTPrWTAzF9QWzJrM/JmTOTb/53kzJlNXU1XuapbNga4e7l67g1sefppbHt12yN3fAJNrqzl+TgOLZzccbP1a1JgLYnOn14/ZnYopJbbt7eDJlkJQ28fa7bngtnHHfrp6EhG5MWDLD3ZLzmJpU8OohPPd7V1c/Zu1XP2btewtGqv27CWz+OhLT+R5xzeV/D1LbdueA1x9xzq+c9f6Pscwf0Y9737BUv7w2QtpmDT6Y3FHGsy+kFK6ot/2z5G7u/JIweyXwNSU0nP7bX8pcBPPBLOvkmsVm9dvv2XAYxwmmBWzxUyDKXzAHQxqW/bw8OY9rG1tO+IXbE1VsHD2FBbMnMz8mfXMn5n7kluQf5w3o37EXX4pJX63LjeQ/6cPbKG9q+9vyOceN5s3rljIK0+fO6qD97t7elnX2sZDW57pCr13wy52t3cdsm9VwJkLZ/KCE5p4/rJmli+aSW11ab/sU0qs3d52MIjd9WQrB7oODY6NDXU89/hGaqur2NPexZ4DXexp784/dtFW4jFAWdXkA9+0fMibXl/L1PoaqiPo6umlqzfR1d1Ld28vXT2Jrp5euvOPXb1Ffy56rv9EzMMRAXOmTTp4DS+Y9UxoWzBzCgtmTWbqGHwxjaWtuw9wyyPbuOWRbdyxZvsh/8cAFs2ewoUnz+HCk+fwnKWzK35sandPL1t2H2DmlFqmjfEvSbv2d/LVXz/JNXes69Oi99yljXzsZSfy7CWzx7Q+Wazd3sZVv3qCH96zqc//oxOPmcr7X3g8l541v+SfYUcykmD2NPCjYXZlXguclaEr8/8D/gy7MlUGbR3dPLJ1b5/WtUe27Bm0K7RY09S6PmEt9+dnQlxjQ92ArVtP7znAD1c9xX+ufIq1/QbpHjN9Eq8/+1jesGJhWacD6elNrN68m18/vp3fPL6de9bvHDAcNNRVc97SRs4/oYkXLGvihDlTh9Wit6+jmzufaOX2x7Zx+2MtA7ZmVFcFZy+ayQUnNnPBiXM4bf70I7YOdPX0Hmzp6h/adg+wbU/RvpNqqvPhqm/AKn6cdvAx9+fp+cf62qqSd++klOjpTXT35kJaIby17utk0652Nu3MTcuwaWc7T+1qZ9PO9j5312U1Y3Itc6fXM7uh7uDPrIY6Zk+pZfbUScyeUry9tuJCTE9v4r6Nu7j1kW3c/Mg2Ht6y55B9aqqCZy+ZzYUnz+HFJ8/h+OaGcdcdV2472jr5yq+e4Jv/s75P2H3BsiY++tITOXvRrDLWLuf+jbv48u1P8PPVW/v8Ev7sJbP4wAXH8+KT5pSl638kwewWoC6l9Px+22/Llz/S4P9Pkptio//g/08BV/DM4P+3Ad8AlhXu/Mzv9w7gGhz8rzFWaDlanQ9qG3fsZ9OuA2ze1U7L3qF/yU2qqcqHtFz30twZk3lw025ue3RbnwXda6uDi045hjeuWMgLljVRM4a/vWXV3tnDb9ft4DePt/Drx7fzyNa9A+43d3r9wZD2vBMamTOtfsD9Uko8snVvrlXs0RZWrt9BV8+hn0vzZtTng1gzzzuhaUTjZY42B7p62JIPa5t27T8Y2jbvamfTrna27DpwyBQ0QzV1Ug2zGmoPBrZciKtj9tTc46yGOhob6pgxuZbJddXU11YzuTb3WKpuwd3tXfzqsRZueSQX6ncMMHShsaGOF52UaxV7wYlNZZlLbCJq2dvBl29/gm/ftb7PL7UvPqmZj770RM44duaY1ielxK8e386Xb3uCO59s7fPcRafM4QMXHM+KMrfqjSSYfQT4v+TumHwyv20JubFfH08pfeEIZc8C7gXekVL6Rn5bDfB7YE1K6ZL8tsJ0GX+bUvp0Ufn/Bo5JKZ0+2AEazDRWOrp72Lr7QL51op3N+cC2eXfuS27zrvYBu9uO5MRjpvLGFQt57fIFNI6zdUFb9nZwx5rt/GZNrkVt656B7449ee60g+PTTp47nZXrd3D7o7kuym0DhN266irOPW52Loyd1MyyYbbAaXA9vYltewvBrf3gtb1tbwc72jrZ2dbJjv2dfe7KK6W66irqa6sOCWz1tVVMrq3Oba+ppr6u8FxV0T7V7DnQxW2PtnDPYQbuP2vBdC48+RguPHkOZyyY4Y0Ro+jpPQf4t9ue4Lt3b+jTsn7RKcfwnhccx/HNU2maOnAPQil09/Tyswe38uXbnuCholbSmqrgsrMW8P4LlnLiMdNG5b2HaiTBrIHcBLPtPDPB7GeAaeQmmN2X328x8AS52fqvLCr/PeBichPMrgU+CLwaeF5KaVXRfoU7M/83sAp4E7lVAC5LKf14sAM0mKlSpJTYub/r4Jfc5sLP7vY+rW7TJtVwyVnzedOKhZxx7IwJEToKd7gVuj3verJ1SOO7ljROORjEzlva6GS4Faa7p5dd7V3saOs8GNhaix/3dx7y3FCGBJRKQ101z1/WxIUnz+FFJ83hmOkDt9Zq9Gze1c6/3LqG76/ceEgLeF1NFfNn1OduSpnxzDjHkYzbPdDVw3+u3MhVv36yz/CHKXXV/NG5i3j3849j/szJJTm2UhnpkkyL6Lsk083klmRaV7TPEnLB69MppU8VbZ8MfA54M7klme4nd9PAbf3eo5pc9+Z76bsk0w+yHKDBTONJR3cPNVVVE34Nuq6eXu7buCsf1Fq4/6ndfVo0JtdW87zjG7ngpGZeuKyZJUf50loT0f7O7oNhbUdbJ7vbu+jo6qW9q4f2rh4OFB47ezjQb/vB57p6ae/soaO7h/bO3Lb+DWNLGqfw4vzA/XOPq/yB+0eLjTv28y+3ruE/73lqSHPBNU2ddHCcbvHY3WPzdxjPmlJLRLBrfyffunM9X/+fdX3uum9sqOMdz1vC5c9dXLFzVrqIuaSy23OgizufaOXJljbOOHYGK5bM8gtUQ5ZSoqsnHQxwAcyxVayibdndzu+f2l3Ui3Dg4DjH4Yzbra/NjdvduvtAn7tCF86ezPtesJQ/OGchk+sq+7PlcMHMfgJJY2Z6fS0XnzZ38B2lI4gI6mqCupoqbwIZJ+bNyM2fN5CO7h625Id5FMY4FsJb4e+d/brED3T19lmS7pR50/ngi47nlc+aW5E3TQ2FwUySJJXNpJpqljQ1HHYoQ0qJ1rbO/M1WfcNbTVUVb3z2Ql64bPytPHA4BjNJklSxIoKmqZNomjqJMxeO7bQb5TC+2/skSZImEIOZJElShTCYSZIkVQiDmSRJUoUwmEmSJFUIg5kkSVKFMJhJkiRVCIOZJElShTCYSZIkVQiDmSRJUoUwmEmSJFUIg5kkSVKFMJhJkiRVCIOZJElShTCYSZIkVYhIKZW7DiURES3A+nLXYxQ0AdvLXYky8xx4Dgo8D54D8BwUeB7G9zlYnFJq7r9xwgSziSoiVqaUVpS7HuXkOfAcFHgePAfgOSjwPEzMc2BXpiRJUoUwmEmSJFUIg1nlu6rcFagAngPPQYHnwXMAnoMCz8MEPAeOMZMkSaoQtphJkiRVCIOZJElShTCYjbGI+IOI+GFErI+I9oh4NCL+LiKmZSibDvNz1ljUvZQi4kWHOZZdGcrWR8TfR8SW/Dm8MyJeOBb1LqWIuO0I/6Y/H6TsuLwWIuLYiPin/L/Z/nydlwyw36yI+PeI2B4RbRHx3xFxesb3qIqIKyJiXUQciIj7I+L1pT6W4cpyDiLiJRHx7Yh4In+NPxER/xYRczK+x7rDXB+vGY1jGqqM52DJEa7zmRneo6KvA8h8Hr5+hPPwSIb3qPRrIdN34kT+TOivptwVOAr9ObAB+N/AU8By4FPAiyPieSml3kHKfx34Sr9tj5W4jmPpT4HfFf29O0OZrwGvAv4CeBL4EPCLiHhuSum+0ldx1PwxML3ftucCXwRuzFD+64y/a+EE4I3APcCvgZf13yEigtzxHwf8CbATuAK4NSLOSik9Nch7fIbc/7O/yr/PHwL/GRGvTin9rFQHMgKDngPgA8BU4LPkrvFlwKeBiyPijJTSvgzv8wtyny3FHh1mnUstyzko+DsO/f+wN8N7VPp1ANnOw2eAL/fbtgT4D7J9TkBlXwuDficeBZ8JfaWU/BnDH6B5gG1vAxJw4SBlE/DZch9Dic7Di/LHc9EQy52ZL/fOom015D5kbiz3cZXgvHwN6ABmT8RrAagq+vN78sexpN8+l+W3v7ho2wxgB/D/Bnn9Ofnz9+l+228GHij38Q/hHAz0OfHC/L7vyvAe64Bvl/tYR3gOluS3v2cYr1/x10HW83CYcn+d3/e0CXAtDPqdONE/E/r/2JU5xlJKLQNsLrQYLRjLuoxTlwJdwLWFDSmlbuB75FoTJpWrYiMVEZOBNwA/TintKHd9RkMavEUYcv/Gm1NKtxaV2w38mNwH9JFcDNQB3+63/dvA6RFx3BCqOyqynIOJ/jmR8ToYiYq/DmBE5+FtwD0ppdWlrE85ZLzWJ/RnQn8Gs8pwQf7x4Qz7fjAiOvLjEW6JiBeMZsXGwHcioiciWiPiuxGxaJD9TwPWppT299u+mtx/vhNGpZZj43XANOAbGfefaNdCwWnAgwNsXw0sioipg5TtANYMUBbg1JFXr2yG8jkBcEn+2uiIiLsqZUzRMPxdRHRHxO6IuDHjuKIJex1ExPnkPueyfk7A+LsW+l/rR9VngsGszCJiAXAl8N8ppZWD7P5tcuOSLgLeBzQCt0TEi0a1kqNjN/AFcs33F5IbA3ARcOcgA5xnkxtf0N+OoufHq7cB24D/yrDvRLoW+hvs33jWIGV3pXxfxQBlx+X1kR8I/SVyX1Q/ylDkx+TG4lwMvAU4AFwfEW8dtUqWXge5MZTvB15MbozQ6cD/RMQpg5SdkNdB3tvI9Rr8R8b9x9W1cJjvxKPqM8HB/2WUT/k3kBvw/s7B9k8pXV70119HxA3kfov4LPD8UankKEkp3QvcW7Tp9oj4FfBbcjcEfOIwRYPcWIOBto9bETGfXMj6x3zX7BFNpGthACP5N55w10dE1JD7El4AnJ/x+viTfq9xPXAXuYH0/bt0KlJKaQu5myAKfh25u5VXkxvEfaRgMeGuA4D8UI03Aj9JKW3PUmY8XQtH+E48qj4TbDErk4ioJ3eXyVLg4jT4XSWHSCntBX4KPLvE1SuLlNIqcncVHul4djDwbzizip4fj95K7v/jULonDppg18Jg/8YD/eZcXHZW/i6ugcqOq+sjIgrXxEXAa1JKDwzndVJKPcB/AsdGxLwSVnFMpZQ2Ar9h8Ot8Ql0HRS4DZjLMzwmo3GthkO/Eo+ozwWBWBhFRC/wQOBd4ZUrp9yN5OQb+bWC8Gux4VgPHRcSUfttPBTo5dBzBePE24P6U0v0jeI2Jci2sJjcupL9TgQ3pyFNFrAYmAccPUBbgoZFXb0x9GXgT8IcppZtH+FqFL6bxfo1kuc4n2nVQ8HZgOzDSKR4q6lrI8J14VH0mGMzGWP434O8ALwEuSyndNYLXmk5uPq+7S1S9soqIFcCJHPl4bgRqyd29WChXQ+7L66aUUseoVnIU5I/7NEbwW/AEuxZuBBZERGEAcOH4LmHweZt+Ti6gv6Xf9rcCD6aU1payoqMpIgpjMN+ZUsoyruxIr1VD7v/MhpTS1lLUrxzyNwedz+DX+YS5Dgoi4hhy85x9N6XUNYLXqahrIeN34lH1meAYs7H3L+T+U3wOaIuI84qeeyql9FRELAaeAK5MKV0JEBF/DpwE3ApsBhaTGww7l0MvuIoXEd8B1gKrgF3kJhW8AtgE/FN+n0POQ0rpvoi4FvhS/restcAHyU08OO7OQ97byI2p+G7/JybitRARf5D/4zn5x1dERAvQklK6ndwH7Z3AtyPiL3hmMskA/k+/1+oGvpFSejdASmlbRPwDcEVE7CV3fb2J3A0mg91WP2YGOwcR8ZfAx4Crgcf7fU60pJSeKHqtPucgIv6I3LH+DNgIHENuEuZzgD8axcMakgzn4AvkGg/uBFrIXfNXAL3A3/Z7rXF5HUCm/w8FbyH3nX3YX+DG6bUw6HciR8FnQh/lnkjtaPshN9lfOszPp/L7LCn+e37bJcAd5Jqxu4BWchfrueU+pmGehyuAB8jdndlF7kPjKmBe0T6HnIf89snkZsffSu4Oo7uBF5X7mIZ5HmrJfen8+DDPT7hr4QjX/21F+8wmF0p2APvJTQZ55mFe6+v9tlWTu3lkPbk7+x4A/qDcxz2UcwDcdoR9+h9vn23AecAtwNP562M38N/kxu2U/diHcA7eRW4+q53kfnHZSu6Xl5MmynWQ5TwU7Xc/8PsMrzWurgUyfCfm95vQnwnFP5GvtCRJksrMMWaSJEkVwmAmSZJUIQxmkiRJFcJgJkmSVCEMZpIkSRXCYCZJ/3+7dSwAAAAAMMjfeho7iiKACTEDAJgQMwCAiQAY0pyCdyEj7gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Get radius data\n",
    "r = df['r']\n",
    "\n",
    "# Get correlations for T=4\n",
    "c = df[corrTs[4]]\n",
    "\n",
    "# PLot the correlation vs the radius\n",
    "plt.plot(r, c);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the function `optimize.curve_fit()` from `scipy` module, fit the correlation verus the radius for each temperature, by a decreasing exponential function:\n",
    "$$\n",
    "f(r) = A \\; e^{-\\frac{r-r_0}{\\xi}}\n",
    "$$\n",
    "where $A$, $r_0$ and $\\xi$ are parameters to be fitted. \n",
    "\n",
    "**HINT:** When using `optimize.curve_fit()` function, it can be useful to give initial value of free parameters (argument `p0`), and also to restrict the parameters ranges during the minimization (argument `bounds`). \n",
    "\n",
    "Plot $\\xi$ vesus $T$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}