{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Siegert neuron integration\n", "\n", "#### Alexander van Meegen, 2020-12-03\n", "\n", "This notebook describes how NEST handles the numerical integration of the 'Siegert' function.\n", "\n", "For an alternative approach, which was implemented in NEST before, see Appendix A.1 in ([Hahne et al., 2017](https://www.frontiersin.org/articles/10.3389/fninf.2017.00034/full#h10)). The current approach seems to be faster and more stable, in particular in the noise-free limit.\n", "\n", "Let's start with some imports:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy.special import erf, erfcx\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "\n", "We want to determine the firing rate of an integrate-and-fire neuron with exponentially decaying post–synaptic currents driven by a mean input $\\mu$ and white noise of strength $\\sigma$. For small synaptic time constant $\\tau_{\\mathrm{s}}$ compared to the membrane time constant $\\tau_{\\mathrm{m}}$, the firing rate is given by the 'Siegert' ([Fourcaud and Brunel, 2002](https://doi.org/10.1162/089976602320264015))\n", "\n", "$$ \\phi(\\mu,\\sigma)\t=\t\\left(\\tau_{\\mathrm{ref}}+\\tau_{\\mathrm{m}}\\sqrt{\\pi}I(\\tilde{V}_{\\mathrm{th}},\\tilde{V}_{\\mathrm{r}})\\right)^{-1} $$\n", "\n", "with the refractory period $\\tau_{\\mathrm{ref}}$ and the integral\n", "\n", "$$ I(\\tilde{V}_{\\mathrm{th}},\\tilde{V}_{\\mathrm{r}})\t=\t\\int_{\\tilde{V}_{\\mathrm{r}}}^{\\tilde{V}_{\\mathrm{th}}}e^{s^{2}}(1+\\mathrm{erf}(s))ds $$\n", "\n", "involving the shifted and scaled threshold voltage $\\tilde{V}_{\\mathrm{th}}=\\frac{V_{\\mathrm{th}}-\\mu}{\\sigma}+\\frac{\\alpha}{2}\\sqrt{\\frac{\\tau_{\\mathrm{s}}}{\\tau_{\\mathrm{m}}}}$, the shifted and scaled reset voltage $\\tilde{V}_{\\mathrm{r}}=\\frac{V_{\\mathrm{r}}-\\mu}{\\sigma}+\\frac{\\alpha}{2}\\sqrt{\\frac{\\tau_{\\mathrm{s}}}{\\tau_{\\mathrm{m}}}}$, and the constant $\\alpha=\\sqrt{2}\\left|\\zeta(1/2)\\right|$ where $\\zeta(x)$ denotes the Riemann zeta function. \n", "\n", "Numerically, the integral in $I(\\tilde{V}_{\\mathrm{th}},\\tilde{V}_{\\mathrm{r}})$ is problematic due to the interplay of $e^{s^{2}}$ and $\\mathrm{erf}(s)$ in the integrand. Already for moderate values of $s$, it causes numerical problems (note the order of magnitude):" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAENCAYAAAAL98L+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnVElEQVR4nO3de3wU9dn38c8VlIOCHMQDcm5VFGurmIKHuz5WeRShBYu2kIpFBSIUKGpVoCIeUCoKyikgKKhYBb09UBQQEG9EW9QEUA4iFBFuUxAQtRDCmev5YxefTdjAkt3s7Cbf9+u1L3Z+MzvzzW/JXvnNzM6YuyMiInJIRtABREQktagwiIhIESoMIiJShAqDiIgUocIgIiJFqDCIiEgRKgwiIlKECoOIiBRxXNABROJlZtcAPwXOBv7o7vsCjiSS1jRikLTn7nPc/XGgAKgcdB6RdKfCIGnFzI43sxfM7J9m9pGZ1Qu39wLmuPvOYsv/1cxuDyLrkZhZMzNbamY7zOxP4baYsprZx2Z23tHaREpLhUHSzc+B7e5+KXCxu28ys57AVcA5Zlb70IJmdgrwB2BCeLqPmeWZ2R4zey6A7JHuARa4ew13H10861EMBx6KoU2kVFQYJN0sBnaY2TSgLYC7P+XuN7j7SHf/LmLZm4FZ7r4rPL0ReBiYHOvGzOwBM3sgIclD6zt0XK8xsDJi1s0UzXokM4BfHhotHaFNpFRUGCRlmVkPM/vMzP5jZrPN7FSgkrsPAHoS+qv7SK4F3js04e6vu/t0YFsZZD3DzF4zs61m9uWh3UPheevNrL+ZLQN2mtm7wC+BsWZWYGZnF88aft2PzOwtM/sm3Afzwj/HbkIF8uqIn+2wNpHS0llJkpLM7C/A9UB74H+BcYT+2j/JzBoCxwODj7Ka84HVZZkTwMwygDeBvwNZQAPgHTNb7e5zwotlAe2Ab9x9l5ktAP7m7s+E1xEt6xRgKtCB0M97UcS8VcDPii0frU3kmKkwSMoJjwwGARe4+9pw2yRgnLtfeAyrqgXsSHzCw/wcOMXdD+3jX2dmTwOdgUOFYbS7f3WEddTi8Kw/BioRGiXtBv4RMW8HUHy3UbQ2kWOmwiCp6CpCp51+bGaH2gxYeozr+Q6ocawbN7O3gP8KT1YNt90env7A3X9V7CWNgTPM7PuItkrA+xHTRyoKJWW9EbgXGGxmfwfudvdvw/NqAN8XWz5am8gx0zEGSUV1gDfcvVbEo6a7X3GM61lG6Etvx8Tdf3Vou8CjwKMROYoXBQh96H9ZLG8Nd28budpjzeru77r7VUBzQruIbo6YfS7wabF1RGsTOWYqDJKKlhA6w6YFgJmdZGYdLGL4EKNZwP85NGFmx5lZVcK7Z8ysasRZQvH4GNgePsBczcwqmdlPzOzncWTtaGZnhX/mGkBt4JPwvCqEjjfMi1j+sDaR0lJhkJTj7osInZP/mpkVAJ8BbfzYb1A+BWhrZtXC04OAXcAAoEv4+aAE5D0A/Bq4APgS+AZ4BqgZR9b/InSW0g5CReNRd383PK89oe9AbIx4fbQ2kVKxY/9dE0kfZjYU2OLuI4POcjSxZjWzj4Bu7r7iSG0ipaXCICIiRWhXkoiIFKHCICIiRagwiIhIESoMIiJSREp+87lu3brepEmToGOIiKSVxYsXf+Pup8S7npQsDE2aNCEvLy/oGCIiacXMNiRiPdqVJCIiRagwiIhIESoMIiJShAqDiIgUUSaFwcxONLPFZvarI7WJiEjqiakwmNlkM9tiZiuKtbcxs9VmttbMBkTM6g+8Umw10dpERCTFxDpieA5oE9lgZpWAHEI3MW8OZJlZczNrTegyyZsjlj2sTUREUlNM32Nw94Vm1qRYc0tgrbuvAzCzaYRuWl4dOJFQsdhlZrOAXxZvc/eDkSszs2wgG6BRo0al/oFERCqiTz75JGHriucLbvUpeh/bfKCVu/cBMLObgW/CBeDeKG1FuPtEYCJAZmamrgUuInIMhg8fnrB1xVMYot1m8YcPdHd/7rCZUdpERCQ+//73v3n55ZcTtr54zkrKBxpGTDcAdFtBEZEky8nJ4cCBAwlbXzyFIRc4y8yamllloDMwIzGxREQkFoWFhUyYMIHrrrsuYeuM9XTVqcAioJmZ5ZtZN3ffD/QB5gCrgFfcfWXCkomIyFFNmTKFb7/9ljvuuCNh60zJez5nZma6rq4qInJkBw8e5LzzzuPEE08kNzeXjIyMxe6eGe96U/Ky2yIicnRz5szh888/529/+xtm0c4HKh1dK0lEJE0NHz6c+vXr89vf/jah61VhEBFJQ0uWLOHdd9+lX79+VK5cOaHrVmEQEUlDI0aMoEaNGmRnZyd83SoMIiJpZsOGDbz88stkZ2dTs2bNhK9fhUFEJM2MGjUKM6Nfv35lsn4VBhGRNPL999/z9NNP06lTJxo2bHj0F5SCCoOISBqZMGECBQUF3HXXXWW2DRUGEZE0sXfvXkaNGkXr1q254IILymw7+oKbiEiamDp1Kps2beLZZ58t0+1oxCAikgbcneHDh3P++edz9dVXl+m2NGIQEUkDc+bMYcWKFTz//PMJvfxFNBoxiIikgeHDh3PGGWfQuXPnMt+WCoOISIpbsmQJ8+fPL5PLX0SjwiAikuJGjBhB9erVy+TyF9GoMIiIpLAvv/zyh8tf1KpVKynbTPjBZzM7F+gH1AXmu/t4M2sEjAW+Ada4+6OJ3q6ISHn0+OOPk5GRwZ133pm0bcZ6a8/JZrbFzFYUa29jZqvNbK2ZDQBw91Xu3hP4HXDoTkJnAzPd/VageQLzi4iUW19//TWTJ0+ma9eu1K9fP2nbjXVX0nNAm8gGM6sE5ADXEvqwzzKz5uF57YEPgPnhxZcCnc3sXeB/4o8tIlL+jRw5kn379nHPPfckdbsxFQZ3Xwh8W6y5JbDW3de5+15gGtAhvPwMd78UuDG87C3A/e5+JdAu2jbMLNvM8swsb+vWraX4UUREyo/vv/+ecePGccMNN3DWWWclddvxHGOoD3wVMZ0PtDKzK4COQBVgVnje28ADZvZ7YH20lbn7RGAiQGZmpseRS0Qk7Y0bN44dO3YwYMCApG87nsIQ7at37u4LgAXFGlcAN8SxLRGRCqOwsJCRI0dyzTXXcOGFFyZ9+/GcrpoPRF4MvAGwMb44IiLy7LPPsnXrVgYOHBjI9uMpDLnAWWbW1MwqA52BGYmJJSJSMe3bt4/HH3+cSy65hMsvvzyQDLGerjoVWAQ0M7N8M+vm7vuBPsAcYBXwiruvLLuoIiLl37Rp09iwYQMDBw4s84vllcTcU+84b2Zmpufl5QUdQ0QkqQ4ePMj5559PRkYGn376KRkZx7ZTx8wWu3vm0Zc8Ml12W0QkRbz55pt89tln/O1vfzvmopBIulaSiEgKcHf++te/0rRpUzp16hRoFo0YRERSwDvvvMNHH33E+PHjOe64YD+aNWIQEQmYu/Pggw/SoEEDbrnllqDjaMQgIhK0BQsW8I9//IOxY8dSpUqVoONoxCAiErSHHnqIevXq0a1bt6CjABoxiIgE6v3332fBggWMHDmSqlWrBh0H0IhBRCRQQ4YM4dRTT6VHjx5BR/mBCoOISEAWLVrEvHnzuPvuuznhhBOCjvMDFQYRkYAMGTKEunXr0rNnz6CjFKHCICISgNzcXGbPns2f//xnqlevHnScIlQYREQCMGTIEOrUqUPv3r2DjnIYFQYRkSRbunQpb775JnfccQc1atQIOs5hVBhERJLs4YcfpmbNmvTt2zfoKFGpMIiIJNHSpUt5/fXX6devHzVr1gw6TlQqDCIiSTR48GBq1arFHXfcEXSUEpVJYTCzc83sKTN71cx6hduuM7OnzezvZnZ1WWxXRCSVffTRR7z11lvcfffd1KpVK+g4JYq5MJjZZDPbYmYrirW3MbPVZrbWzAYAuPsqd+8J/A7IDLdNd/cewM1AsBcbFxEJwODBg6lbty5/+tOfgo5yRMcyYngOaBPZYGaVgBzgWqA5kGVmzcPz2gMfAPOLrWdQ+DUiIhXGBx98wNy5c+nfv3/KfW+huJgLg7svBL4t1twSWOvu69x9LzAN6BBefoa7XwrcCGAhw4DZ7r4kIelFRNKAuzNo0CBOP/10/vjHPwYd56jivbpqfeCriOl8oJWZXQF0BKoAs8Lz+gKtgZpmdqa7PxW5IjPLBrIBGjVqFGcsEZHU8e677/Lee+8xatSolLomUkniLQwWpc3dfQGwoFjjaGB0SSty94nARIDMzEyPM5eISEpwd+677z4aNGhAdnZ20HFiEm9hyAcaRkw3ADbGuU4RkXLj7bffZtGiRYwfPz5l7rdwNPGerpoLnGVmTc2sMtAZmBF/LBGR9HdotNCkSRNuvfXWoOPE7FhOV50KLAKamVm+mXVz9/1AH2AOsAp4xd1Xlk1UEZH08ve//53FixczePBgKleuHHScmJl76u3Oz8zM9Ly8vKBjiIiU2sGDB7ngggvYvXs3n332GccdV/Z3Ujazxe6eGe96dM9nEZEy8NJLL7F8+XJeeumlpBSFRNKIQUQkwfbu3UuzZs2oVasWixcvJiMjOZel04hBRCRFTZgwgfXr1zN79uykFYVESr/EIiIpbMeOHQwZMoQrrriCa665Jug4paLCICKSQE8++SRbt27l0UcfxSzad4BTnwqDiEiCbN26lccff5yOHTvSqlWroOOUmgqDiEiCPPLIIxQWFvLwww8HHSUuKgwiIgmwfv16xo8fzy233MK5554bdJy4qDCIiCTA/fffT0ZGBg888EDQUeKmwiAiEqfly5fzwgsv0LdvXxo0aBB0nLipMIiIxGnAgAGcdNJJDBgwIOgoCaEvuImIxOGdd95h1qxZPPbYY9SpUyfoOAmhS2KIiJTSgQMHuOiii/jPf/7DqlWrAr/fgi6JISISsClTpvDpp58yderUwItCIukYg4hIKezcuZN7772XVq1a0alTp6DjJJRGDCIipTBixAg2bdrEq6++mraXviiJRgwiIsdo06ZNPPbYY9xwww1ceumlQcdJuISPGMzsOqAdcCqQ4+5zzSwDGAKcBOS5+/OJ3q6ISLLcd9997N27l0cffTToKGUiphGDmU02sy1mtqJYexszW21ma81sAIC7T3f3HsDNwKEdbx2A+sA+ID9h6UVEkmzZsmVMnjyZPn368OMf/zjoOGUi1l1JzwFtIhvMrBKQA1wLNAeyzKx5xCKDwvMBmgGL3P1OoFc8gUVEgnT33XdTq1YtBg0aFHSUMhNTYXD3hcC3xZpbAmvdfZ277wWmAR0sZBgw292XhJfNB74LPz8QbRtmlm1meWaWt3Xr1mP+QUREytrs2bOZO3cugwcPLjdfZosmnoPP9YGvIqbzw219gdbADWbWMzzvdeAaMxsDLIy2Mnef6O6Z7p55yimnxBFLRCTx9u7dy+23387ZZ5/NH//4x6DjlKl4Dj5HOz/L3X00MLpYYyHQLY5tiYgEasyYMaxZs4aZM2dSuXLloOOUqXhGDPlAw4jpBsDG+OKIiKSer7/+mgcffJC2bdvStm3boOOUuXgKQy5wlpk1NbPKQGdgRmJiiYikjr/85S/s3r2bJ598MugoSRHr6apTgUVAMzPLN7Nu7r4f6APMAVYBr7j7yrKLKiKSfLm5uTz77LP069ePs88+O+g4SaGrq4qIlODgwYNcdtllfPnll6xZs4aTTjop6EhHpKurioiUsRdffJEPP/yQyZMnp3xRSCSNGEREotixYwfNmjWjQYMGfPjhh2RkpP6l5TRiEBEpQ0OHDmXTpk288cYbaVEUEqli/bQiIjFYs2YNTzzxBH/4wx9o1apV0HGSToVBRCSCu9O7d2+qVavGsGHDgo4TCO1KEhGJ8Morr/DOO+8wduxYTj/99KDjBEIHn0VEwrZv384555xDvXr1+Pjjj6lUqVLQkY6JDj6LiCTY4MGD+frrr5k+fXraFYVE0jEGERHgk08+YcyYMdx22220bNky6DiBUmEQkQrv4MGD9OrVi5NPPpmhQ4cGHSdw2pUkIhXepEmT+PDDD3n++eepXbt20HECpxGDiFRoW7dupX///vziF7/gpptuCjpOSlBhEJEKrX///uzYsYPx48djFu3+YxWPCoOIVFjz58/n2Wef5c9//jPnnXde0HFShr7HICIVUmFhIT/96U8xM5YtW0a1atWCjhQ3fY9BRCQODz74IF988QXvvvtuuSgKiZTwwmBm1wHtgFOBHHefa2YnAuOAvcACd38x0dsVEYnV0qVLGTFiBN26deOXv/xl0HFSTqy39pxsZlvMbEWx9jZmttrM1prZAAB3n+7uPYCbgU7hRTsCr4bb2ycuvojIsdm/fz/du3enbt26PP7440HHSUmxHnx+DmgT2WBmlYAc4FqgOZBlZs0jFhkUng/QAPgq/PxAacOKiMRr5MiRLFmyhLFjx+o7CyWIqTC4+0Lg22LNLYG17r7O3fcC04AOFjIMmO3uS8LL5hMqDjFvU0Qk0datW8fgwYPp0KED119/fdBxUlY8H9L1+f+jAAh9+NcH+gKtgRvMrGd43uvA9WY2Hngz2srMLNvM8swsb+vWrXHEEhE5nLtz2223cfzxx5OTk6PvLBxBPAefo/Wqu/toYHSxxp3ALUdambtPBCZC6HTVOHKJiBzmmWee4Z133mH8+PHUr18/6DgpLZ4RQz7QMGK6AbAxvjgiIom3YcMG7rzzTq688kqys7ODjpPy4ikMucBZZtbUzCoDnYEZiYklIpIY7k63bt2A0MXyMjJ0mPNoYj1ddSqwCGhmZvlm1s3d9wN9gDnAKuAVd19ZdlFFRI7dhAkTmD9/PsOHD6dJkyZBx0kLuiSGiJRbX375Jeeffz6XXHIJc+fOLfcHnBN1SQyNqUSkXDp48CC33norGRkZTJo0qdwXhUTStZJEpFwaN24cCxYs4Omnn6ZRo0ZBx0krGjGISLnzxRdf0L9/f6655pofDjxL7FQYRKRc2b9/PzfddBPHH388Tz/9tHYhlYJ2JYlIuTJ06FAWLVrESy+9RMOGDY/+AjmMRgwiUm58+OGHPPTQQ9x4441kZWUFHSdtqTCISLlQUFBAly5dqF+/Pjk5OUd/gZRIu5JEpFy4/fbbWbduHQsWLKBmzZpBx0lrGjGISNp74403mDRpEgMGDODyyy8POk7aU2EQkbS2ceNGevToQYsWLXjggQeCjlMuqDCISNo6ePAgN998M4WFhbz44otUrlw56Ejlgo4xiEjaGjZsGPPmzWPChAmcc845QccpNzRiEJG09MEHH3DffffRuXNnevToEXScckWFQUTSzrZt28jKyqJp06ZMmDBB325OMO1KEpG04u507dqVLVu2sGjRIk466aSgI5U7KgwiklaefPJJZs6cyZgxY2jRokXQccol7UoSkbTx8ccf079/fzp27Ejv3r2DjlNuJbwwmNmPzGySmb1arP1EM1tsZr9K9DZFpPzbtm0bv/vd76hfv75uvFPGYr3n82Qz22JmK4q1tzGz1Wa21swGALj7OnePdgH0/sAr8UcWkYrmwIED/P73v2fTpk28+uqr1KpVK+hI5VqsI4bngDaRDWZWCcgBrgWaA1lm1jzai82sNfAZsLnUSUWkwrr//vuZO3cuOTk5ZGbGfUtjOYqYDj67+0Iza1KsuSWw1t3XAZjZNKADoQJQ3C+BEwkVkF1mNsvdD0YuYGbZQDag2/CJyA9mzJjBI488Qvfu3enevXvQcSqEeI4x1Ae+ipjOB+qb2clm9hRwoZkNBHD3e939duAl4OniRSG8zER3z3T3zFNOOSWOWCJSXvzrX//ipptu4qKLLmLMmDFBx6kw4jldNdqRH3f3bUDPaC9w9+fi2J6IVCA7d+6kY8eOHHfccbz22mtUrVo16EgVRjyFIR+IvG9eA2BjfHFEREJfYuvevTsrV67k7bffpnHjxkFHqlDi2ZWUC5xlZk3NrDLQGZiRmFgiUpH99a9/Zdq0aTzyyCNcffXVQcepcGI9XXUqsAhoZmb5ZtbN3fcDfYA5wCrgFXdfWXZRRaQimD59Ovfeey9ZWVkMGDAg6DgVkrl70BkOk5mZ6Xl5eUHHEJEkW7ZsGZdeeinNmzfnvffeo1q1akFHSitmttjd4z6fV5fEEJGUsGXLFtq3b0/NmjWZPn26ikKAdBE9EQnc3r17uf7669m8eTPvv/8+Z5xxRtCRKjQVBhEJlLvTq1cvPvjgA6ZNm6ZvNqcA7UoSkUA98sgjTJ48mcGDB9OpU6eg4wgqDCISoClTpnDfffdx00038cADDwQdR8JUGEQkEPPnz6dbt25cddVVPPPMM7qMdgpRYRCRpFu+fDkdO3bknHPO4bXXXqNy5cpBR5IIKgwiklT5+flce+211KhRg1mzZlGzZs2gI0kxOitJRJLmu+++o23btmzfvp3333+fhg0bHv1FknQqDCKSFDt37qRdu3asXr2amTNn8rOf/SzoSFICFQYRKXN79uyhY8eOfPTRR/z3f/83rVu3DjqSHIEKg4iUqQMHDtClSxfmzp3LpEmT6NixY9CR5Ch08FlEyoy707NnT1599VVGjBjBrbfeGnQkiYEKg4iUCXfnrrvu4plnnmHQoEHceeedQUeSGKkwiEjCuTv33HMPTzzxBH379uWhhx4KOpIcAxUGEUkod2fgwIEMHz6c3r17M2rUKH2rOc0k5eCzmTUCxgLfAGvc/dFkbFdEksvdGTRoEMOGDaNXr16MGTNGRSENlXrEYGaTzWyLma0o1t7GzFab2VozO3RfvrOBme5+K9A8jrwikqLcnfvvv5+hQ4eSnZ3N2LFjVRTSVDy7kp4D2kQ2mFklIAe4llAByDKz5sBSoLOZvQv8TxzbFJEU5O4MGDCAIUOG0L17d8aPH09GhvZUp6tSv3PuvhD4tlhzS2Ctu69z973ANKADcAtwv7tfCbQr7TZFJPUcPHiQ3r1789hjj9GrVy8mTJigopDmEv3u1Qe+ipjOD7e9DfzJzJ4C1kd7oZllm1memeVt3bo1wbFEpCzs37+fm2++mfHjx9O/f39ycnJUFMqBRB98jrZD0d19BXDDkV7o7hOBiQCZmZme4FwikmB79uwhKyuLN954g6FDhzJw4MCgI0mCJLow5AORl0tsAGxM8DZEJGA7duzg+uuvZ968eYwePZq+ffsGHUkSKNGFIRc4y8yaAv8GOgO/T/A2RCRAmzZtom3btixfvpznnnuOrl27Bh1JEiye01WnAouAZmaWb2bd3H0/0AeYA6wCXnH3lYmJKiJBW7VqFRdffDFr165l5syZKgrlVKlHDO6eVUL7LGBWqROJSEp6//33ad++PVWrVuW9996jRYsWQUeSMqLTB0TkqKZNm0br1q05/fTTWbRokYpCOafCICIlOnjwIH/5y1/IysqiVatW/OMf/6BJkyZBx5Iyphv1iEhU27dv58Ybb+Stt94iOzubMWPGULly5aBjSRKoMIjIYdauXUv79u1Zs2YNOTk59OrVS9c9qkBUGESkiFmzZtGlSxfMjLlz53LllVcGHUmSTMcYRAQIXd5i4MCBtGvXjoYNG5Kbm6uiUEFpxCAibNy4kaysLBYuXEiPHj0YNWoU1apVCzqWBESFQaSCmzNnDjfddBM7d+7khRdeoEuXLkFHkoBpV5JIBVVYWEjfvn1p06YNp5xyCrm5uSoKAqgwiFRIixcv5qKLLmLs2LH069ePvLw8mjfXzRUlRIVBpALZt28fDz/8MBdffDE7duxg3rx5jBw5UscTpAgdYxCpIHJzc+nevTvLli2jU6dOjBs3jjp16gQdS1KQRgwi5VxBQQF33HEHF198Md988w3Tp09n2rRpKgpSIhUGkXLK3Zk+fTo/+clPGDlyJLfddhufffYZHTp0CDqapDjtShIph1asWMHtt9/O/PnzOe+88/jggw+47LLLgo4laUIjBpFyZNu2bfTp04ef/exnLFmyhLFjx/LJJ5+oKMgx0YhBpBzYvn07Tz75JE888QQFBQX06tWLBx98kJNPPjnoaJKGklIYzOw6oB1wKpDj7nOTsV2R8q6wsJCcnByGDRvGtm3b+M1vfsOQIUM477zzgo4maSyeez5PNrMtZraiWHsbM1ttZmvNbACAu0939x7AzUCnuBKLCN9//z2PPvooP/7xj7nnnnvIzMwkNzeX119/XUVB4hbPiOE5YCww5VCDmVUCcoD/C+QDuWY2w90/Cy8yKDxfREohPz+fkSNHMmHCBAoKCrj66qt5+eWXufzyy4OOJuVIqQuDuy80sybFmlsCa919HYCZTQM6mNkq4FFgtrsvibY+M8sGsgEaNWpU2lgi5Y67s2DBAp566ilef/113J1OnTpx1113ceGFFwYdT8qhRB9jqA98FTGdD7QC+gKtgZpmdqa7P1X8he4+EZgIkJmZ6QnOJZJ2tm7dypQpU5g4cSJr1qyhdu3a9OnTh379+um+y1KmEl0Yot37z919NDA6wdsSKXe2b9/O9OnTmTp1KvPmzePAgQNcdtllDBo0iBtuuEHXNJKkSHRhyAcaRkw3ADYmeBsi5cqWLVuYPXs2M2bMYNasWezevZvGjRtz11130aVLF37yk58EHVEqmEQXhlzgLDNrCvwb6Az8PsHbEElr+/fvZ+nSpbz99tvMnDmTjz/+GHenXr16dO/enaysLC655BLMog3ARcpeqQuDmU0FrgDqmlk+cL+7TzKzPsAcoBIw2d1XJiSpSJravXs3eXl5vPfeeyxcuJB//vOfFBQUYGb8/Oc/58EHH6Rdu3ZceOGFKgaSEuI5KymrhPZZwKxSJxJJY4WFhXz66acsWbKEJUuWsHjxYlauXMn+/fsBOP/88+natSuXX345V1xxBaeeemrAiUUOp0tiiBwjd2fz5s2sXr2azz///Id/P//8c9avX4976KS6unXrctFFF9G2bVtatmzJL37xC12iQtKCCoNIBHdn+/btbN68mY0bN7JhwwY2bNjA//7v/xZ5vmfPnh9ec8IJJ9CsWTNatWpF165dueCCC2jRogUNGjTQriFJSyoMUq7t3buX77777rDHtm3b2Lx5c9RH5If+IaeffjqNGzfmggsuoEOHDjRu3Jizzz6bc845hwYNGpCRoQsVS/mhwiAp58CBA+zcuZOCggJ27NhBQUFBkUdJbd9///1hBaCwsLDE7VSqVIlTTz2V0047jdNOO41zzz33h+ennXYa9erVo1GjRjRs2JCqVasmsQdEgqXCkAZGjx5Ny5Ytufjii5O2TXdn37597N2794fHnj17ijzftWsXu3fvZteuXT88IqePNC/asoWFhRQUFLBr166Yc1apUoXq1atTvXp1atasSe3atTnzzDOpXbv2ER916tTh5JNP1l/6IlGoMKSBfv36AXDaaadx1VVXkZmZya5du6hevTq7du1iz549RT60S/owj3Xenj172LdvX0KyV61alapVq1KtWjWqVatW5PmJJ55I3bp1i7TXqFHjhw/66tWrF5kuPq969eocf/zxCckpIv+fHTqDIpVkZmZ6Xl5e0DFSRiwHMDMyMqhSpQpVqlShcuXKRR6xth3rsoc+4It/4B+arlKliv4iF0kiM1vs7pnxrkcjhjRw4okncuaZZ/LrX/+aVq1a0bx5c+rVq8fOnTs54YQTqFKlCpUqVQo6poiUEyoMaaJ169YMGTKkSJsuqCYiZUHjfBERKUKFQUREilBhEBGRIlQYRESkCBUGEREpQoVBRESKUGEQEZEiVBhERKSIlLwkhpntAFYHnSMGdYFvgg4RA+VMLOVMrHTImQ4ZAZq5e414V5Kq33xenYjrfZQ1M8tTzsRRzsRSzsRJh4wQypmI9WhXkoiIFKHCICIiRaRqYZgYdIAYKWdiKWdiKWfipENGSFDOlDz4LCIiwUnVEYOIiAQksMJgZr81s5VmdtDMMovNG2hma81stZldU8Lr65jZPDP7V/jf2knI/LKZfRJ+rDezT0pYbr2ZLQ8vl/Rb0ZnZA2b274isbUtYrk24j9ea2YAAcj5uZp+b2TIze8PMapWwXCD9ebT+sZDR4fnLzKxFsrKFt9/QzP7HzFaFf5f6RVnmCjP7T8T/hcHJzBiR44jvYdB9Gc7QLKKfPjGz7WZ2e7FlAulPM5tsZlvMbEVEW0yfgaX6PXf3QB7AuUAzYAGQGdHeHPgUqAI0Bb4AKkV5/WPAgPDzAcCwJOcfAQwuYd56oG6AffsAcNdRlqkU7tsfAZXDfd48yTmvBo4LPx9W0nsYRH/G0j9AW2A2YMDFwEdJzlgPaBF+XgNYEyXjFcBbycxVmvcw6L4s4f3/GmicCv0JXA60AFZEtB31M7C0v+eBjRjcfZW7R/sSWwdgmrvvcfcvgbVAyxKWez78/HngujIJGoWFbsL8O2BqsrZZBloCa919nbvvBaYR6tOkcfe57r4/PPkh0CCZ2z+KWPqnAzDFQz4EaplZvWQFdPdN7r4k/HwHsAqon6ztJ1igfRnFVcAX7r4hwAw/cPeFwLfFmmP5DCzV73kqHmOoD3wVMZ1P9P/sp7n7Jgj9ggCnJiHbIb8ANrv7v0qY78BcM1tsZtlJzBWpT3hIPrmEIWas/ZwstxL6izGaIPozlv5JmT40sybAhcBHUWZfYmafmtlsMzsvucl+cLT3MGX6MqwzJf/hlwr9CbF9BpaqX8v0m89m9g5wepRZ97r730t6WZS2pJ06FWPmLI48WrjM3Tea2anAPDP7PFzxk5ITGA8MIdRvQwjt9rq1+CqivDbh/RxLf5rZvcB+4MUSVlPm/RlFLP0T6P/VH0KYVQdeA2539+3FZi8htDukIHysaTpwVpIjwtHfw5ToSwAzqwy0BwZGmZ0q/RmrUvVrmRYGd29dipflAw0jphsAG6Mst9nM6rn7pvCQc0tpMhZ3tMxmdhzQEbjoCOvYGP53i5m9QWg4l9APslj71syeBt6KMivWfo5LDP3ZFfgVcJWHd4pGWUeZ92cUsfRPUvrwSMzseEJF4UV3f734/MhC4e6zzGycmdV196Re9yeG9zDwvoxwLbDE3TcXn5Eq/RkWy2dgqfo1FXclzQA6m1kVM2tKqBp/XMJyXcPPuwIljUASrTXwubvnR5tpZieaWY1DzwkdYF0RbdmyUmzf7G9K2H4ucJaZNQ3/hdSZUJ8mjZm1AfoD7d29sIRlgurPWPpnBvCH8Bk1FwP/OTS0T4bwsa5JwCp3f6KEZU4PL4eZtST0O78tWRnD243lPQy0L4spcY9AKvRnhFg+A0v3e57so+sRR8t/Q6ia7QE2A3Mi5t1L6Ej6auDaiPZnCJ/BBJwMzAf+Ff63TpJyPwf0LNZ2BjAr/PxHhI78fwqsJLTLJNl9+wKwHFgW/k9Qr3jO8HRbQmeyfBFQzrWE9n9+En48lUr9Ga1/gJ6H3n9Cw/Sc8PzlRJxdl6R8/0Vot8CyiD5sWyxjn3C/fUroAP+lAbzPUd/DVOrLiKwnEPqgrxnRFnh/EipUm4B94c/NbiV9Bibi91zffBYRkSJScVeSiIgESIVBRESKUGEQEZEiVBhERKQIFQYRESlChUFERIpQYRARkSJUGEREpIj/BzxRsP5GJB1tAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "s = np.linspace(-10, 10, 1001)\n", "plt.plot(s, np.exp(s**2) * (1 + erf(s)), c=\"black\")\n", "plt.xlim(s[0], s[-1])\n", "plt.yscale(\"log\")\n", "plt.title(r\"$e^{s^{2}}(1+\\mathrm{erf}(s))$\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The main trick here is to use the [scaled complementary error function](https://en.wikipedia.org/wiki/Error_function#Complementary_error_function)\n", "\n", "$$\\mathrm{erf}(s)=1-e^{-s^{2}}\\mathrm{erfcx}(s)$$\n", "\n", "to extract the leading exponential contribution. For positive $s$, we have $0\\le\\mathrm{erfcx}(s)\\le1$, i.e. the exponential contribution is in the prefactor $e^{-s^{2}}$ which nicely cancels with the $e^{s^{2}}$ in the integrand. In the following, we separate three cases according to the sign of $\\tilde{V}_{\\mathrm{th}}$ and $\\tilde{V}_{\\mathrm{r}}$ because for a negative arguments, the integrand simplifies to $e^{s^{2}}(1+\\mathrm{erf}(-s))=\\mathrm{erfcx}(s)$. Eventually, only integrals of $\\mathrm{erfcx}(s)$ for positive $s\\ge0$ need to be solved numerically which are certainly better behaved:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEKCAYAAAAW8vJGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhQklEQVR4nO3deXRV9d3v8fc3IYEIYVaBIISQxEAgIDPRAhWsQkHpI0NoRa22drLVetfFsU+nddUuqtV7aXU5PX0QyyBSBEEeAyhFijIaJYQQiISxhqECoiQMv/tHTmJOmAJJzt7nnM9rrb2Ss7PP3t/zW3I+/n77t/c25xwiIiKVYrwuQERE/EXBICIiQRQMIiISRMEgIiJBFAwiIhJEwSAiIkEUDCIiEkTBICIiQRQMIjWY2dVmttHMjprZL+q4ryfM7P5abLfGzDLrciyR+mK68lkkmJm9DBxxzv2yjvu5HPgISHXOfXWBbScAE51zt9blmCL1QT0GkQAzaxT4tTOQXw+7vBNYfKFQCFgAfNPM2tfDcUXqRMEgEc/MOpjZG2a238w+rT48ZGY7zOxBM/sYOGZmy4FvAtPM7AszSzezq8xsXuD9B81sWuC9Xc3skJn1qXacA2Y2LLD7kcCKasdKMbO3AtscNrPcyr85544D64FvNXBziFyQgkEimpnFAAuBPCAJGA7cb2Y3VttsEvBtoKVz7npgJXCvc64ZsB14CygBkgP7mAXgnNsOPAi8ZmaXAf8F/NU5915gvz2BwmrHmQ68DVwZWH5To9wCoFddP7NIXSkYJNL1By53zv3OOVfunCsGXgRyqm3zf51zu84x5DMA6AD8b+fcMefccefc+5V/dM69CBQBHwLtgUervbclcLTa665ALBAb2M+qGsc6GniPiKcUDBLpOgMdzOzzygV4hIr/Y6+06zzvvwoocc6dPM82LwI9gP/nnCurtv7fQGK1198DbgH2mtnLZta6xn4Sgc/P92FEQkHBIJFuF/Cpc65ltSXROTeq2jbnm5q3C+hU7cR0EDNrBjwDvAz8psaX/cdAetVBnFvunBsOdKdiyOjOGrvrRsWQl4inFAwS6dYARwInmBPMLNbMephZ/4t4/z7gSTNramZNzOzaan9/FljvnPsBsAh4vtrfFgNDAczsP8wszcyMip5BKyqmshL4e2OgL5CLiMcUDBLRnHOngDFAb+BT4ADwEtDiIt+fCuwEdgMTAczsFuAm4MeBzR8A+pjZ9wKvpwOjzCwBuI6KGUpHqQiMJ51zy6sd6mbgPefc3kv6oCL1SBe4iTQgM3scKHXOPXOB7T4E7nbObQpJYSLnoWAQEZEgGkoSEZEgCgYREQmiYBARkSBnnZvtF23btnXJyclelyEiEjbWr19/wDl3eV324etgSE5OZt26dV6XISISNsyspK770FCSiIgEUTCIiEgQBYOIiARRMIiISBAFg4iIBAnZrCQzawr8BSin4mZhr4Xq2CIiUnt16jGY2StmVmpmm2qsv8nMCs1sm5k9FFj9H8Bc59wPqbiTpIiI+FBdh5L+SsVth6uYWSzwZyoehN4dmGRm3YGOfP2krFO12fmxY8fqWJ6IiFysOgWDc+4fwKEaqwcA25xzxc65cioenH4LFfex73ih45rZPWa2zszW/etf/6pLeSIicgka4uRzEsHP0N0dWDcPuNXMngMWnuvNzrkXnHP9nHP9Gjdu3ADliYjI+TTEyWc7yzrnnDsGfP9idnTixIn6qUhERGqtIXoMu4Grqr3uCFzS4wpPnjxZLwWJiEjtNUSPYS2QZmZdgD1ADvDdi9mBmY0BxsTHxzdAeSIicj51na46E1gNXG1mu83sbufcSeBe4H+AAmCOcy7/YvbrnFvonLtHjx0VEQk9Xz/zOSYmxp06dQqzs522EBGRmsxsvXOuX1324etbYjjnOHz4sNdliIhEFV8HA8D+/fu9LkFEJKr4MhjMbIyZvQBQWlrqdTkiIlHFl8FQefIZFAwiIqHmy2Co7rPPPvO6BBGRqOL7YND9kkREQsvXwdCoUSP27NnjdRkiIlElZA/quRiVVz43btxYwSAiEmK+7DFUnnxu0qSJgkFEJMR8GQyV4uPjFQwiIiHm62CIi4vj4MGDlJWVeV2KiEjU8H0wAOzde0l37RYRkUvg62CovO22hpNERELHl8FQeUuM48ePAwoGEZFQ8mUwVM5KatOmDaBgEBEJJV8GQ6XY2FgSEhIUDCIiIeTrYADo0KGDgkFEJIR8HwxJSUns3r3b6zJERKKG74Ohc+fOlJSUeF2GiEjU8H0wpKSksGfPHl3kJiISIr4MhsrpqocPH6ZLly4459RrEBEJEV8GQ+V01RYtWpCSkgLAp59+6nFVIiLRwZfBUF2XLl0ABYOISKj4Phg6dOhAfHy8gkFEJER8HwwxMTEkJydTXFzsdSkiIlHB98EAFcNJ6jGIiIRGWARDSkqKgkFEJETCIhi6dOnCoUOHOHz4sNeliIhEPF8GQ/XrGABSU1MBKCoq8rIsEZGo4MtgqH4dA0BGRgYAW7Zs8bIsEZGo4MtgqKlr167ExsYqGEREQiAsgiE+Pp7U1FQKCgq8LkVEJOKFRTBAxXCSegwiIg0vrIKhqKiIkydPel2KiEhEC6tgOHHihK6AFhFpYGETDN26dQM0M0lEpKGFTTBUTlndvHmzx5WIiES2sAmGFi1a0KlTJ/Ly8rwuRUQkooVNMAD07t2bjz76yOsyREQimi+DoeYtMSr17t2brVu38uWXX3pUmYhI5PNlMNS8JUal3r17c/r0aTZt2uRRZSIikc+XwXAuvXv3BtBwkohIAwqrYEhOTqZ58+YKBhGRBhRWwWBm9OrVi40bN3pdiohIxAqrYADo378/GzdupLy83OtSREQiUtgFw+DBgykrK9NwkohIAwm7YBg0aBAAq1ev9rgSEZHIFHbB0LFjRzp27MgHH3zgdSkiIhEp7IIBKnoN6jGIiDSMsAyGwYMHU1JSwr59+7wuRUQk4oRlMFSeZ9BwkohI/QvLYOjTpw+NGzdm5cqVXpciIhJxwjIYmjRpQnZ2NsuXL/e6FBGRiBOWwQBw/fXXk5eXx4EDB7wuRUQkooRtMAwfPhyAd9991+NKREQiiy+D4VzPY6iuf//+JCYmajhJRKSe+TIYzvU8huoaNWrEkCFDWLZsWQgrExGJfL4MhtoaPnw4RUVF7Ny50+tSREQiRlgHw0033QTAokWLPK5ERCRyhHUwZGRk0LVrVxYuXOh1KSIiESOsg8HMGDNmDMuXL+fYsWNelyMiEhHCOhgAxowZQ1lZGbm5uV6XIiISEcI+GK677jqaN2+u4SQRkXoS9sEQHx/PyJEjWbBgASdPnvS6HBGRsBf2wQAwceJEDhw4oGsaRETqQUQEw8iRI2nevDmzZs3yuhQRkbAXEcHQpEkTvvOd7zBv3jyOHz/udTkiImEtIoIBYNKkSRw5coQlS5Z4XYqISFiLmGC4/vrradu2LTNnzvS6FBGRsBYxwRAXF8fEiRN58803OXTokNfliIiErYgJBoAf/OAHlJWVMWPGDK9LEREJWxEVDL1796Z///688MILOOe8LkdEJCxFVDAA/PCHPyQ/P58PPvjA61JERMJSxAVDTk4OTZs25cUXX/S6FBGRsBRxwZCYmMj3vvc9Zs6cyf79+70uR0Qk7ERcMADcd999HD9+nOeff97rUkREwk5EBkP37t0ZOXIk06ZN05XQIiIXKSKDAeCBBx6gtLRUF7yJiFykkAWDmaWY2ctmNjcUxxs+fDg9e/bkj3/8I6dPnw7FIUVEIkKtgsHMXjGzUjPbVGP9TWZWaGbbzOyh8+3DOVfsnLu7LsVeDDPj4YcfZvPmzcydG5IsEhGJCFabC8HMbAjwBTDdOdcjsC4W2ArcAOwG1gKTgFjgiRq7uMs5Vxp431zn3LjaFNevXz+3bt26Wn6UM506dYqePXtiZnzyySfExETsyJmICABmtt45168u+6jVN6Vz7h9AzRsQDQC2BXoC5cAs4Bbn3CfOudE1ltLaFmRm95jZOjNbV9fpprGxsfznf/6neg0iIhehLv8LnQTsqvZ6d2DdWZlZGzN7HrjGzB4+13bOuRecc/2cc/0uv/zyOpRXYfz48XTr1o3f/va3OtcgIlILdQkGO8u6c45LOecOOud+7Jzr6pyrOdTUYGJjY/n1r3/N5s2bdXM9EZFaqEsw7Aauqva6I7C3buU0jPHjx9OvXz8eeeQRvvzyS6/LERHxtboEw1ogzcy6mFk8kAMsqI+izGyMmb1w+PDh+tgdMTExPP300+zZs4ennnqqXvYpIhKpajtddSawGrjazHab2d3OuZPAvcD/AAXAHOdcfn0U5Zxb6Jy7p0WLFvWxOwC+8Y1vcOutt/KHP/yBffv21dt+RUQiTa2mq3qlrtNVa9q+fTvdunUjJyeH6dOn19t+RUT8ImTTVSNF165defDBB3n11VdZvny51+WIiPhSVPUYAL766iuysrIwMz7++GOaNGlSr/sXEfFSxPYY6vvkc3UJCQk899xzFBUV8cQTIZs1KyISNnwZDA1x8rm6ESNGcNttt/HEE0+wadOmC79BRCSK+DIYQuHpp5+mZcuWTJ48mfLycq/LERHxjagNhssvv5yXXnqJjz76iN/85jdelyMi4htRGwwAN998M3fffTd/+MMfWLVqldfliIj4gi9nJZnZGGBMamrqD4uKihr0WEePHqVXr14459iwYQOtWrVq0OOJiDSkiJ2V1NAnn6tLTEzkb3/7G3v27OGOO+7QHVhFJOr5MhhCbdCgQTz11FMsXLiQqVOnel2OiIinFAwB9957LxMmTOCRRx7hvffe87ocERHPKBgCzIyXXnqJ9PR0Jk6cSElJidcliYh4QsFQTWJiIn//+98pKytj9OjRHDlyxOuSRERCzpfB0JC3xLiQjIwM5s6dS0FBATk5OZw8eTLkNYiIeMmXwRDKWUlnM2LECP785z/z9ttv88tf/hI/TukVEWkojbwuwK9+9KMfsXXrVp5++mnatWvHo48+6nVJIiIhoWA4j6lTp1JaWspjjz1G69at+clPfuJ1SSIiDU7BcB4xMTG88sorfP755/zsZz+jZcuWTJo0yeuyREQalC/PMfhJXFwcc+bM4Rvf+Aa333478+fP97okEZEGpWCohYSEBBYsWEDfvn0ZP348c+fO9bokEZEG48tg8HK66rm0aNGCd955hwEDBpCTk8Ps2bO9LklEpEH4Mhi8nq56Ls2bN2fJkiVkZ2fz3e9+lxkzZnhdkohIvfNlMPhZYmIib7/9NkOHDmXy5Mk8++yzXpckIlKvFAyXoGnTpixatIjvfOc73H///UyZMkW36xaRiKFguEQJCQm8/vrr/PSnP2Xq1Knccccdena0iEQEXcdQB7GxsUybNo0OHTrw2GOPsXfvXubMmUObNm28Lk1E5JKpx1BHZsajjz7K9OnTef/99xk4cCCbN2/2uiwRkUumYKgnkydPZsWKFXzxxRcMGjSIt956y+uSREQuiS+DwY/XMdTGoEGDWLduHWlpadx88808/vjjOiktImHHl8Hg1+sYaqNjx46sXLmSnJwcHn30UUaPHs2BAwe8LktEpNZ8GQzh7rLLLuO1117jL3/5C8uWLeOaa65h1apVXpclIlIrCoYGYmb85Cc/YfXq1cTHxzN06FCmTp2qoSUR8T0FQwPr06cPGzZsYOzYsUyZMoURI0awc+dOr8sSETknBUMItGjRgtdff52XX36ZtWvX0rNnT2bMmKFHhoqILykYQsTMuOuuu8jLyyMrK4vJkyczYcIEDh486HVpIiJBFAwhlpKSwnvvvceTTz7Jm2++Sbdu3Zg1a5Z6DyLiGwoGD8TGxvLggw+ybt06kpOTmTRpEqNHj9a5BxHxBQWDh7Kysli9ejXPPPMMK1asoHv37jz77LOcOnXK69JEJIopGDwWGxvLfffdR35+PkOGDOH+++9n4MCBrF692uvSRCRK+TIYwvWWGHXRuXNnFi1axMyZM9m3bx/Z2dncfvvt7Nu3z+vSRCTK+DIYwvmWGHVhZuTk5FBYWMjDDz/M7NmzSU9PZ+rUqXrWg4iEjC+DIdo1a9aMxx9/nPz8fL75zW8yZcoUevTowRtvvKHZSyLS4BQMPpaamsqCBQtYvHgxcXFxjBs3juzsbFauXOl1aSISwRQMYWDkyJHk5eXx8ssvs2vXLoYMGcKYMWPIz8/3ujQRiUAKhjDRqFEj7rrrLoqKinjyySdZuXIlWVlZ3HnnnWzbts3r8kQkgigYwkxCQgIPPvggxcXFPPDAA8yZM4eMjAy+//3vKyBEpF4oGMJU69atmTp1KsXFxdx3333Mnj1bASEi9ULBEObatWvHU089dUZA3H777WzatMnr8kQkDCkYIkTNgJg3bx49e/bk29/+NitWrNA0VxGpNQVDhKkMiJ07d/L73/+etWvXMmzYMAYNGsQbb7yh+zCJyAUpGCJU69ateeyxxygpKeG5557j4MGDjBs3joyMDKZNm8bRo0e9LlFEfErBEOESEhL48Y9/TGFhIa+//jqtW7fm5z//OUlJSfziF79g69atXpcoIj6jYIgSsbGxjBs3jg8//JAPPviAm2++meeff56rr76am266iUWLFnH69GmvyxQRH1AwRKGBAwcyY8YMdu7cye9+9zs+/vhjRo8eXXXDvtLSUq9LFBEPKRiiWLt27fjVr35FSUkJs2bNon379kyZMoWOHTsyfvx43nnnHfUiRKKQL4MhGp/H4KW4uDgmTpzIypUryc/P59577+Xdd9/lxhtvJCUlhd///vfs2bPH6zJFJETMz/Pb+/Xr59atW+d1GVGprKyM+fPn8+KLL7Js2TJiYmIYNWoUd9xxB6NHj6ZJkyZelygiZ2Fm651z/eqyD1/2GMR7jRs3ZuLEiSxdupRt27bx0EMPsX79esaPH0/79u350Y9+xPvvv68L50QikHoMUmunTp1i+fLlTJ8+nXnz5vHll1/SpUsXbrvtNiZPnkxaWprXJYpEvfroMSgY5JJ88cUX/P3vf+fVV19l6dKlOOcYNGgQOTk5jB8/ng4dOnhdokhUUjCIL+zZs4e//e1vvPbaa+Tl5WFmXHfddUyYMIFx48bRrl07r0sUiRoKBvGdwsJC5syZw5w5c9i0aRNmxtChQ5kwYQK33norV1xxhdclikQ0BYP42ubNm5kzZw6zZ89my5YtxMTEMGzYMMaOHcstt9xCp06dvC5RJOIoGCQsOOfIz89n9uzZvPHGGxQUFADQp08fxo4dy9ixY+nRowdm5nGlIuFPwSBhqbCwkDfffJP58+fzwQcf4JwjJSWlKiSys7OJjY31ukyRsKRgkLC3b98+Fi5cyPz581m2bBnl5eW0bduWUaNGMWrUKL71rW/RqlUrr8sUCRsKBokoR48eZcmSJcyfP58lS5Zw6NAhYmJiyM7OrgqKrKwsDTmJnIeCQSLWqVOnWLNmDYsXL2bx4sVs2LABgKSkJEaOHMmoUaMYMWIEiYmJHlcq4i8KBoka+/btY8mSJSxevJh33nmHI0eOEBcXx+DBg7nhhhu44YYb6Nu3L40aNfK6VBFPKRgkKp04cYJ//vOfLF68mNzcXDZu3AhAixYtuP7667nhhhsYMWIEqampGnaSqKNgEAH279/P8uXLyc3NJTc3l507dwLQuXPnqpAYPnw4bdu29bhSkYanYBCpwTnHtm3byM3NZenSpSxfvpzK53pkZWUxbNgwhg4dypAhQxQUEpEUDCIXcPLkSdavX09ubi4rVqxg1apVfPXVVwD06NGjKiiGDh3K5Zdf7nG1InWnYBC5SOXl5axbt4733nuPFStW8P777/Pll18C0L17d4YNG8awYcMYMmQIV155pcfVilw8BYNIHZ04cYL169cHBcUXX3wBQGpqKtdee23VkpGRQUyMnm0l/qZgEKlnJ0+eZMOGDaxYsYJ//vOfrFq1iv379wPQqlUrBg8eXBUU/fv357LLLvO4YpFgCgaRBlZ5MnvVqlVVS+VNABs1akSfPn249tpryc7OJjs7Ww8oEs8pGEQ8cPDgQVavXl0VFGvXruX48eMAdOzYkYEDBzJgwAAGDhxI3759adasmccVSzQJq2Aws7HAt4ErgD8759650HsUDBIOysvL2bhxI6tXr2bNmjV8+OGHFBcXAxATE0NmZmZVUAwYMIDMzExdoS0NJmTBYGavAKOBUudcj2rrbwKeBWKBl5xzT9ZiX62APzrn7r7QtgoGCVcHDhxgzZo1VUGxZs0aDh06BMBll11G3759q4Kif//+dO7cWVdpS70IZTAMAb4AplcGg5nFAluBG4DdwFpgEhUh8USNXdzlnCsNvO8p4DXn3IYLHVfBIJHCOcf27duDgmLjxo2UlZUB0Lp1a/r06UOfPn3o27cvffr0oWvXrgoLuWghHUoys2TgrWrBMBj4jXPuxsDrhwGcczVDofL9BjwJ5Drnlp7nOPcA9wB06tSpb0lJSa0/jEg4KS8vJy8vj/Xr17NhwwbWr1/PJ598wokTJ4CKez9dc801VUHRt29f0tLSNGVWzqs+gqEuA51JwK5qr3cDA8+z/c+BEUALM0t1zj1/to2ccy8AL0BFj6EO9Yn4Wnx8PP3796d///5V68rLy9m0aVNVUGzYsIFp06ZV9SyaNWvGNddcUxUUvXr1IiMjg/j4eK8+hkSguvQYxgM3Oud+EHg9GRjgnPt5fRWnoSSRiovwCgoKgsLio48+qrpiOy4uju7du9OrVy969epFVlYWvXr10i0+opTXPYbdwFXVXncE9talGBE5U1xcHFlZWWRlZXHnnXcCFQ8yKiwsJC8vr2rJzc1l+vTpVe9r3759VVhULunp6ZoRJRdUlx5DIypOPg8H9lBx8vm7zrn8OhdlNgYYk5qa+sOioqK67k4kauzfv5+8vDw+/vjjqsDYvHlz1XmLxo0bk5mZWRUUmZmZ9OjRgyuvvFInuiNEKGclzQSGAW2Bz4BfO+deNrNRwDNUzER6xTn3f+pSTE0aShKpu/LycrZs2RLUu8jLy6u61QdUzIrq0aNHVVBkZmaSmZmpW5OHobC6wO1SKBhEGs5nn31Gfn4+mzZtCvpZ+fwKgCuvvPKMsMjMzKRly5beFS7npWAQkXrlnGPv3r1nhEV+fj7Hjh2r2i4pKakqLDIyMujWrRsZGRnqYfiAgkFEQuL06dPs3LnzjB5GQUFB1X2iANq0aRMUFJVLcnIysbGxHn6C6BGxwaCTzyLhoTIwCgoK2LJlS9VSUFAQdA6jcePGpKenB4VFRkYGV199NU2bNvXwE0SeiA2GSuoxiISvgwcPUlhYWBUUlaFRXFzM6dOnq7br1KkTGRkZpKenk56eTlpaGmlpaXTu3FlTay+BgkFEwk5ZWRnbtm0LCouCggKKioo4evRo1XZxcXGkpKSQlpYWFBjp6ekkJSXp1iDn4PUFbiIiF63yWorMzMyg9c45PvvsM4qKiigqKmLr1q1Vvy9dujToXEaTJk1ITU09IzDS0tJ0TUY9UI9BRHzv9OnT7NmzJygwKn8WFxdXXcAHkJiYSGpqKqmpqaSkpNC1a9eqn1dddVXEnwSP2KEknXwWkdo6efIkJSUlZ/Q0iouL2bFjR1BoNGrUiOTk5KCwqPw9JSUlIp62F7HBUEk9BhGpi1OnTrFr1y6Ki4vZvn37GT8///zzoO2vuOKKs4ZG165dadeuXVgMUSkYRETq4N///jfbt28/a2js2rWL6t+PCQkJdOnSheTk5KClcl2bNm18ERwKBhGRBlJWVkZJSUlQWOzYsYNPP/2UHTt2nNHbaNq06RmhUX0JVXBoVpKISAOpvCgvPT39rH///PPPKSkpYceOHWcsq1atOiM4mjVrdt7gaN26tS96HKBgEBG5JC1btqRly5b06tXrrH8/W3BU9jb+8Y9/cOTIkaDtmzVrRqdOnc65JCUlhexJfb4MhmqzkrwuRUTkktQmOKoHxs6dO6uWDRs2UFpaGrS9mdG+ffvzhkfr1q3rpXadYxAR8aGvvvqK3bt3BwVGzaX6RX9QcZ7j2LFjOscgIhKJEhISqq7qPhvnHAcOHAgKipKSEv70pz/V+djqMYiIRJD6mJWku1CJiEgQBYOIiARRMIiISBBfBoOZjTGzF6o/lFxERELDl8HgnFvonLunRYsWXpciIhJ1fBkMIiLiHQWDiIgEUTCIiEgQX1/gZmZHgUKv6/CJtsABr4vwCbXF19QWX1NbVLjaOZdYlx34/ZYYhXW9gi9SmNk6tUUFtcXX1BZfU1tUMLM63y5CQ0kiIhJEwSAiIkH8HgwveF2Aj6gtvqa2+Jra4mtqiwp1bgdfn3wWEZHQ83uPQUREQkzBICIiQXwZDGZ2k5kVmtk2M3vI63pCycyuMrN3zazAzPLN7L7A+tZmlmtmRYGfrbyuNVTMLNbMNprZW4HXUdkWZtbSzOaa2ZbAfx+Do7gtfhn497HJzGaaWZNoaQsze8XMSs1sU7V15/zsZvZw4Lu00MxurM0xfBcMZhYL/BkYCXQHJplZd2+rCqmTwP9yznUDBgE/C3z+h4Blzrk0YFngdbS4Dyio9jpa2+JZYIlzLgPoRUWbRF1bmFkS8Augn3OuBxAL5BA9bfFX4KYa68762QPfHTlAZuA9fwl8x56X74IBGABsc84VO+fKgVnALR7XFDLOuX3OuQ2B349S8Y8/iYo2+O/AZv8NjPWkwBAzs47At4GXqq2OurYws+bAEOBlAOdcuXPuc6KwLQIaAQlm1gi4DNhLlLSFc+4fwKEaq8/12W8BZjnnypxznwLbqPiOPS8/BkMSsKva692BdVHHzJKBa4APgSudc/ugIjyAKzwsLZSeAaYAp6uti8a2SAH2A/8VGFZ7ycyaEoVt4ZzbA/wR2AnsAw47594hCtuimnN99kv6PvVjMNhZ1kXdnFozawa8AdzvnDvidT1eMLPRQKlzbr3XtfhAI6AP8Jxz7hrgGJE7VHJegfHzW4AuQAegqZnd5m1VvnVJ36d+DIbdwFXVXnekopsYNcwsjopQeM05Ny+w+jMzax/4e3ug1Kv6Quha4GYz20HFkOL1ZjaD6GyL3cBu59yHgddzqQiKaGyLEcCnzrn9zrkTwDwgm+hsi0rn+uyX9H3qx2BYC6SZWRczi6fixMkCj2sKGTMzKsaRC5xzT1f70wLgjsDvdwBvhrq2UHPOPeyc6+icS6biv4PlzrnbiM62+Bewy8yuDqwaDmwmCtuCiiGkQWZ2WeDfy3AqzsVFY1tUOtdnXwDkmFljM+sCpAFrLrg355zvFmAUsBXYDjzqdT0h/uzXUdHV+xj4KLCMAtpQMdugKPCztde1hrhdhgFvBX6PyrYAegPrAv9tzAdaRXFb/BbYAmwCXgUaR0tbADOpOLdygooewd3n++zAo4Hv0kJgZG2OoVtiiIhIED8OJYmIiIcUDCIiEkTBICIiQRQMIiISRMEgIiJBFAwiIhJEwSAiIkH+P/g4cgB8BLzwAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "s = np.linspace(0, 100, 1001)\n", "plt.plot(s, erfcx(s), c=\"black\")\n", "plt.xlim(s[0], s[-1])\n", "plt.yscale(\"log\")\n", "plt.title(r\"$\\mathrm{erfcx}(s)$\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mathematical Reformulation\n", "\n", "### Strong Inhibition\n", "\n", "We have to consider three different cases; let us start with strong inhibitory input such that $0<\\tilde{V}_{\\mathrm{r}}<\\tilde{V}_{\\mathrm{th}}$ or equivalently $\\muV_{\\mathrm{th}}+\\frac{\\alpha}{2}\\sigma\\sqrt{\\frac{\\tau_{\\mathrm{s}}}{\\tau_{\\mathrm{m}}}}$. In this regime, we can change variables $s\\to-s$ to make the domain of integration positive again. Using $\\mathrm{erf}(-s)=-\\mathrm{erf}(s)$ as well as $\\mathrm{erfcx}(s)$, we get\n", "\n", "$$I(\\tilde{V}_{\\mathrm{th}},\\tilde{V}_{\\mathrm{r}})=\\int_{|\\tilde{V}_{\\mathrm{th}}|}^{|\\tilde{V}_{\\mathrm{r}}|}\\mathrm{erfcx}(s)ds.$$\n", "\n", "In particular, there is no exponential contribution involved in this regime. Thus, we get\n", "\n", "$$\\phi(\\mu,\\sigma)=\\frac{1}{\\tau_{\\mathrm{ref}}+\\tau_{\\mathrm{m}}\\sqrt{\\pi}\\int_{|\\tilde{V}_{\\mathrm{th}}|}^{|\\tilde{V}_{\\mathrm{r}}|}\\mathrm{erfcx}(s)ds}$$\n", "\n", "as a numerically safe expression for $\\tilde{V}_{\\mathrm{r}}<\\tilde{V}_{\\mathrm{th}}<0$.\n", "\n", "### Intermediate Regime\n", "\n", "In the intermediate regime, we have $\\tilde{V}_{\\mathrm{r}}\\le0\\le\\tilde{V}_{\\mathrm{th}}$ or $V_{\\mathrm{r}}+\\frac{\\alpha}{2}\\sigma\\sqrt{\\frac{\\tau_{\\mathrm{s}}}{\\tau_{\\mathrm{m}}}}\\le\\mu\\le V_{\\mathrm{th}}+\\frac{\\alpha}{2}\\sigma\\sqrt{\\frac{\\tau_{\\mathrm{s}}}{\\tau_{\\mathrm{m}}}}$. Thus, we split the integral at zero and use the previous steps for the respective parts to get\n", "\n", "$$I(\\tilde{V}_{\\mathrm{th}},\\tilde{V}_{\\mathrm{r}})=2e^{\\tilde{V}_{\\mathrm{th}}^{2}}D(\\tilde{V}_{\\mathrm{th}})+\\int_{\\tilde{V}_{\\mathrm{th}}}^{|\\tilde{V}_{\\mathrm{r}}|}\\mathrm{erfcx}(s)ds.$$\n", "\n", "Note that the sign of the second integral depends on whether $\\left|\\tilde{V}_{\\mathrm{r}}\\right|>\\tilde{V}_{\\mathrm{th}}$ (+) or not (-). Again, we extract the leading contribution $e^{\\tilde{V}_{\\mathrm{th}}^{2}}$ from the denominator and arrive at\n", "\n", "$$\\phi(\\mu,\\sigma)\t=\t\\frac{e^{-\\tilde{V}_{\\mathrm{th}}^{2}}}{e^{-\\tilde{V}_{\\mathrm{th}}^{2}}\\tau_{\\mathrm{ref}}+\\tau_{\\mathrm{m}}\\sqrt{\\pi}\\left(2D(\\tilde{V}_{\\mathrm{th}})+e^{-\\tilde{V}_{\\mathrm{th}}^{2}}\\int_{\\tilde{V}_{\\mathrm{th}}}^{|\\tilde{V}_{\\mathrm{r}}|}\\mathrm{erfcx}(s)ds\\right)}$$\n", "\n", "as a numerically safe expressions for $\\tilde{V}_{\\mathrm{r}}\\le0\\le\\tilde{V}_{\\mathrm{th}}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Noise-free Limit\n", "\n", "Even the noise-free limit $\\sigma\\ll\\mu$, where the implementation from ([Hahne et al., 2017](https://www.frontiersin.org/articles/10.3389/fninf.2017.00034/full#h10)) eventually breaks, works flawlessly. In this limit, $\\left|\\tilde{V}_{\\mathrm{th}}\\right|\\gg1$ as long as $\\mu\\neq V_{\\mathrm{th}}$; thus, we get both in the 'strong inhibition' and in the 'itermediate' regime $\\phi(\\mu,\\sigma)\\sim e^{-\\tilde{V}_{\\mathrm{th}}^{2}}\\approx0$ for $\\tilde{V}_{\\mathrm{th}}\\ge0$. Accordingly, the only interesting case is the 'strong excitation' regime $\\tilde{V}_{\\mathrm{r}}<\\tilde{V}_{\\mathrm{th}}<0$. Since also $\\left|\\tilde{V}_{\\mathrm{r}}\\right|\\gg1$, the integrand $\\mathrm{erfcx}(s)$ is only evaluated at $s\\gg1$. Using the only the first term of the [asymptotic expansion](https://en.wikipedia.org/wiki/Error_function#Asymptotic_expansion)\n", "\n", "$$\\mathrm{erfcx}(s)=\\frac{1}{s\\sqrt{\\pi}}\\sum_{n=0}^{\\infty}(-1)^{n}\\frac{(2n-1)!!}{(2s^{2})^{n}}$$\n", "\n", "leads to the analytically solvable integral\n", "\n", "$$I(\\tilde{V}_{\\mathrm{th}},\\tilde{V}_{\\mathrm{r}})=\\int_{|\\tilde{V}_{\\mathrm{th}}|}^{|\\tilde{V}_{\\mathrm{r}}|}\\mathrm{erfcx}(s)ds\\approx\\frac{1}{\\sqrt{\\pi}}\\int_{|\\tilde{V}_{\\mathrm{th}}|}^{|\\tilde{V}_{\\mathrm{r}}|}\\frac{1}{s}ds=\\frac{1}{\\sqrt{\\pi}}\\log\\frac{\\left|\\tilde{V}_{\\mathrm{r}}\\right|}{\\left|\\tilde{V}_{\\mathrm{th}}\\right|}.$$\n", "\n", "Inserting this into $\\phi(\\mu,\\sigma)$ and using $\\tilde{V}_{\\mathrm{th}}\\approx\\frac{V_{\\mathrm{th}}-\\mu}{\\sigma}, \\tilde{V}_{\\mathrm{r}}\\approx\\frac{V_{\\mathrm{r}}-\\mu}{\\sigma}$ yields\n", "\n", "$$\\phi(\\mu,\\sigma)\\approx\\begin{cases}\n", "0 & \\mu\\le V_{\\mathrm{th}}\\\\\n", "\\frac{1}{\\tau_{\\mathrm{ref}}+\\tau_{\\mathrm{m}}\\log\\frac{\\mu-V_{\\mathrm{r}}}{\\mu-V_{\\mathrm{th}}}} & \\mu>V_{\\mathrm{th}}\\end{cases}$$\n", "\n", "[as it should](https://neuronaldynamics.epfl.ch/online/Ch8.S3.html). Thus, as long as the numerical solution of the integral $\\frac{1}{\\sqrt{\\pi}}\\int_{|\\tilde{V}_{\\mathrm{th}}|}^{|\\tilde{V}_{\\mathrm{r}}|}\\frac{1}{s}ds$ is precise, the deterministic limit is also numerically safe.\n", "\n", "### Relevance of Noise-free Limit\n", "\n", "Let us briefly estimate for which values the noise-free limit becomes relevant. We have $\\left|\\tilde{V}_{\\mathrm{r}}\\right|>\\left|\\tilde{V}_{\\mathrm{th}}\\right|\\gg1$, thus the integrand $\\mathrm{erfcx}(s)$ is only evaluated for arguments $s>\\left|\\tilde{V}_{\\mathrm{th}}\\right|\\gg1$. Looking at the difference between $\\mathrm{erfcx}(s)$ and the first order asymptotics shown below, we see that the absolute difference to the asymptotics is only $O(10^{-7})$ for moderate values $\\left|\\tilde{V}_{\\mathrm{th}}\\right|\\approx O(100)$. Since we saw above that the noise free limit is equivalent to the first order asymptotics, we can conclude that it is certainly relevant for $\\left|\\tilde{V}_{\\mathrm{th}}\\right|\\approx\\frac{\\mu-V_{\\mathrm{th}}}{\\sigma}\\approx O(100)$; e.g. for $\\mu-V_{\\mathrm{th}}\\approx10$mV a noise strength of $\\sigma\\approx0.1$mV corresponds to the noise-free limit." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEOCAYAAABmVAtTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA0bUlEQVR4nO3dd1yV9fvH8dfFEhxoKk5UyoETcc/UcuTOzMy0Ybn6alpZlraXaX1N/VlaWq7MHJmZZjly5d57ooi5Byoiggh8fn+AfBEBQcZ9c7iej8d55LnndZ9PnPe51+cWYwxKKaUUgJPVBSillLIPDQWllFLxNBSUUkrF01BQSikVT0NBKaVUPA0FpZRS8TQUlFJKxdNQUEopFU9DIZsTkf0i0szqOpIjItNE5DOr68gusro9RcRXRHaKSKiIDLLT8kRkhIi8lorptohIlfSsS/2PhkI2ISJBIhIuItcTvEoYY6oYY1anY5ktMrjUHOt+Ps/E86SnPe/TW8BqY0w+Y8w4uyxPRLyA54GJqZh8FPDJ/a5L3UlDIXvpYIzJm+B1JqWJRcQlqwrL6HVlZe05XBlgf3oXkqC9MmR5QE/gT2NMeCqmXQg8IiLFM2C9OZ6GQjaX+Jdm3Pu3RWQPECYiLnHvT8ft0h8WkeYiMgMoDSyK2+t4K4llVxKR1SJyNe6wRsck1p14XTVEZEfcuuYA7gmmLyEiv4rIRRE5nvjwQlLLSzR+qIgci1v2ARF5ItH4pLZziIj8mmi6r0VkbIJ1DhGRPSISJiKTRaSoiPwVt5y/ReSBRDUOi1v/FRGZKiLuyX2eKX2GSc2TRHuWEpH5cZ9ZsIh8k9y2JvG/x73WvxJ4BPgmbv0Vkpg/2TZLor3uWl4K9ZcVkcsiUjPBei7J/w6dtQHWJKrlIRH5I266EBFZDmCMiQC2A62S+gxUGhlj9JUNXkAQ0OJew+Pe7wJKAR6AL3ASKBE33gcom9Iy48a5AkeBdwA34FEgFPBNYV1uwAng9bj5uwC3gM+I/QGyHfggbrqHgEDgseSWl0RNTwEl4pb1NBAGFI8bl+R2AsXjpisQN9wFuADUSrDOTUBRoGTcuB1ADSAXsBL4MFGN++JqLAisBz5Lpi1S+xkmbr8Wcf92BnYDY4A8xAZs45Ta9D7acDXQO5n/B1Jss6TaK+Hykqs/wfL7AAeB3MBSYFSCcReBOonqWQcMiFuuO9AowbhxwGir/04d4aV7CtnLgrhffFdFZEEK040zxpw0sbve0cR+uVUWEVdjTJAx5lgq1lUfyAuMNMZEGmNWAn8Az6SwrvrEfhGNNcbcMsbMA7bGTVcH8DLGfBK3vEDge6BbCsu7gzHmF2PMGWNMjDFmDhAA1I0bneR2GmPOAv8QGygArYFLxpjtCRb9tTHmvDHmNLAW2GyM2WmMuQn8RmxAJPRNXI2XgeFJfCZp/QyTU5fYEBxijAkzxkQYY9Ylt62ZsP7UtFmy7ZVC/QAYY74ntg03Exve7yaYtwCxAZZQWWIDwTluWesTjAuNm0elk4ZC9tLJGFMg7tUphelO3v6HMeYo8BrwEXBBRGaLSIlUrKsEcNIYE5Ng2Alif00nua64eU4bYxL2x34i7r9lgBIJQu0qsb9gi6awvDuIyPMisivB/FWBwnDP7ZwOPBv372eBGYkWfT7Bv8OTeJ83hRpPELvdSUntZ5icUsAJY0xUwoFpaNP0rj81bZZseyVXfyLfE9uOX8eF8G1XgHyJpu0BPA6ciTvMVzDBuHzA1ZQ2RqWOhoJjuuMhGcaYn40xjYn9IzfAF0lNl8gZoJSIJPx/pDRwOoV1nQVKiogkmgdivzyOJwi1Aib2CpW2KdV+m4iUIfYL5BWgkDGmALGHceLXlcJ2LgD8RKQq0B6Ymfxmp0qpBP8uTexnlVTtqfkMU2qDk0DpxOdWIMVtTev6U5KaNruv+gFEJC8wFpgMfJToS34PcMc5DmPMSmNMc6AyUJ3Yk9G3VSL2UJVKJw0FByex140/KiK5gAhif/lGx40+T+xx4qRsJvZY/Fsi4hp3ArADMDuF1W0EooBBEnvSuTP/O7yzBbgWd2LSQ0ScRaSqiNRJ5abkIfYL6GLcdr1I7C/Me26niT0ROQ/4GdhijPk3letMzgAR8Y77EnsHmBM3PPHnmZrPMKU22EJs0I4UkTxxJ7Qb3aNNE7qfNky8/vS0WZL1Jxj/f8B2Y0xvYDHwXYJxfwJNb78Rkc4iUj7uB0c+4AFiz2cQ9znUApansi6VAg0Fx5cLGAlcAs4BRYj9IgMYAbwXd2jgzYQzGWMigY7EXgVyCZgAPG+MOZTciuLm6UzsL7grxJ4Mnh83LprYLyR/4HjcMn8A8qdmI4wxB4CviA2e80A1Yk/ypmY7IfYQUjXuPnR0P34GlhF70jWQ2BPpkOjzTOVnmFIb3P7MygH/AqeI/Uzvta23509zGyazfn/ur82Sqx8ReZzY8zsvx00+GKgpIj3i3v8ItBURj7j3jYm9GimU2MAYGXeOhLhtXG3ucYm2Sh258/CvUo5JREoDh4Bixphr6VhOELFX1/ydUbWppInI58AFY8zYe0y3GehljNmXJYU5OL1BSDm8uGPqg4HZ6QkElbWMMXft/SQzXb3MriUn0VBQDk1E8hB7uOkEsYcrlFIp0MNHSiml4umJZqWUUvE0FJRSSsWz9TmFwoULGx8fH6vLUEqpbGP79u2XjDFe9zu/rUPBx8eHbdu2WV2GUkplGyJy4t5TJU8PHymllIpny1AQkQ4iMikkJMTqUpRSKkexZSgYYxYZY/rmz5+qu+mVUkplEFufU1BK5Ty3bt3i1KlTREREWF2Krbm7u+Pt7Y2rq2uGLldDQSllK6dOnSJfvnz4+PhwZy/s6jZjDMHBwZw6dYoHH3wwQ5dty8NHSqmcKyIigkKFCmkgpEBEKFSoUKbsTdkyFG6faL5yVU80K5UTaSDcW2Z9RrYMhdsnms+FC/N3nEL7Z1JKqaxhy1C4zdXZicFzd9Plu43sO617DUopldlsHQrliuTlyyf9CLoURodv1jH27yNWl6SUUncYN24clSpVokePHveeOJHw8HCaNm1KdHRST1OFyMhImjRpQlRUVHrLTDVbhwJA1zqlWPlmM3o29MG3aD4AbkXHEB2jh5SUUtYxxhATE8OECRP4888/mTlzZpqXMWXKFDp37oyzs3OS493c3GjevDlz5sxJcnxmsGUoJL6jOb+HKx92qEKbasUBmPRPIB2+XsfWoMtWlqmUcmA//fQTdevWxd/fn379+hEdHU1QUBCVKlWif//+1KxZk169ehEYGEjHjh0ZM2YMP/74I35+flSvXp3nnnsOgK1bt+Ln50dERARhYWFUqVKFfftinxw6c+ZMHn/88fh1Tp8+nVq1auHn58fDDz8MQKdOne4rcO6XrR+yU7t2bZNUh3hL9p3j40X7ORsSQSf/EgxrW4minu4WVKiUymgHDx6kUqVKALz22mvs2rUrQ5fv7+/P2LFj71nDW2+9xfz583F1daV///7Ur1+fJk2a8NBDD7Fhwwbq168P/K/jzvPnz9O5c2fWr19P4cKFuXz5MgULFgTgvffeIyIigvDwcLy9vRk2bBiRkZGULl2ac+fOARAaGkq9evXYtWsXbm5uXL16lQIFChAdHU2xYsW4ePFiknXe/qxuE5Htxpja9/v5ZMub11pXLUaTCoWZsOoYk/4JZPmB83zeuRqP+5e0ujSllANYsWIF27dvp06dOkDssf8iRYrQpEkTypQpEx8ICa1cuZIuXbpQuHBhgPhAAPjggw+oU6cO7u7ujBs3DoBLly5RoECB+GmcnZ0JDw/njTfe4IUXXqB27drxw93c3AgNDSVfvnyZtcnxsmUoAOR2c+HNx3zpUsubzxYfoEyhPABERcfg4mzLo2JKqTS61y/6zGKM4YUXXmDEiBF3DA8KCiJPnjzJzpPcvQOXL1/m+vXr3Lp1i4iICPLkyYOHh8cdN5/lzp2bffv2sWjRIvr27Uvv3r3p378/ADdv3sTdPWuOhmT7b0+fwnn44YU6+JcqAMD7v++jz4/bOHn5hrWFKaWyrebNmzNv3jwuXLgAxH6pnziR8mMKmjdvzty5cwkODo6f57a+ffvy6aef0qNHD95++20AHnjgAaKjo+ODISAggDx58tCtWzfat28fPzw4OBgvL68M7+MoOdk+FBLzKZSH9Ucv0Xz0GkYvO0x4ZNKXeimlVHIqV67MZ599RqtWrfDz86Nly5acPXs2xXmqVKnCu+++S9OmTalevTqDBw8G4Mcff8TFxYXu3bszdOhQtm7dysqVKwFo1aoV69atA2D48OH4+vpSs2ZNjh8/Hr+XsGrVKtq2bZuJW3snW55oFpEOQIdy5cr1CQgISPP850Ii+PzPgyzcfYaSBTwY94w/tcoUvPeMSinLJXXy1FHt3LmT0aNHM2PGjGSn6dy5MyNGjMDX1/eucZlxotmWewrpfZ5CsfzujHumBnP61qeIZy6K5fcAIEbvbVBK2UiNGjV45JFHUrx5rVOnTkkGQmax5Z7Cbcldkno/jDH0nLqVckXy8mqL8ni6Z83xOaVU2uSkPYX0yjF7CpkhMjqGEgXcmbL+OI+OWsO87ad0z0EppRLJMaGQy8WZEZ39+H1AI0oV9ODNX3bz5Hcb9ColpZRKIMeEwm1+3gX49eWGjHqqOreiY3ggjxuAds+tlFLkwFAAcHISutTyZtErjcmby4XIqBie/HYD0zcEERUdY3V5SillmRwZCrfdvvswJPwWHm7OfLhwP+2/XsfmwGCLK1NKKWvk6FC4zStfLn7qVY9ve9QkNCKKpydtYuCsnYRG3LK6NKWUylK27Psowc1rWblO2lQrTjPfIny75hjrj14it5stPx6llMo0ttxTSO/Na+nh4ebM4JYV+KVfA5ydhKs3Iuk8YT2rDl3I8lqUUtZ56aWXKFKkCFWrVr1rXL9+/Vi/fv1dw5999lnWrFmDiKTqZUe2DAU7cHKKbbBz1yK4euMWL07bSq9pWzkRHGZxZUqprNCzZ0+WLFmS5LjNmzff1X322bNnqVixIk2bNsUYk6qXHWko3EPFYp4sea0Jw9pUZFNgMC1H/8N/lx7SG9+UcnBNmjS545kItx08eJAKFSrc9QjN2bNn07179/j3R48excvLCx8fH/z9/SlYsCBly5bl2rVrQNJPWbMDDYVUcHNxol/Tsqx8sxnt/IoTcP56/J6EUipn+euvv2jdujXBwcEsXbo0fviJEyd46KGH4t+XK1eOxo0bM2PGDHbt2oWfnx8LFizA09OT0NBQvvjiCzZu3MiePXtYtGiRFZuSJD2TmgZFPd0Z87Q/t+LuZQi8eJ2PFh3g3baV8C2W+U9EUionenrixruGtfcrznMNfAiPjKbn1C13je9Sy5unapficlgk//lp+x3j5vRrkK56li5dytSpUwkKCuLNN9+kXr16nDlzBj8/v7um3b9/f/w5iUOHDsV3bJfcU9bsQPcU7oNr3JPdTgTfYM+pq7Qdt5aPF+0nJFwvYVXKkd24cYOrV69SokQJatWqxVNPPcUvv/zCr7/+SpcuXe6YNjw8nIiICB544AFOnjxJoUKFcHOL7UHh9lPWGjVqRN++fZkwYYIVm5Mk3VNIh0cqFmHVG80Yteww0zYEsXDXGYa1rUSXWt5Wl6aUw0jpl72Hm3OK4wvmcUv3nkFCq1at4pFHHol//+yzz9KzZ08aN26Mp6fnHdMeOHAgvgfTxL2ZBgQEUL58ebp168aBAwfueCyn1XRPIZ0eyOPG8CeqseiVxvgUzsPBs9esLkkplQGeeeYZGjRowOHDh/H29mby5Mnx5xNue+ihh4iOjqZZs2Z3zZ/w0JGHhwc7duzg0KFDQPJPWbODHPM8haxgjCEyOoZcLs5sOHaJBTtP81brihTOm8vq0pTKNuz8PIWaNWuyefPmO56XvG7dOho2bIiTU9b/xtbnKdiciJDLJfYytcPnQpm/4zSPjFrN1PXHtaM9pRzAjh077ggEgMaNG1sSCJnFllsiIh1EZFJISIjVpdy3Fxs9yJLXmuBfqgAfLzpAu3Ha0Z5Syv5sGQpWdnORkcoVycuPL9Vl4nO1CIuM0vMNSinb06uPMpmI8FiVYjSt4IVL3A1vC3ae5tSVG/R++CHcXZ3vsQSllMo6ttxTcETurs64xN3fsDXoMqOWHaHVmH9YcfC8xZUpZT92vgDGLjLrM9JQsMDwJ6rxU696uDoLvaZv48WpWzh+STvaUwrA3d2d4OBgDYYUGGMIDg7G3d09w5eth48s0rh8YZa81oTpG4IY+3cAxy9d58HCeawuSynLeXt7c+rUKS5evGh1Kbbm7u6Ot3fG3yir9ynYwNUbkRTIHXv7+5R1xymU142O1UvYtr91pZR96X0KDuB2IMTEGP7ce5ZXZ+/i6Umb9GolpVSW01CwEScnYU6/Bnz+RDUCzofSbtxaPvh9HyE3tKM9pVTW0FCwGWcnoXu90qx6sxnP1i/DrC3/cvpquNVlKaVyCA0FmyqQ241PHq/K+rcfpXKJ2N4Xx60IYPuJKxZXppRyZBoKNlfEM/aSs5DwW/y8+V+e/HYDb8zdzYVQ+3S1q5RyHBoK2UR+D1dWvNGUl5uWZeHu0zQftYYf1gbGPwVOKaUygoZCNpInlwtD21Rk6WtNqFnmAcb+HcCVG5FWl6WUciB681o29JBXXqa9WIdTV8Ipks8dYwxjlh+ha51SeD+Q2+rylFLZmO4pZFMiQqmCsQEQcOE6k9YG0mL0Gv7v7wAibkVbXJ1SKrvKslAQkTwiMl1EvheRHlm13pygQtF8rHijGc0rFmXM30doOWYNy/af075jlFJplq5QEJEpInJBRPYlGt5aRA6LyFERGRo3uDMwzxjTB+iYnvWqu5Us4MH4HjX5uXc9PFyd+XjRAW5G6UlopVTapHdPYRrQOuEAEXEGxgNtgMrAMyJSGfAGTsZNpsc3MknDcoVZPOhhfupdD3dXZyJuRfPNygCu34yyujSlVDaQrlAwxvwDXE40uC5w1BgTaIyJBGYDjwOniA2GdK9XpczV2Sm+x9W1AZcYtewIj45azW87T+khJaVUijLjy7kk/9sjgNgwKAnMB54UkW+BRcnNLCJ9RWSbiGzTrnPTr2XlovzWvyHF8rvz+pzdPPXdRvadzr7PvlZKZa7MCIWk+ns2xpgwY8yLxpj/GGNmJjezMWaSMaa2Maa2l5dXJpSX89Qo/QAL+jdiZOdqBF4K4+NF+3WPQSmVpMy4T+EUUCrBe2/gzP0s6MSJE7z33nsUKVKEokWLUqRIkfh/FyxYECcnPQqVWk5OQre6pWlTtTgh4bcQES6G3mT5gfM8XacUzk767AalVAY8ZEdEfIA/jDFV4967AEeA5sBpYCvQ3RizPw3L7AB0cHZ27mOMISbm7qtonJ2d8fLyuiMokgqP2//OjMfWZXffrj7GF0sOUaWEJx93rEJtn4JWl6SUSqf0PmQnXaEgIrOAZkBh4DzwoTFmsoi0BcYCzsAUY8zw+1l+7dq1zZYtWwgODubChQtcuHCB8+fP3/HfxP8OC0v6Wceenp7JhkbiYQUKFMgRTz0zxvDHnrN8/udBzoZE0LlGSYa2qRjfCZ9SKvuxNBQy2/08jjMsLCzF0EgYKsk9HNzV1RUvLy+KFi1K8eLFKVu2LOXKlYt/+fj44ObmllGbabkbkVGMX3WU7/85Tkf/Eox6qrrVJSml7pNDhsLtw0flypXrExAQkGnriYqKIjg4ONnQuHDhAqdPn+bo0aNcv349fj4nJyfKlClzR1Dcfj300EPZ9lBV0KUwPNycKerpzpHzoZwLiaBJBT3Zr1R24pChcNv97ClkBmMMFy9e5OjRo3e9AgICuHr1avy0IoK3t3eSgVG2bFny5Mlj3YakweC5u5i/4zStKhfl/faV4/tZUkrZm4aCDVy+fDnJwDh69CiJ77UoXbo01atXx9/fP/7l4+NjuyupbkZF88Pa43yz8igxxvBy07L8p1lZ3F2drS5NKZUCDQWbCwkJ4dixY/F7Ffv372fXrl0cPnw4/qoqT09P/Pz87giKKlWq2OIw1NmQcD7/8xCLdp9hyGO+DHiknNUlKaVS4JChkFXnFKx048YN9u3bx+7du9m1axe7du1iz5498ecunJ2dqVixYnxI1K5dmzp16lh2+GlzYDDVvPOT282FHf9ewdPdhXJF8llSi1IqeQ4ZCrc5wp5CWsTExBAYGBgfErcD49SpU0BsUNSoUYOGDRvGv0qVKnWPpWa8x8evZ//pEHo29OHVFuXJ5+6a5TUopZKmoZADXLp0iS1btrBhwwY2bNjA5s2buXHjBgDe3t53hIS/vz+urpn7JR18/Sb/XXqYOdtOUihPLoa2qUjnGiVx0ruilbKchkIOFBUVxZ49e+JDYsOGDZw4cQIADw8PGjZsSPPmzWnRogU1a9bE2TlzTg7vPnmVDxfuZ9fJq3z3bE1aVy2eKetRSqWehoIC4PTp02zcuJF169axcuVK9u7dC0CBAgVo1qwZzZs3p3nz5lSsWDFD79aOiTEsO3COVpWL4eQkbDwWjG+xfBTM4zg39ymVnThkKOSEE82Z7fz586xatYoVK1bw999/ExQUBECJEiVo3rw5rVu3pnXr1hQsmHH9HUXciqbRyJVExRjeaFWB7nVL4+Jsr0ttlXJ0DhkKt+meQsYJDAxkxYoV8a9Lly7h5OREw4YNadeuHe3ataNq1arp3os4cj6UjxbuZ8OxYCoVj+1or+6D2tGeUllFQ0GlWXR0NFu3bmXx4sUsXryYnTt3AlCqVCnatWtHhw4daN68Obly5bqv5Rtj+GvfOT774wBnQiJY+loTfIvp5atKZQUNBZVuZ86c4c8//2Tx4sUsX76csLAwPD09ad++PZ07d6Z169b3dX9EeGQ0yw6c43H/kgBsOHaJ2mUK4uaih5SUyiwaCipD3bx5kxUrVjB//nwWLFhAcHAw7u7utG7dms6dO9OpUyfy5Uv7r/6zIeE8/MUqShfMzQcdKtPMt0gmVK+UcshQ0BPN9hAVFcXatWuZP38+8+fP58yZM3h4eNCxY0d69OjBY489lqYuxFcfvsAniw4QeCmMFpWK8kH7ypQupB3tKZWRHDIUbtM9BfuIiYlh48aN/Pzzz8yZM4fg4GAKFixI165d6dGjBw0bNkxVp36RUTFMWX+ccSsCcBJh/dBHye+hd0QrlVE0FFSWu3XrFsuWLWPmzJksWLCA8PBwypQpQ/fu3enRowdVqlS55zLOhUSwJegyHauXAGBb0GVqlXkgRzzxTqnMpKGgLHX9+nUWLFjAzJkzWb58OdHR0VSvXp1nn32W559/niJF7n3uYMvxy3SduJGGZQvxUccqVCiqVyopdb80FJRtnD9/nrlz5zJz5kw2b96Mq6srjz/+OL1796Zly5bJHl6Kio5h1pZ/GbXsCNdvRvFCAx9ea1keT+1oT6k001BQtnTw4EF++OEHpk+fTnBwMGXKlKFXr1689NJLlCxZMsl5LodFMmrZYWZt+ZeyXnlZ9loT7WRPqTRyyFDQq48cx82bN1mwYAHff/89K1aswMnJibZt29KnTx/atm2Li4vLXfPsOx3C6avhPFalGNExhsPnQqlcwtOC6pXKfhwyFG7TPQXHcuzYMSZPnszUqVM5d+4cJUuW5OWXX6Zv377JnnuYu+0kb/+6h661SvFWa18K5b2/u6yVyinSGwp6a6nKMmXLluXzzz/n33//Zf78+VSqVIn333+fUqVK8fzzz7N169a75mlTtRi9Gz/IrztO0WzUaqauP05UdIwF1SuVM+iegrLUwYMHGT9+PNOnT+f69evUr1+fgQMH0qVLlztujDt6IZSPFx1gbcAlWlcpxnfP1bKwaqXsSw8fKYdw7do1pk2bxjfffENAQADFihVjwIAB/Oc//6FQoUJAbEd7S/efx9PDhYZlCxN2M4prEbcont/D4uqVsg89fKQcgqenJ4MGDeLQoUMsWbKEGjVq8P7771O6dGkGDhxIYGAgIkLrqsVoWLYwABNWH+XRUWsYv+ooN6OiLd4CpRyDhoKyFScnJx577DH+/PNP9u7dS9euXZk4cSLly5ena9eud5x36FanNE0qFOa/Sw/z2Jh/WHnovIWVK+UYNBSUbVWtWpWpU6cSFBTEkCFDWLZsGXXr1qVp06b88ccflCzgzsTnajOjV12cnYSXpm1j9LLDVpetVLZmy3MKep+CSkpoaCg//PADY8aM4eTJk1SsWJEhQ4bw7LPPgpML0zcE0aSCF77F8nE5LJJcLk7kyXX3fRBKOTKHPKdgjFlkjOmbP39+q0tRNpIvXz5ef/11jh07xsyZM3F3d6dXr16UL1+eyd9P5Lm6JeKf8Pbhwv00/2oNC3efwY4/fJSyK1uGglIpcXV1pXv37uzYsYPFixdTokQJ+vfvT9myZRk7diw3btygZ8MyFMrrxqBZO+k2aROHzl2zumylsgUNBZVtiQht27Zlw4YNrFixAl9fX15//XV8fHxYPmsSPz3vx/AnqnL4fCjtxq3j912nrS5ZKdvTUFDZnojw6KOPsnLlStatW0etWrUYNmwYZR96kKNLprGgtz8vNPChUbnYS1kvXb9JTIweUlIqKRoKyqE0atSIv/76i61bt9KkSRM++ugj/CqWI2LjTzhFhhETY+g9fRtPTFjPzn+vWF2uUrajoaAcUu3atVmwYAG7d++mTZs2jBw5kgcffJCPP/6Yrv5FOBsSwRMTNjDkl91cDL1pdblK2YaGgnJofn5+zJkzh71799KyZUs++eRjXulQj3aynZcalGLBrtM8Omq17jUoFUdDQeUIVapUYd68eezYsYNGjRrx0btD+aZPC572PEbLSl5UKh77vIaQG7csrlQpa2koqBylRo0aLFq0iE2bNuHv78/woYOYPbg9P0z8lquhN2g7bi39Z27n9NVwq0tVyhIaCipHqlevHsuWLeOff/6hfPnyDBw4EL9qVSkr51l56ALNv1rN1ysCiLilHe2pnEVDQeVoDz/8MKtXr2b58uWULF6UGUO7E73wI8rnieSr5UdoNeYfzoVEWF2mUlnGlqEgIh1EZFJISIjVpagcQERo0aIFGzZsYNGiRbjH3OCPYZ3w3D6N0u4RFPWMfQRo2M0oiytVKvPZMhS07yNlBRGhffv27Nq1i+nTpxNyZAszX21Lq1at+HvdFhp/sZIRfx3kuoaDcmC2DAWlrOTs7Mzzzz/PkSNHGDNmDDt37uSxVi1xPn+QiWsCaf7Van7fdVo72lMOSUNBqWTkypWL1157jWPHjvHOG4M4OO0dLvz8FhFXzvPq7F08PWmTPvFNORwNBaXuIX/+/Hz66accO3aMFzs048DYFwld8R3XgvYTEXYdQK9SUg5DQ0GpVCpWrBjjx4/n4MEDtCybh6Uj+1C2bFneHz2RRiNXMnPzCaK1oz2VzWkoKJVG5cqVY/bs2Wzfvp3q1aszcvhnXA46wLu/7ePxb9ax/cRlq0tU6r5pKCh1n2rWrMnff//Ngh+/w33jRC4u/JKDx0/z5LcbGTZ/r9XlKXVfNBSUSgcRoV27duzds4fRr3bn+ty3CNk4l/XLFhIUFIQxhqjoGKvLVCrVNBSUygAuLi7069ePo4f280qT0mz7cTi+vr489/aXtBy9mrUBF60uUalU0VBQKgPly5ePTz/9lCNHjtC9e3d+nTuLY4FBPDd5C32mb+Xk5RtWl6hUijQUlMoE3t7eTJ06lY0LplHu6ByurJnO8r0neeS/K5myLtDq8pRKloaCUpnI39+fFcuWMPuDF/FYNYqQA2v5evR/2bx5MzExRu+KVrajoaBUJhMR2rRpw97N/zC8XVnObviN+vXr06L/p3Sd8A9HL1y3ukSl4mkoKJVFXFxc6Nu3LwEBAbzzzjvs3LqZzUfP03L0Kj5asIfQCH3qm7KehoJSWSxfvnwMHz6cnb+Op/a5xVzbvZypG0/Q4LMlLD9wzuryVA6XZaEgIg+JyGQRmZdV61TKznx8fJg/azq/v9uVAlt/IPjfAN4a/CqbNm3Scw3KMqkKBRGZIiIXRGRfouGtReSwiBwVkaEpLcMYE2iM6ZWeYpVyRI0bN2bn378x4tFCnN+zlgYNGtDola94dcZGroRFWl2eymFSu6cwDWidcICIOAPjgTZAZeAZEaksItVE5I9EryIZWrVSDsbJyYmePV/gyJEjvPvuewQcOcKCvRep/+lfTF5zRDvaU1lGUrubKiI+wB/GmKpx7xsAHxljHot7PwzAGDPiHsuZZ4zpkpp11q5d22zbti1V9SnlSE6cOMHA90awJaoM7mX8KO4exQ99mlKlZAGrS1M2JyLbjTG173f+9JxTKAmcTPD+VNywJIlIIRH5DqhxO0CSma6viGwTkW0XL2rXACpnKlOmDAtnfMe8AQ+Tb+8v/Hv2As/3eIaNGzdaXZpycOkJBUliWLK7HcaYYGPMy8aYsintTRhjJhljahtjant5eaWjPKWyv8aNG7N70RSGN8zF2SO7adiwIXUGfsPI37cTGaUd7amMl55QOAWUSvDeGziTvnKUUok5OTnR84XYZ0a/894HnDp/ie82nqP2+/NZuufkvRegVBqkJxS2AuVF5EERcQO6AQszoigR6SAik0JCQjJicUo5hLx58zL804/ZPOolKp5fSXDwZfr9vIc2ny/gXEi41eUpB5HaS1JnARsBXxE5JSK9jDFRwCvAUuAgMNcYsz8jijLGLDLG9M2fP39GLE4ph1K6dGmWTP2Kn5+rQp6jy9l78jJPdGzP9u3brS5NOYBUX31kBb36SKmUxcTEMHnqNN57ZxgXLwVTa9C3DH7yYbo18kUkqdN+ytGl9+ojW4aCiHQAOpQrV65PQECA1eUoZXshISG88+kXLLhcAtciD1LKNYyJ/VpQ2bug1aWpLOaQoXCb7ikolTYHDh7ipc+ncLpQLZzdPHi0tAvf9G1FbjcXq0tTWcTK+xSUUjZTuVJFNs34kq8eyYvTic0s2RnEU12eRPe4VWppKCjlgLo+3o5DMz5kUMVw1q5eRVX/WjR6ezqbDutV4ypltjx8pOcUlMo4586d45UPvmSzW02c8uSn9gORfNe/LV6eHlaXpjKBnlNQSqXK6g2bGTD+D66XqI1TzC161irM+92a4uSkVyk5Ej2noJRKlWYN67F3xscM8r1O9MXjTFi4jhdfepGzZ89aXZqyEQ0FpXIQJycn3ujdg33f9OM5nwhmz5pFxZoNeOyDmZy4oD0IKJuGgnZzoVTm8vT0ZNTIz9i/fz/+zTtxMMyDpl+uYOC3i4i4FWV1ecpCtgwF7eZCqaxRrlw51vw0lpFN8uB8MYBFJ5yoNmQ2M/7WLjNyKluGglIqa3V//DEO/TCYJx44TUR4OIO/ms4bb7yB7q3nPBoKSikAXF1dGfN2XzZ/3JEOpaIYM2YMvg1b0+OLOVyP0GdF5xQaCkqpO5QoVpQpkyawbds2ilR7mPVX8uI3dB6j563Gzpewq4yhoaCUSlLNmjXZPesL+pULI/pGCOO2hVFj8FTW7Q20ujSViWwZCnr1kVL2ICIM692VPaOepY4c43JMbjoP/JBRo0YRGamHlByRLUNBrz5Syl7ye+bjlxGD+K1PDep4XmPIkCFUafUM70/9i5gYPaTkSGwZCkope6pV1ZfFC39n8eLFRJf0Z8bhGKq+Po3FmzLkoYvKBjQUlFJp1rZtWw5MHkLz3Ce5HuPKgAVBNB82lX8vXLG6NJVOGgpKqfvi7u7O5A9eZsXgppQMPcjRqII07vwSc+fO1auUsjENBaVUupR/sBTrx7/J2FYPUDAsiKeffpo6XQfy899brS5N3Qdbdp2tz1NQKnuKjo5m0vc/8PkOg3PBUpSIOsvkge2pVKaY1aXlGPo8BaWU7Zw5f5GeI3/ikFMZBGhW9CbfDnqS3O5uVpfm8PR5Ckop2ylR1ItlY15n8pNlyH0tiDWXPWnQ9ik2btxodWnqHjQUlFKZpmXDWhyY9Dqv+V4n+Mh2GjZsSMv/fMLmA8etLk0lQw8fKaWyxPXr1/lo+BfMDa+KuLpTK+9VJg9+igfy5ba6NIeih4+UUtlC3rx5GTXiU+a9VJ0Hrh1jR3hharz3G5/++JdewmojGgpKqSxV168iuyYO5k1/kIhQvt8dTvunehAYqB3t2YGGglLKEq90a8f+sT3pWvAka5YspHLlyjz1zjecvnTV6tJyNFueU9D7FJTKWc6cOcOAdz9nh1criAznyQqufNmvEy7O+rs1rRzynIL2kqpUzlKiRAl+m/oNnzfLj+uNS8z/NxeVX53KLyv1QpOsZstQUErlTD3aNuXQhH50KBxMhLgxeP5BBgx6lStXtKO9rKKhoJSyFRcXF75+83k2vNuapnKA78Z/QwXfSvT/aiYRkbesLs/haSgopWzJu5gXP309gh07dlCmUQf+vFiAKoN/YtLv/1hdmkPTUFBK2Vr16tXZ8utEnve5QZQRPt8YSu1XJ7LjcJDVpTkkDQWllO05OTnxyctPsWN4Z/zkBBddvOgw4jd9VnQm0FBQSmUbhQp4snBEf35+rhKVr+9iyJAhVKtVjy9n6l3RGUVDQSmV7TTyr8TyedNZvHgxUT4NmbA3Br+BE1m946DVpWV7GgpKqWyrbdu27J4zikYeZwlxLcQLsw7T9p3vuXA5xOrSsi0NBaVUtpY3twczP+zNnwPqUjTiJAdiSlCn35f6rOj7pKGglHIIVcuVYcvXA/moUR7yn97I008/zcNtnmDBmu1Wl5at2DIURKSDiEwKCdFdQKVU2vTs0Iyda5fz7bffcjxPZV5dfJqmb3zL8TMXrC4tW7Blh3i36UN2lFLpceLMBV4YNZfjLqUhMpw23lGMe/Vp3FxdrC4t0zhkh3hKKZURypQowurRr/B/bYvjfjOYJZcK4N/tTTZt2mR1abaloaCUcnidmtXh4PiXeaZ0GCHb/6BBgwY81ftV9h87aXVptqOHj5RSOUpoaCifDR/OjHPFcC3ozcMFrzPxjWfI45HL6tIyhB4+UkqpNMiXLx9fjBzJTwOaky/8HOuue1F1yM/835xlVpdmCxoKSqkc6dE61dg7YQD9KkVjDIzZeYvmLwzmxIkTVpdmKQ0FpVSOJSIMe6Eje754moa5TrLx1++pVKkSr340iivXwqwuzxIaCkqpHM8zb25+/vhlDh3YT9sOHZl/qRg13v2VT6YsJCYmxuryspSGglJKxSldujTz5szmnWZFcYqOZMoRZ6oOnMjyzXutLi3LaCgopVQi/+nSkv2jn+MRz4tcdytE73mB9HvrI0JDQ60uLdNpKCilVBI83HMx9Z2eLH+1ERXC9zPpvx/j6+vLf7+fRXS04x5S0lBQSqkU+PqUZPk377Bp0yaKl63M14dyUWngD8xbsdnq0jKFhoJSSqVCvXr12LxqKZ28I7jpmpc3ll3g4dcncPTkOatLy1AaCkoplUouLs6Me70764a1olz0v/zr6s2jo1YxdsJEoqOjrS4vQ2goKKVUGpUqVpgVowYwoaM3Xue28PqAl6lTpw7zlv5jdWnplmWhICKdROR7EfldRFpl1XqVUiqztGtck22zRzNnzhwumry8uSqUOoMmsDsg+94VnapQEJEpInJBRPYlGt5aRA6LyFERGZrSMowxC4wxfYCewNP3XbFSStmIiNC1a1d2rvqDanKKC27F6fjdNp75ZArXb0RYXV6apXZPYRrQOuEAEXEGxgNtgMrAMyJSWUSqicgfiV5FEsz6Xtx8SinlMAoX8GTRiH781KMi+cPPsvFGUaoNnsZfS5ZYXVqapCoUjDH/AJcTDa4LHDXGBBpjIoHZwOPGmL3GmPaJXhck1hfAX8aYHRm7GUopZQ8P16jEngmvMKAqOB1dS9s2bXi8Uyc27z1idWmpkp5zCiWBhE+oOBU3LDkDgRZAFxF5ObmJRKSviGwTkW0XL15MR3lKKWWdIc+248DiKYwYMYK1J27Qdfp+OrwzieAQe98VnZ5QkCSGJfvEHmPMOGNMLWPMy8aY71KYbpIxprYxpraXl1c6ylNKKWvlypWLoUOHsnL2JLxunmZvTElqvreA9yfNt21He+kJhVNAqQTvvYEz6StHKaUcj7+vD9u+foV363vgHB3BjMBcVO/7FQcOHLC6tLukJxS2AuVF5EERcQO6AQszoigR6SAik0JCQjJicUopZQt9Oj3KwbE9aVEgmHPbllC9enVeHfwmJ88HW11avFQ9o1lEZgHNgMLAeeBDY8xkEWkLjAWcgSnGmOEZWZw+o1kp5aguXbrEu+++y5y9IeSv05FODwpfvdIVFxfndC03vc9oTlUoWEVDQSnl6H75exPv/LqLW/lL4XrtFJ928qNbqwb3vbz0hoItu7nQw0dKqZziqRb1OfRNHzqXCCXSJS9vr7hEy0FfYtXVl7YMBWPMImNM3/z581tdilJKZTpnZ2dGD+rGxvfb4GtOsX7+VCpUqMDoceOJuBmZpbXYMhSUUionKuFVkGVf/ofty+dTu3Zthv91mMqDZ/Dd/BVZVoOGglJK2UylSpVYtmwZg7o8Soy4MHJLBLUGjmfHocBMX7ctTzSLSAegQ7ly5foEBARYXY5SSlnmckgoL/13FjsjCoMxtM5/jnFv9yZXrlxJTu+QJ5r1nIJSSsUqmD8fCz7ry+znq1Iw/DSTR31E1apVmb9wcaasz5ahoJRS6k4N/Cqwc8Ig/vr1Z5xdXPjPz7upOmACq7dn7F3RGgpKKZWNtGzZkp07d9GiSnFCc3nxwqwjtB02kQtXMuYSfluGgt6noJRSyfNwz8XM91/kj/71KBp5mgPGmzofLmL01F/SvWxbhoKeU1BKqXurVq40W8YN4KPG+XALO8+b/Z5P9zJtGQpKKaVSr2f7Jhya9Brffj023cvSUFBKKQfg7OxMv3790r0cDQWllFLxNBSUUkrFs2Uo6NVHSillDVuGgl59pJRS1rBlKCillLKGhoJSSql4GgpKKaXiaSgopZSK52J1AUm5/TwFIEJE9ltcTn4gIy+Dut/lpWW+1Eyb0jTJjUvL8MLApXvUkBXs0H5Z2XYpjc9u7WeHtkvrfJn1t5fcuKSG+d5j/Skzxtj2BWyzQQ2T7LC8tMyXmmlTmia5cWkZboe2s0v7ZWXbOVL72aHtsrr90jouM9pODx/d2yKbLC8t86Vm2pSmSW5cWofbgR3aLyvbLqXx2a397NB2aZ0vs/72khuX4W1ny8dx3iYi20w6HiunrKNtl71p+2Vf6W07u+8pTLK6AHXftO2yN22/7CtdbWfrPQWllFJZy+57CkoppbKQhoJSSql4GgpKKaXiZctQEJGHRGSyiMyzuhaVOiKSR0Smi8j3ItLD6npU2ujfXPYlIp3i/u5+F5FW95o+y0NBRKaIyAUR2ZdoeGsROSwiR0VkaErLMMYEGmN6ZW6l6l7S2JadgXnGmD5AxywvVt0lLe2nf3P2ksa2WxD3d9cTePpey7ZiT2Ea0DrhABFxBsYDbYDKwDMiUllEqonIH4leRbK+ZJWMaaSyLQFv4GTcZNFZWKNK3jRS337KXqaR9rZ7L258irK87yNjzD8i4pNocF3gqDEmEEBEZgOPG2NGAO2zuESVSmlpS+AUscGwi2x62NLRpLH9DmRxeSoFaWk7ETkIjAT+MsbsuNey7fLHWZL//YqE2C+QkslNLCKFROQ7oIaIDMvs4lSaJNeW84EnReRb7Nutgkqm/fRvLltI7m9vINAC6CIiL99rIXbpJVWSGJbsXXXGmGDgnhunLJFkWxpjwoAXs7oYlWbJtZ/+zdlfcm03DhiX2oXYZU/hFFAqwXtv4IxFtaj00bbM3rT9sq8MaTu7hMJWoLyIPCgibkA3YKHFNan7o22ZvWn7ZV8Z0nZWXJI6C9gI+IrIKRHpZYyJAl4BlgIHgbnGGKsfrqPuQdsye9P2y74ys+20QzyllFLx7HL4SCmllA1oKCillIqnoaCUUiqehoJSSql4GgpKKaXiaSgopZSKp6GglFIqnoaCUkqpeBoKSiml4v0/mYGMTEYNZicAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "s = np.linspace(0.1, 100, 1000)\n", "plt.plot(s, erfcx(s), c=\"black\", label=r\"$\\mathrm{erfcx}(s)$\")\n", "plt.plot(s, 1 / (np.sqrt(np.pi) * s), ls=\"--\", label=r\"$1/\\sqrt{\\pi}s$\")\n", "plt.xlim(s[0], s[-1])\n", "plt.xscale(\"log\")\n", "plt.yscale(\"log\")\n", "plt.legend()\n", "plt.title(r\"First order asymptotics of $\\mathrm{erfcx}(s)$\")\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEOCAYAAACNY7BQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvpUlEQVR4nO3deXQUVfbA8e8NSYQBHEYdGFlkB0UgAUJABgkgMAmSoDgaFBQRWZRFcUFcQcUFfzMIiAooyCDK4kpAEBxQwmYkhoRVJLIIbqDiirLe3x/dMJ02CZ10J1VJ7uecPif9qurVrXrVfVP1ql6LqmKMMcacEuZ0AMYYY9zFEoMxxpgcLDEYY4zJwRKDMcaYHCwxGGOMycESgzHGmBwsMRhjjMnBEoMxxpgcynRiEJE9ItLF7XWWVCLSWEQ2isjPIjIikOkislVEOhZ3rIEQkVkiMs7pOEqS4mzPMx1vTtYnIk+IyO0BzPeRiFwczLpCoUwkBhH5QEQOichZTsfiqwwkkVHAB6paWVUnBzJdVS9W1Q8KuqIysC+LXUH3aW7zF7Y9C+lMx5sj9YnIX4EbgGkBzP4v4JHCritUSn1iEJE6wKWAAknORuMMEQkPpKygdQSgNrA1iOnBrr9QQr2u4oy9jAv4eMqPT3uFpD7gRmCJqv4WwLwpQCcROT8E6y08VS3VL+AhYC0wAVjsN20PcC+wDTgEvASU9067B/gC+BnYAVzms9xFwAfAD3gOnCS/Ort4/1aggc+0WcA4798vAyeB34BfgFHe8urAG8BBYDcwIp9ty3Nebxz3AJuAI0B4HmVn2pYc8+cSQ67LAyuBE8Dv3u1r5LdcrtP99l9u8f6hXfLal4HGms+6WgAZ3nXNB+adar/C7P9c4hkNfOatfxtwpc+03LbzbuANvzqeASb6rPNu7zp/BWYA1YCl3nr+C/wlwOP/D/v0DPsvr+PZtz1rAW9699d3wJQzfdZCdbwVor3+UF8+8dcHvgda+qznW6CjT2x9/WKpByz2zvcj8J7PtPeAfo5+bzq58mLZQMgGbgVaAceAan4HwxZvg5+DJ4GMAxoD+4Dq3vnqAPW9f0d467wPiAQ6ew/oxrl8EPJMDP7zet+HAR/jSWaR3oNnF/CPXLYr33m9dWd6t61CbmUBbkuOOvxiONPyHwA359M2f5jOHxODb7z5tUuOfVmIWP3XFQnsBUZ6l/0nnuNnXGH3fy4xXY3nSyQMSMbzZX5+XtvpnfYrUMVbHg4cAFr5rPNDPMmghndaBp4EdxaeL6gxZzr+82iLfPdfXm1wqgwoB2QBTwMVgfJA+/zaNJTHW2Hay7e+vOL3qX8gsB34E7AM+JfPtINAa7941gBDvfWWB/7uM20yMMHR700nV17kG+c58I4B53nffwKM9Dtoh/i8747nP7gG3g9VFyDCr85Lga+BMJ+yucDYXD5MBU0MbYDP/dZ3L/BSLtuW77zeum/ym56jLMBtucl/3QVYPs8Pal7T+WNi8I03v3bJsS8LEav/ujoAXwLiU7aO/yWGAu//AI7XTKDnGbZzKTDQ+3cPYJvfPujj8/4N4Hmf98OBt890/OfRFvnuv7zagP8lhkvwfEGG+03Pc1tDebwVpr3ImRhyjd9v/hRgM56zjrN8yo8BF/rN+xUwAojMpZ7HgJkFOXZC/SrtfQz9gOWq+q33/aveMl/7fP7ei+c/l2zgdmAscEBE5olIde881YF9qnrSb7kaIYi3NlBdRH449cLz31G1Qs67L5flfMsC2Zbc6ijI8sE6vf4ztMuZFHRbqwNfqPeT6jP/KYXd/6eJyA0ikumzfFM8/8Tkt53/Afp6/+6L5xKOr298/v4tl/eV/Ob/w/GfR7jBtnUtYK+qHvctLGCbBhNDsO2Va/x+XsDThs+o6hGf8kNAZb95++D5J+BLEZkhIuf4TKuM51KZY0ptYhCRCsA1QJyIfC0iX+O5LBAlIlE+s9by+fsCPP8loqqvqmp7PAeUAuO983wJ1BKRML/lvsgljMN4Ti1P+ZvfdPV7vw/YrapVfF6VVbV7LnUHMq9//f5lgWxLbnUUZPlg5Vh/Pu2SX5xQ8G39CqghIuI3/ymF3f8AiEhtPF8kw4BzVbUKnss6cobtfBtoLiJN8ZwxvJL/Zp9Rrsd/LvEHe6zsAy7IrSM+n231F8zxFlR75Rc/gIhUAibi6dcZ6/dFvwlPH8X/VqS6UlUvA5oAUXg6qE+5CM9lK8eU2sQAXIGn86gJEO19XQSsxnPr2ClDRaSmtyHvA+Z771/u7L299Xc8/2md8M6fhuc67ygRifDeo52Ip2PSXyZwnYiUE5F4IM5v+jd4rnWe8hHwk4jcIyIVvMs1FZHWudRdkHnzUpBtKYrlC+QM7eK/L4ONdT1wHBghIuEi0guI9Zke7P6viOeL6KB32/rj+W8z3+1U1d+B1/Gc/X6kqp8HuL68/OH495nmu08D2X/5tcFHeJLtkyJSUUTKi8jfz9Cm/oI53oJtr1zj95k+CfhYVW8G3gGm+kxbgs9nX0R6iUhD7z8dlYG/4PmuwLsfWuHpgHZMaU4M/fBcP/xcVb8+9QKmAH18Mv+rwHI8HVG78HQ+nwU8ieeOga+Bqng+NKjqUTy3vSZ4pz8H3KCqn+QSw214Dtwf8Jw6vu03/QngAe+p7V2qesI7fzSeuya+BV4E/uxfcUHmzUsBtyXkyxdCnu2C374MNlbv/L3w/Cd3CE/n8Js+04Pa/6q6Dfg3ngT0DdAMT+fvmbYTPJeTmvHHy0iFkdvxf8rpfYrneviZ9l+ebeCzvxoAnwP78ezTM22rbx2FPt5C0F55xY+I9ATigSHe2e8AWopIH+/72UB371UM8PR9rsLTcb4EeFJVV3qnJeF5dsL3zK3YSc5LqMYYtxORC/DcSPE3Vf0piHr24Olc/W+oYjO5E5HHgQOqOvEM86UBA1R1S7EElgd78MaYEsR7ff0OYF4wScEUL1XN9Swol/naFHUsgbDEYEwJISIV8Vx22ovn0oUxRaLYLiWJSD3gfuDPqvrPYlmpMcaYAguq81lEZorIARHZ4lceLyI7RCRbREYDqOouVR0QzPqMMcYUvWDvSpqF3ymtiJQDnsVz50AT4FoRaRLkeowxxhSToPoYVDVVPKOX+ooFslV1F4CIzMPzhN+2gtZ/3nnnaZ06/tUbY4zJz8cff/ytqv61sMsXRedzDXI+Wr4faCMi5+IZA6SFiNyrqk/ktrCIDAIGAVxwwQWkp6cXQYjGGFN6icjeM8+Vt6JIDJJLmarqd/zvAZA8qep0EfkKSIyMjGwV8uiMMcbkqyiefN5PzvFXapJz/BVjjDEuVhSJYQPQUETqikgk0BvPcLQBU9VFqjroz38OeHQHY4wxIRLs7apz8Yz10lhE9ovIAO+wtMPw/FjFdmCBqhbo5/FEJFFEpv/444/BhGeMMaYQXD1WUkxMjFrnszHGFIyIfKyqMYVd3pWjq9oZgzHGOMeVieFUH8NZZ52Fm89ojDGmNHLlIHoikohn7HMaNmxIUlISiYmJtG/fnoiICIejM8aY0s3VfQy1a9fWpk2bsmLFCo4cOUKVKlXo3r07SUlJxMfHY3ctGWPMHwXbx+DKxHDqjKFBgwYDd+7cyS+//MJ7773HokWLWLx4MQcPHiQ8PJy4uLjTZxN169Z1OmxjjHGFUpkYTsntrqQTJ06QlpbGokWLSElJYds2zxBMTZs2PZ0kYmNjCQtzZfeJMcYUuTKXGPxlZ2ezaNEiFi1aRGpqKidOnKBatWr06NGDxMREunTpQsWKFYspYmOMcV6ZTwy+Dh06xLvvvktKSgpLly7lxx9/pHz58nTp0oWkpCR69OjB+eefX4QRG2OM80plYvDvYyiMY8eOsXr1alJSUkhJSWH37t0AtG7d+vQlp+bNmyOS25h/xhhTcpXKxHBKqJ58VlW2bt16ul8iLS0NVeWCCy44nSQ6duxIZGRkCKI2xhhnWWIohG+++YZ33nmHlJQUli9fzm+//UblypWJj48nMTGR7t27c+6554Z8vcYYUxwsMQTpt99+Y+XKlaSkpLBo0SK++uorwsLCaN++PYmJiSQlJdGoUaMijcEYY0LJEkMInTx5koyMjNP9EllZWQA0b96c5ORkkpOTqV+/frHFY4wxhVEqE0MoOp9DYe/evSxcuJD58+ezbt06AFq1akVycjLXXHMNtWvXdiw2Y4zJS6lMDKe4adjtzz//nNdee4358+ezYcMGANq2bUtycjJXX301NWrUcDhCY4zxsMTggF27drFgwQLmz59PZmYmIkL79u3p3bs3V111FdWqVXM6RGNMGWaJwWE7duxg/vz5zJ8/n23bthEWFkanTp1ITk6mV69edneTMabYlZjEICIVgeeAo8AHqvrKmZYpCYnB15YtW04niZ07dxIeHk6XLl1ITk7miiuuoEqVKk6HaIwpAxz9BTcRmSkiB0Rki195vIjsEJFsERntLe4FvK6qA4GkYNbrVk2bNuXRRx9lx44dZGRkcOedd/LJJ5/Qv39/qlWrRlJSEnPnzuXw4cNOh2qMMXkKdgjSWUC8b4GIlAOeBRKAJsC1ItIEqAns8852Isj1upqI0KJFC5588kl27dpFWloaQ4cOJSMjg+uuu47zzz+fgQMHsnbtWvuFOmOM6wSVGFQ1FfjerzgWyFbVXap6FJgH9AT240kOQa+3JBERYmNjmTBhAp9//jnvv/8+V155JXPnzqV9+/Y0atSIcePGsXfvXqdDNcYYoGi+oGvwvzMD8CSEGsCbwFUi8jywKK+FRWSQiKSLSPrBgweLIDznhIWF0bFjR2bNmsXXX3/NrFmzqFmzJg8++CB16tThsssuY/bs2fz6669Oh2qMKcOKIjHkNlypquqvqtpfVW/Jr+NZVacDDwMZpXlQu0qVKtGvXz/ef/99du/ezcMPP8yePXvo168ff/vb3+jfvz+rVq3i5MmTTodqjCljiiIx7Adq+byvCXxZkApUdZGqDiorv+lcp04dHnroIbKzs0lNTSU5OZk33niDjh07Ur9+fcaOHXt62HBjjClqRZEYNgANRaSuiEQCvYGUglQgIokiMv3HH38sgvDcS0S49NJLefHFF/nqq694+eWXadCgAY888gj169fnH//4B2+99RbHjh1zOlRjTCkW7O2qc4H1QGMR2S8iA1T1ODAMWAZsBxao6tbgQy1bKlasSN++fXnvvffYs2cPY8aMYdu2bfTq1YvatWvz4IMPWoe1MaZI2JPPJcjx48dZsmQJ06ZNY+nSpQAkJCQwePBgunfvTnh4uMMRGmPcwNEH3IpKWb2UdCbh4eEkJSXxzjvvsHv3bu6//342btxIz549qVOnDmPHjmX//v1Oh2mMKeHsjKGEO3bsGIsXL2batGksX74cESExMZERI0bQqVMn+01rY8ogO2Mo4yIiIrjyyit59913+eyzz7jnnntYu3Ytl112Gc2aNWP69Ok2BIcxpkBcmRjK2u2qoVK3bl0ef/xxPv/8c2bOnElERASDBw+mZs2ajBo1ij179jgdojGmBHBlYjDBqVChAv379ycjI4PU1FS6dOnChAkTqF+/Pr169eL999+3MZqMMXlyZWKwS0mhceq5iAULFrB7927uueceUlNT6dy5My1btuSVV16xZyKMMX/gysRgl5JCr1atWjz++OPs27ePF154gSNHjtC3b1/q1avHv//9b3766SenQzTGuIQrE4MpOhUqVODmm29my5YtLF68mAYNGnDXXXdRq1Yt7r77brvd1RjjzsRgl5KKXlhYGJdffjnvv/8+GzZsoHv37jz99NPUrVuXG264ge3btzsdojHGIa5MDHYpqXjFxMQwd+5csrOzGTZsGG+++SYXX3wx11xzDVlZWU6HZ4wpZq5MDMYZderU4emnn2bPnj3cd999LFu2jOjoaJKSkvjoo4+cDs8YU0wsMZg/OO+8807/qtyjjz7K2rVradOmDd26dSM1NdXp8IwxRcwSg8lTlSpVeOCBB9i7dy//93//x6ZNm4iLi6NDhw4sX77cnoUwppRyZWKwzmd3qVSpEnfddRe7d+9m8uTJ7N69m3/84x+0adOGlJQUSxDGlDKuTAzW+exOFSpUYPjw4WRnZzN9+nS+/fZbevbsSYsWLSxBGFOKuDIxGHc766yzGDhwIJ9++imzZ8/m8OHD9OzZk7Zt27JixQqnwzPGBMkSgym08PBwrr/+erZt23b650i7dOlC586dWb9+vdPhGWMKyRKDCVp4eDgDBgzg008/ZdKkSWzdupV27dqRmJhoz0EYUwIVW2IQkXoiMkNEXi+udZriVb58eUaMGMFnn33G448/zpo1a4iOjqZ3797s2LHD6fCMMQEKKDGIyEwROSAiW/zK40Vkh4hki8jo/OpQ1V2qOiCYYE3JUKlSJe69997TPz+6ePFimjRpwoABA/jiiy+cDs8YcwaBnjHMAuJ9C0SkHPAskAA0Aa4VkSYi0kxEFvu9qoY0alMiVKlShXHjxrFr1y6GDx/OnDlzaNiwIQ8++CA///yz0+EZY/IQUGJQ1VTge7/iWCDbeyZwFJgH9FTVzaraw+91INCARGSQiKSLSPrBgwcD3hDjXlWrVmXixIls376dpKQkxo0bR4MGDXj++eft9yCMcaFg+hhqAPt83u/3luVKRM4VkalACxG5N6/5VHW6qsaoasxf//rXIMIzblOvXj3mzZtHWloajRs35tZbb6VZs2YsXLjQnoEwxkWCSQySS1men25V/U5Vh6hqfVV9It+K7cnnUi02NpZVq1bx9ttvo6pcccUVdOzYkQ0bNjgdmjGG4BLDfqCWz/uawJfBhWPKChGhZ8+ebNmyhWeffZbt27cTGxvLjTfeyFdffeV0eMaUacEkhg1AQxGpKyKRQG8gJRRB2ZAYZUdERAS33nor2dnZjBo1ildffZVGjRoxfvx4jhw54nR4xpRJgd6uOhdYDzQWkf0iMkBVjwPDgGXAdmCBqm4NRVB2KansOfvssxk/fjxbt26lY8eOjB49mqZNm7Jo0SLrfzCmmAV6V9K1qnq+qkaoak1VneEtX6Kqjbz9Bo8VbaimLGjYsCGLFi1i6dKlhIeHk5SUREJCAp988onToRlTZrhySAy7lGTi4+PZtGkTEyZMYP369TRr1ow77riDH374wenQjCn1XJkYjAFP/8PIkSPZuXMn/fv3Z+LEiTRq1IgXX3yRkydPOh2eMaWWKxOD9TEYX1WrVmX69Omkp6fTuHFjBg4cyN///ncyMzOdDs2YUsmVicEuJZnctGzZktTUVGbPns1nn31Gq1atGDlypA2vYUyIuTIx2BmDyYuIcP3117Njxw4GDhzIpEmTuPDCC3nttdfs7iVjQsSVicHOGMyZ/OUvf2Hq1KmsX7+eqlWrcs0115CQkEB2drbToRlT4rkyMRgTqDZt2rBhwwYmTZrEunXraNq0KY888og9HGdMECwxmBIvPDycESNG8Mknn9CzZ0/GjBlDs2bN+O9//+t0aMaUSK5MDNbHYAqjevXqzJ8/n2XLlqGqdO3aleuuu44DBwIe9d0Yg0sTg/UxmGB069aNzZs3M2bMGF5//XUuuugiZs+ebZ3TxgTIlYnBmGCVL1+esWPHkpmZSePGjenXrx/x8fHs2bPH6dCMcT1LDKZUa9KkCWvWrOGZZ55h3bp1XHzxxUycOJETJ044HZoxruXKxGB9DCaUwsLCGDZsGFu3biUuLo6RI0fSrl07Nm/e7HRoxriSKxOD9TGYonDBBRfwzjvv8Morr7Br1y5atmzJQw89ZLe2GuPHlYnBmKIiIlx33XVs376d3r178+ijjxIdHc3atWudDs0Y17DEYMqk8847j5dffpklS5Zw+PBhLr30UoYPH86vv/7qdGjGOM4SgynTEhIS2Lp1K8OGDWPKlCk0b96c1NRUp8MyxlHFmhhE5AoReUFEFopIt+JctzF5qVSpEpMnT2bVqlUAxMXFcdttt9nZgymzAk4MIjJTRA6IyBa/8ngR2SEi2SIyOr86VPVtVR0I3AgkFypiY4pIhw4d2LRpE8OHD2fy5MlERUWxevVqp8MyptgV5IxhFhDvWyAi5YBngQSgCXCtiDQRkWYistjvVdVn0Qe8yxnjKhUrVmTy5Mm8//77nDx5kri4OG6//XYOHz7sdGjGFJuAE4OqpgLf+xXHAtmquktVjwLzgJ6qullVe/i9DojHeGCpqmaEbjOMCa2OHTuyadMmbr31ViZNmkRUVBRr1qxxOixjikWwfQw1gH0+7/d7y/IyHOgC/FNEhuQ2g4gMEpF0EUk/ePBgkOEZU3iVKlViypQprFy5kuPHj9OhQwfuuOMOO3swpV6wiUFyKctzpDJVnayqrVR1iKpOzWOe6cDDQEZkZGSQ4RkTvE6dOrF582aGDBnC008/TXR0NOvWrXM6LGOKTLCJYT9Qy+d9TeDLIOs0xnUqVarEc889x4oVKzh69Cjt27fnrrvu4vfff3c6NGNCLtjEsAFoKCJ1RSQS6A2kBBuUDYlh3Kpz585s3ryZQYMG8e9//5vWrVuTmZnpdFjGhFRBbledC6wHGovIfhEZoKrHgWHAMmA7sEBVtwYblA2iZ9yscuXKTJ06lSVLlvDtt98SGxvLE088YSO2mlJD3PzjJTExMZqenu50GMbk6bvvvuOWW27htddeo127dsyePZv69es7HZYp40TkY1WNKezyrhwSw84YTElx7rnnMn/+fObMmcPWrVuJiopi+vTp9mtxpkRzZWKwPgZTkogIffr0YfPmzbRt25bBgwfTo0cPvvrqK6dDM6ZQXJkY7IzBlES1atVi+fLlTJo0iZUrV9KsWTNef/11p8MypsBcmRjsjMGUVGFhYYwYMYKNGzdSp04drr76aq6//np++OEHp0MzJmCuTAzGlHQXXngh69evZ8yYMcydO5fmzZuzcuVKp8MyJiCuTAx2KcmUBhEREYwdO5Z169ZRoUIFLrvsMkaOHGkPxRnXc2VisEtJpjSJjY1l48aNDB06lIkTJ9K6dWs2b97sdFjG5MmVicGY0uZPf/oTU6ZMYcmSJRw8eJDWrVszadIkTp486XRoxvyBKxODXUoypVVCQgKbNm2iW7du3H777SQkJNhtrcZ1XJkY7FKSKc2qVq3KwoULef7551m9ejXNmjVj4cKFTodlzGmuTAzGlHYiwpAhQ8jIyOCCCy7giiuuYPDgwfY708YVLDEY46ALL7yQDz/8kFGjRvHCCy/QsmVLbHww4zRLDMY4LDIykvHjx7NixQoOHz7MJZdcYqO1Gke5MjFY57Mpizp16sSmTZu48sorue++++jcuTOff/6502GZMsiVicE6n01Z9Ze//IX58+cza9YsMjIyaN68OfPmzXM6LFPGuDIxGFOWiQj9+vUjMzOTJk2acO2113LDDTfw008/OR2aKSMsMRjjUvXr1yc1NZWxY8fy6quvEhUVxdq1a50Oy5QBlhiMcbHw8HDGjBnD6tWrCQsLo0OHDowdO5bjx487HZopxYotMYjIRSIyVUReF5Fbimu9xpQGl1xyCZmZmfTt25eHH36Yjh07snfvXqfDMqVUQIlBRGaKyAER2eJXHi8iO0QkW0RG51eHqm5X1SHANUChf4vUmLKqcuXK/Oc//+GVV15h06ZNREVFsWDBAqfDMqVQoGcMs4B43wIRKQc8CyQATYBrRaSJiDQTkcV+r6reZZKANcCKkG2BMWXMddddR2ZmJhdeeCHJycncfPPN9sS0CamAEoOqpgLf+xXHAtmquktVjwLzgJ6qullVe/i9DnjrSVHVdkCfvNYlIoNEJF1E0g8ePFi4rTKmlKtXrx6rV6/m/vvvZ+bMmbRs2ZKMjAynwzKlRDB9DDWAfT7v93vLciUiHUVksohMA5bkNZ+qTgceBjIiIyODCM+Y0i0iIoJx48axYsUKfv31V9q2bcuECRNsKG8TtGASg+RSpnnNrKofqOoIVR2sqs/mV7E94GZM4Dp16kRWVhbdu3fnzjvv5PLLL+ebb75xOixTggWTGPYDtXze1wS+DC4cDxsSw5iCOffcc3nrrbd4/vnn+eCDD2jevDnvvvuu02GZEiqYxLABaCgidUUkEugNpIQmLGNMQZ0ayjs9PZ1q1aqRkJDAHXfcwZEjR5wOzZQwgd6uOhdYDzQWkf0iMkBVjwPDgGXAdmCBqm4NRVB2KcmYwrv44otJS0tj2LBhPP3007Rt25ZPPvnE6bBMCSKqeXYLOEZEEoHEBg0aDNy5c6fT4RhTYi1atIj+/fvz22+/MWnSJAYMGIBIbt2DpjQRkY9VtdDPi7lySAw7YzAmNBITE9m0aROXXHIJAwcO5JprruHQoUNOh2VczpWJwTqfjQmd6tWrs3z5cp588knefvttoqOjWbNmjdNhGRdzZWKwMwZjQissLIx77rmHtWvXEhERQVxcnA3GZ/LkysRgjCkasbGxbNy4kT59+thgfCZPrkwMdinJmKJTuXJlZs+ezZw5c2wwPpMrVyYGu5RkTNHr06cPmZmZNG7cmOTkZAYOHMjhw4edDsu4gCsTgzGmeNSrV481a9YwevRoZsyYQUxMDJs2bXI6LOMwVyYGu5RkTPGJiIjgiSeeYPny5Rw6dIjY2Fiee+453PiMkykerkwMdinJmOLXpUsXsrKy6NSpE0OHDuWqq67i++/9R9s3ZYErE4MxxhlVq1blnXfe4V//+heLFy+2Zx7KKEsMxpgcwsLCuPPOO1m3bh2RkZHExcXx6KOPcuLECadDM8XEEoMxJlcxMTFkZGTQu3dvHnroIbp06cIXX3zhdFimGLgyMVjnszHucPbZZzNnzhxeeuklPvroI6Kioli8eLHTYZki5srEYJ3PxriHiHDjjTeSkZFBrVq1SExM5LbbbrPfeSjFXJkYjDHu07hxY9avX8+IESOYPHkybdu2ZceOHU6HZYqAJQZjTMDKly/PpEmTSElJYd++fbRq1Yr//Oc/9sxDKWOJwRhTYImJiWRlZRETE8ONN97I9ddfz88//+x0WCZEijUxiEhFEflYRHoU53qNMaFXo0YNVqxYwcMPP8zcuXNp0aIF6enpTodlQiDQ33yeKSIHRGSLX3m8iOwQkWwRGR1AVfcANoyjMaVEuXLleOihh/jggw84cuQI7dq1Y8KECZw8edLp0EwQAj1jmAXE+xaISDngWSABaAJcKyJNRKSZiCz2e1UVkS7ANuCbEMZvjHGBSy+9lKysLC6//HLuvPNOevTowYEDB5wOyxRSQIlBVVMB/0FTYoFsVd2lqkeBeUBPVd2sqj38XgeATkBb4DpgoIjkum4RGSQi6SKSfvDgwUJvmDGmeJ1zzjm8+eabTJkyhZUrVxIVFcWKFSucDssUQjB9DDWAfT7v93vLcqWq96vq7cCrwAuqmuu5pqpOV9UYVY3561//GkR4xpjiJiIMHTqUtLQ0qlSpQteuXbnvvvs4duyY06GZAggmMUguZWe8Z01VZ6lqvo9O2pPPxpRsUVFRpKenc9NNN/HEE0/QoUMH9uzZ43RYJkDBJIb9QC2f9zWBL4MLxxhTWlSsWJEXX3yRefPmsW3bNqKjo3nttdecDssEIJjEsAFoKCJ1RSQS6A2khCIoGxLDmNIjOTmZjRs30rhxY6655hoGDx5sPyHqcoHerjoXWA80FpH9IjJAVY8Dw4BlwHZggapuDUVQdinJmNLl1E+Ijho1iunTpxMbG8uWLVvOvKBxRKB3JV2rqueraoSq1lTVGd7yJaraSFXrq+pjRRuqMaYki4iIYPz48SxbtoyDBw/SunVrpk2bZsNpuJArh8SwS0nGlF7dunVj06ZNdOjQgSFDhnD11Vdz6NAhp8MyPlyZGIwxpVu1atVYunQpTz31FAsXLiQ6Opq1a9c6HZbxcmVisD4GY0q/sLAw7r77btauXUt4eDhxcXE89thj9hOiLuDKxGCXkowpO2JjY8nIyODqq6/mgQceoGvXrnz5pd357iRXJgY7YzCmbPnzn//Mq6++yowZM0hLSyMqKoolS5Y4HVaZ5crEYGcMxpQ9IsJNN91Eeno61atXPz0g39GjR50OrcxxZWIwxpRdF110EWlpaQwdOpQJEybw97//nc8++8zpsMoUSwzGGNcpX748U6ZM4c033yQ7O5sWLVowb948p8MqM1yZGKyPwRgDcOWVV5KZmUmzZs249tprufnmm/n111+dDqvUc2VisD4GY8wptWvXZtWqVdx3333MnDmT1q1bs3nzZqfDKtVcmRiMMcZXeHg4jz32GMuXL+fQoUPExsYydepUG06jiFhiMMaUGF26dCEzM5O4uDhuueUWrr76an744Qenwyp1XJkYrI/BGJOXatWqsWTJkhzDaXz44YdOh1WquDIxWB+DMSY/p4bTWLNmDSJC+/btGT9+PCdP5vqLwaaAXJkYjDEmEG3atGHjxo306tWL0aNHEx8fzzfffON0WCWeJQZjTIlWpUoV5s+fz7Rp01i9ejVRUVG89957TodVolliMMaUeCLCoEGD2LBhA+eeey7dunXj3nvv5dixY06HViIVW2IQkY4islpEpopIx+JarzGm7GjatCkbNmxg4MCBPPnkk3To0IE9e/Y4HVaJE+hvPs8UkQMissWvPF5EdohItoiMPkM1CvwClAf2Fy5cY4zJ35/+9CemT5/OvHnz2LZtG9HR0bzxxhtOh1WiBHrGMAuI9y0QkXLAs0AC0AS4VkSaiEgzEVns96oKrFbVBOAe4OHQbYIxxvxRcnIyGzdupFGjRvzzn//klltu4bfffnM6rBIhoMSgqqnA937FsUC2qu5S1aPAPKCnqm5W1R5+rwOqeuo+skPAWSHbAmOMyUO9evVYs2YNd999N1OnTqVNmzZs377d6bBcL5g+hhrAPp/3+71luRKRXiIyDXgZmJLPfINEJF1E0g8ePBhEeMYYA5GRkTz11FMsXbqUr7/+mlatWjFz5kwbTiMfwSQGyaUszz2tqm+q6mBVTVbVD/KZbzqeS00ZkZGRQYRnjDH/Ex8fT1ZWFpdccgkDBgygT58+/PTTT06H5UrBJIb9QC2f9zUB+6FWY4xrnX/++Sxfvpxx48axYMECWrZsSXp6utNhuU4wiWED0FBE6opIJNAbSAlFUDYkhjGmqJQrV47777+fDz74gKNHj9KuXTsmTJhgw2n4CPR21bnAeqCxiOwXkQGqehwYBiwDtgMLVHVrKIKyQfSMMUWtffv2ZGZmnv5t6cTERKxf00Pc3AETExOjdppnjClKqspzzz3HHXfcwXnnnccrr7xCx44dnQ4rKCLysarGFHZ5Vw6JYWcMxpjiIiIMHTqUtLQ0KleuTOfOnRkzZgzHjx93OjTHuDIxWB+DMaa4RUdHk56eTr9+/XjkkUfo3Lkz+/btO/OCpZArE4OdMRhjnFCpUiVeeuklXn75ZTZu3Eh0dDQpKSG5p6ZEcWVisDMGY4yT+vbtS0ZGBnXq1KFnz57cdtttHDlyxOmwio0rE4MxxjitYcOGrFu3jttuu43JkydzySWX8OmnnzodVrFwZWKwS0nGGDc466yzmDhxIgsXLmTv3r20bNmSl19+2emwipwrE4NdSjLGuElSUhJZWVm0bNmSG264gX79+vHLL784HVaRcWViMMYYt6lZsyYrV65kzJgxzJkzh1atWpGZmel0WEXClYnBLiUZY9woPDycsWPHsmLFCn755RfatGnDlClTSt1Ira5MDHYpyRjjZh07diQrK4uuXbsyfPhwrrzySr7/3v8na0ouVyYGY4xxu/POO49FixYxYcIElixZQnR0NGvWrHE6rJCwxGCMMYUkIowcOZJ169YRGRlJXFwc48aN48SJE06HFhRLDMYYE6SYmBgyMjLo3bs3Dz74IF27duXLL0vuz9O4MjFY57MxpqQ5++yzmTNnDjNnziQtLY3o6GiWLl3qdFiF4srEYJ3PxpiSSETo378/6enp/O1vf6N79+7cfffdHD161OnQCsSVicEYY0qyiy66iLS0NG699Vb+9a9/0b59e3bt2uV0WAGzxGCMMUWgQoUKPPvss7zxxhvs3LmTFi1aMH/+fKfDCoglBmOMKUK9evUiMzOTiy++mN69ezNw4EAOHz7sdFj5KrbEICJhIvKYiDwjIv2Ka73GGOO02rVrs2rVKu69915mzJhB69at2bJli9Nh5SmgxCAiM0XkgIhs8SuPF5EdIpItIqPPUE1PoAZwDNhfuHCNMaZkioiI4PHHH2fZsmV89913tG7dmmnTprlyOI1AzxhmAfG+BSJSDngWSACaANeKSBMRaSYii/1eVYHGwHpVvQO4JXSbYIwxJUfXrl3JysqiQ4cODBkyhOTkZH744Qenw8ohoMSgqqmA/0AgsUC2qu5S1aPAPKCnqm5W1R5+rwN4zhIOeZfN87FAERkkIukikn7w4MGCb5ExxrhctWrVWLp0KePHj+ett96iRYsWfPjhh06HdVowfQw1AN9fyt7vLcvLm8A/ROQZIDWvmVR1OvAwkBEZGRlEeMYY415hYWGMGjWK1atXA3DppZfy1FNPcfLkSYcjCy4xSC5leV4sU9XDqjpAVYer6rP5VWwPuBljyoq2bduyceNGrrjiCu655x4SEhL45ptvHI0pmMSwH6jl874mEJLBQWxIDGNMWVKlShUWLFjA1KlTSU1NJSoqiv/+97+OxRNMYtgANBSRuiISCfQGUkITljHGlC0iwuDBg/noo48455xz6NatG/fddx/Hjh0r9lgCvV11LrAeaCwi+0VkgKoeB4YBy4DtwAJV3RqKoOxSkjGmrGrWrBnp6ekMGDCAJ554gri4OPbu3VusMYgb76EVkUQgsUGDBgN37tzpdDjGGOOIefPmMWjQIMqVK8eMGTPo1atXQMuJyMeqGlPY9bpySAw7YzDGGOjduzeZmZk0bNiQq666iqFDh/L7778X+XpdmRis89kYYzzq1avHmjVruOuuu3juuedo06YN27dvL9J1ujIx2BmDMcb8T2RkJP/3f//HkiVL+PLLL4mJieGll14qsuE0XJkYjDHG/FFCQgJZWVm0adOGm266ib59+/LTTz+FfD2uTAx2KckYY3JXvXp13nvvPR599FHmzZtHy5YtSU9PD+k6XJkY7FKSMcbkrVy5cjzwwAOsWrWKI0eO0K5dOyZOnBiyS0uuTAzGGGPOrH379mRlZdG9e3dGjhxJUlIS3377bdD1ujIx2KUkY4wJzDnnnMNbb73FM888w/Lly4mKigq6TlcmBruUZIwxgRMRhg0bRlpaGpUqVQq6PlcmBmOMMQUXHR1NVlZW0PVYYjDGmFKkfPnyQddhicEYY0wOrkwM1vlsjDHOcWVisM5nY4xxjisTgzHGGOdYYjDGGJODJQZjjDE5WGIwxhiTgyt/2vMUEfkZ2OF0HMCfgVDdIlXYugqyXCDz5jdPQaflNf95QPADtwQnlG0XTH2hbL/CTi9IuRvaDuyzF8i03Mobq2rlM8SRN1V17QtIdzoGbxzTna6rIMsFMm9+8xR0Wl7zu6H9Qtl2bmm/wk4vSLkb2i7U7eeGtjvTPIWZVhTtZ5eSArPIBXUVZLlA5s1vnoJOC+X+CbVQx+aG9ivs9IKWu4F99s48LeTt5/ZLSemqGuN0HKZwrP1KLmu7ki3Y9nP7GcN0pwMwQbH2K7ms7Uq2oNrP1WcMxhhjip/bzxiMMcYUM0sMxhhjcrDEYIwxJocSmxhEpJ6IzBCR152OxZyZiFQUkf+IyAsi0sfpeEzB2Oet5BKRK7yfu4Ui0i2QZRxJDCIyU0QOiMgWv/J4EdkhItkiMjq/OlR1l6oOKNpITX4K2I69gNdVdSCQVOzBmj8oSPvZ581dCth2b3s/dzcCyYHU79QZwywg3rdARMoBzwIJQBPgWhFpIiLNRGSx36tq8YdscjGLANsRqAns8852ohhjNHmbReDtZ9xlFgVvuwe8088oPDQxFoyqpopIHb/iWCBbVXcBiMg8oKeqPgH0KOYQTQAK0o7AfjzJIZMSfAmzNClg+20r5vBMPgrSdiKyHXgSWKqqGYHU76YPaA3+9x8leL5IauQ1s4icKyJTgRYicm9RB2cCllc7vglcJSLP4+4hGMq6XNvPPm8lQl6fveFAF+CfIjIkkIocOWPIg+RSlufTd6r6HRDQRppilWs7quqvQP/iDsYUWF7tZ58398ur7SYDkwtSkZvOGPYDtXze1wS+dCgWU3jWjiWbtV/JFbK2c1Ni2AA0FJG6IhIJ9AZSHI7JFJy1Y8lm7VdyhaztnLpddS6wHmgsIvtFZICqHgeGAcuA7cACVd3qRHwmMNaOJZu1X8lV1G1ng+gZY4zJwU2XkowxxriAJQZjjDE5WGIwxhiTgyUGY4wxOVhiMMYYk4MlBmOMMTlYYjDGGJODJQZjjDE5WGIwxhiTw/8DrP6W0pQHu1EAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "s = np.linspace(0.1, 100, 1000)\n", "plt.plot(s, 1 / (np.sqrt(np.pi) * s) - erfcx(s), c=\"black\")\n", "plt.xlim(s[0], s[-1])\n", "plt.xscale(\"log\")\n", "plt.yscale(\"log\")\n", "plt.title(r\"Absolute error of first order asymptotics of $\\mathrm{erfcx}(s)$\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "-----------------------------\n", "### License\n", "\n", "This file is part of NEST. Copyright (C) 2004 The NEST Initiative\n", "\n", "NEST is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 2 of the License, or (at your option) any later version.\n", "\n", "NEST is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.5" } }, "nbformat": 4, "nbformat_minor": 1 }