{ "cells": [ { "cell_type": "markdown", "id": "70913a1d", "metadata": {}, "source": [ "### Fit $y(x)=\\frac{A}{1+Bx}$" ] }, { "cell_type": "code", "execution_count": 1, "id": "08c5b3cf", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import math\n", "from scipy.stats import chi2" ] }, { "cell_type": "code", "execution_count": 2, "id": "5ce9fede", "metadata": {}, "outputs": [], "source": [ "# the chisq\n", "def getChisq(y0,y,e2):\n", " return ((y-y0)*(y-y0)/e2).sum()\n", "\n", "# the function\n", "def f(x,p):\n", " return p[0]/(1+p[1]*x)\n", "\n", "# These are the derivatives\n", "# i = 0 gives df/dp0\n", "# i = 1 gives df/dp1\n", "def dydp(i, p, x):\n", " temp = 1 + p[1]*x\n", " if i == 0:\n", " return 1./temp\n", " elif i == 1:\n", " return -p[0]*x/(temp*temp)\n", " else:\n", " print(\"Illegal call to dydp with i = \", i)\n", " return -1" ] }, { "cell_type": "code", "execution_count": 3, "id": "2de6cd85", "metadata": {}, "outputs": [], "source": [ "# The dataset : x, y, and the errors in y\n", "# Note: y was generated as y=2/(1+x) + random_number_between(-0.5 and 0.5)\n", "x = np.array([0., 0.5, 1., 1.5, 1.8, 2.1])\n", "y = np.array([1.77, 1.59, 0.83, 0.56, 0.34, 1.11])\n", "N = len(x) # number of points\n", "e2 = np.full(N, 1.)/12. # variance of y" ] }, { "cell_type": "code", "execution_count": 4, "id": "81c817e1", "metadata": {}, "outputs": [], "source": [ "# number of iterations\n", "niter = 10 # a lot!!\n", "\n", "# The inverse covariance matrix of the y's\n", "W = np.diag(1./e2)\n", "\n", "# we will fit as y = p[0] / (1 + p[1]*x)\n", "p = np.array([12., 10.]) # initial guess (way off!!!!)\n", "y0 = f(x,p)" ] }, { "cell_type": "code", "execution_count": 5, "id": "d87ccfbc", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " \n", "-- At ---\n", "[[ 1. 0.16666667 0.09090909 0.0625 0.05263158 0.04545455]\n", " [-0. -0.16666667 -0.09917355 -0.0703125 -0.0598338 -0.05206612]]\n", " \n", "-- A ---\n", "[[ 1. -0. ]\n", " [ 0.16666667 -0.16666667]\n", " [ 0.09090909 -0.09917355]\n", " [ 0.0625 -0.0703125 ]\n", " [ 0.05263158 -0.0598338 ]\n", " [ 0.04545455 -0.05206612]]\n", " \n", "-- W ---\n", "[[12. 0. 0. 0. 0. 0.]\n", " [ 0. 12. 0. 0. 0. 0.]\n", " [ 0. 0. 12. 0. 0. 0.]\n", " [ 0. 0. 0. 12. 0. 0.]\n", " [ 0. 0. 0. 0. 12. 0.]\n", " [ 0. 0. 0. 0. 0. 12.]]\n", " \n", "-- dy ---\n", "[[-10.23 ]\n", " [ -0.41 ]\n", " [ -0.26090909]\n", " [ -0.19 ]\n", " [ -0.29157895]\n", " [ 0.56454545]]\n", " \n", "Iteration no. 1\n", "Chisq = 15.302731583955243\n", "A = p[0] = 1.7691492997948757\n", "B = p[1] = 2.175730103746452\n", "Covariance:\n", "[[0.08332243 0.07966512]\n", " [0.07966512 1.78214122]]\n", " \n", "Iteration no. 2\n", "Chisq = 336.0119810810809\n", "A = p[0] = 1.7967300735822131\n", "B = p[1] = -0.26941802164603246\n", "Covariance:\n", "[[0.08207281 0.15432027]\n", " [0.15432027 0.95699781]]\n", " \n", "Iteration no. 3\n", "Chisq = 43.58890636881732\n", "A = p[0] = 1.1148650427694777\n", "B = p[1] = -0.16716553241190188\n", "Covariance:\n", "[[0.02409395 0.00344169]\n", " [0.00344169 0.00062615]]\n", " \n", "Iteration no. 4\n", "Chisq = 16.482413819807512\n", "A = p[0] = 1.3956871593105673\n", "B = p[1] = 0.1308860859432799\n", "Covariance:\n", "[[0.03481601 0.01229534]\n", " [0.01229534 0.00573888]]\n", " \n", "Iteration no. 5\n", "Chisq = 7.66429797056098\n", "A = p[0] = 1.7510115409801505\n", "B = p[1] = 0.584555797763564\n", "Covariance:\n", "[[0.05590782 0.03186749]\n", " [0.03186749 0.02678354]]\n", " \n", "Iteration no. 6\n", "Chisq = 6.410809864219388\n", "A = p[0] = 1.8506032711652856\n", "B = p[1] = 0.8664517292972878\n", "Covariance:\n", "[[0.07133768 0.05704988]\n", " [0.05704988 0.07996813]]\n", " \n", "Iteration no. 7\n", "Chisq = 6.385058573379272\n", "A = p[0] = 1.8472731978151837\n", "B = p[1] = 0.9051073348295428\n", "Covariance:\n", "[[0.0758692 0.07231552]\n", " [0.07231552 0.13516539]]\n", " \n", "Iteration no. 8\n", "Chisq = 6.385028685203338\n", "A = p[0] = 1.8457950049412064\n", "B = p[1] = 0.903169790129984\n", "Covariance:\n", "[[0.07632048 0.07490154]\n", " [0.07490154 0.14640934]]\n", " \n", "Iteration no. 9\n", "Chisq = 6.385028540993686\n", "A = p[0] = 1.845872135980709\n", "B = p[1] = 0.903320916984551\n", "Covariance:\n", "[[0.07629864 0.07483874]\n", " [0.07483874 0.14609163]]\n", " \n", "Iteration no. 10\n", "Chisq = 6.385028540126873\n", "A = p[0] = 1.8458661304818023\n", "B = p[1] = 0.9033091988170872\n", "Covariance:\n", "[[0.07630034 0.07484519]\n", " [0.07484519 0.14612245]]\n" ] } ], "source": [ "for iteration in range(niter):\n", " At = np.array( [dydp(0,p,x), dydp(1,p,x)] ) # 2xN matrix\n", " A = (At.T).copy() # Nx2 matrix\n", " dy = (np.array( [(y-y0),] )).T # Nx1 column vector\n", " # debug: look at these matrices on the 1st iteration\n", " if iteration == 0:\n", " print(\" \")\n", " print(\"-- At ---\")\n", " print(At)\n", " print(\" \")\n", " print(\"-- A ---\")\n", " print(A)\n", " print(\" \")\n", " print(\"-- W ---\")\n", " print(W)\n", " print(\" \")\n", " print(\"-- dy ---\")\n", " print(dy)\n", "\n", " # find the matrix to be inverted, and invert it\n", " temp = np.matmul(At, W) # 2xN * NxN = 2xN\n", " temp2 = np.matmul(temp, A) # 2xN * Nx2 = 2x2\n", " temp3 = np.linalg.inv(temp2) # 2x2 ... this is the covariance matrix\n", "\n", " # multiply again\n", " temp4 = np.matmul(temp3, At) # 2x2 * 2xN = 2xN\n", " temp5 = np.matmul(temp4, W) # 2xN * NxN = 2xN\n", " dpar = np.matmul(temp5, dy) # 2xN * Nx1 = 2x1 column vector\n", "\n", " # the new values of the parameters\n", " p[0] = p[0] + dpar[0][0]\n", " p[1] = p[1] + dpar[1][0]\n", "\n", " # the currently fitted value of y\n", " y0 = f(x,p)\n", "\n", " # the chisq\n", " chisq = getChisq(y, y0, e2)\n", "\n", " # output stuff\n", " print(\" \")\n", " print(\"Iteration no. \", iteration+1)\n", " print(\"Chisq = \", chisq)\n", " print(\"A = p[0] = \", p[0])\n", " print(\"B = p[1] = \", p[1])\n", " print(\"Covariance:\")\n", " print(temp3)\n", "\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "4e334265", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "FIT RESULTS\n", "-----------\n", "A = 1.85 ± 0.28\n", "B = 0.90 ± 0.38\n", "chisq = 6.39 for 4 dof\n", "prob = 0.17\n" ] } ], "source": [ "# Nicely formatted output \n", "print(\"FIT RESULTS\")\n", "print(\"-----------\")\n", "print(\"A = \" + \"%.2f \" %p[0] + u\"\\u00B1\" + \n", " \" %.2f\" %math.sqrt(temp3[0][0]))\n", "print(\"B = \" + \"%.2f \" %p[1] + u\"\\u00B1\" + \n", " \" %.2f\" %math.sqrt(temp3[1][1]))\n", "print(\"chisq = \" + \"%.2f \" %chisq+ \"for \" + \n", " str(len(x)-2) + \" dof\")\n", "print(\"prob = %.2f\" %(1-chi2.cdf(chisq, len(x)-2)))" ] }, { "cell_type": "code", "execution_count": 7, "id": "74962f80", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAApEUlEQVR4nO3deXhV1d328e8ihEAIEAENIENQBAEVhMggVqJVG5GpqFibIiCWxwGtYyvFR7AtTx1eFVsqfVGpoCmi5VUBRaBPSRxBwQZlhkJABEGGACEQCaz3j5VwMpOQc84+w/25rn2Rs/fOPr8s47mz9rCWsdYiIiLRq47XBYiIiLcUBCIiUU5BICIS5RQEIiJRTkEgIhLl6npdQE01b97cJicnB+W9jhw5QsOGDYPyXqFObeGjtvBRW/iEelusXLlyr7X27Iq2hV0QJCcns2LFiqC8V2ZmJqmpqUF5r1CntvBRW/ioLXxCvS2MMdsq26ZTQyIiUU5BICIS5cLu1JCIBJ+1lmPHjpVbX1BQwNGjRz2oKPQEqy2MMcTFxWGM8dsxFQQiUqnNmzczbdo0li1bBlDuw+fo0aM0aNDAi9JCTrDa4uTJk8TExHD55Zdzzz330LZt21ofU0EgIhXauXMnd999N2PGjGHSpEk0atSo3D6HDx+ucH00CmZb5Obm8u6773LnnXfy+uuv07Rp01odT9cIRKRC8+fPJy0tjVtuuUUf9iEmMTGRkSNH0rt3bz744INaH09BICIVWrVqFX369PG6DKlCnz59yM7OrvVxwi4Idu3yugKR6HDs2DHi4+O9LkOqEB8fT0FBQa2PE3ZBsHMnZGV5XYVIdPDnnSnif/767xN2QVCvHtx5J/zwg9eViMiAAQM466yz/PJXaVlTp04lJSWFuLg4Ro0aVeW+OTk5p2pp0aIF48aNo7Cw8NT21NRU6tevT0JCAgkJCXTq1Mnv9e7fv5+f/vSnNGzYkHbt2vH3v/+9wv0KCgoYM2YM7dq1o1GjRnTv3p2FCxfW6Ofxt7ALgrZtYf16eOYZrysRiW45OTl8+umnGGOYN2+e34/fqlUrHnvsMW6//fbT7nv33XdzzjnnsGvXLrKzs8nKyuLFF18stc/UqVPJy8sjLy+PDRs2VLuOSZMmMWnSpNPu99BDD1GvXj12795NRkYGd911F2vWrCm3X2FhIW3atCErK4uDBw/yhz/8geHDh5OTk1Ojn8efwi4ImjSBm2+GP/wB/vMfr6sRiV6zZs3isssuY9SoUcycOdPvxx82bBhDhw6lWbNmp91369atDB8+nPr169OiRQvS0tIq/BAOlCNHjjBv3jx+//vfk5CQwBVXXMHgwYN57bXXyu3bsGFDJk2aRHJyMnXq1GHgwIG0b9+elStXevbzhF0QAEyZArGxcPfdoCmXRbwxa9Yshg8fTnp6OosWLWL37t2V7jtw4EASExMrXAYOHFjrWu6//37eeOMN8vPz+fbbb1m4cCFpaWml9hk/fjzNmzenX79+ZGZm1vo9S9q4cSN169alY8eOp9Z169atWh/eu3fvZuPGjXTt2vXUuur8PP4UlkHQqhVMngyLF8OcOV5XIxJ9Pv74Y7Zt28awYcPo2bMn559/fqXnxAEWLFhAbm5uhcuCBQtqXc+VV17JmjVraNy4Ma1btyYlJYWhQ4ee2v7UU0+xZcsWvv32W8aOHcugQYP4TxWnFEoG15NPPsmTTz5ZZXDl5eWVe9aiSZMmHD58uMq6jx8/Tnp6OiNHjuTCCy+s9s/jb2EZBOB6AykpcP/9sH+/19WIRJeZM2dy3XXXnTpt8/Of/zwgp4eq4+TJk6SlpTFs2DCOHDnC3r17OXDgAL/5zW9O7dO7d28aNWpEXFwcI0eOpF+/frz//vuVHrNkcD366KM8+uijVQZXQkJCuQ/9Q4cOVfkg3smTJxkxYgT16tVj6tSpNfp5/C1sgyAmBl56Cfbtc2EgIsFx9OhR3nzzTbKysujQoQMtWrTg+eefZ9WqVaxatarC77n++utP3bFTdrn++utrVc/+/fvZvn0748aNIy4ujmbNmjF69OgqP+iNMVg/nlfu2LEjhYWFbNq06dS6VatWlTrdU5K1ljFjxrB7927mzp1LbGxsrX6e2grbIADo3h3Gj4fXXoP33vO6GpHo8M477xATE8PatWv55JNPyM7OZt26dfzoRz9i1qxZFX7PwoULT92xU3Ype+tkscLCQo4dO8aJEyc4ceIEx44dq/AWyubNm9O+fXumTZtGYWEhubm5zJw5k0suuQRw4/IsWrTo1PdnZGTw4YcfVvuce3XuGmrYsCGDBg3i8ccf58iRI3zyySe8++67jBgxosL977rrLtatW8f8+fPLDVR3up8nIKy1YbX07NnTllRQYO3FF1vbqpW1Bw5Yv1q6dKl/DxjG1BY+0dIWo0ePttnZ2eXW/+QnP7EPPvigtdbaQ4cOnVo/Z84cm5SUZI8fP+6X9584caIFSi0TJ048tT0tLc1OnjzZWmvtv//9b9u/f3+bmJhomzVrZm+++Wb73XffWWut3bNnj01JSbEJCQm2SZMmtnfv3nbx4sVVvndaWppt2LBhhUtaWlqF35OTk2OHDBli4+PjbZs2bWxGRka5Y06ePNnm5ORYwMbFxZU67uuvv35q36p+npI++ugje99991WrPYEVtpLPVc8/2Gu6lA0Ca61dscLamBhrR4+uVntUW7T8D18dagufaGmLyoKgpJJBEO28aAt/BUFYnxoq1rMn/PrX8Le/gR8G4hMRiSoREQQAjz8OnTvDL38JBw54XY2ISPiImCCoXx9mznSjk+pBMxGR6ouYIAC47DKYNAneeAOqeLZFRERKiKggAHc7ab9+rldQYgwnERGpRMQFQUyMe67AWrjtNjhxwuuKRCJLcnIyDRo0ICEhgbZt23LDDTfwzTff+PU9qjukc7F169Zx9dVX06RJEzp06MDbb79dq+MFuuZf/OIXtGzZksaNG9OxY0defvnlUts1DLUftG8PU6fCRx/B0097XY1I5Jk/fz55eXls3LiRpKQk7r33Xr8e/5577qnWkM7gHjwbMmQIAwcOZP/+/UyfPp1f/OIXbNy48YyOV1J1h6CG6g9DDW4AvJycHA4dOsS8efN47LHHSo0+qmGo/WTECBg+3N1N9MUXZ3CA1FS6a+wKkSrVr1+fm266ibVr1/rtmEeOHGHu3LnVGtIZYP369ezcuZMHHniAmJgYrr76avr163dq/5oe70xrru4w1ABdu3YlLi4OcMNdGGNKDYKnYaj9xBiYNs2NVHrLLZCb63VFIpEnPz+fOXPmVDnJfU2HoK7NkM7FrLWsXr3ab8c7nTN5j7vvvpv4+HguvPBCWrZsyYABA05ti5hhqI0xbYwxS40xa40xa4wxv6pgH2OM+ZMxZrMx5itjTA9/1tC0qRum+ptv4PbbdUupiL8MHTqUxMREWrduzZIlS3jkkUcq3bemQ1Dn5eXRuHHjUuuqGtK5U6dOnHPOOTzzzDMcP36cxYsXk5WVRX5+/hkdr6ZDUBe/R02HoX7xxRc5fPgwH330EcOGDTvVQ4DIGoa6EHjIWtsF6APcY4zpUmaf64ELipaxwDR/F9GnDzz1FLz9Nrzwgr+PLhKd3nnnHXJzc/n++++ZOnUq/fv357vvvvPLsRMSEjh06FCpdVUN6RwbG8s777zDe++9R4sWLXj22WcZPnw4rVu3PqPj1XQI6uL3qOkw1AAxMTFcccUV7Nixg2nT3MdfRA1Dba3dZa39sujrw8A64Nwyuw0BZhUNhbEMSDTGtPR3LQ88AIMHwyOPwPLl/j66SPSKiYlh2LBhxMTE8PHHH1e4T02HoK7pkM4Al1xyCVlZWezbt49FixaxZcsWevXqdcbHq6navkdhYeGpawReDENdN2BHLsEYkwxcCpT9GD4XKHnf2Y6idbvKfP9YXI+BpKSkM5pm7o476rJ8eQpDhsD06Sto3LjqW7G65+Zy4sQJv09pF67y8vLUFkWipS327t3LkSNHyv2la60lPz+fw4cPU1hYyBtvvMGBAwdo27ZthadC3nzzzSrfp6LvGTRoEOPHj2fq1Kl8/fXXvPvuuyxZsqTSUy2rV6+mQ4cOnDx5kpdffpmdO3dy4403ntq/pscr9tBDD1VaY1kDBw6s1nt8//33ZGVlkZaWRoMGDVi6dCmzZ89mxowZHD58mLi4OJKTk5kyZQr33XcfeXl5vPLKK3Tu3LncsfLz8/nuu+9q//tY2Wh0/lqABGAlMKyCbQuAK0q8/l8gparjVTT6aHUtX25tbKy1N9xg7YkTp9m5f397oFu3M36vSBMtI25WR7S0RWWjj7Zr187Wr1/fNmzY0CYkJNiuXbuWGkLZH/bt21etIZ2LPfzwwzYxMfHUMNGbNm2q0fHKHrumQ1BbW/Uw1CXr3bNnj73yyittkyZNbKNGjexFF11kp0+fXupYETUMNRALLAIerGT7/wVuLfF6A9CyqmPWJgistXbqVPdT//d/n2ZHBUEp0fLhVx3R0hYahrpmNAx1BYwxBngFWGetfa6S3eYBtxXdPdQHOGit3VXJvn5x990wejT8/vfuArKISLQL5DWCfsAI4GtjTHbRut8CbQGstX8F3gcGAJuBfGB0AOsB3PMFL74Ia9a4ISiWL4cuZe9lEhHq1avHsWPHvC5DqnDs2LFS8x2fqYAFgbX2Y8CcZh8L3BOoGipTvz7MnesmtBk6FD7/HBITg12FSGjr3LkzK1eupHfv3l6XIpX48ssv6eKHv2Qj9sni02nd2oXB1q2Qnq7B6UTKGjRoEG+//TaZmZmc0P8gIaWwsJAPPviAJUuW+OWJ46DcPhqqrrgC/vxnuOsuePhheP55rysSCR3nnXcef/zjH3nhhRd44oknOOuss3CX/nzy8/OJj4/3qMLQEqy2sNayf/9+2rdvz/PPP0+rVq1qfcyoDgKAO++EdetgyhS44AJ3MVlEnJSUFF577TX2799f7ulcgGXLllU5zlA0CVZbGGNo0qQJiX48nx31QQDw3HOwZQvcdx+cdx4EcGwnkbDUtGlTmjZtWm59Tk4OycnJwS8oBIVzW0TtNYKSYmJg9my4+GI3dPXXX3tdkYhI8CgIiiQkwPz50KgR3HAD7Coo/9dP1NLcDCIRTUFQQuvWsGAB7NsHA1f/kcMndBFMRCKfgqCMSy+Ft96CVXkdGJHzHAUFXlckIhJYCoIKDBgAMzo9RVZeH267Tc8YiEhk011DlbitxWK2H4rjv998iLPPds8bmCqfkxYRCU/qEVRh3Dmv8cgj8Je/uEHqREQikXoEp/HUU7BnD0ycCM2awT1BHxlJRCSwFAQVyciAZctoUlCAaZ/My7+bTG5uOuPGQVwc3HGH1wWKiPiPgqCsjAwYOxYKCtzQqdu2Ufeusbz1Igz5IZ2xY6FePTeEtYhIJNA1grImTID8/NLr8vOJnTiBuXPhxz92E9vMmeNNeSIi/qYgKGv79krXN2gA77zjRi1NT9cMZyISGRQEZbVtW+X6hg3d08e9esEtt8C8eUGsTUQkABQEZU2eDGXHFI+Pd+uLNGoECxdC9+5w443wj38Et0QREX9SEJSVng7Tp0NcHBagXTv3Oj291G5NmsA//wm9e7ueQUaGJ9WKiNSa7hqqSHo6vPQSB3NzSczOrnS3xo3hgw9g8GAYMQKOHYMxY4JXpoiIP6hHUEsJCfDee/CTn7jnC/7yF68rEhGpGQWBHxTfTTRkCIwbB08+CdZ6XZWISPUoCPwkLs4NX33rrTB+PDz4IJw86XVVIiKnp2sEfhQbC6+/DmefDVOmwO7d8Oqr7klkEZFQpSDwszp1XAi0bOl6Bnv3wty57pZTEZFQpFNDAWAMPPoozJgB//oXXH21G8FURCQUKQgCaPRoNwzF6tVw+eWwYYPXFYmIlKcgCLBBg1yv4NAh6NsXMjO9rkhEpDQFQRD07QvLl0OLFnDttfC3v3ldUQ0Uz82wahUkJ+sRapEIpCAIkvbt4dNP4aqr4Pbb3TWEkL+9tIK5GRg7VmEgEmEUBEGUmOieQv6v/3JTYN58Mxw54nVVVahkbgYmTPCmHhEJCAVBkMXGwrRp8Nxz7kJy376wZYvXVVWiirkZRCRyKAg8YAw88IAbynrHDkhJgUWLvK6qAqeZm0FEIoOCwEM/+Ql88QW0bg0DBrjTRSE1RlE15mYQkfCnIPDY+efDZ5+56wWPPgrDh0NentdVFanm3AwiEt40xEQIaNgQZs92p4h+8xtYu9YNYNeli9eVUe25GUQkfAWsR2CMmWGM2WOMWV3J9lRjzEFjTHbR8nigagkHxsDDD7trBd9/D5ddBjNnel2ViESDQJ4aehVIO80+H1lruxctvwtgLWHjmmsgOxt69YJRo9wwFSF9i6mIhL2ABYG19kNgf6COH8latXLzIT/+uOsV9OoFa9Z4XZWIRCqvrxH0NcasAnYCD1trK/y4M8aMBcYCJCUlkRmEAXu65+Zy4sSJoLxXZa66Cho3PovJkzvTs2cM9967mQEDdmFMcOsIhbYIJXl5eWqLImoLn7BuC2ttwBYgGVhdybbGQELR1wOATdU5Zs+ePW1Q9O9vD3TrFpz3Oo2dO6398Y+tBWuHDrV2z54gFxBCbREKli5d6nUJIUNt4RPqbQGssJV8rnp2+6i19pC1Nq/o6/eBWGNMc6/qKSczk+wpU7yuAnCT3CxeDM8+C++/Dxdf7B5GExHxB8+CwBjTwhh3ksMY06uoln1e1RPq6tRx8yB/8YWbCnPAALjnnvJDAYmI1FQgbx+dDXwGdDLG7DDGjDHG3GmMubNol5uA1UXXCP4E/Kyo+yJVuOQSFwYPPggvvgg9erjXIiJnKmAXi621t55m+1RgaqDeP5LVr+9OEw0YACNHQp8+8MgjMGmS2yYiUhMaYiKM/fjHbhrM0aPdOEXdu7s5D0REakJBEOYSE+Hll90TyUePwhVXuNNGunYgItWlIIgQ113negd33gnPP++uJWRleV2ViIQDBUEEadTIXUBeutQNZ52a6qbF3LvX68pEpJzUVLeEAAVBBEpNha+/dsNav/YadOoEr7wSBnMki4gnFAQRKj4e/vhHN4Bdly5wxx1w5ZUuIERESlIQRLiuXd21ghkzYP16uPRS+PWvQ2jyGxHxnIIgCtSp424x3bDBDW39zDPQsaM7baTTRSKiIIgizZq5W00//dTNk3zbbXD55bB8udeViYiXFARRqG9fWLYMXn0Vtm1zTyaPHAk7d3pdmYh4QUEQperUcR/+Gze6u4veeMOdLvqf/3EPpolI9FAQRLlGjdzdRWvXumkyJ0yACy5wF5dPnPC6OhEJBgWBAHD++fDOO5CZ6a4fjBnjnk6ePx/s0tCZm0FE/E9BIKX07w+ffQb/+AcUFsLgwW7dmjWNvS5NRAJEQSDlGAM33ujGLpo2zV1HGDeuB8OGuXUiElkUBFKp2Fg3iN3mzXD77Vv55z/d6aJbbnHXFEQkMigI5LQSEmDEiG3k5MBvf+vmTb7oIvj5z93TyiIS3hQEUm1Nm8If/gBbt8JvfgPvvuuGsBgxAjZt8ro6ETlTCgKpsebN3S2nW7e6SXDmzoXOnV0grFnjdXUiUlMKAjlj55zjxi3asgV+9St4+213ymjoUA1bIRJOFARSay1awLPPuuEqJk6EDz90w1ZcfTUsWeImyRGR0KUgEL9p1gwmTYLt210wbNjgptC87DJ3+khPKouEJgWB+F1Cgrt2sGULvPQSHDwIN93kxjL605/g8GGvKxSRkioNAmPM+8aY5CDWIhEmLs7NjLZ+Pbz1ljuF9KtfQZs28MgjrucgIt6rqkfwN2CxMWaCMSY2WAVJ5ImJcT2CTz5xw1+npcHzz8N557mH05Yt87pCkehWaRBYa98CegCNgRXGmIeNMQ8WL0GrUCJK795uyOstW9zpo0WL3PwIffvC3/8OBQVeVygSfU53jeAH4AgQBzQqs4icsbZt4emn4Ztv3HWDvXshPd2dNho/HnJyvK5QJIAyMlxXOCsLkpPdaw9VdY0gDcgG4oEe1tqJ1tonipdgFSiRrVEjuPded4fR4sXQr58LiPPOg0GD3HAWuttIIkpGBowd6+v+btvmXnsYBlX1CCYAN1trH7XW5gerIIlOderAtde6h9JyctwEOV98ATfc4CbKefpp+P57r6sEUlPpfv/9Xlch4WzCBMgv85Gan+/We6SqawQ/stZqwAAJujZt4Pe/d3cVzZkD7dq5sY3OPddddFYvQcJaZbfLeXgbnZ4jkJBVrx4MHw5Ll7oxjO69151SveEGFw6PPQb/+Y/XVYrUUNu2NVsfBAoCCQtdurinlb/91j2l3K2bG/iuQwe46ip4/fXyvW2RkDR5MsTHl14XH+/We0RBIGGlXj0YNgzee89dY5s82d15NGIEtGwJv/yl6zWcPOl1pSKVSE+H6dPdE5fgurfTp7v1HlEQSNhq3dpNlLNxozt9NHQozJ4Nqanujrzx4zUstoSo9HQ3MmP//u7uCA9DABQEEgHq1HEf/jNnwu7d7sG0iy92Q2RfdBFceik89xzs3Ol1pSKhSUEgEaVhQ7j1Vnfq6Ntv4YUX3NzLDz3k7ka69lqYMQP27/e6UpHQEbAgMMbMMMbsMcasrmS7Mcb8yRiz2RjzlTGmR6BqkeiUlAT33Qeff+4GvpswwQ1tMWaM2zZgAPztb3DggNeVingrkD2CV4G0KrZfD1xQtIwFpgWwFolynTrB734HmzfDihVunKN16+D2210o3HCDO7WUm+t1pSLBF7AgsNZ+CFTVAR8CzLLOMiDRGNMyUPWIABgDPXvCU0+53sHnn8P997uLyqNGuek3Bw2CWbN0+kiiR10P3/tc4JsSr3cUrdtVdkdjzFhcr4GkpCQyMzODUR95eXlBe69QF8ltMWAAXH89rF/fiMzMc8jKOpsFC+pTp46lW7dcrrhiL/367SUpqYDuubmcOHEiYtuipiL596KmatoW3Yu6n9mh0H7W2oAtQDKwupJtC4ArSrz+XyDldMfs2bOnDZalS5cG7b1CXTS1xcmT1i5fbu1vf2ttly7WulmXre3Rw9on2r1iP+p4kz150usqQ0M0/V6cTo3bon9/twQJsMJW8rnq5V1D3wJtSrxuXbROxFPGQK9e7mG1NWvcyKhPPw0NGsCkbaP40ca3OP98d50hKwsKC72uWKR2vAyCecBtRXcP9QEOWmvLnRYS8VrHjm5qzY8/hl19b+SF1k/QpQu8+KJ7fuHss91MazNnwp49XlcrUnMBu0ZgjJkNpALNjTE7gIlALIC19q/A+8AAYDOQD4wOVC0i/pJU7wC3NXub+xZM5PBhN4fCwoVuRNQ333T7XHaZu+4wYACkpLgH3kRCWcCCwFp762m2W+CeQL2/SKA1agQ33ugWayE72wXCe++5W1WfeML1Fq6/3oXCddfBWWd5XbVIeV7eNSQSMYxxQ1lceql7cG3vXtdbeO89WLDA3Y5ap47rIVx7rVv69nWD6Il4TZ1WkQBo3hx+/nM3++CePfDJJ27+hJgYePJJd22haVP3INuUKe6itLt5TiT41CMQCbCYGLj8crc88QQcPOhGS/3nP2HJEnc6CaBVK7jmGtdbuOYaaNHC27oleigIRIKsSRM3ZPbQoe719u0uEJYscaeSZs1y6zt3dj2H1FQ3WnFSkjf1SuRTEIh4rG1bNxDemDFuQp3sbNdbyMpyM69NKxqFS8EggaIgEAkhdepAjx5u+fWv3cNq//43ZGa6pWQwXHhh6WDQqSQ5UwoCkerKyIBly2hSUOCmQJs8OeAzS9Wt655LuOwy91Bb2WDIyIC//tXt26ED9OvnWy68UM8wSPUoCESqIyMDxo6FggIMuAmTx45124I4zWBFwZCd7ULhk0/cNYaZM92+TZu6C9TFwZCS4obJ8JvUVDdwWna2Hw8qXlAQiFTHhAmQn196XX6+W+/hfLN167oP+JQUePhhdwvqpk0uFIqXBQvcvrGxbgjukr2Gc87xrHQJIQoCkerYvr1m6z1ijBsbqWNHGF00aMvevfDpp75g+POf4dln3bbkZOjd2w2y17u3uzbh116DhAUFgUh1tG3rTgdVtD7ENW8Ogwe7BaCgAFaudOHw+efw2WcwZ47bVrcuXHyxC4XigNC1hgAJhXkIiigIRKpj8mR3TaDk6aH4eLc+zMTF+R5wK/bddy4Uli93y9//7rsI3bixuyZR3GtISXEPvxlvypcAUBCIVEfxdYAxY7AFBZh27YJy11CwtGhRutdw8qSbh6FkODzzjG/uhaQk6HH8SbrU/Zp+b7tTSm3bulNTEn4UBCLVlZ4OL73EwdxcEiP8Tpk6ddwDbJ07w8iRbt3Ro+7W1ZUr4csvYeVbZ7N4/2ieHea2N2vmAqFnT9+/7dsrHMKBgkBEqqVBgzKnlLaOYdf+o2x/abkLhqKAePZZOH7c7ZKY6HtA7tJL4ZJLoFMndweThA4FgYicsQZ1Ck5dWC5WUACrV/uC4csv3Z1KBQVue7160KULdOvmgqFbN7c0b+7NzyAKAhHxs7g4d1qoZ0/fuuPH3TWHVavgq6/cv4sX+x5+A2jZsnw4dOyo3kMwKAhEJOBiY+Gii9xS8vr6nj0uGIrD4auv4F//gh9+cNvr1YOuXd33de3qW9q10y2t/qQgEBHPnHOOm3vhmmt8644fh/XrfeGwapULh9de8+3TsKG7kF0cDF26uH/btlVAnAkFgYiElNhY91DbxReX7j3k5sLatW42t+Kl7Omlhg19oVByad1aAVEVBYGIhIXExPIPwgEcOFA+ID74AF591bdPfLy7W6lTJ/ekdPHXHTu68Ih2CgIRCWtnneUbRK+k/ft9wbBhgzvdtHy5G06j5PzQbdr4wqHkv+eeGz3PQCgIRCQiNW0KP/qRW0o6ehQ2b3bBUBwQGza4U0yHD/v2a9jQ13Po1AkuuMDN+XDBBS58IomCQESiSoMGvmsQJVkLu3aVDof1692gfG+8UboX0bSpC4Xi5YIL4NChxlx0kXvCOtx6EgoCERHch3erVm656qrS244eha1bXU9i0yb37+bNbgTX2bOLQ6IH99wDTZr4wqFsWJx9dmiGhIJAROQ0GjRwdyN16VJ+W0GBC4m33/6aBg0uPhUUn38Ob77pBvAr1qiRG3+pfXs47zzf18VLfHzwfqaSFAQiIrUQF+cuLvftu4/U1NLbfvgBcnJ8PYjNm11obNrkbn09erT0/klJlQdFmzZuvohAUBCISGQq/lT2cAKYevV8M8aVZa17snrrVtiyxf1bvHz2metNnDjh2z8mxj0wVzIckpPdU9bt2rlTWjExZ1angkBExAPGuB5AUhL06VN+e2EhfPONLxxKhsWCBbB7d+n969Z1D84VB0PZpSoKAhGREFS3ru8v/4ocOeKmzN62rfzyr3/Bzp2lr09U+V7+K1tEokZGBixbRpOCAnd+IoJmawsXxeMtde5c8fbjx2HHDl84jBpV+bEUBCJSMxkZbv7mggI3b/G2be41KAxCSGxs6R5FVUGgYZhEpGYmTID8/NLr8vPdeglLCgIRqZnt22u2XkKegkBEaqZt25qtl5AX0CAwxqQZYzYYYzYbYx6tYPsoY8z3xpjsouWOQNYjUmuZmWRPmeJ1Fd6aPLn8I7Dx8W69hKWAXSw2xsQAfwGuBXYAXxhj5llr15bZdY61dlyg6hARPyu+IDxmDLagANOune4aCnOBvGuoF7DZWrsFwBjzBjAEKBsEIhJu0tPhpZc4mJtLYna219VILQXy1NC5wDclXu8oWlfWjcaYr4wx/zDGtAlgPSIiUgGvnyOYD8y21hYYY/4LmAlcXXYnY8xYYCxAUlISmUEaOyQvLy9o7xXq1BY+agune24uJ06cCNm26J6bC0C2Pi9OK5BB8C1Q8i/81kXrTrHW7ivx8mXg6YoOZK2dDkwHSElJsallh/gLkMzMTIL1XqFObeGjtiiSmEhubm7otkViIkDQ6gvn34tAnhr6ArjAGNPeGFMP+Bkwr+QOxpiWJV4OBtYFsB4REalAwHoE1tpCY8w4YBEQA8yw1q4xxvwOWGGtnQfcZ4wZDBQC+4FRgapHREQqFtBrBNba94H3y6x7vMTX44HxgaxBRESqpieLRUSinIJARCTKKQhERKKcgkBEJMopCEREopyCQEQkyikIRCTyFM2pTFaWm1M5I8PrikKagkBEIkuJOZUB35zKCoNKKQhEJLJoTuUaUxCISGTRnMo1piAQkciiOZVrTEEgIpFFcyrXmIJARCJLejpMnw5xce51u3buteZUrpTXM5SJiPhf0ZzKAITprGHBpB6BiEiUU49ARM5MZibZmZmkel2H1Jp6BCIiUU5BICIS5RQEIiJRTkEgIhLlFAQiIlFOQSAiEuUUBCIiUU5BICIS5RQEIiJRTkEgIhLlFAQiIlFOQSAiEuUUBCIiUU5BICIS5RQEIiJRTkEgIhLlFAQiIlFOQSAiEuUUBCIiUU5zFotIZMrM9LqCsBHQHoExJs0Ys8EYs9kY82gF2+OMMXOKti83xiQHsh4RESkvYEFgjIkB/gJcD3QBbjXGdCmz2xjggLW2A/A88FSg6hERkYoFskfQC9hsrd1irf0BeAMYUmafIcDMoq//AfzYGGMCWJOIiJQRyCA4F/imxOsdResq3MdaWwgcBJoFsCYRESkjLC4WG2PGAmMBkpKSyAzSRaC8vLygvVeoU1v4qC181BY+4dwWgQyCb4E2JV63LlpX0T47jDF1gSbAvrIHstZOB6YDpKSk2NTU1EDUW05mZibBeq9Qp7bwUVv4qC18wrktAnlq6AvgAmNMe2NMPeBnwLwy+8wDRhZ9fRPwL2utDWBNIiJSRsB6BNbaQmPMOGAREAPMsNauMcb8DlhhrZ0HvAK8ZozZDOzHhYWIiARRQK8RWGvfB94vs+7xEl8fA24OZA0iIlI1E25nYowx3wPbgvR2zYG9QXqvUKe28FFb+KgtfEK9LdpZa8+uaEPYBUEwGWNWWGtTvK4jFKgtfNQWPmoLn3BuCw06JyIS5RQEIiJRTkFQteleFxBC1BY+agsftYVP2LaFrhGIiEQ59QhERKKcgkBEJMopCNAEOiVVoy1GGWO+N8ZkFy13eFFnoBljZhhj9hhjVley3Rhj/lTUTl8ZY3oEu8ZgqUZbpBpjDpb4nXi8ov0igTGmjTFmqTFmrTFmjTHmVxXsE36/G9baqF5ww1/8BzgPqAesArqU2edu4K9FX/8MmON13R62xShgqte1BqEtrgR6AKsr2T4AWAgYoA+w3OuaPWyLVGCB13UGqS1aAj2Kvm4EbKzg/5Gw+91Qj0AT6JRUnbaICtbaD3HjX1VmCDDLOsuARGNMy+BUF1zVaIuoYa3dZa39sujrw8A6ys+zEna/GwoCTaBTUnXaAuDGoi7vP4wxbSrYHg2q21bRoq8xZpUxZqExpqvXxQRD0SniS4HlZTaF3e+GgkBqaj6QbK29BFiCr6ck0etL3Dg23YA/A+94W07gGWMSgLnA/dbaQ17XU1sKgppNoENVE+hEgNO2hbV2n7W2oOjly0DPINUWaqrzexMVrLWHrLV5RV+/D8QaY5p7XFbAGGNicSGQYa39fxXsEna/GwoCTaBT0mnbosy5zsG4c6TRaB5wW9EdIn2Ag9baXV4X5QVjTIvia2bGmF64z5VI/EOJop/zFWCdtfa5SnYLu9+NsJizOJCsJtA5pZptcZ8xZjBQiGuLUZ4VHEDGmNm4u2GaG2N2ABOBWABr7V9x82wMADYD+cBobyoNvGq0xU3AXcaYQuAo8LMI/UMJoB8wAvjaGJNdtO63QFsI398NDTEhIhLldGpIRCTKKQhERKKcgkBEJMopCEREopyCQEQkyikIRGqhaDTKrcaYpkWvzyp6nexxaSLVpiAQqQVr7TfANODJolVPAtOttTmeFSVSQ3qOQKSWioYcWAnMAH4JdLfWHve2KpHqi/oni0Vqy1p73BjzCPABcJ1CQMKNTg2J+Mf1wC7gIq8LEakpBYFILRljugPX4majeiDUJyERKUtBIFILRaNRTsONS78deAb4P95WJVIzCgKR2vklsN1au6To9YtAZ2NMfw9rEqkR3TUkIhLl1CMQEYlyCgIRkSinIBARiXIKAhGRKKcgEBGJcgoCEZEopyAQEYly/x/wlKhk6yyHVQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Now plot it\n", "a = plt.subplot(111)\n", "xmin = x.min()-0.2\n", "xmax = x.max()+0.2\n", "a.errorbar(x, y, yerr=np.sqrt(e2), linestyle='none', color=\"red\", marker='o')\n", "a.set_xlim(xmin, xmax)\n", "xpl = np.linspace(xmin, xmax, 1000)\n", "ypl = f(xpl, p)\n", "a.plot(xpl, ypl, color='blue', linestyle='solid')\n", "a.set_xlabel('X')\n", "a.set_ylabel('Y')\n", "# txtRes = \"A = \" + str(round(p[0],2)) + \"$\\pm$\" \n", "txtRes = \"A = \" + \"%.2f\" %p[0] + \" $\\pm$ \" \n", "txtRes = txtRes + \"%.2f\" %math.sqrt(temp3[0][0]) + \"\\n\"\n", "txtRes = txtRes + \"B = \" + \"%.2f\" %p[1] + \" $\\pm$ \" \n", "txtRes = txtRes + \"%.2f\" %math.sqrt(temp3[1][1])\n", "props = dict(boxstyle='round', facecolor='white', alpha=0.8)\n", "a.text(0.96, 0.96, txtRes,\n", " transform=a.transAxes,\n", " bbox=props, fontsize='large',\n", " horizontalalignment='right',\n", " verticalalignment='top' )\n", "a.grid()" ] }, { "cell_type": "code", "execution_count": 8, "id": "101cbf40", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABC7ElEQVR4nO2dd3hUVfrHPyedJCQhhCT0UKSJgoAKChrRZQFRXEQsLIKyoiIW1LX8QGRXcFFUQFEQsYAiioUqRSwBQVERgwUUUUOXHkihJOH8/ngzTAhJSDIzmUzm/TzPfWbm3su5Z06G+73nbcdYa1EURVH8lwBvd0BRFEXxLioEiqIofo4KgaIoip+jQqAoiuLnqBAoiqL4OUHe7kBZiYuLs0lJSRVyraysLCIiIirkWpUdHQsnOhZOdCycVPax+O677/ZZa2sVdcznhCApKYm1a9dWyLVSUlJITk6ukGtVdnQsnOhYONGxcFLZx8IYs6W4Y2oaUhRF8XNUCBRFUfwcnzMNKYpS8VhrOXr06Gn7jx07xpEjR7zQo8pHRY2FMYbQ0FCMMW5rU4VAUZRi2bx5M1OmTGHNmjUAp918jhw5QrVq1bzRtUpHRY3FiRMnCAwM5KKLLuKuu+6iQYMGLrepQqAoSpHs3LmToUOHMnjwYEaPHk316tVPOycjI6PI/f5IRY5Feno68+fP54477uCtt94iNjbWpfbUR6AoSpEsXLiQ7t27c/311+vNvpIRExPDwIEDufDCC1m6dKnL7akQKIpSJOvXr6djx47e7oZSAh07diQ1NdXldlQIFEUpkqNHjxIeHu7tbiglEB4ezrFjx1xux+eE4IcfZFMUxfO4MzJFcT/u+vv4nBAEB8OmTd7uhaIoAD179qRGjRpueSotzOTJk+nQoQOhoaEMGjSoxHPT0tJO9iUxMZFhw4aRm5t78nhycjJhYWFERkYSGRlJ8+bN3d7fAwcO8I9//IOIiAgaNmzI22+/XeR5x44dY/DgwTRs2JDq1avTtm1blixZUqbv4258Tgjq14fq1cGDY6IoSilIS0vjyy+/xBjDggUL3N5+nTp1GDlyJLfeeusZzx06dCjx8fHs2rWL1NRUVqxYwUsvvXTKOZMnTyYzM5PMzEx+/fXXUvdj9OjRjB49+oznPfDAA4SEhLB7925mzZrFnXfeyc8//3zaebm5udSvX58VK1Zw6NAhxowZQ79+/UhLSyvT93EnPicExsD+/bB4sbd7oij+zcyZMzn//PMZNGgQM2bMcHv7ffr04ZprrqFmzZpnPPfPP/+kX79+hIWFkZiYSPfu3Yu8CXuKrKwsFixYwBNPPEFkZCSdO3fm6quv5s033zzt3IiICEaPHk1SUhIBAQH06tWLRo0a8d1333nt+/icEAC8+ioMGACa0Kgo3mPmzJn069eP/v37s2zZMnbv3l3sub169SImJqbIrVevXi735b777uOdd94hOzubHTt2sGTJErp3737KOY8++ihxcXFcfPHFpKSkuHzNgmzatImgoCCaNWt2cl+bNm1KdfPevXs3mzZt4uyzzz65rzTfx534pBAMGgT/+Q9kZXm7J4rin6xatYotW7bQp08f2rdvT5MmTYq1iQMsWrSI9PT0IrdFixa53J9LLrmEn3/+maioKOrVq0eHDh245pprTh5/6qmn+OOPP9ixYwdDhgzhqquu4vfffy+2vYLCNW7cOMaNG1eicGVmZp6WaxEdHU1GRkaJ/c7JyaF///4MHDiQFi1alPr7uBufFILGjaFdO/jtN7DW271RFP9jxowZdOvW7aTZ5qabbvKIeag0nDhxgu7du9OnTx+ysrLYt28fBw8e5OGHHz55zoUXXkj16tUJDQ1l4MCBXHzxxSwuwb5cULgeeeQRHnnkkRKFKzIy8rSb/uHDh0tMxDtx4gQDBgwgJCSEyZMnl+n7uBufFAKAwECYMAFmzfJ2TxTFvzhy5Ahz5sxhxYoVNG3alMTERCZMmMD69etZv359kf+mR48eJyN2Cm89evRwqT8HDhxg69atDBs2jNDQUGrWrMktt9xS4o3eGIN141Nks2bNyM3N5bfffju5b/369aeYewpirWXw4MHs3r2bDz74gODgYJe+j6v4rBCEhMDatfDZZ97uiaL4F/PmzSMwMJANGzawevVqUlNT2bhxI126dGHmzJlF/pslS5acjNgpvBUOnXSQm5vL0aNHycvLIy8vj6NHjxYZQhkXF0ejRo2YMmUKubm5pKenM2PGDM4991xA6vIsW7bs5L+fNWsWK1euLLXNvTRRQxEREVx11VWMGjWKrKwsVq9ezfz58xkwYECR5995551s3LiRhQsXnlao7kzfxxP4rBAEBMBbb0HfvnD4sLd7oyj+w4wZM7jlllto0KABCQkJJCYmnox1nzVrltvi3ceMGUO1atUYN24cb731FtWqVWPMmDEnj/fo0YMnn3wSgA8//JClS5dSq1YtmjZtSnBwMBMmTADEDj9y5Ehq1apFXFwcL7zwAvPmzTvFsVuY8sxgnnvuOY4cOUJ8fDw33ngjU6ZMOWVG4Ojvli1bePnll0lNTSUxMfFku7MKmDdK+j6ewLhzelQRtGzZwU6fvpawMPl86JDkFPztbxJa6k4q+9JzFYmOhRN/GYtbb72Ve++9lzZt2hR7jlYfdeKNsVi1ahXvvfcekyZNOuO5xpjvrLUdijrmszMCB2lpcOWV8N573u6JoiiKb+LzQnD22WIe0ocSRVGU8uHzQhAUBA88IO8zM73bF0VRFF/E54UAxDewdy/cf7/mFSiKopSVKiEEIMllM2fCunXe7omiKIpv4TEhMMbUN8Z8bozZYIz52RhzbxHnJBtjDhljUvO3UeW9Xq9eEk4aUGWkTVEqJ0lJSVSrVo3IyEgaNGjAlVdeybZt29x6jdKWdHawceNGunbtSnR0NE2bNmXu3LkutefpPv/zn/+kdu3aREVF0axZM6ZPn37K8apUhjoXeMBa2wroCNxljGlVxHlfWGvb5m//Le/FgoIgKQl27YKtW8vbiqIopWHhwoVkZmayadMmEhISuPvuu93a/l133VWqks4giWe9e/emV69eHDhwgGnTpvHPf/6TTQUWLilLewUpbQlqKH0ZapACeGlpaRw+fJgFCxYwcuTIU6qPVpky1NbaXdbadfnvM4CNQF1PXQ/EVzBrFnTqBMePe/JKiqIAhIWF0bdvXzZs2OC2NrOysvjggw9KVdIZ4JdffmHnzp0MHz6cwMBAunbtysUXX3zy/LK2V94+l7YMNcDZZ59NaGgoIOUujDGnFMGr6DLUQR5ruQDGmCTgPODrIg53MsasB3YCD1prT/u2xpghwBCA+PgE0tNTik0ea9WqBsZE8ckn2wgPP+FSvzMzM91ertZX0bFw4i9jsW/fPrKysk4rpmatJTs7m4yMDDIyMnjrrbdo3759sZU2r7vuOtasWVPksY4dO/JeoSSg9evXExQURO3atU+22aJFC1atWlXkNbLyyxBnZGScXLoxNzeX1NRUMjIyytxeQRwrr53pvPJcY/jw4bz99tscOXKENm3a0KVLl5Pn3nHHHbz55pu0b9+e9PR0PvroI0aOHHlaW9nZ2fz1118u/x49LgTGmEjgA+A+a23hYhDrgIbW2kxjTE9gHnBW4TastdOAaSCZxTExyScziwvTvTt06QJBQY3o1Mm1bGN/ySAtDToWTvxlLGbOnElERMRp2bLGGG666SaCgoLIysqiVq1aLFu2rNis2qVLl5bputZaoqKiTmkvPj6eI0eOFHmNdu3aER8fz9SpUxk+fDiff/45q1at4rLLLqN69eplbq9Xr16sWrUKgKNHjwIwZcoUADp37lxk9VFrLdWrVy/1NQCmT5/Oyy+/zFdffUVKSgpxcXEni89169aNmTNnUrduXfLy8hg4cCA33njjaWsUh4eHk5iY6PLv0aOuVWNMMCICs6y1HxY+bq09bK3NzH+/GAg2xsS5et2ICFi1SqqTKorifubNm0d6ejp79+5l8uTJXHrppfz1119uaTsyMpLDhQqIlVTSOTg4mHnz5vHRRx+RmJjIs88+S79+/ahXr1652itrCWrHNcpahhogMDCQzp07s3379pNiU6XKUBuRrleBjdba54o5JzH/PIwxF+T3Z787rr98uQhBvqAriuIBAgMD6dOnD4GBgSefogtT1gJuZS3pDHDuueeyYsUK9u/fz7Jly/jjjz+44IILyt1eWXH1Grm5uSd9BN4oQ+1J09DFwADgR2NMav6+/wMaAFhrpwJ9gTuNMbnAEeAG66YqeI88IlVJd++Ghg3d0aKiKIWx1jJ//nwOHjxIy5YtizynuDLTxREREUGfPn0YNWoU06dPJzU1lfnz5/Pll18W+29++OEHmjVrxokTJ3jppZfYtWsXgwYNKnd7DkobMVSwDPWZrrFnzx4+++wzevXqRbVq1fjkk0+YPXs2s2fPBk4tQ/3ggw+SmZnpu2WorbWrrLXGWntugfDQxdbaqfkigLV2srX2bGttG2ttR2vtmf8ypSQmBurUgV9+ge3b3dWqoigAV111FZGRkdStW5cRI0YwY8YMtz5hv/TSS6Uq6ezgzTffpHbt2sTHx/Ppp5+yfPnyk1E5pWmvIOVdRKekMtQF+2uMYcqUKdSrV48aNWrw4IMPMnHiRK6++uqTbVV0GeoKiRryFkFB8PjjcOyYZBy7u0y1ovgjaWlpJ997qvRybGws8+bNK/Z44VnG+PHjGT9+fLnbK6nt0lLSNQq2WatWLVasWFFiW23btq3QyLQqLQQAvXvDgQOQkQFRUd7ujaIoSuWjygtBt27iK/j1V+jQQWcFilJaQkJCToZPKpWTo0ePnrLecXnxi8o8UVHw7rswcqS3e6IovkPLli1PKXugVD7WrVtHq1ZFVe4pG34hBACbNsGCBXDkiLd7oii+wVVXXcXcuXNJSUkhLy/P291RCpCbm8vSpUtZvnw53bt3d7m9Km8acvDww7K+8fbtcNZpucuKohSmcePG/O9//2PSpEn85z//oUaNGqdltmZnZxMeHu6lHlYuKmosrLUcOHCARo0aMWHCBOrUqeNym34jBGFhEBICP/0EO3fCpZd6u0eKUvnp0KEDb775JgcOHDgtOxdgzZo1dOzY0Qs9q3xU1FgYY4iOjiYmJsZtbfqNEICsVTBpEvzwg8wM9EFGUUpHbGwssbGxp+1PS0sjKSmp4jtUCfHlsfArIQAYOlTWKzh4UIVAURQF/FAImjSBBg3g558lmsgDuTCKoig+hd9EDRUkOBg+/hjOOUdyDBRFUfwZnxOCgwfd007r1lC/PuzY4Z72FEVRfBWfE4Ldu+Grr1xv55xz4KmnRAjyFyFSFEXxS3xOCEJD4YknwB1rYISGwr59MGAAZGa63p6iKIov4nNCULcu5OVJglhOjuvtZWRIxnE5Cw4qiqL4PD4nBCEh8H//J1E/Eye63l67djB3rqxfoPW1FEXxR3xOCECygm+6SQrJLV/uenvx8ZJs9vrrknWsKIriT/ikEADccw+ce674CzZvdk+bw4eLyUlRFMWf8FkhCAqCceMkO/iBByA93bX2YmJksfs+fcRvoCiK4i/4rBCAmHTGj4c9e+DRRyE317X2LrgAYmNh7Vr3zTIURVEqOz4tBCD5ACNGwLffyhO9q0RGwujRcNllcPy4zw+PoijKGakStYZ69ZKFZ95+W9YauOYa19obOFDaM+aEW/qnKIpSmakyj7z33AMXXih+g++/d62tNm1ETI4ckUxmRVGUqkyVEYKgIPjf/6BOHXEep6W51l5wMPzySxTNm8PKlW7poqIoSqWkyggBSFnp55+HwEC4917Yv9+19pKSsjn3XJkZWOuePiqKolQ2qpQQANSrJ07jffvg/vtdyxaOjMxlwgQpaZGW5npUkqIoSmWkygkBSInpsWNhwwaJKMrLK3sbzYYkc+Ej9wFQs6asbHbTTTozUBSl6lElhQAgOVl8BStWwDPPuHYDDwqStQuCgrRKqaIoVQ+PhY8aY+oDM4EEwALTrLWTCp1jgElATyAbGGStXeeuPtxwg0T9vPkmREfDHXeUv63bbxcR+P57STwLC3NXLxVFUbyLJ2cEucAD1tpWQEfgLmNMq0Ln9ADOyt+GAFPc3Yl77oHevWH6dHjrLdfaioyEP/6QGkfffuue/imKongbjwmBtXaX4+neWpsBbATqFjqtNzDTCmuAGGNMbXf2wxgpW3355VK2ev5819pLTJTX339Xf4GiKFWDCsksNsYkAecBXxc6VBfYVuDz9vx9uwr9+yHIjIH4+ATS01Mwpmx9uPtuw8GD5zB2bA1yczfQufPeEs9PyknH2jx270457diECRJBtHy5+A0CqqynxUlmZiYpKSne7kalQMfCiY6FE18eC48LgTEmEvgAuM9ae7g8bVhrpwHTAFq27GBjYpLLZaOfNAnuugvGjz+bxES4+OLizw0JjuF4TjoJCclFHs/LgxkzYP16+PRTWfayKpOSkkJycrK3u1Ep0LFwomPhxJfHwqPPssaYYEQEZllrPyzilB1A/QKf6+Xv8wjVqol5qEkT+Pe/4csvy99WYCAkJIgguGP9ZEVRFG/hMSHIjwh6FdhorX2umNMWADcboSNwyFq7q5hz3UJUFLz0EjRqBA8+6JoYXHmlhKZu3CjRSSe0Rp2iKD6IJ2cEFwMDgK7GmNT8racx5g5jjCOQczHwB7AZeAUY6sH+nCQ62n1iEBYmZqHu3cXspCiK4mt4zEdgrV0FlOjStdZawCu3T4cYDB0qYvDMM3DRReVrq3p1aNlS2szKgogI9/ZVURTFk/hBvEvxFJ4ZlNfhb4xkMffoAV9/raubKYriW/i1EIBTDM46SxauX7y4/G1FRsLOndC2LTz7rNu6qCiK4lH8XgjAKQbt2sGoUTBnTvnbatZMFrVp0ACys93WRUVRFI+hQpBPRISEll56KTz9NIzf9c9yZQ4HBsLw4RAXB2vWwIcfVoEM5ORk2t53n7d7oSiKh1AhKEBoKDz1FPTsCU/s+heP/zW83CGhUVESjXTttTB7tnv7qSiK4k6qxOL17iQoCEaPhnrffshLewfyxygxF4WElL2t7t0l4axWLTh8WMRBURSlsqEzgiIICIDx9Z7nsYRJLF0Kw4bJjbysGAO9ekmuwZIlcMstEl6qKIpSmVAhKAZj4N741xkzBn78EQYPloig8hAeLtVK333XteQ1RVEUT6BCUAQ1lswi4sc1xP64ngcnJ7FkwCz27ZMn+o0by9fm3/8OH3wAx4/D9u06M1AUpfKgQlCIGktm0XDsEAJyjmGA0L+2kPz2EJbfMouQELjtNvjss/K1HR8v0UTTp0NSEvzwgzt7riiKUj5UCApR98URBB49NQEg8Gg2beeM4PXXJfHsoYfg5ZfLV2QuKAhat4ZWrSTPIDfXTR1XFEUpJz4nBLm5nq3yGbJ7a7H74+Jg6lS46ip45RXJRC5P0lijRpJ5fOiQ+AxefFErlyqK4j18Tgi2boXHHvNcktbxhAYl7g8NlXDS4cNhxQq49VbYUY4VFIwRM9HixRKVtHSpK71WFEUpPz4nBLVrw3XXwd69siDM9Olw4ID72t9x11jywsJP2ZcXFs6Ou8ae/GwM9O8Pzz8v6xDcfDN88035rte3L0yeLG3u3AlHjrjSe0VRlLLjc0IQFyd1/7t0kSfxqVNh3ToxGbljlnCwR3+2jJjGieBQLHAssSFbRkzjYI/+p53bsaMsVxkbK0/106eX3cRjjLRTowYsXAj168uroihKReFzQuAgMhKGDIGff5YZwuHD8mR9//2Qk+Na2wd79CfrnI4cOKcNPy1KK1IEHDRoIGLQrZuI0n33QXp62a8ZHAwNG0Lz5iJoOjNQFKWi8FkhcNCqFTRuDJddJiGZISFw8KDcSHfsqJiCb+Hh8MQT8Mgj8O23Yjb68ceyt1OrlpibIiJg1SpxRv/yi/v7qyiKUhCfFwIHjhpBy5ZJOelDh6BfPykvXREYI/b+116TCqS33QZvv10+IYqJkWikqVNhyhSpV6QoiuIpqowQOAgIgIQE6NpVntK7dROH7u+/w8cfe/6m2rIlvPWWLHv53HNw772wf3/Z26lXD955B/72N3FEr1vnXqe4oiiKgyonBA6qVYN//xuuvx46dYJPPpGwz7Q0z5uLoqIkT+Chh+C77+DGG8XUU1bi4yVK6sgR6NNHqpn6/NoGiqJUOqqsEDgwRiJypk4VMWjSBPbskRv1kiWevW6/fjBzpkQV3XcfjB8Px46Vva3oaPjPf6S9778Xs9GhQ27vsqIofkqVFwIHAQGy+lj79nDBBVLnJzVVInys9ZzJqEkTiSq68UapPnrzzfDbb2Vvp3VrcYgfOAAPPAAtWkjxOo8zaxasWUP0+vXijZ81qwIuqihKReI3QlCQuDiJ6pk6Vd6vXg29e3suQic0VG7ezz8vwjNggOQclKfOUI0acP75knuwe7eYjTxWnmLWLInRPSYF+NiyRT6rGChKlcIvhQCcJR7atJGtUSNZQCY9XTKWj59w/+JtF10ks4KuXUWEBg4s3+ygTRsYOVLCZD/6SArhLVrk9u7CiBGnF1PKzpb9iqJUGfxWCAqSnCzF3/7+d4nlf+QR6LFpokeuFRMDTz4JTz8tZTJcnR0EBkL16uIzyMhwszN5a9EF+IrdryiKT6JCUIDISDjnHBgzBu6t/S7WylP34sWyoIw76doV5sw5dXbw669lb6dBA3j1VahbVyKT/vUvePBBNwlCg6IL8BW7X1EUn0SFoBDGiL+gf70vCAyUcNNRo6T+j7tDNwvPDm6+GSZOLF9p6+rVoWZNiYjatEn8B9a62OexYyVtuiDh4bJfUZQqg8eEwBjzmjFmjzHmp2KOJxtjDhljUvO3UZ7qS3kxyFrFixc7K54uWiSJXu5cUKZrV3jvPbj6aklG69cPVq4sezuBgZJd/dBDkoA2Z45kWf9U5F+gFPTvD9OmQagU4KNhQ/ncv/jaS4qi+B6enBG8AXQ/wzlfWGvb5m//9WBfyo0x0KMHnHeeVDz9/nv48EOnPd5ds4ToaPHBTp8uD9333y8Jcbt3l72tatUku3rPHvEd7NvnQnRR//7QsSOH2rSR6ZGKgKJUOTwmBNbalUCVKooQGQnz5kFKimT97tgh98VPP3XfNdq2lejMYcPEgX3ddTJLKE9FVUeUUk6O9Ll7dxEYRVGUgrg/RrJsdDLGrAd2Ag9aa38u6iRjzBBgCEBCQgIpKSke71jb9HTy8vJKvNaRI6GEhbXg+PEt7NyZTk6O6GpYmOuB/T16QNu2YUydehYTJ9bk/fezue22zXToUD5tzcszhIU14dCho3z88XZCQuDw4SCios5s4yrNWPgTmZmZOhb56Fg48eWxMNaDxWuMMUnAImtt6yKORQEnrLWZxpiewCRr7VlnarNDhw527dq17u9sYZKTSU9PJyY19YynZmZKUbvnnhOn8pw5EtrpLlatkra3bhXz1PDh5Q/cyc2VSKiffhIn+LJl0maJlGEs/IGUlBSSk5O93Y1KgY6Fk8o+FsaY76y1HYo65rWoIWvtYWttZv77xUCwMSbOW/05jZQUUidOLNWpkZGS5DVggJhycnIkMW3tWvcsMNO5s5h47rlHnMD9+sGkSSJAZSUoSHIlGjYUJ3V2NuzaJeU29u51va+KovgeXhMCY0yiMcbkv78gvy/lKNhcebjiCnj5Zblxh4ZKoblx49zjUA4OlvDSDz4Qs9Gbb8K118L775cvgqlxYylkV726iMD114sw6NoHiuJ/eDJ8dDbwFdDcGLPdGDPYGHOHMeaO/FP6Aj/l+wieB26wnrRTVSDVq4sYLFwId9whT9q//gqvv16+HIGCxMXB44/DG2+IeWjcOJkhfPpp+QQnNFQc36NHy4xm5UoJDnroIfjzT9f6qiiKb1Css9gYsxgYaq1NK0/D1tobz3B8MjC5PG37CldcIa8ZGZKD9fLLcPHFUhtI5kLlp3VrCen/4gtZq/nhh2XfPfdI7kBZadZMtuPHJTJq4kTp58CBEGRd7KyiKJWakmYErwMfG2NGGGOCK6pDVZHq1eXJff16WWN5715ZPc3VIp7GwCWXyJKYI0dKzsGQIeJM3ry5fG2GhIjzeN48KVuRkgLP/Hkt3X6bQUaGa/1VFKVyUqwQWGvfA9oBUcBaY8yDxpj7HVuF9bAK0aqVOJU7dRJn8vbtzvUQjh4tf7tBQXDNNTB3Ltx1lziUb7wRHn20/OadhATZYmIgmnSSgrfx++/S35UrdWEcRalKnMlHcBzIAkKB6oU2pZzExMiT9vTpYp9fvVocwK5GZ4aFwS23wPz5MGiQhJ326yezhbS08rUZGAg31fqEKQ1GcuCA9PvKK6X0hjqWFcUFkpNlqwSU5CPoDjwHLADaWWtddHMqhYmKkmqnx4+LczkmRlYgy8iQ5S2rl1NuY2JkZnDTTRJdNGcOfPyxZBb/61/lz0GIjpZt4kTp8+efQ0SElMKYOBE6FBmhrChKZaekzOIRwHXFZfsq7qN9eylsl5UlT+433CBiMGeOPJGXlxo1xHn8z386BWHZMll3YdAgCSEtD23bymtuLnz1lSxctmOHtJeeLjOFs86YGqgoSmWhJB9BFxWBiiUiAs4+W6KLHnlEZgd798qax/tdyLCIjYV77xWT0fXXw2eficnogQdcqEyK+Ca6dJEifNHR8PXXIjxt20r2sqIovoGuR1AJuegiyT+49FLJHn7xRfjkE9dLX8fFSdG5RYvgttukkuqgQXD77bBmTfkT3wICxIwVHy9t/fvf0t5338FVV+nKlopS2VEhqMRUqwZ9+8KPP8Ktt0qkzjvvSI0gVxLTYmLkhr1okWQ/b90q1U4HDHBdcOrUkZt/rVrS3+PHpRz2b7/B4cPwwgtiRlIUpfLg7eqjSilo1UpeGzeGFSvEJp+RIWUnsrLkxl4ewsPFf9Cvn/goZs4Uk1Tt2mJC6t27/A5rY8Rc9OSTIix//gnLl4vpKDcX7rxTZiBHjojpSlEU76EzAh8iJEQylL/7TrKI09PlJl7K2ngltnvNNbJK2jPPiBBMnAg9e8oymp+MTOHrceW/SFCQLKN54YVSPK9xY4k4euYZSEyURDtFUbyHzgh8kNBQCQGNjZUn61q1xPwC8MsvEooaUA6JDwx0hjb/8ouYoebOFYE4//zWDBoE55/vWnmMJk3k1VpIShIh27lT2vzgA8mInjnTtWgpRVHKhs4IfJjISBgzBu6+W+oLLV0qzmBXIoEctGghhegWLpTcg19/jWLoUDEZvftu+UpgF8QYMXkNHy5O7P37JXR22zbYuFEipiZPlnBXRVE8i84IqgABAVIOYvx4qT1Uv77MEN57T2z8N99c/rbj4sSx3LPnGr7//hLee0+u88ILkqDWt6+Ihqv9j4oS/0FenvR9yxYRuc6dJVEtJkZEoXNnOVdRFPehQlCFCAmRNQpAInYmTBDH7IED4rjNyCi/Yzkk5ARXXw1XXw0bNogZZ8kSKU7XqpUIQrduUubCFQIDnRnMH34ohfS++Qb27RNBe/xxCUc1Rr5jzZquXU9RFDUNVVmio6WsxNy5MkP49VepZzR/vuttt2oFjz0mpqgHH5SCef/9r8wQxo8X/4I7cPhC4uOhaVNZle2ssyQh7pVXZP+nn8q5ubnuWQBIUfwRnRFUcRxP19WrS6G4Fi3kKfvAAXnt0qX8jtnq1aUcxvXXS8G8998X4Xn3XVnb4KqrRHzKOwspSEiIrOUAcsM/cECuffiwlLn49FN49VX49lsxZymKUnpUCPyEunVhyhSxwe/dKwlkixeLHyEhQc4pbzSQMXDeebIdOiS2/IUL4dln5Sm+SxcxKXXqJKGkrmKMCM39+cXQHetCN28uyXd168qiPX/9JTWWXF0ESFGqOioEfkZgoMTuz54NX34pkUd79shCOa1bi2PYFaKjJSS0Xz8JBV2wQHwJn38u9vyePWVzZ1G6atVk5tGjB+TkSKb0rl0SiZSaKt/3v/8Vk9bdd7vvuopSVVAh8FOCg6WWEchTfO3aEr2zd6+YfL7/Hi64wLV4/qZN5an97rtlzYUFC2Q1tTfflHyC7t2lEmqdOu75TiDfq0YNue6JE2I62rVLTEaHDolYxMbK8euvl8+K4u+oEChER0sUUE6OmFPmzBEn8KhRYud31bQSHOxMVDt4UOoZLV0qxfRefBHOPVdE4Yor3FtuIiBAZjyRkVLR9ehRKXWxdq040uvUkWsbIzOGu++W6q+K4m+oECgnCQ6WCKNhw+SG3Ly5zBBWroS1a5vz+ONihnGFGjXguutk27lTbshLlkgpi2eflTIU3brJbKW8dY6KIyxMtthYiZ7KyBCfwsaNks184YWSpb1jB7z1lswa6tZ1bx8UBZAFy9esgWPHJMV+7Fjo399r3VEhUE4jNFSWvAQpaieF7iI4fFj2/fmnhHVGRrp2nTp1pAz2oEHiT1i6VLbRo8WpfMEFcPnlIgruiDwqSECAzIRAbv7Ll0tF13XrZDnOF16QonthYZLHsHKl5C+4+p0VhVmzYMgQEQGQ7MkhQ+S9l8RA8wiUEomIgHHj4OWX19GunUQdDR8Ojz7q3jWLmzaVmciCBfDGG3DjjVJy4oknxI8wdKiYr1xZoKckHDOFWrVktrJ8uZjKvvlGTGVTp8L27ZBzcTIbb/6EJ57QvAWlnIwYcXod+exsry7coTMCpVQYI9E3CQkScrp3r9j7Dx+WWe2wYRI+6ioBARK91Lq1lJz45RdJIPvkE/jf/+Cpp2QFNMdMITHR9WsWRcEn/2HDZNayfTvEZ8CazLZ8O1vyGGJjpYpqdLSU8FaUM7J1a9n2VwAqBEqZMAYuu0ze5+TAF1+IrT07W5K8jhyRKJ127cpXAbXwtVq2lG3oUDEfffaZJI+NHy9bs2YiCF26SLKcq9csDocwBAfDc/XG8NMLfdmyRfq0erX4M665RkJkb7tNHOP33eeZvig+ToMGYg4qar+XUCFQyk1wMHTtKjfDjAwRgNGjxZTy/vvyu87Lc09JaWMk9+CssyTXIS1N7PYrV0pG8SuviFmnc2cRhg4dXK97VBIFZwwTJshKbDt2yFjs3i0lPX78UXwbl18ODz8smd3WioCGhHiub0olZ+xY8QkUNA+Fh8t+L6FCoLiMMVIRNCpKwjT79oV69SRR7amnxPH71FPuvWZSkmw33ywL9KxaJaKwbJmUuQgLkyigLl1kDej4ePdevzAhIbLFxEj57NxcKZS3caM4xXfulLDV48dFPF9/XfwgR4/KONWvrxnQfoPDITx4sDiMGzbUqCGlahEV5ayAmp0tEUf798vNLiRE6hB17uzeeP2YGOjVS7bjx2UFN8dsYcUKOadpUylx0amT+Bg8/UQeFOQUx2efleS27GwxA/fqJTOoL76ATZtkvYdFiyS57a+/ZFZx4YUSvaVUUfr3l2ksSJial/GYEBhjXgN6AXusta2LOG6ASUBPIBsYZK1d56n+KBVPeLhEHFkrWb0//iiRc8bIQ1BoqNz0WrRw39NwSIjzhv/QQ9L+V1/JNnu2ZDWHhclKa506yWyhXj33XLskAgIkAsvh8wARrYgI8X8cPSoO8U8+Ed9HaqrkcaSmyrgNGOBZU5fi33hyRvAGMBmYWczxHsBZ+duFwJT8V6WKYYw8tXfpIk+8+/ZJxNHCheJTeOYZca5a617zSEG/ws03yxP5d99JjaUvv5QnchCzTKdOkrfQvr37E9mKIyTE2T+QWcOll0oE0o4dYk564w0Jqb3iCknGe+cdWQ9i0iQ1JSnuw2NCYK1daYxJKuGU3sBMa60F1hhjYowxta21uzzVJ8X7VK8uW6NGYjuPjJSn8z17RBhSUiRmPyLC/dcODxcx6tJFPm/bJoKwZo3cbOfMkSf3Fi2kT+efL2akinoSDwiQcNju3Z37hg+XvIZNm0QoP/tM/A4//SQRSg89JJFaH3wg5+/eLWW4dc1npSx400dQF9hW4PP2/H2nCYExZggwBCAhIYGUCrKpZWZmVti1KjueGov27eVJODQUYmLiiY+vSXr6RtLTYdGiegQFnaBXr51uvy7IE7mjBlJOjuHXX6NYv74G69fHMGtWFDNmBBAUdIIWLQ7Tps1B2rRJp86xTIJNHrt3p3ikT0URHu5M3hs0SAQhLU0yvAMD6xMUFMiiRWkEBsI997QnNvY448b9iDHwxRdx1Kt3hEaNsjzSN/0/4qSsY9E2PR2A1EowfsZ6MD0yf0awqBgfwSJgnLV2Vf7nT4GHrbVrS2qzQ4cOdu3aEk9xGykpKSQnJ1fItSo7FTkWWVmSsHbDDeJ0fewxmUWsWCFO5oqo/3PkiNjnv/lGon1++UVuwBEBR+gY/h0tB3WmbVspbV0ZQkFzcqTPy5aJqHbsKK+9e8M//gETJ4qg3H235DtceaV7rqv/R5yUeSwc51aQEBhjvrPWdijqmDdnBDuA+gU+18vfp/g5ERGyrVkjEUcZGWIXf+wxiUi65x65yW3b5rkcnGrVnE5nEGf3unWw6emlrDh8Lp9Olv2hoSJOjoV5zjnHM2atMxEcLNt11zn35eTA9OkSofj11yIUCxeKqDqS7zp1Euf0gAFy/qZNkqQXHFzx30HxHt4UggXAMGPMO4iT+JD6B5TC1KwpW1IS/PCDiEJOjkQB3XOPLGbfq5fTdOKOFdCKIjpaMqpvnz2J4znpfP9cKqmpsm7D999LXsCrr4ptvnlz8S20ayev7i6YV1qCgyVstiALFkiE0h9/yKzr3HPFP7N6tTjy+/aF116Dm26S/TNnikh4MelVqQA8GT46G0gG4owx24HHgWAAa+1UYDESOroZCR+9xVN9UaoGrVo539evL0/pHTpIFFJKiizF+cYbFXPTqlFDhMFRbiMrS8I8HcLw/vuyCA+IY/ycc+Sm27o1NG7suVIYZyIw0Dnjio2V+k3WSihrYCD8+99i6vrsMzGJjRolM56wMHGsjx0rYbhNm0pJkf37Q9we7aVUPJ6MGrrxDMctcJenrq9UbRo2hJEj5f2xYzIj2LxZTDV79sCHH0p0zbPPem6WUJCICLHLd+won48fF3NWaqpsK1bI07jj3LPPdgpD69bemzWA3MRDQ0Vcr7/euf/vf5fIKWMk7HbjRsmY3rhRZmZz58ITT1zE5s3y9/j6a/mut91WOfwmSunRzGLF5wkNFefnlVfKjerQIckR+P13eWoNCJBEtqgo5zoLniYkRMxCbdvKZ2vFp/HDDxL6+cMPYoI5cUKON2ggswbH1qRJxQhYSQQGSv0mB926yZabK1VnGzaEO+7YxK+/NuO33+T7LFok4a9RUeKg/uILEUFjxOEeFHS6uUrxPioESpUiKEh8Ck8+KTffzEwRgx07xIS0Z4+YOV5+WZ52KyrgxRi52TdoID4NkAS3DRucwvDVV/DRR3IsNFR8Da1aObORGzasHPkBQUGS/9G2LdSuvZOEhGaA5Dz07y8CfOKEzBqqVROzXXS0rGHxxx9icgoLk3wRgHvvldfcXO+Ln7+iw65UWYxxJrB9+qk4SQ8flgrAX3whN6PWreWm+8orcoNu1qzi+hceLj6ODvkBfdZKNvGPP8LPP4sJZt48ySYGuam2aCGi0KqVbPXqec/fUJiQkFPLdQwcKK+5uSLI114rUWBffSX7HQl83brJDOKaa8SfMnu2CN6774r4OcxtiudQIVD8BseaxfHxYqZJT5en1lWrJDM3KUmOZWWJiaN371MXvqmxZBYRP64hMucYkb2S2HHXWA72cF/FSGMkR6JuXWd2cV6eJI9t3OjcPvjA6YiOjHSKQ4sWMouoX79yzBwcBAXJ1r79qfsnThT/zq5dUoyvXTv5PsuXi+jdeaeIxKRJ8nfr2VPE5P775d+vWCHfOSGhwr9SlUOFQPFLAgIkaiY2Vp46r75aROHwYYn4efVVuTGFh+evIfL2LAasHEJAjqwzG/rXFhqOlXVm3SkGhQkMFH9BkyZOk1JurmQVb9ggwrBhg8wacnLkeGios4ZR8+Yyy2naVL5LZSIgQG741arJZ8eyvSDf5fXXJffh++/lOwcGiolv9WpnVvhDD8kKj8ZIPak77hAnd16ezK7q1q08M6bKjAqBouA0IdWpI9FIgwfLzWXvXrnJ/t/yEQRy6jqzgUezqfviCI8KQVEEBTlv9L17y76cHBGHTZtkUZxNm8QcNneuHDdGZgrNmp261apV/tDPZkOSScpJJ+31VLd8r4IEB59eFXbCBGeo66FD8PTT4g/68ksxOa1bJ5FLderI58suk7Uh/vUv8Q/9738iFK1bSxvHj5+6wJA/o0KgKEVQu7a8JibmOzVf2QpFVGMJ2b2Vd9+V9wVDLyua4GDnzd0xc7DWuVrab7/J68aNUuraQUyMCIpj1tGkieQ5VNYbpCPUNTRUFvhxUKuWzOSOHxcR2LdPymlUqwaffy7f/Y03JGQ3LEyc89deC0uWwN/+JsL5zjuy+l2dOjIDCQjwn9mECoGinAFjKHad2ZzEBqxaJU/kXbvKTWbaNHF69ulT8X0tiDEiZImJUt7aQWam5Fw4Zg6bN8P8+WKGcZCQcKo4NG0qPpTKvCZCQIDTD+SocOsgLg4+/liEIi1NZhSDBkk59OXLnclynTqJWWnhQnjgATFLNW8ukU4rVsiMIjJS/t5BQS4m0lWCYnMOVAgUpTQUs85syPixrL5JHM9Hj4op6aefxAm6d6/cLEaNktXHevb0Wu9PwRH66chxAAn33LVLwjt//925ffut0/dgjJhrHOLQ+cBlNA76icCjlVsg4NSZBEhm+LnnOo9feaWYkvLyxOeSkyN/r02bRP/nzpWQ427dxBz10kvi7P7rLxnPZcvgww8bcOmlcq3sbBkTX5lRqBAoSmkosM6sPXYMU2CdWYPcWEBMSuvXy9N1ZqbcRLKzJZdh7155f/vtEnPvcGrm5Hj/RhoQ4IxYcqzXAGIi2b79VHH4/XcJv30173EATBcxpzjWkS64xcT4RvkJx9KiDjp3lg3ExDZwoAQU7N0rghkZKf6ZVatk7F57DT7+uB4PPiiZ4yNGyAzkzz/Fsf3225LDMny4tHnwoLRRWYr7qRAoSmnJX2f2UHo6MampxZ5mjETohIdLOGpqqpgkMjPFHHP++WK73rtXbPf33SertF1yiZxz8KA8eVeGG2hQkPOmfvnlzv3HjgGDB7Mhqybf9niatDQxuaxdm38sn+hoicoqLBB16vhO8pgxcjMvGKbqyLIGEYrhw+GGG9awbdsl5OaKn+Xqq6VmkzESAbVzp4QFR0RIMMLevWKSCgkRx3dYmDM09pdfREQLhi97Eh/5UyiKbxMSIqGqnTrB4sXOxexjYyWqpUEDeWJMSZGbwvTpEt2ybZs8gXfu7P1ZQ0FCQ6FZ+O80D/6O84Y8fXL/iRNiLnEIg2NbvdpZawnkSbh+fdkaNHC+1qsn4ukrJhVwCkVY2ImTNaMKrvdgLYwZI7PEXbtkBtixozwcrFwp53z0kQhE167yAPGPf4hgvvuutH3bbdCmDTz4oJz/yScyXu5KgFQhUBQvEBAgpoELLpAN5MZQr56IwznnSF7DvHliVpg7V/Z/9ZU8Ld5+e+V8og4IkKf9OnXgootOPebI6nasrpaWJkL31Vfy3R2Ehso4OITBIRT160t0kC+JBIhQONaLcHDNNaeeM2GCs07W/v2SExEcLDMGEN9NUJCMVXi4RDz16SNrSYSEyIPCzTc7hWLyZBn/du3kc8HxLYpK+FNSFP8kJMRZdA7kSbJDB7np160rfoZ162TWcN11kmD15pvyBD5unPybQ4ckYqYy3iyjok79fg7y8mQ2tHWrCMO2bfL+zz85GZHlwFElteBsok4dGZ+EhMopjqUhMFA2hzPbYXZy8PzzIhS5uRKY8NRTcu7atTJ+cXEyhp98IsJz991iZoqPl/Fr0gQgoFhjo48Om6JUfYwRG7vDedu4sdThycpy1k2KipInyD175PxHH5UbwyuvyE1x9Wq5STRv7t3vUhKBgeJkr10bLrzw1GN5eZILUVAktm2TJ+QvvpAbY8F2EhNFFBzi4HitW9d3HNfF4SjVAeJnKsj48fLqCD5YsEDMdD//LH6ngQPhjTeqF5sdokKgKD6EMWJSioyUG/yLL8rM4ehRsUEPHiwmpawsMQeMGSNrHzz2mDxBvvCCLKl5xRXSXmVfVCYw0GlqKlx8LjdXRGLnTik9sWOH8/3KlTKDKkh4+KnCUKeOmJ4c7VcmH0x5ccws6tRx7ouLg2HD4I03iv9LqxAoio9jjLNmz7Bhzv3Hj0uZicOHxb+wZ4+YWoKCxPForQRC3XKLLE1prZRoaNlSZiKVnaAg50298BMyiDPeIQwFxWLbNvmeR4+een6NGjIrSUiQ18TEU1+joyu3aLqCCoGiVFEcPgcHLVuKiSU7W26Cu3bJzCAhQUIZd+8WIbn3XnFGZmdLElXfvmJasvklNnzlZhgeLhnRRS2EY63MGAoKxF9/OSOevvrqdKEIC3NmahclGNb6yMAUgQqBovgRjmglh2lpzhzZn5MjN8b335cbXFiYZEgvXy5P27GxkgPx+ONij27XDv7KieW7w62pkyWhj76EMZIhXLPm6c5rEKE4dMgpDrt2nfq6adPppqeAgEuoVUuEIT5exjE+/tQtLq5yOrQrYZcURalogoPlxnXttc59HTrArbfKzOD4cXFEXnaZOKj37oVvDrZn6PYRTP9Nnrq//16clI88IjfYjAynT8PXMEacyzExsuZBURw9KrMoh0D8/vsWMjKS2LVLRHPlylOT68BZ/twhFLVqnfre8VrR/goVAkVRiqVgfR7HutAgM4iMeqtJCh1E485vnCyj8dtvIhwnTkgy1BtvSB2emBipW7R5s/gjAgIqv6P6TISFSdZ0w4byeffuNBISkk4et1b8M3v2FL1t2SJjkpl5etvR0afPJmrVkhmFY6tRw30LEKkQKIpSZoKDITY0m06RqcTkh6aed56Yjo4dc9rX69QRETh8WGrvpKRIjHxAAEyZIuan2bPl89dfy8yjYK0jV2g2JBmATdNS3NNgGXGE/0ZHS6nv4sjOLloodu+WmdfGjaeboUBEIDZWBKJmTadQFH4tjWCoECiK4jYCA511lrp3dy65CVIKe88emWFkZ0u5jRo1JEEqL0/q8aSni1M7JETqLxkD//mP/Pt168QsVZTz15cJD3fWYCqOnBzJF9m7V9ZaKPi6f7+Yp376SepUFSYgQISiJFQIFEUpO7NmwZo1RB87Jnew/EqsJREcLKGeDh5+2Pn++HFx2v71lzzJZmSI09VRzhvgySdl35gx0tbo0VJr57bb5Pg330j0Tv36bv2mlYLgYGfEUkk4BGPfvtMFo2Ctp8KoECiKUjZmzZK1GY4dw4AYux0LDp9BDIojJMSZE+Dg1VflNSdHBOHttyVprl49sasHBkpS2d69Yo9/6CFxZt93n7TX/ddJ9Iv9hPxq0qxaJbOJiqro6Q1KEgwVAkVR3MeIEacu0APyecSIcgtBSTgKtnXqdOr+5ctFAI4dk23ePDmvVi0xkVQPzCLQ5pxcB+K++yQKqn9/EZdhw6Ty6xVXiE/js88kLLYqC0VxqBAoilI2tm4t234PYoxzecqC6yUALD3v/7DAseRHyMyU6CWH83bLFonEccwotm6VleQeekjWMN65U4q2PfaYVPbctw+WLhXRSEyUf2eM+6J2vI1HhcAY0x2YBAQC06214wodHwSMB3bk75psrZ3uyT4pikukpJCakkKyt/vhTYpZv5kGDSq+L2fA4BSKghU9W7SQFeKsFf9EVpYsXRkdLbOKY8dkBhIeLkKxbp0sTVmnjtz816+XCdCLL0q01ObNMGdOU26/XQTm0CFxfNetWzkTyArjsS4aYwKBF4G/AduBb40xC6y1Gwqd+q61dthpDSiKUjkpZv1mxo71Xp/KScG1jAsWtWva1JkzkZcnq8f17y8iYa1z+cq6dcWx/fPPsGxZ4sl1Bj7+WITj7bfFef3FF/Dhh1I+OjpaFhv680+JpHK06c2cCk9q1QXAZmvtHwDGmHeA3kBhIVAUxZcoYf3mqogjJNaROAYiAD16OD8nJ8MVV6yiXbtkcnMlXLN+fQmFPX5ctuxscXIfPw4LF8Jbb4kDt1o18b/Pny/7Q0Ml32LDBrjrLml/xw7xa5QUYuoKnhSCusC2Ap+3AxcWcd61xphLgE3AcGvttsInGGOGAEMAEhISSElJcX9viyAzM7PCrlXZ0bFwomMB1K1L2xYtyMvL48cXXpB9lWxM2qanA5BaQf3Kysrk+++d12rXTsxDIE/+l17qnE307RvIJZeEERaWRV4eNGxYk+TkGA4e/J0TJ+DLLxuzalU811yzBoAXXmjOunWxvPHGV/lrIDdm27YIRo/+EYDVq+M4ejSQyy/fDcDBg8GEhp4gPDyvVH33tvVqITDbWnvMGHM7MAPoWvgka+00YBpAhw4dbHJycoV0LiUlhYq6VmVHx8KJjkU+MTGkp6dX3rHIX0DYF+4Xf/+74119rJXZxvHjkJubTE6OFPXbtg2aNk3myBGZkZw4AcYkY4yYog4fhr/9rSXGwHPPSajt1Knio5g4seTre1IIdgAFUzvq4XQKA2Ct3V/g43TgaRRFUfyYotY4LhwR5ViL+MQJMRl99pmEzFarJhFNQ4eKGSoyUgRBlvs8caK4a3pSCL4FzjLGNEIE4AbgpoInGGNqW2t35X+8Gtjowf4oiqJUKQICnM7u6tWd+wcPPvW8iy8GYw4XUd5O8JgQWGtzjTHDgGVI+Ohr1tqfjTH/BdZaaxcA9xhjrgZygQPAIE/1R1EURSkaj/oIrLWLgcWF9o0q8P5R4FFP9kFRFEUpmQBvd0BRFEXxLioEiqIofo4KgaIoip+jQqAoStUjf70EVqyQdNxZs7zdo0qNCoGiKFWLAuslAM71ElQMikWFQFGUqkVJ6yUoRaJCoChK1aISrZfgK6gQKIpStShuXYRKuF5CZUGFQFGUqsXYsVI3uiA+ul5CRaFCoChK1aJ/f5g2TQrwgCwkMG1alV0vwR14uwy1oiiK++nfH155Rd5XsnUSKiMqBIqilA9dv7nKoKYhRVEUP0eFQFEUxc9RIVAURfFzVAgURVH8HBUCRVEUP0eFQFEUxc9RIVAURfFzVAgURVH8HBUCRVEUP0eFQFEUxc9RIVAURfFzVAgURVH8HBUCRVEUP0eFQFEUxc9RIVAURfFzdD0CRVGqJrogTanx6IzAGNPdGPOrMWazMeaRIo6HGmPezT/+tTEmyZP9URRFUU7HY0JgjAkEXgR6AK2AG40xrQqdNhg4aK1tCkwAnvJUfxRFUZSi8eSM4AJgs7X2D2vtceAdoHehc3oDM/Lfvw9cbowxHuyToiiKUghP+gjqAtsKfN4OXFjcOdbaXGPMIaAmsK/gScaYIcAQgISEBFIqyPaXmZlZYdeq7OhYONGxcKJj4cSXx8InnMXW2mnANIAOHTrY5OTkCrluSkoKFXWtyo6OhRMdCyc6Fk58eSw8aRraAdQv8Lle/r4izzHGBAHRwH4P9klRFEUphCeF4FvgLGNMI2NMCHADsKDQOQuAgfnv+wKfWWutB/ukKIqiFMJjpqF8m/8wYBkQCLxmrf3ZGPNfYK21dgHwKvCmMWYzcAARC0VRFKUC8aiPwFq7GFhcaN+oAu+PAtd5sg+KoihKyRhfs8QYY/YCWyrocnEUimDyY3QsnOhYONGxcFLZx6KhtbZWUQd8TggqEmPMWmttB2/3ozKgY+FEx8KJjoUTXx4LLTqnKIri56gQKIqi+DkqBCUzzdsdqEToWDjRsXCiY+HEZ8dCfQSKoih+js4IFEVR/BwVAkVRFD9HhQBdQKcgpRiLQcaYvcaY1PztX97op6cxxrxmjNljjPmpmOPGGPN8/jj9YIxpV9F9rChKMRbJxphDBX4To4o6rypgjKlvjPncGLPBGPOzMebeIs7xvd+GtdavN6T8xe9AYyAEWA+0KnTOUGBq/vsbgHe93W8vjsUgYLK3+1oBY3EJ0A74qZjjPYElgAE6Al97u89eHItkYJG3+1lBY1EbaJf/vjqwqYj/Iz7329AZgS6gU5DSjIVfYK1didS/Ko7ewEwrrAFijDG1K6Z3FUspxsJvsNbustauy3+fAWxE1lUpiM/9NlQIil5Ap/Af9pQFdADHAjpVjdKMBcC1+VPe940x9Ys47g+Udqz8hU7GmPXGmCXGmLO93ZmKIN9EfB7wdaFDPvfbUCFQyspCIMlaey6wHOdMSfFf1iF1bNoALwDzvNsdz2OMiQQ+AO6z1h72dn9cRYVAF9ApyBnHwlq731p7LP/jdKB9BfWtslGa341fYK09bK3NzH+/GAg2xsR5uVsewxgTjIjALGvth0Wc4nO/DRUCXUCnIGcci0K2zqsRG6k/sgC4OT9CpCNwyFq7y9ud8gbGmESHz8wYcwFyX6mKD0rkf89XgY3W2ueKOc3nfhs+sWaxJ7G6gM5JSjkW9xhjrgZykbEY5LUOexBjzGwkGibOGLMdeBwIBrDWTkXW2egJbAaygVu801PPU4qx6AvcaYzJBY4AN1TRByWAi4EBwI/GmNT8ff8HNADf/W1oiQlFURQ/R01DiqIofo4KgaIoip+jQqAoiuLnqBAoiqL4OSoEiqIofo4KgaK4QH41yj+NMbH5n2vkf07yctcUpdSoECiKC1hrtwFTgHH5u8YB06y1aV7rlKKUEc0jUBQXyS858B3wGnAb0NZam+PdXilK6fH7zGJFcRVrbY4x5t/AUqCbioDia6hpSFHcQw9gF9Da2x1RlLKiQqAoLmKMaQv8DVmNanhlX4REUQqjQqAoLpBfjXIKUpd+KzAeeMa7vVKUsqFCoCiucRuw1Vq7PP/zS0BLY8ylXuyTopQJjRpSFEXxc3RGoCiK4ueoECiKovg5KgSKoih+jgqBoiiKn6NCoCiK4ueoECiKovg5KgSKoih+zv8DeNaNaq0crfEAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Cute addition to the simple plot. Add a 1 sigma band\n", "# Code is not very pythonesque...\n", "# Note: temp3 is the covariance matrix\n", "dy = np.empty(len(xpl))\n", "for ii in range(len(xpl)):\n", " xx = xpl[ii]\n", " Bt = np.array( [dydp(0,p,xx), dydp(1,p,xx)] ) # 2x1 matrix\n", " B = (Bt.T).copy() # 1x2 matrix\n", " blah = np.matmul(B, temp3) # 1x2 * 2x2 = 1x2\n", " dy2 = np.matmul(blah, Bt) # 1x2 * 2x1 = 1x1\n", " dy[ii] = math.sqrt(dy2) \n", "\n", "a2 = plt.subplot(111)\n", "a2.errorbar(x, y, yerr=np.sqrt(e2), linestyle='none', color=\"red\", marker='o')\n", "a2.set_xlim(xmin, xmax)\n", "a2.plot(xpl, ypl, color='blue', linestyle='solid')\n", "a2.plot(xpl, ypl-dy, color='blue', linestyle='dotted')\n", "a2.plot(xpl, ypl+dy, color='blue', linestyle='dotted')\n", "a2.fill_between(xpl, ypl+dy, ypl-dy, color='blue', alpha=.2)\n", "a2.set_xlabel('X')\n", "a2.set_ylabel('Y')\n", "a2.text(0.96, 0.96, txtRes,\n", " transform=a2.transAxes,\n", " bbox=props, fontsize='large',\n", " horizontalalignment='right',\n", " verticalalignment='top' )\n", "a2.grid()" ] }, { "cell_type": "code", "execution_count": 9, "id": "5cfe4121", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABW7ElEQVR4nO2dd3gURRvAf5MKJLQECKElSO9ViqCAoIReREBRARE+RRQURBAVVFAUFQuKIr0GFIRQgoiQ0DuhdwiQEFqA9HbJfH9MgASSS7vLXZL5Pc8+d7c7N/tO9rLvzrxNSCnRaDQajSY9bCwtgEaj0WisG60oNBqNRmMUrSg0Go1GYxStKDQajUZjFK0oNBqNRmMUO0sLYGpKlSolPT09U+2LiorCycnJMgJZgII2XtBjLggUtPFC7o754MGDt6WUpdM6lu8UhaenJwcOHEi1z8/Pj7Zt21pGIAtQ0MYLeswFgYI2XsjdMQshLqd3TC89aTQajcYoWlFoNBqNxij5bukpOyQkJGAwGCwthsmIi4sjJibG0mJkCzs7O+zt7S0thkajSUGBVRRJSUksXLgQHx8fgoOD89XNKSYmhsKFC1tajGyRkJBApUqV6NmzJy+//DJCCEuLpNEUeAqsopgxYwYBAQFMnjyZWrVq5asbUkREBEWLFrW0GNkiKSmJEydO8M033xAVFcWwYcMsLZJGU+ApkDaK2NhYVq5cybfffkvt2rXzlZLI69jY2FCvXj2mTZuGt7c3CQkJlhZJoynwFEhFcfbsWTw8PHBxcbG0KJp0KFu2LK6urgQGBlpaFI2mwFMgFUVcXFyeXcMvSBQpUoTY2FhLi6HRFHgKpKIA9HJTHkBfI43GOiiwiiIj2rZtS8mSJYmLizN53zNmzKBp06Y4OjoyaNAgo20DAwPp3LkzJUuWpGzZsowYMSKVK2/btm0pVKgQzs7OODs7U6NGDZPLe+fOHXr16oWTkxMeHh4sXbo0zXZxcXEMGTIEDw8PihYtSsOGDfH19c3SeDQajfWhFUUaBAYGsn37doQQ+Pj4mLz/cuXK8fHHH/P6669n2Hb48OGUKVOGkJAQAgIC8Pf359dff03VZsaMGURGRhIZGcmZM2cyLcekSZOYNGlShu3efvttHBwcuHHjBkuWLOGtt97ixIkTj7UzGAxUrFgRf39/wsLCmDx5Mn379k1lZ8jMeDQajXWhFUUaLFy4kBYtWjBo0CAWLFhg8v579+5Nz549cXV1zbDtpUuX6Nu3L4UKFaJs2bJ4eXmleZM2F1FRUaxcuZIvvvgCZ2dnWrduTffu3Vm0aNFjbZ2cnJg0aRKenp7Y2NjQtWtXKleuzMGDB61mPBqNJutoRZEGCxcuZMCAAQwYMIB//vmHGzdupNu2a9eulChRIs2ta9euOZZl1KhReHt7Ex0dTXBwML6+vnh5eaVqM378eEqVKkWrVq3w8/PL8TlTcvbsWezs7KhevfqDfQ0aNMjUzf3GjRucPXuWOnXqPNiXmfFoNBrrwmKKQghRUQixVQhxUghxQggxMo02bYUQYUKIgOTtU3PLtWPHDi5fvkzfvn1p0qQJVapUSXdNHmDdunXcu3cvzW3dunU5lueZZ57hxIkTFCtWjAoVKtC0aVN69uz54PjXX3/NxYsXCQ4OZtiwYXTr1o2LFy+m219KxTZ16lSmTp1qVLFFRkZSrFixVPuKFy9ORESEUbkTEhIYMGAAAwcOpGbNmpkej0ajsT4sOaMwAKOllLWBFsDbQojaabTbLqVsmLx9bm6hFixYwPPPP0+pUqUAePnll82y/JQZkpKS8PLyonfv3kRFRXH79m3u3r3Lhx9++KBN8+bNKVq0KI6OjgwcOJBWrVqxadOmdPtMqdjGjRvHuHHjjCo2Z2dnwsPDU+0LDw83GvmdlJTEq6++ioODAzNmzMjSeDQajfVhMUUhpQyRUh5Kfh8BnALKW0oeUDmSVqxYgb+/P2XLlqVs2bJMnz6dI0eOcOTIkTS/06lTpwceR49unTp1ypE8d+7c4cqVK4wYMQJHR0dcXV0ZPHgwGzZsSPc7QgiklDk6b0qqV6+OwWDg3LlzD/YdOXIk1XJSSqSUDBkyhBs3brBy5cpUObSyMx6NRmN5rCLXkxDCE2gE7E3jcEshxBHgGjBGSvnY4rgQYhgwDMDNze2xdfrIyMhU+06fPk1UVNRjyyd//vknNjY27N69GwcHhwf7Bw4cyOzZs/nyyy8fE27FihVGx5bWEo3BYMBgMBATE0NsbCy3bt3Czs4OO7vUl8PR0RFPT09++OEH3n33XSIjI5kzZw61atUiIiKCe/fuceDAAVq3bo2dnR0rV65k27ZtfPHFFxkuDQGMHj06XRlT0q1bN8aPH8+MGTM4duwYa9as4d9//03ze6NGjeL48eP4+PhgMBhStcloPI8SFRXF/v37CQ0NzXAsj17jgkBBG3NBGy9Y0ZillBbdAGfgINA7jWPFAOfk952Bcxn116RJE/koW7duTfV537598n//+99j7Tp27Cjff//9x/YvX75curm5yYSEhMeOZYeJEydKINU2ceLEB8e9vLzklClTpJRSHj58WLZp00aWKFFCurq6yhdffFFev35dSinlzZs3ZdOmTaWzs7MsXry4bN68udy0aZMMDw9P99xeXl7Syckpzc3LyyvN74SGhsoePXrIIkWKyIoVK8olS5Y81ueUKVNkYGCgBKSjo2OqfhcvXvygrbHxPMrAgQPl0aNHM/U3ffQaFwQK2pgL2nilzN0xAwdkevfp9A7kxgbYA/8A72eyfSBQylibnCiK/IIxRZGX0IrCOAVtzAVtvFJaj6KwpNeTAOYAp6SU36fTpmxyO4QQzVA2lYzXITQajUZjMixpo2gFvAocE0IEJO/7CKgEIKX8DegDvCWEMAAxQP9kzafRaDSaXMJiikJKuQMwmvVNSjkDmGGsjUaj0WjMi47M1mg0Go1RtKLQaDQajVG0ongET09PChcujLOzMyVLlqRLly5cvXrVpOfIbNru+5w6dYpnn32W4sWLU7VqVf7++2+j/WUU22FumV955RXc3d0pVqwY1atXZ/bs2amO61TjGk3eQiuKNFi7di2RkZGEhITg5ubGO++8Y9L+M5u2G1RwXo8ePejatSt37txh1qxZvPLKK5w9ezbd/t5///1MJe3LbJrxrMo8fvx4AgMDCQ8Px8fHh48//jhVBlmdalyjyVtoRWGEQoUK0adPH06ePGmyPrOSthtUFPm1a9d47733sLW15dlnn6VVq1YP2qfVX6dOndLtLzdkrlOnDo6OjoBKKSKE4MKFCw+O61TjGk3eQisKI0RHR7N8+XJatGiRbpusphnPSdru+0gpOX78eLr91atXz6Q33uzIPHz4cIoUKULNmjVxd3enc+fOD47pVOMaTd7CKnI9WRs9e/bEzs6OqKgoSpcuzT///JNu26ymEs9q2u4aNWpQpkwZpk2bxnvvvcfWrVvx9/enXbt26fZXrFixdPvr2rUrO3bsACA2NhaAH374AYDWrVunOZ7spBr/9ddf+fnnn9m9ezd+fn4PZhigUo3PmjWLYsWKkZiYyMCBA3WqcY3GitEzijRYvXo19+7dIzY2lhkzZtCmTRuuX79ukr6zmrbb3t6e1atXs379esqWLct3331H3759qVChQrr9RUREpNtfVtOMZ0fm+9ja2tK6dWuCgoKYOXMmoFONazR5Ea0ojGBra0vv3r2xtbV98BT+KFlNM57VtN0A9evXx9/fn9DQUP755x8uXrxIs2bN0u3v2LFjRvvLKtmROSUGg+GBjUKnGtdo8h5aURhBSsmaNWu4e/cutWrVSrONr68vkZGRaW6+vr6PtXdycqJ37958+umnREVFsXPnTtasWcOrr76arhxHjx4lNjaW6Ohovv32W0JCQhg0aFC6/W3YsMFof/fJrNdTVmS+efMm3t7eREZGkpiYyD///MOyZcto3749AKVKlaJy5crMnDkTg8HAvXv3WLBgAfXr189QDo1GYxm0okiDbt264ezsTLFixZgwYQILFiww6RP6r7/+SkxMDGXKlOGll15i5syZqfrv1KlTqtoXixYtwt3dnTJlyvDff//x77//plrzf7S/77//Pl15s1toyZjMKeUVQjBz5kwqVKhAyZIlGTNmDD/88APdu3d/0NeqVavYuHEjpUuXpmrVqtjb2zN9+vTs/TE1Go3Z0cbsRwgMDDT7OVxcXFi9enW6xx+diUybNo1p06Zluj9jRua0ZjmZwZjMKfssXbo0/v7+Rvtq2LChdRRj0Wg0mULPKDQajUZjlAKpKBwcHB64hmqsl9jY2FQlaTUajWUokIqiSpUqXLp0iaioKEuLokmHu3fvcv36dTw8PCwtikZT4CmQisLZ2Zn27dvz6aefcuvWLUuLo3mEGzduMHHiRDp37kyhQoUsLY5GU+ApsMbs8ePHM336dPr164ezszP29vaWFslkREdHU6RIEUuLkS3i4+OJioqia9eujBw50tLiaDQaCrCisLe3Z+zYsbz33nvcuHEjX6W53rNnj9H8VNaMnZ0dZcuWxc6uwP40NRqro8D/N9rb2z9Ih5FfCAwMxNPT09JiaAookbGRhISFcDPiJrcjb3M74jZ3ou8QFh1GeGw44THhRMVHERMfQ0xCDHGGOAyJBgxJBhKTEhFCYCNssBE22Nna4WDrgKOdI5FhkXic88DZ0RlnR2eKFy5OSaeSlChcAldnV0oXLU2ZomUoU7QMRRzz5ozaWinwikKj0WSN8Jhwzt08x9nrZzl/6zxX7lzhcuhlrty5wrV714iITTuOx0bYUKxwMYoVKoaToxOF7QtT2KEwjnaOODk4YWdrh62NLaBygiXKRBKTEokzxBEWE8bt6NvcDLxJVFwUEbERRMZFpitj8cLFKVeiHOVKlKNiyYpUcqmEZylPnij1BNXdqlO2eFmEEGb5++RHLKYohBAVgYWAGyCBWVLKHx9pI4Afgc5ANDBISnkot2XVaAoicQlxHL92nIArARy/dpwT105wPPg4IWEhqdqVLV6WSi6VqFe+Hl51vHAv7o57cXfcirlRqmgpSjmXwtXJFSdHpxzdnP38/Gjbtu2Dz4ZEA2ExYdyNvktoZCi3Im5xM+ImNyNucu3eNa7du0bwvWD+Pfkv18KuIaV88F1nR2equVWjVtla1C5Xm9rutalbvi5VSlfBxqZA+vgYxZIzCgMwWkp5SAhRFDgohPhXSpmySlAnoFry1hyYmfyq0WhMSGJSIievnWTPxT3subiHA5cPcDLkJIZEZbsr4lCE2u61eb7289R0r0kNtxpUd6vOE6WfoLBDYYvIbGdrh6uzK67OrlQtU9Vo23hDPFfvXOXCrQucu3mOczfOcebGGXZe2MnSfQ/L+hZxKEK98vVoWLEhjT0a09SjKfXK18PeLv84u2QHiykKKWUIEJL8PkIIcQooD6RUFD2AhVI9CuwRQpQQQrgnf1ej0WQTQ6KBQ1cO4XfGD/+z/uw4v4PwGJVK3sXJhWaVm9G1flcaVWpEo4qNqFyqcp5+0nawc6BKmSpUKVOF5+s8n+pYRGwEp0JOcTz4OEeCjnA06CjLDyzn922/A1DYoTBNPZrSqmorWlVpxVNVn8LFycUSw7AYIuV0zGJCCOEJbAPqSinDU+xfB0yVUu5I/vwf8KGU8sAj3x8GDANwc3Nr4u3tnar/yMhInJ2dzToGa6KgjRf0mDPDtYhr7AvZx8GQgxy+cZioBBVwWqlYJRqUaUC9MvWoXao25ZzLWeX6fW5eYyklIZEhnA49zanQUxy/dZxzd86RKBMRCJ4o+QQNyjSgoVtDGpZpSFFH47VZsktujrldu3YHpZRN0zpmcWO2EMIZWAmMSqkksoKUchYwC6Bp06Yy5TomPL62md8paOMFPea0SEpKYu+lvfgc8cEnwIeTIWqy7uHqwUstXqJDrQ60qd6GssXL5pLEOcPS1zg6Lpp9gfvYfm47/mf98b3gy6ozq7ARNjR/ojledbzoUr8LjSs1NpmitfSY72NRRSGEsEcpiSVSylVpNAkGKqb4XCF5n0ajSYPEpER2nNvBXwf/YtXhVVy7dw07WzueqfYMQ58ZSue6nanmVs0qZwzWThHHIrSt0Za2NdryCZ8Qb4hn36V9/HvyX3yP+zJp7SQm+kykXIlydK3fle4NutO+VnsK2ef97AKW9HoSwBzglJTy+3Sa+QAjhBDeKCN2mLZPaDSPc+TqERbvWczSfUu5du8ahewL0bleZ15o/AKd63WmRJESlhYx3+Fg50Draq1pXa01n/X4jJvhN/E97svaI2tZuncps7bNwtnRme4NutO/WX861umIg13eTHJpyRlFK+BV4JgQIiB530dAJQAp5W/ABpRr7HmUe+zg3BdTo7FOwuLC+HHzj8zdOZejQUexs7WjU91OfN/3e7rU64JzoYJls7E0ZYqVYeBTAxn41EDiEuLYemYrKw+tZNWhVSzdtxRXZ1f6Ne3Hqy1epfkTzfPUrM6SXk87AKN/qWRvp7dzRyKNxvqRUrL93HZm+s1k5cGVJCQl0NSjKTNenkG/pv0oVbSUpUXUAI72jnjV9cKrrhe/vPwL/578l8V7FjN351x+9fuVqmWqMqT1EIa0HkLpoqUtLW6GWNyYrdFoMiY6Lpql+5by85afORp0lBJFStCtWjcm9p9I/Qq63rg142DnQJf6XehSvwvhMeGsOrSK+bvmM37VeCb6TOTFJi/yzrPv0PwJ6w0Ry7uO0RpNAeBWxC0m+UzCY5wHQxcOBeCP1/4g+Jtg3mn6jlYSeYxihYsxqNUg/D7w4+TnJ3mzzZusPbqWFl+1oNmUZizdu/RBkKM1oWcUGo0Vcjn0Mt9s/Ia5O+cSmxBLtwbdGP3caJ6p/kyurm3HxsLVq3D9OoSEqNfQULh7V23h4RAdrbaYGEhMhKQktQkB9vZqc3QEZ2e1FS0KLi5QqpTa3NygfHmoUEF9zsNxfVmilnstfuz/I1N6TmHh7oX8tOUnBswewKdrPuWjzh/xaotXLS3iA7Si0GisiIu3LvKV71fM3zUfgeC1lq8x+vnR1HKvZbZzxsfDyZNw6hScPQtnzsCFC3D5Mty4kfZ3iheHkiXVq5MTFC6sbv52dupGb2MDUkJCgtri4uDOHdVnRIR6Hx39eL+OjuDpCVWqQNWqUL061KkD9fPxxMm5kDPD2w3nzTZv4nPEhy/WfcGQBUP4Yt0X9K7Sm1atW1k8hYhWFBqNFXD1zlU+X/c583bOw87GjjefeZOxXmOp6FIx4y9ngfh4CAiAvXvhwAH1/uRJuF+ORQjw8FA36q5d1ftKlaBcOXB3V0//Li5ga5tzWaKj1ewkJASCg9V25QpcvKgU1fbtSqncp0yZFrRuDc2bQ7Nmasuj9bnSxMbGhp6NetKjYQ98j/vy+drP+X7f92y4soFvX/yWzvU6W8xTSisKjcaC3I64zZcbvuRXv1+RSN5u9zYfen1IuRLlTNJ/WBjs3Anbtqkb74EDSlmAuuk3bgxdukCDBlC7tnqKL5xLOf6KFFFbxXR0oZRKiRw/DkeOwMaNYQQEFGJVcmiuvT00bQpt2kDbtvDMM7knuzkRQtC5Xmc61e3El0u+ZOHphXT9uSvP1X6O7/t+T93ydXNdJq0oNBoLEJcQx89bfuaL9V8QGRvJwKcGMrHbRDxcPXLUb1ISHDwIGzeqbc8etc/eHp58Et59Vz2RN2+ubALW7MovhJrJlCsHzz8PTz55irZt3bh9W82Itm9XCvDbb2HqVKUkOnSA7t2hZ09l78jLCCFoVaEVH/T/gN/8f2OSzyQaft6QsR3HMrHbRBztHXNNFq0oNJpcRErJmoA1jPlzDBduXaBLvS580+cbaperne0+Y2Jg82ZYswbWroWbN9VN9sknYcIEaNdOKYb8skxTqpSaBXXpoj5HRSmlsWED+Piov8Gbb0L79tCvH/TqpewpeRUHOwfebf8uA5oPYMyfY/jK9ytWB6xm/uD5NKvcLFdkKCD+BRqN5Tl/8zydf+pMr1974WjnyMaRG1n37rpsKYmEBNi924VXXoHSpdVT9J9/KqWweLEyQu/dC59/rvblFyWRFk5O4OUFP/0Ely7BoUPwwQdw/jwMGQJly8LLL8OOHWo5K6/i6uzKvMHz8B3pS2RcJC2/askX674gKSnJ7OfWMwqNxszEJsQy1XcqU32n4mDnwPR+0xnRbgR2tln795NS3QTnz4dlyyA0tD4lS6qbYJ8+ap3eIW+mEjIZQkCjRmr78ku1DLdwodqWLVPeU++/Dy+9lHf/Vl51vTg28Rgjlo3g0zWfsvfiXhYNWURJJ/NNm/SMQqMxI9vPbqfBZw34bO1n9G7cm9NfnGZUh1FZUhJhYfDjj+om17Qp/PGHWlaZMuUY16/DrFlqDT+v3vjMhRDq7/XTT8qjatYsFecxaJBywf3uOxUnkhcpXqQ4C19fyC8v/8Kmk5toMrkJp0JOme18WlFoNGYgIjaCt5e8zTPTniHeEM+m9zaxdOjSLHkzHT+u1trLl4dRo5SxduZM5Qm0fDk89VSoVg6ZxMkJhg6FY8eUkb9WLRgzBmrUgCVLlME/ryGEYHi74Wz7YBvR8dE8/c3THAg8kPEXs4FWFBqNifE/40/diXWZ6T+TUR1Gcfyz4zxX+7lMfVdKZZj28oJ69WDBAujbV7m17tunFEdeNsxaGiGgY0f47z/YskUZxl95RRn+d+2ytHTZo0WVFmwfux1nR2ee/e5Ztp/dbvJzaEWh0ZiI2IRYRq8YTbvv2uFo58jOD3cyvd90nBydMvxuYiKsWKHiGp57TgXCTZkCQUEwdy40aWJ++Qsa7drB/v2waBHcugWtW8OHH6oo8rxGNbdq7PxwJ+VLlKfrjK4cCzpm0v61otBoTMDRoKM0ndyU7//9nrfavMXhTw/TskrLDL+XmKiMrPXqKVfOmBiYPVuluvjoI3B1zQXhCzA2NmpGceIEvPEGfPONml0cOWJpybJO+ZLl2fTeJpwdnenycxeu3btmsr61otBocoCUkt/8fqPZlGaERoXiO9KXXwb8kuEsQkpYtQrq1lVeSzY2yu5w4oRy6XTMvVgqDSpR4axZsG6dikNp2RJWr7a0VFmnoktF1r+znrtRd3npj5dM5jqrFYVGk03uRd+j7+99eWvJW7Sp3oaATwPwquuV4ff8/NSN6IUX1Oc//4SjR5UtwhQ5lDTZp0sXNZuoXx9691beZnmNhpUa8vNLP7Pt7DZ+2PyDSfrUikKjyQZHg47SZHITVges5usXvsZ3pC9uxdyMfufiRRUl3K6dsj3MmaO8cPr0KTiptfMCbm7K0N2zp/I2mzLF0hJlnYFPDaRHwx589PdHXA69nOP+9M9To8ki3vu8aflVS2ITYvEf489Yr7HYGLnTR0Yqe0OtWvDvv+rGc+4cvP66SsutsT6KFFEzvVdegY8/VnajvIQQgp9f+hmJZMr6nGs6iyoKIcRcIcRNIcTxdI63FUKECSECkrdPc1tGjeY+hkQDo1eM5qU/XqJxpcYc/PggT1V9yuh31q5VWVm/+koZq8+cUUojP2Q5ze/Y2iqPMy8v5Zbs52dpibJGRZeKDH16KPN2zSPoTlCO+rL0jGI+kNGi7nYpZcPk7fNckEmjeYzwmHC6zejG9/9+z9vt3ua/0f9RtnjZdNsHBak17u7doVgxlWdo4UIVPKfJO9jbq5nFE0/A4MFqdpiXGNVhFIZEA8sPLM9RPxZVFFLKbcAdS8qg0WRE4O1AWn3dis2nNjPr1VnMeHkGDnZph0RLqWwPtWurCOCvvlL5mVq1ymWhNSbD2RnmzVMuy+PHW1qarFG1TFUaV2rMnwf+zFE/Qlo4naIQwhNYJ6V8rBqHEKItsBIIAq4BY6SUJ9JoNwwYBuDm5tbE29s71fHIyEicnZ1NLbrVUtDGC+Yb8/Fbx/nY/2MSZSKTnp5Ek7LpR77duuXIt99WZ98+Vxo0uMfYsacpV858yYQyM2ab2Fgc7t3DPiwM+3v3sI2OxsZgQMTHIxITkQ4OJCVvCcWKEe/iQryLC4mFC1tdsQpL/66nT6/Ghg3uLF26l9KlcycqzxRjnh0wm2Unl+HbzxcH2/RzvrRr1+6glLJpmgellBbdAE/geDrHigHOye87A+cy6q9JkybyUbZu3frYvvxMQRuvlOYZ85rDa2ShtwrJah9Vk2dCzhhtu3SplMWLS1mkiJQ//yxlYqLJxXmMB2NOSpIyMFDK1aulnDRJyn79pGzWTMrSpaVUk5ysb8WLS9mihZSDB0s5fbqUAQHqPBbE0r/rS5ektLWVcvTo3DunKca8dM9SyRvIY0HHjLYDDsh07qtW7XMhpQxP8X6DEOJXIUQpKeVtS8qlyf/M9JvJiKUjaOLRhPXvrqd00dJptouIgBEjlP3hqafUa5UqZhYuIQEOHKCitzdMn65qnYaGqmNCqAX1J55Q/p2enqogQ+nSKrFR0aIqms/RUblcxcerFKoxMXD7Nly/rrZLl+DUKVUNaN481XfZsiq/SM+eqqB2ActI6Omphr54MUybZnUTrnTxLOUJqLrs2S2jatWKQghRFrghpZRCiGYom0qohcXS5GOklHyy+hOmbJhC1/pd8R7mnW6U9b59qq5BYCBMnKjcKM3m7nr5MqxfD5s2KSf/iAiqAFSvrizmzZpBw4YqF4hTxrmlskRQkMpUuGkT+Pqq5Eiurmrww4apcxYQunSBlStVZt+8MmwbkXNTtEUVhRBiGdAWKCWECAImAvYAUsrfgD7AW0IIAxAD9E+eImk0JicpKYmR3iOZsXUGQ58eyq8Dfk2zboSUKt33qFHg7g7+/iqhnMm5ckW53KxYobQSQOXKKufHc8+x08aGVr16meHEj1ChgiriMGgQGAwqGGTBAlUY45dfVJKkL7/M+0WqM0Hbtup11668oyhiE5SdLKuFslJiUUUhpXwpg+MzgBm5JI6mAJOYlMjQhUOZt3Meo58fzbQ+0xBprC1ER8P//qeWHzp3Vg/XLi4mFCQiQj2yLljw0HG/SRP4+msV1l2t2oOmCZZw7Lezg06d1Hbnjooe/PFH+OsvmDxZ/XHycR6SChXU6/XrlpUjK5y7eQ6AqqWrZrsPS8dRaDQWJ94Qz8t/vMy8nfOY2G1iukoiMFDlaFqyBD77TAXTmUxJnDihorrc3ZXDflCQKnh9/rwqRjF2bColYRW4uKgycUeOqNqjb7+t8pHkxTzdmcTeHooXV+acvMLRoKMUdihMJddK2e7Dqm0UGo25MSQaePmPl1l5aCXT+kxjTMcxabbbsUM90CcmKvuuV8a5/zJGSmV3+OEHVUmnUCG17j9kiLKM5xVraZ06yobx009qPa53bzUrKlTI0pKZHCkhKkr5BOQFpJSsO7qOdjXaYWuT/ZmenlFoCiyGRAOvzHmFlYdWMr3f9HSVxPz5qka1iwvs2WMCJXG/SlHDhtCtm8rr8dVXcPWqyhnRqlXeURL3EQJGjoTff1cG7/791V01n3H3rjLTlE7bCc7qOHTlEJduX6JXo5zZsrSi0BRIDIkGBswewPL9y/mmzzeM6jDqsTZJSTBunFoJeuYZpSSqV8/BSZOSlHG6bl2V+Ck+XtkiLl5UJ8oPxuBhw5Tv6Jo1ypCTz9i/X73WzZ6Xaa4z/d/pODs680LjF3LUj1YUmgJHUlISr819jRUHVvDti9/yQccPHmsTH6+cfL7+Gt56Sz0k56hW9ebNyoX1ftGJFSuUj+Vrr6mF7/zEqFHQogWMHp33kiNlgJ+fsuc/ZTwXpFVw8dZFvPd7M+yZYZR0ylmhda0oNAUKKSUjvUeybN8yvur9FaOfH/1Ym7Aw5dSzaJFy5PnllxzER5w9q9yjnntOWUAXLFDG3xdfzL/eQba2aint1i2lYfMJSUng7a1ml6YOVTE193/nhewL8f5z7+e4P60oNAWKz9d+zoytMxjz/BjGdRr32PEbN9SNYNs2dU+fMCGb5oKICPjwQ7VGsXOn8g46c0bNIPKrgkhJ69bKqLN2raUlMRmbNinPt2HDLC1Jxqw+vJp1R9fxWffPKF8y5ymLtdeTpsAw/d/pTFo7icGtBvNNn28eOx4UpIzWV68qZ6Tnn8/miXx81HrVtWuqOtFXX0GZMjkTPq9hZ6eUxeHDlpbEJEgJX3yhspj07GlpaYxzJfQKQxcNpUGFBrz77Lsm6VMrCk2BYPn+5by/4n1eaPwCs16d9VicxMWLSkmEhqonx2xFWt+6Be+8A8uXq7DdVaugeXPTDCAvUqJEvrFRrF6torF//12lybJW4g3x9JvVj3hDPCv+twJ7O9PYv9JVFEKI8PSO3W8ChEgpc+IHotGYne1nt/Pa3Nd4utrTLHljyWOpDM6cUXWs4+JUGqWmaSdaNs6qVWpNIiJCPXqOHVvgkualSUKCpSXIMWFhyj5fu7aaIForSUlJDFkwhD0X97DifyuoXtZ0t2ZjM4oLUspGxr4shMgf80pNvuVUyCl6/NKDyqUqs/rt1Tjap34cvD+TSExUOZuy7PYYFaXuIrNnKw0zf74KQNPA6dNQo4alpcgx77wDwcHK1GTNNc7HrhzL4j2LmdxzMi82fdGkfRszZmfG8TZnzrkajRm5HnadTj92wsHOAd93fXFxSp1vIzBQzSRiYpT3apaVxKFD0LixKmk3bpy6k2gloYiJUe6/eSXgIB0WLFDebxMmWO8q4v2Mx99t+o4R7UbwUeePTH6OdPWjlPJiRl/OTBuNxhLEJcTR69de3Iq4hf8H/lQuXTnV8cBAlQk0IkIpiSxnAp07F4YPV0Fy//2nNI7mIT4+KoNijx6WliTb7NypVhPbt1cp5K0RKSWjV4xm+ubpvPH0G/zQ/4c085TllGy5xwohjplaEI3GVEgpGb5kOHsu7mHh6wtp6pna6HDjBnToAOHhSkk0bpyFzmNjVYbUIUPg6achIEAribSYPVulWr2flzuPcfq0yu3l6amC6a0xJjLeEM+QBUOYvnk6I9uPZNars3KUz8kYxozZvdM7BJQ1izQajQmY/u905u6cyyddP+GFJqlXR8PDVTBdSIiaCGRJSVy7pp6QDxxQS02TJxeMmIis8t9/SgN/8w3Y5L1QrfPn4dlnlejr1uUwIt9MhEaG8sLMF/A/68/EbhOZ2G2iWWYS9zFmmlkOLAHSyuyV/9JCavIF64+uZ8xfY+jTpA+Tuk1KdSw2ViU2PXZMrYy0aJGFjo8cUeU/792Dv/+2fmd6S2EwwJgx4OGhrMB5jAsX1FJTfLxK12Ftmd0BjgUdo/fM3ly5c4XFQxYzoMUAs5/TmKI4CnwrpTz+6AEhRAfziaTRZI+Lty4yYPYAGlVsxILBC7BJ8TSbmKgKw/33n6pr3alTFjpet05lQy1RQuUbb9DA5LLnGyZPVstxf/6Z59KMHz6sfhf3i/hZox1+/s75DF86nGKFirFl9BZaVW2VK+c1Ni8cBaQXS5EL9Rc1mswTnxhP39/7IoRg5VsrKeJYJNXxsWPVRGD6dHj11Sx0PHOmWm6qWVOVI9VKIn127VKK4tVXVQGjPMSWLcqc4uiongUaGQ0MyH0iYyMZMn8Ig+cPpnnl5gR8GpBrSgKMKAop5XYp5ZV0jh0wn0gaTdb56cBPHLx8kPmD5+NZyjPVsd9/h++/VwXYRo3KQqdff608mzp3VkEW5cqZUuT8xaVLyvrr4QE//2xpaTKNlCrp4/PPQ8WKytOpZk1LS5WaXed30fDzhszbNY8JnSew+f3NlC2eu2bidBWFECLD1FeZaZPB9+cKIW4KIR5b3ko+LoQQPwkhzgshjgohsmJ61BQQ5myfw/rz6xnfaTw9GqZ2x/znH6UgOndWheQyhZTKcX7cOFVxbtUq608Xaknu3YMuXdTC/vr1qlZoHiAuTlWfHTFC/T527XpYE9saiEuI44+AP3j6m6dJTErEb4wfk3tNNptnkzGM2SjGCSGMVYYVwEhgVg7OPx+YASxM53gnoFry1hyYmfyq0QAQcCWA4UuH07RsU77o+UWqY8ePq2zedeuq9NCZiqqVUlVq+/lnGDpULT1pz6b0CQtTJf/On1da2doex9PhwgVVO+rgQfjoI5V1xZoctHac28GwRcM4FXKK11u9zvR+0ylWuJjF5DH2r+MPdMvg+//m5ORSym1CCE8jTXoAC6WUEtgjhCghhHCXUobk5Lya/EFMfAwDZg/A1cmVj1t9nOpJKyxMrYQ4OSlbdKZqHEsJH3yglMR776nU4HmtJGluEhYGHTuqu+1ff+WZeJK//lJhMDY21ufAFhYdxrhV4/jN/zc8XD2Y2nYqHw740NJiGY3MHpybgqRDeeBqis9Byfu0otDwwV8fcDLkJJve24T9zYcRUUlJquxDYCBs3ZqF5YSJE5VyGDFCK4mMuH5dLTcdPao8nPJABPb95H7z56t0HN7eKqDOGkhKSmLh7oWMWzWOWxG3eK/De3ze43MO7LEOc3CGk3EhxBPAj0ALVEzFbmCUlPKSmWXLNMm2kmEAbm5u+Pn5pToeGRn52L78TEEY7/ar2/ll2y+8WPNF7G/apxrz/Pme+Ph48s475zAYgsnMn8Jj4UIqz5tHSOfOnOnVSxmvrRxLXecily9T/8MPsQ8L48QXX3CnRAky9UfOITkZ7+HDJfj665rcuuXIgAFXGDgwkMBASWCgSUXMFmdCz/DTgZ84efsktUvV5rOOn1HDtQYH9hywnv9lKaXRDdgDvIpSKnbAK8DejL6X2Q3wBI6nc+x34KUUn88A7sb6a9KkiXyUrVu3PrYvP5Pfxxt4O1CWeLeEbPpFUxmXECelfDjm1aulBCkHDZIyKSmTHf78s/rSwIFSJiaaRWZzYJHrvHGjlCVKSFmmjJT79+fqqbMz3nv3pBw+XF3eatWk3L3b9HJll6A7QXLwvMFSDBXS7X03OX/nfJn4yO8vN68xcECmc1/NjPmmiJRykZTSkLwtJvcis32A15K9n1oAYVLbJwo0SUlJDJo3iMSkRLyHeeNg97DmQ2AgDBwITz6pbNCZWjlav14Zr7t3V1lgzWTRXLJkCZ6entjY2ODp6cmSJUvMch6zISV8+aWKSKtUCfbsyWbhjtxjzRpVQ2LmTHWJAwKyGI1vJsJjwpnw9wSqfVyNJXuX8P5z73PmizMMfGpgqiBRayIzfiC+QohxgDdq6akfsEEI4QIgpbyT3ZMLIZYBbYFSQoggYCJgn9zvb8AGoDNwHogGrMFuorEgv2/7Hb8zfsx+bTZVylR5sD8xUfDKK+p+tnx5JoOCAwKU60vDhrB0qdm8m5YsWcKwYcOIjo4G4PLlywxLLrw8YID50y/kmDt3lPV39WrlLvzHH1btLnzlCrz/PqxcCfXrK4N1s2aWlgpiE2KZtW0Wk9dP5lbELV5q9hJTek55LLOxVZLeVEM+XO65ZGS7mNH3c3vTS0/5d7znb5yXTm87yee+f04mPbKuNGBAoAQply7NZGeXL0vp7i5lhQpSBgebXtgUeHh4SNRDVqrNw8MjR/3mynX291d/Izs7Kb//PgvreaYno/HGxko5ebKUhQurbcoUKePjc0c2Y8QnxMvf/X+XFT6oIHkD2XZaW7nv4r5Mfddalp4ynFFIKfOAutPkdxKTEhk8fzB2NnbMHTg3VabMTZtg6dJKvPGGeuDNkJgY5aUTFaXyNZg54vrKlTQTHKS73yqIj1fpOKZMgSeegN27rXapSUq1gvjeeyqc44UXVCR+pUqWlcuQaGDZvmV8vu5zzt88T/PKzZk/eD7P1nzWrJlezYEVF/bTaB7y/b/fs/3cdhYMXkAFl4f+rjduwCuvgIdHND/+mMnlkHfeUctO69Zlo2JR1qlUqRKXL19Oc79VEhAAgwapjLmvvQYzZmQyECX3CQiA0aNVrqYaNdRDw3PPWVamBEMCi/Ys4ssNX3Lh1gXqV6iPzwgfutbvmucUxH2s03Ki0aTg7PWzfLz6Y3o37s2rLR9m9JNSpWIKD4eJE09QpIiRTu4zf74yWn/8sYoDyAWmTJlCkUeEK1KkCFOmTMmV82ea+Hj4/HPlDXD9ulrcX7DAapXEJ5+oeiJHjsBPP6n08ZZUEnEJcfyx7Q+qf1ydIQuGULxwcf4e/jeHPzlMtwbd8qySgAxmFEKNrIKU8qqxdhqNuZBS8taStyhsX5hfXv4l1T+bt7dKwzR1Knh6Rmfc2ZEj8NZbqirNpEnmE/oR7husJ0yYwJUrV6hUqRJTpkyxLkP29u0q8dHJkyof+08/gaurxcS5e1cVDEpKSt8RrUEDZbSeMMGyxYXCosP4fdvv/LD5B0LCQmhWuRkzXp5B53qd87RySIlRRSGllEKIDYD55+caTRos2buELae3MHPAzFQZM0NCVLK/Fi3U0sOOHRl0FBGhUl+7uMCyZbmev2nAgAHWpRjuc+eOysE+Z47K/Lp2rSrQZEFWrlQJ+r77zri3cp8+ls1mfiP8Bj9s/oFf/X4lPCacDrU6sPD1hbSv1T7fKIj7ZMZGcUgI8aSUcr/ZpdFoUnAr4hbvLX+PFk+0YNgzDxMVSwlvvKFs0gsWZDLZ39ixKhOcvz+UKWM+ofMKBgPMmqXWb8LC1N/n008t4vYqJfj6qlRRhQtDaKjK03X/mLXdc48FHWP65uks2buEhMQE+jTuw4edPqSJRxNLi2Y2MvMv1hwYIIS4DEShssZKKWV9s0qmKfCM9B5JWEwYf7z2R6pApNmzYcMGtTpSvXomOtqyBX77Ta1TPP20+QTOK/j7w7vvqjxNbdvCjz+qgAMLEBoK3bop/bRnjwqYjIuD06ehdWvrURJSSjad2MR3/37Hvyf/pbBDYYa0HsKo9qOoXjYzP8K8TWYURUezS6HRPILvMV+W7VvGZ90/o275hzUpb95UD79t26qlpwyJjFTTj6pVVS7pgsy5czB+vFrbqVRJJfN74YVcvxvfuqVSwD/1FJw5o4Lix42DRYuU15LBoGYXW7aoff37PyxLmtszjMjYSBbtWcSMLTM4GXIS9+LufNnrS/7X5n+4OLnkniAWJkOvJynlZaAEKuV4N6BE8j6NxiwkGBJ4f8X7VCtTjXGdxqU6NnasCn/49ddMZtv48EOV22PePDLnFpUPuXVLuQTXrg0bNypD/qlTaoE/F++6V68qe3mbNqo+taOjWgXctk1l+V23TnktnT0Lly8rEXv0UC7Q94vm5Za4F29dZPSK0VQYW4HhS4ZTyL4Q8wfPJ3BqIOM7jy9QSgIylz12JDAUWJW8a7EQYpaUMu/UO9TkKWZtn8Xp66dZ8/aaVLmctm1TNonx46FWrUx05OenNMr776t1jIJGRIRaVpo2TWnXoUNVKvWyuVtGMzERoqNVqY+bN1XpcWdndaxqVWUmWbtWTXJcXVVJ0tKllehr1sDevWqJKinJvHImJSXxz4l/mOk/k3VH12FrY0ufxn1459l3aFmlZb4zUGeFzCw9DQGaSymjAIQQX6NSjWtFoTE5tyNu88nqT3i25rN0a/CwblZcHAwbphxzJkzIREdJScodysNDRRgXJGJjVSa8L7+E27fV2s7UqZnUrqZj61ZlL69dW01o2rZVogGcOKEcrp56Clq2VPUh7s8Qd+xQdSOaN1fR1osXwzffmK8C3e2I28zdOZff/H/j0u1LlClaho86fcRbbd+ifMny5jlpHiMzikIAiSk+Jybv02hMzoTVEwiPDeen/j+leoL76iu1nu3rm0nHnBUr4NAhNQUpXNh8AlsTcXFqiW3KFAgKgg4dlJJsnvvVg3fsgM8+Uzbzw4dVvXJXV7h4UdnN3d3VzGHmTHVtPTxUvF9EhIr3K1RI6fnLl5WyqFTJtCUvpJTsPL+TWdtmseLACuIMcbSp3oaven9Fr0a9Us1kNZlTFPOAvUKIv5M/9wTmmE0iTYFl/6X9/LH9D0a2H0md8nUe7L9yRT0Q9++vyjNnSEKCWs+oWxesMXbBHOzZo2wOwcHqEX3hQouWJg0IUDf/3r2hfXs1C7xzR01qOnRQk5x791Qg+LRpKkuIg4NaptqzR9kiHB0z6dWWBUIjQ1m0ZxGzts3iVMgpihYqyhtPv8Fbbd5K9ZvTpCYzSQG/F0L4AfcXeQdLKQ+bVSpNgUNKycjlI3Er5sZn3T9LdeyTT9Tr119nsrMFC1R2uDVrcj2wziwsWKDuuMZSaVSrphTj/Pnqzmzh9fQqVVQp7Rs3wM1NLTlFRCgd1r27alOihBL57Fm1UiiEMm6vX29a8aWUbD+3nVnbZvHXwb+IM8TR4okWzBk4h35P9sPJ0XpTplsLmUoKKKU8BBwysyyaAsy6o+vYfWE3v7/6O8UKF3uwPyBAuUiOGZPJbKBxceoxtUUL5aCfVzEYVCRhTIyKd+jSRVmA07uDuroqjyYr4amnVKqo//1PTfCcnZWtIjLyYZtdu5ReGzs2tf0hUwGUmSD4bjALdi1g3q55nL95nuKFi/PG028w7Jlh1K+gw8Cygs4eq7E4SUlJTPh7AlXLVGXwUw9rU0mpPGVKloSPPspkZ/PmKT/MOXMs/lSdLTZvVnfOnj1V/EdsrLI3lCplacmyRPHiKiByzx61xNSzp/JUjohQnk+vvaY8od5+27QZQ+IN8aw/up45O+bge9yXJJlE2xpt+bTrp7zQ+AWKOBZQF+kcohWFxuIs3beUY8HHWDZ0GfZ29g/2b9yo7pvTp6tligyREr79Vq1vdOhgNnlNTng4FCumZkMbNyorcFSUCkHv0+dhHYj4eHXXtaIUJFFRqjhgzZqPB70XKqQ8ne4TE6Oyuru4qEmfKavOBVwJYP6u+SzZu4TbkbcpV6Ic4zqN4/VWr6eqhKjJHllSFEKIUkBocjUkjSbHxCbEMuHvCTTxaELfpn0f7E9IUF4v1aqpVOKZYs8elc/pk0/yxmxi2zbl9hMYqJbJxo1TIcujRil3n/h4tdB/756aJQ0ZAo0aqag0CyvC8+dViMrcuSpV1PDhaWdHkVJ5Na1YoZYOmzRRS0umUBK3Im6xZO8S5u+cz5GgIzjYOdC9QXcGthyIV10v7Gz1c7CpSPcvKYRoAUwF7gBfAIuAUoCNEOI1KaX1LIhq8iw//fcTV+5cYf7g+anyOS1YoCJzV69W3jCZYulS9Rh7P6OcNWIwqLUXR0cVHNC5s7qTjhihos7GjFGzonXr1Pvq1ZUr0LJlqmybra16Xzn3C08aDErE339XqTZsbVUGkBEjoFWrtL8jhAqqmz9flSbPKbEJsaw7uo5Fuxex4fgGDIkGnvR8kl9e/oX+zfoXuIjp3MKYyp0BfAQUB7YAnaSUe4QQNYFlgFYUmhxxN+ouX274ks71OtOu5kNXzsRE5eHUpMlDD5kMkVLdxTp2VMs41kZoqFIIS5eqVzs7Jef96VKVKmpd5okn1B31gw+gfHk1rerbV/kIb96slqHu3VPtr+ZOmZjgYLUK9scf6n358irAe9gwFQ+REX37ZtzGGPdjHr7d+y07Vu0gLCYM9+LujGo/ioFPDUyVC0xjHowpCjsp5SYAIcTnUso9AFLK06YKZRdCeAE/ArbAbCnl1EeODwKmAcHJu2ZIKWeb5OQai/PL1l8IiwljSs/Uld5WrlRLG3/9lfkVpCKBgWq5JlNh2xZg0yaVcykgQNkYgoOV0hg0SOW0CA9XswwvL1VY6T4VKqj9w4apR/l331WWYDOTmKhE/uMP8PFR7qsdO8IvvygHLFN5Jhnj7PWzLNm7hCV7l3Dh1gUK2Raiz5N9eK3lazxb81lsbfKB63MewdjlTplZJeaRYzm2UQghbIFfgOeAIGC/EMJHSnnykabLpZQjcno+jXURFRfFD//9QOd6nWlYqeGD/ffXtGvUUJ4ymaX48ePqTcqbrLUQH6+CAz7+WCmJixdVWPK6dSrseMIEFRg4Y4ayTxw9qr6XmKjyYLzwgsrPNG6cepQ3I8HBymFszhw1iSldWqXKevNNNdkxNzfCb+C9z5vFexZz4PIBhBC0q9GOT7p+Qumo0nR+rrP5hdA8hjFF0UAIEY5K11E4+T3JnwuZ4NzNgPNSyosAQghvoAfwqKLQ5ENmbZtFaGQoEzqnngFs3KgeuufNy1qsXNGzZ5UfbW7czbKKg4NaLvrzT1Xjc9s2FWhQubKaZdyPHh8xQvmUXr6swpptbeGVV1QQgpRqZmEGDAaVGmXOHKW7EhOVvfzbb1X21kzbiLJJeEw4awLWsHTfUjad2ESSTKJxpcZ89+J39G/Wn3IlygHgZ8ocHposISzlwCSE6AN4SSnfSP78Kir54IgUbQYBXwG3gLPAe2nV7xZCDAOGAbi5uTXx9vZOdTwyMhLn++kqCwDWPt74xHheXvMyFYpW4IfnfniwX0p4551G3L7tyKJFe7G3z/xvs86oUTgYDByeMcMMEuecMlu24LFoEVf79eO6lxduGzdS8tAh3DZv5sTnn3O7dWvKrV6N88WLnHv3XWQm1nZyep2Dggrj61uWf/4pS2ioIyVLxtOpUwhduoRQrlxstvvNDPGJ8ey9tpf/Av9jd/Bu4hPjcSviRofKHXiu8nN4FPd47DvW/rs2B7k55nbt2h2UUjZN86CU0iIb0Adll7j/+VWUDSJlG1fAMfn9/4AtGfXbpEkT+Shbt259bF9+xtrH+/N/P0veQG45tSXV/q1bpQQpZ8zIep9hNWpI2amTaQQ0B4GBUrq4SLlypfocEiLlp59KOXiwlEOGSFmzppQ9ekjp75/pLrNznaOipFy4UMo2bdTf2sZGyq5dpfz7bynj47PcXZZIMCTIzSc3yyHzh8ji7xSXvIEs/V5p+faSt+XOcztlYmKi0e9b++/aHOTmmIEDMp37qiUdjYOBiik+V+Ch0RoAKWVoio+zgW9yQS6NGYlLiGOq71RaV21N2xptUx379luVF+j117Per0hKUlMSS3Hjhlpa8vZWdoY+fVIf9/BQ7q6LF6u8TUeOKFvE4sXqeGhoJnOUZB0pYfdutZy3fLmKjn7iCZVkdtAgKFfOLKcFVNT9rgu78N7vzZ8H/uRmxE2cHZ3p1agXA5oPoH2t9jreIQ9gySu0H6gmhKiMUhD9gZdTNhBCuEspQ5I/dgdO5a6IGlMzf9d8gu8FM3/w/FRpxIOD1Tr5uHHZywoeXbEiRU/l8s8jNFQlNPL2VkbnpCTl4ppe4YRx45ThulMnlUr1nXdUznQpM5k7PWuEhKg8WXPnqhTtTk5Kfw0erILjzFXfQUrJoSuH8N7nzfIDy7l65yqFHQrTtV5X+j3Zj871OlPYoYCkfs8nWExRSCkNQogRwD8o99i5UsoTQojPUVMgH+BdIUR3wIAK/BtkKXk1OUdKyY///UgTjya0r9U+1bGFC9V9dvDgdL6cAZFVq+K2ZQucPq3ySZiLc+eUv6iPjyq6kJSkwsc/+gj69XtY3DkthFDFhK5eVV5PKfebiJgYlTR3wQLl3pqUpIr7jR0LL75oPAFtWpy5foZC9oXwcPVASplulTcpJUeuHmHFgRWsOLCCC7cuYG9rj1ddL6b2nkr3Bt1xLlSw7Av5CYvO+aSUG4ANj+z7NMX78cD43JZLYx78z/pzKuQU8wbNS3XDkVLd2Fq3VlG82eG6lxdVli1Td0QfHxNJjErKt307bNigXFzPnVP7GzRQs4OePVVajazc7FMqCRNwf2lp/ny1tBQerk4xbpwKuahRI2v93Ym68yDCedm+ZVQrUw0PV490lUTQnSDaf9+eszfOYmtjS/ua7RnXaRy9G/fWkdL5BL04qMk1ft36KyWLlKTfk/1S7d+/Xy2NjBmT/b4T7qeYHT9eVXV7993sRWiHhcGBA0o5+Pmp/FFxcSo1SLt2MHKkSnfq8bhXTm4THFyIzz5TZo7z56FIEbW0NHCgSsaX1aUl32O+DF04FK+6Xgx6ahBNPJoQfC+Y8Z3Gk5iUmG6AW7kS5WhUqRGjnx9N70a9KVU0b2W61WSMVhSaXCHoThCrDq/ivQ7vPbY+vXixChF48cUcnuS999Sj9SefqDxKgwerVN1Vqqi7aEpiY9US0KlTarnq+HGlsU6fVsdtbNRM4e23VRBfu3aP92EB7t5VCfYWLIDdu1sghFIKH32klERWlpb8zvhxM/wmfZ/sS7whnp3nd7L67dUUti/Mt5u+pW75ukgpcbQ3Hr9hY2OD9zBvo200eRutKDS5wm/+v5EkkxjeNnUqWINBLZd066ZqGOQIR0e1QL9/vwpcmzlTvYK6g5YpoxTEnTtqMT8l5cqpPEqvvKJSmzZrZgKBTENsrFr5WrJEBcTFx0OdOjBs2AU+/rhKlleyboTfoPevvXF2dKZ8yfKcv3mecZ3Gce7mOQyJBm7H3aZZ5WasOrSKeuXrERoZyvhV42lUqRG9G/fGrZibUXuFJv+hFYXG7MQmxDJr+yy61e9G5dKps576+6tkqi+9ZMITPvmkcveZNk2Fel+/rtxXb95US0guLmpzd1dFnGvUyGTBi9wjKUn9bZYsUTmvwsKUnnvzTbW01KgR+PtfpWLFjGstbDi2gaV7l1KvfD3ebPMmfx/6G6+6XnzS9RMSkxJ5ff7rLNi9gG9f/Jbf/H9j6d6lfN/3exztHbkZcZNtZ7fRskpLXJxcmLxuMj+//LNWEgUMrSg0ZmfZvmXcirjFyPYjHzu2YoVy2+zUyQwnLltWBQrkEaRU4RXLlql8gUFBKntH794qy8ezz2Y9Gd/mk5uZsWUGA58ayJ6Le1i8dzG7L+ymkL3KwmNrY4ujnSO/bP2Fwa0GM67TOKb0Ukka31r8Fp90/YR3lr3DLy//gqOdI/su7eNWxC1KFy1t6uFrrBitKDRmZ6bfTOqUq5MqlTioG+O6dUpJZCd2Ir9w/vxD5XD6tFIGzz+vzCw9euTMNLLrwi6qlK5Cvyf70fKJlszaNgsnRyd2XdjFnwf+JDo+mpJFSnLx1kWu3btGuRLlHhiuhRDY2tjSq1Evvt74NasOreJ/z/xPK4kCiFYUGrNy9vpZ9gfu57sXv3tsueLMGbh2TaWvLmgEB6tA7qVLlUkF4JlnlFNVnz6mK5FdyaUSl0Mvq/euldh5YSdVSldhSs8p+J/159LtS3z9wtcE3wtm5/mdvNj0RWxtbNl3aR+GRANuxdzo/2R/nvR8ktdavkajSo1MI5gmT6EVhcasLD+wHCHEYy6xAHv3qtf0qqPlN27eVPaG5cuV962U0LixMqX0769KT5iabg26sXD3Qt5e8jZBd4PwdPWkiEMRyhQtw3d9v0vVNmUBIA9XDxpXagyAna0dNcpmMRhDk6/QikJjVpbvX07rqq0pX7L8Y8eOHlVLTtWrW0CwXOLu3YdZPv77Txmpa9eGzz5Tld+yGgyXVVydXVnyxhL+O/0fZYqWobpbdcatGodnKU/uRd/jePBxFu9ZTNFCRanlXuvB99yKufFm2zfNK5wmz6AVhcZsnLx2khPXTvBT/5/SPH7ihHI6ykrdibzAnTuq1veff6rqpQaDCuUYP17NHIxl+cgO4THhrD+6njrl61C/Qv3HjruXcOeVFq8AEBYdRvDdYFycXFh7ZC1zdsyhUaVGvPPsO6YVSpOv0IpCYzb+OvgXQgheaPJCmsevXjVvWqbc5O5dpRxWrHioHDw9VQzgiy+qEA1TepSGRoay4cIGph2dxuZTm4k3xPNBxw/4pk/aCZa3nt7K3J1z2XdpHx96fYiDnQPdGnRL99poNCnRikJjNv46+BetqrR6UKHsUYKDoX37NA/lCW7eVMph5UrYsuWhcnj/faUcmjQxrXK4EnqF1QGrWX14NdvObSMxKREPVw9GtBtBr0a9aFmlZbrfjTfE07BiQ6b0nEIlV5XO3MHOzKXrNPkGrSg0ZuH8zfMcCz7GD/1+SLdNbKxVZMXIEkFBsGqVUg73k8dWqWK+mcOZ62dYdWgVKw+t5ODlgwDUdq/Nh14f4iE9GNpraKaC3zrW7UjHugXQvUxjErSi0JiFvw//DUCvRr3SbZOYaNqbqjmQEk6eVJlB1qyBffvU/rp1VUqp3r1VCQpTjUNKycHLB1l9eDV/H/6bkyGqhHzzys2Z2nsqvRr1onpZZf338/PTEdKaXEErCo1Z8AnwoVGlRg+WOdKiZElV+8faMBhg166HyuHCBbX/ySdVOYkXXjCtp5Yh0cC2s9tYdXgVawLWEHQ3CFsbW56u9jRvtnmTXo16UcHFDL6zGk0m0YpCY3JCI0PZdWEXH3f52Gi7KlUe3oQtzb178M8/quTEhg1KgTk4KBvKmDHQvbtpS4ZGx0Wz6eQm1gSsweeID3ei7lDYoTAda3dkSs8pdKnfBVdnV9OdUKPJAVpRaEzOjvM7SJJJPF/neaPtGjZUufuio3PfVnE/r5Kvr9p27VJLYS4u0KWLUgwdO2a9IpwxbkfcZu3Rtaw+vJpNJzcRmxBLiSIl6Fq/K70b9aZjnY4UccxjRhtNgUArCo3JOXj5IDbChkYVjad7ePFF+O03FW8wcKD55bpxQ7mubtqktuvX1f5GjeDDD5WCaN7ctHEdl0MvsyZgDX8f/pttZ7eRJJOo6FKRoU8PpUfDHjxT7Rns7exNd0KNxgxoRaExOUeuHqFm2ZoZPh23aaNSWLz5pir90LOnaeW4fv1hoTo/P2WUBjVreO45NWPw8lLZxk2FlJKAqwGsPryaNQFrOBJ0BFCeSuM7jadX4140rtRYG6E1eQqLKgohhBfwI2ALzJZSTn3kuCOwEGgChAL9pJSBuS2nJmucCjlFg4oNMmxna6vsAl27KgPxjBlKaWTnHhofb8OBA3DwoCpyt2PHQ/uHszM8/bSqH92hg1ryMuWsIcGQwI7zOx7EOFy5cwUbYcNTVZ5iWp9pdG/Q/YGnkkaTF7GYohBC2AK/AM8BQcB+IYSPlPJkimZDgLtSyqpCiP7A18Dj2eU0VoMh0cDF2xfp27RvptqXKqVyIPXrB8OHq6WoNm1UosBmzVTZawcHtUVEqCC3mzdVVPfZs2o7fRpOnGhNYqLqs3Rp9f233oLWrdWsxd7EqzvhMeH4HvfFJ8CHDcc3cC/6HoXsC/F87eeZ1H0SXet31em4NfkGS84omgHnpZQXAYQQ3kAPIKWi6AFMSn7/FzBDCCGklDI3BdVknhvhN0hMSqSiS+brczo5qQjnP/5Q2VXnzIGff874e7a2ULmyclWtV+8qPXt60KSJio42x8pO0J0gfI74sCZgDVvPbCUhMYFSzqXo2bAn3Rp0o2Odjjg5Opn+xBqNhRGWuucKIfoAXlLKN5I/vwo0l1KOSNHmeHKboOTPF5Lb3H6kr2HAMAA3N7cm3t6pC71HRkbi7OxszuFYFZYc75nQM7y58U0mt5lMqwrZyx9uMAguXHDm/Hln4uJsMBgE8fE2FCmSSIkS8ZQokUCpUvG4u8dgb69+v+YYc5JM4tydc+wK3sXuoN2cu3sOgIpFK/JUhadoVaEVtUvVxtbGMlkN9e86/5ObY27Xrt1BKWXTtI7lC2O2lHIWMAugadOmsm3btqmO+/n58ei+/Iwlx2tz1gY2QosmLWhbK/sydOiQtfamGnNMfAxbTm9h7ZG1rD26lmv3rj2wNwxpN4QeDXtQ0906Mhnq33X+x1rGbElFEQykXJ+okLwvrTZBQgg7oDjKqK2xUhISEwCwt807Lp9Bd4LYcHwD64+uZ/OpzUTHR+Ps6IxXXS+61e9G53qdKVXURCXnNJo8iCUVxX6gmhCiMkoh9AdefqSNDzAQ2A30AbZo+4R1U7SQilALiwmzsCTpk5iUyL5L+1h/dD3rjq574MLq4erB4FaD6d6gO22qt8HR3tHCkmo01oHFFIWU0iCEGAH8g3KPnSulPCGE+Bw4IKX0AeYAi4QQ54E7KGWisWLKFVd5LoLvPjo5tCy3I26z6eQmfI/7svH4Rm5H3sbWxpZWVVvx9Qtf06VeF2qXq63jGzSaNLCojUJKuQHY8Mi+T1O8jwVezG25NNnHvYQ7xQsX5+CVgxaVw5BoYO+lvfxz4h82Ht/IgcsHkFLi6uxKp7qd6Fq/K8/Xfp6STiUtKqdGkxfIF8ZsjfVga2NLh1od8D3mi5Qy157QpZRcuHmB/07/x6YTm9h8ajNhMWHYCBtaPNGCz7p/Rsc6HWni0cRiXkoaTV5FKwqNyenRsAcrD61kpt9MhrcbbrbzXA69jN8ZP/zO+OEb4MuNpTcAqFCyAn2a9MGrrhfta7bXswaNJodoRaExOf2f7M/y/csZsWwERQsV5dWWr+a4z6SkJE6FnGLnhZ3sOLeD7ee2ExgaCICLkwt1XesysddEnq35LNXdqmtbg0ZjQrSi0Jgcezt7/nrrL7r81IXX5r7G5lOb+aTrJ1QpXSVTN/DouGjO3DjDiWsnOHT5EIeuHOLw1cOEx4QDUKZoGVpVbcX7z71P2xptqVOuDtu2bbMKf3ONJj+iFYXGLBSyL8S6d9Yxef1kpv0zjYW7F1KmaBlaVmlJ/Qr1cbB1wNbGFiEEtyNvE3IvhJCwEC7dvsTlO5e57wVdyL4QDSs2ZEDzATSv3JxWVVtlWuFoNBrToBWFxmwUdijMlF5TGNJ6CJtObGL3xd3subiHNQFrUrUrZF8I9+LuuBd3p2WVlrze+nVqlq1Jbffa1ChbAztb/TPVaCyJ/g/UmJ0nSj/Bm23f5M22bwLK3pCYlEiiTCQpKYnCDoX1DEGjsWK0otDkOjY2NtjY2GBP3knzodEUZGwsLYBGo9ForButKDQajUZjFK0oNBqNRmMUrSg0Go1GYxStKDQajUZjFK0oNBqNRmMUrSg0Go1GYxStKDQajUZjFK0oNBqNRmMUrSg0Go1GYxStKDQajUZjFK0oNBqNRmMUiygKIYSLEOJfIcS55Nc0a1UKIRKFEAHJm09uy6nRaDQay80oxgH/SSmrAf8lf06LGCllw+Ste+6Jp9FoNJr7WEpR9AAWJL9fAPS0kBwajUajyQBxv+Rkrp5UiHtSyhLJ7wVw9/7nR9oZgADAAEyVUq5Op79hwDAANze3Jt7e3qmOR0ZG4uzsbLoBWDkFbbygx1wQKGjjhdwdc7t27Q5KKZumdcxshYuEEJuBsmkcmpDyg5RSCiHS01YeUspgIcQTwBYhxDEp5YVHG0kpZwGzAJo2bSrbtm2b6rifnx+P7svPFLTxgh5zQaCgjResZ8xmUxRSyg7pHRNC3BBCuEspQ4QQ7sDNdPoITn69KITwAxoBjykKjUaj0ZgPS9kofICBye8HAmsebSCEKCmEcEx+XwpoBZzMNQk1Go1GA1hOUUwFnhNCnAM6JH9GCNFUCDE7uU0t4IAQ4giwFWWj0IpCo9FochmzLT0ZQ0oZCrRPY/8B4I3k97uAerksmkaj0WgeQUdmazQajcYoWlFoNBqNxihaUWg0Go3GKFpRaDQajcYoWlFoNBqNxihaUWg0Go3GKFpRaDQajcYoWlFoNBqNxihaUWg0Go3GKFpRaDQajcYoWlFoNBqNxihaUWg0Go3GKFpRaDQajcYoWlFoNBqNxihaUWg0Go3GKFpRaDQajcYoWlFoNBqNxihaUWg0Go3GKFpRaDQajcYoFlEUQogXhRAnhBBJQoimRtp5CSHOCCHOCyHG5aaMGo1Go1FYakZxHOgNbEuvgRDCFvgF6ATUBl4SQtTOHfE0Go1Gcx87S5xUSnkKQAhhrFkz4LyU8mJyW2+gB3DS7AJqNBqN5gEWURSZpDxwNcXnIKB5Wg2FEMOAYQBubm74+fmlOh4ZGfnYvvxMQRsv6DEXBAraeMF6xmw2RSGE2AyUTePQBCnlGlOeS0o5C5iVfN5b7dq1u/xIk1LAbVOe08opaOMFPeaCQEEbL+TumD3SO2A2RSGl7JDDLoKBiik+V0jel9F5Sz+6TwhxQEqZrtE8v1HQxgt6zAWBgjZesJ4xW7N77H6gmhCishDCAegP+FhYJo1GoylwWMo9tpcQIghoCawXQvyTvL+cEGIDgJTSAIwA/gFOASuklCcsIa9Go9EUZCzl9fQ38Hca+68BnVN83gBsMMEpZ5mgj7xEQRsv6DEXBAraeMFKxiyklJaWQaPRaDRWjDXbKDQajUZjBWhFodFoNBqj5BtFkVFeKCHEICHELSFEQPL2hiXkNCVCiLlCiJtCiOPpHBdCiJ+S/yZHhRCNc1tGU5KJ8bYVQoSluMaf5raMpkQIUVEIsVUIcTI5N9rINNrkt2ucmTHnt+tcSAixTwhxJHnMn6XRxlEIsTz5Ou8VQnjmqpBSyjy/AbbABeAJwAE4AtR+pM0gYIalZTXxuJ8BGgPH0zneGfAFBNAC2Gtpmc083rbAOkvLacLxugONk98XBc6m8bvOb9c4M2POb9dZAM7J7+2BvUCLR9oMB35Lft8fWJ6bMuaXGcWDvFBSynjgfl6ofI2Uchtwx0iTHsBCqdgDlBBCuOeOdKYnE+PNV0gpQ6SUh5LfR6DcxMs/0iy/XePMjDlfkXztIpM/2idvj3oZ9QAWJL//C2gvMkiWZ0ryi6JIKy9UWj+uF5Kn538JISqmcTy/kdm/S36iZfIU3lcIUcfSwpiK5KWGRqinzZTk22tsZMyQz66zEMJWCBEA3AT+lVKme52lijELA1xzS778oigyw1rAU0pZH/iXh9pZk384BHhIKRsAPwOrLSuOaRBCOAMrgVFSynBLy5MbZDDmfHedpZSJUsqGqFRFzYQQdS0sUiryi6LIMC+UlDJUShmX/HE20CSXZLMk2cqXlVeRUobfn8JLFaxpL4QoZWGxcoQQwh51w1wipVyVRpN8d40zGnN+vM73kVLeA7YCXo8cenCdhRB2QHEgNLfkyi+KIsO8UI+s23ZHrX3md3yA15I9Y1oAYVLKEEsLZS6EEGXvr9sKIZqhft+59s9kapLHMgc4JaX8Pp1m+eoaZ2bM+fA6lxZClEh+Xxh4Djj9SDMfYGDy+z7AFpls2c4NrLkeRaaRUhqEEPfzQtkCc6WUJ4QQnwMHpJQ+wLtCiO6AAWUQHWQxgU2EEGIZygOkVHLurIkoQxhSyt9Q6U86A+eBaGCwZSQ1DZkYbx/gLSGEAYgB+ufmP5MZaAW8ChxLXr8G+AioBPnzGpO5Mee36+wOLBCqqqcNKq/dukfuX3OARUKI86j7V//cFFCn8NBoNBqNUfLL0pNGo9FozIRWFBqNRqMxilYUGo1GozGKVhQajUajMYpWFBqNRqMxilYUGo0RhBA9hRBSCFEzl843XwhxSQjxZvLnNLOGCiGeTs6wmmYmXY3GlGhFodEY5yVgR/JrtkgOhsvK/9oHyfECAEOAu1LKqsB04GsAKeV2UpQN1mjMiVYUGk06JOcbao26Wacb4CSEeF8IcTx5G5W8z1Oo+igLgeOkTrOBECJQCPGNEOJYci2Cqul0b9GsoRoNaEWh0RijB7BRSnkWCBVCPJYfLHnfYKA5qh7EUCFEo+TD1YBfpZR1pJSX0+g/TEpZD5gB/JCODBbNGqrRgFYUGo0xXkLVNiH5Na3lp9bA31LKqOREdauAp5OPXU6uEZEey1K8tjSBvBqNWcgXuZ40GlMjhHABngXqCSEkKoeYFEJ8kIW8QlEZHJfpvE/J/ayhQZbIGqrRgJ5RaDTp0QdYJKX0kFJ6SikrApd4OFu4z3agpxCiiBDCCeiVvC8z9EvxujudNhbNGqrRgJ5RaDTp8RLJHkYpWJm8f9v9HVLKQ0KI+cC+5F2zpZSH77uxZkBJIcRRII70vaosmjVUowGdPVajsQhCiECgqZTy9iP75wPrpJR/ZaIPz+S2VlUNTZP/0EtPGo11EQZ8cT/gLj2EEE+jyvveNtZOozEFekah0Wg0GqPoGYVGo9FojKIVhUaj0WiMohWFRqPRaIyiFYVGo9FojKIVhUaj0WiM8n/sbxmoBoh4RAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Now a scan of the chisq\n", "# We calculate the chisq for various values of p[0] and p[1]\n", "# Note: temp3[0][0] and temp3[1][1] are the diagonal\n", "# elements of the covariance matrix, ie, the square of the errors\n", "# on p[0] and [p1],\n", "# We will scan out to 5 sigma\n", "ax2 = plt.subplot(111)\n", "p0min = p[0] - 5*np.sqrt(temp3[0][0])\n", "p0max = p[0] + 5*np.sqrt(temp3[0][0])\n", "p1min = p[1] - 5*np.sqrt(temp3[1][1])\n", "p1max = p[1] + 5*np.sqrt(temp3[1][1])\n", "nscan = 100\n", "p0 = np.linspace(p0min, p0max, nscan)\n", "p1 = np.linspace(p1min, p1max, nscan)\n", "ax2.set_xlim(p0min,p0max)\n", "ax2.set_ylim(p1min,p1max)\n", "z = np.zeros( shape=(nscan,nscan) ) # an emptyt array for now\n", "for i0 in range(0, nscan):\n", " this_p0 = p0[i0]\n", " for i1 in range(0, nscan):\n", " this_p1 = p1[i1]\n", " yy = f(x, [this_p0, this_p1])\n", " z[i0][i1] = getChisq(y,yy,e2) - chisq\n", "\n", "# \"z\" is a map of chisq-chisq_at_minimum\n", "# Correspond to the 68% 95% 99% coverage for 2 parameters\n", "# (PDG Table 38.2)\n", "CS = ax2.contour(p0, p1, z, [2.30, 5.99, 9.21], colors=['red','blue','darkgreen'])\n", "fmt = {}\n", "strs = [ '68%', '95%', '99%' ]\n", "for l,s in zip( CS.levels, strs ):\n", " fmt[l] = s\n", "ax2.clabel(CS, inline=True, fmt=fmt, fontsize=10)\n", " \n", "# put a point at the best fit\n", "ax2.plot(p[0], p[1], 'ko')\n", "\n", "# label the axes\n", "ax2.set_xlabel('A or p[0]')\n", "ax2.set_ylabel('B or p[1]')\n", "\n", "# The resut\n", "ax2.text(1-0.96, 0.96, txtRes,\n", " transform=ax2.transAxes,\n", " bbox=props, fontsize='large',\n", " horizontalalignment='left',\n", " verticalalignment='top' )\n", "\n", "# put a grid\n", "ax2.grid()" ] }, { "cell_type": "code", "execution_count": null, "id": "b6e3d0c5", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.5" } }, "nbformat": 4, "nbformat_minor": 5 }