{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# MC for Phys106 Week 1 - Introduction to error analysis\n", "\n", "Create raw data for $g$ measurement, twenty sets of the periods obtained by measuring 20 oscillations for each of 10 pendulum lengths." ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " \n", "lengths (m)\n", " [0.1 0.2 0.2 0.3 0.3 0.4 0.4 0.5 0.5 0.6]\n", "Tideal (s)\n", " [0.8 0.9 1. 1.1 1.2 1.3 1.3 1.4 1.5 1.6]\n", "sigma = 0.025 s\n", " \n", "Length time measurements:\n", " \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\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.020\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" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgAAAAGECAYAAABAsZipAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XuUZWdd5vHv0wkQC6KA3SCTpLsYDCOgEDJFCCugiUBMlOHiMEgsMSJYi5uiXJZgIyjYLEdcoAzXEkJAD2FUCBe5ZhSMXALpYExIIhBiutMmkoYECBQYkvzmj72LnK7U5VSnzqmqs7+ftc7aZ7/73ef8aq/q3k+9+5aqQpIkdcuW9S5AkiSNngFAkqQOMgBIktRBBgBJkjrIACBJUgcZACRJ6iADgKShS/LEJF5zLG0gBgCpI5KcmaSSvGWRZX/SLvu7Be2HJ3lFkkuTfDfJV5N8IslpSYb6/0eSK5O8YJjfIXWZAUDqlquAX0py5/mGJIcCTwH29ndMclfgM8CvA68CpoCHA28Hfh/YPqKaJQ2BAUDqlouALwNP6mv7BeB7wCcW9H0lcG/goVX1tqq6pKq+XFVvA44F/mOpL0nyq0n2JJlrRxXuuWD5fZK8L8l/JPlOks8neUzf8k8AO4BXtSMT1bb/aJKzkuxrRyQuSfLUg9wWUqcZAKTueSvNX/Xzfh14G/CDY/Tt8P6TgV5V7Vv4AVX1var63mIfnuShwJnALHAM8AHg5Qu63QX4MPBo4EHAu4H3JPmJdvkvAvva9e7VvgAOAz4PPAZ4APDnwJuTPHKAn1tSn/gsAKkbkpwJbKUZ7r8aeCBwA7AHOJpmZ7u1qh6T5B7AV4HnVdVrVvk97wS2VdWj+9reAjytqrLMeucBf1dVf9TOXwm8rqr+dIXvexfw7ap6+mrqlLrOEQCpY6rqeuBsmr/8Twc+UVV7F3Rbckc9gPvRnDvQ74D5JHduTzy8NMn1Sb5Nc47BsucVJDkkyc4kFyX5erveL660nqTbOnS9C5C0Ls6gOZnv28BLF1m+H7ieZme+WoOEhz8FTgFeQHNOwhzwDuCOK6z3AuD5wHOBi2nqfyVwj4OoU+o0RwCkbvp74EaaQwLvXbiwqm4B/i8wneTIhcuTHJbksCU++1Lg+AVtC+cfDryjqt5dVRfRHO+/z4I+NwKHLLLeB6rqL6vqQuArwH2XqEPSMgwAUgdVc/LPA4F7V9V/LtHt92guDfxskqcmeUCSH0/yFOAC4MeWWO+1wKOSvDjJ0Ul+A3jCgj5fAp6Q5NgkPwX8Fc0Jfv2uBB6R5IgkW/vWe2SSh7cnDL6O5koFSatkAJA6qqpuqKpvLbP8epq/3M8Efpdmp/9p4GnAK1hw34C+9c5r+zyT5rLDXwT+YEG35wHXAv9EczXAee37fi8FjqL5K39/2/ZHwOfadc4FvgP0VvhRJS3CqwAkSeogRwAkSeogA4AkSR1kAJAkqYMMAJIkdZABQJKkDhrrOwFu3bq1Jicn17sMSZJG5oILLvhaVW1bqd9YB4DJyUl279693mVIkjQySfYM0s9DAJIkdZABQJKkDjIASJLUQQYASZI6yAAgSVIHGQAkSeogA4AkSR1kAJAkqYMMAJIkdZABQJKkUev1YHIStmxppr3eyEsY61sBS5K04fR6MDMDc3PN/J49zTzA9PTIynAEQJKkUdq589ad/7y5uaZ9hAwAkiSN0t69q2sfEgOAJEmjtH376tqHxAAgSdIo7doFExMHtk1MNO0jZACQJGmUpqdhdhZ27ICkmc7OjvQEQPAqAEmSRm96euQ7/IUcAZAkqYMMAJIkdZABQJKkDjIASJLUQQYASZI6yAAgSVIHGQAkSeogA4AkSR1kAJAkqYMMAJIkdZABQJKkDjIASJLUQQYASZI6yAAgSVIHGQAkSeogA4AkSR1kAJAkqYMMAJIkddBIA0CSo5J8PMllSS5J8txF+kwnuah9fTrJg/qWXZnk4iQXJtk9ytolSRonh474+24Cnl9Vn09yOHBBknOq6tK+Pv8G/ExVXZ/kVGAWeGjf8pOq6msjrFmSpLEz0gBQVdcA17Tvb0hyGXAEcGlfn0/3rXIecOQoa5QkqQvW7RyAJJPAg4HPLtPtacCH++YL+FiSC5LMLPG5M0l2J9m9f//+tSpXkqSxMupDAAAkuQvwbuC3q+pbS/Q5iSYAPLyv+YSqujrJPYBzkvxrVZ3bv15VzdIcNmBqaqqG8gNIkrTJjXwEIMkdaHb+vap6zxJ9Hgi8BXhcVX19vr2qrm6n1wJnA8cNv2JJksbPqK8CCPBW4LKqevUSfbYD7wGeUlVf6mu/c3viIEnuDJwMfGH4VUuSNH5GfQjgBOApwMVJLmzbfg/YDlBVbwJeCvwo8IYmL3BTVU0B9wTObtsOBd5ZVR8ZbfmSJI2HUV8F8EkgK/R5OvD0RdqvAB502zUkSdJqeSdASZI6yAAgSdrcej2YnIQtW5ppr7feFW0K63IZoCRJa6LXg5kZmJtr5vfsaeYBpqfXr65NwBEASdLmtXPnrTv/eXNzTbuWZQCQJG1ee/eurl0/YACQJG1e27evrl0/YACQJG1eu3bBxMSBbRMTTbuWZQCQJG1e09MwOws7dkDSTGdnPQFwAF4FIEna3Kan3eEfBEcAJEnqIAOAJEkdZACQJKmDDACSJHWQAUCSpA4yAEiS1EEGAEmSOsgAIElSBxkAJEnqIAOAJEkdZACQJKmDDACSJHWQAUCSpA4yAEiS1EEGAEmSOsgAIElSBxkAJEnqIAOAJEkdZACQJKmDDACSJHWQAUCSdKteDyYnYcuWZtrrrXdFGpJD17sASdIG0evBzAzMzTXze/Y08wDT0+tXl4bCEQBJUmPnzlt3/vPm5pp2jZ2RBoAkRyX5eJLLklyS5LmL9EmS1ya5PMlFSY7tW3Z6ki+3r9NHWbskjb29e1fXrk1t1CMANwHPr6r7AccDz05y/wV9TgWObl8zwBsBktwdeBnwUOA44GVJ7jaqwiVp7G3fvrp2bWojDQBVdU1Vfb59fwNwGXDEgm6PA95RjfOAuya5F/BzwDlVdV1VXQ+cA5wywvIlabzt2gUTEwe2TUw07Ro763YOQJJJ4MHAZxcsOgK4qm9+X9u2VPvCz51JsjvJ7v37969lyZI03qanYXYWduyApJnOznoC4Jhal6sAktwFeDfw21X1rYWLF1mllmk/sKFqFpgFmJqaus1ySdIypqfd4XfEyEcAktyBZuffq6r3LNJlH3BU3/yRwNXLtEuSpFUa9VUAAd4KXFZVr16i2/uBX22vBjge+GZVXQN8FDg5yd3ak/9ObtskSdIqjfoQwAnAU4CLk1zYtv0esB2gqt4EfAj4eeByYA54arvsuiSvAM5v13t5VV03wtolSRobIw0AVfVJFj+W39+ngGcvsewM4IwhlCZJUqd4J0BJkjrIACBJUgcZACRJ6iADgCRJHWQAkCSpgwwAkiR1kAFAkqQOMgBIktRBBgBJkjrIACBJUgcZACRJ6iADgCRJHWQAkCSpgwwAkiR1kAFAkqQOMgBIktRBBgBJkjrIACBJUgcZACRJ6iADgCRJHWQAkCSpgwwAkjQsvR5MTsKWLc2011vviqQfOHS9C5CksdTrwcwMzM0183v2NPMA09PrV5fUcgRAkoZh585bd/7z5uaadmkDMABI0jDs3bu6dmnEDACSNAzbt6+uXRoxA4AkDcOuXTAxcWDbxETTLm0ABgBJGobpaZidhR07IGmms7OeAKgNw6sAJGlYpqfd4WvDcgRAkqQOMgBIktRBBgBJkjpopOcAJDkDeAxwbVX95CLLXwjMHzA7FLgfsK2qrktyJXADcDNwU1VNjaZqSZLGz6hHAM4ETllqYVW9qqqOqapjgBcD/1hV1/V1Oald7s5fkqTbYaQBoKrOBa5bsWPjNOCsIZYjSVJnbchzAJJM0IwUvLuvuYCPJbkgycwy684k2Z1k9/79+4ddqiRJm9KGDADA/wA+tWD4/4SqOhY4FXh2kp9ebMWqmq2qqaqa2rZt2yhqlSRp09moAeDJLBj+r6qr2+m1wNnAcetQlyRJY2HDBYAkPwL8DPC+vrY7Jzl8/j1wMvCF9alQkqTNb9SXAZ4FnAhsTbIPeBlwB4CqelPb7QnAx6rqO32r3hM4Owk0Nb+zqj4yqrolSRo3Iw0AVXXaAH3OpLlcsL/tCuBBw6lKkqTu2XCHACRJ0vAZACRJ6iADgCRJHWQAkCSpgwwAkiR1kAFAkqQOMgBIktRBBgBJkjrooAJAe2veQ9a6GEmSNBoDBYAkW5L8cpIPJrkW+FfgmiSXJHlVkqOHW6YkSVpLg44AfBy4D/Bi4Meq6qiqugfwCOA84I+T/MqQapQkSWts0GcBPKqqvr+wsaquA94NvDvJHda0MkmSNDQDjQAstvM/mD6SJGljWDEAJHl0kr9Ickw7PzP8siRJ0jANcgjgWcBTgZckuTtwzHBLkiRJwzbIIYD9VfWNqnoBcDLwkCHXJEmShmyQAPDB+TdV9SLgHcMrR5IkjcKKAaCq3jf/Psn/As5s378kyXuSHDu88iSp1evB5CRs2dJMe731rkja1FZ7J8Dfr6obkjwc+Dng7cAb174sSerT68HMDOzZA1XNdGbGECDdDqsNADe3018A3tiODtxxbUuSpAV27oS5uQPb5uaadkkHZbUB4N+TvBl4EvChJHc6iM+QpNXZu3d17ZJWtNqd95OAjwKnVNU3gLsDL1zzqiSp3/btq2uXtKJBHwYUgKqaq6r3VNWX2/lrqupj/X0kac3t2gUTEwe2TUw07ZIOysAPA0rym0kOiNtJ7pjkZ5O8HTh97cuTJGB6GmZnYccOSJrp7GzTLumgpKpW7pQcBvw6MA3cG/gGcBhwCPAx4PVVdeEQ6zwoU1NTtXv37vUuQ5KkkUlyQVVNrdRvoKcBVtX3gDcAb2if+rcV+G57HoAkSdpkBn0c8A+0T/27Zgi1SJKkEfESPkmSOsgAIElSBxkAJEnqoBUDQJJHJ/mLJMe08zPDL0uSJA3TICcBPgt4KvCSJHcHjhluSZIkadgGOQSwv6q+UVUvAE4GHnKwX5bkjCTXJvnCEstPTPLNJBe2r5f2LTslyReTXJ7kRQdbgyRJGiwAfHD+TVW9CHjH7fi+M4FTVujzT1V1TPt6OUCSQ4DXA6cC9wdOS3L/21GHJEmdtuIhgKp6X5Ln9bctnK+qVw/yZVV1bpLJ1RTYOg64vKquaL//XcDjgEsP4rMkSeq8Qa8COLx9TQHPBI5oX8+g+Yt8LT0syb8k+XCSB7RtRwBX9fXZ17ZJkqSDMOitgP8QIMnHgGOr6oZ2/g+Av1nDej4P7Kiqbyf5eeC9wNHAYk8aXPQhBu1VCjMA231UqCRJi1rtfQC2Azf2zd8ITK5VMVX1rar6dvv+Q8Adkmyl+Yv/qL6uRwJXL/EZs1U1VVVT27ZtW6vSJEkaK6t9FsBfAp9LcnY7/3jg7WtVTJIfA75aVZXkOJqA8nWapw8eneTewL8DTwZ+ea2+V5KkrllVAKiqXUk+DDyCZgj+qVX1z4Oun+Qs4ERga5J9wMuAO7Sf/SbgicAzk9wEfBd4cjXPK74pyXOAj9I8gviMqrpkNbVLkqRbrfppgMDNwC00AeCW1axYVaetsPx1wOuWWPYh4EOr+T5JkrS4VZ0DkOS5QA/YCtwD+KskvzmMwiRJ0vCsdgTgacBDq+o7AEn+N/AZ4P+sdWGSJGl4VnsVQGgOAcy7mcUv0ZMkSRvYakcA3gZ8dsFVAG9d25IkSdKwrfYqgFcn+UfgBJq//Fd1FYAkSdoYVn0VQFVdAFwwhFokSdKIDBQAknyyqh6e5AYOvAVvgKqqHx5KdZIkaSgGOgmw3fkHeEBV/XDf63B3/tIm1OvB5CRs2dJMe731rkjSiA18FUB7R76zV+woaWPr9WBmBvbsgapmOjNjCJA6ZrWXAZ6X5CFDqUTSaOzcCXNzB7bNzTXtkjpjtScBngQ8I8mVwHe49RyAB651YZKGZO/e1bVLGkurDQCnDqUKSaOzfXsz7L9Yu6TOWO0hgL00TwI8var20FwRcM81r0rS8OzaBRMTB7ZNTDTtkjpjtQHgDcDDgPmn+t0AvH5NK5I0XNPTMDsLO3ZA0kxnZ5t2SZ2x2kMAD62qY5P8M0BVXZ/kjkOoS9IwTU+7w5c6brUjAN9PcgjtzYCSbANuWfOqJEnSUK02ALyW5l4A90iyC/gk8Mo1r0qSJA3Vah8G1EtyAfBImksAH19Vlw2lMkmSNDSDPgvgMOAZwI8DFwNvrqqbhlmYJEkankEPAbwdmKLZ+Z8K/OnQKpIkSUM36CGA+1fVTwEkeSvwueGVJEmShm3QEYDvz79x6F+SpM1v0BGAByX5Vvs+wA+18/PPAvCRwJIkbSIDBYCqOmTYhUiSpNFZ7X0AJEnSGDAASJLUQQYASZI6yAAgSVIHGQAkSeogA4AkSR1kAJAkqYMMAJIkddBIA0CSM5Jcm+QLSyyfTnJR+/p0kgf1LbsyycVJLkyye3RVS5I0fkY9AnAmcMoyy/8N+JmqeiDwCmB2wfKTquqYqpoaUn2SJHXCoM8CWBNVdW6SyWWWf7pv9jzgyGHXJElSF23kcwCeBny4b76AjyW5IMnMOtUkSdJY2JABIMlJNAHgd/uaT6iqY4FTgWcn+ekl1p1JsjvJ7v3794+gWnVerweTk7BlSzPt9da7Ikla0YYLAEkeCLwFeFxVfX2+vaqubqfXAmcDxy22flXNVtVUVU1t27ZtFCWry3o9mJmBPXugqpnOzBgCJG14GyoAJNkOvAd4SlV9qa/9zkkOn38PnAwseiWBNFI7d8Lc3IFtc3NNuyRtYCM9CTDJWcCJwNYk+4CXAXcAqKo3AS8FfhR4QxKAm9oz/u8JnN22HQq8s6o+MsrapUXt3bu6dknaIEZ9FcBpKyx/OvD0RdqvAB502zWkdbZ9ezPsv1i7JG1gG+oQgLTp7NoFExMHtk1MNO2StIEZAKTbY3oaZmdhxw5ImunsbNMuSRvYSA8BSGNpetodvqRNxxEASZI6yAAgSVIHGQAkSeogA4AkSR1kAJAkqYMMAJIkdZABQJKkDjIASJLUQQYASZI6yAAgSVIHGQAkSeogA4AkSR1kAJAkqYMMAJIkdZABQJKkDjIASJLUQQYASZI6yAAgSVIHGQAkSeogA4AkSR1kAJAkqYMMAJIkdZABQJKkDjIASJLUQQYASZI6yAAgSVIHGQAkSeogA4A2ll4PJidhy5Zm2uutd0WSNJYOXe8CpB/o9WBmBubmmvk9e5p5gOnp9atLksbQSEcAkpyR5NokX1hieZK8NsnlSS5KcmzfstOTfLl9nT66qjUyO3feuvOfNzfXtEuS1tSoDwGcCZyyzPJTgaPb1wzwRoAkdwdeBjwUOA54WZK7DbVSjd7evatrlyQdtJEGgKo6F7humS6PA95RjfOAuya5F/BzwDlVdV1VXQ+cw/JBQpvR9u2ra5ckHbSNdhLgEcBVffP72ral2jVOdu2CiYkD2yYmmnZJ0praaAEgi7TVMu23/YBkJsnuJLv379+/psVpyKanYXYWduyApJnOznoCoCQNwUYLAPuAo/rmjwSuXqb9Nqpqtqqmqmpq27ZtQytUQzI9DVdeCbfc0kzd+UvSUGy0APB+4FfbqwGOB75ZVdcAHwVOTnK39uS/k9s2SZJ0EEZ6H4AkZwEnAluT7KM5s/8OAFX1JuBDwM8DlwNzwFPbZdcleQVwfvtRL6+q5U4mlCRJyxhpAKiq01ZYXsCzl1h2BnDGMOqSJKlrNtohAEmSNAIGAEmSOsgAIElSBxkAJEnqIAOAJEkdZACQJKmDDACSJHWQAUCSpA4yAEiS1EEGAEmSOsgAIElSBxkAJEnqIAOAJEkdZACQJKmDDACSJHWQAUCSpA4yAEiS1EEGAEmSOsgAIElSBxkAJEnqIAOAJEkdZACQJKmDDACSJHWQAWCc9XowOQlbtjTTXm+9K5IkbRCHrncBGpJeD2ZmYG6umd+zp5kHmJ5ev7okSRuCIwDjaufOW3f+8+bmmnZJUucZAMbV3r2ra5ckdYoBYFxt3766dklSpxgAxtWuXTAxcWDbxETTLknqPAPAuJqehtlZ2LEDkmY6O+sJgJIkwKsAxtv0tDt8SdKiHAGQJKmDDACSJHXQyANAklOSfDHJ5UletMjy1yS5sH19Kck3+pbd3Lfs/aOtXJKk8THScwCSHAK8Hng0sA84P8n7q+rS+T5V9Tt9/X8TeHDfR3y3qo4ZVb2SJI2rUY8AHAdcXlVXVNWNwLuAxy3T/zTgrJFUJklSh4w6ABwBXNU3v69tu40kO4B7A//Q13xYkt1Jzkvy+CXWm2n77N6/f/9a1S1J0lgZdQDIIm21RN8nA39bVTf3tW2vqingl4E/S3Kf23xY1WxVTVXV1LZt225/xZIkjaFRB4B9wFF980cCVy/R98ksGP6vqqvb6RXAJzjw/ABJkjSgUQeA84Gjk9w7yR1pdvK3OZs/yX8D7gZ8pq/tbknu1L7fCpwAXLpwXUmStLKRXgVQVTcleQ7wUeAQ4IyquiTJy4HdVTUfBk4D3lVV/YcH7ge8OcktNMHlj/uvHpAkSYPLgfvY8TI1NVW7d+9e7zIkSRqZJBe058styzsBSpLUQQYASZI6yAAgSVIHGQAkSeogA4AkSR1kAJAkqYMMAJIkdZABQJKkDjIASJLUQQYASZI6yAAgSVIHGQAkSeogA4AkSR1kAJAkqYMMAJIkdZABYFC9HkxOwpYtzbTXW++KJEk6aIeudwGbQq8HMzMwN9fM79nTzANMT69fXZIkHSRHAAaxc+etO/95c3NNuyRJm5ABYBB7966uXZKkDc4AMIjt21fXLknSBmcAGMSuXTAxcWDbxETTLknSJmQAGMT0NMzOwo4dkDTT2VlPAJQkbVpeBTCo6Wl3+JKkseEIgCRJHWQAkCSpgwwAkiR1kAFAkqQOMgBIktRBBgBJkjrIACBJUgcZACRJ6iADgCRJHWQAkCSpg1JV613D0CTZD+xZ7zo2gK3A19a7iA5wO4+G23k03M6jMYztvKOqtq3UaawDgBpJdlfV1HrXMe7czqPhdh4Nt/NorOd29hCAJEkdZACQJKmDDADdMLveBXSE23k03M6j4XYejXXbzp4DIElSBzkCIElSBxkAxkSSU5J8McnlSV60yPLnJbk0yUVJ/j7JjvWocxwMsK2fkeTiJBcm+WSS+69HnZvdStu5r98Tk1QSz1g/CAP8Pv9akv3t7/OFSZ6+HnVudoP8Pid5Uvv/9CVJ3jn0mjwEsPklOQT4EvBoYB9wPnBaVV3a1+ck4LNVNZfkmcCJVfVL61LwJjbgtv7hqvpW+/6xwLOq6pT1qHezGmQ7t/0OBz4I3BF4TlXtHnWtm9mAv8+/BkxV1XPWpcgxMOB2Phr4a+Bnq+r6JPeoqmuHWZcjAOPhOODyqrqiqm4E3gU8rr9DVX28quba2fOAI0dc47gYZFt/q2/2zoApe/VW3M6tVwB/AnxvlMWNkUG3s26fQbbzbwCvr6rrAYa98wcDwLg4Ariqb35f27aUpwEfHmpF42ugbZ3k2Um+QrNz+q0R1TZOVtzOSR4MHFVVfzfKwsbMoP93/M/28OHfJjlqNKWNlUG2832B+yb5VJLzkgx91NAAMB6ySNuif3Um+RVgCnjVUCsaXwNt66p6fVXdB/hd4CVDr2r8LLudk2wBXgM8f2QVjadBfp8/AExW1QOB/we8fehVjZ9BtvOhwNHAicBpwFuS3HWYRRkAxsM+oD+VHwlcvbBTkkcBO4HHVtV/jqi2cTPQtu7zLuDxQ61oPK20nQ8HfhL4RJIrgeOB93si4Kqt+PtcVV/v+//iL4D/PqLaxskg/2/sA95XVd+vqn8DvkgTCIbGADAezgeOTnLvJHcEngy8v79DO1z6Zpqd/9CPLY2xQbZ1/z/aXwC+PML6xsWy27mqvllVW6tqsqomac5reawnAa7aIL/P9+qbfSxw2QjrGxcrbmfgvcBJAEm20hwSuGKYRR06zA/XaFTVTUmeA3wUOAQ4o6ouSfJyYHdVvZ9myP8uwN8kAdhbVY9dt6I3qQG39XPa0ZbvA9cDp69fxZvTgNtZt9OA2/m32qtZbgKuA35t3QrepAbczh8FTk5yKXAz8MKq+vow6/IyQEmSOshDAJIkdZABQJKkDjIASJLUQQYASZI6yAAgSVIHGQCkDSLJze3T1r6Q5G+STKzR5357gD6fGMZNdIbxuUnumuRZffMnJhnodsBJ/izJT6/iu7Yl+cjB1CltdAYAaeP4blUdU1U/CdwIPGO9C9qg7go8a8VeCyS5O3B8VZ076DpVtR+4JskJq/0+aaMzAEgb0z8BPw7N8xuSfK4dHXhz+2hRknw7ya4k/9I+POSebfu9k3wmyflJXjH/gQv/Uk7yuvZRrwfoHzFI8sQkZ7bvz0zyxiQfT3JFkp9JckaSy+b7LCfJyW1dn29HOO7Stl+Z5A/b9ouT/ETbvi3JOW37m5Psae+Q9sfAfdrtMf9Mi7u0D6r51yS9tHe7WuCJwA/+mm+/95VtTbuTHJvko0m+kqQ/fL0XmF7p55M2GwOAtMEkORQ4Fbg4yf2AXwJOqKpjaO4QNr8zujNwXlU9CDiX5nGiAH8OvLGqHgL8xxqXdzfgZ4HfoXlIzGuABwA/leSYZX6mrTQPRXpUVR0L7Aae19fla237G4EXtG0vA/6hbT8b2N62vwj4Sjta8sK27cHAbwP3B/4rsNhf7CcAFyxou6qqHkYTuM6kCQnHAy/v67MbeMRSP5u0WRkApI3jh5JcSLPD2Qu8FXgkzcNXzm+XPZJmBwfNYYL5v+gvACbb9ycAZ7Xv/3KNa/xANbcPvRj4alVdXFW3AJf0ff9ijqfZOX+q/TlOB3b0LX9PO+3/OR5O8zAlquojNLdVXsrnqmpfW8uFS9RyL2D/grb5WwpfDHy2qm5oh/2/1/cktmuB/7LMd0ubks8CkDaO77Z/5f9AO5T99qp68SL9v1+33sv7Zg7897zYPb5v4sDQf9gSddQyfeafCndL3/v5+eX+PwmtP/ltAAABcElEQVRwTlWdtsTy+c/q/zkWG8ZfSn8tC7fFvO9ycD/PYe260lhxBEDa2P4eeGKSe0BzIluSHSus8ymap43Bgceu9wD3T3KnJD9CM5qwmK8muV+SLcATbkft/c4DTkgyf17DRJL7rrDOJ4Entf1Ppjn8AHADzeOAV+sy2vMqVum+wBcOYj1pQzMASBtYVV1Kc+z8Y0kuAs6hGcpeznOBZyc5H/iRvs+6Cvhr4CKgB/zzEuu/iObQwj8A19yuH+DW795P8xS5s9qf4zzgJ1ZY7Q9pno72eZpzIq4BbmifkPap9nLJVy37CQf6IHDiamuneUTrBw9iPWlD82mAkjakJHcCbm4fpfowmhMblzzRcMDP/CTwmKr6xirWORd4XFUtdw6CtOl4DoCkjWo78NftoYgbufUqh9vj+e3nDhQAkmwDXu3OX+PIEQBJkjrIcwAkSeogA4AkSR1kAJAkqYMMAJIkdZABQJKkDjIASJLUQf8flHdSR7tEb7AAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.optimize import leastsq\n", "%matplotlib inline\n", "#\n", "Debug = False\n", "printData = True\n", "#\n", "nLen = 10\n", "nMeasPerLen = 20\n", "nOscPerMeas = 20\n", "Tdata = np.zeros((nLen, nMeasPerLen))\n", "lenBot = 0.15\n", "lenTop = 0.60\n", "lengths = np.linspace(lenBot, lenTop, nLen)\n", "g = 9.81\n", "Tideal = 2*np.pi*np.sqrt(lengths/g)\n", "sigma = 0.5/nOscPerMeas\n", "print(\" \")\n", "print(\"lengths (m)\\n\",lengths)\n", "print(\"Tideal (s)\\n\",Tideal)\n", "print(\"sigma = {:5.3f} s\".format(sigma))\n", "for n in range(0, nLen):\n", " Tdata[n, 0:nMeasPerLen] = np.random.normal(Tideal[n], sigma, nMeasPerLen)\n", "#\n", "if (printData):\n", " print(\" \")\n", " print(\"Length time measurements:\")\n", " np.set_printoptions(precision = 3)\n", " for n in range(0, nLen):\n", " print(\" \")\n", " print(\"lengths[{:d}] = {:5.4f}\".format(n, lengths[n]))\n", " for i in range(0, nMeasPerLen):\n", " print(\"Tdata[{:d}, {:d}] = {:5.4f}\".\\\n", " format(n, i, Tdata[n, i]))\n", "\n", "Tmean = np.zeros(nLen)\n", "for n in range(0, nLen):\n", " Tmean[n] = np.cumsum(Tdata, 1)[n, nMeasPerLen - 1]/nMeasPerLen\n", "#\n", "if (Debug):\n", " print(\" \")\n", " print(\"Raw pendulum data:\")\n", " np.set_printoptions(precision = 5)\n", " for n in range(0, nPoints):\n", " print(\" \")\n", " print(\"Length {:2.3f} m\".format(lengths[n]))\n", " print(\"Periods (s):\\n\",Tdata[n, 0:nData])\n", " #\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], \\\n", " Tmean[n]**2))\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.xlim(0, 1.2*np.max(lengths))\n", "#plt.ylim(0, 1.2*np.max(Taverage**2))\n", "plt.title('MC data', fontsize = 14)\n", "plt.xlabel('Pendulum length (m)')\n", "plt.ylabel('Period$^2$ (s$^2$)')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculate standard deviation of period measurements (taken over nOsc scillations), then values of $\\Delta T$ for each pendulum length." ] }, { "cell_type": "code", "execution_count": 39, "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.0254 0.0239 0.0218 0.0245 0.0264]\n", "Average standard deviation 0.0227 s\n", "Input standard deviation 0.025 s\n", "Errors 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": [ "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", "#\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(\"Average standard deviation {:5.4f} s\".format(np.sum(sDev)/nLen))\n", "print(\"Input standard deviation {:5.3f} s\".format(sigma))\n", "print(\"Errors in T in seconds\\n\",Terror)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculate the values of $\\Delta T^2$ and add error bars to graph. Add also your estimated errors in the measurements of the length of the pendulum." ] }, { "cell_type": "code", "execution_count": 40, "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": "iVBORw0KGgoAAAANSUhEUgAAAgAAAAGECAYAAABAsZipAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XuUZWdd5vHvQy4gIRqwi8skaRoxKBdJyBQBVtBUBGJQJLhkMBExMmiPQBTkMiLDWNWwcDmDCxwVgQiZgAvCEkMgDJcQlVQIGEg1trkChkhIk0gaEjpRUOjwmz/2bnJSqa461amzT9XZ389aZ9XZ7373Ob/aqfR5zrsvb6oKSZLUL/cadwGSJKl7BgBJknrIACBJUg8ZACRJ6iEDgCRJPWQAkCSphwwAkiT1kAFA0pKSnJ2k2sd3k9yc5BNJXpzkoFW8zkz7GptGWa+k1TEASFrO3wAPAbYAJwEfArYBn0xyyBjrknQPGQAkLec/qupfquqrVbWjqt4IzADHAv8dIMmvJLksye3tKMH7khzertsCfKJ9rV3tSMDZ7bqTk3wyya1JbklyQZJHdvz7Sb1lAJC0KlV1JfAx4BfbpoOBWeBo4BnAJuCcdt0NA/0eTTOa8JJ2+RDgj4HjaELFbuBDSQ4e7W8gCeDAcRcgaUO6GngqQFWdNdB+XZIXAtckOaKqdia5pV13c1V9fW/Hqjp38AWTPB+4jSYQXDLS6iU5AiBpvwQogCTHJvlgkuuT3A4stH02L/sCycOTvCfJl5LcBnyN5t+kZbeTtDYMAJL2x6Novu0fAlwAfAt4HvB44OS2z0pD+R8CpoD/BjwBeBywZ4jtJK0BA4CkVUnyGJoP+b8GfpzmmP+rq+riqvo88MBFm3yn/XnAwGv8MPBI4A+q6m+q6hrgUDwsKXXG/9kkLefeSR5M82VhCngK8GpgO/BHwH2B/wDOSPJmmg/11y16jetpDhf8XJIPAd8GbgW+DvxGkhuAw4E30IwASOqAIwCSlvNU4CbgK8DfAs+kuQ/AT1XVv1XVLuB04Fk0JwbOAi8bfIGq+mrb/nqa4/x/VlXfA34JeCxwJfBm4H/ShAlJHUhVjbsGSZLUMUcAJEnqIQOAJEk9ZACQJKmHDACSJPVQpwEgyZHtdKLXJLkqyUuW6DOTZHeSHe3j9wfWnZzkC0muTfKqLmuXJGmSdH0fgD3Ay6vqc0kOBbYnubCqrl7U75NV9YzBhiQH0Fwq9DRgJ3BZkvOX2Pb7Nm3aVFu2bFnb30CSpHVs+/btX6+qqZX6dRoAquommmuKqarbk1xDcwOQfX6IDzgOuLaqrgNI8l7glOW23bJlCwsLC/taLUnSxEly/TD9xnYOQDtP+OOAzyyx+klJ/jHJR5M8um07nGZq0b12tm2LX3drkoUkC7t27VrjqiVJmgxjCQBJ7gecC7y0qm5btPpzwEOr6mjgT4EP7N1siZe6212MqurMqpququmpqRVHQCRJ6qXOA0CSg2g+/N9dVe9fvL6qbquqf22ffwQ4KMkmmm/8Rw50PQK4sYOSJUmaOF1fBRDgHcA1VfXGffR5cNuPJMfR1PgN4DLgqCQPS3IwcCpwfjeVS5I0Wbq+CuB4mjnDr0iyo217NbAZoKreCjwbeGGSPTSzhp1azYQFe5KcQTP3+AHAWVV1Vcf1S5I0ESZ6MqDp6enyKgBJUp8k2V5V0yv1806AkiT1kAFAkqQeMgBIktRDBgBJknrIACBJUg8ZACRJ6iEDgCRJXZuZgTHPVmsAkCRpHK4fatK+kTEASJLUpZkZ2LFjxW6jZgCQJKkrc3MwPw+7dzfLSfOYm+u8FAOAJEldmZuDKjjhhGa5qnkYACRJUhcMAJIkde2ii2B2dqwlGAAkSRqHMQz7DzIASJLUQwYASZJ6yAAgSVIPGQAkSeohA4AkST1kAJAkqYcMAJIk9ZABQJKkHjIASJLUQwYASZJ6yAAgSVIPGQAkSeohA4AkST1kAJAkqYcMAJIk9ZABQJKkHjIASJLUQ50GgCRHJvlEkmuSXJXkJUv0eW6Sy9vHp5McPbDuy0muSLIjyUKXtUuSNEkO7Pj99gAvr6rPJTkU2J7kwqq6eqDPPwMnVNWtSZ4OnAk8YWD9iVX19Q5rliRp4nQaAKrqJuCm9vntSa4BDgeuHujz6YFNLgWO6LJGSZL6YGznACTZAjwO+Mwy3V4AfHRguYCPJ9meZOvoqpMkabJ1fQgAgCT3A84FXlpVt+2jz4k0AeDJA83HV9WNSR4IXJjk81V18aLttgJbATZv3jyS+iVJ2ug6HwFIchDNh/+7q+r9++jzWODtwClV9Y297VV1Y/vzZuA84LjF21bVmVU1XVXTU1NTo/gVJEna8Lq+CiDAO4BrquqN++izGXg/8Lyq+uJA+yHtiYMkOQQ4Cbhy9FVLkjR5uj4EcDzwPOCKJDvatlcDmwGq6q3A7wM/DPx5kxfYU1XTwIOA89q2A4H3VNXHui1fkqTJ0PVVAJcAWaHPrwO/vkT7dcDRd99CkiStlncClCSphwwAkqSNb25u3BVsOAYASdLGNjMD27aNu4oNxwAgSdq4ZmZgx447n2toBgBJ0sY0Nwfz87B7d7M8Pw+JhwOGZACQJG1Mc3NQBSec0CxXNQ8DwFAMAJIk9ZABQJK0sV10EczOjruKDccAIEna+Bz2XzUDgCRJPWQAkCSphwwAkiT1kAFAkqQeMgBIktRDBgBJknrIACBJUg8ZACRJ6iEDgCRJPWQAkCSphwwAkiT1kAFAkqQeMgBIktRDBgBJknrIACBJUg8ZACRJ6iEDgCRJPWQAkCSphwwAkiT1kAFAkqQeMgBIktRDBgBJ0p1mZmDLlnFXoQ4YACRJd3X99eOuQB3oNAAkOTLJJ5Jck+SqJC9Zok+S/EmSa5NcnuTYgXWnJ/mn9nF6l7VLUi/s2NH8nJkZaxkava5HAPYAL6+qRwJPBF6c5FGL+jwdOKp9bAXeApDkAcAs8ATgOGA2yf27KlySJtrcHCSwe3ezPD/fLM/NjbMqjVCnAaCqbqqqz7XPbweuAQ5f1O0U4F3VuBQ4LMlDgJ8BLqyqW6rqVuBC4OQOy5ekyTU3B1VwwgnNclXzMABMrLGdA5BkC/A44DOLVh0O3DCwvLNt21f74tfdmmQhycKuXbvWsmRJmnwXXTTuCtSRsQSAJPcDzgVeWlW3LV69xCa1TPtdG6rOrKrpqpqempq658VKUt/Mzo67AnWg8wCQ5CCaD/93V9X7l+iyEzhyYPkI4MZl2iVJa8lh/17o+iqAAO8ArqmqN+6j2/nAr7ZXAzwR2F1VNwEXACcluX978t9JbZskSVqlAzt+v+OB5wFXJGmvNeHVwGaAqnor8BHgZ4FrgW8Bz2/X3ZLkdcBl7XavrapbOqxdkqSJ0WkAqKpLWPpY/mCfAl68j3VnAWeNoDRJknrFOwFKktRDBgBJknrIACBJUg8ZACRJ6iEDgCRJPWQAkCSphwwAkiT1kAFAkqQeMgBIktRDBgBJknrIACBJUg8ZACRJ6iEDgCRJPWQAkCSphwwAkiT1kAFAkqQeMgBIktRDBgBJknrIACBJUg8ZACRJ6iEDgCSN2tzcuCuQ7sYAIEmjNDMD27aNuwrpbgwAkjQqMzOwY8e4q5CWZACQpFGYm4P5edi9u1lOmoeHA7ROGAAkaRTm5qAKTjihWa5qHgYArRMGAEmSesgAIEmjdNFFMDs77iqkuzEASNKoOeyvdcgAIElSDxkAJEnqIQOAJEk9dGCXb5bkLOAZwM1V9Zgl1r8SeO5AbY8EpqrqliRfBm4H7gD2VNV0N1VLkjR5uh4BOBs4eV8rq+oNVXVMVR0D/B4wX1W3DHQ5sV3vh78kSfdApwGgqi4GblmxY+M04JwRliNJUm+ty3MAktyXZqTg3IHmAj6eZHuSreOpTJKkydDpOQCr8PPApxYN/x9fVTcmeSBwYZLPtyMKd9GGg60Amzdv7qZaSZI2mHU5AgCcyqLh/6q6sf15M3AecNxSG1bVmVU1XVXTU1NTIy9UkqSNaN0FgCQ/BJwAfHCg7ZAkh+59DpwEXDmeCiVJ2vi6vgzwHGAG2JRkJzALHARQVW9tu/0C8PGq+reBTR8EnJcEmprfU1Uf66puSZImTacBoKpOG6LP2TSXCw62XQccPZqqJEnqn3V3CECSJI2eAUCSpB4yAEiS1EMGAEmSesgAIElSDxkAJEnqIQOAJEk9tF8BoL0z3wFrXYwkSerGUAEgyb2S/HKSDye5Gfg8cFOSq5K8IclRoy1TkiStpWFHAD4BPBz4PeDBVXVkVT0Q+EngUuAPk/zKiGqUJElrbNhbAT+1qr67uLGdrvdc4NwkB61pZZIkaWSGGgFY6sN/f/pIkqT1YcUAkORpSf4iyTHt8tbRlyVJkkZpmEMALwKeD7wmyQOAY0ZbkiRJGrVhDgHsqqpvVtUrgJOAx4+4JkmSNGLDBIAP731SVa8C3jW6ciRJUhdWDABV9cG9z5P8F+Ds9vlrkrw/ybGjK0+SJI3Cau8E+D+r6vYkTwZ+Bngn8Ja1L0uSljAzA1u2jLsKaSKsNgDc0f78OeAt7ejAwWtbkiQt4/rrx12BNBFWGwC+muRtwHOAjyS59368hiSt3swM7Nhx53NJ98hqP7yfA1wAnFxV3wQeALxyzauSpEFzczA/D7t3N8vz85A07ZL2S6pq5U5JaoWOw/Tp2vT0dC0sLIy7DElrZWam+fBfX//USOtKku1VNb1Sv6EnA0ryW0k2L3qTg5P8dJJ3AqfvT6GSNLSLLhp3BdLEGHYyoJOB/wqck+RhwDeB+wAHAB8H3lRVO0ZToiQNmJ0ddwXSRBjqEMBdNmhm/dsEfLs9D2Dd8hCAJKlvhj0EMOwIwPe1s/7dtF9VSZKkdcFL+CRJ6iEDgCRJPWQAkCSph1YMAEmeluQvkhzTLm8dfVmSJGmUhjkJ8EXA84HXJHkAcMxoS5IkSaM2zCGAXVX1zap6BXAS8PgR1yRJkkZsmADw4b1PqupVwLv2982SnJXk5iRX7mP9TJLdSXa0j98fWHdyki8kuTbJq/a3BkmSNMQhgKr6YJKXDbYtXq6qNw75fmcDf8byIeKTVfWMRe93APBm4GnATuCyJOdX1dVDvq8kSRow7I2ADm1//hjNIYDz2+WfBy4e9s2q6uIkW4btP+A44Nqqug4gyXuBUwADgCRJ+2GoAFBV2wCSfBw4tqpub5fngPetcU1PSvKPwI3AK6rqKuBw4IaBPjuBJyy1cXuVwlaAzZs3L9VFkqTeW+19ADYD3xlY/g6wZc2qgc8BD62qo4E/BT7QtmeJvktOYlBVZ1bVdFVNT01NrWFpkiRNjtXOBfCXwGeTnNcuPwt451oVU1W3DTz/SJI/T7KJ5hv/kQNdj6AZIZAkSfthVQGgql6f5KPAT9J8A39+Vf3DWhWT5MHA16qqkhxHM0LxDZrph49qpyL+KnAq8Mtr9b6SJPXNqmcDBO4AvkcTAL63mg2TnAPMAJuS7ARmgYMAquqtwLOBFybZA3wbOLWa+Yr3JDkDuAA4ADirPTdAkiTthzSfr0N2Tl4C/AZwLs1x+V8AzqyqPx1NeffM9PR0LSwsjLsMSZI6k2R7VU2v1G+1IwAvAJ5QVf/Wvsn/Av6e5oQ9SZK0Qaz2KoDQHALY6w6WPkNfkiStY6sdAfi/wGcWXQXwjrUtSZIkjdpqrwJ4Y5J54Hiab/5rehWAJEnqxqqvAqiq7cD2EdQiSZI6MlQASHJJVT05ye3c9Q58AaqqfnAk1UmSpJEYdi6AJycJ8Oiq+sqIa5IkSSM29FUA7Q15zluxo6SNY25u3BVIGpPVXgZ4aZLHj6QSSd2amYFt28ZdhaQxWW0AOJEmBHwpyeVJrkhy+SgKkzRCMzOwY8edzyX1zmqvAnj6SKqQ1J25OZifv3N5fh4SmJ31kIDUI6sdAfgKzUyAp1fV9TRXBDxozauSNDpzc1DVPODO5374S72y2gDw58CTgNPa5duBN69pRZIkaeRWewjgCVV1bJJ/AKiqW5McPIK6JHVhdnbcFUgak9WOAHw3yQG0NwNKMgV8b82rktQNh/2l3lptAPgTmnsBPDDJ64FLgD9Y86okSdJIrXYyoHcn2Q48heY2wM+qqmtGUpkkSRqZYecCuA/wm8CPAlcAb6uqPaMsTJIkjc6whwDeCUzTfPg/HfijkVUkSZJGbthDAI+qqp8ASPIO4LOjK0mSJI3asCMA3937xKF/SZI2vmFHAI5Oclv7PMAPtMuhmSjwB0dSnSRJGomhAkBVHTDqQiRJUndWex8ASZI0AQwAkiT1kAFAkqQeMgBIktRDBgBJknrIACBJUg8ZACRJ6iEDgCRJPWQAkCSphzoNAEnOSnJzkiv3sf65SS5vH59OcvTAui8nuSLJjiQL3VUtSdLk6XoE4Gzg5GXW/zNwQlU9FngdcOai9SdW1TFVNT2i+iRJ6oVhJwNaE1V1cZIty6z/9MDipcARo65JkqQ+Ws/nALwA+OjAcgEfT7I9ydZ9bZRka5KFJAu7du0aeZHSXczNjbsCSRrKugwASU6kCQC/O9B8fFUdCzwdeHGSn1pq26o6s6qmq2p6amqqg2ql1swMbNs27iokaSjrLgAkeSzwduCUqvrG3vaqurH9eTNwHnDceCqUljAzAzt2jLsKSRraugoASTYD7weeV1VfHGg/JMmhe58DJwFLXkkgdW5uDubnYffuZjlpHh4OkLSOpaq6e7PkHGAG2AR8DZgFDgKoqrcmeTvwi8D17SZ7qmo6yY/QfOuH5sTF91TV61d6v+np6VpY8IpBdWRmpgkCHf4/JUmLJdk+zNVynQaArhkA1LnEACBprIYNAOvqEIC04c3OjrsCSRqKAUBaSx73l7RBGAAkSeohA4AkST1kAJAkqYcMAJIk9ZABQJKkHjIASJLUQwYASZJ6yAAgSVIPGQAkSeohA4AkST1kAJAkqYcMAJIk9ZABQJKkHjIASJLUQwYASZJ6yAAgSVIPGQAkSeohA4AkST1kAJAkqYcMAJIk9ZABQJKkHjIASJLUQwYASZJ6yAAgSVIPGQAkSeohA4AkST1kAND6NTc37gokaWIZALR+bds27gokaWIZALT+zMzAYYfd+VyStOY6DQBJzkpyc5Ir97E+Sf4kybVJLk9y7MC605P8U/s4vbuq1am5OZifh927m+X5eUg8HCBJa6zrEYCzgZOXWf904Kj2sRV4C0CSBwCzwBOA44DZJPcfaaUaj7k5qIITTmiWq5qHAUCS1lSnAaCqLgZuWabLKcC7qnEpcFiShwA/A1xYVbdU1a3AhSwfJLTRXXTRuCuQpIm23s4BOBy4YWB5Z9u2r3ZNstnZcVcgSRNrvQWALNFWy7Tf/QWSrUkWkizs2rVrTYtTxxz2l6SRWW8BYCdw5MDyEcCNy7TfTVWdWVXTVTU9NTU1skIlSdrI1lsAOB/41fZqgCcCu6vqJuAC4KQk929P/jupbZMkSfvhwC7fLMk5wAywKclOmjP7DwKoqrcCHwF+FrgW+Bbw/HbdLUleB1zWvtRrq2q5kwklSdIyOg0AVXXaCusLePE+1p0FnDWKuiRJ6pv1dghAkiR1wAAgSVIPGQAkSeohA4AkST1kAJAkqYcMAJIk9ZABQJKkHjIASJLUQwYASZJ6yAAgSVIPGQAkSeohA4AkST1kAJAkqYcMAJIk9ZABQJKkHjIASJLUQwYASZJ6yAAgSVIPGQAkSeohA4AkST1kAJAkqYcMAJIk9ZABQJKkHjIA9MXc3LgrkCStIwaAPpiZgW3bxl2FJGkdMQBMupkZ2LHjzueSJGEAmGxzczA/D7t3N8vz85B4OECSRKpq3DWMzPT0dC0sLIy7jPUhgQn+by1JaiTZXlXTK/VzBECSpB4yAPTF7Oy4K5AkrSMGgL7wuL8kaYABQJKkHuo8ACQ5OckXklyb5FVLrH9Tkh3t44tJvjmw7o6Bded3W7kkSZPjwC7fLMkBwJuBpwE7gcuSnF9VV+/tU1W/M9D/t4DHDbzEt6vqmK7qlSRpUnU9AnAccG1VXVdV3wHeC5yyTP/TgHM6qUySpB7pOgAcDtwwsLyzbbubJA8FHgb83UDzfZIsJLk0ybP2sd3Wts/Crl271qpuSZImStcBIEu07evuNKcCf11Vdwy0bW5vbvDLwB8nefjdXqzqzKqarqrpqampe16xJEkTqOsAsBM4cmD5CODGffQ9lUXD/1V1Y/vzOuAi7np+gCRJGlLXAeAy4KgkD0tyMM2H/N3O5k/yY8D9gb8faLt/knu3zzcBxwNXL95WkiStrNOrAKpqT5IzgAuAA4CzquqqJK8FFqpqbxg4DXhv3XWigkcCb0vyPZrg8oeDVw9IkqThORmQJEkTxMmAJEnSPhkAJEnqIQOAJEk9ZACQJKmHDACSJPWQAUCSpB4yAEiS1EMGAEmSesgAIElSDxkAJEnqIQOAJEk9ZACQJKmHDACSJPWQAUCSpB4yAEiS1EMGgP0xNzfuCiRJukcMAPtj27ZxVyBJ0j1iAFitww5rfs7MjLUMSZLuCQPAsObmIIHdu5vl+flm2cMBkqQNKFU17hpGZnp6uhYWFtb+hROY4P0mSdq4kmyvqumV+jkCIElSDxkA9sfs7LgrkCTpHjEA7A+P+0uSNjgDgCRJPWQAkCSphwwAkiT1kAFAkqQeMgBIktRDBgBJknrIACBJUg8ZACRJ6iEDgCRJPWQAkCSphwwAkiT10ERPB5xkF3D9uOtYJzYBXx93ET3gfh499/HouY+7Mar9/NCqmlqp00QHAN0pycIw80PrnnE/j577ePTcx90Y9372EIAkST1kAJAkqYcMAP1x5rgL6An38+i5j0fPfdyNse5nzwGQJKmHHAGQJKmHDAATJsnJSb6Q5Nokr1pi/cuSXJ3k8iR/m+Sh46hzIxtiH/9mkiuS7EhySZJHjaPOjW6l/TzQ79lJKolnra/SEH/Lv5ZkV/u3vCPJr4+jzo1smL/jJM9p/12+Ksl7OqvNQwCTI8kBwBeBpwE7gcuA06rq6oE+JwKfqapvJXkhMFNVvzSWgjegIffxD1bVbe3zZwIvqqqTx1HvRjXMfm77HQp8GDgYOKOqFrqudaMa8m/514DpqjpjLEVucEPu46OAvwJ+uqpuTfLAqrq5i/ocAZgsxwHXVtV1VfUd4L3AKYMdquoTVfWtdvFS4IiOa9zohtnHtw0sHgKYsldvxf3ceh3wv4F/77K4CTHsPtb+G2Yf/wbw5qq6FaCrD38wAEyaw4EbBpZ3tm378gLgoyOtaPIMtY+TvDjJl2g+nH67o9omyYr7OcnjgCOr6v91WdgEGfbfi19sDxn+dZIjuyltYgyzjx8BPCLJp5JcmqSz0UIDwGTJEm1LfvtM8ivANPCGkVY0eYbax1X15qp6OPC7wGtGXtXkWXY/J7kX8Cbg5Z1VNHmG+Vv+ELClqh4L/A3wzpFXNVmG2ccHAkcBM8BpwNuTHDbiugADwKTZCQwm9COAGxd3SvJU4H8Az6yq/+iotkkx1D4e8F7gWSOtaDKttJ8PBR4DXJTky8ATgfM9EXBVVvxbrqpvDPwb8RfAf+6otkkxzL8XO4EPVtV3q+qfgS/QBIKRMwBMlsuAo5I8LMnBwKnA+YMd2mHTt9F8+Hd2rGmCDLOPB//n/Tngnzqsb1Isu5+randVbaqqLVW1heZ8lmd6EuCqDPO3/JCBxWcC13RY3yRYcR8DHwBOBEiyieaQwHVdFHdgF2+iblTVniRnABcABwBnVdVVSV4LLFTV+TRD/vcD3pcE4CtV9cyxFb3BDLmPz2hHWb4L3AqcPr6KN6Yh97PugSH38W+3V7LsAW4Bfm1sBW9AQ+7jC4CTklwN3AG8sqq+0UV9XgYoSVIPeQhAkqQeMgBIktRDBgBJknrIACBJUg8ZACRJ6iEDgLROJLmjnXHtyiTvS3LfNXrdfx2iz0WjuInOKF43yWFJXjSwPJNkqNsBJ/njJD+1iveaSvKx/alTWu8MANL68e2qOqaqHgN8B/jNcRe0Th0GvGjFXoskeQDwxKq6eNhtqmoXcFOS41f7ftJ6ZwCQ1qdPAj8KzbwNST7bjg68rZ1ilCT/muT1Sf6xnUTkQW37w5L8fZLLkrxu7wsu/qac5M/a6V7vYnDEIMmzk5zdPj87yVuSfCLJdUlOSHJWkmv29llOkpPauj7XjnDcr23/cpJtbfsVSX68bZ9KcmHb/rYk17d3SvtD4OHt/tg7l8X92slqPp/k3WnvcrXIs4Hvf5tv3/cP2poWkhyb5IIkX0oyGL4+ADx3pd9P2mgMANI6k+RA4OnAFUkeCfwScHxVHUNzp7C9H0aHAJdW1dHAxTTTigL8H+AtVfV44F/WuLz7Az8N/A7NRDFvAh4N/ESSY5b5nTbRTIr01Ko6FlgAXjbQ5ett+1uAV7Rts8Dfte3nAZvb9lcBX2pHS17Ztj0OeCnwKOBHgKW+sR8PbF/UdkNVPYkmcJ1NExKeCLx2oM8C8JP7+t2kjcoAIK0fP5BkB80HzleAdwBPoZmA5bJ23VNoPuCgOUyw9xv9dmBL+/x44Jz2+V+ucY0fqub2oVcAX6uqK6rqe8BVA++/lCfSfDh/qv09TgceOrD+/e3Pwd/jyTSTKVFVH6O5rfK+fLaqdra17NhHLQ8Bdi1q23tL4SuAz1TV7e2w/78PzMh2M/CflnlvaUNyLgBp/fh2+y3/+9qh7HdW1e8t0f+7dee9vO/grv8/L3WP7z3cNfTfZx911DJ99s4M972B53uXl/v3JMCFVXXaPtbvfa3B32OpYfx9Gaxl8b7Y69vs3+9zn3ZbaaI4AiCtb38LPDvJA6E5kS3JQ1fY5lM0s47BXY9dXw88Ksm9k/wQzWjCUr6W5JFJ7gX8wj2ofdClwPFJ9p7XcN8kj1hhm0uA57T9T6I5/ABwO810wKt1De15Fav0CODK/dhOWtcMANI6VlVX0xw7/3iSy4ELaYayl/MS4MVJLgN+aOC1bgD+CrgceDfwD/vY/lU0hxb+DrjpHv0Cd773LpqZ5M5pf49LgR9fYbNtNLOkfY7mnIibgNvbmdI+1V4u+YZlX+GuPgzjjyL6AAAAiUlEQVTMrLZ2mqlaP7wf20nrmrMBSlqXktwbuKOdUvVJNCc27vNEwyFf8xLgGVX1zVVsczFwSlUtdw6CtOF4DoCk9Woz8FftoYjvcOdVDvfEy9vXHSoAJJkC3uiHvyaRIwCSJPWQ5wBIktRDBgBJknrIACBJUg8ZACRJ6iEDgCRJPWQAkCSph/4/n8pAk6uOnJMAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "T2error = 2*Tmean*Terror\n", "#\n", "lenError = np.zeros(nLen)\n", "lenError[0:nLen] = 0.002\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 = '',\\\n", " marker = '+') \n", "plt.title('Data', fontsize = 14)\n", "plt.xlabel('Pendulum length (m)')\n", "plt.ylabel('Period$^2$ (s$^2$)')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Carry out a $\\chi^2$ fit to pendulum data." ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.optimize import leastsq\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define the fit and error functions, carry out the fit, test whether it has worked and print out and plot the results if it has.\n", "\n", "The routine used to carry out the fit is [optimize.leastsq](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html). This allows the definition of an error function (called a cost function in the documentation) that makes it possible to cope with errors in both $x$ and $y$." ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " \n", "Fit quality:\n", "chisq per point = \n", " [0.5 0.1 0.5 2.1 0.2 0.6 0.5 0.2 0.9 0. ]\n", "chisq = 5.570 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": "iVBORw0KGgoAAAANSUhEUgAAAfUAAAGHCAYAAACposvbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd8VFX+//HXh6IYwbUANiTBXfTn6krEWMDG+nXtiq6oYATBElGxrF3RdS1Y0JVdK0YsqKOoqICIIBaaiBIQWBXbuoCxIiqggVByfn+ciYY4IRMyd+7Mnffz8ciDzMxN5nMN5O2593zOMeccIiIikv2ahF2AiIiIpIZCXUREJCIU6iIiIhGhUBcREYkIhbqIiEhEKNRFREQiQqEuIiISEQp1EckIZtbFzN4ys8lm9pSZNQ+7JpFso1AXkUyxEDjYOXcQ8BnQPeR6RLKOQl1E1svM3jezbut5fYGZHdLY93HOfemcWxF/uAaoquP9djazd81suZldkEyNIrlCoS4SkngYroiH049mNt3M+ptZUv8uUxWm9XHO7eqcm5Su9zSzDsARwNg6DrkcmOSca+WcuyuMGkUylUJdJFzHOOdaAfnArcAVwEPhlhQeM9sMGA70ds6tquOwfOD99FUlkj0U6iIZwDm31Dk3BjgZOM3MdgMwsyvN7L/x0fwHZnZ8/PnHgfbAi2b2k5ldvr7jazOzfmb2Yo3Hn5rZMzUef25mhfHPF5jZIXW9Z1yhmc0zs6Vm9rSZtajjfT8ys/+Y2Tbxx7vFv2YXM2sGPAX8wzn3UR1f/zrwZ+CeeA07NaBGkchTqItkEOfcO0A5cED8qf/GP/8dcD3whJlt65zrDSzCj/RbOucGr+/4BG81GTjAzJrEX28O7AdgZjsCLYF5tWqr6z0BTgIOBzoAuwN96zjFQuAn4Jh4iA8HbnbOzQd6AfsAfzezSWZ2coL/PgcDU4EB8Ro+bkCNIpGnUBfJPF8CWwI4556NTyCrcs49DXwC7F3XFyZ7vHPuM2A5PmQPAiYAX5jZ/4s/nuqcSzhRrQ53xd/3e+DF+PdNVN+K+Hv9CbgGPyHujvhrjzvnWjvnusU/nm7A+4sI0CzsAkTkN7YHvgcwsz7AxUBB/LWWQOu6vrCBx08GugF/iH/+Iz7Qu8QfN8TXNT6vALZbz7HvAYOBNsDezrm1DXwvEamDRuoiGcTM9sKH+jQzywceBAYAWznnNscHosUPd7W+tr7ja6sO9QPin0/Gh/pB1B3qro7nG+ID/GX6G+OX3VMtFTWKZCWFukgGMLPNzOxoYATwhHPuP8Cm+IBaHD+mH7BbjS/7BtixxuP6jq9tMn7S2SbOuXL8verDga2Ad+v4mtrvuSH6xP98pJHfpy6pqFEkKynURcL1opktBz4HBgJ3Av0AnHMfAP8E3sIH1Z+AN2t87S3ANfEe90uTOH4d8UlmP+HDHOfcMvxKbm+u55L4Ou/Z0JM1s72Bc/GTAXdp6NcnqVE1imQzc05XqkQkeGa2MTAbGIqfvDfDOXdvuFWJRItG6iKSLtfjryDcg2+XO8rMNgq3JJFo0UhdRAIXnwD4ClDonFsYXwp2ArDEOdcl3OpEokOhLiIiEhG6/C4iIhIRCnUREZGIUKiLiIhERNYtE9u6dWtXUFAQdhkiIiKBWbECPvsMVq6ErbeGb76Z9Z1zrk19X5d1oV5QUEBZWVnYZYiIiKScc1BaChddBJttBqNGwWGHgZktTObrdfldREQkA3z/PfToAf37w0EHwbx5PtAbQqEuIiISsqlTobAQXnwR7rgDxo3zl90bSqEuIiISkrVr4YYboFs32GgjmD4dLrkEmmxgOmfdPfVEVq9eTXl5OStXrgy7lMC1aNGCdu3a0bx587BLERGRRvj8czj1VJgyxf95333QqlXjvmckQr28vJxWrVpRUFCAWV1bR2c/5xxLliyhvLycDh06hF2OiIhsoFGj4PTTYfVqeOwx6N07Nd83EpffV65cyVZbbRXpQAcwM7baaqucuCIhIhJFK1bAeefB8cfDjjvC7NmpC3SISKgDkQ/0arlyniIiUfP++7D33v4y+6WX+vvnHTum9j0iE+pha9q0KYWFhb98LFiwgLKyMi644AIAJk2axPTp00OuUkRE0q2693yvveDbb2H8eLj9dj8xLtUicU89E2yyySbMmTNnnecKCgooKioCfKi3bNmSrl27hlGeiIiE4Icf4Kyz4Lnn4NBDYfhw2Gab4N5PI/UATZo0iaOPPpoFCxYwdOhQhgwZQmFhIVOnTg27NBERCdi0adCpE4weDYMHw8svBxvoEMGR+kUXQa0Bc6MVFsK//rX+Y1asWEFhYSEAHTp04IUXXvjltYKCAvr370/Lli259NJLU1uciIhklLVrYdAguP566NDB3zvfa6/0vHfkQj0siS6/i4hIbikvh+Ji33teXOwnxW22WfreP3KhXt+IWkREJAijR/ve88rK1PaeN4TuqadJq1atWL58edhliIhIiq1YAQMGwHHHQUEBvPtuOIEOCvW0OeaYY3jhhRc0UU5EJEI++AD22Qfuvdev2f7WW6nvPW+IyF1+D8tPP/30m+e6detGt27dANhpp52YN29emqsSEZEgOAcPPugnZ7ds6XdVO+KIsKvSSF1ERKRBfvgBTjoJzj4b9t/f73ueCYEOCnUREZGkvfmmb3MeNcr3no8fH3zveUMo1EVEROqxdi3cdBMcdBA0a+bD/bLLNnzf86DonrqIiMh6lJf7/c4nT4ZTToH7709v73lDKNRFRETqULP3fPhw36qWyZtlZtiFgzTq1s1/iIiI1LJyJZx/vu89z8/3+5736ZPZgQ65HOopVr316q677kqnTp248847qaqqWu/XLFiwgCeffDJNFYqISDI++MDve37PPXDxxb73fKedwq4qObkZ6rEYzJjhb5AUFPjHjVS99vv777/PxIkTGTduHNdff/16v0ahLiKSOap7z4uK4Ouv4aWX4J//hI03Druy5OVeqMdiUFLib5AALFzoH6cg2Ku1bduW0tJS7rnnHpxzLFiwgAMOOIDOnTvTuXNnpk+fDsCVV17J1KlTKSwsZMiQIXUeJyIiwfrxRzj5ZB8H++0Hc+fCkUeGXVXD5c5Euer75zNm/Bro1Soq4Iwz/P+iTZqUkrfbcccdqaqq4ttvv6Vt27ZMnDiRFi1a8Mknn9CrVy/Kysq49dZbueOOOxg7dmy8jIqEx4mISHCmT/ez2r/4Am67DS69NPNa1ZKVO6FerXag1/d8IzjnAFi9ejUDBgxgzpw5NG3alI8//jjh8ckeJyIijbd2Ldx6K1x3nZ8M9+ab/l56NsudUK8egRcU+EvuteXnp2yUDvDZZ5/RtGlT2rZty/XXX8/WW2/N3LlzqaqqokWLFgm/ZsiQIUkdJyIijfPFF773fNIk6NULhg7N3N7zhsjSCwyNMGgQ5OWt+1xenn8+RRYvXkz//v0ZMGAAZsbSpUvZdtttadKkCY8//jhr164Ffrsda13HiYhI6owZA7vvDjNnwqOP+ilVUQh0yMVQLy6G0tJfpzPm5/vHxcWN+rYrVqz4paXtkEMO4dBDD+W6664D4Nxzz2X48OHsu+++fPzxx2y66aYA7L777jRr1oxOnToxZMiQOo8TEZHGq+49797d/+qfNQtOOy3ze88bwqrv+2aLoqIiV3vy2Pz589lll10a9o2qJ86l8JJ7umzQ+YqI5LD586FnT7+j2t/+Brfckl2tamY2yzlXVN9xuXNPvbYsDHMREWkY5+Dhh+GCC2DTTX3veTa2qiUr9y6/i4hITvjxRz86P/NM6NIle3vPG0KhLiIikTN9ut/3/PnnfdvaK6/AttuGXVXwIhPq2TY3YEPlynmKiGyItWt9M9OBB/oFZKZNgyuuyN7FZBoqEqfZokULlixZEvnAc86xZMkS9a+LiCTwxRfwl7/ANdfAiSfCu+/CPvuEXVV6RWKiXLt27SgvL2fx4sVhlxK4Fi1a0K5du7DLEBHJKC++CP36wYoV8Mgj0WtVS1YkQr158+Z06NAh7DJERCTNVq6Eyy+Hu+/299BHjICddw67qvAEdvndzHYwszfMbL6ZvW9mFyY4ppuZLTWzOfGPvwdVj4iIRMuHH8K++/pAv+giv19XLgc6BDtSXwNc4pybbWatgFlmNtE590Gt46Y6544OsA4REYmQmr3neXkwdiwcdVTYVWWGwEbqzrmvnHOz458vB+YD2wf1fiIiEn0//ug3YKnuPZ83T4FeU1pmv5tZAbAH8HaCl7uY2Vwze9nMdq3j60vMrMzMynJhMpyIiPzWW2/BHnvAyJF+mddc6T1viMBD3cxaAs8BFznnltV6eTaQ75zrBNwNjEr0PZxzpc65IudcUZs2bYItWEREMsratXDzzXDAAf7xtGlw5ZW503veEIH+JzGz5vhAjznnnq/9unNumXPup/jn44DmZtY6yJpERCR7fPml7z0fOBB69IA5c/zkOEksyNnvBjwEzHfO3VnHMdvEj8PM9o7XsySomkREJHuMHev3PX/7bT8x7qmn4He/C7uqzBbk7Pf9gN7Af8xsTvy5q4H2AM65oUAP4BwzWwOsAHq6qC8LJyIi61VZ6XvP77pLvecNFVioO+emAetdz8c5dw9wT1A1iIhIdvnwQz+7fc4cuPBCuO227Nr3PGyRWFFORESym3N+edfzz/e95y++CEdrBZMG09xBEREJ1dKlfnR+xhl+A5a5cxXoG0qhLiIi6RWLQUEBNGlC5bYF/P0PMUaO9G1rEyfCdtuFXWD20uV3ERFJn1gMSkqgogKAjb9eyC1WwoDroONVxSEXl/00UhcRkfQZOPCXQK+W5yro+MjAkAqKFoW6iIikjVu0KPELdT0vDaJQFxGRwFVW+u1RF7r2iQ9oX8fz0iAKdRERCdRHH/mlXf/9b5hy2CBcXt66B+TlwaBB4RQXMQp1EREJRHXveefO8PnnMGYM9BlfjJWWQn4+mPk/S0uhWJPkUkGz30VEJOWWLoX+/f0Sr3/+Mzz+OGy/ffzF4mKFeEA0UhcRkZR6+22/7/mzz/qr6hMn1gh0CZRCXUREUqKqCm69Ffbf338+dSpcfTU0bRp2ZblDl99FRKTRvvoKeveG116Dk06CBx6AzTcPu6rco1AXEZFGGTcOTjsNfv4Zhg2D00/3c+Ak/XT5XURENkhlJfztb3DUUX699lmz/KYsCvTwaKQuIiIN9tFH0LOn3/f8/PNh8GBo0SLsqkShLiIiSXMOhg+HAQN8iI8ZA8ccE3ZVUk2X30VEJCnLlvn28n79YK+9/L7nCvTMolAXEZF6vf02FBbCM8/ATTfBq6+q9zwTKdRFRKROVVVw222/9p5PmeJ3T1XveWbSPXUREUnoq6+gTx8/Kj/xRL9Eu3rPM5tCXUREfqNm7/mDD6pVLVvo8ruIiPyishIuvnjd3vMzz1SgZwuN1EVEBICPP4ZevWD2bN+ydvvt6j3PNgp1EZEc5xw89hicdx5svDGMHg3HHht2VbIhdPldRCSHLVsGp54KfftCURHMm6dAz2YKdRGRXBSLUbltAS1/14Sbnyzg+R4xXntNvefZTpffRURyTNUTMdaeXsLGqysAyGch+eNKYAR+yTjJWhqpi4jkkK++gsVnDaR5PNB/UVHhV5WRrKZQFxHJES+/DJ06QZuVixIfsKiO5yVrKNRFRCKushIuuQSOPBK22QbWbNs+8YHt63hesoZCXUQkwj7+GLp2hTvv9L3n77wDG90+CPLy1j0wLw8GDQqnSEkZhbqISARV95537gwLFsCoUXD33fHFZIqL/ULu+fl+qbj8fP9Yk+Synma/i4hEzLJlcO65EIvBQQfBE09Au3a1DiouVohHkEbqIiIR8s47sMce8NRTcMMN8NprCQJdIkuhLiISAVVVMHgw7LcfrFnj9z2/9lrte55rdPldRCTLff213/d84kQ44QS/VeoWW4RdlYRBI3URkSw2frzvPZ82zc91e/ZZBXouU6iLiGShVavg0kvhiCOgbVsoK4OzztK+57lOl99FRLLMJ5/4fc9nzfKz3O+4AzbZJOyqJBMo1EVEsshjj/kg32gjeOEFOO64sCuSTKLL7yIiWWDZMujdG047DfbcE+bOVaDLbynURUQy3MyZfmW4J5/0veevvw477BB2VZKJFOoiIhmqqsrfL+/a1U+MmzxZveeyfrqnLiKSgb7+2l9qf+UV+OtfYdgwtapJ/TRSFxHJMBMm+N7zKVNg6FAYOVKBLslRqIuIhCUWg4ICaNIECgpYMzzGZZfB4Yf/2nt+9tnqPZfkBRbqZraDmb1hZvPN7H0zuzDBMWZmd5nZp2Y2z8w6B1WPiEhGicWgpAQWLvT7pC5cyOp+JXxxR4xzz/Ubs+y6a9hFSrYJ8p76GuAS59xsM2sFzDKzic65D2occwTQMf6xD3B//E8RkWgbOBAqKtZ5ahNXwbA2A8m7V1uiyoYJbKTunPvKOTc7/vlyYD6wfa3DugOPOW8GsLmZbRtUTSIiGWPRooRP532X+HmRZKTlnrqZFQB7AG/Xeml74PMaj8v5bfCLiERO5dbtE7/Qvo7nRZIQeKibWUvgOeAi59yy2i8n+BKX4HuUmFmZmZUtXrw4iDJFRNKiqgr++U8489tBVFjeui/m5cGgQeEUJpEQaKibWXN8oMecc88nOKQcqLkuUjvgy9oHOedKnXNFzrmiNm3aBFOsiEjAvvkGjjzS7672c/diqoaWQn6+n96en+/3Ti3W/XTZcIFNlDMzAx4C5jvn7qzjsDHAADMbgZ8gt9Q591VQNYmIhGXCBOjTx6/hfv/91a1qxVCiEJfUCXL2+35Ab+A/ZjYn/tzVQHsA59xQYBxwJPApUAH0C7AeEZG0W7XKT3S/4w7fovbaa7DbbmFXJVEVWKg756aR+J55zWMccF5QNYiIhOnTT/2+52VlcM45/l669j2XIGntdxGRAMRi0L8/NG8Ozz8Pxx8fdkWSC7RMrIhICi1f7jdiOfVUKCyEOXMU6JI+CnURkRSZNcvve/7EE3DddfDGG2o7l/RSqIuINFJ173mXLrBypQ/zf/wDmukGp6SZ/sqJiDTCN99A374wfry/zD5sGGy5ZdhVSa7SSF1EZAO98orf93zSJN97/txzCnQJl0JdRKSBVq2CK66Aww6D1q1h5kw/0137nkvYdPldRKQB/vtf33s+c6ZfFe7OO/2S7SKZQKEuIpKkWMwvItO0KYwcCSecEHZFIuvS5XcRkXrU7D3v1AnmzlWgS2ZSqIuIrMesWbDnnuo9l+ygUBcRSaCqyt8v79IFVqyA119X77lkPoW6iOS2WAwKCqBJE/9nLMa338JRR8Ell/g/58yBgw4Ku1CR+un/OUUkd8ViUFICFRX+8cKFrD2jhOs2hjcqi7nvPrWqSXYxv/tp9igqKnJlZWVhlyEiUVBQAAsX/ubpL5rl8/3sBfzpT+kvSSQRM5vlnCuq7ziN1EUkdy1alPDp7dYuYnsFumQh3VMXkdxVxzR20/R2yVIKdRHJST/9BEPbD+Jnai0Hl5cHgwaFU5RIIynURSTnzJ7t9z0/781iXj6+FNc+38+Gy8+H0lIoLg67RJENonvqIpIznIN//ctvxtK2re89P+igYkAhLtGgUBeRnPDtt9CvH4wbB927w0MPwVZbhV2VSGrp8ruIRN6rr/o12197De65B154QYEu0aRQF5HIWr0arrwSDj0UttgC3nkHzjtPi8lIdOnyu4hE0mef+X3P33nHLxo3ZIj2PZfoU6iLSOQ89RScfbbf9/zZZ6FHj7ArEkkPXX4Xkcj46Sc/Ge6UU2D33f1GLAp0ySUKdRGJhHff9fueDx8O114Lkyb5tnORXKJQF5GsVt17vu++fqT++utwww3a91xyk/7ai0jWWrwY+vb1vefHHut7z1u3DrsqkfBopC4iWem11/x98+re81GjFOgiCnURySqrV8NVV8Ff/qLec5HadPldRLLG//7ne8/ffhvOOsv3nm+6adhViWQOhbqIZIURI3zvuRk88wyceGLYFYlkHl1+F5GM9vPPcPrpfoS+224wd64CXaQuCnURyVjVveePPgrXXAOTJ6v3XGR9FOoiknGcg3//2/eeL1/uZ7jfeKN6z0Xqo38iIpJRFi/2S72+9BIccww8/LBa1USSpZG6iIQvFoOCAlyTJqzcpoCtxse4+24YPVqBLtIQGqmLSLhiMVxJCVZRgQE7uIU80qKEJlsAVhx2dSJZRSN1EQnV6ssHYhUV6zzXZGUFDBwYUkUi2UuhLiKhefppaPrlosQvLqrjeRGpk0JdRNLu55/hjDOgZ0/4ZqP2iQ9qX8fzIlInhbqIpNWcOb73/JFH/BX2Ng8Ogry8dQ/Ky4NBg8IpUCSLKdRFJC2cg7vugn32+bX3/KaboFmfYigt9avKmPk/S0uhWJPkRBpKs99FJHD19p4XFyvERVJAI3URCdTrr0OnTjBxoh+pq/dcJDgKdREJxOrV/p75IYfAZpv57VLPP1/7nosESZffRSTl/vc/OOUUmDHDz3L/97+177lIOgQ2Ujezh83sWzN7r47Xu5nZUjObE//4e1C1iEj6PP00FBbCBx/4PdCHDVOgi6RLkJffHwUOr+eYqc65wvjHDQHWIiIB+/lnOPNM33v+xz/61rWTTw67KpHcElioO+emAN8H9f1FJHNU954//DBcfTVMmQIdOoRdlUjuCXuiXBczm2tmL5vZrnUdZGYlZlZmZmWLFy9OZ30ish41e8+XLYNXX/VrxjRvHnZlIrkpzFCfDeQ75zoBdwOj6jrQOVfqnCtyzhW1adMmbQWKSN2++w66d4cLL4RDD4V58+Dgg8OuSiS3hRbqzrllzrmf4p+PA5qbmbpXRbLAG2/43vMJE/zM9jFj1HsukglCC3Uz28bMd6ya2d7xWpaEVY+I1G/1arjmGvi//4NWrXzv+QUXqPdcJFME1qduZk8B3YDWZlYOXAc0B3DODQV6AOeY2RpgBdDTOeeCqkdEGmfBAt97/tZbcPrp/l66WtVEMku9oW5mA4CYc+6Hhnxj51yvel6/B7inId9TRMLxzDNQUuInxj31lG9bE5HMk8zl922AmWb2jJkdXn3JXESi7+ef4ayzfL/5Lrv41jUFukjmqjfUnXPXAB2Bh4C+wCdmdrOZ/T7g2kQkRHPnQlERPPQQXHWVes9FskFSE+Xi97q/jn+sAbYARprZ4ABrE5F0isWgoADXpAnLtizgzj1jLF3qd1e7+Wb1notkg2TuqV8AnAZ8BwwDLnPOrTazJsAnwOXBligigYvF/E3zigoM2OyHhTzQpIRVf4fN/k/7nItki2RG6q2BvzrnDnPOPeucWw3gnKsCjg60OhFJj4EDoaJinadaVFWw2a0DQypIRDZEvSN151ydu6c55+anthwRSbc1a6DpwkUknAG7aFG6yxGRRgh77XcRCdGCBXDggbCQ9okPaF/H8yKSkRTqIjnq2Wf9vufvvQdfnjsI8vLWPSAvz+/OIiJZQ6EukmMqKvycuJNOgp139r3nXe8thtJSyM/3a77m5/vHxZokJ5JNAlsmVkQyz7x5fvGYDz+EK6+EG26o0apWXKwQF8lyGqmL5ADn4N57Ye+94Ycf4JVX4JZb1HsuEjUaqYtE3JIlfgOWMWPgyCPh0UehTZuwqxKRIGikLhJhkyb5fc9ffhmGDIGxYxXoIlGmUBeJoDVr4Npr4eCD/faoM2bARRdp33ORqNPld5GIWbjQ73s+fTr07Qt33w0tW4ZdlYikg0JdJEJGjvRbpa5d65dzP+WUsCsSkXTS5XeRCKiogLPPhhNPhJ128r3nCnSR3KNQF8ly8+b5fc9LS+GKK2DaNNhxx7CrEpEwKNRFslSi3vNbb1XvuUgu0z11kSy0ZAmccQaMHg1HHOF7z9u2DbsqEQmbRuoiWWbyZN97Pm4c3Hmn7z1XoIsIKNRFssaaNXDddb73PC/P957/7W/QRP+KRSROvw5EMkksBgUFPqkLCvxjfO95t25+A5bevWH2bOjcOcxCRSQT6Z66SKaIxfyeqBUV/vHChVBSwtvvwOGPFbN2LTzxhDZSE5G6aaQukikGDvw10KtVVLD1XQPp2BHefVeBLiLrp5G6SKZYtCjh0/ksYto02GijNNcjIllHI3WRTNG+fcKnLb+9Al1EkqJQF8kQP109iJVN89Z9Mi8PBg0KpyARyToKdZEMMGUK7HJjMWe5UpZtkY8zg/x8v/arbqSLSJIU6iIhqu49//OfoUULuOidYjb7fgFWVQULFijQRaRBNFFOJCSLFvnMnjYNTjvN73veqlXYVYlINlOoi4TguefgzDP9SF295yKSKrr8LpJGFRXQvz/06AEdO/p9zxXoIpIqCnWRNHnvPb9N6gMPwGWX+cvuv/992FWJSJQo1EUC5hzcfz/stRd89x1MmACDB2sxGRFJPYW6SIC+/x5OOAHOPRcOOgjmzoVDDw27KhGJKoW6SECmTPH7no8dC3fc4fc/33rrsKsSkShTqIuk2Jo18I9//Np7Pn06XHKJ9j0XkeCppU0khRYtglNPhalToU8fuOce9Z6LSPoo1EVS5Pnnfe/56tXw+OM+3EVE0kkXBEUaacUKOOccPyHu97/3+54r0EUkDAp1kUZ47z3fqjZ0KFx6Kbz5JvzhD2FXJSK5SqEusgGc80G+116weDGMHw+3367ecxEJl0JdpIG+/94v83rOOb73fN48OOywsKsSEVGoizTI1KlQWAgvvqjecxHJPAp1kSSsWQPXXw/duvlL7Oo9F5FMpF9JIrXFYlBQ4BO7oIAld8U4+GC/oMwpp/jZ7UVFYRcpIvJbgYW6mT1sZt+a2Xt1vG5mdpeZfWpm88ysc1C1iCQtFoOSEli40M+GW7iQTS4s4Q/vxHjsMd9/rsVkRCRTBTlSfxQ4fD2vHwF0jH+UAPcHWItIcgYO9Jue15BHBQ9sNZDevUOqSUQkSYGFunNuCvD9eg7pDjzmvBnA5ma2bVD1iCRl0aKETzf/KvHzIiKZJMx76tsDn9d4XB5/TiQUzsHyLdonfrF9Hc+LiGSQMEPdEjznEh5oVmJmZWb37RILAAARmUlEQVRWtnjx4oDLklxU3Xt+9veDWNkkb90X8/Jg0KBwChMRaYAwQ70c2KHG43bAl4kOdM6VOueKnHNFbdq0SUtxkjuqe8/HjIE9Bhez0fBSyM8HM/9naSkUF4ddpohIvcLcpW0MMMDMRgD7AEudc1+FWI/kmLVr/QD8+uuhQwffe77XXgDFcKpCXESyT2ChbmZPAd2A1mZWDlwHNAdwzg0FxgFHAp8CFUC/oGoRqa283A++p0zxf953H2y2WdhViYg0TmCh7pzrVc/rDjgvqPcXqcuoUXDGGVBZCY89hlrVRCQytKKc5IwVK+C88+D44/2Cce++q0AXkWhRqEtOeP992Htvf5n9kkvgrbegY8ewqxIRSa0wJ8qJBM45ePBBuOgiaNnS76p2xBFhVyUiEgyN1CWyfvgBTjwRzj4b9t/f73uuQBeRKFOoSyS9+abvPR89GgYPhvHjYZttwq5KRCRYCnWJlLVr4cYb4cADoVkzH+6XXaZ9z0UkN+ieukRGeTmceipMnuz3Pb//fvWei0huUahLJIweDaef7nvPhw/3rWqWaHcBEZEI00VJyWorVsCAAXDccb73fPZs6NNHgS4iuUmhLlnrgw9gn33g3nvh4ov92u077RR2VSIi4dHld8k6tXvPX3oJjjwy7KpERMKnkbpklR9+gJNO8r3n++0Hc+cq0EVEqinUJWtMn+57z0eNgttugwkTYNttw65KRCRzKNQlM8VifuZbkya4/AJeODG2Tu/55Zer91xEpDbdU5fME4tBSQlUVABgixZy6KIS/tUF+owvVu+5iEgdNNaRzDNw4C+BXm1TKjjvi4EKdBGR9VCoS8ZxixYlfN4+T/y8iIh4CnXJKPPnw1fN2id+sX0dz4uICKBQlwzhHAwbBnvuCTe0GMSajfPWPSAvDwYNCqc4EZEsoVCX0P34I5x8Mpx1FnTtCtd9VEyzh0ohP9+v95qfD6WlUFwcdqkiIhlNs98lVNOn+x3VvvgCbr21xjapxcUKcRGRBtJIXUKxdq2/mn7ggT7Ep02DK65Q77mISGNopC5p98UXft/zSZOgZ08YOhR+97uwqxIRyX4KdUmrMWOgXz9YuRIeeQROO03bpIqIpIoudkparFwJF1wA3bv7zrTZs6FvXwW6iEgqKdQlcPPn+33P777bb5c6YwbsvHPYVYmIRI8uv0tgnIOHH/Yj9Lw8GDsWjjoq7KpERKJLI3UJxI8/+klwZ54JXbrAvHkKdBGRoCnUJeWq9z1/7jm45RZ45RXtey4ikg4KdUmZmr3nZr73/Mor1XsuIpIuuqcuKfHFF9C7N7zxhl/y9YEH1HsuIpJuCnVptLFjfXvaihV+Ypxa1UREwqELo7LBVq6ECy+EY46BHXbwvef9+inQRUTColCXDfLhh7DvvnDXXT7Y1XsuIhI+XX6XBqnde/7ii3D00WFXJSIioJG6NMCPP0KvXr73fJ99YO5cBbqISCZRqEtS3noL9tgDRo6Em2+GiRNhu+3CrkpERGpSqEvdYjFcfgHOmrBd1wKO/SnGtGlw1VXQtGnYxYmISG26py6JxWJUnVVCkxUVAOSzkH/9XIL9F9i3ONzaREQkIY3UJaGKvw38JdCr2YoKGDgwpIpERKQ+CnVZR2Wlb1FrsXhR4gMW1fG8iIiETqEuv6jZe/5jq/aJD2pfx/MiIhI6hbr80nu+557w+ecwZgxsef8g34heU16e37FFREQykkI9xy1dCqecAmec8Wvv+THHAMXFUFoK+fl+3df8fP+4WJPkREQylWa/57AZM/xiMp9/7gfgV1xRq1WtuFghLiKSRTRSz0FVVXDLLbD//v7S+9SpcPXV6j0XEcl2GqnnmC+/hD594LXX4KST/L7nm28edlUiIpIKCvUc8tJLfq/zn3+GYcPg9NO1TaqISJQEevndzA43s4/M7FMzuzLB633NbLGZzYl/nBlkPbmqshIuushvvrLddjBrlp8Yp0AXEYmWwEbqZtYUuBf4C1AOzDSzMc65D2od+rRzbkBQdeS6jz6Cnj1hzhw4/3wYPBhatAi7KhERCUKQI/W9gU+dc58551YBI4DuAb6f1OAcPPIIdO78a+/5XXcp0EVEoizIUN8e+LzG4/L4c7WdYGbzzGykme2Q6BuZWYmZlZlZ2eLFi4OoNVKWLvWdaKefDnvvXaP3XEREIi3IUE90x9bVevwiUOCc2x14FRie6Bs550qdc0XOuaI2bdqkuMxoefttv+/5M8/ATTfBq6/C9on+V0pERCInyFAvB2qOvNsBX9Y8wDm3xDlXGX/4ILBngPVEWlUV3Hqr7z2vqoIpU/yGauo9FxHJHUGG+kygo5l1MLONgJ7AmJoHmNm2NR4eC8wPsJ7I+uorOPRQuOoqOP54Pymua9ewqxIRkXQLbPa7c26NmQ0AJgBNgYedc++b2Q1AmXNuDHCBmR0LrAG+B/oGVU9UjRsHp53me88ffFCtaiIiucycq32bO7MVFRW5srKysMsIXWWlH5kPGQK77w4jRsAuu4RdlYiIBMHMZjnniuo7TivKZaGPP/a95+++CwMGwO23q1VNREQU6lnFORg+3Af5xhvD6NFw7LFhVyUiIplCu7RliWXLfO95v35QVATz5inQRURkXQr1LPD221BY6HvPb7zR77Cm3nMREalNoZ7Bqqrgttt87/natTB5MlxzjXrPRUQkMd1Tz1BffeX3PX/1VejRA0pLYYstwq5KREQymUI9A738su89/+knH+ZnnqnecxERqZ8uv2eQykq4+GI48kjYZhsoK4OzzlKgi4hIcjRSzxAffwy9esHs2eo9FxGRDaNQD5lz8NhjcN55vvd81Cjorl3nRURkA+jye4iWLYNTT4W+fX3v+dy5CnQREdlwCvWQvPOO3/d8xAi44Qbfe96uXdhViYhINlOop1lVFQweDPvtB2vW+H3Pr71WveciItJ4uqeeRl9/7XvPJ06EE07wW6Wq91xERFJFI/U0efllv0XqtGm+9/zZZxXoIiKSWgr1gFVWwiWX+N7zrbdW77mIiARHl98DVLP3/Nxz4Y47YJNNwq5KRESiSqEeAOfg8cd9kG+0EbzwAhx3XNhViYhI1Onye4otWwa9e/u12/fc0/eeK9BFRCQdFOopNHMmdO4MTz3le89ffx122CHsqkREJFco1FOgqsqv1d61K6xa5fc9V++5iIikm+6pN1LN3vO//hWGDVOrmoiIhEMj9UYYPx46dYKpU2HoUBg5UoEuIiLhUahvgFWr4NJL4YgjoG1b33t+9tnqPRcRkXDp8nsDffKJ7z2fNUu95yIiklkU6g1Q3XvevDk8/zwcf3zYFYmIiPxKl9+TsHy57z3v08dvlzp3rgJdREQyj0K9HmVlPsiffBKuvx7eeEO95yIikpkU6nWoqvL3y7t08RPjJk2Cv/9dveciIpK5dE89ga+/9su8vvKKv8w+bBhsuWXYVYmIiKyfRuq1TJjge8+nTIH774fnnlOgi4hIdlCox61aBZddBocfDm3a+HXc+/dX77mIiGQPXX4HPv3U956XlcE558A//6necxERyT45H+pPPOGDXL3nIiKS7XL28vvy5b7vvHdvKCyEOXMU6CIikt1yMtTLyvy+57EYXHed7z1v3z7sqkRERBonp0K9qsrfL+/aFVau9GH+j39As5y/CSEiIlGQM3H2zTe+93zCBPWei4hINOXESP2VV3zv+eTJ6j0XEZHoinSor1oFl18Ohx0GrVur91xERKItspffa/aen3023Hkn5OWFXZWIiEhwIhnqsZgfkTdrBiNHwgknhF2RiIhI8CJ1+X35cj8Z7tRTfe/53LkKdBERyR2RCfVZs3zv+RNPqPdcRERyU9aHenXveZcuvvf89dfVey4iIrkpq6Pvm2+gb18YPx6OO873nm+1VdhViYiIhCNrR+oTJ/re8zfegPvu85uxKNBFRCSXBRrqZna4mX1kZp+a2ZUJXt/YzJ6Ov/62mRXU9z2dgyuugEMP9SE+c6bfZU295yIikusCC3UzawrcCxwB/BHoZWZ/rHXYGcAPzrk/AEOA2+r7vh9+CIMH+97zmTPhT39KdeUiIiLZKciR+t7Ap865z5xzq4ARQPdax3QHhsc/Hwn8n9n6x9yVlfDsszB0qBaTERERqSnIUN8e+LzG4/L4cwmPcc6tAZYC670z/sc/Qo8eKaxSREQkIoKc/Z5oxO024BjMrAQoiT+sNLP3GllbNmoNfBd2ESHQeecWnXdu0XknLz+Zg4IM9XJghxqP2wFf1nFMuZk1A34HfF/7GznnSoFSADMrc84VBVJxBtN55xadd27ReeeWIM87yMvvM4GOZtbBzDYCegJjah0zBjgt/nkP4HXn3G9G6iIiIlK/wEbqzrk1ZjYAmAA0BR52zr1vZjcAZc65McBDwONm9il+hN4zqHpERESiLtAV5Zxz44BxtZ77e43PVwInNvDblqagtGyk884tOu/covPOLYGdt+lqt4iISDRk7TKxIiIisq6MDfUglpjNBkmc94FmNtvM1phZZDr2kzjvi83sAzObZ2avmVlS7R2ZLonz7m9m/zGzOWY2LcGqjFmpvvOucVwPM3NmFokZ0kn8vPua2eL4z3uOmZ0ZRp2plszP28xOiv8bf9/Mnkx3jUFI4uc9pMbP+mMz+7HRb+qcy7gP/MS6/wI7AhsBc4E/1jrmXGBo/POewNNh152m8y4AdgceA3qEXXMaz/vPQF7883Ny6Oe9WY3PjwXGh113Os47flwrYAowAygKu+40/bz7AveEXWsI590ReBfYIv64bdh1p+O8ax1/Pn5CeaPeN1NH6oEsMZsF6j1v59wC59w8oCqMAgOSzHm/4ZyriD+cgV/3INslc97LajzclASLM2WhZP59A9wIDAZWprO4ACV73lGTzHmfBdzrnPsBwDn3bZprDEJDf969gKca+6aZGuqBLDGbBZI57yhq6HmfAbwcaEXpkdR5m9l5ZvZffMBdkKbaglTveZvZHsAOzrmx6SwsYMn+PT8hfptppJntkOD1bJPMee8E7GRmb5rZDDM7PG3VBSfp32vx24kdgNcb+6aZGuopW2I2y0TxnJKR9Hmb2alAEXB7oBWlR1Ln7Zy71zn3e+AK4JrAqwrees/bzJrgd228JG0VpUcyP+8XgQLn3O7Aq/x6NTKbJXPezfCX4LvhR6zDzGzzgOsKWkN+n/cERjrn1jb2TTM11BuyxCzrW2I2yyRz3lGU1Hmb2SHAQOBY51xlmmoLUkN/3iOA4wKtKD3qO+9WwG7AJDNbAOwLjInAZLl6f97OuSU1/m4/COyZptqClOzv89HOudXOuf8BH+FDPps15N93T1Jw6R0yN9RzdYnZZM47iuo97/jl2AfwgR6F+22Q3HnX/MV2FPBJGusLynrP2zm31DnX2jlX4JwrwM+hONY5VxZOuSmTzM972xoPjwXmp7G+oCTze20UfjIsZtYafzn+s7RWmXpJ/T43s52BLYC3UvGmGRnq8Xvk1UvMzgeecfElZs3s2PhhDwFbxZeYvRiosy0mWyRz3ma2l5mV41fie8DM3g+v4tRI8ud9O9ASeDbe/pH1/7OT5HkPiLf4zMH/PT+tjm+XNZI878hJ8rwviP+85+LnT/QNp9rUSfK8JwBLzOwD4A3gMufcknAqTo0G/D3vBYxI1aBUK8qJiIhEREaO1EVERKThFOoiIiIRoVAXERGJCIW6iIhIRCjURUREIkKhLiIiEhEKdRERkYhQqIvIesUXPJpnZi3MbNP44ii7hV2XiPyWFp8RkXqZ2U1AC2AToNw5d0vIJYlIAgp1EalXfO3qmfi9zbumYjcpEUk9XX4XkWRsiV97vxV+xC4iGUgjdRGpV3wDnRFAB2Bb59yAkEsSkQSahV2AiGQ2M+sDrHHOPWlmTYHpZnawc+71sGsTkXVppC4iIhIRuqcuIiISEQp1ERGRiFCoi4iIRIRCXUREJCIU6iIiIhGhUBcREYkIhbqIiEhEKNRFREQi4v8Dco6IwdjRRBkAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "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 = 1)\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", " xPlot = np.zeros(2)\n", " xPlot[1] = 1.1*Tmean[nLen - 1]**2\n", " plt.plot(xPlot, fitLine(pFinal, xPlot), color = 'b',\\\n", " linestyle = '-', label = \"Fit\")\n", " plt.errorbar(lengths, Tmean**2, xerr = lenError, yerr = T2error, \n", " color ='r', marker = 'o', 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": [ "Calculate value of $g$ and estimate its error." ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " \n", "Acceleration due to gravity 9.89 +/- 0.08 m/s^2\n" ] } ], "source": [ "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": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.6.6" } }, "nbformat": 4, "nbformat_minor": 2 }