play.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "import pathlib\n",
    "import pandas as pd\n",
    "import os\n",
    "from sklearn.linear_model import LinearRegression\n",
    "import seaborn as sns\n",
    "from sklearn.model_selection import train_test_split"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def calc_step(self, d):\n",
    "    return np.sqrt(4 * d * 2)\n",
    "\n",
    "def convert_d_from_nm2_ns_to_cm2_s(d: float):\n",
    "    return d * 1e-14 / 1e-9\n",
    "\n",
    "def convert_d_from_cm2_s_to_nm2_ns(d: float):\n",
    "    return d * 1e-9 / 1e-14\n",
    "\n",
    "def calc_d_from_step(step):\n",
    "        \"\"\"calculates the d value based on the step nm distance\"\"\"\n",
    "        return (step**2) / (4 * 2)\n",
    "\n",
    "\n",
    "def calc_step_rads(d):\n",
    "    return np.sqrt(2 * d * 2)\n",
    "\n",
    "\n",
    "def convert_d_from_rads2_s_to_rads2_ns(d: float):\n",
    "    \"\"\"converts time scale from s to ns\"\"\"\n",
    "    return d * 1e-9\n",
    "\n",
    "def convert_d_from_rads2_ns_to_rads2_s(d: float):\n",
    "    \"\"\"converts time scale from s to ns\"\"\"\n",
    "    return d / 1e-9\n",
    "\n",
    "def calc_rot_d_from_step(step):\n",
    "    \"\"\"calculate the d value for rotation based on the rotation mean in radians\n",
    "    Rotation diffusion coefficient: D = rad^2 / (2*t)\n",
    "    \"\"\"\n",
    "    return (step**2) / (2 * 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "file_str = \"src\\\\grana_model\\\\res\\log\\\\01072022_displacement_data.csv\"\n",
    "\n",
    "filename = f\"{pathlib.Path.cwd()}/{file_str}\"\n",
    "\n",
    "if os.path.exists(filename):\n",
    "    df = pd.read_csv(filename)\n",
    "else:\n",
    "    print(f\"{filename} does not exist\")\n",
    "\n",
    "df[\"displacement_sq\"] = df[\"displacement\"]**2\n",
    "df[\"rotation_sq\"] = df[\"rot_from_origin\"]**2\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "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>displacement</th>\n",
       "      <th>rot_from_origin</th>\n",
       "      <th>mass</th>\n",
       "      <th>rotation_scalar</th>\n",
       "      <th>diffusion_scalar</th>\n",
       "      <th>x</th>\n",
       "      <th>y</th>\n",
       "      <th>theta</th>\n",
       "      <th>displacement_sq</th>\n",
       "      <th>rotation_sq</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000</td>\n",
       "      <td>-3.417132</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>0.001785</td>\n",
       "      <td>1220.0</td>\n",
       "      <td>200.131083</td>\n",
       "      <td>295.020177</td>\n",
       "      <td>2.139292</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>11.676792</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2.0</td>\n",
       "      <td>0.011</td>\n",
       "      <td>-3.416959</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>0.001785</td>\n",
       "      <td>1220.0</td>\n",
       "      <td>200.139894</td>\n",
       "      <td>295.013698</td>\n",
       "      <td>2.139465</td>\n",
       "      <td>0.000121</td>\n",
       "      <td>11.675607</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4.0</td>\n",
       "      <td>0.017</td>\n",
       "      <td>-3.418452</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>0.001785</td>\n",
       "      <td>1220.0</td>\n",
       "      <td>200.148089</td>\n",
       "      <td>295.017132</td>\n",
       "      <td>2.137973</td>\n",
       "      <td>0.000289</td>\n",
       "      <td>11.685811</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>6.0</td>\n",
       "      <td>0.023</td>\n",
       "      <td>-3.413262</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>0.001785</td>\n",
       "      <td>1220.0</td>\n",
       "      <td>200.153368</td>\n",
       "      <td>295.025093</td>\n",
       "      <td>2.143163</td>\n",
       "      <td>0.000529</td>\n",
       "      <td>11.650355</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>8.0</td>\n",
       "      <td>0.045</td>\n",
       "      <td>-3.409736</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>0.001785</td>\n",
       "      <td>1220.0</td>\n",
       "      <td>200.173892</td>\n",
       "      <td>295.033322</td>\n",
       "      <td>2.146689</td>\n",
       "      <td>0.002025</td>\n",
       "      <td>11.626298</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   time  displacement  rot_from_origin    mass  rotation_scalar  \\\n",
       "0   0.0         0.000        -3.417132  1000.0         0.001785   \n",
       "1   2.0         0.011        -3.416959  1000.0         0.001785   \n",
       "2   4.0         0.017        -3.418452  1000.0         0.001785   \n",
       "3   6.0         0.023        -3.413262  1000.0         0.001785   \n",
       "4   8.0         0.045        -3.409736  1000.0         0.001785   \n",
       "\n",
       "   diffusion_scalar           x           y     theta  displacement_sq  \\\n",
       "0            1220.0  200.131083  295.020177  2.139292         0.000000   \n",
       "1            1220.0  200.139894  295.013698  2.139465         0.000121   \n",
       "2            1220.0  200.148089  295.017132  2.137973         0.000289   \n",
       "3            1220.0  200.153368  295.025093  2.143163         0.000529   \n",
       "4            1220.0  200.173892  295.033322  2.146689         0.002025   \n",
       "\n",
       "   rotation_sq  \n",
       "0    11.676792  \n",
       "1    11.675607  \n",
       "2    11.685811  \n",
       "3    11.650355  \n",
       "4    11.626298  "
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "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>displacement_sq</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2.0</td>\n",
       "      <td>0.000121</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4.0</td>\n",
       "      <td>0.000289</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>6.0</td>\n",
       "      <td>0.000529</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>8.0</td>\n",
       "      <td>0.002025</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9999995</th>\n",
       "      <td>1990.0</td>\n",
       "      <td>21.687649</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9999996</th>\n",
       "      <td>1992.0</td>\n",
       "      <td>21.399876</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9999997</th>\n",
       "      <td>1994.0</td>\n",
       "      <td>21.114025</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9999998</th>\n",
       "      <td>1996.0</td>\n",
       "      <td>20.912329</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9999999</th>\n",
       "      <td>1998.0</td>\n",
       "      <td>20.675209</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>10000000 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "           time  displacement_sq\n",
       "0           0.0         0.000000\n",
       "1           2.0         0.000121\n",
       "2           4.0         0.000289\n",
       "3           6.0         0.000529\n",
       "4           8.0         0.002025\n",
       "...         ...              ...\n",
       "9999995  1990.0        21.687649\n",
       "9999996  1992.0        21.399876\n",
       "9999997  1994.0        21.114025\n",
       "9999998  1996.0        20.912329\n",
       "9999999  1998.0        20.675209\n",
       "\n",
       "[10000000 rows x 2 columns]"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dflm = df.drop(columns=[\"rot_from_origin\", \"displacement\", \"mass\", \"rotation_scalar\", \"diffusion_scalar\", \"x\", \"y\", \"theta\", \"rotation_sq\"]).copy()\n",
    "dflm = dflm.dropna()\n",
    "dflm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Intercept: [-3.42339188], coef: [[0.04771214]]\n"
     ]
    }
   ],
   "source": [
    "# linear regression on data\n",
    "y = dflm[[\"displacement_sq\"]]\n",
    "x = dflm[[\"time\"]]\n",
    "\n",
    "x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=26)\n",
    "\n",
    "\n",
    "slr = LinearRegression()\n",
    "slr.fit(x_train, y_train)\n",
    "\n",
    "print(f\"Intercept: {slr.intercept_}, coef: {slr.coef_}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3000000,)\n"
     ]
    },
    {
     "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>y_test</th>\n",
       "      <th>yhat</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>NaN</td>\n",
       "      <td>80.072849</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>NaN</td>\n",
       "      <td>53.067779</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>NaN</td>\n",
       "      <td>31.883590</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.000529</td>\n",
       "      <td>11.939916</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>NaN</td>\n",
       "      <td>20.718950</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2999995</th>\n",
       "      <td>129.072321</td>\n",
       "      <td>33.696651</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2999996</th>\n",
       "      <td>127.938721</td>\n",
       "      <td>37.227349</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2999997</th>\n",
       "      <td>NaN</td>\n",
       "      <td>77.305545</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2999998</th>\n",
       "      <td>NaN</td>\n",
       "      <td>65.186662</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2999999</th>\n",
       "      <td>NaN</td>\n",
       "      <td>54.212870</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>3000000 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "             y_test       yhat\n",
       "0               NaN  80.072849\n",
       "1               NaN  53.067779\n",
       "2               NaN  31.883590\n",
       "3          0.000529  11.939916\n",
       "4               NaN  20.718950\n",
       "...             ...        ...\n",
       "2999995  129.072321  33.696651\n",
       "2999996  127.938721  37.227349\n",
       "2999997         NaN  77.305545\n",
       "2999998         NaN  65.186662\n",
       "2999999         NaN  54.212870\n",
       "\n",
       "[3000000 rows x 2 columns]"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# predict values\n",
    "yhat = slr.predict(x_test).flatten()\n",
    "print(yhat.shape)\n",
    "slr_diff = pd.DataFrame({\"y_test\": y_test[\"displacement_sq\"], \"yhat\": yhat}, index=[i for i in range(yhat.shape[0])])\n",
    "slr_diff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjyUlEQVR4nO3dd3hUdfr+8fcTCBA6KihFDDYUG0jUtaGIigUEC0VZBEHBRrFEYF1X+C67iqjYUEREUVCKYkBcRcW1rqhg6GVBZJEYQaWJBAjJ5/dHhmPCLwEmU86U+3VdXuQ5mcncnoSbkzOnmHMOERFJLCl+BxARkfBTuYuIJCCVu4hIAlK5i4gkIJW7iEgCUrmLiCSgA5a7mY03s41mtqTYskPM7AMzWxX4s06xzw0xs9VmttLM2kYquIiIlO1gttxfBi7bZ9lgYI5z7jhgTmDGzJoBXYGTAs951swqhC2tiIgclAOWu3PuU2DTPos7ABMCH08AOhZbPtk5t8s59z2wGjgzPFFFRORgVSzn8w53zuUCOOdyzaxeYHlDYG6xx60PLPv/mFkfoA9AtWrVWp5wwgnljBIfluVuo6Cw9LOBK6RYmZ8TkeRySsNaB/3Y+fPn/+Kcq1va58pb7mWxUpaV2lrOubHAWICMjAw3b968MEeJLU0Gv1P6iqBopanaRQRg3sNXHvRjzex/ZX2uvEfLbDCz+oEvXh/YGFi+Hjiy2OMaAT+W8zUSSoPaaaUub1g7rczPiYiUV3nLfSbQI/BxD2BGseVdzayymTUBjgO+Di1iYshs25S01JLvLaelViCzbVMy2zb1KZWIJKoD7pYxs9eBC4HDzGw98CDwMDDVzHoD64BOAM65pWY2FVgG7AHucM4VRCh7XOnYouith5GzV/Ljljwa1E6j9Ql1vVm7ZkQknCwWLvmbDPvc95WVncOQ6YvJy9e/fSLyh7XB7XOf75zLKO1z4X5DVcqQlZ3D0JlL2ZKXD0CKgQ6QEZFIUblHQVZ2DpnTFpJfrM1V7CISSbq2TBSMnL2yRLGLiESayj0KftyS53cEEUkyKvcoqJWW6ncEEUkyKvcosNLO2xURiSCVexRs2ZHvdwQRSTIq9ygo6/ICdapqd42IRIbKPQrKuvTALp3AJCIRouPco2DfSw/USkvFDDZrd42IRIi23KOkY4uGfDH4IkZ1ac6uPYUqdhGJKJV7lI2cvVLXkxERjv1lHatGduDK5Z9F5Otrt0yUZGXnMHL2SnJ0QpNI0nv0nVFct2QOAIft2BKR11C5R4GuACkiAI22buDzMb29uV/7TN5udkFEXkvlHgXaFSMiQ/49nr5fT/fmE+96g7xKVSL2eir3KNCuGJHkddjvm5n3THdv/uultzOxxRURf12Ve4RlZefoLksiSer2L6dy36evePOpAyazrUr1qLy2yj3CRs5eqWIXSTI1d25n0ZNdvXnk+d0ZfU6XqGZQuUeYLvcrklxuWPAu/5w92pvPuONVfq5eJ+o5VO4R1qB2mva5iySBKvk7WfH4dd48LqMDw9vc4lsencQUYaVdV0ZEEku75Z+WKPbz+o7ztdhBW+4RV/y6MtqCF0ksqQX5zB3dg0PztgEw/aTW3N3uHp9TFVG5i4iUwwVr5jNh2oPefEmv0ayqe5SPiUpSuUeYzk4VSSwphQW8+1I/mv6yDoCPm7SkZ6ehMXfLNZV7hOnsVJHEcfr65UyflOnNHbs/xoIGTX1MVDaVe4RpP7tIAnCO1ybfzznrFgGw5PBjaN9jFM5i95gUlXsE6exUkfh34sY1vPtSf2/u1mU4X6Q39y/QQVK5R5DOThWJb0/PGEH7FUXXW8+tfijn3/oieyrER23GR8o4pbNTReLTUZt/5JOxfby5b8e/MLvpOT4mCp7KPYJ0dqpI/Bn6wRh6fjsLgN0pFTll4BR2pVb2OVXwYvfdgASgs1NF4sfhv/3C2hHtvGLPvLw/x2dmxWWxg7bcI65KaooOhRSJcQM/n8TAL1735pMHTmV75ao+Jgqdyj1CdPKSSOyrlfcbC5+63puHt+7FuDOv8TFR+KjcI0QnL4nEth7z32bYh8978+n9JrGpai0fE4WXyj1CdKSMSGyqujuPZaM6efPoP3Vi5AU9fEwUGSGVu5ndBdxM0Xk6i4GbgKrAFCAdWAt0ds5tDillHNKRMiKx5+olHzHqnce9+ezbXiK3Zl0fE0VOuY+WMbOGQH8gwzl3MlAB6AoMBuY4544D5gTmpKMjZURiR9runawd0c4r9tdOa0v6oFkJW+wQ+qGQFYE0M6tI0Rb7j0AHYELg8xOAjiG+Rlzq2KIh17ZsSEpsXShOJOnc8tV0lo/640YaF908hr9c1s/HRNFR7t0yzrkcM3sUWAfkAe875943s8Odc7mBx+SaWb3Snm9mfYA+AI0bNy5vjJiVlZ3Dm/NzKNT1B0R8UbFgD6sf7VhiWfqgWf6E8UEou2XqULSV3gRoAFQzsz8f7POdc2OdcxnOuYy6dRPvVyMdLSPin2sXzylR7Ddd92BSFTuE9obqxcD3zrmfAcxsOnAOsMHM6ge22usDG8OQM+7oaBmR6DNXyPePXFViWZP7Zsb0pXkjJZT/43XAn8ysqpkZ0AZYDswE9h5X1AOYEVrE+NSgdprfEUSSykWrvy5R7HdfeRfpg2YlZbFDaPvcvzKzN4BvgT1ANjAWqA5MNbPeFP0D0Knsr5K4Mts2ZeCUBX7HEEl8zrH2kfYlFh17b1bcXJo3UkL6J80596Bz7gTn3MnOue7OuV3OuV+dc22cc8cF/twUrrDxpGOLhn5HEEl4LdcvK1Hsf2/dm/RBs5K+2EFnqEbMX7MW+x1BJKF98/SfqbtjizefcPcb7Eyt4l+gGJOcO6OiYNLcdX5HEElITX9ey9oR7bxiH3PWtaQPmqVi34e23CMgKztHt9cTiYCZEwZy6k+rvfnUAZPZVqW6j4lil8o9Au57Y6HfEUQSSqOtG/h8TG9vfvOk1tzT7h4fE8U+lXsE7C7QdrtIuDw/fThtV8315rNuf5kNNQ7zMVF8ULmLSEw69PctzH/mj5Pe/9P4VG64/p8+JoovKvcwy8rO8TuCSNz7v/ef48bsd7y59S3P8/0hOrw4GCr3MBs5e6XfEUTiVrVdO1j6RGdv/r5OfVr3ecHHRPFL5R5mukGHSPnc+Z/J3PvZRG9uf+MoFtc/zsdE8U3lHmYVzChwekNV5GBV2pPPfx+72pt3p1Tk+Mws/wIlCJV7mKnYRQ7eDQve5Z+zR3tzty7D+SK9uX+BEojKPYz0ZqrIwUkpLGDNyA4llqXf9zaYbl0WLrr8QBgNe3up3xFEYt7lKz4vUex3XDWo6EYaKvaw0pZ7GG3eke93BJHYVcqleY/OnEFhim4kHwnacheRiDv7fwtLFPtfL72d9EGzVOwRpC13EYmoZY9fS9X8Xd7c9J7p7KpYycdEyUFb7mFUOy3V7wgiMeOkn1azdkQ7r9hHnXsD6YNmqdijRFvuYdTutPpM1HXcRfhg3G0c9+sP3nzywKlsr1zVx0TJR1vuYZKVncOb83UopCS39E05rB3Rziv2ic0vJ33QLBW7D7TlHiYjZ68kL7/A7xgivnllygO0WpvtzRl3vsov1er4mCi5qdzDRNeUkWRVd/smvhl9ozd/eMwZ3Hzdgz4mElC5h02KQaGuPCBJZuQ7T9BpyYfefH7fcfxQ+wgfE8leKvcwyMrOUbFLUqm5czuLnuzqzUvrHc2VNz3lYyLZl8o9DHQNd0kmmZ9M4I6507z5spueZkW9Jj4mktKo3MNA+9slGVTO38XKx6/15s1VatBiwOs+JpL9UbmHga7hLomu1zcz+NtHf9wRqdMND/PNkSf7mEgOROUeBip2SVQVCgv4TpfmjUs6iSkM6lTVZQck8XRY+u8SxX7zNQ/o0rxxRFvuYaANd0kopVyat8l9M3GmbcF4ou9WGGzN03XcJTFc+N28EsWeefkA0gfNUrHHIW25h0GD2mk6Ykbi3toR7UrMx937FvkVtMsxXumf4zDIbNuUtFTddEDi0+k5y0sU+8MX9CR90CwVe5zTlnsYdGzREIC7pi7Q/neJK1+O7kH97b96c7O7prGjUpqPiSRctOUeRip2iRfH/rKOtSPaecU+LqMD6YNmqdgTSEhb7mZWGxgHnAw4oBewEpgCpANrgc7Ouc2hvE6sy8rOIXPaQr9jiByUNyZmkpGz3Jub93+NLWk1fUwkkRDqlvuTwHvOuROA04DlwGBgjnPuOGBOYE5oQ2cuJV9XDpMY12DbRtaOaOcV+8wTW5E+aJaKPUGVe8vdzGoCrYCeAM653cBuM+sAXBh42ATgY2BQKCFjWVZ2Dlt0KKTEuGdmjKDdis+8+ezbXiK3Zl0fE0mkhbJb5mjgZ+AlMzsNmA8MAA53zuUCOOdyzaxeaU82sz5AH4DGjRuHEMNfQ6Yv8juCSJnq7NhK9tPdvPnrRs3o3O0RHxNJtIRS7hWB04F+zrmvzOxJgtgF45wbC4wFyMjIiMt9GlnZOeTlF/odQ6RUD8x5gd7zZnhzm97P8d1hR/qYSKIplHJfD6x3zn0VmN+gqNw3mFn9wFZ7fWBjqCFjla7jLrGo6u48lo3q5M3ra9bjvNvG+5hI/FDuN1Sdcz8BP5hZ08CiNsAyYCbQI7CsBzCjlKcnhB91VqrEmNvmTitR7B27P6ZiT1KhnsTUD5hkZpWANcBNFP2DMdXMegPrgE77eX5cS0tNYYd2y0gMSC3IZ9WjV5dYlj5olk9pJBaEVO7OuQVARimfahPK140HWdk5KnaJCZ0Xvs8j7/1x/9IenYbxydEtfUwksUCXHygn7W8Xv5kr5PtHriqxTDfSkL10+YFy0v528dMlq+aWKPb+7e/VjTSkBG25l5Mu8yu+KOVGGsdkzqAgRVcllZK05V5OmW2bkqKNJImiM39YUqLYh7bpQ/qgWSp2KZW23EOgy8lItCx8ogu1dv3uzU3vfpNdqZV9TCSxTlvu5aQ3VCUaTty4hrUj2nnF/vTZXUgfNEvFLgekLfdy0huqEmmvTv4r5/9vgTefOmAy26pU9y+QxBWVeznpDVWJlMN+38y8Z7p785RTLmHQFQN8TCTxSOVeTumHqtwl/G7/cir3ffqKN2trXcpL5V4OWdk5fPHdJr9jSAKpuXM7i57s6s0jz+/O6HO6+JhI4p3KvRz0ZqqE0/UL3uOh2c948xl3vMrP1ev4mEgSgcq9HLQ7RsKhSv5OVjx+nTePy+jA8Da3+JhIEonKPUhZ2TkYRXcDFymvK5d/xuiZI7z5vL7jWF/7CB8TSaJRuQdp5OyVKnYpt9SCfL58tieH7dgKwPSTWnN3u3t8TiWJSOUeJB3fLuXVas18Xpn2oDdf0ms0q+oe5WMiSWQq9yDp+HYJVkphAe++1I+mv6wD4OMmLenZaaiu4CgRpcsPBCmzbdMDP0gk4PT1y1kzsoNX7B27P0bPzsNU7BJx2nIXiQTneG3y/ZyzbhEASw4/hvY9RuFM21MSHSr3IGRl5zBk+mK/Y0iMO2Hj97z3Uj9v7tZlOF+kN/cvkCQllXsQRs5eSV5+gd8xJIY9PWME7Vd8BkBu9UM5/9YX2VNBf80k+vRTFwQdKSNlabw5l0/H/nECUt+Of2F203N8TCTJTuUeBB0pI6UZ+sEYen47C4DdKRU5ZeAUXW9dfKd3d4KQ2bapVph46v32K2tHtPOKPfPy/hyfmaVil5igrgpCxxYNKfQ7hMSEAZ+/xtfP9vDmkwdOZdqpl/qYSKQk7ZYRCUKtvN9Y+NT13jy8dS/GnXmNj4lESqdyFzlIPea/zbAPn/fm0/tNYlPVWj4mEimbyj1I1SpV4PfdOhwymVTdnceyUZ28+ZmzO/Noqxt9TCRyYCr3IDU/spbuwpRErl7yEaPeedybz7ltPD/WrOdjIpGDo3I/SEVnpy4iL19vqSaDSnvyyX7qeqrl7wTg9VMvZcjl/X1OJXLwVO4HISs7h8xpC8kv1JXck0Gb1V/x4pt/9+aLbh7DmkMb+ZhIJHgq94MwdOZSFXsSSCks4N8v9OWoLT8B8N7xZ3Pr1ff7nEqkfFTuB2FLXr7fESTCzlq3mCmvD/HmK3s8wdIjjvUxkUhoVO4HkJWd43cEiSTneHNiJi1/XAHANw2b0anbCF1vXeKeyv0ARs5e6XcEiZCTf1rNrAkDvbnL9Q/xVeNT/AskEkYq9wPQhcISkHOMfesfXLpqLgBr6jTg4pufozClgs/BRMIn5HI3swrAPCDHOdfOzA4BpgDpwFqgs3Nuc6iv4wftkkk8x/z6A3PG3ebNva79Gx8de6aPiUQiIxwXDhsALC82DwbmOOeOA+YE5riTlZ1D5hsL/Y4hYfTwu095xf5bpTSOv+ctFbskrJC23M2sEXAl8A/g7sDiDsCFgY8nAB8Dg0J5HT8Me3sp+QU6/DERNNi2kf8818ubB7a7h6yTWvuYSCTyQt0t8wRwH1Cj2LLDnXO5AM65XDMr9VxtM+sD9AFo3LhxiDHCKys7h807dPhjIrjvk5e5fe4b3tzsrmnsqJTmYyKR6Ch3uZtZO2Cjc26+mV0Y7POdc2OBsQAZGRkxtYmsI2Ti3yE7tvLt0928+cGL+zKhZXsfE4lEVyhb7ucCV5nZFUAVoKaZTQQ2mFn9wFZ7fWBjOIJGk+6VGt9u+Wo693883ptP6/86W9Nq7OcZIomn3G+oOueGOOcaOefSga7AR865PwMzgb23qOkBzAg5ZZQ1qK1f2+NRjV2/s3ZEO6/YR517A+mDZqnYJSlF4jZ7DwOXmNkq4JLAHFcy2zYlLVXHPMeTzgvfZ/ETXbz5zNsn8OR5N/iYSMRfYTmJyTn3MUVHxeCc+xVoE46v65eOLRoCMHDKAn+DyAFVzt/F0lGdqOiKLsX8Usv2DLu4r8+pRPynM1Qlbl228gvGZD3kza36vMC6OvV9TCQSO1TuZdARM7GrYsEevhjTi8O3F90Ra+aJreh/1X0+pxKJLSr3MuiImdh03vfZTJz6gDdfdtPTrKjXxMdEIrFJ5V6GKqkpuqVeDDFXyDsvD6DZxu8B+Pyo0/hzl+G6NK9IGVTuZdipYo8ZLXJW8NbEe735mm4j+bbRiT4mEol9KvcyxNQps8nKOV6Z+jdarc0GYHnddK646SmcReIIXpHEonKXmHT8z2t5f/yd3nxjp2F8enRLHxOJxBeVu8ScUW8/ytXLPgbg52q1Ofu2l9lTQT+qIsHQ35hSZGXnkGJQqH0zUdVoy098/vzN3nx7h8H864TzfEwkEr9U7vvIys5hyPTFKvYoe2DOC/SeV3QZokKMZndPY2dqFZ9TicQvlfs+Rs5eSV5+gd8xkkbd7Zv4ZvSN3jy47Z1Mbn6Zj4lEEoPKfR86eSl67vzPZO79bKI3nzJwCr9VruZjIpHEoXLfR+2qqboLU4TV3LmdRU929eYRF/TguT918jGRSOJRuReTlZ3DVhV7RP3523cY/sFz3pxx56v8Uq2Oj4lEEpPKvZihM5ei81IjI233TpaPus6bnz/zGh5q3Ws/zxCRUKjci9mSp632SLhq2Sc89fZIbz731vHk1Cr1vukiEiYqd4mY1IJ8vnmmO7V3bgdg2skXk3nlQH9DiSQJlbtExIXffcPLbwzz5ja9n+O7w470MZFIclG5S1ilFBbwwYt3cMym9QB8eMwZ3Hzt33RpXpEoU7kHZGXn+B0h7p3xwxKmvTbYm9vfOIrF9Y/zMZFI8lK5B+i2eiFwjimvDeas9UsByK7flGu6j9SleUV8pHIP0Jmp5XPShu945+UB3nx913/y5VGn+phIREDl7mlQO40cFXxQnn3rn1zx3/8AsL5mPS7o+wIFKRV8TiUioHL3pB+qcj9Y6Zty+PiFvt7c5+r7ef/4s31MJCL7UrlT9GbqF99t8jtGXPj7+8/SPftfAORVrEzzAa+zq2Iln1OJyL5U7hRddkD274htvzD3uZ7efM8Vd/HmKW38CyQi+6VyR5cdOJC7P32V/l9O8eaTBk7l98pVfUwkIgeicpcy1dmxleynu3nz/110C+PP6OBjIhE5WCp3KVXvb7J44KNx3ty8/2tsSavpYyIRCYbKHaiamsKOfF3sF6Darh0sfaKzNz95TldGnf9nHxOJSHkkfblnZeeo2AOuW/whj/7rCW8+6/aX2VDjMP8CiUi5JX25D5m+yO8Ivqu8ZzeLnuhC5YKiN5ZfbXEFD1x6u8+pRCQUSV3uWdk55CX5Vvul//2SsW/9w5svvOV51h7S0MdEIhIOSV3uyXx8e4XCAj55/mYabfsZgHeanssdHYf4nEpEwqXc5W5mRwKvAEcAhcBY59yTZnYIMAVIB9YCnZ1zm0OPGn7Jenz72f9byOuT7/fmK3o+xbLDj/YxkYiEWyhb7nuAe5xz35pZDWC+mX0A9ATmOOceNrPBwGBgUOhRJVTmCnnr1XtonrsKgLlHnkzX6x/SjTREElC5y905lwvkBj7+zcyWAw2BDsCFgYdNAD4mBss92W7OcWruf5n5yt3efF23EcxrdJKPiUQkksKyz93M0oEWwFfA4YHixzmXa2al3ubezPoAfQAaN24cjhhBSZqjZJxj/BvDuGjNPABWHXokbXs9Q6EuzSuS0EIudzOrDrwJDHTObbOD/BXfOTcWGAuQkZHhQs0RrGQ4SubYX9bx4Yt/HNLY87oH+fiYM3xMJCLRElK5m1kqRcU+yTk3PbB4g5nVD2y11wc2hhpSgvfoO6O4bskcADal1eSsOyaQXyHV51QiEi2hHC1jwIvAcufc48U+NRPoATwc+HNGSAklKI22buDzMb29uV/7TN5udoGPiUTED6FsuZ8LdAcWm9mCwLK/UFTqU82sN7AO6BRSQjloQ/49nr5fT/fmE+96g7xKVXxMJCJ+CeVomc+Bsnawx/RdHBLtSJnDft/MvGe6e/NfL72diS2u8DGRiPgtKc9QHTl7pd8Rwua2udMY9MkEbz51wGS2VanuYyIRiQVJWe6JcCPsmju3s+jJrt488vzujD6ni4+JRCSWJF25/zVrsd8RQnb9gvd4aPYz3nzGHa/yc/U6PiYSkViTdOU+ae46vyOUW5X8nax4/DpvHpfRgeFtbvExkYjEqqQr96ifLRUmVy7/jNEzR3jzeX3Hsb72ET4mEpFYlnTlHm9SC/L58tmeHLZjKwDTT2rN3e3u8TmViMS6pCr3eNvf3mrNfF6Z9qA3X9JrNKvqHuVjIhGJF0lT7lnZOUyMk/3tKYUFvPtSP5r+UpT34yYt6dlpqC7NKyIHLWnKPV7uunT6+uVMn5TpzR27P8aCBk19TCQi8Shpyj3m77rkHJOm3M+5/yu6FPGSw4+hfY9ROEvxOZiIxKOkKfdYdsLG73nvpX7e3K3LcL5Ib+5fIBGJe0lT7ikGhTF4HORTMx/hquWfApBb/VDOv/VF9lRImm+LiERI0rRIrBV74825fDr2jxOQ+nb8C7ObnuNjIhFJJElT7rFk6Adj6PntLAB2p1TklIFT2JVa2edUIpJIVO5RVO+3X/n62R7enHl5f6adeqmPiUQkUSVFucfC9dsHfP4ad33xmjefPHAq2ytX9TGRiCSypCj3u6cs8O21a+X9xsKnrvfm4a17Me7Ma3zLIyLJIeHLPSs7h0KfXrvH/LcZ9uHz3nx6v0lsqlrLpzQikkwSvtyHvR39M1Or7s5j2ag/bh37zNmdebTVjVHPISLJK+HLffOO6J6Z2nHpv3li1mPefM5t4/mxZr2oZhARSehyj+YbqZX25DP/6RuosbvoFn6vn3opQy7vH7XXFxEpLqHLPVo3wm6z+itefPPv3nzRzWNYc2ijqLy2iEhpErrcf4zwjbBTCgv46IVbSd+SC8B7x5/NrVffH9HXFBE5GAld7qkVjN0FkbnuwFnrFjPl9SHefGWPJ1h6xLEReS0RkWAldLlHpNid482JmbT8cQUA3zRsRqduI3QjDRGJKQlb7pG4pd7JP61m1oSB3tzl+of4qvEpYX8dEZFQJWy5h/WWes4x9q1/cOmquQCsqdOAi29+jsKUCuF7DRGRMErIcg/nIZDH/PoDc8bd5s29rv0bHx17Zti+vohIJCRkuYfrWjL/fO9pblg4G4BtlaqS0W8SuyumhuVri4hEUkKWe6jXkqm/7We+fO4mbx7Y7h6yTmod4lcVEYmehCz3UGR+MoE75k7z5mZ3TWNHpTQfE4mIBC/hyj198Dvlet4hO7by7dPdvPnBi/syoWX7cMUSEYmqhCr3s/7xQbmed8tX07n/4/HefFr/19maViNcsUREoi6hyn3Db7uDenyNXb+z+Iku3jzq3Bt48rwbwh1LRCTqEqLcy7MrpvPC93nkvae8+czbJ7CxxqHhjCUi4puIlbuZXQY8CVQAxjnnHg73a5Sn1Cvn72LpqE5UdEXH1LzUsj3DLu4b7mgiIr6KSLmbWQVgNHAJsB74xsxmOueWhes1ylPsl638gjFZD3lzqz4vsK5O/XBFEhGJGZHacj8TWO2cWwNgZpOBDkDYyj0YFQv28PmYXhyxfRMAM09sRf+r7vMjiohIVJhz4b9yopldB1zmnLs5MHcHznLO3VnsMX2APoGxKRDUnTUqHXFsy70fF+zYSoUYvPG0cgVHuYKjXMGJl1y7f1o9P4inH+Wcq1vaJyJV7p2AtvuU+5nOuX5hf7Girz/POZcRia8dCuUKjnIFR7mCk2y5UsL9BQPWA0cWmxsBP0botUREZB+RKvdvgOPMrImZVQK6AjMj9FoiIrKPiLyh6pzbY2Z3ArMpOhRyvHNuaSReK2BsBL92KJQrOMoVHOUKTlLlisg+dxER8VekdsuIiIiPVO4iIgkorsvdzC4zs5VmttrMBkf5tY80s3+b2XIzW2pmAwLLh5pZjpktCPx3RbHnDAlkXWlmbSOYba2ZLQ68/rzAskPM7AMzWxX4s040c5lZ02LrZIGZbTOzgX6sLzMbb2YbzWxJsWVBrx8zaxlYz6vN7CkzswjkGmlmK8xskZm9ZWa1A8vTzSyv2HobE+VcQX/fwp1rP9mmFMu11swWBJZHZZ3tpxui+zPmnIvL/yh6o/Y74GigErAQaBbF168PnB74uAbwX6AZMBS4t5THNwtkrAw0CWSvEKFsa4HD9ln2CDA48PFgYES0c+3zvfsJOMqP9QW0Ak4HloSyfoCvgbMBA94FLo9ArkuBioGPRxTLlV78cft8nWjkCvr7Fu5cZWXb5/OPAX+L5jqj7G6I6s9YPG+5e5c4cM7tBvZe4iAqnHO5zrlvAx//BiwHGu7nKR2Ayc65Xc6574HVFP0/REsHYELg4wlARx9ztQG+c879bz+PiVgu59ynwKZSXu+g14+Z1QdqOue+dEV/C18p9pyw5XLOve+c2xMY51J0zkiZopVrP6K2vg6ULbCV2xl4fX9fI9zZ9tMNUf0Zi+dybwj8UGxez/7LNWLMLB1oAXwVWHRn4Nfo8cV+9YpmXge8b2bzregyDwCHO+dyoeiHD6jnQ669ulLyL5zf6wuCXz8NAx9HKx9AL4q23vZqYmbZZvaJmZ0fWBbNXMF83/xYX+cDG5xzq4oti+o626cbovozFs/lXtq+p6gf12lm1YE3gYHOuW3Ac8AxQHMgl6JfCyG6ec91zp0OXA7cYWat9vPYqK5HKzqp7Spg741qY2F97U9ZOaK93u4H9gCTAotygcbOuRbA3cBrZlYzirmC/b758f28npIbEVFdZ6V0Q5kPLeP1Q8oVz+Xu+yUOzCyVom/eJOfcdADn3AbnXIFzrhB4gT92JUQtr3Pux8CfG4G3Ahk2BH7N2/tr6MZo5wq4HPjWObchkNH39RUQ7PpZT8ldJBHLZ2Y9gHZAt8Cv5wR+hf818PF8ivbTHh+tXOX4vkVtfQGYWUXgGmBKscxRW2eldQNR/hmL53L39RIHgf15LwLLnXOPF1te/ALxVwN738WfCXQ1s8pm1gQ4jqI3S8Kdq5qZ1dj7MUVvyC0JvH6PwMN6ADOimauYEltTfq+vYoJaP4Ffq38zsz8FfhZuLPacsLGim94MAq5yzu0otryuFd03ATM7OpBrTRRzBfV9i1auYi4GVjjnvN0a0VpnZXUD0f4ZK+87wrHwH3AFRe9EfwfcH+XXPo+iX5EWAQsC/10BvAosDiyfCdQv9pz7A1lXEoYjBcrIdTRF77wvBJbuXS/AocAcYFXgz0OimSvwOlWBX4FaxZZFfX1R9I9LLpBP0dZR7/KsHyCDolL7DniGwBnfYc61mqL9sXt/xsYEHntt4Pu7EPgWaB/lXEF/38Kdq6xsgeUvA7fu89iorDPK7oao/ozp8gMiIgkonnfLiIhIGVTuIiIJSOUuIpKAVO4iIglI5S4ikoBU7iIiCUjlLiKSgP4f3nRjzwuCrRAAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(x_test, y_test)\n",
    "plt.plot(x_test, yhat, 'Red')\n",
    "plt.ylim([0, 100])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "meas D = 2.85e-04\n",
      "exp D = 1.88e-04\n",
      "step_nm = 0.04771 nm\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x17f3e942460>]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEjCAYAAAA1ymrVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABKrElEQVR4nO3dd3gU5fbA8e8JvVdBQu+9Q+gIgggogggoNhC9FvR6UUGKBfWionBtP9tFxIr0jqiIgAiICAqh907oNUD6+f0xA3eB3U0C2U07n+fZJ7uz78yenU1Ozs687zuiqhhjjMk8QlI7AGOMMcFlid8YYzIZS/zGGJPJWOI3xphMxhK/McZkMpb4jTEmk7HEn46IyJciMsK930pEtqTANneLSPvrjy79E8cXInJSRFamdjz+pKdYAUTkPhGZn9JtzbWxxJ9Oqepvqlo1teNIS0Skr4gsvY5NtARuAUqpalhyti8ii0XkkSuWtRGR/Vcsu1VElojIWRE5KiK/isgd3rafyD9lv7GmNao6XlU7pHRbc20s8RvzP2WB3ap6LhAbF5EewBTga6AUUBx4GehyDZvzG6uIZL3WOFNaWorFOCzxp2EiUl9E/nKrw0lATo/nLqsmRWSwiBxw224RkXbu8ldEZKqITHKf+0tE6vp4vTAR+V1ETolIhIh8KCLZPZ6vKSI/i8gJETksIsPc5SEiMkREdojIcRGZLCKF3efKiYiKyEMiss89NPG4iDQWkXD3tT68Io5+IrLJbfuTiJT1eE7d9be5z3/kHvaoDnwKNBORSBE55eM9horIbPc9bBeRf7jLHwbGeqz/avI+Lf9ERIB3gH+r6lhVPa2qCar6q6r+I5nbuirWi78P7u/BIeALEckhIu+JyEH39p6I5HC3cbH98yJyxP28u4lIZxHZ6u6fYX5iKCAiX7vfWvaIyIsiEuI+11dElonIuyJyAnjFy7eZDu7v6WkR+dj95vOIx/qebb1+5snZZ+ZylvjTKDfhzgS+AQrjVIp3+WhbFXgKaKyq+YBbgd0eTbq66xcGvgNmikg2L5uKB54BigLNgHZAf/c18gELgB+BUKAS8Iu73tNAN+Am97mTwEdXbLsJUBm4G3gPeAFoD9QEeonITe7rdAOGAd2BG4DfgAlXbOt2oDFQF+gF3Kqqm4DHgd9VNa+qFvS2r9xt7Xfj7AG8ISLtVPXzK9Yf7mP9a1UVKA1Mvd4N+Yn1RpzPuCzwKM4+bgrUw9lXYcCLHpu6EaeYKInzzeMz4H6gIdAKeFlEKvgI4/+AAkAFnM/9QeAhj+ebADuBYsDrniuKSFGc/TAUKAJsAZon8rav+swTaW/8UVW7pcEb0Bo4CIjHsuXACPd+G2C/e78ScAQnkWa7YjuvACs8HocAEUAr9/FuoL2PGAYAM9z7vYG/fbTbBLTzeFwCiAWyAuUABUp6PH8cuNvj8TRggHv/B+DhK+I9D5R1HyvQ0uP5ycAQ935fYKmffVoa559bPo9lbwJfJnF9n88Di904T3ncIj0+oxZu7DmTuv1EPpsr27YBYjy3D+wAOns8vhXn8NDF9heALO7jfG58TTzarwa6eXntLEA0UMNj2WPAYo/Y9vqKF+efxO8ezwmwD3jEx3vz+Znb7dpuVvGnXaHAAXV/0117vDVU1e04SfoV4IiITBSRUI8m+zzaJvC/ivcyIlJFROaKyCEROQO8gVP9g5M0d/iItSwwwz1scwrnH0E8zjHsiw573L/g5XFej22977GtEziJoaRH+0Me9897rJuYUOCEqp71WLbnim1fj6dVteDFG06VetFx92eJFHotb46qapTH41Au/53Zw+Wf+3FVjXfvX3B/+vpcPBUFsnvZtud+3IdvoVz+O6k4v5P+XOtnbrywxJ92RQAlrziWWcZXY1X9TlVb4iROBd7yeLr0xTvucdhSON8mrvQJsBmorKr5cQ65XHz9fUBFHy+/D+jkmfRUNaeqHvD7Dn1v67ErtpVLVZcnYd3Eppo9CBR2D1tdVAa4ljiTawvOe/N6uC6FXPn+D+L8PlxUBu+fe3Idw/lGd+W2Pfejv88iAud3ELh0/qOU7+YmpVniT7t+B+KAp0Ukq4h0xzlGexURqSoiN7sn7qJwKrV4jyYNRaS7OL0rBuB8TV/hZVP5gDNApIhUA57weG4ucKOIDHBPGuYTkSbuc58Cr188CSsiN4hI12t8358CQ0WkprutAiLSM4nrHgZKiccJaU+qug/ncNmbIpJTROoADwPjkxGfuOteuiVlJbeqfRZ4SZwT3fnFOSneUkTGJOP1k2MC8KL7eRTFOY7/7fVu1P2WMBnnM8/nfu7PJmPb3wO13ZPJWYEncc43mCCxxJ9GqWoMzgnOvjgnS+8GpvtongMYiVOJHcI5oebZI2OWu/5J4AGgu6rGetnOQOBe4CzOib5JHvGcxek33sV9jW1AW/fp94HZwHwROYvzT6UJ10BVZ+B8W5noHm5aD3RK4uoLgQ3AIRE55qNNb5zzDgeBGcBwVf05GSE2x/nHeukmSeyuqKpTcT6Hfu7rHwZG4Hw+gTACWAWEA+uAv9xlKeGfwDmcE7hLcToNjEvKiqp6DOgJvI1zCKyGG2d0CsVmEiGXH0I2GY2IvAJUUtX7UzsWY7xxDz/uB+5T1UWpHU9mYBW/MSboxBnBXNA9PHnxXJK3w48mACzxG2NSQzOcXmLHcA4fdlPVC/5XMSnFDvUYY0wmYxW/SVNE5AcR6ZPacXgSkYYiMvyKsRHpVlrcxya4rOI3ASMiu3EGccXhdC/diDNB2Rh3IFma53Zr/Q1nfEMBoLWqnvJ4vg/OlBWVcbrCfgcMU9U4jzahwEpVTXN91UWkL86I2ZapHYsJHqv4TaB1UWf+oLI4XU4HA5+nbkhJIyKlcKaQGIozhcYiYPYVffdz44yNKIrThbUdTrdYT51x5jgyJk2wxG+CQp3ZKGfj9GPvIyK1vLUTj3ntL87SKCKj3VkZd4mIzz794sxfP1CcWT9PizMjaU73uWTNRinO7KLfAy+pM5umquq/cAaATRSRLO77+kSdayPEuCOVx+PMy+OpMzAvCTEWdafMOOXG85vb1fHK93lxxtOsHsuSvN8uthUfM5q6+2SjOLO5HhCRK/+RmXTOEr8JKlVdidNnu1USV2mCM91BUZwBP5+7Q/x96QV0BMoDdXAGwF2U5NkoVfWEqtZV1ctGo6rqEFXt5jHHzZVa4wwiA0CcWVBbA56DxHzF+BzOvrkB5xDZMBKfhsKXRPeb+p7R9HOcaTPyAbVwBsaZDMQSv0kNB3GmD06KPar6mZtov8KZ5Ky4n/YfqOpBVT0BzMGZkviiWOB1d9TyRJyk+L6qnlXVDTgJu07y3sr/iMhDQCNgtMfi1sDaKyaG8xVjrPv+yqpqrPtN4loTf3L3m6dYoIaI5FfVk6r61zXGYNIoS/wmNZTEmXUzKS7Nyqiq5927/mZm9DeL47XORpkoca4jMBJnsjrP6SIuHeZJQoyjgO04U1/sFJEh1xLLla+RxP3m6S6cuPeIc4GUZtcRh0mDLPGboBKRxjiJ/3qujZumiEhHnMNGXVR13RVPd8Y5V5Ao95vHc6paAWdQ07PiXkntChcvt5jbY9m1TnJ21TcKVf1TVbvizPk0E2dCNpOBWOI3QeHORnk7ziGWb70kyHRJRG7GOaF7l3v+wvO58kAOVd2cxG3dLiKV3GPxZ3C6wF51LkFVj+JMgXy/iGQRkX74njI7MZfNaCoi2UXkPhEp4B4SuxiHyUAs8ZtAm+PO2LkP51KA73D5JfrSu5dw+vfPc3vGRIrID+5zt3H1YR5/KuNc3jISZ1ruj1V1sY+2/wAG4cxuWROnt9G18Daj6QPAbnd21MdxToCbDMQGcBkTICIyD/hQVZOT/I0JOKv4jQmcxTiDvoxJU6ziN8aYTMYqfmOMyWSSdMm41Fa0aFEtV65caodhjDHpyurVq4+p6g1XLk8Xib9cuXKsWrUqtcMwxph0RUT2eFtuh3qMMSaTscRvjDFpUXwsXDgVkE1b4jfGmLTm0Dr4rC3MfAIC0PPSEr8xxqQV8bGw+C0Y0wbOHoZ694HfWcivTbo4uWuMMRneofVOhX8oHGr3hE5vQ+6kzl6ePJb4jTEmNcXHwtJ34de3IVdBuHs8VL89oC9pid8YY1LL4Q1OlR+xFmr1gM6jAlble7LEb4wxwRYfC0vfg1/fcqr8Xt9AjTsua3LkTBQFcmcjR9YsKf7ydnLXGGOC6fAGGNsOFo1wkn3/P65K+qv3nOC2/1vKm/OSdCmHZLOK3xhjgiE+Dpa96/TayVnAa5Wvqny7Yg+vzd1IaMFc9A4rE5BQLPEbY0ygHd7oHstfAzW7Q+fRkKfIZU3Ox8QxbPo6Zq45yM3VivFur3oUyJ0tIOFY4jfGmEC5qsr/Gmp0varZjqORPPHtarYdiWRghyr0b1OJkJCU779/kR3jT4Y2bdpQtWpV6tSpQ7Vq1Xjqqac4derUdW93/Pjx1KlThzp16tC8eXPWrl3rtd0vv/xCgwYNqFevHi1btmT79u0ALF68mAIFClCvXj3q1avHa6+9luwY3nzzTSpVqkTVqlX56aefrnr+jjvuoFatWsnerjGZ1pFN8Hl7WDjC6Z755B9ek/4P6yLo+uEyjkXG8HW/MJ66uXJAkz7gHFNK67eGDRtqWnDTTTfpn3/+qaqq0dHR+uyzz2rr1q2ve7vLli3TEydOqKrqvHnzNCwszGu7ypUr68aNG1VV9aOPPtI+ffqoquqiRYv0tttuu+bX37Bhg9apU0ejoqJ0586dWqFCBY2Li7v0/LRp07R3795as2bNa34NYzKNuFjVJaNVXyuq+lYF1fUzvDaLjo3XV2dv0LKD5+odHy7VAyfPp3gowCr1klPTXcW/e/duqlWrxiOPPEKtWrW47777WLBgAS1atKBy5cqsXLkSgHPnztGvXz8aN25M/fr1mTVr1qX1W7VqRYMGDWjQoAHLlzvXqI6IiKB169bUq1ePWrVq8dtvv/mNI3v27Lz99tvs3bvXZ4WeVM2bN6dQoUIANG3alP3793ttJyKcOXMGgNOnTxMaGprotr/99lvCwsKoV68ejz32GPHx8Ve1mTVrFvfccw85cuSgfPnyVKpU6dJ+jIyM5J133uHFF1+81rdnTOZxscr/5TWo2tmp8mt2u6rZvhPn6fnpcsYt20Xf5uWY/FhTQgvmClqY6fIY//bt25kyZQpjxoyhcePGfPfddyxdupTZs2fzxhtvMHPmTF5//XVuvvlmxo0bx6lTpwgLC6N9+/YUK1aMn3/+mZw5c7Jt2zZ69+7NqlWr+O6777j11lt54YUXiI+P5/z584nGkSVLFurWrcvmzZupW7fuZc/dfffdbNmy5ap1nn32WR588EGf2/z888/p1KmT1+fGjh1L586dyZUrF/nz52fFihWXnvv999+pW7cuoaGhjB49mpo1a7Jp0yYmTZrEsmXLyJYtG/3792f8+PFXvf6BAwdo2rTppcelSpXiwIEDALz00ks899xz5M6dO9H9YUymFR8Hy9+HxSMhRz7o+SXUvNNr0x/WRfD8tHAAPr2/AR1rlQhioI50mfjLly9P7dq1AahZsybt2rVDRKhduza7d+8GYP78+cyePZvRo0cDEBUVxd69ewkNDeWpp55izZo1ZMmSha1btwLQuHFj+vXrR2xsLN26daNevXpJikV9zJw3adKkZL+vRYsW8fnnn7N06VKvz7/77rvMmzePJk2aMGrUKJ599lnGjh1LgwYN2LNnD3nz5mXevHl069aNbdu28csvv7B69WoaN24MwIULFyhWrFiS3oOIsGbNGrZv38677757ab8aY65wZLPTY+fgX1CjG9z2H8hT9KpmUbHxvDFvE1//voe6pQrw4b0NKF04dQqqdJn4c+TIcel+SEjIpcchISHExcUBTjKbNm0aVatWvWzdV155heLFi7N27VoSEhLImTMnAK1bt2bJkiV8//33PPDAAwwaNMhvZQ4QHx/PunXrqF69+lXPJbfiDw8P55FHHuGHH36gSJEiVz1/9OhR1q5dS5MmTS5tv2PHjgDkz5//UrvOnTvTv39/jh07hqrSp08f3nzzzcu2NWPGDF599VXA+RZRqlQp9u3bd+n5/fv3Exoayu+//87q1aspV64ccXFxHDlyhDZt2rB48WK/+8WYTCE+DpZ/AIvfdKr8Hl9Are5em+48Gsk/J/zNhoNneLhleQZ3rEb2rKl4pN3bgf+0dvM8ubtr167LTjL26dNHp0yZctVzQ4cO1SeffFITEhJUVfWvv/5SVdUBAwbo6NGjVVV13Lhx6uwC1d27d2tsbKyqqr777rv6r3/966oTJZ4nd2NiYnTQoEEpcnJ3z549WrFiRV22bJnPNrGxsVqkSBHdsmWLqqqOHTtWu3fvrqqqERERl97nH3/8oaVLl9aEhATdsGGDVqpUSQ8fPqyqqsePH9fdu3dfte3169dfdnK3fPnyl53cVb16vxuTqR3epDqmrerw/KoT71c9e8Rrs4SEBJ2yap9Wf+kHrfvqTzp/w6GghomPk7vpsuJPipdeeokBAwZQp04dVJVy5coxd+5c+vfvz1133cWUKVNo27YtefLkAZwukaNGjSJbtmzkzZuXr7/+2ut277vvPnLkyEF0dDTt27e/dNL4erz22mscP36c/v37A5A1a9ZL1xju3LkzY8eOJTQ0lM8++4y77rqLkJAQChUqxLhx4wCYOnUqn3zyCVmzZiVXrlxMnDgREaFGjRqMGDGCDh06kJCQQLZs2fjoo48oW7bsZa9fs2ZNevXqRY0aNciaNSsfffQRWbKk/PwgxqR78XHw+//Bojcge16/Vf7ZqFhemrmemWsO0qR8Yd6/pz43FsgZ5IC9Ew3A1V1SWqNGjdQutm6MSVVHtzjH8g+shup3wG3vQN4bvDZdu+8UT0/8m30nzjOgfRWebFuJLIHum++FiKxW1UZXLs+wFb8xxqSIhHhYfrHKzwM9xjnTLni5MlZCgjJ26U7e/nELxfLlYNJjzWhcLvDTLCdXuuvHHyiBGpW7efNmmjVrRo4cOS71MPLG16jc06dP06VLF+rWrUvNmjX54osvkh2Djco15hod3Qqfd4AFw6FKB6dffq27vCb9PcfPce/YFbwxbzPtqhdj3r9apcmkD5b4LzN+/HjCw8MJDw8nR44cdO169fDq5CpcuDAffPABAwcO9NvuiSeeYPz48axZs4Z7772XESNGAPDRRx9Ro0YN1q5dy+LFi3nuueeIiYlJ8utv3LiRiRMnsmHDBn788Uf69+9/2SCu6dOnkzdv3mt7c8ZkVAnxsOx9+LQlnNgJd33uzKaZ9+ru0AkJyrilu+j43m9sOHCGkd1r8+n9DSmYO3sqBJ406TLxexuN+ueff1KnTh2ioqI4d+4cNWvWZP369SxevJjWrVtz5513UqNGDR5//HESEhL8bj8lR+UWK1aMxo0bky2b/1n2fI3KFRHOnj2LqhIZGUnhwoXJmjWrz/1wJRuVa0wyHd0K426Fn1+Gyrc4VX7tHl6r/B1HI+n13995be5GmlYozPxnW3NPWBkkABdIT0np7hi/v9God9xxBy+++CIXLlzg/vvvp1atWixevJiVK1eyceNGypYtS8eOHZk+fTo9evTw+zqBGJXrj69RuU899RR33HEHoaGhnD17lkmTJhESEmKjco1JaQnx8PuHsPB1yJ7bqfJ9HNaJT1DG/raTd37eSo6sIfynZ126NyiZ5hP+Reku8fsbjfryyy/TuHFjcubMyQcffHBpnbCwMCpUqABA7969Wbp0aaKJH1J2VG5ifI3K/emnn6hXrx4LFy5kx44d3HLLLbRq1cpG5RqTko5tg5n9Yf9KqHa702MnX3GvTbcePsugqeGs3XeKDjWKM6JbLYrlTxvdNJMq3SV+9TEaFeDEiRNERkYSGxtLVFTUpT76V/4XTsp/5ZQclZsYf6Nyv/jiC4YMGYKIUKlSJcqXL8/mzZttVK4xKSEhHlZ87EydnC0XdB/r87BObHwC//11Bx/8sp08ObLwQe/6dKlTIt1U+ZfxNqorrd08R+76G43apUsXHT9+vI4YMUKffPJJVXWmLM6ZM6fu3LlT4+PjtUOHDjp16tSrRrgFalTuRcOHD9dRo0Z5fc7fqNzHH39chw8frqqqhw4d0tDQUD169KiNyjXmeh3dqvpZe2f07Xe9Vc/4HlW74cBp7fz+Ei07eK72H79aj56NCmKg146MMnLX12jUX3/9laxZs3LvvfcSHx9P8+bNWbhwISEhITRr1owhQ4awbt26Syd6vQnEqNxDhw7RqFEjzpw5Q0hICO+99x4bN24kf/78SRqV+9JLL9G3b19q166NqvLWW29RtGhRihYtaqNyjbkWnlV+1pzQ/TOo3dNrlR8Tl8CHi7bz8aLtFMydPdVm00xpGX7k7uLFixk9ejRz585N4aiMMenOse0wqz/s+wOq3ga3v+vzWH74/lMMmhLOlsNnubN+SV6+vQaF8qTdLpre2MhdY0zmlRAPKz6Bhf92qvw7x0CdXl6r/KjYeN7/ZRtjluykaN7sfN6nEe2qe//nkF5l+MTfpk0b2rRpk9phGGNSy7HtMOtJ2LcCqnSCLu9Bvhu9Nl295yTPT13LjqPnuLtRaYbdVp0CufyPwUmP0uUALm/ee++9JF01K71YvXo1tWvXplKlSjz99NM+u5YmdzqGZ5555tJF2atUqULBggUB5yIwF5fXq1ePnDlzMnPmTAA+/PBDKlWqhIhw7NixFH+vxgREQjz8/hF82gKObnaq/N4TvCb9CzHxjJi7kR6fLicqNoGv+4XxVo86GTLpA4Ht1QM8A2wA1gMTgJxAYeBnYJv7s1Bi20nKxdbLli2rR48evY7z32lL48aNdfny5ZqQkKAdO3bUefPmXdXmei+S/sEHH+hDDz101fLjx49roUKF9Ny5c6rqXMtg165dGW4fmwzs2HbVsR2cHjvj71Y9E+Gz6Yodx/Smtxdq2cFz9YUZ4XrmQkwQAw0sgn2xdREpCTwNNFLVWkAW4B5gCPCLqlYGfnEfJ9m5c+e47bbbqFu3LrVq1WLSpEl88MEHHDx4kLZt29K2bVvAufRis2bNaNCgAT179iQyMhKAcuXKMXjwYMLCwggLC7s0GZo3ffv25emnn6Z58+ZUqFCBqVOnAs4J45tuuolevXpRpUoVhgwZwvjx4wkLC6N27drs2LEjubvrMhEREZw5c4ZmzZohIjz44IOXqm9P1zsdw4QJE+jdu/dVy6dOnUqnTp0ujeitX78+5cqV8xvzl19+Sffu3enYsSOVK1fm+eefv/Rc3rx5GTx4MA0bNqR9+/asXLmSNm3aUKFCBWbPnp2EPWJMEiUkwO8fwyct4OgmuPO/Pqv8c9FxDJ+1nrvHrCBB4bt/NGFEt9rky5lBq3xP3v4bpMQNKAnsw6nwswJzgQ7AFqCE26YEsCWxbXlW/FOnTtVHHnnk0uNTp06p6uUV/9GjR7VVq1YaGRmpqqojR47UV1999VK7ESNGqKrqV199pbfddpvP/5Z9+vTRHj16aHx8vG7YsEErVqyoqs7YgAIFCujBgwc1KipKQ0ND9eWXX1ZV1ffee8/r1bsWLlyodevWverWrFmzq9r++eef2q5du0uPlyxZ4jXOJ598Ur/55ptLj/v163fpamQDBgzQ6dOn++yjv3v3br3xxhuv6tOvqtq2bVudM2fOVcv9VfxffPGFli9fXk+dOqUXLlzQMmXK6N69e1VVFbj0jaVbt256yy23aExMjK5Zs0br1q3rdXvGJNux7aqf3+pW+b1UTx/02XTFjmPa6q2FWm7IXB0+a72ei44NYqDBQ7D78avqAREZDewFLgDzVXW+iBRX1Qi3TYSIXD3PACAijwKPApQpU+bS8tq1azNw4EAGDx7M7bffTqtWra5ad8WKFWzcuJEWLVoAEBMTQ7NmzS49f7HK7d27N88884zf99GtWzdCQkKoUaMGhw8fvrS8cePGlCjh9OetWLEiHTp0uBTfokWLrtpO27ZtWbNmjd/Xukh9TLOQ1HZJmY5h4sSJ9OjR46o+/REREaxbt45bb701SbF6ateuHQUKFACc8RZ79uyhdOnSZM+e/dJI5Nq1a5MjRw6yZctG7dq1bboIc/0SEmDlf2HBq5A1O3T7FOre47XHzoWYeN7+aTNfLNtN2SK5mfRoM8LKp82pkwMpYIlfRAoBXYHywClgiojcn9T1VXUMMAacfvwXl1epUoXVq1czb948hg4dSocOHXj55ZevXJdbbrmFCRMm+IrN631vPC/s7plok3LBd0+LFi3y+k8md+7cLF++/LJlpUqVYv/+/ZceX5xm4UrXMx3DxIkT+eijj67a5uTJk7nzzjsTnU3UG899kiVLlkv7IVu2bJf2c1L2lTFJdnwHzHoK9i6HyrdCl/chv/cBVqt2n2DglLXsPn6evs3L8XzHquTOnuE7NnoVyF497YFdqnpUVWOB6UBz4LCIlABwfx5JzkYPHjxI7ty5uf/++xk4cCB//fUXAPny5ePs2bMANG3alGXLll06fn/+/Hm2bt16aRsXJ1mbNGnSZd8EAulixX/l7cqkD1CiRAny5cvHihUrUFW+/vprr9cGuOOOO5g4cSLR0dHs2rWLbdu2ERYWxhNPPMHBgwfZvXs3S5cupUqVKpcl/S1btnDy5Emv793XcX9j0pSEBKdf/ict4PAG6PYJ3DvJa9KPinV67PT87+/EJSgT/tGUV+6omWmTPgS2H/9eoKmI5MY51NMOWAWcA/oAI92fyZoXYd26dQwaNIiQkBCyZcvGJ598AsCjjz5Kp06dKFGiBIsWLeLLL7+kd+/eREdHAzBixAiqVKkCQHR0NE2aNCEhIcHnt4LU9sknn9C3b18uXLhAp06d6NSpEwCzZ89m1apVvPbaa9c8HcOECRO45557rvq2s3v3bvbt28dNN9102fIPPviAt99+m0OHDlGnTp1LU00Ykyouq/I7uFX+1d+IAf7ae5KBk9ey89g57m9ahqGdqpMnR+ZN+BcFdMoGEXkVuBuIA/4GHgHyApOBMjj/HHqq6gl/20nJi62XK1eOVatWUbRo0RTZnjEmSBISYOUYWPAKZMkOHd+Eevf6HH377oKtfLZkJyUK5OLtHnVoUSnz/c2nypQNqjocGH7F4mic6t8YY5LmxE6nyt+zDCrdAnd84LPKX7vvFM9NWcv2I5H0DivDsM7VMkcXzWTIdN95vPUief3115kyZcply3r27MkLL7wQpKiMMV4lJMCfnzlVfkg26PoR1LvPa5UfHRfPB79s49Nfd1IsXw6+6hfGTVVuCH7M6UCGn53TGJNOndjlVvlLnSq/y/tQoKTXpuv2n2bglLVsOXyWXo1K8eLtNchvVb7NzmmMSScSEuDPsbBgOIRkhTs+hPr3+54vf+E2Plq8g6J5s/NF38a0reZ1aJDxYInfGJN2nNgFs/8Ju3+DSu3dKr+U16YbDp5m4JRwNkWcoXuDkgy/vSYFcluVnxSW+I0xqS8hAVZ9Dj8Ph5Asfqv82PgEPl60g/9buI1CebLz2YONuKVGxpovP9As8RtjUtfJ3c6x/N2/QcV2To8dH1X+5kNneG7yWjYcPEPXeqG80qVmursqVlpgid8YkzquqvL/D+o/4LXKj4tP4L9LdvLegq0UyJWNT+9vSMda3i+mYhJnid8YE3yXVfk3O0nfR5W/7fBZnpuylvD9p7m9Tgle61qLwlblXxdL/MaY4ElIgNXjYP7LICHQ5QNo8KDPKv+z33bx7s9byZszKx/d24Db6nifgM0kjyV+Y0xwnNzjXPt2929Qoa1T5Rcs7bXp9iORDJyyljX7TtGp1o38u1stiubN4bWtST5L/MaYwLqqyn8fGvTxWuUnJCjjlu1i1E9byJU9Cx/0rk+XOiUSnT7dJI8lfmNM4JzcA7Ofgl1LEq3y9x4/z8Cpa1m56wTtqxfnje61KJYvZ5ADzhws8RtjUp4qrBoHP7sXSbr9PWjY12uVr6qM/2Mvb8zbRBYRRvesy10NSlqVH0CW+I0xKevUXqfHzq5foUIbt8ov47VpxOkLPD81nN+2HaNlpaK83aMOoQVzBTfeTMgSvzEmZajC6i9g/kvO49vfhYYP+azyZ/x9gOGzNxAXr/y7a03ub1rWqvwgscRvjLl+p/Y6c+zsXAzlb4KuH/qs8o9FRjNs+jrmbzxMo7KFGN2zLuWK5gluvJmcJX5jzLVThdVfulW+wm3vQKN+Xqt8gB/XRzBsxnoio+IY2qkaj7SqQJYQq/KDzRK/MebanNrnVvmLoHxrZ2K1QmW9Nj15LoZX5mxg1pqD1CqZn3d61aNK8XxBDthcZInfGJM8qvDXV/DTi6AJiVb58zccYtiM9Zw6H8O/2lXmqZsrkS1LSJCDNp4s8Rtjku7UPpjzNOxYmKwqv3qJ/HzVrzE1QwsEOWDjjSV+Y0zirqry/wMN+0HI1ZW7qjInPILX5mzg1PlYBrSvTP82lcie1ar8tMISvzHGv9P7YfbTsOMXKNfK6bFTqJzXpgdOXeDFGetYtOUodUsV4Ot+TagRmj+48ZpEWeI3xninCn9/Az+9AAnx0Hk0NHrYa5Ufn6B8tXw3o+dvAeCl22vQt3k567GTRlniN8ZcLRlV/saDZxg6PZy1+0/TtuoN/LtbLUoVyh3ceE2yWOI3xvzPZVV+nN8qPyo2nvd/2cZnS3ZSIFc23r+nHnfUDbXRt+mAJX5jjOP0AafHzvYFULalU+UXLu+16e87jjN0eji7j5+nR8NSvNC5ul37Nh2xxG9MZqcKf38LPw1zqvxOo6DxI16r/HPRcYz8YTPfrNhD2SK5Gf9IE1pUKpoKQZvrYYnfmMzs9AGY8y/Y/jOUbeFW+RW8Nl2+/RjPTwvnwKkLPNyyPAM7VCVX9ixBDtikBEv8xmRGqrBmPPw4DBJiodPb0PgfXqv8yOg43py3ifF/7KV80TxMeawZjcoVToWgTUqxxG9MZnPmoFPlb5ufaJW/dNsxBk8L5+DpC/yjVXme61CVnNmsyk/vLPEbk1mowprv4MehEB8DHd+CsEe9Vvlno2J5Y95mJqzcS4WieZj6eDMalrUqP6OwxG9MZuBZ5Zdp7lT5RSp6bfrbtqMMmbaOiNMXeLR1BZ69pYpV+RmMJX5jMjJVWDsBfhiSxCp/ExNW7qPCDXmY8nhzGpYtlApBm0CzxG9MRnUmwq3yf0q0yv9j53GenbyWiNMXeOymCjzT3qr8jMwSvzEZzcUq/8chEBcDHUdC2GNeq/yYuATeW7CVT37dQdnCuZn6RHMalLEqP6MLaOIXkYLAWKAWoEA/YAswCSgH7AZ6qerJQMZhTKZxJgLmDoCtP0KZZtD1I59V/o6jkQyYuIZ1B05zd6PSvNylBnlyWC2YGQT6U34f+FFVe4hIdiA3MAz4RVVHisgQYAgwOMBxGJOxqUL4JPjheafKv/VNaPIYhFx9uEZVmfjnPl6bs5Ec2UL49P4GdKxVIhWCNqklYIlfRPIDrYG+AKoaA8SISFegjdvsK2AxlviNuXaeVX7pptDtY59V/olzMQyeFs7PGw/TslJRRvesy40FcgY3XpPqAlnxVwCOAl+ISF1gNfAvoLiqRgCoaoSIFAtgDMZkXJdV+dFw6xvQ5HGvVT7Akq1HeW7KWk6fj+XF26rTr0V5Qmy+/EwpkIk/K9AA+Keq/iEi7+Mc1kkSEXkUeBSgTJkygYnQmPTq7CGY+wxsmQelm0DXj6FoJa9No2LjefvHLYxbtosqxfPy1UNhdlWsTC6QiX8/sF9V/3AfT8VJ/IdFpIRb7ZcAjnhbWVXHAGMAGjVqpAGM05j0QxXCJ7tVfhR0eB2aPuGzyt9+JJJ/TvibTRFn6Nu8HEM6VbNumiZwiV9VD4nIPhGpqqpbgHbARvfWBxjp/pwVqBiMyVA8q/xSYc6x/KKVvTZVVSav2scrszeSO3sWvujbmLbV7KiqcQS6V88/gfFuj56dwENACDBZRB4G9gI9AxyDMembKqybAvMGuVX+CGja32eVfyYqlmHT1zE3PIIWlYrwbq96FMtvJ3DN/ySa+N3eOTeo6o4rltdR1XB/66rqGqCRl6faJSdIYzKts4fdKv/7RKt8gL/2nuTpCX8TcTqKQbdW5YmbKtoJXHMVv4lfRHoB7wFHRCQb0FdV/3Sf/hLn5K0xJqWpwrqp8MMgiL2QaJWfkKB8umQH/5m/lRvz52TyY81snh3jU2IV/zCgoXsiNgz4RkSGqep0wMoIYwLh7GH4/lnYPBdKNYZun/it8o+cieKZyWtYtv04t9UpwRt31qZArmxBDNikN4kl/iwefe5XikhbYK6IlMKZgsEYk1JUYf00mDcQYs7DLf+GZk/6rPIBFm05wsDJazkXE8fI7rW5u3FpRKwmM/4llvjPikjFi8f33cq/DTATqBnY0IzJRCKPOMfyL1b5XT+GG6r4bB4VG8+on7bw+dJdVLsxH5PubUqlYvmCGLBJzxJL/E9wxSEdVT0rIh2BXgGLypjM4qoq/zVo9pTfKn/d/tM8O3kN245E8mCzsgzrXN365ptk8Zv4VXWtj+WxwPiARGRMZhF5xDmWv2kOlGzk9Ni5oarP5rHxCXy8aAf/t3AbRfJm58uHGtOmqvXNN8mXpH78ItJXVb8McCzGZA6qsGE6fD8QYs5B+1edKj+L7z/H7UcieW7yGtbuP80ddUN5rWtNCubOHsSgTUaSlH78LwNVcLpvGmOuR+RR+P4Zt8pv6PTY8VPlJyQoXy7fzVs/biZX9ix8eG99bq8TGsSATUaUWD/+MUA+4L7ghGNMBnVZlR8J7V+BZv/0W+XvP3meQVPC+X3ncW6uVoyR3WvbCFyTIhKr+O8FwlQ1IRjBGJMhRR51j+XPhtAGTpVfrJrP5qrKlNX7eW3ORlTVummaFJdY4r8dZ16drldO2WCMSYL1050eO9Fnod1waP603yr/6Nlohk5fx4JNhwkrV5j/9KpL6cK5gxiwyQwS69WzWETuAb4FmgUnJGMygMijMO852DgLQuu7VX51n81VlTnhEbwyewOR0XG80Lk6/VqWJ4vNs2MCINGTu6q6XkTuCkYwxmQIG2bA98+5Vf7L0Pxffqv8Q6ejeHHmehZsOkydUgUY3bMuVYrbYCwTOEnqzqmqBwMdiDHp3rljTsLfODPJVf6kP/fx+rxNxMQlMKxzNfq1KE/WLCHBi9lkSsmaj19EBgDfquqxwIRjTDq1YaZb5Z9JUpW/78R5hkwPZ9n244SVL8xbd9WhfNE8wYvXZGpJTvwiUgfnqllZgdEBi8iY9OTcMefk7YYZUKKeU+UXr+GzeXyC8tXy3Yz6aQtZQoQR3Wpxb1gZmzPfBFVyKv6HgeeBR7DEb4xz4nbusxB1Gm5+CVoM8Fvl7zwayaCp4azec5I2VW/gjTtrE1owV/DiNcaV1CkbcgCdgVpAGxFpqapLAxqZMWnVueNulT/dqfL7zPFb5SckKF/97oy+zZE1C+/0qsud9Utav3yTapJa8d8F/KSq0SLyBU71b4nfZD6XVfkvulW+74ueeI6+bVP1Bt66qw7FbfStSWVJTfwPAwPd+/OA90Qkr6pGBiYsY9KYc8edyyCunwYl6kKf2VDc9yUpVJXJq/bx77mbbPStSXOSMklbQeCgqv4NoKrxIvIhEAYsDGx4xqQBG2c7Uy5cOAVtX4SWA/xW+UfORDFk+joWbj5C0wqFGdXDRt+atCUpA7hOAQ9csezdQAVkTJpx/oRzLH/9NLixDjwwE26s5XeV2WsP8tLM9UTFxjO8Sw36NCtnPXZMmpPUk7tZgNuAcp7rqOo7gQnLmFS2aY5zKcQLp6DtC9DyGb9V/olzMbw0az3fh0dQr3RB/tOrLhVvyBu8eI1JhqQe458DRAHrAJup02Rc50/AvEGwfmqSq/wFGw8zZPo6Tl+IYdCtVXmsdQUbfWvStKQm/lKqWiegkRiT2jbNdav8E9BmGLR61m+VfyYqltfmbGTq6v1UL5Gfbx4Oo3qJ/EEM2Jhrk9TE/4OIdFDV+QGNxpjUcP4E/PA8rJsCN9aGB6Y7P/1Ytv0Yg6as5dCZKJ5qW4mn21Ume1ar8k36kNTEvwKYISIhQCwggKqqlTcmfdv8PcwZ4Fb5Q6HVc36r/PMxcYz8YTNf/76HCjfkYdoTzalfplDw4jUmBSQ18f8HZz7+daqqAYzHmOA4fwJ+GAzrJie5yl+95wTPTV7LnhPnebhleQbdWpWc2bIEKWBjUk5SE/82YL0lfZMhbJ4HcwfA+eNOld/yWcia3WfzqNh43l2wlc+W7CS0YC4m/KMpTSsUCV68xqSwpCb+CGCxiPwARF9caN05Tbpy/gT8OATCJ0Hx2nDfVCjhv8/Cmn2nGDRlLduORNI7rAwv3FadvDmSNZu5MWlOUn+Dd7m37O7NmPRlyw8w519OlX/TEOdYfiJV/nsLtjFmyQ6K58/JV/3CuKnKDUEM2JjASeoVuF4NdCDGBMRlVX4tuG+KM9eOH3/vPcnAKWvZcfQcvcNKM7RzdfLn9H3C15j0Jqkjd6vgTNJWjstH7t4cmLCMSQFbfnB67Jw/BjcNhlYDEz+W//NWPvttJzfmz8nX/cJobVW+yYCSeqhnCvApMBaID1w4xqSACyfhx6GwdgIUqwn3TU60yv9r70kGXaryyzCsczXyWZVvMqikJv44Vf0koJEYkxK2/Ogcyz93FFo/D60HJVrlv/PzVsb+tpMSBXLxzcNhtKpsVb7J2JI8V4+I9AdmcHmvnhMBicqY5Lqyyr93EoTW87vK6j0nGTR1LTuPnuPeJmUY2smqfJM5JDXx93F/DvJYpkCFxFZ0Z/ZcBRxQ1dtFpDAwCed8wW6gl6qeTGrAxlxl609OlR95xKnwWz+faJX/n/lbGLt0F6EFcvHtw01oWbloEAM2JnUltVdP+et4jX8Bm4CL0zsMAX5R1ZEiMsR9PPg6tm8yqwun3Cr/OyhWA3pPgND6fldZvecEg6aEs/PYOe5rUoahna1fvsl8knIFrtxAZVVd67GsDBCvqgcSWbcUzjz+rwPPuou7Am3c+18Bi7HEb5Jr63yY87RT5bcaCDc9D1lz+Gx+Icap8j9f5lT54x9pQotKVuWbzCkppU4sMF1E6qjqOXfZWGAY4DfxA+8BzwP5PJYVV9UIAFWNEJFi3lYUkUeBRwHKlCmThDBNpnDhFPw0DNaMd6r8e76Dkg38rrJq9wkGTQ1n17Fz3N+0DEM6WZVvMrekXHoxVkRmAHcD49xq/wZVXeVvPRG5HTiiqqtFpE1yA1PVMcAYgEaNGtkcQQa2/Qyzn4bIw0mu8kfP38K4ZbsoWTAX3z3ShOZW5RuT5JO7Y4HPgHHAg8AXSVinBXCHiHQGcgL5ReRb4LCIlHCr/RLAkWuI22QmF07BTy/Amm/hhupwz/hEq/w/d5/gebfKf6BpWYZ0qkYeq/KNAZJ+cneziFwcwdsbaJmEdYYCQwHcin+gqt4vIqNwegmNdH/OurbQTaZwWZX/nDMCN5Eqf9RPW/hiuVvl/6MJzStalW+Mp+SUQJ/jVP7h19n9ciQwWUQeBvYCPa9jWyajijrtHMv/+1u4oVqSqvyVu04waOpa9hw/b1W+MX4k569iMvA+8FpyX0RVF+P03kFVjwPtkrsNk4lsW+D02Dkb4cyV32aI3yr/fEwcb/+4ha9+302pQs58+c0q2nz5xviS5MSvqueBAgGMxWR2V1b5d38DJRv6XWXFzuM8PzWcvSfO07d5OZ7vWJXc2a3KN8Yf+wsxacP2Bc6x/LMR0PIZZ878bDl9Nj8XHcfbP27mq9/3ULZIbiY+alfFMiapLPGb1BV1Gua/CH99DUWrwsMLoJT/Kn/5jmMMnhbO/pMXeKhFOQbdalW+Mclhfy0m9Wz/BWb/06nyWwxwrn+bSJU/8ofNfLNiD+WK5GbSo80IK184ePEak0FY4jfBF3UG5r/gVvlV4OGfoVQjv6ss336M56eFc+DUBR5uWZ6BHaqSK3uWIAVsTMZiid8E1/Zf3GP5B6HFv6DNML9VfmR0HG/O28T4P/ZSvmgepjzWjEblrMo35npY4jfBEXXGPZb/lVPl95sPpRv7XWXpNudY/sHTF/hHq/I816EqObNZlW/M9bLEbwJvx0Knyj9zAJo/DW1f8Fvln42K5Y15m5mwci8VbsjD1Meb0bCsVfnGpBRL/CZwos7Azy/B6i+hSOUkVfnLdxxj0JRwIk5f4LHWFXjmlipW5RuTwizxm8DYscjpsXOpyh8G2XL5bB4VG8/on5yrYpUvmoepTzSnQZlCQQzYmMzDEr9JWdFnYf5LsPoLKFIJ+v0EpcP8rrLh4GmembSGrYcjeaBpWYZ2rmb98o0JIPvrMinnYpV/ej80/6d7LN93lR+foPx3yQ7e/XkrhXJn58uHGtOmqtfr8hhjUpAlfnP9os/Czy/DqnFOlf/w/ESr/L3Hz/Ps5DWs2nOS22qXYES3WhTK4/sC6caYlGOJ31yfnb/CrKfg9D5o9hTc/KLfKl9VmbxqH6/N2UhIiPDe3fXoWi8UEQli0MZkbpb4zbW5ssrv9xOUaeJ3laNnoxk6PZwFm47QvGIRRvesS2hB3/8kjDGBYYnfJN/OX2H2U3AqaVU+wI/rI3hhxnrORsfx0u01eKh5OUJCrMo3JjVY4jdJFx3pVvmfQ+GK0O9HKNPU7yonzsUwfPYG5qw9SK2S+ZnQqx5ViucLUsDGGG8s8Zuk2bUEZj3pVPlNn3Sq/Oy5/a7y4/pDvDhzHacvxPLcLVV4vE1FsmUJCVLAxhhfLPEb/6IjYcFw+HMsFK4AD/0AZZv5XeWkW+XPXnuQmqH5+ebhJlQvkT9IARtjEmOJ3/i26ze3yt8LTfvDzS8lWuX/tOEQL8xYz6nzMTzTvgr921qVb0xaY4nfXC06Eha8An9+luQq/9T5GF6ZvYGZaw5So0R+vu4XRo1Qq/KNSYss8ZvLeVb5TZ6Adi8nWuX/vPEww2as4+S5GAa0r8yTbStZlW9MGmaJ3zhizjlV/soxUKg8PDQPyjb3u8qJc06VP3vtQaqXyM+XDzWmZmiB4MRrjLlmlvgN7F7qVPkn9yS5yp+3LoKXZ63n9IVYnmlfhSfaVCR7VqvyjUkPLPFnZjHnYMGrsPK/UKgc9P0eyrXwu8rRs9EMn72eeesOUbtkAb59pAnVbrRj+cakJ5b4M6vdy2BWfzi5G5o87lb5eXw2V1VmrTnIq3M2cC46nuc7VuXRVhXIasfyjUl3LPFnNjHn4JfX4I9P3Sp/XqJV/sFTF3hhxjoWbTlK/TIFGdWjDpWK2ehbY9IrS/yZye5l7rH8XRD2GLQf7rfKT0hQxq/cy1s/bCY+QXn59hr0aV6OLDbHjjHpmiX+zOBSlf9fKFTWPZbf0u8qO45GMnTaOlbuPkGrykV5487alC7s/4SvMSZ9sMSf0e1ZDjP7u1X+o9D+Fb9Vfmx8Ap/9tpP3FmwjZ9YQRvWoQ4+GpWy+fGMyEEv8GVXM+f8dyy9YBvrMhfKt/K6y/sBpnp8azsaIM3SqdSOvdq1JsXw5gxSwMSZYLPFnRHt+d3rsnNgJjf/hVPk58vpsHhufwP8t3M5Hi7ZTOE92Pr2/AR1rlQhevMaYoLLEn5HEnIeF/4YVn0DB0tBnDpRv7XeVrYfP8uzkNaw/cIbu9UsyvEtNCuTOFqSAjTGpwRJ/RrF3hXMs/8QOaPwItH/Vb5Ufn6B8vnQno+dvJW+OrFblG5OJWOJP72LOw8IRsOLjJFf5e4+fZ+CUtazcfYJbahTnjTtrc0O+HEEK2BiT2izxp2eeVX6jh+GW1/xW+arKhJX7GPH9RrKI8J+edeneoKT12DEmkwlY4heR0sDXwI1AAjBGVd8XkcLAJKAcsBvopaonAxVHhhR7wanyf/8ICpSGB2dDhZv8rnL4TBSDp4WzeMtRWlQqwts96lKyoP8LpBtjMqZAVvxxwHOq+peI5ANWi8jPQF/gF1UdKSJDgCHA4ADGkbHs/cPpsXN8u1vlvwo5/E+fMGftQV6cuZ7ouHhe6VKDB5uVI8RG3xqTaQUs8atqBBDh3j8rIpuAkkBXoI3b7CtgMZb4E3dVlT8LKrTxu8qp8zG8NGsDc9YepG7pgrzTqy4Vb/B9KMgYkzkE5Ri/iJQD6gN/AMXdfwqoaoSIFPOxzqPAowBlypQJRphp176VMPMJt8rv5x7L91/l/7j+EC/PWs+JczE8d4szX77NpGmMgSAkfhHJC0wDBqjqmaSeSFTVMcAYgEaNGmngIkzDYi/AotedKj9/ySRV+YfPRDF81gZ+3HCI6iXyM65vY2qVtKtiGWP+J6CJX0Sy4ST98ao63V18WERKuNV+CeBIIGNIt/atdHrsHN8GDR+CDv/2W+WrKhP/3Mcb8zYRHZfA8x2r8o9WFezat8aYqwSyV48AnwObVPUdj6dmA32Ake7PWYGKIV26ssp/YCZUbOt3lZ1HIxk6fR1/7DpB0wqFebN7HcoX9T0RmzEmcwtkxd8CeABYJyJr3GXDcBL+ZBF5GNgL9AxgDOnLvj+dHjvHtkLDvnDLvyGn78saxsYnMGbJTt7/ZRs5soYwsntt7m5c2vrlG2P8CmSvnqWArwzULlCvmy7FRrlV/oeQLxQemAEVb/a7Svj+Uwyeto5NF2fSvKMmxfLbTJrGmMTZyN3Utn+V02Pn2FZo0Ac6jPBb5Z+PiePdn7fy+dJdFM2bg0/vb0jHWjcGMWBjTHpniT+1xEbB4jdg+f85Vf7906GS/y9Cv207yrAZ69h34gL3NinD4I7VKJDLZtI0xiSPJf7UsH+1W+VvgQYPulW+7y6XJ8/FMOL7TUz7az8ViuZh0qNNaVKhSBADNsZkJJb4gyk2Cha/Ccs/gHwl4P5pUKm9z+aqypzwCF6dvYHTF2J5qm0lnrq5EjmzZQli0MaYjMYSf7B4Vvn1H4BbX/db5R84dYGXZq5n4eYj1C1VgG8faUL1Er6P/RtjTFJZ4g+0uGinyl/2vlPl3zcNKvuu8hMSlG9W7OHtHzeToPDS7TXo27wcWWxSNWNMCrHEH0gHVjujb49uTlKVv/XwWYZMC+evvadoXeUGXu9Wi9KFcwcxYGNMZmCJPxA8q/y8N8J9U6HyLT6bR8fF8/GiHXy8eDt5c2Tl3bvr0q2eXSDFGBMYlvhT2mVV/v3Q4XXIVdBn89V7TjJkWjjbjkTSrV4oL91egyJ57TKIxpjAscSfUuKi4de3YOl7kLd4olV+ZHQco37czNcr9hBaIBdfPNSYtlW9zlBtjDEpyhJ/Sjjwl1vlb4J69zvH8v1U+Qs3H+aFGes5dCaKvs3LMbBDVfLksI/CGBMclm2ux2VVfjG4dwpU6eCz+bHIaF6ds5E5aw9StXg+Pr6vAfXLFApevMYYgyX+a3fwb6fKP7IR6t0Ht77hs8pXVab/dYB/f7+R89HxPHdLFR67qSLZs9pc+caY4LPEn1xx0fDr27D0XbfKnwxVbvXZPOL0BYZNX8eiLUdpVLYQI++qQ6Vidt1bY0zqscSfHJ5Vft17oeMbkMv7oRpVZfKqfYyYu4m4BGV4lxr0aVaOEBuIZYxJZZb4kyIuBpa8Db+9A3luSLTKP3DqAkOmhfPbtmM0rVCYt+6qQ9kidkUsY0zaYIk/MQfXuFX+BqjbGzq+6bfK/27lXt6ctxlV5d/danFfWBmr8o0xaYolfl/iYmDJKPjtP06V33sSVO3os/m+E+cZPC2c5TuO07JSUd7sXtumWzDGpEmW+L3xrPLr3AOdRvqs8hMSlG//2MPIHzYTIsKb3Wtzj1331hiThlni9xQXA7+Ndqr83EWg90So2sln893HzjF4Wjh/7DrBTVVu4M3utQktmCuIARtjTPJZ4r8oItyp8g+vc6r8jm9C7sJem8YnKF8u382onzaTLUsIo3rUoUfDUlblG2PSBUv8V1b590yAap19Nt95NJLnp4azas9J2lUrxut31ubGAjmDGLAxxlyfzJ34L6vy74aOI/1W+eOW7mL0/C3kzJaFd3rV5c76NnWyMSb9yZyJPy4Glr7j9NpJQpW//chZBk0N5++9p+hQozgjutWiWH6r8o0x6VPmS/yH1sGMJ5wqv3Yv6PSWzyo/Lj6Bz37bxbsLtpInexY+6F2fLnVKWJVvjEnXMk/ij491juMvGQW5CsPd46H67T6bbzh4mqHT1xG+/zSdat3Ia11rcUM+u0CKMSb9yxyJ/9A6mPmE87N2T+j0ts8q/3xMHO8t2MbnS3dRKHd2Prq3AbfVKRHkgI0xJnAyduKPj3Xm11nydpKq/EVbjvDijPUcOHWB3mGlGdKxOgVyZwtiwMYYE3gZN/EfWu9W+eFQqwd0HuWzyj9yNorX5mxkbngElYrlZfJjzQgr772tMcakdxkv8cfHOnPl//q2c2GUu7+F6l28No2LT2D8H3v5z/wtRMUm8OwtVXjspgrkyJoluDEbY0wQZazEf3iDU+VHrIVad0GnUZCniNemK3Ye55XZG9h86CwtKhXhta61qHiDXSDFGJPxZYzEHx/rXPf217ecKr/XN1DjDq9NI05f4I15m5mz9iAlC+bik/sa0LHWjdZF0xiTaaT/xO9Z5dfsDp1He63yz0TFMubXnYxdupMEhafbVeaJmyqSK7sd1jHGZC7pN/HHx8Gyd2HxW5CzAPT6Gmp0vapZdFw8367Yy4cLt3HyfCxd6oby/K1Vba58Y0ymlT4T/+GNbpW/xmeVHxkdx3d/7GHsb7s4cjaalpWKMrhjNWqXKpA6MRtjTBqRKolfRDoC7wNZgLGqOjJJK15Z5ff8Cmp2u6zJ3uPnmfjnXsb/sZfTF2JpXrEI795djxaViqb02zDGmHQp6IlfRLIAHwG3APuBP0Vktqpu9LviZVX+nW6V7yTzY5HRLNx8hFlrDrBs+3FCBG6pUZwn2lSiXumCgX1DxhiTzqRGxR8GbFfVnQAiMhHoCvhO/JGH0TE3QY58nOnyObuKtWf31nOE79/I6r0nCd9/ClUoWTAXz91ShZ6NStsc+cYY40NqJP6SwD6Px/uBJn7XOHOQ72Pa8fK5hzgxJRewDICc2UKoGVqAAe2q0K56MWqG5rdumcYYk4jUSPzeMrNe1UjkUeBR92F0lxE/r4efr1pxCzA9RcO7JkWBY6kdhB8W3/VJy/Gl5djA4rte1xtfWW8LUyPx7wdKezwuBRy8spGqjgHGAIjIKlVtFJzwks/iuz4W37VLy7GBxXe9AhVfSEpvMAn+BCqLSHkRyQ7cA8xOhTiMMSZTCnrFr6pxIvIU8BNOd85xqroh2HEYY0xmlSr9+FV1HjAvGauMCVQsKcTiuz4W37VLy7GBxXe9AhKfqF51XtUYY0wGlhrH+I0xxqSiNJ34RaSjiGwRke0iMiSVYigtIotEZJOIbBCRf7nLXxGRAyKyxr119lhnqBvzFhG5NQgx7haRdW4cq9xlhUXkZxHZ5v4slBrxiUhVj320RkTOiMiA1Nx/IjJORI6IyHqPZcneXyLS0N3v20XkA0mhQSQ+4hslIptFJFxEZohIQXd5ORG54LEfP02l+JL9eQY5vkkese0WkTXu8qDuPz/5JLi/f6qaJm84J353ABWA7MBaoEYqxFECaODezwdsBWoArwADvbSv4caaAyjvvocsAY5xN1D0imVvA0Pc+0OAt1Irvis+00M4fYtTbf8BrYEGwPrr2V/ASqAZztiUH4BOAYyvA5DVvf+WR3zlPNtdsZ1gxpfszzOY8V3x/H+Al1Nj/+E7nwT19y8tV/yXpnZQ1Rjg4tQOQaWqEar6l3v/LLAJZ/SxL12Biaoaraq7gO047yXYugJfufe/Arp5LE+t+NoBO1R1j582AY9PVZcAJ7y8bpL3l4iUAPKr6u/q/BV+7bFOisenqvNVNc59uAJn/ItPwY7PjzSx/y5yq+JewAR/2whUfH7ySVB//9Jy4vc2tYO/hBtwIlIOqA/84S56yv3qPc7jq1lqxK3AfBFZLc6IZ4DiqhoBzi8bUCwV47voHi7/g0sr+w+Sv79KuveDHSdAP5wK76LyIvK3iPwqIq3cZakRX3I+z9Taf62Aw6q6zWNZquy/K/JJUH//0nLiT9LUDsEiInmBacAAVT0DfAJUBOoBEThfHyF14m6hqg2ATsCTItLaT9tU2a/iDNa7A5jiLkpL+88fX/Gk1n58AYgDxruLIoAyqlofeBb4TkTyp0J8yf08U+tz7s3lxUeq7D8v+cRnUx9xXFd8aTnxJ2lqh2AQkWw4H9J4VZ0OoKqHVTVeVROAz/jf4Yigx62qB92fR4AZbiyH3a+DF7+2Hkmt+FydgL9U9bAba5rZf67k7q/9XH64JeBxikgf4HbgPvfrPe4hgOPu/dU4x4CrBDu+a/g8U2P/ZQW6A5M84g76/vOWTwjy719aTvxpYmoH95jg58AmVX3HY3kJj2Z3Ahd7EMwG7hGRHCJSHqiMcxImUPHlEZF8F+/jnARc78bRx23WB5iVGvF5uKzSSiv7z0Oy9pf7dfysiDR1f0ce9FgnxYlz8aLBwB2qet5j+Q3iXOMCEangxrczFeJL1ucZ7Phc7YHNqnrpEEmw95+vfEKwf/+u9yx1IG9AZ5yz3juAF1IphpY4X6HCgTXurTPwDbDOXT4bKOGxzgtuzFtIoZ4KfuKrgHPWfy2w4eJ+AooAvwDb3J+FUyM+9/VyA8eBAh7LUm3/4fwDigBicSqnh69lfwGNcBLcDuBD3AGRAYpvO86x3ou/g5+6be9yP/e1wF9Al1SKL9mfZzDjc5d/CTx+Rdug7j9855Og/v7ZyF1jjMlk0vKhHmOMMQFgid8YYzIZS/zGGJPJWOI3xphMxhK/McZkMpb4jTEmk7HEbzIFESkoIv09HoeKyNQAvVY3EXn5Gtdd4DklrzGBYP34TabgTog1V1VrBeG1luOMsD12Dev2AUqp6uspH5kxDqv4TWYxEqgozsU2RolzAY71ACLSV0RmisgcEdklIk+JyLPujI0rRKSw266iiPzozoL6m4hUu/JFRKQKEH0x6YvIl+5FMpaLyE4R6eEuLyEiS9x41nvMCjkbZ3oLYwLGEr/JLIbgXAugnqoO8vJ8LeBenMnFXgfOqzNj4+8486CAc+Hrf6pqQ2Ag8LGX7bTAGfrvqQTOUP3bcf4B4b7WT6paD6iLM3QfVT0J5BCRItfwHo1JkqypHYAxacQidS6McVZETgNz3OXrgDruNLrNgSnyvyvc5fCynRLA0SuWzVRn1sqNIlLcXfYnMM6dqXGmqq7xaH8ECMWZ38iYFGcVvzGOaI/7CR6PE3AKpBDglPuN4eKtupftXABy+tm2wKWrRLUGDgDfiMiDHm1yutsxJiAs8ZvM4izONU6viToXy9glIj3BmV5XROp6aboJqJTY9kSkLHBEVT/Dmaa3wcXtAjfiXEfZmICwxG8yBXUutrHMPZE66ho3cx/wsIhcnALb2zWglwD1xeN4kA9tgDUi8jfO1MDvu8sbAiv0f9fXNSbFWXdOY1KYiLwPzFHVBde47mxV/SXlIzPGYRW/MSnvDZyLz1yL9Zb0TaBZxW+MMZmMVfzGGJPJWOI3xphMxhK/McZkMpb4jTEmk7HEb4wxmcz/A7RX9uuB0VdTAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "selection = df[['time', 'displacement_sq']].groupby(\"time\").mean()\n",
    "\n",
    "lm_time = np.arange(0, 2000)\n",
    "\n",
    "# line_y = np.linspace(0, np.max(selection[\"displacement_sq\"]), 2000)\n",
    "\n",
    "lm = np.poly1d(np.polyfit(df[\"time\"], df[\"displacement_sq\"], 1))\n",
    "   \n",
    "def predict_line(m, b, n):\n",
    "    return [m*i + b for i in range(n)]\n",
    "\n",
    "line_y = predict_line(m=lm.coef[0], b=0, n=len(lm_time))   \n",
    "\n",
    "d = calc_d_from_step(step=lm.coef[0])\n",
    "d_text = f\"meas D = {d:.2e}\"\n",
    "exp_text = f\"exp D = {convert_d_from_cm2_s_to_nm2_ns(d=1.88e-9):.2e}\" # 1.88e-9 cm2/s\n",
    "    \n",
    "step_text = f\"step_nm = {round(lm.coef[0], 5)} nm\"\n",
    "\n",
    "print(d_text)\n",
    "print(exp_text)\n",
    "print(step_text)\n",
    "\n",
    "# plot the displacement_sq against the fit\n",
    "\n",
    "plt.suptitle(\"displacement of LHCII from origin\")\n",
    "plt.title(\"D in nm^2/ns units\")\n",
    "plt.xlabel(\"time (ns)\")\n",
    "plt.ylabel(\"<nm>^2\")\n",
    "plt.xlim(0, np.max(df[\"time\"]) * 1.05)\n",
    "plt.ylim(0, np.max(selection[\"displacement_sq\"]) * 1.05)\n",
    "plt.annotate(text=d_text, xy=(20, np.max(selection[\"displacement_sq\"]) * 0.95))\n",
    "plt.annotate(text=step_text, xy=(20, np.max(selection[\"displacement_sq\"]) * 0.85))\n",
    "plt.annotate(text=exp_text, xy=(20, np.max(selection[\"displacement_sq\"]) * 0.90))\n",
    "plt.plot(selection[\"displacement_sq\"])\n",
    "plt.plot(lm_time, line_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "meas D = 2.85e-04\n",
      "exp D = 1.88e-04\n",
      "step_nm = 0.04771 nm\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x17f3eb8b7c0>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEjCAYAAAA1ymrVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA7BElEQVR4nO3dd3hUZfbA8e9JgNCbgCa0gJHepISiKCyKgCiIgKKuYlnbuv4sIMXuYnctrK6uYpelg4KLWGEVBRGQFHpvCb0G0nN+f9wbHEISJmUyk8z5PM88mbn1zL2TM+/c+95zRVUxxhgTPEL8HYAxxpiSZYnfGGOCjCV+Y4wJMpb4jTEmyFjiN8aYIGOJ3xhjgowl/lJERD4SkQnu854isr4YlrlNRC4renSlnzg+FJHDIrLM3/HkpzTFCiAiN4rIN8U9rSkcS/yllKr+pKrN/R1HIBGRkSKyuAiLuBi4HGigqtEFWb6ILBKRO3IM6yUiu3IMu0JEfhSR4yKyX0T+JyJX57b8s3wp5xtroFHVyarat7inNYVjid+YPzQGtqnqCV8sXESGAjOAT4AGwLnAE8BVhVhcvrGKSLnCxlncAikW47DEH8BE5EIRWem2DqcBFT3GndaaFJExIrLbnXa9iPRxhz8lIjNFZJo7bqWItM9jfdEiskREjohIooi8KSIVPMa3FpFvReSQiOwVkfHu8BARGSsim0XkoIhMF5Ha7rhIEVERuVVEdrqHJu4WkS4iEuuu680ccdwmImvdab8WkcYe49Sdf6M7/i33sEdL4B2gu4gkiciRPN5jhIjMdd/DJhH5izv8dmCSx/xPF2xv5U9EBHgV+LuqTlLVo6qapar/U9W/FHBZZ8Sa/XlwPwd7gA9FJExEXheRBPfxuoiEucvInv4REdnn7u/BIjJARDa422d8PjHUEJFP3F8t20XkMREJcceNFJGfReQ1ETkEPJXLr5m+7uf0qIj8y/3lc4fH/J7T5rrPC7LNzOks8QcoN+F+DnwK1MZpKV6bx7TNgfuALqpaDbgC2OYxySB3/trAf4DPRaR8LovKBB4E6gDdgT7Ave46qgHfAQuACCAK+N6d735gMHCpO+4w8FaOZXcFLgCuA14HHgUuA1oDw0XkUnc9g4HxwBCgLvATMCXHsgYCXYD2wHDgClVdC9wNLFHVqqpaM7dt5S5rlxvnUOA5Eemjqu/nmP/JPOYvrOZAQ2BmUReUT6zn4ezjxsCdONu4G9ABZ1tFA495LOo8nMZEfZxfHu8BNwGdgJ7AEyLSNI8w/gnUAJri7PebgVs9xncFtgD1gGc9ZxSROjjbYRxwDrAe6HGWt33GPj/L9CY/qmqPAHwAlwAJgHgM+wWY4D7vBexyn0cB+3ASafkcy3kKWOrxOgRIBHq6r7cBl+URwwPAHPf5COD3PKZbC/TxeB0OpAPlgEhAgfoe4w8C13m8ngU84D7/Crg9R7wngcbuawUu9hg/HRjrPh8JLM5nmzbE+XKr5jHseeAjL+fPczywyI3ziMcjyWMfXeTGXtHb5Z9l3+SctheQ5rl8YDMwwOP1FTiHh7KnTwZC3dfV3Pi6eky/Ahicy7pDgVSglcewu4BFHrHtyCtenC+JJR7jBNgJ3JHHe8tzn9ujcA9r8QeuCGC3up901/bcJlTVTThJ+ilgn4hMFZEIj0l2ekybxR8t3tOISDMR+VJE9ojIMeA5nNY/OElzcx6xNgbmuIdtjuB8EWTiHMPOttfjeXIur6t6LOsNj2UdwkkM9T2m3+Px/KTHvGcTARxS1eMew7bnWHZR3K+qNbMfOK3UbAfdv+HFtK7c7FfVFI/XEZz+mdnO6fv9oKpmus+T3b957RdPdYAKuSzbczvuJG8RnP6ZVJzPZH4Ku89NLizxB65EoH6OY5mN8ppYVf+jqhfjJE4FXvQY3TD7iXsctgHOr4mc3gbWAReoanWcQy7Z698JnJ/H6ncC/T2TnqpWVNXd+b7DvJd1V45lVVLVX7yY92ylZhOA2u5hq2yNgMLEWVDrcd5brofriknO95+A83nI1ojc93tBHcD5RZdz2Z7bMb99kYjzGQROnf9okPfkprhZ4g9cS4AM4H4RKSciQ3CO0Z5BRJqLyJ/cE3cpOC21TI9JOonIEHF6VzyA8zN9aS6LqgYcA5JEpAVwj8e4L4HzROQB96RhNRHp6o57B3g2+ySsiNQVkUGFfN/vAONEpLW7rBoiMszLefcCDcTjhLQnVd2Jc7jseRGpKCLtgNuByQWIT9x5Tz28mclt1T4EPC7Oie7q4pwUv1hE3i3A+gtiCvCYuz/q4BzH/6yoC3V/JUzH2efV3P3+UAGW/V+grXsyuRzwV5zzDaaEWOIPUKqahnOCcyTOydLrgNl5TB4GvIDTEtuDc0LNs0fGF+78h4E/A0NUNT2X5YwCbgCO45zom+YRz3GcfuNXuevYCPR2R78BzAW+EZHjOF8qXSkEVZ2D82tlqnu4KR7o7+XsPwCrgT0iciCPaUbgnHdIAOYAT6rqtwUIsQfOF+uph3jZXVFVZ+Lsh9vc9e8FJuDsH1+YACwHYoE4YKU7rDj8DTiBcwJ3MU6ngQ+8mVFVDwDDgJdwDoG1cuNMLabYzFnI6YeQTVkjIk8BUap6k79jMSY37uHHXcCNqrrQ3/EEA2vxG2NKnDhXMNd0D09mn0vK7fCj8QFL/MYYf+iO00vsAM7hw8Gqmpz/LKa42KEeY4wJMtbiNwFFRL4SkVv8HYcnEekkIk/muDai1ArEbWxKlrX4jc+IyDaci7gycLqXrsEpUPaueyFZwHO7tf6Ec31DDeASVT3iMf4WnJIVF+B0hf0PMF5VMzymiQCWqWrA9VUXkZE4V8xe7O9YTMmxFr/xtavUqR/UGKfL6Rjgff+G5B0RaYBTQmIcTgmNhcDcHH33K+NcG1EHpwtrH5xusZ4G4NQ4MiYgWOI3JUKdapRzcfqx3yIibXKbTjzq2mdXaRSRV9yqjFtFJM8+/eLUrx8lTtXPo+JUJK3ojitQNUpxqov+F3hcnWqaqqr/h3MB2FQRCXXf19vq3Bshzb1SeTJOXR5PA4D5XsRYxy2ZccSN5ye3q2PO95ld8bScxzCvt1v2tJJHRVN3m6wRp5rrbhHJ+UVmSjlL/KZEqeoynD7bPb2cpStOuYM6OBf8vO9e4p+X4UA/oAnQDucCuGxeV6NU1UOq2l5VT7saVVXHqupgjxo3OV2CcxEZAOJUQb0E8LxILK8YH8bZNnVxDpGN5+xlKPJy1u2meVc0fR+nbEY1oA3OhXGmDLHEb/whAad8sDe2q+p7bqL9GKfI2bn5TD9RVRNU9RAwD6ckcbZ04Fn3quWpOEnxDVU9rqqrcRJ2u4K9lT+IyK1AZ+AVj8GXADE5CsPlFWO6+/4aq2q6+0uisIm/oNvNUzrQSkSqq+phVV1ZyBhMgLLEb/yhPk7VTW+cqsqoqifdp/lVZsyvimNhq1GelTj3EXgBp1idZ7mIU4d5vIjxZWATTumLLSIytjCx5FyHl9vN07U4cW8X5wYp3YsQhwlAlvhNiRKRLjiJvyj3xg0oItIP57DRVaoal2P0AJxzBWfl/vJ4WFWb4lzU9JC4d1LLIft2i5U9hhW2yNkZvyhU9TdVHYRT8+lznIJspgyxxG9KhFuNciDOIZbPckmQpZKI/AnnhO617vkLz3FNgDBVXeflsgaKSJR7LP4YThfYM84lqOp+nBLIN4lIqIjcRt4ls8/mtIqmIlJBRG4UkRruIbHsOEwZYonf+No8t2LnTpxbAb7K6bfoK+0ex+nfP9/tGZMkIl+5467kzMM8+bkA5/aWSThluf+lqovymPYvwGic6patcXobFUZuFU3/DGxzq6PejXMC3JQhdgGXMT4iIvOBN1W1IMnfGJ+zFr8xvrMI56IvYwKKtfiNMSbIWIvfGGOCjFe3jPO3OnXqaGRkpL/DMMaYUmXFihUHVLVuzuGlIvFHRkayfPlyf4dhjDGliohsz224HeoxxpggY4nfGGOCjCV+Y4wJMpb4jTEmyFjiN8aYIGOJ3xhjgowlfmOMCTKW+I0xJgDtO5ZCaoZvKmJb4jfGmACzYvshrvznYp6f79WtHArMEr8xxgQIVeXTJdu4/t2lVK4QyojoRj5ZT6ko2WCMMWXdybQMxs+O4/NVCfypRT1eG96BGpXL+2RdlviNMcbPNu9P4p7PVrBxXxKj+jbj3l5RhISIz9Znid8YY/zoq7hERs+MpUK5ED65LZqeF5xRTLPYWeI3xhg/SMvI4oWv1vHBz1tp37Amb9/YkYialUpk3Zb4jTGmhO08dJL7/rOSmF1HGdkjknEDWhBWLrTE1m+J3xhjStBXcYk8MisWgHdu6ki/NuElHoMlfmOMKQEp6Zk8N38tnyzZTvsGNXjzho40rF3ZL7FY4jfGGB/bsj+Jv035ndUJx7j94iaM6deCCuX8dxmVJX5jjPERVWXWyt088UU8FcqF8N7Nnbm81bn+DssSvzHG+MLxlHQe/zyez1cl0LVJbd64/kLOq1HR32EBlviNMabYxew8wv1Tf2fnoZM8dHkz/to7ilAfXpBVUJb4jTGmmGRlKZMWb+GlBeupVy2MaXd1p0tkbX+HdQYr0ubq1asXzZs3p127drRo0YL77ruPI0eOFHm569ato3v37oSFhfHKK6/kOd33339Px44d6dChAxdffDGbNm0C4OjRo1x11VW0b9+e1q1b8+GHHxY4hueff56oqCiaN2/O119/fcb4q6++mjZt2hR4ucaYP2w/eIIbJi3lufnr6NOyHvP/r2dAJn2wxH+ayZMnExsbS2xsLGFhYQwaNKjIy6xduzYTJ05k1KhR+U53zz33MHnyZFatWsUNN9zAhAkTAHjrrbdo1aoVMTExLFq0iIcffpi0tDSv179mzRqmTp3K6tWrWbBgAffeey+ZmX/U+J49ezZVq1Yt3JszxpCVpXyweCv9Xv+J1buP8cKQtrxzUydqVq7g79DyVCoT/2effUZ0dDQdOnTgrrvuIjMzk99++4127dqRkpLCiRMnaN26NfHx8SxatIhLLrmEa665hlatWnH33XeTlZWV7/IrVKjASy+9xI4dO4iJiSlSrPXq1aNLly6UL59/lT0R4dixY4DTyo+IiDg1/Pjx46gqSUlJ1K5dm3LlyuW5HXL64osvuP766wkLC6NJkyZERUWxbNkyAJKSknj11Vd57LHHivQejQlWm/cnMfzfS3jmyzV0a1qbbx66hOujGyESOMfzc1PqjvGvXbuWadOm8fPPP1O+fHnuvfdeJk+ezM0338zVV1/NY489RnJyMjfddBNt2rRh0aJFLFu2jDVr1tC4cWP69evH7NmzGTp0aL7rCQ0NpX379qxbt4727dufNu66665j/fr1Z8zz0EMPcfPNNxfqfU2aNIkBAwZQqVIlqlevztKlSwG47777uPrqq4mIiOD48eNMmzaNkJCQfLeDp927d9OtW7dTrxs0aMDu3bsBePzxx3n44YepXNk/F5EYU1plZimTftrCq99uIKxcCP8Y1p4hHesHfMLPVuoS//fff8+KFSvo0qULAMnJydSrVw+AJ554gi5dulCxYkUmTpx4ap7o6GiaNm0KwIgRI1i8ePFZEz84fXBzM23atKK+jTO89tprzJ8/n65du/Lyyy/z0EMPMWnSJL7++ms6dOjADz/8wObNm7n88svp2bNnvtvhbO9BRFi1ahWbNm3itddeY9u2bcX+fowpqzbsPc7ombHE7DxC31bnMmFwG+pVD4xumt4qdYlfVbnlllt4/vnnzxh36NAhkpKSSE9PJyUlhSpVqgCc8S3szbdyZmYmcXFxtGzZ8oxxxd3i379/PzExMXTt2vXU8vv16wfAhx9+yNixYxERoqKiaNKkCevWrctzO8yZM4enn34acH5FNGjQgJ07d54av2vXLiIiIliyZAkrVqwgMjKSjIwM9u3bR69evVi0aFGB4zcmGKRnZvHv/21m4vebqBIWysQRF3JVu/BS08o/jaoG/KNTp06abfXq1RoVFaV79+5VVdWDBw/qtm3bVFX1qquu0smTJ+uECRP0r3/9q6qqLly4UCtWrKhbtmzRzMxM7du3r86cOVNzuvTSS/W3335TVdW0tDQdPXq0XnLJJWdMV1hPPvmkvvzyy7mOS09P13POOUfXr1+vqqqTJk3SIUOGqKrq3XffrU8++aSqqu7Zs0cjIiJ0//79+W4HT/Hx8dquXTtNSUnRLVu2aJMmTTQjI+O0abZu3aqtW7currdqTJmzevdRHfDGj9p4zJd67+QVuv94ir9D8gqwXHPJqaWuxd+qVSsmTJhA3759ycrKonz58rz11lv873//o1y5ctxwww1kZmbSo0cPfvjhB0JCQujevTtjx44lLi7u1Ine3Nx4442EhYWRmprKZZddxhdffFHkePfs2UPnzp05duwYISEhvP7666xZs4bq1aszYMAAJk2aREREBO+99x7XXnstISEh1KpViw8++ABwjsOPHDmStm3boqq8+OKL1KlThzp16uS6HRo3bnza+lu3bs3w4cNp1aoV5cqV46233iI0tOTKvxpTmqVlZPHmwk38a+Emalau4LdqmsVNNI/j2IGkc+fOunz58kLNu2jRIl555RW+/PLLYo7KGFOWxe46wugZsazfe5xrLqzPEwNbUatK4HbRzI2IrFDVzjmHl7oWvzHG+FJKeiZvfL+Rd3/cQp2qFXj/ls70aen/wmrFqcwn/l69etGrVy9/h2GMKQVWbD/MIzNj2Lz/BNd1bsj4K1tSo1L+1+CURqXyAq7cvP7665w8edLfYRSbFStW0LZtW6Kiorj//vvz7Fpa0HIMDz74IB06dKBDhw40a9aMmjVrArBw4cJTwzt06EDFihX5/PPPAXjzzTeJiopCRDhw4ECxv1dj/C05LZMJX65h6Du/kJKexSe3RfPi0HZlMukDvu3VAzwIrAbigSlARaA28C2w0f1b62zL8ezVk5fGjRvr/v37i3D+O7B06dJFf/nlF83KytJ+/frp/Pnzz5hm9erVp/XYadq06Wk9dmbNmqUjRozIs8fOxIkT9dZbbz1j+MGDB7VWrVp64sQJVVVduXKlbt26tcxtY2NUVZduPqCXvvSDNh7zpT46J1aPJaf5O6RiQx69enzW4heR+sD9QGdVbQOEAtcDY4HvVfUC4Hv3tddOnDjBlVdeSfv27WnTpg3Tpk1j4sSJJCQk0Lt3b3r37g3AN998Q/fu3enYsSPDhg0jKSkJgMjISMaMGUN0dDTR0dGniqHlZuTIkdx///306NGDpk2bMnPmTMA5YXzppZcyfPhwmjVrxtixY5k8eTLR0dG0bduWzZs3F3RznSYxMZFjx47RvXt3RISbb775VOvbU1HLMUyZMoURI0acMXzmzJn079//1BW9F154IZGRkfnG/NFHHzFkyBD69evHBRdcwCOPPHJqXNWqVRkzZgydOnXisssuY9myZfTq1YumTZsyd+5cL7aIMcXvRGoGT34Rz3XvLiVL4T9/6cqEwW2pVrGMtvI9+PpQTzmgkoiUAyoDCcAg4GN3/MfA4IIscMGCBURERBATE0N8fDz9+vXj/vvvJyIigoULF7Jw4UIOHDjAhAkT+O6771i5ciWdO3fm1VdfPbWM6tWrs2zZMu677z4eeOCBfNeXmJjI4sWL+fLLLxk79o/vqJiYGN544w3i4uL49NNP2bBhA8uWLeOOO+7gn//85xnLyXkoJfvRo0ePM6bdvXs3DRo0OPXas8xCzukaNmyY63RnK8ewfft2tm7dyp/+9Kczxk2dOjXXL4SzWbVqFdOmTSMuLo5p06adunDsxIkT9OrVixUrVlCtWjUee+wxvv32W+bMmcMTTzxR4PUYU1S/bjlI/zd+4pOl2xnZI5IFD/Skx/l1/B1WifHZyV1V3S0irwA7gGTgG1X9RkTOVdVEd5pEETmzzgAgIncCdwI0atTo1PC2bdsyatQoxowZw8CBA+nZs+cZ8y5dupQ1a9Zw0UUXAZCWlkb37t1Pjc9OaiNGjODBBx/M930MHjyYkJAQWrVqxd69e08N79KlC+HhTn/e888/n759+56Kb+HChWcsp3fv3qxatSrfdWXTPMoseDudN+UYpk6dytChQ8/o05+YmEhcXBxXXHGFV7F66tOnDzVq1ACc6y22b99Ow4YNqVChwqkrkdu2bUtYWBjly5enbdu2Vi7ClKjktExe+nodH/68jcbnVGband2JbhKYpZN9yWeJX0Rq4bTumwBHgBkicpO386vqu8C74PTjzx7erFkzVqxYwfz58xk3bhx9+/Y9o9Woqlx++eVMmTIlr9hyfZ6bsLCw05ab2/CQkJBTr0NCQsjIyDhjOQsXLsz1S6Zy5cr88ssvpw1r0KABu3btOvU6u8xCTkUpxzB16lTeeuutM5Y5ffp0rrnmmrNWE82N5zYJDQ09tR3Kly9/ajt7s62M8YXl2w4xakYM2w6eZGSPSB7p15zKFcp8x8Zc+fJQz2XAVlXdr6rpwGygB7BXRMIB3L/7CrLQhIQEKleuzE033cSoUaNYuXIlANWqVeP48eMAdOvWjZ9//vnU8fuTJ0+yYcOGU8vILrI2bdq0034J+FJ2iz/nI2fSBwgPD6datWosXboUVeWTTz7J9d4AV199NVOnTiU1NZWtW7eyceNGoqOjueeee0hISGDbtm0sXryYZs2anZb0169fz+HDh3N973kd9zemtEpJd3rsDPv3EjKylCl/6cZTV7cO2qQPvu3HvwPoJiKVcQ719AGWAyeAW4AX3L8FqosQFxfH6NGjCQkJoXz58rz99tsA3HnnnfTv35/w8HAWLlzIRx99xIgRI0hNTQVgwoQJNGvWDIDU1FS6du1KVlZWnr8K/O3tt99m5MiRJCcn079/f/r37w/A3LlzWb58Oc8880yhyzFMmTKF66+//oxfO9u2bWPnzp1ceumlpw2fOHEiL730Env27KFdu3anSk0YE+hW7jjMqOkxbDlwgpu6NWJc/5ZUCQvehJ/NpyUbRORp4DogA/gduAOoCkwHGuF8OQxT1UP5LacoJRtyioyMZPny5dSpEzwncowJNinpmbz23Qbe+3EL4TUq8dLQdlwUFXz/834p2aCqTwJP5hicitP6N8aYYhez8wgPz4hh074kRkQ3YvyAFkHRRbMggu43T269SJ599llmzJhx2rBhw4bx6KOPllBUxpiiSs3IZOL3G3nnf1uoVy2Mj2+L5tJmdf0dVkAq89U5jTFlX9yuo4yaEcP6vccZ3rkBjw1sRXVr5Vt1TmNM2ZOWkcWbP2zkrUWbqVO1Ah+O7ELvFrleGmQ8WOI3xpRKqxOOMmpGLGsTjzGkY32eHNiaGpWtle8NS/zGmFIlPTOLfy3czD9/2EitKhV47+bOXN6qbNXL9zVL/MaYUmPdnmM8PD2G1QnHGNQhgqeual3q7ooVCCzxG2MCXkZmFv/+cQuvf7eBGpXK885NnejX5jx/h1VqWeI3xgS0jXuP8/CMGGJ3HWVgu3CeGdSG2tbKLxJL/MaYgJSRmcV7P23ltW83ULViOd66oSNXtgv3d1hlgiV+Y0zA2bQviVEzYli18wj925zH3we3oU7VsLPPaLxiid8YEzCyspQPft7Ky1+vp1KFUCaOuJCr2oWftXy6KRhL/MaYgLDj4ElGzYxh2dZDXNbyXJ4b0oZ61Sr6O6wyyRK/McavVJXJv+7guflrCRXhlWHtubZjfWvl+5AlfmOM3yQeTeaRmbH8tPEAF0fV4aWh7YioWcnfYZV5lviNMSVOVZnz+26enLuajEzl74Nac1O3xtbKLyGW+I0xJepAUirjZ8fxzZq9dG5ci1eGtSeyThV/hxVULPEbY0rMgvhExs+JJyklg3H9W3BHz6aEhlgrv6RZ4jfG+NzhE2k8NW81X6xKoE396rw6vAPNzq3m77CCliV+Y4xPfbN6D+PnxHPkZBr/1+cC7vtTFOVDQ/wdVlCzxG+M8QnPVn7L8Op8fFsXWkfU8HdYBkv8xphipqrMi03kmXmrOXIynQcuu4B7e0VRoZy18gOFJX5jTLHZfSSZx+bEsXD9fto3qMEnt3WlVUR1f4dlcrDEb4wpssws5eNftvHKN+sBeHxgK0b2iLQeOwHKEr8xpkjWJBxj3OxYYnYdpXfzuvx9cBsa1Krs77BMPizxG2MKJSU9kze+38h7P26hRqXyvHF9B65uH2FX35YClviNMQW2ZPNBxs2OZdvBkwzt1IBHB7S0e9+WIpb4jTFeO5GawQtfrePTpdtpfE5lJt/RlYui6vg7LFNAlviNMV75ZdMBHpkVy+4jydx+cRNG9W1OpQqh/g7LFIIlfmNMvpJSM3h+/lom/7qDJnWqMOOu7nSOrO3vsEwRWOI3xuRp8cYDjJkVS8LRZP7SswkP921OxfLWyi/tLPEbY85wPCWd5+avY8qyHTStU4WZd3enU2Nr5ZcVlviNMaf5aeN+xs6KI/FoMnde0pSHLm9mrfwyxhK/MQbIbuWvZcqynTStW4UZd/egU+Na/g7L+IAlfmMMv245yEPTY0g8msxdlzblwcuslV+WWeI3JoilZWTx+ncbePt/m2lcuzIz7+lBx0bWyi/rfJr4RaQmMAloAyhwG7AemAZEAtuA4ap62JdxGGPOtHl/Eg9MXUXc7qNc17khT1zViiph1hYMBr4ukP0GsEBVWwDtgbXAWOB7Vb0A+N59bYwpIarKlGU7GDhxMTsPn+Sdmzry4tB2lvSDiM/2tIhUBy4BRgKoahqQJiKDgF7uZB8Di4AxvorDGPOHQyfSGDMrlm/X7OXiqDq8Mqw959Wo6O+wTAnz5Vd8U2A/8KGItAdWAP8HnKuqiQCqmigi9XwYgzHG9eOG/Tw8I4ajJ9N57MqW3HZRE0KsXn5Q8mXiLwd0BP6mqr+KyBsU4LCOiNwJ3AnQqFEj30RoTBBISc/kpQXr+eDnrTQ7tyof3xptd8UKcr5M/LuAXar6q/t6Jk7i3ysi4W5rPxzYl9vMqvou8C5A586d1YdxGlNmbdqXxN+m/M7axGOM7BHJ2P4trJum8V3iV9U9IrJTRJqr6nqgD7DGfdwCvOD+/cJXMRgTrFSV6ct38tTcNVSuEMqHI7vQu4UdVTUOX5/G/xswWUQqAFuAW3F6Ek0XkduBHcAwH8dgTFA5lpLO+NlxfBmbyEVR5/Da8A7Uq24ncM0fzpr43d45dVV1c47h7VQ1Nr95VXUV0DmXUX0KEqQxxjsrdxzm/im/k3g0hdFXNOeeS8+3E7jmDPn24xeR4cA6YJaIrBaRLh6jP/JlYMYY72VlKf9atIlh7yxBFabf1Z2/9o6ypG9ydbYW/3igk3siNhr4VETGq+pswD5RxgSAfcdSeHD6Kn7edJAr24Xz3DVtqVGpvL/DMgHsbIk/1KPP/TIR6Q18KSINcEowGGP8aOH6fYyaHsOJtAxeGNKW67o0RMTaZCZ/Z0v8x0Xk/Ozj+27LvxfwOdDat6EZY/KSkp7Jy1+v5/3FW2lxXjWm3dCNqHrV/B2WKSXOlvjvIcchHVU9LiL9gOE+i8oYk6e4XUd5aPoqNu5L4ubujRk/oKX1zTcFkm/iV9WYPIanA5N9EpExJlfpmVn8a+Fm/vnDRs6pWoGPbu1Cr+bWN98UnFf9+EVkpKp+5ONYjDF52LQviYenryJm11Gubh/BM4NaU7NyBX+HZUopb/rxPwE0w7pvGlPisrKUj37ZxosL1lGpQihv3nAhA9tF+DssU8rlm/hF5F2gGnBjyYRjjMm26/BJRs+IZcmWg/ypRT1eGNLWrsA1xeJsLf4bgGhVzSqJYIwxTp2dGSt28cy8NaiqddM0xe5siX8gTl2dQTlLNhhjit/+46mMmx3Hd2v3Eh1Zm38Mb0/D2pX9HZYpY87Wq2eRiFwPfAZ0L5mQjAk+qsq82ESemruapNQMHh3QktsubkKolVwwPnDWk7uqGi8i15ZEMMYEoz1HU3js83i+W7uXdg1q8Mqw9jQ71y7GMr7jVXdOVU3wdSDGBBtVZdpvO3l2/lrSMrIYP6AFt13UhHKh+dZONKbIClSPX0QeAD5T1QO+CceY4LDz0EnGzo7l500HiW5SmxevbUeTOlX8HZYJEl4nfhFph3PXrHLAKz6LyJgyLDNL+fiXbbz89XpCQ4QJg9twQ3QjK59sSlRBWvy3A48Ad2CJ35gC27I/idEzY1mx/TC9mtfluWvaElGzkr/DMkHI25INYcAAoA3QS0QuVtXFPo3MmDIiK0v5eIlz9W1YuVBeHd6eay6sb/3yjd942+K/FvhaVVNF5EOc1r8lfmPOwvPq217N6/Lite04166+NX7mbeK/HRjlPp8PvC4iVVU1yTdhGVO6qSrTl+/k71+utatvTcDxpkhbTSBBVX8HUNVMEXkTiAZ+8G14xpQ++46lMHZ2HD+s20e3prV5eahdfWsCizcXcB0B/pxj2Gu+CsiY0mxuTAKPfx5PSnomT17Vilu6R1qPHRNwvD25GwpcCUR6zqOqr/omLGNKl0Mn0nj8i3j+G5tIh4Y1+cfw9pxft6q/wzImV94e458HpABxgFXqNMbDd2v2MnZ2HEeT0xh9RXPuuqSpXX1rApq3ib+BqrbzaSTGlDLHUtJ5Zt4aZq7YRcvw6nx6ezQtw6v7OyxjzsrbxP+ViPRV1W98Go0xpcTPmw4wekYMe46lcF/vKO7vcwEVylkr35QO3ib+pcAcEQkB0gEBVFWteWOCysm0DF74ah2fLNlO07pVmHVPDy5sVMvfYRlTIN4m/n/g1OOPU1X1YTzGBKwV2w/x8PQYth86ye0XN2H0Fc2pWD7U32EZU2DeJv6NQLwlfROMUtIzee27Dbz34xYialZiyl+60a3pOf4Oy5hC8zbxJwKLROQrIDV7oHXnNGXdqp1HGD0jho37khgR3YhHr2xJ1bACVTM3JuB4+wne6j4quA9jyrSU9Exe/24j7/64mXOrV+Tj26K5tFldf4dlTLHw9g5cT/s6EGMCxe87DjNqRgyb959gRHRDxg1oSfWK5f0dljHFxtsrd5vhFGmL5PQrd//km7CMKXkp6Zm89u0G3vtpC+dVr8gnt0VzibXyTRnk7aGeGcA7wCQg03fhGOMfK3ccZvSpVn4jxg9oQTVr5ZsyytvEn6Gqb/s0EmP8ICU9k1e/3cCkn7YQXqMSn94eTc8LrJVvyjava/WIyL3AHE7v1XPIJ1EZUwJWbD/M6JkxbNl/ghu6NmJcf2vlm+DgbeK/xf072mOYAk3PNqNb2XM5sFtVB4pIbWAazvmCbcBwVT3sbcDGFFVKeib/+GY9kxZvJaJGJT67vSsXX1DH32EZU2K87dXTpAjr+D9gLZBd3mEs8L2qviAiY93XY4qwfGO8tmL7IUbPiGXLgRPc2LUR4wZYv3wTfM5aVUpEKotI+xzDGolIfS/mbYBTx3+Sx+BBwMfu84+BwV5Ha0whJadlMuHLNQx9ZwmpGVlMvqMrz17T1pK+CUrefOrTgdki0k5VT7jDJgHjgd1nmfd14BGgmsewc1U1EUBVE0WkXm4zisidwJ0AjRo18iJMY3K3fNshRs+MZeuBE9zUrRFj+1sr3wS3s7b4VTUd56TudeC09oG6qro8v/lEZCCwT1VXFCYwVX1XVTuraue6da2XhSm45LRM/v7lGob9ewnpmVn8546uTBhsrXxjvP0PmAS8B3wA3Ax86MU8FwFXi8gAoCJQXUQ+A/aKSLjb2g8H9hUibmPy9du2QzzitvL/3K0xY/u3oIolfGMAL1r8AKq6Dk5dwTsC+NSLecapagNVjQSuB35Q1ZuAufzRS+gW4ItCxG1MrpLTMnlm3hqGZ7fy/9KVvw9uY0nfGA8F+W94H6flH1vE7pcvANNF5HZgBzCsCMsy5pRlWw8xemYM2w+etFa+MfkoyH/FdOAN4JmCrkRVFwGL3OcHgT4FXYYxeTmZlsFLC9bz8ZJtNKjl1Mvvfr7VyzcmL14nflU9CdTwYSzGFNjSLQd5ZGYsOw6dZGSPSB7p15zKFayVb0x+7D/ElEonUjN4acE6Pl6yncbnVGbqnXZXLGO8ZYnflDq/bD7AmFmx7DqczK0XRTL6CmvlG1MQ9t9iSo0TqRm88NU6Pl26nchzKjPtzu5EN6nt77CMKXUs8ZtS4ZdNB3hkViy7jyRz+8VNGNW3OZUqhPo7LGNKJUv8JqAlpWbw/Py1TP51B03qVGHGXd3pHGmtfGOKwhK/CViLNzrH8hOOJvOXnk14uG9zKpa3Vr4xRWWJ3wSc4ynpPDd/HVOW7aBp3SrMvLs7nRpbK9+Y4mKJ3wSUXzYfYPSMWBKPJnPXJU158PJm1so3pphZ4jcBISU9k1e+du6K1aROFWbe04OOjWr5OyxjyiRL/MbvVicc5cFpq9iwN4k/d2vMuAEtrF++MT5k/13GbzKzlH//uJnXvt1ArcoV+OjWLvRqnut9eYwxxcgSv/GLHQdP8tD0VSzffpgr24YzYXAbalWp4O+wjAkKlvhNiVJVpi/fyTPz1hASIrx+XQcGdYhARPwdmjFBwxK/KTH7j6cybnYs363dR4/zz+GVYe2JqFnJ32EZE3Qs8ZsSsSA+kUfnxHM8NYPHB7bi1h6RhIRYK98Yf7DEb3zq0Ik0npy7mnkxCbSpX50pwzvQ7Nxq/g7LmKBmid/4zIL4PTz2eRxHk9N5+PJm3N3rfMqHenWbZ2OMD1niN8XusNvKnxuTQOuI6nx6e1dahlf3d1jGGJclflOsvl69h0fnxHPkZBoPXtaMe3tbK9+YQGOJ3xSLIyfTeGruaj5flUCr8Op8cls0rSKslW9MILLEb4rs2zV7GT8njsMn0njgsgv4a+8oa+UbE8As8ZtCO3TCaeXPjUmgZXh1Prq1C60javg7LGPMWVjiN4UyPy6RJ76I52hyOg9e1ox7ep1PhXLWyjemNLDEbwpk//FUnpwbz/y4PbStX4PP7uhKi/PsWL4xpYklfuMVVeWLVQk8PW81J1IzeaRfc+7s2ZRydizfmFLHEr85q4QjyTw6J46F6/dzYaOavDy0HVH17OpbY0orS/wmT1lZyuRlO3jxq3VkZilPDGzFLT0iCbUaO8aUapb4Ta42709i3Kw4lm07RM8L6vDcNW1pWLuyv8MyxhQDS/zmNOmZWbz30xZe/24jFcuF8PLQdgzt1MDq5RtThljiN6fE7z7KIzNjWZN4jP5tzuPpQa2pV62iv8MyxhQzS/yG9Mws/vnDJt5auInaVSrwzk0d6dcm3N9hGWN8xBJ/kNuw9zgPTV9F/O5jDLmwPk9e1Zoalcv7OyxjjA9Z4g9SmVnK+4u38Mo3G6gaVs5a+cYEEUv8QWjHwZOMmhHDsm2HuLzVuTx3TVvqVgvzd1jGmBJiiT+IqCpTlu1kwn/XECrCP4a1Z0jH+tZjx5gg47PELyINgU+A84As4F1VfUNEagPTgEhgGzBcVQ/7Kg7j2HsshTGzYlm0fj8XRZ3DS0PbU79mJX+HZYzxA1+2+DOAh1V1pYhUA1aIyLfASOB7VX1BRMYCY4ExPowj6M2LSeCxz+NJzcjkqatacXP3SELs6ltjgpbPEr+qJgKJ7vPjIrIWqA8MAnq5k30MLMISv08cOZnG41+sZl5MAu0b1uTV4e05v25Vf4dljPGzEjnGLyKRwIXAr8C57pcCqpooIvXymOdO4E6ARo0alUSYZcqC+D088UU8h06k8fDlTr18q6RpjIESSPwiUhWYBTygqse8PZGoqu8C7wJ07txZfRdh2bL3WApPfrGaBav30DK8Oh+M7EKb+nZXLGPMH3ya+EWkPE7Sn6yqs93Be0Uk3G3thwP7fBlDsFBVpv62k+fmryU1I4tH+jXnLz2b2r1vjTFn8GWvHgHeB9aq6qseo+YCtwAvuH+/8FUMwWLL/iTGzY7j162H6Na0Ns8PaUeTOlX8HZYxJkD5ssV/EfBnIE5EVrnDxuMk/OkicjuwAxjmwxjKtPTMLN79cQtvfL+RsHIhvDCkLdd1aWj98o0x+fJlr57FQF4ZqI+v1hssYncdYcysONZmV9K8ujX1qlslTWPM2dmVu6XMybQMXvt2A+8v3kqdqmG8c1Mn+rU5z99hGWNKEUv8pchPG/czfk4cOw8lc0PXRozp14IalaySpjGmYCzxlwKHT6Qx4b9rmbVyF03rVGHand3o2vQcf4dljCmlLPEHMFVlXmwiT89dzdHkdO7rHcV9f4qiYvlQf4dmjCnFLPEHqN1Hknn883h+WLeP9g1q8NkdXWkZXt3fYRljygBL/AEmK0v5dOl2XlqwjiyFxwe2YmSPSEKtqJoxpphY4g8gG/YeZ+ysWFbuOMIlzery7OA2NKxd2d9hGWPKGEv8ASA1I5N/LdzMvxZtompYOV67rj2DO9gNUowxvmGJ389WbD/M2FmxbNyXxOAOETw+sBXnVLXbIBpjfMcSv58kpWbw8oJ1fLJ0OxE1KvHhrV3o3TzXCtXGGFOsLPH7wQ/r9vLonHj2HEthZI9IRvVtTpUw2xXGmJJh2aYEHUhK5el5a5gXk0Dzc6vxrxs7cmGjWv4OyxgTZCzxlwBVZfbK3fz9v2s4mZrJw5c3465Lz6dCOauVb4wpeZb4fSzxaDLjZ8excP1+OjeuxQvXtiOqnt331hjjP5b4fURVmb58JxO+XEtGlvLkVa24pXskIXYhljHGzyzx+8DuI8mMnRXLTxsP0K1pbV68th2Nz7E7YhljAoMl/mKkqvxn2Q6en78OVeXvg9twY3Qja+UbYwKKJf5isvPQScbMiuWXzQe5OKoOzw9pa+UWjDEByRJ/EWVlKZ/9up0XvlpHiAjPD2nL9XbfW2NMALPEXwTbDpxgzKxYft16iEub1eX5IW2JqFnJ32EZY0y+LPEXQmaW8tEv23j563WUDw3h5aHtGNqpgbXyjTGlgiX+AtqyP4lHZsayfPth+rSox7PXtOW8GhX9HZYxxnjNEr+XMrOUDxZv5ZVv1lOxfCivDm/PNRda6WRjTOljid8Lm/YdZ/TMWH7fcYS+rc5lwuA21KturXxjTOlkiT8fGZlZvPfTVl77bgNVKoQyccSFXNUu3Fr5xphSzRJ/HlYnHGXc7Dhidx2lf5vzeGZQG+pWsxukGGNKP0v8OZxMy+D17zby/uKt1Kpcgbdu6MiV7cL9HZYxxhQbS/weFq7fx2Nz4tl9JJkR0Q0Z268lNSqX93dYxhhTrCzxA/uOp/DMvDV8GZtIVL2qTL+rO9FNavs7LGOM8YmgTvwZmVlM/nUH//hmPSnpWTx0eTPuurQpYeVC/R2aMcb4TNAm/qVbDvLU3NWs23Oci6LO4ZlBbTi/rt0gxRhT9gVd4k88msxz89cxLyaB+jUr8faNHenX5jzrommMCRpBk/iPpaTz7v+2MGnxFrIU7u9zAfdcej6VKthhHWNMcCnziT81I5PPlu7gzR82cvhkOle1j+CRK5pbrXxjTNAqs4k/KTWD//y6nUk/bWXf8VQujqrDmH4taNughr9DM8YYv/JL4heRfsAbQCgwSVVfKK5l7zh4kqm/7WDyrzs4mpxOj/PP4bXrOnBRVJ3iWoUxxpRqJZ74RSQUeAu4HNgF/CYic1V1TWGXeSAplR/W7eOLVbv5edNBQgQub3Uu9/SKokPDmsUUuTHGlA3+aPFHA5tUdQuAiEwFBgH5Jn5VJS0zi2PJGew+ksy2AyeI3XWUFTsOE7vrCKpQv2YlHr68GcM6N7Qa+cYYkwd/JP76wE6P17uArvnNEL/7KE3GzT9jeMXyIbSOqMEDfZrRp2U9WkdUt26ZxhhzFv5I/LllZj1jIpE7gTvdl6nbXxwYn9vC1gOziy+2wqoDHPB3EPmw+IomkOML5NjA4iuqosbXOLeB/kj8u4CGHq8bAAk5J1LVd4F3AURkuap2LpnwCs7iKxqLr/ACOTaw+IrKV/GFFPcCvfAbcIGINBGRCsD1wFw/xGGMMUGpxFv8qpohIvcBX+N05/xAVVeXdBzGGBOs/NKPX1XnA2eerc3bu76KpZhYfEVj8RVeIMcGFl9R+SQ+UT3jvKoxxpgyzB/H+I0xxvhRQCd+EeknIutFZJOIjPVTDA1FZKGIrBWR1SLyf+7wp0Rkt4isch8DPOYZ58a8XkSuKIEYt4lInBvHcndYbRH5VkQ2un9r+SM+EWnusY1WicgxEXnAn9tPRD4QkX0iEu8xrMDbS0Q6udt9k4hMlGK6iCSP+F4WkXUiEisic0Skpjs8UkSSPbbjO36Kr8D7s4Tjm+YR2zYRWeUOL9Htl08+KdnPn6oG5APnxO9moClQAYgBWvkhjnCgo/u8GrABaAU8BYzKZfpWbqxhQBP3PYT6OMZtQJ0cw14CxrrPxwIv+iu+HPt0D07fYr9tP+ASoCMQX5TtBSwDuuNcm/IV0N+H8fUFyrnPX/SIL9JzuhzLKcn4Crw/SzK+HOP/ATzhj+1H3vmkRD9/gdziP1XaQVXTgOzSDiVKVRNVdaX7/DiwFufq47wMAqaqaqqqbgU24byXkjYI+Nh9/jEw2GO4v+LrA2xW1e35TOPz+FT1R+BQLuv1enuJSDhQXVWXqPNf+InHPMUen6p+o6oZ7sulONe/5Kmk48tHQGy/bG6reDgwJb9l+Cq+fPJJiX7+Ajnx51baIb+E63MiEglcCPzqDrrP/en9gcdPM3/ErcA3IrJCnCueAc5V1URwPmxAPT/Gl+16Tv+HC5TtBwXfXvXd5yUdJ8BtOC28bE1E5HcR+Z+I9HSH+SO+guxPf22/nsBeVd3oMcwv2y9HPinRz18gJ36vSjuUFBGpCswCHlDVY8DbwPlAByAR5+cj+Cfui1S1I9Af+KuIXJLPtH7ZruJcrHc1MMMdFEjbLz95xeOv7fgokAFMdgclAo1U9ULgIeA/IlLdD/EVdH/6az+P4PTGh1+2Xy75JM9J84ijSPEFcuL3qrRDSRCR8jg7abKqzgZQ1b2qmqmqWcB7/HE4osTjVtUE9+8+YI4by17352D2z9Z9/orP1R9Yqap73VgDZvu5Crq9dnH64RafxykitwADgRvdn/e4hwAOus9X4BwDblbS8RVif/pj+5UDhgDTPOIu8e2XWz6hhD9/gZz4A6K0g3tM8H1graq+6jE83GOya4DsHgRzgetFJExEmgAX4JyE8VV8VUSkWvZznJOA8W4ct7iT3QJ84Y/4PJzW0gqU7eehQNvL/Tl+XES6uZ+Rmz3mKXbi3LxoDHC1qp70GF5XnHtcICJN3fi2+CG+Au3Pko7PdRmwTlVPHSIp6e2XVz6hpD9/RT1L7csHMADnrPdm4FE/xXAxzk+oWGCV+xgAfArEucPnAuEe8zzqxryeYuqpkE98TXHO+scAq7O3E3AO8D2w0f1b2x/xueurDBwEangM89v2w/kCSgTScVpOtxdmewGdcRLcZuBN3AsifRTfJpxjvdmfwXfcaa9193sMsBK4yk/xFXh/lmR87vCPgLtzTFui24+880mJfv7syl1jjAkygXyoxxhjjA9Y4jfGmCBjid8YY4KMJX5jjAkylviNMSbIWOI3xpggY4nfBAURqSki93q8jhCRmT5a12AReaKQ837nWZLXGF+wfvwmKLgFsb5U1TYlsK5fcK6wPVCIeW8BGqjqs8UfmTEOa/GbYPECcL44N9t4WZwbcMQDiMhIEflcROaJyFYRuU9EHnIrNi4VkdrudOeLyAK3CupPItIi50pEpBmQmp30ReQj9yYZv4jIFhEZ6g4PF5Ef3XjiPapCzsUpb2GMz1jiN8FiLM69ADqo6uhcxrcBbsApLvYscFKdio1LcOqggHPj67+paidgFPCvXJZzEc6l/57CcS7VH4jzBYS7rq9VtQPQHufSfVT1MBAmIucU4j0a45Vy/g7AmACxUJ0bYxwXkaPAPHd4HNDOLaPbA5ghf9zhLiyX5YQD+3MM+1ydqpVrRORcd9hvwAdupcbPVXWVx/T7gAic+kbGFDtr8RvjSPV4nuXxOgungRQCHHF/MWQ/WuaynGSgYj7LFjh1l6hLgN3ApyJys8c0Fd3lGOMTlvhNsDiOc4/TQlHnZhlbRWQYOOV1RaR9LpOuBaLOtjwRaQzsU9X3cMr0dsxeLnAezn2UjfEJS/wmKKhzs42f3ROpLxdyMTcCt4tIdgns3O4B/SNwoXgcD8pDL2CViPyOUxr4DXd4J2Cp/nF/XWOKnXXnNKaYicgbwDxV/a6Q885V1e+LPzJjHNbiN6b4PYdz85nCiLekb3zNWvzGGBNkrMVvjDFBxhK/McYEGUv8xhgTZCzxG2NMkLHEb4wxQeb/AXhsgun7JlApAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "selection = df[['time', 'displacement_sq']].groupby(\"time\").mean()\n",
    "\n",
    "lm_time = np.arange(0, 2000)\n",
    "\n",
    "# line_y = np.linspace(0, np.max(selection[\"displacement_sq\"]), 2000)\n",
    "\n",
    "lm = np.poly1d(np.polyfit(df[\"time\"], df[\"displacement_sq\"], 1))\n",
    "   \n",
    "def predict_line(m, b, n):\n",
    "    return [m*i + b for i in range(n)]\n",
    "\n",
    "line_y = predict_line(m=lm.coef[0], b=0, n=len(lm_time))   \n",
    "\n",
    "d = calc_d_from_step(step=lm.coef[0])\n",
    "d_text = f\"meas D = {d:.2e}\"\n",
    "exp_text = f\"exp D = {convert_d_from_cm2_s_to_nm2_ns(d=1.88e-9):.2e}\" # 1.88e-9 cm2/s\n",
    "    \n",
    "step_text = f\"step_nm = {round(lm.coef[0], 5)} nm\"\n",
    "\n",
    "print(d_text)\n",
    "print(exp_text)\n",
    "print(step_text)\n",
    "\n",
    "# plot the displacement_sq against the fit\n",
    "\n",
    "plt.suptitle(\"displacement of LHCII from origin\")\n",
    "plt.title(\"D in nm^2/ns units\")\n",
    "plt.xlabel(\"time (ns)\")\n",
    "plt.ylabel(\"<nm>^2\")\n",
    "plt.xlim(0, np.max(df[\"time\"]) * 1.05)\n",
    "plt.ylim(0, np.max(selection[\"displacement_sq\"]) * 1.05)\n",
    "# plt.annotate(text=d_text, xy=(20, np.max(selection[\"displacement_sq\"]) * 0.95))\n",
    "plt.annotate(text=step_text, xy=(20, np.max(selection[\"displacement_sq\"]) * 0.85))\n",
    "plt.annotate(text=exp_text, xy=(20, np.max(selection[\"displacement_sq\"]) * 0.90))\n",
    "plt.plot(selection[\"displacement_sq\"])\n",
    "# plt.plot(lm_time, line_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "meas D = 4.33e-11\n",
      "exp D = 2.00e-06\n",
      "step_rads = 1e-05 rads/ns\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x213840800a0>]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEjCAYAAAA41BqSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAz8klEQVR4nO3deZgU5bXH8e+ZQSABUdkURWQTUWEGEBBFEdQoiLhECBCMgCiuuCWKxrglmFVNxHhzL3HDiErclyiJGhAXVBZhZIvKpghBAUFAtpk+94+qbnpmunt6lp4Z2t/nefqZ7rffeutUdc/p6uqqU+buiIhI9smp6QBERCQzlOBFRLKUEryISJZSghcRyVJK8CIiWUoJXkQkSynBZyEza21mbmZ1ajqWdJjZz83sgRqY77lm9rmZbTWzrtU9//LYm2IFCONsW9V9pXyU4KVSzOx2M3usHP37mtnq+DZ3/7W7X1T10ZXpLuBKd2/o7h+WfDL8kGyfoH2Umb2doH2lmZ0a97iFmT1oZmvNbIuZLTWzO8ysQcnx01iPKWOtbcI4l1d1XykfJXhJaW/5FlBBhwGLMjGwmTUGZgHfA45z932BHwD7A+0qMGTSWGvTa1SbYhEl+BplZjea2bJw626xmZ0b99woM3vbzO4ys6/NbIWZDYh7vo2ZzQynfd3M7k+2BWhm+8VtSX5hZhPMLDdJ39vN7Gkze8zMvgFGmdnBZvaimW00s0/N7OKwb3/g58DQ8Gv2grB9tJktCWNbbmaXhO0NgFeBg8P+W8Oxi229mtlZZrbIzDaZ2QwzOzLuuZVm9jMzKzCzzWY21czqJ1mWHDP7hZmtMrMvzezRcF3UM7OtQC6wwMyWpfmSlcd1wBbgfHdfCeDun7v71e5ekO4gyWIN18N4MysAtplZnTTW2/XhetsWvh8ONLNX495DB6SI4+Lwtd8YvhcOjnvOzewKM/sE+CSuLfrtpImZvWRm35jZ7PD993aJ6aN9Hwnfy/8I43rfzCrygSgowde0ZcCJwH7AHcBjZtYi7vljgf8ATYHfAw+amYXPPQ58ADQBbgd+kmI+k4FCoD3QFTgNSLVL5GzgaYKtzSnAE8Bq4GBgMPBrMzvF3acBvwamhl+z88PpvwTOBBoBo4E/mlk3d98GDADWhP0buvua+BmbWYdwftcAzYBXgJfMrG5ctx8B/YE2QB4wKslyjApv/YC2QEPgz+6+090bhn3y3T0TCeRU4Fl3j1RmkDJiHQ4MJHid2lL2ejuP4FtEB2AQwYftzwneXznAVYliMLOTgd8QrPcWwCrgyRLdziF4vx6VYIj7gW3AQcDI8JbKcIL/hwOAT4E7y+gvSSjB1yB3f8rd17h7xN2nEmz99Izrssrd/+ruRQRJugVwoJm1AnoAt7r7Lnd/G3gx0TzM7ECCpHqNu29z9y+BPwLDUoQ2y92fD5NTU+AEYLy773D3+cADpPhAcfd/uPsyD7wJ/IvggywdQ4F/uPtr7r6bYN/z94Dj4/pMDNfbRuAloEuSsUYA97j7cnffCtwEDLOq2Y3QK9xSjt2AVnHPNwHWVsF8UpkYfivYTnrr7T53X+fuXwBvAe+7+4fuvhN4juDDP5ERwEPuPi/sexNwnJm1juvzG3ffGMYSE35TPA+4zd2/dffFBO/lVJ519w/cvZBgA6NLmWtCElKCr0FmdoGZzY9LEJ0IEmrUf6N33P3b8G5Dgi3pjXFtAJ8nmc1hwD7A2rj5/B/QPEVo8WNF57Ulrm0VcEiK5RpgZu+FX+c3AWeUWK5UDg7HByD8kPm8xPz+G3f/W4J1UuZY4f06wIFpxpLKe+6+f/wN+Czu+Q0EH8iZVPJ1Kmu9rYu7vz3B47TWY/hhuaHE2Mnef80I1vnnafSNSvf1lTIowdcQMzsM+CtwJdAkTBALAUs1XWgt0NjMvh/XdmiSvp8DO4GmccmokbsfnWL8+BKja8J57RvX1gr4IkFfzKwe8AzBFuSB4XK9wp7lKqt86RqCD6XoeEawbF8knSLNscK4Cyme2DLldeBcM8vk/1jJ16mq1ltJJcduQPANJX7sZK/rVwTrvGVcW7L3qlQxJfia04Dgn+IrCH6YJNiCL5O7rwLmALebWV0zO45gn2qivmsJdpHcbWaNwh8e25nZSWnO63PgXeA3ZlbfzPKAMQRfnSFIlq3jElldoF64XIUW/DB8WtyQ64AmZrZfkln+HRhoZqeY2T7ATwk+oN5NJ94SngCuteAH6Ybs+b2gsBxj1A2XO3pL+ON0AvcQ/AYxOfwwx8wOMbN7wnVY1apyvZX0ODDazLqEH+C/Jti9s7KsCcPdi88SvFe/b2YdgQuqICZJgxJ8DQn3Rd5NcCjdOqAz8E45hhgBHEfwVXkCMJXgHzqRCwgS72Lga4IfUMuz+2A40JpgS+45gv2pr4XPPRX+3WBm88JdOVcRJJyvgR8T9/uAuy8lSLzLw11GsaMxwuf/A5wP3AesJ/jgGuTuu8oRb9RDwN+AmcAKYAcwrpxjLCLYfRG9jU5novD3geOB3cD7ZrYFeAPYTPDDYZWq4vVWcuw3gFsIvpmtJTjMM9VvOCVdSXAgwX8JXo8nSP5elSpkuuBHdjCzqcBSd7+tpmMRScXMfgcc5O5lHU0jlaQt+L2UmfUId7XkWHA8+tnA8zUclkgpZtbRzPIs0JNgF99zNR3Xd4HOOtt7HUSwb7MJwTHql+0Np7DLd9K+BLtlDiY4R+Ju4IUajeg7QrtoRESylHbRSK0QnjJfI/tkbS+rvpmIqSKjJKAEL1XOgron28NaIpvM7F0zuzTVMeHuPsDdyzrDsVpYUP/nE9tTAbLUYX1m9mMze7wm4kskviJjWM9lQk3HJDVPCV4yZVBYQfEw4LfAeODBTM+0irbCtxEcZrgfQd2Ue83s+BJ9ziA4gUuk1lKCl4xy983u/iJBrZSRZpbwZC4Lqh9eFN5PWUkzwbSJKiumqtSZG4693syWExTsio/5NndfGtYIep+gbstxcdPnEBTtmha3e2ekmX0WjnlzXN+eZjbHgkqK68zsniTLUKrGvJWjymK0r5mNJThH4oZwt81L4fPjLagkusXM/mNmpyRbn5I9lOClWrj7BwRH+6RbdCxVJc1EYpUVwzNVU1XqvJig2mVXoDtBhcyEzOx7BIXd4mux9wSWu/v6uLYTgCOAU4BbbU+p3nuBe929EcEJQn8vY7lTKbPKortPIjjL+PfhbptBZnYEwclGPcJvVacDKysRh+wllOClOq0BGqfZN2ElzRT94ysrllWp80fAn8L+GwlK4Sbzv8AC4J9xbQMpvXvmDnff7u4Lwv7R0sm7gfZm1tTdt7r7e2UueXIVrbJYRFA+4igz28fdV7p7JmrgSy2jBC/V6RBgY5p9k1XSTKZYhUJLXanz4BL94ytOxo/xh3C6H3nx44kT7X9PVgFxDEH99aUWXOzizBTLUJYKVVl0908J6sTfDnxpZk+WLBEh2UkJXqqFmfUgSPClrmVaRWIJ2Mqu1LmW4hUN4+u4R8e4g6CO/mnu/k1c+0EE3ybmpRWU+yfuPpygPPPvgKctvCZrCduAWHXQcD4VVerkFnd/3N1PIPjR28NYJMspwUtGWVDB8kyCKwA95u4fVcNsy6rU+XfgKjNracFl6m4sEfNNBEXSfuDuG0qMfQYwrcQWfVJmdr6ZNQvrs28Km4sSdF0AHG1Bxcb6BFvbFbWO4ApP0RiOMLOTLagEuYOgaFqiGCTLKMFLprxkQQXFz4GbCcrnplWJsbLSqNT5V4J96gsItsSfLTHErwm26j+xPdeO/Xn4XHkPj+wPLLLguqr3AsPcfUeCmD8GfklQR/4TKvdN50GC/e2bzOx5gv3vvyWoMvlfgm8TP08+uWQLlSoQSVN4jP1/gXbuvrmm4xEpi7bgRdLXGLhFyV32FtqCFxHJUtqCFxHJUrWqel7Tpk29devWNR2GiMheY+7cuevdvVmi52pVgm/dujVz5syp6TBERPYaZpbwRD3QLhoRkaylBC8ikqWU4EVEspQSvIhIllKCFxHJUkrwIiJZSgleRCRLKcGLiJSTu1MUcSpb6iUScb7dVUhRJDMlY2rViU4bt+3ib7NWBg/MYldnMIPoIzMwd8Axi35CBSvHwjaAHAezsN2cHCd2f09/wyMRCiNOJFKExeYFORaOb4a5k2NGjgXzyInOw5wcLBg/1hZOi5GDQxhjTmzcsB/BmGYe+5TNsSA+j0BRJEJhpAiPgBPE7B4se1BaPFjuOjmQm5NDbg7skxPEmBudV9iWmxOuvbj1FruPxe5Hxwzaw/te4n4wcThNfHv0PiX6E9c3cQx7XsLk8w0nCF+LuPZi45ZsT3zfPYK7EwE8Etx3D9oj7mHXYD1HwrGDacDDf2on/JvwFs6DYByi07rjROLGD5axTk6weBH32PjR+5HwtS6KOJFIhKKIUxSJFL8fJpvCoqJg+nB5zYP7Fi67l7F+rFgfirXn5ASvq4Xvw+j7JHi9Y1Nhsfcm4WsYF0P4+rrveT9HY4jed4/O27Ho/1H4P2Lh/9+eC/M6kUjwXiqKRLDwNY3OF0r339NOsT7R/0uiryF78kz8utrzOsZNG/YI/u/DsW3PfC1uXe1ZB2EuCtf/foMmcGiPQVS1WpXgD9i8mGHTuiV8cXJMRdGkakT/Cff2r68Rj30EBv8pBvH/NfEf1XvSWbQtWZ+wPRwr2j9+2nCzoFh7yflGt7QSzSs6bWxrrNj9uHjciH5EeJLlMjPcgw0bN8Os9LKAhctTMs648TxcFgOznCBKA8gpFqeZhRtmQcxBsk8QZ1y72Z71YBabCRBsQObk5LBvw0ZkQq1K8N6gGbt7XVzsMzfZC1XyjRO7bxa+MYh7kyUYJ/o2y8kJV3L4Ivme8YMXfs8/jzvBFgJGxKPjhH0s2lb8DRN9LrphEszDS40fibsffRPl5Bg5llP8TRHet3D6Iocid4oiwVZekUPEIRIuS2EkQmEkbv3Er6ti/9Bx/2Ce6h8i/p8s0WthsbE8+h8DROIThlvc/OKTR+I+8Ykk4uEyR5zCCBRGnMJIkDbMordgneUEX4kwywm/HcXfzwkfQ074HrC4W3T9Y9HniL0WZjnhaxOMH70fTBf2zwluObGEkBMbI5gfYDm4w66i4FtknZwccnOM3Do55FgOdXKC+eTm5lAnx6ibW4d99smhbm4OdesEt33CaUQSqVUJPme/Q/j+gF/WdBgiIllhb/+WmhF9+/bliCOOIC8vj44dO3LllVeyadOmKht/9uzZ5Obm8vTTTyd8fsyYMeTn55OXl8fgwYPZunUrAC+88AJ5eXl06dKF7t278/bb5buq24YNG+jXrx8NGzbkyiuvLPbczTffzKGHHkrDhg0rtlAiUusowScxZcoUCgoKKCgooF69epx99tlVMm5RURHjx4/n9NNPT9rnj3/8IwsWLKCgoIBWrVrx5z//GYBTTjmFBQsWMH/+fB566CEuuuiics27fv36/OpXv+Kuu+4q9dygQYP44IMPyrcwIlKr1doEv3LlSjp27MhFF11Ep06dGDFiBK+//jq9e/fm8MMPjyWjbdu2ceGFF9KjRw+6du3KCy+8EJv+xBNPpFu3bnTr1o13330XgLVr19KnTx+6dOlCp06deOutt1LGUbduXX7/+9/z2WefsWDBgkov13333cd5551H8+bNk/Zp1Cj4wcXd2b59e+yHo4YNG8bub9u2Le4HJXjsscfo2bMnXbp04ZJLLqGoqKjUuA0aNOCEE06gfv36pZ7r1asXLVq0qNSyiUjtUmsTPMCnn37K1VdfTUFBAUuXLuXxxx/n7bff5q677uLXv/41AHfeeScnn3wys2fPZvr06Vx//fVs27aN5s2b89prrzFv3jymTp3KVVddBcDjjz/O6aefzvz581mwYAFdunQpM47c3Fzy8/NZunRpqeeGDh1Kly5dSt0effTRUn2/+OILnnvuOS699NIy5zl69GgOOuggli5dyrhx42Ltzz33HB07dmTgwIE89NBDACxZsoSpU6fyzjvvMH/+fHJzc5kyZUqZ8xCR7FarfmQtqU2bNnTu3BmAo48+mlNOOQUzo3PnzqxcuRKAf/3rX7z44oux3Q47duzgs88+4+CDD+bKK6+MJbyPP/4YgB49enDhhReye/duzjnnnLQSPESP3S1t6tSpaS/PNddcw+9+9ztyc3PL7Pvwww9TVFTEuHHjmDp1KqNHjwbg3HPP5dxzz2XmzJnccsstvP7667zxxhvMnTuXHj16ALB9+/aU3xBE5LuhVif4evXqxe7n5OTEHufk5FBYWAgEifeZZ57hiCOOKDbt7bffzoEHHsiCBQuIRCKx3RJ9+vRh5syZ/OMf/+AnP/kJ119/PRdccEHKOIqKivjoo4848sgjSz03dOhQ/vOf/5Rqv+6660qNO2fOHIYNGwbA+vXreeWVV6hTpw7nnHNOwvnm5uYydOhQ/vCHP8QSfFSfPn1YtmwZ69evx90ZOXIkv/nNb4r1ee6557jjjjsAeOCBB+jevXvK5RSR7FKrE3w6Tj/9dO677z7uu+8+zIwPP/yQrl27snnzZlq2bElOTg6TJ0+O7ZNetWoVhxxyCBdffDHbtm1j3rx5KRP87t27Y0eY5OXllXq+PFvwK1asiN0fNWoUZ555Zqnk7u4sW7aM9u3b4+689NJLdOzYEQh2WbVr1w4zY968eezatYsmTZpwyimncPbZZ3PttdfSvHlzNm7cyJYtW2Jb+yLy3VSr98Gn45ZbbmH37t3k5eXRqVMnbrnlFgAuv/xyJk+eTK9evfj4449p0KABADNmzKBLly507dqVZ555hquvvjrhuCNGjIiNuW3bttiPt5lyxhlnsGbNmtjWeOfOnencuTNr167l1ltvBeCZZ56hU6dOdOnShSuuuIKpU6diZhx11FFMmDCB0047jby8PH7wgx+wdu3ahPNp3bo11113HY888ggtW7Zk8eLFANxwww20bNmSb7/9lpYtW3L77bdndHlFJPOsssVyqlL37t1dF90WEUmfmc1194T7X/f6LXgREUnsO5fgM3WW6pQpU8jLyyMvL4/jjz8+6THzK1as4Nhjj+Xwww9n6NCh7Nq1Cwj2vV911VW0b9+evLw85s2bV+4YfvOb39C+fXuOOOII/vnPf8bad+3axdixY+nQoQMdO3bkmWeeqdhCisheJaMJ3sz2N7OnzWypmS0xs+MyOb90ZeIs1TZt2vDmm29SUFDALbfcwtixYxP2Gz9+PNdeey2ffPIJBxxwAA8++CAAr776Kp988gmffPIJkyZN4rLLLivX/BcvXsyTTz7JokWLmDZtGpdffnnsh+U777yT5s2b8/HHH7N48WJOOumkyi2siOwVMr0Ffy8wzd07AvnAkvJMnOjszNmzZ5OXl8eOHTvYtm0bRx99NAsXLmTGjBn06dOHc889l6OOOopLL72USCSScvyqPEv1+OOP54ADDgCCs0JXr15dqo+78+9//5vBgwcDMHLkSJ5//nkgqDNzwQUXYGb06tWLTZs2xX4oTecs1RdeeIFhw4ZRr1492rRpQ/v27WNn+z700EPcdNNNQHCIadOmTSu1rCKyd8hYgjezRkAf4EEAd9/l7pvSnT7Z2Zk9evTgrLPO4he/+AU33HAD559/Pp06dQLggw8+4O677+ajjz5i2bJlPPvss2XOp6rOUo334IMPMmDAgFLtGzZsYP/996dOneDo1JYtW/LFF18AwVmuhx56aKxv9Ll0z1JNNn1099Mtt9xCt27dGDJkCOvWrStzvYjI3i+Tx8G3Bb4CHjazfGAucLW7b4vvZGZjgbEArVq1irWnOjvz1ltvpUePHtSvX5+JEyfGpunZsydt27YFYPjw4bz99tuxreVUquIs1ajp06fz4IMPJqz0mGg+0XoyyZ5L9yzVZNMXFhayevVqevfuzT333MM999zDz372M/72t7+Ve9lEZO+SyQRfB+gGjHP3983sXuBG4Jb4Tu4+CZgEwWGSce0Jz84E2LhxI1u3bmX37t3s2LEjdox7fPGtRI8TqaqzVAEKCgq46KKLePXVV2nSpEmp55s2bcqmTZsoLCykTp06rF69moMPPhgItrg///zzWN/oc+mepZps+iZNmvD9738/dsLTkCFDYvv9RSTLJb+uZOVuwEHAyrjHJwL/SDXNMccc41GLFi3y9u3b+7p169zdfcOGDb5y5Up3dx80aJBPmTLFJ0yY4FdccYW7u0+fPt3r16/vy5cv96KiIj/ttNP86aef9pJOOukknz17tru779q1y6+//nrv06dPqX7ltWrVKm/Xrp2/8847KfsNHjzYn3jiCXd3v+SSS/z+++93d/eXX37Z+/fv75FIxGfNmuU9evQocz3EW7hwoefl5fmOHTt8+fLl3qZNGy8sLHR396FDh/obb7zh7u4PP/ywDx48uNLLKyK1AzDHk+XhZE9UxQ14CzgivH878IdU/eMTvLv7k08+6fn5+d65c2fv1q2bz5o1yydPnuznnnuuu7sXFhZ6z549/Y033vDp06d7v379/Ec/+pEfeeSRfskll3hRUVGplXHSSSd5hw4dvHPnzt6hQwe//PLL/euvv67kKnYfM2aM77///p6fn+/5+fkevywDBgzwL774wt3dly1b5j169PB27dr54MGDfceOHe7uHolE/PLLL/e2bdt6p06dYh9CydZDIhMmTPC2bdt6hw4d/JVXXom1r1y50k888UTv3Lmzn3zyyb5q1apKL6+I1A6pEnxGz2Q1sy7AA0BdYDkw2t2/Tta/Mmeyzpgxg7vuuouXX365QtOLiOyNUp3JmtFiY+4+H1AJQxGRGrDXV5OM6tu3L3379q3pMEREao29rlTBn/70J7799tuaDgMIPlQyVRztwgsvpHnz5rFj/Mtr7ty5dO7cmfbt23PVVVfFDqN85JFHaNasWeyY/gceeKDKYp4xYwZnnnlmmf369+8fO/5fRDJHCb6E6IVEatqoUaOYNm1ahae/7LLLmDRpUqz8QfxYQ4cOZf78+cyfP7/MC3e7e5lnBJfH9u3b2bhxI4ccckiVjSkiidXaBL9t2zYGDhxIfn4+nTp1YurUqUycOJE1a9bQr18/+vXrBwSX7DvuuONiZ2lu3boVCOqejx8/np49e9KzZ08+/fTTpPMaNWoU1113Hf369WP8+PF88MEHHH/88XTt2pXjjz8+diz89u3bGTZsGHl5eQwdOpTt27cDwbH0o0aNolOnTnTu3Jk//vGPlV7+Pn360Lhx41Lty5Yto3///hxzzDGceOKJCc/AXbt2Ld988w3HHXccZsYFF1wQK4mQjpUrV3LkkUdy+eWX061bNz7//HMuu+wyunfvztFHH81tt90W6ztt2jQ6duzICSecUOzM4TfffDP2LaFr165s2bIFCLbyo7vSWrduzW233Ua3bt3o3LlzbFmSTSsi5ZTs8JqauMUfWvj000/7RRddFHu8adMmd3c/7LDD/KuvvnJ396+++spPPPFE37p1q7u7//a3v/U77rgj1m/ChAnu7j558mQfOHBg0sOMRo4c6QMHDowdN75582bfvXu3u7u/9tpr/sMf/tDd3e+++24fPXq0u7svWLDAc3Nzffbs2T5nzhw/9dRTY+MlOuzyscceix1CGX8777zzksa1YsUKP/roo4u1nXzyyf7xxx+7u/t7773n/fr1KzXd7Nmz/ZRTTok9njlzZmz5H374YT/ooIO8c+fOft555/lnn32WcL5mVuxwzA0bNrh7cGjqSSed5AsWLPDt27d7y5Yt/eOPP/ZIJOJDhgyJzefMM8/0t99+293dt2zZEluf48aNix2Tf9hhh/nEiRPd3f3+++/3MWPGpJxWREojxWGStXYLvnPnzrz++uuMHz+et956i/32269Un/fee4/FixfTu3dvunTpwuTJk1m1alXs+eHDh8f+zpo1K+X8hgwZErsY9ubNmxkyZAidOnXi2muvZdGiRQDMnDmT888/HyBWGhigbdu2LF++nHHjxjFt2jQaNWpUavwRI0bEdovE355++um018nWrVt59913GTJkSKzwWKIrN3mKkgiDBg1i5cqVFBQUcOqppzJy5MiE8zrssMPo1atX7PHf//53unXrRteuXVm0aBGLFy9m6dKltGnThsMPPxwzi60bgN69e3PdddcxceJENm3aFKu/884773DCCSfE+v3whz8E4JhjjoldSD3ZtCJSPrU2wXfo0CH2Q+FNN93EL3/5y1J93J0f/OAHsWS5ePHiYqfhx5cqKKtsQbTcAQSFufr168fChQt56aWX2LFjR8pxDjjgABYsWEDfvn25//77E+7XnjJlSsLCZenUyomKRCLsv//+xT4glixZQlFRUWy8W2+9lZYtWxarZhlfEqFJkyaxi5dffPHFzJ07t8z1sWLFCu666y7eeOMNCgoKGDhwYGydJFuvN954Iw888ADbt2+nV69eLF26lOXLl3PooYdSt27dWL9oLLm5ubHfPxJNKyLlV2s3jdasWUPjxo05//zzadiwIY888ggA++67L1u2bKFp06b06tWLK664gk8//ZT27dvz7bffsnr1ajp06AAExcJuvPFGpk6dynHHpV+KfvPmzbEfAaPzhWC/+JQpU2LJv6CgAID169dTt25dzjvvPNq1a8eoUaNKjTlixAhGjBhRsZURatSoEW3atOGpp55iyJAhuDsFBQXk5+czf/78Yn333Xdf3nvvPY499lgeffRRxo0bBwT751u0aAHAiy++mLAGT0nffPMNDRo0YL/99mPdunW8+uqr9O3bl44dO7JixQqWLVtGu3bteOKJJ2LTLFu2LHZd2VmzZrF06VK++OIL+vfvX+b8Ek0bvfC4iKSv1ib4jz76iOuvv56cnBz22Wcf/vKXvwAwduxYBgwYQIsWLZg+fTqPPPIIw4cPZ+fOnQBMmDAhluB37tzJscceSyQSKZZ8ynLDDTcwcuRI7rnnHk4++eRY+2WXXcbo0aPJy8ujS5cu9OzZEwhK9Y4ePTp2tEmiAmnlNXz4cGbMmMH69etp2bIld9xxB2PGjGHKlClcdtllTJgwgd27dzNs2DDy8/NLTf+Xv/yFUaNGsX37dgYMGBArXzxx4kRefPFF6tSpQ+PGjYt9gCWTn59P165dOfroo2nbti29e/cGoH79+kyaNImBAwfStGlTTjjhBBYuXAgERztNnz6d3NxcjjrqKAYMGMDgwYO57777ypxfomlFpPyy9qLbrVu3Zs6cObq4RS2xc+dOevfunbHzBkS+q3TRbalx9erVU3IXqWa1dhdNZUWPyIh355138tRTTxVrGzJkCDfffHM1RSUiUn2ydheNiMh3gXbRiIh8BynBi4hkKSV4EZEspQQvIpKllOBFRLKUEryISJZSghcRyVJJE7yZ5ZrZJWb2KzPrXeK5X2Q+NBERqYxUW/D/B5wEbAAmmtk9cc/9MKNRiYhIpaVK8D3d/cfu/ifgWKChmT1rZvWA1MXVRUSkxqWqRRO7KoO7FwJjzexW4N9Aw3QGN7OVwBagCChMdjqtiIhUvVQJfo6Z9Xf3adEGd/+lma0B/lKOefRz9/UVjlBERCok6S4adz8/PrnHtT/g7vtkNiwREamsMg+TNLPcSozvwL/MbK6ZjU0y/lgzm2Nmc7766qtKzEpEROKlTPBmti/wQiXG7+3u3YABwBVm1qdkB3ef5O7d3b17s2bNKjErERGJl+o4+BbA68Ckig7u7mvCv18CzwE9KzqWiIiUT6ot+LeA37r7ixUZ2MwahN8AMLMGwGnAwoqMJSIi5ZfqKJqvgUMqMfaBwHNmFp3P44l+tBURkcxIleD7An83M3f3+8s7sLsvB/IrGpiIiFROqsMktwFnAV2rLxwREakqqbbgcfci4KJqikVERKqQygWLiGSptBK8me1jZvPNrEemAxIRkaqR7hb82QTFxy7OYCwiIlKF0k3wY4ALgb5m9v0MxiMiIlUknVo0hwLN3f094HlgaKaDEhGRyktnC3408Gh4/2GCrXkREanlyio2ZsD5wN8A3H0JkGtmR1RDbCIiUgllbcHvC1zj7hvj2i5Hl+wTEan1yjrR6RvglRJtH2Y0IhERqRIpEzyAmXUArgcOi+/v7idnMC4REamkMhM88BTwv8BfCS6eLSIie4F0Enyhu5fnItsiIlILpHOY5EtmdrmZtTCzxtFbxiMTEZFKSWcLfmT49/q4NgfaVn04IiJSVcpM8O7epjoCERGRqpU0wZvZye7+bzP7YaLn3f3ZzIUlIiKVlWoL/iTg38CgBM85oAQvIlKLJU3w7n5b+Hd09YUjIiJVJZ0Tna5L0LwZmOvu86s8IhERqRLpHCbZHbgUOCS8jQX6An81sxsyF5qIiFRGOgm+CdDN3X/q7j8lSPjNgD7AqLImNrNcM/vQzF6uVKQiIlIu6ST4VsCuuMe7gcPcfTuwM43prwaWVCA2ERGphHROdHoceM/MXggfDwKeMLMGwOJUE5pZS2AgcCeQaF++iIhkSDonOv3KzF4BTiCoA3+pu88Jnx5RxuR/Am4gqCsvIiLVKNWJTo3c/Zuw7syK8BZ9rnGJi4Akmv5M4Et3n2tmfVP0G0vwwy2tWrUqX/QiIpJUqi34x4EzgbkEJzZFGenVoukNnGVmZwD1gUZm9pi7nx/fyd0nAZMAunfv7qWHERGRikh1otOZ4TVZT3L3z8o7sLvfBNwEEG7B/6xkchcRkcxJeRSNuzvwXDXFIiIiVSidwyTfM7MelZmJu89w9zMrM4aIiJRPOodJ9gMuMbNVwDbCffDunpfRyEREpFLSSfADMh6FiIhUuTJ30bj7KmB/ghOcBgH7h20iIlKLlZngzexqYArQPLw9ZmbjMh2YiIhUTjq7aMYAx7r7NgAz+x0wC7gvk4GJiEjlpHMUjQFFcY+LwjYREanF0tmCfxh438yix8OfAzyYsYhERKRKpFNs7B4zm8GeYmOj3f3DTAcmIiKVk84WPO4+D5iX4VhERKQKpbMPXkRE9kJK8CIiWSqtXTQQ1IeP719WPXgREalZZSZ4M7sE+CWwnT114dOpBy8iIjUonS34nwFHu/v6TAcjIiJVJ5198MuAbzMdiIiIVK10tuBvAt41s/eBndFGd78qY1GJiEilpZPg/w/4N/AREMlsOCIiUlXSSfCF7n5dxiMREZEqlc4++OlmNtbMWphZ4+gt45GJiEilpLMF/+Pw701xbTpMUkSklkun2Fib6ghERESqVjonOu0DXAb0CZtmAP/n7rszGJeIiFRSOrto/gLsA/xP+PgnYdtFmQpKREQqL50E38Pd8+Me/9vMFpQ1kZnVB2YC9cL5PO3ut1UsTBERKa90jqIpMrN20Qdm1pbil/BLZidwcvjh0AXob2a9KhSliIiUW7q1aKab2XKCKzodBowuayJ3d2Br+HCf8ObJpxARkaqUMsGbWS6QDxwOHEGQ4Je6+85U05WYfi7QHrjf3d9P0GcsMBagVatW5QpeRESSS7mLxt2LgLPcfae7F7j7gnSTe3R6d+8CtAR6mlmnBH0muXt3d+/erFmz8sYvIiJJpLOL5l0z+zMwFdgWbQyv05oWd98UXri7P7CwvEGKiEj5pZPgjw///jKuzYGTU01kZs2A3WFy/x5wKvC7CkUpIiLlljTBm9nV7n4vcIu7v12BsVsAk8P98DnA39395QrGKSIi5ZRqC340cC8wEehW3oHdvQDoWsG4RESkklIl+CVmthJoZmYFce1GcBRkXkYjExGRSkma4N19uJkdBPwTOKv6QhIRkaqQ8kdWd/8vwXHwIiKyl0mnVIGIiOyFlOBFRLKUEryISJZK54IfzYDxwFFA/Wi7u6c80UlERGpWOlvwU4AlQBvgDmAlMDuDMYmISBVIJ8E3cfcHCcoOvOnuFwKq6y4iUsulU4smeu3VtWY2EFhDUB1SRERqsXQS/AQz2w/4KXAf0Ai4JpNBiYhI5aWT4L92983AZqAfgJn1zmhUIiJSaensg78vzTYREalFUpULPo6gFnwzM7su7qlGQG6mAxMRkcpJtYumLtAw7LNvXPs3wOBMBiUiIpWXqprkm8CbZvaIu68yswbuvi1ZfxERqV3S2Qd/sJktJjjZCTPLN7P/yWxYIiJSWekk+D8BpwMbANx9AdAngzGJiEgVSKvYmLt/XqKpKAOxiIhIFUrnOPjPzex4wM2sLnAV4e4aERGpvdLZgr8UuAI4BFgNdAkfi4hILVbmFry7rwdGVEMsIiJShdKtB38x0Dq+f1hVMtV0hwKPAgcBEWCSu99bmWBFRCR96eyDfwF4C3id8v24Wgj81N3nmdm+wFwze83dF1cgThERKad0Evz33X18eQd297XA2vD+FjNbQrAfXwleRKQapPMj68tmdkZlZmJmrYGuwPuVGUdERNKXqtjYFsABA35uZjsJLv5hgLt7o3RmYGYNgWeAa9z9mwTPjwXGArRq1arcCyAiIomlqkWzb7Ln0mVm+xAk9ynu/myS+UwCJgF0797dKztPEREJlLmLxszeSKctQR8DHgSWuPs9FQtPREQqKtUumvpAA6CpmR1AsGsGgnrwB6cxdm/gJ8BHZjY/bPu5u79S8XBFRCRdqY6iuYTg2qsHA3PZk+C/Ae4va2B3fztuGhERqWap9sHfC9xrZuPcXZfoExHZy6TcB29m3wNmlmhrZWaHZDQqERGptLJ+ZC0EnjWzBnFtDwAtMheSiIhUhZQJ3t13A88BQyHYegeaufucaohNREQqIZ0zWR8ARof3LwAezlw4IiJSVdIpF7zUzDCzDsBw4ITMhyUiIpWV1iX7CE5YegAocPevMxiPiIhUkXQT/N+BfIJELyIie4F0ygXj7t8C+2U4FhERqULpbsGLiMheRgleRCRLKcGLiGQpJXgRkSylBC8ikqWU4EVEspQSvIhIllKCFxHJUkrwIiJZSgleRCRLKcGLiGQpJXgRkSylBC8ikqWU4EVEslTGEryZPWRmX5rZwkzNQ0REksvkFvwjQP8Mji8iIilkLMG7+0xgY6bGFxGR1Gp8H7yZjTWzOWY256uvvqrpcEREskaNJ3h3n+Tu3d29e7NmzWo6HBGRrFHjCV5ERDJDCV5EJEtl8jDJJ4BZwBFmttrMxmRqXiIiUlqdTA3s7sMzNbaIiJRNu2hERLKUEryISJZSghcRyVJK8CIiWUoJXkQkSynBi4hkKSV4EZEspQQvIpKllOBFRLKUEryISJZSghcRyVJK8CIiWUoJXkQkSynBi4hkKSV4EZEspQQvIpKllOBFRLKUEryISJZSghcRyVJK8CIiWUoJXkQkSynBi4hkqYwmeDPrb2b/MbNPzezGTM5LRESKy1iCN7Nc4H5gAHAUMNzMjsrU/EREpLhMbsH3BD519+Xuvgt4Ejg7g/MTEZE4mUzwhwCfxz1eHbaJiEg1qJPBsS1Bm5fqZDYWGBs+3GlmCzMYU2U0BdbXdBApKL7KUXyVo/gqpzLxHZbsiUwm+NXAoXGPWwJrSnZy90nAJAAzm+Pu3TMYU4XV5thA8VWW4qscxVc5mYovk7toZgOHm1kbM6sLDANezOD8REQkTsa24N290MyuBP4J5AIPufuiTM1PRESKy+QuGtz9FeCVckwyKVOxVIHaHBsovspSfJWj+ConI/GZe6nfPUVEJAuoVIGISJaqFQm+NpQ0MLNDzWy6mS0xs0VmdnXYfruZfWFm88PbGXHT3BTG/B8zO70aYlxpZh+FccwJ2xqb2Wtm9kn494CaiM/MjohbR/PN7Bszu6Ym15+ZPWRmX8YfeluR9WVmx4Tr/VMzm2hmiQ4BrorY/mBmS82swMyeM7P9w/bWZrY9bh3+byZjSxFfuV/Lao5valxsK81sftheE+svWT6p3vefu9fojeAH2GVAW6AusAA4qgbiaAF0C+/vC3xMUGLhduBnCfofFcZaD2gTLkNuhmNcCTQt0fZ74Mbw/o3A72oqvhKv6X8Jjs+tsfUH9AG6AQsrs76AD4DjCM7teBUYkKHYTgPqhPd/Fxdb6/h+Jcap8thSxFfu17I64yvx/N3ArTW4/pLlk2p9/9WGLfhaUdLA3de6+7zw/hZgCanPvD0beNLdd7r7CuBTgmWpbmcDk8P7k4Fz4tprKr5TgGXuvipFn4zH5+4zgY0J5pv2+jKzFkAjd5/lwX/bo3HTVGls7v4vdy8MH75HcO5IUpmKLVl8KVTruisrvnAL90fAE6nGyHB8yfJJtb7/akOCr3UlDcysNdAVeD9sujL82vxQ3FeqmojbgX+Z2VwLzgAGONDd10LwpgKa12B8UcMo/s9VW9YflH99HRLer+44LyTYWotqY2YfmtmbZnZi2FYTsZXntaypdXcisM7dP4lrq7H1VyKfVOv7rzYk+LRKGlQXM2sIPANc4+7fAH8B2gFdgLUEX/2gZuLu7e7dCCp0XmFmfVL0rZH1asFJbWcBT4VNtWn9pZIsnmqP08xuBgqBKWHTWqCVu3cFrgMeN7NGNRBbeV/LmnqNh1N8A6PG1l+CfJK0a5JYKhVjbUjwaZU0qA5mtg/BizHF3Z8FcPd17l7k7hHgr+zZjVDtcbv7mvDvl8BzYSzrwq9x0a+cX9ZUfKEBwDx3XxfGWmvWX6i862s1xXeVZDROMxsJnAmMCL+SE35t3xDen0uwf7ZDdcdWgdeyWuMDMLM6wA+BqXFx18j6S5RPqOb3X21I8LWipEG43+5BYIm73xPX3iKu27lA9Ff7F4FhZlbPzNoAhxP8GJKp+BqY2b7R+wQ/yC0M4xgZdhsJvFAT8cUptvVUW9ZfnHKtr/Br9BYz6xW+Ry6Im6ZKmVl/YDxwlrt/G9fezILrK2BmbcPYlldnbOG8y/VaVnd8oVOBpe4e261RE+svWT6hut9/VfGLcWVvwBkEvzIvA26uoRhOIPjqUwDMD29nAH8DPgrbXwRaxE1zcxjzf6iiX99TxNeW4Ff2BcCi6HoCmgBvAJ+EfxvXRHzh/L4PbAD2i2ursfVH8EGzFthNsCU0piLrC+hOkMyWAX8mPEEwA7F9SrAfNvr++9+w73nha74AmAcMymRsKeIr92tZnfGF7Y8Al5boWxPrL1k+qdb3n85kFRHJUrVhF42IiGSAEryISJZSghcRyVJK8CIiWUoJXkQkSynBi4hkKSV4yRpmtr+ZXR73+GAzezpD8zrHzG6t4LSvx5eJFckUHQcvWSMs6vSyu3eqhnm9S3DG6foKTDsSaOnud1Z9ZCJ7aAtesslvgXYWXNThDxZc6GEhgJmNMrPnzewlM1thZlea2XVhhcH3zKxx2K+dmU0LK3a+ZWYdS87EzDoAO6PJ3cweCS/E8K6ZLTezwWF7CzObGcazMK6K4YsEJR1EMkoJXrLJjQR16Lu4+/UJnu8E/JigSNadwLceVBicRVDjA4KLH49z92OAnwH/k2Cc3gSnvMdrQXB6+pkEHzSE8/qnu3cB8glOV8fdvwbqmVmTCiyjSNrq1HQAItVougcXX9hiZpuBl8L2j4C8sLTr8cBTtueqaPUSjNMC+KpE2/MeVFlcbGYHhm2zgYfCqoLPu/v8uP5fAgcT1O4RyQhtwct3yc64+5G4xxGCjZ0cYFP4DSB6OzLBONuB+inGNohddagP8AXwNzO7IK5P/XAckYxRgpdssoXg+pcV4sEFGVaY2RAISr6aWX6CrkuA9mWNZ2aHAV+6+18JSsd2i44LHERwjV2RjFGCl6zhwUUd3gl/0PxDBYcZAYwxs2hZ5kTXB54JdDUr8+r2fYH5ZvYhQcnae8P2Y4D3fM/1V0UyQodJilSAmd0LvOTur1dw2hfd/Y2qj0xkD23Bi1TMrwkucFIRC5XcpTpoC15EJEtpC15EJEspwYuIZCkleBGRLKUELyKSpZTgRUSy1P8DIfMjZyEkWdcAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "selection = df[['time', 'rotation_sq']].groupby(\"time\").mean()\n",
    "\n",
    "lm_time = np.arange(0, 2000)\n",
    "\n",
    "line_y = np.linspace(0, np.max(selection[\"rotation_sq\"]), 2000)\n",
    "\n",
    "lm = np.poly1d(np.polyfit(df[\"time\"], df[\"rotation_sq\"], 1))\n",
    "   \n",
    "d = calc_rot_d_from_step(step=lm.coef[0])\n",
    "d_text = f\"meas D = {d:.2e}\"\n",
    "exp_text = f\"exp D = {convert_d_from_rads2_s_to_rads2_ns(d=2e3):.2e}\" # 2e3 rads2/s\n",
    "    \n",
    "step_text = f\"step_rads = {round(lm.coef[0], 5)} rads/ns\"\n",
    "\n",
    "print(d_text)\n",
    "print(exp_text)\n",
    "print(step_text)\n",
    "\n",
    "# plot the displacement_sq against the fit\n",
    "\n",
    "plt.suptitle(\"angle rotation of LHCII from origin\")\n",
    "plt.title(\"D in rad2/ns units\")\n",
    "plt.xlabel(\"time (ns)\")\n",
    "plt.ylabel(\"<theta from origin>^2\")\n",
    "plt.xlim(0, np.max(df[\"time\"]) * 1.05)\n",
    "plt.ylim(0, np.max(selection[\"rotation_sq\"]) * 1.05)\n",
    "plt.annotate(text=d_text, xy=(20, np.max(selection[\"rotation_sq\"]) * 0.95))\n",
    "plt.annotate(text=step_text, xy=(20, np.max(selection[\"rotation_sq\"]) * 0.85))\n",
    "plt.annotate(text=exp_text, xy=(20, np.max(selection[\"rotation_sq\"]) * 0.90))\n",
    "plt.plot(selection[\"rotation_sq\"])\n",
    "plt.plot(lm(lm_time))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d0 = 0.00018 # d0 in nm2/ns\n",
    "\n",
    "d1 = convert_d_from_nm2_ns_to_cm2_s(d0) # convert to cm2/s units\n",
    "\n",
    "d2 = convert_d_from_cm2_s_to_nm2_ns(d1) # convert the cm2/s units back to nm2/ns units\n",
    "\n",
    "print(f\"d0==d2 is {d0 == d2}\")\n",
    "\n",
    "step_nm = calc_step(d0)\n",
    "\n",
    "step_cm = calc_step(d1)\n",
    "\n",
    "print(f\"d0: {d0}, d1: {d1}, d2: {d2}\")\n",
    "\n",
    "print(f\"nm/ns step: {step_nm}\")\n",
    "\n",
    "print(f\"cm/s step: {step_cm}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_d = pd.read_csv(\"C:\\\\projects\\\\grana_model\\\\src\\\\grana_model\\\\res\\log\\\\d_calc.csv\")\n",
    "print(df_d.mean())\n",
    "\n",
    "df_d.hist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from src.grana_model.dcalibrator import DCalibrator\n",
    "\n",
    "structure_dict={\n",
    "                \"LHCII\": {\n",
    "                    \"d\": 1.8e-9,  # 1.8e-9 in cm2/s\n",
    "                    \"d_rot\": 2e13,  # 2 x 10^3  rad^2 s^(-1)\n",
    "                    \"simulation_limit\": 100000,\n",
    "                    \"distance_scalar\": \"well\",\n",
    "                    \"diffusion_scalar\": 1220.0,  # average over 250 steps, gave us this number for keeping step_nm equal to calculated step\n",
    "                    \"distance_threshold\": 50.0,\n",
    "                    \"mass\": 1000.0,\n",
    "                    \"mass_scalar\": 1.0,\n",
    "                    \"rotation_scalar\": 180.0,  # average over 250 steps, gives us <d_rot> of 2.00e13\n",
    "                    \"time_per_step\": 2,  # in ns\n",
    "                    \"average_step_over\": 250,\n",
    "                    \"calibrate_rot_d\": False,\n",
    "                    \"calibrate_diff_d\": False,\n",
    "                }\n",
    "}\n",
    "\n",
    "dcal = DCalibrator(structure_dict[\"LHCII\"])\n",
    "\n",
    "d_ns = dcal.convert_d_from_rads2_s_to_rads2_ns(d=2e13)\n",
    "print(d_ns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(\"lhcii_circle_4p5_6_limit_100_coords.csv\")\n",
    "\n",
    "print(df.head())"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.12 ('sim')",
   "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.12"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "f9bc5cb7ec0992daba0695a9823e244b71b39dc358dc9742cfefbb2ae8088a2b"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}