2020-12-01-advection-diffusion.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 5, Evolution in time and space\n",
    "\n",
    "> Advection and diffusion in 1D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
    "import ipywidgets as widgets\n",
    "\n",
    "from IPython.display import HTML\n",
    "from IPython.display import display"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Evolution in time and space: Advection and diffusion in 1D\n",
    "\n",
    "So far we have been looking at dynamics in time, for example how does the temperature of the Earth change over time. But the Earth does not have a single, uniform temperature; rather, at a particular moment in time, different places on Earth are at different temperatures, and those different temperatures change over time due to several mechanisms. \n",
    "\n",
    "\n",
    "In this notebook we will look at two fundamental mechanisms: **advection** and **diffusion**. Let's think about the temperature in the ocean. Since the ocean is a fluid that is in motion, a warm \"parcel\" of water can flow (move) to a new location due to the physical motion of the water itself; this is **advection**.  And even if the water doesn't move, temperature or a higher concentration of some substance dissolved in the fluid can spread out due to molecular mechanisms; this is **diffusion**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook we will restrict ourselves to 1 space dimension.\n",
    "So we will think about the temperature $T$, for example, being a function \n",
    "$$T(t, x)$$\n",
    "of\n",
    "\n",
    "- time: $t$\n",
    "- space: $x$\n",
    "\n",
    "We want to calculate a value of the temperature $T$ at each point $(t, x)$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Temperature profiles and discretization\n",
    "\n",
    "An ordinary differential equation needs an initial value for each variable. Similarly, we will need an initial *function* $T_0(x)$ that gives us the temperature at each position $x$. Let's suppose that the position is restricted to the interval $[0, L_x]$.\n",
    "\n",
    "As usual, to represent the continuous function $T_0(x)$ on a computer, we will need to **discretise** it in some way, i.e. *approximate* the continuous function by a *finite* set of numbers in the computer.\n",
    "\n",
    "The simplest (but far from the only!) discretisation method is to **sample** the function at discrete **grid points** (or **nodes**) $x_i$, for $i = 1, \\dots, N_x$. For simplicity we will take these equally spaced, with spacing $x_{i+1} - x_i =: \\delta x := L_x / N_x$.\n",
    "\n",
    "For example, let's consider an initial temperature profile given by a sine curve:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def t0(x):\n",
    "    return np.sin(2*np.pi * x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and define the grid points as follows:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nx = 20\n",
    "Lx = 1.0\n",
    "\n",
    "deltax = Lx / Nx\n",
    "\n",
    "xs = np.arange(0.025, Lx, deltax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.025, 0.075, 0.125, 0.175, 0.225, 0.275, 0.325, 0.375, 0.425,\n",
       "       0.475, 0.525, 0.575, 0.625, 0.675, 0.725, 0.775, 0.825, 0.875,\n",
       "       0.925, 0.975])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It turns out to be a good idea to take the grid points at the *centre* of each interval, so that we have $N_x$ intervals and $N_x$ grid points, starting at $x_1 = \\delta x/2$ and finishing at $x_N = L_x - \\delta x / 2$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We call such a function of $x$ at a given time a **temperature profile**. Let's draw it as a function and as a heatmap:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fc1addd36a0>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEGCAYAAABCa2PoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOy9d3xU15n//z4zaox6QUJCgCRAVFGFRHEwGANuGMd2NrGJy+7G+nmT7HeTLYmz5LVOsiHJ7n7j+Jddbxy8ibFjbCcbO44hjo0NyIBpoqjQEaigAqj3kaac7x8zjApqc0eaOyOd9+s1r5l77/mc5zya0Txz7imPkFKiUCgUCoWnGPRugEKhUCjGBiqgKBQKhWJEUAFFoVAoFCOCCigKhUKhGBFUQFEoFArFiBCgdwP0JC4uTqakpGjStrW1ERoaOrIN8nGUz+OD8ebzePMXPPf55MmTtVLKiX3Pj+uAkpKSwokTJzRpc3NzWbNmzcg2yMdRPo8PxpvP481f8NxnIURZf+fVLS+FQqFQjAgqoCgUCoViRFABRaFQKBQjwrgeQ1EoFOMDi8VCRUUFZrP5tmuRkZGcP39eh1bpx3B9DgkJITk5mcDAwGHVqwKKQqEY81RUVBAeHk5KSgpCiF7XWlpaCA8P16ll+jAcn6WU1NXVUVFRQWpq6rDq9ZlbXkKIXwshbgohzgxwXQghfi6EKBZCFAohlvS4do8Q4qLz2nPea7VCofAHzGYzsbGxtwUTxcAIIYiNje23VzcQPhNQgB3APYNcvxeY6XzkAL8AEEIYgZec1+cCjwkh5o5qSxXdlOyE91LgTYPjuWQnADa7pMVs4WaLmYqGdqoaO7jRbKapw4LdLofUKxQjjQom7uPu38xnbnlJKQ8IIVIGKbIZeF069ts/KoSIEkIkAilAsZTyKoAQ4m1n2XOj22IFJTvheA7Y2mmxTaCszkjFnleomRxEw4RF2Oz9p0YwCEHEhAAmmfNIKvs3pgbUEx0gob3MUR9A6hYvOqJQKEYCnwkow2AycK3HcYXzXH/nsweqRAiRg6OHQ0JCArm5uZoa09raqlnrr/T1uau+gYrOf6e8PZiGrh6DdiWdEFIydIXmLpD/AEBkoJXJEzqZGtrJhIJ6KMsdXOsl1Ps8NoiMjKSlpaXfazabbcBrYxV3fDabzcP+PPhTQOmv7yUHOd8vUsrtwHaAzMxMqXW16HheXdvUbuFEWT0XT7+MVRqJBCL7Fp79PKYgI6YgI8EBRuxSYrVL2rustHXaHGWKXusub4X2DrjUYCc95BpL73qG+IgQL3k2MOP5fR5LnD9/fsBBaDUoPzghISEsXrx4WGX9KaBUAFN6HCcDVUDQAOcVI0yHVbLn7HXOV7dglxICYsDSBMCJvGIeviOMaUHVJEYEE7tmOiGBxn7rsdjs1LV2UXW9kmvNkvKuRA4fKyMrKw27NHDBuoALx8qZNSmcFWmxRIcGedNNhWJUeOWVV3jppZcAKCwsZMGCBQDcddddvPDCC3o2bcTwp4DyPvB15xhJNtAkpawWQtQAM4UQqUAl8CXgcR3bOeaw2SWnyxv4pMxCsmzuvjBpHQk3X2Ne8AXqOq/wSLQdjCbI2g4DBBOAQKOBSZEhTLrjGZYcz6HT0sWF9lAmB4VTaU2CSesAuHi9hcs3Wlk6LZrstBgCjb40h0ShcI9nnnmGZ555hsrKSlauXEl+fr7eTRpxfCagCCHeAtYAcUKICuB5IBBASvky8AFwH1AMtAN/6bxmFUJ8HfgIMAK/llKe9boDY5Ta1k4+Onudm82dWO3d56fEmMheei/JjWZE4VYmrrODaRos3Db8AXVnueCCrWzdWEbS5EvcmPE0x20rKL7ZCoBdSvJK67l8s4W75yQwJcY00i4qFF7lzJkzZGRk6N2MUcFnAoqU8rEhrkvgawNc+wBHwFGMEFJKTpU3cLi4DmuP2VqxYUHcmT6RabHOra9jtkDaFrZ/73t872vfc99Q6hZI7dYnAJuA6qYODl6qpbKxA4DGdgvvnKogKyWG5WmxGAxqCqhCGz/7+FKv466uToKCgkek7m+uTx+yTFFREfPnzx8Re76GzwQUhe9gttjYc+4GV5y9BIAAg2B+nJEt2dMw9vNlrjWvzED6xMgJfCEzmTOVzRwsrqHTYkdKOFZST0VDB/dmTCI8ZHjbQSgUvsSZM2dYv36967itrY2vfvWrBAUFsWbNGrZs8d8p8+qmtKIXda2dvH28vFcwiY8I5rHsqcyMNvYbTADS04f+ZTYY/emFEGQkR/LkihSm9rjVVdnYwVvHy7neNPwVvAqFr9C3h/Luu+/y6KOP8sorr/D+++/r2DLPUT0UhWOBYsFWrjWaeb/tQbriN0KU4x7vkmnR3DEjbsBAcos9e/awcuVKzU0YTB8WHMDnF0/mRFkDR67UYZeStk4b//vJbjZ2/pR0joFpqnvjN4pxS9/bUt6cNmy327l8+TKzZ892nauoqHCNqRiNA09m8QdUD2W841ztfqle8oeGtXR1tkHlLgKbC7k3YxJ3pk8cMpgAbNq0yaNmDKU3GARZqTE8vGSyYzpyYxHWa7v50/U0TrWld6+yV1u3KHyY4uJikpOTCQ7uHrNJTk6moqICcAQcf0YFlPFOwVYKW5L4oGkVNun4OISJZv7C9gNmT4oYdjVaUym7q58SY+JLy6YQXfc+2C0AfNqylKOt85HWdijY6lE7FIrRJD09nXPneu8K9fDDD/POO+/wN3/zNx7/MNMbdctrnFNQG8y+5kwAiotvkDV7Ag9F7yfS2u5WPdXV1R61wx19dGgQXwr7Le9bV1PZNZHi4hswYwFdMpDPkd/v1gkKha8SGhrKq6++qnczRgTVQxnHFFxrZF/HOtfxvCQrfxHzMZHGNseYhBvk5OR41BZ39SFhk3goaj/TgquZONFx//tk2xwOW9cPoVQoFKOFCijjlKKKJvZduOlYlW4IJDGwlqpTB5lg6HSsdl+4za36tm/f7lF73NYv3EZQYBAPRn3KxROnHOcMgRw35XC8pN6jtigUCm2oW17jkOKbLey9cMNxEJXBpFA7D5m/gz1Rw2p3J6MxbXhQnO0LKNjK/SnlhEeYuRr5MERl8FlxLQFGwZKp0R61SaFQuIcKKOOMioZ2/lx0Helc/J4QEcLnlzxMSOAXSIrJBY27zCYlJXnULk165yr75Khc7vjcav6YX0V5vWPs59OLNYQGBTBr0vjaRVah0BN1y2scUdvayfsFVa6tVKJNgTy0OMm1K7AnOTA8zZ/hqe0Ao4FNC5OYHDXBdX7P2euurVsUCsXoowLKOKG9y8p7pyvptDjmuYcGG/n84mRMQd2d1EceeURz/Z5oR8p2UICBBxclEePc7t5ql7yfX0VDW5dHbVMoFMNDBZRxgM0u2V1QTYvZCji+eB9aNJlIU++9sPy5h3KLkEAjDy2ajCnI0esyW2y8l1+J2WLzqH0KhWJoVEAZ40gp2X/hpuvWjxBw7/xJ/WZDrKur02zHE+1I2440BbJ50WQCjY4VKY3tFv58phr7ADnuFYqxxpo1a9xabJybm8sDDzzgsV0VUMY4hRVNFFU2uY5XzYgjbWJYv2U9WUvi7XUoQ2knRYawcd4k13FpbTtHr3oW9BQKxeCogDKGqWrsIPdijet49qRwMqcNPJXWk7UkXl+HMgztzIRwslJjXMfHSupdibsUikEp2QnvpcCbBsfzCOwR19bWxv3338/ChQuZP38+v/3tb/nBD37AsmXLmD9/Pjk5OUjn9Ms1a9bwzW9+k9WrVzNnzhzy8vJ4+OGHmTlzJt/97ncBKC0tZfbs2Tz11FMsWLCARx99lPb223e42LNnDytWrGDJkiV84QtfoLXV8T/w4YcfMnv2bO644w7effddj/0DHwsoQoh7hBAXhRDFQojn+rn+T0KIfOfjjBDCJoSIcV4rFUIUOa95trHUGKCjy8YHRdWO3O84pgffPTcBIQbemMSTLHKeZqAbLdsr0mKZFtu99f1HZ69TrwbpFYPh3DCV9jJAjtjGox9++CFJSUkUFBRw5swZ7rnnHr7+9a+Tl5fHmTNn6OjoYPfu3a7yQUFBHDhwgGeffZbNmzfz0ksvcebMGXbs2OG6zXvx4kVycnIoLCwkIiKC//7v/+5ls7a2lh/+8Id88sknnDp1iszMTF544QXMZjPPPPMMu3bt4uDBg1y/ft0j327hMwFFCGEEXgLuBeYCjwkh5vYsI6X8DynlIinlIuA7wKdSyp7Lotc6r2d6reE+iJSSPeeuuwbhQwKNPLAwccic7GFh/d8KGw6eaEfTtsEguHd+IhETHBMQuqx2Piiqxmrz711dFaNIwVaw9fmlb/N849GMjAw++eQTvv3tb3Pw4EEiIyPZv38/2dnZZGRksG/fPs6e7c5e/uCDD7p08+bNIzExkeDgYNLS0rh27RoAU6ZMYdWqVQB8+ctf5tChQ71sHj16lHPnzrFq1SoWLVrEa6+9RllZGZcuXSI1NZWZM2cihODLX/6yR77dwmcCCpAFFEspr0opu4C3gc2DlH8MeMsrLfMzTpY1cLWmzXW8YV4CEcPIbnjkyBHNNj3RjrbtCUFGNi1IJMC5DX9NSycHi2s121OMcdrL3Ts/TNLT0zl58iQZGRl85zvf4Qc/+AFf/epX+f3vf09RURHPPPMMZnN30rhbW9wbDIZe290bDAasVsePxb53HPoeSylZv349+fn55Ofnc+7cOX71q1/1W3Yk8KWV8pOBaz2OK4Ds/goKIUzAPcDXe5yWwB4hhAR+KaXs98a6ECIHyAFISEjQPF21tbXV46myo0Fdh51DlVZuTWiaGW3g2tmqXn/YgUhLSxvUp8F8Hkrrqe2R0Ia22ii46Zg+XFoK1y8HkBg2+G8qX32fR5Ox6HNkZCQtLS39XrPZbL2uhU5IxtBx+3+MfUIybQPUMRyqq6uJjo5m8+bNGI1Gdu7ciZSS4OBgqqur+d3vfsfmzZtpaWnBZrPR1tZGS0sL7e3tWK1WVxtvXQsJCaG8vJxPPvmE7OxsXn/9dZYtW9ZLP3/+fA4dOkR+fj7Tp0+nvb2dyspKpk+fztWrVykoKCAtLY3XX3+9l42emM3mYX8efCmg9BcuB5rnuQn4rM/trlVSyiohRDzwsRDigpTywG0VOgLNdoDMzEy5RuNWI7m5uWjVjgolO+k8/Tw7K+YxVUyGSetInLaML2ROGVaCLICf/vSnPPnkkwNeH8znobSe2h4J7Z1Ssquw2pHeuLGIpvI/80D4W4SHTxxw/zKfe5+9wFj0+fz58wNmZbwtY+OiHzvGTHre9jKaMCz6sUeZHQ8fPsyjjz6KwWAgMDCQX/ziF7z33nusXLmSlJQUsrOzCQ4OJjw8HKPRSGhoKOHh4ZhMJgICAly2b10LCwtjzpw5/P73v+fv//7vmTlzJt/4xjcwmUyuMqmpqbz22ms888wzdHZ2AvDDH/6QmTNn8sorr/DFL36RuLg47rjjDs6cOdOvfyEhISxevHhYPvpSQKkApvQ4TgaqBij7Jfrc7pJSVjmfbwoh/oDjFtptAWVM4hxEPFA/nyZrGNBEcNU73JvRidHw+LCrufWB04InWm/ZFkKwYW4Cb5QepqVyNx12Cx/Zl/OIYR/iuHPqsUohrLj1GSjY6rjNNULppTdu3MjGjRt7ncvMzOSHP/zhbWV79gjWrFnTK8DfulZaWorBYODll18eVH/XXXeRl5fX63pLSwv33HMPFy5ccN+RQfClMZQ8YKYQIlUIEYQjaLzft5AQIhK4E/hjj3OhQojwW6+BDcAZr7TaFyjYypW2aM50zHCduivsMJEX/tmtah5/fPjBZyS13rQdEmjkHvOPENIx0+ta1yTyO9JHZNBVMYZI3QIPlcLjdsez+qExLHwmoEgprTjGRD4CzgO/k1KeFUI8K4R4tkfRzwN7pJRtPc4lAIeEEAXAceBPUsoPvdV2vWlvvcEnzY7hpuLiG6SHlDErpMztQcQdO3ZoboMnWm/bTrbls8zkSMNaXHyDz1oW0WAN93jQVaHwJikpKZw541u/m33plhdSyg+AD/qce7nP8Q5gR59zV4GFo9w8n0RKySfm+2i3O7ZSEZY21kWcQgjczrq4dOlSze3wROt126apZMsirnYlcbbTgkUGsKd5OV+YfNl3fmEpFH6I+v/xcy5cb+FK5ONgcEwLXmS6SIihS1PWxXHDwm0EBISwMeIIjkmBUGVN4tSkH+vcMIXCv1EBxY9p77Ly6aUaiMqAyZtYGN1IRUWDI+ti1na37/uePHlSc1s80XrdduoWyNpOfGQk8sZFCIyEyZs40r6UulbPJhcoFOMZn7rlpXCPTy/W0NHlWFcRPmkxd6zYRdrackhJ0VTf008/rbktnmh1se3M9vidBSUcqzFyo9mM1S75+NwNvrhsyqgs+lIoxjqqh+KnXK1p5cL17kVId89JICjAwJtvvqm5Tr20etr+7dtvsWFegmutTnWTmcKKpiFUCoXvUVpayvz583Vtgwoofkin1ca+Czddx3MSw0mJCwXotUWDu+il1dN2cHAwcWHBZKZ078J8qLiWFrNFc50KxXhFBRQ/5LPiWtfGjxOCjNyZHu+65kmSHL20etq+pc1KiSHa1L2BZM9t/xXjj51FO0l5MQXD9w2kvJjCziLPdhrOy8tjwYIFmM1m2tramDdv3m1TfktLS5kzZw7PPPMM8+bNY8OGDXR0OBLj5efns3z5chYsWMDnP/95GhoaAMf44cKFC1mxYgUvvfSSqy6bzcY//dM/sWzZMhYsWMAvf/lLwLH9y+rVq1m1ahXz58/n4MGDHvnVFxVQ/Iyqxg4KrnXfklkzayITnOluAd56S/t+mXpp9bR9SxtgNLBuToLrfPHNVpU7ZZyys2gnObtyKGsqQyIpayojZ1eOR0Fl2bJlPPjgg3z3u9/lW9/6Fl/+8pf7vT11+fJlvva1r3H27FmioqJ45513AHjyySf5t3/7NwoLC8nIyOD73/8+AH/5l3/Jz3/+89s2SP3Vr35FZGQkeXl55OXl8corr1BSUsKbb77Jxo0b+eyzzygoKGDRokWafeoPFVD8CLtd9rrVlTYxlFkJvffeWbFiheb69dLqabundkqMiXlJEa7j3Is3sai0weOOrXu30m7pvX19u6WdrXs920nhX/7lX/j44485ceIE3/rWt/otk5qa6vqSX7p0KaWlpTQ1NdHY2Midd94JwFNPPcWBAwduO//EE0+46tmzZw+vv/46ixYtIjs7m7q6Oi5fvsyyZct49dVX+dGPfkRRUZFHe5P1hwoofkRhZRM1LY5prQEGwZpZ8bfNRrqVjU0Lemn1tN1Xuzp9IiZnj6/FbOV8nU1z3Qr/pLyp/x0TBjo/XOrr62ltbaWlpaXXNvU96TkeaDQaXdvU94eUcsDZiFJK/vM//9O1bX1JSQkbNmxg9erVHDhwgKSkJJ544glef/11j3zqiwoofkJbp5XDV7pzeGSlxhA54fYcJ0VFRZpt6KXV03ZfbUigkTtnTXQdlzTZqVVrU8YVUyP732FioPPDJScnh3/9139ly5YtfPvb3x62LjIykujoaNd4x29+8xvuvPNOoqKiiIyMdCXV2rmz+5bcxo0b+cUvfoHF4phccunSJdra2igrKyM+Pp6nn36av/7rv+bUqVMe+dQXtQ7FTzh4uZZOiyPLYJQpkKUD5IbPycnRbEMvrZ62+9POSginqKKJioYO7BJyL9bwyJLJam3KOGHbum3k7MrpddvLFGhi2zrtO0+8/vrrBAQE8Pjjj2Oz2Vi5ciX79u3jrrvuGpb+tdde49lnn6W9vZ20tDReffVVAF599VX+6q/+CpPJ1Gsn46985SuUlpayZMkSpJRMnDiR9957j9zcXP7jP/4Do9FIRETEiPdQkFKO28fSpUulVvbv369Z6y7X6tvkC3suuh4lNa0Dln3++ec12xlKO5jPntj1VD8a2poWs3zx40vy/2z/SL6w56K8eL1Zsw1/w5ufbW9x7ty5Aa81N9/+3r5R+Iac9rNpUnxPyGk/mybfKHxjNJvndfrzeSD6+9sBJ2Q/36mqh+Lj2O2S/T2msM5MCHOtOemP2NhYzbb00uppeyBtXFgwC6dEcrXEcXzgUg0psaEEBai7xOOBLRlb2JKhtqx3F/Xf4eMUVDRS6xyIDzQKVqdPHLS8J5n29NLqaXsw7fK0WIKdP7lazFbySusHLKtQKFRA8WnMFhtHr3Z/iWWnxRIRcvtAfE9uzVvXgl5aPW0Ppg0JNDIvtnuNz8myBhraujTbUuiL406Nwh3c/ZupgOLDHLlah9nimLYaOSGQxVOihtT44i99X7Y9lHZquIGkKEeuGZtdcuCyWkHvj4SEhFBXV6eCihtIKamrqyMkJGTYGp8aQxFC3AP8/4AR+B8p5U/6XF+DI/Wv884270opfzAcrb9R39ZFYY8V8avT4wgwDh3/q6qqNNvUS6un7aG0QgjWzornzePlSAlXa9oor2tnaqxJs02F90lOTqaiooKamtt/EJjNZre+NMcCw/U5JCSE5OTkYdfrMwFFCGEEXgLWAxVAnhDifSnluT5FD0opH9Co9RsOXq7B7vw1lRw9gekTw4alu3Tpkmabemn1tD0cbXxECHMTIzhb1QzAp5dr2BI9FYNBTSP2FwIDA0lNTe33Wm5uLosXL/Zyi/RltHz2pVteWUCxlPKqlLILeBvY7AWtz1FW18bVmjYAhIA70ycOew2EL63n8Afbw9WunBHnmuFV29LJuepmzTYVirGKz/RQgMnAtR7HFUB2P+VWCCEKgCrgH6WUZ93QIoTIAXIAEhISyM3N1dTY1tZWzdrBsEvJ/mtWmjsdvZOUSAPnTlUx3K7Wjh07NCerGko7mM+e2PVU7y2fjQ02Sp1bsfymooz10wIIGIO9lNH6bPsq481fGD2ffSmg9Pef2XcE7RQwTUrZKoS4D3gPmDlMreOklNuB7QCZmZlS64Bubm6uxwPRt1Gyk4LPfklM7RRiAiMJSrqLpx94iNDg4b9NFy9e1NyuobSD+eyJXU/13vJ5pdXO60dKabl+Gq7vxVT8GSvjm2DhNrfTLfsyo/LZ9mHGm78wej770i2vCmBKj+NkHL0QF1LKZillq/P1B0CgECJuOFqfp2QnnUe/ypFaZ24TSxPLGn9MaNVv3aomMzNTcxP00upp2x1tUICBFSHHoXIXWJo41TablpabcDwHSjzLl6FQjAV8KaDkATOFEKlCiCDgS8D7PQsIISYJ52CCECILR/vrhqP1eQq2crIlhQ67Y7fRCGMrS0LyocC9LbN37dqluQl6afW07a52bsVzxBtvAGCRAXzWughs7W6/TwrFWMRnAoqU0gp8HfgIOA/8Tkp5VgjxrBDiWWexR4EzzjGUnwNfcm4t06/W+15op7WlhlNts13Hq8IKCBB2aHdvy+wNGzZoboNeWj1tu6sVHeWsDu/eofV8Ryo3LDFuv08KxVjEZwIKOG5jSSnTpZTTpZTbnOdellK+7Hz9X1LKeVLKhVLK5VLKw4Np/Ylj1rVYpGOs5PypImaFlDkumNzbMtuXp9/6om23taapTAm6yfSQCo4fvwrAZ60L3X6fFIqxiE8FlPFKfVsXZ8KfBoNjWxVTy1WEAIwmx4CvG5SWlmpuh15aPW27rV24DYwm7gjLp6HBkZyrzDKV8tQfa26DQjFWUAHFBzh8pRZ75AKYvImp4Rb+aX0bmKZB1na3Zw/5+noOX7PttjZ1C2RtJyYimqfXBEJgJEzexCHzKrWth2LcowKKzlQ3dXD5hjMNbVQGd3zxdV6pfR4eKtU0FXX79u2a26KXVk/bmrSpW+ChUgra7idg7jchKoMbzWaKb3qWBlmh8HdUQNERKSUHL3en9Z01KZyEiBBSUlI01+mPWj1te6KdNWM6i6Z2b9j5WXEtdrvqpSjGLyqg6EhJbRuVDR0AGIRg5XRHsqf09HTNdfqjVk/bnmqXpcQQHOj4N2pot7j2+1IoxiMqoOiE3S75rLi7d7IgOZIoUxAAe/bs0VyvP2r1tO2pNiTQSOa0GNe5o1frsNjsmutUKPwZFVB04uKNFmpbHcmaggIMZKV2fylt2rRJc73+qNXT9khoF0+NIsy5PU5rp5WCa42a61Qo/BkVUHTAZpccvVrnOl48NarXfl0nTpzQXLc/avW0PRLaQKOB7LTuHwTHS+tdidEUivGECig6cL66mcZ2C+BIM7tkanSv69XV1Zrr9ketnrZHSjsvKZIok2MdUafFzonSBs31KhT+igooXsZqs/fqnSydFk1IoLFXGb9ZkzFCWj1tj5TWaBCsnB7nOs6/1kB7l1Vz3QqFP6ICipc5U9VMi9nxRWMKMrKonzzxfrUmYwS0etoeSW16Qhhx4Y7NPS02qXopinGHCihexGKzc7yku3eSmRLjygLYE3+dQusJ/tjuvlohBCvSYl3HBdcaae1UvRTF+EEFFC9SWNFIW6djsDYsOIAFyZH9lktKStJswx+1etoeae30iaEkRIQAYLVL8krrNdevUPgbKqB4iU6rjbwet0CyUmMINPb/5/ckNac/avW0PdJaIQQrpnf3Uooqmmg2WzTbUCj8CRVQvER+eSMdXY7eScSEQOZP7r93AvDII49otuOPWj1tj4Y2JdZEYqSjl2KzS/JKVC9FMT7wqYAihLhHCHFRCFEshHiun+tbhBCFzsdhIcTCHtdKhRBFQoh8IYRniyJGGLPFxsny7t5JdmoMRoMYsLwv/eL2hlZP26OhFaL3jK8zlc00tateimLs4zMBRQhhBF4C7gXmAo8JIeb2KVYC3CmlXAD8K9B3is5aKeUiKaVnCc5HmFNlDXRaHNtxRJsCmZsYMWj5urq6Qa+PNa2etkdLOyVmApOjJwBgl5JjJZ79fRQKf8BnAgqQBRRLKa9KKbuAt4HNPQtIKQ9LKW/91D8KJHu5jW7T0WXjdI+tOJZPj8UwSO8EfGNdhTe1etoeLa3osdknwPnqFhraujTbUij8AeErSYGEEI8C90gpv+I8fgLIllJ+fYDy/wjM7lG+BGgAJPBLKWW/CwyEEDlADkBCQsLSt99+W1N7W1tbCdXBgJcAACAASURBVAsLG7LcmVoblxucYyfBgrumBCDE4AFlx44dPP3005raNZrawXz2xK6nel/2+bNKCzfbHf9jU8INZE4KGLS8LzDcz/ZYYbz5C577vHbt2pP93QnypU93f9+y/UY7IcRa4K+BO3qcXiWlrBJCxAMfCyEuSCkP3FahI9BsB8jMzJRr1qzR1Njc3FyG0rZ3WTlzqISUSIcbmxYmMiM+fMi66+rqhqxbD+1gPnti11O9L/uc3tjBb/OuASAEZGROIzYsWFNbvcVwPttjifHmL4yez750y6sCmNLjOBmo6ltICLEA+B9gs5TSdWNaSlnlfL4J/AHHLTRdOVXWiMXmCCYTw4OZPnF4vwg8+eXgj1o9bY+2NilqAqlxoQBICUevqhlfirGLLwWUPGCmECJVCBEEfAl4v2cBIcRU4F3gCSnlpR7nQ4UQ4bdeAxuAM15reT90dNkoqOgeO8lOjRnyVtctjhw5otmuP2r1tO0Nbc91KZdvtlDX2qnZpkLhy/hMQJFSWoGvAx8B54HfSSnPCiGeFUI86yz2L0As8N99pgcnAIeEEAXAceBPUsoPvexCL06VN9BldczsigsLYkb88H8JP/bYY5rt+qNWT9ve0CZEhJA2sbuXckytS1GMUXwmoABIKT+QUqZLKadLKbc5z70spXzZ+forUspo59Rg1/Rg58ywhc7HvFtavTBbbOT3mNmVnRY77N4JwO7duzXb9ketnra9pc1O7e6lXLqheimKsYlPBZSxwqmy7t5JbFgQM93onQB0dmr/svFHrZ62vaWdFBnSayxF7fGlGIv40iyvMYH58k5O7/kUOtsgMJLslZ9DiBS36nj88cc12/dHrZ62vanNTouhpLYNGou4cHEv2ad2Eh0RDQu3QeoWze1QKHwF1UMZSUp2cvrTF+nqbAMgRl5j5uVnoGSnW9Xs2LFDcxP8UaunbW9qEyMnkGLPh8pdyK4mjrXNg/YyOJ7j9mdEofBFVEAZQcynn+d0a4rrODu0CIO9DQq2ulXP0qVLNbfBH7V62va2NrvxR2B37Ot1wZxCgzUcbO1uf0YUCl9EBZQRJL/WRKc9CIAbJVdJDyl3XGgv17FVCl8iyVbItGBHLvrLl29yvG2e44L6jCjGACqgjBBmi41TlmzXcUvZBQzCudDfNNWtuk6ePKm5Hf6o1dO217WmqWSHOpZIFV+5wQVzCo3WMLc/IwqFL6ICyghRcK2Rzvh7wBBIVEAL/3CXc9qw0eQYdHUDT/bE8ketnra9rl24jckT2pgSdJ2sZWnYpYHjHUvc/owoFL6ICigjQKfVxqnyRojKgMmbyIqt5+0jgGkaZG13ewbPm2++qbkt/qjV07bXtalbIGs7yyc2cOpUGQRGcj72b2lK+AvN7VAofAUVUEaAgmtNmC2OHYUjE5cw57GDBGf9X3ioVNN00OBg7ZsH+qNWT9u6aFO3kPzFU8Smb4LZ38AeuYDjal2KYgygAoqHdFntnCrvnSveYBA88MADmuscb1o9bevp81891p1C+FxVM00dKqujwr9RAcVDCit654qf48zG+NZbb2muc7xp9bStp8+5f36vV1ZHlXte4e+ogOIBXVY7J8t69E5SunPFr1ixQnO9402rp229fV6R1r3H17nqZprNqpei8F9UQPGAospG2p29k/CQAOYmdeeKb21t1VzveNPqaVtvn5OjJzA5ytFLsdlVL0Xh36iAohGrXXKitPfYibFHrviioiLNdY83rZ629fZZCMHyHr2Us1Wql6LwX1RA0Uhpk7137yQxotf1nJwczXWPN62etn3B5ykxE0iKCgEcvZSTPX6oKBT+hAooGrDY7FxutLuOl6XEEGDs/afcvn275vrHm1ZP277gsxCiV76UosomWlQvReGHuB1QhBCPCyHeFkLsFEK8KYTwLFVf77rvEUJcFEIUCyGe6+e6EEL83Hm9UAixZLjakaSosgmz1bGtSnhIAPN6jJ3sLNpJyospfP/490l5MYWdRcPfRXa8af213aPh87RYE4mR3b2UE2X991JuaQ3fN2hutz9p/bXd/uqzp2jJh3KnlPJLtw6EEC8Bns2fdNRjBF4C1gMVQJ4Q4n0p5bkexe4FZjof2cAvgOxhakcEq83e65bE0mnRrt7JzqKd5OzKod3SDilQ1lRGzi7HbY0tGYMvcBxvWn9t92j6vDLtId47XQnAmYomlqXEEBYc0L8WD9rtJ1p/bbev+zya65203PIKFkLcL4RYIIS4D5gwQm3JAoqd6Xy7gLeBzX3KbAZelw6OAlFCiMRhakeEM1XNtHZaAQgLDiBjcqTr2ta9W11vJM5Q1m5pZ+veobcmH29af233aPqcEmsiIcLRS3FM+qgfWIsH7fYTrZ62x6rPbZ1WXj9cyqFKC9fq2/urwiOElHLwAkL8XynlP/Y4NgEPA8nANeAPUkqPWyaEeBS4R0r5FefxE0C2lPLrPcrsBn4ipTzkPN4LfBtIGUrbo44cIAcgISFh6dtvvz3sNkop2VtupaVLUlFRgb3qDJtWzCc/P5/m5mamZk3l6N6jJE5NpO5GHV2dXSxeuZhzp8+REJLA3XffzYcffkh6ejoAly5d4p577uGTTz7hhvkGcxfP5fTh0wQFBxGbEEt1eTXL1y2n/Hg5ERERLFq0iAMHDrBo0SJqa2upqKjggQceYPtb24mIjmBK6hQOfnSQFetWUFNdQ31NPTmP5bB7924mTZpEcnIyJ06cYMWKFVy8eJH6+npXmycmTqSlqQVzu5kFWQu4cv4KkUTywAMPsHv3blJTUxFCcPXqVe666y6OHj1KaV0pi1ctJu/TPIwBRhKnJFJRUsGyO5dxo/AGQUFBLF++nH379jF37lw6OjooKSlx1dlEE9PnTKfweCEhphDCI8Opqa5x+RwTE8OsWbM4cuQImZmZVFRUcP36dZfPMRNjmJg4kSN7j/C5jZ/jWsk1mhuaXT4nJycTFxdHfn4+q1ev1vQ+mc1mysvLR+R92r17N60BrUxJncLZU2cJDQ8lOCSY+pp6l940KZWG8OlUVVcxLTmZqJunaa6vcfk8MXEiEdER5H2ax5r713Dl/BXaWtpcPqempjJhwgTOnTvnep+6urpIWJBA3qd5JKcmU32tGpvVxrI7l3H6s9OkxKb0ep8aGhqorq7W/X26pU+cmkhwcDCnDp9i3eZ1nDt9DpvFxpOff3LA/6egoCDCp4dz+vBpUmamUFZchpSS5euWc3TvUWZPnu16n2bNmkVbW5tb79NA/0+evk933303v/7fX5OcmgxA4fFCNj66kdOfncYYaOQLG74w4P9TREQEQZODKDxeyPS507l85jIGg8Hl8/xp8wmcsoCjV24SPzGeeWmTuTM5ACHEgN95A7F27dqTUsrM275fhxFQTkkplwxaaAQQQnwB2NgnKGRJKf+2R5k/AT/uE1C+BaQNpe2PzMxMeeLECbfa2dFl41R5Ax98dprvfXkdgT0G41NeTKGsqew2zbTIaZR+o3TQev1Bm5uby5o1a0bErjfb7YnW2z5LKXnzeDk3mx356pdOi2Z1+sTbtfuBtcO37Y52UJ9H0e5I64er7euvt+yOhn4w7bmvFvPrQyVYbJLS0lL+z0MrmT4xbMj29IcQot+A4kuzvCqAKT2Ok4GqYZYZjnZEmBBkZNWMONZNDegVTAC2rduGKdDU65wp0MS2dUNvTT7etHra9mWf+874KqxopL3Lert2rgft9hOtv7bbV30+VdaIxeboQEQGC9LiQofVHreQUg76AGxACfA+8CPgMSADCBxK684DxwSBq0AqEAQUAPP6lLkf+DMggOXA8eFq+3ssXbpUamX//v39nn+j8A057WfTpPiekNN+Nk2+UfjGsOv0dW1/Pnti11vt9kSrh892u13+5kipfGHPRfnCnovywKWbt2m5F80+D6UdzOfRtDsa+uFoh/pfHgs+t3da5X/tu+z6TL2xa69b7ekLcEL29z3e38leBeA0kARsAP4BeBU4ATQBZ4bSu/MA7gMuAVeArc5zzwLPOl8LHLO5rgBFQOZg2qEeoxFQxjLKZ+9x+Uaz65//v/Zdlu2d1l7Xn3/+ec11D6UdzOfRtDua+sG0Q73HY8HnQ5drXJ+n1w+XyH379nnUroECyrCmDUspq3DcQtpz65xwjOTMGI5+uEgpPwA+6HPu5R6vJfC14WoVCn9l+sQw4sKDqW3pdG1CesfMONd1f9wdQO2moI9ts8VG/rVG1/nstFiqzpd61K6BGM4Yykv9nXQGqssj3B6FQoFjLGV5aozruKBHmgTwz90B1G4K+tg+Vd5Al9Wxs0dsWBAz47UNxA+HIQOKlPJ/Rs26QqEYkBnxYcSFBQGOVAmneyRyy8jI0FyvP2r1tO3PPpstNk6Xd/dOslJjNE0THi6+NMtLoVD0QAhBVo8ZX6evNbpSTYeFaf+V6Y9aPW37s8+nyxtdvZOY0CDS48M9as9QqICiUPgwM+PDiO3RS7mVbvrIkSOa6/RHrZ62/dXng4c+4/S129OTjyYqoCgUPozBIMjqMZaS7+ylPPaY9j1Z/VGrp21/9XnBnffRaXH0TqJNgcxKGN3eCaiAolD4POnx4USbAgHotNg5Xd7I7t27Ndfnj1o9bfujz51WG79994+u42Ve6J2ACigKhc/j6KX0HEtpoK29Q3N9nZ2dfqfV07Y/+lxwrQmz2aGNnBDInEkRQyhGBhVQFAo/YPakcKJMgdBYRGfRC2QYfwPvpUCJ+7kuHn/8cc3t0Eurp21/87mreCen/vwtliRdgQsvkhV8zCu9E1ABRaHwCwwGQVbQMajcBZYm/ueAkc7WSjie43ZQ2bFjh+Z26KXV07Zf+Vyyk8JPf0qHuYPjeVeJsFcy58pXNP3w0IIKKAqFnzCn4jtEGhyzdhInT6SgYxbY2qFgeHk2brF06VLNbdBLq6dtf/K56/S/cLIlFYApyTFkhZ7FaG9z+zOiFRVQFAo/wdBRRlboGQCmT0/gVNtsuuwB0F6uc8sUvkJRXTDtdkeStkXpYcydUOK44KXPiAooCoW/YJrKnJBSIoythIQE0mEPpqAjHUxT3arm5MmTmpugl1ZP2/7is8Vm54Rlpet4dewljMIxbdjdz4hWVEBRKPyFhdswBoSQFXrWdepkxwK65g8vz8Ytnn76ac1N0Eurp21/8bmoson2ifeCIZBwYztzJ1x1XDCaYKF7nxGtqICiUPgLqVsgaztzY22EG9shMJKOxEcpCrzXrWrefPNNzU3QS6unbX/w2WKzc6K0HqIyYPImMmPrCRASTNMga7vjs+MFhrV9vUKh8BFSt2BM3UJWRSN7z98E4ERpAwuSo27LIDoQwcHBms3rpdXTtj/4XFTZRFunc5+3hMXMX7Ufhvl5GElUD0Wh8EPmJkYQHuL4PdjeZaOwomnY2gceeECzXb20etr2dZ+7rHbySupdx5kp0QToEEzARwKKECJGCPGxEOKy8zm6nzJThBD7hRDnhRBnhRB/1+Pa94QQlUKIfOfjPu96oFB4lwCjgWUp3Xt8nSyrx2KzD0v71ltvabarl1ZP277uc2FFI+3OXDnhIQFkTI70yKYn+ERAAZ4D9kopZwJ7ncd9sQL/IKWcgyOf/NeEEHN7XP+ZlHKR86EyNyrGPPOSIggLdvRS2jptFFUOr5eyYsUKzTb10upp25d97rTaOFHWe0dhvXon4DsBZTPwmvP1a8BDfQtIKaullKecr1uA88Bkr7VQofAxAowGMlO6O/MnSxuG1UtpbW3VbFMvrZ62fdnngmtNrkyeERMCmZekX+8EfGdQPkFKWQ2OwCGEiB+ssBAiBVgMHOtx+utCiCeBEzh6Mg39SBFC5AA5AAkJCeTm5mpqcGtrq2atv6J89j1sdsn1CitmqwTgtboSZkQbB9Xs2rWL2NjYAa8P5vNQWk/sjqZ+MO1Q77Gv+txlk3xcZuFWZujF8UYOHigbVr2j9rmWUnrlAXwCnOnnsRlo7FO2YZB6woCTwMM9ziUARhw9rm3Ar4fTpqVLl0qt7N+/X7PWX1E++yanyxvkC3suyhf2XJQv5xbLTott0PKVlZWDXh/M56G0ntgdTf1g2qHeY1/1+ciVWtf7/utDV6XVZh92vZ5+roETsp/vVK/d8pJS3i2lnN/P44/ADSFEIoDz+WZ/dQghAoF3gJ1Synd71H1DSmmTUtqBV4Cs0fdIofAN5if1nvGVf61x0PLbt2/XbEsvrZ62fdFns8Xmyt4JkJ0ai9FLOwoPhq+MobwPPOV8/RTwx74FhBAC+BVwXkr5Qp9riT0OP4+j56NQjAsCjAaWp3XfFjlRVu/KPd8fntyC0Uurp21f9PlUWUOvbIyzJ41+Nsbh4CsB5SfAeiHEZWC98xghRJIQ4taMrVXAE8Bd/UwP/nchRJEQohBYC3zTy+1XKHRlTmKEI18KjqyOp8r6HUIEYM2aNZrt6KXV07av+dzRZeN0j17o8umxXst3MhQ+EVCklHVSynVSypnO53rn+Sop5X3O14eklEJKuUD2mR4spXxCSpnhvPagdA7wKxTjBaNBsGJ6z6yOjbR3Wfst+84772i2o5dWT9u+5vPJsga6rI7eSWxYEOnxvtE7AR8JKAqFwnPS48OJCwsCHKunT5T230sZS7/Wx7K2P317l5X8a93v64o03+mdgAooCsWYwdCnl1JwrZHWztt7KVVVVZpt6KXV07Yv+ZxX2oDF5pgiHhcezIz4MI/qH2lUQFEoxhDTJ4YRH+HYUNBqlxwvqbutzKVLlzTXr5dWT9u+4nOz2UJhj7GTFWmxOOYq+Q4qoCgUYwghBCunx7mOz1Q209Rh6VUmJydHc/16afW07Ss+H71Sh9Xu6J0kRoYwfWKoR3WPBiqgKBRjjJRYE5OjJgCOlfTHrvbupYy1NRljVdtTX9faybnqZtf5VTPifK53AiqgKBRjDiF6j6Wcr26hoa3LdZyYmNifbFjopdXTti/4fPhKHY5NQWBarIkpMSaP6h0tVEBRKMYgU2JMTHV+6dil5PCV7l5KZmam5nr10uppW2+frzeZKb7ZvUnkqhlxgyj0RQUUhWKM0vOL59KNFq43mQHHhoNa0Uurp209fX7//fc5VFzrOk5PCCchIsSjOkcTFVAUijHKpMgQZiZ0Tys9eLkGKSUbNmzQXKdeWj1t6+nzwuWruVbfDoBBCFZO92wrl9FGBRSFYgyzanocBufgbUVDByW1bWNiCu140Eop+ehwvut4XlIE0aFBmuvzBiqgKBRjmOjQIDKSI1zHnxXXUlJSorm+0tJSXbR62tZLe/lmq0sfYBBkp8UMLvABVEBRKMY42amxBAUYoLGI2hP/zp0Rr8F7KVCy0+26fGVNxljX2q7s5LN3/5YVaRVw4UUWBZ0iPCRQczu8hQooCsUYJzQ4gCVBJ6FyF1iaePHjECxtFXA8x+2g4gtrMsa8tmQnhZ/+O40dFo4cLSbYVsOyiv9P0w8Ab6MCikIxDlh6458JFY6FcaFR0eS3zwJbOxRsdauelJQUzW3wRKunbW9rzae/x9HmmQDERIeRHXqGENnk9nulByqgKBTjgCBzCdmhjrxzWVlp5LXNxWwPgvZyt+pJT0/X3AZPtHra9rb2eE0EZrtjP7b1K+NZaHIO7Lv5XumBTwQUIUSMEOJjIcRl53P0AOVKnYm08oUQJ9zVKxTjFtNU5k+4QnSAo5fSaQ/iWNt8ME11q5o9e/ZoboInWj1te1Pb1G4h37LMdXxHWD4BwpH7xN33Sg98IqAAzwF7pZQzgb3O44FY60yu1XMJqjt6hWL8sXAbxoAQVoUVuE7ld8yjYdaP3Kpm06ZNmpvgiVZP297UfnalFlvCejAEkhhYy8xgZ6/EaIKF2zS3w1v4SkDZDLzmfP0a8JCX9QrF2CZ1C2RtZ0a0gclBNRAYiX3yZg503ulWNSdOnBi60Cho9bTtLW11UwcXr7dAVAZM3sTq+BuODSBN0yBru+M99HEC9G6Ak4RbaXullNVCiPgByklgjxBCAr+UUm53U48QIgfIAUhISCA3N1dTg1tbWzVr/RXls78zGaJ3EDjBTlmFFdkIpXln6ay6SLyp+7flYD4fOXKEWbNmabLuiXY0bQ/1HnvDZyklByqt1Hc4doCcHB7JpUkv4loWWQaUDdxGdxmtz7XXAooQ4hNgUj+X3Jm6sEpKWeUMGB8LIS5IKQ+40w5nENoOkJmZKbWm6MzNzfU4vae/oXweO4Sevc7ZKud4SlgQq7OnuVLJDuZzeno6SUlJmmx6oh1N20O9x97w+dKNFiKs1UQARoPgyRXTiDKN3qr40fpce+2Wl5Tybinl/H4efwRuCCESAZzPNweoo8r5fBP4A5DlvDQsvUKhcLByRpxjsSNQ29pFUWXTsHRqHcrIa7usdg5cqnEdL5wSNarBZDTxlTGU94GnnK+fAv7Yt4AQIlQIEX7rNbABODNcvUKh6CYsOIBlKd1beRy5WofZYhtSp6YNj7z2RGk9LWYrAKYgI9mpvr/FykD4SkD5CbBeCHEZWO88RgiRJIT4wFkmATgkhCgAjgN/klJ+OJheoVAMzJKpUURMcGzn0dFl4+jV2/PP98WTW1aeaPW0PZrapnYLJ8saXMerZsQREmjUbE9vfCKgSCnrpJTrpJQznc/1zvNVUsr7nK+vSikXOh/zpJTbhtIrFIqBCTAaWD2zO2dKwbUm6lo7B9V4MpDr6SCwXrZHU/vp5RpXnvhJkSHMS4oYtLyv4xMBRaFQ6MOM+DCSox355+1Ssu/CTeStXLP98Mgjj2i25YlWT9ujpS2tbeNKj0yMa2ZN9Mk88e6gAopCMY4RQrBmVnyvnCkVrQMHFNVDGRmtzS75tMdA/NykCBIjJ2i24yuogKJQjHMmhgezeGqU6/hMrW3AAfq6uqHHWQbCE62etkdDe7q8gfq2LgCCAgzc4cN54t1BBRSFQkF2WgzhIY5laWar5MgAA/QqH4rn2qYOS68JEMvTYgkN9pU15p6hAopCoSA4wMjq9Imu44JrjdxsNt9WTq1D8UwrpWT/hZtYbI7binHhwSyaEtWf1C9RAUWhUAAwMz6MabEmAKSEvf0M0GdkZGiu3xOtnrZHUlt8s5WS2jYAhIB1s+MxGvx7IL4nKqAoFArAMUC/dlY8t77frjeZKazovYI+LCxMc/2eaPW0PVJas8VG7sXugfgFyZEkRfn/QHxPVEBRKBQuokODSI/uXlh3qLiWZrPFdXzkyBHNdXui1dP2SGmPXKmjtdOxIj402MjK6WNjIL4nKqAoFIpepEcbiAl17CXVZbWz73z3ra/HHntMc72eaPW0PRLa6qYOCioaXefXzIr36xXxA6ECikKh6IXRIFg/N4Fba+xKats4X90CwO7duzXX64lWT9ueaq02Ox+fu8Gt4ajUuFBmxnt2+89XUQFFoVDcRlLUBBb2mH306aUa2jqtdHYOvjXLYHii9VSvp/bo1XrqWrvXnKydHe/3K+IHQgUUhULRL6umx7k2j7w1oPz4449rrs8Trad6vbR3P/AwJ8q6txa8Y0Yckc6/6VhEBRSFQtEvQQEG7p7jTH7aWMSlT/+Vn/7dQngvBUp2ul3fjh07PGqPJ3qva0t2Yv1DGj/69oPI8y9CYxFTYkwsSI7U3A5/QAUUhUIxINNiQ5kXcBoqd4GlCevEebS03ITjOW4HlaVLl3rUFk/0XtWW7ITjORy5GUlcUiJYmgiq+gPrQw+M2Vtdt1ABRaFQDMrqhq2EC8cMJYsM5OPm5UhrOxS4k717HFGwlYqOUE62z3Gd+lzoMSIv/LOOjfIOKqAoFIpBCTFfZWPkEYSQBAcHUtaZSH5HOrSXu1XPyZMnPWqHJ3pvas2t1/mwaSVSCoKDA5kafJ2MCcVu/738EZ8IKEKIGCHEx0KIy87n6H7KzBJC5Pd4NAshvuG89j0hRGWPa/d53wuFYoximsqUoJssNZ1nxowEAA61LKY2cK5b1Tz99NMeNcMTvbe0Ukr2dt5Liy0UgPnpUWyIOOKYgm2aqrkN/oJPBBTgOWCvlHImsNd53Asp5UUp5SIp5SJgKdAO/KFHkZ/dui6l/KCvXqFQaGThNjCaWBFWyMRAR7paqwjhzxOex2qzD7uaN99806NmeKL3lvZsVTOXIr4MBsdMrvURxwg3doDR5Pg7jnF8JaBsBl5zvn4NeGiI8uuAK1LKslFtlUKhgNQtkLWdgNAp3Bt5mICgcJi8idrghb32phqK4OBgj5rhid4b2oa2LkfSrKgMmLyJBdFNzAipBNM0yNru+DuOccRg6T691gghGqWUUT2OG6SUt9326nH918ApKeV/OY+/BzwNNAMngH+QUjYMoM0BcgASEhKWvv3225ra3Nra6vFmd/6G8nl8MJTPV5tsFNzsTsCVOSmAKeFD/zYtLy9n6lTtt3080Q+mHcrf4di12iWfVlhp7nR8n4YHCdZMCSDAR3cS9vRzvXbt2pNSyszbLkgpvfIAPgHO9PPYDDT2KdswSD1BQC2Q0ONcAmDE0ePaBvx6OG1aunSp1Mr+/fs1a/0V5fP4YCif7Xa7/FNhlXxhz0X5wp6L8j/3XpK1LeYh633++ec9apcn+sG0Q/k7lF273S4/6PH3+Pknl+SN5g73G+lFPP1cAydkP9+pXksTJqW8e6BrQogbQohEKWW1ECIRuDlIVffi6J3c6FG367UQ4hXAs02DFArFgAghWDcnnpqWTurburDYJH8qquZLy6YSFDBwT2XFihUe2fVEP5ragoomLlxvcR2vnR1PfHiIZnv+jK+MobwPPOV8/RTwx0HKPga81fOEMwjd4vM4ej4KhWKUCA4wcv+CRAKNjls6da1dfHT2+m0JuXrS2trqkU1P9KOlrWrs4MCl7nGk+ZMjmT95bK+GHwxfCSg/AdYLIS4D653HCCGShBCuGVtCCJPz+rt99P8uduh/dgAADNJJREFUhCgSQhQCa4FveqfZCsX4JS4smLtmJ7iOi2+2cuRK/7noAYqKijyy54l+NLTNZgu7C6uw2R1BNCEihLWzJvZbdrzgtVtegyGlrMMxc6vv+Srgvh7H7UBsP+WeGNUGKhSKfpmbFEFNayenyhxzYI6V1BMTFsTsSRG3lc3JyfHIlif6kdZ2Wm38Mb+Ktk7H5ISQQEePLcDoK7/R9WF8e69QKDzmczPiSIkzuY4/PnuD603m28pt377dIzue6EdSa7dLPjxzndoWx7b2RoPggQWJY3oX4eGiAopCofAIg0Fw7/xEV5ZHq13yx/xKGtq6epWLjb3t5oJbeKIfKa2Ukk8v13C1ps117q7Z8UyJMfUnHXeogKJQKDwmJNDI5kVJrrS27V023j1dSUuPfPRr1qzxyIYn+pHSHiupJ7+8O5XvspSYcT0I3xcVUBQKxYgQZQpi86Ik18yv5g4L752uxGxxjDO88847HtXviX4ktKfLG3pNOkhPCGfVDM96XWMNFVAUCsWIkRQ1gfsXJGFw5v2obe3iD86g4s89lHNVzb22mZkWa2LjvIQxn9/EXVRAUSgUI0pqXCgb5nVPJ75else7v36S0t+t15ztEaCqqkpzmzRpS3bCeykce/OL7Hn3OWh0TB9OigrhgQVJ435GV3+ov4hCoRhx5iRGsG5OvONLuHIXN9rg/ZIUOlqva8r2CHDp0iXN7XFb68y6WFAbxL6yeGRXE1TuIq6zgM2LJg+6I8B4Rv1VFArFqLAgOYq7O19ESMdsrwVZ8/nf+rtp6hKasj16cx2KzN/K8eZU9jUvY8XyGQDEG2/waOc/uyYeKG5HBRSFQjFqZBg+4+6IYwghqalpoc4ayW/rN3Cjyf2tULy1DsVml3xyfRKftSwEoKamhcTAWh6J3suEziua2zAe8ImV8gqFYoximsp8rhIobBhmrsAmoc02gf9tfYSNN1qYmRA+7KoSExOHLuSh1myx8UFRNWXWJUATAGvmwoNR+wg2WB25TRQDonooCoVi9HBme5wVUsbD0XsJNnSBIRBL/AZ2F1bz6aUa115YQ5GZeXv6jeEyHO31JjM7j5VTVtcOk9aBIZA5E0p4OHq/I5iMk6yLnqACikKhGD2c2R4xTSM5qJYvTj5PZNr9jqyGwKmyBn5/8hqN7V1DVAS7du3S3IzBtFJKTpU38LsT12jucC7EjMpgxbJ72DipCqOQ4yrroieoW14KhWJ0Sd3i+iKOBR632Pjo7HXX9iVVjWbeOFrGyhlxLJ4SNeDajg0bNmhuwkDahrYuDlZaCbd2rzEJDjSwYe4kZsSnA1/SbHM8onooCoXCq4QEGnlwYRJ3zIxzLYC02CSfXqzht3nXqGrs6Fc3ktOGO602DhfX8sbRMuo6um+5JUSEsCVrGjPix1fa55FC9VAUCoXXEUKwLCWGqTEm9py9Tm2r45ZXdZOZ3+ZdY2ZCGFmpMb0yH5aWlmq2d0vbZbVztqqJvNJ619bzAAYhWJYSTVZqjFqw6AEqoCgUCt1IiAjhsaypHC+t50Rpg2uA/vKNVi7faGVarIkFHCSlbCs5E0vgvR2OgfHhjmWU7ISCrfxFdB2HXztMQUQO5rAFvYpEhwgey54ybtP2jiQ+EYqFEF8QQpwVQtiFEANOxxBC3COEuCiEKBZCPNfjfIwQ4mMhxGXnc7R3Wq5QKDwlwGhg5fQ4nlqRwqxJvacRl105xq4DH/FK2RK+9cFkztcJWo78HfLq4CvtrTY7lYVvcnTvz3jr2my++UEKx+omYi7b7dpCJTwkgI3zJnFncoAKJiOEr/RQzgAPA78cqIAQwgi8hCMFcAWQJ4R4X0p5DngO2Cul/Ikz0DwHfHv0m61QKEaKSFMg92UksnRaNCdKG7h8swV5fS/YLZgJxhY+hQ+b0qAJQv6cS/SyzxEaHECg0YDRILDY7HRabTS2W2jqsCDPHwLLLABioh1rSrBbiKzdxZLl65iXFEGg0cDNS2qDx5HCJwKKlPI8MNTOnVlAsZTyqrPs28Bm4JzzeY2z3GtALiqgKBR+SUJECPcvSKSp/f+1d78hctx1HMffn/vvNdFUL4mxTZprqVoRg21M2lIlqYptnpRCH5SWVopSRCKCCAUL+kDF+ExEpIRSig/sIdraiNFgKWeFmDYtNEljUGOC8WKwqSYt6XHJ/fn6YKZhPffu5nK/nZ3Jfl6w3M7O7N73wy7zndmdmd8Qh4/t5cjMNbw1vYxNm669uMzE+XFONRkV8n9Mvnnx7i2b17O+f4wbBo5zXf9JutY+2aLqO5siip1UVAZJo8DXI+LlJvPuAe6IiC/m0w8AmyNiu6SzEbGiYdkzEdH0ay9JDwMPA6xevfqmkZGRS6r13LlzLFvWWUeCOHNnqFTms4eI6Qucmezh9EQvp8/3cmayl6nogYFVcz5NgismTzHUN8FQ/ySr+i/Q352v67r6Lp4HAxXLW5KlZt66desrEfF/P0+Utoci6Tng/U1mPRoRzxZ5iSaPLbobRsROYCfAxo0b41LHSRgdHV3y+A5148ydoVKZj5/Mrk48Mw59QB9E1yDnPr6TsytvY2JymvNTM8xE0NfTRW93F+8e6GXFYC+9J56Cl74G0+MwSXbrHsxPUNxy8V9UKm9JWpW5tIYSEZ9Z4kuMAWsbpq8G3hnk4F+S1kTEKUlrgNeX+L/MrAreOZrrwKMwfgIG16EN32X58P0seBWwJs9d1BFitmiV+A2loP3A9ZKGgZNkp7Del8/bBXwe2JH/LbLHY2Z10HCmfanPtUWrymHDd0saA24Bfi1pT/74ByTtBoiIKWA7sAc4AvwsIg7nL7ED+Kykv5IdBbaj7AxmZp2uEnsoEfEM8EyTx/8JbGuY3g3sbrLcv4FPt7JGMzObXyX2UMzMrP7cUMzMLAk3FDMzS8INxczMkqjUmfJlk3Qa+PslPn0IeCNhOXXgzJ2h0zJ3Wl5YeuZrImLl7Ac7uqEshaSXm1164HLmzJ2h0zJ3Wl5oXWZ/5WVmZkm4oZiZWRJuKJduZ7sLaANn7gydlrnT8kKLMvs3FDMzS8J7KGZmloQbipmZJeGGsgBJd0j6s6Sj+Xj1s+dL0g/z+Qcl3diOOlMpkPf+POdBSXslbWhHnSktlLlhuU9Ims5HD621IpklbZH0qqTDkn5fdo2pFfhsv0fSryQdyDM/1I46U5H0hKTXJb02x/z0666I8G2OG9AN/A24lmy8uAPAR2Ytsw34DdmIkjcDL7a77hbnvRW4Mr9/Z53zFs3csNzzZFe7vqfddZfwPq8A/gSsy6dXtbvuEjJ/A/h+fn8l8B+gr921LyHzp4AbgdfmmJ983eU9lPltAo5GxLGIuACMAHfNWuYu4CeR2QesyEeNrKMF80bE3og4k0/uIxs5s86KvMcAXwF+weUxGmiRzPcBT0fECYCIqHvuIpkDWC5JwDKyhjJVbpnpRMQLZBnmknzd5YYyv6uAfzRMj+WPLXaZulhsli+QbeHU2YKZJV0F3A08VmJdrVTkff4gcKWkUUmvSHqwtOpao0jmHwE3kA0tfgj4akTMlFNeWyRfd1VigK0KU5PHZh9nXWSZuiicRdJWsoZyW0srar0imX8APBIR09nGa+0VydwD3EQ2cN27gD9K2hcRf2l1cS1SJPPngFeB24HrgN9J+kNEvNXq4tok+brLDWV+Y8DahumrybZeFrtMXRTKIuljwOPAnZGNlllnRTJvBEbyZjIEbJM0FRG/LKfE5Ip+rt+IiLeBtyW9AGwA6tpQimR+CNgR2Q8MRyUdBz4MvFROiaVLvu7yV17z2w9cL2lYUh9wL7Br1jK7gAfzIyZuBt6MiFNlF5rIgnklrQOeBh6o8dZqowUzR8RwRKyPiPXAz4Ev17iZQLHP9bPAJyX1SBoENgNHSq4zpSKZT5APJS5pNfAh4FipVZYr+brLeyjziIgpSduBPWRHiTwREYclfSmf/xjZUT/bgKPAONlWTi0VzPtN4H3Aj/Mt9qmo8ZVaC2a+rBTJHBFHJP0WOAjMAI9HRNPDT+ug4Pv8beBJSYfIvg56JCJqe1l7SU8BW4AhSWPAt4BeaN26y5deMTOzJPyVl5mZJeGGYmZmSbihmJlZEm4oZmaWhBuKmZkl4YZiZmZJuKGYmVkSbihmFZKPuXJQ0oCkK/JxOT7a7rrMivCJjWYVI+k7wADZRRnHIuJ7bS7JrBA3FLOKya81tR+YAG6NiOk2l2RWiL/yMque95IN8LScbE/FrBa8h2JWMZJ2kY0oOAysiYjtbS7JrBBfbdisQvKREaci4qeSuoG9km6PiOfbXZvZQryHYmZmSfg3FDMzS8INxczMknBDMTOzJNxQzMwsCTcUMzNLwg3FzMyScEMxM7Mk/guaVT0rLapzfwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "t = np.arange(0,Lx, 0.001)\n",
    "\n",
    "plt.plot(t, t0(t), label = \"$T_0$\", lw = 3, alpha = 0.5)\n",
    "plt.scatter(xs, t0(xs), label = \"sampled\", color = \"orange\")\n",
    "plt.scatter(xs, np.zeros((xs.size)), color =\"green\", label = \"x nodes\")\n",
    "plt.hlines(0,xmin=0, xmax = Lx, color = \"black\", ls = \"--\", lw = 0.5)\n",
    "for x in xs:\n",
    "    plt.plot([x, x], [0, t0(x)], color = \"black\", ls = \"--\", lw = 0.5)\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"$T_0$\")\n",
    "plt.grid()\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will denote by $T^0_i$ the initial temperature at grid point number $i$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A useful way to think about $T^n_i$ is as the (spatial) average of $T(t_n, x)$ over the interval of positions between neighbouring grid points, so $T_i$ is the average over the interval between $x_i - \\frac{\\delta x}{2}$ and $x_i + \\frac{\\delta x}{2}$. We can thus think of the following **piecewise constant** approximation to the original continuous function:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fc1ae03a2b0>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEGCAYAAABCa2PoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOydeXyU1b3/32cm+76RQAiQBBPWsIYdMaLiUlGqdhGs0lq4vqp9/WztrbetvXLr0vZWvd7FDbSiLWitWirUKipEFoMQICGsYUkgCWFJyL7PzPn9McMkk2SeJ5nJZCbJeb9e88o8zzmf8z3fPJN85zzPOecrpJQoFAqFQuEuBm93QKFQKBSDAxVQFAqFQtEnqICiUCgUij5BBRSFQqFQ9AkqoCgUCoWiT/Dzdge8SVxcnExOTnZJ29DQQGhoaN92yMdRPg8NhprPQ81fcN/n/fv3V0gph3U+P6QDSnJyMrm5uS5ps7OzycrK6tsO+TjK56HBUPN5qPkL7vsshDjb3Xl1y0uhUCgUfYIKKAqFQqHoE1RAUSgUCkWfMKSfoSgUiqFBW1sbpaWlNDc3dymLjIzk2LFjXuiV9+ipz0FBQSQlJeHv79+jdlVAUSgUg57S0lLCw8NJTk5GCOFQVldXR3h4uJd65h164rOUksrKSkpLS0lJSelRuz5zy0sI8UchxCUhxGEn5UII8T9CiFNCiENCiBkdym4RQpywlf1b//VaoVAMBJqbm4mNje0STBTOEUIQGxvb7ajOGT4TUID1wC0a5bcCabbXauAVACGEEXjJVj4RuFcIMdGjPVW0U7QBNiXDRoP1Z9EGAMwWSV1zG5fqmimtauR8dRMXa5upaWrDYpG6eoWir1HBpPf09nfmM7e8pJQ7hBDJGlXuBN6W1v329wghooQQI4Bk4JSU8gyAEOJdW92jnu2xgqINsHc1mBupMwdzttJI6dZ1XB4ZQFXwNMyW7lMjGIQgItiP4c37SDz7e0b7XSHaT0LjWWt7ACkr+tERhULRF/hMQOkBI4GSDseltnPdnZ/jrBEhxGqsIxwSEhLIzs52qTP19fUuawcqP97/Y4x5RvuxbGvFbJqN2ZyExRzZXvHCDjDk6DdoaQMWAYswGOow+pVjNJYiStaA3x/s1V6c9mKf+dBbhuJ1How+R0ZGUldX122Z2Wx2WjZY6Y3Pzc3NPf48DKSA0t3YS2qc7xYp5VpgLUBmZqZ0dbXoUFxd25TbRFJUEmZzIC1NI2mtacQojRixAFWOlYPGYDC0IQxtCGGiqqmaqKBYpCUAi8U2Y6ThYnt9M2AehhCx+BvPETg8Cj+/BmtZMmQlZ3newW4Yitd5MPp87Ngxpw+h1UN5bYKCgpg+fXqP6g6kgFIKjOpwnAScBwKcnFf0MXOjs7h/1tMcK6/DIiUcfxHaagDI3XeKuxaGMSagnBERgcTe9TZB/u2jmaz1WWSv/DMAbWYLlfWtnP/oFkpqJedaR/DV12eZPTsVMID/VBjzKOOGhzMvNZb/3ves1wKKQtFXrFu3jpdeegmAQ4cOMWXKFAAWL17MCy+84M2u9RkDKaB8BDxie0YyB6iRUpYLIS4DaUKIFKAM+C6w3Iv9HHSYLZKD56ooq7qGI+dr2wuG30DCpbeYFHicypbT3B1tAWMIzF4LHYJJZ/yNBoZHBjF84Spm7F1NS1srxxtDGRkQTpkpEYbfAMCJC3WcvFjP+eoI2swW/I2+NIdEoegdq1atYtWqVZSVlTF//nzy8vK83aU+x2cCihDiHSALiBNClAJPAv4AUspXgY+B24BTQCPwfVuZSQjxCPApYAT+KKU80u8ODFIq6lv49MgFLtW2MCE8w35+VEwIc2beSlJ1M+LQrxh2gwVCxsDUZ7p9oL7xro1dG7fVC8z/Fb+6+SyJIwu5eM1K9prncepSPQAWKRkbdhN/3nOWGyckMComxDOOKhT9xOHDh8nIyNCvOADxmYAipbxXp1wCDzsp+xhrwFH0EVJKDpyr4qtTlZhss7UazY3EhgVwXfowxsTatr6OWQGpK1i7Zg1rHl7jtL38i/kkRiR2LUhZASnt+gRgKVBe08TOwgrKqpuoa6kjwBjIBwdKmZ0cw9zUWAwGNQVU4Rr/9Vmhw3FrawsBAYF90vZPbkrXrVNQUMDkyZP7xJ6voe4hKLrQ3GZm86FydhRW2IOJn0FQwx5WzBnTHkw6oJdX5ve7f69Z3lk/IjKYb2UmceOEBI5WHgJASvi66Arv7y+lrrmt5w4pFD5E5xFKQ0MDDzzwAKtWrWLDhoG9DksFFIUDlfUtvLv3HKdtt5wA4iMCuXfOaCJDKjE6GRmkp+t/M9OiO70QgoykSCaMvsToDre6yqqbeGfvOS7U9HwFr0LhK3QeoXz44Yfcc889rFu3jo8++siLPXMfn7nlpfAiRRsg/1eUVDfzUcMdtMbfDFHWb1AzxkSz8Jo4p4HkKlu3bmX+/Pkud0FL7+9n4ZvTR5J7toqc05VYpKShxcxfP9/CzS3Pk87XEDLa6fMbhaIjnW9L9ee0YYvFwsmTJxk/frz9XGlpqX3EYjQ6n8wyEFAjlKGObbV74RXJ36qup7WlAco24197iFszhnNd+jB7MJkXO89pM0uXLtU088PpP9Qs19IvTV+KwSCYnRLDXTNGWqcjVxdgKtnCPy6kcqAhvX2Vvdq6ReHDnDp1iqSkJAID25/ZJCUlUVpaClgDzkBGjVCGGFnrsxxPVOyhpXUUTW2zkNIaOAyiltDWR3i/brxD1TXJa5y2m5uby8yZM52W3zf1Ps1+aekf2/qYw7HZHETD+WDMZuttsPW1Mwn29yfQ7zCi5EGIW2evm70yW9OuQtGfpKenc/So465Qd911F4888gj/+Mc/dL+Y+TpqhDLEqG6udjhuaR1NY+tspBTUF1/EaKglLHArflzsov3d8d85bbe8vFzTbpdA1gv98YrjDsdGYzNhAR/jZ7wMQH3xRZraptBsmo40tzjUzS7O1rSrUHib0NBQ3nzzTV555RVWrBjYt2zVCGWIsWz8MtZkrQEgv6SabZt+aV/tHniNie8nf0awocW6pmRZtoN22ovTnLa7evVqt/qlpR8fN77rSGNTMq3129lScy2HZCSREQATmB07jAUPfGqvtiZ7jVplr1D0E2qEMkQpKK1h2/FL1lXpBn9G+Fdw/sBOazAxhlgfcPeCtWvXutWfXuunPkOAfwB3RH3JidwD1nMGf/aGrGZv0RW3+qJQKFxDjVCGGOmx6Zy6VMcXx223tKIyGB5qYVnzL7CM0F7tnhyS7LxdnWnDk+InafdLQ9+t1tY/v/xf8Y3kc4RHNHMm8i6IymD3qQr8jIIZo6NJj3VvOrNCoeg5KqAMMRYlLeNvB8qQtv2YEyKC+OaMuwjy/xaJMdmgscvso+mPOi1LTOxmFXwHXrrtJc1yLb1TrW2VfVJUNguvXcTf885z7kojAF+euExogB/LM9S2bgpFf6FueQ0hKupbeOrTf9pXv0eH+LNseqJ9V2C9nAcvFjrPS6KnffjjbnfN6ZG+J1o/o4GlUxMZGRVsP7/1yAVe3vOOplahUPQdKqAMERpbTWw6WMaFukoAQgONfHN6EiEB7YPUu+++W7ON4sZip2V62iOXtPfr1NL3VBvgZ+COaYnEhAYAYLJIvjhaQ1VDq6ZeoVD0DSqgDAHMFsmW/HLqmk2A9R/vsmkjiQzxd6jnTpY+dzP89ZXtIH8jy6aNJCTAOuoymQ1syiujuc3sVv8UCoU+KqAMcqSUbD9+ibLqJgCEgFsnDyc+IqhL3crKSpftuKPta9uRIf7cOW0k/kbrQs3qxjb+ebgci5Mc9wrFYCMrK4vc3Nwe18/Ozub22293264KKIOcQ6U1FJTV2I9/deM9pA4L67au3loSrdzuelq9Fetaele0wyODuHnScPsalOKKRvaccS/oKRQKbVRAGcScr24i+8Rl+/H44eEI/zNO6+utBdl6cavL2j/n/1mzXEvvqjYtIZyoiAr78ddFV+yJuxQKTYo2wKZk2Giw/uyDPeIaGhr4xje+wdSpU5k8eTJ/+ctf+M1vfsOsWbOYPHkyq1evRtqmX2ZlZfGTn/yERYsWMWHCBPbt28ddd91FWloaTzzxBADFxcWMHz+eBx54gClTpnDPPffQ2NjYxe7WrVuZN28eM2bM4Fvf+hb19da/gU8++YTx48ezcOFCPvzwQ7f9Ax8LKEKIW4QQJ4QQp4QQ/9ZN+b8KIfJsr8NCCLMQIsZWViyEKLCV9XysN0hpajXzcUG5Nfc71unBN05MYMvJLU41elnkPi53nsNMT/v6wdc1y7X07miLG7cyJrZ96/tPj1zginpIr9DCtmEqjWcB2Wcbj37yySckJiaSn5/P4cOHueWWW3jkkUfYt28fhw8fpqmpiS1b2v8+AwIC2LFjBw899BB33nknL730EocPH2b9+vX227wnTpxg9erVHDp0iIiICF5++WUHmxUVFTz99NN8/vnnHDhwgMzMTF544QWam5tZtWoVmzdvZufOnVy4cMEt367iMwFFCGEEXgJuBSYC9wohJnasI6X8g5RympRyGvAL4EspZcdl0dfbyjP7reM+iJSSrUcv2B/CB/kbuX3qCN2c7GFh3d8K6wnuaD1p2/rMaAQRwdYJCK0mCx8XlGMyD+xdXRUeJP9XYO70Td/caD3vBhkZGXz++ec8/vjj7Ny5k8jISLZv386cOXPIyMhg27ZtHDnSPqPxjjvusOsmTZrEiBEjCAwMJDU1lZKSEgBGjRrFggULALjvvvvYtWuXg809e/Zw9OhRFixYwLRp03jrrbc4e/YshYWFpKSkkJaWhhCC++7T3ry1p/hMQAFmA6eklGeklK3Au8CdGvXvBdQig27Yf7aKM5cb7MdLJiUQEeSvobCSk5Pjsk13tJ62HRxgZOmUEfjZtuG/XNfCzlMVmhrFEKbxXO/O95D09HT2799PRkYGv/jFL/jNb37Dj370I95//30KCgpYtWoVzc3tSeOubnFvMBgctrs3GAyYTNYvi0I45inqfCyl5KabbiIvL4+8vDyOHj3KG2+80W3dvsCXVsqPBEo6HJcCc7qrKIQIAW4BHulwWgJbhRASeE1K2e2NdSHEamA1QEJCgsvTVevr692eKusJKpss7CozcXVCU1q0gZIj5+2/2NTGVKf9Tk11XgawbNgyl7W3ht+qWa6ld0fb0d/QejP5l6zTh4uL4cJJP0aEaX+n8tXr7EkGo8+RkZHU1dV1W2Y2mx3KQoOTMDSVdKlnCU6iwUkbPaG8vJzo6GjuvPNOjEYjGzZsQEpJYGAg5eXlvPfee9x5553U1dVhNptpaGigrq6OxsZGTCaTvY9Xy4KCgjh37hyff/45c+bM4e2332bWrFkO+smTJ7Nr1y7y8vIYO3YsjY2NlJWVMXbsWM6cOUN+fj6pqam8/fbbDjY60tzc3OPPgy8FlO7CpbN5nkuB3Z1udy2QUp4XQsQDnwkhjkspd3Rp0Bpo1gJkZmbKLI2tRrTIzs7GVa1HKNpAy8En2VA6idFiJAy/gRFjZvGtzFEO2RbrWuoID+w+O93zzz/P/fff79RExScVTn3W06bXppMY4Xx7FS29O9qO/l4nJZsPlVvTG1cXUHPun9we/g7h4cOc7l/mc9e5HxiMPh87dsxpVsYuGRun/db6zKTjbS9jCIZpv3Urs+NXX33FPffcg8FgwN/fn1deeYVNmzYxf/58kpOTmTNnDoGBgYSHh2M0GgkNDSU8PJyQkBD8/Pzstq+WhYWFMWHCBN5//31++tOfkpaWxqOPPkpISIi9TkpKCm+99RarVq2ipcWa2uHpp58mLS2NdevW8Z3vfIe4uDgWLlzI4cOHu/UvKCiI6dOn98hHXwoopcCoDsdJwHkndb9Lp9tdUsrztp+XhBB/w3oLrUtAGZTYHiLuuDKZGlMYUEPg+Q+4NaMFo8FxL6vnc563b1/fmasfOGc8ffxp7rnlHpe0yz9crjn9V0vvjrajv0IIlkxM4M/FX1FXtoUmSxufWuZyt2EbYq9t6rFKIay4+hnI/5X1NlcfpZe++eabufnmmx3OZWZm8vTTT3ep23FEkJWV5RDgr5YVFxdjMBh49dVXNfWLFy9m3759DuV1dXXccsstHD9+nL7ElwLKPiBNCJEClGENGl129hNCRALXAfd1OBcKGKSUdbb3S4Df9Euv+5luE1VV7KGtLY76lmvsp0IDvuK9LR/CHsc7f1q5QZYvd30jRXe0nrS9Pm99lyRbpounqW++1pahcjjP16YT6FcInbI9gnaWSsUgxrbxqKJ3+MxDeSmlCeszkU+BY8B7UsojQoiHhBAPdaj6TWCrlLKhw7kEYJcQIh/YC/xDSvlJf/W9P+mccRHAYobGVuvjpvriiwT4ncXfeBYsXb+5F1cXO217/fr1LvfLHa0nbTebmruc8xOlBPpZ07Basz1Ow2wJ7/b3lVed53K/FApPkpyczOHDh73dDQd8aYSClPJj4ONO517tdLweWN/p3Blgqoe75xM8df1TLB3XnndaSsnmN+7hdK11K5XT/g08l3yAIAPdZl389bZfO21bKyc8QEak8/UeetprR1+rWa6ld0e7asYqnlr8lOPJTcmYGgp450oiX5a0MSncj8SAuXxr5EkM38x2qPrs+89q2lYoFO34zAhF0TM6BhOA4xfqOB25HAzWacHTQk4QZGh1mnWxyz/XXvBgyoMua92x2+faqc/g5xfEzRE5WCcFwnlTIgeG/7ZL1flx8122rVAMNVRAGWC8lvua/X1jq4kvCy9DVAaMXMrU6GpKS6usI5PZa7u9B7ziA+f3hffv369p++ljXR8e9lSrZVdP3+falBUwey3xkZHIiyfAPxJGLiWncSaV9Y63vT46/5GmbYVC0Y5P3fJS6FNeX25//+WJyzS1WtdVhA+fzsJ5m0m9/hwkJzvVl9WVOS1buXKlpu2KFueLAfW0Wnb19B7R2h66/mJKEV9fNnKxthmTRfLZ0Yt8Z9Yo+6KvK60qP71C0VPUCGWAcuZyPccvtC9CunFCAgF+BjZu3Ohym97SetP2X959hyWTEuxrdcprmjlUWqOjUih8j+LiYiZPnuzVPqiAMsAIDwinxWRm2/FL9nMTRoSTHBcK4LBFgzO9M/S0IcYQp2Xu2NXTe1obFxZIZnK0/dyuUxXUNbcB2j4rFApHVEAZYDw2/zF2n6qwb/wYHGDkuvR4e7lekpzNyzc7LdPTPpvhfMaTO3b19P2hnZ0cQ3RI+waSV7f9//aob2vqFYOTQxcOsSZ7DT/4+w9Yk72GQxcOudXevn37mDJlCs3NzTQ0NDBp0qQuU36Li4uZMGECq1atYtKkSSxZsoSmJmtivLy8PObOncuUKVP45je/SVVVFWB9fjh16lTmzZvHSy+9ZG/LbDbzr//6r8yaNYspU6bw2mvWZ6/l5eUsWrSIBQsWMHnyZHbu3OmWX51RAWWA8cHhbeSXtN+SyRo3jGBbuluAd97R3i/zxT3Ok2Tpad8vfd9lrZZdPX1/aP2MBm6YkGA/f+pSPacu1at1KEOQQxcO8VzOc1Q1VZEUkURVUxXP5TznVlCZNWsWd9xxB0888QQ///nPue+++7q9PXXy5Ekefvhhjhw5QlRUFB988AEA999/P7///e85dOgQGRkZ/Md//AcA3//+9/mf//mfLhukvvHGG0RGRrJv3z727dvHunXrKCoqYuPGjdx8883s3r2b/Px8pk2b5rJP3aECygDCYpH89cAx+3HqsFDGJTje0pk3b55mG5uOb3JapqfdVbHLaZk7dvX0/aUdFRPCpMQI+3H2iUscuOLeN1PFwOPD4x8SHRRNdHA0BmEgOjia6KBoPjzuXhKqf//3f+ezzz4jNzeXn//8593WSUlJsf+TnzlzJsXFxdTU1FBdXc11110HwAMPPMCOHTu6nP/e975nb2fr1q28/fbbTJs2jTlz5lBZWcnJkyeZNWsWb775Js8++ywFBQVu7U3WHSqgDCAOldXQ1GK9LeNnEGSNi++yBfXVbGyu4C2tN2131i5KH0aIbcRX12yiqnGYy20rBibnas4RGRTpcC4yKJJzNe5tX3/lyhXq6+upq6tz2Ka+Ix2fBxqNRvs29d0hpXS6Bb2Ukv/93/+1b1tfVFTEkiVLWLRoETt27CAxMZHvfe97vP3222751BkVUAYIDS0mvjrdPm13dkoMkcFdc5wUFBS4bMNbWm/a7qwN8jdy3bj2IFLXHENFvfbGl4rBxejI0dQ0O870q2muYXTkaLfaXb16NU899RQrVqzg8ccf77EuMjKS6Oho+/OOP/3pT1x33XVERUURGRlpT6q1YUN7Rsmbb76ZV155hbY26+SSwsJCGhoaOHv2LPHx8axcuZIHH3yQAwcOuOVTZ9Q6lAHCzpMVtLRZmJk4k6gQf2aOie623urVqzXbeeHmF5yW6Wl/NPZHLmu17Orp+1s7LiGcgtIaSquaSAsbR/aJy9w9Y6RHEhIpfI+7xt/FcznPAdaRSU1zDVXNVTw43fWdIt5++238/PxYvnw5ZrOZ+fPns23bNhYvXtwj/VtvvcVDDz1EY2MjqampvPnmmwC8+eab/OAHPyAkJMRhJ+Mf/vCHFBcXM2PGDKSUDBs2jE2bNpGdnc0f/vAHjEYjERERfT5CQUo5ZF8zZ86UrrJ9+3aXtb2l5EqDfGHrCfnC1hPyP7bkyqLL9U7rPvnkk5pt7T+/32Xtax+95hG7enpvaC/XNcsXPyuUP3zlb/KFrSfkiQu1mu0MJvrzs91fHD161GlZbW3Xa5tfni+f3P6k/P6m78sntz8p88vzPdm9fqc7n53R3e8OyJXd/E9Vt7x8HItFst02hRXgdO0u+5qT7oiNjdVs76ef/tRl7cunX/aIXT29N7RxYYFMHRVJYV0hADsKL9NqUnnohwpThk9hTdYa/njnH1mTtYYpw6d4u0sDAhVQfJz80moq6qz38P2NgpFxtZr13cm05y2tN21raeemxmI0WB+K1jWb2FestmFRKLRQAcWHaW4zs+dM+z+xOamxBPibNTVX5627gre03rStpQ3yNxIdetF+vP9sFVUNrS7bUngX650aRW/o7e9MBRQfJudMJc1t1gASGezP9FFRmhkXQf/b+rLxy1zWLoxb6BG7enpvaQEWJqSQGGXNNWO2SHacvKxZX+GbBAUFUVlZqYJKL5BSUllZSVBQUI81PjXLSwhxC/DfgBF4XUr5u07lWcDfgSLbqQ+llL/piXagcaWhlUMdVsQvSo/Dz2jQDSjnz5/XLH907qMua+9J6j6fvLt29fTe0gJMj57GxHHxbNx7DinhzOUGzlU2MjpW7fE1kEhKSqK0tJTLl7t+IWhubu7VP83BQE99DgoKIikpqcft+kxAEUIYgZeAm4BSYJ8Q4iMp5dFOVXdKKW93UTtg2HnyMhbbt6mk6GDGDgsD4Pmvnuex+Y851RUWFmq2u3TjUqf7W+lpf1nwS77K+qrP7erpvaUFeK/kPV7OymLiiAiOnLc+v/ry5GVWRI/GYFDTiAcK/v7+pKSkdFuWnZ3N9OnT+7lH3sVTPvvSLa/ZwCkp5RkpZSvwLnBnP2h9jrOVDZy53ACAEHBd+jD7Goi61jotqe56EC29nrbR3OgRu3p6b2mh3ef518QR4Gf9c6moa+FoufbkCIViKOIzIxRgJFDS4bgUmNNNvXlCiHzgPPAzKeWRXmgRQqwGVgMkJCSQnZ3tUmfr6+td1mphkZLtJSZqW6yjk+RIA0cPnOfqUKu4uFjT7vr16zUTTlVXVzvV62nNZrPLWi27enpvaQFaW1vtemOVmeJK6zOtP5We5aYxfvgNwlGKpz7bvspQ8xc857MvBZTu/jI7P0E7AIyRUtYLIW4DNgFpPdRaT0q5FlgLkJmZKV2dcpqdne32VNkuFG0gf/drxFSMIsY/koDExay8fRmhge2X6UTuCbIynds9ceKEZr8mVU5yWq6nTTiW4LJWy66e3ltasKYAvlo+32Th7Zxi6i4chAtfEHJqN/Pja2DqM92mWx6oeOSz7cMMNX/Bcz770i2vUmBUh+MkrKMQO1LKWillve39x4C/ECKuJ1qfp2gDLXt+RE6FLbdJWw2zqn9L6Pm/OFT7l8x/0WwmMzNTs3zD3Ruclulpn5jwhEfs6um9pQW4I/EO+/sAPwPzgvZC2WZoq+FAw3jq6i7B3tVQpG1HoRgK+NIIZR+QJoRIAcqA7wLLO1YQQgwHLkoppRBiNtaAWAlU62l9iaz1WV1PVuyhqWUczW3W3UYNhno2teUhPngQ4tbZqz027zGWjlvqtO3Nmzczc+ZMp+W/3vZrnlr8lEvaN4recPqtxh27enpvaQF+sO8HjC5u3xRQXt5DfdP1mCwxgB9vNUwjNCAHShyvE0D2ymyn7SoUgxGfGaFIKU3AI8CnwDHgPSnlESHEQ0KIh2zV7gEO256h/A/wXdvWMt1q+9+LnlHdXN3lnMVsoMU03n4c7JePEBawOO50+8Ex7UV8S5Ys0Szfec55hjY9bUGN81193bGrp/eWFqC6zfFaCdlCkH/7Dq2tphRrcLF03ZE4uzhbs22FYrDhSyOUq7exPu507tUO7/8P+L+ean2VZeOXsSZrjcO5L968nUNV1stx7EABr916FiGAkDGwLNter9vRTQcKCwuZP3++S/3yltabtvW0o0NGO440NiVD41k+qi7lnR2tzJ6dypjAqdw18ozDdQJYk71Gd92QQjGY8JkRylDmSkMrh8NXgsGa3ySk7ow1mBhDrA98e0FxcbHL/fCW1pu2e62d+gwYQ1gYlkdVlTU519m20ZxL+a3LfVAoBgsqoHiBx+Y5Lkz86nQFlsgpMHIpo8Pb+NebGqwjk9lru8we2njXRs229dZVaOn1tE+Md/5Q3h27enpvaaEbn1NWwOy1xEREszLLH/wjYeRSdjUv6LKtR+frrFAMdlRA8QLn69onoJXXNHHyoi0NbVQGC7/zNusqnoRlxd1ORc2/mK/Z9tq1azXLtfR62tMNpz1iV0/vLS048TllBSwrJr/hG/hN/AlEZXCxtplTlxzTCXe8zgrFUEAFFC/wzhxXDycAACAASURBVOF3AOvmaztPtqf1HTc8nISIIJKTk51qf7/795pta2n19Hrad0re8YhdPb23tKDt87hrxjJtdJT9ePepCiyW9lHK1eusUAwVVEDxIkUVDZRVNQFgEIL5Y63JntLT011ucyBqvWnbXe2s5BgC/a1/RlWNbfb9vhSKoYgKKF7CYpHsPtU+OpmSFElUSAAAW7dudbndgaj1pm13tUH+RjLHxNjP7TlTSZtZZXZUDE1UQPECS9OXcuJiHRX11mRNAX4GZqe0/1NautT5wsUfTv+hdtsaWj29nva2Ebd5xK6e3lta6JnP00dHEWbbHqe+xUR+iXXtytJ07bYVisGGCiheYNrwGew5U2k/nj46ymG/rtzcXKfa+6bep9m2llZPr6ddkuB8EaA7dvX03tJCz3z2NxqYk9r+hWBv8RWa28zMTHS+Al+hGIyogOIFHt38B6ob2wBrmtkZo6MdysvLy51q9RY2amn19HraR/OcJ6tyx66e3lta6LnPkxIjiQqxriNqabOQW1zFmuw1mm0rFIMNFVD6GZPZwoWqcPvxzDHRBPkbHerorY3QYiBqvWm7r7RGg2D+2Dj7cV5JFW0m9eelGFqoT3w/c/h8La1t1gASEmBk2qioLnX01kZoMRC13rTdl9r0hDDiwq2be7aZJZeqw1xuW6EYiKiA0o+0mS3sLaokNsQ6PTgzOcaeBbAjWlNZJ8VP0rShNw1WS6+nTQ5J9ohdPb23tNA7n4UQzEuNtR/L1jHUt5g021coBhMqoPQjh0qraWgxkxGfQVigH1OSIrutl5iY6LSNl257SdOGllZPr6d9NN358wR37OrpvaWF3vs8dlgoCRFBAEwalsG+4iua7SsUgwkVUPqJFpOZfcVVABRcKmB2Sgz+xu5//VqpOR/++GFNO3ppPbX0etoXC1/0iF09vbe00HufhRDMsy1QLbhUQEFpDbXNbZo2FIrBggoo/UTeuWqaWq35yOvaLjF5ZPejE4C7777badmRS9ppXrS0eno9bXFjsUfs6um9pQXXfE6ODWFEZBCVjZWYLZJ9RWqUohga+FRAEULcIoQ4IYQ4JYT4t27KVwghDtleXwkhpnYoKxZCFAgh8oQQ2osL+pnmNjP7z1XZj4fH1GE0CKf19b41azEQtd607QmtEI4zvg6X1VLTqEYpisGPzwQUIYQReAm4FZgI3CuEmNipWhFwnZRyCvAU0HmKzvVSymlSSu1E4f3MgbNVtLRZt+OIDvEnJrxRs35lZaVm+WDTetO2p7SjYoIJC7buhGCRkq+L3Pv9KBQDAimlT7yAecCnHY5/AfxCo340UNbhuBiI643NmTNnSlfZvn17j+o1tpjk/207KV/YekK+sPWEPFZeo6spKytzuV+e1Gr57I5dd/W+6nPJlQb7dX/xs0J5pb7FlS72Oz39bA8Whpq/UrrvM5Aru/mf6kspgEcCJR2OS4E5GvUfBP7Z4VgCW4UQEnhNStntAgMhxGpgNUBCQoLLtzzq6+t7pD1cYeZklfXZSUSgoNyvjC/3FTIufJxTzfr161m5cmW3ZVsvbtXcDkRLq6fX024+uxmy+96unt5bWnDP5xN1J2isTeVSo3VL+zcqzpI53Jf+5Lqnp5/twcJQ8xc86HN3UcYbL+BbwOsdjr8H/K+TutcDx4DYDucSbT/jgXxgkZ5NT49QGlra5P9+UWj/lnryYq2UUsontz+pqXv//fedll335nUua/X0etqp/zXVI3b19N7SSumez09uf1KWVTXar/9/fXZCVtQ1a2p8gaH2jX2o+Sul50YoPvMMBeuIZFSH4ySgS8o7IcQU4HXgTiml/ca0lPK87ecl4G/AbI/2tgccOFtNm9n67XRYeCBjh/Vs5XRYmOsrrAei1pu2Pa1NjAomJS4UAClhzxk140sxePGlgLIPSBNCpAghAoDvAh91rCCEGA18CHxPSlnY4XyoECL86ntgCXC433reDU2tZvJLq+3Hc1JiEML5zK6O5OTkuGx3IGq9abs/tFfXpQCcvFRHZX2LyzYVCl/GZwKKlNIEPAJ8ivV21ntSyiNCiIeEEA/Zqv07EAu83Gl6cAKwSwiRD+wF/iGl/KSfXXDgwLkqWk3WmV1xYQFcE9/+bfbeyfdqau+913n54wsed1mrp9fT3jvKebk7dvX03tKCez5fvc4JEUGkDmsfpXyt1qUoBik+E1AApJQfSynTpZRjpZTP2M69KqV81fb+h1LKaGmdGmyfHiylPCOlnGp7Tbqq9RbNbWbySjqMTlJjHUYnieHa231s2bLFadnUhKlOy/S0eno97djQsR6xq6f3lhbc87njdZ6T0j5KKbyoRimKwYlPBZTBwoGz7aOT2LAA0uId77U/n/O8pr6lxfk/m+UfLndZq6fX0z59/GmP2NXTe0sL7vnc8ToPjwxyeJai9vhSDEZ8fw7jAKP55AYObv0SWhrAP5I5869FiORetbF8ufY/wcGm9abt/tTOSY2hqKIBqgs4fuIL5hzYQHRENEx9BlJWuNwPhcJXUCOUvqRoAwe/fJHWlgYAYmQJaSdXQdGGXjWzfv16l7swELXetN2f2hGRwSRb8qBsM7K1hq8bJkHjWdi7utefEYXCF1EjFBd5NO9Rooodk2NZLu+ntvEWrMthIDSggDebG6D0QYhbZ6/32LzHNNueOdN5LvJrR1/rslZPr6fNiMzwiF09vbe04J7PHxz9gOzibIdzposnqWvKAkDUJfNMzWGMhjoocfyMZK901CkUAwE1QulDWluTkTIAgMazZ/A3nrMWWBzvtS8dt9RlG08tfsplrbv6B1Me9Ipdb2nBPZ+vJlLriJ84j7/Rmou+rugSzSZbAjCLekivGPiogOIii+MXk70y2/76ZMUX3B81l5URsDICllQe58tRkuwkyE4f41B3xQfa98v379/vtMwdrZ5eT/v0MecPqN2xq6f35u/LHZ9Hho90uO7ZK7PJTh/D30cdZmUEXFN+kftCk9k0PKzLZ+S13Nc021YofBEVUFzkSqvjLJ38kmpa4m8Bgz9RfnU8ttg2bdgYYn3o2oGyujLNtrX2h3JHq6fX01a0VHjErp7em7+vPvd56jOMDG5gVMAFZs9KxSIN7G2a0eUzUl5frtm2QuGLqIDSB7SYzBw4Vw1RGTByKbNjr/BuDhAyBmav7fUMno0bN7rcl4Go9abtftemrIDZa5k7rIoDB86CfyTHYn9MTcK3Xe6HQuErqIDiIiHGEPv7/JIamtusOwpHjpjBhHt3Ejj7OVhW3G0wCQ8I12w7MDDQaZk7Wj29nrajz31pV0/vzd+XR3xOWUHSdw4Qm74Uxj+KJXIKezutS9Hrt0Lhi6iA4iLfHmX9RtlqsnCgQzbG2SkxGAyC22+/3al28/LNmm17Squn19M+m/GsR+zq6b35+/Kkzz+4tz2F8NHztdQ0tWd1fGy+9kxAhcIXUQHFRfKq8wA4VNqeKz4i2J8JIyIAeOedd5xqX9zzombbntLq6fW075e+7xG7enpv/r486XP2PzcxMjoYsGZ17Jh7vvN0Y4ViIKACiovkVefRarKw/2yH0UlyjD1X/Lx585xqNx3fpNm2p7R6ej3tropdHrGrp/fm78vTPs9LbZ9afLS8ltpm6yhFBRTFQEQFFDcoKKum0TY6CQ/yY2JihL2svr7e5XaHmtabtr3tc1J0MCOjrKMUs8VxlKJQDDRUQHERixTkFjs+O7k6OgEoKChwue2hpvWmbW/7LIRgbodRypHz7aMUhWLA0V0ax6HycicF8B/+8oE9teu6Hadlm8nsUF5WVuZUu//8fs22PaXV0+tpX/voNY/Y1dN78/fVHz5bLBb57t6z9s/TtmMXZVmtdtueZKilxB1q/ko5NFIADxjazBaKa9p/dbOSY/AzOv4q165d63L7Q03rTdu+4LMQwiFfSkFZDQ0tZrfaVii8Qa83hxRCLAfuAMyAADZLKbWnu/S87VuA/waMwOtSyt91Khe28tuARmCllPJAT7R9SUFZDQXVhYxlIuFBfkzq8Ozk0IVDfHj8Qz4p+wSy4a7xdzFl+BQH/U8//Wm3m/95WutM31Pty6dfZjWr+8xuf/js7u+rv3zOSMhgRGQQ5TXNmC2S329/n9e/9ROn2nM15xgdOdqp7e4YiNqB2u+B6rO7uDJCuU5K+V0p5Qop5XJgYV90RAhhBF4CbgUmAvcKISZ2qnYrkGZ7rQZe6YW2TzCZLezv8Oxk5pho++jk0IVDPJfzHFVNVUydO5Wqpiqey3mOQxcO6bY71LQDtd+e8rngYgFzOjxLqawNob7F5FSbFJHkcr8Hinag9tvXfe643qmvcSWgBAohviGEmCKEuA0I7qO+zAZOSWs631bgXeDOTnXuBN623cbbA0QJIUb0UNsnHD5fa/9DDwv0I2NkpL3sw+MfEh0UTXRwNCd2niA6OJrooGg+PP6hbrtDTTtQ++1Jn5NjQ0iICALAYhHkdlo931FrEAaX+z1QtAO1377sc0OLibe/KmZXWRslVxp71J/eIOTV5B3OKgjxnJTyZx2OQ4C7gCSgBPiblNLtngkh7gFukVL+0Hb8PWCOlPKRDnW2AL+TUu6yHX8BPA4k62k7tLEa6+iGhISEme+++26P+yil5ItzJupaJQUlBUSWl7F03mTy8vKora1l2nXT2Ll1J6OTR3PxwkVaWlqYd+08cvfmEmYM48Ybb+STTz7h8tjLLDAuoLCwkFtuuYXPP/+cenM9mbMzydmZQ2BgIAnDEzhXfI5rl1xL3pd5REREMG3aNF7Ne5XlqcupqKigtLSU22+/nT+//2diY2NJTUvl082fcsMtN3C+7DwXL1zkvnvuY8uWLQwfPpykpCTeLnybByc9yIkTJ7hy5Yq9zyNGjqCmpobGhkZmzZvF4UOH8Wvz4/bbb2fLli2kpKSQQw7Di4azePFi9uzZw6XaS1ybdS07tu3Az8+PUcmjKDpZxLzr53F873ECAgKYO3cu27ZtozKtkkxTJkVFRfY2Tf4mJk+ZzL6cfYSEhhAZGUl5Wbnd55iYGMaNG8cbR97g/vT7KS0t5cKFC3afE4YnkDgykS8++YKbl97MmZNnqKystPuclJTEkfAjRB6LZNGiRb2+Tunp6WQ3Z5N4LrHX12nHjh3UTKhhUt0k+3XasmULhEBaWhoH9h0gPCKcoKAgLl+8bNeHDE+hKnws+VfyyRw2nahLB6m9ctnu88ikkcTExvDl51/yjW9+g6MFR6mpqbH7nJKSQq5fLrEnY+3XqbW1lfGzx5OzPYeUtBRKikswmUwsWryIndk7iY+It1+niRMn8mntp4wqHdXr65STk0PdpDrSKtN6fZ3i4uLYeGYjD017yH6druqTU5IJCg7iqy+/4s5v3cmBvQdoaW3hntvvsV+n3ebdDDs9zH6dAgICGD15NPt27iNtQhqnjp9CIlm8ZDGfffIZo+NH26/T5bGXmd4yvVfX6erfU25uLg0ZDYwqG8WVK1d6dZ2Cg4N5v/R9Hsl8xH6dbrzxRt7d9C7XpF0DAvbu3ss9y+9hd/ZuhJ/gtutvs1+nbU3bGF403N7niIgI4sfGc3DPQSZOmcjhvMMYjAa7z6mJqfiPmsKe05fwi/Xj2msyuC7JD+uThN5x/fXX75dSZnb5/9qDgHJASjmj1xZ7iRDiW8DNnYLCbCnljzvU+Qfw204B5edAqp62OzIzM2Vubm6v+tnUaubAuSo+3n2QNffdgH+Hh/FrstdQ1VRFdHC0/dzV4zVZazTbHQja7OxssrKy+sRuf/bbHW1/+yylZOPec7x36FOykrOYOSaaRenDumiz12eTtTKrW9tZ67O6PLvpqRZg2ovTyHs0zyVtZ9vuaPvL587+Dlaffz7/Cf64q4g2s2RT3ib+uPJHjB0WhisIIboNKL40y6sUGNXhOAk438M6PdH2CcEBRhZcE8cNo/0cgglYH65WNVdR1VSFRVqoaqqiqrmKu8bfpdvuUNMO1H572ufOM74OlVbT2Grqoh137TiX+z1QtAO1377q84Gz1bSZrQOIAL9mUuNCe9Sf3tCTgDJVCFEkhPhICPGsEOJeIUSGEMK/j/uyD0gTQqQIIQKA7wIfdarzEXC/sDIXqJFSlvdQ26cYuhkmThk+hZ/N+xnRwdGU1pYSHRzNz+b9rEczLIaadqD2uz98HjsslOBA64PTNrO0b+/TUZu/J9/lfg8U7UDtty/6nBYzifzSanu9qJDLLt3q0qW7xSkdX8BBIBFYAjwGvAnkAjXAYT19b15YpwMXAqeBX9nOPQQ8ZHsvsM7mOg0UAJlaWr2XOwsb1WKooYG3fN595rR9oeP/bTspG1tMDuVPPvmkU63egkotrZTaizn1tFq23dHq6d3Ravmrp9Wz7Ss+7zp52f55evurIvneP/+qqdUDJwsbe7QORUp5HustpK1Xz9nWhFzjeijr1s7HwMedzr3a4b0EHu6pVqEYqIyJDeT0RSMVdS32TUgXpsXZy1evXq2h1mYgar1pe6D73NxmJq+kfXQyJzWWQ3nOVO7Rk1teL3V30haoTvZxfxQKBbDuwDrmpsTYj/M7pEkA7VX6P/30p5pt663wf/n0yy5rtWy7o9XTu6PV8ldPq2fbF3w+cK6KVpMFgNiwANLiw9hSvkVT6yq6AUVK+bpHLCsUCk2uiQ8jLiwAsCZyO9ghkVtGRobL7Q5ErTdtD2Sfm9vMHDzXPjqZnRLjmWcnNnxplpdCoeiAEILZHWZ8HSyptqeaDgtzbbrnQNV60/ZA9vnguWr76CQmNID0eM+mllYBRaHwQbKSswBIiw8jtsMo5Wq66ZycHKfaZeOXabatpQVYGOd8NyU9rZZtd7R6ene0Wv7qafVse9Pnnbt2c7Cka3pygGlR0zS1LtPdk/qh8lKzvHqH8tk7HCuvsc/QeWn7SdnUapLHjx93uT09rZbPnrTrSb2WVu8aD1Sf//JZjv1z8+auM9JsttjL1Pb1CsUQ4vmvnre/T48PJzrEuuyrpc3CwXPV1u1BnLB041LNtrW0AL8s+KXLWi3b7mj19O5otfzV0+rZ9pbPLSYzf/nw7/bjWR1GJwDvlbynaddVVEBRKHyQutY6+3uDofOzlCoaGpt6pO2OlpYWzfJGs/Ot+fS0Wrbd0erp3dFq+aun1bPtLZ/zS2pobrZqI4P9mTA8wqFcz2dXUQFFoRgAjB8eTlSIP1QX0FLwAhnGP8GmZCja0Ou2li9f7nI/vKX1pu2B5nPrqQ0c+OfPmZF4Go6/yOzArx1GJ55EBRSFwgcZETbC4dhgEMwO+BrKNkNbDa/vMNJSXwZ7V3cJKiPDR2q2vX79es3yuMA4p2V6Wi3b7mj19O5otfzV0+rZ7nefizZw6MvnaWpuYu++M0RYyphw+oddPiMxATFdtX1ArzM2KhQKz/PO4Xd457BjIlR5+WvqmpZgtoRRGz2M684aCPI/AiUPQtw6e73uskR2ZObMmZrlT0x4wmXthrudj5jc0erp3dFq+aun1bPtSZ/L6srIWp/lcE5ezqW28WYsEmpjYni94Qh/OtcApY6fkTXJazTtuooaoSgUPojJYupyTshmAv0OAxCanECLaTxS+oHF8T77e4fde+D6RtEbLmt/ve3XA07rjr/u2nZHe6n+UpdzLa2jsEhrkrbw1DACjEXWgk6fkUaTZ56hqBGKQuGD7PrBrq4nNyVjbihmfcVkas3WBW8Lw9OZNawBlmXbq63JXqPZ9v79+1m61PkMoYKaApe1O8/t9IhWT++OVstfPa2ebU/6/O3J33bIj9JmtvDH11bT2GydsHFDRCFTQqyLGgkZ4/AZyc5uf9+XqBGKQjFQmPoMRr8gZocesZ/a3zSF1snP9KqZlStXutwFb2m9aXug+FxQVkPjsFvB4E+4sZGJwWesBcYQmNq7z4irqICiUAwUUlbA7LVMjDUTbmwE/0iaRtxDgf+tvWpm48aNLnfBW1pv2h4IPreZLeQWX4GoDBi5lMzYK/gJaR2ZzF5r/ez0A+qWl0IxkEhZgTFlBbNLq/nimPUeem5xFVOSouwZRB+b95hmE4GBgZrlT4x3/pBaT7vxLuf/AN3R6und0Wr5q6fVs+1Jnzte54KyGhpabPu8JUxn8oLtYOz/8YIaoSgUA5CJIyIID7J+H2xsNXOotMZedr5OO/v17bffrll+uuG0y9r8i/ke0erp3dFq+aun1bPtSZ+vXudWk4V9RVfs5zOTo/HzQjABHwkoQogYIcRnQoiTtp/R3dQZJYTYLoQ4JoQ4IoT4fx3K1gghyoQQebbXbf3rgULRv/gZDcxKbl9LsP/sFdrM1gewnacbd+add3TKS5yX62l/v/v3HtHq6d3Ravmrp9Wz7Umfr17nQ6XVNNpy5YQH+ZExMlKzTU/iEwEF+DfgCyllGvCF7bgzJuAxKeUEYC7wsBBiYofy/5JSTrO9VOZGxaBnUmIEYYHWUUpDi5mCshodhZV58+a5bNNbWm/a9mWfW0xmcs867ijsrdEJ+E5AuRN4y/b+LaDLvsxSynIp5QHb+zrgGKC9zFShGMT4GQ1kJrcP5vcXV9lHKVrU19e7bNNbWm/a9mWf80tq7Jk8I4L9mZTovdEJ+M5D+QQpZTlYA4cQIl6rshAiGZgOfN3h9CNCiPuBXKwjmapupAghVgOrARISElyej11fX++xudy+ivLZ9zBbJBdKTTSbJABvVRYxwm+EZp83b95MbGys0/Ibom9wqtfTzg+e7xGtnt4drZa/elo92570eVjNCP66PZermaGnxxvZueOs07Y64rHPdXd72nviBXwOHO7mdSdQ3alulUY7YcB+4K4O5xIAI9YR1zPAH3vSJ5UPpXcon32Tg+eq7HkvXs0+JVvazJr1y8rKNMu1fNbTumPXk3otrd419lWfc05X2K/7H3edkaYO+U70GPD5UKSUN0opJ3fz+jtwUQgxAsD2s+ueAtYyf+ADYIOU8sMObV+UUpqllBZgHTDb8x4pFL7B5ETHGV//76MXNOuvXbtWs/zRvEdd1nbeW6qvtHp6d7Ra/upp9Wx7yufmNjP/mf1X+/GclFiM/bSjsBa+8gzlI+AB2/sHgL93riCEEMAbwDEp5QudyjpuzfpNrCMfhWJI4Gc0MDe1/bbIxeowe+757tC6BaOHt7TetO2LPh84W4XZbP33HR3iz/jhns0V31N8JaD8DrhJCHESuMl2jBAiUQhxdcbWAuB7wOJupgf/pxCiQAhxCLge+Ek/91+h8CoTRkRY86UAZrOBA2e7fYQIQFZWlst2vKX1pm1f87mp1czBkmr78dyxsf2W70QPnwgoUspKKeUNUso0288rtvPnpZS32d7vklIKKeUU2Wl6sJTye1LKDFvZHdL2gF+hGCoYDYJ5Y63fZmNDYjlYUk1ja9cdiwE++OADzbaSQ5KdlulpJ8VP8ohWT++OVstfPa2ebU/4vP9sFa0mC7EhscSGBZAe7xujE/CRgKJQKNwnPT6cuLAAMuIzaDVZyC3ufpSi96350XTnzxT0tC/d9pJHtHp6d7Ra/upp9Wz3tc+NrSbySqzXNSM+g3mpvjM6ARVQFIpBg8E2Sim4ZN2OPb+kmvqWrqOU8+e1t2Z5sfBFp2V62oc/ftgjWj29O1otf/W0erb72ud9xVW0ma1TxItqD3FNfJimvr9RAUWhGESMHRZGk/kCACaLZG9RZZc6hYWFmm0UNxY7LdPTHrl0xGmZO1o9vTtaLX/1tHq2+9Ln2uY2DnV4dmIIPIV1rpLvoAKKQjGIEEIwIrbWfny4rJaapjaHOqtXr3a5fW9pvWnbV3zec7oSk8U6OhkRGURkaLNbbXsCFVAUikFGeEgLI6OCAetK+q/POI5S9NZGaOEtrTdt+4LPlfUtHC1v/6Kw4Jo4fGxwAqiAolAMOv7j+jX2GV8Ax8rrqGpotR+PGDGiO5mdF6c5f6agp81eme0RrZ7eHa2Wv3paPdt95fNXpyuxbgoCY2JDGBUT4pD+11dQAUWhGGTsP7+fUTEhjI4JAcAiJV+dbh+lZGZmauq3XtzqtExP++f8P3tEq6d3R6vlr55Wz3Zf+HyhpplTl9o3iVxwTRxgvc6+hgooCsUgY3PhZqD9Hw9A4cU6LtRY77lv3rxZU/9xufPsD3ra1w++7hGtnt4drZa/elo92+76/NFHH7HrVIX9OD0hnISIIGvbhdptewMVUBSKQcrwyCDSEtqnle48eRkpJUuWLHG5TW9pvWnbmz5PnbuIkiuNABiEYP5Y97Zy8TQqoCgUg5gFY+Mw2J7ellY1UVTRoDuVVQtvab1p21taKSWffpVnP56UGEF0aIDL7fUHKqAoFIOMeyffa38fHRpARlKE/Xj3qQqKioq09aPudVpWXFysqX18weMe0erp3dFq+aun1bPtjvbkpXq73s8gmJMa41De8Tr7CiqgKBSDjMTwRIfjOSmxBPgZoLqAitz/5LqIt2BTMhRt6FY/NnSs07b11lVMTZjqEa2e3h2tlr96Wj3brmrNpzew+8MfMy+1FI6/yLSAA4QH+TvU6XydfQEVUBSKQcbzOc87HIcG+jEjYD+UbYa2Gl78LIi2hlLYu7rboPL08aedtq23rmL5h8s9otXTu6PV8ldPq2fbJW3RBg59+Z9UN7WRs+cUgebLzCr9ly7XqvN19gV8JQWwQqHoI9bnrSe7ONvhnLycS23jEiwymCtB0VxbPI4g/6NQ8iDEretx28nJyS73yx2tN217Upt3Ia9Lki3L5QPUNS7BIuFKcBg1DYd5r6QGyhyvVVayo84XUCMUhWKQERUU1eWckA0E+VvzzsVMS6XZNBGLDABLS5e6w4OGO207PT3d5X65o/WmbU9qg/yCupxrab0GiwwEYNiMeAL9bA/2O10rXwwoPjFCEULEAH8BkoFi4NtSyi57bwshioE6wAyYpJSZvdErFEOBvIfyUU+4YQAAGa9JREFUup7clIy54TR/qhxHlSkCCGBG6GSuS6iEZdkOVZ99/1mnbW/dupX58+c7Lb929LUe0erp3dFmRGa4rNWzraddt3QdS8cttR/XNLbx1h9/jLm1DoBvROWRHmSxFoaM6XKtfA1fGaH8G/CFlDIN+MJ27Izrbcm1Oi5B7Y1eoRh6TH0Go18QC8Ly7afymiZRNa5r8Jgf5/wf4NKlS52WATy1+CmPaPX07mgfTHnQZa2ebT1tx2ACsPt0BeaEm8Dgzwj/CtICz1kLjCEw9RnNtnwBXwkodwJv2d6/BSzrZ71CMbhJWQGz13JNtIGRAZfBPxLLyDvZ0XJdl6ofnf/IaTO5ubmaZlZ8sMIjWj29O9qnj2k/lPekz6/lvmZ/X17TxIkLdRCVASOXsij+onV7+pAxMHut9Rr6OD5xywtIuJq2V0pZLoSId1JPAluFEBJ4TUq5tpd6hBCrgdUACQkJZGdnu9Th+vp6l7UDFeXzQGckRK/HP9jC2VITshqK9x2h5fwJ4kPav1teqL/g1OecnBzGjRvn1MKRkiMe0erp3dFebLzoslbPtp42pziHcfXjkFKyo8zElSbrDpAjwyMpHP4i9mWRZ4GzzvvYWzz1ue63gCKE+Bzo7mnfr3rRzAIp5XlbwPhMCHFcSrmjN/2wBaG1AJmZmVIvRaczsrOzddN7DjaUz4OH0CMXOHLeuh16S1gAi+aMsaeSXV+83qnP6enpJCY6X/8QVRzlEa2e3h2tMc/oslbPtp726uer8GIdEaZyIgCjQXD/vDFEhXhuVbynPtf9dstLSnmjlHJyN6+/AxeFECMAbD8vOWnjvO3nJeBvwGxbUY/0CoXCyvxr4qyLHYGK+lYKymrsZSHGEKc6vXUV4QHhHtHq6d3Ravmrp9Wz3RNtq8nCjsLL9nNTR0V5NJh4El95hvIR8IDt/QPA3ztXEEKECiHCr74HlgCHe6pXKBTthAX6MSu5fSuPnDOVNLeZAfj2qG871elNg9283PkOuO5o9fTuaJ/NcD6rTU+rZ1tP+9j8x8gtvkJdswmAkAAjc1JiNDW+jK8ElN8BNwkhTgI32Y4RQiQKIa7uLZ0A7BJC5AN7gX9IKT/R0isUCufMGB1FRLB1O4+mVjN7bJkd86q7mXZsQ+v2DcCLe5wnq3JHq6d3R/t+6fsua/Vs62n/cWI7+8+2r3BYcE0cQf5GTY0v4xMBRUpZKaW8QUqZZvt5xXb+vJTyNtv7M1LKqbbXJCnlM3p6hULhHD+jgUVp7TlT8ktqqKxv0Qwoeg9yNx3f5BGtnt4d7a6KXS5r9Wzraf9yoMCeJ354ZBCTEiM06/s6PhFQFAqFd7gmPoykaGv+eYuUbDt+yZ5qtjvuvvtul225o/WmbU9piysaqKlvXymfNW6YdZrwAEYFFIViCCOEIGtcvEPOlIaWSKf13Zlq6u40VW/Z9oTWbJF82eFB/MTECEZEBrtsx1dQAUWhGOIMCw9k+uj2/b+GG6+3P6DvTGVlZbfnr/LCzS84LXNHq6d3R/ujsT9yWatn25n24LkqrjS0MjNxJgF+BhZ2SNc8kFEBRaFQMCc1hvAg67K0FrN11ld36OX30MIdrTdt97W2pqnNPgECYG5qLKGBvrLG3D1UQFEoFAT6GVmUPgyAwrpC8kuquVTb3KWe3rqKn376U6dl7mj19O5oXz79sstaPdudtVJKth+/RJvZ+qDqaOUepo3qujv0QEUFFIVCAUBafBhjYq2L/KSEL45fQnZ6Qp+Rob0zrxbuaL1puy+1py7VU1TRAIAQMGpYDUbDwH4Q3xEVUBQKBWB9QH/9uHisW+XBhZpmDpXWONQJCwtzuX13tN603Vfa5jYz2SfaH8RPSYokNLjV5bZ9ERVQFAqFnejQAGYOa78Fs+tUBbXNbfbjnJwcTf2y8c43+nZHq6d3R7swbqHLWj3bHbU5pyupb7GuiA8NNDJ/bJxPJslyBxVQFAqFA7eNSicm1LqXVKvJwrZj7be+7r33Xk3to3MfdVrmjlZP7472nqR7XNbq2b6qLa9pIr+02n4+a1w8Qf5GFVAUCsXg5oOyv3LTxASurrErqmjgWLk1g+CWLVs0tUs3Ok8o5Y5WT++O9pcFv3RZq2d7y5YtmMwWPjt60b5gNCUulLR4662w5796XrPtgYYKKAqFwoFGcyOJUcFM7TD76MvCyzS0mGhp6ZqDviN1ttS13eGOVk/vjrbR3OiyVs92S0sLe85cobLe+qwkwM/A9ePj7Svi9fo90FABRaFQdMuCsXH2zSOvPlBevny5y+25o3VX7y3tjbffRe7Z9q0FF14TR6TtdzoYUQFFoVA4EBNg3T49wM/AjRNsyU+rCyj88ime/39TYVMyFG3oVjsyfKTTdtevX69pV0urp3dHGxeovUrdpX4XbcD0t1SeffwO5LEXobqAUTEhTEly3NZmRNgIzbYHGiqgKBQKB+5IvMP+fkxsKJP8DkLZZmirwTRsEnV1l2Dv6m6Dyoa7uw80ADNnztS0q6XV07ujfWLCEy5ru7VdtAH2ribnUiRxiSOgrYaA83/jptAdXTZ//JfMf9Fse6AxONb7KxSKPuP7+77PmOIx9mPL5f3UNd2AxRJKbYM/+4rnEhqwHVHyIMStc9BeO/panlr8lEt2f73t117RvlH0hlvpcMf81xhSolPaT1TswdQWQX3rBGoay9hXC/+/vbsPrqq8Ezj+/d2bm4RAQiAhISCYWAFfkGB5sdi3UCsqfWHr2BdfqnbaZZzWnbbbndWpO8WZ3e7anZ3W6bY7XWpd3bGItqKCpUWLRp1FKKJAREqMQAIJAQIhIQlJbu797R/nJiQh9+SSe3Jzb+7vM5NJznme37nPL8mc5z73nPM8OZk7+N1Lv4ftA5+c/8GyH/CFee43FKQSG6EYYwbwy8AFnnzaRk7gLUQUX1aAYKiE7tBcCF94sfpI65Gox921a5fr675Z96ZruVt8PLFVLVUjjgXoCg38PYRDSnvwelQFX1aAgL+RTH/NkL+v3Cz3pYtTTVKMUERkKvAMUAocBr6iqs2D6syL1Ol1GfAjVX1URB4G/hbofQz1h6q6GWPMRXts8WMD37G/UAodtbx5dj9v514FQIZcy+2XBCj8cuWA2IcrH4563HvvvTeudsUTP5qx9y2+j4crHgacubo2P34r1S3OFDbZC/O5q2AzuX4g51L4m8oRtyMVJMsI5UFgq6rOAbZGtgdQ1QOqulBVFwKLgA7g+X5VftZbbp2JMR4q/zH4c1g2aS/TAs77vB7J5o8T1tATCsd8mHXr1sXVjHjiExW7r6GV6ry7wOfcyXVj3g5y/efAn+P8Hse5ZOlQVgFPRn5+EnCfRwFuAD5U1dpRbZUxBsruhKVryZg4i1smbyMjMxdmfoGmrPIBc1OBc00gmqysLNeXWXer+4nbLT6e2H+6wv2i/HDt7s25ub3bWTQr/xqY+QUWTGnh8ux6Z2SydK3zexznZPBsomPSCJEzqprfb7tZVae41H8ceEdVfxHZfhi4F2gF3gZ+MPgjs36xq4HVAMXFxYvWr18/oja3tbXFPdldqrGc08NwOR9sCbHnxPkFuBZPz2BWrvPetK6jjtk5s4eMq6urY/bsocsAdpzawXUF10Utd4uPJ7ayvpKKmRUjigUn5xnZs3j9aA+tXc75NDdTqJiVQUaSziQc7//18uXLd6nq4gsKVDUhX8CfgfeG+FoFnBlUt9nlOJlAE1Dcb18x4McZcf0YeDyWNi1atEhH6rXXXhtxbKqynNPDcDmHw2H9w94G/enLB/SnLx/Q/9xarU1nO1VVdc1ra6LGrVkTvUxV9dP/82nXcrf4eGLLf1Y+4lhV1R+9ukY39/t9/PzP1Xq89ZxrzFiL9/8aeFuHOKcm7KK8qn42WpmIHBeRElU9JiIlwAmXQ92CMzo53u/YfT+LyK8B98l3jDEjJiLccGURJ892cbq9m2BI+UPVMb62JPq7eIBly5bF9brxxI9mbFPLRP7ac34KleVXFFGUmz3i10tlyXINZSNwT+Tne4AXXereDjzdf0ekE+r1JZyRjzFmlGRl+PncghICfucjnVNt3WzZ19g3AeJQ2tra4nrNeOJHK7bhzDkamvL6tufPnMz8mZOj1h/vkqVDeQS4UUQ+AG6MbCMiM0Sk744tEcmJlG8YFP/vIlIlInuB5cD3E9NsY9JX4aQsPnNFcd92zYk2ynJWRK1fVeX+vMe3rv2Wa7lbfDyxK0tWjii2tTPIS3sbmDN1HgDFedksnzfN9VjjXVI8h6Kqp3Du3Bq8vwFY2W+7AygYot7XR7WBxpghXTUjj5NtXbxT69wD09xayF8bW7liet4FdVevXu16rLvK73Itd4uPJ3ZFcfROMFpsV0+IF3c30N4VoiR3BtkBZ8SW4U+W9+hjI72zN8bE7ZOXF1Ja6DzIV3m4klf2HaexpfOCemvXrr1gX38VT1S4lrvFxxP7vd3ui3MNjg2HlT+910jTWefJ9zdqX+PzC0rG9SzCsbIOxRgTF59PuGV+Sd8qjz1h5cXd9TS3D1wvvaDggg8XLko88V7Fqiqvf3CSgyfb+/ZdUtTCrKk5Iz7+eGIdijEmbtkBP6sWziDD7zw539EdYsO79Zzttx59PBMwxhvvVeyOQ6fZXXd+Kd8lpVMpyHNfoCudWIdijPFEfk4mN1w1ue/Or9ZzQV54t57OoPMQ5HPPPecaf3XR1a7lbvHxxJbmlMYU+25dM299eKpv/9ziXD5+eQFzC+a6xqcT61CMMZ759sdu53MLZuCLrPvR1NbN85FOZbhRwi9X/tK13C0+ntjvzXW/hlJRUcH7Da0Dppm5tCCHm64uRkS445r4VqIcT6xDMcZ4Zl3VOsoKJ7Li6vO3EzfW7mTD43dz+NkbXVd7/M7m77geu6GhIWpZPLGPVj86dMGh38ILpexY91Ve3vAgnHFuH56Rn83nF8zou6NrXVV8k16OJ9ahGGM8U32qGoArS/K44coi5yRcv4nj7bDxUCnn2hqjrva478Q+92NXV0ctiyf2cMfhC3dGVl3c05TJq7VFaHcL1G+isGsPqxbOJDPj/KmzN2eTJM+hGGPGhyd2P0Hl4cq+7a7GE5zrWoiqELxqPh8/OI+JmZX4h1jtcTjDPccy0tiatpoLbjvWk9vp6r6cc8FygvM7eKIVMnzHmbj/m/ymaeDHaxWlA2PTmY1QjDGeyc/OH7Cd5dvPhMAORJSuU2cJhSfT1rWCnp6JF8SW5pe6Hnu451hGGpvpyxywrSqc6yrnXLAcgK5TZ8nwNTExays+PXtBvHUo59kIxRjjmd337R6444VS6DjIgc4QW/KWEVKACQSybuOmW/6DOcWxL4FbUlIStazy3soRx264fkPfRfvOYIjNVceobftXCLYAcMli+GL+q2T5etJi1cV42AjFGDN6Iqs9zsuu5dYpW8nydYMvQLBoBS/tPcbr1ScJhZ0ZJXc1uK/dvnjxhctv9Hpqz1Mjjj1w9gAAjS2d/HZHHbWnOmD6DeALcOWEQ9w65TWnM0mTVRfjYR2KMWb0RFZ7JOdSLsls4qsz9zP5ss85qxoC79Q28/tdRzjT0c2m6k2uh9q0KXr5Y+8+NuLYbU1v8U5dM8++fYTWc5EHMfOvYdmSm7lpegN+0bRadTEe9pGXMWZ0ld3ZdyIuAO4Ihtiyr7Fv+pKGM508tb2WE60TUVVEhl7lcMUK90kc3USLbW7vprGllNf7PWOSFfCx4qrpXF40F/jaiF8zHdkIxRiTUNkBP18sn8En5hT2PQAZDCn1JyfzzM4jNJw5N2Sc262/wxkc29UTYltNE09tr6UzeH4eruK8bO5ceimXF6XXss9esRGKMSbhRIQlpVOZPTWHl/c10tTWzTVF8znW0skzO48wp3gSS8umDlj58PDhw1GP98DHH3B9vd7Y7p4w+xpa2Hn4NO1dzpQwZRPL8ImwpHQKS8umpv0U9PGwDsUYM2aK87K5fels/nL4NP9Xc3524g+Ot/HB8TYuLchhAW9SWvsQq6cdgheecC6MD7qWUV5cPvQLHPot7HmIr0w5xbYnt7EnbzWdkxYMqDIjJ4fbr5uVtsv2eikpumIR+bKI7BORsIhEvR1DRG4WkQMiUiMiD/bbP1VEXhGRDyLfpySm5caYeGX4fVz/kUJO+55n3vSBtxHXfriDTW9s4de1H+UfN89k/ynh7FvfRQ8OfNL+jg0D59PqCYWp37uO7Vt/xtNHruD7m0vZcWoanbUv9U2hkpudwU1XT6def2ediUeSZYTyHnAr8N/RKoiIH/glzhLAR4GdIrJRVd8HHgS2quojkY7mQcB9DGyMSSpP7/sNb9VvpSc4ka5zMwl2F6LtR0FDQBanA7OoOnIZAL6aTfjytuPzdYOEOXTGx5Kfr0Y1g3Aom3A4OxLrLM97ekILja0AQfxn15FVGCAz+wRPfRimVErHJuFxKCk6FFXdD0S9uyNiKVCjqgcjddcDq4D3I98rIvWeBCqxDsWYlNL7lH1GoJ2MQDWhUC3dXY10h0oJhycxdeFlfXXDYT/h4PnRTK5k0d01aNlhDfX9WHBtKRm+o2RmHCLgO4pM+FRf2cL8haOUUfoRVR3rNvQRkUrgH1T17SHKbgNuVtVvRba/DlynqveLyBlVze9Xt1lVh/zYS0RWA6sBiouLF61fv35EbW1ra2PSpPS6E8RyTg9JlfOZKjTUTXMwg5OdAU52BWgOBujRDMguihomAhODxyjM7KQwK0hRVjdZ/si5zpfZ9xwMJFm+CRJvzsuXL9+lqhdcnkjYCEVE/gxMH6LoIVV9MZZDDLHvontDVV0LrAVYvHixjnQlt8rKyrhXoEs1lnN6SKqcD9U7sxOHOyATyAT15dB27VrOTPsEncEQXT1hwqpkZvgI+H3kZQfIzwkQqHsa/vL3EOqAIM6XPyfygGJF30skVb4JMlo5J6xDUdXPxnmIo8CsftuXAL2LHBwXkRJVPSYiJcCJOF/LGJMMeu/m2vMQdNRBzmyk/Mfklt3JsLOADRE71B1ixjtJcQ0lRjuBOSJSBtTjPMLae2vHRuAe4JHI91hGPMaYVNDvSfuExpqLliy3DX9JRI4Cy4A/iMiWyP4ZIrIZQFV7gPuBLcB+4FlV7V1V5xHgRhH5AOcusEcSnYMxxqS7pBihqOrzwPND7G8AVvbb3gxsHqLeKeCG0WyjMcYYd0kxQjHGGJP6rEMxxhjjCetQjDHGeMI6FGOMMZ5IqiflE01ETgK1IwwvBJo8bE4qsJzTQ7rlnG75Qvw5X6qq0wbvTOsOJR4i8vZQUw+MZ5Zzeki3nNMtXxi9nO0jL2OMMZ6wDsUYY4wnrEMZubVj3YAxYDmnh3TLOd3yhVHK2a6hGGOM8YSNUIwxxnjCOhRjjDGesA5lGCJys4gcEJGayHr1g8tFRH4eKd8rIh8di3Z6JYZ874zkuVdEtolI+Vi000vD5dyv3hIRCUVWD01pseQsIhUisltE9onI64luo9di+N+eLCKbRGRPJOdvjEU7vSIij4vICRF5L0q59+cuVbWvKF+AH/gQuAxnvbg9wFWD6qwE/oizouTHgB1j3e5Rzvd6YErk51tSOd9Yc+5X71Wc2a5vG+t2J+DvnA+8D8yObBeNdbsTkPMPgZ9Efp4GnAYyx7rtceT8KeCjwHtRyj0/d9kIxd1SoEZVD6pqN7AeWDWozirgf9WxHciPrBqZiobNV1W3qWpzZHM7zsqZqSyWvzHA3wHPMT5WA40l5zuADapaB6CqqZ53LDkrkCsiAkzC6VB6EttM76jqGzg5ROP5ucs6FHczgSP9to9G9l1snVRxsbl8E+cdTiobNmcRmQl8CfhVAts1mmL5O88FpohIpYjsEpG7E9a60RFLzr8ArsRZWrwK+K6qhhPTvDHh+bkrKRbYSmIyxL7B91nHUidVxJyLiCzH6VA+MaotGn2x5Pwo8ICqhpw3rykvlpwzgEU4C9dNAN4Ske2qWj3ajRslseR8E7Ab+AzwEeAVEXlTVVtHu3FjxPNzl3Uo7o4Cs/ptX4Lz7uVi66SKmHIRkQXAY8At6qyWmcpiyXkxsD7SmRQCK0WkR1VfSEwTPRfr/3WTqrYD7SLyBlAOpGqHEkvO3wAeUecCQ42IHAKuAP6SmCYmnOfnLvvIy91OYI6IlIlIJvA1YOOgOhuBuyN3THwMaFHVY4luqEeGzVdEZgMbgK+n8LvV/obNWVXLVLVUVUuB3wPfTuHOBGL7v34R+KSIZIhIDnAdsD/B7fRSLDnXEVlKXESKgXnAwYS2MrE8P3fZCMWFqvaIyP3AFpy7RB5X1X0icl+k/Fc4d/2sBGqADpx3OSkpxnx/BBQA/xV5x96jKTxTa4w5jyux5Kyq+0XkT8BeIAw8pqpD3n6aCmL8O/8z8ISIVOF8HPSAqqbstPYi8jRQARSKyFFgDRCA0Tt32dQrxhhjPGEfeRljjPGEdSjGGGM8YR2KMcYYT1iHYowxxhPWoRhjjPGEdSjGGGM8YR2KMcYYT1iHYkwSiay5sldEskVkYmRdjvlj3S5jYmEPNhqTZETkX4BsnEkZj6rqv41xk4yJiXUoxiSZyFxTO4FO4HpVDY1xk4yJiX3kZUzymYqzwFMuzkjFmJRgIxRjkoyIbMRZUbAMKFHV+8e4ScbExGYbNiaJRFZG7FHVdSLiB7aJyGdU9dWxbpsxw7ERijHGGE/YNRRjjDGesA7FGGOMJ6xDMcYY4wnrUIwxxnjCOhRjjDGesA7FGGOMJ6xDMcYY44n/Bzbi67APWaNXAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "t = np.arange(0, Lx, 0.001)\n",
    "\n",
    "plt.plot(t, t0(t), label = \"$T_0$\", lw = 3, alpha = 0.5)\n",
    "plt.scatter(xs, t0(xs), label = \"sampled\", color = \"orange\")\n",
    "plt.scatter(xs, np.zeros((xs.size)), color =\"green\", label = \"x nodes\", alpha = 0.5)\n",
    "plt.hlines(0,xmin=0, xmax = Lx, color = \"black\", ls = \"--\", lw = 0.5)\n",
    "\n",
    "for x in xs:\n",
    "    plt.plot([x, x], [0, t0(x)], color = \"black\", ls = \"--\", lw = 0.5)\n",
    "    plt.xlabel(\"x\")\n",
    "    plt.hlines(t0(x), xmin = x-deltax/2, xmax = x+deltax/2, color = \"green\")\n",
    "    plt.vlines(x-deltax/2, ymin=0, ymax=t0(x), color = \"green\", alpha = 0.5, ls=\"--\", lw = 1)\n",
    "    plt.vlines(x+deltax/2, ymin=0, ymax=t0(x), color = \"green\", alpha = 0.5, ls=\"--\", lw = 1)\n",
    "plt.ylabel(\"$T_0$\")\n",
    "plt.grid()\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Advection\n",
    "\n",
    "Now let's think of this as the temperature in a fluid, and let's suppose that the fluid is moving to the right with a constant, uniform speed $U$. (**Uniform** here means that the speed is the same in all parts of the fluid.) Then the temperature profile should also *move with the fluid*. We call a quantity, such as the temperature, that is carried along with the fluid a **tracer**.\n",
    "\n",
    "\n",
    "If we fix our attention at a single point in space, say the grid point $x_i$, the temperature there will vary over time due to the motion of the fluid. How it varies in time depends on the values at neighbouring grid points, since they determine how much heat will be *transported into* the current cell."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's think about a time-stepping algorithm, and let $\\delta t$ be the duration of a time step. We'll denote by $T^n_i$ the approximate value of $T$ to calculate at position $x_i$ and at the $n$th time step $t_n$, i.e. an approximation of $T(t_n, x_i)$, where $t_n = n \\, \\delta t$. \n",
    "For example, $T^{n+1}_i \\simeq T(t_n + \\delta t, x_i)$ and $T^{n}_{i+1} \\simeq T(t_n, x_i + \\delta x)$.\n",
    "\n",
    "\n",
    "**Note** that the superscript $n$ in these algorithms does not denote a power; it's just a label for the time step. We could write $T_i^{(n)}$ instead, but that is annoying to both write and read, so we suppress the parentheses."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose the fluid is moving to the right with speed $U$. During a time step of duration $\\delta t$, the temperature $T^n_i$ at cell $i$ changes due to the amount of heat that leaves cell $i$, and the amount that enters that cell.\n",
    "Most of the fluid that starts within cell $i$ remains within that cell during the time step (if the time step is short enough). Only fluid that is within a distance $U \\, \\delta t$ from the boundary of the cell will cross over that cell. So a *proportion* $(U \\, \\delta t) / \\delta x$ of the amount in cell $i$ crosses the boundary.\n",
    "\n",
    "\n",
    "Hence, roughly an amount $T^n_i (U \\delta t) / \\delta x$ will leave cell number $i$ and cross into cell $i+1$ (the cell to the right). Similarly, an amount $T^n_{i-1} (U \\delta t) / \\delta x$ will *enter* cell $i$ from the neighbouring cell $i-1$ on the left. Hence we have\n",
    "$$T^{n+1}_i = T^{n}_i + (T^n_{i-1} - T^n_{i})\\, U \\, \\delta t / \\delta x.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus\n",
    "$$\\frac{T^{n+1}_i - T^{n}_i}{\\delta t} =  \\frac{T^n_{i-1} - T^n_{i}}{\\delta x}\\, U.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Taking the continuum limit when $\\delta t \\to 0$ and $\\delta x \\to 0$, we recognise the definition of **partial derivatives** with respect to time and space variables from multivariable calculus. (Note the different indices that change on the two sides of the equation.) \n",
    "Denoting these partial derivatives using $\\partial$, we arrive at the **advection equation**:\n",
    "$$\\frac{\\partial T(t, x)}{\\partial t} = -U \\frac{\\partial T(t, x)}{\\partial x},$$\n",
    "or for short\n",
    "$$\\frac{\\partial T}{\\partial t} = -U \\frac{\\partial T}{\\partial x}.$$\n",
    "Since $T$ is a function of both $x$ and $t$, and this equation involves partial derivatives with respect to both of the independent variables, this is a **partial differential equation** (PDE). It describes how the function $T(t, x)$ changes continuously as a function both of time and space. \n",
    "Although there are some analytical methods to solve PDEs, often it's necessary to use numerical methods. Here we'll look at simple numerical methods to solve such equations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerics for the advection equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's return to the version of the equation in which the value at the *following* time step is isolated:\n",
    "$$T^{n+1}_i = T^{n}_i - \\left( U \\frac{\\delta t}{\\delta x} \\right) (T^n_{i} - T^n_{i-1}).$$\n",
    "In the last term on the right-hand side, we see that we require combinations of values of $T$ at the *same* time step from *different* places, with certain coefficients.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are many approaches to implementing this numerically. The simplest is to directly transcribe the equation for the $i$th entry of the vector.\n",
    "Calling `T` the current vector, i.e. $\\mathbf{T}^n := (T^n_i)_{i=1, \\ldots, N_x}$, and `T′` the new vector at the next time step, we have the following basic expression:\n",
    "\n",
    "> T′[i] = T[i] + δt * U * (T[i-1] - T[i]) / δx\n",
    "\n",
    "But now we realise a problem: What should we do when $i=0$? This will try to access the index -1 of the vector `T`, which does not exist!\n",
    "\n",
    "This illustrates the necessity of choosing **boundary conditions** that specify what happens at the edge of the domain.\n",
    "For simplicity we will choose to use **periodic boundary conditions**. This is a convenient mathematical fiction that allows us to treat all cells as being on the same footing, by wrapping the system around a torus, so that cells $i=0$ and $i=N_x$ are neighbours.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then write this as follows, where we separate out the case $i=0$:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def advection(T, deltat, deltax, U):\n",
    "    N = len(T)\n",
    "    T2 = T.copy()\n",
    "    \n",
    "    for i in range(1, N):\n",
    "        T2[i] = T[i] - deltat * U * (T[i]-T[i-1]) / deltax\n",
    "        \n",
    "    T2[0] = T[0] - deltat * U * (T[0]-T[-1]) / deltax\n",
    "    \n",
    "    return T2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use a function `plot_func` that carries out a generic time evolution. This is basically the Euler method again, except now we apply it to a vector of values at different points in space at each step. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "deltat = 0.001\n",
    "U = 0.2\n",
    "initcondition = t0(xs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.15643447,  0.4539905 ,  0.70710678,  0.89100652,  0.98768834,\n",
       "        0.98768834,  0.89100652,  0.70710678,  0.4539905 ,  0.15643447,\n",
       "       -0.15643447, -0.4539905 , -0.70710678, -0.89100652, -0.98768834,\n",
       "       -0.98768834, -0.89100652, -0.70710678, -0.4539905 , -0.15643447])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "initcondition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_func(function, deltat, U, xs, steps, T0 = initcondition):\n",
    "    T = T0.copy()\n",
    "    deltax = xs[1] - xs[0]\n",
    "    for step in range(steps):\n",
    "        T = function(T, deltat, deltax, U)    \n",
    "    plt.figure(figsize=(8,8))\n",
    "    plt.plot(xs, T)\n",
    "    plt.scatter(xs,T)\n",
    "    plt.title(\"{}\".format(step*deltat))\n",
    "    plt.xlabel(\"t\")\n",
    "    plt.ylabel(\"x(t)\")\n",
    "    plt.ylim([-1,1])\n",
    "    plt.grid()\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4398e9698cf14c4ba20ba4e7869c5561",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=1, description='Timesteps', max=10000, min=1), Output()), _dom_classes=(…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "interact(plot_func, \n",
    "         function = fixed(advection),\n",
    "         steps = widgets.IntSlider(description=\"Timesteps\", value = 1, min = 1, max = 10000, step = 1),\n",
    "         deltat = fixed(deltat), #widgets.FloatSlider(description = \"$\\Delta t$\",value = 0.001, min = 0.0001, max = 0.01, step = 0.0001, readout_format='.4f'),\n",
    "         U = fixed(U),         \n",
    "         xs = fixed(xs),\n",
    "         T0 = fixed(advection(initcondition, deltat, xs[0], U)));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unfortunately this does *not* behave as we expect: instead of preserving the shape of the profile over time, it is decaying. This is due to the way we are approximating. \n",
    "\n",
    "A better way to discretize the spatial derivative is using the following **centered difference**:\n",
    "$$\\frac{\\partial T(t_n, x_i)}{\\partial x} \\simeq \\frac{T^n_{i+1} - T^n_{i-1}}{2 \\delta x}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def advection2(T, deltat, deltax, U):\n",
    "    N = len(T)\n",
    "    T2 = T.copy()\n",
    "    \n",
    "    for i in range(1, N-1):\n",
    "        T2[i] = T[i] - deltat * U * (T[i+1]-T[i-1]) / (2*deltax)\n",
    "        \n",
    "    T2[0] = T[0] - deltat * U * (T[1]-T[-1]) / (2*deltax)\n",
    "    T2[-1] = T[-1] - deltat * U * (T[0]-T[-2]) / (2*deltax)\n",
    "    \n",
    "    return T2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "fcdb59f62019478fb022e250220e9237",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=1, description='Timesteps', max=10000, min=1), Output()), _dom_classes=(…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "interact(plot_func, \n",
    "         function = fixed(advection2),\n",
    "         steps = widgets.IntSlider(description=\"Timesteps\", value = 1, min = 1, max = 10000, step = 1),\n",
    "         deltat = fixed(deltat),#widgets.FloatSlider(description = \"$\\Delta t$\",value = 0.001, min = 0.0001, max = 0.01, step = 0.0001, readout_format='.4f'),\n",
    "         U = fixed(U),         \n",
    "         xs = fixed(xs),\n",
    "         T0 = fixed(advection(initcondition, deltat, xs[0], U)));\n",
    "         \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Diffusion\n",
    "\n",
    "Another key physical process is **diffusion**. This models how temperature or mass spreads out from hot or high concentration regions towards regions where it is cold or where there is a low concentration.\n",
    "The PDE describing this is the **heat equation**, or **diffusion equation**:\n",
    "\n",
    "$$\\frac{\\partial T}{\\partial t} = D \\frac{\\partial^2 T}{\\partial x^2}.$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To obtain a numerical method to solve this equation, we again need to discretise this, in particular the second derivative. One possible discretisation is\n",
    "$$\\frac{\\partial^2 T}{\\partial x^2}(t_n, x_i) \\simeq \\frac{T^n_{i+1} - 2 T^n_i + T^n_{i-1}}{\\delta x^2}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This may again be transcribed directly into code:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def diffusion(T, deltat, deltax, D):\n",
    "    N = len(T)\n",
    "    T2 = T.copy()\n",
    "    \n",
    "    for i in range(1, N-1):\n",
    "        T2[i] = T[i] + deltat * D * (T[i+1]-2*T[i] + T[i-1]) / (deltax**2)\n",
    "        \n",
    "    T2[0] = T[0] + deltat * D *     (T[1]-2*T[0] + T[-1]) / (deltax**2)\n",
    "    T2[N-1] = T[N-1] + deltat * D * (T[0]-2*T[N-1] + T[N-2]) / (deltax**2)\n",
    "    \n",
    "    return T2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "aebebd168b33464980e5643833f5130d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=1, description='Timesteps', max=10000, min=1), Output()), _dom_classes=(…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "interact(plot_func, \n",
    "         function = fixed(diffusion),\n",
    "         steps = widgets.IntSlider(description=\"Timesteps\", value = 1, min = 1, max = 10000, step = 1),\n",
    "         deltat = fixed(deltat),#widgets.FloatSlider(description = \"$\\Delta t$\",value = 0.001, min = 0.0001, max = 0.01, step = 0.0001, readout_format='.4f'),\n",
    "         U = fixed(0.01),         \n",
    "         xs = fixed(xs),\n",
    "         T0 = fixed(advection(initcondition, deltat, xs[0], U)));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Advection--diffusion\n",
    "\n",
    "Finally we can combine both mechanisms, to describe a tracer that is both being advected at a constant speed and diffusing. This basically utilises the composition of the advection and diffusion functions:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def advection_diffusion(T, deltat, deltax, U):\n",
    "    D = 0.01\n",
    "    temp = advection2(T, deltat, deltax, U)\n",
    "    return diffusion(temp, deltat, deltax, D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "59dfbc03caf245f29c6bfcd28102b113",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=1, description='Timesteps', max=10000, min=1), Output()), _dom_classes=(…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "interact(plot_func, \n",
    "         function = fixed(advection_diffusion),\n",
    "         steps = widgets.IntSlider(description=\"Timesteps\", value = 1, min = 1, max = 10000, step = 1),\n",
    "         deltat = fixed(deltat),#widgets.FloatSlider(description = \"$\\Delta t$\",value = 0.001, min = 0.0001, max = 0.01, step = 0.0001, readout_format='.4f'),\n",
    "         U = fixed(0.2),         \n",
    "         xs = fixed(xs),\n",
    "         T0 = fixed(advection(initcondition, deltat, xs[0], U)));"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}