mirror of
https://github.com/vale981/phys512
synced 2025-03-04 09:11:37 -05:00
436 lines
104 KiB
Text
436 lines
104 KiB
Text
{
|
||
"cells": [
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "cd079071",
|
||
"metadata": {},
|
||
"source": [
|
||
"# Integration exercises"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "c0431bce",
|
||
"metadata": {},
|
||
"source": [
|
||
"## Newton-Cotes"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 10,
|
||
"id": "18ed0344",
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"3 0.34075853066724404 0.0519405510314801 0.0022798774922103693\n",
|
||
"5 0.1834653418221377 0.012884199027224597 0.00013458497419382986\n",
|
||
"9 0.0949599423108507 0.003214828113830337 8.295523967971619e-06\n",
|
||
"17 0.04828406569741284 0.0008033195149277361 5.166847063531321e-07\n",
|
||
"33 0.024342886926189244 0.00020080567998093102 3.2265001115305836e-08\n",
|
||
"65 0.012221646395186303 5.0199907898895724e-05 2.0161288194486815e-09\n",
|
||
"129 0.006123373269068644 1.2549882473678053e-05 1.2600120946615334e-10\n",
|
||
"257 0.003064824111058906 3.137464712255067e-06 7.87503395827116e-12\n",
|
||
"513 0.0015331964220766103 7.843658088590999e-07 4.922728891187944e-13\n",
|
||
"1025 0.0007667943025135848 1.960914292054028e-07 3.086420008457935e-14\n",
|
||
"2049 0.00038344617411567583 4.902285577479404e-08 1.9984014443252818e-15\n"
|
||
]
|
||
},
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGiCAYAAAAfnjf+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACGK0lEQVR4nOzdd3QUZd/G8e/upnfSewMSegIBQq/B0EITBCvV3hELFsAKShERHrsCVkSlSIdApPdekpCQQEjvyW6STdl9/4jmeXiJChqYlN/nnD3HnZ3MXhuRvZx75r5VRqPRiBBCCCFEA6dWOoAQQgghRF2QUiOEEEKIRkFKjRBCCCEaBSk1QgghhGgUpNQIIYQQolGQUiOEEEKIRkFKjRBCCCEaBSk1QgghhGgUTJQOcLsYDAbS0tKwtbVFpVIpHUcIIYQQN8BoNFJcXIynpydq9V+fi2kypSYtLQ0fHx+lYwghhBDiH0hJScHb2/sv92kypcbW1hao/qXY2dkpnEYIIYQQN6KoqAgfH5+a7/G/0qBKzYYNG3juuecwGAy8+OKLTJs27YZ/9o8hJzs7Oyk1QgghRANzI5eONJhSU1lZyfTp09m1axf29vaEhYUxevRonJyclI4mhBBCiHqgwdz9dPjwYdq2bYuXlxc2NjYMGTKEbdu2KR1LCCGEEPXEbSs1u3fvJioqCk9PT1QqFWvXrr1un2XLluHv74+FhQXh4eEcPny45rW0tDS8vLxqnnt5eZGamno7ogshhBCiAbhtpUan0xESEsKyZctqfX3VqlVMnz6d2bNnc/z4cUJCQoiMjCQrK+sfvZ9er6eoqOiahxBCCCEar9tWaoYMGcJbb73F6NGja3190aJFPPjgg0yePJk2bdrw8ccfY2VlxZdffgmAp6fnNWdmUlNT8fT0/NP3mzt3Lvb29jUPuZ1bCCGEaNzqxTU15eXlHDt2jIiIiJptarWaiIgIDhw4AEDXrl05e/YsqampaLVaNm/eTGRk5J8ec+bMmRQWFtY8UlJSbvnnEEIIIYRy6sXdTzk5OVRVVeHm5nbNdjc3N2JjYwEwMTFh4cKF9O/fH4PBwAsvvPCXdz6Zm5tjbm5+S3MLIYQQov6oF6XmRo0YMYIRI0YoHUMIIYQQ9VC9GH5ydnZGo9GQmZl5zfbMzEzc3d0VSiWEEEKIhqRelBozMzPCwsKIjo6u2WYwGIiOjqZ79+4KJhNCCCFEQ3Hbhp+0Wi0JCQk1z5OSkjh58iSOjo74+voyffp0Jk6cSOfOnenatSuLFy9Gp9MxefLk2xVRCCGEEA3YbSs1R48epX///jXPp0+fDsDEiRNZvnw548ePJzs7m1mzZpGRkUFoaChbtmy57uJhIYQQQojaqIxGo1HpELdDUVER9vb2FBYW3pIFLau0WjQ2NnV+XCGEEKIpu5nv73pxTU1DZ6yo4GLPXiQMGEhlTs5/t1dVKZhKCCGEaFqk1PxLVQYjD81bi1GvpzS/kH05VWQVlWE0Gsmav4CEAQPJX71a6ZhCCCFEo9eg5qmpjy7n6tius2LfsDfx0OWRuPwoAI7WZrwTcwC/tDSOXCnC/WoBLV1t0WRlcOWBB7Ds2BHPBfNRqVQKfwIhhBCicZBS8y+52lnw2QOdicso4kJGMaQXkZSjI09XznMh99HcP5UrKdYULN2HWgWjiuJ4MC2NXI0F585n0trDDi8HSzLfeJ2qvHycHnwQy/btlP5YQgghRIMjpeZfsjE3YVAbNwa1+e9dWmUVVVzM1BKbUURsRjG2GUVcSC8mT1fOZgs/Eno+gomhiuNfHwPA2kzDlxu3Yqcr4EjnO/C29STY3RaT82fJW74c6149aXbXXUp9RCGEEKJBkFJzC1iYamjvbU97b/uabUajkWytnriMYmLTi7mQUURZejEJWVp0+kreCR1Pi4KrbDitR3++ehHPqSl7GHtsG3FZWkqCe9HK3ZYAZ2ty3p2HibMLDuPGYtKsmVIfUwghhKhXpNTcJiqVCldbC1xtLejd0qVme0WVgeQcHRcyOhGbXkTPjGLiMopJLSgl2q4FBW2HcdXclUPfnwDAjgq+X/sNaoxs9OxIYLAfrdztsDpzjNKz57Du0QPLdm2V+phCCCGEYqTUKMxUo6almy0t3WwZEeJZs72wtKL6rE5GEWbpxVRkFBGXUUylrpSv2g7FU5vDkr2ZsLd6vawXzv5C/4T97DmZTOWUR2jpZkuLZubo/rMU85YtsY8ajspE/nULIYRovORbrp6ytzSla4AjXQMca7YZDEau5pdyIaMXcRnFDM0oIja9mKRcHUft/FB76dhW5sTxn04D4FeUzsc7v6TMzJJNpi1p6WZHCzcbPPdvR5OVie0dg7Bo1UqpjyiEEELUKSk1DYharcLXyQpfJysi2/539fLS8iriM3sSl1FMSEYRFhnFXMzSUq4zZV1gL1RGI5/tTa7Zf+7elYTmJPJVkh79QAMtXW0IUmlxXLUcuw7tcXzgfgU+nRBCCPHvSKlpBCzNNIT4OBDi43DN9oKSci5mRXExU8vkrGIuZmq5mFVMjHcnMq2c+LXckZS9SQB0TzvLrMO/cmH/Cb4tD6KFqw1Bbra03bASBypwu/9eLIKCFPh0QgghxI2RUtOIOViZ0cXfkS7+jtdsLyztS0JWMe6ZWi5maYnPLEaLNytaD6bY1JK9CTnsTahe7mHF1o2YlRYwtcgTY7s8Wrra0LEgmebbVuPQpzc+D09T4qMJIYQQ15FS0wTZW5oS5udImN//lp1wisvu4mKWlr6/n9GJzyhmbd4oHNIvc8LEiZLEXPYn5qK9+BvB5w6zKb2Mz7L9aOFqQ0tXW4Z9+y6Wdja4v/A8bsGBin0+IYQQTZOUGlHD1sKUTr7N6OT7P3PfTA1Hq69kZJaWi5nV8+pku1Sx0t6WeKM1OdpycrR5HL+Yyfgzx9BgZPBHfVA5JdLC1YYhlw/R/lg0qjuG4PPQFFxszFGpVBiNRlkiQgghRJ2SUiP+lo25CaE+DoT+cc3O0NbACErKK0nI0nIxU0tCWj6/WD1P1eUkCixsMerKOZyUR9eTp7C+nMAPv51jRXY0DlamBDtZ8NIXMyh396L0jQU0D/TA08ESVVkpKnNzVBqNkh9XCCFEAyWlRvxjVmYmdPB2oIO3A4R5Q1R7AF4qryIxu3oI62qIHevOhHPRaIMKKCipIDUzFYuSYgxXkpi0+gKoYjE3UTM9dgM9Yvdwcfh9mNx9H81drAloZok6JRkzf3/UZmaKfl4hhBD1m5QaUecszTS087KnnZc9dPSGu3sD1WtiJWZruZhawJ4unuReSSfY1o6kHB36SgOW2WmYVFawJVnLtt9nUPbQ5fDl9nmUm5qz7q1vaO5uS3MXG/xyU7C30GDRvDlqKyslP64QQoh6QkqNuG0sTDW09bSnrac9dPEDYCZQZTByNb+EhPSOHIlNxK1EQ2cdJGRrccwtpsTEnAxLR748cLnmWHMOfEF45gXW9buX3IFRNHe1oaVFFb5nDuIe0gbbLp0V+pRCCCGUIqVGKE6jVuHnZI2fkzW087zmtTxdPxIy70Z3OZNpOkjM1pKYrUNnZkmeuS37Kmw5c+wqAJ2y4nh7/2ccsHVl/t1v0NzFhuau1oSd2oWrhRrvqCE4BPgq8RGFEELcBlJqRL3maG1G10AnugY6XbO97Nk+JOXoeDirmMTsEhKztRhOZnDcvRVpls24mFU9Bw/n4NMd36PRZjP5ZAnpLTrQ3NWazvpMOh/bjmWnTnjffw9uduZyN5YQQjRwUmrqwvonybXsRrNeI1Bb2iqdpkmwMNXQ2sOO1h52/914d0cMhimkFZbSNVtHYpaWhGwtCRk9yEhNJtnOnfyiMjKKynBIOszgU9EcSkxj7FU3rM00NHe14eGdn2FjYUrFpIfwC2mNn5M1ZiZq5T6oEEKIGyal5t/KOIP+6E/8lD0Ei7VbGNv/GNY97wGPDkona5LUahXezazwbmZF3yCX6o2j5wEwsqSCxBwtiVlask+q2GuvIl5jj0atQldexdmUfLxjj2NmqGSyUx8ydmajUasYn3OKqNNbyOjSF/390wh0sSHQxRqHihI0Dg5yhkcIIeoJKTX/lp0XeZ3eRLO9ClOVDqszH8HZj8ArjKKgadiEj0BtYaN0SgHYW/3P5IKdfWDaUADerjRwJU9HQnoh571foSwhEffm/mhzS9HqK7FKu4xDQRYxiel88ssZAFRGA79seBWjRsPaJ97FrWUAgc7WBBi0eJoZsGkeILegCyHEbaYyGo1GpUPcDkVFRdjb21NYWIidnd3f/8BNqiyvpPjUXpolL4cLGzBWVfJNzn8wYsLQXrE4DxgLbm3q/H3FrWM0Gskq1pN48SqZp8+TVGHGSRNHLmVrKU3PYMWWtzCqVIyKmkuVunrCwEnnNjH+4k5iWvXm4IhpBLpYE+hsTesTMbi1aYlXjy5SdoQQ4ibczPe3nKmpIyZmJjTr0g+69ANtNoW7f0S/wRajEexjP4S4heATTkmbKVh0GoHaXOZWqe9UKhVudha4hbWAsBbXvFZWUUXSs4NIiU3iWdNmXMrWkZijwyIWdCYWXDB14rf4bH6Lz6ZZWRHfbXmXQlQMv/M9fNztCXS2oVvaaXyKs2nWtzeBPcOwMpP/HIUQ4t+QMzW3UGVZBblH9+CW+hXEbgRjFevzZlFQ5cOA8CS8h4wCl+DbkkXcHkajkexiPUkZhVwq0HMpW0tu/CV6bPyKSn05L/V8pGbfmYdX0iftNJ+2i2JNi7542FvQ1sbIhOjlqPz9MXn0KQJdbPBysEStlut2hBBNk5ypqSdMLExx6zUAGADFGegPfUvWL0HoDdbYxr8Kl94Bv56Ut5+MaYcoVGYWSkcW/5JKpcLVzgJXOwvC/9g4rA08O5zySgOd80q4lK3lUo4ONd04dcGKdI/mAKQXltHsUhJescfIvJzIRIvuAJibqJlxfh3NC1O5OmwCNv36EuhiQ4CjBXamalQynCWEEICUmtvH1h3ziOeY1KuctD2/YZ/TGeI3w+V9HDzThuTlRnp2SqP5iOHg3OLvjycaHDMTNS1cbWjh+vuF432fBWACUFBSTmK2jivxXpxyqySnuIwWrjZczq1eQsItNQGPwlQ+O3mVA1mnAAjKv8Ki3UtJ8WzB7sffJNC5+q4sf30+3t4umDs7/UkSIYRonKTU3GYmFmb4DhoEDIKiNIzHvib55wCKK50xu7gUls4G/95UdZyMuvVwVGbmSkcWt4GDlRlhfmaE+TWDQaEAPAFUVhm4ml/Klf6OxJ6Lp5VzcwylGi7l6PC+ko3GaKBAX8X3h1NqjvX+b0soz7/CpxEPo+3ai0Bna1pqygjMTsIjpC3u7YLkNnQhRKMkpUZJdp6o+r/IPd3LSdoWg3eROySchuQ9nD1nx6nSKrqE5ND6zsHgGKh0WqEAE40af2dr/CO7Q2R3Rv/Pa4WlfUi+MAG79FyeNHWsvlg5W4u5oQKA4wYbLp/PBKBfynFePPYde5wCeHPQ09W3nztb0+9MNM3srXEeGklAkK9crCyEaNDkQuH6piAFjq/kl/UupJe1oI/dJ7S32gKB/TGGTYagIahM5RoK8ecMBiOpmXkkFVSQlF9GUo4Ok7276LR/A6fsffmk/ciafX/YNBv7ch1P9HuWRAcvPOwtGKBLplfcPio7dcVm9BgCXazxcrDERCMzKwshbj+5ULghc/CBAa8Q1V1P4qadBBSrIVkFl3Zx5UIhu4vLCGlfRIe77oBmfkqnFfWQWq3Cx8MJHw/o88fGEW2BJxhQUcWY3BKScrRcyijiSnY/TNOuonPxgIrqi5Wr4s7if+EA2wsqWFToA4CpRsV/dr0P1tacm/gsHsEBBDjbEGCrxsnOSubeEULUC3KmpiHIT4ZjK9i2xZyL2q50sPqV3nZfQYuBEDYZY8tIVCamSqcUDVy+rpxLOTrSj5xEf+QwFy1diHFoQVKODk1pCT9vfBWAO4e9SYmpJQDj4nfywIUt7A+JIHb0ZAKcrQl0sSYgIxHfdi2x8XCT63eEEP9Koz1Tk5KSwv33309WVhYmJia89tprjBs3TulYt14zf4iYTf+eZfhujMY9Nw/SjZCwg4LYc/xaqKNt6xI63j0QVTNfpdOKBqqZtRlh1mbg1x/G9gfgJX4fzsot4upAV3Lik7grsDWXcnRcytbiqcvFxGgguUzFmhOpAJhX6lm74RWuAk+On4ertxuBLta0L0zBt6II9/BO+LRpIcNZQog616DO1KSnp5OZmUloaCgZGRmEhYURHx+PtbX13/5sgz5TU5vcRDi+ggPRJRwvGIqv2TGinN6BlndA2CSMzQeiMpEhAXFrleoruBJ/mctF5SRUmHEpW0dBYjL3/TQfswo99wyZU7PvUydWM+TyIb4NHsSqdoPxc7KmhYMpUft+xDwgALsJdxPobo+TtZmc3RFC1Gi0Z2o8PDzw8PAAwN3dHWdnZ/Ly8m6o1DQ6Ts1h0Bt07lmCw8Zo7NMTIM8A8Vsoj/2NH/I+wN+7mB5jgzEJ7Alq+b9iUfcszU0Jbt+CYOCOmq0h8NJI8otK+KVAz6VsHUk5WuyK/EkqyyHF0YuKKiMJWVoqEtJ5et9WtIcsGFfcHFQqbC1MmJocQ1BhKvkDhmHXpzcBztb4O1lhbW4ihUcI8afqtNTs3r2b+fPnc+zYMdLT01mzZg2jRo26Zp9ly5Yxf/58MjIyCAkJ4cMPP6Rr1643/V7Hjh2jqqoKHx+fOkrfMJlaWdF6XBQQBTkz4NhykvZdprjShZQUPZqvh4OdJ7QbQ2nzMVg27wTypSBug2Z2VjSzs6peFR0gcg4Agw1G0gpLScrRkRp7ifPG0RTq9Hg7WpFaUEpxWSVuCWcIyEnkZ8vmRKdWT1boXZzF4j0fkukewL5H59Tclu5PSfVkg7Y2Cn1SIUR9UaelRqfTERISwpQpUxgzZsx1r69atYrp06fz8ccfEx4ezuLFi4mMjCQuLg5XV1cAQkNDqaysvO5nt23bhqenJwB5eXk88MADfPbZZ3+aRa/Xo9fra54XFRX9249X/zm3hMi3adGvHPPf9mJISkGVZw/FaRj3L+Wn9W0xMYkjckAGjj2Gg0uQ0olFE6RWq/BuZoV3Myto6QJR1QtKTKF6odAreSWkdtQQf/4Cfp5t6WywISlHh3d6FtblpVQVFvLdoSs1x3t736d0yo7nq36TyenWnwBna1pYGGiefQmvkDZ4tm4uZ3eEaCJu2TU1KpXqujM14eHhdOnShaVLlwJgMBjw8fHhySef5KWXXrqh4+r1egYNGsSDDz7I/fff/6f7zZkzh9dff/267Y3mmpobVamHi9spOLyV7w8MR62qZIrrJExV5eDegTyfezHvMBhrH3+lkwrxlwoKdFw+fYG0zHwu2PtyKUdHUo6Ox757ncCCVGb0fpxzTgEAdMk4zxsHv+SSnQczIl8gwNmaABdr+sTvx8XalGYR/QkMDsDeSu4aFKK+u5lram5bqSkvL8fKyoqffvrpmqIzceJECgoKWLdu3d8e02g0cs899xAcHMycOXP+ct/aztT4+Pg0vVLzP/T5+eQc2IVX/veQuBMMlfya9xpXykMZELCF1gNaQZtRYO2sdFQhbpjRaCQ9JYvkEgNJheUkZetQHdxL150/Em/jzvxOd9fs+8X2uXjqcnmh5yOccWmBo7UZvSsyGHQhhsr2oViMGUuAsw1+TlZYmGoU/FRCiD/UywuFc3JyqKqqws3N7Zrtbm5uxMbG3tAx9u3bx6pVq+jQoQNr164F4Ouvv6Z9+/bX7Wtubo65uayb9L/MmzXDa+gYYAzocjGcW0flag2Uq/Eo3gAbP4NNL5Drfie5jkMJuKMfpvbNlI4txF9SqVR4+rrhCfT4Y+PwNsBD9KsyMDS/tHqywSwtBTnd0V29jN7TFyogT1eOMfE8zc/sZ19OEW+VBP5+TFi8/yMsTE04Pf5RnNsEVZ/tsTXBy9EaEwv5u0WI+qhB3f3Uq1cvDAaD0jEaB2sn1F2nMLoraK8kY3P1cTizGtJPcT7OltMl9rQ6sJiB3ZKg3djqW8VNLZROLcRNMdWoq8uIszUDWrlBn/kAjAB0+kqSc3WkHrXj4n4HSswdCHV14FK2Fm2JnsCsJEyMBl45k0NOYvV1fsOS9vPoqTUcatWDQ2MeIcDFmkBna5pnJ+Pd0g/XQG/UGjnDI4RSblupcXZ2RqPRkJmZec32zMxM3N3db1cMUQsbX3/wfRJ6PAk5Cdj+vAe7czk0N98N54/D+XWUmHhxiOm07BGAd9++oGlQfViI61ibm9DW0562I3rDiN4APEX1cFZuUSlXIpxIi01gXIuOJOVW363ldT4PDUbSK03Y9vtioWqjgTW/vky+oZJ7hr2GTYAvAc7WhOjS8Ndl4RLWEf+wtthayPU7Qtxqt+2byczMjLCwMKKjo2uuqTEYDERHR/PEE0/crhji7zi3IPThFoQYDJARBGd/grO/kJjegfPFgeSsuci4Yw9XX3vTfhz4dJVbxEWjolKpcLa3wvmOnnBHTyL/57Wqp3qRmphCv/wSfA2WXMrWkXE5nTwbR+x1BSSb2GK4Wsjpq4U4n/2V7gm/saZ5bz5tPxIXW3MCnSy5d8+3mHj7YDr+bgK8nfBxtMLcRM7uCFEX6rTUaLVaEhISap4nJSVx8uRJHB0d8fX1Zfr06UycOJHOnTvTtWtXFi9ejE6nY/LkyXUZQ9QBlVoNnqHVj4g3cDu8nzY7L+JWsg902XDkMwyHv2BN4QK8/NR0Gt0JM78OUnBEo6bRqPEN8sMX6FWztS08GUFpWTlbC/UkZutIztVhZojlkj6LDLfqO7Kyi/UYMzNpeWwXlcfVjDS2xaDWoFbBlJS9hGXHk9YzEpOIO2omG/S0NcPETM7wCHGj6vTup5iYGPr373/d9okTJ7J8+XIAli5dWjP5XmhoKEuWLCE8PLyuIvypRrdMglKqKuBSDJz5iZQTyazPfhELVRGTXKegcW0J7cZS2Wo0Jm4tlE4qRL1RVFZBco6OK4lXqfp1LSV5BXzfaRSXsnVo9ZW8emg5PdPP8lH7UaxvXl2XmpUVsXLrW2Q7uLHq0XkEuNgR6GxNYFURPm52uPh5oZaZwkUTUC9u6a5vpNTUvQqtluTtMegvnaBdyRKoKgfg59x3MJrZ0ndABS69h4Kdh8JJhaifjEYjOdpyLh86Qd7pc8Q5+XFG5UBSjg6b2NO8s/s/pFk7MXXQzJqfmXl4JX3STvNVx9HE9RxaPdmgtZq2l0/i0joI/15dsDGXa95E41Evb+kWjY+pjQ0tRw8HhkPZU3BhAyXHNpCZEYSxQoXVwWlw5Dnw74U2YDwmbSOxcHZVOrYQ9YZKpcLF1hyXiG4Q0e1/1s+CyspepF4chkVyOnPsPEnK0XEpR4fdMSNVKjVJ5s04dbWQU1cLaZObxMA9y8iwasbwO17B1dacAGdrIi/tx11dgcWAgfiGtJLrd0SjJ2dqRJ3TpaWSvieGFtoVkHIIgOjCJ4gv7UvvFgdoF9kWggaDuazVI8Q/UVpSxpW8EpILy0nK0aE9dpw2m74jzcSGt0In1Oy3ZNf7tCxMZU74ZA55tEWtgu6GXO49tR5dizaU3jul5pZ3DzsLNBoZzhL1j5ypEYqy9vSixfh7gXsh/zLGMz9TuM4BAyY45m2GnxeBiQVF3qNItx6M/8DemDvKLMZC3ChLKwuCrSwI9v59Q9/mMH0cAHf9fv1OUo6OKu4g7tJFzINaYlNuglZfiUVKEgHJ5zipLWOWzX+vZ1y4ZxmOVWX8NvIhLDuGEuBkTYA1+FqpcPH1lOt3RIMgpUbcWs38UPWZzpg+kHvmDI5Zw+FsGeQncfFsOQe19vjuXU5U5/3QOgpaDQcbF6VTC9Fg2VmY0sHbgQ7eDhD6AgCjqL5+J1urJ/msH1f2eaKvNGGQpxtJOTou52jxKcrAtqKUXSlaLhcmAtAv5TgvHvuOne7B/DB2xu9ndWxoc/Us7t7u+HVpj529nHEV9YcMP9UBg9GAWiX/F3PDjEbIPMv5DQc5ecaGEPOfaWu1HYAKLNhc+g5+wVa0G9ULjaOPwmGFaPwqqwykxiaRejaOeLcWJBWWcylHR/PffmXs4V/Y5teVJR3H1ez/w6ZZ2JeX8Hi/Zyn2CSTA2Zqupem0yYzHqmMoXn174usk1++IuiHDT7eRwWjguZjnaO7QnMdCH5NycyNUKnBvT5tp7WkDGLPCIK4HXFhPyiUzUgoDKDiWSYfL7cCrE7SOojwwCjOvlkonF6JRMtGo8WvbHL+2zf+7fhbA1HBKda8wLD2PtnoNSTk6rqTlk3fEB0NeOmk2zpQV68ku1tM8djcBsVvZduwck45XoVaBVzNLnt23HBNHR4rHT8I70ItAZ2s8HSzRqGVOK1H3pNT8SwfSDrDjyg6S9x0gMT2Wtwa9i7WptdKxGhSVaxC4Tofe03G7come2w5hkn0alV4Facch7Tg//eCE2vQIg/pm4dQ9AtzaykR/QtwGltaWBLfwIvh/N969FoDD/3P9ji46l3gzHfmuLbE1N6FYX0luZh6t444AcKdjb0pMUwG4K/E3hqQcIbbzAAqGjSXQ2Rp/J2v8VKW4+sn1O+Kfk+GnOvDtmvfJ+CGaIqsKzg2z5P3ID/GxlWGTf02bBbEbKD65i6+PPoAKA1NcJ2GuLgHHQLLd78bg1xfXzp1RySKCQtQbf8y/k5SSTcHmzWhT0tgaNoykHB3JuSU8fvg77rhylBWtB/NDcAQAFpV61mx4hTITM96evBBPTycCnK1ppc/B29yIT2hrHFwcFf5kQgky+V4tbmWpyb6cxOp3Z3HGKYu9LdKwN7dnYd+FhHvc+pmSm4qy7CyyDu7BV/szJEZDZRmb81/kkr4b3ZzXE9bDrPpCY78eoJaCI0R9VWUwcvXiFdJOneeyqR2xmurJBssuXuS1NW+jNbNiwtDXa/Z/6sRqhlw+xDfBg9jaJap6skF7U3qd2oFti0Dchg3Gz9kGC1P5776xklJTi1s9T02ZVksBWp6NeZazuWcxM2iY3nUG97S5F5UMk9QtvRYSthP9UzaJGV6MdnwZF9NkALI1HTmrvp8W4X749OsDJmbKZhVC3LCy0jJSLl4hSWNH0u/DWsE/fUZw7BE+bTOMGJ9OAPgVpfPxzoUUm1py19A3UKlVeNpbcteV/QSUZKPvfwcu3boQ4GyNl4MlJjL/ToMmpaYWt2vyvbLKMl7fPwftz4cwrVRjd2d3Xuk7G1ONLEp3K1SWlKC5shtV7K8Qt5EDWUM5rhtLc/P9DHb/FIIHQ+soqvz7o7GUW0+FaKiKyyq4nFvCpRwdmWficF3/HUUVRt4LGUdxWSUAc/d+TGhOAvM73c1O3zAA/EuymX1oOflegcROe/73xUKtCbACV7dmcv1OAyClpha3c0bhnJTLrJj5FFWVlWzploFPcDsW9VuEk6XTLX3fJq+qgoz9e4ndm4Rf+UYCjDsAKDPY8HX2x/g4ZxMx2haTNpFgIbNKC9EYGI1G8nTVMyvnbNpM+bnzHAjqxgmjPUk5OjqmnGb2oeXEO3jzdL9nan7u7X2f0iYvmR8GTUXXtVf1YqGWRvxLsvHr0JpmHjIhaH0hpaYWt3uZhKzkS+w7vZ352hUUVxTjbu3Okv5LaO3U+pa/twAMBrh6BC6s5+KhFLalTcLJJJkJzs+CxgwC+5HheCf2nfpj6eamdFohxC1gMBhJS83m6qHjZBaUcMqtFUk5WpJydMz64TU8dLnM6P0455wCAOiScZ43Dn5Jop0ns6JeqllColvSUZwdrHDr0xP/QC8szeT6ndtJSk0tlFr7Kakwiad2PkVmVgqtUx2YMOkFhgQOvW3vL8BoMJB9/Dj6+P345H0DOXEYjSpWZH9KiaEZo9uuxqNrZ2g9HOw8lY4rhLgN9KV6Us4nkKy2JVlbyaUcHRb7Yhi48zvONPNjXpf7a/b9cttcPEpyeaHXo5xxbo6nvQU9KzPpdfEAtAvBdkQUAc42eDezxFSu36lzUmpqoeSClgVlBSx9YQqWmeWc9ysiZMJYnuj4hEzUp5TsOHTHN7FhiyuFZQ5McZ2IiaoCgEvW91Jo35PmA7tg1zxI4aBCCCVoS/RcLigjOaeEpOxivL94H4u0FN7pNokrWAIwInEPj55Zxz6PdrwVPgkAE7WKdw9/gbm5KRfHTsOldcvqOXgczHBzsEYjMyz/I1JqaqH0Kt3n9uxiy/IP+TksiWLrSvr79Gdu77kyUZ/Cyq4mYnFlM1z4FVIOsS5vNlfLQ+huu4JOzROh9QiMrYaDSytUckGhEE1evq6cpFwdGQeOUrlvD0lWLmzxCCU5R4deX866X2diYjRwf+Sr5Fg6ADDs0j4ePLeBQ617cWr01Jo1tALzr+LXyh9HLze5S/YvSKmphdKlBqCyvJzNKVuZs38O5YZyWps3Z+HQD/Gxk4n66oXiDM6t/Y2Es6X0M3sXe00aAFf17YnRPUnrlkWEjWgHPl1lLhwhxDWMRiMZeTpSdh8gLy6Bo216k5RXSlKOjjtivmdkwm5+atGXL9pFAaA2Glj760xMDVU8GTULG39fAp2t6VCSgZ8+F5eOHQgIbY21uUz8L6WmFvWh1PzhTPYZXlv7LN12mXGpeTmTH3mD7l49/v4Hxe1Tkgdxm+HCenYf8+OMbjCtLbczwP4/YOUMQYO5aj0ct+69MLWxVTqtEKIeq9CXkxKbREpRBQkqa5JytGRfTuf+lXOwLtNy5/C3Mfx+OcK0s79yZ8JvrGnem0/bj8TNzpwAJ0sm7PsBUx8fzMeOx9/LCV9HK8xMmsbZYyk1tahPpQZg188rOP7jatKcStnZNZcZ4c9zT6t75BRkPVReWMCVmL3Y5e/BNeMHKCtAW+XEiuzPMVGVMbnHD5i1HQRBg8HGRem4QogGpERbyhVtBUnZOi7l6LDc8DMeJ/axxS+cdW6hALiU5LNy29tUqtSMjJqLQa1BrYIpqfvomJ1Aes87MBsY8fscPFZ42pqjMW08Z3ik1NSivpUagHP7Y/het5l1aZsAGNNyDK+Ev4KZRmbBrbeqKuDKAdL37Wf7weZYq7K40+nl319UcUzzNCq3IIIjwrAObKVoVCFEw1ZYUkFSro6UhBSM636mJL+IbzqPJilbh668itcOLadH+lk+aj+K9c17AeBUWshX298hy8Gd1Y+8Q4CLHYHO1gRUFePrZo+Ln0eDm3BQSk0t6mOpgepx2JXnV7Lo2CKCkqxxc/VhzqRlOFvKxE/1ndFgQJ98FouULRC7EUPaab7MWo7eaMtox1fw9KyE4KFUNh+KxrczKrnzQQhRB4xGI9nFei4fOEb+6XPEOQdwWtOMpBwddrGneXv3f0i1dmbaoJdqfublwyvpnXaarzqOJrbn0Oo1tKyg7ZUzuLRuiV/PzthZ1M+Z76XU1KK+lpo/bNn3I2eXrECFiv399bw57gPaOLVROpa4CZU5KVzYuJ+r8YVEmr2K2qgH4Jh2DKdKR9I56BIdBreGwH5gaqlsWCFEo1RZUUlqXBJXr2Ry0d6reg2t3BKGf/ceba6e5/VuUzjiXj0JbNvcJBbsWUaGVTMm3/EKzjbmBDhbMfjSAdzVFZgPGIhvSCt8Ha0UXTBUSk0t6nupMVRV8euXi9lzaRdbgy9jYWLBmz3fZHDAYKWjiX+irBASdkDsJtbtDeVqWVv62n1EO6ttYGpFhf8gEjSj8O/fU2Y0FkLcFqW6UlLyS0kqLCc5R4fu6DFab/qWNFM73gydULPfkl3v07IwlTnhkznk0RaVCrpVZnP/6Q3oWrRCe89U/J2tCXS2xtPeAtNbfBZaSk0t6nup+UNRWSEv7n2Jval7URtgis89PDnwRZmorwGrKisjbd9+nIt2YJm0DoqucqmsK5sLZmKvSee+sG+g1VAIHgpOzZWOK4Rogv53wVDD99+gunSRNaFDOVZuRXFZJRFXjvDc8VWcdG7BzF6P1Pzcwj1LcawsZeeIh7EIDSHA2YYJXXxQq+vuphcpNbVoKKUGoMpQxQcnPuDCD+tonmpN3kB3Zt2/BBszWWW6wTMaIeM0l3Yc4PAxW7zVR+hl91XNyxt1b+PoaUPosDZYtugMDeyCPiFE4/LHgqHJZxPI23eAtCpTDnq2qx7WytbyzbpXsKks45EBz3HZzgMnazOOvTaoTjNIqalFQyo1UD1R32ezn0B3KZXoztmYt/RkyYAl+Nr5Kh1N1CFD7mXUCdUXGucnJvFd1geoqWCK6yTM7ewgeDA6r2GYt+qFiZWV0nGFEKJGVZWBtITLpJ6JI8GtOZcKKjDRqHh5aN0u3CylphYNrdRA9XU2u/ev5830pWSVZmFnZseCvgvo7tld6WjiFigvyCV5x16KLiXQufJ9KC8GYFP+i6SUhzKg9UFa9m0DLe8AK0eF0wohxO0hpaYWDbHU/CG7JJtnYp4hNu0MfU66EDJuDBP7PCYT9TVmlXpI3oPxwia+29aFggo37nKajotpEqg05LpEcdV8MAH9OmPXvKXSaYUQ4paRUlOLhlxqAPRVeha+9SDm5/PItylHNbkrr3WfJRP1NQFGg4GcE8dwLtiGKn4TZJ7lQPG9HNeNpYXFXiKDttRcaGz0CJWFN4UQjYqUmlo09FIDoCvI58v5z/OL+0ly7PSEuISwuP9imaivqclPJn7zXs6frKKtejUtLfYAUGaw4ce8xfh5FdNrtC+aFr3BxFzhsEII8e9IqalFYyg1f9ifup8Zu2dQXF5MyzJ3Zg2bS6h3Z6VjCSWU5MHFbRC3ibhTpezIfRQnk2QmOD8LZrbQYgAZDiOwD+2Dpaur0mmFEOKmSampRWMqNQCXiy7z/PonCNlSSbmZkU6PT2Zk2HilYwkFVZaUkLp3H4YrRwkoWA7aDIxGWJH9GSWGZoxu/T0enUMgaAg4twS5JksI0QBIqalFYys1AMkXz7Jq7isUmJSwuVsGk0Km8mTHJ9GoZY2hJs9ggLQTlJzewfptPhTq7ZnqOhETVQUA8Zo7ybLoTXDv5rh06Q6a+rnmixBCNOpSU1JSQuvWrRk3bhwLFiy44Z9rjKUGoDg/l89OfcpXl78DoI93H+b1moetua3CyUR9UpaahEXqDojbDMl7WJc9k6vlofSw/YqOTr9Bi0FUtRhCpW8/zB3lGi0hRP3RqEvNK6+8QkJCAj4+PlJq/sfGSxuZvX82PskmtMpz5p5n36S9T0elY4n6SF/MpW27STqRSZjmUxwqLgCQou/AhvzXaO6UwB3DjRA8BBwDFQ4rhGjqbub72+Q2ZaoTFy9eJDY2lqioKM6ePat0nHplWOAwvEzc2LJlNqYV8NbyJxkyehoPtH1A1o0S1zK3JTBqGIFRgGEiXD0KcZvI2FOOARNMSq7C1mWw9WVwDuaUehquIW1xCw9HbdKg/soQQjQxdfZtt3v3bqKiovD09ESlUrF27drr9lm2bBn+/v5YWFgQHh7O4cOHb+o9ZsyYwdy5c+soceMT6tuZ0a+8TnFrW877FLLw2EKmbZtGujZd6WiivlJrwDccBr1Olzfmct9zXnQaHAABfUBtgjYzh72nW/DL16WUzQ+DNY/AubUYy4qUTi6EENeps//t0ul0hISEMGXKFMaMGXPd66tWrWL69Ol8/PHHhIeHs3jxYiIjI4mLi8P191tNQ0NDqaysvO5nt23bxpEjRwgKCiIoKIj9+/f/bR69Xo9er695XlTUNP4SDm7dmdmzv6PdxZ9578h7HE0/wjuvP0DfqHsY22ey0vFEPWffMhhaBgMPQmkBlcd20XJHMuWFhVjpk+FUMpz6nl1FT6I1a06XHho8+vQDB1mTTAihvFtyTY1KpWLNmjWMGjWqZlt4eDhdunRh6dKlABgMBnx8fHjyySd56aWX/vaYM2fO5JtvvkGj0aDVaqmoqOC5555j1qxZte4/Z84cXn/99eu2N9ZrampzpegKiz6Zjs/RMspMqyia1JaXe7+Gvbm90tFEQ1NVCSkHIW4zhtitfBn7GnqjLaMdX8bT7AK4tUPrMwKtS1/cwrqgMpE78IQQdUPxC4X/f6kpLy/HysqKn3766ZqiM3HiRAoKCli3bt1NHX/58uWcPXv2Ly8Uru1MjY+PT5MqNQDFhXl8segF9pidJ967GDcrN97u9TbhHuFKRxMNlNFopCD2Apf3nqCD+lvUVw+A0cAx7RgOau+npc1h7uh9tfpC48B+YGatdGQhRANW7y4UzsnJoaqqCjc3t2u2u7m5ERsbe0ve09zcHHNzmSLe1t6Rp+d8xoDsM7y872UuF13muZ8fZaRlX56cOBcLEwulI4oGRqVS0ax1G5q1bgPc+/usxtup3JKGqa4UT/UxOLENTnxNhdqe7frZ+Laypc3wnqibeSkdXwjRiDXIWxkmTZqkdIQGRaVS0cG1Az8O/5GFh+ZTumIv6qKLvJBxF48/Op9gx2ClI4qGzMoRQsYTHgKdy8owJltCUiDEbSYlw42kggByD2bQNqENeIZA8FB0nndg1UIW3xRC1K3bUmqcnZ3RaDRkZmZesz0zMxN3d/fbEUEAVqZWvNrjNX5MWsTFHTs56HyVvRvv5qmOT8mt36JOaCwsoNWA6sfgeTjHn6XbzlOY5p1FVaGC9FOQfoq12Z5UqJIY1v0MLl16VN9tZSpnDYUQ/85tKTVmZmaEhYURHR1dc02NwWAgOjqaJ5544nZEEL9TqdWMv38GmSPuJfPYO8SkxLDw2EKOH9jOjFFv4esSoHRE0VioVNgFtycsuH31c+1zcHErutMxaDOdMRg12Md9DBffB1Mrrja7hyKHHvj3746Vh6ey2YUQDVKdlRqtVktCQkLN86SkJE6ePImjoyO+vr5Mnz6diRMn0rlzZ7p27crixYvR6XRMniy3GSvBzd6DJf2X8MvFX/hs20J89hawfP9jdJg+hZHtxqKSxQ5FXbNxgY73Yd3xPqaOKyHn6D7Miu6B+C1QlMqZOCcu6ZvR5fgcugYnQvBgjC0jwa2dDFMJIW5Ind39FBMTQ//+/a/bPnHiRJYvXw7A0qVLmT9/PhkZGYSGhrJkyRLCw2/PXTiNfZmEf+PEiRg2f7iQNDstv4XmMDhgMK92e1Vu/Ra3h9EIGWc4se4ECfFq+lkuxMU0CYCM8iC2Fr1IkF8e3UcFgH9vGaYSoolR/Jbu+khKzV/TaQtZeW4ln8R9QZWxCndzV2a2epYBocOVjiaamqJ0uLgV4rdy6LgzR4tH08JiL5EOC8HUCgL7k2Q+EtfwHlh7+SidVghxi0mpqYWUmhtzNucsM/fMxPlQIa0v20JEEE9Pmoe5Rm6PF7dfpU7L1T0HsMjej3v2KihKpcxgw5dZyzGiYVK7RVi36wVBkeARCjJsKkSjI6WmFlJqbpxWr+U/sx9Bk1TAjrAsLIK9mNd7ntz6LZT1+zBV7pHf2BnjQFV5BROcn615+VD5g5RYt6F9X0+cu/YBMysFwwoh6oqUmlpIqbk5RqORjb99y4LMz8kty8VUbcrjwQ8xqfODaNQyBb5QXmV+BibJOyBuM8aEXaxMW4zW4MzwZm/iZ30eAvpS5jeESp9+2PjJXX1CNFRSamohpeafySvLY87+OexJimHEXg8M7jZMffZd/F1bKB1NiBrG8jJSftvH5eOX6W6yFJPiRABO64ayp/hB2jgepX9EOQQNAc+OIHdTCdFgSKmphZSaf85oNPLthg9I/3Y7JeZV7BhQxAu9ZjI8cLjc+i3qH6MRsi5A/Bb27ijnVGZnetiuoKP1egCqrNzZU/E8vh088R/YG7WVrcKBhRB/RUpNLaTU/HsnT+3hoxPL2K86B0CkfySvdHmZZlaOCicT4s+VZmaiuhyDRfImSNxJSpEf6/PfwFJdwGT3R1AF9IKgwVT4D8LUTYaphKhvpNTUQkpN3ag0VPLFmS/46NRHuGWZ0j3OhZ4PPsigrqOVjibE36ssJ//EPs7uSsasMI5w049qXvoxZwFVJtYM7JmGa3hv8O4Mcv2YEIqTUlMLKTV160z2GX5+bSa2+XDev4igO4fxdKenZdVv0XAYjZBzEeK3UHYuhq+OP4IBNZNdpmClKQRLRzJcJqB17IFvv56YOcgZSSGUIKWmFlJq6l5uXgaff/oq37sepkpjpLl9c+b1mUcrx1ZKRxPippVlZ5FxcD/++l8hYTuUFbKj4CniyvoTar2enh0uQdDg6jlxnJorHVeIJkNKTS2k1Nw6u6/uZta+WeSW5dL9nBMdO/ThwbtnYaK5LeulClH3qioh5SBH18cSl2BDf5sP8DQ7D0BuhS9bta/Q3F9H+OhW4BMO8mddiFtGSk0tpNTcWnllebz744u4bkzHgJGE0c2YPfw9PG1ktWXR8BlzElBd3Abxmzl21pWDxffiZ36U4c3eBgsHaDmIVOvhOHXpiYWTi9JxhWhUpNTUQkrNrVdVVcnKFW+zOzmGo81zsDG14eXwl+XWb9Go6PPzuBKzD8ucw3jnfQOleVQYzPkiayUGNNwX8iV2HbpXD1U5t5SlG4T4l6TU1EJKze2TUpTCzL0zOZV9Cgu9muE5HXj88QU427spHU2IumWogqtHyDv6G1t/86WiEu53fqSmx5wwTEJrE0qbPn44de4FJmbK5hWiAZJSUwspNbfXH7d+n//iR3wzLcnwrGLM87Po7tld6WhC3DL6tETMU3ZA/FZI3sN3GQvJr/Im0mE+LezOQPP+lPsPodK3P1YeHkrHFaJBkFJTCyk1yti9bx2/ffkZ20PTyLer4L7W98mt36JJMJYVk7R9L8mn0ulp8j7mZVcAuFAygJ1Fj9Pa6TQD7qiqvpvKvb0MUwnxJ6TU1EJKjXJ0ei3vn1jMqrhVAHTWBfBYn+l0addP2WBC3C4GA6SfhPit7NtVxcmsHnSx+Z6uNj8CYLT1Zl/lM3h38MWnX080ljbK5hWiHpFSUwspNcrbc3UP7+yYTZ/tFmgMKszuDefRIS9irjFXOpoQt5U25Qrq5N+wuroZLu0iQ+fFz3nvYabSMcXzYTSBvSAoksqAQZg4+yodVwhFSamphZSa+iE1I4nlC16kuCifTT0y8LP3Z3b32XR276x0NCGUUVFG/vG9nN6Vgrogid4WH9S8tCb3TfQaJ/r1yMK9ey/w6iRLN4gmR0pNLaTU1B9Go5FtcZt49/RCskuzwQh3l/bhkXvn4Ggnc3yIJux/Vhgvv7CTL44/iQFT7nd+GDuTLLByJtttHIUOvfDt2xMzh2ZKJxbilpNSUwspNfVPUXkRi44u4uyOrXQ/50SRvYG+Lz/HAP+BSkcTol4ozcwkbf9+mleuh4Ro0BfyW+FDnC0dQjurrfQNuSBLN4hGT0pNLaTU1F9bd33P0ZXfcMo/nwv+xUT6R/JS15dwtnRWOpoQ9UdVBVw5yPENF7hw0Z6e1h/jb34cgKJKV34tepNAvxK6jW6Byq87aEwVDixE3ZBSUwspNfVbUXE+n8V+ydcXvqbKWIWn3oEHXMZw98inUavVSscTot7579INWzh93pY9hVPxND3HaKdXwdweWgwkzSYKx869sHBxVTquEP+YlJpaSKlpGM7nnmf2vlkEbCrALd+CnM52PPHofLxsvJSOJkS9pc/PI+W3/ZhmHcOvYAWU5FJl1PBl1goqjBaMb7scp9Cw6qEq19YyJ45oUKTU1EJKTcNRVl7K55++RtHB86zvlYbB1ownOz7JPa3uQSN3fgjx1wxVkHqcohO72LTTm5JyCya7TEGlqv6r/ozhHvKtutC6ly8u4b3AVCbCFPWblJpaSKlpeBKy4nj7+DyOZh4FoF9RK6YMeoaOrXoqnEyIhkOfnoz51ejqpRuSfuPHjDfJrmzBALsPaW1/EAL7UxE4mHKvflh7y5w4ov6RUlMLKTUNk8Fo4OeLP/PZzg+IiLEHFVhO68PD/Z7FTCOLAwpxU8pLSI7eQ9KJVLpqPsG6NBaAhLIebC14nhYOZ4mM1FXfTeURCnI9m6gHpNTUQkpNw5Z4+RzfL5lDjj6X6M7ZBDoE8nqP1wl1DVU6mhANk9EIGWcgfiuHY3QcSe9DqNVaetqtqH7Z2o39hul4tPHFb0BPNNb2CgcWTZWUmlpIqWn4DAYDWy5u5N2TC8gry0NTpeJuXS8envwWDjaOSscTokErSU/DmPgb1qkbIXEXeSUOfJ/zIRrKmeoxDdPALtDy96Ub3GROHHH7SKmphZSaxqNQX8j8I/O5vDmG0AQH8p0M3DHzJfr49FE6mhCNQ2U5Raf2cir6MlV5qfSzmF/z0ob8VyjChz5d0/Du2Q18wkFjomBY0dhJqamFlJrGZ+P2lZz8dhUHW+WQ7FHCsMBhvNjlRZpZyNTxQtQZoxFyEyB+C1Wx2/niyENUGK2Y4PQ0TqZXwMKeXPc7ybHph1/f7jInjqhzUmpqIaWmcSoszuPTC5/zTey3GIwGArXOTPAZzYSoJ1HJXBxC1Dl9Xg5X9xwgsHIDqoRtUJrH/uIHOKEbTZDlbwxqd7D6QmOZE0fUESk1tZBS07idyT7DnL2zaPNrCQ46U7J7NuOZqfNxt3ZXOpoQjZehCq4e5fTm05y/YE2Y+de0tNwPQEmVPT8XLMDPS0uv0T6oA3vLnDjiH2m0pSYpKYkpU6aQmZmJRqPh4MGDWFtb39DPSqlp/ErLSvjsk1coOH6Bdb3TMbG04JmwZxgfPB61Sm5NFeJWM+ZfqT57E7+N2LNGovMfxcUkkbucZ4CpFQT2I9V2BA4hPbD28VM6rmggGm2p6du3L2+99Ra9e/cmLy8POzs7TExu7AI1KTVNR3xWLG8de4cTWScAiMhtzeSo5+jQMlzhZEI0HRXFxVzdsx/SjhNQuBKK0zAaVXyV/QWlhmbc2fIz3Du2qx6m8uwoc+KIP9UoS825c+d4+umn2bFjxz/6eSk1TYvBaGBV3Cq+2baMAfscqFQbsHv8DqZ1fxxTWb1YiNvr9zlxSk5Hs2mbM/klDkxxnYRGVQnA+arRpJv1oXV3Vzx79QEL+Tta/NfNfH/XWTXevXs3UVFReHp6olKpWLt27XX7LFu2DH9/fywsLAgPD+fw4cM3fPyLFy9iY2NDVFQUnTp14p133qmr6KIRUqvU3N3qbhZH/YcSLwsueen4T+IXjN84nrM5Z5WOJ0TTolKBRwesIp9l7ML7mfhWFzSjP4Q2o8DcjvjCjsSm+pOz9Vt4LxBWjqRq38cUxscqnVw0MHU2uYBOpyMkJIQpU6YwZsyY615ftWoV06dP5+OPPyY8PJzFixcTGRlJXFwcrq7VtwCGhoZSWVl53c9u27aNyspK9uzZw8mTJ3F1dWXw4MF06dKFQYMG1dVHEI1Qy4AOvLLgRzZe/JULJxZyMf8ik9bdx3hdTx6e+jZ21g5KRxSiyTFzcgeneyD0Hqgsp8ve/TgfScJfnQ3FFXAphqsXitiQH4S39W+MvONK9TCVbzeQM63iL9yS4SeVSsWaNWsYNWpUzbbw8HC6dOnC0qVLgerZYX18fHjyySd56aWX/vaYBw4cYM6cOWzduhWA+fOrJ4N6/vnna91fr9ej1+trnhcVFeHj4yPDT01Yflk+7x55l7x1BwhOsSXbw8DIF16ju2d3paMJIf6QkwAXt3Lqt0z2J/eitWU0/ew/rn7N3J69xudxaeFDYEQPTB1lTpym4GaGn27LNJDl5eUcO3aMmTNn1mxTq9VERERw4MCBGzpGly5dyMrKIj8/H3t7e3bv3s3DDz/8p/vPnTuX119//V9nF41HM4tmzOs9j/W6Lzn9w88c9c1i4/aHGNViFDM6z8DeXNa2EUJxzi3AuQUh3aFVfh6VF0shtQgubqWoWM2p7LaorlThd7YTpn5tICiSct9BmPq0QyUXGzd5t6XU5OTkUFVVhZub2zXb3dzciI29sTFTExMT3nnnHfr06YPRaOSOO+5g+PDhf7r/zJkzmT59es3zP87UCDFi8BT69RmD5sxHfB/7PWsT1hJ3dD/jWtzJ2KGPyqR9QtQT5s0cMe86EhgJhio0sUcJ2x5LSWYWFupiSDkEKYfYVaAlozKePqGJBPQJhYA+YGqpdHyhgAa1YMeQIUMYMmTIDe1rbm6Oubn5LU4kGio7Kwdmhs9kSMAQ3oyZRfvocq7s38TMq0eZfv88XK3ktLYQ9Ypag3WbcLq1+X1qhsK74eI2jHFbSdvXjhJDMywvvQdX3wETSwrcR3DVPAL/Pl2w8QtQNru4bW7LuTpnZ2c0Gg2ZmZnXbM/MzMTdXWZ8FcoJdQ3l66jvsA0LJt+2gi2mRxi5diSr41djMBqUjieE+DP23tB5Cqp7V3Hf/EiGjarCtXtPsPOGylIuxhr5bZ8zexZ9Dx/1gug3IOUwxqrrb0YRjcdtOVNjZmZGWFgY0dHRNRcPGwwGoqOjeeKJJ25HBCH+lJWlDc88/QFx2Re4dPhNzuSc4Y39b3Dkl9VMHPsCbVt0VjqiEOIvmFrb4D94EDCoek6crPPYbTmE+5mr+KuPQOYZyDxD+W//4bvcZXi7FdJvhAMmrQaChVxL15jUWanRarUkJCTUPE9KSuLkyZM4Ojri6+vL9OnTmThxIp07d6Zr164sXrwYnU7H5MmT6yqCEP9KsEtrvh7yNd/FfseajZ/idqKE9WdmceCZEUwKm4aJukGN1grRNKlU4NaW4IltCQbQDYOEHRC/hZRTheiqmpGZWYLJ2smgNgHf7ly1HYN1m244BLeSi40buDq7pTsmJob+/ftft33ixIksX74cgKVLlzJ//nwyMjIIDQ1lyZIlhIffnqnrZUZhcTNiE0/w44dvctE6i5NBhQQ3C2ZW91l0cOmgdDQhxD9UVa4n4+AhypOOEVD8HeTEYzTCNzkfUVTlznCfz/AL84eWd4BfTzAxUzqyoJEuk/BvSakRN8tgMLD+4joWnFhIob4Q61ITRuWHMeXhN3B19FQ6nhDi38q7RPnZHWzZaEZGkTOTXKZipi4DILZ8MEmaSFp3tsV/YF+wkZsHlCKlphZSasQ/lVeWx6Kjiyj6cT9+mVake1bQ78knGRYwTG7/FqKRqNIVormyG+K3wMXtbEyZSrK+C91svibM5hfwCsPQPJJc+4E4d+wkw1S3kZSaWkipEf/Wtt0/cvCblWxvn0aBXQXh7uG80u0VAuzldlEhGhWDgezjR0k6EEtL4zqa5VUvpJxW3po1ee/gZJbK+IGHUQUPhsB+YGatbN5GTkpNLaTUiLpQXlnOygsr+fjUx+ir9LRJsaeHfWemTXsTa0tbpeMJIW6F4gy4uI3Y3Yn8FhdOoPkhBjksrn5NY85e40vY+XgRFNkZC8/mikZtjKTU1EJKjahLKcUpzIt5E+8frmJWqeZ8VyNTJ7xCD68eSkcTQtxCVaWl6OP3Y5W2DeI2U5qXz5dZXwFqJrpMxcbdFYIiKfOJxCywK2ozWYDz35JSUwspNaKuGQwGfv71Iw7HbGBLyFVQwRD/IczoPANXa7moUIhGz2ik7Eos57eeoCAlmwEmr4OxCoCdhY9xSd+dXq1O06pfMLQYCFaOCgdumKTU1EJKjbhVtOValp1cxnex32E0GIg86kFw155MvGcmpnJLqBBNR2k+JO6E+K18v6sXeRXejGg2Gx/z06BSU+R6BwnqKPx7dMCxXfvqOXXE35JSUwspNeJWu5B7gaXfvUbg3jLKTQycG2XDywPm0MapjdLRhBC3maGykszDh3Et3oEmcQtkneOkbgT7iifjbXaSkQFfQVAkBA3G4NsTtbkswPlnpNTUQkqNuB0qKstZ+e1cotN2ccYzF7VKzd2t7uax9o9iZynTsQvRZBWkkBS9h3NH9fgbd9DOYgMAlUZTVmZ/hluzAiKGqjBvfwfYyTxY/0tKTS2k1IjbKac0h/lH5rMpaRMORabccdyd4NHDGBf1mMxtI0RTV14CSdVz4qScvML61CewVucy0WVa9YiUewdS7O/CNLALbmFdUJlolE6sKCk1tZBSI5RwIO0Avyx+G/fLkOymo2JEMK90fQUfOx+lowkh6gGjwUDOqVPoYg/jr1sNV48CRn7IWURuZQCDXD8jqKND9VBV8wFg0fS+v6TU1EJKjVCKrrSILz6fzSrzPRSZ6zHXmPNgm6k80PYBLM1l0i4hxP/Q5VAVu50d60pJyXHmPufHsFAXA3BR34dzVXfSJkRDUGR3cG6hcNjbQ0pNLaTUCKUlFybz9qG3OZh+kI5xDrTItqfH5MlE9LhT6WhCiHrIUFGO+upBiN8K8VvZlhjFxbLedLL+ie6234Jjc4wtB5NqMQiP7t3RWFgoHfmWkFJTCyk1oj4wGo1sjP+Vk3M/xrJUTXRYFh17DGJ62HScLJ2UjieEqMcK42NJ3nsCb/12nLLXg6GC7IoAfsxdhJUmn4m916EOjoSWgxrVApxSamohpUbUJ9n5aXz+4zt8b7EbI0bszOx4suVDjOt4LxqNidLxhBD1nb4YEneRtPc0u061w93kAkObzfv9RRV7DC9i4eZNm0HtsG4Z2qDnxJFSUwspNaI+OpN9hjcOvsHF7FhG7vFEbW3BiKdfomPL7kpHE0I0EMbKKvTJx7FI2Q7xWylPjeWLrJUYMOVe58dxaGaElndQ6h2JScvemNo2rO9AKTW1kFIj6qtKQyUrty8la8U2KkwMrO+byYSQ+3g05FGsTK2UjieEaGAqclK5uPUAWQmZ9DN5Gyp0AOwpmsK5kki6BRwmtJ87BN0BzfyVDXsDpNTUQkqNqO8upVzg071L2Fi+FwB3a3ee9pzC8B53K5xMCNFgVZTB5b0Qv411OwK4WhLMEId5BFocAkDn0IVThvvw79oCz+7doB4Of0upqYWUGtFQ7L66m3cOvUPV5RwGH3JH52/FlFcW4W3nrXQ0IUQDZjQYyL9wHtucnZhe2gJXDnJON4CYosdwM41jrNc8aBEBQYOp9OuPib2z0pEBKTW1klIjGpLSylI+/uo1KqMvEO+j5VRoKY+GPMp9be7DVG2qdDwhRGNQmk/a7t84dzAXN/1+Opj9BIDBqGZ59hc4WBUTOagQ69AIcG2t2MXGUmpqIaVGNEQnzu/lw4ufcqTgBACtLJvzuO9k+nUbqXAyIUSjYqiCq0cgfitZJ8+wOuExzFXFTHGdhFplAHsfrjjci9GrM149u2NibXPbokmpqYWUGtFQGY1G1ieuZ+HRhbQ+oiLoqi36Hl489vC7OFg4KB1PCNEIFScnUnDyID66tdXrVFWW8XPuO2RUtKavw+e062D4fZXxyFu+AKeUmlpIqRENXX5JHh8vmI7puWw2dc+k0t2K5zo/x4jmI2SRTCHErVNegvHSbnb/kkpymgN3NnseG00uAJfKunJEP4nWwaV0GBYCXp1AXbcLcEqpqYWUGtFY7Du3gwUJy0goSACgb2UHpvV4nNDWPRROJoRo7IwGA6qscxC/BeK3setcZ86XDqK91Qb62H0B1q7w9Ckwq7vpKKTU1EJKjWhMKgwVfHP+Gz4/8jFDdjpiVqHGbEIXHh7+IpYmlkrHE0I0EaWZGVz57QCO2j24ZK4G5yB4MLpO30NKTS2k1IjGKPHqeb77YDbleYWs75WOh50XL4e/TB/vPkpHE0I0NVUVUJwBDj51elgpNbWQUiMaK6PRyPa4Tcw/u5gMXQYYYXRuJyaPn0mAdyul4wkhxL9yM9/f6tuUSQhxi6hUKu5oNYx1I9cxqe0kAjNssD+cy3czp7Pi5JdUGiqVjiiEELeFlBohGgkr0+q7oV4d+g46Fw0X/IpYcOp97tpwF8cyjykdTwghbjkZfhKiEaqqqmRN/Bo+OLWEAn0BNiUmRGW15/6HZ+Pj3lzpeEIIccNk+EmIJk6jMWFs63H8OupXxgaNpev5Zpifz+P9uY/y3YXvZEhKCNEoSakRohFzsHBgdvfZjH3gObTOag4HZTP38Fzu3ng3J7NOKh1PCCHqlAw/CdFEVFZV8kvCL3xw/AOKyotoe8mWEHVLJj76Op6u/krHE0KIWsnwkxDiOiYaE+4KvotfR//KaJ8RhF50wOx8Ls99NZkf436kylCldEQhhPhXGlSpef/992nbti1t2rThqaeeoomcZBKiTjlaOPLGgLfp9tTDZLbQcNYthzcPvsm9m+7ldPpJpeMJIcQ/1mCGn7Kzs+nWrRvnzp3D1NSUPn36sGDBArp3735DPy/DT0Jcr9JQyaq4VSw9sRRduZZh+9yx9/NhymNv4ep4a1feFUKIG9Foh58qKyspKyujoqKCiooKXF1dlY4kRINmojbh3tb3Vg9JmfbDucgcw4V07v51Ar9c/AWD0aB0RCGEuGF1Vmp2795NVFQUnp6eqFQq1q5de90+y5Ytw9/fHwsLC8LDwzl8+PANH9/FxYUZM2bg6+uLp6cnERERNG8u820IURecLZ15/d4PCXlyEpe6mZFFPrP3z+aBzQ9w4tIhpeMJIcQNqbNSo9PpCAkJYdmyZbW+vmrVKqZPn87s2bM5fvw4ISEhREZGkpWVVbNPaGgo7dq1u+6RlpZGfn4+GzZsIDk5mdTUVPbv38/u3bvrKr4QAojoNZaPHv+JGZ1nYGViRUrCBba//AbvzZ1GQVmB0vGEEOIv3ZJralQqFWvWrGHUqFE128LDw+nSpQtLly4FwGAw4OPjw5NPPslLL730t8dcvXo1MTExNaVp/vz5GI1GXnjhhVr31+v16PX6mudFRUX4+PjINTVC3KCskiw++uhFbA5nc8lDx9nuBp7r/BxRgVGoVCql4wkhmoh6d01NeXk5x44dIyIi4r9vrFYTERHBgQMHbugYPj4+7N+/n7KyMqqqqoiJiSE4OPhP9587dy729vY1Dx+ful0KXYjGztXKldnPfUXbx+4ls6sNeWV5vLL3Fab8OpGjZ2OUjieEENe5LaUmJyeHqqoq3Nzcrtnu5uZGRkbGDR2jW7duDB06lI4dO9KhQweaN2/OiBEj/nT/mTNnUlhYWPNISUn5V59BiKZqcN+7+X78LzzT6RksTSxhfxK73prP/GVPoi3XKh1PCCFqNKi7n95++20uXLjAuXPnWLJkyV+eAjc3N8fOzu6ahxDinzHVmDK1/VTWjVxHoMYbtVFFdMlBRqwdwaZLm2TOKCFEvXBbSo2zszMajYbMzMxrtmdmZuLu7n47Iggh6oCHjQevvf097Z6dhEkLN7JLs3lxz4s8ufIBTpzfp3Q8IUQTd1tKjZmZGWFhYURHR9dsMxgMREdH3/DkeUKI+iOy21h+GfELT3Z8EmujBS67stnxxjssXP0KJRUlSscTQjRRdVZqtFotJ0+e5OTJkwAkJSVx8uRJrly5AsD06dP57LPPWLFiBRcuXODRRx9Fp9MxefLkuooghLiNzDRmPNThIb69YyVqVztKzKv4puhXotZGsTV5qwxJCSFuuzq7pTsmJob+/ftft33ixIksX74cgKVLlzJ//nwyMjIIDQ1lyZIlhIeH18Xb/y1ZJkGIWyv6wmbmn/+AVG0qAMOy2nHPyGfoEHR7/hsXQjRON/P93WDWfvq3pNQIceuVVZbx5dkv2bjzawYccqJSbcDi8QE81O3x6junhBDiJtW7eWqEEE2DhYkFj4U+xoIRy9D5WBDvq+XzxBWMXDuS6CvRMiQlhLil5EyNEOKWMBgM7LoczXvHFpCmS8NCr2Z4UjB3TX2J1i06KR1PCNFAyPBTLaTUCKGM0spSPjv9GbHfrqH5VWuympXj98goprSbgoWJhdLxhBD1nAw/CSHqDUsTS57q9BSPPjQPnZc5h1vn8tGpjxi9bjS/pfymdDwhRCMiZ2qEELeN0Whk++XtvHvkXbJKsmiZYkPHYh/uevQVgv1DlI4nhKiH5EyNEKJeUqlU3OF/B7+O+pVJrSbSKd4Bq+QS3ljxFJ+e/pTyqnKlIwohGjApNUKI287K1Irnwmcw7KVXyAuy4IxvHh+e+JAx68ewJ0mGpIQQ/4wMPwkhFGU0GtmctJkFRxeQXZLN4ENu2Ds4c/djs2ju00bpeEIIhcnwkxCiwVCpVAwNHMr6Ueu533kUrvnmmCYXMXXLFD4/87kMSQkhbpicqRFC1CtHTu/imz2fsNPmHAB+dn480/wRIjoMVziZEEIJMk9NLaTUCNFwGI1GNiZtZOHRhehzChi1xxO9rw0PvDAPP6cApeMJIW4jGX4SQjRoKpWK4YHD+XXUr4y2GIDKCPnaHO7cNI6PT32MvkqvdEQhRD0kZ2qEEPXe8XN7+OT8Z+wvOQGAr5UPjziOJ2rgRIWTCSFuNRl+qoWUGiEaNqPRyJbkLSw4sgC3M2WExTejuJUNU2fMx8fWR+l4QohbRIafhBCNjkqlYkjAEH4d/Sthzh0xqIycME9m1NpRLDu5jLLKMqUjCiEUJmdqhBAN0pmEwyxJ/JSDGYcAaFPmxXifMYwaPA21Wv5/TYjGQs7UCCEavfYtuvLpHZ+xsO9CPCzdaXHEQNKK9bz84QNcLrqsdDwhhAKk1AghGqw/1pL6adhqnNoFU2JexQ6rM4xeN5olx5dQUlGidEQhxG0kw09CiEYjMTue+acWsS91HwC9kr0YGDaCMUMfkSEpIRooGX4SQjRJzV2C+GjgRyzuv5hWei+an9dw+etNPL1qKpcKLykdTwhxi0mpEUI0KiqVioG+A/ni3h8w6dGCBN8SYsqPcuf6O1l0bBG6cp3SEYUQt4gMPwkhGrUrhVd49+i77L66G9MKFcOOetNuyBDGRT0uQ1JCNAAy/CSEEL/ztfdl2cBlLB2wlG5p3jjkqzm3bgMPbXmQxIJEpeMJIeqQlBohRJPQ16cv859fhWnvII63L+ZQ9mHGrh/L/MPvkV+Uo3Q8IUQdkOEnIUSTc7X4Ku8deY9dKbvwzbCkx3lnAkbfwd0jn0alUikdTwjxP2T4SQgh/oK3rTdLBizhPwP/Q2iqCxZlajYf+YnJWycTnx+vdDwhxD8kZ2qEEE1aaZmOL755k29MotGpytCoNNztcyfTOj6Ek4Ob0vGEaPJkle5aSKkRQvyVdG0684/OZ/vl7fQ/5oJbgQWBdw9j/OBHUavkpLYQSpHhJyGEuEkeNh4s6reIZT0X41JihVm5iv8kfsHEzROJzYtVOp4Q4gbImRohhPh/yvQlLN+6hC+Lf6G0shS1Ss0Ei0imRjyDq6On0vGEaFLkTI0QQvwLFuZWPDLiJdaPWs9g/8GYl6pQrT3L509P48cDyzEYDUpHFELUQkqNEEL8CXdrd+b3nc+8zm9Qbq2myKqCN+MWcv+m+zmXe07peEKI/0eGn4QQ4gaUlZfy3fEVfJL4FSWVJaiNKiYU9mDKva/i5uytdDwhGq0GP/w0evRomjVrxtixY697bcOGDQQHB9OyZUs+//xzBdIJIZoiCzNLpnR7hF9H/8rQgKG0vGKN2f6rfPzCg/wYu4oqQ5XSEYVo8urlmZqYmBiKi4tZsWIFP/30U832yspK2rRpw65du7C3tycsLIz9+/fj5OT0t8eUMzVCiLoUve9n9q1czhnPXGL9i2nr1JaXw1+mg0sHpaMJ0ag0+DM1/fr1w9bW9rrthw8fpm3btnh5eWFjY8OQIUPYtm2bAgmFEE3dwJ538tLSnxhx5yPYmNpwLvccT/wwiblvTiIt+7LS8YRokm661OzevZuoqCg8PT1RqVSsXbv2un2WLVuGv78/FhYWhIeHc/jw4brISlpaGl5eXjXPvby8SE1NrZNjCyHEzTIzNef+dg/w6+hfGREYRfg5R8zO5jD3vWl8H/s9lYZKpSMK0aTcdKnR6XSEhISwbNmyWl9ftWoV06dPZ/bs2Rw/fpyQkBAiIyPJysqq2Sc0NJR27dpd90hLS/vnn0QIIRTibOnM273fYeikJyhyUXO4eTbvHHqHCRsmcCzjmNLxhGgyTG72B4YMGcKQIUP+9PVFixbx4IMPMnnyZAA+/vhjNm7cyJdffslLL70EwMmTJ/9RWE9Pz2vOzKSmptK1a9da99Xr9ej1+prnRUVF/+g9hRDiRvXvNoo+XaNoHb+aJSeWEJcfx5IPniHI1I/7HpmNn2dLpSMK0ajV6TU15eXlHDt2jIiIiP++gVpNREQEBw4c+NfH79q1K2fPniU1NRWtVsvmzZuJjIysdd+5c+dib29f8/Dx8fnX7y+EEH9Ho9YwodUENozewFivEbRNtsMiroBnfpjKinMrqDBUKB1RiEarTktNTk4OVVVVuLldu7Ktm5sbGRkZN3yciIgIxo0bx6ZNm/D29q4pRCYmJixcuJD+/fsTGhrKc88996d3Ps2cOZPCwsKaR0pKyj//YEIIcZMcLRyZHfE2PWc8QUZbMxIc81lwdAFj149lT+IupeMJ0Sjd9PDT7bBjx44/fW3EiBGMGDHib49hbm6Oubl5XcYSQoib1jNsCN07RdIjYR2Ljy/mct4ldrw1jxiPT7jvyTcI8AhSOqIQjUadnqlxdnZGo9GQmZl5zfbMzEzc3d3r8q2EEKLBUKvUjG45mvWj1jPBagg2JSYYrxZwz9b7+PT0p+ir9H9/ECHE36rTUmNmZkZYWBjR0dE12wwGA9HR0XTv3r0u30oIIRoce3N7XprwHv1emUFmX0e0qlI+PPEho9eNZvOBH5WOJ0SDd9PDT1qtloSEhJrnSUlJnDx5EkdHR3x9fZk+fToTJ06kc+fOdO3alcWLF6PT6WruhhJCiKauS4f+dG7fj01Jm1h4dCEVSdmcX7WSI2tWM+nlBfg6+CkdUYgG6aaXSYiJiaF///7XbZ84cSLLly8HYOnSpcyfP5+MjAxCQ0NZsmQJ4eHhdRL4n5JlEoQQ9ZGuQsfHX76KYVc8cb7FnGivY1K7SUxrPw1LE0ul4wmhuJv5/q6Xaz/dClJqhBD12anYA3x08Qv25R0CwE/jwSSH0YwZ+jBqdb1c0UaI20JKTS2k1Agh6juj0cjOKzt598i7BBwoJ+iqDXltrZj61DwCHQKVjieEIhr8gpZCCNEUqVQqBvoNZO3ItbT2C6FSbeSAXSJ3rr+ThUcXoqvQKR1RiHpNztQIIUQ9lZB6gQ9i/0PM1RgAOuS6M8xvCBNGPyNDUqLJkDM1QgjRCLTwas2HAz9k2cBlBJj70PqECek/7uTZ/9xDXF6c0vGEqHek1AghRD3Xx7sPP4z6EYce7ci1r2CX9Tnu2nAX7xx6h0J9odLxhKg3ZPhJCCEakLSiVBYeX8S2y9vACBGnPQnt0p/77pqBRlMvV74R4l+R4SchhGikPO28WNhvIZ/f8TldtAF4p5qSsX43k3+6l7M5Z5WOJ4SipNQIIUQDFO4Rzn8eWoX1oBDOtSrlRNl57tl4D3P2zyFbm6V0PCEUIcNPQgjRwGWXZPP+sff59dKv2OpMGHrQA/eBXZl478uYmpgpHU+If0WGn4QQoglxsXLhnd7vsGLwCrpn+mKpV3Ph0B7u3ng3xzOPKx1PiNtGztQIIUQjUlFZzorv5vJT2Q5SzQsAGO43lEeCpuHn2VLZcEL8A7JMQi2k1AghmpL8snyWnFjCz/E/0zbRlpBEe5oN7sKUe1/FVG2qdDwhbpgMPwkhRBPXzKIZs7vP5ruh3xFU6IxppZrNV7Yydv1YDqQdUDqeELeEnKkRQohGrqqqku83fMhnpWvJ0+cBMMSmDw91fpQWfu0UTifEX5Php1pIqRFCNHWF+kL+c/I/rLrwA8P2umGnM8V2XDemjnwBc4250vGEqJUMPwkhhLiOvbk9M8NnsrLfF5hZWVGlMfBl7k+MXjea3Vd3Kx1PiH9NztQIIUQTZDAYWHf0B5Ymf0FWafVkfcOLOvHAkKdp3aKTwumE+C8ZfqqFlBohhLierkLHJ6c/YeP+Hxiy1wWDGswf7suDPZ/AytRK6XhCyPCTEEKIG2Ntas30sOm8H7mEEk8LLruV8Pnlrxm5biTbkrfRRP6/VzQScqZGCCEEUD0kFX1pOwtOLiJNl4ZphYrhCUGMfuBZOrbpqXQ80UTJ8FMtpNQIIcSNKa0s5auzX3Hsxx9pk2hDgU0Ftg8P4NHQx7A1s1U6nmhiZPhJCCHEP2ZpYsljoY8xY9r7lPhbcaRVHl9f+IaoNVGsvbiWKkOV0hGFqJWcqRFCCPGX9qbu5d3D75JclIxPpiXhVzy4Y+rjdO94h9LRRBMgw0+1kFIjhBD/XEVVBSvPryT5g9U4FJtwunkhzaMG8VTHp3CwcFA6nmjEpNTUQkqNEEL8e0lXY/lm+Tv84nGSShMjdmZ2PNHqEca2H4+piZnS8UQjJKWmFlJqhBCi7hzNOMrcw3OJz4+n33FnXPQ29HvwEfqEDVM6mmhk5EJhIYQQt1Rn986sGr6KF1o/g2euJZYFBl7d/xov73mZ7JJspeOJJkrO1AghhPhX0rKS+XLjIn5Ux2DEiLWpNQ+6TODevg9jYWapdDzRwMmZGiGEELeNp6s/r05ewnfDvqO9c3sM2lKylm/jvcfHsev8FqXjiSZESo0QQog60c65Hd8M/YZn/B/GqIFSlZ6nDj/P9JjppGnTlI4nmgAZfhJCCFHnsvPT+eLwx/yQuY4qYxWWagvuqxzIlPEvY2MlfweLGyfDT0IIIRTl0syDlyJf58eoH+ns1hnvZBMqtpxl/rN3szN5pyyUKW4JKTVCCCFumaBmQXwZ+SUTOt5PqaWRWI98nv7taR6NfpTkwmSl44lGpl6WmtGjR9OsWTPGjh17zfaUlBT69etHmzZt6NChA6tXr1YooRBCiBulUqm4a/hjPPnh13QbOhYTtQn7Uvcx+ZtxLFz0OIXFeUpHFI1EvbymJiYmhuLiYlasWMFPP/1Usz09PZ3MzExCQ0PJyMggLCyM+Ph4rK2t//aYck2NEELUD8mFybx7eB42vyTgnmfBlYBKIh96iiEBQ1CpVErHE/VMg7+mpl+/ftjaXr+8vYeHB6GhoQC4u7vj7OxMXp40fCGEaEj87f1ZNvA/9Bx5N1o7I4f8M3hxz4tM3jqZuLw4peOJBuymS83u3buJiorC09MTlUrF2rVrr9tn2bJl+Pv7Y2FhQXh4OIcPH66LrNc4duwYVVVV+Pj41PmxhRBC3FpqtZpRg6fy4kc/M7n7o1hoLDiWeYy3F07jvXkPkp2frnRE0QDddKnR6XSEhISwbNmyWl9ftWoV06dPZ/bs2Rw/fpyQkBAiIyPJysqq2Sc0NJR27dpd90hLu7F5DPLy8njggQf49NNPbza+EEKIesTCxIKHQx5m3ah1RLr0p90lOzQn0nl8+b38cvEXDEaD0hFFA/KvrqlRqVSsWbOGUaNG1WwLDw+nS5cuLF26FACDwYCPjw9PPvkkL7300g0fOyYmhqVLl15zTQ2AXq9n0KBBPPjgg9x///1/+vN6vR69Xl/zvKioCB8fH7mmRggh6rEtMd+xbcf3bG+RDCpo59SOFzvOINQrTOloQiGKXVNTXl7OsWPHiIiI+O8bqNVERERw4MCBf318o9HIpEmTGDBgwF8WGoC5c+dib29f85BhKiGEqP8G97uHd9/4hRldZmBtas35rLP88trLzH1jImk5V5SOJ+q5Oi01OTk5VFVV4ebmds12Nzc3MjIybvg4ERERjBs3jk2bNuHt7V1TiPbt28eqVatYu3YtoaGhhIaGcubMmVqPMXPmTAoLC2seKSkp//yDCSGEuG1M1aZMbDuRX0f9ymiTvtjrTKlMyGTCxgl8e+FbKg2VSkcU9ZSJ0gFqs2PHjlq39+rVC4PhxsZXzc3NMTc3r8tYQgghbiMXKxfmPLCUmKB1fHtyJfmkMu/wPH6++DNP+06jb+hQpSOKeqZOz9Q4Ozuj0WjIzMy8ZntmZibu7u51+VZCCCGaiH7dRvLxQz/yWrfXsDe3pyjhCkfn/oe3Xr2X9GJZKFP8V52WGjMzM8LCwoiOjq7ZZjAYiI6Opnv37nX5VkIIIZoQjVrDXcF3sWHUBgaadMGIkcv6q4xcP4rPz3xOeVW50hFFPXDTw09arZaEhISa50lJSZw8eRJHR0d8fX2ZPn06EydOpHPnznTt2pXFixej0+mYPHlynQYXQgjR9DhYODDzqY840HsbZ5NWUFp8mg+Of8DGs2uZ6nIXwyMeUDqiUNBN39IdExND//79r9s+ceJEli9fDsDSpUuZP38+GRkZhIaGsmTJEsLDw+sk8D91o7eEVVVVUVFRcRuTifrI1NQUjUajdAwhxF8wGo1suLSBhUcXEnTESHCKLYUh9jz45Dx8bOWO18biZm7prpdrP90Kf/dLMRqNZGRkUFBQcPvDiXrJwcEBd3d3WYtGiHquWF/MJ0uex3gshS3hmRQ4G5nUbhLT2k/D0sRS6XjiX5JSU4u/+6Wkp6dTUFCAq6srVlZW8kXWhBmNRkpKSsjKysLBwQEPDw+lIwkhbsDZ5ON8EP8RB9MPAhCS78nIFiO5c+gjqNX1cqlDcQOk1NTir34pVVVVxMfH4+rqipOTk0IJRX2Tm5tLVlYWQUFBMhQlRANhNBrZcWUH7++fT49NGizLNaT2tefxe96guUNzpeOJf6DBr9J9u/1xDY2VlZXCSUR98sefB7nGSoiGQ6VSMchvED+M+BGbLkHk21UQbXmaO9ffybuH36W4vFjpiOIWqpeT7ylFhpzE/5I/D0I0XHbWDjzz1AdcKbhM8fFF7EzZyTfnvyHjx510Co/gnjHPotHIV2BjI2dqhBBCNFq+Dn58MOADPo74mDCtH94ppqT/EsPUn+/nXM45peOJOialRjRoycnJqFQqTp48qXQUIUQ91tOrJx9N+xHLAe04F6zjWOlZ7t54N3P2zyFHm6V0PFFHpNQ0cJMmTUKlUqFSqTA1NSUgIIAXXniBsrKyOjv+qFGj6uRYQgihJEsLKx57eB7zX/iR4YHDMWJk66l1fPT4JD5b+ToVlTIrcUMnpaYRGDx4MOnp6Vy6dIn333+fTz75hNmzZysdq8EoL5e/yIRoSlytXJnbey4rBq+gZ4YvVmVq4g7uYfzG8RzNOKp0PPEvSKn5E0ajkZLySkUeN3uXvbm5Oe7u7vj4+DBq1CgiIiLYvn07UL321ty5cwkICMDS0pKQkBB++umna37+3LlzDB8+HDs7O2xtbenduzeJiYnMmTOHFStWsG7dupqzQTExMQC8+OKLBAUFYWVlRWBgIK+99to1dwnNmTOH0NBQvv76a/z9/bG3t2fChAkUF//3zoPi4mLuvfderK2t8fDw4P3336dfv34888wzNfuoVCrWrl17TV4HB4ea2av/v6qqKqZOnVrzeYODg/nggw+u2eePs09vv/02np6eBAcH39TvWwjROHRy68TcV1fjMLQr50LLuViQwOStk3lh5/MkpcYqHU/8A3Lp958oraiizaytirz3+TcisTL7Z/9qzp49y/79+/Hz8wNg7ty5fPPNN3z88ce0bNmS3bt3c9999+Hi4kLfvn1JTU2lT58+9OvXj507d2JnZ8e+ffuorKxkxowZXLhwgaKiIr766isAHB0dAbC1tWX58uV4enpy5swZHnzwQWxtbXnhhRdqsiQmJrJ27Vo2bNhAfn4+d911F/PmzePtt98GYPr06ezbt4/169fj5ubGrFmzOH78OKGhof/4d2cwGPD29mb16tU4OTmxf/9+HnroITw8PLjrrrtq9ouOjsbOzq6m/AkhmiZTUzOmTpzFnWVP8eGJD1kdv5orv+1j1efnsIvsxLT7ZmGmMVM6prhBUmoagQ0bNmBjY0NlZSV6vR61Ws3SpUvR6/W888477Nixo2aV9MDAQPbu3csnn3xC3759WbZsGfb29vzwww+YmpoCEBQUVHNsS0tL9Ho97u7u17znq6++WvPP/v7+zJgxgx9++OGaUmMwGFi+fDm2trYA3H///URHR/P2229TXFzMihUr+O677xg4cCAAX331FZ6env/qd2Fqasrrr79e8zwgIIADBw7w448/XlNqrK2t+fzzzzEzk7+shBDVC2W+1v01xrQcw49vvYJpVSXbU3awaf1JXujyAn28+ygdUdwAKTV/wtJUw/k3IhV775vRv39/PvroI3Q6He+//z4mJibceeednDt3jpKSEgYNGnTN/uXl5XTs2BGAkydP0rt375pCc6NWrVrFkiVLSExMRKvVUllZed1Mj/7+/jWFBsDDw4OsrOq7DC5dukRFRQVdu3ated3e3r5OhoKWLVvGl19+yZUrVygtLaW8vPy6sz/t27eXQiOEuE5b57bMXvgzP2xcRkHpGnKKLvN49ONEWHXnoU6P0rp5R6Ujir8gpeZPqFSqfzwEdLtZW1vTokULAL788ktCQkL44osvaNeuHQAbN27Ey8vrmp8xNzcHqs/E3KwDBw5w77338vrrrxMZGVlzpmfhwoXX7Pf/i5JKpcJgMNzUe6lUquuuMfqrGX5/+OEHZsyYwcKFC+nevTu2trbMnz+fQ4cOXbOftbX1TeUQQjQdao2Ge0Y8xYjyKXxy+hO+PfcNltuSWP/Lq2wd1ZGHxryMlanMQF8fNYxvbXHD1Go1L7/8MtOnTyc+Ph5zc3OuXLlC3759a92/Q4cOrFixgoqKilrP1piZmVFVVXXNtj+u2XnllVdqtl2+fPmmcgYGBmJqasqRI0fw9fUFoLCwkPj4ePr0+e9pXhcXF9LT02ueX7x4kZKSkj897r59++jRowePPfZYzbbExMSbyiaEEAA2ZjY81/k5hroPYvWJORhKyvi28Fc2rD3IjC4ziPSLlJnH6xm5+6kRGjduHBqNhk8++YQZM2bw7LPPsmLFChITEzl+/DgffvghK1asAOCJJ56gqKiICRMmcPToUS5evMjXX39NXFwcUD2EdPr0aeLi4sjJyaGiooKWLVty5coVfvjhBxITE1myZAlr1qy5qYy2trZMnDiR559/nl27dnHu3DmmTp2KWq2+5i+JAQMGsHTpUk6cOMHRo0d55JFH/nKorGXLlhw9epStW7cSHx/Pa6+9xpEjR/7Bb1EIIaq19u7Aq4t+IuS5qTg5eZBZksnzvz3P80vu4fiFvUrHE/9DSk0jZGJiwhNPPMF7773HzJkzee2115g7dy6tW7dm8ODBbNy4kYCAAACcnJzYuXMnWq2Wvn37EhYWxmeffVZTHB588EGCg4Pp3LkzLi4u7Nu3jxEjRvDss8/yxBNPEBoayv79+3nttdduOueiRYvo3r07w4cPJyIigp49e9K6dWssLCxq9lm4cCE+Pj707t2be+65hxkzZvzlwqMPP/wwY8aMYfz48YSHh5Obm3vNWRshhPgn1Go1QzvdydqRa3ks5DHctNa47S9ix+tzeXf76xSVFykdUQAq481OitJA/dXS5WVlZSQlJREQEHDNF6q4vXQ6HV5eXixcuJCpU6cqHUf+XAgh/lRc0ilWf/QOGeVZ/NYxB0cLR57p9AwjW4xErZLzBXXpr76//z/5zQvFnDhxgu+//75mWOzee+8FYOTIkQonE0KIvxYcEMKr761iwrOvE2AfQF5ZHm/sns2sF8aw74Qyc5wJKTVCYQsWLCAkJISIiAh0Oh179uzB2dlZ6VhCCHFDevn14eeon5nReQZhSc40u1LJ1iWLmLXnNXJLc5WO1+TI8BMyzCBqJ38uhBA3I/lqHN99/CZ77OK56laKraktj3d8nLtajsPURObF+qdk+EkIIYS4zfy9g3n5rW+YO+lTWju2priimG9+/YB3nhjDjgO/KB2vSZBSI4QQQtShUNdQvh/2Pa+Gv0pYoiM2+fDN+veZ8dsMMnQZSsdr1KTUCCGEEHVMo9YwvtV4nnjrUyo7uXOmZTFbk7cyYu0IPjm0lBK9TumIjZKUGiGEEOIWcXfx5cUXP+fbUT/QybUTpZWlnPn+ZxY8cRcb936vdLxGR0qNEEIIcYu1cmzF8sHLebPjLDzzrLAsgndPzOexHY9xuejmlpkRf05KjRBCCHEbqFQqRnUYxyNLvkI9sj3FDkb2pO5h9LrRLN7wBkW6AqUjNnhSakS9NGnSJEaNGvWX+/Tr149nnnnmho8ZExODSqWioKDgX2UTQoh/w9HBlefumccvI36hp1dPNKVVlH5/kA+euId1x1bRRGZauSWk1DRgKpXqLx9z5sxROuI/9sEHH7B8+XKlYwghxC0TYB/ARwM/YnabFzGYqigxqeC1M28xacskYvNilY7XIJkoHUD8c+np6TX/vGrVKmbNmlWzujaAjY1NzT8bjUaqqqowMWkY/8rt7e2VjiCEELecSqVieK976BM6hG+Ofon51R84nnWc8b+OZ7yuJ9PufhXXZp5Kx2ww5EzNnzEaoVynzOMGTz26u7vXPOzt7VGpVDXPY2NjsbW1ZfPmzYSFhWFubs7evXtJTExk5MiRuLm5YWNjQ5cuXdixY8c1x/X39+fNN9/k7rvvxtraGi8vL5YtW3bNPgUFBUybNg0XFxfs7OwYMGAAp06duuYYtZ09+sOZM2cYMGAAlpaWODk58dBDD6HVamte///DTzqdjgceeAAbGxs8PDxYuHDhdb+Pr7/+ms6dO2Nra4u7uzv33HMPWVlZN/S7FEIIJdnZNOOxfs/x6+hfGeI/hMAUS8x/u8J/Zkzlm3NfU2moVDpig9Aw/rddCRUl8I5C7fjlNDCzrpNDvfTSSyxYsIDAwECaNWtGSkoKQ4cO5e2338bc3JyVK1cSFRVFXFwcvr6+NT83f/58Xn75ZV5//XW2bt3K008/TVBQEIMGDQJg3LhxWFpasnnzZuzt7fnkk08YOHAg8fHxODo6cuTIEaqqqgCoqqpi7NixmJqaAtUFJTIyku7du3PkyBGysrKYNm0aTzzxxJ8OOT3//PP89ttvrFu3DldXV15++WWOHz9OaGhozT4VFRW8+eabBAcHk5WVxfTp05k0aRKbNm2qk9+lEELcau7W7rzX9z22GX/gSNK3xHoX8vPR9/g54Rde7Poi3Ty6KR2xXpNS08i98cYbNUUEwNHRkZCQkJrnb775JmvWrGH9+vU88cQTNdt79uzJSy+9BEBQUBD79u3j/fffZ9CgQezdu5fDhw+TlZWFubk5UL0w5dq1a/npp5946KGHcHFxqTnW008/TXp6OkeOHAHgu+++o6ysjJUrV2JtXV3eli5dSlRUFO+++y5ubm7XfAatVssXX3zBN998w8CBAwFYsWIF3t7e1+w3ZcqUmn8ODAxkyZIldOnSBa1We81QnBBC1Hd39JtA324jWXtpHamnl5FQkMDzPz9GRF5r7ntwFs09gpWOWC9JqfkzplbVZ0yUeu860rlz52uea7Va5syZw8aNG0lPT6eyspLS0lKuXLlyzX7du3e/7vnixYsBOHXqFFqtFicnp2v2KS0tJTEx8Zptn376KV988QX79++vKToXLlwgJCSkptBAdYkyGAzExcVdV2oSExMpLy8nPDy8ZpujoyPBwdf+R33s2DHmzJnDqVOnyM/Px2AwAHDlyhXatGnzl78nIYSob8wtLBnfZgKDmw/hPyeWkbtyJza5hbw//1GC7x/N1HZTsarD74vGoF6WmtGjRxMTE8PAgQP56aefrnu9pKSE1q1bM27cOBYsWHBrQqhUdTYEpKT/LQ4AM2bMYPv27SxYsIAWLVpgaWnJ2LFjKS8vv+FjarVaPDw8iImJue41BweHmn/etWsXTz75JN9//z0dOnT4px/hhvwxpBUZGcm3336Li4sLV65cITIy8qY+mxBC1Df25vbM7PYyBwgl+uvPONYyn99Of8q6hHVMD5vOkIAh11yz2JTVywuFn376aVauXPmnr7/99tt06ybjiv/Evn37mDRpEqNHj6Z9+/a4u7uTnJx83X4HDx687nnr1q0B6NSpExkZGZiYmNCiRYtrHs7OzgAkJCQwduxYXn75ZcaMGXPNsVq3bs2pU6fQ6f679sm+fftQq9XXnX0BaN68Oaamphw6dKhmW35+PvHx8TXPY2Njyc3NZd68efTu3ZtWrVrJRcJCiEale7ehvPzhz7w5dD5eNl5klmTy/Rfv8OorYzl56dDfH6AJqJelpl+/ftja2tb62sWLF4mNjWXIkCG3OVXj0LJlS3755RdOnjzJqVOnuOeee2qGaf7Xvn37eO+994iPj2fZsmWsXr2ap59+GoCIiAi6d+/OqFGj2LZtG8nJyezfv59XXnmFo0ePUlpaSlRUFB07duShhx4iIyOj5gFw7733YmFhwcSJEzl79mzNGZ3777//uqEnqL41ferUqTz//PPs3LmTs2fPMmnSJNTq//7x9fX1xczMjA8//JBLly6xfv163nzzzVv0WxRCCGWo1Woi/CJYO3ItjwU9SLtL9jgm6nl19VO8ceAN8svylY6oqJsuNbt37yYqKgpPT09UKhVr1669bp9ly5bh7++PhYUF4eHhHD58uC6yAtXDJ3Pnzq2z4zU1ixYtolmzZvTo0YOoqCgiIyPp1KnTdfs999xzHD16lI4dO/LWW2+xaNEiIiMjgep5FTZt2kSfPn2YPHkyQUFBTJgwgcuXL+Pm5kZmZiaxsbFER0fj6emJh4dHzQPAysqKrVu3kpeXR5cuXRg7diwDBw5k6dKlf5p7/vz59O7dm6ioKCIiIujVqxdhYWE1r7u4uLB8+XJWr15NmzZtmDdv3q0bmhRCCIVZmFjwaPenGP7Kq+hCHLnsVsLq+NUMWzOMr49/1WRvAVcZb3I+5s2bN7Nv3z7CwsIYM2YMa9asuWY+kVWrVvHAAw/w8ccfEx4ezuLFi1m9ejVxcXG4uroCEBoaSmXl9b/wbdu24elZfRt1TEwMS5cuveaamnXr1rF3717mz5/P8uXLOXv27A1/cRUVFWFvb09hYSF2dnbXvFZWVkZSUhIBAQFYWFjczK+jUfL39+eZZ565qSUIGiP5cyGEaCiOZhxl3uF5XMyNY+QeTyoczRjy0NP0CYpQOtq/9lff3//fTV8oPGTIkL8c+lm0aBEPPvggkydPBuDjjz9m48aNfPnllzW3CJ88efJm3xaovq7jhx9+YPXq1Wi1WioqKrCzs2PWrFnX7avX69Hr9TXPi4qK/tF7CiGEEPVdZ/fOrBq+im+3LCWzZBv6ikqe2T2dvqkDea7zc3jbev/9QRqBOr2mpry8nGPHjhER8d9mqFariYiI4MCBA//6+HPnziUlJYXk5GQWLFjAgw8+WGuh+WNfe3v7moePj8+/fn8hhBCivtKoNTww9GnGvv0uJsPbYzBTs+PKDkauHckHm95CV677+4M0cHVaanJycqiqqrruYk83N7eai0RvREREBOPGjWPTpk14e3v/o0I0c+ZMCgsLax4pKSk3fYymKjk5uckPPQkhREMV0LwdL41/j9VRqwn3CMchR0XlioO88fxYNiVsbNSrgNfLeWr+/1pEtZk0adJfvm5ubl4z260QQgjR1LRs1pLPBn3Gj/mLuaLeQb5pCS/ue4lVF3/kpa4v0dqptdIR61ydnqlxdnZGo9GQmZl5zfbMzEzc3d3r8q2EEEII8TdUKhXjxz/L/fOX0unOO7E0seR41nEe+OVu3v7qKfJK85SOWKfqtNSYmZkRFhZGdHR0zTaDwUB0dPR10+4LIYQQ4vZw9w7gke5PsX7Ueob4DyE03h6LLZeYNfsuvr3wLRWGCqUj1ombLjVarZaTJ0/W3MGUlJTEyZMna9YOmj59Op999hkrVqzgwoULPProo+h0upq7oYQQQgihjD9WAR8aMpoqE4h1L2De4XmMWz+OA2n//oYepd30NTVHjx6lf//+Nc+nT58OwMSJE1m+fDnjx48nOzubWbNmkZGRQWhoKFu2bKl1plghhBBC3H6j73sG3YhJ+Kdv48MTH5JYmMjbK58hxKYND9//Or72vkpH/EduevK9hkom3xM3S/5cCCGagkJ9If85vITKz/ZjpddwqH0BvYeNrzergN/M5Hv1cu0nUbf+bDkLIYQQwt7cnhd7vEzXkWMpdTYlzrOQT09/yoi1I9ictLlB3QIupaYRyM7O5tFHH8XX1xdzc3Pc3d2JjIxk3759AKSnp8sCoEIIIf6UWqMh8s6pvPzhzywauLh6FXBdJpsWv8eMBeM5l3lG6Yg3pF7OUyNuzp133kl5eTkrVqwgMDCQzMxMoqOjyf2/9u4/psk7jwP4uy0/a6llMPlxFE0cY3aWloBF12F0oGxHZDJhEZMJ6HZzOt2iM95wruLd4slhhjNm2dw5XTIzJzfJ7jJztzUYmA4UNnDZPx4GHDuRqgRaQAHpc384unEOV1zp0z59v5Imtn2+fd5tPsKHb5/n+V6/DgA8nZ6IiNwil8uRNTML5t+Z8bd/VuLmlQaM2gbwh2nFyEldjo2pGxEZFil2zAlxpmYCgiBgcGRQlNtkpvp6e3tRX1+PPXv2YPHixZg5cyZMJhNeffVV5OXlARj/9VNHRwdkMhk+/vhjZGZmIjw8HPPmzcOFCxdw7tw5pKenQ6VS4YknnsDVq1dd+ykpKcHy5ctRXl6O+++/H2q1GuvWrcPw8LBrm+rqauj1eoSHhyMqKgrZ2dkYGLh9WW6n04ldu3YhISEBoaGhrgPIx4zl+uSTT7B48WIolUoYDAaPLK9BRESTExYUhvV5ZZhfWoLBjBjYp438tAp46xGfPQWcMzUTuHHrBjKOZoiy78ZVjW4fnKVSqaBSqVBTU4P58+e7fRVli8WCqqoqJCYmYs2aNVi1ahUiIiKwb98+KJVKPP3003j99dfx9ttvu8ZYrVaEhYXh1KlT6OjoQGlpKaKiovDGG2+gq6sLRUVFqKioQH5+PhwOB+rr610N2r59+7B371688847SE1NxaFDh5CXl4fvvvsOSUlJrn1s374dlZWVSEpKwvbt21FUVIS2tjYEBbFUiYi8SSaXw/x4AcyPF2BZdzP+cvYv6Pzvf/D9Xz/CluS/Y2XJNjyiNYsdcxzO1Pi5oKAgHD58GEeOHIFGo4HZbEZZWRnOnz9/13GvvPIKcnJyMGfOHLz00ktobm7Gjh07YDabkZqairVr16K2tnbcmJCQEBw6dAgPP/wwcnNzsWvXLrz11ltwOp3o6urCrVu38NRTT2HWrFnQ6/VYv349VCoVAKCyshLbtm3DypUrkZycjD179sBoNKKqquqOXLm5uXjwwQdRXl6OS5cuoa2tzaOfGRERTU5aTBo+yv0IzyIXYSMKhF4exPPWdXi59mX84PhB7Hgu/PN3AuFB4Whc1SjavidjxYoVyM3NRX19PRoaGnDy5ElUVFTgvffem3CNrJSUFNe/x64hpNfrxz1ms9nGjTEYDFAqf5pBWrBgAfr7+9HZ2QmDwYCsrCzo9Xrk5ORg6dKlKCgoQGRkJOx2Oy5fvgyzeXxHbzab0draOmGuuLg4AIDNZsNDDz00iU+EiIg8TSFXYO3zf0aL7t/4V/+XUFz7B6zfW/Hl9/VYnVCIZxduEv0UcM7UTEAmk0EZrBTlJpPJJp03LCwMS5YswY4dO3DmzBmUlJTAYrFMuH1wcPC49/pLjzmdTrf3r1Ao8Pnnn+PkyZPQ6XTYv38/kpOT0d7ePqn38Uu5JpODiIimjkwmQ+rCHPzx939C9bJqZMRlYHZ7KIYO1uPl3fminwLOpkaidDqd6yBdT2ltbcWNGzdc9xsaGqBSqaDVagHcLnaz2Yzy8nJ88803CAkJwYkTJ6BWqxEfH+86xXzM6dOnodPpPJqRiIi844HIB3BwyUHkhD8CuSDDNZkdO8/sRM9N8RbJ5NdPfu769esoLCzEmjVrkJKSgoiICDQ1NaGiogJPPvmkR/c1PDyMtWvX4rXXXkNHRwcsFgtefPFFyOVyNDY2wmq1YunSpZgxYwYaGxtx9epVzJlze2n7rVu3wmKxYPbs2TAajXj//ffR0tKCDz/80KMZiYjIe2QyGdZsq8TF801QOZuhClUhKjxKtDxsavycSqVCRkYG3nzzTVy8eBEjIyPQarV47rnnUFZW5tF9ZWVlISkpCQsXLsTQ0BCKioqwc+dOAIBarUZdXR2qqqpgt9sxc+ZM7N2713XRv02bNqGvrw9btmyBzWaDTqfDp59+Ou7MJyIi8k+zU9IxG+lix+DaTwDX+HFHSUkJent7A2q5BdYFEZH4uPYTERERBRw2NURERCQJPKaG3HL48GGxIxAREd0VZ2qIiIhIEtjU/EyAHDNNbmI9EBH5FzY1+OkqtoODgyInIV8yVg8/v8oxERH5Lh5Tg9uX+NdoNK61jpTKe1uqgKRBEAQMDg7CZrNBo9FAoVCIHYmIiNzApuZHsbGxAHDHIo4UuDQajasuiIjI97Gp+ZFMJkNcXBxmzJiBkZERseOQyIKDgzlDQ0TkZ9jU/B+FQsFfZkRERH6IBwoTERGRJLCpISIiIklgU0NERESSEDDH1IxdSM1ut4uchIiIiNw19nvbnQuiBkxT43A4AABarVbkJERERDRZDocD06dPv+s2MiFArgXvdDpx+fJlRERETMmF9ebNm4dz5875xOvdy1h3x3hqu4met9vt0Gq16OzshFqt/tX9+CJP14K39+ertefutqw91t5UjGHt3d1U1p4gCHA4HIiPj4dcfvejZgJmpkYulyMhIWHKXl+hUHi0GH/L693LWHfHeGq7X3terVb77X9uT9eCt/fnq7Xn7rasPdbeVIxh7d3dVNfer83QjOGBwh6yYcMGn3m9exnr7hhPbefpz8uXePu9BUrtubsta89/98fa81++8t4C5usn8g92ux3Tp09HX1+f3/7FQv6JtUdiYe15DmdqyKeEhobCYrEgNDRU7CgUYFh7JBbWnudwpoaIiIgkgTM1REREJAlsaoiIiEgS2NQQERGRJLCpISIiIklgU0N+JT8/H5GRkSgoKBA7CgWQzs5OLFq0CDqdDikpKTh+/LjYkShA9Pb2Ij09HUajEXPnzsXBgwfFjuTTePYT+ZVTp07B4XDgyJEjqK6uFjsOBYiuri50d3fDaDTiypUrSEtLw4ULFzBt2jSxo5HEjY6OYmhoCEqlEgMDA5g7dy6ampoQFRUldjSfxJka8iuLFi1CRESE2DEowMTFxcFoNAIAYmNjER0djZ6eHnFDUUBQKBRQKpUAgKGhIQiC4NZq1YGKTQ15TV1dHZYtW4b4+HjIZDLU1NTcsc2BAwcwa9YshIWFISMjA2fPnvV+UJIcT9Zec3MzRkdHodVqpzg1SYEnaq+3txcGgwEJCQnYunUroqOjvZTe/7CpIa8ZGBiAwWDAgQMHfvH5Y8eOYfPmzbBYLPj6669hMBiQk5MDm83m5aQkNZ6qvZ6eHqxevRrvvvuuN2KTBHii9jQaDVpbW9He3o6jR4+iu7vbW/H9j0AkAgDCiRMnxj1mMpmEDRs2uO6Pjo4K8fHxwu7du8dtV1tbK6xYscIbMUmC7rX2bt68KWRmZgoffPCBt6KSxPyWn3tjXnjhBeH48eNTGdOvcaaGfMLw8DCam5uRnZ3tekwulyM7OxtfffWViMlI6typPUEQUFJSgsceewzPPPOMWFFJYtypve7ubjgcDgBAX18f6urqkJycLEpef8CmhnzCtWvXMDo6ipiYmHGPx8TE4MqVK6772dnZKCwsxGeffYaEhAQ2PPSbuVN7p0+fxrFjx1BTUwOj0Qij0Yhvv/1WjLgkIe7U3qVLl5CZmQmDwYDMzExs3LgRer1ejLh+IUjsAEST8cUXX4gdgQLQo48+CqfTKXYMCkAmkwktLS1ix/AbnKkhnxAdHQ2FQnHHAXDd3d2IjY0VKRUFAtYeiYW153lsasgnhISEIC0tDVar1fWY0+mE1WrFggULRExGUsfaI7Gw9jyPXz+R1/T396Otrc11v729HS0tLbjvvvuQmJiIzZs3o7i4GOnp6TCZTKiqqsLAwABKS0tFTE1SwNojsbD2vEzs068ocNTW1goA7rgVFxe7ttm/f7+QmJgohISECCaTSWhoaBAvMEkGa4/EwtrzLq79RERERJLAY2qIiIhIEtjUEBERkSSwqSEiIiJJYFNDREREksCmhoiIiCSBTQ0RERFJApsaIiIikgQ2NURERCQJbGqIiIhIEtjUEBERkSSwqSEiIiJJYFNDREREksCmhoiIiCThfx0pArS0/H03AAAAAElFTkSuQmCC",
|
||
"text/plain": [
|
||
"<Figure size 640x480 with 1 Axes>"
|
||
]
|
||
},
|
||
"metadata": {},
|
||
"output_type": "display_data"
|
||
}
|
||
],
|
||
"source": [
|
||
"import numpy as np\n",
|
||
"import matplotlib.pyplot as plt\n",
|
||
"\n",
|
||
"def do_int(n):\n",
|
||
" # integration limits\n",
|
||
" x1 = 0.0\n",
|
||
" x2 = np.pi/2\n",
|
||
" x = np.linspace(x1, x2, n)\n",
|
||
" dx = (x2-x1)/(n-1)\n",
|
||
" \n",
|
||
" # Choose the function to integrate\n",
|
||
" if True:\n",
|
||
" f = np.cos(x)\n",
|
||
" I0 = 1.0\n",
|
||
" else: \n",
|
||
" # Polynomial\n",
|
||
" a = 4\n",
|
||
" I0 = (np.pi/2)**(a+1)/(a+1) # The analytic value of the integral\n",
|
||
" f = x**a # The function to integrate\n",
|
||
" \n",
|
||
" I1 = np.sum(f[:-1])*dx # rectangle rule\n",
|
||
" I2 = np.sum(f[1:-1]*dx) + 0.5*dx*(f[0]+f[-1]) # trapzoidal rule\n",
|
||
" I3 = np.sum(4*f[1:-1:2]*dx) + np.sum(2*f[2:-1:2]*dx) + f[0]*dx + f[-1]*dx #simpson's rule\n",
|
||
" I3 = I3/3.0\n",
|
||
" return abs((I1-I0)/I0), abs((I2-I0)/I0), abs((I3-I0)/I0)\n",
|
||
"\n",
|
||
"# try different numbers of points, generated as 2^n + 1 so the number is odd\n",
|
||
"n_vals = [2**(i+1) + 1 for i in range(11)]\n",
|
||
"err1=np.array([])\n",
|
||
"err2=np.array([])\n",
|
||
"err3=np.array([])\n",
|
||
"for n in n_vals:\n",
|
||
" e1, e2, e3 = do_int(n)\n",
|
||
" print(n, e1, e2, e3)\n",
|
||
" err1 = np.append(err1, e1)\n",
|
||
" err2 = np.append(err2, e2)\n",
|
||
" err3 = np.append(err3, e3)\n",
|
||
" \n",
|
||
"plt.plot(n_vals, err1, label='Rectangular')\n",
|
||
"plt.plot(n_vals, err2, label='Trapezoidal')\n",
|
||
"plt.plot(n_vals, err3, label='Simpson')\n",
|
||
"\n",
|
||
"n = np.array(n_vals)\n",
|
||
"dx = (np.pi/2)/(n-1)\n",
|
||
"# Plot the Euler Maclaurin error formulas\n",
|
||
"plt.plot(n_vals, dx/2, \":\")\n",
|
||
"plt.plot(n_vals, dx**2/12, \":\")\n",
|
||
"plt.plot(n_vals, dx**4/180, \":\")\n",
|
||
"\n",
|
||
"plt.yscale('log')\n",
|
||
"plt.xscale('log')\n",
|
||
"plt.legend()\n",
|
||
"plt.show()"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "e0bb8f63",
|
||
"metadata": {},
|
||
"source": [
|
||
"Notes:\n",
|
||
"- The scalings with the number of points are $1/N$ for rectangular, $1/N^2$ for trapezoidal, and $1/N^4$ for Simpson's rule. (You might have been expecting $1/N^3$ for Simpson's rule but because of the way it was constructed using double intervals, the third order term is antisymmetric over the double interval and cancels.)\n",
|
||
"- If you try integrating polynomials, you should find that the error goes to zero (machine precision) for cubic polynomials and below (Simpson), linear or below (trapezoid) or constant (rectangular).\n",
|
||
"- Some special cases can give surprising results. For example, for $\\int \\sin^2(x) dx$, all the methods can give exact results. To see why, you can rewrite $\\sin^2 x=(1-\\cos(2x))/2$. If you have an odd number of sample points, the cos term averages to zero and the remaining term is a constant, which all three methods will fit exactly. \n",
|
||
"- It is possible to derive analytic expressions for the errors in the different approximations (these are known as [Euler Maclaurin](https://en.wikipedia.org/wiki/Euler–Maclaurin_formula) formulae). For the trapezoidal rule, the leading term is \n",
|
||
"\n",
|
||
"$$\\epsilon \\approx {1\\over 12} (\\Delta x)^2 \\left[f^\\prime(a) - f^\\prime(b)\\right]$$\n",
|
||
"\n",
|
||
"and for Simpson's rule,\n",
|
||
"\n",
|
||
"$$\\epsilon \\approx {1\\over 180} (\\Delta x)^4 \\left[f^{\\prime\\prime\\prime}(a) - f^{\\prime\\prime\\prime}(b)\\right].$$\n",
|
||
"\n",
|
||
"These are plotted as dotted lines in the Figure, you can see that they agree remarkably well with the numerical results. (Note that when plotting we use the fact that the derivative terms are equal to unity for $\\cos(x)$). "
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "a5381b56",
|
||
"metadata": {},
|
||
"source": [
|
||
"## Gaussian quadrature"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 5,
|
||
"id": "7eb969c8",
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"quad gives I = 1.493648265624854 1.6582826951881447e-14\n",
|
||
"1 1 2.0 0.33900332898208696\n",
|
||
"2 3 1.4330626211475785 -0.0405621898217985\n",
|
||
"3 5 1.4986795956600294 0.0033684838331537233\n",
|
||
"4 7 1.4933346224495394 -0.00020998462792936129\n",
|
||
"5 9 1.493663920702629 1.048110062810725e-05\n",
|
||
"6 11 1.4936476141506054 -4.361630937756038e-07\n",
|
||
"7 13 1.493648288869414 1.5562271646873954e-08\n",
|
||
"8 15 1.4936482648990144 -4.859508274756656e-10\n",
|
||
"9 17 1.4936482656450039 1.3490379379714472e-11\n"
|
||
]
|
||
},
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGjCAYAAADUwuRbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABFWUlEQVR4nO3dZ3gU5cLG8f/sppEKIRAIJPReEkoIIYCiUUSkiAoqagAFRQQ0loN6jujxVTxybEioFlBUEJUgKKKiSEkwEAjSewmEJNQ0SN19P3jMeXlBaZtMsrl/17Ufdnayc496ufc188zzGHa73Y6IiIhIJWcxO4CIiIiII6jUiIiIiFNQqRERERGnoFIjIiIiTkGlRkRERJyCSo2IiIg4BZUaERERcQoqNSIiIuIUVGpERETEKajUiIiIiFNQqRERERGnUKlKzdKlS2nRogXNmjXjvffeMzuOiIiIVCBGZVnQsri4mNatW/Pzzz/j5+dHp06dSEhIoGbNmmZHExERkQqg0lypSUpKok2bNtSrVw9vb2/69OnD999/b3YsERERqSBcyutAq1atYvLkySQnJ3Ps2DEWLVrEwIEDz9snLi6OyZMnk56eTmhoKO+++y5dunQBIC0tjXr16pXuW69ePY4ePXrZx7fZbKSlpeHj44NhGA45JxERESlbdrudnJwcgoKCsFj++lpMuZWavLw8QkNDGTFiBIMGDbrg8wULFhAbG8uMGTOIiIjg7bffpnfv3uzatYvatWtf8fEKCgooKCgofX/06FFat259TecgIiIi5khNTaV+/fp/uU+5lZo+ffrQp0+fP/38zTffZOTIkQwfPhyAGTNm8M033/DBBx8wYcIEgoKCzrsyc/To0dKrOBczadIkXnrppQu2p6am4uvrew1nIiIiIuUlOzub4OBgfHx8LrmvKQOFDcM47/ZTYWEhnp6efPHFF+fdkoqJieHMmTMsXryY4uJiWrVqxcqVKy9roPD/v1Lzxz+UrKwslRoREZFKIjs7Gz8/v8v6/S63KzV/5cSJE5SUlBAYGHje9sDAQHbu3AmAi4sLb7zxBr169cJms/HMM8/85ZNP7u7uuLu7l2luERERqTgqRKm5XP3796d///5mxxAREZEKqEI80h0QEIDVaiUjI+O87RkZGdSpU8ekVCIiIlKZVIhS4+bmRqdOnVixYkXpNpvNxooVK4iMjDQxmYiIiFQW5Xb7KTc3l71795a+P3DgACkpKfj7+xMSEkJsbCwxMTF07tyZLl268Pbbb5OXl1f6NJSIiIjIXym3UrNhwwZ69epV+j42Nhb4/QmnOXPmMGTIEI4fP84LL7xAeno6YWFhfPfddxcMHhYRERG5mEqz9tO1upJHwkRERKRiuJLf7woxpkZERETkWqnUiIiIiFNQqRERERGnoFIjIiIiTkGlRkRERJyCSs01stvtjPlkIws3pFJiqxIPkomIiFRIKjXX6McdmXyz5RhPf/Ebvd9exXdb06kiT8mLiIhUKCo116hHswCeu7Ul1T1d2ZuZyyPzkhk4LYG1e0+YHU1ERKRK0eR7jvr+/CJmr9rP+2sOcLawBICopjV5pndLQoOrO/x4IiIiVcGV/H6r1DjY8ZwC4n7eyye/HqKo5Pd/tLe0qcNTvZvTtLZPmR1XRETEGanUXER5L5OQeuosb/+4h0WbjmCzg8WAOzrW5/GbmlOverUyP76IiIgzUKm5CLPWftqdkcO/l+/i++0ZALhZLQztGsKYXk0J8HYvtxwiIiKVkUrNRZi9oOWmw6d5/btdJO4/CYCXm5UHezRmZI9G+Hi4lnseERGRykCl5iLMLjXw+5w2a/ae4PXvdrHlaBYANTxdGdOrKfd1bYCHq9WUXCIiIhWVSs1FVIRS8we73c53W9P59/e72Hc8D4C6fh6Mv7EZd3aqj4tVT9qLiIiASs1FVaRS84fiEhtfbTzK2z/uJi0rH4DGAV48eXML+rStg8VimJxQRETEXCo1F1ERS80f8otK+OTXw8T9vJdTeYUAtK3ny9O9W9KzWQCGoXIjIiJVk0rNRVTkUvOHnPwi3l9zgPdWHyC3oBiAro39eeaWlnQMqWFyOhERkfKnUnMRlaHU/OFkbgHTVu7j43WHKCy2ARDdKpCne7egRR1N4CciIlWHSs1FVKZS84e0M+d458c9LExOxWYHw4Dbw+rxxE3NCfb3NDueiIhImVOpuYjKWGr+sDczlzd/2MW3W9IBcLUa3NslhMduaEYtH03gJyIizkul5iIqc6n5w5YjWby+fCer9/y+Ang1VysjujdkVM8m+FXTBH4iIuJ8VGouwhlKzR8S9v0+gV9K6hkA/Kq5Mvr6JsRENqSamybwExER56FScxHOVGrg9wn8vt+ewb+X72JPZi4AtX3cGR/djMGdg3HVBH4iIuIEVGouwtlKzR9KbHbiNx3lrR93c+T0OQAa1vTkiZua0699kCbwExGRSk2l5iKctdT8oaC4hM9+PczUn/dyIvf3Cfxa1fXlmd4tuL5FLU3gJyIilZJKzUU4e6n5Q15BMR+uPcDMX/aT858J/MIb1uCZW1oS3tDf5HQiIiJXRqXmIqpKqfnD6bxCZvyyjzkJByn4zwR+vVrU4uneLWkd5PznLyIizkGl5iKqWqn5Q3pWPlN+2sOC9amU2H7/V90/NIjYm5rTMMDL5HQiIiJ/TaXmIqpqqfnDgRN5vPnDbpZsTgPAxWIwJDyYcTc2I9DXw+R0IiIiF6dScxFVvdT8YevRLP79/S5W7joOgIerhZhuDRl9XROqe7qZnE5EROR8KjUXoVJzvqQDp3j9u51sOHQaAB8PFx65rgnDoxri6eZicjoREZHfqdRchErNhex2Oz/vyuT173axMz0HgABvd8be0JR7uoTg5qIJ/ERExFwqNRehUvPnbDY7S35L443vd3P41FkAgv2r8UR0cwaE1cOqCfxERMQkKjUXoVJzaYXFNhZsSGXKij0czykAoEWgD0/e3JybWgdqAj8RESl3KjUXoVJz+c4VljAn4SDTV+4lO//3Cfw6hlTnnwPa0raen8npRESkKlGpuQiVmiuXdbaImav28cHaA+QX2bAYENOtIU/e3AJvdw0mFhGRsqdScxEqNVcvIzuf//lmR+kcN3V8PZjYrzW3tK2jW1IiIlKmruT3W4+3yCUF+nrw7j0d+GhEFxrU9CQ9O5/Rn2zkwbkbSP3PwGIRERGzqdTIZevZvBbLH+/J2Bua4mo1+GlnJje99QvTV+6jqMRmdjwREaniVGrkini4Wnny5hYsG9+TiEb+5BfZ+Nd3O+k7ZTUbDp4yO56IiFRhKjVyVZrW9mb+qK68cVco/l5u7M7I5c4ZiUz48jfOnC00O56IiFRBKjVy1QzD4I5O9VkRex13hwcDMH99Kje88QtfJh+hioxBFxGRCkKlRq5ZDS83XrujPQsfiaR5oDen8gp5cuFm7pm9jr2ZuWbHExGRKkKlRhwmvKE/S8f24G+3tMTD1cK6/afo884q3vx+F/lFJWbHExERJ6dSIw7l5mJh9PVN+OGJ6+jVohZFJXam/LSX3m+vYtXu42bHExERJ6ZSI2Ui2N+TD4aFM31oRwJ93Tl08iwPfJDE2M82kZmTb3Y8ERFxQio1UmYMw6BPu7r8GHsdw6MaYjFgyeY0bnzjFz5OPEiJTQOJRUTEcbRMgpSbLUeyeD5+C78dyQIgNLg6rwzUIpkiIvLntEyCVEjt6vux6NEoXurfBm93FzannqH/1DW8vHQ7uQXFZscTEZFKrtKUmtTUVK6//npat25N+/btWbhwodmR5CpYLQYx3Rqy4snr6Nu+LjY7vL/mADe9+QvfbU3X3DYiInLVKs3tp2PHjpGRkUFYWBjp6el06tSJ3bt34+XldVl/r9tPFdPKXZn8Y/FWUk+dAyC6VW1e7N+G+jU8TU4mIiIVgVPefqpbty5hYWEA1KlTh4CAAE6d0lpDld31LWrz/ePXMaZXE1ytBj/uyOSmN1cx8xctkikiIlfGYaVm1apV9OvXj6CgIAzDID4+/oJ94uLiaNiwIR4eHkRERJCUlHRVx0pOTqakpITg4OBrTC0VQTU3K0/3bsm343rQpaE/54pKmLRsJ/3eXUPyIRVXERG5PA4rNXl5eYSGhhIXF3fRzxcsWEBsbCwTJ05k48aNhIaG0rt3bzIzM0v3CQsLo23bthe80tLSSvc5deoUDzzwALNmzXJUdKkgmgX6sODhrrx+Z3tqeLqyMz2HO6Yn8uxXW7RIpoiIXFKZjKkxDINFixYxcODA0m0RERGEh4czdepUAGw2G8HBwYwdO5YJEyZc1vcWFBRw0003MXLkSO6///5L7ltQUFD6Pjs7m+DgYI2pqSRO5RUy6dsdLEw+AkBNLzf+flsrBobVwzAMk9OJiEh5qXBjagoLC0lOTiY6Ovq/B7ZYiI6OJjEx8bK+w263M2zYMG644YZLFhqASZMm4efnV/rSrarKxd/Ljcl3hbJgVFea1vbmZF4hTyzYzND3fmXfcS2SKSIiFyqXUnPixAlKSkoIDAw8b3tgYCDp6emX9R1r165lwYIFxMfHExYWRlhYGFu2bPnT/Z999lmysrJKX6mpqdd0DmKOiMY1+XZcD57u3QJ3FwsJ+07S5+3VvPnDbi2SKSIi53ExO8Dl6t69Ozbb5T8N4+7ujru7exkmkvLi5mJhTK+m9GsfxD8Wb+WX3ceZsmIPSzan8fKAtnRvFmB2RBERqQDK5UpNQEAAVquVjIyM87ZnZGRQp06d8oggTiCkpidzhocTd29Havu4c+BEHve9/yvj52/ieE7Bpb9AREScWrmUGjc3Nzp16sSKFStKt9lsNlasWEFkZGR5RBAnYRgGfdvX5ccnryMmsgGGAYtT0rjhjZXMW3cImxbJFBGpshxWanJzc0lJSSElJQWAAwcOkJKSwuHDhwGIjY1l9uzZzJ07lx07djB69Gjy8vIYPny4oyJIFeLr4cpLA9qyeEwUbev5kpNfzN/jt3LHjAS2p2WbHU9EREzgsEe6V65cSa9evS7YHhMTw5w5cwCYOnUqkydPJj09nbCwMKZMmUJERIQjDn9JWibBeRWX2Ph43SHe+H43uQXFWC0GI6Ia8nh0c7zcK82wMRERuYgr+f2uNGs/XSuVGueXnpXPS0u2sWzr70/UBfl58GL/NtzcRuO2REQqqwo3T41Ieajj58H0+zrxwbDO1K9RjbSsfEZ9nMzIjzZw9Mw5s+OJiEgZU6kRp3NDy0B+eOI6HrmuCS4Wgx+2Z3DTm78we9V+LZIpIuLEVGrEKVVzszKhT0u+GdeDzg1qcLawhFe+3UG/d9ew8fBps+OJiEgZUKkRp9aijg+fPxzJv+5oR/XSRTITeH7RFrLOFpkdT0REHEilRpyexWIwJDyEFbHXcUfH+tjt8Mmvh7nxzZUsTjlKFRkrLyLi9FRqpMqo6e3OG4ND+WxkVxrX8uJEbiHj56dw//tJHDiRZ3Y8ERG5Rio1UuVENqnJsvE9ePKm5ri5WFiz9wS9315F3M97KdGMxCIilZZKjVRJ7i5Wxt7YjO8f70mPZgEUFtuYvHwXg2cmcuikrtqIiFRGKjVSpTUM8OKjEV2YfGd7vN1dSD50mj7vrObTXw9rrI2ISCWjUiNVnmEY3NU5mGXje9ClkT9nC0t4btEWRsxZT2Z2vtnxRETkMqnUiPxHsL8n80d25flbW+FmtfDzruP0fnsV3245ZnY0ERG5DCo1Iv+HxWIwsmdjloztTuu6vpw+W8Sjn2zkiQUpZJ3TvDYiIhWZSo3IRbSo40P8mCjG9GqCxYBFm47S5+1VrN17wuxoIiLyJ1RqRP6Em4uFp3u3ZOEjkTSo6UlaVj5D3/uVl5ZsI7+oxOx4IiLy/6jUiFxCpwb+fDuuB/dGhADw4dqD3PbuGrYcyTI5mYiI/F8qNSKXwcvdhVdvb8eHw8Kp5ePO3sxcbp+2lnd+3EOxVv4WEakQVGpErkCvlrVZ/nhPbm1Xh2Kbnbd+3M0dMxLZfzzX7GgiIlWeSo3IFfL3ciPu3o68PSQMHw8XNqee4dYpq/ko8aAm7BMRMZFKjchVMAyDgR3qsfzxnkQ1rUl+kY0XFm/jgQ+SSM/ShH0iImZQqRG5BkHVq/HxiAgm9muNu4uF1Xt+Xxzz681pZkcTEalyVGpErpHFYjA8qhHfjOtOu3p+ZJ0rYtxnm3js042cOVtodjwRkSpDpUbEQZrW9uGrR7sx/sZmWC0GS387Ru+3V/HL7uNmRxMRqRJUakQcyNVq4YmbmvPl6G40DvAiI7uAmA+S+Ef8Vs4WFpsdT0TEqanUiJSBsODqfDOuBzGRDQD4eN0h+k5Zw6bDp01OJiLivFRqRMpINTcrLw1oy8cPdqGOrwcHTuRx54xE3vx+F0WasE9ExOFUakTKWI9mtVj+eE8GhAVRYrMz5ae93D5tLXsycsyOJiLiVFRqRMqBn6cr79zdgan3dsCvmitbj2bT9901vL/mADabJuwTEXEElRqRcnRb+yC+f6InPZvXorDYxstLt3Pf+79y9Mw5s6OJiFR6KjUi5SzQ14O5w8N5eWBbqrlaSdh3klveWsVXG49omQURkWugUiNiAsMwuL9rA74d34Ow4OrkFBQT+/lmRs/byKk8TdgnInI1VGpETNQowIsvHonkqZub42Ix+G5bOje/tYqfdmaYHU1EpNJRqRExmYvVwmM3NCN+TBTNantzIreAEXM28OxXv5FXoAn7REQul0qNSAXRtp4fS8Z258HujQD4LCmVPu+sZsPBUyYnExGpHFRqRCoQD1cr/7itNZ+OjKBe9WocPnWWwTMT+dd3Oyks1oR9IiJ/RaVGpALq1iSAZY/34I6O9bHZYfrKfQyIW8vO9Gyzo4mIVFgqNSIVlK+HK28MDmXGfR3x93Jjx7Fs+r+7llmr9lGiCftERC6gUiNSwd3Sti7fPd6DG1vWprDExqvf7uSe2etIPXXW7GgiIhWKSo1IJVDbx4P3Yjrz2qB2eLlZSTpwilveXsXn61M1YZ+IyH+o1IhUEoZhcHeXEJaN70l4wxrkFZbwzJe/MerjZE7kFpgdT0TEdCo1IpVMSE1P5o+KZEKflrhaDX7YnkHvt1bx/bZ0s6OJiJhKpUakErJaDB65rglfP9adlnV8OJlXyKiPk3l64WZy8ovMjiciYgqVGpFKrFVdXxY/FsUj1zXBMGBh8hFueXs16/afNDuaiEi5U6kRqeTcXaxM6NOSzx+OJNi/GkfPnOOe2et45Zvt5BeVmB1PRKTcqNSIOInwhv4sG9+Tu8ODsdth9uoDDJi6lm1pWWZHExEpFyo1Ik7E292F1+5oz3sPdCbA241dGTkMjFtL3M97NWGfiDg9lRoRJxTdOpDlj/ekd5tAikrsTF6+iyEzEzl65pzZ0UREyoxKjYiTquntzoz7OvHGXaH4uLuw4dBp+k5ZzY/bM8yOJiJSJlRqRJyYYRjc0ak+34zrQWh9P86cLeKhjzbwP0u3a9VvEXE6KjUiVUBITU8WPtKNEVGNAHhvzQEGz0zU+lEi4lRUakSqCDcXCy/0a82s+zvh6+FCSuoZ+k5ZzXLNRCwiTkKlRqSKublNHb4d34Ow4Opk5xfz8MfJvLRkm25HiUilV+lKzdmzZ2nQoAFPPfWU2VFEKq36NTz5/OFIRvb4/XbUh2sPcueMBA6f1O0oEam8Kl2peeWVV+jatavZMUQqPTcXC8/3bc37MZ2p7unKb0ey6DtlNcu2HDM7mojIValUpWbPnj3s3LmTPn36mB1FxGnc2CqQb8f1oFODGuQUFDP6k428sHirllgQkUrHYaVm1apV9OvXj6CgIAzDID4+/oJ94uLiaNiwIR4eHkRERJCUlHRFx3jqqaeYNGmSgxKLyB+Cqldj/qiuPHJdEwA+SjzEHdMTOHgiz+RkIiKXz2GlJi8vj9DQUOLi4i76+YIFC4iNjWXixIls3LiR0NBQevfuTWZmZuk+YWFhtG3b9oJXWloaixcvpnnz5jRv3txRkUXk/3C1WpjQpyUfDg/H38uNbWnZ3PbuGpZsTjM7mojIZTHsdrvDF4QxDINFixYxcODA0m0RERGEh4czdepUAGw2G8HBwYwdO5YJEyZc8jufffZZ5s2bh9VqJTc3l6KiIp588kleeOGFi+5fUFBAQUFB6fvs7GyCg4PJysrC19f32k5QxMmlZ+Uz7rNNJB08BcC9ESG8cFtrPFytJicTkaomOzsbPz+/y/r9LpcxNYWFhSQnJxMdHf3fA1ssREdHk5iYeFnfMWnSJFJTUzl48CD//ve/GTly5J8Wmj/29/PzK30FBwdf83mIVBV1/Dz4dGQEj/VqimHAp78eZmDcWvYdzzU7mojInyqXUnPixAlKSkoIDAw8b3tgYCDp6WUz8dezzz5LVlZW6Ss1NbVMjiPirFysFp7q3YKPRnShppcbO9Nz6PfuGuI3HTU7mojIRbmYHeBqDBs27JL7uLu74+7uXvZhRJxcj2a1WDa+B+Pmb2Ld/lM8viCFxH0nebF/G6q56XaUiFQc5XKlJiAgAKvVSkbG+asDZ2RkUKdOnfKIICLXoLavB5881JXxNzbDMGDBhlQGxq1lb2aO2dFEREqVS6lxc3OjU6dOrFixonSbzWZjxYoVREZGlkcEEblGVovBEzc155MHIwjwdmdXRg793l3LF8lHzI4mIgI4sNTk5uaSkpJCSkoKAAcOHCAlJYXDhw8DEBsby+zZs5k7dy47duxg9OjR5OXlMXz4cEdFEJFy0K1pAMvG96B70wDOFZXw1MLNPPn5Zs4WFpsdTUSqOIc90r1y5Up69ep1wfaYmBjmzJkDwNSpU5k8eTLp6emEhYUxZcoUIiIiHHH4S7qSR8JE5NJKbHam/byXt37cjc0OTWt7E3dvR1rU8TE7mog4kSv5/S6TeWoqIpUakbKxbv9Jxs/fREZ2AR6uFl7q34bBnYMxDMPsaCLiBCrcPDUi4ry6Nq7Jt+N60LN5LfKLbPztyy08sSCFvALdjhKR8qVSIyLXrKa3O3OGhfPMLS2wWgziU9Lo9+4adhzLNjuaiFQhKjUi4hAWi8Gj1zdl/qiu1PXzYP+JPAbEreXTXw9TRe5yi4jJVGpExKHCG/rzzbge3NCyNoXFNp5btIVx81PIyS8yO5qIODmVGhFxOH8vN957oDPP3doSF4vBks2/347aejTL7Ggi4sRUakSkTFgsBqN6NmHBw5HUq16NgyfPMmhaAh8lHtTtKBEpEyo1IlKmOjWowTfjuhPdKpDCEhsvLN7GmE83kq3bUSLiYCo1IlLmqnu6MfuBTvy9bytcrQbfbkmn75TV/HbkjNnRRMSJqNSISLkwDIOHejRm4SPdqF+jGqmnznHH9AQ+WHNAt6NExCFUakSkXIUFV+ebcT3o3SaQohI7/1y6nYc/TibrrG5Hici1UakRkXLnV82VGfd14sV+rXGzWvh+ewa3TlnNpsOnzY4mIpWYSo2ImMIwDIZFNeLL0d0I8ffk6Jlz3DUjkfdW79ftKBG5Kio1ImKqdvX9WDquO33b1aXYZud/vtnByI82cOZsodnRRKSSUakREdP5ergy9d4OvDywLW4uFn7ckcmt76wm+dAps6OJSCWiUiMiFYJhGNzftQGLHu1GowAv0rLyGTxzHTN+2YfNpttRInJpKjUiUqG0CfJjydju9A8NosRm57VlOxkxdz2n8nQ7SkT+mkqNiFQ43u4uvHN3GJMGtcPdxcLKXce59Z3VJB3Q7SgR+XMqNSJSIRmGwT1dQogfE0XjWl6kZ+dz96xEpv60R7ejROSiVGpEpEJrVdeXJY91Z1CHetjs8O/vdxPzYRIncgvMjiYiFYxKjYhUeF7uLrwxOJTX72yPh6uF1XtOcOs7q0ncd9LsaCJSgajUiEilYBgGgzsH8/Vj3WlW25vMnAKGvreOd37cQ4luR4kIKjUiUsk0D/Rh8WNR3NWpPjY7vPXjbu5//1cyc/LNjiYiJlOpEZFKx9PNhcl3hfLm4FCquVpJ2HeSW99Zw7r9uh0lUpWp1IhIpTWoY32WjO1Oi0AfTuQWMPS9X3l/zQGtHSVSRanUiEil1rS2N/Fjori9Qz1KbHZeXrqdxxekcK6wxOxoIlLOVGpEpNKr5mblzcGhTOzXGqvFYHFKGrdPW8vhk2fNjiYi5UilRkScgmEYDI9qxKcPRRDg7cbO9Bxue3c1K3dlmh1NRMqJSo2IOJWIxjVZOrYHYcHVyc4vZvic9ZqFWKSKUKkREadTx8+DBQ935d6IEOz/mYX44XnJ5OQXmR1NRMqQSo2IOCV3Fyuv3t6Of93RDjerhR+2ZzAgbi17M3PMjiYiZUSlRkSc2pDwED5/JJK6fh7sP57HgKlr+W7rMbNjiUgZUKkREacXFlydJWO707WxP3mFJTwybyOvf7dTyyuIOBmVGhGpEgK83Zn3YAQPdW8EwLSV+xj2YRKn8wpNTiYijqJSIyJVhovVwt9va807d4eVrvbdb+oatqVlmR1NRBxApUZEqpwBYfVY9GgUIf6eHDl9jjumJ7Bo0xGzY4nINVKpEZEqqVVdX5Y81p3rW9Qiv8jGEws28+LX2ygqsZkdTUSukkqNiFRZfp6uvB8TztgbmgIwJ+EgQ2f/SmZOvsnJRORqqNSISJVmtRg8eXMLZt3fCW93F5IOnqLfu2vYePi02dFE5Aqp1IiIADe3qcPix6JoWtubjOwChsxM5JNfD2G367FvkcpCpUZE5D+a1PImfkwUfdrWoajEzvOLtjLhyy3kF5WYHU1ELoNKjYjI/+Ht7sK0oR155pYWWAxYsCGVITMTSTtzzuxoInIJKjUiIv+PYRg8en1T5gzvQnVPVzYfyaLfu2tI3HfS7Ggi8hdUakRE/kTP5rVY8lh3Wtf15WReIfe9/yvvrd6vcTYiFZRKjYjIXwj29+TL0d0Y1KEeJTY7//PNDsbNT+FsYbHZ0UTk/1GpERG5hGpuVt4YHMqL/VrjYjFYsjmNQdMSOHQyz+xoIvJ/qNSIiFwGwzAYFtWIT0d2JcDbnZ3pOfR7dw0/78o0O5qI/IdKjYjIFejSyJ+lY7vTIaQ62fnFjJiznikr9mCzaZyNiNlUakRErlAdPw/mj+rK0IgQ7HZ484fdjPo4mez8IrOjiVRpKjUiIlfB3cXKK7e34/U72uPmYuHHHRkMnLqWPRk5ZkcTqbJUakRErsHg8GAWPhxJkJ8H+0/kMSBuLd9uOWZ2LJEqqVKVmgMHDtCrVy9at25Nu3btyMvTkwciYr7Q4OosGdudyMY1OVtYwqOfbOS1ZTsp0TgbkXJVqUrNsGHD+Oc//8n27dv55ZdfcHd3NzuSiAgANb3d+fjBLozq2RiAGb/sI+aDJE7nFZqcTKTqqDSlZtu2bbi6utKjRw8A/P39cXFxMTmViMh/uVgtPHdrK969pwPVXK2s2XuC295dw9ajWWZHE6kSHFZqVq1aRb9+/QgKCsIwDOLj4y/YJy4ujoYNG+Lh4UFERARJSUmX/f179uzB29ubfv360bFjR1599VVHRRcRcah+oUEsGtONBjU9OXrmHHdMT+CrjUfMjiXi9BxWavLy8ggNDSUuLu6iny9YsIDY2FgmTpzIxo0bCQ0NpXfv3mRm/nfiqrCwMNq2bXvBKy0tjeLiYlavXs20adNITEzkhx9+4IcffnBUfBERh2pZx5evH+tOrxa1KCi2Efv5ZiYu3kpRic3saCJOy7CXwcpshmGwaNEiBg4cWLotIiKC8PBwpk6dCoDNZiM4OJixY8cyYcKES35nYmIiL774IsuXLwdg8uTJADz99NMX3b+goICCgoLS99nZ2QQHB5OVlYWvr+/VnpqIyBWx2ey8s2IP76zYA0B4wxrEDe1IbR8Pk5OJVA7Z2dn4+fld1u93uYypKSwsJDk5mejo6P8e2GIhOjqaxMTEy/qO8PBwMjMzOX36NDabjVWrVtGqVas/3X/SpEn4+fmVvoKDg6/5PERErpTFYvDETc1574HO+Li7sP7gaW6bsobkQ6fNjibidMql1Jw4cYKSkhICAwPP2x4YGEh6evplfYeLiwuvvvoqPXv2pH379jRr1ozbbrvtT/d/9tlnycrKKn2lpqZe0zmIiFyL6NaBLH4sima1vcnMKeDuWYnMW3eIMrhYLlJlVarHh/r06UOfPn0ua193d3c98i0iFUrjWt7Ej4ni6S828+2WdP4ev5XNqWd4eWBbPFytZscTqfTK5UpNQEAAVquVjIyM87ZnZGRQp06d8oggIlIheLm7EHdvR57t0xKLAQuTjzB4ZiJHz5wzO5pIpVcupcbNzY1OnTqxYsWK0m02m40VK1YQGRlZHhFERCoMwzB4+LomfDQighqervx2JIt+764hYe8Js6OJVGoOKzW5ubmkpKSQkpIC/L6kQUpKCocPHwYgNjaW2bNnM3fuXHbs2MHo0aPJy8tj+PDhjoogIlKpdG8WwNePdadtPV9O5RVy3/u/MnvVfo2zEblKDnuke+XKlfTq1euC7TExMcyZMweAqVOnMnnyZNLT0wkLC2PKlClEREQ44vCXdCWPhImIlKf8ohKeX7SVL/8zQd9t7evy+p3t8XSrVMMeRcrElfx+l8k8NRWRSo2IVGR2u5156w7x0pLtFNvstAj0Yeb9nWgY4GV2NBFTVbh5akRE5K8ZhsH9kQ2ZP6ortXzc2ZWRQ7+pa/hpZ8al/1hEAJUaEZEKpXNDf5aO7U6nBjXIyS/mwbkbeOfHPdhsVeKiusg1UakREalgAn09+GxkV+7v2gC7Hd76cTejPt5A1rkis6OJVGgqNSIiFZCbi4WXB7Zl8p3tcXOx8OOOTG6ftpYDJ/LMjiZSYanUiIhUYHd1DubLR7oR5OfB/uN5DIxby5o9ms9G5GJUakREKrh29f2IfyyKDiHVyTpXRMyHSXyceNDsWCIVjkqNiEglUNvn93E2gzrUo8Rm5x+Lt/H3+C0UldjMjiZSYajUiIhUEh6uVt4YHMrfbmmJYcC8dYeJ+SCJM2cLzY4mUiGo1IiIVCKGYTD6+ibMur8zXm5WEvadZGDcWvZm5pgdTcR0KjUiIpXQTa0D+fLRbtSvUY2DJ89ye1wCK3dlmh1LxFQqNSIilVTLOr4sHhNFeMMa5BQUM2LOet5fc0ALYkqVpVIjIlKJ1fR255OHujK4c31sdnh56Xae/WoLhcUaQCxVj0qNiEgl5+Zi4V93tOfvfVthMWD++lTue/9XTuVpALFULSo1IiJOwDAMHurRmPeHhePj7kLSgVP0n7qGXekaQCxVh0qNiIgT6dWiNovGdKNBTU+OnD7HoGlrWbFDK31L1aBSIyLiZJrW9iH+0SgiG9ckr7CEhz7awMxf9mkAsTg9lRoRESdUw8uNjx7swtCIEOx2mLRsJ08u3ExBcYnZ0UTKjEqNiIiTcrVa+J+BbfnngDZYLQZfbTzKPbPWcTynwOxoImVCpUZExIkZhsEDkQ2ZMzwcXw8XNh4+w4Cpa9iWlmV2NBGHU6kREakCejSrRfyYKBoHeJGWlc+d0xP5bmu62bFEHEqlRkSkimhcy5tFj0bRo1kA54pKeGReMlN/2qMBxOI0VGpERKoQP09XPhwWzrBuDQH49/e7GT8/hfwiDSCWyk+lRkSkinGxWnixfxtevb0dLhaDrzenMWRmIhnZ+WZHE7kmKjUiIlXUvREhfPxgBNU9Xdl8JIsBU9ey5YgGEEvlpVIjIlKFRTapyeIxUTSt7U16dj53zUxg6W9pZscSuSoqNSIiVVyDml589Wg3erWoRX6Rjcc+3cRbP+zGZtMAYqlcVGpERARfD1feiwlnZI9GALyzYg+PfbaRc4UaQCyVh0qNiIgAYLUYPN+3Na/f2R5Xq8G3W9K5a2YCx7LOmR1N5LKo1IiIyHkGdw7m05Fdqenlxtaj2fSfupZNh0+bHUvkklRqRETkAuEN/YkfE0XLOj4czylgyKx1xG86anYskb+kUiMiIhcV7O/JF6O7Ed0qkMJiG48vSOH173ZqALFUWCo1IiLyp7zdXZh1fydGX98EgGkr9/HwvGTyCopNTiZyIZUaERH5SxaLwd9uaclbQ0Jxc7Hww/YM7piewJHTZ82OJnIelRoREbkst3eoz/xRXQnwdmdneg4Dpq5lw8FTZscSKaVSIyIil61jSA2+fiyKNkG+nMwr5J7Z61i4IdXsWCKASo2IiFyhoOrVWPhIJH3a1qGoxM7TX/zGq9/uoEQDiMVkKjUiInLFPN1ciLu3I+NubAbArFX7GfnRBnLyi0xOJlWZSo2IiFwVi8Ug9qbmvHtPB9xdLPy0M5NB0xI4fFIDiMUcKjUiInJN+oUGsfCRSAJ93dmTmcuAuDWs23/S7FhSBanUiIjINWtfvzpfP9ad0Pp+nD5bxH3v/cpnSYfNjiVVjEqNiIg4RKCvBwsejqRfaBDFNjvPfrWFF7/eRnGJzexoUkWo1IiIiMN4uFqZcncYT97UHIA5CQcZPmc9Wec0gFjKnkqNiIg4lGEYjL2xGTPu60g1Vyur95zg9mlr2X881+xo4uRUakREpEzc0rYuX4yOJMjPg/3H8xgYt5Y1e06YHUucmEqNiIiUmTZBfix+rDsdQ6qTnV9MzIdJfJR40OxY4qRUakREpEzV8nHn05FdGdSxHiU2Oy8s3sbf47dQpAHE4mAqNSIiUuY8XK28cVcoE/q0xDBg3rrDxHyQxJmzhWZHEyeiUiMiIuXCMAweua4Js+/vjJeblYR9JxkQt5a9mTlmRxMnoVIjIiLlKrp1IF8+2o36Napx6ORZbo9LYOWuTLNjiRNQqRERkXLXso4vi8dE0aWhPzkFxYyYs5731xzAbtdK33L1VGpERMQUNb3dmfdQBEM6B2Ozw8tLt/PCYs1ALFevUpWat956izZt2tC6dWvGjRunRi8iUsm5uVh47Y52/L1vKwwDPl53iJEfbSC3oNjsaFIJVZpSc/z4caZOnUpycjJbtmwhOTmZdevWmR1LRESukWEYPNSjMdOHdsLD1cLPu44zeEYi6Vn5ZkeTSqbSlBqA4uJi8vPzKSoqoqioiNq1a5sdSUREHOSWtnWYPyqSAG83th/LZmDcWranZZsdSyoRh5WaVatW0a9fP4KCgjAMg/j4+Av2iYuLo2HDhnh4eBAREUFSUtJlf3+tWrV46qmnCAkJISgoiOjoaJo0aeKo+CIiUgGEBVdn0aNRNK3tTXp2PnfNSOBnPRkll8lhpSYvL4/Q0FDi4uIu+vmCBQuIjY1l4sSJbNy4kdDQUHr37k1m5n//Yw0LC6Nt27YXvNLS0jh9+jRLly7l4MGDHD16lISEBFatWuWo+CIiUkEE+3vy5ehudGtSk7zCEh6au4F56w6ZHUsqAcNeBqNtDcNg0aJFDBw4sHRbREQE4eHhTJ06FQCbzUZwcDBjx45lwoQJl/zOhQsXsnLlytLSNHnyZOx2O88888xF9y8oKKCgoKD0fXZ2NsHBwWRlZeHr63sNZyciIuWhsNjGc4u28EXyEQBG9WzMhFtaYrEYJieT8pSdnY2fn99l/X6Xy5iawsJCkpOTiY6O/u+BLRaio6NJTEy8rO8IDg4mISGB/Px8SkpKWLlyJS1atPjT/SdNmoSfn1/pKzg4+JrPQ0REyo+bi4XJd7bnyZuaAzBr1X7GfLqR/KISk5NJRVUupebEiROUlJQQGBh43vbAwEDS09Mv6zu6du3KrbfeSocOHWjfvj1NmjShf//+f7r/s88+S1ZWVukrNTX1ms5BRETKn2EYjL2xGe/cHYab1cKyrencPWsdJ3ILLv3HUuW4mB3gSrzyyiu88sorl7Wvu7s77u7uZZxIRETKw4CwetT1q8aojzeQknqG26et5cNh4TSt7WN2NKlAyuVKTUBAAFarlYyMjPO2Z2RkUKdOnfKIICIilVyXRv58NbobDWp6knrqHIOmJZCw74TZsaQCKZdS4+bmRqdOnVixYkXpNpvNxooVK4iMjCyPCCIi4gQa1/Jm0aNRdGpQg+z8YmI+SCodSCzisFKTm5tLSkoKKSkpABw4cICUlBQOHz4MQGxsLLNnz2bu3Lns2LGD0aNHk5eXx/Dhwx0VQUREqgB/Lzc+eSiC29rXpajEzlMLN/PmD7u1dI44bkzNhg0b6NWrV+n72NhYAGJiYpgzZw5Dhgzh+PHjvPDCC6SnpxMWFsZ33313weBhERGRS/FwtTLl7g6E+HsybeU+pqzYQ+qps7x2RzvcXaxmxxOTlMk8NRXRlTznLiIilcf8pMM8H7+VEpudLo38mXV/J6p7upkdSxykws1TIyIiUlbu7hLCnOHh+Li7kHTgFIOmJ3DoZJ7ZscQEKjUiIlLp9WhWiy9GdyPIz4P9x/O4fVoCyYdOmx1LyplKjYiIOIUWdXyIHxNFu3p+nMor5J7Z6/jmt2Nmx5JypFIjIiJOo7avBwse7kp0q9oUFtsY8+lGZvyyT09GVREqNSIi4lQ83VyYeX9nhnVrCMBry3by3KKtFJXYzA0mZU6lRkREnI7VYvBi/zZM7NcaiwGfJR1mxJz15OQXmR1NypBKjYiIOK3hUY2YeX9nqrlaWb3nBHfNSCTtzDmzY0kZUakRERGndlPrQD5/OJJaPu7sTM9hYNxath7NMjuWlAGVGhERcXrt6vsRPyaKFoE+ZOYUMHhmIit2ZFz6D6VSUakREZEqoV71aiwcHUmPZgGcLSxh5EcbmJtw0OxY4kAqNSIiUmX4erjywbBw7g4PxmaHiV9v459LtlNi0yPfzkClRkREqhRXq4VJg9rxzC0tAPhg7QEemZfM2cJik5PJtVKpERGRKscwDB69vinv3tMBNxcLP2zP4O5Z68jMyTc7mlwDlRoREamy+oUG8elDEdTwdOW3I1ncHpfA7owcs2PJVVKpERGRKq1zQ38WPRpFowAvjp45xx3TEliz54TZseQqqNSIiEiV1zDAi69Gd6NLQ39yCooZ9mESn69PNTuWXCGVGhEREaCGlxsfP9SFAWFBFNvsPPPlb0xevhObnoyqNFRqRERE/sPdxcrbQ8IYd0NTAOJ+3sf4BSnkF5WYnEwuh0qNiIjI/2EYBrE3t+D1O9vjYjFYsjmN+977lVN5hWZHk0tQqREREbmIwZ2D+WhEF3w8XNhw6DSDpq3lwIk8s2PJX1CpERER+RPdmgbw1ehu1KtejYMnz3L7tLWsP3jK7FjyJ1RqRERE/kKzQB/ix0QRWt+PM2eLGDr7VxanHDU7llyESo2IiMgl1PJxZ/6oSHq3CaSwxMb4+SnE/bwXu11PRlUkKjUiIiKXoZqblWlDO/FQ90YATF6+i799+RtFJTaTk8kfVGpEREQuk9Vi8PfbWvPygDZYDPh8wxGGfZhE1rkis6MJKjUiIiJX7P7IhrwX0xlPNytr957kzukJHDl91uxYVZ5KjYiIyFW4oWUgnz8cSaCvO3sycxkYl8Dm1DNmx6rSVGpERESuUtt6fsSPiaJlHR9O5BYwZFYiy7elmx2rylKpERERuQZ1/arxxehuXNe8FvlFNh6Zl8z7aw7oySgTqNSIiIhcI293F96P6czQiBDsdnh56XZe/HobxXoyqlyp1IiIiDiAi9XC/wxsy3O3tgRgbuIhRn2cTF5BscnJqg6VGhEREQcxDINRPZswfWhH3F0s/LQzk8EzE8nIzjc7WpWgUiMiIuJgfdrV5bNRXanp5ca2tGwGxq1lx7Fss2M5PZUaERGRMtAxpAaLHo2iSS0vjmXlc9eMRFbtPm52LKemUiMiIlJGQmp68tXoKLo29ie3oJgRc9bz1cYjZsdyWio1IiIiZcjP05W5I7rQPzSIYpud2M83M33lPj3yXQZUakRERMqYu4uVt4eEMbLH74th/uu7nbz49TZKbCo2jqRSIyIiUg4sFoPn+7bm731bAb8/8v3YpxvJLyoxOZnzUKkREREpRw/1aMy793TAzWph2dZ0Hng/iayzWuXbEVRqREREylm/0CDmjAjHx92FpIOnuHNGAmlnzpkdq9JTqRERETFBtyYBLBz931W+B01LYGe65rK5Fio1IiIiJmlZx5evHo2iWW1v0rN/n8smcd9Js2NVWio1IiIiJqpXvRoLH4kkvGENcvKLifkgiaW/pZkdq1JSqRERETFZdU83Pn4wglva1KGwxMbYzzbxwZoDZseqdFRqREREKgAPVytxQzvyQGQD7Hb459LtvPrtDmyay+ayqdSIiIhUEFaLwUv92/C3W1oCMGvVfp74PIXCYpvJySoHlRoREZEKxDAMRl/fhDcHh+JiMVicksbwOUnk5Gsum0tRqREREamABnWszwfDwvFys7J270kGz1xHRna+2bEqNJUaERGRCqpn81oseDiSAG83dhzLZtC0BPZm5podq8JSqREREanA2tbz46vRUTQK8OLomXPcOSOB5EOnzI5VIanUiIiIVHAhNT354pFIQoOrc+ZsEffO/pXl29LNjlXhVMhSc/vtt1OjRg3uvPPOCz5bunQpLVq0oFmzZrz33nsmpBMRESl/Nb3d+WxkBDe2rE1BsY3R85KZt+6Q2bEqlApZasaPH89HH310wfbi4mJiY2P56aef2LRpE5MnT+bkSU0nLSIiVYOnmwsz7+/E3eHB2Ozw9/it/Hv5Lux2zWUDFbTUXH/99fj4+FywPSkpiTZt2lCvXj28vb3p06cP33//vQkJRUREzOFitTBpUDsej24GwNSf9/LMF79RVKK5bK641KxatYp+/foRFBSEYRjEx8dfsE9cXBwNGzbEw8ODiIgIkpKSHJGVtLQ06tWrV/q+Xr16HD161CHfLSIiUlkYhsHj0c2ZNKgdFgMWJh9h5EcbyCsoNjuaqa641OTl5REaGkpcXNxFP1+wYAGxsbFMnDiRjRs3EhoaSu/evcnMzCzdJywsjLZt217wSktz3AJeBQUFZGdnn/cSERFxJvd0CWH2A53xcLWwctdx7pm9jhO5BWbHMo3Llf5Bnz596NOnz59+/uabbzJy5EiGDx8OwIwZM/jmm2/44IMPmDBhAgApKSlXFTYoKOi8KzNHjx6lS5cuF9130qRJvPTSS1d1HBERkcrixlaBfDayKyPmrOe3I1ncMT2BucO70DDAy+xo5c6hY2oKCwtJTk4mOjr6vwewWIiOjiYxMfGav79Lly5s3bqVo0ePkpuby7Jly+jdu/dF93322WfJysoqfaWmpl7z8UVERCqiDiE1+HJ0N4L9q3Ho5FnumJ7A5tQzZscqdw4tNSdOnKCkpITAwMDztgcGBpKefvnP00dHR3PXXXfx7bffUr9+/dJC5OLiwhtvvEGvXr0ICwvjySefpGbNmhf9Dnd3d3x9fc97iYiIOKvGtbz5cnQ32tbz5WReIXfPWsfPuzIv/YdO5IpvP5WHH3/88U8/69+/P/379y/HNCIiIpVDbR8P5o+KZPS8ZFbvOcFDczcwaVA7BncONjtauXDolZqAgACsVisZGRnnbc/IyKBOnTqOPJSIiIhchLe7C+/HhDOoQz1KbHae+eI33l2xp0rMZePQUuPm5kanTp1YsWJF6TabzcaKFSuIjIx05KFERETkT7i5WHhjcCiPXt8EgDd+2M3f47dSYnPuYnPFt59yc3PZu3dv6fsDBw6QkpKCv78/ISEhxMbGEhMTQ+fOnenSpQtvv/02eXl5pU9DiYiISNkzDINnbmlJoK8HLy7Zxie/HiYzp4Apd3egmpvV7HhlwrBf4fWolStX0qtXrwu2x8TEMGfOHACmTp3K5MmTSU9PJywsjClTphAREeGQwFcrOzsbPz8/srKyNGhYRESqlGVbjjF+QQqFxTY6hlTn/Zhwani5mR3rslzJ7/cVl5rKSqVGRESqsqQDp3ho7nqy84tpXMuLucO7EOzvaXasS7qS3+8KufaTiIiIOFaXRv58Mbobdf082H88j0HTE9iWlmV2LIdSqREREakimgf68NWj3WgR6MPxnAKGzFzH2r0nzI7lMCo1IiIiVUhdv2p8/kgkXRv7k1tQzLAPk1ic4hyLQ6vUiIiIVDF+1VyZO6ILfdvXpajEzvj5Kcxata/Sz2WjUiMiIlIFubtYeffuDoyIagTAq9/u5OWlO7BV4rlsVGpERESqKIvF4IV+rXn+1lYAfLD2AGPnbyK/qMTkZFdHpUZERKSKG9mzMe/cHYar1eCb344R80ESWeeKzI51xVRqREREhAFh9Zg7vAve7i78euAUg2ckcizrnNmxrohKjYiIiADQrWkAnz8cSW0fd3Zl5DBoWgK7M3LMjnXZVGpERESkVOsgX756tBtNanlxLCufO6cnkHTglNmxLotKjYiIiJynfg1PvhzdjU4NapCdX8x97//Kt1uOmR3rklRqRERE5ALVPd345KEIbm4dSGGxjTGfbmTO2gNmx/pLKjUiIiJyUR6uVqbf14n7uoZgt8OLS7bz2rKdFXYuG5UaERER+VNWi8HLA9rydO8WAMz4ZR9PLtxMYbHN5GQXUqkRERGRv2QYBmN6NWXyne2xWgwWbTrKg3PXk1tQbHa086jUiIiIyGW5q3Mw78V0xtPNyuo9JxgyM5HMnHyzY5VSqREREZHL1qtFbT4b2ZWaXm5sS8tm0LQE9h/PNTsWoFIjIiIiVyg0uDpfPdqNBjU9OXL6HHdMT2Dj4dNmx1KpERERkSvXoKYXX47uRvv6fpw+W8S9s9fx4/YMUzOp1IiIiMhVCfB257ORXbm+RS3yi2w8sSCFM2cLTcvjYtqRRUREpNLzcndh9gOdeWHxVm5qHUh1TzfTsqjUiIiIyDVxtVqYNKi92TF0+0lEREScg0qNiIiIOAWVGhEREXEKKjUiIiLiFFRqRERExCmo1IiIiIhTUKkRERERp6BSIyIiIk5BpUZEREScgkqNiIiIOAWVGhEREXEKKjUiIiLiFFRqRERExClUmVW67XY7ANnZ2SYnERERkcv1x+/2H7/jf6XKlJqcnBwAgoODTU4iIiIiVyonJwc/P7+/3MewX071cQI2m420tDR8fHwwDMOh352dnU1wcDCpqan4+vo69LsrAmc/P3D+c9T5VX7Ofo46v8qvrM7RbreTk5NDUFAQFstfj5qpMldqLBYL9evXL9Nj+Pr6Ou1/rOD85wfOf446v8rP2c9R51f5lcU5XuoKzR80UFhEREScgkqNiIiIOAWVGgdwd3dn4sSJuLu7mx2lTDj7+YHzn6POr/Jz9nPU+VV+FeEcq8xAYREREXFuulIjIiIiTkGlRkRERJyCSo2IiIg4BZUaERERcQoqNddg1apV9OvXj6CgIAzDID4+3uxIDjVp0iTCw8Px8fGhdu3aDBw4kF27dpkdy2GmT59O+/btSyeKioyMZNmyZWbHKjOvvfYahmHw+OOPmx3FYV588UUMwzjv1bJlS7NjOdTRo0e57777qFmzJtWqVaNdu3Zs2LDB7FgO07Bhwwv+HRqGwZgxY8yO5hAlJSX84x//oFGjRlSrVo0mTZrw8ssvX9Y6RpVFTk4Ojz/+OA0aNKBatWp069aN9evXm5KlyswoXBby8vIIDQ1lxIgRDBo0yOw4DvfLL78wZswYwsPDKS4u5rnnnuPmm29m+/bteHl5mR3vmtWvX5/XXnuNZs2aYbfbmTt3LgMGDGDTpk20adPG7HgOtX79embOnEn79u3NjuJwbdq04ccffyx97+LiPP9bO336NFFRUfTq1Ytly5ZRq1Yt9uzZQ40aNcyO5jDr16+npKSk9P3WrVu56aabuOuuu0xM5Tj/+te/mD59OnPnzqVNmzZs2LCB4cOH4+fnx7hx48yO5xAPPfQQW7du5eOPPyYoKIh58+YRHR3N9u3bqVevXvmGsYtDAPZFixaZHaNMZWZm2gH7L7/8YnaUMlOjRg37e++9Z3YMh8rJybE3a9bM/sMPP9ivu+46+/jx482O5DATJ060h4aGmh2jzPztb3+zd+/e3ewY5Wr8+PH2Jk2a2G02m9lRHKJv3772ESNGnLdt0KBB9qFDh5qUyLHOnj1rt1qt9qVLl563vWPHjvbnn3++3PPo9pNctqysLAD8/f1NTuJ4JSUlzJ8/n7y8PCIjI82O41Bjxoyhb9++REdHmx2lTOzZs4egoCAaN27M0KFDOXz4sNmRHObrr7+mc+fO3HXXXdSuXZsOHTowe/Zss2OVmcLCQubNm8eIESMcvvCwWbp168aKFSvYvXs3AJs3b2bNmjX06dPH5GSOUVxcTElJCR4eHudtr1atGmvWrCn3PM5znVbKlM1m4/HHHycqKoq2bduaHcdhtmzZQmRkJPn5+Xh7e7No0SJat25tdiyHmT9/Phs3bjTt/nZZi4iIYM6cObRo0YJjx47x0ksv0aNHD7Zu3YqPj4/Z8a7Z/v37mT59OrGxsTz33HOsX7+ecePG4ebmRkxMjNnxHC4+Pp4zZ84wbNgws6M4zIQJE8jOzqZly5ZYrVZKSkp45ZVXGDp0qNnRHMLHx4fIyEhefvllWrVqRWBgIJ999hmJiYk0bdq0/AOV+7UhJ4WT33565JFH7A0aNLCnpqaaHcWhCgoK7Hv27LFv2LDBPmHCBHtAQIB927ZtZsdyiMOHD9tr165t37x5c+k2Z7v99P+dPn3a7uvr6zS3EF1dXe2RkZHnbRs7dqy9a9euJiUqWzfffLP9tttuMzuGQ3322Wf2+vXr2z/77DP7b7/9Zv/oo4/s/v7+9jlz5pgdzWH27t1r79mzpx2wW61We3h4uH3o0KH2li1blnsWlRoHceZSM2bMGHv9+vXt+/fvNztKmbvxxhvto0aNMjuGQyxatKj0fzJ/vAC7YRh2q9VqLy4uNjtimejcubN9woQJZsdwiJCQEPuDDz543rZp06bZg4KCTEpUdg4ePGi3WCz2+Ph4s6M4VP369e1Tp049b9vLL79sb9GihUmJyk5ubq49LS3Nbrfb7YMHD7bfeuut5Z5BY2rkT9ntdh577DEWLVrETz/9RKNGjcyOVOZsNhsFBQVmx3CIG2+8kS1btpCSklL66ty5M0OHDiUlJQWr1Wp2RIfLzc1l37591K1b1+woDhEVFXXBNAq7d++mQYMGJiUqOx9++CG1a9emb9++ZkdxqLNnz2KxnP9Ta7VasdlsJiUqO15eXtStW5fTp0+zfPlyBgwYUO4ZNKbmGuTm5rJ3797S9wcOHCAlJQV/f39CQkJMTOYYY8aM4dNPP2Xx4sX4+PiQnp4OgJ+fH9WqVTM53bV79tln6dOnDyEhIeTk5PDpp5+ycuVKli9fbnY0h/Dx8blg/JOXlxc1a9Z0mnFRTz31FP369aNBgwakpaUxceJErFYr99xzj9nRHOKJJ56gW7duvPrqqwwePJikpCRmzZrFrFmzzI7mUDabjQ8//JCYmBineiQfoF+/frzyyiuEhITQpk0bNm3axJtvvsmIESPMjuYwy5cvx26306JFC/bu3cvTTz9Ny5YtGT58ePmHKfdrQ07k559/tgMXvGJiYsyO5hAXOzfA/uGHH5odzSFGjBhhb9Cggd3Nzc1eq1Yt+4033mj//vvvzY5VppxtTM2QIUPsdevWtbu5udnr1atnHzJkiH3v3r1mx3KoJUuW2Nu2bWt3d3e3t2zZ0j5r1iyzIznc8uXL7YB9165dZkdxuOzsbPv48ePtISEhdg8PD3vjxo3tzz//vL2goMDsaA6zYMECe+PGje1ubm72OnXq2MeMGWM/c+aMKVkMu92JpjUUERGRKktjakRERMQpqNSIiIiIU1CpEREREaegUiMiIiJOQaVGREREnIJKjYiIiDgFlRoRERFxCio1IiIi4hRUakRERMQpqNSIiIiIU1CpEREREaegUiMiIiJO4X8BBDeEHpTMo7IAAAAASUVORK5CYII=",
|
||
"text/plain": [
|
||
"<Figure size 640x480 with 1 Axes>"
|
||
]
|
||
},
|
||
"metadata": {},
|
||
"output_type": "display_data"
|
||
}
|
||
],
|
||
"source": [
|
||
"# Use Gaussian quadrature to evalute the integral of polynomials from x=-1 to +1\n",
|
||
"\n",
|
||
"import numpy as np\n",
|
||
"import scipy.integrate\n",
|
||
"import matplotlib.pyplot as plt\n",
|
||
"\n",
|
||
"# Some choices for functions to integrate:\n",
|
||
"# The first one is a random polynomial (the size parameter sets the degree)\n",
|
||
"#func = np.polynomial.Polynomial(np.random.randint(-10,high=10,size=20)) \n",
|
||
"#func = np.cos\n",
|
||
"#func = np.exp\n",
|
||
"func = lambda x: np.exp(-x**2)\n",
|
||
"\n",
|
||
"# Integrate this with quad to get the true value\n",
|
||
"I0, err0 = scipy.integrate.quad(func, -1.0, 1.0)\n",
|
||
"print('quad gives I = ', I0, err0)\n",
|
||
"\n",
|
||
"# then use Gaussian quadrature with different numbers of points\n",
|
||
"n_vec = np.array([])\n",
|
||
"err_vec = np.array([])\n",
|
||
"\n",
|
||
"for npoints in range(1,10):\n",
|
||
" I = 0.0\n",
|
||
" x, w = np.polynomial.legendre.leggauss(npoints)\n",
|
||
" for i in range(npoints):\n",
|
||
" I += w[i] * func(x[i])\n",
|
||
" err = (I-I0)/I0\n",
|
||
" print(npoints, 2*npoints-1, I, err)\n",
|
||
" n_vec = np.append(n_vec, npoints)\n",
|
||
" err_vec = np.append(err_vec, abs(err))\n",
|
||
" \n",
|
||
"plt.plot(n_vec, err_vec)\n",
|
||
"plt.yscale('log')\n",
|
||
"plt.xscale('linear')\n",
|
||
"plt.show()"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "22c1aadf",
|
||
"metadata": {},
|
||
"source": [
|
||
"- You can use the code above to check that the answer is exact for polynomials of degree $2n-1$ and smaller\n",
|
||
"- The scaling of the error with number of points depends on the function you are integrating, but decreases approximately exponentially (rather than power law with the Newton-Cotes methods)."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "1650be55",
|
||
"metadata": {},
|
||
"source": [
|
||
"Let's try the Gauss-Hermite quadrature which has weighting function $W(x)=e^{-x^2}$ for integrals from $-\\infty$ to $+\\infty$. In the next example I try $x^4 e^{-x^2}$ which becomes exact once we go to 3 terms:"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 3,
|
||
"id": "38475475",
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"quad gives I = 1.3293403881791366 1.5859180523983767e-08\n",
|
||
"1 1 0.0 -1.0\n",
|
||
"2 3 0.4431134627263788 -0.6666666666666666\n",
|
||
"3 5 1.3293403881791368 1.6703367090890604e-16\n",
|
||
"4 7 1.3293403881791377 8.351683545445302e-16\n",
|
||
"5 9 1.3293403881791368 1.6703367090890604e-16\n",
|
||
"6 11 1.3293403881791355 -8.351683545445302e-16\n",
|
||
"7 13 1.3293403881791372 5.011010127267181e-16\n",
|
||
"8 15 1.3293403881791364 -1.6703367090890604e-16\n",
|
||
"9 17 1.3293403881791372 5.011010127267181e-16\n"
|
||
]
|
||
},
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGdCAYAAADqsoKGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA8K0lEQVR4nO3dfXDU5333+8/uSloJaSWQhFYISciAJGywJYcHWdhuoFbCkV1uPxw7nI5zqkJDp2eM41TTTCG5C83JJKR3OwynY2qadGziJnZo7/s2SesEu1bsEAcaDETEqW2BQMYY0OoBpNWu0Era/Z0/pF0QEqCH3f3tw/s1s+PZB+3vu4DZD9d1fa/LYhiGIQAAgDhnNbsAAACAcCDUAACAhECoAQAACYFQAwAAEgKhBgAAJARCDQAASAiEGgAAkBAINQAAICGkmF1AtAQCAV28eFEOh0MWi8XscgAAwCQYhqG+vj4VFRXJar31WEzShJqLFy+qpKTE7DIAAMA0nD9/XsXFxbd8TdKEGofDIWnkFyU7O9vkagAAwGS43W6VlJSEvsdvJWlCTXDKKTs7m1ADAECcmczSERYKAwCAhECoAQAACYFQAwAAEgKhBgAAJARCDQAASAiEGgAAkBAINQAAICEQagAAQEIg1AAAgIRAqAEAAAmBUAMAABICoQYAACSEpDnQMlI+vdKvF945o5yM1LG3WWPvZ9lTJnUYFwAAmB5CzQx9euWqfvjrT277OpvVouz0lFDIyb4xBN1wy74uHDkIRAAA3BahZobm5aTryw+Vy311SL03uQ0OB+QPGLrSP6Qr/UNTvobVonEh6HahKPgahz1FViuBCACQ+Ag1M7QgL1ONn6u45WsGhvxjg07/+OBzs1DkGw4oYEg9/UPqmWYgcqTfYiToFjdHOoEIABA/CDVRkJ5qU3qqTc7s9Cn/7MCQ/5ajQLcKRQNDI4EoeH+qLBbJYU/R7FlpunOeQ3//h/fKnmKb8vsAABANhJoYFwxEBdMNRAMjgadngtGh2wUiw5DcA8NyDwzrk8v9Otp2WQ+Wz43ApwQAYOYINQksFIgcUw9EvmF/KPB8498+0C9Pd+mUy0OoAQDELPapwYTsKSNhaHGBQ/eWzJYknXb1mVsUAAC3QKjBbVUUOiRJpwg1AIAYRqjBbVU4R0LNaZdHhmGYXA0AABMj1OC2yvIylWK1qM83rEu9A2aXAwDAhAg1uK20FKsWzs2UJLUwBQUAiFGEGkxKeWgKilADAIhNhBpMSuVoqGlp95hcCQAAEyPUYFIqnFmSpNMdjNQAAGIToQaTcn0HVCBABxQAIPYQajApC/IylZZi1dUhvz69ctXscgAAGCeuQs3jjz+uOXPm6MknnzS7lKRjs1q0aO7IFBSb8AEAYlFchZrnnntOL7/8stllJK3K0XU1tHUDAGJRXIWaNWvWyOFwmF1G0qKtGwAQy8IWag4dOqT169erqKhIFotFBw4cGPeaPXv2qKysTOnp6aqpqdHRo0fDdXlEQait20VbNwAg9oQt1Hi9XlVVVWnPnj0TPr9//341NjZqx44dOnHihKqqqrRu3Tp1dHSEXlNdXa1ly5aNu128eDFcZWIGgh1QZzo98tMBBQCIMSnheqP6+nrV19ff9Pldu3Zp8+bN2rhxoyRp7969ev311/Xiiy9q69atkqTm5uZwlSOfzyefzxe673a7w/beyap4ToYyUm26OuTXuW6vFo4uHAYAIBZEZU3N4OCgjh8/rrq6umsXtlpVV1enI0eOROSaO3fuVE5OTuhWUlISkeskE6vVonInHVAAgNgUlVDT1dUlv98vp9M55nGn06n29vZJv09dXZ2eeuop/fSnP1VxcfEtA9G2bdvU29sbup0/f37a9eOa8oKRKahTrKsBAMSYsE0/RcNbb7016dfa7XbZ7fYIVpOcKgtp6wYAxKaojNTk5+fLZrPJ5XKNedzlcqmwsDAaJSBMaOsGAMSqqISatLQ0LV++XE1NTaHHAoGAmpqaVFtbG40SECbBtu6znV4NDgdMrgYAgGvCNv3k8XjU2toaut/W1qbm5mbl5uaqtLRUjY2Namho0IoVK7Rq1Srt3r1bXq831A2F+DAvJ10Oe4r6fMP6uNsbavMGAMBsYQs1x44d09q1a0P3GxsbJUkNDQ3at2+fNmzYoM7OTm3fvl3t7e2qrq7WwYMHxy0eRmyzWEY6oE580qNTrj5CDQAgZoQt1KxZs0aGcesN2bZs2aItW7aE65IwSYXTMRJq2vuke8yuBgCAEXF19hNiQ3CxMG3dAIBYQqjBlFWGQg0dUACA2EGowZRVjO4q/HG3VwNDfpOrAQBgBKEGUzbXYdfsWakKGCOHWwIAEAsINZgyi8WiioLgJnyEGgBAbCDUYFoqCjnYEgAQWwg1mJYKFgsDAGIMoQbTwmndAIBYQ6jBtAQ7oD653K/+wWGTqwEAgFCDacrLsis/K02S1NrBaA0AwHyEGkxbcF1NSzvragAA5iPUYNqCoeY0IzUAgBhAqMG00QEFAIglhBpMW3Cx8CmmnwAAMYBQg2kLntZ9sXdAfQNDJlcDAEh2hBpMW05Gqgqz0yWxXw0AwHyEGsxI+egU1GnW1QAATEaowYxUBtu6CTUAAJMRajAjobZupp8AACYj1GBGKgoZqQEAxAZCDWakvGBkTU1nn089/YMmVwMASGaEGsxIpj1F82dnSKIDCgBgLkINZqySKSgAQAwg1GDGaOsGAMQCQg1mrJLTugEAMYBQgxm7/mBLwzBMrgYAkKwINZixxQVZslikK/1D6vLQAQUAMAehBjOWnmrTgtxZklhXAwAwT9yEmp6eHq1YsULV1dVatmyZvve975ldEq5Tft0UFAAAZkgxu4DJcjgcOnTokGbNmiWv16tly5bpiSeeUF5entmlQSOLhf/jA5da2KsGAGCSuBmpsdlsmjVrZIrD5/PJMAwWpcYQ2roBAGYLW6g5dOiQ1q9fr6KiIlksFh04cGDca/bs2aOysjKlp6erpqZGR48endI1enp6VFVVpeLiYn31q19Vfn5+mKrHTF2/AR9hEwBghrCFGq/Xq6qqKu3Zs2fC5/fv36/Gxkbt2LFDJ06cUFVVldatW6eOjo7Qa4LrZW68Xbx4UZI0e/ZsnTx5Um1tbXrllVfkcrnCVT5m6I78TNmsFvUNDMvl9pldDgAgCYVtTU19fb3q6+tv+vyuXbu0efNmbdy4UZK0d+9evf7663rxxRe1detWSVJzc/OkruV0OlVVVaVf/vKXevLJJyd8jc/nk8937cvV7XZP8pNgOuwpNt2Rn6nWDo9aXH0qzEk3uyQAQJKJypqawcFBHT9+XHV1ddcubLWqrq5OR44cmdR7uFwu9fWNrNfo7e3VoUOHVFlZedPX79y5Uzk5OaFbSUnJzD4EbquCdTUAABNFJdR0dXXJ7/fL6XSOedzpdKq9vX1S73Hu3Dk9+OCDqqqq0oMPPqhnn31Wd999901fv23bNvX29oZu58+fn9FnwO2VF9DWDQAwT9y0dK9atWrS01OSZLfbZbfbI1cQxrm2WJi2bgBA9EVlpCY/P182m23cwl6Xy6XCwsJolIAoCE4/tbr6FAjQAQUAiK6ohJq0tDQtX75cTU1NoccCgYCamppUW1sbjRIQBQvyMpVms8o76NeFnqtmlwMASDJhm37yeDxqbW0N3W9ra1Nzc7Nyc3NVWlqqxsZGNTQ0aMWKFVq1apV2794tr9cb6oZC/Eu1WbVwbqY+au/T6Y4+lYyeBwUAQDSELdQcO3ZMa9euDd1vbGyUJDU0NGjfvn3asGGDOjs7tX37drW3t6u6uloHDx4ct3gY8a3C6dBH7X1qaffo95fwewsAiJ6whZo1a9bcdifZLVu2aMuWLeG6JGIQbd0AALPEzdlPiA+h07o7CDUAgOgi1CCsKkdDzWmXR346oAAAUUSoQViV5M6SPcUq33BA5y/3m10OACCJEGoQVjarReWj62paWFcDAIgiQg3CrqIgOAVFqAEARA+hBmFXwXEJAAATEGoQdrR1AwDMQKhB2AVP6z7b6dWQP2ByNQCAZEGoQdjNn52hzDSbBv0Bnev2ml0OACBJEGoQdlarRYuDm/CxrgYAECWEGkREZbCtu511NQCA6CDUICIqgjsLc1wCACBKCDWIiGCoYaQGABAthBpERDDUfNzdL9+w3+RqAADJgFCDiHBm2+VIT5E/YKitiw4oAEDkEWoQERaLJXRiN1NQAIBoINQgYsqDi4Vp6wYARAGhBhFTyWndAIAoItQgYkJt3YQaAEAUEGoQMcHTus9d7tfVQTqgAACRRahBxORn2ZWbmSbDkM50sq4GABBZhBpEVHnByLqaU0xBAQAijFCDiKocnYJisTAAINIINYgo2roBANFCqEFEsQEfACBaCDWIqIrRvWou9FyVxzdscjUAgERGqEFEzZ6VpgKHXRL71QAAIotQg4irYF0NACAKCDWIuHInbd0AgMhLMbuAqSgrK1N2drasVqvmzJmjt99+2+ySMAmhxcKEGgBABMVVqJGkw4cPKysry+wyMAW0dQMAooHpJ0RcsAOq3T2g3qtDJlcDAEhUYQs1hw4d0vr161VUVCSLxaIDBw6Me82ePXtUVlam9PR01dTU6OjRo1O6hsVi0Wc/+1mtXLlSP/zhD8NUOSLNkZ6qopx0SXRAAQAiJ2zTT16vV1VVVdq0aZOeeOKJcc/v379fjY2N2rt3r2pqarR7926tW7dOLS0tKigokCRVV1dreHj8XiZvvvmmioqK9O6772r+/Pm6dOmS6urqdPfdd+uee+4J10dABFUUOnSxd0Atrj6tKMs1uxwAQAIKW6ipr69XfX39TZ/ftWuXNm/erI0bN0qS9u7dq9dff10vvviitm7dKklqbm6+5TXmz58vSZo3b54efvhhnThx4qahxufzyefzhe673e6pfByEWYXToXdaOllXAwCImKisqRkcHNTx48dVV1d37cJWq+rq6nTkyJFJvYfX61Vf38jUhcfj0c9//nMtXbr0pq/fuXOncnJyQreSkpKZfQjMCKd1AwAiLSqhpqurS36/X06nc8zjTqdT7e3tk3oPl8ulBx54QFVVVbrvvvv0R3/0R1q5cuVNX79t2zb19vaGbufPn5/RZ8DMBE/rJtQAACIlblq6Fy5cqJMnT0769Xa7XXa7PYIVYSoWj47UdHkG1e3xKS+L3xsAQHhFZaQmPz9fNptNLpdrzOMul0uFhYXRKAEmm5WWotLcWZKkU6yrAQBEQFRCTVpampYvX66mpqbQY4FAQE1NTaqtrY1GCYgBwf1qTncwBQUACL+wTT95PB61traG7re1tam5uVm5ubkqLS1VY2OjGhoatGLFCq1atUq7d++W1+sNdUMh8VU4HXrrww61tBNqAADhF7ZQc+zYMa1duzZ0v7GxUZLU0NCgffv2acOGDers7NT27dvV3t6u6upqHTx4cNziYSQuTusGAESSxTAMw+wiosHtdisnJ0e9vb3Kzs42u5yk9F8Xe/XI37+rnIxUNW//nCwWi9klAQBi3FS+vzn7CVGzaG6WrBap9+qQOvt8t/8BAACmgFCDqElPtaksL1MSHVAAgPAj1CCqgutqWtiEDwAQZoQaRFWorZtQAwAIM0INoqqikJEaAEBkEGoQVde3dSdJ4x0AIEoINYiqsrxMpVgt8viGdbF3wOxyAAAJhFCDqEpLsWrh3GAHFFNQAIDwIdQg6spDU1CEGgBA+BBqEHWVwbbudvaqAQCED6EGUcdp3QCASCDUIOqu74AKBOiAAgCEB6EGUbcgL1NpKVZdHfLr0ytXzS4HAJAgCDWIOpvVokVzR6ag2IQPABAuhBqYonJ0XQ1t3QCAcCHUwBS0dQMAwo1QA1OE2rpdtHUDAMKDUANTBDugznR6NOwPmFwNACAREGpgiuI5GcpItWlwOKBzl/vNLgcAkAAINTCF1WpReXATPtbVAADCgFAD05QXcFwCACB8CDUwTWXhaFs3xyUAAMKAUAPT0NYNAAgnQg1ME2zrPtvp1eAwHVAAgJkh1MA083LS5bCnaDhg6ONur9nlAADiHKEGprFYrnVAtbQzBQUAmBlCDUxVwboaAECYEGpgqvLQcQmEGgDAzBBqYKrK0EgNe9UAAGYmbkJNS0uLqqurQ7eMjAwdOHDA7LIwQxWja2o+7vZqYMhvcjUAgHiWYnYBk1VZWanm5mZJksfjUVlZmT73uc+ZWxRmbK7DrtmzUtXTP6QznR4tLcoxuyQAQJyKm5Ga6/3kJz/RQw89pMzMTLNLwQxZLBZVFDAFBQCYubCFmkOHDmn9+vUqKiqSxWKZcGpoz549KisrU3p6umpqanT06NFpXetf/uVftGHDhhlWjFhRMXpcAouFAQAzEbZQ4/V6VVVVpT179kz4/P79+9XY2KgdO3boxIkTqqqq0rp169TR0RF6TXV1tZYtWzbudvHixdBr3G63Dh8+rIcffjhcpcNktHUDAMIhbGtq6uvrVV9ff9Pnd+3apc2bN2vjxo2SpL179+r111/Xiy++qK1bt0pSaM3Mrfz4xz/W5z//eaWnp9/ydT6fTz6fL3Tf7XZP4lPADKHTugk1AIAZiMqamsHBQR0/flx1dXXXLmy1qq6uTkeOHJnSe0126mnnzp3KyckJ3UpKSqZcN6Ij2AF1/vJV9Q8Om1wNACBeRSXUdHV1ye/3y+l0jnnc6XSqvb190u/T29uro0ePat26dbd97bZt29Tb2xu6nT9/fsp1IzrysuzKz0qTJLV2sFgYADA9cdPSLUk5OTlyuVyTeq3dbpfdbo9wRQiXCqdDXZ5utbT36Z7i2WaXAwCIQ1EZqcnPz5fNZhsXSFwulwoLC6NRAmJcaLEwIzUAgGmKSqhJS0vT8uXL1dTUFHosEAioqalJtbW10SgBMS4YajitGwAwXWGbfvJ4PGptbQ3db2trU3Nzs3Jzc1VaWqrGxkY1NDRoxYoVWrVqlXbv3i2v1xvqhkJyCy4Wpq0bADBdYQs1x44d09q1a0P3GxsbJUkNDQ3at2+fNmzYoM7OTm3fvl3t7e2qrq7WwYMHxy0eRnIKntZ9sXdA7oEhZaenmlwRACDeWAzDMMwuIhrcbrdycnLU29ur7Oxss8vBBO77dpPa3QP6X//Pai1fMMfscgAAMWAq399xefYTElM5U1AAgBkg1CBmVDrZWRgAMH2EGsSMa2dA0dYNAJg6Qg1iRkUhIzUAgOkj1CBmlBeMrKnp7PPpinfQ5GoAAPGGUIOYkWlP0fzZGZKkU4zWAACmiFCDmFI5OgV1iuMSAABTRKhBTKGtGwAwXYQaxJRKzoACAEwToQYxJdjWfcrVpyTZ7BoAECaEGsSUxQVZslikK/1D6vLQAQUAmDxCDWJKeqpNC3JnSWJdDQBgagg1iDnlHJcAAJgGQg1iTmVoXQ1t3QCAySPUIObQ1g0AmA5CDWJO5XVnQNEBBQCYLEINYs4d+ZmyWS3qGxiWy+0zuxwAQJwg1CDm2FNsuiM/UxKLhQEAk0eoQUyqYF0NAGCKCDWISeUFHJcAAJgaQg1iEqd1AwCmilCDmHT99FMgQAcUAOD2CDWISQvyMpVms6p/0K8LPVfNLgcAEAcINYhJqTarFs4d6YA63cG6GgDA7RFqELMqgmdAtbOuBgBwe4QaxCzaugEAU0GoQczitG4AwFQQahCzgqd1t3Z45KcDCgBwG3EVav7u7/5OS5cu1bJly/SDH/zA7HIQYSW5s2RPsco3HNAnl/vNLgcAEOPiJtS8//77euWVV3T8+HG99957ev7559XT02N2WYggm9Wi8tF1NaeYggIA3EbchJoPP/xQtbW1Sk9PV0ZGhqqqqnTw4EGzy0KEVYwel8BiYQDA7YQt1Bw6dEjr169XUVGRLBaLDhw4MO41e/bsUVlZmdLT01VTU6OjR49O+v2XLVumd955Rz09Pbpy5YreeecdXbhwIVzlI0ZVFAYXC9PWDQC4tZRwvZHX61VVVZU2bdqkJ554Ytzz+/fvV2Njo/bu3auamhrt3r1b69atU0tLiwoKCiRJ1dXVGh4eHvezb775pu666y59+ctf1u///u8rJydH9913n2w2W7jKR4yirRsAMFkWwzDC3lZisVj02muv6bHHHgs9VlNTo5UrV+r555+XJAUCAZWUlOjZZ5/V1q1bp3yNL33pS3r88cf1yCOPTPi8z+eTz+cL3Xe73SopKVFvb6+ys7OnfD2Y4/zlfj34P95Wqs2iD/7f/0OptriZMQUAhIHb7VZOTs6kvr+j8g0xODio48ePq66u7tqFrVbV1dXpyJEjk36fjo4OSVJLS4uOHj2qdevW3fS1O3fuVE5OTuhWUlIy/Q8A08yfnaHMNJuG/IbOdXvNLgcAEMOiEmq6urrk9/vldDrHPO50OtXe3j7p93n00Ud111136Ytf/KJeeuklpaTcfPZs27Zt6u3tDd3Onz8/7fphHqvVosUclwAAmISwramJhqmM6tjtdtnt9ghWg2ipdGbp5PkenXL16RHNM7scAECMispITX5+vmw2m1wu15jHXS6XCgsLo1EC4ljwYEtO6wYA3EpUQk1aWpqWL1+upqam0GOBQEBNTU2qra2NRgmIY9dO6ybUAABuLmzTTx6PR62traH7bW1tam5uVm5urkpLS9XY2KiGhgatWLFCq1at0u7du+X1erVx48ZwlYAEFQw1H3f3yzfslz2FVn4AwHhhCzXHjh3T2rVrQ/cbGxslSQ0NDdq3b582bNigzs5Obd++Xe3t7aqurtbBgwfHLR4GbuTMtsuRnqK+gWGd7fTqznm05AMAxovIPjWxaCp97og9T75wWMfOXdH/939V69Hq+WaXAwCIkpjbpwaYqfLRKSgOtgQA3AyhBnGhMnRaN3vVAAAmRqhBXAi1dTNSAwC4CUIN4kLwtO5zl/t1ddBvcjUAgFhEqEFcyM+yKzczTYYhnelkCgoAMB6hBnGjvGBkXQ2b8AEAJkKoQdyoHJ2COsVxCQCACRBqEDdCbd2M1AAAJkCoQdyoDO1Vw5oaAMB4hBrEjYrRvWou9FyVxzdscjUAgFhDqEHcmD0rTQUOuyT2qwEAjEeoQVy5tgkfU1AAgLEINYgr5aNTUC2M1AAAbkCoQVyp5GBLAMBNEGoQVzitGwBwM4QaxJVgB5TL7VPv1SGTqwEAxBJCDeKKIz1VRTnpkuiAAgCMRahB3Ame2M1iYQDA9Qg1iDu0dQMAJkKoQdzhtG4AwEQINYg7wdO6T3NaNwDgOoQaxJ3FoyM1XZ5BdXt8JlcDAIgVhBrEnVlpKSrNnSWJE7sBANcQahCXgvvVMAUFAAgi1CAuBTugWCwMAAgi1CAu0dYNALgRoQZx6frTug3DMLkaAEAsINQgLi2amyWrReq9OqTOPjqgAACEGsSp9FSbyvIyJXFcAgBgREyGmscff1xz5szRk08+OaXnkFyC62po6wYASDEaap577jm9/PLLU34OySXU1s1IDQBAMRpq1qxZI4fDMeXnkFw4rRsAcL0ph5pDhw5p/fr1KioqksVi0YEDB8a9Zs+ePSorK1N6erpqamp09OjRcNQKjHF9WzcdUACAKYcar9erqqoq7dmzZ8Ln9+/fr8bGRu3YsUMnTpxQVVWV1q1bp46OjtBrqqurtWzZsnG3ixcvTv+TIOmU5WUqxWqRxzesi70DZpcDADBZylR/oL6+XvX19Td9fteuXdq8ebM2btwoSdq7d69ef/11vfjii9q6daskqbm5eXrVToHP55PPd63V1+12R/yaiK60FKsWzs3UKZdHp1x9mj87w+ySAAAmCuuamsHBQR0/flx1dXXXLmC1qq6uTkeOHAnnpW5r586dysnJCd1KSkqien1ER3mwA4rjEgAg6YU11HR1dcnv98vpdI553Ol0qr29fdLvU1dXp6eeeko//elPVVxcPCYQ3eq5623btk29vb2h2/nz56f3oRDTKmnrBgCMmvL0UzS89dZb03ruena7XXa7PVwlIUZxWjcAICisIzX5+fmy2WxyuVxjHne5XCosLAznpQBJYzugAgE6oAAgmYU11KSlpWn58uVqamoKPRYIBNTU1KTa2tpwXgqQJC3Iy1RailVXh/z69MpVs8sBAJhoytNPHo9Hra2tofttbW1qbm5Wbm6uSktL1djYqIaGBq1YsUKrVq3S7t275fV6Q91QQDjZrBYtmpulDy+51eLqU2neLLNLAgCYZMqh5tixY1q7dm3ofmNjoySpoaFB+/bt04YNG9TZ2ant27ervb1d1dXVOnjw4LjFw0C4VDpHQs0pV58+dxd/zgAgWU051KxZs+a2u7du2bJFW7ZsmXZRwFSE2ro5LgEAklpMnv0ETAVt3QAAiVCDBBDsgDrT6dGwP2ByNQAAsxBqEPeK52QoI9WmweGAzl3uN7scAIBJCDWIe1arReXBTfhYVwMASYtQg4RQXjAyBdXSzroaAEhWhBokhMrCkZGaUxyXAABJi1CDhMBp3QAAQg0SQrCtu63Lq8FhOqAAIBkRapAQ5uWky2FP0XDAUFuX1+xyAAAmINQgIVgs1zqg2FkYAJIToQYJI7gJH23dAJCcCDVIGMHFwi2EGgBISoQaJIzK0EgNe9UAQDIi1CBhVIyuqfm426uBIb/J1QAAoo1Qg4Qx12HX7FmpChgjh1sCAJILoQYJw2KxqGL0uAQ6oAAg+RBqkFAqgsclsK4GAJIOoQYJhbZuAEhehBoklNBp3YQaAEg6hBoklGAH1PnLV9U/OGxyNQCAaCLUIKHkZdmVn5Umif1qACDZEGqQcILrauiAAoDkQqhBwiHUAEByItQg4VwLNUw/AUAyIdQg4QQXC9PWDQDJhVCDhBM8rfti74DcA0MmVwMAiBZCDRJOTkaqCrPTJdEBBQDJhFCDhFTuDB6XwBQUACSLmAw1jz/+uObMmaMnn3xyzOM9PT1asWKFqqurtWzZMn3ve98zqULEuko6oAAg6cRkqHnuuef08ssvj3vc4XDo0KFDam5u1q9//Wt9+9vfVnd3twkVItbR1g0AyScmQ82aNWvkcDjGPW6z2TRr1ixJks/nk2EYMgwj2uUhDlQU0tYNAMlmyqHm0KFDWr9+vYqKimSxWHTgwIFxr9mzZ4/KysqUnp6umpoaHT16NBy1ShqZgqqqqlJxcbG++tWvKj8/P2zvjcRRXjCypqazz6cr3kGTqwEARMOUQ43X61VVVZX27Nkz4fP79+9XY2OjduzYoRMnTqiqqkrr1q1TR0dH6DXBNTE33i5evHjb68+ePVsnT55UW1ubXnnlFblcrql+BCSBTHuK5s/OkMQUFAAki5Sp/kB9fb3q6+tv+vyuXbu0efNmbdy4UZK0d+9evf7663rxxRe1detWSVJzc/P0qr2O0+lUVVWVfvnLX45bUCyNTE/5fL7QfbfbPeNrIr5UFjp0oeeqTnV4VLMwz+xyAAARFtY1NYODgzp+/Ljq6uquXcBqVV1dnY4cOTLj93e5XOrrG/lXd29vrw4dOqTKysoJX7tz507l5OSEbiUlJTO+PuJLqK27nZEaAEgGYQ01XV1d8vv9cjqdYx53Op1qb2+f9PvU1dXpqaee0k9/+lMVFxeHAtG5c+f04IMPqqqqSg8++KCeffZZ3X333RO+x7Zt29Tb2xu6nT9/fvofDHGJtm4ASC5Tnn6KhrfeemvCx1etWjXpqSu73S673R7GqhBvrm/rNgxDFovF5IoAAJEU1pGa/Px82Wy2cYt3XS6XCgsLw3kp4LYWF2TJYpGu9A+py0MHFAAkurCGmrS0NC1fvlxNTU2hxwKBgJqamlRbWxvOSwG3lZ5q04LckX2NOLEbABLflKefPB6PWltbQ/fb2trU3Nys3NxclZaWqrGxUQ0NDVqxYoVWrVql3bt3y+v1hrqhgGgqdzr0cXe/Wlx9Wr2YPY0AIJFNOdQcO3ZMa9euDd1vbGyUJDU0NGjfvn3asGGDOjs7tX37drW3t6u6uloHDx4ct3gYiIZKp0P/8YGLnYUBIAlMOdSsWbPmtkcTbNmyRVu2bJl2UUC4cFo3ACSPmDz7CQiXysKxHVAAgMRFqEFCuyM/UzarRX0Dw2p3D5hdDgAgggg1SGj2FJvuyM+UxIndAJDoCDVIeBWj62po6waAxEaoQcIrLxhZV9PCGVAAkNAINUh4ocXCHUw/AUAiI9Qg4V0//RQI0AEFAImKUIOEtyAvU2k2q/oH/brQc9XscgAAEUKoQcJLtVm1cG6wA4p1NQCQqAg1SAoVzuAmfKyrAYBERahBUqCtGwASH6EGSaF8dKSmhVADAAmLUIOkUDkaalo7PPLTAQUACYlQg6RQkjtL9hSrfMMBfXK53+xyAAARQKhBUrBZLSofXVdDBxQAJCZCDZJGxehxCac4LgEAEhKhBkmjguMSACChEWqQNGjrBoDERqhB0gie1n2m06Mhf8DkagAA4UaoQdKYPztDmWk2DfkNnev2ml0OACDMCDVIGlarRYuDm/C1s64GABINoQZJpZK2bgBIWIQaJJVrB1sSagAg0RBqkFQINQCQuAg1SCrBUPNxd798w36TqwEAhBOhBknFmW2XIz1F/oChs510QAFAIiHUIKlYLJbQid1MQQFAYiHUIOmUE2oAICERapB0rrV1s1cNACSSmAw1jz/+uObMmaMnn3xy3HNlZWW65557VF1drbVr15pQHeIdHVAAkJhiMtQ899xzevnll2/6/OHDh9Xc3Ky33347ilUhUQRP6/7kcr+uDtIBBQCJIiZDzZo1a+RwOMwuAwkqP8uu3Mw0GYbU2sEUFAAkiimHmkOHDmn9+vUqKiqSxWLRgQMHxr1mz549KisrU3p6umpqanT06NFw1CpppHvls5/9rFauXKkf/vCHYXtfJJfyAo5LAIBEkzLVH/B6vaqqqtKmTZv0xBNPjHt+//79amxs1N69e1VTU6Pdu3dr3bp1amlpUUFBgSSpurpaw8PD4372zTffVFFR0S2v/+6772r+/Pm6dOmS6urqdPfdd+uee+6Z6sdAkqssdOjXbZd1qoNQAwCJYsqhpr6+XvX19Td9fteuXdq8ebM2btwoSdq7d69ef/11vfjii9q6daskqbm5eXrVSpo/f74kad68eXr44Yd14sSJCUONz+eTz+cL3Xe73dO+JhJPqK27nVADAIliyqHmVgYHB3X8+HFt27Yt9JjValVdXZ2OHDky4/f3er0KBAJyOBzyeDz6+c9/ri984QsTvnbnzp36xje+MeNrIjFd24Av8dfUfHqlX0fbLisj1ab7FuZpTmaa2SVhCgIBQx9ccus353tkt1mVl5Wm3My00NqwWWk2WSwWs8sEYkJYQ01XV5f8fr+cTueYx51Opz766KNJv09dXZ1Onjwpr9er4uJi/eu//qtqa2vlcrn0+OOPS5L8fr82b96slStXTvge27ZtU2NjY+i+2+1WSUnJND4VElHF6F41F3quyuMbVpY9rP8rmOqKd1BHznbrV61d+lVrlz7u7g89Z7FIy4pytHpxnh5YnK+VZblKT7WZWC0m8kl3v94d/f07fKZLV/qHbvra9FSr8jJHAs6NgSdv9LHrn5+Vljh/1oEbxeSf7rfeemvCxxcuXKiTJ09O6j3sdrvsdns4y0ICmT0rTQUOuzr6fDrt6tO9pXPMLmnarg76dezcZb3b2qXDrd363cVeGca1521Wi+4pzpHXN6xTLo/ev9Cr9y/06h9/cVZpKVYtL52jB8rztXpRnu6en6MUW0w2RSa0Lo9Ph89063Brl95t7dKnV66OeT4zzaYVZbmyWKTL3kF1ewbV5fHJNxzQwFBAF3qu6kLP1Zu8+1gZqbbR4JM2GnTsofCTm2kfDUHXwhGhF/EkrKEmPz9fNptNLpdrzOMul0uFhYXhvBQwYxVOhzr6fDoVZ6Fm2B/Q+xd6dfhMt9493aXj565o0B8Y85oKZ5buX5yv+xflq2ZhrhzpqZKkDvfAyM+1dulwa5cu9g7oyNluHTnbLUlypKeodmHeyM8uzteiuZlMbUSA1zesox9f1q9Od+lXZ7r14aWxa/5SbRbdWzpH9y/K1wPlebqneLZSbwibhmGof9Cvy96RgHPZO6ju0cBz2etTt2f0vteny55BdXkHNTgc0NUh/5RC0Kw0Wyjw5GeODUK5140EBUeJCEEwU1hDTVpampYvX66mpiY99thjkqRAIKCmpiZt2bIlnJcCZqzcmaV3W7tifl2NYRg60+nVr0b/Ff+fZ7vVNzC2e7AoJz0URFYvylNBdvqE71WQna7H7p2vx+6dL8Mw1Nbl1a/OdOtXp0emOdwDw3rzA5fe/GDkHybObLvuX5yvB0bf23mT98WtDfkDOnm+R79qHZkWPPHJFQ0HjDGvuXNeth5YnKfVi/O1qixXmbeZErVYLMq0pyjTnqKS3Fm3rcEwDHkH/er2+MaEny7PoC57B8eGo9HHBv0B9Q/61X/5qs5fnlwIyrKnXBd2Rqe+stLGjgZdN1VmTyEEIXymHGo8Ho9aW1tD99va2tTc3Kzc3FyVlpaqsbFRDQ0NWrFihVatWqXdu3fL6/WGuqGAWBHLp3W39w6MrIk5M7KuwuX2jXk+Oz1Fqxfl6/7yfN2/KE935E99RMVisWjh3CwtnJul//u+BfIHDP3Xxd7QWo73Pr4il9un/33igv73iQuSpMUFWbp/0chIzn2L8pQ9OgKEsQzD0CmXJ/Rr+euz3fLesHt18ZwMPVier9WLRoJoXlZkp8stFouy7CnKsqdoQV7mbV9vGIb6fMO67BkZ7QmO/NwYfkYeG3l+OGDI4xuWxzesTy733/YakuSwp4QCTvGcWVoyz6E7C7O1ZJ5DhdnpjBTGsCF/QG1dXn14ya2P2vv00SW3KpwObXv4TtNqshiGYdz+Zde88847E5651NDQoH379kmSnn/+ef3t3/6t2tvbVV1drb//+79XTU1NWAqeLrfbrZycHPX29io7O9vUWhAbjp+7ov/zhcNyZtv166/VmVqLe2BI/3lmdHHvme5xOx2npVi1smxOaNRkaVGObNbI/mU/MOTX8XNXQguOf3th7Fodq0W6p3h2aBTnMwtmJ/W/ui/0XA39Wv2qtVtdnrFBdM6sVK0Ojnotyldp3u1HV+KJYRhyDwyrezTwBEeAQiND3rHTYpe9g/IHbv31M3tWqpYUOrSkMFt3zhv5b4XToYy05P1zZpYuj08fXerTR+1ufTj639Muz7ip72Xzs/Xvzz4Y1mtP5ft7yqEmXhFqcKO+gSHd/ddvSpJObv+8cmZFb9TBNzwSGA63jqxv+e2nPQrcEBjunp8TmlJavmCO6WsVevuHrnVVnenS2U7vmOfTU61aWZYbCl53zcuWNcLBy0w9/YM6ElyfdKZbbV1jfz0yUm1aeUeuHlg8MrJ1Z2Fi/3pMVSBgyD0wFJoK6/b41NbtDX1xnun0Thh6LBbpjrxMLZkXDDvZWlLoUPGcDEZ1wsA37NeZDq8+ah8Zffnw0kiIuTGkB2Wm2VRZ6Bj5fZiXraVF2fpMmNcoEmomQKjBRFbvbNLF3gH965/VamVZbsSuE9xr5NrUzmUNDI39F87CuZm6f9FIiKldmBfVkDUdF0dHJoILjzv7xo9M1I5OVd2/KF8L8mbF9ZfOwJBf7318ObQuZqIus6riHD2wOF+rF+fr3tLkHrmaKd+wX60dnpFRgUvXvmC7vYMTvt5hT1FloeO6sONQZWF2Qm3XEE6GYaijzxeaOvrwklsfXerTmU7PuPVe0kiYLMvLDI2cBacJi+dkRDysE2omQKjBRP74paN6p6VT33p8mZ6uWRC29zUMQ+e6+0NrYg6f6VbPDXuNzHXYR74AR7/4i2ZnhO360WYYhlo7rq0h+c+zl+XxjV3MPH92xugXfp5WL8rXXEdsb7ngDxh6/0JvaErp2LkrGhwe32W2etHIyNT1XWaInM4+38gowqXRUYT2PrV29GnIP/FXWWnurJEv4nnZunP0vwtyZyXVqNnAkF+nXZ7RXy93aDTsZvsfZaenjPn1WlLoUIXTcdvF65FCqJkAoQYT+fZPP9R3D51VQ+0CfePRZTN6r84+nw6fubam4saW2Sx7iu5beG16ZnFBVlyPXNzKsD+gk59eCwQnPrky7ktnSaEj9Gux6o7bd/tEWrDL7PCZLr17eqTLzH1Dl9m8UJdZnu5flH/TLjNE15A/oLOd3jHrPT685B63wD4oI9WmikKH7hod1QmOPsT66OjtGIahCz1Xr619GV2829bl1UTLl6wWaeHcLC0JTh+NhpiinNhaoE2omQChBhP512Pn9dX/+VvVLszTq39635R+1uMb1tG27tB0xEc3nCOVarPoM6VzQutiqoqTd2O7/sFhHW27HNpb54Mb9mVJsVp0b+ns0K9Vdcn4fVkiweUeCIXQX7V2qd09MOb57PQU1S7KC00pLZxGlxnMc9k7GBrVCa4RaWnvk++GEbegopz0kRGK66awyvIyY/L/W69vWC2uvmuf7VKfPmx3j9vuIWjOrNTR4DIydXTXvGwtLsgyfa3eZBBqJkCowUR++2mP/tvzv1J+VpqO/ffP3fK1Q/6Ams/36N3RPV1+80nPuLnnu+Zl64HykS/mlWVz2JL+Jro9vtFFxyNh4sb238w0m1bdMTqqVZ6vSqcjLGHCPTCkX5+9HNrz51ZdZvcvytey+ZHvMkN0+QMj+zNdH3Y+vNR3080I01KsqnBmhUZ07hpdEJsbpTPUAgFD56/0h0aggjWfu9yvib69U6wWLS64bvRldBpprsMet4GcUDMBQg0m0j84rLu2vyFJOv7f68bsFRIIGGpx9YWmUH7ddln9N+w1Upo7KzQdUbsw8nuNJKrzl/tDQePwmW5dvmExaH5WWmjtyurFeSqeM7l2aN+wXyfO9YxMKbV26bef9o7pqLFc12X2QIx0mcEcvVeH1NI+tmW5pb1v3P/zQQUO+3XrTkZGdhbNzVJayvRHddwDozWMrhX68NKta5jrsF8XtMJTQywi1EyAUIOb+b3/8bY+udyvVzffp5LcjNEv124dOdOlLs/YL9fczDStHp2OuH9x/qR2csXUBAKGPmrvC4Wco22XdXVo7F/qZXmzQnu+1F538niwyyz4sxN2meVnhoLofQvzNHsWp5ZjYtePkox0CY1MYZ3rnnhjwVSbRYvmZo1ZnzLRKIk/YOjj69rXg23TNx0tsllV7rz2vnfOy1ZloUP5SfKPKELNBAg1uJkvff89vfVhhxzpKePmozNSbapZmBtqtV5S6EiqrolYMDgc0G8+uTJynENrl5rP94wbbVlalK2inAy99/HlcR0dcx320C7I8d5lhthw43qWYDt0n2/i9Sy5mWlaUuiQMztdZzo9t1zXMy8nfVwouiM/Ntf1RAuhZgKEGtzM7rdOafdbpyWN7DVSXTK6YHVRnu4tnZNwQ7nxrm9gSEfbrp1K3nLDMRfBLrPVi0bW45QncJcZYsdEnUcfXnLr45t0HqWnWlVZODp9dV3rNCOH4xFqJkCowc30DQzpfx7/VKW5s7TqDvYaiTcdfQM63Notl3tAK8rmTHiiNWCWq4N+ne4YGdVxuQe0aHQR74K8TBahTxKhZgKEGgAA4s9Uvr/55wwAAEgIhBoAAJAQCDUAACAhEGoAAEBCINQAAICEQKgBAAAJgVADAAASAqEGAAAkBEINAABICIQaAACQEAg1AAAgIRBqAABAQiDUAACAhJBidgHREjyM3O12m1wJAACYrOD3dvB7/FaSJtT09fVJkkpKSkyuBAAATFVfX59ycnJu+RqLMZnokwACgYAuXrwoh8Mhi8US1vd2u90qKSnR+fPnlZ2dHdb3jgWJ/vmkxP+MfL74l+ifkc8X/yL1GQ3DUF9fn4qKimS13nrVTNKM1FitVhUXF0f0GtnZ2Qn7h1VK/M8nJf5n5PPFv0T/jHy++BeJz3i7EZogFgoDAICEQKgBAAAJgVATBna7XTt27JDdbje7lIhI9M8nJf5n5PPFv0T/jHy++BcLnzFpFgoDAIDExkgNAABICIQaAACQEAg1AAAgIRBqAABAQiDUzMChQ4e0fv16FRUVyWKx6MCBA2aXFFY7d+7UypUr5XA4VFBQoMcee0wtLS1mlxU2L7zwgu65557QRlG1tbX62c9+ZnZZEfOd73xHFotFX/nKV8wuJWz++q//WhaLZcxtyZIlZpcVVhcuXNAXv/hF5eXlKSMjQ3fffbeOHTtmdllhU1ZWNu730GKx6JlnnjG7tLDw+/36q7/6K91xxx3KyMjQokWL9M1vfnNS5xjFi76+Pn3lK1/RggULlJGRodWrV+u9994zpZak2VE4Erxer6qqqrRp0yY98cQTZpcTdr/4xS/0zDPPaOXKlRoeHtbXvvY1ff7zn9cHH3ygzMxMs8ubseLiYn3nO99ReXm5DMPQ97//fT366KP6zW9+o6VLl5pdXli99957+sd//Efdc889ZpcSdkuXLtVbb70Vup+Skjh/rV25ckX333+/1q5dq5/97GeaO3euTp8+rTlz5phdWti899578vv9ofu/+93v9LnPfU5PPfWUiVWFz9/8zd/ohRde0Pe//30tXbpUx44d08aNG5WTk6Mvf/nLZpcXFl/60pf0u9/9Tv/8z/+soqIi/eAHP1BdXZ0++OADzZ8/P7rFGAgLScZrr71mdhkR1dHRYUgyfvGLX5hdSsTMmTPH+Kd/+iezywirvr4+o7y83PiP//gP47Of/azx3HPPmV1S2OzYscOoqqoyu4yI+cu//EvjgQceMLuMqHruueeMRYsWGYFAwOxSwuKRRx4xNm3aNOaxJ554wnj66adNqii8+vv7DZvNZvz7v//7mMc/85nPGF//+tejXg/TT5i03t5eSVJubq7JlYSf3+/Xj370I3m9XtXW1ppdTlg988wzeuSRR1RXV2d2KRFx+vRpFRUVaeHChXr66af1ySefmF1S2PzkJz/RihUr9NRTT6mgoED33nuvvve975ldVsQMDg7qBz/4gTZt2hT2g4fNsnr1ajU1NenUqVOSpJMnT+rdd99VfX29yZWFx/DwsPx+v9LT08c8npGRoXfffTfq9STOOC0iKhAI6Ctf+Yruv/9+LVu2zOxywub9999XbW2tBgYGlJWVpddee0133XWX2WWFzY9+9COdOHHCtPntSKupqdG+fftUWVmpS5cu6Rvf+IYefPBB/e53v5PD4TC7vBk7e/asXnjhBTU2NuprX/ua3nvvPX35y19WWlqaGhoazC4v7A4cOKCenh798R//sdmlhM3WrVvldru1ZMkS2Ww2+f1+fetb39LTTz9tdmlh4XA4VFtbq29+85u688475XQ69eqrr+rIkSNavHhx9AuK+thQglKCTz/92Z/9mbFgwQLj/PnzZpcSVj6fzzh9+rRx7NgxY+vWrUZ+fr7xX//1X2aXFRaffPKJUVBQYJw8eTL0WKJNP93oypUrRnZ2dsJMIaamphq1tbVjHnv22WeN++67z6SKIuvzn/+88Qd/8AdmlxFWr776qlFcXGy8+uqrxm9/+1vj5ZdfNnJzc419+/aZXVrYtLa2Gr/3e79nSDJsNpuxcuVK4+mnnzaWLFkS9VoINWGSyKHmmWeeMYqLi42zZ8+aXUrEPfTQQ8af/umfml1GWLz22muhv2SCN0mGxWIxbDabMTw8bHaJEbFixQpj69atZpcRFqWlpcaf/MmfjHnsH/7hH4yioiKTKoqcjz/+2LBarcaBAwfMLiWsiouLjeeff37MY9/85jeNyspKkyqKHI/HY1y8eNEwDMP4whe+YDz88MNRr4E1NbgpwzC0ZcsWvfbaa/r5z3+uO+64w+ySIi4QCMjn85ldRlg89NBDev/999Xc3By6rVixQk8//bSam5tls9nMLjHsPB6Pzpw5o3nz5pldSljcf//947ZROHXqlBYsWGBSRZHz0ksvqaCgQI888ojZpYRVf3+/rNaxX7U2m02BQMCkiiInMzNT8+bN05UrV/TGG2/o0UcfjXoNrKmZAY/Ho9bW1tD9trY2NTc3Kzc3V6WlpSZWFh7PPPOMXnnlFf34xz+Ww+FQe3u7JCknJ0cZGRkmVzdz27ZtU319vUpLS9XX16dXXnlF77zzjt544w2zSwsLh8Mxbv1TZmam8vLyEmZd1F/8xV9o/fr1WrBggS5evKgdO3bIZrPpD//wD80uLSz+/M//XKtXr9a3v/1tfeELX9DRo0f13e9+V9/97nfNLi2sAoGAXnrpJTU0NCRUS74krV+/Xt/61rdUWlqqpUuX6je/+Y127dqlTZs2mV1a2LzxxhsyDEOVlZVqbW3VV7/6VS1ZskQbN26MfjFRHxtKIG+//bYhadytoaHB7NLCYqLPJsl46aWXzC4tLDZt2mQsWLDASEtLM+bOnWs89NBDxptvvml2WRGVaGtqNmzYYMybN89IS0sz5s+fb2zYsMFobW01u6yw+rd/+zdj2bJlht1uN5YsWWJ897vfNbuksHvjjTcMSUZLS4vZpYSd2+02nnvuOaO0tNRIT083Fi5caHz96183fD6f2aWFzf79+42FCxcaaWlpRmFhofHMM88YPT09ptRiMYwE2tYQAAAkLdbUAACAhECoAQAACYFQAwAAEgKhBgAAJARCDQAASAiEGgAAkBAINQAAICEQagAAQEIg1AAAgIRAqAEAAAmBUAMAABICoQYAACSE/x/m7yWIbJkCFgAAAABJRU5ErkJggg==",
|
||
"text/plain": [
|
||
"<Figure size 640x480 with 1 Axes>"
|
||
]
|
||
},
|
||
"metadata": {},
|
||
"output_type": "display_data"
|
||
}
|
||
],
|
||
"source": [
|
||
"# Gauss-Hermite quadrature\n",
|
||
"\n",
|
||
"import numpy as np\n",
|
||
"import scipy.integrate\n",
|
||
"import matplotlib.pyplot as plt\n",
|
||
"\n",
|
||
"# Function to integrate\n",
|
||
"func = lambda x: x**4 * np.exp(-x**2)\n",
|
||
"\n",
|
||
"# Integrate this with quad to get true value\n",
|
||
"I0, err0 = scipy.integrate.quad(func, -np.inf , np.inf)\n",
|
||
"print('quad gives I = ', I0, err0)\n",
|
||
"\n",
|
||
"# Now try the Gauss-Hermite quadrature with different numbers of points and plot the error\n",
|
||
"n_vec = np.array([])\n",
|
||
"err_vec = np.array([])\n",
|
||
"\n",
|
||
"for npoints in range(1,10):\n",
|
||
" I = 0.0\n",
|
||
" x, w = np.polynomial.hermite.hermgauss(npoints)\n",
|
||
" for i in range(npoints):\n",
|
||
" # note that on the next line we include only the part of the function multiplying \n",
|
||
" # the exp(-x^2), because the exp(-x^2) is already taken into account in the weights\n",
|
||
" I += w[i] * func(x[i]) / np.exp(-x[i]**2) \n",
|
||
" err = (I-I0)/I0\n",
|
||
" print(npoints, 2*npoints-1, I, err)\n",
|
||
" n_vec = np.append(n_vec, npoints)\n",
|
||
" err_vec = np.append(err_vec, abs(err))\n",
|
||
" \n",
|
||
"plt.plot(n_vec, err_vec)\n",
|
||
"plt.yscale('log')\n",
|
||
"plt.xscale('linear')\n",
|
||
"plt.show()"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "f7749401",
|
||
"metadata": {},
|
||
"source": [
|
||
"## Maxwell-Boltzmann distribution"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "0d023c71",
|
||
"metadata": {},
|
||
"source": [
|
||
"First summarize some results for the Maxwell-Boltzmann distribution from the Wikipedia page: \n",
|
||
"\n",
|
||
"In 1D, the distribution of velocities is\n",
|
||
"\n",
|
||
"$$f(v) = \\left({m\\over 2\\pi k_BT}\\right)^{3/2} e^{-mv^2/2k_BT} ,$$\n",
|
||
"\n",
|
||
"which is normalized so that\n",
|
||
"\n",
|
||
"$$\\int_0^\\infty d^3v f(v) =1.$$\n",
|
||
"\n",
|
||
"The mean speed is\n",
|
||
"\n",
|
||
"$$\\langle v\\rangle = \\int_0^\\infty d^3v\\ v f(v).$$"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "563e6fa9",
|
||
"metadata": {},
|
||
"source": [
|
||
"To get this into a simpler form, change integration variable to $x=(m/2k_BT)^{1/2} v$ and use the spherical volume element $d^3 v= 4\\pi v^2 dv$. This gives\n",
|
||
"\n",
|
||
"$$\\langle v\\rangle = \\left({k_BT\\over \\pi m}\\right)^{1/2}\\ 4\\sqrt{2}\\ I\\hspace{1cm} I = \\int_0^\\infty dx\\ x^3\\ e^{-x^2}.$$\n",
|
||
"\n",
|
||
"We can evaluate the integral $I$ numerically (aiming for an error of 0.1%). (Spoiler: the analytic answer is 1/2).\n",
|
||
"\n",
|
||
"For Gaussian quadrature, you might think of using the Gauss-Hermite polynomials since we have an $e^{-x^2}$ in the integral, but the limits are from zero to infinity, not -infinity to infinity. One way to do it is to write the integral as\n",
|
||
"\n",
|
||
"$${1\\over 2}\\int_{-\\infty}^\\infty dx\\ \\left|x\\right|^3\\ e^{-x^2}.$$\n",
|
||
"\n",
|
||
"Note that this is not the same as having a cubic polynomial, so the usual rule about how many terms we need does not apply.\n",
|
||
"\n",
|
||
"The other way to do it is to rewrite\n",
|
||
"\n",
|
||
"$$\\int_0^\\infty dx\\ x^3\\ e^{-x^2}= \\int_0^\\infty dy\\ {y\\over 2} e^{-y}$$\n",
|
||
"\n",
|
||
"and use Laguerre polynomials."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 158,
|
||
"id": "94aa2a9a",
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"quad gives 0.5 with error 3.36729e-05; number of evaluations=75\n",
|
||
"Simpson gives 0.499296 with error 0.00140896; number of points=21\n",
|
||
"Gauss-Hermite gives 0.500992 with error 0.0019831; number of terms=15\n",
|
||
"Gauss-Laguerre gives 0.5 with error 0; number of terms=1\n"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"import numpy as np\n",
|
||
"import scipy.integrate\n",
|
||
"\n",
|
||
"func = lambda x: x**3 * np.exp(-x**2)\n",
|
||
"I0 = 0.5 # the analytic result\n",
|
||
"\n",
|
||
"# First use quad\n",
|
||
"# We'll set the relative error we are looking for to 1e-3 and also ask for \n",
|
||
"# full output which gives us information such as how many function evaluations were done\n",
|
||
"I1, err1, info = scipy.integrate.quad(func,0.0,np.inf, full_output=True, epsrel=1e-3)\n",
|
||
"print(\"quad gives %lg with error %lg; number of evaluations=%d\" % (I1, err1, info['neval']))\n",
|
||
"\n",
|
||
"# Now Simpson's rule\n",
|
||
"# Sample the function first\n",
|
||
"npoints = 21\n",
|
||
"xp = np.linspace(0.0,7.0,npoints) # integrate to x=7\n",
|
||
"fp = func(xp)\n",
|
||
"I2 = scipy.integrate.simpson(fp, xp)\n",
|
||
"print(\"Simpson gives %lg with error %lg; number of points=%d\" %(I2, abs(I2-I0)/I0, npoints))\n",
|
||
"\n",
|
||
"# Now use Gaussian quadrature\n",
|
||
"npoints = 15\n",
|
||
"x, w = np.polynomial.hermite.hermgauss(npoints)\n",
|
||
"I3 = 0.5 * np.sum(w * abs(x)**3)\n",
|
||
"print(\"Gauss-Hermite gives %lg with error %lg; number of terms=%d\" %(I3, abs(I3-I0)/I0, npoints))\n",
|
||
"\n",
|
||
"# Gauss-Laguerre\n",
|
||
"npoints = 1\n",
|
||
"x, w = np.polynomial.laguerre.laggauss(npoints)\n",
|
||
"I4 = 0.5 * np.sum(w * x)\n",
|
||
"print(\"Gauss-Laguerre gives %lg with error %lg; number of terms=%d\" %(I4, abs(I4-I0)/I0, npoints))\n",
|
||
"\n"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "4f9b0ea2",
|
||
"metadata": {},
|
||
"source": [
|
||
"In the Simpson's method, we integrate to $x=7$. You can experiment with this and see how much it changes the answer. At $x=7$, the integrand becomes comparable to the machine precision, so it seems a reasonable place to stop."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"id": "da407d88",
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": []
|
||
}
|
||
],
|
||
"metadata": {
|
||
"kernelspec": {
|
||
"display_name": "Python 3 (ipykernel)",
|
||
"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.11.5"
|
||
}
|
||
},
|
||
"nbformat": 4,
|
||
"nbformat_minor": 5
|
||
}
|