derivation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 387,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from uncertainties import ufloat,umath\n",
    "from matplotlib.patches import Ellipse\n",
    "\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Estimation of Errors in Ellipse Fitting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ellipse Fitting\n",
    "An ellipse is a special case of a general conic section described by\n",
    "\n",
    "$$F(x,y) = ax^2 + bxy + cy^2 + dx + ey + f = 0$$\n",
    "\n",
    "with the constraint $b^2 - 4ac < 0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 388,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generating a noisy ellipse\n",
    "def generate(a,b,center=(0,0),n=100,phi=0,weight=0.2):\n",
    "    t = np.linspace(0,2*np.pi,n)[:-1]\n",
    "    x = a*np.cos(t)*np.cos(phi) - b*np.sin(t)*np.sin(phi) + center[0]\n",
    "    y = a*np.cos(t)*np.sin(phi) + b*np.sin(t)*np.cos(phi) + center[1]\n",
    "    x += np.random.randn(len(t))*weight\n",
    "    y += np.random.randn(len(t))*weight\n",
    "    return x,y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 389,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.607823991300505 1.4667002536011795 [2.7913013  2.73931612] 1.9007356427967286\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEGCAYAAABrQF4qAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAf8ElEQVR4nO3df4wc9Znn8ffD2ME/sM3GdlgL+84oG+UUNhiHMSRL9pQBccLYm9hJduWTMreRWA341msibmVlTjeOMJezbHZMYmWP21wSbQJJHJRgxJL49pCm90JyB3gMNsFmjwWHHINATGDjQMABM8/9Ud1MT7u7p39UdX2r6vOSSv2rpvpb011Pf+up7w9zd0REJL/OSbsAIiKSLAV6EZGcU6AXEck5BXoRkZxToBcRybk5aRegXcuWLfPVq1enXQwRkaAcOXLkl+6+vN5rmQv0q1evZnx8PO1iiIgExcx+0eg1pW5ERHJOgV5EJOcU6EVEci5zOXoRkVpvvfUWExMTnD59Ou2iJG7evHmsXLmSuXPntvw3CvQiknkTExMsWrSI1atXY2ZpFycx7s7LL7/MxMQEF110Uct/p9SNiJxt714olWY+VypFzwfo9OnTLF26NNdBHsDMWLp0adtnLgr0InK2devgT/5kOtiXStHjdevSLVcTeQ/yFZ3sp1I3InK2gQG4++4ouG/dCnfcET0eGEi7ZNIB1ehFpL6BgSjI33prdKsg31RfXx+XXnopF198MWvWrGF0dJSpqammf/Pss8/yne98J/GyKdCLSH2lUlSTHxmJbmtz9nFI41pAQu85f/58jh49yvHjx3nggQc4dOgQt9xyS9O/6VWgx90ztVx22WUuUmh79riPjc18bmwsej4uY2Puy5ZNv0/t48De58SJEz1/z1oLFy6c8fiZZ57xd7/73T41NeU///nP/aMf/aivXbvW165d6z/96U/d3f2KK67wxYsX+5o1a3zfvn0N16tVb3+BcW8QNxMPzEAf8Bhwf53XPgtMAkfLy5/Ntj0Feim8XgThXvyYVG932TL3kZHm+9GkTG0F+nbesw21gd7dfcmSJf7iiy/6b37zG3/jjTfc3f2pp57yShwrlUq+YcOGd9ZvtF6tdgN9Ly7G3gQ8CSxu8Pr33H1bD8ohkg+9uFC6Y0f9900iT199LWBkpPF7VFoCVfa10hLo7ruTe8+YvPXWW2zbto2jR4/S19fHU0891dV67Uo0R29mK4ENwNeSfB+RrLX77lqeLpS2ei2g+gdu586ZQT+p9+zCyZMn6evr4z3veQ+33347F1xwAceOHWN8fJw333yz7t+0ul67kr4Y+yVgB9Ds0vOnzOxxM/u+ma2qt4KZDZnZuJmNT05OJlFOyboMtvvuSi8ulPZCda18167pQN4s2Hf7A9fue3ZgcnKSG2+8kW3btmFmnDp1ihUrVnDOOedw55138vbbbwOwaNEiXn311Xf+rtF6XWuU0+l2ATYC/7V8/2PUz9EvBc4t378BGJttu8rRS0MJ5F1b0st8dmXbSeToe70fnbxng8+4rRx9Qvt5zjnn+Jo1a/wDH/iAX3LJJX7bbbf522+/7e5Rvv2DH/ygX3LJJb5jx4538vlvvvmmDwwM+CWXXOL79u1ruF6tYC7GAruBCeBZ4EXgdeCuJuv3Aadm264CfYF0ckCOjERf65GRZMtWW6ZetFCpSCogJ7kfcZS5SfnavhibccEE+hlv0rhGv6Lq/mbgodm2pUBfIO0GnrRq9Gm/d5yS2o84fkTibHWTccEHemAX8HGfrvUfB44BJeBfzbYtBfqCaTXw9LpWXU8aZxNJaGc/2qmpJ/hjqEAfQKCPc1GgL6BWAk8a+eXa9+omiKVd/ur3bGc/2v2BTejH8MSJEz41NRXrNkM1NTWlQC85k4WUSBxnEyGckXRahnbPuhL4LE+ePOmTk5O5D/ZTU1M+OTnpJ0+ePOu1ZoFeo1dKuKqbwVU663TTdjophw/PLFOlvffhw9H9vXujZp7VZS6VotcrHZNCGC1ytv1opJXORwl/litXrmRiYoJ3ml+fOgXnngvz5k2vdPo0/Pa3sGRJ1++XpsoMU21p9AsQ6qIafYGEks7oVjs15Szm+VupqeelCWrAUOpGJGWtBMM00lTdBuCQA2oW0n4xUqAXCUGz2npaAbPb9w39rCuLZ0gdUqAXSdtstcs0A2Zea7553a8GFOglGaHX5kIRcnqjIm813yz8z2PWLNBrhinpXLsDiWVphMk4y9qsNUsI8jJAWrXbboPh4Zn/8+Hh6PkiavQLEOqiGn1g2jk9DrGW1eisZGgovLImIcTPpFXNziizvF8dQqkbSVQ7p/2h5U2bBYTQypqELKffqj+rPXvcR0dnfk6jo+4LF+b786uiQC/J6SQYhpYPbrYPoZVVZqp8doOD7mZRcK99viCfnwK9JKOT0+NQa8n1AnqoZZWZKp/d4ODMz6tSwy/I56dAL8nodNKI0PKm9QJ6qGWVmWo/u0oNvhL0C/T5KdBLGELMBzcK6END4ZVVZqr97EZHo/TN4KD7ggXTaZzq9XP8+TUL9Ba9nh39/f0+Pj6edjEkL1oZcEzCVP3ZVZr2Dg/DmTPTTX9DGwAvQWZ2xN37676mQC8imacfbAV6EZG8axboE+8Za2Z9ZvaYmd1f57Vzzex7Zva0mT1sZquTLo/kWJZ63or0UC+GQLgJeLLBa9cD/+zuvwfcDuzpQXkkr9odkkGkIBIN9Ga2EtgAfK3BKp8Avlm+/33gajOzJMskOVY9S9POnYW7GCfSSNI1+i8BO4CpBq9fCDwH4O5ngFPA0tqVzGzIzMbNbPydqcJE6qme1m7rVgV5aU9O03+JBXoz2wi85O5Hut2Wu3/V3fvdvX/58uUxlE66FuoBkceRGKV3cpr+S7JGfyXwcTN7FjgAXGVmd9Ws8zywCsDM5gBLgJcTLJPEJcQDonoC6l27ptM4CvbSqpym/xIL9O4+7O4r3X01sAUYc/fP1Kx2H/Cn5fufLq+TrfaeRRXiARH6uO+SDTlM/83p9Rua2S6irrr3AV8H7jSzp4FXiH4QJCuqD4iRkfQPiHodYwYG0i+XZEtt+i8H36GeBHp3/wfgH8r3d1Y9fxr4416UQRKQwwNCCq46/Vf5PodwttolTSUonVE+XPIop+k/BXrpTN4OiFBbEUlv7dhxds19YCDz4+Uo0Etn8nZAhNiKSCQmCvQiEGYrIglXxs4AFehFKnLYrE4SkrEzQAX6PMpYbSMY6lUrrerkDDDF41KBPo8yVtsIgloRSbvaPQNM87hsNMdgqIvmjG1RvQmvpbEQ57OVsHVyjCV4XKLJwQtqZCT6iEdG0i6JSG/06ge70aTyrQTuhI7LZoFeqZu8Ur5ZiqhX6ZFO+5FUjsurr4b9+2cel0nm6xv9AoS6FKpG32ntpJvahkjWhZi23LPHfXR0ujxjY+4LF7rPmzf9uMuyohp9RnVaO8lbr1WRdoTYTHbduujsenh4ujx9fdHyxS8m32+j0S9AqEuhavTuYdZORELWzTETd46/enuVcg0Oui9YED2OMV+PLsZmXCgXVdUyRULXbdoy7rRn7d8PDkbH8uBg7JU4BfosC6lGr9y/hC6Oykgcx1y9mvw110wH+cWL3ZcsifVYUqDPqhADa0g/PCJJ6fYsuvZYrQT5a66JHg8NRcG++vjp8uxYgT6rQk2VhJJKEklCXJWZ6px8JcjXVtxiPJYV6KWxdn9MVKOXPIv7LLo6Jx/H9ppoFugTa15pZvPM7BEzO2Zmx83sljrrfNbMJs3saHn5s6TKIw2004RT48FI3sXZNLlUgnvugcFBOHQoepxWU+dGvwDdLoAB55XvzwUeBj5cs85nga+0s13V6KvEldpptZYeaipJJDQpXF8jjRp9+b1fKz+cW148qfcrpLi6e7fawSRvs0qJJCWwTouJ9ow1sz4zOwq8BDzg7g/XWe1TZva4mX3fzFY12M6QmY2b2fjk5GSSRc6WuGZF0rg4Ip1pNMY8hFUpalTVj3MBzgdKwO/XPL8UOLd8/wZgbLZtKXVTRzetYEJswimSFQEdP6Q91o27/6oc6K+tef5ld/9t+eHXgMt6UZ5c6bY2HtgppkimZGWu4Ua/AN0uwHLg/PL9+cCDwMaadVZU3d8MPDTbdlWjrxJQbUKk0ALoW0JKNfoVQMnMHgcOE+Xo7zezXWb28fI628tNL48B24la4UirVBsXSUd1br5yVj04CKOjQV7jsuiHIDv6+/t9fHw87WKISJFVWrgND8Pu3WffppC+MbMj7t5f7zWNRx+CFGeHF5EOVM6ed+6E9eung/vNNwd5Vq1AH4I0Z4cXkc4MDESB/c47Z/Y/CbBviQJ9CLJy5V5EpmWo/4kCfShCnP5MROrL2LhPCvShyFDtQKTwMtbiTa1uQlBdOxgYOPuxiMgs1OomdBmrHYhItsxJuwCFsXdv1IqmuoZeKkXBvN4V+oEB1eZF8qTdGBAj1eh7pZsmlGpnLxKuVo/PNJtRNxobIdQl02PddDoNn8a0EQlXO8dnglNxojljA9Lp4Eeaq1UkXO0cnwkNgNYs0Ct1E7dmp3HdNKFUO3uRcLV6fKbVjLrRL0CoS/A1+kancaOj3aVfVKMXCVcrx2fCKViUuumxeh96NxNrK0cvEq5Wj89uYkALmgV6dZhKys6d0WncyEjURbobKTbLEpFZBHJ8NuswpUCfhEqzqa1bozzcbD1cA/miiORCQY8n9YztpU4GO9IwxSLxafV4KlD/FAX6uHUynIGGKRaJT6vHU5EqWI2S990uwDzgEeAYcBy4pc465wLfA54GHgZWz7bdTFyM7VQAEwyL5EYrx1OOWrORUjv63wJXufsa4FLgWjP7cM061wP/7O6/B9wO7EmwPGHTMMUi8Wn1eCpK/5RGvwBxLsAC4FHgiprn/x74SPn+HOCXlC8QN1pyWaNX80mR+AQyJEGvkVbPWDPrM7OjwEvAA+7+cM0qFwLPlX9wzgCngKV1tjNkZuNmNj45OZlkkdOhYYpF4tPq8ZSxWaK60ZPmlWZ2PnAQ+At3f6Lq+SeAa919ovz4GaJa/y8bbSsTzStFJHw5a4bZrHllT8ajd/dfmVkJuBZ4ouql54FVwISZzQGWAC/3okwiUnAFmgcisdSNmS0v1+Qxs/nANcA/1qx2H/Cn5fufBsa8F6cYIiIFkmSNfgXwTTPrI/pBudvd7zezXUQXDe4Dvg7caWZPA68AWxIsj4hIISUW6N39cWBtned3Vt0/DfxxUmUQERH1jBURSVcPhmJQoBcRSVMPhmLoSasbERFpoHpsnlZHvG2TavQiImlLeCgGBXoRkbQlPNaVAr2ISJp6MBTDrIHezP7CzH4ntnfMmgJNTiAiKejBWFet1OgvAA6b2d1mdq2ZWWzvngVFmpxARHpvx46zc/IDA7GOtzNroHf3/wS8j6gX62eBfzKz/2Jm742tFCHT7E8iknEt5ejL48+8WF7OAL8DfN/MipG/KMrkBCKSS63k6G8ysyPAXuCnwAfdfStwGfCphMsXBs3+JCIZ1kqHqXcDn3T3X1Q/6e5TZrYxmWIFpPqKeGUIU6VvRCRDWsnRf6E2yFe99mT8RQqMZn8SEYha2t1ww8wz+lIpei7wVng9mWEqTkHPMJWzGWtEpEqpBJs2gRkcPBg9V/045TP8ZjNMqcNUnNQUUyS/Bgbg3nvBHTZuhA0bggnys1Ggj5OaYork28AA3HQTvP46vPEGbN+eieNbgT5uaoopkl+lEnz5y7BgAcyfD/v3Z6IVXpJzxq4ys5KZnTCz42Z2U511PmZmp8zsaHnZWW9bmaKmmCL5VJ2jv/9++OEPozTO5s3BH+dJjkd/BvgP7v6omS0CjpjZA+5+oma9B909H8001RRTJL8OH4YtW6Klcjzfey8cOBC9FvAxnuScsS8AL5Tvv2pmTwIXArWBPjtma1XTrClmwF8CEWlBvZZzlQpd4HqSozez1UQThT9c5+WPmNkxMztkZhc3+PshMxs3s/HJyckki9rcbK1qejA4kYhIuxJvR29m5wH/C/iiu99T89piYMrdXzOz64Avu/v7mm0v9Xb0leCe0JRfIiKdSK0dvZnNBX4AfLs2yAO4+6/d/bXy/R8Bc81sWZJlmtVs48+rVY2IZEySrW6MaGjjJ919X4N1frcyvr2ZXV4uz8tJlakl9dIzGzfCnDnTj++4AwYHYXQ0+KvtIiJJtrq5EhgEfmZmR8vP/UfgXwC4+38DPg1sNbMzwBvAFk97TIZ6M7Lfeivs3h29vns3DA9Ht7feqlY1IhK8JFvd/ARoOhuVu38F+EpSZehYdXpmZARuvhnWroU/+iP45CejIF8J7mvXqlWNiAQtyRp9dtV2eqo0obr55ungX92EUkFeRAKmIRBqNZqRfd8+9XgVyYrZGlUUjAJ9rXqdnoaHo0HKaoO/gr1ImDSS7AxK3dSq17npzBn4u79Tj1eRrKjXqKLAjSY08UhSNAmJSPp27py+rrZrV9qlSZQmHkmDTh1F0qWRZN+hQA/JXLjRJCQi6WnUqKKgwV6BHpKrfWu4BJF0NBtJtoCUo69IYrAyDYAmIj2iHP1sKima6tp39fOd0KmjiARCgR6iFM3mzdFckCMj0e3mzd2lbnTqKCKBUDv6CvdoLkiIbrtNaWV4NhoRyRfV6CGqZd97L2zfHqVutm+PHqv2LZJPBRsiQYEepmvf1W1uq58XkXwpWD8XBXrQhVORPKtXe4doyPGC9HNRoAddOBXJs0a19y1bCtPPRe3oRST/6vVpgVz1cyl2O/o0L7oU7IKPSLBqe6lDodK1SU4OvsrMSmZ2wsyOm9lNddYxM9tvZk+b2eNm9qHYC5LmRZeCXfARCVbtAGcHDhQrXevuiSzACuBD5fuLgKeAD9Sscx1wiGhu2Q8DD8+23csuu8zbNjbmvmyZ+8hIdDs21v42OpXme4vI9DFYOfZqH+cEMO4N4mZiNXp3f8HdHy3ffxV4EriwZrVPAN8ql/Mh4HwzWxF7YdIcXEwDm4mkS40tepOjN7PVwFrg4ZqXLgSeq3o8wdk/BpjZkJmNm9n45ORk+wVIc1xqjYktkq4dO86uYA0MFKqfTOKB3szOA34AfM7df93JNtz9q+7e7+79y5cvb++P02wjr/b5IhKARAO9mc0lCvLfdvd76qzyPLCq6vHK8nPxSfO0TaeMIhKAxNrRm5kB3wRecffPNVhnA7CN6KLsFcB+d7+82XbVjl5E5GxptaO/EhgErjKzo+XlOjO70cxuLK/zI+Ak8DTw34F/n2B5RKRo1JcFSHCYYnf/CVGzyWbrOPDnSZVBRAqu0pelkkKtvm5WIPnvGSsixVBbe9+7Fx57bObgZZs3R48L1sxZgV5EIllPc9T2RJ8zB/7yL+H975/uy/Lmm9FgZgWjGaZEJJL1NEelVVv1QGV/9Vdwyy1w5gzMnw/velfapUyFavQiEqkOlFkdo722J/ratVGQf/31qHZ/8GAh+7Io0IvItKwP2VHbE310NErhVM8cV8C+LAr0IjItiSE7epX7r+2JPjwMP/oRfOELM3umQ6GGPwAFehGpSGrIjl4N113bE/3MmShHf+ZM9LjAPdM1w1Rc9u6NvrjVp7qlUvSlKljtQTKq2+9ws7+vBPuczOYUomLPMNUrmmREsq7bUR6bHQNZz/1nnAJ9XPLQYkGkG82OAQ3XnSoF+la1ckFJtRYpunrHgIbrTp0CfataSc20U2vppiVC1nswSvY1+g7ecMPZx4CG605fozkGQ106mjM2Ls3mf213Xspu5rEsyByYErB638HFi92XLNH3MiU0mTM29cDd7pJqoHePgjxEt9X27Dn7Cz02Fj3fSDcTh2vScUlb7XdwaKj9Y0Bio0AflySCa6MfjqT/ViQO+g4Go1mgV46+VTfcAJs2zbygtGlT9HynummJoFYMkjZ9BzNDgb4dZs0ft6OblghqxSDN9OJivb6DmaJA36q/+Zvpke8qbYQPHoye70Q3LRHUikGa6UXnPX0Hs6VRTqfbBfgG8BLwRIPXPwacAo6Wl52tbDfYi7GSjE4ucosu1hcQKeXo/xa4dpZ1HnT3S8vLrgTLEg/lJHtPQ0t0prrj0po1Z7+ufheFkligd/cfA68ktf2eU04yHRpaojPVlZLx8ajhgH4sCyvtHP1HzOyYmR0ys4sbrWRmQ2Y2bmbjk5OTvSzftE5zkurF2r04hpbo9nPI0udYWyk5eDBqOLBpk34si6pRTieOBVhN4xz9YuC88v3rgH9qZZup5+jbpV6s3Ysj39zt55Clz7HRdY2rr9b1pRwjrQ5TzQJ9nXWfBZbNtl7mAr17OhfG8nIRs9UA28r+dvs5ZPkCZ5bLLi0JMtADv8v0xCeXA/+v8rjZkslA79771jpZqoE20+oPVqv72+3nkMVWV3n5LkhTqQR64LvAC8BbwARwPXAjcGP59W3AceAY8BDwB61sN5OBPq3aVNFqcbPtb0g1+l6eceXl7E6aSq1Gn8SSuUCfdm0qizXQbjTa39By9END0WiPtaM/Dg11tj0pvGaBPu1WN/kXZw/Cdlt+FK3df7P97fZziLsn6JYtM1vCbNoUPd6ypbPtiTTT6Bcg1CVzNfo4tVOrTPtMoteyuL9jY+4LFkRnIAsWhF1WCR6q0edEO52HijYWSVb3N7peNX0rkoBKq5fM6O/v9/Hx8bSLka6dO6POQyMjUYcYyZ5SaTpds3077N8fBft771VHJumImR1x9/56r6lGnzVFy7vn1YEDUZA/eHBm79XPfz47PXAlMxToQ1d9AbbStX14GM47T+PtZNl73xsF9+pU08GDcOmlGsRNYqfUTeiqxy05fBjmzIHdu6fz0aVS9PyOHWmXVOJS+cy3bo3O2jQujbRAqZtQdDIwVvUF2NdemxnkK68ryJ8tS4OQ1YpjEDeRKgr0vdTp2Oo68NuX5XHsdR1G4tao3WWoS+bb0XfSjb5oQxnEJYv/tyz2B5AgoHb0AWm3dq4JTzqXxTOhrPYHkKAp0Pdau6flOvA7l8UUyI4dZ/8g6TqMdEmtbnqpunZeaTGj2X6Sof+1FIxa3YRCtfPe0f9a5B2q0YuI5IBq9NK5LLdHFxFAgV5mk+X26CICJBjozewbZvaSmT3R4HUzs/1m9rSZPW5mH0qqLNKFw4ejsXWqh0YeHlauWyRDkqzR/y1wbZPX1wPvKy9DwB0JlkU6tW5dNOzC+vVRe/T166PHqtGLZEZigd7dfwy80mSVTwDfKnfqegg438xWJFUe6dDAQFSDv+su+MM/jG6Hh9VEUSRD0szRXwg8V/V4ovzcWcxsyMzGzWx8cnKyJ4WTslIpqsF/5jPw4IPR7e7d2eh8JCJARi7GuvtX3b3f3fuXL1+ednHiFXqrlkqO/tChqIfpoUPFztGH/nmJ1JFmoH8eWFX1eGX5uWIJvVVLJUdfPdZOkXP0oX9eIvU0Gu0sjgVYDTzR4LUNwCHAgA8Dj7SyzcyPXllPr0ZZ3LPn7G2PjUXPN7J+vfvo6MznRkej57Osk/9F9XpZGxVTco8mo1cmGeS/C7wAvEWUf78euBG4sfy6AX8NPAP8DOhvZbu5DPTuUdCA6DYpnQyBm9dhc7vdr158XiJtSCXQJ7XkMtD3soao8fCndbpfef1/SKYp0IcsjRpzJ7XRvNZg292vvJ7hSOY1C/SZaHWTa70eZbGTMdqzOK57KzrZL42KKVnU6Bcg1CV3NfpeUo5+Wl73SwoL1egF6Kw2mtcabF73S6QOjUcvIpIDGo9eRKTAFOhFRHJOgV5EJOcU6EVEck6BXkQk5zLX6sbMJoFfdLmZZcAvYyhOSPK4T5DP/dI+ZUPW9ulfunvdcdwzF+jjYGbjjZohZVUe9wnyuV/ap2zI0z4pdSMiknMK9CIiOVfUQP/VtAuQgDzuE+Rzv7RP2ZCbfSpkjl5EpEiKWqMXESkMBXoRkZwrVKA3s2vN7P+a2dNm9vm0yxMHM/uGmb1kZk+kXZa4mNkqMyuZ2QkzO25mN6Vdpm6Z2Twze8TMjpX36Za0yxQXM+szs8fM7P60yxIXM3vWzH5mZkfNLPPD5RYmR29mfcBTwDVEk5UfBv6tu59ItWBdMrN/DbwGfMvdfz/t8sTBzFYAK9z9UTNbBBwBNmX5szIzAxa6+2tmNhf4CXCTuz+UctG6ZmY3A/3AYnffmHZ54mBmzwL97p6lDlMNFalGfznwtLufdPc3gQPAJ1IuU9fc/cfAK2mXI07u/oK7P1q+/yrwJHBhuqXqTnkSoNfKD+eWl8zXssxsJbAB+FraZZHGihToLwSeq3o8QcaDRxGY2WpgLfBwykXpWjnFcRR4CXjA3TO/T8CXgB3AVMrliJsD/9PMjpjZUNqF6VaRAr1kjJmdB/wA+Jy7/zrt8nTL3d9290uBlcDlZpbpVJuZbQRecvcjaZclAR919w8B64E/L6dIM6tIgf55YFXV45Xl5yRA5Tz2D4Bvu/s9aZcnTu7+K6AEXJtyUbp1JfDxcj77AHCVmd2VbpHi4e7Pl29fAg4SpX4zq0iB/jDwPjO7yMzeBWwB7ku5TFJH+cLl14En3X1f2uWJg5ktN7Pzy/fnEzUK+MdUC9Uldx9295XuvproeBpz98+kXKyumdnCciMAzGwh8G+ATLdqK0ygd/czwDbg74ku7t3t7sfTLVX3zOy7wP8B3m9mE2Z2fdplisGVwCBRDfFoebku7UJ1aQVQMrPHiSodD7h7bpoj5swFwE/M7BjwCPBDd/8fKZepK4VpXikiUlSFqdGLiBSVAr2ISM4p0IuI5JwCvYhIzinQi4jknAK9iEjOKdCLiOScAr1Im8xsnZk9Xh5jfmF5fPlMj1sj+aYOUyIdMLP/DMwD5gMT7r475SKJNKRAL9KB8nhJh4HTwB+4+9spF0mkIaVuRDqzFDgPWERUsxcJlmr0Ih0ws/uIhua9iGjaw20pF0mkoTlpF0Aka8zs3wFvuft3ynMR/28zu8rdx9Ium0g9qtGLiOSccvQiIjmnQC8iknMK9CIiOadALyKScwr0IiI5p0AvIpJzCvQiIjn3/wEVLDkzzY0VdwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# generate a noisy ellipse with random parameters\n",
    "minor = np.random.uniform(1,5)\n",
    "major = np.random.uniform(minor,8)\n",
    "center = np.random.uniform(-5,5,size=2)\n",
    "angle = np.random.uniform(0,np.pi)\n",
    "print(major,minor,center,angle)\n",
    "\n",
    "x,y = generate(minor,major,center,100,angle)\n",
    "\n",
    "# plot of generated ellipse\n",
    "fig,ax = plt.subplots()\n",
    "ax.plot(x,y,'xr',label='Data')\n",
    "ax.legend()\n",
    "ax.set_xlabel('x')\n",
    "ax.set_ylabel('y')\n",
    "ax.set_aspect('equal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fitzgibbon, Pilu, & Fisher (1996) proposed that the ellipse-specific fitting problem can be reformulated as \n",
    "\n",
    "$$\\min_{\\bf{a}} \\lVert{\\bf{D}a}\\rVert^2 \\text{ subject to } \\bf{a}^T\\bf{Ca} = 1$$\n",
    "\n",
    "where $\\bf{a} = \\begin{pmatrix}a,b,c,d,e,f\\end{pmatrix}^T$, the design matrix $\\bf{D}$ representing the least squares minimization \n",
    "\n",
    "$$\\bf{D} = \\begin{pmatrix}x_1^2 & x_1y_1 & y_1^2 & x_1 & y_1 & 1 \\\\ \\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots \\\\ x_i^2 & x_iy_i & y_i^2 & x_i & y_i & 1 \\\\ \\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots \\\\ x_N^2 & x_Ny_N & y_N^2 & x_N & y_N & 1\\end{pmatrix}$$\n",
    "and the constraint matrix $\\bf{C}$ \n",
    "\n",
    "$$\\bf{C} = \\begin{pmatrix}0 & 0 & 2 & 0 & 0 & 0 \\\\ 0 & -1 & 0 & 0 & 0 & 0 \\\\ 2 & 0 & 0 & 0 & 0 & 0 \\\\ 0 & 0 & 0 & 0 & 0 & 0 \\\\ 0 & 0 & 0 & 0 & 0 & 0 \\\\ 0 & 0 & 0 & 0 & 0 & 0\\end{pmatrix}$$\n",
    "\n",
    "The minimization problem can be reduced by using Lagrange multipliers, giving the two important relations,\n",
    "\n",
    "$$\\begin{align*}\n",
    "\\bf{Sa} &= \\lambda \\bf{Ca}\\\\\n",
    "\\bf{a}^T\\bf{Ca} &= 1\n",
    "\\end{align*}$$\n",
    "\n",
    "where $\\bf{S} = \\bf{D}^T\\bf{D}$ or $S_{ab} = \\sum_{i=1}^N x_i^ay_i^b$. There will be 6 solutions for $\\lambda$, and the minimal positive value corresponds to the ellipse solution. \n",
    "\n",
    "Halir & Flusser (1998) improves and numerically simplifies the method further. First, the design matrix $\\bf{D}$ is decomposed as $\\bf{D} = (\\bf{D}_1 | \\bf{D}_2)$ where $\\bf{D}_1$ and $\\bf{D}_2$ are the quadratic and linear parts respectively, expressed as \n",
    "$$\\bf{D}_1 = \\begin{pmatrix}x_1^2 & x_1y_1 & y_1^2 \\\\ \\vdots & \\vdots & \\vdots \\\\ x_i^2 & x_iy_i & y_i^2 \\\\\\vdots & \\vdots & \\vdots \\\\ x_N^2 & x_Ny_N & y_N^2\\end{pmatrix}, \\quad\\bf{D}_2 = \\begin{pmatrix}x_1 & y_1 & 1 \\\\ \\vdots & \\vdots & \\vdots \\\\ x_i & y_i & 1 \\\\\\vdots & \\vdots & \\vdots \\\\ x_N & y_N & 1\\end{pmatrix}$$\n",
    "\n",
    "The scatter matrix can also be decomposed as \n",
    "\n",
    "$$\\bf{S} = \\left(\\begin{array}{c|c}\\bf{S}_1 & \\bf{S}_2 \\\\\\hline \\bf{S}_2^T & \\bf{S}_3\\end{array}\\right) \\quad\\text{where}\\quad \\begin{matrix}\\bf{S}_1 = \\bf{D}_1^T\\bf{D}_1 \\\\ \\bf{S}_2 = \\bf{D}_1^T\\bf{D}_2 \\\\ \\bf{S}_3 = \\bf{D}_2^T\\bf{D}_2\\end{matrix}$$\n",
    "\n",
    "Similarly, the constraint matrix and vector of coefficients can be split as \n",
    "\n",
    "$$\\begin{align*}\n",
    "\\bf{C} &= \\left(\\begin{array}{c|c}\\bf{C}_1 & \\bf{0} \\\\\\hline \\bf{0} & \\bf{0}\\end{array}\\right) \\quad\\text{where}\\quad \\bf{C}_1 = \\begin{pmatrix}0 & 0 & 2 \\\\ 0 & -1 & 0 \\\\ 2 & 0 & 0\\end{pmatrix} \\\\\n",
    "\\bf{a} &= \\begin{pmatrix}\\bf{a}_1 \\\\\\hline \\bf{a}_2 \\end{pmatrix} \\quad\\text{where}\\quad \\bf{a}_1 = \\begin{pmatrix}a \\\\ b \\\\ c\\end{pmatrix}, \\bf{a}_2 = \\begin{pmatrix}d \\\\ e \\\\ f\\end{pmatrix} \n",
    "\\end{align*}$$\n",
    "\n",
    "Rewriting the first relation, \n",
    "\n",
    "$$\\begin{align*}\n",
    "\\bf{S}_1\\bf{a}_1 + \\bf{S}_2\\bf{a}_2 &= \\lambda\\bf{C}_1\\bf{a}_1 \\\\\n",
    "\\bf{S}_2^T\\bf{a}_1 + \\bf{S}_3\\bf{S}_2 &= \\bf{0}\n",
    "\\end{align*}$$\n",
    "\n",
    "from which we can obtain that $\\bf{a}_2 = -\\bf{S}_3^{-1}\\bf{S}_2^T\\bf{a}_1$. The second relation on the other hand would simply be $\\bf{a}_1^T\\bf{C}_1\\bf{a}_1 = 1$, which gives the final set of equations\n",
    "\n",
    "$$\\begin{align*}\n",
    "\\bf{Ma}_1 &= \\lambda\\bf{a}_1 \\\\\n",
    "\\bf{a}_1^T\\bf{C}_1\\bf{a}_1 &= 1 \\\\\n",
    "\\bf{a}_2 &= -\\bf{S}_3^{-1}\\bf{S}_2^T\\bf{a}_1 \\\\\n",
    "\\bf{a} &= \\begin{pmatrix}\\bf{a}_1 \\\\\\hline \\bf{a}_2\\end{pmatrix}\n",
    "\\end{align*}$$\n",
    "\n",
    "where $\\bf{M} = \\bf{C}_1^{-1} (\\bf{S}_1 - \\bf{S}_2\\bf{S}_3^{-1}\\bf{S}_2^T)$ is the reduced scatter matrix. The minimal non-negative eigenvalue $\\lambda$ would give the appropriate eigenvector $\\bf{a}_1$ of matrix $\\bf{M}$, for which $\\bf{a}_2$ and consequently the best fit parameters $\\bf{a}$ can be obtained."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 390,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fitting an ellipse \n",
    "# based on Fitzgibbon, Pilu, & Fisher (1996) and Halir & Flusser (1998)\n",
    "def fit(x,y):\n",
    "    D1 = np.vstack([x**2,x*y,y**2]).T\n",
    "    D2 = np.vstack([x,y,np.ones_like(x)]).T\n",
    "    S1,S2,S3 = D1.T @ D1, D1.T @ D2, D2.T @ D2\n",
    "    C1 = np.array([[0,0,2],[0,-1,0],[2,0,0]])\n",
    "    M = np.linalg.inv(C1) @ (S1 - S2 @ np.linalg.inv(S3) @ S2.T)\n",
    "    vec = np.linalg.eig(M)[1]\n",
    "    cond = 4*(vec[0]*vec[2]) - vec[1]**2\n",
    "    a1 = vec[:,np.nonzero(cond > 0)[0]]\n",
    "    return np.vstack([a1,np.linalg.inv(-S3) @ S2.T @ a1]).flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 391,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.38132347,  0.35651649, -0.8529293 ,  1.1310537 ,  3.72028399,\n",
       "       -4.50970767])"
      ]
     },
     "execution_count": 391,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "params = fit(x,y)\n",
    "params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimation of Errors\n",
    "To estimate the error for the fit parameters $\\bf{a}$, we start by defining the auxilary dataset $z_{ij}$\n",
    "$$z_{ij} = \\begin{pmatrix}x_j^2 & x_jy_j & y_j^2 & x_j & y_j & \\alpha_j\\end{pmatrix}^T$$\n",
    "where $\\alpha_j = 1$. We then obtain the goodness of fit of the data using the chi-squared and equating it to the degrees of freedom i.e. $N-6$, giving \n",
    "$$\\chi^2 = \\sum_{j=1}^N \\frac{\\left(\\sum_{i=1}^6 a_iz_{ij} - 1\\right)^2}{\\sigma_j^2} = N - 6$$\n",
    "where $\\sigma_j$ is the uncertainty for data point $(x_j,y_j)$. The best fit parameters $\\bf{a}$ will minimize the distribution. After derivation based on McDonald, 2014, the uncertainties for the parameters can be obtained using\n",
    "$$\\sigma_i = \\sqrt{\\frac{\\sum_{j=1}^N \\left(\\sum_{k=1}^6 a_kz_{kj} - 1\\right)^2}{(N-6)\\sum_{j=1}^N z_{ij}^2\n",
    "}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 392,
   "metadata": {},
   "outputs": [],
   "source": [
    "# estimating the errors using chi-squared\n",
    "def errors(x,y,coeffs):\n",
    "    z = np.vstack((x**2,x*y,y**2,x,y,np.ones_like(x)))\n",
    "    numerator = np.sum(((coeffs @ z)-1)**2)\n",
    "    denominator = (len(x)-6)*np.sum(z**2,axis=1)\n",
    "    unc = np.sqrt(numerator/denominator)\n",
    "    return tuple(ufloat(i,j) for i,j in zip(coeffs,unc))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 393,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.3813234733723383+/-0.007463327642510613,\n",
       " 0.35651649404035146+/-0.010273363681267213,\n",
       " -0.852929304299267+/-0.009995413138221472,\n",
       " 1.1310537012738004+/-0.03402667857190656,\n",
       " 3.7202839927133007+/-0.03735161370496613,\n",
       " -4.5097076661287+/-0.11284004340463694)"
      ]
     },
     "execution_count": 393,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "params = errors(x,y,params)\n",
    "params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Converting Fit Coefficients to Ellipse Parameters\n",
    "The fit parameters obtained from the fit (as well as their estimated errors) correspond to the polynomial parameters. Therefore, they have to be converted into the usual ellipse parameters i.e. center coordinates, semi major and semi minor axes, and the rotation angle, in order for the result to have a direct physical significance. This gives the equation of the ellipse as\n",
    "$$\\frac{1}{a^2} \\left[(x-x_0) \\cos{\\phi} + (y-y_0)\\sin{\\phi}\\right]^2 + \\frac{1}{b^2} \\left[(x-x_0) \\sin{\\phi} - (y-y_0)\\cos{\\phi}\\right]^2 = 1$$\n",
    "or parametrically,\n",
    "$$\\left.\\begin{matrix}x = x_0 + a\\cos{t}\\cos{\\phi} - b\\sin{t}\\sin{\\phi} \\\\ y = y_0 + a\\cos{t}\\sin{\\phi} + b\\sin{t}\\cos{\\phi}\\end{matrix}\\right\\} 0 \\leq t \\leq 2\\pi$$\n",
    "where $(x_0,y_0)$ is the center of the ellipse, $a$ is the semi major axis, $b$ is the semi minor axis, and $\\phi$ is the rotation angle of the ellipse's major axis with respect to the horizontal (counterclockwise). The errors for these parameters are propagated using the `uncertainties` package. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 394,
   "metadata": {},
   "outputs": [],
   "source": [
    "# converting the coefficients to ellipse parameters\n",
    "def convert(coeffs):\n",
    "    a,b,c,d,e,f = coeffs\n",
    "    b /= 2\n",
    "    d /= 2\n",
    "    e /= 2\n",
    "    x0 = (c*d - b*e) / (b**2 - a*c)\n",
    "    y0 = (a*e - b*d) / (b**2 - a*c)\n",
    "    center = (x0, y0)\n",
    "    numerator = 2 * (a*e**2 + c*d**2 + f*b**2 - 2*b*d*e - a*c*f)\n",
    "    denominator1 = (b*b-a*c)*((c-a)*umath.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))\n",
    "    denominator2 = (b*b-a*c)*((a-c)*umath.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))\n",
    "    major = umath.sqrt(numerator/denominator1) if numerator/denominator1 > 0 else 0\n",
    "    minor = umath.sqrt(numerator/denominator2) if numerator/denominator2 > 0 else 0\n",
    "    phi = .5*umath.atan((2*b)/(a-c))\n",
    "    major, minor, phi = (major, minor, phi) if major > minor else (minor, major, np.pi/2+phi)\n",
    "    return center, major, minor, phi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 395,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((2.773541669303805+/-0.09466162130269165,\n",
       "  2.760543764290714+/-0.05361930769043788),\n",
       " 2.61210485340965+/-0.14555611751219827,\n",
       " 1.550347820215571+/-0.07933262615860551,\n",
       " 0.32365321552125015+/-0.009408231946422919)"
      ]
     },
     "execution_count": 395,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "center, a, b, phi = convert(params)\n",
    "center, a, b, phi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 396,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generating fit ellipse line\n",
    "def line(coeffs,n=100):\n",
    "    t = np.linspace(0,2*np.pi,n)\n",
    "    center,major,minor,phi = convert(coeffs)\n",
    "    x = major*np.cos(t)*np.cos(phi) - minor*np.sin(t)*np.sin(phi) + center[0]\n",
    "    y = major*np.cos(t)*np.sin(phi) + minor*np.sin(t)*np.cos(phi) + center[1]\n",
    "    return x,y\n",
    "\n",
    "# alternative using matplotlib artists\n",
    "def artist(coeffs,*args,**kwargs):\n",
    "    center,major,minor,phi = convert(coeffs)\n",
    "    return Ellipse(xy=(center[0],center[1]),width=2*major,height=2*minor,\n",
    "                     angle=np.rad2deg(phi),*args,**kwargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 397,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEGCAYAAABrQF4qAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA2nUlEQVR4nO3dd5hTVfrA8e+R3qWLICKKiA2RYldGxQIsAquIuqPs8ltExbqKsAoqKAgK9gprA0URxYLrAjqDBRRm6FVEQB1AGYYifWaY9/fHyUgISSbt5qa8n+e5zyQ3NzfvneS+OTn3FCMiKKWUSl1HuB2AUkopZ2miV0qpFKeJXimlUpwmeqWUSnGa6JVSKsWVdzuAcNWrV0+aNWvmdhhKKZVQ5s+fv0VE6vt7LOkSfbNmzcjNzXU7DKWUSijGmJ8DPaZVN0opleI00SulVIrTRK+UUiku6ero/SkqKiIvL499+/a5HYrjKleuTJMmTahQoYLboSilkkRKJPq8vDxq1KhBs2bNMMa4HY5jRISCggLy8vI47rjj3A5HKZUkUqLqZt++fdStWzelkzyAMYa6deumxS8X5bLRoyE7+9B12dl2vUo6KZHogZRP8qXS5TiVy9q3h169Dib77Gx7v317d+NSEUmJqhulVIxlZMDkyTa533ILvPSSvZ+R4XZkKgIpU6J3W7ly5TjjjDM45ZRTaN26NWPGjKGkpCToc9avX88777wTpwiVClNGhk3yw4fbv5rkk1b6JXqH6h6rVKnCokWLWL58OTNnzuTzzz/nkUceCfocTfQqoWVn25L8kCH2r+95Ewt6LSA+RCSplrZt24qvFStWHLYuoKwskXr17F9/9yNUrVq1Q+7/9NNPUqdOHSkpKZF169bJ+eefL23atJE2bdrI7NmzRUTkrLPOkpo1a0rr1q1l7NixAbfzFdbxqtQzatThn9esLLs+Vhw6T1x7nTQA5EqAvOl4YgbKAQuBaX4e6wPkA4s8y/+Vtb+oE73IwQ/TkCEx+1D5JnoRkVq1aslvv/0mu3fvlr1794qIyOrVq6X0GLKzs6VLly5/bh9oO1+a6NNcPJJjPL5MvPcbyvkYz5iSULBEH4+LsXcCK4GaAR5/T0QGxCGOg7zrHocMcbzusaioiAEDBrBo0SLKlSvH6tWro9pOpbl4XCgdOND/6zpxroR6Ppa2BCo91tKWQJMnxz6mFONoHb0xpgnQBRjv5OuELQ51j2vXrqVcuXI0aNCAp556ioYNG7J48WJyc3MpLCz0+5xQt1N+pFtdbypdKA31fPT+ghs69NCkr4Jy+mLs08BAIFjzk78aY5YYY6YYY47xt4Expp8xJtcYk5ufnx9dRN6lgGHDDn5wYpjs8/Pz6d+/PwMGDMAYw44dO2jUqBFHHHEEEyZM4MCBAwDUqFGDnTt3/vm8QNupEKRbu+94XCiNh3DPx1T6gounQHU60S5AV+BFz+2O+K+jrwtU8ty+Gcgqa79R19E7VM93xBFHSOvWreXkk0+W008/XZ544gk5cOCAiNj69tNOO01OP/10GThw4J/1+YWFhZKRkSGnn366jB07NuB2vrSOPgAHrr2EJN51x07V0btRBx7ua7r1HicB3LgYC4wE8oD1wG/AHmBikO3LATvK2m9MLsYmubQ53kgSz5Ah9mM9ZIizsfnG5NDF0X37RNasEfnuO5FPPxVZu1akcMQTMrD3evnHP0R69RJ5/HH7mrefmyutW4uceKJI8+b2+a+8IlK7tkidOjak7GyRn38WOfpokcaNRY4/XmTECLvtbVf9Kr0rfiD/7LJBhg+3+8w98mKZOmyJzJ0r8uuvIsXFER5ILL5EtIVOUK4k+kNeJHCJvpHX7R7A92XtSxN9Gh1vuCe2m6W9CF57/36bxOfPt/enThXp00ckI0PkuONENm4UmT7d3m7fXuTKK0VmzhQpKbHJfdw4kXffFfn+e/v8FStEFiwQWbXKJnMRkb17RQoKRPLzRTZvtl8cRUUieXkiv/wi8uOPIhs22G2//FLk7QeWyUvV/iXju04VqVdPJj24TLp1E2nXzn4x/PSTyLffipxxhkjnziJ33CGSk2Ofv2GDjS3o/yeaJK2tboJKqEQPDAO6ycFS/3JgMZANnFTWvjTRp9nxhppAE6G0F+DXxNatNjmOGycybZpd17GjSIUKIscea0vlIiJffGG3mTFD5IcfbEJ2RRm/inbtEsnNFfn4Y5Enu2TJvBdzZO9ekaOOEqlWTaTNCX/IE53t/33pUpFt2zxP1GoXR7me6GO5aKJPv+MNqTrG7dJeVpYU1j1KFv7zBRlf/U659apf5euvRbZvF6lZU6RDB5GbbhL54AO7+W+/+SRyt+P3fs1wkrHPF+q2T76WubU6Se7L80RE5IYbRKpXF2nRQuS220RkyBDZREPZM+gRhw8k/WiiTzFpdbwJWgosLral1TfeEBnQ/VfZVOdk+fLJBdKqlUhmp00yttoDsmbCHBEJUp3hLRF+kUQaQxnvUen/6n+jFonUqyeDz/tKqpudcsmZBTJypK1OUtHTRJ9i0uZ4EyH5eezcaevLR4ywiXvUKFtKve46kTFdvpT8qd8cHntpaTzU0rrbX2rR/Koo61eXz3u349Ov5JMa18ug69ZLSYnIq6+KXHON/ZuXF+VxiCTOL6Q40kSfYtLmeF08WbduPXiR8b77bN3zBReIPPCAvYgaUim9VDhfWG60GopWKF9QZbyXGzeKvP66yPXX2xZCGzfaJeKkn0CFhHjRRB8Hpe3oS5d169bJOeecIyIi69atk7fffjtmr5UIx5uqRowQad3a1iv37m3XbdhgW69EJZRk6EaJPtovUwcSamGh/Tt5sm0aev75Is88I7JjR5g7cvsXUpxpoo+DQJ2bRA4fvCxaiXC8qWDdOpEXXxTp1k3k2mvtug8/FJkz52CyialgpXW3SqDRvq7Dv7r27bMtlW66yf7KWrRI5L//DaM9fzL+QoqQJvo48JfoS9f5DkccrUQ43mRUWGjboc+YYe+fe65IZqbIO++IbNni8IuXVbp0s045iUq+WVkibduKNGsmMnKkbeoZdOMkOa5Y0EQfB95VN927dxeRg4k+ZUv0SXLBq7hY5O9/F6lb1zZzfOutOAeQDPXFSVbyzckRGTDAlvjnzBFZuNBng2T4n8dYsESfkjNMPfwwGHNwmT/fLt7rHn7Ybnv00QfXtW1r1/Xrd+i2GzeW/ZqlM0wtWrSIqVOnOnVoiSXcgcTiNMLknj3w4Ydwww12nKxy5eDSS2HhQpg7FzIzQ9hJLGPNyTl0lMXSURhzcsLflxOScIC0du3gueegUiXYsAE6d4arroLcXM8GTzwBgwcf+j8fPNiuT0eBvgESdUnUEn2wqpuULdGLhPfz2MFS1v79B7v4N2ggcskltv5906YynhjoV0m/fulRIkzmkq/Xe7dnj8hzz4lc3n6LlDw+SrZ+/HXyHleE0Kob5wVL9Lm5uXLhhRfG7LUS4XgPEc7P/hjWm5aUiMyda3tc1qsn8uijdv3OnWHsJFiiS4c63iSpfvPL+70aNUpkzJg/72dkiFzecp0sr9I2td8/L5ro4yBYovcdjjhaiXC8f4okGUZZH/zzz7an5f79ImefLTJ8uG1BE7Fgx5Bkdddpp/S9y8wUMcYmexHZPz1bnq72b6lXaYc8x21p8f5pok8xCXO8kfzsj7CUXFwsMmGCHdmxbl2RF16IQfze/CX0dCjRp4LS9y4z89D3a8wY+b3OSbLylmdkZ91jZcLg5eKZIiIlaaJPMQlzvJFOGhHGF8OaNXbkx5ISkX/8Q2TKFAfGRvGX0JO57jqd+L53mZmHJn3P+7XunTlyVvkcOavVdlm2zOWYHaKJPsUk7fGG8cXw8ccil19uz9Wnn3YwpkAJvV+/5K27The+792YMbb6JjNTpGrVP6txSh34Ikte7fm5NGxoL9ynmmCJ3tjHk0e7du0k9882VNbKlSs56aSTMMa4FFX8iAirVq2iVatWbocSc7//DitXQseOMGAAdOgA11wDVao4+KKjR9vmoN5zj2Zn26aPAwc6+MIqat7vXWnT3sGDobj4YNNfP5OH79gBtWrBG29Ap07QuLE74ceaMWa+iLTz91j5eAfjhMqVK1NQUEDdunVTOtmLCAUFBVSuXNntUGJq5UoYOxamTLHzPXfsCM8/H6cX95fMMzJ00ulk4P3e+fZVgIN9FXzey1q17N8NG+DMM+GFF+Dqq+MQr4tSokRfVFREXl4e+/btcymq+KlcuTJNmjShQoUKbocSFRFYtw6aN4fbboMGDeDWW6F+fbcjU+lk7ly4/noYMQKuvdbtaKITrETveKI3xpQDcoENItLV57FKwFtAW6AAuFZE1gfbn79Er5LHgQMwdartoLhvHyxYYHuuxoRWw6gIFBRAtWqQlwcNG0KNGm5HFJlgiT4eQyDcCawM8FhfYJuInAA8BYyKQzzKBSUl9u/DD9tqmkGDYpzkIfwhGZQC6taFypXhzTfh/PPh55/djij2HE30xpgmQBdgfIBNrgLe9NyeAlxiUrmSPQ0VFcG4cXDiibbENGQIzJkDPXrEOMnDwTFkevWCoUMDXoxTyp9hw6BPHzjnnMQZhihWnC7RPw0MBEoCPN4Y+BVARIqBHUBd342MMf2MMbnGmNz8/HyHQlWxlptrE/z778OECdCkCVSs6PCLZmTYK7rDh9u/muRViIyBu4tG8/KtS9i/3+sBBwbeizfHEr0xpiuwWUTmR7svEXlVRNqJSLv6erUuMQQY3bF45BO8/rodKbJ5c5g4EWbMsKWkuEjCkRhVAmnfnm7PXML5Rdk8/DDMezE3Nar/AjWwj3YBRgJ5wHrgN2APMNFnm+nAOZ7b5YEteC4QB1r8dZhSLvDprFLyZZZMqvF/cvzRuyUjw8/44C7EpL1ZVUQ8n5tPrp0oDczvMu/FHLcjCglu94wFOgLT/Ky/DXjZc7s3MLmsfWmiTyCeE2LX/cNkd91j5C/n5LubU5N5JEaVWDzj53xy7URp0EAkGTqjB0v0cWlHb4zpCNwrIl2NMcM8AX1ijKkMTADaAFuB3iKyNti+tHll4li9Gu7vtgLzwyo+HLLIXs1SKtmVtta65RZ46SVyh39O677tSPSuK243r0REZomnDb2IDBWRTzy394nINSJygoh0KCvJq8Tx+ONwXodCzvn1fd4ZvEzrw1VqKE3ykyfbgsvkybQbciW/fTiHnj1t349klJJTCSpnFBbCK6/YoUQuq53DynKnMXDahVQeMfRgs0ZN9iqZBZj2scn6bylfHvr3t726k40mehWS7Gxo3Ro++QT++APO3JFNvSkvJ+48qOGK03y2KsENHHh4k9yMDMz9A3n9dVi8GJ56yp3QopESg5opZy1aZDuSPPOMnYDZGFJvMDDf0Q69f8IrhR0m4eOPYetWtyMJn5bolV/FxfDcc7bq/Ywz4IcfoHt3T5JPRdqrVoWgaVM45RR49PJv2D991qEPJvAvQE306jBz59qx4D/8EC66yK5LsZGR/dNetSoE5ctDzq6TeLTH/KQZV0kTfSqKsL659CLThAnwr39BVhacfLJDMSYi7VWrQmAMvDylPq9Uup0FPR8N/Regm9eBAjWwT9RFO0yFIIIeot9+K9Kmjcjvv8cpxkSjvWpVmCZOFHnuik/lsEnlA3H4M4bbPWNjuWiiD5G/Ca/92L1b5O67RRo1EvnwwzjHmEi0V60Kl+cc2z1oWOgJO8TzMhKa6NOVpxt3oNJGcbHIhg0iffuKbNkS59iUckK8vrA9CXvlG9/LCSeI7J+eHXriLuO8jFSwRK919KkqSH3z/v22Dr5vXzj6aBg/3k6+oFTSi9fkM56OVSfddBbNm8N/fuoYWj+S0vPykkvg2WcPrbN3sr4+0DdAoi5pVaKPtHQSpC5w1SqRM84Q6dlTJD/fmbCVcpWD1SP+5OSIHH20rQYNaNQokTFjDsaTlSVSrZpI5coH70cZK1qiT1KRlk4CdOMmJ4c5c+Dmm2HKFKhXz9nwlXJFnJvJtmtnp8cMqn17++t68OCD8ZQrZ5fHHnO+30agb4BEXdKqRC8Sk9LJjh0i11+f5hdbVfqI5pyJ8Fd0SYnI/Pn2b8D9lcaVmSlStaq9H8P6evRibJKL4sOQkyNy/PEi/fqV8dMyFNoyRSW6aJswRvj8AwdEjjvOnm9B95eZac/lzMyYVzFpok9mUX4YevYUmTw5xrFoW3OVqGJRGInwnHvsMZF//tNPHKX769TpYJKvWVOkVq2Ynkua6JNVhIl1zx6Re+4R+e03B2OK04UupVwRwa/ojRtFjjxS5I8/5PBztTTJd+pk7/frZ5O99/kT5a9jTfTJKoLSyS+/iLRtK9K7dwyqagJxqB2wUgkhisLMFVeIvPuuz35Kq2s6dTq84BbDak9N9Gli3z6R5s1FnnjCz0WhQML9MtESvUplUVZPPv+8yE03ea3wrpOPYH/hcCXRA5WBecBiYDnwiJ9t+gD5wCLP8n9l7VcT/eFKSkSys+3tsMeqCeeDrXX0KtVFWcf/4492OJGSEjnYVj4z09GSfCm3Er0BqntuVwDmAmf7bNMHeD6c/Wqi9zJqlBTOyJa+fUVOO82rbjDcD1GopXRtdaNUmW6/XWTHp1/FvVAULNE71mHK89q7PHcreJYknG0xce069Wy6dhF+X76FOXOgRm6E3b1D7WASYJo1v7NNKZWmnn0Waq74PmCnRTc42jPWGFPOGLMI2AzMFJG5fjb7qzFmiTFmijHmmAD76WeMyTXG5Obn5zsZctIQgUqdLuQv/Rsz9cdTqT46ilmRdBx2pSLjZ4z52c8t4IZJXROrUBSoqB/LBTgSyAZO9VlfF6jkuX0zkFXWvrTqRmTFCpFzz/VU1YhE1wpG692Vipyf8+eH2h2keaM9cQ8Ft8e6EZHtnkR/hc/6AhHZ77k7Hmgbj3iS2TffQMeO0L8/1KhB9KXxIOPiKKXK4Geu4eaTR7FpexX27HE7OC+BvgGiXYD6wJGe21WAb4CuPts08rrdA/i+rP2mc4l+507bfHLmTM8KLY0rlRh8flX/9a8ia9fGNwRcKtE3ArKNMUuAHGwd/TRjzDBjTDfPNncYY5YbYxYDd2Bb4Sg/cnOhWjVYvhwuvdSzUkvjSrnDu26+9Fd1ZiaMGQPZ2UyZAscd526I3oz9Ikge7dq1k9zcXLfDiKvXX4cHH4T58+Goo9yORin155DhgwfDyJGH/Z1+3xdsbtSazMz4hWSMmS8i7fw9Vj5+YaiARo+2TSK9r9JnZ0NODq/UGsijj0JWliZ5pRJG6a/nv/wFeva0Sb7013WbNqx5ZjPLjyauiT4YnXgkEQSYYETatWfFCnu3ZUt3Q1RK+cjIgHvugQkTDu1/kpHBkdd0YscOd8Pzpok+Efi5cv/SDd+yunEGzzwDJ5zgdoBKqcMEafFWqxZs3+5eaL400ScKr96pL7Z/jTHTWlK1qttBKaX8Kq2jnzwZhg07WFDzJPvLLoP333c5Ri+a6BOFp3Twbo/3GDG9LTOHfccxfvsJK6VcV0aLt/377XW1RKGtbhKBV+lgTqUMaq6ax6n3d3F2smCllGOWLoXrroNly+L3msFa3WiJPhHk5DBn6P8Ymp3BuefCqf/ooO3hlUpiu3fbfi+JQhN9vPgZ/IjsbBg9miVXDKTHo2057zyvx3RUSKWS1p49HH6NLUgOcJom+ngJ0IRyQ7Pz6NzZDm16+eUBnuviB0QpVQY/52errbMZ0uLdQ7cLkAPCHlY8EoHGRkjUJanHuvEzwUdhocgXX4T4PB3TRqnE4+f83F+3kez936zA2zowFSc6Z2wC8Qx+dODBoXLrrSKrV4f4PJ2rVanE5XN+vnLPSunbN8C20QwrHkSwRK9VN7EWrJrFq4PFI2Oqs+irHTRtGuJ+Q50FSikVfz7n59oKJ/kf1MytSX4CfQMk6pLwJfpA1Sxjxvy5/r33RJo22Cu/1WkVeslcS/RKJS6f8/PqC3+XSZMCbONQFSxaoo8jP8MZMHkyFBf/2S7+hBPg4+mVaTjlhdCaUJbRC08p5SI/5+eZua9yxn6fmVNdHFZcO0w5ZehQ+zNuyBD75mPHvnjySXv3iHC+YoOMbqlNMJVyWYKcn8E6TGmid0LpN/wtt9h6uMmTkY4Z9Oplhxp+7jmf7RPkg6JUSnD5fNqyBXr3hi++cPylDqE9Y+MpQDXLuHtX8eOP8MQTfp7jZvtapVJNqOeTQ/1T5s+3NbWJRBN9rAWoh/tt6RbefRcqV/bznED1+tqyRqnwhXo+OVTAmj2bQ3u5J4JAV2mjXYDKwDxgMbAceMTPNpWA94A1wFygWVn7TfhWNz527xaZOzfEjR1qX6tUWgrlfHKgNdvtt4v8739R7yZsuNFhCjBAdc/tCp5EfrbPNrcCL3tu9wbeK2u/yZbob7lF5IYbQthQm08qFTvhnE8pUsByJdEf8iJQFVgAnOWzfjpwjud2eWALngvEgZZkSvTZ2SJNmohs317GhjrEgVKxE875FOMC1rJlIiNGRLWLiAVL9I7W0RtjyhljFgGbgZki4tOwlMbArwAiUgzsAOr62U8/Y0yuMSY3Pz/fyZBj6vHH4YUX7LRiQbnYvlaplBPq+eRA/5QZM2D9+shDd0pcmlcaY44EpgK3i8gyr/XLgCtEJM9z/ydsqX9LoH0lRfNKQAQKC6FSJbcjUUr55UAzzAsvtE/t2jVGMYbB9eaVIrIdyAau8HloA3AMgDGmPFALKIhHTE5avBguvRQqVnQ7EqVUQAMHHt4SJ4p5IHbuhLVr7bmfaBxL9MaY+p6SPMaYKkAnYJXPZp8AN3luXw1kSTx+YjjowAH45z/tNGLGuB2NUipeatSw1TZ+m1C7zMkSfSMg2xizBMjB1tFPM8YMM8Z082zzH6CuMWYNcA8wyMF44uLll+3MMn37uh2JUiqeRo6E3393Owr/dAiEGPviC2jUCE45xe1IlFLxsm0bNGsGGze6N1es63X06WL2bLjgAk3ySqWbd9+Fyy6LMMnHYapQTfQxkpcH3bpBErX+VErFyEcfQf/+ET45DmNdlY/ZntLcQw/BzTdDkyZuR6KUirdp06B8pNnUe2werxFvYznWlZboY+Dnn+0bff/9bkeilIq34cPhxx+jbGXn8FShmuhj4NhjYenSEHrAKqVSyubNMHasbYARFYfnktVEH6U1a+wb3aCB25EopeJt/Hjo2RNq145iJ3GYKrTMRG+Mud0YE81hJLcyroiPGWObViml0s/XX8Ndd0W5kziMdVVmO3pjzKPYIYQXAK8B093svRr3dvTe37YZGYfc33xKBi1bwqpV0LBh/EJSSiUGkcTpAR9VO3oReRBoge3F2gf40RgzwhhzfEyjTFRBZquZNQuuv16TvFLpZtcuOOss2L3b7UhCE1IdvacE/5tnKQZqA1OMMbFr0Z/IAlwR79ULnn/e5diUUnH3wgu2J2z16m5HEpoyW34aY+4EbsROCjIeuE9EiowxRwA/As5Pq+423yviGRlMKchg0ya4/Xa3g1NKxdPOnfba3KxZbkcSulCa+NcBeorIz94rRaTEGOPCqMtx5ltHn5EBvXrx6rE/0Pe+Om5Hp5SKs4ICuOceOPlktyMJnQ5qVhY/kxNsmDyb0/qcyYaCKlSpEr9QlFIuGj2azcs2U7lXN2p2vdCuy862A90cf3zE49jHSrCLsToEQln8vHnLap3HP27l8CTvwIw1SqkE0b49Ax7YzpmTsxj0+QG7rnt32+xm6lRXQyuLJvoIXH65XQ5TOjiRn6aYSqnk9mVJBrn19/Lmrv7Q9QnbtrJiRZvkYzxkQaxpz9gwLV4Mt94a4MEgTTGVUslLBO6+G556qQpV7roZ9uyBvXvhjjuS4vzWRB+mzz4rYy5YhwcnUkrFnzHw4YfQrUY2PPOMnUauShV49tmYj0vjBCfnjD3GGJNtjFlhjFnuaabpu01HY8wOY8wizzLUqXhiZfr0ANU2pRwenEgpFV9r18LgwXDCr9mYHt1t1p82zZb6RKBHj4Q/z52soy8G/iUiC4wxNYD5xpiZIrLCZ7tvRCQpmmkWF8Mff8BFFwXYIEBTTK2+USo5HTgAffrAVVdhG1X07m2X0vP5o49sq5ucnIQ+xx1L9CKyCdjkub3TGLMSaAz4JvqkUX7saBaObQ9VA7SqCTY4UQJ/CJRS/j39tP17111AOT8t50oLdAkuLnX0xphmQBtgrp+HzzHGLDbGfG6M8TvbqjGmnzEm1xiTm+/iXH1P/9yD73qMDjzl18CBh7/pGRnatFKpJLVrF7zxBpQr53Yk0XG8w5QxpjrwFfCYiHzo81hNoEREdhljOgPPiEiLYPuLe4cpLyefDBPvzOHMBzs7NuWXUsp9RUWwZAm0bet2JKGLavTKKF+4AvAB8LZvkgcQkT9EZJfn9n+BCsaYek7GVKYA48/vfWws69bBqX9vr61qlEpxw4fDww+7HUXsONnqxmCHNl4pImMDbHOUZzuMMR088RQ4FVNI/M3I3rUrq7fVo1UrqDjb06omM9OObJTgV9uVUuGZPh3+8x949VW3I4kdJ0v05wGZwMVezSc7G2P6G2P6e7a5GlhmjFkMPAv0dnNSE8B/p6fhw2n95r+Ye93T9v7gwfD55/ZrP8ZTfiml3CNiT+tJk2IwD2wCcbLVzbdA0LlXROR5IPFGdPfu9DRkCNxzD9P2XUqLh26i5dVXwsiRB+vm27TRVjVKpYD9+6GkxA4/XD7FBofRnrH++On09MT008n7650wYcKhdfPaqkaplHDHHXZu7lRL8qCDmh0uQKendcU/ccLSZw6ZfERL8UolqDBHkn39dTvR97x5cYwxjrRE78tPpycZNJjN2yvS4J2n7Vd+aR2+1s0rlZj8Narw7vPi5bffYNAgO5ZNjRpxjjNOtETvy8+3vRQVk/Xscqpc4Rn7QHu8KpXYvBtVBOnzsn8/HHWUbTPfsKFLscaBluhDsO+OgTS40qfnRFl18wHa4zM6PeZTV8p1ZYwku2ULnHEGrFyZ2kkeNNGHZP58O7BRWML46aiUckCQkWT37oVu3ezAk61auRhjnGiihzJL3wUFULdumPvUSUiUco93owo/19UeegiOOw4efdTlOONEEz2EVPo2QXsEBKCTkCjljgAjycq8HPbsgQcegNdegyPSJAM6PqhZrDk2qFlpcvdz4WblSlt987e/xW6fSqn4e+ghO5HIhAluRxJ7rg1qljRKL5B6l7691rdqFUWSD/DTUSkVX489Bu+/b4eoSjea6MFW0fToYeeCHDLE/u3R48+qm1mz7BhmYQk2CYlSKq6WLYOJE+HLL6FBA7ejiT9tR19K5GBFvDH2vkdhoe1UERZ/TS+1N61ScbdqFZx6KixaBJUquR2NO7RED7aU/dFHdrCL4cPt348++rP0Xbu2bXmjlEouL78MnTvDnj0+ST7N+rloooeDpW/vNrde6489Fuq5Ox2KUipMzz0HI0bAzJlQtarPg2nWz0UTPZR54bRBA5gxw+UYlVIhKyyEb7+1A5Ud/4Gf0jtAz55p089FEz2EdOF06FDYsMGl+JRSISkpgccft5N6v/ceNGtG4NJ7795p089FEz3YKhrfN9lnLJtZs2D16viGpZQKXVER3HQTfPaZT0eoQL3UIeAQCakm9RN9jC66tGhhr9678dpKqeBE4NprYetWO+frkUf6bODbSx3Sqp+Lk5ODH2OMyTbGrDDGLDfG3OlnG2OMedYYs8YYs8QYc2bMA4nRRZezz45gUoI0u+CjlBv27bMtoksbyx124RUOH+Ds3XfTq5+LiDiyAI2AMz23awCrgZN9tukMfI6dW/ZsYG5Z+23btq2ELStLpF49kSFD7N+srLB3sWePSGFh+C8di9dWSvm3YoXI8ceLLFwYZKPSc7D03PO9nyKAXAmQNx0r0YvIJhFZ4Lm9E1gJNPbZ7CrgLU+c3wNHGmNiP/d6DAYXq1IFpkyBbdvi/9pKqcPNnAkdO8KDD9px5QPSXurOlei9F6AZ8AtQ02f9NOB8r/tfAu38PL8fkAvkNm3aNPyvuhiVqi+5ROSzz9x5baXUQSUlIj17inz1lduRJA7cKNGXMsZUBz4A7hKRPyLZh4i8KiLtRKRd/fr1w3tyDAcXO/98yMpy57WVUlBcDI88YmeH+uADuPBCtyNKDo4memNMBWySf1tEPvSzyQbgGK/7TTzrYieGP9v++lc7+l3IIzvrT0alYqagALp2hdmzoUIFt6NJLo6NR2+MMcCbwFYRuSvANl2AAdiLsmcBz4pIh2D7dWw8+hCIwPLlcMopEU5EopSKSFERnHaaTfQjR2qi98et8ejPAzKBi40xizxLZ2NMf2NMf882/wXWAmuAccCtDsYTNWPscAha86JUfIjYzooVKthhSJ58Mswkr31ZAAeHKRaRb7HNJoNtI8BtTsXghPx8O1H4+vXpMw2ZUm7YuRP69bMzvM2eDU2bRrCT0r4spVWo3tfN0oimqjCdcoodtviLL9yORKnUlZcHHTpAtWrw3Xf2b5l8S++jR8PChYcOXtajh72fZs2cNdFH4O67YexYt6NQKsYSoJpDBDZtgoYN7eBk48fbPiwh8e2JXr483HsvtGx5sC9LYaEdzCzN6AxTEbj+erjkErejUCrGXK7myM+3VTVFRTBtGlx1VZg78B687JZb7FAHTz5p22MWF9tvjIoVHYk90WmJPgIVK9rl1VfdjkSpGAo0ymMcqjmysqB1azjxRNs+PmK+PdHbtLFJfs8eW7qfOjUt+7Jooo9QlSoweDCsW+d2JErFUJyH7Ni92+bg2rVh0iQYNSrKeV19By8bM8ZW4XjPHJeOfVkCdZlN1CWiQc0c8tBDIr17ux2FUjHkxJAdo0Ydvp+sLJk74E1p0ULkzTejf4nSfR4S85gxIsbYv/4eTzG4OQRCKhs4EFassBePlEp6Tg3Z4XORVLKyGdxlCX95uzePPQY33hiD2OHwnujFxbaOvrjY3k/jnumO9Yx1ips9Y/0pKbHt6WXUaEyH9of+1M3Oth8qr5mqlEpYo0fbpBzpZzjY8z3J/qdr/83x743g7X5f0enOk2nQIPaHka6C9Yx1vSom3CWRqm5KPfywyNO3rU6LMa+VCijIuO9btojcePoCac4a2T1omLtxpii06sZZvXrBo++14OdnP06bWeWVOkyAVjs51TM4pcV+6qyey+L7J1F1/LNp1+rFbZroQxWkM0mrVnDffXDdc+dS1O82nWREpS+vVjs/9nqAlUdlcNLmr/m0pCtP/bcl1R9/UIfrdoEm+lCVMf/rvffCPZcto/wrL4Q2q3w0vRAToAejSnOBPoM338wfL0zg/nO+5pyXMln49nJqLP+e9lP/rcN1uylQnU6iLq7W0QdreuZ57LMRi2TGDCm7jj6aeSzTZA5MlcD8fQZr1pSSmrXk7JO3S58+Ihvf/1Y/l3FEkDp61xN3uIvrF2OHDLH/tiFDDl3vaSucnS1y1FEieXliP+CjRgXeVzRtlnWKQuU2r8/g97Uuk3+0/EaKZ2bJjh0+2wQ7B1TMaKKPlRCT6+jRIqeeKrJtWwj7DPTFEYponqtUDGy4a7TcyBtydI0d8uabIgcOuB1R+gqW6LWOPlQ33wzdux/amaR7d7vex733wk032e7dQfl21w7n4lQ0z1UqSlu3QuGMWSx4bSGNzm3Oqgqnc+Mx2TpHQ4LStyUcvvMHBphP0Bib7GvWtA1wSjvmHSKaXog66bgKxsGL9bt3w4gR0LJ5Id/1eoquH/2Tx2dfQI0pr+tnMIFpog/VK68cHPmutI3w1Kl2fQCVKsE339jWZiUlPg9GM3G4TjqugimjhVikfv8dWrSApUthTr83uWjqXfoZTBaB6nSiXYDXgM3AsgCPdwR2AIs8y9BQ9puwF2MD+OMPkXPPFbnxRpGSEodjS0UBBsTSC3xliNHF+qIiO+jYq6/a+6tXxzBGFVO4VEf/BnBFGdt8IyJneJZhDsYSGxHUi9eoATNn2kkUjLEz6KgwOFQ6TXneww23bn344yFU5bz2Gpx0EowbB61a2XUtWjgQq3JeoG+AWCxAM4KX6KeFu0/XSvQxaLs+b57IRReJbN3qTIgpS5uShs/7f1arlkjNmiF9dnfvFvnqK3t7xIiDt1Xiw63mlSEk+gJgMfA5cEqQ/fQDcoHcpk2bOvefCibSKgSv5x04IHLXXSKnNd8p6we95FCgKSrapqTRVgElUxWSv0JJabIP8GW5fbtN7A0bitxwg1YzJqNETfQ1geqe252BH0PZp+t19OHyOelKvsySsdUekMvaFbgcWBKJRYk+2l9kydQbOdCX0iWXHPZluX27/fuvf4n87W8iy5bFMU4VUwmZ6P1sux6oV9Z2SZfoRfwmqqIikR077OQ3xcUOvGYylUCDCTXBhnK80X5hJHMVkk/s817MkRtusKu2b9cSfCpIyEQPHMXBiU86AL+U3g+2JGWiF/Fb9ZCfL3Lxxbag9fvvMX69ZCqBBhPqF1aoxxttFVAy9kb2/C+KZmaLiMi4f62UY4/4WZ64+Ue9XpRCXEn0wCRgE1AE5AF9gf5Af8/jA4Dlnjr674FzQ9lvUib6ICXB4mKRBx8UadpUZNeu+L1uSirreBOpRB/HX1z5Q5+VEf+3Rpo0EVm+3F5wLZqZnXy/7lRQrpXonViSLtGHWNL85Rf79/XXRTZujOHrJ2MJNBqBjjfR6uj79Tu8JUzNmnZ9DBw4YJfsbHsdtk8fkQULYrJrlaCCJXrtGeu0EHuxHnOMbWO/fr1t9jx+vJ829+F2bU+38XCCHW+0vYlj3Ru5d2/bsaJ7d9vTunt3e79378j257Fxox2ioEULmDULOnSAn36C11+HNm2i2rVKZoG+ARJ1SboSfQSWLBHp0EHk+ed9HginVJkqdfShSsbjzcoSqVrV/gKpWjXiWAsLRfbvF5k7V6R2bfujYO5cvcCabtASfXI57TSYMwf69rUFxuuus6WyQHNy+p2yMN3Gw0nW4y392XbYz7eynzZvHtx9NzRtCv/7H5x5Jvz6qx1+qUOHgGPuqXQU6BsgUZd0KNF727lTZPhwkbp1RW67zdPuOd3q3VNRaZ18rVr+e68GsHSpyObNIgsXipx4oshDD4msWhWHeFXCQ0v0yat6dXjwQVi5EmrVgorffcXy57LYeNfo9Kh3T1XvvmuL3FOn2qGmp0619wcNOuw9zZ/6LY9d8TWnngqdO8OKFfY6zqpV8PDD0LKlO4egkocm+kTnuQBbvz48dmk2VTKvZlbGI5z6ygD6tF3Ckp4Pa7JPRscfb5O7d1XT1KlwxhnINb1Y+GoOjzwCS/8zjz19b2dT5ea8/LK9WH/RRfY7QatmVKhKOywljXbt2klubq7bYcSP9yQjOTlQvjyMHMm2/3zIKysuYMZ72/iy9zg2ZQ6kUSM9+ZPV/v1QrhwsXAg9u+yjUsEmruqwkVt/uIvjPxjt/zqMUl6MMfNFpJ2/x7REH0+RzPzjfQF21y4YORImT6Z2twsYNAiyFtbG3D+QO++0TeqGDbOlvrTn4CxLsVJQAG+9BddcAw0b2u/xli1hxleV+fHfrzPm+/M5fsCVmuRV1DTRx1OkY6t7jy1+yy1+T/zJk2HSJMjPh3vusetyc+GPP2J8DMkiAcex37EDPv4Y7rgDtm2zpfePP4YuXWD1ajjnHDv9ZKvfsjEvp1H/B+W8QFdpE3VJ+lY3kXSjj7Dr/c0328Yc114r8sEH0YWdlFweAmLvXpE5c+ztF18UqV5dpFMnkccfFykINHhpMvYHUAkBHQIhwYTTPDLKE3/zZpFXXrFjjYuIPPaYyJNPivzwQ4SxJ5s4NkXdts3+HTdO5KyzbB+os84S2bPHPrZ3bwg7SZVRR1XcaaJPJOGWMmN84s+YYXtONmokcuGFdt2qVSIbNkS0u8TmYIn+wAH7d/x4kauvFmncWOTMM+262bNFZs1yYJA6pYLQRJ8oEuhn+YEDIj//bG8/+qhInToizZuL3HST7U6/b1+Sd6GP0f+6pMT+P4qKbMe1q68WadFCpHdv+/hLL4m89ZbImjVJ/v9SSS9YoteLsfGUQN30jzjCdp0HeOABexH300/h4ouhYkV48UWoXRvOPx9uvRXy8mwTwJ074x5qZML8X+/fbzsg7dkDv/wCt91m26vXqWOvh5YrB0VF0KOHbe7+1lv2ef37Q2ambRavTVtVotJ29CqgggJYuhSWLIHrr7ctQzp1gqOOglat4K674NJLbcuRY46BZs3sl0OiJrxff7VNTzdutMd2660wcaLtebxpk/3i++ADe3yTJsHJJ9seqA0auB25UmUL1o6+fLyDUcmjbl3oOG80Hdu3h3oZ1Ktnm2uumfg9K7M20rRpT4qK7BC469fDunW2qeA778C//23X1a1rE+cDD8CaNXabOnWgalX7xVChgm12WKkSVK5s+4N5E4EDB+z6fftg715bst650ybm7dth9my7jx074IILbHL++99tQt+40Q7PO3Gi7WOwciU0agRNmkBJCVxyCZx9Nhx7rI2l1J13xu//rJTTNNGr4Erbo3uqQcp9nU3Le3vRcvJkONFu8tFH9q8IFBba23/5C6xda0vOu3fbdYsX22qQggKbtCdOhGrV4NxzbdXJvn1w333w+OO249CaNTYZn3aa/VUxdCiMG2eTfo0a8NVXsHmz/aI58kg7FlDbtrZa6rLLoH59OPpoaNzYvv64cYcfXqNGTv7zlEoMjlXdGGNeA7oCm0XkVD+PG+AZoDOwB+gjIgvK2q9W3cTZ6NF/DrvALbfYTD14MBQXw8CBMX+5khKbqPfvty97xBGJWxWkVCJxawiEN4Argjx+JdDCs/QDXnIwFhWp9u1tkr/yStsz98or7X2Hepge4flEVqpkL4Bqklcqeo4lehH5GtgaZJOrgLc8LYO+B440xugP6USTkWFL8BMn2grwiRPtfR1/Ramk4WbzysbAr1738zzrDmOM6WeMyTXG5Obn58clOOWRnW1L8H/7G3zzjf07cqSOv6JUEkmKdvQi8qqItBORdvXr13c7nNhK9FEWc3JsCf7zz+0gW59/bu8n+hR9Tkn090spP9xM9BuAY7zuN/GsSy8JOMriIUrr6CdPtu0TJ092tI4+4SX6+6WUP4G6zMZiAZoBywI81gX4HDDA2cC8UPaZ1EMgBBKvURYjGTfnyitFxow5dN2YMXZ9MotmDCGXR8VUyh/cGOsGmARsAoqw9e99gf5Af8/jBngB+AlYCrQLZb8pmehF4jPKYiTjvyTQ+DwxFe1x6QTtKsG4kuidWlIy0cezhBjH8fATXqTHlar/D5XUNNEnMjdKzJGURlO1BBvucaXqLxyV9IIl+qRodZPS4j2iZXa27d0azjR1kTwnGURyXAk0AqlSIQv0DZCoS8qV6ONJ6+gPStXjUmkLLdErILLSaKqWYFP1uJTyQ8ejV0qpFODWoGZKKaUSgCZ6pZRKcZrolVIqxWmiV0qpFKeJXimlUlzStboxxuQDP0e5m3rAlhiEk0hS8ZggNY9Ljyk5JNsxHSsifsdxT7pEHwvGmNxAzZCSVSoeE6TmcekxJYdUOiatulFKqRSniV4ppVJcuib6V90OwAGpeEyQmselx5QcUuaY0rKOXiml0km6luiVUiptaKJXSqkUl1aJ3hhzhTHmB2PMGmPMILfjiQVjzGvGmM3GmGVuxxIrxphjjDHZxpgVxpjlxpg73Y4pWsaYysaYecaYxZ5jesTtmGLFGFPOGLPQGDPN7VhixRiz3hiz1BizyBiT9MPlpk0dvTGmHLAa6ISdrDwHuE5EVrgaWJSMMRcCu4C3RORUt+OJBWNMI6CRiCwwxtQA5gPdk/m9MsYYoJqI7DLGVAC+Be4Uke9dDi1qxph7gHZATRHp6nY8sWCMWQ+0E5Fk6jAVUDqV6DsAa0RkrYgUAu8CV7kcU9RE5Gtgq9txxJKIbBKRBZ7bO4GVQGN3o4qOZxKgXZ67FTxL0peyjDFNgC7AeLdjUYGlU6JvDPzqdT+PJE8e6cAY0wxoA8x1OZSoeao4FgGbgZkikvTHBDwNDARKXI4j1gSYYYyZb4zp53Yw0UqnRK+SjDGmOvABcJeI/OF2PNESkQMicgbQBOhgjEnqqjZjTFdgs4jMdzsWB5wvImcCVwK3eapIk1Y6JfoNwDFe95t41qkE5KnH/gB4W0Q+dDueWBKR7UA2cIXLoUTrPKCbpz77XeBiY8xEd0OKDRHZ4Pm7GZiKrfpNWumU6HOAFsaY44wxFYHewCcux6T88Fy4/A+wUkTGuh1PLBhj6htjjvTcroJtFLDK1aCiJCKDRaSJiDTDnk9ZIvI3l8OKmjGmmqcRAMaYasBlQFK3akubRC8ixcAAYDr24t5kEVnublTRM8ZMAr4DWhpj8owxfd2OKQbOAzKxJcRFnqWz20FFqRGQbYxZgi10zBSRlGmOmGIaAt8aYxYD84DPROR/LscUlbRpXqmUUukqbUr0SimVrjTRK6VUitNEr5RSKU4TvVJKpThN9EopleI00SulVIrTRK+UUilOE71SYTLGtDfGLPGMMV/NM758Uo9bo1KbdphSKgLGmEeBykAVIE9ERrocklIBaaJXKgKe8ZJygH3AuSJywOWQlApIq26UikxdoDpQA1uyVyphaYleqQgYYz7BDs17HHbawwEuh6RUQOXdDkCpZGOMuREoEpF3PHMRzzHGXCwiWW7HppQ/WqJXSqkUp3X0SimV4jTRK6VUitNEr5RSKU4TvVJKpThN9EopleI00SulVIrTRK+UUinu/wHs465PT06pmwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot of fit\n",
    "fig,ax = plt.subplots()\n",
    "ax.plot(x,y,'xr',label='Data')\n",
    "ax.plot(*line([i.n for i in params]),'--b',lw=1,label='Fit')\n",
    "ax.legend()\n",
    "ax.set_xlabel('x')\n",
    "ax.set_ylabel('y')\n",
    "ax.set_aspect('equal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Confidence interval/area\n",
    "\n",
    "To visualize the accuracy of the fit, the region in which data points may validly lie within a given range of statistical confidence has to be establish, which is referred to as the confidence interval/area. This explanation is based on O’Leary & Zsombor-Murray, 2004.\n",
    "\n",
    "The singular value $\\sigma$ is the two-norm distance between the fit and the original data \n",
    "$$\\sigma = \\lVert\\bf{r}\\rVert_2 = \\sqrt{\\sum_{i=1}^N r_i^2}$$\n",
    "with the fit having a residual standard deviation of $\\delta_r = \\sigma/\\sqrt{N}$. The fitting procedure shown earlier gives a mean residual $\\bar{r}$ of zero. Assuming the residual has a normal distribution, the conic equation can be rewritten as \n",
    "\n",
    "$$ax^2+bxy+cy^2+dx+ey+f = \\aleph(0,\\delta_r)$$ \n",
    "\n",
    "Note that $\\aleph$ here refers to the normal distribution. In usual applications, the confidence interval is defined as a multiple $\\mu$ of the standard deviation. This gives an expression for the confidence interval for the ellipse fit as\n",
    "$$ax^2+bxy+cy^2+dx+ey+f \\pm \\frac{\\mu\\sigma}{\\sqrt{N}} = 0$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 398,
   "metadata": {},
   "outputs": [],
   "source": [
    "# obtaining the confidence interval/area for ellipse fit\n",
    "def confidence_area(x,y,coeffs,f=1): # f here is the multiple of sigma\n",
    "    c = coeffs\n",
    "    res = c[0]*x**2 + c[1]*x*y + c[2]*y**2 + c[3]*x + c[4]*y + c[5]\n",
    "    c_up = np.array([c[0],c[1],c[2],c[3],c[4],c[5] + f*np.std(res)])\n",
    "    c_do = np.array([c[0],c[1],c[2],c[3],c[4],c[5] - f*np.std(res)])\n",
    "    if convert(c_do) > convert(c_up):\n",
    "        c_do, c_up = c_up,c_do\n",
    "    return c_up,c_do"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 399,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([-0.38132347,  0.35651649, -0.8529293 ,  1.1310537 ,  3.72028399,\n",
       "        -3.65143635]),\n",
       " array([-0.38132347,  0.35651649, -0.8529293 ,  1.1310537 ,  3.72028399,\n",
       "        -5.36797898]))"
      ]
     },
     "execution_count": 399,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c_up, c_do = confidence_area(x,y,[i.n for i in params],f=2)\n",
    "c_up, c_do"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 400,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAEJCAYAAAC0U81tAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABJcklEQVR4nO2deXgT5fbHv2+6b1BKWyiUHUE2WSyoiHjjCoKKohWvyNUrIni9iogK/gQVVLQXlOtVcUERcEFEQUQFUeIu2iKLIIIsFcrSlZYudEve3x+n06ZpkmaZySTp+TzPPEnaycybZOY7Z857FiGlBMMwDBM4GPQeAMMwDOMeLNwMwzABBgs3wzBMgMHCzTAME2CwcDMMwwQYLNwMwzABRqiWGxdCZAMoBWAGUCulTNNyfwzDMC0BTYW7DqOUssAH+2EYhmkR+EK4XSYxMVF27dpV72EwDMPozrZt2wqklEn2/qe1cEsAXwghJIBXpZSvOVu5a9euyMrK0nhIDMMw/o8Q4i9H/9NauEdIKY8JIZIBbBZC/CGl/NZmcFMATAGAzp07azwchmGYwEfTqBIp5bG6xzwAawEMs7POa1LKNCllWlKS3bsChmEYxgrNhFsIESOEiFOeA7gCwG6t9scwDNNS0NJV0g7AWiGEsp93pZQb3d1ITU0NcnJyUFlZqfb4/IrIyEikpqYiLCxM76EwDOPnaCbcUspDAAZ6u52cnBzExcWha9euqLsIBB1SShQWFiInJwfdunXTezgMw/g5fp85WVlZibZt2wataAOAEAJt27YN+rsKhmHUwe+FG0BQi7ZCS/iMDMOoQ0AIN8MwDNMAC7cLhISEYNCgQejXrx8GDhyIRYsWwWKxOH1PdnY23n33XR+NkGGYlkRwCXdGBmAyNf6byUR/94KoqCjs2LEDe/bswebNm/H555/jiSeecPoeFm6GYbQiuIR76FAgPb1BvE0mej10qGq7SE5OxmuvvYYXX3wRUkpkZ2fjoosuwpAhQzBkyBD8+OOPAIBZs2bhu+++w6BBg/D88887XI9hGMZd/KrIlNcYjcDq1STW06YBS5bQa6NR1d10794dZrMZeXl5SE5OxubNmxEZGYk///wTN998M7KysvDMM89g4cKF2LBhAwCgoqLC7noMY0ttbcMipeMFsP93g4GWkBBarJ8rr3kuPLAJLuEGSKSnTQPmzwfmzFFdtG2pqanBPffcgx07diAkJAT79+/3aj0m+DCbgerqxoJcWwvU1DR9bTY3iLKWKGKuPIaFARERtISHNzw3BNc9edAQfMJtMpGlPWcOPRqNqov3oUOHEBISguTkZDzxxBNo164ddu7cCYvFgsjISLvvef75511ajwlMLBagqgqorGz6WFur9+iaYrHQ0hyhoU3FPDwciIykR0Yfgku4FZ+24h4xGhu/VoH8/HxMnToV99xzD4QQKCkpQWpqKgwGA5YvXw6z2QwAiIuLQ2lpaf37HK3HBA5SOhbnmhq9R6cNyt1AeXnT/4WGAtHRtMTE0COLuW8ILuHOzGws0orPOzPTK+E+c+YMBg0ahJqaGoSGhuLWW2/FjBkzAAB33303xo8fjxUrVmDUqFGIiYkBAJxzzjkICQnBwIEDcdtttzlcj/FfKiuBigoSrfJy4MwZ16zUlkJtLXD6NC0KYWEk4rGxtERHsz9dC4T0hUPNRdLS0qTthN3evXvRp08fnUbkW1rSZ/U3qqsbi3RFBfmbA412yzNQ3m8oytIaDJXYLBNi9mQi9x8P+Xw8BgOJd6tWQHw8EBXl8yEELEKIbY769AaXxc0wLmA2Nxbo8vLgcXWU9xuK7rPSceiZ1ShLMyI2y1T/Wg8sFqCsjJbjx8lHHh9PS2ysLkMKCli4maBHShJo5ba+vNw3kRt6UJZmxKFnVqP7rHTk3zANSWuW1Iu4P1BVBeTm0hIaCrRuDbRpA8TFcQSLO7BwM0FJVVWDUJeWBqbbw1PK0ozIv2EaOiydj+OT56gq2mq6YmprgcJCWgwGcqe0aUPWOIu4c/jrYYICsxk4dQo4cgTYvZuWI0eA4mL/FO12yzMQm9W4PENslgntlntXnkHZTtKaJTg+eQ6S1ixpsh9vUFwxyjYVV0x5P++yky0W+q0OHwZ27QKOHqWLL2MfFm4mYCkvJ7/pH38AO3cChw4B+fmBccJrJYDWPu0TU+fVu03UEm9rV0zKK3Mb+dOt8ebCZDYDeXl08T1woHHUCkOwcDP+j1XxsLIyssYOvWHCqdkZOHEiMH3Wrgqgu8TsyWy0HWU/MXsy1Rh2/TYVV0z+DdPsjlmtC1NJCfDnn8CePSTm/nj3pAcs3C6glHVVluzsbAwfPhwAVwHUGimB8r5DYb4hHQeXmrBvH1DxmQmdZnpvneqNKwLoLrn/eKjJdsrSjKqGArriilH7wlRZSRfs335jNwrAk5MuoZR1tUap7qcI99///ncdRha8lJUBRUXkt65NMSL26dXo/mA6ovwwUsJTbAWwNM3o95/J2hVTlmZEaZrRoShrMUmquFHy8oC2bYGOHSnpp6XBFreHxNYFodqWb2XcwKZ+ekUFkL/ahBP3Z2DfPvJXK3U+tLBOHaHlxKH19jo8fCt+fPhj7L+FfNF5MzOwa+UOfPMN8OuvtN7PPwNr1wLr1gEff0zfR04O8NlnwMaNwHffAeLllxD1iwlHjlCERmUlEJPp/XjtfQ8Jm1ah2Hi9S64YLSdJAfqsu3cDJ060vIxWtrhdQEl5B4Bu3bph7dq19f+zLd/KuMHQoZDp6Tj1ymoc721E2A+Ok0V8aZ16k8QiJflli4uBrl1JeH/6ie4ciouB2bOBkyeB+/59IaotR9B6kQHTpgHXXmvEgx0HIezDEohuQJ8+wJAhJNK//95QsvWqq4CCAtpmbS2l4SeedwGunHUjRoXvRpk5CuWlZkwwnMQjLwzFU0/RvpOSyDqdOJFiqM1moH1752F39r6HeNNHTb6HMju/hTuWuTdYLDRBXVBAny8hQbVN+zUBJ9yPPw5YN59RMuTTrBJDH3uM1uvQga7GAJ0E27YBU6YAr7/esO6xY7SeM+y5ShjvOH0aKOhiRO2Tq9FtcjoSnLhAfCUCCs0lsUhJdwOHD9OSmAhcdhkwcyYJdWgocPbZVJxSSopN7t6dHlu1olv8TzeHIyamcR2Pl95pA6BNo7GMH0+LNYMG0dLAEBSdtRQHZ3WuH++Bp2m842PpGM/Pp7R+ANjxny/wv50jcaoiEp06Ac89B7Ta8yN2bClC4j/GomdP+gzeJPM4myTV4jerrqbfIi8PSE0N/qxMrlXiArGxsSgrK7P7t6+//lo1i9sfPquW1NaSZVRQ0HhyKeWVufV+0BNT5zV5n171N1JemYu4pc/jk1Ev4Zfek3DoEHDllSSa115LFnW3bsCFFwIXXURGQo8Ni2EYPFCXWiHNfY/W4+k+Kx2/Pf4h9iaOxDmnvkbtI3Mw5+w12FvYDsePk+FjNAKbNgEX73oBI9ZMx8nJjzrdrj/Rpg1Z4BEReo/Ec7hWiYbYlm9lmlJaShZfcXHTsD1XXCD2BM/e7bm3lJVR2NnevcD+H/Lw4oEN2H7NQrz+WTd0qT6KAed3QqdOVIt606am709JAQyDB+pSK8QdV5Ji/Q6YNR7tFUv6P6vxRFo7AA0lXMvLgZ8+zsWb26/FqfA7MfedJzEmzYTDXYxITPTvqn+nTpHLqlMnuiMKNli4vcS2fOv999+v95D8AiWd2VlCjK9dINZISZmV27eTD3n2bBLjjRuBcxKO4pZ9T6LkqcVof9FIvHKVCd1nDcGh9NWI+TLTqfWvR60QT75HZxEfSsXh2CwTPj6YjkMvr8Zf3bogJHMsuj18DUZHZuOMiMGwYcCll9Idh6doeTdlsQB//UWGQ+fO1OknWOCoEhewdZNY/y0sLAxbtmzBzp07WbRBVquStpyT4zze1hfJIgoWC2XhffstvV6wALj7bpr36NOHLjTjx9P8x5N938N5z02AuGhkk3G5kljiywgYwLPv0ZWID+vttm0LxI+6AIefXY1vbnwJL74I9O5NE50A8NJLwKef0u/vDlplkFpTVER3URUVqm1Sd9jH7UcE8mctKaFoCXdPXC2pqaEY3+efBzZsoInBESOABx6gkLmICM9u9xVxcWRRN/d/tfDUWrW10G1fe8LGjcAXX9CFcNw44P776a7Gle/XV9+XEDRxmZys+qY1gX3cjCZISb7EkycpLM0f2LED+PFHWtq2Bf77X5pk+/vfgXbtGtbzpuWnMzeDL90/noYtahHxMWoULWVldDwAwF13UXTN9dfTRK4jtKxmaI2UlHVZWkrjCWTXCbtKGNewSpaxWMh3nb3MhIrHM3QV7TNngC+/pAWgTnVSklW9aBH9bdCgxqLtLc7cDL50/3iaVq5lWnxsLNCzJz1/7DHqQTllCrB4sZP3aJyoY0txMc1r2OujGSiwxc24hlWyTM5ZRkT8ZN+680XonpKM8sgjlIgyYABZdQDw9NOO36fG2JqzqH0VAWO9bV9Yq57QsSNwzz0k3EVF5Lq66y5KIrrmGqDTexmwhIYiZdmC+u+vNi4ePaePxYHFGzT9LNXVwP79dJGJi9NsN5rBFjfTLLW1wLFeRhx8ejXiJqcj8X+OrTutJptKSynl+777KGHEYKCTf/164MUXgUsuaby+vXRtS2goek6/2qux+dKidgVfW6ue0Om9DPTMMSEsjHzf330HXD+mEt9mRaHjkrk4cfvseldPyrIFODZtvk++T2XCOhDLxvLkpB/hb5+1upqiBgoKGmpBuJLkodZkk9JFvFUrSnjp25fCz0aMaD4zztEE3InbZyNl2QK/bOvlLlpMMmqB9bhi9mTCEhqK469/hoIZTyG5dTW2zv4YE8/ZhbgD23UZu8EA9OhBx5k/wZOTjFuYzZQFmJfXOGHG1SQPb27fpaRb2A0bKErh6qvpdnvtWvKXuoqzeOrQ0mK/dC24i6/Tyj3F+rcoGT4abT9/G5i+EGHXnIfSL37E+prReC3rNsy9KhM9dBi3Ynn36EE9MAMBdpUw9UhJFvbu3fRoK9qudlbx5Pa9uJhOHimB//yHkkBef51EG3BPtBXsxVMHgmvBVXxRe1stlN8i8bOVKBw9ESnLFiDllbkY8dRofB49HtMv/Q3/t3Ekyjb/qMv4pAQOHqTjMBAIOIt72zZ1t3fuuc2vc/ToUUyaNAm5ubkQQmDKlCm477771B2IzhQVUZU1Rwkzrlp37oTDSUkXiTVrgG++ocp1PXsCS5eq85lsRbo2Lr7RRJgvMzVbOra/Rcnw0eiwdD7MkdE4uHgDBqcZsf4XE3o+ko6pW7Yhtk9n3HKLb0P2pKT2d926Ua0TfybghFsPQkNDsWjRIgwZMgSlpaU499xzcfnll6Nv3756D81rysoow7G50ChXoyVcEfgzZyie9qyzgBUrKCrk/vupu7da2LuA9Jw+Fsemzfd710KwYftb1MbFo9PimSjvPRgROQfr16sYRr/HHT9+jn9/fxe+/54qgaak+G6sUlLmb0iI//m8rWFXiQukpKRgyJAhAKioVJ8+fXDs2DEAwM6dOzFy5Ej07dsXBoMBQgjMnTtXz+G6RGUluSb27VM3ntXZ7XtODrlBxo6lxgBC0OtJk9QVbcD+BeTA4g0wKJ0ZbMbGaIf1b6FEjhydvhCnrpiAg4vWNXK5laUZEXHvXViyBBg+HPj8c9+PVxHvmhrf79tV2OJ2k+zsbGzfvh3nnXceKisrcdNNN2HFihUYNmwY5syZg8rKSjxhXTDcz6ipIZdIYaHvGuwWFFCFtldfpeL977xDj1ri63hqd1FSwYVovAD0u5jNgdcA2RHWv4XtBRWA3buekBDgttvo+c6d5E575BEgKso3Y66tJbdJr17+WQVRc+EWQoQAyAJwTEo5Vuv9aUlZWRnGjx+PxYsXo1WrVtiwYQOGDBmCYcOGAaBKgRs3boTww19aSkpFPnnSN22eLBbg+++Bt98md8w77wDz53u+Pb1qcruCwUB1T8LCqAGB8qgs1q8NBteFwGKhxWxuWCwWEpWqKgrXrKqixZ+tQ2s8uaD27k3f2513Ugy/r2qNlJWRkdOxo2/25w6+sLjvA7AXgB97jJqnpqYG48ePxy233ILr69L0du/ejQEDBtSv8+uvv9a7VPyJ0lIqYVpZqf2+lMJCL70EbN0K3HordYfx9lrmTTsxNRCCxDkykh6tn3sS8eIKBgMtoS6cpRYLCXl1Nf3OFRW0+EsNGW+IjKTGDsuX0+PLL/tu3ydPUs6Av4UJaircQohUAGMAPAVghpb70hIpJe644w706dMHM2Y0fIy2bdtiy5YtAID9+/fjo48+qu/+7g+YzTTxWFCg/b5qaoBPPiEL+803gcmTKZRPrZsPX9a5NhiA6GgKSVQew8P985ZZwWAggYuMbDypZrGQgJeXkwVZVtbQgDmQEIJcJzffTBejX3/1rg64Oxw+TMlfWl2gPUFri3sxgIcAqFYNwJXwPbX54YcfsHLlSgwYMKC+afDTTz+Nm2++GevXr0f//v2RmJiI9957D23btvX9AO1QVESRG744SXftAh59lIrVP/64+hONClrU5RCC/KbWIh0Z6d8i7Q4GA1mMsbENhbbOnKF45eLiwKtRHREBVPz3dTz9ya2YNj0SY+ucr1q6zcxm8nf37u0/x4Vmwi2EGAsgT0q5TQjxNyfrTQEwBQA6d+6s1XC8YsSIEXBUGuCTTz7x8WicU11NbpGSEvW2ac+/HPHz19j0UQW6TrsKHToA8+bZNrBVH7U6vUdFkVXaqhUJmrNO58FIVBQtKSl0vCgiHigd+JIv7InN64y47PmvIUQEbmqnvdusvJyyiZtrLO4rtDxkLwRwjRAiG8AqAJcIId62XUlK+ZqUMk1KmZaUlKThcIIbJetxzx51RRtoXDhKSuDbV/Zg3L+74IOcC1BbSxEjvhBtVzM3bQkLo9rc3boB55xDt72pqSTcLU20bQkPp8m+Xr2AgQOBLl3ozsOfKUszIvw/T2OzxYiwTZ/4LIkqN9d/JoE1s7illLMBzAaAOot7ppRyolb7a8lUVFBvPa1uexX/coeHb8Wxa6fhs3eHY8698eg/sZs2O7SDu3U54uJoQqlVK9+FkAU6oaF0EU5MJAszL48aZfhjWGJZmhEpN12G2UuvwVuXv43cYiMu03ifFgtZ3f7gGAiIOG4ppV+G2KmJJ1UaLRYKV7ItBqU2R44A/3vfiPjET/HRikFYPnkOTky0XxlQK5oLIxOCRDo+nhZXIjEYx8TE0B1Kaio1zcjP969JTWu3Wb9VSzHjpxuRmhqOs8/Wdr8FBTRXEBGh7X6awyc3ilLKrz2N4Y6MjERhYaFHwhYoSClRWFiISDf6aVVUUANU22JQarN0KXD77cDg1gfxVsFYvyvO1KoV3d4PHEh1ThITWbTVJCyM/LoDBgCdOvlHuy9bt1mrhXOxxHIXHrq3EoWF2u5bSjKW9MbvD/HU1FTk5OQgPz9f76FoSmRkJFJTU11aNy+Pwvy0EmyzGdiyhWpfDxsGTOrxPdKeug6HnvWP4kwxMUBCAhUCCgtzsmJGBjB0KDWdVDCZgMxM4CFOc3cHg4F84QkJDSV/9cKe22zQc8D9H5hgsYzWfP9FRWR16zkX4PeNFJgGamuB7Gz1Jx+t2bMHeOYZColbuJD8xP6QtWgw0ARjYqIbJ4zJBKSnUyNKo7Hpa8ZjKivJeNDyWPSEU6eAr78GrrtO2/20akVF0rTEWSMFFu4AobRU+8I3v/9OVfr+/W9gzBj/iFmNjgaSksi69ug2XRHradOAJUtYtFXm9GmaGK+u1nskRGkpMGEC8N8LVmHQqHaaGhu9emnbr9KZcLfwYCj/R/Gp7d+vjWhLSd1mNm4E+vQBPvyQqvfpKdoGA1nWffrQkpjohW/VaCTRnj+fHlm0VaVVK/qNtEq6cpe4OGD2bOCB765F8sO3q9771Bo93UUs3H6M0on6xInm17XXHDc2y4R2yzMcvufECeBf/wJWraIQJyGa7+WoJeHhNAF2zjkqxhObTGRpz5lDjyb/mFQNJkJDqe1Xp07+cZc2YgQw/JIofPLPteg+Kx0przhubg14du4A5CbSK9KGhdtPKS4m10VZmWvre9Jd/cUXafLxrbcoKUUvIiOBrl2B/v1pAky1yAVrn/a8efSYns7irRHJycDZZ+sfKgcADz8M9LtlMP669t5Grevs4cm5A9DdalGR6kN3CfZx+xkWC036eBJE40p39RMngMWLKagiIUFfCyk6mupyx8drNA6OKtEFs5n6N+qdQv/rWzuxckk53rttI5I/dF6UzJVzxx7R0eQq0gL2cQcINTXkGvE08tFec1wFKcl/feutdKC1bq2faMfF0Yx8nz406ajZOB56qKlP22hk0daYkBCKqXfW+stT94SrxGaZMHHllTiWNBBrBzVfIsHZueOMigrHfVq1hIXbTzhzBvjjD0o19vSgdtbBPDcX2LSJutDcdps+SSqxsVRhrVcv/+7nx3iPwUB+b0e/s6fuCVeJ2ZOJI8++hynTY/DSS0DpuQ0lEqxRzrXYLBPavf0cCq66FcmrXkDnp+6qH5crvm5fw64SP6CkhEL9zGZ6bdtc1fa1PRy9Z8WkzfgqfxAeeMCHH8iGyEjqIuIvkQeM77BYqCSqPXHz1D3h7v6/+AK4/HL7cyexWSb0eGAcIASO3zkXHV6fB1FbAxkahuN3zkXKsgXNjisujowRtWFXiR+Tm0tNexXRBho3DWhuRlzBNpuseJARdw/LwqNv9MDFF2v9KewTGkrRKn37smi3VBTL214HGU/cE+7ejRoMJNq7d9vfTlmaEaeumABIieh9OyCqKiFDw1A88hp0XDLHpYtJWVnj89cXsHDrhJRUvCknx/7/3T2obburb9gA7C3vgpXr4pBm95rtHG98kAYD1Xru35+SZ/whRIzRDyEoasi2g4wz154jPHGxmM3AzJnUWMTedo7836sovvhaJH62EqcuT0fehHuR+NlK5E58wKWLiZS+n4hl4dYBsxn480/nk5CeHNQA8NVXwE8/AddcQ9Ejbdp4NkZPfZCJiSTYHTr4R0Eixj8IDQW6d2+4iHtaX92Tu9HwcMoEXreOXivGh7KdLnMnoe1nK1Heewjiv/kYyatecPu880U/V2tYuH1MVRVNQjq7QntyUJvN1KB38WJyS7jTTdwe7p4gUVEUw9ulSzOFn5gWS0wMlYkFnNdXbw5PXCzjxtFdaE1Ng1ECACXDRyPxs5WQ4REovOoWOmmkRKnV8e+KePtauP2+OmAwUVpKEzXNZVu52zQAAJ5/nqz4FSs8t7JtcaXHoxBkXbdrxy4RpnmSk8kn3Fx9dWd40sKua1fKGfj1V+C88+h86jHzOhjOVMASHgFLWASisvfh4MK1AFBf06S5806BhTtIOXWKIkdcCeJx56A+fJhOhttuc62BgDuV/po7QeLiaPLRjTLiDIMuXag4lScTerbRU+6UGB43jsJuFURNNQzmGhy/fU79doqunNDoXHP1YsKukiDEHdF2hy1bgDvvpNR4VxsIuOq7duauCQkhC6ZXLxZtxn1CQmjy2hO8cbFcdx3wt7/R84RNqyDDwuuNEgAub8ceZrNv+1FyHLfGFBVRDW01v+Z2yzPw35zxWPljD2RkAP36uVey0pX4WUeWeeKhTLR68iH2YzNeYbFQiJ6vm+/edx/w1NgfMfzZa93Kk3CFs85SN7GM47h1orBQfUtbSrKa+25ahNXTf6gXbXeyzlyZ3LENLzQYgMQbjGj7LIs24z0GA82NaI1tWGt0NLDrg/0oNl7vkdXuDF+6S1i4NaKggCxtNamspDIbX8OI4c/diAueHedySJQ17oYaRkVRXZG2bdX4FAxDtG2rvavN1jU4KO4A/txdhaIrJzRaryzN6HWDBRbuACc/n7qCqMmpU8DUqXSgDxzoeVEcd0MNlVKd7Mtm1EYIivTQEtuw1ks2PYzCC8Zo0ivVl9mTLNwqk59PGZFq88wzVKF03jyKk/Y0QcfVyR3r4vgGPkoYjWjTRvvjy9rIOefmfpixyLWm3P4MT06qSF5e47RaNTh6lA7u0NAGq9eTIlTuEBsLdOvWNEWZYbQgO5vmg9TEenJdOT9Kho9Gmy0fYspFe3HrI51V7xeZkEDnjVrw5KQPyM1VX7QPHQKmTAGyshq7KrwJiWqOpCQK82PRZnyFWglj1ii+7aR3nkP3Wek4cftstP7xcxybNh/bt5xC/sbANRABTsBRhYICx8WimsNR2N2Rb7Jx9+bbce+9DbGnCt5knTmjUyfyaTOML4mLI3eJxaLeNhVjpuf0q3HqkusblWeN2liEqr1fAvCg+pqfwBa3l5SUeOfTdpQQkxV2PmbMAK66SqWBOsFgoI4lLNqMHhgM2jTWKEszInfiDCR+trLRBH50xwQcPT9d/R36ELa4vaC8nNwZ3kwTWM96598wDTmrfsALN32NK6Zp1MjOhvBwEu2oKJ/sjmHsEh1NDbLVxFHJhkceCfwoKba4PaSqihogqHF7p8x6ly59H1eb10Ge3c/7jbpATAyF+rFoM3qjtpA6C3stKgKOHVN3f76GhdsDamupEl9zVf5cJTbLhKrV63BZ7FY8hf/D1XGuhfZ5Q3w8TUJyFiTjD6htPDibwN+4Edi8Wd39Ab6tjsmuEjeREjh4UL3OzoplsG3uWkypaIOLkq9TNbTPHm3aUNgSl2Fl/IWICHW352wCv+q/2kSy+NIIYovbTf76i+oJu4ujVmDRy17ExAE7UDNsBK66St3QPnsosaYs2ow/IYTv/M5nzjTelzdt+qzxpd+chdsNTp70PFHAXvRIp4f/jgklr6MmsWOjW0WlboJaB5RCQgKVY2XRZvwRta1uwL4oX9/1V4w+trT+tadt+mxh4fZDiou9m9Cw1wrszkGZQEICZs2yL6ZqHVAAFfRh0WZaGvbOoWtfH4P2I3rUr+NJH0t7+FK42cftAlVV6lT6s20FNrh3KqYMddxU1zZU0FHt7OZITKRONSzaTEvD3jk0sv1+TDS3wfk26zXXps8ZoaG+bY7NFnczSEk1tV2t/OXMvaHEla6+6k18/s4pXB1nQmys8+15WgVQISGBRZtp2dieQ3+dblPftFjB06JtCr6OC2fhboYTJyjRxlUcuTcsoaHoPisdXz+wHtO23o7Qu/7pUgdpbw6oVq3YPcIw1udQ6w+WoiDf0qicrLulju2hhX/eGSzcTigrI+F26z0O/GWG2lrsfmINpi67AFOmAP0nDm42esSbAyoqCujenUWbadnYnkO7576P2w3LEb+j4RxSo2ibry1uLuvqALOZmvBWV3v2/pRX5tb7y05MnQeA6pp8/jkwwar5hrOu6wBc7shuTVgYdazh5BrGr8nIoCLzRjq+9+8HpMn13qmu4Oz8UmsfANWuj49XbXMAuKyrR/z1l+eibc+9sXUrxY9OaNwxyWnkiG3fR6D5FktKwSgWbcbvGToUSE8HTHTsh37XNGrK25BY23No7Vpg5VHv25TZ0txcldpoJtxCiEghxC9CiJ1CiD1CiCe02pfaFBZSqzBPsOfekA8+hDmzqu0W0VErFEmha1cq2MMwfo/RCKxeDaSnQ86Zi04zmx77aobEAsAPP6h/fsTEUFSJL9HS4q4CcImUciCAQQBGCSHOd/4W/amq8q4hgq2/7NRAI9ITvsJ9g7/D2Wfbf4+3kSMK7dtrk8rLMJphNALTpkE8af/YV9OwkRLYuZN6tqpJ69bqbs8VNBNuSSjJ4WF1i/841O3gbuifPWxvzQ4eBDr2aYUxz13q8D3ehiIBdNXv0MGjITOMfphMwJIlODPT8bGvlmFz6hSQmqp+g+KgEm4AEEKECCF2AMgDsFlK+bOW+/OWkyfdC/1rjpISoHdv4MknHUd3qBGKZDBw/REmADGZyMe9ejWKZzg+9tUwbADKaVi2TI2BNxAWpo9rUlPhllKapZSDAKQCGCaE6G+7jhBiihAiSwiRlZ+fr+VwnFJdTcKtFrW1wN13A1u3Ol9PjVCkzp19H0fKMF6TmUk+bqMRZ87YP/bVMGwU3nqLauiriR7WNuDDcEAhxFwAFVLKhY7W0TMc8PBhoKhIve2tWAH8/DPw4ovaWsJt2lC8NsMEKorv2Z6LUq1wvpoa4Ior6DqRlKTGqAktwgAVnIUDajYXKoRIAlAjpSwWQkQBuBzAs1rtzxvKy9UV7ZMngeXLadFStMPDydpmmECmvNzxvJJajbEzM8mdqKZoC6FNr0xX0DKIJQXAciFECMgls1pKuUHD/XmMN1Ek9khOBl55BU3qIahNly6+D0NiGLVRu9ekPXbuBC65RN1tKt3p9UCz015KuQvAYK22rxZFRepOSG7dSv7ykSPV26Y94uP1u9ozjJqUlGi/j2nTvIsWs0dCgrrbc4cWnTlpsajbNLS2Fli4UPursBDaW/MM4wuqq4HKSm33sX07sGGDe2VXm8vYDA1l4daNkyc9T2u3x/r11LDgwgvV26Y9kpM5ioQJDnzhJnn7bffP8+YyNpOS9A2/bbHCXV0N5Oaqu83MTODee7X9QUNDgZQU7bbPML5CSiAvT9t9nDwJ7NgBjBrl3vucZWwKoe4kpye02KmtY8fIVaIWUgILFqi3PUd07OjbThsMoxWnT1OJCS3JygKuusqzJBlHXXHatNG/iFuLtLirqz0vImWP8nLgppuAigr1tmmPyEhyxTBMMKC1tQ0AY8cCM2Z49l5HGZvJySoO0EOaFW4hxL+FEEFVuigvjyxkV2luouLjjykJRuvU13btOK2dCQ4qK8ni1pIvvwTWrPHsnHGUsZm024SYGPXH6i6uWNztAGQKIVYLIUYJEdjSYbEABQXuvcfZREVtLfDuu8Ctt2owWCv0nsVmGDXR2tqurQVefhno1Mmz9zsqRZF42PVSFFrSrI9bSvmoEGIOgCsA3A7gRSHEagBvSCkPaj1AtSksdD+e01m39bJi4JprgH79NBluPYmJ+gX7M4ya1NbSeaglmzaRoTNsmGfvt5exWXWBEVEDPKtMqDYuSYGkgiYn65ZaAG0ArBFCuNaGwo/w9Epvr7SklHQbNmWKumO0h96z2AyjFidOqBsYYI/cXGDqVHVdi/7kqnTFx32fEGIbgAwAPwAYIKWcBuBcAOM1Hp+qnD7tebC/vYmKXbuoAqDWtGlDdUkYJtCprNTeTVJdDfzzn0Ca3fJMnhER4R+TkgquhAMmALheSvmX9R+llBYhxFhthqUNnh4w1hMVZWlGlKYZ0X1WOh4bsAOXXdZR3UHawZ8OGIbxBrXrAtlSXQ38/e+Uwdy1q3rb7djRf6xtwAWLW0r5mK1oW/1vr/pD0obKSs9rItibqNg7fw02ZSZg9OiG9bxtbGqPsDDfNyJlGC0oKdEukqTd8gx0fuoufPifg+jUiUQ7NsuEzk/d5dX5B1AxKX9rCdhipru86dFgr9t6Yb+LMXV6VKM2SGo3NgX0K9TOMGoiJZCTo932y/sNRfnGb/HW2njMueInxGaZ0OOBcWiz+X2vzj/AP+sC+ayRgito1UjBWaF2Tzl1yv5VWBFr2+gTT9GyUDvD+Iq8PO3dJPjuW5yctRijxCZASsiwcBxcuNar8y8xkcon64GzRgotwuJ2VqjdE8rKgHHj7E90qtXYFCCfWlyc5+NkGH+gslLdKpz22L0b2Jc8EoNv7Y+QygqEVJ1B3oR7vTr/QkLIt+2PtAjhVrveb1YWxW1HRjb9n1qNTQESba5LwgQyUlJbQC3D/6qrgXnzgIJvdiP5vf/CHBkNc0QUkle94NX51769/zYqYeH2gK1bgQsuaPp3NRubAuzfZgKf48e1r+Hz+utA19h8THv7QkAIHFi8AQf++ykgJXrMvM6j8y8iguK2/RU/vZ6oR3U1cOaMettrtzwDg6OvRp+L+9T/TWleCsBhx3ZPbtm0rn3CMFpSWkplVbXkzBlqyv32Be+juMcEFF05of5cO7hoHRI2rfLo/EtN9a/wP1uCfnIyPx84ckS97UX+/DXOeuRGHH6WBNo2xltNBg7031s1hnFGbS3w++/UXV0rqqrIlRgSoq7IJiRQY2G90aXLu7+gtpvkveN/Q/bg7Vg2a7BqkSP2CAlh0WYClyNHtBVtAFi8mCK71Cw5EREBdO6s3va0IqilwWKh2zVPaLc8A+X9htYLcrvlGbCEhuLol39Du8uHIL8nRY6cHnqp6qIN2J/4ZJhAIDdX3Xr39vjpJ+Dbb4FVq9TbphCUuBMIAQFBPTlZWur5bLZtMo0lNBSdFs/E74XtcI7cXj97Hf17plcz145g4WYCkeJibRNtADqvH38ceOIJdcNl27cPnCzloBduT7HtOZeybAGOTl+IfnkmXPXS2PrZ64OL1nkVOeIIbgbMBBpnzlDon9bExQGvvqpuEanY2MDq5RrUwu1pJUAF22Sa/Ftm4IXBbyGp+nh9cL915IiaBMLtGsMoVFcDf/6pfbnWhQspHFfNAlIhITQZ6c9RJLawcDvBNpnmz2c+wLO7RjVJrilLM9otvM4wLQGzGThwQPvJyA0byLfdv7+62+3cOfDKJgft5KSUZAV4im2YX21cPE4/vxkH+t2LE1P715d21SKihGECBYsFOHjQ81wJ2yAAoCEvwtoY2rePokhefVVdP3RCQmC2BAxai7uqyr2GwLbYlnI11Nbiz+GTkBRJjnOtXCQMEyhYLMChQ97NJblaUbNdO2DBAiq6phbh4YER+mePoE3AKS4mS0BNXn4ZaNsWuOmm5td11ZJwRKdO3ECB8V8U0VYjT8JZRc3SUiAjA3jkESAqyvt9KRgMQO/e/p2d3CKrA3rr37bH3Xe7JtqANrW5GcYfUNwjaiW3OaqoWVMDPPgg1exROzy2a1f/Fu3mYOG2orkONhs3ul6e0jac0F1/uDf+eYbRCouFJiLV7GTjqKLmk0+SP/v++9WN+OjQwf862rhL0Ap3VZX772nOSv7kE/eKwXtTm1uLOwaG8QYlesQbn7YtzipqXnEFibeaobEJCYEVr+0IFm4rytKMKDZejx4PjGtkJQNkjRsM7jVk8KY2Nws3409oIdqA/X6uD478GZvXluHCC9V1kcTE6NfNRm2CVrg9pejKCRDm2norGUC91W0wuB6p4m1t7qoq7ZMZGMYVqqooHK+sTP1t2/ZzXb4cWLezO866/2pV9xMRAfTsSZOSwUDQxnF7gwwJhTkyGu1WLkTyqhfq+9ZNb+t6zKc9S8Ld2txVVerOpDOMu5SXk6VdW6v9vjZuBNauBV57jXo9qkVYGHDWWcFVbTOIPor3KFbywUXrEJdlQoel82HtGYmNdd0Kthfyp6TIu0pFBQs3ox9FRUB2tnf5EK5SUwOMGAEMGtQQButtSC1AFnbPnsFX+ydIbhya4snBpljJAOp90zIsHAmbqHbkihXAp5+qOUrnqF1LnGFc5fhxKhjlC9FeswaYOZMMo/btG/7ubUitEED37oEd9ueIoBVuT1Cu4ta+6YML1yLe9BGFBbajWsO+4vRp9nMzvkVJrDlxQvt9SQksXQqsXAk82+NVu/M/xcbrPQqpVWprB2vfVhZuG5z5pn0t3Gaz+rP4DOOIqipg/37tmyAo/PknsGUL8MYbQOLwXnat66IrJ7gdUisEVfsLxBokrtKihLu5BBug6Sw30FD9r08fYORIdffXHIWFru+PYTzl1Clg716ajNSa2loqzdqrF7kfExMdJ6wBcCukVhHtQE+waQ7NhFsI0UkIYRJC/C6E2COEuE+rfbmKtz6z1FRgzBjX/X5qpL0XF7sXO84w7mCxAH/9Re4RXxxnVVXArFnAu+/Svq0jPWwT1gC4FVKr+LSDXbQBbS3uWgAPSCn7AjgfwL+EEH013F8j7MVrepuGDgB33OF68So19iclkJfn8uoM4zJnzpCVXVDgm/0VFwP/+heF5y1a1PQctU1YS9i0yqHb0hYhqHJgfLz2n8Mf0Ey4pZQnpJS/1j0vBbAXQEet9meLo4wrb9LQAara99tvrq/v7f4A4ORJ38TRMi2HvDwSbV9l6EpJIXmXXQY89RSJtzX2EtbiTR812Y69piVKyF+wTkTawyc+biFEVwCDAfzsi/0BjoXbmzR0ABgwANi1y/X1vd0fQLeUvpjlZ4Kfmhq6Yzx61DehfgB1rbnrLhLuCRPs3w07CwpwhsFAlnarVlqM3H/RPAFHCBEL4EMA06WUTWqKCSGmAJgCAJ1VrGpuT7htu9p40sVm0CDXRVSN/Snk51Mx+UBrscT4DwUF1IHdl3Mm778PvPkm8OyzztPNPUlYUyxtNTu9BwqaNlIQQoQB2ABgk5TyuebWV7ORQmkphTZZo0Ymljuovb+2bdVtksq0DCoraQJSi1ojzsjOBmbPBv7zH5rYVxODgdLY1Wxj5m84a6SgmXALIQSA5QCKpJTTXXmPmsJdU+OeS8Md1q+ndODbbtNm+85oab48xnMsFso7OHHCd24RgOZktm4Fxo0j617NsqwA+cd79gzOjEhr9OqAcyGAWwFcIoTYUbdcpeH+GhEWpv4Bo9CtG6W+69H1LTubmywwzVNWRpOPx4/79jj96SfgH/9oaLTgyjnoTr5DdDTQp0/wi3ZzaBlV8r2UUkgpz5FSDqpbPtNqf/bQqkBT//4USqV2T0tXqK2lmFs/ahXK+BE1NeQW2bfP9zXdv/kGmDcPePppYNIk19/nar5DfDz1ibSNSGmJBHXmpNp96hSEAG68kW4J9aC83PUWakwLICMD5i9NOHEC2L2bJiHdzdD1hlOn6E7wvPOo7si55zpf39bCjtmTiaLLbkTP6VfX5zucuH12o4iS9u0peiRY6ml7S1B/DVpOXPzjH1SGUq8iULm5vqspwfgvUgIlvYZCpqfj9HoTLBbfNqbeuhW45Rbgxx/JUHKljrathW0JDUXyh6+gbOBwdFg6HyXDRyNl2QKU9xtaXyyqo88yQAKDoK7H3bo1WcdauRVWrqTb0Tvv1Gb7zXH4MH2+lpItxjSmpITC+yo7GRG7gDJ082+YhqQ1SzwKOXWXpUup8cHjjwPDhrn+PuuMYmW8eeOnIvnDV1A6+CK0/fxtHJ2+EJXnG9GrR3BHjnhKUFvcoaHUZ04rRowAPvhAv/6QUpKfnS3vlkV5OYW6HjjQcOypkaHrKkpHnJEjgffec0+0FazHWzJ8NBK+/ACFoycibvt3KBw9ER3eWoA+J00s2g4IauEGtLVGu3WjTMoPPtBuH65w6BCFJzLBTWkplUL944+m5X7VyNBVcBTlkbjsP3jjDWDaNLrb69XL84xF6/G22fIhii67Ea1//BzHJ89B/E+fQzwyG+E7nWdNtmSC2lUCUKWwnBzttn/vvdrFi7uD0q2kbVu9R8KoTUkJTYQ7SqBRM0MXaPBBK++PzTIh5eFJ+Fvb3QhJIBehdacad7Edb21cPDotnomj0xfCMGMGQm42QtyUDqxe7flOgpygF+7wcPKRaZU11qUL0LkzzeZfus23mZm2ZGdTmGKHDjz7HuhISS6wkyfpN3WGN42p7WX3Ag2dZ3KuuwdHPvgFWLgCd9S2xrBh3h9btuM11Nbi+IyFSGpVi6hUAKlGEu3MTMCorZ8+UNE05d1d1MyctKaggGJbtaKyEhg/Hnju1l9x49IrG1kq3lg+nhIdTXWJg61BakvAYqHmGbm5VLtaa2yPUevXezccxOwNF6J39xo8+v45EEKbMcTFkduR47Mbo1fmpN/Qpo16Fqg9/1/ibhPmDFqPuWuH4Pd5a7yqv60GFRWUNceTloFDZSW59HbtAo4c8Y1oA45rxn/4Y3v869MxmHHpTqwuvBRx2zz3mTtCCArz69WLRdtdWoRwh4So13/OUZbXpePi0L07sCL74mZn99VoadYcZjNNWh45wh10/BWLhSaV9+8H9uwhK1uP30qJ8EhZOh+vDXoZR/Ijcfu6cVi36ACGPDsBh5913nnGE8LDKQvSG195S6ZFCDegnt/XkYVSPtSIxx4D7uhpQsIHrzqd3fe0pZkngp+fT6JQUMBp8v5CVRVZ17/9RpPKvmwIbe8YSnrnOeSs2ILz2mVj6Xe9EWbajNMZr8Aw8iIArtfGdpU2bYC+fbUN1Q12Woxwh4UBycnqbMtRzGzibhNSZ0/COTEHkTnWcY88T1uaeSr4Sv2KvXu5a7xemM0N1vXu3WRd69HVyPYYSnrnObR9fg7+HvsxrryzC975XyEu/PV/Td5nr/OMuxgMNJnfvbt2BeBaCi1iclLBbCYrx9vbUUUwbbPUlBn6pQeM+Phj4PXXgfZ/OI4qSXllLjosnY/jk+fgxNR5Xu3bHeLjybeoVS0XhrBYqEpeURGF9OlVHsGW2CwTUh++BRlnv4kDmcVYcN9x5N48o/6OVItIqKgoEmw+5lzH2eRk0IcDWhMSQj41Tws0tVueAUtoKFKWLWgUg9pz+lgcWLyh/kC/6VyycB97DFi0yH4XD9uEidJmun0oWFv7xyfP8Wjis7iYljZt6C6Es9PUQxHrU6dIrP1xfuHzM0a8IHfh7J9/wawbcpF/y4xGt97NdZ5xByGoc1NKCoeoqkmL+yqTkz2fwS7vNxQdl8zFidtnNyQmLFuAY9PmN/L/CQHMnAncfz/5lW1PXnuNUV2d/FEzQ+7UKSr/uXcvhaD5i0UYaNTUkFV96BCwcyeVISgq8ly0tZi8lhLYsYMeQ/74DS+ap2Hp5K244MsnVZ10tCYmhmpnd+zIoq02LcriBugASkmhaAt3KUsz4sDiT9B9VjpCS4uduipCQqhd08qVZOE//DDq42A9TZhQO0NOoaKCkndycqi6W2KiD2PAMzKAoUMbJ1qYTJR88ZD2SUueYLFQQtfp07Q0lyDjLvYyF5XXnrBrF/Dii3RxfvvuHzD1/XE4tGg1Tqh4DFljMJBYqzWnxDSlRfm4FaSkSAtPY2Xd8U2XlQF3302hT7NmeTcp48uemdHR5AuPj9euIQUAEun0uvRmo7Hpaz+hooJEurSUFq1PGzXmMgAgK4tcdlOmAGPGAB3f0fYYat2aMom5qbX36NJz0hN8JdwAVVjbt8/9E9CTE6q8HHjwQeDyy4HrrvNi0DoREdEg4jExUD+DThHradOAJUt0F+3aWvrNKirosbxcnwgQxUA4PfRSnLjj/1wW2927qbP65ZcDV15Jrhyt76DCwoBOnWjehFGHFp85aY+YGPeLs3vim263PAPt9pqweDFwzTVUEhPffeuz7iRqUFVF4Wv79pEP98CBhqJHqvjFjUYS7fnz6dEXop2RAZhMMJvpc+TmAsffMeHkjIz6z3j8OE0w6iHa1nMZ0Xuz0OOBcc2GgSp3dw8/TN1ojEZyW2gt2omJQL9+LNq+pMUKN0Cz3e6UfXXmm3aE4q9M2GVCSAjw1bK/8M8HErA3+SIvR68PZjOJ2bFjJOQ7dgC//05RNPn5JB41NW5u1GQiS3vOHHo0uTFZVifATbaX0XBhlJIuPiUlQF4ecPQokJMyFLXj03FgqQn79gHF60xIuicdJb217xrTHLYGwsGFawEh0OOBcU3i/qUEvv8e2LyZjJEbbwTWrQNuukn70LvISHIBdunCcdm+psW6ShTMZoqq0LI2hLV7JfGDJZg36ge8/EUvzJsHXHCBevvxpQ+8OQwG8nNGRDQ8Ks9DQhoW8bXnPm4pActXJhhuTseZt1ajargR4msT4iano+Cl1ShNM6Kqin5be4e5Wn5ktXH0O6a88RRaZX5VP7eyfj3w7rvkurr7buAiH9kCBgOF1bZvr4HbjKmH47idEBJCiQF//KHdhJNt7PV1U3vhrNFksVRV0RhCVfgl1I5G8AaLhQonNdcdqP26TJxZsBqVSUaE7AUMHYyIeW41Ij7LRFEHsijbLM3AmX5DUT6MhExKIHKrCZG76IIU++RqdJ+YjuI6AT7wzGqUnWUESpzvW42YeC1wdJGN+nMntt2cgWPv/YreaSbk5hpx333A+ef7RkCFoHrvHTpwUSi9adGuEoXoaJpY0Qp7sdf9+wM9ewIffghMnapOswdPU+n15OSkh1AymCzjioo6f3NfI45MeAhlZTQxWHzWUHSckY7Q702orATCfjChy4MNPl5P23apERPvi4JhsVkmFM5cgHE9d8O44UF8etEz6D4rHfcPNuGCC3wj2q1bU0x2ly4s2v4AC3cdSUnqVRC0prkJzQkTgIsvpq7x777r/f582XvQVzR3QfJEgL1JgrKmvN9Q9Jw+FknvPNdou5bQUK/Fu7wc+Plnmlt5e+Rr6DWiHdavB6Y82UXVok/OiI6msqs9e2ocFsq4RYv3cVtjsVDW2+nT6m3TVb9zTg6wbRtw7bVUyS8x0f1tKH/3R7+tGtiLn3fWCMDZ51ZzPiDpnefQafFMFI6eiNY/fo4Tt89uVBbBXX7/HVizhtz9550HLFjge19yeDhFXWlhzDCuwXHcbqB0Ti9pxj+qFaWlwA03AFddBdx+OzVjdVWcPBWxQKC5wl56T8h2mTsJiZ+tROngixB5eK/b3/mpU2RdjxoFrFpFcwNjxza+gPuCkBDKLE5O5olHvWHhdhMpqe5EcbE++y8oAF59lSyuhQuBQYNcs6T9RcTUxt8vSMp4Krv1Qdz271Bw1a34a94Kl977888UQLNtG7nM5sxRZ6LaXYQgsU5J4dA+f4GF2wOkpCL3erb/ys6mk+nQIboLuPPYY+j05jy3ysAGA/58QVJEW3GPlAwfjbafv42SC8cgd+KMJmMO3fUrPu7xAH75hQqRff01TcgajfpUaVQiRVJSOE3d32Dh9hApSTyLivQdx/79wPOPncKpg6dw+/l7cO+eqTj27Nt+YW22dOyV+k165zmkvvR/sISF4+CidSg914i4bSZ8fv8XeNwwD917heGKK6j8gR7WNUCx2ElJlITGUSL+CQu3l2RnU2U1X2JtZcZmmdDt4XR8cuX/8MHuvnj0X0WonPU4zjz6FKIvvdC3A2OaYO+OwPzNd/j1vX3YsrMtfom8GH8Y+uDL+9ZDDDtP1z6LoaF0F5eUpN9Fg3ENrlXiJV27kmXiS6xbTMXsycTJf87GNV/8G0/cW4iq8/6GNVcuxdVzz8UTTwDbt7ecfpK+iJt2l9x/PISyNCNychqy7xduvQjvRk3GsPMEtpb2RdGNdyHlGv1EOzycchUGDCC3CIt2YMPC7SKpqUCPHp4f8O4KjnXssuFMWZPwshtn9cSHn0aia1fg/ffJV7l9uzqJPP6Mp303tUApPrVsGTB+PHDHHcAvv9DfHn4YeP0WE+7bfScsk6d43fTCUyIjyfDo358sbW5oEBzwddcN4uMpIeHwYZpQcgdP0tGbS8mOj6fEHYWdO4F33qETdcwYYNw498YYCFhf0PSIVc/KArZuJYGOjQVefpkq411wASWqWPdt1KLphatER1MtEa7YF5ywj9sDpKSSnydPuvc+d5NjPEmmqakBfviBypJOngysWEEdWi66CDj77OCxuDxptOwu+fnUPea330gEJ0wAnnySYquHDSO3g6OJPT0iYYSguP/kZHpkAhuenNSI06fJ+nanXrOrgqNW7PLu3cBXXwHffksXnDVrqKypEOTzDMQkCy2yQ8+coYvdnj10gevbl7JYe/cmgT7/fLKs/ZGwsIaWcxzSFzywcGtITQ2Jd2lp8+u6IzhaWGxKKv2aNeSXra2l5J5HHyXfvS+K7nuLtxc0Kalpwv79wJ9/UoZqWRlw223kYurXD7j6ahJrKf37wta6Nf2erVv79zgZz2Dh9gH5+eQ+cWR9+1v2n5TAiRPkF7/ySuCLL6gBTYcOVFBo8mTqHVhYSBE1/iIM7lzQsrOp2cPRo9T44dFHgY8/pqxUpXDSjTfS55MyMCIt2LpuObBw+wizmcQwL69peJ4/Z/8p1NSQ2P35J5CWRvUy7ryTLNIOHYCJE8l98NFH5ENNSSHR82U9DSmpal5BAV0sExOBbt2Al14igc7NpUnb558HXnuN3B+dO1NU0JgxlM4diH7+Vq0o9pqt65YDC7ePqaoiC0/PdHlPcHRxkdt3YqdxOmJiSKxffJFalR0/ThXkMjLIWt+6lYQlPh5YtIjS9DdsIMswPJws+65dgQ8+IPExm+n1iBFkCR87RhUaQ0Ko9eRXX1G98rIyWp58krbzz3+SYCclUeTM6NH0/qgompjr1InSuAOdsDD6HImJ/u/CYtRHlw44Qog3AYwFkCel7K/VfvyRiAjqqlNWRlZgRYXeI3KNiJyDaL/saRxcuLbendNj5nU4dflN6NmzYb177mn63gcfpAtVSQkV54qIAOLi6HuorqbFYCCxPn6c3hMS0thaV9qaKXWfe/UCJk2isLvY2IZ6Gt9+23T/116r2tegKyEhFMKXkECfma1rxh6aWdxCiJEAygCscFW4g8XitqWwkKxJt5vo+pjYLBN6PDAOEAJ5E+5F8qoXAClxcNE6rouiIQYD3aUkJJBLhMWaAXR0lQghugLY0NKFGyAXQGEh+b+b68OoJ7FZJvS8bwxCqs7AHBGFA//9lEVbAwwGciu1aUOPgeh3Z7SFmwX7AUo1tqQkiv/Oy9OvWUOzKCYfm36qEhJClnV8PFnWLNaMp+h+6AghpgghsoQQWfn5+XoPxye0akWhaAM3ZaDzQVOjwvV6FkxSfNoyJBTHJ8+BDAlFj5nX6VJjw5doWbgqIoIu1j17AgMH0mRsfDyLNuMduh8+UsrXpJRpUsq0pKQkvYfjU0IvGIqkf6XjnEITOncG2u5St2CSu4KUsGlVvU/7xNR5OLhoHSAl/T1A8ESE1SxcpVjVnTtTYaf+/ek5h/ExaqK7cLdojEZg9WoYJqQj6X9z0fWhdNS8vRqRo4yqJIO4K0hhuUdxfMpj9T7tsjQjjk95DGG5R70fjI/wRISb6yLvDCEo+qNDB0qVHziQqkgmJXEIH6MdWoYDvgfgbwAShRA5AB6TUr6h1f4CFqORgpbnzwfmzEH0GCO6AOjShcIJi4tpqapyf9PuVtLLnfQgus9Kx5neg+vDAZVysoGCp9UDm6vEaE1kJLm74uJo4R6NjK/RTLillDdrte2gwmQCliyhLrFLlpCQG0k0lPjl1FQqgqSIuDtx4e4Ikt4lU9XCnc+sEJtlQtKaJTg+eQ6S1ixBaZoRZWlGCEElUqOjgZgYEmpONWf0hqNK9MRkAtLTqc23ItjWr62IiqIlJYWSWYqLKSqlrIxCDR3hSJAc4Yno+RvufmbrujG1I4wwXGLEWf9OR9WK1YgcbWTfNON3sHDrSWZmY5Gu83kjM7OJcFsTHk6p3cnJVLujooLqdyip4UqijyfF/N0VPX/Dnc8cHk5WdFJOJmreXo2eo4zk9uhnBNqvRlRmJnBV4Hx2puXAtUqCkOpqEnLDwgyUnj0UBQOMMJvpf84KW/lbBUNPsK23IgTQZocJcX9k4sy/H0JEBPmoo6MDoxog03LhIlMMqqrIT15ZScJeVdXwqBwCgVDB0BHh4agXZevHiAgOw2MCExZuxilKEShrMVee19To10HeYCCrWFnCwpq+VgSbE1qYYINT3hmnKGVXY2Pt/19KqupnNtNEqPLc9rXyXEEI+4u9/1mLtCLQLMYMYx8WbqZZhGgQVYZh9IdtGoZhmACDhZthGCbAYOFmGIYJMFi4GYZhAgwWboZhmACDhZthGCbAYOFmGIYJMFi4GYZhAgwWboZhmACDhZthGCbAYOFmGIYJMFi4GYZhAgy/KusqhMgH8JePd5sIoMDH+1QTHr++8Pj1JZjH30VKmWTvH34l3HoghMhyVPM2EODx6wuPX19a6vjZVcIwDBNgsHAzDMMEGCzcwGt6D8BLePz6wuPXlxY5/hbv42YYhgk02OJmGIYJMFqscAshRgkh9gkhDgghZuk9HncRQrwphMgTQuzWeyzuIoToJIQwCSF+F0LsEULcp/eY3EEIESmE+EUIsbNu/E/oPSZPEEKECCG2CyE26D0WdxFCZAshfhNC7BBCZOk9HncRQsQLIdYIIf4QQuwVQlzg1vtboqtECBECYD+AywHkAMgEcLOU8nddB+YGQoiRAMoArJBS9td7PO4ghEgBkCKl/FUIEQdgG4BxgfL9CyEEgBgpZZkQIgzA9wDuk1Ju1XlobiGEmAEgDUArKeVYvcfjDkKIbABpUsqAjOEWQiwH8J2UcqkQIhxAtJSy2NX3t1SLexiAA1LKQ1LKagCrAFyr85jcQkr5LYAivcfhCVLKE1LKX+uelwLYC6CjvqNyHUmU1b0Mq1sCygISQqQCGANgqd5jaWkIIVoDGAngDQCQUla7I9pAyxXujgCOWr3OQQAJRzAhhOgKYDCAn3UeilvUuRl2AMgDsFlKGVDjB7AYwEMALDqPw1MkgC+EENuEEFP0HoybdAOQD2BZnatqqRAixp0NtFThZvwAIUQsgA8BTJdSntZ7PO4gpTRLKQcBSAUwTAgRMO4qIcRYAHlSym16j8ULRkgphwAYDeBfda7DQCEUwBAAS6SUgwGUA3Brnq2lCvcxAJ2sXqfW/Y3xEXW+4Q8BvCOl/Ejv8XhK3S2uCcAonYfiDhcCuKbOT7wKwCVCiLf1HZJ7SCmP1T3mAVgLcn8GCjkAcqzu0taAhNxlWqpwZwI4SwjRrW5iYAKA9TqPqcVQN7n3BoC9Usrn9B6PuwghkoQQ8XXPo0CT3H/oOig3kFLOllKmSim7go79LVLKiToPy2WEEDF1k9qoczFcASBgoquklCcBHBVC9K7706UA3JqYD1V9VAGAlLJWCHEPgE0AQgC8KaXco/Ow3EII8R6AvwFIFELkAHhMSvmGvqNymQsB3Argtzo/MQA8IqX8TL8huUUKgOV10UkGAKullAEXUhfAtAOwlq7/CAXwrpRyo75Dcpt/A3inznA8BOB2d97cIsMBGYZhApmW6iphGIYJWFi4GYZhAgwWboZhmACDhZthGCbAYOFmGIYJMFi4GYZhAgwWboZhmACDhZthAAghhgohdtXV2o6pq7MdMPVHmJYFJ+AwTB1CiCcBRAKIAtWSWKDzkBjGLizcDFNHXfpxJoBKAMOllGadh8QwdmFXCcM00BZALIA4kOXNMH4JW9wMU4cQYj2ozGk3UGu1e3QeEsPYpUVWB2QYW4QQkwDUSCnfrav696MQ4hIp5Ra9x8YwtrDFzTAME2Cwj5thGCbAYOFmGIYJMFi4GYZhAgwWboZhmACDhZthGCbAYOFmGIYJMFi4GYZhAgwWboZhmADj/wHmFI67vugq4gAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot of fit and confidence area\n",
    "fig,ax = plt.subplots()\n",
    "ax.plot(x,y,'xr',label='Data')\n",
    "ax.plot(*line([i.n for i in params]),'--b',lw=1,label='Fit')\n",
    "ax.add_patch(artist(c_up,ec='none',fc='b',alpha=0.2,label=r'2$\\sigma$'))\n",
    "ax.add_patch(artist(c_do,ec='none',fc='white'))\n",
    "ax.legend()\n",
    "ax.set_xlabel('x')\n",
    "ax.set_ylabel('y')\n",
    "ax.set_aspect('equal')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.4 64-bit",
   "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.10.4"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "26de051ba29f2982a8de78e945f0abaf191376122a1563185a90213a26c5da77"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}