Exercise 3 - Monte-Carlo Error Propagation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Practice Problem - Monte-Carlo Error Propagation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Part 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "You have likely encountered the concept of propagation of uncertainty before (see [the usual rules here](http://en.wikipedia.org/wiki/Propagation_of_uncertainty#Example_formulas)). The idea is that given measurements with uncertainties, we can find the uncertainty on the final result of an equation.\n",
    "\n",
    "For example, let us consider the following equation:\n",
    "\n",
    "$$F = \\frac{G~M_1~M_2}{r^2}$$\n",
    "\n",
    "which gives the gravitational force between two masses $M_1$ and $M_2$ separated by a distance $r$.\n",
    "\n",
    "Let us now imagine that we have two masses:\n",
    "\n",
    "$$M_1=40\\times10^4\\pm0.05\\times10^4\\rm{kg}$$\n",
    "\n",
    "and\n",
    "\n",
    "$$M_2=30\\times10^4\\pm0.1\\times10^4\\rm{kg}$$\n",
    "\n",
    "separated by a distance:\n",
    "\n",
    "$$r=3.2\\pm0.01~\\rm{m}$$\n",
    "\n",
    "where the uncertaintes are the standard deviations of Gaussian distributions which could be e.g. measurement errors.\n",
    "\n",
    "We also know:\n",
    "\n",
    "$$G = 6.67384\\times10^{-11}~\\rm{m}^3~\\rm{kg}^{-1}~\\rm{s}^{-2}$$\n",
    "\n",
    "(exact value, no uncertainty)\n",
    "\n",
    "Use the [standard error propagation rules](http://en.wikipedia.org/wiki/Propagation_of_uncertainty#Example_formulas) to determine the resulting force and uncertainty in your script (you can just derive the equation by hand and implement it in a single line in your code).\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Resulting Force F = 7.82e-01 +- 5.63e-03\n"
     ]
    }
   ],
   "source": [
    "# Computing with standard error propagation rules: \n",
    "M1 = 40e4\n",
    "dM1 = 0.05e4\n",
    "M2 = 30e4\n",
    "dM2 = 0.1e4\n",
    "r = 3.2\n",
    "dr = 0.01\n",
    "G = 6.67384e-11\n",
    "\n",
    "F_oldschool = (G*M1*M2)/r**2\n",
    "dF_oldschool = ( (dM1*G*M2/r**2)**2 + (dM2*G*M1/r**2)**2 + (-2*G*M1*M2*dr/r**3)**2 )**0.5\n",
    "\n",
    "print(\"Resulting Force F = {0:.2e} +- {1:.2e}\".format(F_oldschool, dF_oldschool))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Now, we can try using a **Monte-Carlo** technique instead. The idea behind Monte-Carlo techniques is to generate many possible solutions using random numbers and using these to look at the overall results. In the above case, you can propagate uncertainties with a Monte-Carlo method by doing the following:\n",
    "\n",
    "* randomly sample values of $M_1$, $M_2$, and $r$, 1000000 times, using the means and standard deviations given above\n",
    "\n",
    "* compute the gravitational force for each set of values\n",
    "\n",
    "You should do this with Numpy arrays, and **without any loops**. You should then get an array of 1000000 different values for the forces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# These are the parameters\n",
    "# d.. values are uncertainties, standard deviations \n",
    "M1 = 40e4\n",
    "dM1 = 0.05e4\n",
    "M2 = 30e4\n",
    "dM2 = 0.1e4\n",
    "r = 3.2\n",
    "dr = 0.01\n",
    "G = 6.67384e-11\n",
    "number = 1000000 #number of values we want to generate\n",
    "\n",
    "# Ranomly sampling values\n",
    "M1ar = np.random.normal(M1, dM1, number)\n",
    "M2ar = np.random.normal(M2, dM2, number)\n",
    "rar = np.random.normal(r, dr, number)\n",
    "Gar = np.ones(number)*G\n",
    "\n",
    "# Now we compute the force with the formula\n",
    "F = (M1ar*M2ar*Gar)/rar**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Make a plot of the normalized histogram of these values of the force, and then overplot a Gaussian function with the mean and standard deviation derived with the standard error propagation rules. Make sure that you pick the range of x values in the plot wisely, so that the two distributions can be seen. Make sure there are also a sensible number of bins in the histogram so that you can compare the shape of the histogram and the Gaussian function. The two distributions should agree pretty well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgsAAAFkCAYAAACuFXjcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xl4HNWZ7/Hva3nfhMF4A8vYgB2xYxE2G2IgQAhDEsgE\nouCQBJLJ3CzD+N6ZJMwkN2SZyQwQTBa4CRMCmEWEsIYlGGJ2DBgs78grtmVbi2kv8iLLsqRz/+iW\npaputdVSV1cvv0+efuC8dbr7PWkhvX3q1ClzziEiIiLSlT5hJyAiIiLZTcWCiIiIJKViQURERJJS\nsSAiIiJJqVgQERGRpFQsiIiISFIqFkRERCQpFQsiIiKSlIoFERERSUrFgoiIiCSVUrFgZn3M7Gdm\n9qGZNZrZWjP7YYJ+PzWzmlifl8zsuPSlLCIiIpmU6szCD4BvAt8CPgZ8D/iemX2nvYOZfR/4Tqzf\nmcBeYK6Z9U9LxiIiIpJRlsqNpMzsGaDOOfeNTrHHgEbn3HWxdg1wq3Nudqw9HKgHvuKcezSdyYuI\niEjwUp1ZmA9cZGbHA5jZqcA04PlYeyIwBpjX/gTn3C7gXeCcdCQsIiIimdU3xf7/BQwHVppZK9Fi\n49+dc4/Ejo8BHNGZhM7qY8fimNkRwKXABqApxXxEREQK2UDgGGCuc25bUG+SarFwDfAl4IvAB8Bp\nwK/MrMY590CS5xnRIiKRS4GHUsxDREREOlwLPBzUi6daLNwC/Kdz7s+x9gozOwa4CXgAqCNaGIzG\nO7swCljUxWtuAHjwwQcpLS1NMZ3cMmvWLGbPnh12GoErlHFC4YxV48wvGmf+qKqqYubMmRD7WxqU\nVIuFwcTPELQRW/vgnFtvZnXARcBSOLjA8Szgzi5eswmgtLSUqVOnpphObikuLs77MULhjBMKZ6wa\nZ37ROPNSoKfxUy0WngH+3cw2ASuAqcAs4A+d+twB/NDM1hKtdH4GbAae7nW2IiIiknGpFgvfIfrH\n/06ipxZqgP8XiwHgnLvFzAYDvwcOA94ALnPONaclYxEREcmolIoF59xe4H/HHsn63Qzc3OOsRERE\nJGvo3hAZVF5eHnYKGVEo44TCGavGmV80TklVSjs4BpKA2VRg4cKFCwtpIYqIiEivVVZWUlZWBlDm\nnKsM6n00syAiIiJJqVgQERGRpFQsiIiISFIqFkRERCQpFQsiIiKSlIoFERERSUrFgoiIiCSlYkFE\nRESSUrEgIiIiSalYEBERkaRULIiIiEhSKhZEREQkKRULIiIikpSKBREREUlKxYKIiIgkpWJBRERE\nklKxICIiIkmpWBAREZGkVCyIiIhIUioWREREJCkVCyIiIpKUigURERFJSsWCiIiIJKViQURERJJS\nsSAiIiJJ9Q07ARHJD/sO7OPN6jdZs30NA4oGUDaujFNHn4qZhZ2aiPRSSsWCma0HJiQ4dKdz7rtm\nNgC4HbgGGADMBb7lnNva60xFJCu1tLVw2/zbuHX+rWzft91z7JTRp3Dbxbdx8bEXh5SdiKRDqqch\nzgDGdHpcDDjg0djxO4DLgc8D5wPjgMfTkqmIZJ1IY4QZ983gpnk3xRUKAEvrl3LJg5fwo5d/hHMu\nhAxFJB1Smllwzm3r3DazK4B1zrk3zGw4cD3wRefca7HjXwOqzOxM59yCdCUtIuHbvm87n5zzSZbU\nLzlk35+/8XOaWpq45eJbdFpCJAf1eIGjmfUDrgXuiYXOIFp8zGvv45xbBVQD5/QiRxHJMq1trVz9\n56sTFgplY8uYUBx/tvK2t2/jd+//LhPpiUia9WaB45VAMXB/rD0aaHbO7fL1qyd6ykJE8sS/PvOv\nzFs/zxMbPXA0/33Gf3PyiJNpda38af2fmL1iNm20Hexz4ws3UrynmOnHTqekpCTTaYtID/WmWLge\n+Ktzru4Q/YzouoakZs2aRXFxsSdWXl5OeXl5zzMUkbR7aclLzK6c7Z2X3AX1d9Tz1Z1f9XY+Bbiq\no3mg7QDX/ulaBswZxOqVK1UwiKSgoqKCiooKT6yhoSEj792jYsHMSoBPAp/rFK4D+pvZcN/swiii\nswtJzZ49m6lTp/YkHRHJEOcc33/9+95Coa0I/nw37Dwt/glLgTGz4dwHO2JjYf8p+4hEIioWRFKQ\n6At0ZWUlZWVlgb93T2cWridaADzfKbYQaAEuAp4EMLPJQAnwdi9yFJEs8eiKR1m0fZE3+MZNsOn6\nrp807x44bhGMWtERmwG7D+wOJEcRSb+UFzhadCnzV4H7nHMHT0bGZhPuAW43sxlmVgbcC7ylKyFE\ncl9rWys/ee0n3uCOifDGvx3iif3hWd/CxkHwwLoH0pugiASmJ1dDfBIYT7QQ8JsFPAs8BrwK1BDd\nc0FEctyjKx6lKlLlDb54K7QMOvSTq6fDii94QhUfVhBpjKQxQxEJSsrFgnPuJedckXNubYJj+51z\n33XOjXTODXPOfUG7N4rkPuccP3/j595g3WRYeWX3X+SVn0Bbx6+cxtZGXUopkiN0IykROaR56+fx\nwUcfeIOv/gO4FH6FREph2Zc8oTvfu5P9LfvTkKGIBEnFgogc0m8W/MYb+AhYOSP1F3r7f3uadXvq\n+NOKP/U4LxHJDBULIpLUhp0beGbVM97gAohuoZKiutNhg/cyr7hCRESyjooFEUnqd+//DtdpX7Uh\nfYfAoW8H0bW3vaci3q95n2X1y3rxgiISNBULItKl1rZWHljqvcTxivFXQHMvXnTNdPBtCn/v4kQX\nV4lItlCxICJdemXDK9TsrvHErixJ4QqIRNr6xs1MPLD0AZpbe1OBiEiQVCyISJfmLJnjaZ8+5nSO\nG35c71/YtwlkpDHCc6uf6/3rikggVCyISEJ7mvfweNXjnth1p16XnhffDqcffron9MiKR9Lz2iKS\ndioWRCShp1c+TeOBxoPtIiui/KT03QX2sqMu87SfXf0se5v3pu31RSR9VCyISEL+WYVLj7uU0UNH\np+31Lxx7IX2s046OBxp5fs3zSZ4hImFRsSAicfY27+WFtS94YlefcHVa32PEgBFcOPFCT+zRDx5N\n63uISHqoWBCROHPXzWVfy76D7SIr4oopV6T9ffwFyHOrn2NP8560v4+I9I6KBRGJ80TVE572jGNm\ncPigw9P+PleWXkmRFR1s72vZx4vrXkz7+4hI76hYEBGP5tZmnl39rCd2VelVgbzXyMEjmXHMDE/s\nmdXPJO4sIqFRsSAiHq+sf4WG/Q2e2Oc+9rnA3u+Kyd7TG8+tfo421xbY+4lI6lQsiIiHf1bh7KPP\nZtywcYG9n38txEeNH7Fgy4LA3k9EUqdiQUQ8XljnvQrC/80/3SaNmMQJR57gicXd5VJEQqViQUQO\nWrt9LWu3r/XELjvusi56p8/fHf93nvaza57toqeIhEHFgogc9Nc1f/W0xwwdw2ljTgv8ff2nIpbW\nL6W6oTrw9xWR7lGxICIH/XWtt1i49NhLMbPA3/eco8+JuzTzpXUvBf6+ItI9KhZEBIB9B/bx6oZX\nPbFMnIIAKOpTxCcnfdITe+lDFQsi2aJv2AmISOZVV1cTiUQ8sflb53t2bexDH0btGUVlZaWnX1VV\nVSA5XTzpYh5d0bHd87z182hzbZ77R4hIOFQsiBSY6upqpkwppamp0XvgUuCcjmbbpjYuPMd774Yg\nXTzpYk870hhhcd1ipo6dmrEcRCQxlewiBSYSicQKhQeBhR2Picd7O679pvf4wcfPAslrwmETOP5w\nbw5atyCSHVQsiBSsUmBq9DF4PIxZ4z384Zc7jnseEwPLyD+7oHULItlBpyFEBI55zdtuHgw1Hw/0\nLROtfZjoK0Te2PgG89+bz8CigZ74yJEjKSkpCTQ/EemgYkFEYOLL3nb1edDaP6A3qwX6MHPmzPhD\nA4HvcXDOs7mtmWnl02Cdr9vAwaxaVaWCQSRDUj4NYWbjzOwBM4uYWaOZLTGzqb4+PzWzmtjxl8zs\nuPSlLCJpd8wr3vb6CwJ8s51AG3FrJlgITQthy8ne7hOu9/V7kKamxrirOUQkOCnNLJjZYcBbwDyi\na6cjwPHAjk59vg98B/gKsB74OTDXzEqdc81pyltE0mVoLRy50hvbEGSx0K59zYTPhsth/LKO9jGr\nEvcTkYxJdWbhB0C1c+7rzrmFzrmNzrm/OefWd+pzI/Az59wzzrnlwHXAOCC4e9yKSM9N9M0qNA2H\n2hD/OG/8hLd91ALo15i4r4hkRKrFwhXA+2b2qJnVm1mlmX29/aCZTQTGEJ15AMA5twt4F88V3CKS\nNfynIDaeD20hLmeqngZtRR3togNw9Dvh5SMiKRcLk4D/BawCLgF+B/zazNpXKo0BHFDve1597JiI\nZBv/4saMnIJIonlY/MzGhNcS9xWRjEi1WOgDLHTO/cg5t8Q5dzfwP0QLiGSMaBEhItlk+CY4/ENv\nbH3mdm3s0gbfqQj/pZ0iklGpzjXWAv6Lo6uAq2L/Xke0MBiNd3ZhFLAo2QvPmjWL4uJiT6y8vJzy\n8vIUUxSRbit5y9vedxjUnxJOLp1tPB+m3dbRPvod6NsELQO7fo5InquoqKCiosITa2hoyMh7p1os\nvAVM8cWmABsBnHPrzawOuAhYCmBmw4GzgDuTvfDs2bOZOlUrnkUyquRNb3vTNHBZsLFr9XngDCw2\nIdl3P4x7LxoXKVCJvkBXVlZSVlYW+Hun+lthNnC2md1kZsea2ZeArwO/7dTnDuCHZnaFmZ0MzAE2\nA0+nJWMRSZ/xvpmF6mnh5OHXdBjUneqN6VSESGhSKhacc+8DVwLlwDLg34EbnXOPdOpzC/Ab4PdE\nr4IYBFymPRZEssyAPTB6qTdWPT2cXBLxX0KpRY4ioUl5vtE597xz7hTn3GDn3InOuT8m6HOzc25c\nrM+lzrm16UlXRNLm6GXQp62j3doPas4ILx8//yLHo98Baw0nF5EClwUnJ0UkFCWLve2aM6BlUDi5\nJLLpXG97wB4YtTycXEQKnIoFkUI1fom3nS3rFdrtHQ3bJ3lj4+eHk4tIgVOxIFKI+hA9DdFZNq1X\naOefXVCxIBIKFQsihWgM0L/JG/P/Yc4GKhZEsoKKBZFCVOJrR6ZA45GhpJKUv1g4/EMYsi2cXEQK\nmIoFkUI03tfOxlMQAFtPgv1DvbHxSxP3FZHAqFgQKTDOufiZhWxb3NjOFcHms70xFQsiGadiQaTA\n1O2rg2G+YDauV2gXt25hSeJ+IhIYFQsiBWbZTt9VEPtGwLbJ4STTHZvP8bbHVUFROKmIFCoVCyIF\nZvkO38ZGW84kerPYLOU/DdG3GcaGk4pIoVKxIFJg4oqFzWeFk0h3NR0GW0/0xvwLNEUkUCoWRArI\ngdYDrGxY6Q1uyfJiAeLXLRwdThoihUrFgkgBWbZ1Gfvb9nuDW84MJ5lU+E9FHBVOGiKFSsWCSAF5\nd/O73sD2Y6FxZDjJpMJf0BwG2/ZrcyaRTFGxIFJA3t3iKxayfb1Cu49KoXmIJ7Ri54qQkhEpPCoW\nRApIXLGQC+sVILo5U80ZntCKHSoWRDJFxYJIgdjZtJOVEd/ixlyZWYC4UxGaWRDJHBULIgXivS3v\neQMt/aDutHCS6YkExYJzLqRkRAqLigWRAhF3CqJuMrQOCCeZnvAVC7sO7GLdjnUhJSNSWFQsiBSI\n+PUKJ4WTSE81jIc9oz2hBVsWhJSMSGFRsSBSAJxz8ZdNbs6xYgGLm11QsSCSGSoWRArAhp0b+Kjx\nI29wy8nhJNMbKhZEQqFiQaQAxP1RbQS25+Ceyb5iobK2kgOtB0JKRqRwqFgQKQDv1fiuhNgCWX2n\nya749lrY37qfZVuXddFZRNJFxYJIAXi/5n1vYEs4efTavsNhW4knpFMRIsFTsSCS59pcG5W1ld5g\nbTi5pMWWEzxNFQsiwVOxIJLn1mxbw+7m3d5gTTi5pMWWEz1NFQsiwVOxIJLnFtYu9LSPGHAE7O6i\ncy7w7Q/xwUcfsKd5T0jJiBSGlIoFM/uxmbX5Hh90Oj7AzO40s4iZ7Tazx8xsVPrTFpHu8q9XOOGw\nE7romSPqJkNbR9PhWFy3OLx8RApAT2YWlgOjgTGxx/ROx+4ALgc+D5wPjAMe72WOItIL/pmF0uLS\nkDJJk5aBsNUbilvAKSJp1bcHz2lxzn3kD5rZcOB64IvOuddisa8BVWZ2pnNOJxZFMizR4sbSw3K8\nWIDoAs0xHU1/QSQi6dWTmYXjzWyLma0zswfNbHwsXka0+JjX3tE5twqoBs7pfaoikqrV21bHnc/P\n+ZkFiFugubBGxYJIkFItFt4BvgpcCvwjMBF43cyGEK3zm51zu3zPqcfzHUBEMsU/PT9u2DiOHHhk\nSNmkke/Sz5WRlVrkKBKglIoF59xc59zjzrnlzrmXgE8DI4CrkzzNAN10XiQE/m/cZWPLQsokzeqg\nyIoONrXIUSRYPVmzcJBzrsHMVgPHAX8D+pvZcN/swiiiswtJzZo1i+LiYk+svLyc8vLy3qQoUtDe\nr/XOLJwx7owueuaYFpg0bBJrdq05GHq/5n2ml0xP8iSR3FZRUUFFRYUn1tDQkJH37lWxYGZDgWOB\n+4GFQAtwEfBk7PhkoAR4+1CvNXv2bKZOndqbdESkk9a2VhbVLvLEysaWQZ7M1pcWl3qKBS1ylHyX\n6At0ZWUlZWXBzximus/CrWZ2vplNMLNziRYFLcAjsdmEe4DbzWyGmZUB9wJv6UoIkcxbvW01ew/s\n9cTKxuXJaQjiF2pqkaNIcFKdWTgaeBg4AvgIeBM42zm3LXZ8FtAKPAYMAF4Avp2eVEUkFf7FjUcN\nO4oxQ8dQk9N7PXfwXwLavshxaP+hIWUkkr9SKhacc0kXEDjn9gPfjT1EJET+afm8Wa8Qc/zw4ymy\nIlpdKxBd5LiodhHnTTgv5MxE8o/uDSGSp/wzC3lzJUTMwKKBnDjKe1MprVsQCYaKBZE81NrWyqI6\n7+LGfJtZgPgCSMWCSDBULIjkoZWRlTQeaPTE8mlxYzt/AaRFjiLBULEgkof837DHDx/PqCH5dwNY\n/8yCdnIUCYaKBZE8FLdeIQ9nFQBOGX1K3E6O/r0lRKT3VCyI5KG4KyHG5t96BYBB/QZpkaNIBqhY\nEMkzLW0t8Ts35unMAsQXQioWRNJPxYJInlkZWcm+ln2eWL5dNtmZvxDSIkeR9FOxIJJn/OsVSopL\nOHJIHtyWuguJFjnu3r87pGxE8pOKBZE84/9mnY/7K3SWaJGjblctkl4qFkTyjP+21Pl8CgK0yFEk\nE1QsiOSRlrYWltQt8cTyvVgALXIUCZqKBZE8knBxYx5fCdFOixxFgqViQSSP+P9IlhSXMHLwyJCy\nyRzt5CgSLBULInmksrbS0546dmpImWSWFjmKBEvFgkge8Z+rL4T1CtDFIkedihBJGxULInmita01\n7tt0ocwsgG5XLRKkvmEnICKpq66uJhKJeGLrd69n74G9nli/SD8qd3tPTVRVVQWeXxjKxpZx7+J7\nD7ZVLIikj4oFkRxTXV3NlCmlNDU1eg+cDHy+U3sXXHLuJZlMLVT+KyJWRlayt3kvQ/oPCSkjkfyh\n0xAiOSYSicQKhQeBhR2Pcdd6O9ae7z1+8PGzTKabMaeOPtWzyLHNtWmRo0iaaGZBJGeVAp3WJIzd\n7D1ce4H3+EH5cRoi0emUiUMnsnb32oPtpxY8xaDIoITPHzlyJCUlJYHlJ5JPVCyI5ANrg7HetQnU\n5OuVELVAH2bOnBl/6LPA6R3N2x66jdueui3hqwwcOJhVq6pUMIh0g4oFkXwwYh0M8N1psTZfr4TY\nCbQRPQ1T6j1U+wicfmtHe9yxwKMJXqOKpqaZRCIRFQsi3aBiQSQfjPOt/N8zGnaPCyeXjPGdhgGo\nbQI6FQsj10O/KXBAixxFekMLHEXygf8URO1UwEJJJVR1p0Fbp19rfdpgzJKu+4tIt6hYEMkHY30z\nC3m7XuEQDgyGiO/UhP//GxFJmYoFkZznuphZKFD+Qsl/ikZEUqZiQSTXjVgPg3Z6Y7UFOrMA8WPX\nzIJIr/WqWDCzm8yszcxu7xQbYGZ3mlnEzHab2WNmNqr3qYpIQv5ZhcYjoGF8OLlkA//MwpEfQL/G\nxH1FpFt6XCyY2ceBbwD+1UN3AJcT3Xj2fGAc8HhP30dEDiHheoUCXNzYLtEix9Fa5CjSGz0qFsxs\nKNGLnL9O9KLn9vhw4HpglnPuNefcIuBrwDQzOzMN+YqIn9YreB0YApGPeWNatyDSKz2dWbgTeMY5\n97IvfgbRvRvmtQecc6uAauCcHr6XiHTJxf8hLOT1Cu20bkEkrVIuFszsi8BpwE0JDo8Gmp1zu3zx\nemBM6umJSFLFm2DwNm+sUC+b7ExXRIikVUo7OJrZ0UTXJFzsnDuQylMBl6zDrFmzKC4u9sTKy8sp\nLy9PJUWRwuL/xrxvBOw8JpRUsop/ZuHID6DvPmhJfFMpkVxQUVFBRUWFJ9bQ0JCR9051u+cy4Ehg\noZm1r6AqAs43s+8AnwIGmNlw3+zCKKKzC12aPXs2U6cW+LlWkVRp58bE6k4DZ2Cx7yh9WqM7OW4+\nO9y8RHoh0RfoyspKysqCn01M9TTE34CTiZ6GODX2eJ/oYsf2fz8AXNT+BDObDJQAb6chXxHpLG69\nggpuAJqHxi9y1LoFkR5LaWbBObcX+KBzzMz2Atucc1Wx9j3A7Wa2A9gN/Bp4yzm3ID0pi0iU0zbP\nydSUwZFVHW2tWxDpsXTs4OhfizALeBZ4DHgVqCG654KIpNOwj2DoVm9MMwsddEWESNr0+hbVzrkL\nfe39wHdjDxEJytiV3nbTcNhxbDi5ZCP/LMuoFVrkKNJDujeESK4aV+Vt150OTv9JH1R3enSRY7s+\nrTB6aXj5iOQw/WYRyVX+mQWtV/BqHgqRKd6Y1i2I9IiKBZFcNdY3s6D1CvG0bkEkLVQsiOSiocDw\nj7wxbfMcTzs5iqSFigWRXDTW124eAtuODyWVrBa3k+MK6NsUTi4iOUzFgkgu8hcLtaeDKwollaxW\n61vkWNSiRY4iPaBiQSQXxRULWq+QUPMw2DbZG9O6BZGUqVgQyUXjfG2tV+ia1i2I9JqKBZEcs23/\nNij2BTWz0DVdESHSayoWRHLMBzs/8Aaah8BHpeEkkwvidnJcDn33h5OLSI5SsSCSY6p2JthfQYsb\nu1Z3urdd1AKj1oaTi0iOUrEgkmNW7FzhDdScEU4iuWL/cIj4Fjn6t8oWkaRULIjkEOccVQ2+P3Rb\nPh5OMrkkbt2CigWRVKhYEMkhm3dtji5w7EwzC4cWd0WEigWRVKhYEMkh79e87w00FcP248JJJpf4\nZxZGrQMt8xDpNhULIjkkrlioOQOwhH2lk9pEixzDSUUkF6lYEMkh79W85w1ovUL37C+Ov3eGf2Mr\nEemSigWRHOGc62JmQbrFv27Bv2W2iHRJxYJIjvhwx4fsaNrhDdZoZqHb/OsWNLMg0m0qFkRyRNys\nwt4R0DA+nGRyUdxOjnCg7UA4uYjkGBULIjkifr3CCWhxYwr898/oC2t3aSdHke5QsSCSI+LXK5wQ\nTiK5an8xbPNeZhq3wZWIJKRiQSQHtLa1srDWd7dEFQup861biLspl4gkpGJBJAes3raaPc17vEEV\nC6nzXWoad58NEUlIxYJIDohbr9AA7BkZSi45bcuZnubaXWvZ27w3pGREcoeKBZEcEL9eIZw8cl7t\nVGjr2Oe5jTYqaytDTEgkN6hYEMkBcTMLKhZ65sAQ2HqSJ/TulndDSkYkd6hYEMlyB1oPsLhusTeo\nYqHnNp/laapYEDm0lIoFM/tHM1tiZg2xx3wz+1Sn4wPM7E4zi5jZbjN7zMx0uxaRXli2dRlNLU3e\noIqFntviLRYWbFkQUiIiuSPVmYVNwPeBstjjZeBpMyuNHb8DuBz4PHA+0Q1VH09PqiKF6Z3N73ja\nE4ZMgH0hJZMPfIscqxuqqdtTF1IyIrkhpWLBOfecc+4F59za2OOHwB7gbDMbDlwPzHLOveacWwR8\nDZhmZmcme10R6Zp/mvykESd10VO65aNS2D/YE3p3s05FiCTT4zULZtbHzL4IDAbeJjrT0BeY197H\nObcKqAbO6WWeIgXLP7OgYqGXXFHcHhU6FSGSXMrFgpmdZGa7gf3AXcCVzrmVwBig2Tm3y/eU+tgx\nEUnR9n3bWb1ttSd28mEnh5RNHtlyoqepRY4iyfXtwXNWAqcChxFdmzDHzM5P0t8Ad6gXnTVrFsXF\nxZ5YeXk55eXlPUhRJD/4v/EO7DuQ44Yf10Vv6bbN3tmZ92reo8210cd0gZhkr4qKCioqKjyxhoaG\njLx3ysWCc64F+DDWrIytR7gReBTob2bDfbMLo4jOLiQ1e/Zspk6deqhuIgXFfy69bGwZ/fr0Cymb\nPLLFOzuza/8uVkVWUXpkaRdPEAlfoi/QlZWVlJWVdfGM9ElHGd0HGAAsBFqAi9oPmNlkoITomgYR\nSdE7W7zrFc466qwuekpKdh8JvhOmOhUh0rVU91n4DzObbmYTYmsXfgF8AngwNptwD3C7mc0wszLg\nXuAt55xWD4mkyDkXN7Nw9tFnh5RNHtrsbeqKCJGupXoaYjQwBxhL9FY2S4FLnHMvx47PAlqBx4jO\nNrwAfDs9qYoUljXb17CjaYcndtbRZxFZFwkpozyzBeh0UcSCGn2nEelKSsWCc+7rhzi+H/hu7CEi\nveC/ZHLs0LGMHz6eCCoW0mKLt7m0fin7DuxjUL9B4eQjksW09FckS/mnxc86+izMLKRs8lAN9On0\nK7ClrUV3oBTpgooFkSzlX9x49lFar5BWzTBp2CRPyD+bIyJRKhZEslDjgUaW1i/1xM46WldCpNvJ\nI7yXUM7fPD+kTESym4oFkSxUWVtJS1vLwXYf68MZ484IMaP8dMrhp3ja8zfNx7lD7iEnUnBULIhk\nobj7QYw6iaH9h4aUTf46bcRpnnbdnjrW71wfUjYi2UvFgkgWemvTW5621isEY/yQ8Rw5+EhPbP4m\nnYoQ8VOxIJJlnHO8Ve0tFqaVTAspm/xmZpw7/lxPzP//vYioWBDJOmu3r+Wjxo88sekl00PKJv/5\niwUtchSfxAL6AAAcpklEQVSJp2JBJMv4T0GMGTqGiYdNDCmb/DdtvHfWZln9MhqaMnMnP5FcoWJB\nJMu8Wf2mpz1t/DRtxhSgsnFl9C/qf7DtcLqplIiPigWRLOOfWdApiGAN7DuQsrHeW/xq3YKIl4oF\nkSwSaYywMrLSE/NPk0v6ad2CSHIqFkSyiP+yvcH9BnPamNO66C3p4i/I3tn8Dq1trSFlI5J9VCyI\nZBH/9PeZR51Jv6J+IWVTOPwzC3ua97Bs67KQshHJPioWRLJI3HqF8VqvkAmjh47m2BHHemLanEmk\ng4oFkSzR1NLEezXveWLajClz4jZn2qRFjiLt+oadgEihqa6uJhKJxMUXb19Mc2vzwbZhDNo2iMpd\nlZ5+VVVVgedYiKaNn8YDSx842H594+s453TZqggqFkQyqrq6milTSmlqaow/OA24uKPp6hwzzp6R\nqdQK3nkTzvO0N+/azPqd65k0YlJIGYlkDxULIhkUiURihcKDQKn34IQbgU4bMm36e+CmBK/yPPCj\noFIsKJ1naZxzjOg/gh3NOw7G5rw2h8+UfCbhc0eOHElJSUngOYpkAxULIqEoBaZ2NPu0QMlSb5cN\nf+/tc5BOQ/ReLdCHmTNnesNfAE7saP7k/p/wk6d+kvAVBg4czKpVVSoYpCCoWBDJBmMWw8Bd3tjG\nT4STS0HYCbQRN8Oz8U9w4i0d7QnjgGcSPL+KpqaZRCIRFQtSEFQsiGSDCa9525EpsGdMOLkUFN8M\nz4Z+QKdiYUQNFB8BDRMynZhIVtGlkyLZ4JhXve0NM8LIQj46ERoP98aOeS1xX5EComJBJGzWChNe\n98ZULITD9YGN53tj/lkfkQKkYkEkbInWK2zQeoXQ+NeKaGZBRMWCSOj8f4wik2HP2HBykfhZncPX\nwbAtoaQiki1ULIiETesVskv9ybDvMG9MswtS4FQsiIRJ6xWyjyuCau9ujlq3IIUupWLBzG4yswVm\ntsvM6s3sSTOb7OszwMzuNLOIme02s8fMbFR60xbJE2OWwMAGb0z7K4TPv2Zk4svh5CGSJVKdWTgP\n+A1wFvBJoB/wopkN6tTnDuBy4PPA+cA44PHepyqSh/ynICKTYfe4UFKRTjZc4G0fsRaKN4aTi0gW\nSKlYcM592jn3gHOuyjm3DPgqUAKUAZjZcOB6YJZz7jXn3CLga8A0MzszvamL5IFJL3nbOgWRHepO\ng8YjvLFJfwsnF5Es0Ns1C4cBDtgea5cR3RVyXnsH59wqoBo4p5fvJZJfivbHL5z78JPh5CJerg98\neJE3duxLifuKFIAeFwsWvcn7HcCbzrkPYuExQLNzznfROPWxYyLSbvx86Levo+0M1l8YXj7i5S/c\nJs4DawsnF5GQ9ebeEHcBJwDTu9HXiM5AdGnWrFkUFxd7YuXl5ZSXl/c4QZGs5v+mWlMG+45I3Fcy\n78OLve0hERi9NHqKQiQEFRUVVFRUeGINDQ1d9E6vHhULZvZb4NPAec65mk6H6oD+ZjbcN7swiujs\nQpdmz57N1KmJbscrkqeOfdHb9v9xknDtPAa2HxvdlKndpJdULEhoEn2BrqyspKysLPD3Tvk0RKxQ\n+CxwgXOu2nd4IdACXNSp/2SiiyDf7kWeIvll0E4YW+mNrVOxkHX8n4kWOUqBSnWfhbuAa4EvAXvN\nbHTsMRAgNptwD3C7mc0wszLgXuAt59yCNOcukrsmLQDrdGaueTBsOje8fCQx/7qFCa9D36ZwchEJ\nUaozC/8IDAdeBWo6Pa7u1GcW8CzwWKd+n+9lniL5ZdK73vbG86F1QDi5SNfWXxhdeNquX1N0YapI\ngUl1n4U+zrmiBI85nfrsd8591zk30jk3zDn3Befc1vSnLpLDjvUVC1qvkJ2aRkDNGd6Yf28MkQKg\ne0OIZNoRwGG13pjWK2Qv/2dz/F/DyUMkRCoWRDLteF979xjYelIoqUg3rP2Utz1mCQxPenGXSN5R\nsSCSaZN97TWfJroViWSlzefAvhHe2PFvhZOLSEhULIhk0N6WvTDBF1xzeSi5SDe19YW1l3pjx78Z\nTi4iIVGxIJJB7370LhR1CrT203qFXLDm0972JN/nKJLnVCyIZNCb9b5vpBs+Ac3DwklGum/tp7yX\nUPZvgmNCy0Yk41QsiGRIm2vjza2+YkGnIHJD45Gw+SxvzL9QVSSPqVgQyZBFtYvYtn+bN7haxULO\n8Bd2k8G5pPfHE8kbKhZEMuS5Nc95A9uOh+36epoz/OsWDoeNezeGk4tIhqlYEMmQv6z6izegWYXc\nUns67B7rCb1e93pIyYhklooFkQzYuHMjC2sXeoOrrwgnGekhg1Xez+zlupdDykUks1QsiGTAkyuf\n9AYai6M3j5LcUnWVp7lsxzK27NoSUjIimaNiQSQDnqh6whtY9YnoZj+SWzZcAE3FntBTK58KKRmR\nzFGxIBKw+j31vFntu2Sy6oJwkpHeae0fdyriiZVPdNFZJH+oWBAJ2NOrnsbR6RK7/cCHZ3XZX7Kc\n71TEaxteI9IYCSkZkcxQsSASsLhTEGuAlgGh5CJpsO5SONDx+bW6Vp5Z9UyICYkET8WCSIB2Nu1k\n3vp53mBVOLlImhwYDGumeUKPVz0eUjIimaFiQSRAT1Q9QUtby8F2/z79ozMLktt8a05eXPci2/dt\nDykZkeCpWBAJ0MPLHva0zx11LjSHlIykz+rzoKMG5EDbAR7/QLMLkr9ULIgEpHZ3LS+v927ac9lR\nl4WUjaTV/mGw2ht6ePnDifuK5AEVCyIB+dOKP3mughjWfxjTR08PMSNJq6Xe5msbXmPzrs3h5CIS\nMBULIgHxn4K4qvQqBhYNDCkbSbu1MLTv0INNh+OR5Y+EmJBIcFQsiARg7fa1vFfznif2pZO/FFI2\nEogWuGjsRZ6Qv0AUyRcqFkQC8NDShzztUUNGceHEC0PKRoLyqaM+5WkvqltE1Ue6Nlbyj4oFkTRr\nc23ct+Q+T+zqE66mbx/dCyLflI0sY+xQ722r71t8XzjJiARIxYJImr28/mU27NzgiX3ltK+Ek4wE\nqsiK4k4v3bfkPg60HggpI5FgqFgQSbM/VP7B0z5l9CmUjS0LKRsJ2g2n3+Bpb927lWdXPxtSNiLB\nSLlYMLPzzOwvZrbFzNrM7DMJ+vzUzGrMrNHMXjKz49KTrkh229a4jSdXPumJ3XD6DZhZSBlJ0EqP\nLOXc8ed6Yn9Y9Icueovkpp7MLAwBFgPfhs630osys+8D3wG+CZwJ7AXmmln/XuQpkhMeWvYQza0d\nWzT2L+rPtSdfG2JGkglfP/3rnvYLa1/QnguSV1IuFpxzLzjn/q9z7ikg0delG4GfOeeecc4tB64D\nxgGf612qItnNORd3CuLKj13JEYOPCCkjyZQvnPgFhvUfdrDd5tq4d9G9IWYkkl5pXZ5tZhOBMcDB\n2+w553aZ2bvAOcCj6Xw/kTBUV1cTiUTi4gsjC1m2dZkndv6w86msrDzYrqrSZXX5aGj/oZSfVM7d\nlXcfjN1deTc/mP4D+hX1CzEzkfRI97VcY4iemqj3xetjx0RyWnV1NVOmlNLU1Bh/8GrghE7tbfDt\ny76d4GSd5IvOxd95Q87jbjqKhc27NvPLv/6SS8ZdkvC5I0eOpKSkJPAcRdIhUxd+G/qVKXkgEonE\nCoUHgdKOA8W18LHPAG0dsQX/Aq7c9wrPAz8KPE8JWi3Qh5kzZ3rDNwDjO5o3PXkTN/3xpoSvMHDg\nYFatqlLBIDkh3cVCHdHCYDTe2YVRwKJkT5w1axbFxcWeWHl5OeXl/l+2ItmgFJja0fz4D6BPp0Jh\n/1BY/CNguO95Og2RH3YSLQx9ReO7c2H8v3W0S4BxD0DNCd6nU0VT00wikYiKBem2iooKKioqPLGG\nhoaMvHdaiwXn3HozqwMuInZPNjMbDpwF3JnsubNnz2bq1KnJuohkp36NMPV/vLHFX4X9/kJB8o+v\naPzgZNh1Jwzf0hE760V4cmbcM0VSlegLdGVlJWVlwe/j0pN9FoaY2almdlosNCnWbp98uwP4oZld\nYWYnA3OAzcDT6UlZJMucfg8M3u6NLfhOOLlIuNr6wXvf8sZOegSGbUncXyRH9GSfhTOInlJYSHQd\nwi+BSuAnAM65W4DfAL8H3gUGAZc555oTvppILitqhmm3emNrLoNtU8LJR8K38B+gZUBHu+gAnHtb\nePmIpEFP9ll4zTnXxzlX5Htc36nPzc65cc65wc65S51za9ObtkiWOPkhKN7kjb2ReEGbFIjGkbDo\nem/sjN/DkK3h5COSBro3hEhPWStM/y9vbON0qD4vnHwke7z1PWjttCSs3z44e3Z4+Yj0kooFkZ46\n8c8wcrU39qZmFQTYeQws9S1qPPO3MGh7wu4i2U7FgkhP9DkAF/j2S6g7NbpeQQSihaPrtCP+gD1a\nuyA5S8WCSE+c/hc4wrcU57X/S+LbpUhB2jYZll/jjZ19h66MkJykYkEkVf2AGXd7Y5vPhKorQ0lH\nstirN0NbUUe73z6YcXNY2Yj0mIoFkVSdDQzz3Uhq3i/QrILE2TYFKm/wxk7/Ixz5YTj5iPSQigWR\nFNTtqwP/xQ7rLob1F4aSj+SAV2+G5sEd7T5tcImujJDcomJBJAW/+uBX0L9TwBn87b+67C/CnrHw\nzixv7Pj5nltKiGQ7FQsi3fTK+ld4seZFb3DhP0Ct7mkih/DW92D3WG/sU9DYkuBW5yJZSMWCSDc0\nHmjkm89+0xc8HOb9RzgJSW7ZPxzm3u6NFcPdq+9O3F8ky6hYEOmGH778Q9ZsX+MNvvwfsO+IcBKS\n3LP8GvjwIk/ooXUPMX/T/JASEuk+FQsih/Bm9Zvc8c4d3uCmk2HhN8JJSHKUwXN3QkvHopc22rju\nyevY07wnxLxEDk3FgkgSO5t2ct2T1+FwHcEW4Okfgyvq8nkiCW2bEr06opN1O9bxry/+azj5iHST\nigWRLjjnuP7p61m/c733wMtAZGIoOUkemP+vsOkUT+h3C3/Hn5b/KaSERA5NxYJIF3717q94cuWT\nntipI06Ft0NKSPJDW1948qfQ7A3f8JcbWLF1RTg5iRyCigWRBOauncu/vPgvntgRg47gP8v+k85n\nJER6ZPt4+Ks3tPfAXq569Cp2Nu0MJyeRJFQsiPgsq1/GF/78BVpdqyf+wJUPMGbQmJCykryzCD5b\n8llPaPW21Xz2kc/S1NIUUlIiialYEOlkw84NXP7w5exu3u2J/+j8H3HZ8br9tKTX90/6PmVjyzyx\n1ze+zpef/DKtba1dPEsk81QsiMRUN1Rz4f0XsmnXJk/8mhOv4WbdKVACMKBoAE9c8wRjh3p3d3zs\ng8e44S83qGCQrNE37AREMq26uppIxHvXyC2NW/jW299ic+NmT/zUEady44QbWbxoMQBVVVUZy1Py\nX1VVFaWU8supv+Qb87/B3pa9B4/dv+R+6iJ1/PS0n9K3T/yv6pEjR1JSUpLJdKWAqViQglJdXc2U\nKaU0NXXak38c8CVgqK9zPSy5ZQnnNp6bwQylMNQCfZg5c2ZHaCIwE+i0fcfcLXOZ+8pceIy4qycG\nDhzMqlVVKhgkI3QaQgpKJBKJFQoPAguh9Bb46qD4QmHrJJjzEjQujPY7+PhZhjOW/LQTaOPgzyEL\nYf1CeOQOzw6PAEwGbjgWDnumoy8P0tTUGDdDJhIUzSxIYSo6Di6eA2f/Kv5Y/UnwwEuwN9GVDzoN\nIelUCnS6a+maqVBRCl/8LPTrdEXE6HXwzevgL3+AqqsynqWIZhak8IwDvnFd4kJh/QXwxzdhjy6R\nlJCsuyRarDb6blI2aAdc83m44hswcHfi54oERMWCFIwd+3Zw+4rb4evAmNXxHRZfBw++APuLM56b\niEf1dPifBbD1hPhjZX+A71wFp0a3JBfJBBULkvf2Nu/l1rduZdKvJ/HQhw/F/9S3DIC/3A1P3Qet\n/RO9hEjm7ZgE98yH5VfHHxu6Ha6EmW/M5JlVz6hokMBpzYLkrS27tnDne3fyu/d/x46mHYk71Z8M\nT86ButMym5xId+wvhscegbWXwae/A/33eg6vbFjJZx75DFPHTuW7Z36Xa068hkH9BoWUrOQzzSxk\nUEVFRdgpZESY49zTvIcHlz7Ipx78FCV3lPCLN3+RuFBo6Qcv/wzufr+XhUJhfKYaZ5gMFn8V7loO\nKz+TsEdlbSVfe/prHHX7UXz7uW/z6oZXk27opN9FkqrAigUz+7aZrTezfWb2jpl9PKj3yhWF8oOb\nyXG2tLWwqHYRt751KxfNuYjD//twvvzkl5m7bi5tri3xk5YDd/0ZXv9hGk47FMZnqnFmgZ3HwCNP\nw8PPQGRCwi47mnZw1/t3ccH9FzDu9nF85amvcP/i+6luqPb00+8iSVUgpyHM7Brgl8A/AAuAWcBc\nM5vsnNOFwdIjkcYIq7etZv7q+SyuWcyKnSuo2lnF/rb93Xr+tFHTuKT/Jfz45h8D44NNViQoq/8O\n1o6Ck85i9DWjqT9Qn7Db1r1bmbNkDnOWzAFg9MDRfKz4Y0wpnsK6+nX8bcnfmHbCNJ22kG4Jas3C\nLOD3zrk5AGb2j8DlwPXALQG9p+QQ5xx7mvewa/8uGvY3RP/ZFP3ntn3bqNldQ+3uWmr2RP+5sWEj\n2/dtT/2NmoGlwDvwVuQt3uKtdA9FJPPaPoKlfahfXg/HAx8Hjkv+lPqmeuqb6nmt/jXYCRc/dTE8\nBeOGjWPSiEmMGzaO0UNGRx9DRzNqyCiKBxQzbMAwhg8YzrD+wxg2YBiD+g7CzDIxSskiaS8WzKwf\nUAb8Z3vMOefM7G/AOel+v1xw86s3s3b7WiprK7n2iWsPrlx2xP6Z4XbQ77Gsehlld5fR3Np88LG/\nZb+33dq92YAeaS2CD8+GpZfByhlwoPM3p+eBHwX33iIZEdsBsu1BWFUKq4ARm+CkF+GEeTB2Vbdf\nqWZ3DTW7a7rdv8iKGNh3IP2L+tO/qD8D+g44+O/9i/ozoGgA/Yr60cf6eB6GedtmCfuYGUbqxUii\nAmZh7ULKHy/v+jlpep+k/XvwHqnY/mEPvkT1QBAzCyOJ7m7unxurB6Yk6D8Q8vsGPY/Pe5zlW5fD\nVnh47sNhpxO8HbBjYRdXHwRl15Gw5WTYdApsObFTgbDB17H9l+Lz9Gw3xvaZifbnbwYe6sXze6K3\nr9GT53ceZxjvn87nJ3uN7n6eYY+h/fnrO0I7gDfOjT6G18PRS+GoD2BsFQzc4316Ex3/KaSolVb2\nxv6X9bbCIy8+EnYWweo4sT8wyLexdF+fa2ZjgS3AOc65dzvFbwGmO+fO9fX/Eqn9thURERGva51z\ngX0bDWJmIQK0AqN98VHEzzYAzAWuJfoVsCnBcREREUlsIHAM0b+lgUn7zAKAmb0DvOucuzHWNqAa\n+LVz7ta0v6GIiIgEJqirIW4H7jezhXRcOjkYuC+g9xMREZGABFIsOOceNbORwE+Jno5YDFzqnPso\niPcTERGR4ARyGkJERETyh+4NISIiIkmpWBAREZGkAikWUrmJlJm9YmZtCR7PdNH/97Hj/xRE7qlI\n9zjN7N4Ex5/PzGi6FsTnaWalZva0me00sz1m9q6ZHR38aLoWwOfZZmatCfr8n8yMqMvc0z3OIWb2\nWzPbZGaNZrbCzL6ZmdF0LYBxjjKz+8xsi5ntNbPnzewQmywHL5Vxxvr/s5mtjH1W1WZ2u5kN6M1r\nZkK6x2lm55nZX2KfZ5uZJb6lZwgCGOtNZrbAzHaZWb2ZPWlmk1NKyjmX1gdwDdH9Eq4DPgb8HtgO\njOyi/2FE92Bof5wAHAC+nKDv54BFwCbgn9Kde9jjBO4FngOO7NSvOA/HeSzR/Th+AZwCTAT+rqvX\nzOFxjvI9vgq0ABPybJx3A6uB84AS4BuxPn+XZ+N8G3gVmEr0jgy/I7o/zKAcGueXgH2x55UAnyS6\nid5tPX3NHB7np4guwv8c0b2BPhPW+DIw1ueBLwOlwMnAs6n+7AYx0HeAX3VqG9E9VL/Xzef/M9GN\nzwf54kcR3auhlOgep2EXC2kfJ9Fi4Ymwf1gzMM4K4P6wxxb0OBP0eQp4Kd/GCSwD/t3X733gp/ky\nTqLFQRvwMd9r1gPX58o4gd/4fwaB24DX0/X/Xa6M03esjewpFgIda+z4yNiYp3c3r7SehrCOm0jN\na4+5aGap3ETqeqDCObev0+saMAe4xTkX+k0kghpnzIzYNNFKM7vLzA5PS9I9EMQ4Y5/l5cAaM3sh\nNtZ3zOyz6c2++wL+PNvfYxTwaeAPvcu25wIc53zgM2Y2LvY+FxD94xrojnJdCWicAwAHHLwDWuw1\n9wPT05B2yno4zvlAWfu0tplNIvpz+VwvXjNQQYwzW2VwrIcR/Xnu9l2o0r1mIdlNpMYc6slmdiZw\nIvG/UH8ANDvnfpuOJNMgqHH+lejU04XA94BPAM/H/sCGIYhxjgKGAt8nOjV2MfAk8ISZnZeGnHsi\nqM+zs68Cu4iONSxBjfO7xO6qZWbNRD/XbzvnwrofeBDjXEl0ZvMXZnaYmfU3s+8DRwNj05J16lIe\np3OuAvgx8Gbss1oDvOKc+++evmYGBDHObBX4WGN/T+4A3nTOfdDdxILawdHPoNP9kbt2A7DcObfw\n4BPNyoB/Ak4PKLd06vE4IbqZVafmCjNbBqwDZgCvpCvJNOjNONsL1Kecc7+O/ftSMzsX+EfgjfSl\n2Wu9+jx9vgY86JxrTktm6dXbcf4TcBbRdSfVwPnAXWZW45x7Oa2Z9k6Px+mcazGzq4B7iH4bayH6\nbS/0BcgJdDlOM5sB/BvR/9YWAMcBvzazWufcz3vymiEKYpzZKp1jvYvompxpqSSQ7mIh1ZtIHWRm\ng4gu0Pih79B0ogv+NnX6gl0E3G5m/+ycm9TbpHsgiHHGcc6tN7MI0Q8/jGIhiHFGiP6i9Z9OqiLF\nH940CvTzjM2YTAa+0Ls0ey3t4zSzgcB/AJ91zr0QCy83s9OBfwHCKBYC+Tydc4uAqWY2DOjvnNtm\n0fvgvJeWrFPXk3H+FJjjnLs31l5hZkOJLlL9eQ9fM2jpHOfviY4zWwU6VjP7LdFTFOc552pTSSyt\npyGccweAhcBFnZKzWHv+IZ5+DdCf+NtVzyG6Yv7UTo8a4Bbg0rQknqKAxhnHopcSHgGk9KGmSxDj\njL3me8AUX//JwMZeptwjGfg8bwAWOueW9zLVXglonP1iD/+3nlZC2scl6M/TObc7VigcD5xBdOFq\nxvVwnIOJLmzrrK39ub38/y4QaR6nhXha95CCHGusUPgscIFzrronyaV7JefVRC/j6HzZxzbgyNjx\nOcB/JnjeG8DD3XyPbLgaIq3jBIYQLYDOAiYQ/eF4n+g37n75Ms7Ysc8RvTTo60Qvo/wO0Ayck0/j\njB0fDuwBvhHmz2vAn+crwFKia2yOIbo+oxH4hzwb59/HxjiR6C/d9cCjufR5Ej23vZNoUXQM0TVD\nazqP+VCvmUfjHEL0y+dpRP+4/nOsPT4PP9O7gB1EL28e3ekxsNt5BTTYbxG9hnMf0WuTz+h07GXg\nj77+xxP9JnJhN1//Q0IuFtI9TqL3JH8BqCP6h/RD4P+F+R9okJ8n0T8oq4G9QCUhXpMf8Di/QbRY\nGBb2+IIaJ9Ep0nuI7n+yF/gAuDEPx/ldomsymogWCjcDfXNpnERne37U6b+9DcCvgeHdfc18GSfR\nwq8t9pl3fvwxU2PK4FgTjbMVuK67OelGUiIiIpKU7g0hIiIiSalYEBERkaRULIiIiEhSKhZEREQk\nKRULIiIikpSKBREREUlKxYKIiIgkpWJBREREklKxICIiIkmpWBAREZGkVCyIiIhIUv8fnI2fNWLX\n9iMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1048e7320>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.mlab as mlab\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "\n",
    "x = np.linspace(0.75, 0.82, 10e3)\n",
    "plt.hist(F, normed=True, bins=20)\n",
    "plt.plot(x,mlab.normpdf(x, F_oldschool, dF_oldschool), linewidth=3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Part 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Now repeat the experiment above with the following values:\n",
    "\n",
    "$$M_1=40\\times10^4\\pm2\\times10^4\\rm{kg}$$\n",
    "$$M_2=30\\times10^4\\pm10\\times10^4\\rm{kg}$$\n",
    "$$r=3.2\\pm1.0~\\rm{m}$$\n",
    "\n",
    "and as above, produce a plot."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Resulting Force F = 7.82e-01 +- 5.55e-01\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgsAAAFkCAYAAACuFXjcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xd4VGXax/HvkwQINSChKUZUEGNB6UUQVuzYXRejUdeu\nWFGxvbtrQUVdFbGgYEEQjQ17Q3FFUZqEphBQBAxdIr2EtOf94wwhJ2WYGWbmTCa/z3XNBeee58y5\nD8NM7pzzFGOtRURERKQqCV4nICIiIrFNxYKIiIj4pWJBRERE/FKxICIiIn6pWBARERG/VCyIiIiI\nXyoWRERExC8VCyIiIuKXigURERHxS8WCiIiI+BV0sWCM6WOM+dgYs8oYU2KMOSuAffoZY7KNMfnG\nmF+NMZeFlq6IiIhEWyhXFuoDc4EbgL0uLGGMaQN8CnwDHAOMAF42xpwUwrFFREQkysy+LCRljCkB\nzrHWfuynzWPAadbaDmViWUCKtfb0kA8uIiIiURGNPgs9gEnlYhOBnlE4toiIiOyjpCgcoyWwrlxs\nHdDIGFPHWrur/A7GmKbAKcByID/iGYqIiMSPZKANMNFa+1c4XjAaxUJljO/Pqu6BnAK8EaVcRERE\n4tHFwJvheKFoFAtrgRblYs2BLdbagir2WQ4wfvx40tPTI5hadAwePJjhw4d7nUbY6Hy8lZOTQ2Zm\nJjAeKP/5GAyUP5ccILNafp6q23uzN/F0PvF0LhBf57PnO8L5WRoO0SgWpgGnlYud7ItXJR8gPT2d\nTp06RSqvqElJSYmL89hN5xMr0oHyeadUEvO1roafp+r73lQuns4nns4F4u98fMJ2Gz+UeRbqG2OO\nMcYc6wsd4ts+0Pf8MGPM2DK7vAgcaox5zBjT3hgzCPg78NQ+Zy8iIiIRF8poiC7AHCAbp8/Bk8Bs\n4AHf8y2BA3c3ttYuBwYAJ+LMzzAYuNJaW36EhIiIiMSgoG9DWGu/w0+RYa29vIp9Ogd7LBEREfGe\n1oaIgoyMDK9TCCudTyyLp3OJt/cmvs4nns4F4u98wm2fZnCMFGNMJyA7Ozs7HjuciOyT2bNn07lz\nZ5w7gYF8PmYDndHnSaRm2PMdQWdr7exwvKauLIiIiIhfKhZERETELxULIiIi4peKBREREfFLxYKI\niIj4pWJBRERE/FKxICIiIn6pWBARERG/VCyIiIiIXyoWRERExC8VCyIiIuJX0KtOikj45ebmkpeX\nF1DbnJycCGcjIuKmYkHEY7m5ubRvn05+/g6vUxERqZSKBRGP5eXl+QqF8UB6AHt8Dvw7skmJiJSh\nYkEkZqQT2JLTug0hItGlYkGkhgi2r0NqaippaWkRykZEqhMVCyJxbw2QQGZmZlB7JSfXY/HiHBUM\nIqJiQST+bQJKCLxPBEAO+fmZ5OXlqVgQERULIjVHoH0iRETcNCmTiIiI+KViQURERPxSsSAiIiJ+\nqVgQERERv1QsiIiIiF8qFkRERMQvFQsiIiLil4oFERER8UvFgoiIiPilYkFERET8UrEgIiIifqlY\nEBEREb9ULIiIiIhfKhZERETELxULIiIi4peKBREREfFLxYKIiIj4pWJBRERE/FKxICIiIn6pWBAR\nERG/VCyIiIiIXyoWRERExC8VCyIiIuKXigURERHxS8WCiIiI+KViQURERPxSsSAiIiJ+qVgQERER\nv1QsiIiIiF8qFkRERMSvkIoFY8wNxphlxpidxpjpxpiue2l/qzFmkTFmhzEm1xjzlDGmTmgpi4iI\nSDQFXSwYYwYCTwL3AR2BecBEY0xqFe0vAob52h8OXAEMBB4OMWcRERGJolCuLAwGRllrx1lrFwHX\nATtwioDK9AR+sNa+ba3NtdZOArKAbiFlLCIiIlEVVLFgjKkFdAa+2R2z1lpgEk5RUJmpQOfdtyqM\nMYcApwOfhZKwiIiIRFdSkO1TgURgXbn4OqB9ZTtYa7N8tyh+MMYY3/4vWmsfCzZZERERib5gi4Wq\nGMBW+oQx/YB7cW5XzATaAs8YY9ZYax/y96KDBw8mJSXFFcvIyCAjIyMcOYuIiFRrWVlZZGVluWKb\nN28O+3GCLRbygGKgRbl4cypebdjtQWCctXaMb3uBMaYBMArwWywMHz6cTp06BZmiiIhIzVDZL9Cz\nZ8+mc+fOYT1OUH0WrLWFQDbQf3fMd2uhP07fhMrUA0rKxUp8u5pgji8iIiLRF8ptiKeAscaYbJzb\nCoNxCoLXAIwx44CV1tp7fe0/AQYbY+YCM4B2OFcbPvJ1jhQREZEYFnSxYK19x9dh8UGc2xFzgVOs\ntet9TVoDRWV2GYpzJWEocACwHvgY+Nc+5C0iIiJRElIHR2vtSGBkFc+dUG57d6EwNJRjiYh3cnJy\nAm6bmppKWlpaBLMREa+EazSEiMSVNUACmZmZAe+RnFyPxYtzVDCIxCEVCyJSiU04dw/HA+kBtM8h\nPz+TvLw8FQsicUjFgoj4kQ5o+LJITaclqkVERMQvFQsiIiLil4oFERER8UvFgoiIiPilYkFERET8\nUrEgIiIifqlYEBEREb9ULIiIiIhfKhZERETELxULIiIi4peKBREREfFLxYKIiIj4pWJBRERE/FKx\nICIiIn5piWoRCV6dzZC6GJI3Qq2dULAGtkN+cb7XmYlIBKhYEJG9SyiCQydC+4/h0K+gyfJKmx3/\nxfEcPftoTj30VC448gI6tuyIMSa6uYpI2KlYEJGq1d4O3YZB1xcgZcVemxfbYuauncvctXN59MdH\n6dyqM3f0uoMLjriAxITEKCQsIpGgPgsiUpEpga7AzefAifcGVChUJntNNhkTMug0uhOTl08Oa4oi\nEj0qFkTEremvcMVQGAA02FB1u5JE2LI/bDgEtjcGW3XT+evm87exf+OyDy9jy64tYU9ZRCJLtyFE\nZI9jxsIZ1zudFiuzqiv8MhCWngR5h0Nxbd8Ts6F2Z1769CWWmCW89ctb/LH5jwq7j5s3jil/TOGt\nv79FtwO6Re48RCSsdGVBRJwOjKfeCuf+s/JCYfEZMHomvDQTpt0O6zqUKRR8CqBT0048euKj/H7z\n77x3wXt0aNGhwkst27SM48ccT9bPWZE5FxEJOxULIjVd0k648GzoMaLic38dCK99C1mfwOquAb9k\nYkIi5x9xPtnXZDPy9JE0qtPI9fyu4l1c9P5FPPjdg1jr5/6FiMQEFQsiNVntrZB5Ghz2ecXnfgJe\nzILl/UJ++aSEJK7vej3zrpvHcQceV+H5+ybfx12T7lLBIBLjVCyI1FS1t8KlJ0Gb79zxotrw4TXw\nGVBYNyyHatO4DZP/OZlBXQZVeO6/U//L4ImDVTCIxDAVCyI1UeIuuPBcaD3DHd/ZGMZ+C3OPD/sh\nkxKSeH7A8zx72rMY3BM1jZgxgvsm3xf2Y4pIeKhYEKlpTDGcdwkc8o07vr2ZUyis6BXRw9/Y7UbG\nnjOWBOP++hn6/VBe+OmFiB5bREKjYkGkpjnpTjjyXXdseyqM+Q7WHhuVFC455hLGnzu+QsFww+c3\n8NGij6KSg4gETsWCSE1yzDjo9ZQ7tqsBvPEF5KVHNZWMozN4/vTnXTGLJfODTHLW50Q1FxHxT8WC\nSE2x/wI48xp3rLgWvPUhrO7iSUrXdbmO/xz/H1dsW8E2znn7HDbnb/YkJxGpSMWCSE1QDxg4BJJ2\nueOfvgjL+nuS0m7397ufSzpc4or9+tevZH6QSYkt8SgrESlLxYJI3LNwNpCyzh2ecSPMucKTjMoy\nxjDqjFF0atXJFf/0108ZMb2SiaJEJOpULIjEu66ToH252PK+MPGpSpt7oW6tunww8ANS66W64ndN\nuos5a+Z4lJWI7KZiQSSeNf8ZTnnTHdveDN57C0pqeZNTFdJS0sg6P8s1B0NhSSEXvX8R2wu2e5iZ\niKhYEIlXibvg/IsgqdAd/2AsbGvpTU57ceIhJzKk1xBXbFHeIm7/6naPMhIRULEgEr/6PAItfnHH\npg2GJad5k0+Ahp4wlC77u0dnjMoexbfLvvUoIxFRsSASj1rMc4qFsta0h0nDvMknCLUTa/PmeW9S\nv1Z9V/yqT67S7QgRj6hYEIk3CYVwzuWQWLQnVgx8dB8U1/EsrWC0a9qOR0981BVbunEp//rfvzzK\nSKRmU7EgEm96PQmtyo0g+AFYW35IRGwb1HUQvdN6u2IjZoxg2oppHmUkUnOpWBCJJ42XQ98H3LE/\nD4DvPclmnySYBF4+82XqJO65GmKxXP/Z9RSVFPnZU0TCTcWCSDw5ZTDUyt+zXZIAH13t3Iaohtqn\ntufBvz3ois1bN48XZ73oUUYiNVOS1wmISJi0/RLSP3THZl0Pq9pGLYWcnOAWgEpNTSUtLc1vm9t6\n3sb4+eP5+c+fS2P/+t+/uOCIC2jRoEVIeYpIcFQsiMSDxF1w2k3u2PZm8L+hwOdRSGANkEBmZmZQ\neyUn12Px4hy/BUNSQhLPnf4cfV/rWxrbvGszd39zN2POHhNqwiISBN2GEIkHvZ6Epkvcsa8fg/wm\nUUpgE1ACjAeyA3yMJz9/B3l5eXt99eMPOp7MDu5C5LW5rzF1xdQwnoOIVEVXFkSquwZroM/D7tiK\nHjDvMg+SSQc67bVVKB4/8XE+WvQRWwu2lsZu/fJWpl81nQSj33tEIkmfMJHq7m//gdo79mxbA58/\nDza+Pt6tGrbigX7ukR4/rf6Jdxa841FGIjVHfH2biNQ0zX+Bjq+6Y7OvhDWR+e3eazd2u5HDmh7m\nit3zzT3sKtrlUUYiNYOKBZHq7KQ7IaFkz3ZBffj2warbV3O1Emvx2ImPuWLLNy3nuZnPeZSRSM2g\nYkGkujp4JrT7wh2begdsa+VNPlFydvuz6ZPWxxV7aMpD/LXjL48yEol/IRULxpgbjDHLjDE7jTHT\njTFd99I+xRjzvDFmtW+fRcaYU0NLWUQwwMlPu2NbWzrFQpwzxvDEyU+4YpvyN/HQ9w95lJFI/Au6\nWDDGDASeBO4DOgLzgInGmNQq2tcCJgFpwHlAe+BqYFWIOYvIUUCrxe7Y5AegoIEn6URbtwO6ceFR\nF7piI2eNJHdzrkcZicS3UIZODgZGWWvHARhjrgMGAFcAj1fS/kqgMdDDWrt70ll9oiVu5ebmBjR3\nwG7BznpIQjH0Kxdbnw5zrgjudaq5R054hPdz3qeguACAguICHvr+IUafOdrjzETiT1DFgu8qQWfg\nkd0xa601xkwCelax25nANGCkMeZsYD3wJvCYtbakin1EqqXc3Fzat08nP3/H3huH6pgfoGm52KRh\nUFKzpk05uMnBXNPpGp77aU/nxlfnvMqdx91J2/2iN8W1SE0Q7G2IVCARWFcuvg5oWcU+hwAX+I51\nGjAUuB24N8hji8S8vLw8X6EQzEyGQwM/QGIB9P3AHVvVFRafte/JV0P39rmX5KTk0u1iW8yD38Xv\naBARr4TrVxED2CqeS8ApJq6x1lpgjjHmAOAOwG+PpMGDB5OSkuKKZWRkkJGRse8Zi0RUMDMZBnEb\nouOr0LjcLY5vH8T5CNY8rRq24sauN/LEtD0dHsfPH889ve8hvVm6h5mJREdWVhZZWVmu2ObNm8N+\nnGCLhTycxW7LL/XWnIpXG3ZbAxT4CoXdcoCWxpgka22VC9MPHz6cTp3ic3IZkaAl5cPx5err3F6w\n5BRv8okRdx53Jy9mv8i2gm0AWCz3f3c/b//9bY8zE4m8yn6Bnj17Np07dw7rcYK6DWGtLcS5btp/\nd8wYY3zbVa3o8iNQ/gZie2CNv0JBRMrpPBoalRtE9O1QaupVhd2a1W/Grd1vdcXeWfAOc9fO9Sgj\nkfgTyjwLTwHXGGMuNcYcDrwI1ANeAzDGjDPGPFKm/QtAU2PMCGNMO2PMAOAeQFOuiQSq1g7oPcwd\nW9YZlp3gTT4x5vZet9M4ubErpnkXRMIn6GLBWvsOTgfFB4E5QAfgFGvtel+T1pTp7GitXQmcDHTF\nmZPhaWA44J6zVUSq1vEVaLjWHfv2em9yiUGNkxtzR0/3hFQTciaw4M8FHmUkEl9CmsHRWjvSWtvG\nWlvXWtvTWjurzHMnWGuvKNd+hrW2l7W2nrW2nbX2sXJ9GESkKokFcFy5KUx+B3I7epJOrLqx242k\n1HF3iB72w7AqWotIMLQ2hEis6/A6pKx0x77zJpVYlpKcws3db3bFsn7JYsmGJR5lJBI/VCyIxLKE\nIuj9qDv2R3vNgVqFW7rfQv1a9Uu3S2wJj/7wqJ89RCQQKhZEYtkR70LTcr8Zf3+2N7lUA03rNWVQ\n10Gu2Nh5Y7VmhMg+UrEgEqtMCfR5xB1b3Rl+P9qbfKqJ23re5prVsaikiMd/rGzZGhEJlIoFkVh1\n2KfQ4hd3bMq91PR5FfamZYOWXN3palfs5dkvs2brGo8yEqn+VCyIxCQLfR52h/48Ahad40061cyQ\nXkOolVCrdHtX8S6env60hxmJVG8qFkRi0SHfQOuZ7tgP94DVRzYQB6YcyD+P/acr9mL2i2zZtcWb\nhESqOX3ziMSi8vMqbDgEfrnQm1yqqSG9hmDK3LLZsmsLo7NHe5iRSPWlYkEk1rSYB4d+7Y79eCeU\nhGuR2JqhXdN2nJt+riv29PSnKSgu8CgjkepLxYJIrOn1pHt7W3OYd5k3uVRzQ3oNcW2v2rqKrJ+z\nqmgtIlVRsSASSxqthKPK/TCbeSMUJVfeXvzq0boHfdL6uGL/nfpfNNu8SHBULIjEku7PQGKZldsL\n68IsLRi1L8pfXViwfgFfLPnCo2xEqicVCyKxos426DzKHZtzOexI9SafODHgsAGkp6a7Yv+d+l+P\nshGpnlQsiMSKTh9CcpmhfdbA9MHe5RMnEkwCd/RyL189eflkflr1k0cZiVQ/KhZEYkEC0ONNdyzn\nXNjQ1pN04s3FR19MqwatXDFdXRAJnIoFkVhwJJCyzh2bOqTSphK8Okl1uKX7La7YhJwJLNu4zKOM\nRKoXFQsinrPQq1wo9zhY2cOTbOLVtV2upWHthqXbJbaEZ2c+62FGItWHigURrx28AFqVi029o9Km\nErrGyY25ouMVrtjLs1/WFNAiAVCxIOK1XuWG8f3VDhaf6U0uce7m7je7poDeWrCV1+a+5l1CItWE\nigURLzX9FdrNc8em3QY20Zt84twhTQ7h7MPPdsVGzBhBcUmxRxmJVA+abF7ES92ec2/v2A/mXepN\nLh7JyckJuG1qaippaWn7dLzBPQbz4aIPS7eXblzKp79+WqGIEJE9VCyIeKXOFjh2jDs2+2oorOdN\nPlG3BkggMzMz4D2Sk+uxeHHOPhUMfdL60LFlR+asnVMae3rG0yoWRPzQbQgRrxz7mjNr424lCfBT\nTZraeRNQAowHsgN4jCc/fwd5eXn7dFRjDLf2uNUVm7x8MnPXzt2n1xWJZyoWRLxgSqBbuWF7i/rB\n5oM8Scdb6UCnAB7pVb1A0AYeOZCWDVq6Yk9Pfzpsry8Sb1QsiHih7ZfQdIk7NuNCb3Kpgeok1WFQ\nl0GuWNYvWazdttajjERim4oFES90f8a9vRb4o5MnqdRU13W5jjqJdUq3C4oLeOGnFzzMSCR2qVgQ\nibbURdB2ojs2Aygz/l8ir1n9ZmR2cHeufGHWC+QX5XuUkUjsUrEgEm0Vhks2gJ+9SaWmK79exPod\n68n6OcujbERil4oFkWiqs9kZBVFW9t+gyJNsaryjWxzNiYec6Io9PeNprLUeZSQSm1QsiERTxzFQ\ne/ue7ZJEmNXfu3yEW7u7h1HOXzefKblTPMpGJDapWBCJFlNccbhkzrmwOdWbfASA09qdRtv92rpi\nWo1SxE3Fgki0tPsC9lvqjs242ZtcpFSCSeCGrje4Yh/kfMCKzSs8ykgk9qhYEImW8sMl1xwLub29\nyUVcLj/2curXql+6XWyLeWGWhlGK7KZiQSQami2EQ792x2bcjIZLxoaU5BQuO+YyV+yl2S9pGKWI\nj4oFkWgoP1xyeyr8kuFNLlKpG7vd6NrO25HHW7+85VE2IrFFxYJIpCVvgmPGumPZ10BRsjf5SKXS\nm6VXGEb57MxnNYxSBBULIpHX8VWovWPPdkkizKpJq0tWHzd1u8m1PXvNbKatnOZRNiKxQ8WCSCSZ\n4oq3IBaeD1tae5OP+DWg3QDaNG7jij0z45nKG4vUICoWRCLpsM+gyTJ3TMMlY1ZiQiI3dnX3XZiQ\nM4HVW1d7lJFIbFCxIBJJ5SdhWt0JVvTyJhcJyBUdr6BerXql20UlRbw460UPMxLxnooFkUhpthAO\nneSOabhkzGtStwmZR7tXoxyVPYpdRbs8ykjEeyoWRCKl/FWF7c1gwUBvcpGglB9G+ef2P3l34bse\nZSPiPRULIpGQvAWOGeeOabhktXF0i6Pp16afK6b1IqQmU7EgEgkdP6o4XPInDZesTsoPo5y5aiYz\nVs7wKBsRb6lYEAk3A3R7xx1b+HfYeoAn6Uhozmp/Fmkpaa6Yri5ITaViQSTcDgOalBtqp+GS1U5S\nQhKDugxyxd5Z8A5rt631KCMR76hYEAm37uW2V3eGFT09SUX2zVWdriI5aU8/k8KSQkZnj/YwIxFv\nqFgQCaMlW5bAIeWCGi5ZbTWt15SLjrrIFXtx1osUFBd4lJGIN1QsiITR28vedge2NYdfNFyyOrup\nu7uj45pta5iwcIJH2Yh4I8nrBETixYadG/hs5WfuYPa1UFzHm4TiVE5OTsBtU1NTSUtL23tDP45t\neSy903rzQ+4PpbFnZz5LxtFaYlxqDhULImHyyuxX2FVSZpa/4iSYdZ13CcWdNUACmZmZe225W3Jy\nPRYvztnnguHmbje7ioVpK6eRvTqbzvt33qfXFakuQroNYYy5wRizzBiz0xgz3RjTNcD9LjTGlBhj\n3g/luCKxqqikiOd/et4dXHgBbN3fm4Ti0iagBBgPZAfwGE9+/g7y8vL2+cjnHH4OBzR0D33VMEqp\nSYIuFowxA4EngfuAjsA8YKIxJnUv+x0E/Bf4PoQ8RWLax4s/5o/Nf7iDM26qvLHso3SgUwCP9LAd\nsVZiLa7v4p5UK+uXLP7c/mfYjiESy0K5sjAYGGWtHWetXQRcB+wArqhqB2NMAs6vA/8BllXVTqS6\nembGM+7AqiNhZQ9vkpGIuLrz1dROrF26XVBcoGGUUmMEVSwYY2oBnYFvdsestRaYBPgbSH4f8Ke1\ndkwoSYrEsnlr5/HdH9+5gzMuRMMl40vz+s3JOMrdqXHkTyM1jFJqhGCvLKQCicC6cvF1QMvKdjDG\nHAdcDlwVdHYi1UCFqwpbgQUneZKLRNbN3d0zcWoYpdQU4RoNYQBbIWhMA+B14Gpr7cZgX3Tw4MGk\npKS4YhkZGWRkaMiSxIa8HXm88fMb7uAsoLiWJ/lIZHVq1anCMMoRM0ZoGKV4Jisri6ysLFds8+bN\nYT9OsMVCHlAMtCgXb07Fqw0AhwIHAZ8YY3Zfk00AMMYUAO2ttVX2YRg+fDidOnUKMkWR6Hkp+yV2\nFe8ZLplkkijKLvIwI4m0W7rf4ioWZqyawYyVM+jeuvw83yKRV9kv0LNnz6Zz5/AO6w3qNoS1thBn\nTFL/3TFfEdAfmFrJLjnA0cCxwDG+x8fA/3x/XxFS1iIxoLC4kJGzRrpiJ+9/MmzzKCGJinMOP4cD\nGx3oio2YMcKjbESiI5TREE8B1xhjLjXGHA68CNQDXgMwxowzxjwCYK0tsNYuLPvAGSy91VqbY63V\nr2BSbX2w6ANWblnpil14yIUeZSPRkpSQxA1db3DF3l34Lqu2rPIoI5HIC7pYsNa+A9wOPAjMAToA\np1hr1/uatKaKzo4i8aR8x8aerXtyZOMjPcpGounqzldTN6lu6XZRSREvznrRw4xEIiukGRyttSOt\ntW2stXWttT2ttbPKPHeCtbbKORestZdba88L5bgisSJ7dTY/rvjRFSvfU17i13519+OSDpe4YqOy\nR5FflO9RRiKRpVUnRUJQfqrf/Rvuz/np53uUjXihfHG4fsd6sn7OqqK1SPWmYkEkSOu2rSPrF/cP\nhUFdBlErUcMla5Ijmx/JiYec6IqNmDECZ546kfiiYkEkSKOzR7tm7auTWIdrOl/jYUbilZu7ua8u\nzFs3jym5UzzKRiRyVCyIBKGguIAXZr3gimUcnUGz+s08yki8NOCwARza5FBXTMMoJR6pWBAJwoSF\nE1izbY0rVv63S6k5EkwCN3Vzry764aIPWb5puTcJiUSIigWRIJT/rbFPWh86turoUTYSCy7veDkN\najco3S6xJTw/83kPMxIJPxULIgGasXIGM1bNcMU0XFIa1WnE5cde7oq9POdlthds9ygjkfBTsSAS\noKemP+Xabt2oNeccfo5H2UgsuanbTZgyS5Jvyt/E6/Nf9zAjkfBSsSASgOWblvPewvdcsRu73khS\nQrgWbpXqrF3Tdpze7nRX7JkZz2gYpcQNFQsiAXh2xrOU2JLS7fq16mu4pLjc0v0W13ZOXg5fL/3a\no2xEwkvFgshebNm1hZdmv+SKXdnxSprUbeJRRhKLTjzkRNJT012xp6Y9VUVrkepFxYLIXrwy+xW2\nFmwt3TYYbulxi589pCYyxlS4ujDx94n8vO5njzISCR/dcBXxo6ikqMJwyXPTz+WQJod4lJEEKycn\nJ6j2qamppKWlhXSsS4+5lH99+y/yduSVxp6a/hRjzh4T0uuJxAoVCyJ+vJ/zPn9s/sMVu63HbR5l\nI8FZAySQmZkZ1F7JyfVYvDgnpIKhbq263ND1Bh747oHS2Bvz3+DhEx5m/4b7B/16IrFCtyFEqmCt\n5clpT7pi3Q/oTq8De3mUkQRnE1ACjAeyA3yMJz9/B3l5eZW+YiAGdR1EclJy6XZhSSHPzXwu5NcT\niQUqFkSqMG3lNGaumumK3dbzNowxVewhsSkd6BTgI72K1whc8/rNubTDpa7YC7NeYFvBtn1+bRGv\nqFgQqUL5nuxpKWmcl36eR9lIdXJbT/etqk35mxgzR/0WpPpSsSBSiaUbl/LBog9csVu636JJmCQg\n7VPbc+ZhZ7piw6cPp7ik2KOMRPaNigWRSoyYPsI1CVPD2g25qtNVHmYk1c0dve5wbS/btKxCASpS\nXahYECmFB81uAAAgAElEQVRnw84NvDLnFVfs6k5X06hOI48ykuqoT1ofuuzfxRV7YuoTmgJaqiUV\nCyLlPD/zebYX7lkxMNEkanVJCZoxhjt6uq8uzFg1g6krpnqUkUjoVCyIlLGjcAfPzHzGFRt41EAO\nanyQRxlJdXb+EedzUIr7/84T057wKBuR0KlYECnj1TmvumbfA7jruLs8ykaqu6SEJG7tcasr9tGi\nj1ict9ijjERCo2JBxKewuJAnprp/6zut7Wl0aNHBo4wkHlzZ8UpS6qSUblssj//4uIcZiQRPxYKI\nzzsL3qkwtfPdve/2KBuJFw3rNOT6Lte7Yq/Pf52VW1Z6lJFI8FQsiOBM7fzYj4+5Yj1a96BPWh+P\nMpJ4cmuPWytMAf3k1Cf97CESW1QsiABfLPmCn/90LyV893F3a2pnCYsWDVpwZccrXbHRs0dX6B8j\nEqtULIhAhasK6anpnNn+zCpaiwTvjl53kGgSS7d3FO7gmRnP+NlDJHaoWJAab9qKaXz/x/eu2JBe\nQ0gw+nhI+LRp3IaLjr7IFXt25rNs3bXVo4xEAqdvQ6nxHv3xUdf2AQ0P4OIOF3uUjcSz8h1mN+Vv\n4sVZL3qUjUjgVCxIjTZ37Vw+XvyxK3Zbz9uonVjbo4wknh3R7AjOOfwcV+yp6U+RX5TvUUYigdES\nelKjPfT9Q67t/erux9WdrvYoG4kVOTk5AbdNTU0lLS0t4Pb39L6HDxd9WLq9dttaxs4dy7Vdrg0q\nR5FoUrEgNdYvf/7ChJwJrthtPW6jYZ2GHmUk3lsDJJCZmRnwHsnJ9Vi8OCfggqHbAd3of3B/vln2\nTWns8amPc2WnK7UEusQs3YaQGuvhKQ+7thsnN+bGbjd6lI3Ehk1ACTAeyA7gMZ78/B3k5QU3BPKe\n3ve4tpduXMob89/Y1+RFIkZlrNRIi/IW8fYvb7tit3S/hZTklCr2kJolHegUsVc/4eAT6HZAN2au\nmlkaG/r9UC7ucLGuLkhM0pUFqZEenvIwFlu63bB2Q27pfouHGUlNYozhvr73uWK/b/xdVxckZqlY\nkBpnyYYlvPnzm67YTd1uokndJh5lJDXRaW1Po9sB3Vyxod8PpaikyKOMRKqmYkFqnEemPEKJLSnd\nrl+rPoN7DvYwI6mJjDHc3/d+V+z3jb8zfv54bxIS8UPFgtQoSzcuZdy8ca7YoK6DSK2X6lFGUpOd\n2vbUClcXHvr+IV1dkJijYkFqlAe+e4BiW1y6XTepLrf3vN3DjKQmq+rqwuvzXvcmIZEqqFiQGiNn\nfU6FS7zXdbmOFg1aeJSRiHN1ofsB3V2xh6Y8RGFxoUcZiVSkYkFqjPsm31ehr0L5ufpFos0Yw/39\n7nfFlm5cyuvzdXVBYoeKBakR5qyZw7sL33XFbu1xK83rN/coI5E9Tjn0lApXFx747gGtGSExQ8WC\n1Aj//vbfru3GyY25o9cdHmUj4maM4YF+D7hiuZtztSKlxAwVCxL3pq2Yxme/feaKDek1hMbJjT3K\nSKSikw89mb4H9XXFHp7yMFt2bfEoI5E9VCxI3Pu///2fa7t5/ebc3P1mj7IRqZwxhmH9h7lieTvy\neGraUx5lJLKHigWJa5OWTuLb5d+6Yvf0vocGtRt4lJFI1Xoe2JOz25/tij057Un+3P6nRxmJOFQs\nSNwqsSXc+fWdrljrRq25rst1HmUksncPn/AwBlO6va1gG49MecTDjES06qTEsTfmv8GctXNcsf8c\n/x+Sk5I9ykjiVU5OTlDtU1NTSUtLq/S5I5sfyaXHXMrYeWNLYy/MeoFbe9xKm8Zt9iVNkZCpWJC4\ntLNwZ4W+Cump6Vze8XKPMpL4tAZIIDMzM6i9kpPrsXhxTpUFw/397ifrlywKigsAKCgu4L7J9zH2\nnLGVtheJNN2GkLj0zIxnWLFlhSv2+EmPk5Sg+ljCaRNQAowHsgN8jCc/fwd5eXlVvmqbxm24vsv1\nrtjr815n9prZETgHkb0LqVgwxtxgjFlmjNlpjJlujOnqp+1VxpjvjTEbfI+v/bUX2Vd5O/J45Af3\nPd5+bfoxoN0AjzKS+JcOdArwkR7QK97b515XR1yLZfDEwVhrw5q5SCCCLhaMMQOBJ4H7gI7APGCi\nMaaqZfv6Am8C/YAewArgK2NMq1ASFtmbod8NrTA2/YmTnsAYU8UeIrGnef3m3Nv7Xlfs+z++54NF\nH3iUkdRkoVyTHQyMstaOAzDGXAcMAK4AHi/f2Fp7SdltY8xVwPlAf5xrdyJhs2TDEkbOGumKXXT0\nRXTev3PIr5mbm+v3knFZwXZ0E/FncM/BjMoexR+b/yiNDfl6CAPaDaBOUh0PM5OaJqhiwRhTC+gM\nlF7jtdZaY8wkoGeAL1MfqAVsCObYIoEYPHEwRSVFpdu1E2vz8AkPh/x6ubm5tG+fTn7+jnCkJxKU\n5KRkHj/pcQa+N7A0tnTjUp6d+aymK5eoCvY2RCqQCKwrF18HtAzwNR4DVgGTgjy2iF+f//Y5n/76\nqSt2c7eb92m4WV5enq9QCLQD29CQjyVSmQuOuIDjDjzOFRv6/VDWb1/vUUZSE4Wra7gB9trrxhhz\nN/APoK+1tmBv7QcPHkxKSoorlpGRQUZGRqh5SpzaVbSLW7+81RVrUb8F/+777yr2CNbuDmx7o9sQ\nEl7GGIafMpxuL3crjW3ZtYX/fPsfXjjjBQ8zk1iQlZVFVlaWK7Z58+awHyfYYiEPKAZalIs3p+LV\nBhdjzB3AnUB/a+2CQA42fPhwOnUK5Ataarqnpz/Nbxt+c8UeO/ExGtVp5FFGIuHT9YCuZHbIZPz8\nPd28RmWP4qpOV+1Tfxyp/ir7BXr27Nl07hze/xdB3Yaw1hbiXGvtvztmnC7m/YGpVe1njBkC/B9w\nirV2TlXtREKxeutqhn7vvvzfo3UPLjnmkir2EKl+hvUfRt2kuqXbFsugzwdRYks8zEpqilBuQzwF\njDXGZAMzcUZH1ANeAzDGjANWWmvv9W3fCTwIZAC5xpjdVyW2WWu371v6InDn13eyvXDPfyWD4ZlT\nnyHBaM4xiV3BjJzZPT30v47/l2tm0pmrZvLy7Je5pvM1kUhRpFTQxYK19h3fnAoP4tyOmItzxWB3\nb5vWQFGZXa7HGf3wXrmXesD3GiIhm7x8Mm/8/IYrdkXHK+h6gOb9klgV/BTRu6eHvr3n7YydN5Zf\n//q19Lm7J93NuYefS7P6zSKQq4gjpA6O1tqRwMgqnjuh3PbBoRxDZG/yi/K59tNrXbGUOik80l8r\n9EksKztFdCCzOeaQn59JXl4eaWlpPH/685z0+kmlz27M38jdk+7mlbNfiVC+IlobQqqxYVOGuX7D\nAnjohIdoXr+5RxmJBCPQKaLdBcWJh5zIwCMHumKvzn2VqSuq7DYmss9ULEi1lLM+h2E/DHPFuh/Q\nvcLiOyLx6MmTn3StGwFw9SdXs6tol0cZSbxTsSDVTokt4dpPr6WwpLA0lmgSGX3maBITEj3MTCQ6\nDmh0AA/0e8AVW7h+IY9M0S04iQwVC1LtvDrnVabkTnHFbu95Ox1adPAoI5Hou6nbTXRs2dEVe+SH\nR5i/br5HGUk8U7Eg1crKLSu54yv3nPgHNz6Y+/rd51FGIt6olViLV856hUSz52paUUkRV358pWt9\nFJFwULEg1Ya1lis/vpLNu9xTmb4w4AXq1arnUVYi3unYqiN3HXeXKzZr9SyGTxvuUUYSr1QsSLXx\n0uyX+Or3r1yxSzpcwiltT/EoIxHv/bvvvzk89XBX7D+T/1NhpJDIvgjXQlIiEbV803Ju/+p2V2z/\nhvsz4tQRHmUkEl3+Zny8q/1dXJF3Bda3nl9+UT6Z72fy4xU/UiuxVrRSlDimYkFiXokt4fKPLmdb\nwTZX/JWzXqFJ3SYeZSUSLQHO+Hgq0GPP5k+rf+Kh7x/igb89UOUuIoFSsSAx7+npTzN5+WRX7KqO\nV3Fq21O9SUgkqgKc8fGbndD2Ykj9ozT00JSHOLXtqfQ8sGekk5Q4p2JBYlr26mzunnS3K5aWksaT\npzzpUUYiXtk942MVCoEJE+Cq7pBYDDhX5TI/yGTutXNpWKdhVLKU+KQOjhKztu7ayoUTLnRNvmQw\njDl7DI3qNPIwM5EYtaYzTL7OFVq6cSm3fHmLRwlJvFCxIDHrxi9uZMmGJa7Y3b3v5oSDT6hiDxHh\nh8sg1x0aM3cMr8973Zt8JC6oWJCYNH7+eMbNG+eKdT+ge4UpbkWkHJsI70P9pPqu8HWfXcfC9Qs9\nSkqqO/VZkJiz4M8FXPep+1JqozqNyDo/S8PARAKxCS5rehkj140sDe0o3MEZ485gXJ9x1EuqOIlZ\namoqaWlp0cxSqhEVCxJTNudv5ty3z2V74XZXfPQZozm4ycEeZSVSnThDLUdePxLOALrseWbZtmX0\neawPvF9xr+TkeixenKOCQSql2xASM0psCZd+eCm/bfjNFb+q41UMPGqgR1mJVDdlhlp+ORXWtHc/\n3QHoeheQXeYxnvz8HeTl5UU5V6kuVCxIzBg2ZRgfL/7YFeuyfxeePf1ZjzISqc7SoagnvPMZ5Jcb\nPXTaE9BmM85QzE74nb9BBBULEiM+WfwJ//72365Yar1UJvxjAslJyR5lJRIHNh4KH41xxxKK4R9/\nhyZLvclJqh0VC+K5OWvmkDEho3Ree4AEk8Dbf3+btBTdPxXZZznnwff3umP1NkDGWVBnizc5SbWi\nDo7iqZVbVnJG1hkVOjQO6z9M8ymIhNO3Q6H5Ajj8oz2x5gvgvIvh7X9Dif/FqsrT6ImaRcWCeGZb\nwTbOzDqT1VtXu+KXHXMZQ3oN8SgrkThlE+D91+HKXtDilz3x9p/C6UXwqdn7YlVlaPREzaJiQTxR\nUFzABe9ewNy1c13xvgf1ZfSZozHGeJSZSBwraAhZH8M1XaHeX3viXb6ErcB3e1msqlQO+fmZ5OXl\nqVioIVQsSNQVlxRz6QeX8uWSL13xdvu14/2B71M7sbZHmYnUAJsOhrc+hEtPhKRde+J/A7blQPbF\nnqUmsUsdHCWqrLUM+mwQby942xXfr+5+fHbRZ+xXdz+PMhOpQXJ7w3tZUFLuR8CAYXDEu97kJDFN\nVxYkaqy13PPNPYyePdoVr1+rPp9f9DntmrbzKDORGmjRufD5c3DGoD2xhBI4/yIoSXKe34tgOkSC\nOkVWZyoWJCqstdz7zb089uNjrnjtxNp8dOFHdG/d3aPMRGqwWddDw9XQ96E9scQiuGAgvD0Bfj2z\nih2dKaWD6RAJ6hRZnalYkIiz1jLk6yE8Oe1JVzzBJPDW+W/R/5D+HmUmInz7INSdDt0m7YklFjqT\nNr39Pvw2oJKdykwpHfDsj+oUWZ2pWJCIstYyeOJgRswY4YobDK+c9Qrnpu/9UqeIRJKBLy6FhEmu\nRadIKoALz3WGWy6oam2WdJzpoiXeqViQiCksLuSaT6/htbmvueIGw6tnv8o/j/2nJ3mJSDk2AT4D\nEs6GTmUmbUoshL9nQPImyL42LIfSxE/Vk4oFiYhtBdv4x7v/4IslX7jiCSaBseeMJbNDcPc6RSTC\nLPDJ/0FCYzh27J64sXDmdc68DFPuAUKdAyX4fg7q4xA7VCxI2P25/U8GvDmAWatnueKJJpHx543n\nwqMu9CgzEfHLJsJHr8KuhtD9Ofdz/f8PGi+Hz553uisELdh+DurjEEtULEhYzVs7j3PePoflm5a7\n4nWT6vLW39/irPZneZOYiATGJsAXz8DOptDvAfdznV+C/ZbAOwNhZ6gHUD+H6kjFgoTNuwve5Z8f\n/ZMdhTtc8aZ1m/JJxif0PLCnR5mJSHAMTL4fdu4Hp93ifurgb+GqHMgC8iKfifo4xAYVC7LPikqK\nuO/b+3jkh0cqPHdQykFMzJxI+9T2HmQmIvtkxs2wdX8491KoVeZSQtO1cDXwyUT4JVJXCdTHIZao\nWJB9smLzCi5+/2Km5E6p8FzP1j15f+D7tGzQ0oPMRCQsFv4dNrWBjLOg4Zo98TrA3++Fg1bAxKeg\nKDnMB1Yfh1iiYkFC9kHOB1z58ZVszN9Y4bmrOl7Fc6c/R52kOh5kJiJhtboLvDTTKRhazXE/1/UF\naD3dmY9h/ZEROHhwfRw0BXVkqFiQoG3K38TtE2/n1bmvVnguKSGJEaeO4Pou12uZaZF4sqU1vPoD\nnH4jdBzjfq7VHLi2E/zvIZh2mzOqIuo0BXUkqViQoHy06COu/+x61mxbU+G5No3bkHV+Fj1a9/Ag\nMxGJuMJ6ztDKP2rBgNFQq8xzSQVw8p1w+IdOm7+i3U8p9Cmop0yZQnp6YPvU1CsRKhYkICu3rOT2\nr27nnQXvVPr8P478B6POGEXj5MZRzkxEom7u8bB6NFxwMDRb5n4ubSoMOhqm3gHf/x8U1o9ycsHc\ntlAnykCpWBC/dhbu5ImpT/Doj49WGBIJ0LB2Q4afMpwrOl6h2w4iNcmfwKg34G/vQ68nnZked0ss\nhD7D4Og3YOJwyInVNWBC60QZzJUIiI+rESoWpFLFJcW89ctb/N///o8/Nv9RaZvT2p7GqDNGcWDK\ngVHOTkRiQlEd+Pq/sPhsOOefsN/v7ucb58LA8yH3OJh0AuR6kmUAAr0aUXP7RahYEJcSW8KEhRO4\n/7v7Wbh+YaVtmiQ3YcSpI8jskKmrCSICub3hhXnQ90Ho+RQkFrmfT/sRrvgRFgPfLoa11XUGx5rb\nL0LFggDOCpHvLnyXx398nHnr5lXaJsEkcE2na3jwbw/SrH6zKGcoIjGtsD5MegzmXQan3wAHT67Y\npj3Q/iL47XX44W74ow+hL0zlpcj2i6hTJ5kJE96jVatWAbWPRnGhYqGG27hzIy/NfolnZjzDqq2r\nqmzX96C+PH3q0xzb8tgoZici1c76I2Ds/+Cot6H/vdBkWcU27b5wHit6wPTBsOgcKK4d/VyjItir\nEVPYtes2zjjjjICPEI3bHCoWaiBrLVNyp/DqnFd5d+G7lXZc3O2YFsfw4N8e5MzDzoyLWw65ubnk\n5QU+oX2wE7yICICBXy6EnPOg82g4fig0+LNiswOnw4EDYVtzmHMlZF8Nmw6OfrpREejViBz2tdNl\nJL63VCzUIL9v+J13FrzDmLlj+G3Db37bHtHsCB7o9wDnpZ9HgkmIUoaRlZubS/v26eTnV10ciUgY\nFdeGmTfC3H9C9yugx7tQ2UjKBn86oyd6Pwp/HA8/Z8BCuw8rW8aDyHa6DJaKhTi3OG8xE3Im8N7C\n95izds5e2/dO681tPW7jrPZnkZjgxSxskZOXl+crFILpnPQ58O/IJSVSExQ0gClnw/R3oeOd0Ott\naFzJKCtjoc13zuP0RPgdyPkQfmsF2wK7f1/zVHabIwcIb/GgYiHObNm1hW+XfcvE3ycy8feJLN24\ndK/7JJpEzj/ifG7veTvdDugWhSy9FkznJN2GEAmbQmDmQJj1EBzxHnR5Edp8X3nbxGI4DDhsKDAU\nVneGXwfAktOcv5fUqny/Giu4NTSCpWKhmluzdQ3TVk5j6oqpTF0xlZ9W/0RRSdHedwTaN23PlR2v\n5JJjLtHKkCISPSW14JcM59FsIXQeBceMg7qbqt5n/2zn0e9BKKgHK3rBH32d2xaruziFiESMioUo\nyMrKIiMjY59ew1rLii0rmL9uPvPWzmP+n/P5adVPLNtUSU9jP5rVa8Z56edx6TGX0rN1z5A6LYbj\nfGJLFhAv5xNP5wLxeT7isv4I+HKEM+yy7RdwdBYc9gnUyq96n9o74NBJzgOgJAHy9ofVwJosWL0T\n/jwKdqUEkUi8/V8Lr5CKBWPMDcAdQEtgHnCTtfYnP+0vAB4E2gC/Andba78I5djVUTA/XDfnb+b3\njb+zZMMSlmxYwu8bfue3Db/x858/synfT9XtR8sGLTnv8PO44MgL6JPWZ5/7IqhYiGXxdC4Qn+cz\n0OskYlNRMiw613nU2QLt74L2L8Kh9SF5u/99E0qg+UpoDhz7BPCEE9/aCvIOh/Xpzp9/HQab2sDm\nNCiqW+5F4u3/WngFXSwYYwYCTwLXADOBwcBEY8xh1toKY9KMMT2BN4G7gM+Ai4APjTEdrbWVTxEY\nZyyWTfmb2LBzA3/t+Iu129ayausqVm9dzeqtq0v/vnLLSjbs3LDPx6udWJveab05+ZCTOaXtKXRo\n0SFuRjSISA2wqxHM7w3zX4TEb+DAHXDYZ9D2S2i+IPDXabjGeRz8bcXntjeDTQ1gM7D5SZi5BA4a\nA9ubw/YWvj+bO0WMhHRlYTAwylo7DsAYcx0wALgCeLyS9rcAX1hrn/Jt32eMORm4ERgUwvEjosSW\nUFBc4Pexo3AH2wu2s61gG9sLfX/6tsvGNuVv4q+df7Fh5wY27NxA3q95NHmsScRyb1i7IT1a96Bn\n6570OrAXvdN6U792tFd6ExGJgOJasPxvzuOrJ6Deekj7wRkxcdD30GKec2UhWPXXO48DAN6EZcA5\nV1Rst6sB5Dd2HrsKIB/YdS/kH+yLpUBBfedKRWFdZxnvIt+fhXWhcCUUAYXrofgvp79GcS3nz5JE\nqssMlkEVC8aYWkBn4JHdMWutNcZMAnpWsVtPnCsRZU0Ezt7b8W7+4mYa5jSkxJYE/CguKd5rm6KS\nIgpLCl2FQKCdAkNi994kUCl1UujQogMdWnTgmBbH0L11d45sdmTcDXMMVDCTLGmCJZE4sKPZntsV\nALW2Q8t50GoU7D8OWh0KzZZDQnF4jldnm/NIWVkmODGEFzq18nBx2eIhCYqLnJGQxWdBSQMnZhOc\nwsImVPL4y/kZY68Fm+LE8rZBmG/0B3tlIRVIBNaVi6/DmfW7Mi2raO+v+30ywI+zf4zhVcqCkI/T\n8SZABkPz+s1pVrsZzZOb06JuC1rVbUVa/TRSk1OdRhudx8+LfuZnfiYhIYGSksCr62Dbl91n5cqV\nvPHGGxE9RiDy8vIYMuRuCgv9dISq1Oe4h0SuBKo6nx+r2KcqwbYP9zEqOxevc9qXfcqeT6zktC/t\nV0bhGLF43hE8RiGwAlix+2rqPyHhQGi0DhqvhiZrnD8br4KG6yF5W+WvE+T3dPgUUvlQjiqm3zdU\ncTFi1p6/7umOEbZ7KMbawH/tNca0wjmDntbaGWXijwO9rbW9KtlnF3CptfbtMrFBwL+stftXcZyL\nqPrbW0RERPbuYmvtm+F4oWCvLOQBxUCLcvHmVLx6sNvaINuDc43nYmA5Tr0nIiIigUnGGX0Yyv2S\nSgV1ZQHAGDMdmGGtvcW3bXBuFjxjrf1vJe3fAupaa88uE/sRmGetjZkOjiIiIlK5UEZDPAWMNcZk\ns2foZD3gNQBjzDhgpbX2Xl/7EcB3xpjbcIZOZuB0krx631IXERGRaAi6WLDWvmOMScWZZKkFMBc4\nxVq73tekNc5Akd3tpxljMoCHfY/fgLNryhwLIiIi1V3QtyFERESkZtG0fiIiIuKXigURERHxKyaK\nBWPMvcaYH40x240xAS2OYIwZY4wpKff4PNK5BiKU8/Ht96AxZrUxZocx5mtjTNtI5hkoY0wTY8wb\nxpjNxpiNxpiXjTF+55M2xkwu994UG2NGRivnMnncYIxZZozZaYyZbozpupf2Fxhjcnzt5xljTotW\nroEI5nyMMZeV+bff/T7siGa+VTHG9DHGfGyMWeXL66wA9ulnjMk2xuQbY341xlwWjVwDEez5GGP6\nVvL9VWyMaR6tnP3kdo8xZqYxZosxZp0x5gNjzGEB7BeTn51QzifGPzvX+f59N/seU40xVUwPWbrP\nPr83MVEsALWAd4AXgtzvC5xOli19j1hZMizo8zHG3IWzXsa1QDdgO84CXbUjkmFw3gTSgf4464Ac\nD4zayz4WGM2e96cVcGcEc6zA7Fn07D6gI84KqRN9HXQra7970bOXgGOBD3EWPTsiOhn7F+z5+Gxm\nz+ejJXBQpPMMUH2cztE3EMCE6MaYNsCnwDfAMTijrF42xpwUuRSDEtT5+FigHXvem1bW2j8jk15Q\n+gDPAt2BE3G+z74yxpRfprFUjH92gj4fn1j97KzAWZixs+/xP+AjY0x6ZY3D9t5Ya2PmAVwGbAiw\n7Rjgfa9zDuP5rAYGl9luBOwE/uHxORyOM1N5xzKxU3BGvLT0s9+3wFMe5z4dGFFm2+DMt3tnFe3f\nAj4uF5sGjPT6/1KI5xPw/z+Pz6sEOGsvbR4D5peLZQGfe51/iOfTF2eCu0Ze5xvA+aT6zqm3nzYx\n/dkJ4XyqxWenTL5/AZdH8r2JlSsLoernu6y0yBgz0hizn9cJhcIYczBO5frN7pi1dgswg6oX6IqW\nnsBGa+2cMrFJOL8Vdd/LvhcbY9YbY342xjwSQCUfNmbPomdl/00tTu7+Fj2bVC420U/7qAnxfAAa\nGGOWG2NyjTGx8pteKHoQo+/NPjDAXN+tx6+MMRWmy48RjXE+7/5uqcbsZ6cSgZwPVIPPjjEmwRhz\nIc5cR9OqaBaW9yaUSZlixRfABJyFRQ8FhgGfG2N6+r5Eq5OWOP95g11wKxpaAq5Lo9baYuP0xfCX\n2xvAHzhXTDrgLF9+GPD3COVZXrQWPYuWUM5nMc7S8fOBFGAIMNUYc6S1topVamJWVe9NI2NMHWvt\nLg9y2hdrcG45zgLq4ExSN9kY081aO9fTzMowxhjgaeAH639unFj+7JQK4nxi+rNjjDkKpzhIBrYC\n51prF1XRPCzvTcSKBWPMMJz7KlWxQLq19tdQXt9a+06ZzQXGmJ+B34F+OJfAwyrS51PVYQnrAtdl\nXjjA8/H3EvjJzVr7cpnNBcaYtcAkY8zB1tplQSUbXsH+m0bsPQiTKvOz1k7HuXXhNDRmGs7yfdfg\n9Huo7navvRfL70+lfN8TZb8rphtjDsWZETdmOm4CI4EjgONC2DcWPzsBnU81+Owswum70xg4Hxhn\njDneT8FQXtDvTSSvLDyB06/An6XhOpi1dpkxJg9oSwSKBSJ7Pmtx3rwWuCvA5sCcSvfYd4Gez1pf\nHomcZ+gAAANASURBVKWMMYlAE/wvBlbeDJxzbItzNSjSorXoWbSEcj4u1toiY8wcnPeguqnqvdli\nrS3wIJ9ImEloP5QjwhjzHHA60Mdau2YvzWP5swMEfT4usfbZsdYWsefnzWxjTDfgFuD6SpqH5b2J\nWLFgrf0Lp9NFVBhjWgNNcS7vhV0kz8dX6KzFGW0wH8AY0winT8DzETpmQOfjq6gbG2M6lum30B/n\nB/+MqvesoCNOJRuR96c8a22hcdYv6Q98DKWXIPsDz1Sx27RKnj+Jqu8FRk2I5+NijEkAjgJiYohx\nkKYB5Yd7nUwMvDdhdCxR+nzsje8H69lAX2ttbgC7xOxnB0I6n/L7x/pnJwHndlZlwvPeeN2L09e9\n4ECcSyr/wRmucozvUb9Mm0U4a0qAM0zpcZwfpgf5/iFm4VwmqlXdzse3fSfOD+8zgaNxhrf8BtSO\ngfP53Pfv2xXnN5/FwOtlnt/f92/fxbd9CPAvoJPv/TkLWAL8L8p5/wNnRMmlOKM6Rvn+jZv5nh8H\nPFKmfU+gALgNpx/A/ThLpB/h9XsQ4vn8G+dL4WCcYi0LZ0ju4TFwLvV9n4ljcXqm3+rbPtD3/DBg\nbJn2bYBtOKMi2gODfO/ViV6fS4jnc4vvc3EocCTOffRCoF8MnMtIYCPOkMMWZR7JZdqMrS6fnRDP\nJ5Y/Ow8DvX3frUf5/m8VASf4no/I95qnJ13mZMbgXGIt/zi+TJti4FLf35OBL3Eur+TjXI55YfeX\nptePYM+nTOx+nA6BO3B6q7b1+lx8eTUGxuMUPhtxxuvWK/P8QWXPD2cxscnAet+5LPb9h27gQe6D\ngOU4P2Sn4StofM/9D3i1XPvzcQq5nThXeU7x+t8/1PPBWSF2ma/tauAToIPX5+DLrS/OD9Xyn5FX\nfc+PoVxx6dsn23c+vwGXeH0eoZ4PToe533B+AK3HGeVyvBe5V3IulZ2H6/uqOn12QjmfGP/svIzz\nM28nzs/Ar/AVCpF8b7SQlIiIiPhV3edZEBERkQhTsSAiIiJ+qVgQERERv1QsiIiIiF8qFkRERMQv\nFQsiIiLil4oFERER8UvFgoiIiPilYkFERET8UrEgIiIifqlYEBEREb/+H4ePqwMvj4tlAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10447ae80>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Computing with standard error propagation rules: \n",
    "M1 = 40e4\n",
    "dM1 = 2e4\n",
    "M2 = 30e4\n",
    "dM2 = 10e4\n",
    "r = 3.2\n",
    "dr = 1.0\n",
    "G = 6.67384e-11\n",
    "\n",
    "F_oldschool = (G*M1*M2)/r**2\n",
    "dF_oldschool = ( (dM1*G*M2/r**2)**2 + (dM2*G*M1/r**2)**2 + (-2*G*M1*M2*dr/r**3)**2 )**0.5\n",
    "\n",
    "print(\"Resulting Force F = {0:.2e} +- {1:.2e}\".format(F_oldschool, dF_oldschool))\n",
    "\n",
    "\n",
    "number = 1000000 #number of values we want to generate\n",
    "\n",
    "# Ranomly sampling values\n",
    "M1ar = np.random.normal(M1, dM1, number)\n",
    "M2ar = np.random.normal(M2, dM2, number)\n",
    "rar = np.random.normal(r, dr, number)\n",
    "Gar = np.ones(number)*G\n",
    "\n",
    "# Now we compute the force with the formula\n",
    "F = (M1ar*M2ar*Gar)/rar**2\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.mlab as mlab\n",
    "%matplotlib inline\n",
    "x = np.linspace(-2, 3, 10e3)\n",
    "plt.hist(F, range=[-1., 5.], normed=True, bins=50)\n",
    "plt.plot(x,mlab.normpdf(x, F_oldschool, dF_oldschool), linewidth=3)\n",
    "plt.axis((-1.5, 3, 0, 1))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "source": [
    "In this case, which method do you think is more accurate? Why? What do you think are the advantages of using a Monte-Carlo technique?\n",
    "\n",
    "The Monte-Carlo technique is applicable in a much greater variety of circumstances."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Computing with standard error propagation rules: \n",
    "\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "from ipywidgets import interact, interactive\n",
    "from IPython.display import clear_output, display, HTML\n",
    "import numpy as np\n",
    "#from scipy import integrate\n",
    "from matplotlib import pyplot as plt\n",
    "import matplotlib.mlab as mlab\n",
    "M1 = 40e4\n",
    "dM1 = 0.1e4\n",
    "M2 = 30e4\n",
    "dM2 = 0.3e4\n",
    "r = 3.2\n",
    "dr = 0.4\n",
    "G = 6.67384e-11\n",
    "\n",
    "F_oldschool = (G*M1*M2)/r**2\n",
    "dF_oldschool = ( (dM1*G*M2/r**2)**2 + (dM2*G*M1/r**2)**2 + (-2*G*M1*M2*dr/r**3)**2 )**0.5\n",
    "    \n",
    "def MonteCarlo(N=10000, dM1=1, dM2=1, dr=):\n",
    "    # Ranomly sampling values\n",
    "    M1ar = np.random.normal(M1, dM1, N)\n",
    "    M2ar = np.random.normal(M2, dM2, N)\n",
    "    rar = np.random.normal(r, dr, number)\n",
    "    Gar = np.ones(number)*G\n",
    "    \n",
    "    # Now we compute the force with the formula\n",
    "    F = (M1ar*M2ar*Gar)/rar**2\n",
    "    \n",
    "    x = np.linspace(-2, 3, 10e3)\n",
    "    plt.hist(F, normed=True, bins=10)\n",
    "    plt.plot(x,mlab.normpdf(x, F_oldschool, dF_oldschool), linewidth=3)\n",
    "    plt.axis((-1.5, 3, 0, 1))\n",
    "    plt.show()\n",
    "    \n",
    "w = interactive()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}