{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Lecture 14 " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Table of contents\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "\n", "In this lecture, we will introduce various types of experimental error and solve some error analysis problems.\n", "\n", "## The simple pendulum\n", "### Theory\n", "\n", "Figure 1 shows a schematic diagram of a simple pendulum of length $l$ \n", "with a bob of mass $m$. The bob is at an angle $\\theta$, and an arc-length $s$, \n", "from its equilibrium position. The combination of the gravitiational force, $mg$, \n", "and the tension in the string, $T_s$, ensure that the pendulum moves along the arc indicated in the figure. The restoring force, $F$, is the force that acts towards the pendulum's equilibrium position.\n", "\n", "![Simple pendulum](SimplePendulum.png)\n", "\n", "**Figure 1:** *Simple pendulum illustrating the forces acting on the bob and the definitions of the arc $s$ and angle $\\theta$.*\n", "\n", "The force $F$ is given by:\n", "\n", "$$\n", "F = - mg \\sin \\theta.\n", "$$\n", "\n", "Note that $F$ acts towards the equilibrium position, in the negative $s$ direction in Figure 1.\n", "\n", "Newton's second law states that $F = ma$, where $a = \\frac{d^2s}{dt^2} = \\ddot s$.\n", "Since $s = l \\theta$, we have:\n", "\n", "$$\n", "m \\ddot s = m l \\frac{{{d^2}\\theta }}{{d{t^2}}} = m l\\ddot \\theta.\n", "$$\n", "\n", "Equating this to the restoring force gives $-mg \\sin \\theta = ml \\ddot \\theta$.\n", "Hence:\n", "\n", "$$\n", "\\ddot \\theta = - \\frac{g}{l}\\theta.\n", "$$\n", "\n", "The acceleration of the bob is proportional to its displacement from equilibrium and directed towards its equilibrium position. This is an example of simple harmonic motion (SHM), a phenomenon that is encountered in many areas of physics. The motion of the bob can be described by the equation $\\theta = \\theta_0 \\sin(\\omega t + \\phi)$. Here, $\\omega$ is the angular frequency of the bob's motion, $\\theta_0$ the amplitude and the value of $\\phi$ is determined by the *initial conditions* (the position and speed of the pendulum at $t = 0$). \n", "\n", "The angular frequency is given by:\n", "\n", "$$\n", "\\omega = \\sqrt {\\frac{g}{l}},\n", "$$\n", "\n", "and this in turn is related to the period of the pendulum's oscillations by:\n", "\n", "\\begin{align}\n", "T &= \\frac{2 \\pi}{\\omega} \\\\ \n", " &= 2 \\pi \\sqrt {\\frac{l}{g}}.\n", "\\end{align}\n", "\n", "## Measurement of $g$\n", "\n", "The period of oscillation of a simple pendulum with a constant mass and variable length is determined as a function of the pendulum's length." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "### Raw data\n", "\n", "Five measurements of 20 oscillations for each of 10 pendulum lengths gave the results below. " ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "#\n", "nLen = 10\n", "nMeasPerLen = 20\n", "lengths = np.zeros(nLen)\n", "Tmean = np.zeros(nLen)\n", "Tdata = np.zeros((nLen, nMeasPerLen))\n", "#\n", "# Measured lengths and times\n", "lengths[0] = 0.1500\n", "Tdata[0, 0] = 0.7876\n", "Tdata[0, 1] = 0.7671\n", "Tdata[0, 2] = 0.7690\n", "Tdata[0, 3] = 0.8073\n", "Tdata[0, 4] = 0.7555\n", "Tdata[0, 5] = 0.7742\n", "Tdata[0, 6] = 0.7558\n", "Tdata[0, 7] = 0.8065\n", "Tdata[0, 8] = 0.7632\n", "Tdata[0, 9] = 0.7500\n", "Tdata[0, 10] = 0.7797\n", "Tdata[0, 11] = 0.8078\n", "Tdata[0, 12] = 0.7855\n", "Tdata[0, 13] = 0.7870\n", "Tdata[0, 14] = 0.7854\n", "Tdata[0, 15] = 0.7876\n", "Tdata[0, 16] = 0.7629\n", "Tdata[0, 17] = 0.7256\n", "Tdata[0, 18] = 0.7934\n", "Tdata[0, 19] = 0.7944\n", " \n", "lengths[1] = 0.2000\n", "Tdata[1, 0] = 0.8908\n", "Tdata[1, 1] = 0.9079\n", "Tdata[1, 2] = 0.9031\n", "Tdata[1, 3] = 0.9162\n", "Tdata[1, 4] = 0.9197\n", "Tdata[1, 5] = 0.9046\n", "Tdata[1, 6] = 0.9005\n", "Tdata[1, 7] = 0.8523\n", "Tdata[1, 8] = 0.8949\n", "Tdata[1, 9] = 0.9214\n", "Tdata[1, 10] = 0.8598\n", "Tdata[1, 11] = 0.8714\n", "Tdata[1, 12] = 0.8703\n", "Tdata[1, 13] = 0.9229\n", "Tdata[1, 14] = 0.9219\n", "Tdata[1, 15] = 0.8925\n", "Tdata[1, 16] = 0.9243\n", "Tdata[1, 17] = 0.9010\n", "Tdata[1, 18] = 0.8974\n", "Tdata[1, 19] = 0.9066\n", " \n", "lengths[2] = 0.2500\n", "Tdata[2, 0] = 1.0283\n", "Tdata[2, 1] = 1.0112\n", "Tdata[2, 2] = 1.0018\n", "Tdata[2, 3] = 0.9988\n", "Tdata[2, 4] = 1.0060\n", "Tdata[2, 5] = 0.9832\n", "Tdata[2, 6] = 0.9920\n", "Tdata[2, 7] = 1.0337\n", "Tdata[2, 8] = 1.0353\n", "Tdata[2, 9] = 1.0129\n", "Tdata[2, 10] = 0.9599\n", "Tdata[2, 11] = 0.9925\n", "Tdata[2, 12] = 1.0186\n", "Tdata[2, 13] = 1.0311\n", "Tdata[2, 14] = 1.0189\n", "Tdata[2, 15] = 0.9973\n", "Tdata[2, 16] = 1.0422\n", "Tdata[2, 17] = 1.0100\n", "Tdata[2, 18] = 1.0355\n", "Tdata[2, 19] = 0.9847\n", " \n", "lengths[3] = 0.3000\n", "Tdata[3, 0] = 1.1011\n", "Tdata[3, 1] = 1.0804\n", "Tdata[3, 2] = 1.1011\n", "Tdata[3, 3] = 1.1370\n", "Tdata[3, 4] = 1.1124\n", "Tdata[3, 5] = 1.1406\n", "Tdata[3, 6] = 1.1440\n", "Tdata[3, 7] = 1.1132\n", "Tdata[3, 8] = 1.0770\n", "Tdata[3, 9] = 1.1217\n", "Tdata[3, 10] = 1.0653\n", "Tdata[3, 11] = 1.1078\n", "Tdata[3, 12] = 1.0762\n", "Tdata[3, 13] = 1.0891\n", "Tdata[3, 14] = 1.1007\n", "Tdata[3, 15] = 1.1104\n", "Tdata[3, 16] = 1.1148\n", "Tdata[3, 17] = 1.1340\n", "Tdata[3, 18] = 1.1263\n", "Tdata[3, 19] = 1.1383\n", " \n", "lengths[4] = 0.3500\n", "Tdata[4, 0] = 1.2023\n", "Tdata[4, 1] = 1.2088\n", "Tdata[4, 2] = 1.1898\n", "Tdata[4, 3] = 1.1578\n", "Tdata[4, 4] = 1.1580\n", "Tdata[4, 5] = 1.1691\n", "Tdata[4, 6] = 1.1709\n", "Tdata[4, 7] = 1.2005\n", "Tdata[4, 8] = 1.1617\n", "Tdata[4, 9] = 1.2121\n", "Tdata[4, 10] = 1.1792\n", "Tdata[4, 11] = 1.1794\n", "Tdata[4, 12] = 1.1899\n", "Tdata[4, 13] = 1.1773\n", "Tdata[4, 14] = 1.1880\n", "Tdata[4, 15] = 1.1941\n", "Tdata[4, 16] = 1.1906\n", "Tdata[4, 17] = 1.2101\n", "Tdata[4, 18] = 1.1644\n", "Tdata[4, 19] = 1.2042\n", " \n", "lengths[5] = 0.4000\n", "Tdata[5, 0] = 1.2866\n", "Tdata[5, 1] = 1.2700\n", "Tdata[5, 2] = 1.3042\n", "Tdata[5, 3] = 1.2549\n", "Tdata[5, 4] = 1.2391\n", "Tdata[5, 5] = 1.2739\n", "Tdata[5, 6] = 1.2845\n", "Tdata[5, 7] = 1.2883\n", "Tdata[5, 8] = 1.2333\n", "Tdata[5, 9] = 1.2882\n", "Tdata[5, 10] = 1.2963\n", "Tdata[5, 11] = 1.2802\n", "Tdata[5, 12] = 1.2337\n", "Tdata[5, 13] = 1.2428\n", "Tdata[5, 14] = 1.2219\n", "Tdata[5, 15] = 1.2674\n", "Tdata[5, 16] = 1.2369\n", "Tdata[5, 17] = 1.2852\n", "Tdata[5, 18] = 1.2316\n", "Tdata[5, 19] = 1.2530\n", " \n", "lengths[6] = 0.4500\n", "Tdata[6, 0] = 1.3222\n", "Tdata[6, 1] = 1.3378\n", "Tdata[6, 2] = 1.3694\n", "Tdata[6, 3] = 1.3523\n", "Tdata[6, 4] = 1.3741\n", "Tdata[6, 5] = 1.3856\n", "Tdata[6, 6] = 1.3410\n", "Tdata[6, 7] = 1.3486\n", "Tdata[6, 8] = 1.3504\n", "Tdata[6, 9] = 1.3355\n", "Tdata[6, 10] = 1.3103\n", "Tdata[6, 11] = 1.3260\n", "Tdata[6, 12] = 1.3859\n", "Tdata[6, 13] = 1.3772\n", "Tdata[6, 14] = 1.3549\n", "Tdata[6, 15] = 1.3051\n", "Tdata[6, 16] = 1.3625\n", "Tdata[6, 17] = 1.3229\n", "Tdata[6, 18] = 1.3649\n", "Tdata[6, 19] = 1.3589\n", " \n", "lengths[7] = 0.5000\n", "Tdata[7, 0] = 1.3949\n", "Tdata[7, 1] = 1.4244\n", "Tdata[7, 2] = 1.4159\n", "Tdata[7, 3] = 1.4001\n", "Tdata[7, 4] = 1.4074\n", "Tdata[7, 5] = 1.4798\n", "Tdata[7, 6] = 1.4318\n", "Tdata[7, 7] = 1.4391\n", "Tdata[7, 8] = 1.4176\n", "Tdata[7, 9] = 1.4158\n", "Tdata[7, 10] = 1.4333\n", "Tdata[7, 11] = 1.3937\n", "Tdata[7, 12] = 1.4353\n", "Tdata[7, 13] = 1.4086\n", "Tdata[7, 14] = 1.4065\n", "Tdata[7, 15] = 1.3917\n", "Tdata[7, 16] = 1.4250\n", "Tdata[7, 17] = 1.4085\n", "Tdata[7, 18] = 1.4167\n", "Tdata[7, 19] = 1.4566\n", " \n", "lengths[8] = 0.5500\n", "Tdata[8, 0] = 1.4738\n", "Tdata[8, 1] = 1.5150\n", "Tdata[8, 2] = 1.5008\n", "Tdata[8, 3] = 1.4711\n", "Tdata[8, 4] = 1.4941\n", "Tdata[8, 5] = 1.4934\n", "Tdata[8, 6] = 1.4623\n", "Tdata[8, 7] = 1.4486\n", "Tdata[8, 8] = 1.4486\n", "Tdata[8, 9] = 1.4887\n", "Tdata[8, 10] = 1.4855\n", "Tdata[8, 11] = 1.4628\n", "Tdata[8, 12] = 1.4511\n", "Tdata[8, 13] = 1.4383\n", "Tdata[8, 14] = 1.4956\n", "Tdata[8, 15] = 1.4849\n", "Tdata[8, 16] = 1.5078\n", "Tdata[8, 17] = 1.5277\n", "Tdata[8, 18] = 1.4606\n", "Tdata[8, 19] = 1.4962\n", " \n", "lengths[9] = 0.6000\n", "Tdata[9, 0] = 1.5181\n", "Tdata[9, 1] = 1.5678\n", "Tdata[9, 2] = 1.5354\n", "Tdata[9, 3] = 1.5145\n", "Tdata[9, 4] = 1.5495\n", "Tdata[9, 5] = 1.6030\n", "Tdata[9, 6] = 1.5545\n", "Tdata[9, 7] = 1.5595\n", "Tdata[9, 8] = 1.5181\n", "Tdata[9, 9] = 1.5540\n", "Tdata[9, 10] = 1.5795\n", "Tdata[9, 11] = 1.5783\n", "Tdata[9, 12] = 1.5721\n", "Tdata[9, 13] = 1.5580\n", "Tdata[9, 14] = 1.5535\n", "Tdata[9, 15] = 1.5230\n", "Tdata[9, 16] = 1.5730\n", "Tdata[9, 17] = 1.5690\n", "Tdata[9, 18] = 1.4998\n", "Tdata[9, 19] = 1.5381" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From these measurements, the value of $\\overline T$ for each set of 20 oscillations was determined." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " \n", "Raw length (m) period (s^2) data\n", " \n", "Data set 0 lengths 0.15\n", "Tdata\n", " [0.7876 0.7671 0.769 0.8073 0.7555 0.7742 0.7558 0.8065 0.7632 0.75\n", " 0.7797 0.8078 0.7855 0.787 0.7854 0.7876 0.7629 0.7256 0.7934 0.7944]\n", " \n", "Data set 1 lengths 0.2\n", "Tdata\n", " [0.8908 0.9079 0.9031 0.9162 0.9197 0.9046 0.9005 0.8523 0.8949 0.9214\n", " 0.8598 0.8714 0.8703 0.9229 0.9219 0.8925 0.9243 0.901 0.8974 0.9066]\n", " \n", "Data set 2 lengths 0.25\n", "Tdata\n", " [1.0283 1.0112 1.0018 0.9988 1.006 0.9832 0.992 1.0337 1.0353 1.0129\n", " 0.9599 0.9925 1.0186 1.0311 1.0189 0.9973 1.0422 1.01 1.0355 0.9847]\n", " \n", "Data set 3 lengths 0.3\n", "Tdata\n", " [1.1011 1.0804 1.1011 1.137 1.1124 1.1406 1.144 1.1132 1.077 1.1217\n", " 1.0653 1.1078 1.0762 1.0891 1.1007 1.1104 1.1148 1.134 1.1263 1.1383]\n", " \n", "Data set 4 lengths 0.35\n", "Tdata\n", " [1.2023 1.2088 1.1898 1.1578 1.158 1.1691 1.1709 1.2005 1.1617 1.2121\n", " 1.1792 1.1794 1.1899 1.1773 1.188 1.1941 1.1906 1.2101 1.1644 1.2042]\n", " \n", "Data set 5 lengths 0.4\n", "Tdata\n", " [1.2866 1.27 1.3042 1.2549 1.2391 1.2739 1.2845 1.2883 1.2333 1.2882\n", " 1.2963 1.2802 1.2337 1.2428 1.2219 1.2674 1.2369 1.2852 1.2316 1.253 ]\n", " \n", "Data set 6 lengths 0.45\n", "Tdata\n", " [1.3222 1.3378 1.3694 1.3523 1.3741 1.3856 1.341 1.3486 1.3504 1.3355\n", " 1.3103 1.326 1.3859 1.3772 1.3549 1.3051 1.3625 1.3229 1.3649 1.3589]\n", " \n", "Data set 7 lengths 0.5\n", "Tdata\n", " [1.3949 1.4244 1.4159 1.4001 1.4074 1.4798 1.4318 1.4391 1.4176 1.4158\n", " 1.4333 1.3937 1.4353 1.4086 1.4065 1.3917 1.425 1.4085 1.4167 1.4566]\n", " \n", "Data set 8 lengths 0.55\n", "Tdata\n", " [1.4738 1.515 1.5008 1.4711 1.4941 1.4934 1.4623 1.4486 1.4486 1.4887\n", " 1.4855 1.4628 1.4511 1.4383 1.4956 1.4849 1.5078 1.5277 1.4606 1.4962]\n", " \n", "Data set 9 lengths 0.6\n", "Tdata\n", " [1.5181 1.5678 1.5354 1.5145 1.5495 1.603 1.5545 1.5595 1.5181 1.554\n", " 1.5795 1.5783 1.5721 1.558 1.5535 1.523 1.573 1.569 1.4998 1.5381]\n", " \n", "Average periods and squared average periods \n", "Length (m)\tAverage period (s)\tAverage period squared (s^2)\n", "0.150\t\t0.777\t\t\t0.604\n", "0.200\t\t0.899\t\t\t0.808\n", "0.250\t\t1.010\t\t\t1.019\n", "0.300\t\t1.110\t\t\t1.231\n", "0.350\t\t1.185\t\t\t1.405\n", "0.400\t\t1.264\t\t\t1.597\n", "0.450\t\t1.349\t\t\t1.821\n", "0.500\t\t1.420\t\t\t2.017\n", "0.550\t\t1.480\t\t\t2.191\n", "0.600\t\t1.551\t\t\t2.405\n" ] } ], "source": [ "for n in range(0, nLen):\n", " Tmean[n] = np.cumsum(Tdata, 1)[n, nMeasPerLen - 1]/nMeasPerLen\n", "print(\" \")\n", "print(\"Raw length (m) period (s^2) data\")\n", "np.set_printoptions(precision = 4)\n", "for n in range(0, nLen):\n", " print(\" \")\n", " print(\"Data set\",n,\"lengths\",lengths[n])\n", " print(\"Tdata\\n\",Tdata[n, 0:nMeasPerLen])\n", "#\n", "print(\" \")\n", "print(\"Average periods and squared average periods \")\n", "print(\"Length (m)\\tAverage period (s)\\tAverage period squared (s^2)\")\n", "for n in range(0, nLen):\n", " print(\"{:2.3f}\\t\\t{:2.3f}\\t\\t\\t{:2.3f}\".format(lengths[n], Tmean[n], Tmean[n]**2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Data analysis\n", "\n", "The best way to display your data and extract the value of $g$ from it is to plot it on a graph. If appropriately presented, the data will lie on a straight line, which can be described by the equation $y = mx + c$. In this case, a little manipulation of your results is needed to get the data into the required form. We have:\n", "\n", "$$\n", "T = 2\\pi \\sqrt {\\frac{l}{g}} \\\\\n", "\\Rightarrow {T^2} = \\frac{4\\pi^2}{g} l\n", "$$\n", "\n", "Written in this way, we see that a plot of $l$ (on the $x$ axis) against $T^2$ (on the $y$ axis) will have a gradient $m = \\frac{4\\pi^2}{g}$ and we would expect the intercept to be $c = 0$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Exercise 3\n", "\n", "Plot your data using the cell below. Add a straight line to the plot and adjust the values of the gradient and the intercept until you get what looks like the best fit to your data." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfwAAAGECAYAAADTI5K/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAgPElEQVR4nO3dfbQkdX3n8fdnBoheJA8648MBZi4aPBETQDJiEnyAjbJgHognJkJuMEs0Nz6uedwYJ9Ek7uQka47JcUXxRll0t8XdBPAhIhATIxqCMijhUQ1BZhxhw/iEyGXFge/+UTWh53Jnphu6+z7U+3VOn+r61q+qv7fOnfncqq7uSlUhSZJWtzVL3YAkSRo/A1+SpA4w8CVJ6gADX5KkDjDwJUnqAANfkqQOMPAlSeoAA1/SopKcl6Tax3eS3JHkY0lemeTAIbZzYruNdePsV9K+GfiS9uWjwBOAaeBk4EPAHwKfSHLwEvYlaUgGvqR9+XZV/d+q+nJVXVNVbwZOBI4D/gtAkl9MclWSu9qzAH+V5NB22TTwsXZbO9sj/fPaZack+USSryf5WpJLkzxl0j+g1BUGvqShVNX1wCXAz7alg4A3AMcAPwmsA85vl32pb9xTac4WvKadPxj4C+B4mj8i7gQ+lOSgsf4AUkcdsNQNSFqRbgSeC1BV5/bVb0nycuCmJIdV1Y4kX2uX3VFVX9k9sKou6N9gkrOAb9L8AfDJsXYvdZBH+JIeigAFkOS4JB9Isi3JXcDWdsyGfW4geVKS9yb51yTfBP6N5v+kfa4n6aEx8CU9FEfRHM0fDFwKzANnAk8HTmnH7O/U/IeA9cCvAs8AngbsGmA9SQ+BgS9pKEl+kCbU/xr4AZr37F9XVZdX1eeAxy5Y5d52urZvG48BngL8cVV9tKpuAg7BtxmlsfEfl6R9+a4kj6c5OFgP/DjwOuBq4M+AKeDbwKuSnE0T4m9csI1tNKf/fyLJh4B7gK8DXwF+JcmXgEOBN9Ec4UsaA4/wJe3Lc4Hbge3A3wE/TfM5/GdX1d1VtRP4JeBnaC7kewPwG/0bqKovt/UtNO/Tv7Wq7gdeBBwNXA+cDfw+zR8PksYgVbXUPUiSpDHzCF+SpA4w8CVJ6gADX5KkDjDwJUnqAANfkqQOWNWfw1+3bl1NT08vdRuSJE3M1Vdf/ZWqWr+wvqoDf3p6mq1bt+5/oCRJq0SSbYvVPaUvSVIHGPiSJHWAgS9JUgcY+JIkdYCBL0lSBxj4kiR1gIEvSVIHGPiSJHWAgS9JUgcY+JIkTVqvB9PTsGZNM+31xv6Sq/qrdSVJWnZ6PZidhfn5Zn7btmYeYGZmbC/rEb4kSZO0efMDYb/b/HxTHyMDX5KkSdq+fbj6iBj4kiRN0oYNw9VHxMCXJGmStmyBqak9a1NTTX2MDHxJkiZpZgbm5mDjRkia6dzcWC/YA6/SlyRp8mZmxh7wC3mEL0lSBxj4kiR1gIEvSVIHGPiSJHWAgS9JUgcY+JIkdYCBL0lSBxj4kiR1gIEvSVIHGPiSJHWAgS9JUgcY+JIkdYCBL0lSBxj4kiR1gIEvSVIHGPiSJHWAgS9JUgcY+JIkdcBEAz/J4Uk+luSmJDckec0iY2aSXNs+rkhyTN+yW5Ncl+SaJFsn2bskSSvZARN+vV3Ab1bVZ5IcAlyd5G+r6sa+MV8EnlNVX09yKjAHPKNv+UlV9ZUJ9ixJ0oo30cCvqtuB29vndyW5CTgUuLFvzBV9q1wJHDbJHiVJWo2W7D38JNPA04BP7WPYS4CP9M0XcFmSq5PM7mW7s0m2Jtm6c+fOkfUrSdJKNulT+gAkeRRwAfBrVfXNvYw5iSbwn9lXPqGqbkvyWOBvk3yuqi7vX6+q5mjeBmDTpk01lh9AkqQVZuJH+EkOpAn7XlVduJcxRwPvBE6rqq/urlfVbe30DuAi4PjxdyxJ0so36av0A7wLuKmq3ryXMRuAC4Ezq+oLffWD2wv9SHIwcDJw/fi7liRp5Zv0Kf0TgDOB65Jc09ZeB2wAqKpzgNcDjwHe1vx9wK6q2gQ8DriorR0AvLeqLplo95IkrVCTvkr/k0D2M+alwEsXqd8CHPPgNSRJ0v74TXuSJHWAgS9JWtl6PZiehjVrmmmvt9QdLUtL8rE8SZJGoteD2VmYn2/mt21r5gFmZpaur2XII3xJ0sq1efMDYb/b/HxT1x4MfEnSyrV9+3D1DjPwJUkr14YNw9U7zMCXJK1cW7bA1NSetamppq49GPiSpJVrZgbm5mDjRkia6dycF+wtwqv0JUkr28yMAT8Aj/AlSeoAA1+SpA4w8CVJ6gADX5KkDjDwJUnqAANfkqQOMPAlSeoAA1+SpA4w8CVJ6gADX5KkDjDwJUnqAANfkqQOMPAlSeoAA1+SpA4w8CVJ6gADX5KkDjDwJUnqAANfkqQOMPAlSeoAA1+SpA4w8CVJD+j1YHoa1qxppr3eUnekETlgqRuQJC0TvR7MzsL8fDO/bVszDzAzs3R9aSQ8wpckNTZvfiDsd5ufb+pa8SYa+EkOT/KxJDcluSHJaxYZkyRvSXJzkmuTHNe37JQkn2+XvXaSvUvSqrd9+3B1rSiTPsLfBfxmVT0F+BHglUmOWjDmVODI9jELvB0gyVrg7Hb5UcAZi6wrSXqoNmwYrq4VZaKBX1W3V9Vn2ud3ATcBhy4YdhrwnmpcCXxvkicAxwM3V9UtVXUv8L52rCRpFLZsgampPWtTU01dK96SvYefZBp4GvCpBYsOBb7UN7+jre2tvnC7s0m2Jtm6c+fOkfYsSavazAzMzcHGjZA007k5L9hbJZbkKv0kjwIuAH6tqr65cPEiq9Q+6nsWquaAOYBNmzY9aLkkaR9mZgz4VWrigZ/kQJqw71XVhYsM2QEc3jd/GHAbcNBe6pIkaT8mfZV+gHcBN1XVm/cy7IPAi9ur9X8EuLOqbgeuAo5MckSSg4DT27GSJGk/Jn2EfwJwJnBdkmva2uuADQBVdQ5wMfB84GZgHjirXbYryauAS4G1wLlVdcNEu5ckaYWaaOBX1SdZ/L34/jEFvHIvyy6m+YNAkiQNwW/akySpAwx8SZI6wMCXJKkDDHxJkjrAwJckqQMMfEmSOsDAlySpAwx8SZI6wMCXJKkDDHxJkjrAwJckqQMMfEmSOsDAlySpAwx8SZI6wMCXJKkDDHxJkjrAwJckqQMMfEmSOsDAlySpAwx8SZI6wMCXJKkDDHxJGpdeD6anYc2aZtrrLXVH6rADlroBSVqVej2YnYX5+WZ+27ZmHmBmZun6Umd5hC9J47B58wNhv9v8fFOXloCBL0njsH37cHVpzAx8SRqHDRuGq0tjZuBL0jhs2QJTU3vWpqaaurQEDHxJGoeZGZibg40bIWmmc3NesKcl41X6kjQuMzMGvJYNj/AlSeoAA1+SpA4w8CVJ6oCJvoef5FzgJ4E7quoHF1n+28DuN7wOAJ4CrK+qryW5FbgLuA/YVVWbJtO1JEkr36SP8M8DTtnbwqp6U1UdW1XHAr8LfLyqvtY35KR2uWEvSdIQJhr4VXU58LX9DmycAZw/xnYkSeqMZfkefpIpmjMBF/SVC7gsydVJZvex7mySrUm27ty5c9ytSpK0IizLwAd+CvjHBafzT6iq44BTgVcmefZiK1bVXFVtqqpN69evn0SvkiQte8s18E9nwen8qrqtnd4BXAQcvwR9SZK0Ii27wE/yPcBzgA/01Q5Ocsju58DJwPVL06EkSSvPpD+Wdz5wIrAuyQ7gDcCBAFV1TjvsBcBlVXV336qPAy5KAk3P762qSybVtyRJK91EA7+qzhhgzHk0H9/rr90CHDOeriRJWv2W3Sl9SZI0ega+JEkdYOBLktQBBr4kSR1g4EuS1AEGviRJHWDgS5LUAQa+JEkd8JACv/2q27WjbkaSJI3HQIGfZE2SX0jy4SR3AJ8Dbk9yQ5I3JTlyvG1KkqSHY9Aj/I8BTwJ+F3h8VR1eVY8FngVcCfxJkl8cU4+SJOlhGvS79J9bVd9ZWGzvV38BcEGSA0famSRJGpmBjvAXC/uHMkaSJC2N/QZ+kucl+cskx7bzs2PvSpIkjdQgp/RfAZwF/F6SRwPHjrUjSZI0coOc0t9ZVd+oqt8CTgaePuaeJEnSiA0S+B/e/aSqXgu8Z3ztSJKkcdhv4FfVB3Y/T/JzwHnt899LcmGS48bXniS1ej2YnoY1a5ppr7fUHUkryrDftPf7VXVXkmcC/xF4N/D20bclSX16PZidhW3boKqZzs4a+tIQhg38+9rpTwBvb4/+DxptS5K0wObNMD+/Z21+vqlLGsiwgf/lJO8Afh64OMl3PYRtSNJwtm8fri7pQYYN658HLgVOqapvAI8GfnvUTUnSHjZsGK4u6UEGvXlOAKpqvqourKp/aedvr6rL+sdI0sht2QJTU3vWpqaauqSBDHzznCSvTrLHn9NJDkryH5K8G/il0bcnScDMDMzNwcaNkDTTubmmLmkgqar9D0oeAfwyMAMcAXwDeASwFrgMOLuqrhlblw/Rpk2bauvWrUvdhiRJE5Pk6qratLA+0N3yqur/AW8D3tbeFW8dcE/7Pr4kSVrmBr097r9r74p3+xh6kSRJY+JH6iRJ6gADX5KkDjDwJUnqgP0GfpLnJfnLJMe287Nj70qSJI3UIBftvQI4C/i9JI8Gjh1rR5IkaeQGOaW/s6q+UVW/BZwMPP2hvliSc5PckeT6vSw/McmdSa5pH6/vW3ZKks8nuTnJax9qD5IkddEggf/h3U+q6rXAex7G650HnLKfMZ+oqmPbxx8BJFkLnA2cChwFnJHkqIfRhyRJnbLfU/pV9YEkv9FfWzhfVW8e5MWq6vIk00N12DgeuLmqbmlf/33AacCND2FbkiR1zqBX6R/SPjYBLwcObR8vozniHqUfTfLPST6S5Klt7VDgS31jdrQ1SZI0gEG/WvcPAZJcBhxXVXe1838A/NUI+/kMsLGqvpXk+cD7gSOBxe7Et+hNANpPEcwCbPDWmZIkAcN/Dn8DcG/f/L3A9KiaqapvVtW32ucXAwcmWUdzRH9439DDgNv2so25qtpUVZvWr18/qtYkSVrRhv0u/f8JfDrJRe38zwDvHlUzSR4P/FtVVZLjaf4g+SrN3fmOTHIE8GXgdOAXRvW6kiStdkMFflVtSfIR4Fk0p9TPqqrPDrp+kvOBE4F1SXYAbwAObLd9DvBC4OVJdgH3AKdXc//eXUleBVxKc0vec6vqhmF6lySpy4a+Wx5wH3A/TeDfP8yKVXXGfpa/FXjrXpZdDFw8zOtJkqTGUO/hJ3kN0APWAY8F/leSV4+jMUmSNDrDHuG/BHhGVd0NkORPgX8C/vuoG5MkSaMz7FX6oTmlv9t9LP6ROUmStIwMe4T/P4BPLbhK/10j7UiSJI3csFfpvznJx4ETaI7sh7pKX5IkLY2hr9KvqquBq8fQiyRJGpOBAj/JJ6vqmUnuYs+vtA1QVfXdY+lOkiSNxEAX7bVhH+CpVfXdfY9DDHtpBer1YHoa1qxppr3eUnckacwGvkq//ca7i/Y7UNLy1uvB7Cxs2wZVzXR21tCXVrlhP5Z3ZZKnj6UTSZOxeTPMz+9Zm59v6pJWrWEv2jsJeFmSW4G7eeA9/KNH3ZikMdm+fbi6pFVh2MA/dSxdSJqcDRua0/iL1SWtWsOe0t9Oc6e8X6qqbTRX7D9u5F1JGp8tW2Bqas/a1FRTl7RqDRv4bwN+FNh917u7gLNH2pGk8ZqZgbk52LgRkmY6N9fUJa1aw57Sf0ZVHZfkswBV9fUkB42hL0njNDNjwEsdM+wR/neSrKX98p0k64H7R96VJEkaqWED/y00n8V/bJItwCeBPx55V5IkaaSGvXlOL8nVwI/TfCTvZ6rqprF0JkmSRmbQ79J/BPAy4PuB64B3VNWucTYmSZJGZ9BT+u8GNtGE/anAn42tI0mSNHKDntI/qqp+CCDJu4BPj68lSZI0aoMe4X9n9xNP5UuStPIMeoR/TJJvts8DPLKd3/1d+t4iV5KkZWygwK+qteNuRJIkjc+wn8OXJEkrkIEvSVIHGPiSJHWAgS9JUgcY+JIkdYCBL0lSBxj4kiR1gIEvSVIHTDTwk5yb5I4k1+9l+UySa9vHFUmO6Vt2a5LrklyTZOvkupYkaeWb9BH+ecAp+1j+ReA5VXU08EZgbsHyk6rq2KraNKb+JElalQb9Lv2RqKrLk0zvY/kVfbNXAoeNvSlJkjpgOb+H/xLgI33zBVyW5Ooks0vUkyRJK9KyDPwkJ9EE/u/0lU+oquOAU4FXJnn2XtadTbI1ydadO3dOoFt1Xq8H09OwZk0z7fWWuiNJepBlF/hJjgbeCZxWVV/dXa+q29rpHcBFwPGLrV9Vc1W1qao2rV+/fhItq8t6PZidhW3boKqZzs4a+pKWnWUV+Ek2ABcCZ1bVF/rqByc5ZPdz4GRg0Sv9pYnavBnm5/eszc83dUlaRiZ60V6S84ETgXVJdgBvAA4EqKpzgNcDjwHelgRgV3tF/uOAi9raAcB7q+qSSfYuLWr79uHqkrREJn2V/hn7Wf5S4KWL1G8BjnnwGtIS27ChOY2/WF2SlpFldUpfWnG2bIGpqT1rU1NNXZKWEQNfejhmZmBuDjZuhKSZzs01dUlaRiZ6Sl9alWZmDHhJy55H+JIkdYCBL0lSBxj4kiR1gIEvSVIHGPiSJHWAgS9JUgcY+JIkdYCBL0lSBxj4kiR1gIEvSVIHGPiSJHWAgS9JUgcY+JIkdYCBL0lSBxj4kiR1gIEvSVIHGPiSJHWAgS9JUgcY+JIkdYCBL0lSBxj4kiR1gIEvSVIHGPiSJHWAgS9JUgcY+JIkdYCBL0lSBxj4kiR1gIGv5aXXg+lpWLOmmfZ6S92RJK0KByx1A9K/6/Vgdhbm55v5bduaeYCZmaXrS5JWgYke4Sc5N8kdSa7fy/IkeUuSm5Ncm+S4vmWnJPl8u+y1k+taE7N58wNhv9v8fFOXJD0skz6lfx5wyj6Wnwoc2T5mgbcDJFkLnN0uPwo4I8lRY+1Uk7d9+3B1SdLAJhr4VXU58LV9DDkNeE81rgS+N8kTgOOBm6vqlqq6F3hfO1aryYYNw9UlSQNbbhftHQp8qW9+R1vbW12ryZYtMDW1Z21qqqlLkh6W5Rb4WaRW+6g/eAPJbJKtSbbu3LlzpM1pzGZmYG4ONm6EpJnOzXnBniSNwHK7Sn8HcHjf/GHAbcBBe6k/SFXNAXMAmzZtWvSPAi1jMzMGvCSNwXI7wv8g8OL2av0fAe6sqtuBq4AjkxyR5CDg9HasJEkawESP8JOcD5wIrEuyA3gDcCBAVZ0DXAw8H7gZmAfOapftSvIq4FJgLXBuVd0wyd4lSVrJJhr4VXXGfpYX8Mq9LLuY5g8CSZI0pOV2Sl+SJI2BgS9JUgcY+JIkdYCBL0lSBxj4kiR1gIEvSVIHGPiSJHWAgS9JUgcY+JIkdYCBL0lSBxj4kiR1gIEvSVIHGPiSJHWAgS9JUgcY+JIkdYCBL0lSBxj4kiR1gIEvSVIHGPiSJHWAgS9JUgcY+JIkdYCBL0lSBxj4kiR1gIG/mvV6MD0Na9Y0015vqTuSJC2RA5a6AY1JrwezszA/38xv29bMA8zMLF1fkqQl4RH+arV58wNhv9v8fFOXJHWOgb9abd8+XF2StKoZ+KvVhg3D1SVJq5qBv1pt2QJTU3vWpqaauiSpcwz81WpmBubmYONGSJrp3JwX7ElSR3mV/mo2M2PAS5IAj/AlSeoEA1+SpA6YeOAnOSXJ55PcnOS1iyz/7STXtI/rk9yX5NHtsluTXNcu2zrp3iVJWqkm+h5+krXA2cDzgB3AVUk+WFU37h5TVW8C3tSO/yng16vqa32bOamqvjLBtiVJWvEmfYR/PHBzVd1SVfcC7wNO28f4M4DzJ9KZJEmr2KQD/1DgS33zO9ragySZAk4BLugrF3BZkquTzO5lvdkkW5Ns3blz54jaliRpZZt04GeRWu1l7E8B/7jgdP4JVXUccCrwyiTPftDGquaqalNVbVq/fv3D71iSpFVg0oG/Azi8b/4w4La9jD2dBafzq+q2dnoHcBHNWwSSJGk/Jh34VwFHJjkiyUE0of7BhYOSfA/wHOADfbWDkxyy+zlwMnD9RLqWJGmFm+hV+lW1K8mrgEuBtcC5VXVDkpe1y89ph74AuKyq7u5b/XHARUmg6fu9VXXJ5LqXJGnlStXe3kJf+TZt2lRbt/pxfUlSdyS5uqo2Laz7TXuSJHWAgS9JUgcY+JIkdYCBL0lSBxj4kiR1gIEvSVIHGPiSJHWAgS9JUgcY+JIkdYCBL0lSBxj4kiR1gIEvSVIHGPiSJHWAgS9JUgcY+JIkdYCBP6heD6anYc2aZtrrLXVHkiQN7IClbmBF6PVgdhbm55v5bduaeYCZmaXrS5KkAXmEP4jNmx8I+93m55u6JEkrgIE/iO3bh6tLkrTMGPiD2LBhuLokScuMgT+ILVtgamrP2tRUU5ckaQUw8AcxMwNzc7BxIyTNdG7OC/YkSSuGV+kPambGgJckrVge4UuS1AEGviRJHWDgS5LUAQa+JEkdYOBLktQBBr4kSR1g4EuS1AEGviRJHWDgS5LUAQa+JEkdkKpa6h7GJslOYNtS97EMrAO+stRNdID7eTLcz5Phfp6MceznjVW1fmFxVQe+Gkm2VtWmpe5jtXM/T4b7eTLcz5Mxyf3sKX1JkjrAwJckqQMM/G6YW+oGOsL9PBnu58lwP0/GxPaz7+FLktQBHuFLktQBBv4qkeSUJJ9PcnOS1y6yfCbJte3jiiTHLEWfq8EA+/q0dj9fk2RrkmcuRZ8r3f72c9+4pye5L8kLJ9nfajHA7/OJSe5sf5+vSfL6pehzpRvk97nd19ckuSHJx0feRFX5WOEPYC3wr8ATgYOAfwaOWjDmx4Dva5+fCnxqqfteiY8B9/WjeODtsqOBzy113yvtMch+7hv398DFwAuXuu+V9hjw9/lE4G+WuteV/BhwP38vcCOwoZ1/7Kj78Ah/dTgeuLmqbqmqe4H3Aaf1D6iqK6rq6+3slcBhE+5xtRhkX3+r2n+xwMGAF8oMb7/7ufVq4ALgjkk2t4oMup/18Ayyn38BuLCqtgNU1ch/pw381eFQ4Et98zva2t68BPjIWDtavQba10lekORzwIeBX55Qb6vJfvdzkkOBFwDnTLCv1WbQ/zt+NMk/J/lIkqdOprVVZZD9/GTg+5L8Q5Krk7x41E0cMOoNaklkkdqiR5VJTqIJfN9XfmgG2tdVdRFwUZJnA28EnjvuxlaZQfbzXwC/U1X3JYsN1wAG2c+fofmq1m8leT7wfuDIcTe2ygyynw8Afhj4ceCRwD8lubKqvjCqJgz81WEHcHjf/GHAbQsHJTkaeCdwalV9dUK9rTYD7evdquryJE9Ksq6q/F7ywQ2ynzcB72vDfh3w/CS7qur9E+lwddjvfq6qb/Y9vzjJ2/x9Htogv887gK9U1d3A3UkuB44BRhb4ntJfHa4CjkxyRJKDgNOBD/YPSLIBuBA4c5R/MXbQIPv6+9OmUJLjaC7S8Q+s4ex3P1fVEVU1XVXTwF8DrzDshzbI7/Pj+36fj6fJDX+fh7Pf/Qx8AHhWkgOSTAHPAG4aZRMe4a8CVbUryauAS2muBj23qm5I8rJ2+TnA64HHAG9r/+3uKm+MMbQB9/XPAi9O8h3gHuBFfRfxaQAD7mc9TAPu5xcCL0+yi+b3+XR/n4czyH6uqpuSXAJcC9wPvLOqrh9lH37TniRJHeApfUmSOsDAlySpAwx8SZI6wMCXJKkDDHxJkjrAwJeWifaOb9ckuT7JX7WfxX242/xPSd66nzHTSUb68Z8xb/fEJD/WN3/eIHfKS/LIJB9PsnaI13pVkrMeaq/ScmLgS8vHPVV1bFX9IHAv8LKlbmiZOpHm7o/D+mWam5PcN8Q65wL/+SG8lrTsGPjS8vQJ4PuTHJzk3CRXJflsktPg34/cL0xySZJ/SfLfdq+Y5KwkX2jvp31CX32PI+Ek31r4ogvPCCT5myQn7h6f5E/bG3t8NMnx7Y0+bkny0/v6YZKsTfKm9ue4NsmvtvUT2238dZLPJen1favb89vaJ5O8pe1lmuYPoV9vz4Y8q32JZye5ou1lb0f7MzTfZrb7dT+e5P+0++pPkswk+XSS65I8CaCq5oFb22+Yk1Y0A19aZpIcAJwKXAdsBv6+qp4OnAS8KcnB7dBjgRcBPwS8KMnhSZ4A/CFN0D8POGqErR0M/ENV/TBwF/Bf29d4AfBH+1n3JcCd7c/xdOBXkhzRLnsa8Gttr08ETkjyCOAdNPd9eCawHqCqbqW5O96ft2dDPtFu4wk0N4T6SeBPFr54+3WmT2zX3+0Y4DU0++9M4MlVdTzN/SZe3TduK/AspBXOr9aVlo9HJrmmff4J4F3AFcBPJ/mttv4IYEP7/O+q6k6AJDcCG2luIvMPVbWzrf9vmttujsK9wCXt8+uAb1fVd5JcB0zvZ92TgaP7jr6/h+aOa/cCn66qHW2/17Tb+hZwS1V9sR1/PjC7j+2/v6ruB25M8rhFlq8DvrGgdlVV3d6+7r8Cl/X9bCf1jbsD+IH9/HzSsmfgS8vHPVV1bH+hPb39s1X1+QX1ZwDf7ivdxwP/nvf2fdm7aM/qtds9aF9jWo/oe/6dvu9Qv3/361fV/e1ZiX0J8OqqunTBz3HiXn6OYe9327+Nxda9hz1/loXr3N83fz97/t/4iHZ9aUXzlL60vF0KvLrvfe2n7Wf8p4ATkzwmyYHAz/Utu5XmftsApwEHLrL+rcCxSdYkORwY1XvXl9LcgOVAgCRP7ntrYjGfA57YvmcPzVsXu90FHDLMi1fV14G17VsFw3oyMPJPG0iTZuBLy9sbaYL52vYjbm/c1+D2FPUfAP8EfBT4TN/ivwSek+TTNLfevHuRTfwj8EWa09p/tmD9h+OdwI3AZ9qf4x3s4wxjVd0DvAK4JMkngX8D7mwXfwh4wYKL9gZxGc37/MM6gWZfSiuad8uTtCwleVRVfas9u3E28C9V9ecPY3tPA36jqs4c5zrScuURvqTl6lfai/huoLnI7x0PZ2NV9VngY8N88Q7NxX6//3BeV1ouPMKXJKkDPMKXJKkDDHxJkjrAwJckqQMMfEmSOsDAlySpAwx8SZI64P8DMJsp3IKa27QAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# \n", "#\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "#\n", "plt.figure(figsize = (8, 6))\n", "#\n", "# Plot the data\n", "plt.plot(lengths, Tmean**2, marker = 'o', color = 'r', linestyle = '')\n", "plt.title('Data', fontsize = 14)\n", "plt.xlabel('Pendulum length (m)')\n", "plt.ylabel('Period$^2$ (s$^2$)')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfYAAAGECAYAAADEAQJ2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA3G0lEQVR4nO3dd5gUVdrG4d8LgghmwYAIY1w/dSUNYEDFnFDXsCiSBAQUUVR0DZgVF8MaWAOishgGMWAGBZUsOQyYVtcEYgIVQRgkvt8fp9BxJPRAd1eH576uuaa7q7r7rWmYZ86pOueYuyMiIiK5oULcBYiIiEjyKNhFRERyiIJdREQkhyjYRUREcoiCXUREJIco2EVERHKIgl1EfmNmbmZnbWCfAWb2ejlft1n02tU3rcKNY2bvm9lNG9ins5nNMbPVZnaTmZ1nZovTVKJI0ijYRZIgCjuPvlaY2TwzG2lmF5lZpXK+VpwhuAvwWlRHQVRHYQx1pJWZbQc8CNwF7ArcDTwL7FFqn5vM7P14KhRJnIJdJHneJgRjAXAcISBvBsaaWbUY60qYu3/n7sviriMGdYDNgNfd/Vt3X+zuS919XtyFiZSXgl0keZZFwfi1uxe7+z1AM6AB8I81O5lZazObYma/RC37581s12hbATAy2nV+1GIeEG07wczGmtkCM/vJzIaZ2f+tqxgz+7/o+TtH96ua2XIze6PUPp3M7H+l7pfuiv8i+j4lenxUmdfvbmZfR/X8x8yqJvAzqmtmk8ysxMymmlmDUq+3g5k9Y2ZzzWypmX1gZu3LvOcoM3vIzG43sx+in9/dZlah1D47mtkr0WvMNrMO6yvIzM4DZkR3P4+OtaB0V3y0z43A/qV6Zs5L4HhF0k7BLpJC7v4+8CZwZqmHKxNCoi7QHKgOPBNt+6rUvvsTegC6R/erAfcBjQl/MCwEXjOzyut474+A76N9AQ6NntPUzDaLHmsGjFpH+Y2j7ydEdZxRatthwAHAMcDZwOml6lyffwJXE/7Y+REoMjOLtlUBphN+JvsD9wOPmNnRZV6jFbASOAToBlwa1bDGAGCvqLa/AW0JvSjr8mx0jBCOeRfC51B2n38BH0fbd4keE8k4m214FxHZRB8SQgYAd+9fatvnZnYh8JGZ1XL3uWb2U7Rtnrv/UOp5g0u/aNSaXUQIo3HreO/RwJHAIEKIvwCcCDQCJgBHAFet47nzo+8/uvt3ZbYtAi5095VR7c8DRxOCe32ud/eRUf23RHXvCsx1968J57jX6GdmRwEtgXdKPf6hu98Q3f7EzDpF7/2Mme0THV9Td383ep92wOfrKsjdl5rZj2uOec2x/v73xm/7LAZWruVnIZJR1GIXST0DflttycwaRF3Fs83sF2BqtKn2el/EbE8zG2hmn5nZIkJrvMIGnjeK31vszQjd/KOBZma2NyFUR5X3gAjhurLU/W+AHRN43qwyz2HN88ysopn1NLNZZvZjFKRn8Ofjm1Xmfun3/j9gNTB5zUZ3n13qvURynlrsIqm3H1GLMbqIbhjhQrs2wDxCV/xYQhf9+rwGfA10ib6vJPQGrO95o4CHohAvjO5XI7SCfwA+jVrK5bWizH0nsYZC6eet+WNnzfOuAHoQuvTfAxYDt/PnPxjW996GSJ5Ti10khczsAML52xeih/YlBPm17j7G3f/Ln4NrefS9YqnX2YHQGr3d3d+Ozp9vxQb+OC91nr0nIcTnEVrthwLHsv7W+p/qSLGmwGvu/pS7FwOfAfuU8zU+Ivxea7TmATOrDdRMQn3LSd/PQmSjKdhFkmdzM9vZzGqaWV0zu5wQnNMI46IB5gDLgG5mtoeZnQzcWuZ1ZhNaoSebWQ0z2xJYQGhhdzKzvczsCKAvodW+IaOB1kRX27v7l4Tz52ew/mCfBywFjjezncxsmwTea1N8AhxtZk3NbF/gAWD38ryAu39MuFjxETM72MzqES6mW5qE+r4E6kSnUqqb2eZJeE2RpFOwiyTPMcC3hPB+BziVMI79cHdfAuDu84F2hKu1PyRcHX956ReJusZvBHoRWtsPuPtqwpXfBwLvEyZTuZ7wR8KGjCS0NEeVemzUWh77g+gc+iXA+YRz1K8k8F6b4jbCufE3gDHAEqBoI17nPMJQvRGE0xcDCaG8qQYDQwmf7XzC6QyRjGPuvuG9REREJCuoxS4iIpJDFOwiIiI5RMEuIiKSQxTsIiIiOUTBLiIikkNyYua56tWre0FBQdxliIiIpMW0adN+cPcaa9uWE8FeUFDA1KlTN7yjiIhIDjCz2evapq54ERGRHKJgFxERySEKdhERkRyiYBcREckhCnYREZEcomAXERHJIQp2ERGRHKJgFxERySEKdhERkRyiYBcREckhCnYREZEcomAXERHJIQp2ERGRHKJgFxERySEKdhERkRyiYBcREckhCnYREZEcomAXERHJIQp2ERGRHKJgFxERySEKdhERkRyiYBcREckhCnYREZEcomAXERHJIQp2ERGRdfj5Z3jrrbirKB8Fu4iIyFpMmAD16sGZZ8KCBXFXkzgFu4iISCmrVsHtt8Nhh4EZDB8O220Xd1WJ2yzuAkRERDLFN99AmzYwYgScfTY88ghss03cVZWPgl1ERAQYMgTOOw9KSuDxx6F9+9BizzbqihcRkby2bBlcdhk0bw41a8LUqdChQ3aGOqjFLiIieeyTT+Ccc2DGDLj4YrjzTqhSJe6qNo2CXURE8o47PPEEdOsGm28Or7wCp54ad1XJoa54ERHJK4sWQevW4Rx6YSHMnJk7oQ4KdhERySNTpkD9+jBoENxyC7zzDtSqFXdVyaVgFxGRnLd6Ndx1FxxyCKxcCWPGwPXXQ8WKcVeWfDrHLiIiOe3776FdOxg2DM44Ax57LLsmnCkvtdhFRCRnDRsGBx4Io0dD377wwgu5HeqgYBcRkRy0fDlceSWccALUqBHOrXfpkr1j08sjrcFuZruZ2Ugz+8jMPjCz7mvZp5mZLTSz4ujrhnTWKCIi2e2zz6BpU7j77hDmkyfDAQfEXVX6pPsc+0qgh7tPN7OtgGlm9pa7f1hmv7Hu3jzNtYmISJYrKoILLwwXxb3wQliZLd+ktcXu7t+6+/To9i/AR8Cu6axBRERyz+LFYZ731q3DOfXi4vwMdYjxHLuZFQD1gUlr2Xywmc00szfMbP91PL+zmU01s6nz589PZakiIpLBZsyAhg3hySfDELZRo6BOnbirik8swW5mWwKDgUvdfVGZzdOBOu5eF/g38PLaXsPd+7l7obsX1qhRI6X1iohI5nGH++6Dgw6CJUvCUqu33AKb5flA7rQHu5lVIoR6kbu/WHa7uy9y98XR7aFAJTOrnuYyRUQkg82fH1Zju+wyOP740PXerFncVWWGdF8Vb8DjwEfufs869tk52g8za0yo8cf0VSkiIplsxAioWxfefhv69AkLuFRX8+836e6wOBRoA7xnZsXRY9cCtQHcvS9wFnChma0ElgLnuLunuU4REckwK1bAjTdC796wzz7wxhsh4OWP0hrs7j4OWO/0AO7+APBAeioSEZFs8OWX0LIlTJwIHTvC/fdDtWpxV5WZ8vwSAxERyXTPPw+dOoWL5QYNgrPPjruizKYpZUVEJCMtWRICvUUL2HffcIGcQn3DFOwiIhK/oiIoKIAKFaCggFn/HEJhITz+OFx9NYwdC7vvHneR2UFd8SIiEq+iIujcGUpKcOCh2SfT49qj2W6bpQwfvgXHHBN3gdlFLXYREYlXz55QUsKPbM/pvEQ3HuQoRjBzq6YK9Y2gFruIiMRrzhzGcBitKOJ7duIeLqM791Ph67gLy05qsYuISGxWroSbtr6HIxlJFX5lAgdzGfdRAYfateMuLysp2EVEJBZffQVHHQU3L7yU1hUHMZ0GNGR62Fi1KvTqFW+BWUrBLiIiaffSS2HWuBkz4Kmn4IknYKs6O4BZWJqtXz9o1SruMrOSzrGLiEjaLF0KPXrAww+HpVYHDYK99gJopSBPErXYRUQkLT74ABo3DqHeoweMH78m1CWZ1GIXEZGUcodHH4VLL4WttgqLt5xwQtxV5S612EVEJGUWLAhTwnbpAk2bwsyZCvVUU7CLiEhKjB8P9erByy/DHXfAm2/CzjvHXVXuU7CLiEhSrVoVRqodfjhUrAjjxsE//hGmgZfU0zl2ERFJmq+/htatYdSosH76ww/DNtvEXVV+UbCLiEhSvPYatG8fhrT17w/nnReGpUt6qWNEREQ2ya+/QvfucOqpsNtuMH16CHiFejwU7CIistE+/hgOPhj69AnhPnEi/OUvcVeV39QVLyIi5eYOAwZAt26wxRahG75587irElCLXUREymnhwjD7a4cOYSa5mTMV6plEwS4iIgmbNAnq14fnnoPbboO334Zdd427KilNwS4iIhu0enWYZKZp0zBOfcwY6NkzjFOXzKJz7CIisl7ffQdt2oTW+VlnhRVVt9su7qpkXRTsIiKyTm++CW3bwuLFIdDPP1/D2DKduuJFRORPli+HK66AE0+EnXaCqVOhUyeFejZQsIuISFBUBAUF/M/24ZCtZvGvf0HXrjB5Muy3X9zFSaIU7CIiEkK9c2eemn0YDZjG58tr8WLlc3jwkCK22CLu4qQ8FOwiIsIv19xO25KHactT1KOYYupx+vJnw6XvklV08ZyISJ6bNg3O+eplPmcPbuQmruM2NmNV2DhnTrzFSbmpxS4ikqdWr4Z77glzvf9asRojOIqbuPn3UAeoXTu+AmWjKNhFRPLQvHlhGtgePeCkk6D4gXc5ourUP+5UtSr06hVPgbLRFOwiInnm7behbl0YMQIefBBeegl2uODvYaB6nTphTFudOuF+q1ZxlyvlpHPsIiJ5YsUKuP56uPNO2HdfGDYMDjyw1A6tWinIc4CCXUQkD3zxBbRsGRZx6dQJ7r0XqlWLuypJBQW7iEiOGzQIunQJPezPPgstWsRdkaSSzrGLiOSoJUugY8fQUt9vPyguVqjnAwW7iEgOKi6Ghg3hP/+Ba68Ny6wWFMRdlaSDgl1EJIe4w7//DU2awKJF4Qr4Xr2gUqW4K5N00Tl2EZEc8cMP0KEDvPYanHxyaK3XqBF3VZJuarGLiOSAUaPC2PRhw+C++0K4K9Tzk4JdRCSLrVwZxqYfdVQYvjZhAnTvrnXT85m64kVEstTs2WE+mXffhfPOC+fWt9wy7qokbgp2EZEsNHgwnH8+rFoVllI/99y4K5JMoa54EZEssnQpXHABnHUW7L03zJihUJc/UrCLiGSJ99+HRo3gkUfgyith3DjYc8+4q5JMo2AXEckkRUVhJpkKFcL3oiLcoW/fEOrz54cr3++8EypXjrtYyURpDXYz283MRprZR2b2gZl1X8s+ZmZ9zOxTM5tlZg3SWaOISGyKiqBz53BVnDvMns1Pna7irCZzuPBCOPxwmDULjjsu7kIlk6W7xb4S6OHu/wccBFxkZvuV2edEYO/oqzPwcHpLFBGJSc+eUFLy291xHEq9peN5dcou3HUXvPEG7LRTjPVJVkhrsLv7t+4+Pbr9C/ARsGuZ3U4DnvRgIrCtme2SzjpFRGIxZw4Aq6jALVzPEYymMssZz6FccUXonRfZkNj+mZhZAVAfmFRm067AV6Xuz+XP4S8ikntq12Yuu3I073Ajt9CSZ5hOAxrVmRd3ZZJFYgl2M9sSGAxc6u6Lym5ey1N8La/R2cymmtnU+fPnp6JMEZG0euWMJ6jLTKZSyADa8RRt2LrqqrCKi0iC0h7sZlaJEOpF7v7iWnaZC+xW6n4t4JuyO7l7P3cvdPfCGpoQWUSy2K+/wsUXw9/uPYI6dYzpNU+hnT2F1akD/fqF6eVEEpTWmefMzIDHgY/c/Z517PYq0M3MBgFNgIXu/m26ahQRSaePPoJzzglXu196KfTuvT2bbz4y7rIki6V7StlDgTbAe2ZWHD12LVAbwN37AkOBk4BPgRKgfZprFBFJOXfo3x8uuQSqVoXXXw9LrYpsqrQGu7uPY+3n0Evv48BF6alIRCT9Fi6ELl3g2WfDqmxPPQU1a8ZdleQKDZ4QEUmjiROhXj144QW4/XYYPlyhLsmlYBcRSYPVq+Gf/4SmTUM3/NixcM01ULFi3JVJrtGyrSIiKfbNN9C2LbzzDrRoERZx2XbbuKuSXKVgFxFJoaFDoV07WLIEHn0UOnYEW++VRiKbRl3xIiIpsGwZXH55uNK9Zk2YNg3OP1+hLqmnFruISJJ98gm0bAnTp0O3bnDXXVClStxVSb5QsIuIJNGTT0LXrrD55vDyy3DaaXFXJPlGXfEiIkmwaBG0aRPOpzdsCDNnKtQlHgp2EZFNNGUKNGgAAwfCzTfDiBFQq1bcVUm+UrCLiGyk1avh7rvhkENg+XIYNQpuuEFj0yVeOscuIrIRvv8+dLsPGwannw6PPQbbbx93VSJqsYuIlNvw4VC3LoweDQ8/DIMHK9QlcyjYRUTKKiqCggKoUCF8LyoCQnf7VVfB8cfDDjuEc+sXXKCx6ZJZ1BUvIlJaURF07gwlJeH+7NnQuTOffb8l5z57GpMnh5XZ7rknLLcqkmkU7CIipfXs+XuoRwaWnMYFVxxFxW3g+efhrLNiqk0kAeqKFxEpbc6c324uphrt6U8rBvJXn0VxsUJdMp+CXUSktNq1AZhBPRoyjSdox3XcyujabalTJ+baRBKgYBcRKcVv68X9la7gICaymC15h6O5tWpvNrv9lrhLE0mIgl1EJDJ/PpwyqBWXrriL47cYw0zqcWSdL6BfP2jVKu7yRBKii+dERICRI0N2//gj9OkD3bodi9n8uMsSKTe12EUkr61YAdddB0cfDVtvDZMmwcUXa2y6ZC+12EUkb335JZx7LkyYAB06hJZ6tWpxVyWyaRTsIpKXnn8eOnUKC7kMHAgtW8ZdkUhyqCteRPJKSUmYWK5FC/jLX6C4WKEuuUXBLiJ54733oLAQHn00zPk+bhzssUfcVYkkl4JdRHKeOzz0EDRqBAsWhNXZeveGSpXirkwk+XSOXURy2k8/QceO8PLLcMIJ8MQTsOOOcVclkjpqsYtIzhozJqybPmQI/Otf4btCXXKdgl1Ecs7KlXDTTXDkkVClCowfD5dfHpZXF8l16ooXkZzy1VdhBrmxY6FNG3jwQdhqq7irEkkfBbuI5IyXXw4TzaxYAU8+GYJdJN9sVMeUmVUzs4rJLkZEZGMsXQoXXQSnnx6Gr02frlCX/JVQsJtZBTM718yGmNk84L/At2b2gZndZWZ7p7ZMEZG1+/BDaNIkDGe7/PJwPn1v/UaSPJZoi30ksCdwDbCzu+/m7jsChwETgd5m1jpFNYqI/Il7WE21sBC++w6GDg1XvleuHHdlIvFK9Bz7Me6+ouyD7v4TMBgYbGaa6kFE0uLnn8O0sM8/D8ccE86n77JL3FWJZIaEWuxrC/WN2UdEZFONHw/16sFLL4XZ44YNU6iLlLbBYDezY83sUTOrF93vnPKqRETKWLUKevWCww8P49HHjQvzvWtsusgfJfJfoitwJdDazI4C6qW0IhERgKIiKCiAChX4plZjjj3wO667Dv7+d5gxI1wwJyJ/lkiwz3f3n939CuA4oFGKaxKRfFdUFE6iz57N634SB349lEkfbsXjnSYwcCBss03cBYpkrkSCfciaG+5+NfBk6soREQF69mRZyUq6cx+n8Dq1mMs0GtJheEvM4i5OJLNtMNjd/ZU1t83s78CA6PZ1ZvaimTVIXXkiko8+nl2Fg5hIH7pzMX2YyEHsy8cwZ07cpYlkvPJednK9u/9iZk2B44EngIeTX5aI5CN3GDAAGto0vmI3XuUU+tCdKiwLO9SuHWt9ItmgvMG+Kvp+MvBw1JrXdBAisskWLQqLt7RvD432/YWZVQ7iFF7/fYeqVcNl8SKyXuUN9q/N7BGgBTDUzDbfiNcQEfmDyZOhfn147jm49VZ4+72d2fWxm6FOHTAL3/v1C8kvIutl7p74zmZVgROA99z9f2a2C/BXdx+eqgITUVhY6FOnTo2zBBHZCKtXw913Q8+eULMmDBwIhx4ad1Uimc/Mprl74dq2JTSlrJmZByXAi2sed/dvgW9L75OMgkUk9333HbRtC2+9BWeeCY8+CtttF3dVItkv4UVgzOxiM/vDlStmVtnMjjKzJ4B2yS9PRHLRsGFQty6MHQt9+4Y53xXqIsmR6CIwJwAdgGfMbHfgZ6AKUBEYDtzr7sWpKFBEcsfy5aHb/e67Yf/9YcSI8F1EkiehYHf3X4GHgIeiVdyqA0vd/efyvJmZ9QeaA/Pc/YC1bG8GvAJ8ET30orvfUp73EJHM9Omn0LIlTJ0KF14YlljdYou4qxLJPYm22H8TreL27Ua+3wDgAdY/e91Yd2++ka8vIhmoqAguuAA22wwGD4Yzzoi7IpHcldahau4+Bvgpne8pIvFZvBjatYPWrcNSqzNnKtRFUi0Tx6AfbGYzzewNM1vn2Tcz62xmU81s6vz589NZn4gkYPp0aNAAnn4abrgBRo7UxHEi6ZBpwT4dqOPudYF/Ay+va0d37+fuhe5eWKNGjXTVJyIb4A733gsHHQQlJeECuZtvDt3wIpJ6Gwx2MzvWzB41s3rR/c6pKsbdF7n74uj2UKCSmVVP1fuJSHLNmwfNm8Pll8OJJ4au9yOOiLsqkfySyN/QXYH2wHVmtj1QL1XFmNnOwPfu7mbWmPCHx4+pej8RSZ533gnn0hcsgAcegK5d0RKrIjFIJNjnR8ParjCz3kCjjX0zM3sGaAZUN7O5wI1AJQB37wucBVxoZiuBpcA5ms1OJLOtWAE33gi9e8Nf/gJvvhkmnxGReCQS7EPW3HD3q83s4o19M3dvuYHtDxCGw4lIFvjiCzj3XJg4Ec4/H+67D6pVi7sqkfy2wWB391fM7PLSj5W97+73JLswEclszz4LnaMrbgYNgrPPjrceEQkSvU51q+j7Xwhd8a9G908BxiS7KBHJXEuWQPfu8Pjj4cr3gQNh993jrkpE1khouJu73+zuNxOmkm3g7j3cvQfQEKiVygJFJEZFRVBQABUqQEEBM28fQmEh9O8P11wDY8Yo1EUyTXlHltYGlpe6vxwoSFo1IpI5iopCX3tJCQ48OLs5V/Q8mu23LeGtt6py9NFxFygia1PeYH8KmGxmL0X3/wY8kdSKRCQz9OwJJSX8yPZ0oD+vchonMYQBW15PjaOnx12diKxDuYLd3XuZ2RvAYYAD7d19RkoqE5F4zZnDaA6nFUXMY0fu5VK6cz/2tQani2SyjZlSdhWwutSXiOSYlSvhxq3v5ShGUJUSJnIQl3I/BprwXSTDlSvYzaw7UES4iG5H4OlNGdcuIplnzhxo1gxuWdidNhWfYRoNaUDUMVe1KvTqFWt9IrJ+5T3H3hFo4u5LAMzsDmACYcEWEclyL74IHTuGFvvTT0MrDHruAHOWhJZ6r17QqlXcZYrIepQ32I3QFb/GqugxEcliS5eGhVv69oXCQnjmGdhrL4BWCnKRLFPeYP8PMKnMVfGPJ7UiEUmrDz4Is8Z98AFceSXcdhtUrhx3VSKyscp7Vfw9ZjYaOJTQUtdV8SJZyh369YNLL4Wttw6Ltxx/fNxVicimKm+LHXefBkxLQS0ikiYLFkCnTjB4MBx7LDz5JOy8c9xViUgyJHRVvJmNi77/YmaLSn39YmaLUluiiCTTuHFhWdVXXoE77wwtdYW6SO5IqMXu7k3NzID93X1OimsSkRRYtQpuvx1uuilM//7uu9C4cdxViUiyJTyO3d0deGmDO4pIxvn6azjmGLjhhnCh3IwZCnWRXFXemecmmlmjlFQiIinx6qtw4IEwZQoMGBDWdtl667irEpFUKW+wH0kI98/MbJaZvWdms1JRmIhsml9/hUsugdNOgzp1YNo0aNcOTDNPiOS08l4Vf2JKqhCRpPrvf+Gcc2DmTOjeHe64AzbfPO6qRCQdyttin0NY2a2du88mrPC2U9KrEpGN4g6PPw4NG4bz6q+9Bvfdp1AXySflDfaHgIOBltH9X4AHk1qRiGyUhQuhZUs4/3xo0iS01ps3j7sqEUm38gZ7E3e/CPgVwN0XAJp8UiRmkyZB/frwwgthnZa33oKaNeOuSkTiUN5gX2FmFQld8JhZDbQmu0hqFRWFgecVKoTvRUW/bVq9Gnr3hqZNw+0xY+Daa6FixdiqFZGYlffiuT6Esew7mlkv4CzguqRXJSJBURF07gwlJeH+7NnhPvDtUa1o2xbefhv+/vcw7/u228ZXqohkhvIuAlNkZtOAowmLwPzN3T9KSWUiAj17/h7qa5SU8Mblb9HOW7F4cQj088/XMDYRCRIKdjOrAlwA7AW8Bzzi7itTWZiIAHP+OIPzMipzLbdzz7we/PWvMGgQ7LdfTLWJSEZK9Bz7E0AhIdRPBO5OWUUi8rvatX+7+T/24hDGcw896Lrlk0yapFAXkT9LNNj3c/fW7v4I4bz64SmsSUTW6NULqlblKVrTgOl8we68VPlsHuxbkS22iLs4EclEiQb7ijU31AUvkj6/nNqKNvXfpy1PUZ8ZzNz1ZP7W/1Ro1Sru0kQkQyV68VzdUuuuG7BFdN8IC79pSQmRJJs2LUwL+/nnu3PTTXDddYdTseKEuMsSkQyX6HrsGhUrkiarV8O998I118BOO8GoUXDYYXFXJSLZorzj2EUkhb7/Hs47D958E/72tzDv+/bbx12ViGST8s48JyIp8tZbULcujBwJDz4IL76oUBeR8lOwi8RsxQq46io47rgQ5FOmQNeumnBGRDaOuuJFYvT552FFtsmTw0yx994LVavGXZWIZDMFu0hMBg2CLl1Cy/y558J87yIim0pd8SJptmQJdOgQWuoHHBDWTVeoi0iyKNhF0qi4GBo2hAEDwvouo0dDnTpxVyUiuUTBLpIG7tCnDzRpAosWhaVWb7sNNtPJMBFJMv1aEUmxH36A9u3h9deheXP4z3+gevW4qxKRXKUWu0gKjRwZxqYPHw733w+vvqpQF5HUUrCLpMDKlXD99XD00bDlljBxIlxyicami0jqqSteJMlmz4Zzz4Xx40MXfJ8+IdxFRNJBwS6SRC+8AJ06wapVMHBgGNImIpJO6ooXSYKSkjDZzN//DnvvDTNmKNRFJB4KdpFN9N570KgR9OsH//gHjBsHe+4Zd1Uikq8U7CKJKCqCggKoUCF8LyrCHR5+GBo3hh9/hGHD4I47oHLluIsVkXyW1mA3s/5mNs/M3l/HdjOzPmb2qZnNMrMG6axPZK2KisIKLbNnh5lmZs/mp05XcWbjOXTtCkccEaaFPe64uAsVEUl/i30AcMJ6tp8I7B19dQYeTkNNIuvXs2c4iR4ZS1PqLR3P61N35u67YehQ2GmnGOsTESklrcHu7mOAn9azy2nAkx5MBLY1s13SU53IOsyZA8AqKnAL19OMUVRmOeM5lB49Qu+8iEimyLRfSbsCX5W6Pzd6TCQ+tWvzFbU4ihHcyC205Bmm04DCOvPjrkxE5E8yLdjXNi+Xr3VHs85mNtXMps6fr1+wkjqvnPEE9ShmGg15grY8TRu2rroKevWKuzQRkT/JtGCfC+xW6n4t4Ju17eju/dy90N0La9SokZbiJL/8+it06wZ/u/cICgpges1TaGtPh3VW+/WDVq3iLlFE5E8ybea5V4FuZjYIaAIsdPdvY65J8tBHH8HZZ4cx6pdfDrffvgObbz4y7rJERDYorcFuZs8AzYDqZjYXuBGoBODufYGhwEnAp0AJ0D6d9Ym4w+OPhwVbttwShgyBk06KuyoRkcSlNdjdfb2TbLq7AxelqRyRP/j55zAt7HPPhVXZnnoKdtGYDBHJMpl2jl0kFhMmQP36MHgw/POfYf10hbqIZCMFu+S1Vavg9tvhsMPC/XHj4OqrNTZdRLJXpl08J5I233wDbdrAiBHQogU88ghsu23cVYmIbBoFu+SlIUPgvPPCTLGPPQYdOoCtbRYFEZEsow5HySvLlsFll0Hz5lCzJkydCh07KtRFJHeoxS5545NP4JxzYMYMuPhiuPNOqFIl7qpERJJLwS45zx2efBIuugg23xxeeQVOPTXuqkREUkNd8ZLTFi2C1q3D+fTCwrBuukJdRHKZgl1y1pQpYWz6oEFwyy3wzjtQq1bcVYmIpJaCXXLO6tVw111wyCGwYgWMHg3XXw8VK8ZdmYhI6ukcu+SU77+Hdu1g2DA444wwlG277eKuSkQkfdRil5wxfDgceGBoofftCy+8oFAXkfyjYJest3w5/OMfcPzxUKNGOLfepYvGpotIflKwS/YoKoKCgjCRe0EBFBXx2WfQtGk4p96lC0yeDAccEHehIiLx0Tl2yQ5FRdC5c5gDFmD2bIo6vMOFFVpQsUolXngBzjwz3hJFRDKBWuySHXr2/C3UF1ON8/gPrZf350AvprhYoS4isoaCXbLDnDkAzKAeDZnGk7Tlem5h1LJDqFMn5tpERDKIgl2ygu9Wm/vozkFMZAnVGMFR3MKNbFZn17hLExHJKDrHLhlv/nxov/27DJmzK6fwKv3pQHV+hKpVoVevuMsTEckoarFLRhsxAurWhbc+3JU+bafwSu1LqG4/QZ060K8ftGoVd4kiIhlFLXbJSCtWwI03Qu/esM8+8MYbULduI+DLuEsTEcloCnbJOF9+CS1bwsSJ0LEj3H8/VKsWd1UiItlBwS4Z5fnnoVOnsIb6oEFw9tlxVyQikl10jl0yQklJCPQWLWDffaG4WKEuIrIxFOwSu1mzoLAQHn8crr4axo6F3XePuyoRkeykYJfYuMODD0LjxrBgQVid7Z//hEqV4q5MRCR76Ry7xOLHH8OFca+8AieeCAMGwI47xl2ViEj2U4td0m7MGKhXD4YOhXvugddfV6iLiCSLgl3SZuVKuOkmOPJIqFIFJkyAyy4Lq7CKiEhyqCte0uKrr8IkcWPHQtu28MADsNVWcVclIpJ7FOySci+9FM6nr1gBTz0FrVvHXZGISO5SJ6ikzNKl0LUrnHEG7LEHzJihUBcRSTUFu6TEBx+EYWwPPww9esD48bDXXnFXJSKS+9QVL0nlDo8+CpdeGs6hv/EGnHBC3FWJiOQPtdglaRYsCFPCdukCTZvCzJkKdRGRdFOwS1KMHx/Gpr/8MtxxB7z5Juy8c9xViYjkHwW7bJJVq6BXLzj8cKhYEcaNg3/8Q2PTRUTiol+/Uj5FRVBQABUq8HWtJhx74Hdcd13ogp8xA5o0ibtAEZH8pmCXxBUVQefOMHs2r/nJ1P16CJM+3Ir+nSZQVATbbBN3gSIiomCXxPXsya8lq+jOfZzKa9RiLtNoSPvhLTGLuzgREQENd5Ny+Hh2Fc5hAsXUpzv30ZurqcIymKNUFxHJFAp22SD3sKxqN5vOFl7CazSnOUN+36F27dhqExGRP1JXvKzXwoVh8ZYOHaDxvouYWeWgP4Z61arhsngREckICnZZp0mToH59eO45uO02ePu9ndn1sZuhTh0wC9/79QvJLyIiGUFd8fInq1fDXXfBdddBzZowZgwccki0sVUrBbmISAZTsMsffPcdtGkDb78NZ54Z5n3fbru4qxIRkUSpK15+8+abcOCB8O67oYf9+ecV6iIi2UbBLixfDldcASeeCDvtBFOnQqdOaGy6iEgWSnuwm9kJZvaxmX1qZlevZXszM1toZsXR1w3prjGffPppOH/+r39B164weTLst1/cVYmIyMZK6zl2M6sIPAgcC8wFppjZq+7+YZldx7p783TWlo+efhouvBAqVYIXX4TTT4+7IhER2VTpbrE3Bj5198/dfTkwCDgtzTXkvV9+gbZtw0Vy9epBcbFCXUQkV6Q72HcFvip1f270WFkHm9lMM3vDzPZf2wuZWWczm2pmU+fPn5+KWnPStGnQoEFYz+XGG2HkSE0cJyKSS9Id7Gu7HMvL3J8O1HH3usC/gZfX9kLu3s/dC929sEaNGsmtMgetXg333AMHHwy//gojRsBNN8FmGvAoIpJT0h3sc4HdSt2vBXxTegd3X+Tui6PbQ4FKZlY9fSXmnnnzoHlz6NEDTjopdL0fcUTcVYmISCqkO9inAHub2e5mVhk4B3i19A5mtrNZGGhlZo2jGn9Mc5054+23oW7d0EJ/8EF46SXYYYe4qxIRkVRJa0esu680s27AMKAi0N/dPzCzC6LtfYGzgAvNbCWwFDjH3ct218sGrFgBN9wAd9wB++4Lw4aFyWdERCS3WS5kZmFhoU+dOjXuMjLGF19Ay5ZhEZdOneDee6FatbirEhGRZDGzae5euLZtunQqxwwaBF26hFnjnn0WWrSIuyIREUknTSmbI5YsgY4dQ0t9v/3CBXIKdRGR/KNgzwHFxdCwIfznP3DttWGZ1YKCuKsSEZE4KNizmDv8+9/QpAksWhSugO/VK0wRKyIi+Unn2LPUDz9Ahw7w2mtw8smhta55ekRERC32LDRqVBibPmwY3HdfCHeFuoiIgII9q6xcCddfD0cdFYavTZgA3btr3XQREfmduuKzxOzZ0KoVvPsunHdeOLe+5ZZxVyUiIplGwZ4FBg+G88+HVavCqmznnht3RSIikqnUFZ/Bli6FCy6As86CvfeGGTMU6iIisn4K9gz1/vvQqBE88ghceSWMGwd77hl3VSIikunUFZ9h3EOYX3YZbL11uPL9uOPirkpERLKFWuwZ5KefQrf7hRfC4YfDrFkKdRERKR8Fe4YYNw7q1YNXX4W77oI33oCddoq7KhERyTYK9pitWgW33AJHHAGVK8P48XDFFVBBn4yIiGwEnWOP0dy50Lo1jB4dxqg/9FA4ry4iIrKxFOwxefVVaN8eli2DAQOgbVvNICciIptOHb5p9uuvcPHFcNppUKcOTJ8O7dop1EVEJDkU7Gn00UdhidUHHoBLLw1zve+zT9xViYhILlFXfBq4Q//+cMklULUqvP56WGpVREQk2dRiT7GFC6FlyzDX+0EHwcyZCnUREUkdBXsKTZwYxqa/8ALcfjsMHw41a8ZdlYiI5DIFewqsXg29e0PTpqEbfuxYuOYaqFgx7spERCTX6Rx7kn37LbRpA++8Ay1ahHnft9027qpERCRfKNiTaOjQMHRtyRJ49FHo2FHD2EREJL3UFZ8Ey5bB5ZeHi+Jq1oRp08LFcgp1ERFJN7XYN9Enn4Sr3qdPh27dwgIuVarEXZWIiOQrBfsmePJJ6NoVNt8cXn45zCYnIiISJ3XFb4RffgkXyLVrBw0bhrHpCnUREckECvZymjoV6teHgQPh5pthxAioVSvuqkRERAIFe4JWr4a774aDD4bly2HUKLjhBo1NFxGRzKJz7An4/vvQ7T5sGJx+Ojz2GGy/fdxViYiI/Jla7BswfDjUrQujR8PDD8PgwQp1ERHJXAr2dVi+HK66Co4/HnbYAaZMgQsu0Nh0ERHJbOqKX4vPPw9j0ydPhi5d4J57wnKrIiIimU7BXsYzz4Qwr1gRnn8ezjor7opEREQSp674yOLF0L49nHsu/PWvUFysUBcRkeyjYAdmzAgTzTzxBFx3XbhQrk6duKsSEREpv7wOdne4/3446KDQYn/nHbj1VthMJyhERCRL5W2EzZ8fut6HDIFTToH+/aF69birEhER2TR52WIfOTKMTX/rLejTB155RaEuIiK5Ia+CfeXKcA796KNh661h0iS4+GKNTRcRkdyRN13xX34ZrnifMAE6dAgt9WrV4q5KREQkufIi2J9/Hjp1Cgu5DBwYJp8RERHJRTndFV9SAp07Q4sW8Je/hLHpCnUREcllORvs770HhYXw6KNhzvdx42CPPeKuSkREJLVyLtjd4aGHoFEjWLAgrM7WuzdUqhR3ZSIiIqmX9mA3sxPM7GMz+9TMrl7LdjOzPtH2WWbWINHX/uknOPNMuOgiOPJImDkTjj02ufWLiIhksrQGu5lVBB4ETgT2A1qa2X5ldjsR2Dv66gw8nMhrjxkTxqa//jr8619h4pkdd0xi8SIiIlkg3VfFNwY+dffPAcxsEHAa8GGpfU4DnnR3Byaa2bZmtou7f7uuF/3mm9BC32MPGD8+nFsXERHJR+nuit8V+KrU/bnRY+Xd5w++/TaMUZ8+XaEuIiL5Ld0t9rXN8eYbsQ9m1pnQVQ+w7Omn7f2nn97E6rJPdeCHuIuIgY47v+i484uOOzHrXIM03cE+F9it1P1awDcbsQ/u3g/oB2BmU90979rqOu78ouPOLzru/JLM4053V/wUYG8z293MKgPnAK+W2edVoG10dfxBwML1nV8XERGR36W1xe7uK82sGzAMqAj0d/cPzOyCaHtfYChwEvApUAK0T2eNIiIi2Sztc8W7+1BCeJd+rG+p2w5cVM6X7ZeE0rKRjju/6Ljzi447vyTtuC3kqIiIiOSCnJtSVkREJJ9lVbCncjraTJbAce9rZhPMbJmZXRFHjamQwHG3ij7nWWY23szqxlFnsiVw3KdFx1xsZlPNrGkcdSbbho671H6NzGyVmZ2VzvpSJYHPu5mZLYw+72IzuyGOOpMtkc87OvZiM/vAzEanu8ZUSODzvrLUZ/1+9G99+3K9ibtnxRfhYrvPgD2AysBMYL8y+5wEvEEYC38QMCnuutN03DsCjYBewBVx15zG4z4E2C66fWIefd5b8vtptAOB/8ZddzqOu9R+IwjX6ZwVd91p+rybAa/HXWsMx70tYVbS2tH9HeOuOx3HXWb/U4AR5X2fbGqx/zYdrbsvB9ZMR1vab9PRuvtEYFsz2yXdhSbZBo/b3ee5+xRgRRwFpkgixz3e3RdEdycS5jzIdokc92KP/tcD1VjLBE5ZKJH/3wAXA4OBeeksLoUSPe5ck8hxnwu86O5zIPyeS3ONqVDez7sl8Ex53ySbgj0l09FmgVw8pkSU97g7Enprsl1Cx21mp5vZf4EhQIc01ZZKGzxuM9sVOB3oS+5I9N/5wWY208zeMLP901NaSiVy3PsA25nZKDObZmZt01Zd6iT8e83MqgInEP6QLZe0D3fbBEmbjjbL5OIxJSLh4zazIwnBngvnmhM6bnd/CXjJzA4HbgWOSXVhKZbIcd8HXOXuq8zWtntWSuS4pwN13H2xmZ0EvExY/TKbJXLcmwENgaOBLYAJZjbR3T9JdXEpVJ7f56cA77r7T+V9k2wK9qRNR5tlcvGYEpHQcZvZgcBjwInu/mOaakulcn3e7j7GzPY0s+runs3zaydy3IXAoCjUqwMnmdlKd385LRWmxgaP290Xlbo91MweypPPey7wg7svAZaY2RigLpDNwV6e/9/nsBHd8EBWXTy3GfA5sDu/X3Swf5l9TuaPF89NjrvudBx3qX1vIncunkvk865NmKHwkLjrTfNx78XvF881AL5ecz9bv8rz7zzafwC5cfFcIp/3zqU+78bAnHz4vIH/A96J9q0KvA8cEHftqT7uaL9tgJ+AahvzPlnTYvc8nY42keM2s52BqcDWwGozu5RwpeWidb1upkvw874B2AF4KGrFrfQsXzwiweM+k7CewgpgKXC2R78NslWCx51zEjzus4ALzWwl4fM+Jx8+b3f/yMzeBGYBq4HH3P39+KredOX4d346MNxDb0W5aeY5ERGRHJJNV8WLiIjIBijYRUREcoiCXUREJIco2EVERHKIgl1ERCSHKNhF0ixarWnNyk3PR1NHbuprnmdmD2xgnwIzS/pwoRS+bjMzO6TU/QGJrOhmZluY2Wgzq1iO9+pmZlk/PFYEFOwicVjq7vXc/QBgOXBB3AVlqGaEFfzKqwNh8ZBV5XhOf+CSjXgvkYyjYBeJ11hgLzOrZmb9zWyKmc0ws9Pgt5b4i2b2ppn9z8zuXPNEM2tvZp9E61QfWurxP7RszWxx2Tct28I3s9fNrNma/c3sjmjhjbfNrHG0EMfnZnbq+g7GzCqa2V3Rccwysy7R482i13jBzP5rZkUWzSpkZidFj40zsz5RLQWEP3gui3o3Dove4nAzGx/Vsq7WeyvglVLvO9rMnot+Vr3NrJWZTTaz98xsTwB3LwG+NLPG6zs+kWygYBeJiZltRlhH/j2gJ2Hd5UbAkcBdZlYt2rUecDbwV+BsM9vNwnLENxMC/VhgvySWVg0Y5e4NgV+A26L3OB24ZQPP7QgsjI6jEdDJzHaPttUHLo1q3QM41MyqAI8Q5vpvCtQAcPcvCau43Rv1boyNXmMXwmI/zYHeZd/czCoDe0TPX6Mu0J3w82sD7OPujQlrDFxcar+pwGGIZLmsmVJWJIdsYWbF0e2xwOPAeOBUM7sierwKYS58gHfcfSGAmX0I1CEsgjLK3edHjz9LWOYyGZYDb0a33wOWufsKM3sPKNjAc48DDizVmt6GsBLZcsLaDXOjeouj11oMfO7uX0T7PwN0Xs/rv+zuq4EPzWyntWyvDvxc5rEp7v5t9L6fAcNLHduRpfabB+y7geMTyXgKdpH0W+ru9Uo/EHVLn+nuH5d5vAmwrNRDq/j9/+265oNeSdQbF71u5fXtE6lS6vaKUnORr17z/u6+OuplWB8DLnb3YWWOo9k6jqO866+Wfo21PXcpfzyWss9ZXer+av74O7BK9HyRrKaueJHMMAy4uNR55/ob2H8S0MzMdjCzSsDfS237krCONcBpQKW1PP9LoJ6ZVTCz3QirhiXDMMKCJZUAzGyfUqcU1ua/wB7ROXUIpxzW+AXYqjxv7u4LgIpRF3957UNYQUwkqynYRTLDrYQAnhUNHbt1fTtHXcs3AROAt4HppTY/ChxhZpOBJsDaVoh6F/iC0B19d5nnb4rHgA+B6dFxPMJ6egbdfSnQFXjTzMYB3wMLo82vAaeXuXguEcMJ5+HL61DCz1Ikq2l1NxGJlZlt6e6Lo96KB4H/ufu9m/B69YHL3b1NKp8jkqnUYheRuHWKLqb7gHCx3SOb8mLuPgMYWZ4JaggX3V2/Ke8rkinUYhcREckharGLiIjkEAW7iIhIDlGwi4iI5BAFu4iISA5RsIuIiOQQBbuIiEgO+X9Bch3SaOCsQAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# \n", "#\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "#\n", "gradient = 4.0\n", "intercept = 0.0\n", "#\n", "plt.figure(figsize = (8, 6))\n", "#\n", "# Plot the data\n", "plt.plot(lengths, Tmean**2, marker = 'o', color = 'r', linestyle = '')\n", "#\n", "# Draw a straight line\n", "xPoints = np.zeros(2)\n", "xPoints[1] = 1.1*np.max(lengths)\n", "plt.plot(xPoints, gradient*xPoints + intercept, color = 'b', linestyle = '-')\n", "#\n", "plt.xlim(0, 1.2*np.max(lengths))\n", "plt.ylim(0, 1.2*np.max(Tmean**2))\n", "plt.title('Data with hand fit', fontsize = 14)\n", "plt.xlabel('Pendulum length (m)')\n", "plt.ylabel('Period$^2$ (s$^2$)')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Exercise 4\n", "\n", "Now use the value of the gradient from the graph above to determine $g$." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " \n", "Acceleration due to gravity 9.870 m/s^2.\n" ] } ], "source": [ "# \n", "#\n", "gravAccel = (4*np.pi**2)/gradient\n", "print(\" \")\n", "print(\"Acceleration due to gravity {:5.3f} m/s^2.\".format(gravAccel))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Errors\n", "\n", "You should observe that the line of best fit does not go exactly through all of your points. This is because in any experiment we cannot measure things perfectly; our measurements are subject to uncertainties or errors. We must be able to estimate these uncertainties and understand how they affect our overall result. In the lab book, the theory underpinning these uncertainties and how you should handle them is laid out, but we will some of the relevant aspects here.\n", "\n", "In this experiment you took multiple readings of the period. Taking multiple readings is one way that we can estimate and ultimately reduce what are known as *statistical* errors. These are random uncertainties that influence our experimental results, such as how accurately we can start and stop stopwatch. These are random in the sense that we are just as likely to push the stop button on the watch a little too early as we are to press it too late. So we are just as likely to get a measurement that it smaller than the true value as one that is larger. Further, we are less likely to get a measurement that is too long (or too short) by a large amount, than to get a measurement that is close to the true value. (In the second half of the experiment, when we are dealing with radioactive decay, we will meet another type of statistical uncertainty associated with counting only a finite number of decays.)\n", "\n", "In addition to statistical errors, there are also *systematic* errors. These are not random, but produce systematic shifts in your measurements. Here is an example:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Exercise 5\n", "\n", "If your stopwatch is 10% slow, all your measured times would be 10% too short. What would then happen to your average times, to your graph, to the gradient you measure from it and to your final measurement of $g$?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Quantifying errors\n", "\n", "Because you took multiple readings, you can determine the error on the measured period for those readings. This error can be used to estimate the uncertainty in the $T^2$ values in your graph, the errors in the gradient of the straight line and hence the precision of the results for $g$.\n", "\n", "The first step is to calculate the standard deviation $\\sigma_{N-1}$ (also known as the standard error on a single reading). This is a measure of the error on any one of your measurements of the period and is calculated using the equation:\n", "\n", "$$\n", "\\sigma = \\sqrt {\\frac{{\\sum\\limits_{i = 1}^N {{{({T_i} - \\overline T )}^2}} }}{{N - 1}}}.\n", "$$\n", "\n", "The standard deviation tells us the instrinsic spread of our individual measurements, due to the precision of the methods we are using, and does not decrease as we collect more data.\n", "However, when we take more measurements of the same quantity and average them, the precision of that mean value does increase. The uncertainty in value of the mean is given by:\n", "\n", "$$\n", "\\Delta T = \\frac{{{\\sigma _{N - 1}}}}{{\\sqrt N }}.\n", "$$\n", "\n", "This is the error on the mean, and as you can see, the more measurements you take the smaller the uncertainty will be.\n", "\n", "While you can do all these calculations by hand, it is often more efficient to use the statistics mode on your calculator or to use an Excel spreadsheet or Jupyter Notebook to automate the calculation. This is a good habit to get into for future experiments!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Exercise 6\n", "\n", "Calculate the standard deviation of the period measurements, $T$ (taken over 20 oscillations), then the values of $\\Delta T$ for each pendulum length, using the cell below. " ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " \n", "Standard deviations in seconds\n", " [0.0213 0.0213 0.0214 0.0234 0.0178 0.0255 0.0239 0.0218 0.0245 0.0264]\n", "Errort in T in seconds\n", " [0.0048 0.0048 0.0048 0.0052 0.004 0.0057 0.0053 0.0049 0.0055 0.0059]\n" ] } ], "source": [ "# \n", "#\n", "sDev = np.zeros(nLen)\n", "for n in range(0, nLen):\n", " sDev[n] = np.sqrt(np.sum((Tdata[n, 0:nMeasPerLen] - \\\n", " Tmean[n])**2)/(nMeasPerLen - 1))\n", "Terror = sDev/np.sqrt(nMeasPerLen)\n", "#\n", "np.set_printoptions(precision = 4)\n", "print(\" \")\n", "print(\"Standard deviations in seconds\\n\",sDev)\n", "print(\"Errort in T in seconds\\n\",Terror)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "You now need the error on $T^2$. Following the rules laid out in the handbook, which you tested in your homework, you can calculate $\\Delta T^2$ via the equation:\n", "\n", "$$\n", "\\frac{{\\Delta {T^2}}}{{{T^2}}} = 2\\frac{{\\Delta T}}{T}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Exercise 7\n", "\n", "Calculate the values of $\\Delta T^2$ in the cell below and use these to add error bars to your graph. Add also your estimated errors in the measurements of the length of the pendulum." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " \n", "Error in T^2 in s^2\n", " [0.0074 0.0085 0.0097 0.0116 0.0094 0.0144 0.0144 0.0139 0.0162 0.0183]\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfYAAAGECAYAAADEAQJ2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA1PUlEQVR4nO3dd5gUVdrG4d8rgijGFYwoGNdVV0ABEyrmsJhdFBEQEBAUs65rYjDnwIoiCmJAUcQsCCpZchgw7fqZxYiKZIFh3u+PU2gzDtAN3V0dnvu65prururut6Zhnjmn6pxj7o6IiIgUhvXiLkBERETSR8EuIiJSQBTsIiIiBUTBLiIiUkAU7CIiIgVEwS4iIlJAFOwi8jszczM7Yw379DOz11N83abRa9dctwrXjpm9b2Yla9ino5l9ZWblZlZiZuea2YIslSiSNgp2kTSIws6jr2Vm9qOZjTCzC8ysaoqvFWcIbgu8FtVRN6qjYQx1ZJWZbQH0BO4CtgfuBp4Ddk7Yp8TM3o+nQpHkKdhF0udtQjDWBY4hBGR3YIyZ1YixrqS5+/fuviTuOmJQB1gfeN3dv3P3Be6+2N1/jLswkVQp2EXSZ0kUjN+4e6m73ws0BfYFrlqxk5mdY2aTzWx+1LIfaGbbR9vqAiOiXWdHLeZ+0bbjzGyMmc0xs1/MbKiZ/W1VxZjZ36LnbxPd38jMlprZkIR9OpjZ/yXcT+yK/zz6Pjl6fGSF17/YzL6J6nnczDZK4mdUz8wmmtkiM5tiZvsmvN6WZvasmc0ys8Vm9oGZta3wniPN7CEzu9XMfop+fneb2XoJ+2xlZq9Er/GlmbVbXUFmdi4wPbr7WXSsdRO74qN9ugF7JfTMnJvE8YpknYJdJIPc/X3gTeD0hIerEUKiHtAMqAk8G237OmHfvQg9ABdH92sA9wONCX8wzAVeM7Nqq3jvj4Afon0BDo6e08TM1o8eawqMXEX5jaPvx0V1nJaw7RBgb+Ao4Ezg1IQ6V+c24GrCHzs/A/3NzKJt1YFphJ/JXsADwCNmdmSF12gJlAEHARcCl0Q1rNAP2DWq7RSgNaEXZVWei44RwjFvS/gcKu5zD/C/aPu20WMiOWf9Ne8iIuvoQ0LIAODufRO2fWZmnYGPzKy2u88ys1+ibT+6+08JzxuU+KJRa3YeIYzGruK9RwGHAwMIIf4CcDzQCBgPHAb8axXPnR19/9ndv6+wbR7Q2d3LotoHAkcSgnt1rnf3EVH9N0Z1bw/McvdvCOe4V+htZkcALYB3Eh7/0N1viG5/bGYdovd+1sx2j46vibu/G71PG+CzVRXk7ovN7OcVx7ziWP/4e+P3fRYAZZX8LERyilrsIplnwO+rLZnZvlFX8ZdmNh+YEm3acbUvYraLmT1jZp+a2TxCa3y9NTxvJH+02JsSuvlHAU3NbDdCqI5M9YAI4VqWcP9bYKsknjezwnNY8Twzq2Jm15rZTDP7OQrS0/jz8c2scD/xvf8GlAOTVmx09y8T3kuk4KnFLpJ5exK1GKOL6IYSLrRrBfxI6IofQ+iiX53XgG+ATtH3MkJvwOqeNxJ4KArxhtH9GoRW8E/AJ1FLOVXLKtx3kmsoJD5vxR87K553BXA5oUv/PWABcCt//oNhde9tiBQ5tdhFMsjM9iacv30hemgPQpBf4+6j3f2//Dm4lkbfqyS8zpaE1uit7v52dP58E9bwx3nCefZrCSH+I6HVfjBwNKtvrf+pjgxrArzm7k+5eynwKbB7iq/xEeH3WqMVD5jZjsB2aahvKdn7WYisNQW7SPpsYGbbmNl2ZlbPzC4jBOdUwrhogK+AJcCFZrazmf0DuKnC63xJaIX+w8xqmdnGwBxCC7uDme1qZocBvQit9jUZBZxDdLW9u39BOH9+GqsP9h+BxcCxZra1mW2WxHuti4+BI82siZntATwI7JTKC7j7/wgXKz5iZgeaWX3CxXSL01DfF0Cd6FRKTTPbIA2vKZJ2CnaR9DkK+I4Q3u8AJxHGsR/q7gsB3H020IZwtfaHhKvjL0t8kahrvBtwC6G1/aC7lxOu/N4HeJ8wmcr1hD8S1mQEoaU5MuGxkZU8tpLoHPpFwHmEc9SvJPFe6+JmwrnxIcBoYCHQfy1e51zCUL3hhNMXzxBCeV0NAgYTPtvZhNMZIjnH3H3Ne4mIiEheUItdRESkgCjYRURECoiCXUREpIAo2EVERAqIgl1ERKSAFMTMczVr1vS6devGXYaIiEhWTJ069Sd3r1XZtoII9rp16zJlypQ17ygiIlIAzOzLVW1TV7yIiEgBUbCLiIgUEAW7iIhIAVGwi4iIFBAFu4iISAFRsIuIiBQQBbuIiEgBUbCLiIgUEAW7iIhIAVGwi4iIFBAFu4iISAFRsIuIiBQQBbuIiEgBUbCLiIgUEAW7iIhIAVGwi4iIFBAFu4iISAFRsIuIiBQQBbuIiEgBUbCLiIgUEAW7iIhIAVGwi4iIFBAFu4iISAFRsIuIiBQQBbuIiMgq/PorvPVW3FWkRsEuIiJSifHjoX59OP10mDMn7mqSp2AXERFJsHw53HorHHIImMGwYbDFFnFXlbz14y5AREQkV3z7LbRqBcOHw5lnwiOPwGabxV1VahTsIiIiwBtvwLnnwqJF0KcPtG0bWuz5Rl3xIiJS1JYsgUsvhWbNYLvtYMoUaNcuP0Md1GIXEZEi9vHHcNZZMH06dO0Kd94J1avHXdW6UbCLiEjRcYcnnoALL4QNNoBXXoGTToq7qvRQV7yIiBSVefPgnHPCOfSGDWHGjMIJdVCwi4hIEZk8GRo0gAED4MYb4Z13oHbtuKtKLwW7iIgUvPJyuOsuOOggKCuD0aPh+uuhSpW4K0s/nWMXEZGC9sMP0KYNDB0Kp50Gjz2WXxPOpEotdhERKVhDh8I++8CoUdCrF7zwQmGHOijYRUSkAC1dCldeCccdB7VqhXPrnTrl79j0VGQ12M1sBzMbYWYfmdkHZnZxJfs0NbO5ZlYafd2QzRpFRCS/ffopNGkCd98dwnzSJNh777iryp5sn2MvAy5392lmtgkw1czecvcPK+w3xt2bZbk2ERHJc/37Q+fO4aK4F14IK7MVm6y22N39O3efFt2eD3wEbJ/NGkREpPAsWBDmeT/nnHBOvbS0OEMdYjzHbmZ1gQbAxEo2H2hmM8xsiJnttYrndzSzKWY2Zfbs2ZksVUREctj06bDffvDkk2EI28iRUKdO3FXFJ5ZgN7ONgUHAJe4+r8LmaUAdd68H/Ad4ubLXcPfe7t7Q3RvWqlUro/WKiEjucYf774cDDoCFC8NSqzfeCOsX+UDurAe7mVUlhHp/d3+x4nZ3n+fuC6Lbg4GqZlYzy2WKiEgOmz07rMZ26aVw7LGh671p07iryg3ZviregD7AR+5+7yr22SbaDzNrTKjx5+xVKSIiuWz4cKhXD95+G3r0CAu41FTz73fZ7rA4GGgFvGdmpdFj1wA7Arh7L+AMoLOZlQGLgbPc3bNcp4iI5Jhly6BbN7j9dth9dxgyJAS8rCyrwe7uY4HVTg/g7g8CD2anIhERyQdffAEtWsCECdC+PTzwANSoEXdVuanILzEQEZFcN3AgdOgQLpYbMADOPDPuinKbppQVEZGctHBhCPTmzWGPPcIFcgr1NVOwi4hI7igpAWDmTGjYEPr0gauvhjFjYKed4i0tX6grXkREckPTpjBqFD1rlXD55WEVtmHD4Kij4i4sv6jFLiIiOeHnaV8CcOGFcMQRMGOGQn1tKNhFRCReJSVgxpbzvwDAMQYPMbZ6qCTOqvKWgl1ERGJTVgYllFBlPWe3Db8OD7qHr+h8u6RGwS4iIrH4+uvQ5d69e1iVbdoPteMuqSDo4jkREcm6l14KE80sWwZPPRWCHQhTy8k6UYtdRESyZvFi6NIFTjsNdt45LLn6e6iDut/TQMEuIiJZ8cEH0LgxPPwwXH45jBsHu+4ad1WFR13xIiKSUe7w6KNwySWwySZh8Zbjjou7qsKlFruIiGTMnDlhSthOnaBJkzA2XaGeWQp2ERHJiHHjoH59ePlluOMOePNN2GabuKsqfAp2ERFJq+XL4ZZb4NBDoUoVGDsWrroK1lPiZIXOsYuISNp88024yn3kyLB++sMPw2abxV1VcVGwi4hIWrz2GrRtG4a09e0L554LZnFXVXzUMSIiIuvkt9/g4ovhpJNghx1g2rQQ8Ar1eCjYRURkrf3vf3DggdCjRwj3CRPgr3+Nu6ripq54ERFJmTv06xeWWN1ww9AN36xZ3FUJqMUuIiIpmjsXWraEdu3CTHIzZijUc4mCXUREkjZxIjRoAM8/DzffDG+/DdtvH3dVkkjBLiIia1ReHiaZadIkjFMfPRquvTaMU5fconPsIiKyWt9/D61ahdb5GWdA796wxRZxVyWromAXEZFVevNNaN0aFiwIgX7eeRrGluvUFS8iIn+ydClccQUcfzxsvTVMmQIdOijU84GCXUREVvJz1xIOOgjuuQe6dIFJk2DPPeOuSpKlYBcRkd899bdb2PLB7nz2Gbz4IvTsGcapS/5QsIuICPPnh3Pprf97LQClpXDqqfHWJGtHwS4iUuSmToXH65Tw5FOGE06i71jHwgn1kpJ4i5OUKdhFRIpUeTnce2+Y6/2uGiWMGulw2GFho3v4UrDnHQ13ExEpQj/+GJZVHTIETj4Z+vSBLbckLKSuS9/zmlrsIiJF5u23oV49GD48XBz30ktRqK/QrVtstcm6U7CLiBSJZcvg6qvhmGPCzHGTJoXhbH9qoKv7Pa+pK15EpAh8/jm0aBEWcenQAe67D2rUiLsqyQQFu4hIgRswADp1Ci3z556D5s3jrkgySV3xIiIFauFCaN8+tNT33DOMTVeoFz4Fu4hIASothf32g8cfh2uuCcus1q0bd1WSDQp2EZEC4g7/+Q/svz/MmxeugL/lFqhaNe7KJFt0jl1EpED89BO0awevvQb/+EdordeqFXdVkm1qsYuIFICRI8PY9KFD4f77Q7gr1IuTgl1EJI+VlcH118MRR4Tha+PHw8UXa/K4YqaueBGRPPXll9CyJbz7bpge9j//gY03jrsqiZuCXUQkDw0aBOedB8uXQ//+cPbZcVckuUJd8SIieWTxYjj/fDjjDNhtN5g+XaEuK1Owi4jkifffh0aN4JFH4MorYexY2GWXuKuSXKNgFxHJRQkLsbhDr14h1GfPDle+33knVKsWX3mSu7Ia7Ga2g5mNMLOPzOwDM7u4kn3MzHqY2SdmNtPM9s1mjSIisWvaFLp3B+CXX0K3e+fOcOihMHNmWJ1NZFWy3WIvAy53978BBwAXmNmeFfY5Htgt+uoIPJzdEkVEYlZaCoSu9vr14dVX4a67YMgQ2HrrWCuTPJDVYHf379x9WnR7PvARsH2F3U4GnvRgArC5mW2bzTpFRGJRUhIGoM+dC0CTQ4yvvja+OLeEK66A9XTyVJIQ2z8TM6sLNAAmVti0PfB1wv1Z/Dn8RUQKT0kJs752mh7mAJzT0pk319n+0ZJ465K8Ekuwm9nGwCDgEnefV3FzJU/xSl6jo5lNMbMps2fPzkSZIiJZ9corYVrYKVPC/aeegk03jbcmyT9ZD3Yzq0oI9f7u/mIlu8wCdki4Xxv4tuJO7t7b3Ru6e8NamhBZRPLYb79B165wyilQpw5MmwZ066ZpYWWtZPuqeAP6AB+5+72r2O1VoHV0dfwBwFx3/y5rRYqIZNFHH4UlVh98EC65JMz1vvvurDTcTSQV2Z5S9mCgFfCemZVGj10D7Ajg7r2AwcAJwCfAIqBtlmsUEck4d+jbFy66CDbaCF5/PSy1KrKushrs7j6Wys+hJ+7jwAXZqUhEJPvmzoVOneC558KqbE89BdttF3dVUig0eEJEJIsmTAhj0194AW69FYYNU6hLeinYRUSyoLwcbrsNmjQJ3fBjxsC//w1VqsRdmRQaLdsqIpJh334LrVvDO+9A8+ZhEZfNN4+7KilUCnYRkQwaPBjatIGFC+HRR6F9ezSMTTJKXfEiIhmwZAlcdlm40n277WDqVDjvPIW6ZJ5a7CIiafbxx9CiRZho5sILwwIu1avHXZUUCwW7iEgaPfkkdOkCG2wAL78MJ58cd0VSbNQVLyKSBvPmQatW4Xz6fvvBjBkKdYmHgl1EZB1Nngz77gvPPAPdu8Pw4VC7dtxVSbFSsIuIrKXycrj7bjjoIFi6FEaOhBtu0Nh0iZfOsYuIrIUffgjd7kOHwqmnwmOPwV/+EndVImqxi4ikbNiwsG76qFHw8MMwaJBCXXKHgl1EZFUqLJ26dCn8619w7LGw5Zbh3Pr552tsuuQWBbuISGWaNg1XwkU+/RQOOQTuvDOszDZ5Muy9d3zliayKzrGLiFSmtPT3m888E1rmVarAwIFwxhnxlSWyJmqxi4gkKikJfetz54b7Zpzd0rhvsxJKSxXqkvsU7CIiiUpKwJ3p+7YHYD1zrr/OafN5CXXqxFuaSDIU7CIiCdzhgQfggPcfA8JSqzfdBOvrxKXkCf1TFRGJzJ4NbdvCG2/AiSfCoj26cfjhcVclkhoFu4gIMGIEtGwJP/8MPXqEVdnMSuIuSyRl6ooXkaK2bBlcdx0ceSRsuilMnAhdu2psuuQvtdhFpGh98QWcfTaMHw/t2oWWeo0acVclsm4U7CJSlAYOhA4dwkIuzzwDLVrEXZFIeqgrXkSKyqJF0LEjNG8Of/1rmIdGoS6FRMEuIkXjvfegYUN49NEw5/vYsbDzznFXJZJeCnYRKXju8NBD0KgRzJkTVme7/XaoWjXuykTST+fYRaSg/fILtG8PL78Mxx0HTzwBW20Vd1UimaMWu4gUrNGjw7rpb7wB99wTvivUpdAp2EWk4JSVhSnfDz8cqleHcePgsstgPf3GkyKgrngRKShffx1mkBszBlq1gp49YZNN4q5KJHsU7CJSMF5+OUw0s2wZPPlkCHaRYrNWHVNmVsPMqqS7GBGRtbF4MVxwAZx6ahi+Nm2aQl2KV1LBbmbrmdnZZvaGmf0I/Bf4zsw+MLO7zGy3zJYpIlK5Dz+E/fcPw9kuuyycT99Nv5GkiCXbYh8B7AL8G9jG3Xdw962AQ4AJwO1mdk6GahQR+RN36N07TDjz/fcweHC48r1atbgrE4lXsufYj3L3ZRUfdPdfgEHAIDPTVA8ikhW//hqmhR04EI46KpxP33bbuKsSyQ1JtdgrC/W12UdEZF2NGwf168NLL4XZ44YOVaiLJFpjsJvZ0Wb2qJnVj+53zHhVIiIVLF8Ot9wChx4axqOPHRvme9fYdJGVJfNfogtwJXCOmR0B1M9oRSIiiUpK+PZbOPpouO46+Oc/Yfr0cMGciPxZMsE+291/dfcrgGOARhmuSUTkD927s88+MHEi9OkT1k7fbLO4ixLJXckE+xsrbrj71cCTmStHRCRYsgQu3qAXALVrw9SpYfIZs5gLE8lxawx2d39lxW0z+yfQL7p9nZm9aGb7Zq48ESlGP11YwgbVjQeWdgagdIaxx98sTAAvIquV6mUn17v7fDNrAhwLPAE8nP6yRKQYuUO/flC3Xwk1t3Re2+vqPza4K9hFkpBqsC+Pvv8DeDhqzWs6CBFZZ/PmhcVb2raFRo1gxgw48f3b4i5LJO+kGuzfmNkjQHNgsJltsBavISKykkmToEEDeP55uOkmePtt2H77aGO3brHWJpJvUg3l5sBQ4Dh3/xX4C2EonIhIysrL4c474eCDwxrqo0aFIW1VEpeYUve7SEqSmlLWzMyDRcCLKx539++A7xL3yUyZIlJovv8eWreGt96C00+HRx+FLbaIuyqR/Jf0IjBm1tXMdkx80MyqmdkRZvYE0Cb95YlIIRo6FOrVgzFjoFevMOe7Ql0kPZJdBOY4oB3wrJntBPwKVAeqAMOA+9y9NBMFikjhWLoUrr0W7r4b9toLhg8P30UkfZIKdnf/DXgIeChaxa0msDg6z540M+sLNAN+dPe9K9neFHgF+Dx66EV3vzGV9xCR3PTJJ9CiBUyZAp07hyVWN9ww7qpECk+yLfbfRau4fbeW79cPeJDVz143xt2breXri0gO6t8fzj8f1l8fBg2C006LuyKRwpXVoWruPhr4JZvvKSLxWbAA2rSBc84JS63OmKFQF8m0XByDfqCZzTCzIWa2yrNvZtbRzKaY2ZTZs2dnsz4RScK0abDvvvD003DDDTBiBOy445qfJyLrJteCfRpQx93rAf8BXl7Vju7e290bunvDWrVqZas+EVkDd7jvPjjgAFi0KFwg17176IYXkcxbY7Cb2dFm9qiZ1Y/ud8xUMe4+z90XRLcHA1XNrGam3k9E0uvHH6FZM7jsMjj++ND1fthhcVclUlyS+Ru6C9AWuM7M/gLUz1QxZrYN8IO7u5k1Jvzh8XOm3k9E0uedd8K59Dlz4MEHoUsXLbEqEodkgn12NKztCjO7HWi0tm9mZs8CTYGaZjYL6AZUBXD3XsAZQGczKwMWA2dpNjuR3LZsWZjO/fbb4a9/hTffDJPPiEg8kgn2N1bccPerzazr2r6Zu7dYw/YHCcPhRCQPfP45nH02TJgA550H998PNWrEXZVIcVtjsLv7K2Z2WeJjFe+7+73pLkxEcttzz0HH6IqbAQPgzDPjrUdEgmSvU90k+v5XQlf8q9H9E4HR6S5KRHLXwoVw8cXQp0+48v2ZZ2CnneKuSkRWSGq4m7t3d/fuhKlk93X3y939cmA/oHYmCxSRHBAtnTpjBjRsCH37wr//DaNHK9RFck2q49h3BJYm3F8K1E1bNSKSm7p358EHYf/9Ye7csNTqrbdC1apxFyYiFaU6ZcRTwCQzeym6fwrwRForEpGc8vOmO7El0LUrnHAC9OsHmhNKJHel1GJ391sIY9rnEOZ8b+vut2WiMBGJWUkJmLHl/C8AcIw3Bhu1epbEWZWIrMHaTCm7HChP+BKRAlNWBt28hCrrObvvFk0l4R6+ovPtIpKbUgp2M7sY6E+4iG4r4Ol1GdcuIrnnq6+gaVO48UZo1QqmTo27IhFJRarn2NsD+7v7QgAzuwMYT1iwRUTy3IsvQvv2ocX+9NPQsmW0oVu3WOsSkeSl2hVvhK74FZZHj4lIHlu8GDp3htNPh113henTE0Id1P0ukkdSbbE/DkyscFV8n7RWJCJZ9cEHYda4Dz6AK6+Em2+GatXirkpE1lZKwe7u95rZKOBgQku9rbtPz0hlIpJR7tC7N1xyCWy6aVi85dhj465KRNZVqi123H0qoMtpRPLYnDnQoQMMGgRHHw1PPgnbbBN3VSKSDkmdYzezsdH3+WY2L+FrvpnNy2yJIpJOY8eGZVVfeQXuvDO01BXqIoUjqRa7uzcxMwP2cvevMlyTiGTA8uVhGtiSEqhbF959Fxo3jrsqEUm3pK+Kd3cHXlrjjiKSc775Bo46Cm64IVwoN326Ql2kUKU63G2CmTXKSCUikhGvvgr77AOTJ4d53vv3DxfLiUhhSjXYDyeE+6dmNtPM3jOzmZkoTETWzW+/wUUXwcknQ506YQa5Nm3ANPOESEFL9ar44zNShYik1X//C2edFdZPv/hiuOMO2GCDuKsSkWxItcX+FXAI0MbdvwQc2DrtVYnIWnGHPn1gv/3CefXXXoP771eoixSTVIP9IeBAoEV0fz7QM60VichamTsXWrSA886D/fcPrfVmzeKuSkSyLdVg39/dLwB+A3D3OYAmnxSJ2cSJ0KABvPAC3HILvPUWbLdd3FWJSBxSDfZlZlaF0AWPmdVCa7KLZEclC7GUl8Ptt0OTJuH26NFwzTVQpUr2yxOR3JBqsPcgjGXfysxuAcYCt6a9KhFZWdOm0L37Sg99912Y2/3f/4ZTT4XSUjjooFiqE5EckuoiMP3NbCpwJGERmFPc/aOMVCYifygtXenukCFh6NqCBWEhl/PO0zA2EQmSCnYzqw6cD+wKvAc84u5lmSxMRAjd74kt9Si9J9KNbf5ewoABsOee8ZQmIrkp2a74J4CGhFA/Hrg7YxWJyB9KSsIYNncA9tvXMZzZXUqYOFGhLiJ/lmxX/J7u/ncAM+sDTMpcSSJS0VNPQSvg88/hpZfglFPirkhEclWywb5sxQ13LzOdzBPJivnzoUsXePppKNuxGzPGwg47xF2ViOSyZIO9XsK66wZsGN03wsJvWlJCJM2mTg3Twn72WeiRb31diYaxicgaJbseu36diGRJeTncd18Yxrb11jByJBxySNxViUi+SHURGBHJoB9+gHPPhTffDOfR+/SBv/wl7qpEJJ+kOkGNiGTIW29BvXowYgT07AkvvqhQF5HUKdhFYrZsGfzrX3DMMSHIJ08OF8zpGlURWRvqiheJ0WefhRXZJk2Cjh3DufWNNoq7KhHJZwp2kZgMGACdOoWW+fPPwz//GXdFIlII1BUvkmULF0K7dqGlvvfeYd10hbqIpIuCXSSLSkthv/2gXz+49loYNQrq1Im7KhEpJAp2kSxwhx49YP/9Yd48ePttuPlmWF8nw0QkzfRrRSTDfvoJ2raF11+HZs3g8cehZs24qxKRQqUWu0gGjRgRxqYPGwYPPACvvqpQF5HMUrCLZEBZGVx/PRx5JGy8MUyYABddpLHpIpJ56ooXSbMvv4Szz4Zx40IXfI8eIdxFRLJBwS6SRi+8AB06wPLl8MwzYUibiEg2qSteJA0WLQqTzfzzn7DbbjB9ukJdROKhYBdZR++9B40aQe/ecNVVMHYs7LJL3FWJSLFSsIukoqTk95vu8PDD0Lgx/PwzDB0Kd9wB1arFV56ISFaD3cz6mtmPZvb+KrabmfUws0/MbKaZ7ZvN+kTWqHt3AH75BU4/PazCdthhYVrYY46JuTYREbLfYu8HHLea7ccDu0VfHYGHs1CTSHI23xyAMWOgfv0w4czdd8PgwbD11rFWJiLyu6wGu7uPBn5ZzS4nA096MAHY3My2zU51IqtQUhIGoM+dC8AhhxpffW18cW4Jl18O6+mElojkkFz7lbQ98HXC/VnRYyLxKSnh66+cww51AM5p6cyb62zXuyTeukREKpFrwV7ZvFxe6Y5mHc1siplNmT17dobLkmL2yiuh633q1HD/6adh001jLUlEZJVyLdhnATsk3K8NfFvZju7e290bunvDWrVqZaU4KS6//QYXXginnAJ168K0aUC3bjFXJSKyerkW7K8CraOr4w8A5rr7d3EXJcXno4/CMLaePeGyy8L0sLvvzkrD3UREclFWp5Q1s2eBpkBNM5sFdAOqArh7L2AwcALwCbAIaJvN+kTcoU+fsGDLxhvDG2/ACSfEXZWISPKyGuzuvtpJNt3dgQuyVI7ISn79NUwL+/zzYVW2p56CbTUmQ0TyTK51xYvEYvx4aNAABg2C224L66cr1EUkHynYpagtXw633gqHHBLujx0LV1+tsekikr+0bKsUrW+/hVatYPhwaN4cHnnk98nlRETyloJditIbb8C554blVh97DNq1C5PLiYjkO3U4SlFZsgQuvRSaNYPttoMpU6B9e4W6iBQOtdilaHz8MZx1FkyfDl27wp13QvXqcVclIpJeCnYpeO7w5JNwwQWwwQZhitiTToq7KhGRzFBXvBS0efPgnHPC+fSGDcO66Qp1ESlkCnYpWJMnh7HpAwbAjTfCO+9A7dpxVyUiklkKdik45eVw111w0EGwbBmMGgXXXw9VqsRdmYhI5ukcuxSUH36ANm1g6FA47bQwlG2LLeKuSkQke9Ril4IxbBjss09ooffqBS+8oFAXkeKjYJe8t3QpXHUVHHss1KoVzq136qSx6SJSnBTskn8S1kT/9FNo0iScU+/UCSZNgr33jq80EZG4Kdgl/3TvDkD//uGq9//7v9Dt3qsXbLRRzLWJiMRMwS75pWlTIIxLP+eccE69tBROPz3OokREcoeCXfJDSUk4aT5qFAD9njAcY/QRJdSpE2tlIiI5RcEuecG7lXD/fc4G1RyAkSMc3FnvxpI4yxIRyTkKdsl5s2fDiSeGVdmOPTY8FvXIi4hIBQp2yWnDh0O9evDWW9CjR1jAhW7d4i5LRCRnKdglJy1bBtdcA0cdBZtuGoaxde0ajU1PGO4mIiIr05SyknO++AJatIAJE6B9e3jgAahRI+6qRETyg4JdcsrAgdChQ1hDfcAAOPPMuCsSEckv6oqXnLBoUQj05s1hjz3C2HSFuohI6hTsEruZM6FhQ+jTB66+GsaMgZ12irsqEZH8pGCX2LhDz57QuDHMmRNWZ7vtNqhaNe7KRETyl86xSyx+/jlcGPfKK3D88dCvH2y1VdxViYjkP7XYJetGj4b69WHwYLj3Xnj9dYW6iEi6KNgla8rKwhD0ww+H6tVh/Pgwm9x6+lcoIpI26oqXrPj6a2jZMlwY17o1PPggbLJJ3FWJiBQeBbtk3EsvhfPpy5bBU0+F5VZFRCQz1AkqGbN4MXTpAqedBjvvDNOnK9RFRDJNwS4Z8cEHYRjbww/D5ZfDuHGw665xVyUiUvjUFS9p5Q6PPgqXXBLOoQ8ZAscdF3dVIiLFQy12SZs5c8KUsJ06QZMmMGOGQl1EJNsU7JIW48aFsekvvwx33AFvvgnbbBN3VSIixUfBLutk+XK45RY49FCoUgXGjoWrrtLYdBGRuOjXr6ydkhK++QaOPhquuy50wU+fDvvvH3dhIiLFTcEua6d7d+rVg4kToW9f6N8fNtss7qJERETBLin57Te4eINeANSuDVOnQtu2YBZzYSIiAijYJQU/XVhC9Q2NB5Z2BqB0hrHH3yxMAC8iIjlBwS5r5A6PPw51Hi+h5pbO66/5HxvcFewiIjlEwS6rNXduWLylXbswk9yMGdCsWdxViYjIqijYZZUmToQGDeD55+Hmm+Htt2H77aON3brFWpuIiFROwS5/Ul4eJplp0iSMUx89Gq69NoxT/52630VEcpLmipeVfP89tGoVWuennx7mfd9ii7irEhGRZKnFLr97803YZx94913o3RsGDlSoi4jkGwW7sHQpXHEFHH88bL01TJkCHTpobLqISD7KerCb2XFm9j8z+8TMrq5ke1Mzm2tmpdHXDdmusZh88gkcdBDccw906QKTJsGee8ZdlYiIrK2snmM3sypAT+BoYBYw2cxedfcPK+w6xt01qCrDnn4aOneGqlXhxRfh1FPjrkhERNZVtlvsjYFP3P0zd18KDABOznINRW/+fGjdOlwkV78+lJYq1EVECkW2g3174OuE+7Oixyo60MxmmNkQM9urshcys45mNsXMpsyePTsTtRakqVNh333Doi3dusGIEbDjjnFXJSIi6ZLtYK/sciyvcH8aUMfd6wH/AV6u7IXcvbe7N3T3hrVq1UpvlQWovBzuvRcOPDAs5DJ8eBiKvr4GPIqIFJRsB/ssYIeE+7WBbxN3cPd57r4guj0YqGpmNbNXYuH58ccwDezll8MJJ4Su98MOi7sqERHJhGwH+2RgNzPbycyqAWcBrybuYGbbmIWBVmbWOKrx5yzXWTDefhvq1Qst9J494aWXYMst465KREQyJasdse5eZmYXAkOBKkBfd//AzM6PtvcCzgA6m1kZsBg4y90rdtfLGixbBjfcEKaG3WMPGDo0TD4jIiKFzQohMxs2bOhTpkyJu4yc8fnn0KJFWMSlQwe47z6oUSPuqkREJF3MbKq7N6xsmy6dKjADBkCnTmHWuOeeg+bN465IRESySVPKFoiFC6F9+9BS33PPcIGcQl1EpPgo2AtAaSnstx88/jhcc01YZrVu3birEhGROCjY85g7/Oc/sP/+MG9euAL+llvCFLEiIlKcdI49T/30E7RrB6+9Bv/4R2ita54eERFRiz0PjRwZxqYPHQr33x/CXaEuIiKgYM8rZWVw/fVwxBFh+Nr48XDxxVo3XURE/qCu+Dzx5ZfQsiW8+y6ce244t77xxnFXJSIiuUbBngcGDYLzzoPly8OqbGefHXdFIiKSq9QVn8MWL4bzz4czzoDddoPp0xXqIiKyegr2HPX++9CoETzyCFx5JYwdC7vsEndVIiKS69QVn2PcQ5hfeilsumm48v2YY+KuSkRE8oVa7Dnkl19Ct3vnznDooTBzpkJdRERSo2DPEWPHQv368OqrcNddMGQIbL113FWJiEi+UbDHbPlyuPFGOOwwqFYNxo2DK66A9fTJiIjIWtA59hjNmgXnnAOjRoUx6g89FM6ri4iIrC0Fe0xefRXatoUlS6BfP2jdWjPIiYjIulOHb5b99ht07Qonnwx16sC0adCmjUJdRETSQ8GeRR99FJZYffBBuOSSMNf77rvHXZWIiBQSdcVngTv07QsXXQQbbQSvvx6WWhUREUk3tdgzbO5caNEizPV+wAEwY4ZCXUREMkfBnkETJoSx6S+8ALfeCsOGwXbbxV2ViIgUMgV7BpSXw+23Q5MmoRt+zBj497+hSpW4KxMRkUKnc+xp9t130KoVvPMONG8e5n3ffPO4qxIRkWKhYE+jwYPD0LWFC+HRR6F9ew1jExGR7FJXfBosWQKXXRYuittuO5g6NVwsp1AXEZFsU4t9HX38cbjqfdo0uPDCsIBL9epxVyUiIsVKwb4OnnwSunSBDTaAl18Os8mJiIjESV3xa2H+/HCBXJs2sN9+YWy6Ql1ERHKBgj1FU6ZAgwbwzDPQvTsMHw61a8ddlYiISKBgT1J5Odx9Nxx4ICxdCiNHwg03aGy6iIjkFp1jT8IPP4Ru96FD4dRT4bHH4C9/ibsqERGRP1OLfQ2GDYN69WDUKHj4YRg0SKEuIiK5S8G+CkuXwr/+BcceC1tuCZMnw/nna2y6iIjkNnXFV+Kzz8LY9EmToFMnuPfesNyqiIhIrlOwV/DssyHMq1SBgQPhjDPirkhERCR56oqPLFgAbdvC2WfD3/8OpaUKdRERyT8KdmD69DDRzBNPwHXXhQvl6tSJuyoREZHUFXWwu8MDD8ABB4QW+zvvwE03wfo6QSEiInmqaCNs9uzQ9f7GG3DiidC3L9SsGXdVIiIi66YoW+wjRoSx6W+9BT16wCuvKNRFRKQwFFWwl5WFc+hHHgmbbgoTJ0LXrhqbLiIihaNouuK/+CJc8T5+PLRrF1rqNWrEXZWIiEh6FUWwDxwIHTqEhVyeeSZMPiMiIlKICrorftEi6NgRmjeHv/41jE1XqIuISCEr2GB/7z1o2BAefTTM+T52LOy8c9xViYiIZFbBBbs7PPQQNGoEc+aE1dluvx2qVo27MhERkczLerCb2XFm9j8z+8TMrq5ku5lZj2j7TDPbN9nX/uUXOP10uOACOPxwmDEDjj46vfWLiIjksqwGu5lVAXoCxwN7Ai3MbM8Kux0P7BZ9dQQeTua1R48OY9Nffx3uuSdMPLPVVmksXkREJA9k+6r4xsAn7v4ZgJkNAE4GPkzY52TgSXd3YIKZbW5m27r7d6t60W+/DS30nXeGcePCuXUREZFilO2u+O2BrxPuz4oeS3WflXz3XRijPm2aQl1ERIpbtlvslc3x5muxD2bWkdBVD7Dk6aft/aefXsfq8k9N4Ke4i4iBjru46LiLi447OatcgzTbwT4L2CHhfm3g27XYB3fvDfQGMLMp7l50bXUdd3HRcRcXHXdxSedxZ7srfjKwm5ntZGbVgLOAVyvs8yrQOro6/gBg7urOr4uIiMgfstpid/cyM7sQGApUAfq6+wdmdn60vRcwGDgB+ARYBLTNZo0iIiL5LOtzxbv7YEJ4Jz7WK+G2Axek+LK901BaPtJxFxcdd3HRcReXtB23hRwVERGRQlBwU8qKiIgUs7wK9kxOR5vLkjjuPcxsvJktMbMr4qgxE5I47pbR5zzTzMaZWb046ky3JI775OiYS81sipk1iaPOdFvTcSfs18jMlpvZGdmsL1OS+Lybmtnc6PMuNbMb4qgz3ZL5vKNjLzWzD8xsVLZrzIQkPu8rEz7r96N/639J6U3cPS++CBfbfQrsDFQDZgB7VtjnBGAIYSz8AcDEuOvO0nFvBTQCbgGuiLvmLB73QcAW0e3ji+jz3pg/TqPtA/w37rqzcdwJ+w0nXKdzRtx1Z+nzbgq8HnetMRz35oRZSXeM7m8Vd93ZOO4K+58IDE/1ffKpxf77dLTuvhRYMR1tot+no3X3CcDmZrZttgtNszUet7v/6O6TgWVxFJghyRz3OHefE92dQJjzIN8lc9wLPPpfD9Sgkgmc8lAy/78BugKDgB+zWVwGJXvchSaZ4z4beNHdv4Lwey7LNWZCqp93C+DZVN8kn4I9I9PR5oFCPKZkpHrc7Qm9NfkuqeM2s1PN7L/AG0C7LNWWSWs8bjPbHjgV6EXhSPbf+YFmNsPMhpjZXtkpLaOSOe7dgS3MbKSZTTWz1lmrLnOS/r1mZhsBxxH+kE1J1oe7rYO0TUebZwrxmJKR9HGb2eGEYC+Ec81JHbe7vwS8ZGaHAjcBR2W6sAxL5rjvB/7l7svNKts9LyVz3NOAOu6+wMxOAF4mrH6Zz5I57vWB/YAjgQ2B8WY2wd0/znRxGZTK7/MTgXfd/ZdU3ySfgj1t09HmmUI8pmQkddxmtg/wGHC8u/+cpdoyKaXP291Hm9kuZlbT3fN5fu1kjrshMCAK9ZrACWZW5u4vZ6XCzFjjcbv7vITbg83soSL5vGcBP7n7QmChmY0G6gH5HOyp/P8+i7Xohgfy6uK59YHPgJ3446KDvSrs8w9WvnhuUtx1Z+O4E/YtoXAunkvm896RMEPhQXHXm+Xj3pU/Lp7bF/hmxf18/Url33m0fz8K4+K5ZD7vbRI+78bAV8XweQN/A96J9t0IeB/YO+7aM33c0X6bAb8ANdbmffKmxe5FOh1tMsdtZtsAU4BNgXIzu4RwpeW8Vb1urkvy874B2BJ4KGrFlXmeLx6R5HGfTlhPYRmwGDjTo98G+SrJ4y44SR73GUBnMysjfN5nFcPn7e4fmdmbwEygHHjM3d+Pr+p1l8K/81OBYR56K1KmmedEREQKSD5dFS8iIiJroGAXEREpIAp2ERGRAqJgFxERKSAKdhERkQKiYBfJsmi1phUrNw2Mpo5c19c818weXMM+dc0s7cOFMvi6Tc3soIT7/ZJZ0c3MNjSzUWZWJYX3utDM8n54rAgo2EXisNjd67v73sBS4Py4C8pRTQkr+KWqHWHxkOUpPKcvcNFavJdIzlGwi8RrDLCrmdUws75mNtnMppvZyfB7S/xFM3vTzP7PzO5c8UQza2tmH0frVB+c8PhKLVszW1DxTSu28M3sdTNrumJ/M7sjWnjjbTNrHC3E8ZmZnbS6gzGzKmZ2V3QcM82sU/R40+g1XjCz/5pZf4tmFTKzE6LHxppZj6iWuoQ/eC6NejcOid7iUDMbF9WyqtZ7S+CVhPcdZWbPRz+r282spZlNMrP3zGwXAHdfBHxhZo1Xd3wi+UDBLhITM1ufsI78e8C1hHWXGwGHA3eZWY1o1/rAmcDfgTPNbAcLyxF3JwT60cCeaSytBjDS3fcD5gM3R+9xKnDjGp7bHpgbHUcjoIOZ7RRtawBcEtW6M3CwmVUHHiHM9d8EqAXg7l8QVnG7L+rdGBO9xraExX6aAbdXfHMzqwbsHD1/hXrAxYSfXytgd3dvTFhjoGvCflOAQxDJc3kzpaxIAdnQzEqj22OAPsA44CQzuyJ6vDphLnyAd9x9LoCZfQjUISyCMtLdZ0ePP0dY5jIdlgJvRrffA5a4+zIzew+ou4bnHgPsk9Ca3oywEtlSwtoNs6J6S6PXWgB85u6fR/s/C3Rczeu/7O7lwIdmtnUl22sCv1Z4bLK7fxe976fAsIRjOzxhvx+BPdZwfCI5T8Eukn2L3b1+4gNRt/Tp7v6/Co/vDyxJeGg5f/y/XdV80GVEvXHR61Zb3T6R6gm3lyXMRV6+4v3dvTzqZVgdA7q6+9AKx9F0FceR6vqria9R2XMXs/KxVHxOecL9clb+HVg9er5IXlNXvEhuGAp0TTjv3GAN+08EmprZlmZWFfhnwrYvCOtYA5wMVK3k+V8A9c1sPTPbgbBqWDoMJSxYUhXAzHZPOKVQmf8CO0fn1CGcclhhPrBJKm/u7nOAKlEXf6p2J6wgJpLXFOwiueEmQgDPjIaO3bS6naOu5RJgPPA2MC1h86PAYWY2CdgfqGyFqHeBzwnd0XdXeP66eAz4EJgWHccjrKZn0N0XA12AN81sLPADMDfa/BpwaoWL55IxjHAePlUHE36WInlNq7uJSKzMbGN3XxD1VvQE/s/d71uH12sAXOburTL5HJFcpRa7iMStQ3Qx3QeEi+0eWZcXc/fpwIhUJqghXHR3/bq8r0iuUItdRESkgKjFLiIiUkAU7CIiIgVEwS4iIlJAFOwiIiIFRMEuIiJSQBTsIiIiBeT/AYWdkIXEalb3AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# \n", "#\n", "T2error = np.zeros(nLen)\n", "T2error = 2*Tmean*Terror\n", "#\n", "lenError = np.zeros(nLen)\n", "lenError[0:nLen] = 0.002\n", "np.set_printoptions(precision = 4)\n", "print(\" \")\n", "print(\"Error in T^2 in s^2\\n\",T2error)\n", "#\n", "plt.figure(figsize = (8, 6))\n", "#\n", "# Plot the data\n", "plt.errorbar(lengths, Tmean**2, xerr = lenError, yerr = T2error, color = 'r', \\\n", " linestyle = '', marker = '+') \n", "#\n", "# Draw a straight line\n", "xPoints = np.zeros(2)\n", "xPoints[1] = 1.1*np.max(lengths)\n", "plt.plot(xPoints, gradient*xPoints + intercept, color = 'b', linestyle = '-')\n", "#\n", "plt.xlim(0, 1.2*np.max(lengths))\n", "plt.ylim(0, 1.2*np.max(Tmean**2))\n", "plt.title('Data with hand fit', fontsize = 14)\n", "plt.xlabel('Pendulum length (m)')\n", "plt.ylabel('Period$^2$ (s$^2$)')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Exercise 8\n", "\n", "For how many of your points does the line go through the error bar? You should expect that 2/3 of the points lie within 1 error bar of the line. The probability that a point is 3 error bars from the line should be less than 0.3% if everything is working out correctly." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "Now that you have a set of data with error bars you can conduct a fit rather than estimate a line of best fit by eye. The method you will use makes use of a statistic called $\\chi^2$, which, when minimised, tells you the line of best fit and gives you the uncertainties on the patameters describing the line. (Here, we use a linear relationship between $T^2$ and $l$, so the parameters are the gradient and intercept of the line.) We can also use the $\\chi^2$ value to tell us if the model we are using is a good fit to the data.\n", "\n", "The $\\chi^2$ statistic for a linear hypothesis, with errors on the $x$ and $y$ measurements, is given by the formula:\n", "\n", "$$\n", "{\\chi ^2} = \\sum\\limits_{i = 1}^N {\\frac{{{{\\left( {{y_i} - (m{x_i} + c)} \\right)}^2}}}{{\\Delta {y_i}^2 + {m^2}\\Delta {x_i}^2}}} \n", "$$\n", "\n", "To determine the quality of the fit, we look at the calue of the statistic $\\chi^2$ per degree of freedom ($\\chi^2$/NDF). The number of degrees of freedom (NDF) is the number of data points minus the number of free parameters in the fit, which in the case of a straight line is 2 (the gradient $m$ and the intercept $c$). A good fit is indicated by values of $\\chi^2$/NDF between 0.25 and 4. Any higher and the error bars are too small or the model is wrong, any lower and the error bars are too large and we cannot reliably judge how well model agrees with the data." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Exercise 9\n", "\n", "Use the cells below to carry out a $\\chi^2$ fit to your pendulum data. (Your values of $l$, $\\Delta l$, $T^2$ and $\\Delta T^2$ should already be available as NumPy arrays in the cells above. Running the cells here will allow you to perform the fit to your data.)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " \n", "Fit quality:\n", "chisq per point = \n", " [0.49 0.08 0.46 2.12 0.16 0.64 0.45 0.23 0.91 0.02]\n", "chisq = 5.568, chisq/NDF = 0.696.\n", " \n", "Parameters returned by fit:\n", "Intercept = 0.0129 +- 0.0114\n", "Gradient = 3.9923 +- 0.0323\n", " \n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfEAAAGHCAYAAACgSWuhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAtPklEQVR4nO3deZyd493H8c9PhNhaKpSK9Xm0j/JIaMTeppuiFKW2EmJfYqmltNSuVJU+KKmidmopImKJIAuCSUxSW1HViqUiKonsyVzPH9dJjZFlJpkz91k+79drXplzzj1zfndP4tvruu/fdUVKCUmSVH2WKLoASZK0aAxxSZKqlCEuSVKVMsQlSapShrgkSVXKEJckqUoZ4pIkVSlDXFJFiIgtI+LpiBgaEbdFROeia5IqnSEuqVL8A/hWSukbwBvALgXXI1U8Q1zSAkXEixHRewGvvxkR31nc90kpvZNSmlZ6OBtoms/7fSUino+IyRFxbGtqlGqVIS4VpBR+00ph9FFEPBURR0REq/5dtld4LkxKacOU0hMd9Z4RsS6wAzBwPof8FHgipbRCSumyImqUKoUhLhVr55TSCsDawIXAKcC1xZZUnIj4HHADsH9KaeZ8DlsbeLHjqpIqlyEuVYCU0sSU0gBgL+CAiNgIICJOjYi/lUbrL0XEbqXnbwLWAu6PiI8j4qcLOr6liOgbEfc3e/x6RNzR7PFbEdGj9P2bEfGd+b1nSY+IGBsREyPiTxHRZT7v+9eI+EtErFZ6vFHpZzaIiCWB24CzUkp/nc/PPwZ8E7iiVMOX21CjVHMMcamCpJSeBcYB25ae+lvp+88DZwM3R8TqKaX9gX+SR/LLp5QuWtDx83irocC2EbFE6fXOwNYAEbEesDwwtkVt83tPgD2B7YF1gY2BA+dzij2Aj4GdS6F9A/DLlNLLwD7A5sAZEfFEROw1j/99vgUMB/qVani1DTVKNccQlyrPO8AXAFJKd5Zu+GpKKf0JeA3oNb8fbO3xKaU3gMnkUP0G8DDwdkT8T+nx8JTSPG8sm4/LSu/7IXB/6ffOq75ppff6X+B08g1sF5deuyml1DWl1Lv09ac2vL9Ul5YsugBJn7EG8CFARPQBTgDWKb22PNB1fj/YxuOHAr2B/y59/xE5wLcsPW6L95p9PxX40gKOfQG4CFgF6JVSmtPG95JU4khcqiARsRk5xEdExNrAH4B+wMoppRXJARilw1OLn13Y8S3NDfFtS98PJYf4N5h/iKf5PN8WL5Gn3c8tTaO3t/aoUaoKhrhUASLicxGxE3A7cHNK6S/AcuRAGl86pi+wUbMf+xewXrPHCzu+paHkm8SWSSmNI19r3h5YGXh+Pj/T8j0XRZ/Sn39czN8zP+1Ro1QVDHGpWPdHxGTgLeA04BKgL0BK6SXgN8DT5GD6X+DJZj97AXB6qcf8pFYc/ymlm8I+Joc3KaVJ5JXSnlzAFPen3rOtJxsRvYCjyDfvbdDWn2+lxapRqiaRkjNPksovIpYGRgP9yTfbjUwp/a7YqqTq5khcUkc5mzxDcAW5fe37EbFUsSVJ1c2RuKSyK92w9wjQI6X0j9LSqg8DE1JKWxZbnVS9DHFJkqqU0+mSJFUpQ1ySpCpliEuSVKWqbtnVrl27pnXWWafoMiRJ6jCjRo36IKW0Ssvnqy7E11lnHRoaGoouQ5KkDhMR/5jX806nS5JUpQxxSZKqlCEuSVKVqrpr4vMya9Ysxo0bx/Tp04supey6dOlCt27d6Ny5c9GlSJIKVhMhPm7cOFZYYQXWWWcdIua3dXL1SykxYcIExo0bx7rrrlt0OZKkgtXEdPr06dNZeeWVazrAASKClVdeuS5mHCRJC1cTIQ7UfIDPVS/nKUlauJoJ8aJ16tSJHj16/OfrzTffZKuttgLgzTff5NZbby24QklSramJa+KVYJlllqGxsfFTzz311FPAJyG+7777FlCZJKlWORIvo+WXXx6AU089leHDh9OjRw8uvfTSgquSJNWKmhuJH388tBgQL7YePeC3v13wMdOmTaNHjx4ArLvuutxzzz3/ee3CCy/k4osvZuDAge1bmCSprtVciBdlXtPpkiSVU82F+MJGzJIk1QqviXeAFVZYgcmTJxddhiSpxhjiHWDjjTdmySWXpHv37t7YJklqNzU3nV6Ujz/+eL7Pde7cmSFDhnR0SZKkGudIXJKkKmWIS5JUpQxxSZKqlCEuSVKVMsQlSapS9RvivXvnL0mSqlT9hng7m7sV6YYbbkj37t255JJLaGpqWuDPuEWpJGlxGOLtZO7a6S+++CKDBw9m0KBBnH322Qv8GUNckipHStDQUHQVbWOIl8Gqq67K1VdfzRVXXEFKiTfffJNtt92WTTfdlE033fQ/+4y33KJ0fsdJksrro49gr72gVy947rmiq2m9+lmxreX177k7jrV8/okn2uXt1ltvPZqamnj//fdZddVVGTx4MF26dOG1115jn332oaGh4TNblE6dOnWex0mSyuepp2DffeHtt+HCC+FrXyu6otarnxAvQEoJgFmzZtGvXz8aGxvp1KkTr7766jyPb+1xkqTFN2dODu0zz4S114Ynn8wj8WpSPyHecoQ9dwTeTiPvlt544w06derEqquuytlnn80Xv/hFxowZQ1NTE126dJnnz1x66aWtOk6StHjefhv22y9HwD77QP/+8LnPFV1V23lNvAzGjx/PEUccQb9+/YgIJk6cyOqrr84SSyzBTTfdxJw5c4DPblE6v+MkSe1nwADYeON87fv66+GWW6ozwMEQbzfTpk37T4vZd77zHbbbbjvOPPNMAI466ihuuOEGtthiC1599VWWW2454LNblM7vOEnS4ps+HY45BnbZJU+fjx4NBxwAEUVXtuhi7nXbatGzZ8/U8mavl19+mQ022KBtv6jM0+nltEjnK0l17OWXYe+9YexY+MlP4IILYOmli66q9SJiVEqpZ8vn6+eaeEtVGN6SpLZJCa67Do49FpZbDh54AHbcseiq2o/T6ZKkmvTRR3n0fcghsOWWMGZMbQU4GOKSpBr01FPQowf8+c+5jeyRR2D11Yuuqv3VTIhX27X9RVUv5ylJi2LOHDj/fPj612GJJWDECDjllPx9LaqJ0+rSpQsTJkyo+YBLKTFhwgT7xyVpHt5+G777XTj9dNhzT3j+edh886KrKq+auLGtW7dujBs3jvHjxxddStl16dKFbt26FV2GJFWU+++Hvn1h2jT44x+rv3WstWoixDt37sy6665bdBmSpA42fXqeLr/ssnwN/Pbb4StfKbqqjlO26fSIWDMiHo+IlyPixYg4bh7H9I6IiRHRWPo6o1z1SJJqyyuvwBZb5AA//ngYObK+AhzKOxKfDZyYUhodESsAoyJicErppRbHDU8p7VTGOiRJNaR57/eyy8LAgfD97xddVTHKNhJPKb2bUhpd+n4y8DKwRrneT5JU+z76KG9YMrf3e+zY+g1w6KC70yNiHWAT4Jl5vLxlRIyJiAcjYsP5/PxhEdEQEQ31cPOaJOmznn4aNtkE7rorL5taq73fbVH2EI+I5YG7geNTSpNavDwaWDul1B24HLh3Xr8jpXR1SqlnSqnnKqusUtZ6JUmVZc4c+OUvYdtt8+MRI+DUU2u397styvo/QUR0Jgf4LSmlP7d8PaU0KaX0cen7QUDniOhazpokSdXjnXdy7/dpp8Eee0BjY76ZTVk5704P4Frg5ZTSJfM5ZrXScUREr1I9E8pVkySpegwcmPf9fuaZfCPbbbfB5z9fdFWVpZx3p28N7A/8JSIaS8/9HFgLIKXUH9gDODIiZgPTgL1TrS+7JklaoBkz4Kc/rd/e77YoW4inlEYAC1wvJ6V0BXBFuWqQJFWXV17Jd583NsJxx8GvflVd+353tJpYsU2SVN1SysulHnNM7v2+/37YyRVEFsp7+yRJhZo4MY++Dz4437Q2ZowB3lqGuCSpMCNH5uved92V28geeQS+9KWiq6oehrgkqcM1feObXLDe1WyzTX48YgT87GfQqVOxdVUbr4lLkjrUO+9AnyfPZMic3uy1F/z+97aOLSpH4pKkDvPAA9C9Ozw9ZzOuXaafvd+LyRCXJJXdjBl5u9CddoI11oBRmx3JQb1eIBbYiKyFcTpdklRWf/0r7L137v0+9tjc+92ly41Fl1UTDHFJUlmkBNdfD/36wTLLwIABsPPORVdVW5xOlyS1u4kTYd994aCDYPPNc++3Ad7+DHFJUrt65pm87/edd8L558Pgwfk6uNqfIS5JahdNTXDhhbDNNvn74cPh5z+397ucvCYuSVps774L++8PQ4bAnnvm3u8VVyy6qtpniEuSFsugQXDAATB1KlxzTb4ObutYx3A6XZK0SGbMgJ/8BL7//bzeeUND3sTEAO84jsQlSW326qu59/v55/P2oRddBF26FF1V/THEJUmtlhLccEPu/e7Sxd7vojmdLklqlUmT4Mc/hr59YbPN7P2uBIa4JGmhnnkm7/t9xx1w3nnw6KP2flcCQ1ySNF9NTXmt87m938OGwWmn2ftdKbwmLkmap3ffhT598qj7Rz+Cq6+297vSGOKSpM+Y2/s9ZQr84Q+2jlUqp9MlSf8xYwaccMInvd+jRsEhhxjglcqRuCQJyL3f++wDo0fnFrJf/9re70pniEtSnUsJbrwRjj4all4a7rsPfvCDoqtSazidLkl1bNIk2G8/OPDA3Ps9dqwBXk0McUmqU88+m/f9/tOf4Nxz7f2uRoa4JNWZpqa81vnWW8xm9tvvMXQonH66vd/VyBCXpDry3nvwve/BKafArp0G0rj05my9ddFVaVEZ4pJUJx58EDbeGJ58Mi/ccseyB7JSTCy6LC0GQ1ySatyMGXDiibDjjrDaannf70MPhdikR14QXVXLFjNJqmEL7P1+4okiS1M7MMQlqQalBDfdBEcdlXu/770Xdtml6KrU3pxOl6QaM2kS7L9/Xvu8Z8+877cBXpsMcUmqIXN7v2+7Dc45B4YMgW7diq5K5WKIS1IN+E/v99Ywe3be9/sXv7D3u9Z5TVySqtx77+V9vwcPht13z1uHrrRS0VWpIzgSl6Qq9tBD0L07jBiRe7/vvNMAryeGuCRVoZkz4aSTYIcd4ItfbNb77b7fdcXpdEmqMq+9lnu/R43KLWQXXwzLLFN0VSqCIS5JVeTGG3NwL7UU3HMP7Lpr0RWpSE6nS1IVmDz5k97vr30t934b4DLEJanCNTTk3u9bb8293489BmuuWXRVqgSGuCRVqKamfL17yy3zjWxDh9r7rU/zmrgkVaD33stT5488Yu+35s+RuCRVmIcfzr3fw4ZB//72fmv+DHFJqhAzZ8LJJ8P228Oqq+Zr4Ycfbu+35q9sIR4Ra0bE4xHxckS8GBHHzeOYiIjLIuL1iBgbEZuWqx5Jqki9e0Pv3rz2Gmy1Vb4GftRReSOTDTcsujhVunJeE58NnJhSGh0RKwCjImJwSumlZsfsAKxf+tocuKr0pyTVh8ZGbpq5F0dtCp07w5//DLvtVnRRqhZlG4mnlN5NKY0ufT8ZeBlYo8VhuwA3pmwksGJErF6umiSpkkyeDH2m9qfPtN+zySa599sAV1t0yDXxiFgH2AR4psVLawBvNXs8js8GvSTVnIYG2HRTuGXO3px9Njz+uL3faruyh3hELA/cDRyfUprU8uV5/Eiax+84LCIaIqJh/Pjx5ShTkjpEUxP85jf5+veMGbn3+4wz7P3WoilriEdEZ3KA35JS+vM8DhkHNP//nt2Ad1oelFK6OqXUM6XUc5VVVilPsZJUZv/6F+y4Y959bKedoLERttmm6KpUzcp5d3oA1wIvp5Qumc9hA4A+pbvUtwAmppTeLVdNklSUhx+GjTfOI++rroK774YvfKHoqlTtynl3+tbA/sBfIqKx9NzPgbUAUkr9gUHAjsDrwFSgbxnrkaQON3MmnHZabh3bcEMYMgQ22qjoqlQryhbiKaURzPuad/NjEnB0uWqQpCK9/nre97uhAY48Ml8Ld99vtSfXTpekMrjlFjjiCHu/VV4uuypJ7Wjy5LxxyX775e1DGxsNcJWPIS5J7WTUqNz7ffPNcNZZed/vtdYquirVMkNckhbT3N7vLbeE6dPzwi1nnglLesFSZeZfMUlaDP/6Fxx4IDz0UJ42v+YaW8fUcRyJS9IieuSRvO/3E0/Y+61iGOKS1EYzZ8Ipp8D3vgddu8Jzz+U70d33Wx3N6XRJaoO//S33fs8N7t/8BpZdtuiqVK8McUlqpVtuyYu2dOoEd90Fu+9edEWqd06nS9JCNO/97t497/ttgKsSGOKStACjRsHXvpZ7v888M7eP2futSmGIS9I8NDXBJZfk3u9p0/LCLWedZe+3Kot/HSWphfffz9PnDz0Eu+6ae79XXrnoqqTPciQuSb175y9g8OC87/fjj8OVV+bNSwxwVSpH4pLU2MistCSnnwIXXQRf/WoO8//936ILkxbMkbikuve3OeuwzZSHuOgiOPzw3ANugKsaOBKXVNduvRWOmPYknaKJO++EPfYouiKp9RyJS6pLH3+cNy758Y9h4y2Wo/H1FQxwVR1DXFLdGT067/t9001wxhl5A5O11y66KqntDHFJdSMluPRS2GILmDo1936ffba936pe/tWVVBfefx/69oVBg2CXXeDaa20dU/VzJC6p5j36aF7zfMgQ+N3v4J57DHDVBkNcUs2aNQtOPRW22w5WWgmefRaOOsp9v1U7nE6XVJPeeCPv+/3ss3DYYflauPt+q9YY4pJqzm235UVbOnXC3m/VNKfTJdWMjz/ON6/tu29e/7yx0QBXbTPEJdWE55/P+37fcAP84hf2fqs+GOKSqlpK8Nvf5t7vKVNy7/c559j7rfrgX3NJVWv8+Lx06qBB8IMf5N7vrl2LrkrqOI7EJVWlIUPyde8hQ+CKK+Deew1w1R9DXFJVmTULfvYz+O53P+n9Pvpoe79Vn5xOl1Q1/v733Pv9zDNw6KG593u55YquSiqOIS6pKtx+e+79joA77oAf/ajoiqTiOZ0uqaJNmQIHHZRH4BttBGPGGODSXIa4pIo1t/f7+uvh9NNh6FB7v6XmDHFJFScl+L//y73fkyfnO9DPPdfeb6kl/0lIqijjx+elUx94AHbeGa67ztYxaX4ciUuqGEOG5H2/H30ULr8c7rvPAJcWxBCXVLhZs+DnP8+935//fG4h69fP3m9pYQxxSYX6+9/h61+HCy6Ag5e8gYaGPBqXtHCGuKTC/OlP0KMHvPwy/GmZA/nDsse5eIvUBoa4pA43ZQocfDDsvTd89at53+89e72ZE11Sq3l3uqQO1diYw/vVV+G00+DMM6FzZ/IG4JLaxJG4pA6RElx2GWy++Se93+edVwpwSYvEkbiksrP3WyoPR+KSyuqxx/Ld5oMH55G4vd9S+zHEJZXFrFn5mvd3vgOf+1ze9/uYY+z9ltqT0+mS2t3f/w777gsjR+a70P/v/9z3WyqHso3EI+K6iHg/Il6Yz+u9I2JiRDSWvs4oVy2SOs7c3u+XXsp7gF9zjQEulUs5p9OvB7ZfyDHDU0o9Sl/nlLEWSWU2ZQoccsine7/32qvoqqTaVrYQTykNAz4s1++XVDkaG/O+39ddl9dAHzYM1l236Kqk2lf0jW1bRsSYiHgwIjac30ERcVhENEREw/jx4zuyPkkL0Lz3e9KkvPvY+efb+y11lCJDfDSwdkqpO3A5cO/8DkwpXZ1S6plS6rnKKqt0VH2SFuCDD2CXXeC442C77WDsWPjWt4quSqovhYV4SmlSSunj0veDgM4RYfeoVAUefzz3fj/8cL7zfMAAe7+lIhQW4hGxWkTuGI2IXqVaJhRVj6SFmzULTj8dvv1tWGGFvO/3scfa+y0VpWx94hFxG9Ab6BoR44Azgc4AKaX+wB7AkRExG5gG7J1SSuWqR9LiefPN3Pv99NNw0EH5WritY1KxFhriEdEPuCWl9O+2/OKU0j4Lef0K4Iq2/E5JxbjjDjjssHwj22235TYyScVrzXT6asBzEXFHRGw/dwpcUu2bMgUOPTT3e2+wwSfbiEqqDAsN8ZTS6cD6wLXAgcBrEfHLiPivMtcmqUBjxkDPnnDttfZ+S5WqVTe2la5Vv1f6mg2sBNwVEReVsTZJBUgJrrgi935PnJh3H7P3W6pMrbkmfixwAPABcA1wckppVkQsAbwG/LS8JUrqKB98kG9au/9+2HFHuP56cGkGqXK15u70rsAPU0r/aP5kSqkpInYqT1mSOtrjj8N+336HD9LK/Pa3S9s6JlWB1lwTP6NlgDd77eX2L0lSR5o9G37xi9z7vXxMYeTy3+W44wxwqRoUvXa6pAK9+SZ8/etw3nnQty+Mmrg+m0weVnRZklqpbIu9SKpsd96Z28dSgltvhX0WuLKDpErkSFyqM1On5oVb9twT/ud/4PnnDXCpWhniUh0ZOzb3fl9zDZx6KgwfDuutV3RVkhaVIS7VgZTgd7+DXr3g3/+GRx6BCy6w91uqdl4Tl2rchAm593vAAHu/pVrjSFyqYU88kff9fvBBuPRSGDjQAJdqiSEu1aC5vd/f+lbeLnTkSDj+eHu/pVrjdLpUY/7xj7zv91NPwYEHwuWXw/LLF12VpHIwxKUactddufd7zhy45ZYc5pJql9PpUg2YOhUOPxx+9CP48pfzvt8GuFT7DHGpys3t/b76ajjlFBgxwt5vqV4Y4lKVmlfv94UX2vst1ROviUtVaMIEOPhguO8+2GGH3Pu96qpFVyWpozkSl6rM0KG593vQoE96vw1wqT4Z4lKVmD0bzjwz934vu+wnvd9L+K9YqltOp0tV4B//gB//GJ58Eg44AK64wt5vSYa4VHl6985/PvEEAHffDYccknu/b745h7kkgSEuVZ7GRiD3fv/kJ7l1bLPN4Lbb4L/+q9jSJFUWQ1yqQH+Z81X23gxeein3fp9zDiy1VNFVSao0hrhUQVKCq375ESecACtOyL3f3/1u0VVJqlSGuFQhPvww937fey9svz3ccIOtY5IWzOYUqQIMG5Z7vx94AH7zm/ynAS5pYQxxqUBze7+/+U3o0gWefhpOOMHeb0mt43S6VJB//jO3i40YkXu/L78cVlih6KokVRNDXCrA3N7v2bPt/Za06Jy0kzrQ1KlwxBGwxx6w/vq5JdwAl7SoDHGpg7zwQt429Pe/h5NPztPoLt4iaXEY4lKZpQRXXZVXXfvgA3j4YbjoIhdvkbT4DHGpjD78EHbfHY46Cr7xDRgzBrbbruiqJNUKQ1wqk7m93wMHwsUX5/2/v/jFoquSVEsMcamdzZ4NZ531Se/3U0/BiSfa+y2p/dliJrWjf/4T9tsPhg+HPn3yvt/2fksqF0Ncaif33JPXPp81C266KYe5JJWTE3zSYpo2Ld+49sMf5pax5583wCV1DENcWgwvvJBbx666Ck46CZ58Ev77v4uuSlK9MMSlRZAS9O+fA3z8eHjoIfj1r+39ltSxDHGpjT78MC+beuSRufd77Fj43veKrkpSPTLEpTYYPhx69ID777f3W1LxDHGpFWbPhrPPht69Yeml7f2WVBlsMZMW4q238k5jw4fnu86vvNLeb0mVoWzjiIi4LiLej4gX5vN6RMRlEfF6RIyNiE3LVYu0qO65Jy+d+vzzcOONuf/bAJdUKco5GXg9sP0CXt8BWL/0dRhwVRlrkdpk2udX46ilr+GHP4T11oPRo2H//YuuSpI+rWwhnlIaBny4gEN2AW5M2UhgxYhYvVz1SK314ovQ6+MhXDXzEE46KV//Xn/9oquSpM8q8racNYC3mj0eV3pOKkRK8PvfQ8+e8H7XDe39llTxigzxmMdzaZ4HRhwWEQ0R0TB+/Pgyl6V6NLf3+4gj4Otfz/t+2/stqdIVGeLjgDWbPe4GvDOvA1NKV6eUeqaUeq6yyiodUpzqx4gRufd7wAC46CJ48EFYbbWiq5KkhSsyxAcAfUp3qW8BTEwpvVtgPaozc+bAOefkVdeWWipf+z75ZHu/JVWPsvWJR8RtQG+ga0SMA84EOgOklPoDg4AdgdeBqUDfctUitTRuXO79HjYs/3nllfC5zxVdlSS1TdlCPKW0z0JeT8DR5Xp/aX7uvTfv+z1jRu79tnVMUrVy4lB1Y9o0OPpo2G03WGedvICLAS6pmhniqgsvvgi9euVp8xNPhKeftvdbUvVz7XTVtJTgD3+A44/Py6U++CBsv6B1BCWpijgSV83697/hRz+Cww+HbbbJvd8GuKRaYoirJj35ZO79vu++3Pv90EP2fkuqPYa4asqcOXDuuXnVtSWXzGFu77ekWuU1cdWMcePyft9Dh8K++8JVV9n7Lam2GeKqCffdBwcdlHu/b7ght47FvFbnl6Qa4iSjqtq0adCvH+y6a+79Hj0a+vQxwCXVB0NcVeull2DzzeF3v4MTTshrn3/5y0VXJUkdx+l0VZ3mvd/LLw8PPAA77lh0VZLU8RyJq6r8+9+w556593vrrXPvtwEuqV4Z4qoaTz2Ve7/vvRd+9St4+GFYffWiq5Kk4hjiqnhz5sB553269/unP7X3W5K8Jq6K1rz3e599oH9/e78laS7HMqpY9210Gt3XnUhDA1x/PdxyiwEuSc0Z4qo406fDMcfAri+ez9pNbzJqFBxwgL3fktSSIa6K8vLLuff7iivgJz+Bp6d25ytfKboqSapMXhNXRUgJrr0Wjj3W3m9Jai1H4ircRx/BXnvBoYfCVlvZ+y1JrWWIq1Bze7/vuQcuvBAeecTeb0lqLUNchZgzB84/P/d+L7EEjBgBp5xi77cktYXXxNXh3n47934/8UTu/b7qKvj854uuSpKqjyGuDjVgAPTtm9vI/vhHW8ckaXE4eakOMX16vvN8l11grbXyvt8HHmiAS9LiMMRVdnN7vy+/PG8fOnIk9n5LUjtwOl1lkxJcd10egS+7LAwcCN//ftFVSVLtcCSusvjoI9h7bzjkENhySxg71gCXpPZmiKvdze39vvtuuOACe78lqVwMcbWb5r3fEbn3+9RT7f2WpHLxmrjaxdtvw/77w+OP5yVUf/97e78lqdwMcS22gQNzu9i0aflGNlvHJKljONGpRTZ9Ohx3HOy8M6y5Zu797tvXAJekjmKIa5G88gpssQVcdlkOcnu/JanjOZ2uNmnZ+33//bDTTkVXJUn1yZG4Wu2jj/KGJYcckldgGzPGAJekIhniapWnn4ZNNoG77oJf/hIGD4YvfanoqiSpvhniWqA5c/KCLdtumx+PGAE/+xl06lRsXZIkr4lrAd55J/d+P/aYvd+SVIkciWueBg6EjTeGkY9N4dpl+nHbbQa4JFUaQ1yfMmPGJ73f3brBqM2O5KBeL9j7LUkVyOl0/ccrr+S7zxsbcwvZr34FXbrcWHRZkqT5MMRFSvDHP8Ixx8Ayy8CAAXkkLkmqbE6n17mJE2HffeHggz/p/TbAJak6GOJ1bOTIvO/3nXfmLUQHD4Y11ii6KklSaxnidaipKfd+b7NNnkofPhx+/nN7vyWp2nhNvM688w706QNDhsCee+be7xVXLLoqSdKiMMTryAMP5L2+p0yBa66Bgw5y21BJqmZlnU6PiO0j4q8R8XpEnDqP13tHxMSIaCx9nVHOeurVjBlw/PF5s5IvfQlGjco3shngklTdyjYSj4hOwO+A7wLjgOciYkBK6aUWhw5PKbkXVpn89a+w99659/uYY+Cii6BLl6KrkiS1h3KOxHsBr6eU3kgpzQRuB3Yp4/upmbm935tuCm+9lXu/L7vMAJekWlLOEF8DeKvZ43Gl51raMiLGRMSDEbHhvH5RRBwWEQ0R0TB+/Phy1FpTJk6EH/84X/Pu1cveb0mqVeUM8XldcU0tHo8G1k4pdQcuB+6d1y9KKV2dUuqZUuq5yiqrtG+VNeaZZ/K+33fcAeedB48+au+3JNWqcob4OGDNZo+7Ae80PyClNCml9HHp+0FA54joWsaaalZTE1x4Ye79bmqCYcPgtNPs/ZakWlbOEH8OWD8i1o2IpYC9gQHND4iI1SLyPdIR0atUz4Qy1lST3n0XttsOfvYz2G23fBPbVlsVXZUkqdzKdnd6Sml2RPQDHgY6AdellF6MiCNKr/cH9gCOjIjZwDRg75RSyyl3LcCgQXDAAbn3+w9/sHVMkupJVFtm9uzZMzU0NBRdRuFmzMgj70svhY03httvhw02KLoqSVI5RMSolFLPls+7YlsVevXV3Pv9/PPQrx/8+te2jklSPTLEq0hKcMMNObiXXhruuw9+8IOiq5IkFcVdzKrEpEm597tvX+jZE8aONcAlqd4Z4lXgmWfyvt933AHnnpt3ILP3W5JkiFewpib41a9y7/ecOTB0KJx+ur3fkqTMa+IV6t13877fjz4Ke+wBV18NK61UdFWSpEpiiFegBx/Mvd8ff5zD+5BD7P2WJH2W0+kVZMYMOOEE2HFHWG01aGiAQw81wCVJ8+ZIvEK8+irssw+MHm3vtySpdQzxgqUEN94IRx+de7/vvRd2cdd1SVIrOJ1eoEmTYL/94MADc+/3mDEGuCSp9Qzxgjz7bN73+/bb4Zxzcu93t25FVyVJqiaGeAdraoKLLoKtt4bZs/O+37/4hb3fkqS285p4B3rvvdz7PXgw7L573jrU3m9J0qJyJN5BHnwwbxk6YkTu/b7zTgNckrR4DPEymzEDTjwx935/8Yv2fkuS2o/T6WXUvPf7qKPg4othmWWKrkqSVCsM8TJICW66KQf3UkvBPffArrsWXZUkqdY4nd7OJk2C/ffPa59/7Wu599sAlySVgyHejp57DjbdFG67Lfd+P/YYrLlm0VVJkmqVId4OmpryWudbbQUzZ+Z9v+39liSVm9fEF1Pz3u8f/hCuucbWMUlSx3Akvhgeegi6d4fhw6F/f7jrLgNcktRxDPFFMHMmnHQS7LADrLpq7v0+/HB7vyVJHcvp9DZ67bXc+z1qlL3fkqRiGeJtMLf3u3Nn+POfYbfdiq5IklTPnE5vhcmTc+93nz55+9AxYwxwSVLxDPGFaGjIwX3rrXD22fD44/Z+S5IqgyE+H01N+Xr3llvmG9meeALOOMPeb0lS5fCa+Dy8915eNvWRR/K0+TXXwBe+UHRVkiR9miPxFh5+OPd+DxsGV10Fd99tgEuSKpMhXjJzJpx8Mmy/PayySl4H/Ygj7P2WJFUup9OB11/Pvd8NDXDkkfCb39j7LUmqfHUf4jffnIPb3m9JUrWp2+n0yZNz3/f++0OPHtDYaIBLkqpLXYZ4Q0Pe9/uWW+DMM3Pv91prFV2VJEltU1ch3tSUr3dvtRVMn57D+6yzYMm6v6ggSapGdRNf//pX7v1++GF7vyVJtaEuRuKPPJJ7v4cOtfdbklQ7ajrEZ86En/4Uvvc96NrV3m9JUm2p2en05r3fhx8Ol1wCyy5bdFWSJLWfmgzxW27JI+4ll4S77oLddy+6IkmS2l9NTadPnpxvXttvv9z7PWaMAS5Jql01E+KjRuXe75tvtvdbklQfqj7E5/Z+b7ll7v1+7DF7vyVJ9aGqo+5f/4IDD4SHHoJdd8293yuvXHRVkiR1jKodiQ8enHu/H38crrwyb15igEuS6klZQzwito+Iv0bE6xFx6jxej4i4rPT62IjYdGG/MyU45RTYbrsc2s89l3chs/dbklRvyjadHhGdgN8B3wXGAc9FxICU0kvNDtsBWL/0tTlwVenP+XrlFRg92t5vSZLKORLvBbyeUnojpTQTuB3YpcUxuwA3pmwksGJErL6gXzpjBtx5J/Tvb4BLkupbOUN8DeCtZo/HlZ5r6zGf8tWvwh57tEt9kiRVtXLenT6vq9RpEY4hIg4DDis9nBERLyxmbdWoK/BB0UUUwPOuL553ffG8W2/teT1ZzhAfB6zZ7HE34J1FOIaU0tXA1QAR0ZBS6tm+pVY+z7u+eN71xfOuL+153uWcTn8OWD8i1o2IpYC9gQEtjhkA9Cndpb4FMDGl9G4Za5IkqWaUbSSeUpodEf2Ah4FOwHUppRcj4ojS6/2BQcCOwOvAVKBvueqRJKnWlHXFtpTSIHJQN3+uf7PvE3B0G3/t1e1QWjXyvOuL511fPO/60m7nHTlHJUlStanaZVclSap3FRvi5ViytRq04rz/JyKejogZEXFSETWWQyvO+8elz3lsRDwVEd2LqLO9teK8dymdc2NENETENkXU2d4Wdt7NjtssIuZERE2sDtGKz7t3REwsfd6NEXFGEXW2t9Z83qVzb4yIFyNiaEfXWA6t+LxPbvZZv1D6u/6FNr1JSqnivsg3wv0NWA9YChgDfLXFMTsCD5J7zbcAnim67g4671WBzYDzgZOKrrkDz3srYKXS9zvU0ee9PJ9c9toYeKXoujvivJsd9xj5vpo9iq67gz7v3sDAomst4LxXBF4C1io9XrXoujvivFscvzPwWFvfp1JH4mVZsrUKLPS8U0rvp5SeA2YVUWCZtOa8n0op/bv0cCR5TYFq15rz/jiV/oUDyzGPxZCqUGv+fQMcA9wNvN+RxZVRa8+71rTmvPcF/pxS+ifk/851cI3l0NbPex/gtra+SaWGeFmWbK0CtXhOrdHW8z6YPAtT7Vp13hGxW0S8AjwAHNRBtZXTQs87ItYAdgP6Uzta+/d8y4gYExEPRsSGHVNaWbXmvL8MrBQRT0TEqIjo02HVlU+r/7sWEcsC25P/T2ublLXFbDG025KtVaYWz6k1Wn3eEfFNcojXwrXhVp13Suke4J6I+DpwLvCdchdWZq05798Cp6SU5kTt7DPcmvMeDaydUvo4InYE7iXv8ljNWnPeSwJfA74NLAM8HREjU0qvlru4MmrLf893Bp5MKX3Y1jep1BBvtyVbq0wtnlNrtOq8I2Jj4Bpgh5TShA6qrZza9HmnlIZFxH9FRNeUUjWvN92a8+4J3F4K8K7AjhExO6V0b4dUWB4LPe+U0qRm3w+KiCvr5PMeB3yQUpoCTImIYUB3oJpDvC3/vvdmEabSgYq9sW1J4A1gXT65IWDDFsd8n0/f2PZs0XV3xHk3O/YsaufGttZ83muRV/bbquh6O/i8/5tPbmzbFHh77uNq/WrL3/PS8ddTGze2tebzXq3Z590L+Gc9fN7ABsCQ0rHLAi8AGxVde7nPu3Tc54EPgeUW5X0qciSe6nTJ1tacd0SsBjQAnwOaIuJ48h2Pk+b3eytdKz/vM4CVgStLo7PZqco3Tmjlee9O3l9gFjAN2CuV/uVXq1aed81p5XnvARwZEbPJn/fe9fB5p5RejoiHgLFAE3BNSqmqd6tsw9/z3YBHUp6FaDNXbJMkqUpV6t3pkiRpIQxxSZKqlCEuSVKVMsQlSapShrgkSVXKEJckqUoZ4pIkVSlDXNIClfb0HhsRXSJiudJ+zxsVXZckF3uR1AoRcR7Qhbw5xbiU0gUFlyQJQ1xSK0TEUsBzwHTy+vVzCi5JEk6nS2qdLwDLAyuQR+SSKoAjcUkLFREDgNvJOzKtnlLqV3BJkqjc/cQlVYiI6EPeNe7WiOgEPBUR30opPVZ0bVK9cyQuSVKV8pq4JElVyhCXJKlKGeKSJFUpQ1ySpCpliEuSVKUMcUmSqpQhLklSlTLEJUmqUv8PYid5wZvVGNUAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# \n", "\n", "'''\n", "Define the fit and error functions, carry out the fit, test whether it has worked and print out \n", "and plot the results if it has. The routine used to carry out the fit is optimize.leastsq\n", "(see: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html).\n", "This allows the definition of an error function (called a cost function in the documentation)\n", "that makes it possible to cope with errors in both x and y.\n", "'''\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.optimize import leastsq\n", "%matplotlib inline\n", "def fitLine(p, x):\n", " '''\n", " Straight line\n", " '''\n", " f = p[0] + p[1]*x\n", " return f\n", "def fitError(p, x, y, xerr, yerr):\n", " '''\n", " Error function for straight line fit\n", " '''\n", " e = (y - fitLine(p, x))/(np.sqrt(yerr**2 + p[1]**2*xerr**2))\n", " return e\n", "#\n", "# Set initial values of fit parameters, run fit\n", "pInit = [1.0, 1.0]\n", "out = leastsq(fitError, pInit, args=(lengths, Tmean**2, lenError, T2error), \\\n", " full_output = 1)\n", "#\n", "# Get output\n", "pFinal = out[0]\n", "covar = out[1]\n", "#\n", "# Test if fit failed (i.e. returned nan = not a number)\n", "if np.isnan(pFinal[0]):\n", " print(\" \")\n", " print(\"Fit failed\")\n", " print(\" \")\n", " print(\"pfinal\\n\",pFinal)\n", " print(\" \")\n", " print(\"covar\\n\",covar)\n", "else:\n", " #\n", " # calculate chis**2 per point, summed chi**2 and chi**2/NDF\n", " chiarr = fitError(pFinal, lengths, Tmean**2, lenError, T2error)**2\n", " chisq = np.sum(chiarr)\n", " NDF = nLen - 2\n", " redchisq = chisq/NDF\n", " #\n", " print(\" \")\n", " print(\"Fit quality:\")\n", " np.set_printoptions(precision = 2)\n", " print(\"chisq per point = \\n\",chiarr)\n", " print(\"chisq = {:5.3f}, chisq/NDF = {:5.3f}.\".format(chisq, redchisq))\n", " #\n", " cf = pFinal[0]\n", " mf = pFinal[1]\n", " #\n", " cErr = np.sqrt(covar[0, 0])\n", " mErr = np.sqrt(covar[1, 1])\n", " #\n", " print(\" \")\n", " print(\"Parameters returned by fit:\")\n", " print(\"Intercept = {:5.4f} +- {:5.4f}\".format(cf, cErr))\n", " print(\"Gradient = {:5.4f} +- {:5.4f}\".format(mf, mErr))\n", " print(\" \")\n", " #\n", " # Plot data\n", " fig = plt.figure(figsize = (8, 6))\n", " plt.plot(xPoints, fitLine(pFinal, xPoints), color = 'b',\\\n", " linestyle = '-', label = \"Fit\")\n", " plt.errorbar(lengths, Tmean**2, xerr = lenError, yerr = T2error, fmt='r', \\\n", " linestyle = '', label = \"Data\") \n", " plt.title('Data with $\\chi^2$ fit')\n", " plt.xlabel('x')\n", " plt.ylabel('y')\n", " plt.xlim(0.0, 1.2*np.max(lengths))\n", " plt.ylim(0.0, 1.2*np.max(Tmean**2))\n", " plt.legend(loc = 2)\n", " plt.savefig(\"LineFitPlot.png\")\n", " plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Exercise 10\n", "\n", "The fit tells you the gradient and intercept of the line of best fit and the errors on these parameters. Use these values to calculate a new value of $g$ and estimate its error." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " \n", "Acceleration due to gravity 9.89 +/- 0.08 m/s^2\n" ] } ], "source": [ "# \n", "\n", "gravAccel = (4*np.pi**2)/mf\n", "errAccel = (4*np.pi**2)/mf**2*mErr\n", "print(\" \")\n", "print(\"Acceleration due to gravity {:5.2f} +/- {:5.2f} m/s^2\".\\\n", " format(gravAccel, errAccel))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Exercise 11\n", "\n", "The fit tells you the value of the statistic $\\chi^2$/NDF. What does this tells you about the quality of your fit to the pendulum data?\n", "\n", "How does your value of $g$ compare with the accepted value of 9.81 ms$^{-2}$? Is it consistent?\n", "\n", "What do you think are the main reasons for any discrepancy?\n", "\n", "What do you think are the main systematic and random errors?\n", "\n", "How could the experimental method be changed to improve your estimate of $g$?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Group exercise\n", "\n", "Your measurement of $g$ and its error will be used in the group statistics exercise, so please provide your result on the group sheets at the front of the optics lab. " ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.8.8" } }, "nbformat": 4, "nbformat_minor": 4 }