{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Heat advection and diffusion\n", "\n", "## General equation\n", "\n", "If we add advection and internal heat generation to our heat equation, we get\n", "\n", "\\begin{equation}\n", " \\rho C_p \\left(\\frac{\\partial T}{\\partial t} + v_x\\frac{\\partial T}{\\partial x}\\right) = k \\frac{\\partial^2 T}{\\partial x^2} + H\n", "\\end{equation}\n", "\n", "If we rearrange the terms to have only the time derivative on the left-hand side, we get\n", "\n", "\\begin{equation}\n", " \\frac{\\partial T}{\\partial t} = \\frac{k}{\\rho C_p} \\frac{\\partial^2 T}{\\partial x^2} + \\frac{H}{\\rho C_p} - v_x\\frac{\\partial T}{\\partial x}\n", "\\end{equation}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise - Advection and diffusion\n", "\n", "- Write the finite difference approximation (i.e., *discretize*) the advection-diffusion equation. Use the forward difference in time and the central difference in space. Write the expression in form\n", "\n", " \\begin{equation}\n", " T_i^{q+1} = \\ldots~T^q~\\mathrm{terms~only~on~this~side}~\\ldots\n", " \\end{equation}\n", "\n", "- Modify the intrusion code below (or use the one you modified previously):\n", " - Add variables that store the internal heat production $H$ and the advection velocity $v_x$ (what are suitable values?)\n", " - Modify the line within the time loop that calculates the next value of $T$: include the terms for heat generation and advection\n", "- Start with $H=0$ and $v_x = 0$. Gradually change these values.\n", " - How do the results change as you give non-zero values for these terms?\n", " - What happens if $v_x > 0$ compared to $v_x < 0$?\n", " - What if $v_x$ has a large (positive or negative) value?\n", " \n", "#### Extra tasks\n", "\n", "- Remove the initial hot intrusion. Instead, add spatially variable internal heat generation. For this, you will need to introduce new numpy array (of size `nx`), fill it with values of internal heat generation, and use it (`H[ix]` instead of `H`) when you calculate the new temperature value.\n", "- Calculating the heat flow\n", " - Heat flow is defined as $q=k\\frac{dT}{dx}$. Write an approximation for $q$ by replacing the derivative with a finite difference expression.\n", " - Modify the code so that you create an array (of size `nt`) that stores surface heat flow values each time step. Within the time loop, calculate each time step the heat flow between two topmost grid points.\n", " - After the time loop, plot the recorded surface heat flow history." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "dx is 1206.896551724138 m\n", "Calculating 10 time steps\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAD8CAYAAACPWyg8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VOXZ//HPNTNZ2Jew74RNFgEBNbgWq4K2Fa0bKpXWKrXi0qe1j/rYWp9uv+4qro1bpQXRuoGKj4JaFUvEsMi+xEgkgAohIBXINvfvjzkTJmGykWVmku/79ZpXztxnmSu8wrnmvu/rnGPOOUREROrCF+sAREQk8Sh5iIhInSl5iIhInSl5iIhInSl5iIhInSl5iIhInSl5iIhInSl5iIhInSl5iIhInQViHUBj6dKlixswYECswxARSSgrVqzY45zrWtN2zTZ5DBgwgOzs7FiHISKSUMwsrzbbadhKRETqTMlDRETqTMlDRETqrNnOeYhIy1ZSUkJ+fj6HDx+OdShxKTU1lT59+pCUlHRM+yt5iEizlJ+fT7t27RgwYABmFutw4opzjoKCAvLz8xk4cOAxHSNhhq3MbIqZbTazHDO7PdbxiEh8O3z4MGlpaUocUZgZaWlp9eqVJUTyMDM/8CBwHjACuMLMRjTGZ63IK+TBt3NYkVfYGIcXkSakxFG1+v7bJMqw1UlAjnMuF8DM5gNTgQ0N+SEr8gq58tEsikuDpCT5mHttBuP7d2rIjxARaRYSoucB9Aa2R7zP99oqMLOZZpZtZtm7d++u84dk5RZQXBrEAYdLgize8NkxBywiLVtBQQFjx45l7Nix9OjRg969e5e/Ly4uLt/u/PPPZ9++fdUe66677mLJkiUA3HvvvRw8eLBO+w8YMIA9e/bU47c5WqL0PKL1r9xRDc5lApkAEyZMOGp9TTLS00hJ8lFcGiToYG5WHhMHdeHMoTVeqS8iUkFaWhqrV68G4O6776Zt27bceuutR223aNGiGo/1y1/+snz53nvvZfr06bRu3brW+zeGROl55AN9I973AXY29IeM79+Juddm8JNzhzH7irH07tSa7z65nD+/sZmyYJ1zkYgkmFjMeYZ7Bdu2bWP48OFcd911jBw5knPPPZdDhw4B8N3vfpfnnnuO2bNns3PnTiZNmsSkSZMq7A9w4YUXMn78eEaOHElmZmajxp0oPY8PgSFmNhDYAUwDrmyMDxrfv1P5PMc5w3vwi4XruP+tHLK3FXLfFWPp1i61MT5WRBrR/768ng07v6x2mwOHS9j02QGCDnwGx/VoR7vUqq+BGNGrPb/41sgGjXPr1q08/fTTPProo1x22WU8//zzTJ8+vXz9zTffzF/+8hfefvttunTpctT+TzzxBJ07d+bQoUOceOKJXHzxxaSlpTVojGEJ0fNwzpUCNwKvAxuBZ51z6xv7c1sl+/nDJWP44yWjWbW9kPPvW8q/P27YcUMRiQ9fHi4lPMAQdKH3TW3gwIGMHTsWgPHjx7Nt27Y67T979mzGjBlDRkYG27dvZ+vWrY0QZUii9Dxwzi0CYjK4d+mEvozu05Efzl3B9Mc+YNqJfenVsRUTB3VRNZZIAqhND2FFXiFXPZZFSWmQpICP+6ad0OT/v1NSUsqX/X5/+bBVbfzrX/9iyZIlLFu2jNatW/O1r32tUa+uT5jkEWvDerTj5RtP4/p/rGDe8lDhV0ogh3nXqZxXpDkIz3lm5RaQkZ4Wt/+v27Vrx4EDB44attq/fz+dOnWidevWbNq0iaysrEaNIyGGreJFm5QAGemdy0u/ikqDvLAyP6YxiUjDGd+/E7MmDY7bxAEwc+ZMzjvvvPIJ87ApU6ZQWlrK6NGj+fnPf05GRkajxmHONc8qogkTJrjGeBhUuGtbXBrEeRNrt583nGtPH6irWUXiyMaNGxk+fHisw4hr0f6NzGyFc25CTftq2KqOIru2o/t0YG7Wp/xm0UaWb9vLny4ZQ4fWx3aHShGRRKLkcQwiy3lPG9yFJ9/fxm8XbeQb97/HQ1eNY3SfjjGOUESkcWnOo57MjGtOG8iz10/EObjk4WXMWbaN5jocKJJI9P+wavX9t1HyaCDj+nXilZtO47QhXbhrwXqueuwD7lm8WXfnFYmR1NRUCgoKlECiCD/PIzX12C961oR5AwsGHb9YuJ6/Z+UBkBzw8bTKeUWanJ4kWL2qniSoCfMY8fmMHh1S8VnoKtXi0iCZ737MI9PHqxpLpAklJSUd81PypGYatmoEGelpJAd8+C1Uyvv6+s/5ybMfcbC46W93ICLSGNTzaASR5bwnDezM+zl7uO/NrazdsZ+Hp49jcLd2sQ5RRKReNOfRRJZu3cMt81dxsLiM3357FBed0CfWIYmIHKW2cx4atmoipw3pwqJbTuf4Ph34r2c+4o4X1nC4pCzWYYmIHBMNWzWh7u1TmXftyfxl8RYe+tfH/PvjAs4Z0Z3zRvVUNZaIJBT1PJpYwO/jv6ccx53nDyev4CCPvfcJ0zKX6XoQEUkoSh4xUlwWxOdV7paUOX7z6gaKSjWMJSKJQckjRiLLef0+Y+Wn+7jskWVs33sw1qGJiNRI1VYxtCKvsPzBM7sPFPHTf36Ez2f8+dIxnD2ie6zDE5EWSFeYJ4DIu/MCDO/ZjhvmruTaOdn84Ix0bp08jCS/OociEn8a7cxkZneb2Q4zW+29zo9Yd4eZ5ZjZZjObHNE+xWvLMbPbI9oHmtkHZrbVzJ4xs+TGijuW+qe14fkfnsL0jH789d1crsjMYtf+2j/DWESkqTTasJWZ3Q38xzn3p0rtI4CngZOAXsASYKi3egtwDpAPfAhc4ZzbYGbPAi845+ab2SPAR865h6v7/EQYtqrOgtU7uOOFtfh9xplDu3Lq4C6M7tOBlICPJL+P5Iifyf7Qy+fTvbNEpH7iedhqKjDfOVcEfGJmOYQSCUCOcy4XwMzmA1PNbCNwFnClt81TwN1Atckj0U0d2xsDbpm/mlfW7OKVNbtq3MfvM5L9PpL8RnLAT7LfKiSZCskmcGS7JL8dSUp+H0kR20QeL8k7XnKl4yX5fZWS2tHbBXymG0OKNCONnTxuNLOrgWzgJ865QqA3kBWxTb7XBrC9UvvJQBqwzzlXGmX7CsxsJjAToF+/fg31O8TM9sJDmFH+rPQLxvTi68O7U1IWpLg0SElZkKLSICVlrvx9sbeuuCxISfhneZujuLSMQyVl7D9UEtF+5HjF4eOVBRv0dzEjlGQiklNSIJzsKiefaEnu6J7WkeOEj2sk++ue5JJ86rWJ1FW9koeZLQF6RFl1J6Gewa8A5/38M3ANEO1/qSP6/IurZvujG53LBDIhNGxVQ/hxL1zOW1IaJCng4zsTBzTZlejOufIkUlJaKSlVl7yiJqPI5FVp/yhJ7qvisir2P/KzoUdbk/xWMXlF9rwiklz1SSlKkovao6uY5GpKnn4lNolD9Uoezrmza7OdmT0KvOK9zQf6RqzuA+z0lqO17wE6mlnA631Ebt+sRd6dNyM9rUlvYWJmoZNcwAcpTfaxtVZadiRpFZcFoye5qD2vaMmv5iQXPu5/ikor9tCiJMvSYMNmNp9RdZI5qudVfZILtR89nBl92NJLctX0EJP8Go5sqRpt2MrMejrnwgP1FwHrvOWFwDwz+wuhCfMhwHJCPYwhZjYQ2AFMA650zjkzexu4BJgPzAAWNFbc8aZyOa+EBPw+An5oleyPdShHCQZdpR5alCTjJbdoPa+oyauK4cXI5FdUEuTA4dIoPcSKn9/Qqh1erDIpVTecWHWSi3bco5KcikiaRGPOefzBzMYSGmLaBvwAwDm33que2gCUArOcc2UAZnYj8DrgB55wzq33jnUbMN/Mfg2sAh5vxLhF6sXnM1J9flKT4i+xOecoDVY1zBi9JxVORLXupUVNlI6Dh0qOSpSVj9vAnTYCPqtxeLH6nlcNycs7XvVzdhHbRSQ5fyMVkURefNyYXzx1hbmIxI2yYPQEdnTychSXlVFc6qoYoqy+l1Zdkqtqfq+krGHPlWaU95AqJ5kaqyUjC0QCFhq29Pv44kARTy//lKBzJAd8zL02o84JJJ5LdUVEovL7jFbJfloRf722YNBREqzYS6u6R1YxyZWUuopDlBFzdUdVOkZJfl8VlR75nIgkFx76jDYcWVIaJCu3oNF6H0oeIiK14PMZKT4/KYH4S2zh4cjF6z/nhnkrMSAp4CMjPa3RPlPJQ0QkwZkZSX5jwsBQL+PsEd25/sxBjTrnobvuiYg0M2cO7droVZpKHiIiUmdKHiIizcw7W3Y3+qOtlTxERJqJNfn7AViy4XOueiyrUROIkoeISDOx0ksWjiOluo1FyUNEpJkY502Sq1RXRERqbXSfDoBKdUVE5BioVFdEROKSkoeISDOjUl0REak1leqKiEidqVRXRETqTKW6IiJSZyrVFRGRY6ZSXRERiUv1Sh5mdqmZrTezoJlNqLTuDjPLMbPNZjY5on2K15ZjZrdHtA80sw/MbKuZPWNmyV57ivc+x1s/oD4xi4g0d01RqlvfOY91wLeBv0Y2mtkIYBowEugFLDGzod7qB4FzgHzgQzNb6JzbAPweuMc5N9/MHgG+Dzzs/Sx0zg02s2nedpfXM24RkUbhXOj54kWlQYpKghSVllFc6r0vDVJUUkZRaTCiray8vbgsvE9ke9A7XlmFddGOebC4DAiV6r63dTdzr82Iz2eYO+c2QugRiJVMBeY754qAT8wsBzjJW5fjnMv19psPTDWzjcBZwJXeNk8BdxNKHlO9ZYDngAfMzJxzrj6xi0jzE36Wd+QJ9ciJuqzCCb2q9gon5fAJu6zifhX2rXSyLy4N1vv38BmkJvlJDvhICfhICfhDP5NCy8l+H23aBCqsSw742PTZAVbmFVYo1Y3L5FGN3kBWxPt8rw1ge6X2k4E0YJ9zrjTK9r3D+zjnSs1sv7f9nsYJXUSOVWlZ5Im18rfno0+0lb851/lEH+Wben2/Vppx1Em5/H1SaLlDqySvrWJ7cuSJPuAjJclPit9Xvr7iMf0V2o98jo+A/9hmFFbkFXLVY1mUlAZjX6prZkuAHlFW3emcW1DVblHaHNHnWFw121d3rKM/1GwmMBOgX79+VYQm0jyVBV35CTby5H24JMpJudrhk+gn+grbVXHMsmD9BwSifdtO9nsn4oCPtikB0tpEnqCrP9FXPKFXf6JP9vtI8lu00ZSEML5/J+Zem0FWbgEZ6WmNWnFVY/Jwzp19DMfNB/pGvO8D7PSWo7XvATqaWcDrfURuHz5WvpkFgA7A3ipizQQyASZMmKBhLWkyzrmjhjqKaznUUZuTcm2OWVLWACduv6/CN+DwCTvc1irJT8dWSUdO6JW+PVc4gScdOWFXd8wKJ3S/L2FP3PFifP9OjV6mC403bLUQmGdmfyE0YT4EWE6oFzHEzAYCOwhNql/pnHNm9jZwCTAfmAEsiDjWDGCZt/4tzXdIJOccJWXuqG/VRw1pVDHxeKS9LsMnFcfDG2KcO+CzKoc0wu3tUgM1fqtOLh8mOfKtuqpjVt7P59OJW2qnXsnDzC4C7ge6Aq+a2Wrn3GTn3HozexbYAJQCs5xzZd4+NwKvA37gCefceu9wtwHzzezXwCrgca/9ceDv3qT7XkIJR+LYim17eXXdLoZ3b8/Arm1q/lZdUhYxIRnlJF3hhH70sEtRA5y4zSC1uiGNgJ+OrZMrfaOutF0Nwycp1Zzok/3HPs4tEgvWXL/ET5gwwWVnZ8c6jBZnRV4h0zKX1XkIpfKQRugE7D+qPWr1ib/ifhWGSZIqnrCjtScHfAR8iTvOLdKQzGyFc25CTdvp3lbSoLJyCypMmnZsncQ9l42lc5vkoyY+wyf0RJ6gFGmp1E+WBpWRnkZywIffQmP4+w+W8KtXNpCS5OO4Hu1J79qW3h1b0aVtCu1Sk0gOaIJUJBFp2Eoa3Iq8wvJSwZKyIDc9vYovD5XwqwtHcdmEvjUfQERiprbDVkoe0uh2HyjiR8+s4v2cAi4Z34dfTR1Fq2R/rMMSkSg05yFxo2u7FOZcczL3vbmV+9/aytr8/cyaNIjthYca/UImEWkc6nlIk3p3y25mzVvJgcOlGJCS5GvUm7eJSN3UtuehCXNpUmcM7cr0k/sDoXvMFJUEWbp1d2yDEpE6U/KQJnf2iO6kJvkwQgnkpdU7yCv4KtZhiUgdKHlIkwvfvO3WycO447zj2PtVCd+cvZT/W7cr1qGJSC1pzkNiLr/wILPmreKj7fv43qkDuOO84SQH9L1GJBY05yEJo0+n1vzzBxP53qkDePL9bVz212Xs2Hco1mGJSDXU85C4smjtLv77uTUE/MYNXxtESZlTOa9IE9J1HpKQzj++JyN6tue7Ty7nt4s2qZxXJE5p2ErizoAubbjohNBTiB1wuCTIko2fxzYoEalAyUPi0mlDupKa5CP8bKK5WXm8p+tBROKGkofEpXA570/OHcbsaWPp0SGVq59Yzj2LtzTIc7JFpH405yFxK/JZzGeP6M7PXlrHfW9uJTtvL/defgJd26XEOEKRlks9D0kIrZMD/PnSMfzh4tFkbyvkG7Pf44PcgliHJdJiqechCcPMuOzEvhzfpwM3zF3JFY9mccVJ/ejZIZWJg7qoGkukCdWr52Fml5rZejMLmtmEiPYBZnbIzFZ7r0ci1o03s7VmlmNms817jJyZdTazxWa21fvZyWs3b7scM1tjZuPqE7MkvuE927PwxlPJSE9j7gef8qc3tnDlo1msyCuMdWgiLUZ9h63WAd8G3o2y7mPn3FjvdX1E+8PATGCI95ritd8OvOmcGwK86b0HOC9i25ne/tLCtUtN4tTBaYQfYFtUGuSFlfkxjUmkJalX8nDObXTOba7t9mbWE2jvnFvmQpe2zwEu9FZPBZ7ylp+q1D7HhWQBHb3jSAuXkd6FFK+c14D5yz/l8aWf0FzvmiASTxpzzmOgma0CvgR+5px7D+gNRH49zPfaALo753YBOOd2mVk3r703sD3KProFawsXLufNyi3g+F4d+PsHefzqlQ0s/6SAP1wyhg6tkmIdokizVWPyMLMlQI8oq+50zi2oYrddQD/nXIGZjQdeMrORUD7KEKmmr4m13sfMZhIa2qJfv341HFaag8hy3tOHduHxpZ/wu9c28a37l/LQVeMY1btDjCMUaZ5qHLZyzp3tnBsV5VVV4sA5V+ScK/CWVwAfA0MJ9Rr6RGzaB9jpLX8eHo7yfn7htecDfavYp/LnZjrnJjjnJnTt2rWmX02aGTPj2tPTeeYHGZSUBfn2Q//mH1l5GsYSaQSNcp2HmXU1M7+3nE5osjvXG5Y6YGYZXpXV1UA4CS0EZnjLMyq1X+1VXWUA+8PDWyLRjO/fmVdvPp2Jg9L42Uvr+M7joSvTVY0l0nDqW6p7kZnlAxOBV83sdW/VGcAaM/sIeA643jm311v3Q+AxIIdQj+Q1r/13wDlmthU4x3sPsAjI9bZ/FLihPjFLy9C5TTJPfvdErjy5H0tz9nDfm1tVzivSgOo1Ye6cexF4MUr788DzVeyTDYyK0l4AfD1KuwNm1SdOaZl8PqN3x1b4DIIuVM772Hu5jO8/PtahiSQ83Z5EmrWM9DSSAz78Bj6D19Z9xk//+RGHistiHZpIQtPtSaRZiyznPWlgZ97bspv7385hTf5+Hpo+jkFd28Y6RJGEpMfQSovzzpbd/NczqykqKeO33z6eqWN717yTSAtR28fQathKWpwzh3bl1ZtPY3jP9twyfzV3vriWwyUaxhKpCw1bSYvUs0Mrnp6ZwZ/e2Mxf38ll2cd7OHt4DyaP6qG784rUgnoe0mIl+X3ccd5wbj/vOHL3HCTzvVymZS5TOa9ILSh5SItXFnTlz0ovKXP8v0UbKS4NxjYokTin5CEtXmQ5r99nZOcVcnnmMnbsOxTr0ETilqqtRIAVeYVk5RaQkZ7GZ/sPc9vzawj4jXsuG8uk47rVfACRZqK21VaaMBeh4t15AUb0as8Nc1fyvb99yA+/NoifnDOUgF8ddZEwJQ+RKAZ2acOLN5zC/768gYf/9TEr8gqZeXo6mz8/QEZ6miqypMXTsJVIDV5clc9tz62luCyIzyA54GPutRlKINIs6SJBkQZy0Ql9mJ4RerhY0EFRSZBlH++JcVQisaXkIVIL3xjdi9RA6L+LAxZv+Jw9/ymKbVAiMaTkIVIL4/t3Yu51Gfx08lB+eGY6mz47wDdmv8fyT/bWvLNIM6Q5D5FjsGHnl8yat5JP9x7kp5OHMfP0dHzhKw1FEpjmPEQa0Yhe7Vl446lMGdmD3722ievmZLPvYHGswxJpMirVFTlG7VKTeODKEzg5qzO/emUD35i9lBvPGszer4pVzivNnpKHSD2YGVdPHMCYPh25ds6H3PHCWgxISVI5rzRvGrYSaQBj+nbkihND5bwOOFwS5J0tX8Q2KJFGVK/kYWZ/NLNNZrbGzF40s44R6+4wsxwz22xmkyPap3htOWZ2e0T7QDP7wMy2mtkzZpbstad473O89QPqE7NIYzlzWDdSk3yEp82f/XA763bsj2lMIo2lvj2PxcAo59xoYAtwB4CZjQCmASOBKcBDZuY3Mz/wIHAeMAK4wtsW4PfAPc65IUAh8H2v/ftAoXNuMHCPt51I3Ak/L/3WycP4zYWjAOPbD/+beR98SnOtapSWq17Jwzn3hnOu1HubBfTxlqcC851zRc65T4Ac4CTvleOcy3XOFQPzgalmZsBZwHPe/k8BF0Yc6ylv+Tng6972InFnfP9OzJo0mKsy+vPqzaeRkZ7G/7y4lv96ZjVfFZXWfACRBNGQcx7XAK95y72B7RHr8r22qtrTgH0RiSjcXuFY3vr93vZHMbOZZpZtZtm7d++u9y8kUh9pbVP423dP5NZzh7Lwo51c8MBSNn92INZhiTSIGqutzGwJ0CPKqjudcwu8be4ESoG54d2ibO+InqxcNdtXd6yjG53LBDIhdJFgtG1EmpLPZ9x41hDG9e/EzU+vZuqDS7n29HRaJflVzisJrcbk4Zw7u7r1ZjYD+CbwdXdkYDcf6BuxWR9gp7ccrX0P0NHMAl7vInL78LHyzSwAdAB0TwhJKKcM6sKiW07je08u54G3cgBIDfiYe53KeSUx1bfaagpwG3CBc+5gxKqFwDSvUmogMARYDnwIDPEqq5IJTaov9JLO28Al3v4zgAURx5rhLV8CvOU0+ygJqFu7VKaM6lnelT5cGuS1tbtiGpPIsarvnMcDQDtgsZmtNrNHAJxz64FngQ3A/wGznHNlXq/iRuB1YCPwrLcthJLQj80sh9CcxuNe++NAmtf+Y6C8vFck0ZwyqAspST7Ct8H6R1YeL3+0s/qdROKQbowo0sTCz0sf3K0tme/msiKvkKsn9ufObwwnJeCPdXjSwukZ5iJxKvJ56Wcd140/vr6ZzHdzWfXpPh66ahx9O7eOcYQiNdPtSURiKMnv43/OH07md8aTV/AV589+jzfWfxbrsERqpJ6HSBw4d2QPXu3ZnlnzVjLz7yv41uieDO3RjlMGdVE1lsQl9TxE4kTfzq355/UTOf/4Hry8Zhd/fmMLVz6axYq8wliHJnIUJQ+ROJIS8DOyV4fyct6i0iDzlufFNCaRaJQ8ROJMRnpaeTmvAc+v2MGfXt9MaVkw1qGJlNOch0icCd+dNyu3gHH9OrJg9U4eeDuH7Ly9zJ52At3ap8Y6RBFd5yGSCJ5fkc/PXlpHm5QAs6eN5ZTBXWIdkjRTtb3OQ8NWIgng4vF9WHDjqXRsncT0xz9g9ptbCQab5xc/SQxKHiIJYmj3diyYdSoXjOnFXxZv4aKH3uePr29SNZbEhJKHSAJpkxLgnsvHcv2Z6XyUv58H3/6YKzJVzitNT8lDJMGYGe1Sk8pvrlhcFmT2m1v1qFtpUkoeIgkoIz2N5IAPv4HP4J0tu7luTjb7DhbHOjRpIVRtJZKgwnfnzRjYmbU79vObRRvp1i6VB68ax9i+HWMdniSo2lZbKXmINBOrt+9j1tyVfHHgMHeeP5wZpwzALNpTnEWqplJdkRZmbN+OvHrzaZwxpCt3v7yBWfNW8uXhkliHJc2UrjAXaUY6tk7m0asn8Oh7ufzh9c2szHuHySN7cMHY3ro7rzQo9TxEmhmfz/jBmYP43wtG8tmXRTy1LI9pmctYsW1vrEOTZkTJQ6SZ2n+opLyct6TM8fMF6/mqqDS2QUmzUa/kYWZ/NLNNZrbGzF40s45e+wAzO2Rmq73XIxH7jDeztWaWY2azzZvRM7POZrbYzLZ6Pzt57eZtl+N9zrj6xCzSUkSW8wZ8xoZdXzL1wffZ+vmBWIcmzUB9ex6LgVHOudHAFuCOiHUfO+fGeq/rI9ofBmYCQ7zXFK/9duBN59wQ4E3vPcB5EdvO9PYXkRqE787743OH8cwPJjL32pPZd7CYCx54nxdW5sc6PElw9Uoezrk3nHPhfnAW0Ke67c2sJ9DeObfMhWqE5wAXequnAk95y09Vap/jQrKAjt5xRKQG4/t3YtakwYzv34lTB3dh0c2nM7pPB3787Efc/vwaDpeUxTpESVANOedxDfBaxPuBZrbKzN4xs9O9tt5A5FeefK8NoLtzbheA97NbxD7bq9inAjObaWbZZpa9e/fu+v02Is1Qt/apzL32ZGZNGsT8D7dz4YPv8/JHO3jw7RzdH0vqpMZSXTNbAvSIsupO59wCb5s7gVJgrrduF9DPOVdgZuOBl8xsJBDtiqWarlKs9T7OuUwgE0IXCdZwXJEWKeD38dPJxzFhQGdumreSm55ejQEpST7mXpuhkl6plRqTh3Pu7OrWm9kM4JvA172hKJxzRUCRt7zCzD4GhhLqNUQObfUBdnrLn5tZT+fcLm9Y6guvPR/oW8U+InKMJg3rxvSM/jzyTi4OKCoJ8n7ObiUPqZX6VltNAW4DLnDOHYxo72pmfm85ndBkd643HHXAzDK8KqurgQXebguBGd7yjErtV3tVVxnA/vDwlojUzzkjepCa5MMIdecXrN7J9r0Ha9pNpH73tjKzHCAFKPCaspxz15vZxcAvCQ1llQG/cM697O0zAfgb0IrQHMlNzjlnZmmYIMqjAAAPYElEQVTAs0A/4FPgUufcXi/JPECoKusg8D3nXI03rdK9rURqJ3yDRTN4+F8fY8CfLxvLOSO6xzo0iQHdGFHJQ6TOPi04yKx5K1m7Yz8zz0jnp5OHkeTXtcQtiW6MKCJ11i+tNf+8fiLfyehP5ru5TMvMYtf+Q7EOS+KQeh4iEtXCj3Zyx/NrSA74mDVpMEWlQTLS0zSh3szVtuehu+qKSFQXjOnFyF7tuebJD/n1qxtVzisVaNhKRKo0qGtbvj0udE2uAw6XBHlz4+exDUrigpKHiFTrtCFdy8t5AeZ+8Cn//nhPTGOS2FPyEJFqhW+weOvkYdxz2Vi6tE1m+mMfcP+bWwkGm+ecqdRMcx4iUqPx/TuVz3OcO7I7//PiWv68eAsf5hVy7+Vj6dwmOcYRSlNTz0NE6qRNSoB7Lx/Lby4aRVZuAeff9x7Zekphi6Oeh4jUmZlx1cn9GdOnI7PmreTyzCyuOrkf3dunkJHeRdVYLYCSh4gcs1G9O/DyTacxc042c5blAZASyGHedSrnbe40bCUi9dI+NYnTh3Qpr8YqKg3y0qodMY1JGp+Sh4jUW0Z6F1KSfPgs9ACeecvzeOrf22iud7AQDVuJSAMIl/Nm5RYwqld75izL4xcL17P8k7387uLjaZeaFOsQpYEpeYhIg4gs5z19SFcy38vlj69vZsOuL3nwynGM6NU+xhFKQ9KwlYg0OJ/PuP7MQTx9XQYHi0u56KH3mb/8Uw1jNSO6q66INKo9/yniR/NXszRnD2cO7cIJ/Tpx+pCuqsaKU3qeh4jEhS5tU3jqmpO4fEJf3tmyh3uXbOXKR7NYkVcY69CkHpQ8RKTR+X1Gv7TW+Lx63qLSIE8szY1tUFIvSh4i0iQy0tNIDvjwG/gMXl37GXe8sIbDJWWxDk2OQb2Th5n9yszWmNlqM3vDzHp57WZms80sx1s/LmKfGWa21XvNiGgfb2ZrvX1mm5l57Z3NbLG3/WIz02CpSIIJl/P++NxhzJ+ZwQ1fG8TTy7dz0UP/5pM9X8U6PKmjek+Ym1l759yX3vLNwAjn3PVmdj5wE3A+cDJwn3PuZDPrDGQDEwg9X2YFMN45V2hmy4FbgCxgETDbOfeamf0B2Ouc+52Z3Q50cs7dVl1cmjAXiX9vb/qC/3p2NaVljt9fPJpvjO4Z65BavCabMA8nDk8bQgkBYCowx4VkAR3NrCcwGVjsnNvrnCsEFgNTvHXtnXPLXCijzQEujDjWU97yUxHtIpLAJh3XjVdvPp0h3dsya95K7l64nqJSDWMlgga5SNDMfgNcDewHJnnNvYHtEZvle23VtedHaQfo7pzbBeCc22Vm3RoibhGJvd4dW/HMzIn8/v828fjST1i6dTdnDe/O5JE9VM4bx2rV8zCzJWa2LsprKoBz7k7nXF9gLnBjeLcoh3LH0F5rZjbTzLLNLHv37t112VVEYig54OPn3xzBf08eRs7ur8h8N5dpmctUzhvHapU8nHNnO+dGRXktqLTpPOBibzkf6Buxrg+ws4b2PlHaAT73hrXwfn5RRZyZzrkJzrkJXbt2rc2vJiJxxEF5OW9JmeP3r22kpCwY05gkuoaothoS8fYCYJO3vBC42qu6ygD2e0NPrwPnmlknr2rqXOB1b90BM8vwqqyuBhZEHCtclTUjol1EmpHIcl6/z1i+rZBpmVns2n8o1qFJJQ1RbfU8MAwIAnnA9c65HV4CeACYAhwEvuecy/b2uQb4H+8Qv3HOPem1TwD+BrQCXgNucs45M0sDngX6AZ8Clzrnqn3upaqtRBLTirxCsnILyEhPY8e+Q9zx/BpSkvzce/lYzhiqEYXGVttqK93bSkTi2se7/8MN/1jJli8OcNOkwdxy9lD8vmhTpNIQaps8dEt2EYlrg7q25aVZp3LXgnXMfiuHD7cV8v3TB7L5swNkpKepIitGlDxEJO61Svbzx0vHcOLAzvzsxbUsyy3AZ6EqrbnX6nnpsaB7W4lIwrhsQl+uPLk/AEEHRSVBln28J8ZRtUxKHiKSUL41phepgdCpywFLNn7B3q+KYxtUC6TkISIJZXz/Tsy9LoOfTh7KD85IZ8POL/nG7PfI3lZtAaY0MFVbiUhCW7djPzfMXcmOfYe4bcowrjs9He+G3HIM9CRBEWkRRvXuwCs3n8Y5w7vz20WbuG7OCvYfLIl1WM2eeh4i0iw453jy/W38dtFGenRI5aazBrPnP8Uq560jXechIi2KmXHNaQMZ268jM+dkc9vzazEgJUnlvI1Bw1Yi0qyM69eJaSf2A0LVWIdLgry7Jeq9VKUelDxEpNmZdFw3UpN85c95eDY7nw07v6x2H6kbJQ8RaXbCz0u/dfIwfj11FEHnuOih95m//FOa6zxvU9Och4g0S+P7dyqf55hyfA9+NH81t7+wluWf7OXXF42idbJOf/WhnoeINHtd2qbw1DUn8aOzh/Di6h1MfeB9tn5+INZhJTSV6opIi7J06x5umb+Kg8VlzDxjIMkBv8p5I+h5HkoeIlKFz788zHefXM7GXaHeR2rAx9zrVM4LusJcRKRK3duncv7xPcursQ6XBvm/dbtiGlOiUfIQkRbplEFdSEnyEX4o4T+yPmXRWiWQ2lLyEJEWKVzO+5Nzh/HX74xnWI923DB3JXcvXE9xaTDW4cW9eiUPM/uVma0xs9Vm9oaZ9fLav2Zm+7321WZ2V8Q+U8xss5nlmNntEe0DzewDM9tqZs+YWbLXnuK9z/HWD6hPzCIiYeP7d2LWpMFMHtmDZ38wkWtOHcjf/r2NS/+6jO17D8Y6vLhW357HH51zo51zY4FXgLsi1r3nnBvrvX4JYGZ+4EHgPGAEcIWZjfC2/z1wj3NuCFAIfN9r/z5Q6JwbDNzjbSci0qCSAz7u+tYIHpk+jtwv/sM371/Kmxs/j3VYcateycM5F3m9fxtCt5KpzklAjnMu1zlXDMwHplro5vtnAc952z0FXOgtT/Xe463/uulm/SLSSKaM6skrN59Gn06t+P5T2dwyfxX3v7WVFXmFsQ4trtR7zsPMfmNm24GrqNjzmGhmH5nZa2Y20mvrDWyP2Cbfa0sD9jnnSiu1V9jHW7/f215EpFH0T2vD8z88hXNHdGfB6p38+Y0tXPVolhJIhBqTh5ktMbN1UV5TAZxzdzrn+gJzgRu93VYC/Z1zY4D7gZfCh4vyEa6a9ur2iRbrTDPLNrPs3bt31/SriYhUKTXJz5i+HctPQCVlQbJyC2IaUzypMXk45852zo2K8lpQadN5wMXePl865/7jLS8CksysC6EeRd+IffoAO4E9QEczC1RqJ3Ifb30HIOrDip1zmc65Cc65CV27dq3xlxcRqU5GehopST78BkkBHxnpGvQIq9edwcxsiHNuq/f2AmCT194D+Nw558zsJEJJqgDYBwwxs4HADmAacKW33dvAJYTmQWYA4eS00Hu/zFv/lmuul8WLSFwJl/Nm5RboFiaV1Pe2kr8zs2FAEMgDrvfaLwF+aGalwCFgmnfCLzWzG4HXAT/whHNuvbfPbcB8M/s1sAp43Gt/HPi7meUQ6nFMq2fMIiK1Fnl3XjlC97YSEZFyureViIg0GiUPERGpMyUPERGpMyUPERGpMyUPERGps2ZbbWVmuwmVDx+LLoQuXEwUiRRvIsUKiRVvIsUKiRVvIsUK9Yu3v3Ouxqusm23yqA8zy65NqVq8SKR4EylWSKx4EylWSKx4EylWaJp4NWwlIiJ1puQhIiJ1puQRXWasA6ijRIo3kWKFxIo3kWKFxIo3kWKFJohXcx4iIlJn6nmIiEidKXlUYmZTzGyzmeWY2e1xEM8TZvaFma2LaOtsZovNbKv3s5PXbmY224t9jZmNi0G8fc3sbTPbaGbrzeyWeI3ZzFLNbLn3xMv1Zva/XvtAM/vAi/UZM0v22lO89zne+gFNFWtEzH4zW2VmryRArNvMbK2ZrTazbK8t7v4OvM/vaGbPmdkm7293YhzHOsz7Nw2/vjSzHzV5vM45vbwXodvEfwykA8nAR8CIGMd0BjAOWBfR9gfgdm/5duD33vL5wGuEnr6YAXwQg3h7AuO85XbAFmBEPMbsfWZbbzkJ+MCL4VlCjxEAeAT4obd8A/CItzwNeCYG/74/JvTgtVe89/Ec6zagS6W2uPs78D7/KeBabzkZ6BivsVaK2w98BvRv6nhj8gvH6wuYCLwe8f4O4I44iGtApeSxGejpLfcENnvLfwWuiLZdDGNfAJwT7zEDrQk9PvlkQhdXBSr/TRB6Ds1EbzngbWdNGGMf4E3gLOAV72QQl7F6nxstecTd3wHQHvik8r9PPMYaJfZzgfdjEa+GrSrqDWyPeJ/vtcWb7s65XQDez25ee1zF7w2VnEDoG31cxuwNA60GvgAWE+p57nPOlUaJpzxWb/1+oCmfS3ov8N+EHr6G99nxGiuAA94wsxVmNtNri8e/g3RgN/CkNyT4mJm1idNYK5sGPO0tN2m8Sh4VWZS2RCpHi5v4zawt8DzwI+fcl9VtGqWtyWJ2zpU558YS+lZ/EjC8mnhiFquZfRP4wjm3IrK5mnji4W/hVOfcOOA8YJaZnVHNtrGMN0BoaPhh59wJwFeEhn2qEg//tnjzWxcA/6xp0yht9Y5XyaOifKBvxPs+wM4YxVKdz82sJ4D38wuvPS7iN7MkQoljrnPuBa85rmN2zu0D/kVoTLijmYUf0RwZT3ms3voOhB6N3BROBS4ws23AfEJDV/fGaawAOOd2ej+/AF4klJzj8e8gH8h3zn3gvX+OUDKJx1gjnQesdM597r1v0niVPCr6EBjiVbAkE+oSLoxxTNEsBGZ4yzMIzSuE26/2qisygP3hbmxTMTMj9Nz5jc65v0SsiruYzayrmXX0llsBZwMbgbeBS6qINfw7XAK85bxB5MbmnLvDOdfHOTeA0N/lW865q+IxVgAza2Nm7cLLhMbm1xGHfwfOuc+A7WY2zGv6OrAhHmOt5AqODFmF42q6eGMxyRPPL0KVCVsIjX3fGQfxPA3sAkoIfYP4PqGx6zeBrd7Pzt62Bjzoxb4WmBCDeE8j1CVeA6z2XufHY8zAaGCVF+s64C6vPR1YDuQQGhJI8dpTvfc53vr0GP1NfI0j1VZxGasX10fea334/1I8/h14nz8WyPb+Fl4COsVrrF4MrYECoENEW5PGqyvMRUSkzjRsJSIidabkISIidabkISIidabkISIidabkISIidabkISIidabkISIidabkISIidfb/ARp3t0E+n1mkAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAD8CAYAAACPWyg8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VdW58PHfczKShMwDQyAhEGYRSarBWVHE6hVbsXW4lbcOXFut9qpv1Q5eb7X37XRrtXUoVVu1WJzFOiEgigMRE0EghCEEAiEkZCIJQ8bzvH+cnXACYQgZTk7yfD+f8zl7r732OWtrkoe11rPXFlXFGGOM6QyXrxtgjDHG/1jwMMYY02kWPIwxxnSaBQ9jjDGdZsHDGGNMp1nwMMYY02kWPIwxxnSaBQ9jjDGdZsHDGGNMpwX6ugE9JT4+XlNTU33dDGOM8Su5ubkVqppwvHr9NnikpqaSk5Pj62YYY4xfEZGiE6lnw1bGGGM6zYKHMcaYTrPgYYwxptP67ZyHMcZ0pKmpieLiYurr633dFJ8KDQ0lOTmZoKCgkzrfgocxZkApLi5m8ODBpKamIiK+bo5PqCqVlZUUFxczatSok/oMvxm2EpFZIrJJRApE5D5ft8cY45/q6+uJi4sbsIEDQESIi4vrUu/LL4KHiAQAjwOXAhOBa0VkYk98V25RNY8vLyC3qLonPt4Y0wcM5MDRqqv/Dfxl2Op0oEBVCwFEZCEwG9jQnV+SW1TNtX/NpqnZTXCgixdvPoOM1Nju/ApjjOkX/KLnAQwHdnrtFztl7YjIPBHJEZGc8vLyTn9JdmElTc1uFGhodnPT8zk8/PYGVm2rosVtz3o3xnRdZWUlU6dOZerUqQwZMoThw4e37Tc2NrbVe+WVV5gwYQIXXHABOTk53HHHHSf9nampqVRUVHRH89v4S8+jo/7VEX/NVXU+MB8gMzOz03/ts9LiCAly0dTsxuUSUuPCeX5lEU9/uo3Y8GBmjE/k4olJnJOewKDggM5fhTFmwIuLi2PNmjUAPPjgg0RERHDPPfccUe+ZZ57hiSee4IILLgAgMzOzV9t5PP4SPIqBEV77yUBJd39JRkoMC27OIruwkqy0ODJSYqirb+LjzeUs2VDG+3mlvJJbTGiQi3PSE5g5MYkZE5KIDQ/u7qYYY/qQ3KLqdn8Xetovf/lLPv30U7Zt28YVV1zBZZddxu9//3vefvttHnzwQXbs2EFhYSE7duzgxz/+cVuv5Morr2Tnzp3U19dz5513Mm/evB5ro78Ejy+BdBEZBewCrgGu64kvykiJaffDMTg0iMunDOPyKcNobHazalsVSzaU8sGGMpZsKMMlkJkay8yJSVw8MYmUuPCeaJYxpgf897/y2FBSe8w6dfVNbCytw63gEhg/ZDCDQ49+b8TEYZH8179N6lK7HnjgAT788EN+//vfk5mZyUcffdTu+MaNG1m+fDl1dXWMGzeOH/zgBwQFBfHss88SGxvLwYMH+cY3vsFVV11FXFxcl9pyNH4RPFS1WURuBxYDAcCzqprX2+0IDnRxdno8Z6fH8+AVk8grqeWDPE8gefidfB5+J59xSYOZOckTSE4ZHmVZHcb4udr6ZlqnPN3q2T9W8OgNl112GSEhIYSEhJCYmEhZWRnJyck89thjvPHGGwDs3LmTLVu2DOzgAaCq7wLv+rodrUSEycOjmDw8irtmjmNH5QGW5JfxQV4pjy8v4E8fFjA0KpSLJiQxc1ISwQEucoqqe63ba4w5vhPpIeQWVXP9054szKBAF49ec5rPf4dDQkLatgMCAmhubuajjz5i6dKlrFy5krCwMM4///wevYveb4JHXzcyLoybzh7FTWePomp/Ix9u3MOSDaW8mlvMC9mHVjgOdAkPXTmZq6YlExzoL8luxgxcHc2F9kU1NTXExMQQFhbGxo0byc7O7tHvs+DRA2LDg5mTkcycjGTqm1r46evreH31LgCa3cr9r6/jv97KY/KwSE4bGcO0kTGcNjKaoVGhNsxlTB90+FxoXzRr1iyeeuoppkyZwrhx48jKyurR7xPV/nn/QmZmpvaVh0F5d3sDA1zcOWMMew82s3pHNWuLa2hodgOQFBnSFkhOGxnDKcOjCA2ylGBjulN+fj4TJkzwdTP6hI7+W4hIrqoeNy/Yeh694Fjd3qYWN/m7a1m9Yy9f7ahm9Y69vLe+FPAMcU0cFslpI6KZlhLDaSNiGBE7yHonxhifs+DRS47W7Q0KcDElOZopydHMPTMVgIp9DazesZfVTjB5JbeY51Z65k3iI4KZOqK1dxLNqcnRhIfY/0ZjTO+yvzp9UHxECBc7940ANLe42Vy2j9U7q/mqaC+rd1azNL8M8OSdjxsSyWkjo9uGvNLiw613YswxqOqA/x3p6pSFBQ8/EBjgYuKwSCYOi+T6M1IA2HugkdU797b1UP71dQkvfrEDgKhBQZ6eidNDEYG1xTV9OlPEmN4SGhpKZWXlgF6WvfV5HqGhoSf9GTZh3k+43crW8n1t8yard+xl8546vP/3ugT+7dRhZKXFMSo+nLSEcBIiQgbsL5AZmOxJgh5He5LgiU6YW/Dox2rrm3j47Q28klPctopkgEvarRAcERLIqPjwtmAyKj6ctPgIRiWEE2FzKcYMOJZtZYgMDeK73xjJW1+XtN0d+8JNZzAkMpRtFfvZVrGfwvJ9FFbs56sd1fxrbUm7nkri4JAOg8qImDC7wdGYAc56HgPAia4IWt/UQlHlAbZVeALKtvL9nveK/VTtP/ScgQCXMCJmEGkJEe16LWnxESRF2jCYMf7Meh6mzYneHRsaFMC4IYMZN2TwEcf2HmhsCyitvZat5fv4fGsF9U3utnphwQGkxrUGk3BGJYQzKt4TZKIG+XYxOWNM97HgYU5IdFgw00YGM21k+yDkdiultfXthsC2Vexn3a4a3l23G+8HMMZHBHv1VCKcobBwRsaFERJod9Ib408seJgucbmEYdGDGBY9iLPGxLc71tDcws6qAxSWt86veN4/3FjOyznFhz5DIDkm7IghsP2NzWwpq2P66HhLMTamj7HgYXpMSGAAYxIHMybxyGGwmoNNbG+dtHd6Ldsq9vPl9ioONLa0qytsZsLQwYwfEtkWqIZFh5IcM4ihUYPsDntjfMB+64xPRA0K4tQR0Zw6Irpduaqyp66BPyzZxMtfelKMFU+w+WJbFaW19e1SjQGiw4IYFuUJKsOjQ70CzCCGRw8icXAILpdN4hvTnSx4mD5FREiKDOU7mSNZtOZQivFj104jIyWG5hY3e+oaKNl7kF3Oq2TvQUr21lNcfYAvtlVSV9/c7jODAoQhUaEMi/IEE+/eS+u+9V6M6RxL1TV91ommGB+utr6J3Xvr2wJMyWFBxnovxhydpeoav3eyD+CJDA0ickhQhynHwBG9lxIn0JTsPWi9F2NOUI/9tIvIg8AtQLlT9FPnOeSIyP3ATUALcIeqLnbKZwGPAgHA06r6a6d8FLAQiAW+Ar6nqofuWjOmEwIDXG1//I/2z6uOei+tPZeTnXtJGBxCgPVeTD/R0/9UekRVf+9dICITgWuAScAwYKmIjHUOPw5cDBQDX4rIW6q6AfiN81kLReQpPIHnyR5uuxnAutp7WbWtktoT7L0caGymZO9BZkxIOiLd2Zi+yhf97NnAQlVtALaJSAFwunOsQFULAURkITBbRPKBC4HrnDrPAQ9iwcP40In0Xurqm9qCyon0Xp79bDsRIQGMjA1v670Mjzms9xJhcy+mb+jp4HG7iNwA5AB3q2o1MBzI9qpT7JQB7Dys/AwgDtirqs0d1G9HROYB8wBGjhzZXddgzEkZHBrEuGP0Xlrcyu8Wb2T+ikLcCgKMSRxMXHjwMedehkZ55lpaA8rhczBhwTb3Ynpel37KRGQpMKSDQz/D0zN4CE+a/kPA/wI34vkdOZwCHS3Tqseof2Sh6nxgPniyrY7TfGN8KsAlXDxxCH//fHtbSvIvLp/YLkmgde5l194D7PIaGttVfZDsrZWU1tZz2NQLMWFBnh5L1KEey6EeTCjx4dZ7MV3XpeChqhedSD0R+SvwtrNbDIzwOpwMlDjbHZVXANEiEuj0PrzrG+PXMlJiWHBz1lFTkk9k7qWsde6luv3wWFHlAT7fWsm+hva9l+AAF0O9ssTahsiiw9p6NKFBttaYObaezLYaqqq7nd1vAeud7beAF0XkD3gmzNOBVXh6GOlOZtUuPJPq16mqishyYA6ejKu5wKKearcxve1kU5LBM/fSOnT1jdQjj6sqtfXNh3osh03wf7qlgrK6eg6/3SsuPLit1+Kdktzag4kLD7al9we4nhwc/a2ITMUzxLQd+A8AVc0TkZeBDUAzcJuqtgCIyO3AYjypus+qap7zWfcCC0XkYWA18EwPttuYfkNEiBoURNSgICYMjeywTlOLm9IaZ0isprUH49nfWr6PFVvKj1hvLCTQdcR9Lt5zMEOiQo/ovZzsTZ+mb7I7zI0xx6Sq1Bxsatdj2XVYBtmeuoYjei/xESFt97sEuoT31pfiViUowMWLN59BRmqsby7IHJPdYW6M6RYiQnRYMNFhwUwaFtVhnYbmFspqGtoFlNYAs7msjqLKAzQ7M/sNzW6uf+YLTk2OZtyQwYxNGtz2bg8M8x8WPIwxXRYSGMDIuDBGxoV1eDx3exXXPf0FTS1uXCKcMyaByv0NvP7VrnYT+kOjQtsFk3FJg0lPirAJ/D7IgocxpsdlpMby4i1HZpWpKiU19WwurWNTWR2bS+vYWFrHysJKGps9jzcWgdS4cMYmRTAuaTBjh3iCSmp8OEEBHWX4m95gcx7GmD6nucVNUdWBQ0GlzBNUtlfsb7uvJTjARVpCeLteyrghgxkePcjuY+kCm/MwxvitwAAXoxMiGJ0QwaWnDG0rr29qYWv5PjaX1bGp1POes72aRWsO3foVHhxAuhNMWnspY4dEkBARYunF3ciChzHGb4QGBTBpWNQRE/d19U1sLmsNKp6eytL8Ml7KObTiUWx48BFDX+nOJL2lEXeeBQ9jjN8bHBrU4c2WFfsa2g19bSqt47XDJunjw4OpOtCIKoQEuVhwc5YFkBNgwcMY02/FR4QQPyaEM72WulfVthTiTaX7+NfXJVTs9zweqL7JzSdbyi14nABLVTDGDCgiQnJMGBeOT+IH54/moSsnExrkaluB9fXcXWyv2O/TNvoDCx7GmAGtdXHKey4Zxy8un0BtQxP/9udPWZZf5uum9WmWqmuMMV52Vh3g1n/kkldSy50z0rlzRvqASv090VRd63kYY4yXEbFhvPaDM/n2tOE8umwLNz33JTUHmnzdrD7HJsyNMeYwoUEB/O/Vp3LaiGj++18buOLxT7lzRjq7a+otnddhwcMYYzogInxveioRoYH850tfc9fLX+MSCA60dF6w4GGMMe2oKpvL9vH++lLezyslf3dt2zG3QlOzm+zCSgsevm6AMcb4mqqytriG9/NKWby+lMKK/YhAZkoMP79sAsOjB/GfL69pe9Z8Vlqcr5vscxY8jDEDUotbydlexXvrS/kgr5SSmnoCXcL00XHcePYoZk5KInFwaFv9xMhQW8LEiwUPY8yA0djs5vOtFSzOK+WDvDIq9zcSEuji3LEJ3D1zHDMmJBIdFtzhuV151nx/ZMHDGNOvHWxs4ePN5SzOK2Vpfhl19c2EBwdw4YQkZk0awvnjEggPsT+FndWl/2IicjXwIDABOF1Vc7yO3Q/cBLQAd6jqYqd8FvAoEAA8raq/dspHAQuBWOAr4Huq2igiIcDzQAZQCXxXVbd3pd3GmP4rt6iajzfvQRA2ldbx0eY91De5iQ4LYtakIVx6yhDOHB1vTyfsoq6G2/XAt4G/eBeKyETgGmASMAxYKiJjncOPAxcDxcCXIvKWqm4AfgM8oqoLReQpPIHnSee9WlXHiMg1Tr3vdrHdxph+omJfA/m7a8nfXcunBRV8srmC1nUzYsKC+E7mCGZNGsLpo2IJtCcPdpsuBQ9VzQc6esDKbGChqjYA20SkADjdOVagqoXOeQuB2SKSD1wIXOfUeQ5Pj+ZJ57MedMpfBf4sIqL9dV0VY0yHmlrcFJbvbwsUG3bXsrG0jvK6hrY6ESGBbYHDJXDT2aO4/cJ03zS4n+upgb7hQLbXfrFTBrDzsPIzgDhgr6o2d1B/eOs5qtosIjVO/Yqeaboxxteq9jey0QkQ+bvryN9dS8GefTS2eJ5rHhzgIj0pgvPGJjB+yGAmDo1k/NBItlXs5/qns9tSaqePjj/ON5mTddzgISJLgSEdHPqZqi462mkdlCkdr6Wlx6h/rM868ktF5gHzAEaOHHmUphlj+ormFjfbKvaTX1rX1qPI311LWe2h3kTC4BAmDI3knLHxTBgSyYShkaQlhBPUwRBUbHgwC27OspTaXnDc4KGqF53E5xYDI7z2k4HWhwx3VF4BRItIoNP78K7f+lnFIhIIRAFVR2nrfGA+eFbVPYl2G2N6SM2BJqcnUcvGUk+PYnNZHQ3Nnt5EUIAwOiGCs0bHM2FoJOOHDmbC0EjiI0I69T2WUts7emrY6i3gRRH5A54J83RgFZ5eRLqTWbULz6T6daqqIrIcmIMn42ousMjrs+YCK53jH9p8hzF9V4tb2V55aG5iozPsVFJT31YnLjyYCUMjuWF6CuOd3sSYxAiCA21C2190NVX3W8CfgATgHRFZo6qXqGqeiLwMbACagdtUtcU553ZgMZ5U3WdVNc/5uHuBhSLyMLAaeMYpfwZ4wZl0r8ITcIwxfcAnW8p5d10p4cEB7GtoJr+0jk2ltdQ3eXoTAS5hdEI43xgV6+lNOPMTCYNDOkq0MX7EHgZljDkuVaW8roG8klrySmrIK6ll9Y5qSmu9M50COGV4NBOGRjLBGXIakxhh91P4mRN9GJTdVmmMacftVoqqDrQFibySWjaU1FCxr7GtTkpcGJGDgiirbfBkwgj84PzR3HaBpcUOFBY8jBnAGppb2FK2jw1ePYr83bXsb2wBINAlpCcN5ryxiUwaFsmkYZFMGBZJZGgQuUXV7dJis9IsLXYgseBhzABRV9/kBAnP/RN5JbUU7KmjqcUzdB0eHMCEoZHMyUhm4rBIJg2LIj0pgpDAjoedMlJiLC12ALPgYUw/tKe2vt38xIbdtRRVHmg7Hh8RzMRhUZw/LoGJQz09itS4cFyuzk1iW1rswGXBwxg/kltU3e5f+h3PT9RSse/QRPbI2DAmDYvkaq8eRaJlO5kusuBhjJ/ILarmur9m09jsxuUS0hMj2Fl1oN38xJhEz5Idh89PGNPdLHgY04e53cqa4r18mL+Hl77c0XY3dotbqa1v4qqMZCdQHHt+wpjuZsHDmD5mX0Mzn24pZ2n+Hj7atIeKfY0EuISxSRFUH2jCrUpwoIs/XTvN5huMz1jwMKYP2Fl1gGX5ZSzbuIcvCqtobHETGRrI+eMSmTEhkfPGJhAdFnzEnIcxvmLBwxgfaHErq3dUs2zjHpbll7G5bB8AaQnhzD0zhRkTkshIiTli5VjLbjJ9hQUPY3pJbX0Tn2yuYFl+Gcs37aH6QBOBLuH0UbF8J3MEMyYkMSo+3NfNNOaEWPAwpge0Di+lxoVRWtvAhxvL+KKwima3Eh0WxAXOcNQ56QlEDbJsKON/LHgY080+L6jghmdX0ew+tOhoemIEN5+TxowJiUwbGUNAJ2/GM6avseBhTDcqrannxy+taQscAtxyTho/vWyCbxtmTDezJ68Y003ySmq48vHPqD3YRHCAECAQEuTikskdPcXZGP9mPQ9jusGHG8u4/cXVRA0K4vUfnsXBphZLqTX9mgUPY7rob59t46G3NzBpWBTPzM0kMTIUwIKG6dcseBhzkppb3Dz09gaeW1nEzIlJ/PGaqYQF26+UGRjsJ92Yk/DJlnIeWJTHtor9zDs3jftmje/0cubG+LMuTZiLyNUikicibhHJ9CpPFZGDIrLGeT3ldSxDRNaJSIGIPCbOutAiEisiS0Rki/Me45SLU69ARNaKyLSutNmYrlq5tYIbnlnFtor9BAUIl0waYoHDDDhdzbZaD3wbWNHBsa2qOtV53epV/iQwD0h3XrOc8vuAZaqaDixz9gEu9ao7zznfGJ/JLqyi9Q4Ot1vJLqz0aXuM8YUuBQ9VzVfVTSdaX0SGApGqulJVFXgeuNI5PBt4ztl+7rDy59UjG4h2PscYnzh3bAKtHQ3Ps7vjfNsgY3ygJ+/zGCUiq0XkYxE5xykbDhR71Sl2ygCSVHU3gPOe6HXOzqOcY0yvy0iJ4f5vem76u+PCdMuqMgPScYOHiCwVkfUdvGYf47TdwEhVPQ24C3hRRCLx3HB7OO2grF0TTvQcEZknIjkiklNeXn6cjzXm5H3/zFQSB4fw1Y69vm6KMT5x3GwrVb2osx+qqg1Ag7OdKyJbgbF4eg3JXlWTgRJnu0xEhqrqbmdYao9TXgyMOMo5h3/vfGA+QGZm5vGCkjEnLTDAxbemDeeZT7ZRsa+B+IgQXzfJmF7VI8NWIpIgIgHOdhqeye5CZziqTkSynCyrG4BFzmlvAXOd7bmHld/gZF1lATWtw1vG+NKcack0u5X/++pacouqfd0cY3pVV1N1vyUixcB04B0RWewcOhdYKyJfA68Ct6pqlXPsB8DTQAGwFXjPKf81cLGIbAEudvYB3gUKnfp/BX7YlTYb011q65sRgeUb93DdX7MtgJgBpUs3CarqG8AbHZS/Brx2lHNygMkdlFcCMzooV+C2rrTTmJ6QXViJ4JmAa2h2887aEps8NwOGraprzEnKSosjONDVlra78MudfF5Q4dtGGdNLLHgYc5IyUmJYcHMWd88cx1P/Po3kmEHc8OwqXs0tPv7Jxvg5W9vKmC7ISIlpG6qaPjqeHy7I5Z5XvmZn1QF+fFE6zuo7xvQ71vMwpptEDQrib//ndOZkJPPosi3c/fLXNDa7fd0sY3qE9TyM6UbBgS5+N2cKKbFh/O+SzWwqq2XG+CTOG5dok+mmX7HgYUw3ExF+NCOdZrebR5cVkFdSx1MrCvnnLVkWQEy/YcNWxvSQ4MCAtkysxmY3v31/ow1jmX7DgocxPaQ1lTdAIECEL7ZVceXjn7G5rM7XTTOmy8RzD17/k5mZqTk5Ob5uhhngcouqyS6sJCstjsp9Ddz/+jrqGpr5ySXjuPGsUfYQKdPniEiuqmYer57NeRjTg7xTeQGmpcRw32vrePidfJbl7+H33zmV4dGDfNhCY06ODVsZ04viI0L46w0Z/PaqKawt3susR1bwWm4xudureHx5ga2PZfyG9TyM6WUiwne+MYLpo+O46+U13P3K120T68GBLhbcbFlZpu+znocxPjIiNoyF86Zz3tgE3Apu9WRl2TPRjT+w4GGMDwW4hDtmpBMS6PlVdCt8sqWcPXX1Pm6ZMcdmwcMYH8tIieHFW7K4e+ZY/v2MFL7asZeZj6zgjdXF9NdsSOP/LFXXmD6mYM8+fvLq13y1Yy8Xjk/kV9+azNAoy8gyveNEU3Wt52FMHzMmMYJXbj2TX1w+kc+3VjDzDytYuGqH9UJMn2I9D2P6sKLK/dz72lqyC6s4e0w8158xksKK/WSlxVlGlukRdpOgMf1ASlw4L96cxYurdvDw2xv4tKACAUKCLKXX+JYNWxnTx7lcwr9npTD3rFTA88z0+iY3767b7dN2mYGtS8FDRH4nIhtFZK2IvCEi0V7H7heRAhHZJCKXeJXPcsoKROQ+r/JRIvKFiGwRkZdEJNgpD3H2C5zjqV1pszH+aubEIYQGuWhdDeu5ldt54qMCmltspV7T+7ra81gCTFbVKcBm4H4AEZkIXANMAmYBT4hIgIgEAI8DlwITgWudugC/AR5R1XSgGrjJKb8JqFbVMcAjTj1jBpzWZ6bfc8k4np6byUXjk/jt+5u48onPyCup8XXzzADTpeChqh+oarOzmw0kO9uzgYWq2qCq24AC4HTnVaCqharaCCwEZovnQc8XAq865z8HXOn1Wc85268CM8QeDG0GqIyUGG67YAwXTUjiqe9l8OT10yitaeCKP3/G7xZvpL6pxddNNANEd8553Ai852wPB3Z6HSt2yo5WHgfs9QpEreXtPss5XuPUP4KIzBORHBHJKS8v7/IFGdPXXXrKUJbedS5XTh3O48u3ctljn5BbVOXrZpkB4LjBQ0SWisj6Dl6zver8DGgGFrQWdfBRehLlx/qsIwtV56tqpqpmJiQkHO2SjOlXosOC+d/vnMpzN55OfZObOU+t5MG38visoMJW6jU95ripuqp60bGOi8hc4HJghh66aaQYGOFVLRkocbY7Kq8AokUk0OldeNdv/axiEQkEogD7p5UxhzlvbAKL//Ncfvv+Rv7++Xae+3w7IrZSr+kZXc22mgXcC1yhqge8Dr0FXONkSo0C0oFVwJdAupNZFYxnUv0tJ+gsB+Y4588FFnl91lxnew7wofbXOxuN6aKIkEB+OXsy154+EsWz0GJ9k5uPN+3xddNMP9PVOY8/A4OBJSKyRkSeAlDVPOBlYAPwPnCbqrY4vYrbgcVAPvCyUxc8QeguESnAM6fxjFP+DBDnlN8FtKX3GmM6NicjmdDAQ2m9z2cX8f76Up+2yfQvtjyJMf1U6/PTEweH8LfPtrNhdy3fPGUID14xicTBob5unumjTnR5EgsexgwATS1u5q8o5NFlWxgUFMADl0/k29OGY1nv5nC2qq4xpk1QgIvbLhjDu3ecw5jECO5+5Wvm/u1LiqsPHP9kYzpgPQ9jBhi3W3khu4jfvL8RAa49YyTRg4KYPjreMrKMraprjOmYyyXMPTOVC8cncvuLX/H0J9sACAks4MVbLKXXnBgbtjJmgBoRG8bMSUltGVkNzW4eXbqZJlto0ZwACx7GDGBZafGEBLkIEHAJrNhSwZWPf8b6XbbQojk2m/MwZoBrTenNSoujvK6en7+ZR/WBRm49L40fXZhOaFCAr5toepHNeRhjTkhGSky7eY7pafE8/M4GHl++lffWl/Lbq6aQmRrrwxaavsiGrYwx7USFBfG7q0/l+RtPp6HJzdV/Wcl/LVrPvobm459sBgwbtjLGHNX+hmb6tuG+AAATQklEQVR+t3gTz63czrCoQdx4dir1TW6y0uIsK6ufsjvMLXgY021yi6q445+r2bW3HoCQQJel9fZTdoe5MabbZKTE8p3MEe3Sep//fLsvm2R8zIKHMeaEnJ2eQEiQC5eACCz6uoRbX8hlT229r5tmfMCGrYwxJ6w1rfcbqTHkFu3lkaWbCQ108fPLJ3J1RrIttNgP2JyHBQ9jelxh+T7ue20dq7ZXcU56PP/zrVMYERvm62aZLrA5D2NMj0tLiGDhvCwemj2Jr4qqueSPK/jbZ9tocffPf5SaQ6znYYzpFrv2HuSnr6/j483ljE2K4Oz0eC47ZZhlZPkZ63kYY3rV8OhB/P373+COC8ewuWwfz366ne/+ZSWrtlX6ummmB1jwMMZ0GxEhJCgAlzNv3uxWblvwFeuKbaHF/qZLwUNEficiG0VkrYi8ISLRTnmqiBwUkTXO6ymvczJEZJ2IFIjIY+KkZ4hIrIgsEZEtznuMUy5OvQLne6Z1pc3GmJ6VlRZHcKBnpd6gAKHJrVz5xGf8+r2N1De1+Lp5ppt0teexBJisqlOAzcD9Xse2qupU53WrV/mTwDwg3XnNcsrvA5apajqwzNkHuNSr7jznfGNMH5WREsOCm7O4a+Y4Fs6bzsf3XMCcack89fFWLn30E74otGGs/qBLwUNVP1DV1tXSsoHkY9UXkaFApKquVM9M/fPAlc7h2cBzzvZzh5U/rx7ZQLTzOcaYPiojJYbbLhhDRkoMUWFB/GbOFBbcfAbNbjffnZ/NL95cT119k6+babqgO+c8bgTe89ofJSKrReRjETnHKRsOFHvVKXbKAJJUdTeA857odc7Oo5zTjojME5EcEckpLy/v2tUYY7rVWWPiWfzjc7nxrFH844siLnlkBcs37SG3qJrHlxeQW1Tt6yaaTjju8zxEZCkwpINDP1PVRU6dnwHNwALn2G5gpKpWikgG8KaITAI6uv30eLnCJ3yOqs4H5oMnVfc4n2uM6WVhwYE88G8TufzUodz76lq+/7cvCRBBUYIDXSy42RZb9BfH7Xmo6kWqOrmDV2vgmAtcDlzvDEWhqg2qWuls5wJbgbF4eg3eQ1vJQImzXdY6HOW873HKi4ERRznHGOOHpo2M4e07zmZ6WiwtqrgVGpvdZBdW+Lpp5gR1NdtqFnAvcIWqHvAqTxCRAGc7Dc9kd6EzHFUnIllOltUNwCLntLeAuc723MPKb3CyrrKAmtbhLWOM/woJDOCeS8YTEuj5M+RW+HhTOWW20KJf6NId5iJSAIQArekT2ap6q4hcBfwSz1BWC/Bfqvov55xM4O/AIDxzJD9SVRWROOBlYCSwA7haVaucIPNnPFlZB4Dvq+pxbx23O8yN8Q+5RdV8vrWC8roGXvpyJ8GBLn5+2QTPEvC20GKvs4URLXgY43e2V+zn3tfW8sW2Ks4cHcevvz2FkXG20GJvsuVJjDF+JzU+nH/eksWvvjWZtcU1XPLHFTzzqS202BdZz8MY0yftrjnIz95Yz4cb9zB1RDT/58xUdu09aM9P72E2bGXBwxi/p6q89XUJP39jHXUNLQgQEmQpvT3Jhq2MMX5PRJg9dTg3nJkKeG7wqm9y8+bqXT5tl7HgYYzxAxeOTyLUeX46wD+yi/ifd/M52GgLLfrKce8wN8YYX2tdbDG7sJIpyVG8u66U+SsK+SCvlP/37SlMHx3n6yYOODbnYYzxS59vreD+19dRVHmA684YyX2XjicyNMjXzfJ7NudhjOnXzhwdz/t3nsst54xi4aodzPzDCj7cWObrZg0Y1vMwxvi9NTv3cu+ra9lUVsfsqcO4cuowNuyus7Tek2CpuhY8jBlQGpvdPPnRVh77cDMtbiyt9yTZsJUxZkAJDnRx50Xp3DA9FTiU1rtkQ6lP29VfWfAwxvQrl08ZRmiQq+1BQM9/XsSLX+zAbUucdCsbtjLG9Du5RdVkF1aSGhfGP7J3sLKwkqy0WH797Smkxof7unl9ms15WPAwxuBZ4uSlL3fyq3fyaXK7ufvicdx49igCXLbce0dszsMYY/AscXLN6SNZctd5nD0mgV+9m8+3n/iMjaW1vm6aX7OehzFmwFBV3l67mwffyqPmYBPfmjacETGDOGtMgmVkOWzYyoKHMeYoqvY38uOXVrNis+eZ6cGBLv55i6X0gg1bGWPMUcWGB3PGqLi2jKzGZjf/7918DjQ2+7Rd/sSChzFmQMpKiyMkyEWAQIBLyCmqZtYfP+HzggpfN80vdDl4iMhDIrJWRNaIyAciMswpFxF5TEQKnOPTvM6ZKyJbnNdcr/IMEVnnnPOYiIhTHisiS5z6S0TE+pbGmC5pXan3rpnjePk/prNwXhYugeue/oL7XltLzcEmXzexT+vynIeIRKpqrbN9BzBRVW8VkW8CPwK+CZwBPKqqZ4hILJADZOK5CTQXyFDVahFZBdwJZAPvAo+p6nsi8lugSlV/LSL3ATGqeu+x2mVzHsaYzqpvauGRpZv564pC4iNCePjKycycNMTXzepVvTbn0Ro4HOF4AgLAbOB59cgGokVkKHAJsERVq1S1GlgCzHKORarqSvVEtOeBK70+6zln+zmvcmOM6TahQQHcf+kE3rztLGLDg5n3Qi63v/gVFfsafN20PqdbHgYlIr8CbgBqgAuc4uHATq9qxU7ZscqLOygHSFLV3QCqultEEruj3cYY05EpydH860dn89RHW/nThwV8WlDB3OkpBAe6yEqLt6wsTrDnISJLRWR9B6/ZAKr6M1UdASwAbm89rYOP0pMoP2EiMk9EckQkp7y8vDOnGmNMO0EBLn40I5137jibxMEhPLqsgN8t3sx1f80mt6ja183zuRMKHqp6kapO7uC16LCqLwJXOdvFwAivY8lAyXHKkzsoByhzhrVw3vccpZ3zVTVTVTMTEhJO5NKMMeaY0pMGc8Wpw9r+ddvQ7ObJjwoG/EKL3ZFtle61ewWw0dl+C7jBybrKAmqcoafFwEwRiXGypmYCi51jdSKS5WRZ3QAs8vqs1qysuV7lxhjT46aPjickyIVLwCWwNH8P18zPprB8n6+b5jPdkW31GjAOcANFwK2qussJAH8GZgEHgO+rao5zzo3AT52P+JWq/s0pzwT+DgwC3gN+pKoqInHAy8BIYAdwtapWHatdlm1ljOlOrSv1Zo2KZWv5fh56ZwMNzW7+86Kx3HLOKAID+sdtc7Y8iQUPY0wPKqut5xdvrueDDWVMHh7Jb686lYnDIn3drC6z4GHBwxjTw1SV99aX8sCi9ew90MSt543mrDFxfLVjr98+P92ChwUPY0wvqd7fyEPvbOD1r3YhgIhnsUV/fH66LYxojDG9JCY8mD98ZypzMoajgFuhocnNJ5v77y0DFjyMMaabXHt6CqGBnuenK7BgVRGfbOmfAcSGrYwxphu1ZmVFhQbx7GfbKKzYz9UZyfz8solEhQX5unnHZXMeFjyMMT5W39TCY8u28JcVhcSGB/PQ7MnMmty3F1q0OQ9jjPGx0KAAfjJrPItuO4uEiBBu/UcuP1yQy566el83rcu6ZWFEY4wxRzd5eBSLbj+L+SsKeXTZFj4rqOSByyeSGhdG9rYqv0zrtWErY4zpRQV79nHfa2vJKarG5SyY1ZfSem3Yyhhj+qAxiRG8/B/TmTEhEbceSutdudW/Hn9rwcMYY3qZyyX88PwxhAR6/gQr8M7a3Wz1o4UWbdjKGGN8xJPWW8HBJjcvrCziYFMLP74onVvOSSPIRwstWqquBQ9jjB/ZU1fPg2/l8e66UiYNi+Q3V01h8vCoXm+HzXkYY4wfSRwcyhPXZ/DUv0+jrLaB2Y9/xm/f30h9U4uvm9YhS9U1xpg+ZNbkoUxPi+fhdzbwxEdbeT+vlJvOGsXeg019KqXXhq2MMaaPWrG5nLtfWUN5XSMAoYEuFtzSsym9NmxljDF+7tyxCVx/ekrbfn2zm3+u2uHDFh1iwcMYY/qwc8YmEOo8P12AV3OLuevlNew90OjTdtmchzHG9GEZKTEsuDmL7MJKpo2M5rOCSp76eCsrNpfzy9mT+eYpQ33Sri71PETkIRFZKyJrROQDERnmlJ8vIjVO+RoRecDrnFkisklECkTkPq/yUSLyhYhsEZGXRCTYKQ9x9guc46ldabMxxvibjJQYbrtgDNNHx3PPJeNYdPtZDIkK5YcLvuI/XshhT23vL7TY1WGr36nqFFWdCrwNPOB17BNVneq8fgkgIgHA48ClwETgWhGZ6NT/DfCIqqYD1cBNTvlNQLWqjgEeceoZY8yANWlYFG/+8Czuu3Q8H20q56I/fMzLOTvpzQSoLgUPVa312g3Hc5f9sZwOFKhqoao2AguB2SIiwIXAq06954Arne3Zzj7O8RlOfWOMGbACA1zcet5o3rvzHMYPieQnr67le8+s4r11u3l8eQG5RdU9+/1d/QAR+RVwA1ADXOB1aLqIfA2UAPeoah4wHNjpVacYOAOIA/aqarNX+XBnu+0cVW0WkRqnvn+tImaMMT0gLSGChfOyWLBqB796ewOfFlQgQEhQz67Ue9yeh4gsFZH1HbxmA6jqz1R1BLAAuN057SsgRVVPBf4EvNn6cR18hR6j/FjndNTWeSKSIyI55eX987nBxhhzOJdL+F5WCnPPSgU8fyCbmt1kF1b23Hcer4KqXqSqkzt4LTqs6ovAVc45taq6z9l+FwgSkXg8PYoRXuck4+mZVADRIhJ4WDne5zjHo4Cqo7R1vqpmqmpmQkLCcS/eGGP6k5kThxAa5CJAICjQRVZaXI99V5eGrUQkXVW3OLtXABud8iFAmaqqiJyOJ0hVAnuBdBEZBewCrgGuc+otB+bgmQeZC7QGp7ec/ZXO8Q+1v94Wb4wxXeCd1tvTS5l0dc7j1yIyDnADRcCtTvkc4Aci0gwcBK5x/uA3i8jtwGIgAHjWmQsBuBdYKCIPA6uBZ5zyZ4AXRKQAT4/jmi622Rhj+q2MlJheWf/K1rYyxhjTxta2MsYY02MseBhjjOk0Cx7GGGM6zYKHMcaYTrPgYYwxptP6bbaViJTjSR8+GfH0n+VP7Fr6nv5yHWDX0ld15VpSVPW4d1n32+DRFSKScyKpav7ArqXv6S/XAXYtfVVvXIsNWxljjOk0Cx7GGGM6zYJHx+b7ugHdyK6l7+kv1wF2LX1Vj1+LzXkYY4zpNOt5GGOM6TQLHocRkVkisklECkTkPl+353hE5FkR2SMi673KYkVkiYhscd5jnHIRkceca1srItN81/L2RGSEiCwXkXwRyRORO51yf7yWUBFZJSJfO9fy3075KBH5wrmWl0Qk2CkPcfYLnOOpvmz/4UQkQERWi8jbzr6/Xsd2EVknImtEJMcp87ufLwARiRaRV0Vko/M7M723r8WChxcRCQAeBy4FJgLXishE37bquP4OzDqs7D5gmaqmA8ucffBcV7rzmgc82UttPBHNwN2qOgHIAm5z/tv747U0ABc6T9KcCswSkSzgN8AjzrVUAzc59W8CqlV1DPCIU68vuRPI99r31+sAuEBVp3qlsfrjzxfAo8D7qjoeOBXP/5/evRZVtZfzAqYDi7327wfu93W7TqDdqcB6r/1NwFBneyiwydn+C3BtR/X62gvPw8Au9vdrAcLwPJb5DDw3bQUe/rOG5/k2053tQKee+LrtTnuS8fwhuhB4G89jof3uOpw2bQfiDyvzu58vIBLYdvh/296+Fut5tDcc2Om1X+yU+ZskVd0N4LwnOuV+cX3OcMdpwBf46bU4Qz1rgD3AEmArsFdVm50q3u1tuxbneA3Qc88P7Zw/Aj/B88A38LTLH68DPI/2/kBEckVknlPmjz9faUA58DdnOPFpEQmnl6/Fgkd70kFZf0pH6/PXJyIRwGvAj1W19lhVOyjrM9eiqi2qOhXPv9xPByZ0VM1575PXIiKXA3tUNde7uIOqffo6vJylqtPwDOPcJiLnHqNuX76WQGAa8KSqngbs59AQVUd65FoseLRXDIzw2k8GSnzUlq4oE5GhAM77Hqe8T1+fiAThCRwLVPV1p9gvr6WVqu4FPsIzjxMtIq2PfvZub9u1OMej8Dxy2dfOAq4Qke3AQjxDV3/E/64DAFUtcd73AG/gCer++PNVDBSr6hfO/qt4gkmvXosFj/a+BNKdbJJgPM9Lf8vHbToZbwFzne25eOYPWstvcLIvsoCa1m6ur4mI4Hlefb6q/sHrkD9eS4KIRDvbg4CL8ExoLgfmONUOv5bWa5wDfKjO4LQvqer9qpqsqql4fhc+VNXr8bPrABCRcBEZ3LoNzATW44c/X6paCuwUkXFO0QxgA719Lb6e/OlrL+CbwGY8Y9Q/83V7TqC9/wR2A014/oVxE55x5mXAFuc91qkreLLJtgLrgExft9/rOs7G05VeC6xxXt/002uZAqx2rmU98IBTngasAgqAV4AQpzzU2S9wjqf5+ho6uKbzgbf99TqcNn/tvPJaf7f98efLad9UIMf5GXsTiOnta7E7zI0xxnSaDVsZY4zpNAsexhhjOs2ChzHGmE6z4GGMMabTLHgYY4zpNAsexhhjOs2ChzHGmE6z4GGMMabT/j+tVt+qAXhlTQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "yr2sec = 60*60*24*365.25 # num of seconds in a year\n", "\n", "# Set material properties\n", "rho = 2900 # kg m^-3\n", "Cp = 900 # J kg^-1 K^-1\n", "k = 2.5 # W m^-1 K^-1\n", "\n", "# Set geometry and dimensions\n", "L = 35e3 # m\n", "nx = 30 # -\n", "dx = L / (nx-1) # m\n", "print(\"dx is\", dx, \"m\")\n", "\n", "silltop = 5e3 # m\n", "sillbott = 10e3 # m\n", "\n", "# Set total time and time step\n", "t_total = 0.1e6 * yr2sec # s\n", "dt = 1e4 * yr2sec # s\n", "nt = int(np.floor(t_total / dt)) # -\n", "print(\"Calculating\", nt, \"time steps\")\n", "\n", "# Set boundary temperature values and intrusion temperature\n", "Tsurf = 0.0 # deg C\n", "Tmoho = 600.0 # deg C\n", "Tintrusion = 700.0 # deg C\n", "\n", "# Create arrays to hold temperature fields\n", "Tnew = np.zeros(nx)\n", "Tprev = np.zeros(nx)\n", "\n", "# Create coordinates of the grid points\n", "x = np.linspace(0, L, nx)\n", "\n", "# Generate initial temperature field\n", "for ix in range(nx):\n", " if (x[ix] < sillbott) and (x[ix] > silltop):\n", " Tprev[ix] = Tintrusion\n", " else:\n", " Tprev[ix] = x[ix] * (Tmoho - Tsurf) / L\n", "\n", "# Plot initial temperature field\n", "plt.plot(Tprev, -x, '.-', label='T initial')\n", "plt.legend()\n", "plt.show()\n", "\n", "# Start the loop over time steps\n", "curtime = 0\n", "while curtime < t_total:\n", " curtime = curtime + dt\n", " \n", " # boundary conditions\n", " Tnew[0] = Tsurf\n", " Tnew[nx-1] = Tmoho\n", " \n", " # internal grid points\n", " for ix in range(1, nx-1):\n", " Tnew[ix] = (k * (Tprev[ix-1] - 2.0*Tprev[ix] + Tprev[ix+1]) / dx**2) * dt / (rho * Cp) + Tprev[ix]\n", " \n", " # copy Tnew to Tprev\n", " Tprev[:] = Tnew[:]\n", " \n", "\n", "# Plot the final temperature field\n", "plt.plot(Tnew, -x, '.-', label='T final')\n", "plt.legend()\n", "plt.show()\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.8.2" } }, "nbformat": 4, "nbformat_minor": 2 }