"
]
},
{
"cell_type": "markdown",
"metadata": {
"lines_to_next_cell": 0,
"nbgrader": {
"grade": false,
"grade_id": "cell-0a4af3607542b85a",
"locked": true,
"schema_version": 1,
"solution": false
}
},
"source": [
"# Introduction"
]
},
{
"cell_type": "markdown",
"metadata": {
"lines_to_next_cell": 0,
"nbgrader": {
"grade": false,
"grade_id": "cell-6f2c0b04c4f4fdf6",
"locked": true,
"schema_version": 1,
"solution": false
}
},
"source": [
"Let us recall that $J_{EC}$ and $J_{WC}$ represent the flux of mass [M/T] towards the center. If these are positive, the mass in the center will increase. If there is a source of mass in the center, the mass will increase as well. Therefore, the steady-state equation writes:\n",
"\n",
"\\begin{equation}\n",
"J_{EC}+J_{WC} + Q = 0\n",
"\\end{equation}\n",
"\n",
"With Q being written in terms of mg/s (Q: [M/T]) . Therefore mass balance equation with a source term writes:\n",
"\n",
"\\begin{equation}\n",
"\\left( -2 \\frac{D\\theta}{\\Delta x} c_{\\left[i\\right]} + \\frac{D\\theta}{\\Delta x}c_{\\left[i-1\\right]} + \\frac{D\\theta}{\\Delta x} c_{\\left[i+1\\right]} \\right) \\Delta y \\Delta z + q \\Delta x \\Delta y \\Delta z = 0\n",
"\\end{equation}\n",
"\n",
"with $q$ being the volumetric source ([M/T/L$^3$])\n",
"\n",
"We can rearrange the latter equation as:\n",
"\n",
"\\begin{equation}\n",
"\\left( 2\\frac{D\\theta}{\\Delta x} c_{\\left[i\\right]} - \\frac{D\\theta}{\\Delta x}c_{\\left[i-1\\right]} - \\frac{D\\theta}{\\Delta x} c_{\\left[i+1\\right]} \\right) \\Delta y \\Delta z = q \\Delta x \\Delta y \\Delta z\n",
"\\end{equation}\n",
"\n",
"We can also divide everything by $\\Delta x \\Delta y \\Delta z$:\n",
"\n",
"\\begin{equation}\n",
"\\frac{1}{\\Delta x} \\left( 2\\frac{D\\theta}{\\Delta x} c_{\\left[i\\right]} - \\frac{D\\theta}{\\Delta x}c_{\\left[i-1\\right]} - \\frac{D\\theta}{\\Delta x} c_{\\left[i+1\\right]} \\right) = q\n",
"\\end{equation}\n",
"\n",
"which gives\n",
"\\begin{equation}\n",
"2\\frac{D\\theta}{\\Delta x^2} c_{\\left[i\\right]} - \\frac{D\\theta}{\\Delta x^2}c_{\\left[i-1\\right]} - \\frac{D\\theta}{\\Delta x^2} c_{\\left[i+1\\right]} = q\n",
"\\end{equation}"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-98a865dfd995366b",
"locked": true,
"schema_version": 1,
"solution": false
}
},
"outputs": [],
"source": [
"import matplotlib.animation as animation\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from IPython.display import display\n",
"from IPython.display import HTML\n",
"from numpy.testing import assert_allclose\n",
"\n",
"# the next two imports are to show animations and video!"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-bd0b49f6ae1acef3",
"locked": true,
"schema_version": 1,
"solution": false
}
},
"outputs": [],
"source": [
"# this function deals with harmonic averaging when diffusion is not the same everywhere.\n",
"# It doesn't change anything when diffusion is homogeneous but \n",
"# you can try to see how it affects the behavior.\n",
"# To modify a read-only cell, create a new cell below and copy these lines\n",
"# into that cell\n",
"\n",
"\n",
"def avg(Di, Dj):\n",
" \"\"\"\n",
" Computes the harmonic average between two values Di and Dj\n",
" Returns 0 if either of them is zero\n",
" \"\"\"\n",
" if (Di * Dj) == 0:\n",
" return 0\n",
" else:\n",
" return 2 / (1 / Di + 1 / Dj)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-54d67b8f8547140a",
"locked": true,
"schema_version": 1,
"solution": false
}
},
"outputs": [],
"source": [
"def Build_1D_Inhomo_Matrix_Source(bc_left, bc_right, n, D, Width, poro, Q):\n",
" \"\"\"\n",
" Constructs a coefficient matrix and an array with varying diffusion coefficient and a source term\n",
"\n",
" Parameters\n",
" ----------\n",
" bc_left: (float)\n",
" left boundary condition (mg/L)\n",
" bc_right: (float)\n",
" right boundary condition (mg/L)\n",
" n: (int)\n",
" number of cells\n",
" D: (float vector)\n",
" values of the diffusion coefficient ((dm)^2/day)\n",
" Width: (float)\n",
" Total phyiscal width of the domain (dm)\n",
" poro: (float)\n",
" porosity value\n",
" Q: (float vector)\n",
" volumetric source term (mg/L/day)\n",
"\n",
" Returns\n",
" -------\n",
"\n",
" A: nxn float array\n",
" b: vector of length n\n",
"\n",
" that solve the\n",
" discretized 1D diffusion problem Ax = b\n",
" \"\"\"\n",
" Matrix = np.zeros((n, n))\n",
" RHS = np.zeros(n)\n",
" dx = Width / (n - 1)\n",
" coef = poro / dx / dx\n",
" for i in range(n):\n",
" if i == 0:\n",
" RHS[i] = bc_left\n",
" Matrix[i, i] = 1\n",
" elif i == n - 1:\n",
" RHS[i] = bc_right\n",
" Matrix[i, i] = 1\n",
" else:\n",
" RHS[i] = Q[i]\n",
" East = coef * avg(D[i], D[i + 1])\n",
" West = coef * avg(D[i], D[i - 1])\n",
" Matrix[i, i] = East + West\n",
" Matrix[i, i + 1] = -East\n",
" Matrix[i, i - 1] = -West\n",
" return Matrix, RHS"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Initial conditions"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-b41a29613812db59",
"locked": true,
"schema_version": 1,
"solution": false
}
},
"source": [
"We will use similar conditions than in the previous assignment. However, it is a good practise to use units which are representative of the current problem, otherwise we have to deal with very large or very small numbers.\n",
"\n",
"The diffusion coefficient of solutes in pure water is usually 2$\\times$10$^{-9}$ m$^2$/s. Concentrations are usually expressed in mg/L, corresponding to mg/dm$^3$. To avoid future unit problems, let us define\n",
"\n",
"- one day as the time unit\n",
"- one dm as the length unit (leading to Liters being the adequate volume unit)\n",
"- mg as the mass unit\n",
"\n",
"In these units, the diffusion is expressed in dm$^2$/day, the width in dm, and the rate in mg/L/day.\n",
"We will here change the values of the parameters used previously.."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"lines_to_next_cell": 2,
"nbgrader": {
"grade": false,
"grade_id": "cell-b9bd7eaf90619a64",
"locked": true,
"schema_version": 1,
"solution": false
}
},
"outputs": [],
"source": [
"c_left = 1 # mg/L\n",
"c_right = 0 # mg/L\n",
"n = 51\n",
"Diff = 2e-9 * 100 * 24 * 3600 # dmĀ²/day\n",
"D = Diff * np.ones(n)\n",
"Q = np.zeros(n)\n",
"Q0 = 5e-2 # mg/L/day\n",
"Q[int(n / 4) : int(n / 2)] = Q0 # mg/L/day\n",
"Width = 2 # dm\n",
"poro = 0.4\n",
"nTstp = 201\n",
"dt = 0.5 # days\n",
"\n",
"c_init = np.zeros((nTstp, n))\n",
"c_init[:, 0] = 1 # Boundary condition"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# If you want to change the parameters value, you are welcome to do so here.\n",
"# But know that your results will be tested based on the parameters given above\n",
"# We encourage you to try different things\n",
"# For example, you change the diffusion at certain places, change the source term, ...\n",
"# Trying out these things will help you build physical intuition\n",
"# So we leave this cell for you to change these parameters.\n",
"# As before, make a writeable copy of that cell to make changes\n",
"#\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-e983be851ed2d326",
"locked": true,
"schema_version": 1,
"solution": false
}
},
"source": [
"## Save these into a dict for later use"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"lines_to_next_cell": 2,
"nbgrader": {
"grade": false,
"grade_id": "cell-b19ecc741595f838",
"locked": true,
"schema_version": 1,
"solution": false
}
},
"outputs": [],
"source": [
"init_dict = dict(\n",
" c_left=c_left,\n",
" c_right=c_right,\n",
" c_init=c_init,\n",
" n=n,\n",
" Diff=Diff,\n",
" D=D,\n",
" Q=Q,\n",
" Width=Width,\n",
" poro=poro,\n",
" nTstp=nTstp,\n",
" dt=dt,\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-8ce4321aef65caa0",
"locked": true,
"schema_version": 1,
"solution": false
}
},
"source": [
"# Steady state"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-1a17c575606df6ff",
"locked": true,
"schema_version": 1,
"solution": false
}
},
"outputs": [
{
"data": {
"text/plain": [
"Text(0, 0.5, 'Concentration (mg/L)')"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd8VGXa//HPlQQyAVIEAoQaqnQIJKjgo2JZAQviimCBhOACyu6i/nTXsuuqj7u66u6KZRXEJIAorr2sHQsWkITeBQEBE0goqaTn/v0xhzwxpkzKmZNMrvfrNS9mTv3mZMg1c8597luMMSillFIAfk4HUEop1XRoUVBKKVVOi4JSSqlyWhSUUkqV06KglFKqnBYFpZRS5bQoKKWUKqdFQSmlVDktCkoppcoFOB2grjp27GgiIyOdjqGUUs3K+vXrjxljwmtbrtkVhcjISFJSUpyOoZRSzYqI/OjJcnr6SCmlVDktCkoppcppUVBKKVVOi4JSSqlyWhSUUkqV06KglFKqnBYFpZRS5ZrdfQrKt+UVlnAku4Cj1uNIViFlxhDWphVntGlNWFArQq3nHdsF0jpAP9co1Zi0KCjHlJSWsflwJl/uzuDLPcfYl55LTmGJx+sHBvgxqucZnN2nA2f3ac+IHmG4WvnbmFgp36dFQXlVenYBX+zO4MvvM/hqTwbZBSX4CUT1PINfj+5O5xAXXUID6RzsonOoi84hLgL8hKz8Yk6eKiLzVDGZp4o4eaqYPUdzWbvvOE+s+h7zKbQO8GNUzzDOH9CJqdHd6dgu0OkfV6lmR4wxTmeok+joaKPdXDQ/O1KzWbz6B97bkkZJmaFzSCDnDwjn/AGdOLdfR0LbtKr3trNOFbPuwAm+23ectfuPs+2nbFr7+3H5iAhiz4lkRI+wRvxJlGqeRGS9MSa61uW0KCi7GGP4Zu9xFq3+ga/2HKNNa3+mx/RkanR3BnYJRkRs2e/e9FyWrznAa+sPk1dUysgeYcSNjWTSsAi9BqFaLC0KyjHGGN7bksazX/zAjrRswoMDiRsbyY1n9WrQN4K6yiko5vX1h1m25kf2HcujS4iLey4bxBXDI2wrSEo1VVoUlCP2H8vjnje2smbfcfqGt2XOeX24KqobgQHOXQAuKzN8tfcYj3+0m60/ZXFOnw48OHkI/TsHO5ZJKW/ToqC8qri0jMWr97Fw1R4C/f24a9JArovpiZ9f0/lEXlpmeHndQR77aDd5hSXEn9ub31/Un3aB2t5C+T4tCsprNh48yV2vb2X30RwmDevCX64YQucQl9OxqnUir4hHP9zFyuRDdA4J5M+XD+by4V2djqWUrbQoKNsVFJfyyAe7WLrmAJ2DXfzvVUO5ZHBnp2N5bMPBk9z39ja2/ZTNtOgePDB5iN7noHyWp0XBtqYYItJDRD4XkZ0isl1EFlSxjIjIkyKyV0S2iMgou/KoxnUkq4Bpi9eS9O0BZp7di09uP69ZFQSAUT3P4O355/Lb8f14JeUQ1zz3LYdOnHI6llKOsrN9Xgnw/4wxg4CzgfkiMrjSMhOB/tZjDvCsjXlUI1n/40muePpr9h7NYdGM0TwweSjBLu+1KmpM/n7CHZeeyZKZ0fx4/BSXP/U1n+9OdzqWUo6xrSgYY9KMMRus5znATqBbpcUmA8uM21ogTEQi7MqkGu6V5INct3gtbVr78+b8cVw6pIvTkRrFxYM78+5vzyUi1EV8UjL/+uR7ysqa16lVpRqDV+7kEZFIIAr4rtKsbsChCq8P88vCgYjMEZEUEUnJyMiwK6aqQXFpGX95ext/fH0rZ/Vpz9vzxzHAx5p0RnZsy5u3jGNKVDcWrtpD/NJkcuvQF5NSvsD2oiAi7YDXgVuNMdmVZ1exyi8+nhljFhtjoo0x0eHh4XbEVDXIyi9mxgvfsXTNj9x0bm8S42IIa9Pa6Vi2CGrtzz+mjuChq4by9Z5j3LDkOzJPFTkdSymvsbUoiEgr3AVhhTHmjSoWOQz0qPC6O5BqZyZVN1mn3AVh/Y8n+ee1I/jT5YMJ8PftriJEhBvP7sWzN45mZ2o20xev5VhuodOxlPIKO1sfCfACsNMY889qFnsHmGm1QjobyDLGpNmVSdXNybwirl+yll1pOTx342iuHtXd6UhedcngzrwQF82B43lcu2gNaVn5TkdSynZ2fuQbB8wALhSRTdZjkojME5F51jLvA/uAvcDzwC025lF1cCKviOuXfMee9FwWzRzNRYOaV3PTxvI//cNZFn8W6dmFTH1uDQePa5NV5dv05jX1C8dzC7lhyXfsP5bH8zOjOW+AXsfZcjiTmQnrcAX48+JNZ9GvUzunIylVJ47fvKaap4ycQq57fi0HjufxQmyMFgTL8O5hrJxzNiVlZUxbtIa96TlOR1LKFloUVLljue6CcOhEPglxMZzbv6PTkZqUgV1C+M/ccxARYhOSSc8ucDqSUo1Oi4IC3P0Y/WZZCodPniJxVgxj+2pBqEqf8HYkxsVw8lQRcYl6H4PyPVoUFGVlhtv/s4lNhzJ5YloUZ/fp4HSkJm1Y91CeuWEUu4/mcPOL6ykuLXM6klKNRouC4tGPdvP+1iPcM3EQE4b6RrcVdht/ZicenjKMr/Yc467Xt9LcGmwoVR0dXaSFe3ndQZ778gduOKsnN/1Pb6fjNCvXxvQgNSufJz7dQ7cwF7f/6kynIynVYFoUWrCv9mTwp7e2cf6AcB64coiOW1wPCy7qT1pmAU9+tpeIsCCuG9PT6UhKNYgWhRZq95EcbnlxA/07tePp66N8vusKu4gID00ZytGcAv701jY6hwRy4cCWeaOf8g36l6AFSs8uID4pmaDW/iTExTTbsRCailb+fjxz/SgGR4Qwf8VGNh/KdDqSUvWmRaGFOVVUwk3LUjiRV0RCXAxdw4KcjuQT2gYGkBAXQ8fg1sQnJfPj8TynIylVL1oUWpDSMsOClZvY+lMWT10XxdBuoU5H8inhwYEkzRpDmTHEJqzjuPasqpohLQotyMPv7+STHUf582WDubiZjafcXPQNb8eS2GjSsgqYvTSF/KJSpyMpVSdaFFqI5WsOsOTr/cSNjST+XG16aqfRvdqzcHoUmw9n8ruXN1Kqw3qqZkSLQgvw+a50/vLOdi4a2Ik/Xz7Y6TgtwoShXXjgyiF8uvMof3lnm97cppoNbZLq43akZvPblzYwKCKEJ6+Lwt9P70XwlpnnRPJTZj6LvtxHRGgQ88f3czqSUrXSouDDjmS5m56GBLUiIS6GtoH66/a2P146kCNZBTz20W4iQl0tbvQ61fzoXwkflVdYwuylyeQUFPPqvLF0DnE5HalF8vMTHr1mOOnZhfzhtS10CnZpl+SqSdNrCj6otMzw+5c3sjMtm6dvGMXgriFOR2rRAgP8eW7GaPqGt2Pei+vZkZrtdCSlqqVFwQf973s7WLUrnQcmD2X8mZ2cjqOA0KBWJMXH0C4wgFlJ6/gpM9/pSEpVSYuCj0n8Zj9J3x7gpnN7M+PsXk7HURVEhAaRFB/DqcJS4hLWkXWq2OlISv2CFgUf8smOozz43g4uHdKZuycNcjqOqsLALiEsmjmaA8fzmLM8hcISvblNNS1aFHzE1sNZ/P7ljQzvFsoT07TpaVM2tm9HHp86gu/2n+D//WczZXpzm2pCtPWRD0jNzGf20mTat23N87HRBLX2dzqSqsXkkd1IyyrgkQ920TUsiHv0m51qIrQoNHM5BcXEJyWTX1TK8pvPolOwNj1tLuae14e0zHwWr95HRKiLWeO0+xHlPC0KzVhJaRm/fWkje9JzSZoVw5ldgp2OpOpARLjviiEcyS7gwfd20CXExcRhEU7HUi2cXlNopowx3PfOdr78PoOHrhrK//QPdzqSqgd/P2Hh9CiieoSx4JVNpBw44XQk1cJpUWimlny1n5e+O8jc8/vouMDNnKuVPy/ExtA9LIjZS1PYm57rdCTVgmlRaIY+3JbG3z7YyaRhXfjjpQOdjqMawRltW7M0fgyt/IXYhHWk5xQ4HUm1UFoUmplNhzK59ZVNjOgexj+vHYmfNj31GT3atyExbgwnTxUxKzGZ3MISpyOpFkiLQjNy6MQpblqaTHhwIEtio3G10qanvmZY91CeuWEUu47kcMuKDRSXljkdSbUwtRYFEekkIlNEZL6IxIvIGBHRYuJlWfnupqeFJWUkxsXQsV2g05GUTcaf2Ym/TRnK6u8zuPuNrTpAj/Kqapukish44C6gPbARSAdcwFVAXxF5DfiHMUa7fLRZcWkZ81dsYP+xPJbNHkO/Ttr01NdNi+lJamYBC1ftoWtYELdfMsDpSKqFqOk+hUnAb4wxByvPEJEA4HLgEuB1m7Ip3E1P//TmNr7ee4zHrhnO2L7aF39LcevF/UnLyufJVXuICHVpKzPlFdUWBWPMnTWsN9kYo8XAC577ch+vpBzidxf2Y2p0D6fjKC8SEf46ZRhHswv501vb6BwSyIUDOzsdS/m4+l4b+FejplBV+u+WNP7+4S6uHNFVTx+0UK38/fj3DaMYFBHM/BUb2Xwo0+lIysfVtyhoO0ibrf/xJLf9ZxPRvc7g0WuGI6KHvKVqGxhAQlwMHdq1Jj4pmR+P5zkdSfmw+hYFbQ5ho4PHTzFnWQoRoS4Wz9Smpwo6BbtYGj+GUmOITVjH8dxCpyMpH1VtURCRrSKypYrHVqDWE5sikiAi6SKyrZr5F4hIlohssh73NeDn8BlZp4qZlbSOkjJDYlwM7du2djqSaiL6hrfjhdho0rIKuGlZCvlFOkCPanw1tT66vIHbTgKeBpbVsMxXxpiG7sdnFJWUMe/F9Rw8cYoXZ59Fn/B2TkdSTczoXu1ZOH0kN6/YwO9XbuS5G0frgEqqUdV0+mgxcDUQZIz5sfKjtg0bY1YD2uWjh4wx3P3GVtbsO86j1wznrD4dnI6kmqgJQyP4y+WD+WTHUf7yzja9uU01qpqKQixwErhfRDaIyLMiMllEGvPj6zkisllEPhCRIY243Wbn6c/28vqGw9x6cX+mRHV3Oo5q4uLG9WbueX14ce1Bnv3yB6fjKB9S030KR3CfAkqyurU4C5gI/EFE8oGPjTGPNmDfG4BexphcEZkEvAX0r2pBEZkDzAHo2dP3buB5e9NP/OOT75kS1Y0FF1V5CJT6hT9OGEhqVgGPfribiFCXfphQjcKj1kfGmDJjzBpjzH3GmHHAdOCnhuzYGJNtjMm1nr8PtBKRKm/XNcYsNsZEG2Oiw8N9azCZ5AMnuPPVLYzp3Z5Hfj1Mm54qj/n5CY9PHc7Zfdrzh9e28M3eY05HUj6g1uE4ReQpftkENQtIaciORaQLcNQYY0RkDO4Cdbwh22xu9h/LY86yFLqfEcTiGaMJDNCmp6puAgP8WTQjmqnPfcu85ev5z7xzGBQR4nQs1Yx58k0hEBgJ7LEew3F3kjdbRJ6obiUReRlYA5wpIodFZLaIzBORedYi1wDbRGQz8CQw3bSgK2Yn84qIT0oGICEuhrA22vRU1U9oUCuSZo2hbWAAcYnrSM3MdzqSasaktr/DIvIZ8CtjTIn1OgD4GHdneFuNMYNtT1lBdHS0SUlp0JcUxxWWlDJjyTo2Hc7kpZvOIjqyvdORlA/YdSSbqc+uISLMxavzxhIa1MrpSKoJEZH1xpjo2pbz5JtCN6Bthddtga7GmFJAb6usI2MMf3htC+sOnODxqSO0IKhGM7BLCItmjC4/LVlYoje3qbrzpCg8CmwSkUQRScI9tsLjItIW+NTOcL7oX5/u4e1Nqdx56ZlcOaKr03GUjxnbryOPTx3Bd/tPcMerWygrazFnZFUjqfVCszHmBRF5HxiDuyO8e4wxqdbsmrrXVpW8vv4wT67aw9TR3bnlgr5Ox1E+avLIbqRlFfDIB7uICHVxz6RBTkdSzUitRcESjrsFkj9wtohgjHnDvli+Z80Px7nrjS2M7duBv07RpqfKXnPP60NaZj6LV++ja6iLuHG9nY6kmglPmqQm4G5xtB04PYq4AbQoeGhvei5zl6fQq0Nbnr1xNK0DdIhrZS8R4b4rhnAku4AH3ttBl1AXE4ZGOB1LNQOefFM429stjHzJ8dxC4pOSaR3gR2JcjLYIUV7j7ycsnB7F9c+vZcHKTay4KVAbNqhaefKRdY2IaFGoh4LiUn6zLIWj2QU8PzOaHu3bOB1JtTCuVv4siY2hW1gQs5emsDc91+lIqonzpCgsxV0Ydp8eT0FEttgdrLkrKzPc8epmNhzM5IlpI4nqeYbTkVQL1b5ta5JmjaGVvxCbsI70nAKnI6kmzJOikADMACYAV+AeZ+EKO0P5gsc/3s17W9K4e+JAJg7Tc7nKWT07tCEhLoaTp4qYlZhMbmGJ05FUE+VJUThojHnHGLO/LuMptGSvJB/k31/8wPVn9WTOeX2cjqMUAMO7h/HM9aPYdSSHW1ZsoLi0rPaVVIvjSVHYJSIvich1InL16YftyZqpr/cc4943t3HegHAevHKINj1VTcr4gZ3425ShrP4+g3ve2KoD9Khf8KT1URDu7ix+VWGaNkmtwvdHc7j5xfX069SOZ66PIsBfm56qpmdaTE9SMwtYuGoPEWFB3H7JAKcjqSbEkzuaZ3kjSHOXkVPIrMRkXK39eSEuhmCXNj1VTdetF/cnLSufJ1ftoWuoi+ljfG/wKlU/1X6UFZE/iUi1jZpF5EIRudyeWM1LflEpNy1L4UReES/ERtMtLMjpSErVSET465RhnD8gnHvf2sbnu9KdjqSaiJrOb2wF3hWRVSLymIj8QUTuE5HlIrIVdwuk77wTs+kqKzPc9somthzOZOH0kQzvHuZ0JKU80srfj3/fMIpBEcHcsmIDWw5nOh1JNQHVFgVjzNvW0JvzcHdx4Q9kAy8CY4wxtxljMrwTs+l65MNdfLj9CH+6bDC/GtLF6ThK1UnbwAAS4mLo0K418UnJHDx+yulIymG1Xgk1xuwxxiQZYx42xjxhjPnIGKNDOwEvrv2Rxav3MfOcXsSPi3Q6jlL10inYxdL4MZSUGWIT13Eir8jpSMpB2jymnr7Ync5f3tnO+DPDue/ywdr0VDVrfcPb8UJsNKmZ+cxemkx+kQ7Q01JpUaiHnWnZ/PaljQzoHMxT14/SpqfKJ4zu1Z6F06PYdCiTBSs3UqoD9LRI+tesjo5mFxCflEzbQH8S4qJpF+jpkBRKNX0Thnbh/iuG8PGOo9z/zna9ua0F8mQ8hXDgN0BkxeWNMfH2xWqa8gpLiE9KJiu/mP/MPYeIUG16qnxP7NhIUjPzWbR6H13DgrhZRwlsUTz5mPs28BXu8Zhb7InG0jLDgpUb2ZmWzfMzoxnaLdTpSErZ5o8TBpKaVcDfP3QP6XlVVDenIykv8aQotDHG/NH2JE3cQ//dwac703ngyiFcNKiz03GUspWfn/D41OFk5BRw52ub6RQcyNh+HZ2OpbzAk2sK74nIJNuTNGFJ3+wn8ZsDxI/rTezYSKfjKOUVgQH+LJoRTe+ObZm7fD0707KdjqS8wJOisAB3YSgQkRzr0WLeHat2HuXB93Zw8aDO3HvZIKfjKOVVoUGtSJo1hraBAcxKTCY1U29R8nWe3LwWbIzxM8a4rOfBxpgQb4Rz2rafsvjdyxsZ3DWEJ68bib+f3ougWp6uYUEkxceQV1hCXOI6svKLnY6kbORRk1QRuVJEHrceLaITvLQs9008YUGtSIiNoU1rbXqqWq6BXUJYNGM0+4/lMWdZCoUlLbbNic+rtSiIyCO4TyHtsB4LrGk+K7ewhPikFPIKS0mYFUOnEJfTkZRy3Nh+HXnsmhF8t/8Ed7y6hTK9uc0nefLxdxIw0hhTBiAiS4GNwF12BnNKSWkZv3tpA98fzSEhLoaBXVrEmTKlPHJVVDfSKjRVvWeSXmfzNZ6eEwkDTljPfbaBvjGGB97dwee7M/ib1de8Uurn5p3fh9TMfBav3kdEqItZ43o7HUk1Ik+KwsPARhH5HBDgPOBuW1M5JOGbAyxf+yNzz+vD9WfpSFRKVUVEuP/KIRzJLuDB93YQEepiwtAIp2OpRuJJ66OXgbNxj8n8BnCOMWal3cG87ePtR3jovzuYOLQLf5ww0Ok4SjVp/n7Ck9OjiOoRxoKVm0g5cKL2lVSzUNNwnAOtf0cBEcBh4BDQ1ZrmM7YczmTByk0M7x7GP68diZ82PVWqVkGt/VkSG0PXsCBuWpbC3vRcpyOpRiDV9YIoIouNMXOs00aVGWPMhfZGq1p0dLRJSUlptO39lJnPVc98Q2t/P96aP47w4MBG27ZSLcHB46e4+tlvcLXy541bxtIpWFvrNUUist4YE13bcjUNxznHejrRGDO+4gN3i6RmL7ugmPjEZAqKS0mcFaMFQal66NmhDQlxMRzPLSI+KZncwhKnI6kG8OTmtW89nNasFJeWMX/FBn7IyOXZG0YzoHOw05GUaraGdw/j3zeMYmdaDvNXbKC4tMzpSKqearqm0EVERgNBIhIlIqOsxwVAG68ltIExhvve3sZXe47x1ylDObe/9v6oVEONH9iJv00ZypffZ3DPG1t1gJ5mqqYmqZcCcUB34J8VpucA99S2YRFJAC4H0o0xQ6uYL8BC3KeiTgFxxpgNHidvgMWr9/HyukPcckFfpsVo01OlGsu0mJ6kZhawcNUeIsKCuP2SAU5HUnVUbVEwxiwFlorIr40xr9dj20nA08CyauZPBPpbj7OAZ61/bfX+1jQe/mAXlw+P4I5fnWn37pRqcW69uD9pWfk8uWoPXUNdTB+jH7yak1pvXjPGvC4ilwFDAFeF6Q/Wst5qEYmsYZHJwDLj/o65VkTCRCTCGJPmUfJ62HDwJLe9sonRvc7g8akjtOmpUjYQEf46ZRhHswu5961tdA5xMX5gJ6djKQ950iHec8A04He472ieCvRqhH13w33fw2mHrWm2OHTiFL9ZmkLnEBeLZ4zG1crfrl0p1eK18vfj3zeMYlBEMLes2MCWw5lOR1Ie8qT10VhjzEzgpDHmAeAcoEcj7Luqj+lVXpkSkTkikiIiKRkZGfXa2Z70HPz9hMRZMXRop01PlbJb28AAEuJi6NCuNfFJyRw8fsrpSMoDnhSFAuvfUyLSFSgGGqMHrMP8vLh0B1KrWtAYs9gYE22MiQ4Pr18ndRcO7MzqP4ynb3i7eq2vlKq7TsEulsaPoaTMEJu4jhN5RU5HUrXwpCi8KyJhwGPABuAA8HIj7PsdYKa4nQ1k2Xk9AdBTRko5oG94O5bMjCY10z1wVX6RDtDTlNVYFETED1hljMm0WiD1AgYaY+6rbcMi8jKwBjhTRA6LyGwRmSci86xF3gf2AXuB54FbGvKDKKWarujI9iycPpJNhzJZsHIjpTpAT5NVbd9H5QuIrDHGnOOlPLVq7L6PlFLek/TNfu5/dwczzu7Fg5OH4L5dSXlDg/s+quBjEfm16G9PKdVAceN6M+e8Pixf+yOLVu9zOo6qgieD7NwOtAVKRKQAd6shY4zRcSqVUnV214SBpGUV8MgH7iE9J4+0rSW6qgdPbl7TnuKUUo3Gz094fOpw0rMLuOPVzYS3C2RsP+1/rKnw5Oa1VZ5MU0opTwUG+LN4ZjS9O7Zl7vL17EzLdjqSstTUS6pLRNoDHUXkDBFpbz0iga7eCqiU8k2hQa1ImjWGNoH+zEpMJjUz3+lIipq/KcwF1gMDrX9PP94GnrE/mlLK13UNCyJp1hhyC0uIS1xHVn6x05FavJpGXltojOkN3GGM6WOM6W09RhhjnvZiRqWUDxsUEcKiGaPZfyyPuctTKCzRm9ucVOs1BWPMUyIyVkSuF5GZpx/eCKeUahnG9evIY9eMYO2+E9z56hbK9OY2x9Ta+khElgN9gU3A6RJuqH6cBKWUqrOrorqRmpXPox/uJiLUxd2TBjkdqUXy5D6FaGCw0bH1lFI2u/n8vqRlFrBo9T66hgUROzbS6UgtjidFYRvQBbC1szqllBIR7r9yCEeyC7j/3e10DglkwtAIp2O1KJ50c9ER2CEiH4nIO6cfdgdTSrVM/n7Ck9OjGNkjjAUrN5Fy4ITTkVoUTzrEO7+q6caYL21JVAvtEE+pluFEXhG/fvZbTp4q4vWbx+pYKA3UaB3iWX/8DwCtrOfJuMdVUEop27Rv25qkWTEE+AmxCetIzymofSXVYJ50c/Eb4DVgkTWpG/CWnaGUUgqgV4e2vBAbw/HcIuKTksktLHE6ks/z5JrCfGAckA1gjNkDdLIzlFJKnTaiRxjP3BDFzrQc5q/YQHFpmdORfJonRaHQGFM+sKqIBOC+T0EppbziwoGdeeiqoXz5fQb3vrkVbSFvH0+apH4pIvcAQSJyCe5hM9+1N5ZSSv3cdWN6kpaZz5Of7SUiNIjbLhngdCSf5ElRuAuYDWzF3Une+8ASO0MppVRVbrtkAKlZBSxctYeuYS6mxfR0OpLP8aQoBAEJxpjnAUTE35p2ys5gSilVmYjw8NXDOJpdwD1vbqNTsIvxA/USZ2Py5JrCKtxF4LQg4FN74iilVM1a+fvx7I2jGdglmFtWbGDL4UynI/kUT4qCyxiTe/qF9byNfZGUUqpm7QIDSIyLoX3b1sQnJXPwuJ64aCyeFIU8ERl1+oWIjAZ0iCSllKM6hbhYGh9DcakhNnEdJ/KKal9J1cqTonAr8KqIfCUiXwGvAL+1N5ZSStWuX6dglsRG81NmPjctTaagWAfoaShPurlIxj0k5824m6MOMsastzuYUkp5IiayPQunjWTjoUx+//JGSnWAngbx5JsCQAwwHIgCrtOR15RSTcnEYRH8+bLBfLzjKA+8u11vbmsAHXlNKeUT4s/tTVpWPs9/tZ+uYUHMO7+v05GaJR15TSnlM+6eOIi0rAIe+WAXEaEuJo/s5nSkZseT00enR15TSqkmzc9P+Me1Izird3vueHUz3/5wzOlIzY6OvKaU8imBAf4snhlN745tmbtsPbuOZDsdqVnRkdeUUj4pNTOfKf/+BkF4c/5YIkKDal/JhzX2yGu7gGDrsdOpgqCUUp7qGhZE0qwx5BaWEJeQTFZ+sdORmgVPRl67FlgHTAWuBb4TkWvsDqaUUg01KCKERTNG80NGLnOXp1BYoje31cZXBIshAAARCUlEQVSTawr3AjHGmFhjzExgDPBne2MppVTjGNevI49NHc7afSe489UtlOnNbTXypEmqnzEmvcLr43h+05tSSjluSlR30rIKePTD3USEubh74iCnIzVZnhSFD0XkI+Bl6/U04AP7IimlVOO7+fy+pGbms+jLfXQNDSJ2bKTTkZqkWouCMeZOEbkaOBcQYLEx5k3bkymlVCMSER64cihHswu5/93tdA5xMWGo3oJVWbWngUSkn4iMAzDGvGGMud0YcxtwXET0/nGlVLPj7yc8OT2KEd3DWLByI+t/POF0pCanpmsDTwA5VUw/Zc2rlYhMEJHdIrJXRO6qYn6ciGSIyCbrcZNnsZVSqn6CWvvzQmw0EaEuZi9N4YeM3NpXakFqKgqRxpgtlScaY1KAyNo2bI3l/AwwERiMu3fVwVUs+ooxZqT1WOJZbKWUqr8O7QJZGj8GfxFiE9aRnlPgdKQmo6ai4Kphnie3Bo4B9hpj9hljioCVwOS6hFNKKbv06tCWhLgYjucWMTsphbzCEqcjNQk1FYVkEflN5YkiMhvwZJCdbsChCq8PW9Mq+7WIbBGR10SkhwfbVUqpRjGiRxjP3BDF9tQs5r+0geLSMqcjOa6monArMEtEvhCRf1iPL4GbgAUebFuqmFb5rpF3cZ+mGg58CiytckMic0QkRURSMjIyPNi1Ukp55sKBnfnrlGF8sTuDP725rcUP0FNtk1RjzFFgrIiMB4Zak/9rjPnMw20fBip+8u8OpFbax/EKL58H/l5NlsXAYnB3iOfh/pVSyiPXjelJWmY+T362l4gwF7dePMDpSI7x5D6Fz4HP67HtZKC/iPQGfgKmA9dXXEBEIowxadbLK4Gd9diPUko12G2XDCA1q4AnPt1D19Agro1pmWezPbmjuV6MMSUi8lvgI8AfSDDGbBeRB4EUY8w7wO9F5EqgBDgBxNmVRymlaiIiPHz1MI5mF3D3m1sJDwlk/JmdnI7ldbWOp9DU6HgKSik75RaWMG3RGvYfy+OVOecwrHuo05EaRaONp6CUUi1Ju8AAEuNiOKNNa2YlJXPoxCmnI3mVFgWllKqkU4iLpfExFJeWEZuwjpN5RU5H8hotCkopVYV+nYJZEhvN4cx8Zi9NpqC4ZQzQo0VBKaWqERPZnoXTRrLxUCYLVm6ktAUM0KNFQSmlajBxWAR/vmwwH20/yoPvbvf5m9tsa5KqlFK+Iv7c3qRl5fP8V/vpGhbE3PN9d/QALQpKKeWBuycOIi2rgIc/2EWXUBeTR1bVlVvzp0VBKaU84Ocn/OPaEWTkFHLHq5sJDw5kbN+OTsdqdHpNQSmlPBQY4M/iGdFEdmjL3GXr2XUk2+lIjU6LglJK1UFom1YkxY+hTaA/cQnJpGbmOx2pUWlRUEqpOuoWFkRi3BhyC0uYlZhMVn6x05EajRYFpZSqh8FdQ3juxtH8kJHL3OUpFJb4xs1tWhSUUqqezu3fkUevGc7afSe489UtlPnAzW3a+kgppRrg6lHdScsq4LGPdhMR5uLuiYOcjtQgWhSUUqqBbrmgL2lZ+Sz6ch9dQ4OIHRvpdKR606KglFINJCI8cOVQjmQVcv+72+kc4mLC0C5Ox6oXvaaglFKNwN9PeOq6KEb2CGPByo2s//GE05HqRYuCUko1kqDW/iyZGU1EqIvZS1P4ISPX6Uh1pkVBKaUaUYd2gSyNH4O/CHGJ68jIKXQ6Up1oUVBKqUbWq0NbEuJiOJZTRHxSMnmFJU5H8pgWBaWUssGIHmE8fX0U21OzmP/SBkpKy5yO5BEtCkopZZOLBnXmoauG8cXuDO59c1uzGKBHm6QqpZSNrj+rJ2lZ+Tz12V66hgWx4OL+TkeqkRYFpZSy2e2XDCA1s4B/ffo9EWEuro3u4XSkamlRUEopm4kIj/x6GOk5Bdz9xlY6BQdywZmdnI5VJb2moJRSXtDK349nbxzNmZ2DuWXFBrb9lOV0pCppUVBKKS9pFxhA0qwYzmjTmrjEZA6dOOV0pF/QoqCUUl7UKcTF0vgYikvLiE1cx8m8Iqcj/YwWBaWU8rJ+nYJZEhvN4ZP53LQshYLipjNAjxYFpZRyQExke56YNpINB0+yYOVGSpvIAD1aFJRSyiGThkXw58sG89H2o/zvezuaxM1t2iRVKaUcFH9ub1Iz81ny9X66hrmYc15fR/NoUVBKKYfdM2kQadkF/O39XXQOcTF5ZDfHsmhRUEoph/n5Cf+YOoKMnELueHUz4cGBjO3b0ZksjuxVKaXUz7ha+fP8jGgiO7Rl7vL17D6S40gOLQpKKdVEhLZpRVL8GIJa+ROXuI60rHyvZ9CioJRSTUi3sCASZ8WQU1DCrMRksguKvbp/W4uCiEwQkd0isldE7qpifqCIvGLN/05EIu3Mo5RSzcGQrqE8d+No9qbnMm/5eopKvDdAj21FQUT8gWeAicBg4DoRGVxpsdnASWNMP+BfwN/tyqOUUs3Juf078ug1w/n2h+P84bXNlHnp5jY7vymMAfYaY/YZY4qAlcDkSstMBpZaz18DLhIRsTGTUko1G1eP6s6dl57JW5tSefSj3V7Zp51NUrsBhyq8PgycVd0yxpgSEckCOgDHbMyllFLNxi0X9CU1M5/nvvyBrmEuZp4Taev+7CwKVX3ir/z9x5NlEJE5wByAnj17NjyZUko1EyLCA1cOIbughG5hQbbvz86icBioOOZcdyC1mmUOi0gAEAqcqLwhY8xiYDFAdHS0852DKKWUFwX4+/HUdVFe2Zed1xSSgf4i0ltEWgPTgXcqLfMOEGs9vwb4zDSFHqGUUqqFsu2bgnWN4LfAR4A/kGCM2S4iDwIpxph3gBeA5SKyF/c3hOl25VFKKVU7W/s+Msa8D7xfadp9FZ4XAFPtzKCUUspzekezUkqpcloUlFJKldOioJRSqpwWBaWUUuW0KCillConze22ABHJAH6s5+odaZpdaDTVXNB0s2muutFcdeOLuXoZY8JrW6jZFYWGEJEUY0y00zkqa6q5oOlm01x1o7nqpiXn0tNHSimlymlRUEopVa6lFYXFTgeoRlPNBU03m+aqG81VNy02V4u6pqCUUqpmLe2bglJKqRr4TFEQkQkisltE9orIXVXMDxSRV6z534lIZIV5d1vTd4vIpV7OdbuI7BCRLSKySkR6VZhXKiKbrEflbsftzhUnIhkV9n9ThXmxIrLHesRWXtfmXP+qkOl7EcmsMM/O45UgIukisq2a+SIiT1q5t4jIqArz7DxeteW6wcqzRUS+FZERFeYdEJGt1vFK8XKuC0Qkq8Lv674K82p8D9ic684KmbZZ76n21jxbjpeI9BCRz0Vkp4hsF5EFVSzjvfeXMabZP3B3zf0D0AdoDWwGBlda5hbgOev5dOAV6/lga/lAoLe1HX8v5hoPtLGe33w6l/U618HjFQc8XcW67YF91r9nWM/P8FauSsv/DneX7LYeL2vb5wGjgG3VzJ8EfIB7NMGzge/sPl4e5hp7en/AxNO5rNcHgI4OHa8LgPca+h5o7FyVlr0C9xgvth4vIAIYZT0PBr6v4v+j195fvvJNYQyw1xizzxhTBKwEJldaZjKw1Hr+GnCRiIg1faUxptAYsx/Ya23PK7mMMZ8bY05ZL9fiHqHObp4cr+pcCnxijDlhjDkJfAJMcCjXdcDLjbTvGhljVlPFqIAVTAaWGbe1QJiIRGDv8ao1lzHmW2u/4L33lyfHqzoNeW82di6vvL+MMWnGmA3W8xxgJ+7x6yvy2vvLV4pCN+BQhdeH+eVBLV/GGFMCZAEdPFzXzlwVzcb9aeA0l4ikiMhaEbmqkTLVJdevra+qr4nI6aFVm8Txsk6z9QY+qzDZruPlieqy23m86qry+8sAH4vIenGPg+5t54jIZhH5QESGWNOaxPESkTa4/7i+XmGy7cdL3Ke1o4DvKs3y2vvL1kF2vEiqmFa5WVV1y3iybn15vG0RuRGIBs6vMLmnMSZVRPoAn4nIVmPMD17K9S7wsjGmUETm4f6WdaGH69qZ67TpwGvGmNIK0+w6Xp5w4v3lMREZj7sonFth8jjreHUCPhGRXdYnaW/YgLvbhVwRmQS8BfSniRwv3KeOvjHGVPxWYevxEpF2uIvQrcaY7Mqzq1jFlveXr3xTOAz0qPC6O5Ba3TIiEgCE4v4a6cm6duZCRC4G7gWuNMYUnp5ujEm1/t0HfIH7E4RXchljjlfI8jww2tN17cxVwXQqfbW38Xh5orrsdh4vj4jIcGAJMNkYc/z09ArHKx14k8Y7bVorY0y2MSbXev4+0EpEOtIEjpelpvdXox8vEWmFuyCsMMa8UcUi3nt/NfZFEyceuL/x7MN9OuH0xakhlZaZz88vNP/Hej6En19o3kfjXWj2JFcU7gtr/StNPwMItJ53BPbQSBfcPMwVUeH5FGCt+b8LW/utfGdYz9t7K5e13Jm4L/qJN45XhX1EUv2F08v4+YXAdXYfLw9z9cR9nWxspeltgeAKz78FJngxV5fTvz/cf1wPWsfOo/eAXbms+ac/MLb1xvGyfu5lwBM1LOO191ejHWinH7ivzn+P+w/svda0B3F/+gZwAa9a/0HWAX0qrHuvtd5uYKKXc30KHAU2WY93rOljga3Wf4qtwGwv53oY2G7t/3NgYIV1463juBeY5c1c1uv7gUcqrWf38XoZSAOKcX86mw3MA+ZZ8wV4xsq9FYj20vGqLdcS4GSF91eKNb2Pdaw2W7/ne72c67cV3l9rqVC0qnoPeCuXtUwc7sYnFdez7XjhPqVngC0Vfk+TnHp/6R3NSimlyvnKNQWllFKNQIuCUkqpcloUlFJKldOioJRSqpwWBaWUUuW0KCillCqnRUGpOhKReSIysw7LR4nIkmrmHbDu5PV0W5eLyAOeLq9UXel9CkrZTEReBR4yxmyuYt4B3DciHfNwW4K736Bx5v9611Wq0eg3BeWzRCTG6uXVJSJtrQFMhlax3BXiHnhpo4h8KiKdrelPnh78RUQuFZHVIuInIveLyB3W9N/L/w2StLKKbQcDw08XBBHpICIfW/tahNWhmYhEisguEVliDe6yQkQuFpFvrMFTxgAY96e4L4DLbTloqsXToqB8ljEmGXgHeAh4FHjRGFPViFtfA2cbY6Jw99//B2v6XcA0q4fRJ3F3IVBWad27gChjzHDc3RJUFg1U3OdfgK+tfb2Du2+i0/oBC4HhwEDgetxdINwB3FNhuRTgf2r40ZWqN1/pOlup6jwIJAMFwO+rWaY78Io1aElr3J2KYYw5JSK/AVYDt5mqu+HeAqwQkbdwd/9cWQSQUeH1ecDV1vb/KyInK8zbb4zZCiAi24FVxhgjIltxd+J2WjrQtfofWan6028Kyte1B9rhHubQBSAifz09Dq+1zFO4hx4dBsw9vZxlGHCc6v8IX4a7o7LRwHqrW/aK8ittD6rv776wwvOyCq/L+PkHOJe1XaUanRYF5esWA38GVgB/BzDG3GuMGWmMGWktEwr8ZD0vH/jcGt3t/+Hu3nyiiJxVccMi4gf0MMZ8jvuUUxjuAlTRTtynhU5bDdxgrT8Rd3fHdTWAn5+SUqrRaFFQPstqNlpijHkJeASIEZELq1j0fuBVEfkKOGatK8ALwB3GPbjKbGCJiFT81O8PvGid3tkI/MsYk1lxw8aYXUCodcEZ4AHgPBHZAPwK9zgCdTUe+G891lOqVtokVSmbichtQI4xpsp7Feq4rc7AS8aYixqeTKlf0m8KStnvWX5+vaAheuI+paWULfSbglJKqXL6TUEppVQ5LQpKKaXKaVFQSilVTouCUkqpcloUlFJKlfv/5Rlrw3ZsWPQAAAAASUVORK5CYII=\n",
"text/plain": [
"