{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "\n", "In this lecture, we will briefly discuss the calculation of the error of some \"experimental\" measurements and look at comparing different theoretical predictions to experimental data. We will also briefly revise some Monte Carlo methods.\n", "\n", "## The simple pendulum\n", "\n", "Figure 1 shows a schematic diagram of a pendulum of length $l$ with a bob of mass $m$. \n", "The bob is at an angle $\\theta$, and an arc-length $s$, from its equilibrium position. The combination of the gravitiational force, $mg$, \n", "and the tension in the supporting (massless) rod, $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 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$ (the phase constant) 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}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The pendulum at large amplitude\n", "\n", "If large ampitude oscillations are considered, the pendulum is much more difficult to describe mathematically than is the case for small amplitude oscillations. In fact, there is no algebraic expression for $T$ that is valid at large angles, and measurements can deviate noticeably from the predictions of the above simple harmonic motion (SHM) formula. [This paper](https://arxiv.org/pdf/physics/0510206.pdf), by Lima and Arun, gives more accurate (but still approximate) results for the period at large amplitude. Their formula (referred to as the LA formula in the following) is:\n", "\n", "$$\n", "T = - 2\\pi \\sqrt {\\frac{l}{g}} \\frac{{\\log \\left[ {\\cos \\frac{{{\\theta _0}}}{2}} \\right]}}{{1 - \\cos \\frac{{{\\theta _0}}}{2}}}.\n", "$$\n", "\n", "Here, $\\theta_0$, the maximum displacement of the pendulum, is assumed to occur at time $t = 0$.\n", "\n", "The SHM and LA expectations for the period of a pendulum of length 0.5 m are plotted below as a function of amplitude." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfgAAAGECAYAAADayDLFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA+oklEQVR4nO3dd5hURfb/8fchSzCsGEEFQVYkKiisEb+6RlbdFbMYF8wBdV1zXF1z+IkBVhHDCgpiABGzgAIKKCJBFBGVVSQoKAxImPP7o+5oOw4zPUP33O7bn9fz9DPdN57qnukzVbdulbk7IiIikiw14g5AREREMk8JXkREJIGU4EVERBJICV5ERCSBlOBFREQSSAleREQkgZTgpWCYWTczczNrvJ7HmWZm12UorHTPuaWZvWpmy80s9ntbzWyumV0SdxwAZnadmX0XfbanxB1PaaXfq2y+d2Z2iZnNzcaxJf8owUvOMbOB0Ze1m9lqM5tjZneYWYP1PPQ4YCtgcQbCrG6XAFsDHQllqBZR8pxWxqpdgQeqK451MbO2wLXAmYT35el4I0rLb9676Pe8R4zxSELVijsAkXV4HegJ1Ab2Ah4GGgBnVeVgZlbb3VcB8zMWYfVqCUx298/iDgTA3RfGHUOkZfTzec+TUbty6L2ThFMNXnLVz+4+392/dvengP8CRwBYcKmZfW5mK8zsYzM7sWRHM2sW1YqOM7M3zWwFcEZZTfRm9rdo/5/N7Gszu9LMLGX95mb2QnSeL83stIoCN7MW0T7zoyb1D8yse6lt/mZmU6Pjfm9mo81si3Ucby5wOHBSFP/AaPnvan5lNAe7mfU2syFRLHNS36tom63N7L9mttjMisxsipntGzV3Xwu0SWlROWUd59nWzJ4zs5+ixzAza5qy/rro0sax0ef2k5k9X9HlEjNrZ2avp7xPA81so5JjAs9FmxaXd+nCzG4xs1nRceaa2W1mVq+M+E6O1i8zs0fNrI6ZnR39biw2s7vMrEbKfnOjfZ+M9plvFTS/p7539mtz+pDo/Z2bGk+p/U4xs2Wlll0anXOZmT0ONCzjfKea2QwzW2lmn5pZn9QySHLpQ5Z8sYJQmwf4F3A6cA6wE/BvoJ+ZHVpqn38TmkJ3Ap4vfUAz6wQMAYYB7YDLgMuBc1M2G0ioJe5P+AfjJKBZBbE2BF4G/gx0AJ4FhpnZjtF5twQGA48BrYG9gSfKOd6uhBaNZwjN0BdUcP7SrgFeiGJ5GhhgZttFsTQARkdl+ivhfbgh2u9p4E5gVnTeMpvAo3+Inge2AP4P2JdwOeH51H+WonMcE53nAGBn4KZ1BW1m9YFRwDJgt2i/3YEB0SZ3AL2i5yXxrcty4DTC+302cCxwZaltmhH+keoOHAkcRXjfdo3i/TtwXhRHqouAmcAuhH+Ibjazv5UTS6pdo5+9ovh3LWfb3zCzowl/C9dG554VxZK6TS/gZsLvQGvgYuCfhPdAks7d9dAjpx6EpDoi5fVuwCJCcmlASPZ7ldrnHmBk9LwZ4MDFpbbpFi1vHL3+L/BmqW2uA+ZFz1tF2++Rsn47YC1wXSXLNAG4Knq+S3Tc7Sqx/whgYKllDvQotWwucEmpbf6d8roWUAScGL3uBfxU8p6Ucd7rgGllLP/lPIR/ZNYCzVLWbw8UA/unHGclsFHKNlcCs8spcy9gKdCojM+wZfS6R/gaq/Tv2Jmp547iW1EqvqHAQqBOyrK3gb6l3ofXSh37YeCdcj6Tsj6j0p/j79534BRgWcrrccB/Sm3zOjA35fVXQM9S21wIzKjse6ZH/j1Ug5dcdVDU7LgSGA+MIdSedgLqAaOi9cuiZsuzgBaljjGpgnO0Bt4ttewdoImZbRitLwbeL1np7l8C35R3UDNrEDUBzzCzH6L4OgPbRpt8RPginmZmz5rZWWa2WQWxro+pJU/cfQ0haW0eLdoZmOrui9bj+K2Bb9x9bsp55hDep51StvvS3ZemvP4mJY51HXequ/+Usmwc4TPZqexdymZmPczsnZLmbOBufv08SnxVKr7vgE899N1IXVY65vFlvK5UfFXUeh3nBiD6ndqG0LqV+rdyC7//W5EEUic7yVVjgN7AakLyWA1gZs2j9X8h1E5SrS71enkF5zBC7aksHq2vijuAgwg93z8j1JgfB+oAuPtaMzsA6Epo+j0d+LeZ7ePuH1XiPGXFWLuM7Uq/L86vl+eqWsZUFb2P6cSxPsctl5l1JVwSuR7oAywBDiN8TqnKiq+sZTXTPfd6KCa9z7Y8Je/tmYR/jKTAqAYvuarI3We7+5clyT0yA/iZ0Lw9u9Tjy0qeYwawZ6llexKa6H8iXFetQcp1UTPblnB9uTx7Ao+7+7PuPhWYR6kakwfj3f366PjfEK5PV8ZCUq47W+ikV9lb6D4A2pfT2W0VFSe0GYRWj2YpsWxPeJ9mVDKe0sftYGaNUpbtTvhMZlbiOHsA/3P3G919ooc7EbZbj7hK61rG68rEt5rfv8cLgS1K9WHoWGqbmes4NwDu/h3wP6BFGX8rsysRn+Qp1eAlr7j7T2Z2B3BH9OU3htCprStQ7O79K3G4O4GJUW/spwiJ9mLgiuhcs8xsFKGJszfhGu1d0c/yfAr81cxeIHx5X0u4rAD8UqPcH3iF0OS7M6EptbLJ8E3gHDMbR7gGfjPhOndlPEXoXPi8mV1O+GekHfCTu79FuF68nZntQmgx+cndfy51jNcJlx3+a2bnE2qe9xH+eXizkvGk+i+h1v24mV0DbAL0A4ZVMkF9SvgH5ARCE/aBwHHrEVdpXaP3biihj8BJwAmV2H8usJ+ZjSbcPfID4Vr/H4ArzGxwdNzS98rfS3hvJkbb9wC6AN+nbHMdcJ+ZLQFGEloBdgGauPu/KxGj5CHV4CUfXU344roEmA68Ruj1/EVlDuLuHxB6Sh8JTCNcm7wF6Juy2SnRcd8EhhMS4twKDn0RsAAYS+hNPyF6XmIpoVY5gtCEfydwo7s/WZn4Cf+MzCF8uQ8ldO5aUJkDuPtyYB9CTW844f28nl+bwJ8lJIY3CLXK3yVGd3fCHQYLo1jeIow3cES0rkrcvYiQjDck9IN4gZCgK7xVsdRxhgO3EzpiTiV0CrymqnGV4S6gPfAhoVf7Ne4+tBL7X0y48+Dr6Bi4+0xCv5LeKTHfnLqTuz9N+Du4KdqvXRRL6jYPE96vnoR/wsZGx6zU34rkJ1uPvz8RkYIW3bfe191LX88XiZ1q8CIiIgmkBC8iIpJAaqIXERFJINXgRUREEkgJXkREJIESdR9848aNvVmzZhk73mfff8YOf9ghY8eLW5LKk6SygMqTy5JUFlB5cllVyjJ58uRF7l72UNdxD4afyUenTp08kzr1y+zx4pak8iSpLO4qTy5LUlncVZ5cVpWyAJNck82IiIgUDiV4ERGRBFKCFxERSaBEdbIry+rVq5k3bx4rV1Z2Dg64reNtzJxZmUmhclu2ylOvXj2aNm1K7dqVnc1SRESyJfEJft68eTRq1IhmzZrx25kXK+YLndabtc5SZNUvG+VxdxYvXsy8efNo3rx5xTuIiEi1SHwT/cqVK9l0000rndwlPWbGpptuWqUWEhERyZ7EJ3hAyT3L9P6KiOSegkjwcWvYsOE6111wwQU0adKE4uLidW5z3HHH0b59e+6+++5shFemU045haFDKzOltYiI5JLEX4PPZcXFxTz33HNss802jBkzhm7duv1um/nz5zNu3Di+/PLLtI+7Zs0aatXSRysiUshUg4/RW2+9Rdu2bTnrrLMYNGhQmdsccMABLFiwgI4dOzJ27FimTJlC165dad++PX/961/54YcfAOjWrRtXXHEF++yzD/feey/dunWjT58+7L333rRu3ZqJEydywSkXsMMOO3DVVVcBMHfuXNq2bfvLue644w6uu+6638Vwww03sOuuu9K2bVt69+6NawZCEZGcV1DVvAsvhClT0t++aPV21K/gzq+OHeGee6oWz6BBgzjuuOM4/PDDueKKK1i9evXvbjV78cUX6d69O1OiwNu3b899993HPvvswzXXXMP111/PPVEAS5YsYfTo0QAMHz6cOnXqMGbMGO69914OP/xwBr06iK47dKVFixb06dMn7TjPPfdcrrnmGgB69uzJiBEj+Mtf/lK1QouISLVQDT4mq1atYuTIkRxxxBFsuOGGdOnShVdffbXcfZYuXcqSJUvYZ599ADj55JMZM2bML+uPOeaY32x/2GGHAdCuXTvatGnDZltsRt26ddl+++35+uuv0471rbfeokuXLrRr144333yT6dOnp72viIgEw4bB8uXVd76CqsFXtqY9Y+GX7LTZTlmJZdSoUSxdupR27doBUFRURP369Tn00EOrfMwGDRr85nXdunUBqFGjxi/PS16XXKdP7dxX1q1uK1eu5Oyzz2bSpElss802XHfddbolTkSkkt57D448Eq6/HqIG0axTDT4mgwYN4uGHH2bu3LnMnTuXL774gldffZWioqJ17rPRRhuxySabMHbsWACeeOKJX2rzVbHFFluwYMECFi9ezM8//8yIESN+t01JMm/cuDHLli1Tz3oRkUoqLobzz4ctt4RKXB1dbwVVg49LUVERTZs2/eX12WefzSuvvEK/fv1+WdagQQP23HNPhg8f/rum9lSPPfYYZ555JkVFRWy//fY8+uijVY6rdu3aXHPNNXTp0oXmzZuz4447/m6bjTfemF69etGuXTuaNWvGrrvuWuXziYgUoieegPffh4EDoVGjajzxuuaRzcdHWfPBz5gxo9Lz65aYvmB6lffNRdksz/q8z1WRpDmg3VWeXJaksrirPNXtxx/dt9zSfbfd3NeuLX/bTM8Hrxq8iIhIltx0E8yfD88/DzWq+aK4rsGLiIhkwezZcPfdcPLJ0KVL9Z9fCV5ERCQLLroI6tSBf/87nvOriV5ERCTDXnkFhg+HW2+FrbaKJwbV4EVERDJo9eowcmrLlnDBBfHFoRq8iIhIBvXtC598Ai++CCljjFU71eCrwU033USbNm1o3749HTt25L333qNbt25MmjTpl21SJ355++23MTMeeeSRX9Z/+OGHmBl33HFHtccvIiLpWbAArrsODjwQunePNxbV4LNs/PjxjBgxgg8++IC6deuyaNEiVq1aVeF+7dq14+mnn+b0008HYPDgwXTo0CHb4YqIyHq48kooKgpDo5vFG4tq8Fn27bff0rhx41/Ggm/cuDFbb711hfttu+22rFy5ku+++w53Z9SoURx88MHZDldERKro/ffhkUfgvPOgjIFBq11h1eArOV/sdquLoHb98jeqYL7YAw44gBtuuIFWrVqx//77c8wxx/wyfvwJJ5zABhtsAITZ5WqUGgWhR48eDBkyhJ133plddtnlNxPGiIhI7li7Fs45B7bYIjTR5wLV4LOsYcOGTJ48mf79+7PZZptxzDHHMHDgQAD++9//MmXKFKZMmcLIkSN/t+/RRx/NkCFDfpk3XkREctPDD8OkSXDnnbDhhnFHExRWDb6S88V+uXBGRqaLrVmzJt26daNbt260a9eOxx57LK39ttxyS2rXrs1rr73Gvffey7hx49Y7FhERyaxFi+Dyy2GffSCX6mKFleBjMGvWLGrUqMEOO+wAwJQpU9huu+2YNm1aWvvfcMMNLFiwgJo1a2YzTBERqaLLL4cff4T774+/Y10qJfgsW7ZsGeeddx5LliyhVq1atGzZkv79+9OjR4+09t99992zHKGIiFTVhAmhef7ii6FNm7ij+S0l+Czr1KlTmU3rb7/99m9eN2vW7JdafUlzfmnX5UrPDRER+aVj3dZbw7XXxh3N7ynBi4iIVEH//vDBBzBoEDRqFHc0v6de9CIiIpW0cCFccQXsuy8cc0zc0ZRNCV5ERKSSLrsMli0L487nUse6VAWR4N097hASTe+viBSS8eNhwIAw3/tO638nddYkPsHXq1ePxYsXKwllibuzePFi6tWrF3coIiJZt2YNnH02NGkCV18ddzTlS3wnu6ZNmzJv3jwWLlxY6X3n/zQfW5SjbS9VkK3y1KtXj6ZNm2b8uCIiueahh8KI5888Aw0bxh1N+RKf4GvXrk3z5s2rtG/P/j2Z1HtSxRvmiaSVR0SkOn33HVx1Fey/P6Q5lEmsEt9ELyIikgmXXBKmgr3vvtztWJdKCV5ERKQCr78OTz4Zes/nwlSw6VCCFxERKceKFXDWWdCyZbj3PV8k/hq8iIjI+rj5Zpg9G157DfLphqGs1eDNbICZLTCzdU6bZmbdzGyKmU03s9Epy+ea2cfROvUKExGRWMyYAbfeCieeGDrX5ZNs1uAHAn2Bx8taaWYbAw8AB7n7V2a2ealN9nX3RVmMT0REZJ2Ki+HMM8PtcHfeGXc0lZe1BO/uY8ysWTmbHA8Mc/evou0XZCsWERGRyho4EMaODdPBbl66CpoHLJsjvEUJfoS7ty1j3T1AbaAN0Ai4190fj9Z9AfwAONDP3fuXc47eQG+AepvW69Tm5sxNyDtz0UxaN26dsePFLUnlSVJZQOXJZUkqC6g86Vr94yZMv3YoG2z9Oa0uPgOrkf3RUKtSlslnTJ7s7p3LXOnuWXsAzYBp61jXF5gANAAaA58BraJ1W0c/Nwc+AvZO53ydOnXyTOrUL7PHi1uSypOksrirPLksSWVxV3nS1bOne+3a7tOnZ+XwZapKWYBJvo6cGOdtcvOAUe6+3MO19jFABwB3/yb6uQB4DtgttihFRKSgvPEGPPEEXHppbk8mU5E4E/wLwF5mVsvM6gNdgJlm1sDMGgGYWQPgAGCdPfFFREQyZeXK0LGuRQu48sq4o1k/WetkZ2aDgG5AYzObB1xLuOaOuz/k7jPNbBQwFSgGHnb3aWa2PfCchXEAawFPufuobMUpIiJSouSe91dfhQ02iDua9ZPNXvTHpbHN7cDtpZbNIWqqFxERqS6ffAK33AInnAB//nPc0aw/DVUrIiIFr7gYzjgj3PN+111xR5MZGqpWREQKXv/+MGZM/t7zXhbV4EVEpKB9/XXoMb/ffnDaaXFHkzlK8CIiUrDcw0xxa9eGWnw+zPOeLjXRi4hIwRo0CF56Ce6+G7bfPu5oMks1eBERKUgLF8L550OXLnDeeXFHk3lK8CIiUpAuuAB+/BEeeQRq1ow7msxTghcRkYIzfHhonr/qKmiTuTnKcooSvIiIFJSlS0PHurZt4bLL4o4me9TJTkRECso//wnffgvDhkGdOnFHkz2qwYuISMF4+23o1w/69IHdEj5PqRK8iIgUhKIi6NUr3A53ww1xR5N9aqIXEZGCcN11Yaa4N9+E+vXjjib7VIMXEZHEmzgR7rwz1OD33TfuaKqHEryIiCTaypVwyimw1VZw221xR1N91EQvIiKJdu21MGMGjBoFG28cdzTVRzV4ERFJrPHj4Y47QtP8gQfGHU31UoIXEZFEWrEiNM03bRqSfKFRE72IiCTSlVfCp5/C66/DhhvGHU31Uw1eREQSZ+xYuOceOPts2G+/uKOJhxK8iIgkyvLlcOqp0Lw53Hpr3NHER030IiKSKJddBp9/DqNHQ8OGcUcTH9XgRUQkMd56C/r2DXO977133NHESwleREQS4aef4LTTYIcd4Oab444mfmqiFxGRRPjHP+DLL+GddwpjrPmKqAYvIiJ579VXwzSwF18Mu+8edzS5QQleRETy2g8/wOmnw447FsY0sOlSE72IiOS1c86B+fPDsLQbbBB3NLlDCV5ERPLW9+8fyKBBcOON0Llz3NHkFjXRi4hIXvr6a/jqqcv405/Cve/yW0rwIiKSd4qL4eSTwYtr8sQTUEvt0b+jBC8iInnnnnvCoDbbHH0nLVrEHU1uUoIXEZG88vHHcPnlcNhhsOkeL8QdTs5SghcRkbzx889w4omw8cbwn/+AWdwR5S5dtRARkbxx9dUwdSqMGAGbbx53NLlNNXgREckLo0fDHXfAGWfAoYfGHU3uU4IXEZGct3QpnHQStGwJd94ZdzT5QU30IiKS8849F/73Pxg3Dho0iDua/KAavIiI5LTBg+HJJ8P19912izua/KEELyIiOeuLL8I19z/9Ca68Mu5o8osSvIiI5KTVq+H448OtcE89pdHqKktvl4iI5KTrroMJE0ITfbNmcUeTf1SDFxGRnPPWW/Dvf4d53o85Ju5o8pMSvIiI5JRFi8Joda1awb33xh1N/lITvYiI5Ax3OO20kORfekm3xK0PJXgREckZ998Pw4eH2eI6dow7mvymJnoREckJU6fCJZeEYWjPPz/uaPKfEryIiMSuqAiOPRY22QQefVSzxGVC1hK8mQ0wswVmNq2cbbqZ2RQzm25mo1OWH2Rms8xstpldlq0YRUQkN/TpA598Ekas22yzuKNJhmzW4AcCB61rpZltDDwAHObubYCjouU1gfuBg4GdgOPMbKcsxikiIjEaOhT694d//hP22y/uaJIjawne3ccA35ezyfHAMHf/Ktp+QbR8N2C2u89x91XAYODwbMUpIiLxmTMH/v73MMb8DTfEHU2ymLtn7+BmzYAR7t62jHX3ALWBNkAj4F53f9zMegAHufvfo+16Al3c/dx1nKM30Bug3qb1OrW5uU3G4p+5aCatG7fO2PHilqTyJKksoPLksiSVBXKrPMWrazPrtkf4eVFTWl95AnUbf1vpY+RSedZXVcoy+YzJk929c5kr3T1rD6AZMG0d6/oCE4AGQGPgM6AVoan+4ZTtegL3pXO+Tp06eSZ16pfZ48UtSeVJUlncVZ5clqSyuOdWec491x3cn3++6sfIpfKsr6qUBZjk68iJcd4HPw9Y5O7LgeVmNgboEC3fJmW7psA3McQnIiJZMmQI9O0LF10Eh+sibFbEeZvcC8BeZlbLzOoDXYCZwERgBzNrbmZ1gGOBF2OMU0REMmj27DDGfJcuYbx5yY6s1eDNbBDQDWhsZvOAawnX3HH3h9x9ppmNAqYCxYRm+WnRvucCrwA1gQHuPj1bcYqISPVZuRKOOipM/fr001CnTtwRJVfWEry7H5fGNrcDt5exfCQwMhtxiYhIfPr0gSlTwnC0220XdzTJppHsRESkWgwaBA89BP/4B3TvHnc0yacELyIiWTdrFvTuDbvvDjfdFHc0hUEJXkREsmrFCjj6aKhbN1x3r1077ogKg6aLFRGRrDr//DBT3MiR0LRp3NEUDtXgRUQkax5/HB5+GC6/HA4+OO5oCosSvIiIZMWUKXDGGdCtm8aZj4MSvIiIZNwPP8Df/gabbhquu9fSBeFqp7dcREQyqrgYTjwR5s2DMWNg883jjqgwKcGLiEhG3Xhj6FD3wAPQtWvc0RQuNdGLiEjGjBwJ118PJ58MZ54ZdzSFTQleREQyYs4cOOEE6NABHnwQzOKOqLApwYuIyHorKgqd6szg2Wdhgw3ijkh0DV5ERNaLO5x1VhjM5qWXYPvt445IQDV4ERFZTw89FAa0ufZaDWaTS5TgRUSkysaPhwsugEMOgauvjjsaSaUELyIiVfLtt9CjRxhf/sknoYYySk7RNXgREam0lStDp7qlS2HcONhkk7gjktKU4EVEpFJKOtVNmABDh0L79nFHJGVRg4qIiFTKfffBwIFwzTVw5JFxRyProgQvIiJpe+MNuOgiOPzw0GtecpcSvIiIpGXOHDj6aNhxR3jiCXWqy3X6eEREpELLloVauzu88AI0ahR3RFIRdbITEZFyFRfDSSfBjBkwahS0aBF3RJIOJXgRESnXv/4Fzz0Hd98Nf/5z3NFIutRELyIi6/T886Ez3cknhxHrJH9UKsGbWQMzq5mtYEREJHd8/DH07Am77RbGm9f0r/ml3ARvZjXM7Hgze8nMFgCfAN+a2XQzu93MdqieMEVEpDp99x107w4bbgjDhkG9enFHJJVVUQ3+LaAFcDmwpbtv4+6bA3sBE4BbzOzELMcoIiLVaMUKOOIIWLQIhg+HJk3ijkiqoqJOdvu7++rSC939e+BZ4Fkzq52VyEREpNq5w2mnhWFohw2DXXaJOyKpqnJr8CXJ3cy6mtkvdz2aWSMz65K6jYiI5L/rr4fBg+GWW+Cvf407Glkf6XayexBYlvJ6ebRMREQS4qmnQoI/9VS49NK4o5H1lW6CN3f3khfuXozuoRcRSYxx40LT/N57q8d8UqSb4OeY2flmVjt6XADMyWZgIiJSPebODZ3qmjYN193r1Ik7IsmEdBP8mcDuwP+AeUAXoHe2ghIRkerx44/hdrjVq2HECNh007gjkkxJq5nd3RcAx2Y5FhERqUZr1sAxx8CsWWGM+R13jDsiyaS0avBm1srM3jCzadHr9mZ2VXZDExGRbHGHPn1CYn/gAdhvv7gjkkxLt4n+P4TBblYDuPtUVKMXEclbd90FffvCxRdDr15xRyPZkG6Cr+/u75datibTwYiISPY98wxccgn06AG33RZ3NJIt6Sb4RWbWAnAAM+sBfJu1qEREJCt++qwjPXvCHnvAE09ADc0pmljp3st+DtAf2NHM/gd8AWgMehGRPPLJJ/D5A3fSsjm88IImkEm6dHvRzwH2N7MGQA13/ym7YYmISCbNnw8HHwxWcw0vv6zb4QpBur3oLzCzDYEi4G4z+8DMDshuaCIikgnLloV73RcsgJbnXkjz5nFHJNUh3asvp7n7j8ABwObAqcAtWYtKREQyouRe9w8/hKefhgbNZsYdklSTtMeij34eAjzq7h+lLBMRkRzkDuecAyNHwv33h1q8FI50E/xkM3uVkOBfiaaOLc5eWCIisr7+/W/o3x8uuwzOPDPuaKS6pduL/nSgIzDH3YvMbFNCM72IiOSgxx6DK6+E44+Hm26KOxqJQ7k1eDNrBmF6WHf/wN2XRK8Xu/tUC5pmP0wREUnXiBFw+ulh+NkBA3Sve6GqqAZ/u5nVAF4AJgMLgXpAS2BfYD/gWsIMcyIiErN334WjjoKOHeG556Bu3bgjkriUm+Dd/Sgz2wk4ATgN2Ipwq9xMYCRwk7uvLGtfMxsAdAcWuHvbMtZ3I/zj8EW0aJi73xCtmwv8BKwF1rh758oWTESk0EybFjrSbbstvPwyNGoUd0QSpwqvwbv7DODKKhx7INAXeLycbca6+7r6de7r7ouqcF4RkYIzdy4ceCDUrw+vvAKbbRZ3RBK3dDvZVZq7jym5hi8iItmzcGFI7kVFMGYMNGsWd0SSC8zds3fwkOBHlNNE/yzh+v03wCXuPj1a9wXwA2Fym37u3r+cc/QGegPU27RepzY3t8lY/DMXzaR149YZO17cklSeJJUFVJ5clutlWbuyPp/e9SArvmlBqwvPoWHLj8rdPtfLU1lJKk9VyjL5jMmT13kZ292z9gCaAdPWsW5DoGH0/BDgs5R1W0c/Nwc+AvZO53ydOnXyTOrUL7PHi1uSypOksrirPLksl8uycqX7fvu516zpPnx4evvkcnmqIknlqUpZgEm+jpxYbhO9me1S3np3/yD9/zN+t++PKc9HmtkDZtbY3Re5+zfR8gVm9hywGzCmqucSEUmatWuhZ0944w0YOFCj1MnvVXQN/s7oZz2gM6E2bUB74D1gz6qe2My2BL5zdzez3Qj35C9OnbEuen4AcENVzyMikjTucN55MGQI3H47nHxy3BFJLqroNrl9AcxsMNDb3T+OXrcFLilvXzMbBHQDGpvZPML98rWj4z4E9ADOMrM1wArg2CjZbwE8Z2Yl8T3l7qOqXEIRkQRxD0PPPvggXHopXFLuN7EUsnR70e9YktwB3H2amXUsbwd3P66C9X0Jt9GVXj4H6JBmXCIiBeWmm+C22+Css+AWzekp5Ug3wc80s4eBJwk9208kDHYjIiLV5J574Oqrw7X3vn3BNKenlCPdBH8qcBZwQfR6DPBgViISEZHfeeQR6NMH/vY3jS8v6Ukrwbv7SjO7H3idUIOf5e6rsxqZiIgAMHgw9OoFBx0ETz0FtbI2RJkkSVq/JtGgNI8Bcwm96Lcxs5PdXbeuiYhk0fDhoUl+r73g2Wc1eYykL93/A+8EDnD3WQBm1goYBHTKVmAiIoXujTfCzHA77xwSff36cUck+STdqzi1S5I7gLt/SnTLm4iIZN64cXDYYdCqFYwaBRtuGHdEkm/SrcFPMrNHgCei1ycQ5ocXEZEMmzwZDjkEmjSB116DP/wh7ogkH6Wb4M8CzgHOJ1yDHwM8kK2gREQK1QcfwP77wyabwOuvwxZbxB2R5Kt0e9H/DNwVPUREJAumTAnJfcMN4a23YNtt445I8llFk8084+5Hm9nHhNvjfsPd22ctMhGRAjJ1akjuDRvC229rTndZfxXV4EsGttE8RSIiWTJtGuy3H2ywQai5N28ed0SSBBVNNvOtmdUEHnH3/aspJhGRgjFjBvzf/0GdOiG5t2gRd0SSFBXeJufua4EiM9uoGuIRESkYn3wSknutWiG5t2wZd0SSJOn2ol8JfGxmrwHLSxa6+/lZiUpEJOFmzYJ99w3P33wz3O8ukknpJviXooeIiKynzz4Lyb24OHSo23HHuCOSJEr3NrnHzGwDYNvUEe1ERKRySprl16wJzfKtW8cdkSRVWkPVmtlfgCnAqOh1RzN7MYtxiYgkzrRpsM8+oeb+5pvQpk3cEUmSpTsW/XXAbsASAHefAuhGDhGRNH34IXTrFjrUjR4NbdvGHZEkXboJfo27Ly217HcD34iIyO9NnBia5evXD8n9j3+MOyIpBOkm+GlmdjxQ08x2MLP7gHFZjEtEJBHGjft1bPkxY3QrnFSfdBP8eUAb4GfgKWApcGGWYhIRSYTRo+GAA8KEMaNHa/hZqV4VjUVfDzgTaAl8DPzJ3ddUR2AiIvns9dfDfO7NmsEbb8BWW8UdkRSaimrwjwGdCcn9YOCOrEckIpLnRo6E7t1Dc/zbbyu5Szwqug9+J3dvB2BmjwDvZz8kEZH8NWwYHHts6CX/2muw6aZxRySFqqIa/OqSJ2qaFxEp36OPwlFHQefOoVleyV3iVFENvoOZ/Rg9N2CD6LUB7u4bZjU6EZE8cffdcNFFoVPdsGHQoEHcEUmhq2i62JrVFYiISD5yh2uugX/9C3r0gCefhLp1445KJP3JZkREpJTiYrjgAujbF04/Hfr1g5qqFkmOSPc+eBERSbF6NZx8ckjuF18M//mPkrvkFtXgRUQqaeVKOPpoGD4cbroJLr8czOKOSuS3lOBFRCrhp5/CADajR8P998PZZ8cdkUjZlOBFRNL03Xdw6KEwZUroTHf88XFHJLJuSvAiImmYPRsOPBDmz4cXXgiJXiSXKcGLiFRg4sSQ0N3hrbdgt93ijkikYupFLyJSjqXTdqdbN2jYEN59V8ld8ocSvIjIOgwcCLPvv4s//jHM696qVdwRiaRPCV5EpBR3uPlmOPVUaPTHSYweDVtuGXdUIpWja/AiIinWroXzz4cHHoATToAZu19Io0YT4g5LpNJUgxcRiaxYEQaweeAB+Mc/4PHHoUYtTaQp+Uk1eBERYMECOPxweO+9MDPchRfGHZHI+lGCF5GCN3NmuA1u/nx49ln461/jjkhk/SnBi0hBe+MNOPJIqFcvDD+7665xRySSGboGLyIFa8AAOOggaNo0NM0ruUuSKMGLSMEpLoYrrghzuO+7bxjAZrvt4o5KJLPURC8iBWXFCjjlFHjmGejVK8wIV7t23FGJZJ4SvIgUjIULQ0/58ePhttvgkks0j7sklxK8iBSEqVNDcp8/H4YODR3rRJJM1+BFJPGefx523x1WrYIxY5TcpTBkLcGb2QAzW2Bm09axvpuZLTWzKdHjmpR1B5nZLDObbWaXZStGEUk2d/jXv8J97W3ahGlf1VNeCkU2m+gHAn2Bx8vZZqy7d09dYGY1gfuBPwPzgIlm9qK7z8hWoCKSPEVFYbKYZ54JY8r/5z+wwQZxRyVSfbJWg3f3McD3Vdh1N2C2u89x91XAYODwjAYnIon29dew554wZAjceis88YSSuxQec/fsHdysGTDC3duWsa4b8Cyhlv4NcIm7TzezHsBB7v73aLueQBd3P3cd5+gN9Aaot2m9Tm1ubpOx+Gcumknrxq0zdry4Jak8SSoLqDyZtOzz9nz+0G0Ur6rH9n+/ko3avbtex9Nnk9uSVJ6qlGXyGZMnu3vnMle6e9YeQDNg2jrWbQg0jJ4fAnwWPT8KeDhlu57Afemcr1OnTp5Jnfpl9nhxS1J5klQWd5UnUx591L1OHfcWLdynT8/MMfXZ5LYklacqZQEm+TpyYmy96N39R3dfFj0fCdQ2s8aEGv02KZs2JdTwRUTKtGoVnHtuuOa+117w/vuw005xRyUSr9gSvJltaRaGmDCz3aJYFgMTgR3MrLmZ1QGOBV6MK04RyW3ffBOGm73/frjoInj5ZfjDH+KOSiR+WetFb2aDgG5AYzObB1wL1AZw94eAHsBZZrYGWAEcGzU3rDGzc4FXgJrAAHefnq04RSR/jR0LRx8NP/4IgwfDMcfEHZFI7shagnf34ypY35dwG11Z60YCI7MRl4jkP3fo2zfU2Js3h9deg7a/68orUtg0kp2I5JWiIujZE84/Hw45JAxeo+Qu8ntK8CKSNz7/HP70J3jqKbjxRnjuOdhoo7ijEslNmmxGRPLCSy/BiSeG2d9GjoSDDoo7IpHcphq8iOS01avhssuge3do1gwmT1ZyF0mHavAikrPmzYNjj4V334UzzoB77oF69eKOSiQ/KMGLSE56+eXQme7nn8M19+PKvS9HREpTE72I5JQ1a+CKK0IP+SZNYNIkJXeRqlANXkRyxv/+F5L52LHQu3doktcscCJVowQvIjnhlVdCL/kVK+DJJ8Mc7iJSdWqiF5FYrVoFl14KBx8MW24ZmuSV3EXWn2rwIhKbTz+F448Pt76deSbceSfUrx93VCLJoAQvItXOHR59FM47L9z29txzcMQRcUclkixqoheRavXDD2EGuNNPh65dYepUJXeRbFCCF5FqM2YMdOgAzz8Pt94aZoFr0iTuqESSSQleRLJu9Wq4+mrYd1+oWxfGjQsd62roG0gka3QNXkSyatYsOOkkeP99OPVU+H//Dxo2jDsqkeTT/88ikhXFxSGZd+wIs2fD00/DgAFK7iLVRTV4Ecm4Vd9vwZ//DG++GYacffhh2GqruKMSKSxK8CKSMe7w2GMw/fqnqV8L/vOf0FveLO7IRAqPEryIZMSCBWH8+BdegIY7zGLqK51o3jzuqEQKl67Bi8h6GzYM2rSBUaPCaHStLjpTyV0kZkrwIlJlCxfCscfCkUfCttvCBx/ARReB1fC4QxMpeErwIlJp7vDUU9C6dRhm9sYbYcIE2GmnuCMTkRK6Bi8ilfK//4WJYUaMgC5dwq1vSuwiuUc1eBFJi3voFb/TTvDGG3DXXfDuu0ruIrlKNXgRqdCcOdCrV7ivfd99Q6Jv0SLuqESkPKrBi8g6rVkTaurt2sHEidCvH7z+upK7SD5QDV5EyvT++3DGGTBlChx6KDz4IGyzTdxRiUi6VIMXkd9YuhTOOSfM1b5gAQwdCsOHK7mL5BsleBEBQie6p5+GHXeEhx6C886DmTPDPe4aalYk/6iJXkSYMwfOPhteeQU6dQo19s6d445KRNaHavAiBeznn+Hmm8Mws+PGwb33wnvvKbmLJIFq8CIF6qWX4MILw1ztRx4ZknuTJnFHJSKZohq8SIGZPRu6dw+PGjXCBDFDhyq5iySNErxIgVi+HK64IjTHjx4Nt98OH38MBx4Yd2Qikg1qohdJuJLe8ZdcEsaR79kTbr0Vttoq7shEJJtUgxdJsI8+CkPLHnccbL45vPMOPP64krtIIVCCF0mgb76B00+HnXcOzfAPPRSGmt1jj7gjE5HqoiZ6kQRZvhzuuANuuw1Wr4Y+feCqq2CTTeKOTESqmxK8SAKsXRua3q+6KtTee/SAW27RpDAihUxN9CJ57s03w8A0p50Wxot/5x0YMkTJXaTQKcGL5Klp0+Avf4H99oMffoBBg2D8eF1nF5FACV4kz8yZE251a98exo4NTfGffALHHqtJYUTkV7oGL5In5s+Hf/0L+veHmjXhH/+Af/4T/vCHuCMTkVykBC+S45YsCb3i770XVq2Cv/8drr4att467shEJJcpwYvkqKIiuO++0AS/ZEkYrOaGG6Bly7gjE5F8oAQvkmOKisLANLfeCgsWwKGHwk03QYcOcUcmIvlEnexEcsTy5XDnndC8OVx88a+d6EaMUHIXkcrLWoI3swFmtsDMplWw3a5mttbMeqQsm2tmH5vZFDOblK0YRXLB8uXw3WsnsP32YUKYksT+2muw555xRyci+SqbNfiBwEHlbWBmNYFbgVfKWL2vu3d0985ZiE0kdiU19u23h3lD+9ChQxikRoldRDIha9fg3X2MmTWrYLPzgGeBXbMVh0iuWbIEHngA7rkHFi6EAw6ALzuexqu3Dog7NBFJEHP37B08JPgR7t62jHVNgKeA/wMeibYbGq37AvgBcKCfu/cv5xy9gd4A9Tat16nNzW0yFv/MRTNp3bh1xo4XtySVJx/LsmpJYxa8cTwLx/yN4pUN2bDtO2x1yAAatpial+UpT5LKk6SygMqTy6pSlslnTJ68zpZud8/aA2gGTFvHuiFA1+j5QKBHyrqto5+bAx8Be6dzvk6dOnkmdeqX2ePFLUnlyaeyfPqpe69e7nXquNeo4X788e5Tpvx2m3wqTzqSVJ4klcVd5cllVSkLMMnXkRPjvE2uMzDYwtiajYFDzGyNuz/v7t8AuPsCM3sO2A0YE1+oIpU3eXK41W3oUKhTJ8zPfskl4Zq7iEi2xZbg3b15yXMzG0hoon/ezBoANdz9p+j5AcANMYUpUinFxTByJNx9d5jlbcMN4bLL4IILYIst4o5ORApJ1hK8mQ0CugGNzWwecC1QG8DdHypn1y2A56KafS3gKXcfla04RTJh+XJ47LEwnOynn0KTJmEEujPPhI02ijs6ESlE2exFf1wltj0l5fkcQMN6SF74+mvo2zdMALNkCey6Kzz1FPToAbVrxx2diBQyDVUrUgXvvx+a4YcMAXf429+gTx/40580ZauI5AYleJE0FRXB4MHhHvbJk8P19QsvhHPPhWbN4o5OROS3lOBFKjBrVpj8ZeDA0Azfpk1olj/pJGjUKO7oRETKpgQvUobVq+HFF+HBB+GNN8L19COPhLPPDsPIqhleRHKdErxIijlz4NFHYcAA+OYb2HZbuPlmOO003eYmIvlFCV4KXlERDBsGjzwCb78NNWrAgQdCv35w8MFQs2bcEYqIVJ4SvBQkd5g4MdTUBw2CH3+EFi3gppvCtfWmTeOOUERk/SjBS0H59tuQ0B99FKZNgw02gKOOCk3we++ta+sikhxK8JJ4P/4YmuD/+98wfGxxMXTpEprgjz023O4mIpI0SvCSSKtWwcsvh6Q+fDisXBkmebnySjjhBPjjH+OOUEQku5TgJTHWroWxY8NgNEOGwPffQ+PGYRa3E06Arl3VBC8ihUMJXvLamjWh5/vQofDcc7BgAdSvD0ccEZL6n/+sMeFFpDApwUveWb0alk7vyt//Ds8/D4sXQ4MGcOihYZKXgw+Ghg3jjlJEJF5K8JIXfvoJXnstjC734ovwww99+a4R/OUvIakfeGCouYuISKAELznrq69CB7nhw+Gtt0LHuY03hsMOg3cb9mHanXdTr17cUYqI5CYleMkZxcVh8JmSpD51aljeqhWcd16ore+xB9SqBZ37j1VyFxEphxK8xGrevND0/uqr4efixWFo2D32gNtvD0ldt7SJiFSeErxUq6IiGDMmJPRXX4Xp08PyLbcMneQOOCB0kvvDH+KNU0Qk3ynBS1b9/DO8/364le3tt+Hdd8OyunXD0LCnnBKSert2ukddRCSTlOAlo1auhAkTYPTokNAnTAjLzEISP+eckND32ku93kVEskkJXtbLwoUhiY8fH2rn770Xaug1akDHjnDWWbDPPiGhq9ldRKT6KMFL2lavho8/Dsl8/PiQ2D//PKyrVQt23jn0dt9nH9hzz3BLm4iIxEMJXsq0Zg3MnAkffAAffgiTJ4fHihVh/VZbwZ/+BGecEX526hSmXhURkdygBC8UFYXe7B9+GBL6Bx+EmvrKlWF9/frQoQP07h2SedeusO226hQnIpLLlOALSPHqOkyZEpJ56mPOHHAP22y0EeyyS+gMt8suodm9Vatwb7qIiOQPJfiEKS4Og8d89hnMnh0en30GM2bAZ7PHsnOUyGvVCol7l12gZ09o2zY8b9ZMNXMRkSRQgs8z7mG0t6++Co+vv4a5c39N5HPmhF7sJerWhe23D7eoLWs1gLt79qJNm5Dc69SJrRgiIpJlSvA5wh2WLYP58+G77377mDfv12T+1Ve/dnQrUa8etGwZhnTt3j08b9kSdtgBmjQJt6wBdO7fj2OO6VX9hRMRkWqnBJ8ha9eGxFtU9NvHjz/CDz/AkiVl//z++18TeenEDaG5fMstQ6e29u3DcK7bbvvrY5ttYLPN1KwuIiK/ZV7SuyoBOnfu7JMmTcrIsV5scSEb/+9d6tduAB5q2O7gpDx3KF4La4t/7aRWESNc//7lURvq1A7N5bXrhJ8lr+vUgdq1M5e8J387mU5bdcrMwWKWpLKAypPLklQWUHli17Ej3HNPmas69+/MpN6Vy2FmNtndO5e1TjX4dahdG6zmGurWASwk2V8e/Pq8Rk2oWaPsnzVqQK2aIYmXJPSaNcP+IiIi2aQEvw4Hf3JP9N/U23GHkjFnJKg8SSoLqDy5LEllAZWnkNSIOwARERHJPCV4ERGRBFKCFxERSSAleBERkQRSghcREUkgJXgREZEEUoIXERFJICV4ERGRBFKCFxERSSAleBERkQRSghcREUkgJXgREZEEUoIXERFJoETNB29mC4EvM3jIxsCiDB4vbkkqT5LKAipPLktSWUDlyWVVKct27r5ZWSsSleAzzcwmuXvnuOPIlCSVJ0llAZUnlyWpLKDy5LJMl0VN9CIiIgmkBC8iIpJASvDl6x93ABmWpPIkqSyg8uSyJJUFVJ5cltGy6Bq8iIhIAqkGLyIikkAFn+DN7CAzm2Vms83ssjLWm5n9v2j9VDPbJY4405VGeU6IyjHVzMaZWYc44kxXReVJ2W5XM1trZj2qM77KSKcsZtbNzKaY2XQzG13dMVZGGr9rG5nZcDP7KCrPqXHEmQ4zG2BmC8xs2jrW59v3QEXlybfvgXLLk7JdPnwPVFiWjH0PuHvBPoCawOfA9kAd4CNgp1LbHAK8DBjQFXgv7rjXszy7A5tEzw/O9/KkbPcmMBLoEXfc6/HZbAzMALaNXm8ed9zrWZ4rgFuj55sB3wN14o59HeXZG9gFmLaO9XnzPZBmefLmeyCd8kTb5Pz3QJqfTca+Bwq9Br8bMNvd57j7KmAwcHipbQ4HHvdgArCxmW1V3YGmqcLyuPs4d/8hejkBaFrNMVZGOp8PwHnAs8CC6gyuktIpy/HAMHf/CsDd8708DjQyMwMaEhL8muoNMz3uPoYQ37rk0/dAheXJs++BdD4fyI/vgXTKkrHvgUJP8E2Ar1Nez4uWVXabXFHZWE8n1EpyVYXlMbMmwF+Bh6oxrqpI57NpBWxiZm+b2WQzO6naoqu8dMrTF2gNfAN8DFzg7sXVE17G5dP3QGXl+vdAhfLoeyAdGfseqJXBoPKRlbGs9G0F6WyTK9KO1cz2Jfxh75nViNZPOuW5B/inu68NFcWclU5ZagGdgP2ADYDxZjbB3T/NdnBVkE55DgSmAP8HtABeM7Ox7v5jlmPLhnz6HkhbnnwPpOMe8uN7IB0Z+x4o9AQ/D9gm5XVTQm2jstvkirRiNbP2wMPAwe6+uJpiq4p0ytMZGBz9UTcGDjGzNe7+fLVEmL50f9cWuftyYLmZjQE6ALmY4NMpz6nALR4uJM42sy+AHYH3qyfEjMqn74G05NH3QDry5XsgHRn7Hij0JvqJwA5m1tzM6gDHAi+W2uZF4KSoF21XYKm7f1vdgaapwvKY2bbAMKBnjtYMU1VYHndv7u7N3L0ZMBQ4O0f/qNP5XXsB2MvMaplZfaALMLOa40xXOuX5ilALwcy2AP4IzKnWKDMnn74HKpRn3wMVyqPvgXRk7HugoGvw7r7GzM4FXiH0wBzg7tPN7Mxo/UOEHpmHALOBIkKtJCelWZ5rgE2BB6L/dtd4jk7UkGZ58kI6ZXH3mWY2CpgKFAMPu3u5twXFJc3P5kZgoJl9TGji/qe75+SsX2Y2COgGNDazecC1QG3Iv+8BSKs8efM9AGmVJ29UVJZMfg9oJDsREZEEKvQmehERkURSghcREUkgJXgREZEEUoIXERFJICV4ERGRBFKCF8kjZvZXM3Mz2zGDx+xmZiOi54dZNDOcmR1hZjtV4Xhvm1mlbrkys3vMbO/KnqvUMZZFPzeLbjMSKWhK8CL55TjgHcLAMhnn7i+6+y3RyyOASif4yjKzPwBdo0k4Sq+rWdnjuftC4Fsz2yMT8YnkKyV4kTxhZg2BPQhjhx+bsrybmY02s2fM7FMzu8XCfN/vm9nHZtYi2m6gmT1kZmOj7bqXcY5TzKyvme0OHAbcHs1L3SK1Zm5mjc1sbvR8AzMbbGFu8acJ42eXHO8AMxtvZh+Y2ZCoDKX1AEal7DPXzK4xs3eAo8ysl5lNtDCv/LPR6F5Eo+iNj9bdWOqYzwMnVPpNFkkQJXiR/HEEMCoaWvR7M9slZV0H4AKgHdATaOXuuxHGGj8vZbtmwD7AocBDZlavrBO5+zjC8Kz/cPeO7v55OXGdBRS5e3vgJsJEGZhZY+AqYH933wWYBFxUxv57AJNLLVvp7nu6+2DC1Jm7unsHwpCdp0fb3As86O67AvNL7T8J2KucmEUSTwleJH8cR5h3nejncSnrJrr7t+7+M/A58Gq0/GNCUi/xjLsXu/tnhHHhM3Etf2/gSQB3n0oYYhOgK6GJ/10zmwKcDGxXxv5bAQtLLXs65XnbqNXhY0KtvE20fA9gUPT8iVL7LwC2rnRJRBKkoMeiF8kXZrYpYdrVtmbmhPHf3cwujTb5OWXz4pTXxfz277z02NSVGat6Db9WCkrX/Ms6jgGvuftxZaxLtaKM4y1PeT4QOMLdPzKzUwjjeJd33pL4VlRwXpFEUw1eJD/0AB539+2iWbO2Ab6g8vN4H2VmNaLr8tsDs8rZ9iegUcrruUTN71E8JcYQXe82s7ZA+2j5BGAPM2sZratvZq3KOM9MoGU5cTQidJqrzW+vq7/Lr30RSl9vbwXk5EQ9ItVFCV4kPxwHPFdq2bPA8ZU8zixgNPAycKa7ryxn28HAP8zsw+gfgjuAs8xsHGHO7RIPAg3NbCpwKdF871Fv9lOAQdG6CZR9SeAlflsrL+1q4D3gNeCTlOUXAOeY2URgo1L77BsdV6RgaTY5kQJhZgOBEe4+NO5YSot6zHd39yUZOt4Y4HB3/yETxxPJR6rBi0guuBjYNhMHMrPNgLuU3KXQqQYvIiKSQKrBi4iIJJASvIiISAIpwYuIiCSQEryIiEgCKcGLiIgkkBK8iIhIAv1/N4PPN2oxOycAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "#\n", "def periodLima(L, thetaZero):\n", " #\n", " T = -2*np.pi*np.sqrt(L/g)*np.log(np.cos(thetaZero/2))/(1 - np.cos(thetaZero/2))\n", " #\n", " return T\n", "#\n", "g = 9.81 # m/s**2\n", "L = 0.5 # m\n", "periodSHM = 2*np.pi*np.sqrt(L/g)\n", "thetaBot = 0.0001\n", "thetaTop = np.pi/2\n", "#thetaTop = np.pi\n", "nTheta = 240\n", "thetaArr = np.linspace(thetaBot, thetaTop, nTheta)\n", "periodArr = periodLima(L, thetaArr)\n", "#\n", "plt.figure(figsize = (8, 6))\n", "plt.title('Period as function of amplitude', fontsize = 14)\n", "plt.xlabel('Amplitude (rad)')\n", "plt.ylabel('Period (sec)')\n", "plt.plot(thetaArr, periodArr, label = \"LA formula\", color = 'b')\n", "plt.plot(thetaArr, periodSHM*np.ones(nTheta), label = \"SHM\", color = 'r')\n", "plt.grid(color = 'g')\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Large angle pendulum experimental study\n", "\n", "We will look at data taken using a pendulum of length $l = 0.5$ m (see below for how the \"measurements\" were really made!). The data are a set of measurements of the period for a range of angular amplitudes. We will investigate whether a significant deviation from the \"simple\" form for $T$ can be identified. The data are available as three files (which you can look at using MS excel):\n", "\n", "* _thetaTable.csv_ contains the values of the amplitude of the pendulum's oscillations in radians. This is just a single column of data.\n", "\n", "* _periodTableSmall.csv_ contains a table formed of $10$ measurements of the period for each of the amplitudes in _thetaTable.csv_. The rows in _periodTableSmall.csv_ correspond to the rows in _thetaTable.csv_. E.g. the third row (counting from zero!) contains the measured periods for the amplitude $0.394$ radians.\n", "\n", "* _periodTableLarge.csv_ is the same as _periodTableSmall.csv_, except that there are $1000$ measurements for each amplitude of the pendulum." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "!\"C:\\Program Files\\Microsoft Office\\root\\Office16\\EXCEL.EXE\" thetaTable.csv" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "!\"C:\\Program Files\\Microsoft Office\\root\\Office16\\EXCEL.EXE\" periodTableSmall.csv" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "!\"C:\\Program Files\\Microsoft Office\\root\\Office16\\EXCEL.EXE\" periodTableLarge.csv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Low precision experiment\n", "We will read in _thetaTable.csv_ and _periodTableSmall.csv_ in the following cell, then check that they have the expected numbers of rows and columns (using the `array.shape` syntax, which gives us access to the length of 1D arrays, the number of rows and columns for 2D arrays, etc.)." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " \n", "Shape of thetaArr = (16,)\n", "Shape of periodArrSm = (16, 10)\n", "nTheta = 16 and nMeas = 10\n", " \n", "thetaArr =\n", " [0.1 0.198 0.296 0.394 0.492 0.59 0.688 0.786 0.884 0.982 1.081 1.179\n", " 1.277 1.375 1.473 1.571]\n", " \n", "periodArr =\n", " [[0.591 1.789 1.382 1.464 0.703 1.228 1.505 1.482 1.944 1.416]\n", " [2.099 0.75 2.135 0.863 1.935 2.659 0.743 1.278 1.377 1.27 ]\n", " [1.895 0.796 0.682 1.112 2.059 1.19 1.244 1.811 1.625 1.547]\n", " [0.974 1.159 1.413 1.623 1.171 1.006 1.996 0.705 2.115 1.739]\n", " [1.355 2.282 1.99 1.163 1.876 1.06 1.445 1.83 1.228 1.055]\n", " [1.272 1.759 1.541 1.655 1.322 1.422 1.285 1.597 1.594 2.151]\n", " [1.47 1.306 2.062 1.889 1.078 1.764 1.773 0.83 1.226 0.583]\n", " [2.197 1.705 0.886 1.138 1.592 0.42 1.19 0.595 2.127 1.243]\n", " [2.161 1.667 2.033 2.044 1.974 1.115 2.041 0.291 0.736 1.802]\n", " [1.412 1.146 1.2 1.249 0.391 2.078 1.098 1.193 1.279 1.459]\n", " [0.953 1.445 1.088 2.064 1.031 2.737 1.232 2.29 1.313 1.082]\n", " [1.144 1.513 0.802 1.711 0.77 1.297 1.363 0.853 1.238 1.552]\n", " [1.693 1.682 2.373 1.325 1.157 2.177 2.753 0.718 1.561 1.254]\n", " [2.11 1.619 2.506 1.877 1.536 1.814 0.941 1.271 2.136 1.844]\n", " [1.591 1.535 2.081 1.998 1.032 1.27 1.384 2.306 1.323 1.84 ]\n", " [1.88 2.343 0.893 1.439 1.508 2.369 1.636 1.474 0.874 2.362]]\n" ] } ], "source": [ "thetaArr = np.loadtxt(\"thetaTable.csv\", delimiter = ',')\n", "periodArrSm = np.loadtxt(\"periodTableSmall.csv\", delimiter = ',')\n", "np.set_printoptions(precision = 3)\n", "print(\" \")\n", "print(\"Shape of thetaArr =\",thetaArr.shape)\n", "print(\"Shape of periodArrSm =\",periodArrSm.shape)\n", "#\n", "nTheta = thetaArr.shape[0]\n", "nMeas = periodArrSm.shape[1]\n", "print(\"nTheta =\",nTheta,\"and nMeas =\",nMeas)\n", "print(\" \")\n", "print(\"thetaArr =\\n\",thetaArr)\n", "print(\" \")\n", "print(\"periodArr =\\n\",periodArrSm)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now calculate the mean value of the period for each angle:\n", "\n", "$$ \\overline {T} = \\frac{1}{N}\\sum_{n = 1}^N T_n.$$\n", "\n", "The *standard error* of the mean is given by:\n", "\n", "$$ \\delta \\overline {T} = \\frac{1}{\\sqrt {N}} \\sqrt{ \\sum_{n = 1}^N \\frac{(T_n - \\overline {T})^2}{N - 1}}.$$\n", "\n", "In this formula, the factor\n", "\n", "$$s = \\sqrt{ \\sum_{n = 1}^N \\frac{(T_n - \\overline {T})^2}{N - 1}}$$\n", "\n", "describes the spread of $T$ values we obtain from our experiment. It is the sample standard deviation of the $T$ distribution. (While the precision with which we know the average of this distribution increases as we make more measurements by the factor $1/\\sqrt N$, \n", "the value of $s$ is a property of the apparatus and experimental procedure we are using; it doesn't change - within errors - with the number of measurements we make. We can only decrease $s$ by using better apparatus or improving our experimental procedure.)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " \n", "Number\t Amplitude (rad)\t Mean period (s)\n", "0\t 0.10\t\t\t 1.35 +- 0.13\n", "1\t 0.20\t\t\t 1.51 +- 0.21\n", "2\t 0.30\t\t\t 1.40 +- 0.15\n", "3\t 0.39\t\t\t 1.39 +- 0.15\n", "4\t 0.49\t\t\t 1.53 +- 0.14\n", "5\t 0.59\t\t\t 1.56 +- 0.08\n", "6\t 0.69\t\t\t 1.40 +- 0.15\n", "7\t 0.79\t\t\t 1.31 +- 0.19\n", "8\t 0.88\t\t\t 1.59 +- 0.20\n", "9\t 0.98\t\t\t 1.25 +- 0.13\n", "10\t 1.08\t\t\t 1.52 +- 0.20\n", "11\t 1.18\t\t\t 1.22 +- 0.10\n", "12\t 1.28\t\t\t 1.67 +- 0.19\n", "13\t 1.37\t\t\t 1.77 +- 0.14\n", "14\t 1.47\t\t\t 1.64 +- 0.13\n", "15\t 1.57\t\t\t 1.68 +- 0.18\n" ] } ], "source": [ "#\n", "# Calculate means for each amplitude\n", "meanPeriods = np.mean(periodArrSm, axis = 1)\n", "#\n", "# Calculate sample standard deviation using numpy function, hence get standard error of mean.\n", "errPeriods = np.std(periodArrSm, axis = 1, ddof = 1)/np.sqrt(nMeas)\n", "print(\" \")\n", "print(\"Number\\t Amplitude (rad)\\t Mean period (s)\")\n", "for n in range(0, nTheta):\n", " print(f\"{n:d}\\t {thetaArr[n]:.2f}\\t\\t\\t {meanPeriods[n]:.2f} +- {errPeriods[n]:.2f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now plot our measurements and compare them with the theoretical SHM and Lima curves. " ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " \n", "Chi2 for simple approximation = 23.01, chi2/NDF = 1.53\n", " \n", "Chi2 for Lima/Arun approximation = 19.15, chi2/NDF = 1.28\n", " \n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfEAAAGECAYAAAAm3RkPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA5IklEQVR4nO3deZgU1b3/8feXTRZZFHBFWVTAIMM24ALCuMZ9iQtBRdEYBJdorjcxRiMTlxgj3uvNzwQh0YCJcTeu0aAiokFEiDiAiAuCoggCKjuyfH9/nJ5hZuiZ6YGu7q6ez+t56unuquqqc6ab/nBqOcfcHREREYmfetkugIiIiOwYhbiIiEhMKcRFRERiSiEuIiISUwpxERGRmFKIi4iIxJRCXCTPmFmhmbmZdajFe4oS72kTYdFStiN1EKmLFOIiETGz8YkgcjPbZGYLzGy0mTXLdtnykZlNNrN7sl0OkUxqkO0CiOS5l4GhQEPgSODPQDNgZDYLJSL5QS1xkWhtdPcv3f0zd/878CBwBoAFPzezj81svZnNNrMLSt9oZh0SrfizzOwlM1tnZu+Z2XHld2BmJ5jZ+2a2wcxeBzpXWj7MzNZUmlft4fNU3lO6jpmdmNj/OjN7xsxamtnZZvahmX1rZn81sybV/ZFSqENrM3vIzBYn/lZzzezicsvHA4OAK8od/ehgZvXN7D4z+yTxvg8Tf3P99kle0BdZJLPWE1rlALcCPwKuAL4H3A6MNbOTK73nNuD3QA/gbeBhM9sVwMz2A54CXgJ6Av8P+F2kNahoF+Ba4HzgGKAQeBy4CDiL8B+WU4DLq9pAinVoDPwnsa1uwP8R/lbHJJZfDbwJ/AXYOzF9RviN+xw4FzgYuAH4JXAxInlAh9NFMsTM+gHnAa8kzov/F3C8u7+eWOWTxDpXAM+Xe+v/uvuziW38EriQEHZvEA7Lfwr8xMNACO+bWWfglgxUCcJvyBXuPj9Rvr8DPwX2dPfliXlPA0cBd1WxjRrr4O6fA3eWe884MzsaGAK84u7fmtl3wDp3/7LceluAm8q9XmhmvRPvu29HKy2SKxTiItE6IXFYugGhBf40cBWh5d0YeNHMyo9C1BBYWGkbJeWef5F43CPxeDAwzSuOZPRmeoqeko2lAZ6wFPiyNMDLzfteNduosQ5mVh/4BTAY2JdwBKARMLmmAprZCOBSoD3QhPA3XlTT+0TiQCEuEq0pwHBgE/CFu28CMLOOieWnElqh5W2q6rW7u5nBtlNhlkIZtiZZr2GyFXfgPZsrvXa2L79T/am7VOrw34TD9lcDs4E1wG/Y9p+Z5Bs2GwzcnXj/VGAV4UjHmSnsUyTnKcRForXO3T9KMv89YCPQ3t0n7cT23wPOMjMr15I9rNI6XwFNzayFu69KzOtZw3Z35D07KpU6DACedfe/QrgokHDx2zfl1vkOqJ/kfW+5e9mtZ2Z2QBrLLpJVurBNJAvcfTUwGhhtZpeY2YFm1tPMRpjZ8Fps6l6gA3C3mXUxs7OBEZXWeQtYC9ye2M9ZVHOh2U68Z0elUocPgGPMbICZdQXuATpWWmch0C9xVXqbxBXoHwC9E1fQH2RmvyJcxS6SFxTiItnzK6CYcKh3LuHq7LOAT1LdgLt/CvwAOAF4l3BR2S8qrbOScPX4cYRD0cMT+65uu7V+z45KpQ6EK/mnAy8QTlGsJdyuV95oQmv8PcKRhP2BscCjwN8JV/Z3oOoL7ERixypeSyIiIiJxoZa4iIhITCnERUREYkohLiIiElMKcRERkZhSiIuIiMRU7Dp7adS8kRd0Kch2MdLmw5UfctDuB2W7GGmTT/XJp7qA6pPL8qkuoPpEYebMmcvdvW3l+bEL8Ya7N2TGjBnZLkbaFI4rZMZw1ScX5VNdQPXJZflUF1B9omBmSfv71+F0ERGRmFKIi4iIxJRCXEREJKZid05cRERqZ9OmTSxevJgNGzZkZH+/6/k75s2bl5F9ZUIm69O4cWPatWtHw4Y1jRYcKMRFRPLc4sWLad68OR06dCAxHn2k/Cvn4LYHR76fTMlUfdydFStWsHjxYjp2rDxIX3I6nC4ikuc2bNhA69atMxLgsuPMjNatW9fqiIlCXESkDlCAx0NtPyeFuIiIRG7p0qWcd955dOrUiT59+nD44Yfzj3/8Iy3bLioqyqv+Q2pDIS4iIpFyd8444wwGDhzIggULmDlzJg8//DCLFy/OdtFqbcuWLdkuQgUKcRERidSkSZNo1KgRI0aMKJvXvn17rrrqKiCcs7/44ovp3r07vXr14tVXX612/vr16/nhD39IQUEBgwcPZv369Un3e/PNN9O3b18OOeQQhg8fjrsDoeV+zTXXcMQRR3DIIYcwffp0AIqLixk6dChHH300Bx10EH/6058AmP7v6Rx11FGcd955dO/evcpynX766TzwwAMAjB07lvPPPz/df8rt6Op0EZG65JprYNas9G6zZ0+4++4qF8+dO5fevXtXufwPf/gDALNnz+b999/n+OOP54MPPqhy/pgxY2jatCklJSWUlJRUue0rr7ySm266CYChQ4fy3HPPceqppwKwdu1apk6dypQpU7jkkkuYM2cOACUlJUybNo21a9fSq1cvTj75ZACmT5/OnDlz6NixI3fddVfSco0bN47+/fuXrTNt2rTU/4Y7SC1xERHJqCuuuIIePXrQt29fAN544w2GDh0KQNeuXWnfvj0ffPBBlfOnTJnCBRdcAEBBQQEFBckHxXr11Vc59NBD6d69O5MmTWLu3Llly4YMGQLAwIEDWbVqFd988w0QWtNNmjShTZs2HHXUUWWt9H79+pXd9lVVufbcc09uvvlmjjrqKO666y523333dP7ZklJLXESkLqmmxRyVbt268cQTT5S9/sMf/sDy5cspLCwEKDvMXVlV86Hmq7g3bNjA5ZdfzowZM9hvv/0oLi6ucOtW5feXvk4636FZs2YplWv27Nm0bt2aL774otrypYta4iIiEqmjjz6aDRs2MGbMmLJ569atK3s+cOBAHnzwQQA++OADPv30U7p06ZLS/Dlz5lBSUrLdPksDu02bNqxZs4bHH3+8wvJHHnkECK3qli1b0rJlSwCefvppNmzYwIoVK5g8eXLZ0YLyqirX9OnTeeGFF3jnnXcYPXo0n3zyyY79wWpBLXEREYmUmfHUU0/x05/+lN/97ne0bduWZs2acccddwBw+eWXM2LECLp3706DBg0YP348u+yyS5XzR44cycUXX0xBQQE9e/akX79+2+2zVatW/PjHP6Z79+506NBhuzDebbfdOOKII1i1ahX3339/2fx+/fpx8skn8+mnn/KrX/2KffbZZ7ttJysXwI9//GP+8pe/sM8++3DXXXdxySWXMGnSJMaOHQtQ4cK+dFGIi4hI5Pbee28efvjhpMsaN25cFoSpzG/SpEmV2yrv1ltv5dZbb0267KyzzuL222/fbn7nzp0ZN25chXn9+vdj2BnDaizXu+++W/b8tNNO47TTTgOiCe9SOpwuIiISU2qJi4hInTJ58uSk84uLizNajnRQS1xERCSmFOIiIiIxpRAXERGJKYW4iIhITCnERUQkcrfddhvdunUru7f7rbfeAuDSSy/lvffeS8s+dt1117RsJ050dbqIiETqzTff5LnnnuM///kPu+yyC8uXL+e7774D4M9//nOWS5demzdvpkGDzEWrWuIiIhKpJUuW0KZNG3bZZRcgdIVa2hNaUVERM2bMAEJL+rrrrqNPnz4ce+yxTJ8+naKiIjp16sQzzzwDwPjx4zn99NM54YQT6NKlC7/+9a+T7vPOO++kb9++FBQUMGrUqKTrjBw5ksLCQrp161ZhnQ4dOnDdddfRr18/+vXrx6IFiwAYNmwYI0aM4Mgjj6Rz584899xzZWU655xzOPXUUzn++ONZuXIlZ5xxBgUFBRx22GGUlJSwefNm+vbtW3Z72/XXX88NN9ywk39ZtcRFROqeoqLt5517Llx+OaxbByedtP3yYcPCtHw5nH12xWVV3Hdd6vjjj+fmm2+mc+fOHHvssQwePJhBgwZtt97atWspKirijjvu4Mwzz+TGG2/kpZde4r333uOiiy4q6wGtdFjQpk2b0rdvX04++eSywVQAJk6cyIcffsj06dNxd0477TSmTJnCwIEDK+zvtttuY/fdd2fLli0cc8wxlJSUlI2I1qJFC6ZPn84DDzzAHb+6gxMnngjAwoULee211/j444856qij+Oijj4BwtKGkpITdd9+dq666il69evHUU08xadIkLrzwQmbNmsX48eM5++yz+f3vf8+LL75YdkphZ6glLiIikdp1112ZOXMm48aNo23btgwePDhpt6WNGjXihBNOAKB79+4MGjSIhg0b0r17dxYuXFi23nHHHUfr1q1p0qQJP/jBD3jjjTcqbGfixIlMnDiRXr160bt3b95//30+/PDD7fb36KOP0rt3b3r16sXcuXMrnJsvHap0yJAhzJoxq2z+ueeeS7169TjooIPo1KkT77//flmZSoceLT9U6dFHH82KFSv49ttv6datG0OHDuXUU0/l/vvvp1GjRrX/Y1ailriISF1TXcu5adPql7dpU2PLO5n69etTVFREUVER3bt3Z8KECQwbNqzCOg0bNiwbCrRevXplh9/r1avH5s2by9arahjRUu7O9ddfz2WXXVZleT755BNGjx7N22+/zW677cawYcOqHKq0quflX9c0VGnperNnz6ZVq1YsXbq0yrLVRmQtcTPbz8xeNbN5ZjbXzK5Oso6Z2e/N7CMzKzGz3lGVR0REsmP+/PkVWsKzZs2iffv2O7y9l156iZUrV7J+/Xqeeuop+vfvX2H597//fe6//37WrFkDwOeff86yZcsqrLNq1SqaNWtGy5YtWbp0KS+88EKF5aVDlT7yyCP0KOxRNv+xxx5j69atfPzxxyxYsIAuXbpsV77yQ5VOnjyZNm3a0KJFC5588klWrFjBlClT+MlPfsI333yzw3+DUlG2xDcD17r7f8ysOTDTzF5y9/L3EpwIHJSYDgXGJB5FRCRPrFmzhquuuopvvvmGBg0acOCBB243UlhtDBgwgKFDh/LRRx9x3nnnVTgfDuEc/Lx58zj88MOBcDj/b3/7G3vssUfZOj169KBXr15069aNTp06bfcfgY0bN3LooYeydetWbr7n5rL5Xbp0YdCgQSxdupR7772Xxo0bb1e+4uLisqFSmzZtyoQJE1i+fDm/+MUveOWVV9hvv/248sorufrqq5kwYcIO/x2A0OzPxAQ8DRxXad5YYEi51/OBvavbTtP9m3o+6TO2T7aLkFb5VJ98qou76pPLoq7Le++9F+n2K5u7bG5k2/7LX/7iV1xxRWTbd3dv3769f/XVV2WvS+tz0UUX+WOPPRbpvt2Tf17ADE+SiRm5sM3MOgC9gMqX4u0LfFbu9eLEPBEREamBeZIT8GndgdmuwGvAbe7+ZKVlzwO3u/sbidevAD9395mV1hsODAewVtan9x35c+p83vJ5HNzm4GwXI23yqT75VBdQfXJZ1HX5Xc/fsVfHvSLbfmUbNm+gcYPtDzPHVabr8+UnX/LzWT+vMG/mZTNnunvhdisna56nawIaAv8C/quK5TqcnkeHBN3zqz75VBd31SeX6XB6bst0fXLicLqF6+nvA+a5+/9UsdozwIWJq9QPA7519yVRlUlERCSfRHl1en9gKDDbzGYl5v0S2B/A3e8F/gmcBHwErAMujrA8IiIieSWyEPdwnttqWMeBK6Iqg4iI7JjSnll3oF8XySB1uyoiIhlXXFzM6NGjq1z+1FNPpW2I0nymEBcRkaTKdVeecQrx1CjERUQkqUWL0ru92267jS5dunDssccyf/58AP70pz/Rt29fevTowVlnncW6deuYOnUqzzzzDD/72c/o2bMnH3/8cdL1RCEuIiIZMHPmTB5++GHeeecdnnzySd5++20AfvCDH/D222/z7rvvcvDBB3PfffdxxBFHcNppp3HnnXcya9YsDjjggKTriUJcRETKKS4GM3jttfDaLEzFxTu33ddff50zzzyTpk2b0qJFi7KxwefMmcORRx5J9+7defDBB5k7d27S96e6Xl2jEBcRkTLFxeAOgwaF1+5h2tkQh+2H8QQYNmwY99xzD7Nnz2bUqFEVhgPdkfXqGoW4iIhEbuDAgfzjH/9g/fr1rF69mmeffRaA1atXs/fee7Np06ay4TsBmjdvzurVq8teV7VeXacQFxGRpHZiyO/t9O7dm8GDB9OzZ0/OOussjjzySABuueUWDj30UI477ji6du1atv4Pf/hD7rzzTnr16sXHH39c5Xp1XZQ9tomISIx16JDe7d1www3ccMMN280fOXLkdvP69+9f4RazkSNHJl2vrlOIi4jIdtRTWzzocLqIiEhMKcRFRERiSiEuIlIHhPGmJNfV9nNSiIuI5LnGjRuzYsUKBXmOc3dWrFhB48aNU36PLmwTEclz7dq1Y/HixXz11VcZ2d+Xq7/Ellc7EnWsZLI+jRs3pl27dimvrxAXEclzDRs2pGPHjhnb39BxQ5kxfEbG9he1XK6PDqeLiIjElEJcREQkphTiIiIiMaUQFxERiSmFuIiISEwpxEVEpFaKisIk2acQFxERiSmFuIiISEwpxEVERGJKIS4iIhJTCnEREZGYUoiLiIjElEJcREQkphTiIiIiMaUQFxERiSmFuIiISEwpxEVERGJKIS4iIhJTCnERkTyiwUnqFoW4iIhITCnERUSkzor7kQuFuIiISEwpxEVERGJKIS4iIhJTCnEREZGYUoiLiIjElEJcREQkphTiIiIiMaUQFxERiSmFuIiISEwpxEVERGJKIS4iIhJTkYW4md1vZsvMbE4Vy1ua2bNm9q6ZzTWzi6Mqi4iISD6KsiU+HjihmuVXAO+5ew+gCLjLzBpFWB4REZG8ElmIu/sUYGV1qwDNzcyAXRPrbo6qPCIiIvmmQRb3fQ/wDPAF0BwY7O5bs1geERGRWMlmiH8fmAUcDRwAvGRmr7v7qsormtlwYDiAtTIKxxVmspyRmrd8nuqTo/KpLqD65LJ01mX+krEAFI67LC3b25F9xOmzSeXvldP1cffIJqADMKeKZc8DR5Z7PQnoV9M2m+7f1PNJn7F9sl2EtMqn+uRTXdxVn1yWzroMGhSmKNW0jzh9Nqn8vXKhPsAMT5KJ2bzF7FPgGAAz2xPoAizIYnlEpJaKisIkItkR2eF0M3uIcNV5GzNbDIwCGgK4+73ALcB4M5sNGHCduy+PqjwiIiL5JrIQd/chNSz/Ajg+qv2LiIjkO/XYJiIiElMKcRERkZhSiIuISK0tXJjtEggoxEVEZAcsWpTtEggoxEVERGJLIS4iIikpLgYzeO218NosTMXF0exP/RDUTCEuIiIpKS4Gdxg0KLx2D1NUIS41U4iLiIjElEJcRERqrX37bJdAQCEuIiI7oEOHbJdAQCEuIiISWwpxERGRmFKIi4iIxJRCXEREJKYU4iIiIjGlEBcREYkphbiIiEhMKcRFRERiqkG2C1BXlXXqf142SyEisO3f4+TJ2SyFSO2pJS4iIhJTCnEREZGYUoiLiIjElEJcREQkphTiIiIiMaUQFxERiSmFuIiISEwpxEVERGJKIS4iInXawoXZLsGOU4iLiEid1XTzKhYtynYpdpxCXEQkz8S5ZZkxs2fDyJE88eY+2S7JTlGIi4jkmTi3LCP13Xfw8MMwcCDFBU9g946h6dY1AJiFqbg4u0WsLYW4iIjkt08/hRtvhP32gyFD4IsvKB7dHF++gkGDwiruYVKIi4hIxhUXh5bka6+F13FtWabN1q0wcSKccQZ07Ai33w6HHQYvvggffADXXgutW2e7lDtNQ5FKLGioSJHqFReHqagoBLl7lguULV9/DePHw5gx8OGH0LYtXHcdXHYZtG+f9C1VzI4FhbiIiMTfzJnwxz/CQw/B+vXQv3/4X81ZZ8Euu1T71g4dMlLCSCjERUTyTJxblrWyfj08+mgI7+nToVkzuPBCGDkSevTIdukyIqUQN7N6QA9gH2A9MNfdl0ZZMBER2TFxblmmZMECuPdeuO8+WLkSunaF3/8+BHjLltkuXUZVG+JmdgBwHXAs8CHwFdAY6Gxm64CxwAR33xp1QUVEpA7bsgVeeCG0ul98EerVgzPPhMsvDxcCmGW7hFlRU0v8VmAMcJl7xcskzGwP4DxgKDAhmuKJiEhd1vK7r+C394WW96JFsPfeMGoUXHop7LtvtouXddWGuLsPqWbZMuDudBdIRETquC1b4NVXuWHeeKYs6wJv3gRHHQWjR8Ppp0PDhtkuYc5I9Zz4OcCL7r7azG4EegO3uvt/Ii2dSIbpVjaRLJo/HyZMgL/+FRYv5vD6LTmev3HLe2fDwQdnu3Q5KdXOXn6VCPABwPcJh8/HRFcsERGpE1auDPd0H3ZYuEDtjjugoAAeeYQfHPFlWEcBXqVUQ3xL4vFkYIy7Pw00iqZIIiKS1zZtgueeg3POCee4L78c1q0Lh8s//5zivs9jg8/lpdcbA+p9rjqp3if+uZmNJVylfoeZ7YK6bBURkdqYNQseeAAefBCWLQu9qY0cCRddBD17ll1hrt7nUpdqiJ8LnACMdvdvzGxv4GfRFUtERPLC0qUhtCdMgJKScFHaqaeG4D7xRF2ktpNSDfEC4CV3X514vQb4NpoiiYhInDXatBUeeywE94svhqvN+/WDP/wBBg+u1cAjdab3uR2UaoiPIVyRXmptknkiIlIHJL17wx3eegsmTODFB0pg3bnhPu6f/Sz0pLaDF6flfe9zOynVELfynb24+1YzU7/rIiJ13WefhVvCJkwIQ3w2acIb3Vty0q0Pw9FHQ/362S5hXkv14rQFZvYTM2uYmK4GFkRZMBERyVGrV4fgPvbYcLz7hhtgr71CX+ZffslNP+oIxx2nAM+AVFvTI4DfAzcCDrwCDK/uDWZ2P3AKsMzdD6linSJCr28NgeXuPijF8oiISCatWgXPPhvOdb/4ImzcCJ06hS5Qhw4NzyXjUgrxRBerP6zltscD9wAPJFtoZq2APwInuPunib7YRUQkV3z7LTzzDDz+OPzrXyG4990XRowI93gfcUSdHXgkV6Ta7WpnwoVse7r7IWZWAJzm7rdW9R53n2JmHarZ7HnAk+7+aWL9ZakXW0REIvHNNyG4H3sMJk6E776Ddu3C/dznnBN6VqunbkJyRaqH0/9EuC98LIC7l5jZ3wmjnO2ozkBDM5sMNAf+z92rarUPJ3H43loZheMKd2K3uWH+krEAbF0+Ly/qU2peRPUp/XsVjrss7duuaj9R1SVboqhPpj6XZNJZn2zWA7Jfl+ZrNzPo3W84dubXHDpvNQ23OEt2b8QrA1vxcp/dmNuhGV7vdZjzOsypeXvpqk8mPpdU9pHTvwXuXuMEvJ14fKfcvFkpvK8DMKeKZfcA04BmQBvCeOWda9pm0/2bej4YNChMfcb2yXZR0iqq+pT+vaJWfj/6bGo2aJB7+/Zp32xK0lmfTH2/qpKVuqxY4X7ffe4nnujesKE7hA/z2mvdp01z37p1h8uQrvpk4nNJZR+58FsAzPAkmZhqS3y5mR1AuKgNMzsbWLKT/39YTLiYbS2w1symAD2AD3ZyuyKSQYsWZbsEkrIVK+Cpp8Kh8ldegc2bw43Y11wTDpUXFuocd8ykGuJXAOOArmb2OfAJcMFO7vtp4J7E/eaNgEOB/93JbYqISHnLl8M//hGCe9Kk0Htax47wX/8VgrtPHwV3jKV6dfoC4FgzawbU823dr1bJzB4CioA2ZrYYGEW4lQx3v9fd55nZi0AJsBX4s7uncLZFRLKtuBh+/ettr0szYNQojTSVC1p+9xWMfTJcVf7qqyG4Dzgg9J52zjnQq5eCO0+kenX61cBfgNXAn8ysN/ALd59Y1XvcfUhN23X3O4E7UyyriOQIjTKVY9xh/nx47jn+d9bzFHw7Bd7cCgceCD//eQjucqOESf5I9XD6Je7+f2b2fWAP4GJCqFcZ4iIiEqGNG2HKFHj++TA298cfA9CiWXf+vv/1DH3mHCgoUHDnuZT7Tk88ngT8xd3fNdM3Q0Q0ylRGffkl/POfIbgnToQ1a2CXXeCYY+Daa+Hkk/nRhfsDMLRHlssqGZFqiM80s4lAR+B6M2tOOI8tInWcRpmK0Nat8M47oaX9/PPw9tth/r77wvnnw8knh0FGmjXLbjkla1IN8R8BPYEF7r7OzFoTDqmL5J2FC7NdAqnT1qyBl18Oof3887BkSTgkfuihcOutIbh79NBhcgFqCHEz6+DuC919K/Cf0vnuvgJYkTikvq+7L464nCIZo/ueJdP2Wb4R/t//C6H96quhq9MWLeD734dTToETT4S2bbNdTMlBNbXE7zSzeoR7umcCXwGNgQOBo4BjCLeOKcRFRFK1eTNMnVp2mPyZ994DfgKdO8OVV4bgHjAAGjbMdkklx1Ub4u5+jpl9DzgfuATYG1gHzAP+Cdzm7hsiL6VIxJLf9zyD4i9037OkyZdfhsPk//wnvPBCGGikYUMYOJC7uq3i2tsmwUEHZbuUEjM1nhN39/eAGzJQFpGsSXbfc+G4QoqHz8hyySS2Vq0KX6aXXw5dnM6dG+bvsQeccUZobR93HLRowUPjCrlWAS47INUL20QkTYqKwuPkydkshaTdxo0wbdq20J4+PfSU1rgxHHkkXHhhuBWsVy8N5SlpoxAXqUT3PUtKtm6FWbNCYL/8Mrz+OqxfHwK6Xz/4xS9CaB9+eAhykQgoxPOcWn21p/ueJSn30CtaaUt70iRYuTIs+9734NJL4dhjYdAgaNkyu2WVOqOmW8x6V7fc3f9T3XIRkVj78ssQ1qWt7U8/DfPbtYPTTgst7aOPhn32yW45pc6qqSV+V+KxMVAIvEvogrUAeAsYEF3RJC4y1dpXJywSudKL0V55JUxzEgMrtmoVwvq660Jr+6CD1NmK5ISabjE7CsDMHgaGu/vsxOtDgP+Ovngi26gTFkm7lSth6lR+vODf/HvpgbD7ZdsuRhswAC64YNvFaPXrZ7u0IttJ9Zx419IAB3D3OWbWM5oiiYhEwD0cznnjDfj3v8Nj4ravwdaAC3wTf/zlgtDS1sVoEhOphvg8M/sz8DfAgQsIHb6IRCp5JywwapQ6YZEabN4MJSUhrEuD+4svwrIWLeCII2DIEBgwgJNv7AtvALfdltUii9RWqiF+MTASuDrxegowJpISiZSTrBMWkaTWrIG33toW2tOmhXkA++8frhofMCBM3bpB/frhP4lF2zah/yRK3KQU4u6+wcz+ALxMaInPd/dNkZasDli4ENpkuxAiMdX6203w+OPbQnvWrHA+2wwKCkLnKgMGQP/+IcST0H8SJe5SCnEzKwImAAsJV6fvZ2YXufuUyEpWByxapBCvDXXCUodt3Qrz528L7Dfe4F8LFgDnQJMmYZjO668PgX344bpPW+qMVA+n3wUc7+7zAcysM/AQ0CeqgolUpk5Y6pClS+Htt2HGjPA4bdq2jlXatoUBA/jf3t/x0589AT17QqNGO71L/SdR4ijVEG9YGuAA7v6BmWmMvB1Q+UKtmZfNwC7TOTipw77+GmbODGFdGtyffRaW1asHBx8cBgwpPZ994IFgxoPjCvlpv35pK4b+kyhxlGqIzzCz+4C/Jl6fTxhfXGqp8jm4PmMLmaGRsqSuWLMG3nmnYmB/9NG25QceGIK6sBD69g33Z++6a/bKK1ml7qJrlmqIjwSuAH5COCc+BfhjVIUSkTywYUO4xav8YfF588L5bYD99gtBfckl4bFPH9htt+yWWSRmUr06fSPwP4lJ0kTn4CQqGR/4ZvPm0HFK+cCePRs2JW5iads2BPXZZ4fHwkLYc88MFU6kanFv7dc0AMqj7n6umc0m3FpWgbsXRFayOqBDB1iT7UKI1Na6dSGwS0q46qN36bJ6BjR/J7S8IVwZXlgI1167LbD32099jYtEoKaWeGnnLqdEXRARyTFbt4b7IEtKKk4fflh2Q/VJ9Zrx4a69YOTIbeexDzggXJAmIpGraQCUJWZWH7jP3Y/NUJlEJNNWrQqHv8uH9ezZsHp1WG4WwrmgAM47LzwWFHDSJR1xq8dknWgTyYoaz4m7+xYzW2dmLd3920wUKtsyfj5RJFO2bAlXg1duXZcf57VVqxDSF11UFtZ065b0KnHXEXKRrEr16vQNwGwzewlYWzrT3X8SSalEZKe12LQCXt0W1BNemQfXNIf168MK9etDly5w2GEwfPi2wG7XTuevRWIi1RB/PjGJSC7ZsiWct37//W3T/Pk8+eb77L5pGRydWK9NG9a0qQ8jLg1B3aNH6ERFw22KxFqqt5hNMLMmwP7le24TyTc5ewplzZrQd3j5sH7//XCR2caN29Zr0wa6dOHN1qfyadOujByTaF3vuSdX/KkvM4br5HW+y9nvsEQi1QFQTgVGA42AjmbWE7jZ3U+LsGwidYs7LF6cPKw//3zbevXqhYvMunaFE04Ij127hkPjbcKQOncWhVVHHp/5aohI5qR6OL0Y6AdMBnD3WWbWMaIyieS1XbasY9/1H8Fj87c7DM7atdtWbNEihPMxx2wL6a5dQ4Dvskv2KiAiOSPVEN/s7t9axYtdNPKuSDJbt8KSJbBgwfbTJ5/wryVLwnrnJtZv3z6E85FHVmxV77WXLjATkWqlGuJzzOw8oL6ZHUToQ31qdMUSyXGrV8Mnn1QI5wrPy5+nNgs9lnXqBCeeyJ8ndeKLJgdw00MHw0EHQdOm2auHiMRaqiF+FXADsBH4O/Av4NaoCiWSdVu2sNeKjfDqq9uH9IIF8NVXFddv0SIc5u7WDU49NQR26bT//hXGu/5bUXi8qUfmqiMi+ammvtMbAyOAA4HZwOHuvjkTBROJjDt8+224iOyzz8Jj6fTZZyGwFy3iuU2bKLtHq379cNi7Uyc488xtAd2xY3jcbTcd+haRjKupJT4B2AS8DpwIHAxcE3GZRHacO3zzTcVwThbU5S8ggxDAe+8dOjrp0wfOOYdbP3uQG4fdH0J6v/2gQaoHrkREMqOmX6XvuXt3ADO7D5gefZFEquAOK1ZWHcylz9etq/i+evVgn31CQB9ySLgta7/9wut27cLzvfaChg0rvO2pcRO58ZhjMlhBEZHaqSnEN5U+cffNpsOFkm6lLeelS7efli0re/7Qu0vZ/bsvoc2Giu+vX39bQPfoAaecUjGc27ULAa1WtIjkoZp+2XqY2arEcwOaJF4b4O7eItLSSTxt2QIrVlQbyhXmbdq0/Tbq1YO2bWHPPWHPPZnd4iC+brQng69tV7EVvddeIchFROqgmoYi1a9jXbdhA3z9Naxcue2x0vMb531Nq+++goKl/GvhPBjZKNwrXVnDhmWhzF57lXUHut20xx7QunWFcP5NUXgcfE1Gai0iEgs6xlgHmG+Fb1YlDeCaArpsxKtk6tWDVq3oun53VjVoDR078lrLJfygaMT2obznnmGIS52SERFJG4V4Ltq6NVyctXp1GPhi9eqKz2vx+MwXq2m2+VvYLUnLuFTTpuEWqd13D9MBB0DfvuF5+fmlz0sfW7SAevW4oChsZvLT8Jtxhfxg+C0Z+TOJiNR1CvHquMN334VDyhs3pvXx13M30mjrehqNng/39q4YwGvWpF7Gxo1h112hefNtj61ahfPFzZvz8sRdWd1gNy66JkkYlz5XP9wiIrGkEE/iljln0n/F01AvTd3Dm4WpWbPQem3QgH4rv8Kpx5b166DJonD+t1ev0OPXxo3w4othXunUoAFcemkYDGPRIrj11jC//OHp3/wGjjgCpk6FX/4Sli+nU+nt0E8Dd98NPXvCyy/D5ZdvX86xY0Of3c8+C3fdtf3yv/41XFT2yCMwZkzZ7LtnJZ4sfzw8jh8fpsr++c/Q6v/jH+HRR7dfXjqG4ujR8NxzFRbdUdKE6wpeCC9uuQVeeaXie1u3hieeCM+vvx7efLPi8nbt4G9/C8+vuQZmzaq4vHNnGDcuPB8+HD74gLFL5sPfi8K8nj3D3w/gggvCrWzlHX443H57eH7WWeHCvvKOOQZ+9atEXU5kl63roajc8lNOgf/+7/C8qIjtnHtu+MzWrYOTTtp++bBhYVq+HM4+e9tnUrqpkSPD42efwdCh27//2mtDT3Pz58Nll22//MYb4dhjw9/tmmvKZpftZ2ql715l5b97tybp7HEHvnsVPp/HHw8juO3ody+M7ZT0u0eTJvBCxN+9wsTzxHevgjR+9zjxxO1PkaX5uweVPpuRI2Hw4LR/98pU/t2rLNXvXkwpxJOY1vokOq+ZyZ571QvnfevVC2F56KHhHmMz+O1vty0rXX7WWTBkSOhI5LLLti0rVe7LPL97+DJv2XMmffbuHpZfeeW2L/PcudsXrGtXOPjgEPK6ZUpEpM4z93gNRtasfTNfu2htzSvuhNL/jJY2DKPcx5rzCpkxfEbk+8lEXSZPhsJx0dQnE/WoLM51SbaPKOqTjc+lVDrrk816QHTftWxRfdLPzGa6e2Hl+fWSrZymHd5vZsvMbE4N6/U1sy1mdnZUZREREclHkYU4MB44oboVzKw+cAdhVDSJyMKF2S6BVKbPRETSIbIQd/cpwMoaVrsKeAJYFlU5JFwHJ7lFn4mIpEPWro4ys32BMwljPfatYd3hwHAAa2UUjtvutEBazV8yFoDCcUmulEzzPrYunxdpfcJ++kS+j43L96Zw3GnMi6g+mfhMKou2LtF/JlDx7xVFfbLxuZRKZ32yWQ+I7ruWLapPBrl7ZBPQAZhTxbLHgMMSz8cDZ6eyzab7N/WoDRoUpkzso8/YPpFsf9Qo93Cje8Vp1Kj072vQoLBt9+jqk4nPpLJ01yXTn0nlv1cUn002PpdS6axPNuvhHt2/m2xRfdIPmOFJMjHKc+I1KQQeNrOFwNnAH83sjCyWJ68UF4eIGDQovC6NjOLibJaqbtNnIiLplrXD6e7esfS5mY0HnnP3p7JVHqm94mL49a+3vQ79zsyg+Iv0B1O2bv0REcllUd5i9hDwJtDFzBab2Y/MbISZjYhqn5Jc+/bRbDdZy7LP2EK1LFMQ1WciInVLZC1xdx9Si3WHRVUOgQ4dsl0CqUyfiYikg/rulLRQy7JuysRpjmz3piaSy7J5YZvkEbUsRUQyTyEuIiISUwpxERGRmFKIi4iIxJRCXEREJKYU4lXQKFMiIpLrFOJV0ChTIiKS6xTiIiIiMaUQL6e4OPT//dpr4bVZmNSNqIiI5CKFeDkaZUpEROJEIS4iIhJTCvEqqC9wERHJdQrxKqgvcBERyXUKcRERkZjSUKQiUudlaphTDasq6aaWuIiISEwpxEXylLoOFsl/CnGRPKWug0Xyn0JcREQkphTiInlEXQeL1C0KcZE8oq6DReoWhbiIiEhMKcSzZPJk3Ssq0VLXwSL5TyEukqfUdbBI/lOIi4iIxJRCXEREJKYU4iIiIjGlEBcREYkphbiIiEhMKcRFRERiSiEuIiISUwpxERGRmFKIi4iIxJRCXERy3sKF2S6BSG5SiItIzlu0KNslEMlNCnEREZGYUoiLSE4qLgYzeO218NosTBobXWSbBtkugERLw53mHn0mqSkuDlNRUQhy9ywXSCQHKcRlpymURESyQ4fTRSTntW+f7RKI5CaFuIjkvA4dsl0CkdykEBcRySDd8y7ppBAXEckg3fMu6aQQFxERiSmFuIhIxHTPu0RFIS4iErHi4nCf+6BB4bV7mBTisrN0n3gSuu9ZRETiILKWuJndb2bLzGxOFcvPN7OSxDTVzHpEVRYRkVyhe94lnaI8nD4eOKGa5Z8Ag9y9ALgFGBdhWUREcoLueZd0iuxwurtPMbMO1SyfWu7lNKBdVGURERHJR7lyYduPgBeyXQgREZE4yfqFbWZ2FCHEB1SzznBgOIC1MgrHFWaodNGbt3ye6pOj4lyX+UvGAlA47rKyeXGtT7K6QDzrk091qY7qk0HuHtkEdADmVLO8APgY6JzqNpvu39TzSZ+xfbJdhLTKp/rEuS6DBoWpvLjWJ1ld3ONZn3yqS3VUn/QDZniSTMza4XQz2x94Ehjq7h9kqxwiIiJxFdnhdDN7CCgC2pjZYmAU0BDA3e8FbgJaA380M4DN7p6jxytERERyT5RXpw+pYfmlwKVR7V9ERCTf5crV6SIiIlJLCnEREZGYUoiLiIjElEJcREQkphTiIiIiMaUQFxERiSmFuIiISEwpxEVERGJKIS4iIhJTCnEREZGYUoiLiIjElEJcREQkphTiIiIiMaUQFxERiSmFuIiISEwpxEVERGJKIS4iIhJTCnEREZGYapDtAohI+k2enO0SpI/qIlI1tcRFRERiSiEuIiISUwpxERGRmFKIi4iIxJRCXEREJKYU4iIiIjGlEBcREYkphbiIiEhMKcRFRERiSiEuIiISUwpxERGRmFKIi4iIxJRCXEREJKYU4iIiIjGlEBcREYkphbiIiEhMKcRFRERiSiEuIiISUwpxERGRmFKIi4iIxJRCXEREJKYU4iIiIjGlEBcREYkphbiIiEhMKcRFRERiSiEuIiISUwpxERGRmFKIi4iIxJRCXEREJKYiC3Ezu9/MlpnZnCqWm5n93sw+MrMSM+sdVVlERETyUZQt8fHACdUsPxE4KDENB8ZEWBYREZG8E1mIu/sUYGU1q5wOPODBNKCVme0dVXlERETyTTbPie8LfFbu9eLEPBEREUlBgyzu25LM86Qrmg0nHHLHWhmF4wqjLFdGzVs+T/XJUflUF1B9clk+1QVUn4xy98gmoAMwp4plY4Eh5V7PB/auaZtN92/q+aTP2D7ZLkJa5VN98qku7qpPLsunurirPlEAZniSTMzm4fRngAsTV6kfBnzr7kuyWB4REZFYiexwupk9BBQBbcxsMTAKaAjg7vcC/wROAj4C1gEXR1UWERGRfBRZiLv7kBqWO3BFVPsXERHJd+qxTUREJKYsNIjjw8xWEy6CyxdtgOXZLkQa5VN98qkuoPrksnyqC6g+UWjv7m0rz8zmLWY7ar675+i1/rVnZjNUn9yUT3UB1SeX5VNdQPXJJB1OFxERiSmFuIiISEzFMcTHZbsAaab65K58qguoPrksn+oCqk/GxO7CNhEREQni2BIXERERcjjEzewEM5tvZh+Z2S+SLDcz+31ieYmZ9c5GOVORQl3OT9ShxMymmlmPbJQzVTXVp9x6fc1si5mdncny1VYq9TGzIjObZWZzzey1TJexNlL4vrU0s2fN7N1EfXK2t0Qzu9/MlpnZnCqWx+l3oKa6xO13oNr6lFsvLr8DNdYnJ38HknWonu0JqA98DHQCGgHvAt+rtM5JwAuE0dAOA97Kdrl3oi5HALslnp+Yq3VJtT7l1ptE6F737GyXeyc/n1bAe8D+idd7ZLvcO1mfXwJ3JJ63BVYCjbJd9irqMxDoTdUDKcXidyDFusTmdyCV+iTWicXvQIqfT07+DuRqS7wf8JG7L3D374CHgdMrrXM68IAH04BWZrZ3pguaghrr4u5T3f3rxMtpQLsMl7E2UvlsAK4CngCWZbJwOyCV+pwHPOnunwK4ey7XKZX6ONDczAzYlRDimzNbzNS4+xRC+aoSl9+BGusSs9+BVD4biM/vQCr1ycnfgVwN8X2Bz8q9XpyYV9t1ckFty/kjQssiV9VYHzPbFzgTuDeD5dpRqXw+nYHdzGyymc00swszVrraS6U+9wAHA18As4Gr3X1rZoqXdnH5HaitXP8dqFHMfgdSkZO/A7naY5slmVf5MvpU1skFKZfTzI4i/OMdEGmJdk4q9bkbuM7dt4TGXk5LpT4NgD7AMUAT4E0zm+buH0RduB2QSn2+D8wCjgYOAF4ys9fdfVXEZYtCXH4HUhaT34FU3E18fgdSkZO/A7ka4ouB/cq9bkdoNdR2nVyQUjnNrAD4M3Ciu6/IUNl2RCr1KQQeTvzDbQOcZGab3f2pjJSwdlL9ri1397XAWjObAvQAcjHEU6nPxcBvPZzY+8jMPgG6AtMzU8S0isvvQEpi9DuQijj9DqQiJ38HcvVw+tvAQWbW0cwaAT8Enqm0zjPAhYmrUw8DvnX3JZkuaApqrIuZ7Q88CQzN9v/qUlBjfdy9o7t3cPcOwOPA5Tn8DzeV79rTwJFm1sDMmgKHAvMyXM5UpVKfTwmtCcxsT6ALsCCjpUyfuPwO1ChmvwM1itnvQCpy8ncgJ1vi7r7ZzK4E/kW4uvF+d59rZiMSy+8lXO14EvARsI7Qusg5KdblJqA18MfE/1o3e452tp9ifWIjlfq4+zwzexEoAbYCf3b3am+ryZYUP59bgPFmNptwOPo6d8/2CE1JmdlDQBHQxswWA6OAhhCv3wFIqS6x+R2AlOoTKzXVJ1d/B9Rjm4iISEzl6uF0ERERqYFCXEREJKYU4iIiIjGlEBcREYkphbiIiEhMKcRFcpCZnWlmbmZd07jNIjN7LvH8NEuMcGZmZ5jZ93Zge5PNrFa3QJnZ3WY2sLb7qrSNNYnHtolbfkTqLIW4SG4aArxB6Kwl7dz9GXf/beLlGUCtQ7y2zGx34LDEQBOVl9Wv7fbc/StgiZn1T0f5ROJIIS6SY8xsV6A/of/sH5abX2Rmr5nZo2b2gZn91sIY1NPNbLaZHZBYb7yZ3WtmryfWOyXJPoaZ2T1mdgRwGnBnYpzkA8q3sM2sjZktTDxvYmYPWxjv+hFC/9Gl2zvezN40s/+Y2WOJOlR2NvBiufcsNLObzOwN4Bwz+7GZvW1hnPMnEr1ikeh97s3EslsqbfMp4Pxa/5FF8oRCXCT3nAG8mOh6c6WZ9S63rAdwNdAdGAp0dvd+hP62ryq3XgdgEHAycK+ZNU62I3efSui69Gfu3tPdP66mXCOBde5eANxGGAwCM2sD3Agc6+69gRnAfyV5f39gZqV5G9x9gLs/TBjmsa+79yB0Z/mjxDr/B4xx977Al5XePwM4spoyi+Q1hbhI7hlCGAecxOOQcsvedvcl7r4R+BiYmJg/mxDcpR51963u/iGhX/R0nFsfCPwNwN1LCN1PAhxGOBz/bzObBVwEtE/y/r2BryrNe6Tc80MSRw9mE1rX3RLz+wMPJZ7/tdL7lwH71LomInkiJ/tOF6mrzKw1YYjQQ8zMCf2fu5n9PLHKxnKrby33eisV/z1X7k+5Nv0rb2bbf/Art+CTbceAl9x9SJJl5a1Psr215Z6PB85w93fNbBihH+vq9ltavvU17Fckb6klLpJbzgYecPf2iRGg9gM+ofZjS59jZvUS58k7AfOrWXc10Lzc64UkDpUnylNqConzz2Z2CFCQmD8N6G9mByaWNTWzzkn2Mw84sJpyNCdcqNaQiue5/822awMqn//uDGR9EAqRbFGIi+SWIcA/Ks17AjivltuZD7wGvACMcPcN1az7MPAzM3snEfqjgZFmNpUwDnSpMcCuZlYC/JzE+OOJq8SHAQ8llk0j+eH756nYuq7sV8BbwEvA++XmXw1cYWZvAy0rveeoxHZF6iSNYiaSZ8xsPPCcuz+e7bJUlrgS/RR3/yZN25sCnO7uX6djeyJxo5a4iGTStcD+6diQmbUF/kcBLnWZWuIiIiIxpZa4iIhITCnERUREYkohLiIiElMKcRERkZhSiIuIiMSUQlxERCSm/j8LYG324prdawAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "#\n", "# Expected period and chi2 for SHM\n", "ySHM = 2*np.pi*np.sqrt(L/g)*np.ones(nTheta)\n", "chi2SHM = (ySHM - meanPeriods)**2/errPeriods**2\n", "chi2sumSHM = np.sum(chi2SHM)\n", "print(\" \")\n", "print(f\"Chi2 for simple approximation = {chi2sumSHM:.2f}, chi2/NDF = {chi2sumSHM/(nTheta - 1):.2f}\")\n", "#\n", "# Expected period and chii2 for Lima formula\n", "chi2 = (periodLima(L, thetaArr) - meanPeriods)**2/errPeriods**2 \n", "chi2sum = np.sum(chi2)\n", "print(\" \")\n", "print(f\"Chi2 for Lima/Arun approximation = {chi2sum:.2f}, chi2/NDF = {chi2sum/(nTheta - 1):.2f}\")\n", "#\n", "print(\" \")\n", "plt.figure(figsize = (8, 6))\n", "plt.title('Pendulum data', fontsize = 14)\n", "plt.xlabel('Amplitude (rad)')\n", "plt.ylabel('Period (secs)')\n", "plt.errorbar(thetaArr, meanPeriods, yerr = errPeriods, \n", " linestyle = '', marker = '+', color = 'b', label = \"data\")\n", "plt.plot(thetaArr, periodLima(L, thetaArr), color = 'r', label = \"Good approx.\")\n", "plt.plot(thetaArr, ySHM, color = 'r', linestyle = '--', label = \"Simple approx\")\n", "plt.xlim(0.0, 1.1*thetaArr[nTheta - 1])\n", "plt.ylim(0.8*np.min(meanPeriods), 1.2*np.max(meanPeriods))\n", "plt.legend()\n", "plt.grid(color = 'g')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Does the difference between the two $\\chi^2$ values calculated above allow us to state that the pendulum is not following SHM? In this case, no, as acceptable $\\chi^2$/NDF values lie in the range 0.2 to 5.0 (some people use 0.25 to 4) so both of the above values are OK: we cannot say the pendulum is not obeying the SHM prediction." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Increase precision\n", "\n", "Repeat the above calculation using the large dataset in the file _periodTableLarge.csv_. Can we now conclude that there is a deviation from SHM at large amplitudes? " ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " \n", "Shape of periodArrLg = (16, 1000)\n", "nTheta = 16 and nMeas = 1000\n", " \n", "Number\t Amplitude (rad)\t Mean period (s)\n", "0\t 0.10\t\t\t 1.40 +- 0.02\n", "1\t 0.20\t\t\t 1.45 +- 0.02\n", "2\t 0.30\t\t\t 1.45 +- 0.02\n", "3\t 0.39\t\t\t 1.44 +- 0.02\n", "4\t 0.49\t\t\t 1.46 +- 0.02\n", "5\t 0.59\t\t\t 1.44 +- 0.02\n", "6\t 0.69\t\t\t 1.45 +- 0.02\n", "7\t 0.79\t\t\t 1.48 +- 0.02\n", "8\t 0.88\t\t\t 1.49 +- 0.02\n", "9\t 0.98\t\t\t 1.53 +- 0.02\n", "10\t 1.08\t\t\t 1.54 +- 0.02\n", "11\t 1.18\t\t\t 1.53 +- 0.02\n", "12\t 1.28\t\t\t 1.56 +- 0.02\n", "13\t 1.37\t\t\t 1.61 +- 0.02\n", "14\t 1.47\t\t\t 1.66 +- 0.02\n", "15\t 1.57\t\t\t 1.69 +- 0.02\n", " \n", "Chi2 for simple approximation = 956.91, chi2/NDF = 63.79\n", " \n", "Chi2 for Lima/Arun approximation = 16.83, chi2/NDF = 1.12\n", " \n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfEAAAGECAYAAAAm3RkPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA7rElEQVR4nO3dd5xU1fnH8c8jRYoIClZU1kaR3jSKyoINoyKKUUEx2EBsGEtMNIaN5WeMYEmwgIpoNBq70ViwIRo0FEMVQSkqNlyQXqQ8vz/OLGzfWZi7M3f4vl+vee3M3Dt3nrNlvnvuPfcec3dEREQkfnZIdwEiIiKydRTiIiIiMaUQFxERiSmFuIiISEwpxEVERGJKIS4iIhJTCnGRLGNmnczMzSynEq/JTbymUYSlJW1r2iCyPVKIi0TEzEYngsjNbL2ZzTOzoWZWN921ZSMzG2tmw9Ndh0hVqp7uAkSy3NtAP6AGcBTwMFAXGJTOokQkO6gnLhKtde7+vbt/7e7/AJ4EegFY8Fszm2tma8xsupmdW/BCM8tJ9OJ7m9lbZrbazD41s+MKv4GZ9TCzz8xsrZl9ADQttry/ma0s9ly5u8+TeU3BOmZ2YuL9V5vZv8ysvpmdYWafm9kyM/u7mdUu75uURBsamtlTZrYw8b2aaWbnF1o+GugKXFZo70eOmVUzs0fMbH7idZ8nvuf67JOsoF9kkaq1htArB7gVuBC4DDgEuB0YYWYnFXvNbcBfgbbAROBpM9sJwMz2BV4C3gLaAX8D/hJpC4raEbgGOAc4BugEPAf8GuhN+IflZODSsjaQZBtqAZ8kttUSuJfwvTomsXww8BHwKLBX4vY14TPuG+BMoAVwI3ADcD4iWUC700WqiJkdCvQF3kkcF78aON7dP0isMj+xzmXAvwu99G53fyWxjRuA8whh9yFht/xXwJUeJkL4zMyaArdUQZMgfIZc5u6zE/X9A/gNsIe75yeeexnoBgwrYxsVtsHdvwHuLPSakWbWHegDvOPuy8zsZ2C1u39faL2NwB8LPV5gZh0Sr3tkaxstkikU4iLR6pHYLV2d0AN/GbiC0POuBbxhZoVnIaoBLCi2jWmF7n+b+Lp74msL4GMvOpPRR6kpPSnrCgI84Qfg+4IAL/TcIeVso8I2mFk14HfAWUBjwh6AmsDYigo0s0uAi4AmQG3C9/jLil4nEgcKcZFojQMGAOuBb919PYCZ7Z9YfgqhF1rY+rIeu7ubGWw5FGZJ1LCplPVqlLbiVrxmQ7HHTsn6nfIP3SXThmsJu+0HA9OBlcD/seWfmdI3bHYWcE/i9eOB5YQ9Hacl8Z4iGU8hLhKt1e7+RSnPfwqsA5q4+7vbsP1Pgd5mZoV6sr8ots6PQB0z29ndlyeea1fBdrfmNVsrmTYcCbzi7n+HMCiQMPhtaaF1fgaqlfK6/7r75lPPzOzAFNYuklYa2CaSBu6+AhgKDDWzC8zsIDNrZ2aXmNmASmzqQSAHuMfMmpnZGcAlxdb5L7AKuD3xPr0pZ6DZNrxmayXThjnAMWZ2pJk1B4YD+xdbZwFwaGJUeqPECPQ5QIfECPqDzewmwih2kaygEBdJn5uAPMKu3pmE0dm9gfnJbsDdvwJOB3oAUwmDyn5XbJ0lhNHjxxF2RQ9IvHd52630a7ZWMm0gjOSfALxOOESxinC6XmFDCb3xTwl7EvYDRgDPAP8gjOzPoewBdiKxY0XHkoiIiEhcqCcuIiISUwpxERGRmFKIi4iIxJRCXEREJKYU4iIiIjEVu4u91KxX09s0a5PuMlLm8yWfc/CuB6e7jJTJpvZkU1tA7clk2dQWUHuiMHny5Hx3363487EL8Rq71mDSpEnpLiNlOo3sxKQBak8myqa2gNqTybKpLaD2RMHMSr3ev3ani4iIxJRCXEREJKYU4iIiIjEVu2PiIiJSOevXr2fhwoWsXbu2St7vL+3+wqxZs6rkvapCVbanVq1a7LPPPtSoUdFswYFCXEQkyy1cuJB69eqRk5NDYj76SPmPTovdWkT+PlWlqtrj7ixevJiFCxey//7FJ+krnXani4hkubVr19KwYcMqCXDZemZGw4YNK7XHRCEuIrIdUIDHQ2V/TgpxERGJ3A8//EDfvn054IAD6NixI4cffjgvvvhiSradm5ubVdcPqQyFuIiIRMrd6dWrF0cffTTz5s1j8uTJPP300yxcuDDdpVXaxo0b011CEQpxERGJ1LvvvkvNmjW55JJLNj/XpEkTrrjiCiAcsz///PNp3bo17du357333iv3+TVr1nD22WfTpk0bzjrrLNasWVPq+95888107tyZVq1aMWDAANwdCD33q666iiOOOIJWrVoxYcIEAPLy8ujXrx/du3fn4IMP5qGHHgJgwn8m0K1bN/r27Uvr1q3LrOvUU0/l8ccfB2DEiBGcc845qf5WlqDR6SIi25OrroIpU1K7zXbt4J57ylw8c+ZMOnToUOby++67D4Dp06fz2WefcfzxxzNnzpwyn3/ggQeoU6cO06ZNY9q0aWVu+/LLL+ePf/wjAP369ePVV1/llFNOAWDVqlWMHz+ecePGccEFFzBjxgwApk2bxscff8yqVato3749J510EgATJkxgxowZ7L///gwbNqzUukaOHEmXLl02r/Pxxx8n/z3cSuqJi4hIlbrsssto27YtnTt3BuDDDz+kX79+ADRv3pwmTZowZ86cMp8fN24c5557LgBt2rShTZvSJ8V67733OOyww2jdujXvvvsuM2fO3LysT58+ABx99NEsX76cpUuXAqE3Xbt2bRo1akS3bt0299IPPfTQzad9lVXXHnvswc0330y3bt0YNmwYu+66ayq/baVST1xEZHtSTo85Ki1btuT555/f/Pi+++4jPz+fTp06AWzezV1cWc9DxaO4165dy6WXXsqkSZPYd999ycvLK3LqVvHXFzwu9XmHunXrJlXX9OnTadiwId9++2259aWKeuIiIhKp7t27s3btWh544IHNz61evXrz/aOPPponn3wSgDlz5vDVV1/RrFmzpJ6fMWMG06ZNK/GeBYHdqFEjVq5cyXPPPVdk+T//+U8g9Krr169P/fr1AXj55ZdZu3YtixcvZuzYsZv3FhRWVl0TJkzg9ddf53//+x9Dhw5l/vz5W/cNqwT1xEVEJFJmxksvvcRvfvMb/vKXv7DbbrtRt25d7rjjDgAuvfRSLrnkElq3bk316tUZPXo0O+64Y5nPDxo0iPPPP582bdrQrl07Dj300BLv2aBBAy6++GJat25NTk5OiTDeZZddOOKII1i+fDmjRo3a/Pyhhx7KSSedxFdffcVNN93E3nvvXWLbpdUFcPHFF/Poo4+y9957M2zYMC644ALeffddRowYAVBkYF+qKMRFRCRye+21F08//XSpy2rVqrU5CJN5vnbt2mVuq7Bbb72VW2+9tdRlvXv35vbbby/xfNOmTRk5cmSR5w7tcij9e/WvsK6pU6duvt+zZ0969uwJRBPeBbQ7XUREJKbUExcRke3K2LFjS30+Ly+vSutIBfXERUREYkohLiIiElMKcRERkZiKLMTNbJSZLTKzGWUsr29mr5jZVDObaWbnR1WLiIhINoqyJz4a6FHO8suAT929LZALDDOzmhHWIyIiaXLbbbfRsmXLzed2//e//wXgoosu4tNPP03Je+y0004p2U6cRDY63d3HmVlOeasA9Sxc424nYAmwIap6REQkPT766CNeffVVPvnkE3bccUfy8/P5+eefAXj44YfTXF1qbdiwgerVq+7Er3SeYjYc+BfwLVAPOMvdN5W2opkNAAYAWAOj08hOVVZk1Gblz1J7MlQ2tQXUnkwWdVv+0u4v+I9lX+871dZuWMunP27pXU/4bAI77rwjc5fP3bJSDVj641L69+rPtXnX0qpdKzrldKLPBX34eNzH7Fx/ZwbfOJi7br6L7xZ+x/W3Xk/3Ht158ekXeee1d/h53c9889U3nHT6SVx63aUAbPJNm9931PBRvPGvN1i/bj3H/PIYLr/+8hJ13nzdzcyYMoO1a9dy/MnHb17nuI7H0ePUHkz4T5j85Nbh4YIxN1xxAzvuuCNfzP6CxT8u5rc3/5bc43N58ekXGffWONatW8ea1Wu4+5G7uWnwTSz8ciG16tQib2geBzY7kL6/7Mu1Q67l0C6Hcvetd7PDDjsw+IbBJer6fsX39BvZL7lvtrtHdgNygBllLDsDuBsw4CBgPrBzRduss18dzyYdR3RMdwkplU3tyaa2uKs9mSzqtnz66adFn+jateTtvvvCslWrSl/+6KNh+Y8/llxWzMxFM4s8XrFihbdt29YPPvhgHzRokI8dO7ZQKV194sSJ7u4O+Guvvebu7r169fLjjjvOf/75Z58yZYq3bdvW3d0fffRR33PPPT0/P99Xr17tLVu23Pz6unXrurv7m2++6RdffLFv2rTJN27c6CeddJK///77JepcvHixu7tv2LDBu3bt6lOnTnV39yZNmvitt97q7u6PPfaYdz0utPHXv/61n3DCCb5x40afM2eON27c2NesWeOPPvqoN27cePP2Lr/8cs/Ly3N393feeWdz7TNmzPDmzZv7mDFjvF27dr5u3boSNbmX8vMK35tJXkompnN0+vnAC4n6vkiEePM01iMiIhHYaaedmDx5MiNHjmS33XbjrLPOKvWypTVr1qRHjzCUqnXr1nTt2pUaNWrQunVrFixYsHm94447joYNG1K7dm1OP/10PvzwwyLbGTNmDGPGjKF9+/Z06NCBzz77jM8//7zE+z3zzDN06NCB9u3bM3PmzCLH5gumKu3Tpw9TJk3Z/PyZZ57JDjvswMEHH8wBBxzAZ599trmmgqlHC09V2r17dxYvXsyyZcto2bIl/fr145RTTmHUqFHUrLntw8DSuTv9K+AY4AMz2wNoBsxLYz0iItuHMq5YBkCdOuUvb9So/OVlqFatGrm5ueTm5tK6dWsee+wx+vfvX2SdGjVqbJ4KdIcddmDHHXfcfH/Dhi1DpsqaRrSAu/P73/+egQMHllnP/PnzGTp0KBMnTmSXXXahf//+ZU5VWtb9wo8rmqq0YL3p06fToEEDfvjhhzJrq4woTzF7CvgIaGZmC83sQjO7xMwKrgR/C3CEmU0H3gGud/f8qOoREZH0mD17dpGe8JQpU2jSpMlWb++tt95iyZIlrFmzhpdeeokuXboUWX7CCScwatQoVq5cCcA333zDokWLiqyzfPly6tatS/369fnhhx94/fXXiywvmKr0n//8J207td38/LPPPsumTZuYO3cu8+bNo1mzZiXqKzxV6dixY2nUqBE777wzL7zwAosXL2bcuHFceeWVLF26dKu/BwWiHJ3ep4Ll3wLHR/X+IiKSGVauXMkVV1zB0qVLqV69OgcddFCJmcIq48gjj6Rfv3588cUX9O3bl06dig4KPP7445k1axaHH344EHbnP/HEE+y+++6b12nbti3t27enZcuWHHDAASX+EVi3bh2HHXYYmzZt4ubhN29+vlmzZnTt2pUffviBBx98kFq1apWoLy8vb/NUqXXq1OGxxx4jPz+f3/3ud7zzzjvsu+++XH755QwePJjHHntsq78PoAlQREQkYh07dmT8+PGlLis8GUlBzxlKTkZSeNnuu+/O8OHDS2yr8DqDBw9m8OCSI78LK+24fIHLLruMIUOGABQZad+lSxfuvvvuIuv279+/yKGBXXfdlZdffrnENufMmbP5/pVXXllubcnSZVdFRERiSj1xERGJjeK93igUHglfWHk993RRT1xERCSmFOIiIiIxpRAXERGJKYW4iIiUkJsbbpLZFOIiIlLl8vLyGDp0aJnLX3rppZRNUZrNFOIiIlKqMgZpVwmFeHIU4iIiUqovv0zt9m677TaaNWvGsccey+zZswF46KGH6Ny5M23btqV3796sXr2a8ePH869//YvrrruOdu3aMXfu3FLXE4W4iIhUgcmTJ/P000/zv//9jxdeeIGJEycCcPrppzNx4kSmTp1KixYteOSRRzjiiCPo2bMnd955J1OmTOHAAw8sdT1RiIuISCF5eWAG778fHpuFW7GroFbaBx98wGmnnUadOnXYeeed6dmzJwAzZszgqKOOonXr1jz55JPMnDmz1Ncnu972RiEuIiKb5eWBO3TtGh67h9u2hjiUnMYTwhXYhg8fzvTp0xkyZEiR6UC3Zr3tjUJcREQid/TRR/Piiy+yZs0aVqxYwSuvvALAihUr2GuvvVi/fv3m6TsB6tWrx4oVKzY/Lmu97Z1CXERESrUNU36X0KFDB8466yzatWtH7969OeqoowC45ZZbOOywwzjuuONo3rz55vXPPvts7rzzTtq3b8/cuXPLXG97pwlQRESkVDk5qd3ejTfeyI033lji+UGDBpV4rkuXLkVOMRs0aFCp623vFOIiIlJCoWm+JYNpd7qIiEhMKcRFRERiSiEuIrIdcPd0lyBJqOzPSSEuIpLlatWqxeLFixXkGc7dWbx4MbVq1Ur6NRrYJiKS5fbZZx8WLlzIjz/+WCXv9/2K77H8khd2iauqbE+tWrXYZ599kl5fIS4ikuVq1KjB/vvvX2Xv129kPyYNmFRl7xe1TG6PdqeLiIjElEJcREQkphTiIiIiMaUQFxERiSmFuIiISEwpxEVERGJKIS4iIhJTCnEREZGYUoiLiIjElEJcREQkphTiIiIiMaUQFxERianIQtzMRpnZIjObUcby68xsSuI2w8w2mtmuUdUjIiJSXG5uuMVVlD3x0UCPsha6+53u3s7d2wG/B9539yUR1iMiIpJVIgtxdx8HJBvKfYCnoqpFREQkG5m7R7dxsxzgVXdvVc46dYCFwEFl9cTNbAAwAMAaWMcOd3SIoNr0mJU/ixaNWqS7jJTJpvZkU1tA7clk2dQWiFd7Zg8bAUCzawaWuU4mtGfywMmT3b1TiQXuHtkNyAFmVLDOWcAryW6zzn51PJt0HNEx3SWkVDa1J5va4q72ZLJsaot7vNrTtat7kyblr5MJ7QEmeSmZmAmj089Gu9JFRCRNvvwy3RVsvbSGuJnVB7oCL6ezDhERkTiK8hSzp4CPgGZmttDMLjSzS8zskkKrnQaMcfdVUdUhIiJSXF4emMH774fHZuGWl5fOqiqvelQbdvc+SawzmnAqmoiISJXJywu33NwQ5BGO8Y5UJhwTFxERka2gEBcRke1akybprmDrKcRFRGS7VXPjGjrv+XW6y9hqkR0TFxERyWjvvMOjkwayrEZD8I/DyLaYUU9cRES2L4sXw/nnw7HHAvDQ/rfHMsBBIS4iItsLd3jqKWjRAv7+d/jd72i8ZDp3T+2e7sq2mnani4hIRiqYInTs2BRs7MsvYdAgeP116NwZ3noL2rZNwYbTSz1xERHJXhs3wj33QMuWMG4c3H03fPRRVgQ4qCcuIiLZaupUuPhimDgRTjwRHngg3ueTlUI9cRERyS5r1sDvfw+dOsGCBfCPf8C//511AQ7qiYuISAZbsKCSL3j3XRg4EL74Avr3h6FDoWHDCCrLDOqJi4hIxkp6mtAlS+CCC+CYY8Io9LffhkcfzeoAB4W4iIjEmTs8/XQ4bezxx+H662HatBDm2wGFuIiIZJSkpwn96is45RTo0wf22w8mTYI//xnq1KniitNHIS4iIhklLy90sLt2DY/dw21ziG/cCPfeC4ccAu+9B3fdFU4ba9cuPQWnkQa2iYhIfEyfDhddBBMmQI8e4bSxnJx0V5U26omLiEjG2nxW2Nq1cOON0KEDzJsHTz4Jr722XQc4qCcuIiIZLCeHcN3VAQPg88/h17+GYcOyftR5stQTFxGRjLTT+p+4bvZF0K1bOA4+ZgyMHq0AL0Q9cRERySwbN8Kjj/L4xBupv34xXHddGNW2HY06T5ZCXEREMsd778FvfgNTp7LrEUfA8Degfft0V5WxtDtdRETSb+5cOP106N4dli4NF3D58EMFeAUU4iIikj7LloXd5S1ahGPet90Gs2bBWWeFK7xIubQ7XUREqt6GDfDII3DTTZCfHyYrue022GuvdFcWK+qJi4hI1Xr77XC+9yWXQPPm4XKpo0YpwLeCQlxERKrGnDnQsyccdxysXAnPPRcukN6hQ7oriy2FuIiIVEpubrgl7aef4OqroWXLcOGWP/8ZPv0UevfWce9tpGPiIiISjQ0bYMQIGDIkzPd90UVwyy2wxx7prixrqCcuIiKp9+ab0LYtXH45tGkDn3wCI0cqwFNMPXEREam0BQvKWDBrFvf+7XOY0QMOPBBeeikcB9du80ioJy4iIpX25ZfFnli8GK68Elq3pu0XK2HoUJg5E049VQEeIfXERURk661fD/ffD3/6U7hwy8CB9DrwP7x9zTXprmy7oJ64iIgkJS8vdKrffz88NgOrWYO8q36Cjh1h6lS4/36W1quR1jq3JwpxERFJSl4euMOFHacA4BjetBl5r3QKl0xt1Sqt9W2PFOIiIpKc6dOhd28enpyYlOSee8JzJ5+s495pElmIm9koM1tkZjPKWSfXzKaY2Uwzez+qWkREZBsUTEjSpg28/Tajmwzh4H3XwuDBULNmuqvbrkXZEx8N9ChroZk1AO4Herp7S+BXEdYiIiKVNWcOnHtuuNLaa6/BjTfC/PmMzslj7wNqpbs6IcLR6e4+zsxyylmlL/CCu3+VWH9RVLWIiEglzJsXrqz2+ONQq1aYKvS666BRo3RXJsWYu0e38RDir7p7idEOZnYPUANoCdQD7nX3x8vYzgBgAIA1sI4d7siei+XPyp9Fi0Yt0l1GymRTe7KpLaD2ZLJMacuei9dx4Wvfc8r4fDZUM54/ejce67EnS3YuOtp89rARADS7ZmCp28mU9qRKJrRn8sDJk929U4kF7h7ZDcgBZpSxbDjwMVAXaAR8DjStaJt19qvj2aTjiI7pLiGlsqk92dQWd7Unk6W9LV9/7T5okHuNGu41a7pfcYX7N99s9ebS3p4Uy4T2AJO8lExM58VeFgL57r4KWGVm44C2wJw01iQisv347ju4/fYwSYk7XHgh3HAD7LtvuiuTJKXzFLOXgaPMrLqZ1QEOA2alsR4Rke3DokVwzTVwwAHhamvnnRcGsT3wgAI8ZiLriZvZU0Au0MjMFgJDCMfAcfcH3X2Wmb0BTAM2AQ+7e5mno4mIyDbKz4c774Thw2HtWujXD266KUxUIrEU5ej0PkmscydwZ1Q1iIgIYS7vu+6Ce++FVaugT58wx3fTpumuTLaRJkAREclWy5aFq6rddRcsXw5nnhnC+5BD0l2ZpIhCXEQk26xYAX/9a5gOdOlSOO20cOHzNm3SXZmkmEJcRCRb/PRTGJx2111hfu9TTgnh3SF7rq0hRSnERUTibsGCsNv84YfDMe8TTwzze3funO7KJGIKcRGRuJo0Kewyf/ZZ2GEH6NuXC2Zew7zVbRir/N4uaCpSEZE42bQJ/v1vyM0NPe3XXw/nfM+fD489xryddNx7e6KeuIhIHKxbB088AcOGhalB99033L/oIth553RXJ2miEBcRyWRLlsCDD4bR5j/8AO3ahTA/80yoUaPCl0t2U4iLiGSi+fPDYLVHHgmD1U44IUwH2r07mJX70gULqqRCyQAKcRGRTDJxYhis9txzUK0a9O0LV19dqXO8v/wywvokoyjERUTSzDY5vPpquK75uHHhGPe118KVV0LjxukuTzKYRqeLiKTL2rXw8MM886dPw4VZ5s8PF2r5+mu4445KBXheXtjL/v774bFZuOXlRVK5ZIikeuJmtgNhru+9gTXATHf/IcrCRESy1pIl4cpqf/sb/PADP+9bG558En71q60erJaXF265uSHI3VNZsGSqckPczA4ErgeOBT4HfgRqAU3NbDUwAnjM3TdFXaiISOzNm7dlsNrq1dCjB1x3Hed8fh2T+vZNd3USQxX1xG8FHgAGuhf9v87Mdgf6Av2Ax6IpT0Qkerm54evYsRFsfMMGeO21cJrYG29A9epwzjnhAi2tWoV1vih/tHllNWmS0s1JBis3xMubE9zdFwH3pLogEZGs8M034VrmDz8MCxfCXnvBH/4AAwdGPlgtJyfSzUsGSfaY+K+AN9x9hZn9AegA3Orun0RanYhInGzcCGPGwIgRYbT5xo3h/O6//hVOPlkXZ5GUS/YUs5vc/VkzOxI4ARhK2M1+WGSViYjExfffw6hR8NBD4Uoru+8eLsxy8cVwwAHprk6yWLIhvjHx9STgAXd/2czyoilJRKTqVfoqZ5s2wXvvhWPdL70Ujn137x5ODevVC2rWLPGSSI+9E/22JfMkG+LfmNkIwij1O8xsR3SOuYhUgaoIPqjEVc7y82H06LDL/IsvYNddYfBgGDAAmjaNskSREpIN8TOBHsBQd19qZnsB10VXlohIBnGHDz4Ive7nn4eff4YjjwwnZvfuDbVqpbtC2U4l25tuA7zl7p8nHq8ElkVTkohI1ajwKmc//QT33gstW0LXruFUsYEDYcaMEOrnnKMAl7RKNsQfIAR3gVWJ50REYisvL3Syu3YNj93BNzl5J3wE/fvD3nvDVVeFa5mPGgXffhtGmrdsudXvqRnGJJWS3Z1uhS/24u6bzEyTp4hIlaiy4Lv//nCse9o02GmnEOQDB4Y5vFNEM4xJKiUbxPPM7Eq29L4vBeZFU5KISFGRBd+GDfDuu/z+syfoQnO47EZo3z4EeZ8+UK9eRG8skhrJ7k6/BDgC+AZYSDg/fEBURYmIRMYdJk8Oc3Tvuy+ccAJd8v/F4XvOhwkTwrIBA1Ia4JphTKKSVIi7+yJ3P9vdd3f3Pdy9b+KyqyIikUh58M2fD7fdBoccAp06wfDhcPjh8PzznH7E9wxt9hB07hzeJMVKPfbuCnHZdkmFuJk1NbN3zGxG4nGbxOVXRUQikZLgW7w4nBZ25JHhyml/+APstlvYXf799/DCC3D66fy8g0aYSzwle0z8IcJ54SMA3H2amf2DMMuZiEjmWLMmXLf8ySfDKWHr14fe9//9H/TtW+oUX1V5lTPNMCaplGyI13H3CVZ0N9OGCOoRESmhwuDbtCnsd3/iCXjuOVi+PMwaduWVcO650LZtJLvJt4ZmGJNUSjbE883sQMABzOwM4LvIqhIRKaTM4Js2LQT3U0+F6T532ilcQe3cc6FbN6hWrSrLFKlyyYb4ZcBIoLmZfQPMB86NrCoRkbJ8/TX84x889bdP4Zu2UL069OgBQ4fCKadAnTrprlCkyiQV4u4+DzjWzOoCO7j7imjLEhHZYqcNS+GR50Ov+/33wZ3VB9SF++6DM8+ERo3SXWJSNMOYpFqyo9MHm9nOwGrgbjP7xMyOj7Y0EdmuLVkCjz/O2F1O49VJe8JFF4XLnublwRdfcOH1zeHSS2MT4CJRSHZ3+gXufq+ZnQDsDpwPPAqMiawyEdn+LFwY5uZ+8cXQ4964ERo3Dpc+PffccH53wQC1d9JaqUhGSPra6YmvvwQedfepZuUP9TSzUcDJwCJ3b1XK8lzgZcLxdYAX3P3mJOsRkQyQkrm+P/sshPaLL8LEieG55s3ht7+FXr1CcO+Q7MUlRbYvyYb4ZDMbA+wP/N7M6gGbKnjNaGA48Hg563zg7icnWYOIZINNm2DSpC3BPXt2eL5z53Au92mnhRAXkQolG+IXAu2Aee6+2swaEnapl8ndx5lZzraVJyJZYf36sHv8pZfC7ZtvwulfublwxRVw6qmwzz5pLlIkfqzQDKMlF5rluPuCcpYb0NjdF5b1euDVcnanP0+YUOVb4Fp3n1nGdgaQmHDFGljHDnd0KLPmuJmVP4sWjVqku4yUyab2ZFNbIJr2zB42AoBm1wwssWzHnzdx+MxldJuylCOnLaP+6o2srWF81LI+Y9s14IM29Vled+tnNM6mn082tQXUnihMHjh5srt3KrHA3cu8Ac8SgvY8oCVhUNt+QHfgFmA8cFw5r88BZpSxbGdgp8T9XwKfl1dLwa3OfnU8m3Qc0THdJaRUNrUnm9riHk17unZ1b9Kk0BOLF7uPHu3eq5d77drhcue77OJ+3nnuL77ovmpVyt47m34+2dQWd7UnCsAkLyUTy/032N1/ZWaHAOcAFwB7EU4zmwW8Btzm7mu35r8Kd19e6P5rZna/mTVy9/yt2Z6IpMeXXxJmBCs+ovyCC8Lx7aOPhho10l2mSFaqcF+Wu38K3JjqNzazPYEf3N3N7FDCOeuLU/0+ItujlIwaL8v69fDxx/Dmm4yY/AbNmRSOaxeMKD/ttKKngolIZLb+gFQFzOwpIBdoZGYLgSFADQB3fxA4AxhkZhuANcDZiV0GIpJp5s2DN98Mt3ffJW/F1fyJWymYyNBw+AyG1IS8zuktVWR7ElmIu3ufCpYPJ5yCJiKZZsUKeO+9LcE9d254vkkT6NOHvONbk3fMUnJ7NSi4CqqIpEFkIS4iMbJpE/zvf1tCe/x42LAhTCbSrRsMHgzHHw9Nm2o3uUgGKTfEzazcc7nc/ZPUliMiqbJgQQUrfPstvPVWCO233oL8xJjSdu3gmmvghBPgiCNgxx3L3UyFc32LSGQq6okPS3ytBXQCphIuwdoG+C9wZHSlici2+PLLoo9rrt8UwnrMmBDc06eHBbvvHqbyPOEEOO442GOPSr1PmXN9i0jkKjrFrBuAmT0NDHD36YnHrYBroy9PJIh0tHU2+/TTENhjxvDuO1Ng/fFQsyYceST8+c8huNu00bXJRWIq2WPizQsCHMDdZ5hZu2hKEpGt4k7e5fn86f7dNj9lLQ8BDmFIwxo0O+oT+lw9KvxHVLdu2soUkdRJNsRnmdnDwBOAA+cSLvgiIumyYQNMmQIffADjxsGHH5KXn08esKTG7jRc/wM+8qEwIK3J5XQaOZo+J52U0hK0Z0QkvZIN8fOBQcDgxONxwAORVCQipVuzBv773xDaH3wAH30EK1eGZQccACedBEcdBUcdxekXHxz+Si++OK0li0i0kgpxd19rZvcBbxN64rPdfX2klYkUU+Fo62yzdCn85z9bQnvixHC1NDNo1QrOO29zaNO4cdHXmkaNi2wPkgrxxIxjjwELCKPT9zWzX7v7uMgqEymm+GjrrPPtt1sC+4MPwuhxd6hePcy1/ZvfhMDu0gV22aXCzWnUuEj2S3Z3+jDgeHefDWBmTYGngI5RFSaSrXJzAXfGPvxF0dAuuCpa3bpw+OGQlxdC+7DDwkVXRESKSTbEaxQEOIC7zzEzTUu0DQpOmaJv1bxPXAcg5eXBn/605XHBxcKGDAnLYmPJEpg8GSZN4uaZE2m5/CNo+n1Y1rBhCOtLLw1f27XTrF8ikpRkQ3ySmT0C/D3x+BxgcjQliWyRlxduubnE5xrdy5bBJ5/ApElbbvPmbV68f+2D+aTBMRx3c+J4dvPmKT9PO67/tIlI5SQb4oOAy4ArCcfExwH3R1WUxEvce/vbZOXKcM3xwoE9Z86W5Tk5YVrOgQPD1w4d6NerAQDHDUxLxSKSRZIdnb4OuCtxE0mLqhhtXe4/JKtXw9SpRQN71qwtuwf22ScE9Xnnha8dO0KjRqW+z3Y30l5EIlHRBCjPuPuZZjadcGpZEe7eJrLKtgMLFkDpH/Gpf59sUKWjrdet45AFq+CBB7YE9syZsHFjWL7HHmHE+JlnbgnsPfdMevNZP9JeRKpERT3xgou7nBx1IdujL7+smhBXYJTDHRYtCgE9YwbXzJlB0xWTod50Hl+/Hrg0DDzr3Bl69gyB3akT7L23puQUkbSraAKU78ysGvCIux9bRTVJDMWit794cQjrRGBv/rp48eZVulbfhc936gCXXc31Pz3LHTe8C/vtl5LAzpqR9iKSMSo8Ju7uG81stZnVd/dlVVFUNiv+QT554CRsYOo/yKs6MKLu7Vdq0NyyZWH2rsJBPXMmfP/9lnXq1QtXPTvttPC1ZUto2ZKeZ+/Jgi+NBX+Gd0a+ndID8bEcaS8iGS3Z0elrgelm9hawquBJd78ykqqyWPEP8o4jOjFpwKTI3ycrA2PVqhDWxXvWCxduWadOHTjkkDBfdiKoadUqDEIrrXdtOvwgIvGRbIj/O3ET2az03v4k8r5NYW/fPezunjsXPv+8aGDPn79lvR13hBYtoGvXIj1rcnIycq5sXddcRFIh2VPMHjOz2sB+ha/cJtumqj7Io3qf0nr7nUZ2Iq+yexY2bQq957lzS78tK3QUp3p1aNYMDj0Uzj9/S8/6gAPCsm1oS+T/kBSi65qLSCokOwHKKcBQoCawv5m1A252954R1pb1cnJgZRW9T9qtWxd6zqWF9Pz5YXmB6tVD0QcdFK4hfuCB4XbQQeFWs2bKy0vZPyQiIlUo2a5LHnAoMBbA3aeY2f4R1SQx1KQJsGwZzb5aDc8+WzSkv/gi9LQLH5jfaacQzIccAqecsiWoDzwQ9t13m3rVIiLbi2Q/KTe4+zIrOhAoG4dKSWncYfly+OabUm8jJn/DnmsXQIPFPAnAmeF1u+8eQrlr16IhfdBBsNtuGXueddSHObbLy9OKSCSSDfEZZtYXqGZmBxOuoT4+urKyX8EHeaeRVfM+Zdq4MZx6VUZA8803oRe9alXJ1+66KzRuzNIajZldryM9rzqQ6+aN4M5BL4Rj1PXqRdGkyGXE4QcRkSQkG+JXADcC64B/AG8Ct0ZVlKTAhg3w00+Qnw/ffReCuLSA/v77MLCssOrVwxXJGjeG1q3hxBPD/cK3vfeG2rUBuD43vKznb+G9kc9A27ZV21YRke1URddOrwVcAhwETAcOd/cNVVGYFLJuXTjNKj8/fC1+v7THP/1U+rbq198SxK1alQznxo3Dru5KnJal3cMiIulRUU/8MWA98AFwItACuCrimrKPe5gBa+XKsFu64Ovy5ZwwYQms+1vpYVxwv7Rd2QXq1g3X9m7UKHzdf/+ijxs2hL322hLQdetWXbtjSP+QiEicVBTih7h7awAzewSYEH1JGSA/H5YuLRm6pX1NZp1yQvg2IAwxAHbZpWjwtmq15XHhUC54vOuuUKtWFXxDREQkE1UU4usL7rj7BsvQ0cQpd845MGZMxevtsANUq7bltueeoSfcoAFMnhye22WXELjVqoWTkLt2hTVr4L77oFo1ZiydQ6u924Xj0NdeG063mj0bBg4M/wAUvgboH/4A3brBlCnQt2/Jev7v/+CII2D8eLjhhpLL77kH2rWDt9+GW0sZ0jBiRLiQyiuvwLBhJZf//e/h9K9//jNM0Vncc8+Fr6NHh1txr70WLoN6//3wzDMllxd0g4cOhVdfLbqsdm14/fVw/5Zb4J13ii5v2BCefz7c//3v4aOPii7fZx944olw/6qrwvewsKZNYWRilOGAATBnDiO+mw3/yA3PtWsXvn8A555b9NKuEM5nv/32cL937yKTqgBwzDFw003h/oknht+Bwk4+Ofz8Ycuk5oWdeSZcemnYo/PLX5Zc3r9/uOXnwxlnlFw+aFD4+vXX0K9fyeXXXFP0d6+4P/wBjj02fN+uuqrk8jT87hX5+Tz3XPg7i+vvXqfE/cTvXhEx/N0r8rMZNAjOOivzf/diqqIQb2tmyxP3DaideGyAu/vOkVaXLldfDV99VTSgq1ULv2gDB4bw/tWvSr6uog/S44/f8sv81FMArFu9A9SoEWlzREQkO5nHbGaMuk3q+qovyzlGHDOdRkYzAUq6ZFN7sqktoPZksmxqC6g9UTCzye7eqfjzmTczRAbIzS19r5KIiEgmUYiLiIjElEJcREQkpiILcTMbZWaLzGxGBet1NrONZlbKSDAREREpS5Q98dFAj/JWMLNqwB2Ey7hmlAUL0l2BiIhI+SILcXcfByypYLUrgOeBRVHVsbUKn54tIiKSidI2abOZNQZOA7oDnStYdwAwAMAaGJ1Glhhln1KzvxsBdIz8fQBm5c+qkvepKtnUnmxqC6g9mSyb2gJqT5Vy98huQA4wo4xlzwK/SNwfDZyRzDbr7FfHozJkiHu40HnR25Ahkb2ldxzRMbqNp0E2tSeb2uKu9mSybGqLu9oTBWCSl5KJ6Ryd3gl42swWAGcA95tZrzTWQ15eiO2uXcPjghjPy0tnVSIiIqVL2+50d9+/4L6ZjQZedfeX0lWPiIhI3EQW4mb2FJALNDKzhcAQoAaAuz8Y1fumSpMm6a5ARESkfJGFuLv3qcS6/aOqY2vl5KS7AhERkfLpim0iIiIxlbZj4pmsYGphERGRTKaeuIiISEwpxEVERGJKIS4iIhJTCnEREZGYUoiLiIjElEJcREQkphTiIiIiMaUQFxERiSmFuIiISEwpxEVERGJKIS4iIhJTCnEREZGYUoiLiIjElEJcREQkphTiIiIiMaUQFxERiSmFuIiISEwpxEVERGJKIS4iIhJTCnEREZGYUoiLiIjElEJcREQkphTiIiIiMaUQFxERiSmFuIiISEwpxEVERGJKIS4iIhJTCnEREZGYUoiLiIjElEJcREQkphTiIiIiMaUQFxERiSmFuIiISExFFuJmNsrMFpnZjDKWn2pm08xsiplNMrMjo6pFREQkG0XZEx8N9Chn+TtAW3dvB1wAPBxhLSIiIlknshB393HAknKWr3R3TzysC3hZ64qIiEhJtiVHI9i4WQ7wqru3KmP5acDtwO7ASe7+URnrDQAGAFgD69jhjg7RFJwGs/Jn0aJRi3SXkTLZ1J5saguoPZksm9oCak8UJg+cPNndO5VY4O6R3YAcYEYS6x0NvJ3MNuvsV8ezSccRHdNdQkplU3uyqS3uak8my6a2uKs9UQAmeSmZmBGj0z3sej/QzBqluxYREZG4SFuIm9lBZmaJ+x2AmsDidNUjIiISN9Wj2rCZPQXkAo3MbCEwBKgB4O4PAr2B88xsPbAGOCuxy0BERESSEFmIu3ufCpbfAdwR1fuLiIhku4w4Ji4iIiKVpxAXERGJKYW4iIhITCnERUREYkohLiIiElMKcRERkZhSiIuIiMSUQlxERCSmFOIiIiIxpRAXERGJKYW4iIhITCnERUREYkohLiIiElMKcRERkZhSiIuIiMSUQlxERCSmFOIiIiIxpRAXERGJKYW4iIhITCnERUREYkohLiIiElMKcRERkZhSiIuIiMSUQlxERCSmFOIiIiIxpRAXERGJKYW4iIhITCnERUREYkohLiIiElMKcRERkZhSiIuIiMSUQlxERCSmFOIiIiIxpRAXERGJKYW4iIhITEUW4mY2yswWmdmMMpafY2bTErfxZtY2qlpERESyUZQ98dFAj3KWzwe6unsb4BZgZIS1iIiIZJ3qUW3Y3ceZWU45y8cXevgxsE9UtYiIiGQjc/foNh5C/FV3b1XBetcCzd39ojKWDwAGAFgD69jhjg6pLjVtZuXPokWjFukuI2WyqT3Z1BZQezJZNrUF1J4oTB44ebK7dyqxwN0juwE5wIwK1ukGzAIaJrPNOvvV8WzScUTHdJeQUtnUnmxqi7vak8myqS3uak8UgEleSiZGtjs9GWbWBngYONHdF6ezFhERkbhJ2ylmZrYf8ALQz93npKsOERGRuIqsJ25mTwG5QCMzWwgMAWoAuPuDwB+BhsD9ZgawwUvb3y8iIiKlinJ0ep8Kll8ElDqQTURERCqmK7aJiIjElEJcREQkphTiIiIiMaUQFxERiSmFuIiISEwpxEVERGJKIS4iIhJTCnEREZGYUoiLiIjElEJcREQkphTiIiIiMaUQFxERiSmFuIiISEwpxEVERGJKIS4iIhJTCnEREZGYUoiLiIjElEJcREQkphTiIiIiMaUQFxERiSmFuIiISEwpxEVERGJKIS4iIhJTCnEREZGYUoiLiIjElEJcREQkphTiIiIiMaUQFxERiSmFuIiISEwpxEVERGLK3D3dNVSKma0AZqe7jhRqBOSnu4gUyqb2ZFNbQO3JZNnUFlB7otDE3Xcr/mT1dFSyjWa7e6d0F5EqZjZJ7clM2dQWUHsyWTa1BdSeqqTd6SIiIjGlEBcREYmpOIb4yHQXkGJqT+bKpraA2pPJsqktoPZUmdgNbBMREZEgjj1xERERIYND3Mx6mNlsM/vCzH5XynIzs78mlk8zsw7pqDMZSbTlnEQbppnZeDNrm446k1VRewqt19nMNprZGVVZX2Ul0x4zyzWzKWY208zer+oaKyOJ37f6ZvaKmU1NtOf8dNSZDDMbZWaLzGxGGcvj9DlQUVvi9jlQbnsKrReXz4EK25ORnwPunnE3oBowFzgAqAlMBQ4pts4vgdcBA34B/DfddW9DW44AdkncPzFT25Jsewqt9y7wGnBGuuvexp9PA+BTYL/E493TXfc2tucG4I7E/d2AJUDNdNdeRnuOBjoAM8pYHovPgSTbEpvPgWTak1gnFp8DSf58MvJzIFN74ocCX7j7PHf/GXgaOLXYOqcCj3vwMdDAzPaq6kKTUGFb3H28u/+UePgxsE8V11gZyfxsAK4AngcWVWVxWyGZ9vQFXnD3rwDcPZPblEx7HKhnZgbsRAjxDVVbZnLcfRyhvrLE5XOgwrbE7HMgmZ8NxOdzIJn2ZOTnQKaGeGPg60KPFyaeq+w6maCydV5I6FlkqgrbY2aNgdOAB6uwrq2VzM+nKbCLmY01s8lmdl6VVVd5ybRnONAC+BaYDgx2901VU17KxeVzoLIy/XOgQjH7HEhGRn4OZOoV26yU54oPo09mnUyQdJ1m1o3wx3tkpBVtm2Tacw9wvbtvDJ29jJZMe6oDHYFjgNrAR2b2sbvPibq4rZBMe04ApgDdgQOBt8zsA3dfHnFtUYjL50DSYvI5kIx7iM/nQDIy8nMgU0N8IbBvocf7EHoNlV0nEyRVp5m1AR4GTnT3xVVU29ZIpj2dgKcTf7iNgF+a2QZ3f6lKKqycZH/X8t19FbDKzMYBbYFMDPFk2nM+8GcPB/a+MLP5QHNgQtWUmFJx+RxISow+B5IRp8+BZGTk50Cm7k6fCBxsZvubWU3gbOBfxdb5F3BeYnTqL4Bl7v5dVReahArbYmb7AS8A/dL9X10SKmyPu+/v7jnungM8B1yawX+4yfyuvQwcZWbVzawOcBgwq4rrTFYy7fmK0JvAzPYAmgHzqrTK1InL50CFYvY5UKGYfQ4kIyM/BzKyJ+7uG8zscuBNwujGUe4+08wuSSx/kDDa8ZfAF8BqQu8i4yTZlj8CDYH7E/+1bvAMvdh+ku2JjWTa4+6zzOwNYBqwCXjY3cs9rSZdkvz53AKMNrPphN3R17t7umdoKpWZPQXkAo3MbCEwBKgB8focgKTaEpvPAUiqPbFSUXsy9XNAV2wTERGJqUzdnS4iIiIVUIiLiIjElEJcREQkphTiIiIiMaUQFxERiSmFuEgGMrPTzMzNrHkKt5lrZq8m7ve0xAxnZtbLzA7Ziu2NNbNKnQJlZveY2dGVfa9i21iZ+Lpb4pQfke2WQlwkM/UBPiRcrCXl3P1f7v7nxMNeQKVDvLLMbFfgF4mJJoovq1bZ7bn7j8B3ZtYlFfWJxJFCXCTDmNlOQBfC9bPPLvR8rpm9b2bPmNkcM/uzhTmoJ5jZdDM7MLHeaDN70Mw+SKx3cinv0d/MhpvZEUBP4M7EPMkHFu5hm1kjM1uQuF/bzJ62MN/1PwnXjy7Y3vFm9pGZfWJmzybaUNwZwBuFXrPAzP5oZh8CvzKzi81sooV5zp9PXBWLxNXnPkosu6XYNl8Czqn0N1kkSyjERTJPL+CNxKU3l5hZh0LL2gKDgdZAP6Cpux9KuN72FYXWywG6AicBD5pZrdLeyN3HEy5dep27t3P3ueXUNQhY7e5tgNsIk0FgZo2APwDHunsHYBJwdSmv7wJMLvbcWnc/0t2fJkzz2Nnd2xIuZ3lhYp17gQfcvTPwfbHXTwKOKqdmkaymEBfJPH0I84CT+Nqn0LKJ7v6du68D5gJjEs9PJwR3gWfcfZO7f064Lnoqjq0fDTwB4O7TCJefBPgFYXf8f8xsCvBroEkpr98L+LHYc/8sdL9VYu/BdELvumXi+S7AU4n7fy/2+kXA3pVuiUiWyMhrp4tsr8ysIWGK0FZm5oTrn7uZ/TaxyrpCq28q9HgTRf+ei19PuTLXV97Aln/wi/fgS9uOAW+5e59SlhW2ppTtrSp0fzTQy92nmll/wnWsy3vfgvrWVPC+IllLPXGRzHIG8Li7N0nMALUvMJ/Kzy39KzPbIXGc/ABgdjnrrgDqFXq8gMSu8kQ9BcaROP5sZq2ANonnPwa6mNlBiWV1zKxpKe8zCzionDrqEQaq1aDoce7/sGVsQPHj302BtE9CIZIuCnGRzNIHeLHYc88DfSu5ndnA+8DrwCXuvracdZ8GrjOz/yVCfygwyMzGE+aBLvAAsJOZTQN+S2L+8cQo8f7AU4llH1P67vt/U7R3XdxNwH+Bt4DPCj0/GLjMzCYC9Yu9pltiuyLbJc1iJpJlzGw08Kq7P5fuWopLjEQ/2d2Xpmh744BT3f2nVGxPJG7UExeRqnQNsF8qNmRmuwF3KcBle6aeuIiISEypJy4iIhJTCnEREZGYUoiLiIjElEJcREQkphTiIiIiMaUQFxERian/B+OJzB5hf/LLAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "periodArrLg = np.loadtxt(\"periodTableLarge.csv\", delimiter = ',')\n", "print(\" \")\n", "print(\"Shape of periodArrLg =\",periodArrLg.shape)\n", "#\n", "nMeas = periodArrLg.shape[1]\n", "print(\"nTheta =\",nTheta,\"and nMeas =\",nMeas)\n", "#\n", "meanPeriods = np.mean(periodArrLg, axis = 1)\n", "errPeriods = np.std(periodArrLg, axis = 1, ddof = 1)/np.sqrt(nMeas)\n", "print(\" \")\n", "print(\"Number\\t Amplitude (rad)\\t Mean period (s)\")\n", "for n in range(0, nTheta):\n", " print(f\"{n:d}\\t {thetaArr[n]:.2f}\\t\\t\\t {meanPeriods[n]:.2f} +- {errPeriods[n]:.2f}\")#\n", "g = 9.81\n", "ySHM = 2*np.pi*np.sqrt(L/g)*np.ones(nTheta)\n", "chi2SHM = (ySHM - meanPeriods)**2/errPeriods**2\n", "chi2sumSHM = np.sum(chi2SHM)\n", "print(\" \")\n", "print(f\"Chi2 for simple approximation = {chi2sumSHM:.2f}, chi2/NDF = {chi2sumSHM/(nTheta - 1):.2f}\")\n", "#\n", "chi2full = (periodLima(L, thetaArr) - meanPeriods)**2/errPeriods**2 \n", "chi2sumFull = np.sum(chi2full)\n", "print(\" \")\n", "print(f\"Chi2 for Lima/Arun approximation = {chi2sumFull:.2f}, chi2/NDF = {chi2sumFull/(nTheta - 1):.2f}\")\n", "#\n", "print(\" \")\n", "plt.figure(figsize = (8, 6))\n", "plt.title('Pendulum data', fontsize = 14)\n", "plt.xlabel('Amplitude (rad)')\n", "plt.ylabel('Period (secs)')\n", "plt.errorbar(thetaArr, meanPeriods, yerr = errPeriods, \n", " linestyle = '', marker = '+', color = 'b', label = \"data\")\n", "plt.plot(thetaArr, periodLima(L, thetaArr), color = 'r', label = \"Good approx.\")\n", "plt.plot(thetaArr, ySHM, color = 'r', linestyle = '--', label = \"Simple approx\")\n", "plt.xlim(0.0, 1.1*thetaArr[nTheta - 1])\n", "plt.ylim(0.8*np.min(meanPeriods), 1.1*np.max(meanPeriods))\n", "plt.legend()\n", "plt.grid(color = 'g')\n", "plt.savefig(\"periodDataLargeNew.png\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The SHM formula results in a $\\chi^2$/NDF which is outside the range 0.2 to 5.0, so this formula does not describe the data. The $\\chi^2$/NDF value for the Lima/Arum formula is inside the acceptable range, so this does describe the data.\n", "\n", "These results demonstrate the importance of precision in physics experiments!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Monte Carlo revisited\n", "\n", "The \"experimental data\" used above wasn't actually a set of measurements of the motion of a pendulum. (It is difficult to precisely measure the behaviour of a pendulum for large amplitudes!) It was generated using the code in the cell below." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "Debug = False\n", "#\n", "# Range of amplitudes (angles) considered \n", "thetaBot = 0.1 # rad\n", "thetaTop = np.pi/2 # rad\n", "nTheta = 16\n", "thetaArr = np.linspace(thetaBot, thetaTop, nTheta) # rad\n", "#\n", "# Pendulum length \n", "L = 0.5 # m\n", "#\n", "# Assumed measurement precision\n", "periodSigma = 0.5 # seconds\n", "#\n", "nMeas = 10\n", "periodTabSm = np.zeros((nTheta, nMeas))\n", "#\n", "periodLimaValues = np.zeros(nTheta)\n", "#\n", "for n in range(0, nTheta):\n", " periodLimaValues[n] = periodLima(L, thetaArr[n]) # seconds\n", " periodTabSm[n, 0:nMeas] = np.random.normal(periodLimaValues[n], scale = periodSigma, size = nMeas) # seconds\n", "if Debug: print(periodTab)\n", "#\n", "np.savetxt(\"thetaTable.csv\", thetaArr, delimiter = ',', )\n", "np.savetxt(\"periodTableSmall.csv\", periodTabSm, delimiter = ',')\n", "#\n", "nMeas = 1000\n", "periodTabLg = np.zeros((nTheta, nMeas))\n", "#\n", "for n in range(0, nTheta):\n", " periodTabLg[n, 0:nMeas] = np.random.normal(periodLimaValues[n], periodSigma, nMeas) # seconds\n", "if Debug: print(periodTab)\n", "#\n", "np.savetxt(\"periodTableLarge.csv\", periodTabLg, delimiter = ',')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Summary\n", "\n", "I hope that revising calculating means and their standard errors, fitting techniques (last week) and the calculation and use of $\\chi^2$ values is useful for your practicals!" ] } ], "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 }