{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Finite Volume Burgers\n",
"\n",
"Copyright (C) 2010-2020 Luke Olson \n",
"Copyright (C) 2020 Andreas Kloeckner\n",
"\n",
"\n",
"MIT License\n",
"Permission is hereby granted, free of charge, to any person obtaining a copy\n",
"of this software and associated documentation files (the \"Software\"), to deal\n",
"in the Software without restriction, including without limitation the rights\n",
"to use, copy, modify, merge, publish, distribute, sublicense, and/or sell\n",
"copies of the Software, and to permit persons to whom the Software is\n",
"furnished to do so, subject to the following conditions:\n",
"\n",
"The above copyright notice and this permission notice shall be included in\n",
"all copies or substantial portions of the Software.\n",
"\n",
"THE SOFTWARE IS PROVIDED \"AS IS\", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR\n",
"IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,\n",
"FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE\n",
"AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER\n",
"LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,\n",
"OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN\n",
"THE SOFTWARE.\n",
""
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import scipy as sp\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from matplotlib import animation\n",
"from IPython.display import HTML"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def gaussian(x):\n",
" u = sp.exp(-100 * (x - 0.25)**2)\n",
" return u\n",
"\n",
"\n",
"def step(x):\n",
" u = np.zeros(x.shape)\n",
" for j in range(len(x)):\n",
" if (x[j] >= 0.6) and (x[j] <= 0.8):\n",
" u[j] = 1.0\n",
"\n",
" return u\n",
"\n",
"def g1(x):\n",
" return 1+gaussian(x)\n",
"\n",
"def g2(x):\n",
" return 1+gaussian(x) + step(x)\n",
"\n",
"g = g1"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"tsteps = 250\n",
" dx = 0.00609756\n",
"lambda = 0.95\n"
]
},
{
"data": {
"text/plain": [
"[]"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXRdZ33u8e/v6GiwrHmyNVm240nybMs4hAAhCZCkTE2BNrchbVYgDaGsTveWdrWFe0tX7+2l0NuZpjSktBCgTS71pUkhEEJwZtmx5UF2YseTLMmSNVuzdN77xznHuIplyfI5Z+999HzW0rKsvbX3b+vYj/Z533e/rznnEBGR4At5XYCIiCSGAl1EJE0o0EVE0oQCXUQkTSjQRUTSRNirE5eVlbnly5d7dXoRkUDas2fPeedc+eW2eRboy5cvp6mpyavTi4gEkpmdmmmbmlxERNKEAl1EJE0o0EVE0oQCXUQkTSjQRUTSxKyBbma1ZvYjMztsZofM7Ncus4+Z2V+Y2TEzazazbckpV0REZjKXYYuTwG855/aaWT6wx8yecs4dvmSf24HVsY+dwN/G/hQRkRSZ9Q7dOdfunNsb+3wQaAGqp+32QeBrLupFoMjMKhNerVy14fFJ/u+rreza38boxJTX5YhIEl3Vg0VmthzYCrw0bVM1cOaSv7fGvtY+7fvvB+4HWLZs2dVVKlftH58/yRe+d5QLY5MAFOSE+f2faeCjO2o9rkxEkmHOnaJmlgc8Bvy6c25gPidzzj3knGt0zjWWl1/2yVVJkO82t/G5XYfYuqyIb//KW/nGJ3ZSX1nA7zzezNNHznldnogkwZwC3cwyiYb5151zj19ml7PApbd9NbGviQf2n+njt769n8a6Yv7+nkbesqKEG64r46v37qChqoBPf+NVjnYMel2miCTYXEa5GPAPQItz7ksz7LYLuCc22uV6oN851z7DvpJEzjk+u+sQxblZPHRPIzmZGRe35WaF+Ydf2kF2Zgaf/+7hKxxFRIJoLnfobwM+BtxsZvtiH3eY2QNm9kBsnyeAN4BjwN8DDyanXJnND1o62X+mj99492pKFme9afuSghwevOk6dh87z/PHz3tQoYgky6ydos653YDNso8DPpWoomR+IhHHF79/lBVli/m5bTUz7nf39XV85Scn+NPvHeWxT5YSfRMmIkGnJ0XTyPcOdXCkY5Bfv3U14YyZX9qczAw+fcsq9p7u4yev6y5dJF0o0NPIo6+coaowh/dtqpp1349sr6VkcRbfeuXMrPuKSDAo0NNEe/8IP3m9iw9vryEjNHsTSlY4xAe3VPHU4XP0Do2noEIRSTYFepp4fO9ZnIMPb5/7Q0Mf2V7L+FSEf9unEaYi6UCBngacc/xL0xl2rihhWWnunL+voaqADdUF/Mue1iRWJyKpokBPA3tP93Kye5iPNF79I/0f2V7LobYBjnTM6+FfEfERBXoa+P7hc4RDxnvWL7nq771941IAnjqk6QBEgk6BngZ+cPgc168spSAn86q/tyI/hy21RfygRYEuEnQK9IB7o+sCx7uGuLW+Yt7HeHfDEva39nNuYDSBlYlIqinQA+6HLZ0A3FJ/9c0tcbfGvjd+LBEJJgV6wD3Vco51S/OpLZn76Jbp1izJo7ZkkZpdRAJOgR5gfcPjNJ3suXiHPV9mxq31S9h97Dwj41rVSCSoFOgB9uIb3UQc3LT22hcLeeeacsYnIzSd6klAZSLiBQV6gD1/vJvcrAw21RRd87F2LC8hHDKeO9adgMpExAsK9AB77th53rKihKzwtb+Mi7PDbKkt4gXNkS4SWAr0gDo3MMrxriFuuK40Yce8YVUZB8720z8ykbBjikjqKNAD6oXj0aaRG64rS9gxb7iulIiDl0+oHV0kiBToAfXcsfMULsqkobIgYcfcuqyI7HCI546p2UUkiBToAfX88W7eurKU0BzmPp+r7HAGO5aXaK1RkYBSoAdQW98IZ/tG2LmyJOHH3rmihNfOXaBvWIteiASNAj2Amk71AtBYl/hA3768GIhOySsiwaJAD6C9p3pZlJlBfWV+wo+9pbaIjJCx55QCXSRoFOgB1HSqhy21RYQzEv/y5WaFaagsoOmkAl0kaBToATM0NklL+yCNsaaRZNheV8z+1j4mpiJJO4eIJJ4CPWD2n+ljKuLYXpe8QG9cXszoRITDbVqWTiRIFOgBs+dUL2awdVly79Dj5xKR4FCgB0zTqV7WVORTuOjql5ubq8rCRVQXLVKgiwSMAj1AnHPsO9PH1mXXPrvibLYuK+JVDV0UCRQFeoCc7B6mf2SCLbXJD/QttUW09Y/SOah1RkWCQoEeIM2tfQAJmf98NptjvzSaz/Qn/VwikhgK9ADZd6aPnMwQa5bkJf1c66sKCNlPf4mIiP/NGuhm9rCZdZrZwRm2F5rZ/zOz/WZ2yMzuTXyZAtDc2s+GqsKkPFA0XW5WmDVL8tnXqjt0kaCYSzI8Atx2he2fAg475zYDNwFfNLOsay9NLjUxFeHg2f6LTSGpsLmmiObWPpxzKTuniMzfrIHunHsWuNKKBw7INzMD8mL7TiamPIl77dwgY5MRNtUUpuycm2uL6Bue4HTPcMrOKSLzl4j37n8F1ANtwAHg15xzl31m3MzuN7MmM2vq6upKwKkXjv2xzslUjHCJ21wb/eWxX80uIoGQiEB/L7APqAK2AH9lZpddRsc595BzrtE511heXp6AUy8cza19FOVmsqwkN2XnXLMkn+xwiP1n1DEqEgSJCPR7gcdd1DHgBLAuAceVSzS39rOxupBoy1ZqZGaEaKgq4MBZ3aGLBEEiAv00cAuAmS0B1gJvJOC4EjM2OcVr5wbZUJ269vO49VUFtLQNEImoY1TE7+YybPFR4AVgrZm1mtl9ZvaAmT0Q2+XzwA1mdgD4IfAZ55wWpUyg1zouMBlxbKjyItALGRyb5EyvOkZF/C482w7Oubtm2d4GvCdhFcmbHGyLNnlsqL5s10RSra+KnvNQ2wB1pYtTfn4RmTs9KRoAB8/2k58TTmmHaNyaJflkhIxDbWpHF/E7BXoAHGwbYENVajtE43IyM1hdkcchLXYh4nsKdJ+bmIrQ0j7gSXNLXENVgQJdJAAU6D53vOsC45MRT0a4xK2vKqRrcExT6Yr4nALd5w6ejd4Zr/dghEvcpR2jIuJfCnSfO3i2n9ysDFaUeTfCpCEW6Fo0WsTfFOg+19I+wLql0ZEmXinIiU45oJEuIv6mQPcx51w00Cu96xCNW6+OURHfU6D7WFv/KAOjk9T7JNBPdQ8zMDrhdSkiMgMFuo+1xO6IGyrzPa7kp52yLbpLF/EtBbqPtbRHw3PtUu/v0Bs00kXE9xToPtbSMUBdaS552bNOuZN0FfnZlOVlKdBFfEyB7mNH2gep98HdOYCZ0VBVqJEuIj6mQPep4fFJTnQPsc4H7edx66sKONZ5gbHJKa9LEZHLUKD71NGOQZzDFyNc4tZXFTAZcbzWccHrUkTkMhToPtXSPghAg68CPTrSRc0uIv6kQPeplvYB8rPD1BQv8rqUi+pKclmclXFx9I2I+IsC3aeiT4jmezIH+kxCIWPt0vyL7x5ExF8U6D7knONIx6Cv2s/j6isLaOkYwDktGi3iNwp0H2rtHeHCmD8e+Z+uoaqAwdFJWntHvC5FRKZRoPvQ4Vgb9bql/hmyGBf/JXNY7egivqNA96GW9gHMYK0PA33d0nzMUMeoiA8p0H2opX2AFaWLyc3y/pH/6XKzwiwvXaxAF/EhBboPtbT7s0M0rr5SI11E/EiB7jODoxOc7hmm3keP/E9Xv7SA0z3DDGpudBFfUaD7zGvnone+fr5Dj0+le7RDd+kifqJA95nDsaYMPyw7NxONdBHxJwW6z7S0D1CQE6aqMMfrUmZUWZhD4aJMdYyK+IwC3Wda2georyzw1SP/05kZ9ZX5F99NiIg/KNB9JBJxHPXpI//T1VcWcLRjgKmIpgAQ8YtZA93MHjazTjM7eIV9bjKzfWZ2yMx+nNgSF45TPcMMj0/5asrcmTRUFjA6EeFk95DXpYhIzFzu0B8Bbptpo5kVAX8DfMA5tx74SGJKW3iOxNqkg3KHDnBYa4yK+Masge6cexboucIu/wV43Dl3OrZ/Z4JqW3Ba2gcIGaxekud1KbNavSSPcMjUMSriI4loQ18DFJvZM2a2x8zumWlHM7vfzJrMrKmrqysBp04vh9sHWVmeR05mhtelzCo7nMF15XkKdBEfSUSgh4HtwM8A7wX+wMzWXG5H59xDzrlG51xjeXl5Ak6dXuIjXIJCUwCI+EsiAr0V+J5zbsg5dx54FticgOMuKP0jE5ztG/H1I//TNVQV0DEwSu/QuNeliAiJCfR/A240s7CZ5QI7gZYEHHdBCVKHaFy8VjW7iPjDrPOzmtmjwE1AmZm1Ap8DMgGcc192zrWY2X8AzUAE+IpzbsYhjnJ5R2LzogRhyGLcpVMA3LCqzONqRGTWQHfO3TWHfb4AfCEhFS1QLe0DlCzOoiI/2+tS5qwsL5vy/GzN6SLiE3pS1Cda2gdiqwH595H/y6mvLFDHqIhPKNB9YCriOHouGI/8T1dfmc+xzkHGJyNelyKy4CnQfeDE+SFGJyKBDPSGygImphzHuy54XYrIgqdA94GWiyNcgjNkMa5BUwCI+IYC3QeOdAwQDhmrKvz/yP90K8oWkxUOaeiiiA8o0H2gpX2QVRV5ZIf9/8j/dOGMEGuX5NPSoUAX8ZoC3QfiI1yCKj4FgHOaG13ESwp0j/UOjdPePxrIDtG4hsoCeobG6Rwc87oUkQVNge6xeNvz+qpCjyuZPy0aLeIPCnSPHQ7wCJe4dRrpIuILCnSPHWobYGlBDqV5wXnkf7rCRZlUFy3SSBcRjynQPXa4bYCGquC2n8dFpwBQoIt4SYHuodGJKY51XWB9GgR6Q1VB7InXKa9LEVmwFOgeev3cBaYiLlBT5s6koTKfiIOjHZqoS8QrCnQPHWrrB0ibJhfQSBcRLynQPXS4fYC87DC1xblel3LNaotzWZyVoXZ0EQ8p0D10uG2A+sp8QqFgzYF+OaGQsU4doyKeUqB7JBJxtLQPBPqBoukaKgs4oikARDyjQPfIqZ5hhsan0qJDNK6+soDBsUlae0e8LkVkQVKgeyT+VGU6dIjGxZ92VceoiDcU6B453N4f2DnQZ7J2aT5mmgJAxCsKdI8cbhtgVUUeOZnBmwN9JrlZYVaULlbHqIhHFOgeOZQmj/xPV19ZoMUuRDyiQPdA1+AYnYNjadUhGtdQVcCZnhEGRye8LkVkwVGgeyDeJJGed+jRjtEjmgJAJOUU6B6IjwJJxzv0es2NLuIZBboHDrUNUF20iKLcLK9LSbilBTkU5WaqY1TEAwp0Dxw625+WzS0AZsb6qgIO6Q5dJOUU6CnWPzLBG+eH2FyTPo/8T7exuogjHQOMTWpudJFUUqCn2KGz0SlzN9YUeVxJ8myqKWRiynGkXR2jIqmkQE+x/a3RQN9Unb536Jti7z6aW/s8rkRkYZk10M3sYTPrNLODs+y3w8wmzezDiSsv/TS39rGsJJfixenXIRpXXbSI0sVZF395iUhqzOUO/RHgtivtYGYZwJ8A309ATWmtubX/4h1sujIzNtYUckCBLpJSswa6c+5ZoGeW3T4NPAZ0JqKodNV9YYyzfSNpH+gAm2qKeL1zkOHxSa9LEVkwrrkN3cyqgZ8F/nYO+95vZk1m1tTV1XWtpw6c5nj7eRp3iMZtrikk4uDgWQ1fFEmVRHSK/h/gM865yGw7Oucecs41Oucay8vLE3DqYGlu7ccMNqRxh2jcRnWMiqRcOAHHaAS+aWYAZcAdZjbpnPtOAo6dVppb+7iuPI+87ET82P2tIj+HysKci+9KRCT5rjlZnHMr4p+b2SPAdxXmb+acY39rP+9YU+Z1KSmzqaZQd+giKTSXYYuPAi8Aa82s1czuM7MHzOyB5JeXPtr7Rzl/YYzNC6D9PG5TTREnu4fpH9ZUuiKpMOsdunPurrkezDn3y9dUTRr7aYdo+refx8Wv9cDZfm5cvXDemYh4RU+Kpkhzax/hkF2cXnYh2FQdfTeyX80uIimhQE+R5tZ+1i7NT6s1RGdTmJtJXWmu2tFFUkSBngLOOZpb+xbE+PPpNtUU6YlRkRRRoKfAqe5hBkYn03rK3JlsrimkrX+UrsExr0sRSXsK9BSItyFvXICBvrFaDxiJpIoCPQX2n+knOxxizZJ8r0tJuQ3VhYQM9p1RoIskmwI9Bfac7mVzbRGZGQvvx704O8y6pQXsPd3rdSkiaW/hJUyKjU5McehsP9vrir0uxTPb64rZd7qPyalZp/sRkWugQE+y5tZ+JiOO7csWdqAPjU9x9JyWpBNJJgV6ku05FW1q2LbA79AB9p5Ss4tIMinQk2zPqV5Wli2mJI2XnJtNTfEiyvOzL/5yE5HkUKAnkXOOvad7F/TdOUSXpNu+rJg96hgVSSoFehKd7B6mZ2h8QXeIxm2vK+ZMzwidA6NelyKSthToSdR0MroUqwL9p30IanYRSR4FehK9fKKHotxMVpXneV2K5zZWF5KTGeKlE7OtNy4i86VAT6KXTvTwluUlhELmdSmeywqH2LasWIEukkQK9CRp7x/hdM8wO1eWel2Kb+xcUcqRjgGtYCSSJAr0JHnpjeid6M4VJR5X4h87V5bgHLxyUnfpIsmgQE+Sl070kJ8TXlArFM1mS20RWRkhXjrR7XUpImlJgZ4kL53oZsfyEjLUfn5RTmYGW2qLeFnt6CJJoUBPgs7BUd7oGlJzy2XsXFnCwbYBLoxNel2KSNpRoCfBi/H2c3WIvsn1K0uZijheVrOLSMIp0JNg9+tdFOSEL67WIz+1va6Y7HCI3a8r0EUSTYGeYM45dr9+nhuuK1P7+WXkZGawY3kJzx0773UpImlHgZ5gJ7uHaesf5cbVZV6X4ls3ri7j6LlBzesikmAK9ATb/XoXADeuUqDPJP6z2a27dJGEUqAn2O5j56kpXkRdaa7XpfhWQ2UBxbmZCnSRBFOgJ9DkVITnj3dz46oyzNR+PpNQyLhhVRm7Xz+Pc87rckTShgI9gfad6WNwdFLt53PwjtVldA6OcaRD64yKJIoCPYGePtJJOGS8fXW516X43rvWVgDRn5mIJMasgW5mD5tZp5kdnGH7L5pZs5kdMLPnzWxz4ssMhqePdLJjeQmFizK9LsX3Kgpy2FhdqEAXSaC53KE/Atx2he0ngHc65zYCnwceSkBdgdPaO8yRjkFuXlfhdSmBcfO6Cvae7qVnaNzrUkTSwqyB7px7FphxNiXn3PPOufi6Yi8CNQmqLVB+FLvTvLlegT5Xt9RX4Bw8c1R36SKJkOg29PuAJ2faaGb3m1mTmTV1dXUl+NTeevpIJ3WluawsW+x1KYGxoaqQsrxsNbuIJEjCAt3M3kU00D8z0z7OuYecc43Oucby8vTpOBwam+T5493cvK5CwxWvQihk3LyunB8f7WJscsrrckQCLyGBbmabgK8AH3TOLbhZl54+0snYZITb1i/1upTAuW3DUgbHJnn+2IL7ZyOScNcc6Ga2DHgc+Jhz7rVrLyl4njzYTlleNo3LNf/51XrbqjLys8M8caDd61JEAm8uwxYfBV4A1ppZq5ndZ2YPmNkDsV0+C5QCf2Nm+8ysKYn1+s7w+CRPH+nk9g1LNbviPGSHM3h3wxK+f/gcE1MRr8sRCbTwbDs45+6aZfvHgY8nrKKAeeZoF6MTEW7fqOaW+bp9YyWPv3qW549388416dO3IpJqelL0Gj1xoJ3SxVm8Rc0t8/b21WXkZYd5olnNLiLXQoF+DQZHJ/hByzlu27CUcIZ+lPOVk5nBrfUVPHmwndEJjXYRmS+l0DV48mAHoxMR7ty2IJ+lSqg7t9UwMDqpMeki10CBfg0e29PKirLFbFtW5HUpgfe2VWUsKcjmsT2tXpciElgK9Hk60zPMSyd6uHNrtR4mSoCMkPGhrdU881oXXYNjXpcjEkgK9Hl6fO9ZAH52W7XHlaSPD2+rYSri+Ld9Z70uRSSQFOjzMBVxfLvpDDdcV0pNsZaaS5TVS/LZXFvEN185o5WMROZBgT4PP2w5x9m+Ee55a53XpaSdj11fx7HOC7xwXFMBiFwtBfo8fO2FU1QW5nBr/RKvS0k779tUSXFuJv/4wkmvSxEJHAX6VTrWOcjuY+e5+/o6jT1PgpzMDH5+xzKeOhx9FyQic6dEukqPPH+SrIwQP7+j1utS0tYv7lwGwD+9cMrjSkSCRYF+FToHRvl2Uyt3bqumLC/b63LSVm1JLndsrOSfXzxF//CE1+WIBIYC/Sr8/U/eYHIqwidvus7rUtLep961igtjkzzy/EmvSxEJDAX6HPUOjfP1l07zgc1V1JVqmblkq68s4Nb6JTz83AkujE16XY5IICjQ5+jvnn2D4fEpHnzXKq9LWTB+9eZV9I9M8NXdJ7wuRSQQFOhzcKZnmIefO8GdW6tZsyTf63IWjC21RbynYQlf/vFxOgdHvS5HxPcU6HPwhe8dJWTwX9+71utSFpzfvaOesckIf/bU616XIuJ7CvRZ7DnVy679bXzi7SupKlrkdTkLzoqyxXzsrXV865XTHG4b8LocEV9ToF/B2OQUn3msmarCHH7lnRrZ4pVfu2U1JYuz+MxjzUxq3VGRGSnQr+Cvnz7Gsc4L/PGdG8nLnnX5VUmSotws/scHNnDgbD//oA5SkRkp0Gew/0wff/PMce7cVs1Nayu8LmfBu2PjUt67fglfeuo1jnYMel2OiC8p0C+jd2icB7++lyUFOXz2fQ1elyOAmfFHH9pIwaJMPvnPexgc1ROkItMp0KeZijh+49v76Boc42/v3kZRbpbXJUlMeX42f3XXVk71DPPb/9pMJKI500UupUC/hHOO3//OQZ452sVn39/AphqtFeo3O1eW8pnb1vLkwQ7++IkWLYQhcgn19MU45/jT7x/l0ZdP8+BN13H39Vq8wq8+8faVnO0d4Su7T1CSl8WDN+npXRFQoAMQiTg+/++H+epzJ/mFHbX8Nz1A5Gtmxufev57e4Qn+938cZWR8it989xot1i0L3oIP9KGxSX77sWb+vbmde9+2nD/4mQYFQwCEQsaXPrqZ3KwM/vLpY3T0j/L5D20gJzPD69JEPLOgA72lfYBf/cZeTpwf4nduX8evvGOlwjxAwhkh/uedG6koyOEvfvg6B9sG+Mu7trKqIs/r0kQ8sSA7RYfHJ/lfTx7h/X+5m/6RSf754zt54J3XKcwDyMz4zXev4av37qCjf4Q7/vwn/NlTrzE6MeV1aSIpZ16NEmhsbHRNTU0pPeeFsUm++fJpvvzj45y/MM5HG2v43dvrKV6soYnpoHNwlD/6bgu79rexpCCbB29axUcba1mUpWYYSR9mtsc513jZbeke6FMRR9PJHnbtb+M7r55laHyKG1eV8RvvXsP2uuKkn19S78U3uvnS91/j5ZM95OeE+bltNbx/cyVba4sJhfQuTILtmgLdzB4G3gd0Ouc2XGa7AX8O3AEMA7/snNs7W1HJCvTxyQiH2vrZc6qXppO9vHKyh+6hcbLDId63qYq7r1/G1mUK8nTnnOOVk718/aVTPHmgg/GpCOX52exYXsz2uhIa64ppqCogM2NBtjpKgF1roL8DuAB8bYZAvwP4NNFA3wn8uXNu52xFzTfQ+4cneL1zkO6hcXpiHx39o5zsHuJU9zBn+0aYij1BuKwkl8a6Ym6ur+BdaytYrAm2FqSB0QmebunkR0c7aTrZy9m+EQAyQkZN8SKWly6mrjSXJQU5lC7OomRxFqV5WeRmhcnNymBRVgaLMjPIzQqToTt88diVAn3WhHPOPWtmy6+wyweJhr0DXjSzIjOrdM61z6vaWTz7eheffvTV//S1wkWZLC/NZUttER/aUsW6ygIa64qpKMhJRgkSMAU5mXxoazUf2loNQHv/CE0neznaMXjxRmDv6V4GR2dfuzQcMjLiH2ZkZET/DIWMcMgIWXTb9P71y/0amN4J/6Z9EnEM8aWf31HLx9++MuHHTcQtazVw5pK/t8a+9qZAN7P7gfsBli1bNq+T7VxZwiP37qAsL5uS2N2Uxh7L1agsXMT7Ny/i/Zv/89dHJ6boHhqn+8IYPUPjDI9PMTI+xfDEFCPjk4yMRxibnGLKOaamHFPOEYk4JiOOiHNMxT+fNsfM5d4DT39jPH2f6e+cL/s++k3H0DQIQVGWl52U46a0DcI59xDwEESbXOZzjIr8HCrW6s5bEi8nM4PqokVUa2UqCahE9AidBWov+XtN7GsiIpJCiQj0XcA9FnU90J+s9nMREZnZrE0uZvYocBNQZmatwOeATADn3JeBJ4iOcDlGdNjivckqVkREZjaXUS53zbLdAZ9KWEUiIjIveqpCRCRNKNBFRNKEAl1EJE0o0EVE0oRnsy2aWRdwap7fXgacT2A5QbEQr1vXvDAsxGuG+V13nXOu/HIbPAv0a2FmTTNNTpPOFuJ165oXhoV4zZD461aTi4hImlCgi4ikiaAG+kNeF+CRhXjduuaFYSFeMyT4ugPZhi4iIm8W1Dt0ERGZRoEuIpImfB3oZnabmR01s2Nm9juX2Z5tZt+KbX9plqXyAmEO1/ybZnbYzJrN7IdmVudFnYk223Vfst/PmZkzs8APcZvLNZvZR2Ov9yEz+0aqa0y0Ofz7XmZmPzKzV2P/xu/wos5EMrOHzazTzA7OsN3M7C9iP5NmM9s275M553z5AWQAx4GVQBawH2iYts+DwJdjn/8C8C2v607BNb8LyI19/smgX/Ncrzu2Xz7wLPAi0Oh13Sl4rVcDrwLFsb9XeF13Cq75IeCTsc8bgJNe152A634HsA04OMP2O4AniS4Jez3w0nzP5ec79LcAx5xzbzjnxoFvEl2Q+lIfBP4x9vm/ArfY9FVzg2XWa3bO/cg5Nxz764tEV4gKurm81gCfB/4EGE1lcUkyl2v+BPDXzrleAOdcZ4prTLS5XLMDCmKfFwJtKawvKZxzzwI9V9jlg8DXXNSLQJGZVc7nXH4O9JkWn77sPs65SaAfKE1Jdckxl2u+1Fv9J6QAAAIESURBVH1Ef7MH3azXHXsbWuuc+/dUFpZEc3mt1wBrzOw5M3vRzG5LWXXJMZdr/u/A3bHFdJ4APp2a0jx1tf/vZ5TSRaIlcczsbqAReKfXtSSbmYWALwG/7HEpqRYm2uxyE9F3Ys+a2UbnXJ+nVSXXXcAjzrkvmtlbgX8ysw3OuYjXhQWBn+/Q57L49MV9zCxM9C1ad0qqS445LbhtZrcCvwd8wDk3lqLakmm2684HNgDPmNlJou2MuwLeMTqX17oV2OWcm3DOnQBeIxrwQTWXa74P+DaAc+4FIIfoBFbpbE7/7+fCz4H+CrDazFaYWRbRTs9d0/bZBfxS7PMPA0+7WC9DQM16zWa2Ffg7omEe9DbVuCtet3Ou3zlX5pxb7pxbTrTv4APOuSZvyk2Iufz7/g7Ru3PMrIxoE8wbqSwyweZyzaeBWwDMrJ5ooHeltMrU2wXcExvtcj3Q75xrn9eRvO4BnqV3+A6idyXHgd+Lfe0Pif5nhuiL/S9EF6h+GVjpdc0puOYfAOeAfbGPXV7XnIrrnrbvMwR8lMscX2sj2tR0GDgA/ILXNafgmhuA54iOgNkHvMfrmhNwzY8C7cAE0Xdd9wEPAA9c8jr/dexncuBa/m3r0X8RkTTh5yYXERG5Cgp0EZE0oUAXEUkTCnQRkTShQBcRSRMKdBGRNKFAFxFJE/8fP58+feyTWKEAAAAASUVORK5CYII=\n",
"text/plain": [
"