{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "de4c7922-8621-40b8-8acf-ec4e51c6e377",
   "metadata": {},
   "source": [
    "# Linear Regression\n",
    "\n",
    "In this notebook, we will explore some more operations, such as calculating the logarithm in numpy and doing a linear regression.\n",
    "\n",
    "You have the following isothermic reaction:\n",
    "$A+\\frac{1}{6}B→\\frac{1}{4}C+\\frac{1}{2}D$\n",
    "\n",
    "During a laboratory experiment you measure in a batch reactor with constant volume and the initial concentration of $C_A$ is 25 $mol \\cdot m^{-3}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "13048039-3f7d-49ad-bbde-06353d36dfb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# define a distribution of CC\n",
    "cc = np.linspace(0, 4, 11)\n",
    "time = [0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20]\n",
    "\n",
    "# create a pandas dataframe with the data we simulated \n",
    "df = pd.DataFrame()\n",
    "df['Time'] = time\n",
    "df['Cc'] = cc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c9883c57-fe68-4693-8250-dfbc0a86f8a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a function to calculate CA given CC and CA0\n",
    "def calculate_CA(C_A0, Cc):\n",
    "    C_A = C_A0-(4*Cc)\n",
    "    return C_A\n",
    "\n",
    "# Looking at the stochiometry, we know that CA can be calculated from CC and CA0\n",
    "df['Ca'] = calculate_CA(25, df['Cc'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24a6a7bd-e3a8-472d-bb4e-a7d07e4002a6",
   "metadata": {},
   "source": [
    "## Finding the logarithm of a value in Python\n",
    "\n",
    "Here we will use the numpy library to calculate the logarithm of the values defined."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0ce6a9cf-25ca-4670-b920-b99660123d06",
   "metadata": {},
   "outputs": [],
   "source": [
    "# find the log of CA\n",
    "df['ln(Ca)'] = round(np.log(df['Ca']), 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d886559-d7cd-4521-a6d9-f86446051a98",
   "metadata": {},
   "source": [
    "## Simple linear regression in Python with SciPy library\n",
    "\n",
    "Here we calculate a linear least-squares regression for two sets of measurements.\n",
    "Check the documentation [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.linregress.html).\n",
    "\n",
    "The unction returns:\n",
    "- Slope of the regression line\n",
    "- Intercept of the regression line\n",
    "- The Pearson correlation coefficient. The square of rvalue is equal to the coefficient of determination\n",
    "- The p-value for a hypothesis test whose null hypothesis is that the slope is zero, using Wald Test with t-distribution of the test statistic. See alternative above for alternative hypotheses\n",
    "- Standard error of the estimated slope (gradient), under the assumption of residual normality\n",
    "- Standard error of the estimated intercept, under the assumption of residual normality"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ceb896e9-f581-4b1c-8507-a1ffb05d2cab",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Time</th>\n",
       "      <th>Cc</th>\n",
       "      <th>Ca</th>\n",
       "      <th>ln(Ca)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>25.0</td>\n",
       "      <td>3.22</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>0.4</td>\n",
       "      <td>23.4</td>\n",
       "      <td>3.15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4</td>\n",
       "      <td>0.8</td>\n",
       "      <td>21.8</td>\n",
       "      <td>3.08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>6</td>\n",
       "      <td>1.2</td>\n",
       "      <td>20.2</td>\n",
       "      <td>3.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>8</td>\n",
       "      <td>1.6</td>\n",
       "      <td>18.6</td>\n",
       "      <td>2.92</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>10</td>\n",
       "      <td>2.0</td>\n",
       "      <td>17.0</td>\n",
       "      <td>2.83</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>12</td>\n",
       "      <td>2.4</td>\n",
       "      <td>15.4</td>\n",
       "      <td>2.73</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>14</td>\n",
       "      <td>2.8</td>\n",
       "      <td>13.8</td>\n",
       "      <td>2.62</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>16</td>\n",
       "      <td>3.2</td>\n",
       "      <td>12.2</td>\n",
       "      <td>2.50</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>18</td>\n",
       "      <td>3.6</td>\n",
       "      <td>10.6</td>\n",
       "      <td>2.36</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>20</td>\n",
       "      <td>4.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>2.20</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    Time   Cc    Ca  ln(Ca)\n",
       "0      0  0.0  25.0    3.22\n",
       "1      2  0.4  23.4    3.15\n",
       "2      4  0.8  21.8    3.08\n",
       "3      6  1.2  20.2    3.01\n",
       "4      8  1.6  18.6    2.92\n",
       "5     10  2.0  17.0    2.83\n",
       "6     12  2.4  15.4    2.73\n",
       "7     14  2.8  13.8    2.62\n",
       "8     16  3.2  12.2    2.50\n",
       "9     18  3.6  10.6    2.36\n",
       "10    20  4.0   9.0    2.20"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "fce433ff-7735-4752-bce7-8f4746d90c0a",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.\n",
      "C:\\Users\\fiacac\\AppData\\Local\\Temp/ipykernel_15964/2368782444.py:17: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure.\n",
      "  fig.show()\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAEWCAYAAABWn/G6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABD/0lEQVR4nO3deXwU9f3H8deHJEAg3KeEm0AU5AiEwwMFRUHbgrZYsBTxKtqKtdaitYflZy8BxaNaKyqeVdoqxbRFAgrBGwFBIkhIgCCJnIEAgYRcn98fM8Fl2c0Be8Ln+Xjkkd2Z78y8d3ZmPruzc4iqYowxxoRKvXAHMMYYc3axwmOMMSakrPAYY4wJKSs8xhhjQsoKjzHGmJCywmOMMSakIq7wiMhwEckKdw4AEZkhIq+GO4cJLxH5m4j8NgjjFRF5QUQOiMingR7/2URE/iAi+0RkV7izmJqFrfCISK6IjPLurqrvq2pyODKFg7/5YCKHqt6uqr8PwqgvBq4AOqrqkNMdmYh0FREVkdjTjxZ4IvKiiPwhCOPtDNwD9FbV9oEe/9lORDJE5NZAjjPivvGES6SurJFCRGKCNN6Azvcoex+7ALmqeqSuA0bZ6wy2zkCBqu4Jd5BACdb6FjFUNSx/QC4wykf3EUCeV7tfAOuBg8A/gIYe/b8NrAMKgY+Afh79fglsAQ4DG4FrPfrdCHwIPAoUAH/wkWUG8KrH82HuNAqBz4ERHv1uAr50p7UVuM2jX2vgv+5w+4H3cYr+K0AlUAwUAff6mVc3uuM8DGwDJrndY4CHgX1u/zsABWJ9zWMfr+dfwC53vr4H9PHo9yLwNLAIOAKMAjoAbwJ73Rw/9Wg/BFgNHAJ2A3P8vJYRQB5wnzvtV9x5UfVeFQD/BFp6DHMDsN3t91vP1+W+pjeAV91p3wo0A54HdgL5wB+AGLd9ErDCfc37gH+43cVdFva448kEzveYF3/wyPMjIMd9L9OADh79FLgdyHbf76cA8TEfbgFKgAr3vf+/Wo77Dnfc23yM8yu3TZH7d4E73wa5/Se5/ft4ZFjoPm4APAZ87f49BjTw8x76nIduv3OBpW7+LOD7bvepQBlQ6mb7j9v9Pvc9Ouy2v9zPNJsBL+Mse9uB3+AsN6Nw1p9Kd7wvVrPM3eu+vzuBa4Crgc1u1l95tK9peaxuvbkaZ1tz2H1dv/BYhz/wyqVA0imubzPcHK+608oEegH3u69xB3Cl1/zzt07cCHyAsy054E7rKrffH3GW0RJ3/j5JNetKrbf/p1M8TuePuhWeT903oSXOxv12t1+K++KH4myEp7jtG7j9r3OHqwdMcN/QczxmdjlwJxALxPvIMgN3Qw0kugvh1e74rnCft3H7fwvo4b4plwJHgYFuvz8DfwPi3L/huBsjf/PBI0Nj981Ndp+fwzcbjtuBTUAnd94sp26F52agCd9sdNZ59HsRZ8W6yH29jYA1wANAfaA7TrEb7bb/GJjsPk4Ahvl5PSPc+T7TnW48cBfwCdDR7fYM8LrbvjfOAn+xO92HcTZgnoWnDGdDUs8d37/dcTQG2uIsP7e57V8Hfu22bQhc7HYf7b6+5u57eB7fLCsv4hYe4DKcje1AN+tfgPe8Nib/dcfTGWejMcbPvLgRj41RLce91H2vfS2vXT3ff7fby8A97uO5OBvTH3v0u9t9/KD7HrQF2uB8wPq9n9z+5mFjnA3eTTjrVIr7enp7z0f3ebLbvoNH/h5+pvky8BbO8toVp2Dc4mubUc0y9wDO+vcj9315zR1fH5zi1c1t73d5rMV6sxMY7j5uwTfbgBPea4/307Pw1GV9m4FTDEa78/plnILxa4/XuM1jWtWtEzfirEM/wtmO/hjnw0fVNioDuNVjXH7XlVpv/+vSOJB/1K3w/NDj+Szgb+7jp/FaOXA+NV3qZ5rrgHEeM/urGjLO4JvCcx/wilf/dGCKn2EXAnd5rNRvVS1ktZkPHv0b43xy/h5eGxtgGW4Rdp9fSR0Kj9e4mrvDNvNYEV726D/Ue37hfLp6wX38HvB/QOsa5ukInE+9nt9av8Tjky5OcS3DWaEe4MSVvpE7vGfh8dw4twOOec4r4Hpgufv4ZZwNcEevXJfhbMyGAfW8+r3IN4XneWCWR78EN2tX97nibojd5/8EfulnXtzIiYWnNuO+rJp525WTC88tQJrHfL4VmO8+3843G8YtwNUew43G2Q3oazr+5uEE4H2vbs8Av/Oej+7zJJwPjqOAuGpeV4z7nvf26HYbkOGxTNVUeIr55hN+E3c+DfVoswa4pqblsRbrzVdutqbVvdce76dn4anL+jYDWOrR7zs4H9C8X2Nzal4nbgRyvNYxBdq7zzM4sfD4XVdq+xctv/F4HqlyFGeFBGcf+T0iUlj1h/PpvwOAiNwgIus8+p2Ps9uryo46ZOgCXOc1rYtxFkpE5CoR+URE9rv9rvaY1myc3SdLRGSriPzS30TcI6iK3L9fqbP/fwLOt5udIvI/ETnXbd7B6zVsr+2LEZEYEXlIRLaIyCGcIgX+508XoIPX6/8VzkINzgauF7BJRFaJyLermfxeVS3xGve/Pcb7Jc7X+3ber1FVj+J80/TknTMOZ15Vje8ZnE954OxuEeBTEdkgIje7412GsxvhKWCPiMwVkaY+snfAYz6rapGbJ9Gjjb/ltSa1GXddlllwdokNF5FzcDbg/wQuEpGuOLtf1vmatvu4g59x+pyHOPN+qNcyMgnw+YO/quYAP8PZiO4Rkfki4muarXHeU+98iT7a+lOgqhXu42L3/26P/sWcuF3xuTzWYr35Hs66v11EVojIBXXIWJf1zVf+fT5eYwI1rxPgscy661jVsCepw7riV7QUHn92AH9U1eYef41U9XUR6QI8C0wDWqlqc+ALnBWmitZxWq94Tauxqj4kIg1w9sU+DLRzp7WoalqqelhV71HV7sBY4OcicrmvDOocQZXg/v3J7ZauqlfgFLlN7usC52t9J4/BO3tlPoLz6aWK5wbgB8A4nE+bzXA+LYP/+bMD56u75+tvoqpXuxmzVfV6nIV5JvCGiDQ+eTae/JrdcV/lNe6GqprvvsaOVQ1FJB5oVc34duB8umvtMa6mqtrHzblLVX+kqh1wPpn+VUSS3H5PqOognN17vYDpPrJ/jbMiV+Vp7ObJ9/Na66I2465umT2pn7txP4qzS/k9VT2Es5GZivMJvNLXtHGWpa99TsT/PNwBrPB6HxNU9cfV5HtNVS92p604y463fTjfOLzzBWKe+1Ld8ljteqOqq1R1HM56sBCn0IPXuigivopxrde3U3g9fteJWvD1vtVmXfEr3IUnTkQaevzV9UidZ4HbRWSoe05EYxH5log0wdlFpTj7chGRm3C+8ZyqV4HviMho91NPQxEZISIdcfbBNnCnVS4iV+Hs9sKd9rdFJElEBGc/bgXOj6HgfGrp7m+iItJORMa5G6FjOF+nq4b9J/BTEekoIi1wfhD1tA6YKCJxIpIKjPfo18QdXwHOCvGnGl7/p8BhEblPROLdeXC+iAx2c/5QRNq4G7JCd5hKfyPz8jfgj+6HBUSkjYiMc/u9gTPfLxSR+jifjsX3aEBVdwJLgEdEpKmI1BORHiJyqTvu69z3DJwfUhWoFJHB7nIUh7ORKPGT/3XgJhEZ4H7g+BOwUlVza/laq3O6497rZvZenlbgfABb4T7P8HpeNe3fuPO+Nc4uTp/nsPmbhzi/bfUSkcnuMhfnztfz3LYnLOsikiwil7mvtYRvDhI4gfsp/p84y0gTdzn5ub98AVDd8uh3vRGR+iIySUSaqWoZzm+zVa/nc6CP+942xFmOq1Pt+lYXNa0TteD9vtV2XfEr3IVnEc7CVvU3oy4Dq+pqnB/EnsRZAXJw9leiqhuBR3B+9N4N9MU5iu2UqOoOnE86v8JZwXfgVPl6qnoY+CnOynEA51NRmsfgPYF3cIrGx8BfVXW52+/POCt8oYj8wsek6+GsZF/jHH1zKc6Pf+AU3nSchfozYIHXsL/FOeDhAM7vL6959HsZZ3dFPs5ROJ/U8PorcI4gHIDzI+Y+4DmcT30AY4ANIlIEPA5MVNViH6Py5XGc+bVERA67WYa6092A82l9Ps63nyKc3wWOVTO+G3A+DGzEee1v4O4SBQYDK92caTi/w20FmuLMzwN8cwTdbB/z4R2c+fqmm6cHMLGWr7NapztudxfJH4EP3eVpmNtrBc4G8z0/z8E5ymk1ztGjmTjLk79zbnzOQ3c9uNLN/DXON6uqg0jA+Q2rt5ttodv9IZxlaRfOt4T7/UzzTpyN3FacI7BeA+bVMEtOld/lkZrXm8lArrsb7nacXY2o6mac33rfwTkq8YPqAtRifaur6taJmjwOjBfnROcnqGZdEZFficjbNY2w6qgFcwYQZ7/9NpwfasvDHCcoRCQB5xtVT1XdFuY4xphTEO5vPMbUSES+IyKN3N2ND+N8Is8NbypjzKmywmOiwTi+ObGxJ85uPPuqbkyUsl1txhhjQsq+8RhjjAmpM+ZCg61bt9auXbsGbfxHjhyhcWN/p6VEjmjJCdGT1XIGVrTkhOjJejo516xZs09V2wQ4UvVqc3mDaPgbNGiQBtPy5cuDOv5AiZacqtGT1XIGVrTkVI2erKeTE1itId5e2642Y4wxIWWFxxhjTEhZ4THGGBNSZ8zBBcYYc7YrKysjLy+PkpKSk/o1bNiQjh07EhcXF4ZkJ7LCY4wxZ4i8vDyaNGlC165dca5J7FBVCgoKyMvLo1u3bmFM6AjqrjYRGSMiWSKSIz7uQSMil4jIZyJSLiLjvfotdi8m+N9gZjTGmGi1cG0+Fz20jMz8g1z00DL2HzpCq1atTig6ACJCq1atfH4TCoegFR4RicG5UdBVOPdsuF5Eens1+wrnatKvcbLZOFd6NcYY42Xh2nzuX5BJfqFzEfj8wmIKj5ZSWFzms713MQqnYH7jGYJzO9WtqlqKc1n7cZ4NVDVXVdfj+x4c7wKHg5jPGGOi1uz0LIrLKk7oVqmw+2BkfKupTjB/40nkxFu55vHNPS0CQkSm4txJkXbt2pGRkRHI0Z+gqKgoqOMPlGjJCdGT1XIGVrTkhMjOOrHT4eP3H24XD/f0LadpfaVF/UoOH/b9mb2kpCQiXk9UH1ygqnOBuQCpqak6YsSIoE0rIyODYI4/UKIlJ0RPVssZWNGSEyI7668fWnZ8N9s9fct5JDOWud1g/zGhXasmJ7VXVRo2bEhKSkqoo54kmLva8jlejwHoSPDukW6MMWeV6aOTiY+LOaFb3qFy4rUY9brrQNVRbQ0bNgxlRL+C+Y1nFdBTRLrhFJyJOLeENsYYc5quSUkEnN964DCJzeNp3fYctLSYTZs2ndS+6jyeSBC0wqOq5SIyDUgHYoB5qrpBRB7EuShdmogMBv4NtAC+IyL/p6p9AETkfeBcIEFE8oBbVDU9WHmNMSbaXJOSyDUpiWRkZHDnpBHhjlNrQf2NR1UXAYu8uj3g8XgVzi44X8MOD2Y2Y4wx4WHXaguA2267jcaNG7Ns2bITus+ZM4fevXvTr18/Lr/8crZv317juLZt28bQoUNJSkpiwoQJlJaWntSmtLSUm266ib59+9K/f/8TjlLJysqib9++JCUl8dOf/vSkfb3GGBNuVnhOgapSWemcevSHP/yBwsJCVq5cyR133MH69euPt0tJSWH16tWsX7+e8ePHc++999Y47vvuu4+7776bnJwcWrRowfPPP39Sm2effRaAzMxMli5dyj333HM8z2OPPcazzz5LdnY22dnZLF68OBAv2RhjAsYKTy3t2rWL5ORkbrjhBs4//3x27NjBSy+9xIYNG3jttdc4//zzSUtL40c/+hE7djinL40cOZJGjRoBMGzYMPLy8qqdhqqybNkyxo93rh40ZcoUFi5ceFK7jRs3ctlllwHQtm1bmjdvzurVq9m5cydHjhxh2LBhiAg33HCDz+GNMSacrPDUQXZ2Nj/5yU/YsGEDXbp0YcqUKbz++uvExDiHNPbs2ZOVK1fSqVOnk4Z9/vnnueqqq44/HzBgwEltCgoKaN68ObGxzk9vHTt2JD//5CPQ+/fvT1paGuXl5Wzbto01a9awY8cO8vPzadPmmzvY+hveGGPCKapPIA21Ll26MGzYsDoP9+qrr7J69WpWrFhxvNu6detOOcfNN9/Ml19+SWpqKl26dOHCCy88XvyMMSbSWeGpg8aNG9d5mHfeeYc//vGPrFixggYNGlTbtlWrVhQWFlJeXk5sbCx5eXkkJiae1C42NpZHH330+PMLL7yQXr160aJFC/bu3Xu8u7/hjTEmnGxXWxCtXbuW2267jbS0NNq2bVtjexFh5MiRvPHGGwC89NJLjBs37qR2R48e5ciRIwAsXbqU2NhYevfuzTnnnEPjxo355JNPUFVefvlln8MbY0w4WeEJounTp1NUVMR1113HgAEDGDt27PF+vn7jAZg5cyZz5swhKSmJgoICbrnlFgDS0tJ44AHnFKg9e/YwcOBAzjvvPGbOnMkrr7xyfPif/exn3HrrrSQlJdGjR48TflcyxphIYLvaaql9+/Z88cUXdRrmnXfe8dvP32883bt359NPPz2p+9ixY48Xrq5du5KVleVz+OTk5DrnNMaYULJvPMYYY0LKCo8xxpiQssJjjDEmpKzwGGOMCSkrPMYYY0LKCk8dPPHEE5x33nlMmjQpZNN88cUXmTZt2ikPv3jxYpKTk0lKSuKhhx7y2ebYsWNMmDCBpKQkhg4dSm5u7vF+f/7zn0lKSiI5OZn09ODdDqm6DJ78vZ5bbrmF/v37069fP8aPH09RUVHQshpjTo8Vnjr461//ytKlS/n73/9eq/bl5eVBTlS9iooK7rjjDt5++202btzI66+/zsaNG09q9/zzz9OiRQtycnK4++67ue+++wDnYqTz589nw4YNLF68mJ/85CdUVFTUOUdGRgY33nhjtW38Zajt63n00Uf5/PPPWb9+PZ07d+bJJ5+sc05jTGhY4amlOXPmsHXrVq666ioeffRR9u/fzzXXXEO/fv0YNmzY8dshzJgxg8mTJ3PRRRcxefJkZsyYwZQpUxg+fDhdunRhwYIF3HvvvfTt25cxY8ZQVlYGOOfm7Nu3D4DVq1czYsSIkzL85z//YejQoaSkpDBq1Ch2795dbeZPP/2UpKQkunfvTv369Zk4cSJvvfXWSe3eeustpkyZAsD48eN59913UVXeeustJk6cSIMGDejWrRtJSUl8+umnrFq1in79+lFSUsKRI0fo06fPaZ875C9DbV9P06ZNAecK38XFxYjIaeUxxgSPFZ5a+vnPf06HDh1Yvnw5d999N7/73e9ISUlh/fr1/OlPf+KGG2443nbjxo288847vP766wBs2bKFZcuWkZaWxg9/+ENGjhxJZmYm8fHx/O9//6t1hosvvphPPvmEtWvXMnHiRGbNmgU4herWW289qX1+fv4JV8r2d7Vqz3axsbE0a9aMgoICv8MPHjyYsWPH8pvf/IZ7772XH/7wh5x//vm1fh2++MtQl9dz00030b59ezZt2sSdd955WnmMMcET1CsXiMgY4HEgBnhOVR/y6n8J8BjQD5ioqm949JsC/MZ9+gdVfSmYWevqgw8+4M033wTgsssuo6CggEOHDgHOVQbi4+OPt73qqquIi4ujb9++VFRUMGbMGAD69u3r97cMX/Ly8pgwYQI7d+6ktLSUbt26AZCamspzzz0XoFdWOw888ACDBw+mYcOGPPHEEz7bDB06lGPHjlFUVMT+/fuPXyZo5syZjB49OuCZXnjhBSoqKrjzzjv5xz/+wU033RTwaRhjTl/QvvGISAzwFHAV0Bu4XkR6ezX7CrgReM1r2JbA74ChwBDgdyLSIlhZA837KtZVV6WuV68ecXFxx3cD1atX7/jvQLGxscfvIlpSUuJzvHfeeSfTpk0jMzOTZ555xm+7KomJicdvSgf+r1bt2a68vJyDBw/SqlWraocvKCigqKiIw4cP+82xcuVK1q1bx3PPPcfYsWNZt24d69at81l0/GWo6+uJiYlh4sSJxz8UGGMiTzB3tQ0BclR1q6qWAvOBEy6VrKq5qroeqPQadjSwVFX3q+oBYCkwJohZ62z48OHHDzLIyMigdevWx39nOBVdu3ZlzZo1AH43mgcPHjy+oX3ppZq/AA4ePJjs7Gy2bdtGaWkp8+fPP+FCpVXGjh17fHxvvPEGl112GSLC2LFjmT9/PseOHWPbtm1kZ2czZMgQAG677TZ+//vfM2nSJJ8HAtSVvwy1eT2qSk5ODuD8xpOWlsa555572pmMMcERzF1ticAOj+d5ON9gTnXYkz6qi8hUYCpAu3btyMjIOKWgtVFUVERJSQkffvghzZo1Y9SoUcyaNYsXX3yRBg0acM8995CRkUFubi7x8fHHs3g/r6io8Nlv7NixTJ06lUaNGjFgwAAKCwvJyMhg06ZN5Ofnk5GRwfe+9z2+853vkJCQwMCBA4+3ycrKIi0t7fjVsD3nw9SpU7nkkkuorKzkqquuYu/evWRkZDBv3jySk5O56KKL6NmzJ//6179ITEykadOm/Pa3vz0+jsGDB9O9e3diYmK44447eP/990lPT6ewsJAOHTrQrl07pk2bxpw5cxg4cKDPebdu3Tp27dp10vvjmdVfhn379vHwww8fP3Ta1+tZvnw5P/3pTzl69CiqSo8ePbj77rsDtjx4z9NIZTkDL1qyRkvO41Q1KH/AeJzfdaqeTwae9NP2RWC8x/NfAL/xeP5b4BfVTW/QoEEaTMuXLw/q+AMlWnKqRk9WyxlY0ZJTNXqynk5OYLUGqQ74+wvmrrZ8oJPH845ut2APa4wxJoIFs/CsAnqKSDcRqQ9MBNJqOWw6cKWItHAPKrjS7WaMMSbKBa3wqGo5MA2nYHwJ/FNVN4jIgyIyFkBEBotIHnAd8IyIbHCH3Q/8Hqd4rQIedLsZY4yJckE9j0dVFwGLvLo94PF4Fc5uNF/DzgPmBTOfMcaY0LMrFxhjjAkpKzzGGBPlKr2uaxjprPAYY0yUyt13hDtfX8vc9cfCHaVOgvobjzHGmMDbc6iEJ5ZlM//THcTF1GNUp3qoatRcld0KjzHGRImDxWU8s2IL8z7cRnmFcv2Qztx5WRIbP/skaooOWOExxpiIV1JWwUsf5fLXjC0cLC5jbP8O3HNlL7q0ci5IfPLtHSObFR5jjIlQ5RWV/GtNHo+/k82uQyWMSG7D9NHJ9OnQLNzRTosVHmOMiTCqyttf7OLh9Cy27jtCSufmPDZxAMO6t6p54ChghccYYyLIB9n7mJW+ifV5B+nZNoG5kwdxRe92UfUbTk2s8BhjTARYn1fIrMVZfJCzj8Tm8Tx8XX+uTUkkpt6ZU3CqWOExxpgw2rK3iEeWZLEocxctG9fnt9/uzaShnWkYFxPuaEFjhccYY8Jg58Finng3m3+uzqNhbD3uurwntw7vRpOGceGOFnRWeIwxJoQKj5bydMYWXvwol0pVJg/rwrTLkmid0CDc0ULGCo8xxoTA0dJyXvgwl7+t2ELRsXKuTUnk7lG96NSyUbijhZwVHmOMCaKyikrmr9rBE+9ms/fwMUad15ZfjE7m3PZNwx0tbKzwGGNMEFRWKv/N3MkjS7LYXnCUwV1b8PSkgaR2bRnuaGFnhccYYwJIVVmxeS+zFmexcechzm3fhHk3pjIyue0ZdS7O6Qhq4RGRMcDjQAzwnKo+5NW/AfAyMAgoACaoaq6I1AeeAVKBSuAuVc0IZlZjjKmthWvzmZ2exdeFxXRoHs/00clck5LIZ18dYObbm1i5bT+dWsbz2IQBjO3fgXpn4Lk4pyNohUdEYoCngCuAPGCViKSpquf17G4BDqhqkohMBGYCE4AfAahqXxFpC7wtIoNVtTJYeY0xpjYWrs3n/gWZFJdVAJBfWMx9b67n+Q+2kZl/kNYJ9fm/sX24fkhn6sfaLc98CeY3niFAjqpuBRCR+cA4TryQ6jhghvv4DeBJcb6L9gaWAajqHhEpxPn282kQ8xpjTI1mp2cdLzpVjpVX8kX+Qe65ohc3X9yNxg3sV4zqiAbplqkiMh4Yo6q3us8nA0NVdZpHmy/cNnnu8y3AUOC7ON+Urgc6AWuBW1T1Ta9pTAWmArRr127Q/Pnzg/JaAIqKikhISAja+AMlWnJC9GS1nIEVLTnBd9bM/IMcLYdP99Rj3X5nF9qAVsqQNpUM7RKeq0afzjwdOXLkGlVNDXCkakVqWZ4HnAesBrYDHwEV3o1UdS4wFyA1NVVHjBgRtEAZGRkEc/yBEi05IXqyWs7AipaccHLWI8fKufPP73K4pPyEdmv2CbvKG3PflBGEQzTNUwhu4cnH+bZSpaPbzVebPBGJBZoBBep8Dbu7qpGIfARsDmJWY4zxq7S8ktdWbufJ5TkcLimnnkClx86i+LgYpo9ODl/AKBPMwrMK6Cki3XAKzETgB15t0oApwMfAeGCZqqqINMLZDXhERK4Ayr0OSjDGmKCrVOXfa/N4ZMlm8g4UM6x7S5694Vy2Fxz1eVSbqZ2gFR5VLReRaUA6zuHU81R1g4g8CKxW1TTgeeAVEckB9uMUJ4C2QLqIVOIUrcnBymmMMd5UlWWb9vC7D4vJK/qcPh2a8qdr+zK8Z2tEhJTOLazQnIag/sajqouARV7dHvB4XAJc52O4XMC+txpjQm5V7n5mLd7EqtwDtGsk/OX6FL7V9xw7FyeAIvXgAmOMCalNuw4xe3EW727aQ9smDfjjtefT7shWRvXvEO5oZxwrPMaYs9qO/Ud5dOlm/r0un4QGsdw7JpmbLuxGfP0YMjK2hTveGckKjzHmrLT38DGeWp7D31dup54It13Sg9sv7U7zRvXDHe2MZ4XHGHNWOVxSxrPvb+O597dyrLyS76d24q7Le9K+WcNwRztrWOExxpwVSsoqePWT7Ty1PIcDR8v4Vr9zuOeKXnRvEx1XUTiTWOExxpzRKiqVBZ/l8dg72eQXFjO8Z2umj06mX8fm4Y521rLCY4w5I6kqSzbu5uH0LLL3FNG/YzNmje/HRUmtwx3trGeFxxhzxvlkawEzF29i7VeFdG/TmKcnDWTM+e3tRmwRwgqPMeaM8UX+QWanZ7Fi817aN23IzO/15XsDOxIbY/fFiSRWeIwxUS933xEeWbqZ/3z+Nc3i4/jV1edywwVdaRgXE+5oxgcrPMaYqLXnUAlPLMtm/qc7iIupx7SRSfzoku40i48LdzRTDSs8xpioc7C4jGdWbGHeh9sor1CuH9KZOy9Pom0TOxcnGljhMcZEjZKyCl76KJe/ZmzhYHEZ4wZ04OdX9KJLq8bhjmbqwAqPMSbilVdU8sYa51ycXYdKuLRXG+4dk0yfDuG51bQ5PVZ4jDERS1V5+4tdPLwki617j5DSuTmPTRzAsO6twh3NnAYrPMaYiPRhzj5mLt7E+ryD9GybwNzJg7iidzs7F+cMYIXHGBNR1ucVMmtxFh/k7COxeTwPX9efa1MSibEbsZ0xrPAYYyLClr1FzFmymf9l7qRl4/r89tu9mTS0s52LcwYKauERkTHA40AM8JyqPuTVvwHwMjAIKAAmqGquiMQBzwED3Ywvq+qfg5nVGBMeuw6W8Pi7m/nn6jwaxtbjrst7cuvwbjRpePK5OAvX5jM7PYuvC4vp0Dye6aOTuSYlMQypzekIWuERkRjgKeAKIA9YJSJpqrrRo9ktwAFVTRKRicBMYAJwHdBAVfuKSCNgo4i8rqq5wcprjAmtwqOlPL1iCy9+mEulKpOHdWHaZUm0Tmjgs/3CtfncvyCT4rIKAPILi7l/QSaAFZ8oE8xvPEOAHFXdCiAi84FxgGfhGQfMcB+/ATwpzi+HCjQWkVggHigFDgUxqzEmRI6VK08tz+FvK7ZQdKyca1MSuXtULzq1bFTtcLPTs44XnSrFZRXMTs+ywhNlRFWDM2KR8cAYVb3VfT4ZGKqq0zzafOG2yXOfbwGGAgeBV4DLgUbA3ao618c0pgJTAdq1azdo/vz5QXktAEVFRSQkRP4No6IlJ0RPVssZGOWVynt55SzMPsahMmFAmxi+16s+nZrU7gKemfkH/fbrmxic83kifZ5WOZ2cI0eOXKOqqQGOVK1IPbhgCFABdABaAO+LyDtV356quMVoLkBqaqqOGDEiaIEyMjII5vgDJVpyQvRktZynp7JS+W/mTh5ZksX2glJ6tYhh3oShpHZtWafx/PqhZeQXFp/UPbF5PHdOGhGgtCeK1HnqLVpyVgnmtcLzgU4ezzu63Xy2cXerNcM5yOAHwGJVLVPVPcCHQEgrsjHm9KgqKzbv5TtPfsBPX19LfFwM825M5f4hDetcdACmj04m3usIt/i4GKaPTg5UZBMiwfzGswroKSLdcArMRJyC4ikNmAJ8DIwHlqmqishXwGXAKyLSGBgGPBbErMaYAPrsqwPMWryJT7bup1PLeB6bMICx/TtQr56QsevLUxpn1e84dlRb9Ata4VHVchGZBqTjHE49T1U3iMiDwGpVTQOexykuOcB+nOIEztFwL4jIBkCAF1R1fbCyGmMCI3v3YWanZ7Fk425aJ9Tn/8b24fohnakfG5idK9ekJFqhOQME9TceVV0ELPLq9oDH4xKcQ6e9hyvy1d0YE5nyC4t5bOlm3vwsj0b1Y7nnil7cfHE3GjeI1J+RTTjZUmGMOWX7j5Ty1PIcXvl4OwjcfFE3fjIyiZaN64c7molgVniMMXV25Fg5z3+wjbnvbeVoaTnjB3XkrlG9SGweH+5oJgpY4THG1FppeSWvrdzOk8tz2FdUyug+7Zg+Opmktk3CHc1EESs8xpgaVVYqb32ez5ylm9mxv5hh3Vvy7A3nktK5RbijmShkhccY45eqsjxrD7MWZ7Fp12H6dGjKSzf35ZKere2+OOaUWeExxvi0Onc/MxdvYlXuAbq2asRfrk/hW33PoZ7dF8ecploVHhHpCfwZ6A00rOquqt2DlMsYEyabdh3i4fQs3vlyD22aNOAP15zPhMGdiIsJ5oVOzNmktt94XgB+BzwKjARuIriX2zHGhNiO/Ud5dOlm/r0un4QGsdw7JpkbL+xKo/q2Y8QEVm2XqHhVfVdERFW3AzNEZA3wQE0DGmMi276iYzy5LIe/r9xOPRGmXtKdH1/ag+aN7FwcExy1LTzHRKQekO1eBicfiPxrhRtj/DpcUsaz72/jufe3cqy8ku+nduKuy3vSvlnDmgc25jTUtvDchXNfnJ8Cv8fZ3TYlWKGMMcFTUlbB31d+xVPLc9h/pJRv9T2Hn1/Zix5t7LOkCY1qC4+INASaqOoqt1MRcJOItMXuCGpMVKmoVBZ8lsdj72STX1jM8J6tmT46mX4dm4c7mjnL1PSN5wlgMbDAq/tFwJXAj4MRyhgTOKrK0o27mZ2eRfaeIvp1bMas8f24KKl1uKOZs1RNhWeQqk717qiq/xaRPwQpkzEmQD7ZWsDMxZtY+1Uh3ds05ulJAxlzfns7+dOEVU2Fp1E1/exwamMi1IavDzJrcRYrNu+lfdOGPPTdvowf1JFYOxfHRICaCs8eERmiqp96dhSRwcDe4MUyxpyK3H1HmLN0M2mff02z+Dh+dfW53HBBVxp63TLamHCqqfBMB/4pIi8Ca9xuqcANfHO3UGNMmO05VMLLG47x3pIVxMXU446RPZh6SQ+axceFO5oxJ6m28KjqpyIyBLgDuNHtvAEYqqp7ahq5iIwBHse59fVzqvqQV/8GwMvAIKAAmKCquSIyCafoVekHDFTVdbV5UcacLQ4WlzH3vS3M+yCX0vIKfjC0C3delkTbpoE/F2fh2nxmp2fxdWExHZrHM310st2G2pySmg6nngu8DTysqofrMmIRiQGeAq4A8oBVIpKmqhs9mt0CHFDVJBGZCMzEKT5/B/7ujqcvsNCKjjHfKCmr4KWPcvlrxhYOFpcxtn8HLm52gO9ffX5QprdwbT73L8ikuKwCcG51ff+CTAArPqbOavql8XmgP7BIRN4VkftEpH8txz0EyFHVrapaCswHxnm1GQe85D5+A7hcTj7c5np3WGPOeuUVlcz/9CtGzM7gz29vIqVzc/7304t54voU2jYK3oEDs9OzjhedKsVlFcxOzwraNM2ZS1S1dg1FWuGcu3MV0BdYCyxW1X/6aT8eGKOqt7rPJ+Psopvm0eYLt02e+3yL22afR5stwDhV/cLHNKYCUwHatWs3aP784NWnoqIiEhIi/8zuaMkJ0ZM1EnKqKqt3V/Bmdim7jig9mtXjuuT6nNvym4MGgpkzM/+g3359E5vVaVyRMD9rK1qynk7OkSNHrlHV1ABHqlatLzurqgXA68Dr7reS6UBSsIIBiMhQ4KivouNmmgvMBUhNTdURI0YELUtGRgbBHH+gREtOiJ6s4c75Yc4+Zi7exPq8o/Rsm8CD303mit7tTjoXJ5g5f/3QMvILi0/qntg8njsn1W2a4Z6fdREtWaMlZ5VTut65qqqITFPVztU0ywc6eTzv6Hbz1SZPRGKBZjgHGVSZiFPsjDnrrM8rZNbiLD7I2Udi83hmj+/Hdwd2JCYMN2KbPjr5hN94AOLjYpg+OjnkWUz0q+nggvX+egHtahj3KqCniHTDKTATgR94tUnDudjox8B4YJm6+/7cq2F/Hxhew3SMOaNs2VvEnCWb+V/mTlo2rs9vv92bSUM7h/VcnKoDCOyoNhMINX3jaQeMBg54dRfgo+oGVNVy9xYK6TiHU89T1Q0i8iCwWlXTcA5eeEVEcoD9nHhu0CXADlXdWutXY0wU23WwhMff3cw/V+fRMLYed13ek1uHd6NJw8g4F+ealEQrNCYgaio8/wUSfB3KLCIZNY1cVRcBi7y6PeDxuAS4zs+wGcCwmqZhTLQrPFrK0yu28OKHuVSqMnlYF6ZdlkTrhAbhjmZMUNR0Aukt1fTz3m1mjKmDo6XlvPBhLn9bsYWiY+Vcm5LI3aN60alldZdINCb62c3UjQmxsopK/rFqB4+/m83ew8cYdV5bfjE6mXPbNw13NGNCwgqPMSFSWan8N3MnjyzJYnvBUQZ3bcHTkwaS2rVluKMZE1JWeIwJMlXlvex9zFq8iQ1fH+Lc9k2Yd2MqI5Pb2n1xzFnJCo8xQfTZVweYtXgTn2zdT6eW8Tw2YQBj+3egXhjOxTEmUljhMSYIsncfZnZ6Fks27qZ1Qn3+b2wfrh/SmfqxdiM2Y6zwGBNA+YXFPLZ0M29+lkej+rHcc0Uvbr64G40b2KpmTBVbG4wJgP1HSnlqeQ6vfLwdBG6+qBs/GZlEy8b1wx3NmIhjhceY03DkWDnPf7CNue9t5WhpOeMHdeSuUb1IbB4f7mjGRCwrPMacgtLySl5buZ0nl+ewr6iU0X3aMX10Mkltm4Q7mjERzwqPMXVQUamkfZ7PI0s2k3egmGHdW/LsDeeS0rlFuKMZEzWs8BhTC6rKsk17mJ2exaZdh+nToSl/urYvw3u2tnNxjKkjKzzG1GDzgQqeeuZjVuUeoGurRvzl+hS+1fccOxfHmFNkhccYPzbtOsTsxVm8u6mEtk2UP157Pt9P7URcjJ2LY8zpsMJjjJcd+4/y6NLN/HtdPgkNYhnfM47fTx5JfP3w3YjNmDOJFR5jXPuKjvHkshz+vnI79US47ZIe3H5pd9Z9+pEVHWMCyAqPOesdLinj2fe38dz7WzlWXsn3Uztx1+U9ad+sYbijGXNGssJjzlolZRX8feVXPLU8h/1HSvlWv3O454pedG+TEO5oxpzRglp4RGQM8DgQAzynqg959W8AvAwMAgqACaqa6/brBzwDNAUqgcHurbKNOS0VlcqCz/J47J1s8guLGd6zNdNHJ9OvY/MT2i1cm8/s9CwmdjrMrx9axvTRyVyTkhjyvFU5vi4spkPz+LDlMCZQglZ4RCQGeAq4AsgDVolImqpu9Gh2C3BAVZNEZCIwE5ggIrHAq8BkVf1cRFoBZcHKas4OqsrSjbuZnZ5F9p4i+ndsxqzx/bgoqfVJbReuzef+BZkUl1VAJ+fin/cvyAQI6Ub/hByEL4cxgRTM40KHADmqulVVS4H5wDivNuOAl9zHbwCXi3M23pXAelX9HEBVC1S1IohZzRnuk60FfPfpj5j6yhoqVHl60kAW3nGRz6IDMDs96/jGvkpxWQWz07NCETfichgTSKKqwRmxyHhgjKre6j6fDAxV1Wkebb5w2+S5z7cAQ4Ef4ux+awu0Aear6iwf05gKTAVo167doPnz5wfltQAUFRWRkBD5+/6jJSeEJuv2QxW8sbmMzH0VtGggXNMzjos7xBJTw8mfmfkHjz9uFw+7i7/p1zexWbDiVpvDm3eOaHnvoyUnRE/W08k5cuTINaqaGuBI1YrUgwtigYuBwcBR4F0RWaOq73o2UtW5wFyA1NRUHTFiRNACZWRkEMzxB0q05ITgZs3dd4Q5SzeT9vnXNIuP41dX9+SGC7rSMK52h0X/+qFl5Bc61eaevuU8kumsKonN47lzUnAy15TDk68c0fLeR0tOiJ6s0ZKzSjB3teUDnTyed3S7+Wzj/q7TDOcggzzgPVXdp6pHgUXAwCBmNWeIPYdK+M3CTEbNWcHSjbuZNjKJ9+4dydRLetS66ABMH51MvFf7+LgYpo9ODnTkqMhhTCAF8xvPKqCniHTDKTATgR94tUkDpgAfA+OBZaqqIpIO3CsijYBS4FLg0SBmNVHuYHEZc9/bwrwPcimrqOT6IZ258/Ik2jY5tXNxqn64d35LOUximI4m88xhR7WZM0XQCo+qlovINCAd53Dqeaq6QUQeBFarahrwPPCKiOQA+3GKE6p6QETm4BQvBRap6v+CldVEr5KyCl76KJe/ZmzhYHEZ4wZ04OdX9KJLq8anPe5rUhK5JiWRjIyMkO5e85fDmDNFUH/jUdVFOLvJPLs94PG4BLjOz7Cv4hxSbcxJyisqeWONcy7OrkMljEhuw/TRyfTpELof/o0xpyZSDy4wxidV5e0vdvHwkiy27j3CwM7NeWziAIZ1bxXuaMaYWrLCY6LGhzn7mLl4E+vzDtKrXQLP3pDKqPPa2o3YjIkyVnhMxFufV8isxVl8kLOPxObxPHxdf65NSazxXBxjTGSywmMi1pa9RcxZspn/Ze6kZeP6PPDt3kwa1pkGsXaLAmOimRUeE3F2HizmiXez+efqPBrG1uOuy3ty6/BuNGkYF+5oxpgAsMJjIkbh0VKeztjCix/lUqnK5GFdmHZZEq0TGoQ7mjEmgKzwmLA7WlrOCx/m8rcVWyg6Vs61KYncPaoXnVo2Cnc0Y0wQWOExYVNeqbzyyXaeeDebvYePMeq8dkwfnUxy+ybhjmaMCSIrPCbkKiuV/2bu5A8fFLPn6BcM7tqCpycNJLVry3BHM8aEgBUeEzKqyorNe5m1OIuNOw/RqUk9XrgxlRHJbexcHGPOIlZ4TEh89tUBZr69iZXb9tOpZTyPTRhA08LNjDy3bbijGWNCzAqPCars3YeZnZ7Fko27aZ1QnwfH9WHi4M7Uj61HRkZ2uOMZY8LACo8JivzCYh5dupkFn+XRuH4s91zRi5sv7kbjBrbIGXO2s62ACaj9R0p5ankOr3y8HQRuubgbPx6RRMvG9cMdzRgTIazwmIA4cqyc597fxrPvb+VoaTnjB3XkZ6N60aF5fLijGWMijBUec1qOlVfw+sqv+MuyHAqOlDKmT3t+MboXSW3tXBxjjG9WeMwpqahU3lqXz5ylm8k7UMwF3Vtx75hkUjq3CHc0Y0yEC2rhEZExwOM4t75+TlUf8urfAHgZGAQUABNUNVdEugJfAllu009U9fZgZg2nhWvzmZ2exdeFxXRoHs/00clhudVxbXKoKss27WF2ehabdh2mT4em/Onavgzv2Tpg5+JEyvwwxgRH0AqPiMQATwFXAHnAKhFJU9WNHs1uAQ6oapKITARmAhPcfltUdUCw8kWKhWvzuX9BJsVlFYBzNNj9CzIBQrqxrU2OVbn7mfn2JlZvP0DXVo148gcpXH3+OdQL4H1xImV+GGOCp14Qxz0EyFHVrapaCswHxnm1GQe85D5+A7hczrJT2GenZx3fyFYpLqtgdnqWnyFCn+PLnYe45cVVXPe3j/lq/1H+eO35LP35pXy7X4eAFp2achhjzgyiqsEZsch4YIyq3uo+nwwMVdVpHm2+cNvkuc+3AEOBBGADsBk4BPxGVd/3MY2pwFSAdu3aDZo/f35QXgtAUVERCQkJAR9vZv5Bv/36Jjar8/hONaevHAdL4cPd9dhUWI/4WPhW9zhGdYmjQUxgio2vrIGeH4EQrPc+0Cxn4EVL1tPJOXLkyDWqmhrgSNWK1IMLdgKdVbVARAYBC0Wkj6oe8mykqnOBuQCpqak6YsSIoAXKyMggGOP/9UPLyC8sPql7YvN47pxU9+mdak5/OQBuv7QHP760B80aBfZGbL6yBnp+BEKw3vtAs5yBFy1ZoyVnlWDuassHOnk87+h289lGRGKBZkCBqh5T1QIAVV0DbAF6BTFr2EwfnUx83Im3co6Pi2H66OSQ52gYe+LiEFNPmPGd3vzyqnMDXnSqyxEJ88MYEzzB/MazCugpIt1wCsxE4AdebdKAKcDHwHhgmaqqiLQB9qtqhYh0B3oCW4OYNWyqfjAP51FcJWUV7Cs6RozH7zVtmzTgV1efF/If9CNhfhhjgitohUdVy0VkGpCOczj1PFXdICIPAqtVNQ14HnhFRHKA/TjFCeAS4EERKQMqgdtVdX+wsobbNSmJYdmwlldUsmBtPo8t3czXB0sY3rM1944+l74dw/NbSpVwzQ9jTGgE9TceVV0ELPLq9oDH4xLgOh/DvQm8GcxsZzNVZcnG3cxOzyJnTxH9Ozbj4ev6c2FS63BHM8acBSL14AITJB9vKWDm4k2s21FI9zaN+dsPBzK6T3u7EZsxJmSs8Jwlvsg/yKz0LN7bvJdzmjVk5vf68r2BHYmNCebxJcYYczIrPGe43H1HeGTpZv7z+dc0bxTHr68+j8kXdKGh15FjxhgTKlZ4zlB7DpXw+LvZ/GPVDuJi6jFtZBJTL+1O04ahOSzaGGP8scJzhjlSpsxavIl5H26jvEL5wdDOTLssibZNGoY7mjHGAFZ4zhglZRW8+FEuf3nvKEfKtjBuQAd+fkUvurRqHO5oxhhzAis8Ua68opJ/rcnjsXc2s/vQMfq1ieHP119Anw7hPRfHGGP8scITpVSVRZm7eGRJFlv3HWFg5+Y8MTGF4q8yregYYyKaFZ4o9EH2PmYu3kRm/kF6tUvg2RtSGXVeW0SEjK/Cnc4YY6pnhSeKfL6jkFnpm/gwp4DE5vE8cl1/rklJPOEaa8YYE+ms8ESBLXuLeGRJFosyd9GycX0e+HZvJg3rTINYOxfHGBN9rPBEsJ0Hi3n8nWz+tSaPhrH1+Nmontw6vDsJDextM8ZEL9uCRaDCo6U8nbGFFz/KRRVuuKALd4xMonVCg3BHM8aY02aFJ4IcLS3nhQ9z+duKLRQdK+e7KR352aiedGrZKNzRjDEmYKzwRICyikrmr9rBE+9ms/fwMUad147po5NJbt8k3NGMMSbgrPCEUWWl8p/1XzNn6Wa2FxxlcNcWPD1pIKldW4Y7mjHGBI0VnjBQVVZs3susxVls3HmIc9s34YUbBzMiuY3dF8cYc8azwhNia7YfYNbiTazctp9OLeN5bMIAxvbvQD07F8cYc5YIauERkTHA40AM8JyqPuTVvwHwMjAIKAAmqGquR//OwEZghqo+HIyMC9fmMzs9i68Li+nQPJ7po5O5JiUx4NPZvPsws9OzWLpxN60T6vPguD5MHNyZ+rF2IzZjzNklaIVHRGKAp4ArgDxglYikqepGj2a3AAdUNUlEJgIzgQke/ecAbwcr48K1+dy/IJPisgoA8guLuX9BJkDAik/egaM89k42Cz7Lo3H9WO65ohc3X9yNxnYujjHmLBXMrd8QIEdVtwKIyHxgHM43mCrjgBnu4zeAJ0VEVFVF5BpgG3AkWAFnp2cdLzpVissqmJ2eddqFp6DoGE8t38Krn2wHgVsu7saPRyTRsnH90xqvMcZEO1HV4IxYZDwwRlVvdZ9PBoaq6jSPNl+4bfLc51uAoUAJsBTn29IvgCJfu9pEZCowFaBdu3aD5s+fX6eMmfkH/fbrm3jiFZ6LiopISEiocZzF5cqS3DLe3lbGsQq4ODGWa5LiaBUfml1qtc0ZCaIlq+UMrGjJCdGT9XRyjhw5co2qpgY4UrUidX/PDOBRVS2q7igvVZ0LzAVITU3VESNG1Gkiv35oGfmFxSd1T2wez52TThxXRkYG1Y3/WHkFr638iiffz6HgSBlj+rTnF6N7kdQ2tOfi1JQzkkRLVssZWNGSE6Ina7TkrBLMwpMPdPJ43tHt5qtNnojEAs1wDjIYCowXkVlAc6BSREpU9clABpw+OvmE33gA4uNimD46udbjqKhU3lqXz5ylm8k7UMwF3Vtx75hkUjq3CGRUY4w5YwSz8KwCeopIN5wCMxH4gVebNGAK8DEwHlimzr6/4VUNRGQGzq62gBYd+OYAglM5qk1VeffLPcxOzyJr92H6dGjKn67ty/Cere1cHGOMqUbQCo+qlovINCAd53Dqeaq6QUQeBFarahrwPPCKiOQA+3GKU0hdk5JY5wMJVuXuZ+bbm1i9/QBdWzXiyR+kcPX559i5OMYYUwtB/Y1HVRcBi7y6PeDxuAS4roZxzAhKuFPw5c5DzE7PYtmmPbRt0oA/Xns+30/tRFyMnYtjjDG1FakHF0SUrwqO8sz6Ej5Jf58mDWK5b8y53HhhV+Lr243YjDGmrqzwVGPv4WM8uSyb1z79ClS57ZIe/PjSHjRrFBfuaMYYE7Ws8PhwqKSM597bynMfbONYeSUTBndicPxerh1zbrijGWNM1LPC46GkrIJXP9nOU8tzOHC0jG/1O4d7ruhF9zYJZGRkhDueMcacEazwAOUVlSxYm89jSzfz9cEShvdszb2jz6Vvx2Y1D2yMMaZOzvrC81XBUW5+aRU5e4ro37EZD1/XnwuTWoc7ljHGnLHO+sJzTvOGdGnZiF9c2YvRfdrbyZ/GGBNkZ33hiYupx/M3Dg53DGOMOWvYmY/GGGNCygqPMcaYkLLCY4wxJqSs8BhjjAkpKzzGGGNCygqPMcaYkLLCY4wxJqSs8BhjjAkpce40Hf1EZC+wPYiTaA3sC+L4AyVackL0ZLWcgRUtOSF6sp5Ozi6q2iaQYWpyxhSeYBOR1aqaGu4cNYmWnBA9WS1nYEVLToierNGSs4rtajPGGBNSVniMMcaElBWe2psb7gC1FC05IXqyWs7AipacED1ZoyUnYL/xGGOMCTH7xmOMMSakrPAYY4wJKSs8HkSkk4gsF5GNIrJBRO7y0WaEiBwUkXXu3wNhyporIpluhtU++ouIPCEiOSKyXkQGhiFjssd8Wicih0TkZ15twjY/RWSeiOwRkS88urUUkaUiku3+b+Fn2Clum2wRmRKGnLNFZJP73v5bRJr7Gbba5SQEOWeISL7H+3u1n2HHiEiWu7z+Mpg5q8n6D4+cuSKyzs+woZynPrdJkbic1omq2p/7B5wDDHQfNwE2A7292owA/hsBWXOB1tX0vxp4GxBgGLAyzHljgF04J6tFxPwELgEGAl94dJsF/NJ9/Etgpo/hWgJb3f8t3MctQpzzSiDWfTzTV87aLCchyDkD+EUtlo0tQHegPvC593oXiqxe/R8BHoiAeepzmxSJy2ld/uwbjwdV3amqn7mPDwNfAonhTXXKxgEvq+MToLmInBPGPJcDW1Q1mFeXqBNVfQ/Y79V5HPCS+/gl4Bofg44GlqrqflU9ACwFxoQyp6ouUdVy9+knQMdgTb+2/MzP2hgC5KjqVlUtBebjvA9BU11WERHg+8DrwcxQG9VskyJuOa0LKzx+iEhXIAVY6aP3BSLyuYi8LSJ9QpvsOAWWiMgaEZnqo38isMPjeR7hLaIT8b8iR8L8rNJOVXe6j3cB7Xy0ibR5ezPOt1tfalpOQmGau0twnp9dQpE2P4cDu1U120//sMxTr21SNC6nx1nh8UFEEoA3gZ+p6iGv3p/h7C7qD/wFWBjieFUuVtWBwFXAHSJySZhy1EhE6gNjgX/56B0p8/Mk6uyviOjzDUTk10A58Hc/TcK9nDwN9AAGADtxdmFFuuup/ttOyOdpddukaFhOvVnh8SIicThv8N9VdYF3f1U9pKpF7uNFQJyItA5xTFQ13/2/B/g3zu4KT/lAJ4/nHd1u4XAV8Jmq7vbuESnz08Puql2S7v89PtpExLwVkRuBbwOT3I3PSWqxnASVqu5W1QpVrQSe9TP9iJifACISC3wX+Ie/NqGep362SVGznPpihceDu2/3eeBLVZ3jp017tx0iMgRnHhaELiWISGMRaVL1GOeH5i+8mqUBN4hjGHDQ46t5qPn9BBkJ89NLGlB19M8U4C0fbdKBK0Wkhbvr6Eq3W8iIyBjgXmCsqh7106Y2y0lQef2ueK2f6a8CeopIN/fb8USc9yEcRgGbVDXPV89Qz9NqtklRsZz6Fe6jGyLpD7gY5yvremCd+3c1cDtwu9tmGrAB58ibT4ALw5Czuzv9z90sv3a7e+YU4Cmco4UygdQwzdPGOIWkmUe3iJifOMVwJ1CGs//7FqAV8C6QDbwDtHTbpgLPeQx7M5Dj/t0Uhpw5OPvvq5bTv7ltOwCLqltOQpzzFXf5W4+zsTzHO6f7/GqcI7a2BDunv6xu9xerlk2PtuGcp/62SRG3nNblzy6ZY4wxJqRsV5sxxpiQssJjjDEmpKzwGGOMCSkrPMYYY0LKCo8xxpiQssJjjDEmpKzwGONFRFp5XB5/l8dl/YtE5K9BmN6LIrJNRG73eD7eR7seVTkCncGYUIoNdwBjIo2qFuBcWwwRmQEUqerDQZ7sdFV9o4ZcW4ABVnhMtLNvPMbUkjg3rfuv+3iGiLwkIu+LyHYR+a6IzHJvELbYvb4WIjJIRFa4VzJOr8OtKS4RkY9EZKuvbz/GRDMrPMacuh7AZThX3n4VWK6qfYFi4Ftu8fkLMF5VBwHzgD/Wctzn4Fwu5dvAQ4EObkw42a42Y07d26paJiKZOHfRXOx2zwS6AsnA+cBS9zqoMTjXB6uNhepc0XmjiPi614oxUcsKjzGn7hiAqlaKSJl+c+HDSpx1S4ANqnrBqY7bJacX05jIYrvajAmeLKCNiFwAzn1VIuAOq8aEnRUeY4JEVUuB8cBMEfkc55L2F4Y1lDERwG6LYEyYiciLwH9rOpzao32RqiYEN5UxwWPfeIwJv4PA76tOIPWn6gRS4KRbiBsTTewbjzHGmJCybzzGGGNCygqPMcaYkLLCY4wxJqSs8BhjjAmp/weppBLqLrPgrgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df['1/Ca'] = round(1/df['Ca'], 2)\n",
    "\n",
    "# linear least-squares regression\n",
    "m, b, r_value, p_value, std_err = scipy.stats.linregress(df['Time'].iloc[1:], df['1/Ca'].iloc[1:])\n",
    "\n",
    "#plotting the results and annotating the plot\n",
    "fig, ax = plt.subplots()\n",
    "ax.scatter(df['Time'].iloc[1:], df['1/Ca'].iloc[1:])\n",
    "ax.plot(df['Time'].iloc[1:], m*df['Time'].iloc[1:] + b)\n",
    "ax.annotate('r^2: ' + str(\"{:.2f}\".format(r_value**2)), xy=(2.5, 0.105))\n",
    "ax.annotate('formula: ' + str(\"{:.2f}\".format(m)) + 'x + ' + str(\"{:.2f}\".format(b)), xy=(2.5, 0.10))\n",
    "plt.title('Linear least-squares regression for two sets of measurements.')\n",
    "plt.xlabel('Time [h]')\n",
    "plt.ylabel('1/Ca')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "fig.show()"
   ]
  }
 ],
 "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
}