Computational_Seismology_Igel/sp_elasticwave1D.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)\"> SBP-SAT finite difference method for the 1D elastic wave equation in first order form </div>\n",
    "        </div>\n",
    "    </div>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook is based on the paper [Dynamic earthquake rupture simulations on nonplanar faults embedded in 3D geometrically complex, heterogeneous Earth models](https://pangea.stanford.edu/~edunham/publications/Duru_Dunham_FD3d_JCP16.pdf), and on the theory of summation-by-parts (SBP) finite difference methods and weak implementation of boundary conditions using the simultaneous-approximation-term (SAT).\n",
    "\n",
    "\n",
    "##### Authors:\n",
    "* Kenneth Duru\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic Equations ##\n",
    "\n",
    "We consider in this notebook a system of hyperbolic PDE, the elastic wave equation, in 1D. The source-free elastic wave equation in a heterogeneous 1D medium is \n",
    "\n",
    "\\begin{align}\n",
    "\\rho(x)\\frac{\\partial v(x,t)}{\\partial t} -\\frac{\\partial\\sigma(x,t)}{\\partial x} & = 0\\\\\n",
    "\\frac{1}{\\mu(x)}\\frac{\\partial\\sigma(x,t)}{\\partial t} -\\frac{\\partial v(x,t)}{\\partial x} & = 0  \n",
    "\\end{align}\n",
    "\n",
    "with $\\rho(x)$ the density, $\\mu(x)$ the shear modulus and $x \\in [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",
    "We introduce the mechanical energy defined by\n",
    "\\begin{equation}\n",
    "E(t) = \\int_0^L{\\left(\\frac{\\rho(x)}{2} v^2(x, t) + \\frac{1}{2\\mu(x)}\\sigma^2(x, t)\\right) dx},\n",
    "\\end{equation}\n",
    "\n",
    "where $E(t)$ is the sum of the kinetic energy and the strain energy.\n",
    "We have \n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{d E(t)}{dt} = -v(0, t)\\sigma(0, t) + v(L, t)\\sigma(L, t) \\le 0.\n",
    "\\end{equation}\n",
    "\n",
    "From the boundary conditions, it is easy to check that  $v(0, t)\\sigma(0, t) \\ge 0$ and $v(L, t)\\sigma(L, t) \\le 0$, for all $|r_0|, |r_n| \\le 1$. This energy loss through the boundaries is what the numerical method should emulate. \n",
    "\n",
    "1) Discretize the spatial domain $x$ into $N$ discrete nodes with the uniform spatial step $\\Delta{x} = L/(N-1)$, denote the unknown fields at the nodes: $\\mathbf{v}\\left(t\\right) = [v_1\\left(t\\right), v_2\\left(t\\right), \\cdots, v_N\\left(t\\right)]$, and $\\boldsymbol{\\sigma}\\left(t\\right) = [\\sigma_1\\left(t\\right), \\sigma_2\\left(t\\right), \\cdots, \\sigma_N\\left(t\\right)]$.\n",
    "\n",
    "\n",
    "2) At the grid-point $x_j = (j-1)\\Delta{x}$: Approximate the spatial derivative by a finite difference operator $\\partial v/\\partial x\\Big|_{x = x_j} \\approx \\left(\\mathbf{D}\\mathbf{v}\\right)_j $. Here $\\mathbf{D}$ is a finite difference matrix satisfying the summation-by-parts property:\n",
    "\n",
    "\\begin{align}\n",
    "\\mathbf{H}\\mathbf{D} = \\mathbf{Q}, \\quad \\mathbf{Q} + \\mathbf{Q} = \\left(\\boldsymbol{e}_{N}\\boldsymbol{e}_{N}^T -\\boldsymbol{e}_{1}\\boldsymbol{e}_{1}^T\\right), \\quad \\mathbf{H}^T = \\mathbf{H} > 0,\n",
    "\\end{align}\n",
    "\n",
    "where, $\\boldsymbol{e}_{0} = [1, 0, \\dots,  0 ]^T, \\quad  \\boldsymbol{e}_{L} = [ 0, 0, \\dots,  1 ]^T$ and $\\mathbf{H}$ defines a dicrete norm. We consider only diagonal norm SBP operators with $H_{jj} = h_j > 0$, and define the quadrature rule\n",
    "\n",
    "\\begin{equation}\n",
    " \\sum_{i = 1}^{N} f(x_j)h_j \\approx \\int_{0}^{L}f(x) dx.\n",
    "\\end{equation}\n",
    "\n",
    "The second order accurate SBP operator for the first derivative is:\n",
    "\\begin{align}\n",
    "\\left(\\mathbf{D}\\mathbf{v}\\right)_j = \\frac{v_{j+1}-v_{j-1}}{2 \\Delta{x}}, \\quad j = 2, 3, \\cdots N-1, \\quad\n",
    "\\left(\\mathbf{D}\\mathbf{v}\\right)_1 = \\frac{v_{2}-v_{1}}{\\Delta{x}},\\quad\n",
    "\\left(\\mathbf{D}\\mathbf{v}\\right)_N = \\frac{v_{N}-v_{N-1}}{\\Delta{x}}, \\quad j = N.\n",
    "\\end{align}\n",
    "\n",
    "Note that the interior stencils are centered, with second order accuracy and the boundary stencils are one-sided and first order accurate. \n",
    "\n",
    "Higher order SBP operators can be found in the book: High Order Difference Methods for Time Dependent PDE, by B. Gustafsson. In this notebook we implement SBP operators with interior accuracy 2, 4 and 6. The implementation of the spatial derivative operators can be found in the file first_derivative_sbp_operators.py\n",
    "\n",
    "To construct a stable  semi-discrete approximation we replace the spatial derivatives by the SBP operators, and add the boundary conditions as SAT-terms with special penalty weights having:\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{d \\mathbf{v}(t)}{d t} = {\\boldsymbol{\\rho}}^{-1}\n",
    "\\left(\\mathbf{D} \\boldsymbol{\\sigma}(t) - \\underbrace{\\mathbf{H}^{-1}\\left(\\tau_{11}\\boldsymbol{e}_{1}B_0\\left(v_1, \\sigma_1, Z_{s}, r_0\\right) + \\tau_{12}\\boldsymbol{e}_{N}B_L\\left(v_N, \\sigma_N, Z_{s}, r_n\\right)\\right)}_{SAT \\to 0}\\right),\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{d \\boldsymbol{\\sigma}(t)}{d t} = \\boldsymbol{\\mu}\n",
    "\\left(\\mathbf{D} \\mathbf{v}(t) + \\underbrace{\\mathbf{H}^{-1}\\left(\\tau_{21}\\frac{\\boldsymbol{e}_{1}}{Z_{s}}B_0\\left(v_1, \\sigma_1, Z_{s}, r_0\\right) - \\tau_{22}\\frac{\\boldsymbol{e}_{N}}{Z_{s}}B_L\\left(v_N, \\sigma_N, Z_{s}, r_n\\right)\\right)}_{SAT \\to 0}\\right).\n",
    "\\end{align}\n",
    "Here $\\tau_{ij}$ are penalty parameters determined by requiring stability.\n",
    "\n",
    "Approximate the mechanical energy by the above quadrature rule, having \n",
    "\\begin{align}\n",
    "\\mathcal{E}( t) = \\sum_{j}^{N}\\frac{1}{2}\\left(\\rho_jv_j^2 + \\frac{1}{\\mu_j}\\sigma_j^2\\right)h_j > 0.\n",
    "\\end{align}\n",
    "\n",
    "By chosing the penalty parameters $\\tau_{ij} = 1$, the semi-discrete approximation satisfies the energy estimate:\n",
    "\\begin{align}\n",
    "\\frac{d \\mathcal{E}( t)}{d t} = -\\frac{1}{2}\\left(\\left(1-r_0\\right)Zv_1^2 + \\frac{\\left(1+r_0\\right)}{Z}\\sigma_1^2 +\n",
    "\\left(1-r_n\\right)Zv_N^2 + \\frac{\\left(1+r_n\\right)}{Z}\\sigma_N^2\\right) \\le 0.\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "3) The discrete mechanical energy can never grow in time, and thus the semidiscrete numerical approximation is asymptotically stable.\n",
    "\n",
    "4)  Time integration can be performed using any stable time stepping scheme. This notebook implements the fourth order accurate Runge-Kutta method. \n",
    "\n",
    "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",
    "\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 numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import time_integrator\n",
    "import rate\n",
    "import utils\n",
    "import timeit\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
    ]
   },
   "outputs": [],
   "source": [
    "# Initializations\n",
    "L = 10.0         # length of the domain (km)\n",
    "t = 0.0          # initial time\n",
    "tend = 1.45      # final time\n",
    "nx = 501        # grid points in x                                                                                                                       \n",
    "dx = L/(nx-1)    # grid increment in x\n",
    "cs = 3.464       # velocity (km/s) (can be an array)                                                                                                             \n",
    "iplot = 5       # snapshot frequency\n",
    "rho = 2.6702     # density [g/cm^3]\n",
    "mu = rho*cs**2   # shear modulus [GPa]\n",
    "Zs = rho*cs      # shear impedance \n",
    "\n",
    "order = 6        # order of accuracy\n",
    "\n",
    "\n",
    "#Initialize the domain\n",
    "y = np.zeros((nx, 1))\n",
    "\n",
    "# Initial particle velocity perturbation and discretize the domain\n",
    "for j in range(0, nx):\n",
    "    y[j, :] = j*dx                                             # discrete domain\n",
    "\n",
    "\n",
    "\n",
    "# Time stepping parameters\n",
    "cfl = 1.0                         # CFL number\n",
    "dt = (cfl/cs)*dx                  # Time step\n",
    "nt = int(round(tend/dt))          # number of time steps\n",
    "n = 0                             # counter\n",
    "\n",
    "# Boundary condition reflection coefficients \n",
    "r0 = 1                            # r=0:absorbing, r=1:free-surface, r=-1: clamped \n",
    "r1 = 1                            # r=0:absorbing, r=1:free-surface, r=-1: clamped\n",
    "\n",
    "# penalty parameters\n",
    "tau_11 = 1 \n",
    "tau_12 = 1\n",
    "tau_21 = 1 \n",
    "tau_22 = 1\n",
    "\n",
    "# Initialize: particle velocity (v); and shear stress (s)\n",
    "v = np.zeros((nx, 1))\n",
    "s = np.zeros((nx, 1))\n",
    "\n",
    "U = np.zeros((nx, 1))\n",
    "V = np.zeros((nx, 1))\n",
    "U_t = np.zeros((nx, 1))\n",
    "V_t = np.zeros((nx, 1))\n",
    "U_x = np.zeros((nx, 1))\n",
    "V_x = np.zeros((nx, 1))\n",
    "\n",
    "                                \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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "code_folding": [
     0
    ],
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAALICAYAAABiqwZ2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeZydZX3//9fbsMkOErYkkKBBBdRoR9zq0ooVcAm2LliVRS1SRaXFVlza8rO25WvVIi2FIiJYNhGlRAUBU2lrK5YEMBDZYgATCBCQPbIEPr8/zh17ZzKZOclk5iQzr+fjcR7n3Nd2f+6cSfLJleu+7lQVkiRJkjqe0esAJEmSpPWJCbIkSZLUYoIsSZIktZggS5IkSS0myJIkSVKLCbIkSZLUYoIsSRuoJJ9Octowx5iapJJstK7ikqQNXdwHWZLGryRTgVuBjatqeW+jWXeSFDC9qhb0OhZJGx5nkCVpA+SMrySNHBNkSRpEktuSfCLJvCQPJvlmks2ausOS/Lhf+0rynObzGUn+OcklSR5J8t9Jdk5yQpL7k9yY5MWtvrsm+XaSpUluTfKxVt1xSS5IclaSh4DDmrKzWm1+O8n/JHkgyaIkhzXlb0pyTZKHmvLjurz2Y5Nc0K/sK0lObF3/wiQPN/G+ZzXjPKMZ6xdJ7ktyfpLtm7p3NWNs3RwfkOSuJBNb51vUxD43yatb405olpn8oolhbpIpSf6zafKz5tf9Xd1cryStYIIsSUN7J7A/MA14IXDYGvb9LLAD8DjwE+Dq5vgC4MvQSSKB7wI/AyYBrweOTvLG1lgzmz7bAme3T5JkN+AS4B+BicAM4Nqm+lHgkKbfm4A/TnJQF7GfCxzYSl4nNNdzTpItgBOBA6pqK+CVrfP19zHgIOC1wK7A/cBJAFX1zebX5MQkzwK+BnywqpY2fa9qrmV74BzgWyv+gQL8KfBu4EBga+D9wLKqek1T/6Kq2rI5hyR1zQRZkoZ2YlXdWVW/opPEzliDvhdW1dyqegy4EHisqr5RVU8B3wRWzCC/FJhYVZ+rqieqaiHwVeDg1lg/qap/q6qnq+rX/c7zHuCHVXVuVT1ZVfdV1bUAVXVFVV3X9JtHJ/F97VCBV9XtdJL5Fcn079JJQK9sjp8G9knyzKpaUlXzVzPUh4DPVNXiqnocOA54e2uZyEeasa8AvltV32vFcFZzLcur6kvApsBzm+oPAp+tqpuq42dVdd9Q1yVJQzFBlqSh3dX6vAzYcg363t36/OsBjleMtTuwa7M84oEkDwCfBnZqtV80yHmmAL8YqCLJy5L8qFm68SBwJJ0Z7G6cQ2eWFuAPm2Oq6lHgXc1YS5J8P8nzVjPG7sCFreu6AXiK5tqq6gHgW8A+wJf6xX5Mkhua5S0PANu0Yl/tNUvScJggS9LaexTYfMVBkp2HMdYi4Naq2rb12qqqDmy1GWzboUXAs1dTdw4wC5hSVdsApwDpMq5vAa9LMhl4WzNWJ5iqS6vqDcAuwI10ZrxXF9sB/a5ts6q6AyDJDDrLI86ls2yDpvzVwCfpLOvYrqq2BR5sxT7YNUvSWjNBlqS19zNg7yQzmnWxxw1jrP8FHkryySTPbG5A2yfJS7vsfzawX5J3JtkoybOaxBNgK+BXVfVYkn3pzAR3pVkLfAXwdToJ/A0ASXZK8tZmLfLjwCN0ZoUHcgrwN0l2b/pOTDKz+bwZcBad2fLDgUlJPtyKezmwFNgoyV/SWWu8wmnAXyeZno4XNuuYoTNTv0e31ylJbSbIkrSWqupm4HPAD4FbgB8P3mPQsZ4C3kJnffOtwL10EsBtuuz/Szo3qx0D/IrODXMvaqo/DHwuycPAXwLnr2F45wD70Zo9pvP3xzHAnc35XtucZyBfoTODfVkTw5XAy5q6vwMWV9XJzfrk9wKfTzIduJTOjYc3A7cDj7HyMpMvN9dyGfAQnRv8ntnUHQec2SzreOcaXq+kcc4HhUiSJEktziBLkiRJLSbIkiRJUosJsiRJktRigixJkiS1bDR0k7Fvhx12qKlTp/Y6DEmSJI2iuXPn3ltVE/uXmyADU6dOZc6cOb0OQ5IkSaMoye0DlbvEQpIkSWrpaYKcZP8kNyVZkOTYAeqT5MSmfl6Sl7TqbktyXZJrk8xplW+f5PIktzTv243W9UiSJGnD17MEOckE4CTgAGAv4N1J9urX7ABgevM6Aji5X/3vVNWMquprlR0LzK6q6cDs5liSJEnqSi9nkPcFFlTVwqp6AjgPmNmvzUzgG9VxJbBtkl2GGHcmcGbz+UzgoHUZtCRJksa2XibIk4BFrePFTVm3bQq4LMncJEe02uxUVUsAmvcdBzp5kiOSzEkyZ+nSpcO4DEmSJI0lvUyQM0BZrUGbV1XVS+gsw/hIktesycmr6tSq6quqvokTV9ndQ5IkSeNULxPkxcCU1vFk4M5u21TVivd7gAvpLNkAuHvFMozm/Z51HrkkSZLGrF4myFcB05NMS7IJcDAwq1+bWcAhzW4WLwcerKolSbZIshVAki2A3wOub/U5tPl8KHDRSF+IJEmSxo6ePSikqpYnOQq4FJgAnF5V85Mc2dSfAlwMHAgsAJYBhzfddwIuTAKdazinqn7Q1B0PnJ/kA8AvgXeM0iVJkiRpDEhV/2W/409fX1/5JD1JkqTxJcncftsFAz5JT5IkSVqJCbIkSZLUYoIsSZIktZggS5IkSS0myJIkSVKLCbIkSZLUYoIsSZIktZggS5IkSS0myJIkSVKLCbIkSZLUYoIsSZIktZggS5IkSS0myJIkSVKLCbIkSZLUYoIsSZIktZggS5IkSS0bra4iyZ920f/RqvqXdRiPJEmS1FODzSD/GbAlsNUgr2NGOkBJkiRpNK12Bhn416r63GCdk2wxnJMn2R/4CjABOK2qju9Xn6b+QGAZcFhVXZ1kCvANYGfgaeDUqvpK0+c44I+Apc0wn66qi4cTpyRJksaP1SbIVfXnQ3Xups3qJJkAnAS8AVgMXJVkVlX9vNXsAGB683oZcHLzvhw4pkmWtwLmJrm81fcfquqLaxubJEmSxq8hb9JL8vEkW6fja0muTvJ76+Dc+wILqmphVT0BnAfM7NdmJvCN6rgS2DbJLlW1pKquBqiqh4EbgEnrICZJkiSNc93sYvH+qnoI+D1gInA4cPzgXboyCVjUOl7MqknukG2STAVeDPy0VXxUknlJTk+y3UAnT3JEkjlJ5ixdunSgJpIkSRqHukmQ07wfCHy9qn7WKhuOgcaoNWmTZEvg28DRTRIPnWUYzwZmAEuALw108qo6tar6qqpv4sSJaxq7JEmSxqhuEuS5SS6jkyBf2qz5fXodnHsxMKV1PBm4s9s2STamkxyfXVXfWdGgqu6uqqeq6mngq3SWckiSJEldWW2CnGTFDXwfAI4FXlpVy4BN6CyzGK6rgOlJpiXZBDgYmNWvzSzgkGb988uBB6tqSbO7xdeAG6rqy/3i3qV1+Dbg+nUQqyRJksaJwbZ5uzLJYuAHwA+q6gGAqroPuG+4J66q5UmOAi6ls83b6VU1P8mRTf0pwMV0Zq4X0NnmbUVi/irgfcB1Sa5tylZs5/aFJDPoLMW4DfjQcGOVJEnS+JGq/st+W5XJ7nS2Wtufzs1xPwYuAf6jqh4flQhHQV9fX82ZM6fXYUiSJGkUJZlbVX39ywddg1xVt1fVKVV1EPBK4LvAfsB/Jfn+yIQqSZIk9c5gSyxWUlVPAv/evEjivsOSJEkac7p5UMibk1yT5P4kDyV5OMlDVXXHaAQoSZIkjaZuZpBPAH4fuK4GW7AsSZIkjQHd7IO8CLje5FiSJEnjQTczyH8OXJzkP4Df7FzRf/9hSZIkaSzoJkH+G+ARYDM6DwmRJEmSxqxuEuTtq+r3RjwSSZIkaT3QzRrkHyYxQZYkSdK40E2C/BHgB0l+3d7mbaQDkyRJknphyCUWVbXVaAQiSZIkrQ9WO4OcZOehOnfTRpIkSdqQDLbE4uIu+nfTRpIkSdpgDLbE4kVDrDUO4FpkSZIkjSmrTZCrasJoBiJJkiStD7rZxUKSJEkaN0yQJUmSpBYTZEmSJKllyAQ5yReT7D0SJ0+yf5KbkixIcuwA9UlyYlM/L8lLhuqbZPsklye5pXnfbiRilyRJ0tjUzQzyjcCpSX6a5Mgk26yLEyeZAJwEHADsBbw7yV79mh0ATG9eRwAnd9H3WGB2VU0HZjfHkiRJUle6eZLeacBpSZ4LHA7MS/LfwFer6kfDOPe+wIKqWgiQ5DxgJvDzVpuZwDeqqoArk2ybZBdg6iB9ZwKva/qfCVwBfHIYca57Rx/N6844Y5Xid06cyIcnTWLZU09x4HXXrVJ/2M47c9jOO3Pvk0/y9vnzV6n/41135V077siixx7jfTfeuEr9MZMn85YdduCmZcv40M03r1L/2d13Z7/ttuPaRx7h6AULVqn/22nTeOU22/A/Dz7Ip2+9dZX6E57zHGZsuSU/vP9+Pn/77avU/8uee/LczTfnu/fey5cWL16l/l+f9zymbLYZ37znHk6+885V6i/Ye2922HhjzrjrLs64665V6i9+wQvYfMIE/vmOOzh/6dJV6q+YMQOALy5axPfuu2+lumc+4xlc8sIXAvDXt9/O7PvvX6n+WRtvzLf37vxHyqcWLuQnD628w+HkTTflrOc/H4CjFyzg2kceWal+z80359Q99wTgiJtv5uZly1aqn7HllpzwnOcA8N4bbmDx44+vVP+Krbfm7/bYA4A/mD+f+558cqX612+3HX+x++4AHDBvHr9++umV6t/8rGfxiSlTAHjdtdfSnz97/uyBP3v+7Pmz1+bP3uj+7F1x2GFwwgmrjNMrXa1BbmZsn9e87gV+Bvxpk5iurUnAotbx4qasmzaD9d2pqpYANO87DnTyJEckmZNkztIBvlhJkiSNT+lMzg7SIPky8Bbg34GvVdX/tupuqqrnrtWJk3cAb6yqDzbH7wP2raqPttp8H/i7qvpxczwb+HNgj9X1TfJAVW3bGuP+qhp0HXJfX1/NmTNnbS5DkiRJG6gkc6uqr3/5kEssgOuBz1bVsgHq9h1GTIuBKa3jyUD/Of7VtdlkkL53J9mlqpY0yzHuGUaMkiRJGme6WWLxnv7JcTOTS1U9OIxzXwVMTzItySbAwcCsfm1mAYc0u1m8HHiwWTYxWN9ZwKHN50OBi4YRoyRJksaZ1c4gJ9kM2BzYodkqLU3V1sCuwz1xVS1PchRwKTABOL2q5ic5sqk/BbgYOBBYACyjc5Pgavs2Qx8PnJ/kA8AvgXcMN1ZJkiSNH6tdg5zk48DRdJLh9tKHh+jsYPFPIx/e6HANsiRJ0vizxmuQq+orwFeSfLSq/nFEo5MkSZLWE4Mtsfjdqvp34I4kv9+/vqq+M6KRSZIkST0w2C4Wr6WztdtbBqgrwARZkiRJY85gSyz+qnk/fPTCkSRJknpryG3ekvxtkvaDN7ZL8vmRDUuSJEnqjW72QT6gqh5YcVBV99PZek2SJEkac7pJkCck2XTFQZJnApsO0l6SJEnaYHXzqOmzgNlJvk7n5rz3A2eOaFSSJElSjwyZIFfVF5LMA/Zriv66qi4d2bAkSZKk3uhmBhngGmBjOjPI14xcOJIkSVJvdbOLxTuB/wXeDrwT+GmSt490YJIkSVIvdDOD/BngpVV1D0CSicAPgQtGMjBJkiSpF7rZxeIZK5Ljxn1d9pMkSZI2ON3MIP8gyaXAuc3xu4CLRy4kSZIkqXe62cXiz5L8AfAqIMCpVXXhiEcmSZIk9UBXu1hU1beBb49wLJIkSVLPrTZBTvIwnW3dVqkCqqq2HrGoJEmSpB5ZbYJcVVuNZiCSJEnS+qCr3SiS/HaSw5vPOySZNpyTJtk+yeVJbmnet1tNu/2T3JRkQZJjW+V/n+TGJPOSXJhk26Z8apJfJ7m2eZ0ynDglSZI0/nTzoJC/Aj4JfKop2gQ4a5jnPRaYXVXTgdnNcf/zTgBOAg4A9gLenWSvpvpyYJ+qeiFwcys2gF9U1YzmdeQw45QkSdI4080M8tuAtwKPAlTVncBwl1/MBM5sPp8JHDRAm32BBVW1sKqeAM5r+lFVl1XV8qbdlcDkYcYjSZIkAd0lyE9UVdHcsJdki3Vw3p2qaglA877jAG0mAYtax4ubsv7eD1zSOp6W5Jok/5Hk1asLIMkRSeYkmbN06dI1vwJJkiSNSd1s83Z+kn8Btk3yR3QS0q8O1SnJD4GdB6j6TJexZYCylXbVSPIZYDlwdlO0BNitqu5L8lvAvyXZu6oeWmWgqlOBUwH6+voG2q1DkiRJ41A3Dwr5YpI3AA8BzwX+sqou76LffqurS3J3kl2qakmSXYB7Bmi2GJjSOp4M3Nka41DgzcDrmxluqupx4PHm89wkvwD2BOYMFa8kSZIEXSTISf4E+FY3SfEamAUcChzfvF80QJurgOnNjhl3AAcDf9jEtD+dGwdfW1XLWrFOBH5VVU8l2QOYDixch3FLkiRpjOtmDfLWwKVJ/ivJR5LstA7OezzwhiS3AG9ojkmya5KLAZqb8I4CLgVuAM6vqvlN/3+ic6Pg5f22c3sNMC/Jz4ALgCOr6lfrIF5JkiSNE2lWJwzdMHkh8C7gD4DFgy2h2ND09fXVnDmuwpAkSRpPksytqr7+5V09KKRxD3AXcB8D7zohSZIkbfC6eVDIHye5gs4DPXYA/qh5QIckSZI05nSzzdvuwNFVde1IByNJkiT1WjfbvK3yGGhJkiRprFqTNciSJEnSmGeCLEmSJLWYIEuSJEktJsiSJElSiwmyJEmS1GKCLEmSJLWYIEuSJEktJsiSJElSiwmyJEmS1GKCLEmSJLWYIEuSJEktJsiSJElSiwmyJEmS1GKCLEmSJLX0JEFOsn2Sy5Pc0rxvt5p2+ye5KcmCJMe2yo9LckeSa5vXga26TzXtb0ryxtG4HkmSJI0dvZpBPhaYXVXTgdnN8UqSTABOAg4A9gLenWSvVpN/qKoZzevips9ewMHA3sD+wD8340iSJEld6VWCPBM4s/l8JnDQAG32BRZU1cKqegI4r+k31LjnVdXjVXUrsKAZR5IkSepKrxLknapqCUDzvuMAbSYBi1rHi5uyFY5KMi/J6a0lGkP1+Y0kRySZk2TO0qVL1/Y6JEmSNMaMWIKc5IdJrh/gNdQs8G+GGKCsmveTgWcDM4AlwJe66LNyYdWpVdVXVX0TJ07sMiRJkiSNdRuN1MBVtd/q6pLcnWSXqlqSZBfgngGaLQamtI4nA3c2Y9/dGuurwPeG6iNJkiR1o1dLLGYBhzafDwUuGqDNVcD0JNOSbELn5rtZAE1SvcLbgOtb4x6cZNMk04DpwP+OQPySJEkao0ZsBnkIxwPnJ/kA8EvgHQBJdgVOq6oDq2p5kqOAS4EJwOlVNb/p/4UkM+gsn7gN+BBAVc1Pcj7wc2A58JGqemoUr0uSJEkbuFQNuER3XOnr66s5c+b0OgxJkiSNoiRzq6pvlXITZEiyFLh9lE+7A3DvKJ9TI8vvdGzx+xx7/E7HFr/PsaVX3+fuVbXKbg0myD2SZM5A/2LRhsvvdGzx+xx7/E7HFr/PsWV9+z57dZOeJEmStF4yQZYkSZJaTJB759ReB6B1zu90bPH7HHv8TscWv8+xZb36Pl2DLEmSJLU4gyxJkiS1mCBLkiRJLSbIPZBk/yQ3JVmQ5Nhex6O1l2RKkh8luSHJ/CQf73VMWjeSTEhyTZLv9ToWDU+SbZNckOTG5vfqK3odk4YnyZ80f+Zen+TcJJv1OiZ1L8npSe5Jcn2rbPsklye5pXnfrpcxmiCPsiQTgJOAA4C9gHcn2au3UWkYlgPHVNXzgZcDH/H7HDM+DtzQ6yC0TnwF+EFVPQ94EX6vG7Qkk4CPAX1VtQ8wATi4t1FpDZ0B7N+v7FhgdlVNB2Y3xz1jgjz69gUWVNXCqnoCOA+Y2eOYtJaqaklVXd18fpjOX7yTehuVhivJZOBNwGm9jkXDk2Rr4DXA1wCq6omqeqC3UWkd2Ah4ZpKNgM2BO3scj9ZAVf0n8Kt+xTOBM5vPZwIHjWpQ/Zggj75JwKLW8WJMqMaEJFOBFwM/7W0kWgdOAP4ceLrXgWjY9gCWAl9vlsyclmSLXgeltVdVdwBfBH4JLAEerKrLehuV1oGdqmoJdCafgB17GYwJ8ujLAGXutbeBS7Il8G3g6Kp6qNfxaO0leTNwT1XN7XUsWic2Al4CnFxVLwYepcf/davhadamzgSmAbsCWyR5b2+j0lhjgjz6FgNTWseT8b+GNmhJNqaTHJ9dVd/pdTwatlcBb01yG50lUL+b5KzehqRhWAwsrqoV/7NzAZ2EWRuu/YBbq2ppVT0JfAd4ZY9j0vDdnWQXgOb9nl4GY4I8+q4CpieZlmQTOjcWzOpxTFpLSUJnbeMNVfXlXsej4auqT1XV5KqaSuf3579XlbNTG6iqugtYlOS5TdHrgZ/3MCQN3y+BlyfZvPkz+PV44+VYMAs4tPl8KHBRD2Nho16efDyqquVJjgIupXPn7elVNb/HYWntvQp4H3Bdkmubsk9X1cU9jEnSyj4KnN1MSiwEDu9xPBqGqvppkguAq+nsJHQN69ljijW4JOcCrwN2SLIY+CvgeOD8JB+g84+gd/QuQh81LUmSJK3EJRaSJElSiwmyJEmS1GKCLEmSJLWYIEuSJEktJsiSJElSiwmyJEmS1GKCLEljXJKpSX69Yq/u5vj6YYz390nuSvKJdRelJK0/fFCIJI0Pv6iqGetioKr6sySProuxJGl95AyyJG3Akrw0ybwkmyXZIsn8JPusQf89klzTjHNYkn9L8t0ktyY5KsmfNvVXJtl+JK9FktYXJsiStAGrqquAWcDngS8AZ1VVV8snkjwX+DZweDMOwD7AHwL7An8DLKuqFwM/AQ5Zx+FL0nrJJRaStOH7HHAV8BjwsS77TAQuAv6gqua3yn9UVQ8DDyd5EPhuU34d8MJ1FK8krdecQZakDd/2wJbAVsBmXfZ5EFgEvKpf+eOtz0+3jp/GSRVJ44QJsiRt+E4F/gI4G/h/XfZ5AjgIOCTJH45UYJK0IXI2QJI2YEkOAZZX1TlJJgD/k+R3q+rfh+pbVY8meTNwubtSSNL/SVX1OgZJ0ghKMhX4XlV1vbtFF2MeBzxSVV9cV2NK0vrCJRaSNPY9BWyz4kEhw5Xk74H3As46SxqTnEGWJEmSWpxBliRJklpMkCVJkqQWE2RJkiSpxQRZkiRJajFBliRJklpMkCVJkqQWE2RJkiSpxQRZkiRJajFBliRJklpMkCVpA5Xk00lOG+YYU5NUko3WVVyStKHzUdOSNI4lmQrcCmxcVct7G826k6SA6VW1oNexSNrwOIMsSRsgZ3wlaeSYIEvSIJLcluQTSeYleTDJN5Ns1tQdluTH/dpXkuc0n89I8s9JLknySJL/TrJzkhOS3J/kxiQvbvXdNcm3kyxNcmuSj7XqjktyQZKzkjwEHNaUndVq89tJ/ifJA0kWJTmsKX9TkmuSPNSUH9fltR+b5IJ+ZV9JcmLr+hcmebiJ9z2rGecZzVi/SHJfkvOTbN/UvasZY+vm+IAkdyWZ2Drfoib2uUle3Rp3QrPM5BdNDHOTTEnyn02TnzW/7u/q5nolaQUTZEka2juB/YFpwAuBw9aw72eBHYDHgZ8AVzfHFwBfhk4SCXwX+BkwCXg9cHSSN7bGmtn02RY4u32SJLsBlwD/CEwEZgDXNtWPAoc0/d4E/HGSg7qI/VzgwFbyOqG5nnOSbAGcCBxQVVsBr2ydr7+PAQcBrwV2Be4HTgKoqm82vyYnJnkW8DXgg1W1tOl7VXMt2wPnAN9a8Q8U4E+BdwMHAlsD7weWVdVrmvoXVdWWzTkkqWsmyJI0tBOr6s6q+hWdJHbGGvS9sKrmVtVjwIXAY1X1jap6CvgmsGIG+aXAxKr6XFU9UVULga8CB7fG+klV/VtVPV1Vv+53nvcAP6yqc6vqyaq6r6quBaiqK6rquqbfPDqJ72uHCryqbqeTzK9Ipn+XTgJ6ZXP8NLBPkmdW1ZKqmr+aoT4EfKaqFlfV48BxwNtby0Q+0ox9BfDdqvpeK4azmmtZXlVfAjYFnttUfxD4bFXdVB0/q6r7hrouSRqKCbIkDe2u1udlwJZr0Pfu1udfD3C8YqzdgV2b5REPJHkA+DSwU6v9okHOMwX4xUAVSV6W5EfN0o0HgSPpzGB34xw6s7QAf9gcU1WPAu9qxlqS5PtJnreaMXYHLmxd1w3AUzTXVlUPAN8C9gG+1C/2Y5Lc0CxveQDYphX7aq9ZkobDBFmS1t6jwOYrDpLsPIyxFgG3VtW2rddWVXVgq81g2w4tAp69mrpzgFnAlKraBjgFSJdxfQt4XZLJwNuasTrBVF1aVW8AdgFupDPjvbrYDuh3bZtV1R0ASWbQWR5xLp1lGzTlrwY+SWdZx3ZVtS3wYCv2wa5ZktaaCbIkrb2fAXsnmdGsiz1uGGP9L/BQkk8meWZzA9o+SV7aZf+zgf2SvDPJRkme1SSeAFsBv6qqx5LsS2cmuCvNWuArgK/TSeBvAEiyU5K3NmuRHwceoTMrPJBTgL9JsnvTd2KSmc3nzYCz6MyWHw5MSvLhVtzLgaXARkn+ks5a4xVOA/46yfR0vLBZxwydmfo9ur1OSWozQZaktVRVNwOfA34I3AL8ePAeg471FPAWOuubbwXupZMAbtNl/1/SuVntGOBXdG6Ye1FT/WHgc0keBv4SOH8NwzsH2I/W7DGdvz+OAe5szvfa5jwD+QqdGezLmhiuBF7W1P0dsLiqTm7WJ78X+HyS6cCldG48vBm4HXiMlZeZfLm5lsuAh+jc4PfMpu444MxmWcc71/B6JY1zPihEkiRJanEGWZIkSWoxQZYkSZJaTJAlSZKkFhNkSZIkqWWjoZuMfTvssENNnTq112FIkiRpFM2dO/feqprYv9wEGZg6dSpz5szpdRiSJEkaRUluH6i8p0sskuyf5KYkC5IcO0B9kpzY1M9L8pJW3W1JrktybZI5rfLtk1ye5JbmfbvRuh5JkiRt+HqWICeZAJwEHADsBbw7yV79mh0ATG9eRwAn96v/naqaUVV9rfgGzEwAACAASURBVLJjgdlVNR2Y3RxLkiRJXenlDPK+wIKqWlhVTwDnATP7tZkJfKM6rgS2TbLLEOPOBM5sPp8JHLQug5YkSdLY1ssEeRIrPzJ0cVPWbZui89jSuUmOaLXZqaqWADTvOw508iRHJJmTZM7SpUuHcRmSJEkaS3qZIGeAsv7PvR6szauq6iV0lmF8JMlr1uTkVXVqVfVVVd/EiavcvChJkqRxqpcJ8mJgSut4MnBnt22qasX7PcCFdJZsANy9YhlG837POo9ckiRJY1YvE+SrgOlJpiXZBDgYmNWvzSzgkGY3i5cDD1bVkiRbJNkKIMkWwO8B17f6HNp8PhS4aKQvRJIkSWNHz/ZBrqrlSY4CLgUmAKdX1fwkRzb1pwAXAwcCC4BlwOFN952AC5NA5xrOqaofNHXHA+cn+QDwS+Ado3RJkiRJGgNS1X/Z7/jT19dXPihEkiRpfEkyt992wUCPHxQiSZIkrW9MkCVJkqQWE2RJkiSpxQRZkiRJajFBliRJklpMkCVJkqQWE2RJkiSpxQRZkiRJajFBliRJklpMkCVJkqQWE2RJkiSpxQRZkiRJajFBliRJklpMkCVJkqQWE2RJkiSpxQRZkiRJajFBliRJklpMkCVJkqSWnibISfZPclOSBUmOHaA+SU5s6ucleUlTPiXJj5LckGR+ko+3+hyX5I4k1zavA0fzmiRJkrRh26hXJ04yATgJeAOwGLgqyayq+nmr2QHA9Ob1MuDk5n05cExVXZ1kK2Bukstbff+hqr44WtciSZKksWPQBDnJvC7GWFpVr1+Lc+8LLKiqhc25zgNmAu0EeSbwjaoq4Mok2ybZpaqWAEsAqurhJDcAk/r1lSRJktbYUDPIE4DBligEmLWW554ELGodL6YzOzxUm0k0yTFAkqnAi4GfttodleQQYA6dmeb7Vwk8OQI4AmC33XZby0uQJEnSWDPUGuQPVdXtg7xuAz68lufOAGW1Jm2SbAl8Gzi6qh5qik8Gng3MoJNIf2mgk1fVqVXVV1V9EydOXNPYJUmSNEYNmiBX1Y+HGqCbNquxGJjSOp4M3NltmyQb00mOz66q77Tiubuqnqqqp4Gv0lnKIUmSJHWlq10skkxPckGSnydZuOI1zHNfBUxPMi3JJsDBrLpcYxZwSLObxcuBB6tqSZIAXwNuqKov94t1l9bh24DrhxmnJEmSxpFud7H4OvBXwD8AvwMczsDLH7pWVcuTHAVcSmet8+lVNT/JkU39KcDFdNZALwCWNecFeBXwPuC6JNc2ZZ+uqouBLySZQWcpxm3Ah4YTpyRJksaXdDaIGKJRMreqfivJdVX1gqbsv6rq1SMe4Sjo6+urOXPm9DoMSZIkjaImx+3rX97tDPJjSZ4B3NLM+t4B7LguA5QkSZLWB90+Se9oYHPgY8BvAe8FDh2poCRJkqReGXIGOclEOut5N6qqxfzfOmBJkiRpzBl0BjnJB4H5wD8CNyZ566hEJUmSJPXIUDPIRwN7V9XSJHsAZ7P2T86TJEmS1ntDrUF+oqqWAlTVQmDTkQ9JkiRJ6p2hZpAnJzlxdcdV9bGRCUuSJEnqjaES5D/rdzx3pAKRJEmS1geDJshVdeZoBSJJkiStDwZNkJMMekNeVbmrhSRJksaUoZZYvAJYBJwL/BTIiEckSZIk9dBQCfLOwBuAdwN/CHwfOLeq5o90YJIkSVIvDLrNW1U9VVU/qKpDgZcDC4Arknx0VKKTJEmSRlk3j5reFHgTnVnkqcCJwHdGNixJkiSpN4a6Se9MYB/gEuD/q6rrRyUqSZIkqUeGmkF+H/AosCfwseQ39+gFqKraegRjkyRJkkbdUPsgD/UoakmSJGlMMQGWJEmSWgZNkJNcPdQA3bQZpO/+SW5KsiDJsQPUJ8mJTf28JC8Zqm+S7ZNcnuSW5n27tY1PkiRJ489Qa5Cfn2TeIPUBtlmbEyeZAJxEZ5/lxcBVSWZV1c9bzQ4ApjevlwEnAy8bou+xwOyqOr5JnI8FPrk2MUqSJGn8GSpBfl4XYzy1lufeF1hQVQsBkpwHzATaCfJM4BtVVcCVSbZNsgud7eZW13cm8Lqm/5nAFaxvCfLRR/O6M85YpfidEyfy4UmTWPbUUxx43XWr1B+2884ctvPO3Pvkk7x9/qrPavnjXXflXTvuyKLHHuN9N964Sv0xkyfzlh124KZly/jQzTevUv/Z3Xdnv+2249pHHuHoBQtWqf/badN45Tbb8D8PPsinb711lfoTnvMcZmy5JT+8/34+f/vtq9T/y5578tzNN+e7997LlxYvXqX+X5/3PKZsthnfvOceTr7zzlXqL9h7b3bYeGPOuOsuzrjrrlXqL37BC9h8wgT++Y47OH/p0lXqr5gxA4AvLlrE9+67b6W6Zz7jGVzywhcC8Ne3387s++9fqf5ZG2/Mt/feG4BPLVzITx56aKX6yZtuylnPfz4ARy9YwLWPPLJS/Z6bb86pe+4JwBE338zNy5atVD9jyy054TnPAeC9N9zA4scfX6n+FVtvzd/tsQcAfzB/Pvc9+eRK9a/fbjv+YvfdAThg3jx+/fTTK9W/+VnP4hNTpgDwumuvpT9/9vzZA3/2/NnzZ6/Nn73R/dm74rDD4IQTVhmnV4a6SW/VX/F1ZxKdx1ivsJjOLPFQbSYN0XenqloCUFVLkuw40MmTHAEcAbDbbrut5SVIkiRprElncrYHJ07eAbyxqj7YHL8P2LeqPtpq833g76rqx83xbODPgT1W1zfJA1W1bWuM+6tq0HXIfX19NWfOnHV8hZIkSVqfJZlbVX39y3u5i8ViYErreDLQf45/dW0G63t3swyD5v2edRizJEmSxriuEuQkWyR5RvN5zyRvTbLxMM99FTA9ybQkmwAHA7P6tZkFHNLsZvFy4MFm+cRgfWcBhzafDwUuGmackiRJGkeGuklvhf8EXt1smTYbmAO8C3jP2p64qpYnOQq4FJgAnF5V85Mc2dSfAlwMHAgsAJYBhw/Wtxn6eOD8JB8Afgm8Y21jlCRJ0vjT1RrkJFdX1UuSfBR4ZlV9Ick1VfXikQ9x5LkGWZIkafwZ7hrkJHkFnRnj7zdl3c4+S5IkSRuMbhPko4FPARc2yyD2AH40cmFJkiRJvdHVLHBV/QfwHwDNzXr3VtXHRjIwSZIkqRe63cXinCRbJ9mCztPqbkryZyMbmiRJkjT6ul1isVdVPQQcRGdnid2A941YVJIkSVKPdJsgb9zse3wQcFFVPQn05hF8kiRJ0gjqNkH+F+A2YAvgP5PsDjw0UkFJkiRJvdLtTXonAie2im5P8jsjE5IkSZLUO93epLdTkq8luaQ53ov/e5yzJEmSNGZ0u8TiDDqPdd61Ob6Zzt7IkiRJ0pjSbYK8Q1WdDzwNUFXLgadGLCpJkiSpR7pNkB9N8iyanSuSvBx4cMSikiRJknqkq5v0gD8FZgHPTvLfwETg7SMWlSRJktQjQybIzaOlNwNeCzwXCHBTsxeyJEmSNKYMmSBX1dNJvlRVrwDmj0JMkiRJUs90uwb5siR/kCQjGo0kSZLUY2uyBnkLYHmSx+gss6iq2nrEIpMkSZJ6oNsn6W010oFIkiRJ64Nun6Q3u5uybiXZPsnlSW5p3rdbTbv9k9yUZEGSY1vlf5/kxiTzklyYZNumfGqSXye5tnmdsrYxSpIkaXwaNEFOslmS7YEdkmzXJLbbJ5nK/z1Vb20cC8yuqunA7Oa4/7knACcBBwB7Ae9uHnENcDmwT1W9kM5T/T7V6vqLqprRvI4cRoySJEkah4aaQf4QMBd4XvO+4nURneR1bc0Ezmw+nwkcNECbfYEFVbWwqp4Azmv6UVWXNU/zA7gSmDyMWCRJkqTfGDRBrqqvVNU04BNVtUdVTWteL6qqfxrGeXeqqiXNOZYAOw7QZhKwqHW8uCnr7/3AJa3jaUmuSfIfSV49jBglSZI0DnW7i8VdSbaqqoeTfBZ4CfD5qrp6dR2S/BDYeYCqz3R5zoG2lKt+5/gMsBw4uylaAuxWVfcl+S3g35LsXVUPDRDfEcARALvttluXIUmSJGms6zZB/ouq+laS3wbeCHwROBl42eo6VNV+q6tLcneSXapqSZJdgHsGaLYYmNI6ngzc2RrjUODNwOurqppzPg483nyem+QXwJ7AnAHiOxU4FaCvr6/610uSJGl86vZBIU81728CTq6qi4BNhnHeWcChzedD6axp7u8qYHqSaUk2AQ5u+pFkf+CTwFuratmKDkkmNjf3kWQPYDqwcBhxSpIkaZzpNkG+I8m/AO8ELk6y6Rr0HcjxwBuS3AK8oTkmya5JLgZobsI7CrgUuAE4v6pWPOr6n4CtgMv7bef2GmBekp8BFwBHVtWvhhGnJEmSxpk0qxMGb5RsDuwPXFdVtzTLIl5QVZeNdICjoa+vr+bMWWUVhiRJksawJHOrqq9/ebdP0lsGfKd1vITODXGSJEnSmDKcZRKSJEnSmGOCLEmSJLWYIEuSJEktJsiSJElSiwmyJEmS1GKCLEmSJLWYIEuSJEktJsiSJElSiwmyJEmS1GKCLEmSJLWYIEuSJEktJsiSJElSiwmyJEmS1GKCLEmSJLWYIEuSJEktJsiSJElSiwmyJEmS1GKCLEmSJLX0JEFOsn2Sy5Pc0rxvt5p2+ye5KcmCJMe2yo9LckeSa5vXga26TzXtb0ryxtG4HkmSJI0dvZpBPhaYXVXTgdnN8UqSTABOAg4A9gLenWSvVpN/qKoZzevips9ewMHA3sD+wD8340iSJEld6VWCPBM4s/l8JnDQAG32BRZU1cKqegI4r+k31LjnVdXjVXUrsKAZR5IkSepKrxLknapqCUDzvuMAbSYBi1rHi5uyFY5KMi/J6a0lGkP1+Y0kRySZk2TO0qVL1/Y6JEmSNMaMWIKc5IdJrh/gNdQs8G+GGKCsmveTgWcDM4AlwJe66LNyYdWpVdVXVX0TJ07sMiRJkiSNdRuN1MBVtd/q6pLcnWSXqlqSZBfgngGaLQamtI4nA3c2Y9/dGuurwPeG6iNJkiR1o1dLLGYBhzafDwUuGqDNVcD0JNOSbELn5rtZAE1SvcLbgOtb4x6cZNMk04DpwP+OQPySJEkao0ZsBnkIxwPnJ/kA8EvgHQBJdgVOq6oDq2p5kqOAS4EJwOlVNb/p/4UkM+gsn7gN+BBAVc1Pcj7wc2A58JGqemoUr0uSJEkbuFQNuER3XEmyFLh9lE+7A3DvKJ9TI8vvdGzx+xx7/E7HFr/PsaVX3+fuVbXKzWgmyD2SZE5V9fU6Dq07fqdji9/n2ON3Orb4fY4t69v36aOmJUmSpBYTZEmSJKnFBLl3Tu11AFrn/E7HFr/PscfvdGzx+xxb1qvv0zXIkiRJUoszyJIkSVKLCbIkSZLUYoLcA0n2T3JTkgVJju11PFp7SaYk+VGSG5LMT/LxXsekdSPJhCTXJPne0K21PkuybZILktzY/F59Ra9j0vAk+ZPmz9zrk5ybZLNex6TuJTk9yT1Jrm+VbZ/k8iS3NO/b9TJGE+RRlmQCcBJwALAX8O4ke/U2Kg3DcuCYqno+8HLgI36fY8bHgRt6HYTWia8AP6iq5wEvwu91g5ZkEvAxoK+q9qHztN2DexuV1tAZwP79yo4FZlfVdGB2c9wzJsijb19gQVUtrKongPOAmT2OSWupqpZU1dXN54fp/MU7qbdRabiSTAbeBJzW61g0PEm2Bl4DfA2gqp6oqgd6G5XWgY2AZybZCNgcuLPH8WgNVNV/Ar/qVzwTOLP5fCZw0KgG1Y8J8uibBCxqHS/GhGpMSDIVeDHw095GonXgBODPgad7HYiGbQ9gKfD1ZsnMaUm26HVQWntVdQfwReCXwBLgwaq6rLdRaR3YqaqWQGfyCdixl8GYII++DFDmXnsbuCRbAt8Gjq6qh3odj9ZekjcD91TV3F7HonViI+AlwMlV9WLgUXr8X7canmZt6kxgGrArsEWS9/Y2Ko01JsijbzEwpXU8Gf9raIOWZGM6yfHZVfWdXsejYXsV8NYkt9FZAvW7Sc7qbUgahsXA4qpa8T87F9BJmLXh2g+4taqWVtWTwHeAV/Y4Jg3f3Ul2AWje7+llMCbIo+8qYHqSaUk2oXNjwawex6S1lCR01jbeUFVf7nU8Gr6q+lRVTa6qqXR+f/57VTk7tYGqqruARUme2xS9Hvh5D0PS8P0SeHmSzZs/g1+PN16OBbOAQ5vPhwIX9TAWNurlycejqlqe5CjgUjp33p5eVfN7HJbW3quA9wHXJbm2Kft0VV3cw5gkreyjwNnNpMRC4PAex6NhqKqfJrkAuJrOTkLXsJ49pliDS3Iu8DpghySLgb8CjgfOT/IBOv8IekfvIvRR05IkSdJKXGIhSZIktZggS5IkSS0myJIkSVKLCbIkSZLUYoIsSZIktZggS5IkSS0myJI0RiV5XZIHk1zcOv7eMMY7O8mvkrx93UUpSesfE2RJGtv+q6oOXBcDVdV78MmfksYBE2RJGgOSvDTJvCSbJdkiyXxgnyHaX5NkjyTHJTkzyWVJbkvy+0m+kOS6JD9IsvEoXook9ZwJsiSNAVV1FZ3Z3c8DXwDOAq4fqG2SVwKnADOramFT/GzgTcDMpu+PquoFwK+bckkaNzbqdQCSpHXmc8BVwGPAx4BXD9Dm+cCpwO9V1Z2t8kuq6skk1wETgB805dcBU0csYklaDzmDLEljx/bAlsBWwGarabOETgL94n7ljwNU1dPAk1VVTfnTOJkiaZwxQZakseNU4C+As4H/t5o2D9BZMvG3SV43SnFJ0gbFWQFJGgOSHAIsr6pzkkwA/gf4zkBtq+ruJG8BLkny/tGMU5I2BPm//0WTJI0lzQzxJ6rqzetwzDOA71XVBetqTEla37jEQpL+//buPN6uur73/+ttEBxAkBKsJGGQRjHyA9QjcltvxaJlcMBeh8ugKFIxj0Kv1gkc2lpnr1YtD9DcFLlIRbkVQYONIlatFqQlsUwRaSMoCUEIKIOAQsjn98deBxaHfYbknLN39jmv5+ORR/b6ru9a67PX1xPffPd3rzNz3QfsPfyLQiYrydnA8+msYZakGcsZZEmSJKnFGWRJkiSpxYAsSZIktRiQJUmSpBYDsiRJktRiQJYkSZJaDMiSJElSiwFZkiRJajEgS5IkSS0GZEmSJKnFgCxpoCR5X5IvTOL4VUkOnMKSthhJliT5y37XMZqJ3vskleT3elBS+5o/S/LCSZ7jG0leN1U1SeqfrfpdgCRNlyRnAmur6r3DbVX1jP5VNL2qanG/axg22+49QFUdOvw6yeuBP62q5/WvIkmbyxlkSVuMJDP+P9rT8ajx2iZwni32XiWZ0+8aJGkyDMiS+qr5aPukJFcCdyfZKskuSb6SZH2S65P8rzGO/3KSXyS5I8n3kzyjaT8eOBp4Z5JfJ7mgdb0XNte4N8mOrXM9M8mtSR7dbL8hyTVJfpXkwiS7jVHHAUkuSXJ7kivaSwmSfC/Jh5JcDNwDPGWUtl2SLEvyyySrk7yxdY73JTk3yReS3Am8vksNZyb5YPP6wCRrk7wtyS1Jbkpy7Bj1fy/JR5L8e3Mvvzbi3nS9z63rfjbJ8iR3A8eNde+b13OSvDvJT5PclWRlkgVd6tomySeS3JDk5mYZyWNH6Xd7kr1bbXObMd652X5Jksubfpck2WeUe7FNkk8nWdf8+XSSbVr7D2/Oc2dT/yGte/inSZ4OLAH+W/P+b0/ynKb+rVrneUWSy0cbE0n9Y0CWtCU4EngxsAOwEbgAuAKYBxwEvCXJwaMc+w1gIbAz8CPgbICqWtq8/t9VtW1VvbR9UFWtA34IvKLVfBRwblXdn+TlwLuB/wHMBX4AfKlbAUnmAf8EfBDYEXg78JUkc1vdXgscD2wH/HyUti8Ba4FdgFcCH05yUOschwPnNvfp7FHuR9vvAtvTuY/HAacleeIY/Y8B3tBcfwNwSmtf1/vcchTwoea9nMUY977xVjrjfhjwhOa693Tp9zHgqcB+wO817+WvRnaqqt8C5zXnHPZq4F+q6pYkzwLOAN4E/A7wf4Bl7eDb8h7ggOaa+wL7A+8FSLJ/8/7eQWcc/hD42YhargEWAz9s3v8OVXUZcBvwolbX1wD/0OX6kvrMgCxpS3BKVa2pqnuB5wBzq+r9VXVfVV0H/D1wRLcDq+qMqrqrCUjvA/ZNsv0Er/tFmkCVJM01vtjsexPwkaq6pqo2AB8G9htlFvk1wPKqWl5VG6vqImAFnfA37MyqWlVVG6rq/pFtdMLs84CTquo3VXU5cDqdED3sh1X11eYa907g/d0PvL+q7q+q5cCvgaeN0f8fqurqqrob+Evg1WmWS0zgPn+tqi5uavvNBGr7U+C9VXVtdVxRVbe1OzRj8kbgL6rql1V1F51x6Pq/BVrj2TiKh8bzjcD/qap/q6oHqurzwG/pBOGRjqZz326pqvXA3/DQOBwHnFFVFzXv9caq+skE3i/A5+n8b4Vmdv7gVn2StiAGZElbgjWt17sBuzQfS9+e5HY6M7lPGnlQ8zH9R5uPue/koZm8nSZ43XPpfAy+C52ZwKIzUzxcx9+1avglEDozmCPtBrxqRM3PA548ynvs1rYLMBwCh/18xPW6nWMstzXhe9g9wLZj9G+f/+fAo4GdJnifN7W2BcBPx+kzF3gcsLJ1X7/ZtHfzHeCxSZ7b/IfMfsD5zb7dgLeNGKMFdO77SLvw0Cw/zevhfhOpezRfAF6aZFs6s9s/qKqbNvNckqbRFvslD0mzSrVerwGur6qFEzjuKDrLDl5IJ7RtD/yKTpAded5HXrTq9iTfohNWng58qaqGj1kDfKiqJrKUYQ2d2dc3jtGnWy3ttnXAjkm2a4XkXYEbxznHVGqvAd6Vzgz0rYx/n7vVNl6ta4A9gavH6HMrcC/wjKq6cYx+nQtWbUzyj3RmkW8Gvt66l8Pj+aHxzkNnLHYDVjXbuzZt7brHLadLfTcm+SHwJ3RmpD87gfNI6gNnkCVtaf4duDOdL+49tpm93DvJc7r03Y7Ox+S30Zlp/PCI/TcDTxnnel+ks/b2FTz84+4lwLvy0Jf+tk/yqlHOMTwzeHBT72OaL8nNH+faD6qqNcAlwEea4/eh83H+RAL6VHlNkkVJHge8n8567AcY/z53M969Px34QJKF6dgnye+0O1TVRjrLaz7V+qLdvDHWo0NnDP8nnWUS7fH8e2BxM7ucJI9P8uIk23U5x5eA9zZf8tuJzprn4Wdvfw44NslBSR7V1LPXKO9/fpKtR7SfBbwT+P94aHZb0hbGgCxpi9IEspfS+Xj8ejqziKfTmbUc6Sw6H3/fCPwYuHTE/s8Bi5qP1L86yiWX0fny2c1VdUWrjvPpfEHsnGZZwdXAod1O0ITbw+ksBVlPZ5bxHWz6v7FHArvTma08H/jrZj1zr/wDcCbwC+AxwPDTQ8a7z92Md+8/Cfwj8C3gzqb/I55OAZwErAYubcbh24yxjrqq/g24m86SiG+02lfQWYd8Kp3Z79V0eRJI44N01pBfCVxF50uJH2zO8+/AscCngDuAf6Ez2zzSd+jMQP8iya2t9vOb/uc3a70lbYHy0KeJkqTZKsn3gC9U1en9rmWmS/JT4E1V9e1+1yKpO2eQJUnqkSSvoLM++Tv9rkXS6AYqICc5I50H3nf9UkezruyUdB6wf2Xz3EtJkvqumaX/LHBCs75a0hZqoJZYJPlDOs/xPKuq9u6y/zDgz+k8e/S5wN9V1XN7W6UkSZIG2UDNIFfV9+k8i3Q0h9MJz1VVlwI7JHnyGP0lSZKkh5lpz0Gex8MfVr+2aXvEg9iTHE/nV7zy+Mc//tl77dXtKT2SJEmaqVauXHlrVT3ilw/NtICcLm1d15BU1VJgKcDQ0FCtWLFiOuuSJEnSFibJz7u1D9QSiwlYy8N/E9R8HvrtR5IkSdK4ZlpAXgYc0zzN4gDgDn/PvSRJkjbFQC2xSPIl4EBgpyRrgb8GHg1QVUuA5XSeYLEauIfObzuSJEmSJmygAnJVHTnO/gJO6FE5kiRJmoFm2hILSZIkaVIMyJIkSVKLAVmSJElqMSBLkiRJLQZkSZIkqcWALEmSJLUYkCVJkqQWA7IkSZLUYkCWJEmSWgzIkiRJUosBWZIkSWoxIEuSJEktBmRJkiSpxYAsSZIktRiQJUmSpJaBC8hJDklybZLVSU7usn/7JBckuSLJqiTH9qNOSZIkDaaBCshJ5gCnAYcCi4Ajkywa0e0E4MdVtS9wIPC3SbbuaaGSJEkaWAMVkIH9gdVVdV1V3QecAxw+ok8B2yUJsC3wS2BDb8uUJEnSoBq0gDwPWNPaXtu0tZ0KPB1YB1wFvLmqNo48UZLjk6xIsmL9+vXTVa8kSZIGzKAF5HRpqxHbBwOXA7sA+wGnJnnCIw6qWlpVQ1U1NHfu3KmvVJIkSQNp0ALyWmBBa3s+nZnitmOB86pjNXA9sFeP6pMkSdKAG7SAfBmwMMkezRfvjgCWjehzA3AQQJInAU8DrutplZIkSRpYW/W7gE1RVRuSnAhcCMwBzqiqVUkWN/uXAB8AzkxyFZ0lGSdV1a19K1qSJEkDZaACMkBVLQeWj2hb0nq9DvjjXtclSZKkmWHQllhIkiRJ08qALEmSJLUYkCVJkqQWA7IkSZLUYkCWJEmSWgzIkiRJUosBWZIkSWoxIEuSJEktBmRJkiSpxYAsSZIktRiQJUmSpBYDsiRJktRiQJYkSZJaDMiSJElSiwFZkiRJahm4gJzkkCTXJlmd5ORR+hyY5PIkq5L8S69rlCRJ0uDaqt8FbIokc4DTgBcBa4HLkiyrqh+3+uwAfAY4pKpuSLJzf6qVJEnSIBq0GeT9gdVVdV1V3QecAxw+os9RwHlVdQNAVd3S4xolSZI0wAYtIM8D1rS21zZtbU8Fnpjke0lWJjmm24mSHJ9kRZIV69evn6ZyqTaiEQAAIABJREFUJUmSNGgGLSCnS1uN2N4KeDbwYuBg4C+TPPURB1UtraqhqhqaO3fu1FcqSZKkgdSXgJxkTpK/2IxD1wILWtvzgXVd+nyzqu6uqluB7wP7bl6lkiRJmm36EpCr6gEeuXZ4Ii4DFibZI8nWwBHAshF9vgb89yRbJXkc8FzgmkkVLEmSpFmjn0+xuDjJqcD/A+4ebqyqH412QFVtSHIicCEwBzijqlYlWdzsX1JV1yT5JnAlsBE4vaquns43IkmSpJkjVSOX8Pbowsl3uzRXVf1Rr2sZGhqqFStW9PqykiRJ6qMkK6tqaGR732aQq+oF/bq2JEmSNJq+PcUiyfZJPjn8qLUkf5tk+37VI0mSJEF/H/N2BnAX8Ormz53A/+1jPZIkSVJfv6S3Z1W9orX9N0ku71s1kiRJEv2dQb43yfOGN5L8AXBvH+uRJEmS+jqDvBg4q7Xu+FfA6/pYjyRJktSfgJzkUcDTqmrfJE8AqKo7+1GLJEmS1Nav36S3ETixeX2n4ViSJElbin6uQb4oyduTLEiy4/CfPtYjSZIk9XUN8huav09otRXwlD7UIkmSJAH9XYP8mqq6uB/XlyRJkkbTzzXIn+jHtSVJkqSx9HMN8reSvCJJ+liDJEmS9DD9XIP8VuBxwANJfgMEqKp6Qh9rkiRJ0izXz4C8PXA0sEdVvT/JrsCT+1iPJEmS1NclFqcBBwBHNtt3AaeOd1CSQ5Jcm2R1kpPH6PecJA8keeXUlCtJkqTZoJ8B+blVdQLwG4Cq+hWw9VgHJJlDJ1gfCiwCjkyyaJR+HwMunOqiJUmSNLP1MyDf3wTZAkgyF9g4zjH7A6ur6rqqug84Bzi8S78/B74C3DKF9UqSJGkW6GdAPgU4H9g5yYeAfwU+PM4x84A1re21TduDkswD/gRYMtaJkhyfZEWSFevXr9/U2iVJkjRD9e1LelV1dpKVwEF0nmDx8qq6ZpzDuj0SrkZsfxo4qaoeGOsJclW1FFgKMDQ0NPIckiRJmqX6+RQLquonwE824ZC1wILW9nxg3Yg+Q8A5TTjeCTgsyYaq+upkapUkSdLs0NeAvBkuAxYm2QO4ETgCOKrdoar2GH6d5Ezg64ZjSZIkTdRABeSq2pDkRDpPp5gDnFFVq5IsbvaPue5YkiRJGs9ABWSAqloOLB/R1jUYV9Xre1GTJEmSZo5+PsVCkiRJ2uIYkCVJkqQWA7IkSZLUYkCWJEmSWgzIkiRJUosBWZIkSWoxIEuSJEktBmRJkiSpxYAsSZIktRiQJUmSpBYDsiRJktRiQJYkSZJaDMiSJElSiwFZkiRJahm4gJzkkCTXJlmd5OQu+49OcmXz55Ik+/ajTkmSJA2mgQrISeYApwGHAouAI5MsGtHteuD5VbUP8AFgaW+rlCRJ0iAbqIAM7A+srqrrquo+4Bzg8HaHqrqkqn7VbF4KzO9xjZIkSRpggxaQ5wFrWttrm7bRHAd8Y1orkiRJ0oyyVb8L2ETp0lZdOyYvoBOQnzfK/uOB4wF23XXXqapPkiRJA27QZpDXAgta2/OBdSM7JdkHOB04vKpu63aiqlpaVUNVNTR37txpKVaSJEmDZ9AC8mXAwiR7JNkaOAJY1u6QZFfgPOC1VfWffahRkiRJA2ygllhU1YYkJwIXAnOAM6pqVZLFzf4lwF8BvwN8JgnAhqoa6lfNkiRJGiyp6rqEd1YZGhqqFStW9LsMSZIk9VCSld0mUgdtiYUkSZI0rQzIkiRJUosBWZIkSWoxIEuSJEktBmRJkiSpxYAsSZIktRiQJUmSpBYDsiRJktRiQJYkSZJaDMiSJElSiwFZkiRJajEgS5IkSS0GZEmSJKnFgCxJkiS1GJAlSZKkFgOyJEmS1GJAliRJkloMyJIkSVKLAVmSJElqSVX1u4a+S7Ie+Hm/65jBdgJu7XcRmjKO58zieM48junM4nhOr92qau7IRgOypl2SFVU11O86NDUcz5nF8Zx5HNOZxfHsD5dYSJIkSS0GZEmSJKnFgKxeWNrvAjSlHM+ZxfGceRzTmcXx7APXIEuSJEktziBLkiRJLQZkTVqSHZNclOS/mr+fOEq/Q5Jcm2R1kpO77H97kkqy0/RXrbFMdkyTfDzJT5JcmeT8JDv0rnoNm8DPXJKc0uy/MsmzJnqsem9zxzPJgiTfTXJNklVJ3tz76tXNZH5Gm/1zkvxHkq/3rurZwYCsqXAy8M9VtRD452b7YZLMAU4DDgUWAUcmWdTavwB4EXBDTyrWeCY7phcBe1fVPsB/Au/qSdV60Hg/c41DgYXNn+OBz27CseqhyYwnsAF4W1U9HTgAOMHx7L9JjumwNwPXTHOps5IBWVPhcODzzevPAy/v0md/YHVVXVdV9wHnNMcN+xTwTsBF8VuGSY1pVX2rqjY0/S4F5k9zvXqk8X7maLbPqo5LgR2SPHmCx6q3Nns8q+qmqvoRQFXdRSdQzetl8epqMj+jJJkPvBg4vZdFzxYGZE2FJ1XVTQDN3zt36TMPWNPaXtu0keRlwI1VdcV0F6oJm9SYjvAG4BtTXqHGM5HxGa3PRMdWvTOZ8XxQkt2BZwL/NuUValNNdkw/TWdiaeN0FTibbdXvAjQYknwb+N0uu94z0VN0aaskj2vO8cebW5s2z3SN6YhrvIfOx7tnb1p1mgLjjs8YfSZyrHprMuPZ2ZlsC3wFeEtV3TmFtWnzbPaYJnkJcEtVrUxy4JRXJgOyJqaqXjjaviQ3D3+M13z0c0uXbmuBBa3t+cA6YE9gD+CKJMPtP0qyf1X9YsregB5hGsd0+ByvA14CHFQ+T7IfxhyfcfpsPYFj1VuTGU+SPJpOOD67qs6bxjo1cZMZ01cCL0tyGPAY4AlJvlBVr5nGemcVl1hoKiwDXte8fh3wtS59LgMWJtkjydbAEcCyqrqqqnauqt2ranc6/xg8y3Dcd5s9ptD5ZjZwEvCyqrqnB/XqkUYdn5ZlwDHNN+UPAO5oltRM5Fj11maPZzqzD58DrqmqT/a2bI1hs8e0qt5VVfOb/988AviO4XhqOYOsqfBR4B+THEfnKRSvAkiyC3B6VR1WVRuSnAhcCMwBzqiqVX2rWOOZ7JieCmwDXNR8MnBpVS3u9ZuYzUYbnySLm/1LgOXAYcBq4B7g2LGO7cPbUGMy4wn8AfBa4Koklzdt766q5b18D3q4SY6pppm/SU+SJElqcYmFJEmS1GJAliRJkloMyJIkSVKLAVmSJElqMSBLkiRJLQZkSZIkqcWALEkDJskOSf6stb1LknOn4TrvS3JjkveP0WfPJJcn+fVUX1+S+sXnIEvSgEmyO/D1qtp7mq/zPuDXVfWJCfT9dVVtO531SFKvOIMsSYPno8DwzO3Hk+ye5GqAJK9P8tUkFyS5PsmJSd6a5D+SXJpkx6bfnkm+mWRlkh8k2Wu8iyZ5fnPNy5vzbTfN71OS+sJfNS1Jg+dkYO+q2g8enFFu2xt4JvAYOr+i9qSqemaSTwHHAJ8GlgKLq+q/kjwX+AzwR+Nc9+3ACVV1cZJtgd9M0fuRpC2KAVmSZp7vVtVdwF1J7gAuaNqvAvZpwu3vA19OMnzMNhM478XAJ5OcDZxXVWunuG5J2iIYkCVp5vlt6/XG1vZGOv/uPwq4fXgGeqKq6qNJ/gk4DLg0yQur6idTUbAkbUlcgyxJg+cuYLPX/1bVncD1SV4FkI59xzsuyZ5VdVVVfQxYAYy7blmSBpEBWZIGTFXdBlyc5OokH9/M0xwNHJfkCmAVcPgEjnlLc80rgHuBb2zmtSVpi+Zj3iRJXfmYN0mzlTPIkqTR/Bo4fiK/KAS4uXdlSdL0cgZZkiRJanEGWZIkSWoxIEuSJEktBmRJkiSpxYAsSZIktRiQJUmSpBYDsiRJktRiQJYkSZJaDMiSJElSiwFZkiRJajEgS9IEJHlfki9M4vhVSQ6cwpK2GEmWJPnLftchSVNlq34XIEkzTZIzgbVV9d7htqp6Rv8qml5VtXiqz5nk9cCfVtXzpvrckjQeZ5AlzXpJZvxkQToeNV7bBM6zxdyrJHP6XYOkmcmALGlWSvKzJCcluRK4O8lWSXZJ8pUk65Ncn+R/jXH8l5P8IskdSb6f5BlN+/HA0cA7k/w6yQWt672wuca9SXZsneuZSW5N8uhm+w1JrknyqyQXJtltjDoOSHJJktuTXNFexpHke0k+lORi4B7gKaO07ZJkWZJfJlmd5I2tc7wvyblJvpDkTuD1XWo4M8kHm9cHJlmb5G1JbklyU5Jjx6j/9UmuS3JXc8+PTvJ0YAnw35p7eHvrOp9NsjzJ3cALkmyT5BNJbkhyc7Pc47FN/52SfL25N79M8oPh/yBoxv7G5rrXJjlotBolzT4GZEmz2ZHAi4EdgI3ABcAVwDzgIOAtSQ4e5dhvAAuBnYEfAWcDVNXS5vX/rqptq+ql7YOqah3wQ+AVreajgHOr6v4kLwfeDfwPYC7wA+BL3QpIMg/4J+CDwI7A24GvJJnb6vZa4HhgO+Dno7R9CVgL7AK8EvjwiMB4OHBuc5/OHuV+tP0usD2d+3gccFqSJ3ap//HAKcChVbUd8PvA5VV1DbAY+GFzD3doHXYU8KGm9n8FPgY8FdgP+L3mmn/V9H1b877mAk+ic18rydOAE4HnNNc9GPjZBN6XpFnCgCxpNjulqtZU1b3Ac4C5VfX+qrqvqq4D/h44otuBVXVGVd1VVb8F3gfsm2T7CV73i3TCOUnSXOOLzb43AR+pqmuqagPwYWC/UWaRXwMsr6rlVbWxqi4CVgCHtfqcWVWrqmpDVd0/so1OmH0ecFJV/aaqLgdOpxOih/2wqr7aXOPeCby/+4H3V9X9VbUc+DXwtFH6bgT2TvLYqrqpqlaNc+6vVdXFVbUR+C3wRuAvquqXVXUXnfs1PGb3A08Gdmtq+UFVFfAAsA2wKMmjq+pnVfXTCbwvSbOEAVnSbLam9Xo3YJfm4/jbm4/1301n5vFhksxJ8tEkP22WHfys2bXTBK97Lp3lA7sAfwgUnZni4Tr+rlXDL4HQmRkdaTfgVSNqfh6dUNjtPXZr2wUYDpfDfj7iet3OMZbbmvA97B5g25Gdqupu4H/SmS2+Kck/JdlrnHO3a5kLPA5Y2Xr/32zaAT4OrAa+1SzjOLm57mrgLXT+w+aWJOc0YyFJgAFZ0uxWrddrgOuraofWn+2q6rAuxx1FZ9nBC+ksJdi9aU+X8z7yolW3A98CXt2c60vNzOZwHW8aUcdjq+qSLqdaA/zDiL6Pr6qPjvIeu7WtA3ZMsl2rbVfgxnHOMSWq6sKqehGdUP8TOrP2Y12z3X4rcC/wjNb7376qtm3OfVdVva2qngK8FHjr8NKRqvpi84SM3ZpzfmzK35ykgWVAlqSOfwfubL689dhmlnjvJM/p0nc7Oh/v30ZnBvPDI/bfDDxlnOt9ETiGzlrkL7balwDvan3pb/skrxrlHF8AXprk4KbexzRfkps/zrUfVFVrgEuAjzTH70Nn3fBE1hpPSpInJXlZsxb5t3SWYjzQ7L4ZmJ9k69GOb5ZZ/D3wqSQ7N+ecN7xuPMlLkvxes4zlzubcDyR5WpI/SrIN8Bs6IfuB7leRNBsZkCUJqKoH6Mwy7gdcT2d28nQ6M8QjnUVnGcKNwI+BS0fs/xyd9a23J/nqKJdcRudLfjdX1RWtOs6nM5t5TrN842rg0FFqXkNnJvvdwHo6M8rvYNP/bT+Sziz4OuB84K+b9czT7VF0vki3js5SkucDf9bs+w6wCvhFklvHOMdJdJZRXNrcr2/z0Hrnhc32r+l8MfIzVfU9OuuPP0pnjH9B54uW756ydyVp4OWhT/UkSZIkOYMsSZIktQxUQE5yRvPg+atH2Z8kpzQPur8yybN6XaMkSZIG20AFZOBM4JAx9h9KZ83ZQjoPwf9sD2qSJEnSDDJQAbmqvk/nixyjORw4qzouBXZI8uQx+kuSJEkPs1W/C5hi83j4Q+TXNm03jeyY5Hg6s8w8/vGPf/Zee433bHpJkiTNJCtXrry1quaObJ9pATld2ro+pqOqlgJLAYaGhmrFihXTWZckSZK2MEl+3q19oJZYTMBaYEFrez6d52tKkiRJEzLTAvIy4JjmaRYHAHdU1SOWV0iSJEmjGaglFkm+BBwI7JRkLfDXwKMBqmoJsBw4jM5vVboHOLY/lUqSJGlQDVRArqojx9lfwAk9KkeSJEkz0ExbYiFJkiRNigFZkiRJajEgS5IkSS0GZEmSJKnFgCxJkiS1GJAlSZKkFgOyJEmS1GJAliRJkloMyJIkSVKLAVmSJElqMSBLkiRJLQZkSZIkqcWALEmSJLUYkCVJkqQWA7IkSZLUMnABOckhSa5NsjrJyV32b5/kgiRXJFmV5Nh+1ClJkqTBNFABOckc4DTgUGARcGSSRSO6nQD8uKr2BQ4E/jbJ1j0tVJIkSQNroAIysD+wuqquq6r7gHOAw0f0KWC7JAG2BX4JbOhtmZIkSRpUgxaQ5wFrWttrm7a2U4GnA+uAq4A3V9XGkSdKcnySFUlWrF+/frrqlSRJ0oAZtICcLm01Yvtg4HJgF2A/4NQkT3jEQVVLq2qoqobmzp079ZVKkiRpIA1aQF4LLGhtz6czU9x2LHBedawGrgf26lF9kiRJGnCDFpAvAxYm2aP54t0RwLIRfW4ADgJI8iTgacB1Pa1SkiRJA2urfhewKapqQ5ITgQuBOcAZVbUqyeJm/xLgA8CZSa6isyTjpKq6tW9FS5IkaaAMVEAGqKrlwPIRbUtar9cBf9zruiRJkjQzDNoSC0mSJGlaGZAlSZKkFgOyJEmS1GJAliRJkloMyJIkSVKLAVmSJElqMSBLkiRJLQZkSZIkqcWALEmSJLUYkCVJkqQWA7IkSZLUYkCWJEmSWgzIkiRJUosBWZIkSWoxIEuSJEktAxeQkxyS5Nokq5OcPEqfA5NcnmRVkn/pdY2SJEkaXFv1u4BNkWQOcBrwImAtcFmSZVX141afHYDPAIdU1Q1Jdu5PtZIkSRpEgzaDvD+wuqquq6r7gHOAw0f0OQo4r6puAKiqW3pcoyRJkgbYoAXkecCa1vbapq3tqcATk3wvycokx3Q7UZLjk6xIsmL9+vXTVK4kSZIGzaAF5HRpqxHbWwHPBl4MHAz8ZZKnPuKgqqVVNVRVQ3Pnzp36SiVJkjSQ+hKQk8xJ8hebcehaYEFrez6wrkufb1bV3VV1K/B9YN/Nq1SSJEmzTV8CclU9wCPXDk/EZcDCJHsk2Ro4Alg2os/XgP+eZKskjwOeC1wzqYIlSZI0a/TzKRYXJzkV+H/A3cONVfWj0Q6oqg1JTgQuBOYAZ1TVqiSLm/1LquqaJN8ErgQ2AqdX1dXT+UYkSZI0c6Rq5BLeHl04+W6X5qqqP+p1LUNDQ7VixYpeX1aSJEl9lGRlVQ2NbO/bDHJVvaBf15YkSZJG07enWCTZPsknhx+1luRvk2zfr3okSZIk6O9j3s4A7gJe3fy5E/i/faxHkiRJ6uuX9Pasqle0tv8myeV9q0aSJEmivzPI9yZ53vBGkj8A7u1jPZIkSVJfZ5AXA2e11h3/CnhdH+uRJEmS+hOQkzwKeFpV7ZvkCQBVdWc/apEkSZLa+vWb9DYCJzav7zQcS5IkaUvRzzXIFyV5e5IFSXYc/tPHeiRJkqS+rkF+Q/P3Ca22Ap7Sh1okSZIkoL9rkF9TVRf34/qSJEnSaPq5BvkT/bi2JEmSNJZ+rkH+VpJXJEkfa5AkSZIepp9rkN8KPA54IMlvgABVVU/oY02SJEma5foZkLcHjgb2qKr3J9kVeHIf65EkSZL6usTiNOAA4Mhm+y7g1PEOSnJIkmuTrE5y8hj9npPkgSSvnJpyJUmSNBv0MyA/t6pOAH4DUFW/ArYe64Akc+gE60OBRcCRSRaN0u9jwIVTXbQkSZJmtn4G5PubIFsASeYCG8c5Zn9gdVVdV1X3AecAh3fp9+fAV4BbprBeSZIkzQL9DMinAOcDOyf5EPCvwIfHOWYesKa1vbZpe1CSecCfAEvGOlGS45OsSLJi/fr1m1q7JEmSZqi+fUmvqs5OshI4iM4TLF5eVdeMc1i3R8LViO1PAydV1QNjPUGuqpYCSwGGhoZGnkOSJEmzVD+fYkFV/QT4ySYcshZY0NqeD6wb0WcIOKcJxzsBhyXZUFVfnUytkiRJmh36GpA3w2XAwiR7ADcCRwBHtTtU1R7Dr5OcCXzdcCxJkqSJGqiAXFUbkpxI5+kUc4AzqmpVksXN/jHXHUuSJEnjGaiADFBVy4HlI9q6BuOqen0vapIkSdLM0c+nWEiSJElbHAOyJEmS1GJAliRJkloMyJIkSVKLAVmSJElqMSBLkiRJLQZkSZIkqcWALEmSJLUYkCVJkqQWA7IkSZLUYkCWJEmSWgzIkiRJUosBWZIkSWoxIEuSJEktBmRJkiSpZeACcpJDklybZHWSk7vsPzrJlc2fS5Ls2486JUmSNJgGKiAnmQOcBhwKLAKOTLJoRLfrgedX1T7AB4Clva1SkiRJg2ygAjKwP7C6qq6rqvuAc4DD2x2q6pKq+lWzeSkwv8c1SpIkaYANWkCeB6xpba9t2kZzHPCNbjuSHJ9kRZIV69evn8ISJUmSNMgGLSCnS1t17Zi8gE5APqnb/qpaWlVDVTU0d+7cKSxRkiRJg2yrfhewidYCC1rb84F1Izsl2Qc4HTi0qm7rUW2SJEmaAQZtBvkyYGGSPZJsDRwBLGt3SLIrcB7w2qr6zz7UKEmSpAE2UDPIVbUhyYnAhcAc4IyqWpVkcbN/CfBXwO8An0kCsKGqhvpVsyRJkgZLqrou4Z1VhoaGasWKFf0uQ5IkST2UZGW3idRBW2IhSZIkTSsDsiRJktRiQJYkSZJaDMiSJElSiwFZkiRJajEgS5IkSS0GZEmSJKnFgCxJkiS1GJAlSZKkFgOyJEmS1GJAliRJkloMyJIkSVKLAVmSJElqMSBLkiRJLQZkSZIkqcWALEmSJLUYkCVJkqQWA7IkSZLUkqrqdw19l2Q98PN+1zGD7QTc2u8iNGUcz5nF8Zx5HNOZxfGcXrtV1dyRjQZkTbskK6pqqN91aGo4njOL4znzOKYzi+PZHy6xkCRJkloMyJIkSVKLAVm9sLTfBWhKOZ4zi+M58zimM4vj2QeuQZYkSZJanEGWJEmSWgzIkiRJUosBWZOWZMckFyX5r+bvJ47S75Ak1yZZneTkLvvfnqSS7DT9VWsskx3TJB9P8pMkVyY5P8kOvatewybwM5ckpzT7r0zyrIkeq97b3PFMsiDJd5Nck2RVkjf3vnp1M5mf0Wb/nCT/keTrvat6djAgayqcDPxzVS0E/rnZfpgkc4DTgEOBRcCRSRa19i8AXgTc0JOKNZ7JjulFwN5VtQ/wn8C7elK1HjTez1zjUGBh8+d44LObcKx6aDLjCWwA3lZVTwcOAE5wPPtvkmM67M3ANdNc6qxkQNZUOBz4fPP688DLu/TZH1hdVddV1X3AOc1xwz4FvBPwW6NbhkmNaVV9q6o2NP0uBeZPc716pPF+5mi2z6qOS4Edkjx5gseqtzZ7PKvqpqr6EUBV3UUnUM3rZfHqajI/oySZD7wYOL2XRc8WBmRNhSdV1U0Azd87d+kzD1jT2l7btJHkZcCNVXXFdBeqCZvUmI7wBuAbU16hxjOR8Rmtz0THVr0zmfF8UJLdgWcC/zblFWpTTXZMP01nYmnjdBU4m23V7wI0GJJ8G/jdLrveM9FTdGmrJI9rzvHHm1ubNs90jemIa7yHzse7Z29adZoC447PGH0mcqx6azLj2dmZbAt8BXhLVd05hbVp82z2mCZ5CXBLVa1McuCUVyYDsiamql442r4kNw9/jNd89HNLl25rgQWt7fnAOmBPYA/giiTD7T9Ksn9V/WLK3oAeYRrHdPgcrwNeAhxUPnC9H8Ycn3H6bD2BY9VbkxlPkjyaTjg+u6rOm8Y6NXGTGdNXAi9LchjwGOAJSb5QVa+ZxnpnFZdYaCosA17XvH4d8LUufS4DFibZI8nWwBHAsqq6qqp2rqrdq2p3Ov8YPMtw3HebPabQ+WY2cBLwsqq6pwf16pFGHZ+WZcAxzTflDwDuaJbUTORY9dZmj2c6sw+fA66pqk/2tmyNYbPHtKreVVXzm//fPAL4juF4ajmDrKnwUeAfkxxH5ykUrwJIsgtwelUdVlUbkpwIXAjMAc6oqlV9q1jjmeyYngpsA1zUfDJwaVUt7vWbmM1GG58ki5v9S4DlwGHAauAe4Nixju3D21BjMuMJ/AHwWuCqJJc3be+uquW9fA96uEmOqaaZv2pakiRJanGJhSRJktRiQJYkSZJaDMiSJElSiwFZkiRJajEgS5IkSS0GZEmSJKnFgCxJAyzJDkn+rHm9S5Jzp+i8Zya5fviZrElenmRRa//Hk/wiydun4nqStCUxIEvSYNsB+DOAqlpXVa+cwnO/o/llBQAvBx4MyFX1DmBJ16MkacAZkCVpsH0U2DPJ5Um+nORqgCSvT/LVJBc0M8EnJnlrkv9IcmmSHZt+eyb5ZpKVSX6QZK+RF0jy+8DLgI8319mzp+9QknrMgCxJg+1k4KdVtR/wjhH79gaOAvYHPgTcU1XPBH4IHNP0WQr8eVU9G3g78JmRF6iqS4BldGaU96uqn07LO5GkLcRW/S5AkjRtvltVdwF3JbkDuKBpvwrYJ8m2wO8DX04yfMw2vS9TkrYsBmRJmrl+23q9sbW9kc6//48Cbm9mnyVJDZdYSNJguwvYbnMOrKo7geuTvAogHftO9XUkadAYkCVpgFXVbcDFzZfzPr4ZpzgaOC7JFcAq4PBR+p0DvKP5kp9f0pM0o6Wq+l2DJGkLk+RM4OtVNepzlZO8D/h1VX2iV3VJUi84gyxJ6uYO4APDvyhkpCQfB14D3N3TqiSpB5xBliRJklqcQZYkSZJaDMhI1RjzAAAAGElEQVSSJElSiwFZkiRJajEgS5IkSS3/P0lada+ULNXpAAAAAElFTkSuQmCC\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 =  9.381083924999999\n",
      "spatial order  of accuracy =  6\n",
      "number of grid points =  501\n",
      "maximum relative error in particle velocity =  0.0007690689466005808\n",
      "maximum relative error in stress =  0.0007652270810542603\n"
     ]
    }
   ],
   "source": [
    "# Computation and plotting\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, v, '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, s, '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 [s]')\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",
    "\n",
    "t=0   # initial time\n",
    "\n",
    "forcing = 1.0  # forcing function, forcing = 1,  and no forcing function, forcing = 0\n",
    "\n",
    "# type of initial data: Gaussian or Sinusoidal\n",
    "type_0 = 'Gaussian'\n",
    "#type_0 = 'Sinusoidal'\n",
    "\n",
    "\n",
    "if type_0 in ('Sinusoidal'):\n",
    "            forcing = 1.0  # we must use forcing for Sinusoidal initial condition\n",
    "\n",
    "# L2-norm normalizer\n",
    "# Generate  conditions for normalization\n",
    "rate.mms(v, s, U_t, V_t, U_x, V_x, y, 0.65, type_0)\n",
    "A =  (np.linalg.norm(v)) \n",
    "B =  (np.linalg.norm(s))\n",
    "\n",
    "\n",
    "# Loop through time and evolve the wave-fields using ADER time-stepping scheme of N+1 order of accuracy\n",
    "start = timeit.default_timer()\n",
    "\n",
    "# Generate initial conditions\n",
    "rate.mms(v, s, U_t, V_t, U_x, V_x, y, t, type_0)\n",
    "\n",
    "for t in utils.drange (0.0, tend+dt,dt):\n",
    "    n = n+1\n",
    "    \n",
    "    # compute numerical solution \n",
    "    time_integrator.elastic_RK4(v, s, v, s, rho, mu, nx, dx, order, y, t, dt, r0, r1,  tau_11,\\\n",
    "                    tau_21, tau_12, tau_22, type_0, forcing)\n",
    " \n",
    "    # Analytical solution\n",
    "    rate.mms(U, V, U_t, V_t, U_x, V_x, y, t+dt, type_0)\n",
    "    \n",
    "    # compute error and append to the error array\n",
    "    EU.append(np.linalg.norm(U-v)/A)\n",
    "    EV.append(np.linalg.norm(V-s)/B)\n",
    "    \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, v, 'r', y, U, 'k--')\n",
    "        ax1.legend(iter(line1),('Numerical', 'Analytical'))\n",
    "        line2 = ax2.plot(y, s, '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('spatial order  of accuracy = ', order)                                  # print the polynomial degree used\n",
    "print('number of grid points = ', nx)                     # 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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}