Skip to content

Instantly share code, notes, and snippets.

Show Gist options
  • Save adcroft/34d4989dd2dc30919f567810069c6035 to your computer and use it in GitHub Desktop.
Save adcroft/34d4989dd2dc30919f567810069c6035 to your computer and use it in GitHub Desktop.
Examples of quadrature and test of convergence
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def quad_average(y):\n",
" \"\"\"Returns the average value found by quadrature at order n.\n",
" y is a list of values in order from x=-1 to x=1.\"\"\"\n",
" if len(y)==2: # 1, 1\n",
" d = 1./2.\n",
" return d * ( y[0] + y[1] )\n",
" if len(y)==3: # 1/3, 4/3, 1/3\n",
" d = 1./6.\n",
" return d * ( 4. * y[1] + ( y[0] + y[2] ) )\n",
" if len(y)==4: # 1/6, 5/6, 5/6, 1/6\n",
" d = 1. / 12.\n",
" return d * ( 5. * ( y[1] + y[2] ) + ( y[0] + y[3] ) )\n",
" if len(y)==5: # 9/10, 49/90, 64/90, 49/90, 9/90\n",
" d = 1. / 180.\n",
" return d * ( 64.* y[2] + ( 49. * ( y[1] + y[3] ) ) + 9. * ( y[0] + y[4] ) )\n",
" raise Exception('Uncoded order')\n",
"def quad_positions(n=3):\n",
" \"\"\"Returns weights wa and wb so that the element [xa,xb] is sampled at positions\n",
" x=wa(xa+xb*xb).\"\"\"\n",
" if n==2:\n",
" return numpy.array([0.,1.]),numpy.array([1.,0.])\n",
" if n==3:\n",
" return numpy.array([0.,0.5,1.]),numpy.array([1.,0.5,0.])\n",
" if n==4:\n",
" r5 = 0.5 / numpy.sqrt(5.)\n",
" return numpy.array([0.,0.5-r5,0.5+r5,1.]),numpy.array([1.,0.5+r5,0.5-r5,0.])\n",
" if n==5:\n",
" r37 = 0.5 * numpy.sqrt(3./7.)\n",
" return numpy.array([0.,0.5-r37,0.5,0.5+r37,1.]),numpy.array([1.,0.5+r37,0.5,0.5-r37,0.])\n",
" raise Exception('Uncoded order')"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" Analytic = 1.6989760198247865\n",
"Quadrature = 1.7169000409973312\n",
"ratio = 1.0105498965043622\n",
"rel. error = 0.010549896504362164\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xlc1NX++PHXEVAUV8AdFEzcU0TcTXO5pmWamTfLSrH0dsvKltty65Z5637bt5+3a94ys0xLUzNLrpaa+wKKG24oKAgiggLKIjDv3x8zTqggoyIzwPvpYx7z+ZzP+cy85zC+5zNnzud8jIiglFKq8qji7ACUUkqVLU38SilVyWjiV0qpSkYTv1JKVTKa+JVSqpLRxK+UUpWMJn6llKpkSkz8xhhPY8xWY8xOY8xeY8zrRdSpZoz5zhgTY4zZYowJKLTtJVv5AWPMbaUbvlJKqavlyBF/LjBARDoBwcAQY0yPS+o8DJwWkZbAh8DbAMaYdsAYoD0wBPjUGONWWsErpZS6eu4lVRDrqb1nbasettulp/uOAKbalhcC040xxlY+X0RygVhjTAzQDdh0pef09fWVgIAAB1+CUkqpyMjIUyJS35G6JSZ+ANtReiTQEvi3iGy5pEpTIB5ARPKNMemAj618c6F6CbayKwoICCAiIsKR0JRSSgHGmKOO1nXox10RKRCRYMAP6GaM6XDpcxa12xXKL2OMmWSMiTDGRKSkpDgSllJKqWtwVaN6ROQMsAZrf31hCYA/gDHGHagDpBUut/EDEot57JkiEioiofXrO/RtRSml1DVwZFRPfWNMXdtydWAQsP+SakuBcbble4BVtt8GlgJjbKN+AoEgYGtpBa+UUurqOdLH3xj4ytbPXwX4XkSWGWOmAREishT4Avja9uNtGtaRPIjIXmPM90A0kA88LiIF1xJoXl4eCQkJ5OTkXMvuqoLz9PTEz88PDw8PZ4eilMszrjgff2hoqFz6425sbCy1atXCx8cH64AhpaxEhNTUVDIzMwkMDHR2OEo5hTEmUkRCHalbbs7czcnJ0aSvimSMwcfHR78NqvJp/UcQu/bisti11vIbpNwkfkCTviqWvjdUudU0BBaM/yP5x661rjcNuWFP6dA4fqWUUjdIYF8YPZvz8x8iusloOiX/gBk921p+g5SrI35n++STT2jbti1jx45lyZIlTJs27Yr1n3vuOVatWlVG0SmlyitL81tYaG4jOHYm0mXCDU36oEf8V+XTTz9l+fLlBAYG0qtXL5YuXXrF+k888QQTJ05kwIABZRShUqo8ilizhNuyf2Z/60dpEzkLWvTVI35X8Oijj3LkyBGGDx/O22+/TbVq1fD19QVgxIgRzJkzB4DPPvuMsWPHAtC8eXNSU1M5ceKE0+JWSrk2y+HfabXuSd6s8TxBY96C0bMv7vO/AcrlEf/rP+0lOjGjVB+zXZPavHZn+2K3z5gxg/DwcFavXs1PP/1ESMgfP7zMnDmT3r17ExgYyPvvv8/mzX9MTxQSEsKGDRsYNWpUqcarlKoYDuxYy7TcJxgzYhRuVYy9z5/j22/YUX+5TPzOlpSUROFpJRo2bMi0adPo378/ixcvxtvb276tQYMGJCYWOUuFUqqSK7AITx7ri/jCsI5N/tgQeGO7espl4r/SkXlZqF69Ounp6ReV7d69Gx8fn8uSfE5ODtWrVy/L8JRS5cSyXYkcOnmW6fd3th7tlxHt478Gbdu2JSYmxr6+detWli9fzo4dO3jvvfeIjY21bzt48CAdOlw6malSqrIrsAgf/3aI1g1rcXuHxmX63Jr4r0Hfvn3ZsWMHIkJubi4TJ05k1qxZNGnShPfff58JEyYgIuTl5RETE0NoqENnUSulKpGlO49zJOUcUwYFUaUMj/ahnHb1OEtcXJx9edCgQfz2228MGjSInTt32suHDx/O8OHDAVi2bBn33HMP7u7azEqpP+QXWPjktxjaNq7Nbe0blfnz6xH/Nfr73/9OVlbWFevk5+fz7LPPllFESqnyYklUIrGnnHO0D3rEf80aNmxoP7IvzujRo8soGqVUeZFXYOGT3w7RvkltBrdr6JQY9IhfKaXK0OLtxzmWlsXTg1o5bXJBTfxKKVVGcvIK+Pi3Q3T0q8PAtg2cFocmfqWUKiNzNsVx/Ew2Lwxp49SpxDXxK6VUGTiTdZ7pq2Lo16o+vVv6OjUWTfxKKVUG/r06hszcfF66vY2zQyk58Rtj/I0xq40x+4wxe40xTxVR52/GmCjbbY8xpsAY423bFmeM2W3bFnH5M9wATriUmVJKFSc+LYuvNh7lnhA/2jSq7exwHDrizweeFZG2QA/gcWNMu8IVRORdEQkWkWDgJeB3EUkrVKW/bXvZnMJ6Ay5l9o9//IOPP/7Yvv7yyy/z+uuv07dvX4KDg+nQoQPr1q27vriVUhXSeysOUKUKPDO4lbNDARxI/CKSJCLbbcuZwD6g6RV2uQ+YVzrhXaML05ouGA+r3rTeX+elzB5++GG++uorACwWC/Pnz6dmzZrcdtttREVFsXPnToKDg0sjeqVUBbI7IZ0foxJ5uE8gjeu4xoSNV3UClzEmAOgMbClmew1gCDC5ULEAK4wxAnwmIjOvKdKrFdgXQh+Gte9A3+eve4rTgIAAfHx82LFjB8nJyXTu3JmuXbsyYcIE8vLyuOuuuzTxK6UuIiL865d9eHtV5S/9bnJ2OHYO/7hrjKkJ/ABMEZHiroJyJ7Dhkm6e3iISAgzF2k1UZAY2xkwyxkQYYyJSUlIcDat4sWsh4gtr0o/4olSuZvPII48we/ZsvvzySyZMmEDfvn1Zu3YtTZs25cEHH7RfhUsppQDWHEhh05FUnhzQktqeHs4Ox86hxG+M8cCa9OeKyKIrVB3DJd08IpJouz8JLAa6FbWjiMwUkVARCS18kZNrcqFPf/RsGPByqV3KbOTIkYSHh7Nt2zZuu+02jh49SoMGDZg4cSIPP/ww27dvv764lVIVRoFF+L/l+wjwqcH93Zs7O5yLlNjVY6xnGXwB7BORD65Qrw7QD3igUJkXUEVEMm3Lg4Fp1x11SY5vv7hPv5QuZVa1alX69+9P3bp1cXNzY82aNbz77rt4eHhQs2ZNPeJXStktjIznYPJZPh0bQlV31xo570gff2/gQWC3MSbKVvZ3oBmAiMywlY0EVojIuUL7NgQW285Qcwe+FZHw0gj8ivpMubysFC5lZrFY2Lx5MwsWLABg3LhxjBs37roeUylV8ZzLzeeDlQfp3KwuQzuU/bTLJSkx8YvIeqDEc4tFZDYw+5KyI0Cna4zNpURHRzNs2DBGjhxJUFCQs8NRSrmw6atjSM7I5dOxIU6dmqE4Oi2zg9q1a8eRI0ecHYZSysUdSTnL5+uOMCrEjy7NvZ0dTpFcq+NJKaXKMRHh9Z+i8XR344WhrZ0dTrE08SulVCn5dd9Jfj+YwlODgmhQy9PZ4RRLE79SSpWCnLwCpi3bS1CDmozrFeDscK5I+/iVUqoUfPb7EeLTsvn2ke54uLn2MbVrR+eCFi9ejDGG/fv3X9fjjB8/noULF16xzr/+9a+L1nv16nVdz6mUujHi07L4dE0Md3RsTC8nz7XvCE38V2nevHn06dOH+fPn3/DnujTxb9y48YY/55Xk5+c79fmVclVv/ryPKsbw8u1tnR2KQzTxX4WzZ8+yYcMGvvjii4sS/5o1a7j11lu55557aNOmDWPHjkVEAJg2bRpdu3alQ4cOTJo0yV5+wW+//cbIkSPt6ytXruTuu+/mxRdfJDs7m+DgYMaOHQtAzZo17fXeeecdbr75Zjp16sSLL754Waw//fQT3bt3p3PnzgwaNIjk5GQsFgsBAQGcOXPGXq9ly5YkJyeTkpLCqFGj6Nq1K127dmXDhg0ATJ06lUmTJjF48GAeeugh4uLiuOWWWwgJCSEkJMT+YWSxWHjsscdo3749w4YN4/bbb7d/o4mMjKRfv3506dKF2267jaSkpOv6OyjlStYeTCF87wkmD2hJk7quMftmiUTE5W5dunSRS0VHR1+03q9fv8tu//73v0VE5Ny5c0Vu//LLL0VEJCUl5bJtjvj6669lwoQJIiLSs2dPiYyMFBGR1atXS+3atSU+Pl4KCgqkR48esm7dOhERSU1Nte//wAMPyNKlS0VEZNy4cbJgwQKxWCzSunVrOXnypIiI3HffffY6Xl5eFz3/hfVffvlFevbsKefOnbvsOS5IS0sTi8UiIiL//e9/5ZlnnhERkSeffFJmzZolIiKbN2+WgQMH2p/3QsxHjx6VNm3aiIjIa6+9JiEhIZKVlWVv2+zsbBEROXjwoFz4Wy1YsECGDh0qBQUFkpSUJHXr1pUFCxbI+fPnpWfPnvbXN3/+fAkLC3Oova/Wpe8RpW603LwC6f/eaun3zirJyct3aixAhDiYY/XH3aswb948pkyxTgcxZswY5s2bR0iI9eIu3bp1w8/PD4Dg4GDi4uLo06cPq1ev5p133iErK4u0tDTat2/PnXfeaX9MYwwPPvgg33zzDWFhYWzatKnEOX9+/fVXwsLCqFGjBgDe3pefJJKQkMC9995LUlIS58+fJzAwEIB7772XadOmERYWxvz587n33nvtjxkdHW3fPyMjg8zMTACGDx9O9erWI5m8vDwmT55MVFQUbm5uHDx4EID169czevRoqlSpQqNGjejfvz8ABw4cYM+ePfzpT38CoKCggMaNGzvc5kq5splrD3Mk5RyzxodSzd3N2eE4rNwm/jVr1hS7rUaNGlfc7uvre8XtRUlNTWXVqlXs2bMHYwwFBQUYY3jnnXcAqFatmr2um5sb+fn55OTk8NhjjxEREYG/vz9Tp04lJyfnsscOCwvjzjvvxNPTk9GjR+PufuU/i4iUeBr4E088wTPPPMPw4cNZs2YNU6dOBaBnz57ExMSQkpLCkiVLeOWVVwBrV82mTZvsCb4wLy8v+/KHH35Iw4YN2blzJxaLBU9PT3tMxcXavn17Nm3adMV4lSpvDqec5ZPfYrjj5sYMaNPQ2eFcFe3jd9DChQt56KGHOHr0KHFxccTHxxMYGMj69euL3edCkvf19eXs2bPFjuJp0qQJTZo04Y033mD8+PH2cg8PD/Ly8i6rP3jwYGbNmkVWVhYAaWlpl9VJT0+naVPrhdIuXDkMrN8wRo4cyTPPPEPbtm3x8fGxP+b06dPt9aKioihKeno6jRs3pkqVKnz99dcUFBQA0KdPH3744QcsFgvJycn2D9bWrVuTkpJiT/x5eXns3bu3yMdWqrywWISXFu3G06MKrw1vV/IOLkYTv4PmzZt30Y+wAKNGjeLbb78tdp+6desyceJEbr75Zu666y66du1abN2xY8fi7+9Pu3Z/vIkmTZpEx44d7T/uXjBkyBCGDx9OaGgowcHBvPfee5c93tSpUxk9ejS33HILvr4XDy+79957+eabb+zdPACffPIJERERdOzYkXbt2jFjxoxLHxKAxx57jK+++ooePXpw8OBB+7eBUaNG4efnR4cOHfjLX/5C9+7dqVOnDlWrVmXhwoW88MILdOrUieDgYKePTlLqes3fFs/W2DRevqOtS5+hWxxT3Fd0ZwoNDZWIiIiLyvbt20fbtuVjqNS1mDx5Mp07d+bhhx92dijX7OzZs9SsWZPU1FS6devGhg0baNSo7KakrejvEeUakjNyGPTB73RoUodvJ3Z3mdk3jTGRIhLqSN1y28dfkXTp0gUvLy/ef/99Z4dyXYYNG8aZM2c4f/48//jHP8o06StVVl77cS+5+Rb+dffNLpP0r5YmfhcQGRnp7BBKxdX+YK5UeRO+5wThe0/w/JDWBPp6lbyDi9I+fqWUckBGTh6v/riHto1rM/GWFs4O57po4ldKKQe8vXw/p87m8tbdN7v8JGwlKd/RK6VUGdgam8bcLccI6x1IJ/+6zg7nupWY+I0x/saY1caYfcaYvcaYp4qoc6sxJt0YE2W7vVpo2xBjzAFjTIwx5vJJZZRSyoWdy83nuQU78feuzjN/auXscEqFI0f8+cCzItIW6AE8bowp6oyFdSISbLtNAzDGuAH/BoYC7YD7itm3XEhISGDEiBEEBQXRokULJk+eTG5ubqk89po1axg2bNhV7RMXF3fF8wiUUtfvzV/2EX86i/dHB+NVrWKMhykx8YtIkohsty1nAvuApg4+fjcgRkSOiMh5YD4w4lqDdSYR4e677+auu+7i0KFDHDp0iOzsbJ5//vkb+rxXmgr5WhP/hbNtnUGndlblyeoDJ/l2yzEm3dKCboGueeH0a3FVffzGmACgM7CliM09jTE7jTHLjTHtbWVNgfhCdRJw/EPDpaxatQpPT0/CwsIA63w8H374IXPmzOHs2bPMnj2byZMn2+sPGzbMPrzxr3/9K6GhobRv357XXnvNXic8PJw2bdrQp08fFi1aZC93dCrkF198kXXr1hEcHMyHH354xRhq1qzJq6++Svfu3dm0aZNDUyXr1M6qMjt97jwvLNxFq4Y1ebqCdPHYOTqNJ1ATiATuLmJbbaCmbfl24JBteTTweaF6DwL/r5jHnwREABHNmjW7bMrRS6fcHb98vCw+tFhERM4XnJfxy8fL0hjrdMZZeVkyfvl4WX5kuYiIZORmyPjl42Vl3EoREUnLTpPxy8fL6mOrRUQkJSulxClPP/74Y5kyZcpl5cHBwbJjxw758ssv5fHHH7eX33HHHbJ6tfXxL0ybnJ+fL/369ZOdO3dKdna2+Pn5ycGDB8Viscjo0aPljjvuEBHHp0JevXq1fR8RuWIMgHz33XfW9nJwquTyNrWzTsusStPjcyPlppd+lt0JZ5wdikMo7WmZjTEewA/AXBFZdOl2EckotPyLMeZTY4wv1iN8/0JV/YDEYj6AZgIzwTplgyNxlSUpZkZMcWDKi++//56ZM2eSn59PUlIS0dHRWCwWAgMDCQoKAuCBBx5g5syZ9n0cmQr5ari5uTFq1CjA8amSdWpnVVkt3ZnIsl1JPDe4FR2a1nF2OKWuxMRvrNnuC2CfiHxQTJ1GQLKIiDGmG9YupFTgDBBkjAkEjgNjgPtLI/Avh3xpX/ao4nHRenX36het16pa66L1ep71Llr3rV7yNTLbt2/PDz/8cFFZRkYGycnJtG7dmj179mCxWOzbLszMGRsby3vvvce2bduoV68e48ePt2+70unejkyFfCl3d/ciYwDw9PTEzc06X7g4OFWyTu2sKqPkjBz+sWQPwf51ebTfTc4O54ZwpI+/N9YumgGFhmvebox51BjzqK3OPcAeY8xO4BNgjO3bRz4wGfgf1h+FvxeRcjkn78CBA8nKyrJfJKWgoIBnn32WyZMnU716dQICAoiKisJisRAfH8/WrVsB64eDl5cXderUITk5meXLlwPQpk0bYmNjOXz4MGCd/bM4xU2FXKtWLfsRNVBsDJdydKpkndpZVTYiwvMLd5GbX8AHf+6Eezk/Uas4jozqWS8iRkQ6yh/DNX8RkRkiMsNWZ7qItBeRTiLSQ0Q2Ftr/FxFpJSI3icibN/LF3EjGGBYvXszChQsJCgrCx8eHKlWq8PLLLwPQu3dvAgMDufnmm3nuuefsV+bq1KkTnTt3pn379kyYMIHevXsD1iPwmTNncscdd9CnTx+aN29e7HMXNxVyx44dcXd3p1OnTnz44YfFxnApR6dK1qmdVWUzd8sxfj+YwktD29Kifs2SdyindFrma7Rx40buu+8+Fi1aRJcuXZwdToVxPVM7u9p7RJUvB5MzufP/radboDdfhXWjSpXyNfOmTstcBnr16sXRo0edHUaFo1M7K2fIPl/A43O3U8vTnff/3KncJf2rpYlfuRSd2lk5w7Rlezl08ixfP9ytXF5R62qVq18uXLFbSrkGfW+oa/XTzkTmbY3nr7fexC1B9Z0dTpkoN4nf09OT1NRU/Q+uLiMipKamFjvMVaniHEvN4qVFuwlpVrfCTMDmiHLT1ePn50dCQgIpKSnODkW5IE9PT/z8/JwdhipHzudbmDxvO1UMfHJf53I/x/7VKDeJ38PDw37mqFJKXa93/7efXQnpzHggBL96NZwdTpmqPB9xSills2p/Mv9dF8tDPZszpEPlm/5DE79SqlKJT8vi2e930qZRLf5+e+U870MTv1Kq0sjJK+CvcyPJtwj/eaALnh5uzg7JKcpNH79SSl0PEeGVJXvYczyDL8aFEujrVfJOFZQe8SulKoW5W46xMDKBJwcGMbBtQ2eH41Sa+JVSFV7k0dO8/tNe+reuz5SBQc4Ox+k08SulKrSTmTk8NjeSxnWq89G9nSv8PDyO0D5+pVSFlVdgYfLcHaRn57H4sW7UqeHh7JBcgiZ+pVSF9a9f9rE1Lo2PxwTTtnFtZ4fjMrSrRylVIX237RhfbogjrHcAI4KbOjscl6KJXylV4Ww6nMrLi/dwS5AvL1fSk7SuRBO/UqpCiTt1jr/OjaS5Tw2m3x9SYa+bez1KbBFjjL8xZrUxZp8xZq8x5qki6ow1xuyy3TYaYzoV2hZnjNltu0h7xKX7KqVUaUnPymPCV9swwKzxXalTXX/MLYojP+7mA8+KyHZjTC0g0hizUkSiC9WJBfqJyGljzFBgJtC90Pb+InKq9MJWSqmL5RVYeOzbSOLTspj7SA+a+1TeM3NLUmLiF5EkIMm2nGmM2Qc0BaIL1dlYaJfNgE6MrpQqMyLCa0v3siEmlXfv6Ui3QG9nh+TSrqrzyxgTAHQGtlyh2sPA8kLrAqwwxkQaYyZdbYBKKVWSLzfE8e2WYzza7yZGh/o7OxyX5/A4fmNMTeAHYIqIZBRTpz/WxN+nUHFvEUk0xjQAVhpj9ovI2iL2nQRMAmjWrNlVvASlVGX2a3Qyb/wczeB2DXn+ttbODqdccOiI3xjjgTXpzxWRRcXU6Qh8DowQkdQL5SKSaLs/CSwGuhW1v4jMFJFQEQmtX79yXPBYKXV9Io+m8fi327m5aR0+GhOs0zE4yJFRPQb4AtgnIh8UU6cZsAh4UEQOFir3sv0gjDHGCxgM7CmNwJVSlVvMyUwe/iqCxnU8mTW+KzWq6kQEjnKkpXoDDwK7jTFRtrK/A80ARGQG8CrgA3xq/ZwgX0RCgYbAYluZO/CtiISX6itQSlU6J9JzGDdrG+5VqjBnQnd8alZzdkjliiOjetYDV/z+JCKPAI8UUX4E6HT5HkopdW3Ss/MY/+VWzmSd57u/9KSZT+W6UHpp0O9GSqlyIyevgElzIjiccpZZ47vSoWkdZ4dULmniV0qVCwUW4Znvo9gSa51t85YgHQRyrXQSC6WUyxMRXv1xD7/sPsErd7TV2TavkyZ+pZRLExH+9cs+5tpO0HrklhbODqnc08SvlHJpH/56iP+ui2V8rwBeGKInaJUGTfxKKZc14/fDfPLbIf4c6serw9phGxqurpMmfqWUS5qzKY63lu9neKcm/N/dHfWs3FKkiV8p5XK+j4jn1R/38qd2DXn/z51w06RfqjTxK6VcytKdibz4wy76tqrP9Ps746FX0Cp1Oo5fKeUyfow6ztPfRdE1wJvPHuhCNXc3Z4dUIelHqVLKJSzansDT30XRLdCbL8O6Ur2qJv0bRY/4lVJOtyAinud/2EWvm3z4/CFN+jeaJn6llFPN33qMlxbvpk9LX/77UCieHpr0bzTt6lFKOc3cLUd5cdFu+gbV16RfhvSIXynlFHM2xfHqj3sZ0KYB/3kgRH/ILUOa+JVSZUpE+HTNYd793wH+1K4h0+/vrEm/jGniV0qVGRHhzZ/38fn6WO4KbsK7ozvpOH0n0MSvlCoT+QUWXlq0mwWRCYzvFcCrw9rpNAxOoolfKXXD5eQV8OS8HayITmbKoCCeGhikE645UYnfsYwx/saY1caYfcaYvcaYp4qoY4wxnxhjYowxu4wxIYW2jTPGHLLdxpX2C1BKubazuflMmL2NFdHJTL2zHVMGtdKk72SOHPHnA8+KyHZjTC0g0hizUkSiC9UZCgTZbt2B/wDdjTHewGtAKCC2fZeKyOlSfRVKKZd06mwuD8/exp7EDD68txMjO/s5OySFA0f8IpIkIttty5nAPuDS656NAOaI1WagrjGmMXAbsFJE0mzJfiUwpFRfgVLKJR1JOcvdn27kQHImnz3QRZO+C7mqPn5jTADQGdhyyaamQHyh9QRbWXHlRT32JGASQLNmza4mLKWUi4mIS+ORORG4GcO8iT3o3Kyes0NShTg8jsoYUxP4AZgiIhmXbi5iF7lC+eWFIjNFJFREQuvXr+9oWEopF/PL7iTu/3wL9WpUZdFjvTTpuyCHEr8xxgNr0p8rIouKqJIA+Bda9wMSr1CulKpgRITP1x3h8W+3c3PTOvzw11409/FydliqCI6M6jHAF8A+EfmgmGpLgYdso3t6AOkikgT8DxhsjKlnjKkHDLaVKaUqkAKL8PpP0bzx8z6GtG/E3Ee64+1V1dlhqWI40sffG3gQ2G2MibKV/R1oBiAiM4BfgNuBGCALCLNtSzPG/BPYZttvmoiklV74Silny8zJ48l5O1h9IIWH+wTy8u1t9cQsF1di4heR9RTdV1+4jgCPF7NtFjDrmqJTSrm0o6nneOSrCI6cOsc/R7TnwZ4Bzg5JOUDP3FVKXZONh0/x2NztiMDXE7rRq6Wvs0NSDtLEr5S6anO3HOW1H/cS4OvF5w+FEuCrP+KWJ5r4lVIOyyuw8M9l0czZdJT+revz8X2dqe3p4eyw1FXSxK+UckhKZi5PzNvO5iNpTOrbgheGtMFNf8QtlzTxK6VKtP3YaR77Zjuns87z/uhOjOqi0y+UZ5r4lVLFEhG+2XyUacuiaVynOose60X7JnWcHZa6Tpr4lVJFyskr4O+Ld7No+3H6t67PR/d2pk4N7c+vCDTxK6Uucyw1i0e/iWTfiQymDAriyQFBelJWBaKJXyl1keW7k3j+h11UMYZZ47vSv3UDZ4ekSpkmfqUUYO3aefPnfXy9+Sid/Osy/b7O+HvXcHZY6gbQxK+U4kjKWR7/dgf7kjKY1LcFzw1uTVV3h2dtV+WMJn6lKrnFOxJ4efEeqrlXYdb4UAa0aejskNQNpolfqUrqXG4+U5fuZUFkAt3hbG2dAAAWvUlEQVQCvPn4vmAa16nu7LBUGdDEr1QltP3YaZ7+LopjaVk8MaAlTw0Mwt1Nu3YqC038SlUieQUWpq+KYfrqGBrV9uS7ST3pFujt7LBUGdPEr1QlEXvqHE9/F0VU/BnuDmnK1OHtdYK1SkoTv1IVnIjw3bZ4pi2LxsOtCtPv78ywjk2cHZZyIk38SlVgJ9JzeGnRLlYfSKF3Sx/eG91Jf8BVmviVqohEhAWRCfxzWTT5BcJrd7ZjXM8AnXZBAQ4kfmPMLGAYcFJEOhSx/W/A2EKP1xaob7vQehyQCRQA+SISWlqBK6WKVvgov1uAN++O7khzH71ClvqDI0f8s4HpwJyiNorIu8C7AMaYO4GnRSStUJX+InLqOuNUSpVARFgYmcC0ZdHkFVj0KF8Vq8TELyJrjTEBDj7efcC86wlIKXX14tOyeHnJHtYetB7lv3NPR70OripWqfXxG2NqAEOAyYWKBVhhjBHgMxGZeYX9JwGTAJo1a1ZaYSlVoeUXWJi1IZYPVh7EzRim3tmOh/QoX5WgNH/cvRPYcEk3T28RSTTGNABWGmP2i8jaona2fSjMBAgNDZVSjEupCml3QjovLtrF3sQMBrVtyLQR7WlSV0fsqJKVZuIfwyXdPCKSaLs/aYxZDHQDikz8SinHZJ3P54MVB5m1IRbfmtX4z9gQhnRohDF6lK8cUyqTcxhj6gD9gB8LlXkZY2pdWAYGA3tK4/muJCw8jCUxSwDIs+QRFh7GT4d/AiA7P5uw8DDCY8MByDyfSVh4GL8e/RWA0zmnCQsPY038GgBOZZ8iLDyM9cfXA3Di3AnCwsPYlLgJgPjMeMLCw9h2YhsAsemxhIWHEXUyCoBDpw8RFh7GnlPWl70/bT9h4WHsT9sPwJ5TewgLD+PQ6UMARJ2MIiw8jNj0WAC2ndhGWHgY8ZnxAGxK3ERYeBgnzp0AYP3x9YSFh3Eq2/rb+Zr4NYSFh3E65zQAvx79lbDwMDLPZwIQHhtOWHgY2fnZAPx0+CfCwsPIs+QBsCRmCWHhYfa2XHhwIY+seMS+Pn//fB799VH7+jfR3/DEb0/Y12fvmc3Tq5+2r3+++3P+9vvf7Oszds7gxXUv2ten75jOK+tfsa9/FPkRUzdOta+/t+093tj8hn397a1v8/bWt+3rb2x+g/e2vWdfn7pxKh9FfmRff2X9K0zfMd2+/uK6F5mxc4Z9/W+//43Pd39uX3969dPM3jPbvv7Eb0/wTfQ39vVHf32U+fvn29cfWfEICw8utK+XxXvvk40/M+j93/liSwTN2n/FOw/UZOjNjYnLiNP3XgV6791ojgznnAfcCvgaYxKA1wAPABG58EpGAitE5FyhXRsCi21HIe7AtyISXnqhK1V5HEvLYv+JTNZuPURQ7RBeuSuEBUdX41XNzdmhqXLIkVE99zlQZzbWYZ+Fy44Ana41sGsV+1YsH9n+XRD/53h4DOS8EPtWLG/Z/l2QMD4BxkPB2QJi34plKlMveszjfz0O90JeWh6xb8XyEi9dtP3EsyfgTjh/4jyxb8UyhSkXb3/lBB0GdSDnWA6xb8XyKI9etD3lXykE9QoiKyaL2LdiCSPsou2pH6XiH+zPuehzxL4Vy5i3xly8PTAV39a+ZEZlEvt+LCPfGnnR9jPtz1DLvxbpW9OJ/U8sQ98aetH29NB0fH19ObP+DLGzY7n1rVsv2p7VJ4saNWqQtiqN/d/v59Y3Lt7OQOvdqfBT7Fi2g1tf/2N79erVrd8FgZSlKWz+bTO38sd2Hx8f6GNdTl6YzKZNmy7a7ufnBz2sy0nfJhEVFcVyltu3t2rVCrpalxNnJ7Lm4BqWsMS+PTg4GDpblxM+S2Bzwmbm88dRu3tPd/g/6/Kx6cfYkbqD2YXeyg0GNoB21uWj7x9lRvYMZvDHkVvcsDh4zrp8I957sRPj2VL/ANN/2sLJnw7iXy+B7NpLeMf2EvW9Z72rKO+98WvGUxb0zF2lXNTpc+d54+docpvVZmCb+uyPqouHTp2sSoERcb0BNKGhoRIREeHsMJRyigMnMvnnsmjWx5wiqEFNpo3oQM+bfJwdlnJxxphIR2dH0CN+pVxE2rnzfLDyAN9uOUYtTw9eu7MdD/Rorkf5qtRp4lfKyfIKLMzZdJSPfz3IufMFPNijOVMGtaKeV1Vnh6YqKE38SjmJiLAyOpm3wvdzJOUctwT58o9h7WjVsJazQ1MVnCZ+pZxg+7HT/N8v+9gWd5oW9b34/KFQBrZtoCdhqTKhiV+pMhR76hzv/m8/v+w+gW/Narw5sgP3hvrrhc5VmdLEr1QZSMnMZfqqQ8zdcoyq7lWYMiiIibe0wKua/hdUZU/fdUrdQOlZeXy29jBfbojjfIGFMV39eWpQEA1qeTo7NFWJaeJX6gY4l5vP7I1xzPj9MJk5+dzZqQlPDwqiRf2azg5NKU38SpWmnLwCvt1yjE/XxHDq7HkGtW3AM39qTbsmtZ0dmlJ2mviVKgU5eQXM33qM//x+mOSMXHrd5MPMh1oT0qyes0NT6jKa+JW6Djl5BczbeowZtoTfLdCbD/8cTK+Wvs4OTaliaeJXqiTrP4KmIRDY1150/tAadm5ZxeNH+3IyM5fugd58dG9nnVNHlQua+JUqSdMQWDAeRs8ms3FPVi3/gb47/8YH55+gRYAXH4/RhK/KF038SpUksC/pw/6Lx9wHmZs/kNGyghkNX2XKkFF0b6EJX5U/mviVuoLEM9nMXHuE+dvO85j050n3H0gOeYqXRkxydmhKXTNN/EoVYW9iOp+vi+WnnYkAPNcqmYknfoduz9Mw4guIHXRRn79S5YkmfqVsRIR1h04xc+0R1secwquqG+N6BfBo8+PUX/5PuPcra7IPvMXe56/JX5VHjlxsfRYwDDgpIh2K2H4r8CMQaytaJCLTbNuGAB8DbsDnIvLWpfsr5Wy5+QUs25nEf9cdYf+JTBrUqsYLQ9pwf/dm1KnuAetXXJzkA/ta149v18SvyqUSL71ojOkLnAXmXCHxPyciwy4pdwMOAn8CEoBtwH0iEl1SUHrpRVUWUjJzmbvlKN9sPsaps7m0bliLiX1bMLxTE6q662yZqnwp1UsvishaY0zANcTRDYgRkSO2oOYDI4ASE79SN9Ke4+nM2hDLsp1JnC+w0L91fcJ6B3JLkK/Oh68qhdLq4+9pjNkJJGI9+t8LNAXiC9VJALoX9wDGmEnAJIBmzZqVUlhKWZ3Pt7Ai+gRfbYxjW9xpalR1475u/ozrFaATp6lKpzQS/3aguYicNcbcDiwBgoCiDp2K7VcSkZnATLB29ZRCXEqRlJ7NvC3HmLctnpTMXPy9q/PKHW35c1d/ant6ODs8pZziuhO/iGQUWv7FGPOpMcYX6xG+f6Gqfli/ESh1Q4kIGw+n8vWmo6zcl4xFhAGtG/BAz+b0C6pPlSranaMqt+tO/MaYRkCyiIgxphtQBUgFzgBBxphA4DgwBrj/ep9PqeKcOpvLD5EJzN8WT+ypc3h7VWXiLS0Y270Z/t41nB2eUi7DkeGc84BbAV9jTALwGuABICIzgHuAvxpj8oFsYIxYhwrlG2MmA//DOpxzlq3vX6lSY7EIGw6fYv7WeFZEnyCvQOgW4M2TA1sytENjPD3cnB2iUi6nxOGczqDDOVVJktKzWbT9OPO3HSM+LZu6NTwYFeLHfd38admglrPDU6rMlepwTqVcRU5eAb/uS+b7iATWH0rBItCjhTfPDW7Nbe0b6dG9Ug7SxK9cmoiwNzGDBRHxLIlKJD07jyZ1PHm8f0vu6eJHcx8vZ4eoVLmjiV+5pMQz2SyJOs7i7cc5dPIsVd2rcFv7Rvw51I9eN/nipiNzlLpmmviVyzibm8/y3Uks3nGcTUdSEYEuzevx5sgODLu5CXVq6Lh7pUqDJn7lVLn5Bfx+IIWlOxP5dV8yOXkWmnnX4KmBQYzs3FS7cpS6ATTxqzJXYBE2H0llaVQiy/ckkZGTj7dXVe7p4sfIzk0JaVZP58xR6gbSxK/KhMUiRBw9zS+7k/h5dxIpmbl4VXXjtvaNuDO4CX1a+uLhpjNiKlUWNPGrG8ZiEbYfO82yXUks35NEckYuVd2rcGur+owIbsqANg2oXlWHYCpV1jTxq1JVYBG2xaURvucE4XtOcCIjx57s7+jYmIFtG1Kzmr7tlHIm/R+ortv5fAsbD5/if3tPsGJvMqnnzlPNvQp9W9XnpY5tGNCmAbV0JkylXIYmfnVNMnPy+P1gCiujk1m1/ySZOfl4VXVjQNuGDO3QiH6t6uOlR/ZKuST9n6kclpSeza/RyayITmbzkVTyCgRvr6oMad+IoTc3otdNvjptglLlgCZ+VSyLRdh9PJ1V+0+yav9Jdh9PByDQ14uw3oH8qV1DQprV07NolSpnNPGri5zNzWf9oRRbsk/h1NlcjIHO/nV5fkhrBrdryE31a+o4e6XKMU38lZyIcDD5LGsOnOT3gylsi0sjr0Co5elOv1b1Gdi2Af1aNcDbq6qzQ1VKlRJN/JVQRk4eG2NOseZACr8fTCEpPQeANo1qMaF3IP3bNKBL83p6QpVSFZQm/kogv8DCzoQzrD14ivUxp4iKP0OBRahVzZ3eLX15amB9+rWuT+M61Z0dqlKqDGjir4BEhMMp59h0+BTrDp1i0+FUMnPzMQY6+tXlsVtvok9LX0L0qF6pSkkTfwWRlJ7NhphUNsacYuPhVE5kWLtvmtatzrBOjbklqD69bvKhbg3tq1eqsnPkYuuzgGHASRHpUMT2scALttWzwF9FZKdtWxyQCRQA+Y5eD1KV7GRGDpuOpLL5SBqbj6QSe+ocAN5eVel5kw+9b/Kld0sfmnnX0BE4SqmLOHLEPxuYDswpZnss0E9EThtjhgIzge6FtvcXkVPXFaUiOSOHLbHWJL/5SCpHUqyJvlY1d7oFejO2ezN63eRLm0a1qKLj6pVSV1Bi4heRtcaYgCts31hodTPgd/1hVW4iwrG0LLbEprE1No1tcWkcTc0CoKYt0Y/p6k/PFr60a1JbT6BSSl2V0u7jfxhYXmhdgBXGGAE+E5GZxe1ojJkETAJo1qxZKYfl2vILLOxLyiTiaBoRcaeJOJpGckYuAPVqeNA1wJsHezSnW6A37RrXxl1/kFVKXYdSS/zGmP5YE3+fQsW9RSTRGNMAWGmM2S8ia4va3/ahMBMgNDRUSisuV5SenUdU/Bkij54m8mgaO46dIet8AWD9MbZ7oA/dAr3pFuhNy/o1tetGKVWqSiXxG2M6Ap8DQ0Uk9UK5iCTa7k8aYxYD3YAiE39FZbEIMSln2X70NNuPnWb7sTPEnDwLQBUDbRrV5p4ufoQGeBPavB5N6upYeqXUjXXdid8Y0wxYBDwoIgcLlXsBVUQk07Y8GJh2vc/n6k5m5hB17Aw7E84QFX+GXfHpZObmA1C3hged/esyolMTOjerRyf/OjpPvVKqzDkynHMecCvga4xJAF4DPABEZAbwKuADfGobNnhh2GZDYLGtzB34VkTCb8BrcJqMnDz2HE9nV0I6uxLOsDM+neNnsgFwq2Jo06gWw4ObEOxfl5Dm9Wjh66VDK5VSTmdEXK87PTQ0VCIiIpwdxkXO5uYTnZhhS/Rn2HU83T6kEsDfuzod/erS2b8uwf51ad+kjl5PVilVZowxkY6eK6Vn7hYhPSuPvUnp7D2ewZ7EdHYfTyf21DkufEY2qu3JzX51uLtzU272q0vHpnWop7NXKqXKiUqd+EWEpPQc9iZmEJ2Ywd7EdKKTMkg4nW2v06SOJ+2b1uGu4KZ0aFqbDk3q0KC2pxOjVkqp61MxEv/6j6BpCAT2/aMsdi0c3w59pgCQk1fAweRM9iVlsC/Jer//RCbp2XkAGGO9slSwf13Gdm9Ouya1ad+kNr41qznjFSml1A1TMRJ/0xBYMB5Gz6ag+S0k71yB7/K/sLTVv1g1N5L9JzKJO3UOi62rpkZVN1o3qsUdHRvTtlEt2jWpQ5tGtfTi4EqpSqFiZLrAvuSP+pJzc8byTcEgxpiVjMt7ks2RdWjunUHrRrUY1rEJ7RrXok2j2jTzrqEnRSmlKq2KkfgB95v6sdn7Lh5PncPeoEd5qd8kWjaoSY2qFeYlKqVUqag4k77EruW27J+h7/O0P76Ajnm7NOkrpVQRKkbij11r7+NnwMvW+wXjreVKKaUuUjES//Ht1mR/YVRPYF/r+vHtzoxKKaVcUsXoC7EN2bxIYN+Lh3cqpZQCKsoRv1JKKYdp4ldKqUpGE79SSlUymviVUqqS0cSvlFKVjEvOx2+MSQGOXuPuvsCpUgyntGhcV0fjujoa19WpiHE1F5H6jlR0ycR/PYwxEY5ejKAsaVxXR+O6OhrX1anscWlXj1JKVTKa+JVSqpKpiIl/prMDKIbGdXU0rqujcV2dSh1XhevjV0opdWUV8YhfKaXUFZTbxG+M8TTGbDXG7DTG7DXGvG4rDzTGbDHGHDLGfGeMqeoicc02xsQaY6Jst+CyjKtQfG7GmB3GmGW2dae21xXicnp7GWPijDG7bc8fYSvzNsastLXXSmNMPReJa6ox5nih9rq9rOOyxVHXGLPQGLPfGLPPGNPT2W1WTExOby9jTOtCzx9ljMkwxkwpi/Yqt4kfyAUGiEgnIBgYYozpAbwNfCgiQcBp4GEXiQvgbyISbLtFlXFcFzwF7Cu07uz2uuDSuMA12qu/7fkvDLF7EfjN1l6/2dZdIS6w/h0vtNcvTorrYyBcRNoAnbD+TZ3dZkXFBE5uLxE5cOH5gS5AFrCYMmivcpv4xeqsbdXDdhNgALDQVv4VcJeLxOV0xhg/4A7gc9u6wcntVVRcLm4E1nYCJ7WXqzLG1Ab6Al8AiMh5ETmDE9vsCjG5moHAYRE5Shm0V7lN/GDvHogCTgIrgcPAGRHJt1VJAJo6Oy4R2WLb9KYxZpcx5kNjTLWyjgv4CHgesNjWfXCB9ioirguc3V4CrDDGRBpjJtnKGopIEoDtvoGLxAUw2dZes5zRBQW0AFKAL23ddp8bY7xwbpsVFxM4v70KGwPMsy3f8PYq14lfRApsX5P8gG5A26KqlW1Ul8dljOkAvAS0AboC3sALZRmTMWYYcFJEIgsXF1G1TNurmLjAye1l01tEQoChwOPGGFe5sk9Rcf0HuAlr92IS8L4T4nIHQoD/iEhn4BzO6wq7oLiYXKG9ALD9rjYcWFBWz1muE/8Ftq9ua4AeQF1jzIUri/kBiS4Q1xARSbJ1A+UCX2L9oCpLvYHhxpg4YD7WLp6PcH57XRaXMeYbF2gvRCTRdn8Sa99rNyDZGNMYwHZ/0hXiEpFk2wGHBfgvTmgvrN8YEwp9w12INek6s82KjMlF2uuCocB2EUm2rd/w9iq3id8YU98YU9e2XB0YhPVHm9XAPbZq44AfXSCu/YX+kAZrn92esoxLRF4SET8RCcD6tXKViIzFye1VTFwPOLu9jDFexphaF5aBwbYYlmJtJ3DO+6vIuC60l81Iyri9AETkBBBvjGltKxoIROPENisuJldor0Lu449uHiiL9hKRcnkDOgI7gF1Y/2iv2spbAFuBGKxfnaq5SFyrgN22sm+Amk5su1uBZa7QXleIy6ntZWuXnbbbXuBlW7kP1pEWh2z33i4S19e29tqFNXE0dtLfMBiIsMWxBKjnAm1WVEyu0l41gFSgTqGyG95eeuauUkpVMuW2q0cppdS10cSvlFKVjCZ+pZSqZDTxK6VUJaOJXymlKhlN/EopVclo4ldKqUpGE79SSlUy/x8m8zRddmMyJAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"def f(x):\n",
" \"\"\"Example (integrable) function\"\"\"\n",
"# return 2.*x + 10*x**9\n",
"# return 2.*x + 3*x**2\n",
"# return x\n",
"# return numpy.sin(x*numpy.pi/180.)\n",
"# return numpy.cos(x*numpy.pi/180.)\n",
" return 1. / numpy.cos(x*numpy.pi/180.) # secant\n",
"def F(x): # Integral of f(x)\n",
" \"\"\"Analytic integral of f(x)\"\"\"\n",
"# return x**2 + x**10\n",
"# return x**2 + x**3 \n",
"# return 0.5*x**2\n",
"# return -numpy.cos(x*numpy.pi/180.)\n",
"# return numpy.sin(x*numpy.pi/180.)*180./numpy.pi\n",
" return numpy.log( 1./numpy.cos(x*numpy.pi/180.) + numpy.tan(x*numpy.pi/180.) )*(180./numpy.pi)\n",
"x0,dx = 30, 40\n",
"x = numpy.array([x0, x0+dx])\n",
"xp = numpy.linspace(x[0],x[1],50) # x for plotting function\n",
"plt.plot( xp, f(xp), '-', label='f(x)');\n",
"\n",
"a,b = quad_positions(n=3)\n",
"xs = a*x[0] + b*x[1] # Sample position to evalute f(x)\n",
"ys = f(xs) # Values to pass to quad_average()\n",
"plt.plot( xs, ys, 'x', label='ys');\n",
"\n",
"truth = ( F(x[1]) - F(x[0]) ) / dx\n",
"print(' Analytic =',truth)\n",
"print('Quadrature =', quad_average(ys))\n",
"print('ratio =',quad_average( ys ) / truth )\n",
"print('rel. error =', ( quad_average( ys ) - truth ) / truth )\n",
"plt.plot(xp, truth+0*xp, 'k--', label='Analytic average');\n",
"plt.plot(xp, quad_average( ys )+0*xp, ':', label='Quadrature average');\n",
"\n",
"plt.legend();"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.010549896504362164"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def relerror(x,n=3):\n",
" a,b = quad_positions(n=n)\n",
" xs = a*[x0] + b*x[1]\n",
" ys = f(xs)\n",
" truth = ( F(x[1]) - F(x[0]) ) / ( x[1] - x[0] )\n",
" ival = quad_average(ys)\n",
" return numpy.abs( ( ival - truth ) / truth )\n",
"relerror(x)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEcCAYAAADHiMP9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsnXlYlVXXh+/NJIigCDgjioqzqTjliAMerRyy3jL9UnN61cwsZ3PKIUvtzcrMypzKRi1Ls8AhxVnQtFQUUVFQGWUeBM7Z3x8PnFCZOQjovq/rXPE8Zw/rCD2/s/faay0hpUShUCgUCrPSNkChUCgUZQMlCAqFQqEAlCAoFAqFIhMlCAqFQqEAlCAoFAqFIhMlCAqFQqEAlCAoHgOEEJ5CiNBi9F8nhJhvSpsUirKIEgRFuUAIESyESBFCJAohwoQQm4QQlUpgnlFCiMPZ70kpJ0gpl5h6LoWirKEEQVGeGCClrAS0BtoAc0rZnoeGEMKiIPcKO4ZCkR0lCIpyh5QyDPBGEwYAhBAVhBCrhBA3hBDhmds8Njn1F0LMFkJcEUIkCCEuCCGezbzfFFgHPJm5EonNvL9JCLE08+cAIcQz2cayEEJECSHaZl53EkIcFULECiHOCiE8c/scQohaQojtQohIIcQ1IcSUbO8tEkJsE0J8LYSIB0blcq+CEGK1EOJW5mu1EKJC5hieQohQIcQsIUQYsLFo/+KKxwUlCIpyhxCiDtAfCMp2+z3AHU0kGgK1gQW5DHEF6AZUBt4GvhZC1JRSBgATgGNSykpSyio59P0WeCnbtQ6IklKeFkLUBn4DlgJVgenAdiGEcw6fwQzYCZzNtLU3MFUIocvWbBCwDagCbM3l3ltAp8zP/QTQAZiXbYwamba4AuNz+fdQKAAlCIryxQ4hRAIQAkQACwGEEAIYB7whpbwjpUwA3gGG5jSIlPJHKeUtKaVBSvk9cBntQVoQvgEGCiEqZl4Py7wH8H/Abinl7syx9wD+wFM5jNMecJZSLpZSpkkprwJf3GfzMSnljsyxUnK5NxxYLKWMkFJGogncy9nGMAALpZR3s42hUOSI2lNUlCcGSyn3CiF6oD2EnYBYwBmoCJzStAEAAZjnNIgQYgTwJlAv81alzLHyRUoZJIQIAAYIIXYCA9H8GaB9C/+PEGJAti6WwJ85DOUK1MralsrEHDiU7Tokh37336sFXM92fT3zXhaRUsrU3D6PQpEdJQiKcoeU8qAQYhOwChgMRAEpQHMp5c28+gohXNG+ifdG+7atF0KcQRMQgIKk/83aNjIDLkgps7auQoCvpJTjCjBGCHBNStkojzY52XL/vVto4nI+87pu5r28xlAockRtGSnKK6sBLyFEaymlAe0h/4EQohqAEKL2ffvxWdiiPSQjM9u9ArTI9n44UEcIYZXH3N8BfYGJ/LtdBPA12spBJ4QwF0JYZzp26+QwxkkgPtPha5PZvoUQon2BPv2/fAvME0I4CyGc0PwmXxdyDIUCUIKgKKdk7pdvAbICxmahOZmPZ57A2Qs0zqHfBeB94Bjaw78lcCRbk/1o37bDhBBRucx9O7N/Z+D7bPdD0Jy+c9EEJwSYQQ7/n0kp9cAANGfwNbRVzno0R3dhWIrmp/gb+Ac4nXlPoSg0QhXIUSgUCgWoFYJCoVAoMlGCoFAoFApACYJCoVAoMlGCoFAoFApACYJCoVAoMilXgWlOTk6yXr16pW2GQqFQlCtOnToVJaV8IKfW/ZQrQahXrx7+/v6lbYZCoVCUK4QQ1/NvpbaMFAqFQpGJEgSFQqFQAEoQFAqFQpFJufIh5ER6ejqhoaGkpj6aGX6tra2pU6cOlpaWpW2KQqF4xCn3ghAaGoqdnR316tUjWy78RwIpJdHR0YSGhlK/fv3SNkehUDzilPsto9TUVBwdHR85MQAQQuDo6PjIrn4UCkXZotwLAvBIikEWj/JnUygUZYtSFQQhhJsQ4kshxLbStKO02LNnDx4eHrRs2RIPDw/2799f2iYpFIrHmCILghBigxAiQghx7r77/YQQl4QQQUKI2XmNIaW8KqUcU1QbyjtOTk7s3LmTf/75h82bN/Pyyy/n30mhUChKiOKsEDYB/bLfEEKYA58A/YFmwEtCiGZCiJZCiF33vaoVY+5icep6DJ/8GcSp6zEmGS84OJimTZsybtw4mjdvTt++fUlJScm3X5s2bahVS6uH3rx5c1JTU7l7965JbFIoFIrCUuRTRlJKXyFEvftudwCCpJRXAYQQ3wGDpJTLgWeKMo8QYjwwHqBu3bp5tn1753ku3IrPs01CajoXwxIwSDAT0KSGHXbWuR/pbFbLnoUDmudr5+XLl/n222/54osveOGFF9i+fTu3b99m69atD7Tt3r07H3300T33tm/fTps2bahQoUK+cykUCkVJYOpjp7XR6shmEQp0zK2xEMIRWAa0EULMyRSOe5BSfg58DtCuXbti1/uMT83AkDmKQWrXeQlCQalfvz6tW7cGwMPDg+DgYObNm8eMGTPy7Xv+/HlmzZqFj49Pse1QKBSPBmlpaRw9ehRvb29iY2P59NNPS3xOUwtCTkdicn2ISymjgQmmmrwg3+RPXY9h+PrjpGcYsLQw48OhbfBwdSj23Nm/2Zubm5OSksLKlSvzXSGEhoby7LPPsmXLFho0aFBsOxQKRflm586dfPHFF/z5558kJiZiYWGBp6cnBoMBM7OSPQdkakEIBVyyXdcBbpl4jmLh4erA1rGdOH41mk5ujiYRg9yYMWNGniuE2NhYnn76aZYvX06XLl1KzA6FQlE2SUhIYP/+/Xh7e7Ns2TIcHBwICAjg/PnzvPzyy+h0Onr27Im9vf1DscfUguAHNBJC1AduAkOBYSaeo9h4uDqUqBAUlDVr1hAUFMSSJUtYsmQJAD4+PlSrVmr+doVCUcLcvn2bjRs34u3tzdGjR8nIyMDW1pZhw4bRtWtX3nzzTWbOnFkqtgkpi7YtL4T4FvAEnIBwYKGU8kshxFPAasAc2CClXGYiW2nXrp28vx5CQEAATZs2NdUUZZLH4TMqFI8qYWFh+Pj40KBBA7p06cLFixdp2rQpbdq0QafTodPp6Ny5M1ZWViVmgxDilJSyXX7tinPK6KVc7u8Gdhd1XIVCoSjPSCn5888/8fb2xtvbm7NnzwIwYcIEunTpQuPGjQkPDy/UTsCp6zEPZZu73Ce3UygUitJESsnly5e5cuUK/fv3RwjBhAkTCA4OpkuXLixfvhydTscTTzwBaOloCioG8anprPe9ypo/g5ASKliasXVspxITBSUICoVCUUji4+PZt2+fcRUQHByMg4MDkZGRmJub89NPP+Hq6oqdnV2hx05ITWdvQDi//R2Gb2AkaXqD8b30DAPHr0YrQVAoFIrSwmAwcOrUKVq3bo2lpSWLFy/m/fffp1KlSvTu3ZuZM2ei0+kwNzcHoEWLFoUaPyE1nX0BEfz2z20OBkaSlmGgZmVrXn7SlQbOlVi867zxqHwnN8eS+IiAEgSFQqHIkVu3buHj44O3tzd79uwhOjqagwcP0r17d8aPH8+AAQPo3LlzkYpXnboeg29gJELAhVvxHMgUgRr21vxfR1eeblWDNi4OmJlpoV2Na9gpH4JCoVA8LFJTU0lJScHBwQE/Pz86dOgAQPXq1XnqqafQ6XS0bNkSAHd3d9zd3Qs9R3TiXTYcucanB64YMyY4VLRkeMe6PNOq5j0ikJ2HdVReCUIpcvLkScaPHw9ojqlFixbx7LPPlrJVCsXjgZSSS5cuGf0ABw4cYNKkSaxatYrWrVuzYsUKvLy8aNWqVbEihEPuJON9PgyfC+H4B98xCgFo+dTGdK3P5F6NTPCJio8ShFKkRYsW+Pv7Y2Fhwe3bt3niiScYMGAAFhbq16JQlATp6enGLZ527dpx+vRpABo1asSYMWMYOHAgAJaWlgXKQ5YTUkou3I7H53w4PhfCCbitJdxsUsOOyb0a4eJgw/xfzhl9Ak82cDLBJzMNj+eTJ+QkBB+Cet3ApUOxhwsODqZ///507dqVo0ePUrt2bX755RdsbGzy7FexYkXjz6mpqao6mkJhYvR6Pf7+/sZVQExMDBcuXABg2LBhjBs3Dp1OV+ya5Rl6A/7XYzJFIIzQmBSEgPauVZn3dFO8mlXH1dHW2N7NudJD8QkUlkdLEH6fDWH/5N3mbjyEnwNpAGEG1VtAhTzyhNRoCf3fzXfqoqa/PnHiBKNHj+b69et89dVXanWgUJiItWvXMn/+fO7cuYMQAg8PD4YMGUJGRgYWFhZMmzatyGOfuh7D4cuRWFmYcSUyiX0B4cQkp2NlYUa3hk681qshvZtWx6lSzunsy0r6nPt5/J4+qXGaGID239S4vAWhgBQ1/XXHjh05f/48AQEBjBw5kv79+2NtbV1sexSKx4XU1FQOHTpkXAX8+OOPNGnShDp16jBgwAB0Oh1eXl44ORV/ayY8PpVNR67xme9Voy+gopU5fZtVp2/zGvRwd8a2Qvl9rJZfy3OiAN/kCTkJmweCPg3MreC59SbZNipq+ussmjZtiq2tLefOnaNdu3xTjigUjz1BQUG89tprHDx4kJSUFKysrOjatStJSUkADBw40OgTKCoGg+TcrTj2BUSw/2IE/9yMu+d9MwETejRgSu+y4RQuLo+WIBQElw4w8leT+hByI7/019euXcPFxQULCwuuX7/OpUuXqFevXonZo1CUV2JjY42RwR07dmTMmDE4Ojpy/fp1xo4di06nw9PTE1tb2/wHy4ekuxkcuhzF/ovh7L8YSVTiXYSAtnUdmKFrTM3K1sz9+R+jU7hLw7LjFC4uj58ggCYCJSgEBeXw4cO8++67WFpaYmZmxtq1a02yrFUoHhXeffdddu7cyYkTJ9Dr9djb21O7dm0AHBwcjA7i4nIjOpn9F8PZdzGCE1fvkKY3YGdtQQ93Z3o1qYZn42pUtf03G6mro22ZdAoXlyKnvy4NVPprheLRJTQ0FB8fH4KDg1m8eDEAvXv3JiEhwZgmumPHjkWKDL6fDL2B0zdi2XcxnP0BEVyOSATAzdmW3k2q0atJddrVc8DSvGQrlD0sSjz9tUKhUBSXU6dO8c033+Dt7c358+cBqFu3LvPmzcPKygpvb2+TnbyLTU7jYGAk+wIiOBgYSVxKOpbmgo71HRnaoS69mlSjvlPxt5zKM6UqCEKIpsDraEV29kkpS76KtEKhKBWklAQEBODt7c3LL7+Mk5MTvr6+rFmzhu7duzNq1Ch0Oh0tWrQwxuQURwyklARFJLLvYgT7AyLwv65FCTvaWuHVrDq9m1SjayMn7KyLv+J4VCjyv7YQYgPwDBAhpWyR7X4/4EO0imnrpZS5Hv2RUgYAE4QQZsAXRbVFoVCUTRITE/n999+NR0JDQ0MBcHNzY9CgQYwdO5b//ve/9wRpFoWsAjIerlVIy5DsvxjBvovhhNxJAaB5LXte7dmQXk2q8USdKjnmC1IUb4WwCVgDbMm6IYQwBz4BvIBQwE8I8SuaOCy/r/9oKWWEEGIgMDtzLIVCUY7JyMjAz88Pa2tr2rRpQ3h4OC+88AKVK1emT58+LFiwAJ1OR926dQGKVC/gfnzOh/HqN6dJ1//rD7W2NKNrQycm9mhIzybO1Kycd9YAhUZxSmj6CiHq3Xe7AxAkpbwKIIT4DhgkpVyOtprIaZxfgV+FEL8B3xTVHoVCUTqEhIQYVwB79+4lNjaWoUOH8u2339KgQQNOnDhB27ZtTeYLSMsw4H/9DgcvRXLgUiSXwhOM7wlgSNvaLHu2JdaW5iaZ73HC1D6E2kBItutQoGNujYUQnsAQoAK51GEWQowHxgPGbxUKhaL0SElJ4eLFi7Rp0waAfv36ceHCBWrXrs2QIUPQ6XT06dPH2D4rjXRxuBWbwoFLkRy4FMHRK9Ek3s3A0lzQvl5VRri58r1fCBl6LS5gWEdXJQZFxNSCkNPGXK7nWqWUB4ADeQ0opfwc+By0Y6fFsK3McuPGDZo1a8aiRYuYPn16aZujUNyDlJLz588bVwG+vr5UqFCBqKgoLC0t+eSTT3B2dqZZs2YmS9CYlmHAP/gOBwI1EQgM146F1q5iw8DWtfB0d6ZzQycqZaaJGNS69iMZF/CwMbUghAIu2a7rALdMPMcjxxtvvEH//v1L2wyFwkh0dDT29vZYWlqydOlSFixYAECzZs2YNGkSOp3O+PD39PQ0yZw3Y1M4cCmCA5ciORoURVKaHktzQYf6VfmPhwuejZ1pWK1SjqJTVpPFmYI9e/Zw4MABli1bVuJzmVoQ/IBGQoj6wE1gKDDMxHMUmzMRZ/AP96dd9Xa0rta62OMVNf01wI4dO3BzczNJyL1CUVQyMjI4ceKEcRXg5+fHnj176N27N4MHD6ZWrVr07dsXFxeX/AcrIGkZBvyC7xhFICs4rHYVGwa3qY1n42p0buBYrpPFFYd//vmHmTNn8scff+Dm5sasWbOwty9+Is68KM6x028BT8BJCBEKLJRSfimEmAx4o50s2iClPG8SSwvAeyff4+Kdi3m2SUxL5FLMJSQSgaCxQ2MqWVXKtX2Tqk2Y1WFWvnMXJf11UlIS7733Hnv27GHVqlX5f0CFwoRkpYEODAykQ4cOxMXFYWZmRocOHViwYIGxRkDLli2NpSOLS/ZVwJGgKJLT9FiZm9GhflVebK+tAho457wKeFyIjIxk7ty5bNiwAXt7e95//31effXVexJolhTFOWX0Ui73d5OLg7gskJCegMx0a0gkCekJeQpCQSlK+uuFCxfyxhtvUKlS8edXKPIjKSmJgwcPGlcBOp2ODz/8EDc3N4YPH46npyd9+vTBwcF0Wy93M/T4B8fkuAp4Vq0CcmXHjh28/vrrzJs3j6pVqz60eR+p30JBvsmfiTjDOJ9xpBvSsTSz5N1u75pk26go6a9PnDjBtm3bmDlzJrGxsZiZmWFtbc3kyZOLbY9CkZ1hw4axfft20tLSsLa2pkePHrRv3x7QooE/+eQTk80VGpOceSIokqNX1CogP/R6PRs3bmTXrl389NNPODs7ExwcXCrbyI+UIBSE1tVa80XfL0zqQ8iN/NJfHzp0yPjzokWLqFSpkhIDRbGIiopiz549eHt7c/nyZQ4fPowQAjc3NyZPnoxOp6Nbt24F8m8VlLsZevyuZa4CAiMJylwF1HGwYUjb2ni6V+NJtQrIkT/++IMZM2Zw7tw5nnzySe7cuYOTk9ODYnDjhJayv373Es3U/Fj+hlpXa12iQqBQPGx27NjBO++8g7+/P1JKHBwc8PLyIjk5GVtbW5YuXWrS+X7/5zbb/7pJbFIaF27HG1cBHd2qMrS9C56Nq9HA2VatAnLh9u3bjBo1Ch8fH9zc3Pjxxx957rnncv73urATfhwBUoKFtVbPpYRE4bEUBFNTr149zp07Z7wuSizBokWLTGiR4lEmODjY6AdYsmQJzZs3x2AwYGlpyaJFi9DpdLRr1w5zc9MFZ6Wm6zl+NZqDgZF4nw/jVmwqoAUe9W1enRfaufBkA0cqWqlHSl5kZGRgbm5O1apViY6O5oMPPmDixIk5O4ylhFOb4PeZ/5b91adpKwUlCArF40t4eDjLli3D29ubwMBAAFxcXAgNDaV58+YMGTKEIUOGmGw+KSVXIpM4GBjJwcBITlyN5m6GASsLM2pVsUGgRZyaCWhVpwq9m1Y32dyPIomJiaxcuZLvv/+ev/76CxsbG/z8/HJfQcXegF9fg6sHoGZriLwI+nSt7G+9biVmpxIEhaKMIaXk77//xtvbGxcXF1566SVsbGzYvHkzXbp0MQaGNW7c2KRbMgmp6RwJ0lYBvoGR3IzVMoU2cLZlWMe69HB3pmN9Ry7cjmf4+uPGEpKd3BxNZsOjRkZGBhs3bmTBggWEhYXxwgsvkJiYiI2NTc6/Oynh1Ebwma9dP/MBeLwCoX4PpeyvEgSFooywbds2du7ciY+PD2FhYQCMGjWKl156CXt7e2OqCFNhMEgu3I43rgJOX48hwyCxtTKnS0MnJvVsQPdGzrhUvTc1tYerA1vHdlKpIvIhIiKCXr16cf78eTp37szPP/9Mp06dcu8Qc11bFVw7CPV7wMCPwcFVe+8hlf1VgqBQlALp6ekcP36cf/75h0mTJgHw+eefc/r0aby8vNDpdPTt25datWoZ+5hCDO4kpXHociQHL0XiezmKqMS7gFYvYFx3N3q4O9O2rgNWFnmXjnyUU0UUl+joaBwdHXF2dqZt27a8/fbbDBkyJPfVnMGgrQr2aOlBeGY1eIyCUnDIK0FQKB4SoaGh7Nq1C29vb/bt20dCQgJWVla8/PLL2NnZsXXrVqpWrWpSZ3CG3sDZ0FgOXtJWAX/fjENKcKhoSbdGzvRwd6abuxPV7KxNNufjys2bN5k3bx7bt2/n0qVL1KxZky1btuTdKSY4c1XgC249YeBHUKX0sjorQVAoSojExEQOHDhAp06dcHJyYseOHbz22mu4urry0ksvodPp6NWrl7FIjLOzc5HnyqoY1snNkVpVrPENjMQ3MIpDlyOJT83ATECbug5M7e1Oj8bOtKxdGXNVNcwkJCQksHLlSlatWoVer2fKlCn5x3lcPw5HVsOVP8HcEgZ8BG1HlMqqIDtKEEqR4OBgmjZtSuPGjQHo1KkT69atK2WrFEXFYDDw999/88cff+Dt7c2RI0dIT09n8+bNjBgxgqFDh+Ll5YW7u7tJncHHr0Yx4ks/0vXa0cSsHPHV7SvQr0UNerhXo2tDJypXVLWDTU1cXBxNmzbl9u3bvPTSSyxbtsyYAypX/v4Rfh6vHSUVZvDcemg24OEYnA9KEEqZBg0acObMmdI2Q1FEIiIiSEhIoEGDBoSEhBiLxrRq1YqpU6ei0+no2rUrAE5OTjg5ORV7TiklwdHJHLwUge/lKHwDI8kw/FsqpIe7M3OeakLj6nYqMKwEkFLyzz//0KpVKypXrszrr79Oz5498y8ElJ4Kh1bBoff/jStAQHRgidtcUB5LQUj+6y+ST/pRsUN7Kmb+D1wcipP+WlG+SEtL49ixY8bAsNOnT/P888/z448/4urqyo8//kiXLl2oWbOmSedNvJvB0aAo7Ujo5Uhj8fh6jhXxaladvQHhGAwSSwszpvRuRJMaJZsm+XHlzJkzTJ8+nf3793PmzBlatWrFrFn551Dj6gHY9QbcuQoN+sD1ww8lrqCwPFKCEPbOO9wNyDv9tT4xkbsXL2rnfYWgQpMmmOeRbbRC0ybUmDs337mLkv4a4Nq1a7Rp0wZ7e3uWLl1Kt25l549DoREeHk716lrglaenJ8eOHcPCwoInn3ySpUuX8vTTTxvbPv/88yaZM/uRUN/ASE5lHgmtaGVO5waOjO/mRnd3Z1wdtZw32X0I6vSP6QkNDWXevHls2bIFBwcHVq9eTZMmTfLvmBQFPvPg7LdQ1Q1G/AJunhBy8qHEFRSWR0oQCoIhPl4TAwApMcTH5ykIBaUo6a9r1qzJjRs3cHR05NSpUwwePJjz58+XeBEMRd4kJCTw559/GlcBYWFh3LlzBysrK2bMmIEQgl69epn89xSdeJdDmVtAvpcjiUpMA6BZTXvGdnOju7sT7Vyr5ngkVB0DLTlSUlJo3bo1CQkJTJ8+nblz51KlSpW8O0kJZ7ZqYnA3EbrPgG7TwDJz1+AhxRUUlkdKEAryTT75r7+48cpoZHo6wtKSWqtWmmTbqCjprytUqGDs5+HhQYMGDQgMDKRdu3bFtkdRcAwGA1JKzM3NWb9+PZMmTSI9PR1bW1t69uzJ1KlTycjIwMrKimeffdZk86brDfx1I5aDgRH4BkZx7pY6ElpWyMjI4Ndff+XZZ5/FxsaGzz77DA8PD+rVq5d/56jLsHOqti3k0gkGfAjVCrCaKAOUqiAIITyBJcB54Dsp5YGSnrNimzbU3bjBpD6E3Mgv/XVkZKTx3PnVq1e5fPkybm5uJWaP4l/Cw8Px8fHB29sbHx8ftmzZQr9+/fDw8ODNN99Ep9PRuXNnk1epCrmTjO9lbRvoaFA0CXczMDcTtK1bhTf7uNPd3ZkW6khoqSGl5LfffmPmzJkEBARw4MABevTowXPPPZd/54y7cOh/cPh/2kpgwIfQZgSY5R3kV5YoTgnNDcAzQISUskW2+/2AD9FKaK6XUr6bxzASSASsgdCi2lJYKrZpU6JCUFB8fX1ZsGABFhYWmJubs27duodaHelx5Pbt2zz11FPGk13Ozs707dsXR0ctH0+bNm2MJ4WKQ9aefhuXKtzVGzIjgyO5GpkEaBXDnnmiJj3cnXmygROVbdSR0NLmr7/+MjqM3d3d2bFjB927dy9Y52uHNKdx9GVo+R/QvQOVqpWswSWAkFLm3yqnjkJ0R3uYb8kSBCGEORAIeKE94P2Al9DEYfl9Q4wGoqSUBiFEdeB/Usrhec3Zrl076e/vf8+9gIAAmjZtWqTPUF54HD6jqZFSEhQUZPQDNG3alBUrVmAwGBg0aBBPPvkk/fr1o3Xr1piZ8BuclJIdf91ixraz9xwFrZCZBK67uzM93J1UxbAyRnp6OvXr1+fu3bssWrSI8ePHFyxVSPIdLRHdma+hiis88z9o2KfkDS4kQohTUsp896KLU1PZVwhR777bHYAgKeXVTCO+AwZJKZejrSZyIwYo+QrSiseCBQsW8PXXX3Pt2jUA3NzcjGfEzczM2Llzp0nni01O43BQlDE6OCw+1fieAJ73qMOSwS2wtjRdSgpF8YmPj2fdunVMnToVKysrfv75Z9zd3alcuXL+nW+cgONr4cp+SE+Grm9A95lgVTH/vmUYU/sQagMh2a5DgY65NRZCDAF0QBVgTS5txgPjAerWLb0cH4qyh8Fg4PTp08Z4gG3btiGEICYmhhYtWjBt2jR0Oh0NGzY06bxafqA4fDOzhP4dGotBgr21BV0bOeFatSIbjgSTodfSQw/tUFeJQRkiPT2dL774gkWLFhEZGUnLli3p37+/scZ0vvz9A/z838xIYwGD18ETQ0vW6IeEqQUhpzVwrntSUsqfgJ/yGlBK+TnwOWhbRsWyTvFIcOzYMT7++GP27NlDVFQUoO39R0dH4+TkxMcff2zyOW8tbzuOAAAgAElEQVTHpRgF4PDlKOJTMxACnqhThcm9GtHD3Zkn6lTGwlzbfurTrIaKCyhjSCnZtWsXM2fO5OLFi3Tv3p3du3cX/FRfahwceE9bGRgfa2YQf7OkTH7omFoQQgGXbNd1gFsmnuMBpJSP7H5sUX08jwp3797l8OHDeHt7M3LkSJo3b05YWBj79u2jX79+6HQ6vLy8jIFjpiI1Xc/Ja3eMgWGXMwvHV7evgK55DXo0dqZrQyeqVLTKsb+KCyh7SClZuHAhUkp++eUXBgwYULDnhsEAf3+vpadOioTG/bWtoocUaXzz5k22b9/OlClTSnQeML0g+AGNhBD1gZvAUGCYiee4B2tra2P+8UdNFKSUREdHY239eJ1DT0hIYNOmTfzxxx8cOHCA5ORkLC0tadWqFc2bN2fAgAEMGjTI5M7goIjEzNQQUfeUjOxYvyovtHOhu7sz7tWVM7g8cePGDZYuXcry5ctxdHTkl19+oUaNGgWvLXH7LOyeASEnoHY7GP4D1Grz0CKN/fz8GDRoEAkJCTz77LO4uLjk36kYFOeU0beAJ+AEhAMLpZRfCiGeAlajnSzaIKVcZiJbczxllJ6eTmhoKKmpqbn0Kt9YW1tTp04dk1bKKmvExcWxf/9+LC0teeaZZ0hKSqJq1arUrVsXnU6HTqejZ8+eVDJBRPk98yanc+RKljM4kltx2t9QA2dburs7093dmU71HbGxUvv/5Y34+HiWL1/OBx98gBCCbdu23ZNiJF+S78D+peC/ASo6gtdieOKlhxpT8P333zNq1Chq1KjBr7/+SsuWLYs8VkFPGRVZEEqDnARBUT45ffo0u3fvxtvbm2PHjqHX6+nVqxf79u0DICwsjBo1aphsvlPXYzh2JQrHShWIiL/LwcAIzoRozmC7ChZ0aeiUKQJO1HEo3ydFHmeklHz66adGh/HLL7/M0qVLC34gxaCH01tg32LNZ9BhPHjOBpt8UlWYmGXLljFv3jy6du3KTz/9VKxaGfAQjp0qFIXh1q1bnDx5ksGDBwMwf/58du/ejYeHB7NmzUKn0/Hkk08a25tKDMLiUvnqWDCfHrxCtrAAnqhTmVd7NqSHuzOtXaoYncGK8o0Qgr1799K8eXPef/992rZtW/DOIX6wezrcPgOuXaD/CqjRIv9+JUDDhg0ZPXo0a9euNXm0fF6oFYKiREhNTTU6g729vfnnn38ALWVEtWrVCAwMxMHBodjffB6YN12PX/AdY0zApfCEe943EzDJsyHTdY1NOq+i9PD392f27Nl88sknNG7cmKSkJCpWrFhwX09iBOxdpCWjs6sJfZdCi+ceevWy0NBQTp06xaBBgx5473T4aU5HnKZd9Xa0rta60GOrFYLioSKl5NKlS1SvXh0HBwc2b97MhAkTsLS0pGvXrrz77rvodDpjgRh3d3eTzXslMpGDgZov4MS1aFLTDViZm9G+vgND2jbB2a4Cc3/+h/QMLS6gZ5Pyl1JA8SDXr19n7ty5fPPNNzg7O3Pt2jUaN26Mra1twQbQZ4DfevjzHS24rMvrWlbSCnYla3gOZK2e09LS7imrCnDy9knG+YxDIqlgXoEv+n5RJFEoCEoQFEUmNjaWffv2GVcBN27c4Msvv2T06NEMGjSI2rVr4+npaXpncEo6R4OiMpPERXEzVisW4+Zky9D2denh7kxHt6pUtPr3z9vV0VbFBTxCLFiwgBUrViCEYO7cucyaNatw6ciDD2unhyIuQINe2vaQU6OSMzgPvvvuO1555RVq1KiBj4/PPWKQmpHK28fexoBWYS3dkI5/uL8SBEXpo9friYmJwcnJiejoaGrUqEFGRgZ2dnb07t2bOXPmoNPpAM0H8MwzeWUrKcS8BsnfobH4BmoicCYkFr1BYlfBgs4NHZnUswHdGznjUjV3Z7CKCyj/ZGRkYGGhPbKSk5N58cUXWbp0aeGOYsbf0nIPndsGlevCi1uhydOlUtw+Ky5iyZIldOvWje3bt9+zhZqSkcKU/VO4kXADCzMLpJRYmlnSrnrJpcdXgqDIk5s3bxpXAHv37qVbt27s2LEDR0dHVq5ciYeHB506dTL5sdiwuFR8L2uRwUeCoohNTkcIaFm7MpM8G9A90xlsqZzBjzxZgWSzZs1i7dq19O7dm5UrVxYuHiQjTYswPrgCDBnQY7a2RVSKuYeEEOj1ekaPHs2nn36KldW/QY4pGSm8tv81Tt4+ydIuS3G1d8U/3L/IPoSCogRBcQ/Zv4UNGzaMb7/9FtC+8WcFhGUxdepUk82bmzPY2a4CvZtUp7u7E90aOVPVNufIYMWjycmTJ5k+fTqHDh2iadOmxr/NQolB0D74fSZEB0Hjp7TU1FXrl5DF+RMaGkpYWBjt2rVjyZIlCCHu+TwpGSm8tu81ToadZGnXpQxsMBCgRIUgCyUIjzlSSi5evGhcBfj5+RESEoKNjQ19+/alTZs26HQ6WrZsadIIXSklO8/eYseZW8QkpxFwO/4BZ3B3d2ea1LBTkcGPKa+99hpr1qyhWrVqrFu3jjFjxhgFoUDEXAfvuXBxl1bPeNiP4N635AwuACdOnGDw4MHY2dlx4cKFBz5Pcnoyr+1/Db8wP5Z1XcaABgMeqn1KEB5jdu3axaRJkwgJ0RLUNm7cmOHDh5OcnIyNjQ2jRo0y6XzZncF7LoQbawYLoH/LGjzvUYdObo73OIMVjxexsbHY2dlhbm5Oy5YtmTdvHjNnzrzH0Zov6Slw5COtcpkwg94L4MnJYFG6Gfa/+eYbRo8eTe3atdmxY0eOYjB5/2ROhZ8qFTEAJQiPBXq9Hj8/P+MqYN68eTz11FPUqlWL9u3bM2/ePHQ6Ha6uriabMysy2MlOiwz2DYzkr2zO4BqVrYlOTEOixQY0r1WZXk1Mm6BOUX5IS0tj3bp1LF68mBUrVjB69GjGjx9f8AFCTmpVy8zMtXQTsdeh+bNaTEHlOiVneAEwGAwsWLCAZcuW0aNHD7Zt22Y8fp1FdjF4p+s7PO1WiDQbJkQJwiNMXFwc48aNY+/evcTExCCEoH379sYMqm3btmX79u0mnTM8PpUtx4L59MC9kcGt6lRmYg/NGdymbhX+Do1j+PrjxtiATm6OJrVDUT6QUvLzzz8za9YsgoKC6N27Nx4eHoUbJOQkbH5Gq2kMWuWykTuhfgHLXz4E/v77b8aOHcsnn3xyj/MYNDGYtG8Sf0X8VapiAEoQHhlSUlLw9fXF29sbBwcH5s+fj729PVevXmXw4MHodDr69OljrB1sKlLT9fgHx3AwMKJQkcEerg5sHdtJxQY85owZM4aNGzfSrFkzdu/eTb9+/QrnM0qJ0Y6RZokBAtq8XCbEICQkBCEEderUYdu2bVhaWj7w2ZLTk5m4dyJnIs/wbrd36V+/fylZq6EEoZyzadMmvv32W3x9fUlNTaVChQoMHapVbxJCYOpUH1pkcJJ2GuhyJMev/hsZ3K6eA7PbNsG5UgXe2pF/ZLCKDXg8uXbtGk5OTtjZ2TF06FA6derE6NGjC+cw1mfAqY1alHHKHRCZGWnNrcCtR8kYXghOnDjBoEGDaNq0KX/++ecDq4IzEWc4euso+27s40rsFd7r9h796vcrJWv/RQlCOeLOnTvs3buXw4cPs3r1aszMzDh69Cg3btzgv//9L/369aN79+5UrGjas9X5RQZ3d3d6wBlcz0lFBivuJSYmhnfeeYePPvqIWbNmsXjxYvr2LcKpnyv74Y+5EBmg1SPQvQMZqQ+lPkFB2Lp1K2PGjKF27dqsXbv2gffPRJxhrM9Y7uq1Vc2rrV8tE2IAShDKPEFBQXz11VfGI6EGg4EqVaowbdo0XF1d+eSTT0weFKYigxWmJC0tjU8//ZTFixcTExPDK6+8wn//+9/CDxR1GXzmQeAf4FAPXvwamjzzb5RxKQuBwWBg/vz5vPPOO3h6erJt27Yct2gP3zxsFAOBwMKs7DyGy44lCkCr8OTt7U23bt1o0qQJAQEBLF26lA4dOjB//nx0Oh3t27c3Lq9NJQZhcalazeDL90YGt1KRwYpiMmHCBDZu3IiXlxcrV67kiSeeKNwAKTFwcCWc/AwsbLRiNR0nlPox0vtJTk7ml19+Ydy4caxZs+aBbSKAyORIdl3dBYAZZliZW5VoKorCUqrpr4UQ3YDhaMLUTErZOa/2j2L66/T0dPbu3Ws8Enrx4kUAVq5cyfTp00lNTSU5OZmqVauadN6smsFZvoDAcK1mcDW7CsZqYV0bOqnIYEWROH78OLVr18bFxYVz584RGhqKTqcrnMP4Hj9BDLQdAb3mQaWyla02JCQER0dHKlasSFxcHPb29jl+zpuJNxnnM46olCheb/s6KRkpJZ6KIosSr5gmhNgAPANESClbZLvfD/gQrYTmeinluwUYazBQXUr5WV7tHgVBkFJy7tw54uLi6Nq1KykpKcaHfY8ePYwlI5s2bWrSCN1TwXf47dxtDAbJ1ajkf2sGZ0YGd2/kTI/GzjSuriKDFUXn6tWrzJkzhx9++IGJEyfmuIdeIHLyE9RsZVpjTcCxY8cYPHgwAwYMYP369bm2uxp7lXF7xpGakcqnfT6llfPD/SwPox7CJmANsCXbpObAJ4AXEAr4CSF+RROH5ff1Hy2ljMj8eRgwthi2lGmio6PZs2cP3t7e+Pj4cOvWLTw8PPD398fGxoZDhw7RvHlzbGxsTDpvVs3gn06Hsjcgwni/VhVrXuqQc5pohaIo3Llzh2XLlvHxxx9jaWnJwoULmT59euEHys9PUIb4+uuvGTNmDC4uLkybNi3XdheiLzBhzwTMhBkbdBtoXLXsFmcq8pNASukrhKh33+0OQJCU8iqAEOI7YJCUcjnaauIBhBB1gTgpZXwu748HxgMFr4taymRkZHD27FljgM2YMWP45ZdfcHBwwMvLC51Od8/pinbtTLOHqDdIzobGGovGZ9UMtrL4d9/fTMDwjq682rOhSeZUKADefvttPv74Y0aPHs3ixYupVatW4QYoJ34C0JzHc+fO5b333svTeQxwKvwUk/dNxt7Kni/6fkFd+7L9DCuWDyFTEHZlbRkJIZ4H+kkpx2Zevwx0lFJOzmOMtwFvKeXR/OYry1tGwcHB+Pj44O3tzb59+4iLi+PmzZvUqlULPz8/9Ho97du3x9zc3KTz3o5LMWYIPRwURVxKpjO4ThV6NNIKxxukZMSGk8a4gK1jO6kTQIpiIaVk27Zt1K9fn3bt2hEeHk54eDitWhVyK6Sc+AmyExoayhNPPMHzzz/PmjVrcj3YcSj0EG8ceINalWrxudfn1LA1TZ3wolBaJTRzWtflqThSyoUmtuGhkJSUBICtrS3fffcdL730EgAuLi785z//QafTUblyZQDat29vsnlT0/WcyHIGB0ZyOUJzBle3r0DfZtWNzmCH+5zBKipYYSqOHj3KtGnTOH78OK+88gobNmygevXqVK9eyFxU5cRPkEVERATOzs7UqVOHs2fPUrt27Vz9bd7B3sw+NJtGVRqxzmsdVa1NeyikxJBSFvkF1APOZbt+Eu3bftb1HGBOcebI/vLw8JClhcFgkGfPnpUrVqyQvXv3llZWVnLdunVSSilv3bolP/jgA3nhwgVpMBhMPu+lsHj5he8V+X/rj0v3t3ZL11m7ZKO3dsvhXxyXnx+8Ii/ejjf5vArF/QQFBcnnn39eArJmzZpy/fr1MiMjo/ADRQZKufUFKRfaS7m6lZQXfpWyjP/9HjlyRFarVk2+9957+bbdHrhdttrcSo7YPULG340v9tyXLl2SI0eOlElJSUUeA/CXBXjGmnqF4Ac0EkLUB24CQ9EcxuWSrGIxycnJuLu7c/PmTQCaN2/O5MmT6dBBC4SpWbOmSYvFxCancThIKxp/6HIUt+NSAWjgbMuwjnXp7u5Mp/qO2FiZdvtJociLH374gd9//523336badOmFbyYfRblyE+QnS1btjBu3Djq1q3LwIED82y7+fxmVvmvokvtLnzg+QE2FsU7KHL79m10Oh1JSUmEhYXh5uZWrPHyo8iCIIT4FvAEnIQQocBCKeWXQojJgDfayaINUsrzJrH0IZCens7x48eNMQE1atRg586dVKxYkeHDh9O4cWP69u1LnTqmTaeboTdwNjTOGBNwNtMZbGdtQdeGTkzprcUF1K5i2lNICkVe3L17l08++QQ3NzcGDx7M1KlTGTVqFDVr1izcQOXQTwBa2vi5c+eyYsUKevXqxY8//phrPJCUkrVn17Lu7Dr6uvbl3W7vYmlevKDRuLg4+vfvT2RkJAcOHChxMQCKt2X0sF8luWU0d+5caW9vLwFpZmYmO3fuLFesWFFi892MSZbfnrguJ37tL1su/EO6ztol68/eJQetOSzf97kk/YOjZXqGvsTmVyhyw2AwyO+//17Wr19fAnL8+PFFHyxon5RrOmrbQxuflvLWWdMZWsKcOnVKmpubywkTJsi0tLRc2+kNern8xHLZYlMLOf/wfJmhL8I22n2kpqZKT09PaWFhIb29vYs9HqW0ZVTmSUxM5MCBA3h7e3PgwAGOHz+Ora0tNWvW5MUXX0Sn09G7d2+qVKli0nlT0/UcvxptzA8UlOkMrmFvTb8WNYzO4CoVVWSwovQ4efIkU6ZM4cSJE7Rq1QofHx+8vLwKP1A5iie4n8TERCpVqkTbtm05c+YMzZs3z9V5nGHIYNHRRfxy5RdebvYyM9rNMElg57Vr17hw4QIbN24sWgLAIvLYCMLx48d56623OHToEOnp6djY2ODp6UlUVBS2trZMnpzrydgiIaUkMDzRuA104tod0jIMWFmY0bF+VYa2d6G7uzONqlVSkcGKMsPVq1e5ceMGGzZsYMSIEYU/Jl1O/QRZHDlyhOeee47PP/+cgQMH0qJFi1zbpunTmOU7i7039jKp9SQmtJpgsv+XmzRpwuXLl7G3twcgfu9e0oKCqNixIxXbtDHJHDnx2AhChQoViIyM5PXXX0en09G1a1esra1NOkdM0r3O4LB4zRncqFolXu7kSnd3ZzrWr4q1pXIGK8oG0dHRLFmyxBht++KLLzJgwIDCO4z1GXB6E+xfVq78BNnZvHkz48ePx9XVlcaN844mPnHrBIuPL+ZGwg1mtZ/F/zX7P5PY8O6775KYmMiSJUuMYpB4+DA3X5sCUiKsram7cUOJicJjIwht2rTh77//NumYmjM4loOXIjl4OYq/Q2OREuytLejayInujTRncC3lDFaUMVJTU1mzZg1Lly4lISGBN954A9CKKhVaDMpZPMH96PV65syZw8qVK+nduzc//vgjDg65x+ocCj3Eq/teRSKxMLOghVPuq4jCsGnTJubMmcOwYcOQUiKEQEpJ+PJ3ITOAWKank3zSTwlCWeFmbIoxKOxwUBQJqRmYCXjCpQpTejWiu7szT9SpjIVKE60oo+zdu5dx48YRHBzMU089xYoVK2jevHnhBgk5Ced3wM1TEHK83PkJsrNr1y5WrlzJpEmTWL16dZ4p5UMTQpl7eC4yM95WSol/uH+xM5b+9ttvjB07Fi8vLzZu3IiZmfb8uLN5M2lXroCFhbZCsLSkYgfTBbrejxKEfEhJ03P8WrRRBK5EahHKNStb81SLmnR3d6ZLQ0flDFaUefR6Pebm5lSqVIkqVaqwZ88e+vTpU/iBgvbCNy+CIUO7bj9WWxWUEz9BFllxRgMHDuTPP//E09Mzz/YXoi8wae8k0vRpWJlZoZd6LM0si13P4Pjx4/znP/+hdevWbN++3VhHIeXMGSJWvY+dVx+qvjKaZD8/KnZor3wIDxMpJZfCE4z5gU4Ga87gChZmdHRzNGYJbaicwYpywuXLl5k9ezbOzs6sW7eOTp06cfr06cL//erTwe9L2LvwXzEQ5mBfq9yJweHDhxk5ciS//PILLVq0yFcMjt48yhsH3qByhcps0G0gPi0e/3B/k9QzuH79OvXq1WP37t3Y2dkBkBETQ+gbb2JZowY1ly3D3N6eim1LTgiyUIKA5gw+ZHQGRxIer5W3c69eiRGZzuAOyhmsKGdERUWxZMkS1q5di7W1NW+99ZbxvUKJgZRw6XfYMx+ig6Bma4gI0ETB3ErzG5QjNm3axPjx46lfv36OVc3u59crv7LwyEIaVGnA2j5rqVZRc5QXVwgMBgNmZma8+OKLDBkyxLhVJQ0Gbs2ejT4qCtdvvsE807n8MHgsBSFDb+CvkH/TRP99Mw4pobKNJV0bOdGjkTPd3J2oWVk5gxXlk99++43hw4eTkJDAuHHjWLRoETVqFCHb5u2z4P2WVsDeyR2G/QCN+kKoX5kpal9Q9Ho9s2fPZtWqVfTp04cffvghT+exlJIvz33Jh6c/pGPNjqz2XE0lq0omsSU2Npa+ffsya9YsnnvuuXv8FtFffknSQV+qz5+HTUvTOKwLymMhCKeux+B9PgwBBEcncTQomoS7mjO4tUsVXu/diB7uzrSqUwVzM7UNpCifGAwG4uLicHBwoEWLFvTs2ZNly5bRrFmzwg8Wfxv2L4UzW8HGAZ5aBR6jICsdg0uHciMEWXz66aesWrWKyZMn88EHHxjrkueE3qBn+cnlfH/pe552e5olnZcUOxVFFqmpqQwaNIgzZ84Yj5ZmkXzqFJGrP8SuXz8chj38NHClWlO5sBSlHsKp6zEM/fwY6XrtczrZWtEnM010lwZOVK5oml+yQlGa+Pr6Mn36dBwcHPD29i76QGlJcPRjOPKhtiXUcQJ0mwY2po3cf5hkHeFMS0tj586dPPfcc3m2T81IZZbvLPaH7OeVFq8wte1UzIRpTg3q9XpefPFFtm/fzjfffGNMmw+QcecO1wY/i7Cxpv727ZhXMs1qBApeD+GRPxt5/Go0eoMmBmYCXulaj3efa8VTLWsqMVCUewIDAxk8eDA9evTg1q1bDB8+nCJ9yTMY4Mw38LEHHFiubQu9ehL6LinXYuDr60vnzp2Jjo7GysoqXzGITY1lnM84/gz5k9kdZvOmx5smEwMpJVOmTGH79u3873//u0cMpMHArRkz0cfGUmf1apOKQWF45LeMOrk5YmVhZqwW1snNqbRNUihMwq5du3j22WextrZm2bJlTJ06lYoVKxZ+oGuHwOctzV9Qqy38ZxPU7WRyex82X375JRMnTsTNzY24uLhcy1xmcSvxFhP2TuBmwk1W9VhF33qmzyFkY2PDjBkzjIGAWUR/9hlJR45Q4+23sW7a1OTzFpiCZMArK6+iZjv1D74j1+y/LP2D7xSpv0JRVkhOTpaXLl2SUkoZHx8vp02bJsPCwoo2WFSQlN8O0zKRvt9MyrM/SKkv/xl2MzIy5JtvvikB2bdvXxkTE5Nvn4DoANnz+57yyW+elH63/UxuU0pKivHn+4tZJR47Li80bSZDp88osUJXFDDbaak/5AvzKs2KaQpFaaLX6+XXX38t69atK5s0aVK0SmVZJEVL+ftsKd+uKuWyWlIeXCllWrLpjC1l5syZIwH52muvyfT09HzbH7t1THbc2lH2/qG3vHznssnt+fXXX2WdOnVkQEDAA++lR0TIS126yqD+T0l9YqLJ586ioIJQqltGQohmwCIgGtgnpdxWmvYoFGWRgwcPMm3aNE6dOkXbtm1ZtWpV4bOQAmSkgd96OPge3I3XEtB5zgW7QtZCLuNMnTqVxo0bM3LkyHzb7rq6i/lH5lPPvh6f9vmUGrZFOJqbB0ePHuWFF16gRYsWDxTWkno9N6fPwJCYSN0NX2JW2BxSJUFBVCOnF7ABiCBbTeXM+/2AS0AQMDufMaYB3TJ//jW/OdUKQfG48fvvv0tAuri4yK+++krqi7KlYzBIeWGnlB+21raHNg+SMuyc6Y0tRQ4ePChfeOGFPAvZZMdgMMgN/2yQLTa1kK/88YqMuxtncpsuXLggHRwcZMOGDWV4ePg97yWdPi2DR4yUFxo3kTHbfzL53PfDQ1ghbALWAFuybgghzIFPAC8gFPATQvyKVk5z+X39RwNfAQuFEAOBvD0+CsVjQmRkJBcuXKBHjx54eXmxdu1aRo0ahY1NEQIlb53RAsuuHwanxjB8GzTsU+4S0OXF+vXrmThxIg0aNCAyMpJatWrl2V5v0LPSfyVbA7bSr14/lnVdhpW5aXOR3bx5E51Oh5WVFd7e3lSr9m8a8OS//uLGiJHI9HQwN8eqfj2Tzl0ciiwIUkpfIUS9+253AIKklFcBhBDfAYOklMuBZ3IZ6tVMIfmpqLYoFI8CKSkprF69muXLl2Nra8v169exsrJi4sSJhR8s/hbsWwJnv4WKVeHp96HtKDB/dA4W6vV6pk+fzurVq+nbty/ff/99vpUO7+rvMufQHPZc38OIZiOY1m6ayY6VZsfe3p4OHTrw1ltvPVALOcFnjyYGmZRkOuvCYuq/jtpASLbrUKBjbo0zBWUuYAuszKXNeGA8QN26dU1kpkJRdjAYDGzdupW33nqLkJAQBg4cyHvvvVegPDsPcDcRjn4ERz4CqYcuU7TAMuvKpje8lJkwYQLr169nypQpvP/++3lGHgPE3Y1jyv4pnI44zfR20xnZPH8fQ2FJSUlBr9djZ2fHtm0PukQNyckk7N2rXZiZlXg668JiakHIaR2aa5SMlDKYzId9Hm0+Bz4HLVK5OMYpFGWRY8eOMWLECDw8PNiyZUu+mTdz5PoxOL4Wgg9Dyh1oPgT6LNTqFDyivPrqq3To0IFx48bl2zYsKYwJeyZwI+EGK7qvoH/9/ia3R6/XM2zYMMLCwjh06NADAiWl5NZbb5F+8ybV5sxBpqaWeDrrwmJqQQgFXLJd1wFumXgOhaLcExAQwIkTJxg1ahRdunRh79699OzZ01gYpVCc+Ax+nwVIzTfwzGpo94rJbS4LHDhwgH379rFkyRJat25N69b5ZxwNjAlk4t6JJKcns67POjrUNH0OJiklr776Kjt27ODDDz/McbVy5+7uSV8AACAASURBVMsvSfj9D5ynvYnjyBEmt8EUmHrzzA9oJISoL4SwAoYCv5p4DoWi3BIREcHEiRNp2bIlM2bMIClJK7jUu3fvwotB5CWtUM3vM/l3IW6mrRAeQb744gu8vLzYvn078fHxBepz8vZJRv6ubQ1t7r+5RMQAYMmSJXz22WfMnj2bKVOmPPB+4uEjRPzvA+z698Nx7NgSscEkFOQoUk4v4FvgNpCOtjIYk3n/KSAQuAK8VdTxc3qpY6eK8kpSUpJctmyZtLOzk+bm5nLy5MkyIiKiaIMlREi58w0pFzlI+U4dKXe+KeWSatr1kupS3jhhWuNLmfT0dPn6669LQPbv31/GxsYWqN/vV3+Xbba0kYN+HiRvJ94uMfs2btwoATly5MgcI43v3rghL3boKK8MGCj1SUklZkdeoCKVFYqyw8WLF6WFhYUcPHiwvHjxYtEGSUuR0vd9KZfV1h7+u6ZJmRipvXfjhJS+qx45MZBSyiFDhkhAvvHGGwWO0N50bpNssamFHLF7hIxNLZiAFJUrV67IiRMn5hgDoU9KklcGDJQXO3SUd69fL1E78qKggvDIp79WKEqL/fv3s2fPHpYv10Jwrl69+sARxAJhMMC57f/f3p3H2VT+ARz/nNmHWTAYhJkaZd+JZMu+7yTaUKkfIpXIEkUKFdmybyFCWUqUfck+RiGMpZkxGLPvc+fe+/39cSeJGXNnuffO8rxfr/NqzpnnnPudp+t+7znPBnumQEwwPNUR2n4MpZ7K5Yjzpk2bNhEdHc1rZjxqMYqRWadmsebCGtr6tGV6s+k421tmec9r167h6+ub4aM+EeHm6NHE7dpNhUWLcGvW1CJxmMPc6a9t/q0/K5u6Q1Dyg/Pnz0vnzp0FEB8fH4mMzMGkijeOiixqaRphvLCpyNX9uRdoHrZ3715ZtWqV2eX97/jLwrMLZcgvQ6TGyhry6bFPRW/IwXxPmfjzzz+lWLFi8uGHH2ZYJnzpUrlQuYrcXbzYYnGYC/XISFGsKzw8XIYOHSp2dnbi6ekpn3/++X9muczaxQJFvhtoSgSzqoj4ry0QM5Ga45tvvhEHBwepU6eOWZPT+d/xl/pr6kuNlTWkxsoaMvX3qRabNVREJCgoSMqXLy9ly5aV69evp1sm7vBhuVC1mgSPHGXRWMxlbkIoOMMWFSUP2LJlC8OGDWPSpEmULJmNtTcSI+HADNMkdPZO8Nx4eGY4OGVjnYN8Rq/XM3r0aObOnUvHjh357rvvMh1sBrAnaA8phhQANDS8i3qjWWhqjsjISDp06EBsbCyHDh3C19f3oTK64GBujn4X50qVKDdtqsVisQSVEBQlmwwGA2vWrGHLli38+OOPeHl5ce3aNdyys9qVPgVOLIGDMyAlDuq+BM99CO65O/tmXmUwGOjSpQu7du1i9OjRzJgxw6wZXc/dPcfmK5sBsMMOJ3snGnhn/qg8O0SEvn37EhgYyK5du6hVq9ZDZYyJiYQMHwFA+Xlz88YMplmgEoKiZMOePXt47733OHv2LA0bNuTu3bt4e3tnPRmIwIWt8NtHEHUD/Fqblq30rm6RuPMqe3t7mjZtSt++fRkyZIhZ5+y+sZsPD39IKddSTGw8kZvxN2ng3YA6pTMfrJYdmqbx4YcfEh0dne5ochHh1oQJpFy5QoVFi3DKj1PtmPNcKa9sqg1BsbU7d+5Ip06d7jUYr1+/PntTUouIBJ8UWdrO1E4wv7HIlV9zN9h8YM+ePXL48OEsnWM0GmXJuSVSY2UNefGnFyUiKcJC0f37eseOHcu0XF5qRH4QZrYh5P40f4pSABkMBgCKFStGWFgYM2fO5K+//qJ///5ZH2Ec9TdsGgxLW0PkNeg6B4YeMk1LXYgsXLiQdu3aMWHCBFMPFzOkGlOZ/Ptk5pyZQ0ffjixtv5QSLiUsGufkyZNp3LgxR44cybBM/JEjhH3xJe4d8vhI5MyYkzXyyqbuEBRri4+PlylTpkjlypUlPm2Jw2z3GkmMEtk1QeTjkqYRxXs+EUmOzcVo84fU1FQZNmyYANKpUyeJiTFvcZqYlJh73Uq/PvO1GIyW73W1cOFCAWTw4MEZ/n9PCQ6WS/+MRLbgMpg5gep2qijZp9frZdmyZVKuXDkBpFevXtlfzF6vEzm2SOQzX5GPPEW2DBWJDsndgPOJ+Ph4adu2rQDy7rvvmj3yODg2WLr90E3qrK4jP1750cJRmmzevFk0TZMuXbpk2P3VkJAgV7t1l78aPm3TkciZMTchqEZlRXlAREQErVq14ty5czRq1IiNGzfy7LPPZu0iwSfg+iGwswf/byHiCvg2g3ZToZxlGj3zA1dXV0qUKMGyZcsYPHiwWecE3A3g7b1vozfqWdx2MQ3LWH79gBs3bjBgwAAaNWrEhg0b0u3+KiLcmjCRlMuXqbA4nzYiP0AlBEVJExkZSYkSJShRogR169Zl/Pjx9O3bN+v9yINPwMoukNY3Ho/y8MJ38FSHArV0ZVbs3bsXPz8/fHx8WL9+vdl1uuvGLsYfHk/pIqWZ33o+j3s+buFITXx9fVm4cCHdunWjSJGHx4Ak+vsTvmAhCYcOUWr0aNyaNbNKXJamGpWVQu/WrVu8/vrr+Pj4cPPmTTRNY+XKlfTr1y/rySDmJvz07r/JAA3qvwqVOxbaZLBgwQLatWvH2LFjAcyqUxFh6R9Lee/Ae1TzqsbaTmutkgyCgoI4c+YMAIMGDcLL6+Gl3v9ZEznh0CGws6NIA8uMe7AFdYegFFoJCQnMmjWLmTNnotPpGD58eLrfBs2SEgdH5sDReWBMBTsH0xgDeyd4okXuBp5PpKamMmrUKBYsWEDnzp1ZtGiReecZUvnk2Cf8EPgDHR/vyCfPfmKxCeruFxERQfv27UlISODKlSs4O6f/mrG/7Pp3TWRNI/HkSYrUyzurnuWESghKoRQXF0e1atUICQmhb9++TJ8+HT8/v6xfyGgA/zWwdxokhEGN3tD6I4i/AzcOmdoNKlhmUZa8LDo6mj59+twbwPfZZ5+ZNfI4JiWGd/e/y/HbxxlaayjD6gyzytQPiYmJdO3alevXr7N79+4Mk4E+KorYXb+YdvLgmsg5pRKCUqj88ccf1KxZE3d3d0aMGEHTpk1p0qRJ9i4W+BvsnghhF6BCI3hhPZRPe3xQ3KdQJoJ/2NvbExMTw4oVK3j11VfNOic4Lphhe4YRHBfMtKbT6ObXzbJBptHr9Tz//PMcP36c77//nubNm6dbTnQ6bo4chTEikjKTJ2OIiclzayLnmDldkXJjA54AlgGbHnXsUZvqdqpk17lz56R9+/YCyJkzZ3J2sdvnRVb3NI0wnl1L5M8fRPLAjJZ5waFDh+6N1zC3S6mIacbS5t81lybrmsiJWycsFV66vv76awFk4cKFGZYxGo1yc/x4uVC5ikRv22bF6HIHuTkOAVgOhAF/PnC8A3AJCATGmnmthz78VUJQLOXmzZsyZMgQsbOzk2LFismXX34pycnJ2btY7G2RrSNEJhcTmV5B5MhckdRsXqsAmjdvntjb28v777+fpfN2Xt8p9VbXk46bO8q16GsWii5jOp1OtmzZ8sgy4cuWy4XKVeTO7NlWiip35XZCaA7Uuz8hAPaY1k1+AnACAoBqQE1gxwNbaXnEh79KCIolJCUlSenSpcXR0VFGjx4tERHZnPMmJUFk/wyRaeVEppQQ+fkDkQTLzp+Tn+h0OnnrrbcEkK5du0psrHmjr41GoywOWCw1VtaQl35+SSKTcrCQUDZ89913Zg02jN2zVy5UqSrBb48UYz5dkyJXE4Lpevg+kBCeAXbdtz8OGGfGdbKUEIA3gFPAqYoVK1qswpSCQa/Xy5YtW+5NM7Bp0ya5evVq9i5mMIj4rxP5oqrp8dD6AaaFa5R7IiIipHXr1gLImDFjzH5MpNPrZMLhCVJjZQ0Zc2CMJOute6e1ceNG0TRNhg0b9shySRcvysW69eRa7z5iSEy0UnS5zxoJoQ+w9L79l4B5jzjfC/gm7a5iXEbHHrWpOwQlI0ajUXbu3Ck1atQQQH79NYczh147KPJNM1MiWNRC5HrWZuQsLK5duyblypWTlStXmn1OdHK0DP5lsNRYWUPm+c+z+opi+/btEycnJ2nSpIkkJCRkWC41LEwut3xOLjdvIbrbd6wYYe4zNyHkpJdRen3BMpyyUEQigDczO6YoWRUQEMD777/Pr7/+ip+fH5s2baJ169bZu1h4IPw6CS79ZBph3GsJ1OgDWZ3RtIDz9/enTp06PP744wQGBuLq6mrWebbqSfSPgIAAunfvjp+fH9u3b89w3IkxOZng4cMxREfju/ZbHL1LWzVOW8nJuzwEqHDffnkgNGfhKErW6PV6unXrxunTp5k9ezYXLlygd+/eWe+7nhABP4+BBY3g+kFoPQlGnIJa/VQyuI+IMHfuXBo2bMjChQsBzE4GZ8POMvCngUQkRbC47WKrJwOAUaNG4e7uzq5duyhRIv1ps0WEWx+OJ/ncHzw2cwYu1apZOUobMuc2QtJ/ZOQAXAMe599G5ermXi87m3pkpIiIxMbGyowZM+71Fjpx4oRERmazQTI1WeTwbJFPK5h6D20fJRKXvx8PWIpOp5OhQ4cKIN26dTO78VhEZOc1U0+iTps7yfXo65YLMhN3796VixcvPrJM2Nx5eXahm+wil3sZrQduAamY7gyGpB3vBFzG1AYw3pxr5WRTCaFwS01NlUWLFom3t7cAsn379uxfzGgU+WOzyFc1TO0E3/YRufPoD4rCLCIiQlq1aiWAjB071uxV4u7vSfTyzy9bvSeRyL9rWqSkpGRaNnrHDrlQuYrc/GCs1ds2LClXE0Je2VRCKJyMRqP89NNPUq1aNQHk2WefNWtJwwwFHRdZ0saUCBY0EQncm3vBFlAHDx4UNzc3Wb16tdnn6PQ6GX9o/L2eRCn6zD+Qc5tOp5NOnTqJnZ2dHDhw4JFlE8+elYs1a8n1gQPFYEbyyE/MTQhq6golX5g8eTI6nY7NmzfTs2fP7E1JfWEb3A4wtRG4lYFu86DOANOaBUq6bty4ga+vL82aNePGjRvpzv6ZnpiUGEbvH82J2yd4q/ZbvFX7LavMSXQ/EeH111/n559/ZtGiRRlOSQGQGhpK8LDhOHh7U37uXOycnKwYaR5iTtbIK5u6Qyg8goODZejQoRIWFiYiIkFBQWbd8qfr0m7TgLKPPEzbD8NEkuNyMdqCx2g0yuzZs8XBwUF27tyZpXODYoKky5YuUmd1HdkWaLtpHsaNGyeATJ48+ZHl9HHxplXP6jeQ5CtXrBSddWHmHYLqPqHkKXFxcUycOJGnnnqKFStW3FvYvEKFCjhl9VtbSjzs/xw2DACj3nRMswevx8HZLZcjLzhSU1N58803GTVqFF26dKFp06Zmn3s27CwDfx5IVEoUS9ouoatfVwtGmrHQ0FDmz5/P0KFDmTRpUoblxGAg9P33SQkM5LGvvsK5UiUrRpn3qEdGSp4gIixevJhJkyYRFhbGCy+8wKeffoqvr2/WL6bXwemVcHAGJNwFn2fh5ikw6E3rE/gWjNWtLCEiIoI+ffqwf/9+xo0bx9SpU7Ezs9vtzus7mXB4AmWKlmFBmwX4ePhYONqMlStXjtOnT/P4448/8lFV2KwviN+3D++JE3BrZn7iK6hUQlDyBE3T+O2336hcuTLbt2/n6aezMXW00Qh/boK9UyH6b9MH/wvfmaakDj5RqNcnMNeOHTv4/fffWbNmDS+++KJZ54gIS/5Ywlz/udQrXY85z82hmEsxC0eavn379nHmzBlGjx5NpUy+7Udv2kTkihUUHziQEgMHWinCPM6c50p5ZVNtCAXLmTNnpF27dnL+/HkRMXUPzFZXP6NR5NIuU4+hjzxEFj4rcvlXNSV1Ftw/8d+1a+bPOKrT6+TDQx9KjZU15IODH9ikJ9E//P39xd3dXapXr/7IKSlEROKPHZcL1WvI30NeE2NqqpUitB1UG4KSVwUHB/PKK69Qv359Tp8+zfXr1wEoWrRo1nuiBB2HFZ1gXV/QxUPvZfDGQXiyTaFdwzgrRITZs2fzxBNPcP78eQAef9y8tYtjUmIY+ttQtl3dxv9q/4/pTafjZG+b3jnXr1+nY8eOFCtWjF9++eWRS6Hqbtwg5O23cfLx4bGvvkRzUA9K7jEna+SVTd0h5H8ff/yxuLi4iLOzs4wZM0aio6Ozd6E7F0TW9TfdEcyoJHJ8sUhqweo7bmkpKSny2muvCSA9e/aUuDjze17turZLmq5vKrVX1bZpTyIRkbCwMHnyySelePHicuHChUeW1UdHS2D7DnKpUWNJCQqyUoS2hxqHoOQVBoPh3nq68fHx9OrVi08//RQfn2w0OkYHwb7pELAenN2h1QRo/D9wKprLURds4eHh9OnThwMHDjB+/Hg+/vhjsxuP111cx/QT0wFwsnOignuFTM6wrP379xMaGsru3bupWrVqhuUkNZWQkaPQ3byJz4rlOFWwbdx5kUoIisWICDt27GDMmDHMmTOHdu3a8dlnn2VvgFJCBByaBSeXAho8MwyavQtF0p+gTHm0OXPmcOzYMb799lsGmtmgKiKsvrCaL059ce+YQQycunOKOqXrWCrUTPXt25cWLVpQunTGM5KKCLenTiPx2DHKTp9OkQYNrBhh/qHaEBSLOHPmDK1bt6Zbt24YjUYc0p7TZjkZ/DOWYE5tOP6NafbREaeh/TSVDLIhKSkJgEmTJnHixAmzk0FiaiIfHPyAWadmUd+7Ps72zthr9jjaOdLA2/ofrkajkbfeeotffvkF4JHJACBqzRqiN2zA6/XXKdazhzVCzJ/Mea6UVzbVhpA/vPPOOwJIyZIlZd68eaLT6bJ+kdQUkWOLRGb4/btamZp8LtuMRqN8+eWX4ufnd2/0t7mCYoOk19ZeUnNlTVlybokYjUbxv+MvS84tEf87/haK+NHGjBkjgEydOjXTsnH798uFqtUkePjwfLsEZk6hJrdTrCkmJkZS07rvLV26VMaOHZu9BmODQSRgg8hXNU2JYHknkaATuRxt4ZKSkiJDhgwRQHr37i3x8fFmn3s45LA0WddEmqxrIodCDlkwSvN99dVXAshbb72VaTflpL8uyV/16svVnj3FkElX1IJMJQTFKnQ6ncyfP19KlSol33zzTfYv9OBYggVqLEFuuHv3rjRv3lwAmThxYpamrV5ybonUXFlTem3tJUGxeaNHzvr16wWQXr16Zbp+c2p4uFx5rpVcbtpMdLdvWynCvMnchKAalZVsERG2b9/OmDFjuHTpEi1btqRhw4bZu1jQcfhtMgQdheK+prEE1XuplcpywTvvvMPx48dZt24dL7zwglnnJKQmMOHwBH4L+o2Ovh2Z3GQyRRwz7tdvTQcPHqR58+asXbv2Xs+19BhTUggZNhx9ZCQ+a9bg6O1txSjzMXOyRm5swBPAMmDTfceqAt8Am4C3MrtGdu8Qzm76XH4e2UbObvo8W+crD3vjjTcEkMqVK8u2bduyN8I4YIPI3AZqLIEF/PPtOSwsTE6cMP+R27Xoa9Lth25Sa1UtWfnnyjyzSMw/cRiNxkxHIRuNRgl59z25ULmKxOz8xRrh5Xnk8oppy4Ew7ltCM+14B+ASEAiMNfNam9I5Zgcsy+zc7CSEs5s+l7PVq8iflatIQLUqsnf6QIm/c1U9isiGGzdu3GsX+O2332T+/PnZazAWEfltyr/TUU8pIRK4L/cCLcSMRqPMmjVLWrVqleXpwvcF7ZPGaxtLs/XN5FhoDhYgymWBgYHSuHFjuXTpklnl7y5YYFoCc+FCC0eWf5ibEMx9ZLQSmAes/ueApmn2wHygLaZlNU9qmrYNsAemP3D+YBEJS+/CmqZ1A8amXT/XhR76lYp6U8axM0CZlacJWtmZmKIQ56mhK+6IVtId13KlKeZbiTKVG1K6alPsPcqoqQ/SxMTEMH36dGbPns0777zD9OnTad26Na1bt87wnIC9Gwk9/Cvlmraldqt+//5CBA7MgENf/PdY6Gnwa2m5P6IQ0Ol0vPnmm6xYsYLevXuj1+vNmjLcKEa+CfiGhQELqeZVja9afkU5t3JWiDhzYWFhtG/fnujoaIxGY6blY3/5hbtzvsajW1e8hg61QoQFi1kJQUQOaprm+8Dhp4FAEbkGoGnad0B3EZkOdDE3ABHZBmzTNO0nYJ2555mrXLO2pP62AgcDGOzgZqtyaA72GO5E4RCRhGewjuIXIrCTCOAi0Wznrj1EeUKip4axhDOOpTxwK1eOkk9UoVy1JriVrwFu3gX+GXdqaiqLFi1iypQphIeH8/LLL/O///0v0/OOrP8S96lLqGAE/cbDBHyNKSkY9PDTaDizCiq1Nc0+akhVU1Lngrt379K7d28OHTrExIkTmTx5slkjj2N1sXx46EMOhBygm183JjaeiIuDixUizlxcXBydOnUiNDSUvXv3UqVKlQzLJvr7E7N9O9Hfb8K1bl3KTp1q9RXaCoKcNCo/BgTftx8CNMqosKZpXsA0oK6maeNEZLqmaS2BXoAz8HMG570BvAFQsWLFLAdZu/cYAoDgQ79SrllbOvUe81CZlKQEQv86Stj534m99hcpobfQwmJwiUyh2J/JFElJBsIwcJZgviPOFWI9BZ2nHXi54FqqOJ4VKuD9RHXKVH4ae68nwKNcvl+acfjw4SxevJjnnnuOWbNmUa9evYfKGPSpBJ7ZS/Dh3aScDaD4pduUiDHcVwCCD/9K7aad4ftBcGWXaYRxq4kQclJNSZ1Lnn/+eU6ePMn69evp37+/WecERgUyav8obsbd5MNGH9K/cv888yGq0+no3bs3Z8+eZevWrTRu3DjDson+/gS98iqi04Gm4fXm0MK7BGYOaabHS2YUNN0h7BCRGmn7fYH2IvJa2v5LwNMiMsIyoUKDBg3k1KlTlrp8ukSEqLvB3LzwO5F/nSbhRiCG0Ds4hMdTNCqV4jGCw313sno7iPYQEj0Evac9Dl5FcPP2wquCD+WfrINbxRrgWQE8y4O9o1X/FnOcPHmS0qVL4+Pjw8WLF7l27RqdOnW690GRFB/NX4e3c+fYATh3kVJXIymSYjo32t2OyKe8kZLF8fn1AvZGMGqgnzWaOjfXwa2z0GkWNBxiw7+wYBERNE3jjz/+ICkpyex1JHbf2M2EIxMo4lCEL1p+QX3v+haONGuio6Pp0qULQ4YMYdCgQRmW0wUFETJyJCkX/zIdsLOj1MiRlBz6hpUizR80TTstIpkOKc9JQngGmCwi7dP2xwGkPTKyCFskhMzoUpIIvf4HYZfOEHPlHClBN+B2BM7hiXhG6XFP+m/5BBeI8RBSPIxonva4lCyKp3dpvH0q4e1XC4dST0AxH1PCcHC22t9x48YNxo8fz7p163jttddYsmQJAJG3bnDpwFaiThzF6fxVvIMT7iXA295OxFetgFv9Bvg170L5J+vde0wRsHcjN+fP4fHzkVyta0eX6hHQZzlU6WS1v6kgExFmzZpFUFAQc+fONfs8g9HAXP+5LPtzGbVK1eLLFl/iXTRvdcnU6/U4ODj8Z1LEB6XeuUP4goVEb95senRrMIAImpMTFVcsp0jdulaOOm+zRkJwAC4DrYGbwElggIicz2bMmcqLCeFRRIToyFuEXvYnPPAcCVcvog8JweFONEUjUigeY8TxvqcrRg2i3YUED0HvYcTB04GiJd3wKluWcr5P4l6+KlpxHyhW0XSX4ZTzvuHR0dFMnz6dOXPmoGkar73cnw5VSqNd+BO3iyGUCtcBoLOH2z5upFb3o8TTTajcvDte3pnMVhpyih3vvYjfWY2g5+vTfsq3OY5XgZSUFN58801WrlxJ3759Wbdu3b25oh4lJiWGMQfHcDT0KH2f6svYp8fabP2CjHzxxRfs2LGD7du34+b28LrX+shIIhYvIWrdOkSE4n374vXmUFJv3iTxxEmKPN1QJYN05GpC0DRtPdASKAncAT4SkWWapnUCZmPqWbRcRKblKOpM5LeEkBldajK3bpznTuA5Yq7+RfKNQCT0Ds534/CM1OGZ8N/ySU4Q4ymkuBsRdyMunvZ4lvSgdLnyePs+hWMpPyhW4d+E4eLx78npLCGpT07m9Vf6s+r7rbT28eZtNw8qpaZNU+2qcadSCexqV6NM45ZUbdIF1yL3XU/EtHC9PgUMuvu2VNOxC9vg0Ez0riXZdaoUvgERRI4fRLOXHm7DUcwXFhZGr169OHLkCB999BGTJk0yq/H4UuQlRu4bSVhiGOMbjaf3U72tEG3WrF27lhdffJG+ffuyfv36/9wdGOLiiFyxgsiVqzAmJ+PZvTslhw3DqfxjNow4/8j1O4S8oKAlhMxERd0mNNCfiMALxN8IRB/0N/a3wikSkUiJKANO+n/LGoGYf+4u3I3Yuxko6mFPiVLFKFeqGDeuXuBmmAsODkYO33SgZLzQPN6RhGQ9t/V6SpVxJbqCA87lXfAp60wlTw0Hg/6BD3udaQH7f37GjPeOgwuJfddzeOQYvIMTkNmTqNPWvBGzyn8ZDAZq167N1atXWbVqFf369cv8JOCnaz8x+ehkPJw9+KrlV9QqVcvCkWbd7t276dy5M02bNmXnzp24uJh6OhkTE4lcu5aIpcswxsTg3qEDpd4egfMTT9g44vxFJYQCTqdP4VbwRe5cOUfM9Usk/f03EnIT5zvReESmUDzuv/9fBTiXlMTMu2GcSUqidWkP3m71GO5l7XiyjBOPebii2TuDvbOpsdvB2dQd9J/Nwem/+/859sA5l3fBn1sAI2j20Go8kU/24c8+XSkSr6fYqoVUqt3CJvWW3+3YsQNvb2+zpgnRG/V8efpL1lxYQ73S9fii5ReUdC1phSiz5tSpU7Rs2ZJKlSpx4MABPD09Mep0RG/8nvBvvsEQHk7RFs0pPXIkLtWq2TrcfMnchKDmMsqnhQfumgAAFlBJREFUnByc8Xm8Dj6Pp78wSXTMHUKvBhAeeJ7Ly1ax5fdr7IyLw8venuFt6/DVzyfNeu6cLcV94a+fTHcRaWMMSpTxwWfJUm6/+Cq33xyGx8bNlK5Q2TKv/wCjGDn5yyruHPgVrybNqdayF07ORXBycMbBzsEiXS0zHJiXRSLCzJkz8fT0ZOjQoXTpYt4Qn4ikCN4/+D4nb59kYNWBvNvgXRzt8l6vNgBXV1fq1avHhg0b8ChalOjNW7g7fx760FsUadCAUl/PoUg6XZ6V3KfuEAqB0W++xMLFa3m5RAleLlUCt/lTc/QhZZZ02iwAzh/eSupbY4ko7UKjLb/i5mneN1YRITEplpjwm8SG3yIh6g5JkeGkREegi45EHxODMTYW4hKwi0/CISEZpwQdLkkGiiYYcTY8fE29nWmwot4eDPYahn/+62CH0V7DaG+H0cG0ib094mCHONjf23BwMG2ODmiODuDoiObggERE4XPkOnZGSHWAoKmDqfVcX0q6lqSoY1GzE1BKSgpvvPEGq1evZuDAgaxZs8asc8+Hn2fU/lFEJUfx0TMf0dWvq1mvZ23x8fEULWqqD6PBQPzu3dz9ei6669dxqVGDUqNGUfTZJnlmbER+ph4ZFWI6nY6FCxdSvnx5evfuTXJyMge3rsBw6UyOv7HmhmOb5uM+cR6hFYuQUq4kWjlvHEuUQB+d9qEeG49dfCL2aR/qzkl6XJOMuKQ++rrJjpBcxIGUIg7oizpjKOqKuBfF/mYY5a7HYYeprSWkSgnsalVFUlNNm14PqalIqh70aVuqAU2vRzMY0fQG7FINaT8bsTcY0QxG7PWCvcGInUGwN4CDQXAwgN19/6QEuFwOFne0J7i0hrO9MyVdS+Ll6oWXi9e9n0u6lPz3uKsXxjgjA/oO4OjRo0yePJlJkyaZ9cH4w5UfmHpsKl6uXsx+bjbVvPLmI5bY2FhatmxJmzZtmNipE3dnzyHlr79wfrISJd9+G/c2bVQiyEUqIRRCIsKWLVv44IMPuHr1Kq+++iorVqywdVjp2vHe8/jtOPfQ8WRHSHK1J6WII6lFnTC4uSBuRdA83LHz8MDRsxhOxUrgWrwkRUqUxs2rLB5eZXEv4Y2dc/rjNgL2bkTe/ggHg+luQPt6ikWT4tnfvoNRU3BIuyvR7OywMxiJq/wY11r4cb6WJ3eM0YQnhxORFEFUchRyXwO9IdlA4IRADLEGao+oTc3WNR9OIK5pCcTFlEDO3T3Hl6e/5I/wP2hUthEzm8+kuEtxi/2NOZGSkkLnzp05sH8/S559lkZ3wnCsUIFSI4bj0bkz2iOmtVayR7UhFDKnT59m5MiRHDlyhOrVq7Nz507at29v67AyZO/hgUEDewGDBn93q0e7KUtxdHHN9deq3aofAV+bptCwxh1SnTb9Cfja7t7rVa/blpitW3Ha+D3uiw9S190dz65dKdZvNC5VqqA36olKjiI8KZyI5AjCk8LZcH0DJauWxNXXlfDkcK5GX+V40nFidbGPfG17zZ63ar+VZ5OB0WjkpR492LNnD9PLlOVZNEpOmUKxXj3RHPNmG0dhou4QCohNmzYxfPhwPvnkEwYNGmS5BuNcYu1v7XmBiJB0+jRRGzcS98suRKfDpXYtivfrh3uHDnwxfz6NGjWiRYuMe2DpDDoikiLuJY7wpHB23djFsVvHAFNCGF53OK/VfM1af5bZki9f5u1+/VgSEMB7FSsydsoUivfvj51L3phMryBTj4wKuKioKKZNm0bp0qUZM2aMqdE1MZGiRYvaOjSz5VZPnPzIEB1NzLZtRG3YSNyVK0yOCGdrZCRDX3iBb9ZlbdLfs2FneX3366QaU3G0c2RJuyXUKZ1+7zNb0P39N3fnzSd2xw72per4y8+PLzdtxsH94ZHIimWohFBApaSksGDBAj755BOio6MZMWIEc+bMsXVYSjbdvn2bHu3bc/zcOd72LsNQT09ca9WiWL++eHbqhJ2ZCf5s2FlO3TlFA+8GeSYZpN66ZZpvaMsWwoHKgwfjNWQw9sWK2Tq0QkclhAJo//79DBkyhGvXrtGuXTtmzpxJrVp5b9SpYp7Q0FCeeeYZwsPDWb16NT1atyZm23aiv99IypVA7IoUwaNrV4r164tr9eq2Dtds+ogIIhYvJmr9d4gI/vXqMmTjRjZv3kynTmpyQ1tQjcoFyD+zPrq5ueHm5sYvv/ySpxuMFfOUKVOGbt26MWjQoHtrTZR4+SWKv/QiSf5nid64kZgffyR6wwZcqlenWL9+eHTujL1b3nwsaIiNJWL5ciJXr0GSk/Hs2YMbzzzD688/T9WqVWnatKmtQ1Qyoe4Q8rCrV68ybtw4PDw8WLp0KfDv/PdK/iQizJ49m549e+Lr65tpeUNsrOmuYeNGUi5fRitSBM/OnSnWrx+uNWtYPmAzGBMTiVzzLRHLlmGMjcWjU0dKDh/B3/pUnn32Wdzd3Tl69ChlypSxdaiFlrpDyMciIyOZOnUq8+bNw9HRkbFjx95LBCoZ5F/Jycm8/vrrfPvtt0RGRvLJJ59keo69hwclXhxI8YEDSA4IIGrj92lLRX6PS7VqpruGLp1JuXLFqtM/J/r7k/D7MQzR0cT+/DOG8HDcWrak1Mi3calalejoaNrXrYumaezatUslg3xC3SHkMbt27aJ///7ExsYyePBgPv74Y8qWLWvrsJQcun37Nj179uTYsWNMnTqVDz/8MNvJ3RAXZ0oKGzaScukSODubRlcbjWBvT4lBr+JUseLD13/o9TL5/UP7pv/o/v6biCVLTa8JOFetSplJE/+TiESEzz77jDZt2pg1EZ9iWapROR8REaKjoylevDghISEMGzaMqVOnUrNmTVuHpuSCy5cv06ZNGyIiIlizZg29evXKleuKCMl//MGtjz8m5U+LrUv1aHZ2lBr5NiWHDgVMveBCQkLw8/OzTTxKutQjo3ziyJEjvPfee7i6urJnzx7Kly/P1q1bbR2WkovKlStH7dq1+fjjj6mbi49zNE3DtVYtyowfT9Crg5DUVDRHR8p9+QWuD04T/eAXvwf2H/5emHH55AsXCB3zAaLXozk6UiRtHWej0cjLL7/Mnj17uHz5MiVKlMjBX6fYgtUSgqZpTwDjAU8R6ZN2rCXwCXAe+E5E9lsrHlsLDAxk7NixbN68mbJlyzJ16lRbh6TkIhFhyZIlDBgwADc3N7Zv326x1ypSty4VV66wWhuCU4UKOJQu/Z/XExFGjRrFxo0bmTlzpkoG+ZWIZLoBy4Ew4M8HjncALgGBwFgzr7Xpvp9bADuBlUClzM6tX7++FAQ7d+4UR0dHKVq0qEyZMkXi4+NtHZKSi5KSkmTgwIECyJw5c2wdjlV8+umnAsg777wjRqPR1uEoDwBOiTmfz2YVguZAvfsTAqZ1lK8CTwBOQABQDagJ7HhgK33fefcnBLu0/3oDazOLIz8nhKSkJLl06ZKIiMTFxck777wjoaGhNo5KyW23bt2SRo0aCSDTpk0rFB+O27ZtE0AGDBggBoPB1uEo6cjVhGC6Hr4PJIRngF337Y8DxplxnU3pHHNK73ja794ATgGnKlasaLkasxCj0Sjr168XX19fefLJJyU1NdXWISkWEhAQIOXLl5ciRYrIli1bbB2O1SQlJcnUqVMlJSXF1qEoGTA3Idjl4GnTY0DwffshacfSpWmal6Zp3wB1NU0bl3asl6Zpi4A1wLz0zhORxSLSQEQalCpVKgfhWt/hw4dp3LgxL7zwAp6enixYsCDPz0KqZJ+Hhwfe3t4cOXKEnj172jociwsICCA6OhoXFxfGjx+Pk5OTrUNScignn07pdaLOsA+riEQAbz5wbAuwJQcx5Fl79+6ldevWlCtXjhUrVvDSSy9hrxb+KHBEhE2bNtG7d298fX05efJkoRg8eOnSJVq3bk3Tpk358ccfbR2OkktycocQAlS4b788EJqzcPK38PBw9u3bB0CLFi2YP38+ly9f5tVXX1XJoABKSkpi4MCB9OvXjy1bTN9rCkMyCA0NpX379tjb2/PFF1/YOhwlF+XkDuEk8KSmaY8DN4H+wIBciSqfSU5OZu7cuUybNg1HR0eCg4NxcXHhf//7n61DUyzk1q1bdO/enVOnTjF9+nR69+5t65CsIiYmho4dOxIREcH+/fvVALQCxqw7BE3T1gO/A5U1TQvRNG2IiOiB4cAu4CKwUURsNFzSNoxGI+vXr6dKlSqMGTOGpk2bsn//flzUClAF2pkzZ2jYsCEXLlxgy5YtjB07tlDcGQAMHz6cixcvsmXLFurXr2/rcJRcZtYdgoi8kMHxn4GfczWifOTUqVMMGDCAOnXqsGzZMlq3bm3rkBQriI+Px9XVlZ9++onatWvbOhyr+uyzz+jXrx9t27a1dSiKBai5jLLo8uXLHDlyhEGDBgHw22+/0apVK+zsctIco+R1IsKhQ4do3rw5AKmpqTgWkkXhRYTvv/+e3r17q7awfMrcuYzUp5iZwsPDGTFiBNWrV+e9994jLi4OgDZt2qhkUMAlJSUxYMAAWrRowfHjxwEKTTIAmD59Os8//zxr1661dSiKhalPskwkJyczY8YM/Pz8WLhwIa+99hoXLlzA3d3d1qEpVhAaGkqLFi3YsGEDn3/+OU+nTeRWWCxfvpzx48fz4osv8uKLL9o6HMXC1CipTISEhDBhwgTat2/P559/TrUHZ5FUCqzTp0/TrVs3YmJi+PHHH+nWrZutQ7KqHTt28MYbb9C+fXuWL1+u7oQLAZUQ0nHgwAF27tzJZ599RqVKlbh48aLqXlcInT59GkdHR44ePUqtWrVsHY5VxcfH8+qrr1K3bl02bdpUqB6RFWaqUfk+ly5dYsyYMWzbto3y5cvj7+9PyZIlLfZ6St4jIly8ePHenWBcXFyhfTx44sQJfH19KV26tK1DUXJINSpnQVRUFMOHD6d69ers27ePTz/9lMuXL6tkUMgkJibSv39/GjZsyI0bNwAKXTK4efMma9asAeDpp59WyaCQUY+M0nz//fcMHTqUjz76SP0jKIRu3rxJjx49OH36NJ9//jk+Pj62DsnqoqKi6NChA3///Tft2rXD29vb1iEpVlYoE4LRaGTdunVs2LCBH3/8keLFi3P16lXc3NxsHZpiA6dOnaJ79+7ExsaydetWunbtauuQrC4pKYnu3btz6dIlfv75Z5UMCqlC98ho//79NGzYkJdeeolbt24RFhYGoJJBIbZ27VqcnJw4evRooUwGBoOBgQMHcujQIVavXk2bNm1sHZJiI4UmIYSHh9OtWzeee+457t69y5o1azhx4gRly5a1dWiKDRiNRkJDTZPzzpgxg5MnT1KzZk0bR2Ubv/76Kz/88AOzZ8+mf//+tg5HsaFC88jI09OT27dvM336dEaOHImrq6utQ1JsJDExkVdeeYUTJ04QEBBAsWLFCnUHgg4dOnDq1Ck1WZ1SeBKCo6Mjx48fLzSzUirpCwkJoXv37vj7+zNjxgw8PT1tHZLNrFq1Cj8/P5o2baqSgQIUooQAhWPxEiVjJ06coEePHsTFxbFt2za6dOli65BsZuvWrQwePJgePXrQtGlTW4ej5BGFKiEohdvkyZNxdnZm9+7d1KhRw9bh2MyRI0fo378/DRo0YPXq1bYOR8lDVEJQCjSj0UhCQgLu7u6sWbMGo9FIqVKlbB2WzZw/f54uXbpQsWJFfvrpJ4oWLWrrkJQ8xGoJQdO0J4DxgKeI9Ek71gwYmBZHNRFpYq14lIIvISGBV155hbCwMPbs2YOXl5etQ7K5xYsX4+Liwq5duwp1Q7qSPnOX0FyuaVqYpml/PnC8g6ZplzRNC9Q0beyjriEi10RkyAPHDonIm8AOYFVWg1eUjISEhNCsWTN++OEHevTogYODuhkG+Oqrrzh+/Di+vr62DkXJg8wdh7AS6HD/AU3T7IH5QEegGvCCpmnVNE2rqWnajge2zOaCGACsz2LsipKu48eP07BhQwIDA9m+fTujR48u1B0KkpKSGDx4MEFBQdjZ2VGxYkVbh6TkUWbPdqppmi+wQ0RqpO0/A0wWkfZp++MARGR6JtfZ9M8jo7T9isBEEXk9g/JvAG+k7VYGbgMx6RT1fOD4g/slgfBHxZbLHnx9S59vTvlHlcnod+YeT6+cNes8p/Wd1WtYqr4z+p05xwryezyn9f2o3xeG+vYRkcwbz0TErA3wBf68b78PsPS+/ZeAeY843wv4BrgKjLvv+BSgSRbiWGzO8XT2T5n7GrmxZRSnpc43p/yjyphbr+bWt7XrPKf1ndVrWKq+H1GXmR4ryO/xnNb3o35fWOrbnC0nD1bTuwfP8HZDRCKAN9M5/lEWX3e7mcczKmctOX39rJ5vTvlHlTG3XjM6nt/rO6vXsFR9Z/Q7c49ZkzXf4zmt70f9vrDUd6as/sjIVjRNOyVmLBCh5B5V59al6tu6CmJ952Ryu5PAk5qmPa5pmhPQH9iWO2FZxGJbB1AIqTq3LlXf1lXg6tusOwRN09YDLTE1otwBPhKRZZqmdQJmA/bAchGZZsFYFUVRFAvKV2sqK4qiKJZTaNZDUBRFUR5NJQRFURQFUAkBAE3TemiatkTTtK2aprWzdTwFnaZpT2iatkzTtE22jqWg0jStqKZpq9Le1wNtHU9hUBDe1/k+IeTSPEs/immk9KvA8xYMN9+z1LxWSuayWPe9gE1p7+tuVg+2gMhKnReE93W+Twjk7jxLE9LOUzK2EsvOa6VkbCVm1j1QHghOK2awYowFzUrMr/N8L99PASkiB9MGzd3vaSBQRK4BaJr2HdA9bdDcQ8tkaaaZzz4DdorIGctGnL/lRn0r2ZOVugdCMCWFsxSML342kcU6v2Dd6HJfQX2jPMa/347A9I/jsUeUHwG0AfpomvbQ9BpKprJU35qmeWma9g1Q958R7kq2ZVT3W4DemqYtxPZTLhQ06dZ5QXhf5/s7hAxkdZ6lr4GvLRdOgZcr81op2ZJu3YtIAjDI2sEUEhnVeb5/XxfUO4QQoMJ9++WBUBvFUhio+rYdVffWV2DrvKAmhPw2z1J+p+rbdlTdW1+BrfN8nxDS5ln6HaisaVqIpmlDREQPDAd2AReBjSJy3pZxFhSqvm1H1b31FbY6V3MZKYqiKEABuENQFEVRcodKCIqiKAqgEoKiKIqSRiUERVEUBVAJQVEURUmjEoKiKIoCqISgKIqipFEJQVEURQFUQlCUbNE0raemaaJpWhVbx6IouUUlBEXJnheAU5jmsVGUAkElBEXJIk3T3IAWwBBMieGf4/s0TWub9vNUTdPUlOpKvlJQ10NQFEvqAfwmIuc0TUvQNK1e2kp7HwEfpy0TWhe1lrGSz6g7BEXJuheAjWk/b0zbR0QOYlo8ZTTQX0TUWsZKvqISgqJkgaZpXpjW1P0l7dAG4HnNpCZQFkgRkThbxago2aUSgqJkTR/gZxFJARCR68BtoCmwFtNi6wmaprW3XYiKkj1qPQRFyQJN0/YDtYDY+w4/BgQCb4vIr5qmNQc+F5FnbBCiomSbSgiKoigKoB4ZKYqiKGlUQlAURVEAlRAURVGUNCohKIqiKIBKCIqiKEoalRAURVEUQCUERVEUJY1KCIqiKAoA/wddqjT0RLvTxAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"x0,dx = 30., numpy.array( [20,15,10,5,2,1,0.5,.25,.125,.1,.05,.01] )\n",
"n = len(dx)\n",
"y,yt,d = numpy.zeros((n,4)), numpy.zeros(n), numpy.zeros(n)\n",
"for i in range(n):\n",
" x = numpy.array( [x0, x0+dx[i]] )\n",
" for k in range(4):\n",
" y[i,k] = relerror(x,n=k+2)\n",
"plt.loglog(dx, y[:,0], '.-', label='n=2');\n",
"plt.loglog(dx, y[:,1], '.-', label='n=3');\n",
"plt.loglog(dx, y[:,2], '.-', label='n=4');\n",
"plt.loglog(dx, y[:,3], '.-', label='n=5');\n",
"plt.loglog(dx, 3*y[0,0]*(dx/dx[0])**(2), 'k--');\n",
"plt.loglog(dx, 3*y[0,1]*(dx/dx[0])**(4), 'k--');\n",
"plt.loglog(dx, 3*y[0,2]*(dx/dx[0])**(6), 'k--');\n",
"plt.loglog(dx, 3*y[0,3]*(dx/dx[0])**(8), 'k--');\n",
"plt.ylim(1e-17,1);\n",
"plt.legend(); plt.xlabel('$\\Delta x$'); plt.title('Relative error');"
]
}
],
"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.7.1"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment