{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "c4aa3266-0663-4fa1-929f-58eabafebdfb",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import curve_fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "81ce782c-4fb1-4a31-906e-3b1b5f04d9dd",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArcAAAIhCAYAAABUopIpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA6JUlEQVR4nO3de3RU1d3/8c8QkgmZJFxNQiBA1Li4C4jSBipQ5SJgpdRLQRHEKha8BLwUpH0MWsNFDVQiWKwPIJSqFLloRUEQkAIlKKDirf6MFIEQxJSQgdzI/v3Bk5FJJpAEyBk279das5azzz7nfM/J1vVxz54zLmOMEQAAAGCBOk4XAAAAAJwrhFsAAABYg3ALAAAAaxBuAQAAYA3CLQAAAKxBuAUAAIA1CLcAAACwBuEWAAAA1iDcAgAAwBqEWwDVNn/+fLlcrkpf69evd7rE01q/fv0FUefbb7+t1NTUszrGyJEj1apVqxrv36tXL/Xq1atG+6alpWn58uU1PvfpFBcXq3Xr1po6dapf+2uvvaZ27dqpXr16crlc2rlzp1JTU+Vyuc5LHacTaJwNHz5cgwcPrvVagItJXacLAHDhmjdvnlq3bl2hvW3btg5UU3VdunTRli1bgr7Ot99+Wy+88MJZB1ynpKWl6eabbz4vYW727NnKzc3VAw884Gs7dOiQhg8frv79+2v27Nlyu9264oor9Jvf/Eb9+/c/5zXURGpqqlq3bq1169bp5z//udPlAFYi3AKosfbt26tr165Ol1FlxcXFcrlcio6O1k9+8hOny0ENlZSU6JlnntGoUaPk8Xh87V999ZWKi4t1xx13qGfPnr72iIgINW/e3IlSK7jsssvUv39/TZ06lXALnCcsSwBw3rz66qtyuVzKyMjwa3/iiScUEhKiNWvWSJK+/fZbuVwuTZ8+XU8//bRatGih8PBwde3aVWvXrq1w3H//+98aNmyYYmJi5Ha71aZNG73wwgt+fco+El64cKEefvhhNWvWTG63W19//XXAj4tHjhypyMhIffHFF+rXr588Ho+aNm3q+9h769at6tGjhzwej6644gotWLCgQl3Z2dkaPXq0mjdvrrCwMCUmJmry5MkqKSnx9Sm71meffVbp6elKTExUZGSkfvrTn2rr1q1+9ZRd06lLPr799ltJ0gsvvKBrr71WMTEx8ng86tChg6ZPn67i4uJq/IV+ZIzR9OnT1bJlS4WHh6tLly5atWpVhX4FBQV6+OGH1alTJ9WvX1+NGjXST3/6U61YscKvn8vlktfr1YIFC3y1ly1vOHTokMaMGaO2bdsqMjJSMTEx+vnPf64PPvigSrWuXLlS+/bt0/Dhw31tI0eOVI8ePSRJt912m9/5yi9L2LRpk0JDQ/XII4/4Hbdsuc3LL7/sa6vKWJOkL774Qv3791dERISaNGmi++67T0ePHg1Y//Dhw/Xee+/p//2//1el6wVQTQYAqmnevHlGktm6daspLi72e5WUlPj1ve+++0xYWJjJzMw0xhizdu1aU6dOHfP73//e1ycrK8tIMgkJCaZHjx5m6dKlZsmSJebqq682oaGhZvPmzb6+u3fvNvXr1zcdOnQwr7zyilm9erV5+OGHTZ06dUxqaqqv3/vvv28kmWbNmpmbb77ZrFy50rz11lvm8OHDvm3vv/++r/+IESNMWFiYadOmjfnTn/5k1qxZY+666y4jyUycONFcccUV5uWXXzbvvvuuGTRokJFktm/f7tv/wIEDJiEhwbRs2dL8+c9/Nu+995556qmnjNvtNiNHjqxwra1atTL9+/c3y5cvN8uXLzcdOnQwDRs2NP/973+NMcZ8/fXX5uabbzaSzJYtW3yvgoICY4wx48aNM3PmzDHvvPOOWbdunZkxY4Zp0qSJueuuu/zu/4gRI0zLli3P+Dd94oknjCRz9913m1WrVpm5c+eaZs2ambi4ONOzZ09fv//+979m5MiRZuHChWbdunXmnXfeMY888oipU6eOWbBgga/fli1bTL169cyAAQN8te/evdsYY8wXX3xhfvvb35pXX33VrF+/3rz11lvm7rvvNnXq1PH7m1Rm1KhRJiYmxq/t66+/Ni+88IKRZNLS0vzOV3Ztp5o6daqRZFasWGGMMebTTz81ERER5o477vD1qepYy87ONjExMaZZs2Zm3rx55u233za33367adGiRYVxZowxBw8eNJLM888/f8ZrBVB9hFsA1VYWbgO9QkJC/PoWFBSYzp07m8TERPPZZ5+Z2NhY07NnT78QXBb44uPjzfHjx33teXl5plGjRub666/3tfXr1880b97cHDlyxO88999/vwkPDzc//PCDMebHcHvttddWqL+ycCvJLF261NdWXFxsLrnkEiPJfPTRR772w4cPm5CQEDN+/Hhf2+jRo01kZKTZs2eP37meffZZI8kXtMqutUOHDn73YNu2bUaS+dvf/uZrGzt2bIVQFsiJEydMcXGxeeWVV0xISIjvHpRd15nCbW5urgkPDze//OUv/dr/+c9/Gkl+4ba8kpISU1xcbO6++27TuXNnv20ej8eMGDHijPWXHeO6666rUEMgbdq0Mf3796/QXvZ3XbJkiV97oHBbWlpqBgwYYBo0aGA+/fRT07ZtW9O6dWuTn5/v61PVsfa73/3OuFwus3PnTr9+ffr0CRhujTGmWbNm5rbbbjvjtQKoPpYlAKixV155RZmZmX6vf/3rX3593G63Xn/9dR0+fFhdunSRMUZ/+9vfFBISUuF4Q4YMUXh4uO99VFSUbrzxRm3cuFEnTpxQQUGB1q5dq1/+8peKiIhQSUmJ7zVgwAAVFBT4fbQvSb/61a+qfD0ul0sDBgzwva9bt64uv/xyNW3aVJ07d/a1N2rUSDExMdqzZ4+v7a233lLv3r0VHx/vV9cNN9wgSdqwYYPfuQYOHOh3Dzp27ChJfsc8nR07dugXv/iFGjdurJCQEIWGhurOO+/UiRMn9NVXX1X5miVpy5YtKigo0O233+7XnpycrJYtW1bov2TJEnXv3l2RkZGqW7euQkND9fLLL+vzzz+v8jlffPFFdenSReHh4b5jrF27tkrH2L9/v2JiYqp8rkBcLpdeeeUVRUVFqWvXrsrKytLrr7/uW8NbnbH2/vvvq127drryyiv9zjFs2LBKzx8TE6N9+/ad1TUACIxwC6DG2rRpo65du/q9rrrqqgr9Lr/8cv3sZz/zBaimTZsGPF5cXFzAtqKiIuXn5+vw4cMqKSnRrFmzFBoa6vcqC6Xff/+93/6VnSuQiIgIv3AtSWFhYWrUqFGFvmFhYSooKPC9P3jwoN58880KdbVr1y5gXY0bN/Z773a7JUnHjx8/Y53/+c9/9LOf/Uz79u3Tn/70J33wwQfKzMz0rQWtyjFOdfjwYUmV3/9TvfHGG7r11lvVrFkzLVq0SFu2bFFmZqZGjRrldz9OJz09Xb/97W/VrVs3LV26VFu3blVmZqb69+9fpdqPHz9e4e9UE40bN9YvfvELFRQUqH///urQoYNvW3XG2uHDh6t0704VHh5e7b8TgKrhaQkAzru//OUv+sc//qFrrrlGGRkZuu2229StW7cK/bKzswO2hYWFKTIyUqGhoQoJCdHw4cM1duzYgOdKTEz0e19bzzdt0qSJOnbsqKeffjrg9vj4+HN2ruXLl8vr9eqNN97wm1nduXNnjY5XFrQru/+nPid30aJFSkxM1GuvveZ3bwsLC6t8vkWLFqlXr16aM2eOX3tlX8Aqr0mTJvrhhx+qfL7KrFmzRnPmzNE111yjZcuWaenSpb6Z/oYNG1Z5rDVu3LjSe1eZH3744ayePwygcoRbAOfVJ598ogcffFB33nmnXnrpJSUnJ+u2227Tjh071LBhQ7++b7zxhp555hnfrNzRo0f15ptv6mc/+5lCQkIUERGh3r17a8eOHerYsaPCwsKcuKSABg0apLfffluXXXZZheuqqVNnc+vVq+drLwuVZdulk087eOmll2p0np/85CcKDw/XX//6V79lHJs3b9aePXv8QpjL5VJYWJhfsM3Ozq7wtISy+gLNTrpcLr/aJenjjz/Wli1blJCQcMZ6W7dufdZPGjhw4IDvkWFr1qzRkCFDdPfdd6tLly5KTEys1ljr3bu3pk+frl27dvktTVi8eHHA/iUlJdq7d6/fEhgA5w7LEgDU2KeffqqtW7dWeB06dEiS5PV6deuttyoxMVGzZ89WWFiYXn/9df33v//VXXfdVeF4ISEh6tOnj28W7brrrlNeXp4mT57s6/OnP/3J97H8/PnztX79er355puaMWOGo88NffLJJxUaGqrk5GTNmTNH69at09tvv63Zs2dr0KBB+u6776p9zLKPyadNm6Z//etf2r59u4qKitSnTx+FhYVp6NChWrVqlZYtW6Z+/fopNze3RrU3bNhQjzzyiJYtW6bf/OY3evfdd/WXv/xFt956a4WP1gcNGqQvv/xSY8aM0bp167RgwQL16NEj4PKPDh06+P4+27dv15dffuk7xurVq/XEE09o3bp1mjNnjvr161dh1r0yvXr10ieffKJjx47V6HpPnDihoUOHyuVyafHixQoJCdH8+fNVv3593XbbbSoqKpJU9bGWkpKiJk2aaODAgZo/f75WrVqlO+64Q1988UXA83/88cc6duyYevfuXaP6AZyB099oA3DhOd3TEiSZl156yRhjzB133GEiIiJ8Twoos2TJEiPJzJgxwxjz4xMEpk2bZiZPnmyaN29uwsLCTOfOnc27775b4fxZWVlm1KhRplmzZiY0NNRccsklJjk52fzxj3/09ansm/Onbiv/tASPx1Ohb8+ePU27du0qtLds2dIMHDjQr+3QoUPmwQcfNImJiSY0NNQ0atTIXHXVVWbSpEm+b+GXXeszzzxT4ZiSzBNPPOF7X1hYaH7zm9+YSy65xLhcLiPJZGVlGWOMefPNN82VV15pwsPDTbNmzcyjjz5qVq1aFfC6qvIosNLSUjNlyhSTkJBgwsLCTMeOHc2bb75pevbsWeFpCVOnTjWtWrUybrfbtGnTxrz00ksBn0iwc+dO0717dxMREeH31IXCwkLzyCOPmGbNmpnw8HDTpUsXs3z58irX+vXXXxuXy2Vef/11v/aqPi1h0qRJpk6dOmbt2rV+/TZv3mzq1q1rHnroIV9bVcaaMcZ89tlnpk+fPiY8PNw0atTI3H333WbFihUBn5bwhz/8wTRp0sT3WDcA55bLGGNqN04DgL9vv/1WiYmJeuaZZyo8WB8I5MYbb1RJSUnAH5oIZidOnNDll1+uYcOGVbo+G8DZYVkCAOCCM2XKFL333nvKzMx0upRqWbRokfLz8/Xoo486XQpgLcItAOCC0759e82bN++0TyQIRqWlpfrrX/+qBg0aOF0KYC2WJQAAAMAazNwCAADAGoRbAAAAWINwCwAAAGvwC2U6ucB///79ioqKqrWf6gQAAEDVGWN09OhRxcfHq06dyudnCbeS9u/fX6WffAQAAICz9u7dq+bNm1e6nXArKSoqStLJmxUdHe1wNcGtsLBQ6enpGj9+fIXfhsfFiTGBQBgXKI8xgfKqOyby8vKUkJDgy22VIdxKvqUI0dHRhNszKCwsVHh4uKKjo/mPEyQxJhAY4wLlMSZQXk3HxJmWkPKFMgAAAFiDcAsAAABrEG4BAABgDdbcVpExRiUlJTpx4oTTpTiqqKhIHo9HhYWFqu4vN4eGhiokJOQ8VQYAAEC4rZKioiIdOHBAx44dc7oUxxlj1L17d3333XfVfiawy+VS8+bNFRkZeZ6qAwAAFzvC7RmUlpYqKytLISEhio+PV1hY2EX9Qw+lpaX6/vvv1aRJk9M+QLk8Y4wOHTqk7777TklJSczgAgCA84JwewZFRUUqLS1VQkKCIiIinC7HcaWlpapbt67Cw8OrFW4l6ZJLLtG3336r4uJiwi0AADgv+EJZFVU3yKGii3nGGwAA1A4SWy3xeiWX6+TL63W6GgAAADsRbgEAAGANwi0AAACsQbi12MiRI+VyueRyuRQaGqrY2Fj16dNH//u//6vS0tIqH2f+/Plq0KDB+SsUAADgHCHcWq5///46cOCAvv32W61atUq9e/fWQw89pEGDBqmkpMTp8gAAAM4pwu154PUGflV1+7nkdrsVFxenZs2aqUuXLnr88ce1YsUKrVq1SvPnz5ckpaenq0OHDvJ4PEpISNCYMWOUn58vSVq/fr3uuusuHTlyRC6XSyEhIXruueckSYsWLVLXrl0VFRWluLg4DRs2TDk5OefnQgAAAKqAcHseREZWfMXG/rg9NjZwn9ry85//XFdeeaXeeOMNSScfc/b888/r008/1YIFC7Ru3To99thjkqTk5GTNnDlT0dHROnDggPbt26f77rtP0slnAD/11FPatWuXli9frqysLI0cObL2LgQAAKAcfsThItW6dWt9/PHHkqSUlBRfe2Jiop566in99re/1ezZsxUWFqb69evL5XIpLi5OpaWlys7OliSNGjXKt9+ll16q559/Xtdcc43y8/P5iV0AAOAIwu158H+f6Pvxen+cvT14UPJ4arem8owxvh9VeP/995WWlqbPPvtMeXl5KikpUUFBgbxerzynKXTHjh1KTU3Vzp079cMPP/i+pPaf//xHbdu2rZXrAAAAOBXLEs4Djyfwq6rba8Pnn3+uxMRE7dmzRwMGDFD79u21dOlSffjhh3rhhRckScXFxZXu7/V61bdvX0VGRmrRokXKzMzUsmXLJJ1crgAAAOAEZm4vQuvWrdMnn3yicePGafv27SopKdFzzz3n+4nh119/3a9/WFiYTpw44df2xRdf6Pvvv9fUqVOVkJAgSdq+fXvtXAAAADivvN4fvw+Un+/8J87Vwcyt5QoLC5Wdna19+/bpo48+Ulpamm666SYNGjRId955py677DKVlJRo1qxZ+uabb7Rw4UK9+OKLfsdo1aqV8vPztXbtWn3//fc6fvy4WrRoobCwMN9+K1eu1FNPPeXQVQIAAJxEuLXcO++8o6ZNm6pVq1bq37+/3n//fT3//PNasWKFQkJC1KlTJ6Wnp2vatGlq3769/vrXv2rKlCl+x0hOTtZ9992n2267TbGxsZo9e7YuueQSzZ8/X0uWLFHbtm01depUPfvssw5dJQAAwEksS7DY/Pnzfc+yPZ1x48Zp3Lhxfm3Dhw/3ez9nzhzNmTPH72kJQ4cO1dChQ/36GWPOrmgAAICzQLitJR6PRO4DAAA4v1iWAAAAAGsQbgEAAGANR8Ptxo0bdeONNyo+Pl4ul0vLly/3226MUWpqquLj41WvXj316tVLu3fv9utTWFioBx54QE2aNJHH49EvfvELfffdd7V4FQAAAAgWjoZbr9erK6+8UhkZGQG3T58+Xenp6crIyFBmZqbi4uLUp08fHT161NcnJSVFy5Yt06uvvqpNmzYpPz9fgwYNqvBc1rPFF6XOHvcQAACcb45+oeyGG27QDTfcEHCbMUYzZ87UpEmTNGTIEEnSggULFBsbq8WLF2v06NE6cuSIXn75ZS1cuFDXX3+9JGnRokVKSEjQe++9p379+gU8dmFhoQoLC33v8/LyArZLUmlpqYwxys/Pl9vtPutrvtCVBVRjjO/ndquqsLBQxhidOHGiwn3Ghavsb8nfFKdiXKA8xsSF5eSfyf1//1youuchMVZ3TFS1n8sEyXSay+XSsmXLNHjwYEnSN998o8suu0wfffSROnfu7Ot30003qUGDBlqwYIHWrVun6667Tj/88IMaNmzo63PllVdq8ODBmjx5csBzpaamBtw2YcIEhYeHV2i//PLLlZSUpEaNGik0NFQul+ssr/biY4xRbm6uvv32W+3atcvpcgAAwGkUFYUqLe1xSdLjj6cpLKzY4YqkgoICTZ06VUeOHFF0dHSl/YL2UWBlz1KNjY31a4+NjdWePXt8fcLCwvyCbVmfsv0DmThxosaPH+97n5eXp4SEBI0fPz7gzTLG6NChQzp69KiKiopqfE02MMYoLy9P0dHR1Q75brdb119/faWz9bgwFRYWasaMGRo3bhyfbsCHcYHyGBPBy+sN3JaWdvKf77vv4Up/fvdsfpa3umMiLy9PU6dOPWO/oA23ZcoHKGPMGUPVmfq43e6AN7GydklKSEjQiRMnVFzs/P+5OKmoqEhz587Vvffeq7CwsGrtGxYWpjp1eECHrU737w8uXowLlMeYCD4BPrT206JF5X+vc/H5f1XHRFXHTdCG27i4OEknZ2ebNm3qa8/JyfHN5sbFxamoqEi5ubl+s7c5OTlKTk4+5zWFhIQoJCTknB/3QuJyueT1evmPEwAACEpBO42WmJiouLg4rVmzxtdWVFSkDRs2+ILrVVddpdDQUL8+Bw4c0Keffnpewi0AAIBt8vMrvg4e/HH7wYOB++TnO1fz6Tg6c5ufn6+vv/7a9z4rK0s7d+5Uo0aN1KJFC6WkpCgtLU1JSUlKSkpSWlqaIiIiNGzYMElS/fr1dffdd+vhhx9W48aN1ahRIz3yyCPq0KGD7+kJAAAAqNyZ1s16PGe3tra2ORput2/frt69e/vel33Ja8SIEZo/f74ee+wxHT9+XGPGjFFubq66deum1atXKyoqyrfPjBkzVLduXd166606fvy4rrvuOs2fP/+iXz4AAABwMXI03Pbq1eu0D/Z3uVxKTU1VampqpX3Cw8M1a9YszZo16zxUCAAAgAtJ0K65BQAAAKqLcAsAAABrEG4BAABgDcItAAAArEG4BQAAgDUItwAAALAG4RYAAADWINwCAADAGoRbAAAAWMPRXygDAABA8PF4pNP8iGxQY+YWAAAA1iDcAgAAwBqEWwAAAFiDcAsAAABrEG4BAABgDcItAAAArEG4BQAAgDUItwAAALAG4RYAAADWINwCAADAGoRbAAAAWINwCwAAAGsQbgEAAGANwi0AAACsQbgFAACANQi3AAAAsAbhFgAAANYg3AIAAMAahFsAAABYg3ALAAAAaxBuAQAAYA3CLQAAAKxBuAUAAIA1CLcAAACwBuEWAAAA1iDcAgAAwBqEWwAAAFiDcAsAAABrEG4BAABgDcItAAAArEG4BQAAgDUItwAAALAG4RYAAADWINwCAADAGoRbAAAAWINwCwAAAGsQbgEAAGANwi0AAACsQbgFAACANQi3AAAAsAbhFgAAANYg3AIAAMAahFsAAABYg3ALAAAAaxBuAQAAYA3CLQAAAKxBuAUAAFbyeiWX6+TL63W6GtQWwi0AAACsQbgFAACANeo6XQAAAEB1VHWJwan9qrqPx1P9ehBcCLcAAOCCEhlZ/X1iY6vWz5jqHxvBhWUJAAAAsAYztwAA4IKSn1+1fl7vjzO2Bw+y5OBiQbgFAAAXlJqEVI+HcHuxYFkCAAAArEG4BQAAgDUItwAAALAGa24BAICVPB4e7XUxYuYWAAAA1iDcAgAAwBqEWwAAAFiDcAsAAABrEG4BAABgDcItAAAArEG4BQAAgDUItwAAALAG4RYAAADWINwCAADAGoRbAAAAWINwCwAAAGsQbgEAAGANwi0AAACsEdThtqSkRL///e+VmJioevXq6dJLL9WTTz6p0tJSXx9jjFJTUxUfH6969eqpV69e2r17t4NVAwAAwClBHW6nTZumF198URkZGfr88881ffp0PfPMM5o1a5avz/Tp05Wenq6MjAxlZmYqLi5Offr00dGjRx2sHAAAAE4I6nC7ZcsW3XTTTRo4cKBatWqlm2++WX379tX27dslnZy1nTlzpiZNmqQhQ4aoffv2WrBggY4dO6bFixc7XD0AAABqW12nCzidHj166MUXX9RXX32lK664Qrt27dKmTZs0c+ZMSVJWVpays7PVt29f3z5ut1s9e/bU5s2bNXr06IDHLSwsVGFhoe99Xl5ewHZUVHZ/uE8ow5hAIIwLlMeYQHnVHRNV7ecyxpgaV3WeGWP0+OOPa9q0aQoJCdGJEyf09NNPa+LEiZKkzZs3q3v37tq3b5/i4+N9+917773as2eP3n333YDHTU1N1eTJkyu0T5gwQeHh4efnYgAAAFBjBQUFmjp1qo4cOaLo6OhK+wX1zO1rr72mRYsWafHixWrXrp127typlJQUxcfHa8SIEb5+LpfLbz9jTIW2U02cOFHjx4/3vc/Ly1NCQoLGjx9/2puFk//XNGPGDI0bN05ut9vpchAEGBMIhHGB8hgTKK+6YyIvL09Tp049Y7+gDrePPvqoJkyYoF//+teSpA4dOmjPnj2aMmWKRowYobi4OElSdna2mjZt6tsvJydHsbGxlR7X7XYHvImVtaMi7hXKY0wgEMYFymNMoLyqjomqjpug/kLZsWPHVKeOf4khISG+R4ElJiYqLi5Oa9as8W0vKirShg0blJycXKu1AgAAwHlBPXN744036umnn1aLFi3Url077dixQ+np6Ro1apSkk8sRUlJSlJaWpqSkJCUlJSktLU0REREaNmyYw9UDAACgtgV1uJ01a5b+8Ic/aMyYMcrJyVF8fLxGjx6t//mf//H1eeyxx3T8+HGNGTNGubm56tatm1avXq2oqCgHKwcAAIATgjrcRkVFaebMmb5HfwXicrmUmpqq1NTUWqsLAAAAwSmo19wCAAAA1UG4BQAAgDUItwAAALAG4RYAAADWINwCAADAGoRbAAAAWINwCwAAAGsQbgEAAGANwi0AAACsQbgFAACANQi3AAAAsAbhFgAAANYg3AIAAMAahFsAAABYg3ALAAAAaxBuAQAAYA3CLQAAAKxBuAUAAIA1CLcAAACwBuEWAAAA1iDcAgAAwBqEWwAAAFiDcAsAAABrEG4BAABgDcItAAAArEG4BQAAgDUItwAAALAG4RYAAADWINwCAADAGoRbAAAAWINwCwAAAGsQbgEAAGANwi0AAACsQbgFAACANQi3AAAAsAbhFgAAANYg3AIAAMAahFsAAABYg3ALAAAAaxBuAQAAYA3CLQAAAKxBuAUA4DzxeiWX6+TL63W6GuDiQLgFAACANeo6XQAAABeaqs7CntqvujO3Hk/1+gM4iXALAEA1RUZWf5/Y2Or1N6b65wDAsgQAAABYhJlbAACqKT+/av283h9nbA8eZKkBUBsItwAAVFNNQqrHQ7gFagPLEgAAAGANZm4BADhPPB6+GAbUNmZuAQAAYA3CLQAAAKxBuAUAAIA1CLcAAACwBuEWAAAA1iDcAgAAwBqEWwAAAFiDcAsAAABrEG4BAABgDcItAAAArEG4BQAAgDUItwAAALAG4RYAAADWINwCAADAGoRbAAAAWINwCwAAAGsQbgEAAGANwi0AAACsQbgFAACANQi3AAAAsAbhFgAAANYg3AIAAMAahFsAAABYg3ALAAAAaxBuAQAAYA3CLQAAAKxBuAUAAIA1CLcAAACwBuEWAAAA1iDcAgAAwBqEWwAAAFiDcAsAAABrBH243bdvn+644w41btxYERER6tSpkz788EPfdmOMUlNTFR8fr3r16qlXr17avXu3gxUDAADAKUEdbnNzc9W9e3eFhoZq1apV+uyzz/Tcc8+pQYMGvj7Tp09Xenq6MjIylJmZqbi4OPXp00dHjx51rnAAAAA4oq7TBZzOtGnTlJCQoHnz5vnaWrVq5ftnY4xmzpypSZMmaciQIZKkBQsWKDY2VosXL9bo0aNru2QAAAA4KKjD7cqVK9WvXz/dcsst2rBhg5o1a6YxY8bonnvukSRlZWUpOztbffv29e3jdrvVs2dPbd68udJwW1hYqMLCQt/7vLy8gO2oqOz+cJ9QhjGBQBgXKI8xgfKqOyaq2s9ljDE1ruo8Cw8PlySNHz9et9xyi7Zt26aUlBT9+c9/1p133qnNmzere/fu2rdvn+Lj43373XvvvdqzZ4/efffdgMdNTU3V5MmTK7RPmDDBd04AAAAEj4KCAk2dOlVHjhxRdHR0pf2COtyGhYWpa9eu2rx5s6/twQcfVGZmprZs2eILt/v371fTpk19fe655x7t3btX77zzTsDjBpq5TUhIUE5OzmlvFk7euxkzZmjcuHFyu91Ol4MgwJhAIIwLlMeYQHnVHRN5eXmKiYk5Y7gN6mUJTZs2Vdu2bf3a2rRpo6VLl0qS4uLiJEnZ2dl+4TYnJ0exsbGVHtftdge8iZW1oyLuFcpjTCAQxgXKY0ygvKqOiaqOm6B+WkL37t315Zdf+rV99dVXatmypSQpMTFRcXFxWrNmjW97UVGRNmzYoOTk5FqtFQAAAM4L6pnbcePGKTk5WWlpabr11lu1bds2zZ07V3PnzpUkuVwupaSkKC0tTUlJSUpKSlJaWpoiIiI0bNgwh6sHAABAbQvqcHv11Vdr2bJlmjhxop588kklJiZq5syZuv322319HnvsMR0/flxjxoxRbm6uunXrptWrVysqKsrBygEAAOCEoA63kjRo0CANGjSo0u0ul0upqalKTU2tvaIAAAAQlIJ6zS0AAABQHYRbAAAAWINwCwAAAGtUO9yOHDlSGzduPB+1AAAAAGel2uH26NGj6tu3r++xW/v27TsfdQEAAADVVu1wu3TpUu3bt0/333+/lixZolatWumGG27Q3//+dxUXF5+PGgEAAIAqqdGa28aNG+uhhx7Sjh07tG3bNl1++eUaPny44uPjNW7cOP373/8+13UCAAAAZ3RWXyg7cOCAVq9erdWrVyskJEQDBgzQ7t271bZtW82YMeNc1QgAAABUSbXDbXFxsZYuXapBgwapZcuWWrJkicaNG6cDBw5owYIFWr16tRYuXKgnn3zyfNQLAAAAVKrav1DWtGlTlZaWaujQodq2bZs6depUoU+/fv3UoEGDc1AeAAAAUHXVDrczZszQLbfcovDw8Er7NGzYUFlZWWdVGAAAAFBd1Q63w4cPPx91AAAAAGeNXygDAACANQi3AAAAsAbhFgAAANYg3AIAAMAahFsAAABYg3ALAAAAaxBuAQAAYA3CLQBcJLxeyeU6+fJ6na4GAM4Pwi0AAACsQbgFAACANQi3AAAAsAbhFgAAANao63QBAIBz40xfEjt1e1W+UObxnF09AOAEwi0AWCIysup9Y2PP3MeYmtcCAE5hWQIAAACswcwtAFgiP//0273eH2dsDx5k2QEAOxFuAcAS1QmrHg/hFoCdWJYAAAAAaxBuAQAAYA3CLQAAAKxBuAUAAIA1+EIZAFwkPB6eXQvAfszcAgAAwBqEWwAAAFiDcAsAAABrEG4BAABgDcItAAAArEG4BQAAgDUItwAAALAG4RYAAADWINwCAADAGoRbAAAAWINwCwAAAGsQbgEAAGANwi0AAACsQbgFAACANQi3ACDJ65VcrpMvr9fpagAANUW4BQAAgDUItwAAALBGXacLAIDzqapLDE7tV51lCR5P9eoBAJxfhFsAVouMrP4+sbFV72tM9Y8PADh/WJYAAAAAazBzC8Bq+flV6+f1/jhje/Agyw0A4EJFuAVgtZqEVI+HcAsAFyqWJQAAAMAahFsAAABYg2UJAKCTyxB48gEAXPiYuQUAAIA1CLcAAACwBuEWAAAA1iDcAgAAwBqEWwAAAFiDcAsAAABrEG4BAABgDcItAAAArEG4BQAAgDUItwAAALAG4RbAOef1Si7XyZfX63Q1AICLCeEWAAAA1iDcAgAAwBp1nS4AwIUj0BKDwkKpqChUXq9UUlKxX1WXJXg8Z18fAACEWwBVFhkZqNUt6XGlpQXeJza2asc2poZFAQBwCpYlAAAAwBrM3AKosvz8im2FhYV67rnn9PDDD8vtdks6uRShbMb24EGWHAAAag/hFkCVBQqpdetKYWHF8nik/8u2FfYh3AIAagvLEgAAAGANZm4BnHMeD18QAwA4g5lbAAAAWINwCwAAAGsQbgEAAGCNCyrcTpkyRS6XSykpKb42Y4xSU1MVHx+vevXqqVevXtq9e7dzRQIAAMAxF0y4zczM1Ny5c9WxY0e/9unTpys9PV0ZGRnKzMxUXFyc+vTpo6NHjzpUKQAAAJxyQYTb/Px83X777XrppZfUsGFDX7sxRjNnztSkSZM0ZMgQtW/fXgsWLNCxY8e0ePFiBysGAACAEy6IR4GNHTtWAwcO1PXXX68//vGPvvasrCxlZ2erb9++vja3262ePXtq8+bNGj16dMDjFRYWqrCw0Pc+Ly8vYDsqKrs/3CeUYUwgEMYFymNMoLzqjomq9gv6cPvqq6/qo48+UmZmZoVt2dnZkqTYst/5/D+xsbHas2dPpcecMmWKJk+eXKE9PT1d4eHhZ1nxxWHGjBlOl4Agw5hAIIwLlMeYQHlVHRMFBQVV6hfU4Xbv3r166KGHtHr16tOGTpfL5ffeGFOh7VQTJ07U+PHjfe/z8vKUkJCg8ePHKzo6+uwLt1hhYaFmzJihcePGyR3ot1Zx0WFMIBDGBcpjTKC86o6JvLw8TZ069Yz9gjrcfvjhh8rJydFVV13laztx4oQ2btyojIwMffnll5JOzuA2bdrU1ycnJ6fCbO6p3G53wJtYWTsq4l6hPMYEAmFcoDzGBMqr6pio6rgJ6i+UXXfddfrkk0+0c+dO36tr1666/fbbtXPnTl166aWKi4vTmjVrfPsUFRVpw4YNSk5OdrByAAAAOCGoZ26joqLUvn17vzaPx6PGjRv72lNSUpSWlqakpCQlJSUpLS1NERERGjZsmBMlAwAAwEFBHW6r4rHHHtPx48c1ZswY5ebmqlu3blq9erWioqKcLg0AAAC17IILt+vXr/d773K5lJqaqtTUVEfqAQAAQPAI6jW3AAAAQHUQbgEAAGANwi0AAACsQbgFAACANQi3gEW8XsnlOvnyep2uBgCA2ke4BQAAgDUItwAAALDGBfecWwCVLzk4tf1MyxI8nnNXDwAAwYJwC1yAIiPP3Cc29vTbjTk3tQAAEExYlgAAAABrMHMLXIDy8wO3e70/ztgePMjSAwDAxYdwC1yAqhJaPR7CLQDg4sOyBAAAAFiDcAsAAABrEG4BAABgDdbcAhbxeHjEFwDg4sbMLQAAAKxBuAUAAIA1CLcAAACwBuEWAAAA1iDcAgAAwBqEWwAAAFiDcAsAAABrEG4BAABgDcItLlper+RynXx5vU5XAwAAzgXCLQAAAKxBuAUAAIA1CLcAAACwBuEWAAAA1iDcAgAAwBqEWwAAAFiDcAsAAABrEG4BAABgDcItAAAArEG4BQAAgDUItwAAALAG4RYAAADWqOt0AUB5Xm/tn6e2zlnG46nd8wEAcLEg3CLoREbW/jljY2v3fMbU7vkAALhYsCwBAAAA1mDmFkEnP792zuP1/jhje/AgSwUAALAB4RZBx4mQ6fEQbgEAsAHLEgAAAGANwi0AAACsQbgFAACANQi3AAAAsAbhFgAAANYg3AIAAMAahFsAAABYg3ALAAAAaxBuAQAAYA3CLQAAAKxBuAUAAIA1CLcAAACwBuEWAAAA1qjrdAGAUzweyRinqwAAAOcSM7cAAACwBuEWAAAA1iDcAgAAwBqEWwAAAFiDcAsAAABrEG4BAABgDcItAAAArEG4BQAAgDUItwAAALAG4RYAAADWINwCAADAGoRbAAAAWINwCwAAAGsQbgEAAGANwi0AAACsQbgFAACANQi3AAAAsAbhFgAAANYg3AIAAMAahFsAAABYg3ALAAAAaxBuAQAAYA3CLRzl9Uou18mX1+t0NQAA4EJHuAUAAIA1CLcAAACwBuEWAAAA1gjqcDtlyhRdffXVioqKUkxMjAYPHqwvv/zSr48xRqmpqYqPj1e9evXUq1cv7d6926GKAQAA4KSgDrcbNmzQ2LFjtXXrVq1Zs0YlJSXq27evvKd882j69OlKT09XRkaGMjMzFRcXpz59+ujo0aMOVo5AvN7Ar6puBwAAOJO6ThdwOu+8847f+3nz5ikmJkYffvihrr32WhljNHPmTE2aNElDhgyRJC1YsECxsbFavHixRo8e7UTZqERk5Om3x8YGbjfm3NcCAADsFNThtrwjR45Ikho1aiRJysrKUnZ2tvr27evr43a71bNnT23evLnScFtYWKjCwkLf+7y8vIDtqKjs/tTsPrnP6pwITmc3JmArxgXKY0ygvOqOiar2cxlzYcyLGWN00003KTc3Vx988IEkafPmzerevbv27dun+Ph4X997771Xe/bs0bvvvhvwWKmpqZo8eXKF9gkTJig8PPz8XABUVBQasO3ZZx+VJD3yyDMKCyuu0CdQGwAAuLgUFBRo6tSpOnLkiKKjoyvtd8GE27Fjx+of//iHNm3apObNm0v6Mdzu379fTZs29fW95557tHfv3grLGsoEmrlNSEhQTk7OaW8WTt67GTNmaNy4cXK7azYTeyqvV2rc+ORxDh8ulMdz1odELTvXYwJ2YFygPMYEyqvumMjLy1NMTMwZw+0FsSzhgQce0MqVK7Vx40ZfsJWkuLg4SVJ2drZfuM3JyVFsZQs4dXLpQqCbWFk7KjpX96qkpPwxz/qQcAj//iAQxgXKY0ygvKqOiaqOm6B+WoIxRvfff7/eeOMNrVu3TomJiX7bExMTFRcXpzVr1vjaioqKtGHDBiUnJ9d2uQAAAHBYUM/cjh07VosXL9aKFSsUFRWl7OxsSVL9+vVVr149uVwupaSkKC0tTUlJSUpKSlJaWpoiIiI0bNgwh6sHAABAbQvqcDtnzhxJUq9evfza582bp5EjR0qSHnvsMR0/flxjxoxRbm6uunXrptWrVysqKqqWqwUAAIDTgjrcVuW7bi6XS6mpqUpNTT3/BQEAACCoBXW4hf08Hn6kAQAAnDtB/YUyAAAAoDoItwAAALAG4RYAAADWINwCAADAGoRbAAAAWINwCwAAAGsQbgEAAGANwi0AAACsQbgFAACANQi3FwivV3K5Tr68XqerAQAACE6EWwAAAFiDcAsAAABrEG4BAABgDcItAAAArEG4BQAAgDUItwAAALAG4RYAAADWqOt0ARer6j6r9tT+NXnOrcdT/X0AAAAuNIRbh0RG1nzf2Njq72NMzc8HAABwoWBZAgAAAKzBzK1D8vOr19/r/XHG9uBBlhkAAAAEQrh1yNmEU4+HcAsAABAIyxIAAABgDcItAAAArEG4BQAAgDUItwAAALAG4RYAAADWINwCAADAGoRbAAAAWIPn3F4gPB5+QhcAAOBMmLkFAACANQi3AAAAsAbhFgAAANYg3AIAAMAahFsAAABYg3ALAAAAaxBuAQAAYA3CLQAAAKxBuAUAAIA1CLcAAACwBuEWAAAA1iDcAgAAwBqEWwAAAFiDcAsAAABrEG4BAABgDcItAAAArEG4BQAAgDUItwAAALAG4RYAAADWqOt0AcHAGCNJysvLc7iS4FdYWKiCggLl5eXJ7XY7XQ6CAGMCgTAuUB5jAuVVd0yU5bSy3FYZlzlTj4vAd999p4SEBKfLAAAAwBns3btXzZs3r3Q74VZSaWmp9u/fr6ioKLlcLqfLCWp5eXlKSEjQ3r17FR0d7XQ5CAKMCQTCuEB5jAmUV90xYYzR0aNHFR8frzp1Kl9Zy7IESXXq1Dnt/wGgoujoaP7jBD+MCQTCuEB5jAmUV50xUb9+/TP24QtlAAAAsAbhFgAAANYg3KJa3G63nnjiCb7pCh/GBAJhXKA8xgTKO19jgi+UAQAAwBrM3AIAAMAahFsAAABYg3ALAAAAaxBuAQAAYA3CLapkypQpuvrqqxUVFaWYmBgNHjxYX375pdNlIYhMmTJFLpdLKSkpTpcCB+3bt0933HGHGjdurIiICHXq1Ekffvih02XBISUlJfr973+vxMRE1atXT5deeqmefPJJlZaWOl0aasnGjRt14403Kj4+Xi6XS8uXL/fbboxRamqq4uPjVa9ePfXq1Uu7d+8+q3MSblElGzZs0NixY7V161atWbNGJSUl6tu3r7xer9OlIQhkZmZq7ty56tixo9OlwEG5ubnq3r27QkNDtWrVKn322Wd67rnn1KBBA6dLg0OmTZumF198URkZGfr88881ffp0PfPMM5o1a5bTpaGWeL1eXXnllcrIyAi4ffr06UpPT1dGRoYyMzMVFxenPn366OjRozU+J48CQ40cOnRIMTEx2rBhg6699lqny4GD8vPz1aVLF82ePVt//OMf1alTJ82cOdPpsuCACRMm6J///Kc++OADp0tBkBg0aJBiY2P18ssv+9p+9atfKSIiQgsXLnSwMjjB5XJp2bJlGjx4sKSTs7bx8fFKSUnR7373O0lSYWGhYmNjNW3aNI0ePbpG52HmFjVy5MgRSVKjRo0crgROGzt2rAYOHKjrr7/e6VLgsJUrV6pr16665ZZbFBMTo86dO+ull15yuiw4qEePHlq7dq2++uorSdKuXbu0adMmDRgwwOHKEAyysrKUnZ2tvn37+trcbrd69uypzZs31/i4dc9Fcbi4GGM0fvx49ejRQ+3bt3e6HDjo1Vdf1UcffaTMzEynS0EQ+OabbzRnzhyNHz9ejz/+uLZt26YHH3xQbrdbd955p9PlwQG/+93vdOTIEbVu3VohISE6ceKEnn76aQ0dOtTp0hAEsrOzJUmxsbF+7bGxsdqzZ0+Nj0u4RbXdf//9+vjjj7Vp0yanS4GD9u7dq4ceekirV69WeHi40+UgCJSWlqpr165KS0uTJHXu3Fm7d+/WnDlzCLcXqddee02LFi3S4sWL1a5dO+3cuVMpKSmKj4/XiBEjnC4PQcLlcvm9N8ZUaKsOwi2q5YEHHtDKlSu1ceNGNW/e3Oly4KAPP/xQOTk5uuqqq3xtJ06c0MaNG5WRkaHCwkKFhIQ4WCFqW9OmTdW2bVu/tjZt2mjp0qUOVQSnPfroo5owYYJ+/etfS5I6dOigPXv2aMqUKYRbKC4uTtLJGdymTZv62nNycirM5lYHa25RJcYY3X///XrjjTe0bt06JSYmOl0SHHbdddfpk08+0c6dO32vrl276vbbb9fOnTsJtheh7t27V3hE4FdffaWWLVs6VBGcduzYMdWp4x81QkJCeBQYJEmJiYmKi4vTmjVrfG1FRUXasGGDkpOTa3xcZm5RJWPHjtXixYu1YsUKRUVF+dbJ1K9fX/Xq1XO4OjghKiqqwpprj8ejxo0bsxb7IjVu3DglJycrLS1Nt956q7Zt26a5c+dq7ty5TpcGh9x44416+umn1aJFC7Vr1047duxQenq6Ro0a5XRpqCX5+fn6+uuvfe+zsrK0c+dONWrUSC1atFBKSorS0tKUlJSkpKQkpaWlKSIiQsOGDavxOXkUGKqksrUv8+bN08iRI2u3GAStXr168Siwi9xbb72liRMn6t///rcSExM1fvx43XPPPU6XBYccPXpUf/jDH7Rs2TLl5OQoPj5eQ4cO1f/8z/8oLCzM6fJQC9avX6/evXtXaB8xYoTmz58vY4wmT56sP//5z8rNzVW3bt30wgsvnNUkCeEWAAAA1mDNLQAAAKxBuAUAAIA1CLcAAACwBuEWAAAA1iDcAgAAwBqEWwAAAFiDcAsAAABrEG4BAABgDcItAAAArEG4BQAAgDUItwAAALAG4RYALHHo0CHFxcUpLS3N1/avf/1LYWFhWr16tYOVAUDtcRljjNNFAADOjbfffluDBw/W5s2b1bp1a3Xu3FkDBw7UzJkznS4NAGoF4RYALDN27Fi99957uvrqq7Vr1y5lZmYqPDzc6bIAoFYQbgHAMsePH1f79u21d+9ebd++XR07dnS6JACoNay5BQDLfPPNN9q/f79KS0u1Z88ep8sBgFrFzC0AWKSoqEjXXHONOnXqpNatWys9PV2ffPKJYmNjnS4NAGoF4RYALPLoo4/q73//u3bt2qXIyEj17t1bUVFReuutt5wuDQBqBcsSAMAS69ev18yZM7Vw4UJFR0erTp06WrhwoTZt2qQ5c+Y4XR4A1ApmbgEAAGANZm4BAABgDcItAAAArEG4BQAAgDUItwAAALAG4RYAAADWINwCAADAGoRbAAAAWINwCwAAAGsQbgEAAGANwi0AAACsQbgFAACANf4/yMc8ia1kKoEAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xdata, xerror, ydata, yerror = np.loadtxt(\"fitdata.csv\", delimiter = ',', unpack = True)\n",
    "\n",
    "xdata[7] = 7.02\n",
    "ydata[7] = 68.8\n",
    "\n",
    "plt.figure(figsize = (8, 6))\n",
    "plt.title(\"Experimental data (fixed)\")\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"y\")\n",
    "plt.errorbar(xdata, ydata, xerr=xerror, yerr=yerror, color = 'b', linestyle = '', label = 'Data')\n",
    "plt.grid(color = 'grey')\n",
    "plt.legend(loc = 2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "6cea6244-36a9-47d5-93cc-0366941fda80",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def straight_line(xdata, c, m):\n",
    "    \"\"\"\n",
    "    Function for a straight line\n",
    "    - params   array of function parameters\n",
    "    - xdata    array of x data points\n",
    "\n",
    "    \"\"\"\n",
    "    \n",
    "    f = c + m*xdata\n",
    "    return f\n",
    "\n",
    "def straight_line_diff(xdata, c):\n",
    "    \"\"\"\n",
    "    Differential of function for a straight line\n",
    "    - params   array of function parameters\n",
    "    - xdata    array of x data points\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    df = c\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "d33ab8fd-bf5d-4bb9-813e-b1939ca80089",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.091 9.957]\n",
      "[[ 2.734 -0.457]\n",
      " [-0.457  0.088]]\n",
      "{'fvec': array([ 0.346,  1.701, -0.706, -0.693, -2.205,  0.351,  2.1  ,  1.321,\n",
      "       -0.81 , -1.112]), 'nfev': 7, 'fjac': array([[-13.728,   0.099,   0.061,   0.237,   0.33 ,   0.315,   0.273,\n",
      "          0.568,   0.513,   0.242],\n",
      "       [ -2.295,  -0.882,   0.16 ,   0.323,   0.274,   0.089,  -0.059,\n",
      "         -0.437,  -0.587,  -0.333]]), 'ipvt': array([2, 1], dtype=int32), 'qtf': array([ 3.189e-09, -1.874e-07])}\n",
      "Both actual and predicted relative reductions in the sum of squares\n",
      "  are at most 0.000000\n",
      "1\n"
     ]
    }
   ],
   "source": [
    "sigma = yerror # np.sqrt(yerror**2 + straight_line_diff(params, xdata)**2 * xerror**2)\n",
    "\n",
    "init_params = [0.0, 10.0]\n",
    "popt, pcov, info, msg, ierr = curve_fit(straight_line, xdata, ydata, sigma = sigma, p0 = init_params, full_output = True, method = \"lm\")\n",
    "print(popt)\n",
    "print(pcov)\n",
    "print(info)\n",
    "print(msg)\n",
    "print(ierr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b685555f-983b-41fa-b5af-9d3b1b9701f4",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fit succeeded\n"
     ]
    }
   ],
   "source": [
    "if ierr not in [1,2,3,4]:\n",
    "    print (\"ERROR: Fit failed with message {}\".format(msg))\n",
    "    print (\"Please check the data and initial parameter estimates\")\n",
    "else:\n",
    "    print (\"Fit succeeded\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "d2efa852-6bf7-4d66-9af5-6e0967269d7d",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "=== Fit quality ===\n",
      "chisq per point = \n",
      " [0.251 4.919 1.643 0.529 4.374 0.136 6.612 1.571 0.788 3.584]\n",
      "chisq =  24.408, ndf = 8, chisq/NDF =  3.0509, chisq prob = 0.0019575\n",
      "\n",
      "WARNING: chi2 probability for given degrees of freedom less than 0.05 .  Please check the data and initial parameter estimates\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import chi2 as stats_chi2\n",
    "\n",
    "# Get fitted parameters \n",
    "final_params = popt\n",
    "c = final_params[0]\n",
    "m = final_params[1]\n",
    "nparams = len(final_params)\n",
    "\n",
    "# Calculate chi2\n",
    "yfit = straight_line(xdata, c, m)\n",
    "sigma = np.sqrt(yerror**2 + straight_line_diff(xdata, c)**2 * xerror**2)\n",
    "chi2_array = (ydata - yfit)**2 / sigma\n",
    "chi2 = sum(chi2_array)\n",
    "npoints = len(xdata)\n",
    "reduced_chi2 = chi2 / (npoints - nparams) \n",
    "chi2_prob = stats_chi2.sf(chi2, (npoints - nparams))\n",
    "\n",
    "# Print chi2\n",
    "np.set_printoptions(precision = 3)\n",
    "print(\"\\n=== Fit quality ===\")\n",
    "print(\"chisq per point = \\n\",chi2_array)\n",
    "print(\"chisq = {:7.5g}, ndf = {}, chisq/NDF = {:7.5g}, chisq prob = {:7.5g}\\n\".format(chi2, npoints-nparams, reduced_chi2, chi2_prob))\n",
    "\n",
    "if reduced_chi2 < 0.25 or reduced_chi2 > 4:\n",
    "    print(\"WARNING: chi2/ndf suspiciously small or large.  Please check the data and initial parameter estimates\")\n",
    "\n",
    "if chi2_prob < 0.05:\n",
    "    print(\"WARNING: chi2 probability for given degrees of freedom less than 0.05 .  Please check the data and initial parameter estimates\")      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "09720caf-a3e6-43e7-abb9-75cce7a83270",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def prediction_band(x, func, popt, pcov):\n",
    "    # Central prediction\n",
    "    y = func(x, *popt)\n",
    "    \n",
    "    # Jacobian matrix for each x\n",
    "    # Use numerical derivative (can use autograd for complex cases)\n",
    "    eps = np.sqrt(np.finfo(float).eps)\n",
    "    J = []\n",
    "    for i, p in enumerate(popt):\n",
    "        dp = np.zeros_like(popt)\n",
    "        dp[i] = eps\n",
    "        J.append((func(x, *(popt+dp)) - func(x, *popt)) / eps)\n",
    "    J = np.array(J)  # shape: (n_params, len(x))\n",
    "    \n",
    "    # Variance of prediction\n",
    "    var = np.einsum('ij,jk,ik->i', J.T, pcov, J.T)\n",
    "    return y, np.sqrt(var)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "816be94e-7ca4-40b4-a75c-f065d4b83d90",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "yfit = straight_line(xdata, c, m)\n",
    "param_errors = np.sqrt(np.diag(pcov))\n",
    "x_fit = np.linspace(min(xdata), max(xdata), 500)\n",
    "y_fit, y_err = prediction_band(x_fit, straight_line, popt, pcov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "601d5634-f3e1-4626-8735-c072893e3631",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsAAAAIqCAYAAAAq8PwLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAC5aklEQVR4nOzdeXjU1fXH8ffMJBMggbAkEJA9hF0FZBWRWBVXwLpXrGLdKrgAtRa1P0VrAW0FlApqbV0qamlRAS2IS8EiyA6uyBJ2yMaSkHXW3x+XDGQyQICZzJLP63nySL6z3eschpObc8+1eL1eLyIiIiIitYQ13AMQEREREalJSoBFREREpFZRAiwiIiIitYoSYBERERGpVZQAi4iIiEitogRYRERERGoVJcAiIiIiUqsoARYRERGRWkUJsIiIiIjUKkqARURq2BtvvIHFYgn49fDDD7N9+3YsFgtvvPGG7zHLli1jwoQJHDp0KGzjPhMTJkzAYrFUujZjxoxKcxQRqSlx4R6AiEht9frrr9O5c+dK11q0aEGzZs1Yvnw56enpvuvLli3jqaeeYuTIkTRs2LCGRxoaM2bMICUlhZEjR4Z7KCJSyygBFhEJk+7du9O7d++At/Xv3z/or1dSUkK9evWC/rwiItFGJRAiIhHGvwRiwoQJ/Pa3vwWgXbt2vnKJxYsXH/c5Ro4cSVJSEt9++y1Dhgyhfv36XHzxxQA4HA6eeeYZOnfuTEJCAqmpqdxxxx3k5eVVeo4vvviCzMxMmjRpQt26dWndujXXXXcdJSUlACxevDjgOAKVcPhr27Yt33//PUuWLPHNp23btqf0/0lE5HRpBVhEJEzcbjcul6vStbi4qh/Ld911FwcOHGD69Om8//77NG/eHICuXbue8PkdDgfDhg3j3nvvZfz48bhcLjweD8OHD+d///sfjzzyCOeffz47duzgySefJDMzk9WrV1O3bl22b9/OVVddxaBBg/j73/9Ow4YN2bNnDwsXLsThcJzxSvIHH3zA9ddfT3JyMjNmzAAgISHhjJ5TRKS6lACLiIRJoDIHp9NZ5VrLli1p3bo1AD179qz2SqnT6eSJJ57gjjvu8F177733WLhwIXPmzOHaa6/1XT/33HPp06cPb7zxBvfddx9r1qyhrKyMP/3pT5x77rm++91yyy3Vnd4J9ezZk7p169KgQYOQlHuIiJyIEmARkTB566236NKlS6VrgVaAz8R1111X6fuPPvqIhg0bMnTo0Eqrzz169CAtLY3Fixdz33330aNHD+x2O/fccw+jRo1i0KBBtG/fPqhjExEJFyXAIiJh0qVLl+NugguGevXq0aBBg0rXcnJyOHToEHa7PeBj8vPzAUhPT+ezzz7jueeeY/To0RQXF9O+fXsefPBBHnrooZCNWUSkJigBFhGJUf59dwFSUlJo0qQJCxcuDPiY+vXr+/48aNAgBg0ahNvtZvXq1UyfPp0xY8bQrFkzbr75ZurUqQNAeXl5peeoSKJFRCKVukCIiESBig1ipaWlZ/Q8V199Nfv378ftdtO7d+8qX506daryGJvNRr9+/XjppZcAWLt2LYCvFvmbb76pdP958+ZVaywJCQlnPB8RkdOhFWARkShw9tlnA/DCCy9w++23Ex8fT6dOnSqt2FbHzTffzKxZs7jyyit56KGH6Nu3L/Hx8ezevZv//ve/DB8+nJ///Oe8/PLLfPHFF1x11VW0bt2asrIy/v73vwNwySWXAJCWlsYll1zCpEmTaNSoEW3atOHzzz/n/fffr/ac3nvvPf75z3/Svn176tSp45uniEgoKQEWEYkCmZmZPProo7z55pv89a9/xePx8N///pfMzMxTeh6bzca8efN44YUX+Mc//sGkSZOIi4ujZcuWDB482JeA9ujRg0WLFvHkk0+SnZ1NUlIS3bt3Z968eQwZMsT3fP/4xz944IEH+N3vfofb7Wbo0KG8++671aptfuqpp9i3bx933303hw8fpk2bNmzfvv2U5iMicjosXq/XG+5BiIiIiIjUFNUAi4iIiEitogRYRERERGoVJcAiIiIiUqsoARYRERGRWkUJsIiIiIjUKkqARURERKRWUR/gavJ4POzdu5f69esHPF5URERERMLL6/Vy+PBhWrRogdV6/HVeJcDVtHfvXlq1ahXuYYiIiIjISezatYuWLVse93YlwNVUcdzorl27aNCgQZhHE9nKy8uZMmUK48aNIyEhIdzDkQigmJBAFBfiTzEh/k41JgoLC2nVqtVJj4lXAlxNFWUPDRo0UAJ8EuXl5dSpU4cGDRroA0wAxYQEprgQf4oJ8Xe6MXGyclUdhVxNhYWFJCcnU1BQoAT4JDweD/n5+aSkpJyw/kZqD8WEBKK4EH+KCfF3qjFR3XxN0SVBZ7FYSE5O1mZB8VFMSCCKC/GnmBB/oYoJJcASdA6Hg8mTJ+NwOMI9FIkQigkJRHEh/hQT4i9UMaEEWERERERqFSXAIiIiIlKrKAEWERERkVpFXSCqSV0gqs/r9eJwOLDb7drIIIBiQgJTXIg/xYT4O9WYqG6+pj7ANcDpdOJ2u8M9jBrj8Xg4ePAgjRo1UhsbAY7GREpKinp7io/X66WgoICUlBQlOwIoJqSqUMWEEuAQKiwsJD8/n/Ly8nAPpUZVBGtBQYE+wASoHBN16tQhJSVFv0kRnE4nM2fOZPz48frBSADFhFQVqphQAhwihYWF7Nmzh6SkJFJSUoiPj681yaAamYs/j8dDXl4eDRs29P3dAJQEi4hIWCgBDpH8/HySkpJo2bJlrUl8K3g8HuLi4qhTp44SYAFMTMTHx1O/fn0aNGjA7t27yc/PVwIsIiJhoewkBJxOJ+Xl5bX6NJvaOm85voqYqDjVp7y8HKfTGeZRSbjZ7fZwD0EijGJC/IUiJtQFoppOpQtEWVkZ27Zto23bttStW7eGRigSPUpLS9m+fTvt2rWjTp064R6OiIjEiOrma1oBDqHaugrq9XopKytDP1tJBf+YqK1/N6Qyj8fDli1b8Hg84R6KRAjFhPgLVUwoAZag83q9HDhwQAmw+CgmJBCn08msWbNUCiM+ignxF6qYUAIsIiIiIrWKEmARERERqVWUAEtIxMUFr8NeZmYmFouFxYsXB+05peYFMyYkNlgsFlJTU1UTLj6KCfEXqphQAixBZ7Vaadq06Ul7ALdt2xaLxXLCr2nTph338YsXL2bChAlKjKNAdWNCahe73c6oUaPU9kp8FBPRpbi42PfvdXFxcUheI1QxoSUZCTqv10tJSQn16tWr1k9sGRkZNG3aNOBtZ511Fq1bt6ZTp07Uq1ev0m2LFy/mqaeeAswqsUSuU40JqR3cbjcbNmzg3HPPxWazhXs4EgEUE+IvVDGhBFiCzuv1UlBQQN26dauV7Dz22GOMHDnyuLffcMMNQRydhMOpxoTUDi6Xi/nz59OtWzclOwIoJqSqUMWEfh8pIiIiIrWKEmCJeIE2wVksFl/5w1NPPVWpbvhEq8kiIiIiKoGQkEhISAjp8w8cOJCdO3eya9cuWrVqRevWrX23dezYMaSvLacn1DEh0cdisZCenq6yGPFRTIi/UMWEEmAJOqvVSpMmTUL6GkuXLmXChAk89dRT/OpXv2LChAkhfT05MzURExJ97HY7t956a7iHIRFEMSH+QhUTEV8C8eWXXzJ06FBatGiBxWLhww8/rHS71+tlwoQJtGjRgrp165KZmcn3339f6T7l5eU88MADpKSkkJiYyLBhw9i9e3cNzsKP1wvFxdH5VY2jbL1eL4cPH672sbd33HFHwBZo6uwQO041JqR2cLlcLF68GJfLFe6hSIRQTIi/UMVExK8AFxcXc+6553LHHXdw3XXXVbn9ueeeY8qUKbzxxht07NiRZ555hksvvZSffvqJ+vXrAzBmzBjmz5/Pe++9R5MmTfjNb37D1VdfzZo1a8Kzy7SkBJKSav51g6GoCBITT3iXimQnMTHxjNqgnX322ac9TIkspxoTUju43W6WLFnCgAEDdFCKAIoJqSpUMRHx0XXFFVdwxRVXBLzN6/Uybdo0Hn/8ca699loA3nzzTZo1a8Y777zDvffeS0FBAX/729/4xz/+wSWXXALA22+/TatWrfjss8+47LLLAj53eXk55eXlvu8LCwurXLdarcTHx+N0OvF4PL77VvyU4vF4Kl2vWNn0eDyRv/R+HB6PB44zpwrHrvIde73i/v73GT9+PHfccUeV64Feu+L5rFar775er7fKdf/nsVqtAcfiP/bTvR5o7Me7fqIxnur1aJvTse+V1+vF4XD4nvN4f59sNhtxcXE4HI5Kzx8XF4fNZqtyPT4+HqvVWunvb8V1i8WCw+GodN1ut+P1enE6nZWuJyQk4PF4Kl23WCzY7Xbcbnel1YiK6y6XC7fbXWn+mtPx51Th2HlF+5xi8X2qyTkd+zqxMqdjr8fanI5VXl5OXFxcyOZU8d+Tzcn//8HxRHwCfCLbtm0jOzubIUOG+K4lJCQwePBgli1bxr333suaNWtwOp2V7tOiRQu6d+/OsmXLjpsAT5o0yddl4FhTpkyhTp06APTs2ZNhw4axYMEC1q1b57vPRRddRGpqKgUFBRw6dMh3PTk5mcTERPJLSnBv3uy73qhRI+rUqUN2dnalNz8lJQWbzUZOTk6lMTRr1gy3201+fr7vmsViIS0tjbKyMg4ePOi7brPZaNq0KcXFxb4kHkygNGnShMOHD1NUVOS7XrduXRo2bMihQ4coLS31XU9KSqJ+/focLC2l/PDhqnPKzw/464mcnJxKc0pNTcVms5Gdne37ECgoKMDr9eJ2u8nLy6s0p+bNm/v+Ah04cIDs7Gzi4uJo2rSp7y9oUVER2dnZJCQk0KRJE4qKijh8zBjr1atHw4YNKSwspKSkxHe9fv36Zk4HD1b6C3O8OTVu3Jg6deqccE7HSktLO+6cysvLOXDggO96xZxKSkooKCjwXY+1OeXk5JCQkEBiYiLl5eW8+uqrvtODjvf3afDgwWRmZjJ79my2bt3quz506FB69erFa6+9Vmk8I0aMoEOHDkyZMqXSB+R9991HcnIykydPrjSn8ePHU1BQwMyZM33X7HY7jz76KFlZWcyaNavS/5dRo0axYcMG5s+f77uenp7OrbfeytKlS1myZInvuuZ04jm1atUKgKlTp8bMnGLxfarpOVWIpTnF4vuUnp7Oz3/+c9/3zz//PHa7Pehzeumll4CjnxMnm9OxnycnYvFGUVGexWLhgw8+4JprrgFg2bJlDBw4kD179tCiRQvf/e655x527NjBJ598wjvvvMMdd9xR5SeCIUOG0K5dO1555ZWArxVoBbhVq1bk5ubSoEED4Pg/uTmdTnbt2kWbNm18yXLF+CNhFa46189kBbHi190V/5+ON/b27duzY8cO/va3vx13BdhqtZKZmcmSJUv4/PPPfXXBVquVJ598kqeffponnniCJ598MqRzOtn1QGOP9PepJufk8XgoLCykQYMGvp/cs7KyaNmypa87RLSvhMTi6k6o5+R2u/noo4+49NJLfSvC0T6nWHyfanJOTqeTTz/9lKuvvhqLxRITczr2ejS/T4WFhVXm5HQ6adasGQA7d+4kMTEx4JwSExNPe05FRUV8+umnvs+Jk81p//79NG3alIKCgoB5iO91j3tLFPFfhvd6vSetMzzZfRISEgK2bQp0/dhf5VU8N5jAtlqrFjsEunai68cbZ6DrFQnNmV4/1TH6X2/YsGHA+x37usc+tuL7E70n/v8/K/5ssVgqXQ/VnE52PRrfp5NdD9acbDYbjRo1qnLdbref9O9TheOdA3+868druxbousViCXjdarUGvG6z2QLuH4iLiwtYo6Y5Bb5utVorrSAdK1rnBLH3PkHNzSkhIaFSTMTCnI4Vze9TcnJywOeucGw7Un/HJranOqekpKSAnxOnMqdAorUUFTC/jgWq/Jo2NzfX9xNJWloaDoejUlmA/30kuDweD4cOHaqyqhhsdevWBahUqiGRqaZiQqKL0+lk3rx5VVahpPZSTIi/UMVEVCfA7dq1Iy0tjU8//dR3zeFwsGTJEs4//3wAzjvvPOLj4yvdZ9++fXz33Xe++0jwHVuXGirt27cHTCmMWuZEvpqICYkuHo+HdevW6Qcj8VFMRK6ioqIqX8fuUcrJyQl4n2P3GZ2OUMVExJdAFBUVsWXLFt/327ZtY/369TRu3JjWrVszZswYJk6cSEZGBhkZGUycOJF69epxyy23AGbJ/s477+Q3v/kNTZo0oXHjxjz88MOcffbZvq4QEp2GDBlCo0aNWLp0Ka1bt6Z9+/bExcVx+eWXM378+HAPT0REJGYknqQFamJi4knvE0kiPgFevXo1F110ke/7cePGAXD77bfzxhtv8Mgjj1BaWsqoUaM4ePAg/fr1Y9GiRb4ewGB2DsbFxXHjjTdSWlrKxRdfzBtvvBGeHsASNA0aNGDRokU88cQTrFixguXLl+PxeGjbtm24hyYiIiIRLKq6QIRTYWEhycnJJ91VCFBWVsa2bdto165dpS4QtYXX66WoqIikpKSTbkaU2sE/Jmr73xExXC4XS5cu5YILLtChBwIoJqJNcXExSUcO9ioqKgrJCvCpxkR18zVFlwSdxWKptAIvopiQQOLi4nytDUVAMSFVhSomonoTnEQmj8fD/v37tYlBfBQTEojD4eDtt9+u0s9Tai/FhPgLVUwoAZaQqO5RhFJ7KCbEn9frZevWrSc8Al1qF8WE+AtVTCgBFhEREZEz43feQqRTAiwiIiIip664+Oifd+wI3zhOgzbBSdBZLBaSk5PVAUJ8FBMSSFxcHEOHDtVuf/FRTESRXbvgiy+Ofh+i1rKhiglFmASdxWKJqmbYEnqKCQnEZrPRq1evcA9DIohiIgo4nbBhA6xaBWd4ylt1hComVAIhQefxeMjNzdWOf/FRTEggDoeDGTNmaMe/+CgmItyBA/DJJ7BkCdStCzVw8FSoYkIrwBISLpcr3EOQCKOYEH9er5e8vDzt+BcfxUSE8nph0yZYvhzy803im5AAZWUAhPI4o1DFhFaARURERCSwkhKz4rtggUl4O3Y0yS+QmJCA96GHKE1MJHHv3jAP9NRoBVhEREREqtq9G5YtMx0eWraEI8ceA1BYCC+9ZFaFAd59F669NjzjPA1KgCXoLBYLjRs31o5/8VFMSCDx8fGMGDGC+Pj4cA9FIoRiIkI4nfDNN7ByJTgckJFRucvD6tXw4otw6BDExcGQIfDYYyEZSqhiQgmwBJ3FYqFOnVBWBEm0UUxIIFarlQ4dOoR7GBJBFBMR4MABs6r744+QkmJWfiuUlcHrr5tyCIBWrWDcOHC7Q9YGLVQxoRpgCTqPx8O+ffu04198FBMSSHl5OZMmTdIx2eKjmAijio1u8+eb5LdtW2jS5OjtmzbB2LFHk99hw2DKFEhPD+mwQhUTWgGWkNAOXvGnmJBA1O5K/CkmwqCkxPT1Xb8e7Haz0a2iZM3thn/9C957DzwekxQ/9BD06FFjwwtFTGgFWOSIkpISHn74Ydq1a0d8fDwWi4WRI0cCMHLkSCwWC2+88UZQXzNUz1vTMjMzsVgsLF68ONxDERGRU7FnD3z0EaxYAamppuShIvnduxd+9zt45x2T/A4aZGp/azD5DRWtAIdTaakpLo8WdrtpfB0k27Zt47PPPmPlypWsXLmS77//HrfbzR/+8Ad+//vfB+11quvuu+/mnXfeoV69evTo0YOEhAQ6dux4wsesX7+eDz/8kB49enDNNdfUzEBFRETOlMt1dKNbWZnZ6FZx3LDXaw68+NvfoLwcEhPh17+GwYPDO+YgUgIcLqWlMHcuHDwY7pFUX6NGMHz4SZNgi8VCamrqSXf8v/DCC7zwwgvBHOFpO3jwIO+99x716tVj48aNtGrVqtLtzZs3p1OnTiQnJ1e6vn79ep566iluv/12JcAnUN2YkNolPj6e++67Tzv+xUcxUUMOHjQb3X74wZQ0nHVW5dv+8hdTEgFw9tkwZoxZHQ4kNxesVgjRexaqmFACHC4OhwmyunUhGnbHl5WZ8Toc1VoFtlVjN2hKSgpXX301ffv2pU+fPrz22mvMmTMnGKM9ZZs3b8bj8dC9e/cqyS/ApEmTmDRpUhhGFjuqExNSu1gsFpKTk/WDkfgoJkLM64UtW0xv37w8aNOmcg7y9dcm+S0sNAntbbfB0KEmwfXncJj+wPXrw89+ZlaQQyBUMaEEONzq1DG/WogGpaXVupvX6yU7O5u0tLQTBqx/mcN77713RsM7E6VH5lY3iCUeclR1Y0JqF4fDweTJkxk/fjwJR06WktpNMRFCpaWm3KFio1tGxtHEtqTElDt8+qn5vm1b+M1vTIIcSH6+aZeWkQH9+0PTpiEbdqhiQpvgJCplZWVxzz33kJ6eTp06dbBYLAG/xowZc8Ln2b59OxaLhczMTACWLFlS6fHbt28HAm9Wa9u2LXfccQcAb775ZqXHVTzfqdi0aRM33XQTTZs2pW7duvTs2ZO///3vAe9bWlrKu+++y80330ynTp1ISkoiKSmJHj168Mwzz1BcXBzwcW3btvXN6+uvv+aKK66gUaNGJCYmMmjQIL744ovjji8/P59Ro0Zx1llnUadOHTp16sQf/vAHnE7nKc9VRERq0N698PHHJgGu6O1bkfz++KMpcfj0U7P57dpr4fnnAye/TqdZQXY44KKL4IorQpr8hpJWgCXqfPPNNwwePJhDhw5hs9no3r07LpeLH3744ZRbbdWpU4eBAwdSUFDAd999R4MGDTj77LMr3X48ffr0wW63s3nzZpo2bUrGMb/+OfY5qmPz5s2MGTOGsrIyunXrxv79+1m/fj133nkn69ev58UXX6x0/zVr1nDLLbcQFxdHWloaXbp0oaCggO+//54NGzbwwQcfsHTp0uOuaH/00UeMGzeOBg0akJ6ezpYtW1i6dCmXXXYZn376aZUEPjs7m4EDB5KVlUVcXBzdu3enuLiYJ554gpUrV6rFmYhIJHK54NtvTeJbWgodOhzd6OZ0mtZmc+aYDg+pqabPb/fugZ9r/35TNpGeDuefD2lpNTePENAKsESdO++8k0OHDtGjRw+2bt3K+vXr+e677/j2229pc+Qn1ltvvZVVq1bxm9/85oTPlZaWxtKlS5k+fToAPXv2ZOnSpb6vtBP8Bf/Xv/7FY0eOfrziiisqPa7i+arrueeeo1evXuzatYs1a9awfft2/vWvfxEfH8/06dP5+OOPK92/VatWzJ49m4MHD7Jr1y5WrVrFpk2b2LVrF9dffz1r167lueeeO+7rjRs3jqeffpqcnBxWr15NXl4eI0aMwOVyMX78+Cr3HzVqFFlZWfTq1YusrCzWrVvHpk2b+Pzzz1myZAnLK86CFxGRyHDoEHz2GXz+uannTU8/mvzu2gWPPGL6+3o8ZjX3xRcDJ78uF2zdahLowYPh6qujPvkFJcASAhaLJWS1nitXrmT16tUAvPXWW76EF6Bbt26+rhLz588/7oa2SGSz2XjnnXdIPWaX7fXXX88DDzwAUCWZbdOmDTfccANJSUmVrqelpfHWW29ht9uZNWvWcV/v8ssvZ/z48b6NafHx8UybNo2EhARWrFjBwWO6k2zZsoUPP/wQMP/Pj/1/+rOf/YynnnrqpGUQoYwJiV52u53x48djt9vDPRSJEIqJIKjY6DZvHnz3nSllqPi3xeMxPX/HjjVJbf36ps/v2LGB9yMdPAibN5uSiWHDoG9fUz9cg0IVE0qAJSTcbndInvfLL78E4MILLwxYZjBs2DBatWpFQUEBS5YsCckYQuHaa68NuNo8atQoAL766qsqdb0ej4e5c+cyevRorrjiCgYNGsQFF1zApZdeisViYfPmzZSUlAR8vbvuuqvKtZSUFNq2bQuYGusKixYtwuv1cuGFF9KtW7eAz1WdD6ZQxYREL6/XS0FBgUpoxEcxcYbKyuCrr+A//4GiIrNJraKUb/9+eOopePVVU8Pbs6dZ9R04sOrzuFywfTscPmwOv7j66sqt0mpQqGJCCbAEndfrJS8vLyQfYDt27ACga9euAW+3WCx06dIFgB9//DHorx8qFWP21759exISEnC73WzdutV3/dChQ1xwwQVcc801zJgxg4ULF7J06VK++uorvvrqK9+Z6QeP02c6/Thntzc9spmhqKjId23Tpk0nHGP9+vU56yQfjKGMCYleTqeTmTNnaiOl+CgmzsC+fWaj27Jlprdvq1ZHN7otXQoPPADr1pkV3HvvhQkTzP38FRaaFeRmzUwLtAEDwtquNVQxoQRYokpZWRlgViuPpyKJO14nhEjU9Di7aCsOkAA4fPiw7/q4ceNYvnw5nTp1Ys6cOezZs4fy8nK8Xi9er9eXkB7vAyPxOK33rEc+LI9NVCuS4dTjNUEHmjVrdtzbREQkhNxu09ps3jxT25uRAQ0amNuKi2HKFHjuObMi3KEDTJsGV1119LjjY59n+3bT3mzAALPq27p1DU+m5qgLhESVBkf+Uu/bt++499m7d2+l+0aDvLy8gNcrVk7BrLQCuFwuZs+eDcDcuXPp1KlTpce4XC6ys7ODNraKOuPjjREgNzc3aK8nIiLVVFBgTnT7/ntzWmuLFkdv++47mDrVdG6wWuGGG+Cmm45uhDvW4cOwZ4+p9e3f3/QBjvE9G0qAJSRCtdmp4tfw3333XcDb3W63r/TheL+yD6ZgzfN45Rrbtm2jvLwcq9XqK1vIy8ujuLiYxo0bV0l+wfy/CWa9bceOHQHYuHFjwNuLiorYvXv3SZ9HG+AkEG12En+KiWrweiEryyS/2dlmpbai7aXTCW+/DR9+aO6XlgbjxkHnzlWfx+OB3btNzW/fvtC7d0QezhWKmFAJhASd1WqlefPmvl+nB9Pll1+OzWZjxYoVrKo4p/wYs2fPZt++fTRo0IABAwYE/fX9VfTZLa3mKXnHM2fOHHJycqpcnzFjBgADBw70lS1UvGZhYWHA1z1R+7PTMWTIEMBsQPzhhx+q3P7aa6/hcDhO+ByhjAmJXgkJCTz66KM68Ut8FBPVULHR7eOPzQpwRsbR5Hf7dpPsfvCBSX6HDIEXXgic/BYXw6ZNkJxsSiIuvDAik99QxYT+NZKg83q9lJWVhWTDU8uWLbntttsAuP3229myZYvvtlWrVvlOfhszZsxx61yDqX379r7XPl7Hhepwu92MGDGC/Px837UPPvjA10/4t7/9re96w4YN6datGy6Xi7Fjx/qST7fbzbPPPss///nPoP603KFDB4YPH47X6+X222+vtNq7ePFiJkyYQHx8/AmfI5QxIdHL4/GwZcsWPB5PuIciEUIxcRLZ2abDw/LlZgNb69amvMHjMUnvuHGwY4dJah9/HO6//2hyXKFi1Tc7G3r1guHDTW1whP6WLlQxoQRYgs7r9XLgwIGTJjtfffUVKSkpvq/33nsPgEmTJlW6vmvXrkqPmzZtGn379uXHH3+kc+fO9OjRg27dutG3b19yc3O56qqr+P3vfx+y+R2rV69eZGRksG3bNlq3bs35559PZmbmSY9g9vfb3/6W1atX06pVK3r37k27du249tprcTgcjBo1iqFDh1a6/6RJk7BYLLzyyis0b96cPn36kJaWxvjx43n88cdp3rx5EGdpVqLbtm3L6tWrad++Pb169aJTp05cdNFFXHDBBSddba9uTEjt4nQ6mTVrlnb8i49i4jjcbvjmG7PRbccOc6hFxT6XvDz4v/+D118/WsowfTr061f1eUpKzKpvvXrmGOOLLjK9gCNYqGJCNcDhdqSrQcQLwTidTif79++vcr2kpKTSaqp/PWuDBg1YsmQJ06ZN491332XTpk3YbDb69u3LHXfcwd133+074CHUrFYrH3/8MY899hhffvklK1euPK36244dO7Jy5Up+//vfs3jxYgoLCzn33HMZPXp0wJ69Q4cOZcGCBTz99NOsW7eOn376iW7dujFt2jRGjBjBW2+9FYzp+bRo0YKVK1fyxBNPMHfuXH744Qdat27N008/zfjx47n00kuD+noiInJEQQF8/bXZ6JacbFZrwZQ4LFkCr7xiyhnq1IG77oJLL626muv1mjZpJSVw7rkmOU5Orvm5RBCLV0sy1VJYWEhycjIFBQUn7S5QVlbGtm3baNeuHXWO1zuvtBTmzjWnrESLRo3Mr0r8f53ix+PxkJ2dTVpammo+BagaE9X6OyIxr7y8nMmTJzN+/HjVfAqgmKjE64Vt20xf3337zIluFf/+Hj4MM2ea/r4AnTqZ09yO7QJRobQUdu6ElBST+HbqdLQ/cBQ41Ziobr6mFeBwqVvXJJMn2TwUUez2kya/FeICtVmRWk0xIf4q+lyrQ4hUUEwcUVYGa9aYL5sNOnY8mrSuW2c2th04YG67+Wa4/nrz52N5vZCTYw626N7dJL+NGtX8XM5QqGJCK8DVFPQVYJFaTH9HRESOIyfHdHnYssWs6FaUKpSXw1tvwfz55vuzzjKb3jIyqj5HebnpCNGokUl8u3SpmiDHqOrma9GzBi5Rw+v1UlxcrA1P4qOYkEDcbjdr164Nat9qiW61Oibcbvj2W7PRbft2U+tbkfxu3WqS3Yrk98orzYlugZLfvDzz+M6dzW+au3eP6uQ3VDGhBFiCzuv1UlBQoGRHfBQTEojL5WL+/Pm4XK5wD0UiRK2NicJC+Pxz+PRTs4EtIwPi401SPHs2PPywOea4cWN48kn49a/Bvx7W4YDNm81jLrkELr/c1P1GuVDFhIryRERERMLBf6Nb69amRRmYPr1TpkDFKZznnw+jRh1tf3as/HzYv9/UCvfvD02b1twcopQSYBEREZGaVl4Oa9fC6tVmg1vFRjevFz77DF57zXRwqFsX7r3X9Oz13wjmdJq+wPXqmdvPOcesHMtJKQGWkKj17WukCsWE+LNYLKSnp2vHv/jUmpjIzTUtzLZuhbQ0aNjQXC8ogJdeMn1/Abp2Ne3NmjWr+hwHDpjnSU+HAQMgyAcgRYpQxYS6QFSTukCIBI/+johIreR2w48/mgS3sBDatj26YrtqFbz4okmC4+JgxAi45pqqG9hcLrPqm5BgjjLu0aNqPXAtpi4QEjZer5fDhw9rw5P4KCYkEJfLxeLFi2vfhic5rpiOiYqNbosWmTKHio1uZWUwYwb84Q8m+W3dGp5/Hq67rmrye+iQaY/WsiUMHWpanMV48huqmFACLEGnZEf8KSYkELfbzZIlS2pnyysJKGZjYts208JswwbTvzctzVz/6Sd46CFYuNB8P3y42fjWrl3lx7vdprVZYSFccAFcfbVJgmuBUMWEaoBFREREQqFio9uaNeb7io1uLpdpbzZ7Nng8pl3ZQw/BuedWfY7CQtizxxyF3L+/+a+cMSXAIiIiIsGWm2vam23ebDaxVRxDvGePWeXdvNl8P3iw6fKQlFT58W437N5tEuT+/aF376Mt0uSMKQGWkKinv6TiRzEh/qxWKz179sRqVTWeGDEREx6P2ei2fLlZvU1PN7W+Xi8sWAB//7s5tCIxEe67Dy68sOpzFBWZgy9atjTJb7t2VVug1RKhigl1gagmdYEQCR79HRGRmHT4sOnw8O23UL++Wfm1WODgQdPhoaIU4txzTcmD/0ltHo9ZIXY6zX169666MiwnpC4QtVBxcTEWiwWLxUJxcXHYxuHxeDh06BAejydsY5DIopiQQJxOJ/PmzcPpdIZ7KBIhojomtm83G93Wrz+60c1iMSvB999vkt/4eLjrLnjqqarJb3ExbNpkTnq78kpTGqHkN2QxoQRYQqKkpCTcQ5AIo5gQfx6Ph3Xr1ukHI/GJyphwOMyq78cfm+OIMzJMrW5JCbzwAkyaZFaG27WDqVNh2DCzEa5Cxarvvn3Qs6e5PSOj1pY8+AtVTKgGWEREROR05OWZjW6bNlXe6PbDD2ajW26uSWSvvRZuuaXqMcUlJbBzJzRtCpmZR7tESMgpAZawadu2LTt27PB9b7FYSEpKIjk5mc6dO9OvXz9uueUWunbtGrTXPHToENOmTaNhw4aMGTMmaM8rIiK1iP9Gt/btwW43tbvvvAPvv282vTVtao4y7tat8uO9XsjONpvdzjnHHGhRcRyy1AglwBJ0FouF+vXrV/vc7oyMDJo2bQqYzVH5+fl89tlnfPbZZ/zxj3/kuuuu45VXXqFJkyZnPLZDhw7x1FNP0aZNGyXANehUY0JqB5vNxuDBg7H5n3YltVZUxMThw7BihdnolpQEHTqYVd6dO80Jbtu2mftdfDHcfXfV1mVlZeYo48aN4fLLoVOnqie+iU+oYkIJsARdRbJTXY899hgjR46sdC0/P59Zs2bxzDPPMGfOHL7//nu+/vprkpOTgzxaqQmnGhNSO8TFxZGZmRnuYUgEifiY2L7dlDzs2WOOLK5Xz6wGz58Pb75pVoDr14fRo+H88ys/1us1JREFBdC1q2lv1rhxWKYRTUIVEyo0kaDzeDzs37//jArWU1JSeOihh1i9ejXNmzdn48aNWrGNYsGICYk9DoeDt99+G4fDEe6hSISI2JhwOMyqr/9Gt/x8ePJJeO01k/yedx5Mn141+S0vhy1bzJ+HDDFfSn6rJVQxoQRYQqK8vDwoz9OmTRtmzJgBwNtvv82uXbt8t2VlZfHss8+SmZlJq1atSEhIIDU1lcsvv5yPP/64ynONHDmSdkfOV9+xY4evZVzFV4XS0lLeffddbr75Zjp16kRSUhJJSUn06NGDZ555Jqwt5qJZsGJCYofX62Xr1q2oHb1UiMiYyMszB1h8+aVZ3W3b1pQs/O9/8MADsGGDqf/99a/hiSeqJrZ5eWblOCMDhg+Hs8+GOP0CvrpCFRN6ByTiDRs2jBYtWrB3714WLVrEnXfeCcDEiRP529/+RlJSEi1atOCcc85hz549fPLJJ3zyySdMnjyZ3/3ud77n6dixI71792b16tUkJCTQu3fvgK+3Zs0abrnlFuLi4khLS6NLly4UFBTw/fffs2HDBj744AOWLl1K3bp1a2T+IiISOsXFxSQd6bdbVFREYmKiucHjgY0bzUa3goKjG92KiuCVV2DJEnO/jAyz0a1ly8pP7HCYWt/ERFMP3L171S4QEjZKgKNEdVYdj71PdVcpfX/RI5jVamXAgAHMmTOHVatW+RLg6667jrvvvpu+fftWWsH93//+x0033cTjjz/O9ddfT3p6OmBqjW+55RbatWtHWloaS5cuDfh6rVq1Yvbs2VxxxRW+D0WA7OxsHnjgAf7973/z3HPP8eSTT4Zw1iIiEjZFRUdPdDt2o9uGDaa3b36+aVd2443my39Fd/9+s/KbkWFqfdPSwjMPOS4lwFEi6RRPg2nWrFm17heKXzNZLBaSk5ODuuO/VatWAOTm5vquXXHFFQHvO2jQIP7whz9w11138c9//pPHHnvslF6rTZs2tGnTpsr1tLQ03nrrLebNm8esWbOUAJ+CUMSERL+4uDiGDh1KnH4dLEdEREzs2GE2uu3eDa1amRVchwP+8Q+YO9fcp3lzs+rbuXPlxzqdphtEQgJcdJFpcWa31/wcYkioYkKfOhJ0Fosl6CvLFc93+PDhStfz8vJ45513WLFiBbm5uZSVlQFQUFAAwIYNG07r9TweD/Pnz2fRokVkZWVRVFTk+2HBYrGwefNmSkpKqOff3kYCCkVMSPSz2Wz06tUr3MOQCBLMmKjub0Ir/fb0yy9Nf1+32yS/NhvWjRtJmD4d65E9KM5LLiH+nnugTp3KT3TwIOTkmFKJAQOgRYugzKO2C9XnhBLgKFFUVHTS+xQXF/tWfnNycsKWcHg8HvLz80lJScEapBNtKubfoEED37VFixZx4403+pLdQA4cOHDKr3Xo0CGuvPJKli9ffsL7HTx4UAlwNYUiJiT6ORwOXnvtNe666y7sWiUTghsTp/qbU4BmV17p+7MVGAc8c+TPOcCvgP989hneBx88+iCXC3btMmUQgwZBjx5Vk2M5baH6nNC/RFEiMTGxWl+ne/9gc7lcQX2+nTt3AvgOzDh06BA333wzBQUF3HbbbXz99dccPHgQt9uN1+vl008/BcDpdJ7ya40bN47ly5fTqVMn5syZw549eygvL8fr9eL1ejnrrLNO+7lrs2DHhEQ/r9dLXl5eZO34l7CKlJhoDXwB/AlIAD4EugP/8b9jQQFs3mxqfIcONfW+Sn6DKlQxoRVgiXgej8e3Gtu3b18AFixYwMGDBxkwYABvvPFGldrSY9ulnQqXy8Xs2bMBmDt3Lp06dapye3Z29mk9t4iI1JyT/ua0qAhWraJ47VqaHdnTkfPmmySvXIn973/HUlqKt04dHCNHculFF7H92H9n3G6z6muxwMCB0KsXqDNQVFECLBHvww8/JDs7m/j4eIYMGQLA9u3bARgwYEDAjVXHq/092SasvLw8iouLady4cZXkF+C7777D7Xaf4gxERKSmnfC3nMdudGvdGoDGQOOZM4n7+mtzn86dsYwbR0JaGgnHPvbwYXPccatWpta3TRuTCEtUUQIsQWexWGjcuHFQdvzv2LGD+++/H4DbbrvNV35Q0YM3JyenymP279/P3/72t4DPV/G40tLSE95eWFhIaWlplV6/zz333GnMQoIZExI74uPjGTFiBPHqjSpHhDwmHA5Ytw7WrDGruBkZ4HRyGfB3MMmvzQa/+AVcd535cwW32xyB7HZDv37Qu7fpECEhFaqYUA2wBJ3FYqFOnTpnlOzk5+fz4osv0rt3b/bt20fXrl2ZMmWK7/ZBgwYBMHv2bD777DPf9X379nHdddcdt940NTWV+vXrk5uby48//ljl9oYNG9KtWzdcLhdjx471Hb3odrt59tln+ec//6nNOqchGDEhscdqtdKhQwdtjBSfkMZEfj588ok50S0pyZzo5nJh/9vfWAi0ADxnnQV//rPp7Xts8ltUBJs2QXIyXHUVXHihkt8aEqqY0KeOBJ3H42Hfvn14PJ5q3X/ixIlccMEFXHDBBfTp04d27dqRmprKQw89RH5+PjfccAP/+9//KnWAOO+887j++utxOp1ceumlZGRk0LNnT1q3bs3atWuZPHlywNeyWCzccMMNAPTq1Ys+ffqQmZlJZmam7z6TJk3CYrHwyiuv0Lx5c/r06UNaWhrjx4/n8ccfp3nz5qf/P6eWOtWYkNqhvLycSZMm6Zhs8QlJTHg8prXZvHkmiW3fHho1MpvXxowh/pNPAHgRKH32WThyeJLvsbt3m/ZmvXubo4zT01XyUINC9TmhEggJiVPZrbl582Y2b94MmLY1DRs25JJLLqFfv36MGDGCLl26BHzcrFmz6NKlC//4xz/YsWMHTZo04frrr2fChAns27fvuK/3wgsvUL9+febOncuGDRuqdHMYOnQoCxYs4Omnn2bdunX89NNPdOvWjWnTpjFixAjeeuutas9Njgr3rm6JTBW/ZRGpENSYKC6GFSvgm2/MJrWMDJPU/vOf8N57ppyhcWN46CEe7Nmz8mNLSsxGt9RUc5Rxhw7m9DepcaH4nFACLGFTsZHtdNntdp5++mmefvrpKrd16tTpuAlXUlIS06ZNY9q0acd97ssuu4zLLrss4G1nOm4REakBO3fCV1+ZJLZ1a1OysG8fTJ0KGzea+wwcCKNGQf36Rx/n9Zr7lZSYk9z69TOlDxJTlACLiIhI7HA4YP16WL3arPB27GhWbj/5BP72Nygrg3r14N57ITOzcjlDWZnpENGkianz7dRJq74xSglwDElMTIyIXzNbLBZSU1O14Ul8FBMSSHx8PPfdd5+6QIjPGcdEfj4sX25WeJs2NeUNhw7BX/4CK1ea+5x9NowZY0obKni9ps738GHo3h369jWPlbAL1eeEEmAJCduxu2dFUExIVRaLheTkZP1gJD6nHRMeD/z0k0l+Dx40G93sdpP0Tp9uTmyLi4Nf/tJsZDt2Vbe83Kz6JifDkCHQpUvlDhASVqH6nNC6vgSd1+slOzs7IlajJTIoJiQQh8PB5MmTtRFOfE4rJoqLYfFiWLgQnE6z0c3tNqu+zzxjkt82beD55+HnP6+c/OblwfbtptThmmvM6q+S34gSqs8JrQCLiIhIdAq00W3jRpgyBbKzTX3vNdfAiBFmRbiCw2FWfZOS4JJLoFs3UClOraIEWERERKKL02lOdDt2o5vXC2+/Df/+tymJSEmBsWNNze+x9u83tcIdOpijjJs1C88cJKyUAIuIiEj0qNjo9tNPZiNb48bmsIopU2DLFnOfzEy45x6zwlvB6TQrxnXqwEUXmcRYJ3vWWhavivKqpbCwkOTkZAoKCiqdSBZIWVkZ27Zto23bttStW7eGRhg5vF4vXq8Xi8WizS0CVI2J0tJStm/fTrt27ahTp064hydh4vV6cTgc2O12fVYIcJKY8N/o1ratKVv4z3/g9ddNWUNSkunre8EFlR978KDp8pCeDv37Q4sWNTYnOTOn+jlR3XxNK8AhUHFetdvtDvNIwsftdhMXp/CSo46NiYq/G8E+212ii9frpaCggJSUFCXAApwgJoqLTUeH9etND9+MDDhwAF580ZRCAPToAQ89ZHr4VnC5zKpvfDwMGgQ9e0JCQk1OSc5QqD4n9K9PCMTHx2Oz2SgtLQ33UMLC6/WSl5enHf/i4x8TpaWl2Gw29X+t5ZxOJzNnzqxyHLnUXgFjYudOmDcPVq2CtDSzert8OTz4oEl+7XZT7jBhQuXkt6AANm+G5s1h6FCz8qvkN+qE6nNCS3QhYLFYqFevHgUFBTRu3Fj9T0WO4Xa7KSgooF69elr1E5HjC7TRrazMHGX83/+a+7RvD+PGmQ4QFdxu0xXCYjFHHffqBbWwHFFOTAlwiDRt2pTt27ezY8cOGjduTEJCQq35x97j8eByuSgrK9OvuAUwMeF0OikoKODgwYN4PB6aNm0a7mGJSKQ6cMAkvsdudPvuO5g2DXJzTS/f666Dm2+u3L6ssBD27IFWrUyHhzZtKh91LHKEEuAQsdvttGzZkvz8fPbt2xfu4dQor9dLYWEhRUVFtSbplxOriIni4mKSkpJIS0vDrt3XAooDqczjwR4XBx9/bI4wbtvWJLtvvAEffGBanTVrZtqbde1a6XHs3m1Wf/v1g969TU9giQmh+JyI+i4QLpeLCRMmMGvWLLKzs2nevDkjR47k97//vW/10ev18tRTT/Hqq69y8OBB+vXrx0svvUS3bt2q/Tqn0gUi0BhdLtcpPUYkFsXFxWlzpIgE5r/RrXlzU/87ZQps22buc8klcNdd5vYKRUUm+W3e3Kz6tm+vVd9arNZ0gXj22Wd5+eWXefPNN+nWrRurV6/mjjvuIDk5mYceegiA5557jilTpvDGG2/QsWNHnnnmGS699FJ++ukn6tevH/Ix1rZ/9D0eD1lZWbRv314lEAJUjgmRCvqsEJ+dO2HZMjw7dpDVpAntU1KwzpsHb71laoEbNID77zcb2Sp4PLB3L5SXw3nnQZ8+UAP/pkvNCtXnRNR/4ixfvpzhw4dz1VVX0bZtW66//nqGDBnC6tWrAbP6O23aNB5//HGuvfZaunfvzptvvklJSQnvvPNOmEcfm5xOJ7NmzdLObvFRTEggigvB6TSrvh99BHl5ODt0YP4PP8ATT8Df/mZu790bpk+vnPyWlJgOD/XqwZVXmoMvlPzGpFB9TkT9suQFF1zAyy+/zKZNm+jYsSMbNmxg6dKlTJs2DYBt27aRnZ3NkCFDfI9JSEhg8ODBLFu2jHvvvTfg85aXl1NeXu77vrCwsMp1q9VKfHw8TqcTj8fju6/NZiMuLg6Hw1GpFVhcXBw2m63K9fj4eKxWa6XXq7husVhwOByVrtvtdrxeb5VgSEhI8G02qmCxWLDb7bjd7kplGBXXXS5XpX7FwZjTsfOIlTlVjD2W3qeanFPF85eXl8fMnGLxfarpOVU4dl7RPqdYfJ9CNqf9+7GvWYN740ZcKSnQqBHuxYu57+WXsZaX401IwHXHHXiGDAGLBavHQ7zFgnPvXjwlJaYGuHdvbI0bE2e1RsacYvF9ipA5Vfz3ZHPy/39wPFGfAP/ud7+joKCAzp07Y7PZcLvd/PGPf+QXv/gFANnZ2QA08zvru1mzZuzYseO4zztp0iSeeuqpKtenTJniO7mqZ8+eDBs2jAULFrCuohE3MHjwYDIzM5k9ezZbt271XR86dCi9evXitddeIy8vz3d9xIgRdOjQgSlTplR6Q++77z6Sk5OZPHlypTGMHz+egoICZs6c6btmt9t59NFHycrKYtasWb7rqampjBo1ig0bNjB//nzf9fT0dG699VaWLl3KkiVLfNeDNadj/3/Fypxi8X2q6TlNnTo15uYEsfc+1dScWrVqBZi4iJU5xeL7FLI5xcUxKiGBDQ0b8ulPP3Hlxx9z9nffAeDJyGDVrbey0OuFb74xc2rYkGFOJwuAdcXFpi/wqlWRNadYfJ/CPKeXXnoJOPo5cbI5Hft5ciJRvwnuvffe47e//S1/+tOf6NatG+vXr2fMmDFMmTKF22+/nWXLljFw4ED27t1L8+bNfY+7++672bVrFwsXLgz4vIFWgFu1akVubq6vqFo/uQWek8Ph4K233uLuu++uchhGtM6pYuyx9D7V5JxKSkp48803uf3220lISIiJOcXi+1TTc3K5XLz66qvcfvvtvl3e0T6nWHyfgjqnQ4dgzRr49lssdetiP+ss3OvXY33xRSz79+O1Wll10UX0uO8+rHFxuL1e0/khNxdrURHxXbviPO88PMnJkTOnWHyfImhOhw8f9v37YbfbTzqn/fv307Rp05Nugov6BLhVq1aMHz+e0aNH+64988wzvP3222zcuJGsrCzS09NZu3YtPXv29N1n+PDhNGzYkDfffLNar3MmXSBERERqvSMb3di50/TpjY83m9wqVhVbtDCHWnTsePQx5eWwYwckJ5sa4C5dQIdLyQlUN1+L+k1wJSUlVXYF2mw2308d7dq1Iy0tjU8//dR3u8PhYMmSJZx//vk1Otbawu12s3bt2ko/EUrtppiQQBQXtYTfRjc6doScHJPsViS/V1wB06bh7tCBtfn5uD0ec9/t26FTJ7jmGujeXclvLRSqz4morwEeOnQof/zjH2ndujXdunVj3bp1TJkyhV/96leAWaIfM2YMEydOJCMjg4yMDCZOnEi9evW45ZZbwjz62ORyuZg/fz7dunXTMdACKCYkMMVFLZCfD8uXHz3RLTnZHGjxzjvgckHDhvDgg6bTA+Byu5m/cyfdDhzAVr8+XHyxSXyPPe1NapVQfU5EfQI8ffp0/u///o9Ro0aRm5tLixYtuPfee3niiSd893nkkUcoLS1l1KhRvoMwFi1aVCM9gEVERGodj8ckvcuXm2ON27aFgwfh8cfhhx/Mffr3h9GjTVJc4cAB89+2beGCC8ypbyIhEPUJcP369Zk2bZqv7VkgFouFCRMmMGHChBobl4iISK1UXAwrVsCGDaZPb0YG/Pe/8OqrUFoKdevC3Xeb1d2KE9ucTlMbXLHSO2SI+vpKSEV9AiyRx2KxkJ6ejkVHUcoRigkJRHERg3buhK++gl27zEY3jwcmT4avvza3d+kCY8dCWtrRxxw8aGqC27fHct55pC9diiUhITzjl4gTqs+JqO8CUVPUBUJEROQ4nE5Ytw5Wrwa32yS/69bBiy/CoUMQFwe33AI///nRjWwul0mUbTZzlHGPHnCkz77I6ao1XSAk8rhcLhYvXlypX6DUbooJCURxESPy82HhQvjyS0hMNKu7r74KTz9tkt9WreBPf4Lrrz+a/BYUwJYt5r5Dh5p64Dp1FBNSRahiQgmwBJ3b7WbJkiVqbSQ+igkJRHER5Twe+PFHmDcPNm2Cdu1MMjx2LCxYYO4zdChMmQLp6eZ7t9u0Njt4EAYMMLe3bu17SsWE+AtVTKgGWERERE5NUZHp7Vux0a19e/j3v+G990xi3LgxPPQQHHMAFYcPw5490LKlWfFt2/boJjiRGqYEWERERKrv2I1urVubcobx480qMMCgQfDrXx/t4uDxmMTX5YI+fcxXYmL4xi+CEmAJAavVSs+ePauc0Ce1l2JCAlFcRBmHA9avP7rRLSMDPvsM/vY3c2RxYqJJfC+88OjKbnGxSZSbN4d+/aBDhxOu+iomxF+oYkJdIKpJXSBERKTWys+HZcvM4RZNm5ok9i9/gVWrzO1nnw1jxpjT3sCs+u7bZ/r+nn029O0L+rdTaoC6QEjYOJ1O5s2bh9PpDPdQJEIoJiQQxUUU8HjMyW1z58LmzabWd9MmeOABk/zGxcGvfgV/+MPR5LekxNy3bl244gr42c+qnfwqJsRfqGJCCbAEncfjYd26dXg8nnAPRSKEYkICUVxEuKIic4LbJ5+YkoezzoJXXoGJE6Gw0GximzoVrrkGrFbwes2q7549ZtV3+HDo3NncVk2KCfEXqphQDbCIiIhUtmOHKXnYvdv08d25E558ErKzTfnDz38OI0YcPbq4vNy0N2vcGC67zCS+FT1/RSKQEmARERExHI7KJ7q1bQv/+hfMmWPKIVJTTZ/f7t2PPiY31/T17dLFtDdr0iRswxepLiXAEnQ2m43Bgwdj00//coRiQgJRXESYvDxYvtzU+DZtakogHn0Utm41t190Edxzz9EWZg6HWSmuXx+GDIGuXU1N8BlQTIi/UMWEukBUk7pAiIhITPJ4YONGk/wWFJiSh88+gzfeMElu/fowahQMHHj0Mfn5cOCAaYXWv79JmEUigLpASNg4HA7efvttHA5HuIciEUIxIYEoLiJAURF88cXRjW6NGplNbq++apLfnj3hxRePJr9Op1kRdjggM9N0eQhi8quYEH+higmVQEjQeb1etm7din65IBUUExKI4iLM/De6rVsHM2aYpNhuhzvugCuvPHpwxcGDkJMD6ekwYIA53CLIFBPiL1QxoQRYRESkNnE4YO1aWLPGlD+0aGHamy1ebG7v0AHGjYOWLc33LpfpAmG3m1PeevSAhIRwjV4kKJQAi4iI1BZ5efDVV+agimbNzOrv00+bml6rFa6/Hm6++ehmtoIC2LsX2rUzq74VSbFIlFMCLEEXFxfH0KFDiTvD3cASOxQTEojiogZ5PPDjj2ajW2GhKXn45z/NCW9eL6SlmVXfzp3N/d1us+prtZr63169zMluIaaYEH+higl1gagmdYEQEZGodPgwrFgB33xjOjqUlsKUKaYGGEwLszvvPJrgHj58tC54wADTC1gkSqgLhISNw+FgxowZ2sUrPooJCURxUQO2b4d588wGt+bNzQrwb35jkt/kZHj8cbj/fpP8ejxm1Tc/H/r1g2HDajz5VUyIv1DFhH7HIEHn9XrJy8vTLl7xUUxIIIqLECovN0lvxUa35GSYNAm+/dbc3qePSXwbNTLfFxfDrl0mSR4wANq3P9r9oQYpJsRfqGJCCbCIiEgsyc017c02bzY9ejdsMF0eSkqgTh1T7jBkiElwPR7Ytw/Kykydb9++pkxCJMYpARYREYkFHg/88AN8/bXZ6Na0Kfz1r6brA0CnTjB2rGl7BiYh3rnT3O+ii8ypblZVRkrtoE1w1aRNcNXn8XjIysqiffv2WPVhKigmJDDFRRAVFprE97vvzAruvn3wwgvmuGKr1bQ2u+EGsNlM14fsbHPgRbdupt63YcNwzwBQTEhVpxoT1c3XlABXkxJgERGJOF4vbNtmNrft3Wvamc2eDR99ZG4/6yzT3iwjw3xfVmY2wDVubBLfzp1NUiwSI9QFQsKmvLycSZMmUV5eHu6hSIRQTEggioszVFZmEt+PPzbHFNtspqtDRfJ75ZUwbdrR5Dc31yS/XbvC8OFm9TfCkl/FhPgLVUyoBlhCQi1sxJ9iQgJRXJymnBxT27t1q6nh/ewzePddc4BFo0bw4INw3nnmvg6HSXzr1zeb37p2PXrSWwRSTIi/UMRE5P4NEBERkcrc7qMb3Q4fhsREmDwZNm40tw8YAKNHQ8WvfvPzYf9+6NgR+vc3ybKIKAEWERGJChUb3b791iS427fDa6+Zk93q1oV77oGf/cy0N3M6TYeHunVNh4dzzoH4+HDPQCRiaBNcNWkTXPV5PB7y8/NJSUnRLl4BFBMSmOKimio2ui1bZro7NGwIf/+7SYbBlDSMHQvNmpnvDx409b7t25sV4ebNwzb0U6WYEH+nGhPVzde0AixBZ7FYSE5OxhKGU4QkMikmJBDFRTWUlZnT3NasMRvWCgvhj3+EggJTxztiBFxzjbnN5TKrvvHxMGgQ9OgBCQnhnsEpUUyIv1DFhH68kqBzOBxMnjxZGxnERzEhgSguTiI723R4WLYM6tUz3R2eecYkv61bw5//DNddZ5LfggJz8lvz5jB0qGlxFmXJLygmpKpQxYRWgEVERCKJ220OtFixAoqLzfdPPGHKH8C0MPvlL8FuN7ft2mWun3++6fxQt274xi4SJZQAi4iIRIqCAtPb9/vvISkJVq40B1t4PJCSAg89BOeea+57+DDs2QMtW5oOD23bmg1wInJSSoBFRETCzes1PX2XLTM9fuPj4fnnTVkDwIUXwq9/bZJij8ckvi4X9OljvhITwzt+kSijLhDVpC4Q1ef1enE4HNjtdm1kEEAxIYEpLo4oLYXVq2HdOlPP++238Prr5gCLxEST+A4ebO5bXGxKHpo3N6u+6ekxteqrmBB/pxoT6gIhYeP1eikoKCAlJUUfYAIoJiQwxQWwd69Z9c3KMsnu66+bjg9gevc+9BCkppoV4n37oKQEevaEvn2PHnYRQxQT4i9UMaEuEBJ0TqeTmTNn4nQ6wz0UiRCKCQmkVseFy2VWfOfPh927zWltjz5qkt/4eLjzTnj6aZP8lpbCpk1Qpw5ccYU57CIGk1+o5TEhAYUqJrQCLCIiUpMOHjSHWPzwg0lqP/oIPv/c3NauHYwbB23amFXf7GzT+7d7d9ParFGj8I5dJEYoARYREakJXq/Z1LZ8OeTlmZXd554zp7ZZLHDttXDLLWYFuLwcduyA5GS47DLo0sXUB4tIUCgBlpCw2+3hHoJEGMWEBFJr4qKkBFatgvXrwWo1K8Dvvw9eL57UVAbn5bF0zhyKbrqJxEOH4MAB6NzZbHRLSQn36GtUrYkJqbZQxIS6QFSTukCIiMhp2bMHvvoKtm8337/8MmzbZv588cUU//KXJI0cCUDRxIkkNmliWpt1725Wg0Wk2tQFQsLG4/GQlZVF+/btsVq1z1IUExJYrMdF8aFD5kS31auhtJS4H3/E/t57WJxOvPXrU37PPbj79aO4rOzoY5o3Nz1/mzUzbdBOcvxrYoz1/431mJBTF6qYUHRJ0DmdTmbNmqVdvOKjmJBAYjouDhwgqVEjkgYNotPYsSx/7DES/vEPLE4n/wFaHD5M3eefJ+nGG2l2222+hzW74w6S0tNJSkqq1lesiemYkNMSqpjQCrCIiEiweDymZdnXXwNwEzATaASUAL8BXg7f6ETkCCXAIiIiwVBcDCtXwoYN4PHgHDiQuK++AsCdng4PPMCfW7Tgz263Oc3NZoNzzqG4Y0eatWkDQE5OTsyVNYhEIiXAEnQWi4XU1FSd4iM+igkJJKbiYudOc6Lbrl1QUACvvEJcfr7p+HDjjdhuvJF6cXFw+LA5+KJ1axgwwPT7LS72PU1iYmKtToBjKiYkKEIVE+oCUU3qAiEiIlU4neZEt9WroawMFi82p7sBNG8OY8eadmYej0l83W4491zo3dscfQwUFxf76nmLiopqdQIscqaqm69pE5wEndvtZu3atbjd7nAPRSKEYkICifq4yM+HhQvhyy/NUcYvvHA0+b3sMpg2zSS/xcWmLjg5Ga66ynR5OCbJTUxMxOv14vV6a33yG/UxIUEXqphQAixB53K5mD9/Pi6XK9xDkQihmJBAojYuPB5zjPHcubBxo6n5ffppUwaRnAz/938wejQkJJgewNnZ0LMnDB8O6enm1DcJKGpjQkImVDGhGmAREZHqKioyHR6+/dYcZfzOO/D99+a2vn3h/vuhYUNz286dkJoKmZnQsaOpBxaRiKAEWEREpDq2bzcb3XbvNn9+4w2T6NapA3fdBZdeau6XnW02u519NvTrZxJiEYkoSoAl6CwWC+np6drFKz6KCQkkauKivNxsdFuzxiS28+bB8uXmts6dzUa35s3N/XbsMAnvZZeZ22y2sA492kRNTEiNCVVMqAtENakLhIhILZSba1Z9N2+GnBx4/XU4cMAktr/4BVx3nflzXp653rkz9O8PKSnhHrlIraQuEBI2LpeLxYsXaxOD+CgmJJCIjgu3G777zqz2btwIX3wBzz9vktyWLeFPf4IbbzQb4jZvBpcLLr7YrPwq+T1tER0TEhahigklwBJ0brebJUuWqI2N+CgmJJCIjYvCQvj8c1i0CLZtg7/8xbQ7A7j6apg6FTp0MMlwVha0bw/DhplOD/Hx4R17lIvYmJCwCVVMqAZYREQEwOs1Ce3y5Waj27p1MGeOWQ1u3BgefBB69TKrvVu3mjZnF14IPXqA3R7u0YvIKVACLCIiUlZmTnNbu9as7L73Hvz0k7lt4EC47z5o0MAcc7xvH7RrZ2p9W7YM77hF5LQoAZags1qt9OzZE6t6XsoRigkJJGLiYt8+s9Ft61aT9L7zjkmI69WDe+81fXw9HlMOYbXC+efDeeeZ9mcSVBETExIxQhUT6gJRTeoCISISY1wus9FtxQrT7WH+fFi1ytzWvTuMGQNNm5rWZ3v2mNXeAQOgbdtwjlpETkBdICRsnE4n8+bNw+l0hnsoEiEUExJIWOPi0CH49FP47DNzktvzz5vkNy4O7rgD/vAH081h1y7T4qxPH7PRTclvSOmzQvyFKiaUAEvQeTwe1q1bh8fjCfdQJEIoJiSQsMSF12vals2da+p9Fy6E6dNNbW+bNiYR/vnPzaEWmzebut+rroLBgyExsebGWUvps0L8hSomVAMsIiK1Q0mJWeVdv950eXj7bXNsscUC11wDI0aYNmZ795r79ugBffuaJFhEYooSYBERiX27d5uNbllZpuZ33jyzsS0lxRxlfPbZZuPbpk3m2oUXQqdOZtObiMQcJcASdDabjcGDB2Oz2cI9FIkQigkJpEbiwumEb76BlStNPe/s2abbA5juDvfcY0obcnJMGUS3btCvn+n7KzVOnxXiL1QxoS4Q1aQuECIiUWb/fnOoxcaNsGGDSX4dDkhKglGj4IILzPc7dpgyh379oGtXUPIlErXUBULCxuFw8Pbbb+NwOMI9FIkQigkJJGRx4fHAjz+ajW4rV5pa37ffNslujx5m09sFF0B+PmzfDhkZMHy4KYNQ8htW+qwQf6GKCZVASNB5vV62bt2KfrkgFRQTEkhI4qKoyCS9GzaYld933zV9fO12GDkSrrzSHG28dSvUrWvKIM45x2x+k7DTZ4X4C1VMKAEWEZHYsH27KXnYvNn09/3f/8z19u1h3Dho3RoOHjT1vunp5lCL5s3DOmQRCQ8lwCIiEt3Ky2HdOlizxhxl/N575vAKqxWuuw5uvtn8eft2899Bg6BnT0hICPfIRSRMlABL0MXFxTF06FDi4hReYigmJJCgxEVurmlv9uOPsHQpLFhgDrto1sy0N+vaFQoLTW/fNm2gf3+zEiwRSZ8V4i9UMaEuENWkLhAiIhHE7YYffoCvvza9e//9b7PCC3DJJXDXXVCnjun/63abzW+9e0O9euEctYiEmLpASNg4HA5mzJihXbzio5iQQE47LgoKTI3vwoXmvy+8YJLfBg3gscfgwQfNKvCmTZCcbI4yHjRIyW8U0GeF+AtVTOh3DBJ0Xq+XvLw87eIVH8WEBHLKceH1mpPcli8/2ubs++/Nbb17wwMPQMOGsGcPlJZCr17mKOP69UM2BwkufVaIv1DFhBJgERGJfKWlsHq12ey2bh38619QXGw2st15J1x22dGjjFNTTXuzjh11lLGIBKQEWEREItvevWaj23ffwSefmLpfMAnu2LHQooVpbVZYCN27mxPdGjUK75hFJKJpE1w1aRNc9Xk8HrKysmjfvj1Wrb4IigkJ7KRx4XTCt9/CqlWwfr05yvjAAbOqe9NNcOON4HKZo4yTk02Hhy5ddJpbFNNnhfg71Ziobr6mBLialACLiNSgAwfMqu+338IXX5jNbmBWe8eNM6u/+fnmfh07muQ3NTW8YxaRsFMXCAmb8vJyJk2aRHl5ebiHIhFCMSGBBIwLj8ccYTx3Lnz+OcyYcTT5veIKmDYN2rUzRxk7HHDRRXD55Up+Y4Q+K8RfqGJCNcASEmphI/4UExJIpbgoLoaVK80mt//9Dz76yJQ4NGxoWpv17q2jjGsBfVaIv1DEhBJgEREJvx07TMnDhg3w4YfmSGMwpQ2jR0NSko4yFpGgiYkSiD179nDrrbfSpEkT6tWrR48ePVizZo3vdq/Xy4QJE2jRogV169YlMzOT7yt6R4qISHitXg3z55sV3xdfNMlv3bpm1ffRR8Figc2boWlTGDrUJMVKfkXkDET9JriDBw/Ss2dPLrroIu677z6aNm3K1q1badu2Lenp6QA8++yz/PGPf+SNN96gY8eOPPPMM3z55Zf89NNP1K9mg3Rtgqs+j8dDfn4+KSkp2sUrgGJCAvNkZ5O/eDEpP/6I9eOPoWLhoksX096sadOjRxmfey706aPT3GKcPivE36nGRK3pAjF+/Hi++uor/ve//wW83ev10qJFC8aMGcPvfvc7wBRUN2vWjGeffZZ77723Wq+jBLj6vF4vDocDu92OxWIJ93AkAigmpBK3G374Ae/XX+NasYK4OXOwHDoEcXFwyy3w85+bQy127TI1vgMGQPv2ZiVYYpo+K8TfqcZEdfO1qK8BnjdvHpdddhk33HADS5Ys4ayzzmLUqFHcfffdAGzbto3s7GyGDBnie0xCQgKDBw9m2bJlx02Ay8vLK+04LCwsrHLdarUSHx+P0+nE4/H47muz2YiLi8PhcFQ6ui8uLg6bzVblenx8PFartcoOx/j4eCwWS5Xib7vdjtfrxel0VrqekJCAx+OpdN1isWC323G73bhcrirXXS4Xbrfbdz0YcyovL2fq1KmMHz++yv/XaJ1Txdhj6X2qyTkVFxczdepUxo4dS506dWJiTrH4PtXInAoKYNUq4jdswLtoEfFLlgDgadUK15gxxKWnY9m3D0dpKZx9Npx3HjRogB3w+o09YuYUi+9TmOZ07L8fNpstJuZ07PVYeZ9qck6HDx/2/fuRkJBw0jlVt1tE1CfAWVlZzJw5k3HjxvHYY4+xcuVKHnzwQRISErjtttvIzs4GoFmzZpUe16xZM3bs2HHc5500aRJPPfVUletTpkyhTp06APTs2ZNhw4axYMEC1q1b57vP4MGDyczMZPbs2WzdutV3fejQofTq1YvXXnuNvLw83/URI0bQoUMHpkyZUukNve+++0hOTmby5MmVxjB+/HgKCgqYOXOm75rdbufRRx8lKyuLWbNm+a6npqYyatQoNmzYwPz5833X09PTufXWW1m6dClLjvwDFMw5Hfv/K1bmFIvvU03PaerUqTE3J4i99ynUc2qxeze3z5+PPScHgK/79+fziy/Gdfgw9/30E8kpKUzOyYEVK8xXFMwpFt+ncM2pQizNKRbfp5qa00svvQSYfz+qM6eK+51M1JdA2O12evfuzbJly3zXHnzwQVatWsXy5ctZtmwZAwcOZO/evTQ/pl3O3Xffza5du1i4cGHA5w20AtyqVStyc3N9S+r6yU0rwJpT9eakFeBaPieHA8fy5XjXr8f22WfYPvkEi8eDt0kT/nHllVwzfDgJ+/dDYSHxXbti6dcPR1JSZM8pFt+nCJiTVoA1J/85FRYWntIK8P79+2natGnsl0A0b96crl27VrrWpUsX5syZA0BaWhoA2dnZlRLg3NzcKqvCx0pISCAhwC7jQNfj4+MDPofdbj+l64Fe73jXLRZLwOtWqzXgdZvNhi3A8aBxcXHExVUNA81Jczre9dOZU8X1hIQE3xiifU6x+D6FZE579sCyZdhXrYL33zcHWAAMGoTj7rvZlpVFwvbtJDRoAEOGQNeuEBdHoJlGzJxOcj0q36eTXNecNKdwz8k//zqVOQUS9SvAt9xyC7t27aq0CW7s2LGsWLGCZcuW+TbBjR07lkceeQQwDZWbNm2qTXAhok0M4k8xUQs5nfDNN6aEYfFimDcPysshMRF+/WsYPBhvfj6O/HzsHTtiGTDAdH2QWk2fFeJPm+COY+zYsZx//vlMnDiRG2+8kZUrV/Lqq6/y6quvAuangTFjxjBx4kQyMjLIyMhg4sSJ1KtXj1tuuSXMo49NXq+XgoICUlJS9AEmgGKi1tm/H5YvN/19P/zQHG4BZlPbmDHQqBFs3YrXbqegd29SBgzAcmRvhdRu+qwQf6GKiahvstenTx8++OAD3n33Xbp3784f/vAHpk2bxogRI3z3eeSRRxgzZgyjRo2id+/e7Nmzh0WLFlW7B7CcGqfTycyZM6vUFkntpZioJTwe+OEHmDsX5syBKVNM8hsfD3feCX/4A9jtsGULtGyJ84ormPnllziV6MgR+qwQf6GKiahfAQa4+uqrufrqq497u8ViYcKECUyYMKHmBiUiUpscPmzKHVatgo8/hq++MtfbtoXf/AZatoSdO00v3/PPN+3NlPiKSJjERAIsIiJhtG2bKXn46iv4978hN9cktz//OYwYYWp/N22CVq3MMcZt25rbq9mvU0Qk2M4oAV6wYAGXX3656nSkiuPt5pTaSzERg8rKzPHFq1bBggXw2WemDCI11Rxl3LUr7N1rEt0+fcxXgPZmIsdSTIi/UMTEGXWBsFqttGjRgltvvZXbb7+dLl26BHNsEUVdIEREjpGdDcuWma85c2D7dnP9oovgnnvAaoUdO0xnhwEDICNDJQ8iEnLVzdfOaBNct27d2Lt3L3/605/o3r07/fv35+WXX+bQoUNn8rQS5TweD1u2bKnUAFtqN8VEDHG5YP16+OADePttmD7dJL/168Pvfme6PBw+DLt3wznnwPDh0LFjwORXcSH+FBPiL1QxcUYJ8Lfffsvq1asZPXo0TZo0YeXKlYwePZrmzZtz8803s3DhQqK8zbCcBqfTyaxZs7SLV3wUEzHi4EFYtMis+L70kkmCHQ7o2RNefBF694bNm02yO2QIXHIJNGx43KdTXIg/xYT4C1VMnHEbtF69evHiiy+yd+9e3n//fYYNG4bH42H27NlcddVVtGzZkvHjx/Pjjz8GY7wiIlLTPB746SdzmMW778LUqfD996al2b33woQJ4PWakoeOHc2qb/fuEODUKBGRSBC0PsBxcXFcc801fPDBB+zdu5cXXniBHj16sG/fPl+JRL9+/VQiISISTYqLYckSc4zxjBmm7KG4GDp0gGnTzEpvVpZZCc7MhMsvN5vgREQiWEgOwmjSpAkPPPAAK1euZPLkydhsNrxeL6tWrWL06NG0aNGCO++8k23btoXi5SXMLBYLqamp6g4iPoqJKLVjh1n1/ec/zaEWK1eazW033gjPPWc6OmzZYtqbDR1qevvGx1f76RUX4k8xIf5CFRNn1AXieL7//nvefPNNZs2aRXZ2Nl6vl5SUFG655RZycnKYO3cuZWVl1KtXjwULFjBo0KBgDyHo1AVCRGqN8nJYtw6+/hrmz4f//teUOKSlwbhxpqPDzp2mxKFXL1MDrKOMRSQC1EgXiGPt37+f6dOn07t3b8455xz+/Oc/k5uby+WXX86//vUv9uzZw7Rp03j33XfZvXs3o0ePpqSkhEceeSRYQ5AI4Xa7Wbt2LW63O9xDkQihmIgiOTnmJLfZs+H55+GLL0zyO2QIvPCCOdFt82ZT5nDVVabF2Wkmv4oL8aeYEH+hiokzOgjD5XLx0Ucf8eabb7JgwQKcTider5eMjAxGjhzJyJEjad68eZXHNW7cmOnTp7No0SI2bNhwJkOQCORyuZg/fz7dunXDpk0wgmIiKrjdZmPb8uUmAV6wwLQ8S06G++83h1js3WtqfSsOtUhMPKOXVFyIP8WE+AtVTJxRAtyiRQv279+P1+slMTGRX/ziF/zqV7+qdklD8+bN2bJly5kMQUREztShQ6bcYelSc5TxTz+Z6336mOQ3IcGs+jZtalqbdeigQy1EJKqdUQKcn5/PgAED+NWvfsVNN91Ekt8RlyfzwgsvqCOEiEi4eL0msV22zPT3nTsXSkpMScOdd8Kll0JuLuTnw9lnQ//+ZkVYRCTKnVECvHHjRjp27Hjajz/33HPP5OUlQlksFtLT07WLV3wUExGouBhWrYKvvjItztauNdc7dYKxY6FJE9PhITkZLrsMOncOel9fxYX4U0yIv1DFREi6QMQidYEQkZixc6dZ9f3sM5P8Hjxoktubb4brrzffHzhgDrUYMABSUsI9YhGRaqnxLhAiFVwuF4sXL8blcoV7KBIhFBMRorwcVqwwRxm//DL87W8m2T3rLNPX99prTe/fYw+1CGHyq7gQf4oJ8ReqmFACLEHndrtZsmSJ2tiIj2IiAuTmms4O77wDf/6zOd0N4MorzYluqamwdatpc3Yah1qcDsWF+FNMiL9QxcQZ1QCLiEiEO7a92bx58Mkn5lrjxvDAA9CjB+zaZbo6nH++SXx1qIWIxDglwCIisaqivdn//mcOtqhoO3n++TBqlDnWeNMmc5TxgAHQtm04RysiUmOUAEvQWa1WevbsidWqChsxFBM17Nj2Zv/5jznOuKwM6taFe++FwYNh3z5TE1xxqMWRNpbFxcW+lpZFRUUknuFhFyeiuBB/ignxF6qYUBeIalIXCBGJChXtzSoOtag4bbNrV9PerEED0wUiNdX09c3IMCvBvofXXAIsIhJs6gIhYeN0Opk3bx5OpzPcQ5EIoZioITt3mtXeN9+EP/3JJL9xcXD77fDHP5r77N4N3bvD8OGm528YV9oUF+JPMSH+QhUTSoAl6DweD+vWrcPj8YR7KBIhFBMhVl5uan3//W/4y1/gjTegsBBat4bnnzddHbKyzH0vvdR8NWwYzhEDigupSjEh/kIVE6oBFhGJZjk5psPD55/Dv/5lvgezwvvLX8Lhwyb5rTjUomnT8I5XRCQCKAEWEYlGLpdpb7ZsGXz4oTnVzeMxB1c89BB062ZKIux2s+mtRw+w2ykuLj7h0x57+8nuC6hGWESikhJgCTqbzcbgwYOx2WzhHopECMVEkB08aEoeliyBf/4Ttm831wcPNl0ePB7T8qxNG9PyrGVL30MrNrhVR7NmzU56nzPZR624EH+KCfEXqphQF4hqUhcIEQk7j8e0N6s41OLjj82xxYmJcN99cMEFZpOb221WfPv0Ma3PjmGxWII6JP0TIiKRRF0gJGwcDgdvv/02Docj3EORCKGYCIKioqMrvn/+M3zwgUl+zz0Xpk+H3r1NcpycDFddBYMGVUl+zdMUnfArp6KGGMjJyTnp/c+E4kL8KSbEX6hiQiUQEnRer5etW7dqZUh8FBNnaPt2s+q7YIFJfIuKID7etDe76iqz8e3QIbPq27ev6fV7HKdSs5uYmBjSGl/FhfhTTIi/UMWEEmARkUhVVgbr1sFXX5kWZytWmOvt2sFvfgPNmsHWrdCoEVx2GXTuHNa+viIi0UIJsIhIJNq3z6z6fvaZaW+Wnw8WC1x7LdxyCxQUmJXhzp1Ne7MmTcI9YhGRqKEEWIIuLi6OoUOHEhen8BJDMXEKXC749lvT3uyDD+CLL8DrNf17x441/Xx37jT1vT/7GZx9timHiEKKC/GnmBB/oYoJdYGoJnWBEJGQ27/frPouXgyzZ8OuXeb6xRfD3XebTW/79kF6uln1bdEi6EMoLi72tUorKipSn18RiSrqAiFh43A4mDFjhnbxio9i4iQ8HvjxR7Pi+9e/wosvmuS3fn0YPx7uvx9yc03Zw8CBZuNbCJJfMBvfvF4vXq835Mmv4kL8KSbEX6hiQr9jkKDzer3k5eVpF6/4KCZOoLAQVq6E//3PrPpu3Giun3cePPCAOclt0yZzmMWAAWYDXIxQXIg/xYT4C1VMKAEWEQkHrxe2bTO1vgsWwNy5UFxsEt5f/Qouv9yUOxw8aHr89u0Lp3CKm4iIHJ8SYBGRmlZaCqtXm/Zmc+aYPwNkZJiNbikp5lCLlBTIzDQb39TeTEQkaLQJrpq0Ca76PB4PWVlZtG/fHqv+0RYUE5Xs3m1WfT//3PT2PXDAJLc33mi+Dh4017p2hf79oXHjcI84ZBQX4k8xIf5ONSaqm68pAa4mJcAickYcDtiwwXR5+OAD0+kBoHlzs+qbng47dkBioil36N4d1ApKROSUqAuEhE15eTmTJk2ivLw83EORCFHrYyI3FxYuhHfegeeeO5r8XnYZTJt29ES3Nm1g2DBzpHEtSH5rfVxIFYoJ8ReqmIj9T1gJC7WwEX+1MibcbvjhB7PqO3++SYJdLkhOhgcfhF69zKEWNhsMGgQ9e0JCQrhHXaNqZVzICSkmxF8oYkIJsIhIKBw6BF9/DUuXmqOMN20y1/v2NX194+LMtVat4PzzzeqviIjUCCXAIiLB5PGYDg7Ll5sV33nzTNeHOnXgrrvgkktg714oL4c+fUxCrNPWRERqlDbBVZM2wVWfx+MhPz+flJQU7eIVoBbFRFERrFplujy8/z6sXWuud+5sNro1amQ2uqWmmg4PGRm1ur1ZrYkLqTbFhPg71Ziobr6mFWAJOovFQnJyMhaLJdxDkQhRK2Ji+3az6vvZZyb5PXTI1Pb+4hdw3XWwf7+p9+3WzSS/jRqFe8RhVyviQk6JYkL8hSom9OOVBJ3D4WDy5MnayCA+MR0TZWXmQIv334eXX4a//90kvy1bwp/+BD//OWRlmQ1xl1wCl16q5PeImI4LOS2KCfEXqpjQCrCIyOnau9eUO3z5pTnUYt8+c/3qq+H226GkBLZsgQ4dzEa3tLTwjldERAAlwCIip87phG++MV0ePv4YFi0yK7yNG5v2Zueee7S92YUX1sr2ZiIikUwJsIhINRUXF5OUlGT+PG4c9ebONQdYAAwcCPfdZza1qb2ZiEhEUxeIalIXiOrzer04HA7sdrs2MggQIzHhdlO8di1JfftyF/BqQgKW8nKoVw/uvRcGDzYlEGVlZgVY7c1OKibiQoJKMSH+TjUmdBSyhI3X66WgoAD9bCUVoj4mDh0y3R0+/JAPgb+CSX7PPhumT4cBA0zv37p14YorTDKs5Pekoj4uJOgUE+IvVDGhBFiCzul0MnPmTJxOZ7iHIhEi0mOiuLg48FdREcXr11M8ezZlb79N3enTGQ6UA4dvvpnixx+nuLyc4s2bKW7fnuKLL6a4ZUuKS0srPY8EFulxITVPMSH+QhUTqgEWkVqvoq43kERgCnDPke+/BUYA3773Hrz3XrWeX6tZIiKRRSvAIiLH0R9Yj0l+PcCfgD6YJFhERKKXVoAlJOx2e7iHIBEmkmOiqKjI/KGsDNavhzVriF+0iPgvvsDi9eJp0oTy++/n9rZteeSOOwDI+egjEgcMUHuzMxTJcSHhoZgQf6GICXWBqCZ1gRCJcXv3mqOMv/4a5syBbdvM9cxMuOcesFgo3rqVpN//HjBJc6I2uomIRBR1gZCw8Xg8bNmyBY/HE+6hSISI6JhwOmHNGpg7F959F1580SS/SUnwyCMwdiwUFpoWZz16hHu0MSWi40LCQjEh/kIVE0qAJeicTiezZs3SLl7xidiYyMuDhQtN8jtjhln5dThMojt9OvTpYw61iI837c0uvDDcI44pERsXEjaKCfEXqphQDbCI1D5uN/z4oyl3WLrUJL5FRWC3w8iRcOWVcOAAbN8OnTubE90aNyYRdXQQEYkFSoBFpHY5dAhWrIBVq+A//zF1vwDt28O4cdCihSmBqFsXfvYzc9hFfHxYhywiIsGlBFiCzmKxkJqaqmMsxSciYsLjMae1ff21SYD//W9TAmG1wnXXwc03Q2mpuU/79mbVt0WL8I23FoiIuJCIopgQf6GKCXWBqCZ1gRCJYkVFZsV3zRpYtMgca+z1QrNmZpNb586we7cpjejVC3r3hjp1wj1qERE5ReoCIWHjdrtZu3Ytbrc73EORCBG2mPB6TTnDvHnw0Ufwl7/Ap5+a65dcAi+8AG3bmlXf5GS4+moYOFDJbw3RZ4X4U0yIv1DFhBJgCTqXy8X8+fNxuVzhHopEiLDERGkpLFtmkt9//9t0ddixAxo0gMcegwceMO3N9uyBc86BYcNM6YN+9Vpj9Fkh/hQT4i9UMaEaYBGJPbt3m81t69aZFmc//GCu9+5tEt/ERNiyxaz6XnYZdOliaoFFRKRWUAIsIrHD4fAdZczy5fD++1BcbI4rvvNOk+weOGDKIjp2hAEDIDU13KMWEZEapgRYgs5isZCenq5dvOJTIzGRk2OS3g0bYMEC0+kBTKI7diykpZm+vvHx5kCLHj1M318JG31WiD/FhPgLVUyoC0Q1qQuESIRyueC77452efjXv8wqr9UKN90EN94IJSWm1rdNG7Pq26pVuEctIiIhoC4QEjYul4vFixdrE4P4hCwmDhwwbc0WLIDZs+GVV8y1Fi3guedMArxvn+n327ev6fKg5Ddi6LNC/CkmxF+oYkIJsASd2+1myZIlamMjPoFiori4GIvFgsViobi4+NSe0OMxRxl/+KFpazZjhuntC3DFFTBtGrRubdqb1atnrl14ofmzRAx9Vog/xYT4C1VMqAZYRKJLYaGp7/3mG/jqK5g/35RBNGwIDz5oOj3k5cHBg9C9O/Tvb24TERE5QgmwiEQHrxe2bjVHGX//vVn93bjR3Na/P4webVZ4t2wxbc4uuQS6dYM4fcyJiEhl+pdBgs5qtdKzZ0+s6qtaKwUqZ3A6nXTp0oXS0lJfHdex9ztpCURxMaxdS+LGjWaj25w55qCLunXh7rvh4ouhoMAkyOnpcP75puuDRDR9Vog/xYT4C1VMqAtENakLhEj1hKp9URMgv08f0+0BzOEVY8eaPr67d5tr550HPXvqKGMRkVpKXSAkbJxOJ/PmzcPpdIZ7KBIjLge+BZP8xsXBbbfBxInmWONNm6BxY7jqKtPiTMlv1NBnhfhTTIi/UMWEEmAJOo/Hw7p16/B4POEeioRBUVFRla/9+/fz2GOPsX//ft+1nJwc32NycnIqP2bzZorefpui//s/HJmZLACag2lh9qc/wXXXQW4u7N0LvXrB0KHQrl24piynSZ8V4k8xIf5CFROqARaRoEpMTKxyLS4uDrvdTmJiIgkJCQEfk5iYCE6n6e6wapVpYfbeeybJBZPk3nab+fPmzdCokTnauHNnc+iFiIhINSkBFpHIkJdnjjL+8UfT6WHePNPvt3FjeOghU9u7fz/k50OnTqbcISUl3KMWEZEopARYgs5mszF48GBsNlu4hyIRIlBMJCYm4vV6we02bc1WroSsLNPhYdMmc6dBg+DXvzbtzbKywG6HzEw491yIjw/PZCRo9Fkh/hQT4i9UMaEuENWkLhAiIXDggDnU4vvvTenD7NlQXm76+N57LwweDEVFsGcPtGlj2pu1bBnuUYuISIRSFwgJG4fDwdtvv43D4Qj3UCRCVImJiqOM580z5Q7vvgv/+IdJfs8+G1580RxdvGePKY3o2xeuvlrJb4zRZ4X4U0yIv1DFhEogJOi8Xi9bt25Fv1yQCpViorDQlDt8+63ZzPbOO+YQi4r2ZsOGgcNhbktNNYdcZGRAiPoLS/jos0L8KSbEX6hiIuZWgCdNmoTFYmHMmDG+a16vlwkTJtCiRQvq1q1LZmYm33//ffgGKVJbZWWZVd/ly2H+fJg50yS/bdvClClwzTVmo9uOHeYY4+HDoWNHJb8iIhJUMZUAr1q1ildffZVzzjmn0vXnnnuOKVOm8Je//IVVq1aRlpbGpZdeyuHDh8M0UpFapuKo408+gfXr4aWX4L//NYntz38Ozz8PZ50FW7aAywWXXAKXXgoNG4Zz1CIiEqNiJgEuKipixIgR/PWvf6VRo0a+616vl2nTpvH4449z7bXX0r17d958801KSkp45513wjji2BUXF8fQoUOJi1OFjQA7dhD3yScMtVqJX7wYpk6FnBxT3vDHP8Idd0BJCWzdaja6DRsGPXqYkgiJafqsEH+KCfEXqpiImS4Qt99+O40bN2bq1KlkZmbSo0cPpk2bRlZWFunp6axdu5aePXv67j98+HAaNmzIm2++GfD5ysvLKS8v931fWFhIq1atyM3N9e0qtFqtxMfH43Q6K51QYrPZiIuLw+FwVKpZiYuLw2azVbkeHx+P1Wqt9HoV1y0WS5XCb7vdjtfrrXIsYEJCAh6Pp9J1i8WC3W7H7XbjcrmqXHe5XLjdbt91zUlzCtqcysuxf/cd7rVrce/ZQ9zs2VizssyNF12E6667cNepYza6eb1Ye/Qgvm9fnDZb5M4pFt8nzUlz0pw0pxia0/79+2natOlJu0DExI9Y7733HmvXrmXVqlVVbsvOzgagWbNmla43a9aMHTt2HPc5J02axFNPPVXl+pQpU6hTpw4APXv2ZNiwYSxYsIB169b57jN48GAyMzOZPXs2W7du9V0fOnQovXr14rXXXiMvL893fcSIEXTo0IEpU6ZUekPvu+8+kpOTmTx5cqUxjB8/noKCAmbOnOm7ZrfbefTRR8nKymLWrFm+66mpqYwaNYoNGzYwf/583/X09HRuvfVWli5dypIlS3zXgzWn5ORkRo0aFVNzisX3KaRzSkhglMXC3jVrSPvgA6wuFyV167L6hhu48PrrWbpjB0v27z86p0OHGFanDgvmzYvcOcXi+xTGObVu3Zpnn3220j960T6nWHyfanpO8fHxPPzwwyxbtixm5hSL71Okzmnq1KlUR9SvAO/atYvevXuzaNEizj33XIBKK8DLli1j4MCB7N27l+bNm/sed/fdd7Nr1y4WLlwY8Hm1Anz6cyovL2fq1KmMHz++yv/XaJ1Txdhj6X0KyZyKi+G772DtWix5edg//BA2bABgS3o6zX73OxJSUrDn5+MqLMTdpQv06QMNGkTunGLxfYqQOTmdTiZPnszYsWN9R2RH+5xi8X2qyTkd+++HzWaLiTkdez1W3qeanFNhYSFTp071fU5oBfiINWvWkJuby3nnnee75na7+fLLL/nLX/7CTz/9BJiV4GMT4Nzc3CqrwsdKSEjwfSCf7Hr8cU6kstvtp3Q90Osd77rFYgl43Wq1Brxus9kCnqISFxcXsK5Gc9Kcjnf9uHPKzydh+XLTvmzbNtPX9/BhsNtx3n47s1q0YHxyMvZt26BBA+Iuv5y4Ll3AWnkrQkTNKRbfpwibEwT+XI3mOcXi+6Q5aU7hnpP/58SpzCmQqE+AL774Yr799ttK1+644w46d+7M7373O9q3b09aWhqffvqprwbY4XCwZMkSnn322XAMWSS2uFxHjzLOzYVFi+DLL81tHTrAuHF4mjc3K8E7dkCnTjBggNkEJyIiEgZRnwDXr1+f7t27V7qWmJhIkyZNfNfHjBnDxIkTycjIICMjg4kTJ1KvXj1uueWWcAw55sXHxzNixIjj/vQnMeTAAdPT98cfTfL7xhvm5DarFa6/Hm66CaxW4nfuZERyMvH9+sF558EJVgOl9tBnhfhTTIi/UMVE1CfA1fHII49QWlrKqFGjOHjwIP369WPRokXUr18/3EOLSVarlQ4dOoR7GBJKHg9s3GiOMc7NPXqwhdcLaWkwdix06QJFRbBrF9ZWregwcCC0bh3ukUsE0WeF+FNMiL9QxUTM9AE+1uLFi5k2bZrve4vFwoQJE9i3bx9lZWUsWbKkyqqxBE95eTmTJk2qUowvMaKgAD77DBYuNCUNr7xiTnfzemHIEJg2DTp3hr17Tb/fPn0ov+wyJs2apZiQSvRZIf4UE+IvVDFRK1aApeb5786UGOD1mg1uX38N2dmm28O775oa4ORkGD0a+veH8nLYtAkaN4bBg81Rxk6nYkICUlyIP8WE+AtFTCgBFpGTKyqCVavMRrbiYvjnP6Fi82mfPnD//dCoEezfD/n5ZgX4/PNNEiwiIhJhlACLyPF5vbB9u6nx3bULdu6EN980SXBCAtx5J1x2Gbjd5ijjOnXgoovgnHNAm1hERCRCRf1BGDWlsLCQ5OTkkzZWFvB4POTn55OSkoLVGpNl5rVDaSmsWQPr1kFJCXz0ESxdam7r2BHGjYMWLUyv3927oW1bGDgQzjqrylMpJiQQxYX4U0yIv1ONierma1oBlqCzWCwkJydjsVjCPRQJoLi4mKSkJACKiopITEyseqddu8yq744dptXZX/9qyhusVrj5ZrjhBrBYTOLrdEK/fqYUol69gK+pmJBAFBfiTzEh/kIVE/rxSoLO4XAwefJkbWSIRuXlsGKFaWm2cycsWQKTJ5vk96yz4E9/Mgmw02k2xNWrB1dcARdeeNzkFxQTEpjiQvwpJsRfqGJCK8AiYuzbZ1Z9t26FsjKz6rtrl7ntyivhjjtM3W9eHhw8CN27m64PDRuGddgiIiKnSgmwSIwqLi4+6fXi4mKzmvvdd7B2LRQVEb9hA/Fz5mBxu01nhwcfNKe3OZ2wZYtZ6b34YpMABzj/XUREJNLpXy+RGFVR53sizZo18/25LfAP4IIj388Brps+HRo0MIdf7NsH7dub9mbNm4dgxCIiIjVDXSCqSV0gqs/r9eJwOLDb7drIEEan8v9+JPAiUB8oBB4A3gK8H35oNrq53dCrF/TubVqdnSLFhASiuBB/ignxd6oxoS4QEjZer5eCggJSUlL0ARZGRUVFAa8X795Ns86dAch96ikaffQRcatWAeDu0oW40aOZ0bQpM0pKzIluzZubVd927Uznh9OgmJBAFBfiTzEh/kIVE+oCIUHndDqZOXMmTqcz3EOp1RITEyt/1a1L4s6dJC5eDMAVQJMpU0zyGxcHt9+ObeJE6rVuTeKhQyTm58O558KwYab04Qw+eBQTEojiQvwpJsRfqGJCK8AitUFhoWlv9t134HIxA7gPTG1v69bmUIv27c1Gt6wsU/c7ZAh07Qo2W5gHLyIiElxKgEVimddrOjesWGE2sTkc1H3pJZP8As6rriL+jjvAbjetzXJyoEMHc6Jb06ZhHbqIiEioKAGWkLDb7eEeghQVwerVsGGDWcVduxZmz8bq8bALs/Ft3u23E2+zwfbt5pS3QYOgZ0/T7zfIFBMSiOJC/CkmxF8oYkJdIKpJXSAkani9JqFdvtx0cIiLg1deMRvawJza9utfQ1ISlJSYE9/OOstsdGvbNpwjFxEROSPqAiFh4/F4yMrKon379lit2mdZo0pLYc0aWLfOfJ+VBW+8YY44Tkw0ie/gwSZJ3rcPiovNim+/flC/fsiGpZiQQBQX4k8xIf5CFROKLgk6p9PJrFmztIu3pu3aBfPnm5VfqxXeesus/JaXwznnwIsvmuTX4TB1wTYbXHYZ/OxnIU1+QTEhgSkuxJ9iQvyFKia0AiwS7crLYf16s/LrcsH+/fDMM6bzQ3w83HYbDB1qkuIDByA3Fzp1ggEDIDU13KMXERGpcUqARaLZvn1mxXfrVtO67P334fPPzW3t2pn2Zm3amJPctm83q74XXmjKHrTRREREaiklwBJ0FouF1NRUneITSk4nfPON6fJQWmrKGp580rQxs1jg2mvhllvMCnBRkSmPaNXKtDdr3brGh6uYkEAUF+JPMSH+QhUT6gJRTeoCIREjNxe+/tp0dUhOhkWLzMqvx2N6944dC926Hd3oVlJiTnTr189shBMREYlR1c3XtAlOgs7tdrN27Vrcbne4hxJbXC6z6jt3rtnEFhcHkyfDv/9tkt+f/cxsdOvWzdQFb95sVoCvuAIyM8Oa/ComJBDFhfhTTIi/UMWEEmAJOpfLxfz583G5XOEeSuw4cMCs9H76qVnZ3bgRfvc70+asfn0YPx7GjIF69cwmuO3bzUa34cOhc2ezAS6MFBMSiOJC/CkmxF+oYkI1wCInUFxcTFJSEgBFRUUk1vQqqsdjkt2vvzZJcP36MHOm6foA0KsXPPggNG5sNrrt2GE2t2VmmrKH+PiaHa+IiEgUUAIsEqkKCmDFCvjuO1O+kJMDTz1lNrXZ7fCrX5nyBovl6Ea3tm3NiW4tW4Z79CIiIhFLCbAEncViIT09Xbt4T5fHY2p8ly+HvDyzuvvWW7Bkibm9QwfT3qxlS1MOsWcPlJVB377mq1698I4/AMWEBKK4EH+KCfEXqphQF4hqUheI2qnGSyAOH4ZVq8xmt4QEU8/7wguQn2/qeG+4AW66yWyAKy83tb4pKeZQi44dzWqwiIhILaUuEBI2LpeLxYsXaxPDqfB6zWEW8+aZE90aNTKb3p54wiS/aWmm48OIESb5zc83yW+XLjBsmNnwFsHJr2JCAlFciD/FhPgLVUwoAZagc7vdLFmyRG1sqqukBL78Ej7+GA4dMhvXnnwSPvzQJMaXXWZWgTt3Nq3Qtm41q78XXQRDhpgSiQinmJBAFBfiTzEh/kIVE6oBFgmnHTtMre+uXWaV94sv4O23TaKbnAwPPGDqesGUR+zebTa6DRwIZ50V1qGLiIhEKyXAIuFQVgbr1sHataZ9WXIyPPec6fgAJum9/35o2NBsitu71xx33K8f9OkTkRvdREREooUSYAk6q9VKz549sYb58IWItWePWfXdtg2aNTM9fV95xZRC1KkDd90Fl15qanrLyswqcWoqXHKJ6QARwbW+x6OYkEAUF+JPMSH+QhUT6gJRTeoCUTsFtQuEw2GS3TVrTA1v48bw6qvw1Vfm9s6dYexYaN7cfJ+XBwcPQteu0L+/2RgnIiIix6UuEBI2TqeTefPm4XQ6wz2UyJGdDf/5j9nsVrcuFBaaZPerr8Bmg1tvhUmTTPJbsdHN6YSf/cysBkd58quYkEAUF+JPMSH+QhUTSoAl6DweD+vWrcPj8YR7KOHndJpa37lzISvLbFz7979Nl4cDB8z3zz0HN95oEuHCQti82RxyMXQo9Oxp2p5FOcWEBKK4EH+KCfEXqpiI/n9ZRSJVXp45ynjjRlPuYLHAI4+YTg4AV10FI0eaAy88HlMb7HKZQy169zYrxSIiIhJ0SoAlKhUXF9f461T7Nd1u+PFHWL0aCgrgrLOIX7iQ+H/9C4vbjadRIxz33Ye7Rw/T5/fQIdMGLTXVtDdLTwePhxCfOSciIlJrKQGWoLPZbAwePBibzRay16jYmFaTmjVrdsqPaQ+8BQw88v2/gF8fPMiBiRNP+thY2p9aEzEh0UdxIf4UE+IvVDGhLhDVpC4QkcUSBa3AfgVMA+oDBcD9wNun8Hj91RQRETk11c3XtAIsQedwOJg9ezY33ngjdrs9JK9RVFQUkuf1V1xc7Fv5zcnJCdwGraDAlDv88IM5oKJuXRJefZW41asBcHftSvzo0bycmsrLYDa67dsH7dqZ9mZpaTUyl3CqiZiQ6KO4EH+KCfEXqphQAixB5/V62bp1a0hXMM+oH+8ZvGal1/V6YcsWc6hFTo5JaL//Hl580STFcXHwy19iGzaMejab2ei2e7epEc7MNBvd6tSp8XmEQ03EhEQfxYX4U0yIv1DFhBJgkdNRVASrVsGGDaaLQ6tW8Prr8Mkn5vY2bWDcOJMUA5SWws6dZrW3f3+z0S0KyjhERERikRJgkVPh9ZojjJcvh717Tb/e3bvhD38wZQ0WCwwfbg62qPhVTW6u6fRw9tkm+U1ODusUREREajslwBJ0cXFxDB06lLgYOMChkpISc4zx+vXm0Ip27WDOHJg925Q3pKTAmDFwzjnm/k4n7NgBiYnmNLdu3czjaqGYjQk5I4oL8aeYEH+higl1gagmdYGonYqLi30t14reeIPEnBxzelthIUydak5tAxg8GO69FyrasxUUmBXh9HQ4//xasdFNREQk3Kqbr+koZAk6h8PBjBkzcDgc4R7KmSsrO/rnvDzIyIClS81K7+bNJuH97W/hN78xf/Z4TK3vgQMm8b3ySiW/xFhMSNAoLsSfYkL8hSom9DsGCTqv10teXl707+Ldswf++1/ft5YGDeCPfzRlEADnnmsS4SZNzPelpabkoXlzc5xx+/ba6HZEzMSEBJXiQvwpJsRfqGJCCbCIP4fD1PmuWWNKHYCfA3UffhgOHzab226/Ha66CqxHfomSk2PKHs45xyS/KpMRERGJWEqARY6VnQ1ff236+6akQL16vA6MBJP8tm9v2pu1bm3uX7HRLSmp1m90ExERiRbaBFdN2gRXfR6Ph6ysLNq3b4/VGiVl5k4nfPut6e1bUmL6+G7ahGfKFKx5ebgB9zXXYP/lLyE+3jymoMC0QuvQAQYOhCMnxklVURkTEnKKC/GnmBB/pxoT1c3XlABXkxLgGJaXZ1Z9f/oJGjeGhg3hnXfg/ffB6yULuA34ZPZsEuvUOXqim8cDvXrBeefVmhPdREREIpm6QEjYlJeXM2nSJMrLy8M9lBNzucyq79y5sGkTtG1rNrI9/LDp7+v14rzoInoAX1U8pqTE3Dc52dQAn3++kt9qiJqYkBqluBB/ignxF6qYUA2whETEt7A5cABWrIAffjDJbHo6fPQRvPWWKYdo0ABGj8bRsyeHKzpB5OaatmjnnmtOdNNvAk5JxMeEhIXiQvwpJsRfKGJCCbDULh4PbNxoSh4OHDC1vkVFMGECbNhg7tO7NzzwADRqVLkPsNcLQ4ZA167a6CYiIhLFlABL7VFQYFZ9v/vOHE/csSN8+SW8/DIUF0NCAvzqV3D55b7+vYnl5Xife86sEGujm4iISEzQJrhq0ia46vN4POTn55OSkhIZu3g9HtPWbPlys+GtVStz7eWXTQIMJhkeO9Ycc1zxmIqNbuedZ74SEsI3hygXcTEhEUFxIf4UE+LvVGOiuvmaVoAl6CwWC8nJyVgi4RS0w4dNa7NvvjEJbEaG2fj2wguQn28OsrjpJrjhBog78tehpMQcZ9y8udnk1q6dTnQ7QxEVExIxFBfiTzEh/kIVE/rxSoLO4XAwefLk8G5k8Hph61bT4WHNGlO60LQpvP46/N//meS3RQt49ln4xS+OJr/Z2eYI5HPPhWHDdJxxkERETEjEUVyIP8WE+AtVTGgFWGJPcTGsXm02tdlsZtV3xw54/nnYtcvc5/LLTb1vRQszpxO2b4f69c1Gty5dtNFNREQkRikBltjh9ZpEd/lyk+i2bAn16sGHH8KsWabvb8OGpsNDnz5HH3fokFn51UY3ERGRWkEJsMSG0lJYuxbWrTOJcMeOpszhj380vX7B9O4dPdr0/QWzwa1iRfiCC8ypbtroJiIiEvPUBaKa1AWi+rxeLw6HA7vdXjMbGXbtMqu+27ebut4GDeCLL+DVV01iXLcu3H03XHzx0XrekhLzuGM3uknI1HhMSFRQXIg/xYT4O9WYUBcICRuv10tBQQEpKSmh/QArL4f1680mN6fT1PqWlMDkySYhBlPLO3YspKUdfVx2tukOcc45OtGthtRYTEhUUVyIP8WE+AtVTKgLhASd0+lk5syZOJ3O0L3Ivn3w8cemj29SkunWsGGDqe9dvtxsYPvlL2HixKPJr9MJmzebVeAhQ8yKsJLfGlEjMSFRR3Eh/hQT4i9UMaEVYIkuTqfp6btqlTmmOD396KEW//mPuU+rVjBunLmtgja6iYiIyBFKgCV65ObC11/Dpk2QkmJObdu8GaZMMb17AYYOhdtuO7qZTRvdRERExI8SYAkJu90evCdzueC778yq7+HDZsOazQbvvWe+PB5o3Bgeegh69jz6uIoT3Vq00Ea3CBDUmJCYobgQf4oJ8ReKmFAXiGpSF4gwyc83q74//QSNGkFqKuzdC1OnmmsAgwbBr39tDrGoULHRrXt3bXQTERGpJaqbr2kTnASdx+Nhy5YteDye038St9us+s6daxLdNm1M2cPChWal96efIDERfvMbePjho8mvNrpFpKDEhMQcxYX4U0yIv1DFhBJgCTqn08msWbNOf8fmwYPw6afwySemvKGivdkzz8CMGab92dlnw4svwuDBR3v7HjoEW7eaUodhw8zqr44zjghnHBMSkxQX4k8xIf5CFROqAZbI4fGYld0VK0zpQ+vW5hCLFSvgL3+BggKIizPtzYYPB6v16OO00U1ERESqSQmwRIbCQlPr+/33UK+eOcq4tBSmTzerwQBt25r2Zm3bHn2cNrqJiIjIKVICLEFnsVhITU2t3oktXi9s2WIOr8jJMau+9erBjz+ajW7Z2abE4Zpr4NZbIT7+6GMrNrqde642ukW4U4oJqTUUF+JPMSH+QhUT6gJRTeoCEQJFRaa12YYNpmShRQtTzvDee/Dvf5s/p6aao4y7dz/6OKcTtm83CW///tC169FyCBEREam11AVCwsbtdrN27VrcbnfgO3i9kJVlOjysXm1OZWvZ0hxm8cgjMHu2SX4vushsdDs2+Q200U3Jb8Q7aUxIraS4EH+KCfEXqphQ5iBB53K5mD9/Pi6Xq+qNJSXwv//BRx+ZZDYjw7Qz++gjs9K7ZYtpafbII+b7xETzOI8HduwwG+EuuACuvBKaNq3RecnpO2FMSK2luBB/ignxF6qYUA2w1JwdO0yt765dptyhQQPYv9+s8q5bZ+7Tsyc8+CA0aXL0cSUl5jFpaTBwoDa6iYiIyBlRAiyhV1oKa9fC+vXmgIuMDNOf96uvTF/fw4fBboeRI+Gqq4729QXIzqY4P5+kxx4DoKioiMSwTEJERERihRJgCTqLxUJ6errZsbl7t1n13b7drOA2bAjFxfDqq/Df/5oHdOhg2pu1bHn0SY7d6HbxxWGYhQRTpZgQOUJxIf4UE+IvVDGhLhDVpC4Qp6i83Kz4rl0LDodpbxYXZ443njoV8vLM5rXrr4ebbqrc3uzQIdPiLD0dBg6kOCmJpKQk4MgKcKLWgEVERKQqdYGQsHHt3s3iv/8d15Ilpqdv+/am88Prr8Pjj5vkNy0NJk2q3NvXf6PbVVeZDhES9VwuF4sXL9bGFqlEcSH+FBPiL1QxEfUJ8KRJk+jTpw/169enadOmXHPNNfz000+V7uP1epkwYQItWrSgbt26ZGZm8v3334dpxDHM6YQ1a3DPn8+S3Fzc7dpB48amlOE3v4EPPjCJ8JAhMG0adOly9LElJbB5symRuOoqGDBAxxnHELfbzZIlS9TaSCpRXIg/xYT4C1VMRH0N8JIlSxg9ejR9+vTB5XLx+OOPM2TIEH744Qffr8qfe+45pkyZwhtvvEHHjh155plnuPTSS/npp5+oX79+mGcQI3JzzVHGmzZBo0bmms0GH34Ib70FLhckJ8Po0ebwimNVnOh2zjkUn322qfstLvbdXHycPx9LZREiIiJSXVGfAC9cuLDS96+//jpNmzZlzZo1XHjhhXi9XqZNm8bjjz/OtddeC8Cbb75Js2bNeOedd7j33nvDMezY4XKZut5Vq0wS264d2Gw0+O474p980twG0KcP3H//0eQYjm50q1/frAp36UJS3IlDstlxSiJUyi4iIiLVFfUJsL+CggIAGjduDMC2bdvIzs5myJAhvvskJCQwePBgli1bdtwEuLy8nPLyct/3hYWFVa5brVbi4+NxOp14PB7ffW02G3FxcTgcjkqJWVxcHDabrcr1+Ph4rFZrpderuG6xWHA4HJWu2+12vF4vTqez0vWEhAQ8Hk+l6xaLBbvdjtvtrlQ/U3Hd5XJV+rXCKc3pwAHiVq3CtnkzjgYN8B6p9fUuXsz9r76KtbQUb0ICrl/9Cs+ll4LFQrzXiwVwHDwIOTkmYe7fH3vLlmeUxJaXlwdnTsTg+xQBc3I6nZxzzjk4nc6YmVMsvk/hmFNFXMTSnGLxfaqpOVV8Vlit1piZ07HXNadTn9Ox/35UZ07+/w+OJ6YSYK/Xy7hx47jgggvofuT43OzsbKDqymGzZs3YsWPHcZ9r0qRJPPXUU1WuT5kyhTp16gDQs2dPhg0bxoIFC1hXcZADMHjwYDIzM5k9ezZbt271XR86dCi9evXitddeIy8vz3d9xIgRdOjQgSlTplR6Q++77z6Sk5OZPHlypTGMHz+egoICZs6c6btmt9t59NFHycrKYtasWb7rqampjBo1ig0bNjB//nzf9fT0dG699VaWLl3KkiVLfNdPeU4JCfTq0IHXtm7l8JYtXPXxx3SvqK/u2JFXrriCnEaN4JtvzJw6dSI5L4/JBw6Y+/zwA/zwg29Ojx3p91sxp4cffphvv/2WAQMGAPDwww/TokUL7rnnHtavX89//vMfACZPnhy8OcXi+xQhc/rmm29ibk4Qe+9TTc5p48aNfHPk8yFW5hSL71NNzyk+Pp7FixfH1Jxi8X2qiTlNnz4dh8Ph+5w42ZymTp1KdcRUG7TRo0fz8ccfs3TpUloe6Sm7bNkyBg4cyN69e2nevLnvvnfffTe7du2qUkJRIdAKcKtWrcjNzfW11aiVP7nl5uJduRI2boT69Ylr1gybzYZz7Vripk/Hsn8/XquVjZdfToc778RjPWafZUkJ8bt3Y2neHEefPtC2re/QixPN6fDhw77/5/v37ycpKUk/YUfZnEpLS/n000+59NJLsdvtMTGnWHyfanpObrebjz76iEsvvZT4I91gon1Osfg+1fQK8KeffsrVV1+NxWKJiTkdez1W3qeanFNRUZHv34/4+PiTzmn//v00bdr0pG3QYmYF+IEHHmDevHl8+eWXvuQXIC0tDTArwccmwLm5ucetJwUTBAkBuhAEul7xwe3Pbref0vVAr3e86xaLJeB1q9Ua8LrNZsNms1W5HhcXR1yAutsqc/J44Mcfsa9YAfn5pq9v3bqm3+/f/058xU+FZ52F86GHmF1SwnirlYSK18zJMe3NevSA/v1JCBCUJ5pThYSEBN//vzOe0xEx9T4dEWlzio+P55tvvuHKK6/0jSHa5xSL71NNz8npdPriwn9M0ToniL33CWp2TsfGRKzMqUIsvU8VQj2nY//9OHZcpzKnQKK+DZrX6+X+++/n/fff54svvqBdu3aVbm/Xrh1paWl8+umnvmsOh4MlS5Zw/vnn1/Rwo1NBAXz2GSxcaBLejh1N8rt1qznBrSL5vfJKmDYNb0bG0cc6naa9WUX7s4svNl0eRERERMIk6leAR48ezTvvvMPcuXOpX7++r+Y3OTmZunXrYrFYGDNmDBMnTiQjI4OMjAwmTpxIvXr1uOWWW8I8+gjn8cCWLaa9WU6OWfWtVw/cbnj/fXj3XdMFolEjePBBOO8887iKX4UUFpoWZx06wMCBOtRCREREIkLUJ8AVBdyZmZmVrr/++uuMHDkSgEceeYTS0lJGjRrFwYMH6devH4sWLVIP4BMpKjKtzTZsMAdSdOxoji7OzjZHGf/4o7nfgAGmt+8xq7o2r5fBdetiO3TIJL69e+tQi1rOZrMxePDggL8+k9pLcSH+FBPiL1QxEVOb4EKpumdLRz2vF7KyzKrv3r3QsiUkJZnrn38Of/0rlJaaEoh77oGf/cy3kQ0wJ7rt3AnNm8P555s2Z8fefhqKi4tJSkoCoKioSIdeiIiISEDVzdeivgZYgqi4GL78Ej7+GA4dgowMk/wWFMCkSfDiiyb57doVXnjB1PMem9zm5PD/7d15XNTV3gfwz4gwrCKyKCMCChmobImPK6ImJi7lnog3MsFSS5GriZZLVohyc8/cuGoadh8fw8iw9HoVTEVxy5Rr7l4UcUUQlGWY8/wxD/Mwiyw6MOh83q/XvF5y5pzf78zwfY1fzpwFN26gtH17bJXLUeri8tzJL6A85U0IASEEk98XVGlpKbZu3aq1apeMG+OCNDEmSFNdxcQLPwWC9OTaNeDIESA7G2jZ8v+nNBw/rkx8Hz4EGjcGwsOBIUOUxxxXKCsDrl9XJsv9+kF4eOByQgJPZyMVIQQuX77MmCA1jAvSxJggTXUVExwBNnZPngCHDgG7dgF37yrn+jZpAhQXA6tXAwsWKJPfVq2AhARg+HD15Dc/X7lQzs0NGDwY8PFRf56IiPSirKwMn332Gdzd3SGVSuHl5YWVK1fWuH1hYSGio6Mhk8lgbm4Of39/fP/991r1hBBYsWIFvLy8IJVK4ezsjIkTJyIvL0+r7rJlyzBs2DC0bt0aEolEaz1OhRs3biA6OhrBwcFo2rQpJBIJNm3apLOuXC7HkiVL0KFDB1hZWaF58+YIDQ3F4cOH1erNnz8fEonkqQ9dr62h2LVrF9555x34+Pio9rWtrXv37mHq1KmqeKh4nx5UHDQF4MCBA099fzIyMtSut2LFCnTp0gUODg6QSqVwdXXF6NGjca7icCsN169fx3vvvQeZTAapVIqWLVti6NChtX4dhsIRYGOWna2c63v1KiCTAba2yvILF4AlS5RzgAHgzTeBd94BKu/Rp1AAN24od3zo3l25A8T/nZBHRET6N2nSJGzZsgWff/45OnXqhF9//RVTp07Fo0eP1E7RfJphw4YhMzMT8fHxaNu2LZKSkhAWFgaFQqG2K9L06dOxbNkyTJ8+HX379kVWVhbmzp2LzMxMHDlyRG1/1zVr1sDKygp9+vRROyVM06VLl/Ddd9/B398fAwYMwLZt255aNyUlBWfPnsWsWbPQp08fPHjwAPHx8QgODsahQ4fwX//1XwCAyMhI9O/fX6t9VFQULl++rPO5hiI5ORkZGRkICAiAVCrFiRMnatU+JycHQUFBaNy4MebMmYNXXnkF9+7dw/79+3VOFYiLi0Pv3r3VyipOzK1w//59hIaGws/PD3Z2drhy5Qri4+PRuXNnnDhxAq+++qqq7tmzZ9GrVy+0adMGf/vb3+Di4oJbt27h119/rdXrMChBNZKfny8AiPz8fEN35fkVFwuRkSHE118LsWyZED/8IERKihDJyUKEhQnRqJEQgBAODkJ8/rnyucqPf/xDiMWLhdiyRYjLl4VQKNQuL5fLxYkTJ4RcLjfM66MGhzFBujAuau7s2bNCIpGIuLg4tfKoqChhYWEh7t+/X2X7n3/+WQAQSUlJauUhISFCJpOpfgc3btwQJiYm4qOPPlKrl5SUJACIdevWqZWXl5er/t2+fXsRHBys8/6V62VmZgoAYuPGjVr1ioqKhImJiQgPD1crz8nJEQDElClTqnydV69eFRKJRIwdO7bKek8THBwsIiIinqltbVR+PyZPnixqm4699dZbomXLluLBgwdV1tu/f78AILZv3/5M/czKyhIAxJw5c1RlCoVC+Pv7C39/f1FcXPxM162N2n5O1DRf4xQIY5OTo1zklp6unLPbpo1ybm9ODjBzpnJvX4UC6NlTOffXz0+9/Z07ypFfX1/lyHCbNloL3UxMTPDaa69xGxtSYUyQLi97XJw/fx5hYWFo3ry56ivld955R+u42prYuXMnhBAYN26cWvm4cePw5MkT/PLLL1W2T05OhrW1NUaOHKnVPicnB0ePHgUAZGRkoLy8HAMGDFCrN2jQIADAjh071Morn9RZlZrWqzj+tmnTpmrlTZo0QaNGjWBezTeNf//73yGEQGRkZI3uZyg1fT90uXbtGlJSUhAVFQU7Ozs99kqbo6MjAKid9Jaeno7Tp08jOjq6xqeuPY+6+pxgAmwsSkuVC9pSUpTblHl4KA+wEEJ5wtvUqcqpD1ZWwF//CkyfrkyQK5SVKef6lpcDffsqH0/ZXqS0tBSrV6/mKl5SYUyQLi9zXPz+++/o1KkTMjIysGDBAuzevRsLFy5ESUmJ2ut1d3eHu7t7tdc7e/YsHB0d0aJFC7VyX19f1fPVtff29tY6slazfUXfNBObinmqZ86cqbavz0MIgR49emDz5s3YuXMnCgoKcO3aNURFRcHW1hZRUVFPbatQKLBp0yZ4enoiODi4TvtpSAcPHoQQAjKZDGFhYbC2toa5uTl69eqFI0eO6GwzefJkNG7cGE2aNMEbb7yB33777anXLy8vR0lJCc6fP4/IyEg4OTmp/eGVnp4OALCxscGAAQNgbm4Oa2trDBo0COfPn9fvi0XdfU5wDrAxuH1bucPDxYuAo6NylwcAyMsDVq5UJsaAclR36lRlncry84Fbt5RJc7dugMYHsCYhBO7evctVvKTCmCBdXua4iImJQePGjXHs2DHVKBoAhIeHq9XTTEif5v79+2jWrJlWuZWVFczMzHD//v1q27dp00arvOKaFe3btWsHADh06JDanNHDhw9DCFHtfZ6XEAI9e/ZE165dMXz4cCgUCgCAq6sr/vWvf8HT0/Opbffs2YPs7GwsXLiwxvcqrzi5tFKZEAJyuVytvKa/p/pw8+ZNAMq52r1798aOHTtQVFSEzz77DH369MHRo0dVf9jY2tpi6tSp6NWrF+zt7XHp0iUkJCSgV69e+Pnnn/HGG29oXd/Kykr1LUXbtm1x4MABtGrVSuv+48aNw8iRI/Hzzz/j1q1b+PTTTxEUFIQzZ87A2dlZb6+3rj4nGs5vlPSvrAw4e1Z5olthoXK6QsVCtowMYNUq5XHFpqbKRW6DBytPe6tQeaFbt25c6EZEVAOPHz9GWloaxo8fr5b86nLp0qUaX7eqnQJqsotATdr7+fmhZ8+eSEhIwKuvvoqQkBBkZWXhgw8+gImJyXN9dV9T6enpOHr0KObPn4+goCAUFBRg1apVCAkJwZ49exAQEKCzXWJiIho3bqw6BbY6aWlpWgvDKu7/7bffqpVdvXoV7u7uOpPm+k6OK/4ocHFxwY4dO1RTA7p27QpPT08sXrwYW7duBQAEBASovV9BQUEYOnQofHx88PHHH+tMgA8fPozS0lJcvnwZS5cuRe/evbFv3z60b99e7f5du3bFhg0bVO06dOiAgIAAfP311/jiiy/q5sXrERPgl9Xdu8DRo8D580CzZspDLQDlSW0bNgD//Kfy59atgZgY5TZmlT15opwq0aKF8rhjHXN9iYhIW15eHsrLy+Hi4qK3a9rb2+P06dNa5UVFRSgtLdU5OqzZXtfobcWWWZXbb9++He+++y5GjRoFADAzM8O0adPwz3/+Ew8fPnz2F1ED58+fx/79+7Fw4ULExsaqykNDQ9GuXTvExMRg//79Wu3u3buHlJQUDBw4UGuayNN07NgRmZmZamXvv/8+ZDIZ5s2bp1Yuk8kA6E6aK5Lj+mJvbw8A6Nu3r9q8WGdnZ/j5+eHkyZNVtm/atCkGDRqENWvW4MmTJ7CwsFB7/rXXXgMAdOnSBW+++SY8PT0xe/Zs/Pjjj2r310ye/f394ezsXO39GwomwC8RtSODV62CVXEx4O4OVMzlysoCli5VTomQSIBhw4AxY5QjwJXduaPc+9fHB+jS5f+3R6shU1NThIeHq22VQ8aNMUG6vKxx0axZM5iYmODGjRt6u6aPjw++//575ObmqiV4f/zxBwDtLa10td+2bRvkcrnaiKWu9k5OTkhNTcWdO3eQm5sLNzc3WFhYYPXq1RgxYoTeXpMuFXvOdu7cWa3c1NQUfn5+SEtL09luy5YtKC0trdXiNxsbGwQGBmqV2dvba5VX0JU0VyTH9aVieoMuQogajdJXTCeo7psDGxsbeHl54cKFC3q9f23U1ecEF8G9rMrKlKO+Uqny31u2ALNnK5NfJycgLg6IiFBPfisvdAsJUS50q2XyCyhXt3p6etbLV2X0YmBMkC4va1xYWFggODgY27dvx7179/RyzbfeegsSiQSbN29WK9+0aRMsLCyq3fN26NChKCws1NrFYfPmzZDJZFoJJ6BMhH19fWFra4s1a9agqKgIH3744fO/mCpUjJofO3ZMrbykpAQnT5586qh6YmIiZDIZQkND67R/FUlz5YdZ5T3y60Hnzp3h4uKCPXv2qE3HyMnJwe+//44uXbpU2T4vLw+7du2Cv79/tbtq3Lt3D3/88Yfa3OvQ0FBYWlpi9+7danVPnjyJ3Nzcau9fW3X1OcER4JeVjY1ylDc7W3moxeXLyvI+fYAJEwBLS/X6BQXAzZvKhW7du1e70K0qJSUlWLJkCWJiYuplixRq+BgTpMvLHBdLlixBjx490LlzZ8TGxsLT0xO3b99GSkoK1q5dCxsbGwBQJRbVzQVu3749xo8fj3nz5sHExASdOnXCnj17sG7dOnzxxRdqUxgWLFiABQsWYN++fardEEJDQxESEoKJEyeioKAAnp6e2LZtG3755Rds3bpV7av09evXAwA8PDzw8OFD7N69G4mJiYiLi1N9PV7h+PHjuHbtGgCgoKAAQgj8z//8DwCgU6dOcKs0va6i/MqVK6q2Fd9aVowsd+rUCS4uLpg/fz4eP36Mnj17Ij8/HytXrsTVq1exZcsWrffm6NGjOHfuHGbPnv3CbKl3/fp11Ujy5f/7/7ni/XF3d1eNQF+/fh0eHh6IiIhAYmIiAGVCuHTpUowaNQpvvfUWJk6ciKKiInz++ecwMzPDrFmzVPcZM2YMXF1dERgYCAcHB1y8eBFfffUVbt++rXYSX35+PkJCQjBmzBi88sorsLCwwIULF7B8+XKUlJSoTQlp2rQpFixYgOnTp+Pdd99FWFgYcnNzMWfOHLi6umLSpEl6fa/q7HOidtsRG68X4SCMwsJCAUAAEIXr1gkRFSWEmZnyUAsbGyFiY7UPtdi5U3kgxvLlQhw8KMSTJ8/dj+LiYjF//vx62SCbXgyMCdLlZY+LrKwsMXLkSGFvby/MzMyEq6urePfdd9Ver5ubm3Bzc6vR9UpLS8W8efOEq6urMDMzE23bthUrVqzQqjdv3jwBQOzfv1+t/NGjR2LKlCmiRYsWwszMTPj6+opt27ZptV+7dq3w9vYWlpaWwtraWgQFBYmdO3fq7FNERITq/x3Nh+ZBF0+rVzkVKS4uFrGxsWLmzJmqPjg5OYlevXqJ1NRUnX2IiooSEolEXL58uZp3sHr1dRDGxo0bn/peVL7/1atXtcoq7Ny5U3Tq1EmYm5sLW1tb8eabb4pz586p1Vm4cKHw9/cXtra2wsTERDg6OoqhQ4eKY8eOqdUrLi4WkZGRwtvbW1hbW4vGjRsLFxcXMXbsWK1rVli/fr3o0KGDMDMzE/b29iI8PFxkZ2c/93ujqbafEzXN1yRCvIT7z9SBgoIC2NraIj8/H02esv+toVXMAZYB+E+7djDJylI+8dprwJQpysVwlRUXA9euAc2bK3d58PDQy0K3kpISxMfHIzY29qUb1aFnw5ggXRgXpIkxQZpqGxM1zdc4BeIlMxLAGkCZ/JqZAe+9B4SGaie2d+8q9wGuWOimceoOERER0cuKI8A19CKMAD8+fhzmnTqhEYByd3eYfPwxoLlgQC4Hrl8HLCyAzp2BDh2URyHrkUKhwL179+Dg4PDSLW6hZ8OYIF0YF6SJMUGaahsTHAE2QsLbG4sBlAL4a2wsrDS3ZqlY6NamjXKhmx5PaqlMIpHA1ta2Rhuzk3FgTJAujAvSxJggTXUVE/zzqgErKiqq9WMWgHkAisrLUVRcrHw8foyiS5dQlJODIj8/FPXujaImTbTa6ktpaSni4+P1fm43vbgYE6QL44I0MSZIU13FBBPgBsza2rpWj+bNm6vaNp88GdajRikfo0fDOiYG1h9/DOv+/WHt4KCzPREREQB89913CAgIgLm5ORwcHDBmzBhkZ2fXqK0QAitWrICXlxekUimcnZ0xceJE5OXladWVSCRqD3Nzc8yfPx8JCQlq9Xr16qVVt/IjNzdXL6/7eRUWFiI6OhoymQzm5ubw9/fH999/X6O2P/zwA8LCwuDp6QkLCwu4u7sjPDwcFy9e1KpbWlqKuXPnonXr1jAzM4ObmxtmzZqFJ0+e6Lz22bNnMXLkSDg6OkIqlcLd3V1ru7Jz585h0qRJ6Nq1K6ysrCCRSHDgwAGd14uMjESHDh3QtGlTWFhYoG3btpgxY4be9r2uD5wCQURERCorV67ElClTEBkZifj4eNy4cQNz5sxBUFAQTp06BTs7uyrbT58+HcuWLcP06dPRt29fZGVlYe7cucjMzMSRI0e0TvQaMWIE/vrXvwJQJnbffvstwsPD1eqsXr0aBQUFamWPHz9G//790bFjxxoff1zXhg0bhszMTMTHx6Nt27ZISkpCWFgYFAoFxowZU2XbRYsWoUWLFvjkk0/Qpk0bZGdnq/ZezsjIQPv27VV1w8LCkJqairlz56JTp044cuQIvvjiC5w7dw4pKSlq192/fz8GDhyIoKAgrFmzBg4ODvjPf/6DU6dOqdU7fvw4du7ciYCAALz++uv46aefntrXoqIiTJgwAZ6enjA3N8fx48fx5ZdfIjU1FadOnar3w0GeyXNuz2Y0DLEPcGFhYa0et2/fVu0jeHvGDFG4fLkoPHRIFD58WKP2+vKy7+1JtceYIF0YFw1PcXGxsLW1FYMHD1YrP3z4sAAgZs+eXWX7GzduCBMTE/HRRx+plSclJQkAYt26dWrlAMTkyZPV7l/TmNi0aZMAIDZs2FBtXV2gY6/i5/Hzzz8LACIpKUmtPCQkRMhkMiGXy6tsf/v2ba2ymzdvClNTUzF+/HhV2ZEjRwQA8dVXX6nVjYuLEwDEnj17VGVFRUXC2dlZDBw4UCgUiirvX15ervr39u3bde4lXZXVq1cLAGLfvn01blMTdbUPMKdANGBWVla1fqjaenvDasQIWHXrBitb21q1fV5mZmaIjY19Mf4CpHrBmCBdjD0u5s+fD4lEgjNnzmDkyJGwtbVFs2bNEBMTA7lcjj///BP9+/eHjY0N3N3dsXjx4jrv09mzZ5Gfn48BAwaolXft2hXNmjXTOkpZU0ZGBsrLy7XaDxo0CACqbV+bmEhMTIS1tTXefvvtauvWh+TkZFhbW2PkyJFq5ePGjUNOTg6OHj1aZXsnJyetMplMBhcXF7XpJ4cOHQKAGr3H27dvx61btzBjxoxqF5E9764bjo6OAIDGet5Zqq4+J5gAv6wGDwY0d4GoJ0II5OfnQ3CHPfo/jAnShXGhNGrUKPj5+WHHjh2IiorC0qVLMW3aNAwZMgQDBw5EcnIy+vTpg5kzZ+KHH35Qa1sxN1ZfKhYa6TpwQCqV4uLFiyguLq51e1NTU1WyrykpKQkWFhaQSqUIDAzEmjVrqo2Jixcv4uDBgxg9enSDWcNy9uxZeHt7ayWAvr6+qudr68qVK7h+/bra9IenvccVP1d+j9PT0wEA5eXl6NGjB8zMzGBnZ4ewsDDk5OTUuj+a5HI5ioqKcOjQIcyZMwc9evRA9+7dn/u6ldXV5wQT4JeVhYXBbl1WVoZvvvkGZWVlBusDNSyMCdKFcaE0YcIEfPrpp+jbty8WLVoEf39/rFq1CnFxcfjoo4/Qt29frFu3Do6Ojvjuu+/U2pqYmMDExERvfXn11VfRqFEj1ShjhcuXL+PWrVtQKBQ6F7NVaNeuHQBotT98+DCEELh//75a+ZgxY7Bq1Srs2bMHSUlJcHR0RHR0ND755JMq+5mYmAgAGD9+fI1el0KhgFwuV3voKi8vL6/R9XS5f/8+mmmeuAqoyjRfe3XkcjnGjx8Pa2trTJs2TVX+tPf4t99+07rPzZs3AQDDhw9H9+7d8euvvyI+Ph579+5FcHAwHj9+XKs+VZaRkQFTU1NYW1ujR48eaNOmDVJTU/Uaj0DdfU4wASYiIjKgiq+uK3h7e0MikSA0NFRV1rhxY3h6euL69etqdfft26dK5vShWbNmCA8Px7fffou1a9fiwYMHOHPmDMLDw1WJTVVflfv5+aFnz55ISEjA9u3b8fDhQxw+fBgffPABTExMtNp+9913GDNmDIKCgjB8+HD8+OOPaNu2Lf72t7/h7t27Ou8hl8uxefNmtG/fHl26dKnR61qwYAFMTU3VHoAyga5c5uHhoWrzLMlxVaPxtRmpF0Jg/PjxOHjwIL799lu0atVK9VxoaCg8PT0xc+ZM7N27Fw8fPsQvv/yC2bNna73HCoUCAPD2229j0aJF6N27N95//30kJibi0qVLSEpKqnGfNPn4+CAzMxNpaWlYvnw5Tp06hZCQkOdKqusTE2AiIiID0hw1NDMzg6WlJczNzbXKq5p+oC/ffPMN3n77bUyaNAn29vYICAiAl5cXBg4cCKlUCnt7+yrbb9++Hd27d8eoUaNgZ2eH3r17Y9iwYfD390fLli2rvb+vry/kcjmOHz+u8/nU1FTk5uYiMjKyxq9pwoQJyMzMVHsAwLx589TKKu98oJk0V06OdbG3t9c5yvvgwQMA2r/npxFCIDIyElu3bsWmTZvw1ltvqT1vZmaG3bt3w9XVFf369YOdnR1GjBiB2bNnw87OTu09rvhdvfHGG2rXeOONNyCRSHDy5Mka9UkXKysrBAYGomfPnpgyZQqSk5Nx9OhRrF279pmvWZ+4DRrVCWNd1EJPx5ggXRgXDY+VlRW2bNmCFStWIDs7GzKZDA4ODvDy8kK3bt2qXeTk5OSE1NRU3LlzB7m5uXBzc4OFhQVWr16NESNGVHv/6kaaExMTYWZmhr/85S81fk0ymQwyHeti3N3dERgYqLPNhAkT1Ebndc2LrszHxwfbtm2DXC5Xe4/++OMPAECHDh2q7WdF8rtx40YkJiZi7NixOut5enriyJEjuHnzJh48eAAPDw/k5+dj6tSp6Nmzp6qer69vlfsQ6/O46cDAQDRq1AgXLlzQ2zUr1MXnBEeASe+kUilmzZpV7YcFGQ/GBOnCuGjY7Ozs4OvrCwcHB6SkpODPP//E1KlTa9zeyckJvr6+sLW1xZo1a1BUVIQPP/ywyjZSqRSPHj2CqakpOnbsqPV8bm4uUlNTMWTIkGpHop+XTCZDYGCg6uHj41Nl/aFDh6KwsFBrp4vNmzdDJpOhc+fOVbYXQiAqKgobN27E2rVrMW7cuGr72LJlS/j4+MDS0hIJCQmwsrJSmxc9dOhQSCQS7N69W63d7t27IYSo8RSSmkhLS4NCoYCnp6fergnU3ecER4BJ7xQKBa5cuYI2bdro9a9LenExJkgXxsXze/3115GWlqbXecA7duxATk4OvL29UVxcjAMHDmD58uX44IMPtL6Or0h2Ll26pCpbv349AMDDwwMPHz7E7t27kZiYqDrUoUJCQgKysrLw+uuvw8XFBXfu3MGGDRuwd+9ezJs3Dw4ODlp927x5M+Ryea2mP9SX0NBQhISEYOLEiSgoKICnpye2bduGX375BVu3blVbHDZ+/Hhs3rwZly9fhpubGwBgypQpSExMxHvvvQcfHx9kZGSo6kulUgQEBKh+Xrx4MVq0aAFXV1fcvn0b//3f/42dO3diy5YtalMgvLy8MHnyZKxevRo2NjYIDQ3FhQsX8OmnnyIgIACjRo1S1X38+DFSU1MBQHXvtLQ03Lt3D1ZWVqo56bt27cL69evx5ptvws3NDWVlZTh+/DiWLVsGT09Pvf9u6uxz4pl2JTZChjgI40XFze1JE2OCdDH2uJg3b54AIO7evatWHhERIaysrLTqBwcHi/bt22uV6fu/8uTkZOHv7y+srKyEhYWFCAwMFImJiToPUnBzcxNubm5qZWvXrhXe3t7C0tJSWFtbi6CgILFz506ttikpKaJHjx7C0dFRNG7cWNjY2Iju3buL4cOHPzUm2rZtK9zd3as91KEmoOeDMIQQ4tGjR2LKlCmiRYsWwszMTPj6+opt27Zp1YuIiBAAxNWrV1Vlbm5uqsOsNB+a7/Fnn30mPDw8hFQqFU2bNhX9+/cX6enpOvskl8tFfHy88PT0FKampsLZ2VlMnDhR5OXlqdW7evVqje7/73//W4wYMUK4ubkJc3NzYW5uLry8vMSMGTPE/fv3n/Wte6q6OghDIoSRb8BYQwUFBbC1tUV+fj6aNGli6O40aCUlJYiPj0dsbCy/2iQAjAnSjXFBmhgTpKm2MVHTfI3fORERERGRUWECTHonkUjg6Oio19OJ6MXGmCBdGBekiTFBmuoqJjgFooY4BYKIiIioYeMUCDKY8vJynDx58rmOlKSXC2OCdGFckCbGBGmqq5hgAkx6J5fL8dNPP+l1Wx56sTEmSBfGBWliTJCmuooJJsBEREREZFSYABMRERGRUWECTHonkUjg4eHBVbykwpggXRgXpIkxQZrqKia4C0QNcRcIIiIiooaNu0CQwcjlchw4cICLGEiFMUG6MC5IE2OCNNVVTDABJr0rLy9HWloat7EhFcYE6cK4IE2MCdJUVzHBBJiIiIiIjAoTYCIiIiIyKkyASe8aNWqEgIAANGrE8CIlxgTpwrggTYwJ0lRXMcFdIGqIu0AQERERNWzcBYIMpqysDCkpKSgrKzN0V6iBYEyQLowL0sSYIE11FRNMgEnvFAoFTp06BYVCYeiuUAPBmCBdGBekiTFBmuoqJpgAExEREZFRaWzoDrwoKqZKFxQUGLgnDV9JSQmKi4tRUFAAqVRq6O5QA8CYIF0YF6SJMUGaahsTFXladUvcuAiuhm7cuIFWrVoZuhtEREREVI3s7Gy4uLg89XkmwDWkUCiQk5MDGxsbSCQSQ3enQSsoKECrVq2QnZ3NHTMIAGOCdGNckCbGBGmqbUwIIfDo0SPIZLIqt07jFIgaatSoUZV/SZC2Jk2a8AOM1DAmSBfGBWliTJCm2sSEra1ttXW4CI6IiIiIjAoTYCIiIiIyKkyASe+kUinmzZvHFbykwpggXRgXpIkxQZrqKia4CI6IiIiIjApHgImIiIjIqDABJiIiIiKjwgSYiIiIiIwKE2AiIiIiMipMgElvFi5ciE6dOsHGxgZOTk4YMmQI/vzzT0N3ixqQhQsXQiKRIDo62tBdIQO6efMmxo4dC3t7e1haWsLf3x8nTpwwdLfIQORyOT799FO0bt0aFhYWaNOmDRYsWACFQmHorlE9SU9Px+DBgyGTySCRSLBz506154UQmD9/PmQyGSwsLNCrVy+cO3fuue7JBJj0Ji0tDZMnT0ZGRgb27t0LuVyOfv36oaioyNBdowYgMzMT69atg6+vr6G7QgaUl5eH7t27w9TUFLt370ZWVha++uorNG3a1NBdIwNZtGgR1qxZg1WrVuHf//43Fi9ejISEBKxcudLQXaN6UlRUBD8/P6xatUrn84sXL8aSJUuwatUqZGZmokWLFggJCcGjR4+e+Z7cBo3qzN27d+Hk5IS0tDT07NnT0N0hAyosLMRrr72G1atX44svvoC/vz+WLVtm6G6RAcTGxuLQoUM4ePCgobtCDcSgQYPQvHlzJCYmqsqGDx8OS0tLbNmyxYA9I0OQSCRITk7GkCFDAChHf2UyGaKjozFz5kwAQElJCZo3b45Fixbh/ffff6b7cASY6kx+fj4AoFmzZgbuCRna5MmTMXDgQPTt29fQXSEDS0lJQWBgIEaOHAknJycEBARg/fr1hu4WGVCPHj2wb98+XLhwAQDw+++/47fffsOAAQMM3DNqCK5evYrc3Fz069dPVSaVShEcHIzDhw8/83Ub66NzRJqEEIiJiUGPHj3QoUMHQ3eHDOj777/HyZMnkZmZaeiuUANw5coVfPPNN4iJicHs2bNx7NgxTJkyBVKpFO+8846hu0cGMHPmTOTn58PLywsmJiYoLy/Hl19+ibCwMEN3jRqA3NxcAEDz5s3Vyps3b47r168/83WZAFOd+PDDD3HmzBn89ttvhu4KGVB2djamTp2KPXv2wNzc3NDdoQZAoVAgMDAQcXFxAICAgACcO3cO33zzDRNgI/WPf/wDW7duRVJSEtq3b4/Tp08jOjoaMpkMERERhu4eNRASiUTtZyGEVlltMAEmvfvoo4+QkpKC9PR0uLi4GLo7ZEAnTpzAnTt30LFjR1VZeXk50tPTsWrVKpSUlMDExMSAPaT65uzsjHbt2qmVeXt7Y8eOHQbqERnajBkzEBsbi9GjRwMAfHx8cP36dSxcuJAJMKFFixYAlCPBzs7OqvI7d+5ojQrXBucAk94IIfDhhx/ihx9+wL/+9S+0bt3a0F0iA3v99dfxxx9/4PTp06pHYGAgwsPDcfr0aSa/Rqh79+5a2yNeuHABbm5uBuoRGdrjx4/RqJF6OmJiYsJt0AgA0Lp1a7Ro0QJ79+5VlZWWliItLQ3dunV75utyBJj0ZvLkyUhKSsKPP/4IGxsb1bwdW1tbWFhYGLh3ZAg2NjZac8CtrKxgb2/PueFGatq0aejWrRvi4uIwatQoHDt2DOvWrcO6desM3TUykMGDB+PLL7+Eq6sr2rdvj1OnTmHJkiV47733DN01qieFhYW4dOmS6uerV6/i9OnTaNasGVxdXREdHY24uDi88soreOWVVxAXFwdLS0uMGTPmme/JbdBIb542F2fjxo14991367cz1GD16tWL26AZuV27dmHWrFm4ePEiWrdujZiYGERFRRm6W2Qgjx49wpw5c5CcnIw7d+5AJpMhLCwMc+fOhZmZmaG7R/XgwIED6N27t1Z5REQENm3aBCEEPvvsM6xduxZ5eXno3Lkzvv766+caSGECTERERERGhXOAiYiIiMioMAEmIiIiIqPCBJiIiIiIjAoTYCIiIiIyKkyAiYiIiMioMAEmIiIiIqPCBJiIiIiIjAoTYCIiIiIyKkyAiYiIiMioMAEmIiIiIqPCBJiIiIiIjAoTYCIiIiIyKkyAiYiIiMioMAEmIjISkZGRkEgkCAkJgRBC6/m5c+dCIpHAx8cHJSUlBughEVH9kAhdn4JERPTSKSwshJ+fH65cuYKlS5ciOjpa9dzRo0fRvXt3mJiY4NixY/Dz8zNcR4mI6hhHgImIjIS1tTW2bNkCExMTzJo1C+fOnQMAPH78GH/5y19QXl6Ozz//nMkvEb30mAATERmRbt264eOPP0ZxcTHGjh2L0tJSxMTE4OLFi+jZsyemT59u6C4SEdU5ToEgIjIyZWVl6Ny5M06dOoWQkBDs3bsXTZo0wZkzZ+Dm5mbo7hER1TkmwERERigrKwsdO3ZEcXExAGDTpk2IiIgwcK+IiOoHE2AiIiNUWloKHx8fXLhwAba2trhx4wasra0N3S0ionrBOcBEREbok08+wYULF9CoUSPk5+dj2rRphu4SEVG9YQJMRGRk0tPTsWTJElhaWmLv3r1o2rQpNmzYgJ9++snQXSMiqhdMgImIjEhBQQEiIiKgUCiQkJCAPn364OuvvwagPCjj7t27Bu4hEVHdYwJMRGREpkyZgmvXrqFfv36YNGkSAGDMmDF4++23cefOHUyYMMHAPSQiqntcBEdEZCSSk5MxbNgw2NnZ4ezZs5DJZKrn8vLy0KFDB+Tk5ODvf/87xo0bZ8CeEhHVLSbARERG4Pbt2+jQoQPu3buHbdu2YfTo0Vp19uzZg/79+8Pa2hpnzpyBu7t7/XeUiKgeMAEmIiIiIqPCOcBEREREZFSYABMRERGRUWECTERERERGhQkwERERERkVJsBEREREZFSYABMRERGRUWECTERERERGhQkwERERERkVJsBEREREZFSYABMRERGRUWECTERERERGhQkwERERERkVJsBEREREZFT+F4tmGgUEcb2QAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize = (8, 6))\n",
    "plt.title('Fit result')\n",
    "plt.xlabel('x', fontsize=16)\n",
    "plt.ylabel('y', fontsize=16)\n",
    "plt.grid(color = 'grey', linestyle=\"--\")\n",
    "\n",
    "# Plot data\n",
    "plt.errorbar(xdata, ydata, xerr = xerror, yerr = yerror, fmt='k', linestyle = '', label = \"Data\") \n",
    "\n",
    "# Plot fit \n",
    "plt.plot(xdata, yfit, color = 'r', linestyle = '-', label = \"Fit\")\n",
    "\n",
    "# Plot fit error band\n",
    "plt.fill_between(x_fit, y_fit - y_err, y_fit + y_err,\n",
    "                 color='red', alpha=0.3, label=\"1σ fit band\")\n",
    "\n",
    "# Add legend and fit params\n",
    "plt.legend(loc = 2, fontsize=16)    \n",
    "\n",
    "text = \"c: {:7.5g} +- {:7.5g}\\n\".format(final_params[0], param_errors[0])\n",
    "text += \"m: {:7.5g} +- {:7.5g}\\n\".format(final_params[1], param_errors[1])\n",
    "plt.text(0.95, 0.0, text, transform = fig.axes[0].transAxes, ha = \"right\", va = \"bottom\", fontsize=12)        \n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "a47bfce9-c0fb-4b99-a2de-3a1545c0588e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.odr import ODR, Model, RealData\n",
    "\n",
    "def curve_fit_with_xerr(func, xdata, ydata, p0=None, xerr=None, yerr=None):\n",
    "    \"\"\"\n",
    "    Fit function with uncertainties in both x and y using ODR.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    func : callable\n",
    "        The model function, must be callable as func(x, *params).\n",
    "    xdata, ydata : array-like\n",
    "        Data points.\n",
    "    p0 : list or tuple, optional\n",
    "        Initial guess for parameters.\n",
    "    xerr, yerr : array-like, optional\n",
    "        Standard deviations of x and y.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    popt : array\n",
    "        Optimal parameters.\n",
    "    perr : array\n",
    "        Standard deviations of the parameters.\n",
    "    \"\"\"\n",
    "    # ODR model expects parameter array first\n",
    "    def odr_model(beta, x):\n",
    "        return func(x, *beta)\n",
    "\n",
    "    # Default errors\n",
    "    sx = xerr if xerr is not None else np.zeros_like(xdata)\n",
    "    sy = yerr if yerr is not None else np.ones_like(ydata)\n",
    "    \n",
    "    # Build data & model\n",
    "    data = RealData(xdata, ydata, sx=sx, sy=sy)\n",
    "    model = Model(odr_model)\n",
    "    \n",
    "    # Initial guess\n",
    "    if p0 is None:\n",
    "        p0 = np.ones(func.__code__.co_argcount - 1)  # crude default\n",
    "    \n",
    "    # Run ODR\n",
    "    odr = ODR(data, model, beta0=p0)\n",
    "    output = odr.run()\n",
    "    \n",
    "    popt = output.beta\n",
    "    perr = output.sd_beta\n",
    "    return popt, perr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "ccbe0609-a059-40dc-86dd-8838b8c12bf3",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-1.535 10.243]\n",
      "[1.578 0.289]\n",
      "{'fvec': array([ 0.346,  1.701, -0.706, -0.693, -2.205,  0.351,  2.1  ,  1.321,\n",
      "       -0.81 , -1.112]), 'nfev': 7, 'fjac': array([[-13.728,   0.099,   0.061,   0.237,   0.33 ,   0.315,   0.273,\n",
      "          0.568,   0.513,   0.242],\n",
      "       [ -2.295,  -0.882,   0.16 ,   0.323,   0.274,   0.089,  -0.059,\n",
      "         -0.437,  -0.587,  -0.333]]), 'ipvt': array([2, 1], dtype=int32), 'qtf': array([ 3.189e-09, -1.874e-07])}\n",
      "Both actual and predicted relative reductions in the sum of squares\n",
      "  are at most 0.000000\n",
      "1\n"
     ]
    }
   ],
   "source": [
    "popt, pcov = curve_fit_with_xerr(straight_line, xdata, ydata, p0 = init_params, xerr = xerror, yerr = yerror)\n",
    "print(popt)\n",
    "print(pcov)\n",
    "print(info)\n",
    "print(msg)\n",
    "print(ierr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "c6ef1f36-dc7c-4a80-b3a8-f046e7bab750",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.535, 10.243])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "popt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "a559e5d8-b5a9-49df-8d3b-cc30c52a355d",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.578, 0.289])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pcov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6fb75db8-02af-4249-a6c9-2e914b033485",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}