Skip to content

Instantly share code, notes, and snippets.

@shane5ul
Last active October 19, 2020 20:29
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save shane5ul/06d000b066636dbaf111187bb4264a2f to your computer and use it in GitHub Desktop.
Save shane5ul/06d000b066636dbaf111187bb4264a2f to your computer and use it in GitHub Desktop.
Evaluating integrals of functions with smoothness properties in log-log space using trapezoidal rule
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"plt.style.use('myjournal')\n",
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Problem\n",
"\n",
"Consider the problem described in this [StackOverFlow post](https://scicomp.stackexchange.com/questions/20901/integral-in-log-log-space). You have a function with certain smoothness properties that are apparent on a log-log plot. This is often accompanied by a large domain of integration. It seems worthwhile to \"integrate in logspace\", whatever that means.\n",
"\n",
"For example, consider the following integral,\n",
"$$\\int_{x_\\min}^{x_\\max} \\dfrac{1}{1+ax^2} \\, dx = \\dfrac{\\tan^{-1}(\\sqrt{a} x_2) - \\tan^{-1}(\\sqrt{a} x_1)}{\\sqrt{a}}$$\n",
"\n",
"Suppose $a = 0.01$ and the limits of integration span 5 decades."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"fun1 = lambda x: 1.0/(1.0 + 0.01*x**2)\n",
"n = 100\n",
"xmin = 1e-2\n",
"xmax = 1e3\n",
"a = 0.01\n",
"xi = np.geomspace(xmin, xmax)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAvUAAAEoCAYAAADYCg7hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdd3xcZ5X/8c+ZGVXLkm1Zttx7EpcktmPSQ0wgkAC20xfCJtQUyi4su+zCD370LSz8gKUEksAuJMASEpI4DoSE9F4cp9lO4t6bLMuyLdlqc35/zIwsK5I90pQ7o/m+X695jXTn3jtHiZ97zzxznucxd0dERERERPJXKOgAREREREQkNUrqRURERETynJJ6EREREZE8p6ReRERERCTPKakXEREREclzSupFRERERPKcknoREUmKmZ1jZm/my3lFJLPM7Odm9n+DjkNilNRLRpnZBjN7VwbPP9/MtmTq/CKFqqe26+5PuPvx6TxnOs4rUsiCvM+6+/Xu/q1Mvbf0jZJ6EREREZE8p6Ress7MSszsh2a2Lf74oZmVdHn9n81se/y1T5iZm9nUfrzPdDN71Mz2mtkKM1vY5bVqM1tiZvvM7AUz+7aZPZmuv1FkIOreYxfvIfwnM3vVzBrN7DYzK82V84oUqizeZ39lZt+O/zzfzLaY2T+a2a74+T+azr9Ljk5JvQThy8DpwGzgZOBU4CsAZnYB8HngXcBUYH5/3sDMioAlwAPACODvgN+aWeIr/p8CTUAt8OH4Q0T67grgAmAScBLwkRw/r0ghyPh9the1QBUwBvg48FMzG5rG88tRKKmXIHwI+Ka773L3OuAbwFXx164A/sfdV7h7M/D1fr7H6UAF8B/u3uruDwP3Ah80szBwKfA1d29295XAr1P4e0QK2Y/cfZu77yH2QXp2jp9XpBBk4z7bk7b4+7a5+5+BA4DGy2SJknoJwmhgY5ffN8a3JV7b3OW1zp/NbLyZHUg8kniPze4e7fY+Y4AaINLb+4hIn+zo8nMzsQ/TmNl9Xdrrh9J1XhFJSjbusz2pd/f2Lr+r7WZRJOgApCBtAyYAK+K/j49vA9gOjO2y77jED+6+ieQvDtuAcWYW6pLYjwdWAXVAe/x9VnV/HxFJnbtfGHQMIgUsG/dZyTHqqZdsKDKz0sQD+F/gK2ZWY2bDga8Cv4nv+wfgo/FBruVAUvPfdj1//D2eJ9ZD8M9mVmRm84EFwO/dvQO4E/i6mZWb2QnA1en8g0UGiO5tNx0dQUXd2qs6l0RSl/X7rJlZRv4S6Tcl9ZINfwYOdnmUAkuBV4HXgGXAtwHc/T7gR8AjwBrg2fg5Wo5y/jHdzn+QWM/DAuBCYDdwA3C1u78RP+YzxAbz7ABuJXYBPNp7iBSi7m336zl6TpFCF8R9dkq6/whJjbl70DGI9MrMpgPLgZJudXrpfp/vALXurllwRESkYGTrPiuZp556yTlmdnF8jt2hwHeAJem+0JjZCWZ2ksWcSmzqrbvS+R4iIiK5KBv3Wck+JfWSi64DdgFrgQ7gkxl4j8HE6uqbgNuA/wcszsD7iIiI5Jps3Gcly1R+IyIiIiKS59RTLyIiIiKS55TUi4iIiIjkOc0P3AfDhw/3iRMnBh2GSM548cUXd7t7TdBx9EZtVuRIudxm1V5FjtTX9qqkvg8mTpzI0qVLgw5DJGeY2cZj7xUctVmRI+Vym1V7FTlSX9urym9ERERERPJcQSb1ZjbFzJ40s1Vm9pKZzQs6JhERERGR/irIpB74OfBrdz8O+Gfgt2ZmAcckMiCZ2Vgz+7GZPWNmzWbmZjYxyWNDZvYlM9tgZofM7BUzuzSzEYsULrVXkfyVF0l9Xy4yZjbOzO4ws0Yz22dmd5rZ+C6v1wCnA78CcPe/AgackvE/RKQwTQWuABqAJ/p47LeArwM/AS4EngVuN7P3pjNAEemk9iqSp/IiqSfJi4yZlQMPAycAHwauAqYBj5jZoPhu44Ht7t7W5dAN8e0ikn6Pu/tId38vcHuyB5nZCOCfgP9w9++5+yPufh3wCPAfGYpVpNCpvYrkqXxJ6pO9yFwDTAYucve73X0xsBCYQGxJZBHJMneP9vPQ9wDFwG+6bf8NcKKZTUopMBF5C7VXkfyVF1Na9uEisxB41t3XdDl2vZk9BSwCvg9sAkaZWVGX3vqJ8e0puX/FDjbsbkr1NGnRlxECxuGdezsuMeQgZLFaJTPDLPYcMggd8WxEwkY4ZERCISIhoygSojgcojhilETClBaFKS8OU1YUpqI0QlE4Xz5fShbNBFqANd22r4g/zwDW9/fk2xsPcuszG/nQ6RMYM6Ssv6cRkZiMtteOqPOLJ9b19/Ae9ec+2fWYxH3R4tuPuDeSuD/G742h+L0xlLg3xp4T98aicIiSSIjSoth9sbQ4RGVpESWRUOf7iBxLXiT1fTATWNzD9hXA5QDuXmdmzwMfAW42s/OJtb8XU33zO5dt4f4VO1M9TUEqLw5TWVrEkPIihleUMLKylNqqEmqryphaU8G0kRVUDyrWxa2wDAP2urt3276ny+v99uTq3dzw6Fo27mnmp1fOTeVUIpLh9toRdf79vjdSOUVeKgoblaVFVJUXUVNRQs3gEkYMLmXcsDIm11QwefggRg8pIxzSvVEGXlI/jFjdfXd7gKFdfr8e+LWZfQFoBj7Uw4UIADO7FrgWYPz4o5fdnz+jlgnVg466Tzb08qf0sm+Xn3t5zeOvJH6PuuMe2x712PZo1Im60xF/rT3qRKNOW0eU9vhza3uU1o4oLW1RDrV10NzaQXNrO02tiZ872LHvELC/x1iHlhcxY3Qlp0+q5owp1Zw0dgjFEfXwy5H60mYBtu89mOmQRKQXybbXkMG1b5+ctvftz33S6f2+mDhfNH5vdI/9DE5HNHavjEZj98j2qNPR4bRHo7R1HL4/trTH7o0H2zo42NrB/kPttHZEqW9qpb6plXV1PVcClBaFmDNuKG+bNIzTJg1j7vihlBWH+/FfRfLdQEvqk+Luq4Ezk9z3JuAmgHnz5h31KnDZKWNTD64AuTtNrR00Hmxjb3Mru/a3sLPxEDv3tbCloZk1dQdYs/MADc1tPLWmnqfW1MNfoawozDunj+Bv3jaOs6YMJ6SeioGmARhiZtbtQ3eix29PD8ck3WZnjakC4EBLe1qCFSlwGW2vkXCI//Pe6emKNW8caosl9w3NrdTtb6Fufws79x1iQ30z6+oOsH53E7v2t/DMunqeWVcPxL75vnDWKK6YN5ZTJw3TN9wFZKAl9Q0c2SOf0FsPflLMbAGwYOrUqf09hRyFmVFREqGiJMKYIWXM7GEfd2fHvkMs27iXZ+MXrzW7DnDvq9u599XtjBlSxmWnjOVjZ0+iqqwo63+DZMQKoASYwpF1ujPizytTOXl5vCfrYFtHKqcRkZiMttdCVVoUG4NWM7iE40YO7nGf+gMtvLChgRc27OG59fUs37qPPy7bwh+XbWH8sHKuO3cyH3jbeJXoFICBltSvgB5zwhmkcEFx9yXAknnz5l3T33NIasyMUVVlvO+kMt530igAtjQ0c+eyrfxh6Wa2NBzkvx5azS3PbOBz7zqOK08br8G3+e8vQBvwIeAbXbb/LbDc3fs96A5iN0uAQ239nexDRLrIaHuV3lVXlHDBrFoumFULwPrdTdzx4mb++OJWNu1p5st3Led/n9/ENxfNYu74nvo9ZaAYaFnPPcDpZtZZeBdfpOqs+Gv9YmYLzOymxsbGlAOU9Bk7tJy/f+c0Hv/CO/jdJ07jtEnDaGhu42v3rOCCHz7OE6vrgg5R4szsMjO7jMOLvF0Y33Zul33azeyXid/dfRexGau+ZGafN7P5ZvYz4DzgS6nGVBIfj9HarqRepKtcbK+SvEnDB/GF95zAU188j59cOYfRVaUs37qPS254mi/c/opKDgewvOmpj19g4MiLTB1Q5+6PxbfdDHwGWGxmXyE2ruVbwGbgxv6+t3rqc1soZJw5dThnTKnm/hU7+ff7XmdtXRNX/fJ5/uFdx/F3501VvX3wuq8vcUP8+TFgfvzncPzR1ZeBA8BngVrgTeAKd7831YASg6xb2lV+I9JNzrVX6btwyHj/SaM574QR/OThNdz8xDpuf3ELb+7cz68+eirDBhUHHaKkmfVlBHiQzKy3QB9z9/ld9hsP/ABITFX5EPA5d9+Qagzz5s3zpUuXpnoaybDW9ig/f2wtP3hwFe7w7hkj+f7fzKaiJG8+w+YNM3vR3ecFHUdvjtZm2zuiTP3yfYRDxtp/0yr2Uhhyuc3qHptZa3Yd4GO/eoFNe5qZUjOIWz9+GqO1RkdO62t7zZvyG3e3Xh7zu+23yd0vdfdKdx/s7helmtCr/Ca/FEdC/P07p/E/H3kblaURHli5k4t++hQb63NjYTDJDeFQbJGYjmhsyjkRkYFs6ogK7rj+DE6oHczauiYu//kzrKs7EHRYkkZ5k9QHyd2XuPu1VVVVQYcifTD/+BHc85mzOW5kBWt2HeDKm59jR+OhoMOSHGFmFMcHU7d1qK5eRAa+EZWl3HbtGcwdP4Stew/ygZueZfeBlqDDkjRRUi8D2sThg7jzU2d1XsCu/u/n2NvcGnRYkiOKlNSLSIGpKi/iN584jbdNHMqu/S188Y+v9mkxLsldSuplwKsoifDfH3kb00ZUsGpnrKbwYKsGR0psCXaAtg7d0ESkcJQXR/ivD8yhsjTCg6/v4nfPbwo6JEkDJfVJUE19/htSXswtHz+V0VWlLNu0l0/99kX1zgqReE99u/4tiEiBGT2kjH+9+EQAvnXvStaqvj7vKalPgmrqB4ZRVWXc8vHTGFpexCNv1vHjh1YHHZIErCg+1WmbBsqKSAFacPJoLp4zhkNtUT73+5e1bkeeU1IvBWXqiApu+NApmMFPHlnD0g17gg5JAqSeehEpdN9YNJMxQ8p4bWsjNzy6JuhwJAVK6pOg8puB5Ywp1Vz39ilEHT5328vsP9QWdEgSkIhq6kWkwFWWFvG9y08G4JdPrGef7ol5S0l9ElR+M/B8/vzjmDWmki0NB/na4hVBhyMBKQrFe+qj6qkXkcJ1xpRqzphczf6Wdm59ZmPQ4Ug/KamXglQcCfHDv5lDaVGIO1/aypJXtgUdkgQgHK+pb1dPvYgUuE+9YwoA//3kes0Ql6eU1EvBmjqigi+/bwYAX7l7OQ1Nmr++0CSmtGzXQFkRKXBnTx3OSWOrqG9q5Q9LNwcdjvSDknopaH972njOnFJN48E2fvDgqqDDkSxL9NR3qPxGRAqcmfGp+bHe+pseX6dpn/OQkvokaKDswGVmfHXBDEIGv3l2I2/s2Bd0SJJFkVBiRVn11IuIvHtGLVNqBrF170HueVllqflGSX0SNFB2YDuhtpIPnTaBqMM3l6zUctkFJDH7TYfKb0RECIWM68+N9db/7LG1RHVtzCtK6kWIzYZTVVbE02vreWDlzqDDkSxJlN/oa2YRkZiL5oxhzJAy1uw6wGOr6oIOR/pASb0IMHRQMZ8//zgA/vVPr3OoTSP/C0EkpJ56EZGuisIhPnjqOACWvKoSnHyipF4k7kOnjee4kRVs2tPMr57eEHQ4kgXheE29knoRkcMuPHEUAH9duZPWdn2TmS+U1IvERcKhzikub3xsLU0t7QFHJJmmnnoRkbeaUlPBCbWD2X+onafW7g46HEmSkvokaPabwvH2acM5ZcJQGprbuEWr6g14Yc1TLyLSowtnxXrr//zq9oAjkWQpqU+CZr8pHGbGZ985DYCbn1in3voBTj31IiI9e99JtQA8sHKnJhPIE0rqRbo5Z9pw5o4fwp6mVm59Vr31A1li9hv11IuIHGnqiMFMG1FB48E2nl5bH3Q4kgQl9SLdmBmffVdsJpybHl9Hc6t66weqiFaUFRHp1XvjA2bve00lOPlASb1ID94+bTizx8V761VbP2AlZr9RT72IyFslkvr7V+xQCU4eUFIv0gMz43PvitXWq7d+4ArHr4CqqRcReavjRlYwpWYQDc1tPLduT9DhyDEoqRfpxbnH1XDyuCHUN7Vyx4tbgg5HMiCS6KnvUFIvItKdmXX21v9JJTg5T0m9SC/MjOvfPhmAXz65Xr25A1BYs9+IiBxVYmrLB1bsIKprZU5TUi9yFO+eWcu4YWVsrG/mwdd3Bh2OpFnnQFnXjUpEpCfTRw1mVFUp9U2trN51IOhw5CiU1CdBi08VrnDI+NhZkwD4xRPrAo4mP5nZODO7w8wazWyfmd1pZuOTPHa8mf3azDaZ2UEzW2Vm3zazQemILaSeepG3yOU2K9lnZpw6aRgAz6/X1Ja5TEl9ErT4VGG7fN44BpdGeGFDAy9v3ht0OHnFzMqBh4ETgA8DVwHTgEeOdZOPv/4g8Hbg/wLvBX4B/CPw3+mIT4tPiRwp19usBCOR1D+3XoNlc1kk6ABEcl1FSYQrTxvPjY+t4xdPrOMnV84NOqR8cg0wGTje3dcAmNmrwGrgOuD7Rzn2LGLJxHvc/YH4tkfMbBjwT2ZW7u7NqQQXMi0+JdJNTrdZCcZpnT31e3B3LH7tlNyinnqRJHzkzIlEQsZ9y3ewpUH3pD5YCDybSA4A3H098BSw6BjHFsef93XbvpfYtSvlu4oWnxJ5i5xusxKMKTUVDBtUzK79LWys1z0wVympF0nCqKoy3n/SKDqizq+e2hB0OPlkJrC8h+0rgBnHOPZBYr2D3zGzGWZWYWbnAZ8Ffu7uTakGFw4nkvpUzyQyYOR0m5VgmBmnTjzcWy+5SUm9SJI+cU5sesvbXthMU4sWo0rSMKChh+17gKFHO9DdDwFnE7tOrQD2Aw8B9wKf6e04M7vWzJaa2dK6urqjBhc29dSLdJPVNtuX9irBUl197lNSL5KkWWOqmDdhKPtb2rn75a1BhzPgmVkpcBswgthgvXOBLwB/A/y0t+Pc/SZ3n+fu82pqao76HofnqU9T0CIFrD9tti/tVYLVOQPOBs2Ak6s0UFakD646YwJLNzZw6zMbufLU8RosdGwN9Ny711tvYFcfB+YDU919bXzb42bWCNxkZj9391dSCS6smnqR7nK6zUpwpo+qZHBJhM17DrJt70FGDykLOiTppmB76s3s/8bnz42a2UVBxyP54YJZtVQPKuaNHft5ceOx7m9C7Cv4mT1snwGsPMaxJwINXZKDhOfjz9NTjE2LT4m8VU63WQlOOGTMmxj7vPfCBpXg5KKCTeqBvwIXAI8HHYjkj5JImA+cOg6AW57ZGHA0eeEe4HQzm5zYYGYTiU19d88xjt0BDDWzqd22nxZ/TrkGKhyKXQI1T71Ip5xusxKsUydVA6qrz1U5m9Sb2Vgz+7GZPWNmzWbm8QtLT/v2efU7d3/W3bVEqPTZladNIGRw3/Lt1O1vCTqcXHczsAFYbGaLzGwhsBjYDNyY2MnMJphZu5l9tcuxvyI20O7PZvZhM3uHmX0B+B7wIrEp9lISjl8B2zuU1IvE5XSblWCdOkkz4OSynE3qganAFcRq+J7obadUVr8T6Y8xQ8p45/SRtHU4t72wKehwclp8CrvzgFXArcBvgfXAee5+oMuuBoTpck1y9w3A6cDLwLeBPxNbGOcm4Hx3T7kQvrOnXuU3IkDut1kJ1oljqigtCrFm1wF2H1CnVq7J5YGyj7v7SAAz+wTw7l72S2X1O5F+uer0Cfx15U5+99wmrj93CpFwLn8+Dpa7bwIuPcY+G+hhYRp3X0nsw31GJP63qfxG5LBcbrMSrOJIiLnjh/L02nqWbtjDBbNGBR2SdJGzmUgfPtEfc/U7M7vazF6OPz6d/mil0Jw9dTgTq8vZ1niIh97YFXQ40k+Jnvp2JfUiIknRfPW5K2eT+j445up37n6Lu8+OP3qd31okWaGQ8benTwDgd8+pBCdfJWa/iSqpFxFJypzxsRlwXtvSGHAk0t1ASOr7tfqdmX3dzLYAZwC/MLMtZja2h/202p306JK5YykOh3h8dR1bGpqDDkf6IRRfZ0A99SIiyZk5uhKA17fvU4dIjhkISX2/uPvX3X2su5e4+/D4z1t62E+r3UmPhg0q5oJZtbjDH5a+5Z+O5AH11IuI9M3wihJGVpbQ1NrBhvqmoMORLgZCUp/K6ndJMbMFZnZTY6O+apIjJeas/8MLm2nv0MQO+Saxoqx66kVEkjdzdBUAK7btCzgS6WogJPWprH6XFHdf4u7XVlVVpeN0MoCcMbmaidXl7Nh3iMdWqTwr3ySSes1+IyKSvEQJjpL63DIQkvpUVr8TSYmZ8YFTY+uc/e/zmwOORvoqoqReRKTPDif1qmDIJTmd1JvZZWZ2GXBKfNOF8W3ndtktqdXvUoxD5TfSq0vnjiUSMh5+Yyc7Gg8FHY70QUhJvYhInyXKb1Zu24dr8b6ckdNJPXB7/HF9/Pcb4r9/I7FDH1a/6zeV38jR1Awu4fwZI4k63L5UvfX5JNJZU6/xECIiyRo7tIzK0gj1Ta3s3KeVZXNFTif17m69POZ322+Tu1/q7pXuPtjdL4qvdieSFR+Ml+DctnSzZlLJI5019fpfJiKSNDNjhkpwck5OJ/W5QuU3cixnTx3O2KFlbGk4yJNrdgcdjiTp8EBZ9dSLiPSFZsDJPUrqk6DyGzmWUMi4Yl5sesvbX9Sc9fmic0pLddWLiPSJBsvmHiX1Imly6SljMYP7V+xgb3Nr0OFIEiKh2CUwqoFeIiJ90jlYdrt66nOFkvokqPxGkjFmSBlnTx1Oa3uUe17ZFnQ4koRw/AqoxadERPpmSs0gSiIhNu85SOPBtqDDEZTUJ0XlN5KsRAnOHzQLTl4Ix3vqNaWliEjfRMIhTqgdDMSmtpTgKakXSaPzZ4ykqqyI5Vv3qc4wD2jxKRGR/pvROVhW97tcoKReJI1Ki8JcNHs0ALcv1YDZXKfFp0RE+i8xWFY99blBSX0SVFMvfXF5vATn7pe30tLeEXA0cjSHF59SUi8i0leHZ8BRUp8LlNQnQTX10hezxlQxY1Qle5vbeHDlrqDDkaNITGmpBcNERPruhNpKQgZr6g5wqE2dWEFTUi+SAVfMGwvEVpiV3BW2WFLf1qHFp0RE+qqsOMyUmgo6os6bO/YHHU7BU1IvkgGLZo+hOBziydV17Gg8FHQ40otwON5Tr456EZF+mT4qVoLz5k4l9UFTUi+SAUMHFfOuGSOIOtz5kgbM5irNfiMikpppIyoAWLPrQMCRiJL6JGigrPTHZafESnDueHELrhVLc1JYSb2ISEqmjYzNVb9KPfWBU1KfBA2Ulf54+7QahleUsK6uiZc27w06HOlBJL74VHtUNfUiIv0xbWSsp371TvXUB01JvUiGRMIhLpk7Boj11kvuiXfUE3XNgCMi0h8ThpVTHA6xde9Bmlragw6noCmpF8mgS+fGSnCWvLJN033lIDM7XFevEikRkT6LhENMrhkEqK4+aErqRTLo+NrBnDS2iv2H2nlg5c6gwwmEmY0zszvMrNHM9pnZnWY2vg/HTzez281st5kdNLM3zeyz6YpPdfUiR8r1Niu5Z2p8sOxqJfWBUlIvkmFdB8wWGjMrBx4GTgA+DFwFTAMeMbNBSRw/D3gOKAE+AbwX+H9AOF0xalVZkcPyoc1K7pk2IjZYdvUuDZYNUiToAPKBmS0AFkydOjXoUCQPLThpNN++9/XOOetrq0qDDimbrgEmA8e7+xoAM3sVWA1cB3y/twPNLATcAjzk7hd3eemRdAbY2VPfoaRehDxos5J7EoNl12iwbKDUU58EzX4jqeg6Z/0flxVcb/1C4NlEcgDg7uuBp4BFxzh2PjCdoyQR6RAJawYckS5yvs1K7jkuntSvUk99oJTUi2RBogTnj8sKbs76mcDyHravAGYc49iz48+lZvasmbWZ2S4z+5GZlaUrQNXUixwh59us5J4J1YOIhIwtDQdpbtUMOEFRUi+SBedMq2F4RTHr6pp4ZUtBLWI2DGjoYfseYOgxjh0df74NeAA4H/hPYnW6v0tXgKqpFzlCzrdZyT1F4RCThg/CHdbVNQUdTsFSUi+SBUXhEItmx+as/2MBDpjtp8T16Tfu/lV3f9Tdvwd8A7jIzKb3dJCZXWtmS81saV1d3THfRD31ImnT5zbb1/Yqueu4kRosGzQl9SJZkpiz/p5XttHSXjBz1jfQc+9eb72BXdXHn//abfsD8ec5PR3k7je5+zx3n1dTU3PMAMPqqRfpKqtttq/tVXJXYlrLVRosGxgl9SJZMmN0JdNHVdJ4sI2HX98VdDjZsoJYjW53M4CVSRx7NGkZ2Xq4p14DZUXIgzYruSkxA85qJfWBUVIvkkWXzo2X4CzbGnAkWXMPcLqZTU5sMLOJwFnx147mPqAFeE+37RfEn5emI0DV1IscIefbrOSmxFz1a1R+E5h+JfVmNtPM/tHMbjWzp81shZm9YWbPmdnvzexfzOyUdAcrku8WzR5DOGQ8+uYu6g+0BB1ONtwMbAAWm9kiM1sILAY2AzcmdjKzCWbWbmZfTWxz93rg34HrzezfzOxdZvZF4KvAr7tOuZeKcCh2GVRNvQiQB21WctOk4YMIh4xNe5o51FYwJaY5Jemk3mI+ZGYriM1X+y5gJ/An4CfAD4g1/K3EprV6ML409CfNLK9XkjOzBWZ2U2NjQc1aIhlQM7iEc4+roT3qLH55W9DhZJy7NwHnAauAW4HfAuuB89y963e0RmzFye7XpG8C/wxcAfwZ+CTwXWIL5KRFRANlRTrlQ5uV3FQcCTGxupyow9o6leAEIakVZc1sErHGvQn4GPCCux+1Ns7MDJgH/B1wnZld5e6vpRhvINx9CbBk3rx5uihJyi6dO5aH39jFH5dt4WNnTwo6nIxz903ApcfYZwOxJKH7die2kE3GFrPRQFmRI+V6m5XcNW3EYNbWNbFm1wFmjtaCndl2zJ56M5tN7Cu3j7j7le7+3LESeog1bHd/wd2vBi4G/tPM5qccsUiee+f0EVSWRlixbR9v7NgXdDgFTz31IiLpcZwGywYqmfKbi4AFqdTCxZeYXgC8w8w0OFcKWmlRmAUnx9ZouatwBszmrM6e+g4l9SIiqZgan6t+1U4Nlg3CMRNsd/+6u6c8os/d2939a4hH0mYAACAASURBVMn08osMdJfE56y/66Wt6iEOWCSsnnoRkXSYFp+rfvUu9dQHod+95mY2ycz+wcyOtWy0iHQzd/wQJlaXs2t/C0+t2R10OAUtMftNm+apFxFJyaThgwgZbNrTXEiLLOaMVEphvgl8D/hSYkM80f+pmZ2WcmQiA5iZdfbW37lsS8DRFLYild+IiKRFaVGYccPK6Yg6G+ubgw6n4KSS1G8FzgF+lNgQr53/DHC+mb0zxdhEBrSL58QWovrLih0caGkPOJrCpRVlRUTSZ2pNrARnrUpwsi6VpH4vEHX3I7oZ47PefBtYmFJkGWRmQ83sXjNbZWavmNkDZjY16LiksIwbVs6pk4ZxqC3Kfa9tDzqcglUUjpffqKdeRCRlU+J19WuU1GddKkn9jcCvzewZM/uWmZ1nZqVdXi9OMbZMcuCH7n6cu58M3Av8IuCYpABdOjfWW3+nZsEJTGKgbLt66kVEUtbZU68FqLIulaT+ZuBpYktHfxR4EGgws+fM7HGg9GgHH4uZjTWzH8c/NDSbmZvZxF72HWdmd5hZo5ntM7M7zWx8b+d2973u/mCXTU8DPZ5bJJMuPHEUJZEQz6yrZ0uD6g+DEAmpp15EJF2mjBgEwBol9VmXSlK/0d0/6u5XuPtYYAbwBWA7MAT4VIqxTSW2zHQD8ERvO5lZOfAwcALwYeAqYBrwiJkNSvK9PgcsTilakX6oLC3i3TNrAbj7JfXWB0GLT4mIpM+Uzpr6JqK6rmZVKkn9EctDu/sb7v4Td78IuBz4eiqBAY+7+0h3fy9w+1H2uwaYDFzk7ne7+2Ji9fwTgOuO9SZm9rX48V861r4imXBJlxKc2Arrkk2d5TcdKr8REUnVkPJihlcUc7Ctg+37DgUdTkFJJan/jZn9l5mVdd1oZjOBk1MLC/qwSNVC4NmuK97GZ+F5ClgUj+lqM3s5/vh0l1i/ArwXuNDdVfsggThn6nCGV5SwbncTr2xpDDqcgqOBsiIi6TVFM+AEot9JvbsvA34MfNfMxnV56Wrgf4HqFGNL1kxgeQ/bVxArCcLdb3H32fHHT6Gzh34B8G537zWTMrNrzWypmS2tq6vLQPhS6CLhEItmjwY0Z30QEuU3GigrIpIeiRlwNFg2u1Lpqcfd17j7Z9x9c5fNXyNWfvMvKUWWvGHE6u672wP0uNpt/NuErxP74PFYvAd/aU/7uvtN7j7P3efV1NSkKWSRIyVKcJa8so3WdiWX2RTunP1GPfUiIumQmAFH01pmVyTdJ3T3Q8Cd6T5vOrn7CrqNCRAJ0oxRlRw/cjBv7tzPo2/u6hw8K5lXFJ/9RivKioikh3rqg3HUnnozi5jZR9L1Zmb22XSdq4sGeu6R760Hv8/MbIGZ3dTYqHpnyQwz6+ytv0uz4GSVBsqKiKTXlJr4tJa7mgKOpLAcNal393bggJn9sNvCUn1iZkPM7A7gjf6e4yhWEKur724GsDIdb+DuS9z92qqqqnScTqRHi2aPwQween0Xe5tbgw6nYHQOlFX5jYhIWoyuKqOsKMzuAy00NrcFHU7BOGZNvbvfAdwDPG5mf29mPdap98TMRpnZd4DHge+6+/39D7VX9wCnm9nkLu87ETgr/lrK1FMv2VBbVcrZU4fT2hHl3le3Bx1OwQhrnnoRkbQKhYzJNVqEKtuSGijr7g8D7wJGA2vis8H8yMyuMbPLzex8M3uPmX3QzD5tZj81s+XAq8BB4Ax3f66vwZnZZWZ2GXBKfNOF8W3ndtntZmADsNjMFpnZQmILSW0Gbuzre/ZEPfWSLRfPUQlOtiVmv2lT+Y2ISNpMVV191vVloKy5+xfN7JvA+4DzgWuBiUAV4MBeYD3wJLFVWh9391TqCLovOnVD/PkxYD6AuzeZ2XnAD4BbiQ2AfQj4nLvrX5LklffMrKW8eDkvbmxgw+4mJg5PdlFk6a9E+Y0GyoqIpI/mqs++pJJ6M/sZ8HEzu9DdHyKWbB9tlde0cPekZqhx903ApZmKw8wWAAumTp2aqbcQAWBQSYQLZtZy50tbueulrfzD+ccFHdKA1zlQVvPUi4ikjXrqsy/Zeer3EVtQqrOo3Mx+mpGIcpDKbySbLpk7FoiV4Lir9zjTOhefUk+9iEjaTNFc9VmXbFI/E3gEqO+ybVL6wxGRM6ZUM7KyhE17mnlxY1pmZZWjiCTmqddAWRGRtJk4vJyQwaY9zbS0dwQdTkFINqn/HPAVYoNkd5rZvcCE+MDUiZkKLldo9hvJpnDIuGh2bMDsnRowm3GJ8hsNlBURSZ+SSJjxw8qJOmzY3Rx0OAUh2dlv1gDTgYuIDUatjP9+F7DWzBrM7BEz+4GZfcTM5phZsh8Ycp7KbyTbEiU4976yTT0cGaaBsiIimaG6+uxKOvF297Z4cvtP7v524H5gCnA58CNidfcXA/8NLAV2m9ntZnaFmYUzELvIgHV87WBmjKpk36F2Hn59V9DhDGidNfUqvxERSSvV1WdXKr3pO919vbvf6e5fc/dF7j4RqCY23eW/AW3AvwKvmtkJqYcrUjgumasSnGzoXFFW5TciImk1Jd5Tv1pJfVb0O6l394/0sr3B3R929++5+5XuPg34APC1/r5X0FRTL0FYOHs0IYNH39zFnqZUlnsIlpmNM7M7zKzRzPaZ2Z1mNr4f5/mimbmZPZnO+JTUixwp19us5I/jRg4GYPXO/QFHUhgyXvduZvcT67XP24xYNfUShBGDSzlnWg1tHc69r24LOpx+MbNy4GHgBODDwFXANOARM0t6ZS0zm0xssH7aa5GKI7HLYGu7knqRfGizkj+mxXvq19U10a6Ok4zLxmDWh4mtOvt0Ft5LZEDpLMFZlrclONcAk4GL3P1ud18MLAQmANf14Tw/A34LvJ7uAIs0+41IVznfZiV/DCqJMGZIGa0dUTbUawacTMt4Uu/u33H3E939lky/l8hA8+4ZtVSURHh5817W5efsAQuBZ+MzaAHg7uuBp4BFyZzAzK4E5gJfykSAnT31mv1GBPKgzUp+OW5kvK5eJTgZN2CmnRQZiMqKw1wwqxaAu/NzwOxMYHkP21cAM451sJkNBX4A/LO770lzbAAUhxPlN5o6VIQ8aLOSX46rjdXVr9qZlx1TeUVJfRI0UFaCdMmcWAnOXS9vxT3vepOHAT0ti7sHGJrE8d8FVgG/SvYNzexaM1tqZkvr6uqOuX+ip75NPfUikOU229f2KvnnuBGJpF499ZmmpD4JGigrQTp9cjWjqkrZvOcgSzf2dK8dmMzsHOBq4JPeh08z7n6Tu89z93k1NTXH3L8orIGyIunQnzbb1/Yq+ScxA46S+sxTUi+S40IhY9HsvB0w20DPvXu99QZ2dSPwS2CLmQ0xsyFABAjHfy9JR4CHe+qV1IuQB21W8svUERWYwfrdTeo8yTAl9SJ5IDELzp9e3cahtryq/V5BrEa3uxnAymMcOx24nlgikXicBZwe//mT6QhQPfUiR8j5Niv5paw4zLih5bRHnQ31TUGHM6ApqRfJA8eNHMzM0ZXsO9TOw2/k1bTP9wCnx+esBsDMJhK70d9zjGPf0cPjFWKD+N4B3JGOADsHyqqnXgTyoM1K/knMgPPmDpXgZJKSepE8cfGcvCzBuRnYACw2s0VmthBYDGwm9lU9AGY2wczazeyriW3u/mj3B7AXaIz/viUdAR6e0jKajwORRdIt59us5B+tLJsdSuqToNlvJBcsnD2akMGjb+5iT1Nr0OEkxd2bgPOIzYZxK7HFaNYD57l71/nNDAgTwDUpHDJCBu7QEVVSL4UtH9qs5J/Dg2U1rWUmRYIOIB+4+xJgybx5864JOhYpXCMGl3LOtBoeW1XHva9u4+ozJgYdUlLcfRNw6TH22UAsSTjWueanJ6ojFUdCHGqL0toRJRJWjiKFLR/arOSXafHym1W71FOfSbp7ieSRxIDZPCvByXmJwbJt7eqpFxFJtyk1FYQMNtY359tkD3lFSb1IHnn3jFoGFYd5efNe1tXpa8x0KYnX1bd06GYjIpJupUVhJlQPoiPqrKvTDDiZoqReJI+UFYe5YNYoAO5+Sb316dLZU69VZUVEMmLaiFgJzmqV4GSMknqRPJMowbnr5a2arSVNOmfA0Vz1IiIZcXytVpbNNCX1Innm9MnV1FaWsnnPQZZuPNYCj5KMwz31SupFRDJhmmbAyTgl9SJ5JhwyFs0ZDWjAbLoUa1VZEZGMSixApbnqM0dJvUgeumTOWAD+9Oo2zSSQBkURrSorIpJJk4YPIhwyNu5p5mCr7luZoKQ+CVp8SnLN8bWDmTGqkn2H2nn0zV1Bh5P3StRTLyKSUSWRMBOry3GHNbtUgpMJSuqT4O5L3P3aqqqqoEMR6XTxHM1Zny4lRfEpLZXUi4hkzPRRlQCs3K5O0kxQUi+SpxbNHk3I4JE3d9HQ1Bp0OHmttCgMoK+ERUQy6MQxsc7R17Yqqc8EJfUieWpEZSlnTR1OW4dz72vbgw4nryWS+pZ2JfUiIplyOKnfF3AkA5OSepE81jln/bItAUeS38ri5TfqqRcRyZyZ8aT+9e37NIVwBiipF8lj75lZS3lxmGWb9rJht5be7q9ET71mEhIRyZyqsiImVJfT2h5ltearTzsl9SJ5rLw4wgUzawG46yUNmO2vskRNfZt6jkREMmlWvLd+uerq005JvUieuzhegnP3y1tx94CjyU+dA2XVUy8iklEaLJs5BZvUm9ltZvaqmb1kZs+b2TuDjkmkP86cMpyRlSVsrG9m2aaGoMPJS50DZZXUi4hklJL6zCnYpB64zt1Pcvc5wHXA7WZWyP89JE+FQ8ai2fEBsyrB6ZfOgbJK6kVEMmrW6MODZds1WDatcjaJNbOxZvZjM3vGzJrNzM1sYi/7jjOzO8ys0cz2mdmdZjb+aOd3971dftWqUpLXEgtR3fvqdq2K2g9lxZqnXkQkG6rKixg/rJyW9iirtbJsWuVsUg9MBa4AGoAnetvJzMqBh4ETgA8DVwHTgEfMbNDR3sDMfmBm64A/Ape6u7IhyUvTR1VyQu1g9ja38cibu4IOJ+90zn6jD0QiIhmnEpzMyOWk/nF3H+nu7wVuP8p+1wCTgYvc/W53XwwsBCYQK6vplbv/g7tPBj4E/KeZFacpdpGsOzxnvUpw+koryoqIZI9mwMmMnE3q+9BrvhB41t3XdDl2PfAUsAjAzK42s5fjj0/38F5/AYYCJ6YeuUgwFp48BjN4+I1d7G1uDTqcvFKmeepFRLJGPfWZkbNJfR/MBJb3sH0FMAPA3W9x99nxx0/NrMzMJiV2NLMzgGpgXVYiFsmA2qpSzpoynNaOKH96bXvQ4eQVLT4lIpI9s8ZUAhosm24DIakfRqzuvrs9xHrfe1IG/M7MlpvZy8D3iNXUv+U8ZnatmS01s6V1dXVpC1okExIDZlWC0zdlmqdeRCRrhpQXM25YGYfaoqyp02DZdBkISX2fufsedz/D3WfFe+/PcveHe9n3Jnef5+7zampqsh2qSJ9cMKuWsqIwSzc2sKm+Oehw8kapprQUEcmqzhKcLSrBSZeBkNQ30HOPfG89+H1mZgvM7KbGRv3Dk9w2qCTCe2aOBGIrzEpyDi8+pa+BRUSyQYNl028gJPUriNXVdzcDWJmON3D3Je5+bVWVprOX3Hfx3LFAbCEqdw84mvzQOU+9eupFRLJCg2XTbyAk9fcAp5vZ5MSG+CJVZ8VfEykoZ02pZnhFCet3N/Hy5r3HPiDD+rM4XPy4eWZ2k5m9EV+AbpOZ/bbrIPd0KY8n9U0t7ek+tUjeyYc2K/nvpDFDMIPl2/ZpkoI0yemk3swuM7PLgFPimy6Mbzu3y243AxuAxWa2yMwWAouBzcCNaYpD5TeSNyLhEItmjwZivfVBSmVxOOADxL6F+xFwIfBFYC6w1MzGpTPOsqIw4ZDR0h6lpV03Fylc+dJmJf9VlRcxvbaS1vYoyzalpVq64OV0Uk9s0anbgevjv98Q//0biR3cvQk4D1gF3Ar8FlgPnOfuaRlSrfIbyTeJWXCWvLKN1mBXSe334nDAd+KD2G9w98fc/XfABcTG0FyTziDNjMrSCAD7D6m3XgpaXrRZGRjOmFINwDNr6wOOZGDI6aTe3a2Xx/xu+21y90vdvdLdB7v7Re6+IV1xqKde8s3M0ZUcN7KChuY2HlsV6FSsx1wcrjfu/pbA3X0jUAeMSXOcVJYVAbDvYFu6Ty2ST/KmzUr+OzOe1D+tpD4tcjqpzxXqqZd8Y2ZcPCcxYHZLkKEcc3G4vjCz6cAI4PUU43qLweqpF4E8arOS/06dNIxwyHhl816NaUoDJfUiA9Si2aMxgwdf30VjcL3P/VkcrkdmFgF+TqzX75eph3akytJ4T/0h9dRLQcubNiv5b3BpEbPGVNEedV7YsCfocPKeknqRAWr0kDJOn1RNa3uU+17bHnQ46fAT4Ezgb3ta/Tmhv6tAdyb1B9VbJJImx2yzWrVdzlRdfdooqU+CauolX108N1bGemdws+CkZXE4M/sP4FrgY+7+wNH27e8q0JVlifIb9dRLQctqm9Wq7dKZ1K9TUp8qJfVJUE295KsLZ9VSEgnx/Po9bN7THEQIKS8OZ2ZfBv4F+Ht3vzWNsR1hsMpvRCCP2qwMDPMmDKMobCzf2khjs66/qVBSLzKADS4t4t0zawG455VtQYSQ0uJwZvb3wLeBL7v7TzIUI6DyG5G4vGmzMjCUFYeZM24oUYfn1qu3PhVK6kUGuEvic9bfuWwL7p7tt09qcTgzm2Bm7Wb21S7bPgD8EPgL8LCZnd7l0edZOI5F5TciQB61WRk4TlcJTlooqU+Cauoln509bTjVg4pZW9fEa1uz+2+4D4vDGRDmyGvSBfHtFwDPdHvckO5YD5ffqKdeClc+tVkZODRYNj0iQQeQD9x9CbBk3rx5WhFP8k5ROMSCk0fzq6c3cOeyrZw0dkhW39/dNwGXHmOfDcSSga7bPgJ8JFNxdZdYUVaLT0mhy5c2KwPHnPFDKImEeGPHfuoPtFBdURJ0SHlJPfUiBeCS+Cw4S17ZRltHNOBoclPnirIqvxERyaqSSJh5E2OTLj27TvPV95eSepECcOKYKqbUDKK+qZUnVmsu6J4kBsru1ewLIiJZd+aU4QA8/MaugCPJX0rqk6Caesl3ZsYlc8cCcNdLgcyCk/NGVsa+7t2571DAkYiIFJ4LZ8Vmart/xQ4OtXUEHE1+UlKfBM1TLwPBotmjAXhgxQ7N8NKDYYOKKQ6H2HeoneZWDZYVEcmmyTUVnDS2igMt7eqt7ycl9SIFYuzQck6dNIyW9ih/Wb4j6HByjpkxsirWW7+jUb31IiLZtvDkWOfT4pcDWwU9rympFykgiTnr73pJF8ye1FaWArBDJTgiIlm34OTRmMEjb9TRqJnI+kxJvUgBufDEURRHQjyzrp7tjQeDDifnjIwn9aqrFxHJvpGVpZw5pZrWjih/Wb496HDyjpJ6kQJSVVbE+dNH4g53a8DsW4yqiiX121V+IyISiEUnx75RXvyy7lF9paRepMBc1FmCswV3Dzia3NLZU6+kXkQkEBecWNv5jbLGN/WNkvokaEpLGUjOPa6GoeVFrNp5gJXb9wUdTk6prVJNvYhIkCpLizjv+BG4w72vqre+L5TUJ0FTWspAUhwJsSA+w8BdyzRgtqtE+Y16h0REgnPRnNg96m7NgtMnSupFCtDF8RKcxa9so70jGnA0uSNRfrNNSb2ISGDmHz+CwSURlm/dx+v6RjlpSupFCtDscUOYNHwQdftbeHptfdDh5IxRVWWUFYWp29/CnqbWoMMRESlIpUVhLj0ltgr6fz24OuBo8oeSepECZGZcNFtz1ncXDhkzRlcCsGKbxtCIiATlU/OnUBIJ8ZcVO1i+VdfjZCipFylQiRKcvyzfQVNLe8DR5I5Z8aR++VZ95SsiEpQRlaVcdfoEAH7w11UBR5MflNSLFKjx1eWcMmEoB9s6uH/FjqDDyRkzx8QGxC9XT72ISKCunz+F8uIwD72xi5c2NQQdTs5TUi9SwC6eoxKc7maNjiX1K/R1r4hIoIZXlPDhMycC8H311h+TkvokaJ56Gajef9IoisMhnlqzm52amx2AaSMrKI6E2FDfzL5DbUGHIyJS0K57+2QGl0R4YvVunl+/J+hwcpqS+iRonnoZqIaUF/OOE2qIOtyjJbkBKAqHmF47GIDlW/RBXkQkSEPKi/nY2ZMA+Pf7Xtc0zEehpF6kwCVKcO5UCU6n0yZXA/Dn5dsDjkRERD5+ziRGDC7hpU17+e4DbwYdTs5SUi9S4N5xwgiqyop4ffs+3tihGV/g8AedJa9sp6W9I+BoREQKW2VpET/+4BzCIePGx9bxgCZ36JGSepECVxIJ876TRgFw1zL11gNMH1XJjFGVNB5s46HXdwUdjohIwTttcjX//J7jAfjH219hY31TwBHlHiX1IsIl8Z7pxS9voyPqAUeTGxKrGf7xxS0BRyIiIgDXvn0y754xkv2H2vnkb5ZxqE3fpHalpF5EOGXCUMYNK2PHvkM8u64+6HBywqLZo4mEjEdX1Wk1QxGRHGBmfPfyk5lQXc7K7fv42K9eoLFZs5QlKKkXEcyMi2fHB8yqBAeIzY985Wnj6Yg6n/39SzS3atVdEZGgVZUVceNVpzC8ooSn19Zzyc+eUilOnJJ6EQHg4rmxcpO/LN/OwVZ9pQnwpQunM21EBWvrmvjWva8HHY6IiAAn1FZy96fP5PiRg1lb18RFP32KFzZoDvuCTurN7KNm5mZ2UdCxiARt0vBBzB43hKbWDh5YqZkFAMqKw/zog3MojoT43+c38fV7VqiGU0QkB4wdWs4dnzyD+cfX0NDcxgduepav3P0adftbgg4tMAWb1JvZROAa4NlgIxHJHZfMjZXg3JXGOevNbJyZ3WFmjWa2z8zuNLPxSR5bambfNbPtZnbQzJ4xs7enLbgkTB9Vyb9eNItIyPjV0xtY8OMnVWMvA1q+t1kpHINLi/jF1fO45pxJuDu/eXYT5373EX7w11XsL8AVwXMyqTezsWb24/jFoDnemz6xl337fPExsxDwC+DvgML9SCfSzftPig0OfWL17rT0dphZOfAwcALwYeAqYBrwiJkNSuIUvyT24furwPuB7cD9ZjY75eD64PJ547jrU2cxuWYQq3cd4P0/fpLLfvY0tz67kT1NrdkMRSSjBkqblcIRCYf48vtmcP/n3s67po+gubWD/3poNW/71wf59O+Wcf+KHQWz3kgk6AB6MRW4AngReAJ4d087dbn4tBC7+DjwbWIXn5PcvbeRE58HnnL3F80s3bGL5K1hg4qZf3wND76+iyWvbOtcmjsF1wCTgePdfQ2Amb0KrAauA77f24FmdjJwJfAxd/+f+LbHgBXAN4GFqQbXFyeOreJPf3cO3/nLG/z+hU0s3djA0o0NfOOeFUwfVckJtYM5vnYwJ9RWMnVEBdUVxRSFc7LfRORoBkyblcIybeRgfvHht/Hcunq+/9dVPLd+D396dTt/enU7g0sizJ0wlDnjhzBn/FBOHFPF0PIiBloOmKtJ/ePuPhLAzD5BL0k9/bj4mNks4FJAXweK9ODiOWN58PVd3PXS1nQk9QuBZxPtE8Dd15vZU8AijpIgxI9tA27rcmy7mf0e+KKZlbh7Vr9pKysO8/WFM/mn9xzPX1fu4O6XtvHkmt28trWR13ooyRlaXsTwihKqK4oZXlFCVVkR5cVhyoojlBeHGdTl57KiMJGwURQOUdT5fPjnSDhEUcgIhYywxZ5DBuGQEbLYIxzfNtBuVJJVA6rNSuE5bXI1t113Blsamrn31e3c8/I2Vm7fx2Or6nhsVV3nfoNLI0ysHsT46nJGVZZSHb9WVw8qZnBpERUlEQaXRigrDlNaFKY0ErsO57KcTOrdPZrkrse8+JjZ1cR65gFuBqLARGB1/MZXC9xkZmPd/Sdp+hNE8tY7p49gcGmE17Y2smbXfqaOGJzK6WYCi3vYvgK4PIlj17t7cw/HFhP7Rm9FKsH1V0VJhIvnjOXiOWNpPNjGG9v38ebO/byxYz9v7tjPxvom6ptaaWhuo6G5jdVZXpTWjFjibwZGLNHHMIOQGQbx7Udus/gHgsTPsd1i+1jnue2I9znafp17dvmM0fXjRvcPH0e+1nV7t/2S+MzS2web3g7t7Zy9bu/1TEeLqY/79/H8p02u5v+8d3ofj3qLAdlmpfCMHVrO9edO4fpzp7B170Fe2tTAS5v2smxTA6t27Gf/ofZeO2R6EwklOlnizyEjEjLC4cOdLYlrb/frqVm366QZf7z+jLR+UMjJpL4PjnnxcfdbgFu6vf6zxA9m9ijwQ3e/u6c3MLNrgWsBxo9PapyQSF4rLQrzvhNH8fsXNnP/ip2pJvXDgIYetu8BhqZwbOL1t8h2m60qK+K0ydWcNrn6iO0dUWdPUyv1TS3s3t/K7gMt7DvURnNrB82tHRxsbY8/x39v66A9GqWt3WmLRmnriNLe4bTGn9s6orR1OFF3OqKx52jU6XAn6nT+7A7u0O5OrCJRCsWIytJ0nCarbVb3WMmGMUPKGDOkjPefNBoAd6e+qZWN9U1srG+mbn8L9U2x63RDUysHWtrZf6idAy2x6/ShttijPeq0Rzti30floHxP6lO5+CTF3W8CbgKYN2+e7pBSED5xzmQunzeWuePT0oyyKlfabDhk1AwuoWZwSez7wCzxeJKfSPxj2yDqjnd5HQcn9nNim8c2xl4n9rv74Y8GHv/QcPi94sd0/txlvy77dDniiGN7fuXI17zbB5Pu79/jf4NePsz0tn9v+nr+/pyr1/37/A5QWZp/t/Rcaa9SWMyM4RUlDK8o4ZQJPfYP9SjWuRLrYEn83BF1olFoj0Zj19n4NbQj6oevoX7k9dTjF4RwKL2lkvl3BUgzd59/rH3MbAGwYOrUh3ysRAAAC4JJREFUqZkPSCQHTB1Rka5TNdDzB+zePpB3P3ZCL8fC4d4/6cLMCFv6bxZSMNRmRXqRGOuUq3I3suSkcvFJmrsvcfdrq6qq0nVKkUKxgliZXHczgJVJHDspPstV92NbgTVvPUREUqQ2K5Kn8j2pT+XiIyKZdw9wuplNTmyIrzlxVvy1o1kCFNFlcJ6ZRYC/AR7QLBoiGaE2K5Kn8j2pT+XikzQzW2BmNzU2ahVJkT66GdgALDazRWa2kNjg9s3AjYmdzGyCmbWb2VcT29z9JWJT4/3QzD5hZu8Efg9MAr6Wxb9BpJCozYrkqZxN6s3sMjO7DDglvunC+LZzu+yW1MUnVSq/Eemf+AJw5wGrgFuB3wLrgfPc/UCXXQ0I89Zr0keB/yG2qNyfgHHABe6+LMOhixQktVmR/JXLA2Vv7/b7DfHnx4D5ELv4mNl5wA+IXXwMeAj4XLeLj4gExN03EVvw7Wj7bKCHabnd/SCxdSY+3/01EckMtVmR/JSzSb27JzV1QzIXn1Rp9hsRERERyWU5W36TS1R+IyIiIiK5zLyvK2IUMDOrAzYGHUcfDQd2Bx1EDtN/n94l899mgrvXZCOY/kiyzerfQObpv3F2/P/27j1GlrJO4/j3OQEUWW94SxQQj5pdIXgJGmXdKAqKEkWNuxiDdyV4VzaYLGFViP4ju14iSrztioJGgqtyWBMvIKJRiDGgCEaJiopIonjgqKCA+vOPqoG2T09Ngz1TXTXfT9I509Vv93nr6el6f/NWddWgP7OOsaNjNt0W/nm1qB+5JN+uqkf33Y9lZT6r2yzZbJb17JMZbwxz3nhmvjqz6bYe+Xj4jSRJkjRwFvWSJEnSwFnUj9+H+u7AkjOf1W2WbDbLevbJjDeGOW88M1+d2XRbeD4eUy9JkiQNnDP1kiRJ0sBZ1EuSJEkDZ1EvAJK8OckVSf6S5Nl992cZmEm3JGcmuTTJJUm+leSQvvs0ryR7J/l0kh1JfpvkM0n26btfyy7JwUlqxu36qXb3TPKRJNcmuSHJuUkOmPF6d07yX0muSfKHJBcmecLGrVH/kuyV5JR23W9s89x3Rru5skqyJcnxSX6a5I9Jvptk5lXXkxyd5AdJbkrywySvXPwayrFkNnPpdkfGWIt6rfgy8DTga313ZImYSbdjqurhVfUo4BjgrCRLv01JchfgK8A/AS8GXgg8FDg/yR599m1AXg8cNHE7dOWBJAHOofnsvA54LrArTb57Tb3O/wBHA28BngFcA3wxySPXewWWyEOAI4HrgK93tJs3q7cBJwLvA54OXETz2Tx8slGSo4EPAv9H816dBZya5FV/5/poZ44ls5lLt9s/xlaVtwHcgL2AU4ALgRuBAvZdpe3ewKeBHcBvgc8A+8z5/3wVeHbf67tMGQ05k436HQIOBrYDW/pezzlyeAPwZ+AhE8seBPwJ+Pe++7fMt/Z9LuDQjjbPats8aWLZ3dvfj/dOLHtE2+6lE8t2AX4IbOt7XTcw0y0TP79i1udy3qyA+wI3ASdNPf884NKp5/4K+NhUu/+lucLlrn3nssHvgeNrjzkNPZeN+D2ad4xd+lk13Wqu2ZxNPgtpRt0Wnk+Sdyf5Cc1s33Or6i/r1PdFOgK4qKp+tLKgqq4EvkFTkOrvcwTwy6o6f2VBVe2gmb1/1lS7W4AzJ9r9CfgUcFiSO21Md/s152dm3qwOA3YDzph6/hnAAUke1N4/CLjPjHanA/cC/uX2rMMIOHbMx5y69T7GWtQPx9eq6n5VdTjNbtLVHA1spfmr93NVdTbNgPBAmt03Y2ZG3RaeT1UdW1VbgaOAk5Pstk59X6T9gctmLL8c2G+D+zJUn0jy5yS/SfLJqe8jdOW7T5J/mGh3ZVXdOKPdbjQDpBrzZrU/zUz9j2a0g9t+v/dv/51+n6bbbRaOHfMxp269j7EW9QNxO2ZA15yFTPKiJN9pb69ZfG/7sciMxmg986mqLwD3BHb6MuQS2pNmJmXadpp10Op2AO+kOUzkyTTHbx8KXJjkvm2brnzhtozXarfnIjo8EvNmtSdwfbX769dox4zX3JTZO77OxzG22zKMsRb147PmLGRVfbyqHtne3r+hvVsOztR2WzOfJLtP7MonyUE0u+1/siE9VC+q6pKqOq6qzqmqC6rqPTRfdLsfzZdnpTFzfJ2PY2y3dRtjLerH5w7NQiY5MckvaI61/EiSX8w4U8VYzJXRJstk0jz57A58MsllSb4D/DfN8X6znrdsrmP2Z2G19VaHqroYuAJ4TLuoK9+Vx+dpt33GY5vVvFldB9yjPQPRWu2Y8Zpm383xdT6Osd3WbYzdZaHd1GBV1Yk0p0FTy0xWV1XbaTbEQ3Q5tx1TPGk/4Psb3JcxWTnk43LgqTMe3w/4eVX9fqLdc5LcZepY8f2Am9n5uPDNbN6sLgfuBDyYv81vZXb0+xPtoPkcXNPRTgvgWDKbuazujo6xztSPj7OQazOjbmPPZxvwuCRbVxa0F/t5fPuYbockjwb+EfhWu2gb8IAkT5xoczfgmfxtvufQnL/+3yba7QI8D/hSVd20zl0fknmz+gLNWXKOmnr+C4DL2uN2oTnl3rWrtNtOc2yvdjb2beOimFO3dcvHmfrxcRZybWbUbez5fBh4LXB2kv+kmWF+G3AVzcV4tIoknwCuBC4GrgceBRwPXA28t222jaZoPCPJm2gGqeOBACevvFZVXZLkTOA9SXZtX/dVNNcMmC42Ry3Jv7Y/Htj++/QkvwZ+3X53Ya6squpXSd4FHJ/kdzTv0/NovtR8xES7W5K8meZiU1cD57ZtXga8rqpuXs/1HbCxbxsXxZy6rVs+ztSPj7OQazOjbqPOp6puoClgrqA5L/dKofrkiUNDNNtlNMXhR4EvAm+kuWjKY6vqWrj1DBDPoLla5KnAZ2ku9vWkqrpq6vVe2r7W24HP01yQ5WntcfqbyVnt7ZXt/VPb+ydNtJk3qxPaNm+geY8eDxxZVf8/2aiqPkDzh8GRbbvnA6/dxF/unMeot40LZE7d1i2f7HzmKy2ridmcQ2g2/q8Gbp3NadvsAXwX+AMwOQt5V+DhYy9azKib+UjSztw2zsecuvWdj0X9gCRZ7c26oKoOnmi3D/Bu4Ck0u7zPA95YVT9d7z72zYy6mY8k7cxt43zMqVvf+VjUS5IkSQPnMfWSJEnSwFnUS5IkSQNnUS9JkiQNnEW9JEmSNHAW9ZIkSdLAWdRLkiRJA2dRL0mSJA2cRb0kSZI0cBb1kiRJ0sBZ1EuSJEkDt0vfHZAkjU+SA4EXAAXsC7wCOAa4B/AA4K1V9ePeOihJI2NRL0laqCQPBV4CvL6qKslpwEXtsgBfBy4B3tlTFyVpdCzqJUmLdizwpqqq9v4ewHVV9c0kewPvAk7rq3OSNEa5bZsrLTd350vDkOSBVfWziftXA6dV1Qk9dkvSKhxfx8GZeg2Cu/Ol4Zgq6B8G3B84v78eSVqN4+t4ePYbDcWxwH/M2p0P/Bx350vL6hDgZuCbKwuSbO2vO5KmOL6OhEW9huIdVXXDxP1/Bs4FqKqrquq4qvpNP12TtCLJ7klOTnJAu+gpwKVVdWP7+BbguN46KGma4+tIePiNBsHd+dJgHE5TtF+c5BZgK7Bj4vETgNP76JiknTm+jocz9Roid+dLy+sCml31BwIvBx4H/DjJB5OcAlxUVRf22D9Jq3N8HTDPfqOll2R34CTg9Kr6XpKzgftX1WPax7cA76uqV/fZT0mShsTxdVw8/EZD4O58SZIWz/F1RJyp19JLcm/gZGDlizonAqcCf6TZTbitqr7cT+8kSRomx9dxsaiXJEmSBs4vykqSJEkDZ1EvSZIkDZxFvSRJkjRwFvWSJEnSwFnUS5IkSQNnUS9JkiQNnEW9JEmSNHAW9ZIkSdLAWdRLkiRJA2dRL0mSJA3cXwGZP1FnEKykjAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 864x288 with 3 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# using the variable ax for single a Axes\n",
"fig, ax = plt.subplots(1,3, figsize=(12,4))\n",
"ax[0].plot(xi, fun1(xi), lw=2)\n",
"ax[0].set_xscale('log'); ax[0].set_yscale('log')\n",
"ax[0].set_title('Log-Log')\n",
"ax[0].set_xlabel('$x$'); \n",
"ax[0].set_ylabel('$f(x)$'); \n",
"\n",
"ax[1].plot(xi, fun1(xi), lw=2)\n",
"ax[1].set_title('Lin-Lin')\n",
"ax[1].set_xlabel('$x$');\n",
"\n",
"ax[2].plot(xi, fun1(xi), lw=2)\n",
"ax[2].set_xscale('log');\n",
"ax[2].set_title('Log-Lin');\n",
"ax[2].set_xlabel('$x$');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Methods\n",
"\n",
"The integrand is smooth and appears to follow powerlaws over subdomains in log-log space. The method outlined in the [StackOverFlow post](https://scicomp.stackexchange.com/questions/20901/integral-in-log-log-space) and implemented in this Matlab [FileExchange](https://www.mathworks.com/matlabcentral/fileexchange/54650-integration-over-log-space) works something like this.\n",
"\n",
"Over each subinterval assume $\\log y$ is linear in $\\log x$.\n",
"$$\\log y = m \\log x + n$$\n",
"\n",
"This is equivalent to assuming $y = e^n x^m$. The values of $m$ and $n$ are:\n",
"$$m = \\dfrac{\\log y_2 - \\log y_1}{\\log x_2 - \\log x_1}; \\quad \\quad n = \\log y_1 - m \\log x_1$$\n",
"The exact formula for the integral is ($m \\neq 1$):\n",
"$$\\int_{x_1}^{x_2} y(x)\\,dx = \\dfrac{e^n}{m+1} \\left(x_2^{m+1} - x_1^{m+1} \\right) = \\dfrac{x_1 y_1}{m+1} \\left[ \\left(\\dfrac{x_2}{x_1}\\right)^{m+1} - 1 \\right]$$\n",
"For $m = 1$,\n",
"$$\\int_{x_1}^{x_2} y(x)\\,dx = e^n \\log \\left(\\dfrac{x_2}{x_1}\\right) = \\dfrac{y_1}{x_1^m} \\log \\left(\\dfrac{x_2}{x_1}\\right)$$\n",
"\n",
"As shown shortly, so long as $\\log y$ varies slowly with $\\log x$ [nothing approaching steep vertical lines] this is a decent method.\n",
"\n",
"I am going to compare this `LogTrapz1` method with three other methods.\n",
"\n",
"`LinTrapz1` is the regular trapezoidal method with uniform step size.\n",
"\n",
"`LinTrapz2` modifies the above to constant **log-scale** step size. This implies uneven step size $\\Delta x$. Note that for this setting the lower limit $x_\\min > 0$. In principle, we can use $x_\\min = \\epsilon$ and compute the correction to the integral between $x \\in [0, \\epsilon]$ using a linear trapezoidal method.\n",
"\n",
"One problem with `LogTrapz1` for integrating functions with large slopes $d \\log y/d \\log x$ (in stress relaxation example shown shortly) is that at large $x$, the function decays rapidly (steep negative slope) in log-log space. For example, a function like: $G(t) = g_1 \\exp(-t/\\tau_1) + g_2 \\exp(-t/\\tau_2)$.\n",
"$$I = \\int_{0}^{T} G(t)\\, dt = g_1 \\tau_1 \\left(1 - e^{-T/\\tau_1}\\right) + g_2 \\tau_2 \\left(1 - e^{-T/\\tau_2}\\right)$$\n",
"\n",
"Using `LogTrapz1` leads to values of $m$ that are negative and moderately large. Since $n$ is like the intercept, the formula $n_i = \\log y_i - m \\log x_i$ can lead to very large postive values for $n_i$ when $x_i$ is large. This causes overflow issues.\n",
"\n",
"`LogTrapz2` tries to rescue this idea by asserting that $\\log y$ ($\\log G$) is linear in $x$ ($t$) and not $\\log t$. The subsequent derivation is as follows.\n",
"\n",
"$$y = a e^{bx} \\quad \\text{ or } \\quad \\log y = \\log a + bx$$\n",
"It leads to linear interpolation for each interval:\n",
"$$b = \\dfrac{\\log y_2 - \\log y_1}{x_2 - x_1}; \\quad \\quad a = \\log y_1 e^{-bx_1}$$\n",
"The exact formula for the integral is ($m \\neq 1$):\n",
"$$\\int_{x_1}^{x_2} y(x) dx = \\dfrac{y_1}{b} \\left( e^{b(x_2 - x_1)} - 1 \\right)$$\n",
"\n",
"For our problem, at large $t$, $b$ is nonzero, negative and typically small ($\\mathcal{O}(1/t)$). The exponential thus acts on a small -ve values, which don't particularly cause under/overflow problems."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def LinTrapz1(n, xmin, xmax, func):\n",
" \"\"\"Trapezoidal Linear x-scale\"\"\"\n",
" x = np.linspace(xmin, xmax, n)\n",
" y = func(x)\n",
"\n",
" return 0.5*((x[1:]-x[:-1])*(y[1:]+y[:-1])).sum()\n",
"\n",
"def LinTrapz2(n, xmin, xmax, func):\n",
" \"\"\"Trapezoidal Log x-scale; xmin > 0\"\"\"\n",
"\n",
" eps0 = 1e-6\n",
" xmin = max(eps0, xmin)\n",
" x = np.geomspace(xmin, xmax, n)\n",
" y = func(x)\n",
" return 0.5*((x[1:]-x[:-1])*(y[1:]+y[:-1])).sum()\n",
"\n",
"def LogTrapz1(n, xmin, xmax, func):\n",
" \"\"\"assumes log(y) is linear in log(x); xmin > 0\"\"\"\n",
"\n",
" eps0 = 1e-6\n",
" xmin = max(eps0, xmin)\n",
" \n",
" x = np.geomspace(xmin, xmax, n)\n",
" y = func(x)\n",
" \n",
" X = np.log(x)\n",
" Y = np.log(y)\n",
" \n",
" m = np.diff(Y) / np.diff(X)\n",
" n = Y[:-1] - m * X[:-1]\n",
" \n",
" Intg = 0.0\n",
" \n",
" for i in range(len(m)):\n",
" if np.abs(m[i] + 1.0) < eps0:\n",
" print(\"triggered\", i,m[i])\n",
" Intg += np.exp(n[i]) * np.log(x[i+1]/x[i])\n",
" else:\n",
" Intg += np.exp(n[i])/(m[i]+1) * ( x[i+1]**(m[i]+1) - x[i]**(m[i]+1) );\n",
"# Intg += x[i]*y[i]/(m[i]+1.0) * ( (x[i+1]/x[i])**(m[i]+1.0) );\n",
" \n",
" return Intg\n",
"\n",
"def LogTrapz2(n, xmin, xmax, func):\n",
" \"\"\"assumes log(y) is linear in x; xmin > 0\"\"\"\n",
"\n",
" # assumes log(y) is linear in x\n",
" \n",
" eps0 = 1e-6\n",
" xmin = max(eps0, xmin)\n",
" \n",
" x = np.geomspace(xmin, xmax, n)\n",
" y = func(x)\n",
" Y = np.log(y)\n",
" \n",
" b = np.diff(Y) / np.diff(x)\n",
" #a = y[:-1] * np.exp(-b * x[:-1])\n",
" \n",
" Intg = 0.0\n",
" \n",
" for i in range(len(b)):\n",
" Intg += y[i]/b[i] * ( np.exp(b[i]*(x[i+1] - x[i])) - 1.0)\n",
" \n",
" return Intg\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Experiments\n",
"\n",
"Let's start with the first problem. The exact integral is given by,\n",
"\n",
"$$\\int_{x_\\min}^{x_\\max} \\dfrac{1}{1+ax^2} \\, dx = \\dfrac{\\tan^{-1}(\\sqrt{a} x_2) - \\tan^{-1}(\\sqrt{a} x_1)}{\\sqrt{a}}$$"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def Iexact2(a, t1, t2):\n",
" sqrt_a = np.sqrt(a)\n",
" Iexact = np.arctan(sqrt_a*t2) - np.arctan(sqrt_a*t1)\n",
" return Iexact/sqrt_a"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"15.597966604415644"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a = 0.01\n",
"tmin = 1e-2\n",
"tmax = 1e3\n",
"\n",
"Iexact2(a, tmin, tmax)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/sachin/.local/lib/python3.6/site-packages/scipy/integrate/quadrature.py:251: AccuracyWarning: maxiter (50) exceeded. Latest difference = 1.160097e-05\n",
" AccuracyWarning)\n"
]
},
{
"data": {
"text/plain": [
"(15.597929302702289, 1.1600965324376489e-05)"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from scipy.integrate import quadrature\n",
"quadrature(fun1, tmin, tmax)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let us vary the number of gridpoints (between $2^4$ and $2^7$) and compute the absolute error by comparing with the exact value."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"n = 2**np.arange(4,8)\n",
"err = np.zeros((len(n),4))\n",
"tmin = 1e-2\n",
"tmax = 1e3\n",
"Iex = Iexact2(0.01, tmin, tmax) \n",
"\n",
"for i, n1 in enumerate(n):\n",
" err[i,0] = np.abs(LinTrapz1(n1, tmin, tmax, fun1) - Iex)\n",
" err[i,1] = np.abs(LinTrapz2(n1, tmin, tmax, fun1) - Iex)\n",
" err[i,2] = np.abs(LogTrapz1(n1, tmin, tmax, fun1) - Iex)\n",
" err[i,3] = np.abs(LogTrapz2(n1, tmin, tmax, fun1) - Iex)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x7faafb4137b8>"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAEqCAYAAAAh0uZRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdeXzU1bn48c/Jvq8kEEhCCJCwLxIWWYMIAoparW21Ktyq2FurtfXeW9pqxVp/be+v1dtfb7WCWm3dilWxLqCiJIACsoiyhkAIWSD7vicz5/fHTCaZZAaSySSTyTzv18uXme/3zMzJDJlnzvk+5zxKa40QQgjR37xc3QEhhBCeQQKOEEKIASEBRwghxICQgCOEEGJASMARQggxICTgCCGEGBAScIQQQgwICThCCCEGhAQcIexQSiW7ug89MZj6qZQar5S6SSm1USk1y9X9EYOLBBwhbDB/iM/rQbvRSqnvDECX7D1/j/o5gNYAF4Angf8A179GYvCQgCOEbd8HXrtcI631eSBYKTXZGU+qlJqhlPrDJc6HKaUmdDrUrZ+Xe4z+pLV+Umu9D4gHcs3HnPoaCfclAUcIM6XUj5VS/0cp9SOgQPd8o8FXgB864fkfAh4Foi/R7FtAnbn9dLr0s4eP0d72lj512P7jKuAbwBOdDjvlNRLuTQKOEIBSKhz4NvAvIBTY2dP7aq2bAH+lVFhf+qC1/gPwzmWaJWitC8w/r6FLP3v4GO36a8SxBvgTMLL9gLNeI+HefFzdASEGibnAEa31PqXUz4Df9PL+R4D5wPb2A0qpscA9l7jPPq311p4+gXkq7VSnQ7Md6Ge/UkrdBPwceADIBB7vdLrbayQ8iwQc4fGUUvOAB4EL5g/MIK21oUubGwADsAg4CqwEntBanzQ3uQCMp9OHqdb6LLDBiV29EfifTre79dOZlFLjgXXAHuAWYLvWeov5nL3X4y3gLTsP2e01Ep5FAo7weOZRTSPwR631UaWU1bUGpdRo4KTW+rRS6lfA74BqIK9Tsyogtb/6qJTyBnzNU1PtvHv5GJOAOzsdWqiUCuh0e7fW+n1z22BMgWOp1rpMKXU/cNx8rievhy39+hqJwU8CjhAmE4ET5p/bOp8wZ1mhlBoO1GqtK4B3u9w/EKjvfMDJU2rLgY+7HGuz1dAerfUJOo24lFIbtdYb7TS/CThuDjY+wBit9XHz4/Tk9bCl22skPIsEHOHxzB+cZZ2mp4qVUqFa61rz+YmAPzAT2GU+tkZr3flDNgoo6vy4Tp5Sm2cjOFj108ligMPmn9OBL5RSK4BPgBQu/3rY0u01Ep5FstSEMCUMfN7pdiYwp9PtFcB1mP5eAszXeUq6PMY04LO+dMI8bXUXsFQp9Zg5cw6lVARQaeMuXftp9zEc8BowSim1GhgD1ALDzEG5J6+HLX1+jYR7Uz1faiDE0GLeeuUeTB/mW7TWX5qPRwH/obX+eS8e6wWt9ff6qZ/rgfe01he6HO91P7vc/1at9WUXtzpLf75Gwj3ICEd4MiNQCJS3BxsA8zWJCqVUTE8eRCk1h+7XV5xpZNdgA73vp437D2Sw6e/XSLgBCTjCY2mtv9RaP661/r2N008B37zcY5gvqC8D/uHs/pkffyzw9SWa9KifrtTfr5FwHzKlJkQfKKXigBqttWRf2SGvkWgnAUcIIcSAkLToLpRSa4A1oaGh96SkpLi6O0II4VYOHTpUprW2eV1RRjh2pKWl6YMHD7q6G0II4VaUUoe01mm2zknSgBBCiAEhAacLpdQapdSm6upqV3dFCCGGFAk4XWit39Varw8Pd3SBthBCCFsk4AghhBgQEnCEEEIMCAk4QgghBoQEnC4kaUAIIfqHBJwuJGlACCH6hwQcIYQQA0ICjhBCiAEhAccDGIyyfZEQwvVk804P8PK+87z4eS7pqTGkp8Yyd0wUAb7eru6WEMLDSMDxADuzSjhXVs+5snr++lkugb7ezB8bbQlACVFBru6iEMIDSMDpor08wbhx41zdFadoMxj5usA6xbux1cAnp0r45FQJcJxxsSGkp8SwdEIsaUmR+PvI6EcI4XxSnsCOoVSeoKnVwBfnKsjIKiUjq4ScMvuFF4P8vFkwbphl9DMqInAAeyqEcHeXKk8gAceOoRRwujpfXk9GVik7s0rYe7ac5jaj3bapw0MtwSctKRJfb8kzEULYJwHHAUM54HTW1Gpgb045mVmlfHqqhLyKBrttQ/19WDBuGEsnxLAkJZYR4QED2FMhhDuQgOMATwk4nWmtOVfWMfrZf66ClkuMfibGhbHUPPq5IjECHxn9COHxPD7gKKUeAe4AxgE3aa23Xu4+nhhwumpoaWPv2XJ2ZpWw81QphVWNdtuGBfiwaHwM6akxLEmNITZURj9CeCIJOErNA0qAF4D/kYDTe1przpbWWUY/X5yroNVg/9/OlFFhLE2NJT01hhkJkXh7qQHsrRDCVdwu4Cil4oGfAmnAdCAQGKO1zrXRNgF4ClgOKGAH8KDWOs9G2wwk4DhFXXMbn58pY2dWKZlZJVyobrLbNiLIl8Xto5+UGKJD/Aewp0KIgXSpgDNY1+GMA74FHAJ2AytsNVJKBQGfAs3AWkADvwZ2KqWmaa3t5/+KPgnx92HF5BGsmDwCrTWni+vYmVVCRlYJB3Mraeu0nU5VQyv/+uoC//rqAkrBtFHhpKfGsnRCLNNGheMlox8hPMJgDTi7tNbDAZRSd2Mn4AD3AMlAqtb6jLn910A2cC/w5AD01eMppUgdEUrqiFC+v2QsNU2tptHPKdP0W0lts6Wt1vBVQTVfFVTzx0+yiQr2Y0mKafSzeHwMkcF+LvxNhBD9aVAGHK21/dQoa9cD+9qDjfm+55RSnwE3IAHHJcICfFk5JY6VU+LQWnPyYi07s0rIzCrlUF6l1WaiFfUtvP1lIW9/WYiXghkJEabRT2osk0eGyehHiCFkUAacXpgMvGPj+HHglt4+mFJqPbAeIDExsW89E4Bp9DNpZBiTRoZx39JxVDe0svtMqXnXg1LK6jpGP0YNh/OqOJxXxZMfn2ZYiD9LUmJYOiGGReNiCA/ydeFvIoToq0GZNNCZeUptMzaSBpRSLcCTWusNXY7/GtigtfYx394I3A3EALVAEzBPa11g73klaaD/GY2a4xdqyMgqYWdWCUfyq7BXScHbS3FFomn0k54aw6S4MJSS0Y8Qg407Jg04ldZ6I7CxJ22H2uadg5mXl2JqfDhT48O5f9l4Kutb2JVdSmZWKRmnS6mob7G0NRg1B3IrOZBbyf/9MIvhYebRT2osC8YPIyxARj9CDHbuPsIpBrZqre/tcvxp4BatdYyjzysjHNcyGjVfF1abRz+lfF1Qhb1/qj5eilmjI1k6wTT6SR0eKqMfIVzE7dbhdHaZgPMp4Ke1XtjleAam322JA8/XPsK5Jzs72+F+C+cqr2tmV3YpO0+Vsiu7lKqGVrtt48IDLFNvC8YNI8TfIwbyQgwKQzngPAj8HkjRWueYjyVhSoveoLX+g6PPKyOcwctg1BzJryIjq4SMrFKOFlbbbevrrZidFMXS1FiWTohhbEyIjH6E6EduGXCUUt80/7gM+D7wA6AUKNVaZ5rbBANfAY3Aw5gWfj4OhALTtNZ1DjyvjHDcTEltE7tOl7Ezq4Rdp0upbWqz23ZURCBLJ8SQnhLL/HHRBPnJ6EcIZ3LXgGOvY5la6/RO7RKx3trmE0xb2+T25fllhOOe2gxGvsyvYucp0+jnxMUau239vL2YmxxlXvcTw5hhwTL6EaKP3DLguJoEnKGhqLqJzNOm4LM7u4y6Zvujn8SoIFO5hQmxXJkcTYCvlNoWorck4PSCTKkNXa0GIwdzK8k4XULGqVKyimvttvX38eLKsdFcNSGWG2aMIjxQ0q6F6AkJOA6QEc7Qd6Gq0bzjQQl7zpTR0GKw2S7Yz5tb5yTyvYVjGBkROMC9FMK9SMBxgAQcz9LcZjCNfszrfs6UdM838fFSrJk+knsWJTNpZJgLeinE4CcBpxdkSk0A5Fc0sDOrhJf3ned0cffgs2j8MO5dPJYF46Il0UCITiTgOEBGOAJMlU4zskp5dtdZ9uVUdDs/eWQY6xcnc+3UOHy8vVzQQyEGFwk4DpCAI7r6Kr+KTbty2HbsYrdNRkdFBHLXwjF8e3YCwbKzgfBgEnAcIAFH2HO+vJ7n95xjy8F8mlqtSzeFB/py+7xE1s5PIjY0wEU9FMJ1JOD0glzDET1VUd/C3/ee56W9uVY7WwP4+Xhx08xR3LM4mbExIa7poBAuIAHHATLCET3V2GLgn4cLeG53DufLG7qdXz5pOPcuTiYtKcoFvRNiYEnAcYAEHNFbBqPmo+NFPLsrhyP5Vd3OX5EYwfrFY1k+aTjeUjpbDFEScBwgAUc4SmtTsbhNu86y42RJt/NjhgVz96Ix3HxFvGyfI4YcCTi9INdwhDOdKall064ctn55gRaDdYLBsBA/1l6ZxB1XjiYiyM9FPRTCuSTgOEBGOMKZSmqa+Ovnuby873y38gmBvt58e3YCdy0cQ0JUkIt6KIRzSMBxgAQc0R/qmtt4/Ys8XthzjgvVTVbnvBSsnhrHvYvHMjU+3EU9FKJvJOA4QAKO6E+tBiPvf32Rv2Se5VRR912r54+NZv3iZJakxMjWOcKtSMBxgAQcMRC01uzOLmPTrhz2nCnrdn7CiFDuWZTMmukj8fORrXPE4CcBxwEScMRAO1ZYzebdObz39UUMXfbOGREWwF0Lx/CdOQmEBkhtHjF4ScBxgAQc4SoFlQ08v+cc/ziQ361GT6i/D7fNS+R7C8YwPEy2zhGDjwScXpC0aDFYVDe08vL+8/z1s1zK6pqtzvl6K26YMYr1i5NJGR7qoh4K0Z0EHAfICEcMFk2tBrZ+Wcim3TnklNZ3O780NYb1i8cyLzlKEgyEy0nAccCQCjjnP4czO2DCtTDyCpAPJbdkNGp2nCxm064cDp6v7HZ+enw46xePZeWUEbJ1jnAZCTgOGFIB51/3w+G/mX4OG2UKPBOug9HzwVsuQLujQ+dNW+d8dKKYrn/CiVFB3L1oDLfMSiDQT7bOEQNLAo4DhkzAMRrg9ynQ0D3lloAISF1lCj5jrwI/WeXubnJK63huzzn+eaiAljbrrXMig3y548ok1l45mugQfxf1UHgajw44SqmxwEtALFAP3KO1vmwkGTIBx9AGp7fDqfcgaxs0dd/FGACfQBi3zBR8Uq6BINlK352U1TXzt89z+du+81Q1tFqd8/fx4pa0eO5emEzSsGAX9VB4Ck8POB8DW7TWm5VSy4H/BSboy/ziQybgdGZohfOfwan3Tf/VFNpup7whaaEp+Ey4FsJHDWw/hcMaWtrYciCf5/aco6Cy0eqcUrBy8gjWL05mZmKki3oohjq3CjhKqXjgp0AaMB0IBMZorXNttE0AngKWAwrYATyotc4zn48BcoAorXWr+dhp4LbLjXKGZMDpTGu4cBhOvmcKPmVZ9tuOvMIUeCaugZjUgeujcFibwci2Y0Vs2pXD0cLqbufnJEVx75JklqbG4iUJBsKJ3C3gpAP/AA4B3sAKbAQcpVQQ8BXQDDwMaODXQBAwTWtdr5SaBbymtU7pdL+PgL9ord+6VD+GfMDpqiwbTr5rCj6Fl/i9o8fDxOtMo5+RV4CXbLcymGmt2ZtTzrOZOWSeLu12flxsCOsXJXPDzJH4+0iCgeg7dws4Xlpro/nnu4HN2A44PwKeBFK11mfMx8YA2cB/aa2flIDjoJoLHdNuubvB2Ga7XWhcR8Zb0kLJeBvkThXVsGlXDv86coG2LlvnxIb6s25BEt+dO5rwQHkfhePcKuB0dpmA8wkQoLVe0OV4JoDWeolMqTlBYyWc/ghOvQtnPoHWBtvtAsIhZaUp+IxbBn5ycXqwuljdyF8/y+XV/XnUNVt/mQj28+Y7cxL53sIxjIoIdFEPhTsbqgGnCHhHa31vl+NPA7dorWPMtz8BXu+UNPA0kGIraUAptR5YD5CYmDjr/Pnz/fBbubHWRji7syPjrbHCdjufAFOa9YTrTGnXkvE2KNU0tfLa/jxe+OwcxTXWW+f4eCnWTB/JPYuSmTQyzEU9FO5oqAacFuBJrfWGLsd/DWzQWvuYb4/HlBY9DGgA1mutv7jcc8sI5zIMbZC31xR8Tr0P1fm22ylv0wLT9oy3iISB7ae4rJY2I+8cKWTz7hxOF9d1O79o/DDuXTyWBeOiZesccVkeHXAceE7ZvLO3tIaLX5mCz8n3oPSk/bZx02HCGlPiQcwE2WZnENFak5FVyrO7zrIvp/vodfLIMNYvTubaqXH4eEuyiLBtqAacYmDr5abUHCUjnD4oP9sRfAouMZiMGtuR8TYqTTLeBpGvC6p4dlcO245epEt+AaMiArlr4Ri+PTuBYH+HvteJIWyoBpxPAT+t9cIuxzMw/V5LHHxOGeE4U20RZH1gCj7ndoGx1Xa7kBEwYbVp2i1pMfj4DWw/hU155Q08tyeHLQfzaWq13jonPNCX2+clsnZ+ErGhUptHmAzVgPMg8HtMCQA55mNJmNKiN2it/9CX55YRTj9orILsj02jn+yPobX7VvsA+IdDygpzxtvV4B8ysP0U3VTUt/D3ved5aW8uFfUtVuf8vL246YpR3LM4mbEx8l55OrcLOEqpb5p/XAZ8H/gBUAqUaq0zzW2CMS38bKRj4efjQCimhZ/dr3727LllhDMQWpsgJ8OUbp21DRrKbbfz9oexSzsy3oKHDWg3hbWmVgP/PFTA5t05nC/vniK/fNJw7l2cTFqSZCZ6KncMOPY6lam1Tu/ULhHrrW0+wbS1TW5f+yAjnAFkNEDevo7rPtV5ttspL0i8siPjLXL0wPZTWBiMmo+OF/HsrhyO5HffEHbR+GH88TsziQqWqVFP43YBZzCQgOMiWkPR0Y7gU3LcftsR00zBZ+J1EDtJMt5cQGvNgVxTbZ4dJ0usziVGBfHCujTGxUoJbE8iAacXZEptkKnIMa3zOfke5O/HNHNqQ+SYjg1G4+dIxpsLnCmpZdOuHN44VGApChca4MPT372CReP7lDQq3IgEHAfICGcQqivplPGWCYYW2+2CY80Zb9fBmMXgI8XHBtKHx4t48PUjNLYaAPD2Ujx2/WRunydToJ5AAo4DJOAMck01cOZjU/DJ/hhaam238wvtyHgbvxz8ZXpnIBwrrObulw5SVNNkObZufhIPXztRFo0OcRJwekGm1NxQW7Npjc/Jd00joPru2/AD4O0HyenmjLfVECLTPP2puKaJu186aFWPJz01hj/dOpPQANmReqiSgOMAGeG4KaMBCg6Ya/u8B5W5dhoqSJzXkXQQmTSAnfQcjS0GHnrjCB8cLbIcSxkewvNrZ5MQFeTCnon+IgHHAUMp4BiMBryUF23FxfgMH+45GzBqDcXHzbV93jVlv9kzfEpH8Bk+RTLenMho1Dz58Wn+d+cZy7HoYD823TmLWaNlvc5QIwHHAUMp4GzJ2sJrh1/g8SfyUXGxRK9YRfjyFQRMm4bypGyuylw49YFp5JO3F7TRdruI0R3BJ2EueEklTGd463ABG948SovB9Lr7eXvx39+cxo0zR7m4Z8KZJOD0wlC8hvPvO/4dw45d/Hir9QdsW3Q4EVcvJ3L5SoLnzEb5edAivfoy0/WeU++bavwYmm23CxpmLiy32nT9RwrL9cmB3Aru/fshq+1xHrhqHA9enYKXl4wqhwIJOA4YKiOcVmMry7YsY96ecm7NNBJgZ+9MQ5A/gUsWEbNyDSGLFuIV5EHz6821cGaHKfic/hCaa2y38wmAMUtMW+ykrISwuIHt5xCRV97AXS8dILukY/epa6fF8YdbphPgK6NJdycBxwFDJeAANBuaOVB0gF1nd1Cc+TEpRytJy9aENdpuXzJvHLF/+G8mRE3wnOs97dpaIHeXKd066wOoK7bfduQVpmy31FUwfLJc9+mFmqZWfvjql+w63ZFROD0hgs13zCI2THaedmcScBwwlAJOZ1prTlScIDP3U87t3sbwg7nMOa2J6fSl/k9rvNg9xYvhQcNJT0hnacJSJucaCR4zFt+RI13X+YFmNMKFw6bAk7UNSk7YbxueaAo8qatg9AIpr9ADbQYjj793gpf2dpRyjwsP4Lm1aUweGe7Cnom+kIDTC0PxGs6lXKy7SEb+To7vfY+gz75mZraBx27zpj6w49u6Mmo2/8lIWIOmaewoYlZeS8zK6/AbN86zRkAV5+D0dlMAOv85GNtst/MPM5VVSF0N46+GwMiB7aeb+dveXB579wQGc6W3ID9v/vidmSyfNNy1HRMOkYDjgKE6wrmU2pZaPrvwGRn5Gewq2EWtefX+hHzNr142dGvfFBdJ6NXLGXXtNzwv462xynTdJ2ubaaeD5mrb7ZQ3jJ5vnnpbCVHJA9tPN5F5upQfvnKY2mZTEFcKfr5qIncvGuNZX2qGAAk4DvDEgNNZq7GVIyVH+DTvU87t2caSj0uYcl7jYyeTuCkiCN8l8xl9/XcIXbBgYDvraoZW04gnaxtkvQ9VdsorAMRM7Jh6k7LaVrKLa/neSwfIr+i4uPjttAQev3EKfj7yOrkLCTgO8PSA05nWmrNVZ9mdtZ2LH7/PiEPnmXlW28x4O5vgw4HHbiI9Pp15I+cR6BM48B12Ja2h5GTHdZ/CS/wbCo6BlGtMo5/kdEm5Bsrrmrn374c4eL7ScmxechR/uX0WEUFyXcwdSMBxgAQc+8oay9h1dgc5O7YSsvc4M0+3WTLe/r7Ui3fnmb6N+nv7c2XclVyfE8nE0LGMuOZ6fCI97HpGbZEp1TprG+TshLYm2+18AkxBpz3lOnTEQPZyUGluM/Czt47y1uFCy7Exw4J5fm0ayVLCetCTgOMACTg909jWyL78zzie8SZk7mfrjBaKI63n3H+/uY3EMjAqqJ2UQOSKlYy57tv4jfKwFeYtDaay2lkfmJIP7G0yCjBqlnnqbbVHFpfTWvN0xln+74dZlmPhgb48890rmD9OyowPZhJwHDDUAo7Wmq+++orw8HCio6MJDQ11+sVYg9HA0bKjZORnkJGfwdnqs4yo0Py/Z7snHABUJUUTeFU642+4naCUVM+6OGw0QuGhjqm30pP220Ykdqz3Gb0AvD1np+UPjl7kJ1uO0NRqunjo46V4/MYp3Don0cU9E/ZIwOmFoZoWXVNTw+HDhy23fX19iY6OJjo6msjISHx8fJz+nHk1eew+9j7VW7cy8nA+4ws09i791sQGoxfPYfJPHyc0NNrpfRn0KnIgazuc3ga5n4G2HaTxDzelWqeuhnHLPCLl+uuCKu5+6SAltR3bD929cAw/Wz0Rb9kOZ9CRgOOAoTbCOXfuHOfPn7d5TilFRESEJQAFBjr/Qn9VUxWfH/2Awm1vE77/JBPPGbplvBVHwI9/4M+cuLmkJ6STHp9OXIgHbh/TWAlnPjGNfrI/tr/VjpcPJF7ZMfqJGjOw/RxAF6sbufulgxy/0PFaLJsQyx9vnUmIv/O/LAnHScBxwFALOG1tbVRWVlJeXk55eTmtrXY2VQNiY2OZNGlSv/WlxdDCwTOZZH/wOr57DjPxdBMBrfDuHMXfl1nvpfWds8OZXxhM7Ko1TFj5HbxDPOyicVsL5LWnXH/Qw5Tr1aZrQEMs5bqhpY0HXz/CRyc6thuaMCKU59fNZlSEh2VDDmIScBww1AJOZ1pramtrLcGnrq7O6nxycjKJidZz5HV1dfj7++Pr69zrB1prTl38mqPbXiGDU+z2OWd1/uHXDEzLNf0bbfGBsimjCFu+nCk3rCNwmIetRNfatL2OJeX6kP22wbFdUq6HxmasRqPmvz/M4i+ZZy3HhoX4s/nOWcxMHPrTi+5AAo4DhnLA6aq5udkSfCorK5k1axbBwdZrQg4cOEB9fb0l6SA6OpqgoCCnX+gvqi+yJB0cO7ePvzzVjLeNf6JGBUXjovBdupDJ3/g3osdMcGo/3EJtkXmrnW2m7LdLplwv7ZRy7f6B+o2D+fz87aO0Gkz/OPx8vPj9LdO5froH7fU3SHl8wFFKPQLcAYwDbtJab73cfTwp4HRmMBjw8vKyCiRNTU3s27evW9uAgABL8ImIiMDLyVM4tc21HNi9hZLt7xL9xRniS+xcSAeK4oMo+vmdzE+7kdFho53aD7fQUt8p5frDy6Rcp3VKuZ7otinX+3PKufflQ1Q1dEwPP3j1eH60bLxnZTwOMhJwlJoHlAAvAP8jAad3amtryc7OpqbGzsVrwMvLi8jISKKjo4mLi3P6H3ybsY0jh7dz/r0tBH7+NaPzmq0y3hr84a4feWPwViSHJ5t2uY5bzJTYafj4eE4aMdDLlOvR1rtcu1nKdW5ZPd976QA5pfWWY9dPH8l/f3Oa1NZxkUEfcJRS8cBPgTRgOhAIjNFa59pomwA8BSwHFLADeFBrfYmrqZb7ZiABx2EtLS1UVFRQXl5ORUUFBkP3EUdAQABz587t12+YWmvO5hzi5Dt/Q+3az+jTNeyboPh/N1h/wMzJMnL3x1A6K4moFauZufJ2goM8cNv79pTr9l2ue5RyfTUERgxsPx1U3dDKfa8eZs+ZMsuxmYkRbLojjZhQfxf2zDO5Q8BJB/4BHAK8gRXYCDhKqSDgK6AZeBjQwK+BIGCa1rqeS5CA4zxGo5Hq6mrLtZ/GRtPeNqNGjWL8+PFWbYuKiqiurras+fH2du43z5KSXPaf2cnHTV+y98Jemgymaxk//JeBxcc7/n3X+0Ph1OEELktnxpp/I3aYB069NVZC9g5T8Dmz49Ip1+27XKesHPQp160GIxv/dZxX9nd87xwVEcjz69KYMCLMhT3zPO4QcLy01kbzz3cDm7EdcH4EPAmkaq3PmI+NAbKB/9JaP3mZ58lAAk6/aGhooLy8nPDwcMLCrP/Av/76ayoqKgDTmp/IyEiioqL6Zc1PY1sj+y/uJyNvJ+n/9SbDy21/m2/xhtzUMNTieUy44U7GJ13hefP+bS1w/jNzyvU2qL7EJEHspI7rPiOvGJQp11prXvw8l8ffO4G5tA7Bft786baZXDXB/RMl3MWgDzidXSbgfAIEaK0XdDmeCaC1XrUi8LcAACAASURBVKKUuhP4ifnUZq31nzu1y0ACzoAyGAzs2bMHe//OgoKCLIkH4eHhTv3QN7S1cnznm1z44C3C9p0kstJ2wTSjgufuiGXE1atJT0jniuFX4OvlXtcy+kxrKD7esd7nwmH7bYNjTbV9UlfDmCWDLuV656kS7n/tS+rMtXW8FPzi2kl8b0GS532pcIGhFHCKgHe01vd2Of40cIvWOuYyj52BBJwBdbk1P535+Pgwc+bMbinZzurH+YMZnP3XK/ju+ZKYiw2WcwYF9zzgTV2Q6cMo1C+UxXELWc5EZl95M2H+HjglU3PROuXa0Gy7nU8gjDWnXI+/ZtCkXJ8qquGuFw9SWNVRW+e2uYk8dv1kfL0H3+hsKBlKAacFeFJrvaHL8V8DG7TWNve4UEptBO4GYoBaoAmYp7Uu6NJuPbAeIDExcZa9rWCE45qamiyJB5WVlRiNHfvbeHt7s2DBAqv0aoPBQFNTk9PX/JRnH+PE1hcxZHxOmbGGx75j/dipBZrH/26gKFJRMDOOsKuWccWy75IQ6YHXfVrq4exOU/A5vR0ayuy3taRcr3L5Ltdldc2s/9tBDudVWY4tGBfN07fNIjzIw0awA8jjA04vn39Ibt45GBkMBqqqqiyjn7CwMCZPnmzVprS0lOPHj/frmp+WpgYOVhwhIz+Dnfk7Kaov4vZPDVy/3/pvozYAsieF4b1wLhNWfYepo+fi7eVhqbdGQ5eU61P220YkmhIOUlZC0kLwGfiMsaZWAz9982veOXLBcix5WDAvrJtN0jApeNcfhlLAKQa2Ojql1hsypTawtNYYDIZuu1afOnWKoqIiq2Pe3t6WNT/R0dH4+TmvEqTWmqzKLPKe+BWxH3+Ff4vtmtptXpA1xo+GeZNIWP1N5k5bRZDv4LqWMSDKz3ZMvV0q5dovBMZeZZ56WwHBA1fTRmvNnz49w5Mfn7Yciwjy5S+3z2JesgfuTN7PnBpwlFLXaq3ft3MuSGvdYOtcLx7/UgHnU8BPa72wy/EMTL/Lkr48t/mxZIQziJw+fZri4mKba37ahYaGEh0dTWxsLEFBzvvQNzY3U5j5Iec/eAP/vV8TUt1is93fl3rx4fwA5sTNYWnCUhbHL2ZEsAdW7GyoMKVaZ31g2u3aXso1CuJnmxIPUlYN2G4H7319gYe2fEVzm+lLhK+34okbp/Kt2Qn9/tyexNkB502t9c3mn+/SWj/f6dw9wF6t9bE+dPZSAedB4PdAitY6x3wsCVNa9Aat9R8cfd6uZIQzeNhb89PVhAkTGDGifz7otdZUfXWY7PdeoS1zL5H5HdcFfrTem4vR1h+YPzgQReSsuUy95lYmxc3wvOwoyy7X5ho/lbn220YkmgJPyjX9PvV2JN9UW6esriMJ4t7Fyfx05QS8pLaOUzg74GzVWt9o/vktrfVNnc75A49orR92oJPfNP+4DPg+8AOgFCjVWmea2wRjWvjZSMfCz8eBUEwLP+2nQPWSBJzBSWtNY2OjJfhUV1dbUq7nz59vNb2mtebEiROWWj8BAQFO60dzYQGn33uV4sN7efpaRXZlx2h4eIXmT+Yqp02+cGp8AC1XTid59beYnXoVAT7O64db0BrKTnckHeTvB217qnIgpt4Kq0y1dU5e7BiBLZ80nP/59gyCpbZOnzk74PwReEtrnamUeltr/Y0u5x/XWj/iQCftdSRTa53eqV0i1lvbfIJpa5vc3j6nnX7IlJobaWtro6KigoaGBpKSkqzOda1yGhwcbLnuExYW5tRRR2FdIRn5GWTmZzLsX3u5Y0f3NT9GBdnx3lSkJRO74lrmzbuZYYEDdy1j0KgvhzMfm4LP5abeEuaYRj5Onnqrb27jR69/yY6TJZZjk+LCeH5dGnHhUlunL5wdcFKBbcC3gZ91GeFMBr6ntX6oD/0dFGSE4/4uVeXUx8fHqsS2M+v8VHx1iOwtz8PuLwgrsb/b0oUoOL5gFAF3fIsl8UtIiUyRqbeeTL2lroTRC8Gnb8kiBqPmt9tOsnl3Rw2m2FB/Nt+ZxvQE99hHbjByepaaUuqHwB+BCuB5oAzT1v/fBJZorY873l3XkhHO0NHU1GSZequqqrJa89OZUoqRI0d22wOur7TWNJ7NJvvd16jbuZOI7GK8uvy5vTdb8berTanVI4NHsiRhCenxS5g9Yg6+brZzc59pDaVZpsCTtR0KvrjE1FsojLvKlHLdx6m317/I4+Gtx2gz74fj7+PFU9+eweqpHlje3An6JS3a/MH8K0y7OwOcAh7QWu9w6AEHGRnhDC0Gg8GqxHZLi3XGma0qp42Njfj7+zttzU9beTm529+i+KP3CPnyDH4tRjbe5sWJ0daPv2GLAR+8qZmTwqiVNzJ/+hoiAjzwG3f71FvWNtPUW0utnYbtU28rTdd+Yib0eurt87Nl/PvLh6lu7Kit8x8rUrhv6TjPG3X2Ub+uw1FKhQE+WuuKPj3QICMBZ+jSWlNXV2cJPrW1tcyePdtmldOmpqZ+WfNjbG6mdPenHEhsIfPiHvYU7qGutY6gJs1zfzTg0+mLfc4IReH0kYQvu5q0xbeQHJ7seR+C7RuNtq/5qbrELiARo83Bp3dTb+fK6vneiwc4V9YxDfqNmaP47c1T8ffxsAW+feBWCz9dTabUPE9LSwu+vr49qnLavuYnOjqakJAQp33wtxpaOVRyiJPvvsy8Jz+x264sFE5PDsN70TwmL/82M+Pn4OPlYZlVDk29tWe9XXqhZ1VDC//+8mH25pRbjs0aHcmzd8xiWIjU1ukJCTgOkBGOZ6uurubkyZM0NTXZbePn52cJPtHR0U4LPi0FBZz74J+U79hO2PE8vA22/0bLQuFnP45kQfwi0hPSWTBqAWF+HrjRaH05ZH9kCkBnPu3z1Furwcgv3znGa1/kW47FRwbywrrZpAwP7adfYuiQgOMACTjiUmt+OuvPKqeG2lqKPt1Gwba38d9/DP/GjnTrzyYq/nhjx1SPj/LhaiYyZ+Rc5s29iYRQD1xB39upt9RVpgA0eoHV1JvWmuf3nOOJD07S/paH+vvwp9tmkp4a28+/hHuTgOMACTiiq9bWVqvEg7Y204e/rSqnJSUl1NXVOXXNj25tpfrgfnLe+weGXfv4Z7ovH4yz/jZ/37sGlhzTFEbBmSmRBKQvYsbSbzN1+HTP22hUa9Pmolnb4PSHDk297ThRzAOvf0lDi2khr5eCR9dMZu38pAH6JdyPBJxekGs4oie01tTU1FBeXm4pHtdZ5yqnvr6+lgqnUVFR3TYodfT5tdHI6epsdubvJDM/k5Olx9j8RwOhXWYBawLheGoArVfOYNw13+TKsUs9c6PRnk69KS+In2PZ6+1Eaxx3/+0gF6o7Xtg75o3m0TWT8JHaOt1IwHGAjHCEoy5V5VQpRXh4uOW6T2BgoNOm4ooKsjj/iw0EHj6Nb6vtb/Kt3nA8yYuKtLHE3HAzi1Ov8cyNRtta4PyejgWnVZcorx2ZREPS1fw2ZwyvFSfQiukLw6Lxw/jf264gPNDD1ktdhgQcB0jAEY4yGo2WInO21vx0FhgYyNSpU527y3VTE5V7Msl9/w3U54cIrO6e+GBQcPePvKkPVEyImsCS+CUsTVjKxOiJeCkP+9ZuNfW2HfK/wLRNY3eNXsF80jqFTwxXkGGcTnTsSF5YO5vEaA8cMdohAccBEnCEM9ha89OZvSqnBoPBKWt+tNFIw7Gj5L6/hfqdmYTmmdJ9jyfCY9+1ntpLvqhZdcIfFs1m4lXfZO7ohZ630ShAfZlp6i1rG5z9FFps7wls0IrDejyfe89m6fVrmTZjjksrnA4WEnAcIAFH9Ifm5marEttRUVF2q5yGhYVZpt6Cg4OdMvXWUlBI4fatnPQq5p2EYr4o+oI2oyn54bs7Ddywz/R50OAHR8f6UDsnlfgVN7Bo0irP3Gi0rdmU9daDqbe6oARCpl5n2my0S9abJ5GA0wuSNCAGitFopLW1FX9/6wWFtqqc+vv7W5XY9vZ2TsZZXUsdn1/4nMyCTJZteIe4su6F7gwKTiVA4YxRRCxbzrzZN3rmRqNaQ8lJy4JTXXAAZWfqDf8w6zILQVED21cXkoDjABnhCFc5fvw4paWlds97eXlZbbfTNWA5QmtN3b595H6whZaMzwgqtbd4Egqi4a0bY4hfdA1LE5aSNiINP28P/DZfV0rFV+9zfOc/mNl6mBBlZ5Gw8oKEuR1lFmJSh/TUmwQcB0jAEa7U2tpqmXqrqKiwrPnpqj+qnGqtaTlzhvxtb1O540OCsy/QtVpV5yqnQT5BLBi1gKUx81k4dhmRAZFO7c9gV9vUyk9e3U/Tmd0s8zrM1d6HiVdl9u8QmdRRZiFx/pCbepOA4wAJOGKw0FpbldhuaGiwnLNV5TQ7O5uIiAinrflpKy2lZMd2Lnz4Dv6HTlIWBg+st85kG16p+cNmA8fGKIpnJjLs6lUsmH6dx2w0ajBqnnj/JC98dg7QpKgCbgz6mnXRpwgqOYy9rDf8w2DcMvOC0+VDYupNAo4DJOCIwap9u52mpibGjRtnda5zldOua36ckXptbGyksSCP4yHVpgWnBZnk1+az+gsj6z6xXvtzdgRkT4ogMH0R0xfcSNqI2UO+xs8r+8/zy3eOYzDX1gn09ebPNyRwlfcR07WfszvtZr2hvGDMErjjbbeecpOA4wAJOMIdXarKaWBgoCX4hIeHO6XOj9aac9XnyPvZTxm+85jdduWh8HWKHy1zpzL26m+wMPmqITv1tie7jH9/5RC1TR3ToD9dOYHvL0lGGVogd495r7ftUN0l6238CvjuGwPcY+eSgNMLkqUm3FldXR1lZWU21/x05u3tTUJCAklJSU577pb8fEo+ep/iD98j4FgOXkbbny3vzlG8crUv02Omszh+Menx6YyNGDukpt7OlNRx10sHOF/eMf158xXx/J+bpnTU1tEaSk50LDgtOAjX/gFm3+WiXjuHBBwHyAhHuLvOa34qKiq6ldi2VeXUVm0gRxhqaqjelUne9rdQew/jV9+x28Jjt3lxvEuV0xtPhzNy+pVMX3DjkCmvXVnfwr0vH+KLcx21KeckRfGXO2YRFWwjUaCuFHz8IcC9S0xIwHGABBwxlBiNRqqqqiyJB01NTXarnLa1tTl1zY9ua6Phyy8p2L6V6v2f8/T3R/FlxVGM5p2bA5s0z5urnHaeektediMLx15FVID7XkhvaTPy8NajbDlYYDmWGBXEC+vSGBc7NGvrSMBxgAQcMVRprWloaCAoKOiyVU77Y80PQGVTJbsLd5OZn0nLxzv5wZuN3do0+8DRJEXJzESGLV/J/KnXMi5inNtNvWmt2bQrh99uP9VRWyfAh6e/ewWLxse4tnP9QAKOAyTgCE9TUVHBiRMn7K75AQgJCbEEn9DQUKd8+NcePkjOC3+GLlNvXZ0ZAcdmRsJ31rAkfonbLTj98HgRD75+hMZW024O3l6KjddP5o55o13cM+eSgOMACTjCExmNRkudn65rfjrrjyqnuq2NhsNfUvDhVuoyMggqrOjWpnOV0/YFp4vjF7M4frFbTL0dK6zm7pcOUlTTsSvBuvlJPHztxCFTW8ejA45SKhL4O5ACNALFwA+01mcudT8JOEJgVWK7qqrKUuPHVpXTsrIympqaLHV++qrl/HlKPvqA4o/ft2S9/fF6Lz6bbP3BfM92AxF1UDIzkZirV3Ll1NWMjxg/aKfeimuauOdvB/m6oNpyLD01hj/dOpPQAPdPlvD0gBMBpGmtd5hvPwDcpLVOv9T9JOAIYa2trc1SYnvEiBFERERYne9c5TQoKMgy9RYWFtbnNT+GmhqqMzM4nRpMZtVBMvIzKKgrwMuou1U5PRMHpyeFE7BkITPm38jsuDmDbuqtscXAQ28c4YOjHZu0pgwP4fm1s0mIcu/aOoM+4Cil4oGfAmnAdCAQGKO1zrXRNgF4ClgOKGAH8KDW+hIl+6zunwb8U2uddKl2EnCE6LlLVTn18fGxKrHt69v3b/HtC04PfvoqU3/+it12ZaHwdYovzfOmMXbZDSwccxXRgdF9fn5nMBo1T358mv/d2THZEh3sx0c/Xkx0iHOSM1zBHQJOOvAP4BDgDazARsBRSgUBXwHNwMOYNij6NRAETNNa1/fguV4GyrXWP7pUOwk4QvScwWCguLjYUuen65qfzsLCwpg0aRIBAc4p7taSm0vJx9uspt5safCHe+/3IXXkNNLj01kcv3hQlFl463ABG948SovByL2Lk/nZ6oku7U9fuUPA8dLalJSvlLob2IztgPMj4Ekgtf0ajFJqDJAN/JfW+snLPM+jwDXA1Vpr21dDzSTgCOEYg8FgteanubnZ6rytKqftAcopU2+7MijY9jZ6n3XW24kE2Hi79Wamk1tiWBg2kxkLbmRO3FyXTb0dzK1gy8F8fnPTNLy9Bue1p54a9AGns8sEnE+AAK31gi7HMwG01kuUUncCPzGf2qy1/rO5zcPAGmCF1rqay5CAI0Tfaa2pr6+3BJ+amhqGDRvGlClTrNqVlpZy8uRJp675ac96K/xwK3U7M9iXFsKzU4osC04BvvupgRv2a8vUW9O8qYy96gYWJl/lmRVOnWAoBZwi4B2t9b1djj8N3KK1trmKyjyyWc1lgo1Saj2wHiAxMXGWvU0QhRCOaWlpwWAwdMtis1Xl1NlrfrTBQE1bHbsLd7Mrfxd7Cvfw2J8riS+3btfka1pwWjwzgWFXr2T+lNWDYurNXQylgNMCPKm13tDl+K+BDVrrbsU/lFKTgWPAWaB9X/A2ey9IOxnhCDFwvvzyS6qr7U88+Pn5WRIPIiMjnVLnp6WpnlM/+aFp6q3B/oLT7Dg4PSkM1lzN3CkrmRM3B39v972o398uFXD6/q4Nclrr45iy2Xqk027R/dcpIYSVmTNn0tDQYNlotPOaHzCNjIqKiigqKnJalVO/gGCmPf1XdGureertHeoyMgi8YL3gdPxFGH+xhh+NfYeXiv9FoE8gV8ZdSXpCOoviF8nUWy+4W8CpBGwV0YgynxNCuKmgoCCCgoJISEiwWvNTXl5Oa2urpV1UlPWOAlprcnNziYyMJDw8vNdTX8rXl+C5c0iZOwd+Cc3nzlFmznrzP34OL6PmYiRcND9tY1sjn+Z/yvGvP6H2EyNFMxMYtuwa5k9dTWpkqky9XYK7Tal9CvhprRd2OZ6B6XdZ4qx+yJSaEIOD1pra2lpLxtuECROszneucursNT+Gqiqqd2VwviqXHRNb2VWwi9yaXIBuVU6z4yBrUhgBSxYy/crrmTtynkdOvQ2lazgPAr8HUrTWOeZjSZjSojdorf/grH5IwBHCPVyqymnXEtvOGH3kVueSWZDJ8A1PMzq7xmabsjA4kuJL09wpjLvqBhYlL/OYqTe3CDhKqW+af1wGfB/4AVAKlGqtM81tgjEt/GykY+Hn40AopoWfdoqF96ofUvFTCDdSVVVFSUmJzTU/nQUEBBAfH098fLxTnrdj6u0D/I/bX3Da5At/Xe5F6VXTWJKwhCXxS5gQNWHITr25S8Cx15HMzvueKaUSsd7a5hNMW9vkOrM/MsIRwr3YWvPTla0qp0ajse8LTquqqN6VScGHW2HvYXy7ZL09dqsXx5M6nmN40HBW+80i7YrVzImbS4CPc3ZdGAzcIuAMFjLCEWJoaGlpsSqxbTAYSEtLIyQkxKrdgQMH8PLysky9hYSE9Gn00Z71duHDd6jbuROqarj7AW9avDqu97RXOa0Mga9SfGmcO5lxS01TbzFB7l2UTQKOA2SEI8TQ0V7np2sWm60qp35+fpbgExkZ2ecS221lZTSE+fFZ4Wdk5Gewp3APk7+q4sdbrfebe3mpF/uWjmDHLTvcerrNo9fhCCGEl5dXt3IKALW1tSiluq35uXjxIhcvXkQpZbXdjiMbjvoMG0YYsGrMKlaNWUWbsY1j1b+nNeg1q6m3Q+MUi+IXuXWwuRwZ4XQhU2pCeJa2tjarqbfOa346c3aVU93aSsOhw1z48F+UnzjMc/82kjunrGVx/GKnPL6ryJSaA2RKTQjPo7W2KrFdX99R8cRWldP2AOWsOj9DgUypCSFEDyilCA8PJzw8nOTkZJqamizBZ9iw7uto8vPzqaysRClFWFiY09f8DDUScIQQwo6AgABGjRrFqFGjup1rr/sDppFRdXU11dXV5OTkEBAQYAk+ERERfU67Hiok4HQhm3cKIXpCa01SUpLNNT9NTU0UFhZSWFiIt7c3kZGRjB8/vs81ftydXMOxQ67hCCF6ytaan85sVTlt/+wdalNvcg2nn1RXV1NWVkZLi/1aGkIIEz8/P4YNG0Z4eLiru+J0fn5+jBgxghEjRmA0Gqmurqa8vJyysjKampqIiorqNq1WVlZGdna2U9f8DHYScBzU1NREcXEx8fHxBAYGDrlvKUI4k9aaxsZGCgoK8Pf3d2g9i7vw8vIiMjKSyMhIxo4dS2NjI7ZmksrLy63W/IwdO5aEhAQX9HjgyJUsB5WWlhITEyPZKEL0gFKKoKAghg0bRmlpqau7M2Daf+/g4OBu57pe94mOjh6obrmMjHC66GnSQFNTk1OqDgrhSUJDQykvL3d1NwaF2bNnW9b81NfXExgY6Oou9TsJOF1ord8F3k1LS7vnUu3a2tqcUlddCE/i4+NDW1ubq7sxKHRe8+MpZEqtD2QqTYjekb8ZzyYBRwghxICQgCOEEGJASMARFi+++CJKKc6cOWPz/Lp160hKSnLosd99911uu+02UlJS8PLyIj093fGOim4u9945m7yfwhEScLpQSq1RSm2qrq52dVcGnUceeYS3337boftu3bqVI0eOMG/ePKfVlBeuI++ncISkWXXR0yw1TzR27FiH77t582bLSuuFCxc6q0vCReT9FI6QEY7osa5Tarm5uSilePbZZ/nlL39JXFwcERERrFmzhoKCAqv7ym65rtXa2srDDz9MUlISfn5+JCUl8fDDD3crNpaTk8Pq1asJCgoiNjaWhx56iE2bNqGUIjc319JO3k/hCBnhOFnShvdd3QWL3N9eOyDP85vf/Ib58+fzwgsvUFJSwkMPPcTtt99ORkbGgDy/M0x9aaqru2BxdO1Rpz/m2rVr2bJlCz//+c9ZuHAhn3/+OU888QQ5OTm8+uqrgGkDyuXLl9Pc3MwzzzxDTEwMzz33HP/85z+d3h/hmSTgiD5LSkqyfGiBaduf//zP/+TChQuMHDnShT0TAMeOHeO1117j0UcfZePGjQCsWLECHx8fHnnkETZs2MC0adN48cUXycnJYf/+/cyZMweAVatWMWPGDPLy8lz4G4ihQsbFos9Wr15tdXvqVNNoQT6kBoddu3YBcPvtt1sdb7+dmZkJwL59+0hMTLQEGzAt1Lz55psHqKdiqBvyIxyl1D+AiYABaAV+prX+pL+eb6CmsQaTqKgoq9vtRaaamppc0R2H9Mc01mBRUVEBQFxcnNXx9r0A289fvHiR2NjYbvcfPnx4P/dQeApPGOHcq7WeprWeCdwLvKGU8oTfWwig4wtBUVGR1fH22+3n4+LiKCkp6Xb/4uLifu6h8BQu/+BVSsUrpf6klNqrlGpQSmmlVJKdtglKqX8qpaqVUjVKqbeUUomXenytdVWnm56zS54QZosXLwbg9ddftzr+yiuvAFgWbc6bN4+8vDy++OILSxutNW+++ebAdFQMeYNhSm0c8C3gELAbWGGrkVIqCPgUaAbWAhr4NbBTKTVNa11v7wmUUk8BN2AKODdrrY1O/Q2GmO3bt3crvdDXHW3Pnz/PgQMHAFPhKS8vL0v20+zZsxk9enSfHl+Y2Hvvbr31VjZu3EhbWxvz589n7969PP7449x6662Wa27r1q3jd7/7HTfddBNPPPGEJUutsrISsE6FlvdTOERr7dL/AK9OP9+NKZAk2Wj3I0zXYcZ1OjYGaAN+0sPnWgkcAPwu13bWrFn6Uk6cOHHJ8+7or3/9qza//t3+mzx5sl67dq0ePXq0pf25c+c0oDdv3mz1ODt37tSA3rlzZ48e+69//evA/IJD2OXeu+bmZv2LX/xCJyYmah8fH52YmKh/8Ytf6JaWFqvHOXPmjF61apUOCAjQw4YN0w888ID+7W9/qwFdVVXVo+e73Ps5FP92RAfgoLbzuaq0jdKnrqKUuhvYDIzRWud2OfcJEKC1XtDleCaA1nqJUupO4CfmU5u11n+28RxngG9rrQ9dqi9paWn64MGDds+fPHmSiRMnXv6XEsLNXXfddZw8eZKzZ8865fHkb2doU0od0lqn2To3GKbUemoy8I6N48eBWwC01n8D/tZ+QikVCIzQWp8z374SiAZybD2BUmo9sB4gMfGSl4aEGJKefPJJQkJCGD9+PLW1tbzxxhu8//77PPPMM67umhgC3CngRAGVNo5XAJF27hMIvKqUCsU09VaP6RqOrcdBa70J2ASmEU6feyyEm/H39+epp54iLy8Pg8FAamoqzz33HHfddZeruyaGAHcKOL2mta4AruzNfZRSa4A148aN659OCTGI3Xfffdx3332u7oYYolyeFt0Lldgeydgb+ThEa/2u1nq9J9UZF0KIgeBOAec4pus4XU0CTjjrSaQejhBC9A93Cjj/AuYppZLbD5gXiC4wn3MKGeEIIUT/GBTXcJRS3zT/OMv8/1VKqVKgVGudaT62Gfgh8I5S6mFMOf+PA/nAs07si1zDEUKIfjAoAg7wRpfbT5v/nwmkA2it65VSVwFPAX8HFPAJ8KDWus5ZHdFS8VMIIfrFoAg4WmvVw3Z5QL/ulS4jHCGE6B/udA1nQMg1HCGE6B8ScITFiy++iFKKM2fO2Dy/bt06kpKSev24NTU1/OpXv2L+/PlER0cTERHB/Pnz2bp1ax97LNpd7r1zJnk/haMk4HQhadH2PfLII7z99tu9vl9eXh5PP/00S5YsoqF68AAACmVJREFU4eWXX+Yf//gHKSkpfOMb3+DPf+623Z0Y5OT9FI4aFNdwBhNJGrBv7NixDt1vzJgx5OTkEBQUZDl2zTXXkJ+fz+9+9ztZ2e5m5P0UjpIRjuixrlNqubm5KKV49tln+eUvf0lcXBwRERGsWbOGgoICS7vg4GCrD6d2aWlpXLhwYSC67vFaW1t5+OGHSUpKws/Pj6SkJB5++GFaW1ut2uXk5LB69WqCgoKIjY3loYceYtOmTSilyM3NBeT9FI6TEY6zbRxEyQYbB2Za8De/+Q3z58/nhRdeoKSkhIceeojbb7+djIyMS95v165dTJgwYUD62BOlf/pfyno4JRRxyy3EPf4rq2MXH/klVW90zfC3bdh99xFz/w+tjjUeO07gFFubafTd2rVr2bJlCz//+c9ZuHAhn3/+OU888QQ5OTm8+uqrALS0tLB8+XKam5t55plnLAXY2gurXc5gez/F4CMBpwtJi+69pKQky4cWQGlpKf/5n//JhQsXGDlypM37bNq0iX379vHyyy8PVDc91rFjx3jttdd49NFH2bhxIwArVqzAx8eHRx55hA0bNjBt2jRefPFFcnJy2L9/P3PmzAFg1apVzJgxg7y8vEs+h7yfoidkSq0LSYvuvdWrV1vdbi9ZbO9DKiMjgwceeIA777yT7373u/3eP0+3a9cuAG6//Xar4+23MzNNm3ns27ePxMRES7ABUEpx882XXvom76foKRnhONsATWMNJlFRUVa3/f39AWhqaurW9sCBA1x//fVcddVVPPfccwPSv56Kuf+H3aa5eiPu8V91m2brjf6aTquoqAAgLi7O6viIESOszl+8eJHY2Nhu9x8+fLjdxx7M76cYfGSEIwbM0aNHueaaa5gxYwZvvvkmvr6+ru6SR2j/QlBUVGR1vP12+/m4uDhKSkq63b+4uNjm48r7KXpLAo4YENnZ2Sxfvpzk5GTee+89AgMDXd0lj7F48WIAXn/9davjr7zyCgDp6ekAzJs3j7y8PL744gtLG601b775ZrfHlPdTOEKm1LqQpAHYvn27ZbqlXV+uaZWUlLB8+XJaWlp47LHHOHHCunzRzJkzLdNwom/svXe33norGzdupK2tjfnz57N3714ef/xxbr31Vss1t3Xr1vG73/2Om266iSeeeMKSpVZZaapv6OVl+n4q76dwlAScLmThJ9x///3djk2ePJm0tDSHHu/EiROcP38egOuuu67b+XPnzjm0ZY7ozt57d/jwYZKTk3nhhRf49a9/zciRI/npT3/Ko48+amnn5+fHRx99xP3338/3v/99QkJCuO2225g7dy4bNmywfOmQ91M4SmmtXd2HQSktLU0fPHjQ7vmTJ08yceLEAeyREK5x3XXXcfLkSc6ePeuUx5O/naFNKXVIa23z26mMcIQQFk8++SQhISGMHz+e2tpa3njjDd5//32eeeYZV3dNDAEScIQQFv7+/jz11FPk5eVhMBhITU3lueee46677nJ118QQIAFHCGFx3333yeabot9IWrQQQogBIQGni97Uw5GECyF6R/5mPJsEnC56upear68vjY2NA9QrIYaGxsZG2ZHAg0nAcVBsbCyFhYU0NDTItzYhLkNrTUNDA4WFhTb3axOeQZIGHBQWFgbAhQsXuhWxEkJ05+vry/Dhwy1/O8LzSMDpg7CwMPnjEUKIHvKYKTWl1L8ppbRS6kZX90UIITyRRwQcpVQScA+wz7U9EUIIz+XygKOUildK/UkptVcp1WAehSTZaZuglPqnUqpaKVWjlHpLKZV4mcf3Ap4D7geanf4LCCGE6BGXBxxgHPAtoBLYba+RUioI+BSYAKwF7gDGAzuVUsGXePyfAJ9prQ85rcdCCCF6bTAkDezSWg///+3dS6hVZRiH8edFaxJdzEYRFVIGGkqNjMqSSjQoByXNbFIWNnHYwEiIGkSZ40MjJw0rlQZdUCND0kGWp0FGORAKjHIQigW+Dfa2Dvucrfuy9rfWdj0/OHBct+/lXL6/a6+zvxcgIl4A1vc57kVgGXBPZv7UPf474CTwErCr94SIuBd4Blg7gbolSUOo/Q4nMy8OeOjTwJFLYdM99xfgMLAJICK2RMS33Y9XgIeBO4GTEXEKWAPMRMTojeslSSNpwh3OoFYCHy+wfRbYDJCZe4A9Pfv/W1c9Ig4CuzPzo4UGiIitwNbuPy9ExIkxa26jG4ErrwvULHXXXGL8qseo6nrjXGeUc4c95xbg9yHHaLu7++2YpsC5mc5znl5/AEuqGCAzZ4AZgIg41q+JkPqLiJnM3HrlI5uj7ppLjF/1GFVdb5zrjHLusOc4DwwvImb67ZumwBlbZj5adw0tsK/uAkZQd80lxq96jKquN851Rjm37u91G/T9Gtf+DGcIf7LwnUy/Ox/VIDOn7he67ppLjF/1GFVdb5zrjHJu3d/rNrjc13iaAmeWznOcXiuAHyYwXt/bQkmt4TxQoWkKnL3AmohYdmlD9w2iD3b3Var7PEdSizkPVCuasLR+RDzb/fQx4GVgG3AGOJOZh7rHXAccB84DO4AE3gCuB1Zl5l+l65YkDa4pgdOviENzH/R3l7F5D3gCCOALYHtmnpp0jZKk8TQicCRJV79peobTKBHxWkT8GBEXbXkgtUtELImI/d054HhEfBoRd9VdV9MZOKP7DNgAfFl3IZKKSzqrlizPzNXAfjqr0usyWhM4VbdByMwjmflzidolja/KOSAzz2bm53NO+ZrOuo26jNYEDpNvgyCp2SY5B2xn4bUeNUeblraZWBsESVNhInNARLzePX6q1hCsQ2vucKpsgyBp+kxiDoiIHcCTwMbMPFdVrVer1gTOEFYCC7UlmKWzjI6kq9tAc0D3zuYpYH1mTltLjloYOPMN1AYhInZGxGngAeD9iDgdEbcVqlHS5FxxDoiIlcBOYClwqNv08VixCqdUm57hVCozd9L5gZPUMpk5S2e1Ew3BO5z5bIMgtZtzwIQYOPOVboMgqVmcAybEwJmvaBsESY3jHDAhrVq80zYIUrs5B9SrbYFjGwSpxZwD6tWqwJEk1cdnOJKkIgwcSVIRBo4kqQgDR5JUhIEjSSrCwJEkFWHgSJKKMHAkSUUYOJKkIgwcSVIRBo4kqQgDR5JUhIEjSSpicd0FSLq8iNgEPA6sBp6n0+p4c3f3Q8BbmflJTeVJA7M9gdRgEXEt8HZmbo+Io8AF4ENgV2ZmRLwKbMvM22stVBqAL6lJzfYI8FVEBLAM+C0z383//6cYwNLaqpOG4EtqUrOdAM4Cq+i8lLa7Z/9qOu2QpcbzDkdqsMz8NTPPA+uAc8A3l/ZFxDXABmBfTeVJQzFwpOmwDjicmX/P2bYRuAH4ICIWRcQd9ZQmDcbAkRouIhYBa4EDPbu2AAcy8xSdv2K7v3Bp0lAMHKn57gNuAg72bF8O7I2IxcBzwP7CdUlD8Y8GpOa7FfgeONqz/U0678tZAbyTmf+ULkwahu/DkSQV4UtqkqQiDBxJUhEGjiSpCANHklSEgSNJKsLAkSQVYeBIkoowcCRJRRg4kqQi/gXJ9ueixsCNOwAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.plot(n, err[:,0], label='Lin1')\n",
"plt.plot(n, err[:,1], label='Lin2')\n",
"plt.plot(n, err[:,2], label='Log1')\n",
"plt.plot(n, err[:,3], '--', label='Log2')\n",
"\n",
"plt.plot(n, 50.0/n**2, '--', c='gray', alpha=0.5)\n",
"\n",
"plt.xlabel(r'$n$')\n",
"plt.ylabel(r'$\\epsilon$')\n",
"\n",
"plt.title('$f(x) = 1/(1+ax^2)$')\n",
"plt.xscale('log')\n",
"plt.yscale('log')\n",
"\n",
"plt.xlim(10, 200)\n",
"plt.ylim(1e-4, 50)\n",
"\n",
"plt.legend(ncol=2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this case both the `Log` methods outperform the `Lin` methods, and essentially converge. There is some funny business with `Lin1`, but in general it is not the best method."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let us work the the second problem. Recall that for $G(t) = g_1 \\exp(-t/\\tau_1) + g_2 \\exp(-t/\\tau_2)$, the integral,\n",
"\n",
"$$I(T) = \\int_{0}^{T} G(t)\\, dt = g_1 \\tau_1 \\left(1 - e^{-T/\\tau_1}\\right) + g_2 \\tau_2 \\left(1 - e^{-T/\\tau_2}\\right)$$\n",
"\n",
"Thus,\n",
"$$I = \\int_{t_\\min}^{t_\\max} G(t)\\, dt = I(t_\\max) - I(t_\\min)$$"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"def G(t, g=[0.5, 0.5], tau=[0.01, 100]):\n",
" G = np.zeros(len(t))\n",
" for i in range(len(g)):\n",
" G += g[i] * np.exp(-t/tau[i])\n",
" return G\n",
"\n",
"def Iexact(T=1000.0, g=[0.5, 0.5], tau=[0.01, 100]):\n",
" Iexact = 0.0\n",
" for i in range(len(g)):\n",
" Iexact += g[i] * tau[i] * (1.0 - np.exp(-T/tau[i]))\n",
" return Iexact"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"tmin = 1e-4\n",
"tmax = 1e+4\n",
"g = [0.5, 0.5]\n",
"tau = [0.01, 100]\n",
"\n",
"n = 100\n",
"ti = np.geomspace(tmin, tmax, n)\n",
"Gi = G(ti, g, tau)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAvwAAAEoCAYAAAAkHyUrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdeZxcdZX38c+pql7Tnc7WWUjSCdkgIQSQABEEAgICCkFBXFF8ZNMZHx2dcXRGEZWZZ3x01EeRcVDHBXHDYURxYwv7EgNIIIEskJAEQtLZ093ppbrO88e91WmaTnd1uqpuLd/361WvTt2699aphl/1qV+de37m7oiIiIiISGmKRR2AiIiIiIjkjhJ+EREREZESpoRfRERERKSEKeEXERERESlhSvhFREREREqYEn4RERERkRKmhF9ERIbNzE41s9XFcl4RyS0z+66ZfT7qOCSghF8iY2YbzOysHJ5/sZltztX5RcpVf2PX3R909yOyec5snFeknEX5d9bdr3H3L+fquWVolPCLiIiIiJQwJfxSUMysysy+aWavhLdvmllVr8c/bWZbwseuMDM3s1mH8Dxzzew+M9ttZivN7MJej401s9+Z2V4z+4uZXW9mD2XrNYqUor4zfeHM4t+b2Qoz22NmvzSz6kI5r0i5yuPf2R+Z2fXhvxeb2WYz+5SZbQvP/6Fsvi4ZmBJ+KTT/DCwCjgWOAU4EPgdgZucCnwTOAmYBiw/lCcysAvgdcCcwHvgYcIuZpcsGvgO0AhOBD4Y3ERm6S4FzgcOBBcDlBX5ekXKQ87+zBzERaAAmAx8GvmNmo7N4fhmAEn4pNO8DvuTu29y9GfgicFn42KXAD919pbu3Adcd4nMsAuqAf3P3Tne/F7gDeI+ZxYGLgS+4e5u7rwJ+PIzXI1LOvuXur7j7ToIP2ccW+HlFykE+/s72pyt83i53/wPQAuj6nDxRwi+F5jDgpV73Xwq3pR/b1Ouxnn+bWZOZtaRvGTzHJndP9XmeyUAjkDjY84jIkLza699tBB+0MbM/9hqv78vWeUUkI/n4O9ufHe6e7HVfYzePElEHINLHK8A0YGV4vyncBrAFmNJr36npf7j7RjJ/43gFmGpmsV5JfxOwBmgGkuHzrOn7PCIyfO5+XtQxiJSxfPydlQKjGX6JWoWZVadvwM+Bz5lZo5mNA64Ffhru+yvgQ+EFt7VARv19e58/fI5lBDMLnzazCjNbDFwA/MLdu4HbgOvMrNbMjgQ+kM0XLFIi+o7dbEwgVfQZr5qUEhm+vP+dNTPLySuRQ6aEX6L2B2B/r1s1sBxYATwDPAlcD+DufwS+BSwF1gGPhefoGOD8k/ucfz/BjMUFwHnAduBG4APu/nx4zN8SXFj0KnAzwZvjQM8hUo76jt3rCvScIuUuir+zM7P9ImR4zN2jjkHkkJjZXOBZoKpPXWC2n+crwER3V7ceEREpG/n6Oyu5pxl+KSpm9vawh/Bo4CvA77L9JmRmR5rZAgucSNA+7H+y+RwiIiKFKB9/ZyX/lPBLsbka2Aa8AHQDH8nBc9QT1PG3Ar8E/h24PQfPIyIiUmjy8XdW8kwlPSIiIiIiJUwz/CIiIiIiJUwJv4iIiIhICVOP4ywZN26cT58+PeowRArGE088sd3dG6OO42A0ZkVeq5DHrMaryGsNdbwq4c+S6dOns3z58qjDECkYZvbS4HtFR2NW5LUKecxqvIq81lDHq0p6RERERERKmBJ+EYmMmU0xs2+b2aNm1mZmbmbTMzw2ZmafNbMNZtZuZk+b2cW5jVikfGm8ihQvJfwiEqVZwKXALuDBIR77ZeA64AbgPIIl4G81s/OzGaCI9NB4FSlSquEXkSg94O4TAMzsCuCcTA4ys/HA3wP/5u5fCzcvNbNZwL8Bf8hFsCJlTuNVpEhphl9EIuPuqUM89C1AJfDTPtt/ChxtZocPKzAReR2NV5HipYS/H2b2eTNbY2YpM7so6nhE5HWOAjqAdX22rwx/zstvOCIyAI1XkYippKd/dwG3AP+VtROu2sqLzS09980OPGbY67YBxMwwAwMs/W8L9o6ZEbMD+8RjRsyMeOzALRH+rIjHSMSMikSMyniMqkSMykSMqkSc6ooY1RVxqhIxrG8AIoVrDLDb3b3P9p29Hj9k67a1cM9zW5nRWMfZ8yYM51QikuPxCvCdpesYWVPBpJHVTBlTQ9OYWmorleKIpBXlaDCzKcA/AguBY4Aa4HB339DPvlOBbwBnE+TOdwOfcPeNBzu/uz8WHpu1mH/z1Mv8/pktWTtftplBTUWc2soEI6qCn/VVCeqqE9RVJaivDv5dXxXcHxHeairj1IQfGKoScSriRiIeIx5+EInFLPzAcuCDTX/PPWh8GWw0rOdc6YfS/w17YujzIar3h6bgQ1SwTR9+SpOZXQVcBdDU1HTQ/Z5/dS//54/P89ajJynhF4lIpuO1M5nia3eupu/HiUkN1cwaX8e8SSOZP7mBhdNHM6mhJpchixSsokz4OdAp4AmCTgH9XjhkZrXAvQRfJX4QcOB6gouFFrh7a37ChTfPHc/k0cEbTe9JjvQ//XX3HfdgXw+3p7elPNijOxXc73YnlXJSvf7d1e2k3OnqTpHsDn52pZyuZIrO7hQdyW7au1J0dAU/O7tTtHV209bZzfYDX0SUrZhBIhY78E1J3EjEYlTGe31TUnHgW5Kail4flKoTjKxO0FBTwajaSsaOqKSxvorx9dWMrEnow0R27AJGmZn1mTVMzxTu7OcY3P0m4CaAhQsX9p1tFJHcyOl47U45n3jzHLbs2c8re9rZvLONTbva2LKnnS172nlw7faefaeOqeHU2Y28+cjxvGn2OKoS8WG/OJFiUKwJf6adAq4EZgBHuPu6cP8VwFrgauDreYgVgHe8YUq+nuqQJLtT7O8KEv7WjiStHd20dCTZ195Fa2eSlvYk+zqS7GsP/t3amaSto5u2rm72dybpTKboSKaCDxip4MNGKgUp954PK/153Re8/e2T0XHez4cn77nf8+HJ0/eDiNLxpX92+4EPVZ3dKegePL6hGFEZZ/LoGqaNHcGMxhHMGV/P3EkjmTOhjkRcl9QMwUqgCpjJa+uC07XAq/IekYgcTE7Ha01lnI+fNfs127pTzqadbazeuo9Vr+zlr5t28+TGXWzauZ+fPb6Rnz2+kfrqBG9bMIn3nNjEgimjhhOCSMEryoR/CJ0CLgQeSyf74bHrzexhYAnwdTP7APDJ8OHvuft3shttcUjEY9THY9RXV0QdSuTcg29PusOfXd3Bz2T4LUlnMhV+wOmmI5miPfygtL8z/SEpyd72Lvbs72JXayc7WjvZvq+DrXvbae3sZs3WFtZsfe3XKDUVcY6Z2sAbZ4zj9CMaWTC5gVhM3wQM4E9AF/A+4Iu9tr8feNbd10cSlYj0J+/jNR4zpo8bwfRxI3jLUROB4EPAsy/v4b7Vzfxp5as8t2UvP1+2iZ8v28SJ08fw0TNmcvqcRn0LKyWpKBP+ITgKuL2f7SuBdwK4+0+An+QzKClsZkYiblkfHO7O3v1JNu1qY8OOVl7Y1srqrXtZ+cpeXtrRxmMv7uSxF3fyjbvXML6+ircumMQlx0/hqMMashxJYTGzS8J/Hh/+PM/MmoFmd78/3CcJ/NjdPwzg7tvM7OvAZ81sH/Ak8C7gTIIP+iKSA8U8XuMx45ipozhm6ig+ftZs1m3bxy+WbeJXyzexbMNOlv1wJyfPHMvn3jqPeYeNzFdYInlR6gn/GILawb52AqMPdpCZXQdcATQC883sBmCRu2/us19GFxSJQPBBoqG2gobaBuZPfm0Sv7O1k2Xrd/Lwuu3c89xWXtnTzg8f3sAPH97AwmmjueLUGbzlqAmlOvN0a5/7N4Y/7wcWh/+Oh7fe/hloAT4OTARWA5e6+x25CVNEKKHxOmt8PZ972zw+ftZsbnl8I/9x3ws88sIOLrjhIT66eCYfO3M2lQmVWkppKPWE/5C4+3UES4APtp8uAJSsGDOiknPnT+Tc+RP50pKjWLF5D7c9uZnbnnqZ5S/tYvlLT7BgSgOfPW8ub5w5Nupws8rdB/0U098+7t5NcBH+9bmIS0RerxTHa311BdecPpN3nzCVb969lh8/uoFv37uOpau38d33H8+U0bVRhygybKX+0XUX/c/kH2zmXyRyZsHXzl9cMp/HPvtmvnDBPBrrq1ixeQ/v+d5jfPrXT7OnrSvqMEVESsqo2kquu/AofnX1G5k6poZnX97LhTc8zOMv7og6NJFhK/WEfyVBHX9f81AXDykCI6oSfOiUw7n/Hxbzd2fNoTIe41fLN3P+tx5kxebdUYcnIlJyTpg+hjv+9lROm9PIztZOPvBfy7j3+a1RhyUyLKWe8P8WWGRmM9IbzGw6cEr4mEhRqK1M8PGzZvOHj5/KMVMaeHn3fi75j0f5zVMvRx2aiEjJaait4IeXn8B7T2qiI5niqp88wZ9Xvhp1WCKHrGgTfjO7JOwW0LtTwCVmdnqv3b4HbABuN7MlZnYhQdeeTcB/5jVgkSyYNb6OX13zRt6/qInO7hSf+OVf+cmjG6IOS0Sk5MRjxr9cNJ+rT5tBMuV87OdP8ZcN/a4RJlLwijbhJ+gUcCtwTXj/xvB+T4/fcCXdM4E1wM3ALcB64Ex313qyUpSqEnGuv+ho/vn8uQBce/tKfvzIhmiDEhEpQWbGZ847kvee1ERnMsWHf/QXXmxW+iDFp2gTfne3g9wW99lvo7tf7O4j3b3e3S9y9w3RRC2SPVeeNoN/eft8AK773Ur+9OyWiCMSESk9ZsaXl8znrLkT2Nue5KO3PMn+ziwvwy6SY0Wb8IsIvO+kafz9OXNwh4//4q88+/KeqEMSESk58ZjxjXcdw4xxI3j+1X1ce/uzUYckMiRK+EWK3N+cMYt3LZxKRzLFx37+FC0dyahDEhEpOfXVFdz4/jdQXRHj1ic261tVKSpK+EWKnJnxxSVHceTEetZvb9XMk4hIjhw5cSSfPS+4fupzv1nJ7rbOiCMSyYwSfpESUF0R54b3HkdNRZzbnnyZpau3RR2SiEhJumzRNE6YPprtLR18+Y7nog5HJCNK+EVKxKzx9Xzy7DkAfOH2lbR36aIyEZFsi8WMr1y8gMpEjP9+cjNPbtwVdUgig1LCL1JCLj9lOnMm1LFxZxs33vdC1OGIiJSkGY11XPGmwwG4/o5VuHvEEYkMTAm/SAmpiMe4/qKjAbjpgRdo3tcRcUQiIqXpI4tnMq6ukic37ub3z+gCXilsSvhFSsyJh4/hrLkTaO9K8d37NcsvIpIL9dUVfPLsIwD42p9Xk+xORRyRyMEp4RcpQZ84azYAP33sJbbtbY84GhGR0nTpwik0jallw4427lihWX4pXEr4RUrQ/MkNvOWoCXQkU3z3/hejDkdEpCQl4jE+ungmADcsXUcqpVp+KUxK+EVK1MfODGb5b12+SYtxiYjkyDveMIXJo2pYt62FPz77atThiPRLCb9IiZo/uYGF00azryPJ/zz1ctThiIiUpMpEjKtPnwHADx7SN6pSmJTwi5SwD5w8HYCbH92gtnEiIjlyyfFTqK9O8OTG3azYvDvqcEReRwm/SAk796iJNNZXsWZrC4+9uDPqcERESlJtZYJ3nzAVgB89siHaYET6oYRfpIRVJmK8J/wj9OsnNkccjYhI6bps0XTM4I6nt7C9RWugSGFRwi9S4i46bjIAf175Ku1d3RFHIyJSmprG1vLmI8fT2Z3SBIsUHCX8IiVuRmMdR09uoKUjydLnt0UdjohIyXr3CU1A0B1N101JIVHCL1IGLjzmMAB++/QrEUciIlK6Fh/RyLi6Kl5obuWpTbp4VwqHEn6RMvC2YyZhBvc8v4197V1RhyMiUpIS8RjveENQRnnrcpX1SOFQwi9SBiY11HDi9DF0JlPcq7IeEZGcueT4KQDc8fQr7O/UdVNSGJTwi5SJs+dNAOC+1c0RRyIiUrrmTKhnwZQG9nUkWbpaEyxSGJTwi5SJxUeMB+D+Nc10pwrjYjIzm2pmvzazPWa218xuM7OmDI9tMrMfm9lGM9tvZmvM7HozG5HruEXKlcZsZi5YEFw3dccKXTclhUEJv0iZmNk4gqljatjZ2lkQK0GaWS1wL3Ak8EHgMmA2sHSwBCB8/G7gNODzwPnA94FPAf+Vw7BFypbGbObeumASAPc+v43WjmTE0YhAIuoARCQ/zIwzjhjPTx59iaWrmzmuaXTUIV0JzACOcPd1AGa2AlgLXA18fYBjTyFINN7i7neG25aa2Rjg782s1t3bche6SFnSmM3QYaNqOH7aaJ54aRd3P7eVJcdOjjokKXOa4RcpI2eEZT33FUZd6YXAY+nEAcDd1wMPA0sGObYy/Lm3z/bdBO9rlq0gRaSHxuwQvC2c5b9jxZaIIxFRwi9SVhbNGEtVIsaKzXvYEf3S70cBz/azfSUwb5Bj7yaYVfyKmc0zszozOxP4OPBdd2/NbqgigsbskLz16KAd8v1rmlXWI5Ery4TfzD4fXiyUMrOL+jw208weCh9/yswWRhWnSLbVVMY5flpQyvOXDTsjjoYxwK5+tu8EBqw3cvd24E0E72ErgX3APcAdwN8e7Dgzu8rMlpvZ8uZmdSsSGaK8jtliH6/jR1bzhqbRdCZTPLCm+OKX0lKWCT9wF3Au8EA/j30X+LG7zwE+DdxiZiX3VaOUrxMPHwPA4+sjT/gPmZlVA78ExhNcOHg68A/Au4DvHOw4d7/J3Re6+8LGxsa8xCoihzZmS2G8ptsh37lqa8SRSLkryITfzKaY2bfN7FEzazMzN7PpB9l3yC3C3P0xd3+xn3M1AouAH4X73UVQV3j8MF+SSME4cXqQ8BfADP8u+p8VPNgsYm8fBhYD57v7T939AXf/GkHHj2vM7JisRioioDE7ZOmE/97nt9HVnYo4GilnBZnwA7OASwneQB482E7DaRF2EE3AFnfv6rVtQ7hdpCQc1zSaRMxY9cpe9rV3DX5A7qwkqAnuax6wapBjjwZ2ufsLfbYvC3/OHWZsIvJ6GrNDNLOxjpmNI9izv6sQJlmkjBVqwv+Au09w9/OBWwfYL90i7CJ3/427307QRWAaQYswEemjpjLO0VMaSDk88dJgk3I59VtgkZnNSG8Iv8k7JXxsIK8Co81sVp/tJ4U/X85SjCJygMbsITh73kQA7lJZj0SoIBN+d8/0e69BW4SZ2QfM7K/h7W8GOd9GYJKZVfTaNj3cLlIy0nX8y6Kt4/8ewTdot5vZEjO7ELgd2AT8Z3onM5tmZkkzu7bXsT8iuOjvD2b2QTM7w8z+Afga8ATBe4CIZJfG7CFIl/Xc89w23AtjlXMpPwWZ8A/BoC3C3P0n7n5seDvoxXzhvs0EXy9eDmBmZxPU8D+RzaBFonZSAST8YRu+M4E1wM3ALcB64Ex3b+m1qwFxer1fufsGgutt/gpcD/yB4Bu/m4CzhzBpICIZ0pg9NMdOHcXo2go27mxj/faS6z4qRaLYV9o9pBZhZnYdcAXQCMw3sxuARe6+GbgG+HE489AGvM8P8pHczK4CrgJoalKZvxSP45uChH/Fy3vo6k5REY/ms7+7bwQuHmSfDfSzKI+7ryK41kdE8kRjdujiMeO0OY3c/tdXWLq6mRmNdVGHJGWo2Gf4D4m7X+fuU9y9yt3Hhf/eHD621t1Pdvc54bcCywY4T9G3DJPy1FBbwbSxtXQmU6zd2jL4ASIicsgWHxHkCAWyyrmUoWJP+IfTIkykrM2f3ADAsy/viTgSEZHSdtrsRszg8Rd30tapVXcl/4o94R9OizCRsnZ0mPA/o4RfRCSnxtZVsWDKKDq7Uzz6wo6ow5EyVOwJ/3BahImUtXTCv0IJv4hIzp0RlvUsVVmPRKBgE34zu8TMLuHAKrfnhdtO77VbRi3CROT15h8WJPzPbdmrFSBFRHLstDlBwv/Q2u0RRyLlqJC79PRdcOvG8Of9BMtz4+6tZnYm8A2CFmEG3AN8ok+LMBHpI33h7ks72li7tYV5h42MOiQRkZK1YHIDI6sTbNjRxsYdbTSNrY06JCkjBTvD7+52kNviPvttdPeL3X2ku9e7+0VhSzARGYQu3BURyY9EPMYps8YB8OC65oijkXJTsAm/iOSeLtwVEcmfU2cHZT0PrlFZj+SXEn6RMpau41+1ZW/EkYiIlL5TZwcz/A+/sJ2krp2SPFLCL1LG5kwIVnxcu3UfB1lQWkREsmTqmFoOHzeCfe1Jnt6sb1Ylf5Twi5SxxvoqRlYn2NuepHlfR9ThiIiUvPQs/4NrVccv+aOEX6SMmRlzJtQDsHabGluJiOTam8ILdx9ZpwW4JH+U8IuUudlhWc+arfsijkREpPQtmjmWmMGTG3fR0pGMOhwpE0r4Rcrc7PGa4RcRyZeR1RUcM3UUyZSzbL1m+SU/lPCLlLnZvS7cFRGR3EuX9Ty0Vgm/5IcSfpEyl67hX7O1RZ16RETyIJ3wP7xO/fglP5Twi5S58fVV1Fcn2LO/i+0tnVGHIyJS8o5rGk1NRZzVW/exbW971OFIGVDCL1LmXtOpR2U9IiI5V5mIcdKMMQA8pFl+yQMl/CLC7PFhHb8u3BURyYueOn4l/JIHSvhFhFlhwv9CsxJ+EZF8OHV2IwAPrd2u66ck55TwiwhNY2oBeGlHW8SRiIiUhzkT6misr2Lbvg59uyo5p4RfRJg2dgQAG3cq4RcRyQcz6ynreXCtynokt5Twi0jPDP/mXW10p/TVsohIPhzox98ccSRS6pTwiwg1lXEmjKyiq9t5Zff+qMMRESkLb5odJPyPvbiTjmR3xNFIKVPCLyIATBsTlPWojl9EJD8mjKzmyIn17O/q5okNu6IOR0qYEn4RAaBpbHjh7s7WiCMRESkfp88JuvXcv0ZlPZI7SvhFBIBpYR3/Rs3wi4jkTTrhv2+1En7JHSX8IgLAtHH5L+kxs6lm9msz22Nme83sNjNrGsLxc83sVjPbbmb7zWy1mX08lzGLlDON2exbOH0MtZVxVm/dx5Y9uoZKckMJv4gAB2b4X8pTa04zqwXuBY4EPghcBswGlprZiAyOXwg8DlQBVwDnA/8OxHMVs0g505jNjcpEjJNnBhfvPqCyHsmRRNQBiEhhmJau4d/RirtjZrl+yiuBGcAR7r4OwMxWAGuBq4GvH+xAM4sBPwHucfe393poae7CFSl7GrM5cvoRjdz93FbuW93Mu07I+AsTkYxphl9EABhVW8nI6gRtnd1sb+nMx1NeCDyWThwA3H098DCwZJBjFwNzGSDBEJGs05jNkcVhHf+Da7erPafkhBJ+EekxfVx6xd28dOo5Cni2n+0rgXmDHPum8Ge1mT1mZl1mts3MvmVmNVmNUkTSNGZzZOqYWuZOGklLR5JHX9gRdThSgpTwi0iPqelOPfmp4x8D9Nd4eicwepBjDwt//hK4Ezgb+L8EdcE/O9hBZnaVmS03s+XNzaqVFRmivI7ZchuvZ8+bAMCdq7ZGHImUorJM+M3s82a2xsxSZnZRn8fODd9gVoSzEMdEFadIvk0eFUy0vbK7PeJIBpV+7/qpu1/r7ve5+9eALwIXmdnc/g5y95vcfaG7L2xsbMxbsCIy9DFbbuP1nDDhv3vVVlIpjzgaKTVlmfADdwHnAg/03mhmo4FbgA+6+wLgH8L7ImVhUkM1QL5aw+2i/1nBg80i9pb+zvuuPtvvDH8eN4y4RKR/GrM5dNRhI5k8qoZt+zp4evPuqMORElOQCb+ZTTGzb5vZo2bWZmZuZtMPsu+QewK7+2Pu/mI/D80Edrj7ynC/B4EmM3vDMF+SSFGY1BDM8G/Jzwz/SoKa4L7mAasyOHYgqUOKSEQGojGbQ2bWU9bz55Uq65HsOqSE38yOMrNPmdnNZvaIma00s+fN7HEz+4WZ/aOZHT+MuGYBlxLMGDw4QBzD6gncj7XAWDM7JTz/hUA9MP0QziVSdA4bFczwv7InLwn/b4FFZjYjvSH8YH9K+NhA/gh0AG/ps/3c8Ofy7IQoIr1ozObYufMnAvC7p19RWY9kVcYJvwXeZ2YrCVpwnQVsBX4P3AB8A7gdeJngavy7wxX0PmJmQ11U4wF3n+Du5wO3DrBfuifwRe7+G3e/naBt2DSCnsBD4u57gEuAfzWzJ4BzCGYtkkM9l0gxOqynhj8vJT3fAzYAt5vZkvAD9u3AJuA/0zuZ2TQzS5rZtelt7r4D+D/ANWb2r2Z2lpl9BrgW+HHvtoEikjUaszl24vQxTB5Vw8u797Nsw86ow5ESktHCW2Z2OHAzsBH4X8Bf3H3Ar98sWLVnIfAx4Gozu8zdn8nk+QY7dy/99gQ2s3RP4K+b2QeAT4YPf8/dvzPIcy8lXAjEzKqAVxn8q0qRkjB2RCWViRh79nfR1pmktjJ3a/O5e6uZnUkwWXAzYMA9wCfcvaXXrkawEmffCYovAfuAjwJ/D2wBvgp8OWdBi5Qxjdnci8WMJccexo33vcD/PPkyi2aMjTokKRGD/jU3s2MJWmddPpRP4O7uwF+AD4QfGG40s6+4+32HGmw/jiKYXehrJfDOMI6fEKzulxEzm+TuW8K7nwfuPdjrNrOrgKsAmpq0Mp4UPzNjUkM1L+1o45Xd7cwaX5fT53P3jcDFg+yzgSCB6LvdCRbx0UI+InmiMZt773jDZG687wX+8MwWvrjkKKorhlokIfJ6mZT0XARcMJyv28KV+C4AzgiX186WQ+oJbGbXmdlm4I3A981ss5lNCR/+Ung9wjqC2v0PH+w85dYyTMpDnjv1iIhIL7PG13P05Ab2dST588pXow5HSsSgybe7X+fuHcN9IndPuvsXhlCukzPha5ri7lXuPi789+bwsSvd/Uh3n+Xu73d39caSsnJYfjv1iIhIH5eeMBWA/3poPcEXIyLDc8iz7Wb2kX62TTGzCcMLaUiG0xNYRPqRvnD35fxcuCsiIn1c8oYpjK6t4OnNe/jLBqUzMnxD6dIz1cwqe226oJ/d9hFcoPv3w44sM8PpCSwi/Zg0SiU9IiJRqqmMc9miaQDc9EB/ywaJDM1QZvgfA3ab2VIz+yJQ1+cDAO6+x92/BNxnZp/KZqAHMZyewJg/TC0AACAASURBVCLSj56Snvz04hcRkX584OTpVCVi3P3cVp7bsjfqcKTIDSXhPxX4J4Llsa8m6LW/x8weMLMvhz11awHcfTkwrNIeM7vEzC4B0gt4nRduO73Xbhn1BBaRzKVn+PPUi19ERPoxrq6K95wYdAD8/G+e1UJcMiwZJ/zu/qK7f9PdL3H3icAy4OMEC21dAdwJ7ApX2/0lQQec4bg1vF0T3r8xvP/FXjG1AmcCawh6At8CrAfO7NMTWEQylK7h37KnXReLiYhE6JPnzKGxvorlL+3i1ic2RR2OFLHhtMjcFbalfI+7TwLmA58iSLhjwN8MJzB3t4PcFvfZb6O7X+zuI9293t0vCnsAi8ghGFldQV1VgrbObvbs74o6HBGRsjWyuoLPv20eAP/y++dYt21fxBFJsRpOwv/D3nfcfZW73+Du73b3d7r7imHGJiIRmRj24n91r+r4RUSidMGCSZwzbwJ725O8//vL2LyrLeqQpAgdcsLv7r/KZiAiUjga66oAaN437CU4ipqjkiYRiZaZ8f/efRwnTh/Dq3vbefuNj/DHZ7ao5FKGZMCE38wSZnZ5tp7MzD6erXOJSO6MH1neCb9hUYcgItKjpjLO9y9fyImHj6F5XwcfueVJzv/WQ3zrnrXc+/xW1m7dx9a97ext76K9q5uu7pQu8pXXSAz0oLsnzazFzL4JfMbdD+n7fTMbBXyfoKuOiBQ4zfCLiBSWkdUV/OLKRdyybCNf/dPzPLdl76DtOs0gETPiMSMRi5GIGxXxGJXxGBVxozIRoyoRp7oiRnVFnJqKOHVVCWqr4tRXV9BQU8GomgpG1VYyrq6S8fXVjB9ZRXVFPE+vWrJlwIQfwN1/bWY7gQfM7KfAze6e0bJvZjYJ+ARwHnCluz8+rGhFJC8a64OEf1uZJvyxcIK/WzNkIlJAYjHjskXTuHThFB5au5371zSzfnsrL+/az972JPs7k3SlnO7w5g5d3U5XtwOprMUxrq6SKaNrOXzcCGY2jmDOhHrmT25gUkM1ZvqGtBANmvADuPu9ZnYWQR/+dWa2HngEeAbYHd5iwJjwNg84naAX/3eAN4YtNEWkCJR7SU88zPiV8ItIIapKxHnz3Am8ee7Blzxyd1IOyVSKVAq6UimS3U5Xd4rOZCr42Z2ioytFe1c37ckUbR1JWju7aWnvoqUjyd72JLtaO9nV1klzSyfb93WwdW8721s62d7SyV837X7NczbWV3HC9NG8ceY4zjiikSmja3P9q5AMZZTwA7j7XuAzZvYl4K3A2cBVwHSgAXCCxH898BDBzP4D7t6Z5ZhFJMca64IuPeWa8FfEg8ubglkxEZHiY2bEDeKxoPymhuyU4XSnnG372tm4o40NO1pZu7WF51/dxzMv76F5Xwd/eOZV/vDMqwDMnTSSJccexkXHTu7p/ibRyDjhT3P3Ng4siiUiJShd0tPcUp4JfyIezPAnU9n7ClxEpBTEY8akhhomNdRw0oyxPdvdnRe3t/L4izt5cG0zD6xp7rnO4Kt/Xs258ydy5akzOHbqqAijL1+DJvwWFGNNcvdX8hCPiBSAnhr+Mu3Dny7p0Qy/iEhmzIyZjXXMbKzjvSc10ZHs5v7Vzdz25Mvc9dxWfr9iC79fsYWz5o7nH889ktkT6qMOuaxkMsMfA24ws7nA3cA9wFJ335PewczOcPelOYpRRPJsVE0FFXFjb3uS9q7usuvIkC7pUQ2/iMihqUrEOeeoiZxz1ES27NnPjx7ZwE8eeYm7n9vGfaubufK0GXz8zbPL7u9LVAZdeMvdu4FLgKeAjwLXA78zs7N77XZ9bsITkSjEYsa4sDXn9jIs60mEM/zJbpX0iIgM16SGGj573lwe+PQZvPekJrrd+Y/7XmDJDQ+zduu+qMMrC5mutDsdqARmu/t8dz/N3e/q9fh4M1tmZteZ2elmVpH1SEUkr8q5Nacu2hURyb7G+ir+9e1H898fOZkZ40aweus+LrjhIf74zJaoQyt5mSb8nyHoo//iAPu8AbgWuBfYZWbnDDc4EYnO+Prybc2pi3ZFRHLnDU2j+d3H3sQ73jCZ9q4UH/3Zk3z/wYFSTBmuTBP+kYMstrULqAPOAn4E3Obudw4zNhGJUGM5J/yx4K0xqRl+EZGcGFGV4N/feQz/eO6RuMP1v3+OG+9bF3VYJSvThH/MII/vdvd2d7/X3T8MrDGz04YZm4hEqLGunEt6wi49muEXEckZM+Mji2fy7+88BjP4v39azQ8eWh91WCUp04S/2swG6ujz0T73vw5cfkgRiUhBKOsZ/nQNf1Iz/CIiuXbx8VP4yjsWAPDlO1Zx58pXI46o9GSa8D8CXHqwB919XZ/7bcDIYcQlIhFrrC/f1XYrYqrhFxHJp0tPmMo/nnskAJ/45V95bsveiCMqLZkm/D8C/tXMhrI8WsPQwxGRQlHOq+2qS4+ISP5dc/oM3n7cZNo6u/noLU/S2pGMOqSSkVHC7+7PA38C7jCzQWfuw/Kf8cOMTUQi1NOlJ4er7ZrZVDP7tZntMbO9ZnabmTUdwnk+Y2ZuZg9lI66eLj3qwy/yGoU6ZqU0mBn/5x1Hc+TEetZvb+X636+KOqSSkekMP8CnCMp0/mJmpw6y75XA84cclYhELj3Dv72lE/fsz3SbWS1BG98jgQ8ClwGzgaVmNmII55kBfA7Ylq3Yemb4tdKuSI9CHrNSOqor4nzz3cdSmYjx82WbuGvV1qhDKgkZJ/zu3gq8FUgC95nZnWZ2qZk1pvcxsxFm9knga8A3sh6tiORNdUWcmoo4nd0pWju7c/EUVwIzgIvc/TfufjtwITANuHoI5/kP4BbguWwFdqCkRzP8Ir0U7JiV0nLkxJF8+i1HAPCF25+lrVOlPcM1lBl+3H0TcDLwM4Ke+z8HXjWz3Wb2MkE//q8A/+Tuj2U7WBHJrzEjKgHY1dqZi9NfCDzW+6J/d18PPAwsyeQEZvZegkX/PpvNwOIxI2bgDt2a5RdJK9gxK6XnQ6cczvzJI3llTzvfukf9+YdrSAk/gLvvcffLgBOB7wMvAJXhuX4PnOHu/y+rUYpIJEaPqABgZ24S/qOAZ/vZvhKYN9jBZjaa4JvET7v7zizHpll+kdcr6DErpSUeM768ZD5m8P0HX2TdtpaoQypqQ07409x9ubtf7e5z3L3W3Se5+9vdXRfgiJSI0bXBDP/Otpwk/GMIvhXsaycwOoPjvwqsIegilnXphL8jqYRfJFTQY1ZKz3FNo3nXwqkkU87X71oddThF7ZATfhEpfWNzW9JzyMLGAR8APuJDuKLYzK4ys+Vmtry5uXnAfSvUqUckaw5lzA5lvErp+sRZc6hKxPjDM6/yzOY9UYdTtMou4Tez0WZ2h5mtMbOnw4uPZ/V6/Hwze9LM/mpmK83smijjFYnS6DDhz1FJzy76nxU82Cxib/8J/ADYbGajwjVCEkA8vF/V30HufpO7L3T3hY2Njf3t0iM9w59UDb9IWl7H7FDGq5SuiQ3VfPDk6QB89U7N8h+qskv4AQe+GZYiHQPcQXAtAmYWI7gg+XJ3PxY4G/h3MzsssmhFIjSmNqcJ/0qCmuC+5gGDNV+eC1xDkGSkb6cAi8J/f2S4waUT/k6V9IikFfSYldL1kdNnUleV4IE1zTy5cbDPltKfgkz4zWyKmX3bzB41s7ZwcY7pB9l3SIuAuPtud7+716ZHgPS5LfyZXlG4HtgH6EoRKUvpGf5duanh/y2wKOzJDUA4zk8JHxvIGf3cnia4oPAM4NfDDa4yESb8KukRSSvoMSula/SISi574zQAbrr/xYijKU4FmfADs4BLCT71P3iwnbK0CMgngNsB3L0beCfw32b2EvAEcI277z3E1yFS1MbktqTne8AG4HYzW2JmFxKMxU0EX/8DYGbTzCxpZtemt7n7fX1vwG5gT3h/83CDS9fwq0uPSI+CHrNS2j508nQq4zH+vOpV1m9vjTqcolOoCf8D7j7B3c8Hbh1gv2EtAmJmXwiP/2x4PwH8M3Cpu08j+KrxPw5l2XCRUpDu0rOrtSvr5w4X8zuToGvHzQQL8awHznT33t+qGRAnz+9X6Rn+rqRq+EWg8MeslLbxI6u56LjDcA/adMrQJKIOoD/unumUWr+LgJhZehGQr5vZB4BPhg9/z92/A2BmnwPOB85x97bw8WOBw9x9aXiuZ83sWeAkYONwX5dIsemZ4c9NSQ/uvhG4eJB9NnCg3G6g/RZnJ6pATw1/d05WGRYpSoU8ZqX0XXXaDH61fDO/fmIzn37LkTTUVkQdUtEo9k/fgy4C4u4/cfdjw1s62f8CcAFBst+7x9Mm4DAzmx/uNwVYEJ5PpOykF94qtLac+XDgol3N8IuIFIJZ4+t506xxdCRT3PaUqsCGotgT/iEvAmJmRwHXAWOB+8P2m8sB3H0rcAXwMzN7GvgT8E/u3m/3AfUIllLXU9LT1kmqzNpTViW00q6ISKF570lBlfXPHt/IEJZhKXsFWdKTS+6+kgG+anT3XwC/yPBcNwE3ASxcuFD/10nJqYjHGFmdYG97kj37u3q69pQDteUUESk8Z8+bwLi6KtZua2H5S7s4YfqYqEMqCsU+wz+cRUBEJAO5ruMvVJVxzfCLiBSainiMSxdOAYJZfslMsSf8w1kEREQy0NOLv8zq+CvUh19EpCC9+4SgrOdPz75Ka0cy4miKQ7En/MNZBEREMpDj1XYLVnqGv0MlPSIiBaVpbC0Lp41mf1c3d656NepwikLBJvxmdomZXQIcH246L9x2eq/dMloEREQOXY5X2y1YPSvtKuEXESk4S46bDMD/PPVKxJEUh4JN+AkW3LoVuCa8f2N4/4vpHYawCIiIHKIDq+1mf/GtQlalhF9EpGC97ehJJGLGQ2ub2bavPepwCl7BJvzubge5Le6z30Z3v9jdR7p7vbtfFC76ISJZ0Ls1ZzmpVA2/iEjBGj2iksVHjCfl8Lunt0QdTsEr2IRfRArD2HCGf0dLmSX86Rr+LiX8IiKFaMmxhwHwh2eU8A9GCb+IDKhca/h7Snq6uyOORERE+nPGkeOpTMR44qVdbNursp6BKOEXkQGNrq0AyjDhr1ANv4hIIaurSnDa7EYA/rxS3XoGooRfRAY0siZI+PfsL6+LdtWWU0Sk8J03fyIAf3xWCf9AlPCLyIAawoR/7/7yWtykqiIOqIZfRKSQnTV3AomY8fj6nWW3XsxQKOEXkQEdSPi7cPeIo8mfdA1/R1I1/CIihaqhtoKTZ42jO+XcvWpr1OEULCX8IjKg6oo4lYkYnd0p2stotrsqEc7wq6RHRKSgnT1vAgD3PK+E/2CU8IvIoEZWh7P87eVTx39ghl8Jv4hIITvzyPEAPLh2u76VPQgl/CIyqIaaBFBeF+5WhzX87V364yEiUsgmj6ph7qSRtHV28/iLO6MOpyAp4ReRQTWUYaeedFtOzfCLiBS+N4ez/Pc+vy3iSAqTEn4RGVRPwt9WRgm/LtoVESkaZ84NEv57nt9aVg0mMqWEX0QGle7FX041/AdKejTDLyJS6I6ZMoqxIyrZtHM/67a1RB1OwVHCLyKDKseSnuqEavhFRIpFPGacPidYdff+Nc0RR1N4lPCLyKDKMuEPa/g1wy8iUhxOU8J/UEr4RWRQZZnwV2qGX0SkmJw6exxm8Pj6nezv1Ht3b0r4RWRQPX349ycjjiR/atSWU0SkqIytq2L+YQ10JlM8vn5H1OEUFCX8IjKokWU4w18RjxGPGcmU09Wtsh4RkWJw2pxxgMp6+lLCLyKDSpf07C2jhB80yy8iUmxOnxO051TC/1pK+EVkUA05astpZlPN7NdmtsfM9prZbWbWlMFxC83sJjN73szazGyjmd1iZodnM750a879SvhFgMIfsyLHNY2irirBi82tvLx7f9ThFAwl/CIyqJE1CSC7JT1mVgvcCxwJfBC4DJgNLDWzEYMc/m7gKOBbwHnAZ4A3AMvNbGq2YqypDN4idfGXSHGMWZGKeIxFM8YC8PDa7RFHUzgSUQcgIoUvR116rgRmAEe4+zoAM1sBrAWuBr4+wLFfcffXfF9rZg8D68PzXpuNAGsrgrfI1g4l/CIUwZgVgaBbz93PbeXBddu59AR9ngTN8ItIBuqqEsRjRltndzYvYL0QeCydOAC4+3rgYWDJQAf2TRzCbS8BzcDkbAVYU5ku6Smf7kQiAyj4MSsC8KbZwYW7D6/bTirlEUdTGJTwi8igzIyR1cFsdxYv3D0KeLaf7SuBeUM9mZnNBcYDzw0zrh61YcLfppIeESiCMSsCMGPcCA5rqGZnayfPvbo36nAKghJ+EclIDlpzjgF29bN9JzB6KCcyswTwXYLZwh8MsN9VZrbczJY3Nw/ewUEJv8hr5HXMDnW8iqSZWc8s/0Oq4weU8ItIhgp8td0bgJOB97t7fwkJAO5+k7svdPeFjY2Ng560pjL4VkMX7Ypk3aBjdqjjVaS3U2aFCf86JfyghF9EMpSDhH8X/c8KHmwWsV9m9m/AVcD/cvc7sxQbACPCGf7WTtXwi1AEY1Yk7U1hwr9s/U6tpUKZJvxm9kszW2FmT5nZMjN7c7g9bmZ/7XVbaWZuZguijlkkaiOr0734s5b8riSoCe5rHrAqkxOY2T8D/wj8b3e/OVuBpdVVBTP8Ldl7zSLFrODHrEja2Loq5k0aSUcyxZMvZfx5tGSVZcIPXO3uC9z9OIJWYreaWczdu9392PQNuB5Y4e4rog1XJHo5qOH/LbDIzGakN5jZdOCU8LEBmdn/Jhij/+zuN2QrqN5GVKXbcirhF6EIxqxIb+k6/gdV1lOYCb+ZTTGzb5vZo+GKfB6+qfS375BX/XP33b3uNgyw6xXA94f8AkRKUM9qu9lL+L8HbABuN7MlZnYhcDuwCfjP9E5mNs3MkmZ2ba9t7wa+CfwJuNfMFvW6DblbyMH0zPCrD78IFMGYFektXcf/sBL+gl14axZwKfAE8CBwTn879Vr1r4Ng1T8nmD1YamYL3L31YE9gZt8g6BvcAFzs7qk+j88E3ghcMuxXI1ICsl3D7+6tZnYm8A3gZsCAe4BPuHtLr10NiPPaCYpzw+3nhrfe7gcWZyNGzfCLHFAMY1aktxOnj6EyHuOZl/ewu62TUbWVUYcUmUJN+B9w9wkAZnYFB0n4Gcaqf+7+d8Dfmdm5wP81s1PcvbPXLh8G/nugjh8i5aQ+7MO/L4v17O6+Ebh4kH02ECQKvbddDlyetUAOoi58zS26aFcEKPwxK9JbTWWc46eN5tEXd/DICzs4/+hJUYcUmYIs6ek72z6AQVf9M7MP9LoI92/6ea4/EXQdODq9zcziBN8YqJxHJJRO+FvKaLa7riro0qOLdkVEilNPHX+Z9+MvyIR/CAZd9c/df9LrQtzvmFmNmR2e3tHM3giMBV7sdfxbgRZ3vz+HsYsUlQMdawqyD39OpDsT7Suj1ywiUkrS7TkfXNuMu0ccTXQKtaQnU4ey6l8N8DMzqweSQCtBDX/v81zBAKt1ppnZVQS9hGlqGvA6YZGid+AC1vKZ7U53JspiK1IREcmj+ZMbGFVbweZd+3lpRxvTx42IOqRIFHvCP2TuvpPgYtyB9rkww3PdBNwEsHDhwvL92ChloS4HNfyFrmftgcJcXVhERAYRjxlvmjWOO1Zs4YG1zWWb8Bd7SU9WVv0TkcHVVwXJbzmtOtvTilQlPSIiReu02Y0APLCmfOv4iz3hH/aqfyKSmZ6ONWU0w19dEaMibrR3pehIqhe/iEgxOnVOUMf/6Avb6erOtC9MaSn2hH9Yq/6JSOZGpDvWdCTL5sInM+tV1lM+H3RERErJpIYaZo+vo7Wzm6c27h78gBJUsAm/mV1iZpcAx4ebzgu3nd5rt4xW/ROR4atKxKmMx+jqdjqS5TNDcmDBsc5B9hQRkUJ1aljWc9/qbRFHEo2CTfiBW8PbNeH9G8P7X0zvEK6keyawhmDVv1uA9cCZfVb9E5EsqCvDXvxjRgQrM+5oUcIvIlKszjxyPAD3Pl+eCX/Bdulxdxt8r8xW/ROR7KirSrCztZOW9iTj6qqiDicv0gn/zlYl/CIixerEw8cwojLO86/u4+Xd+5k8qibqkPKqkGf4RaTAlGMv/rF14Qy/En4RkaJVmYj1lPUsLcNZfiX8IpKxcuzFr5IeEZHSUM5lPUr4RSRj9WU4wz9mRFC6tLO1I+JIRERkOBYfGczwP7xuO/s7y6vVshJ+EclYeoa/tYwS/nFhSc92zfCLiBS18fXVLJjSQEcyxUPrymsRLiX8IpKxdA3/vjJK+CeOrAZgy579EUciIiLDde78iQD84ZktEUeSX0r4RSRj5bja7mFhJ4cte9ojjkRERIbrrUdPAuCuVVtp7yqfsh4l/CKSsbrKdA1/V8SR5M+EkdWYwda97STLdEl2EZFSMW3sCOZPHklLR5IH15ZPWY8SfhHJWDnO8FcmYoyrqyLlsG2fLtwVESl254ez/L9f8UrEkeSPEn4RyVg51vADHNagOn4RkVKRLuu5+7ltZdOtRwm/iGSsvgxn+AGaxo4AYP32togjERGR4Zo2dgTHNY2ipSPJHWUyy6+EX0QyVldVAZRXH36A2ePrAFi7bV/EkYiISDa858QmAH62bGPEkeSHEn4RyVhPDX+ZJvzrtrZEHImIiGTDBQsOo746wVMbd7Pqlb1Rh5NzSvhFJGN1ZbjSLsDsCekZfiX8IiKloKYyztuPmwzAz5a9FHE0uaeEX0QyVq41/NPGjqAyHmPTrjZ2t2nFXRGRUvC+k6YB8OsnNrNtX2mvtaKEX0QyVq4z/BXxGMc1jcIdHntxZ9ThiIhIFhwxsZ5z5k2gvSvFDfeuizqcnFLCLyIZq62MYwZtnd10pzzqcPLq5JnjAHjsxR0RRyIiItnyqXOOwAx+vmwjm3aWbic2JfwikjEz67XabnnN8r9x5lgA7l/TjHt5fdgRESlVR0ys56JjJ9PV7Vz/+1Ul+/6uhF9EhiSbnXrMbKqZ/drM9pjZXjO7zcyaMjy22sy+amZbzGy/mT1qZqcNO6iDOK5pFOPrq1i/vbWslmMX6a2YxqxIpj51zhzqqxL8eeVWbn6sNC/gVcIvIkPSU8c/zAt3zawWuBc4EvggcBkwG1hqZiMyOMUPgCuBa4G3AVuAP5vZscMK7CAq4jE+ePJ0AG564MWSnQUSOZhiG7MimZoyupZ/u3gBANff8RyPl2DpphJ+ERmSAzP8XcM91ZXADOAid/+Nu98OXAhMA64e6EAzOwZ4L/B37v49d78HuBTYCHxpuIEdzHtPbGJEZZyH1m3nxvteyNXTiBSqohuzIpl664JJXLZoGp3dKS77wTJue3JzSU3sKOEXkSFJz/DvG35rzguBx9y9pzWCu68HHgaWZHBsF/DLXscmgV8AbzGzquEG15/RIyr5+ruOxQy++ufVXPmT5Ty4tpk9+4f94UekGBTdmBUZii9cMI8PvjFI+j/5q6e5+D8e4TdPvcwru/dHHdqwJaIOQESKS7oXf2tH93BPdRRwez/bVwLvzODY9e7et6XCSqASmBX+O+vectRE/uWio/nyHau4a9VW7lq1FQg+CDXUVFBTGScRMxJxwzDMwADMchGOSEYWThvN5982b7inKcoxK5KpRDzGF5fMZ87Eer7659U8uXE3T278KxB0qRszopLayjgV8RjxmOX0vf1XVy+iKhHP2vmU8IvIkBzoxT/sWe0xwK5+tu8ERg/j2PTjr2NmVwFXATQ1ZXSdYb/ee1ITZ80dzw8eXs+y9TtZ+cpeWjqSZde5SIrH2BGV2ThNXsdstsaryFC976RpXHTsZG5dvol7Vzfz9Kbd7NnfRVtn/mb6s11NpIRfRIbkmtNn8t6TpjFtTG3UoQyZu98E3ASwcOHCYb2djh9ZzWfPmwtAKuXs60iyp62L9mQ3yW4nmUrhDh4877BjFxmOkTUVUYcwZNkcryJDNaIqweWnHM7lpxyOu9Pa2c2Olg72d3XTlXRSHtxypTKe3ap7JfwiMiQzGuuydapd9D8reLCZwL7HTjvIsXBg1jAvYjGjoaaChiJMqkSGoGTGrMhQmBl1VYmeb7iLkS7aFZGorCSo6+1rHrAqg2MPD9sE9j22EyjtNdJFoqExK1KklPCLSFR+CywysxnpDWY2HTglfGwgvwMq6HWhoJklgHcBd7p7R7aDFRGNWZFipYRfRKLyPWADcLuZLTGzCwk6gGwC/jO9k5lNM7OkmV2b3ubuTxG09/ummV1hZm8maO93OPCFPL4GkXKiMStSpJTwi0gk3L0VOBNYA9wM3AKsB85095ZeuxoQ5/XvVx8CfghcD/wemAqc6+5P5jh0kbKkMStSvIr36gMRKXruvhG4eJB9NhC2su+zfT/wyfAmInmgMStSnDTDLyIiIiJSwpTwi4iIiIiUMNOCMNlhZs3AS1HHkQPjgO1RB1Eg9Ls4IJPfxTR3b8xHMIciwzGr/+a5p99xfhT1mNXf2LKg38UBWR+vSvhlQGa23N0XRh1HIdDv4oBy+V2Uy+uMkn7H+aHfc2HSf5cD9Ls4IBe/C5X0iIiIiIiUMCX8IiIiIiIlTAm/DOamqAMoIPpdHFAuv4tyeZ1R0u84P/R7Lkz673KAfhcHZP13oRp+EREREZESphl+EREREZESpoRfRERERKSEKeGXYTOz0WZ2h5mtMbOnzexOM5sVdVz5ZmafD38HKTO7KOp4CoGZfcjMvJh+H2Y21cx+bWZ7zGyvmd1mZk1Rx1VIzGxx+N+17213n/1Gm9n3zWy7mbWa2d1mdnQ/56s2s6+a2RYz229mj5rZaf3sFzOzz5rZBjNrD99vLs7la80XM5tiZt8OX3tb+Puc3s9+Wf9d2f9v715j5CrLAI7/H9tCBC9YEEQBy0UxrVoMkECA0CJSIXJRQT6gGA3ljhDURExFDIkfDCAgYkA0hOHY5QAACRtJREFUUSFKECOFREARisGCKIWGGrkkLchNSMo1QCny+OE9w4zTmWXbzu7ZOfP/JSe7e+bZmfc8szvPs+++c07Ewoj4V0SsiYgHIuLEPnFHRMSy6v4eiYhFETFtY49d/Vlf26yx61qfGmvDr0FI4MLM/HBmzgVuAK6oeUx1+CPwaeD2ugcyFVTNykLgznpHMn4RsRnwZ+AjwJeBLwEfAm6NiM3rHNsU9TVg747twNYNERHA9ZTfidOAzwMzKLncrut+fkb5WTkb+AzwJHBTROzWFXcucA5wCXAw5Wfrmog4ZKBHVY9dgC8AzwJ/GSNuoLmKiIXAZcC1lOfqGuDSiDipK25BFXN3dX8XAYuA76/ncWr9WF/brLEd1rvGZqbbCGzAdsCPgKXAy5QXkVl9YrcHfgs8D7wA/A7YYT0eaw9gVd3HXFdOgNuAI+o+vjpzQZlM+BOw+zDlAzgd+C+wS8e+HYHXgTPrHt9U2YB51c/IgWPEHF7FzO/Y925gNXBxx765VdxXOvZNBx4AFnfs2xpYA3yv63FuAZbXnZMB5PRtHZ8f1+t3cNC5qr73aeAXXXE/p1zlc0bHvmXAkq64s4HXgPfVnb+anzvr6yTmZZhqykTlYkNqrDP8o2Ncs0cDmuE8A7huo0Y7OSYzJ1PdoHNxJnBHZv5jwkY8MQ4D7szMh1s7MnMlcAelgdX4HQY8kZm3tnZk5vOUWf/Du+LWAld3xL0O/AZYEBGbVrsXAJsAV3Y9zpXAxyJix4EfwSTKzDfGETboXO0NvLdH3K+ALYF9oSxzA3brEzeDMuM/yqyvvVlj22qvsdM3YNAaTrdn5jYAEXEccFCfuIXATsCuraYnIpYDDwEnABeM9SAR8d3q+48f0Lgn0qTkZEgMLBcR8VHK8o111hUPgTn0LqYrgKMmeSzD4KqI2Ap4DrgJ+FZmPlrdNge4v8f3rACOjYh3ZOZLVdzKzHy5R9wmlEK5oopbAzzcIw5gNrByI49nqht0ruZUX3c/T51xt/aLy8yVEfFyFTfKrK+9WWPbaq+xzvCPiHHOHsE4Zjgj4tiIuLfaTmnFRcQi4BDg4B4FacoZZE6G3YBzsR8wC3goIlYBewGXR8SpAxvwxJlJmYHpthp4zySPZSp7HjifsvTkAMp68QOBpRGxdRUzVi6hnc+3ipvZ8fG5rP6fPUZckw06V62P3fc53rjWvlHIfV/W196ssW1Tocba8KvbWLNyswEy85eZuVu1/RjenHk4FDio+rd9k7xlTkbIeH4+fpKZ22bmrMycRXlD0fGZecnkDVMTKTOXZeY3MvP6zFySmRdS3ky3DeWNvJLWZX3tzRrbNmE11oZf3dZ7hjMi5lDOBrElsKSamfj7hI1w8o0rJxFxTkQ8RlkXe0VEPNbjbCTDrukz4M/S+zj6HbcqmXkP8CCwZ7VrrFy2bh9P3OqOuC2qs/+MFddkg85V6znovs/xxrX2jULuB8H62ps1tm3Caqxr+LXRMnMF0F1YRk5mnkN5YVaHzJxX9xjWQ2v9c7fZwD8neSzDqrWMZAW916nOBh6t1u+34j4bEZt1LVWYTTkDzMMdcZsCO/P/a9NbM4Cj8PwMOlettfpzKKf3HE/c0lZQdVrAzRiN3NfC+tpmje1tvDXWGX51c4ZzXeakrem5WAzsFRE7tXZUTc0+1W3qIyL2AHYF/lbtWgx8ICL274h5F2VpQmcur6ec6eWojrjpwNHAzZm5ptp9I+UMNcd0PfQXgfurda5NN+hcLaWcfrNX3GrKumGqN2Lf1yduLfCHDT+kkdL0188NZV7aJiwXzvCrmzOc6zInbU3PxU+BU4HrqjfJJeUNqf+mXJxIQERcRTnLyz2UM/R8AjgLeBy4uApbTGkor4yIb1KK1VmU2coftO4rM5dFxNXAhRExo7rfkyjXPzimI+7piLgAOCsiXqwe+2jKm4YPm7ijnTwRcWT16e7Vx4Mj4hngmeq9EgPNVWaujYjvUC609TjlvN4HAF8FTsvM1zqG923ghoi4DPg15TlfBFyUmU8NNhON1fTXzw1lXtomLhfre/EAt+Hf6HNRl+q2MygXGdqpY98syizO1+seuzkxF5Nw/DtQrij6AvAi8PteuRjljdK4L6ecrWct5Q+iy4Ftu+JmUi7itJpysZlbgLk97u/tlFPvPQW8CtwFzOsRN43SZD5COe3kcuDIuvMxwLxmn+22icwV5XR/D1ZxDwEn94n7HGWmfw3wKOXCW9PqzttU2kb99dO8TN1cRHVnGgEds0efBE4ETgbenD2qYjanvKC/QikWrRnOdwIfz/a620YwJ23mQpI2jK+fvZmXtrpzYcM/QiKi35O9JDve9BEROwA/BD5F+ff7LcAZmblqosc42cxJm7mQpA3j62dv5qWt7lzY8EuSJEkN5ll6JEmSpAaz4ZckSZIazIZfkiRJajAbfkmSJKnBbPglSZKkBrPhlyRJkhrMhl+SVJuI2DkinoiID9Y9FklqKht+SVKdDgVmAv+peyCS1FQ2/JKkOu0H3JWZr9Y9EElqKht+NY5LBKShsi9we92DkPTWrK/Dy4ZfTeQSAWkKi4ijI+LGiLgL2BqYX319St1jkzQm6+uQisysewzSQEXEtcBWmbl/3WOR1F9EnABcDGyRma/UPR5JY7O+Di9n+NVELhGQhsN84G6bfWloWF+HlA2/GsElAtJQmgcsqXsQkvqzvjaDS3rUKC4RkIZDRMwB7gcWZObNdY9H0tisr8PNGX41jUsEpOEwH3gd+CtARGwREdvXOyRJY7C+DjEbfjXNPFwiIA2D/YB7M/Ol6uvTKX8ASJqa5mF9HVo2/GqMaonANviCJA2DacAqgIjYE3glM5+sdUSSerK+Dj8bfjWJSwSk4XEu8P6IOI8yc3hevcORNAbr65CbXvcApAHqtUTg8hrHI6mPzLwP2KfucUgaF+vrkHOGX03iEgFJkgbP+jrkPC2nGiMi5gKXAkspl/0+PzPfqHdUkiQNN+vr8LPhlyRJkhrMJT2SJElSg9nwS5IkSQ1mwy9JkiQ1mA2/JEmS1GA2/JIkSVKD2fBLkiRJDWbDL0mSJDWYDb8kSZLUYDb8kiRJUoPZ8EuSJEkN9j/C26846DvprQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 864x288 with 3 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(1,3, figsize=(12,4))\n",
"ax[0].plot(ti, Gi, lw=2)\n",
"ax[0].set_xscale('log'); ax[0].set_yscale('log')\n",
"ax[0].set_title('Log-Log')\n",
"ax[0].set_xlabel('$t$'); \n",
"ax[0].set_ylabel('$G(t)$'); \n",
"\n",
"ax[1].plot(ti, Gi, lw=2)\n",
"ax[1].set_title('Lin-Lin')\n",
"ax[1].set_xlabel('$t$');\n",
"\n",
"ax[2].plot(ti, Gi, lw=2)\n",
"ax[2].set_xscale('log');\n",
"ax[2].set_title('Log-Lin');\n",
"ax[2].set_xlabel('$t$');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note the steep fall in $G(t)$ at large $t$. In a certain sense, we should expect this to cause trouble for `LogTrapz1`. On the other hand the last plot looks reasonably smooth (Log-Lin). Thus, we should expect `LogTrapz2` to work ok. "
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"50.004900249193746\n"
]
}
],
"source": [
"print(Iexact(tmax, g, tau) - Iexact(tmin, g, tau))"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(49.99994987558861, 5.385833148352503e-07)"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# notice some loss of precision already\n",
"quadrature(G, tmin, tmax)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/usr/lib/python3/dist-packages/ipykernel_launcher.py:39: RuntimeWarning: overflow encountered in exp\n",
"/usr/lib/python3/dist-packages/ipykernel_launcher.py:39: RuntimeWarning: invalid value encountered in double_scalars\n"
]
}
],
"source": [
"tmin = 1e-4\n",
"n = 2**np.arange(4,8)\n",
"err = np.zeros((len(n),4))\n",
"Iex = Iexact(tmax, g, tau) - Iexact(tmin, g, tau) \n",
"\n",
"for i, n1 in enumerate(n):\n",
" err[i,0] = np.abs(LinTrapz1(n1, tmin, tmax, G) - Iex)\n",
" err[i,1] = np.abs(LinTrapz2(n1, tmin, tmax, G) - Iex)\n",
" err[i,2] = np.abs(LogTrapz1(n1, tmin, tmax, G) - Iex)\n",
" err[i,3] = np.abs(LogTrapz2(n1, tmin, tmax, G) - Iex)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x7fab03b5bd68>"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAEYCAYAAABoYED3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXhc1Z3n//cp7arSalmbbVmWJcs7Jghsa/ECGALDj04g5NckZuskkAkhYZJJ4knHDTRhMsx0E35P00ADT6A7BGejswAdJx1iybZsecED3m3JsixbuyzJUmlX1fn9cUulKqlKlspSLdL39Tz3iXXrVNUph/JH59xzv0dprRFCCCGmmynQHRBCCDE7SOAIIYTwCwkcIYQQfiGBI4QQwi8kcIQQQviFBI4QQgi/CA90B4JVSkqKzs7ODnQ3hBAipHz00UetWuu5nh6TwPEiOzubw4cPB7obQggRUpRSF7w9JlNqQggh/EICRwghhF9I4AghhPALCRwhhBB+IYEjhBDCLyRwQlBTZx82u1T5FkKEFlkWHYLueXkfrdZ+8tIs5KfFszQ9jvz0OJamxzE3LgqlVKC7KMS06uzspLm5mcHBwUB3ZdaIiIggNTWV+Ph4n19DAifEdPUNUtfRC8Dxuk6O13W6PZ4UG8GStDhHCMWT7wgjS5T8Xy1mhs7OTpqampg3bx4xMTHyC5YfaK3p7e2lrq4OwOfQkX+FQkzjlT5S46Jo7ur3+Hh7zyAHzrdx4Hyb2/n5STHOkVB+ejz5aXHkzDUTESazqiK0NDc3M2/ePGJjYwPdlVlDKUVsbCzz5s2jvr5eAme2yEuL4+Df3kp79wBnmro409jF6cYuzjR2crbJirV/yOPzLrX3cqm9lz+fanaeiwhTLJ5rcY6ChkdFmQnR8lujCFqDg4PExMQEuhuzUkxMzDVNY0rghKgkcyTrcuawLmeO85zWmkvtvZxp7OJM00gQVbd0M+RhkcGgTXPaEViu4qLDyU+LY8lwCKXFsTQ9noTYiGn/XEJMhPxCFBjX+vcugTODKKVYkBzLguRYbl2e5jw/MGSnutXqMhoyjuFrQaN19Q1x+EI7hy+0u51Pj492GQkZR26qhajwsGn9XEKImUECZxaIDDexND2epenx/JXL+c6+Qc46Quhs00gYXen1PGRu7OyjsbOPsrMtznNhJkX2nFiWuixQWJoex4KkWEwm+S1UCDFixgeOUioJ+CmwBOgFmoCvaa2rAtqxIBAfHUFBdjIF2cnOc1prmjr7Od3Y6RwJnW7soqrFysCQfcxr2Oyacy3dnGvp5oNjDc7zsZFh5KXFsTRtZDSUnx5HiiXKL59NiFDx1ltv8cgjj1BZWUlubu6Yxx9++GFKS0upqamZ9Gu/99577Nixg8OHD1NVVcWGDRsoLS299k77aMYHDqCBF7XWfwZQSn0DeAPYFMhOBSulFOkJ0aQnRLMpP9V5fshmp+ZytzEaGp6aa+qitq0H7eEe1J4BG59c7OCTix1u51MskUb4uNw/tCQtjphImZYTwpPt27fzzW9+06fn/va3v+Xjjz9m3bp19PX1TXHPJi8oA0cpNR/4HlAAXAfEAIu01jUe2i4AfgxsARTwZ+BJrXUtgNa6w3Fu2D7gW9PZ/5koPMxEbmocualxsHrkfM/AEGebrJxp7HS7PnS5e8Dj67RaB2itukx51WXnOaVgYXLsmPuHsufEEi7LtsUst3jxYp+f+/rrr2MyGd+h4uLiqeqSz4IycIBc4PPAR8Ae4DZPjZRSscBfgH7gIYzRzA+BXUqp1Vrrbg9PexL43XR0ejaKjQxnzYJE1ixIdDvf0tXvcl3ImJ4722Sld9A25jW0hprLPdRc7uFPJ5uc5yPDTeSlWpzXhZY4pufS42XZtpg9Rk+p1dTUsGjRIl599VXq6up4/fXX6e3tpaSkhFdeeYX58+c7nzscNsEiWANnt9Y6DUAp9WW8BA7wFSAHyB++JqOUOgpUAo8BL7g2Vko95Wj/6DT1WzjMjYtiblwURbkpznN2u6a2rWdkJNRkjIpqWrvxVBpuYMjOifpOTtS7V1MYvWx7SZqxdDvJHDndH0sEqextHwS6C041/+u/+OV9fvSjH1FYWMhPfvITmpub+fa3v83WrVsDeo3maoIycLTWY69Oe3Y3UOG6AEBrfV4pVQ78FS6Bo5T6AXAncJvWumcq+ysmxmRSZKeYyU4x8+mV6c7zfYM2qpqtY+4faur0XE3B27LtuXFRRhA5puaWpMeRl2rBLGV9xAyUnZ3NO++84/y5paWF73znO9TX15OZmRnAnnkX6t/EFXieHjsB3Df8g2NkMxw2V/zUNzFB0RFhrJyXwMp5CW7nO3oGnKOhs01dzkDq6vNcTaGlq5+Wrn72VrW6nV+QHOMMouHVcjkpFiLDg2u6QYjJuPPOO91+XrVqFQC1tbUSONMkGWj3cL4NSAJQSq0AngbOAWWOuf8hrXXB6CcppR7FMd2WlZU1PT0WE5YY67maQmNnn0sIWTnT1Ellk5V+D8u2AS629XKxzb2sT7hJsSjFzBJHJYXhUdGC5FjC5P6hkOOvaaxgkpyc7PZzVJRxy0EwrEbzJtQD56q01icwVq9NpO1rwGsABQUFsuFMEFJKkZEQQ0ZCjNuybZvj+pAziJqM5dvVrd0e9w4asmsqm61UNlv5gJH7h6IjTOSlDo+GLI4giictXrZ9EOJahXrgtOMYyYzibeQjZqgwx4hl0ajrQ/1DNqpbup1TcsNhdLHNc1mfvkE7x+qucKzOfeY1Pjrcec9QvixUEMInoR44JzCu44y2HDjp576IIBQVHsayjHiWZbiXU7f2D1HZNDItNxxELV62fejsG+JQTTuHatx/j0mNi3ILIFmoIHy1c+dO0tPT3c4lJCR4aT0xFy5c4NChQwBcvnwZk8nEr3/9awBuvPFGFi5ceE2vP1mh/q34PfAPSqkcrXU1gFIqGygCtgWwXyLIWaLCuT4rieuz3AfIl639nG2yuk3LjbdQobmrn+aufvZUui9UyHLcyDo8LScLFcTVPPHEE2POrVixgoKCMZebJ2zXrl088sgjbufuu89YT/Xmm2/y8MMP+/zavlDaU12SIKCU+pzjj7cAXwW+BrQALVrrMkcbM/AJRo20H2Dc+PksEAes1lpbfX3/goICffjwYd8/gJgxtNY0XOlzC6CzTV3jLlTwZHihQn76yGgoP00WKkzWqVOnWLZsWaC7MWtd7e9fKfWRp0VZENwjnF+N+vllx/+W4aiDprXuVkrdjFHa5qcYiwM+xCht43PYCOFKKUVmYgyZiTFs9rhQodNtWu78BBYqvO9hocJwEOWlWchNtZCZECMVt8WMErSBo7We6MqyWuDeae6OEGO4L1QYOe+6UGG42OmZpi4utU9uoUJsZBiL5xrhk5tqYfFcC3lpFhYmS405EZqCNnCECFWuCxVc9x9yXahw2uU+olar54UKPQM2j0EUEabInmN2BpFrIEVHSNVtEbwkcITwk4ksVDjd2EVVcxeVzVY6ejxvhDdoG5mac6UUzE+KIXeuhby0OHLnWljsCKOEGNkeXASeBI4QATbHEsV6SxTrF7tXVLjcPUBVs9V5nGuxUtlkpbHT853kWo9UVdh1psXtsblxUeS5jogcU3Vz4+SGVuE/EjhCBCGlFCmWKFIsUW6lfQC6+gY519LtEkZdVDVbqW3r8Vh1G0bqzO07d9ntfHx0uNu0XF5qHLmpFuYlyoIFMfUkcIQIMXHRER73IOobtFFz2QiiyiYrVS1WzjVbqW7t9rg9OBg3tB6p7eBIrfvOrNERJnJShkNoJJAWzjHLvUTCZxI4QswQ0RFhLE2PZ2m6e1UFm11zsa3HCKLhUZEjjKz9nm9o7Ru0c7Khk5MN7nsRhZsUWXNiR03PxbE41UxspPxzIsYn/4UIMcOFuexDdOvyNOd5rTVNnf1UOqbkXK8VtVo9bxE+ZNdUt3RT3dLNH080uT02LzFmzMq53LkWqTcnnCRwhJillFKkJ0STnhBNSd5ct8fauweoarG6BVFVs5W6Ds/3EgHUdfRS19FL2Vn3BQsplki3+4mGrxNJBe7ZRwJHCDFGkjmSG83J3JjtvudKd/8Q1S3dVLV0uV0runC5x2N1BYBW6wCt1jYOnG9zOx8XFU6OYxT01N3LiY+enUu333rrLR555BEqKyvJzc0d8/jDDz9MaWkpNTU1k3rdzs5OXnzxRXbu3MmZM2ew2WwsX76c7373u3zmM5+Zot5PjgSOEGLCzFHhrJqfwKr57lWMB4bsXBhesDBqes5bvbmu/iE+udjBqfpOnr93lT+6H5K2b9/ON7/5zUk/r7a2lpdffplHHnmE7du3YzKZ2LFjB5/97Gd56aWXePzxx6eht+OTwBFCXLPIcBN5aXHkpcVxh8t5u11T19E75jpRVbOVTkcF7kUpZinVM47Fixf79LxFixZRXV1NbGys89ztt9/OxYsXef755yVwhBAzi8mkWJAcy4LkWG5e6r5goaWrn6pmKwO2iVfcno1GT6nV1NSwaNEiXn31Verq6nj99dfp7e2lpKSEV155hfnz5wNgNps9vl5BQQFlZWX+6r4bCRwhhN8ppUiNjyY1PnpqXvDpa9uobEo9feXqbabAj370IwoLC/nJT35Cc3Mz3/72t9m6dSulpaXjPm/37t0sXbrUL30cTQInFJ39I0QnwvwCMEmxRiFmo+zsbN555x3nzy0tLXznO9+hvr6ezMxMj8957bXXqKio4O233/ZXN91I4IQarWHnNmirhtg5sOTTkH8H5GyGKEugeyeE8JM777zT7edVq4yFF7W1tR4Dp7S0lG984xs8+OCDfPGLX/RLH0eTwAk1rZVG2AD0XIaPf2YcYVGQs3EkgOI9/4YjxIzkp2msYJKc7L5kPSoqCoC+vrHFXQ8dOsTdd9/NzTffzBtvvOGX/nkiS0NCjSkMPvUgmFPdz9v6ofJP8MG34IVl8C8bofR5aDhqjIqEELPSsWPHuP3221mzZg3vvvsuERGBu99JRjihZs5iuPufwG6H+iNw5j/gzE5oPuHeruFj4yj9nxA/H/IdI5/sEgiPCkzfhRB+VVlZyZYtW8jJyeH9998nJiYmoP2RwAlVJpOxaGB+Adzyd9BeYwTPmf+AC+VgdynK2HkJDr1hHJEWyL0FltwBebeBeY7XtxBC+M/OnTtJT093O5eQ4Pvqu+bmZrZs2cLAwADPPPMMJ0+edHv8+uuvd07D+YsEzkyRlA3rvmocfVeg6s9w5g/GNFufy/z2gBVO/s44lAkWrDNGPvl3QEpewLovxGz3xBNPjDm3YsUKCgoKfHq9kydPcuHCBQDuuuuuMY+fP3+e7Oxsn17bV0rL/L5HBQUF+vDhw4HuxrWzDULtfsfo5wNjJOTNnFxH+NwJ82+CMPl9RASfU6dOsWzZskB3Y9a62t+/UuojrbXHlJR/UWa6sAhYtME4bn8OWs4Y025nd8LFg4DLLxyXq2DfPxlHTLIx5ZZ/hzEFFxUXsI8ghJgZJHBmE6UgdalxlHwLrC1Q+Udj6u3cX2CwZ6Rtbxsc/blxhEUaiw2Gp94S5gfuMwghQpYETgiyazsmNQUr2i1z4fqtxjHYB+d3j4x+uhpG2tkG4NyHxvEf/x3SVxnTbvl3QMYaI8iEEOIqJHBC0Nc+/Bq9g70UzSuiMLOQZcnLCLvWEjcR0bDkNuOw240l1Wcdq94aj7m3bTxmHGXPQ1yG42bTO41pu4gpqo0lhJhxZNGAF8G6aKB3qJeiHUUM2ged5xKiEliXsY7CzEIKMwtJN6eP8wo+6LjoCJ8/GKMgl/d2E2GGxZuNkU/e7cYISogpJosGAksWDcwiJy+fdAsbgCv9V/hjzR/5Y80fAchJyKEws5D1mespSCsgNiLW00tNXOICuOkrxtHXaVzvOfMH4/pPb/tIu8FuOP2+caBgwU0jq95SlsjUmxCznIxwvAjWEQ5Aa28rFQ0V7K/fz776fbT2tnptG2GK4PrU61mfuZ6izCLyk/On5voPgG0ILh10VDv4g7HKzZukRSPXfbLWy5Jr4TMZ4QTWtYxwJHC8CObAcaW15mz7WWf4fNT0EQP2Aa/tk6OTndNv6zPXkxqb6rXtpLVWGsFz5g9wsQK0l421ohMdS64/Dbm3QnQQ7WUigp4ETmBJ4EyDUAmc0fqG+jjSdIR99fsory+nqmOcUQeQm5hLUaax+OBTaZ8iOnyKLvp3XzaqHJz9A1R9aFQ48MQUDtnFxuhnyachaeHUvL+YsSRwAksCZxqEauCM1tzT7Bz9VDRU0NbX5rVtpCmSG9JucI5+liQtQU3FdZehfqjZ4xj97DRqu3mTumLkuk/m9UbNOCFcSOAElgTONJgpgePKru2caTtDeX05++v3c6T5CEOuRT5HSYlJcYbPuox1pMSkXHsntDaWVJ/5g3Htp+Fj720tabDkdseS640QeY2LH8SMMNMC56233uKRRx6hsrKS3NzcaX+/9957jx07dnD48GGqqqrYsGHDVbeldiWr1MSEmJSJZXOWsWzOMr686sv0DPZwuOmwcwRUfaXarX1rbyu/P/d7fn/u9wAsTV7qXHxwfer1RIZFTr4TSkHGauPY9D3orB9Zcl1dZuzrM8zaBEf+zTjCo42pt9wtkLfF2KZBCDFpv/3tb/n4449Zt26dx83appMEziwWGxHLhvkb2DB/AwAN1gb2N4xMv13pd99F8XTbaU63nebN428SHRZNQXqB896fnIQc36bf4jOh4G+Mo98K1aVG+JzdCT0uq++G+owK2FV/hp3fM1a95W0xAii7WEY/QkzQ66+/jskxVV1cXOzX95bAEU4ZlgzuybuHe/LuwWa3cartlLH4oK6coy1HGdIj0299tj721u1lb91eANJi05zhszZjLUnRSZPvQJQFlt1lHHYbXDo8Umqn5bR72/bzcPA145DRj5jBBgcHeeaZZ3j77bepr68nMzOTrVu38tRTT7nt3lldXc3Xv/51SktLsVgsPPDAA+Tn5/PYY4+5bUVgCuB1UQkc4VGYKYyVKStZmbKSR1c/inXAyqHGQ+yr38f+hv1c6Lzg1r6pp4nfVP2G31T9BoVi+Zzlzus/a+auISJsktvamsIga61xbHkGOmqN0U3lfxpTb4PdI21l9DPrrfrXVYHugtOxh45dvdEkPPTQQ/zyl7/k+9//PsXFxezbt4/nnnuO6upq3nnnHQAGBgbYsmUL/f39vPLKK8ydO5c33niDX//611Pal2slgSMmxBJpYXPWZjZnbQbgUtclI3zq93Og4QBdg13OthrNicsnOHH5BK8fe53Y8FhuSr+J9ZnrKcwsZGH8wslPvyVmjUy9DfUbe/xU/qcRMjL6ETPU8ePH2bFjB0899RRPP/00ALfddhvh4eFs376dbdu2sXr1at566y2qq6s5cOAAN910EwB33HEHa9asoba2NoCfwJ0EjvDJ/Lj5fD7/83w+//MM2Yc43nrcufjgaOtR7C43ffYM9VB6qZTSS6UAZJozjcUH84q4Kf0mEqImeeNneBTkbDKO25+T0Y+YsXbv3g3A1q1b3c5v3bqV7du3U1ZWxurVq6moqCArK8sZNgBKKe69916OHj3q1z6PRwJHXLNwUzhrUtewJnUN/3XNf6VzoJODDQfZV7+PffX7qLPWubWv767n3cp3ebfyXUzKxMqUlc7rP6tSVhFumuR/ljL6mfWmehorWLS1GffNZWRkuJ1PT093e7yhoYHU1LFVQ9LS0qa5h5MjgSOmXHxkPLcuvJVbF96K1pqLXRed4XOw8SDdLiMQu7ZztOUoR1uO8uonr2KJsLA2Y63z+s+CuAWTe3NPo5/h8JHRjwgxycnJADQ2NrJ48cgvRI2NjW6PZ2RkcPLkyTHPb2pq8kMvJ27GB45SajHwr0Aq0A18RWs9s+7oDGJKKbLis8iKz+Kvl/41g/ZBjrYcdV7/Od56HO2yzbV10MqHtR/yYe2HACyIW+C2+s0cYZ5cBxKz4MYvGYeMfkSI2bDBuGXh5z//OX/7t3/rPP+zn/0MgE2bNgGwbt063nzzTQ4ePOicVtNa8+677/q3w1cx4wMHeBX4V63160qpLcDPlFJLtZRYCIgIUwQ3pN3ADWk38MT1T9DR18GBxgPOEVBjd6Nb+4tdF/nFmV/wizO/IFwZU3dF84p8q3w9FaOfvNuMIIqImYK/DSFG7Ny50zlVNiwhIYH777+fp59+mqGhIQoLC9m/fz/PPvss999/P6tWGavzHn74YZ5//nnuuecennvuOecqtfZ2Y/sQ16XQFy5c4NChQwBcvnwZk8nkXM124403snDh9NUzDLrSNkqp+cD3gALgOiAGWKS1rvHQdgHwY2ALoIA/A09qrWsdj88FqoFkrfWg49xZ4AtXG+XMxNI2wU5rzfnO887FB4caD9E71Ou1fXJ0MkWZRRTNK2J95nqSo5N9f3PX0U/lf0LrGe9tZfQTUDO1tI0nK1as4MiRI/z93/89P/3pT5334TzwwANj7sM5d+4cTzzxBLt27cJisfCFL3yBzMxMtm3bRkdHBwkJCVd9vzfffJOHH3543P7OqFpqSqlNwC+Aj4Aw4DY8BI5SKhb4BOgHfgBo4IdALLBaa92tlLoB2KG1XuLyvD8Br2qt/328fkjgBN6AbYBPWj6hvK6cffX7ONV2ymtb13t/iucVs3ru6skvPnA13uhnNBn9+NVMC5zpdNddd3Hq1CnOnTs3Za8502qp7dZapwEopb6METiefAXIAfK11lWO9keBSuAx4AU/9FVMo8iwSG5Mv5Eb02/kyRuepLW31Vn5YH/9ftr7R3YbHX3vjyXCYuz7M6+QoswiMi2Zk3tzb9d+PI1+5NqPCAIvvPACFouFvLw8urq6+NWvfsUHH3zAK6+8EuiuOQXdCMeVI3Bex/MI50MgWmtdNOp8GYDWeqNMqc1cdm3n1OVT7K3by776fXzS8gk2bfPaflHCIooyiyieV8wNaTdc274/MvoJKBnhePbP//zPvPTSS9TW1mKz2cjPz+cb3/gGX/rSl6b0fWbUlJqrqwROI/A7rfVjo86/DNyntZ7r+PlD4OcuiwZeBpZ4WjSglHoUeBQgKyvrhgsXLoxuIoLU8L0/e+v2Ul5fPmbxgauosCgK0gqc02+LEhb5vu+PXPvxOwmcwJqtgTMAvKC13jbq/A+BbVrrcMfPeRjLolOAHuBRrfXBq723jHBCl9aa81fOO0c/hxoPjbvtdro53bn4YG3GWuIj431/8/YLI6vbZPQzLSRwAmumXcOZUlrrSqAw0P0Q/qOUIicxh5zEHB5c8SC9Q7181PQR5XXllNeXc/7Kebf2jd2NzsoHYSqM1XNXO6ffls1ZNrml10kL3a/9XNg3UnZHrv2IWS6URzhNwG+vNqXmKxnhzFz11nrK68vZV2fs+2MdtHptmxSV5Kz7VphZeG27nk5q9JMNubcaR3aJsXWDAGSEE2izdUrtL0Ck1rp41PlSjM+18VreWwJndhi0D3Ks5Zhz+u3E5RPjtl+avNQ5/ebTtgvDrjb6cWWKgKx1IwGUtsLYOXWWksAJrNkaOE8C/4CxAKDacS4bY1n0Nq31P17Le0vgzE6Xey8bu57W7aO8vpy2vjavbWPDY1mbsdYZQPPj5vv+xpMZ/VjSHeFzM+RshthruOE1BEngBNaMCxyl1Occf7wF+CrwNaAFaNFalznamDFu/Oxl5MbPZ4E4jBs/vc+TTIAEjrBrO2fazlBeX055XTkfN3/stuvpaAvjFzrDpyCtgNgIH4t/Dg3AxQpHAH0ITce9t1UmmHfDyOgn83pj87oZTAInsGZi4HjrVJnWepNLuyzcS9t8iFHapuZa+yCBI0azDlg52HjQufhg9LYLroZrxg0HUG5iru9Lrzsb4NxfjAA69xfo6/DeNibJGPXk3gq5t0Bcuve2IUoCJ7BmXOAEAwkcMR6tNRc6LzhHP4caD9Fn6/PaPjU2laLMIgrnFbI+Y/3kN50bZrdB3RE496ERQJcOA+N8h9NWGVNvubfCgnUQHunb+waRmRY4w7XNKisryc3Nndb36uzs5MUXX2Tnzp2cOXMGm83G8uXL+e53v8tnPvOZCb2GLIsWws+UUmQnZJOdkM0Xl32Rfls/R5qOOEc/VR1Vbu2be5r5TdVv+E3Vb5ybzhVnFlM4r5CVc1YSNtFpMFMYLLjRODZtg542qN5lTL1VfQjWUTe8Nh0zjvL/DyItsGiDMfJZfAskL5qivw0RKmpra3n55Zd55JFH2L59OyaTiR07dvDZz36Wl156iccff3xa319GOF7ICEdci8buxpG6bw376Rro8to2ISqB9RnrKcwspGheEamxY3dunBCtoenEyOKD2gqwD3pvn7zYZel16Gw4JyMc33V3d6OUIjbW/f/rW265hcrKSmpra6/6GjLCESLIpJvTuSfvHu7Ju4ch+xDHW4877/051nrMbdO5K/1X2Fmzk501OwHIS8rj1qxb+dqar03uTZWC9JXGUfwk9FuhZs/I0uuOUaWa2s7BwXNw8F8gLAoWrh8JoLlLZ/XS62AyODjIM888w9tvv+3cnmDr1q1jtieorq7m61//OqWlpVgsFh544AHy8/N57LHHOH/+PNnZ2ZjNnjcwLCgooKysbNo/iwSOENMs3GRsHLcmdQ2Pr3mcjr4OKhoqnHXfWntb3dpXtleSYc7w8mqTEGWB/DuMQ2toq3ZMvf0Zzu8G172GbP1QXWocf/oBxM8bmXrL2QQxidfen2l0amnwjHiWnfa+jYYvHnroIX75y1/y/e9/n+LiYvbt28dzzz1HdXU177zzDgADAwNs2bKF/v5+XnnlFecGbMMbq13N7t27Wbp06ZT22xMJHCH8LDE6kU8v+jSfXvRptNacbT/rHP181PwRQ/YhijKLrv5Ck6GUUS5nzmJY+ygM9hlFR4eXXreM+keysw6O/JtxqDCYf+PIyreMNWCaRLkf4bPjx4+zY8cOnnrqKZ5++mkAbrvtNsLDw9m+fTvbtm1j9erVvPXWW1RXV3PgwAHnFtN33HEHa9asueo02WuvvUZFRQVvv/32dH8cCRwhAl2NnbgAAByfSURBVEkpRX5yPvnJ+fzNyr+hZ7CHg40HWT5n+fS+cUQ0LN5sHLc/B1fqRla+nSuF/isjbbXNuC/oYgXs+iHEzoHFjpVvi28Gi4/XnMRV7d69G4CtW7e6nd+6dSvbt2+nrKyM1atXU1FRQVZWljNswPhv69577+Xo0aNeX7+0tJRvfOMbPPjgg3zxi1+cng/hQgJHiCASGxHLpgWb/P/GCfPgUw8ah20I6j4aWXxQ/39xW3rdcxmO/co4ADKuM6becm+FBTeBr+V+rsFUT2MFi7Y2o9JFRob7FGt6errb4w0NDaSmjg3+tLQ0r6996NAh7r77bm6++WbeeOONqeryuCRwhBDuwsIha61x3Py30N0K53aNjIC6W9zbN3xiHHtfgMg4yNk4Mv2WmBWYzzBDJCcbZYsaGxtZvHikgnhjY6Pb4xkZGZw8eXLM85uamjy+7rFjx7j99ttZs2YN7777rtvig+kkE7FCiPGZU2D1ffDZV+HbZ+Gx3XDL38HCIjCN+p11oAtOvw/vPwkvroKXboSd/8MIqsFez68vvNqwYQMAP//5z93O/+xnPwNg06ZNAKxbt47a2loOHhzZ6ktrzbvvvjvmNSsrK9myZQs5OTm8//77xMT4bx8mGeEIISbOZDKm0DKug5JvQ1+nseJtePHBlVEXqFvPGkfFy8aePwuLRpZep+TJ0msXO3fudE6VDUtISOD+++/n6aefZmhoiMLCQvbv38+zzz7L/fffz6pVqwB4+OGHef7557nnnnt47rnnnKvU2tvbATA5Fnk0NzezZcsWBgYGeOaZZ8aMiq6//nqioqKm7TNK4AghfBcdD8vuMg6tobVyZOqtZi8MuZT7GeozHjv3Ifzxf0BCljHtNrz8OkRuPJ0uTzzxxJhzK1as4MiRI+Tk5PCTn/yEH/7wh2RmZvK9732Pp556ytkuMjKSP/3pTzzxxBN89atfxWKx8IUvfIG1a9eybds2EhKMUkonT57kwgXjfqy77rprzPsN368zXaTSgBdSaUCIazTYCxfKR8rujLfnz5PHIXHBhF52plUamE533XUXp06d4ty5c1P2mlJpQAgRfCJiRqbPADpqR248rS4zrveAUdVggmEjvHvhhRewWCzk5eXR1dXFr371Kz744ANeeeWVQHfNSQJHCOEfiVlQ8Ihx2Abh0iEjfMxyH89UiIqK4sc//jG1tbXYbDby8/N54403+NKXvhTorjlJ4Agh/C8sAhYWGoeYEo8//vi0V3u+VrIsWgghhF9I4AghQo4sdgqMa/17l8ARQoSUiIgIenvlJtJA6O3tvaaqBBI4QoiQkpqaSl1dHT09PTLS8ROtNT09PdTV1Xms2TZRsmhACBFS4uPjAaivr2dwcJwdTcWUioiIIC0tzfn37wsJHCFEyImPj7+mf/hEYEx6Sk0p9V/GeWx216YQQgjhlS/XcP5m+A9KqdF3FH1RKbXy2rokhBBiJvIlcMJc/jx6tPNvwF/73h0hhBAzlS+Bc0EptdHxZ7fa4lrrfty2BhRCCCEMvgTOy8CbSqkbGRUuSqkVgFzHEUIIMcakV6lprc8opV4AKoA2pdT/AlqBXOBzwMbxni+EEGJ28mlZtNb6JaXUBeDvge86Tp8G/lprfWKqOieEEGLm8Pk+HK31e8B7Sql4IFxr3TZ13RJCCDHTXPONn1rrzqnoiBBCiJlNaqkJIYTwCwkcIYQQfiGBI4QQwi8kcIQQQviFBI4QQgi/kMARQgjhFxI4Qggh/EICRwghhF9I4AghhPCLGR84SqkkpdT7SqmzSqlPlFJ/UkrlBrpfQggx28z4wMHYQuFFrfUSrfV1wPvAGwHukxBCzDpBEThKqflKqX9SSu1XSvUopbRSKttL2wVKqV8rpa4opTqVUv+ulMry9tpa6w6t9Z9dTu0DPL62EEKI6RMUgYOxl87ngXZgj7dGSqlY4C/AUuAh4AEgD9illDJP8L2eBH53Tb0VQggxaddcLXqK7NZapwEopb4M3Oal3VeAHCBfa13laH8UqAQeA14Y702UUk85nv/oFPVbCCHEBAXFCEdrbZ9g07uBiuGwcTz3PFAO/BWAUupBpdTHjuPx4XZKqR8AdwJ3aK17pq73QgghJiJYRjgTtQLP02EngPsAtNb/Bvyb64OOkc2dwG1a6yvT3UkhhBBjhVrgJGNc5xmtDUjy9ASl1ArgaeAcUKaUAhjSWhd4aPsojum2rCyv6xCEEEL4INQCZ9K01icANcG2rwGvARQUFOjp7JcQQsw2QXENZxLa8TyS8TbyEUIIESRCLXBOYFzHGW05cNLPfRFCCDEJoRY4vwfWKaVyhk84bhAtcjwmhBAiSAXNNRyl1Occf7zB8b93KKVagBatdZnj3OvA14HfOZY5a+BZ4CLwL/7srxBCiMkJmsABfjXq55cd/1sGbALQWncrpW4Gfgz8FGMxwIfAk1prq5/6KYQQwgdBEzha64muJKsF7p3m7gghhJhioXYNRwghRIiSwAlRAwMDtLW1YbPZAt0VIYSYkKCZUhOT09raytmzZzGZTCQkJJCcnExSUhJmsxlHNQUhhAgqEjghqq2tDQC73U57ezvt7cZ9r5GRkSQnJzsDKCIiIpDdFEIIJwmcEGU2m+np6aGnx73w9cDAAI2NjTQ2NgIQFxdHdnY2c+bMCUQ3hRDCSQInRC1atIhFixbR19dHe3s7bW1ttLe3MzQ05Nauq6vL4/MHBgaIjIz0R1eFEAKQwAl50dHRZGRkkJGRgdaarq4u2traaGtro7OzE6UUiYmJbs+x2WxUVFQQFRXlnHpLSkoiLCwsQJ9CCDEbSODMIEop4uPjiY+PJzs7m6GhIaxW65gg6ejowG6309vbS11dHXV1dSil3BYfWCwWWXwghJhSEjgzWHh4+JjRDUB/fz8mkwm7fWSjVa01HR0ddHR0ABAREUFycjIpKSnMnTvXb30WQsxcEjizUGZmJunp6Vy5csU5/dbd3e3WZnBwkKamJmw2mwSOEGJKSODMUiaTyXntZvHixfT397stPhgcHAQgOTl5zHMrKyvp6+sjKSmJ5ORkYmJiZPpNCHFVEjgCgKioKNLT00lPT0drjdVqpa2tbUzgaK1paWlhYGCAy5cvA8bCheHwSUpKIjxc/rMSQowl/zKIMZRSxMXFERcXN+axnp4eBgYG3M719fXR0NBAQ0ODc+HCcADFxcXJ6EcIAUgtNTFJZrOZtWvXkpeXR0pKypgVcFprrly5Qk1NDUeOHHFOzQkhhIxwQpDWOqCjhpiYGObNm8e8efOw2+10dnY6r/243mgaFxc35ubSrq4umpqaSE5OJiEhQe79EWIWkcAJQXX/7VvYe3uI27wZy8aNRGRkBKwvJpOJxMRE5/LrgYEBZ203s9k8pn1rayuXLl3i0qVLboVHk5OTiY2Nlek3IWYwCZwQY+/txbprF7q/n+6y3QBELVtG3OZNWDZtInrlSpQpcDOlkZGRpKWlkZaW5vHx4SKj4F549Ny5c0RFRbktPpDCo0LMLBI4Iab344/R/f1u5/pPnaL/1ClaX36FsJQULBs3ELd5M+bCQkyxsQHqqWfZ2dnOe39GFx7t7+93Kzy6dOlS0tPTA9FNIcQ0UFrrQPchKBUUFOjDhw8HuhseDVy6hHVXKdbSUroPHgQvF+ZVZCSxa9di2byJuE2biMjM9HNPx9fX1+e89uOp8GhBQQEWi8XtXFNTEwkJCURHR/uzq0KICVJKfaS1LvD4mASOZ8EcOK5s1m66y8uxlpZiLSvD5tgnx5Oo/Hxn+ESvXh3QqbfRtNZ0dnY6bz7t7+9n3bp1btd0+vr6qKioAIyFC8PXfhITE2XxgRBBQgLHB6ESOK60zUbv0aNYS8uw7tpF/9mzXtuGzZmDZeNGLJs2Yi4sIswy9gJ/INntdkyjArGhoYEzZ86MaetaeDQ5OVl2PRUigCRwfBCKgTPaYF0dXaWlWHeV0nPgANrb1FtEBLE33YRl82biNm8iYt48P/d0Yi5fvkxdXZ2z2rU3kZGRpKenk5OT48feCSFAAscnMyFwXNm7u7Hu22dc+ykrw+YoS+NJVF4els2bsWzeRMzq1aggm66y2+3jFh4FmDdvHnl5eW7n+vv7iYiIGDNyEkJMHQkcH8y0wHGl7Xb6jh2ja9curKVl9J8+7bVtWHIylg0bsGzejLmokLBRF/GDgafCoytXriQlJcWt3SeffEJnZyeJiYnOpddSeFSIqSWB44OZHDijDdbXG1NvpaX0VBxAj6qV5hQRgfnGG7Fs2oTl5s1Ezp/v345OwHDh0djYWLeFBDabjfLy8jFTcdHR0W67nkrhUSGujQSOD2ZT4Liyd3fTXVHhHP3YWlu9to3MXWxUO9i8mZjrrgu6qTdX3d3dHDt2jL6+Pq9tXAuPzp8/X8JHCB9I4PhgtgaOK22303fiBNZdu+gqLaX/5CmvbcMSE7FsHJ56KyLMQ6XpQNNaO+/9aWtro6OjA5vNNqZdWFgYRUVFbtd6Al2/TohQIYHjAwmcsQYbG437fXaV0l1RMabigVN4OLE3Fhijn02biMzK8m9HJ8i18GhbWxtWqxWAOXPmsGrVKre2LS0t1NTUOKffpPCoEJ5J4PhAAmd89p4euisqnBUPhlpavLaNXLzYWestZs0aVJBOVQ0XHo2MjCQpKcntsTNnztDQ0OD8ebho6XDtNyk8KoRBAscHEjgTZ0y9nXSMfnbRd/Kk17ZhCQmYN2wgbvMmzCUlQTn15smBAwfo7e31+rgUHhXCIIHjAwkc3w02NTmrHXTv3z/+1NsNNxjldjZvJnLhQv92dBJsNhsdHR3OpdejC4+6WrZsmddq2ULMdBI4PpDAmRr23l73qbfmZq9tIxctMm443bSR2E99Kmin3oAxiw9cC48WFha6bTynteb06dPO8jtSeFTMZBI4PpDAmXpaa/pOnnQuPOg7ftxrW1N8PJaSEiOASooJS0jwY08nZ7jw6HDR0aVLl7o93tnZyZEjR5w/S+FRMZNJ4PhAAmf6DTY1Y91dZqx627cP7e0embAwYq+/HvPGDVg2bCBqyZKQukBfU1NDTU2Nx8ek8KiYaSRwfCCB41/2vj56Dhxw3nA65NiEzZPwtDTMJcVYNmzAXBic5XZc9fb2cvnyZef023iFR1NSUli5cqUfeyfE1JLA8YEETuBorek/fdoIn12l9B075r1xeLgx+tlQEhKjH7vdTkdHh7P22+jCo4sWLWLhqMUTly9fJiwsjPj4eCk8KoKeBI4PJHCCx1BLC9a95XTv2Y21fB/2K1e8tg1PS8OyoQRzSUlIjH76+/vddj1dvXo1caOWih86dIju7m7CwsKchUeTk5OJiYkJUK+F8E4CxwcSOMFJDw0Zm8zt3k337j3j3vMzPPqxbNyAuWQDUUvygnr0M/xddO1jf38/+/fv99h+uPDo8OIDqf0mgoEEjg8kcELD8OjHuruM7vJ92Ds7vbYNT0/HUlIcMqMfMJZf19bW0tbWNqHCo6tWrZLgEQElgeMDCZzQ4zr6se7ePW6xUcLDif3UpxzTb6Ex+unt7XVe+/FUeDQ6Opp169a5nRscHMRutxMVFeXP7opZTAIHUEo9AvwE+KzW+rdXay+BE/qGWlqw7tmLdc/uCY5+SjBvKMG8fn3Qj348FR7NzMxkyZIlbu0uXbpEVVUVZrPZWXYnMTFRFh+IaTPrA0cplQ28AyjgeQmc2UcPDdH7ySdYd+/BumcSo58NG4jKC+7RDxiFR7XWY0YyR48epa2tze3ccOHR4QCSwqNiKgV14Cil5gPfAwqA64AYYJHWusZD2wXAj4EtGOHxZ+BJrXXtOK9vAv7keI9/BF6UwBGDzc1079mLdc8eusvLsXd1eW3rPvopJMxi9mNPfae15vjx47S1tTHe9zwqKork5GTmzZuHJchHdiL4BXvgbAJ+AXwEhAG34SFwlFKxwCdAP/ADQAM/BGKB1Vpr9xsaRp7334E4rfVTSqlSJHDEKG6jn9276T81zugnIsLl2k9JSIx+XAuPtrW1ea16fd11143ZlkE2nhOTFeyBY9Ja2x1//jLwOp4D55vAC0C+1rrKcW4RUAl8V2v9gofXXul4vQ1a60EJHDERkxr9ZGQYNd82lBC7bn1IjH5cC4+2t7djs9k87nI6NDTEgQMH3Pb9kcKj4mqCOnBcXSVwPgSitdZFo86XAWitNyqlHgS+5XjodcAO/B3GqAggHegE/l5r/dJ4fZHAEeAY/Xz8sePaz54Jjn42YNlQQmRubtCPDoYLj/b19Y3ZUqGlpYUTJ064nYuNjXWGjxQeFZ7MlMBpBH6ntX5s1PmXgfu01nMn8PqljDPCUUo9CjwKkJWVdcOFCxd8+RhiBnOOfnbvpnvfvhk3+nFVVVXFpUuXvD6ulHKOfubMmYPZHFqfT0yPmRI4A8ALWutto87/ENimtb7q3W4ypSamkh4cdL/2c/q098YREcZmc44ACpXRT09Pj3P67cqVK14Lj6amprJ8+XI/91AEo/ECZ1bdkqy13hToPoiZQ0VEEFtQQGxBAanf+m8MNjXTvXcP1t17xo5+Bgfpqaigp6KC5v/zfwjPzMBSYky9mdetwxSEowOlFGazGbPZzIIFC7DZbFy5csV57ce18GhycvKY558/fx6tNcnJyVJ4VAChFTjtQJKH88mOx4QIqIi0VBLvvZfEe+81Rj+u135GjX6G6hvo+MUv6PjFL9xHPxs3ELl4cVCOfsLCwpy128C98OjowNFaU1dXx9DQELW1tVJ4VAChNaX2FyBSa1086nwpxufYOJV9kSk1MZUGm5ro3uMy+rFavbYNz8zAUlyCuaQ4JKoeeDJ6l9PRhguPpqSkeBwdidA1U67hPAn8A7BEa13tOJeNsSx6m9b6H6eyLxI4Yrq4jX5276b/zBnvjcPDiVlznTOAopctQ4XA1JTNZnOOfsYrPJqWlsayZcv83DsxnYI+cJRSn3P88Rbgq8DXgBagRWtd5mhjxrjxs5eRGz+fBeIwbvz0/iujDyRwhL9MZvQTNmcOluIizMUlmIsKCQ+B0cF4hUeXLVs2Zjm2CG2hEDjeOlHmeqFfKZWFe2mbDzFK29RMdZ8kcEQgOFe+7dlL956r7PejFNErVhjbbZeUELN6NSoEtiaw2+1cuXKF9vZ25s+fT2RkZKC7JKZQ0AdOMJLAEcFgqLWV7n37jADauxdbu/f1Maa4OMzr1xsBVFxMREaGH3sqhEECxwcSOCLYaLudvhMnjaXXe/bS+/HH4OW+GICovFzMxSVYSoqJueEGTLInjvADCRwfSOCIYGfr7KR7f4UzgIYaG722VdHRxK69CYsjgCIWLgzKpdci9Eng+EACR4QSrTUDVVWOqbc99Bw6jB4c9No+YsECY7vt4mJib1obcmV3RPCSwPGBBI4IZfaeHnoOHXIuPhgYry7gcNHRkmJjy4UlS2T0I3wmgeMDCRwxkwxcvEj33r1GAFVUoHt6vLYNnzsXc3GxEUCFhYQlJvqxpyLUSeD4QAJHzFR6YICeI//Xee1n3BtPTSZiVq3CXFKCpbiI6FWrULIlgRiHBI4PJHDEbDHY1Ex3ebkRQOX7sF+54rVtWEIC5qJC48bT4iIiUlP92FMRCiRwfCCBI2YjbbPRd+yY876f3qNHYZx/I6Ly8x2LD0qI/dT1KLmJc9aTwPGBBI4QMNTeTs/+/Vj37MW6dw+2llavbU2xscSuW+dc/Ra5YIEfeyqChQSODyRwhHCntab/7Fmj7tuevfQcOQLjLL2OXLjQuPZTUkzsTTdhki0JZgUJHB9I4AgxPpu1m56DB7Du2UP3nr0MjrcddWQksQU3OCsfhMKOp8I3Ejg+kMARYuK01gxeuOCceus5cBDtZUsCgPD0dOe1H/P6dYTFx/uxt2I6SeD4QAJHCN/Z+/vpOXyY7r3G6rf+yirvjcPCyHn/PaIWLfJfB8W0GS9wgr+WuRAi5JiiorAUFWEpKoLvfZfBhgase/fSvWcv3fv3Y+/qcrYNS04iMjs7cJ0VfiOBI4SYdhEZGSTddx9J992HHhqi9+hR57WfqLw8uZ4zS8iUmhcypSaEf2ibTaoXzCDjTakF/+boQogZTcJm9pDAEUII4RcSOEIIIfxCAkcIIYRfSOAIIYTwCwkcIYQQfiGBI4QQwi/kPhwvlFJdwDhbIc44CYD3nbf8xx/9mMr3uNbX8uX5k3nORNtOpF0K4H1/gplHvhO+ydNaJ3h8RGsth4cDOBzoPvj5874W6D74qx9T+R7X+lq+PH8yz5lo24m0k+/EzO2Hv74TMqUmhr0X6A44+KMfU/ke1/pavjx/Ms+ZaNtg+f8/mATL38mM+U7IlJoXSqnD2kt5BiFmI/lOiGslIxzvXgt0B4QIMvKdENdERjhCCCH8QkY4Qggh/EICRwghhF9I4PhIKbVdKXVWKWVXSn0m0P0RIpCUUklKqfcd34lPlFJ/UkrlBrpfIrhI4PjuP4FPA7sD3REhgoAGXtRaL9FaXwe8D7wR4D6JIDNrAkcpNV8p9U9Kqf1KqR6llFZKZXtpu0Ap9Wul1BWlVKdS6t+VUlmubbTWFVrran/0XYjpMJXfCa11h9b6zy5P2Qd4fC0xe82awAFygc8D7cAeb42UUrHAX4ClwEPAA0AesEspZfZDP4Xwl+n8TjwJ/G5KeytCXnigO+BHu7XWaQBKqS8Dt3lp9xUgB8jXWlc52h8FKoHHgBf80Fch/GFavhNKqacc7R+dpn6LEDVrRjhaa/sEm94NVAx/sRzPPQ+UA381HX0TIhCm4zuhlPoBcCdwh9a6Z6r6KmaGWRM4k7ACOO7h/AlguZ/7IkQwmNB3wjGy+X+A27TWwVBlWQQZCZyxkjHmtEdrA5KGf1BKPa2UugSsB95QSl1SSs33Ux+F8KerfieUUiuAp4E5QJlS6mOl1GG/9VCEhNl0DWdKaa2fxviCCTHraa1PACrQ/RDBTUY4Y7XjMpJx4e23PCFmOvlOiCkhgTPWCYw569GWAyf93BchgoF8J8SUkMAZ6/fAOqVUzvAJx81wRY7HhJht5DshpsSs2p5AKfU5xx9vAb4KfA1oAVq01mWONmbgE6AX+AFGyY5ngThgtdba6u9+CzFd5Dsh/Gm2BY63D1umtd7k0i4L+DGwBeNC6IfAk1rrmunuoxD+JN8J4U+zKnCEEEIEjlzDEUII4RcSOEIIIfxCAkcIIYRfSOAIIYTwCwkcIYQQfiGBI4QQwi8kcIQQQviFBI4QQgi/kMARQgjhFxI4Qggh/EICRwghhF9I4AghhPALCRwhhBB+ER7oDgghxqeU+ivgVuA64CGMrZ3vczxcDPxPrfV/BKh7QkyYbE8gRBBTSkUC/1tr/aRS6hDQD/wGeEFrrZVS24Cvaa2zAtpRISZAptSECG4bgb1KKQXkAI1a63/UI78pKmBOwHonxCTIlJoQwe040AGsxphKe3HU49dhbP8sRNCTEY4QQUxr3aC17gU2Az3AweHHlFIRwKeB9wLUPSEmRQJHiNCwGSjXWg+4nLsDiAd2KKXClFILA9M1ISZGAkeIIKeUCgM2ALtGPfQgsEtrXYOxiu1Tfu6aEJMigSNE8LseSARKR51fAvxeKRUO/L/A+37ulxCTIosGhAh+mcAx4NCo889h3JezHPgHrfWgvzsmxGTIfThCCCH8QqbUhBBC+IUEjhBCCL+QwBFCCOEXEjhCCCH8QgJHCCGEX0jgCCGE8AsJHCGEEH4hgSOEEMIvJHCEEEL4hQSOEEIIv/j/AVYnNEoqY51uAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.plot(n, err[:,0], label='Lin1')\n",
"plt.plot(n, err[:,1], label='Lin2')\n",
"plt.plot(n, err[:,2], label='Log1')\n",
"plt.plot(n, err[:,3], label='Log2')\n",
"\n",
"plt.plot(n, 1.0/n**2, '--', c='gray', alpha=0.5)\n",
"\n",
"plt.xlabel(r'$n$')\n",
"plt.ylabel(r'$\\epsilon$')\n",
"\n",
"plt.xscale('log')\n",
"plt.yscale('log')\n",
"\n",
"\n",
"plt.xlim(10, 500)\n",
"#plt.ylim(1e-4, 50)\n",
"\n",
"plt.legend(loc='upper right')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Thus, overall `LogTrapz2` seems like a good choice to integrate integrands that decay smoothly or rapidly over large domains."
]
}
],
"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.9"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment