Phase_Transition.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0dd8f032-65ec-484c-b668-8291b9ef6547",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import curve_fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "187759c4-e738-41fe-9fa1-121805697a45",
   "metadata": {},
   "outputs": [],
   "source": [
    "N_avg = np.genfromtxt(\"Functions/Result/N_avg.csv\", delimiter = ',')/100\n",
    "N_std = np.genfromtxt(\"Functions/Result/N_std.csv\", delimiter = ',')/100\n",
    "N_std[N_std == 0] = 1e-5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0aba9975-d770-433c-971e-98b51af833d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def func(x, a, b):\n",
    "    return a*x + b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a64e42e5-c0d9-468f-ba7c-da538844ca16",
   "metadata": {},
   "outputs": [],
   "source": [
    "xdata = np.linspace(1,50,50)/100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b4d6da3b-4e0f-4b6e-906c-cdbe3e9105bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = []\n",
    "b = []\n",
    "ea = []\n",
    "eb = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "11e82829-3e79-465f-9f77-ee26158d1ca3",
   "metadata": {},
   "outputs": [],
   "source": [
    "ab, pcov = curve_fit(func, xdata[1:15], N_avg[1:15], sigma = N_std[1:15])\n",
    "a.append(ab[0])\n",
    "b.append(ab[1])\n",
    "ea.append(np.sqrt(pcov[0,0]))\n",
    "eb.append(np.sqrt(pcov[1,1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "de1f4148-c682-4c51-aa3e-1ba9d193a4dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "ab, pcov = curve_fit(func, xdata[14:22], N_avg[14:22], sigma = N_std[14:22])\n",
    "a.append(ab[0])\n",
    "b.append(ab[1])\n",
    "ea.append(np.sqrt(pcov[0,0]))\n",
    "eb.append(np.sqrt(pcov[1,1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "97952723-e59c-41a4-b978-31ecaaea49a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "ab, pcov = curve_fit(func, xdata[21:25], N_avg[21:25], sigma = N_std[21:25])\n",
    "a.append(ab[0])\n",
    "b.append(ab[1])\n",
    "ea.append(np.sqrt(pcov[0,0]))\n",
    "eb.append(np.sqrt(pcov[1,1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "8105f733-f509-4970-98e1-07479a9ab85a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGBCAYAAACOxg66AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeXxU1fnH8c8DCigJGmQRIbRqEX7RqLUuLKJQpVURbN1aFFsVRKnSCljErQaslaq4VSyiYkGxohQV64IoBBCCClaNRlEQWYobGjTDHnh+f8xAh5DATGYmd5L5vl+v+5qZe++595njCA/nnHuOuTsiIiIimaRe0AGIiIiI1DQlQCIiIpJxlACJiIhIxlECJCIiIhlHCZCIiIhkHCVAIiIiknH2CjqAdLL//vv7j370o6DDqNPWrVtH48aNgw6jTlMdp57qOLVUv6mXKXW8aNGiNe7evLJjSoCitGzZkoULFwYdRp1WWFhIt27dgg6jTlMdp57qOLVUv6mXKXVsZsurOqYuMBEREck4SoBEREQk4ygBEhERkYyjBEhEREQyjhIgERERyThKgERERCTjKAESERGRjKMESERERDKOEiARERHJOEqAREREJOMoARIREZG4hUIhioqKCIVCQYdSLUqAREREBICCggLMLKYtOzubzp07k52dHXOZgoKCoL/iDkqAREREJOMoARIREREg3ALk7nvcysrKyM3NJSsri9zcXMrKymIqV1kLUFBdaXvV6N1EREQkUKFQiA8++IBjjz2WrKysal0jKyuLkpISiouLyc/Pr/Z1QqEQeXl5lJaWkpOTQ0lJSbWvFS+1AImIiGSI7QnHH//4R/Ly8hJqdcnKyqJTp06VJizdunWLeRzRypUrCYVCrFy5MubxRN26dUugFsKUAImIiNRysQ5e3p5wbNiwIa6EI97By4WFhSntSissLEy4ztQFJiIiIoFIVldadagFSEREpJaLdfDyunXrOOiggxg8eHDMrS1VDV5Olt11paWSEiAREZEMMXr0aLp06ULv3r1rPOFIN+oCExERyQCrV6/mnnvuYeHChSxfvjzocAKnFiAREZEMcMMNN3DZZZdx8MEHBx1KWkjrBMjMTjOzxWa2xMyGV3L8LDN7z8zeMbOFZnZirGVFREQyxdtvv83LL7/M9ddfH3QoaSNtEyAzqw+MAU4H8oA+ZpZX4bTXgKPc/WjgUuDhOMqKiIjUee7O4MGDGTFiBE2aNAk6nLSRtgkQcDywxN0/dffNwJPAWdEnuHvI3T3ysTHgsZYVERHJBM888wylpaX069cv6FDSSjoPgm4NrIz6vAo4oeJJZvZL4DagBdAznrKR8gOAAQDNmzdPyuRKUrVQKKQ6TjHVceqpjlNL9Zs8mzdvZtCgQQwePJi5c+fu2K86Tu8EyCrZ57vscH8GeMbMTgJuAU6NtWyk/DhgHED79u09GdNrS9UKCwuTMoW5VE11nHqq49RS/SbP6NGj+fGPf8w111yz037VcXonQKuA3KjPbYDVVZ3s7nPM7FAzaxZvWRERkbrm66+/ZtSoUTu1/Mj/pPMYoLeAdmZ2sJk1AH4NTIs+wcx+ZGYWeX8M0AD4JpayIiIidVlBQQF9+vShQ4cOQYeSltK2Bcjdy83sKmA6UB8Y7+4fmNkVkeNjgXOA35jZFmAD8KvIoOhKywbyRURERGpYSUkJTz/9NB9++GHQoaSttE2AANz9ReDFCvvGRr3/K/DXWMuKiIhkgqFDh3L99ddzwAEHBB1K2krnLjARERGJ08svv8zSpUv53e9+F3QoaU0JkIiISB1RXl7O0KFDueOOO2jQoEHQ4aQ1JUAiIiJ1xLhx42jZsiW9e/cOOpS0l9ZjgERERCQ2a9euZcSIEUyfPp3IA9KyG2oBEhERqQNuvfVWevXqxdFHHx10KLWCWoBERERquaVLl/Loo4/y/vvvBx1KraEWIBERkVpu2LBhDBkyhAMPPDDoUGoNtQCJiIjUYrNnz2bRokU8/vjjQYdSq6gFSEREpJbatm0bQ4YMYdSoUeyzzz5Bh1OrKAESERGppSZOnEiDBg341a9+FXQotY66wERERGqhdevWccMNN/Cvf/1Lj71Xg1qAREREaqHbb7+dk08+mY4dOwYdSq2kFiAREZFaZtWqVdx///385z//CTqUWkstQCIiIrXMddddx8CBA2nbtm3QodRaagESERGpJUKhEJMnT+bVV1/lk08+CTqcWk0tQCIiIgEqKCjAzGLasrOz6d+/P1988QXZ2dkxlSkoKAj6K6YlJUAiIiIBKigowN33uM2fP59GjRoBkJWVxfz582MqpwSockqAREREAhRrC1Dnzp3ZuHEjEO4K69y5s1qAEqAESEREJECxtgCNGjWKnj17Mn/+fMrKymIqoxagqmkQtIiISJr78ssvueOOOygqKqJdu3ZBh1MnqAVIREQkzf3pT3/iN7/5jZKfJFILkIiISBorLi7m2Wef5aOPPgo6lDpFLUAiIiJpyt0ZMmQIN910Ezk5OUGHU6coARIREUlTL774IqtWreLyyy8POpQ6R11gIiIiaWjLli0MHTqU0aNHs/feewcdTp2jFiAREZE0NHbsWNq2bcsZZ5wRdCh1klqARERE0kxpaSm33HILM2fOxMyCDqdOUguQiIhImhk5ciRnn302RxxxRNCh1FlqARIREUkjH3/8MY899hglJSVBh1KnqQVIREQkjQwbNoxhw4bRokWLoEOp09QCJCIikiZmzZrFe++9x5NPPhl0KHWeWoBERETSwNatWxk8eDB//etfadSoUdDh1HlKgERERNLAP/7xD7Kzszn33HODDiUjqAtMREQkYGVlZdx0000899xzeuy9hqgFSEREJGCjRo3i1FNP5bjjjgs6lIyhFiAREZEALV++nLFjx/Luu+8GHUpGUQuQiIhIgK677joGDRpEmzZtgg4lo6gFSEREJCALFixg7ty5PPTQQ0GHknHUAiQiIhIAd2fw4MHceuutNG7cOOhwMo4SIBERkQA8+eSTbNmyhb59+wYdSkZSF5iIiEgN27BhA8OHD+fxxx+nXj21RQQhrWvdzE4zs8VmtsTMhldy/EIzey+yzTezo6KOfWZmxWb2jpktrNnIRUREqnbXXXdx3HHH0bVr16BDyVhp2wJkZvWBMUAPYBXwlplNc/fo5XGXASe7e6mZnQ6MA06IOt7d3dfUWNAiIiJ78MUXX3D33Xfz5ptvBh1KRkvnFqDjgSXu/qm7bwaeBM6KPsHd57t7aeTjAkDPEIqISFq78cYbufTSSznkkEOCDiWjpW0LENAaWBn1eRU7t+5U1A94KeqzA6+YmQMPuvu4ygqZ2QBgAEDz5s0pLCxMJGbZg1AopDpOMdVx6qmOU6su1++SJUuYOnUqEydODPQ71uU6jlU6J0CVLYbilZ5o1p1wAnRi1O4u7r7azFoAM8zsI3efs8sFw4nROID27dt7t27dEg5cqlZYWIjqOLVUx6mnOk6tulq/7s7IkSP5y1/+wplnnhloLHW1juMRcwJkZn9K8F4T3f2zOM5fBeRGfW4DrK4kriOBh4HT3f2b7fvdfXXk9Ssze4Zwl9ouCZCIiEhNmDZtGl999RX9+/cPOhQhvhagAsItMNVZptaB14HP4ijzFtDOzA4G/gv8Grgg+gQzawtMBS5y94+j9jcG6rl7WeT9z4CR1YhbREQkYZs3b+aaa67h/vvvZ6+90rnzJXPE+19hMPBcnGWaAoviLIO7l5vZVcB0oD4w3t0/MLMrIsfHAn8CDgAeMDOAcnc/FmgJPBPZtxfwhLu/HG8MIiIiyfDAAw/Qrl07fv7znwcdikTEmwCtcffl8RQws1Cc99jB3V8EXqywb2zU+/7ALm2J7v4pcFTF/SIiIjXtm2++4S9/+UvGDzpON/EkQJ2AJdW4x9pI2ZI9nSgiIlLXjBgxgvPPP5+8vLygQ5EoMc8D5O5vRA8yjqPc1kjZsnjLioikky+//JK+ffuybt26oEORWuKjjz7in//8JwUFBUGHIhVoJJaISAxKSkro2bMn69ev54MPPuDpp5/mBz/4AW3btt3xmpubyz777BN0qJJGrrnmGq677jqaNWsWdChSgRIgEZE9mDlzJn369OH222+nU6dOPP/88+y///6sWLGCuXPnMmnSJFasWMGqVavYb7/9aNu27U6J0fbXtm3b0qxZM9atW0dxcTH5+flkZWUF/fUkRWbMmMHixYuZOnVq0KFIJWJKgMzsPKAL8C7wmLuXRx17wd17pig+EZFATZgwgWHDhjF58uQdE8f95Cc/qXQSuW3btvHll1+yfPlyVqxYwfLly1m6dCmzZs3asW/Dhg2Ul5fj7jRq1Iirr76aww47bEei1Lp1axo0aJBw3AUFBYwYMSLh61Tl5ptvVrfObpSXlzNkyBDuuOMOGjRoQCgUUtKbbtx9txvwB+AL4BFgMfAG0DTq+Pd7ukZt2Q477DCX1Jo1a1bQIdR5quPk2LZtm990001+yCGH+IcffrjTscrq+OSTT3bCc56lZDviiCP83Xff9dLS0qR+z7KyMs/NzXXAc3NzvaysLKnXr47a8hu++eabU/rf/Oabb05Z7LWljhMFLPQq/s6PpQXod8DP3f3dyArtfwNmmtlP3f1bqjcxoohI2tq0aROXXnopn376KUVFRbRo0WKPZWJ9xDkUCpGXl0dpaSk5OTmUlJTs1CJQXl7O6tWrd7QgRb9ecMEFLF++nHr16u3StRbd3daqVSvq168fUzzFxcWUlobXlC4tLaW4uJhOnTrFVFakNoslAWrl7u9C+Iku4HdmNhqYZWY/pYr1uUREaovddRe1bNky4etHdxdlZWVRUlJSZXfIXnvttSOpOfHEE3e5lruzdu3aHYnR9uToP//5z459a9as4aCDDqp0DNL2940bNwYgPz+fnJwcQqEQOTk55OfnJ/x9M0VBQUGl3YDDhw/nq6++Yvz48cCek14JSFVNQ9s3YClwcCX77wHeA9bv6Rq1ZVMXWOplSrNrkFTH1ffJJ5/4YYcd5n/84x9969atlZ5TVlbmQFp0FVVl48aN/vvf/17dMwH49NNP/YADDvDVq1fvtL+srMznz5+fNr+b2lzH8SDBLrDXgIuBmyskTleb2X3AEdXIu0RE0sq8efM455xzGDFiBJdffnml52z/lzxAXl5e2v5LvmHDhtx7773ce++9lR53d9asWbNLF1v0+++//57c3NwqW5Fyc3MrvXamu/baa7n66qtp1arVTvuzsrLUtZhmYkmArqrqPHf/vZndmdyQRERq1uTJkxk0aBATJ07ktNNOq/K8ujJexsxo3rw5zZs359hjj630nPXr17Ny5cqdkqTCwsId7//73//StGnT3Y5FysnJIbImY0Z4/fXXeeONN5gwYULQoUgM9pgAuftmYPP2z2ZWDzjR3edEjq9IXXgiIqnj7tx2222MHTuWV199lSOPPHK359fEeJl0eXx93333pX379rRv377S41u3buWLL77YqdVo8eLFzJgxY8e+8vLy3SZIBx10EHvvvXeSv2Ewtm3bxuDBg7nttts0GWYtUZ2JEPcBZhFeoV1EpFbasmULAwcO5O233+a8887jqKPiWz955cqVZGdnx3RuPHPmVDWwtqLt3XErV64kNze3xrvj6tevT+vWrWndunWVrWDffffdTq1IK1as4MUXX9yRIH355Ze0bNlypwSpZcuWfPzxx+y777506NCBJk2a1Nh3SsSkSZOoX78+ffr0CToUiVF1Z4LOnDZNEalzvvvuO84991waNWrEnDlzyMrKYvTo0TGXLywsrHQixJpUG7rj9ttvP/bbbz+mTJlSZavWqlWrWLVqFfPmzdtp/9ixY2O+T2UJZk1OPLhu3Tquv/56Jk+enFFdfrVddRMgPfouIrXS8uXL6dmzJ927d+eee+6Jeb6cdFObHl+PtVWrqKiIn/3sZ4RCIbKyspg+fTodOnTYZbB2dIvSt99+y8SJE5k9e/aOrrUWLVowcuRI1q1bR9OmTfnwww9TmgSNHj2aLl260Llz55TdQ5JPa4GJSMZYuHAhZ511FsOGDeMPf/hD0OEkZPt8QtnZ2WnzNFpVrS7VGdcUCoXo0qVLTOcuW7aMZcuWVXps/fr1MXdVnnzyyTFPaLndf//7X+677z4WLVoUVzkJnhIgEckIzz33HP379+fhhx/mrLPOCjqcKlUnWYj1L3iIfw2vmh6UHQqFmDBhAr/97W+rndRVnHhw1qxZfPPNN1U+8r9+/Xpyc3PZe++96d+//y4Dt3Nzc6tcn+2GG25gwIAB/OAHP6hWrBIcJUAiUqe5O/feey933HEHL730UpWPfaeLWLuLakp1u69eeeWVao1JysrK4vDDD0+oRauy2bYPPfRQjj/++ErPD4VCrFy5cqek6LXXXtvxftCgQQwZMmSXcosWLeKVV15h8eLF1Y5VgqNB0CJSZ5WXlzN48GBmzZrF/Pnz9a/0aqhu91Ws42FStap8PBMPZmVl8X//93/83//9X8zXd3cGDx7MyJEj42qBk/RRrxpltgCTkh2IiEgyhUIhfvGLX7B48WLmzZun5KeaCgoKYl5OqKysjPnz51NWVhZzmXRq7apMKBSiqKiIUCi00/6pU6fy/fffc8kllwQUmSQq7gTI3Te7+0WpCEZEJBlWr17NSSedRKtWrXjhhRfYb7/9gg4pI2xvdUmHAdm7U1BQgJnFtGVnZ9O5c2eys7N32n/uuefy7rvvstdee+1SJt2TOgmrTguQiEjaevfdd+nYsSPnn38+48aNqzMzDUvyxNqqNX/+/B3JXFZWFvPnz8fduf322+nVq1etbdWSMCVAIlLj4vkXeLzb0UcfzQknnMDw4cM1KZ1UKtbfX+fOnXd0fW0f12RmDBs2jOeff77KckqAaoeEEiAzu8DM5pnZV2a2tZKtPFmBikjdEe+4EqDScSU333xzpdefMmVKUpIp/UVWNyUyrmngwIH84Q9/qNXjmiSs2o/Bm9mNwAjgS2A+UJqsoERE4H/zuQDk5eXtMuFfQUEBf/rTn7j22muZNm0aL7zwAj/60Y+qvE5Q62ZJ7RX9NNkHH3zAlClT+OijjwKOSpIhkXmAfgcUAqe5+5bkhCMi8j97Wu9qw4YNXHTRRXz11VfMnz+fAw44oFrXEYnF0KFDufHGG2natGnQoUgSJNIF1gR4SsmPiKTK9vWugF3Wu/ryyy/p3r07jRo1YsaMGVUmP3u6jkgsXnrpJT777DMGDhwYdCiSJIkkQP8BcpMViIhIRdtn9AV26rb68MMP6dSpEz//+c957LHHaNiwYbWuIxKLLVu2MHToUO688049VViHJJIA3QhcYWbHJCsYEckM8c7DAuw0D0teXh7Lli1j5MiR1KtXr9rX0SBoicW4ceNo3bo1PXv2DDoUSaJqjwFy99lm1g9YYGZFwGfA1l1P834JxCcidVB117uaMGECw4YN48knn6R79+7JD0ykgtLSUkaOHMmMGTM0rUIdk8hTYCcA/4hco2tkq8gBJUAikpDtj7w//vjjFBYWxrVmk0gi7r33Xs466yyOPPLIoEORJEvkKbB7Ca8LdhYw193XJickEZH/2bRpE5deeilLly5lwYIFtGjRIuiQJINce+21bNq0KegwJAUSGQN0JHCnuz+v5EdEUuGbb76hR48ebNq0iVmzZin5kRq3zz77sP/++wcdhqRAIgnQV8DmZAUiIhJtyZIldO7cmY4dO/LUU0+xzz77BB2SiNQhiSRA44G+ZpZIN5qIyC7mz59P165dGTx4MLfffjv16mnZQhFJrkSSl9eBMwk/BfYAsIxdnwLD3eckcA8RyTCTJ09m0KBBTJw4kdNOOy3ocESkjkokAXo16v3DhJ/4imaRffUTuIeIZAh3Z9SoUfz9739nxowZHHXUUUGHJCJ1WCIJ0CVJi0JEMtqWLVsYOHAgb7/9NkVFRbRu3TrokESkjktkIsQJyQxERDLTd999x7nnnkvDhg2ZM2eOlqkQkRqhkYUiEpjly5fTpUsX2rdvz7PPPqvkR0RqTMwJkJl1NrOql1uuulz9SNnsapQ9zcwWm9kSMxteyfELzey9yDbfzI6KtayIBGvhwoV07tyZfv368be//Y299tIDpSJSc+JpAZoL/Lwa99g/Uva4eAqZWX1gDHA6kAf0MbO8CqctA0529yOBW4BxcZQVkYA899xznH766YwZM4bBgwdrjSURqXHxJEAGHGBmbePZgNxI2XgdDyxx90/dfTPwJOFlN3Zw9/nuXhr5uABoE2tZEYlfPKu47277xS9+wZo1a/jlL3+p1ddFJBDxtjnfE9niVfER+Vi0BlZGfV4FnLCb8/sBL1WzrIjEIN5V3M0M9/D//uXl5QwePJiZM2fywgsv8MMf/jA1QYqIxCCeBGhEgvf6NM7zK2s1qjSRMrPuhBOgE6tRdgAwAKB58+YUFhbGGabEIxQKqY5TLF3qeMOGDQC89FL43yW33HILmzdv5q9//SufffYZn332WYDRJSZd6riuUv2mnuo4jgTI3RNNgOK1inD32XZtgNUVTzKzIwlPxHi6u38TT1kAdx9HZOxQ+/btvVu3bgkHLlUrLCxEdZxa6VDHoVCIvLzwsLv+/fvTrFkzjj32WMaOHcvee+8daGzJkA51XJepflNPdZzej8G/BbQzs4PNrAHwa2Ba9AmRMUZTgYvc/eN4yopI6hQXF1NaGh6e9/nnn9OpUycefvjhOpH8iEjdkLbPnbp7uZldBUwnvJzGeHf/wMyuiBwfC/wJOAB4IPIUSbm7H1tV2UC+iEgGys/PZ5999iEUCtG0aVPuvPNOPeklImklbRMgAHd/EXixwr6xUe/7A/1jLSsiNWPSpEk73n/22Wea4FBE0k46d4GJSC2zbds2hg0bxujRo5k3bx6Akh8RSUtKgEQkZnuaB6h+/frccccdfPLJJxx22GEAcc0PpHmARKSmKAESkZgVFBTg7rtsX375JSeccAJ9+vRhw4YNlZ4Ty6YESERqihIgEUnIhx9+SMeOHenRoweTJk2iUaNGQYckIrJHSoBEMkCylrCobMvLyyM/P59bbrlFT3qJSK2hBEgkA1TVdVXVBuzxnAkTJtC8eXNeffVVnnvuuYC/oYhIfGrsMXgze8Ddf1dT9xOR6gmFQjteK3uCy90ZMWIEEyZMoLCwcMeMzyIitUlNzgOktnGRNBe9hEVeXh4lJSU7JUGbNm2if//+LF68mAULFtCyZcugQhURSUiNdYG5+8CaupeIVE/0EhalpaUUFxfvOPbtt9/ys5/9jHXr1lFYWKjkR0RqtZQkQGbWxMyapOLaIpI6+fn55OTkAJCTk0N+fj4AS5cupVOnThx33HE8/fTT7LvvvkGGKSKSsKQnQGZ2K7AWKDWzFWb2nJndbGa9k30vEUmurKwsSkpKAHZ0f82fP58TTzyRq6++mjvvvJP69esHHKWISOJSMQboSqAL8AFwNPBj4BjgbLQiu0ja2z7mJysri6eeeoorr7ySCRMmcMYZZwQcmYhI8qSiC6wUeNPdv3f3Oe5+r7v/1t2PSsG9RCQG8c4DBOElLH71q1+xZs0aevbsqSUsRKROSUUCNA64IAXXFZFqimceoM2bN3PZZZdx1FFHsXLlSi1hISJ1Uiq6wH4FtDOz44BngbfdfW0K7iMiSfbdd99x3nnnsddeezF37lyys7ODDklEJCVS0QJUANwJ/BCYAHxjZkvN7OkU3EtEkmTFihWceOKJtGvXjmnTpin5EZE6LektQO7+LOGWHwDMrBnhQdBHJ/teIpIcixYtonfv3gwdOpTBgwdrTS8RqfOSngCZWS93f377Z3dfA7wS2UQkzUybNo1+/frx4IMPcvbZZwcdjohIjUhFF9idZvaamR2ZgmuLSBLdd999XHHFFbzwwgtKfkQko6QiAToc+Dfwmpk9ZGYtUnAPkYwR7yPs3bt3j/ncP/zhD5x//vkcf/zxQX9NEZEalfQEyN3L3f1uoAOwEfjQzK43s4bJvpdIJojnEXZ3B9jt8bKyMnr16sVPf/pTSktLueeeewL+hiIiNS9li6G6+zfuPgjoAfQHPjazPqm6n4iEV3OPfq3o888/5+STT+aAAw7gpZdeYv/996/J8ERE0kYq1gL7iZldYWYPm9m7wAJgb+At4Bozm2Nmucm+r0imC4VC5OXlAZCXl7dLElRcXEzHjh05++yzGT9+PA0aNAgiTBGRtJCKFqC5wIWEF0S9BTjE3XPd/Vx3/wnh8UFTUnBfkYxWXFxMaWkpAKWlpRQXF+849sorr3DKKacwatQobrjhBj3mLiIZLxUzQTdx9/KqDrr77Wb2pxTcVySj5efnk5OTQygUIicnh/z8fAAeeughbrrpJqZOncqJJ54YcJQiIukhKQmQmX0HvAu8A7xrZu8A77v7piqKdEvGfUXkf7KysigpKSE7O5uSkhL23Xdfhg8fzr/+9S/mzp1Lu3btgg5RRCRtJKsFqCtwRGT7BTAGqGdmHxNJioB33H0GgLsvTNJ9RSRKVlYWAPXr1+fXv/41q1evpqioiGbNmgUcmYhIeknKGCB3f8/dn3D364EPgEeALsANhJOsUcCDybiXSKaJdx4ggH333Zenn36aefPm0bx5892er5XcRSQTpWIQ9ABgkLu/5e7PufuvgcHAyym4l0idF888QB999BEHHXQQ119/PVu3bo2pjBIgEclEqUiAviI8G3S0B4HTU3AvEYmYPXs2J510EhdeeCG33nor9eqlbJovEZFaLxVPgd0JPGlmF7v7W5F9RwOacU0kRR5//HGGDBnCE088wV57peJ/axGRuiXpf1K6+8Nm1hh41czWAF8TToA0375Ikrk7I0eO5NFHH2XWrFkcfvjhFBYWBh2WiEjaS8k/Fd39XjN7lPDj7s2BEncvSsW9RDLV5s2bueyyy/jwww9ZsGABBx54YNAhiYjUGgklQGb2O6Af0Ab4FpgNTHb3We7+PTAt8RBFpKLS0lLOPvts9t9/fwoLC9l3332DDklEpFap9ijJyGzO9wP5wOrItS4l3PU13cxaJidEkdop3sfX49maNm3Kxgpq5UUAACAASURBVI0bmTJlipIfEZFqSOQxkcuBj4C27v5jd28PtAT+ABwPFJmZZl+TjBXP4+vuDrDHc4qKijjwwAMZM2YMRUVF1K9fP+BvKSJSOyWSALUEHnL3L7bvcPdSd78f6Ag0A0YkGJ9IRti+cnvFFdyjPf300/Tq1YuHH36Y3/3udzUVmohInZRIAlQKNKrsgLsvBsYBv0zg+iIZIRQKkZeXB0BeXt4uSZC7c/vttzNkyBBmzJhBz549gwhTRKROSSQBegG4yMwqTYKApcB+CVxfJCMUFxdTWloKhAc3FxcX7zi2ZcsWrrjiCp544gmKioo4+uijgwpTRKROSSQBugbIAV40sx9Wcvx0wougishu5Ofnk5OTA0BOTg75+fkAfP/995x55pmsXLmSuXPn0qZNmyDDFBGpU6qdALn7GqAH8APgEzObY2Z3mdlIM5sHHAMMSlKcInVWVlYWJSUlAJSUlJCVlcXKlSs58cQTOfTQQ5k2bRrZ2dkBRykiUrcktFiQu78P5AHXAw0JPwF2I9AJqA/cYmajzOxCM8s3s73jub6ZnWZmi81siZkNr+R4BzMrMrNNZnZNhWOfmVmxmb1jZgur+x1FakJWVtaO10WLFtGpUycuvvhixowZo6UtRERSIOHVEt19k7vf4e4nAAcApxDuHnuNcOvQUOAx4B2gLNbrmll9YAzhrrQ8oI+Z5VU47Vvg94TXH6tMd3c/2t2PjeMriSRFvPMAAZgZxx57LP/9738ZOnQo9erVq/J8reIuIlJ9SV0u2t3XRmaBvsvd+7r74UAToDPh7rCJcVzueGCJu3/q7puBJ4GzKtzvq8iCq1uS9BVEkibeeYDuu+8+WrVqxRtvvBHT+UqARESqL+Vt6+6+AVgQ2eLRGlgZ9XkVcEI8twZeMTMHHnT3cZWdZGYDgAEAzZs310KSKRYKhVTHFWzdupUHHniARYsWcdddd7F+/fqE6kh1nHqq49RS/aae6rgGEqAEWCX7PI7yXdx9tZm1AGaY2UfuPmeXC4YTo3EA7du3927dulUrWIlNYWEhquP/WbduHRdccAGhUIh33nmH/fffP+Frqo5TT3WcWqrf1FMdJ7kLLMlWAblRn9sQXnMsJu6+OvL6FfAM4S41kbTx+eefc9JJJ9G0aVNeeumlpCQ/IiISm3ROgN4C2pnZwWbWAPg1Ma4ub2aNzSx7+3vgZ8D7KYtUJE7FxcV07NiRs88+m/Hjx9OgQYOgQxIRyShp2wXm7uVmdhUwnfAj9ePd/QMzuyJyfKyZHQgsJDzQepuZXU34ibFmwDORJ2v2Ap5w95eD+B4iFb3yyiv07duXe++9lz59+gQdjohIRkrbBAjA3V8EXqywb2zU+y8Id41V9D1wVGqjE4nfQw89xE033cS//vUvunbtGnQ4IiIZK527wEQCEe/8PfFsAwYM4LzzzlPyIyISsGonQGb2sZldG+mGEqkz4p2/B9jt8fXr13P++efTpUsXvv76a/72t78F/A1FRCSRFqAtwG3ACjN71szONDO1KElGCYVCO71W9PXXX3PKKadQr149Xn31VZo1a1aT4YmISBUSWQz1cMIzPE8AugPPASvN7FYzOzRJ8YmkrVAoRF5eeHWWvLy8XZKgxYsX06lTJ7p3786kSZNo1KhREGGKiEglEl0MdYG7Xwa0AvoDy4DrgI/NbKaZXWBmDZMQp0jaKS4uprS0FIDS0lKKi4t3HJs9ezYnnXQS119/Pbfeeiv16qlxVEQknSTlT2V3X+/uj7r7iUAHwut2dSO8COpqM7vbzNom414i6SI/P5+cnBwAcnJyyM/PB+Dxxx/nvPPOY9KkSVx66aVBhigiIlVI2mPwkdXbewP9gNMIL1sxC9hEeCHUAWZ2gbs/l6x7igQpKyuLkpISsrOzKSkpoXHjxowYMYJHH32UWbNmcfjhhwcdooiIVCHhBMjMOhBOei4CWgBfAXcCD7n70sg5PwKeAm4nPFZIpE7IysoCoEGDBlx88cV8+OGHLFiwgAMP1MORIiLpLJHH4C81s3nAB8AQ4F3gPKCNuw/fnvwAuPsS4D5Ag6Ml7cU7DxBAw4YNmThxIm+99RatWrXa7fkFBQXBfkEREUloDNDDwMHAKOBQd/+5u//L3curOL+E8JggkbQWzzxAS5cupUOHDgwZMoTy8vKYyigBEhEJXiJdYOcA09x9aywnu/ubwJsJ3E8krSxYsIBf/vKX3HjjjVx55ZVBhyMiInFIpAWoF3BsVQfN7HgzG5/A9UXS1pQpU+jVqxcPP/ywkh8RkVookQToYnY/pudg4LcJXF8k7bg7d9xxB4MHD+aVV16hZ8+eQYckIiLVkMrV4BsTXi5DpE4oLy/nqquuoqioiKKiItq0aRN0SCIiUk1xJUCRyQx/GLWrg5mdVMmpTYGBwJLqhyYSn4KCAkaMGJGy6x966KG0a9eO119/nezs7JTdR0REUi/eLrBLgELCExw6cEPkfcVtKvBj4K/JClRkT+J5equsrAyAsrKyPZ67YsUK8vPz6dGjB88//7ySHxGROiDeLrBngc8AA8YD44CiCuc4EALecveViQYokmwVFzEtKSnZMaFhRW+//Ta9e/dm8ODBDBkyZMe8PyIiUrvFlQC5+7uEJzzEzH4A/Mvd309FYCKpUtkipp06ddrlvH//+99ccskljB07lnPOOaemwxQRkRSq9iBod0/dYAuRFNq+iGkoFNppEdNof/vb37jtttv497//zQknnBBAlCIikkoxJ0Bm9pvI28fc3aM+75a7T6xWZCIpUnER0+jur61btzJ06FBeeeUV5s2bx8EHHxxgpCIikirxtAD9g/D4nieBzVGfdzcowgElQJJ2tic90cnPunXruPDCC/n++++ZN28eOTk5QYUnIiIpFk8C1B3A3TdHfxZJF9V5DL6qQc1NmzbdZd/NN9+sdbxEROqImBMgd5+9u88iQSsoKKhWgvL+++9z5pln0q9fP2688UY96SUikgGSOhO0me0FnEV4IsTn3f2LZF5fJNlmzJjBhRdeyN13382FF14YdDgiIlJDqr0WmJndbmZvRX024FXgKeBBoNjMdrdWmEigHnnkES666CKmTJmi5EdEJMMkshjqacDcqM+9gJOAO4ALIvuGJ3B9kZTYtm0b119/Pbfddhtz5szhpJMqW81FRETqskS6wHKBT6I+9wKWuftwADM7HNA/qyWtbNy4kYsvvpiVK1eyYMECmjVrFnRIIiISgERagBoAW6M+dyfcBbbdp0CrBK4vklRff/01p5xyCgCvvfaakh8RkQyWSAK0EugIO1p7DgGinwxrQXhNMJHALV68mE6dOtGtWzeeeOIJGjVqFHRIIiISoES6wJ4EbjKzFsDhwPfAi1HHfwwsTeD6IkkxZ84czjvvPP7yl7/Qr1+/oMMREZE0kEgL0G2EZ4PuRHjG59+4+1oAM9sP6A28lmiAIomYNGkS5557LpMmTVLyIyIiOySyGOomoF9kq6iM8Pif9dW9vkgi3J0///nPPPLII8ycOZMjjjgi6JBERCSNJHUixO3cfRvwXSquLbInmzdvZsCAAbz//vsUFRXRqpXG4ouIyM4S6QITSZqCggLMLClbw4YNmTBhAosWLeKggw7CzLSGl4iI7CShBMjMOpnZJDN708yWmtmnFTYNgpaYFBQU4O4xbWVlZQCUlZXttP/TTz+lQ4cOXH311ZSXl+90TAmQiIhES2QpjN8ArwPnAI2AFcDyCtuKJMQoskMoFCIvLw+AvLw8QqHwTAtvvPEGXbp04corr+Tuu++mfv36QYYpIiJpLpExQDcAi4FT3X11kuIR2a3i4mJKS0sBKC0tpbi4mNWrV3PFFVfw6KOPcuaZZwYcoYiI1AaJJEA/AP6o5EdqUn5+Pjk5OYRCIXJycpg5cyZ///vfmT59Osccc0zQ4YmISC2RSAK0CmiYrEBEYpGVlUVJSQnZ2dn06NGDyZMnU1RURG5ubtChiYhILZLIIOixwIVmpsEWUqO2bdsGwOrVq3n99deV/IiISNwSSYAWARuBN83sUjPrbmYnVdySFKfUcfE8Br/ffvsB8PLLL7PffvvFVEZPgYmISLREEqDXgBMIr/n1MOGV4GdFbYWR12ozs9PMbLGZLTGz4ZUc72BmRWa2ycyuiaespJdYHoN/++23ad26Nbfffjvbtm2L+bF5PQYvIiIVJTIG6JKkRVGJSNfaGKAH4fFGb5nZNHcviTrtW+D3wC+qUVZqkRdeeIGLL76Yv//975x77rlBhyMiIrVcImuBTUhmIJU4Hlji7p8CmNmTwFnAjiTG3b8CvjKznvGWldpjzJgx/PnPf+b555+nY8eOQYcjIiJ1QFLWAjOzhkAz4Gt335yMawKtgZVRn1cR7nJLalkzGwAMAGjevDmFhYVxByqxC4VCMdfx1q1bGTt2LG+++SajR49m48aN+u8Tg3jqWKpHdZxaqt/UUx0nmACZ2THAncCJQH3CXU4zzawF8E/gNnd/tbqXr2SfJ7usu48DxgG0b9/eu3XrFuMtpDoKCwuJpY7XrVtH3759Wbt2Le+88w45OTmpD66OiLWOpfpUx6ml+k091XFiS2EcDcwFDgUmRh+LdE3tA/w2gdhWAdHPN7cBYp10MZGyErAvvviCbt260aRJE6ZPn67kR0REki6Rp8BGEk4qDgeGs2ury2uEx+JU11tAOzM72MwaAL8GptVAWQnQBx98QMeOHenVqxf/+Mc/aNCgQdAhiYhIHZRIF1hXwl1cocgYoIpWAAdV9+LuXm5mVwHTCXevjXf3D8zsisjxsWZ2ILAQaAJsM7OrgTx3/76ystWNRWrGq6++ygUXXMBdd91F3759gw5HRETqsERagBoB3+3meJMErg2Au7/o7oe5+6Hufmtk31h3Hxt5/4W7t3H3Ju6+f+T991WVleSLZwJDM6N79+5VHuvRowdff/01F110kSYwFBGRlEokAVoK/GQ3x3+KHjuv82KZwHD7VlZWBkBZWdmOfVu3buW6667j0EMP5aOPPtIEhiIiUiMSSYCeAC4ys1Oj9jmAmQ0FTgMeS+D6UoeEQiHy8vIAyMvLIxQKsXHjRi644AJmz55NUVER7du3DzhKERHJFIkkQHcCCwiPs5lDOPm528z+C9wOzAAeSDhCqROKi4spLS0FoLS0lLlz53Lqqafi7rz22ms0b9484AhFRCSTVDsBikx42AO4BthAeGHUw4A1wDDgTHfflowgpfbLz8/f8Th7dnY2V111FV27duWf//wnjRo1Cjg6ERHJNAlNhOju5cDdkU2kSllZWZSUlJCdnU15eTnDhw/nsssuCzosERHJUIlMhNjZzG41s6fNbHrk9RYz02JNGSSep8Cys7MB+PrrrxkwYEBMZTQIWkREUiHuBMjMmpjZC4Rngb4OOIdwV9g5wA3APDObZmbZSY1U0lIsT4Ft27aNW265hbZt2/LII4/E/NSYngITEZFUqU4X2BTgVOB14BHgPeB7wvP+HAn0B84EJgNnJCdMqa02b97M5ZdfTnFxMQsWLGDx4sVBhyQiIhJfC5CZ/Zxw8jPa3U9y9wnu/h93Xxp5neDuXYG7gJ+bWY9UBC21w9q1azn99NP59ttvmT17Nq1atQo6JBERESD+LrA+wHLCT3ntzjDCS2FcUJ2gpPZbtmwZnTt3Jj8/n6lTp9K4ceOgQxIREdkh3gToJ8Cz7u67Oyny+PuzwLHVDUxqrzfeeIMuXbowcOBA7rnnHurXrx90SCIiIjuJNwFqDcQ6iGMx0CbO60stN3XqVM4880wefPBBBg0aFHQ4IiIilYp3EHQToCzGc8uArDivL7WUu3PXXXdx9913M336dI455pigQxIREalSvAlQPSLrfcVxvtRx5eXlDBo0iHnz5lFUVERubm7QIYmIiOxWdR6DP8PMDozhvN2tFC91RFlZGb/61a/Ytm0br7/+Ok2aNAk6JBERkT2qTgJ0AbE/3RVPa5HUMqtWraJnz5507NiR+++/n7333jvokERERGISbwLUPSVRSK3zn//8h969e/P73/+ea665BjMLOiQREZGYxTVGx91nx7ulKnBJTDxreFW2HXPMMaxatYphw4ZRr149reElIiK1igYpZ6hY1vDavpWVhR/8Kysr4/777+fAAw+kqKhIa3iJiEitVZ0xQJJBQqEQeXl5ALRq1YqDDjqIefPmccghhwQcmYiISPWpBUh2q7i4mNLSUgA2bNjAfffdp+RHRERqPSVAslstWrRg8+bNQLgFqEuXLgFHJCIikjglQBkq1kHQP/rRj3YkQKtWrSI7OzumchoDJCIi6UwJUIba0yDoGTNm0Lx5cx577LGYB0trELSIiNQWSoBkF+PHj+fCCy9kypQp9O3bN+hwREREkk5PgckO27Zt46abbmLy5MnMmTOH9u3bBx2SiIhISigBEgA2btzIJZdcwvLlyykqKqJ58+ZBhyQiIpIy6gIT1qxZw6mnnsrWrVt57bXXlPyIiEidpwQow33yySd06tSJrl278uSTT7LPPvsEHZKIiEjKKQHKYHPnzqVr164MGzaM2267jXr19HMQEZHMoDFAGeqJJ57g6quvZtKkSfTo0SPocERERGqUEqAM4+7ceuutPPTQQ8ycOZMjjjgi6JBERERqnBKgDLJ582Yuv/xyiouLWbBgAa1atQo6JBERkUAoAcoQa9eu5ZxzziErK4vZs2fTuHHjoEMSEREJjEa9ZoBly5bRuXNn8vPzmTp1qpIfERHJeEqA6rg333yTLl26MHDgQO655x7q168fdEgiIiKBUxdYHfbMM88wYMAAxo8fT69evYIOR0REJG2oBaiWKSgowMxi2s4++2zWrFlD7969Yy6jVdxFRCQTKAGqZQoKCnD3KrctW7YwcOBA8vLydjzllZubS1lZ2W7Lbd+UAImISCZQAlSHlJWV0bt3b5YuXcq9995LWVkZAKWlpRQXFwccnYiISPpQAlRHrFq1iq5du9KmTRv+/e9/07FjR3JycgDIyckhPz8/4AhFRETSR1onQGZ2mpktNrMlZja8kuNmZvdFjr9nZsdEHfvMzIrN7B0zW1izkdesd955h06dOnHBBRfw4IMPsvfee5OVlUVJSQkAJSUlZGVlBRyliIhI+kjbBMjM6gNjgNOBPKCPmeVVOO10oF1kGwD8vcLx7u5+tLsfm+p4a0plg6B//OMfs2rVKq699lrq1au3Y392djYA2dnZGgQtIiISJW0TIOB4YIm7f+rum4EngbMqnHMWMNHDFgD7m1mdXt8hehD0mDFjOPDAA5k/f35MA5w1CFpERCQsnecBag2sjPq8CjghhnNaA58DDrxiZg486O7jKruJmQ0g3HpE8+bNKSwsTErwqbRt2zbGjh3LggULGD16NJs2baoVcQOEQqFaE2ttpTpOPdVxaql+U091nN4JkFWyz+M4p4u7rzazFsAMM/vI3efscnI4MRoH0L59e+/WrVsCIafe+vXr6du3L99++y3vvPMOTZs2DTqkuBQWFpLudVzbqY5TT3WcWqrf1FMdp3cX2CogN+pzG2B1rOe4+/bXr4BnCHep1Wpffvkl3bt3Jysri+nTp9e65EdERCRdpHMC9BbQzswONrMGwK+BaRXOmQb8JvI0WEfgO3f/3Mwam1k2gJk1Bn4GvF+TwSdbSUkJHTt25IwzzmDChAk0bNgw6JBERERqrbTtAnP3cjO7CpgO1AfGu/sHZnZF5PhY4EXgDGAJsB64JFK8JfCMmUH4Oz7h7i/X8FdImtdee40+ffowevRoLrrooqDDERERqfXSNgECcPcXCSc50fvGRr134MpKyn0KHJXyAGvAo48+yvDhw3nqqacyvr9WREQkWdI6Acpk7s5NN93EP//5T2bPnk2HDh2CDklERKTOUAKUhjZu3Mill17KsmXLWLBgAc2bNw86JBERkTolnQdBZ6wtW7bQtm1bZs6cqeRHREQkBdQClIays7MZNWpU0GGIiIjUWWoBEhERkYyjBKiGVLaIaTI3reElIiISOyVANSR6EdPdbWVlZeTmhie3zs3NpaysTIuYioiIJJkSoDRTXFxMaWkpAKWlpRQXFwcckYiISN2jBCjN5Ofnk5OTA0BOTg75+fkBRyQiIlL3KAGqIbGOAcrOzmblypUArFy5kuzsbI0BEhERSTI9Bl9DCgoKlKSIiIikCbUAiYiISMZRAiQiIiIZRwmQiIiIZBwlQCIiIpJxlACJiIhIxlECJCIiIhlHCZCIiIhkHCVAIiIiknGUAImIiEjGUQIkIiIiGUcJkIiIiGQcJUAiIiKScZQAiYiISMZRAiQiIiIZRwmQiIiIZBwlQCIiIpJxlACJiIhIxlECJCIiIhlHCZCIiIhkHCVAIiIiknGUAImIiEjGUQIkIiIiGUcJkIiIiGQcJUAiIiKScZQAiYiISMZRAiQiIiIZRwmQiIiIZBwlQCIiIpJxlACJiIhIxknrBMjMTjOzxWa2xMyGV3LczOy+yPH3zOyYWMuKiIhI5krbBMjM6gNjgNOBPKCPmeVVOO10oF1kGwD8PY6yIiIikqHSNgECjgeWuPun7r4ZeBI4q8I5ZwETPWwBsL+ZtYqxrIiIiGSodE6AWgMroz6viuyL5ZxYyoqIiEiG2ivoAHbDKtnnMZ4TS9nwBcwGEO4+A9hkZu/HHKFURzNgTdBB1HGq49RTHaeW6jf1MqWOf1DVgXROgFYBuVGf2wCrYzynQQxlAXD3ccA4ADNb6O7HJha27I7qOPVUx6mnOk4t1W/qqY7TuwvsLaCdmR1sZg2AXwPTKpwzDfhN5GmwjsB37v55jGVFREQkQ6VtC5C7l5vZVcB0oD4w3t0/MLMrIsfHAi8CZwBLgPXAJbsrG8DXEBERkTSUtgkQgLu/SDjJid43Nuq9A1fGWjYG4+KNUeKmOk491XHqqY5TS/WbehlfxxbOIUREREQyRzqPARIRERFJiYxIgLSkRuolWMefmVmxmb1jZgtrNvLaI4Y67mBmRWa2ycyuiaeshCVYx/odxyCGOr4w8mfEe2Y238yOirWshCVYx5nzO3b3Or0RHgS9FDiE8OPx7wJ5Fc45A3iJ8PxBHYE3Yi2rLbE6jhz7DGgW9PdI5y3GOm4BHAfcClwTT1ltidVx5Jh+x8mp485ATuT96frzuObqOPI5Y37HmdACpCU1Ui+ROpbY7LGO3f0rd38L2BJvWQESq2OJTSx1PN/dSyMfFxCexy2msgIkVscZJRMSIC2pkXqJ1DGEZ+l+xcwWRWbmll0l8lvU7zg2idaTfsd7Fm8d9yPcclydspkqkTqGDPodp/Vj8ElSI0tqZLhE6higi7uvNrMWwAwz+8jd5yQ1wtovkd+ifsexSbSe9Dves3iWKepO+C/nE+Mtm+ESqWPIoN9xJrQAJbKkRixlJbE6xt23v34FPEO4CVd2lshvUb/j2CRUT/odxySmOjazI4GHgbPc/Zt4ykpCdZxRv+NMSIC0pEbqVbuOzayxmWUDmFlj4GeAFqTdVSK/Rf2OY1PtetLvOGZ7rGMzawtMBS5y94/jKStAAnWcab/jOt8F5lpSI+USqWOgJfCMmUH49/iEu79cw18h7cVSx2Z2ILAQaAJsM7OrCT/98b1+x3uWSB0TXllbv+M9iPHPij8BBwAPROqz3N2P1Z/HsUmkjsmwP481E7SIiIhknEzoAhMRERHZiRIgERERyThKgERERCTjKAESERGRjKMESERERDKOEiARERHJOEqARERkBzNraGYPmdmnZhYys08i8x2J1Cl1fiJEERGJy17AF4RnAf4UOBKYbmafu/vkQCMTSSK1AIlI0plZNzNzM7s46FhqCzMriNTZ9u3EPZdKPndf5+43ufsSd9/m7u8ALwBdKsTbsUK8BUHEK1JdSoBE0lxUMrF922pmpWb2vplNMLPTLDJ3fTozs6Mjf8n/MOhY0txg4CJgcdCBAJjZXoRXC3+vwqElhOMcXONBiSSBusBEao9/El5TzYBsoD3wC+A3wKtmdp67rw0wvmhzgH2ALVH7jgZuBgqBz2o+pFrjWXf/LBUXNrOGwN67OWWDu2+tsO8+4DtgYvROd18DPB5JaO9OYpgiNUIJkEjt8ba7Px69w8yGALcDQwgnSKcHEVhF7r4N2Bh0HLKLR4ALd3O8O+EEFQAzG0249een7r45taGJ1Cx1gYnUYu6+1d2HAq8Dp0WPG4k8zXO9mX1gZhvNbK2ZPW9mP46+hpldHOla+6mZXWNmS81sk5l9bGa/rXhPM2sU6cpabGbrI9ctNrM7os7ZaQxQZHzIo5HDs6K68/5hZmdH3vev7DtG4l8STzefmTWPXHtNpMtwtZndXBu6CisT9d/oFDP7k5ktN7MNZvaGmXWMnHOymb1uZuvM7HMzu6niddy9r7vbbrbCqHveQ3gg9CmR1h6ROkUtQCJ1wyOE/6XeE3jdzPYGXgY6A48B9wP7AZcB88zsJHdfWOEafyHcbfUgsAkYCPzDzJa4+7yo88YAlxLuErkbqA+0A366m/imAq2AAZH7fBjZvxR4i/BTR/2Ah6MLRf5yzwNucHePpSLM7AjgVWBtJL6Nke9dAHxDuC5qq1GE6/teoAEwlPATWr8l/BsYB0wCzgdGmtmyiq2GsTCz+wj/9+zu7l8nK3iRdKIESKRu2D5A9bDI61VAN+A0d5++/SQzewB4H7gzcjxaQ+C47V0dZjaF8GPQVwHRCdAvgZfcfZfWoaq4+3tmVkQ4AZoR3dIQudejwHVmlufuJVGH+gFbgX/Ech8zawQ8R3gA8enuvj6y/wlgReR6tTkBqg90jPpvVEL4+04BOrn7W5H9jwDLgSuBuBIgM/sBMIhwErwsqtFsrrunRRerSDKoC0ykbvg+8tok8toX+Aj+v737C5GyCuM4/n2wIl0FXRUlKjOKkmrDRbC/RKkZ7cZvFwAABQdJREFUGglCLBsZYRfVRX9uWoIuSgwqKsmbjDai7I9IVLiFkaEhgtWGa2h54eZu2p8NK23LCit6unjOC6/jzOzM7uLs7vw+MMy88573vGcZ2Hk453nOsMvMpmUPYtbgI+AaMxtf0Mfz+TwPd/8e2E/M7uT1A5ekmZbh0g44EaAAYGYNQAsRbP1QYT/3AecBK7PgB8Dd+4CvgLOHa8A1sq4gF2dHev40C34AUptOTv7sBuTuB9Ny2JnuPjH3UPAjY4oCIJGxIQt8skBoDnAx8FORx0piJmFaQR89Rfr9BZha8N6DwBRgb8oXesnMlpnZoP+fuHsvsWy1Ii3fQSzjTKJgWWwAtwKfuPuBIucmAEcq6cTMOs2spdKbVtt+CE74jNz9aHrZW6TtUU7+7EQk0RKYyNjQlJ6zvWMM2EtUh5VSmNtRWP6cOSFx2N03pdLnJcB1wEJi5maHmS0cQrXQi8BbwC3A26nPH4lN+AZkZlOAeUTZduG5RmA2FSylmdk44DLgiwrvW1X7ISr1GZV6X0RKUAAkMjZkS0dZsNANTAe2pZL0YeXuR4jcktdTZdWTQBuwjAhiil42QLebgMPAXWb2JbHz8FPu/m+Fw5pLBGvHipy7jfh/txHAzGYCzxFl3w1EftByYnnvAJEPtSvlv9xIBDcvA9cSyeT7iSTxQ8Xau/vOCscsIjWiJTCRUczMxpnZM0QF2OZctdZ6YCYlZoDMbMYQ7jc5/16qztqdDhvLXH6sXBt3/4eYoVlMbJgIUdlUqay8f0G+3N3MLgIeB7YC23L9dgPnApOBO4HvUq7QSqAzl/uykwh6NhKbTzYSuTdry7QXkRFOM0Aio0ezmd2eXud3gp4FbCFmOTJrgUXA02Z2A/HF/xvxhb+AKA2/fhBjmAT0mVkHEfQcJpaW7iVyTt4rc+3nwH/AI2m56g+g190/y7VpBx4CWoHt7t5dxdia03hmAR1m9gFwPnA3kfuzIjcbdiFROXeGu/9OJAzn+9mdO86SqN/NjlOF3E2l2ovIyKcZIJHRo5XY0+dVYslpEbCdKPde7O79WcM0m7IUeIBYCltF7InTQiTSPjHIMfxJLB3NJgKVdcTvQXUA88tVa7n7IWK2ZHy6bgMROOXbfA18nA6rmf2BmAHaTSzDTQeeJarh3iDK+/tybe8ALgW+MbP3zawpd66ZgnweM7vZzLamDRX7gQ+BfaXai8jIZxXuLSYickqY2WbgSuAsd/+rwmsmEDNca9y9rYp7ZRs/TnX3pem9n4El7t6Zjq8mkrJbgZ3uftzMNgD73H11YfvBSrtlP0oEVN8Cv1aR/3TKWfxI6mTgHKALWOXuj9V0UCJV0AyQiIwYZnYBkQP0WqXBT3I5Udq/t4J7LDezOalsvwGYAexK504nSvzz5hIVc3uA08zsYaJEv6tE+6HqSve7Ypj7HW7ziHF21XogIoOhHCARqTkzm0/sXXQ/8DewpsousgToPWVbhauIHKlGYp+jN4HVEEuHKal8S5rhaErnlxOzMgfT2Iz4cdqT2rt7j5m9kPq7p4q/YT3xm26ZAYO5GttHLMNmiu0jJTJiaQlMRGrOzF4h8nJ6gDZ3f6fK69uJSq6J7n582AcoImOOAiARERGpO8oBEhERkbqjAEhERETqjgIgERERqTsKgERERKTuKAASERGRuqMASEREROqOAiARERGpOwqAREREpO4oABIREZG68z9fmw3MUzlRqgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 648x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "\n",
    "fig = plt.figure(figsize = (9,6))\n",
    "ax = plt.axes()\n",
    "ax.errorbar(xdata, N_avg, N_std, capsize = 12, elinewidth = 1, c = 'k', marker = '.', markersize = 5, linewidth = 0)\n",
    "ax.plot(xdata[1:15], func(xdata[1:15], a[0], b[0]), c='k', linewidth = 1)\n",
    "ax.plot(xdata[14:22], func(xdata[14:22], a[1], b[1]), c='k', linewidth = 1)\n",
    "ax.plot(xdata[21:25], func(xdata[21:25], a[2], b[2]), c='k', linewidth = 1)\n",
    "ax.grid(True)\n",
    "ax.set_xlabel(r\"Density $\\rho_{stat.}$ [m$^{-2}$]\", size = 18)\n",
    "ax.set_ylabel(r\"Density $\\rho_{dyn.}$ [m$^{-2}$]\", size = 18)\n",
    "ax.set_xlim(0)\n",
    "ax.set_ylim(0)\n",
    "\n",
    "fig.savefig(\"phase_transition.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ec3c0dae-8786-4388-8dd4-b68d1d9f63c3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1.54 ± 0.035) rho + (0.011 ± 0.0031)\n",
      "(-0.212 ± 0.001) rho + (0.273 ± 0.0002)\n",
      "(3.156 ± 0.317) rho + (-0.481 ± 0.0758)\n"
     ]
    }
   ],
   "source": [
    "for i in range(3):\n",
    "    print(f\"({round(a[i],3)} ± {round(ea[i],3)}) rho + ({round(b[i],3)} ± {round(eb[i],4)})\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "293dc563-31b2-428e-a1fa-ff64b693fad7",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}