eteq/centroid_testing.ipynb

## centroid_testing.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ed0aba5c-d7d4-4f43-823e-83672b6fadf9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from astropy import modeling\n",
    "from astropy import units as u\n",
    "\n",
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "99da5749-8e26-40fe-bc74-ccc09854572d",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f98ba8269d0>]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA16klEQVR4nO3de3yU5Z3///fMJJkJ5AAYyAGjgBXEomChsKntt7ablVZLcbvbUusKy1r7qML+WrPdVqqC2ta425Wyjy2Wb1W0u18VelDWrTxQy6P8/FqjVA61HhA5KAgmnCSTA8kkM9f3j+SeSSSBTHLfc8995/V8PPKwmdwzc10d7vDmuj7XdQWMMUYAAAAuCbrdAAAAMLwRRgAAgKsIIwAAwFWEEQAA4CrCCAAAcBVhBAAAuIowAgAAXEUYAQAArspxuwEDkUgkdPjwYRUWFioQCLjdHAAAMADGGDU1NamiokLBYP/jH54II4cPH1ZlZaXbzQAAAINw8OBBnXvuuf3+3BNhpLCwUFJXZ4qKilxuDQAAGIhoNKrKysrk3+P98UQYsaZmioqKCCMAAHjM2UosKGAFAACuIowAAABXEUYAAICrCCMAAMBVhBEAAOAqwggAAHAVYQQAALiKMAIAAFxFGAEAAK5KO4w8//zzmjdvnioqKhQIBLRhw4azPmfLli362Mc+pnA4rI985CN65JFHBtFUAADgR2mHkZaWFk2fPl2rV68e0PX79+/X1Vdfrc985jPauXOnvv3tb+vrX/+6nnnmmbQbCwAA/Cfts2k+//nP6/Of//yAr1+zZo0mTpyo++67T5I0depUvfDCC/rJT36iuXPnpvv2AADAZxyvGamrq1N1dXWvx+bOnau6urp+n9Pe3q5oNNrrC4B9WmOdeuD5fTra1O52UwDA+TBSX1+v0tLSXo+VlpYqGo3q1KlTfT6ntrZWxcXFya/KykqnmwkMK4+8+I5+tPFN/dszb7ndFADIztU0y5YtU2NjY/Lr4MGDbjcJ8JU332+SJO04+IHLLQGAQdSMpKusrEwNDQ29HmtoaFBRUZHy8/P7fE44HFY4HHa6acCwtf9YsyRpz5FmtbR3amTY8V8FANAvx0dGqqqqtHnz5l6PPffcc6qqqnL6rQH0wRij/UdbJEkJI73xPjVZANyVdhhpbm7Wzp07tXPnTkldS3d37typAwcOSOqaYlm4cGHy+m9+85vat2+fvvvd72rXrl26//779ctf/lK33HKLPT0AkJaGaLtaYvHk96++1+hiawBgEGHklVde0WWXXabLLrtMklRTU6PLLrtMy5cvlyS9//77yWAiSRMnTtTTTz+t5557TtOnT9d9992nBx98kGW9gEv2dU/RWF47RBgB4K60J4qvuOIKGWP6/Xlfu6teccUV2rFjR7pvBcAB+7qnaEbkhdQai+vV90662yAAw15WrqYB4BwrjMz9aFnX98da1NTW4WaTAAxzhBFgmLFW0syaMFrjR+XLGOn1wxSxAnAPYQQYZvYd6xoZmVRSoGnjiyRRNwLAXYQRYBhp74zr4IlWSdKksSN16bmjJLGiBoC7CCPAMHLgeKsSRhqZF9K4wrAuGV8sSfozIyMAXEQYAYaR5BTN2AIFAoFkGNl/rEWNpyhiBeAOwggwjFgraSaNHSlJGj0yT+eO7jqW4XVGRwC4hDACDCP7jnatpJlUUpB87NJzmaoB4C7CCDCM7O+eppnYPTIiSZeMHyVJepUwAsAlhBFgGEkt602FkeTICCtqALiEMAIMEydbYzrREpOUqhmRpGkVXWHkwIlWnWyNudI2AMMbYQQYJvZ2F6+WF0c0Ii91LFXxiFydf84ISdJrh9iJFUDmEUaAYSJZL9JjisZiLfF99dDJTDYJACQRRoBhI7mSZuzpYYS6EQBuIowAw0Ryj5Eey3ot06yREcIIABcQRoBhYv+x3hue9WSFkUMnTyWLXAEgUwgjwDAQTxjtP97/yEhRJFfjCsOSpMMnT2W0bQBAGAGGgcMnTynWmVBeTlDju7d//7DRI/IkSSdbOaMGQGYRRoBhYG938eqEc0YoFAz0eU3xiFxJ0slTTNMAyCzCCDAMJOtF+piisYzK7w4jjIwAyDDCCDAMWCtpJvZRvGoZ1T0y0niKMAIgswgjwDCw75h1Wu+ZwkhXzQhhBECmEUaAYWC/tcfI2P6naYqT0zTUjADILMII4HPtnXEdbmyT1PdW8JZiakYAuIQwAvhcY3e4CAak0d11IX0ZlVxNQxgBkFmEEcDnrBqQovxcBQJ9L+uVpFH53TUjjIwAyDDCCOBzVhixpmH6M4p9RgC4hDAC+NxAwwg1IwDcQhgBfC7dkZH2zoTaOuKOtwsALIQRwOeiPWpGzqQgnJPcKp7REQCZRBgBfK7xVKekrpN5zyQQCKS2hKduBEAGEUYAnxvoNI3U47A8RkYAZBBhBPC5dMIIh+UBcANhBPC5tMJI9/k0UTY+A5BBhBHA56JtgxgZoWYEQAYRRgCfi1IzAiDLEUYAn0ttB59z1muTG58xTQMggwgjgM8NpoCV82kAZBJhBPCxjnhCrbGu3VTTKWClZgRAJhFGAB/ruSqm8CybnknUjABwB2EE8DFriqYwktrq/UzYZwSAGwgjgI8li1cHMCoipaZpGilgBZBBhBHAx9IpXpVSIyPN7Z3qiCccaxcA9EQYAXws3TDS82RfRkcAZAphBPCxdDY8k6RQMKCiSNd+JNSNAMgUwgjgY9G2TkkDDyMSdSMAMo8wAvhYOruvWkZ1L+9tZK8RABlCGAF8zNpJNZ2RkWKW9wLIMMII4GPpFrBKPXZhJYwAyBDCCOBjqWmadEZGugtYqRkBkCGEEcDHom2DGBnJ7y5gbaVmBEBmEEYAHxvMyIhVwMrICIBMIYwAPjaYmhEKWAFkGmEE8Kl4wqhpCPuMMDICIFMII4BPNbWlwkR6YaR7nxFqRgBkCGEE8ClrimZEXki5oYHf6tZheYyMAMiUQYWR1atXa8KECYpEIpozZ462bt16xutXrVqlKVOmKD8/X5WVlbrlllvU1tY2qAYDGJjoqfSnaCSpOLkDa4cSCWN7uwDgw9IOI+vXr1dNTY1WrFih7du3a/r06Zo7d66OHDnS5/WPPfaYbr31Vq1YsUJvvvmmHnroIa1fv17f//73h9x4AP1LrqSJpBlGusOLMUrWnACAk9IOIytXrtSNN96oxYsX6+KLL9aaNWs0YsQIrV27ts/rX3zxRV1++eX62te+pgkTJujKK6/Utddee9bRFABDM5iVNJIUzglpRF6o12sAgJPSCiOxWEzbtm1TdXV16gWCQVVXV6uurq7P53ziE5/Qtm3bkuFj37592rhxo6666qohNBvA2QxmjxFLqm6EIlYAzhv4UZ6Sjh07png8rtLS0l6Pl5aWateuXX0+52tf+5qOHTumT37ykzLGqLOzU9/85jfPOE3T3t6u9vb25PfRaDSdZgLQ4EdGJKl4RJ4ON7ax1wiAjHB8Nc2WLVt0zz336P7779f27dv1xBNP6Omnn9YPfvCDfp9TW1ur4uLi5FdlZaXTzQR8ZzBbwVtYUQMgk9IaGSkpKVEoFFJDQ0OvxxsaGlRWVtbnc+644w5df/31+vrXvy5JuuSSS9TS0qJvfOMbuu222xQMnp6Hli1bppqamuT30WiUQAKkKTVNk9ZtLom9RgBkVlojI3l5eZo5c6Y2b96cfCyRSGjz5s2qqqrq8zmtra2nBY5QqKs4zpi+lw2Gw2EVFRX1+gKQniFN07AlPIAMSvufTDU1NVq0aJFmzZql2bNna9WqVWppadHixYslSQsXLtT48eNVW1srSZo3b55Wrlypyy67THPmzNGePXt0xx13aN68eclQAsB+0SHVjDBNAyBz0g4jCxYs0NGjR7V8+XLV19drxowZ2rRpU7Ko9cCBA71GQm6//XYFAgHdfvvtOnTokMaOHat58+bpRz/6kX29AHCaoYyMjMrvPp+GkREAGZD+ZLKkpUuXaunSpX3+bMuWLb3fICdHK1as0IoVKwbzVgAGaSgjI8maEZb2AsgAzqYBfMqWfUYYGQGQAYQRwIeMMYq2De5sGomaEQCZRRgBfKi5vVPx7kPuqBkBkO0II4APWVM0eTlBRXLTX7Vm1YxET3X0uwQfAOxCGAF8aCgraaRUGInFEzrVEbetXQDQF8II4EPRU131IkWRQS2YU35uSHmhrl8PTNUAcBphBPChoY6MBAKBVBErYQSAwwgjgA8NZY8RS+qwPPYaAeAswgjgQ0MdGZF6HpbHyAgAZxFGAB+yI4wkD8tjrxEADiOMAD4UbbMjjLDXCIDMIIwAPjSUreAto0ZQMwIgMwgjgA/ZEkbyqRkBkBmEEcCH7Cxg/aCVkREAziKMAD5kRxixRlWaug/cAwCnEEYAH7Jjn5HC7t1bm9sJIwCcRRgBfMYYk9oOfghhpCDMyAiAzCCMAD7T1pFQLJ6QZM/ICGEEgNMII4DPWPUioWBAI/NCg36dgrA1TcNqGgDOIowAPtOzeDUQCAz6dayRkbaOhDq6R1oAwAmEEcBn7FhJI6VGRiSpmakaAA4ijAA+Y62kKYrknOXKM8sJBZWf2zXNw4oaAE4ijAA+Y8fuqxZrqsY66wYAnEAYAXzGrmkaSSqw9hphmgaAgwgjgM/YGUYKw2x8BsB5hBHAZ2wNIxE2PgPgPMII4DNWcLCCxFBYK2qaGBkB4CDCCOAz1iZlhUNcTSOlakaaKGAF4CDCCOAzVn2HHWGkkAJWABlAGAF8JjVNY0MYoYAVQAYQRgCfsUYxrFN3h6KAw/IAZABhBPCZaDKM2DFNw2oaAM4jjAA+Y2sBKyf3AsgAwgjgIx3xhNo6uk7YtXc1DSMjAJxDGAF8pOeql5E2TNNYh+1RwArASYQRwEes0JCfG1JuaOi3t1UEy9JeAE4ijAA+Yk2nFNgwRSOlpnqYpgHgJMII4CPWTqmFNkzRSKlQE4sn1N4Zt+U1AeDDCCOAj9i5+6okjcxLvQ5TNQCcQhgBfMQKI3ZN04SCgdRheYQRAA4hjAA+Ym14VmjD7quWAraEB+AwwgjgI802F7D2fK0oJ/cCcAhhBPARq4DVjq3gLZzcC8BphBHAR6yplCI7R0aYpgHgMMII4CNOTNOw1wgApxFGAB9JndhrXwGrVQzLyAgApxBGAB+x88ReC4flAXAaYQTwEbv3GZHUY58RVtMAcAZhBPCRpuQ+Iw6spmGaBoBDCCOAj1gFrIURG2tGWNoLwGGEEcBHmhyYprGCDTUjAJxCGAF8or0zrlhnQpLNBaxWzQjTNAAcQhgBfKLnNErP03aHqiBZM0IBKwBnEEYAn7CmUUbmhRQKBmx73SKW9gJwGGEE8AlrtYudxatSagO15rZOGWNsfW0AkAgjgG80ObAVfM/X60wYtXUkbH1tAJAII4BvOHFir9Q17RPonvVpom4EgAMGFUZWr16tCRMmKBKJaM6cOdq6desZrz958qSWLFmi8vJyhcNhTZ48WRs3bhxUgwH0LTVNY28YCQQCqZN7qRsB4IC0f2utX79eNTU1WrNmjebMmaNVq1Zp7ty5euuttzRu3LjTro/FYvqrv/orjRs3Tr/+9a81fvx4vfvuuxo1apQd7QfQzakwInXt6NrU1kkRKwBHpP1ba+XKlbrxxhu1ePFiSdKaNWv09NNPa+3atbr11ltPu37t2rU6ceKEXnzxReXmdhXCTZgwYWitBnCaZM2IzdM0UndRbGMbW8IDcERa0zSxWEzbtm1TdXV16gWCQVVXV6uurq7P5zz11FOqqqrSkiVLVFpaqmnTpumee+5RPB7v933a29sVjUZ7fQE4syYHtoK3cHIvACelFUaOHTumeDyu0tLSXo+Xlpaqvr6+z+fs27dPv/71rxWPx7Vx40bdcccduu+++/TDH/6w3/epra1VcXFx8quysjKdZgLDkrUpmRMjI5zcC8BJjq+mSSQSGjdunH7+859r5syZWrBggW677TatWbOm3+csW7ZMjY2Nya+DBw863UzA81IjI05M03ByLwDnpPVbq6SkRKFQSA0NDb0eb2hoUFlZWZ/PKS8vV25urkKhUPKxqVOnqr6+XrFYTHl5eac9JxwOKxwOp9M0YNhrzkQYYZoGgAPSGhnJy8vTzJkztXnz5uRjiURCmzdvVlVVVZ/Pufzyy7Vnzx4lEqnNknbv3q3y8vI+gwiAwUkVsNpfM5I8uZeREQAOSHuapqamRg888IB+8Ytf6M0339RNN92klpaW5OqahQsXatmyZcnrb7rpJp04cULf+ta3tHv3bj399NO65557tGTJEvt6ASAZFJwYGUnVjBBGANgv7d9aCxYs0NGjR7V8+XLV19drxowZ2rRpU7Ko9cCBAwoGUxmnsrJSzzzzjG655RZdeumlGj9+vL71rW/pe9/7nn29AJAqYHUwjFAzAsAJg/qttXTpUi1durTPn23ZsuW0x6qqqvTSSy8N5q0ADFCygNWRfUZYTQPAOZxNA/iAMaZHAasTNSMUsAJwDmEE8IH2zoQ6E0aSU9M03QWshBEADiCMAD4Q7Z4+CQSkEbmhs1ydPvYZAeAkwgjgA809zqUJBgO2v34BNSMAHEQYAXwgeWKvA8WrPV+3ub1TxhhH3gPA8EUYAXzAyUPyer5uwkitsf4PuQSAwSCMAD6Q3H3VgeJVSYrkBhXqnv6hbgSA3QgjgA9YtRxOnNgrSYFAgL1GADiGMAL4QLODW8Fb2BIegFMII4APOHlir4Ut4QE4hTAC+IB1SJ5T0zSSVBRh4zMAziCMAD7g9GoaKVUcy5bwAOxGGAF8oDkDIyPJmhGmaQDYjDAC+EByNY2DNSOspgHgFMII4APW1EmRkwWsTNMAcAhhBPCB1DSNczUjhSztBeAQwgjgA00ZWNprFceytBeA3QgjgA9komaEAlYATiGMAB5njHH81F4pFXQoYAVgN8II4HGtsbgSput/O7nPSCEFrAAcQhgBPM4aFQkFA4rkOndLF4apGQHgDMII4HE9T+wNBAKOvU9qnxHCCAB7EUYAj8vEShqpxz4j7Z1KWPNCAGADwgjgcZnYCv7Dr98SY3QEgH0II4DHZWpkJJIbUl4o2Os9AcAOhBHA45ozcGKvpedUDQDYhTACeFxThqZper4HIyMA7EQYATzOWk3j9DRNz/eIsvEZABsRRgCPs6ZpnNwK3lLUPRUUPUUYAWAfwgjgcckC1gxM0xTnd4cRpmkA2IgwAnhc8lyaDBSwFuV3T9MwMgLARoQRwOMyWcCanKahZgSAjQgjgMclt4PPRM1IPjUjAOxHGAE8rjlDm55JPWpGTlEzAsA+hBHA45I1I+EM1owwTQPARoQRwOOaXFja28g0DQAbEUYAD0skTI/VNJmcpiGMALAPYQTwsJ6n52ZkNQ37jABwAGEE8DBriiYvFFQkN+T4+/WcpjHGOP5+AIYHwgjgYdEMLuuVUgWs8YRRayyekfcE4H+EEcDDGlu7wohVy+G0/NyQckMBSayoAWAfwgjgYdaqlqIMhZFAINDjsDzqRgDYgzACeJgVRjI1MiKlgg/LewHYhTACeJgrYSTCYXkA7EUYATwsmgwjmSlglXou7yWMALAHYQTwMGu/D6ZpAHgZYQTwMHemaShgBWAvwgjgYW6EkWKmaQDYjDACeJg7q2koYAVgL8II4GGZ3mdE4uReAPYjjAAe5uY+I0zTALALYQTwMFdrRihgBWATwgjgUW0dccU6E5Lc2fSMaRoAdiGMAB5lhYFQMKCCMJueAfAuwgjgUcni1UiOAoFAxt7XGoVpbu9UImEy9r4A/IswAniUG/UiklTYPU1jjNTUTt0IgKEbVBhZvXq1JkyYoEgkojlz5mjr1q0Det66desUCAR0zTXXDOZtAfTQ2OpOGAnnhBTJ7frVwV4jAOyQdhhZv369ampqtGLFCm3fvl3Tp0/X3LlzdeTIkTM+75133tF3vvMdfepTnxp0YwGkuLHHiIW9RgDYKe0wsnLlSt14441avHixLr74Yq1Zs0YjRozQ2rVr+31OPB7Xddddp7vuukuTJk0aUoMBdHFrmqbne1LECsAOaYWRWCymbdu2qbq6OvUCwaCqq6tVV1fX7/PuvvtujRs3TjfccMOA3qe9vV3RaLTXF4De3AwjyRU1jIwAsEFaYeTYsWOKx+MqLS3t9Xhpaanq6+v7fM4LL7yghx56SA888MCA36e2tlbFxcXJr8rKynSaCQwLroaRiHU+DQWsAIbO0dU0TU1Nuv766/XAAw+opKRkwM9btmyZGhsbk18HDx50sJWAN0WzYWSEaRoANkhrp6SSkhKFQiE1NDT0eryhoUFlZWWnXb9371698847mjdvXvKxRKJrx8icnBy99dZbuuCCC057XjgcVjgcTqdpwLCTFTUjTNMAsEFaIyN5eXmaOXOmNm/enHwskUho8+bNqqqqOu36iy66SH/+85+1c+fO5NcXv/hFfeYzn9HOnTuZfgGGwN1pGlbTALBP2ntI19TUaNGiRZo1a5Zmz56tVatWqaWlRYsXL5YkLVy4UOPHj1dtba0ikYimTZvW6/mjRo2SpNMeB5AedwtYu2tG2qgZATB0aYeRBQsW6OjRo1q+fLnq6+s1Y8YMbdq0KVnUeuDAAQWDbOwKOM3NfUaYpgFgp0GdrrV06VItXbq0z59t2bLljM995JFHBvOWAD6EaRoAfsEQBuBBbR1xtXd2FYMXj2A1DQBvI4wAHmRNjwQDUkHeoAY4h8QaGWGfEQB2IIwAHtSzXiQYDGT8/dkOHoCdCCOAB7lZLyKlVtO0xuLqiCdcaQMA/yCMAB6UHBmJuBNGCnu8LytqAAwVYQTwILdHRkLBgArD7DUCwB6EEcCD3A4jUmpFDct7AQwVYQTwIDc3PLMUJk/uJYwAGBrCCOBB2TQywooaAENFGAE8KBvCSGpLeGpGAAwNYQTwoGgWhBG2hAdgF8II4EHZMDKSOrmXMAJgaAgjgAdlQxjh5F4AdiGMAB6UDWGEaRoAdiGMAB6UFWEkuZqGAlYAQ0MYATymvTOuto6u82DcHRlhnxEA9iCMAB5jjYoEAqmNx9zAyb0A7EIYATzGGokoDOcoGAy41o4iClgB2IQwAnhMsl5khHtTNFLPMNIpY4yrbQHgbYQRwGOyoXhVStWMxOIJtXcmXG0LAG8jjAAeky1hpCCcI2uWiOW9AIaCMAJ4TGNrdoSRQCBA3QgAWxBGAI9p7D6Yzu0wIqU2PmNFDYChIIwAHmNNiRRlQRjh5F4AdiCMAB6TLTUjUuqwPGpGAAwFYQTwmKwKI0zTALABYQTwmGg2hhFGRgAMAWEE8JhsGhmxNl5jmgbAUBBGAI/JpjCSOiyPAlYAg0cYATzGqs/IijDCYXkAbEAYATykI55QaywuKTvCCCf3ArADYQTwkJ61GYUR98OIVcBKzQiAoSCMAB5i/aVfGMlRyDoYxkWjR+ZJkk40x1xuCQAvI4wAHpJNxauSNLYwLEk62twuY4zLrQHgVYQRwEOyLYyc0z0y0hE3TNUAGDTCCOAh2bThmSRFckPJ5b3Hmttdbg0AryKMAB6SbSMjUmqq5kgTYQTA4BBGAA9pbO0+sTcLVtJYknUjhBEAg0QYATwkOTIyInvCSElBVxg5xooaAINEGAE8JJunaRgZATBYhBHAQ6wwUkQYAeAjhBHAQ452r1gZW5DncktSrGmao6ymATBIhBHAQ45Eu8NIYcTllqRYIyPHGBkBMEiEEcAjjDHJqZDSorDLrUkZy8gIgCEijAAecbK1Q7F4QlJqNCIbjOtuy/HmdsUTbAkPIH2EEcAjGpraJEmjR+QqnBNyuTUpY0bmKRCQEkb6oJXlvQDSRxgBPMKqFyktyp56EUnKCQU1ZkRXQS0ragAMBmEE8IiGaNfISDZN0VhY3gtgKAgjgEccacrOkRGpx/JewgiAQSCMAB5xpHtkJJtW0liSy3tZUQNgEAgjgEc0dNeMjMuiPUYsTNMAGArCCOARR5qyd2SkpHtHWPYaATAYhBHAIxqycPdVC9M0AIaCMAJ4QLbuvmoZW9AVkJimATAYhBHAA7J191ULNSMAhoIwAnhAtu6+arFqRj5o7VBHd2gCgIEaVBhZvXq1JkyYoEgkojlz5mjr1q39XvvAAw/oU5/6lEaPHq3Ro0erurr6jNcDOF227r5qGT0iT6FgQJJ0vJkt4QGkJ+0wsn79etXU1GjFihXavn27pk+frrlz5+rIkSN9Xr9lyxZde+21+v3vf6+6ujpVVlbqyiuv1KFDh4bceGC4yObdVyUpGAykVtQwVQMgTWmHkZUrV+rGG2/U4sWLdfHFF2vNmjUaMWKE1q5d2+f1jz76qG6++WbNmDFDF110kR588EElEglt3rx5yI0Hhots3n3VktyFtbnN5ZYA8Jq0wkgsFtO2bdtUXV2deoFgUNXV1aqrqxvQa7S2tqqjo0Njxozp95r29nZFo9FeX8BwZu2+Oi5LR0akHst7m5imAZCetMLIsWPHFI/HVVpa2uvx0tJS1dfXD+g1vve976mioqJXoPmw2tpaFRcXJ78qKyvTaSbgO14YGRmbHBlhmgZAejK6mubee+/VunXr9OSTTyoS6f+X6rJly9TY2Jj8OnjwYAZbCWSfBg+NjFAzAiBdOelcXFJSolAopIaGhl6PNzQ0qKys7IzP/bd/+zfde++9+t3vfqdLL730jNeGw2GFw9n7SxfINGtkZFwWj4xwci+AwUprZCQvL08zZ87sVXxqFaNWVVX1+7x//dd/1Q9+8ANt2rRJs2bNGnxrgWHIGNNjaW/2hvTkyAjTNADSlNbIiCTV1NRo0aJFmjVrlmbPnq1Vq1appaVFixcvliQtXLhQ48ePV21trSTpX/7lX7R8+XI99thjmjBhQrK2pKCgQAUFBTZ2BfCnbN991ZIqYCWMAEhP2mFkwYIFOnr0qJYvX676+nrNmDFDmzZtSha1HjhwQMFgasDlZz/7mWKxmP72b/+21+usWLFCd95559BaDwwD1hRNtu6+amGaBsBgpR1GJGnp0qVaunRpnz/bsmVLr+/feeedwbwFgG6p4tXsrReRUiMjTe2dauuIK5KbvcEJQHbhbBogy6WKV7N3ikaSiiI5ysvp+pXC6AiAdBBGgCznlZGRQCDAXiMABoUwAmS5o03Zv5LGUsJeIwAGgTACZDlrZCSbd1+1WCMjxxgZAZAGwgiQ5ZI1I1m8rNfCLqwABoMwAmS5ZM2IJ0ZG8iQRRgCkhzACZLGeu68yMgLArwgjQBZrPJXafTXbl/ZKPXZhpWYEQBoII0AWa4h6Y/dVSwlLewEMAmEEyGJHmryxx4il5zSNMcbl1gDwCsIIkMWskREvTNFIqZGRto6EWmJxl1sDwCsII0AW89rIyMhwjkbmdU0nUcQKYKAII0AWs1bSeGH3VQu7sAJIF2EEyGKpkRHvhJHk+TSEEQADRBgBslhDcmTEG9M0klQxKl+SdOBEq8stAeAVhBEgi3lp91XLheMKJElvNzS53BIAXkEYAbKUMcZT59JYLiwtlCTtPkIYATAwhBEgSzWe6lCs0zu7r1oml3aNjOw50qxEgr1GAJwdYQTIUtaoyCiP7L5qOW/MCOWFgmrrSOi9D0653RwAHkAYAbLUex90FYCWeaheRJJyQkFNGjtSkrSbuhEAA0AYAbLUm+93/UU+pazQ5ZakbzJ1IwDSQBgBstQb70clSReXF7nckvRZdSNvNzS73BIAXkAYAbLUm4e7wshUD4aRj4zrGhl5m5ERAANAGAGyUGusU/uPt0jyZhhhRQ2AdBBGgCy0q75JxkhjC8Ma66E9RiznnzNSeTldK2oOfsBOrADOjDACZKE3PVwvIkmhYEAXjO0aHdlN3QiAsyCMAFnoDQ/Xi1iS28JTNwLgLAgjQBZKjoxUeDeMsKIGwEARRoAsk0gY7arvGk24uNx7e4xYkmfUsPEZgLMgjABZ5t0TrWqNxRXJDWpiSYHbzRk0a+OzPUeaFWdFDYAzIIwAWcaqF5lSWqhQMOByawbvvDEjlJcTVHtnIrm1PQD0hTACZBk/1ItIrKgBMHCEESDLWNvAe3kljcUqYqVuBMCZEEaALOP1PUZ66lk3AgD9IYwAWeSDlpjeb2yTJF3kgzBi7TXCyAiAMyGMAFnEGhU5/5wRKgjnuNyaobuQFTUABoAwAmSRZL1ImfdHRaSuFTXh7hU1B0+wogZA3wgjQBaxlvV6fSWNpeeKmrepGwHQD8IIkEX8tJLGwooaAGdDGAGyRHtnPLnqxC8jI1KqbuRtwgiAfhBGgCyx50izOhNGRZEcVRRH3G6ObVIrapimAdA3wgiQJXrWiwQC3t0G/sOsKafdDU063tzucmsAZCPCCJAl3ny/axrDT/UiklQ5ZoQuPbdYnQmj//nTYbebAyALEUaALPHaoUZJ/th59cO+dNl4SdITOw653BIA2YgwAmSBd4+36I/vnpAkzZl4jsutsd+86RXKCQb06nuN2nOEQlYAvRFGgCzwf156V8ZIn548VuedM8Lt5tjunIKwrpgyVpL0xHZGRwD0RhgBXHYqFtcvX3lPkrSw6nyXW+OcL33sXEnSkzsOKcHW8AB6IIwALnvqT4fUeKpDlWPydcWUcW43xzGfvWicCiM5er+xTS/tO+52cwBkEcII4CJjjP6z7l1J0t/NOV+hoH+W9H5YJDekL1xaIUn6DVM1AHogjAAu2n7gpF4/HFU4J6ivzKp0uzmO+5uPda2q2fTa+2qNdbrcGgDZgjACuOi/6t6R1LXaZPTIPHcbkwEzzx+t88aMUEssrmdfb3C7OQCyBGEEcMmx5nZt/HO9JH8XrvYUCAT01917jvxm+3sutwZAtiCMAC5Z/8eDisUTml45SpeeO8rt5mTMl7qnav6w55gaom0utwZANiCMAC7ojCf0f17qKlxdNExGRSznnzNSM88frYSRaje+yTJfAIQRINMSCaMfP/OW3m9s05iRebrqknK3m5Rx/99fXqhQMKANOw/rtg2vyRgCCTCcEUaADIp1JlTzy53638/vkyTd8leTFckNudyqzPv05LFa+ZXpCgSkx7ce0F3/8waBBBjGBhVGVq9erQkTJigSiWjOnDnaunXrGa//1a9+pYsuukiRSESXXHKJNm7cOKjGAl7WeKpDi9Zu1Yadh5UTDOjHf3uprv+L4TVF09P8GeP1r39zqSTpkRff0b2bdhFIgGEq7TCyfv161dTUaMWKFdq+fbumT5+uuXPn6siRI31e/+KLL+raa6/VDTfcoB07duiaa67RNddco9dee23IjQe84uCJVn15zYuq23dcBeEcrf37j+vLw2BfkbP58qxK/fCaaZKk//3/79Nd//OGDp5odblVADItYNL8p8icOXP08Y9/XD/96U8lSYlEQpWVlfrHf/xH3Xrrraddv2DBArW0tOi3v/1t8rG/+Iu/0IwZM7RmzZoBvWc0GlVxcbEaGxtVVOS/49XhP20dcW179wP937eP6f++fVSvH45KksYVhvXw4o/roxXFLrcwuzz0wn794LdvJL+fNr5In59Wrs9NK9PEc0Yq6OOdaQE/G+jf3znpvGgsFtO2bdu0bNmy5GPBYFDV1dWqq6vr8zl1dXWqqanp9djcuXO1YcOGft+nvb1d7e3tye+j0Wg6zRywh17Yr/c+4F9hSOkvmhtjZHpc05kwiicS3f81ao3F1djaoZOnYjrZ2qEPWmPqiPd+sY9PGK1VX71M40flO9sJD7rhkxM1rjCsx14+oJf3H9drh6J67VBUP37mLQUD0piRYZUU5KmkIKyi/BzlBIPKCQWU2/3fYCCgQECyIksgcHp46eMhAD38w+UTVTnGnVPD0wojx44dUzweV2lpaa/HS0tLtWvXrj6fU19f3+f19fX1/b5PbW2t7rrrrnSaNihPv3pY2w+cdPx9MDyNKwzrkxeW6H9dOFaXf6REYwvDbjcpq82bXqF50yt0rLldz73RoI1/fl91e4+rM2F0rLldx5rbJTW53UzAt+ZNr/BGGMmUZcuW9RpNiUajqqy0f379b2aeq6oLzrH9deFtAfX+J7T1L+pAj29yggGFgoHkfyO5IY0ekadRI3JVnJ+r0SPzVFEc6fNf6DizkoKwrp19nq6dfZ464gmdaInpWHO7jjfHdLylXU1tneqIG3XGu0amOuIJJYySw1o9/mdSalwLQH9KiyKuvXdaYaSkpEShUEgNDb3PlGhoaFBZWVmfzykrK0vrekkKh8MKh53/V+R1c4bvSgbAC3JDQZUWRVz9JQnAeWmtpsnLy9PMmTO1efPm5GOJREKbN29WVVVVn8+pqqrqdb0kPffcc/1eDwAAhpe0p2lqamq0aNEizZo1S7Nnz9aqVavU0tKixYsXS5IWLlyo8ePHq7a2VpL0rW99S5/+9Kd133336eqrr9a6dev0yiuv6Oc//7m9PQEAAJ6UdhhZsGCBjh49quXLl6u+vl4zZszQpk2bkkWqBw4cUDCYGnD5xCc+occee0y33367vv/97+vCCy/Uhg0bNG3aNPt6AQAAPCvtfUbcwD4jAAB4z0D//uZsGgAA4CrCCAAAcBVhBAAAuIowAgAAXEUYAQAAriKMAAAAVxFGAACAqwgjAADAVYQRAADgqrS3g3eDtUlsNBp1uSUAAGCgrL+3z7bZuyfCSFNTkySpsrLS5ZYAAIB0NTU1qbi4uN+fe+JsmkQiocOHD6uwsFCBQMC2141Go6qsrNTBgwd9e+aN3/tI/7zP732kf97n9z462T9jjJqamlRRUdHrEN0P88TISDAY1LnnnuvY6xcVFfnyD1hPfu8j/fM+v/eR/nmf3/voVP/ONCJioYAVAAC4ijACAABcNazDSDgc1ooVKxQOh91uimP83kf6531+7yP98z6/9zEb+ueJAlYAAOBfw3pkBAAAuI8wAgAAXEUYAQAAriKMAAAAV3k2jDz//POaN2+eKioqFAgEtGHDhrM+Z/Xq1Zo6dary8/M1ZcoU/ed//mevnz/yyCMKBAK9viKRSK9rjDFavny5ysvLlZ+fr+rqar399tt2dk2Se/37+7//+9Ou+dznPmdn15Kc6KMknTx5UkuWLFF5ebnC4bAmT56sjRs3nvY6EyZMUCQS0Zw5c7R161a7upXkVv/uvPPO0z7Diy66yM6uSXKmf1dcccVpbQ8EArr66quT12TqHpTc62Om7kOn/oyuWrVKU6ZMUX5+viorK3XLLbeora3ttNdx+h6U3Oujl+/Djo4O3X333brgggsUiUQ0ffp0bdq0qc/Xse0zNB61ceNGc9ttt5knnnjCSDJPPvnkGa+///77TWFhoVm3bp3Zu3evefzxx01BQYF56qmnktc8/PDDpqioyLz//vvJr/r6+l6vc++995ri4mKzYcMG86c//cl88YtfNBMnTjSnTp3yRf8WLVpkPve5z/W65sSJE7b2zeJEH9vb282sWbPMVVddZV544QWzf/9+s2XLFrNz587kNevWrTN5eXlm7dq15vXXXzc33nijGTVqlGloaPBF/1asWGE++tGP9voMjx49amvfnOrf8ePHe7X7tddeM6FQyDz88MPJazJ1D7rZx0zdh07079FHHzXhcNg8+uijZv/+/eaZZ54x5eXl5pZbbklek6l70M0+evk+/O53v2sqKirM008/bfbu3Wvuv/9+E4lEzPbt25PX2P0ZejaM9DSQD6Cqqsp85zvf6fVYTU2Nufzyy5PfP/zww6a4uLjf10gkEqasrMz8+Mc/Tj528uRJEw6HzeOPPz6otg9EpvpnTNcvwfnz5w+ypYNnVx9/9rOfmUmTJplYLNbv68yePdssWbIk+X08HjcVFRWmtrZ2cI0fgEz2b8WKFWb69OlDaW7a7Orfh/3kJz8xhYWFprm52Rjj3j1oTOb6aIw796Fd/VuyZIn57Gc/e8Zr3LgHjclsH718H5aXl5uf/vSnva750pe+ZK677rrk93Z/hp6dpklXe3v7aVMS+fn52rp1qzo6OpKPNTc36/zzz1dlZaXmz5+v119/Pfmz/fv3q76+XtXV1cnHiouLNWfOHNXV1TnfiTOwo3+WLVu2aNy4cZoyZYpuuukmHT9+3PH2D8RA+vjUU0+pqqpKS5YsUWlpqaZNm6Z77rlH8XhckhSLxbRt27Zen2EwGFR1dbUnPsOz9c/y9ttvq6KiQpMmTdJ1112nAwcOZKwf/Rnon9GeHnroIX31q1/VyJEjJWX3PSjZ00dLNt6HA+nfJz7xCW3bti05ZL9v3z5t3LhRV111laTsvgcle/po8ep92N81L7zwgiSHPsNBRZgsowGkwWXLlpmysjLzyiuvmEQiYf74xz+a0tJSI8kcPnzYGGPMiy++aH7xi1+YHTt2mC1btpgvfOELpqioyBw8eNAYY8wf/vCHXtdbvvzlL5uvfOUrjvTNmMz1zxhjHn/8cfPf//3f5tVXXzVPPvmkmTp1qvn4xz9uOjs7HeufMfb1ccqUKSYcDpt/+Id/MK+88opZt26dGTNmjLnzzjuNMcYcOnTISDIvvvhir9f+53/+ZzN79mxH+mZM5vpnTNew7S9/+Uvzpz/9yWzatMlUVVWZ8847z0Sj0azvX08vv/yykWRefvnl5GNu3YPGZK6PxrhzH9rZv3//9383ubm5Jicnx0gy3/zmN5M/c+seNCZzfTTG2/fhtddeay6++GKze/duE4/HzbPPPmvy8/NNXl6eMcaZz3DYhJHW1lazePFik5OTY0KhkKmoqDDf/e53jaTT6iYssVjMXHDBBeb22283xmR3GLGjf33Zu3evkWR+97vfDaULZ2VXHy+88EJTWVnZ65f2fffdZ8rKyowx2R1G7OhfXz744ANTVFRkHnzwQVv60hcn/ox+4xvfMJdcckmvx7I9jNjRx75k4j60q3+///3vTWlpqXnggQfMq6++ap544glTWVlp7r77bmNM9ocRO/rYFy/dh0eOHDHz5883wWDQhEIhM3nyZHPzzTebSCRijHHmMxw20zT5+flau3atWltb9c477+jAgQOaMGGCCgsLNXbs2D6fk5ubq8suu0x79uyRJJWVlUmSGhoael3X0NCQ/Jlb7OhfXyZNmqSSkpIzXpMpA+ljeXm5Jk+erFAolHze1KlTVV9fr1gsppKSEoVCIc9+hmfrX19GjRqlyZMnu/4ZpvNntKWlRevWrdMNN9zQ6/Fsvgcle/rYl2y5DwfSvzvuuEPXX3+9vv71r+uSSy7RX//1X+uee+5RbW2tEolEVt+Dkj197IuX7sOxY8dqw4YNamlp0bvvvqtdu3apoKBAkyZNkiRHPsNhE0Ysubm5OvfccxUKhbRu3Tp94QtfUDDY9/8N8Xhcf/7zn1VeXi5JmjhxosrKyrR58+bkNdFoVC+//LKqqqoy0v6zGUr/+vLee+/p+PHjZ7wm087Ux8svv1x79uzp9Qth9+7dKi8vV15envLy8jRz5sxen2EikdDmzZs98RmerX99aW5u1t69e7PmMxzIn9Ff/epXam9v19/93d/1etwL96A0tD72JdvuwzP1r7W19bS+WuHZGOOJe1AaWh/74sX7MBKJaPz48ers7NRvfvMbzZ8/X5Kc+QwHNZ6SBZqamsyOHTvMjh07jCSzcuVKs2PHDvPuu+8aY4y59dZbzfXXX5+8/q233jL/9V//ZXbv3m1efvlls2DBAjNmzBizf//+5DV33XWXeeaZZ8zevXvNtm3bzFe/+lUTiUTM66+/nrzm3nvvNaNGjUrO586fP9+RZYVu9K+pqcl85zvfMXV1dWb//v3md7/7nfnYxz5mLrzwQtPW1mZr/5zq44EDB0xhYaFZunSpeeutt8xvf/tbM27cOPPDH/4wec26detMOBw2jzzyiHnjjTfMN77xDTNq1Kh+p7O81r9/+qd/Mlu2bDH79+83f/jDH0x1dbUpKSkxR44cyfr+WT75yU+aBQsW9Pm+mboH3epjJu9DJ/q3YsUKU1hYaB5//HGzb98+8+yzz5oLLrig1zRapu5BN/vo5fvwpZdeMr/5zW/M3r17zfPPP28++9nPmokTJ5oPPvggeY3dn6Fnw8jvf/97I+m0r0WLFhljupbGffrTn05e/8Ybb5gZM2aY/Px8U1RUZObPn2927drV6zW//e1vm/POO8/k5eWZ0tJSc9VVV/VaV21M19LCO+64w5SWlppwOGz+8i//0rz11lu+6F9ra6u58sorzdixY01ubq45//zzzY033ujILwin+mhMV6HunDlzTDgcNpMmTTI/+tGPTiv8+4//+I/k/xezZ882L730km/6t2DBAlNeXm7y8vLM+PHjzYIFC8yePXs8079du3YZSebZZ5/t830zdQ8a404fM3kfOtG/jo4Oc+edd5oLLrjARCIRU1lZaW6++eZef5EZk5l70M0+evk+3LJli5k6daoJh8PmnHPOMddff705dOjQae9t52cYMKafMSUAAIAMGHY1IwAAILsQRgAAgKsIIwAAwFWEEQAA4CrCCAAAcBVhBAAAuIowAgAAXEUYAQAAriKMAAAAVxFGAACAqwgjAADAVYQRAADgqv8HF9t9MbIdZTQAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "line = modeling.models.Gaussian1D(mean=1.967, stddev=.001)\n",
    "x = np.linspace(1.95, 1.99, 100)\n",
    "plt.plot(x, line(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "bc285af3-5ec9-4ba1-9669-63e06a23334e",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f98ba78de90>]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABndklEQVR4nO3dd3xb5b0/8M/R9t4zceLsQSaBpAFSlkuAsjq4gdJAw2gZubclpSNtIZS2hHtbKL0tNC0Q1i0EaIFC4RcgKWElZCdAdmI7zvJe8pIs6fn9cfQcSbZkS7JkjXzer5dfpbZknxNbR9/zfMejCCEEiIiIiGJEF+sDICIiotMbgxEiIiKKKQYjREREFFMMRoiIiCimGIwQERFRTDEYISIiophiMEJEREQxxWCEiIiIYsoQ6wMIhsvlwsmTJ5GRkQFFUWJ9OERERBQEIQSsVitKS0uh0wVe/0iIYOTkyZMoKyuL9WEQERFRGI4dO4aRI0cG/HpCBCMZGRkA1JPJzMyM8dEQERFRMNrb21FWVqa9jweSEMGITM1kZmYyGCEiIkowg5VYsICViIiIYorBCBEREcUUgxEiIiKKKQYjREREFFMMRoiIiCimGIwQERFRTDEYISIiophiMEJEREQxxWCEiIiIYorBCBEREcUUgxEiIiKKKQYjREREFFMMRohOQ21dvXjs/cM42dod60MhImIwQnQ6em5TNX77zgH8ecORWB8KERGDEaLT0ecn2gAAp9p6YnwkREQMRohOS/trrQCA5k5bjI+EiIjBCNFpp8PmQE1zFwCgqdMe46MhImIwQnTaOVDbrv13cweDESKKPQYjRKeZfaes2n9bbQ7YHM4YHg0REYMRotPOfq+VEQBo6eyN0ZEQEakYjBCdZrxXRgCgiUWsRBRjDEaITiMul8ABdyeN2aC+/JtZxEpEMcZghOg0cqK1Gx02B0x6HWaOzAbAYISIYi/kYOTDDz/ElVdeidLSUiiKgtdff33Q52zYsAFnnnkmzGYzxo8fj2eeeSaMQyWiodp3Sq0XGV+YjoJMMwCgiR01RBRjIQcjnZ2dmDlzJh577LGgHl9VVYWvfvWruPDCC7Fr1y784Ac/wK233op33nkn5IMloqGR9SJTSjKRl2YCwJURIoo9Q6hPuOyyy3DZZZcF/fhVq1ZhzJgxePjhhwEAU6ZMwccff4zf//73WLhwYag/noiGQHbSTCnJQIfNAYCDz4go9qJeM7Jp0yZUVFT4fG7hwoXYtGlTwOfYbDa0t7f7fBDR0Mkx8L4rI+ymIaLYinowUltbi6KiIp/PFRUVob29Hd3d/rcvX7lyJbKysrSPsrKyaB8mUdLrsjtQ3dQJAJhcnIHcNLVmhGkaIoq1uOymWb58Odra2rSPY8eOxfqQiBLegVorhAAKMszISzcj170ywgJWIoq1kGtGQlVcXIy6ujqfz9XV1SEzMxMpKSl+n2M2m2E2m6N9aESnFe8UDQDkpbuDEa6MEFGMRX1lZP78+Vi/fr3P59577z3Mnz8/2j+aiLzItt4pxRkAoK2MtHX3otfpitlxERGFHIx0dHRg165d2LVrFwC1dXfXrl2oqakBoKZYbrzxRu3xt99+OyorK/HjH/8Y+/fvx+OPP46XX34Zd999d2TOgIiCst/d1ju5RA1GclJNUBT1ay1dXB0hotgJORjZtm0bZs+ejdmzZwMAli1bhtmzZ+O+++4DAJw6dUoLTABgzJgxeOutt/Dee+9h5syZePjhh/Hkk0+yrZdoGAkhsE9r61XTNHqdgpxUzhohotgLuWbkggsugBAi4Nf9TVe94IILsHPnzlB/FBFFyInWblh7HDDqFYzNT9c+n5tmQnOnHc0sYiWiGIrLbhoiiiyZohlXkA6TwfOy1zpquDJCRDHEYIToNCAnr051p2gkjoQnonjAYIToNLCvT/GqxJURIooHDEaITgOyeHVycaCVEY6EJ6LYYTBClOTsDheqGz1j4L3lMk1DRHGAwQhRkmvr7oVLAIoC5Kf7TjbOdf//RnbTEFEMMRghSnLtPb0AgHSzATqd4vM1FrASUTxgMEKU5Kw9DgBApsXY72tM0xBRPGAwQpTk2rvVlZEMS/8Zh3KzvJYuO5yuwMMMiYiiicEIUZLTVkZS+q+MyHHwQgCt3J+GiGKEwQhRkpM1I5l+VkaMeh2y3EEKUzVEFCsMRoiSnFULRvqvjACeIlYOPiOiWGEwQpTk2rsDp2kAFrESUewxGCFKcjJN46+AFeBIeCKKPQYjREluoNZewNNR08zBZ0QUIwxGiJLcQK29gHeahvvTEFFsMBghSnIDtfYCQG6aOhKeaRoiihUGI0RJbrCaEa2bhmkaIooRBiNESW6wmhF20xBRrDEYIUpysmZksNZepmmIKFYYjBAlMadLwGpTV0YCpWny09WakZYuO1zcn4aIYoDBCFES63AHIkDgYCQnTV0xcbqEVl9CRDScGIwQJTGZojEbdDAb9H4fYzbokWFWAxWmaogoFhiMECWxwdp6pdx0FrESUewwGCFKYoO19Uq5bO8lohhiMEKUxAZr65Xy2N5LRDHEYIQoiQ3W1itxJDwRxRKDEaIkFnyahiPhiSh2GIwQJbFQ0zSsGSGiWGAwQpTEtDRNkAWsrBkholhgMEKUxEJt7WWahohigcEIURILtmYk310zwgJWIooFBiNESSzYmhHvoWdCcH8aIhpeDEaIkphcGclMGXhlRBaw9jo9G+sREQ0XBiNESUwWsGYMsjJiMeqRalL3rmm0MlVDRMOLwQhREgs2TQMAo3JTAQDVTZ1RPSYior4YjBAlKSFE0AWsADC+MB0AcKiuI6rHRUTUF4MRoiRlc7jQ61SLUQdr7QWACYUZAIBD9QxGiGh4MRghSlKyXkSnAGnuepCBTCxyr4wwGCGiYcZghChJtbvrRTIsRiiKMujjJ7iDkcN1Vrb3EtGwYjBClKSCbeuVRuelwaBT0Gl34mRbTzQPjYjIB4MRoiSltfWaB68XAQCjXocx+WkAgEN11qgdFxFRXwxGiJKUZ1+a4FZGAK9UDetGiGgYMRghSlKett7gVkYAYLzsqGF7LxENIwYjREkqlIFnkqejhmkaIho+DEaIkpRnFHwIaRqvlRF21BDRcAkrGHnsscdQXl4Oi8WCefPmYcuWLQM+/tFHH8WkSZOQkpKCsrIy3H333ejpYbU+UTR5akaCXxkpz0+FXqfAanOgrp171BDR8Ag5GHnppZewbNkyrFixAjt27MDMmTOxcOFC1NfX+338Cy+8gJ/+9KdYsWIF9u3bh6eeegovvfQSfvaznw354IkoMK21N4SVEbNBj9F56h41TNUQ0XAJORh55JFHcNttt2HJkiWYOnUqVq1ahdTUVKxevdrv4zdu3Ihzzz0X3/rWt1BeXo5LLrkE119//aCrKUQ0NDJNE0rNCABM4B41RDTMQgpG7HY7tm/fjoqKCs830OlQUVGBTZs2+X3OOeecg+3bt2vBR2VlJd5++21cfvnlAX+OzWZDe3u7zwcRhSac1l6Ae9QQ0fAL6SrV2NgIp9OJoqIin88XFRVh//79fp/zrW99C42NjTjvvPMghIDD4cDtt98+YJpm5cqV+OUvfxnKoRFRH+G09gLes0aYpiGi4RH1bpoNGzbgwQcfxOOPP44dO3bg1VdfxVtvvYVf/epXAZ+zfPlytLW1aR/Hjh2L9mESJZ1wWnsBz8rIQXbUENEwCWllJD8/H3q9HnV1dT6fr6urQ3Fxsd/n3HvvvVi8eDFuvfVWAMD06dPR2dmJ7373u/j5z38Ona5/PGQ2m2E2m0M5NCLqI5zWXgAYW5AGnQK0dfeiocOGwgxLNA6PiEgT0sqIyWTCnDlzsH79eu1zLpcL69evx/z58/0+p6urq1/Aoder25nzrosoOhxOFzrtTgChtfYCgMWox6hctaPmMItYiWgYhJymWbZsGZ544gk8++yz2LdvH+644w50dnZiyZIlAIAbb7wRy5cv1x5/5ZVX4s9//jPWrFmDqqoqvPfee7j33ntx5ZVXakEJEUVWh82h/XeoKyOA11h4FrES0TAI+Sq1aNEiNDQ04L777kNtbS1mzZqFtWvXakWtNTU1Pishv/jFL6AoCn7xi1/gxIkTKCgowJVXXonf/OY3kTsLIvLR3q0GIylGPYz60EvDJhSlY92+Os4aIaJhoYgEyJW0t7cjKysLbW1tyMzMjPXhEMW9L0604Yo/foyiTDM2/6xi8Cf08eqO41j28m7MG5OLl77nPwVLRDSYYN+/uTcNURIKt61XmljENA0RDR8GI0RJyNPWG3q9CACMK0iHogDNnXY0dXCPGiKKLgYjREnI09Yb3spIikmPkTkpALg6QkTRx2CEKAm1h7Fjb18cC09Ew4XBCFESsoaxY29fcsO8w3XsqCGi6GIwQpSEZGtvuGkaQJ3ECgDVTV0ROSYiokAYjBAlIW1lJMQde71luVM8nV4D1IiIooHBCFESGmprLwCkmtRApoPBCBFFGYMRoiQ01NZeAEgzq8/ttDMYIaLoYjBClITatQLW8FdG0mUwYnNG5JiIiAJhMEKUhGQB61BqRtLM6kaWrBkhomhjMEKUhKwRqBlJc9eM2BwuOJyuiBwXEZE/DEaIkowQwjP0bCjBiNmzqsJUDRFFE4MRoiTT3euE06Vuxj2UNI3JoINJr14iOljESkRRxGCEKMnIehG9TkGKUT+k75XqrhvpYt0IEUURgxGiJOM9Cl5RlCF9rzTOGiGiYcBghCjJRGLgmcT2XiIaDgxGiJJMJNp6Ja29lzUjRBRFDEaIkoy2MmIe+sqINoWVaRoiiiIGI0RJRmvrjcTKiInBCBFFH4MRoiRjjcAoeEmujHSwZoSIoojBCFGSkTUjkSlgdbf2smaEiKKIwQhRkvF00ww9TZNqZmsvEUUfgxGiJNPaZQcA5KaZhvy90lnASkTDgMEIUZJp7lSDkZwIBCNpJrlzL2tGiCh6GIwQJZmWTjVNk5sagWBEroywZoSIoojBCFGSae6SKyOcM0JEiYHBCFESEUKgpTNyNSNs7SWi4cBghCiJWG0OOFwCAJATgTSNbO3lyggRRRODEaIk0uquF0kx6mEx6of8/eTKCOeMEFE0MRghSiLNEWzrBTzj4DlnhIiiicEIURJp6Yxc8SrgWRnp6XXB4XRF5HsSEfXFYIQoiWgzRiJQLwIAaWZPqqfTziJWIooOBiNESaQlwmkas0EPo14BwLoRIooeBiNESSTSKyMAkGrirBEiii4GI0RJpKVL7aaJZDCSzlkjRBRlDEaIkohn4FlkClgBT90IV0aIKFoYjBAlEc8oeKZpiChxMBghSiLaykgU0jTcLI+IooXBCFESaYnCyohM07BmhIiihcEIUZIQQmgFrJFq7QW8RsIzTUNEUcJghChJtPc44HRvkpedGsECVtaMEFGUMRghShKyXiTNpIfZMPRN8qQ0tvYSUZQxGCFKEtHopAGAdLb2ElGUMRghShKeGSORDUbS2E1DRFHGYIQoSURjFDzAmhEiir6wgpHHHnsM5eXlsFgsmDdvHrZs2TLg41tbW3HXXXehpKQEZrMZEydOxNtvvx3WARORf61R6KQBvFZGWDNCRFFiCPUJL730EpYtW4ZVq1Zh3rx5ePTRR7Fw4UIcOHAAhYWF/R5vt9vxla98BYWFhfj73/+OESNG4OjRo8jOzo7E8RORm6wZiWQnDeA9Z4QrI0QUHSEHI4888ghuu+02LFmyBACwatUqvPXWW1i9ejV++tOf9nv86tWr0dzcjI0bN8JoVC+S5eXlQztqIuonGtNXAc8E1i7WjBBRlISUprHb7di+fTsqKio830CnQ0VFBTZt2uT3OW+88Qbmz5+Pu+66C0VFRZg2bRoefPBBOJ2Bl3xtNhva29t9PohoYFrNSITTNHJvGrb2ElG0hBSMNDY2wul0oqioyOfzRUVFqK2t9fucyspK/P3vf4fT6cTbb7+Ne++9Fw8//DB+/etfB/w5K1euRFZWlvZRVlYWymESnZbkKPhI14xoe9MwTUNEURL1bhqXy4XCwkL89a9/xZw5c7Bo0SL8/Oc/x6pVqwI+Z/ny5Whra9M+jh07Fu3DJEp4UeumcdeMdPc6tQmvRESRFFLNSH5+PvR6Perq6nw+X1dXh+LiYr/PKSkpgdFohF7vmQg5ZcoU1NbWwm63w2Tqf+E0m80wm82hHBrRaS/a3TSAWjeSYYlsgSwRUUgrIyaTCXPmzMH69eu1z7lcLqxfvx7z58/3+5xzzz0Xhw8fhsvl0j538OBBlJSU+A1EiCh0Lpfw7Ngb4W4as0EHvU4BwPZeIoqOkNM0y5YtwxNPPIFnn30W+/btwx133IHOzk6tu+bGG2/E8uXLtcffcccdaG5uxve//30cPHgQb731Fh588EHcddddkTsLotNce08vZAYlO8JpGkVRkGZiey8RRU/Irb2LFi1CQ0MD7rvvPtTW1mLWrFlYu3atVtRaU1MDnc4T45SVleGdd97B3XffjRkzZmDEiBH4/ve/j5/85CeROwui05ysF8kwG2AyRL4ULN1sQHuPg0WsRBQVIQcjALB06VIsXbrU79c2bNjQ73Pz58/Hp59+Gs6PIqIgtERpkzyJ+9MQUTRxbxqiJNDcqRavRisYSeVIeCKKIgYjRElAmzES4eJVKd3d3ss0DRFFA4MRoiTQEqUZI1KaNoWVwQgRRR6DEaIk0DxMNSPcn4aIooHBCFES0DbJi1owIlt7WTNCRJHHYIQoCWgFrNFK03B/GiKKIgYjREnAs0lelApYTUzTEFH0MBghSgKeUfDRbe1lmoaIooHBCFES0LppolQzwtZeIoomBiNECc7pEmjtHp6aEbb2ElE0MBghSnBt3b0Q2iZ50akZYWsvEUUTgxGiBCc3ycu0GGDUR+clLYeecRw8EUUDgxGiBOfppIlOigbwnjPClREiijwGI0QJLtrFqwCQzjkjRBRFDEaIEly023oB75oRJ1wuEbWfQ0SnJwYjRAku2tNXAU/NCAB09bJuhIgii8EIUYKL9vRVALAYddAp6n8zVUNEkcZghCjBNQ9DzYiiKJw1QkRRw2CEKMFpO/ZGMU0DeIpYu9jeS0QRxmCEKMFpBaxRXBkBgFQT23uJKDoYjBAluJau6BewAmzvJaLoYTBClOBkzUg0C1gBT3tvJ0fCE1GEMRghSmAOpwttUd4kT9KCEdaMEFGEMRghSmByt15FAbJSorwy4q4ZYZqGiCKNwQhRAmt114tkmA0wRGmTPImtvUQULQxGiBKYtccdjFiiuyoCsICViKKHwQhRApOrFBkWwyCPHDpPAStrRogoshiMECWwjp7hC0ZSWTNCRFHCYIQogVndgYFMoUQT0zREFC0MRogSmNW9MpI+DDUjnDNCRNHCYIQogQ1nmibNLNM0rBkhoshiMEKUwDpsntbeaEszMU1DRNHBYIQogXUMY80I54wQUbQwGCFKYO1azcjwFbB2sbWXiCKMwQhRAvPUjES/gDVV1ozYHRBCRP3nEdHpg8EIUQIbzjSN/BlCcHWEiCKLwQhRAhvObpoUox46Rf1vFrESUSQxGCFKYHJvmuFYGVEUxdNRw5URIoogBiNECcw6jHvTAF51I1wZIaIIYjBClKCEEJ6akWEKRtjeS0TRwGCEKEF12Z2QTS0Z5uh30wBAfroZAFDb1jMsP4+ITg8MRogSlNyXRq9TYDEOz0t5bH4aAKCysXNYfh4RnR4YjBAlKG0UvMUARVGG5WeOkcFIQ8ew/DwiOj0wGCFKUNqOvcPQSSONLUgHAFRxZYSIIojBCFGCGs6BZ5JcGalq7OQUViKKGAYjRAnKOowDz6RRuanQ6xR02Z2oa7cN288louQWVjDy2GOPoby8HBaLBfPmzcOWLVuCet6aNWugKAquueaacH4sEXkZzn1pJJNBh7KcFABAZSPrRogoMkIORl566SUsW7YMK1aswI4dOzBz5kwsXLgQ9fX1Az6vuroa99xzDxYsWBD2wRKRhzUGaRrAN1VDRBQJIQcjjzzyCG677TYsWbIEU6dOxapVq5CamorVq1cHfI7T6cQNN9yAX/7ylxg7duyQDpiIVHJlZLgGnkmyiLWygcEIEUVGSMGI3W7H9u3bUVFR4fkGOh0qKiqwadOmgM974IEHUFhYiFtuuSWon2Oz2dDe3u7zQUS+5L40GVwZIaIEF1Iw0tjYCKfTiaKiIp/PFxUVoba21u9zPv74Yzz11FN44okngv45K1euRFZWlvZRVlYWymESnRY6hnlfGmksgxEiirCodtNYrVYsXrwYTzzxBPLz84N+3vLly9HW1qZ9HDt2LIpHSZSYYlUzItM0Nc1dsDtcw/qziSg5hXQVy8/Ph16vR11dnc/n6+rqUFxc3O/xR44cQXV1Na688krtcy6XevEyGAw4cOAAxo0b1+95ZrMZZrM5lEMjOu14akaGr5sGAIoyzUgx6tHd68Sxli6McwcnREThCmllxGQyYc6cOVi/fr32OZfLhfXr12P+/Pn9Hj958mR8/vnn2LVrl/Zx1VVX4cILL8SuXbuYfiEaAlkzMtwrI4qieOpGWMRKRBEQ8lVs2bJluOmmm3DWWWdh7ty5ePTRR9HZ2YklS5YAAG688UaMGDECK1euhMViwbRp03yen52dDQD9Pk9EoZE1I5nDXDMCAGML0rD3VLt71kjRoI8nIhpIyFexRYsWoaGhAffddx9qa2sxa9YsrF27VitqrampgU7Hwa5E0Rar1l6ARaxEFFlhXcWWLl2KpUuX+v3ahg0bBnzuM888E86PJKI+YrFRnsRZI0QUSVzCIEpALpdAhz12KyOyZqSSKyNEFAEMRogSUFevE3LT3Mxh7qYBgDEFajDSYLVphbREROFiMEKUgGS9iEGnwGwY/pdxpsWI/HS1/Z51I0Q0VAxGiBKQ1tZrMUBRlJgcA4tYiShSGIwQJaBYTV/1ptWNsIiViIaIwQhRAuqIYSeNNLaARaxEFBkMRogSkGfg2fAXr0qe3Xs7YnYMRJQcGIwQJSDvmpFYkbNGqho6IWRrDxFRGBiMECWgWA48k0blpkKnAJ12J+qttpgdBxElPgYjRAlIpmliuTJiMuhQlpsKgEWsRDQ0DEaIEpAsYM2IYTACsL2XiCKDwQhRApJpmowYpmkAYEy+3KOGRaxEFD4GI0QJqCMO5owAnrHwXBkhoqFgMEKUgLShZzFs7QWAMXlqMHK0uSumx0FEiY3BCFEC6nC39sa6ZqQgQ92fpqmD3TREFD4GI0QJKF5qRvLSTQCAlq5eOJyumB4LESUuBiNECSgeWnsBICfVBJ17n77mLntMj4WIEheDEaIEFA970wCAXqcgN01dHWnqYDBCROFhMEKUYFwugQ67nDMS2wJWAMhLk3UjDEaIKDwMRogSTKfdAbkVTKwLWAFP3Ugji1iJKEwMRogSjKwXMegUmA2xfwnnpasrIwxGiChcsb+SEVFItHoRiwGKosT4aIA8WTPSyTQNEYWHwQhRgpEDz+IhRQMA+emygJUrI0QUHgYjRAnGqnXSxL54FQDy01nASkRDw2CEKMF0xMnAM0mrGWGahojCxGCEKMF02NRR8LEeeCblMU1DREPEYIQowWij4OMkGMlP83TTCNlzTEQUAgYjRAnGGifTVyW5MtLT60KX3RnjoyGiRMRghCjBxMu+NFKqSQ+LUb2UsIiViMLBYIQowcRbAauiKFpHTWMn60aIKHQMRogSjNVdwBoP+9JIeWzvJaIhYDBClGDirWYEAPLT2FFDROFjMEKUYOKtZgTgZnlENDQMRogSTLzVjADem+UxTUNEoWMwQpRgPHNG4qhmhJvlEdEQMBghSjDxmKYpyJAFrEzTEFHoGIwQJRCXS3iCkXhK06Sxm4aIwsdghCiBdNod2n/Hyzh4wGt/Gs4ZIaIwMBghSiCyXsSoV2A2xM/LVwYjzZ12OF3cn4aIQhM/VzMiGpR3ikZRlBgfjUduqhqMuATQ0sVUDRGFhsEIUQLRBp7FUYoGAAx6HXJS1e4e1o0QUagYjBAlEM/KSPy09UqekfCsGyGi0DAYIUog1h65L018rYwAQL6cwspZI0QUIgYjRAkkHqevSlwZIaJwMRghSiDxOPBM8myWx5URIgoNgxGiBBKPO/ZKnv1puDJCRKEJKxh57LHHUF5eDovFgnnz5mHLli0BH/vEE09gwYIFyMnJQU5ODioqKgZ8PBEFFo/70kienXu5MkJEoQk5GHnppZewbNkyrFixAjt27MDMmTOxcOFC1NfX+338hg0bcP311+P999/Hpk2bUFZWhksuuQQnTpwY8sETnW46bPFbwKqNhOcUViIKUcjByCOPPILbbrsNS5YswdSpU7Fq1SqkpqZi9erVfh//t7/9DXfeeSdmzZqFyZMn48knn4TL5cL69euHfPBEp5t43JdGKshgzQgRhSekYMRut2P79u2oqKjwfAOdDhUVFdi0aVNQ36Orqwu9vb3Izc0N+BibzYb29nafDyICGqzqqkOOu1g0nng2y+PKCBGFJqRgpLGxEU6nE0VFRT6fLyoqQm1tbVDf4yc/+QlKS0t9Apq+Vq5ciaysLO2jrKwslMMkSlo1zV0AgFG5qTE+kv5kzUin3YluuzPGR0NEiWRYu2keeughrFmzBq+99hosFkvAxy1fvhxtbW3ax7Fjx4bxKIniU0+vE3Xt6qrD6DgMRtLNBpjcm/exboSIQhFS4jk/Px96vR51dXU+n6+rq0NxcfGAz/3d736Hhx56COvWrcOMGTMGfKzZbIbZbA7l0IiS3jH3qkiG2YDs1PjrplEUBflpJpxs60Fjhx0jc+IvYCKi+BTSyojJZMKcOXN8ik9lMer8+fMDPu9//ud/8Ktf/Qpr167FWWedFf7REp3GZIqmLDc1rnbs9cYprEQUjpBL8pctW4abbroJZ511FubOnYtHH30UnZ2dWLJkCQDgxhtvxIgRI7By5UoAwH//93/jvvvuwwsvvIDy8nKttiQ9PR3p6ekRPBWi5BbP9SKS3J+GHTVEFIqQg5FFixahoaEB9913H2prazFr1iysXbtWK2qtqamBTudZcPnzn/8Mu92Ob37zmz7fZ8WKFbj//vuHdvREpxEZjIzOi99gRJvCypoRIgpBWMMKli5diqVLl/r92oYNG3z+f3V1dTg/goj6qGnypGniVR5XRogoDNybhihBJESahrNGiCgMDEaIEoAQIiGCEe5PQ0ThYDBClAAarDbYHC7oFGBETkqsDycg7txLROFgMEKUAI66V0VKs1Ng1MfvyzbPPaa+qZMrI0QUvPi9qhGRRhavxnOKBgAKMtSVkeZOO1wuEeOjIaJEwWCEKAEkQr0IAOSkqisjTpdAW3dvjI+GiBIFgxGiBCBHwY+K4xkjAGAy6JCVoo6q5/40RBQsBiNECSBRVkYAT0dNg5V1I0QUHAYjRAngaAIFI3LWCDtqiChYDEaI4ly33YkGq/rGPjo3LcZHM7j8DDlrhMEIEQWHwQhRnDvWoq6KZFoMyEo1xvhoBleYYQEA1FsZjBBRcBiMEMU5ra03zotXJdneW9/OYISIgsNghCjOJVK9CAAUymDE2hPjIyGiRMFghCjOaW29CVAvAgCFmWqapoFpGiIKEoMRojiXSG29gPfKCIOR4SaEwO/fO4i/fngk1odCFBJDrA+AiAaWqMFIc6cdvU5XXO+lk2ze/OwU/rD+EADgm3PKkOveK4go3vEqQRTHXC7hlaZJjGAkJ9UEg04BwPbe4dRpc+DBt/Zp/3/fqfYYHg1RaBiMEMWxeqsNNocLep2C0mxLrA8nKDqdgvx0dtQMt8c3HEZtu6doeO9JBiOUOBiMEMUxmaIZkZ0CQwKlOwozh1430tPrxJu7T3LDvSBUN3biiQ+rAACzR2UD4MoIJZbEuboRnYYSrV5EikR77/99ehT/+eJOPL7hcKQOK2n9+q29sDtdWDAhH3ecPw4AsJfBCCUQFrASxTEZjJQlWDBSIKewDiFNI9MMcugb+ff+gXqs21cPg07BiivPgMWo3mMeaeiA3eGCycB7Top//CslimM1TZ0AgNEJMn1VikR7b5X73Js6uftvIHaHCw+8uRcAsOTccowvTMeI7BRkWgzodQocru+I8RESBYfBSJx76uMqfPV/P2JXwmkqYdM07pqRhiGkaaob1WCkmcFIQP/36VFUNXaiIMOM/7p4AgBAURRMLskEwFQNJQ4GI330Ol1wukSsD0Pz8tZj2HOyHR8dagj4mCMNHfjj+kPotDmG8choONQ0dwNIwGBkiJvltXbZ0dKlFq4yGAlsw0H1unD7+eOQYfFsojjVHYywiJUSBYMRLw1WG+b86j0seWYrXHESkDS4V0SODpA3/907B/Dwewfxjx3Hh+uwaBh02R3ailii1YwUDnGzvCr3qgigBibxdIMQT2RdjeygkRiMUKJhMOJlW3Uz2nsc+PBgA17beSLWh4NepwstXepdoVyu90fmhb840TYsx0XD45h7VSQ71YisFOMgj44vMk3T2GELK7D3DkZcAmzv9aPBakNjhw2KAkwuzvD52hSvYEQIBnKx9PPXPseFv9uA1i6u8A2EwYiXSq8L4ENr98PaE9sLYHOnHfI6cixAMOJyCS1Q2V9rHa5Do2Fw1F3AmWgpGgDITzdDUQCHS6A5jItwtddrEQCaO1kz1Zdc9RiTl4ZUk29j5ISidOh1Clq6en0GodHwauqw4cUtNahq7MSnlU2xPpy4xmDEy5EGT+V5g9WGP/07tvMNvHc9DbQyUmftgc3hAgAcqLVyOTuJJGrxKgAY9Trkpqr7ooSTqqnsE4w0dfCusi9ZnDqlNLPf1yxGPcbmq7s8M1UTO+v21UFekg/WsbNpIAxGvFQ2qBfAb84ZCQBY/UmVT4ASih01LUOOhBu8Omjq2m3o6XX2e4x3LYnN4UJ1U2e/x1Bikr/bRGvrlQqGMPis798xi1j7k0GGrA/pa2qpTNVwxTRW1n5Rq/33wTr+HgbCYMRNCKEFHrctGIsLJxWg1ynwwJt7Q8652h0u3PTUFix+arPP6kao+j73eEv/1ZGjfS7a+0/DC09bVy9OtXXH+jAi7qh7ZWR0blqMjyQ8BWHOGhFCoLpRPXd5dx9OqifZyeLVKSUZfr8+he29MdXe04tPDntuSA9xZWRADEbcGjvssPY4oCjqnei9V0yFUa/gg4MN+Pf++pC+1/GWLlhtDvQ6BXbUtAzhmHwv4v46avp+bn/t6XfhWbx6My5++IMhBX7Rdri+A5tDXCmTA88SrZNGku29of5eGjps6LA5oFOAmWXZAIBmpml89PQ6tVTW1JIsv4+ZkuAdNd12Z8zr9obi/f31sDtdyElVi88rGzvgcLpifFTxi8GIW6V7VWRkToqaby1Ix83njQEAPPCvvbA5+qdIAvFeYt51rDXsY+p7EfdXNyKDkbLcFACn35Ls8ZYufHa8DV12Z9zeAQohcONTm3HdE59id5B/Dw6nC8db1NWeRE3TeAafhRaMyFWRETkpKM5SAxpOYfV1qK4DTpdATqoRRe5/577kiklVYye67Ik1g0gIgW+u2ogLf/cB2roSMyB5Z4+aorlu7iikmvTodQpUc2uDgBiMuMm7jLH56drn/vOiCSjIMONoUxf+tftU0N9LXkwBYOeQVkbUC7DZvbeEv2BEBj6XnlEM4PRbGdlc2az9d6COo1g70tCJk209EAJ4dlN1UM851dYDh0vAZNChONMS3QOMknA3y6tqVG8MyvPStCJY1oz42ntKbeOfWpoJRVH8PqYww4L8dBOEUIvbE8mh+g7sOdmOxg4bPjnSGOvDCVlPrxMbDqgD6S6bVowJher7yiHWjQTEYMTtiHtWx7gCTzCSbjbgihklAEJ7k/eu4/jseFvYHS5ylLZcqu77ZiuE0DYRW+gORo63dKM9gZc2Q7WlyisY8VNTEw+2VXuO8V+7TwU12l9b8cpJgU7n/80m3hWGuVlelVe9SG4agxF/5ArolGL/xauSJ1WTWG+CGw97ApCNCRiMfHSoEV12J0qzLJg+IgsTitRVqqF01LR22fGP7cfRm6SpHgYjbtrKSIFvsaAsoKtqDP6NznsprsvuDLuKWi5vnzkqB0D/lZHmTjusNrXOZdqILJS4l7QPJshd0Km27iG/yWyu8tRhHG+OzyLWrdWe1TG704WXth4b9DlHm+UGeYlZvAp40jShFrDKGSPl+WnITVeDEaZpfGltvQE6aaREncT6yRHP63rjkcSbzyG7aBZOK4aiKNrKyMH68K/N97yyGz98ZTee3VgdiUOMOwxG3GTNSN9gpFwLRoKPaGXqJM2kBxB+3YhM08wZ7QlGvDt7ZLdFcaYFFqNem8K4LwGCkaYOGxb+/kMs+sumsCdE1rX3+AR+8boysv2oujJy1cxSAOrmZoMVsskVr0ScMSJ5p2lC+R1XeQUjee6VkZYAwUh7Ty9++eaeIaVDE40QwtPW62fGiLd47aj5+FBjwIJuh9PlMxahsqETtW2JM7it1+nCun11ADzp84nulZFw0zSVDR1Yt09tpHh3b10EjjL+MBgBYHM4ccxdLOidpgGAMe5gpKa5K6h0S69X4eHCaeofYjgXSpvDqY3AnlmWBZ0C9PS6fIoBj/bZXn5yAt0FbXWP3j9U3xH2Zmqb3SkaGfTFY81IvVUNmBQFuPeKqchLM+FUW492sQok0WeMAJ40TU+vC9YgN3F0uYQWzPdN0/gLaN7+7BSe/qQav33nQISOOv4db+mGtccBo17pd73qSwYj+0+1x81+W4frrVi8ejMWP7XFb8pyz8l2WHscyLAYcIY72EqkVM3myma0dfciL82Es8pzAagTcQE10A4nzfKM12rI9qMtAYt6jzV34Ycv7+43wTgRMBiBehfqdAmkmw3a3ZxUmpUCk0GHXqfAiZbB0wAnWrrhdAlYjDqtjiOclRE5cdKoV1CQbkZJltot452qkW9Y5e6lfLkysj8BgpEdNa3af4cbPG1xp2gun67W9bR09aIjznYu3u5O0UwqykBBhhnXzS0DADy78eiAz9NmjCRwMJJi0iPDrI4pD7Zu5FS7OlHYoFMwIjtFC0bsTpff361cDdtfaz1t9mCRqxwTCjNgMgx8CR9bkAaTXodOu1P7t2rtsuPVHcexLkZ32H/5oBJCqL/Ttz7r3xggC1bnjcnDggkFABIrVbN2j3pOl5xRBL273mtEdgrS3B01fWdDDaatqxevbFM3QU0z6eF0CXx02P8u7v+9dj/+seM4/hjj6eHhYDACtdsBUF+4fSvTdToF5e43hKog/ojkXV15XppW63GoviPkfnm5AqLu8aFoy/X+gpFR7uOT+eEDtda4uQsKZMdRz2pRuJX+spOmYmqRtpGcv8FwsSTrRc523yHdMG80dAqwqbIp4HmrhclyX5rErRkBgILM0Dpq5B3dqLxUGPQ6pJoMsBjVy5S/+qJTrT3a1xpPk1kk+4KsFwHUsfwTi9W78sfeP4xvP7kZc369Dste3o1bn9umbbI5XE61deP1XZ5NSL3/W9roHhR27vg8nDMuDwCw6UhTQgSbLpfAu3vUIE/ejAKAoigYL+tG/BSx1lt78P7+er/X7Ze21aC714nJxRm44UujAQDv7+8fjHTYHNqK67ajzf2+Hu8YjEAdRgN4ilX7kqmaqiBGw8uL6ei8VBRkmDEyJwVCqF01ofAORuT3A/oGI57ARx6nvAs6HsQqTl+tXXb88s09YY/AD5bd4cJnXjsMhxOMNHXYcMh9IT27PFebs3JsCEWsf/3wCOb+Zh02RfAuTF4UzipXA9PS7BRcMlW9SD0XoM23qdOOTrsTiuKZH5Oo5EpjsLNGZL3IGK/C3bw09Xv4K2I90er5fSda+2q4PMGI/8mrfcmOm5e3HcfHhxvhdLeMA8CGA4EHOkbjzf/pT6rR6xSYUpIJnQLsrGn1WSmwOZzY6u4+O3d8Ps4qz4FRr+BEa7ffoY/x5tOqJtRbbcgwG3DOuHyfr3k6avr/nX7/xV1Y8sxWPLR2v8/nHU6Xtop687ljcOGkQgDABwf7By7v7qlFT6+aAjra1IX6BNsgkcEIPHvSjA2Qf/UUsQazMuKbOpnlbssNtW5E5lLlSO2ygVZG3F8z6HVabnJfGPNG/vphJZ7+pBr3/fOLkJ8bir2n2mF3ePKm4ew2LC9Yk4oykJtmQlmO+m8Qbt1IbVsPHn73IOqtNtz1wg6cbPUf1HTaHPj3/jq/+wT11WV3YI97ZLfMHQPAjeeodzev7jih1QV5k7/XkkwLzAZ9yOcSTwpCbO/1Ll6VcgcoYj3ptQ1AvM7YEUJg/b66gEW4odobZPGqdPWsEcgwGzB7VDaWXzYZG+65AD+6ZBIA4MNDgWsxfvT3zzDnV+9F7E2trbsXL2yuAQD8eOEknDtefbP+566T2mN2HG2FzeFCfroZEwrTkWoyYLZ7hTneUzUHaq1Y+sJOAGq9YN8U2sQiOWvE92avtq0Hm9wFu3/9sFL7NwLUYtUTrd3ITTPhqlmlOKs8B+lmAxo77Pj8hO8Nrve/IwBsO5pYRd0MRhC4k0bS2nuDiMw9RaXqc+QLaadXjUQw5J1kgXtlRAYc8s3W2tOr3Sl61xVMLpYFa6G/wX/kvjBtOtIU1uZmwZKBmVy2PFzfEXJR16fuFM28seqbvAzWwu2o+cP6Q9rux82ddtzxtx39pu7WW3vwjT9vxM3PbMOLW2r8fRsfu2pa4XQJlGZZMCLbs8Ixf2weJhalo7vXib9vP97veTXNnlRFogt18JlcWRzjJxjpuzLicgmfLot43Yjs1R0ncMuz27D81c+H/L3ae3q11b9AG+T1dd6EfHz+y4V47c5z8b3zx6E8Pw1fnqjWYmyubPIbWDd22PDqjuNo6rQPGLCE4v8+PYoOmwOTijJwwaQCXD1rBAA1VSNXYTa560XOGZenpcxlqiaei1gP1lnxrSc+RXOnHdNGZOLer07t9xi5MnKoT3vv//tCrTFJdRfi3/vPL/DRITUNs/rjKgDAt+eNgsWoh1Gvw4IJahDnvU1JY4cNH7tns8ivb61OrFTNaR+MqBvkqRfAQJXpY/JlJXQQaRq5MpKvvpHIlZFdx1pDWvaUKyP5GeqFWAYj8q5Z/m9emgkZFqP2PLl0G+pdYkunHV+cVCNtl4DfwrJIkcWrV80sRapJD7vTFXL1txx2NneMOxjJCT9NU9nQgZe3qbM/HvmPmchKMWL3sVY88OZe7THHmrtw7apN2ipOMLl2WS/ivSoCqPnjxe7c7xt+cuZaJ02C14sA3sFIiGmafO80jf/BZ40dNvQ6Pa+peE3TvP25+lr694H6IY9llzcZpVkWZLun04ZjYlE6ijMtsDlcft+03t1TB5kF+OJEaClmf3p6nXj6k2oAwPfOHwtFUbDwjCKYDTpUNnTiixPq9UrOFzl3fJ72XJnuCLdu5K3PTuE7T2+J2t5Vh9yBSFOnHWeUZuL/bpmHrFRjv8fJ9t6+HTXyWrvsKxPxtdkj4HQJ3Pl/O/D37cex7WgLjHoF33ZfLwDgwslqqsY7xfb256fgdAnMHJmFa89Si+S3VXNlJKE0d9rR1t0LRfG9AHqTgcWJlu4B96hxOF3ayoVM05xRmgmjXkFTpz2kN8qGDv8rI/VWG7rtzoCtn9rKSIgX5k+ONML7df7G7pOBHzxEsnh1zugc7QUayvG2dfVqaSgZjIx0//uEU8D6yHsH4XQJXDS5EF8/cyQevW4WFAX42+YavLLtGA7UWvGNP2/E0aYuyGGowby5ynqRs931It4qphYBAD4/0davTa+mT2FyItMGnwWRpnE4XVoa0vu1mBMgGJH1IrJj4WBdR9wVbnfbndodq93hwsdDXGUIpXh1IIqiaHfQHx7sXwwp79YB9EsHhOPVHSfQ2GFDaZYFV7rn7WRYjNrr4LWdJ9Bhc2h7N3nXW8wqy4bFqENTpz2sCaarPjiCDQca8H+fDtzBFo7D9R24/onNaOywY2pJJv5267yAQWJplkXrqJE3X7VtPVo65aszSvDQN6bj7PIcWG0O3PPKbgDAlTNKUei1JcQF7lWt3cfbtABLpmiumjVCu97sOdkWd92FAzntgxE5ebU0S90gz5+CdDPSzQa4xMA1CSdb++8nYjHqMbVU3VVz57HgI9VGq3rhzXffWWanGrU2yeMtXQEndE52r4xUN4W2OdYn7gvmVTNLtcKymigUjNW39+BEa7e2I6tcyQnlrnZrdTOEUNNncpaFXBk53tId0t3TFyfa8C/3nck97jz6hZMK8YOLJwIAfvH6F/iPv2xCvdWGSUUZeODqadp5DMThdHkFXbn9vl6SlYJxBWlwCWj5YqkmCdp6JW0kfBBpmhOt3XC4BMx99uPR0jR9umVOuVM000ZkwWTQobvXGXeD7zYeadTSfwCwfl9oO4D3FalgBICWqvnwoG+A1NbV61PEvfdke8AZS1urm3H/G3sGvNY4XQJ//fAIAOCWBWNh1Hvedq5xp2re/OwkPj3SBIdLoCw3xWenapNBp3WjfXI4tGBOCM8bv5yKGilCCHzv+W1o7LBhyiCBCODuqOkzFl6ump01OgclWSkwG/T4y+KzfF77S84d4/N9CjMtmDZC/f1/cLABx5q7sP1oCxQFuHJGCUqyUjAiOwUuoaaKEwWDkUHqRQD1j0iujgw0Fl629Y7OTfXZT2S2VsTaGvRx9V0ZURRFu1Ouae7C0Ub/b1j56Wbkp5shRPD7IAghtHqRr80egfnuHO2bn0V+dWSHu15kYlEG0s0GTCoKPa20pdq3XgQARroLWDtsDrSGsMunHJZ19axSn4LA/7xoPC6cVACbw4W27l7MHpWNl773JUwfoQaWg62M7K+1otPuRIbZgEnF/rseznMX8PW9wGozRk6zNI28MSjPS/N5/WhTWLt8gxFZZDwqNxXj3SnWeEvVyKmZchz4+gDtm8EKtXh1IOeNz4eiAAfqrD61N+/tq4PDJTChMB1pJj26e50BO+x++eYePLOxesBVh3f31KK6qQtZKUZcd3aZz9fOn1iA7FQjGqw2/H7dQQDAuX26UABoxa6hFrHKLTMA9TwrI9gp2Nxpx5GGTigK8Pwtc7UVvIFM1Np71b/Tt9zBiJyVBKjB9+rvnI2ROSm4cmYppo/M6vd9LnJ31bx/oF5bxT5nXJ62giJXR/yl4GTX5FB2lI+GsIKRxx57DOXl5bBYLJg3bx62bNky4ONfeeUVTJ48GRaLBdOnT8fbb78d1sFGw2D1IlIwdSPajJE+6Z7Zo7IBhDb8TCtg9RrC5j1rxLMy0v/uWasbCXKY2NGmLhxv6YZRr2De2FxcPVO9W3ljVzSCkVYAwJnuEfeTwkgryTHSMkUDqCtQ8t8q2LvjTyub8MHBBhh0CpZ9ZaLP13Q6BY8umo0FE/LxtdkjtLueIveLvcFqG/BNRW6Od+boHC2N0Ne5foKRLrtD+90nRZrGvTJi7XEM2oHkr3gVCFzAKtM0pVkWbeBfPAUjQgj8e7869+HHl052d0HYfNraQ+FwurTzi8TKSE6aCTNGZgMAPjzkSdWs/cLzBimDns/9jCZo6+7VusXeG2CA2gvuYu/FXxqNNPfqrmQy6PBV9xux/F7njO8fjMgi1s2VTYNupeCtus9sqP8XwdUR71V1OYJhMBO9ilhPtXVju3v11DsYAdT3o49+fCH+eP1sv9/nAnfdyIcHG/DaTrXuTF63AU+d2nY/HTV/WH8IT39Sjev/+qnP2P1YCzkYeemll7Bs2TKsWLECO3bswMyZM7Fw4ULU1/tffty4cSOuv/563HLLLdi5cyeuueYaXHPNNfjii+i2jwZLRsrjBlgZAYAxcvDZAIWW1Y2yXsT3TUQWse492T5gzYnUbXdqub78QMGIVjPS/7i1SaxBXpg/cr8ZnjkqB6kmg9qWptfhQJ017HbJw/VWtPa5kwU89SJyIJw81uMt3UHlNztsDnzhvmjNG5Pn87VQiliFEPgfd0//dXPL/P47ZqUa8fwt8/D7RbOQalIvovnpJigK4HAJNPs5P2nrUTnsrH+9iPSlcXnQKepFTb6xyhRNdqpRG+SWyDJTDFqL42AFhP7aegEgL13WjPg+Xw48K81O0Vaf9ke5o6aqsTPoYus9J9tR125DqkmPBRPy8eWJ6pvs+kG2AgjkUH0HbA4XUk16jI7QnkVfdteNyJXRDptD6565bHoxprlXAv3VjWxzp0sB9U3P32j35k67tprxjTkj/R7DNbNH+Pz/+WPz+j3mjNIsZFgMsNo87fLBkNdkeT/gXQszVNpO74UD38h6G+/V3vv252pgdHZ5DoqzLP0e23cAp7eZI7ORm2aCtceBw/UdMOl12vYj6vdUg5EdNS0+wVt7Ty9edm/U2d3rxJKnt8ZNl1LIwcgjjzyC2267DUuWLMHUqVOxatUqpKamYvXq1X4f/4c//AGXXnopfvSjH2HKlCn41a9+hTPPPBN/+tOfhnzwkTDYjBFpTMHgs0b6tvVKo3JTkZtmgt3pwt4gXkjyRW026LQ6EcDTvnqorkPLl5f7eROdEuIeNR+774pkQVtWihHnT1LzyeGsjry3tw5f+f2HuHbVJp/gy3vY2Znu1aKcNBOK3EWOwdzVbj/aAqc7r1ya7TsQLJT23uc2HcWOmlZYjDr810UTgjovQJ3lIodw1QWoGxFCaCsjfTtpvGVajJjpDlTl6kjf2TGJTlGUoNt7PZ00vuee487DN/epGZEzRkqzUzDRHYxEc8fqzZVNWPj7D3H1Y58ENWdGTsNcMCEfFqMeF09WizUHWkUYyIYD6ut03phcnzTWUMi6kY8PNcDpEvj3/nrYHS6MyU/DpKIMLS3pr6NG7g0FqB14//ZTD/POnlo4XQJnlGYGbBCYMypHa32X2yb0pdcp+JI7SPkkhDdPuTJyydRi6BTgixPtEdvDStvpPcB5+ePdUSM76fquigRDr1Nwvvt3BwAXTi7wuXmZUJiOTIsBXXYn9nmNeXh56zF02p0YX5iOCyYVoLvXiZuf2YqNIdbiRENIwYjdbsf27dtRUVHh+QY6HSoqKrBp0ya/z9m0aZPP4wFg4cKFAR8PADabDe3t7T4f0dDrVb0/UM0I4J2mGWBlpM9EVElRFG11ZP2++kE33JP1InIUvCTfoGQeMMNsQI6fFjLZUbPvVPugxZwOp0u7czlvgueP++pZasX7G7tPhlQQeqK1G/e8shtCqHdyT35UpX1NDjvLSTX6XJg8qZrBf89yP5q55f3vnoIdfLbpSBMe+JfatnvPJZN8KtWDMVgdxPGWbtS122DQKZjpXgYPpG/dSDLs1tuX9u81SEeNJxjxvTGQwV+n3ekTBMiakdJsT5qmsrEzqNXHUB2u78B3n98Ou9Plk54YiCxWlUHIhZMLoVPUFctwur7ed7dyytbOSJhVlo0MswEtXb344kSblqK5bFoxFEXRgpE9fopYZbpUzgt6d2//FIgs0PzqjMBvuDqdgv9wt6NWTA18bvJm6YkPK4NenZKjFs4cna2tpEaqkFVbGRnkvcNbaZYF6WYDHC6B3cfboCjAZdNCD0YA4IJJ3tdr39UlnU7RdnuX7xdOl9A23LvlvDFY9e05uHBSAXp6Xbj52a0hFwdHWkjBSGNjI5xOJ4qKinw+X1RUhNpa/7/g2trakB4PACtXrkRWVpb2UVZWFvCxQ3G0qQsOl0CqSe9Tve+PHE9d125Dp590gtMltPRAeX7/NxL5h/Gn9w9j3oPr8LPXPsfHhxr95j/91YsAnjcoWZ0/Oj/V71Le+MJ0mA06tPc4Bp0a+/mJNlh7HMi0GLQLD6BeQNNMehxv6fbZ1G4gvU4X/uvFnWjr7tVyqH/89yEtOJApmtmjcnyOe0oI+f6tVer3mDumf/pDjk4faBT+seYu3Pm37XC6BK6ZVYpbzhsT8LGBFGntqv7v9OWLf9qILKSYBp6g6l03IoQYsBYoUXk6agIHI7uOtWqpqr6vn8wUAwzulQDZ3tvT69T2oinNSkFxpgUZFgOcLqGtdnrrdboGDKp7ep34yd8/wwNv7u03KbXBasOSZ7b4TMvdPUj9V117Dz4/ob7ZyOAhN82kXQe8B1YFo62rV8v/y5HgkWDU63COe6bHu3trtT1P5Bvk2IJ0pLqLWL2LP73TpT+7fDIANdXj3VXT1GHTbnS+Osjd/10XjsMzS87Gfw6wSnntnDLMGJmFlq5efOfpLX73Kuqr2qso+rLpahojUqkauTIyWL2hN+89agDg7NG5flM0wTh/YgEyLQYUZZpxkZ8AVa7KyhED7+2txfGWbuSkGvG12SNgMeqxavEcXDS5UA1InoltQBKX3TTLly9HW1ub9nHs2LGo/BzvTpqB8nOAWj8gC+n6FkUB6l2a3emCSa/Tdtj19q25o3Dd2WXITjWiscOOFzbX4NtPbcb5v93QbyR4oGBkRE4KvFdnA3VbmAw6Ldc7WNGsnHtwzrh8n0LLFJMel7g3enozyJkjv3/vILYfbUGG2YB/3DEfXxqbi55eF37pHh4mO2lkd5E0Kcgal55ep3Y+c8cMsDIS4K6zy+7Abc9tQ0tXL6aPyMJD35gx6O/dn8JBRpzLPXNm+KmC72v2qGykGPVo7LDjQJ01qQaeSQWDpGkO11vxnae3QAigYkqh9u8rKYrSb9aI7P5IMeqRnWqEoigBi1iPNHRg9gPv4bvPbw+4KvnIewfx0rZjWP1JFS743QY8v6kaTpdAt92JW5/bhmPN3Ridl4rvnFMOANh9vHXAc5bBxsyR2T6v44unqDdm60Js8f3osJpGGV+Y7tP2GgkyVfPUx1Xo7nViZE6K1jqq1ynapFc5FBFQ60VkuvTCSYUYlZsKm8Pl0yb8zp46OF0C00Zk+q3J8mbQ63DBpMKA4xUA9Zr05E1nYUR2CqqbunDbc9sGTJcJIXyaCuTGdTtqWnGqbfC6soHYHd6r6sEHI4BnLDww8IrRYLJTTVj7gy/jzaXn+f13k3UjW6tbIITA6o+rAaibdcrHmw16/PnbZ6JiSiEyLP13rR9OIQUj+fn50Ov1qKvzzXnW1dWhuLjY73OKi4tDejwAmM1mZGZm+nxEgyfnF9wfkyxMrfbT3ivfRMpyU/x2T+SkmfDQN2Zg688r8Pwtc3H93FFIM+lxorVbG4EsNXqlabwZ9TqfOomB7p6DbSeWxavnTehfwX6VezjRvz47NWgF+wcHG/D4BnWWwEPfmIHReWn49TXTYNApWLevDu/trdOORXbSSJO83kQGunv97Hgb7E5134q+RcKAp2bkeEt3v04XIQTueWU39tdakZ9uxl9vnDPghW8gcpBXXYA3V7kyE0yqxWzQa11BHx9q1C5wydBJIw2UpjnR2o3FT21Ba1cvZpVl4w/X+e8e6DuF1TtFIwNKmY8/0KeI9cmPKtFhc+C9vXV47P3+W6tvP9qMJz6qBKC+xtu6e3HvP/fgij9+jO8+vw27j7UiO9WIp79ztrbKMdjKiCxSrZjie8cq//+nR5pCGkglgxt/d8BD9WV3elZusnbpGcU+QbpWxHrck5qS9SLzxqhj27/iHl7mnarRUjTTSyN2rIUZFjx789nItBiw/WgLfvjy7oBdbc2ddlh7HFAU9bVYlGnRVqbeGWKqpqa5E06XQJpJr62UBmtCofp3qqZoAr8PBqM0OyVgmnnGyCyY9Do0WG146/NT2FLdDKNeweL5o30eZzbo8fgNc/DqHedqI+tjIaRgxGQyYc6cOVi/fr32OZfLhfXr12P+/Pl+nzN//nyfxwPAe++9F/DxwymYGSPeBmrvrQpQL9KXurdAAVZ+fTqucuf5+qZBAq2MAL5vcAMFI7OCaCfutDm0fWIW+AlGzpuQj5xUIxo7bD6tf33Vt/dg2Uu7AAA3zBulRfvjCzNw64KxAICfvfa5z7Azb+ML06HXKWjr7kXtAMPEZPpj7pgcvysaJVkW6HUK7A6XVncj/eXDSrz9eS2MegWrvn2m39WrYMkXf6CVEZmWGpkT3M+QdSMfHGzACXcgk1Rpmkz/NTZNHTYsfmozTrX1YHxhOp7+ztn9Wj+l3D7BiNbW6xWc+1sZaem0a62PAPDouoM+3QPddifueeUzCAF848yRWLfsfDxw9RnISjFi36l2fHSoESaDDk/eeBbGFqRjpnu1q7qpy2+3GKCu4MmpqxdN9k1RjytIx+i8VNidLq1wfDAul8AH7uJV7zqBSCnLTfWp4ZLpDMlfEausF5nnDqQvcQcj/95fD4fT5U7RqP8Gg6VoQjW+MAOrFs+BUa/grc9P4b/f2e/3cdVeG07KGw/55j/UFt/D9Z7Gh1BXV+eNzVXTd5MKQ65XC4XFqNdmlNz7utq9esWMUm08gTeTQRfzG6CQ0zTLli3DE088gWeffRb79u3DHXfcgc7OTixZsgQAcOONN2L58uXa47///e9j7dq1ePjhh7F//37cf//92LZtG5YuXRq5swhTsJ00kqzyr/RTh3E0QFviQGRHSd8dfbUde9P7D9HxDUYC/yy5Qd++U+3otvtfytxc1YRep8DInBS/d/FGvQ7fOFNtx/PeSbKvn732BZo67ZhcnIF7r/DdIOq/Lh6P0iyLFmDJYWfezAa9VpE+UKpG7kdzdoAOFYPeM7nTu4i1p9eJx913xPdfdcaAHS7BkHf6dQMUsAKeQWyD8a4bkRN8izKid5EabjLtUtPchY8ONeD9A/VYv68OS57ZisqGTpRmWfDczQMPjcrpM2tEdpOVegWV2sqI19/QS9uOoafXhaklmbh2zki4BPD9Nbu0v8f/eWc/qho7UZxpwX1XToVBr8ON88vx/j0X4Pq5o1CSZcH/XjdL+5vJTjVpq3K7/czeANSpqz29LpRmWbSZP5KiKF5dNcGlaj470YamTjvSzYaAf/tDJVt8izLNmF3mu3Ip39D2nGyDyyXQZXfgM/e5yw6XOaNzkJtmQmtXL7ZWt2Dtnlq4hBrIRONN7pxx+fifb84AAPzlg0q/3T7Vfq7Jl7qDka3VzUPaq6ayMfTiVWnGyGy8d/eX8b8BZohE0lnu0QIt7kGQN58beo3ccAk5GFm0aBF+97vf4b777sOsWbOwa9curF27VitSrampwalTngKhc845By+88AL++te/YubMmfj73/+O119/HdOmTYvcWYRp/rg8LJiQr91RDUaujPir5NY2yAvhhSfTFZ8db4Pda2T0QCsjZUGujJRmWVCYYYbDJXxyvd7kbIEFE/IDRvfXzxsFQL3jkUvj3vadase6fXXQKcAfr5/dL/WRajLgvivP0P5/3xSNNGmQIlanS2gFfN7DzvqSRazedSPv7q1De48DI7JTcP3ZowI+N1ja4DM/qzhqYaX6+wt2ZWRycQby0kzaxmSj+kzwTXRyZaSqsROLn9qCJU9vxS3PbsNnx9uQm2bCc7fM69em3ZcnTaP+2570uzKipnNPtHbD2tMLh9OF5zepk0G/c245Hrh6GiYWpaPBasPdL+3CxiON2uZtD31juk9rZG6aCSu/Ph2bll+MS/t0O8iVvUCpGlkPcvGUIr+vK9kx8v6BwTvrAOB9d4pmwYR8n1HqkbTo7FHITTPhtgVj+/3tjStIR4pRj067E5WNndhxtBUO927U8m/coNdpKaR399Zqm78NpSZiMF+bPRIXuleK/A3v8jdqYWROKmaMzIJL+O/+CdaR+tBuZPsaX9j/piwazvbaiuLs8hy/01zjRVh/2UuXLsXRo0dhs9mwefNmzJs3T/vahg0b8Mwzz/g8/tprr8WBAwdgs9nwxRdf4PLLLx/SQUfKDy+ZhOdvmafdUQ1GLmX661AJNGNkIGPz05CdaoTN4fKZCdIQoGZE/f5qAGIe5O7Zu5040P4EsnL6vPGBl37HFaRj/tg8uASwZmv/QuJVH6h1IpdNLwmYb1x4RpGWUw7UCTDYBM19p9rRYXMgw2zQ3nT88bT3egKnV9w78n7jzBEReZP3bu3tm6+WqyLpZkPQQ8t0OsVn6mSkBlrFi8nFmbhyZikmFqVjcnEGpo3IxMyRWTh/YgGeu3muT3dBIJ40jXqHJ9M0Jdme10BWqlFbGTtYZ8W6fXU40dqN3DQTrppZihSTHo/fcCZSjHp8fLgR33l6KwDg+rlluCCEDhXZru0vGHG5hDZv46Ip/r/n2eW5yLAY0Nxpx4tbAq84StFo6e1ramkmdtz7FS2t6k2vU7RJrF+caMNmd3v9vLF5PsGWTNW8ufuUFhxEOkXTl1yx8lcbV9UkN130fT3JTqH/93n4wYhnZSS8YGS4zPG6+Qunc3A4xWU3TbySLYctXb0++WKXS2j7iQxWM+JNURSt0FR2mgghtE3y/K2MTCvNguKuuxjsjVWmavxt0FfT1IWDdR3QKZ5Ry4Hc8CV1NWHNlhqfra+PNXdpm8zdcf64gM9XFAWPfetMvP1fC7SgpC/v2Sj+yBTNnPLA49UB7yJW9fdxorVby99/c05kWsTl78XhEv32S5E/d2ROSki55PO8tkyPde420vQ6BX+8fjbevft8rP3Bl/Gv/1yAfy49D8/ePFcrjhxM35URmaYZ0WdFZaIW1HZgtXvV41tzR2krduMLM/Dra9RVWbvDhRHZKfjZ5VNCOh9tZeR4a7+C6x01Laht70GG2eB3kiigpj9vd79eVryxRxsZ70+D1aalRC6YGPl6kWBN95rEurlSFq/6rlAumFAAi1GHxg4bXAKYOTIr4p0/fclJzjtq+l/jAt0gXu6uifnkSGPAPXcGIoTQZowEW28YKzlpJvzwKxPx7S+NwlemDq1YNtoYjIQg1WTQ7ry8V0dOtffA7nDBqFdQmh1arn+29mJqBaAOdup2t6v5Wxkpz0/DumXn44kbzxr0ew+0MiI3wZs/Lm/QDZ4umVqM/HQT6q02n11Hn/yoEk6XwIIJ+YO+qZgMugE395JpmiMNHT4BjySLVwfLmWtpGvfKyKvbj0MI4EtjcyP2Jm/U67Q3x75FmaHWi0jnJvHKSCTkugefNXfaIYTwm6YBPCtsr+86gS1VzdDrFHz7S77dA9+YMxLfOaccFqMOv7t2JjIsoY3dP6M0EwadgsYOu7ZCI8ng/CtTiwbs1rrzgnH45pyRcLoE7vrbTnwWoFV4g3tVZNqIzKgWOw5Gvr63VTdrRfHz+gRbKSY9FngNTgxnsmioZpZlQa9TcKqtx6ddVwjhNUTPN2AYnZeGiilFEEK9hoWqqdOOdneXTqCpsvHkPy+egF9fM33Am7h4wGAkRHJ1xHvWiCxeLctJhSHEnK4W2bvrIRrdb26pJn3AzoJxBelBpQBmjMyCTgFOtvX47MoJeC6aV8wYvO3OZNDhWveERLnpVVOHDS+50x+3D7AqEqyROSlINxvQ6xT90mBCCG1lZKB6EcB31ojLJfDK9uMAoE14jBT5xtB3JLwnGAmtW2dkTqp2lxXL9rp4lZOm/r03ddrR1t2LLndRdkmfgVEy5Sr/Xi6bVux3qNT9V52BfQ9cqu1QHQqLUY/J7sLU3cc89VhOl9DaWa+YOfAbsaIoWPn16VgwIV8byS2n73qTI+AviuCgs3DIlZHd7vb6wgz/7fWXeK18DkcwkmoyaAHojqOt2udbunph7VFbp/0V599+vpqO+sf2E4NuU9CXbHwYkZ0S9ngA6o/BSIi09l6vKY9V2nJg6He0M8vUtMuJ1m7Ut/do9SL+UjShSjMbtFHru7xSNUcaOrDvVDsMOgWXnhHc0p0s/PzwYANqmrrw7MZq9PS6MH1E1qBpnmAoiqINA+rbUVPZ2ImmTjtMBt2gg8TksvCpth5sqmxCTXMX0s2GsEcuBxJoJLx3miZU/3vdbPzqmmkR+fdMNnleKyMn3Rvk5aWZ+r0Z9C1GX3JuecDvGc7AO0mrG/Fa0dha3Yx6qw2ZFsOAdViSUa/Dn789B1NLMtHYYe83VbTX6cKHB90tvVGsFwnGuII0WIyet4u+9SLSJWcUY3JxBr4+e0TUUzSSv1SNvFkszbL4DRjOKs/FnNE5sDtdeMadzgvWkYbEqBdJNAxGQuTd3nusuQv3v7EHv3lrH4DQ2nqlDIsRk9x3cztqWrVOmmC3pB6MTNXs9Cq2k5Xu547PHzRFI43KS/Wa1FiJZ91dCndcMG5IF3VvMo3z+s4TPrn4re673Fll2TAbBr4TKUg3w2TQwekS+MP6QwCAK2aUDDqWPVSBRsKHm6YB1KXwxV8aHbF/z2QiC1jbunu1wXD+OnDGF6ZrU4pnjMzS3qgizV9HjXxdLTyjWNupeDDpZgOeXnI2RmSnoLKxE1/+n/dx1ws78PrOE/j3/npYbQ7kppkG3eMo2gx6nTaJFehfLyJlpRix9gdfxiOLZg3Tkan7zgB9gpHGwRsKvvdldXXk+U+PhjSALtT5VBQcBiMhkisj6/fV4/zfvo9nNlajy+7E5OIMfGtueG2jWqFpTYvXjJHIBCOztVkmrdrn/uWuF7kixLa7G9xtvs9uOoq27l6M8RqxHAk3zS+HSa/Dv/fX4/VdnkFVW+SwsyBmLOh0irYqIZfqr41wigYIvN9KuGkaGpjcEFIIT5Gzv/osi1GvTbj8zjnlUQvsZJD/+Yk2OF0CDqdL2/Mk1HbWokwLnllyNspyU9Bhc+Ctz07hBy/twvee3w5A3YMkHvL93nVhXxobnXkn4ZAB554T7domidqohQFuECumFGFsQRqsPQ6sCaKrSTrSEPqeNDQ4BiMhkgVL3b1OuITa+//czXPx/76/IOxcvxx+tqOmZcAZI+GQ3TqfH2+Dw+nCwTorDtZ1wKTXaXvPBOviyYU+o49vWzA2ohfJCUUZ+K+LxwMA7n9jr5bL1YadDVIvInmvSowtSNP+fSNJGwnvtTLiPWOkLIyVEQrMoNch2x2QyAFXgabo/vbaGfj1NdNwTZ+dTCNpXEE60kx6dNmdOFzfgc1VzWjssCM71ehTjBysCUUZ+OCeC/HanefgzgvG+e5fMgy1F8GQwUh+uimu3ohH5aYiL80Eu9OFL06ogapng7zAr0OdTsF33a3MT31c5bdw3h+ujERH9KeuJJlxBWm444JxaO3qxY3zR2NKSeAOkWB5Dz8b5d4gLVJpmnEF6cgwG2C1OXCgzqrtyfDliflBz8GQDHodFp09Cv+7/hAKMsz4+pmRv9h/7/xx+H9f1GLPyXbc9/oerLhqKo63qGPkgw0qyrxWJa6dUxaVu2N/KyNyVSTDbEBmCl9akZabqk74/NwdjPRt65VmjMzGjCinNfQ6BTNGZmNTZRN2H2vV0qCXnlEc9mAynU7B7FE5mD0qBz++dDKONnWiscOGOaPjYxXi0mnF+Ndnp3DJVP/D3GJFUdR/t3X76rCzpgVzRudobb2Dpc6vmT0CD793EKfaevDm7pP4unvidCA2h1NLE46Po4AsGXBlJESKouAnl07Gyq9Pj0ggAgBj8tKQlaIOP/v4sFqwFqmVEZ1O0fap2VnTqnXRhDsZ8dYFY3Dd2WV4+NqZUakkN+p1+O03Z8KgU7B2Ty0ecO/4e0ZpVtAtmLJwTq9T8I0oBEyA134r7d7BiHqRGhHijBEKTm6fdurBprZGm6wb2Xa0GWu/CL47LVij89LiJhABgEyLEc/dPLdfq3Q88K4b8W7rHWzuk8Wo14qc//JB5YAbdQLqfCaXUGt9InWNJhWDkTig3hFlAwDq2mUBa3CFpcGQ+e01W2tQ2dgJk0GHiin+h48NJtNixEPfmKEVs0bD1NJM3HWhmq6RG1qFsieHnDp42bTiqM1mkCPh66092gVsKMWrNLjcPsXWJSHO9Im0WWVq2uL1XSfR0tWLvDRTXNVSnE48IxJa0drVi3Z3W28wHY43zBuNNJMeB+qsWit1IJ5OmjTecEQYg5E40bfqP5JRtwx0ZD71wkkFIQ96Gm53XTjep01z7pjguyLOLs/FumVfxu+unRmNQwPgKTDudQptEyoWr0ZXXp8APVCaZrjIlRG5r9Sl04pDnjNEkTFjpDr8rLZdbekH1Bk0wazeZqUYcYN7ted//31owNWRIyFurkrB4ysnTkQzGOnbFhjJpeRoMRnUdI1ep8CoV0LerXR8YUZUBxKZDDqvtIFaxDqUGSM0OO+VEaNeiVjHWbiKMy3avBkgMV5XySrVZNB2SH51hzroMJStOW5dMAZmgw47a1rxoXsDUX+8V0YoshiMxAk5/EyKVAErAOSlm7XlyhSjHhcH2MAr3kwfmYUXb/sSnrt5HvJi/Mbjj3wjkqk1pmmiS46EB4DiLEvMdzZWFEVbHSnIMA86HZiiS97QyVRLeX7wr8PCDItWC/OHdQcDro5UcmUkahiMxAnv4WcZFkPE7+pli+9FUwqRakqcTo+5Y3LDGtk9HGQ9ihx8xjRNdOWmeVKLgdp6h9uXJ6htvF8/c0RczAI5nclgxOHeSTuUlREA+N75Y2E26LCjplXbXNObEEJbGWFbb+QlzrvSaWD2qBzsr7VGZfn59gvGobvXiWVfmRjx73268h4J323njJFo814ZiXW9iPSteaMxsShDa8+n2Omb6h5o+qo/hRkWfGveKDz9STX+sO4Qzhuf71Ok2thhh9W9QV6ogQ4NjisjceTscvXFFI0ugcnFmfjL4rPialhRotOCkfYenGhV60U4YyR68rxqRkLdHTta9DoF88bmhT1bhCKnLDfFpwsxnB117zh/HMwGHbYdbcEnh5t8viZXRUbmcIO8aOArKI5cMaMUy74yEcsvmxLrQ6EgeNp7bTjmTtFwxkj0eBewxkuahuKHHH4m+dutdzCFmRZc797W4w/rPbUjvU4XPnBvWsgbuujgLVwcMRl0+K+LJ8T6MChIngLWHhavDgPvYCRe0jQUX84clYP39tahJMsS9uaYd1wwDi9sqcHW6ha8vusEqhq7sGZLjTZsb3JxZIZdki+ujBCFqdBrZYRtvdFnMeqRYVHvn/jvTP58ZWohzAYdvjwh/KGMRZkWbdPTu1/ajf9dfwj1Vhvy002468JxuOvCcZE6XPLClRGiMHlqRmw43qyujJSFsTRMwfvV1dNQ3dSJ8YVcKqf+xhdmYNsvKpAyxJqO288fh79vP44OmwNfGpuLG+aNxsIzimEy8P49WhiMEIVJ7k9jd7qw56S6eRvv2KPrmtnR24mXkkMkpksXZ1mw9gcL4HSJkLtyKDwMRojCZDbokZ1qRGtXL6qbmKYhSias/xpeXHMiGoKiDN8WU17AiIhCx2CEaAhkqgZQJ+dmpcT3BoRERPGIwQjREHhvaMhVESKi8DAYIRoCOfgMYL0IEVG4GIwQDUGhz8oIgxEionAwGCEaAt+VEaZpiIjCwWCEaAi4MkJENHQMRoiGgDUjRERDx6FnRENQkGGG2aCDSwiOgiciChODEaIhsBj1+OuNZ8HpciEzAmOoiYhORwxGiIbo/Inh7xBKRESsGSEiIqIYYzBCREREMcVghIiIiGKKwQgRERHFFIMRIiIiiikGI0RERBRTDEaIiIgophiMEBERUUwxGCEiIqKYYjBCREREMcVghIiIiGKKwQgRERHFFIMRIiIiiqmE2LVXCAEAaG9vj/GREBERUbDk+7Z8Hw8kIYIRq9UKACgrK4vxkRAREVGorFYrsrKyAn5dEYOFK3HA5XLh5MmTyMjIgKIoEfu+7e3tKCsrw7Fjx5CZmRmx7xtPkv0ceX6JL9nPkeeX+JL9HKN5fkIIWK1WlJaWQqcLXBmSECsjOp0OI0eOjNr3z8zMTMo/MG/Jfo48v8SX7OfI80t8yX6O0Tq/gVZEJBawEhERUUwxGCEiIqKYOq2DEbPZjBUrVsBsNsf6UKIm2c+R55f4kv0ceX6JL9nPMR7OLyEKWImIiCh5ndYrI0RERBR7DEaIiIgophiMEBERUUwxGCEiIqKYSthg5MMPP8SVV16J0tJSKIqC119/fdDnPPbYY5gyZQpSUlIwadIkPPfccz5ff+aZZ6Aois+HxWLxeYwQAvfddx9KSkqQkpKCiooKHDp0KJKnBiB25/ed73yn32MuvfTSSJ6aJhrnCACtra246667UFJSArPZjIkTJ+Ltt9/u933Ky8thsVgwb948bNmyJVKnpYnV+d1///39foeTJ0+O5KkBiM75XXDBBf2OXVEUfPWrX9UeM1yvQSB25zhcr8No/Y0++uijmDRpElJSUlBWVoa7774bPT09/b5PtF+DQOzOMZFfh729vXjggQcwbtw4WCwWzJw5E2vXrvX7fSL2OxQJ6u233xY///nPxauvvioAiNdee23Axz/++OMiIyNDrFmzRhw5ckS8+OKLIj09XbzxxhvaY55++mmRmZkpTp06pX3U1tb6fJ+HHnpIZGVliddff13s3r1bXHXVVWLMmDGiu7s7Kc7vpptuEpdeeqnPY5qbmyN6blI0ztFms4mzzjpLXH755eLjjz8WVVVVYsOGDWLXrl3aY9asWSNMJpNYvXq12LNnj7jttttEdna2qKurS4rzW7FihTjjjDN8focNDQ0RPbdonV9TU5PPcX/xxRdCr9eLp59+WnvMcL0GY3mOw/U6jMb5/e1vfxNms1n87W9/E1VVVeKdd94RJSUl4u6779YeM1yvwVieYyK/Dn/84x+L0tJS8dZbb4kjR46Ixx9/XFgsFrFjxw7tMZH+HSZsMOItmF/A/PnzxT333OPzuWXLlolzzz1X+/9PP/20yMrKCvg9XC6XKC4uFr/97W+1z7W2tgqz2SxefPHFsI49GMN1fkKoF8Grr746zCMNX6TO8c9//rMYO3assNvtAb/P3LlzxV133aX9f6fTKUpLS8XKlSvDO/ggDOf5rVixQsycOXMohxuySJ1fX7///e9FRkaG6OjoEELE7jUoxPCdoxCxeR1G6vzuuusucdFFFw34mFi8BoUY3nNM5NdhSUmJ+NOf/uTzmK9//evihhtu0P5/pH+HCZumCZXNZuuXkkhJScGWLVvQ29urfa6jowOjR49GWVkZrr76auzZs0f7WlVVFWpra1FRUaF9LisrC/PmzcOmTZuifxIDiMT5SRs2bEBhYSEmTZqEO+64A01NTVE//mAEc45vvPEG5s+fj7vuugtFRUWYNm0aHnzwQTidTgCA3W7H9u3bfX6HOp0OFRUVCfE7HOz8pEOHDqG0tBRjx47FDTfcgJqammE7j0CC/Rv19tRTT+G6665DWloagPh+DQKROUcpHl+HwZzfOeecg+3bt2tL9pWVlXj77bdx+eWXA4jv1yAQmXOUEvV1GOgxH3/8MYAo/Q7DCmHiDIKIBpcvXy6Ki4vFtm3bhMvlElu3bhVFRUUCgDh58qQQQoiNGzeKZ599VuzcuVNs2LBBXHHFFSIzM1McO3ZMCCHEJ5984vN46dprrxX/8R//EZVzE2L4zk8IIV588UXxz3/+U3z22WfitddeE1OmTBFnn322cDgcUTs/ISJ3jpMmTRJms1ncfPPNYtu2bWLNmjUiNzdX3H///UIIIU6cOCEAiI0bN/p87x/96Edi7ty5UTk3IYbv/IRQl21ffvllsXv3brF27Voxf/58MWrUKNHe3h735+dt8+bNAoDYvHmz9rlYvQaFGL5zFCI2r8NInt8f/vAHYTQahcFgEADE7bffrn0tVq9BIYbvHIVI7Nfh9ddfL6ZOnSoOHjwonE6nePfdd0VKSoowmUxCiOj8Dk+bYKSrq0ssWbJEGAwGodfrRWlpqfjxj38sAPSrm5DsdrsYN26c+MUvfiGEiO9gJBLn58+RI0cEALFu3bqhnMKgInWOEyZMEGVlZT4X7YcfflgUFxcLIeI7GInE+fnT0tIiMjMzxZNPPhmRc/EnGn+j3/3ud8X06dN9PhfvwUgkztGf4XgdRur83n//fVFUVCSeeOIJ8dlnn4lXX31VlJWViQceeEAIEf/BSCTO0Z9Eeh3W19eLq6++Wuh0OqHX68XEiRPFnXfeKSwWixAiOr/D0yZNk5KSgtWrV6OrqwvV1dWoqalBeXk5MjIyUFBQ4Pc5RqMRs2fPxuHDhwEAxcXFAIC6ujqfx9XV1Wlfi5VInJ8/Y8eORX5+/oCPGS7BnGNJSQkmTpwIvV6vPW/KlCmora2F3W5Hfn4+9Hp9wv4OBzs/f7KzszFx4sSY/w5D+Rvt7OzEmjVrcMstt/h8Pp5fg0BkztGfeHkdBnN+9957LxYvXoxbb70V06dPx9e+9jU8+OCDWLlyJVwuV1y/BoHInKM/ifQ6LCgowOuvv47Ozk4cPXoU+/fvR3p6OsaOHQsAUfkdnjbBiGQ0GjFy5Ejo9XqsWbMGV1xxBXQ6//8MTqcTn3/+OUpKSgAAY8aMQXFxMdavX689pr29HZs3b8b8+fOH5fgHM5Tz8+f48eNoamoa8DHDbaBzPPfcc3H48GGfC8LBgwdRUlICk8kEk8mEOXPm+PwOXS4X1q9fnxC/w8HOz5+Ojg4cOXIkbn6HwfyNvvLKK7DZbPj2t7/t8/lEeA0CQztHf+LtdTjQ+XV1dfU7Vxk8CyES4jUIDO0c/UnE16HFYsGIESPgcDjwj3/8A1dffTUAROd3GNZ6ShywWq1i586dYufOnQKAeOSRR8TOnTvF0aNHhRBC/PSnPxWLFy/WHn/gwAHx/PPPi4MHD4rNmzeLRYsWidzcXFFVVaU95pe//KV45513xJEjR8T27dvFddddJywWi9izZ4/2mIceekhkZ2dr+dyrr746Km2FsTg/q9Uq7rnnHrFp0yZRVVUl1q1bJ84880wxYcIE0dPTE9Hzi9Y51tTUiIyMDLF06VJx4MAB8a9//UsUFhaKX//619pj1qxZI8xms3jmmWfE3r17xXe/+12RnZ0dMJ2VaOf3wx/+UGzYsEFUVVWJTz75RFRUVIj8/HxRX18f9+cnnXfeeWLRokV+f+5wvQZjdY7D+TqMxvmtWLFCZGRkiBdffFFUVlaKd999V4wbN84njTZcr8FYnmMivw4//fRT8Y9//EMcOXJEfPjhh+Kiiy4SY8aMES0tLdpjIv07TNhg5P333xcA+n3cdNNNQgi1Ne7888/XHr93714xa9YskZKSIjIzM8XVV18t9u/f7/M9f/CDH4hRo0YJk8kkioqKxOWXX+7TVy2E2lp47733iqKiImE2m8XFF18sDhw4kBTn19XVJS655BJRUFAgjEajGD16tLjtttuicoGI1jkKoRbqzps3T5jNZjF27Fjxm9/8pl/h3x//+Eft32Lu3Lni008/TZrzW7RokSgpKREmk0mMGDFCLFq0SBw+fDhhzm///v0CgHj33Xf9/tzheg0KEZtzHM7XYTTOr7e3V9x///1i3LhxwmKxiLKyMnHnnXf6vJEJMTyvwVieYyK/Djds2CCmTJkizGazyMvLE4sXLxYnTpzo97Mj+TtUhAiwpkREREQ0DE67mhEiIiKKLwxGiIiIKKYYjBAREVFMMRghIiKimGIwQkRERDHFYISIiIhiisEIERERxRSDESIiIoopBiNEREQUUwxGiIiIKKYYjBAREVFMMRghIiKimPr/icQl5dv8VdUAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def make_y(x, per_pixel_unc=.05):\n",
    "    return line(x) + np.random.randn(len(x))*per_pixel_unc\n",
    "\n",
    "y = make_y(x)\n",
    "plt.plot(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "49086069-c56b-4f0a-847e-4f1a13b902ed",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "window = (1.962<x)&(x<1.97)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ec2b87f9-028b-4ab1-b020-b1fc115f9c2c",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.9671758916798359, 1.9670225484221437)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def centroid(x, y):\n",
    "    return np.sum(x*y)/np.sum(y)\n",
    "\n",
    "\n",
    "centroid(x,y), centroid(x[window],y[window])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b6607048-5ddd-4840-a97b-5820f45acaec",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.9670001267242403, 9.353676409979682e-05)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi4AAAGdCAYAAAA1/PiZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA18ElEQVR4nO3deXRU9f3/8ddMlskkkIQtCWGNNioKX0FBjNhqbQQrWlCKtU0roA2/arQgdSFasFoRqv2K+i2ouCB+BbV6ClJRrMaViuxYq5bFgqCYYEuTQJbJMp/fHxzu9w6LgM7kfoZ5Ps6Zcy53bibvD9fh9fbez3zGZ4wxAgAAiAN+rwsAAAA4UjQuAAAgbtC4AACAuEHjAgAA4gaNCwAAiBs0LgAAIG7QuAAAgLhB4wIAAOJGstcFfB3hcFg7duxQ+/bt5fP5vC4HAAAcAWOMdu/erfz8fPn9X+/aSVw2Ljt27FCPHj28LgMAAHwN27dvV/fu3b/Wz8Zl49K+fXtJeweemZnpcTUAoqWurk75+fmS9v4PSkZGhscVAYim2tpa9ejRw8nxryMuG5d9t4cyMzOj2rg0NTXp/vvvlyRNmDBBqampUXttAIeXlJTkbGdmZtK4AB5oiyz8JtM8fPH4JYu1tbXKyspSTU1NVBuXuro6tWvXTpK0Z88e/tEE2hjvQcB7sXwfRiO/4/KKS6wkJydrzJgxzjYAAInG9izkigsAa3DFBTi2RSO/WccFAADEDRoXAAAQN2hcXOrq6pSdna3s7GzV1dV5XQ4AAG3O9iy0b9aNx2pqarwuAQAAT9mchTQuLsFgUBs3bnS2AQBINLZnIY2Li9/vV2FhoddlAADgGduzkDkuAAAgbnDFxaW5uVlz5syRJI0fP14pKSkeVwQAQNuyPQtZgM6Fxa8Ab/EeBLzHkv9xJCkpST/84Q+dbQAAEo3tWcgVFwDW4IoLcGzjiguAuNZ78pKIP4ebGj2qBEC84FNFAAAgbtC4uNTX16tbt27q1q2b6uvrvS4HAIA2Z3sWcqvIxRijHTt2ONsAACQa27OQxsUlLS1N69atc7YBAEg0tmchjYtLUlKS+vfv73UZAAB4xvYsZI4LAACIG1xxcWlubtb8+fMlSSUlJdYtcwwAQKzZnoU0Li5NTU0aN26cJGn06NHWnSwAAGLN9iykcXFJSkrShRde6GwDAJBobM9CGheXtLQ0LVmy5PAHAgBwjLI9C5mcCwAA4gaNCwAAiBs0Li719fUqLCxUYWGhlcscAwAQa7ZnIXNcXIwx2rx5s7MNAECisT0LaVxc0tLStGzZMmcbAIBEY3sW0ri4JCUlaciQIV6XAQCAZ2zPQua4AACAuMEVF5eWlhYtXLhQknTJJZcoOZm/HgBAYrE9C+2qxmOhUEiXXXaZJGnPnj3WnSwgnvWebO+CVgD+j+1ZaFc1HvP7/TrnnHOcbQAAEo3tWUjj4hIMBvXmm296XQYAAJ6xPQvta6UAAAAOgcYFAADEDRoXl4aGBvXv31/9+/dXQ0OD1+UAANDmbM9C5ri4hMNhvf/++842AACJxvYspHFxSUtL01/+8hdnGwCARGN7FtK4uCQlJen888/3ugwAADxjexYe9RyXt99+WxdffLHy8/Pl8/m0aNGiiOeNMZo6daq6du2qYDCo4uJibdq0KeKYXbt2qaSkRJmZmcrOztZVV12lPXv2fKOBAACAY99RNy51dXU69dRTNWvWrIM+f/fdd+uBBx7QQw89pBUrVigjI0PDhg1TY2Ojc0xJSYk+/PBDvfrqq3rxxRf19ttva/z48V9/FFHS0tKiJUuWaMmSJWppafG6HAAA2pztWegzxpiv/cM+nxYuXKiRI0dK2nu1JT8/X7/61a90ww03SJJqamqUm5urJ554Qpdffrk+/vhjnXzyyVq1apUGDhwoSVq6dKkuvPBCffbZZ8rPzz/s762trVVWVpZqamqUmZn5dcs/QF1dndq1aydp7zLHGRkZUXttINEdyZL/4aZGbZ/5Q0m8BwGvxDILo5HfUf049JYtW1RZWani4mJnX1ZWlgYPHqzly5dLkpYvX67s7GynaZGk4uJi+f1+rVix4qCvGwqFVFtbG/GIBb/fr4EDB2rgwIFWLnMMAECs2Z6FUZ2cW1lZKUnKzc2N2J+bm+s8V1lZqZycnMgikpPVsWNH55j9TZ8+Xbfffns0Sz2oYDCoVatWxfz3AABgK9uz0L5W6iDKy8tVU1PjPLZv3+51SQAAwANRbVzy8vIkSVVVVRH7q6qqnOfy8vK0c+fOiOdbWlq0a9cu55j9BQIBZWZmRjwAAEDiiWrjUlBQoLy8PFVUVDj7amtrtWLFChUVFUmSioqKVF1drTVr1jjHvP766wqHwxo8eHA0yzlqDQ0NGjJkiIYMGWLlMscAAMSa7Vl41HNc9uzZo82bNzt/3rJli9avX6+OHTuqZ8+emjhxou68804VFhaqoKBAU6ZMUX5+vvPJoz59+uiCCy5QaWmpHnroITU3N+vaa6/V5ZdffkSfKIqlcDisd99919kGACDR2J6FR924rF69Wt/97nedP0+aNEmSNGbMGD3xxBO66aabVFdXp/Hjx6u6ulpnn322li5dGrFs8Pz583Xttdfqe9/7nvx+v0aNGqUHHnggCsP5ZgKBgBYuXOhsAwCQaGzPwm+0jotXYrWOC4DYYR0XANHIb76rCEBcOFTjs3XG8DauBICXaFxcWltb9c4770iSvv3tbyspKcnjigAAaFu2ZyGNi0tjY6Mzf4fL1ACARGR7FtK4uPh8Pp188snONgDv9JmyVP7UtMMfCCCqbM9CGheX9PR0ffjhh16XAQCAZ2zPwrhY8h8AAECicQEAAHGExsWloaFB559/vs4//3wrlzkGACDWbM9C5ri4hMNhvfbaa842AACJxvYspHFxCQQCeuqpp5xtAAASje1ZSOPikpycrJKSEq/LAADAM7ZnIXNcAABA3OCKi0tra6vWrl0rSTrttNOsW+YYAIBYsz0LaVxcGhsbdcYZZ0iyc5ljAABizfYspHFx8fl86tWrl7MNAECisT0LaVxc0tPTtXXrVq/LAADAM7ZnIZNzAQBA3KBxAQAAcYPGxaWxsVEjR47UyJEj1djY6HU5AAC0OduzkDkuLq2trXrhhRecbQAAEo3tWUjj4pKamqo5c+Y42wAAJBrbs5DGxSUlJUWlpaVelwEAgGdsz0LmuAAAgLjBFReXcDisjz/+WJLUp08f+f30dQCAxGJ7FtK4uDQ0NKhv376S7FzmGACAWLM9C2lc9tO5c2evSwAAwFM2ZyGNi0tGRoa+/PJLr8sAAMAztmehXTeuAAAAvgKNCwAAiBs0Li6NjY0qKSlRSUmJlcscAwAQa7ZnIY2LS2trqxYsWKAFCxZYucwxAACxZnsWMjnXJTU1VTNnznS2AQBINLZnIY2LS0pKiiZOnOh1GQAAeMb2LORWEQAAiBtccXEJh8Patm2bJKlnz57WLXMMAECs2Z6FNC4uDQ0NKigokGTnMscAAMSa7VlI47Kf9PR0r0sAAMBTNmchjYtLRkaG6urqvC4DAADP2J6Fdt24AgAA+Ao0LgAAIG7QuLiEQiGVlpaqtLRUoVDI63IAAGhztmchjYtLS0uLHn30UT366KNqaWnxuhwAANqc7VnI5FyXlJQU3Xnnnc42AACJxvYspHFxSU1N1a233up1GQAAeMb2LORWEQAAiBtccXExxuhf//qXJKlz587y+XweVwQAQNuyPQtpXFzq6+uVk5Mjyc5ljgEAiDXbs5BbRQAAIG5wxcUlIyNDxhivywAAwDO2ZyFXXAAAQNygcQEAAHGDxsUlFApp4sSJmjhxopXLHAMAEGu2ZyGNi0tLS4vuv/9+3X///VYucwwAQKzZnoVRb1xaW1s1ZcoUFRQUKBgM6vjjj9dvf/vbiIk+xhhNnTpVXbt2VTAYVHFxsTZt2hTtUo5aSkqKbrnlFt1yyy1WLnMMAECs2Z6FUf9U0e9+9zs9+OCDmjdvnk455RStXr1a48aNU1ZWln75y19Kku6++2498MADmjdvngoKCjRlyhQNGzZMH330kdLS0qJd0hFLTU3VtGnTPPv9AAB4zfYsjHrj8u6772rEiBEaPny4JKl37956+umntXLlSkl7r7bcd999+vWvf60RI0ZIkp588knl5uZq0aJFuvzyy6NdEgAAOEZE/VbRWWedpYqKCm3cuFGS9P7772vZsmX6/ve/L0nasmWLKisrVVxc7PxMVlaWBg8erOXLlx/0NUOhkGprayMesWCMUV1dnerq6qz+DDsAALFiexZG/YrL5MmTVVtbq5NOOklJSUlqbW3VtGnTVFJSIkmqrKyUJOXm5kb8XG5urvPc/qZPn67bb7892qUeoL6+Xu3atZNk5zLHAADEmu1ZGPUrLn/84x81f/58LViwQGvXrtW8efP0+9//XvPmzfvar1leXq6amhrnsX379ihWDAAA4kXUr7jceOONmjx5sjNXpV+/fvr00081ffp0jRkzRnl5eZKkqqoqde3a1fm5qqoq9e/f/6CvGQgEFAgEol3qAdLT07Vnzx5nGwCARGN7Fka9camvr5ffH3khJykpSeFwWJJUUFCgvLw8VVRUOI1KbW2tVqxYoauvvjra5RwVn89n3SUxAF+t9+QlB+zbOmO4B5UAxwbbszDqjcvFF1+sadOmqWfPnjrllFO0bt063Xvvvbryyisl7f0LmThxou68804VFhY6H4fOz8/XyJEjo10OAAA4hkS9cfmf//kfTZkyRddcc4127typ/Px8/b//9/80depU55ibbrpJdXV1Gj9+vKqrq3X22Wdr6dKlnq7hIklNTU3OJODbbrtNqampntYDAEBbsz0LfcbGzzodRm1trbKyslRTU6PMzMyovW5dXZ3VM6mBeHawWzr7Czc1avvMH0qSelz/vPypX+9/ZrhVBHx9sczCaOR31K+4xLPk5GRNmDDB2QYQnw7VJNHQAIdnexbaV5GHAoGA7rvvPq/LAADAM7ZnId8ODQAA4gaNCwAAiBs0Li51dXXy+Xzy+Xyqq6vzuhwAANqc7VlI4wIAAOIGk3Nd0tPTtXPnTmcbAIBEY3sW0ri4+Hw+denSxesyAADwjO1ZyK0iAAAQN7ji4tLU1KR77rlH0t5vubZtmWMAAGLN9iykcXFpbm7Wr3/9a0nSxIkTrTtZAADEmu1ZSOPikpycrJ///OfONgAAicb2LLSvIg8FAgE98sgjXpcBAIBnbM9CJucCAIC4QeMCAADiBo2LS11dnTIyMpSRkWHlMscAAMSa7VnIHJf91NfXe10CAACesjkLaVxcgsGgtmzZ4mwDAJBobM9CGhcXv9+v3r17e10GAACesT0LmeMCAADiBldcXJqbmzVr1ixJUllZmVJSUjyuCACAtmV7FtK4uDQ1Nen666+XJJWWllp3sgAAiDXbs5DGxSUpKUk/+clPnG0AABKN7VlI4+KSlpam+fPne10GAACesT0LmZwLAADiBo0LAACIGzQuLnV1derSpYu6dOli5TLHAADEmu1ZyByX/fzrX//yugQAADxlcxbSuLgEg0H9/e9/d7YBAEg0tmchjYuL3+/XKaec4nUZAAB4xvYsZI4LAACIG1xxcWlubtYTTzwhSRo7dqx1qwUCABBrtmchjYtLU1OTxo8fL0n6yU9+Yt3JAgAg1mzPQhoXl6SkJI0YMcLZBgAg0diehTQuLmlpaVq0aJHXZQAA4Bnbs5DJuQAAIG7QuAAAgLhB4+JSX1+v3r17q3fv3qqvr/e6HAAA2pztWcgcFxdjjD799FNnG8Ch9Z685KD7t84Y3saVAIgm27OQxsUlLS1NK1eudLYBAEg0tmchjYtLUlKSBg0a5HUZAAB4xvYspHEBEFWHuoUEANFA4+LS0tKiZ599VpL0ox/9SMnJ/PUAABKL7VloVzUeC4VC+ulPfypJGjlypHUnCwCAWLM9C+2qxmN+v1/FxcXONgAAicb2LKRxcQkGg3r11Ve9LgMAAM/YnoX2tVIAAACHQOMCAADiBo2LS319vU455RSdcsopVi5zDABArNmehcxxcTHG6KOPPnK2AQBINLZnIY2LS1pamt544w1nGwCARGN7FtK4uCQlJencc8/1ugwAADxjexYyxwUAAMSNmDQun3/+uX7605+qU6dOCgaD6tevn1avXu08b4zR1KlT1bVrVwWDQRUXF2vTpk2xKOWotLS0aNGiRVq0aJFaWlq8LgcAgDZnexZG/VbRf/7zHw0ZMkTf/e539fLLL6tLly7atGmTOnTo4Bxz991364EHHtC8efNUUFCgKVOmaNiwYfroo488vZ8WCoV0ySWXSJL27Nlj3TLHAADEmu1ZGPVqfve736lHjx6aO3eus6+goMDZNsbovvvu069//WuNGDFCkvTkk08qNzdXixYt0uWXXx7tko6Y3+/XWWed5WwDAJBobM/CqFe0ePFiDRw4UKNHj1ZOTo4GDBigRx55xHl+y5YtqqysdL4HQZKysrI0ePBgLV++/KCvGQqFVFtbG/GIhWAwqL/+9a/661//qmAwGJPfAQCAzWzPwqg3Lv/85z/14IMPqrCwUK+88oquvvpq/fKXv9S8efMkSZWVlZKk3NzciJ/Lzc11ntvf9OnTlZWV5Tx69OgR7bIBAEAciHrjEg6Hddppp+muu+7SgAEDNH78eJWWluqhhx762q9ZXl6umpoa57F9+/YoVgwAAOJF1BuXrl276uSTT47Y16dPH23btk2SlJeXJ0mqqqqKOKaqqsp5bn+BQECZmZkRj1hoaGjQoEGDNGjQIDU0NMTkdwAAYDPbszDqk3OHDBmiDRs2ROzbuHGjevXqJWnvRN28vDxVVFSof//+kqTa2lqtWLFCV199dbTLOSrhcNj52HY4HPa0FsAmvScv8boEAG3E9iyMeuNy/fXX66yzztJdd92lyy67TCtXrtScOXM0Z84cSZLP59PEiRN15513qrCw0Pk4dH5+vkaOHBntco5KIBDQiy++6GwDAJBobM/CqDcugwYN0sKFC1VeXq477rhDBQUFuu+++1RSUuIcc9NNN6murk7jx49XdXW1zj77bC1dutTz70RITk7W8OHDPa0BAAAv2Z6FMVlV5qKLLtJFF110yOd9Pp/uuOMO3XHHHbH49QAA4Bhl13J4HmttbdXrr78uSTrvvPOUlJTkcUUAALQt27OQxsWlsbFRQ4cOlbR3meOMjAyPKwIAoG3ZnoU0Li5+v1+nnnqqsw0AQKKxPQtpXFyCwaDWr1/vdRkAAHjG9iy0r5UCAAA4BBoXAAAQN2hcXBoaGnTuuefq3HPPtXKZYwAAYs32LGSOi0s4HNZbb73lbAMAkGhsz0IaF5dAIKA//vGPzjYAAInG9iykcXFJTk7W6NGjvS4DAADP2J6FzHEBAABxgysuLq2trXrvvfckSWeeeaZ1yxwDABBrtmchjYtLY2Ojzj77bEl2LnMMAECs2Z6FNC4uPp9P3/rWt5xtAAASje1ZSOPikp6erk2bNnldBgAAnrE9C5mcCwAA4gaNCwAAiBs0Li6NjY0aPny4hg8frsbGRq/LAQCgzdmehcxxcWltbdVLL73kbAMAkGhsz0IaF5fU1FTNnTvX2QYAINHYnoU0Li4pKSkaO3as12UAAOAZ27OQOS4AACBucMXFpbW1VR988IEkqV+/ftYtcwwAQKzZnoU0Li6NjY0aMGCAJDuXOQYAINZsz0IaFxefz6f8/HxnGwCARGN7FtK4uKSnp+vzzz/3ugwAADxjexYyORcAAMQNGhcAABA3aFxcGhsbNXr0aI0ePdrKZY4BAIg127OQxsWltbVVzz//vJ5//nkrlzkGACDWbM9CJue6pKam6g9/+IOzDQBAorE9C2lcXFJSUlRWVuZ1GQAAeMb2LORWEQAAiBtccXEJh8P65JNPJEnHH3+8/H76OgBAYrE9C2lcXBoaGnTCCSdIsnOZYwAAYs32LKRx2U9WVpbXJQAA4Cmbs5DGxSUjI0PV1dVelwEAgGdsz0K7blwBAAB8Ba64AEgYvScvOWDf1hnDPagEwNfFFReXUCiksWPHauzYsQqFQl6XAwBAm7M9C2lcXFpaWjRv3jzNmzdPLS0tXpcDAECbsz0LuVXkkpKSorvvvtvZBgAg0diehTQuLqmpqbrxxhu9LgMAAM/YnoXcKgIAAHGDKy4u4XBYX3zxhSSpa9eu1i1zDABArNmehTQuLg0NDerevbskO5c5BgAg1mzPQhqX/SQn81cCAEhsNmehvZV5ICMjQ83NzV6XAQCAZ2zPQrtuXAEAAHwFGhcAABA3aFxcQqGQysrKVFZWZuUyxwAAxJrtWUjj4tLS0qLZs2dr9uzZVi5zDABArNmehUzOdUlJSdFtt93mbAMAkGhsz8KYX3GZMWOGfD6fJk6c6OxrbGxUWVmZOnXqpHbt2mnUqFGqqqqKdSmHlZqaqt/85jf6zW9+o9TUVK/LAQCgzdmehTFtXFatWqWHH35Y//Vf/xWx//rrr9ef//xnPffcc3rrrbe0Y8cOXXrppbEsBQAAHANi1rjs2bNHJSUleuSRR9ShQwdnf01NjR577DHde++9Ou+883T66adr7ty5evfdd/Xee+/FqpwjYoxRdXW1qqurZYzxtBYAALxgexbGrHEpKyvT8OHDVVxcHLF/zZo1am5ujth/0kknqWfPnlq+fPlBXysUCqm2tjbiEQv19fXq0KGDOnTooPr6+pj8DgAAbGZ7FsZkcu4zzzyjtWvXatWqVQc8V1lZqdTUVGVnZ0fsz83NVWVl5UFfb/r06br99ttjUSoAAIgjUb/isn37dk2YMEHz589XWlpaVF6zvLxcNTU1zmP79u1Red39paenq6mpSU1NTUpPT4/J7wAAwGa2Z2HUr7isWbNGO3fu1Gmnnebsa21t1dtvv60//OEPeuWVV9TU1KTq6uqIqy5VVVXKy8s76GsGAgEFAoFol3oAn89n5Ue/AABoK7ZnYdQbl+9973v64IMPIvaNGzdOJ510km6++Wb16NFDKSkpqqio0KhRoyRJGzZs0LZt21RUVBTtcgAAwDEk6o1L+/bt1bdv34h9GRkZ6tSpk7P/qquu0qRJk9SxY0dlZmbquuuuU1FRkc4888xol3NUmpqadOutt0qSpk2bZuXn1wEAiCXbs9CTlXNnzpwpv9+vUaNGKRQKadiwYZo9e7YXpURobm7W73//e0myduEdAABiyfYsbJPG5c0334z4c1pammbNmqVZs2a1xa8/YikpKbrhhhucbQAAEo3tWch3Fbmkpqbqnnvu8boMAAA8Y3sW8u3QAAAgbnDFxcUY43yFd3Jysnw+n8cVAQDQtmzPQhoXl/r6erVr107S3u9aysjI8LgioG31nrzE6xLa3KHGvHXG8DauBLCD7VnIrSIAABA3uOLikp6erv/85z/ONgAAicb2LKRxcfH5fAd8+SMAAInE9izkVhEAAIgbXHFxaWpq0l133SVJuuWWW6xbLRAAgFizPQtpXFyam5t1++23S5JuvPFG604WAACxZnsW0ri4JCcn65prrnG2AQBINLZnoX0VeSgQCFj3/UkAALQl27OQybkAACBu0LgAAIC4QePiUldXp5SUFKWkpKiurs7rcgAAaHO2ZyFzXPaz74ulAABIVDZnIY2LSzAY1GeffeZsAwCQaGzPQhoXF7/fr27dunldBgAAnrE9C5njAgAA4gZXXFyampp0//33S5ImTJhg3WqBAADEmu1ZSOPi0tzcrJtuukmSdM0111h3sgAAiDXbs5DGxSU5OVljxoxxtgEASDS2Z6F9FXkoEAjoiSee8LoMAAA8Y3sWMjkXAADEDRoXAAAQN2hcXOrq6pSdna3s7GwrlzkGACDWbM9C5rjsp6amxusSAADwlM1ZSOPiEgwGtXHjRmcbAIBEY3sW0ri4+P1+FRYWel0GAACesT0LmeMCAADiBldcXJqbmzVnzhxJ0vjx45WSkuJxRQAAtC3bs5DGxaWpqUnXXnutJGns2LHWnSwAAGLN9iykcXFJSkrSD3/4Q2cbANx6T15y0P1bZwxv40qA2LE9C2lcXNLS0vTcc895XQYAAJ6xPQuZnAsAAOIGjQsAAIgbNC4u9fX16tatm7p166b6+nqvywEAoM3ZnoXMcXExxmjHjh3ONgAAicb2LKRxcUlLS9O6deucbQAAEo3tWUjj4pKUlKT+/ft7XQYAAJ6xPQuZ4wIAAOIGV1xcmpubNX/+fElSSUmJdasFAgAQa7ZnIY2LS1NTk8aNGydJGj16tHUnCwCAWLM9C2lcXJKSknThhRc62wAAJBrbs5DGxSUtLU1Llhz8u0gAAEgEtmchk3MBAEDcoHEBAABxg8bFpb6+XoWFhSosLLRymWMAAGLN9ixkjouLMUabN292toFjWe/J9t7DtgF/P0hUtmchjYtLWlqali1b5mwDAJBobM9CGheXpKQkDRkyxOsyAADwjO1ZyBwXAAAQN7ji4tLS0qKFCxdKki655BIlJ/PXAwBILLZnYdSvuEyfPl2DBg1S+/btlZOTo5EjR2rDhg0RxzQ2NqqsrEydOnVSu3btNGrUKFVVVUW7lKMWCoV02WWX6bLLLlMoFPK6HAAA2pztWRj1xuWtt95SWVmZ3nvvPb366qtqbm7W0KFDVVdX5xxz/fXX689//rOee+45vfXWW9qxY4cuvfTSaJdy1Px+v8455xydc8458vu5iwYASDy2Z2HUr/8sXbo04s9PPPGEcnJytGbNGn3nO99RTU2NHnvsMS1YsEDnnXeeJGnu3Lnq06eP3nvvPZ155pnRLumIBYNBvfnmm579fgAAvGZ7Fsa8laqpqZEkdezYUZK0Zs0aNTc3q7i42DnmpJNOUs+ePbV8+fKDvkYoFFJtbW3EAwAAJJ6YNi7hcFgTJ07UkCFD1LdvX0lSZWWlUlNTlZ2dHXFsbm6uKisrD/o606dPV1ZWlvPo0aNHLMsGAACWimnjUlZWpr///e965plnvtHrlJeXq6amxnls3749ShVGamhoUP/+/dW/f381NDTE5HcAAGAz27MwZp9xuvbaa/Xiiy/q7bffVvfu3Z39eXl5ampqUnV1dcRVl6qqKuXl5R30tQKBgAKBQKxKdYTDYb3//vvONgAAicb2LIx642KM0XXXXaeFCxfqzTffVEFBQcTzp59+ulJSUlRRUaFRo0ZJkjZs2KBt27apqKgo2uUclbS0NP3lL39xtgEASDS2Z2HUG5eysjItWLBAL7zwgtq3b+/MW8nKylIwGFRWVpauuuoqTZo0SR07dlRmZqauu+46FRUVefqJImnvMsfnn3++pzUAAOAl27Mw6o3Lgw8+KEk699xzI/bPnTtXY8eOlSTNnDlTfr9fo0aNUigU0rBhwzR79uxolwIAAI4xMblVdDhpaWmaNWuWZs2aFe1f/420tLTolVdekSQNGzbMumWOAQCINduz0K5qPBYKhXTRRRdJkvbs2WPdyQIAINZsz0K7qvGY3+/XwIEDnW0AABKN7VlI4+ISDAa1atUqr8sAAMAztmehfa0UAADAIdC4AACAuEHj4tLQ0KAhQ4ZoyJAhVi5zDABArNmehcxxcQmHw3r33XedbeBY0HvyEq9LABBHbM9CGheXQCCghQsXOtsAACQa27OQxsUlOTlZI0eO9LoMAAA8Y3sWMscFAADEDa64uLS2tuqdd96RJH37299WUlKSxxUBANC2bM9CGheXxsZGffe735W0d5njjIwMjysCAKBt2Z6FNC4uPp9PJ598srMNAECisT0LaVxc0tPT9eGHH3pdBgAAnrE9C5mcCwAA4gaNCwAAiBs0Li4NDQ06//zzdf7551u5zDEAALFmexYyx8UlHA7rtddec7YBAEg0tmchjYtLIBDQU0895WwDAJBobM9CGheX5ORklZSUeF0GAACesT0LmeMCAADiBldcXFpbW7V27VpJ0mmnnWbdMscAAMSa7VlI4+LS2NioM844Q5KdyxwDABBrtmchjYuLz+dTr169nG0AABKN7VlI4+KSnp6urVu3el0GAACesT0LmZwLAADiBo0LAACIGzQuLo2NjRo5cqRGjhypxsZGr8sBAKDN2Z6FzHFxaW1t1QsvvOBsA7bqPXnJQfdvnTG8jSsBcKyxPQtpXFxSU1M1Z84cZxsAgERjexbSuLikpKSotLTU6zIAAPCM7VnIHBcAABA3uOLiEg6H9fHHH0uS+vTpI7+fvg7A4R1szhHzjRCvbM9CGheXhoYG9e3bV5KdyxwDABBrtmchjct+Onfu7HUJAAB4yuYspHFxycjI0Jdfful1GcDXdqiPSQPAkbI9C+26cQUAAPAVaFwAAEDc4FaRS2Njo6666ipJ0mOPPaa0tDSPKwIQr1jdGPHK9izkiotLa2urFixYoAULFli5zDEAALFmexZyxcUlNTVVM2fOdLYBAEg0tmchjYtLSkqKJk6c6HUZAAB4xvYs5FYRAACIG1xxcQmHw9q2bZskqWfPntYtcwwAQKzZnoU0Li4NDQ0qKCiQZOcyxwAAxJrtWUjjsp/09HSvSwAAwFM2ZyGNi0tGRobq6uq8LgMAAM/YnoV23bgCAAD4CjQuAAAgbtC4uIRCIZWWlqq0tFShUMjrcgAAaHO2ZyGNi0tLS4seffRRPfroo2ppafG6HAAA2pztWcjkXJeUlBTdeeedzjYARNuhvnzxYPhCRnjB9iykcXFJTU3Vrbfe6nUZAAB4xvYs9PRW0axZs9S7d2+lpaVp8ODBWrlypZflAAAAy3nWuDz77LOaNGmSbrvtNq1du1annnqqhg0bpp07d3pVkowx+vLLL/Xll1/KGONZHQAAeMX2LPSscbn33ntVWlqqcePG6eSTT9ZDDz2k9PR0Pf74416VpPr6euXk5CgnJ0f19fWe1QEAgFdsz0JP5rg0NTVpzZo1Ki8vd/b5/X4VFxdr+fLlBxwfCoUiPpJVU1MjSaqtrY1qXe6VAmtra9Xa2hrV1weiJRyy7x+TaAg3Nf7fdqheMmEPq/FetP+NA45ELLNw33/T3+RKjieNy7/+9S+1trYqNzc3Yn9ubq7+8Y9/HHD89OnTdfvttx+wv0ePHjGrMT8/P2avDeDwPp99hdcleC7rPq8rQKKLVRbu3r1bWVlZX+tn4+JTReXl5Zo0aZLz53A4rF27dqlTp07y+XxR/V21tbXq0aOHtm/frszMzKi+tg0YX/w71sfI+OLfsT7GY318UuzGaIzR7t27v1FD5Enj0rlzZyUlJamqqipif1VVlfLy8g44PhAIKBAIROzLzs6OZYnKzMw8Zv+DlBjfseBYHyPji3/H+hiP9fFJsRnj173Sso8nk3NTU1N1+umnq6KiwtkXDodVUVGhoqIiL0oCAABxwLNbRZMmTdKYMWM0cOBAnXHGGbrvvvtUV1encePGeVUSAACwnGeNy49+9CN9+eWXmjp1qiorK9W/f38tXbr0gAm7bS0QCOi222474NbUsYLxxb9jfYyML/4d62M81scn2T1Gn7FxdRkAAICD4NuhAQBA3KBxAQAAcYPGBQAAxA0aFwAAEDfipnF5++23dfHFFys/P18+n0+LFi067M/MmjVLffr0UTAY1Iknnqgnn3zygGOqq6tVVlamrl27KhAI6IQTTtBLL70Uccznn3+un/70p+rUqZOCwaD69eun1atXO88bYzR16lR17dpVwWBQxcXF2rRpU8Rr7Nq1SyUlJcrMzFR2drauuuoq7dmzx/rxNTc36+abb1a/fv2UkZGh/Px8XXHFFdqxY0fEa/Tu3Vs+ny/iMWPGDOvHJ0ljx449oPYLLrgg4jUOd/5sH+P+49v3uOeee5xjbD2HB6vL5/OprKzMOaaxsVFlZWXq1KmT2rVrp1GjRh2wwOW2bds0fPhwpaenKycnRzfeeKNaWlqsH9+uXbt03XXX6cQTT1QwGFTPnj31y1/+0vnOtn0O9hrPPPNMxDG2jlGSzj333AOe/8UvfhHxe+L1HG7duvWQ78HnnnvOeZ3DnUOvxtfa2qopU6aooKBAwWBQxx9/vH77299GfN9QNHLwiJk48dJLL5lbb73V/OlPfzKSzMKFC7/y+NmzZ5v27dubZ555xnzyySfm6aefNu3atTOLFy92jgmFQmbgwIHmwgsvNMuWLTNbtmwxb775plm/fr1zzK5du0yvXr3M2LFjzYoVK8w///lP88orr5jNmzc7x8yYMcNkZWWZRYsWmffff9/84Ac/MAUFBaahocE55oILLjCnnnqqee+998w777xjvvWtb5kf//jH1o+vurraFBcXm2effdb84x//MMuXLzdnnHGGOf300yPq6dWrl7njjjvMF1984Tz27Nlj/fiMMWbMmDHmggsuiKh9165dEfUc7vzZPkb32L744gvz+OOPG5/PZz755BPnGFvP4c6dOyNqevXVV40k88YbbzjH/OIXvzA9evQwFRUVZvXq1ebMM880Z511lvN8S0uL6du3rykuLjbr1q0zL730kuncubMpLy+3fnwffPCBufTSS83ixYvN5s2bTUVFhSksLDSjRo2KqEeSmTt3bsRruf8NsnmMxhhzzjnnmNLS0ojjampqnOfj+Ry2tLQc8B68/fbbTbt27czu3buP+Bx6Nb5p06aZTp06mRdffNFs2bLFPPfcc6Zdu3bm/vvvd46JRg4eqbhpXNyO5IQVFRWZG264IWLfpEmTzJAhQ5w/P/jgg+a4444zTU1Nh3ydm2++2Zx99tmHfD4cDpu8vDxzzz33OPuqq6tNIBAwTz/9tDHGmI8++shIMqtWrXKOefnll43P5zOff/651eM7mJUrVxpJ5tNPP3X29erVy8ycOfOIft628Y0ZM8aMGDHikM8f7fkzxr4x7m/EiBHmvPPOi9hn6znc34QJE8zxxx9vwuGwMWbv+y0lJcU899xzzjEff/yxkWSWL19ujNn7D77f7zeVlZURvzszM9OEQiGrx3cwf/zjH01qaqppbm4+qprdbBvjOeecYyZMmHDInznWzmH//v3NlVdeedQ1H82x0Rrf8OHDD6j10ksvNSUlJcaY2OTgV4mbW0VHKxQKKS0tLWJfMBjUypUr1dzcLElavHixioqKVFZWptzcXPXt21d33XVXxFd4L168WAMHDtTo0aOVk5OjAQMG6JFHHnGe37JliyorK1VcXOzsy8rK0uDBg7V8+XJJ0vLly5Wdna2BAwc6xxQXF8vv92vFihVWj+9gampq5PP5Dvi+qBkzZqhTp04aMGCA7rnnnohLuLaP780331ROTo5OPPFEXX311fr3v//tPBeL8+fFGPepqqrSkiVLdNVVVx3wnI3n0K2pqUlPPfWUrrzySucLVtesWaPm5uaI9+BJJ52knj17RrwH+/XrF7HA5bBhw1RbW6sPP/zQ6vEdTE1NjTIzM5WcHLmGaFlZmTp37qwzzjhDjz/+eMSl/HgY4/z589W5c2f17dtX5eXlqq+vd547ls7hmjVrtH79+oO+B6N5DqM1vrPOOksVFRXauHGjJOn999/XsmXL9P3vf1+SBzl4VG2OJXQEnWZ5ebnJy8szq1evNuFw2Kxatcrk5uYaSWbHjh3GGGNOPPFEEwgEzJVXXmlWr15tnnnmGdOxY0fzm9/8xnmdQCBgAoGAKS8vN2vXrjUPP/ywSUtLM0888YQxxpi//vWvEa+5z+jRo81ll11mjNl7me2EE044oMYuXbqY2bNnWz2+/TU0NJjTTjvN/OQnP4nY/9///d/mjTfeMO+//7558MEHTXZ2trn++usP+hq2je/pp582L7zwgvnb3/5mFi5caPr06WMGDRpkWlpajDFHf/5sHKPb7373O9OhQ4cDbiPYeg7dnn32WZOUlBTxf2jz5883qampBxw7aNAgc9NNNxljjCktLTVDhw6NeL6urs5IMi+99JLV49vfl19+aXr27GluueWWiP133HGHWbZsmVm7dq2ZMWOGCQQCEZfybR/jww8/bJYuXWr+9re/maeeesp069bNXHLJJc7zx9I5vPrqq02fPn0O2H8057Atx9fa2mpuvvlm4/P5THJysvH5fOauu+5yno9FDn6VY7Zxqa+vN+PGjTPJyckmKSnJ5Ofnm5tuuslIci41FhYWmh49ejgBZczef7zz8vKcP6ekpJiioqKI177uuuvMmWeeaYzxrnFpq/G5NTU1mYsvvtgMGDAg4t7zwTz22GMmOTnZNDY2xs349vnkk0+MJPPaa68ZY2LXuHg1xhNPPNFce+21X1mbMfacQ7ehQ4eaiy66KGKfV41LW43PraamxpxxxhnmggsuOOytiylTppju3bsf8nlbx7hPRUWFkeTM1TpWzmF9fb3Jysoyv//977+yNmO++hy25fiefvpp0717d/P000+bv/3tb+bJJ580HTt2jOn/wH+VY/ZWUTAY1OOPP676+npt3bpV27ZtU+/evdW+fXt16dJFktS1a1edcMIJSkpKcn6uT58+qqysVFNTk3PMySefHPHaffr00bZt2yRJeXl5knTAJxiqqqqc5/Ly8rRz586I51taWrRr1y7nGFvHt09zc7Muu+wyffrpp3r11VcP+zXngwcPVktLi7Zu3RoX43M77rjj1LlzZ23evFlSbM6fV2N85513tGHDBv385z8/bH22nMN9Pv30U7322msH1J6Xl6empiZVV1dH7N//PXiw9+i+52we3z67d+/WBRdcoPbt22vhwoVKSUn5yvoGDx6szz77TKFQ6GuNz4sx7l+/pIj3YbyfQ0l6/vnnVV9fryuuuOKw9X3Tcxit8d14442aPHmyLr/8cvXr108/+9nPdP3112v69OmS2j4Hj9nGZZ+UlBR1795dSUlJeuaZZ3TRRRfJ79877CFDhmjz5s0Kh8PO8Rs3blTXrl2VmprqHLNhw4aI19y4caN69eolSSooKFBeXp4qKiqc52tra7VixQoVFRVJkoqKilRdXa01a9Y4x7z++usKh8POm9PW8Un/17Rs2rRJr732mjp16nTYutavXy+/36+cnBzrx7e/zz77TP/+97/VtWtXSbE9f209xscee0ynn366Tj311MPWZcs53Gfu3LnKycnR8OHDI/affvrpSklJiXgPbtiwQdu2bYt4D37wwQcR/3Dua8D3b/psG5+099+UoUOHKjU1VYsXLz5g3sLBrF+/Xh06dIjKl+S1xRgPVr+kiPdhPJ/DfR577DH94Ac/cBqHrxKtc/hNx1dfX+8cv09SUpLzM22eg0d1fcZDu3fvNuvWrTPr1q0zksy9995r1q1b53yyZfLkyeZnP/uZc/yGDRvM//7v/5qNGzeaFStWmB/96EemY8eOZsuWLc4x27ZtM+3btzfXXnut2bBhg3nxxRdNTk6OufPOO51jVq5caZKTk820adPMpk2bzPz58016erp56qmnnGNmzJhhsrOznXkSI0aMOOjHwAYMGGBWrFhhli1bZgoLCyM+Bmbr+JqamswPfvAD0717d7N+/fqIj+ntm8n/7rvvmpkzZ5r169ebTz75xDz11FOmS5cu5oorrrB+fLt37zY33HCDWb58udmyZYt57bXXzGmnnWYKCwsjbpEc7vzZPMZ9ampqTHp6unnwwQfN/mw+h8bsvcfes2dPc/PNNx9QuzF7Pw7ds2dP8/rrr5vVq1eboqKiiNtn+z5KO3ToULN+/XqzdOlS06VLl4iP0to6vpqaGjN48GDTr18/s3nz5oj34L7L+4sXLzaPPPKI+eCDD8ymTZvM7NmzTXp6upk6dWrEa9k6xs2bN5s77rjDrF692mzZssW88MIL5rjjjjPf+c53jolzuM+mTZuMz+czL7/88gHPHck59Gp8Y8aMMd26dXM+Dv2nP/3JdO7c2bkVa0x0cvBIxU3j8sYbbxhJBzzGjBljjNn7F3vOOec4x3/00Uemf//+JhgMmszMTDNixAjzj3/844DXfffdd83gwYNNIBAwxx13nJk2bVrEvT5jjPnzn/9s+vbtawKBgDnppJPMnDlzIp4Ph8NmypQpJjc31wQCAfO9733PbNiwIeKYf//73+bHP/6xadeuncnMzDTjxo2L+Py+rePbsmXLQeuSa42CNWvWmMGDB5usrCyTlpZm+vTpY+66666I4Ld1fPX19Wbo0KGmS5cuJiUlxfTq1cuUlpZGfOTySM6fzWPc5+GHHzbBYNBUV1cf8Jzt5/CVV14xkg54X+3T0NBgrrnmGtOhQweTnp5uLrnkEvPFF19EHLN161bz/e9/3wSDQdO5c2fzq1/9KuLjxLaO71B1SXIC6OWXXzb9+/c37dq1MxkZGebUU081Dz30kGltbT2i1/J6jNu2bTPf+c53TMeOHU0gEDDf+ta3zI033njAXLp4PYf7lJeXmx49ehxwXow5snPo1fhqa2vNhAkTTM+ePU1aWpo57rjjzK233hrxMfRo5OCR8hnzDT8vBwAA0EaO+TkuAADg2EHjAgAA4gaNCwAAiBs0LgAAIG7QuAAAgLhB4wIAAOIGjQsAAIgbNC4AACBu0LgAAIC4QeMCAADiBo0LAACIGzQuAAAgbvx/h9zSei0tE+wAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "centroids = np.array([centroid(x[window], make_y(x)[window]) for _ in range(1000)])\n",
    "plt.hist(centroids, bins='auto')\n",
    "plt.axvline(line.mean.value, c='k')\n",
    "plt.axvline(line.mean.value + line.stddev.value, c='k', ls=':')\n",
    "plt.axvline(line.mean.value - line.stddev.value, c='k', ls=':')\n",
    "np.mean(centroids), np.std(centroids)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4dac9cd-abd0-4649-8085-192f55de0765",
   "metadata": {},
   "source": [
    "So the \"correct\" uncertainty answer should be ~ .0001 with that window."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b10c426a-3d6b-42d4-ae96-7768d743ca78",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.6941356767952795"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def centroid_uncertainty(x, y, yunc=.05):\n",
    "    N = np.sum(y)\n",
    "    innerterm1 = x/N\n",
    "    innerterm2 = np.sum(y*x*np.log(N))\n",
    "    tosum = (yunc*(innerterm1 + innerterm2))**2\n",
    "    return np.sqrt(np.sum(tosum))\n",
    "centroid_uncertainty(x[window], y[window])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dfd67bee-f489-44a4-90c6-efd607b9b5dc",
   "metadata": {},
   "source": [
    "New derivation... yields\n",
    "\n",
    "$ c \\equiv \\frac{\\sum_i x_i y_i}{\\sum_i y_i}$\n",
    "\n",
    "$ \\Delta c = \\sqrt{\\sum_j  \\frac{\\Delta y_j^2}{(\\sum_i y_i)^2} (x_j - c)^2}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f1ade972-5f5f-4867-bc79-f4d181991350",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9.68704826953507e-05"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def centroid_uncertainty_fixed(x, y, yunc=.05):\n",
    "    N = np.sum(y)\n",
    "    s2 = np.sum(yunc**2*(x-centroid(x,y))**2)*N**-2\n",
    "    return np.sqrt(s2)\n",
    "centroid_uncertainty_fixed(x[window], y[window])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d8f5027-a841-4ff0-8e76-2b85fd23bcfc",
   "metadata": {},
   "source": [
    "Well that looks about right!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4bf4e00-c34a-4ad6-a7fd-5faefefcd775",
   "metadata": {},
   "source": [
    "Now lets try running Ricky's https://github.com/astropy/specutils/pull/1057 (SHA 264b5c55db325889506e75f4fae9d929db44c50d at the time of this writing)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "7b5fe1ad-f4f5-4f14-86c9-05678f947c5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from specutils import Spectrum1D, SpectralRegion\n",
    "from specutils.analysis import location\n",
    "from astropy import units as u\n",
    "from astropy import nddata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "181435f7-97e5-427c-9420-ab17ab3950bb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n",
      "Analytic\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(<Quantity 1.96702255 micron>, <Quantity 9.68704827e-05 micron>)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spec = Spectrum1D(y*u.Jy, x*u.micron, uncertainty=nddata.StdDevUncertainty([0.05]*len(y)))\n",
    "reg = SpectralRegion(1.962*u.micron, 1.97*u.micron)\n",
    "c = location.centroid(spec, reg, analytic_uncertainty=True)\n",
    "c, c.uncertainty"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b82b904c-2e10-48e8-b2a7-8f304d75c2dc",
   "metadata": {},
   "source": [
    "For kicks, lets see whether we also get consistent results if we centroid on the *whole* spectrum, not the windowed version:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "4e9a9f78-e80a-4075-b575-2d71f8549607",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0010084652456220772"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "centroid_uncertainty_fixed(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "6b275742-29fc-4c3f-b379-827fd05b236f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.9670090581464357, 0.0009322443601365744)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA0fklEQVR4nO3de3xU1b338e+eSyYhkCCX3LgGi4KAoBUxwqlUo9SiD1SK8DT2KCicp42WyzkqHEVbi4K0VUQpeCviBa+taL3gBVs8KnITb8UCtggUmqBVEnObJDPr+YPjbqO0gpmwZu/9eb9e83qt7Nmz89uL3zC/rL1mbccYYwQAAJBGQrYDAAAA+DwKFAAAkHYoUAAAQNqhQAEAAGmHAgUAAKQdChQAAJB2KFAAAEDaoUABAABpJ2I7gK8imUxq79696tChgxzHsR0OAAA4BMYYffrppyoqKlIo9K/HSDxZoOzdu1c9evSwHQYAAPgKdu/ere7du//LfTxZoHTo0EHSgRPMycmxHA0A22pra1VUVCTpwB8w2dnZliMCcDDV1dXq0aOH+zn+r3iyQPnssk5OTg4FCgKvsbFRt9xyiyRp2rRpysjIsBzRkRcOh912Tk5OIAsU8gBecijTMxwv3iywurpaubm5qqqqokBB4NXW1qp9+/aSpJqamkB+ONMH9AG84XA+vz05ggLg7yKRiC688EK3jWAiD+A3jKAA8DxGDwBvOJzPb9ZBAQAAaYcCBQAApB0KFMDjamtr1bFjR3Xs2FG1tbW2w4El5AH8hplUgA9UVVXZDgFpgDyAn1CgAB6XlZWlbdu2uW0EE3kAv6FAATwuFAqpb9++tsOAZeQB/IY5KAAAIO0wggJ4XFNTk+644w5J0tSpUxWNRi1HBBvIA/gNC7UBHsciZfSBRB/AG1jqHgiQcDis7373u24bwUQewG8YQQHgeYweAN7ACAqAQOs96+lWvf6D+aNTFAmAr4pv8QAAgLRDgQJ4XF1dnbp166Zu3bqprq7OdjiwhDyA33CJB/A4Y4z27t3rthFM5AH8hgIF8LjMzExt3rzZbSOYyAP4DQUK4HHhcFhDhgyxHQYsIw/gN8xBAQAAaYcRFMDjmpqa9MADD0iSysrKWOI8oMgD+A0FCuBxjY2NmjRpkiRp/PjxfDAFFHkAv6FAATwuHA7r29/+tttGMJEH8BsKFMDjMjMz9fTTrVs5Nd0c7kqwycYGt91/ziqFMoL3LRY/5gGCjUmyAAAg7VCgAACAtEOBAnhcXV2d+vbtq759+7LEeYCRB/Ab5qAAHmeM0fvvv++2EUzkAfyGAgXwuMzMTL3yyituG8FEHsBvKFAAjwuHwxo+fLjtMGAZeQC/oUABgM853K85f94H80enKBIguChQAI9rbm7W448/Lkn6zne+o0iEt3UQkQfwm8P+Fs/LL7+sc889V0VFRXIcRytXrmzxvDFG11xzjQoLC5WVlaXS0lJt3769xT4ff/yxysrKlJOTo44dO+riiy9WTU1Nq04ECKp4PK7zzz9f559/vuLxuO1wYAl5AL857AKltrZWgwcP1uLFiw/6/IIFC7Ro0SItXbpU69atU3Z2tkaNGqWGhr+v9FhWVqY//OEPeuGFF/TUU0/p5Zdf1tSpU7/6WQABFgqFdNppp+m0005TKMTKAUFFHsBvDnsM8Oyzz9bZZ5990OeMMVq4cKGuvvpqjRkzRpJ07733Kj8/XytXrtTEiRP13nvvadWqVdqwYYNOOukkSdKtt96qb3/72/r5z3+uoqKiVpwOEDxZWVn6/e9/bzsMWEYewG9SWmbv2LFDFRUVKi0tdbfl5uZq2LBhWrt2rSRp7dq16tixo1ucSFJpaalCoZDWrVt30OPG43FVV1e3eAAAAP9KaYFSUVEhScrPz2+xPT8/332uoqJCeXl5LZ6PRCLq1KmTu8/nzZs3T7m5ue6jR48eqQwbAACkGU9cqJw9e7aqqqrcx+7du22HBKSN+vp6DRkyREOGDFF9fb3tcGAJeQC/Sen30AoKCiRJlZWVKiwsdLdXVlZqyJAh7j779u1r8brm5mZ9/PHH7us/LxaLKRaLpTJUwDeSyaTeeustt41gIg/gNykdQSkuLlZBQYFWr17tbquurta6detUUlIiSSopKdH+/fu1adMmd5+XXnpJyWRSw4YNS2U4QCBkZmbq+eef1/PPP88S5wFGHsBvDnsEpaamxr0hlXRgYuybb76pTp06qWfPnpo+fbrmzp2rvn37qri4WHPmzFFRUZHGjh0rSerfv7++9a1vacqUKVq6dKmampp06aWXauLEiXyDB/gKwuGwzjzzTNthtNDalVhx+NIxD4DWOOwCZePGjfrmN7/p/jxz5kxJ0oUXXqh77rlHV1xxhWprazV16lTt379fI0aM0KpVq1pU9A888IAuvfRSnXHGGQqFQho3bpwWLVqUgtMBAAB+cNgFysiRI//lrbwdx9F1112n66677p/u06lTJ61YseJwfzWAg2hubtZzzz0nSRo1ahRLnAcUeQC/IYMBj4vH4zrnnHMkHbgEywdTMJEH8BsyGPC4UCjkLnzIEufBRR7AbyhQAI/LysrShg0bbIcBy8gD+A1lNgAASDsUKAAAIO1QoAAeV19fr+HDh2v48OEscR5g5AH8hjkogMclk0m99tprbhvBRB7AbyhQAI+LxWJ6/PHH3TaCiTyA31CgAB4XiUTcW0kguMgD+A1zUAAAQNphBAXwuEQiof/5n/+RJP3bv/2bwuGw5YhgA3kAv6FAATyuoaHBvYFnTU2NsrOzLUcEG8gD+A0FCuBxjuPouOOOc9sIJvIAfkOBAnhcu3bt9Ic//MF2GLCMPIDfMEkWAACkHQoUAACQdihQAI+rr6/XmWeeqTPPPJMlzgOMPIDfMAcF8LhkMqkXX3zRbSOYyAP4DQUK4HGxWEz333+/20YwkQfwGwoUwOMikYjKyspshwHLyAP4DXNQAABA2mEEBfC4RCKhN954Q5J04oknssR5QJEH8BsKFMDjGhoadPLJJ0tiifMgIw/gNxQogMc5jqNevXq5bQQTeQC/oUABPK5du3b64IMPbIcBy8gD+A2TZAEAQNqhQAEAAGmHAgXwuIaGBo0dO1Zjx45VQ0OD7XBgCXkAv2EOCuBxiURCTzzxhNtGMJEH8BsKFMDjMjIydMcdd7htBBN5AL+hQAE8LhqNasqUKbbDgGXkAfyGOSgAACDtMIICeFwymdR7770nSerfv79CIf7uCCLyAH5DgQJ4XH19vQYOHCiJJc6DjDyA31CgAD7QpUsX2yEgDZAH8BMKFMDjsrOz9eGHH9oOA5aRB/AbLlICAIC0Q4ECAADSDgUK4HENDQ0qKytTWVkZS5wHGHkAv6FAATwukUhoxYoVWrFiBUucBxh5AL9hkizgcRkZGbr55pvdNoKJPIDfUKAAHheNRjV9+nTbYcAy8gB+wyUeAACQdhhBATwumUxq165dkqSePXuyxHlAkQfwGwoUwOPq6+tVXFwsiSXOg4w8gN9QoAA+0K5dO9shIA2QB/ATChTA47Kzs1VbW2s7DFhGHsBvuEgJAADSDgUKAABIO1ziATwuHo/r0ksvlSTddtttisVirT5m71lPt/oYOLLaIg8AmxhBATyuublZd911l+666y41NzfbDgeWkAfwG0ZQAI+LRqOaO3eu20YwkQfwGwoUwOMyMjJ01VVX2Q4DlpEH8BsKFABIsVTM4flg/ugURAJ4FwUK4HHGGH300UeSpC5dushxHMsRwQbyAH6T8kmyiURCc+bMUXFxsbKysnT00Ufrpz/9qYwx7j7GGF1zzTUqLCxUVlaWSktLtX379lSHAgRCXV2d8vLylJeXp7q6OtvhwBLyAH6T8gLlxhtv1JIlS3Tbbbfpvffe04033qgFCxbo1ltvdfdZsGCBFi1apKVLl2rdunXKzs7WqFGj1NDQkOpwAACAB6X8Es9rr72mMWPGaPToA9dPe/furQcffFDr16+XdGD0ZOHChbr66qs1ZswYSdK9996r/Px8rVy5UhMnTkx1SICvZWdntxihRDCRB/CblI+gnHrqqVq9erW2bdsmSXrrrbf0yiuv6Oyzz5Yk7dixQxUVFSotLXVfk5ubq2HDhmnt2rUHPWY8Hld1dXWLBwAA8K+Uj6DMmjVL1dXV6tevn8LhsBKJhK6//nqVlZVJkioqKiRJ+fn5LV6Xn5/vPvd58+bN009+8pNUhwoAANJUykdQHnnkET3wwANasWKF3njjDS1fvlw///nPtXz58q98zNmzZ6uqqsp97N69O4URA94Wj8c1ffp0TZ8+XfF43HY4sIQ8gN+kfATl8ssv16xZs9y5JIMGDdLOnTs1b948XXjhhSooKJAkVVZWqrCw0H1dZWWlhgwZctBjxmIx7isB/BPNzc265ZZbJEnXX38975WAIg/gNykvUOrq6hQKtRyYCYfDSiaTkqTi4mIVFBRo9erVbkFSXV2tdevW6Qc/+EGqwwF8LxqN6r//+7/dNoKJPIDfpLxAOffcc3X99derZ8+eGjBggDZv3qybbrpJkydPliQ5jqPp06dr7ty56tu3r4qLizVnzhwVFRVp7NixqQ4H8L2MjAxdf/31tsOAZeQB/CblBcqtt96qOXPm6Ic//KH27dunoqIi/cd//IeuueYad58rrrhCtbW1mjp1qvbv368RI0Zo1apVyszMTHU4AADAgxzjwS/OV1dXKzc3V1VVVcrJybEdDmCVMcZdObRdu3YpWeI8FfeSOZKSjQ3affN3JUk9ZjymUIb3/9g53HvxtEUeAKl2OJ/fKf8WD4Ajq66uTu3bt1f79u1Z4jzAyAP4DQUKAABIO9zNGPC4du3aqaamxm0jmMgD+A0FCuBxjuMoOzvbdhiwjDyA33CJBwAApB0KFMDjGhsbddVVV+mqq65SY2Oj7XBgCXkAv6FAATyuqalJN9xwg2644QY1NTXZDgeWkAfwG+agAB4XiUQ0bdo0t41gIg/gN2Qx4HGxWEwLFy60HQYsIw/gN1ziAQAAaYcCBQAApB0KFMDjamtr5TiOHMdRbW2t7XBgCXkAv6FAAQAAaYdJsoDHtWvXTvv27XPbCCbyAH5DgQJ4nOM46tq1q+0wYBl5AL/hEg8AAEg7jKAAHtfY2Kif/exnkqTLL79cGRkZliOCDeQB/IYCBfC4pqYmXX311ZKk6dOn88EUUOQB/IYCBfC4SCSiSy65xG0jmMgD+A1ZDHhcLBbTnXfeaTsMpFjvWU8f/os6j5UkHXvti/pg/ujUBgQcYUySBQAAaYcCBQAApB0KFMDjamtrlZ2drezsbJY4D7BkY4N23TROu24ap2Rjg+1wgFZjDgrgA3V1dbZDQBowTXHbIQApQ4ECeFxWVpZ27NjhthFMTjRD3f7f3W4b8DoKFMDjQqGQevfubTsMWOY4IUVy822HAaQMc1AAAEDaYQQF8LimpiYtXrxYklReXq5oNGo5IthgEs369I0Da6d0OJE1UOB9FCiAxzU2NmrGjBmSpClTplCgBJRJNOuTlw4s2Nd+8CjL0QCtR4ECeFw4HNb3vvc9t41gckIhtTvuNLcNeB0FCuBxmZmZeuCBB2yHAcucSIa6nnu57TCAlKHMBgAAaYcCBQAApB0KFMDjamtr1bVrV3Xt2pWl7gMs2dig3Yu+p92LvsdS9/AF5qAAPvDRRx/ZDgFpIFlfbTsEIGUoUACPy8rK0rvvvuu2EUxONEOFkxe7bcDrKFAAjwuFQhowYIDtMGCZ44SU0bWX7TCAlGEOCgAASDuMoAAe19TUpHvuuUeSdNFFF7GSbECZRLNq3nlRktR+UKnlaIDWo0ABPK6xsVFTp06VJH3ve9+jQAkok2jWx8/dJknKPm6k3WCAFKBAATwuHA5rzJgxbhvB5IRCyup7itsGvI4CBfC4zMxMrVy50nYYsMyJZCjvvKtthwGkDGU2AABIOxQoAAAg7VCgAB5XV1en3r17q3fv3qqrq7MdDixJNjXoL0sm6y9LJivZxFL38D7moAAeZ4zRzp073TYCykiJ6n1uG/A6ChTA4zIzM7V+/Xq3LUm9Zz1tMyRY4ESiKvj3m9w24HUUKIDHhcNhDR061HYYsMwJhRUrPMZ2GEDKMAcFAACkHUZQAI9rbm7Www8/LEmaMGGCIhHe1kFkkgnVvveyJCm7/zcsRwO0Hv+TAR4Xj8d1wQUXSJLGjh1LgRJQprlJf3vqF5Kkdn1LLEcDtB7/kwEeFwqFVFpa6rYRUI6jzF5D3DbgdRQogMdlZWXphRdesB0GLAtFY8qfONd2GEDK8OcWAABIOxQoAAAg7bRJgbJnzx5dcMEF6ty5s7KysjRo0CBt3LjRfd4Yo2uuuUaFhYXKyspSaWmptm/f3hahAL5XV1enAQMGaMCAASx1H2DJpgbtveuH2nvXD1nqHr6Q8gLlk08+0fDhwxWNRvXss89qy5Yt+sUvfqGjjjrK3WfBggVatGiRli5dqnXr1ik7O1ujRo1SQwNvKuBwGWO0ZcsWbdmyhaXug8xITX/bpaa/7WKpe/hCyifJ3njjjerRo4eWLVvmbisuLnbbxhgtXLhQV199tcaMGSNJuvfee5Wfn6+VK1dq4sSJqQ4J8LXMzEz97ne/c9sIJicSVf7/vcFtA16X8hGUJ598UieddJLGjx+vvLw8nXDCCbrzzjvd53fs2KGKigr3a5GSlJubq2HDhmnt2rUHPWY8Hld1dXWLB4ADwuGwRo4cqZEjRyocDtsOB5Y4obAyex6vzJ7HywmRB/C+lBcof/7zn7VkyRL17dtXzz33nH7wgx/oRz/6kZYvXy5JqqiokCTl5+e3eF1+fr773OfNmzdPubm57qNHjx6pDhsAAKSRlBcoyWRSJ554om644QadcMIJmjp1qqZMmaKlS5d+5WPOnj1bVVVV7mP37t0pjBjwtubmZq1cuVIrV65Uc3Oz7XBgiUkmVLdtreq2rZVJJmyHA7RayuegFBYW6rjjjmuxrX///vr1r38tSSooKJAkVVZWqrCw0N2nsrJSQ4YMOegxY7GYYrFYqkMFfCEej+s73/mOJKmmpoal7gPKNDfpw8evlyT1mPGY5WiA1kv5CMrw4cO1devWFtu2bdumXr16STowYbagoECrV692n6+urta6detUUsL9I4DDFQqFdOqpp+rUU09lqfsgcxzFuvVXrFt/lrqHL6T8T60ZM2bo1FNP1Q033KDzzz9f69ev1x133KE77rhDkuQ4jqZPn665c+eqb9++Ki4u1pw5c1RUVKSxY8emOhzA97KysvTqq6/aDgOWhaIxFVzwM9thACmT8gJl6NChevzxxzV79mxdd911Ki4u1sKFC1VWVubuc8UVV6i2tlZTp07V/v37NWLECK1atYqvSAIAAEltdLPAc845R+ecc84/fd5xHF133XW67rrr2uLXAwAAj+OCNeBx9fX1Gjp0qIYOHar6+nrb4cCSZFNcf10+Q39dPkPJprjtcIBWY7o/4HHJZNK911UymbQcDawxRo0V29024HUUKIDHxWIxPfXUU24bweREour63WvdNuB1FCiAx0UiEY0ePdp2GLDMCYXV7uihtsMAUoY5KAAAIO0wggJ4XCKR0EsvvSRJOv3007lhYECZZEINO9+WJGX2Ot5yNEDrUaAAHtfQ0KCzzjpL0oGl7rOzsy1HBBtMc5P2PTJHEkvdwx8oUACPC4VCGjx4sNtGQDmOonnFbhvwOgoUwOOysrL05ptv2g4DloWiMRVNutV2GEDK8OcWAABIOxQoAAAg7VCgAB5XX1+vkSNHauTIkSx1H2DJprgqVsxSxYpZLHUPX2AOCuBxyWRSa9ascdsIKGMU3/2u2wa8jgIF8LhYLKZHHnnEbSOYnEhUXcbMctuA11GgAB4XiUQ0fvx422HAMicUVna/EbbDAFKGOSgAACDtMIICeFwikdDrr78uSTrllFNY6j6gTDKh+N6tkqRY0bGWowFajwIF8LiGhgaNGHFgaJ+l7oPLNDep8oErJLHUPfyBAgXwOMdx9LWvfc1tI6AcKXJUodsGvI4CBfC4du3aafv27bbDgGWhaKa6Tb3T/bn3rKdbdbwP5o9ubUhAqzBJFgAApB0KFAAAkHYoUACPa2ho0OjRozV69Gg1NDTYDgeWmOZG7Xv0x9r36I9lmhtthwO0GnNQAI9LJBJ65pln3DaCySSTqv/zRrfNPFl4HQUK4HEZGRlatmyZ20YwOeGIOn97utsGvI4sBjwuGo3qoosush0GLHPCEbUfVGo7DCBlmIMCAADSDiMogMclEgm98847kqRBgwax1H1AmWRCTR/ulCRFu/aSEyIP4G0UKIDHNTQ06IQTTpDEUvdBZpqb9Nd7fiTpwFL3TgYFCryNAgXwOMdxVFRU5LYRUI4Ubt/JbQNeR4ECeFy7du20Z88e22HAslA0U93L77UdBpAyTJIFAABphwIFAACkHQoUwOMaGho0fvx4jR8/nqXuA8w0N+rDlfP04cp5LHUPX6BAATwukUjoscce02OPPcZS9wFmkknVbX1VdVtflUkmbYcDtBqTZAGPy8jI0G233ea2EUxOOKJOZ/4/tw14HVkMeFw0GlV5ebntMGCZE46ow4nn2A4DSBku8QAAgLTDCArgcclkUn/6058kSUcffbRCIf7uCCJjkmr+5K+SpMhRhXIc8gDeRoECeFx9fb2OOeYYSSx1H2SmqVF77/wPSZ8tdZ9pOSKgdShQAB/Izc21HQLSgBOjOIV/UKAAHpedna39+/fbDgOWhTIy1XP6w7bDAFKGi5QAACDtUKAAAIC0Q4ECeFw8HtdFF12kiy66SPF43HY4sMQ0N+mjp2/WR0/fLNPcZDscoNUoUACPa25u1vLly7V8+XI1NzfbDgeWmGRCte+uVu27q2WS3PIA3sckWcDjotGoFixY4LYRTE44rI4jJ7ltwOsoUACPy8jI0OWXX247DFjmhKPKHTbOdhhAynCJBwAApB1GUACPSyaT+utfDyxxXlhYyFL3AWVMUomajyVJ4fadWOoenkeBAnhcfX29unfvLoml7oPMNDVqzy8vksRS9/AHChTAByIR3sqQFGJyLPyD/9WANNR71tOHtX+3/1wpSRrw09+nPBZ4QygjU70uf8J2GEDKcJESAACkHQoUAACQdtq8QJk/f74cx9H06dPdbQ0NDSovL1fnzp3Vvn17jRs3TpWVlW0dCuBLprlJf3t+if72/BKWOA8w8gB+06YFyoYNG3T77bfr+OOPb7F9xowZ+u1vf6tHH31Ua9as0d69e3Xeeee1ZSiAb5lkQjWbn1bN5qdZ4jzAyAP4TZtNkq2pqVFZWZnuvPNOzZ07191eVVWlu+++WytWrNDpp58uSVq2bJn69++v119/XaecckpbhQT4khMOK3f4/3XbCCbyAH7TZiMo5eXlGj16tEpLS1ts37Rpk5qamlps79evn3r27Km1a9ce9FjxeFzV1dUtHgAOcMJRdRxRpo4jyuSEuRdPUJEH8Js2GUF56KGH9MYbb2jDhg1feK6iokIZGRnq2LFji+35+fmqqKg46PHmzZunn/zkJ20RKgAASEMpH0HZvXu3pk2bpgceeECZmalZyXD27NmqqqpyH7t3707JcQE/MMYo2VCjZEONjDG2w4El5AH8JuUFyqZNm7Rv3z6deOKJikQiikQiWrNmjRYtWqRIJKL8/Hw1NjZq//79LV5XWVmpgoKCgx4zFospJyenxQPAAaYprt23TNTuWybKNMVthwNLyAP4Tcov8Zxxxhl65513WmybNGmS+vXrpyuvvFI9evRQNBrV6tWrNW7cgVuDb926Vbt27VJJSUmqwwEAAB6U8gKlQ4cOGjhwYItt2dnZ6ty5s7v94osv1syZM9WpUyfl5OTosssuU0lJCd/gAb4CJxpTz/9aeeAH7sUSWOQB/MbKvXhuvvlmhUIhjRs3TvF4XKNGjdIvf/lLG6EAnuc4jhTmtlpBRx7AbxzjwdlU1dXVys3NVVVVFfNRkJYO92Z/aJ1kY4N23/xdSVKPGY8plJGaCfr46j6YP9p2CEhDh/P5TbkNeJxJNGn/y/dJkjp+4/usgRFQ5AH8hpsFAh5nEglVr/+Nqtf/RibBEudBRR7AbxhBATzOCYeVc/J5bhvBRB7AbyhQAI9zwlEd9c3JtsOAZeQB/IZLPAAAIO0wggJ4nDFGSv7vnINQ+MDXTRE45AH8hhEUwONMU1y7fj5Wu34+liXOA4w8gN9QoAAAgLTDJR7A45xoTD2mPeS2EUzkAfyGAgXwOMdx5GS2tx0GLCMP4Ddc4gEAAGmHERTA40yiSVVrH5Ek5ZaczxLnAUUewG8oUACPM4mEql59UJKUc/I4PpgCijyA31CgAB7nhMJqf8Jot41gIg/gNxQogMc5kag6n/UD22HAMvIAfsMkWQAAkHYoUAAAQNqhQAE8LtnYoJ0/G6OdPxujZGOD7XBgCXkAv2EOCuAHn90kDsFGHsBHKFAAj3OiGer2w3vcNoKJPIDfUKAAHuc4IUU6dLEdBiwjD+A3zEEBAABphxEUwONMoknVG5+UJOWc9H9YQTSgyAP4DQUK4HEmkdD+3y+TJHU4YTQfTAFFHsBvKFAAj3NCYWUPPMNtI5jIA/gNBQrgcU4kqi6jZ9gOA5aRB/AbJskCAIC0Q4ECAADSDgUK4HHJxgbtWjhBuxZOYInzACMP4DfMQQF8wMRrbYeANEAewE8oUACPc6IZKppyu9tGMJEH8BsKFMDjHCekaKdutsOAZeQB/IYCBb7Te9bTtkMAALQSBQrgcSbRrJq3VkmS2g/+lpwwb+sgIg/gN2Qw4HEm0ayPX1gqScoeWMoHU0CRB/AbMhjwOCcUUrtjh7ttBBN5AL+hQAE8zolkqOvY2bbDgGXkAfyGMhsAAKQdChQAAJB2KFAAj0s2Negvi/9df1n870o2scR5UJEH8BvmoABeZ6REzcduGwFFHsBnKFAAj3MiURVetMhtI5jIA/gNBQrgcU4orIz8PrbDgGXkAfyGOSgAACDtMIICeJxJNKt2y+8lSdnHjWQF0YAiD+A3ZDDgcSbRrL89s1CS1O7YEXwwBRR5AL8hgwGPc0IhZfU5yW0jmMgD+A0FCuBxTiRDeeN/bDsMWEYewG8oswEAQNqhQAEAAGmHAgXwuGRTg/bcMUV77pjCEucBRh7Ab5iDAnidkZo/+avbRkCRB/AZChTA45xIVPllC9w2gok8gN9QoAAe54TCyux+nO0wYBl5AL9J+RyUefPmaejQoerQoYPy8vI0duxYbd26tcU+DQ0NKi8vV+fOndW+fXuNGzdOlZWVqQ4FAAB4VMoLlDVr1qi8vFyvv/66XnjhBTU1Nemss85SbW2tu8+MGTP029/+Vo8++qjWrFmjvXv36rzzzkt1KEAgmGRCtX98RbV/fEUmmbAdDiwhD+A3Kb/Es2rVqhY/33PPPcrLy9OmTZv0jW98Q1VVVbr77ru1YsUKnX766ZKkZcuWqX///nr99dd1yimnpDokwNdMc5M+emK+JKnHjMfkZIQtRwQbyAP4TZvPQamqqpIkderUSZK0adMmNTU1qbS01N2nX79+6tmzp9auXXvQAiUejysej7s/V1dXt3HUgIc4jmI9BrptBBR5AJ9p0wIlmUxq+vTpGj58uAYOPPDGqaioUEZGhjp27Nhi3/z8fFVUVBz0OPPmzdNPfvKTtgwV8KxQNKaC7823HQYsIw/gN226UFt5ebneffddPfTQQ606zuzZs1VVVeU+du/enaIIAQBAOmqzEZRLL71UTz31lF5++WV1797d3V5QUKDGxkbt37+/xShKZWWlCgoKDnqsWCymWCzWVqECAIA0k/IRFGOMLr30Uj3++ON66aWXVFxc3OL5r3/964pGo1q9erW7bevWrdq1a5dKSkpSHQ7ge8mmuPYuu0x7l12mZFP8y18AXyIP4DcpH0EpLy/XihUr9MQTT6hDhw7uvJLc3FxlZWUpNzdXF198sWbOnKlOnTopJydHl112mUpKSvgGD/BVGKOmfTvcNgKKPIDPpLxAWbJkiSRp5MiRLbYvW7ZMF110kSTp5ptvVigU0rhx4xSPxzVq1Cj98pe/THUoQCA4kajyzv+p20YwkQfwm5QXKOYQKvfMzEwtXrxYixcvTvWvBwLHCYWVVXyC7TBgGXkAv2nTb/EAAAB8FdwsEPA4k0yofscbkqSs4hPlhFhBNIjIA/gNBQrgcaa5SR8+dmAhQ5Y4Dy7yAH5DgQJ4neMoo6Cv20ZAkQfwGQoUwONC0ZgKL7zZdhiwjDyA3zBJFgAApB0KFAAAkHa4xAN4XLIprn0PXy1JypswV6Eo960KonTLg96znrb6+z+YP9rq70frUaAAXmeM4nvec9sIKPIAPkOBAnicE4mq63euctsIJvIAfkOBAnicEwqr3THcCTzoyAP4DZNkAQBA2mEEBfA4k0wo/pc/SJJi3QewxHlAkQfwGwoUwONMc5MqH/xvSSxxHmTkAfyGAgXwOkeKdu7pthFQ5AF8hgIF8LhQNFNFl/zSdhiwjDyA3zBJFgAApB0KFAAAkHa4xAN4XLIprg9//VNJUtdxc6wvcQ47yAP4DQUK4HXGqGHnm24bAUUewGcoUACPcyJRdT7nP902gok8gN9QoAAe54TCaj/gm7bDgGXkAfyGSbIAACDtMIICeJxJJtRY+SdJUkb+0SxxHlDkAfyGAgXwONPcpIp7Z0piifMgIw/gNxQogNc5Ujgnz20joMgD+AwFCuBxoWimuv/gV7bDgGXkAfyGSbIAACDtUKAAAIC0wyUewONMc6M+fHKBJKnr/7lCTiTDckSwgTyA31CgAB5nkknVb3/dbTM/MpjIA/gNBQrgcU44ok6jLnXbCCbyAH5DFgMe54Qj6jDkW7bDgGXkAfyGAgVpp/esp22HAACwjAIF8Dhjkmr6aLckKdqlhxyHL+cFEXkAv6FAATzONDXqr78ql/TZEueZliOCDeQB/IYCBfCBUFaO7RCQBsgD+AkFCuBxoYxM9fjRCtthwDLyAH7DRUoAAJB2KFAAAEDa4RIP4HGmuVEfPXuLJKnL2dNY4jygyIOWbC9X8MH80VZ/vx8wggJ4nEkmVbdljeq2rJFJJm2HA0vIA/gNIyiAxznhiI46fYrbRjCRB/AbshjwOCccUc7QMbbDgGXkAfyGSzwAACDtMIICeJwxSSWqP5QkhXO6ssR5QJEH8BsKFKSc7dnzQWOaGrVn6cWSWOI8yMgD+A0FCuADTjRmOwSkAfIAfkKBAnhcKCNTPWf+2nYYsIw8gN9wkRIAAKQdChQAAJB2uMQDeJxpbtLHLyyRJHU68wdyIlHLEcEG8gB+wwgK4HEmmVDN28+r5u3nZZIJ2+HAEvIAfsMICr6Arwl7ixMOq+O/fd9tI5jIA/gNBQrgcU44qtxTJ9gOA5aRB/Abq5d4Fi9erN69eyszM1PDhg3T+vXrbYYDAADShLURlIcfflgzZ87U0qVLNWzYMC1cuFCjRo3S1q1blZeXZyssSa2/xPHB/NEpigT4csYYJeurJUmhrBw5jmM5IthAHqQX258jqbhUb/uzzNoIyk033aQpU6Zo0qRJOu6447R06VK1a9dOv/rVr2yFBHiSaYrrL7eW6S+3lsk0xW2HA0vIA/iNlRGUxsZGbdq0SbNnz3a3hUIhlZaWau3atV/YPx6PKx7/+xuuqqpKklRdXd0m8SXjda16fVvFdaS09vxxZCUbG/7ejtdJJmkxGjvoA/rAb1r7OZKK/8fb4rPss2MaY758Z2PBnj17jCTz2muvtdh++eWXm5NPPvkL+1977bVGEg8ePHjw4MHDB4/du3d/aa3giW/xzJ49WzNnznR/TiaT+vjjjxWNRtWzZ0/t3r1bOTk5FiO0p7q6Wj169KAP6AP6gD6gD0QfSOndB8YYffrppyoqKvrSfa0UKF26dFE4HFZlZWWL7ZWVlSooKPjC/rFYTLFYy7t0duzY0R0qysnJSbt/hCONPqAPJPpAog8k+kCiD6T07YPc3NxD2s/KJNmMjAx9/etf1+rVq91tyWRSq1evVklJiY2QAABAGrF2iWfmzJm68MILddJJJ+nkk0/WwoULVVtbq0mTJtkKCQAApAlrBcqECRP04Ycf6pprrlFFRYWGDBmiVatWKT8//5CPEYvFdO21137h8k+Q0Af0gUQfSPSBRB9I9IHknz5wjDmU7/oAAAAcOdzNGAAApB0KFAAAkHYoUAAAQNqhQAEAAGnHSoHy8ssv69xzz1VRUZEcx9HKlSu/9DWLFy9W//79lZWVpWOPPVb33nvvF/bZv3+/ysvLVVhYqFgspmOOOUbPPPOM+/ySJUt0/PHHu4vXlJSU6Nlnn03lqR0yW33wj+bPny/HcTR9+vRWns1XY6sPfvzjH8txnBaPfv36pfLUDpnNPNizZ48uuOACde7cWVlZWRo0aJA2btyYqlM7LLb6oXfv3l/IBcdxVF5ensrTOyS2+iCRSGjOnDkqLi5WVlaWjj76aP30pz89tHulpJitPvj00081ffp09erVS1lZWTr11FO1YcOGVJ7aIWuLPhg5cuRB83z06L/frdgYo2uuuUaFhYXKyspSaWmptm/fnurTOyxWvmZcW1urwYMHa/LkyTrvvPO+dP8lS5Zo9uzZuvPOOzV06FCtX79eU6ZM0VFHHaVzzz1X0oEbEJ555pnKy8vTY489pm7dumnnzp3q2LGje5zu3btr/vz56tu3r4wxWr58ucaMGaPNmzdrwIABbXW6B2WrDz6zYcMG3X777Tr++ONTfWqHzGYfDBgwQC+++KL7cyRi5xv3tvrgk08+0fDhw/XNb35Tzz77rLp27art27frqKOOaqtT/Zds9cOGDRuUSCTcn999912deeaZGj9+fMrP8cvY6oMbb7xRS5Ys0fLlyzVgwABt3LhRkyZNUm5urn70ox+11ekelK0+uOSSS/Tuu+/qvvvuU1FRke6//36VlpZqy5Yt6tatW1ud7kG1RR/85je/UWNjo/uav/3tbxo8eHCLPF+wYIEWLVqk5cuXq7i4WHPmzNGoUaO0ZcsWZWZmpv5ED0UK7v3XKpLM448//i/3KSkpMf/1X//VYtvMmTPN8OHD3Z+XLFli+vTpYxobGw/r9x911FHmrrvuOqzXpNqR7oNPP/3U9O3b17zwwgvmtNNOM9OmTfuqoafMkeyDa6+91gwePLg14baJI9kHV155pRkxYkSr4m0rNv9PmDZtmjn66KNNMpk8rJhT7Uj2wejRo83kyZNbbDvvvPNMWVnZ4QeeQkeqD+rq6kw4HDZPPfVUi+0nnniiueqqq75a8CmSqj74vJtvvtl06NDB1NTUGGOMSSaTpqCgwPzsZz9z99m/f7+JxWLmwQcf/Oon0EqemIMSj8e/UMFlZWVp/fr1ampqkiQ9+eSTKikpUXl5ufLz8zVw4EDdcMMNLf46+keJREIPPfSQamtrPbG8fir7oLy8XKNHj1ZpaekRiz8VUtkH27dvV1FRkfr06aOysjLt2rXriJ1Ha6SqD5588kmddNJJGj9+vPLy8nTCCSfozjvvPKLn0hpt8X9CY2Oj7r//fk2ePFmO47T5ObRWqvrg1FNP1erVq7Vt2zZJ0ltvvaVXXnlFZ5999pE7ma8oFX3Q3NysRCJx0OO88sorR+ZEWuFQ+uDz7r77bk2cOFHZ2dmSpB07dqiioqLFZ0Jubq6GDRumtWvXtl3wX8ZaafS/dAgV4uzZs01BQYHZuHGjSSaTZsOGDSY/P99IMnv37jXGGHPssceaWCxmJk+ebDZu3Ggeeugh06lTJ/PjH/+4xbHefvttk52dbcLhsMnNzTVPP/10W53aITuSffDggw+agQMHmvr6emOM8dQISqr64JlnnjGPPPKIeeutt8yqVatMSUmJ6dmzp6murm7LU/xSR7IPYrGYicViZvbs2eaNN94wt99+u8nMzDT33HNPW57iITnS/yd85uGHHzbhcNjs2bMn1ad02I5kHyQSCXPllVcax3FMJBIxjuOYG264oS1P75AcyT4oKSkxp512mtmzZ49pbm429913nwmFQuaYY45py1P8Uqnqg3+0bt06I8msW7fO3fbqq68edP/x48eb888/PyXn8lV4okCpq6szkyZNMpFIxITDYVNUVGSuuOIKI8lUVFQYY4zp27ev6dGjh2lubnZf94tf/MIUFBS0OFY8Hjfbt283GzduNLNmzTJdunQxf/jDH1J+XofjSPXBrl27TF5ennnrrbfc571UoKQyD/7RJ598YnJycjxxqS9VfRCNRk1JSUmLY1922WXmlFNOSd0JfUW2cuGss84y55xzTsrOozWOZB88+OCDpnv37ubBBx80b7/9trn33ntNp06drBerR7IP3n//ffONb3zDSDLhcNgMHTrUlJWVmX79+rXJuR2qVPXBP5o6daoZNGhQi23pWqB44hJPVlaWfvWrX6murk4ffPCBdu3apd69e6tDhw7q2rWrJKmwsFDHHHOMwuGw+7r+/furoqKixeSgjIwMfe1rX9PXv/51zZs3T4MHD9Ytt9xyxM/pcKWiDzZt2qR9+/bpxBNPVCQSUSQS0Zo1a7Ro0SJFIpF/OvSdLlKZB/+oY8eOOuaYY/T+++8fkfNojVT1QWFhoY477rgWx+7fv79nLnWlOhd27typF198UZdccskRPY/WSFUfXH755Zo1a5YmTpyoQYMG6fvf/75mzJihefPmWTmvw5GqPjj66KO1Zs0a1dTUaPfu3e7lkT59+lg5r8NxKH3wmdraWj300EO6+OKLW2wvKCiQJFVWVrbYXllZ6T5ngycKlM9Eo1F1795d4XBYDz30kM455xyFQgdOYfjw4Xr//feVTCbd/bdt26bCwkJlZGT802Mmk0nF4/E2jz1VWtMHZ5xxht555x29+eab7uOkk05SWVmZ3nzzzRZv4HSW6jyoqanRn/70JxUWFh6R+FOhtX0wfPhwbd26tcUxt23bpl69eh25k0iBVOXCsmXLlJeX1+Jrl17R2j6oq6tz9/9MOBxu8Zp0l6o8yM7OVmFhoT755BM999xzGjNmzBE9j9b4V33wmUcffVTxeFwXXHBBi+3FxcUqKCjQ6tWr3W3V1dVat26d3TmaNoZtPv30U7N582azefNmI8ncdNNNZvPmzWbnzp3GGGNmzZplvv/977v7b9261dx3331m27ZtZt26dWbChAmmU6dOZseOHe4+u3btMh06dDCXXnqp2bp1q3nqqadMXl6emTt3rrvPrFmzzJo1a8yOHTvM22+/bWbNmmUcxzHPP//8ETv3z9jqg8+zeYnHVh/853/+p/n9739vduzYYV599VVTWlpqunTpYvbt23fEzv0ztvpg/fr1JhKJmOuvv95s377dPPDAA6Zdu3bm/vvvP2Ln/o9svh8SiYTp2bOnufLKK4/Iuf4ztvrgwgsvNN26dTNPPfWU2bFjh/nNb35junTpYq644oojdu6fsdUHq1atMs8++6z585//bJ5//nkzePBgM2zYsMP+VmgqtEUffGbEiBFmwoQJB/298+fPNx07djRPPPGEefvtt82YMWNMcXGxO1/RBisFyu9+9zsj6QuPCy+80Bhz4A1z2mmnuftv2bLFDBkyxGRlZZmcnBwzZswY88c//vELx33ttdfMsGHDTCwWM3369DHXX399i+uOkydPNr169TIZGRmma9eu5owzzrBSnBhjrw8+z2aBYqsPJkyYYAoLC01GRobp1q2bmTBhgnn//ffb+nQPymYe/Pa3vzUDBw40sVjM9OvXz9xxxx1tear/ks1+eO6554wks3Xr1rY8xS9lqw+qq6vNtGnTTM+ePU1mZqbp06ePueqqq0w8Hm/rU/4CW33w8MMPmz59+piMjAxTUFBgysvLzf79+9v6dA+qrfrgj3/8o5H0Tz/zksmkmTNnjsnPzzexWMycccYZ1t8TjjEWlgsEAAD4Fzw1BwUAAAQDBQoAAEg7FCgAACDtUKAAAIC0Q4ECAADSDgUKAABIOxQoAAAg7VCgAACAtEOBAgAA0g4FCgAASDsUKAAAIO1QoAAAgLTz/wEjVeHktUjDMAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "centroids = np.array([centroid(x, make_y(x)) for _ in range(1000)])\n",
    "plt.hist(centroids, bins='auto')\n",
    "plt.axvline(line.mean.value, c='k')\n",
    "plt.axvline(line.mean.value + line.stddev.value, c='k', ls=':')\n",
    "plt.axvline(line.mean.value - line.stddev.value, c='k', ls=':')\n",
    "np.mean(centroids), np.std(centroids)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb1d4dd8-1ac4-42bc-963c-2f4c2303feee",
   "metadata": {},
   "source": [
    "Huzzah, that is also consistent!"
   ]
  }
 ],
 "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
}