{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Numerical Methods 1\n", "### [Gerard Gorman](http://www.imperial.ac.uk/people/g.gorman), [Matthew Piggott](http://www.imperial.ac.uk/people/m.d.piggott), [Nicolas Barral](http://www.imperial.ac.uk/people/n.barral)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Revision exercises" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Question 1\n", "\n", "Below is a extract of shot data from a seismic survey.\n", "\n", "| Time (ms) | Shot |\n", "|-----|-----|\n", "| 0 | -0.021373 |\n", "| 4 | -0.024578 |\n", "| 8 | -0.023914 |\n", "| 12 | -0.018227 |\n", "| 16 | -0.00781 |\n", "| 20 | 0.005602 |\n", "| 24 | 0.019264 |\n", "| 28 | 0.030235 |\n", "| 32 | 0.036059 |\n", "| 36 | 0.035334 |\n", "\n", " 1. Calculate the Lagrange polynomial for these points. Plot both the Lagrange polynomial and the raw data points.\n", " 2. The full shot is available in the file [shot.txt](data/shot.txt) (in the data folder) - where the sample interval is 4ms as above. Note that the file only contains one column as you can calculate the time column yourself. Use cubic-polynomial splines to re-interpolate the data for a sample interval of 7.07ms. Plot both the original shot data and the interpolated time series." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEKCAYAAAA1qaOTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNUZ//HPw5oECCKiICigooisEsGlij9AlFaLK2VR\nsVUxWKwWN7R1A61a25/SCtRotfxQK2hRcUG0iApqkQBSWVwwYTWBCEpFCCTw/P64F4whkACTuTPJ\n9/16zStzz5y58+RC8uScc8855u6IiIgcqBpRByAiIlWDEoqIiMSEEoqIiMSEEoqIiMSEEoqIiMSE\nEoqIiMSEEoqIiMSEEoqIiMSEEoqIiMREragDiKdDDjnEW7VqFXUYIiJJZd68eV+7e5Py6lWrhNKq\nVSuys7OjDkNEJKmY2YqK1FOXl4iIxIQSioiIxIQSioiIxES1GkMpS1FREatXr6awsDDqUJJaSkoK\nLVq0oHbt2lGHIiIRqfYJZfXq1TRo0IBWrVphZlGHk5TcnfXr17N69Wpat24ddTgiEpFq3+VVWFhI\n48aNlUwOgJnRuHFjtfJEqrlqn1AAJZMY0DUUkWrf5SUiUhny8vIYMGAAkyZNomnTpj96betWWLEC\nVq2CDRtg/XrYuBGKi2H7dtixA+rVg/r1oUEDaNoUmjcPHunpEX1DFaCEkgDq16/Ppk2bog5DRGJo\n9OjRzJ49m5tvHsW5545jwQJYsACWLIE1a8B9/8572GFw/PHQrh1kZED37tC2LdRIgP4mJZQqyN1x\nd2okwv8wkWomJSWVrVt/GE98+unxPP30eCCFTp22cOaZcMwx0Lo1HHkkHHIING4MDRtC7dpQs2bw\nvs2b4fvvg5ZLfn6QhFatgk8/haVLYeJEGDcuqNuwIfToAX36wFlnQZs2sLMXem8tpZjb+cunOjy6\ndu3qpS1ZsmS3snirV6/ebmVTp071bt26eefOnb1Xr16en5/v7u7r1q3z3r17e7t27fzKK6/0I488\n0gsKCjw3N9ePPfZYv+yyy7xdu3a+fPlyz8zM9K5du3q7du38zjvv3HXuli1b+p133uldunTx9u3b\n+9KlS/d6bnf3iRMn+kknneSdOnXyoUOHenFx8W4xJ8K1FInC11+7jx3r3rOne82aXzkMckhzwGvX\nTvO+fQf78uV5Mf3M7dvdlyxxf+op96uvdm/d2j1o97gfd5z7bbe5z53rnpk5zGvUqOHDhg3b788C\nsr0Cv2PN97fdlYQyMjK89FpeS5cu5fjjjwfghhvg449j+5mdO8Mjj+y9TlldXt988w0HHXQQZsYT\nTzzB0qVL+fOf/8zw4cNp3rw5t912G2+88QZ9+/aloKCATZs2cdRRR/HBBx9w8sknA7BhwwYOPvhg\ntm/fTq9evfjLX/5Cx44dadWqFTfeeCPXXXcd48aNY/78+TzxxBN7PHdBQQG33HILU6ZMoXbt2lx7\n7bWcfPLJXH755T+KueS1FKnqiovhtdfgqafg9dehqCjoejr/fPj002FMnZpFnTp12LZtG9dccw3j\ndjYnKtGXX8Ibb8CLL8KMGanA7ndepqSksGXLln06r5nNc/eM8uqpTyRBrV69mrPPPpsOHTrw0EMP\nsXjxYgBmz57NgAEDADjnnHNo1KjRrve0bNlyVzIBmDx5MieeeCJdunRh8eLFLFmyZNdrF154IQBd\nu3Zl+fLlez33jBkzmDdvHieddBKdO3dmxowZ5OTkVN43L5LACgrg/vvh6KOD5DFnDlx33Q/jI/ff\nD2ZryczM5D//+Q+ZmZnk5+fHJbajj4Zf/xr+/W9YtCiH7t0HUaNGGgCpqWkMHjyY3NzcSvt8jaGU\nUF5LIp6uu+46RowYwc9//nPeeecd7r777nLfU69evV3Pc3Nz+dOf/sTcuXNp1KgRV1xxxY/midSt\nWxeAmjVrUlxcvNfzujtDhgzh/vvv379vRqQK+OoreOgheOwx2LIFevUKfmecdx7UKvWbdMqUKbue\njx07Ns6RBk44oRlduqQzd24hKSkpbN1aSHp6eqWOo6iFkqA2btxI8+bNAZgwYcKu8tNOO43JkycD\n8Oabb/LNN9+U+f7//e9/1KtXj4YNG7J27VqmTZtW7mfu6dy9evXihRdeYN26dUDQlbZiRYVWsxZJ\neitXBn/1H3UU/PWv0L8/LF4ctAIuuGD3ZJJI1q6Nb0spgS9F9bF582ZatGix63jEiBHcfffdXHLJ\nJTRq1IiePXvuaqbeddddDBw4kIkTJ3LKKafQtGlTGjRosNsYTKdOnejSpQtt27bliCOO4LTTTis3\njj2d+5BDDuHee++lT58+7Nixg9q1azN27FhatmwZ2wshkkA2bIB77oHx44PjK66AkSODxJIs4t5S\nqsjIfVV5JOpdXvuisLDQi4qK3N39gw8+8E6dOiXMuZPtWoqUZds29zFj3Bs1cq9RI7iDasWKqKOK\nFhW8y0stlCSzcuVK+vfvz44dO6hTpw6PP/54UpxbJBlMmwa//S189lkwRvLww9ChQ9RRJQ8llCTT\npk0bFixYkHTnFklk69bBb34DkyYFkwKnToVzz/1hcqBUjAblRaTacg9mnB9/PEyZEoyZLFoU3Lml\nZLLv1EIRkWpp1SoYOjSYCHjyyfDEE3DCCVFHldzUQhGRaueFF6BjR3jvPRgzBmbPVjKJBSUUEak2\nvv8erroKLrkkGCtZuDAYO9m5IKMcGCWUaqJ+/fp7ff3bb7+Ny1pDIlFZsAC6doUnn4TbboP33w9W\n/ZXYUULZD3l5efTo0SPms07dnR07dsT0nBWlhCJVWVZWME7y3XfBDPc//CFYKl5iSwllP+zcOGfU\nqFEHfK7ly5dz3HHHcfnll9O+fXtWrVrFsGHDyMjI4IQTTuCuu+4CYO7cubsWdHz55ZdJTU1l27Zt\nFBYWclQZU3dzc3M55ZRT6NChA7///e93lW/atIlevXpx4okn0qFDB15++WUARo4cyZdffknnzp25\n+eab91hPJJkUFsLVV8M118D/+T9BF1fPnlFHVYVVZPZjZT2Ac4DPgGXAyDJerwtMCl+fA7QKy7sB\nH4ePhcAFFfm8A50pn5KS4sBuj5SUlAqfo7Tc3Fw3M//www93la1fv97d3YuLi71Hjx6+cOFCLyoq\n8tatW7u7+4033ugZGRk+e/Zsf+edd3zAgAG7nfe8887zCRMmuLv7o48+umvPlaKiIt+4caO7uxcU\nFPjRRx/tO3bs8NzcXD/hhBN2vX9P9fZGM+UlUXz11VfevfsZ3rlznoP7737nXsYWPlJBVHCmfGQt\nFDOrCYwF+gLtgIFm1q5UtSuBb9z9GOBh4MGwfBGQ4e6dCZLSY2ZW6bdA5+TkMGjQINLSguWg09Ji\nsxx0RZadr1WrFkcffTRLly7lo48+YsSIEbz33nvMmjWL008/fbdzvv/++wwcOBCAyy67bFe5u3P7\n7bfTsWNHevfuzZo1a1i7du1u769oPZFElJk5mjlzZrN48SimTIF779XAezxEOQ+lG7DM3XMAzOw5\noB+wpESdfsDd4fMXgEfNzNx9c4k6KQQthUrXrFkz0tPTKSwMloMuLIzNctAVXXb+jDPOYNq0adSu\nXZvevXtzxRVXsH37dh566KEyz2tlzMx65plnKCgoYN68edSuXZtWrVr9aFn7fa0nkkhSU1N/9P+0\nqGg8F144fr82lZJ9F+UYSnNgVYnj1WFZmXXcvRjYCDQGMLPuZrYY+ATIDF/fjZkNNbNsM8suKCg4\n4KAreznovS07f/rpp/PII49wyimn0KRJE9avX89nn31G+/btdzvPaaedxnPPPQcEyWGnjRs3cuih\nh1K7dm1mzpy5axn6Bg0a8N1335VbTyRRucM11+QAP2wqFateBKmYpJ0p7+5zgBPM7HhggplNc/fd\n/oR29ywgC4ItgA/0cyt7Oei9LTvfvXt31q5dyxlnnAFAx44dyc/PL7MlMmbMGAYNGsSDDz5Iv379\ndpUPHjyY8847jw4dOpCRkUHbtm0BaNy4Maeddhrt27enb9++3HrrrWXWE0lEhYVw5ZXw7LPNaNs2\nnc8/j20vglRMZHvKm9kpwN3ufnZ4fBuAu99fos70sM6H4RhJPtDESwVtZm8Dt7j7jzeML6W8PeXl\nwOhaShTWrw+24p09O9h+96OPLqRZs2YMHTqUrKws8vLyfvSHoOy7iu4pH2ULZS7QxsxaA2uAAcCg\nUnWmAkOAD4GLgbfd3cP3rHL3YjNrCbQFlsctchFJCGvWQJ8+8OWXwUrB/fsDRL/9bnUVWUIJk8Fw\nYDpQE3jS3Reb2SiCW9SmAn8HJprZMmADQdIB+Akw0syKgB3Ate7+dfy/CxGJyuefB8lkwwaYPh16\n9Ig6Iol0DMXdXwdeL1V2Z4nnhcAlZbxvIjAxhnGUOQ4hFRdV16lUTwsWwNlnB8/feQdOPDHScCRU\n7WfKp6SksH79ev1CPADuzvr160lJSYk6FKkG3n03aI2kpgbjJkomiSNp7/KKlRYtWrB69WpicUtx\ndZaSkkKLFi2iDkOquDffhH79oHXr4Ln+yyWWap9QateuTevWraMOQ0TK8cYbwd1cbdvCW29BkyZR\nRySlVfsuLxFJfK+9FrRM2rWDGTOUTBKVEoqIJLSpU+GCC4IdFmfMgMaNo45I9kQJRUQS1ksvwUUX\nBQPvb70FjRpFHZHsTbUfQxGRxPTGG8FExYyMYJ5JenrUEUl51EIRkYTz7rtBN1f79jBtmpJJslBC\nEZGEMmcOnHsuHHVUcGvwQQdFHZFUlBKKiCSMhQvhnHPgsMOCMZNDDok6ItkXSigikhCWLoWzzoIG\nDYK7uQ4/POqIZF8poYhI5HJyoHdvqFED/v1vaNky6ohkf+guLxGJ1Nq1warBhYXBQo/HHht1RLK/\nlFBEJDKbNsHPfgZffQUzZ0KHDlFHJAdCCUVEIrFtWzBp8eOPgwmM3btHHZEcKCUUEYm7HTuCPeDf\nfBP+/vfgNmFJfhqUF5G4u+02ePppGD0afvWrqKORWFFCEZG4GjMG/vhHyMyE3/0u6mgklpRQRCRu\nJk+G3/422Nfk0UdBO29XLUooIhIXH3wAl18Op54Kzz4LNWtGHZHEmhKKiFS6nJxgg6wjjoCXXw72\ng5eqRwlFRCrVt98Gc022bw92XtQGWVWXbhsWkUpTVAQXXwxffhks9qhZ8FVbpC0UMzvHzD4zs2Vm\nNrKM1+ua2aTw9Tlm1iosP8vM5pnZJ+HXnvGOXUT2zh2uvTZY6DErC3r0iDoiqWyRJRQzqwmMBfoC\n7YCBZtauVLUrgW/c/RjgYeDBsPxr4Dx37wAMASbGJ2oRqag//xmeeAJuvx2uuCLqaCQeomyhdAOW\nuXuOu28DngP6larTD5gQPn8B6GVm5u4L3P2rsHwxkGpmdeMStYiU68UX4ZZbgi18R4+OOhqJlygT\nSnNgVYnj1WFZmXXcvRjYCJQe0rsImO/uWyspThGpoLy8PLp27cHgwfl06wb/+EewJL1UD0n9T21m\nJxB0g12zlzpDzSzbzLILCgriF5xINXT77aOZP382NWqM4qWXdHtwdRPlXV5rgCNKHLcIy8qqs9rM\nagENgfUAZtYCeBG43N2/3NOHuHsWkAWQkZHhMYteRHZJTU2lsLBw1/H334+nWbPxpKSksGXLlggj\nk3iKsoUyF2hjZq3NrA4wAJhaqs5UgkF3gIuBt93dzewg4DVgpLu/H7eIRaRMOTk5tGkzCEgDIC0t\njcGDB5ObmxttYBJXkSWUcExkODAdWApMdvfFZjbKzH4eVvs70NjMlgEjgJ23Fg8HjgHuNLOPw8eh\ncf4WRCT0yivN+OKLdKCQlJQUCgsLSU9Pp2nTplGHJnFk7tWnFygjI8Ozs7OjDkOkSpk9G3r2hEaN\nLuTCC5uRmTmUrKws8vLymDJlStThSQyY2Tx3zyi3nhKKiOyvVasgIwMaNoQ5c6BRo6gjkspQ0YSi\npVdEZL9s3hwsQ19YCO++q2QiSigish/c4aqrYMECmDoV2raNOiJJBEooIrLP/vQn+Oc/4Q9/0H7w\n8oOkntgoIvH39tswciRccknwVWQnJRQRqbBVq2DAADjuOHjySW3hKz+mhCIiFbJ1a7C3SWEhTJkC\n9etHHZEkGo2hiEiF3HADfPQR/OtfGoSXsqmFIiLl+sc/4G9/g1tvhQsvjDoaSVRKKCKyV/PnQ2Zm\nMBv+3nujjkYSmRKKiOzR+vVw0UXQpElwm3AtdZLLXui/h4iUaft2GDwYvvoKZs2CQ7X8qpRDCUVE\nynTPPTB9ejB20q1b1NFIMlCXl4jsZtq0YC/4K66AoUOjjkaShRKKiPzIqlVw6aXQsSOMG6fJi1Jx\nSigisktREfziF7BtGzz/vPaEl32jMRQR2eX22+HDD4M7uo49NupoJNmohSIiALzySrCKcGZmsF6X\nyL5SQhERVqyAIUOgSxd4+OGoo5FkpYQiUs1t2wb9+0NxMUyeDCkpUUckyUpjKCLV3K23Bos+Pv88\nHHNM1NFIMlMLRaQae/FFeOQRGD48WJpe5EAooYhUUzk58MtfQkZGMBgvcqAiTShmdo6ZfWZmy8xs\nt81EzayumU0KX59jZq3C8sZmNtPMNpnZo/GOWyTZbd0ajJtAMG5St2608UjVEFlCMbOawFigL9AO\nGGhm7UpVuxL4xt2PAR4GHgzLC4E7gJviFK5IlXLTTTBvHjz1FLRuHXU0UlVE2ULpBixz9xx33wY8\nB/QrVacfMCF8/gLQy8zM3b9399kEiUVE9sHzz8OjjwY7MF5wQdTRSFUSZUJpDqwqcbw6LCuzjrsX\nAxuBxnGJTqQKWrYMrrwyWD34wQfLry+yL6r8oLyZDTWzbDPLLigoiDockcgUFsIllwSbZE2eDHXq\nRB2RVDVRJpQ1wBEljluEZWXWMbNaQENg/b58iLtnuXuGu2c0adLkAMIVSW4jRsDHH8OECdCyZdTR\nSFUUZUKZC7Qxs9ZmVgcYAEwtVWcqMCR8fjHwtrt7HGMUqRImTYLx44PB+PPOizoaqaoimynv7sVm\nNhyYDtQEnnT3xWY2Csh296nA34GJZrYM2ECQdAAws+VAOlDHzM4H+rj7knh/HyKJ7vPP4aqr4NRT\n4Q9/iDoaqcqsIn/wm9lp7v5+eWWJLiMjw7Ozs6MOQyRutmyBk0+GNWtgwQI44ojy3yNSmpnNc/eM\n8upVtMvrrxUsE5EEcv318N//wsSJSiZS+fba5WVmpwCnAk3MbESJl9IJuqlEJEE98ww8/jiMHAl9\n+0YdjVQH5Y2h1AHqh/UalCj/H8EguYgkoE8/hWuugdNPh9Gjo45Gqou9JhR3fxd418z+4e4rzKx+\nWL4pLtGJyD7bvDmYb5KWFmzlW0ubVEicVPS/WgMzWwAcDGBmXwND3H1RpUUmIvtl+HBYvBjeeAOa\nl157QqQSVXRQPgsY4e4t3b0lcGNYJiIJZMKEYMHH3/0O+vSJOhqpbiqaUOq5+8ydB+7+DlCvUiIS\nkX2Wl5dHRkYPMjPzOfNMuPvuqCOS6qiiCSXHzO4ws1bh4/dATmUGJiIVd8cdo5k3bzZmo3j2Waip\nezAlAhUdQ/kVcA8wJTyeFZaJSIRSU1MpLPxhF4ctW8Zz+OHjSUlJYcuWLRFGJtVRhVoo7v6Nu//G\n3U8MH9e7+zeVHZyI7F1OTg7duw8C0gBIS0tj8ODB5ObmRhuYVEsVaqGY2bEEuyO2Kvked+9ZOWGJ\nSEUUFDQjOzsdKCQlJYXCwkLS09Np2rRp1KFJNVTRLq/ngb8BTwDbKy8cEamo777bub/JWi67LJMb\nbhhKVlYWeXl5UYcm1VRFE0qxu4+v1EhEpMLcg5nwy5bB229PoUePoHzs2LHRBibVWnlreR0cPn3F\nzK4FXgS27nzd3TdUYmwisgdZWcEs+HvvZVcyEYlaeS2UeYADFh7fHB7vdFRlBCUiezZ/PvzmN3D2\n2XDbbVFHI/KDvd7l5e6t3f0o4Fagk7u3Bp4CFqLFIUXi7ttvg3GTJk3g6aehRpR7roqUUtH/jr93\n9/+Z2U+AngSD8xpTEYkjd/jVr2DlSpg8GQ45JOqIRH6sogll551dPwMed/fXCJa2F5E4eeQRePFF\nePDBYDtfkURT0YSyxsweA34BvG5mdffhvSJygD78EG65Bc4/H37726ijESlbRZNCf2A6cLa7f0uw\njP3NlRaViOzy9dfQv3+whe9TT4FZ+e8RiUKF5qG4+2Z+WMcLd88DNHtKpJLt2AGXXw7r1sEHH8BB\nB0UdkcieaS83kQT2wAMwbRqMHw9du0YdjcjeaRxEJEHNnAl33AEDBwaz4kUSXaQJxczOMbPPzGyZ\nmY0s4/W6ZjYpfH2OmbUq8dptYflnZnZ2POMWqWz5+UEiOfbYYFa8xk0kGUSWUMysJjAW6Au0Awaa\nWbtS1a4EvnH3Y4CHgQfD97YDBgAnAOcA48LziSS94uIgmfzvf/D881C/ftQRiVRMlC2UbsAyd89x\n923Ac0C/UnX6ARPC5y8AvczMwvLn3H2ru+cCy8LziSS9u++Gd94Jxk3at486GpGKizKhNAdWlThe\nHZaVWcfdi4GNQOMKvlck6UybBvfdB1deCUOGRB2NyL6p8oPyZjbUzLLNLLugoCDqcET2aOVKuPRS\n6NgR/vrXqKMR2XdRJpQ1wBEljluEZWXWMbNaQENgfQXfC4C7Z7l7hrtnNGnSJEahi8TW1q1w8cVQ\nVBSMm6SmRh2RyL6LMqHMBdqYWWszq0MwyD61VJ2pwM6G/8XA2+7uYfmA8C6w1kAb4KM4xS0Sc9df\nD3PnwoQJwZ1dIskosomN7l5sZsMJlnSpCTzp7ovNbBSQ7e5Tgb8DE81sGbCBIOkQ1psMLAGKgV+7\nu7YmlqT01FPw2GMwciRccEHU0YjsPwv+4K8eMjIyPDs7O+owRHaZPz9YOfgnP4E33oBaWrtCEpCZ\nzXP3jPLqVflBeZFEtX49XHQRHHposJ2vkokkO/0XFonA9u0weDB89RXMmhXswCiS7JRQRCIwahRM\nnx6MnXTTlFypItTlJRJnr74aJJRf/hKuvjrqaERiRwlFJI6WLQsmL554Iowdq0UfpWpRQhGJk++/\nDwbha9aEf/1Lkxel6tEYikgcuAddXIsWweuvQ6tWUUckEntKKCJxcP/9wZIqDz0EZ2v3Hqmi1OUl\nUslefRV+//vgNuEbb4w6GpHKo4QiUomWLoVBg6BLF3j8cQ3CS9WmhCJSSb79Fvr1CwbfX3pJg/BS\n9WkMRaQSbN8ebOO7fDm8/TYccUS5bxFJekooIpXg9tuDxR6zsoKFH0WqA3V5icTYs8/CH/8I116r\nmfBSvSihiMTQvHnBfvBnnAGPPBJ1NCLxpYQiEiNr1gSD8IceGsw5qV076ohE4ktjKCIx8P33cN55\nsHEjvP9+kFREqhslFJEDkJeXxy9+MYC0tEksXNiUV16Bjh2jjkokGkooIgdg9OjRzJo1GxjFmDHj\n+OlPo45IJDraU15kP6SmplJYWLhbeUpKClu2bIkgIpHKoz3lRSpRTk4OPXsOAtIASE1NY/DgweTm\n5kYbmEiElFBE9sPGjc2YPTsdKKRu3RS2bi0kPT2dpk2bRh2aSGSUUET20dq1hGMla7n00kzmzPkP\nmZmZ5OfnRx2aSKQiGZQ3s4OBSUArYDnQ392/KaPeEOD34eG97j4hLL8PuBxo5O714xGzCMB33wXJ\nZO1amDVrCt26BeVjx46NNjCRBBBVC2UkMMPd2wAzwuMfCZPOXUB3oBtwl5k1Cl9+JSwTiZuiIrjk\nEli4ECZPZlcyEZFAVAmlHzAhfD4BOL+MOmcDb7n7hrD18hZwDoC7/8fd8+ISqQjBFr5XXw3Tp8Nj\nj8HPfhZ1RCKJJ6qEcliJhJAPHFZGnebAqhLHq8Mykbi74w6YMAHuuSdYq0tEdldpYyhm9m+grFte\nflfywN3dzCptMoyZDQWGAhx55JGV9TFShY0fD/fdF7RQ7rgj6mhEElelJRR3772n18xsrZk1c/c8\nM2sGrCuj2hrgzBLHLYB39iOOLCALgomN+/p+qd5eegmGD4dzz4Vx47SFr8jeRNXlNRUYEj4fArxc\nRp3pQB8zaxQOxvcJy0TiYuZMGDAATjoJnnsOammhIpG9iiqhPACcZWZfAL3DY8wsw8yeAHD3DcBo\nYG74GBWWYWZ/NLPVQJqZrTazuyP4HqQKmzMHfv5zOOYYeO01qFcv6ohEEp/W8hIp5b//hTPPhIMP\nhlmzoFmzqCMSiZbW8hLZD198AX36QFoa/PvfSiYi+0K9wiKhlSuhd2/Yvh3eeQdatYo6IpHkooQi\nQrCUSu/ewY6LM2dC27ZRRySSfJRQpNr7+uugm2vNGnjrLejSJeqIRJKTEopUa19/HbRMPv8cpk6F\nU0+NOiKR5KVBeam2diaTzz6Dl1+Gs86KOiKR5KaEItVS6WTSp0/UEYkkPyUUqXaUTEQqhxKKVCt5\necGkRSUTkdjToLxUGytWQK9ekJ8fLKfSs2fUEYlULUooUi189lnQzbVpUzAD/uSTo45IpOpRQpEq\nb+HC4A4uM3j3XejYMeqIRKomjaFIlTZ7djBmUrcuvPeekolIZVJCkSrrX/8KurkOOyxILMcdF3VE\nIlWbEopUSX/5C1xyCXTtCu+/Dy1bRh2RSNWnhCJVyo4dcPPNcP310K9fMADfuHHUUYlUDxqUlypj\nyxb45S9h0iS49tqglVKzZtRRiVQfSihSJaxZA+efD/PmwQMPwC23BHd1iUj8KKFI0vvooyCZfPcd\nvPRSsBe8iMSfxlAkKeXl5dGjRw/Gj8+nR4/gtuAPPlAyEYmSEookpXvuGc17783m2mtHcdJJQSul\nQ4eooxKp3tTlJUklNTWVwsLCEiXjmTVrPEcemcKWLVsii0tE1EKRJDNhQg516w4C0gBIS0tj8ODB\n5ObmRhuYiCihSHLYvh1Gj4YBA5qRlpaOWSEpKSkUFhaSnp5O06ZNow5RpNqLJKGY2cFm9paZfRF+\nbbSHekPCOl+Y2ZCwLM3MXjOzT81ssZk9EN/oJd5yc6FHD7jzThg0CE4/fS3DhmXyn//8h8zMTPLz\n86MOUUTsEzMAAAALZElEQVQAc/f4f6jZH4EN7v6AmY0EGrn7raXqHAxkAxmAA/OArsBWoLu7zzSz\nOsAM4A/uPq28z83IyPDs7OwYfzdSWdxhwgS47jqoUQPGjoXBgzW/RCTezGyeu2eUVy+qLq9+wITw\n+QTg/DLqnA285e4b3P0b4C3gHHff7O4zAdx9GzAfaBGHmCWOvv46WIvrl78M1uP673/h0kuVTEQS\nWVQJ5TB3zwuf5wOHlVGnObCqxPHqsGwXMzsIOI+glVImMxtqZtlmll1QUHBgUUulc4dnnoHjj4ep\nU+HBB2HGDC3uKJIMKu22YTP7N1DWSOnvSh64u5vZPve7mVkt4J/AX9w9Z0/13D0LyIKgy2tfP0fi\nJzcXhg2D6dOhe3d4/HHNLRFJJpXWQnH33u7evozHy8BaM2sGEH5dV8Yp1gBHlDhuEZbtlAV84e6P\nVNb3sNPOWdka/K0cRUXw5z9D+/bBUvN//WvwVclEJLlE1eU1FRgSPh8CvFxGnelAHzNrFN4F1ics\nw8zuBRoCN8QhVm66aTSzZ89m1KhR8fi4asMdXnstSBw33QS9esGSJTB8uFYJFklGUd3l1RiYDBwJ\nrAD6u/sGM8sAMt39qrDer4Dbw7fd5+5PmVkLgrGVTwnu+AJ41N2fKO9z9/Uur91nZQdSUjQr+0At\nWgQjRsBbb8GxxwYtlJ/9TIPuIokooe/ycvf17t7L3duEXWMbwvLsnckkPH7S3Y8JH0+FZavd3dz9\neHfvHD7KTSb7Iycnh4EDB1G3blpYkgYMpkuXXP75Tygj10g5VqyAoUOhUyeYOxceeQQ++QTOPVfJ\nRCTZaab8XjRr1oyGDdMpKgpmZZsV0r17Ovn5TRk0CJo3h9/+Nuimkb1buRIyM6FNm2Buya9/DcuW\nBTsr1qkTdXQiEgtKKOVYu3YtmZnBrOxhwzI5/PB8li0Lump69w4m251wAnTrFgwmryvr9oJqbNmy\n4M6tY46BJ5+Eq64Kyv7yF23NK1LVRDKGEpXKmCm/bh08/TRMnAgffxwMJp99Nlx2WbA3R1pa+eeo\natzh7bdhzBh49VWoVQuuvBJuuw2OPDLq6ERkX1V0DEUJJYYWLQqSyzPPwOrVUL8+nHceXHAB9O0b\nHFdl334b7Of+6KPBtWjSJOjmGjYMmjWLOjoR2V8JPShfVbVvH+xnvmJF8Bf6L34RdI317w+HHBK0\nWJ56KlhWJJmVnJdTXAyvvx58r02bBgmkRo2ge2vlShg1SslEpLpQC6WSFRcHk/RefDF4rFwZdIud\nfDKcdRb06QMnnRR0CyWLoUOv5YknHuP4469h/fpxrF0bjIcMHAiXXw4ZGbpjS6QqUZdXGaJebdgd\n5s+Hl14KlhfJzg7KGjaEM8+En/wETjsNTjwx2CO9tLy8PAYMGMCkSZPiuv+HO3zxBbRvn0pR0e73\nSteuncKmTVt0t5ZIFaUurwRkFqycO3p0sAd6QQFMnhysqvvJJ3DzzXDqqZCeHiSWG24Iuo6ys2Hz\nZhg9Oj4z9r/5BmbPhocfhosuCrqyjjsOiopySEsbRK1awZ0GqanBbokrV+YqmYiIWiiJJD8fPvwQ\nPvgg6Cb7+GMIJuSnAmW3DN57bwuNG8PBB8NBB5W/ZIk7bN0a3J22ejWsWRN8XbkymE+zaBF89dUP\n9Y86Kmg57XyMGTOMxx/Pok6dOmzbto1rrrmGcePGxfIyiEiCqWgLJYl67qu+pk2DO8IuuCA43r4d\ncnLg3XdzGDPmJpYufYnt2zcTzNi/gKKiP3HKKT+83yxIKikpwZhMrVo/JJjNm2HTpuCxY8fun52a\nGiwZ37t3cHNB+/bBbPbDD/9xvXXrgnk5Q4cOJSsri7y8vN1PJiLVkhJKAqtZM5hZ3qZNM+bNS2fJ\nkmDG/rZthVxxRTq/+U1T1qyBDRt+eKxfH7RAiouDhFRcHCSQ+vWhXr0fvjZpEsz0b94cWrSARo0q\nNpA+ZcqUXc/Hjh1bid+9iCQbJZQksXPGfsmWQadOQStCRCQRaAxFRET2Snd5iYhIXCmhiIhITCih\niIhITCihiIhITCihiIhITCihiIhITCihiIhITFSreShmVgCs2M+3HwIkw04mijP2kiVWxRl7yRJr\nZcfZ0t2blFepWiWUA2Fm2RWZ2BM1xRl7yRKr4oy9ZIk1UeJUl5eIiMSEEoqIiMSEEkrFZUUdQAUp\nzthLllgVZ+wlS6wJEafGUEREJCbUQhERkZhQQimHmZ1jZp+Z2TIzGxl1PHtjZsvN7BMz+9jMEmad\nfjN70szWmdmiEmUHm9lbZvZF+LVRlDGGMZUV591mtia8ph+b2U+jjDGM6Qgzm2lmS8xssZldH5Yn\n4jXdU6wJdV3NLMXMPjKzhWGc94Tlrc1sTvjzP8nM6iRonP8ws9wS17NzJPGpy2vPzKwm8DlwFrAa\nmAsMdPclkQa2B2a2HMhw94S6b97MzgA2Af/P3duHZX8ENrj7A2GibuTutyZgnHcDm9z9T1HGVpKZ\nNQOauft8M2sAzAPOB64g8a7pnmLtTwJdVzMzoJ67bzKz2sBs4HpgBDDF3Z8zs78BC919fALGmQm8\n6u4vRBUbqIVSnm7AMnfPcfdtwHNAv4hjSjru/h6woVRxP2BC+HwCwS+ZSO0hzoTj7nnuPj98/h2w\nFGhOYl7TPcWaUDywKTysHT4c6Ans/CUd+TXdS5wJQQll75oDq0ocryYBfxhKcOBNM5tnZkOjDqYc\nh7l7Xvg8HzgsymDKMdzM/ht2iUXejVSSmbUCugBzSPBrWipWSLDramY1zexjYB3wFvAl8K27F4dV\nEuLnv3Sc7r7zet4XXs+HzaxuFLEpoVQtP3H3E4G+wK/DLpyE50G/a8L8lVXKeOBooDOQB/w52nB+\nYGb1gX8BN7j7/0q+lmjXtIxYE+66uvt2d+8MtCDonWgbcUhlKh2nmbUHbiOI9yTgYCCSrk4llL1b\nAxxR4rhFWJaQ3H1N+HUd8CLBD0WiWhv2r+/sZ18XcTxlcve14Q/wDuBxEuSahv3n/wKecfcpYXFC\nXtOyYk3U6wrg7t8CM4FTgIPMrFb4UkL9/JeI85ywa9HdfSvwFBFdTyWUvZsLtAnv9KgDDACmRhxT\nmcysXjjoiZnVA/oAi/b+rkhNBYaEz4cAL0cYyx7t/AUduoAEuKbhwOzfgaXu/n9LvJRw13RPsSba\ndTWzJmZ2UPg8leBGnKUEv7AvDqtFfk33EOenJf6QMIJxnkiup+7yKkd4O+MjQE3gSXe/L+KQymRm\nRxG0SgBqAc8mSqxm9k/gTIIVUdcCdwEvAZOBIwlWgO7v7pEOiO8hzjMJumUcWA5cU2KcIhJm9hNg\nFvAJsCMsvp1gbCLRrumeYh1IAl1XM+tIMOhek+AP7cnuPir8uXqOoBtpAXBp2ApItDjfBpoABnwM\nZJYYvI9ffEooIiISC+ryEhGRmFBCERGRmFBCERGRmFBCERGRmFBCERGRmFBCEakkZnaQmV0bPj/c\nzCJduE+ksum2YZFKEq5d9erOlYtFqrpa5VcRkf30AHB0uJDfF8Dx7t7ezK4gmM1cD2gD/AmoA1wG\nbAV+6u4bzOxoYCzBhLXNwNXu/mn8vw2RilGXl0jlGQl8GS7kd3Op19oDFxIs5ncfsNnduwAfApeH\ndbKA69y9K3ATMC4uUYvsJ7VQRKIxM9wf5Dsz2wi8EpZ/AnQMV+c9FXg+WJ4JgEiWJBepKCUUkWiU\nXA9qR4njHQQ/lzUI9uKIZCtXkf2hLi+RyvMd0GB/3hjuGZJrZpdAsIqsmXWKZXAisaaEIlJJ3H09\n8L6ZLQIe2o9TDAauNLOFwGK0/bQkON02LCIiMaEWioiIxIQSioiIxIQSioiIxIQSioiIxIQSioiI\nxIQSioiIxIQSioiIxIQSioiIxMT/B2Sgweql0v2AAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import scipy.interpolate\n", "%pylab inline\n", "\n", "time = [0,4,8,12,16,20,24,28,32,36]\n", "shot = [-0.021373, -0.024578, -0.023914, -0.018227, -0.00781, 0.005602, 0.019264, 0.030235, 0.036059, 0.035334]\n", "\n", "# Create the Lagrange polynomial for the given points.\n", "lp = scipy.interpolate.lagrange(time, shot)\n", "\n", "# Evaluate this fuction at a high resolution so that we can get a smooth plot. \n", "time2 = numpy.linspace(0, 36, 200)\n", "pylab.plot(time2, lp(time2), 'b', label='Lagrange')\n", "\n", "# Overlay raw data\n", "pylab.plot(time, shot, 'k*', label='raw data')\n", "\n", "# Add a legend\n", "pylab.xlabel('time')\n", "pylab.ylabel('shot')\n", "pylab.legend(loc='best')\n", "\n", "pylab.show()" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXl4FFW6/z8nSWfvhJAECIuEVZQtQNgVdcBdAee6O4rj\nMHq94+C4/cZt1OvuqHMdvS6D47iPuNxxGRX3lV12EYSEPRAgCZCQfTu/P9KnqO5Udzpr0+H9PE+e\nVFedrjpd1f1+613OKaW1RhAEQRCCISLUHRAEQRDCBxENQRAEIWhENARBEISgEdEQBEEQgkZEQxAE\nQQgaEQ1BEAQhaEQ0BEEQhKAR0RAEQRCCRkRDEARBCJqoUHegrUlLS9OZmZmh7oYgCEJYsWLFikKt\ndXpT7TqdaGRmZrJ8+fJQd0MQBCGsUEptD6adhKcEQRCEoBHREARBEIJGREMQBEEImk6X0xAEoe2o\nqakhLy+PysrKUHdFaCNiY2Pp3bs3LperRe8X0RAEwS95eXm43W4yMzNRSoW6O0Ir0VpTVFREXl4e\n/fr1a9E+JDwlCIJfKisrSU1NFcHoJCilSE1NbZXnKKIhCEJARDA6F629niIaYc7KlStZtmxZqLsh\nCMJRgohGmDNmzBjGjx8f6m4IQliTmJgYcPvBgwd55plnOqg3RzYiGoIghAVaa+rr60NybBGNw4ho\nCIJwxLJt2zaOPfZYrrjiCoYNG8bOnTu59tpryc7OZujQodx9990A/PDDD/zyl78E4P333ycuLo7q\n6moqKyvp379/o/1u3bqViRMnMnz4cO68805rfWlpKVOnTmX06NEMHz6c999/H4Bbb72VzZs3k5WV\nxS233OK33dGAlNwKghAUf/jDH1i9enWb7jMrK4snnngiYJucnBxefvllJkyYAMADDzxA165dqaur\nY+rUqaxdu5ZRo0ZZffv+++8ZNmwYP/zwA7W1tY7h2+uvv55rr72WK664gqefftpaHxsby7vvvktS\nUhKFhYVMmDCB6dOn8/DDD7Nu3TrrGLW1tY7tjoaiAfE0BEE4ounbt68lGABvvfUWo0ePZtSoUfz0\n00+sX7+eqKgoBgwYwIYNG1i2bBk33ngj3333Hd9//z0nnnhio30uXLiQSy65BIDLL7/cWq+15vbb\nb2fEiBFMmzaNXbt2sXfv3kbvD7ZdZ0Q8DUEQgqIpj6C9SEhIsJa3bt3KY489xg8//EBKSgpXXnml\nNeZgypQpzJ8/H5fLxbRp07jyyiupq6vj0Ucfddyvk1fw+uuvU1BQwIoVK3C5XGRmZjqOaQi2XWdE\nPA1BEMKGkpISEhISSE5OZu/evcyfP9/aduKJJ/LEE08wceJE0tPTKSoqYuPGjQwbNqzRfiZPnsy8\nefOABgEwFBcX061bN1wuF19//TXbtzfMFu52uzl06FCT7Y4GxNPoJGitj4p4qnB0M3LkSEaNGsWQ\nIUPo06cPkydPtraNHz+evXv3MmXKFABGjBjBnj17HH8Xf/3rX7n00kt55JFHmDFjhrX+sssu49xz\nz2X48OFkZ2czZMgQAFJTU5k8eTLDhg3jzDPP5I9//KNju6MBpbUOdR/alOzsbH00PYTJ/CAqKyuJ\niYkJcW+EzsaGDRs47rjjQt0NoY1xuq5KqRVa6+ym3hvS8JRS6gyl1EalVK5S6tYA7f5DKaWVUk1+\noKOVioqKUHdBEISjgJCJhlIqEngaOBM4HrhEKXW8Qzs3cD2wtGN7GF7U1taGuguCIBwFhNLTGAfk\naq23aK2rgXnADId29wGPAEdHaUILCdVIWUEQji5CKRq9gJ2213medRZKqdFAH631Rx3ZsXBEREMQ\nhI7giC25VUpFAH8Bbgqi7dVKqeVKqeUFBQXt37kjEBENQRA6glCKxi6gj+11b886gxsYBnyjlNoG\nTAA+cEqGa63naq2ztdbZ6enp7djlIxcRDUEQOoJQisYPwCClVD+lVDRwMfCB2ai1LtZap2mtM7XW\nmcASYLrW+uipp20GIhqCEDoyMzMpLCwMuv1LL73Edddd1yb7ffDBB4M+blsQMtHQWtcC1wGfAhuA\nt7TWPyml7lVKTQ9Vv8IJ+xgbEQ2hsxPKqdGPZI4a0QDQWn+stR6stR6gtX7As+4urfUHDm1PFi/D\nG/sPSH5MQmfEaWr0zz77jIkTJzJ69GguuOACSktLG70vPz+fKVOmkJWVxbBhw/j+++8BHKdVh4Y7\n+ttuu42srCyys7NZuXIlp59+OgMGDOC5554D4JtvvmHKlCmcffbZHHvssfznf/6n4+/utddeY9y4\ncWRlZXHNNddQV1cHwIsvvsjgwYMZN24cCxcudPy8RUVFnHbaaQwdOpTZs2d73RjOnDmTMWPGMHTo\nUObOnQs0TNleUVFBVlYWl112md92bYrWulP9jRkzRh8tVFZWakADOicnJ9TdEToh69evP/zi+uu1\nPumktv27/vqAx9+6datWSunFixdrrbUuKCjQJ554oi4tLdVaa/3www/r//7v/270vscee0zff//9\nWmuta2trdUlJidZa66KiImvdSSedpNesWaO11rpv3776mWee0Vpr/Yc//EEPHz5cl5SU6H379ulu\n3bpprbX++uuvdUxMjN68ebOura3V06ZN02+//bb1/oKCAr1+/Xp9zjnn6Orqaq211tdee61++eWX\n9e7du3WfPn30vn37dFVVlZ40aZL+3e9+16jfv//9763P8+GHH2pAFxQUePW9vLxcDx06VBcWFmqt\ntU5ISPDah792dryuqwdguQ7CxsrcU2GMfUCfeBpCZ8U+NfqSJUtYv369NedUdXU1EydObPSesWPH\nctVVV1FTU8PMmTPJysoCGqZVnzt3LrW1teTn57N+/XpGjBgBwPTpDVHx4cOHU1paitvtxu12ExMT\nw8GDBwEYN26c9VCnSy65hAULFnD++edbx/3yyy9ZsWIFY8eOBRpmaujWrRtLly7l5JNPxhTqXHTR\nRWzatKlRv7/77jv+9a9/AXD22WeTkpJibXvyySd59913Adi5cyc5OTmkpqY22kew7VqKiEYYI+Ep\noUM5AqZG11pz6qmn8sYbb3i1Wbp0Kddccw0A9957L9OnT+e7777jo48+4sorr+TGG2/kxBNP9Dut\nOmDN3RYREeE1j1tERIR1g+Y7+aHva601s2bN4qGHHvJa/95777X04wMNobEvvviCxYsXEx8fz8kn\nn+w4FXuw7VrDETtOQ2gaEQ3haGPChAksXLiQ3NxcAMrKyti0aRPjx49n9erVrF69munTp7N9+3a6\nd+/Ob3/7W2bPns3KlSsDTqseLMuWLWPr1q3U19fz5ptvcsIJJ3htnzp1Ku+88w779u0DYP/+/Wzf\nvp3x48fz7bffUlRURE1NDW+//bbj/qdMmcI///lPAObPn8+BAweAhqnYU1JSiI+P5+eff2bJkiXW\ne1wuFzU1NU22ayvE0whjRDSEo4309HReeuklLrnkEqqqqgC4//77GTx4sFe7b775hkcffRSXy0Vi\nYiKvvPIK/fr18zuterCMHTuW6667jtzcXE455RTOO+88r+3HH388999/P6eddhr19fW4XC6efvpp\nJkyYwD333MPEiRPp0qWLFS7z5e677+aSSy5h6NChTJo0iWOOOQaAM844g+eee47jjjuOY4891utJ\nhldffTUjRoxg9OjR/OMf//Dbrq2QqdHDmKKiItLS0gBYs2aNFZsVhLZCpkY/zDfffMNjjz3Ghx9+\nGOqutJqwnRpdaB3iaQiC0NFIeCqMEdEQhI7j5JNP5uSTTw51N0KOeBphjIiG0BF0thD20U5rr6eI\nRhgjoiG0N7GxsRQVFYlwdBK01hQVFREbG9vifUh4Kowx0xOAiIbQPvTu3Zu8vDyO1kcOdEZiY2Pp\n3bt3i98vohHG2IXCLiCC0Fa4XC769esX6m4IRxASngpjJDwlCEJHI6IRxohoCILQ0YhohDG+olFb\nW0t1dXUIeyQIQmdHRCOM8RWNoUOHek20JgiC0NaIaIQxvqLhNNWyIAhCWyKiEcZITkMQhI5GRCOM\nEdEQBKGjEdEIY2RwnyAIHY2IRhgjnoYgCB2NiEYYI6IhCEJHI6IRxohoCILQ0YhohDH+RENmJBUE\nob0Q0Qhj/IlGbW1to7YFBQXs2rWrQ/olCELnRWa5DWP8iUZ1dTUul8urbbdu3QDxQgRBaB3iaYQx\ngURDEAShPRDRCGNENARB6GhENMIYf4P7qqqq/L5HqqwEQWgNIhphjD9PI9BT/EpLS9u1T4IgdG5E\nNMIYf0IRSDSKi4vbtU+CIHRuRDTCGH+i4RuCsldMSb5DEITWIKIRxgTradhfi2gIgtAaRDTCGH8D\n+kQ0BEFoL0Q0wpiWeBo1NTXt3zFBEDotIhphjHgagiB0NCIaYUywnoZdUEQ0BEFoDSEVDaXUGUqp\njUqpXKXUrQ7bb1RKrVdKrVVKfamU6huKfh6p2MUhWE9DwlOCILSGkImGUioSeBo4EzgeuEQpdbxP\ns1VAttZ6BPAO8OeO7eWRTbAltxKeEgShrQilpzEOyNVab9FaVwPzgBn2Blrrr7XW5Z6XS4DeHdzH\nIxrJaQiC0NGEUjR6ATttr/M86/zxG2C+0wal1NVKqeVKqeUFBQVt2MUjG6meEgShowmLRLhS6ldA\nNvCo03at9VytdbbWOjs9Pb1jOxdCxNMQBKGjCeVDmHYBfWyve3vWeaGUmgbcAZyktfY/fetRiIwI\nFwShowmlp/EDMEgp1U8pFQ1cDHxgb6CUGgX8DZiutd4Xgj4e0bTE05DwlCAIrSFkoqG1rgWuAz4F\nNgBvaa1/Ukrdq5Sa7mn2KJAIvK2UWq2U+sDP7o5K/HkX4mkIgtBehPQZ4Vrrj4GPfdbdZVue1uGd\nCiOCLbm1eyH2ZeHIY9WqVbhcLoYNGxbqrgiCI2GRCBecaUl4yt5u+fLljBkzhpycnHbspRAs9fX1\njB49muHDh7fo/RJ6FDoCEY0wprUlt3fffTcrV65kzZo17dhLIVj27t3rd9v+/ft56qmnvJ6NYmfO\nnDmkpKQEfNSvILQFIhphTGs9jYiIhst/6NCh9upip6aurq5Nn7ke6FG8v/71r5kzZw6rV6923P7U\nU09RVlbGjh072qw/guCEiEYY0xJPQ0Sj7cjMzGT06NFBt6+rq2PAgAG8/vrrjtsDXYd169YBcODA\ngYDH2L59e9D9EYSWIKIRxrS25FYpBUBJSUl7dbFTk5eX16zQ3tq1a9myZQtz5sxx3G73NHzDUEbg\nd+/e3eh99uu7Z8+eoPsjCC1BRCOMaa2nYcpvO8LT+P777/nqq6/a/TgdRUuEdufOhllz0tLSHLfb\nRaOiosJrW3x8PABO0+TY14nXKLQ3IhphjD8xqK+vZ+rUqbzxxhsB2xkj1RGexpQpU5g6dWq7H6ej\nOHjwoLUcbNWSSVL7M+z29b7jaYznUV5eji9270JEQ2hvRDTCGH+eRnFxMV999RWXXnop4H+cRllZ\nGUC7V9w0N1m8atUqbrjhBr+VQkcCdq/AnMemMOfZn0jb9+N7TYyIOB1LREPoSEQ0whh/OY1NmzYB\nYCZv9JfTMAaovUeJB4rVOzFjxgyeeOIJdu1qNBXZEYPdODvd/TthhMA3fGiorKxs1NZgzqGTaNhL\ndUU0hPZGRCOMMaIRERHhZYiMEUlJSQGaDk+1t6dhN4bFxcVNtjefa9u2be3VpVZjF8K2Eg37dfAV\nciMWTscy5zQ2NjZg2W5z+fHHH72unSCAiEZYY4xrZGSklxgY42Pu6v2JRkd5GvakbjChHJfLBTRU\nJx2ptEY0/IXr7KLRHE/DeBc9evQIOlTWFDk5OYwYMYK77767TfYndB5ENMIYu2jYhSGQaNjDU6Hw\nNIIxsKY/bXnX3Na0RDSMOAfjadiX6+rqLLF3EoWSkhJiYmJITExss2u5ePFi4PD4EEEwiGiEMfX1\n9URERKCUcvQ0jKg4eRrV1dVey+2J3dPwLSV1wowfaau75vbALhTN9TTAObfjLzxlX3Y61qFDh3C7\n3cTExLSZaBjvxZ/ACUcvIhphTF1dHREREY1yGsbI+Hoa9jCW/U45GNEoKSnhtttua5EhD+Rp1NXV\nWXe1BjOQ7UgWDbtxboloOOUK/Hka9mV/nkZSUhIxMTFtloMwx5HEuuCLiEYYU19fT2RkZNCeRkxM\njBWeClTe6cSf/vQnHn74Yd56661m9zOQp/HEE08wadIkvv76a2ud6e+RLBp249zckltw9rj8eRpN\niUZ7eBpGCNtybi2hcyCiEcaY8JSvp+EvpxEbG+sYGw/G0zAT4bVkIGAgT2PDhg0A/PDDD9Y6Y1CP\nZNEIVOkUzHuaEg1/y/7CU8bTaGvRkOopwRcRjTAm2JyG2WYXDROeio+PD8rQGAFqySyqgTwNIwz2\nifjCQTQCjanwR1O5CX/7NMtxcXF+w1Pt5WkEG3oTjh5ENMKY5noaMTExjTyNrl27BnWnbOY3amtP\nw0zHYUSjrq7O6v+RbLAClccG8x5/noYpAnBKhKekpPgNT7WXp3EkXwMhNIhohDHNrZ6KjY21chrG\n0wj2wT2tEQ27gfQNd5iBaUY8ArU9krD3rSXhKSdjXFVVRVJSUqO2ZjklJcXxfa31NMrLyxvNn2WO\ncyR7e0JoENEIY/x5Gv6qp5w8jZSUlKAm3DNGvSXVNIFCOb6eRiDR0Fozf/78oEaVtzdVVVUkJydb\ny8G+x+DP03C73Y3amuWuXbtSW1vbSKRa62kkJCQwc+ZMr3XiaQj+ENEIY0zJbXOqp3xzGsGKhvEw\nWiIadgPpa/CMAJj927f7GtZPP/2Us846iz//+c/N7kNbU1lZaRn4tkyEG0/DqXqqa9eugPfdv9aa\nsrIyEhISiI2NbbZomIc2ffzxx17rjVhUVVU5jtVYunQps2fPluqqoxARjTDGXnLr5GkYApXcdunS\npUnRqKqqanKG1kAE8jSMaBgjGmgsw48//gg0fjqd1pqbb76ZBQsWNLtvLaWqqor4+HgiIyOblQg3\nQtOc8JS5nkY07GNs7Enylnga/p70ZxcmJ4E79dRTeeGFF8jJyWnW8YTwR0QjjPEXnjKYpGqgkttg\nPA27d9GSqT0qKiqsvtiNWm1trdUPIxB2wfMVDTPrrf1ZFtBg+B5//HFOPPHEZvetpVRVVRETE9Ms\nQ11VVUWXLl2AloenwNugm3MUGxvbItGwT6tuH6Xe1Ih3038j5MLRg4hGGOMvEW7wFQ3f8FRERASJ\niYlNiobxLqKjo4OaBsSXyspK4uPjiY6O9jJqdq/F7NeIRmRkZKNjmUedFhUVea1ftWpVs/vUWior\nK4mNjSU6OrpZ4ammRCMxMREIHJ5yesKf8TTq6uqaNfWHXTT8zaflJBrmSYK+10Lo/IhohDHBehpG\nKHzDU4mJibhcLrTWAQ2NMe49evRoUWLUGFjfO2HjMcTGxjbyNJKSkhp5GoWFhUDj6dXts+F21IOb\n/HkaX3zxBY899pjf9xjR8BeeiouLaySuzfE07O3nzZvHVVddFfCc2EXDfl7Ly8sDhtIiIyMB7/E1\nwtGBiEYYY/c0nIy+U8mt3dNISEiwpiEP5G0Y0cjIyGiRaFRUVDjG3I2R6tGjR6OcRnJyciPRMAbK\nNzxlf1jT/v37rWWtdbslav0J4amnnsott9ziWKoajKfhJERGSFNTUwHnfIM5v2Y/ALNnz+bFF19k\n9erVfj9HINEwD/Fyuubme9SUaOzatYv77ruvRR6qcGQSlGgopSYHs07oWJryNMw6fyW3xtOAwKJh\njEn37t2pqalxDIUFwp+BtYuGr6eRnJzcyNAYA+XraeTn51vLdtGYPn06/fr1axfhMAbeX3jKKdZf\nXV1tJbqdjGhlZaXjPu3jNMA7jBTI0zDXNNBzSexP/bOHC8vLy0lLS7OW7WitrT74CrgvzzzzDHfd\ndRcvvPBCwHZtxerVq1m2bFmHHOtoJVhP46kg1wkdSFM5DbtoRERE4HK5WuVp9OjRAwhuenM7wXga\nNTU11NXVBQxPGUHwHYxmj6vbk/YffvghO3bsYOvWrc3qbzA4CaHdmNvv4A1GaGJjY5vlabQ0PGWu\n9b59+/x+jj179lhCZs6dGQtiPA1fr6msrMwKeTXlaZhz35LpZ5qL1ppRo0Yxfvz4dp/u/2gmoGgo\npSYqpW4C0pVSN9r+7gEiO6SHgl/8TY1u327+R0ZGEhUV5ZjTgOBEo3v37oCzaHz++efs3LnT8f32\npLFTTsPst7Ky0is8VVVVZRmn2tpaDh06ZIV37IassLCQqKgo4LDhsxszMyliW+Jk4O1CYb+D931P\nXFxci0TDhKcCJcLt7aOjo/32xd7PAQMGAIfPqdmn8TR8RcPukTQlGuY74a+015clS5bQp08ffvrp\np6Da27F/ztzc3Ga/XwiOpjyNaCARiALctr8S4Pz27ZrQFP6mRjf4iobd0zADwnxF45NPPuGpp57y\nSp76ehq+4Yrq6mpOO+00hg8f7tjPYDwN084enoLDd9JGYPr06QN4G86ioiL69esHHBYNuwH3J2at\nwal6yiTqofmiobX2CnnZRbylnoa5/v4qnA4ePMiePXvo37+/137NfyNSvh6f3ZtrKjxlvBx72DAQ\nTz75JHl5ebz++utBtbdjHzPiTzTq6uqsKXHCld27d3dYwYcTAUVDa/2t1vq/gQme/48Dj2ut/6K1\nllE9Icae03CK2zt5GoHCU1przjzzTObMmcMzzzxj7aekpITIyEjLiPiKxsaNG4HGuQZDUzkNu6dh\nD0+ZdXDY6PTu3Rto7GlkZmYChw2a3Ui1x7PGnbwCu5A1VzSMSDjlNKqrq63yaKWUo6fhKxrV1dXW\nPp28gd27d5OSkkJdXR2DBw8GDp9Tc32NSPkTjaSkpCY9jebOWWb2vWXLlqDa27F7M9u2bXNs8/DD\nD9OtW7egB4IWFBQwfPhw3nnnnWb3pz1YsmQJvXr14p577glZH4LNabiVUquAn4CflFIrlFLD2rFf\nQhDYcxoGIwLQ/PCU3dDecccdlrEwT4Yztfm+Bs9eveQkXnZPw24Mi4uLiY+Pt0o7KyoqvMJTcNhg\nGeNkPA1j4Gprazl48GAjT8P+WYLxNFauXMkDDzzgt/T4vffe4/HHH7de20XDfKZAoqG1prq62q9o\nmM/tNPajqqqK6OholFIkJCQ4ehq+4Sl7GyfD/s0331jL5557LtBYNMxNgm9fjQAcc8wxAUWjpqbG\nb/GCP5obzrJj9y7txRF2jPH/7rvvrHW5ubler+188sknrFu3jgcffLDZ/WkPPvroIwDeeOONkPUh\nWNGYC9yote6rte4L3ORZJ4QQu6dhMLF98O9pmOoXX0/D/FAvvfRSiouLrbiyEY24uDigsadhv4t0\nCkMEGqfRpUsXYmNjrXb+wlO+omEMtDmeP08jIyMjKE/jpptu4s4772T+/PmNtmmtOe+887j55pst\nMXDK05g+9erVq5FoGLGOjo4OKBr+qqeMICQmJjqW3Pp6GnYBO3DgAPv27WPRokXWOnMnvnDhQrKz\ns4Hmexp9+vShuLjYb6jEHq4L1tMwotGSxPmePXuIi4ujV69efkXD3OD8/PPP1rpTTjmFk046ybF4\nwQwoDTa81hrKy8uZMWMGb775pt825je5bdu2kM0CHaxoJGitredxaq2/ARLapUdC0DTladTX16O1\npra2lqioKGtbfX29o6dh4s+nnnoqAOvWrQMa7hIDeRr2GLfTj9VfTuPAgQOkpKRYYuQUnjLHMkbY\neBTGwBnD1LdvX+CwcTKx9mHDhgXlaZgfo9OYBvtd7/r166mvr6empsZveGrgwIGNDJBdFJoSDZfL\nFVA07AbY7mmYxHd1dXUj0bjllluYPHky69evBxqMc1paGpMmTcLlcuFyufyKhr/r3bt3b+rq6vxO\nn26+T5mZmUF5GuXl5ezfv5+YmBjy8/O9zsH+/fuZPXs2a9as8fv+vXv30r17dzIyMhy/h4cOHbLC\nZeb7VF9fb91U+D6rHg5/n1sy51ptbS133HEHS5cuDar9l19+yQcffMDVV1/tt40RvZqaGtauXdvs\nPrUFwYrGFqXUn5RSmZ6/O4HmBx2FNqUpT8O0sXsa0GCEKioqGnka5ocxbFhD5NF8QUtKSkhOTrZE\nw9fTsIvGnj172LhxIyNHjuS9994Dgvc0AiXCzY/Xt9LHJHm7detGYmKi1Rfz/7jjjiMvLy9g4lBr\nbRm1zZs3Aw0xdSM89uqrbdu2+Q0lGUM9aNAg8vPzvY7ZHNFwymkY0UhNTfW667Unwo1o1NTUeD1k\n68CBA5aXYcTxwIEDVvgJ8Ap7mevrdrtxuVyN7mjN98R4ff4EwRjogQMHUl5ebuXTtm7d6picN+I+\nefJktNZeYv/888/zwgsvcMMNNzgeCxqEoEePHvTs2dNRNMy1NW3t/wFHQTL7OXDggN+xNQsXLnT8\nfn3zzTc8+OCDnHHGGX77bMfcpJWUlFj7mz17Nr/61a+sNgUFBUye3DBELtCgzfYkWNG4CkgH/uX5\nS/esE0KIfWp0g93TMG18RcP8yH09DWNou3fvTnJystddlj085S/GDQ2i8fLLL7N27VqeffZZq32w\nnoa/nEZ+fj5ut9tKmhsDbTyN1NRU3G63l2jExMTQv39/KisrA4YXCgsLLSO9c+dOioqKGDBgABdf\nfDHgnbPZtWuX30kCjdEdOHAgZWVlXmLaXNHwrZ4ygpCWluYV9rGHp+yehuljRkYGBw4csL4jJtxi\nBNvgJBpmuvVAngY0fJ8KCgo455xzvHIlxtMYOHAg0PA9KSoqon///lx22WX4YkJSEydO9OorHJ5f\nbPny5X5vAPbs2RPQ0zCiMWbMGEss7GEwpxl77fux98dwww03cMIJJziGlIyHcfDgQb/emP2z2JP3\n+/fvZ9euXbzwwgu8/vrr1newoKCA7OxskpKSjmzR0Fof0FrP0VqP9vxdr7WWSWdCjCm5tXsagUTD\nbDN30P48jaSkJK8fnm8iPJCnsXv3buvHYsoeW+Jp2JPj0PDjzcjIICGhISrq62mkpaV5iUZpaSlu\nt9sybIHvbGdgAAAgAElEQVRCVCY84XK52LFjBz/88APQ8PwOOCwaUVFRFBcXe3kavuGp+Ph46w7c\nbmTM5/InGsbIN5XT8BWNyspKXC4XkZGRXqJh+pSRkUFpaaklsuaaBiMa8fHxxMXFOeY0IiIiyMjI\nsPb10ksv8dFHH/GHP/zBamcMs6nOKikpsQydObfQIAivvPKKdY3GjRvn1Vc4PML+0KFDfnNU9vCU\n/UbAYL6PkyZNoqCggPr6essb6tKlC5s2bWq0z/z8fGskvpNofPbZZ0DDPF++2PfnNO5k2bJl9OnT\nxxJae0hz27Zt1vfQvL+iooLS0lK6d+9OVlbWkS0aSqnBSqm5SqnPlFJfmb/WHlwpdYZSaqNSKlcp\ndavD9hil1Jue7UuVUpmtPWZnoqmcBgT2NPyJRmJiImlpaZZBbioRfujQIbp160ZSUhJ5eXnWlzkv\nL8+adqQ5OQ2Xy2UJlDFYe/bsoUePHo1EI5CnkZiYaIlGoGS4EYWJEyeyc+dOr3LP6upqdu/eTbdu\n3UhPT+fgwYN+DXxpaSmJiYmWMX3zzTctY2IfcNeaRLj9ukCDqBrRdRINMwbGGGDjcQUSDfM/Pj7e\nazJJw6FDh3C73ZY3WFxcbIW/Nm7caIWh9u7di8vl4phjjgEavkf28RPmc4wePZpZs2ZZBtgk5o2R\nrqqqYuPGjdbU9yYvY6empoaCggJ69uxpnX/fYoTNmzeTnp7OwIEDqauro6ioyDofEyZMICcnh2XL\nltG9e3eef/5567yNGTMG8PY4zXkw3xW7gTfk5ORYn33t2rV8+eWX9O3b1xLMuXPnsmvXLp5++mmg\n4TtuRuFv27bNS3Q2bdpkCVx6ejpjxoxh1apVIZnTK9jw1NvAKuBO4BbbX4tRSkUCTwNnAscDlyil\njvdp9hvggNZ6IPA/wCOtOWZno6U5DeNpOIWnEhISrDEZvqJhT4TfdNNNVlLaPKO6T58+LFmyhP37\n93PcccdRXV1tuf++d+X19fUUFxc38jTsU22Ad3gqIyPDEhi7aMTFxVmlu+aO2tfTyMvLY+3atY5h\nCyMoEydOpLKy0mvuovz8fHbt2kWvXr3o0qWLl2g4eRoJCQn07NkTgHvuuYdLLrnEK+zWFjmN8vJy\nS7iNFweBRcNgrv3Bgwctow/+PQ2n8JS53kZ0iouLraRsZWWlFebZu3cv3bp18xIX+936pk2bvPIh\nb775puUpREdHW9fq559/pq6ujgsvvBBwFo29e/eitSYjI8MSDd9rnZuby4ABA6wQ5969e71Eo7i4\nmDvvvJN9+/bx0EMPcejQIcrKyhg9ejRw2BMwISWTA5k6dSq7d+9uVPywadMmTjvtNBITE1mzZg13\n3XUXO3bssGZBNt8zkzPbs2cPEyZMABpEIzc3l7S0NFwuF7t27fISjTPPPJPKysqQjB+JaroJALVa\n62fb+NjjgFyt9RYApdQ8YAZg/0bMAO7xLL8D/K9SSul2GA5ZuX8/P3mSbGbn1kGUavii2P+bNkqB\n1t7/HfZh73Cz9m/7r4Ga2lqmnHgi7qQkpuTnU1NbS1VVFQM92/tUVLDD9p6I995j9PbtdC0vZ+C6\ndcwAEj7/nBlA5po1uLdsYSaQ+t13DP7pJ/4jIgLefZdphw6xOi+P2rfe4oyKCsbl5RH78cecB/Rd\nvpylb7/NaKDitdcYvmkTverqSI6KYvXy5fQFzunfnw83bKDy1Vc5Hxi2YQM9CgspLC+Ht9+morSU\n87Vm/PbtJH/6KRcAvRcuRO/cyYVak/bll1wIdPv6a6ipYfLOnfxi0CAi3nqLWTExDFq+HN54gwHL\nlvGb+Hh44w3OKi7m4IED8M9/kp2Tw/DKSjK+/prLlKLqxRd55Npr6Zaezu9+9zs++eQTzjn7bDIz\nM0mdP5/LleKX5eXkAV0/+giTeqz6+9/J+vFHTklJoaymhpj164l7+20uBwYtXkz0tm0UVlbCyy8z\ndv16+tbUkPntt1xhu957HnkEd309VwD9FyygevNmKktL4eWXrTap69ZxBdDnq684ads20g4csLaf\nvH17w/fi5ZeZuGkTVwAVzz1HfGoqY9atI66uDl5+GXdJCZcD/b7/noSEBC4HfpGXh/1+e9j69fDK\nK5xdVMQv8vLglVcAmG7i7q+8wsBFi7gciH/nHS6oqCBt0yarHcDodevoUV9Pj88+43Kg64cfcsLW\nrfx+9GhWrFzJwaeegokTGbpqFb2joui/cCG/AhLffZf+P/5ondu6l15iT3o6v7L1r198POq117jW\n7abfwoXw6quUet4/89Ah1icmkvzBB5CWBkpZf5WbN3MRMGbzZroUFfEfAP/3f5CXZ7Xps3IlY8aO\nZUhODmcD9f/+N13XruVMpThDKRYC+vPPmQrorVvZ8ve/czJwilIsj4wkbvlytr76Kn+44QZ+e/XV\nVFZVMR74fyedRMmXX7J53jx6nHACRERQXFJC76IiTnC7qezfn9WvvkpxcTHHA0Xff8+BRYuo/vFH\nBgF60ybqN20iPj+fSaefzpaEBErWrePAzz8zKTOTXfn5HNq8mQM5OaQAGbGxZGdnM3HYMB695x4u\nmTmTqOhovv72W+KSkhg/frxX9KGtUYHsr1Kqq2dxDrAPeBew4gta6xYXLyulzgfO0FrP9ry+HBiv\ntb7O1madp02e5/VmT5tCn31dDVwNcMwxx4xpycCgwg0bSDve19ERBEEID3JjYrh42DCWL1/eovcr\npVZorbObatdUeGoFsByYRUM4apHntfk7ItBaz9VaZ2uts01MsLl0GTCAnCVLyFm8mJzFi8n1/G1Z\nsqThb/Fiti5ZwtbFi9nm+dtu/hYtYsfixexYtIidixezc9Eidi5aRJ7nb9eiRexauJDdixaxe+FC\ndi9cSL7nb4/5W7CAvQsXsnfBAvaZv++/p2DBAgq+/57CBQso/P577rrgAkbHxlK7di1XjBnDlWPH\nctHQoQwDhgEXHX88w4ERnr+iL7/khlNO4YKBA/n4gQcYBbw4Zw5ZwMZ589jyf/9HFvDJww9z+1ln\ncU7v3rB6Nf+YM4eRQM477zAC+OC++2DNGk5MSuKhiy9mODAcWDx3LhcOGcLvTjyRZ//rvxgGnNGr\nFxvfeYdhwN+vv56hwL8feojnrruOoUDNqlX8/M47HA988de/UvbDDxwPvHDjjdw+cyZTe/bkwMKF\nHAe8fscdbPv4Y4YA7z30EGzYwLRevfh/554LP//M+cOGMfuEE+Dnn7nz/POZ0q0bbNzI2QMHct2p\np8LGjVyUlcVg4OLRo7lg5EiGRUfz7A03MAiYd999XD5hAheOGkXRkiUMBAYC1599NgOA526+mQHA\nU9dfz/Vnn83UzEyW/vOf9AcWvfoqz91yC/2Aip9+4pdZWVxxwgmwZQuV69fz88cf0w94+5FHWPjq\nq/QDls2bx7O295SsXs2BFSv48Mkn6Qds/vxz7rniCrKSkmDLFtiyhXOOP57Zv/gFbNnCsnnz6Acs\neOUV2LKFq6dN46whQ2DLFirWraM/8Nwtt/DP++6z+tgf6A+cP2oUUzMzKViyhP7Aa/feC5s3w+bN\n3DRzJif27AmbN3PvrFmMdLth82YuHT+ei8eOpXjlSj5+6il0bi7/kZXF5ZMmoXNzGRQRwRkDBzIA\nWPuvf3HOkCFcNWUK5OZyQo8e3HLeeexZsIABwBv3388VkyZxflYWU/v25fdnnsntF17ImORknv/j\nHxkAPP/HP0JuLteeeipnDBwIublcNWUKZx97LOTmcufFFzM6KQm9aRNs2sRHf/kLl40Zw/9cfTXH\nAQXffEPt6tUMB5655hr2fvYZB7/+mpV//zujgUVPPsm+jz5iLPCvP/6RO089lQt790YvWMB/dO/O\nCcCet95iCvDrfv04Gdj6wgtcO3gwt40dy/1TpnAacDpwllLcMWoUfPQRs1JTeXjyZM4FpgN/GjaM\n84D8//1ftvz5z1wSFcXcadMo+8c/uEQpLouI4BJg4bXX8ivg29/8hiuAxddcwxMjR/LHtDRmA/PP\nO4+5o0Zxb/fufDlzJnOAigcfhL/8hYr77uMW4KvTTqPsjjv4c1UVF110UYvsX3MIGJ7SWvcDUEpd\nCHyitS5RSv0JGA3c18pj7wL62F739qxzapOnlIoCkoF2eb5kVHQ0g8aPb49dtykDt2xh1dtvsyUm\nhi2emPohGuZ3AXAnJbHO1r5qyBByExLY5XZzaMAAVgOb4uNZA0SOGUM9sAYo6NmTTXFxbEtOhpEj\nqRs2jLXAT5GR/Ajo4cNhxAi2ut38qJR1jN1du7KqtpbsXr3IPOccfnrmGfqPHk3ihAn8BKyurWU9\nUD1oECUREawHqgYOJP/AATYAUSNGEJOVxQZgd5cubIuLY2dCAq4RI/gZ2J2czI6EBDYCCWPGwJAh\n7E1JIScyEo49lhWlpZyQlQXHHktZ796sqayEwYNZV11Nes+eMHgwcSNHkrN6NddfdRVXXXUVZWVl\npKam8r+ffspT8+dTUFTEqFGj6DpuHKaSv/dJJ7Hn66/5fMsWtgDxw4dTXVXFj8uWcbBrV7YC9O9P\nZWEh24DKjAxy6+oYmJ4O/foRC/Tq3ZttwPaICBK6dmUboDMzqSsoYBtQ6olNr1y5kmeffZZtgGvw\nYErT08mtqwNPzmibUsR06QL9+pFQXs42YFd0NPTrx/aICIqTk6FfP1y1tWwFCtxuqpOS2Ap0GT0a\nMzH8KaNG8d2//83+Ll3YCkQOHAieyQrLe/RgY00N9O/PTpeLgsRE6N+fAykpFBUVccuzz/L888+z\ncOFCNtbUMLh7d9SAAZSkpfGpJ7nd9xe/IHXcOD75/HNq+/ZlaWEhU4YMIWHECLYAu2JiWFdRQc9e\nvYjt2ZMFu3cTuW8fXceO5ZI77yRi0KCGUtzYWCIGDWLZ8uUwYACfb9nCSSedBAMG0G3iRFbNm8de\nt5v09HSufPBBCgsL+eeKFURERND1hBOIjIxkX7dufL53LzfPmMHAgQP5zW9+wyqg7y9/SVJqKsuB\nDW43PyhFca9eqMmT+cvSpdTW1tJjwADyBw7ke8/nSpoxg7x332Xprl2U791L2uTJLFy4ELRm1Jln\nwllnsW/sWF755BPrd/fvdeuIjo4m/ZpriIqKYu5//qc1d1jO00+zYsUKjjvuOLr87ne8/uyzxCnF\nq8CvL7qI3KgoKzl+zhVXsO6rr3jllVcoGzyYv0VH89dbbwWliAPee/lltiUlkXDuuTz/wAN8MGRI\n0LakpQSbCL/TIxgnAL8A/g60NsfxAzBIKdVPKRUNXAx84NPmAxq8HGiYVfer9shnhBMmsbt79+5W\nl9z6JsLtVTpm4JepbTdlh/Hx8V7JzOLiYg4dOkRSUhKnn346f/vb33jxxRet95sks32ai+rqaqsP\nXbp0ISoqioiICGuivejoaK9EuElomgSnfSqNgoICq+IkMTGR0tJSa5oUU7b7+OOPM3fuXK655hri\n4uJIS0tDKcWll17KokWLyMnJoXfv3l5x4CFDhtC7d28rWdmrVy+Sk5MbJcJ9R2GbZ3zbz29RUZFX\notv0a9euXaxYsQKtNStXrrTOb6DqKXNeTZGCPRFuZjy2j9Mw5b/QMP7mwIED1pQsgRLhpujBVE+Z\nSqTVq1db1xuwkv59+vQhOTmZkSNHkp+fz6JFi6itrWXQoEFWxduhQ4coKioiNTWVwYMHs379en78\n8UdGjRpFYmIiv/nNb6zPYsaX7Nixg7y8PLKysoCGclloGAB4zTXXUFhYaH3/MzMzrcfQZmRk8N57\n71FRUcGPP/7IU089Ra9evejVqxexsbEkJCRY1VNm5Hvfvn2twaPmeNHR0XTt2pXu3buzZcsWcnJy\nmDp1KhdddBEul8say2OS5dDw8C+AY4891io+cbvd1vfLvOeyyy6zZjIwo9F79uxpTYkDDd/Dnj17\nUlxczPbt20lPT/f6nmZlZbFq1Sqr+MD0vz0JVjTMLG5nA89rrT+iYdr0FqO1rgWuAz4FNgBvaa1/\nUkrdq5Sa7mn2ApCqlMoFbgQaleUebZgfqREN35Jb3+qp5pTcmifHwWHjZEokjWjExcU1Eg1TTRMR\nEcHVV19NamoqsbGxxMfHW2WKvnMjGcNl9muqkIyBjIqKIioqisrKSqsqxVQCGQNXUVFBWVmZl2jU\n19dTWVlpldyaz/Lb3/620bkxE/XBYeNqyiuzs7Pp3bu3JXo9e/akS5cu1NTUWIJnphExn8lUTxmU\nUiQlJVFSUuIoGvY6eyNOZjoQM+uw2bc5jjFwdtEwFWVKKeu95nhxcXG89957fPzxx3Tp0oXa2lrr\n+vmW3JaXl1NfX095ebn1OUyll6mo2r9/v3W9zXkBrGnxzznnHJRSzJgxA2gwehEREVY5tBGNQYMG\nUV1dTXV1NeMdPHxTqmom5jNiMWbMGG699Vaqqqp44YUX6N69u1VBdLwtJ2luMFJSUnC5XGzevNmq\nTAKs6sCioiLrnNoZNWoUgHWDkZGRQXFxMfX19YwcOZJXX32VHTt2WJ/77LPPBuCkk07i5ptvJjo6\nmjlz5jTaL8CNN97IsmXLuO2220hMTCQ1NdUax5GRkWEJFkD//v2tc7xmzRp8w+9ZWVls3ryZZcuW\nERkZaQ2kbE+CrZ7apZT6G3Aq8IhSKoY2eL641vpj4GOfdXfZliuBC1p7nM5Et27dgIbRtk4lt8GM\n0zBGLz4+3rojNYbG3OmZH5LxNIyBiY+P9xrHsH//fioqKiwjYic1NdUquU1ISPAysHZPAxqLBhy+\ny83Pz8flcllClpCQQEFBgVWCaB4WZESiqKiImpoar7t+J+zP/xg7dizQ8LS/vLw8MjIyLK8ODnsa\ncLj0silPw/S1rKzMSzRMG7tomDLSuLg4LyE3EyKacxIdHY3b7bZEwz5Ow2yvrq62vLeoqCjLgPvO\nIOvraZj9lZWVNfI0zBieoqIia5wGHJ7za9CgQUDDQL7bb7+dBx54gPj4eGvMhdvtpqCgwAoNnnLK\nKdZxnabZMAMCX3jhBaKjo607eaUUDz30ELfddhuPPvoo06dPZ8yYMbz11ltMmzbNer/5rsyaNYt1\n69bxxRdfWHOqme2FhYUUFRVZ3x87I0eOBA5P+mmeOQIwYsQIXC6XVznzpEmT+OqrrxgxYgSpqamU\nlZU1ukkxREREWN83aPCQioqKrLLxSZMmMXToUCZOnEh0dLQlGj///DOnnXaa176MuM2bN49BgwZZ\n38f2JFjRuBA4A3hMa31QKZVBK8dpCC3DhAVKSkqaPbjPHp6Ki4vzWmc8DfNj8xeeiouL85oMz9yJ\nO4lGWlqaNf1DYmJiI0/D3IGCs2iYu9zCwkJ69OhhfU4TnrLXrZv1cHhQmN0oOqGU4vnnn2fZsmXW\n3W6PHj0sY2BEIz4+nq5du3qJEnh7T+Xl5VRUVDiKRmlpaUBPwxj62NhYIiIivIQoOjramobFfl6d\nPA37viIjI62+GYxABxKNsrIyr/CUGRFuRH7Xrl3U1dVZn2HWrFl8+umnXpPs3XvvvRxzzDEcd9xx\n1vfL7XZb02SkpqZy3HHH8e2333LMMcd4eWcG4zXk5ORwzjnneAkjNPwO7rvvcFr1ggu87y1vvvlm\nSkpKuPHGGzl48CDz58/n17/+tdc5zM/Pp6SkxFE0Jk+ezKRJk3j44YeBhlCTwS4gdowQQmOPPxCZ\nmZmsWLGCjIwMlFLEx8ezbt06y9Ps1auX1dbJ04CGMO0JJ5wQ9DFbQ1CfTGtdTsOcU+Z1PuA897DQ\nrkRFRZGQkOAlGk6D+8yDmfyFp4xx881pmB+nEY0tW7YQGRlp/bCNMTEY0TBiZsd3Qjxf0UhOTrb6\nbkSjsrLSEii7p2HCDWZfZWVl1mhwf6Lh1CdfZs+ezezZsx23GdGIjY1FKWUZSiNWdtEwg8R8RcMI\nnMlRREdHW/1avXo1aWlpHHPMMaxcudI6t/aJB4FGomEfeOnP04iIiPArGsZ4BxIN+zUoLy+3rpN5\nr/kM48eP95oIELDClHaSkpIsD9UY6SlTpuCPJM94g6VLl/Lb3/7Wbzt/ZGVl8cEHDSnSPn36NHqq\nZGpqqjWK3Uk0unTp0pDs9jBu3Dj69OnDaaed5vV7awtMDsP+HQesmyTjacDhSIPBDISsrq72yqu0\nJ2376YUOwcTJA3kaxmD4C08ZI+EvpxETE2O1SUlJsY5hF4309PSAnkYg0Th48KBlmMzxzMhp3/CU\nmULEvq/S0lK/nobJozTlaTSFMWrmDs5XNMzobfAvGk7hKXNeiouLGThwoBW/N/u3exq1tbXU1tZ6\nCYNdNOyJcPNekwj3vTu3exrmSYD2fkJjTyMxMdGashwOi4bT9Q6E2+22xNzJSDvx7rvv8sEHH3jl\nntoKE0Iyy00RFRVFbm4uc+e2/WOEzPX3l8S23/zYk+TQICw33XQTycnJVoK9vRHRCEOSkpKspJy/\nnIa5M/UXnnISDbunAYfzGnbjbr/j7du3rxUn9xeeMjh5Gr6iYcJTpg8mPOXP0/CX02gr0Rg6dCif\nfPIJr732mtdnNB5OMJ5GINGABkNhjIb5HE6z1drPe9euXZsMT9nF12AXjeTkZK+bjUCiAVjzSZlz\n3hLRMAQrGhkZGZx77rntMrrZfg2C7U90dHSbexkAF110ETNnzuT3v/+943b75zc5JDsPPPAABQUF\nHZIEBxGNsMR4Gk5Toxtvwhjeuro66yFMZlt5ebllDJRSREZGNvI0gEb5DcArQZqenm7drTUVnrI/\nJMiIhr16JzY21jERXlJSQmFhoZdoJCYmWs9biIyMtPbjKxrBhKea4vTTT7cMnl00jOdm+mqMeFOi\nER0d7TUh44ABA6zKLXN97KJh5n2yC0NSUpI1MaOvR2Ee4BRINExo0Lef0CAa9iowe77BlLPaz0Ww\n2K9FsEa6PbH3IdT96d69O++++65VueeEKeN1ylsopRrlMtsTEY0wJD4+noqKCsdxGsbwBApPgbcx\ncLlcjp6GMfr2RJwxXl26dPEyPE2Fp+wx9qY8DbtomKStr6cBDQ/zSUtLsz5/oERvW2APT9mrmeCw\np+Gb1LWLhnnONxyeKt3uaZiyVmMAAolGSUmJY+gqGE8DGp8bI2JlZWXWQ7fAWwTt4ZPmCrL92MGE\ng9obe7gz1KIRDK+99ho5OTlHxLkT0QhDTKzfN6cRERFh3Q3aPQ0n0bAbg6ioKEdPw4Sn7Ik4Y1yS\nk5ObJRpAkzkNkwi3V09t3downtk3pwEN8XV7NYn58ZvEbFuLhjlnJSUl1vkN1tMwlVAGc1eZlZVl\nhRWmTp0KOHsadmFwu91UVlZaVWz+chq+ohEdHW1dP7sRN/2EBkGsq6uzRMH+eezhD6exDYEwCdzY\n2NgOvSv2h10Aw0E03G53h4WfmiL4ujDhiMFXNMyddlOiYTciwXgaTiV/5o43ISHBy/A43Xn6/hib\n8jTMIDi7p1FfXw84exqbN2/2qnc3U3kHquhqDXYD2lzR8L3zf+mll1i+fDkjRowAGgb3GSGxV0+Z\na2D3NHwT8k45DXsf7XTp0oXy8nK/4Skz+t6cO/v3xG60mnvHa6YjN9cz1JhxJdC4TF0IjHgaYUgg\nT8NewgreomE3Inbj5nK5KC8vR2vtZdhuvPFGYmNjOf/886115k7V5XIF7WkYL8Hs2zz9zi46/sJT\nBt+cBjRMS2GfJkMpZQmVfeBdWxEVFWUZaN9Kp0CJ8OrqasrKyrzO7ZAhQ7ye/Tx27Fjr2jWVCHcq\n/TXYR4T7ehpw2MPwLe/0Jxr2a2Q3tM314oy3eqTMAuR2u/nkk0/49ttvQ92VsEM8jTDEDLhSSnkJ\nhVLKEhBjMOwPYbIbF19PwyRW7W0mTpzY6AE8xuuora31MhxOBmrkyJFceuml1ohf08bMY+TraVRU\nVFBbW2u1sxsse326XaDso7ahwbvJy8tr89CU/dj2EfC+noZTTgMaRMXpHDlhLxgwOImGeQa3r2iU\nlpZSV1fXKAQFh0XfHnK099NXNOzn3T69RXOriEwC9/LLL2/W+9qT008/PdRdCEtENMIQ8zQ1M4LY\n7mmYOzmn8FQg0TDx8aYMmzHSvXr1ajTQyBeXy8Xrr79uvW5KNMycWKafxmvo2rWrl9dgP66TaEDj\nmH1b4Xa72bdvn+VRmH6ZMlzfkJhdNIL1fMznr6qqssI5volw8C8aZnyH07U0novveTMDGI1oGNG1\nn2szJUhLSE5OZufOnWGRPxACI+GpMMSEp0zJrT2nEUg07EbENzzl5Gk4MX78eO655x6ee+45r9BQ\nMJjjm7mbfEtuzfxGpp1JctsNJngnxX2Nn1PFV1vimyA2fS0sLCQ6OrqRoW6Jp2Gf4ddfIhz8h6f8\nDe4DOO+88wDvKS+gwUtNSEgImNOIjY3lT3/6E++//35Qn8OX3r17N/n9Eo58xNMIQ4xomOd52z0N\nc2fqJBr2hF9LPY3IyEjuvvtu4PCgo2CnL/AVDafqKXs7Ixq+d+j20JM/T6O9RMPs3zc8VVdX5/V5\nDEZc9u/fb5XWNoVdNMyEec1NhPsTjbvuuosrr7zSer67HSfRAFi6dKklxvfee29Qn0HovIinEYbE\nxsZaIQhfT8OIht2YmcF99kGA/jyNYO+GoeGO/+OPP+brr78Oqr0x/v5Ew3f55JNPBg4PbDLYP8fQ\noUO9tplSypY+wbEpfEXDLsRO1Vr2SQ5b42k0JxFuxmk4iUZUVJTfCfcSEhIs4bZ/lnHjxnXIcxqE\n8EA8jTDEGIOKiopGOY1AnoYdf55Gc8MHZ555ZtBtzfMeghWNnj17kpub28ibANiwYQN79+5tVC55\n1VVXsX37du68885mfY5gMQbb9Ml8purq6oCi4S/H4IRdNMykhcEmws2IcKdxGk1h/040d8S3cPQg\nohGGGCNhZh+1V08FKxq+nkaw4anWYqp7wDunYZ8I0d4Hf3e4Q4YMYYjDoy2Tk5N54okn2qq7jTAi\nZRuM77oAABCjSURBVK8eMp6fk2jYDXGwiXAjEBUVFZZoNDen4c/TCIR96pD2/h4I4YuEp8IQYwyc\nxmkY0TAGKljRMFNYtLexMPtPTk72GqHuTzSONG655RZGjhzJRRddZK0z5zKQpwHBe3FNhafM7LrG\n0/DNaVRWVlrP52gOpq++YzgEwY6IRhhiNwb+chrmjtifaJgRuva20P4G29zN+k5DYReNI7nCpm/f\nvqxevdprdLQRi6ZEw/dZJP4wgm+mComJiWl0/cyT8KCxp9HSGwBTPOA7hkMQ7IhohCH2O0tfT8M8\nbMaUpfoTDfvdpF002ttgG8MaSDSOZE/DCd9ZcO20RDSUUlaFXHFxsWNVltvttqYr9xUNQ3OvpZl2\n217SLAi+iGiEIXZj4Ftye/3117Nu3TrrQfe+omEqZ+wGriM9jc4sGk6ehhk0B8GLhmlbVlbGwYMH\nHQcq2o/lTzSaex4nTJgAOE+/LQgGSYSHIb7VMvbwVExMDEOHDrWesOYrGosXL2bHjh1+nyseKtGw\ne0++8zcd6fjOR2VHKUViYiKHDh1qlmi43W4OHTrkVzTsx2orT+PMM89k1apVjcqYBcGOiEYYYjcG\n9vEXdiEwIlFXV0d9fb2VdO7WrVuj6T86MjxlBuYF8jTaenba9sY8iMrfCPmYmJgWi0ZxcbHjNORG\nNOLi4rwquVojGuA9v5QgOCHhqTAkkKdhMKJhpsn2zWnY6UhPw3gRvobQ7l2012SD7YURYX/zKpnr\n0hzRMA9a2rt3r+N+/eVRWhOeEoRgEE8jDPHnaYSDaBh8RcP+fIZwC0898cQTdO/enVNPPdVxe0tE\nw+12s3//fnbv3u04uDEY0TiSq9CE8EU8jTDEbtijoqLCytOYM2cO3bt393pGB3jfpTd32u1Qk5GR\nwZNPPun33JkBes0Ju7ndbnJzc6mpqXGcR8tfmW9HhhqFoxPxNMIQ3/BUW3kaHfGA+hEjRljTiNjx\nfQ5FZ8LM59ScSRTT09M5cOAAAMcee2yj7cbD8PXKJDwltDciGmFIMJ6GWW6OaMTExHgl0zsSpRQn\nn3wy48aNC8nx25OuXbtSWlrqNaCyKewz4o4aNarRdiMa9lH14H1DIZ6G0B6IaIQhvqIRqHqquaIR\nSoKdLTfceP/993nssce8RpE3hZm6vFu3bo6JcJN89w3l+T7/QhDamvAKHgtA86qnTDw9HESjs5KV\nlcVrr73WyCsIxPTp05k1axZvv/2243YjQL7l03bRkOsptAfiaYQh/jwNJ9Ew8XQRjfAiLi6Ol156\nye/2SZMm8de//pXLLrvMa324zOElhC8iGmGIXRyCrZ4KdJdrRCOQsAhHFkop5syZ02i9hKeE9kbC\nU2GOv+qpliTCw63UVWiMXTScph8RhNYiViLM8edpKKWIiIgQ0TjKsItGc3IoghAsYiXCHH/VU9Ag\nFCIaRxfhNppeCD/ESoQ5/qqnQETjaES8C6G9kW9YmOOvegqCFw0ztXd7jwYXOoYFCxZYpdaC0NaI\naIQ5/nIaELxomFp/86hYIbyZPHlyqLsgdGJCEo9QSnVVSn2ulMrx/G/0PEulVJZSarFS6iel1Fql\n1EWh6OuRTlJSUqs9jfT0dEBEQxCEpglVEPtW4Eut9SDgS89rX8qBK7TWQ4EzgCeUUlJD6EN6enpA\nTyOYwX2mNDMzM7N9OikIQqchVOGpGcDJnuWXgW+AP9obaK032ZZ3K6X2AenAwY7pYniQnp4elKcR\nKEE6YMAAXnnlFc4666z266ggCJ2CUHka3bXW+Z7lPUDA6T+VUuOAaGCzn+1XK6WWK6WWFxQUtG1P\nj1CGDRsGNCSxW1tyC3D55Zd7PQhJEATBiXbzNJRSXwA9HDbdYX+htdZKKR1gPxnAq8AsrbVj0F1r\nPReYC5Cdne13X52J7777jt27dwO0OhEuCIIQLO0mGlrraf62KaX2KqUytNb5HlHY56ddEvARcIfW\nekk7dTUsSUlJISWloX6gtYlwQRCEYAlVeOoDYJZneRbwvm8DpVQ08C7witb6nQ7sW9jhz9MIdhoR\nQRCEYAmVaDwMnKqUygGmeV6jlMpWSv3d0+ZCYApwpVJqtecvKzTdPbIxgtCanIYgCEIwhKR6Smtd\nBEx1WL8cmO1Zfg14rYO7Fpb4E4TIyMigHsIkCIIQLDLZUCfACILWutH6YMZpCIIgBIuIRifACILv\niO7IyEhqa2sBmchOEIS2QUSjExDI03BaFgRBaCkiGp0AEQ1BEDoKEY1OgIiGIAgdhYhGJ0BEQxCE\njkJEoxMgoiEIQkchotEJCFQ95bQsCILQUkQ0OgHiaQiC0FGIaHQCRDQEQegoRDQ6AcGIhgzuEwSh\nLRDR6ASIpyEIQkchotEJENEQBKGjENHoBEj1lCAIHYWIRidAPA1BEDoKEY1OgD/RsD/JT0RDEIS2\nQESjE2AEweVyOa73XRYEQWgpIhqdACMIMTExjuuVUo0eBSsIgtASRDQ6ASYB7k80ZIyGIAhthYhG\nJ8A80jU2NtZrvb+wlSAIQksR0egEGNEQT0MQhPZGRKMTUFlZCYinIQhC+yOi0QkwoiGehiAI7Y2I\nRidg/PjxAMycOdNrvYiGIAhtjViTTsCoUaOoqalpJA4SnhIEoa0RT6OT4ORNiKchCEJbI6LRiRFP\nQxCEtkZEoxMjnoYgCG2NiEYnRkRDEIS2RkSjE2NEwz7brSAIQmsQa9KJEdEQBKGtEWvSibHPcisI\ngtAWiGh0YoyHIZ6GIAhthViTToyEpwRBaGvEmnRiJDwlCEJbI6LRiTETGNbV1YW4J4IgdBZENDox\nCQkJoe6CIAidjJCIhlKqq1Lqc6VUjud/SoC2SUqpPKXU/3ZkHzsDiYmJgISnBEFoO0LladwKfKm1\nHgR86Xntj/uA7zqkV50M42lIIlwQhLYiVNZkBvCyZ/llYKZTI6XUGKA78FkH9atTERcXB4inIQhC\n2xEq0eiutc73LO+hQRi8UEpFAI8DNze1M6XU1Uqp5Uqp5QUFBW3b006AiIYgCG1Fu81kp5T6Aujh\nsOkO+wuttVZKaYd2/wV8rLXOa8roaa3nAnMBsrOznfZ1VGJKbk1uQxAEobW0m2horaf526aU2quU\nytBa5yulMoB9Ds0mAicqpf4LSASilVKlWutA+Q/BxsSJE7n99tu57rrrQt0VQRA6CaGaM/sDYBbw\nsOf/+74NtNaXmWWl1JVAtghG84iIiOCBBx4IdTcEQehEhCqn8TBwqlIqB5jmeY1SKlsp9fcQ9UkQ\nBEFoAqV150oBZGdn6+XLl4e6G4IgCGGFUmqF1jq7qXZSwC8IgiAEjYiGIAiCEDQiGoIgCELQiGgI\ngiAIQSOiIQiCIASNiIYgCIIQNJ2u5FYpVQBsb8Uu0oDCNupOWyL9ah7Sr+Yh/WoenbFffbXW6U01\n6nSi0VqUUsuDqVXuaKRfzUP61TykX83jaO6XhKcEQRCEoBHREARBEIJGRKMxc0PdAT9Iv5qH9Kt5\nSL+ax1HbL8lpCIIgCEEjnoYgCIIQNCIaHpRSZyilNiqlcpVSHf7cDqXUNqXUj0qp1Uqp5Z51XZVS\nnyulcjz/UzzrlVLqSU9f1yqlRrdhP/6hlNqnlFpnW9fsfiilZnna5yilZrVTv+5RSu3ynLPVSqmz\nbNtu8/Rro1LqdNv6Nr3OSqk+SqmvlVLrlVI/KaWu96wP6TkL0K+QnjOlVKxSaplSao2nX//tWd9P\nKbXUc4w3lVLRnvUxnte5nu2ZTfW3jfv1klJqq+18ZXnWd9h337PPSKXUKqXUh57XoTtfWuuj/g+I\nBDYD/YFoYA1wfAf3YRuQ5rPuz8CtnuVbgUc8y2cB8wEFTACWtmE/pgCjgXUt7QfQFdji+Z/iWU5p\nh37dA9zs0PZ4zzWMAfp5rm1ke1xnIAMY7Vl2A5s8xw/pOQvQr5CeM8/nTvQsu4ClnvPwFnCxZ/1z\nwLWe5f8CnvMsXwy8Gai/7dCvl4DzHdp32Hffs98bgX8CH3peh+x8iafRwDggV2u9RWtdDcwDZoS4\nT9DQh5c9yy8DM23rX9ENLAG6qIbH5rYarfV3wP5W9uN04HOt9X6t9QHgc+CMduiXP2YA87TWVVrr\nrUAuDde4za+z1jpfa73Ss3wI2AD0IsTnLEC//NEh58zzuUs9L12ePw38AnjHs973fJnz+A4wVSml\nAvS3rfvljw777iulegNnA3/3vFaE8HyJaDTQC9hpe51H4B9Ye6CBz5RSK5RSV3vWddda53uW9wDd\nPcsd3d/m9qMj+3edJzzwDxMCClW/PKGAUTTcpR4x58ynXxDic+YJtawG9tFgVDcDB7XWtQ7HsI7v\n2V4MpHZEv7TW5nw94Dlf/6OUivHtl8/x2+M6PgH8P6De8zqVEJ4vEY0jhxO01qOBM4HfKaWm2Dfq\nBh8z5KVuR0o/PDwLDACygHzg8VB1RCmVCPwf8AetdYl9WyjPmUO/Qn7OtNZ1WussoDcNd7tDOroP\nTvj2Syk1DP5/e3cPIlcVBXD8f2RFJYhRSBFR0AQLyRIVTeFHEaw0BgvBShQ12EQsBMVIIJ0QokVS\nxEILwQ8EVxTUUk0hKgRETVZQE41NiixElEBgEfakuHczM2HZebsz895G/j8Y5s6b2bln7rw3Z+fe\nN/fyKiW+bZQup1fajCkidgJzmflDm/Uux6RRnAZu7rt9U93Wmsw8Xa/ngE8pB9OZxW6nej1XH952\nvCuNo5X4MvNMPdAXgLfpfd1uNa6IuJLywfxBZn5SN3feZkvFtVbarMbyD3AEuJfSvTO1RB0X66/3\nXwecbSmuh2o3X2bmPPAO7bfX/cCjEfEXpWvwQeAQXbbXagZC/m8XYIoyYHUrvcG+LS3Wvw64tq/8\nHaUf9HUGB1MP1PIjDA7CHR1zPLcwOOC8ojgo/5GdogwEXl/LN0wgro195RcpfbYAWxgc9PuTMqA7\n9ve5vvZ3gYOXbO+0zZaJq9M2AzYA62v5GuAbYCcww+DA7u5afp7Bgd2Plot3AnFt7GvPg8D+Lvb9\n+tzb6Q2Ed9ZeY/ugudwvlLMhfqf0r+5tue5N9Q39GfhlsX5KX+RXwAngy8Wdr+6oh2usx4F7xhjL\nh5Rui/8o/Z67VhMH8CxlsO0k8MyE4nqv1nsM+IzBD8S9Na7fgIcn9T4DD1C6no4BP9XLjq7bbJm4\nOm0zYCvwY61/FtjXdwwcra99Briqbr+63j5Z7980LN4xx/V1ba9Z4H16Z1i1tu/3Pe92ekmjs/by\nF+GSpMYc05AkNWbSkCQ1ZtKQJDVm0pAkNWbSkCQ1ZtKQRhQR6yNidy3fGBEfD/sb6XLlKbfSiOrc\nTl9k5nTHoUgTNzX8IZKG2A9srpPdnQBuz8zpiHiaMvvoOuA24A3Kr6qfBOaBHZn5d0RspvxQbANw\nHnguM39t/2VIw9k9JY1uD/BHlsnuXr7kvmngMcqEd68B5zPzLuB74Kn6mLeAFzLzbuAl4M1WopZW\nwW8a0mQdybKexbmI+Bf4vG4/Dmyts9DeB8yUZQ+AMj+QtCaZNKTJmu8rL/TdXqAcf1dQ1ka4s+3A\npNWwe0oa3TnKkqorlmWNi1MR8ThcXHv6jnEGJ42TSUMaUWaeBb6NiFnKlOgr9QSwKyIWZzleC0sN\nS0vylFtJUmN+05AkNWbSkCQ1ZtKQJDVm0pAkNWbSkCQ1ZtKQJDVm0pAkNWbSkCQ1dgHulttdR8PG\nkQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Get original data\n", "infile = open(\"data/shot.txt\", \"r\")\n", "\n", "shot = []\n", "for line in infile:\n", " shot.append(float(line))\n", "time = range(0, 4*len(shot), 4)\n", "\n", "# Interpolate data with cubic polynomials\n", "poly_coeffs = np.polyfit(time, shot, 3)\n", "poly = np.poly1d(poly_coeffs)\n", "\n", "# sample the new interpolated polynomial with an interval of 7.07\n", "time2 = [x * 7.07 for x in range(int(time[-1]/7.07))]\n", "shot2 = poly(time2)\n", "\n", "# Plot both series\n", "pylab.plot(time, shot, 'k', label='raw data')\n", "pylab.plot(time2, shot2, 'r', label='re-sampled data')\n", "\n", "pylab.xlabel('time')\n", "pylab.ylabel('shot')\n", "pylab.legend(loc='best')\n", "\n", "pylab.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Question 2\n", "\n", "1. Calculate the upper triangular form for the matrix:\n", "\\begin{align*}\n", "A = \\begin{bmatrix}\n", "−5 & 3 & 4\\\\\n", "10 & −8 & −9\\\\\n", "15 & 1 & 2\n", "\\end{bmatrix}\n", "\\end{align*}\n", "\n", "2. Consider the matrix:\n", "\\begin{align*}\n", "C = \\begin{bmatrix}\n", "−5 & 3 & 4\\\\\n", "10 & −8 & −9\\\\\n", "15 & 1 & 2\n", "\\end{bmatrix}\n", "\\end{align*}\n", "Does matrix $C$ have an inverse? If not, then why not? " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "array([[-5, 3, 4],\n", " [ 0, -2, -1],\n", " [ 0, 0, 9]])\n" ] } ], "source": [ "from pprint import pprint\n", "import numpy.linalg as la\n", "\n", "# Compute the upper triangular matrix by Gaussian elimination\n", "def upper_triangle(A):\n", " # Assuming it is a square matrix of size nxn.\n", " n,m = A.shape\n", " if n != m:\n", " raise(\"Error, matric should be square\")\n", " \n", " # Create a local copy to avoid modifying input matrix\n", " A = A.copy()\n", " \n", " # Loop over each pivot row.\n", " for k in range(n-1):\n", " # Loop over each equation below the pivot row.\n", " for i in range(k+1, n):\n", " # Define the scaling factor outside the innermost\n", " # loop otherwise its value gets changed as you are\n", " # over-writing A\n", " s = (A[i,k]/A[k,k])\n", " # we don't start the following loop from 0 as we can assume\n", " # some entries in the rows are already zero\n", " for j in range(k, n): \n", " A[i,j] = A[i,j] - s*A[k,j]\n", " \n", " return A\n", "\n", "A = np.array([[-5,3,4], [10,-8,-9], [15,1,2]])\n", "U = upper_triangle(A)\n", "pprint(U)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The determinant of C is non null, so C has a unique inverse.\n" ] } ], "source": [ "C = np.array([[-5,3,4], [10,-8,-9], [15,1,2]])\n", "if not np.isclose(np.linalg.det(C), 0):\n", " print(\"The determinant of C is non null, so C has a unique inverse.\")\n", "else :\n", " print(\"The determinant of C is null, so C cannot be inverted.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Question 3\n", "\n", "Consider the function\n", "$$f(x) = \\dfrac{1}{(x − 0.3)^2 + 0.01} - \\dfrac{1}{(x − 0.8)^2 + 0.04}$$.\n", "\n", "\n", " 1. Write a function that computes the second derivative of $f(x)$ using central differencing. The interface to your function should look like *central_diff2(f, x, h)* where $f$ is the function to be differentiated, $x$ is the position at which the derivative should be estimated, and $h (= \\Delta x)$ is the step size.\n", " 2. Use this function to compute the derivative at $x = 0.5$ for decreasing values of $\\Delta x$. Start with $\\Delta x=1.2$ and keep halving until the relative difference between solutions falls below $1.0^{-6}$ Plot the convergence of the method, i.e. plot $\\Delta x$ against the absolute difference between the analytical value and the finite difference approximation." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def central_diff2(f, x, h):\n", " return (f(x+h)-f(x-h))/(2*h)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have $$f'(x) = \\frac{-2(x-0.3)}{((x − 0.3)^2 + 0.01)^2} + \\frac{2(x-0.8)}{((x − 0.8)^2 + 0.04)^2}\\,.$$\n", "\n", "There are two errors involved in this exercise: \n", " - the relative error between two successive approximations is used as a criterion to stop the loop: $\\left|\\frac{f'_{approx, dx}-f'_{approx, dx/2}}{f'_{approx, dx}}\\right|$\n", " - the error with respect to the exact value is stored to plot the convergence: $\\left|f'_{approx, dx}-f'_{exact}\\right|$\n", " \n", "Note that as the discretisation interval tends to 0, the approximate value of the derivative tends to the exact value, and thus both errors become very similar. This means that in a real problem where you cannot compute the exact value, using the difference between the last two values is a good indicator for the convergence." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/usr/local/lib/python3.6/site-packages/matplotlib/axes/_axes.py:545: UserWarning: No labelled objects found. Use label='...' kwarg on individual plots.\n", " warnings.warn(\"No labelled objects found. \"\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEQCAYAAACnaJNPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VPW5x/HPk7ApqyziAgrIIrgggqCgNlStqCBVUQH1\nVquCVFDb20VrF29bq10sRVGU60JFISpWBMRapEZcqGLYZJGlESWIoCCBsCQkee4fGbkhkslMkpMz\nM/m+X6+8Muc3Z855fDrNw/n9zvn9zN0RERGJVVrYAYiISHJR4RARkbiocIiISFxUOEREJC4qHCIi\nEhcVDhERiYsKh4iIxEWFQ0RE4qLCISIicVHhEBGRuNQLO4AgtG7d2jt06BDoOXbv3k3jxo0DPUcy\nU36iU34qpxxFF0R+srOzv3T3NpXtl1KFw8yGAEM6d+7MBx98EOi5srKyyMjICPQcyUz5iU75qZxy\nFF0Q+TGzT2LZL6W6qtx9truPat68edihiIikrJQqHCIiEryUKhxmNsTMJufl5YUdiohIykqpwqGu\nKhGR4KVU4RARqcvGz1tbK+dR4RARCUEQf+QnzF9X48c8lJS6HVdEJFlMmL+OH17QNeo+JSVOYXEJ\n+4tLKCwqYX+xU1hUQmFxCZ/uLGbZxh2l7xeVUFBcAsBXuws5onGDQGNPqcJR9jkOEZFE4+68//F2\npry7AYDBD731jYJQul36u6jEox/w3Xe+0dTrt/MAuP28LpUWpqpKqcLh7rOB2X369Lk57FhERL7m\n7sxfvZVfvryCzXn7DrSv2LQTgK5tm3D6ca1oUM9okJ5G/fQ0GtT7/98Nym3XTzfWfbSaXj1POaj9\n8kfeZeX/XEjjhsH+aU+pwiEikkiKikt45cPNPPLGf1izZRfHtjiM3ww9iav6tOfEX/6DDfdfUuVj\nZ+1YR0b3tt9oD7pogAqHiEiN27e/mBnZuUxekMOn2/fQ5cgm/OWqngzpeQz104O7J+n287oEduyy\nVDhERGrIrn37efa9T3ni7Y/5YlcBp7VvwS8u6c753duSlmYH7RvEH/mgxjTKS/jCYWbfBS4BmgFP\nuPs/Qw5JROQg2/ILeOqdDTy9cAM79xVxTpfWTBh+Gmd1aoWZHfIztfVHPgihFA4zexIYDGx195PL\ntA8CJgDpwOPufr+7zwRmmtkRwJ8BFQ4RSQibduzlfxfkkLnoUwqKSriwx1H8YOAJnNquRdihBSqs\nK44pwETg6a8bzCwdeBi4AMgFFpnZLHdfFdnlF5H3RURCtX7rLiZl5fDy0k0AfLfXsdzyrRPofGST\nkCOrHaEUDndfYGYdyjX3Bda7ew6AmWUCQ81sNXA/8Kq7L67VQEUkaYyftzaQ7p+yx122cQePZK3n\nn6u20LBeGteeeTw3n9uJY1scVuPnTWSJNMZxLLCxzHYu0A8YB5wPNDezzu7+6KE+bGajgFEAbdu2\nJSsrK9Bg8/PzAz9HMlN+olN+KhdPjkrcmTB/Dx2Kc6mfZtRPg/rpkFbB+EI8JszfTcOdG3klp5CV\n20o4vB4M7lSfC46vT7MGX7Bu6RfUzkQfBwvzO5RIheOQ3P1B4MEY9ptsZpuBIU2bNu0d9MphWp0s\nOuUnOuWncrHm6PO8fdzx3BJgDz/M2nvQe/XSjAb10mhYL42G9dJpWL/M63ppNKxf+nBdRe+V3jq7\njj8u2kebpg2566KOjOx3HE0b1Q/kvzkeYX6HEqlwbALal9luF2kTETmk26YvZtayzd9oP7tzK/p2\nbEVBUTEF+0soKCqdwqOgqJiCopLIT+l7O/cWUVBUHHm/9GfX3v3sLzfdxxe7CthTWJwQRSNsiVQ4\nFgFdzKwjpQVjODAyngNoyhGRuqGwqIQ//OMjZi3bTPejmzFxZC/Oe+DNaj2JfShFxSV0vvvVGj9u\nsgvrdtzpQAbQ2sxygV+7+xNmNhZ4jdLbcZ9095VxHleTHIqkuA1f7mbc9CV8uCmP/zrreH5+cXca\n1U8P5Fz1AnzKO5mFdVfViAra5wJzq3FcXXGIpLCXl27i7pdWkJ5mPHptbwadfNSB94KabqO2pvFI\nJonUVVVtuuIQSU17Cou4Z9ZKnv8glz7HH8GEEb2+cQtsUE9iJ/MT3kFJqeswrTkuknpWb97JkIfe\n5oXsXMYO7EzmqDPr3HMTiUZXHCKSkNydZ/79Cb+Zs4rmh9XnmRv7MaBz67DDElKscGiMQyQ15O3d\nz8NLC/hgywrO7dqGB67sSZumDcMOSyJSqnCISPJb/OlXjJu2hM/zirnrohO5+ZxO35iSXMKVUmMc\nZjbEzCbn5eWFHYqIxKmkxJmU9R+ufHQhZvDzfo0Y/a0TVDQSUEoVDg2OiySnL3YV8L2n3ucP//iI\nQScdxSu3ncMJLYJ5NkOqT11VIhKqt9Z9wQ+fW8auffv5/WWnMKJv+woXP5LEoMIhIqHYX1zCX+at\n5dE3/0PnNk149qZ+dDuqadhhSQxSqnDodlyR5PA/s1ayLHcHiz/dwYi+7fnV4JM4rIG6ppKFxjhE\npFb9Y8Vmnnp3A+u25PPQiF7cd/mpKhpJJqWuOEQkcRUVl3D/qx/x+NsfA/DKbedwXKvDQ45KqkKF\nQ0QCt313IZdOfJvcr/5/oaVz//QGUDqJoOaDSi4qHCISqBWb8hg9NZsv8gv407BTubJPezrc+YrW\nuEhiKVU4NDguklhmLtnEz15cTsvGDXhh9Fn0bN8i7JCkBmhwXERqXFFxCb+ds4o7nltKz/YtmDX2\n7IOKhta4SG4pdcUhIuHbll/A2GlLWJizjev7d+DuS7pTv9xKehrTSG4qHCJSY8qOZ/z5yp4M690u\n7JAkAAlfOMysE3A30Nzdh4Udj4gc2ktLcrnzxQ9p1bgBM245i1PbaTwjVYUyxmFmT5rZVjNbUa59\nkJmtMbP1ZnYngLvnuPuNYcQpIpXbX1zCb2av4ofPLeO09i2YNe5sFY0UF9bg+BRgUNkGM0sHHgYu\nAnoAI8ysR+2HJiKx2pZfwHVPvMeT73zMDQM68MxN/WjdRAsupbpQuqrcfYGZdSjX3BdY7+45AGaW\nCQwFVtVudCISiw9z8xg99QO27S7kL1f15PLTNZ5RV5i7h3Pi0sIxx91PjmwPAwa5+02R7euAfsCv\ngXuBC4DH3f2+Co43ChgF0LZt296ZmZmBxp+fn0+TJk0CPUcyU36iS/b8vLNpP0+tLKRZA+O2Xg3p\n0Lzm55pK9hwFLYj8DBw4MNvd+1S2X8IPjrv7NuCWGPabbGabgSFNmzbtnZGREWhcWVlZBH2OZKb8\nRJes+dlfXMK9r6xmyocbOLNTSx4eeTqtAuqaStYc1ZYw85NIDwBuAtqX2W4XaYuZHgAUCc6X+QVc\n+/h7THl3A98f0JFnbuwXWNGQxJZIVxyLgC5m1pHSgjEcGBnPATTliEgwlufuYPTUbLbvLmT81T25\nrJfGM+qysG7HnQ4sBLqZWa6Z3ejuRcBY4DVgNfC8u6+M57i64hCpWePnrWVGdi7DHl1Imhkvjumv\noiGh3VU1ooL2ucDcqh5XVxwiNWd/cQkT5q8DoP8JrXhoRC91TQmQWF1VIpIgvthVwK3TFgNw49kd\nueuiE6mXnkhDohKmlPomqKtKpPp+OmMZZ9z7Ou9/vB2AJ97+mM53v8r4eWtDjkwSha44ROSAae99\nyktLNtG+5WE8em1vLnnwbS24JN+QUoVDYxwiVbNvfzH3zFpJ5qKNnNu1DQ8OP40WhzcIOyxJUOqq\nEqnjPtuxl6sfW0jmoo2MHdiZp64/40DR0IJLcii64hCpwxb+Zxtjpy2moKiER6/tzaCTjzrofS24\nJIeiKw6ROsjdefytHK594j1aHF6fmbcO+EbREKlISl1xiEjl9hQW8bMXP2T2ss+48KS2/PnKnjRt\nVD/ssCSJqHCI1CGfbNvN6KnZrNmyi59c2I0fZJyAmYUdliSZlCocGuMQqdgbH23l9swlmBlTbujL\nt7q2CTskSVIa4xBJcSUlzoPz1/H9vy2i3RGHM2fc2SoaUi0pdcUhIgfbuW8/P3puGa+v3sJlvY7l\n95edwmENan7RJalbVDhEUtTaLbsYPTWbjdv3cM+QHnyvfweNZ0iNUOEQSUFzP9zMj19YxuEN6vHs\nTf3o16lV2CFJCkmpwqHBcanriopL+NM/1/DYmzn0Oq4Fk67pzVHNG4UdlqQYDY6LpIDx89ayfXch\n1z+1iMfezOGafseROepMFQ0JREpdcYjUVRPmr2NGdi5f7Crgj1ecylVntA87JElhKhwiSe7F7FwA\nStx54Zaz6Nm+RcgRSapT4RBJUn9+bQ0T31h/YHtz3j6GPvwOt5/XRZMTSqASvnCYWWPgEaAQyHL3\nZ0MOSSR0W3bu493/fAnAqHM7MXlBjhZckloTyuC4mT1pZlvNbEW59kFmtsbM1pvZnZHmy4EZ7n4z\ncGmtByuSYN7/eDuXPPg2H32+i4kje/Hzi7uHHZLUMWHdVTUFGFS2wczSgYeBi4AewAgz6wG0AzZG\ndiuuxRhFEoq7M+Wdjxn5v/+maaN6zLx1AINPPQbQgktSu+LqqjKzlu6+vbondfcFZtahXHNfYL27\n50TOlQkMBXIpLR5LSbHbh0VitbewmLv+vpyZSz/j/O5t+cvVPWlWZip0jWlIbTJ3P/QbZr9w999F\nXvcAZgL1AQOudvf3qnXi0sIxx91PjmwPAwa5+02R7euAfsDPgInAPuDtisY4zGwUMAqgbdu2vTMz\nM6sTXqXy8/Np0qRJoOdIZspPdPHkZ+ueEh5aUkDurhIu61KfwZ3qk1YHpg7Rdyi6IPIzcODAbHfv\nU9l+0a44Lgd+F3n9J+B2d3/VzPoCfwX6Vz/Myrn7buCGGPabbGabgSFNmzbtnZGREWhcWVlZBH2O\nZKb8RBdrft5Ys5V7M5cC6Tx5Q28Gdjsy8NgShb5D0YWZn1i7qo5x91cB3P19MzssgFg2AWWfWmoX\naROpc0pKnIffWM9fXl/LiUc147Fre3Ncq8PDDksEiD5m0MnMZpnZbKCdmZX91gaxzuQioIuZdTSz\nBsBwYFY8B9CUI5IKdu7bz6ip2Twwby1Dex7D38f0V9GQhBLtimNoue00ADNrC0yqzknNbDqQAbQ2\ns1zg1+7+hJmNBV4D0oEn3X1lnMfVJIeS1DQVuiSDCguHu79ZQfsWSm+brTJ3H1FB+1xgbjWOOxuY\n3adPn5uregyRsMxZ/hk/nbGcxg3rMX3UmZzRoWXYIYkcUpVub43cwZRwzGyImU3Oy8sLOxSRmBUV\nl/D7uasZO20J3Y9uxpxxZ6toSEKr6nMRCXntrDEOSTbb8gu47on3mbwgh/8663im33wmbZtpKnRJ\nbFWaq8rdH6vpQGqCxjgkmSzbuIMxz2SzbXchf76yJ8N6tws7JJGYRL3iMLMLzWxS5O6qWZHXg6J9\nJky64pBEN37eWgDezN3PlY8uJC3NeHFMfxUNSSoVXnGY2V+BrsDTlE77AaXPVtxmZhe5++21EJ9I\nSpkwfx1bdxUwfUUh53RpzYPDe3FE4wZhhyUSl2hdVRe7+zcmwDGz54C1QMIVDnVVSSLbnLcXgOnv\nf8rgTvWZcENf0tMScrhQJKpoXVX7zOyMQ7SfQem8UQlHXVWSiMbPW0uHO1/hrPv+daBtTs5+Hpy/\nLsSoRKou2hXH9cAkM2vK/3dVtQfyIu+JSCXcnSYN65GeZnRs3Zj1W/PZcP8lkXmGNKOtJKdoDwAu\nBvqZ2VHAsZHmTe7+ea1EVgXqqpJEsrugiJ++uJxXlm/m4lOO4o/DenLyr18LOyyRaqv0dlx3/9zM\n9rj7TjNrVhtBVZWeHJdEkfNFPrc8k836rfncddGJjDq3E2amBZckJcT6HEcWcHqZ3yJSgXmrtvCj\n55ZSv14aU2/sx4DOrQ+8pwWXJBXE+wCgbgERqUBxifPX19fy0L/Wc2q75ky6tjfHtghiBQKRcFXp\nyXEROdiOPYXcnrmUN9d+wdV92vM/Q0+iUf30sMMSCURKFQ4NjksYVn6Wxy3PZLMlr4D7Lj+FEX2P\nCzskkUDFO8nhoRcoTxB6jkNq298X53L5I++yv8h5bvSZKhpSJ8R6xWHlfovUaYVFJdz7yir+tvAT\nzuzUkokjT6d1k4ZhhyVSK2ItHFeX+y1SZ23ZuY8fPLuY7E++4uZzOvKzQSdSL72qKxSIJJ9KC4eZ\nHenuawG+/m1m3dx9TdDBiSSaRRu284NnF7O7oIiHRvRiSM9jwg5JpNbF8s+kt8zsqq83zOy/gZeC\nC+lgZtbJzJ4wsxm1dU6R8tydKe98zIjJ/6ZJw3q89IMBKhpSZ8VSODKA68zsBTNbQOlU631jObiZ\nPWlmW81sRbn2QWa2xszWm9md0Y7h7jnufmMs5xMJwt7CYn70/DLumb2KjG5H8vLYAXQ7qmnYYYmE\nJpYpRzab2T+Au4AS4E53z4/x+FOAiZSu6QGAmaUDDwMXUDp54iIzmwWkA/eV+/z33X1rjOcSqVHj\n563litPbMfqZbD76fCf/fUFXbh3YmTRNhS51XCxjHK8DnwEnUzo77hNmtsDdf1zZZ919gZl1KNfc\nF1jv7jmR42cCQ939PmBwfOGLBGfC/HVMeXcDAE9efwYDux0ZbkAiCSKWu6omuvvMyOsdZtaf0quP\nqjoW2FhmOxfoV9HOZtYKuBfoZWZ3RQrMofYbBYwCaNu2LVlZWdUIsXL5+fmBnyOZJXN+StyZk7Mf\ngGb1ihnXqyG2eRVZm1fV2DmSOT+1RTmKLsz8xNJVNbPcdhHw28Ai+ub5twG3xLDfZDPbDAxp2rRp\n74yMjEDjKl1PIdhzJLNkzc/9r67m0TdzDmxv3FXCTxfs5fbzutToBIXJmp/apBxFF2Z+wphyZBOl\nXV5faxdpqzZNqy7VsebzXby2cgv10oxfXNKde2avYsP9l4QdlkjCCeOppUVAFzPraGYNgOHArJo4\nsJkNMbPJeXl5NXE4qUNmL/uM7z78DvkFRUwfdSbXD+gYdkgiCavSwmFmt8fSVsFnpwMLgW5mlmtm\nN0a6usYCrwGrgefdfWV8YYvUjKLiEn43ZxXjpi/hpGOa8cq4szmjQ0sALbokUoFYuqq+B0wo13b9\nIdq+wd1HVNA+F5gbw7njoq4qiceX+QWMnbaYf+ds5/r+Hfj5xd1pUO///y2lRZdEDq3CwmFmI4CR\nQMfIcxZfawpsDzqwqtC06hKrJZ9+xZhnFrNjbyHjr+7JZb3ahR2SSNKIdsXxLrAZaA08UKZ9F7A8\nyKCqSlccUhl3Z/r7G7ln1kraNm/Ii2P6c9IxmoZfJB4VFg53/wT4BDir9sIRCc6+/cX86uUVPP9B\nLt/q2oYJw0+jxeENwg5LJOnE8uT4mcBDQHegAaVTg+x292YBxxY3dVVJRTbt2MuYZ7JZnpvHbd/u\nzO3ndyVdU4eIVEkst+NOBEYA64DDgJsonWsq4WgFQDmUd9Z/yZCH3ubjL3bzv//Vhx99p5uKhkg1\nxPQch7uvB9LdvdjdnwIGBRtW1eg5DinL3ZmU9R+ue+I9WjdpwMtjB3BBj7ZhhyWS9GK5HXdP5EG9\npWb2R0oHzBNyuTMNjsvX8guK+MkLy3h1xedccurR/PGKU2ncMIyJEkRSTyz/T7qO0kIxFvghpdOF\nXBFkUCLVsX5rPqOnfsCGbXu4++Lu3HROR8zUNSVSU2KZ5PCTyMt9wP8EG45I9fxjxef8+IVlNKyX\nxtQb+9L/hNZhhySSclLq2l13VdVdxSXOn/+5hklZ/6Fn+xZMuuZ0jmlxWNhhiaSklCocGuOom34/\ndzWrN+/krXVfMrLfcfx6SA8a1ksPOyyRlBW1cESWef1DLKv9iYRhxaY8Ji/IoUG9NP5wxSlcfcZx\nYYckkvKiFg53Lzazs2srGJF4vPDBRu6euQKAGbecxantWoQckUjdEEtX1ZLIJIcvALu/bnT3vwcW\nlUgUBUXFDJv0Lh9u2nmg7dKJ7wDU+Ep9IvJNsRSORsA24Ntl2hxIuMKhwfHU93nePsY8m82Hm3Yy\n+lud+Ml3utH57le1Up9ILYrldtwbaiOQmqDB8dT275xtjJ22mL2FxTxyzelcfMrRYYckUifFsgJg\nOzN7ycy2Rn5eNDMtXiC1xt15/K0crnn8PZodVp+Xxw44qGhopT6R2hXL1CFPUbom+DGRn9mRNpHA\n7Sks4rbMpfzuldWc3/1IXr51AJ2PbHrQPhrTEKldsYxxtIlMbPi1KWZ2R1ABHYqZfRe4BGgGPOHu\n/6zN80s4Nny5m9FTs1m3dRc/HdSNMd86QVOHiCSAWK44tpnZtWaWHvm5ltLB8piY2ZORLq4V5doH\nmdkaM1tvZndGO4a7z3T3m4FbgKtjPbckr/mrtzBk4tts2bWPv32/Lz/I6KyiIZIgYrni+D6lCzmN\np/RuqneBeAbMp1C6psfTXzdEHix8GLgAyAUWRW75TQfuK39+d98aef0LEnQtEKkZJSXOX+ev48H5\n6zj52GZMuqY37VseHnZYIlJGLE+OX+7ul1b1BO6+wMw6lGvuC6x395zIeTKBoe5+HzD4EHEYcD/w\nqrsvrmosktjy9uznjueW8MaaLxjWux2/++7JNKqvqUNEEo25e/QdzN53977VOklp4Zjj7idHtocB\ng9z9psj2dUA/dx9bwedvA74HLAKWuvujh9hnFDAKoG3btr0zMzOrE3Kl8vPzadKkSaDnSGbx5ufT\nncU8tKSA7fuca7o3YGD7eindNaXvT+WUo+iCyM/AgQOz3b1PZfvF0lX1jplNBJ7j4CfHa+1f/u7+\nIPBgJftMBiYD9OnTxzMyMgKNKSsri6DPkcziyc/MJZv4/fzlND+sAS/c0JvTjzsi2OASgL4/lVOO\nogszP7EUjtMiv39Tps05+EnyeG2idEGor7WLtFWLnhxPLvuLS7j3ldVMeXcDfTu25OGRp9OmacOw\nwxKRSlQ2xpEGTHL352v4vIuALmbWkdKCMRwYWcPnkAS2dec+bp22mEUbvuLGszty50UnUj89IVck\nFpFyov4/1d1LgJ9W5wRmNh1YCHQzs1wzu9HdiyhdivY1YDXwvLuvrM55IvHOdvdRzZs3r+6hJEDZ\nn2xn8ENvs2LTTiYMP41fDu6hoiGSRGLpqnrdzH7MN8c4tsdyAncfUUH7XGBuLMeIlbqqEtv4eWto\n1aQhv5m9imOPOIynb+zLiUc1CzssEYlTLIXj6wfubi3T5kCnmg+nejTJYeLat7+YCfPXA3DeiUfy\nl6tPo/lh9UOOSkSqIpbZcTvWRiA1QVcciWnj9j2MnpoNwA/P78q4b3cmLS11b7UVSXUVdiyb2U/L\nvL6y3Hu/DzKoqtIYR+K5I3MJ5/zxDVZtLl10afzra+n087mMn7c25MhEpKqijUgOL/P6rnLvDQog\nlmozsyFmNjkvLy/sUOq8Encm/msdLy/7jBOPasqbP8kAYMP9l7Dh/ks0o61IEotWOKyC14faTgi6\n4kgMO/ft56ElBfz5n2sZ2vMYXvrBAI5v1TjssESkhkQb4/AKXh9qWwSAtVt2MXpqNp9uK+aeIT34\nXv8OB6YO0YJLIqkhWuHoaWY7Kb26OCzymsh2o8AjqwINjodrzvLP+OmM5TRuWI+f9W3E9QMOvq9C\n3VMiqaHCrip3T3f3Zu7e1N3rRV5/vZ2Q91GqqyocRcUl/H7uasZOW0L3o5sxZ9zZdD1Cs9qKpKpY\nnuMQqdCX+QWMm7aEhTnb+K+zjucXl/SgQb00VocdmIgEJqUKh7qqatfSjTsY80w223cX8sCVPbmi\nd7uwQxKRWpBSEwSpq6r2TH//U656dCHpacaLY/qraIjUISl1xSHB27e/mHtmrSRz0UbO7dqGB4ef\nRovDG4QdlojUIhUOidmmHXsZ80w2y3PzGPftztxxflfSNXWISJ2jwiExeWf9l4ybvoT9RSVMvq43\n3znpqLBDEpGQpFTh0OB4zXN3Ji/I4Q//+IgT2jThset606mN1oEWqcs0OC4Vyi8o4tZpi7nv1Y+4\n6OSjmXnrABUNEUmtKw6pGePnreXS045h9NRsPv5yN3df3J2bzul4YOoQEanbVDjkGybMX8cTb39M\nw3ppTL2xL/1PaB12SCKSQFQ45IDiEucv89YAcMKRTZh0zekc0+KwkKMSkUST8IXDzLoDtwOtgfnu\nPinkkFLSfXNX89iCnAPbyzbuoP/9/+L287pockIROUigg+Nm9qSZbTWzFeXaB5nZGjNbb2Z3RjuG\nu69291uAq4ABQcZbV63YlMec5ZtpkJ7G/ZefAmjBJRGpWNB3VU2h3GqBZpYOPAxcBPQARphZDzM7\nxczmlPs5MvKZS4FXgLkBx1vnzMjO5YpJ7+LuvHDLWQzve1zYIYlIggu0q8rdF5hZh3LNfYH17p4D\nYGaZwFB3vw8YXMFxZgGzzOwVYFpwEdcdhUUl/GbOSp7596f0P6EVD43oRasmDQEtuCQi0Zl7sIv5\nRQrHHHc/ObI9DBjk7jdFtq8D+rn72Ao+nwFcDjQElrv7wxXsNwoYBdC2bdvemZmZNfsfUk5+fj5N\nmiTnMw1f7Svh4aUFrN9RwkUd6zOsS/0anzokmfNTG5SfyilH0QWRn4EDB2a7e5/K9kv4wXF3zwKy\nYthvspltBoY0bdq0d0ZGRqBxZWVlEfQ5gvBezjZ+PG0xewuNR645nYtPOTqQ8yRrfmqL8lM55Si6\nMPMTxpPjm4D2ZbbbRdqqTU+OV8zdeeLtjxn5+Hs0a1SfmbcOCKxoiEhqC+OKYxHQxcw6UlowhgMj\na+LAmqvq0PYUFnHnix8ya9lnfKdHWx64qidNGyXk6r8ikgQCLRxmNh3IAFqbWS7wa3d/wszGAq8B\n6cCT7r4yyDjqsg1f7uaWZ7JZs2UXP7mwG2O+dQJpmgpdRKoh6LuqRlTQPpcAbq1199nA7D59+txc\n08dORvNXb+GO55aSnmb87Ya+nNu1TdghiUgKSPjB8Xioq6pUSYkzYf46Jsxfx0nHNOPRa3vTvuXh\nYYclIilC06qnmLw9+7nxb4uYMH8dV5zejhfH9FfREJEapSuOFLJ6805GT81mc95efvvdk7m233Ga\nCl1EapyAypNlAAALR0lEQVSuOFLEzCWbuOyRd9i3v5jMUWdy3ZnHq2iISCBS6oqjLnrgn2vYta+I\nKe9uoG+Hlky8phdHNm0UdlgiksJSqnDUta6qrbv28dC/1gPw/QEdueviE6mfnlIXkSKSgFLqr0xd\n6qrK/mQ7gx98G4AJw0/jV0N6qGiISK3QX5ok4+7c8NT7XDFpIVt3FQBwe+ZSOtz5CuPnrQ05OhGp\nC1KqqyrV7dtfzM9f+pA31nzBwG5t+OvVvej5m3+y4f5Lwg5NROqQlCocqTzGsXH7HkZPzWbV5p3c\ncX4Xbvt2F00dIiKhSKmuqlQd43hz7RcMfuhtNn61hye+14c7zu96oGho0SURqW0pdcWRakpKnEey\n1vPAvLV0a9uUR6/tTYfWjQ/aR2uCi0htU+FIUDv37ee/n1/GvFVbGHraMdx3+Skc3kD/c4lI+PSX\nKAGt27KL0VOz+WT7Hn41uAc3DOigp8BFJGGkVOFIhcHxV5Zv5iczlnF4g3pMu6kf/Tq1CjskEZGD\naHA8QRQVl/D7uau5ddpiuh3VlDnjzlbREJGElFJXHMlqW34BY6ctYWHONq4783h+ObgHDeqlVE0X\nkRSiwhGypRt3MOaZbLbvLuTPV/ZkWO92YYckIhJVUvyz1swam9kHZjY47Fhq0vT3P+WqRxeSnma8\nOKa/ioaIJIVAC4eZPWlmW81sRbn2QWa2xszWm9mdMRzqZ8DzwURZ+/btL+bOF5dz198/pF+nlswe\nezYnH5t84zIiUjcF3VU1BZgIPP11g5mlAw8DFwC5wCIzmwWkA/eV+/z3gZ7AKiAlFpn4bMdexjyT\nzbLcPG4deAI/uqAb6Zo6RESSSKCFw90XmFmHcs19gfXungNgZpnAUHe/D/hGV5SZZQCNgR7AXjOb\n6+4lQcYdhPHz1tKvY0vGTl9CYVEJj13XmwtPOirssERE4mbuHuwJSgvHHHc/ObI9DBjk7jdFtq8D\n+rn72EqOcz3wpbvPqeD9UcAogLZt2/bOzMysqf+EQ8rPz6dJkyYx7evu3PDaHgw4uokx7rRGHN0k\nKYaXqiye/NRFyk/llKPogsjPwIEDs929T2X7Jc1dVe4+pZL3JwOTAfr06eMZGRmBxpOVlUUs58gv\nKOJnM5YDe7jolKP447CeNGmYNGmvsljzU1cpP5VTjqILMz9h/AXbBLQvs90u0lZtifbk+K9eXsHT\nCz85sD33w8+Z++Hn3H5eF01OKCJJK4zCsQjoYmYdKS0Yw4GRIcQRqNdWfs7fF2+iZeMGTBzRi5GP\nv6cFl0QkJQR9O+50YCHQzcxyzexGdy8CxgKvAauB5919ZU2cLxGmHCkucf702keMnprNCW0aM3vc\n2fTv3Dq0eEREalrQd1WNqKB9LjC3ps8XdlfVV7sLuS1zCW+t+5LhZ7TnnktPolH9dEALLolI6kip\nW3vCvOJYsSmPIRPf5r2c7dx3+Sncf8WpB4oGaMElEUkdKXV7T1hXHDOyc7n7pQ9p2bgBz99yFqe1\nb1Gr5xcRqU264qiGwqISfjlzBT9+YRmnH3cEs8edraIhIilPVxxV9NW+EoZPXsjiT3cw+txO/OTC\nbtRLT6k6LCJySCn1l662rjjey9nGr9/dx0ef7+Lhkadz18XdVTREpM5IqSuOoLk7T72zgXvnrqZN\nI3h69AC6tm0adlgiIrUqpQpHkF1VewqLuOvvH/Ly0s+4oEdbLjt6l4qGiNRJKdW/ElRX1SfbdnP5\nI+8ya9ln/OTCbjx2bW8Or6+p0EWkbkqpK44gvPHRVm7PXEJamjHlhr58q2ubsEMSEQlVShWOmuyq\nKilxHvzXOibMX0f3o5rx2HW9ad/y8OoHKSKS5NRVVc74eWvJ27Ofm57+gL++vo7Leh3L33/QX0VD\nRCQipa44asKE+euYuXQTn+3Yy2+HnsS1Zx6PmcYzRES+psJRxstLS5cF2VtYTOaoM+l9fMuQIxIR\nSTwqHJR2T02Yv+7A9tZdBVwxaaEWXBIROYSUKhxVHRz/4QVd+eEFXSkucU74+VwtuCQiEoUGx8tI\nT9NYhohIZVKqcNQELbgkIhKdCkc5GtMQEYlOhUNEROKS8IXDzDLM7C0ze9TMMsKOR0Skrgu0cJjZ\nk2a21cxWlGsfZGZrzGy9md1ZyWEcyAcaAblBxSoiIrEJ+nbcKcBE4OmvG8wsHXgYuIDSQrDIzGYB\n6cB95T7/feAtd3/TzNoCfwGuCThmERGJItDC4e4LzKxDuea+wHp3zwEws0xgqLvfBwyOcrivgIZB\nxCkiIrEL4wHAY4GNZbZzgX4V7WxmlwMXAi0ovXqpaL9RwKjIZr6ZrYm8bg7kHeIjh2ov3xZtuzXw\nZUXxVENF8VZn/2j7pHp+Yv1MvDlSfip/r7J8lG8r/34QOUrk/Byqrbbzc3xMe7l7oD9AB2BFme1h\nwONltq8DJgZ4/smxtpdvi7YNfFCb8VZn/2j7pHp+gsqR8hP/d6iynB1i/xrPUSLnJ4bvTOD5ifUn\njLuqNgHty2y3i7QFZXYc7eXbKtsOQrzniGX/aPuken5i/Uy8OVJ+Kn8vlv/+2ZW8X9MSOT+Haqvt\n/MTEIpUruBOUjnHMcfeTI9v1gLXAeZQWjEXASHdfGWggNczMPnD3PmHHkaiUn+iUn8opR9GFmZ+g\nb8edDiwEuplZrpnd6O5FwFjgNWA18HyyFY2IyWEHkOCUn+iUn8opR9GFlp/ArzhERCS1JPyT4yIi\nklhUOEREJC4qHCIiEhcVjoCYWWMz+8DMoj0NXyeZWffIpJUzzGxM2PEkGjP7rpn9r5k9Z2bfCTue\nRGNmnczsCTObEXYsiSLy9+Zvke9N4NMyqXCUU0MTMwL8DHg+mCjDUxP5cffV7n4LcBUwIMh4a1sN\n5Wemu98M3AJcHWS8ta2G8pPj7jcGG2n44szV5cCMyPfm0sBj011VBzOzcymdjffpMs+epFP67MmB\niRmBEVQ8MWNPoBWlM/p+6e5zaif64NVEftx9q5ldCowBprr7tNqKP2g1lZ/I5x4AnnX3xbUUfuBq\nOD8z3H1YbcVe2+LM1VDgVXdfambT3H1kkLGFMVdVQvMamJgxsm5IY6AHsNfM5rp7SZBx15aayE/k\nOLOAWWb2CpAyhaOGvj8G3E/pH4KUKRpQ4xOfprR4ckVpEWkHLKUWepJUOGIT18SM7n43gJldT+kV\nR0oUjSjinbgyg9JL64bA3EAjSwxx5QcYB5wPNDezzu7+aJDBJYB4vz+tgHuBXmZ2V6TA1BUV5epB\nYKKZXUItTE2iwhEgd58SdgyJyN2zgKyQw0hY7v4gpX8I5BDcfRul4z8S4e67gRtq63waHI9NbU/M\nmGyUn+iUn+iUn9glRK5UOGKzCOhiZh3NrAEwHJgVckyJRPmJTvmJTvmJXULkSoWjnBSfmLHalJ/o\nlJ/olJ/YJXKudDuuiIjERVccIiISFxUOERGJiwqHiIjERYVDRETiosIhIiJxUeEQEZG4qHCIiEhc\nVDhERCQuKhwiAYqs5udmdmLYsYjUFBUOkWCNAD6I/BZJCSocIgExsyZABnATZQqHmb1hZhdEXv/O\nzB4KJ0KRqtF6HCLBGQq87u7LzCzfzHq7ezbwa+A3ZnYk0ItaWCNapCbpikMkOCOA5yOvn49s4+4L\nAAN+BAx39+JwwhOpGhUOkQCYWUtKl/T8R6TpeeBqK3UKcDRQ6O67wopRpKpUOESCMQyY6+4FAO6e\nA2wGzgWepbQbK9/MBoUXokjVaD0OkQCY2RtAT2BnmeZjgfXAbe4+z8zOBf7g7meFEaNIValwiIhI\nXNRVJSIicVHhEBGRuKhwiIhIXFQ4REQkLiocIiISFxUOERGJiwqHiIjERYVDRETi8n9vdrk2GfNS\n6AAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def f(x):\n", " a = (x-0.3)*(x-0.3) + 0.01\n", " b = (x-0.8)*(x-0.8) + 0.04\n", " return 1./a - 1./b\n", "\n", "def f_prime(x):\n", " a = -2*(x-0.3)\n", " b = (x-0.3)*(x-0.3) + 0.01\n", " c = 2*(x-0.8)\n", " d = (x-0.8)*(x-0.8) + 0.04\n", " return a/(b*b) + c/(d*d)\n", "\n", "x = 0.5\n", "\n", "diff_exact = f_prime(x)\n", "\n", "delta_x = [] # contains the delta_x\n", "error = [] # contains the error between the analytical and the approximate values\n", "dx = 1.2\n", "diff_prev = 1000. # we store the previous approximate derivative to compute the relative error\n", "error_rel = 1000. # relative error that is used as criteria to stop our loop\n", "\n", "# let us stop after 100 iterations for security\n", "for i in range(100) : \n", " diff = central_diff2(f, x, dx)\n", " delta_x.append(dx)\n", " error.append(abs(diff-diff_exact))\n", " error_rel = abs((diff - diff_prev)/diff_prev)\n", " if error_rel < 1.e-6:\n", " break\n", " diff_prev = diff\n", " dx *= 0.5\n", " \n", "pylab.loglog(delta_x, error, '+-')\n", "pylab.xlabel('$\\Delta x$');pylab.ylabel('Error at x=0.5');pylab.grid(True);pylab.legend(loc='best')\n", "pylab.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Question 4\n", "\n", "Consider the integral:\n", "$$\\int_0^{2\\pi} x^2 \\cos(x) dx.$$\n", "\n", "Show how the absolute error between the exact solution and the numerical solution varies with the size of the integration step for the four numerical integration methods: trapezoid rule; Simpson's rule; composite Simpson's rule; and Weddle's rule. Show your result by plotting error against integration step, $dx$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Two successive integrations by parts give the exact integral:\n", "\n", "$$\n", "\\begin{aligned}\n", "\\int_0^{2\\pi} x^2 \\cos(x) dx & = \\left[x^2 \\sin(x) \\right]_0^{2\\pi} - \\int_0^{2\\pi} 2 x \\sin(x) dx \\\\\n", " & = \\int_0^{2\\pi} 2 x (-\\sin(x)) dx \\\\\n", " & = \\left[2 x \\cos(x) \\right]_0^{2\\pi} - \\int_0^{2\\pi} 2 \\cos(x) dx \\\\\n", " & = 4 \\pi - \\left[2 \\sin(x) \\right]_0^{2\\pi} \\\\\n", " & = 4 \\pi\n", "\\end{aligned}$$" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEQCAYAAACa+vIpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4VFX6wPHvSQ/pPaQTakIngYCigoCiIKCCIIK4Fn72\nuqugrAirwOruWlZBUTAJIGBbsWDXSJEOoSShhEB67z2TmfP74wY2uHRyZ1LO53n22czNzD1HJpk3\n957zvq+QUqIoiqIol8vK0hNQFEVR2jYVSBRFUZQrogKJoiiKckVUIFEURVGuiAokiqIoyhVRgURR\nFEW5IiqQKIqiKFdEBRJFURTliqhAoiiKolwRFUgURVGUK2Jj6QmYg7e3twwLC2vx81ZXV+Pk5NTi\n51Vannqv2gb1PrUue/bsKZJS+lzoeR0ikISFhbF79+4WP29CQgIjRoxo8fMqLU+9V22Dep9aFyFE\n+sU8T93aUhRFUa6ICiSKoijKFVGBRFEURbkiKpAoiqIoV6RNLrYLISYB4wBXYIWU8gcLT0lRFKXD\najVXJEKIlUKIAiHEoT8cHyuEOCKESBVCzAGQUn4hpXwAeBCYaon5KoqiKJpWE0iAWGBs8wNCCGvg\nHeAmIBK4UwgR2ewp85q+ryiKojRjMkn2Z5axJ71E97Faza0tKeUmIUTYHw4PAVKllGkAQoh1wEQh\nRAqwBPhWSrnXrBNVFEVppWobjGxNLeLnw/n8lFJAYWU9w8K9WDt7qK7jtppAcg6BQGazx1lADPAY\nMBpwE0J0k1K++8cXCiFmA7MB/Pz8SEhIaPHJVVVV6XJepeWp96ptUO/TpSurN7G/wMi+AiPJxUYa\nTOBgDX19rJkUZkd/n1rd/01beyA5KynlW8BbF3jOcmA5QHR0tNQjW1Zl4bYd6r1qG9T7dGFSSg7n\nVfJzSj4/phSwP7MMgEB3R+6M8WVUhB8x4Z7Y21ibbU6tPZBkA8HNHgc1HVMURekwGhpN7DhRzM8p\nBfyYnE92WS0A/YPdeWZMD0ZH+tHL3wUhhEXm19oDyS6guxCiC1oAmQZMt+yUFEVR9Fda3UDC0QJ+\nSi7gt6OFVNU3Ym9jxTXdvXn0+m6M6uWLr6uDpacJtKJAIoRYC4wAvIUQWcB8KeUKIcSjwPeANbBS\nSplkwWkqiqLoJq2wSrvqSMln98kSTBK8ne0Z368zoyL8GN7NG0c7892yulitJpBIKe88x/GNwEYz\nT0dRFEV3jUYTezPK+Ckln59S8kkrrAagl78LD4/oxuhIP/oFumFlZZlbVher1QQSRVGUjqCyzsDm\nY0X8lJzPr0cKKK0xYGstGBruxd1DQxkV4UewZydLT/OSqECiKIqis+yyWm2XVXI+29OKMRglbo62\nXN/Ll9ERflzbwxsXB1tLT/OyqUCiKIrSwkwmycHs8qZbVgWk5FYA0MXbiXuuCmN0hB9RoR7YWLem\n4iKXTwUSRVGUFlBn0LLKTwWPwsp6rAREh3ry/M29GBXhR1cfZ/NOSkow1ICdvu2LVSBRFEW5TAWV\ndfySUsBPKQVsSS2kzmDCyc6a63r6MDrCjxE9ffF0sjP/xCrz4cA6SPwIOg+A297TdTgVSBRFUS6S\nlJIj+ZX8lKxddSQ2yyq/IzqY0RbIKj+tsQGOfgv71kDqTyCNEBwD4SN0H1oFEkVRlPM4Z1Z5kBvP\njOnBqAg/IjpbLquc3P1a8Dj4CdSWgEtnuPpxGHAXeHc3yxRUIFEURfmDspoGfj2i3bLadKSQyqas\n8uHdWklWeXURHPhYu3WVfxCs7aDXOBgwA7qOBCvzXhGpQKIoSodXWWdg98lStp8oZkdaCQezyzGa\nJN7O9oxrLVnlRgMc+xES18DR78DUCAED4eZ/QJ/boZOnxaamAomiKB1Oea2BXSdK2HGimB0nSjiU\nXY5Jgq21oF+QOw9d15VREb70D3K3fFZ5frIWPA6sh+pCcPKBmAe1W1d+kRd+vRmoQKIoSrtXVtPA\njhMl7EjTgkdybgVSgp21FQOC3Xl0ZDdiwr0YFOLROmpZ1ZTAoc+0AJKzD6xsoMdYGDgDuo0G69aV\nvKgCiaIo7U5xVT07T5Sw40QJ29OKOZxXCYC9jRWDQjx4YlR3Yrp4MTDEHQfbVhA4AExGOP4L7FsN\nRzaCsQH8+sLYJdB3Cjh5W3qG56QCiaIobV5hZb12m6rpiuNofhUAjrbWRIV68MyYzsSEe9E/2M0y\nW3PPp+iYFjwOrIfKXHD0hOh7tVtXnftd0alN0kRtYy1OtiohUVEU5Qz5FXVsTys+fcVxqmquk501\nUWGeTBwQyNBwT/oGumNn0wrLkNSVw6HPtV1XWTtBWEP3MXDTq9otLJsrS2I0SRO/ZPzCsv3L6OnR\nk0XXLGqhiZ9dmwwkQggnYCnQACRIKddYeEqKougop6z29BXH9rRiThbXAOBib0N0mAd3RAczNNyL\nPgGurbd+lckEJ37TgkfKV9BYCz69YMzfoN9UcPG78iGkiV8zfmXZ/mUcKT1CV8dgrnEe0AKTP79W\nE0iEECuB8UCBlLJPs+NjgTfRGlt9IKVcAtwGfCql/EoIsR5QgURR2pHMkprTVxs7ThSTWaIlAbo6\n2DCkiyd3xYQyNNyLyABXrC29q+pCStIgcS3sXwvlmeDgBgOmw8C7IGAQtEAio5SSXzJ/4d3973K4\n5DAhLiG87vwnwt77Hvtuv8HSO1rgP+TcWk0gAWKBt4H4UweEENbAO8AYIAvYJYT4Eq13+8GmpxnN\nO01FUVqSlJKMkhp2pJWw4UA9L2z/5XT2uHsnW4aEefKnq7oQE+5JL/82EDgA6qsgeYO26yp9KyCg\n6/UwZgH0HAe2LZPMKKUkITOBZfuXkVKSQrBLMEv6zqX/ukQqv3wfQkPwnDmjRcY6n1YTSKSUm4QQ\nYX84PARIlVKmAQgh1gET0YJKEJAItNLrWEVRzkZKyYmi6v9ecaSVkFdRB4CLHVzT05vZ14YTE+5J\nD18Xy+dxXCwpIf13LXgkfQGGavDsCqNehH7TwC2wBYeS/Jb1G0sTl5JSkkKQcxB/u2oh1x4wUvTw\nP6msrsbroQfxfvBBrOztW2zccxFSSt0HuVhNgeTrU7e2hBCTgbFSyvubHs8EYoDn0K5e6oAtZ1sj\nEULMBmYD+Pn5Ra1bt67F51tVVYWzs5nLQiuXRb1XliOlJLdacrjEyOESI0dKTZTXa587bvaCnh5W\n9PK0pqenNa6yBheXtvU+2dcV4J/3K/55v+BYl0ejtSMFvsPJ8x9FhWuvFrl1dYqUkqTaJDaWbySz\nIRMvGy/Guo1laGUIHmvXY3f0KA1du1Jx110YAzpf8XgjR47cI6WMvtDzWs0VyaWQUlYDf7rAc5YD\nywGio6PliBEjWnweCQkJ6HFepeWp98p8TCbJsYKq0+sbO0+UUFTVAIC/qwMjIjyJ6eJFTLgn4d5O\nZxQ7bDPvU0MNHP5a27Z7YhMgocu1MGABNhG3EGDnREALDielZHP2ZpYmLiWpOIlA50AWRi9kXPAN\nlK+IpfjdJQhHR3wXLsB98mSElXlv1LT2QJINBDd7HNR0TFGUVsJgNHE4t5JdJ0tOB47SGgOglVe/\ntrsPMeGeDA33IsSzk+Wq5F4pKSFrlxY8kv4D9RXgHgIj5kD/O8EjVIchJVuyt7Bs/zIOFh0kwCmA\nl4a9xIRuEzDsSSTziTtoSEvD9eab8Zs7Bxsfnxafw8Vo7YFkF9BdCNEFLYBMA6ZbdkqK0rGV1TSw\nL6OMPeml7E4vYX9mObUGbc9LsKcjoyL8GBruRUwXT4I9O1l4ti2gIgf2NzWJKj4Gtp0gcqKWMBh6\nNejw17+Ukq05W1mWuIwDRQfo7NSZ+cPmM7HrRKwqq8l/cQHln36GbWAgwe8vx/maa1p8Dpei1QQS\nIcRaYATgLYTIAuZLKVcIIR4Fvkfb/rtSSplkwWkqSocipSStqJo96aXsTS9lT3opxwq0rHFrK0Hv\nAFemDg4mOsyDQSEeBLg7WnjGLcRQp5UpSVyjlS2RJggZBlc/Ab0ngb2LLsNKKfk953eW7l/KgUIt\ngLw47EUmdZ2EjZUNFV9/Tf7iJRjLy/G6/z68H34Yq06WD9atJpBIKe88x/GNwEYzT0dROqQ6g5ED\nWeXsTi85HThO3aZyc7QlKtSDSQMDGRTiQf9gNzrZtZqPkJZRmQfb3oG98VBXBq6BMPxpLe/Dq6tu\nw0op2Za7jaWJS9lfuB9/J3/+OvSv3NrtVmytbWnIyCBzwUKqt27FoV8/QlauwKFXL93mc6na2U+B\noiiXIr+iTrtFdbKUPRmlJGWX02jSdlSF+zgxJtKPqFAPokI9CPd2bjtbcS9VSRpsfVO7fWVq1G5d\nDbobulyna5MoKSXbc7ezbP8y9hXsw6+TH/Ni5nFr91uxs7ZDGgwUrXifonfeQdjY4PfXeXhMm4aw\nbl31wlQgUZQOotFo4nBeJXszmgJHeunpxD97Gyv6B7sz+9pwokI9GBjigafTldV7ahNyD8DWN7TF\ncysbbd3j6sfBM1zXYaWU7MjbwbLEZewt2ItvJ19eiHmB27rfhp219u9es28feS/Op/7YMVzGjMFv\n3gvY+l15GRU9qECiKO1Uea2BfRlNaxsZpSRmlFHdoC2K+7naEx3qyX3DuxAV6kFEZ9fWWdxQL+m/\nw+Z/QeqPYOcCVz0GQx8GF3/dh96Zu5N3Et/RAoijL3OHzOX2Hrdjb60lDhorKih4/XXK1q3Hxs+P\noKXv4HL99brP60qoQKIo7YCUkvTiGnY3rWvsTS/laEElUmqL4hGdXZgcFcSgUA+iwzwJcHNou9tw\nL5eUcOwHLYBkbodOXnD9PBh8Pzh66D78rrxdLE1cyu783fg4+jBnyBwm95h8OoBIKan8/nvyXnkF\nY3EJnnfPxPuxx7F21rcEfEtQgURR2qA6g5FD2eVnBI7iai3pz9XBhkGhHozv15moUA/6B7vjZN+B\nf9WNjdqtqy2vQ0ESuAXDTa9p3Qbt9N/xtDtvN0v3L2VX3i68Hb3/J4AAGLKzyVv4N6p++w2HyEiC\nl72LY5/eus+tpXTgny5FaTsKKutO76LanV7KoexyDEZtUbyLtxMjevoSHaYtinfzaceL4pfCUKdt\n3/39LSg9Cd49YdK70HeyWVrV7snfw9LEpezM24mXgxfPDn6WKT2m4GDz34KNsrGRklWrKXzrLRAC\n3znP4TljBsKmbX00t63ZKkoHYDRJjuRVsqdpfWN3esnpMup2Nlb0D3Lj3uFdiA71ZFCIO17O+hfl\na1PqKmD3Cti2FKoLIDAKbngFet6sS/LgH+3N38vS/UvZkbsDLwcv/hL9F6b0nIKjzZk5NrUHD5E7\n/0Xqk1NwHjEC/7/Owzaw5Qo7mpMKJIpiYZV1htOZ4nszStmXUUZVfSMAPi72RId6MGtYGFGhHvQO\ncOtYi+KXoqoQdiyDnR9AfTmEj4RrVkDYNS1aOPFc9hXsY2niUrbnbsfTwZM/R/+ZO3re8T8BxFhV\nTeFbb1K6eg02Xl4EvvkmLjeMadNrViqQKIoZnVoUPxU09qSXciRfWxS3EtDL35VbBwaezt0I8nBs\n0x8wZlGWAb//W0sibKyHyAlw9ZMQOMgswycWJLI0cSnbcredDiBTekyhk+3/rr9U/vwzeX97mcb8\nfDzunIbPU09h7aJPlrw5qUCiKDo6lSl+KnA0XxR3sbdhQIg7N/U5tSjuhouD/vfu242CFNjyBhz8\nBIQV9J+qBRDv7mYZfn/hfpYlLmNrzlY87D14OupppvacetYAYsjPJ//ll6n88Sfse/Qg6I3XcRyg\nfwtcc1GBRFFaUG55LXtO7aTKKDszU9zbiZG9fBkUol1tdPdVi+KXJWu3toX3yDdaAcWY/4Nhj4Bb\nkFmGP1B4gKX7l7I1WwsgT0U9xbSe084aQKTRSOlHayl84w2k0YjPM0/jdc89CNv29QeDCiSKcpkM\nRhPJORVnXG3klGud/hxsregf5M4D14YTFeLBQLUofmWk1IonbnkdTm4GB3e4bg4MmQ1OXmaZwsHC\ngyzdv5Qt2Vtwt3fnyUFPcmevO88aQADqUlLInf8SdQcO4HT11fi/NB+74OCzPretU4FEUS5ScVU9\ne5stih/IKqPOYAK0vhtRYZ48EOJ+OlPc1lotil8xkxFSvtICSG4iuHTWdmBF3QP25umkeKjoEMv2\nL2NT1ibc7N14YtAT3NnrTpxsz54oaKqpofDtdyiJi8Pa3Z2Af/wD13E3t+u1LhVIFOUsjCbJsYLK\n07ep9mWUcaKoGgBba0HvADemDwklKtSDQaHudHZrJ+XTW4vGBjiwXquDVZyq9T6f8G/oNxVszHNl\nl1ScxLLEZfyW9Rtu9m48PvBxpkdMP2cAAajatIm8BQsxZGfjPmUyvs88g7W7u1nma0ltMpAIISYB\n4wBXYIWU8gcLT0lp4yrqDCQ2u9pIzCijsmkLrrezHYNCPJg6OJioUA/6BrrhYNu6qq+2F9aNtVoZ\n99/fhsoc8O8HU2IhYoKuVXibO1h4kOUHlpOQlYCrnSuPDXyM6b2m42x37iugxsJC8hcvpmLjt9h1\n7Uro6lV0ir5gq/N2w+yBRAixEhgPFEgp+zQ7PhZ4E62B1QdSyiXnOoeU8gvgCyGEB/APQAUS5aJJ\nKTnZtAX3j3WprAT09HdlwoCA01tw23R72LaipgR2Lmfo9rehsVLL/Zj4b+g6yiw5ICZpYlPWJmKT\nYtmTvwcXOxceGfAId0XchYvdubfnSpOJso8/oeCf/0TW1+P9+GN43X8/VnYdoHJyM5a4IokF3gbi\nTx0QQlgD7wBjgCxglxDiS7SgsvgPr79XSlnQ9PW8ptcpyjnVNhg5UmIkOSGVvU27qUpObcF1sGFg\niAc391VbcC2iIke7+tgTC4Zqyr2G4D1pEQQPNsvw9cZ6vkn7hrikONLK0/B38ucv0X/h9h63n/cW\nFkD9sWPkvjif2n376BQTg/9L87Hv0sUs825tzB5IpJSbhBBhfzg8BEiVUqYBCCHWAROllIvRrl7O\nILQ/D5cA30op9+o7Y6WtySlrvgW3lOSciqYtuEcI93Hi+l6+p682VF0qCylK1dY/9q/T2tj2nQJX\nP8GhlAJGmCGIlNeX8/GRj1mTsobiumJ6evRk8TWLuTHsRmytzv+HhKmujqJl71K8YgXWzs50XrwY\nt0kTO/RVa2tZIwkEMps9zgJizvP8x4DRgJsQopuU8t0/PkEIMRuYDeDn50dCQkLLzbZJVVWVLudV\nLl6jSZJRYeJYmYnUMiOppSZK67W8DTsrCHe3YmyYDYEODfTxd8LFDqAUakrJSYGcFItOv8NxrjxO\nSMan+BRuw2RlS57/GDKDJ1Hn6AcpBbr/ThU3FvNrxa9sq9pGg2ygl0MvpvlOo6dDT0SGYGvG1vO+\n3i4lBZeP1mJTWEjt0KFUTr6dHGdn+O033ebcFrSWQHJJpJRvAW9d4DnLgeUA0dHRcsSIES0+j4SE\nBPQ4r3JuRVX1pxs17U0v5UBWOfWN/92CO7yXB1Eh7kSFetKrs8vpLbjqvbIgKeHkFtjyLy0XxN4V\nhj+F9dCHCHT2pXmZQr3ep+TiZGIPxfJDzg8IBDd1uYlZvWfR07PnRb2+saSE/CVLqPjyK+xCQ/GP\n/RCnoUNbfJ5tVWsJJNlA80ydoKZjSgfWvAruvqbgkV5cA/x3C+6MoU1bcEM88HdzuMAZFbMymeDo\nd1oAydoFTr4w+iWIvhcc3HQfXkrJluwtxCbFsjNvJ062TsyImMGMyBn4O11cJ0QpJeWf/4eCV1/F\nWFOD98MP4fV//4eVvUouba61BJJdQHchRBe0ADINmG7ZKSnm1rw17N6MMvZllJ5uDevtbM+gEHem\nDwkhKtSDPmoLbutlNMChz7Q6WIUp4B4K4/6p9UO31T/fxmA08M0JbQE9tSwV306+PB31NJN7TD7v\nDqw/qk87Qd78+dTs2oVjVBSdF7yEfbduOs687bLE9t+1wAjAWwiRBcyXUq4QQjwKfI+2U2ullDLJ\n3HNTzEdKSVpR9entt3vSSzlWUAX8twrubYOCGBTqTlSIJ8Geqgpuq2eohX2rYetbUJ4Bvr3htg+g\n961grf9HTUVDBZ8e/ZQ1yWsoqC2gm3s3Xhn+CjeF3YTtJTSyMjU0ULz8fYrfew/h6Ij/wgW4T56M\nMEMvk7bKEru27jzH8Y3ARjNPRzGTmoZGEjPLzui7UVZjAP7bGnZCfy13o1+wO84duTVsW1NTArtW\nwM73oLoQgmPg5tegx41myQHJq85jVfIqPjv2GdWGamI6x7Dg6gVcHXD1Jf/xUXfkCDl/eZb6o0dx\nHTcOv7lzsPH21mnm7Yf6bVVanJSSrNLa0/029maUkpJbibGpCm43X2duiPQ7vQU33FttwW2TStK0\nLoT7VkNjLXQbA8OfgtCrzBJAjpQcITYplu9OfIdEckPYDdzT+x4ivSIv+VzSaKQkNpbCN97Eys2N\noHeX4aI2Z1w0FUiUK1bfaORQdsXpW1R7MkoprKwHoJOdNQOC3Xl4RFcGNVXBde/UsbJ+253MnVoj\nqZSvtN7n/e6AYY+Cb4TuQ0sp2Za7jdhDsWzL3YajjSPTek1jZuRMApwDLuucDVnZ5M6ZQ83u3biM\nGY3/woXYeHi08MzbNxVIlEtWUFF3Roe/Q9kVNBi1LbjBno5c3dWrqZihBz39XLBRVXDbPpMRjmzU\nAkjmDq2M+zVPa2XcXS5uB9SVMJgMfHfiO+KS4jhSegRvR2+eGPQEU3pMwc3+8naASSkp37CB/L+9\nDKASC6+ACiTKeTUaTRzOqzwjUzyrtBYAOxsr+gW6cc/VYQwK0arg+rqoLbjtSkMN7P9IK6RYkqbt\nwLrpVW0HlhnKuFc1VPHZsc9YlbyK/Jp8wt3CWXjVQsaFj8PO+vKvbBtLS8mb/xKVP/yAY3QUAUuW\nYBdknsZY7ZEKJMoZSqsb2JdZejpw7M8sp9agbcH1c7UnKtSDe64KY1CoB70DXLG3UVtw26WqQtj1\nPux8H2pLIGCQVoW31y1m2YGVX53PmsNr+OTIJ1QZqoj2i+bFYS8yPHA4VuLKrnCrNm0i54UXMJaV\n4/vnZ/D8058Q1urn+EqoQNKBmUyS1MKqM9Y20gq1nhvWVoLIzq5MHRzMoFAPBoW4E+iutuC2e0XH\nYNvbkLgWjPXQ82a46jEIGWaWBfSchhxe2PICG09sxCRNjAkdwz2976GPd58Lv/gCTDU15L/2GmVr\n12HfvTsh77+PQ69eLTBrRQWSDqSyzsD+zPL/tobNKKWyTuu54dHJlqhQD24fFKRtwQ1yo5Od+vHo\nEKSEjG3a+seRjWBtDwOma33QvbubYXjJzrydfJj0IVtzt+Jo48gdPe5gRuQMgl1apjVt7f795Dz7\nHA0ZGXj+6U/4PPmEyk5vQeqTop2SUpJR0qznRkYZR/IqMEntD8sevi6M7xfQVF7EnS7eTupqo6Mx\nNsLhr7QAkr0HHD21PuiD7wdnH92HbzQ18sPJH4hNiiWlJAVPB0/GuY1jztg5uDu0TFdBaTBQ9O57\nFL37Lja+voR8+CFOQ89XD1a5HCqQtBN1BiMHssr/m7uRXkpxU88NZ3sbBoa4M+b67kSHejAgxB1X\n1XOj46qvgsQ12gJ6WTp4hsO4f0H/O8Guk+7D1xhq+PzY56xKXkVOdQ5hrmHMHzafW7rewrbN21os\niNSnnSDnueeoO3gQt4kT8HvhBaxdXVvk3MqZVCBpo3LLa5sCRhl7MkpJyi5v6rkBXbyduK6nz+mE\nv+6+LlirhD+lMg92Ltey0OvKIHgo3LgIet5klja2RbVFrElZw/oj66lsqGSg70CeG/IcI4JHXPEC\nenNSSkrXrqXg1dewsrcn8I3XcR07tsXOr/wvFUjaAIPRRHJOxekF8X3ppeSU1wFgb2NF/2B3Hrg2\nXNuCG+KOl7O696s0U3AYtv0bDnysFVSMuEVbQA8eYpbh08rSiEuO46vjX9FoamRUyChm9Z7FAN8B\nLT6WoaCA3BfmUb15M05XX03nRYuw9fNt8XGUM6lA0goVV9Wz91RNqvRSDmSXUWfQEv4C3BwYFOrB\nA02l0yM6u2JnoxL+lD+QEk5u1tY/jv0ANo4waBYMfQi8uppheMme/D3EJsXyW9Zv2Fvbc1v325gZ\nOZNQ11Bdxqz4/gfyXnwRU309fn+dh8f06Wrdz0xUILEwo0lyNL/yjLWNk816bkQGuDF9SFPPjVB3\nOrvpX4ZbacOMBkjeAL+/Bbn7wckHRs6DwfdBJ0/9hzcZ+SnjJ2IPxXKo+BDu9u481P8hpvWahqeD\nPuMbKyvJf/kVyjdswKFPHwJe/Tv24eG6jKWcnQokZlZeayAxs6yp50Yp+zLKqKrXtuB6O9sxKMSD\naU09N/qqnhvKxaqvhL3xsH0ZlGeCdw+45S3oNxVs9a82UGOoYcPxDcQnxZNVlUWwSzDzYuYxodsE\nHG30++OneudOcubMoTG/AO+HH8b7oQcRtmojibm12UAihHACfgNeklJ+ben5nI2UkhOnem5kaAvj\nRwsqkVLrudHT35VJAwMYFKItiod4dlKX4sqlqciBHe/C7lioL4fQ4XDzP6D7DWCG/hlFtUV8lPIR\nHx/9mPL6cvr59OOZ6GcYGTwSax0X8E0NDRS+8SYlH36IbUgwYWtW4zig5ddclItjicZWK4HxQIGU\nsk+z42OBN9EaW30gpVxygVM9B3ys20QvQ22Dkf1Z/13b2JtRSmlTzw0XBxsGhXgwrl9nokI96K96\nbihXIu+QloF+8BOQJoicBFc9CoFRZhn+eNlx4pPjTy+gjwweyazesxjoO1D3P4aa9wxxnzYVv2ef\nxaqT/tuWlXOzxCdZLPA2EH/qgBDCGngHGANkAbuEEF+iBZXFf3j9vUB/IBmwWIVAKSVFtSa+3J9z\nusRIcm7F6Z4bXX2cGB3x354bXX1Uzw3lCkkJab9qC+jHfwFbJxj8AAx9EDzCzDC8loEemxTLluwt\nZllAP2N81TOk1bJEh8RNQoiwPxweAqRKKdMAhBDrgIlSysVoVy9nEEKMAJyASKBWCLFRSmnSc971\njUaSciq7jq3aAAAgAElEQVROX2nsSS8lv6Ie2IejrdZz48HrwokK9WBgsAceTqrnhtJCGhsg6XMt\ngOQfAmd/GDUfov8Ejvr3zTCYDHx/8nvikuI4XHIYTwdPHhnwCFN7TsXDwTx9OwzZ2eTMmUvNrl1a\nz5AFC7Dx1H/zgHJxhJTS/INqgeTrU7e2hBCTgbFSyvubHs8EYqSUj17gPPcARWdbIxFCzAZmA/j5\n+UWtW7fukudZUmfih5ONHC8zcqLCRGNTqPJ2FHRztyLYsZHefo4Eu1iphL9WrqqqCmdn/cuetyTr\nxmoCcr4nKOtr7BuKqe4UQmbwJPL9rkVa6b+gXGuq5ffK30moTKDMWIafjR/Xu17PYOfB2Ap9xv+f\n90lKHHbswGXdegAq77iDumFDzVJAUoGRI0fukVJGX+h5bfomvZQy9jzfWw4sB4iOjpYjLuMSOLOk\nhjmbf6NPoCt/6uvRVJfKA19X7Y5aQkICl3Nexfza1HtVlqktoO+Jg4ZK6HIdXPU4Tt1G0UsI9K5X\nm1uVy+qU1ad7oA/2H8w9ve9pkRLuF9L8fVI9Q9qO1hJIsoHmZT6Dmo5ZVJCHIwdeukFtwVXMIydR\nW0A/9Ln2uM/t2gJ65/5mGT65OJnYpFh+OPkDADeE3cCs3rPo7dXbLOM3V7V5MznPP696hrQRrSWQ\n7AK6CyG6oAWQacB0y04JhBAqiCj6Mpkg9SctgfDkZrBz0bLPYx4E95YpoX7e4aWJLdlbiEuKY2fe\nTpxsnbgr4i5mRMygs3Nn3cf/H/X15C1cSOlHa7Hv3o2Q5ctxiNC/F7xyZSyx/XctMALwFkJkAfOl\nlCuEEI8C36Pt1FoppUwy99wUxWwMdXDwY60Cb+FhcAmAMX+DqFngcHk9yC9FvbGer49/TXxyPGnl\nafh28uWZqGe4vcftuNi56D7+2dQeOIDXokWU5hfgec89+Dz1pOoZ0kZYYtfWnec4vhHYaObpKIp5\n1ZRo1Xd3LofqAvDrC7e+B71vAxv9d/qV1ZWx/sh6Pjr8ESV1JfTy7MXiaxZzY9iN2JphAf9spMFA\n0XvLKVq2DOHqSkjshzgNHWqRuSiXp7Xc2lKU9q34OGxfCvvWQGMtdButVeDtcp1ZdiBlVGQQnxzP\nhtQN1BnrGB44nFm9ZxHjH2PRagr1J06Q89wc6g4cwHXCLaSOGEEfFUTaHBVIFEVPGTu09Y/D34C1\nLfS9Q2th6xdpluETCxKJTYrll4xfsLGyYXz4eO6OvJtuHt3MMv65SCkpW7eO/L+/irC3J/D1f+F6\n000cS0iw6LyUy6MCiaK0NJMRDn+tJRBm7QIHd7jmaRgyG1z8dR/eaDLyS+YvxCbFcqDwAK52rtzf\n936mR0zH29Fb9/EvxFBQQO68eVRvOtUz5BVs/fwsPS3lCqhAoigtpaFau3W1/R0oPamVLbnpNRh4\nF9g56T58jaGGL1K/YFXyKrKqsghyDmLukLlM6jaJTratoxZVxfc/kDd/PqbaWtUzpB1RgURRrtQf\nW9gGDYExC6HXeLO0sC2sKWTt4bWsP7KeioYK+vv0N0sF3kuheoa0byqQKMrlyk/Wtu8ePNXCdjwM\newxCYswyfGppKnHJcXyT9o3uLWyvRPXOneTOmYshLw/vhx/C+6GHVM+QdkYFEkW5FFJCWoKWgZ76\nE9h2MnsL2+2524lLjmNr9lYcrB24vfvtzIycSYhriO7jXwpTQwOFb75JycqmniEfrVE9Q9opFUgU\n5WKcrsD7NuQfBCdfuH4eRJunha3BZOC7E98RlxTHkdIjeDl48djAx7ijxx24O7jrPv6lOqNnyNSp\n+D37F6yc9F8nUixDBRJFOZ/aMtgTqxVRrMwFn14w4W3odwfY6J91XdlQyadHP2V1ymoKagoIdwtn\nwVULGBc+Dnvr1pf1rfUMiaPwjTdUz5AORAUSRTmb0nQteOyNh4YqLXFwwr+1REIz7DLKqcrRKvAe\n/Yyaxhpi/GOYP2y+WSrwXq7mPUOcR4+i88KFqmdIB3HBQNLUvfDvUso/m2E+imJZ2Xu021fJX4Cw\n0irwDnsUOvczy/BJRUnEJcXxQ7pWgffGsBuZ1XsWkV7mSWC8HFJKyj//D/mLF4OUdF60CLdbJ6lt\nvR3IBQOJlNIohBhujskoikWYTHD0O20BPX0r2LtqwSPmQXAL1H94aWJz1mZik2LZnb8bJ1snZkTM\n4K6IuyxTgfcSGAoKyHtxPlUJCapnSAd2sbe29jX1UP8EqD51UEr5uS6zUhRzMNTC/rXaFt7iVHAL\nhhsXwcCZ4OCq+/D1xnq+Ov4V8cnxnCg/gb+TP3+O/jO3db/NYhV4L0XFxo3kLViIqa4Ov7lz8Jg5\nE2HVOm+7Kfq62EDiABQD1zc7JgEVSJS2p6oQdn0Au96HmmLoPABuXwGRk8Ba/2XD0rpS1h1Zx7rD\n6yipKyHCM4Il1yzhhrAbLFaB91I0lpaSt2Ahld99h0O/fgQsWaySCzu4i/qtkVL+Se+JXAohhBXw\nN8AV2C2ljLPwlJS2oOiYdvtq/zporIMeY7UKvKFXm2UB/WT5SVYlr2LD8Q3UG+u5NuhaZkXOYrD/\n4DaznlD588/kvjgfY0UFPk8+idf99yFs1J6dju6ifgKEEEHAv4Grmw5tBp6QUmZd6oBCiJXAeKBA\nStmn2fGxwJtoja0+kFIuOc9pJqK14y0GLnkOSgciJaRvpc/BlyFhF1jbQ/9p2hqITw8zDC/ZW7CX\n2KRYfsv8DRsrGyZ0ncDdkXcT7t52/oo3VlSQ/8oiyjdswL5XL0JWrsChZ09LT0tpJS72T4kPgY+A\nKU2PZzQdG3MZY8YCbwPxpw407Qx7p+l8WcCupjUZa2DxH15/L9AT+F1K+Z4Q4lPg58uYh9KeGRsh\nZYNWgTdnH662rnDdHBh8Pzj76D58o6mRnzJ+Iu5QHIeKD+Fu787sfrOZ1mtaq6jAeymqNm8hd948\nGouK8HroQXweeghhp38TLqXtuNhA4iOl/LDZ41ghxJOXM6CUcpMQIuwPh4cAqVLKNAAhxDpgopRy\nMdrVyxmaWvQ2ND00Xs48lHaqvlLL/dj+LpRngGdXGP8628sDuXbkjboPX22o5j/H/sPqlNVkV2UT\n6hrKvJh5TOg2AUcbR93Hb0nGqmoKXnuNsvXrsevalbC3/41j376WnpbSCl1sICkWQswA1jY9vhPt\ntlJLCQQymz3OAs5X+e5z4N9CiGuATWd7ghBiNjAbwM/PjwQdGuZUVVXpcl7l0tnXFRGY/TUBOT9g\nY6ymzC2SzD7PU+w1GKqsqKrV970qayzjt8rf2Fq5lVpZS7h9OPf73E9fx75Y5VmxI2+HbmPrwfbo\nUVzj47EuLqFmzGjyJ0wgs7gYdP55V79TbdPFBpJ70dZIXkfbrfU7YLEFeCllDXDfBZ6zHFgOEB0d\nLUfoUKYhISEBPc6rXILcA9oC+qHPQJogciIMewz3oCiaV6DS6706UnKEuKQ4vs34FhMmRoeMZlbv\nWfTzMU8CY0sz1dVR+PrrlMTFYxsSQsCbb9IpKsps46vfqbbpYjPbb5NSTtBxHtlAcLPHQU3HFOV/\nSQmpP2stbE/8BrZOMPgBGPqg1kxK9+ElW3O2EpcUx/bc7TjaODKt1zTuiriLIJe2m4xXm5hIzpy5\nNJw8icf06fj++RmsOrWOhlhK63axme13ol2N6GUX0F0I0QUtgEwDpus4ntIWNdbDwU+0EiaFKeDS\nGUYvgKh7wFH/CrgNxgY2nthIXFIcqWWp+Dj68OSgJ5ncYzJu9m66j68XU0MDRW+/Q/EHH2Dj50fI\nhytxGjbM0tNS2pCLvbW1VQjxNrCeMzPb917qgEKItcAIwLtp0Xy+lHKFEOJR4Hu0nVorpZRJl3pu\npZ2qKYHdK7UuhFX54NcHJr2r1cGy0X/3UHl9OR8f+ZiPDn9EUW0R3T2688rwV7gp7CZsrVt/AuH5\n1CUnkzNnLvVHj+J2+234zZmDtUvrz6pXWpeLDSSnutEsbHZMcmam+0WRUt55juMbgY2Xej6lHStJ\ng21LIXENGGqg6yi49T0IH2GWBMLMikxWpazii9QvqG2s5aqAq3hl+CsM6zyszSQQnos0GCh6/32K\nli7D2sNdlXtXrsjFrJFYAcuklB+bYT6KAhk7YNu/IeVrsLLRen8MewT8eptl+MSCROKT4/k542es\nhBXjuozj7t5308ND/wRGc6hPTSXnuTnUJSXhOm4cfvNewMbDw9LTUtqwi1kjMQkhngVUIFH0YzLC\n4a+19Y+sneDgDsOfgiGzwVX/CrhGk5FfM38lLimOxMJEXOxcuLfPvdzZ6058O/nqPr45aE2nYil8\n8y2snJwIfOMNXMfqn1ujtH8Xe2vrJyHEn/nfNZISXWaldBwN1bBvDWx/B0pPgnso3PQaDJgO9s66\nD19jqGHD8Q2sSl5FZmUmgc6BzBkyh1u73Uon2/azY6khPZ2cOXOp3bdPazr10kvYeLetDHul9brY\nQDK16f8faXZMAm2nWJDSulTmaYvnu1ZAXRkEDYYxC6HXeLCy1n34otoiPkr5iI+Pfkx5fTn9vPvx\n5KAnGRUyCmszjG8u0mSidO1aCv7xT4SNDQF/X4LrhAltfo1HaV0utvpvF70nonQQ+cla/4+DH4PR\nABHjYdhjEHK+QgYtJ7U0lfjkeL5O+5pGUyPXh1zPrN6zGOAzoN19uBqys8l5YR4127fjNHw4nV/+\nG7b+/paeltIOnTeQCCGelVK+2vT1FCnlJ82+t0hK+bzeE1TaASkhLUHLQE/9CWwcYdAsGPoQeHU1\nw/CSI7VHWP/TerZkb8HB2oHbut/GzMiZhLqG6j6+uWmtbz8nf5HW+tZ/4QLcp0xpd4FSaT0udEUy\nDXi16eu5aB0STxkLqECinFtjAyR9ri2g5x8EJ1+4fh5E3wedPHUf3mAy8N2J74hPjudwyWG8HLx4\ndMCjTO05FXcH/RMYLcGQX0Dui3+l+rdNdBoyhM6LFmEXpH+7YKVju1AgEef4+myPFUVTWwZ7YmHH\ne1CZAz69YMLb0HcK2DroPnxlQyWfHv2UNSlryK/Jp6tbV6Z7Tufpm5/G3tpe9/EtQUpJxdffkPfy\ny8j6evyefx6PGXep1reKWVwokMhzfH22x0pHV5oOO97Vyrg3VEGX62DCW9BttFkSCHOqclidsprP\njn5GTWMNMf4xzB82n6sDr2bTb5vabRBpLCkh76UFVP7wA479+9N5yWLsu6hlTcV8LhRI+gshKtCu\nPhybvqbpsf5/WiptQ/Ye7fZV8hcgrLTSJcMehc7mqYCbVJREbFIsP6b/iEBwY5cbmRU5iwivCLOM\nb0kVP/5I3vyXMFVW4vPM03jdey/Cuv3sOlPahvMGEiml+olUzs5kgqPfaR0IM34He1cteMQ8CG76\n35M3SRObsjYRmxTLnvw9ONs6MzNyJndF3IW/U/vfmWQsLyfvlVeo+PIr7CMjCIj9EIce7SPzXml7\nLjaPRFE0hlrYv1bbwlucCm7BcOMiGDgTHFx1H76usY6v0r4iPimekxUn6ezUmb9E/4Xbut+Gs53+\nCYytQdXmzeS+MI/G4mK8H3kE7wf/D2HbtotHKm2bCiTKxakqhF0fwK73oaYYAgbC7SsgchJY6/9j\nVFJXwvrD61l3ZB0ldSVEekXy6rWvMiZ0DDZWHePH2FhVTcHf/07ZJ59g370bQUuX4tjHPPXHFOV8\nOsZvoHL5Co9q+R/714GxHnrcBFc9BqFXmWUB/UT5CeKT4/nq+FfUG+sZETSCu3vfTbRfdIfKi6je\nsZPc55/HkJOD1/334f3YY1jZt8/NA0rb0yYDiRAiBHgLKAGOSimXWHhK7YuUkL5VW/84+h3YOMCA\nO2HoI+Cj/314KSW783cTnxRPQlYCdlZ2TOg2gZmRMwl361hVeUy1tRT863VKV63CNjSE0DVr6DRo\noKWnpShnMHsgEUKsBMYDBVLKPs2OjwXeRGts9cEFgkNf4FMp5WohxHpdJ9yRGA2QvEELILmJ0MkL\nRsyFwfeDk/4F/hpNjfyY/iNxSXEkFSfhYe/BQ/0fYmrPqXg5euk+fmtTs28fuXPm0pCejseMGfg+\n/ZRqfau0Spa4IokF3gbiTx1o6gv/DjAGyAJ2CSG+RAsqi//w+nuB7cCnQoh7gVVmmHP7Vleh5X7s\neBfKM8GrG4x/HfrfCbaOug9fbajms6OfsSZlDTnVOYS5hvHXoX9lQtcJONh0vF3mpoYGiv79b4pX\nrMTW35+Q2A9xGjrU0tNSlHMyeyCRUm4SQoT94fAQIFVKmQYghFgHTJRSLka7ejlDU0n7+U3n+hT4\nUN9Zt1PlWVrw2BMH9RUQejXc/Bp0vxHMkBGdV53HRykf8cnRT6gyVBHlF8WcIXO4Lvg6rETHzMiu\nTUoid84c6o+l4j5lCr7PPYu1c8fYjaa0Xa1ljSQQyGz2OAs4XznY74CXhBDTgZNne4IQYjYwG8DP\nz4+EhIQWmWhzVVVVupxXb86VaQRnfoFP4RaElBT4Xk1W0EQqXbtDLpC7SdfxMxsy+aXiF/ZW7wVg\nQKcBXO91PaH2oZAGm9JafvxW/14ZjTh9+y1OG7/F5OJCxaOPkN+nD0d277b0zMyq1b9Pylm1lkBy\nSaSUh4DJF3jOcmA5QHR0tByhQz/qhIQE9DivLqTUKu/+/hac2AR2zlry4NAH8XMPwU/34SVbsrcQ\nlxTHjrwddLLpxPSI6cyInEGgs/4JjK35vao/dkxrfZucjOstt+D/wvNYu7fPopIX0prfJ+XcWksg\nyQaCmz0OajqmXClDndb7Y9s7UHgYXAK0BlKDZoGj/h9W9cZ6vkn7hvikeI6XH8e3ky9PRz3N7T1u\nx9VO/wTG1kwajZR8+KHW+tbFhcC33sT1hhssPS1FuWStJZDsAroLIbqgBZBpwHTLTqmNqynRug/u\nfA+qC8GvL9z6HvS+DWzsdB++rK6M9UfWs/bwWorriunp0ZNFwxcxNmwsttYqC7v+xAly5z5PbWIi\nLmPG4P/SfGy8Ot7ONKV9sMT237XACMBbCJGFtmi+QgjxKPA92k6tlVLKJHPPrV0oPg7bl2p90Btr\nodsYuOpRrRKvGRL4MioyiE+OZ0PqBuqMdQwPHM6s3rOI8Y/pUAmE5yJNJkpXr6bgX68j7O0JeO01\nXMePU/82SptmiV1bd57j+EZgo5mn0z5ICZk7tPyPw9+AtS30u0MrouirfwVcKSWJhYnEJcXxS8Yv\n2FjZMD58PHdH3k03j266j99WNGRlkfv8C9Ts3InTddfSeeHfsPXztfS0FOWKtZZbW8rlMBkh5Sst\ngGTvBkcPuOYZGDIbXPRePgejycjPGT8TlxTHgaIDuNm7cX/f+5keMR1vR/0TGNsKKSVl6z+m4NVX\nQQg6v/Iybrfdpq5ClHZDBZK2qL4KEtdoC+hl6eDRBW7+BwyYDnZOug9fY6jhP6n/YVXyKrKrsgl2\nCeb5mOeZ2HUinWxV5nVzhtxccuf9leqtW3G6ahidX34Z24AAS09LUVqUCiRtSUWutni+eyXUlUNw\nDNz4CvS8Gaz0bx1TWFPIR4c/Yv2R9VQ2VDLAZwB/jv4zI4NHYm2G8dsSKSXlX2wgf9EiZGMj/vNf\nxH3aNHUVorRLKpC0BflJWgfCg5+ANEKv8VoF3uAhZhn+aOlR4pPi+ebEN5ikiVEho7g78m4G+A4w\ny/htTWNhIbkvzqfq119xjI4iYNEi7EJCLD0tRdGNCiStlZRw/BethPvxX8C2E0T/CYY+BJ76V8CV\nUrItdxtxSXH8nvM7jjaOTOkxhZkRMwl2Db7wCTqoio0byVuwEFNdHb5znsNz5kzV+lZp91QgaW0a\nG+DQp9oVSEESOPvB9X+F6Huhk6fuwxuMBjae2Eh8cjxHS4/i7ejN4wMf546ed+Bm76b7+G1VY2kp\neQsWUvnddzj060fAksXYh3eskvdKx6UCSWtRWwq7P4Sdy6EyF3wjYeJS6DsZbPRvYFReX84nRz9h\nbcpaCmoL6Obejb9d/Tdu7nIzdtb6JzC2ZZU//UTu/JcwVlTg89RTeN13L8JG/WopHYf6abe00pOw\nfRnsXQWGaggfCRPfhq6jzJJAmFWZxeqU1Xx+7HNqG2sZ1nkYC69eyFUBV6mF4QswlpeTv2gR5Ru+\nxD4igpCVK3Do2dPS01IUs1OBxFKydmv5HylfgrCCvlNg2CPg39cswx8oPEBcUhw/ZfyEFVbcHH4z\nd0feTU9P9UF4Mao2byb3hXk0Fhfj/fDDeD/4fwg7deWmdEwqkJiTyQhHvtUCSOZ2sHeDqx6HmP8D\nV/1zC4wmIwlZCcQlxbGvYB8uti7c0/sepveajp+T/gmM7YGxqoqCv/+dsk8+xa5bV8LeeQfHvn0u\n/EJFacdUIDGHhhrY/5GWQFiSBu4hMHYJDJwB9i66D1/bWMuXqV+yKmUV6RXpBDoH8tzg57i1+604\n2eqfwNheVG/fQe7zz2PIy8Pr/vvwfuwxrOz1X79SlNZOBRI9VRVoi+e7VkBtCQQMgimx0OsWsNb/\nn76otoi1h9fy8ZGPKasvo49XH1677jVGh4zGxkq99RfLVFNDwb9ep3T1auxCQwlds5pOAwdaelqK\n0mqoTxM9FBzW8j8OfAzGBi3z/KpHIWSYWRbQj5cdJz45nq+Pf43BZGBE8Ahm9Z7FIN9BagH9EtXs\n3UvO3LkY0jPwmDkT36efwspR/z72itKWqEDSUqTUOg9uexuO/QA2DjDwLhj6CHjrXwFXSsnOvJ3E\nJcWxOXsz9tb2TOo2iZmRMwlzC9N9/PbGVF9P4ZtvUfLhh9gGBBASF4dTTMtVEjAYDGRlZVFXV9di\n52wP3NzcSElJsfQ0OhwHBweCgoKwtb28XkGtPpAIIcKBFwA3KeXkpmNOwFKgAUiQUq6x2ASNBkj6\nj7aAnncAOnnDiOdh8H3gpH8FXIPJwPcnvyc+KZ6UkhQ8HTx5eMDDTO05FU8H/RMY26PagwfJmTOX\nhuPHcZ86Fd+//AVr55ZdS8rKysLFxYWwsDB1ldhMZWUlLi76rxsq/yWlpLi4mKysLLp06XJZ59A1\nkAghVgLjgQIpZZ9mx8cCb6I1sfpASrnkXOeQUqYB9wkhPm12+DbgUynlV0KI9YD5A0ldOcEZ/4E3\nH4aKbPDuAbe8Bf2mgq2D7sNXNlTy2dHPWJ2ymvyafLq4deGlYS8xvut47K3VAvDlkA0NFC5bRvHy\n97Hx9ib4/fdxvma4LmPV1dWpIKK0CkIIvLy8KCwsvOxz6H1FEgu8DcSfOiCEsAbeAcYAWcAuIcSX\naEFl8R9ef6+UsuAs5w0CDjZ9bWzhOZ9fWSbseBf2xNG1oRLCroHxr2udCK2sdB8+tyqX1Smr+ezY\nZ1QbqhniP4QXh73I8MDhWAn9x2+v6g4fJmfOXOoPH8Zt0iT8np+Ltau+PeVVEFFaiyv9WdQ1kEgp\nNwkhwv5weAiQ2nSlgRBiHTBRSrkY7erlYmShBZNEwDyfnjn7tPpXSf/RHve5jd12w4i+5T6zDJ9U\nlERcUhw/pP8AwI1hNzKr9ywivSLNMn57JRsbKf7gAwrfWYq1mxtBS9/B5frrLT0t3RUXFzNq1CgA\n8vLysLa2xsfHB4CdO3diZ+Hkyh07drBu3Tpef/31//leUFAQhw4dwt3d/Zyvv5jnnM8HH3zAoUOH\neOONNy7r9R2NJdZIAoHMZo+zgJhzPVkI4QW8AgwUQsxtCjifA28LIcYBX53jdbOB2QB+fn4kJCRc\n+kylCa/iPQRlfYFH2SEarR3JDbyFrKDx1Dv4UFVVdXnnvUgmaSKpNolfKn4htT4VB+HAdS7XcZ3L\ndXiaPCk4WEABZ7tgU/7obO+VdW4ubrFx2KanUxcdRcW0aeRYWYGO7+kpbm5uVFZW6j7OudjZ2bF5\n82YAFi1ahLOzM48//jgA9fX11NfXI6VESomVGa60TzEajVRWVhIZGcnChQvP+m8kpaSyshLr81RV\nvpjnGI3Gc36/rq6OhoYGi75H5lZXV3f5n2enflj0+h8QBhxq9ngy2rrIqcczgbf1nENUVJS8LBk7\npJzvKuU/I6Xc+paUtWVnfPvXX3+9vPNeQK2hVn585GM5/vPxsk9sHznq41Ey9lCsrKiv0GW8jqD5\ne2VqbJRFK1bKlL795JGYobJ840azzyc5OdnsY57L/Pnz5WuvvSallPLYsWMyIiJCTp8+XUZERMis\nrCz5wAMPyKioKBkZGSkXLFhw+nWBgYHy2WeflX369JFDhgyRx48fl1JKmZeXJ2+99VYZFRUlBw8e\nLLdt2yallPKGG26Q/fv3l/3795cuLi5y9erVsqamRt59992yT58+cuDAgfLbb7+VUkr5448/yokT\nJ0oppSwoKJCjR4+WkZGRcvbs2TIgIECWlpZKKaUcP368HDRokIyMjJTvv//+GXM79ZxTDAaDdHNz\nk0888YTs27ev/P3338943rZt2+SoUaOklFK+//778oknnjjvf097c7afSWC3vIjPWEtckWQDzRta\nBDUda32CBsOd66HbKLC+vG1xl6KkroT1h9ez7sg6SupKiPCMYMk1S7gh7AZsrfQfvyNoSE8nZ+7z\n1O7di/OoUXRe8BI23qq/fHOHDx8mPj6e6OhoAJYsWYKnpyeNjY2MHDmSyZMnExmp3VL19PTk4MGD\nrFy5kqeffpovvviCxx9/nGeffZahQ4dy8uRJxo8fz6FDh/j+++8B7dbZAw88wC233MJbb72Fvb09\nBw8eJCkpiZtuuonU1NQz5jN//nxGjhzJ888/z4YNG1i+fPnp78XFxeHp6UlNTQ3R0dHcfvvteHh4\nnPO/rbz8/9u787gqy/Tx458bOCyKQIpmgrmLCYiCYOTGzzStTEut1BawLLPRsmlobDF1jEnLKVuY\n0kpxyrSislzKlgHN0a+aZU6Ie06iKOIKhrJdvz8OHkERgQMcwOv9ep2Xnuc8z32u59wcLu5nue6T\n9P4j558AAB+ASURBVOnTp0KHrC61P+o8RySSTUAHY0wbrAlkJDDaAXFcnjEQMKja3+a3k7/x/rb3\n+XLPl5wtOEsf/z5Ed44mvHm4npCtKoWFHPtgERn/+AfGxYUWs2biNWRIrfh8py9LYdvBU1XaZucW\nXky9LbBS27Zr186WRAAWL17Me++9R35+PgcPHmTbtm22RDJq1CgA7rnnHiZPngzAd999x44dO2zb\nHz9+nJycHDw8PMjIyCA6OppPP/0ULy8v1q5dS2xsLACBgYE0b978okSyZs0aVq5cCcDQoUNLXB78\n6quv8uWXXwLWS6r37NlTIvYLubq6cscdd1To8yhrf5RVdV/+uxiIAnyNMWnAVBF5zxgzAViF9Uqt\n+SKSUp1x1EYiwubDm1m4bSHJ+5NxdXLltna3cX/n+2nroxMiVaW8Awfwee11Du/YQcNevbjmhRlY\nmjd3dFi1VsOG5++Z2bVrF6+99hobN27Ex8eHe++9t8RNlKUlYhEp9YR9fn4+d999NzNmzLAlInt8\n9913rFmzhv/7v//Dw8ODXr16XfYGTw8PjxIxu7i4UFhYCHDJbS+1P+q86r5qa9Qllq8EVlbne9dW\n+YX5fPe/70hISSDlaAo+bj6M6zKOkZ1G4uuhh1iqkohwIjGRjJmzsOTn0/xv0/G5885aMQoprrIj\nh5pw6tQpGjVqhJeXF+np6axatYpBg86P0j/66CP+8pe/sHjxYnr27AlA//79iY+P54knngBgy5Yt\ndO3aldjYWMLDwxkxYoRt+969e7No0SL69OlDamoqhw8fpn379hw8eNC2Tp8+ffjwww+ZPHkyy5Yt\ns50AP3nyJI0bN8bDw4OUlBQ2bdpU4f1r3bo1mzdvZsCAAXz66aelrnOp/VHn1fo72+uL03mn+XTn\npyxKXcTB0wdp5dWKKddP4bZ2t+HhokPkqpZ3+DDpU6Zwes0PNOjRg//ddhuBI4Y7Oqw6JzQ0lM6d\nO9OpUydatWplSxbnZGZm0qVLFzw8PFi8eDEA8fHxjB8/ngULFtjOq7z++uvMmTOHoKAgvvnGegn7\n3//+dyZOnMi4ceMIDg7GYrEwd+7ci/7ynz59OqNGjeKDDz6gZ8+etGhhnXLh1ltvZd68eXTu3JmA\ngAB69LjkxZ+XNG3aNB566CF8fHzo06dPqeuUtj/x8fEVfq/6zFhPzNdv3bt3lx9//LHK201OTiYq\nKqrMdQ6dPsSHqR+SuDORrLwsQpuFEh0YTVTLKL2BsBqICKeWLePQC3FIbi7NnnySq+4Zzeo1ay7b\nVzUpNTWV6667ztFh2MXeezVKoyVSHKe0n0ljzGYRufRJpyI6Iqkm249tZ2HKQr7+7WsKKWRAqwFE\nd44muGnNzIB4JcrPzCR92jSyv/sej27daPHi33Ft3drRYSlV72kiqUIiwn8O/oeElAQ2pG/Aw8WD\nkZ1Gcs919+DfyN/R4dVrp77+mkPT/0bh6dM0i42lcUw0poyb0ZT90tLSHB2CqiU0kVSB3IJcVuxd\nwb+2/YvdJ3bTzKMZk0InMaLjCLzdvB0dXr2Wf/w4h2e8wKmVK3EPCqLFzBdxa1/9ZfuVUudpIrHD\n6YLTzNs6j8XbF5OZk0nHqzoS1yuOm1vfjKUGbmC80mX9O4n0qc9TcPwETR9/jCZjx2IqOZ+CUqry\nNJFUwv5T+/nXtn/x2YHPyE3LpWeLnsT1iiPymshad2lpfVSQlcXhv7/Iyc8/xy0ggGvnzcO9jp+4\nVqou00RSAVsytrAwZSHf//49zk7OhDUI46kbn6LjVR0dHdoVIyspiUPTppOfmUmTR8bR9NFHMXqj\nmFIOpdefXkZBYQHf/u9b7ll5D/d9dR8bD23kweAHWTV8Fff63qtJpIbkHzvGgSf/Qtr4R3H28qL1\nksU0mzRJk4gd4uLiCAwMpEuXLnTt2pUNGzYAMHbsWLZt2+bg6M5LSEhg2rRpVd5uTEwMiYmJl19R\nXZaOSMqw49gOJiVNIi07DT9PP56OeJrb299OA0sDALZRe75s9ZWIcGrFSg7HxVGQnY3vhAn4PvyQ\nJhA7rV+/nuXLl/PTTz/h5uZGZmYmubm5gHUujvogPz8fFxf9FVcTdERShpaNWnKt17W8EvUKK+5Y\nwejrRtuSiKp+eYcPk/bonzj4l79gadmSNp8m0nTCnzSJVIH09HR8fX1xc7NOy+zr62u7YzwqKopz\nN/B6enoSGxtLYGAg/fv3Z+PGjURFRdG2bVtbscSEhASGDh1KVFQUHTp0YPr06QCcPn2aW2+9lZCQ\nEIKCgvjoo48A+P777+nWrRvBwcE88MADnD17FrCWK4mLiyM0NJTg4GC2b98OWOtjeXp6AvDJJ58Q\nFBRESEhIqXeiJycn07t3b4YMGULnzp3Zt28fQUG2Wb6ZPXt2qaObzZs307dvX8LCwhg4cCDp6el2\nf8ZXlPLUmq/rj0rPR3IZ1TUfyZWusLBQjn30kWwP6y6pIV0lc/4CKczPt6vN2tZXjp6PJCsrS0JC\nQqRDhw4yfvx4SU5Otr3Wt29f2bRpk4iIALKyaL6W22+/XQYMGCC5ubmyZcsWCQkJERGRBQsWSPPm\nzSUzM1P++OMPCQwMlE2bNkliYqKMHTvW1u6JEyckJydH/P39ZceOHSIict9998mrr74qIiKtWrWS\nl156SURE4uPj5cEHH7wo7qCgIElLSxMRuWi+ERFrPzdo0ED27t0rIiK//fabBAYG2l5/+eWXZerU\nqSIiEh0dLZ988onk5uZKZGSkZGRkiIjIkiVLZMyYMRX9SOu8ujYfiVKXlPv776RPeZ4/NmygQUQE\n18z4G66tWjk6rOr11WQ49N+qbbN5MNw885Ive3p6snnzZn744QeSkpK4++67mTlzJjExMSXWc3V1\ntRVpDA4Oxs3NDYvFQnBwMPv27bOtN2DAAJo0aQLAsGHDWLt2LbfccgtPPvkkf/3rXxk8eDC9e/fm\nl19+oU2bNnTsaD23GB0dTXx8PJMmTQJgyJAhAISFhfHZZ59dFHfPnj2JiYnhrrvuYtiwYaXuW0RE\nBG3atCnf5wTs2LGDX3/9lQEDBgDWmROvueaacm+v9ByJqiWkoIBj/3qfI6+9hnFxofn06fjcOQJT\ng9O8XmmcnZ2JiooiKiqK4OBgFi5ceFEisVgstkvanZycbIfCnJycyM/Pt6134WXvxhg6duzITz/9\nxMqVK3nuuee48cYbGTp0aJkxnWvf2dm5RPvnvP3222zYsIEVK1YQFhbG5s2bbQnsnOJl8IuXiYfS\nS8WLCIGBgaxfv77M2NSl1YlEYoxpCzwLeIvIiKJltwO3Al7AeyLyjQNDVHY4u2sXB599jjNbt+IZ\nFUXzaVOvrPlCyhg5VJcdO3bg5OREhw4dAGtp9FZ2jPy+/fZbjh07hoeHB0uXLmX+/PkcPHiQxo0b\nc++99+Lj48O7777LU089xb59+9i9ezft27fn/fffp2/fvuV+nz179tCjRw969OjBV199xf79+y9K\nJMVdffXVZGRkcPToUTw9PVm+fHmJMvgAAQEBHDlyhPXr1xMZGUleXh47d+4kMLD2lvevbao9kRhj\n5gODgQwRCSq2fBDwGtbJrd4VkUt+m0RkL/CgMSax2LKlwFJjzFXAbEATSR0jublkvvMOmW/PxdnT\nkxazZ+N16y16U2cNyM7OZuLEiZw4cQIXFxfat29fYgrbioqIiGD48OGkpaVx77330r17d1atWkVs\nbCxOTk5YLBbeeust3N3dWbBgAXfeeSf5+fmEh4fzyCOPlPt9YmNj2bVrFyLCjTfeSEhISJnrWywW\nnn/+eSIiIvDz86NTp04XrePq6kpiYiKPPfYYJ0+eJD8/n0mTJmkiqYjynEix5wH0AUKBX4stcwb2\nAG0BV+AXoDMQDCy/4NGs2HaJpbT/DyC0rBj0ZHvt88fWrbLntiGyLaCTpP35Sck7erRa36+29ZWj\nT7ZXpQULFsif/vSnKmnr1KlTVdKOqrhafbJdRNYYY1pfsDgC2C3WkQbGmCXAUBF5Eevo5bKM9c/W\nmcBXIvJT1UWsqlNhTg5H3nyTYwsScPH1xf+f/6RRv//n6LCUUnZw1DkSP2B/sedpwCWnNzPGNAHi\ngG7GmKeLEs5EoD/gbYxpLyJvX7DNw8DDYD1OmpycXLV7gPXwQHW0W19Zdu7E6/0PcDlyhD969SJ7\n+DAOOBmogc+wtvWVt7e3bcrYum748OEMHz68SvanoKCg3nwudc2ZM2cq/R2pEyfbReQo8MgFy14H\nXi9jm3nAPLDOkFgds+OVZ4ZEBQXZ2WTMns2JJR9hadmSaxIW0PD662s0htrWV6mpqToTYCl0hkTH\ncXd3p1u3bpXa1lGJ5ADQsthz/6Jlqp7JXr2a9KnTyM/IoHFMDE0fm4hTA60OoFR94qhEsgnoYIxp\ngzWBjARGOygWVQ3yjx/n8N9f5NSyZbi2b0fr1z7E4zJX2Cil6qaauPx3MRAF+Bpj0oCpIvKeMWYC\nsArrFVzzRSSlumNR1U9EyPr6aw7NeIGCU6fwffRRmjwyDietj6VUvVXttw2LyCgRuUZELCLiLyLv\nFS1fKSIdRaSdiMRVdxyq+uUdziBtwkQOPPFnLC1aWIssPjZRk0gtdejQIUaOHEm7du0ICwvjlltu\nYefOnY4Oy+aGG24AYN++fXz44YcV3l7L5Ndcmfw6cbJd1W4iwslPP+XwrJeQ3FyaxcbSOPp+jJbw\nrrVEhDvuuIPo6GiWLFkCwC+//MLhw4dtdbAcbd26dcD5RDJ6dPmPfmuZ/JqlhYyUXXL37+f3Bx4g\n/bkpuAcE0PaLpTR58AFNIrVcUlISFoulxF3lISEh9O7dGxEhNjaWoKAggoODbeXfk5OT6du3L0OH\nDqVt27ZMnjyZRYsWERERQXBwMHv27AGsfwk/8sgjdO/enY4dO7J8+XLAennpmDFjCA4Oplu3biQl\nJQGQkpJCREQEXbt2JTIykl27dgHYSsdPnjyZH374ga5du/Lqq69SUFBAbGws4eHhdOnShblz5160\nf7W1TP7UqVPrZ5n88ty1WNcfemd71SvMz5ejCQmS2rWbbA8Nk2OLF0thQYGjw7qk2tZXjr6z/bXX\nXpNJkyaV+lpiYqL0799f8vPz5dChQ9KyZUs5ePCgJCUlibe3txw8eFDOnDkjLVq0kOeff15ERObM\nmSOPP/64iFjLsw8cOFAKCgpk586d4ufnJzk5OTJ79mxbefbU1FRp2bKl5OTkyIQJE+SDDz4QEbGV\nohcRadiwoYhY++7WW2+1xTd37lyZMWOGiIicOXNGwsLCbGXjz6mtZfJff/11EamdZfJr9Z3tqv45\nu3s36c9NIWfLFhr26c0106dj0bLblTZr4yy2H9tepW12atyJv0b8tVLbrl27llGjRuHs7MzVV19N\n37592bRpE15eXoSHh9tKrLdr146bbroJsJaYPzfCALjrrrtsRSHbtm3L9u3bWbt2LRMnTrTG16kT\nrVq1YufOnURGRhIXF0daWho33XTTZe9l+Oabb9i6davt+P/JkyfZtWtXidLxtbVM/rnS9/WtTL4m\nElVukpfH0XffJfOfb+HUoAEtXpqF1223aZHFOigwMLBSJ2LPHSqCipeVv5TRo0fTo0cPVqxYwYgR\nI3jnnXfo16/fJdcXEd544w0GDhxYZqxaJt9KaqBMviYSVS45v6aQ/uyznN2xg0Y3D6L5c8/hUkb5\nblV+lR052KNfv34888wzzJs3j4cffhiArVu3cvLkSXr37s3cuXOJjo7m2LFjrFmzhpdfftl2TL88\nPvnkE6Kjo/ntt9/Yu3cvAQEB9O7dm0WLFtGvXz927tzJ77//TkBAAHv37qVt27Y89thj7N69m61b\nt5ZIJI0aNSpRNmXgwIG89dZb9OvXD4vFws6dO/Hz8yvxC1bL5J9XE2XyNZGoMhWeOUNmfDxH5y/A\npXFj/N98g0b9+zs6LGUnYwyff/45kyZNYtasWbi7u9O6dWvmzJlDr169WL9+PSEhIRhjeOmll2je\nvHmFEsm1115LREQEp06d4u2338bd3Z1HH32U8ePHExwcjIuLCwkJCbi5ufHxxx/z/vvvY7FY8PX1\nvehkcZcuXXB2diYkJISYmBgef/xx9u3bR2hoKCJC06ZNWbp0aYlttEz+eTVRJt9Yz6fUb927d5dz\nV2lUpdpWv6mq/fHjj6Q/+xy5//sf3iOGc/VTT+Hs5eXosCqltvVVamoq1113naPDqBYxMTEMHjyY\nESNGVHjb2lhrKyEhgR9//JE333zT0aFUq9J+Jo0xm0Wk++W21RGJukhB9mmOvPIPjn+4GIufH9fO\nf4+GRTeHKaXUhTSRqBKyf/iB9Oenkn/oEFfdfx/NJk3SIouqQhISEhwdQpWKiYm56CS9KkkTiQKs\nRRYzZs7i5Bdf4NquHa0+XESDSpaUVkpdWTSRKE59vYpDM2ZQcPIkTcY/gu/48VofSylVbppIrmB5\nGRkcnvECWd9+i3vnzlz77ju419MTwEqp6qOJ5AokIpz87HMOz5qFnDlD0yf/TJMxY7Q+llKqUmp9\n0UZjTFtjzHvGmMQLljc0xvxojBnsqNjqoty0A+x/cCzpzz6LW4cOtPliKb4PPaRJ5ArzxBNPMGfO\nHNvzgQMHMnbsWNvzJ598kldeeaXc7U2bNo3Zs2dftLx4UcHk5GQGDy776/rDDz9Uy4ntS8Wnqka1\nJhJjzHxjTIYx5tcLlg8yxuwwxuw2xkwuqw0R2SsiD5by0l+Bj6sy3vqsICuLY+9/wN4hQ8jZsoWr\nn59Cq/f/hVsFavao+qNnz562Mu2FhYVkZmaSknJ+brl169bZ5gOpC0SkRKkQVbOqe0SSAJS4X98Y\n4wzEAzcDnYFRxpjOxphgY8zyCx7NSmvUGDMA2AZkVG/4dUfhmTOc3bOH7NWrObZoEYdfepm0xx7n\nt2HD2dHjenaGR3A4Lo4GoaG0XfYljUePxjjV+gGpqiY33HCDrfZSSkoKQUFBNGrUiOPHj3P27FlS\nU1MJDQ0F4OWXX7aVbJ86daqtjbi4ODp27EivXr3YsWOHbfnmzZsJCQkhJCSE+Pj4Ut//9OnTPPDA\nA0RERNCtWze++OILwHoXtre3NwCrV6+ma9eudO3alW7dupUokwLW0U5AQAD3338/QUFB7N+/31aO\nHSAxMbHU0c2ePXsYNGgQYWFh9O7du0J37KvSVevxDBFZY4xpfcHiCGC3iOwFMMYsAYaKyItAeQ9T\nRQENsSaiHGPMShGp13+OSF4eeYcOkXfgAHlpaeSmpZGXVvT/A2kUHMkssb5xc8Pi54fFzw+vkC64\n+vvj1jGAhr16apFFRYsWLXBxceH3339n3bp1REZGcuDAAdavX4+3tzfBwcG4urryzTffsGvXLjZu\n3IiIMGTIENasWUPDhg1ZsmQJW7ZsIT8/n9DQUMLCwgAYM2YMb775Jn369CE2NrbU94+Li6Nfv37M\nnz+fEydOEBERQf/+/enRowf9i0rwzJ49m/j4eHr27El2djbu7u4XtbNr1y4WLlzI9ddfX+59f/jh\nh3n77bfp0KEDGzZs4NFHH+Xf//53JT5FdY4jDoz7AfuLPU8DelxqZWNMEyAO6GaMeVpEXhSRZ4te\niwEyS0sixpiHgYfBWtgsOTm5ynbgnOzs7Kprt7AQp1OncM48ivPRzJL/ZmbidOIEptjQXZycKLzq\nKgp8m1DQoSMFkZEUNPG1Pvf1pbBRI7hwxFGQD6tXV028dUyV9lUV8Pb2tv2FfeIfr5BbxVPcunbs\niM+Tfy5znfDwcL7//ntWr17NhAkT8PX1JSkpCW9vb8LDw8nKymL58uWsWrXKVvMpOzub//73v2Rl\nZXHLLbdQUFCAMYZBgwZx9uxZ9u/fz/Hjx20jiGHDhrFixQqysrL4448/yM/PJysri6+//pqlS5fy\n0ksvAZCTk0Nqairt27e3fS7du3fn8ccf56677mLIkCH4+fmViD87O5trr72WwMDAEqOVc//Pyckh\nLy+PrKwszp49i8ViIT09nXXr1jF8+HDb+mfPnr1otHMlOnPmTKW/I7X+DKuIHAVKrXomIgllbDcP\nmAfWWlvVUWepIvWbRISCEyeso4gDaReNKvIOHkSKpgI9x6VpUyz+/lgCArD4++Hq74/Fzx+Lvx+W\n5s31BHkF1MZaW+dqSp12tVDo7Fyl7VtcLZetWRUVFcXPP//M9u3b6dGjBydPnuStt97Cy8uLMWPG\n0KhRIywWC8888wzjxo0rse2cOXNwc3OzvYerq6vtuTHGtrxhw4Y4OTnRqFEjGjRogIuLi22dzz//\nnICAgBLtFq+1NXXqVIYNG8bKlSsZOHAgq1atKlGU0NPTE09PzxL7Wfy9jTFYLNbPwc3NDTc3Nxo2\nbIiPjw9bt26t5Cdbf7m7u192LphLccRvogNAy2LP/YuW1XmFp0+TWyxR5B04YH2eZn1eePp0ifWd\nvb2x+PvjFhCAZ79+55OFvz+WFi1wKmUor+qf5s8845D3veGGG5g9ezZt27bF2dmZxo0bc+LECVJS\nUnjnnXcA69VcU6ZM4Z577sHT05MDBw5gsVjo06cPMTExPP300+Tn57Ns2TLGjRuHj48PPj4+rF27\nll69erFo0aJS33vgwIG88cYbvPHGGxhj+Pnnny/6JbZnzx6Cg4MJDg5m06ZNbN++vdTqtsVdffXV\npKamEhAQwOeff35RMvXy8qJNmzZ88skn3HnnnYgIW7duvWyVXVU2RySSTUAHY0wbrAlkJDDaAXFU\nWGFurvUcxYGD5KWl4bl+HWlffGEbVRQcP15ifePhgau/HxY/fxqEh5dMFH5+ONeyKqfqyhIcHExm\nZiajR48usSw7OxtfX18AbrrpJlJTU4mMjASso4APPviA0NBQ7r77bkJCQmjWrBnh4eG2NhYsWMAD\nDzyAMcY2g+KFpkyZwqRJk+jSpQuFhYW0adPGNrf7OXPmzCEpKQknJycCAwO5+eabL7tPM2fOZPDg\nwTRt2pTu3buTnZ190TqLFi1i/PjxvPDCC+Tl5TFy5EhNJHaq1jLyxpjFWE+M+wKHgaki8p4x5hZg\nDuAMzBeRuGoLgsqXkc9NO0DmG6/bRhX5GRlQ7PMSZ2dc/f1w9StKDv7+1sRxLlE0bqwntmuJ2nho\nq76WkbdHbSwjf6WotWXkRWTUJZavBFZW53tXldMbN+Hq50fDyEjrVVDFksV/UlOJKmNKUKWUuhLo\n2doyuPr70SGpjMsCi107r5RSVyq9I00ppZRdNJEo5SBXwjTXqm6w92dRE4lSDuDu7s7Ro0c1mSiH\nExGOHj1aauWA8tJzJEo5gL+/P2lpaRw5csTRodQqZ86csesXmqocd3d3/P39K729JhKlHMBisdBG\nKy9fJDk5udJ3VyvH0UNbSiml7KKJRCmllF00kSillLJLtZZIqS2MMUeA/5Xykjdw0o7lvkBmKevV\nlEvFWVNtVWSby61b1uv29hPUn76qbDvl3c6efirr9dKW1+d+qmxbVdVPl1unvP3USkSaXjYaEbli\nH8A8e5YDP9bG+GuqrYpsc7l1y3rd3n6qT31V2XbKu509/VSJPqm3/VTZtqqqny63TkW/U5d7XOmH\ntpZV0XJHqcp4KtNWRba53LplvV7X+wmqLqbKtlPe7ezpp7JeL215fe6nyrZVVf10uXWq9LtzRRza\nqi7GmB+lHJUxleNpX9UN2k9105U+IrHXPEcHoMpN+6pu0H6qg3REopRSyi46IlFKKWUXTSRKKaXs\noolEKaWUXTSRVBNjzO3GmHeMMR8ZY25ydDzq0owxbY0x7xljEh0diyrJGNPQGLOw6Lt0j6PjUaXT\nRFIKY8x8Y0yGMebXC5YPMsbsMMbsNsZMLqsNEVkqIg8BjwB3V2e8V7Iq6qu9IvJg9Uaqzqlgnw0D\nEou+S0NqPFhVLppISpcADCq+wBjjDMQDNwOdgVHGmM7GmGBjzPILHs2Kbfpc0XaqeiRQdX2lakYC\n5ewzwB/YX7RaQQ3GqCpA5yMphYisMca0vmBxBLBbRPYCGGOWAENF5EVg8IVtGGMMMBP4SkR+qt6I\nr1xV0VeqZlWkz4A0rMlkC/qHb62lHVN+fpz/ywisP+B+Zaw/EegPjDDGPFKdgamLVKivjDFNjDFv\nA92MMU9Xd3CqVJfqs8+A4caYt6idJVUUOiKpNiLyOvC6o+NQlyciR7Gey1K1jIicBsY4Og5VNh2R\nlN8BoGWx5/5Fy1Tto31V92if1WGaSMpvE9DBGNPGGOMKjAS+dHBMqnTaV3WP9lkdpomkFMaYxcB6\nIMAYk2aMeVBE8oEJwCogFfhYRFIcGafSvqqLtM/qHy3aqJRSyi46IlFKKWUXTSRKKaXsoolEKaWU\nXTSRKKWUsosmEqWUUnbRRKKUUsoumkiUUkrZRROJUkopu2giUaoGFc2cKcaYTo6ORamqoolEqZo1\nCvix6F+l6gVNJErVEGOMJxAFjKVYIjHGJBljBhT9/wVjzBuOiVCpytH5SJSqOUOB70TkF2NMtjEm\nTEQ2A1OBvxVN+9sNnZtc1TE6IlGq5owCPi76/8dFzxGRNYAB/gyMFBGdm1zVKZpIlKoBxpjGQA/g\n66JFHwN3G6tg4BogV0SyHBWjUpWliUSpmjECWCkiZwFEZC+QDvQBFmE97JVtjBnkuBCVqhydj0Sp\nGmCMSQJCgFPFFvsBu4HHRORbY0wfYJaIRDoiRqUqSxOJUkopu+ihLaWUUnbRRKKUUsoumkiUUkrZ\nRROJUkopu2giUUopZRdNJEoppeyiiUQppZRdNJEopZSyy/8HrEznqQ1UnyUAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from math import pi, cos\n", "\n", "def trapezoid_rule(start_point, end_point, f, number_of_bins):\n", " \n", " bin_size = float(end_point - start_point)/number_of_bins\n", " \n", " running_total = 0.\n", " # Loop to create each trapezoid \n", " for i in range(number_of_bins): # note this function takes a slightly different approach to achieve the same thing as above\n", " # Set the start of this bin\n", " bin_start = start_point + (bin_size * i)\n", " # Find the area of the current trapezoid and add it to the running total\n", " running_total += bin_size*float( f(bin_start)+f(bin_start+bin_size) )/2.0\n", " \n", " return running_total\n", "\n", "\n", "def simpsons_rule(start_point, end_point, f, number_of_bins):\n", " \n", " bin_size = float(end_point - start_point)/number_of_bins\n", " \n", " running_total = 0.\n", " for i in range(number_of_bins):\n", " bin_start = start_point + bin_size * (i)\n", " bin_mid = bin_start + bin_size/2\n", " bin_end = bin_start + bin_size\n", " running_total += (bin_size/6)*float(f(bin_start)+4*f(bin_mid)+f(bin_end))\n", " \n", " return running_total\n", "\n", "\n", "def simpsons_composite_rule(start_point, end_point, f, number_of_bins=10):\n", " \n", " bin_size = float(end_point - start_point)/number_of_bins\n", "\n", " running_total = 0.\n", " for i in range(number_of_bins):\n", " bin_start = start_point + bin_size * (i)\n", " bin_quarter = bin_start + bin_size*1./4\n", " bin_mid = bin_start + bin_size*2./4\n", " bin_three_quarter = bin_start + bin_size*3./4\n", " bin_end = bin_start + bin_size\n", " running_total += (bin_size/12)*float(f(bin_start)+4*f(bin_quarter)\n", " +2*f(bin_mid)+4*f(bin_three_quarter)\n", " +f(bin_end))\n", " \n", " return running_total\n", "\n", "\n", "def weddles_rule(start_point, end_point, f, number_of_bins=10):\n", " result_simpson = simpsons_rule(start_point, end_point, f, number_of_bins)\n", " result_composite = simpsons_composite_rule(start_point, end_point, f, number_of_bins)\n", " \n", " return result_composite + float(result_composite - result_simpson)/15\n", "\n", "\n", "def f(x):\n", " return x*x*cos(x)\n", "\n", "int_exact = 4*pi\n", "start_point = 0.\n", "end_point = 2*pi\n", "\n", "dx = []\n", "err_tra = []\n", "err_sim = []\n", "err_sic = []\n", "err_wed = []\n", "\n", "number_of_bins = 1\n", "\n", "while number_of_bins < 1e3:\n", " \n", " int_tra = trapezoid_rule(start_point, end_point, f, number_of_bins)\n", " int_sim = simpsons_rule(start_point, end_point, f, number_of_bins)\n", " int_sic = simpsons_composite_rule(start_point, end_point, f, number_of_bins)\n", " int_wed = weddles_rule(start_point, end_point, f, number_of_bins)\n", " \n", " err_tra_cur = abs(int_tra - int_exact)\n", " err_sim_cur = abs(int_sim - int_exact)\n", " err_sic_cur = abs(int_sic - int_exact)\n", " err_wed_cur = abs(int_wed - int_exact)\n", " \n", " err_tra.append(err_tra_cur)\n", " err_sim.append(err_sim_cur)\n", " err_sic.append(err_sic_cur)\n", " err_wed.append(err_wed_cur)\n", " \n", " dx.append((end_point-start_point)/number_of_bins)\n", " \n", " number_of_bins *= 3\n", "\n", " \n", "pylab.loglog(dx, err_tra, label = \"Trapezoidal rule\")\n", "pylab.loglog(dx, err_sim, label = \"Simpson's rule\")\n", "pylab.loglog(dx, err_sic, label = \"Composite Simpson's rule\")\n", "pylab.loglog(dx, err_wed, label = \"Weddle's rule\")\n", "pylab.xlabel('$\\Delta x$');pylab.ylabel('Error');pylab.grid(True);pylab.legend(loc='best')\n", "pylab.show()\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "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.2" } }, "nbformat": 4, "nbformat_minor": 1 }