Computational_Seismology_Igel/dg_elastic_physicalfluxes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style='background-image: url(\"../../share/images/header.svg\") ; padding: 0px ; background-size: cover ; border-radius: 5px ; height: 250px'>\n",
    "    <div style=\"float: right ; margin: 50px ; padding: 20px ; background: rgba(255 , 255 , 255 , 0.7) ; width: 50% ; height: 150px\">\n",
    "        <div style=\"position: relative ; top: 50% ; transform: translatey(-50%)\">\n",
    "            <div style=\"font-size: xx-large ; font-weight: 900 ; color: rgba(0 , 0 , 0 , 0.8) ; line-height: 100%\">Computational Seismology</div>\n",
    "            <div style=\"font-size: large ; padding-top: 20px ; color: rgba(0 , 0 , 0 , 0.5)\">Discontinuous Galerkin Method with Physics based numerical flux - 1D Elastic Wave Equation </div>\n",
    "        </div>\n",
    "    </div>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook is based on the paper [A new discontinuous Galerkin spectral element method for elastic waves with physically motivated numerical fluxes](https://www.geophysik.uni-muenchen.de/~gabriel/kduru_waves2017.pdf)\n",
    "Published in the [13th International Conference on Mathematical and Numerical Aspects of Wave Propagation](https://cceevents.umn.edu/waves-2017)\n",
    "\n",
    "##### Authors:\n",
    "* Kenneth Duru\n",
    "* Ashim Rijal ([@ashimrijal](https://github.com/ashimrijal))\n",
    "* Sneha Singh\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic Equations ##\n",
    "\n",
    "The source-free elastic wave equation in a heterogeneous 1D medium is \n",
    "\n",
    "\\begin{align}\n",
    "\\rho(x)\\partial_t v(x,t) -\\partial_x \\sigma(x,t) & = 0\\\\\n",
    "\\frac{1}{\\mu(x)}\\partial_t \\sigma(x,t) -\\partial_x v(x,t) & = 0  \n",
    "\\end{align}\n",
    "\n",
    "with $\\rho(x)$ the density, $\\mu(x)$ the shear modulus and $x = [0, L]$. At the boundaries $ x = 0, x = L$ we pose the general well-posed linear boundary conditions\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "B_0(v, \\sigma, Z_{s}, r_0): =\\frac{Z_{s}}{2}\\left({1-r_0}\\right){v} -\\frac{1+r_0}{2} {\\sigma} = 0,  \\quad \\text{at} \\quad x = 0, \\\\\n",
    " B_L(v, \\sigma, Z_{s}, r_n): =\\frac{Z_{s}}{2} \\left({1-r_n}\\right){v} + \\frac{1+r_n}{2}{\\sigma} = 0, \\quad \\text{at} \\quad  x = L.\n",
    " \\end{split}\n",
    "\\end{equation}\n",
    "\n",
    "with the reflection coefficients $r_0$, $r_n$ being real numbers and  $|r_0|, |r_n| \\le 1$. \n",
    "\n",
    "Note that at $x = 0$,  while  $r_0 = -1$ yields a clamped wall, $r_0 = 0$  yields  an absorbing boundary, and  with $r_0 = 1$  we have a free-surface boundary condition. Similarly, at $x = L$, $r_n = -1$ yields a clamped wall, $r_n = 0$ yields an absorbing boundary, and  $r_n = 1$  gives a free-surface boundary condition.\n",
    "\n",
    "1) Discretize the spatial domain $x$ into $K$ elements and denote the ${k}^{th}$ element $e^k = [x_{k}, x_{k+1}]$ and the element width $\\Delta{x}_k = x_{k+1}-x_{k}$. Consider two adjacent elements  $e^k = [x_{k}, x_{k+1}]$ and  $e^{k+1} = [x_{k+1}, x_{k+2}]$ with an interface at $x_{k+1}$. At the interface we pose the physical conditions for a locked interface\n",
    "\n",
    "\\begin{align}\n",
    "\\text{force balance}:  \\quad &\\sigma^{-} = \\sigma^{+} = \\sigma, \\nonumber \\\\\n",
    "\\text{no slip}: \\quad & [\\![ v]\\!]  = 0,\n",
    "\\end{align}\n",
    "\n",
    "where $[\\![ v]\\!] = v^{+} - v^{-}$, and $v^{-}, \\sigma^{-}$ and $v^{+}, \\sigma^{+}$ are the fields in $e^k = [x_{k}, x_{k+1}]$ and  $e^{k+1} = [x_{k+1}, x_{k+2}]$, respectively. \n",
    "\n",
    "2) Within the element derive the weak form of the equation by multiplying both sides by an arbitrary test function and integrating over the element.\n",
    "\n",
    "3) Next map the $e^k = [x_{k}, x_{k+1}]$ to a reference element $\\xi = [-1, 1]$\n",
    "\n",
    "4) Inside the transformed  element  $\\xi \\in [-1, 1]$, approximate the solution  and material parameters by a polynomial interpolant,  and write \n",
    "\\begin{equation}\n",
    "v^k(\\xi, t) = \\sum_{j = 1}^{N+1}v_j^k(t) \\mathcal{L}_j(\\xi), \\quad \\sigma^k(\\xi, t)  = \\sum_{j = 1}^{N+1}\\sigma_j^k(t) \\mathcal{L}_j(\\xi),\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\rho^k(\\xi) = \\sum_{j = 1}^{N+1}\\rho_j^k \\mathcal{L}_j(\\xi), \\quad \\mu^k(\\xi) = \\sum_{j = 1}^{N+1}\\mu_j^k \\mathcal{L}_j(\\xi),\n",
    "\\end{equation}\n",
    "\n",
    "where $ \\mathcal{L}_j$ is the $j$th interpolating polynomial of degree $N$. If we consider  nodal basis then the interpolating polynomials satisfy $ \\mathcal{L}_j(\\xi_i) = \\delta_{ij}$.\n",
    "\n",
    "The interpolating nodes $\\xi_i$, $i = 1, 2, \\dots, N+1$ are the nodes of a Gauss quadrature with\n",
    "\n",
    "\\begin{equation}\n",
    " \\sum_{i = 1}^{N+1} f(\\xi_i)w_i \\approx \\int_{-1}^{1}f(\\xi) d\\xi,\n",
    "\\end{equation}\n",
    "\n",
    "where $w_i$ are quadrature weights.\n",
    "\n",
    "5) At the element boundaries $\\xi = \\pm 1$, we generate $\\widehat{v}^{k}(\\pm 1, t)$ $\\widehat{\\sigma}^{k}(\\pm 1, t)$ by solving a Riemann problem and constraining the solutions against interface and boundary conditions. Then numerical fluctuations $F^k(-1, t)$ and $G^k(1, t)$ are obtained by penalizing hat variables against the incoming characteristics only.\n",
    "\n",
    "6) Finally, the flux fluctuations are appended to the semi-discrete PDE with special penalty weights and we have \n",
    "\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "\\frac{d \\boldsymbol{v}^k( t)}{ d t} &= \\frac{2}{\\Delta{x}_k} W^{-1}({\\boldsymbol{\\rho}}^{k})\\left(Q \\boldsymbol{\\sigma}^k( t) - \\boldsymbol{e}_{1}F^k(-1, t)- \\boldsymbol{e}_{N+1}G^k(1, t)\\right),\n",
    "\\end{split}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "\\frac{d \\boldsymbol{\\sigma}^k( t)}{ d t} &= \\frac{2}{\\Delta{x}_k} W^{-1}(1/{\\boldsymbol{\\mu}^{k}})\\left(Q \\boldsymbol{v}^k( t)  + \\boldsymbol{e}_{1}\\frac{1}{Z_{s}^{k}(-1)}F^k(-1, t)- \\boldsymbol{e}_{N+1}\\frac{1}{Z_{s}^{k}(1)}G^k(1, t)\\right),\n",
    "\\end{split}\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "where \n",
    "\\begin{align}\n",
    "\\boldsymbol{e}_{1} = [ \\mathcal{L}_1(-1), \\mathcal{L}_2(-1), \\dots,  \\mathcal{L}_{N+1}(-1) ]^T, \\quad  \\boldsymbol{e}_{N+1} = [ \\mathcal{L}_1(1), \\mathcal{L}_2(1), \\dots,  \\mathcal{L}_{N+1}(1) ]^T,\n",
    "\\end{align}\n",
    "and\n",
    "\\begin{align}\n",
    "G^k(1, t):= \\frac{Z_{s}^{k}(1)}{2} \\left(v^{k}(1, t)-\\widehat{v}^{k}(1, t) \\right) + \\frac{1}{2}\\left(\\sigma^{k}(1, t)- \\widehat{\\sigma}^{k}(1, t)\\right), \n",
    "\\end{align}\n",
    "\\begin{align}\n",
    "F^{k}(-1, t):= \\frac{Z_{s}^{k}(-1)}{2} \\left(v^{k}(-1, t)-\\widehat{v}^{k}(-1, t) \\right) - \\frac{1}{2}\\left(\\sigma^{k}(-1, t)- \\widehat{\\sigma}^{k}(-1, t)\\right).\n",
    "\\end{align}\n",
    "\n",
    "And the weighted elemental mass matrix $W^N(a)$ and the stiffness matrix $Q^N $ are defined by\n",
    "\n",
    "\\begin{align}\n",
    "W_{ij}(a) = \\sum_{m = 1}^{N+1} w_m \\mathcal{L}_i(\\xi_m)  {\\mathcal{L}_j(\\xi_m)} a(\\xi_m), \\quad Q_{ij} = \\sum_{m = 1}^{N+1} w_m \\mathcal{L}_i(\\xi_m)  {\\mathcal{L}_j^{\\prime}(\\xi_m)}.\n",
    "\\end{align}\n",
    "\n",
    "7) Time extrapolation can be performed using any stable time stepping scheme like Runge-Kutta or ADER scheme.This notebook implements both Runge-Kutta and ADER schemes for solving the free source version of the elastic wave equation in a homogeneous media. To keep the problem simple, we use as spatial initial condition a Gauss function with half-width $\\delta$\n",
    "\n",
    "\\begin{equation}\n",
    "v(x,t=0)  = e^{-1/\\delta^2 (x - x_{o})^2}, \\quad \\sigma(x,t=0) = 0\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**** Exercises****\n",
    "1. Lagrange polynomial is used to interpolate the solution and the material parameters. First use polynomial degree 2 and then 6. Compare the simulation results in terms of accuracy of the solution (third and fourth figures give erros). At the end of simulation, time required to complete the simulation is also printed. Also compare the time required to complete both simulations.\n",
    "\n",
    "2. We use quadrature rules: Gauss-Legendre-Lobatto and Gauss-Legendre. Run simulations once using Lobatto and once using Legendre rules. Compare the difference.\n",
    "\n",
    "3. Now fix the order of polynomial to be 6, for example. Then use degree of freedom 100 and for another simulation 250. What happpens? Also compare the timre required to complete both simulations.\n",
    "\n",
    "4. Experiment with the boundary conditions by changing the reflection coefficients $r_0$ and $r_n$.\n",
    "\n",
    "5. You can also play around with sinusoidal initial solution instead of the Gaussian.\n",
    "\n",
    "6. Change the time-integrator from RK to ADER. Observe if there are changes in the solution or the CFL number. Vary the polynomial order N."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Parameters initialization and plotting the simulation\n",
    "# Import necessary routines\n",
    "import Lagrange\n",
    "import numpy as np\n",
    "import timeintegrate\n",
    "import quadraturerules\n",
    "import specdiff\n",
    "import utils\n",
    "import matplotlib.pyplot as plt\n",
    "import timeit                         # to check the simulation time\n",
    "\n",
    "#plt.switch_backend(\"TkAgg\")          # plots in external window\n",
    "plt.switch_backend(\"nbagg\")           # plots within this notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "code_folding": [
     0,
     16,
     24,
     42,
     47
    ],
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "points per wavelength:  4\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAALICAYAAABiqwZ2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdebhVddn/8fdHEEdUBBRkEDC0NJOMcCiVJ7UENcx8FMuxfHDIlF9amZn5mE2mppg5lKbmPKSSYg48WTY4ADKIOACiHBlFBRUn4P79sdaOdc7ZZ+/NGfY+Z+/P67r2tfda93etda91rbP3zZfvWksRgZmZmZmZJdardAJmZmZmZu2JC2QzMzMzswwXyGZmZmZmGS6QzczMzMwyXCCbmZmZmWW4QDYzMzMzy3CBbGbWQUk6W9LvW7iOAZJCUufWysvMrKOT74NsZla7JA0AXgbWj4hVlc2m9UgKYHBEzK50LmbW8bgH2cysA3KPr5lZ23GBbGZWgKR5ks6UNF3Sckm3S9owjR0n6R8N2oekj6Wfr5f0W0kPSnpH0j8l9ZJ0qaQ3JT0v6dOZZbeRdLekpZJelnRaJnaepLsk3SRpBXBcOu+mTJvPS/qXpLckzZd0XDr/QEnPSFqRzj+vxH0/S9JdDeZdJmlcZv/nSno7zffrTaxnvXRdcyQtk3SHpC3T2BHpOjZLp0dIWiSpZ2Z789PcJ0vaK7PeTukwkzlpDpMl9ZP097TJtPS4H1HK/pqZ5bhANjMr7nDgAGAg8CnguHVc9hygB/AB8G9gSjp9F3AJJEUk8GdgGtAH2BcYK+lLmXWNSpfZArg5uxFJ/YEHgcuBnsAQYGoafhc4Jl3uQOBkSYeUkPutwMhM8dop3Z9bJG0CjANGRERXYM/M9ho6DTgE2AfYBngTuAIgIm5Pj8k4Sd2Ba4ETImJpuuzT6b5sCdwC3Jn7BwrwHeBIYCSwGfANYGVE7J3Gd4mITdNtmJmVzAWymVlx4yJiQUS8QVLEDlmHZe+JiMkR8T5wD/B+RNwYEauB24FcD/JngZ4RcX5EfBgRc4HfAaMz6/p3RNwbEWsi4r0G2/k68GhE3BoRH0XEsoiYChARj0XEjHS56SSF7z7FEo+IV0iK+Vwx/QWSAvSJdHoN8ElJG0XEwoiY2cSqTgR+GBF1EfEBcB5wWGaYyLfSdT8G/Dki7s/kcFO6L6si4mJgA2CHNHwCcE5EvBCJaRGxrNh+mZkV4wLZzKy4RZnPK4FN12HZxZnP7+WZzq1rW2CbdHjEW5LeAs4Gts60n19gO/2AOfkCknaT9Nd06MZy4CSSHuxS3ELSSwvwtXSaiHgXOCJd10JJD0j6eBPr2Ba4J7Nfs4DVpPsWEW8BdwKfBC5ukPsZkmalw1veAjbP5N7kPpuZtYQLZDOz5nsX2Dg3IalXC9Y1H3g5IrbIvLpGxMhMm0K3HZoPbNdE7BZgPNAvIjYHrgJUYl53AsMl9QW+kq4rSSbioYjYH+gNPE/S491UbiMa7NuGEfEagKQhJMMjbiUZtkE6fy/g+yTDOrpFxBbA8kzuhfbZzKzZXCCbmTXfNGAnSUPScbHntWBdTwErJH1f0kbpBWiflPTZEpe/GdhP0uGSOkvqnhaeAF2BNyLifUnDSHqCS5KOBX4M+ANJAT8LQNLWkr6cjkX+AHiHpFc4n6uAn0raNl22p6RR6ecNgZtIesuPB/pIOiWT9ypgKdBZ0rkkY41zfg/8RNJgJT6VjmOGpKd+UKn7aWaW5QLZzKyZIuJF4HzgUeAl4B+Flyi4rtXAwSTjm18GXicpADcvcflXSS5WOwN4g+SCuV3S8CnA+ZLeBs4F7ljH9G4B9iPTe0zy+3EGsCDd3j7pdvK5jKQH++E0hyeA3dLYz4G6iLgyHZ98FHCBpMHAQyQXHr4IvAK8T/1hJpek+/IwsILkAr+N0th5wA3psI7D13F/zazG+UEhZmZmZmYZ7kE2MzMzM8twgWxmZmZmluEC2czMzMwswwWymZmZmVlG5+JNql+PHj1iwIABlU7DzMzMzMpo8uTJr0dEz4bzXSADAwYMYNKkSZVOw8zMzMzKSNIr+eZ3qCEWkvqlj0udJWmmpNPztJGkcZJmS5ouaddK5GpmZmZmHVNH60FeBZwREVMkdQUmS3okIp7LtBkBDE5fuwFXsvaG9GZmZmZmBXWoHuSIWBgRU9LPbwOzgD4Nmo0CbozEE8AWknqXOVUzMzMz66A6VIGcJWkA8GngyQahPtR/FGkdjYtoJI2RNEnSpKVLl7ZVmmZm7ddFF8HWW8O111Y6EzOzdqVDFsiSNgXuBsZGxIqG4TyLNHqedkRcExFDI2Joz56NLl40M6tuEXD55bBkCfz2t5XOxsysXelwBbKk9UmK45sj4k95mtQB/TLTfYEF5cjNzKzDeOYZePVV+MQnYMoUeCXvhdxmZjWpQxXIkgRcC8yKiEuaaDYeOCa9m8XuwPKIWFi2JM3MOoJ774X11ls7vOK++yqbj5lZO9KhCmTgc8DRwBckTU1fIyWdJOmktM0EYC4wG/gdcEqFcjUza7/uuQf22gv22AN22imZNjMzoIPd5i0i/kH+McbZNgF8qzwZmZl1QCtXwrPPwvnnJ9P77Qe//30yLlkFv2LNzGpCR+tBNjOzlpo3L3n/2MeS9+22g3ffBd/Rx8wMcIFsZlZ75s5N3gcNqv+em29mVuNcIJuZ1ZqXX07eBw6s/+4C2cwMcIFsZlZ75s6FTTaB3D3gBwxI3nOFs5lZjXOBbGZWa+bOTYZV5C7I23hj6N3bPchmZikXyGZmtWbu3LXDKnIGDnSBbGaWcoFsZlZLItb2IGcNGuQC2cws5QLZzKyWLF2a3Ac5X4FcVwcffliZvMzM2hEXyGZmtaThLd5yBg2CNWvg1VfLn5OZWTvjAtnMrJbkCuR8Y5CzcTOzGuYC2cysluQK4Nyt3XL8sBAzs/9wgWxmVkvmzk1u6bbxxvXnb7MNdOniAtnMDBfIZma15eWXG48/BlhvvWSYhR8WYmbmAtnMrKbkuwdyzsCBMGdOefMxM2uHXCCbmdWKNWvgtdegX7/88f79k7iZWY1zgWxmViveeANWr4ZevfLHe/VK7pO8alV58zIza2dcIJuZ1YrFi5P3rbfOH+/VK3nS3tKl5cvJzKwd6lzOjUn6TgnN3o2Iq9s8GTOzWrNoUfJeqAcZYOHC5E4XZmY1qtw9yN8FNgW6Fnid0dTCkq6TtETSs03Eh0taLmlq+jq31ffAzKyjKtaDnCuKc4W0mVmNKmsPMvDHiDi/UANJmxQIXw/8BrixQJvHI+KgZuRmZlbdSu1BdoFsZjWurD3IEfG9lrSJiL8Db7RqUmZmtWLx4uRhIJtvnj+e61l2gWxmNa4iF+lJOl3SZkpcK2mKpC+20ur3kDRN0oOSdmqldZqZdXyLFiW9xFL++EYbJcXzwoXlzcvMrJ2p1F0svhERK4AvAj2B44FftMJ6pwDbRsQuwOXAvU01lDRG0iRJk5b6im0zqwWLFzc9/jind2/3IJtZzatUgZzrvhgJ/CEipmXmNVtErIiId9LPE4D1JfVoou01ETE0Iob27NmzpZs2M2v/Fi0qXiD36uUC2cxqXqUK5MmSHiYpkB+S1BVY09KVSuolJf93KGkYyf4ta+l6zcyqwuLFTV+gl+MC2cys7PdB7hwRq4BvAkOAuRGxUlJ3kmEWxZa/FRgO9JBUB/wYWB8gIq4CDgNOlrQKeA8YHRHRJjtjZtaRrF4NS5aU1oPsMchmVuPKfZu3J9LC9i/AXyLiLYCIWEYJPb0RcWSR+G9IbgNnZmZZy5bBmjXFe5B794Z334V33oFNNy1PbmZm7Uy5b/M2FDg9nbxU0tOSfi3pi5I2KGcuZmY1JTdsopQe5Gx7M7MaVPYxyBHxSkRcFRGHAHsCfwb2Ax6X9EC58zEzqwm5p+iVMgYZXCCbWU0r9xCLeiLiI+D/0heS+lQyHzOzqrWuPcgeh2xmNaxSDwo5SNIzkt6UtELS25JWRMRrlcjHzKzqldqD3Lt38u4eZDOrYZXqQb4UOBSY4btMmJmVwaJFsOGG0LVr4Xbdu0OnTi6QzaymVeo+yPOBZ10cm5mVSe4eyE09ZjpnvfWSYRgukM2shlWqB/l7wARJfwM+yM2MiEsqlI+ZWXUr5THTOX5YiJnVuEoVyD8F3gE2BLpUKAczs9qxaBEMGlRa2969YcGCts3HzKwdq1SBvGVEfLFC2zYzqz2LF8Mee5TWtlcvmDKlbfMxM2vHKjUG+VFJLpDNzMph1SpYurT4HSxyevVKHku9enXb5mVm1k5VqkD+FvAXSe9lb/NWoVzMzKrb669DROljkHv3TorjZcvaNi8zs3aqIkMsIqLIfYbMzKzV5C64W5ceZEgeFrLVVm2Tk5lZO1bWHmRJRb+dS2ljZmbrIPeQkHW5iwX4ThZmVrPKPcRiQiu1MTOzUjW3B9kFspnVqHIPsdilyFhjAR6LbGbWmtyDbGa2TspaIEdEp3Juz8zMSArdTTaBTTctrf0mmySPpF64sG3zMjNrpyp1FwszMyuXdXmKXo6fpmdmNcwFsplZtVu0qPTxxzkukM2shrlANjOrds3pQe7d2wWymdWsihTIki6StFMzlrtO0hJJzzYRl6RxkmZLmi5p15Zna2bWwTW3B9ljkM2sRlWqB/l54BpJT0o6SdLmJS53PXBAgfgIYHD6GgNc2aIszcw6uo8+Sp6I15wxyCtWwMqVbZOXmVk7VpECOSJ+HxGfA44BBgDTJd0i6b+KLPd34I0CTUYBN0biCWALSb1bK28zsw5n6dLkvTk9yLD2FnFmZjWkYmOQJXUCPp6+XgemAd+RdFsLVtsHmJ+Zrkvn5dv+GEmTJE1amvsBMTOrNrlxxM0Zg5xd3syshlRqDPIlJMMsRgI/i4jPRMQvI+Jg4NMtWXWeeZGvYURcExFDI2Joz549W7BJM7N2LDeOuLk9yB6HbGY1qNxP0st5FjgnIvINbhvWgvXWAf0y032BBS1Yn5lZx5YrcLfZZt2W89P0zKyGVWqIxdcbFseSJgJExPIWrHc8cEx6N4vdgeUR4e4PM6tdC9I+gnXtQe7ZE9ZbzwWymdWksvYgS9oQ2BjoIakba4dEbAYU7d6QdCswPF2+DvgxsD5ARFwFTCAZtjEbWAkc38q7YGbWsSxYkBS7Xbqs23KdOsFWW7lANrOaVO4hFicCY0mK4SmZ+SuAK4otHBFHFokH8K2WJGhmVlUWLlx7wd268r2QzaxGlbVAjojLgMskfTsiLi/nts3MatKCBes+/jjHj5s2sxpV7iEWX4iI/wNek3Row3hE/Kmc+ZiZVb2FC2HnnZu3bK9e8GzeB5eamVW1cg+x2Af4P+DgPLEAXCCbmbWW1auTHuDm9iD37p08KGTNmuSCPTOzGlHuIRY/Tt998ZyZWVtbujQpklsyxOKjj+DNN6F799bNzcysHavUg0J+JmmLzHQ3SRdUIhczs6qVu8CuJRfpZddjZlYjKvV/ZiMi4q3cRES8SXJ7NjMzay25eyC3pAcZfKGemdWcShXInSRtkJuQtBGwQYH2Zma2rlpaIOd6nl0gm1mNqdSjpm8CJkr6A8nFed8AbqhQLmZm1Sk3NGLrrZu3vHuQzaxGVaRAjogLJU0H9ktn/SQiHqpELmZmVau5T9HL2XRT2Hhjj0E2s5pTqR5kgGdIHhMd6WczM2tNLXmKHoDkh4WYWU2q1F0sDgeeAg4DDgeelHRYJXIxM6taLXmKXk7v3i6QzazmVKoH+YfAZyNiCYCknsCjwF0VysfMrPosWACf+lTL1tGrF8ya1Tr5mJl1EJW6i8V6ueI4tayCuZiZVZ/Vq5On4LVkiAUkBbLHIJtZjalUD/JfJD0E3JpOHwFMqFAuZmbVp6VP0cvp1St5kt4HH8AGvhunmdWGSt3F4ruSvgp8DhBwTUTcU4lczMyqUq7XtzXGIEPSG92/f8vWZWbWQVTsLhYRcTdwd6W2b2ZW1V59NXnv06dl68kt/+qrLpDNrGaUddyvpLclrcjzelvSinLmYmZW1ebMSd63265l6xk0KHl/+eWWrcfMrAMpaw9yRHQt5/bMzGrW3Lmw+ebQrVvL1rPttsn9kHMFt5lZDajYnSMkfV7S8ennHpIGlrjcAZJekDRb0ll54sMlLZc0NX2d29q5m5m1e3PmJL3HUsvWs8EG0K9fUnCbmdWIioxBlvRjYCiwA/AHoAtwE8lFe4WW6wRcAewP1AFPSxofEc81aPp4RBzU6ombmXUUc+e2/B7IOYMGuQfZzGpKpXqQvwJ8GXgXICIWAKUMvxgGzI6IuRHxIXAbMKrNsjQz64hWr07GDLd0/HHOdtu5QDazmlKpAvnDiAggACRtUuJyfYD5mem6dF5De0iaJulBSTvlW5GkMZImSZq0dOnSdcndzKx9e+01+OijtRfYtdSgQclt3t59t3XWZ2bWzlWqQL5D0tXAFpL+h+Qx078rYbl8g+miwfQUYNuI2AW4HLg334oi4pqIGBoRQ3v27LkOqZuZtXOtdQeLnNx6PA7ZzGpERQrkiLgIuIvkPsg7AOdGxOUlLFoH9MtM9wUWNFj3ioh4J/08AVhfUo9WSdzMrCPIFbKt2YOcXa+ZWZWr1EV6/w+4MyIeWcdFnwYGp3e8eA0YDXytwbp7AYsjIiQNI/lHwLJWSNvMrGOYMwc6d07uPtEacj3IHodsZjWiUk/S2wx4SNIbJBfa3RURi4stFBGrJJ0KPAR0Aq6LiJmSTkrjVwGHASdLWgW8B4xOxzubmdWGuXOT+xd3bqWv+G7dknsquwfZzGpERQrkiPhf4H8lfQo4AvibpLqI2K+EZScAExrMuyrz+TfAb1o5ZTOzjmPOnNYbXgHJvZR9qzczqyEVe1BIagmwiGQIxFYVzsXMrDrMndt6F+jlbLede5DNrGZUpECWdLKkx4CJQA/gfyKile5ob2ZWw956C954o3V7kCFZ38svJ/dYNjOrcpUag7wtMDYiplZo+2Zm1Wn27OS9LXqQP/oI5s+HAQNad91mZu1MpW7zdpaLYzOzNjBlSvK+yy6tu97c+nLrNzOrYpUeg2xmZq3pqadgyy1bf4jFLrvA+usn6zczq3IukM3MqslTT8GwYcmdJ1rThhsmRbILZDOrAS6QzcyqxbvvwsyZSYHcFoYNg0mTYM2atlm/mVk74QLZzKxaTJmSFK9tWSC//Ta88ELbrN/MrJ1wgWxmVi1ywx8++9m2WX+u8PYwCzOrci6QzcyqxVNPJbdg26qNnru0ww7QtasLZDOrei6QzcyqRe4Cvbay3npJ77QLZDOrci6QzcyqwZIlMG9e2w2vyPnsZ2HaNHj//bbdjplZBblANjOrBg88kLzvs0/bbmf48OSJeo880rbbMTOrIBfIZmbV4LbbkoeDDB3attvZd9/kQSS33da22zEzqyAXyGZmHd2SJTBxIowe3foPCGlo/fXhsMNg/HhYubJtt2VmViEukM3MOro774TVq+GII8qzvSOOgHfegT//uTzbMzMrMxfIZmYd2apV8Otfw2c+AzvvXJ5t7rNPMpzjoosgojzbNDMrIxfIZmYd2S9/CXPmwLnntv3wipxOneCcc5LHTv/ud+XZpplZGXW4AlnSAZJekDRb0ll54pI0Lo1Pl7RrJfI0M2tzjz8O552XjD0++ODybvvYY2H//eH//T+YPr282zYza2MdqkCW1Am4AhgB7AgcKWnHBs1GAIPT1xjgyrImaWbW1lavhgsvhC98Afr3h9/+tny9xznrrQfXXw+bbQa77w5XXeXhFmZWNTpXOoF1NAyYHRFzASTdBowCnsu0GQXcGBEBPCFpC0m9I2Jh+dNtwqJFDP/iFxvNPnyvvTjloINY+f77jPzxjxvFj9t3X47bf39eX76cw37+80bxk0eO5Ii992b+0qUcffHFyczMD9YZX/kKB++2Gy/U1XHib37TaPlzRo9mvyFDmDp3LmOvuaZR/GfHHMOen/gE/5o1i7NvvLFR/NITTmDIoEE8OnUqF9xxR6P41aecwg59+/Lnp57i4nvvbRT/49ix9OvZk9sff5wr//KXRvG7vv99emy2GddPnMj1Eyc2ik8491w23mADfjthAnf885+N4o/99KcQwUX33sv9kybVi23UpQsPnnsuAD+54w4mNugR6961K3d/73sA/OCmm/j3Cy/Ui/ft3p2bxo4FYOy11zJ13rx68e179+aak08GYMyVV/Liwvqn45ABA7j0G9+ACI667DLqli2rF99jhx34+de/DsBXf/Urlr39dr34vjvvzI8OOwyAET/9Ke99+GG9+EGf+QxnfvnLAAzPc24dvscenPKlL7Hygw8YmefcOm74cI4bPjw5937960bxk/ffnyP23JP5r7/O0Vdc0Sh+xoEHcvBnPsMLCxZwYp7/kj/n0EPZb+edmTpvHmNvuKFR/GejR7Pn9tvzrxdf5Ow8txe79JhjGDJgAI/OmMEF99zTKH71N7/JDttsw58nT+biCRMaxf948sn0696d2//9b67Mc27ddfrp9Ojalev/9jeuf/zxtYH072vCd7+bnHuPPMIdTz7ZaPnHzjkHgIseeID7n3mmXmyjLl148LvfBeAn997LxJkz68W7b7opd59+Orz/Pj+48cbk3HvjjWTscY8e9N11V27q1g2AsWPHMnXq1HrLb7/99lyT/j2PGTOGF198sV58yJAhXHrppQAcddRR1NXV1Yvvscce/Dw9J7761a+yrMG5ue/RR/OjadPg5JMZcfrpvNetG2y8cXK3i06dknPvkEOgc2eGn3NOo0K+6Pfefvut/d772c8axU8eOZIj9tkn+d676KJG8TMOPXTt997llzeKnzN6NPt9+tNMnTMn//fescey54478q/nnuPsPOfmpWPGMGS77Xj0mWe4IM+5efW3v5187z35JBf/6U+N4n8888zke+9vf+PKPOfmXWefTY/NN+f6Rx7h+kcfbRSf8L//y8Ybbshv77+fO7LnZuqxX/4SgIvuvpv7GzwBcaMuXXjwJz8B4Ce33MLEadPqxbt37crd6bn7gz/8gX8//3y9eN8ePbgpPXfHXn01U+fOrRffvk8frjntNADGjBvHi6+9Vi8+ZNAgLj3xRACO+tWvqHv99XrxPT7xCX5+/PEAfPWCC1i2YkW9+L5DhvCjr30NgBE/+hHvffBBvfhBw4ZxZvq9ODz9/s46fO+915576fd/1nH777/23PvpTxvFTz7wwLXn3q9+1Sh+xqGHcvDuuyfn3rhxjeLnHHnk2nPv6qsbxX923HFrz73rr28Uv/TEE9eee7fe2ih+9WmnJefeE0/kP/e++921517uPuoZd/3wh2vPvTz3PZ9w/vlrz72//71R/LELLwTgorvuanzubbBB/XNv6tSk/bBhyW0k24GOViD3AeZnpuuA3Upo0weoV5FIGkPSw0z//v1bPdGC/vUvmDGj8fwZM5KeoKbMmAHpD1mT8fTLMK9nny2cV/pF2KQ8XzD1pAVik045pXD8hBMKx48+unC82BX8hxxSOP6VrxSOH3po4XieL4j/mDkT8vy41Yvn+YL6j+eegzyFX7347bcXjv/xj4Xj115bOF7o3HzuObjsssLxQvL8+NSTp3Cq5wc/KBw/88zC8fRHvEknnVQ4/s1vFo4fdVTheLFzOxvv0gW6d4cePZLXBhsUXratde0KDz6YnH/f+Q4sWwaLF6+Nz5oFN93U9PKlfO/l+UdZvXih771837XrEi927nz724XjY8YUjh97bOH4kUcWjhf7Xho5sm3jf/1r07EZMyBPZ0e9+H33FY7n6WypFy/0vTZjRuHvtRkzIM8/6OvFL7mkcPwXvygcL6RY/IwzCsdPPbVw/H/+p3D8mGMKx0ePLhwv9ps5YsS6xUeMgMcea/uHHZVI0YH+S0zSfwNfiogT0umjgWER8e1MmweAn0fEP9LpicD3ImJyU+sdOnRoTGrQo9im3ngDXnqpecs2579Rm/tfr95W+ZfxtrytpnTpAltvXfmCuBTvvZcUyR98kAwHWbVq7asD/eZYBfj8qG077pgM2yojSZMjotETljpaD3Id0C8z3RdY0Iw2lbXllrBbw45vM7MqsdFGMGBApbMwM2u2DnWRHvA0MFjSQEldgNHA+AZtxgPHpHez2B1Y3q7GH5uZmZlZu9ahepAjYpWkU4GHgE7AdRExU9JJafwqYAIwEpgNrASOr1S+ZmZmZtbxdKgCGSAiJpAUwdl5V2U+B/CtcudlZmZmZtWhQ12k11YkLQVeKfNmewCvF21V23yMivMxKs7HqDAfn+J8jIrzMSrOx6i4ShyjbSOiZ8OZLpArRNKkfFdN2lo+RsX5GBXnY1SYj09xPkbF+RgV52NUXHs6Rh3tIj0zMzMzszblAtnMzMzMLMMFcuU0fqapNeRjVJyPUXE+RoX5+BTnY1Scj1FxPkbFtZtj5DHIZmZmZmYZ7kE2MzMzM8twgWxmZmZmluECuY1JOkDSC5JmSzorT1ySxqXx6ZJ2rUSelSKpn6S/Spolaaak0/O0GS5puaSp6evcSuRaSZLmSZqR7v+kPPGaPY8k7ZA5N6ZKWiFpbIM2NXcOSbpO0hJJz2bmbSnpEUkvpe/dmli24PdWtWjiGP1K0vPp39E9krZoYtmCf5PVooljdJ6k1zJ/TyObWLaWz6PbM8dnnqSpTSxb9edRU7/z7f77KCL8aqMXyeOw5wCDgC7ANGDHBm1GAg8CAnYHnqx03mU+Rr2BXdPPXYEX8xyj4cD9lc61wsdpHtCjQLymz6PMcegELCK58Xt2fs2dQ8DewK7As5l5FwJnpZ/PAn7ZxDEs+L1VLa8mjtEXgc7p51/mO0ZprODfZLW8mjhG5wFnFlmups+jBvGLgXNr9Txq6ne+vX8fuQe5bQ0DZkfE3Ij4ELgNGNWgzSjgxkg8AWwhqXe5E62UiFgYEVPSz28Ds4A+lc2qQ6rp8yhjX2BORJT7yZjtTkT8HXijwexRwA3p5xuAQ/IsWsr3VlXId4wi4uGIWJVOPgH0LXti7UgT51Epavo8ypEk4HDg1rIm1Y4U+J1v199HLpDbVhtFQJ8AACAASURBVB9gfma6jsbFXyltaoKkAcCngSfzhPeQNE3Sg5J2Kmti7UMAD0uaLGlMnrjPo8Romv4hqvVzCGDriFgIyY8WsFWeNj6X1voGyf/M5FPsb7LanZoOQ7muif8a93mU2AtYHBEvNRGvqfOowe98u/4+coHctpRnXsP76pXSpupJ2hS4GxgbESsahKeQ/Jf5LsDlwL3lzq8d+FxE7AqMAL4lae8G8Zo/jyR1Ab4M3Jkn7HOodDV/LgFI+iGwCri5iSbF/iar2ZXAdsAQYCHJEIKGfB4ljqRw73HNnEdFfuebXCzPvLKcRy6Q21Yd0C8z3RdY0Iw2VU3S+iR/NDdHxJ8axiNiRUS8k36eAKwvqUeZ06yoiFiQvi8B7iH5b6esmj+PSH5gpkTE4oYBn0P/sTg39CZ9X5KnTc2fS5KOBQ4Cvh7pQMiGSvibrFoRsTgiVkfEGuB35N93n0dSZ+BQ4Pam2tTKedTE73y7/j5ygdy2ngYGSxqY9m6NBsY3aDMeOCa9C8HuwPLcfznUgnR81rXArIi4pIk2vdJ2SBpGct4uK1+WlSVpE0ldc59JLiJ6tkGzmj6PUk321NT6OZQxHjg2/XwscF+eNqV8b1UtSQcA3we+HBErm2hTyt9k1WpwfcNXyL/vNX0epfYDno+IunzBWjmPCvzOt+/vo3JdxVirL5K7C7xIchXmD9N5JwEnpZ8FXJHGZwBDK51zmY/P50n+u2Q6MDV9jWxwjE4FZpJcvfoEsGel8y7zMRqU7vu09Dj4PGp8jDYmKXg3z8yr6XOI5B8LC4GPSHphvgl0ByYCL6XvW6ZttwEmZJZt9L1Vja8mjtFskjGPue+jqxoeo6b+Jqvx1cQx+mP6PTOdpFjp7fOo/jFK51+f+w7KtK2586jA73y7/j7yo6bNzMzMzDI8xMLMzMzMLMMFspmZmZlZhgtkMzMzM7MMF8hmZmZmZhkukM3MzMzMMlwgm5mZmZlluEA2M6tSkoZLWi5pQmb6/has72ZJb0g6rPWyNDNrf1wgm5lVt8cjYmRrrCgivk7tPQ3NzGqQC2Qzsyog6bOSpkvaMH2E7Uzgk0XaPyNpkKTzJN0g6WFJ8yQdKulCSTMk/UXS+mXcFTOzinOBbGZWBSLiaZLe3QuAC4GbgGfztZW0J3AVMCoi5qaztwMOBEaly/41InYG3kvnm5nVjM6VTsDMzFrN+cDTwPvAacBeedp8ArgG+GJELMjMfzAiPpI0A+gE/CWdPwMY0GYZm5m1Q+5BNjOrHlsCmwJdgQ2baLOQpID+dIP5HwBExBrgo4iIdP4a3JliZjXGBbKZWfW4BvgRcDPwyybavEUyZOJnkoaXKS8zsw7FvQJmZlVA0jHAqoi4RVIn4F/An/K1jYjFkg4GHpT0jXLmaWbWEWjt/6KZmVk1SXuIz4yIg1pxndcD90fEXa21TjOz9sZDLMzMqteHwCdzDwppKUk3A/uQjGE2M6ta7kE2MzMzM8twD7KZmZmZWYYLZDMzMzOzDBfIZmZmZmYZLpDNzMzMzDJcIJuZmZmZZbhANjMzMzPLcIFsZmZmZpbhAtnMzMzMLMMFspmZmZlZhgtkM7MOStLZkn7fwnUMkBSSOrdWXmZmHZ0fNW1mVsMkDQBeBtaPiFWVzab1SApgcETMrnQuZtbxuAfZzKwDco+vmVnbcYFsZlaApHmSzpQ0XdJySbdL2jCNHSfpHw3ah6SPpZ+vl/RbSQ9KekfSPyX1knSppDclPS/p05llt5F0t6Slkl6WdFomdp6kuyTdJGkFcFw676ZMm89L+pektyTNl3RcOv9ASc9IWpHOP6/EfT9L0l0N5l0maVxm/+dKejvN9+tNrGe9dF1zJC2TdIekLdPYEek6NkunR0haJKlnZnvz09wnS9ors95O6TCTOWkOkyX1k/T3tMm09LgfUcr+mpnluEA2MyvucOAAYCDwKeC4dVz2HKAH8AHwb2BKOn0XcAkkRSTwZ2Aa0AfYFxgr6UuZdY1Kl9kCuDm7EUn9gQeBy4GewBBgahp+FzgmXe5A4GRJh5SQ+63AyEzx2indn1skbQKMA0ZERFdgz8z2GjoNOATYB9gGeBO4AiAibk+PyThJ3YFrgRMiYmm67NPpvmwJ3ALcmfsHCvAd4EhgJLAZ8A1gZUTsncZ3iYhN022YmZXMBbKZWXHjImJBRLxBUsQOWYdl74mIyRHxPnAP8H5E3BgRq4HbgVwP8meBnhFxfkR8GBFzgd8BozPr+ndE3BsRayLivQbb+TrwaETcGhEfRcSyiJgKEBGPRcSMdLnpJIXvPsUSj4hXSIr5XDH9BZIC9Il0eg3wSUkbRcTCiJjZxKpOBH4YEXUR8QFwHnBYZpjIt9J1Pwb8OSLuz+RwU7ovqyLiYmADYIc0fAJwTkS8EIlpEbGs2H6ZmRXjAtnMrLhFmc8rgU3XYdnFmc/v5ZnOrWtbYJt0eMRbkt4Czga2zrSfX2A7/YA5+QKSdpP013ToxnLgJJIe7FLcQtJLC/C1dJqIeBc4Il3XQkkPSPp4E+vYFrgns1+zgNWk+xYRbwF3Ap8ELm6Q+xmSZqXDW94CNs/k3uQ+m5m1hAtkM7PmexfYODchqVcL1jUfeDkitsi8ukbEyEybQrcdmg9s10TsFmA80C8iNgeuAlRiXncCwyX1Bb6SritJJuKhiNgf6A08T9Lj3VRuIxrs24YR8RqApCEkwyNuJRm2QTp/L+D7JMM6ukXEFsDyTO6F9tnMrNlcIJuZNd80YCdJQ9Jxsee1YF1PASskfV/SRukFaJ+U9NkSl78Z2E/S4ZI6S+qeFp4AXYE3IuJ9ScNIeoJLko4Ffgz4A0kBPwtA0taSvpyORf4AeIekVzifq4CfSto2XbanpFHp5w2Bm0h6y48H+kg6JZP3KmAp0FnSuSRjjXN+D/xE0mAlPpWOY4akp35QqftpZpblAtnMrJki4kXgfOBR4CXgH4WXKLiu1cDBJOObXwZeJykANy9x+VdJLlY7A3iD5IK5XdLwKcD5kt4GzgXuWMf0bgH2I9N7TPL7cQawIN3ePul28rmMpAf74TSHJ4Dd0tjPgbqIuDIdn3wUcIGkwcBDJBcevgi8ArxP/WEml6T78jCwguQCv43S2HnADemwjsPXcX/NrMb5QSFmZmZmZhnuQTYzMzMzy3CBbGZmZmaW4QLZzMzMzCzDBbKZmZmZWUbn4k2qX48ePWLAgAGVTsPMzMzMymjy5MmvR0TPhvNdIAMDBgxg0qRJlU7DzMzMzMpI0iv55ld0iIWkAyS9IGm2pLPyxCVpXBqfLmnXTGyepBmSpkqalJm/paRHJL2Uvncr1/6YmZmZWcdXsQJZUifgCmAEsCNwpKQdGzQbAQxOX2OAKxvE/ysihkTE0My8s4CJETEYmJhOm5mZmZmVpJI9yMOA2RExNyI+BG4DRjVoMwq4MRJPAFtI6l1kvaOAG9LPNwCHtGbSZmZmZlbdKlkg96H+I0Pr0nmltgmSx5ZOljQm02briFgIkL5vlW/jksZImiRp0tKlS1uwG2ZmZmZWTSpZICvPvIbPvS7U5nMRsSvJMIxvSdp7XTYeEddExNCIGNqzZ6OLF83MzMysRlWyQK4D+mWm+wILSm0TEbn3JcA9JEM2ABbnhmGk70taPXMzMzMzq1qVLJCfBgZLGiipCzAaGN+gzXjgmPRuFrsDyyNioaRNJHUFkLQJ8EXg2cwyx6afjwXua+sdMTMzM7PqUbH7IEfEKkmnAg8BnYDrImKmpJPS+FXABGAkMBtYCRyfLr41cI8kSPbhloj4Sxr7BXCHpG8CrwL/XaZdMjMzM7MqoIiGw35rz9ChQ8MPCjEzMzOrLZImN7hdMFDhB4WYmZmZmbU3LpDNzMzMzDJcIJuZmZmZZbhANjMzMzPLcIFsZmZmZpbhAtnMzMzMLMMFspmZmZlZhgtkMzMzM7MMF8hmZmZmZhkukM3MzMzMMlwgm5mZmZlluEA2MzMzM8twgWxmZmZmluEC2czMzMwswwWymZmZmVmGC2QzMzMzswwXyGZmZmZmGS6QzczMzMwyKlogSzpA0guSZks6K09cksal8emSdk3n95P0V0mzJM2UdHpmmfMkvSZpavoaWc59MjMzM7OOrXOlNiypE3AFsD9QBzwtaXxEPJdpNgIYnL52A65M31cBZ0TEFEldgcmSHsks++uIuKhc+2JmZmZm1aNggSxpegnrWBoR+zZj28OA2RExN93WbcAoIFsgjwJujIgAnpC0haTeEbEQWAgQEW9LmgX0abCsmZmZmdk6K9aD3AkoNERBwPhmbrsPMD8zXUfSO1ysTR/S4hhA0gDg08CTmXanSjoGmETS0/xmo8SlMcAYgP79+zdzF8zMzMys2hQbg3xiRLxS4DUPOKWZ21aeebEubSRtCtwNjI2IFensK4HtgCEkhfTF+TYeEddExNCIGNqzZ891zd3MzMzMqlTBAjki/lFsBaW0aUId0C8z3RdYUGobSeuTFMc3R8SfMvksjojVEbEG+B3JUA4zMzMzs5KUdBcLSYMl3SXpOUlzc68WbvtpYLCkgZK6AKNpPFxjPHBMejeL3YHlEbFQkoBrgVkRcUmDXHtnJr8CPNvCPM3MzMyshpR6F4s/AD8Gfg38F3A8+Yc/lCwiVkk6FXiIZKzzdRExU9JJafwqYALJGOjZwMp0uwCfA44GZkiams47OyImABdKGkIyFGMecGJL8jQzMzOz2qLkBhFFGkmTI+IzkmZExM7pvMcjYq82z7AMhg4dGpMmTap0GmZmZmZWRmmNO7Th/FJ7kN+XtB7wUtrr+xqwVWsmaGZmZmbWHpT6JL2xwMbAacBngKOAY9sqKTMzMzOzSinagyypJ8l43s4RUcfaccBmZmZmZlWnYA+ypBOAmcDlwPOSvlyWrMzMzMzMKqRYD/JYYKeIWCppEHAzzX9ynpmZmZlZu1dsDPKHEbEUICLmAhu0fUpmZmZmZpVTrAe5r6RxTU1HxGltk5aZmZmZWWUUK5C/22B6clslYmZmZmbWHhQskCPihnIlYmZmZmbWHhQskCUVvCAvInxXCzMzMzOrKsWGWOwBzAduBZ4E1OYZmZmZmZlVULECuRewP3Ak8DXgAeDWiJjZ1omZmZmZmVVCwdu8RcTqiPhLRBwL7A7MBh6T9O2yZGdmZmZmVmalPGp6A+BAkl7kAcA44E9tm5aZmZmZWWUUu0jvBuCTwIPA/0bEs2XJyszMzMysQor1IB8NvAtsD5wm/ecaPQEREZu1YW5mZmZmZmVX7D7IxR5FbWZmZmZWVVwAm5mZmZllFCyQJU0ptoJS2hRY9gBJL0iaLemsPHFJGpfGp0vatdiykraU9Iikl9L3bs3Nz8zMzMxqT7ExyJ+QNL1AXMDmzdmwpE7AFST3Wa4DnpY0PiKeyzQbAQxOX7sBVwK7FVn2LGBiRPwiLZzPAr7fnBzNzMzMrPYUK5A/XsI6Vjdz28OA2RExF0DSbcAoIFsgjwJujIgAnpC0haTeJLeba2rZUcDwdPkbgMdobwXy2LEMv/76RrMP79mTU/r0YeXq1YycMaNR/LhevTiuVy9e/+gjDpvZ+FktJ2+zDUdstRXz33+fo59/vlH8jL59ObhHD15YuZITX3yxUfycbbdlv27dmPrOO4ydPbtR/GcDB7Ln5pvzr+XLOfvllxvFL/3Yxxiy6aY8+uabXPDKK43iV2+/PTtsvDF/fv11Lq6raxT/48c/Tr8NN+T2JUu4csGCRvG7dtqJHuuvz/WLFnH9okWN4hN23pmNO3Xit6+9xh1LlzaKPzZkCAAXzZ/P/cuW1YtttN56PPipTwHwk1deYeKbb9aLd19/fe7eaScAfjB3Lv9esaJevO8GG3DTJz4BwNjZs5n6zjv14ttvvDHXbL89AGNefJEXV66sFx+y6aZc+rGPAXDUrFnUffBBvfgem23GzwcNAuCrM2ey7KOP6sX37daNH227LQAjpk/nvTVr6sUP6t6dM/v1A2D41Kk05HPP5x743PO553Mvy+deec+9x447Di69tNF6KqXYRXqNj3jr6UPyGOucOpJe4mJt+hRZduuIWAgQEQslbZVv45LGAGMA+vfv38xdMDMzM7Nqo6RztgIblv4b+FJEnJBOHw0Mi4hvZ9o8APw8Iv6RTk8EvgcMampZSW9FxBaZdbwZEQXHIQ8dOjQmTZrUyntoZmZmZu2ZpMkRMbTh/ErexaIO6JeZ7gs07ONvqk2hZRenwzBI35e0Ys5mZmZmVuVKKpAlbSJpvfTz9pK+LGn9Fm77aWCwpIGSugCjgfEN2owHjknvZrE7sDwdPlFo2fHAsennY4H7WpinmZmZmdWQYhfp5fwd2Cu9ZdpEYBJwBPD15m44IlZJOhV4COgEXBcRMyWdlMavAiYAI4HZwErg+ELLpqv+BXCHpG8CrwL/3dwczczMzKz2lDQGWdKUiNhV0reBjSLiQknPRMSn2z7FtucxyGZmZma1p6VjkCVpD5Ie4wfSeaX2PpuZmZmZdRilFshjgR8A96TDIAYBf227tMzMzMzMKqOkXuCI+BvwN4D0Yr3XI+K0tkzMzMzMzKwSSr2LxS2SNpO0CcnT6l6Q9N22Tc3MzMzMrPxKHWKxY0SsAA4hubNEf+DoNsvKzMzMzKxCSi2Q10/ve3wIcF9EfARU5hF8ZmZmZmZtqNQC+WpgHrAJ8HdJ2wIr2iopMzMzM7NKKfUivXHAuMysVyT9V9ukZGZmZmZWOaVepLe1pGslPZhO78jaxzmbmZmZmVWNUodYXE/yWOdt0ukXSe6NbGZmZmZWVUotkHtExB3AGoCIWAWsbrOszMzMzMwqpNQC+V1J3UnvXCFpd2B5m2VlZmZmZlYhJV2kB3wHGA9sJ+mfQE/gsDbLyszMzMysQooWyOmjpTcE9gF2AAS8kN4L2czMzMysqhQtkCNijaSLI2IPYGYZcjIzMzMzq5hSxyA/LOmrktSm2ZiZmZmZVdi6jEHeBFgl6X2SYRYREZu1WWZmZmZmZhVQ6pP0urZ1ImZmZmZm7UGpT9KbWMq8UknaUtIjkl5K37s10e4ASS9Imi3prMz8X0l6XtJ0SfdI2iKdP0DSe5Kmpq+rmpujmZmZmdWmggWypA0lbQn0kNQtLWy3lDSAtU/Va46zgIkRMRiYmE433HYn4ApgBLAjcGT6iGuAR4BPRsSnSJ7q94PMonMiYkj6OqkFOZqZmZlZDSrWg3wiMBn4ePqee91HUrw21yjghvTzDcAhedoMA2ZHxNyI+BC4LV2OiHg4fZofwBNA3xbkYmZmZmb2HwUL5Ii4LCIGAmdGxKCIGJi+domI37Rgu1tHxMJ0GwuBrfK06QPMz0zXpfMa+gbwYGZ6oKRnJP1N0l4tyNHMzMzMalCpd7FYJKlrRLwt6RxgV+CCiJjS1AKSHgV65Qn9sMRt5rulXDTYxg+BVcDN6ayFQP+IWCbpM8C9knaKiBV58hsDjAHo379/iSmZmZmZWbUrtUD+UUTcKenzwJeAi4Argd2aWiAi9msqJmmxpN4RsVBSb2BJnmZ1QL/MdF9gQWYdxwIHAftGRKTb/AD4IP08WdIcYHtgUp78rgGuARg6dGg0jJuZmZlZbSr1QSGr0/cDgSsj4j6gSwu2Ox44Nv18LMmY5oaeBgZLGiipCzA6XQ5JBwDfB74cEStzC0jqmV7ch6RBwGBgbgvyNDMzM7MaU2qB/Jqkq4HDgQmSNliHZfP5BbC/pJeA/dNpJG0jaQJAehHeqcBDwCzgjojIPer6N0BX4JEGt3PbG5guaRpwF3BSRLzRgjzNzMzMrMYoHZ1QuJG0MXAAMCMiXkqHRewcEQ+3dYLlMHTo0Jg0qdEoDDMzMzOrYpImR8TQhvNLfZLeSuBPmemFJBfEmZmZmZlVlZYMkzAzMzMzqzoukM3MzMzMMlwgm5mZmZlluEA2MzMzM8twgWxmZmZmluEC2czMzMwswwWymZmZmVmGC2QzMzMzswwXyGZmZmZmGS6QzczMzMwyXCCbmZmZmWW4QDYzMzMzy3CBbGZmZmaW4QLZzMzMzCzDBbKZmZmZWYYLZDMzMzOzDBfIZmZmZmYZLpDNzMzMzDIqUiBL2lLSI5JeSt+7NdHuAEkvSJot6azM/PMkvSZpavoamYn9IG3/gqQvlWN/zMzMzKx6VKoH+SxgYkQMBiam0/VI6gRcAYwAdgSOlLRjpsmvI2JI+pqQLrMjMBrYCTgA+G26HjMzMzOzklSqQB4F3JB+vgE4JE+bYcDsiJgbER8Ct6XLFVvvbRHxQUS8DMxO12NmZmZmVpJKFchbR8RCgPR9qzxt+gDzM9N16bycUyVNl3RdZohGsWX+Q9IYSZMkTVq6dGlz98PMzMzMqkybFciSHpX0bJ5XsV7g/6wiz7xI368EtgOGAAuBi0tYpv7MiGsiYmhEDO3Zs2eJKZmZmZlZtevcViuOiP2aiklaLKl3RCyU1BtYkqdZHdAvM90XWJCue3FmXb8D7i+2jJmZmZlZKSo1xGI8cGz6+VjgvjxtngYGSxooqQvJxXfjAdKiOucrwLOZ9Y6WtIGkgcBg4Kk2yN/MzMzMqlSb9SAX8QvgDknfBF4F/htA0jbA7yNiZESsknQq8BDQCbguImamy18oaQjJ8Il5wIkAETFT0h3Ac8Aq4FsRsbqM+2VmZmZmHZwi8g7RrSmSlgKvlHmzPYDXy7zNjsbHqDgfo+J8jArz8SnOx6g4H6PifIyKq8Qx2jYiGl2M5gK5QiRNioihlc6jPfMxKs7HqDgfo8J8fIrzMSrOx6g4H6Pi2tMx8qOmzczMzMwyXCCbmZmZmWW4QK6cayqdQAfgY1Scj1FxPkaF+fgU52NUnI9RcT5GxbWbY+QxyGZmZmZmGe5BNjMzMzPLcIFsZmZmZpbhArmNSTpA0guSZks6K09cksal8emSdq1EnpUiqZ+kv0qaJWmmpNPztBkuabmkqenr3ErkWkmS5kmake7/pDzxmj2PJO2QOTemSlohaWyDNjV3Dkm6TtISSc9m5m0p6RFJL6Xv3ZpYtuD3VrVo4hj9StLz6d/RPZK2aGLZgn+T1aKJY3SepNcyf08jm1i2ls+j2zPHZ56kqU0sW/XnUVO/8+3++ygi/GqjF8kTAOcAg4AuwDRgxwZtRgIPAgJ2B56sdN5lPka9gV3Tz12BF/Mco+HA/ZXOtcLHaR7Qo0C8ps+jzHHoBCwiufF7dn7NnUPA3sCuwLOZeRcCZ6WfzwJ+2cQxLPi9VS2vJo7RF4HO6edf5jtGaazg32S1vJo4RucB/7+9Ow+zpK7vPf7+OIgbCBIGIzOISHBBL6C2SIxJSNAIuIyJy2VRcImEJ5BoFANqTIxbNBpjvKJzJ8hFI0rigo5mFE3UaEQSBsPiiCQjqDMMwqCyCCgM871/VDUUh9ML093n9Ol+v56nnzlV9auq76nf6e7P/M7vVJ80xX6L+nXUs/1vgD9frK+jiX7Pz/efR44gz60DgfVVdXlV3QqcBazoabMC+HA1zgN2TvKQQRc6LFV1VVV9q318I3ApsGy4VY2kRf066jgE+F5VDfovY847VfU14Cc9q1cAH2offwh4Tp9dp/Nza0Hod42q6otVtaVdPA9YPvDC5pEJXkfTsahfR+OSBHgB8LGBFjWPTPJ7fl7/PDIgz61lwIbO8kbuHv6m02ZRSPIw4HHAf/TZ/KtJLkry+SSPGWhh80MBX0xyQZLj+mz3ddQ4gol/ES321xDAg6vqKmh+aQG79Wnja+lOL6V5Z6afqb4nF7oT22kop0/w1rivo8avA1dX1f9MsH1RvY56fs/P659HBuS5lT7reu+rN502C16SHYBPAq+sqht6Nn+L5i3z/YH/A3x60PXNA79WVY8HDgNOSPIbPdsX/esoyfbAs4GP99nsa2j6Fv1rCSDJ64EtwJkTNJnqe3Ih+wCwN3AAcBXNFIJevo4aRzL56PGieR1N8Xt+wt36rBvI68iAPLc2Ant0lpcDm7ahzYKW5N403zRnVtWnerdX1Q1V9bP28Rrg3kl2HXCZQ1VVm9p/rwHOpnnbqWvRv45ofsF8q6qu7t3ga+gOV49PvWn/vaZPm0X/WkpyLPBM4OhqJ0L2msb35IJVVVdX1e1VtRX4e/o/d19HyXbA7wH/OFGbxfI6muD3/Lz+eWRAnlvnA/sk2asd3ToCWN3TZjVwTHsXgoOA68ffclgM2vlZHwQurap3T9Dml9t2JDmQ5nX748FVOVxJHpBkx/HHNB8i+nZPs0X9OmpNOFKz2F9DHauBY9vHxwKf6dNmOj+3FqwkhwInA8+uqpsnaDOd78kFq+fzDb9L/+e+qF9HracC362qjf02LpbX0SS/5+f3z6NBfYpxsX7R3F3gv2k+hfn6dt3xwPHt4wCnttsvAcaGXfOAr89TaN4uuRi4sP06vOcanQiso/n06nnAk4dd94Cv0cPb535Rex18Hd39Gt2fJvDu1Fm3qF9DNP9ZuAq4jWYU5mXALwH/CvxP++8ubdvdgTWdfe/2c2shfk1wjdbTzHkc/3m0svcaTfQ9uRC/JrhG/9D+nLmYJqw8xNfRXa9Ru/6M8Z9BnbaL7nU0ye/5ef3zyD81LUmSJHU4xUKSJEnqMCBLkiRJHQZkSZIkqcOALEmSJHUYkCVJkqQOA7IkSZLUYUCWpAUqycFJrk+yprP8uRkc78wkP0nyvNmrUpLmHwOyJC1sX6+qw2fjQFV1NIvvr6FJWoQMyJK0ACR5YpKLk9y3/RO264DHTtH+v5I8PMkbk3woyReTfD/J7yX56ySXJPlCknsP8KlI0tAZkCVpAaiq82lGd98C/DXwEeDb/domeTKwElhRVZe3q/cGngGsaPf9SlX9L+CWdr0kLRrbDbsASdKseRNwPvBz4I+BX+/T5tHAKuB3qmpTZ/3nq+q2JJcAS4AvtOsvAR42ZxVLiQWnAgAAIABJREFU0jzkCLIkLRy7ADsAOwL3naDNVTQB+nE9638BUFVbgduqqtr1W3EwRdIiY0CWpIVjFfAG4EzgHRO0uY5mysTbkhw8oLokaaQ4KiBJC0CSY4AtVfXRJEuAc4FP9WtbVVcneRbw+SQvHWSdkjQKcue7aJKkhaQdIT6pqp45i8c8A/hcVX1ito4pSfONUywkaeG6FXjs+B8KmakkZwK/STOHWZIWLEeQJUmSpA5HkCVJkqQOA7IkSZLUYUCWJEmSOgzIkiRJUocBWZIkSeowIEuSJEkdBmRJkiSpw4AsSZIkdRiQJUmSpA4DsqSRkuSNST4yg/3XJTl4FkuaN5KsTPKGYdcxkele+ySV5FcGUFL3nN9P8tQZHuPzSY6drZokDc92wy5AkuZKkjOAjVX1Z+Prquoxw6toblXV8cOuYdxiu/YAVXXY+OMkLwZ+v6qeMryKJG0rR5AlzRtJFvx/2tO411TrpnGceXutkiwZdg2SNBMGZElD1b61fXKSi4GbkmyXZPckn0yyOckVSf54kv0/nuRHSa5P8rUkj2nXHwccDfxpkp8l+WznfE9tz3FLkl06x3pckmuT3LtdfmmSS5P8NMk5SfacpI6Dkpyb5LokF3WnEiT5apK3JvkGcDPw8AnW7Z5kdZKfJFmf5OWdY7wxySeSfCTJDcCL+9RwRpK3tI8PTrIxyauTXJPkqiQvmaT+ryb5qyT/2V7Lz/Rcm77XuXPeDyRZk+Qm4GWTXfv28ZIkr0vyvSQ3JrkgyR596rpPkncl+WGSq9tpJPeboN11SR7bWbe07ePd2uVnJrmwbXdukv0muBb3SfKeJJvar/ckuU9n+4r2ODe09R/auYa/n+TRwErgV9vnf12SJ7b1b9c5znOTXDhRn0gaHgOypPngSOAZwM7AVuCzwEXAMuAQ4JVJnj7Bvp8H9gF2A74FnAlQVavax39dVTtU1bO6O1XVJuCbwHM7q48CPlFVtyV5DvA64PeApcDXgY/1KyDJMuCfgbcAuwAnAZ9MsrTT7EXAccCOwA8mWPcxYCOwO/A84G1JDukcYwXwifY6nTnB9ej6ZWAnmuv4MuDUJA+apP0xwEvb828B3tvZ1vc6dxwFvLV9Lh9mkmvfehVNvx8OPLA978192r0DeARwAPAr7XP5895GVfUL4FPtMce9APi3qromyeOB04E/AH4J+L/A6m7w7Xg9cFB7zv2BA4E/A0hyYPv8XkPTD78BfL+nlkuB44Fvts9/56o6H/gx8LRO0xcC/9Dn/JKGzIAsaT54b1VtqKpbgCcCS6vqTVV1a1VdDvw9cES/Havq9Kq6sQ1IbwT2T7LTNM/7UdpAlSTtOT7abvsD4K+q6tKq2gK8DThgglHkFwJrqmpNVW2tqi8Ba2nC37gzqmpdVW2pqtt619GE2acAJ1fVz6vqQuA0mhA97ptV9en2HLdM4/ndBrypqm6rqjXAz4BHTtL+H6rq21V1E/AG4AVpp0tM4zp/pqq+0db282nU9vvAn1XVZdW4qKp+3G3Q9snLgT+pqp9U1Y00/dD3tUCnP1tHcWd/vhz4v1X1H1V1e1V9CPgFTRDudTTNdbumqjYDf8md/fAy4PSq+lL7XK+squ9O4/kCfIjmtUI7Ov/0Tn2S5hEDsqT5YEPn8Z7A7u3b0tcluY5mJPfBvTu1b9O/vX2b+wbuHMnbdZrn/QTN2+C704wEFs1I8Xgdf9ep4SdAaEYwe+0JPL+n5qcAD5ngOfZbtzswHgLH/aDnfP2OMZkft+F73M3ADpO07x7/B8C9gV2neZ3vaW17AN+bos1S4P7ABZ3r+oV2fT9fBu6X5Entf2QOAM5ut+0JvLqnj/ague69dufOUX7ax+PtplP3RD4CPCvJDjSj21+vqqu28ViS5tC8/ZCHpEWlOo83AFdU1T7T2O8ommkHT6UJbTsBP6UJsr3HvftJq65L8kWasPJo4GNVNb7PBuCtVTWdqQwbaEZfXz5Jm361dNdtAnZJsmMnJD8UuHKKY8ym7hzgh9KMQF/L1Ne5X21T1boB2Bv49iRtrgVuAR5TVVdO0q45YdXWJP9EM4p8NfC5zrUc78+3TnUcmr7YE1jXLj+0Xdete8py+tR3ZZJvAr9LMyL9gWkcR9IQOIIsab75T+CGNB/cu187evnYJE/s03ZHmrfJf0wz0vi2nu1XAw+f4nwfpZl7+1zu+nb3SuC1ufNDfzslef4ExxgfGXx6W+992w/JLZ/i3Heoqg3AucBftfvvR/N2/nQC+mx5YZJ9k9wfeBPNfOzbmfo69zPVtT8NeHOSfdLYL8kvdRtU1Vaa6TV/2/mg3bJJ5qND04f/m2aaRLc//x44vh1dTpIHJHlGkh37HONjwJ+1H/LblWbO8/i9tz8IvCTJIUnu1dbzqAme//Ik2/es/zDwp8D/4s7RbUnzjAFZ0rzSBrJn0bw9fgXNKOJpNKOWvT5M8/b3lcB3gPN6tn8Q2Ld9S/3TE5xyNc2Hz66uqos6dZxN8wGxs9ppBd8GDut3gDbcrqCZCrKZZpTxNdzzn7FHAg+jGa08G/iLdj7zoPwDcAbwI+C+wPjdQ6a6zv1Mde3fDfwT8EXghrb93e5OAZwMrAfOa/vhX5hkHnVV/QdwE82UiM931q+lmYf8PprR7/X0uRNI6y00c8gvBi6h+VDiW9rj/CfwEuBvgeuBf6MZbe71ZZoR6B8lubaz/uy2/dntXG9J81DufDdRkrRYJfkq8JGqOm3YtSx0Sb4H/EFV/cuwa5HUnyPIkiQNSJLn0sxP/vKwa5E0sZEKyElOT3PD+74f6mjnlb03zQ32L27veylJ0tC1o/QfAE5o51dLmqdGaopFkt+guY/nh6vqsX22Hw78Ec29R58E/F1VPWmwVUqSJGmUjdQIclV9jeZepBNZQROeq6rOA3ZO8pBJ2kuSJEl3sdDug7yMu96sfmO77m43Yk9yHM2feOUBD3jAEx71qH536ZEkSdJCdcEFF1xbVXf740MLLSCnz7q+c0iqahWwCmBsbKzWrl07l3VJkiRpnknyg37rR2qKxTRs5K5/CWo5d/71I0mSJGlKCy0grwaOae9mcRBwvX/nXpIkSffESE2xSPIx4GBg1yQbgb8A7g1QVSuBNTR3sFgP3Ezz144kSZKkaRupgFxVR06xvYATBlSOJEmSFqCFNsVCkiRJmhEDsiRJktRhQJYkSZI6DMiSJElShwFZkiRJ6jAgS5IkSR0GZEmSJKnDgCxJkiR1GJAlSZKkDgOyJEmS1GFAliRJkjoMyJIkSVKHAVmSJEnqMCBLkiRJHQZkSZIkqWPkAnKSQ5NclmR9klP6bN8pyWeTXJRkXZKXDKNOSZIkjaaRCshJlgCnAocB+wJHJtm3p9kJwHeqan/gYOBvkmw/0EIlSZI0skYqIAMHAuur6vKquhU4C1jR06aAHZME2AH4CbBlsGVKkiRpVI1aQF4GbOgsb2zXdb0PeDSwCbgEeEVVbe09UJLjkqxNsnbz5s1zVa8kSZJGzKgF5PRZVz3LTwcuBHYHDgDel+SBd9upalVVjVXV2NKlS2e/UkmSJI2kUQvIG4E9OsvLaUaKu14CfKoa64ErgEcNqD5JkiSNuFELyOcD+yTZq/3g3RHA6p42PwQOAUjyYOCRwOUDrVKSJEkja7thF3BPVNWWJCcC5wBLgNOral2S49vtK4E3A2ckuYRmSsbJVXXt0IqWJEnSSBmpgAxQVWuANT3rVnYebwJ+Z9B1SZIkaWEYtSkWkiRJ0pwyIEuSJEkdBmRJkiSpw4AsSZIkdRiQJUmSpA4DsiRJktRhQJYkSZI6DMiSJElShwFZkiRJ6jAgS5IkSR0GZEmSJKnDgCxJkiR1GJAlSZKkDgOyJEmS1GFAliRJkjpGLiAnOTTJZUnWJzllgjYHJ7kwybok/zboGiVJkjS6tht2AfdEkiXAqcDTgI3A+UlWV9V3Om12Bt4PHFpVP0yy23CqlSRJ0igatRHkA4H1VXV5Vd0KnAWs6GlzFPCpqvohQFVdM+AaJUmSNMJGLSAvAzZ0lje267oeATwoyVeTXJDkmH4HSnJckrVJ1m7evHmOypUkSdKoGbWAnD7rqmd5O+AJwDOApwNvSPKIu+1UtaqqxqpqbOnSpbNfqSRJkkbSUAJykiVJ/mQbdt0I7NFZXg5s6tPmC1V1U1VdC3wN2H/bKpUkSdJiM5SAXFW3c/e5w9NxPrBPkr2SbA8cAazuafMZ4NeTbJfk/sCTgEtnVLAkSZIWjWHexeIbSd4H/CNw0/jKqvrWRDtU1ZYkJwLnAEuA06tqXZLj2+0rq+rSJF8ALga2AqdV1bfn8olIkiRp4UhV7xTeAZ04+Uqf1VVVvz3oWsbGxmrt2rWDPq0kSZKGKMkFVTXWu35oI8hV9VvDOrckSZI0kaHdxSLJTknePX6rtSR/k2SnYdUjSZIkwXBv83Y6cCPwgvbrBuD/DbEeSZIkaagf0tu7qp7bWf7LJBcOrRpJkiSJ4Y4g35LkKeMLSX4NuGWI9UiSJElDHUE+HvhwZ97xT4Fjh1iPJEmSNJyAnORewCOrav8kDwSoqhuGUYskSZLUNay/pLcVOLF9fIPhWJIkSfPFMOcgfynJSUn2SLLL+NcQ65EkSZKGOgf5pe2/J3TWFfDwIdQiSZIkAcOdg/zCqvrGMM4vSZIkTWSYc5DfNYxzS5IkSZMZ5hzkLyZ5bpIMsQZJkiTpLoY5B/lVwP2B25P8HAhQVfXAIdYkSZKkRW6YAXkn4Ghgr6p6U5KHAg8ZYj2SJEnSUKdYnAocBBzZLt8IvG+qnZIcmuSyJOuTnDJJuycmuT3J82anXEmSJC0GwwzIT6qqE4CfA1TVT4HtJ9shyRKaYH0YsC9wZJJ9J2j3DuCc2S5akiRJC9swA/JtbZAtgCRLga1T7HMgsL6qLq+qW4GzgBV92v0R8EngmlmsV5IkSYvAMAPye4Gzgd2SvBX4d+BtU+yzDNjQWd7YrrtDkmXA7wIrJztQkuOSrE2ydvPmzfe0dkmSJC1QQ/uQXlWdmeQC4BCaO1g8p6ounWK3freEq57l9wAnV9Xtk91BrqpWAasAxsbGeo8hSZKkRWqYd7Ggqr4LfPce7LIR2KOzvBzY1NNmDDirDce7Aocn2VJVn55JrZIkSVochhqQt8H5wD5J9gKuBI4Ajuo2qKq9xh8nOQP4nOFYkiRJ0zVSAbmqtiQ5kebuFEuA06tqXZLj2+2TzjuWJEmSpjJSARmgqtYAa3rW9Q3GVfXiQdQkSZKkhWOYd7GQJEmS5h0DsiRJktRhQJYkSZI6DMiSJElShwFZkiRJ6jAgS5IkSR0GZEmSJKnDgCxJkiR1GJAlSZKkDgOyJEmS1GFAliRJkjoMyJIkSVKHAVmSJEnqMCBLkiRJHSMXkJMcmuSyJOuTnNJn+9FJLm6/zk2y/zDqlCRJ0mgaqYCcZAlwKnAYsC9wZJJ9e5pdAfxmVe0HvBlYNdgqJUmSNMpGKiADBwLrq+ryqroVOAtY0W1QVedW1U/bxfOA5QOuUZIkSSNs1ALyMmBDZ3lju24iLwM+P6cVSZIkaUHZbtgF3EPps676Nkx+iyYgP2WC7ccBxwE89KEPna36JEmSNOJGbQR5I7BHZ3k5sKm3UZL9gNOAFVX1434HqqpVVTVWVWNLly6dk2IlSZI0ekYtIJ8P7JNkryTbA0cAq7sNkjwU+BTwoqr67yHUKEmSpBE2UlMsqmpLkhOBc4AlwOlVtS7J8e32lcCfA78EvD8JwJaqGhtWzZIkSRotqeo7hXdRGRsbq7Vr1w67DEmSJA1Qkgv6DaSO2hQLSZIkaU4ZkCVJkqQOA7IkSZLUYUCWJEmSOgzIkiRJUocBWZIkSeowIEuSJEkdBmRJkiSpw4AsSZIkdRiQJUmSpA4DsiRJktRhQJYkSZI6DMiSJElShwFZkiRJ6jAgS5IkSR0GZEmSJKnDgCxJkiR1GJAlSZKkDgOyJEmS1JGqGnYNQ5dkM/CDYdexgO0KXDvsIjRr7M+Fxf5ceOzThcX+nFt7VtXS3pUGZM25JGuramzYdWh22J8Li/258NinC4v9ORxOsZAkSZI6DMiSJElShwFZg7Bq2AVoVtmfC4v9ufDYpwuL/TkEzkGWJEmSOhxBliRJkjoMyJqxJLsk+VKS/2n/fdAE7Q5NclmS9UlO6bP9pCSVZNe5r1qTmWmfJnlnku8muTjJ2Ul2Hlz1GjeN77kkeW+7/eIkj5/uvhq8be3PJHsk+UqSS5OsS/KKwVevfmbyPdpuX5Lkv5J8bnBVLw4GZM2GU4B/rap9gH9tl+8iyRLgVOAwYF/gyCT7drbvATwN+OFAKtZUZtqnXwIeW1X7Af8NvHYgVesOU33PtQ4D9mm/jgM+cA/21QDNpD+BLcCrq+rRwEHACfbn8M2wT8e9Arh0jktdlAzImg0rgA+1jz8EPKdPmwOB9VV1eVXdCpzV7jfub4E/BZwUPz/MqE+r6otVtaVtdx6wfI7r1d1N9T1Hu/zhapwH7JzkIdPcV4O1zf1ZVVdV1bcAqupGmkC1bJDFq6+ZfI+SZDnwDOC0QRa9WBiQNRseXFVXAbT/7tanzTJgQ2d5Y7uOJM8Grqyqi+a6UE3bjPq0x0uBz896hZrKdPpnojbT7VsNzkz68w5JHgY8DviPWa9Q99RM+/Q9NANLW+eqwMVsu2EXoNGQ5F+AX+6z6fXTPUSfdZXk/u0xfmdba9O2mas+7TnH62ne3j3znlWnWTBl/0zSZjr7arBm0p/NxmQH4JPAK6vqhlmsTdtmm/s0yTOBa6rqgiQHz3plMiBreqrqqRNtS3L1+Nt47Vs/1/RpthHYo7O8HNgE7A3sBVyUZHz9t5IcWFU/mrUnoLuZwz4dP8axwDOBQ8r7SQ7DpP0zRZvtp7GvBmsm/UmSe9OE4zOr6lNzWKembyZ9+jzg2UkOB+4LPDDJR6rqhXNY76LiFAvNhtXAse3jY4HP9GlzPrBPkr2SbA8cAayuqkuqareqelhVPYzmh8HjDcdDt819Cs0ns4GTgWdX1c0DqFd3N2H/dKwGjmk/KX8QcH07pWY6+2qwtrk/04w+fBC4tKrePdiyNYlt7tOqem1VLW9/bx4BfNlwPLscQdZseDvwT0leRnMXiucDJNkdOK2qDq+qLUlOBM4BlgCnV9W6oVWsqcy0T98H3Af4UvvOwHlVdfygn8RiNlH/JDm+3b4SWAMcDqwHbgZeMtm+Q3gaas2kP4FfA14EXJLkwnbd66pqzSCfg+5qhn2qOeZf0pMkSZI6nGIhSZIkdRiQJUmSpA4DsiRJktRhQJYkSZI6DMiSJElShwFZkiRJ6jAgS9IIS7Jzkj9sH++e5BOzdNwzklwxfk/WJM9Jsm9n+zuT/CjJSbNxPkmaTwzIkjTadgb+EKCqNlXV82bx2K9p/1gBwHOAOwJyVb0GWNl3L0kacQZkSRptbwf2TnJhko8n+TZAkhcn+XSSz7YjwScmeVWS/0pyXpJd2nZ7J/lCkguSfD3Jo3pPkOTJwLOBd7bn2Xugz1CSBsyALEmj7RTge1V1APCanm2PBY4CDgTeCtxcVY8Dvgkc07ZZBfxRVT0BOAl4f+8JqupcYDXNiPIBVfW9OXkmkjRPbDfsAiRJc+YrVXUjcGOS64HPtusvAfZLsgPwZODjScb3uc/gy5Sk+cWALEkL1y86j7d2lrfS/Py/F3BdO/osSWo5xUKSRtuNwI7bsmNV3QBckeT5AGnsP9vnkaRRY0CWpBFWVT8GvtF+OO+d23CIo4GXJbkIWAesmKDdWcBr2g/5+SE9SQtaqmrYNUiS5pkkZwCfq6oJ76uc5I3Az6rqXYOqS5IGwRFkSVI/1wNvHv9DIb2SvBN4IXDTQKuSpAFwBFmSJEnqcARZkiRJ6jAgS5IkSR0GZEmSJKnDgCxJkiR1GJAlSZKkDgOyJEmS1GFAliRJkjoMyJIkSVKHAVmSJEnqMCBL0jQkeWOSj8xg/3VJDp7FkuaNJCuTvGHYdUjSbNlu2AVI0kKT5AxgY1X92fi6qnrM8CqaW1V1/GwfM8mLgd+vqqfM9rElaSqOIEta9JIs+MGCNO411bppHGfeXKskS4Zdg6SFyYAsaVFK8v0kJye5GLgpyXZJdk/yySSbk1yR5I8n2f/jSX6U5PokX0vymHb9ccDRwJ8m+VmSz3bO99T2HLck2aVzrMcluTbJvdvllya5NMlPk5yTZM9J6jgoyblJrktyUXcaR5KvJnlrkm8ANwMPn2Dd7klWJ/lJkvVJXt45xhuTfCLJR5LcALy4Tw1nJHlL+/jgJBuTvDrJNUmuSvKSSep/cZLLk9zYXvOjkzwaWAn8ansNr+uc5wNJ1iS5CfitJPdJ8q4kP0xydTvd435t+12TfK69Nj9J8vXx/xC0fX9le97LkhwyUY2SFh8DsqTF7EjgGcDOwFbgs8BFwDLgEOCVSZ4+wb6fB/YBdgO+BZwJUFWr2sd/XVU7VNWzujtV1Sbgm8BzO6uPAj5RVbcleQ7wOuD3gKXA14GP9SsgyTLgn4G3ALsAJwGfTLK00+xFwHHAjsAPJlj3MWAjsDvwPOBtPYFxBfCJ9jqdOcH16PplYCea6/gy4NQkD+pT/wOA9wKHVdWOwJOBC6vqUuB44JvtNdy5s9tRwFvb2v8deAfwCOAA4Ffac/552/bV7fNaCjyY5rpWkkcCJwJPbM/7dOD703hekhYJA7Kkxey9VbWhqm4Bnggsrao3VdWtVXU58PfAEf12rKrTq+rGqvoF8EZg/yQ7TfO8H6UJ5yRJe46Pttv+APirqrq0qrYAbwMOmGAU+YXAmqpaU1Vbq+pLwFrg8E6bM6pqXVVtqarbetfRhNmnACdX1c+r6kLgNJoQPe6bVfXp9hy3TOP53Qa8qapuq6o1wM+AR07Qdivw2CT3q6qrqmrdFMf+TFV9o6q2Ar8AXg78SVX9pKpupLle4312G/AQYM+2lq9XVQG3A/cB9k1y76r6flV9bxrPS9IiYUCWtJht6DzeE9i9fTv+uvZt/dfRjDzeRZIlSd6e5HvttIPvt5t2neZ5P0EzfWB34DeAohkpHq/j7zo1/AQIzchorz2B5/fU/BSaUNjvOfZbtzswHi7H/aDnfP2OMZkft+F73M3ADr2Nquom4H/TjBZfleSfkzxqimN3a1kK3B+4oPP8v9CuB3gnsB74YjuN45T2vOuBV9L8x+aaJGe1fSFJgAFZ0uJWnccbgCuqaufO145VdXif/Y6imXbwVJqpBA9r16fPce9+0qrrgC8CL2iP9bF2ZHO8jj/oqeN+VXVun0NtAP6hp+0DqurtEzzHfus2Absk2bGz7qHAlVMcY1ZU1TlV9TSaUP9dmlH7yc7ZXX8tcAvwmM7z36mqdmiPfWNVvbqqHg48C3jV+NSRqvpoe4eMPdtjvmPWn5ykkWVAlqTGfwI3tB/eul87SvzYJE/s03ZHmrf3f0wzgvm2nu1XAw+f4nwfBY6hmYv80c76lcBrOx/62ynJ8yc4xkeAZyV5elvvfdsPyS2f4tx3qKoNwLnAX7X770czb3g6c41nJMmDkzy7nYv8C5qpGLe3m68GlifZfqL922kWfw/8bZLd2mMuG583nuSZSX6lncZyQ3vs25M8MslvJ7kP8HOakH17/7NIWowMyJIEVNXtNKOMBwBX0IxOnkYzQtzrwzTTEK4EvgOc17P9gzTzW69L8ukJTrma5kN+V1fVRZ06zqYZzTyrnb7xbeCwCWreQDOS/TpgM82I8mu45z/bj6QZBd8EnA38RTufea7di+aDdJtoppL8JvCH7bYvA+uAHyW5dpJjnEwzjeK89nr9C3fOd96nXf4ZzQcj319VX6WZf/x2mj7+Ec0HLV83a89K0sjLne/qSZIkSXIEWZIkSeoYqYCc5PT2xvPfnmB7kry3vdH9xUkeP+gaJUmSNNpGKiADZwCHTrL9MJo5Z/vQ3AT/AwOoSZIkSQvISAXkqvoazQc5JrIC+HA1zgN2TvKQSdpLkiRJd7HdsAuYZcu4603kN7brruptmOQ4mlFmHvCABzzhUY+a6t70kiRJWkguuOCCa6tqae/6hRaQ02dd39t0VNUqYBXA2NhYrV27di7rkiRJ0jyT5Af91o/UFItp2Ajs0VleTnN/TUmSJGlaFlpAXg0c097N4iDg+qq62/QKSZIkaSIjNcUiyceAg4Fdk2wE/gK4N0BVrQTWAIfT/FWlm4GXDKdSSZIkjaqRCshVdeQU2ws4YUDlSJIkaQFaaFMsJEmSpBkxIEuSJEkdBmRJkiSpw4AsSZIkdRiQJUmSpA4DsiRJktRhQJYkSZI6DMiSJElShwFZkiRJ6jAgS5IkSR0GZEmSJKnDgCxJkiR1GJAlSZKkDgOyJEmS1GFAliRJkjpGLiAnOTTJZUnWJzmlz/adknw2yUVJ1iV5yTDqlCRJ0mgaqYCcZAlwKnAYsC9wZJJ9e5qdAHynqvYHDgb+Jsn2Ay1UkiRJI2ukAjJwILC+qi6vqluBs4AVPW0K2DFJgB2AnwBbBlumJEmSRtWoBeRlwIbO8sZ2Xdf7gEcDm4BLgFdU1dbeAyU5LsnaJGs3b948V/VKkiRpxIxaQE6fddWz/HTgQmB34ADgfUkeeLedqlZV1VhVjS1dunT2K5UkSdJIGrWAvBHYo7O8nGakuOslwKeqsR64AnjUgOqTJEnSiBu1gHw+sE+SvdoP3h0BrO5p80PgEIAkDwYeCVw+0ColSZI0srYbdgH3RFVtSXIicA6wBDi9qtYlOb7dvhJ4M3BGkktopmScXFXXDq1oSZIkjZSRCsgAVbUGWNOzbmXn8SbgdwZdlyRJkhaGUZtiIUmSJM0pA7IkSZLUYUCWJEmSOgzIkiRJUocBWZIkSeowIEuSJEkdBmRJkiSpw4AsSZIkdRiQJUmSpA4DsiRJktRhQJYkSZIbWEzSAAAL/klEQVQ6DMiSJElShwFZkiRJ6jAgS5IkSR0GZEmSJKlj5AJykkOTXJZkfZJTJmhzcJILk6xL8m+DrlGSJEmja7thF3BPJFkCnAo8DdgInJ9kdVV9p9NmZ+D9wKFV9cMkuw2nWkmSJI2iURtBPhBYX1WXV9WtwFnAip42RwGfqqofAlTVNQOuUZIkSSNs1ALyMmBDZ3lju67rEcCDknw1yQVJjul3oCTHJVmbZO3mzZvnqFxJkiSNmlELyOmzrnqWtwOeADwDeDrwhiSPuNtOVauqaqyqxpYuXTr7lUqSJGkkDSUgJ1mS5E+2YdeNwB6d5eXApj5tvlBVN1XVtcDXgP23rVJJkiQtNkMJyFV1O3efOzwd5wP7JNkryfbAEcDqnjafAX49yXZJ7g88Cbh0RgVLkiRp0RjmXSy+keR9wD8CN42vrKpvTbRDVW1JciJwDrAEOL2q1iU5vt2+sqouTfIF4GJgK3BaVX17Lp+IJEmSFo5U9U7hHdCJk6/0WV1V9duDrmVsbKzWrl076NNKkiRpiJJcUFVjveuHNoJcVb81rHNLkiRJExnaXSyS7JTk3eO3WkvyN0l2GlY9kiRJEgz3Nm+nAzcCL2i/bgD+3xDrkSRJkob6Ib29q+q5neW/THLh0KqRJEmSGO4I8i1JnjK+kOTXgFuGWI8kSZI01BHk44EPd+Yd/xQ4doj1SJIkScMJyEnuBTyyqvZP8kCAqrphGLVIkiRJXcP6S3pbgRPbxzcYjiVJkjRfDHMO8peSnJRkjyS7jH8NsR5JkiRpqHOQX9r+e0JnXQEPH0ItkiRJEjDcOcgvrKpvDOP8kiRJ0kSGOQf5XcM4tyRJkjSZYc5B/mKS5ybJEGuQJEmS7mKYc5BfBdwfuD3Jz4EAVVUPHGJNkiRJWuSGGZB3Ao4G9qqqNyV5KPCQIdYjSZIkDXWKxanAQcCR7fKNwPum2inJoUkuS7I+ySmTtHtiktuTPG92ypUkSdJiMMyA/KSqOgH4OUBV/RTYfrIdkiyhCdaHAfsCRybZd4J27wDOme2iJUmStLANMyDf1gbZAkiyFNg6xT4HAuur6vKquhU4C1jRp90fAZ8ErpnFeiVJkrQIDDMgvxc4G9gtyVuBfwfeNsU+y4ANneWN7bo7JFkG/C6wcrIDJTkuydokazdv3nxPa5ckSdICNbQP6VXVmUkuAA6huYPFc6rq0il263dLuOpZfg9wclXdPtkd5KpqFbAKYGxsrPcYkiRJWqSGeRcLquq7wHfvwS4bgT06y8uBTT1txoCz2nC8K3B4ki1V9emZ1CpJkqTFYagBeRucD+yTZC/gSuAI4Khug6raa/xxkjOAzxmOJUmSNF0jFZCrakuSE2nuTrEEOL2q1iU5vt0+6bxjSZIkaSojFZABqmoNsKZnXd9gXFUvHkRNkiRJWjiGeRcLSZIkad4xIEuSJEkdBmRJkiSpw4AsSZIkdRiQJUmSpA4DsiRJktRhQJYkSZI6DMiSJElShwFZkiRJ6jAgS5IkSR0GZEmSJKnDgCxJkiR1GJAlSZKkDgOyJEmS1GFAliRJkjpGLiAnOTTJZUnWJzmlz/ajk1zcfp2bZP9h1ClJkqTRNFIBOckS4FTgMGBf4Mgk+/Y0uwL4zaraD3gzsGqwVUqSJGmUjVRABg4E1lfV5VV1K3AWsKLboKrOraqftovnAcsHXKMkSZJG2KgF5GXAhs7yxnbdRF4GfL7fhiTHJVmbZO3mzZtnsURJkiSNslELyOmzrvo2TH6LJiCf3G97Va2qqrGqGlu6dOkslihJkqRRtt2wC7iHNgJ7dJaXA5t6GyXZDzgNOKyqfjyg2iRJkrQAjNoI8vnAPkn2SrI9cASwutsgyUOBTwEvqqr/HkKNkiRJGmEjNYJcVVuSnAicAywBTq+qdUmOb7evBP4c+CXg/UkAtlTV2LBqliRJ0mhJVd8pvIvK2NhYrV27dthlSJIkaYCSXNBvIHXUplhIkiRJc8qALEmSJHUYkCVJkqQOA7IkSZLUYUCWJEmSOgzIkiRJUocBWZIkSeowIEuSJEkdBmRJkiSpw4AsSZIkdRiQJUmSpA4DsiRJktRhQJYkSZI6DMiSJElShwFZkiRJ6jAgS5IkSR0GZEmSJKnDgCxJkiR1pKqGXcPQJdkM/GDYdSxguwLXDrsIzRr7c2GxPxce+3RhsT/n1p5VtbR3pQFZcy7J2qoaG3Ydmh3258Jify489unCYn8Oh1MsJEmSpA4DsiRJktRhQNYgrBp2AZpV9ufCYn8uPPbpwmJ/DoFzkCVJkqQOR5AlSZKkDgOyJEmS1GFA1owl2SXJl5L8T/vvgyZod2iSy5KsT3JKn+0nJakku8591ZrMTPs0yTuTfDfJxUnOTrLz4KrXuGl8zyXJe9vtFyd5/HT31eBta38m2SPJV5JcmmRdklcMvnr1M5Pv0Xb7kiT/leRzg6t6cTAgazacAvxrVe0D/Gu7fBdJlgCnAocB+wJHJtm3s30P4GnADwdSsaYy0z79EvDYqtoP+G/gtQOpWneY6nuudRiwT/t1HPCBe7CvBmgm/QlsAV5dVY8GDgJOsD+Hb4Z9Ou4VwKVzXOqiZEDWbFgBfKh9/CHgOX3aHAisr6rLq+pW4Kx2v3F/C/wp4KdG54cZ9WlVfbGqtrTtzgOWz3G9urupvudolz9cjfOAnZM8ZJr7arC2uT+r6qqq+hZAVd1IE6iWDbJ49TWT71GSLAeeAZw2yKIXCwOyZsODq+oqgPbf3fq0WQZs6CxvbNeR5NnAlVV10VwXqmmbUZ/2eCnw+VmvUFOZTv9M1Ga6favBmUl/3iHJw4DHAf8x6xXqnpppn76HZmBp61wVuJhtN+wCNBqS/Avwy302vX66h+izrpLcvz3G72xrbdo2c9WnPed4Pc3bu2fes+o0C6bsn0naTGdfDdZM+rPZmOwAfBJ4ZVXdMIu1adtsc58meSZwTVVdkOTgWa9MBmRNT1U9daJtSa4efxuvfevnmj7NNgJ7dJaXA5uAvYG9gIuSjK//VpIDq+pHs/YEdDdz2KfjxzgWeCZwSHnD9WGYtH+maLP9NPbVYM2kP0lyb5pwfGZVfWoO69T0zaRPnwc8O8nhwH2BByb5SFW9cA7rXVScYqHZsBo4tn18LPCZPm3OB/ZJsleS7YEjgNVVdUlV7VZVD6uqh9H8MHi84XjotrlPoflkNnAy8OyqunkA9eruJuyfjtXAMe0n5Q8Crm+n1ExnXw3WNvdnmtGHDwKXVtW7B1u2JrHNfVpVr62q5e3vzSOALxuOZ5cjyJoNbwf+KcnLaO5C8XyAJLsDp1XV4VW1JcmJwDnAEuD0qlo3tIo1lZn26fuA+wBfat8ZOK+qjh/0k1jMJuqfJMe321cCa4DDgfXAzcBLJtt3CE9DrZn0J/BrwIuAS5Jc2K57XVWtGeRz0F3NsE81x/xT05IkSVKHUywkSZKkDgOyJEmS1GFAliRJkjoMyJIkSVKHAVmSJEnqMCBLkiRJHQZkSRphSXZO8oft492TfGKWjntGkivG78ma5DlJ9u1sf2eSHyU5aTbOJ0nziQFZkkbbzsAfAlTVpqp63iwe+zXtHysAeA5wR0CuqtcAK/vuJUkjzoAsSaPt7cDeSS5M8vEk3wZI8uIkn07y2XYk+MQkr0ryX0nOS7JL227vJF9IckGSryd5VO8JkjwZeDbwzvY8ew/0GUrSgBmQJWm0nQJ8r6oOAF7Ts+2xwFHAgcBbgZur6nHAN4Fj2jargD+qqicAJwHv7z1BVZ0LrKYZUT6gqr43J89EkuaJ7YZdgCRpznylqm4EbkxyPfDZdv0lwH5JdgCeDHw8yfg+9xl8mZI0vxiQJWnh+kXn8dbO8laan//3Aq5rR58lSS2nWEjSaLsR2HFbdqyqG4ArkjwfII39Z/s8kjRqDMiSNMKq6sfAN9oP571zGw5xNPCyJBcB64AVE7Q7C3hN+yE/P6QnaUFLVQ27BknSPJPkDOBzVTXhfZWTvBH4WVW9a1B1SdIgOIIsSerneuDN438opFeSdwIvBG4aaFWSNACOIEuSJEkdjiBLkiRJHQZkSZIkqcOALEmSJHUYkCVJkqSO/w8xpre2uvIKZgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x720 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total simulation time =  94.009917293\n",
      "ploynomial degree =  4\n",
      "degree of freedom =  400\n",
      "maximum relative error in particle velocity =  0.00026628844582363175\n",
      "maximum relative error in stress =  0.0003766045233679544\n"
     ]
    }
   ],
   "source": [
    "# Simulation start time\n",
    "start = timeit.default_timer()\n",
    "\n",
    "# Tic\n",
    "iplot = 20\n",
    "\n",
    "# Physical domain x = [ax, bx] (km)\n",
    "ax = 0.0                     # (km)\n",
    "bx = 20.0                    # (km)\n",
    "\n",
    "# Choose quadrature rules and the corresponding nodes\n",
    "# We use Gauss-Legendre-Lobatto (Lobatto) or  Gauss-Legendre (Legendre) quadrature rule.\n",
    "\n",
    "#node = 'Lobatto'\n",
    "node = 'Legendre'\n",
    "\n",
    "if node not in ('Lobatto', 'Legendre'):\n",
    "    print('quadrature rule not implemented. choose node = Legendre or node = Lobatto')\n",
    "    exit(-1)\n",
    "\n",
    "# Polynomial degree N: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]\n",
    "N = 4                        # Lagrange polynomial degree\n",
    "NP = N+1                     # quadrature nodes per element\n",
    "\n",
    "if N < 1 or N > 12: \n",
    "    print('polynomial degree not implemented. choose N>= 1 and N <= 12')\n",
    "    exit(-1)\n",
    "\n",
    "# Degrees of freedom to resolve the wavefield    \n",
    "deg_of_freedom = 400 # = num_element*NP\n",
    "\n",
    "# Estimate the number of elements needed for a given polynomial degree and degrees of freedom\n",
    "num_element = round(deg_of_freedom/NP) # = deg_of_freedom/NP\n",
    "\n",
    "# Initialize the mesh\n",
    "y = np.zeros(NP*num_element)\n",
    "\n",
    "# Generate num_element dG elements in the interval [ax, bx]\n",
    "x0 = np.linspace(ax,bx,num_element+1)\n",
    "dx = np.diff(x0)             # element sizes\n",
    "\n",
    "# Generate Gauss quadrature nodes (psi): [-1, 1] and weights (w)\n",
    "if node == 'Legendre':\n",
    "    GL_return = quadraturerules.GL(N)\n",
    "    psi = GL_return['xi']\n",
    "    w = GL_return['weights'];\n",
    "    \n",
    "if node == 'Lobatto':\n",
    "    gll_return = quadraturerules.gll(N)\n",
    "    psi = gll_return['xi']\n",
    "    w = gll_return['weights']\n",
    "\n",
    "# Use the Gauss quadrature nodes (psi) generate the mesh (y)\n",
    "for i in range (1,num_element+1):\n",
    "    for j in range (1,(N+2)):\n",
    "        y[j+(N+1)*(i-1)-1] = dx[i-1]/2.0 * (psi[j-1] + 1.0) +x0[i-1]\n",
    "\n",
    "# Overide with the exact degrees of freedom        \n",
    "deg_of_freedom = len(y);\n",
    "\n",
    "# Generate the spectral difference operator (D) in the reference element: [-1, 1]  \n",
    "D = specdiff.derivative_GL(N, psi, w)\n",
    "\n",
    "# Boundary condition reflection coefficients \n",
    "r0 = 0                            # r=0:absorbing, r=1:free-surface, r=-1: clamped \n",
    "rn = 0                            # r=0:absorbing, r=1:free-surface, r=-1: clamped \n",
    "\n",
    "# Initialize the wave-fields\n",
    "L = 0.5*(bx-ax)\n",
    "delta = 0.01*(bx-ax)\n",
    "x0 = 0.5*(bx+ax)\n",
    "\n",
    "omega = 4.0\n",
    "#u = np.sin(omega*np.pi*y/L)                                          # Sine function\n",
    "u = 1/np.sqrt(2.0*np.pi*delta**2)*np.exp(-(y-x0)**2/(2.0*delta**2))   # Gaussian\n",
    "    \n",
    "u = np.transpose(u)\n",
    "v = np.zeros(len(u))\n",
    "\n",
    "U = np.zeros(len(u))\n",
    "V = np.zeros(len(u))\n",
    "\n",
    "print('points per wavelength: ', round(delta*deg_of_freedom/(2*L)))\n",
    "\n",
    "# Material parameters\n",
    "cs = 3.464                               # shear wave speed (km/s)\n",
    "rho = 0*y  + 2.67                        # density (g/cm^3)\n",
    "mu = 0*y + rho * cs**2                   # shear modulus (GPa)\n",
    "Zs = rho*cs                              # shear impedance   \n",
    "\n",
    "# Time stepping parameters\n",
    "cfl = 0.5                                # CFL number\n",
    "dt = (0.25/(cs*(2*N+1)))*min(dx)         # time-step (s)\n",
    "t = 0.0                                  # initial time\n",
    "Tend = 2                                # final time (s)\n",
    "n = 0                                    # counter\n",
    "\n",
    "# Difference between analyticla and numerical solutions\n",
    "EV = [0]                                 # initialize errors in V (velocity)\n",
    "EU = [0]                                 # initialize errors in U (stress)\n",
    "T = [0]                                  # later append every time steps to this\n",
    "\n",
    "# Initialize animated plot for velocity and stress\n",
    "fig1 = plt.figure(figsize=(10,10))\n",
    "ax1 = fig1.add_subplot(4,1,1)\n",
    "line1 = ax1.plot(y, u, 'r', y, U, 'k--')\n",
    "plt.title('numerical vs exact')\n",
    "plt.xlabel('x[km]')\n",
    "plt.ylabel('velocity [m/s]')\n",
    "\n",
    "ax2 = fig1.add_subplot(4,1,2)\n",
    "line2 = ax2.plot(y, v, 'r', y, V, 'k--')\n",
    "plt.title('numerical vs exact')\n",
    "plt.xlabel('x[km]')\n",
    "plt.ylabel('stress [MPa]')\n",
    "\n",
    "# Initialize error plot (for velocity and stress)\n",
    "ax3 = fig1.add_subplot(4,1,3)\n",
    "line3 = ax3.plot(T, EV, 'r')\n",
    "plt.title('relative error in particle velocity')\n",
    "plt.xlabel('time[t]')\n",
    "ax3.set_ylim([10**-5, 1])\n",
    "plt.ylabel('error')\n",
    "\n",
    "ax4 = fig1.add_subplot(4,1,4)\n",
    "line4 = ax4.plot(T, EU, 'r') \n",
    "plt.ylabel('error')\n",
    "plt.xlabel('time[t]')\n",
    "ax4.set_ylim([10**-5, 1])\n",
    "plt.title('relative error in stress')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.ion()\n",
    "plt.show()\n",
    "\n",
    "A = (np.linalg.norm(1/np.sqrt(2.0*np.pi*delta**2)*0.5*Zs*(np.exp(-(y+cs*(t+1*0.5)-x0)**2/(2.0*delta**2))\\\n",
    "                                               - np.exp(-(y-cs*(t+1*0.5)-x0)**2/(2.0*delta**2)))))\n",
    "B = (np.linalg.norm(u))\n",
    "\n",
    "# Loop through time and evolve the wave-fields using ADER time-stepping scheme of N+1 order of accuracy\n",
    "time_integrator = 'ADER'\n",
    "\n",
    "for t in utils.drange (0.0, Tend+dt,dt):\n",
    "    n = n+1\n",
    "    \n",
    "    # ADER time-integrator\n",
    "    if time_integrator  in ('ADER'):\n",
    "        ADER_Wave_dG_return = timeintegrate.ADER_Wave_dG(u,v,D,NP,num_element,dx,w,psi,t,r0,rn,dt,rho,mu)\n",
    "        u = ADER_Wave_dG_return['Hu']\n",
    "        v = ADER_Wave_dG_return['Hv']\n",
    "     \n",
    "    # Runge-Kutta time-integrator\n",
    "    if time_integrator  in ('RK'):\n",
    "        RK4_Wave_dG_return = timeintegrate.RK4_Wave_dG(u,v,D,NP,num_element,dx,w,psi,t,r0,rn,dt,rho,mu)\n",
    "        u = RK4_Wave_dG_return['Hu']\n",
    "        v = RK4_Wave_dG_return['Hv']\n",
    "\n",
    "    # Analytical sine wave (use it when sine function is choosen above)\n",
    "    #U = 0.5*(np.sin(omega*np.pi/L*(y+cs*(t+1*dt))) + np.sin(omega*np.pi/L*(y-cs*(t+1*dt))))\n",
    "    #V = 0.5*Zs*(np.sin(omega*np.pi/L*(y+cs*(t+1*dt))) - np.sin(omega*np.pi/L*(y-cs*(t+1*dt))))\n",
    "    \n",
    "    # Analytical Gaussian\n",
    "    U =  1/np.sqrt(2.0*np.pi*delta**2)*0.5*(np.exp(-(y+cs*(t+1*dt)-x0)**2/(2.0*delta**2))\\\n",
    "                                            + np.exp(-(y-cs*(t+1*dt)-x0)**2/(2.0*delta**2)))\n",
    "    V =  1/np.sqrt(2.0*np.pi*delta**2)*0.5*Zs*(np.exp(-(y+cs*(t+1*dt)-x0)**2/(2.0*delta**2))\\\n",
    "                                               - np.exp(-(y-cs*(t+1*dt)-x0)**2/(2.0*delta**2)))\n",
    "    EV.append(np.linalg.norm(V-v)/A)\n",
    "    EU.append(np.linalg.norm(U-u)/B)\n",
    "    \n",
    "    T.append(t)\n",
    "\n",
    "    # Updating plots\n",
    "    if n % iplot == 0: \n",
    "        for l in line1:\n",
    "            l.remove()\n",
    "            del l               \n",
    "        for l in line2:\n",
    "            l.remove()\n",
    "            del l\n",
    "        for l in line3:\n",
    "            l.remove()\n",
    "            del l               \n",
    "        for l in line4:\n",
    "            l.remove()\n",
    "            del l \n",
    "\n",
    "        # Display lines\n",
    "        line1 = ax1.plot(y, u, 'r', y, U, 'k--')\n",
    "        ax1.legend(iter(line1),('Numerical', 'Analytical'))\n",
    "        line2 = ax2.plot(y, v, 'r', y, V, 'k--')\n",
    "        ax2.legend(iter(line2),('Numerical', 'Analytical'))\n",
    "        line3 = ax3.plot(T, EU, 'k--')\n",
    "        ax3.set_yscale(\"log\")#, nonposx='clip')\n",
    "        line4 = ax4.plot(T, EV, 'k--')\n",
    "        ax4.set_yscale(\"log\")#, nonposx='clip')\n",
    "        plt.gcf().canvas.draw()\n",
    "       \n",
    "plt.ioff()\n",
    "plt.show()\n",
    "\n",
    "# Simulation end time\n",
    "stop = timeit.default_timer()\n",
    "print('total simulation time = ', stop - start)                   # print the time required for simulation\n",
    "print('ploynomial degree = ', N)                                  # print the polynomial degree used\n",
    "print('degree of freedom = ', deg_of_freedom)                     # print the degree of freedom\n",
    "print('maximum relative error in particle velocity = ', max(EU))  # max. relative error in particle velocity\n",
    "print('maximum relative error in stress = ', max(EV))             # max. relative error in stress"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}