03_2D_diff_equation_MOL.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SciPy solve_ivp example for solving 2D diffusion equation in a rectangular grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import solve_ivp\n",
    "import ipywidgets as widgets\n",
    "import matplotlib.ticker as mticker\n",
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modelled system\n",
    "\n",
    "The modelled system is a rectangular surface of a pond. When a contaminant is spilled in the pond it will diffuse on the surface. The evolution of the contaminant concentration C is described by Fick's second law\n",
    "\n",
    "\\begin{align*}\n",
    "    \\frac{\\partial C}{\\partial t} &= D \\Delta C \\\\\n",
    "    &= D \\left( \\frac{\\partial^2 C}{\\partial x^2} + \\frac{\\partial^2 C}{\\partial y^2} \\right)\n",
    "\\end{align*}\n",
    "\n",
    "where D is the diffusion coefficient of the contaminant in water\n",
    "\n",
    "We well be using the [method of lines](http://hplgit.github.io/prog4comp/doc/pub/p4c-sphinx-Python/._pylight006.html#reduction-of-a-pde-to-a-system-of-odes) approach to transform the PDE into a system of ODEs. This is done by discretizing the spatial dimension.\n",
    "\n",
    "We will follow the central difference scheme in the above link. The change of concentration in an element i is\n",
    "\n",
    "\\begin{equation*}\n",
    "    \\frac{\\mathrm{d}C_{i,j}}{\\mathrm{d}t} = D\\left( \\frac{C(x_{i+1},t) - 2C(x,t) + C(x_{i-1},t)}{(\\Delta x)^2} +\\frac{C(y_{i+1},t) - 2C(y,t) + C(y_{i-1},t)}{(\\Delta y)^2} \\right)\n",
    "\\end{equation*}\n",
    "\n",
    "Finally, we need to decide what is the change of concentration at the boundaries of the rectangle (boundary conditions). For simplicity we will assume Dirichlet boundary condition and assume that the concentration at the boundary is \n",
    "\n",
    "\\begin{equation*}\n",
    "    C_{\\mathrm{boundary}} = C_0\n",
    "\\end{equation*}\n",
    "\n",
    "Now that we have discretized the equation we are left with $(N-2) \\cdot (M-2)$ ordinary differential equations of the form\n",
    "\n",
    "\\begin{align*}\n",
    "\\frac{\\mathrm{d}C_{1,1}}{\\mathrm{d}t} &= D\\left(  \\frac{C_{2,1} - 2C_{1, 1} + C_{0, 1}}{(\\Delta x)^2} +\\frac{C_{1,2} - 2C_{1, 1} + C_{1,0}}{(\\Delta y)^2} \\right) \\\\\n",
    "&\\;\\;\\vdots \\\\\n",
    "\\frac{\\mathrm{d}C_{N-1,M-1}}{\\mathrm{d}t} &= D\\left(  \\frac{C_{N,M-1} - 2C_{N-1, M-1} + C_{N-2, M}}{(\\Delta x)^2} +\\frac{C_{N-1,M} - 2C_{N-1, M-1} + C_{N-1, M-2}}{(\\Delta y)^2} \\right)\n",
    "\\end{align*}\n",
    "\n",
    "where N is the number of horizontal elements in the rectangular grid and M is the number of vertical elements. Note that in the last equations we have omitted the time subscript for simplicity. We do not need to consider the ODEs at the boundary since their concentration does not change due to the choice of the boundary conditions (we'll set the change in concentration to zero in RHS ODE function). If we would have needed Neumann type boundary conditions, it would have resulted to system of $(N+2) \\cdot (M+2)$ ODEs\n",
    "\n",
    "Next, we are going to solve this ODE system in the domain\n",
    "\n",
    "\\begin{align*}\n",
    "t &\\in [0, 100] \\\\\n",
    "\\{x,y\\} &\\in [0, 10]\n",
    "\\end{align*}\n",
    "\n",
    "with initial conditions\n",
    "\n",
    "\\begin{align*}\n",
    "C_0 = 10^{-3} \\\\\n",
    "C_{4.5,4.5 - 5.5,5.5} = 1\n",
    "\\end{align*}\n",
    "\n",
    "and with the Dirichlet boundary condition specified above"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Diffusion coefficient\n",
    "D = 1e-1\n",
    "\n",
    "# Grid resolution and max length\n",
    "dx = 0.25\n",
    "dy = 0.25\n",
    "lx = 10\n",
    "ly = 10\n",
    "\n",
    "# Number of gridpoints\n",
    "N_x = int(lx/dx) + 1\n",
    "N_y = int(ly/dy) + 1\n",
    "\n",
    "# grid point vectors\n",
    "x = np.linspace(0,lx,N_x)\n",
    "y = np.linspace(0,ly,N_y)\n",
    "\n",
    "# Concentration vector of the contaminant\n",
    "C = np.zeros((N_x,N_y)) + 1e-3\n",
    "\n",
    "# Select region between(4.5,4.5)(5.5,5.5)\n",
    "ind_x = np.where((x >= 4.5) & (x <=5.5))[0]\n",
    "ind_y = np.where((y >= 4.5) & (y <=5.5))[0]\n",
    "\n",
    "for ix in ind_x:\n",
    "    for iy in ind_y:\n",
    "        C[ix,iy] = 1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(41, 41)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define RHS function of the ode system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def odefun(c, D, dx, dy):\n",
    "    dcdt = np.zeros(c.shape)\n",
    "    dcdt[1:-1, 1:-1] = D*((c[2:, 1:-1] - 2*c[1:-1, 1:-1] + c[:-2,1:-1] ) / (dx)**2 + c[1:-1, 2:] - 2*c[1:-1, 1:-1] + c[1:-1,:-2] ) / (dy)**2\n",
    "    return dcdt "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-5.684341886080802e-14\n"
     ]
    }
   ],
   "source": [
    "t = odefun(C,D,dx,dy)\n",
    "print(np.sum(t))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run solve_ivp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "`y0` must be 1-dimensional.",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m/var/folders/_7/7fgf24sd5ggbyr5x9kzrwz7c0000gn/T/ipykernel_4262/2380386514.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0msol\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msolve_ivp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;32mlambda\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mc\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0modefun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mc\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mD\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdx\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdy\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m100\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mC\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/scipy_ode_solver/.venv/lib/python3.9/site-packages/scipy/integrate/_ivp/ivp.py\u001b[0m in \u001b[0;36msolve_ivp\u001b[0;34m(fun, t_span, y0, method, t_eval, dense_output, events, vectorized, args, **options)\u001b[0m\n\u001b[1;32m    540\u001b[0m         \u001b[0mmethod\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mMETHODS\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mmethod\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    541\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 542\u001b[0;31m     \u001b[0msolver\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfun\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvectorized\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mvectorized\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0moptions\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    543\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    544\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mt_eval\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/scipy_ode_solver/.venv/lib/python3.9/site-packages/scipy/integrate/_ivp/rk.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, fun, t0, y0, t_bound, max_step, rtol, atol, vectorized, first_step, **extraneous)\u001b[0m\n\u001b[1;32m     87\u001b[0m                  first_step=None, **extraneous):\n\u001b[1;32m     88\u001b[0m         \u001b[0mwarn_extraneous\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mextraneous\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 89\u001b[0;31m         super().__init__(fun, t0, y0, t_bound, vectorized,\n\u001b[0m\u001b[1;32m     90\u001b[0m                          support_complex=True)\n\u001b[1;32m     91\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0my_old\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/scipy_ode_solver/.venv/lib/python3.9/site-packages/scipy/integrate/_ivp/base.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, fun, t0, y0, t_bound, vectorized, support_complex)\u001b[0m\n\u001b[1;32m    117\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mt_old\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    118\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mt0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 119\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_fun\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0my\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcheck_arguments\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfun\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msupport_complex\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    120\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mt_bound\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mt_bound\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    121\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvectorized\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mvectorized\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/scipy_ode_solver/.venv/lib/python3.9/site-packages/scipy/integrate/_ivp/base.py\u001b[0m in \u001b[0;36mcheck_arguments\u001b[0;34m(fun, y0, support_complex)\u001b[0m\n\u001b[1;32m     15\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     16\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0my0\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndim\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 17\u001b[0;31m         \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"`y0` must be 1-dimensional.\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     18\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     19\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mfun_wrapped\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: `y0` must be 1-dimensional."
     ]
    }
   ],
   "source": [
    "sol = solve_ivp(lambda t,c: odefun(c,D,dx,dy), (0, 100), C)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What happened?\n",
    "\n",
    "We ran into error ```ValueError: `y0` must be 1-dimensional.```. The error is due to the fact that our matrix C is 2-dimensional but for historical reasons the solvers operate with one dimensional arrays. The ode solvers behave the same way in many other languages besides Python. We need to reshape our initial condition matrix C into a 1-dimensional array and because the odefun must return a matrix that has the same shape as C we need to also modify the odefun"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def odefun(c, D, dx, dy, Nx, Ny):\n",
    "    conc = c.reshape(Nx,Ny)\n",
    "    dcdt = np.zeros(conc.shape)\n",
    "    dcdt[1:-1, 1:-1] = D*((conc[2:, 1:-1] - 2*conc[1:-1, 1:-1] + conc[:-2,1:-1] ) / (dx)**2 + (conc[1:-1, 2:] - 2*conc[1:-1, 1:-1] + conc[1:-1,:-2]) / (dy)**2)\n",
    "    return dcdt.reshape(Nx*Ny,)\n",
    "C = C.reshape(N_x*N_y,)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1681, 101)\n"
     ]
    }
   ],
   "source": [
    "sol = solve_ivp(lambda t,c: odefun(c,D,dx,dy, N_x, N_y), (0, 100), C, t_eval=np.linspace(0,100, 101))\n",
    "print(sol.y.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d3f4528ea1844dcfa6f35a4a1952dca4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=0, description='Time'), Output()), _dom_classes=('widget-interact',))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6709220b260f4c6e862d793f0c7b47d0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "\n",
       "            <div style=\"display: inline-block;\">\n",
       "                <div class=\"jupyter-widgets widget-label\" style=\"text-align: center;\">\n",
       "                    Figure\n",
       "                </div>\n",
       "                <img src='data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAA+gAAAH0CAYAAACuKActAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAA9hAAAPYQGoP6dpAABWwUlEQVR4nO3deViVdf7/8dcBBBUELBfQWERJKytTMVxyw6UaxyW/5ZaI5VJNm1pjOplLi5pLTrNkagktmpZlTU1quTEuqYU6Or8ylEXNpdLkiCnr+f3hcMYToNxxzrlv6Pm4rvu65L4/98370HU6vPh87vdtczgcDgEAAAAAAFP5mF0AAAAAAAAgoAMAAAAAYAkEdAAAAAAALICADgAAAACABRDQAQAAAACwAAI6AAAAAAAWQEAHAAAAAMACCOgAAAAAAFgAAR0AAAAAAAsgoAMAAAAAYAEEdAAAAAAALICADgAAAACABRDQAQAAAACwAAI6AAAAAAAWQEAHAAAAAMACCOgAAAAAAFgAAR0AAAAAAAsgoAMAAAAAYAEEdAAAAAAALICADgAAAACABRDQAQAAAACwAAI6AAAAAAAWQEAHAAAAAMACCOgAAAAAAFgAAR0AAAAAAAsgoAMAAAAAYAEEdAAAAAAALICADgAAAACABRDQAQAAAACwAAI6AAAAAAAWQEAHAAAAAMACCOgAAAAAAFgAAR0AAAAAAAsgoAMAAAAAYAEEdAAAAAAALICADgBAFWGz2QxvXbt2lSR17dpVNptNmzZtMvU1eNpXX32lu+++Ww0bNlTNmjXVpEkTPfLII/r+++9/9TVPnjyphx9+WE2aNFFAQIAaNmyou+++W2lpaZc9Lz8/X7Nnz9bNN9+swMBA1a1bV127dtV77733q2sBAFRvfmYXAAAAKmbEiBGl9p04cUJr164t93iLFi08XpdVvPfeexoyZIgKCwsVFxenJk2a6Msvv9Rf//pXvfvuu9qyZYuaNWtm6JrffvutbrvtNn3//feKiYlR//79lZmZqffee0+rV6/WypUrNWDAgFLn/fzzz+rZs6e2bdum0NBQ3X777crNzdWGDRu0efNmTZgwQXPnznXXSwcAVBM2h8PhMLsIAADw62zatEndunWTJF3uI/3w4cP6+eefFRkZqdq1a3urPK85duyYYmNj9fPPP+vVV1/VmDFjJElFRUVKSkrSW2+9pbi4OO3YsUM2m61C13Q4HGrTpo12796t4cOHa+nSpfL19ZUkLVq0SGPHjlVQUJDS09MVFhbmcu7jjz+uP//5z7rxxhu1YcMG1atXT9LFGf6uXbsqNzdX//jHP9SnTx83/hQAAFUdS9wBAPgNiIyMVIsWLaplOJekBQsW6Oeff1aPHj2c4VySfH199corrygkJES7du3SunXrKnzNTz/9VLt371ZoaKj+/ve/O8O5JI0ZM0YJCQnKzc3Vn//8Z5fzfvrpJ73yyiuSpFdeecUZziWpTZs2mjhxoiTp+eef/1WvFQBQfRHQAQD4DSjvHvSkpCTZbDYlJyfrwIEDGjRokBo0aKDAwEDFxcXpww8/dI7dsWOH+vbtq/r166tWrVpq37691q9fX+73PH/+vObNm6f4+HiFhoaqZs2aat68uf74xz/q1KlTbn19H3zwgSRp6NChpY4FBQWpb9++kqT333/f8DX79u2roKCgUsdLvtcvr/nPf/5T+fn5ioyMVMeOHcs974svvtCxY8cqXA8AoPojoAMAAKWlpalNmzbau3evEhISdPPNN+vLL7/UgAEDnPdb33bbbTp69KgSEhLUvHlzffHFF7r99tu1ZcuWUtc7duyYbr31Vj3xxBNKT09XXFyc7rzzTuXl5WnOnDlq27atsrOzS51X8oeEadOmVbj2s2fP6uDBg5Kktm3bljmmZP/u3bsrfN2SsVe6Znp6us6dO1fh82JiYnTVVVdJkvbs2VPhegAA1R8BHQAA6C9/+YsmTZqkr7/+WsuXL9e2bdv08ssvy+FwaNy4cbrvvvv02muvKS0tTe+884727Nmjxx9/XIWFhZo+fbrLtRwOh+655x7t27dP999/v7KysvTZZ5/p/fff18GDBzVhwgRlZWVp5MiRbqk9KyvL+e/IyMgyx0REREiSMjMzK3zdkrFXuqbD4XCp4UrnSdI111xjuB4AQPVHQAcAAGrXrp0mT57s0kDtwQcf1FVXXaWjR4+qR48eGj58uMs5Tz/9tCQpNTVVBQUFzv1r167V1q1b1apVKy1cuFB16tRxHvPz89OLL76oli1bauPGjdq/f7/LNSMjI9W8eXOX+7av5OzZs85/BwYGljmmZIm63W43fN0rXfOX173Seb+2HgBA9cdj1gAAgO64445S3c39/PzUpEkTnT59WnfeeWepc66++mpdddVVOn36tE6dOuXsZP7JJ59IkgYOHCg/v9K/avj4+Khz587av3+/tm3bppYtWzqPvfHGG+58WQAAVCnMoAMAgHKXY5fM9JZ3vGR2/MKFC859GRkZkqQpU6bIZrOVuf3973+XJP3www+Vrv3SGfpL7wW/VG5uriQpODjY8HWvdM1fXvdK5/3aegAA1R8z6AAAQD4+l/+b/ZWOX6q4uFiS1KlTJzVt2vSyY2+44YYKX7c8UVFRzn8fPnxYN954Y6kxR44ckSRFR0dX+LrR0dE6ffq0Dh8+XObxkmvabDaXGkq+R3nnSdLRo0cN1wMAqP4I6AAAwK1Kmqf169dPTzzxhMe/X3BwsJo1a6aDBw/qyy+/LDOgf/nll5Kk1q1bV/i6rVu3VlpamvPc8q4ZGxvrcj96yfco77yMjAydPn1aknTLLbdUuB4AQPXHEncAAOBWd9xxhyTp3XfflcPh8Mr3HDBggCRp2bJlpY7l5ubqH//4hyTprrvuMnzNjz76qMzl6iXf65fXvPPOO+Xv76/Dhw9r69at5Z4XHx+vRo0aVbgeAED1R0AHAABu1a9fP8XFxWnnzp0aOXJkmfeZ//TTT1q4cKEKCwtd9icmJqpFixb661//auh7Pv7446pdu7Y+//xzLV682Lm/qKhIDz30kM6cOaO4uDj16tXL5bzvvvtOLVq0UIsWLfTdd9+5HLvjjjt0yy236MyZM3rooYdUVFTkPLZo0SKtX79eQUFBeuyxx1zOq1u3rh588EFJ0kMPPaRTp045j6WlpWn27NmSpD/96U+GXiMAoPpjiTsAAHArHx8frV69Wr/73e+UkpKi9957TzfffLMiIyOVn5+vjIwM7du3T0VFRUpKSnLp9H748GEdOHBAP/74o6Hv2ahRIyUnJ2vIkCEaM2aMXnvtNUVHR2vXrl3KyMhQw4YNtWzZslKd6gsKCnTgwAHnvy9ls9m0fPly3XbbbXrjjTe0ZcsWxcXFKTMzUzt37pSfn5/eeOMNZ/f6S73wwgvauXOntm/frtjYWHXv3l3nzp3T+vXrVVBQoPHjx6tPnz6GXiMAoPpjBh0AALhdo0aN9MUXX2jhwoVq166dDhw4oPfee09btmyRJD3wwANau3atatas6bbveffdd2vHjh266667lJGRoQ8++EBFRUX6wx/+oL1796pZs2aGr9m8eXP9+9//1h/+8AcVFRXpgw8+UGZmpu666y7t2LHDuQz+l2rXrq1NmzZp5syZaty4sf75z39q+/btat++vVauXKl58+ZV9uUCAKohm8NbN4cBAAAAAIByMYMOAAAAAIAFENABoBxvvfWWxo4dq7Zt2yogIEA2m03Jycnljrfb7Ro/fryioqIUEBCg6OhoPfnkk8rNzfVe0QAAAKiyWOIOAOWIjo5Wdna26tWrp8DAQGVnZ2vp0qVKSkoqNfbcuXPq1KmT9uzZo169eumWW27R7t27tW7dOsXFxSk1NdWt99oCAACg+mEGHQDKsWTJEmVlZemHH37QAw88cNmxL774ovbs2aOJEydq7dq1mjVrltauXauJEydq165deumll7xUNQAAAKoqAjoAlKNHjx6Kioq64jiHw6ElS5YoKChIU6ZMcTk2ZcoUBQUFacmSJZ4qEwAAANUEAR0AKik9PV3Hjh1Tx44dFRgY6HIsMDBQHTt2VEZGho4cOWJShQAAAKgK/MwuAACquvT0dElSbGxsmcdjY2O1du1apaenKyIioswxeXl5ysvLc35dXFys06dP6+qrr5bNZnN/0QAAy3I4HDp79qwaNWokHx/rzadduHBB+fn5brmWv78/PVqASxDQAaCScnJyJEkhISFlHg8ODnYZV5aZM2dq+vTp7i8OAFBlHTlyRNdcc43ZZbi4cOGC6teqJXc9nyQsLEyZmZmEdOC/COgAYAGTJk3S+PHjnV/n5OQoMjJS0jhJAabVBQAwQ56kl1SnTh2zCyklPz9fuXLPp1OepJdOnFB+fj4BHfgvAjoAVFLJzHl5M+R2u91lXFkCAgIUEFDWrzoBIqADwG+TlW9xCpRU2UhNEAFK430BAJVUcu95yb3ov3Sle9QBAKhqavx3q4widxQCVDPW6zoBAFVMbGysGjVqpK1bt+rcuXMux86dO6etW7eqSZMm5TaIAwAAACQCOgBUms1m06hRo5Sbm6tnn33W5dizzz6r3NxcjR492qTqAABwPz83bQBc8b4AgHIsWbJEW7ZskSTt27fPuW/Tpk2SpE6dOmnUqFGSpD/+8Y/68MMPNXv2bO3evVutW7dWWlqa1q1bp7i4OD3++ONmvAQAADzCT5Vf4l7ojkKAaoaADgDl2LJli1JSUlz2bd26VVu3bnV+XRLQAwMDtXnzZk2bNk2rVq3Sxo0bFR4ergkTJmjq1KmqVauWV2sHAABA1WNzOBwOs4sAALiy2+3/7fr+lOjiDgC/NXmSZiknJ0fBwcFmF+Oi5PNpnqTK/un5vKQJkiVfJ2AWZtABAAAAGOKOLu4scQdKo0kcAAAAAAAWwAw6AAAAAEPc0YWdIAKUxvsCAAAAgCHu6OJe4I5CgGqGJe4AAAAAAFgAM+gAAAAADGGJO+AZvC8AAAAAGOKOLu6VPR+ojljiDgAAAACABTCDDgAAAMAQZtABzyCgAwAAADCEe9ABz2CJOwAAAAAAFsAfrgAAAAAY4o7noBNEgNJ4XwAAAAAwhCXugGewxB0AAAAAAAvgD1cAAAAADKGLO+AZBHQAAAAAhrDEHfAMlrgDAAAAAGAB/OEKAAAAgCF0cQc8g/cFAAAAAENY4g54BkvcAQAAAACwAP5wBQAAAMAQurgDnkFABwAAAGAIS9wBz2CJOwAAAAAAFsAfrgAAAAAYQhd3wDN4XwAAAAAwhHvQAc9giTsAAAAAABbADDoAAAAAQ2gSB3gG7wsAAAAAhvj5SjVslbyGQ1KRW8oBqg2WuAMAAAAAYAHMoAMAAAAwxM9P8mMGHXA7AjoAAAAAQ2q4YYl7DYd7agGqE5a4AwAAAABgAcygAwAAADDEbUvcAbggoAMAAAAwpIavVKOSa3FrFLunFqA6YYk7AAAAAAAWwAw6AAAAAGN8VfmpvkoukQeqIwI6AAAAAGP8VPmAzhJ3oBSWuAMAAAAAYAHMoAMAAAAwhhl0wCOYQQcAAABgjJ+bNoPeeustjR07Vm3btlVAQIBsNpuSk5Mr+2oAy2AGHQAAAECV8PTTTys7O1v16tVTeHi4srOzzS4JcCtm0AEAAAAY46OLndwrs/2KJLJkyRJlZWXphx9+0AMPPFDplwFYDTPoAAAAAIzx08WQXRm/4jFrPXr0qOQ3BayNGXQAAAAAACyAGXQAAAAAxpg0gw5UdwR0AAAAAMaU3EfuBna73eXrgIAABQQEuOfiQBXDEncAAAAApomIiFBISIhzmzlzptklAaZhBh0AAACAMW5c4n7kyBEFBwc7dzN7jt8yAjoAAAAAY3zltiQRHBzsEtCB3zKWuAMAAAAAYAHMoAMAAAAwxh1N4hzuKASoXgjoAAAAAIzxkylJYsmSJdqyZYskad++fc59mzZtkiR16tRJo0aN8n5hgJsQ0AEAAABUCVu2bFFKSorLvq1bt2rr1q3Or6tLQD906JCSk5O1detWnThxQpIUFhamDh06aOTIkWratKnJFcITbA6Hg8UlAGAxdrtdISEhkp6SRDdbAPhtyZM0Szk5OZZrnlby+ZTTRQqu5FSfvVAK2SxLvk4zORwOPf3003rxxRdVVFQkm81WaoyPj4+efPJJPf/882UeR9XFDDoAAAAAY0xa4v5b8OSTT2r+/PmqX7++xowZo4SEBDVu3FiSdPToUW3cuFELFy7UrFmzdOHCBc2fP9/kiuFOzKADgAUxgw4Av2VVYAY9wU0z6OuZQb9UWlqa4uLi1Lp1a61Zs0ZXX311meN+/PFH3X777UpLS9OXX36p1q1be7lSeAqPWQMAAABgjI/+18n9124kkVJeeeUV+fn56d133y03nEtSvXr19N5776lGjRr629/+5sUK4WksTAEAAABgjDuWuLOOt5TU1FR1795d0dHRVxwbHR2t7t2761//+pfnC4PX8HcrAAAAALCA7777Ti1btqzw+JYtW+rYsWMerAjexgw6AAAAAGOYQfeI4uJi1ahRo8Lja9SooeLiYg9WBG8joAMAAAAwpuQ+8sogV5YSHh6uzMzMCo/PyMhQo0aNPFgRvI0l7gDgJg6HQ++//766deum8PBw1a5dW82bN9fYsWOVkZFhdnkAAMDiOnbsqM8//1x5eXlXHJuXl6fPP/9cHTp08EJl8BYCOgC4yRNPPKGBAwfqwIED6t+/vx555BE1adJEixcvVqtWrbR//36zSwQAwD383LTBxciRI3X69OkKdWb/61//qp9++kmJiYleqAzewtsCANzgxIkTWrBggaKiorR3797/PsP8opdeeknjx4/X/Pnz9frrr5tYJQAAbuKryicJlriX0q1bNz3//PP68ccflZ+fL39//zLHFRQU6NSpU5oxY4Z69Ojh5SrhSQR0AHCDrKwsFRcXq2PHji7hXJL69Omj8ePH64cffjCpOgAAUFU89dRTVxxTo0YNvfDCC16oBt5GQAcAN4iNjZW/v7+2bt0qu92u4OBg57GPP/5YkpSQkGBWeQAAuJc7msRV9vxqrqCgQN98843z94oWLVoY6vCOqomADgBucPXVV2vWrFmaMGGCWrRooX79+ik4OFh79+7Vhg0b9NBDD+nhhx82u0wAANzDHfeQs8S9TD/++KMmTpyo5cuX68KFC879NWvW1ODBg/Xiiy+qXr16JlYITyKgA4CbjBs3To0bN9aoUaO0cOFC5/5OnTpp6NCh8vMr/3+5eXl5Lh1b7Xa7R2sFAADWc/r0acXHxysjI0P169dX9+7d1bBhQ508eVI7d+5USkqKUlNTtWPHDl199dVmlwsPoIs7ALjJjBkzdO+992ry5Mk6cuSIzp49q3/961+6cOGCunbtqo8++qjcc2fOnKmQkBDnFhER4cXKAQAwiC7uHvHss88qMzNTU6ZMUXZ2tj7++GO99tpr+vjjj3X48GFNnjxZmZmZevbZZ80uFR5iczgcDrOLAICq7vPPP1fPnj01btw4zZ8/3+XYiRMnFBMTo8aNGys9Pb3M88uaQb8Y0p+SFODBygEA1pMnaZZycnJceppYgd1uV0hIiHLGSMFlNxiv+LXypZBFsuTrNEtMTIxiYmL0+eeflzsmISFBmZmZysjI8GJl8BZm0AHADT799FNJFx+P8kthYWFq0aKFDh48qNzc3DLPDwgIUHBwsMsGAAB+W44dO6Zbb731smPi4+N17NgxL1UEb2NhCQC4QX5+viSV+yi1H374QT4+PnRfBQBUDz6qfBd2pgpLCQkJUWZm5mXHZGZmlnqkK6oP3hYA4AYdO3aUJM2fP185OTkuxxYuXKijR4+qffv2CghguToAoBrgHnSP6Natm9577z2tWbOmzONr1qzRe++9p65du3q3MHgN96ADgBsUFRWpe/fuSk1NVYMGDdS3b1+FhoYqLS1NGzZsUK1atbRp0ya1a9euQtcrucePe9AB4LeoCtyD/rAUXMmPJ3ueFPJX7kG/1P/7f/9P8fHxys3NVa9evdStWzeFhYXp+PHj2rhxoz777DPVrl1bO3bs0A033GB2ufAA/m4FAG7g6+urdevW6aWXXtLKlSu1bNky5efnq2HDhs7O7tddd53ZZQIA4B7umAEvckch1cv111+vtWvXauTIkfrss8/02WefyeFwyGazSZJiY2O1dOlSwnk1xgw6AFgQM+gA8FtWBWbQJ7hpBn0eM+jl2bp1q/bs2SO73a7g4GC1atXKeUsdqi9m0AEAAADAYjp27Egg/w0ioAMAAAAwhiXugEcQ0AEAAAAY46vKJ4lCdxRSvYwcObJC4xwOh5KTkz1bDExBQAcAAAAAC3jzzTd1uRZhJQ3jCOjVFwEdAAAAgDHuWOJOEiklIyOjzP12u11paWmaMWOGbrnlFs2dO9fLlcFbeFsAAAAAMMb3v1tlrwEXkZGR5R5r2bKlbr/9dt1444368MMP9eijj3qxMniLj9kFAAAAAACurEGDBurTp49eeeUVs0uBhzCDDgAAAMAYlribJjg4WJmZmWaXAQ/hbQEAAADAGAK6KU6dOqXVq1crLCzM7FLgIbwtAAAAAMACpk+fXub+wsJCHT16VB9++KHOnDmjadOmebcweA0BHQAAAIAxPqp8kze6YZUyY8aMyx6vU6eO/vSnP2nKlCleqgjeRkAHAAAAYAxL3D1i48aNZe739fVVaGiomjdvLj8/fnDVGf91AQAAAMACOnfubHYJMBkLSwAAAAAY4+emDYZ99NFHGjlypNllwEN4WwAAAAAwxleVvwe9sudXUz/99JO2b9+u06dPq6ioqNTxDz/8UKtXr1aXLl1ks9kkSSNGjPB2mfAQAjoAAAAAWMCaNWs0cOBAnT9/XpKcAfyXbDab7r//fkmSw+EgoFcjBHQAAAAAxtAkziOmTZumgIAATZw4UREREfL1Lb3MYPXq1Vq9erWSk5O9XyA8jrcFAAAAAGN8VfkkwRL3Uvbv369Ro0bpmWeeKXdMVlaWVq9ercTERC9WBm+hSRwAAAAAWMD58+dVv359s8uAiZhBBwAAAGAMS9w9orx7zn85piLjUDXxtgAAAABgDF3cPaKwsPCKY6ZMmaIpU6Z4oRqYgYAOAAAAABZRXFysnTt3Kjs7W02bNlXbtm3NLglexD3oAAAAAIzxc9MGF6dOnVKHDh3UoUMHDRkyRO3atVPfvn1VUFAgSTp69KgefPBBpaWlmVwpPIWADgAAAMAYArpHPPXUU9q1a5cGDRqkv/zlL7rjjjv08ccf6y9/+YskqXHjxvr444/1+uuvm1wpPIWADgAAAAAWsGbNGnXq1EnLly/XH/7wB3388cdq3ry5VqxYIelig7jOnTtr8+bNJlcKTyGgAwAAADDGR/9rFPdrN5JIKadPn1Z8fLzz65JA/s033zj3RURE6MiRI2aUBy9gYQkAAAAAY3jMmkfExsYqOzvbZV+DBg2Um5vr/Prnn39Wfn6+t0uDl/B3KwAAAACwgAcffFD/+Mc/9J///Me5r0aNGnI4HM6vN23apGbNmplRHryAgA4AAADAGJrEeUTv3r3Vp08fde7cWXPmzNGOHTt05swZSVJaWppGjRql/fv3a9iwYeYWCo/hbQEAAADAmJL7yCt7Dbho2rSpHA6HbDabJk6c6Nxvs9nUtm1b2Ww2DRw4UBMmTDCxSngSAR0AAAAALCApKcllOXuJgIAANWrUSD169FD79u1NqAzeQkAHAAAAYAxN4jzitddeM7sEmIy3BQAAAABjfFX5JMESd6AUmsQBAAAAAGABzKADAAAAMIYl7h7RpEmTCo1zOBzKysrybDEwBW8LAAAAAMbQxd0jfHx8ymwSZ7fb9dNPP0mSwsPD5e/v7+3S4CUEdAAAAACwgEOHDpV77LvvvtMTTzyhjIwMffbZZ16sCt7EPegAAAAAjPFz04YKa9y4sd5++23l5eW5PCMd1QsBHQAAAIAxJV3cK7OxxN0wHx8fde/eXR9++KHZpcBDCOgAAAAAUEWcPHlSdrvd7DLgISwsAQAAAGAMTeK8qqioSMePH9fy5cu1YsUKdezY0eyS4CEEdAAAAADG8Jg1j/D19S2zi/ulrrrqKi1YsMA7BcHreFsAAAAAgAV07dq13Mespaeny9fXV1u3blXz5s1NqA7eQEAHAAAAYAwz6B6xfv36co9duHBBU6ZMUe/evbV9+3aFh4d7sTJ4C03iAAAAABjDY9a8rmbNmpozZ44aNWqkJ5980uxy4CEEdAAAAACoIjp06KDPPvvM7DLgIfzdCgAAAIAhDh/JUcku7A6mCn+VgwcP6uzZs2aXAQ8hoAMAAAAwpMjv4lbZa8BVdnZ2mfsLCwudj1n76KOP1LNnTy9XBm/hbQEAAAAAFhATE3PZx6zZbDZde+21+tvf/ubFquBNBHQAAAAAhjCD7hlJSUllBvQaNWqoUaNGuvnmm9WvXz/ZbDYTqoM38LYAAAAAYEihr02FvpULiYW+Dknlzxb/Fr322mtmlwCT0ZoBAAAAAAALYAYdAAAAgCFFfn4q8qvcDHqRn0NSgXsKAqoJAjoAAAAAQ4p8fVVUySXuRb4E9F9q0qRJhcY5HA5lZWUpMzNT3bt3d36Nqo+ADgAAAAAW4OPjc9ku7iUubRJX0XNQNRDQAQAAABhSLF8VqXIz6MU0iCvl0KFDhsY3adLE8DmwNgI6AAAAAEMK5avCSgb0QgI6UApd3AEAAAAAsABm0AEAAAAYUiRfFVVyrq9IxW6qpvqYPn16hcY5HA5NmzbNs8XAFAR0AAAAAIa4J6BXbol8dTRjxoxyj5U0grPZbAT0aoyADgAAAAAWsHHjxjL32+12paWlacGCBUpISNAjjzzi5crgLQR0AAAAAIYwg+4ZnTt3LvdYnz59NHToULVu3VoDBgzwYlXwJprEAQAAADDkYkCv/AZjmjVrpn79+mnOnDlmlwIPIaADgJt98MEH6tmzp66++mrVrFlTTZo00ZAhQ3TkyBGzSwMAAFVcw4YNdeDAAbPLgIewxB0A3MThcOiBBx7QokWL1LRpUw0ePFh16tTRsWPHtHnzZmVnZysiIsLsMgEAqLQi+aqQJe5el5eXpzVr1ig0NNTsUuAhBHQAcJOXX35ZixYt0kMPPaSXX35Zvr6uS/cKCwtNqgwAAPcqkh+PWfOAlJSUMvcXFhbq6NGjWr58ub799ls9+uijXq4M3mJzlPTrBwD8aufPn1fjxo1Vt25dHThwQH5+lfv7p91uV0hIiKSnJAW4pUYAQFWRJ2mWcnJyFBwcbHYxLko+n77IaaKg4MoF9Fx7seJDMi35Os3i6+srh8Mhm63s1QU+Pj4aPHiwFi9erJo1a3q5OngDM+gA4Abr1q3TTz/9pJEjR6qoqEgfffSRvv32W4WGhqpHjx5q1qyZ2SUCAOA2RfKpdJO3IjfVUp0sXbq0zP2+vr4KDQ1VmzZtFBYW5uWq4E0EdABwg6+++krSxQ/Qm266Sd9++63zmI+Pj8aNG6e5c+eWe35eXp7y8vKcX9vtds8VCwBAJbmjCzsBvbTExESzS4DJ6OIOAG7w/fffS5Lmz5+vkJAQ7dy5U2fPnlVqaqquvfZazZs3T6+88kq558+cOVMhISHOjWZyAAAAvz0EdABwg+Lii41u/P39tXr1asXFxSkoKEi33Xab3n33Xfn4+GjevHnlnj9p0iTl5OQ4Nx7JBgCwskL5umVDaZ988oluv/12XXfddfrd736nbdu2OY/l5+fr5MmTKigoMLFCeBIBHQDc4GJDN6lt27Zq1KiRy7GWLVsqJiZGhw4d0pkzZ8o8PyAgQMHBwS4bAABWVSy//3Zy//VbMXfblvLBBx/o97//vT7//HP9+OOPWrNmjRISErR3715J0qlTp9S4cWMtXLjQ5ErhKQR0AHCD5s2bS1K5zyUt2X/+/HkvVQQAAKqaF154QaGhodq3b59++OEHbdmyRTabTbNnz5YkhYeHq3Pnzvrkk09MrhSeQkAHADfo1q2bJOnrr78udaygoEAHDx5UYGCg6tev7+3SAABwu5ImcZXdfo1du3bpzjvvVGhoqAIDAxUfH6+VK1dW+PxDhw5p2rRp6tu3rxo3biybzabo6OhfVYu7ffPNN7rnnnt03XXXSZLat2+vXr16aevWrc4xzZs31//7f//PrBLhYQR0AHCDpk2bqlevXjp48KCWLFnicmzWrFk6c+aMBgwYUOnnowMAYAVmBfSNGzeqY8eO2rJli+655x498MADOnHihAYNGnTZXi+X+te//qXp06frn//8p6666ir5+FgnEtWuXbvU881jY2N14sQJ59d169Z1NqdF9cNvigDgJn//+9/VoUMHjR49WqtXr1aLFi20e/dubdiwQVFRUZozZ47ZJQIAUGUVFhZq9OjR8vHxUWpqqlq1aiVJeuaZZ9SuXTtNnjxZ//d//6eoqKjLXqdz587avn27br75ZtWqVatUIDZTQkKCNm3apOLiYucfDoKCglyawmVkZDh736D6sc6fiwCgimvatKm+/PJLJSUl6auvvtLLL7+s9PR0/eEPf9DOnTsVFhZmdokAALhFkXzcMINuLIps2LBBhw4d0tChQ53hXLrYqHXy5MnKz89XSkrKFa8TExOj+Ph41apVy+jL9rgXXnhBx44d08MPP+zsW2Oz2ZzHv/76a3300Ufq1KmTWSXCw5hBBwA3ioiI0NKlS80uAwAAj3LHY9IK5TA0ftOmTZKkXr16lTrWu3dvSdLmzZsrVZPZpk+frpYtW2rhwoVatWqVWrdu7Vze3q9fP23cuFGSNGXKFDPLhAcR0AEAAABYXnp6uqSL92T/UlhYmIKCgpxjqqo333xTDodDNptNP/74o9atWyfp4iz6unXr1LlzZ73wwgsuKwhQvRDQAQAAABhS8izzyl3jIrvd7rI/ICBAAQEBpcbn5ORIUrn3XwcHBzvHVFUZGRll7g8ICFC9evXk61u5VQuwPgI6AAAAAEOKK/GYtP9d4+IS94iICJf9U6dO1bRp0yp17aoqMjLS7BJgMgI6AAAAANMcOXJEwcHBzq/Lmj2X/jdzXt4sud1uV926dd1fIOBFBHQAAAAAhvza55i7XuPiDHpwcLBLQC9Pyb3n6enpatOmjcuxEydOKDc3V+3atatUTWarSBf6EiNGjPBgJTALAR0AAACAIYXycUMX92JD47t06aKZM2dq3bp1Gjx4sMuxtWvXOsdUZffdd5+zSdylHI7/dby32WxyOBwE9GqKgA4AAH4VR6vpZpfwq9n2TDW7BAAGJSQkKCYmRsuWLdOjjz7q7GSek5OjF154Qf7+/kpMTHSOP378uHJychQeHl5uYzmrKe9RrXa7XWlpaXr77bfVp08f9evXz8uVwVsI6AAAAAAMcU8Xd2PPQffz89OSJUvUu3dvde7cWYMHD1adOnW0atUqZWdna+7cuYqOjnaOnzRpklJSUrR06VIlJSU59//444964oknnF8XFBToxx9/dBkzd+5c1atX79e+tF/t0j8wlGXMmDHq2rWrxo4d66WK4G0EdAAAAACGuOcedGNL3CWpW7du2rJli6ZOnaoVK1aooKBAN954o2bPnq1BgwZV6Bq5ubml7vU+d+6cy75p06aZEtCvJD4+Xn379tX06dPVq1cvs8uBBxDQAQAAAFQZ7dq106effnrFccnJyUpOTi61Pzo62uWe7qomOjpan3zyidllwEMI6AAAAAAMMWsG/bfO4XAoNTVVtWrVMrsUeAgBHQAAAIAhRfKtdBd3AnppmzdvLnN/YWGhjh49qpSUFO3cufOK96qj6iKgAwAAAIAFdO/e/bLL7202mzp16qT58+d7sSp4EwEdAAAAgCFmdHH/LZg2bVqZAd3X11ehoaGKi4tTu3btTKgM3kJABwAAAGBIkXzccA96kZuqqT6mTJlidgkwmY/ZBQAAAAAAKubf//53qcfEofpgBh0AAACAIe7p4l6586uzc+fO6aefflJRUelVBsnJyVqwYIG6dOkim80mSYqKivJ2ifAQAjoAAAAAQwjonnHw4EHde++92rVr1xWbxcXExDi/Li6mI351QUAHAAAAAAsYN26cvvzyS3Xp0kURERHy9S39R4w9e/Zoz549SkpK8n6B8DgCOgAAAABD3PMcdGbQf2nbtm0aOnSo3njjjXLHzJgxQ3v27NHrr7/uxcrgLQR0AAAAAIa45zFrLMv+pZycHMXGxppdBkxEF3cAAAAAsACHwyEfnytHtJLmcKh+mEEHAAAAYAhN4jxj48aNatKkyWXHJCUlqWvXrt4pCF5HQAcAAABgSJF83BDQWcz7S507d77imMjISEVGRnqhGpiBdwUAAAAAABbADDoAAAAAQwrd0MW9sucD1REBHQAAAIAh7uniXuSmaoDqgyXuAAAAAABYADPoAAAAAAwpdkMX92KWuAOlMIMOAAAAAIAFMIMOAAAAwBCeg+5Z+/fv18KFC7Vz506dOXNGISEhiouL0wMPPKCbbrrJ7PLgQQR0AAAAAIbwHHTP+dvf/qbHH39cxcXFkiSbzSaHw6G0tDQtWrRI8+bN02OPPWZylfAU3hUAAAAAYAGpqal67LHHFBYWptdff11Hjx5VYWGhjh49qqVLl6pRo0YaP368UlNTzS4VHsIMOgAAAABDCuUrX56D7nbz5s1TrVq1tHXrVkVGRjr3h4eHKzExUd26ddP111+v+fPnq3PnziZWCk8hoAMAAAAwxD3PQSeK/NK2bdvUv39/l3B+qYiICA0YMEDr1q3zcmXwFpa4AwAAAIAF5ObmKioq6rJjIiMjdebMGe8UBK/jz1YAAAAADOE56J7RoEEDnThx4rJjvvvuOzVo0MBLFcHbCOgAAAAADOExa54RHx9/xQZwW7ZsUVxcnJcqgrexxB0AAAAALGDcuHHq06ePjh8/XubxEydO6Pe//73Gjx/v5crgLcygAwAAADCkUL7yoYu728XHxys+Pr7c42FhYZo/f74XK4K3EdABAAAAGHJxiXtlu7gT0K/k7NmzstvtCg4OVp06dcwuB17AEncAAAAAsIi8vDw9++yzatasmUJCQhQREaGQkBA1a9ZMM2bMUF5entklwoOYQQcAAABgCE3iPOP8+fPq1q2bdu7cKX9/f11//fVq2LChTp48qYMHD2ratGn65JNPtHHjRtWuXdvscuEBzKADAAAAMKQkoFd2g6vZs2dr165dSkxMVGZmpvbv36/169dr//79ysjI0PDhw7Vr1y7Nnj3b7FLhIQR0AAAAALCAFStWqE2bNkpOTlZ4eLjLsUaNGiklJUVt2rTRypUrTaoQnkZABwAAAGBIsRtmz4uZQS8lKytLCQkJlx3Ts2dPZWVleacgeB33oAMAAAAwpFC+svGYNberXbu2Tp48edkxJ0+e5P7zaowZdAAAAACwgA4dOmjlypX697//XebxvXv3asWKFWrfvr2XK4O3MIMOAAAAwJAi+cqH56C73eTJk7V27VrdeuutGjFihLp166awsDAdP35cGzduVEpKioqKijR58mSzS4WHENABAAAAGHIxoPOYNXdr37693nnnHY0ePVqLFy/WokWLnMdsNptCQ0O1aNEidejQwcQq4UkEdAAAAACwiLvuuku9e/fW6tWrtWfPHtntdgUHB6tVq1bq37+/AgMDzS4RHkRABwAAAGAIM+ieFRgYqGHDhmnYsGFmlwIvI6ADAAAAMIQu7oBnENABAMCvYtsz1ewSAKBamT59eoXGORwOTZs2zbPFwBQEdAAAAACGFMtPRZWMEsVEkVJmzJhR7jGHwyHpYrM4Anr1xbsCAAAAgCFFbljizj3opW3cuLHM/Xa7XWlpaVqwYIESEhL0yCOPeLkyeAsBHQAAAAAsoHPnzuUe69Onj4YOHarWrVtrwIABXqwK3uRjdgEAAAAAqpYi+ahIvpXciCJGNWvWTP369dOcOXPMLgUewrsCAAAAgCGF8nXLBuMaNmyoAwcOmF0GPISADgAAAABVQF5entasWaPQ0FCzS4GHcA86AAAAAEOK5CdbJaNEZbvAV0cpKSll7i8sLNTRo0e1fPlyffvtt3r00Ue9XBm8hXcFAHjQ7Nmz9dRTT0mStm/frvj4eJMrAgCg8or/ex95Za8BV/fdd58cDodsNluZx318fDRs2DDNmjXLy5XBWwjoAOAh+/fv19SpUxUYGKhz586ZXQ4AALC4pUuXlrnf19dXoaGhatOmjcLCwrxcFbyJgA4AHlBQUKARI0aoVatWio2N1VtvvWV2SQAAuM3F2XOeg+5uiYmJZpcAk9EkDgA84Pnnn9d//vMfvf766/L15RcQAED1UvlHrFV+iTxQHTGDDgBulpaWpueff14zZszQ9ddfb3Y5AACgitm/f78WLlyonTt36syZMwoJCVFcXJweeOAB3XTTTWaXBw8ioAOAG+Xl5SkxMVGtWrXSH//4R7PLAQDAIwrlI0ell7izmLcsf/vb3/T444+ruLhYkmSz2eRwOJSWlqZFixZp3rx5euyxx0yuEp7CuwIA3OiZZ55Renq6li5damhpe15enux2u8sGAIBVFcnPLRtcpaam6rHHHlNYWJhef/11HT161PmItaVLl6pRo0YaP368UlNTzS4VHkJABwA32b59u+bOnaunn35aLVu2NHTuzJkzFRIS4twiIiI8VCUAALCqefPmqVatWtq6datGjBih8PBwSVJ4eLgSExO1detW1a5dW/Pnzze5UngKAR0A3KCwsFAjRozQTTfd5HzuuRGTJk1STk6Oczty5IgHqgQAwD1oEucZ27ZtU//+/RUZGVnm8YiICA0YMEBffPGFlyuDt7CuBADcIDc3V+np6ZIkf3//Mse0b99ekvTBBx+of//+LscCAgIUEBDg0RoBAHCXYjc8Zq2YgF5Kbm6uoqKiLjsmMjJSZ86c8U5B8DoCOgC4QUBAgO6///4yj6Wmpio9PV19+/ZV/fr1FR0d7d3iAABAldCgQQOdOHHismO+++47NWjQwEsVwdsI6ADgBrVq1dKSJUvKPJaUlKT09HRNmjRJ8fHxXq4MAAD3K5SvfJhBd7v4+PgrNoDbsmWL4uLivFQRvI2ADgAAAMCQIvnKUckoQUAvbdy4cVq5cqWOHz/ubBB3qRMnTuj3v/+9Bg4caEJ18AYCOgAAAABYQHx8/GVX24WFhdHBvZqjizsAeFhycrIcDgfL2wEA1QZd3L0jNzdXx44d0/nz580uBV5CQAcAAABgCAHds1599VXdcMMNCg4O1jXXXKOgoCDdfPPNWrx4sdmlwcMI6AAAAABgEffdd58efPBBHT58WG3btpXNZlPjxo2VmZmpsWPHaujQoWaXCA8ioAMAAAAwpKjY1y0bXL333ntKSUnR//3f/+m7777Tzp07JUn333+/fvjhB91///1655139Prrr5tcKTyFgA4AAADAkKJCXxVWcisqJKD/0uLFixUWFqbk5GQFBwe7HAsICNDChQt17bXXEtCrMQI6AAAAAFjAV199pR49eqh27dplHvf19VX37t21b98+L1cGb+ExawAAAAAMKSr0k62wclHCUcnzq6Pc3Fw1atTosmP8/f1VWFjopYrgbbwrAAAAABhSVOgjWyWXqDsKWcz7Sw0aNNDJkyfLPV5YWKjPPvtM119/vRergjfxrgAAAAAAC4iNjdU333xTan9+fr7279+vgQMH6uuvv9ZDDz1kQnXwBgI6AAAAAEOK/tvkrbIbXPXt21c7d+5UVlaWy/6ZM2fqpptu0scff6wnnnhCI0eONKdAeBxL3AEAAAAYUljoK1tBZZe4E9B/afjw4QoPD5fNZnPu69atm2rXrq0WLVpoyJAhuuWWW0ysEJ5GQAcAAAAAC7jqqqt0zz33uOz7/PPPTaoGZiCgAwAAADDEUeQnR1Elo0RlzweqId4VAAAAAIwp9L24VfYacDF9+vQKjXM4HJo2bZpni4EpCOgAAAAAYAEzZswo95jD4ZAk2Ww2Ano1RkAHAAAAYAwz6B6xcePGMvfb7XalpaVpwYIFSkhI0COPPOLlyuAtBHQAAAAAxhTZpELblcdd6Rpw0blz53KP9enTR0OHDlXr1q01YMAAL1YFb+I56AAAAABQBTRr1kz9+vXTnDlzzC4FHkJABwAAAGBMoZs2GNawYUMdOHDA7DLgISxxBwAAAGCMOwI2Ad2wvLw8rVmzRqGhoWaXAg8hoAMAAACABaSkpJS5v7CwUEePHtXy5cv17bff6tFHH/VyZfAWAjoAAAAAY5hB94j77rtPDodDNlvZDfR8fHw0bNgwzZo1y8uVwVsI6AAAAACMKZRU4IZrwMXSpUvL3O/r66vQ0FC1adNGYWFhXq4K3kRABwAAAAALSExMNLsEmIwu7gAAAACMKXLTBsNSU1M1ffp0s8uAhzCDDgAAAMAY7kH3mOLiYn3zzTc6ffq0iopK/xUjJSVFycnJ6tKli/Ne9S5duni7THgIAR0AAAAALGDXrl0aMGCAjh07dtlxNptN3bt3d35dXFzs6dLgJQR0AAAAAMYwg+4REydO1OnTp5WUlKSIiAj5+vqWGrNp0yZt2rRJ06ZN836B8DgCOgAAAABjCOgekZaWpsTERC1cuPCy4zZt2qRnnnnGS1XBm2gSBwAAAAAWcPbsWV1zzTVmlwETEdABAAAAGFOk/82i/9rtV3Zx37Vrl+68806FhoYqMDBQ8fHxWrlypaFr5OXlacaMGYqNjVXNmjXVqFEjjRkzRt9//32psT///LPmzZunoUOHqkWLFvLx8ZHNZlNWVtavewFXUNL4rbJjUDWxxB0AAACAMSYtcd+4caN69+6tmjVravDgwapTp45WrVqlQYMG6ciRI5owYcIVr1FcXKx+/fpp7dq1io+P18CBA5Wenq4lS5Zo/fr1+uKLL1S/fn3n+O+//15PPPGEJCkqKkp169bV6dOnjRdfAZmZmapbt+5lxzz++ONKSkryyPeH+ZhBBwAAAGB5hYWFGj16tHx8fJSamqpFixZp3rx52rt3r6699lpNnjxZ2dnZV7xOSkqK1q5dqyFDhmjbtm2aNWuWVq1apb///e/KyMjQ008/7TK+Xr16WrdunU6dOqWsrCzFxcV56iUqMjJSderUueyY4OBgRUZGeqwGmIsZdAAAAADGmDCDvmHDBh06dEgjR45Uq1atnPtDQkI0efJkJSUlKSUl5YrN0xYvXixJmjlzpstS8bFjx2rOnDl6++23tWDBAtWqVUuSFBQUpJ49exor9lcaOXJkhcY5HA4lJyd7thiYgoAOAAAAwJiC/26VvYYBmzZtkiT16tWr1LHevXtLkjZv3nzZa1y4cEE7duxQ8+bNFRUV5XLMZrOpZ8+eevXVV/Xll1/qtttuM1agG7z55ptyOBzlHnc4HLLZbAT0aoyADgAAAMDy0tPTJUmxsbGljoWFhSkoKMg5pjyHDh1ScXFxmde49Nrp6emmBPSMjIxS+xwOh06ePKldu3bpxRdf1I033qi//vWvXq8N3kFABwAAAGBMkX51F3aXa0iy2+0uuwMCAhQQEFBqeE5OjqSLS9rLEhwc7BxTnopc49Jx3lbeveVRUVFq166d+vfvr1atWmnVqlUVaoiHqocmcQAAAACMceNj1iIiIhQSEuLcZs6c6dWXUpU0btxYv/vd7/T666+bXQo8hBl0AAAAAKY5cuSIc+ZaUpmz59L/Zr3Lm9222+1XfERZRa5x6TgrCg8P16FDh8wuAx7CDDoAAAAAYyo7e35JF/jg4GCXrbyAfun94b904sQJ5ebmlntveYmYmBj5+PiUe6/65e5zt4oGDRpo6NChZpcBDyGgAwAAADDGjQG9orp06SJJWrduXalja9eudRlTnlq1aqldu3Y6cOBAqWemOxwOffbZZwoMDFTbtm2NFedmJY3hiopK3+g/btw4lrhXYwR0AAAAAJaXkJCgmJgYLVu2THv27HHuz8nJ0QsvvCB/f38lJiY69x8/flzffPNNqeXsY8aMkSRNmjTJ5ZFmr776qjIyMjRs2DDnM9C9raioSE8++aSCgoIUHh6uevXqad68ec7jP/74o1asWKHjx4+bUh88j3vQAQAAABjzK2bAy7yGAX5+flqyZIl69+6tzp07a/DgwapTp45WrVql7OxszZ07V9HR0c7xkyZNUkpKipYuXaqkpCTn/hEjRmjFihVavny5MjMz1aVLFx08eFDvv/++mjRpoueee67U937iiSf0448/SpL27dvn3BcUFCRJGjVqlDp16mTsBZVhzpw5mj9/vqKjo9W6dWtt375dTz75pJo0aaK77rpLtWvX1ujRozVu3DhNnz690t8P1kNABwAAAGBMSRf3yl7DoG7dumnLli2aOnWqVqxYoYKCAt14442aPXu2Bg0aVKFr+Pj46MMPP9SsWbP05ptv6qWXXtJVV12l+++/X88995zq169f6pz33nuv1JL4VatWOf/dtWtXtwT05ORkNW3aVHv37lWtWrV06tQp3XDDDXr11VedAT0hIUHr1q0joFdTNsel6zoAAJZgt9v/20H2KUllN8sBAFRXeZJmKScnx6W7uRU4P59m50g1K1nbBbs0McSSr9MstWvX1oMPPuiyrH348OFau3atvv/+e0kX70FftmyZTp48aVaZ8CBm0AEAAAAYY8IS99+Chg0bOh/1VuKaa67RTz/95Py6Ro0a5T4mDlUfTeIAAAAAGFPgpg0u7rnnHn366acugbxWrVou3dx3796t8PBwM8qDFxDQAQAAAMACpk6dqujoaN1xxx3atWuXyzGHw6G33npL69ev15133mlShfA0lrgDAAAAMKZIv6rJW6lrwMUNN9yggoICHTt2TLfeeqtq1aolf39/SVLdunVlt9vVtGlTTZs2zdxC4TEEdAAAAADGcA+6R/j4+Mjf39/lcXGSFBYWpkaNGqlHjx56+OGHVadOHXMKhMcR0AEAAADAAg4dOmR2CTAZAR0AAACAMSY9Bx2o7mgSBwAAAMCYQjdtMOzll19WkyZNzC4DHsIMOgAAAABYxJo1a7RhwwadPn3a5fFqJfbs2aPs7GwlJSXJZrPJ4XAoOTnZ+4XCIwjoAAAAAIwpUOXX4vIc9FLmzZunJ598Ujab7bLjbDab3nzzTUkioFczBHQAAAAAxvCYNY9YvHixYmNjtXjxYkVERMjX17fUmAULFmjBggXKysryfoHwOAI6AAAAAFhAVlaWxo8fr86dO5c7JjQ0VJIUGRnppargTQR0AAAAAMbQxd0jCgoKFBgYaHYZMBEBHQAAAIAxhar8Peh0cS+lS5cuioqKuuyY6Ohode3a1TsFwesI6AAAAABgARs2bLjimMTERCUmJnqhGpiBgA4AAADAmAJJl280XrFrAHBBQAcAAABgDF3cPaJJkyYVGudwOJSVlaXMzEx1797d+TWqPgI6AAAAAFiAj4+PHA7HFcdd+pz0ip6DqoGADgAAAMAYmsR5xKFDhwyNb9KkieFzYG2VfVsBACR99913WrBggXr16qXIyEj5+/srLCxMAwcO1I4dO8wuDwAA9yp5zFplNpa4A6UQ0AHADf7yl79o3LhxysjIUK9evTRhwgR16tRJH374oTp06KAVK1aYXSIAAKhizp49q++++05nz541uxR4CUvcAcAN2rVrp02bNqlLly4u+//1r38pISFBDz74oPr376+AgACTKgQAwI3c0YGdLu5lysvL04svvqiUlBRlZGQ498fExCgxMVETJ07k94lqzOagowAAeFTv3r21bt067dq1S23btq3QOXa7XSEhIZKeksSHMAD8tuRJmqWcnBwFBwebXYwL5+dTnxypRiVrK7BLH4dY8nWa5fz58+rWrZt27twpf39/NWvWTA0bNtTJkyd18OBB5efnKy4uThs3blTt2rXNLhcewBJ3APCwGjVqSJL8/Fi0BAAAyjd79mzt2rVLiYmJyszM1P79+7V+/Xrt379fGRkZGj58uHbt2qXZs2ebXSo8hBl0APCgw4cP69prr9VVV12lI0eOyNfXt8xxeXl5ysvLc35tt9sVEREhZtAB4LeoCsyg3+6mGfQ1zKBf6rrrrlOdOnW0c+fOcsfExcUpNzdXX3/9tRcrg7cwgw4AHlJQUKDhw4crLy9Ps2fPLjecS9LMmTMVEhLi3C6GcwAALKqyHdxLNrjIyspSQkLCZcf07NlTWVlZ3ikIXkdABwAPKC4uVlJSklJTUzV69GgNHz78suMnTZqknJwc53bkyBEvVQoAAKyidu3aOnny5GXHnDx5kvvPqzECOgC4WXFxse677z4tW7ZM9957rxYuXHjFcwICAhQcHOyyAQBgWYW62IW9Mhsz6KV06NBBK1eu1L///e8yj+/du1crVqxQ+/btvVwZvIWORQDgRsXFxRo5cqTeeOMNDRkyRMnJyfLx4W+hAIBqpkiSzQ3XgIvJkydr7dq1uvXWWzVixAh169ZNYWFhOn78uDZu3KiUlBQVFRVp8uTJZpcKDyGgA4CbXBrOBw0apDfffPOy950DAABcqn379nrnnXc0evRoLV68WIsWLXIes9lsCg0N1aJFi9ShQwcTq4QnEdABwA1KlrW/8cYbuvvuu/XWW28RzgEA1Zc7lqezxL1Md911l3r37q3Vq1drz549stvtCg4OVqtWrdS/f38FBgaaXSI8iIAOAG4wY8YMpaSkKCgoSNdee62ee+65UmP69++vVq1aeb84AADcjYDuUYGBgRo2bJiGDRtmdinwMgI6ALhByeNOcnNz9fzzz5c5Jjo6moAOAACAchHQAcANkpOTlZycbHYZAAB4R6EkRyWvQZO4UqZPn16hcQ6HQ9OmTfNsMTAFAR0AAACAMe4I1wT0UmbMmFHuMYfj4l9EbDYbAb0aI6ADAAAAgAVs3LixzP12u11paWlasGCBEhIS9Mgjj3i5MngLAR0AAACAMSxx94jOnTuXe6xPnz4aOnSoWrdurQEDBnixKniTj9kFAAAAAKhiCt20wZBmzZqpX79+mjNnjtmlwEMI6AAAAABQRTRs2FAHDhwwuwx4CEvcAQAAABhTKKm4kteo7Pm/QXl5eVqzZo1CQ0PNLgUeQkAHAAAAYEyRKn8POgG9lJSUlDL3FxYW6ujRo1q+fLm+/fZbPfroo16uDN5CQAcAAAAAC7jvvvvkcDhks9nKPO7j46Nhw4Zp1qxZXq4M3kJABwAAAGBMoSrfzYoZ9FKWLl1a5n5fX1+FhoaqTZs2CgsL83JV8CYCOgAAAABjCOgekZiYaHYJMBld3AEAAAAAsABm0AEAAAAYUyBm0AEPIKADAAAAMKZYle/iXtnzgWqIJe4AAAAAAFgAM+gAAAAAjCmUVPaTwCqOGXSgFAI6AAAAAGMI6IBHsMQdAAAAAAALYAYdAAAAgDEFYgYd8AACOgAAAABjikRABzyAJe4AAAAAAFgAM+gAAAAAjGMGHHA7ZtABAAAAALAAAjoAAAAAABZAQAcAAAAAwAII6AAAAAAAWAABHQAAAAAAC6CLOwAAAACDCv67VfYaAC5FQAcAAABgUOF/t8peA8ClWOIOAAAAAIAFMIMOAAAAwCCWuAOeQEAHAAAAYBBL3AFPYIk7AAAAAAAWwAw6AAAAAIMKVfkl6sygA79EQAcAAABgEPegA57AEncAAAAAACyAGXQAAAAABtEkDvAEAjoAAAAAg7gHHfAElrgDAAAAAGABzKADAAAAMIgl7oAnENABAAAAGEQXd8ATWOIOAAAAAIAFMIMOAAAAwCCWuAOeQEAHAAAAYBBd3AFPYIk7AAAAAAAWwAw6AAAAAINY4g54AgEdAAAAgEF0cQc8gSXuAAAAAABYADPoAAAAAAxiiTvgCQR0AAAAAAbRxR3wBJa4AwAAAABgAcygAwAAADCIJe6AJxDQAQAAABhEF3fAE1jiDgAAAACABTCDDgAAAMAgZtABTyCgAwAAADCIe9ABT2CJOwAAAAAAFsAMOgAAAACDeA464AkEdAAAAAAGscQd8ASWuAMAAAAAYAHMoAMAAAAwqECVjxJ0cQd+iYAOAAAAwCCWuAOewBJ3AHCjXbt26c4771RoaKgCAwMVHx+vlStXml0WAADVhjs+a/Py8jRjxgzFxsaqZs2aatSokcaMGaPvv/++3HPefvtttWvXToGBgapbt6769OmjtLS0Mse+9dZbGjt2rNq2bauAgADZbDYlJycbqhG/TcygA4CbbNy4Ub1791bNmjU1ePBg1alTR6tWrdKgQYN05MgRTZgwwewSAQBwE3O6uLvjs7a4uFj9+vXT2rVrFR8fr4EDByo9PV1LlizR+vXr9cUXX6h+/fou5zz//PN6+umnFRUVpQceeEBnz57VO++8ow4dOmj9+vXq2LGjy/inn35a2dnZqlevnsLDw5WdnW34teK3yeZwOBxmFwEAVV1hYaFatGiho0eP6osvvlCrVq0kSTk5OWrXrp2ysrL07bffKioqqkLXs9vtCgkJkfSUpACP1Q0AsKI8SbOUk5Oj4OBgs4tx8b/Pp+mSalbyahckTa3w63TXZ+3SpUt13333aciQIXr77bdls9kkSQsXLtSDDz6oMWPG6NVXX3WOT09P1/XXX6+YmBjt3Lnzv69f2rNnj+Lj4xUTE6P9+/fLx+d/i5M///xzxcbGKioqSrNmzdKkSZO0dOlSJSUlGfsR4TeHJe4A4AYbNmzQoUOHNHToUOcvDJIUEhKiyZMnKz8/XykpKeYVCABAFeeuz9rFixdLkmbOnOkM55I0duxYxcTE6O2339b58+ed+5cuXarCwkL96U9/coZzSWrVqpWGDBmir7/+Wlu2bHH5Hj169KjwH+WBSxHQAcANNm3aJEnq1atXqWO9e/eWJG3evNmbJQEA4EEFbtoqzh2ftRcuXNCOHTvUvHnzUgHaZrOpZ8+eOnfunL788ku3fl+gogjoAOAG6enpkqTY2NhSx8LCwhQUFOQcAwBA1Vfopq3i3PFZe+jQIRUXF5d5jUuvfel10tPTFRQUpLCwsAqNByqDJnEA4AY5OTmS5LL07VLBwcHOMWXJy8tTXl5eqetdvA8RAPDbcvH//dZuFeWOz6eL17Db7S57AwICFBBQuv9KZT9rK3qNS8eV/LtBgwYVHg9UBgEdACxg5syZmj59ehlHXvJ6LQAAazh16lS5QdIs/v7+CgsL04kT7vl8CgoKUkREhMu+qVOnatq0aW65PlDVENABwA1KfoEq7y/odrtddevWLff8SZMmafz48c6vz5w5o6ioKB0+fNhyv5xZld1uV0REhI4cOWK5rsdWxc/MOH5mxvEzMy4nJ0eRkZG66qqrzC6llJo1ayozM1P5+fluuZ7D4XBp1CapzNlzqfKftRW9xqXjSv5tZDxQGQR0AHCDS+9Ba9OmjcuxEydOKDc3V+3atSv3/PKW84WEhPALrUHBwcH8zAziZ2YcPzPj+JkZd+lju6ykZs2aqlmzso9YM66yn7WSFBMTIx8fn3LvGS/rPvfY2Fht375dJ06cKHUf+uXuiwd+DWu+6wGgiunSpYskad26daWOrV271mUMAAAwzh2ftbVq1VK7du104MABZWdnuxxzOBz67LPPFBgYqLZt27r1+wIVRUAHADdISEhQTEyMli1bpj179jj35+Tk6IUXXpC/v78SExPNKxAAgCrO6Gft8ePH9c0335Ranj5mzBhJF28vu7QR36uvvqqMjAwNGzZMtWrVcu4fOXKk/Pz89Pzzz7tca8+ePVq+fLmuu+46derUyd0vF79RLHEHADfw8/PTkiVL1Lt3b3Xu3FmDBw9WnTp1tGrVKmVnZ2vu3LmKjo6u8PUCAgI0derUcu/DQ2n8zIzjZ2YcPzPj+JkZx8+sbEY/aydNmqSUlBQtXbpUSUlJzv0jRozQihUrtHz5cmVmZqpLly46ePCg3n//fTVp0kTPPfecy/e99tprNW3aND399NO6+eabNXDgQJ09e1bvvPOOJGnx4sWlbkdYsmSJtmzZIknat2+fc1/JM9U7deqkUaNGufknhOrA5rD28xsAoErZuXOnpk6dqm3btqmgoEA33nijxo8fr0GDBpldGgAA1UJFP2uTkpLKDOjSxcebzpo1S2+++aaOHDmiq666Sn369NFzzz2nhg0blvl93377bS1YsED/+c9/5O/vr44dO+rZZ59V69atS40t+d7lGTFihJKTkw2/dlR/BHQAAAAAACyAe9ABAAAAALAAAjoAAAAAABZAQAcAAAAAwAII6ABgIbt27dKdd96p0NBQBQYGKj4+XitXrjS7LEv67rvvtGDBAvXq1UuRkZHy9/dXWFiYBg4cqB07dphdXpUye/Zs2Ww22Ww2ffHFF2aXY1kffPCBevbsqauvvlo1a9ZUkyZNNGTIEB05csTs0izH4XDo/fffV7du3RQeHq7atWurefPmGjt2rDIyMswuz1RvvfWWxo4dq7Zt2yogIEA2m+2yzcLsdrvGjx+vqKgoBQQEKDo6Wk8++aRyc3O9VzQAr6FJHABYxMaNG9W7d2/VrFmzzEfHTJgwwewSLeWpp57S7Nmz1bRpU3Xt2lX169dXenq6Vq9eLYfDoWXLltE9vwL279+vtm3bys/PT+fOndP27dsVHx9vdlmW4nA49MADD2jRokVq2rSpevfurTp16ujYsWPavHmz3n77bZ6B/AsTJkzQ/PnzFR4ern79+ik4OFh79+7VunXrFBQUpG3btqlly5Zml2mK6OhoZWdnq169egoMDFR2dnaZXcYl6dy5c+rUqZP27NmjXr166ZZbbtHu3bu1bt06xcXFKTU1VTVr1vT+iwDgOQ4AgOkKCgocTZs2dQQEBDh2797t3H/mzBnHtdde6/D393dkZWWZV6AFrVq1yrFp06ZS+1NTUx01atRw1K1b13HhwgUTKqs68vPzHa1bt3bceuutjnvvvdchybF9+3azy7KcBQsWOCQ5HnroIUdhYWGp4wUFBSZUZV3Hjx93+Pj4OKKiohxnzpxxOTZ//nyHJMfIkSNNqs58n332mfP/5zNnznRIcixdurTMsc8884xDkmPixIku+ydOnOiQ5HjhhRc8XS4AL2OJOwBYwIYNG3To0CENHTpUrVq1cu4PCQnR5MmTlZ+ff9nnqf4W3XXXXerSpUup/bfddpu6deumn376Sfv27TOhsqrj+eef13/+8x+9/vrr8vX1NbscSzp//rymT5+umJgY/fnPfy7z5+Tn52dCZdaVlZWl4uJidezYUSEhIS7H+vTpI0n64YcfzCjNEnr06KGoqKgrjnM4HFqyZImCgoI0ZcoUl2NTpkxRUFCQlixZ4qkyAZiEgA4AFrBp0yZJUq9evUod6927tyRp8+bN3iypSqtRo4YkgtPlpKWl6fnnn9fUqVN1/fXXm12OZa1bt04//fST+vfvr6KiIr3//vuaNWuWFi5cqIMHD5pdniXFxsbK399fW7duld1udzn28ccfS5ISEhLMKK1KSU9P17Fjx9SxY0cFBga6HAsMDFTHjh2VkZFBDwSgmuE3FwCwgPT0dEkXf7H9pbCwMAUFBTnH4PIOHz6szz//XOHh4brxxhvNLseS8vLylJiYqFatWumPf/yj2eVY2ldffSVJ8vX11U033aRvv/3WeczHx0fjxo3T3LlzzSrPkq6++mrNmjVLEyZMUIsWLVzuQd+wYYMeeughPfzww2aXaXmX+1wo2b927Vqlp6crIiLCm6UB8CACOgBYQE5OjiSVWg5aIjg42DkG5SsoKNDw4cOVl5en2bNns2y7HM8884zS09P11Vdf8TO6gu+//16SNH/+fLVu3Vo7d+7Uddddp927d2vMmDGaN2+emjZtqgcffNDkSq1l3Lhxaty4sUaNGqWFCxc693fq1ElDhw5ldUsFVORz4dJxAKoHlrgDAKqF4uJiJSUlKTU1VaNHj9bw4cPNLsmStm/frrlz5+rpp5/+zXbRNqK4uFiS5O/vr9WrVysuLk5BQUG67bbb9O6778rHx0fz5s0zuUrrmTFjhu69915NnjxZR44c0dmzZ/Wvf/1LFy5cUNeuXfXRRx+ZXSIAWBIBHQAsoGSGpLyZELvdXu4sCi6GqPvuu0/Lli3Tvffe6zJjh/8pLCzUiBEjdNNNN+mpp54yu5wqoeR917ZtWzVq1MjlWMuWLRUTE6NDhw7pzJkzJlRnTZ9//rmmTp2qhx9+WE899ZSuueYaBQUFqVOnTvrHP/6hGjVq8NjICqjI58Kl4wBUDwR0ALCAknsMy7rP/MSJE8rNzS33PsTfuuLiYo0cOVIpKSkaMmSIkpOT5ePDx1tZcnNzlZ6erj179sjf3182m825lTwloH379rLZbFq9erW5xVpE8+bNJUmhoaFlHi/Zf/78eS9VZH2ffvqpJKlbt26ljoWFhalFixY6ePCgcnNzvV1alXK5z4VL9/PZAFQv3AAEABbQpUsXzZw5U+vWrdPgwYNdjq1du9Y5Bq5Kwvkbb7yhQYMG6c033+Se6ssICAjQ/fffX+ax1NRUpaenq2/fvqpfv76io6O9W5xFlYTMr7/+utSxgoICHTx4UIGBgapfv763S7Os/Px8SeU/Su2HH36Qj4+P82kLKFtsbKwaNWqkrVu36ty5cy6d3M+dO6etW7eqSZMmNIgDqhmmGADAAhISEhQTE6Nly5Zpz549zv05OTl64YUX5O/vr8TERPMKtKCSZe1vvPGG7r77br311luE8yuoVauWlixZUubWoUMHSdKkSZO0ZMkStWrVytxiLaJp06bq1auXDh48WOqZ07NmzdKZM2c0YMAAmp5domPHjpIuNtb75fLshQsX6ujRo2rfvr0CAgLMKK/KsNlsGjVqlHJzc/Xss8+6HHv22WeVm5ur0aNHm1QdAE+xORwOh9lFAACkjRs3qnfv3qpZs6YGDx6sOnXqaNWqVcrOztbcuXO5Z/MXpk2bpunTpysoKEiPPfZYmQGpf//+BM0KSkpKUkpKirZv3674+Hizy7GUQ4cOqUOHDvr+++/1u9/9Ti1atNDu3bu1YcMGRUVF6YsvvlBYWJjZZVpGUVGRunfvrtTUVDVo0EB9+/ZVaGio0tLStGHDBtWqVUubNm1Su3btzC7VFEuWLNGWLVskSfv27VNaWpo6duyoZs2aSbrY6X7UqFGSLs6Ud+zYUXv37lWvXr3UunVrpaWlad26dYqLi9PmzZtVq1Yt014LAPcjoAOAhezcuVNTp07Vtm3bVFBQoBtvvFHjx4/XoEGDzC7NckoC5eUsXbpUSUlJ3imoiiOgX96RI0f0zDPPaM2aNTp16pTCwsLUt29fPfPMM2rQoIHZ5VlOXl6eXnrpJa1cuVIHDhxQfn6+GjZsqG7dumny5Mm67rrrzC7RNFf6f9eIESOUnJzs/DonJ0fTpk3TqlWrdOLECYWHh+vuu+/W1KlTVadOHS9UDMCbCOgAAAAAAFgA96ADAAAAAGABBHQAAAAAACyAgA4AAAAAgAUQ0AEAAAAAsAACOgAAAAAAFkBABwAAAADAAgjoAAAAAABYAAEdAAAAAAALIKADAAAAAGABBHQAAAAAACyAgA4AAAAAgAUQ0AEAAAAAsAACOgAAAAAAFkBABwAAAADAAgjoAAAAAABYAAEdAAAAAAALIKADAAAAAGABBHQAAAAAACyAgA4AAAAAgAUQ0AEAAAAAsAACOgAAAAAAFkBABwAAAADAAgjoAAAAAABYAAEdAAAAAAALIKADAAAAAGABBHQAAAAAACyAgA4AAAAAgAUQ0AEAAAAAsAACOgAAAAAAFvD/AQXalII2SCn7AAAAAElFTkSuQmCC' width=1000.0/>\n",
       "            </div>\n",
       "        "
      ],
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 5))\n",
    "img2 = ax.imshow(np.random.uniform(size=(10,10)), vmin=-3, vmax=0.5, cmap='jet')\n",
    "cbar = fig.colorbar(img2)\n",
    "cbar.ax.set_ylabel('Contaminant concentration', rotation=270, labelpad=12, fontsize=14)\n",
    "cbar.ax.set_yticks([-3, -2, -1, 0])\n",
    "cbar.ax.set_yticklabels([0.001, 0.01, 0.1, 1])\n",
    "cbar.ax.tick_params(axis='y', labelsize=14)\n",
    "\n",
    "\n",
    "def plot_solution(Time):\n",
    "    ax.clear()\n",
    "    ax.imshow(np.log10(sol.y[:, Time].reshape(N_x, N_y)), interpolation='nearest',\n",
    "                    cmap='jet', vmin=-3, vmax=0.5, origin='lower', extent=[0, lx, 0, ly])\n",
    "    ax.set_title(f'Time: {sol.t[Time]:.2f}', fontsize=16)\n",
    "    ax.tick_params(axis='both', labelsize=14)\n",
    "    # cbar.ax.yaxis.set_ticks_position('left')\n",
    "\n",
    "\n",
    "widgets.interactive(plot_solution, Time=widgets.IntSlider(min=0, max=(len(sol.t)-1), value=0))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close('all')"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "87df4c83b56fb8b58e6fe6434764317a79f3f72ebd797e1bfe209587b31ed721"
  },
  "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}