Skip to content

Instantly share code, notes, and snippets.

@philippkraft
Last active March 6, 2019 10:33
Show Gist options
  • Save philippkraft/aae02d23fbdad62f98a413ab04fe6d83 to your computer and use it in GitHub Desktop.
Save philippkraft/aae02d23fbdad62f98a413ab04fe6d83 to your computer and use it in GitHub Desktop.
Discussion of an overflow smoothing function as introduced by Clark et al 2008
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# The $\\Phi$ Function\n",
"\n",
"The logistic function to smooth threshold or overflow processes given by [Clark et al 2008](https://doi.org/10.1029/2007WR006735) eq. 12h reads as\n",
"\n",
"$\n",
"\\Phi(S, S_{max}, \\rho_{S}, \\varepsilon) = \\frac{1}{1 + \\exp\\left(\\frac{S - S_{max} - \\omega \\varepsilon}{\\omega}\\right)}\n",
"$\n",
"\n",
"Where:\n",
"\n",
"- $\\omega = \\rho_{S} \\cdot S_{max}$\n",
"- $\\varepsilon = 5$ a coefficient \"that ensures that $S$ does not exceed $S_{max}$\" (Clark et al 2008).\n",
"\n",
"The following will show that this is not the case.\n",
"\n",
"We assume the following single storage system:\n",
"\n",
"$\n",
"\\frac{dS}{dt} = Q_{in} - Q_{out}(S)\n",
"$\n",
"\n",
"and understand $Q_{out}(S)$ as an overflow process:\n",
"\n",
"$\n",
"Q_{out}(S) = Q_{in} \\cdot \\left(1 - \\Phi\\left(S, S_{max}, \\rho_{S}, \\varepsilon\\right)\\right)\n",
"$\n",
"\n",
"## Implementation\n",
"\n",
"### The system equations"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"from math import exp\n",
"\n",
"def Phi(S, S_max, rho_s, eps):\n",
" omega = rho_s * S_max\n",
" return 1/(1 + exp((S - S_max - omega * eps)/omega))\n",
"\n",
"def Qout(Qin, S, S_max, rho_s, eps):\n",
" return Qi * (1 - Phi(S, S_max, rho_s, eps))\n",
"\n",
"def dSdt(Qin, S, S_max, rho_s, eps):\n",
" return Qin - Qout(Qi, S, S_max, rho_s, eps)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Solution\n",
"\n",
"We solve it with a simple explicit Euler method, which is ok, since this differential equation has no stiffness:\n",
"\n",
"$\n",
"S(t+1) = S(t) + \\Delta t \\frac{dS}{dt}\\left(t, S\\right)\n",
"$\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def euler(timesteps, S0, Qi, S_max, rho_s, eps):\n",
" S = [S0]\n",
" dt = timesteps[1] - timesteps[0]\n",
" for t in timesteps[1:]:\n",
" old_S = S[-1]\n",
" new_S = old_S + dt * dSdt(Qi, old_S, S_max, rho_s, eps)\n",
" S.append(new_S)\n",
" return S\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Results\n",
"\n",
"We use the following configuration:\n",
"\n",
"$\n",
"S(t_0) = 90 mm, S_{max} = 100 mm \\\\\n",
"Q_{in} = 0.2 \\frac{mm}{min} = 288 \\frac{mm}{day} \\\\\n",
"\\varepsilon = -5, 0, 1, 5, 10, \\rho_S = 0.01 \n",
"$\n",
"\n",
"And plot the results with:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAmwAAAF6CAYAAACk8gQ+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd4HNd97//32d3Z2Y5OgL1TJCVRlERKkWzJVLMtSy6xLd/opzi2r2/cS3KfOHYS3zi/5Dr2ja/zJD/Lciwn7orkFrk32RZkyValClVIikXsBAhgsWjbZmfO74+ZXewCCxIgAewC+L6eZ54zc87ZnYNh+/BMU1prhBBCCCFE/fLVegBCCCGEEOL0JLAJIYQQQtQ5CWxCCCGEEHVOApsQQgghRJ2TwCaEEEIIUecksAkhhBBC1DkJbEIIIYQQdU4CmxBCCCFEnZPAJoQQQghR5wK1HsB0am1t1atWrZrx/YyMjBCNRmd8P3ONHJeJybGpTo7LxOTYVCfHZWJybKqr5+Oyc+fOXq1122T6zqvAtmrVKp544okZ309nZyc7duyY8f3MNXJcJibHpjo5LhOTY1OdHJeJybGprp6Pi1Lq8GT7yilRIYQQQog6J4FNCCGEEKLOSWATQgghhKhzEtiEEEIIIeqcBDYhhBBCiDongU0IIYQQos5JYBNCCCGEqHMS2IQQQggh6pwENiGEEEKIOieBTQghhBCizklgE0IIIYSoc/PqXaJCCCGEEGfDcTR2wcEpONgFjfJBOBas9bBKJLAJIYQQYtYUg5FtOZVloTwwjW6X6mxd6uvYo23VPje67tDX69D76BOjnyvuz/b6We66dnTFOFde0MLNH7ioRkdpPAlsQgghxDyntcYpC0CFsWGpPDRZY/pM2M/GLuhSP+e0n9EUCg6O5eCMCUbnQvkUfr/Cb/jwBXzuesBbD7jrAMGQH3/A8OrdttF1X+k7/AEfPu87Eq3haRvndJDAJoQQQswwrUdnlQqWV+a9YJS3KRQc7HyxzaZQ3q9U2hMGK6fg0N/vcPK3j04YwKaDUpSCTbEMFMOStx4M+fEbwVJgqta/FJSM0WBVXHxn2C4PWz6fOuOYOzs72bHj4mn5+WtJApsQQogFQ3un48pDUEUwyjsUvBDlhqnR8FNeVy1gFb+jWp1tnVtg8vmVG3TGhZ3RwBQwoXFRxA1AY0JS9cCk8Bv+yv5nCFY+v9yrWCuzFtiUUl8GbgZOaa0v8OpuAf4O2ARcprV+wqtfBewG9noff0Rr/Z7ZGqsQQojZ5dhuELLy7iySlbe9sGRj5b1ZqLxD8oDmGfuoG7TK6gt5b1Yq7/Yf+7li/3MNTsXwUgxPlaWfSDgwrq60HfSCUNDvlb5SOa6v4SNg+N2gZExlJunCc/r5RP2azRm2rwK3A18vq3sOeCPwxSr9D2itt87CuIQQQpyGbTsUcjZWzgtEOdsNVF5p5YrrY4KSNSZQTVRvuReDT9bJx/e5Kwo32AT9XuhxSyPoxzD9hGJBjKDX7vXzByvDUHmICgR8pfaxoak426QmEZyEmAmzFti01r/1Zs7K63YDKCV/AIQQ4lycNlRVrDsVYau83cq5IcotR+udwtQuElc+hRH04Q/6xwWmUMyo2A6MaTeC3kyTF7xK7V4ge3zno1z1ipeXgpb8+yEWCqX19N2tccaduYHtx8VTomX1ncBfjDkl+jzwIjAIfFxr/eAE3/ku4F0A7e3tl95zzz0zM/gyw8PDxGKxGd/PXCPHZWJybKpbqMfFvWMPd7EYt25bkE3nCPhMHKuyrz2mr7a9copn+pQffH7wBdxFBSq3yxflV1Xrfd5n1NjtGZyFWqi/ZyZDjk119Xxcrrnmmp1a622T6VuvNx2cBFZorfuUUpcC31dKna+1HhzbUWt9J3AnwLZt2/SOHTtmfHDudQIzv5+5Ro7LxOTYVDeXjotdcMhnC+QzBfIZ2y1zNla2QD5rY2VtrFyxzq23crbb5s1i5bMFty3vwBn/rxwENAHDhxHyY4QCBE0/wYQfwwxgmH63PujHML3ZKdNfKg3TbQuYle3FuslcE1WP5tLvmdkmx6a6+XJc6jKwaa1zQM5b36mUOgBsAJ6o6cCEEHOOdjRWziaXKXiByy4LXlW2s17fUn+3brIXqwdMP8FimDL9BEMBIg1Bd910g1d5WzFcBb1QVlx/9PGHueb6HXM2WAkhplddBjalVBuQ1FrbSqk1wHrgYI2HJYSoAe1ocpkCubQbnrJpi9xIgVzaKtXn0gXy3va4AJazzzybpSAYChAM+zHDATdkJYI0LgoTDAfcJeSVYf/oennw8ma7put0oD+oJKwJIUpm87EedwM7gFal1DHgE0AS+BzQBvxEKfW01vpVwNXA3yulCoANvEdrnZytsQohppdtO6MhK10MWRbJfZon0ofGhK/yPu4s1+kCl8+nCEYCmJGAG7bCASINEYIhf0XYMsNuqDKrBDDD9MvF60KIujabd4neOkHTvVX6fg/43syOSAgxVVprCnmHzHCe3EiB7LBFdqRsGa6yPmyRz9oTfufJnQcJGD43cEUNzHCAWKNJy5JYKYiFIgZmxA1YoWgAM2IQDLttEraEEAtBXZ4SFULMDsd2yAxbZIYsMkN5b/HWRyxy4wJY4bSvuAmG3Mc2hKLu0tgecddjBmbEKIUt0wtiO59+jGuuewV+Q56eLoQQpyOBTYh5RGtNLl2oDF5DedLVAtmQG8SqUT5FKBooha1Ea5hFKxOl7YqyGMiiAfxTfG1NIKQkrAkhxCRIYBNiDrALDunBPCMDOdIDedIDOUaK5WC+VJcZsnCc6hd8haIG4bhBOB6keUmUcDxIOB4k4tWFS2UQMxyQJ7oLIUQdkcAmRA05jiY9kGO4310qAtngaDDLDleZCVMQjgeJNgSJJExal8WIJILjwlc47s6ATXX2SwghRP2QwCbEDNGOJj2YZ7g/x+BRzTO/Pspwf9YLZ1kvoOXRY2bEfH5FJBEk0mCSaA3TsbbRC2VBog0mkQa3DMcNfBLChBBiQZDAJsRZcmzHDWN9WQZ7Mwz1ZRnsc8vhZI6RVK7i9ORR9uE3fMSaTGJNIZae11Rad0uTaINJKGrI6UghhBAVJLAJMQGtNZkhi4FTaQb7sgz1ZRjsdUPZYK87Q1Y+O6YUxJpCxFtCLF7f4K43mUSbQuw58CzX3PByzGhAHkEhhBBiyiSwiQUvO2IxcCpD6lSa1Kk0A91pUt62Neb5YZFEkERriI41DSRaQyRawsS9MtZsTnid2OGkIhQzZuPHEUIIMQ9JYBMLgnY0g31ZkidHSJ4YJtWV9gJapuKCfqUg3hKicVGEjjWLaWwP09AWIdEaIt4cIhD01/CnEEIIsVBJYBPzinY0Q8ksyRMjXjhzy/6TIxTKXt4dbQjS2BFhzcVtNC6K0LgoTGN7hERLWJ4LJoQQou5IYBNzll1wSJ4YoefoEL1Hhug5Okzv8WEKudHTmNEG95lj51+1lOYlUZqXRGlaHMUMy299IYQQc4f8qyXmBNty6Dk6xKnDg24wOzpE8sQIju1e9G+YflqXx9h05WJalkRpXuyGMzMi140JIYSY+ySwibqjtWa4P0fXwQG6Dw7S9dIAPUeHcApuOAvHDVqXx9l6fQuty2O0LY/T0BaWR2EIIYSYtySwiZrTWpM8McLxF1OceLGfkwcHSA/kAQgYPtpWxrnomuW0r0nQvqqBaGNQHo0hhBBiQZHAJmad1pr+k2mOv9jP8Rf7ObEvRWbIvVMz1myy7LwmOtY00L46QcuymLxSSQghxIIngU3Miny2wLE9/Rx+vo8jz/Ux3J8DINZksuL8FpZuaGTphiYSreEaj1QIIYSoPxLYxIwZSmY5+FQPh57t5cS+FI6tMUJ+lm9sZttrmlm2sZlEa0hObwohhBBnIIFNTKuBngwHnjrFgSd7OHVoEICmxVEuunY5Ky5oYfHaBvwBOcUphBBCTIUENnHOsiMW+x7v5sAvHZ6/52EAFq2Mc8Ufri09mFYIIYQQZ08CmzgrjqM58nwfex7u4qVdPTgFTagRrnzjOtZe0ibXogkhhBDTSAKbmJJc2uKF353k2c5jDPVlCcUMLrh6KRuvWMzzB3Zy8Y4VtR6iEEIIMe9IYBOTMtib4an7jrDnkS4KOZsl6xu58o3rWH1R6+g1aQdqO0YhhBBivpLAJk5roCfDzp8dYu8jXeCDDdvb2XLNctpWxGs9NCGEEGLBkMAmqhoZyPHYDw+y++EufD7FBa9YyiWvWkm00az10IQQQogFRwKbqFCwbJ759VF2/uwwdsHhQglqQgghRM1JYBMlx/f285tv7mGwJ8OqLa287E3raGyXR3IIIYQQtSaBTZDPFnj43gM898BxEm1hXvehrSzf3FzrYQkhhBDCI4Ftges5MsTPv/Qcg70ZLrpuOZe/fg1G0F/rYQkhhBCijAS2BUprze7fn+S3d79IKGbwh//zEpasb6z1sIQQQghRhQS2BchxNL/7zj523X+MZRubuOG/n08kEaz1sIQQQggxAQlsC0zBsvnVl1/gwFM9bLl2GS9783p8PlXrYQkhhBDiNCSwLSC25fCzf3uWI88nedmb17H1enmNlBBCCDEXSGBbIOyCw8+/9BxHnk9yzR9vZPPLl9R6SEIIIYSYJF+tByBmntaa+7+5h0O7ern6jzZIWBNCCCHmGAlsC8BTvzzC3ke62H7zai7csazWwxFCCCHEFElgm+eOvNDHw98/wLpti9h+06paD0cIIYQQZ0EC2zyWGc7z66/upqkjyrV/sgml5G5QIYQQYi6SwDZPaa3p/OZesmmLV75zs7y9QAghhJjDJLDNUwef7uHg0z1c/ro1tC6L13o4QgghhDgHEtjmIStn89C399GyNMbW65bXejhCCCGEOEcS2Oahp+47wnB/jqtv3YDPL7/EQgghxFwn/5rPM9kRi2d+dYQ1F7exZJ28zF0IIYSYDySwzTNP33eEfM7msptX13ooQgghhJgmEtjmkXy2wK7OY6y7ZBEtS2O1Ho4QQgghpom8S3Qe2fNwF1bW5qLr5UYDIYQQAgCt0ZaFzufRloWTz4NX6rxV0VZeBlqaiWzfXuvRl0hgmye0o3m28xiLViXoWN1Q6+EIIYRYYHR5MMrlSuVoMJogII0NS9b4uvEha5Lflc+zyLLYo/WUf57oVVexQgKbmG7HXuwn1Z3m+ndsrvVQhBBCzCJdKKDzedTwMFb3KXQ+NxpeioEplx8NQuV1xX7eZ5zyulwObeVxcjk3HJUFMSc/vk7n89P6c6lgEGUYlWWVOl80gjKCqKAxWlfWfvj4CVZvWO/2Gft9477XKH2XP5GY1p/nXElgmydefKybYMjP2ovbaj0UIYRYUHSh4IagXM4ts1k3wGSzpTq33itzWTcs5bz2Yn0+VxmiSuHo9HXYNgCLgP3n8HOoYBBlmqUA4ysGmfK6eMyr9+rM8r5mWX8vOJV9VhnFMFQtMLlhqxi0MIxpe53iC52dtO7YMS3fVUsS2OaBgmVz8MlTrLm4jYC8gkoIsUBpyxoNSNmsN1tUtl4ekPJesKpaPyZYZbNumCrVl+0jn4dC4ewHrRQqFHKDimlWDUz+RKIUjkohqkpg2n/kCBs2nz9BYDJLwahUNzZMyfum69qsBTal1JeBm4FTWusLvLpbgL8DNgGXaa2fKOv/V8A7ARv4kNb6F7M11rnmyHNJ8lmb9dvbaz0UIYSooLV2Z4IyGZxsFieTcUNPJuMGnmzWLTNZNzhlsjjZDDqTdQNTpqxPNktTVxcv3f75sr6jbecUnHy+0eAUCrnhyAyhTBOfaeKPxVGtbfjMIMqsbFcht0+pPhRCBU18IbMUwnym6X5vWb3PNKd1JunZzk6a5sFMkqhuNmfYvgrcDny9rO454I3AF8s7KqU2A38EnA8sAX6llNqgtbZnZ6hzy74nugnHDZad11TroQgh5pBimHLSaXQ6jZPJuEs6g5NJuyErncFJp90Qlc2NhqmyoKSzGZwJ27JwFhd8Y3gzQeEQvlDYDUGhEAD+5iY3LBXbwiE3LIVMfKFwRWgqD1M+MzguNJWCUyAgM0yirs1aYNNa/1YptWpM3W6g2h+S1wP3aK1zwEtKqf3AZcDDMz/SucW2HY4838faSxfJa6iEmIe01qMzUhOGqrS37dWV1ovb5e2VfXCcKY1HmaYbhsJhL1B5YSocwkgkSiGqGJ584RAqFHYDUiiELxx2v6P4uWJbeDSU+UzTvY6pis7OTi6SWSSxANXrNWxLgUfKto95deMopd4FvAugvb2dzs7OGR/c8PDwrOxnMkZOafJZzTBddHZ213Qs9XRc6o0cm+rm5XGxbVQuh8pmvTKHL5ct286isjlULlvWVtbf227JZnnesty6KcxQaaXQwSDaNMeU3noijjZDlXVj+5smOmi67aUyCIYBvmn6j6Ftw8iIu0zBvPw9M03k2FQ3X45LvQa2avPSVf/G0lrfCdwJsG3bNr1jFv7n1dnZyWzsZzIevvcAR3xHeNWbr8IM1/aXs56OS72RY1NdPRyX0gzW0BD28DCOt9jDwzgjI+7s1ciIt6Rx0l5ZrKtoH0HncpPbsVL4IhF80eho2dxcWu8aGGDZ2jXuzFM44s5IRcKj25EIvkjYnbEKl62HQvP61F49/J6pV3Jsqpsvx6VeA9sxoPxx/cuAEzUaS107/HwfHWsbah7WhKgFXSjgjIxgDw3jDA+5QWtoCGd4BGfYC2Be2+j68GgwGxrCHhkByzrzznw+N1SVB6xoFKOx0VuP4It45Zh+/mgUFYngL9ZHo26wOs1s1b7OTtrnwT8yQojpUa//yv8Q+E+l1D/j3nSwHnistkOqPyMDOfqODXPFH66t9VCEODu2TaG/H2dgAHtoCHtgEGdwAHtwCHtw7Pog9ogXurygpdPpM+8jEMAfi+GLxfDF4/ijUYyODnzxmFcfxxeL4Y8X16P44/HR0OUFr/k+cyWEqG+z+ViPu4EdQKtS6hjwCSAJfA5oA36ilHpaa/0qrfXzSqlvAy8ABeD9cofoeCf3DwCwdIPcHSpqR+fzFFIpN1ANDmIPDLiBamAQe3AAxwtc9uDgaDDz1tvTafad5ruVYeBraMCfSOCPx/HHExiLl7jhKhqrDF3F9XgcX9QLYPG4eyegBC0hxBw3m3eJ3jpB070T9P8k8MmZG9Hcd3J/ioDho3VFrNZDEfOA1to9vZhKYfen3LJ86e8fX5dK4ZxhlssXiVSELmPZMkKJBP5EnCPJJGu3XIQ/EceXSOBvaHBntxIN+BsS+LzHOAghxEJXr6dExSScPDBA+5oEfnmch6hC5/MU+lPYyT4KfcmKcjSEpSik+rFTA9gDA6e9lsvX0IC/sQF/YyP+tlbM9evc9eLS0OAGrUQcfyLh9o/FJnw8A8Duzk6a5TotIYQ4Iwlsc1Q+W6D36BCX3riq1kMRs0Q7DvbAAHYySaGvb7TsS1JIFsukW59M4gwMVP8iw8Df2ECgsRF/QyPm6tX4G5sqw1dTY+V2IoEKyF8XQghRK/I38BzVfXAQrWHJusZaD0WcI53PU+jtpdDTU1qsU6cqtgs9PdjJ/tJLnisohb+pCX9zE4HmFsyN5xFtbsHf0kygucWtb2nB39xMoKXFva5LrukSQog5RQLbHHXyQAqloH1NotZDERPQtu2Gra4uzJ07SR454m6fGhPEUqnxH/b53MDV1kagrY3w+efjb2lxA1hLc0UA8zc2ovz+2f8BhRBCzBoJbHNUz5EhmhZHCYbkl7AWtG1T6O2j0HUSq6u7VFpdJymc7MLq7qZw6lRpRqwR6AYwDAJtrQTa2jBWriC87dJSKKtYWlokhAkhhCiRf+3nqJ4jQyzb2FzrYcxbulBwA9ixY1jHj2EdP07+2DGs4yewTp6gcKoHCoWKz6hQCKOjg0BHB9HLLyewuAOjvQNjcQdPHz3KH9x0kzsbJqcjhRBCTJEEtjloZCDHyECethXxWg9lztJaYyeT5A8fcQPZsWNuIDt2HOv4cayTJyuvF/P5CHS0E1y6jOj27QQ6FmN0tBPo6MBYvJhAe/tpw1ihs5NAkzwvTwghxNmRwDYH9R4dBqBNnr92RvbgIPnDh8kfOkT+kFd6287wcEVff1srwaXLCG/dSuLmmzCWLiW4bBnGsmUYHR2nfTyFEEIIMZMksM1BPUeGAGhdJjNs4D7uwjpxktz+feT37yd34GApmNnJ5GhHpTCWLCG4ahUNr3sdwVWrCK5cgbF8OcaSJfKQViGEEHVLAtsc1HN0iIZFYYIL7IXvWmsKJ06Q27/fXfbtJ3fgALkDByreKRloayO4ejXx665zQ9mqlQRXrcJYvhxfMFjDn0AIIYQ4OwvrX/x5oufw0Lx/nIe2LHIHD5J9YTfZ3S+Qe2E32T17Kk5j+ttaMdeto/FNb8Jctw5z/TrMtWvxNzTUcORCCCHE9JPANsdkhy2GklkueMXSWg9l2mjLIrv3RbLP7iL7wgtkX9hNbt8+dD4PgAqHCW3YQOK1NxM67zw3nK1zX4skhBBCLAQS2OaYnmPu9Wtty+fu9WtW9ykyzzxN5plnyDzzDNnnnkdnswD4GxowN2+i6Y//mNCmTYQ2byK4apU8k0wIIcSCJoFtjkkeHwGgZdncuENUa03+0CHSjz1O+rHHSD/1JIUTJwFQhkFo82aa/ttbCG/dSnjLFgJLlshzyoQQQogxJLDNMckTw4RiBuF4/T5iIn/kCCOPPkr60cdIP/aY+8R/3JsBwtsuJfK2txG+6CLMzZvlJgAhhBBiEiSwzTHJkyM0L47W1SyUk82Sfvxx4t/6Nvs/9Smsw0cA8Le2Er3sMiKXXUbk8svcU5t1NG4hhBBirpDANodorUmeGGHD5R21HgqF3l6GfvUrhu/vZOTRR9HZLGHDIHjlFTS/9U+IXnkFwdWrJaAJIYQQ00AC2xwy3J8jn7VpWRKtyf6t7lMM3XcfQ7/4BemdO8FxMFasoPHNbyb2iqt5IpNh8ytfWZOxCSGEEPOZBLY5JHnSveGgeRYDmzMywuAv72Pg3ntJP/44aI25fh2t730v8Ve9EnP9+tFZtM7OWRuXEEIIsZBIYJtDineINi+e2TtEtdZknniC1L3fZ/DnP0en0xgrV9D6gfeTePWrMdeundH9CyGEEKKSBLY5JHlymEgiSCg2M3eIOtksAz/6Ef1f/wa5ffvwRSIkXnMjjW98I+GLL5br0YQQQogakcA2hyRPjMzI6dBCMknya18n9a1vYadSmBs3sviT/5vEjTfii0SmfX9CCCGEmBoJbHOEdjTJkyNsfvmSafvOQl8fff/xZfrvvhudzRK//jqa3vpWItu3y2yaEEIIUUcksM0RQ8kshbxD8+Jzn2Gzh4fp++KdJL/xDXQ+T+Lmm2h9z3sx16yehpEKIYQQYrpJYJsjkieKd4ie/Q0H2rYZuPdeTv3Lv2L39pJ47Wtpfd97MVdLUBNCCCHqmQS2OaL0SI/FZ3dNWfbFFzn5Nx8n++yzhLdupf2OzxPesmU6hyiEEEKIGSKBbY5IdacJJ4KYkandIaoti75//3d67vgC/liMJZ/5JxI33yzXqAkhhBBziAS2OSLVnaapfWqza9aJExz78z8n+8wuEq+5kfaPf5xAc/MMjVAIIYQQM0UC2xzR351mzda2SfcffvAhTnzkI2jLYuk/f5bEa14zg6MTQgghxEzy1XoA4syyIxbZYYvGSc6wJb/xTY6+610E2tpY9d3vSFgTQggh5jiZYZsDUt1pgDOeEtWOw6nPfpbkf3yZ2LXXsvT/fkYefCuEEELMAxLY5oD+LjewnW6GTWtN9yf/kf677qLx1j+i4+MfR/n9szVEIYQQQswgCWxzQKo7jc+vSLSGqrZrrTn1T5+h/667aH7721n00b+Uu0CFEEKISdBaU9AFLNvCctwlb+cxfAZtkclfOz7TJLDNAanuNA1tYXz+6pccJv/jP0h+5Ss03XabhDUhhBB1Z2woytt58k4ey7bc0rEq2qqVU20rfl9vqpfbf3h75f7G7Eujx435qqVXccf1d9TgaFUngW0O6O9OT3g6dOg3v+HUZ//ZfWzH3/y1hDUhhBAA2I5N3sm74cgLSDk754YW210f215cz9m5UqDJ2bmqfcZ+XzF45e2yQFQWjqqFonNh+AyC/iBBXxDDZ2D4jVJdeRn1RemIdVRtK33GFxzXtji2eFrHe64ksNU5x9EM9KRZdWHLuLb8oUOc+IuPEDr/fBb/4z+ifHLTrxBC1JrWuhRmcoWcW3pLtQBU3LYcq7KtSogqBqBxYcvOM5QZQt2jSiGqoAvT8vMUQ4zpN8etm36ToD9I2AhXDT3FUFQMVUF/8LRtpbIsSFXrE/AFJj1B0dnZyY4dO6blWNSSBLY6N9SXwSnocTNs2rI4/pG/BMNg2e2fwxeqfn2bEEIsVI52SkGnWnjK2lnydn60LIzZPsv6nJ0757H7lK8UhopBqLhu+k0Mv0E0EKXZbMbwu8Gp71QfK5auGP1c2WfHBqxq31e+Xuxn+Ax8SiYD6oEEtjpXvEN07CM9eu64g+yzz7L0X/4Fo6OjFkMTQogpKc48ZQtZd7HdMlPIuAGqkCVjZ8gVchXt5f2O9B7hh50/rBrAxgazvJM/p/EWA4wZMN1yzBKPxKvWT9R/bFgy/ea4MFZcD/im/s9zZ2cnO/5gxzn9zKJ+SWCrc8VnsDV2jAa23L599N35JRpe/3oSr35VrYYmhJgnHO1UBqSJ1u3KcDW2Xyls2RN/x9lcxxRQAUKBEKFACG1pGlONo6EoYNJgNhD0Bwn5Q24ZCFVun6Z+oqAV9AdlZknUFQlsda6/O40ZDRCOBQH3f6hd//C/8cViLPrYR2s8OiHEbCmGqkwhQ7qQJlPIjC5WpmJ7yu2FzFmNyfSbhAJu6AkHwm4ACpiE/WESkUT1Nm+9om1Mv2I4K37G8Bmlfc6X65GEmCoJbHUu1VX50veh++4j/dhjdPzdJwg0NdVwZEKIaoqn/UasEdJW2i0L6dJ6+XZFcLIqA1XvQC+f/t6nzzpU+ZWfcCBMJBAhbIQJB9wlFozRFmkrbVdbigEqFAiNC2LlbTKhT8KwAAAgAElEQVQDJcTsOWNgU0o1T+J7HK11ahrGI8ZIdadZcYF7h6i2bXo/9zmCa9bQeMstNR6ZEPOD1pqcnSsFrHRhNFiNFEbIWJmK7bSVrh7Eytome3deQAVGg5IXqiKBCIlgAp/hY8WiFePai0skEKn62WKd4TPkMT9CzCOTmWE74S2n+5PvB1ZMy4hEST5TID2YL82wDf785+T27WfpP39WXjslFrziTNZQfogRa4Th/DDD1vBoWWV9yBpiJD9SqhvJu6HL1vak9hn0BYkaUSJGhIgRIRqIEgvGaI+2EwlESm1RI1q5HSirL9sO+oMT7quzs5MdV+2YpqMlhJjrJhPYdmutLz5dB6XUU9M0HlGmv3v0HaJaa/ru/BLBdWuJv/rVNR6ZEOcub+cZtAc5mDrIYH6QwfxgKXyNK63xoWzIGqLgnHkmK+QPEQvGiBnuEg1GaQ23EjXcsBUJREphKmpEiQaihI1wab08oJVfSyWEELNpMoHtimnqI6YoVRbYMk88QW7vXjr+4e/lAbmiLmityRQyDOYHGcgNMJQfKgWvwdwgQ9YQg7nRIFasL25n7az7Rceqf79P+YgZMeLBuBuuDPfaq9XG6lJdeVvMiBELjvaPG3GiwaiELCHmM8cBOw92DmwLCjlv21sKeRpSL8BBJmj3Plft863rYfv/qPVPWHLGwKa1zhbXlVKG1to6XR8xfVLdaZRP0dAW5uSn78LX0EDDzTfXelhinikGr1Qu5S7ZVGl9IDdQsV4RyPJDZ7xWK27ESZgJEsEE8WCc1Q2rSZiJUn3X4S4uPf9SEsFEqb44GxYOhOUaLCHqhV2AQnY05JTWs1MIQhOHqvGhq/jZ8vYq3z2JWfaLAZ6e5M+pfOA3wR+EddfNrcBWpJT6d+CNSqkR3GvadgG7tNafm6nBLXT9XWkSLSGcZC9Dv/oVzW9/G75wuNbDEnVMa81gfpBkNlkRtk4XwlK5FJYz7v9hJXEjToPZQKPZSIPZwNLY0lLAKgaxUuAqrgcTxIwYft/pr7XsTHayY82O6T0IQswXWo8GnUJuNOgUsmXbo22Lup+Gp45Vbav+uXxl8CptV2nTzvT9XL6AG4iKS8AEv+EFJcPbDkIwNqYtCIHgaT5bvf2Z53Zz0SXby9rH7Kf8u87wd1YtTeWxHlcB7VprSym1FLgI2DIzwxLgzrA1dkQY+tnPwLZpfNObaj0kUQOZQoZkNkkyk3TLbJK+bB/JbJL+bH+pLplJkswlJ7yuK6ACFcFrRXwFjW2NpbpifZPZVFpvMBvO6onrQswLWruBxcp4QccrJ72dHV2sbGWAqpitmiBc2VN7U8NmgN1VGnwBCIS8UBIaDS2B4hKCUOPodkXbmO2K9WBZXXCSocqY9VDUf8KEVS+b1X3OhKn8TfwI0ASc0lofB44DP52RUQm0oxk4lWbZpiYGvvUTQps3Y65ZU+thiWmSLWTpyfTQm+mlJ91TWu/L9I0GMG+Z6Plb4UCY5lAzLaEWOiIdbGreRHOomeZQM02hJppCo8Gr0WwkZsTkFKOYm7T2Qk+GYK4PkgfPLThNZftclAJSqDIcFUOMEa4MSmcMR6ExbV6d1+/RJ5/m8iuuGg1lxTa57nlemEpguxN4QCn1H8CjuKdDByb7YaXUl4GbcQPfBV5dM/AtYBVwCHiL1rpfKbUD+AHwkvfx/9Ja//0UxjrnDfVnKVgO8UCG7LPPsugjH6n1kMQZaK0Ztobd8JXurQhkz/c8zzd+8Y1S25A1NO7zARWgOdxcCl0rEytL682hZlrCLTSZTTSHm2kym4gYkSqjEGIWOTZYaTfkWGkvIGXcsrh9Vm3Zyn5l/2m5EuDhKY7T7wUloxh4wqOBKWBCqKFye2z7lLbLAtoszyRlIr3QtHJW9ylmz1QC2zeBr3ufeR+wRSkV0lqvneTnvwrc7n1H0ceAX2utP62U+pi3XXzf0oNa6wV7hX3xDlFjv/vElMRrbqzlcBY8rTVD1hDdI910jXTRle6ia6TL3U67ZXe6u+psmOk3iakYy+3lrGtcxxWLr6At0kZruJW2sFdG2mg0G+XJ8WL6OI4bePIjYI1APl22nfa2R0a3JwxQ5SFqTNsUT9mVBMJuuKlYIm7QCTePbwuE3bBlRNh78DDnbb5o8gFKZpjEPDGVwHZMa/2J8gqllDnZD2utf6uUWjWm+vXADm/9a0Ano4FtQSsGNt/OBzAvvBBj8eIaj2h+sxyLrpEujg8f58TwCTeUeUt32g1p6UK64jM+5aM13EpHtIP1Teu5atlVLAovojXiBrG2cButkVbiRpwHHnhA3n8oxjtNqGrpfQx29UB+eEzASo8PYONCWLpiVmpSlM8NTWMDlBFxT9vFq4Sriv5naotMS4A6me3kvK07zvrzQsxVSms9uY5KfQ7Yr7X+17PemRvYflx2SjSltW4sa+/XWjd5p0S/h/uEphPAX2itn5/gO98FvAugvb390nvuuedshzdpw8PDxGKxGd3HyZ0OqZc0V//6g6RfcyMjr33tjO5vOszGcTlbjnYYsAfoK/SVlmQhWVpP2Sk0o38WFIq4P06Tv4nGQGOpbPS7602BJhL+BH41uVMe9XxsammuHRflWAQKGfz22CU7pj6L385U6ZstK7P4nanNUDnKj+MLYftD2H7TK0M4PrNK3UR9Rj9f/l1aBWAOXOM4137PzCY5NtXV83G55pprdmqtt02m71Rm2NqB65VSHwWeBJ4BntZaf+csxngmTwIrtdbDSqnXAN8H1lfrqLW+E/f6OrZt26ZnYxajs7NzxmdLfvD0U6hoHz6tueCtbyW8deuM7m86zMZxOR3bsTk5cpLDg4dHl6HDHB08yomRE+PunlwUXsTShqVsjG1kaWxpaVkSW0J7pB3DP30PXK31salXM35ctHZnn3KDkB2E3BDkBtwyPwK5YcgPeeWIO5uVG3LLUru35IbhNI8/qeALuI8kCMbAjEEkBsEOMOMQjI4uRhSCEXf2KRitKJ94dg/brrjaq4uAEcUXCOJjan9xzzfyZ2licmyqmy/HZdJ/7rXWb4HSadDzgQuBy4FzCWzdSqnFWuuTSqnFwClvX4Nl+/2pUuoOpVSr1rr3HPY1p6S60zSlu/A3NhK68MJaD6duaK3py/bx0sBLHBo8xJHBI6Xy6NDRiueJhQNhViVWsbF5I9etvI5lsWWlULY4thjTP+kz+qIWCnk3aFWELa/MDo4Gr/K2in7e+mSeH+UPeiEq7gasYtiKd7h1wWhlvTmmHLseMM95tmr4kOU+aV0IIZjif9S8Nx3kcGfAnpyG/f8QeBvwaa/8gbefDqBba62VUpcBPqBvGvY3J1g5m+H+HO3du4m+7GUL9kXvfZk+DqQOsD+1nwOpAxwYOMCB1AFSuVSpT9AXZHl8OasSq3jFslewMrGSFYkVrEqsojXcKo+xqKVCHrIpyKQgO1C2niqtn/fSHuj+d699TCibzCMV/CaEEmAm3NmrUAKiqyu3x6172xUBa+KXsAshRD2YtTcdKKXuxr3BoFUpdQz4BG5Q+7ZS6p3AEeAWr/ubgfcqpQpABvgjPdmL7eaB4g0H4d6XiNz6hhqPZubl7Tz7U/vZk9zDnuQe9vXv40DqAP25/lKfuBFnbeNarltxHWsb17K2YS0rG1bSEek449P0xTmwMpBOjglbA2dez6TOfNF7IEyzLwx2u/tYhdgiaFnrBqpQMVQ1lIWteFmbtx2QWVIhxMIwa2860FrfOkHTdVX63o77CJAFqRjYIuluIpdeWuPRTK8Ra4S9yb3sTu4uBbT9qf2l68sigQjrm9Zz7YprWdOwhnWN61jbuJZFkUUyW3YutPYCVdINYOnk6HrVst8trfTpv9dsgHCDG7hCjdC6bnQ93OiWFesN3rr73KuH58m1JUIIMdPkTQd1qL87DWhiwTzBOfx2A0c7HEgdYFfPLnb17mJXzy4OpA6U7sZsDjWzqXkTL9v8Mja2bGRT8yaWx5fLs8gmo5CDkV4Y6SkreyDdCyN9VUJYP2h7gi9TEG6CSLP7DKzEUmi/0N0u1oXHBK9wozvLJbObQggxK2btTQdi8lLdacLWALFLLppTs0pZJ8tDxx/iqVNPsatnF8/1PsewNQxAg9nAltYtvHLVK9ncvJmNzRtl1qycY3szX2Xha6Taeq+75Cb4o+cPQqTVC1pNsGiTG7iKwatUtoz2CTXKg0WFEKLOzeabDsQk9R8fIDx0gsiVl9R6KKeVttI8deopHu96nMe7Hue53udwjjr4lZ8NTRu4ac1NbGnbwpbWLaxMrFyY4ayQg6EuEgN74IVBGO6GoS4Y7oKh7tEy3Vv9bkblc8NVtA2irbBkqxvIitvRtsp1Mz4nnqUlhBBiambtTQdicrTWpE5l6Eh3E77kjbUeTgVHO+zu282Dxx/kd8d/x3O9z1HQBQIqwAWtF3BD4gbedPmbuKjtovn/nku74IavwePecnJ8CBvuck9FApcAPOV9Vvkgugji7RBfAksuhlh7lRDW5s6AyWlHIYRY8KYS2J5WSn24/E0H3iM+xDQaSeUpFBSRbA+hzZtrPRwG84P8/vjvefD4gzx0/CGS2SQKxfkt5/P2C97O9o7tbG3bSsSI0NnZyRVLrqj1kM+dY3th7AQMHHPLweOV60Nd468J8wfd4BVrd+92XHml+xyveAe7XjrFlitvgFiHG8okhAkhhJiCen3TwYKV6h4BoKEpgM+szQRmKpvi/qP384vDv+DRE49S0AUazAauXHIlVy29ipctfRnNoeaajG1a2AU3dKUOQ+oI9B8eXR84BkMnYcxbEQiEoWEpJJbA6leMrieWeeUSdzZsgtORycFOWHzRzP9sQggh5qVzfdPBZZzbmw7EGP1d7mMUWta0zep+R6wRfnnol/zspZ/xWNdj2NpmaWwpb938Vq5dcS0Xtl44d553pjWk+6DvgBvE+g9D6tBoOBs8XhnIlM+9M7JxJax6uRfAlrpLg1eeJowJIYQQM23Kr6Sb5jcdiDGSh3rx2zkaz5/5ezkc7fBE1xN8f//3+dWRX5EpZFgRX8E7LngHN6y8gU3Nm+r7RoHsgBvK+g5A8gD07R/dHnsXZazdDWTLL3PLxhXQtNJdb1gG0/jeUCGEEGK6nTGwKaWe1Fqf9nbFyfQRk5M8lCSc7iZ8waSfSTxlA7kBvvvid/nOi9/h+PBxYkaMm9bcxBvWvYEtrVvqK6Rp7Z6m7NkDp3ZDz143mCUPuI+5KFHQuBya18KWW6BlnbvetMqtN8K1+gmEEEKIczaZGbZNSqldp2lXQMM0jWfBG+jLE82cInTeedP+3fv793PXnrv48YEfk7WzbO/Yzgcv/iDXrbiOUCA07fubEq3dU5Wn9kDPbq/c4wa0/NBov+giaN0A5904Gspa1rnBzKjxzyCEEELMkMkEto2T6DPRI9TFFBQsmxHLoMPM44tGp+17d/ft5gvPfIH7j96P6Te5ec3N3LrxVs5rnv5QOCl2Afr2wcld0LULTj4DXc+676AsirZB20bYeqtbLtrklpE5fLODEEIIcZYmE9i+DHxIa/08gFLqdbjvEP2l1vqxmRzcQjNwKgMomhZNzzPM9ib38vmnP8/9R+8nbsR570Xv5daNt9IUapqW758Ux3ZPZR573Atmu6D7eShk3fZACBZthvPfAO0XeMFsE0RbZm+MQgghRJ2bTGBbVhbWrgS+AXwL+KpS6m+01vfO5AAXkr4DpwBoXrvonL6nN9PL7U/dzn/t+y9iwRjv2/o+btt0G4lgYjqGeXrDPXD8CTj6mBvSjj8JlvuoEswGWLwFtr3TLTu2uKc3/VO+90UIIYRYUCbzL+Vg2fqfAP+mtf6oUmoR8ENAAts06dtzHIDWC1ed1ecLToG7dt/FF575ArlCjj/e/Me8e8u7aTBn8BLD/sNw6CE49BCX7/0NdHa59b6AO2N28W2wbDssvRSa18ijMYQQQoizMJnAtl8p9Wbgt8AbgDcCaK1Pyaupplfy2ABmLkds81VT/uyB1AH+1+/+F8/2PsvVy67mI9s+wqqGVdM/yIHj8NJv4dCD7pI64taHmxmObSB81ftg2WXuQ2KD8/z1VEIIIcQsmUxg+3Pc06B3A/dprX8PoJQygNgMjm3BGegvEMn1EVg0+Yfmaq35+gtf51+f/FeiRpTPXP0ZXrXqVdP3aA7bgiOPwP77YN99cOoFtz7c5D5k9ooPuGXbJp7/7W/Z8bId07NfIYQQQpScMbBprbuAG5RSPq21U9Z0DXD/jI1sgdFaM5Q3WRKyJh22BvODfPyhj3P/0fu5dvm1/O0Vf0tLeBou1s+kYO/PYO9P4OADkBsEnwErr4Ab/gHWXgOLzgef79z3JYQQQogzmsqrqZwx278EfjntI1qg0oN5CsqksWVyT9w/mDrI+3/9frpGuvjo9o9y26bbzm1WLZ2EPT+BF34ABzvBsSC+BM7/Q1j/SljzCjDjZ//9QgghhDhrcntenejbewKApuVnfuTGk91P8sHffBDDZ/CVV3+FrYu2nt1OCzl48efw9H+6pzu17b6q6Q/eC5vfAEsvkZsEhBBCiDogga1O9L5wFIDWTctO2++h4w/x4d98mCWxJXzh+i+wLH76/lWdfAaeugue/TZk+iG+GK78AJz/RvdmAQlpQgghRF2RwFYn+g8n8dkhWrZO/A7RR04+wp/d/2esbVzLF2/44tQegGsXYPcP4dF/g6OPgt+EjTe5j91Ycw34/NPwUwghhBBiJkw6sCn3AqnbgDVa679XSq0AOuRtB9Mj1Zcnkk9jtLVWbX/61NN86DcfYnl8+dTCWm4InvgyPHonDB5z37n5qk+5r3wKz+IbD4QQQghx1qYyw3YH4ADXAn8PDAHfA7bPwLgWnKFMgIQxXLXtxPAJPnz/h2kNt/KlV35pcmEtNwSP3Qm/vx0ySVh9NbzmM7DhVTKbJoQQQswxUwlsl2utL1FKPQWgte5XSgVnaFwLSsGySfviLI8PjWvLFDJ86DcfIm/n+cqrvkJruPoM3OiX5d2g9uD/da9PW/9KeMXHYNmlMzR6IYQQQsy0qQQ2SynlBzSAUqoNd8ZNnKP+PUfRyk/TsvGvkPrsE59lb/9e7rjuDtY0rpn4S7R27/j8xd9A8gCsvQ6u/Rv3lVBCCCGEmNOmEtj+P9z3hrYrpT4JvBn4+IyMaoHpedZ9vVPz+vaK+geOPsC39n6Lt21+G1ctO83rqgZPwI/+DPb9wn2Z+m3fhfU3zOSQhRBCCDGLpvLg3LuUUjuB67yqN2itd8/MsBaW5OE+oIG2i0Zn0Abzg/zt7/+W85rO40OXfKj6B7WGp++Cn/812Hl45Sfh8neDf3IP3xVCCCHE3DCVu0T/55iqG5VSVwI7tdZPT++wFpaBngxBy0dk6aJS3eef+jypXIovXP8Fgv4qlwpmB+AH74fdP4KVL4PXfQ5a1s7iqIUQQggxW6ZySnSbt/zI274JeBx4j1LqO1rrf5ruwS0UA8M+YgyVXi31Yv+L3LP3Hm7ZcAubWzaP/0DXs/Ctt0LqiPtuzys+IO/1FEIIIeaxqQS2FuASrfUwgFLqE8B3gauBnYAEtrM07ERZGkmWtj/31OeIGlE+ePEHx3fe+3P47jsg1ABv/4n7QnYhhBBCzGtTmZZZAeTLti1gpdY6A+SmdVQLSLpnACsQpaEtBMDzvc/TebSTt21+Gw3mmLtGd34V7rnVvbHgXQ9IWBNCCCEWiKnMsP0n8IhS6gfe9muBu5VSUeCFaR/ZAtHzzEEAmle4D8P9t2f+jQazgds23VbZ8bEvwU//AtZdD7d8DczYbA9VCCGEEDUylbtE/0Ep9VPg5YAC3qO1fsJrvm3iT4rT6dvXDQRp2biUI4NHeODYA7z7oncTC5YFsp1fc8PahhvhLV+HgDyvWAghhFhIpnql+kHgYeBJIKKUunr6h7Sw9B8fRDkFWras5e49d+NXfm7ZcMtohz0/gR992J1Ze8vXJKwJIYQQC9BUHuvxP4APA8uAp4E/wA1v187M0BaGgf4CkUKKvOHwg/0/4IZVN7Ao4j3eo+s5+N6fwpKt8JZvQMCs7WCFEEIIURNTmWH7MO6L3g9rra8BLgZ6ZmRUC8hgLkjMyPKbo79hyBriLRve4jZkB+Ce/wdCCfijuyEYqe1AhRBCCFEzUwlsWa11FkApZWqt9wDnzcywFga7YJP2JWhoUPzspZ/REe3gkvZL3MaffgQGjrnXrCUW13agQgghhKipqdwlekwp1Qh8H7hPKdUPnJiZYS0M/XuPo30Bom1Bfn/897z1/LfiUz547nuw61vwio/B8stqPUwhhBBC1NikAptyH8H/Ia11Cvg7pdT9QAPw85kc3HzX+4L70veTjT0UdIGbVt8EmRT89C9h6aVw9UdqPEIhhBBC1INJBTattVZKfR+41Nt+YEZHtUAkX+oFEjwSfZaV4ZVsaNoAP/tLyCThrf8F/qlMgAohhBBivprKNWyPKKW2z9hIFqBUdwbDGuE3hce5aulVqN4X4fF/h23vhMUX1Xp4QgghhKgTU5nCuQZ4t1LqMDCC+/BcrbXeMiMjWwAGhyDs9JPXFlctvQo6PwVGBHb8Va2HJoQQQog6MpXAduOMjWKBGnKihAIvEQ6EuVSF4fl74aq/gGhLrYcmhBBCiDoylVdTHVZKXQRc5VU9qLV+ZmaGNf9lU8PkAzHS6jjbO7ZjPnwHBONwxftrPTQhhBBC1JmpvOngw8CfAv/lVX1TKXWn1vpzMzKyea7nmZcAOBo+yhXNV8Pv/hdsfydEmms8MiGEEKK2tNbYjiZvO1gFTc62sWxNvuBg2Q75gkPeKy27vK6yj2U77DmQ5+nCi2M+qyv6jP8+zSUrGvl/X39BrQ9FyVROib4TuFxrPQKglPo/uK+mksB2FvpePAkEONp8iv9+6iVwLNj+p7UelhBCiAXE8UJRruCQK9hucCm428UQYxUccl6ZP0M4yo/pUwxZ+TF9LG+flu2UQlm+bB9520HrafxB9+3D8CsMv49gwOeW3nrQ78MIKLf0+4iaAYJ+H83R+nod5FQCmwLssm3bqxNnof9oCqWb6FqcZsNzP3Rf7t66rtbDEkIIMcO01hQcXRGMimEpN0FdMfTkLNsrnVLQKvY5cizHd088Oe7zFf2t0fCUK7izVtNFKdwQVB6KAj4MvyIY8BP0KzcgBXzEQoFSe/EzbmjyYwQUpt9XGa4CE3zvuDof5pg+j/z+Ia7bsQOfb25HlqkEtq8Ajyql7vW23wB8efqHtDCkkhbBfJLzWpcROPoCvPpTtR6SEEIsGLajyVo2uYIz6TJXVmaLpeUGn/LgVAxDuXFhyS6FrOmYPVIKTC/ImIYfx7JpsAbd7YAPM+AnEgzQ5IUk0yuDXltFnfcdpt+HaZTNPpVmoMbOSPnGhDNFwD+VJ4XNHsOn5nxYg6nddPDPSqlO4OW4M2vv0Fo/NVMDm+8GswaO3cXF2QwEY7Dh1bUekhBCzLqC7ZD1QtGZglLWsnn2sMW+3x6cUtjKV6kvOGefmJSCUMCPafhKpVkehPw+mqLBUpgpD0fmmPBUbAuWBSXT8Jd9djR8Bcd8NuBTuC8icnV2drJjx45p+FUR9WgqNx38H631R4Enq9RN9ju+DNwMnNJaX+DVNQPfAlYBh4C3aK37vddh/SvwGiANvF1r/WS1751rHNshrRJY/me4/NjzsOm1YIRrPSwhhADcU3Z52yGbd8gWbDJ5m4zlLsVglcm74SdT2ra9vg4ZyyZnjX7GbXPI5sv6e+VZnZLbvRsAn4KQ4Sdk+DEDvlJpemUibNAWN8va3eBzptI8TX0o4MfwVwYlIWbDVE6J3gCMDWc3Vqk7na8CtwNfL6v7GPBrrfWnlVIf87Y/6n33em+5HPiCV855qQMncPxBUsFTXDjcDxe8udZDEkLMIZbtkM67QWgkXyCTt0nnbdLl65ZNJl9w+1k2OcupHry803qZfGWQOpsJKL9PEfYCVDjohptw0N1uCBt0JEzChltnem1ufzdshcrCUvnsVak0/Ox89BGu3XFV6TolIRaKMwY2pdR7gfcBa5RSu4rVQAz43VR2prX+rVJq1Zjq1wM7vPWvAZ24ge31wNe11hr3tViNSqnFWuuTU9lnPep9zn3pe64hRcSIwuqrazwiIcR0K4aqdDE0TRCq0rnRUFXe91hXljv2PDw+kJ3FrFQw4CsFo2KgChluWFoUNwgbbiAKG/6ywDXap/S5oL8UwsJlbSFvezYC1D5TETPlPcti4VH6DFc+KqUagCbgU7izX0VDWuvklHfoBrYfl50STWmtG8va+7XWTUqpHwOf1lo/5NX/Gvio1vqJMd/3LuBdAO3t7Zfec889Ux3SlA0PDxOLxc7684M/3c/RwTX0Lv8r3tfQwPMX/PU0jq52zvW4zGdybKqrl+PingKErA3ZgiZra7IFRsuCLmsb06egyY1py9kw1TN9hg9MP5h+hRkAA8edifKrUn3QK0MBSvXBsna3hFBgtK/pB988On1XL79n6pEcm+rq+bhcc801O7XW2ybTdzL/TdkAHNVa3wqglPoT4E3AYaXU351NaJukan/DjPsrUGt9J3AnwLZt2/RsXHB5rhd2/vhHx/AXMqw0j9P2B+9hx7az/656Ihe8TkyOTXXnclxsRzOcKzCcKzCUtRjOFhjKFhjKFRjxlmFv9mq4VGe7Zd5rK9ue7CnAYMBHzAwQNf1EgwGi8QBtZoCY6d6RFw36CZdKty5SWi8u5XUBwoYf/5i72OT3THVyXCYmx6a6+XJcJhPYvghcD6CUuhr4NPBBYCtuUDrXC7C6i6c6lVKLgVNe/TFgeVm/ZcCJc9xXXRgYdAhY3Wwu5GH9DbUejhCzSmtNruAwmLXoGnF45miKoWyB4ZzFYLYwGryylhfG3BBWEcqyFiN5+8w7A0JGMWAFiA5NlxIAACAASURBVATdYNUcDbK8OUI06CdqBoiVtUW9vtGgG8piZoCIGSAWDBAxZ+e0nxBCjDWZwOYvm0X7b8CdWuvvAd9TSj09DWP4IfA23CD4NuAHZfUfUErdg3uzwcB8uH4NIO3EKLCXjfHV0LCs1sMRYsqyls1AxmIwYzHgLYNZi4G0xUCmMLrt9XFD12jgqnikwoPVL4WNBN2wFA8FiIUMEqEAixtCXp1RanMXw+3n1ZWHrrEzV0IIMRdNKrAppQJa6wJwHd71YlP4fIlS6m7cGwxalVLHgE/gBrVvK6XeCRwBbvG6/xT3kR77cR/r8Y6p7Kte5YYy5I1GbNVNQm42EDWitSZj2fSnLVLpvBesCuMDWMYaE8wKDGYt8gXntN8fDbp3BSa8ZUljiHgoXhbA3JB17OA+Lrv4wooAlggZRE1/3T6EUwghamEygetu4AGlVC+QAR4EUEqtAwamsrPidXBVXFelrwbeP5Xvnwt6nzkIgD90Ela8pcajEfOB42gGMhb96XwpgI2Wo+vJkTyp9Gi/04Uun4JE2HBDV8gtFzeESYQD4+qLy2h9YNJhqzP7Ejs2tU/XoRBCiHnrjIFNa/1J7w7NxcAv9ehtpT7ca9nEFJzcfQgI0xA5Bsv/oNbDEXUoa9n0jeTpG87RN5yndzhH34gbuPqG86UglkpbJL3ZsYlu9vb7FI1hg8aIUbpua8uyBv7/9u48OqrzzPP495FKC0IL2oUQGDD7IgwGBWJDYxsbW3bsdibjmKGn43Z37J5wmHTP9ExIMseTzCQ5dg/dYWY8JyEdJp10YtzBScZJd0c2ccCRDTYB21iYTZjFCIQkNlGS0Fb1zh91wcKSjLbSLUq/zzl16t73bk89XKoevXfLTktmTFoy2WmRaZkfKbzSkwNx8SgXEZF40adDms65N3poOzz04cS/2mN14CYwvqAFssb5HY4Mg85QmMY2x8Ezlz4swJraOdd8pSCLDF8pyJraOntcT0oggdzRXqE1OomxY0aRnZZETtqHbZEi7EohlkxmakB3ZBcRiQO6++AwC57vJLn9PNMmlvodigxSa0eIhmAb9cFW6i+1UXeplfpg24evS600BNs439Ie6QHbVnnN8okJRs7oZHJHJ5OXnsKEnDRyR6eQmx5py02PDOd5bWnJiSq+RERGKBVsw6y1bRQWqqfoJl1wEKs6Q2Hqgm2cabzM6YutHxZiXQuyS61cau3eE5aYYOSnp1CQmUJJ9ijmT8gmPyOFC7UnWDJ/ztVCLC89mczUJB12FBGRPlHBNozC4TBtgVwC7hA2Ts8P9UMo7GgItnG68TK1F1upbbxMbWPk/fTFVs40tlIfbO12E9XkQAKFmSkUZKQytSCd227OpSAzlfyMFAoyIu0FmSnkpCX3WIRt336a5XPHDtOnFBGReKOCbRgFT9QTTkwlkTNQMNPvcOJSZyhMbWMrJy+0UHPhMjXnWzh54TI1F1o4deEydcE2Qh+pxkYlJTJ2TCrFWaO4fWoexVmpjB0ziqKsSFtRZiqZo3QumIiI+EcF2zA6/tZ+AEaPuQSBFJ+juTE55zjf3M7xc818cL6FmvOXOXmhhZPee21j6zUFWYLB2KxRlGSPYvHNuRRnjWLsmFTGZqUyNmsUxd6tKlSMiYhILFPBNoxOHzoB3MTYm1L9DiXmXWrt4PjZZo55r67DHz13LD8jhfHZo7j1pmxKskcxPjuN8TlpjM9OY+yYVD1KSEREbngq2IbRxdNBEkJtTCud53coMcE5x5lLrRyua6K6Lkh1XRNHzzZx7GwLZ5vars5nBsVZo5iUN5oHbylmUl46k/LSmJCTRkl2GqlJiT5+ChERkehTwTaM2lqSSQrVUThxqd+hDKsrhVl1XROHvcLscH2QI3VNBLvccywvPZnJ+encNaOASfmjmZg7msn5o5mQo6JMRERGNhVsw6iDHBLdMawofh9JFQo7jp1t5r3Tjbx3+tLV94stHVfnyR2dzNTCdB5eMI6phRlMLUhnWmEGOaOTfYxcREQkdqlgGybtza20J+WQbr+H5DS/wxkS4bDj/YYm3v7gIlWnGtl58DKnXnmJyx0hAJITE5helMG9s4uYVZzJ1IIMphWmk5uuCy5ERET6QwXbMDn61l6wBNKyWvwOZcAutrTz9smLvH3iAm+fvMg7H1y8ekgzPSVAcRp8dtF4ZhdnMmdcFlMK0nXCv4iIyBBQwTZMTrxdBUym4KZ0v0PpszONrbxx9BxvHD3HrmPnOXq2GYjcKmN6USYP3lLM/AnZ3DJ+DJPzRvO7373K8uWzfY5aREQk/qhgGyYXj58HJjNlYew+Q7Qh2MaO9896Rdp5jnkFWkZqgE9MyuEzC0uYPz6b0pIsRqdo1xERERku+tUdJq2NAZLCFyieFTtXiIbCjndOXmT7oXq2H2qg6lQj8GGBtvoTE1g8OZeZYzNJ1DMvRUREfKOCbZh0hrIJhOuwMeN9jaOlvZNtBxt46b0z/K66gYstHSQYLJiQzX9aOZ2lU/OYXZylAk1ERCSGqGAbBs45OgL5pHbuidwFdpgFWzv47cF6fl11hu2H62ntCJM7Opm7ZhSyfHo+S6fmMSZNt9QQERGJVSrYhkHDiRpCgTRSUi4N2zZDYUdldQMv7Knh5f11tHeGKchI4bMLx3Pf3LEsmpijXjQREZEbhAq2YXBkxw4gn4yi6Kf7aEMTW/bU8PO3aqi71MaYtCRWLRrPp+YVs2BCNgkq0kRERG44KtiGQf2BE0A+42bfFJX1h8OO31U38IPXj/Pq4QYSE4zl0/L52qdKuHNmASkBPdZJRETkRqaCbRg0NXSSEGpn+m33DOl62zvDvLCnhk2vHeX9hmbyM1L4yxXTWFU2noLM1CHdloiIiPhHBdsw6GjNIMk1kFawckjW19YZ4qe7a/jOtiOcbmyltCSLDZ+9hfK5Y0kO6MkCIiIi8UYF2zDotHwC7tSgrxANhx2/ePsU618+RG1jKwsmjOHpf1XK0ql5mA9Xn4qIiMjwUMEWZR0trbQn55Iafm9Q69lz4jz/7Vf72VvTyLySLP76M6XcPkWFmoiIyEiggi3Kjr7xGlgio3I6BrR84+UOvvnP+/np7hoKM1P420fm8Ye3jNPVniIiIiOICrYoO/FWFTCP3Ek5/V5228F6vvzzKuqDrfz5H9zM2jun6BmeIiIiI5B+/aPsQk0QgJuXlPV5mfbOMN/6lwP8/Y7jTC/M4Ht/fCulJWOiFaKIiIjEOBVsUdYWTCHgLjF+zqf6NP/pi5f5wk/e4p2TF3n8tkl86b7puo+aiIjICKeCLco6QjmRh74nXr/o2neqkcd+8HtaO0J8Z/UC7ps7dhgiFBERkVingi3KOhMLSLF9153v1cMNfOHHexiTlszmz3+CqYUZwxCdiIiI3AhUsEXRhQ9O0pmUQUag6WPn23aonid/tIebC9L5+z9ZRKGeUiAiIiJdqGCLosO/2waUkFmU1Os8r1Wf5cl/2MPUwnSe+7PFZKX1Pq+IiEgs6ujooKamhtbWVr9D6SYrK4sDBw74GkNqaiolJSUkJQ38N14FWxSdOXQSKGHsnEk9Tn/vdCOf/9FuJueN5sd/+gkVayIickOqqakhIyODiRMnxtwN3YPBIBkZ/p1m5Jzj3Llz1NTUMGlSz/VAX+jBk1HU1BDGwp1MW3p3t2kNwTY+/8PdjElL4kePl5E9OtmHCEVERAavtbWV3NzcmCvWYoGZkZubO+jeR/WwRVF7aybJNJCRW3JNeyjsWPPcW5xvaeeFP/8kBTpnTUREbnAq1no3FLlRwRZFIfJIDNd3a//uq++z69h5/uZfz2POuCwfIhMREZEbiQ6JRklnWwftyfkkply4pv3dmot8e+th7i8dy6cXjPMpOhEREbmRqIctSo7v2oFLCDAq48OHvofCji//vIqc0cl86w/nqvtYRERE+kQFW5Sc2PM2UErOxA+fAfrcmyd47/Ql/veq+boiVEREJA5MnDiRjIwMEhMTCQQC7N69OyrbUcEWJec/aARgctkCAC61dvA3Ww/zyZtzeaBUj5wSERGJF9u2bSMvLy+q21DBFiWXG1MI0MSE+fcC8IPXjnOxpYOvlM/UoVAREYlbX//Ve+w/fWlI1zmrOJP/+qnZ152voqKCdevWAZCSksLOnTuHNA4/qWCLks7ObALUkRhIprGlg++/dpSVswt1VaiIiEiUrF27lsrKSoqKirpNW7p0KcFgsFv7+vXrWbFixYC3aWbcc889mBlPPvkkTzzxxIDX9XFUsEVJR2IBKeGDAPxk1wmCrZ188a5pPkclIiISXX3pCYuW8vJy5s6dy+rVq9mwYcM10yorK/u1rhUrVnDmzJlu7d/85jd56KGHro6//vrrFBcXU19fz913382MGTNYtmzZwD7Ax1DBFgWNtafpTM4iPSFIZyjMj3ee4LYpucwqzvQ7NBERkbi0Y8cOnHPU1tYSCHQvb/rbw/ab3/ymT9stLi4GoKCggIcffphdu3apYLtRHH71t0AxGYUBtu6v43RjK1970L+/OEREROLdli1bmDZtGoFAAOccwWCQzMwPO0r628PWF83NzYTDYTIyMmhububll1/mqaeeGvLtgG6cGxVnDp4AYOysCfz4zROUZI/irpmFPkclIiISv1atWsXGjRspLS1l8eLFVFdXR32bdXV13H777cybN4+ysjLuv/9+7r333qhsSz1sURCsD4ELkXvrH7Dj7w7z7++cSmKCrgwVERGJlrKyMqqqqoZ1m5MnT2bv3r3Dsi31sEVBe0sGyW3n+O2pRJyDB28p9jskERERuYHFRMFmZl80s31m9p6Z/YXX9jUzO2Vm73ivcr/j7KtOckl09fxq72nmjsvi5vx0v0MSERGRG5jvBZuZzQE+D5QB84AHzGyqN/nbzrlbvNe/+BZkP4Q6OulIKiAxcJ69NY08OE+9ayIiIjI4sXAO20zgDedcC4CZvQo87G9IA3dizxuEE5MJp7YAcPcsXWwgIiIig2POOX8DMJsJvAgsAS4DrwC7gXPAY8Alb/w/Oucu9LD8E8ATAIWFhbc+//zzUY+5qamJ9PSeD3PWb3uFhro7CKb+M7/IuYtnlqVFPZ5Y8XF5GemUm54pL71TbnqmvPTOz9xkZWUxZcoUX7Z9PaFQiMTERL/D4MiRIzQ2Nl7Tdscdd+xxzi3sy/K+97A55w6Y2TPAVqAJ2At0At8B/jvgvPe/AR7vYfnvAd8DWLhwoVu+fHnUY96+fTu9beeFX70KwJ5AIffdchPLl4+c+699XF5GOuWmZ8pL75SbnikvvfMzNwcOHCAjI8OXbV9PMBiMidhSU1OZP3/+gJf3/Rw2AOfcJufcAufcMuA8UO2cq3POhZxzYeDviJzjFvNaGpNI7GxhT8I0lk/P9zscERERiQMxUbCZWYH3PgH4NLDZzMZ2meVhYJ8fsfVXZ0c2SR31JCWnsHhyrt/hiIiISByIiYIN+JmZ7Qd+BazxzlX7azOrMrN3gTuAv/Q1wj7qTCwAV8/8CWNITfL/mLmIiIj0XUVFBdOnT2fKlCk8/fTTfodzle/nsAE455b20PZv/YhlMJoa6ulIzibkLrBoYo7f4YiIiEg/hEIh1qxZw9atWykpKWHRokU8+OCDzJo1y+/QYqNgixeRh74X0JjcxhIVbCIiMhL9eh2cGeJHRBXNhfuu39tVUVHBunXrAEhJSWHnzp392syuXbuYMmUKkydPBuDRRx/lxRdfVMEWb07vPwoUcDI1jTUTxvgdjoiIyIiydu1aKisrKSoq6jZt6dKlBIPBbu3r169nxYoVAJw6dYrx48dfnVZSUsKbb74ZvYD7QQXbELpU1wEuzIXieWSkJvkdjoiIyPDrQ09YtJSXlzN37lxWr17Nhg0brplWWVl53eV7ujetmQ1ZfIOhgm0ItbeMJpnzTJ0xcu69JiIiEgt27NiBc47a2loCge7lTV962EpKSjh58uTVaTU1NRQXx8YjJlWwDaFOl0dCuJ7Skk/6HYqIiMiIsmXLFqZNm0YgEMA5RzAYJDMz8+r0vvSwLVq0iOrqao4dO8a4ceN4/vnnee6556IZdp/Fym09bnihUIiOpALCdpa5JVl+hyMiIjKirFq1io0bN1JaWsrixYuprq7u9zoCgQDPPvssK1euZObMmTzyyCPMnh0bR83UwzZEavbuIZyYSmtSkJvz9Zw7ERGR4VRWVkZV1eCvTi0vL6e8vHwIIhpa6mEbIsfeiFxF0poZIDEhNk5QFBERkfiggm2InD1+DoCEiVN8jkRERETijQ6JDpHm84kkWBvjF/2B36GIiIhInFEP2xAJdWST1F7HvEljrz+ziIiISD+oYBsioYQCzNUzOX+036GIiIhInFHBNgRaLp6nPTmHUOAiSYlKqYiIiAwtVRdD4PD2V8ASCKV3+B2KiIiIxCEVbEPgg6rDACSNy/M5EhEREYlHKtiGQGNtGwBjy27zORIREREZjMcff5yCggLmzJnjdyjXUME2BDpa0klqv0DpvAV+hyIiIiKD8Nhjj1FRUeF3GN3oPmxDwLk8EjvrGDdmlN+hiIiI+OqZXc9w8PzBIV3njJwZfKnsS9edr6KignXr1gGQkpLCzp07+72tZcuWcfz48X4vF20q2AYpHA7TkVRAoPMtEvRIKhEREd+sXbuWyspKioqKuk1bunQpwWCwW/v69etZsWLFcIQ3KCrYBun0vncIBdJIDDT7HYqIiIjv+tITFi3l5eXMnTuX1atXs2HDhmumVVZW+hTV0FDBNkj7K18D5hDIS/Y7FBERkRFrx44dOOeora0lEOhe3qiHbYRreL8BgPzZ032OREREZOTasmUL06ZNIxAI4JwjGAySmZl5dfqN3sOmq0QHqb0xQEKonQV33ut3KCIiIiPWqlWr2LhxI6WlpSxevJjq6uoBr2fJkiUcOnSIkpISNm3aNMSRDox62AYp3JlNkjUwtiDX71BERERGrLKyMqqqqga9ns2bNw9BNENPPWyDFE7Ix8L1mOkKUREREYkOFWyD0HqpkfbkPEi64HcoIiIiEsdUsA1C1W9fBkskIaPT71BEREQkjqlgG4Tj70Tu5Jx+U77PkYiIiEg8U8E2CM11kYe+T719qc+RiIiISDxTwTYI4cujSWpvZM6ChX6HIiIiInFMBdsgRB76Xk8gUWkUERGR6FGlMQidgUKws36HISIiInFON84doFMH3qUzKZ1AUpPfoYiIiMgQmThxIhkZGSQmJhIIBNi9e7ffIQEq2Abs3W3bgTkk5yf5HYqIiIgMoW3btpGXl+d3GNdQwTZA596vB6BwzhSfIxEREYkdZ771LdoOHBzSdabMnEHRV75y3fkqKipYt25dZJmUFHbu3DmkcfhJBdsAdVxMxBI7uXXFA36HIiIiIsDatWuprKykqKio27SlS5cSDAa7ta9fv54VK1ZcHTcz7rnnHsyMJ598kieeeCKqMfeVCrYBch3ZJIcayM3N8TsUERGRmNGXnrBoKS8vZ+7cuaxevZoNGzZcM62ysrJP63j99dcpLi6mvr6eu+++mxkzZrBs2bJohNsvKtgG6MpD30VERMR/O3bswDlHbW0tgUD38qavPWzFxcUAFBQU8PDDD7Nr1y4VbDeqztbLtCfnk9J5yO9QREREBNiyZQvTpk0jEAjgnCMYDJKZmXl1el962JqbmwmHw2RkZNDc3MzLL7/MU089Fc2w+0z3YRuAs9X7cQkBEjM7/A5FREREgFWrVrFx40ZKS0tZvHgx1dXV/V5HXV0dt99+O/PmzaOsrIz777+fe++9NwrR9p962Aag5WTkUGj6hFyfIxERERGAsrIyqqqqBrWOyZMns3fv3iGKaGiph20Awhcj79OXftLfQERERGREUME2ANaWSaCjiTm3LvY7FBERERkBVLANhMsj0FmHJSh9IiIiEn2qOAagM1AA6KHvIiIiMjx00UE/1VQfoDM5iyT00HcREREZHuph66e9v3kFgJS8RJ8jERERkZFCBVs/nX+/FoCi2Xrou4iIiAwPFWz91HExAVyIBXff53coIiIiMkLERMFmZl80s31m9p6Z/YXXlmNmW82s2nvP9jtOANcxhuT2s2Tn5fsdioiIiAyxxx9/nIKCAubMmdNtWkVFBdOnT2fKlCk8/fTTwxqX7wWbmc0BPg+UAfOAB8xsKrAOeMU5NxV4xRv3nbMCEkJ66LuIiEg8euyxx6ioqOjWHgqFWLNmDb/+9a/Zv38/mzdvZv/+/cMWVyxcJToTeMM51wJgZq8CDwMPAcu9eX4IbAe+5EN8V7U1N9ORlE9y5xE/wxAREYlZlT89zNmTQ3snhbzx6Sx9ZNp156uoqGDdukj/TkpKCjt37uz3tpYtW8bx48e7te/atYspU6YwefJkAB599FFefPFFZs2a1e9tDEQsFGz7gG+aWS5wGSgHdgOFzrlaAOdcrZkV9LSwmT0BPAFQWFjI9u3boxZow9HDJIZyCac2RXU7N6qmJuWlN8pNz5SX3ik3PVNeeudnbrKysggGgwB0tHcQCoWGdP0d7R1X1/9x1qxZw0svvURhYSEAzc3NhEIhgsEgK1eupKmpeyH5jW98gzvuuOOatqamJsLh8DXbPHLkCEVFRVfbcnNz2b17d5/iAmhtbR3Uv4/vBZtz7oCZPQNsBZqAvUBnP5b/HvA9gIULF7rly5dHI8yI5cvhcXhlaxZR3c4Navv27cpLL5SbnikvvVNueqa89M7P3Bw4cICMjAwA7vyj2b7EAPDAAw+wZMkSVq9ezYYNGwAIBoNkZGSwY8eOPq8nPT2dhISEq58JIDU1laSkpKtto0aNIiUl5Zp5Pk5qairz58/vx6e5lu8FG4BzbhOwCcDMvgXUAHVmNtbrXRsLxMyJY4lJMZE2ERER8ezYsQPnHLW1tQQC3X+nly5d2mNv2Pr161mxYsV1119SUsLJkyevjtfU1FBcXDy4oPshJioPMytwztWb2QTg08ASYBLwOeBp7/1FH0MUERGRGLZlyxamTZtGIBDAOUcwGCQzM/Pq9MrKykGtf9GiRVRXV3Ps2DHGjRvH888/z3PPPTfYsPvM96tEPT8zs/3Ar4A1zrkLRAq1u82sGrjbGxcRERHpZtWqVWzcuJHS0lIWL15MdXX1gNezZMkSDh06RElJCZs2bQIgEAjw7LPPsnLlSmbOnMkjjzzC7NnDd/g3JnrYnHNLe2g7B9zlQzgiIiJygykrK6OqqmrQ69m8eXOv08rLyykvLx/0NgYiVnrYRERERKQXKthEREREYpwKNhERERk055zfIcSsociNCjYREREZlNTUVM6dO6eirQfOOc6dO0dqauqg1hMTFx2IiIjIjaukpISamhoaGhr8DqWb1tbWQRdLg5WamkpJScmg1qGCTURERAYlKSmJSZMm+R1Gj7Zv3z6oJwzECh0SFREREYlxKthEREREYpwKNhEREZEYZ/F0RYeZNQAnhmFTecDZYdjOjUZ56Z1y0zPlpXfKTc+Ul94pNz2L5bzc5JzL78uMcVWwDRcz2+2cW+h3HLFGeemdctMz5aV3yk3PlJfeKTc9i5e86JCoiIiISIxTwSYiIiIS41SwDcz3/A4gRikvvVNueqa89E656Zny0jvlpmdxkRedwyYiIiIS49TDJiIiIhLjVLCJiIiIxDgVbP1gZvea2SEzO2Jm6/yOZ7iZ2Xgz22ZmB8zsPTP7oteeY2Zbzazae8/22s3M/peXr3fNbIG/nyC6zCzRzN42s3/yxieZ2ZteXv7RzJK99hRv/Ig3faKfcUebmY0xsxfM7KC37yzRPgNm9pfe/6N9ZrbZzFJH6j5jZv/XzOrNbF+Xtn7vI2b2OW/+ajP7nB+fZSj1kpf/4f1fetfMfmFmY7pM+7KXl0NmtrJLe9z9dvWUmy7T/srMnJnleePxsc845/TqwwtIBN4HJgPJwF5glt9xDXMOxgILvOEM4DAwC/hrYJ3Xvg54xhsuB34NGLAYeNPvzxDl/PwH4Dngn7zxnwKPesPfBf6dN/wF4Lve8KPAP/ode5Tz8kPgz7zhZGDMSN9ngHHAMWBUl33lsZG6zwDLgAXAvi5t/dpHgBzgqPee7Q1n+/3ZopCXe4CAN/xMl7zM8n6XUoBJ3u9VYrz+dvWUG699PPASkZvo58XTPqMetr4rA444544659qB54GHfI5pWDnnap1zb3nDQeAAkR+eh4j8KOO9/6E3/BDwIxfxBjDGzMYOc9jDwsxKgPuB73vjBtwJvODN8tG8XMnXC8Bd3vxxx8wyiXyxbgJwzrU75y6ifQYgAIwyswCQBtQyQvcZ59zvgPMfae7vPrIS2OqcO++cuwBsBe6NfvTR01NenHMvO+c6vdE3gBJv+CHgeedcm3PuGHCEyO9WXP529bLPAHwb+M9A1ysq42KfUcHWd+OAk13Ga7y2Eck7JDMfeBModM7VQqSoAwq82UZSzjYQ+ZIIe+O5wMUuX6xdP/vVvHjTG73549FkoAH4gXe4+PtmNpoRvs84504B64EPiBRqjcAetM901d99ZETsOx/xOJGeI1BeMLMHgVPOub0fmRQXuVHB1nc9/TU7Iu+JYmbpwM+Av3DOXfq4WXtoi7ucmdkDQL1zbk/X5h5mdX2YFm8CRA5bfMc5Nx9oJnJ4qzcjIjfe+VgPETl0VQyMBu7rYdaRuM9cT2+5GFE5MrOvAp3AT6409TDbiMmLmaUBXwWe6mlyD203XG5UsPVdDZFj41eUAKd9isU3ZpZEpFj7iXPu515z3ZXDVt57vdc+UnJ2G/CgmR0ncrjhTiI9bmO8w11w7We/mhdvehY9d+3Hgxqgxjn3pjf+ApECbqTvMyuAY865BudcB/Bz4JNon+mqv/vISNl38E6OfwBY7byTsVBebibyB9Be77u4BHjLzIqIk9yoYOu73wNTvau4komc+PtLn2Ma1BAyfgAAA/RJREFUVt45M5uAA865v+0y6ZfAlatrPge82KX9j70rdBYDjVcOccQT59yXnXMlzrmJRPaL3zrnVgPbgM94s300L1fy9Rlv/pj9q24wnHNngJNmNt1rugvYzwjfZ4gcCl1sZmne/6sreRnx+0wX/d1HXgLuMbNsrwfzHq8trpjZvcCXgAedcy1dJv0SeNS7ongSMBXYxQj57XLOVTnnCpxzE73v4hoiF8mdIV72Gb+veriRXkSuNDlM5Iqbr/odjw+f/3Yi3cXvAu94r3Ii59K8AlR77zne/Ab8Hy9fVcBCvz/DMORoOR9eJTqZyBfmEWALkOK1p3rjR7zpk/2OO8o5uQXY7e03/4/I1Vgjfp8Bvg4cBPYB/0Dk6r4Ruc8Am4mcy9dB5If2TweyjxA5p+uI9/oTvz9XlPJyhMh5V1e+g7/bZf6venk5BNzXpT3ufrt6ys1Hph/nw6tE42Kf0aOpRERERGKcDomKiIiIxDgVbCIiIiIxTgWbiIiISIxTwSYiIiIS41SwiYiIiMQ4FWwiIiIiMU4Fm4iIiEiMU8EmIjccM8s1s3e81xkzO9VlPNnMdkRhmxPN7LKZvdPP5T5pZl//mOmjvLjbzSxv8JGKSDzSjXNF5IZmZl8Dmpxz66O8nYlEnmIxJ0rrP07kDuxno7F+EbmxqYdNROKOmTV5PWIHzez7ZrbPzH5iZivM7HUzqzazsi7z/5GZ7fJ6ujaaWWIftrHFzJ41s9fM7ISZ3W5mPzKzw2a26SPz3e4N/8LMvmFmlV7P4IroZEBE4o0KNhGJZ1OA/wmUAjOAf0Pkmbh/BXwFwMxmAp8FbnPO3QKEgNV9WPdc4Khz7nbgh8AmIg/lngN82sxSvPnmEHl+4ZXhi865pcAX+rgdERECfgcgIhJFx5xzVQBm9h7winPOmVkVMNGb5y7gVuD3ZgYwCqj/uJWaWSowBtjgNV0GNjnnar3pLUC7N1+Sc67RzNKALODb3jIB4OKQfEoRiXsq2EQknrV1GQ53GQ/z4fefAT90zn25H+udDbzlnAt74/OA7wCYWQlw2isMZwP7uyyzxzkX8sZLgX39+TAiMnLpkKiIjHSvAJ8xswIAM8sxs5uus8xcYG+X8VLgXW94XpfhuV2G5wDv9LKMiMjHUsEmIiOac24/8F+Al83sXWArMPY6i83FK768w56jnHMXvGldC7G5HxnuWrDNQT1sItJHuq2HiEgf6LYeIuIn9bCJiPRNCMjq741zr+fKjXOBJCLn1omIdKMeNhEREZEYpx42ERERkRingk1EREQkxqlgExEREYlxKthEREREYpwKNhEREZEYp4JNREREJMapYBMRERGJcf8fLLWdTIWSyPIAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 720x432 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"timesteps = range(1440) # minutes per day\n",
"\n",
"from matplotlib import pyplot as plt\n",
"plt.figure(figsize=(10, 6))\n",
"def plot_S(timesteps, S0=90, Qi=0.2, S_max=100, rho_s=0.01, eps=5):\n",
" S = euler(timesteps, S0, Qi, S_max, rho_s, eps)\n",
" plt.plot(timesteps, S, label=f'$\\\\varepsilon={eps:0.4g}$')\n",
" plt.ylabel('Storage $S [mm]$')\n",
" plt.xlabel('Time $[min]$')\n",
" plt.legend(loc=0)\n",
" plt.grid()\n",
" \n",
"for eps in [-5, 0, 1, 5, 10]:\n",
" plot_S(timesteps, eps=eps)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Conclusion\n",
"\n",
"As we can see in the figure above, $S$ exceeds $S_{max}=100$ for every setting of $\\varepsilon$ but more so, if $\\varepsilon \\gt 0$. I guess there was a typo in the original publication of Clark et al 2008 and the smoothing function $\\Phi$ should read instead:\n",
"\n",
"$\n",
"\\Phi(S, S_{max}, \\rho_{S}, \\varepsilon) = \\frac{1}{1 + \\exp\\left(\\frac{S - S_{max} + \\omega \\varepsilon}{\\omega}\\right)}\n",
"$\n",
"\n",
"With this correction, $S < S_{max}$ for quite a long time if $\\varepsilon = 5$ and equals the blue line in the figure above.\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment