Blackbody_PowerVsDistance.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Blackbody_PowerVsDistance",
      "version": "0.3.2",
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "[View in Colaboratory](https://colab.research.google.com/github/jrstevenjlab/ColaboratoryAtomicLab/blob/master/Blackbody_PowerVsDistance.ipynb)"
      ]
    },
    {
      "metadata": {
        "id": "4HMW7dTiCIIg",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "The goal of the excercise is to fit the function \n",
        "\n",
        "$$S_{model}(x) = B + C (x - x_0)^{n}$$\n",
        "\n",
        "to data for power as a function of distance that you recorded in lab.  The parameters of the function are:\n",
        "\n",
        "*   B: a background term, the value the data would have at x=infinity,\n",
        "*   C: a scale factor which multiplies $(x-x_0)^n$, which translates the units\n",
        " from $distance^n$ to radiated power.\n",
        "*   $x_0$: maybe the filament wasn't exactly at x=0? If so, this term allows\n",
        " the fit to try and correct for that.\n",
        "*   n: the power in the power law, which we expect to be -2\n",
        "\n",
        "The function given in the lab manual is a limiting form of this more general function, where B=0 and $x_0$=0.  Here we'll treat those as parameters of the model to be determined by the fit.\n",
        "\n",
        "By \"fitting\" the data we mean that the parameters of the model are allowed to vary until the curve defined by this function best match the data you've measured.  This is done by minimizing the $\\chi^2$ which is a metric for comparing the model function to the observed data.  The $\\chi^2$ is defined as\n",
        "\n",
        "$$\\chi^{2} = \\sum_{i=0} \\frac{(S_i - S_{model}(x_i))^2}{\\sigma_{S_{i}}^2}$$\n",
        "\n",
        "where the sum of $i$ is over your observed data points\n",
        "*  $S_i$ are the individual data points for the power (y-axis value), \n",
        "*  $\\sigma_{S_{i}}$ are the errors on the individual data points,\n",
        "*  and $S_{model}(x_i)$ is the value of your model function at the position $x_i$\n",
        "\n",
        "The steps below describe the process that you should follow:\n",
        "\n",
        "1.   Enter your data for power and position in the ordered lists below, replacing the example data already there.\n",
        "2.   Run the fit code by pressing the \"Play\" button, which will cause a few things to happen:\n",
        "> * Your data will be plotted as blue points in the figure below the code \n",
        "> * The `curve_fit` command in the code will find the \"best fit\" values for the paremeters of the model which gives the smallest value for the $\\chi^2$ possible (ie. it \"minimizes\" the $\\chi^2$).\n",
        "> * The function S(x) with these \"best fit\" parameters will be plotted in red overlayed with your data points\n",
        "> * The values for the parameters are output as text below the figure, along with the uncertainties on those parameters.\n",
        "\n",
        "Finally, we need evaluate the \"quality\" of our fit by looking at the $\\chi^{2}/DOF$, where DOF is the **Degrees Of Freedom** of the fit, defined as:\n",
        "\n",
        "$$DOF = # data points - # parameters$$\n",
        "\n",
        "In our case the # parameters = 4 (B, C, $x_0$, and n). and the # data points is given by the number of points you observed in lab.   In general, a \"good\" fit will have a $\\chi^{2}/DOF$ close to 1.  \n",
        "*  If $\\chi^{2}/DOF$ >> 1, then the errors were likely **underestimated** meaning the model function does not describe the data.\n",
        "*  If $\\chi^{2}/DOF$ << 1, then the errors were likely **overestimated** meaning the model describes the data better than it should (statistically).\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "metadata": {
        "id": "tmFu5KTfCFwO",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 506
        },
        "cellView": "code",
        "outputId": "035ce5db-5ab4-408b-fbf6-fb2acc95b0d0"
      },
      "cell_type": "code",
      "source": [
        "# import modules needed for data analysis\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.optimize import curve_fit\n",
        "\n",
        "# function for power vs distance with paramters B, C, x0 and n\n",
        "def Smodel(x, B, C, x0, n):\n",
        "  return B + C * np.power(x-x0, n)\n",
        "\n",
        "##########################\n",
        "# INSERT YOUR DATA HERE! #\n",
        "##########################\n",
        "\n",
        "# enter as a list here for radiation (S) and distance (x)\n",
        "# x = distance (cm)\n",
        "x = [100,90,80,70,60,55,50,45,42,40,38,35,32,30,28,27,26,25,24,22,20,19,18,17,16.5,16,15.5,15,14,13,12,11,10,9,5]\n",
        "\n",
        "# S = converted sensor reading (W/m^2)\n",
        "S = [1.82,2.39,2.73,3.52,4.43,5.23,5.91,6.82,7.84,8.52,9.66,10.91,12.61,13.98,16.14,17.61,19.55,21.14,22.16,26.02,30.45,35,37.39,40.45,43.41,46.93,49.32,53.18,59.09,66.36,79.32,92.73,109.2,134.2,399.43]\n",
        "\n",
        "# uncertainy (error) on converted sensor reading: YOU NEED TO DETERMINE THIS!!! \n",
        "# -choose between one of the options for uncertainties below\n",
        "\n",
        "# fixed uncertainty on power for all datapoints in W/m^2\n",
        "#Serr_value = 2.0 \n",
        "#Serr = np.repeat(Serr_value, len(S)) \n",
        "\n",
        "# scale uncertainty which is a fraction of the power for each datapoint in W/m^2\n",
        "Serr = 0.1*np.array(S) # eg. 10% error\n",
        "\n",
        "# plot the data on a new figure\n",
        "plt.figure() # create new figure for plotting\n",
        "plt.errorbar(x, S, Serr, fmt='o', label=\"data\") # plot S vs x, where Serr is the uncertainty on your power values\n",
        "plt.axis([0, 100, 0, 500])  # set axes [xmin, xmax, Smin, Smax]\n",
        "plt.xlabel('distance (cm)') # set x-axis label\n",
        "plt.ylabel('Power/Area (W/$m^2$)') # set y-axis label \n",
        "\n",
        "# set boundary conditions for fit parameters\n",
        "initialParameters = [0,10000, 0, -2] # start with initial guess for fit parameters [B, C, x0, n]\n",
        "# upper and lower bounds on parameters\n",
        "# 0 < B < inf (background radiation must be positive)\n",
        "# 0 < C < inf (1/x^2 term in fit function must be positive)\n",
        "# 0 < x0 < inf (position offset must be positive to keep return value real)\n",
        "# -10 < n < 0  (power must be negative as S decreases with x)\n",
        "boundParameters=([0,0,0,-10],[np.inf,np.inf,np.inf,0]) \n",
        "\n",
        "# fit the data\n",
        "finalParameters, finalParameterErrors = curve_fit(Smodel, x, S, initialParameters, Serr, True) #\n",
        "\n",
        "# show red curve for fit results and plot on top of the data\n",
        "plt.plot(x, Smodel(x, *finalParameters), 'r--', label='fit: B=%5.3f, C=%5.3f, x0=%5.3f, n=%5.3f' % tuple(finalParameters))\n",
        "plt.legend()\n",
        "plt.show()\n",
        "\n",
        "# output fit parameter values and errors\n",
        "print(\"Model parameters:\")\n",
        "print(\"B   =\", finalParameters[0], \" +/- \",np.sqrt(finalParameterErrors[0,0]))\n",
        "print(\"C   =\", finalParameters[1], \" +/- \",np.sqrt(finalParameterErrors[1,1]))\n",
        "print(\"x_0 =\", finalParameters[2], \" +/- \",np.sqrt(finalParameterErrors[2,2]))\n",
        "print(\"n   =\", finalParameters[3], \" +/- \",np.sqrt(finalParameterErrors[3,3]))\n",
        "\n",
        "# compare model function with obtained parameters to data to determine fit quality\n",
        "Smodel_i = Smodel(x, *finalParameters)\n",
        "deltaS = S - Smodel_i # calculate difference between data and model\n",
        "chisq = np.sum((deltaS/Serr)**2) # calculate chi^2 by summing over difference for all data points\n",
        "dof = len(S) - len(finalParameters) # Degrees of Freedom = # data points - # parameters\n",
        "print(\"\\nFit quality:\")\n",
        "print(\"chisq =\",chisq,\"dof =\",dof)\n"
      ],
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAFcCAYAAABV673oAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xlc1NX+x/HXDIuIgAiCCq5Zt3JJ\nM/ctTa3UFm9pi1m5ZItW3tJcWvWXmam5r5Wm17L0YovleivNUsSS8rpVLrkBKqDIDsPM/P5AJolB\nZZxhGHg/H48eOt/5Lp85kvP2nO/3HIPVarUiIiIiFZbR3QWIiIiIeykMiIiIVHAKAyIiIhWcwoCI\niEgFpzAgIiJSwSkMiIiIVHDe7rpwTEwMI0aM4LrrrgPgH//4B0888QSjR4/GbDYTFhbG1KlT8fX1\nZc2aNSxbtgyj0cgDDzxAv3793FW2iIhIueO2MADQunVrZs+ebXs9btw4+vfvT8+ePZk+fTpRUVH0\n6dOHefPmERUVhY+PD3379qVHjx4EBwe7sXIREZHyo0wNE8TExNCtWzcAunbtSnR0NLt376Zp06YE\nBgbi5+dHixYtiI2NdXOlIiIi5YdbewYOHTrE008/zfnz53n22WfJysrC19cXgNDQUBITE0lKSiIk\nJMR2TEhICImJie4qWUREpNxxWxioX78+zz77LD179uTEiRM89thjmM1m2/vFzZJ8pbMnW61WDAaD\nU2oVEREpz9wWBmrUqEGvXr0AqFu3LtWrV2fPnj1kZ2fj5+fH6dOnCQ8PJzw8nKSkJNtxZ86coXnz\n5pc9v8FgIDExzWX1C4SFBaqNXUxtXDrUzq6nNna9sLBAh4912z0Da9asYfHixQAkJiaSnJzMfffd\nx8aNGwHYtGkTnTp1olmzZuzZs4fU1FQyMjKIjY2lZcuW7ipbRESk3HFbz8Btt93GqFGj+PbbbzGZ\nTIwfP54bb7yRMWPGsHLlSiIiIujTpw8+Pj6MHDmSIUOGYDAYGD58OIGBjqcfERERKcxQnpcwVpeU\na6nbz/XUxqVD7ex6amPX88hhAhERESkbFAZEREQqOIUBERGRCk5hQEREpIJTGBAREang3DodsYhU\nHHl5eTzzzBDq1atPixYtqVIlgFtv7crmzd/QtWv3Yo97663x/P77AYKCqpKbm8t11/2DkSPHYjQW\n/2+ZFSv+zebN3wAGBg8eSrt2He3u98UXq/noo6V8//0WLBYL06dP4fDhg+Tl5XHvvf/krrv62PY9\ncuQQgwcP4JNPPqNWrYhir33gwD7mz59Nbm4uJpOJjh07M2jQ0GJnRM3MzGTixDdIS0vFZDIxaNBQ\n2rRpx8GDf/Duu5MxGKBhw+sYNWpcsTUeP36MqVMnAfmzr44Z8yp16tS9bJsUd9yaNZ/z9ddf4uVl\npGHDfzBy5JgSzehqr/Yr9dNPMbz33jyMRi/atevAwIFPFHo/PT2dCRNeIT09ncqV/Rk/fiJBQVWL\nPe7IkUOMHTuSBx/sz/33P3jFdZTUL7/s4rXXxjJu3Ot06NCpyPuLFs3j1193YbFY6dy5C4888jgA\ns2e/y6+//oKvry+vv/4mERGRtmMWLpzL3r3/Y+7c91xWdwH1DIhIqUhKSsJkMvHqqxPo1etubr21\nKyaTiZUrV1z22Keeepa5c9/jvfeWcvLkCfbv31vsvvHxcXzzzSbmz1/MlCkzmTNnRqGpzgucO3eW\nrVs3217v2fM/vL29WbBgMbNmLWDhwnlYLBYg/4ty7txZ1K5d55J1ZmSkM2HCa7zwwkssWvQh7723\nlIMH/+Crr74o9pj167+ibt16zJmziIkT32HWrGlA/pfEiBEjWbBgCenp6URHbyu2xi++iGLIkKeY\nM2cRvXvfw4oVy6+oTewdl52dzbffbmL+/A9YsGAJx48fZe/e/13yc/+dvdqv1KxZ05g4cQoLFixm\n584d/PnnkULvr1q1gptvvoUFCxZz661d+eijZcUel5WVxYwZU7nlltYlqr+k4uJOsnLlxzRt2szu\n+0eOHOKXX3axYMESFixYzLp1X5GcnER09I/Ex8exZMlHPProIHbu3GE75s8/j7B7d+ktyqeeAZEK\nKuSWJna3Zw57nuwhTwIQOGwoPjHRRfYx3dKStPeWAuC3fCn+M6dxdlfxX9AAc+a8S1zcSSZNmkCN\nGjUJDg7m6NGjHD58iGnTJjNo0BMsXryI0aNfKfYcubm5ZGZmUq1aCL/9doC5c2cUer9jx84EBATS\ntm17fHx8qFatGjVr1uLo0T9p2PDaQvvOnz+bIUOe5o038v/V2qxZc5o1y5/q/Ny5swQFBdl6H9au\nXUPLlq3Yvv3HS37G//53A50738o11+Rfy9vbm9dem0ClSn5s3/4jK1b8u9D+99xzH1WrBnPo0CEA\nUlNTqVo1GJPJREJCPDfe2BiADh068fPPO3nuuRfs1vj88yNt5zx9+hTh4eGFrhMb+7PdNrF3nJ+f\nH7NmLQAgOzub9PR0QkJC7X7e1NTzPPfcUyxYsASz2cywYUOYPXuh3dqDgvxYsGBRkc9/++132l7H\nxZ0kMDCIGjVqAtCuXQd27dpJgwbX2PbZtesnxo17/cK5OzN69L+KPa5Pn75MmzbLFhguJTb2Zz77\nbBUGg5Fjx/6kS5duDB78JK++OoaUlHO2/Xx8fJgxY16hY0NDq/PWW1OZPPlNu+euUiWA3NwccnNz\nsVgsGAxG/Pz82LbtB3r06Glrp4vNnTuToUOHsWSJ63sFQGFARErJs8++wKuvjuHll99g8eL8L4X+\n/R9l//69jBo1FqDYILBo0Vw++WQ58fFxdO7chcjI2gB2u0///e8lBAdXs72uVq0ayclJhcJAbOzP\nVKpUicaNiwaiV18dw549v/Laa/l/sZ8/n8KGDWuZOXP+ZcPAsWPHbF+CBfz9qwDQvn1H2re3P1yx\nbt3XPPhgH9LS0pgyZSbnz6cUmmm1WrUQkpP/WqPl7zUCHDz4OxMnvkGlSn99mRc4eza52DYp7rjl\ny5cSFfUJ/fo9bGvvvwsKqsqDDz7CRx8tJScnh0cfHYzZbLZbe5cuXWjc+JZi2664OuPi4grtk5z8\n1z4Fn6O447y9vfH2vvKvuf3797FixWosFgv9+t3N4MFPMnHiO5c9zs/P75Lv16hRk65du9O3791Y\nLGYGDnyCKlUCSEhIwN//AGvWfEalSpV48cUx1KxZi3XrvqJ58xaXHI5yNg0TiFRQZ3fttftfQa8A\nQNr89+3uU9ArAJD96MDL9gpcrYJhgqior8jJyeXrr4vvdv+7v8+xajKZWLx4EU8+Odzu/hMnvsOi\nRUuZPv0dMjMzWLBgDkOHPnNFXyoGA1gsRYckLmXjxnXUqFGTlSu/YNasBcyYMcXOZyj8If5eI8B1\n113PsmWfcuedvZk9e/olr3nx6Yo77tFHB7Jq1ZfExETzv//9Wuy5eva8i99+28+xY38W+ld+cbWX\nxOUOLX51W8eud/31N+Dn54e/v79jJyhGXNxJtm7dzKpVX/Lpp1/wxRefce7cWaxWK4GBQcyatYBu\n3W5n7tyZpKaeZ926r3j44QFOreFy1DMgIh7DaDTSufOtfPvtf7n22uvtDhNUrx7G8ePHbNsSE89Q\nvXp12+s//vids2eTGTXqeQCSk5N44YUXeOSRwVitVurXb0DNmrWIiIjk6NE/2bXrJ44cOQzA0aN/\n8vLLo5g1awFBQVWL1Fe3bn0OHNjHnXf2tm1LSUkhOzuLI0cO2x0m2LNnN23atAXguuv+QVJSIkFB\nVTl//rxtv6SkRKpXD+PYsaN2a0xJSaF167Z4e3vTtWs3PvtsVaHrFNcm27f/WOS41NTzHDlymObN\nW1Cpkh9t27Znz57d3HST/dVizWYz2dnZWK0W8vLyCA6uZrf2LVu22B0myMhI59tvNxEcXI3hw0dw\n9mxysX92+Z+lOmfPJhEQEGA7d/XqYZc97kp4eXkV2WZvmGDy5HcZOTL/56d//8eK7fEp8Ntv+2nU\nqImtB6Fhw2s5cuQwISEh3HxzCwBat27H8uUfsmvXT6SknGPYsCcwmXKJi4tj9ux3Cw3puILCgIi4\njcFgtHtz36Xs37+PunXrccMNN9odJjh16hQrV37MkCFPcf58ComJidSv/9eYc+PGTfjkk89sr/v2\nvZsZM2YQFbWGDRvW8vbb08jOzub48WPUqhXJf/6zxrbvs88+ySuvjLcbBABuv70njz/+MHfcsZdG\njZpgMpmYNm0Sbdu25667+tj90khOTmL//r106dKNU6cSqFzZH19fX+rVq8/u3b/SrFlzvv/+O+6/\n/0GOHv3Tbo3//vcS8vLy6Ny5C/v27aVOnXqFrtGiRSu7bfLee/OLHJeXl8dbb01g2bJP8Pf358CB\nfdxxR69i/zw++eQjunXrQU5ODp9++hEDBgy0W/ulhgn++c++tt9nZGSQkBBPWFg427f/yOuvFx6H\nb926Ld999w0DBz7Bli3f0qZNO2rVirjscY4qbpigJHf4R0bWYdWqT7BYLFgsFo4cOURERCRt27Yn\nJiaaJk1u4vffD1CnTj26du1ue7omISGet94a7/IgAAoDIuJG1atXJy/PxKuvjuGFF14q9gbCgnsG\nLBYLoaHVefnlN4o9Z82aNbn77j4MH57/ON+oUfmPIe7YsZ2EhPhCXzwX69y5C7GxP/H004PJzc1l\nwICBVKtWze6+ALNmvUu/fg8VehTM39+fd9+dxZQpk8jJycHLy4sePe4s9Iji39177328/fb/8eyz\nT2I25/HSS/k3ND7//EimTp2E1WqhUaMmtGrVBqvVarfG5557kcmT32TVqhW2RwT/XqO9NrF3XEhI\nKIMGPcHzzz+Nl5cX1157HR073kpyclKRP59TpxLYuvU7FixYgtVqZejQx+ne/Q67tV+pUaPGMn58\n/jVuu60HdevWK3Ttvn0f4s03X2PYsCcICAi0fenbO67gJtNTpxLw9vZm8+ZvmTRpKj/+uNX2aKsz\nFNwcevz4MX7//QBRUZ8yY8Y8li9fys03t6BJk5to3botw4blP+549919qFUrgvDwGrz77mSeeWYw\nXl7el7x51tW0aqE4TKuQuZ7auHQ40s6fffYf2rXrUKo3ebnbvHmzGD58hEPHlqWf5T//PMKBA/vo\n1etud5fiVFq1UESklIWH16hQQcBkMpXoX/hlWXZ2Fm3btnd3GWWKhglERBzQsWNnd5dQqnx8fGjd\nuq27y3CKvz/+KeoZEBERqfAUBkRERCo4hQERKdZL87fz0vzt7i5DRFxMYUBESl1mZiZ9+xZ/J/eW\nLd+WYjUiojAgInbF7D9NSnoOyanZvL44hpj9p0vlugkJ8XzzzcZSuZaI5NPTBCJSRMz+0yxas8/2\n+mRihu11m0Y1HDpnRkY6r7wymtzcXNvUtps2rScqaiVeXkbq12/ImDGvMH36Oxw4sI8PP3yf3r3v\n4c0381eoy8vL49VXJxS7aI6IOE49AyJSxNroo8VsP2Z3+5XYuHE911zTkPnzP+C66/4BQFZWFu++\nO4cFC5Zw/Hj+csYPP/wozZu3YNCgoSQnJzFo0FDmzFlE79738Nln/3H4+iJSPPUMiEgR8UmZdrcn\nJGc4fM6jR4/QvHn+3PQ335z/a1BQEOPG5c+7fuzYn5w/n1LomJCQUGbOnMbixYtIS0vl+utvdPj6\nIlI8hQERKSKiuj8nE4t+8dcKreLwOa1WMBoNAFgsVkwmE9OnT2Hp0hWEhlZn9Oh/FTlm8eJFtGnT\nlj59+rJ58zds3/6jw9cXkeJpmEBEiujdrn4x2+vZ3X4lChaOAYiN/ZnMzEy8vLwIDa3O6dOn+O23\nA+Tl5WE0/rWSYUpKCpGRtbFarfz44/eYTCaHry8ixVMYEJEi2jSqwVP3NMbrwr/ka4cF8NQ9jR2+\neRDgzjt7s2/fHkaMeIYTJ44RHBxMq1ZteOKJx/jww/fp3/9RZs+eTr16Dfj999+YPftd7r33PmbM\nmMrIkc/Trdsd/PprLDt37nDWxxSRC7RqoTisLK1CVl65u40LJhyaOqx8L+ri7nauCNTGrnc1qxbq\nngERKVZ5DwEikk/DBCIiIhWcwoCIiEgFpzBQAlq0RUREyiOFARERkQpOYUBERKSCUxgQERGp4BQG\nREREKjiFARERkQpOYUBERKSCUxgQERGp4BQGREREKjiFARERkQpOYUBERKSCUxgQERGp4BQGRERE\nKjiFARERkQpOYUBERKSCUxgQERGp4BQGREREKjiFARERkQpOYUBERKSCUxgQERGp4NweBrKzs+ne\nvTufffYZCQkJPProo/Tv358RI0aQm5sLwJo1a7j//vvp168f//nPf9xcsYiISPni9jCwYMECqlat\nCsDs2bPp378/K1asoF69ekRFRZGZmcm8efNYunQpy5cvZ9myZaSkpJR6nTH7T5OSnkNyajavL44h\nZv/pUq9BRETEFdwaBg4fPsyhQ4fo0qULADExMXTr1g2Arl27Eh0dze7du2natCmBgYH4+fnRokUL\nYmNjS7XOmP2nWbRmH2aLFYCTiRksWrNPgUBERMoFt4aBd955h7Fjx9peZ2Vl4evrC0BoaCiJiYkk\nJSUREhJi2yckJITExMRSrXNt9NFith8r1TpERERcwdtdF/7iiy9o3rw5derUsfu+1Wot0XZ7wsIC\nHart7+KTM+1uT0jOcNo1PFVF//ylQW1cOtTOrqc2LrvcFga2bNnCiRMn2LJlC6dOncLX1xd/f3+y\ns7Px8/Pj9OnThIeHEx4eTlJSku24M2fO0Lx58yu6RmJimlNqjQj152RiRpHttUKrOO0anigsLLBC\nf/7SoDYuHWpn11Mbu97VhC23DRPMnDmT1atXs2rVKvr168ewYcNo3749GzduBGDTpk106tSJZs2a\nsWfPHlJTU8nIyCA2NpaWLVuWaq2929UvZnu9Uq1DRETEFdzWM2DPc889x5gxY1i5ciURERH06dMH\nHx8fRo4cyZAhQzAYDAwfPpzAwNLtamrTqAYAK1dFk+bjT6S/gZ533GTbLiIi4skM1pIMwnsYZ3dJ\nffrsFJ5bNZH0CZPIeuZZp57bE6nbz/XUxqVD7ex6amPX88hhAk90tmoYAMb4k26uRERExHkUBkrg\nbFB1AIxumPRIRETEVRQGSuB8QAhPvvwFaXMWursUERERp1EYKAmDgVxfP3dXISIi4lQKAyUUmnIa\nn+83Q16eu0sRERFxCoWBEur7zYcE97sX46kEd5ciIiLiFAoDJWR7oiAuzs2ViIiIOIfCQAmdDcoP\nA14JCgMiIlI+KAyU0LmCxwvVMyAiIuWEwkAJJWviIRERKWcUBkrINkwQH+/mSkRERJyjTC1U5AnS\nqlQlZfVXmBtc4+5SREREnEJhoKQMBkydbnV3FSIiIk6jYQJHZWSAxeLuKkRERK6awoADAkaOIKxB\nLYwnT7i7FBERkaumMOAAa0gIAF7xerxQREQ8n8KAA8y1IgAwxunxQhER8XwKAw6wRNYGwKjHC0VE\npBxQGHCAOSISAC9NPCQiIuWAHi0sganD2gNgOZsMaEpiEREpHxQGHGCtFkL6q+MxN2rs7lJERESu\nmsKAIwwGsp5/0d1ViIiIOIXuGRAREangFAYcVHnBXKrd2g7jsaPuLkVEROSqKAw4yJByFu8D+/DS\nLIQiIuLhFAYcZIkomGtATxSIiIhnUxhwkCUyf64BhQEREfF0CgMOMl/oGfDSlMQiIuLhFAYcZIm4\nsD6BegZERMTDaZ4BB1mrBpNz172YmjV3dykiIiJXRWHAUQYDqUuWu7sKERGRq6ZhAhERkQpOYeAq\n+G5cT8DI5zHqJkIREfFgCgNXwfvXWCovX4rXkcPuLkVERMRhCgNXwRJ5YeIh9QyIiIgHUxi4CuYL\njxd66fFCERHxYAoDV8E2JXGcwoCIiHiuEj9aeOrUKZYsWcIPP/xAfHw8AJGRkXTq1ImBAwdSq1Yt\npxdZVtkmHkpQGBAREc9Vop6BqKgoBg0aRO3atZkzZw7R0dFER0cze/ZsIiMjGTJkCKtXr3ZVrWWO\nNagq5nr1sVYJcHcpIiIiDitRz8DBgwdZs2YNPj4+hbZfe+21XHvttTz00EO8++67Ti2wrDv70//c\nXYKIiMhVMVitVqu7i3CVxMQ0d5dQroWFBaqNXUxtXDrUzq6nNna9sLBAh4+9omGCnJwccnNzC207\nc+aMwxctT7x+/41Kn3yEITHR3aWIiIg45LJhYOnSpQwdOpSnn36aN954g+zsbABGjRrl8uI8QaU1\nnxM0Yhje+/a4uxQRERGHXPaegY0bN/LJJ58AsG3bNp5++mnGjx/v6ro8RsHEQ17xcZjcXIuIiIgj\nLhsGLBYLeXl5eHt706FDBxo2bMi4ceM4evRoKZRX9plrXXi8UBMPiYiIh7psGBg9ejTnzp0jLCwM\ngJo1a7Jo0SK+/vprlxfnCWxTEisMiIiIh7psGLjllluKbPP19eW+++5zSUGexhwRCWhKYhER8Vwl\nnoEQYO7cuSxfvpzg4GAaNWpE48aNady4Me3atXN2fWVfQACWqsHqGRAREY/lUBiIioqyzTS4f/9+\n9u3bx4cfflgxwwBwbvM2LNXD3F2GiIiIQxwKAw0aNKB27fyx8tq1a3P77bc7tShPY6ldx90liIiI\nOMyhVQufeeYZ3n33XXJycpxdj2fKzMTrj98hPd3dlYiIiJSYQz0DX375JV9//TUff/wx119/PU2a\nNKFJkybce++9V3yOrKwsxo4dS3JyMjk5OQwbNowbbriB0aNHYzabCQsLY+rUqfj6+rJmzRqWLVuG\n0WjkgQceoF+/fo6U7TL+C+dSZfJEUj79DNNt3d1djoiISIk4FAaio6PZvn073t7e/Pbbb+zdu5ed\nO3eWKAxs3ryZJk2aMHToUOLi4hg8eDAtWrSgf//+9OzZk+nTpxMVFUWfPn2YN28eUVFR+Pj40Ldv\nX3r06EFwcLAjpbvExU8UaOIhERHxNA6FgYYNG+Ln54eXlxfNmjWjWbNmJT5Hr169bL9PSEigRo0a\nxMTEMGHCBAC6du3KkiVLaNCgAU2bNiUwMH8BhhYtWhAbG8ttt93mSOkuYbkQBvREgYiIeCKHwkC1\natUYNWoUY8eOpUaNGldVwEMPPcSpU6dYuHAhgwYNwtfXF4DQ0FASExNJSkoiJCTEtn9ISAiJV7go\n0NWs4FQiTa8HoMrZM1QprWuWEaXWxhWY2rh0qJ1dT21cdjkUBsLDw9m/fz/33nsvAQEBtnkGnnzy\nyRKf69NPP+XAgQO89NJLXLyacnErK5dkxeVSWy6zUlXCgNzDf3K+Ai3RqSVJXU9tXDrUzq6nNna9\nqwlbDoWBi1csjI+Pt801UBJ79+4lNDSUWrVqceONN2I2m6lSpQrZ2dn4+flx+vRpwsPDCQ8PJykp\nyXbcmTNnaN68uSNlu46/P5Zq1TRMICIiHqlEjxb27NmTkSNH8vnnn3PmzBkAIiIi6N69OyNGjCjR\nhX/++WeWLFkCQFJSEpmZmbRv356NGzcCsGnTJjp16kSzZs3Ys2cPqampZGRkEBsbS8uWLUt0rdKQ\nNnshae/OcXcZIiIiJWawlqDf3Wq1smfPHrZt28a2bdvIyMigTZs2dOjQgdatW1OpUqUrvnB2djav\nvPIKCQkJZGdn8+yzz9KkSRPGjBlDTk4OERERvP322/j4+LBhwwYWL16MwWBgwIAB3HPPPVd0DXVJ\nuZa6/VxPbVw61M6upzZ2vasZJihRGNi3bx+NGze2vU5PT2fHjh1s27aNn3/+ma+++srhQlzBLT94\nVisYDKV/XTfQ/9yupzYuHWpn11Mbu16p3TPw4osvkp6eTvv27enUqRMdO3ake/fudO+uiXb8/v0h\nVf7vddIWfkBu9zvcXY6IiMgVK1EY2LhxIwkJCURHR7N161amTp1KWFgYnTp1olOnTmVyLL+0WP38\nMKaexxgf7+5SRERESqTETxPUqlWL++67j/vuuw+A77//ng8++ID33nuPAwcOOL1AT2GJzF+4yRh/\n0s2ViIiIlEyJw8DZs2eJjo5m27Zt7Nq1i/DwcNq0aVPipwnKG9uUxHF6vFBERDxLicLAPffcQ2Zm\nJr179+auu+7i9ddfx8/Pz1W1eRRLrQhAUxKLiIjnKVEYePDBB4mOjmb9+vUcPXqU48eP065dO+rV\nq+eq+jyHnx+W6tUVBkRExOOUKAw88sgjPPLII1gsFvbu3cv27dsZP348SUlJNGnShLfffttVdXqE\nzKefhQtrK4iIiHgKh6YjNhqNNGjQgFOnTpGUlMTZs2eJjY11dm0eJ+v5F91dgoiISImVKAzExMSw\nfft2tm/fzrFjx2jZsiUdOnTg8ccfp06dOq6qUURERFyoRGHg7bffpmPHjowaNYoWLVrg4+Pjqro8\nks933+D/3nwyn38RU/uO7i5HRETkipQoDIwfP/6yKwbu3r2bZs2aXVVRnsp47iy+331Dzh29FAZE\nRMRjlGjVwvnz5zNjxgzOnj1b5L1z584xY8YM5s+f77TiPI2lYK4BPVEgIiIepEQ9AwsXLmTJkiXc\nddddREZGUqtWLQDi4+M5deoUgwcPZsGCBS4p1BMUTDxkjNMshCIi4jlKFAaMRiNPPPEEgwYNYs+e\nPSQkJAD5UxQ3bdoULy8vlxTpKWwTDyVofQIREfEcDj1a6OXlRfPmzS97/0CF4+uLJSwcL/UMiIiI\nB3EoDEjxcm/tCiaTu8sQERG5Yk4LA2lpaQQGBjrrdB4rbf777i5BRESkRBwOA4cOHeLcuXMA5Obm\nMnHiRNavX++0wkRERKR0OBQGJk6cyLZt20hKSqJu3bqcOHGCwYMHO7s2j+S1fx+V1q4h985e5DWt\nmPMtiIiIZynRPAMF9uzZw/r167nhhhtYvXo1S5YsISsry9m1eSTv3w9QZerbeO/c4e5SRERErohD\nYcD3wsp8JpMJq9VKkyZNtFDRBeZaBRMP6fFCERHxDA4NEzRo0ICPP/6Yli1bMmjQIBo0aEBaWpqz\na/NIlkhNPCQiIp7FoTAwYcIEzp8/T1BQEGvXriU5OZmnnnrK2bV5JEvNWlgNBk08JCIiHsOhMGAy\nmVi7di0JCQmMGjWKX3/9lWon/z9BAAAgAElEQVTVqjm7Ns/k44OlRk284rQ+gYiIeAaH7hkYP348\nx48fJyYmBoD9+/czduxYpxbmySyRkZBnAqvV3aWIiIhclkNh4MiRI4wbNw4/Pz8A+vfvz5kzZ5xa\nmCdL+fq/nP31ABgM7i5FRETkshwKA97e+aMLhgtfdpmZmWRnZzuvKk9XwRdsEhERz+JQGLjzzjt5\n/PHHOXnyJBMnTqRPnz7cfffdzq7NYxmSk/HZ/C3G48fcXYqIiMhlOXQD4YABA7jpppvYuXMnvr6+\nTJ8+nSZNmji7No/l++P3BA0dSPpb75A19Bl3lyMiInJJDoWB1atXc//993PTTTc5u55ywRxxYa4B\nTTwkIiIewKFhgv/+97+aZOgSLJG1ATDGa+IhEREp+xzqGcjOzua2226jQYMG+Pj4YLVaMRgMfPzx\nx86uzyNZatTE6uWluQZERMQjOBQGhg0bVmRbvLrE/+LlhaVGTYzxCgMiIlL2ORQGWrduDUBOTg4b\nN25k9erVHD58mD59+ji1OE9miYjE+9dYMJv1qKGIiJRpDoWBX3/9lc8++4x169ZhsVh48803uf32\n251dm0dLmzkPq68vGB26LUNERKTUlOib6v3336dXr1688MILhISEsHr1aurWrUvv3r3x8fFxVY0e\nyfyP67HUb6BZCEVEpMwrUc/AzJkzufbaa3n99ddp27Yt8NcshFKUIS0VDAasAYHuLkVERKRYJeoZ\n2LJlC7179+aNN96gR48ezJ8/H5PJ5KraPJrvuq+p3rA2fh8tc3cpIiIil1SiMBAWFsaTTz7Jxo0b\nmTRpEsePHycuLo6nn36a77//3lU1eiRLzZqAJh4SEZGyz+G721q1asXkyZP54Ycf6NKlC/PmzXNm\nXR7vr4mH9HihiIiUbSW6Z6Bnz540atSIjh070qFDB8LDwwkICOChhx7ioYceclWNHslSPQyrtzde\ncZqFUEREyrYShYF169axZ88etm3bxosvvkhGRgZt2rShQ4cOtG7dmkqVKrmqTs/j5YWlVgTGBA0T\niIhI2WawWq1WRw9OT09nx44dbNu2jZ9//pmvvvrKmbVdtcRE966fEHz3HXj/FEPSySTwdmhKhzIt\nLCzQ7W1c3qmNS4fa2fXUxq4XFub4k2tXtWphQEAA3bt3p3v37g4XUJ5lDh+BIT0NLBZ3lyIiIlIs\nrVroQrl39iKn74Pg6+vuUkRERIrllFULC2jVQhEREc/jtFULpSjvnTEEjnyOrEFDyR481N3liIiI\n2OXQMEHz5s05ePAgW7dutT1F0KxZM2fX5vkq+eL9+294HTnk7kpERESK5VAYGD9+PMePHycmJgaA\nffv2MXbsWKcWVh6Ya0UC4BWniYdERKTscigMHDlyhHHjxuHn5wdA//79OXPmjFMLKw+s1atj9fXF\nmKAwICIiZZdDYcD7wjPzBSsWZmZmkp2d7byqygujEUvNCIzqGRARkTLMoRsI77zzTh5//HFOnjzJ\nxIkT2bp1K/379y/xeaZMmcKuXbvIy8vjqaeeomnTpowePRqz2UxYWBhTp07F19eXNWvWsGzZMoxG\nIw888AD9+vVzpGy3MEdG4rNjO+Tm6hFDEREpkxwKAwMGDOCmm25i586d+Pr6Mn36dJo0aVKic+zY\nsYODBw+ycuVKzp07xz//+U/atWtH//796dmzJ9OnTycqKoo+ffowb948oqKi8PHxoW/fvvTo0YPg\n4GBHSi91ub3uwtyoMYbcHKwKAyIiUgY5FAbuvvtuOnbsSMeOHWnVqhW+DnzJtWrViptuugmAoKAg\nsrKyiImJYcKECQB07dqVJUuW0KBBA5o2bUpgYP40iy1atCA2NpbbbrvNkdJLXdZTw91dgoiIyCU5\nFAY+/PBDYmJi2LhxI1OmTCEsLIyOHTsycODAKz6Hl5cX/v7+AERFRdG5c2d+/PFHW7AIDQ0lMTGR\npKQkQkJCbMeFhISQmJh4Rde4mnma5cqojV1PbVw61M6upzYuuxwKA9WrV6d3797ccsst7Ny5k3Xr\n1rFo0aIShYEC33zzDVFRUSxZsoTbb7/dtr249ZNKsq5SWVgUw+vwQSrPnUVu127k3vNPd5fjVFp4\nxPXUxqVD7ex6amPXK/WFil5++WVOnDhBWFgYt9xyCy+88ALXX399ic/zww8/sHDhQj744AMCAwPx\n9/cnOzsbPz8/Tp8+TXh4OOHh4SQlJdmOOXPmDM2bN3ekbLcwZGRQ+eN/Y61cudyFARERKR8cerQw\nMzMTgICAAIKDgwt141+ptLQ0pkyZwqJFi2w3A7Zv356NGzcCsGnTJjp16kSzZs3Ys2cPqampZGRk\nEBsbS8uWLR0p2y3MEbUBTTwkIiJll0M9AzNnzgTg999/Z+fOnYwbN464uDjWr19/xedYt24d586d\n41//+pdt2+TJk3n11VdZuXIlERER9OnTBx8fH0aOHMmQIUMwGAwMHz7cdjOhJ7CGhmL188MYrzAg\nIiJlk8FakkH4C9LT09m1axc7d+4kNjYWq9VK69atefHFF11Ro8PKyvhUtTbNMaank7yvfK1RoDFA\n11Mblw61s+upjV2v1O8ZuPfee2nfvj3t27dn6NChHvPMv7tYIiLx3vYD5ORApUruLkdERKQQh+4Z\n2LhxI7fccgu7du1i7ty5rF271tl1lSt5jRpjanYzhtRUd5ciIiJShEM9A2+99RZnz56ldevWWK1W\n1q1bxy+//MKrr77q7PrKhYy3pri7BBERkWI5FAYOHjzIRx99ZHs9YMAAh9YmEBEREfdzaJjAZDJh\nsVhsr81mM2az2WlFlTeGpCT8PlqG945od5ciIiJShEM9A7feeit9+/alVatWAMTExNCrVy+nFlae\nGE8lEPjic2QNeoL0tu3cXY6IiEghDoWBYcOG0b59e3bv3o3BYOD//u//bIsOSVGWyEgAjAnxbq5E\nRESkqBKFgfT0dObPn8+RI0do2bIlAwcOxNvboTxRoViDq2H198eoWQhFRKQMKtE9A+PHjwfgwQcf\n5PDhw8ydO9cVNZU/BgPmWhF4xZ90dyUiIiJFlOif9XFxcUybNg2Azp07O7RKYUVliaiN9+FDkJUF\nlSu7uxwRERGbEvUMXDwk4OXl5fRiyjNLZCRWgwHj6VPuLkVERKSQEoUBg8FwyddSvPSJk0k6kYil\nfgN3lyIiIlJIiYYJfvnlF7p06WJ7nZycTJcuXbBarRgMBrZs2eLk8soPa1BVd5cgIiJiV4nCwIYN\nG1xVR/mXm4vXkcNY/fzUOyAiImVKicJA5IXn5aXkvA4fIuTWtmQ9PoT0qTPcXY6IiIiNQ9MRS8nZ\nJh7S44UiIlLGKAyUEmtgEJYqAXhp4iERESljFAZKi8GAJTISY4LCgIiIlC0KA6XIEhGJ8dw5yMx0\ndykiIiI2CgOlyBxZGwCvePUOiIhI2aEwUIqyhj5DypoNmCP0VIaIiJQdWnKwFJkbNcbs7iJERET+\nRj0Dpc1qBZPJ3VWIiIjYKAyUIuOfRwi9JpKAsSPdXYqIiIiNwkApsoaFYcxIxytOEw+JiEjZoTBQ\niqwBgViCqmLU0wQiIlKGKAyUMktEBLl/Huel+dvdXYqIiAigMFDqLBGR+Odk4JejiYdERKRsUBgo\nZQUTD4WcT3RzJSIiIvk0z0Apy+l9N2sSvMioHOjuUkRERACFgVJnuq0H63+r4u4yREREbDRMUMpi\n9p8mJT2H5NRsXl8cQ8z+0+4uSUREKjj1DJSimP2nWbRmn+31ycQM2+s2jWq4qywREang1DNQitZG\nHy1m+7FSrUNERORiCgOlKD7J/uOECckZpVyJiIjIXxQGSlFEdX+722uF6oZCERFxH4WBUtS7Xf1i\nttcr3UJEREQuojBQito0qsFT9zTGy2jAaDFTP+kYT915rW4eFBERt1IYKGVtGtUgOKAS9x7eytvX\npNP22mruLklERCo4PVroJltv6cndw9q7uwwRERH1DLhdXh5YLO6uQkREKjCFATeqtOoTQptci88O\nLWcsIiLuo2ECN5h6YXjAsjkL49mzxEx+n0/uNNq2i4iIlCb1DLiRqUMnLIFB3HJgG1it7i5HREQq\nKIUBd/L1JbfH7VQ/f4a6p464uxoREamgFAbcLLfnXQC0+E33DYiIiHvongE3+6FOC9Y/NosTIbXZ\nuDiG3u3qaxIiEREpVQoDbhSz/zSLNh2B6vnTEWtJYxERcQcNE7iRljQWEZGyQGHAjbSksYiIlAVu\nDwN//PEH3bt356OPPgIgISGBRx99lP79+zNixAhyc3MBWLNmDffffz/9+vXjP//5jztLdhotaSwi\nImWBW8NAZmYmb775Ju3atbNtmz17Nv3792fFihXUq1ePqKgoMjMzmTdvHkuXLmX58uUsW7aMlJQU\nN1buHFrSWEREygK3hgFfX1/ef/99wsPDbdtiYmLo1q0bAF27diU6Oprdu3fTtGlTAgMD8fPzo0WL\nFsTGxrqrbKe5eEljgHrZyby0dhodsk+4uTIREalI3Po0gbe3N97ehUvIysrC19cXgNDQUBITE0lK\nSiIkJMS2T0hICImJiaVaq6u0aVSDqC2HAZjYNpTg+T+SveR9nt6RPzyiKYpFRMTVyvSjhdZipugt\nbvvfhYUFOrMcl/Hyyu8ZCL7/brjuOvy+/Iyq1/yT9CpVy/xnKOv1lQdq49KhdnY9tXHZVebCgL+/\nP9nZ2fj5+XH69GnCw8MJDw8nKSnJts+ZM2do3rz5Zc+VmJjmylKdZvJT+fdMJCZnUPmxwfwc9R2n\nMy1k5GbxzORvyuxERGFhgR7Txp5KbVw61M6upzZ2vasJW25/muDv2rdvz8aNGwHYtGkTnTp1olmz\nZuzZs4fU1FQyMjKIjY2lZcuWbq7UNba0vpt3uw0nw9sP+Gsiopj9p91cmYiIlFdu7RnYu3cv77zz\nDnFxcXh7e7Nx40amTZvG2LFjWblyJREREfTp0wcfHx9GjhzJkCFDMBgMDB8+nMDA8tndtDY2wf72\n6GNlsndAREQ8n8F6pQPwHsgTu6SeeGczFjt/JF5GA++P7uqGioqnbj/XUxuXDrWz66mNXa9cDRNU\ndJqISERESpvCQBmjiYhERKS0lbmnCSq6gvsCPvh6P9a8POomH+euG4JoofsFRETERRQGyqCCiYiC\nU9OYEfUyluBqPJMRSq6vnyYhEhERp9MwQRmWElSdrCeHsS3oGjLPp5Gcms3ri2P0mKGIiDiVegbK\nqIIegM27arHQq5Nte8G8A4AeNRQREadQz0AZ9/WvZ+xuXxt9rJQrERGR8kphoIyLT8q0uz0hOaOU\nKxERkfJKYaCM07wDIiLiagoDZZzmHRAREVdTGCjj2jSqwVP3NMbLaACs1E+JY9T66XTIPuHu0kRE\npJxQGPAAbRrVIDigEqFBlXnzrnrcemArgc8+xSszv+Ol+dvdXZ6IiHg4PVroIQoeNTQBmU8N46ct\ne0hLyyLdB15fHEPvdvX1qKGIiDhEYcADfXf/cBZVOWh7rbkHRETkamiYwAOt3RVvf7vmHhAREQco\nDHggzT0gIiLOpDDggYqbe8BssTJ48ndav0BEREpEYcADFTf3QIGCewgUCERE5EooDHigwnMPgI+X\nwe5+uodARESuhJ4m8FBtGtWwPTnwxDub7e6jewhERORKKAyUAxHV/TmZWPSLv2qAL0OnbMZssVI7\nrIrmIhAREbs0TFAOFHcPwdnUHMwWK6D7CEREpHgKA+XA3+8hqB0WQCg5dvddsvaAAoGIiBSiYYJy\n4uJ7CACeeOc7sBbdz2S2aLZCEREpRD0D5VRE9SqXfF89BCIiUkBhoJy63FwEBT0ECgQiIqJhgnKq\nYAhgydoDmMyWYvdbsvZAof1FRKTiUc9AOdamUQ0G977xkvuoh0BERBQGyrmCJw18vC79R71ozT6t\naSAiUkFpmKACKBgCKHiKoDgFcxH8Z/MhUtJziajur4mKREQqAIWBCuJK7yEAOJuWP0dBQTi4+HgR\nESl/FAYqkCvtIfi7JWsP8P5X+9VTICJSTikMVDAFX+Rro49xMjH9io4p6ElQT4GISPmkGwgroDaN\navB/Q1rz1D2NHTpeSyOLiJQv6hmowC7uJUhIzqBqgC9nU+2vaXCx+KR0Xl8cQ3xyJhGhGjoQEfF0\nCgMV3N/XNIjZf9oWDowGg92bDS1WbEsma+hARMTzKQxIIReHg5j9p6/4ZsN13+ylzXUh4ONjO3Zt\n9FHikzJ146GISBmnMCDF+vswQq3QKsQlpttbDJH4NDOhTa4l+5HH+a7v8EIhQr0HIiJlm8KAXNLf\nhxFeXxxjGyK4WCSZWCv5YUhLY230UbvnWht9TGFARKQMUhiQEundrr7doYOe/2zH2XEHMGSkEz9/\nl91jE86k4rvua0ztO2ANrlZkKOH6utX4/fg5DS2IiJQyhQEpEXtDB73b1bNttwYGEVHd327vQZ2k\nY1Qd+ALZfe7jm3+9U2Qo4eJjCoYWDsWdV0AQEXExg9VqtTcEXC4kJqa5u4RyLSws0G4bF3fj4bAb\nven8xzbymjVn7IkQu4HhSoQEVqowaycU18biXGpn11Mbu15YWKDDx6pnQJyuuN6Dlo1qkElnAOLf\n2ezw+f++doJ6D0REro7CgLjE3288/LvihhIc8e2uk7bfF/fkgh51FBEpnsKAuEVxNyI6y/o1P9F5\nXzrmeg3YVqUui746YHtPgUFEpDCFAXELe0MJ19cN5vfjKYVeX/yv/pKIM/sSNPxJrH5+rH3jK7v7\nrF/3C51OWjHXimS7JYRF6/+wvXelcyMoQIhIeaAwIG5zuaEEgGsjq5Z47QSAiAAv0qbMwJCWRnxS\npt194nKMVB3QF6D4wPDlTrpsO4mlehimFi2x1Kuf/4bVSsyBMw5NrqQAISJljcKAlGmXWjvhUr0H\nvXo0IbtRNwAiipsoqbKV9FfH45UQT3y62e714yyVCHhlDABpU2eS/fhgAIJ7dWfDTY9DtdpFjimY\nXMnr998wnjmNNTgYS9VgrNWqseN4Jou+2m/b97JDFpdYDEqhQkScRWFAPIq93oSLew/+Pu8BXGKi\npDubk9XoDqD4wBAR5MP5xf/GmJiIqV0H23Zz7TqcqFrLbo0Jyfnn8fv3EvzfX1jovQ2PzYLq9Yoc\ns2H1Nm5bEYM1IJBtfnVYlF3H9l5BYDCcOU2bG2tgDQgg5mjaFYUKe1wVIhRORDyXwoB4vMsNN1xu\noiQoPjD06taIXDvnTnt/abEBolZoFQByu9+BNbgahpRzGFNSMJxP4URonSL7A5w0BuD/3gIAvnp0\nJoQV3WfDVz/R+4EX8n//6EwIq190n0+30OPo51grV8bqXwVr5crkNW1GzsMDAPh5QwyLfi06uZPX\niWO0vKN1/sbcXAwZ6Vj9KkOlSmA02q35Yn+fW8KZ61G4MmRcSQ+MSEWgMCAVgjMCw98VFyB6t8v/\nl7+pazdMXbsVeq/YHohQf859+wOGtDSOb7V/X8SJ6vXIvvc+DBnpnKhe1+4+JytVo9KGdYW25fS+\nxxYGvv7lDBiqFDlu3X/32cKAz66fCL63p+09q68v1kp+UKkS5z9aSV6LlgBUvbcnGAzg58f6G/qD\nX0iR867f8Ctdv94Fvr5YfXzB1werjy+WiAhyb8+/hvHYUbwP/XHhfd/86134/XZTVRZ9bedJEIuZ\nNo1r5V/fQa4MMK7iib0vnha4PLGNncFjwsCkSZPYvXs3BoOBl19+mZtuusndJUk5cyU3NP59f3BO\ngOjV+TryLhwXcaCYHofwINLGLc3fp7heiRpBJB08jiEzE0NWJmRkQhV/2/snjVWwt+zkxT0WlqCq\n5PS8C0NONuTkYMjOwpCdAznZ+aEAwGrFJ/ZnDDkXJoC66Rm7nzcux0iVGZOLbDe1aWcLA5U2rCXg\ntXF2j1/71ka72zd+uJ67lr+A1dub1CUfkXtnLwCCu3fGmJQI3t5Yvb3B2xu8fcj+5/1kjRgJQOV5\ns/Hd8i3rGw2AyqFFrxl9jA65cfh9vBy8jPnnMXrln8vLi8xhz2GtFgLZ2VReNO+i94xYvfL3MbXr\ngPn6GwDw+e6/GM+f/+s8Xl7gZcRSoyZ5NzUHwHjiOMaEBPDOf99asJ+PD+br/gFAzP9Osmhd0Sde\nDNlZtG5WJ39/AJMpvzfHaLyqsOQMnha4PK3eAhcHri+n3uPQOTwiDOzcuZNjx46xcuVKDh8+zMsv\nv8zKlSvdXZaISwLE5XocLrlP+wZYqwZjrRps9/oR1asUGzQKmBs3IXXZikt/EIOBpBOJYLFATg4R\ny3/h5NnsoterWomUz9fmDz3k5ub/asrFElzNto+pdVvSX52AwZQLplwMuaYLv+YSfzbL7uVPhNYl\nt10HDCYT1uCLPqufH/j4QF4ehuxsDOY8MOVhPHvWtov3b/vx/X4zJ29+1u65E5Iz8DqcSOVli+2+\nnzXgcazVQjBkZhDw1gS7+6RNmWELA1UmvYnP/34tsk9Or7tJXfpxftnLl1Jl5rQi+1j9/Uk6egqA\ndd8ftHutDR9/S/u8JuS1bgNA9Qa18tu64BwXgkHmi6PJHDUWgMChA/Hd/C0YDfkhwmDEajSSd3ML\nUpfn/91aadUnVJk2GauX11/hwmgEg5FzX2+CgACMJ08QNPRxMOS/d/G+maPGsna//a+YDau20uP4\nGnJv7Wq7Kdfv43/j8+PWv0KM0YjVYMAaEkrGG28C4L1nN37LlxZ6v2D/zOdHYg0Lg5wcqrzzlm27\n1WjI399gJLf77eTd0ir/8336MV6nEvLbh/zzrMu8DvArUu/6tbHcuuPkhfMYwNub7AGPA2A4c4ZK\n32zMr6UgfF34fe6tt2Gtkf//tu9XX+b/PBoMF+1rwHzdPzDfcGP+5/tlF8YzZ8Dw1zmsRiPWwCDy\nWuX/+RoT4vE6+Ift/ehEK4t+t38DdEl4RBiIjo6me/fuADRs2JDz58+Tnp5OQECAmysTKTlnDFk4\n0isBVxY0SsRohMqV6d2xof0ejy7XY7pMTXk330LezbfYfa/4HpCqnP9yfZHtKV9vumzJaTPnkTZt\nFhH/juVkctGwUSu0Crk9WnD2h51gNuf/BW425/+XZ8ZSoyYA1oBAUj5djeHCdsxmsJgx5OWR1/xm\n2/kyR4zM76248B5mC5jNmK9paNvH1LY9mc+/eOEaefn7ms1YfXxs+8Rl2l9G5kT1ulirXRSuOnbO\n7x2wWMBiwXDhV0vNv254tYSFYalTFyxm236YzX/1/EB+LdnZGArOYb2wn+WvOgw52Xjv3WM73mCx\n2N7LGvgE8Ul/hcyLnfSpSqWvv8RyUd3eu37Cb/WqIvua69S1hQHj0T+pvNR+SMseOARzWBgGUy7+\nc2fa3ccSGmoLA5WXLcFn10+F3o/712qwc4tMXK4XgWNH2l5bfXxsYcDr6J8E/mu43eulfL4W04Uw\nEPTc0xgyi/4sZ44YScYrbwDgP2s6ldYVfczZdFNzUr7ZCoDvpg0EvvQv23vrirl/qKQ8IgwkJSXR\nuHFj2+uQkBASExMVBqTcKggMl1rcpaS9EgXHQMlDhLvO6/TwAhe66b3o3eGaYs9tDaqKOajqpc/j\n64vpth6XvVzu3fdedh/Tbd0x3db9kvsU36tT1TaUAHD+088ue72MSVO53GTgOQ8PsN1rUhxzw+vy\ne4guZrXmBwmjkYgzP9m/RyYsgKQDf0Il379qev3/8nsuCsKJ1Zr/a8HwB2Dqclt+SLv4fasVg9WC\nOTJ/qMvqV5lz67658L7V9j5WK+b6DWznSp84GcP584D1QtixErnPwImiHVxE+htIfe/D/HNarYWG\nX8wNriF19gLbZzcU7GO1Yr72ukLXIzc3/z2stn1NzVrY9sl+6BFMrdv+9fku7GsJC7ftk9f8ZjJG\nv2x7/7jlKv5fuIhHrFr42muvceutt9p6Bx5++GEmTZpEgwYNLnOkiHi6rb+c5D/fHuTE6TTq1Aik\nX7fr6Hxz0fkdytq5nW3rLyeZ+tGuIttfGnCLanYST6sX4LlpmzmakGp7/dW7lw+f9nhEGJgzZw5h\nYWE89NBDAHTr1o0vv/xSPQMiIiJOcPkHiMuADh06sHFj/l3F+/btIzw8XEFARETESTzinoEWLVrQ\nuHFjHnroIQwGA2+88Ya7SxIRESk3PGKYQERERFzHI4YJRERExHUUBkRERCo4j7hnoCQ0bbHrTJky\nhV27dpGXl8dTTz1F06ZNGT16NGazmbCwMKZOnYqvr+/lTySXlJ2dzV133cWwYcNo166d2tjJ1qxZ\nwwcffIC3tzfPP/88119/vdrYyTIyMhgzZgznz5/HZDIxfPhwwsLCGD9+PADXX389EybYn71RLu2P\nP/5g2LBhDBw4kAEDBpCQkGD353fNmjUsW7YMo9HIAw88QL9+/S553nLVM3DxtMVvvfUWb731lrtL\nKjd27NjBwYMHWblyJR988AGTJk1i9uzZ9O/fnxUrVlCvXj2ioqLcXWa5sGDBAqpWzZ/wRm3sXOfO\nnWPevHmsWLGChQsX8u2336qNXeDzzz+nQYMGLF++nFmzZtn+Pn755Zf59NNPSU9P5/vvv3d3mR4n\nMzOTN998k3bt2tm22fv5zczMZN68eSxdupTly5ezbNkyUlJSLnnuchUGipu2WK5eq1atmDVrFgBB\nQUFkZWURExNDt275q/J17dqV6Ohod5ZYLhw+fJhDhw7RpUsXALWxk0VHR9OuXTsCAgIIDw/nzTff\nVBu7QLVq1WxfPqmpqQQHBxMXF2frqVU7O8bX15f333+f8PC/ZiS09/O7e/dumjZtSmBgIH5+frRo\n0YLY2NhLnrtchYGkpCSqXTTXdcG0xXL1vLy88PfPX/0uKiqKzp07k5WVZetODQ0NVVs7wTvvvMPY\nsWNtr9XGznXy5Emys7N5+umn6d+/P9HR0WpjF+jduzfx8fH06NGDAQMGMHr0aIKC/lqnQO3sGG9v\nb/z8Ci+kZO/nNykpiUvw8fAAAAc4SURBVJCQv5YUv5LvwnJ3z8DF9NSk833zzTdERUWxZMkSbr/9\ndtt2tfXV++KLL2jevDl16tSx+77a2DlSUlKYO3cu8fHxPPbYY4XaVW3sHF9++SUREREsXryY3377\njeHDhxMYGGh7X+3sGsW165W0d7kKA+Hh4SQlJdlenzlzhrCwMDdWVL788MMPLFy4kA8++IDAwED8\n/f3Jzs7Gz8+P06dPF+q6kpLbsmULJ06cYMuWLZw6dQpfX1+1sZOFhoZy88034+3tTd26dalSpQpe\nXl5qYyeLjY2lY8eOANxwww3k5OSQl5dne1/t7Dz2/o6w913YvHnzS56nXA0TaNpi10lLS2PKlCks\nWrSI4Avrx7dv397W3ps2baJTp07uLNHjzZw5k9WrV7Nq1Sr69evHsGHD1MZO1rFjR3bs2IHFYuHc\nuXNkZmaqjV2gXr167N69G4C4uDiqVKlCw4YN+fnnnwG1szPZ+/lt1qwZe/bsITU1lYyMDGJjY2nZ\nsuUlz1PuZiCcNm0aP//8s23a4htuuMHdJZULK1euZM6cOYVWipw8eTKvvvoqOTk5RERE8Pbbb+Nz\n0frr4rg5c+YQGRlJx44dGTNmjNrYiT799FPbEwPPPPMMTZs2VRs7WUZGBi+//DLJycnk5eUxYsQI\nwsLCeP3117FYLDRr1oxx48a5u0yPs3fvXt555x3i4uLw9vamRo0aTJs2jbFjxxb5+d2wYQOLFy/G\nYDAwYMAA7rnnnkueu9yFARERESmZcjVMICIiIiWnMCAiIlLBKQyIiIhUcAoDIiIiFZzCgIiISAWn\nMCBSjo0aNYrPPvuMxMREnn/++Uvu+9VXX2GxWEqpssLGjRvHhg0bHD7+5MmTPPzww1qLRMRBCgMi\nFUBYWBizZ8++5D5z5sxxSxhYv349OTk53HnnnQ6fo3bt2vTp04epU6c6sTKRikPzDIiUIxaLhVde\neYXff/+dyMhIMjMz6d27N61bt6Z///5s3bqVdevWsXjxYvz9/bFarbz99tt8/vnnzJs3j1atWjF3\n7lzWrVvHl19+iY+PD5UqVWLGjBkEBQVx22238dhjj7F161ZOnjzJhAkTaNeuHUePHuW1117DYrFQ\nqVIl3n77bWrUqMHy5ctZv349ZrOZa665hjfeeKPIQiv33HMPkyZNokmTJmRnZzNu3DgSEhIAePHF\nF2ndujU333wzzzzzDN999x0mk4mnn36aVatW8eeffzJ+/Hg6duyIyWSia9eurFmzptAiLSJyeeoZ\n+P/27iYU2i4M4Pj/9hlFM8m4mcLCbpSMsiJlMSU2BuWjQSlFkY2yMIONjLClRMJoEAuFBVkY2U6m\nZsHCQqYGKfJRNMyz0HtnXs/z5nk99ZS5fqv5ONc5557FnKtz7plLiG/k8PCQ09NT1tbWGB0d5fj4\n+EObqakpHA4HCwsL9Pb2cnFxoR0hzM3NodPpeHp6YmZmhsXFRYxGIxsbG1p8YmIis7OzdHR0MD8/\nD8DAwABtbW24XC5qamrY3t7G5/Oxs7ODy+VieXmZlJQUVldXI+ZyeXnJ1dUVJpMJgJmZGVRVxe12\nMzIyorV/fHwkPz8ft9tNcnIye3t7TE9P09nZydLSEgDx8fGYzWYpjSvE//CtChUJEe1OTk4oLCxE\nURSSkpK0+vHvWa1W+vr6sFgsWCwWCgoKPrTR6XS0t7cTExNDIBCIKPhVXFwMQFZWFre3twD4fD7t\n9crKSgCmp6c5OzujubkZeFvQ4+Iiv3KCwSCqqqIoitZPQ0MDALm5uRHb/kVFRQBkZGRgNpsBUFWV\nu7s7rY3RaCQQCHz68xJCvJFkQIhvJBwOawsr8NN7AFpbW6mqqsLj8eBwOKirq6O+vl57PxgM4nQ6\n2dzcJC0tDafTGRH/fkF/f8r477ESEhIoLy/H4XB8ev6KovzyvoXY2NifPhZCfJ0cEwjxjeTl5XF0\ndEQ4HOb+/l6rHPePl5cXxsbGSElJobq6mq6uLq2NoiiEQiGur6/R6/WkpaVxc3PDwcEBz8/P/zmu\n2WzG4/EAsLW1xcTEBGazmf39fR4eHgBwuVx4vd6IOFVVCQaD2vPCwkKtn/Pzc1paWn7r+gOBAEaj\n8bdihBCyMyDEt1JSUsLGxgZ1dXVkZWV9qGEeGxuLXq+nvr6e1NRUAPr7+wEoLS2lpqaGyclJcnJy\nqK2tJTs7m+7ubgYHBykrK/vluHa7HbvdztLSEnFxcQwPD5OZmUlTUxM2m43ExEQMBgNWqzUizmAw\nkJ6ejt/vx2QyYbPZsNvtNDY28vr6Sk9Pz6evPRQK4fV6GRwc/HSMEOKN/JpACPFXbW9vs7u7y/j4\n+Jf6WVlZwe/3MzQ09IdmJkT0kGMCIcRfVVFRQUJCwpf/dGh9fZ3e3t4/ODMhoofsDAghhBBRTnYG\nhBBCiCgnyYAQQggR5SQZEEIIIaKcJANCCCFElJNkQAghhIhykgwIIYQQUe4HfCzQAD3uAxMAAAAA\nSUVORK5CYII=\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7ff6f852de10>"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "stream",
          "text": [
            "Model parameters:\n",
            "B   = 0.4344889373542267  +/-  0.3257912873481228\n",
            "C   = 8033.032912318835  +/-  3788.0965208151465\n",
            "x_0 = -0.0011326268137890247  +/-  0.8856314393812412\n",
            "n   = -1.863839035112852  +/-  0.1314285599991098\n",
            "\n",
            "Fit quality:\n",
            "chisq = 2.4105856961077214 dof = 31\n"
          ],
          "name": "stdout"
        }
      ]
    }
  ]
}