Computational_Seismology_Igel/dg_elastic_nonlinearfriction.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 for Earthquake Rupture Dynamics - non-linear friction laws </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 physicaly motivaled 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",
    "\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",
    "To begin, consider the domain $\\Omega = \\Omega_{-}\\cup  \\Omega_{+}$, with $  \\Omega_{-}:= [0, x_0]$,  $\\Omega_{+}:= [x_0, L]$, $0<x_0<  L$. \n",
    "\n",
    "The sub-domains are governed by the source-free elastic wave equation in a heterogeneous 1D medium is \n",
    "\n",
    "\\begin{align}\n",
    "\\rho^{\\pm}(x)\\partial_t v^{\\pm}(x,t) -\\partial_x \\sigma^{\\pm}(x,t) & = 0\\\\\n",
    "\\frac{1}{\\mu^{\\pm}(x)}\\partial_t \\sigma^{\\pm}(x,t) -\\partial_x v^{\\pm}(x,t) & = 0  \n",
    "\\end{align}\n",
    "\n",
    "Here, the field variables and material parameters in the sub-domains  $  \\Omega_{\\pm}$ with the superscripts $\\pm$:  $v^{\\pm}$,  $\\sigma^{\\pm}$, $\\rho^{\\pm}$,  $\\mu^{\\pm}$, $Z_s^{\\pm}$.\n",
    "\n",
    "The interface at $x = x_0$ is governed by the general non-linear frictional condition. \n",
    "\n",
    "\\begin{align}\n",
    "\\text{force balance}:  \\quad &\\sigma^{-} = \\sigma^{+} = \\sigma, \\nonumber \\\\\n",
    "\\text{fricition law}: \\quad &\\sigma =  \\sigma_{n}\\frac{f\\left(\\left|[\\![ v ]\\!]\\right|,\\psi\\right)}{\\left|[\\![ v ]\\!]\\right|}[\\![ v ]\\!]\n",
    "\\end{align}\n",
    "\n",
    "Here, $\\sigma_{n} > 0$ is the effective normal stress,  $f\\left(\\left|[\\![ v ]\\!]\\right|,\\psi\\right)$ is a non-linear friction coefficient modeling the frictional response of the interface. In general $f\\left(\\left|[\\![ v ]\\!]\\right|,\\psi\\right)\\ge 0$ is a monotically increasing positive function, with $f\\left(0, \\psi\\right) = 0 $. Note that with $f\\left(\\left|[\\![ v ]\\!]\\right|,\\psi\\right)$ we can describe both rate-and-state, and slip-weakening friction laws.\n",
    "\n",
    "At the external 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 define the mechanical energy in each subdomain by\n",
    "\n",
    "\\begin{equation}\n",
    "E^{\\pm}(t) = \\frac{1}{2}\\int_{\\Omega_{\\pm}}{\\left({\\rho^{\\pm}(x)} |v^{\\pm}(x, t)|^2 + \\frac{1}{\\mu^{\\pm}(x)}|\\sigma^{\\pm}(x, t)|^2\\right) dx}.\n",
    "\\end{equation}\n",
    "\n",
    "The elastic wave equation with the frictional interface condition, satisfies the energy balance\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{d \\left(E^-(t)+E^+(t)\\right)}{dt}  = -\\sigma [\\![ v]\\!]  -v^-(0, t)\\sigma^-(0, t) + v^+(L, t)\\sigma^+(L, t),\\le 0.\n",
    "\\end{equation}\n",
    "\n",
    "With the above friction law the first term in the right hand side of the energy rate, $\\sigma [\\![ v]\\!] = f\\left(\\left|[\\![ v ]\\!]\\right|,\\psi\\right){\\left|[\\![ v ]\\!]\\right|} \\ge 0$, is the work done by friction on the fault which is dissipated as heat. \n",
    "\n",
    "### The discontinuous Galerkin spectral element method (DGSEM)###\n",
    "\n",
    "We will design a provably stable DGSEM obeying the energy balance at the discrete level.\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 internal non-frictional interfaces 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^{-}, \\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 transformed hat-variables $\\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. These hat-variables encode the solutions of the IBVP at the element boundaries. For non-linear friction laws we must solve a non-linear algebraic problem to compute hat-variables at the element boundaries. Once the hat-variables are computed we communicate them to the element by penalizing hat variables against the incoming characteristics only. This gives the numerical fluctuations $F^k(-1, t)$ and $G^k(1, t)$ at the element boundaries. \n",
    "\n",
    "6) Finally, the flux fluctuations are appended to the semi-discrete PDE with special penalty weights and we have \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",
    "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(a)$ and the stiffness matrix $Q $ 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 elastic wave equation in a hoterogeneous medium with a nonlinear frictional interface, modeled by a linear friction law (LN), slip-weakening friction law (SW) or the rate-and-state friction law (RS)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**** Exercises****\n",
    "\n",
    "1. In this notebook we use linear, slip-weakening and rate-and-state friction laws. First choose the linear friction law. Initially the friction coefficient, alpha, is set to 1e1000000. What happens if you change that value to 0.1?\n",
    "2. Now change the frinction law to be slip-weakening. Run the simulation keeping the Tau_0 = 81.24+0.1*0.36, initial load, as it is. What happens if we decrease the initial load below 81.24?\n",
    "3. Change the friction law to be rate-and-state and run the simulation. Observe time evolution of slip-rate, stress on fault and slip.\n",
    "4. 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 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",
    "\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,
     21,
     29,
     53,
     58,
     64,
     68,
     110,
     114,
     122,
     130,
     179,
     191
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAGpCAYAAAAA1HueAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3df5hlBX3f8ffHBTQaCBBWXXdXl5itZuOTKLklWNPG+qu7aFib2gTSKA/aTGghlSbRrrVtmv56bJIaJaWQiZJgYqU26sPWrCKimMQ+KLOg6AY3bAjKuqusNAoRE1j59o97NlzHOzNnZnbm3nvm/Xqeee49v+75cJ7hzGfPueecVBWSJEmafI8ZdQBJkiQdHxY7SZKkjrDYSZIkdYTFTpIkqSMsdpIkSR1xwqgDrKYzzjijtmzZMuoYkjrm7rvvxn2LpONt7969X6mq9YtZZk0Vuy1btjAzMzPqGJI6ptfruW+RdNwl+fxil/FUrCRJUkdY7CRJkjrCYidJktQRFjtJkqSOsNhJkiR1hMVOkiSpIyx2kiRJHWGxkyRJ6giLnSRJUkdY7CRJkjrCYidJktQRFjtJkqSOsNhJkiR1hMVOkiSpIyx2kiRJHWGxkyRJ6giLnSRJUkdY7CRJkjrCYidJktQRIy12SbYn2Z/kQJJdQ6YnyeXN9NuTnDVr+roktyV5/+qlliRJGk8jK3ZJ1gFXADuAbcAFSbbNmm0HsLX5mQKunDX9tcAdKxxVkiRpIozyiN3ZwIGququqHgKuBXbOmmcn8I7quxk4NckGgCSbgJcCb1vN0JIkSeNqlMVuI3DPwPDBZlzbed4CvB54ZL6VJJlKMpNk5siRI8tLLEmN6elper0evV4P9y2SxsUoi12GjKs28yR5GXBvVe1daCVVNV1VvarqrV+/fik5JenbTE1NMTMzw8zMDO5bJI2LURa7g8DmgeFNwKGW8zwPOC/J3fRP4b4gye+tXFRJkqTxN8pidwuwNcmZSU4Czgd2z5pnN/Cq5urYc4CvVdXhqnpDVW2qqi3Nch+pqp9e1fSSJElj5oRRrbiqjia5FLgeWAdcXVX7klzcTL8K2AOcCxwAHgQuGlVeSZKkcTeyYgdQVXvol7fBcVcNvC/gkgU+4ybgphWIJ0mSNFF88oQkSVJHWOwkSZI6wmInSZLUERY7SZKkjrDYSZIkdYTFTpIkqSMsdpIkSR1hsZMkSeoIi50kSVJHWOwkSZI6wmInSZLUERY7SZKkjrDYSZIkdYTFTpIkqSMsdpIkSR1hsZMkSeoIi50kSVJHWOwkSZI6wmInSZLUERY7SZKkjrDYSZIkdYTFTpIkqSNGWuySbE+yP8mBJLuGTE+Sy5vptyc5qxm/OclHk9yRZF+S165+ekmSpPEysmKXZB1wBbAD2AZckGTbrNl2AFubnyngymb8UeAXqur7gHOAS4YsK0mStKaM8ojd2cCBqrqrqh4CrgV2zppnJ/CO6rsZODXJhqo6XFW3AlTVA8AdwMbVDC9JkjRuTphrQpKfb7H816vqN5e47o3APQPDB4EfbjHPRuDwQM4twHOATwxbSZIp+kf7eOpTn7rEqJL0raanp5mengbgyJEjI04jSX3zHbF7HfCdwMnz/PzCMtadIeNqMfMk+U7gPcBlVXX/sJVU1XRV9aqqt379+iWHlaRBU1NTzMzMMDMzg/sWSeNiziN2wO9W1X+Yb+EkT1jGug8CmweGNwGH2s6T5ET6pe6dVfXeZeSQJEnqhDmP2FXV6xdauM0887gF2JrkzCQnAecDu2fNsxt4VXN17DnA16rqcJIAbwfuqKo3LyODJElSZyx48USS1yY5pSlXb09ya5KXLHfFVXUUuBS4nv7FD++uqn1JLk5ycTPbHuAu4ADwW8A/b8Y/D3gl8IIkn2p+zl1uJkmSpEk236nYY15dVW9N8g+A9cBFwG8DH1ruyqtqD/3yNjjuqoH3BVwyZLk/Zvj37yRJktasNrc7OVagzgV+u6o+jaVKkiRp7LQpdnuTfIh+sbs+ycnAIysbS5IkSYs1333sTmi+B/ca4NnAXVX1YJLvpn86VpIkSWNkvu/Y3ZzkIPBB4INV9VWAqroPuG81wkmSJKm9OYtdVfWSPI3+81rfkmQj8MfAB4CPVdVfr1JGSZIktTDvd+yq6vNVdVVVvRz4O8D/AV4E/FGSP1iNgJIkSWqnze1OAKiqh4GPND80R/AkSZI0JtrcoPhlSW5L8hdJ7k/yQJL7q+qLqxFQkiRJ7bQ5YvcW4MeBzzQ3DJYkSdIYanMfu3uAz1rqJEmSxlubI3avB/Yk+RjwN1fCVtWbVyyVJEmSFq1NsfvPwF8CjwNOWtk4kiRJWqo2xe70qnrJiieRJEnSsrT5jt2Hk1jsJEmSxlybYncJ8MEk3xi83clKB5MkSdLiLHgqtqpOXo0gkiRJWp45j9glefJCC7eZR5IkSatjvlOxe1os32YeSZIkrYL5TsX+4ALfpQvgd+0kSZLGxJzFrqrWrWYQSZIkLU+bq2IlSZI0ASx2kiRJHWGxkyRJ6ogFi12SX0vy/Sux8iTbk+xPciDJriHTk+TyZvrtSc5qu6wkSdJa0+aI3eeA6SSfSHJxku86HitOsg64AtgBbAMuSLJt1mw7gK3NzxRw5SKWlSRJWlPaPHnibcDbkjwDuAi4PcnHgd+qqo8uY91nAweq6i6AJNcCO4E/GZhnJ/COqirg5iSnJtkAbGmx7Lfbvx+e//xlRJakIdy3SBoTrb5j1xwhe2bz8xXg08DPN4VqqTYC9wwMH2zGtZmnzbLHsk8lmUky8/DDDy8jriQ9avrQIXp799Lbuxf3LZLGxYJH7JK8Gfgx4CPAf6mqTzaT/muS/ctYd4aMq5bztFm2P7JqGpgG6PV6xU03LSKiJA031fwA9Ho93LdIOu4yrO7Mb8FiB3wW+DdV9eCQaWcveo2POghsHhjeBBxqOc9JLZaVJElaU9qciv0ns0tdkhsBqupry1j3LcDWJGcmOQk4H9g9a57dwKuaq2PPAb5WVYdbLitJkrSmzHnELsnjgMcDZyQ5jUdPf54CPGW5K66qo0kuBa4H1gFXV9W+JBc3068C9gDnAgeAB+lfvDHnssvNJEmSNMnmOxX7s8Bl9EvcrQPj76d/q5Flq6o99Mvb4LirBt4XcEnbZSVJktayOYtdVb0VeGuSn6uq31jFTJIkSVqC+U7FvqCqPgJ8McmPz55eVe9d0WSSJElalPlOxf4o/Vuc/NiQaQVY7CRJksbIfKdif6l5vWj14kiSJGmpFrzdSZL/kuTUgeHTkvynlY0lSZKkxWpzH7sdVfXVYwNV9Rf0b0EiSZKkMdKm2K1L8thjA0m+A3jsPPNLkiRpBNo8Uuz3gBuT/Db9iyZeDVyzoqkkSZK0aAsWu6r6lSS3Ay9qRv3Hqrp+ZWNJkiRpsdocsQO4DTiR/hG721YujiRJkpaqzVWxPwF8EngF8BPAJ5K8YqWDSZIkaXHaHLF7I/C3q+pegCTrgQ8Dv7+SwSRJkrQ4ba6KfcyxUte4r+VykiRJWkVtjth9MMn1wLua4Z8E9qxcJEmSJC1Fm6tiX5fkHwHPAwJMV9X7VjyZJEmSFqXVVbFV9R7gPSucRZIkScswZ7FL8gD925t82ySgquqUFUslSZKkRZuz2FXVyasZRJIkScvT6urWJD+S5KLm/RlJzlzZWJIkSVqsNjco/iXgXwFvaEadRP/5sZIkSRojbY7Y/UPgPODrAFV1CPA0rSRJ0phpU+weqqqiuZAiyRNWNpIkSZKWok2xe3eS3wROTfIz9B8n9lvLWWmS05PckOTO5vW0OebbnmR/kgNJdg2M/9Ukn0tye5L3JTl1OXkkSZK6YMFiV1W/Rv+5sO8BngH8u6r6jWWudxdwY1VtBW5shr9FknXAFcAOYBtwQZJtzeQbgGdV1Q8Af8qj3/+TJElasxa8QXGSfwn876q64Tiudyfw/Ob9NcBN9C/QGHQ2cKCq7mpyXNss9ydV9aGB+W4GXnEcs0mSJE2kNqdiTwGuT/JHSS5J8qTjsN4nVdVhgOb1iUPm2QjcMzB8sBk326uBD8y1oiRTSWaSzBw5cmQZkSXpUdPT0/R6PXq9Hu5bJI2LNqdif7mqvh+4BHgK8LEkH15ouSQfTvLZIT87W2bLsDiz1vFG4CjwznnyT1dVr6p669evb7lqSZrf1NQUMzMzzMzM4L5F0rho9azYxr3Al4D7GH6E7VtU1Yvmmpbky0k2VNXhJBuaz57tILB5YHgTcGjgMy4EXga8sLlqV5IkaU1rc4Pif5bkJvoXOZwB/Exz0cJy7AYubN5fCFw3ZJ5bgK1JzkxyEnB+sxxJttP/Tt55VfXgMrNIkiR1Qpsjdk8DLquqTx3H9b6J/m1UXgN8AfjHAEmeArytqs6tqqNJLgWuB9YBV1fVvmb5/w48FrghCcDNVXXxccwnSZI0cRYsdlX1bbciWa6qug944ZDxh4BzB4b3AHuGzPe9xzuTJEnSpGtzVawkSZImgMVOkiSpIyx2kiRJHWGxkyRJ6giLnSRJUkdY7CRJkjrCYidJktQRFjtJkqSOsNhJkiR1hMVOkiSpIyx2kiRJHWGxkyRJ6giLnSRJUkdY7CRJkjrCYidJktQRFjtJkqSOsNhJkiR1hMVOkiSpIyx2kiRJHWGxkyRJ6giLnSRJUkdY7CRJkjrCYidJktQRIyl2SU5PckOSO5vX0+aYb3uS/UkOJNk1ZPovJqkkZ6x8akmSpPE2qiN2u4Abq2orcGMz/C2SrAOuAHYA24ALkmwbmL4ZeDHwhVVJLEmSNOZGVex2Atc0768BXj5knrOBA1V1V1U9BFzbLHfMrwOvB2olg0qSJE2KURW7J1XVYYDm9YlD5tkI3DMwfLAZR5LzgC9W1acXWlGSqSQzSWaOHDmy/OSSBExPT9Pr9ej1erhvkTQuTlipD07yYeDJQya9se1HDBlXSR7ffMZL2nxIVU0D0wC9Xs+je5KOi6mpKaampgDo9XojTiNJfStW7KrqRXNNS/LlJBuq6nCSDcC9Q2Y7CGweGN4EHAKeDpwJfDrJsfG3Jjm7qr503P4DJEmSJsyoTsXuBi5s3l8IXDdknluArUnOTHIScD6wu6o+U1VPrKotVbWFfgE8y1InSZLWulEVuzcBL05yJ/0rW98EkOQpSfYAVNVR4FLgeuAO4N1VtW9EeSVJksbeip2KnU9V3Qe8cMj4Q8C5A8N7gD0LfNaW451PkiRpEvnkCUmSpI6w2EmSJHWExU6SJKkjLHaSJEkdYbGTJEnqCIudJElSR1jsJEmSOsJiJ0mS1BEWO0mSpI6w2EmSJHWExU6SJKkjLHaSJEkdYbGTJEnqCIudJElSR1jsJEmSOsJiJ0mS1BEWO0mSpI5IVY06w6pJ8gCwf9Q5luAM4CujDrFEk5p9UnPD5Gaf1NwAPwTsHXWIJZjkbT6p2Sc1N0xu9knNDfCMqjp5MQucsFJJxtT+quqNOsRiJZmZxNwwudknNTdMbvZJzQ2Q5OuTmH3Ct/lEZp/U3DC52Sc1N/SzL3YZT8VKkiR1hMVOkiSpI9ZasZsedYAlmtTcMLnZJzU3TG72Sc0N8N5RB1iiSd7mk5p9UnPD5Gaf1NywhOxr6uIJSZKkLltrR+wkSZI6y2InSZLUERY7SZKkjrDYSZIkdYTFTpIkqSMsdpIkSR1hsZMkSeoIi50kSVJHWOwkSZI6wmInSZLUERY7SZKkjrDYSZIkdYTFTpIkqSMsdpIkSR1hsZMkSeoIi50kSVJHWOwkSZI6wmInSZLUERY7SZKkjrDYSZIkdYTFTpIkqSMsdpIkSR1hsZMkSeoIi50kSVJHWOwkSZI6wmInSZLUERY7SZKkjrDYSZIkdYTFTpIkqSNOGHWA1XTGGWfUli1bRh1DUsfcfffduG+RdLzt3bv3K1W1fjHLrKlit2XLFmZmZkYdQ1LH9Ho99y2Sjrskn1/sMp6KlSRJ6giLnSRJUkdY7CRJkjrCYidJktQRFjtJkqSOsNhJkiR1hMVOkiSpIyx2kiRJHWGxkyRJ6giLnSRJUkdY7CRJkjrCYidJktQRFjtJkqSOsNhJkiR1hMVOkiSpIyx2kiRJHWGxkyRJ6giLnSRJUkeMtNgl2Z5kf5IDSXYNmZ4klzfTb09y1qzp65LcluT9q5dakiRpPI2s2CVZB1wB7AC2ARck2TZrth3A1uZnCrhy1vTXAnescFRJkqSJMMojdmcDB6rqrqp6CLgW2Dlrnp3AO6rvZuDUJBsAkmwCXgq8bTVDS5IkjatRFruNwD0DwwebcW3neQvweuCR+VaSZCrJTJKZI0eOLC+xJDWmp6fp9Xr0ej3ct0gaF6MsdhkyrtrMk+RlwL1VtXehlVTVdFX1qqq3fv36peSUpG8zNTXFzMwMMzMzuG+RNC5GWewOApsHhjcBh1rO8zzgvCR30z+F+4Ikv7dyUSVJksbfKIvdLcDWJGcmOQk4H9g9a57dwKuaq2PPAb5WVYer6g1VtamqtjTLfaSqfnpV00uSJI2ZE0a14qo6muRS4HpgHXB1Ve1LcnEz/SpgD3AucAB4ELhoVHklSZLG3ciKHUBV7aFf3gbHXTXwvoBLFviMm4CbViCeJEnSRPHJE5IkSR1hsZMkSeoIi50kSVJHWOwkSZI6wmInSZLUERY7SZKkjrDYSZIkdYTFTpIkqSMsdpIkSR1hsZMkSeoIi50kSVJHWOwkSZI6wmInSZLUERY7SZKkjrDYSZIkdYTFTpIkqSMsdpIkSR1hsZMkSeoIi50kSVJHWOwkSZI6wmInSZLUERY7SZKkjrDYSZIkdcRIi12S7Un2JzmQZNeQ6UlyeTP99iRnNeM3J/lokjuS7Evy2tVPL0mSNF5GVuySrAOuAHYA24ALkmybNdsOYGvzMwVc2Yw/CvxCVX0fcA5wyZBlJUmS1pRRHrE7GzhQVXdV1UPAtcDOWfPsBN5RfTcDpybZUFWHq+pWgKp6ALgD2Lia4SVJksbNKIvdRuCegeGDfHs5W3CeJFuA5wCfGLaSJFNJZpLMHDlyZJmRJalvenqaXq9Hr9fDfYukcTHKYpch42ox8yT5TuA9wGVVdf+wlVTVdFX1qqq3fv36JYeVpEFTU1PMzMwwMzOD+xZJ4+KEuSYkOb3F8o9U1VeXuO6DwOaB4U3AobbzJDmRfql7Z1W9d4kZJEmSOmPOYke/QB1i+FGzY9YBT13ium8BtiY5E/gicD7wU7Pm2Q1cmuRa4IeBr1XV4SQB3g7cUVVvXuL6JUmSOmW+YndHVT1nvoWT3LbUFVfV0SSXAtfTL4hXV9W+JBc3068C9gDnAgeAB4GLmsWfB7wS+EySTzXj/nVV7VlqHkmSpEk3X7F7bovl28wzp6aI7Zk17qqB9wVcMmS5P2b+I4mSJElrzpzFrqr+anA4yROBxw1M/8LseSRJkjQ6C14Vm+S8JHcCfw58DLgb+MAK55IkSdIitbndyX+k/3SHP62qM4EXAh9f0VSSJElatDbF7uGqug94TJLHVNVHgWevcC5JkiQt0nwXTxzz1eZGwH8EvDPJvfSf1SpJkqQx0uaI3U7gG8BlwAeBPwN+bCVDSZIkafEWPGJXVV9P8mTgbOD/Adc3p2YlSZI0RtpcFftPgU8CPw68Arg5yatXOpgkSZIWp8137F4HPOfYUbok3w38X+DqlQwmSZKkxWnzHbuDwAMDww8A96xMHEmSJC1VmyN2XwQ+keQ6oOhfTPHJJD8PUFVvXsF8kiRJaqlNsfuz5ueY65rXk49/HEmSJC1Vm6tif3k1gkiSJGl55ix2SXbPt2BVnXf840iSJGmp5jti91z6F0m8C/gEkFVJJEmSpCWZr9g9GXgxcAHwU8AfAO+qqn2rEUySJEmLM+ftTqrqm1X1waq6EDgHOADclOTnVi2dJEmSWpv34okkjwVeSv+o3RbgcuC9Kx9LkiRJizXfxRPXAM8CPgD8clV9dtVSSZIkadHmO2L3SuDrwN8C/kXyN9dOBKiqOmWFs0mSJGkR5ix2VdXmcWOSJEkaE5Y3SZKkjpiz2CW5daGF28yzwPLbk+xPciDJriHTk+TyZvrtSc5qu6wkSdJaM9937L4vye3zTA/wXUtdcZJ1wBX075V3ELglye6q+pOB2XYAW5ufHwauBH645bKSJElrynzF7pktlv/mMtZ9NnCgqu4CSHItsBMYLGc7gXdUVQE3Jzk1yQb6t15ZaNlvt38/PP/5y4gsSUO4b5E0Jua7eOLzK7zujfQfWXbMQfpH5RaaZ2PLZQFIMgVMAfzAYx+7vMSS1Jg+dIjpw4cBePjhh0ecRpL65r1B8Qob9uzZajlPm2X7I6umgWmAXq9X3HTTIiJK0nB/8y9GoNfr4b5F0nGXYXVnfqMsdgeBzQPDm4BDLec5qcWykiRJa8qCtztJ8oQkj2ne/60k5yU58Tis+xZga5Izk5wEnA/snjXPbuBVzdWx5wBfq6rDLZeVJElaU9ocsftD4O8mOQ24EZgBfhL4J8tZcVUdTXIpcD2wDri6qvYlubiZfhWwBzgXOAA8CFw037LLySNJkjTp2hS7VNWDSV4D/EZV/UqS247HyqtqD/3yNjjuqoH3BVzSdllJkqS1rM2TJ5LkufSP0P1BM26U382TJEnSEG2K3WXAG4D3NadKvwf46MrGkiRJ0mIteOStqj4GfAyguYjiK1X1L1Y6mCRJkhanzVWx/zPJKUmeQP/JDvuTvG7lo0mSJGkx2pyK3VZV9wMvp3+xwlOBV65oKkmSJC1am2J3YnPfupcD11XVw8zxlAdJkiSNTpti95vA3cATgD9M8jTg/pUMJUmSpMVrc/HE5cDlA6M+n+Tvr1wkSZIkLUWbiyeelOTtST7QDG8DLlzxZJIkSVqUNqdif4f+o7ue0gz/Kf1720mSJGmMtCl2Z1TVu4FHoP+cVuCbK5pKkiRJi9am2H09yXfTXAmb5BzgayuaSpIkSYvW5pmvPw/sBp6e5OPAeuAVK5pKkiRJizZvsWseIfY44EeBZwAB9jf3spMkSdIYmbfYVdUjSf5bVT0X2LdKmSRJkrQEbb5j96Ek/yhJVjyNJEmSlqztd+yeABxN8lf0T8dWVZ2yoskkSZK0KG2ePHHyagSRJEnS8rR58sSNbcZJkiRptOY8YpfkccDjgTOSnEb/FCzAKTz6FApJkiSNiflOxf4s/UeHPQXYy6PF7n7gihXOJUmSpEWas9hV1VuBtyb5uar6jVXMJEmSpCVoc7uTLyU5GSDJv0ny3iRnrXAuSZIkLVKbYvdvq+qBJD8C/APgGuDK5aw0yelJbkhyZ/N62hzzbU+yP8mBJLsGxv9qks8luT3J+5Kcupw8kiRJXdCm2H2zeX0pcGVVXQectMz17gJurKqtwI3N8LdIso7+d/l2ANuAC5JsaybfADyrqn4A+FPgDcvMI0mSNPHaFLsvJvlN4CeAPUke23K5+eykf+SP5vXlQ+Y5GzhQVXdV1UPAtc1yVNWHqupoM9/NwKZl5pEkSZp4bQraTwDXA9ur6qvA6cDrlrneJ1XVYYDm9YlD5tkI3DMwfLAZN9urgQ/MtaIkU0lmkswcOXJkGZEl6VHT09P0ej16vR7uWySNizZPnngQeO/A8GHg8ELLJfkw8OQhk97YMtuwZ9PWrHW8ETgKvHOuD6mqaWAaoNfr1VzzSdJiTE1NMTU1BUCv1xtxGknqa/Os2CWpqhfNNS3Jl5NsqKrDSTYA9w6Z7SCweWB4E3Bo4DMuBF4GvLCqLGySJGnNW+535ZZqN3Bh8/5C4Loh89wCbE1yZpKTgPOb5UiyHfhXwHnNEUVJkqQ1b1TF7k3Ai5PcCby4GSbJU5LsAWgujriU/vf77gDeXVX7muX/O3AycEOSTyW5arX/AyRJksbNip2KnU9V3Qe8cMj4Q8C5A8N7gD1D5vveFQ0oSZI0gUZ1xE6SJEnHmcVOkiSpIyx2kiRJHWGxkyRJ6giLnSRJUkdY7CRJkjrCYidJktQRFjtJkqSOsNhJkiR1hMVOkiSpIyx2kiRJHWGxkyRJ6giLnSRJUkdY7CRJkjrCYidJktQRFjtJkqSOsNhJkiR1hMVOkiSpIyx2kiRJHWGxkyRJ6giLnSRJUkdY7CRJkjpiJMUuyelJbkhyZ/N62hzzbU+yP8mBJLuGTP/FJJXkjJVPLUmSNN5GdcRuF3BjVW0FbmyGv0WSdcAVwA5gG3BBkm0D0zcDLwa+sCqJJUmSxtyoit1O4Jrm/TXAy4fMczZwoKruqqqHgGub5Y75deD1QK1kUEmSpEkxqmL3pKo6DNC8PnHIPBuBewaGDzbjSHIe8MWq+vRCK0oylWQmycyRI0eWn1ySgOnpaXq9Hr1eD/ctksbFCSv1wUk+DDx5yKQ3tv2IIeMqyeObz3hJmw+pqmlgGqDX63l0T9JxMTU1xdTUFAC9Xm/EaSSpb8WKXVW9aK5pSb6cZENVHU6yAbh3yGwHgc0Dw5uAQ8DTgTOBTyc5Nv7WJGdX1ZeO23+AJEnShBnVqdjdwIXN+wuB64bMcwuwNcmZSU4Czgd2V9VnquqJVbWlqrbQL4BnWeokSdJaN6pi9ybgxUnupH9l65sAkjwlyR6AqjoKXApcD9wBvLuq9o0oryRJ0thbsVOx86mq+4AXDhl/CDh3YHgPsGeBz9pyvPNJkiRNIp88IUmS1BEWO0mSpI6w2EmSJHWExU6SJKkjLHaSJEkdYbGTJEnqCIudJElSR1jsJEmSOsJiJ0mS1BEWO0mSpI6w2EmSJHWExU6SJKkjLHaSJEkdYbGTJEnqCIudJElSR1jsJEmSOsJiJ0mS1BEWO0mSpI5IVY06w6pJ8gCwf9Q5luAM4CujDrFEk5p9UnPD5Gaf1NwAPwTsHXWIJZjkbT6p2Sc1N0xu9knNDfCMqjp5MQucsFJJxtT+quqNOsRiJZmZxNwwudknNTdMbvZJzQ2Q5OuTmH3Ct/lEZp/U3DC52Sc1N/SzL3YZT8VKkiR1hMVOkiSpI9ZasZsedYAlmtTcMLnZJzU3TG72Sc0N8N5RB1iiSd7mk5p9UnPD5Gaf1NywhOxr6uIJSZKkLltrR+wkSZI6y+U3S1sAAAbDSURBVGInSZLUEWui2CXZnmR/kgNJdo06z2IkuTvJZ5J8aimXPa+mJFcnuTfJZwfGnZ7khiR3Nq+njTLjMHPk/vdJvths908lOXeUGYdJsjnJR5PckWRfktc24ydhm8+Vfay3e5LvSvKXSb6R5K+S3NSM/54k9yV5qHndMtKgQyR5XJJPJvl0s81/uRk/1tv8mCTrktyW5P3N8Nj/nh8zJPvYb/Nhf3smZZvPkX0StvmpSX4/yeeafeNzl7LNO1/skqwDrgB2ANuAC5JsG22qRfv7VfXsCbgPz+8A22eN2wXcWFVbgRub4XHzO3x7boBfb7b7s6tqzypnauMo8AtV9X3AOcAlze/2JGzzubLDeG/3+4GnV9V3AKcBvSSvAd4FfLKqTgI+CVw7woxz+WvgBVX1g8Czge1JzmmmjfM2P+a1wB0Dw5Pwe37M7OwwGdt89t+eSdrmw/5ujvs2fyvwwap6JvCD9H9nFr3NO1/sgLOBA1V1V1U9RH+Hu3PEmTqpqv4Q+H+zRu8ErmneXwO8fFVDtTBH7rFXVYer6tbm/QP0dwIbmYxtPlf2sVZ9X24GH09/H1rAc3h0h7sLOGsE8ebVZP/LZvDE5mcirp5Lsgl4KfC2gdFj/3sOc2afVBOxzSdRklOAvwe8HaCqHqqqr7KEbb4Wit1G4J6B4YNMwB+QAQV8KMneJFOjDrMET6qqw9D/Yw48ccR5FuPSJLc3p2rH8pTDMc2pv+cAn2DCtvms7DDm2z3JiUm+Qf8RRbdX1dXAiVX1aYDmdSyf6tOcEvwUcC9wQ1VNxDYH3gK8HnhkYNyk/J4Pyw7jv82H/e2ZlG0+19/Ncd7m3wMcAX67OW3/tiRPYAnbfC0UuwwZNxH/Sm08r6rOon8q+ZIkf2/UgdaIK4Gn0z9ldRj4b6ONM7ck3wm8B7isqu4fdZ7FGJJ97Ld7VT3cnIrdAjwjycQctaiqb1bVs4FNwNlJnsWYb/MkLwPuraqJexbvPNnHeps3Jvlvz7Ds477NT6B/pP/KqnoO8HWWeKp7LRS7g8DmgeFNwKERZVm0qjrUvN4LvI/+qeVJ8uUkGwCa13tHnKeVqvpy80fwEeC3GNPtnuRE+sXonVV17Ca5E7HNh2WflO0OUFWfB24DXgM8nOQHAZrXo6PMtpDmFM9NwPYJ2ObPA85Lcjf9r9K8IMnvMRm/50OzT8A2n+tvzyRs86HZJ2CbHwQODhxF/336RW/R23wtFLtbgK1JzkxyEnA+sHvEmVpJ8oQkJx97D7wE+Oz8S42d3cCFzfsLgetGmKW1Y/8jNf4hY7jdk4T+9zHuqKo3D0wa+20+V/Zx3+5Jnpnkac3704AfAm4FPgW8qZntTfQL31hJsj7Jqc377wBeBHxu3Ld5Vb2hqjZV1Rb6+++PVNVPMwG/53NlH/dtPs/fnrHf5nNlH/dtXlVfAu5J8oxm1AuBP2EJ23wsvwdyPFXV0SSXAtcD64Crq2rfiGO19STgff2/gZwA/M+q+uBoI80tybuA5wNnJDkI/BL9P3Lvbq4c/ALwj0eXcLg5cj8/ybPpn7a/G/jZkQWc2/OAVwKfab43BfCvmYBtztzZLxjz7f79wO82xTTAx6vql5L8LvDJJA8BfwmM4xXsG4BrmjsFPAZ4d1W9P8nvjvk2n8sk/J7P5VfGfJsP/duT5BbGf5vPlX0Sfs9/DnhncxDqLuAimv9XF7PNfaSYJElSR6yFU7GSJElrgsVOkiSpIyx2kiRJHWGxkyRJ6giLnSRJUkdY7CRJkjrCYqfOS7I5yZ8nOb0ZPq0Zftqs+bYk+cbAfdXafv5PJjmQ5P3HM7ek8ea+RePIYqfOq6p76D8ncPDJANPNI6Fm+7PmWZqL+fz/BfzT5aWUNGnct2gcWey0Vvw6cE6Sy4AfocUDoJt/ZX8uyduSfDbJO5O8KMnHk9yZZNyeNShp9blv0Vix2GlNqKqHgdfR3wlfVlUPtVz0e4G3Aj8APBP4Kfo771+k/wgsSWuY+xaNG4ud1pIdwGHgWYtY5s+r6jNV9QiwD7ix+s/h+wyw5fhHlDSB3LdobFjstCY0D39+MXAO8C+TbGi56F8PvH9kYPgR+g+YlrSGuW/RuLHYqfOShP4XnC+rqi8Avwr82mhTSZp07ls0jix2Wgt+BvhCVd3QDP8P4JlJfnSEmSRNPvctGjvpn9KXlGQL8P6qWsz3ZI4t+3zgF6vqZcc5lqQJ575Fq8kjdtKjvgl811JuIkr/X+p/sSKpJE069y1aNR6xkyRJ6giP2EmSJHWExU6SJKkjLHaSJEkdYbGTJEnqiP8P2jIbGrJkzcwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x504 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"
    }
   ],
   "source": [
    "# Initialization and computation\n",
    "# Tic\n",
    "iplot  = 10  \n",
    "\n",
    "# At the moment, keep the length of the two domain the same\n",
    "# Ahysical domain x = [ax, bx] (km)\n",
    "# For the first domain \n",
    "ax_1 = 0.0\n",
    "bx_1 = 30.0\n",
    "\n",
    "# Physical domain x = [ax, bx] (km)\n",
    "# For the second domain\n",
    "ax_2 = bx_1\n",
    "bx_2 = 60.0\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 = 5                        # 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_freedom1 = 301        # for first domain\n",
    "deg_of_freedom2 = 301        # for second domain\n",
    "\n",
    "# estimate the number of elements needed for a given polynomial degree and degrees of freedom\n",
    "num_element1 = round(deg_of_freedom1/NP)       # for first domain\n",
    "num_element2 = round(deg_of_freedom2/NP)       # for second domain\n",
    "\n",
    "# Initialize the mesh\n",
    "y_1 = np.zeros(NP*num_element1)\n",
    "y_2 = np.zeros(NP*num_element2)\n",
    "\n",
    "# Generate num_element dG elements in the interval [ax, bx]\n",
    "x0_1 = np.linspace(ax_1,bx_1,num_element1+1)\n",
    "dx_1 = np.diff(x0_1)                           # element sizes\n",
    "    \n",
    "x0_2 = np.linspace(ax_2,bx_2,num_element2+1)\n",
    "dx_2 = np.diff(x0_2)                           # 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_element1+1):\n",
    "    for j in range (1,(N+2)):\n",
    "        y_1[j+(N+1)*(i-1)-1] = dx_1[i-1]/2.0 * (psi[j-1] + 1.0) +x0_1[i-1]\n",
    "\n",
    "for i in range (1,num_element2+1):\n",
    "    for j in range (1,(N+2)):\n",
    "        y_2[j+(N+1)*(i-1)-1] = dx_2[i-1]/2.0 * (psi[j-1] + 1.0) +x0_2[i-1]\n",
    "\n",
    "deg_of_freedom1 = len(y_1)                     #same for both the domains\n",
    "deg_of_freedom2 = len(y_2)                     #same for both the domains\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.0                          # r=0:absorbing, r=1:free-surface, r=-1: clamped, \n",
    "rn = 0.0                          # r=0:absorbing, r=1:free-surface, r=-1: clamped, \n",
    "\n",
    "# Initialize the wave-fields\n",
    "L_1 = 0.5*(bx_1-ax_1)\n",
    "L_2 = 0.5*(bx_2-ax_2)\n",
    "\n",
    "delta_1 = 0.01*(bx_1-ax_1)\n",
    "delta_2 = 0.01*(bx_2-ax_2)\n",
    "\n",
    "x0_1 = 0.5*(bx_1+ax_1)\n",
    "x0_2 = 0.5*(bx_2+ax_2)\n",
    "\n",
    "# Material parameters\n",
    "c1s = 3.464 + 0.0*np.sin(2.0*np.pi*y_1)             # shear wave speed (km/s)\n",
    "rho1 = 2.67 + 0.0*y_1                               # density (g/cm^3)\n",
    "mu1 = rho1 * c1s**2                                 # shear modulus (GPa)\n",
    "Z1s = rho1*c1s                                      # shear impedance\n",
    "\n",
    "c2s = 3.464 + 0.0*y_2                               # shear wave speed (km/s)\n",
    "rho2 = 2.67 + 0.0*y_2                               # density (g/cm^3)\n",
    "mu2 = rho2 * c2s**2                                 # shear modulus (GPa)\n",
    "Z2s = rho2*c2s                                      # shear impedance\n",
    "\n",
    "# Choose friction law: fric_law\n",
    "# We use linear (LN: T = alpha*v, alpha >=0)\n",
    "# Slip-weakening (SW)\n",
    "# Rate-and-state friction law (RS)\n",
    "\n",
    "fric_law = 'RS'\n",
    "\n",
    "if fric_law not in ('LN', 'SW', 'RS'):\n",
    "    print('friction law not implemented. choose fric_law = SW or fric_law = LN')\n",
    "    exit(-1)\n",
    "\n",
    "if fric_law  in ('LN'):\n",
    "    alpha = 1e1000000             # initial friction coefficient\n",
    "    S = 0.0                       # initial slip (in m)\n",
    "    Tau_0 = 81.24+1*0.36          # initial load (81.24 in MPa), slight increase will unlock the fault\n",
    "    W = 0.0\n",
    "    u_1 = np.exp(-(y_1-x0_1)**2/delta_1)                                                  # Gaussian\n",
    "    u_2 = 1/np.sqrt(2.0*np.pi*delta_2**2)*np.exp(-(y_2-x0_2)**2/(2.0*delta_2**2))         # Gaussian\n",
    "\n",
    "if fric_law  in ('SW'):\n",
    "    alpha = 1e1000000              # initial friction coefficient\n",
    "    S = 0.0                        # initial slip (in m)\n",
    "    Tau_0 = 81.24+0.1*0.36         # initial load (81.24 in MPa), slight increase will unlock the fault\n",
    "    W = 0.0\n",
    "    u_1 = 0.0*np.exp(-(y_1-x0_1)**2/delta_1)                                                # Gaussian\n",
    "    u_2 = 0.0/np.sqrt(2.0*np.pi*delta_2**2)*np.exp(-(y_2-x0_2)**2/(2.0*delta_2**2))         # Gaussian\n",
    "\n",
    "if fric_law  in ('RS'):\n",
    "    alpha = 1e1000000                      # initial friction coefficient                                                                                   \n",
    "    S = 0.0                                # initial slip (in m)                                                                                                                             \n",
    "    Tau_0 = 81.24+0.1*0.36                 # initial load (81.24 in MPa), slight increase will unlock the fault   \n",
    "    W = 0.4367\n",
    "    u_1 = 0.0*np.exp(-(y_1-x0_1)**2/delta_1)                                                # Gaussian\n",
    "    u_2 = 0.0/np.sqrt(2.0*np.pi*delta_2**2)*np.exp(-(y_2-x0_2)**2/(2.0*delta_2**2))         # Gaussian\n",
    "\n",
    "u_1 = np.transpose(u_1)\n",
    "v_1 = np.transpose(u_1)\n",
    "    \n",
    "u_2 = np.transpose(u_2)\n",
    "v_2 = np.transpose(u_2)\n",
    "\n",
    "# Time stepping parameters\n",
    "cfl = 0.5                            # CFL number\n",
    "dt = (cfl/(max(max(c1s),max(c2s))*(2*N+1)))*min(min(dx_1), min(dx_2))         # time-step (s)\n",
    "t = 0.0                              # initial time\n",
    "Tend = 3.0                           # final time (s)\n",
    "n = 0                                # counter\n",
    "\n",
    "Vd = [0]\n",
    "Sd = [Tau_0]\n",
    "T = [t]\n",
    "Slip = [0]\n",
    "\n",
    "# Initialize animated plot   \n",
    "f1, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, sharex='col', sharey='row',figsize=(10,7))\n",
    "line1 = ax1.plot(y_1, u_1,'r')\n",
    "line2 = ax2.plot(y_2, u_2,'r')\n",
    "line3 = ax3.plot(y_1, v_1,'r')\n",
    "line4 = ax4.plot(y_2, v_2,'r')\n",
    "    \n",
    "ax1.set_ylabel('velocity [m/s]')\n",
    "ax3.set_ylabel('stress [Mpa]')\n",
    "ax3.set_xlabel('X [m]')\n",
    "ax4.set_xlabel('X [m]')\n",
    "\n",
    "ax3.set_xlim([ax_1, bx_1])\n",
    "ax4.set_xlim([ax_2,bx_2])\n",
    "f1.subplots_adjust(wspace=0)\n",
    "plt.ion()\n",
    "plt.show()\n",
    "\n",
    "time_integrator = 'ADER'\n",
    "\n",
    "for t in utils.drange (dt,Tend+dt,dt):\n",
    "    n = n+1\n",
    "    \n",
    "    if time_integrator  in ('ADER'):\n",
    "        ADER_Wave_1D_GL_return = timeintegrate.ADER_Wave_1D_GL(u_1,v_1,u_2,v_2,S,W,D,NP,num_element1,\\\n",
    "                                                            num_element2,dx_1,dx_2,w,psi,t,r0,rn,dt,\\\n",
    "                                                             rho1,mu1,rho2,mu2,alpha,Tau_0,fric_law)\n",
    "        u_1 = ADER_Wave_1D_GL_return['Hu_1']\n",
    "        v_1 = ADER_Wave_1D_GL_return['Hv_1']\n",
    "        u_2 = ADER_Wave_1D_GL_return['Hu_2']\n",
    "        v_2 = ADER_Wave_1D_GL_return['Hv_2']\n",
    "        S   = ADER_Wave_1D_GL_return['H_d']\n",
    "        W   = ADER_Wave_1D_GL_return['H_psi']\n",
    "        \n",
    "        \n",
    "    if time_integrator  in ('RK'):\n",
    "        RK4_Wave_1D_GL_return = timeintegrate.ADER_Wave_1D_GL(u_1,v_1,u_2,v_2,S,W,D,NP,num_element1,\\\n",
    "                                                               num_element2,dx_1,dx_2,w,psi,t,r0,rn,dt,\\\n",
    "                                                               rho1,mu1,rho2,mu2,alpha,Tau_0,fric_law)\n",
    "        u_1 = RK4_Wave_1D_GL_return['Hu_1']\n",
    "        v_1 = RK4_Wave_1D_GL_return['Hv_1']\n",
    "        u_2 = RK4_Wave_1D_GL_return['Hu_2']\n",
    "        v_2 = RK4_Wave_1D_GL_return['Hv_2']\n",
    "        S   = RK4_Wave_1D_GL_return['H_d']\n",
    "        W   = RK4_Wave_1D_GL_return['H_psi']\n",
    "        \n",
    "        \n",
    "    \n",
    "    Vd.append(np.abs(u_1[-1]-u_2[0]))\n",
    "    Sd.append(v_2[0]+1.0*Tau_0)\n",
    "    T.append(t)\n",
    "    Slip.append(S)\n",
    "      \n",
    "    # Update plot                                                                                                                                                                                  \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_1, u_1, 'r')\n",
    "        line2 = ax2.plot(y_2, u_2, 'r')\n",
    "        line3 = ax3.plot(y_1, Tau_0+v_1, 'r')\n",
    "        line4 = ax4.plot(y_2, Tau_0+v_2, 'r')\n",
    "        plt.gcf().canvas.draw()\n",
    "\n",
    "plt.ioff()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "code_folding": [
     0
    ],
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAHwCAYAAAC7apkrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdebwkVX3//9d7hmHfBRGQAVkUlbiO4r6AUcEFk7gbXL4qRuNCYuKSGEV/JhqXGHclcUGDGFRUQFARRQWFCAhBRMIi+wgCDjODwDDM5/dH1WWay116Zqr79r3zej4e9ehazqlzuqtr5nNPf6oqVYUkSZKkxryZ7oAkSZI0SgyQJUmSpB4GyJIkSVIPA2RJkiSphwGyJEmS1MMAWZIkSephgCxpVkuyMMnyJPNnui9rK0kl2XMt6z4+yYVd96mPdu+X5JdJliV54xrUe3mSU3uWlyfZfTC9lKS1Y4AsaVZJclmSp4wtV9UVVbV5Vd0xk/0alvHBdFX9tKruNwNdeQtwSlVtUVUfW9udtMfu0g77NakkX0zy3mG0JWl2M0CWJK2NXYHzZ7oTY5JsMNN9kDR3GCBLmjWSfBlYCBzX/jT/liS7taOqG7RlTkny3iQ/a8scl+QeSY5MsjTJL5Ls1rPPvZOclOTGJBcmef4U7W+V5HNJFie5um1nfpKNkixJsk9P2e2T3JLknu3yq5Nc3LZzbJKdJmnjlCSv6lm+MyUhyU/a1ee27+0FSZ6U5Kqe8vdv97EkyflJnt2z7YtJPpnkO21qxBlJ9pji/T673ceSdp/3b9f/EHgy8Im2H/edoO7Lk1zatvPbJC+ZpI07R8Tb/n2mPR7Lkvw4ya6T1Bs77q9McgXww3b915L8LslNSX6S5IHt+kOAlwBvGftetOt3SvKNJL9v+9l3uoikucsAWdKsUVUHA1cAz2p/mv/AJEVfCBwM7AzsAfwc+AKwLXAB8C6AJJsBJwFfAe4JvAj41FhQNYEjgJXAnsBDgacCr6qq24Bj2vpjng/8uKquS7If8L523Y7A5cBX1+L9P6GdfXD7/v+7d3uSBcBxwPfb9/MG4MgkvSkYLwLeDWwDXAz880RttUHvUcChwPbACTR/mGxYVfsBPwVe3/bj/8bV3Qz4GHBAVW0BPAY4p8+3+RLg/wO2a+scOU35JwL3B57WLp8I7EXz/s8eq19Vh7fzH2j7/Kwk82g+r3Npviv7A4cmeRqS1msGyJLmoi9U1SVVdRNNwHRJVf2gqlYCX6MJbgGeCVxWVV+oqpVVdTbwDeC543eYZAfgAODQqrq5qq4DPkITjEMTZPcGyC9u10ET9H2+qs5ug+m3A4/uHcnuyKOAzYH3V9WKqvohcPy4fh1TVf/TfhZHAg+ZZF8vAL5TVSdV1e3Ah4BNaILdfqwC9kmySVUtrqp+0zG+U1U/aT+nf6T5nHaZovxh7fG4BaCqPl9Vy9r6hwEPTrLVJHUfAWxfVe9pP69Lgf9g9TGVtJ4yZ0vSXHRtz/wtEyxv3s7vCuybZEnP9g2AL0+wz12BBcDiJGPr5gFXtvM/BDZJsi/wO5rA85vttp1oRjMBqKrlSW6gGbW8bE3e2DR2Aq6sqlU96y5v2xnzu575P7L6s5hoX5ePLVTVqiRXjtvXhKrq5iQvAP4O+FyS04A3V9Vv+ngPY5/n2Od0Y9uXK6crn+ZOJv8MPI9m1Hvsc9gOuGmCursCO407/vNpRsclrccMkCXNNtXhvq6kSYP40z7L3gZs146+3rVTTQB5NM1o7bXA8VW1rN18DU0wBtyZgnAP4OoJ2rkZ2LRn+V79vJGednZJMq8nSF4I/N8Udaba15+MLaT5q2AXJu7z3VTV94DvJdkEeC/NyOzj+6h652hxks1p0mKumaqpnvkXAwcBT6H5w2Mr4A9AJigLzTH9bVXt1Ue/JK1HTLGQNNtcC3R139zjgfsmOTjJgnZ6xNjFaL2qajFNbu+Hk2yZZF6SPZI8safYV2hSE17C6vSKsfWvSPKQJBsB/wKcUVWXTdCnc4A/T7Jpe/HaK8dtn+r9n0ETYL+lfS9PAp7FWuQ7A0cDz0iyf5vb/GaaPxB+Nl3FJDu0F/ht1tZZDvR7G74DkzwuyYY0uchnVNVko8fjbdG2dwPNHxn/Mm77+M/uf4ClSd6aZJP2gst9kjyiz/YkzVEGyJJmm/cB72jvrPB367KjdoT3qTQ5p9fQpB/8K7DRJFVeCmwI/JpmZPLrNBfdje1vLEDdiSb3eWz9ycA/0eQ3L6a5cHCyPNePACtogrkjuPtFaocBR7Tv/y533KiqFcCzaXKlrwc+Bby0z9SGu6iqC4G/BD7e7utZNBdHruij+jyagPoa4EaaC+le12fTX6G5iPJG4OE0f2z060s0aSFX0xyj08dt/xzwgPaz+1Z77+xn0aTD/Jbmff4nzcizpPVYqrr8tVKSpLWT5IvAVVX1jpnui6T1myPIkiRJUo+RDpCT/E17k/pfJTkqycYz3SdJkiTNbSObYpFkZ+BU4AFVdUt7dfgJVfXFme2ZJEmS5rKRHkGmuQ3dJmkeIbspU9/qR5IkSVpnI3sf5Kq6OsmHaB4rewvw/ar6/vhySQ4BDgHYbLPNHr733nsPt6OSJEmalc4666zrq2r78etHOcViG5pbIr0AWELzeNivV9V/TVZn0aJFdeaZZw6ph5IkSZrNkpxVVYvGrx/lFIun0Dzh6PdVdTtwDPCYGe6TJEmS5rhRDpCvAB7VPk0qwP7ABTPcJ0mSJM1xIxsgt0+k+jpwNnAeTV8Pn9FOSZIkac4b2Yv0AKrqXTSPHJUkSZKGYmRHkCVJkqSZYIAsSZIk9TBAliRJknoYIEuSJEk9DJAlSZKkHgbIkiRJUg8DZEmSJKmHAbIkSZLUwwBZkiRJ6mGALEmSJPUwQJYkSZJ6GCBLkiRJPQyQJUmSpB4GyJIkSVIPA2RJkiSpxwaD2nGSj/VRbGlVvWNQfZAkSZLW1MACZOAg4J3TlHkbYIAsSZKkkTHIAPkjVXXEVAWSbDPA9iVJkqQ1NrAc5Kr693Utk2TrJF9P8pskFyR5dHc9lCRJku5u4BfpJflAki2TLEhycpLrk/xln9U/Cny3qvYGHgxcMLieSpIkScO5i8VTq2op8EzgKuC+wN9PVynJlsATgM8BVNWKqloyyI5KkiRJwwiQF7SvBwJHVdWNfdbbHfg98IUkv0zyn0k2G18oySFJzkxy5u9///uOuixJkqT11TAC5OOS/AZYBJycZHvg1j7qbQA8DPh0VT0UuJnmrhd3UVWHV9Wiqlq0/fbbd9lvSZIkrYcGFiAn2RGgqt4GPBpYVFW3A3+kuQXcdK4CrqqqM9rlr9MEzJIkSdLADHIE+fNJTk/yfpoL7AJQVTdX1e+mq9yWuTLJ/dpV+wO/HlhvJUmSJAZ4H+SqOiDJxsCTgD8DPpTkCuC7NHemuKKP3bwBODLJhsClwCsG1V9JkiQJBvugEKrqVtqAGCDJfYADgE8kuVdVPXKa+ufQ5C5LkiRJQzHQALlXe9u2m4CvttPyYbUtSZIk9WvgAXKS1wDvAW4Bql1dVbX7oNuWJEmS1tQwRpD/DnhgVV0/hLYkSZKkdTKM+yBfQnNrN0mSJGnkDWME+e3Az5KcAdw2trKq3jiEtiVJkqQ1MowA+bPAD4HzgFVDaE+SJElaa8MIkFdW1d8OoR1JkiRpnQ0jB/lHSQ5JsmOSbcemIbQrSZIkrbFhjCC/uH19e8+6ArzNmyRJkkbOwAPkqrrPoNuQJEmSujKwFIskD+uijCRJkjRMgxxB/kKSJwGZoszngIcOsA+SJEnSGhlkgLwVcBZTB8i/H2D7kiRJ0hobWIBcVbsNat+SJEnSoAzjNm+SJEnSrGGALEmSJPUwQJYkSZJ6DDxATuMvk7yzXV6Y5JFrUH9+kl8mOX5wvZQkSZIawxhB/hTwaOBF7fIy4JNrUP9NwAVdd0qSJEmayDAC5H2r6q+BWwGq6g/Ahv1UTHJv4BnAfw6ue5IkSdJqwwiQb08yHyiAJNsDq/qs++/AW9agvCRJkrROhhEgfwz4JnDPJP8MnAq8b7pKSZ4JXFdVZ01T7pAkZyY58/e/97kjkiRJWjepqsE3kuwN7E/zVL2Tq2ranOIk7wMOBlYCGwNbAsdU1V9OVmfRokV15plndtNpSZIkzWlJzqqqRePXD+MuFl+uqt9U1Ser6hNVdUGSL09Xr6reXlX3bp/I90Lgh1MFx5IkSVIXhpFi8cDehTYf+eFDaFeSJElaYxsMasdJ3g78A7BJkqU06RUAK4DD12RfVXUKcEqX/ZM0AFXNtGrV2k9rW3+s7d6pt0+TTf2UGbV9jf/M+11el7qzZV+j1M9RNRv6ORv6CLOjn7OhjwCPexz8xV/MdC/uNLAAuareB7wvyfuq6u2DakcaulWr4Pbbm2nFitXz6zqtXAl33LH6dV3mu9jH2gS1s+UfYs09yeTLU21b0+U1rTuqZkM/Z0MfYXb0czb0cYMN1o8AeUxVvT3JNsBeNBfbja3/yaDbnjGrVsFpp8EjHgEbbzx9ed3VqlVw223NtGLF6vnxU9fbVqzoL+BdNeS7Ds6b1/zDMX9+M63J/Ph1G23U/z7mz2/anmpKpi/T77Qu+0omn2Dq7f2WGbV99RpkgDeq+5KkARp4gJzkVTRPw7s3cA7wKODnwH6DbntGXHQR/L//B6eeCrvvDp/+NDz1qTPdq7WzciXceivccsvkr1NtW5vXW29tgtCuzJvXBIW904Yb3n3dNts06zfcEBYsWLtpXeouWNAEp+OD1bHgT5IkDc3AA2Sa4PgRwOlV9eT2lm/vHkK7M+NlL4MLLoD3vhf+67/ggAPgM5+BV7+6/31UwbXXNsH20qXN6OattzYB02abweabN6+906abNoHlLbfAH//YvN58M9x0UzMtWbJ6fqJp2bK7B6srV67955DAJps008Yb3/11q61ghx3uvn5smiiInS7InWjbBsP4ikuSpLlkGNHDrVV1axKSbFRVv0lyvyG0O3zXXgunnw7vfjf84z/C3/wNPPe5cMghcMUVcNhhzcjgRFauhJ/9DL797Wa65JLB9XOLLZoAdWy65z1hzz0nD2b7eR2/bsECRz4lSdKsNIwA+aokWwPfAk5K8gfgmiG0O3wnntiM/j7zmc3yppvCt74Fr31tM6J82mnwnvfAYx/bBI9Ll8LJJzcB8fHHww03NKOf++8Pr3897L033OMeq0dDV65sRoWXL29ee6c//rEJSjfddHXAuummTQC89darg+Ettpg8SJckSdJQLtL7s3b2sCQ/ArYCvjvodmfE8cfDTjvBQx6yet2GG8J//ic86lHw1rfC4x+/Ok3i2mubMltvDc94Bhx0EDz96U0QK0mSpBkx0AA5yTzgf6tqH4Cq+vEg25tRK1bA978PL3rRxFdfv/rV8OIXw1e/Cuef3+T97rkn7LtvEzQvWDAz/ZYkSdJdDDRArqpVSc5NsrCqrhhkWzPu9NObC90OPHDyMpttBq985fD6JEmSpDU2jBzkHYHzk/wPcPPYyqp69hDaHp6LLmpeH/zgme2HJEmS1skwAuS5e0u3Xldc0dyzduedZ7onkiRJWgfDuEjvzrzjJM+squMH3eaMuPzy5gI9c4klSZJmtXlDbu89Q25veK64AhYunOleSJIkaR0NO0Ceu0+OuPxy2HXXme6FJEmS1tGwA+TXDLm94Vi1Cq680hFkSZKkOWDgOchJNgZeBzwOqCSnAp+uqlsH3fbQXHst3H67AbIkSdIcMIy7WHwJWAZ8vF1+EfBl4HlDaHs4Lr+8eTXFQpIkadYbRoB8v6rqvTnwj5KcO4R2h+eK9hkojiBLkiTNesPIQf5lkkeNLSTZFzhtukpJdknyoyQXJDk/yZsG2st1YYAsSZI0ZwxjBHlf4KVJxh41vRC4IMl5QFXVgyaptxJ4c1WdnWQL4KwkJ1XVr4fQ5zVz+eWw1VbNJEmSpFltGAHy09emUlUtBha388uSXADsDIxegOw9kCVJkuaMgQXISbasqqU0F+jdTVXduAb72g14KHDGBNsOAQ4BWDhTQeoVV8Auu8xM25IkSerUIEeQvwI8EzgLKO76kJACdu9nJ0k2B74BHNoG3HdRVYcDhwMsWrSo1rHPa+f66+EhD5mRpiVJktStgQXIVfXM9vU+a7uPJAtoguMjq+qYrvrWuZtugq23nuleSJIkqQODTLF42FTbq+rsaeoH+BxwQVX9W5d969Qdd8CyZV6gJ0mSNEcMMsXiw1NsK2C/aeo/FjgYOC/JOe26f6iqE7roXGeWtlkfBsiSJElzwiBTLJ68jvVP5a55y6PpppuaVwNkSZKkOWHgDwpJ8rz2PsYkeUeSY5I8dNDtDs1YgGwOsiRJ0pwwjCfp/VN7H+PHAU8DjgA+M4R2h2PJkubVEWRJkqQ5YRgB8h3t6zOAT1fVt4ENh9DucJhiIUmSNKcMI0C+OslngecDJyTZaEjtDocBsiRJ0pwyjED1+cD3gKdX1RJgW+Dvh9DucJiDLEmSNKcM8jZvAFTVH4FjepYXA4sH3e7QmIMsSZI0p8ydVIeZctNNsPHGsOHcSauWJElanxkgr6ubbnL0WJIkaQ4xQF5XN91k/rEkSdIcYoC8rpYscQRZkiRpDjFAXlemWEiSJM0pBsjrygBZkiRpTjFAXlfmIEuSJM0pBsjryhxkSZKkOcUAeV3cfjvccosBsiRJ0hxigLwuxh4zbYAsSZI0Zxggr4uxANkcZEmSpDnDAHldLFnSvDqCLEmSNGeMdICc5OlJLkxycZK3zXR/7sYUC0mSpDlnZAPkJPOBTwIHAA8AXpTkATPbq3EMkCVJkuackQ2QgUcCF1fVpVW1AvgqcNAM9+muDJAlSZLmnA1mugNT2Bm4smf5KmDf8YWSHAIc0i4uT3LhEPrWazv22OP6IbepiW0HeCxGg8didHgsRofHYnR4LEbHTB+LXSdaOcoBciZYV3dbUXU4cPjguzOxJGdW1aKZal+reSxGh8didHgsRofHYnR4LEbHqB6LUU6xuArYpWf53sA1M9QXSZIkrSdGOUD+BbBXkvsk2RB4IXDsDPdJkiRJc9zIplhU1cokrwe+B8wHPl9V589wtyYyY+kduhuPxejwWIwOj8Xo8FiMDo/F6BjJY5Gqu6X1SpIkSeutUU6xkCRJkobOAFmSJEnqYYDcp+kee53Gx9rt/5vkYTPRz/VBH8fiSUluSnJOO71zJvo51yX5fJLrkvxqku2eE0PSx7HwnBiSJLsk+VGSC5Kcn+RNE5Tx3BiCPo+F58YQJNk4yf8kObc9Fu+eoMxInRcje5HeKOl57PWf0tx+7hdJjq2qX/cUOwDYq532BT7NBA820brp81gA/LSqnjn0Dq5fvgh8AvjSJNs9J4bni0x9LMBzYlhWAm+uqrOTbAGcleQk/7+YEf0cC/DcGIbbgP2qanmSBcCpSU6sqtN7yozUeeEIcn/6eez1QcCXqnE6sHWSHYfd0fXA6D+CfD1RVT8BbpyiiOfEkPRxLDQkVbW4qs5u55cBF9A8GbaX58YQ9HksNATtd315u7igncbfJWKkzgsD5P5M9Njr8SdZP2W07vr9nB/d/pRzYpIHDqdrGsdzYrR4TgxZkt2AhwJnjNvkuTFkUxwL8NwYiiTzk5wDXAecVFUjfV6YYtGffh573dejsbXO+vmczwZ2bX/KORD4Fs1PNhouz4nR4TkxZEk2B74BHFpVS8dvnqCK58aATHMsPDeGpKruAB6SZGvgm0n2qare6yZG6rxwBLk//Tz22kdjD8e0n3NVLR37KaeqTgAWJNlueF1Uy3NiRHhODFebY/kN4MiqOmaCIp4bQzLdsfDcGL6qWgKcAjx93KaROi8MkPvTz2OvjwVe2l6F+SjgpqpaPOyOrgemPRZJ7pUk7fwjab7nNwy9p/KcGBGeE8PTfs6fAy6oqn+bpJjnxhD0cyw8N4YjyfbtyDFJNgGeAvxmXLGROi9MsejDZI+9TvJX7fbPACcABwIXA38EXjFT/Z3L+jwWzwVem2QlcAvwwvKRkZ1LchTwJGC7JFcB76K58MJzYsj6OBaeE8PzWOBg4Lw23xLgH4CF4LkxZP0cC8+N4dgROKK9E9U84OiqOn6U4ygfNS1JkiT1MMVCkiRJ6mGALEmSJPUwQJYkSZJ6GCBLkiRJPQyQJUmSpB4GyJIkSVIPA2RJGkFJtk7yup7lnZJ8fQDtHJbk6iTvmaLMHknOSbK86/YlaRR5H2RJGkFJdgOOr6p9BtzOYcDyqvpQH2WXV9Xmg+yPJI0CR5AlaTS9Hxgbuf1gkt2S/AogycuTfCvJcUl+m+T1Sf42yS+TnJ5k27bcHkm+m+SsJD9Nsvd0jSZ5YtvmOe3+thjw+5SkkeOjpiVpNL0N2KeqHgJ3jij32gd4KLAxzaNZ31pVD03yEeClwL8DhwN/VVUXJdkX+BSw3zTt/h3w11V1WpLNgVs7ej+SNGsYIEvS7PSjqloGLEtyE3Bcu/484EFtcPsY4GtJxups1Md+TwP+LcmRwDFVdVXH/ZakkWeALEmz020986t6llfR/Ns+D1gyNgLdr6p6f5LvAAcCpyd5SlX9posOS9JsYQ6yJI2mZcBa5/9W1VLgt0meB5DGg6erl2SPqjqvqv4VOBOYNm9ZkuYaA2RJGkFVdQNwWpJfJfngWu7mJcArk5wLnA8c1EedQ9s2zwVuAU5cy7YladbyNm+StB7zNm+SdHeOIEvS+m05cEg/DwoBrh1etyRp5jiCLEmSJPVwBFmSJEnqYYAsSZIk9TBAliRJknoYIEuSJEk9DJAlSZKkHgbIkiRJUg8DZEmSJKmHAbIkSZLUwwBZkiRJ6mGALKkTSRYmWZ5k/kz3ZW0lqSR7rmXdxye5sOs+9dHu/ZL8MsmyJG8cdvujLMlu7THdoKP9vTbJte33/B7ruK/Lkjyli35J6p4BsqS1Mv4/+Kq6oqo2r6o7ZrJfwzI+mK6qn1bV/WagK28BTqmqLarqY/1WSvLyJKcOsF9DN8igM8kC4N+Ap7bf8xs63PdhSf6rq/1JWncGyJI0u+0KnD+IHc/mXwMGYAdgYwb0WUsaLQbIktZYki8DC4Hj2p+b3zL+5+wkpyR5b5KftWWOS3KPJEcmWZrkF0l269nn3klOSnJjkguTPH+K9rdK8rkki5Nc3bYzP8lGSZYk2aen7PZJbklyz3b51Ukubts5NslOk7RxSpJX9SzfOeKa5Cft6nPb9/aCJE9KclVP+fu3+1iS5Pwkz+7Z9sUkn0zynTY14owke0zxfp/d7mNJu8/7t+t/CDwZ+ETbj/tOUPflSS5t2/ltkpe09T8DPLqtt6SnX59OckKSm4Ent5/ph5Jc0aYXfCbJJm357ZIc3/brxiQ/TTKv3fbW9tgsa4/n/lMcyy8l+X2Sy5O8o2cfL09yatv+H9r+HzDJfu72nezZ/JK2/9cn+ceeOvOSvC3JJUluSHJ0km0n2Pd9gbH0mSXt506Sjya5sv0+n5Xk8T11vpjkvT3Ld/l+9Kx/OvAPwAvafp870fuTNFwGyJLWWFUdDFwBPKv9ufkDkxR9IXAwsDOwB/Bz4AvAtsAFwLsAkmwGnAR8Bbgn8CLgU0keOMl+jwBWAnsCDwWeCryqqm4Djmnrj3k+8OOqui7JfsD72nU7ApcDX12L9/+EdvbB7fv/797taX6OPw74fvt+3gAcmaQ3BeNFwLuBbYCLgX+eqK02ODsKOBTYHjiBJgjcsKr2A34KvL7tx/+Nq7sZ8DHggKraAngMcE5VXQD8FfDztt7WPdVe3PZlC+BU4F+B+wIPofm8dwbe2ZZ9M3BV268daAK9at/n64FHtO0+Dbhsko/z48BWwO7AE4GXAq/o2b4vTXC6HfAB4HNJMn4n03wnHwfcD9gfeOfYHxjAG4HntO3uBPwB+OQE+/4/YOy7uHX7uQP8ov1ctqX57n4tycaTvM8JVdV3gX8B/rvt94PXpL6kwTBAljRIX6iqS6rqJuBE4JKq+kFVrQS+RhPcAjwTuKyqvlBVK6vqbOAbwHPH7zDJDsABwKFVdXNVXQd8hCYYhyZQ6Q2QX9yuA3gJ8PmqOrsNpt9OM4q6W3dvGYBHAZsD76+qFVX1Q+D4cf06pqr+p/0sjqQJtCbyAuA7VXVSVd0OfAjYhCbY7ccqYJ8km1TV4qqaLkXg21V1WlWtAm4DXg38TVXdWFXLaIK5sc/6dpo/NHatqtvbPOwC7gA2Ah6QZEFVXVZVl4xvKE0KxwuAt1fVsqq6DPgwzR9VYy6vqv9oc9uPaNvboc/3PubdVXVLVZ0LnAuMBaGvAf6xqq5qvw+HAc9Nnxf1VdV/VdUN7Xf2w+17nok8dEkdM0CWNEjX9szfMsHy5u38rsC+7U/1S9qf/F8C3GuCfe4KLAAW95T9LM1ILcAPgU2S7JtkV5rA85vttp1oRo0BqKrlwA00o6Jd2gm4sg0yx1w+rp3f9cz/kdWfxUT76u3zKuBK+uhzVd1ME4D+Fc3n9Z0ke09T7cqe+e2BTYGzej7r77brAT5IM/r9/TaN421tuxfTjHgfBlyX5KuZOJVlO2DD3vfHFJ9TVf2xnZ3ss5rMZJ/1rsA3e97bBTTBfV8BeJI3J7kgyU1t/a1o3pOkWc4AWdLaqg73dSVNGsTWPdPmVfXaScreBmzXU3bLqnog3BlAHk0zWvti4Ph25BPgGpqgCLgzBeEewNUTtHMzTXA4ZqJgfTLXALuM5dK2Fk7STj/76u1zgF363VdVfa+q/pRm5PU3wH+MbZqsSs/89TR/yDyw57Peqqo2b/e9rKreXFW7A88C/nYs17iqvlJVj2v7XjSpGuNdTzMKvWvPurX9nKZ6T5O5kib9pPd7t3FVTdt+m2/8Vpp0nW3aNJWbgLH0jzX5/nR5LknqgAGypLV1LU3eaBeOB+6b5OAkC9rpET25oneqqsU0ub0fTrJle6HVHkme2FPsKzQjpy9hdXrF2PpXJHlIko1o0gXOaH/aH+8c4M+TbJrmdm6vHLd9qvd/Bk2A9Jb2vTyJJoBc43xnmmD/GUn2b3Ob30zzB8LPpquYZIc0F/ht1tZZTnruLdsAACAASURBVDNCOtb/eyfZcLL67R8b/wF8JKsvctw5ydPa+Wcm2bMN2pe2+74jzb2Z92s/41tpguy73f6vTZs4GvjnJFu0I/5/C6ztLc/W9Dv5mbbtXdv3s32Sg/qsuwVNHvzvgQ2SvBPYsmf7OcCBSbZNci+aEfWp+r3buD+oJM2gTk7GJP/bx3RyF21JGhnvA97R/jz9d+uyo3aE96k0ua3X0Pwk/q80OZ0TeSnNT/O/prmw6us0I6Rj+xsLUHeiyX0eW38y8E80+c2LaS4cfCET+wiwgiZ4OYImT7jXYcAR7fu/yx03qmoF8GyaXOnrgU8BL62q30zS1qSq6kLgL2kuZrueJtB+VtvGdObRBNTXADfSXIz2unbbD2luWfa7JNdPsY+30qRRnJ5kKfADVufZ7tUuL6e5APNTVXUKzXF7f9vf39Gkv/zDJPt/A82xupTmosCvAJ/v471NZE2/kx8FjqVJEVkGnE5zUWA/vkfz3fo/mrSQW7lresqXafKdL6P5g+6/mdzX2tcbkpzdZ/uSBijN9RTruJPkfODAqYoAx1bVg9a5MUmSJGmAOnn8JvCaqrp8qgJJXjfVdkmSJGkUdDKCLEmSJM0VXY0gA5BkL5ocsAfQPJITgPYKZ0mSJGnkdX3F7BeAT9Nc2ftk4Es0FypIkiRJs0KnKRZJzqqqhyc5r6r+pF3306p6/HR1u7DddtvVbrvtNoymJEmSNMudddZZ11fV9uPXd5piAdza3sfxoiSvp7nZ+z2nqdOZ3XbbjTPPPHNYzUmSJGkWSzLhTSa6TrE4lObJQW8EHk5z786XddyGJEmSNDCdBchJtqd5XOYGVXVVVb2iqv6iqk7vqo2R873vwUEHwdKlM90TSZIkdaSrJ+m9iuaJTB8HfpPk2WtQ92+SnJ/kV0mOSrJx+2jOk5Jc1L5u00U/O7d4MRx7LNxww0z3RJIkSR3pagT5UOCBVfVo4DHA2/uplGRnmnSMRVW1DzCf5rGvbwNOrqq9gJPb5dGz7bbN6403zmw/JEmS1JmuAuQVVfV7gKq6FNhoDepuAGySZAOa/OVrgIOAI9rtRwDP6aif3TJAliRJmnO6uovFvZN8bLLlqnrjRJWq6uokHwKuAG4Bvl9V30+yQ1UtbsssTjLpnTCSHAIcArBw4cIO3soa2KbN/PjDH4bbriRJkgamqwD578ctn9VPpTa3+CDgPsAS4GtJ/nJNGq6qw4HDARYtWjTc52Y7gixJkjTndBIgV9UR05ea0FOA346lZyQ5hiaH+dokO7ajxzsC13XRz86NjSAbIEuSJM0ZnQTISY6dantVTXZXiyuARyXZlCbFYn/gTOBmmvsnv799/XYX/ezcxhvDJpuYYiFJkjSHdJVi8WjgSuAo4Awg/VSqqjOSfB04G1gJ/JImXWJz4Ogkr6QJop/XUT+7t+22jiBLkiTNIV0FyPcC/hR4EfBi4DvAUVV1/nQVq+pdwLvGrb6NZjR59BkgS5IkzSmd3Oatqu6oqu9W1cuARwEXA6ckeUMX+x9p22xjgCxJkjSHdDWCTJKNgGfQjCLvBnwMOKar/Y+sbbeFSy6Z6V5IkiSpI11dpHcEsA9wIvDuqvpVF/udFbbdFn7xi5nuhSRJkjrS1QjywTR3nrgv8Mbkzmv0AlRVbdlRO6PHHGRJkqQ5pav7IHf1yOrZZ5tt4JZb4NZbm9u+SZIkaVZbfwPbrow9Tc97IUuSJM0JnQTISc7uosys5OOmJUmS5pSucpDvn+R/p9geYKuO2hotY4+bdgRZkiRpTugqQN67jzJ3dNTWaHEEWZIkaU7p6iK9y7vYz6xkgCxJkjSneJHeujLFQpIkaU4xQF5XW24J8+Y5gixJkjRHdBogJ/nXftbNKfPmNaPIBsiSJElzQtcjyH86wboDOm5j9Gy7LVx//Uz3QpIkSR3o5CK9JK8FXgfsPu52b1sAp3XRxkjbbTf47W9nuheSJEnqQFe3efsKcCLwPuBtPeuXVdXczz3YYw84+uiZ7oUkSZI60FWKxXxgKfDXwLKeiSTbTlUxyf2SnNMzLU1yaJJtk5yU5KL2dZuO+tq9PfZocpC9k4UkSdKs19UI8llAtfMZt62A3SerWFUXAg8BSDIfuBr4Js1I9MlV9f4kb2uX39pRf7u1557N6yWXwKJFM9sXSZIkrZOuHhRyny72A+wPXFJVlyc5CHhSu/4I4BRGNUDeY4/m1QBZkiRp1utqBBmAJE+YaH1V/aTPXbwQOKqd36GqFrf1Fye55yRtHgIcArBw4cI163BXdm8HyC+5ZGbalyRJUmc6DZCBv++Z3xh4JE36xX7TVUyyIfBs4O1r0mBVHQ4cDrBo0aKapvhgbLYZ3OtecPHFM9K8JEmSutNpgFxVz+pdTrIL8IE+qx8AnF1V17bL1ybZsR093hG4rsOudm/PPR1BliRJmgMG/ajpq4B9+iz7IlanVwAcC7ysnX8Z8O0O+9W9PfYwQJYkSZoDus5B/jir72Yxj+buFOf2UW9TmqfwvaZn9fuBo5O8ErgCeF6Xfe3cHnvAEUfALbfAJpvMdG8kSZK0lrrOQT6zZ34lcFRVTfskvar6I3CPcetuoLmrxewwdqu3iy6CBz1oZvsiSZKktdZ1DvIRXe5vVnnMY5rXE0+cOkD+3e+aIHr5cli4EPbaCzbccDh9lCRJ0rQ6zUFOsleSryf5dZJLx6Yu2xhZu+4Kj340HHXUxNt/9jN4+tNhp53gCU+AAw+EffaBbbeFZzwDPvIROO88qJm5EYckSZIaXadYfAF4F/AR4MnAK7j7k/Xmrhe+EN70JrjgArj//Zt1K1bAP/4jfPjDsOOO8I53wOMfD5tuCpddBj//OfzgB3DCCU35HXaA/feHxz2uub/yVls1I8wbbQR33AG33Qa33nrX1975Kth882baYovV09jy5ps7Yi1JkjSFVIcjlknOqqqHJzmvqv6kXffTqnp8Z41MYdGiRXXmmWdOX3BQFi+GnXduguD3vAeuuQZe8AI49VR4zWvggx9sgtSJXHklnHxyEyz/4Adw7bUTl+vCggWrg+WxabPNmiB8o41g440nnp9uuZ+yG24I8wZ98xRJkqTptbHr3R6D3PUI8q1J5gEXJXk9cDUw4RPw5qQdd4SnPQ3e+174xS/gJz9pRnSPOqoZXZ7KLrvAy1/eTFVw9dXNCPOyZc0o9G23wfz5dw08JwpOAW6+uam3fPndXyebv/nm5rV3VHr8KPWqVd18TgsW9B94b7jh6mnBgonnB7Ftgw0g68+PH5IkabWuA+RDgU2BNwL/H02axcumrDHXHHUU/Mu/wOc/D895Dhx2WHMh3ppI4N73bqZRsnLlxGkda7M8Xdk//KF5XbECbr+9eR0/v2LFYHO2xwLmDTYY3LRgQX/l5s9fPc2bd9fXidZN99p12aSZT+46L0nSLNRJikWSL1fVwUneVFUf7aBfa2XGUyw0fHfcMXnwPFVg3e+2FSuaNlau7Ga6/fbpy3Q1Uj8KxgLmyQLo3vnptg+j7FhQP9VrP2W6qDNTdde0Tj/6LTuIfc62sjPd/qDKznT7Gn2PfSz82Z8NvdlBp1g8PMmuwP9L8iXGXZhXVTd21I50V/PnNw9mmUsPZ1m16u5B+e23N+vHtt1xx+r5fl8HVeeOO5qR/Kpm3fj5idaNallY/avERK/9lBn/uiZlZ7rumtbpR79lB7HP2VZ2ptsfVNmZbl+zQ9WMBMiT6SpA/gzwXWB34CzuGiBXu15SP+bNa6YFC2a6J5IkrZc6uZ1AVX2squ4PfL6qdq+q+/RMBseSJEmaNTq931ZVvbbL/UmSJEnD5g1pJUmSpB4GyJIkSVKPzgPkJLsmeUo7v0mSSR4dJ0mSJI2eTgPkJK8Gvg58tl11b+BbXbYhSZIkDVLXI8h/DTwWWApQVRexPj1qWpIkSbNe1wHybVW1YmwhyQY090GeUpKtk3w9yW+SXJDk0Um2TXJSkova12067qskSZJ0N10HyD9O8g/AJkn+FPgacFwf9T4KfLeq9gYeDFwAvA04uar2Ak5ulyVJkqSB6jpAfhvwe+A84DXACcA7pqqQZEvgCcDnAKpqRVUtAQ4CjmiLHQE8p+O+SpIkSXfT1aOmAaiqVcB/tFO/dqcJqr+Q5ME0j6p+E7BDVS1u97s4yYS5zEkOAQ4BWLhw4Tr0XpIkSer+LhaPbfOF/y/JpUl+m+TSaaptADwM+HRVPRS4mTVIp6iqw6tqUVUt2n777deh95IkSVLHI8g0aRJ/QzMKfEefda4CrqqqM9rlr9MEyNcm2bEdPd4RuK7jvkqSJEl303UO8k1VdWJVXVdVN4xNU1Woqt8BVya5X7tqf+DXwLHAy9p1LwO+3XFfJUmSpLvpZAQ5ycPa2R8l+SBwDHDb2PaqOnuaXbwBODLJhsClwCtogvejk7wSuAJ4Xhd9lSRJkqbSVYrFh8ctL+qZL2C/qSpX1Tnj6ozZfx37JUmSJK2RTgLkqnoyQJLdq+ouF+Ul2b2LNiRJkqRh6DoH+esTrPtax21IkiRJA9NVDvLewAOBrZL8ec+mLYGNu2hDkiRJGoaucpDvBzwT2Bp4Vs/6ZcCrO2pDkiRJGriucpC/DXw7yaOr6udd7FOSJEmaCZ3mIBscS5Ikabbr+iI9SZIkaVYzQJYkSZJ6dHWRHgBJNgL+Atitd99V9Z4u25EkSZIGpdMAGfg2cBNwFj2PmpYkSZJmi64D5HtX1dM73qckSZI0NF3nIP8syZ90vE9JkiRpaLoeQX4c8PIkv6VJsQhQVfWgjtuRJEmSBqLrAPmAjvcnSZIkDVXXDwq5nNWPm34WsHW7TpIkSZoVOg2Qk7wJOBK4Zzv9V5I39FHvsiTnJTknyZntum2TnJTkovZ1my77KkmSJE2k64v0XgnsW1XvrKp3Ao8CXt1n3SdX1UOqalG7/Dbg5KraCzi5XZYkSZIGqusAOcAdPct3tOvWxkHAEe38EcBz1qFfkiRJUl+6vkjvC8AZSb7ZLj8H+Fwf9Qr4fpICPltVhwM7VNVigKpanOSeHfdVkiRJuptOA+Sq+rckp9Dc7i3AK6rql31UfWxVXdMGwScl+U2/bSY5BDgEYOHChWvRa0mSJGm1rkeQqaqzgbPXsM417et17ejzI4Frk+zYjh7vCFw3Sd3DgcMBFi1aVOvUeUmSJK33us5BXmNJNkuyxdg88FTgV8CxwMvaYi8Dvj0zPZQkSdL6pPMR5LWwA/DNJND05ytV9d0kvwCOTvJK4ArgeTPYR0mSJK0nOg2Q2xHgW6pqVZL7AnsDJ1bV7ZPVqapLgQdPsP4GYP8u+ydJkiRNp+sUi58AGyfZmebexa8AvthxG5IkSdLAdH4f5Kr6I/DnwMer6s+AB3TchiRJkjQwnQfISR4NvAT4TrtuFPKcJUmSpL50HSAfCrwd+GZVnZ9kd+BHHbchSZIkDUzXDwr5MfBjgCTzgOur6o1dtiFJkiQNUqcjyEm+kmTL9m4WvwYuTPL3XbYhSZIkDVLXKRYPqKqlwHOAE4CFwMEdtyFJkiQNTNcB8oIkC2gC5G+39z/28c+SJEmaNboOkD8LXAZsBvwkya7A0o7bkCRJkgam64v0PgZ8rGfV5Ume3GUbkiRJ0iB1fZHeDkk+l+TEdvkBwMu6bEOSJEkapK5TLL4IfA/YqV3+P5p7I0uSJEmzQtcB8nZVdTSwCqCqVgJ3dNyGJEmSNDBdB8g3J7kH7Z0rkjwKuKnjNiRJkqSB6fQiPeBvgWOBPZKcBmwPPLfjNiRJkqSB6SxAbh8tvTHwROB+QIAL23shS5IkSbNCZykWVbUK+HBVrayq86vqV/0Gx0nmJ/llkuPb5W2TnJTkovZ1m676KUmSJE2l6xzk7yf5iyRZw3pvAi7oWX4bcHJV7QWc3C5LkiRJA9d1gPy3wNeA25IsTbIsyZRP0ktyb+AZwH/2rD4IOKKdP4Lm0dWSJEnSwHX9JL0t1qLavwNvAXrr7lBVi9t9Lk5yz8kqJzkEOKRdXJ7kwrXow7rYDrh+yG1qYh6L0eGxGB0ei9HhsRgdHovRMdPHYteJVnYaICc5uar2n25dz7ZnAtdV1VlJnrQ2bVbV4cDha1O3C0nOrKpFM9W+VvNYjA6PxejwWIwOj8Xo8FiMjlE9Fp0EyEk2BjYFtmsvqBvLQd6S1U/Vm8hjgWcnOZDmDhhbJvkv4NokO7ajxzsC13XRT0mSJGk6XeUgvwY4C9i7fR2bvg18crJKVfX2qrp3Ve0GvBD4YVX9Jc29lF/WFntZux9JkiRp4DoZQa6qjwIfTfKGqvp4B7t8P3B0klcCVwDP62CfgzJj6R26G4/F6PBYjA6PxejwWIwOj8XoGMljkarqbmfJ84DvVtWyJO8AHga8t6rO7qwRSZIkaYC6vs3bP7XB8eOAp9Hcou3THbchSZIkDUzXAfId7eszgE9X1beBDTtuQ5IkSRqYrgPkq5N8Fng+cEKSjQbQxoxI8vQkFya5OMndnuyXxsfa7f+b5GEz0c/1QR/H4klJbkpyTju9cyb6Odcl+XyS65L8apLtnhND0sex8JwYkiS7JPlRkguSnJ/kTROU8dwYgj6PhefGECTZOMn/JDm3PRbvnqDMSJ0Xnd4HmSYwfjrwoapa0t6i7e87bmPoksynuRvHnwJXAb9IcmxV/bqn2AHAXu20L01qyb7D7utc1+exAPhpVT1z6B1cv3wR+ATwpUm2e04MzxeZ+liA58SwrATeXFVnJ9kCOCvJSf5/MSP6ORbguTEMtwH7VdXyJAuAU5OcWFWn95QZqfOi09HdqvpjVR1TVRe1y4ur6vtdtjFDHglcXFWXVtUK4Ks0j8PudRDwpWqcDmzd/oGgbvVzLDQEVfUT4MYpinhODEkfx0JD0v6/d3Y7vwy4ANh5XDHPjSHo81hoCNrv+vJ2cUE7jb9LxEidF3Mi/WEIdgau7Fm+irufZP2U0brr93N+dPtTzolJHjicrmkcz4nR4jkxZEl2Ax4KnDFuk+fGkE1xLMBzYyiSzE9yDs3D306qqpE+L7pOsZirMsG68X/59FNG666fz/lsYNf2p5wDgW/R/GSj4fKcGB2eE0OWZHPgG8ChVbV0/OYJqnhuDMg0x8JzY0iq6g7gIUm2Br6ZZJ+q6r1uYqTOC0eQ+3MVsEvP8r2Ba9aijNbdtJ9zVS0d+ymnqk4AFiTZbnhdVMtzYkR4TgxXm2P5DeDIqjpmgiKeG0My3bHw3Bi+qloCnEJzzVqvkTovDJD78wtgryT3SbIhzWOxjx1X5ljgpe1VmI8CbqqqxcPu6Hpg2mOR5F5J0s4/kuZ7fsPQeyrPiRHhOTE87ef8OeCCqvq3SYp5bgxBP8fCc2M4kmzfjhyTZBPgKcBvxhUbqfPCFIs+VNXKJK8HvgfMBz5fVecn+at2+2eAE4ADgYuBPwKvmKn+zmV9HovnAq9NshK4BXhhdfnISAGQ5CjgScB2Sa4C3kVz4YXnxJD1cSw8J4bnscDBwHltviXAPwALwXNjyPo5Fp4bw7EjcER7J6p5wNFVdfwox1GdPmpakiRJmu1MsZAkSZJ6GCBLkiRJPQyQJUmSpB4GyJIkSVIPA2RJkiSphwGyJEmS1MMAWZJGUJKtk7yuZ3mnJF8fQDuHJbk6yXumKLNHknOSLO+6fUkaRd4HWZJGUJLdgOOrap8Bt3MYsLyqPtRH2eVVtfkg+yNJo8ARZEkaTe8HxkZuP5hktyS/Akjy8iTfSnJckt8meX2Sv03yyySnJ9m2LbdHku8mOSvJT5PsPV2jSZ7YtnlOu78tBvw+JWnk+KhpSRpNbwP2qaqHwJ0jyr32AR4KbEzzaNa3VtVDk3wEeCnw78DhwF9V1UVJ9gU+Bew3Tbt/B/x1VZ2WZHPg1o7ejyTNGgbIkjQ7/aiqlgHLktwEHNeuPw94UBvcPgb4WpKxOhv1sd/TgH9LciRwTFVd1XG/JWnkGSBL0ux0W8/8qp7lVTT/ts8DloyNQPerqt6f5DvAgcDpSZ5SVb/posOSNFuYgyxJo2kZsNb5v1W1FPhtkucBpPHg6eol2aOqzquqfwXOBKbNW5akucYAWZJGUFXdAJyW5FdJPriWu3kJ8Mok5wLnAwf1UefQts1zgVuAE9eybUmatbzNmyStx7zNmyTdnSPIkrR+Ww4c0s+DQoBrh9ctSZo5jiBLkiRJPRxBliRJknoYIEuSJEk9DJAlSZKkHgbIkiRJUg8DZEmSJKmHAbIkSZLUwwBZkiRJ6mGALEmSJPUwQJYkSZJ6GCBLWm8lWZhkeZL5M92XtZWkkuy5lnUfn+TCrvvUR7v3S/LLJMuSvHEN6r08yak9y8uT7D6YXkpanxkgS1pvJLksyVPGlqvqiqravKrumMl+Dcv4YLqqflpV95uBrrwFOKWqtqiqj63tTtpjd2mH/ZIkwABZkjR8uwLnz3QnJGkyBsiS1gtJvgwsBI5rf5p/S5Ld2lHVDdoypyR5b5KftWWOS3KPJEcmWZrkF0l269nn3klOSnJjkguTPH+K9rdK8rkki5Nc3bYzP8lGSZYk2aen7PZJbklyz3b51Ukubts5NslOk7RxSpJX9SzfmZKQ5Cft6nPb9/aCJE9KclVP+fu3+1iS5Pwkz+7Z9sUkn0zynTY14owke0zxfp/d7mNJu8/7t+t/CDwZ+ETbj/tOUPflSS5t2/ltkpdM0sadI+Jt/z7THo9lSX6cZNfJ+idJUzFAlrReqKqDgSuAZ7U/zX9gkqIvBA4Gdgb2AH4OfAHYFrgAeBdAks2Ak4CvAPcEXgR8KskDJ9nvEcBKYE/gocBTgVdV1W3AMW39Mc8HflxV1yXZD3hfu25H4HLgq2vx/p/Qzj64ff//3bs9yQLgOOD77ft5A3Bkkt4UjBcB7wa2AS4G/nmittqg9yjg0P+/vTuPrrO+7zz+/nqR9wVs412WF9mWDAmLxhhCEjBLsONCwhYDweCSOCTDOZMzyUwznbRlkp6GnKRJSUkgnpYG5pCkWSChhITSNj0hmQPFNjAgyfsq29jYeJMty1j6zh+/e3meu0i6lnXvcyV9Xufo6C6/e+9XenjMx1//nt8PmAA8R/iLSYW7LwZeBO5P1bEx67UjgO8AS9x9FHA58FqBP+adwFeB8anXPFng60REMiggi4hk+gd33+LuR4BfA1vc/V/c/TTwU0K4BVgGbHf3f3D30+6+Dvg5cEv2G5rZRGAJ8Hl3P+7u+4FvE8I4hJAdD8h3pB6DEPoec/d1qTD9P4DL4p3sHrIIGAk86O6n3P3fgGez6nrK3f8j9bt4Eriwg/f6BPArd3/B3d8FvgkMI4TdQrQD55vZMHff6+6FTsf4lbv/LvV7+p+E39P0Al8rIvIeBWQRkUz7Yrdb8twfmbo9A7g0NYXgsJkdJoTZSXnecwYwGNgbG/t9QqcW4N+AYWZ2aWpawIXA06nnphC6xgC4ezNwkNDh7klTgF3u3h57bEfW57wVu32C6HeR773iNbcDuyigZnc/TgjY9xF+X78ys/kF/QThM9Lv0wy8k6pFROSMDEq6ABGREvIefK9dhGkQ1xY4thUYn+q+Zhbl3m5mPyF0a/cBz7r7sdTTewgBG3hvCsI4YHeezzkODI/dzxfWO7IHmG5mA2IhuRLY2MlrOnuvC9J3zMyA6eSvOYe7Pw88b2bDgL8E/jfwwQJe+l632MxGEqbF7Cm8bBGRQB1kEelP9gE9tW7us8BcM7vLzAanvv5T+mK0OHffS5jb+9dmNtrMBpjZbDP7cGzYDwmd0zuJplekH19pZhea2RDgr4CX3X17nppeA24ys+Gpi9fuzXq+s5//ZULA/u+pn+VK4I/oxnxn4CfAR83s6tTc5i8Q/oLwf7t6oZlNTF3gNyL1mmag0GX4lprZFWZWQZiL/LK77+rqRSIi2RSQRaQ/+Rrw5dQ0hy+ezRulOrzXEeYR7yFMP/g6MKSDl6wAKoAG4BDwM8JFd+n3SwfUKYS5z+nH/xX4M8L85r2ECweXk9+3gVOEIPw4uRepPQA8nvr5M1bccPdTwA2EudIHgO8BK9x9fQef1SF33wB8Evjb1Hv9EeHiyFMFvHwAIVDvIUyR+DDwuQI/+oeEiyjfAS4h/GVDROSMmXtP/oujiIhI6ZnZD4Amd/9y0rWISO+nDrKIiIiISIwCsoiIiIhIjKZYiIiIiIjEqIMsIiIiIhLTp9ZBHj9+vFdVVSVdhoiIiIj0AmvXrj3g7hOyH+9TAbmqqoo1a9YkXYaIiIiI9AJmtiPf45piISIiIiISo4AsIiIiIhLTp6ZYiIiIiEh5O3UKNm+G+npoaAjfr74aPvOZpCuLKCCLiIiISI87dQo2boxCcPr7pk1w+nQYYwazZkFdXbK1ZlNAFhEREZFua20NQTgeghsaQhBuawtjBgwIQXjBAvjYx8L32lqYNw+GD0+2/nyKFpDN7DFgGbDf3c9PPXYu8I9AFbAduM3dD+V57fXAQ8BA4O/c/cFi1SkiIiIiXTt5EjZsyAzBDQ1hukQ8CM+ZE8LvzTeH7wsWwNy5MGxYsvWfiWJ2kH8APAw8EXvsS8C/uvuDZval1P0/ib/IzAYC3wWuBZqAV8zsGXdvKGKtIiIiIgK0tIQgnN0R3rIF2tvDmIEDQxBesABuvTXqCM+dC0OHJlt/TyhaQHb335lZVdbDNwJXpm4/Dvw7WQEZWAhsdvetAGb249TrFJBFREREesiJE7B+fW5HeOvWKAgPGgTV1fC+98Htt0cd4epqGDIk2fqLqdRzkCe6+14Ad99rZuflGTMV2BW73wRc2tEbmtkqYBVAZWVlD5YqIiIi0vsdP54bhOvrYds2cA9jBg0K3d+LLoI778wMwhUVydafhHK8SM/yPOYdDXb31cBqgLq6ug7HiYiIiPRlzc3Q2Bh1gtNhePv2KAgPHhwujKurg7vvDkG4tjYE4cGDEy2/9ymsLwAAGQ9JREFUrJQ6IO8zs8mp7vFkYH+eMU3A9Nj9acCeklQnIiIiUuaOHYuCcLwjvCO2aXJFRQjCCxfCypVRR3j2bAXhQpQ6ID8D3A08mPr+yzxjXgGqzWwmsBtYDtxRsgpFREREysDRoyEIZ18st3NnNGbIEJg/Hy6/HD71qehiudmzw7QJ6Z5iLvP2I8IFeePNrAn4C0Iw/omZ3QvsBG5NjZ1CWM5tqbufNrP7gecJy7w95u71xapTREREJElHjuROi6ivh6amaMzQoSEIX3FFFIIXLICZMxWEi8Hc+8603bq6Ol+zZk3SZYiIiIjkOHw4NwQ3NMDu3dGYoUOhpiYzBNfWhiA8cGBytfdVZrbW3XP28dPfOURERER60KFDuSG4vh727o3GDB8egvDixVEIrq2FqioF4XKggCwiIiLSDQcP5q4hXF8Pb70VjRkxIgTh667L7AjPmBF2nZPypIAsIiIi0okDB3I7wg0NsG9fNGbkyBB8r78+syNcWakg3BspIIuIiIgA+/fn7wi//XY0ZtSoEHw/+tEoBC9YANOng+XbyUF6JQVkERER6Tfcc4Nw+vuBA9G40aND8L3hhsypEdOmKQj3BwrIIiIi0ue4h7nA+ZZPe+edaNyYMSH8fvzjmR3hKVMUhPszBWQRERHptdzD6hD5OsKHDkXjxo4NwfeWWzI7wpMnKwhLLgVkERERKXvusGdP/ovlDh+Oxp17bgi/t92WuZbwxIkKwlI4BWQREREpG+5hB7l8HeGjR6Nx48aF4Hv77Zkd4fPOUxCWs6eALCIiIiXnDrt25e8IHzsWjZswIYTfT34ysyM8YUJytUvfp4AsIiIiRdPeDjt35naEGxuhuTkaN3FiCL933x1dLFdbqyAsyVBAFhERkbPW3g47duR2hBsb4fjxaNykSaEDvHJl5oYa48YlV7tINgVkERERKVh7O2zbljstorERTpyIxk2ZEoLvpz6V2RE+99zkahcplAKyiIiI5GhrC0E4e1e59euhpSUaN3Vq6ASvWhXND66pgXPOSa52kbNV8oBsZvOAf4w9NAv4c3f/m9iYK4FfAttSDz3l7l8pWZEiIiL9RFsbbNmSu6HG+vVw8mQ0bvr0EICvvDJzasSYMYmVLlI0JQ/I7r4BuBDAzAYCu4Gn8wx90d2XlbI2ERGRvur06SgIx+cJb9gAra3RuMrKEICvvjqzIzx6dHK1i5Ra0lMsrga2uPuOhOsQERHpE959NwTh7IvlNmyAU6eicTNmhPB73XVRR7imBkaNSq52kXKRdEBeDvyog+cuM7PXgT3AF929Pt8gM1sFrAKorKwsSpEiIiLl5t13YdOm3I7wxo3hubSZM0P4XbIk6gjPnw8jRyZXu0i5M3dP5oPNKgjhd4G778t6bjTQ7u7NZrYUeMjdq7t6z7q6Ol+zZk1xChYREUnAqVMhCGd3hDduDNMmIOwcN3Nm5kYatbUhCI8YkWz9IuXMzNa6e13240l2kJcA67LDMYC7H43dfs7Mvmdm4939QEkrFBERKZHW1hB6s5dP27QpMwjPnh3C7403RmF43jwYPjzZ+kX6kg4DspkVslJhu7sf7uZn304H0yvMbBKwz93dzBYCA4CD3fwcERGRsnHyZAjC2R3hzZvDihIAAwaEILxgAXz841FHeN48GDYs2fpF+oPOOsh7Ul/WyZiBwBlP/DWz4cC1wGdij90H4O6PArcAnzWz00ALsNyTmgsiIiLSDSdPhqXSsjvCmzeHzTYABg6EOXNC+L311mjptHnzYOjQZOsX6c86C8iN7n5RZy82s1e786HufgIYl/XYo7HbDwMPd+e9RURESqmlJQTh7A01tm7NDMLV1XDBBfCJT0RTI+bOhSFDkq1fRHJ1FpAvK+D1hYwRERHp9U6cCNspZ2+osXUrpP+Nc9CgEHovvBDuuCOaGjF3LlRUJFu/iBSuw4Ds7u/tn2Nm5wDT4+PdfV18jIiISF/Q3Jw7NaK+HrZvj4Lw4MEh9F58Mdx1V9QRnjNHQVikL+hyFQsz+ypwD7AFSM8DdmBx8coSEREprubm0BHOvlhu+/ZoTEVFmA+8cCHcc0/UEZ4zJ4RkEembClnm7TZgtruf6nKkiIhImTl2LHdaRH097NwZjamoCGsGL1oE994bdYRnzw7TJkSkfynktH8TGAvsL3ItIiIi3XbkSP6O8K5d0ZghQ0IQvuKKzA01Zs1SEBaRSCF/HHwNeNXM3gRa0w+6+w1Fq0pERKQDhw/n7wjv3h2NGToUamrgQx/K3F1u5sywooSISGcKCciPA18H3gDai1uOiIhIcOhQbghuaIA9e6Ixw4aFILx4cWZHuKpKQVhEuq+QgHzA3b9T9EpERKRfeued3BDc0AB790Zjhg8Pwfeaa6IQnA7CAwYkVrqI9FGFBOS1ZvY14Bkyp1isK1pVIiLS5xw4kBuC6+th375ozIgRIfh+5CNRCF6wACorFYRFpHQKCcjp3fQWxR7TMm8iIpLX22/nhuCGBtgfu9R71KgQfpcuzZwaMX26grCIJK/LgOzuV5WiEBER6T3cQ+DNd7HcgQPRuNGjQ/BdtizzYrlp08AsufpFRDrTYUA2s2Xu/mxnLy5kjIiI9F7uYQpEvovlDh6Mxo0ZE8Lvxz6W2RGeOlVBWER6n846yN8ws91AZ3+0/RWggCwi0su5h4vi8nWEDx2Kxo0dG8LvTTdldoQnT1YQFpG+o7OAvA/4Vhev39SDtYiISJG5h2XS8nWEDx+Oxp1zTgi+t92W2RGeNElBWET6vg4DsrtfWawPNbPtwDGgDTjt7nVZzxvwELAUOAHco1UzREQK5w5NTbmrRjQ0hB3n0saNC+F3+fLM5dMmTlQQFpH+K8mNNa9y9wMdPLcEqE59XQo8kvouIiIx7e2wc2fu1IiGBmhujsZNmBCC7x13ZE6NmDBBQVhEJFu57jx/I/CEuzvwkpmNNbPJ7r63qxeKiPRFbW2wfXtmJ7i+Hhob4cSJaNykSSH8rlwZdYNrakIQFhGRwiQVkB34ZzNz4Pvuvjrr+anArtj9ptRjOQHZzFYBqwAqKyuLU62ISIm0tcGWLZlBuKEB1q+HlpZo3NSpIfx++tNRN7imBs49N7naRUT6ii4DspkNBT4HXEEItr8HHnH3k2fxuR9w9z1mdh7wgpmtd/ffxT82z2s83xulwvVqgLq6urxjRETKzbvvRkE4Pi1iwwZobY3GVVaGAHzVVVFHuLY2LKsmIiLFUUgH+QnCBXV/m7p/O/B/gFu7+6Huvif1fb+ZPQ0sBOIBuQmYHrs/DdjT3c8TEUnKqVOwaVPuHOGNG0NITquqyt1iuaYm7DgnIiKlVUhAnufu74/d/62Zvd7dDzSzEcAAdz+Wun0d8JWsYc8A95vZjwkX5x3R/GMRKWetraH7mz1HeNOmMG0CwsVws2bl7iw3fz6MGJFs/SIiEikkIL9qZovc/SUAM7sU+MNZfOZE4OmwkhuDgB+6+2/M7D4Ad38UeI6wxNtmwjJvK8/i80REekxLSwjC2Uunbd4cVpQAGDAA5swJ4femm6I5wvPmwbBhydYvIiJds7BQRCcDzBqBecDO1EOVQCPQDri7v6+oFZ6Buro6X7NmTdJliEgfcPx4uDAue47w1q1hjWGAgQOhujpz/eDaWpg7F4YOTbZ+ERHpmpmtzd6PAwrrIF9fhHpERMrCsWNhqbTsqRHbt0djBg8Ooffii+Guu6IgXF0NFRWJlS4iIkXSYUA2s9HufpRwgV4Od3+naFWJiPSwI0dyl05raAibbKRVVIT5wIsWwR//cdQZnj07hGQREekfOusg/xBYBqwlLLEWX3rNgVlFrEtEpFsOHcq/q9zu3dGYoUPDChEf/GDm1IhZs2BQuW6fJCIiJdPh/wrcfVnq+8zSlSMiUpgDB3KnRTQ0wFtvRWOGDw9B+OqrM4NwVVWYPywiIpJPZ1MsLu7she6+rufLERGJuMP+/bnTIurr4e23o3EjR4bge/31mRfMVVaGFSVERETORGf/mPjXnTznwOIerkVE+in30PnNXjqtoQEOHozGjR4dAvANN0RLp9XWwrRpYY1hERGRntDZFIurSlmIiPR97mEucL45wocPR+POOSeE35tvjrrBCxbA5MkKwiIiUnxdXo5iZrcCv0ntfPdl4GLgq+7+atGrE5Feqb0ddu3KP0f4WGxdnHHjQvC9/fbMOcITJyoIi4hIcgq5XvvP3P2nZnYF8BHgm8CjhC2gRaQfa2+HHTvyT404fjwaN3FiCL4rVmTOEZ4wIbnaRUREOlJIQG5Lff8o8Ii7/9LMHiheSSJSbk6fhi1bQvBNb6rR2Bi+WlqicVOmhOB7773RtIiamtApFhER6S0KCci7zez7wDXA181sCKDrwkX6oJMnYePGzBDc0BAee/fdaNz06SEAf+hDUUe4pibMHRYREentCgnItxG2m/6mux82s8nAfytuWSJSTM3NsH59bkd4y5YwbQLC8mizZoXgu2xZFILnz4dRo5KtX0REpJi6DMjufgJ4KnZ/L7C3mEWJSM9I7yqX3RGOb688eDBUV8P73w/Ll0fzg+fODTvOiYiI9DfaVFWkl8veTCMehuO7yqW3V77iiqgbXFsLs2eHkCwiIiJByQOymU0HngAmAe3Aand/KGvMlcAvgW2ph55y96+Usk6RcuMelk5LB+B4GD50KBo3alQIvkuWRCG4pgZmzND2yiIiIoVIooN8GviCu68zs1HAWjN7wd0bssa96O7LEqhPJFFtbbBtW2YAbmgIc4abm6Nx48aF8HvbbVEQrq0NK0loDWEREZHuK3lAjs9hTm0+0ghMBbIDskifduoUbNqU2xHesAFaW6Nx6aXTVq7MnBqhNYRFRESKI9E5yGZWBVwEvJzn6cvM7HVgD/BFd6/v4D1WAasAKisri1OoyFk4cSKE3uyO8ObNoVucVlUVgu9112VOjRgzJrHSRURE+qXEArKZjQR+Dnze3Y9mPb0OmOHuzWa2FPgFUJ3vfdx9NbAaoK6uzotYskinjhyJNs+Ih+Ht28P8YQhzgOfMCcH35pujEDxvHowYkWj5IiIikpJIQDazwYRw/KS7P5X9fDwwu/tzZvY9Mxvv7gdKWadItvSKEY2NYU5wPAzv3h2Nq6gIoXfhQrj77mh+8Jw5MGRIcvWLiIhI15JYxcKAvwca3f1bHYyZBOxzdzezhYSd+w6WsEzp59raQuc3HoTTt+MrRowYETrAixdnzg+eORMGaRFFERGRXimJ/4V/ALgLeMPMXks99qdAJYC7PwrcAnzWzE4DLcByd9f0CelxLS3R1srxILxxY+aFcuedF8JvesWI+fPD92nTwo5zIiIi0ncksYrF74FOF6Fy94eBh0tTkfQH77yT2QVO347PDzYLnd+amuhCuXQYPvfcRMsXERGREtI/Akuf0d4OTU35g/Dbb0fjhgyJ5gevWBEF4epqGDYsufpFRESkPCggS69z6lRYIi07CK9fH5ZUSzvnnBB8b7ghc1qEdpQTERGRziggS9k6ejT3ArnGRtiyJXP94MrKEH4//enMIDxhgnaUExERkTOngCyJcoe9e/MH4T17onGDB4cpEBdcALfeGk2LmDsXRo5Mrn4RERHpexSQpSRaW8O0iA0bQgBOf1+/PnSK00aNCsH32mszu8EzZ4aQLCIiIlJsCsjSY9KbaKQDcDwMb9sWLqJLmzo1XCh3112ZQXjyZE2LEBERkWQpIMsZy9cNTt8+ciQaN2xYmAJxySVw550hEM+fr2kRIiIiUt4UkCWvM+0Gz5+fGYLnzYPp07WJhoiIiPQ+Csj9XKHd4KFDQ+i95BK4444oBM+dG+YNi4iIiPQVCsj9QLobnH2BXGdzg9Pd4HRHWN1gERER6S8UkPuQlpawRvDGjblhOLsbnJ4brG6wiIiISCYF5F6mrQ127AghOPtr587QLU6bMiWE33gInjcvbKyhbrCIiIhIfgrIZcgd9u3LH4K3bAlbLaeNHh06vx/4AKxcGW5XV4cgrG6wiIiIyJlTQE7QkSOwaVP+IHzsWDSuoiKE3vnz4YYbQghOB+HzztO6wSIiIiI9KZGAbGbXAw8BA4G/c/cHs5631PNLgRPAPe6+ruSF9oDW1mhecPbXvn3RODOoqgrB9/LLoxA8d264QG7gwMR+BBEREZF+peQB2cwGAt8FrgWagFfM7Bl3b4gNWwJUp74uBR5JfS9LbW2wa1f+ELxjR+YqERMnhtC7bFlmCJ41K1w8JyIiIiLJSqKDvBDY7O5bAczsx8CNQDwg3wg84e4OvGRmY81ssrvvLX25HfvZz+CBB8I6wq2t0eOjRoXQu2gRrFiROSVizJjEyhURERGRAiQRkKcCu2L3m8jtDucbMxXICchmtgpYBVBZWdmjhXZlzBiYMweWLs3sBk+cqHnBIiIiIr1VEgE5X3T0bowJD7qvBlYD1NXV5R1TLNdeG75EREREpO9IYjXcJmB67P40YE83xoiIiIiI9LgkAvIrQLWZzTSzCmA58EzWmGeAFRYsAo6U2/xjEREREembSj7Fwt1Pm9n9wPOEZd4ec/d6M7sv9fyjwHOEJd42E5Z5W1nqOkVERESkfzL3kk7bLSozexvYUeKPHQ8cKPFnSn46FuVDx6J86FiUDx2L8qFjUT6SPhYz3H1C9oN9KiAnwczWuHtd0nWIjkU50bEoHzoW5UPHonzoWJSPcj0WScxBFhEREREpWwrIIiIiIiIxCshnb3XSBch7dCzKh45F+dCxKB86FuVDx6J8lOWx0BxkEREREZEYdZBFRERERGIUkEVEREREYhSQC2Rm15vZBjPbbGZfyvO8mdl3Us//PzO7OIk6+4MCjsWVZnbEzF5Lff15EnX2dWb2mJntN7M3O3he50SJFHAsdE6UiJlNN7PfmlmjmdWb2X/JM0bnRgkUeCx0bpSAmQ01s/8ws9dTx+J/5RlTVudFyXfS643MbCDwXeBaoAl4xcyecfeG2LAlQHXq61LgkdR36UEFHguAF919WckL7F9+ADwMPNHB8zonSucHdH4sQOdEqZwGvuDu68xsFLDWzF7Q/y8SUcixAJ0bpdAKLHb3ZjMbDPzezH7t7i/FxpTVeaEOcmEWApvdfau7nwJ+DNyYNeZG4AkPXgLGmtnkUhfaDxRyLKQE3P13wDudDNE5USIFHAspEXff6+7rUrePAY3A1KxhOjdKoMBjISWQ+m+9OXV3cOore5WIsjovFJALMxXYFbvfRO5JVsgYOXuF/p4vS/1Tzq/NbEFpSpMsOifKi86JEjOzKuAi4OWsp3RulFgnxwJ0bpSEmQ00s9eA/cAL7l7W54WmWBTG8jyW/TefQsbI2Svk97yOsLd6s5ktBX5B+CcbKS2dE+VD50SJmdlI4OfA5939aPbTeV6ic6NIujgWOjdKxN3bgAvNbCzwtJmd7+7x6ybK6rxQB7kwTcD02P1pwJ5ujJGz1+Xv2d2Ppv8px92fAwab2fjSlSgpOifKhM6J0krNsfw58KS7P5VniM6NEunqWOjcKD13Pwz8O3B91lNldV4oIBfmFaDazGaaWQWwHHgma8wzwIrUVZiLgCPuvrfUhfYDXR4LM5tkZpa6vZDw3/nBklcqOifKhM6J0kn9nv8eaHT3b3UwTOdGCRRyLHRulIaZTUh1jjGzYcA1wPqsYWV1XmiKRQHc/bSZ3Q88DwwEHnP3ejO7L/X8o8BzwFJgM3ACWJlUvX1ZgcfiFuCzZnYaaAGWu7aM7HFm9iPgSmC8mTUBf0G48ELnRIkVcCx0TpTOB4C7gDdS8y0B/hSoBJ0bJVbIsdC5URqTgcdTK1ENAH7i7s+Wc47SVtMiIiIiIjGaYiEiIiIiEqOALCIiIiISo4AsIiIiIhKjgCwiIiIiEqOALCIiIiISo4AsIiIiIhKjgCwiUobMbKyZfS52f4qZ/awIn/OAme02s690Mma2mb1mZs09/fkiIuVI6yCLiJQhM6sCnnX384v8OQ8Aze7+zQLGNrv7yGLWIyJSDtRBFhEpTw8C6c7tN8ysyszeBDCze8zsF2b2T2a2zczuN7P/amavmtlLZnZuatxsM/uNma01sxfNbH5XH2pmH0595mup9xtV5J9TRKTsaKtpEZHy9CXgfHe/EN7rKMedD1wEDCVszfon7n6RmX0bWAH8DbAauM/dN5nZpcD3gMVdfO4Xgf/s7n8ws5HAyR76eUREeg0FZBGR3um37n4MOGZmR4B/Sj3+BvC+VLi9HPipmaVfM6SA9/0D8C0zexJ4yt2berhuEZGyp4AsItI7tcZut8futxP+bB8AHE53oAvl7g+a2a+ApcBLZnaNu6/viYJFRHoLzUEWESlPx4Buz/9196PANjO7FcCC93f1OjOb7e5vuPvXgTVAl/OWRUT6GgVkEZEy5O4HgT+Y2Ztm9o1uvs2dwL1m9jpQD9xYwGs+n/rM14EW4Nfd/GwRkV5Ly7yJiPRjWuZNRCSXOsgiIv1bM7CqkI1CgH2lK0tEJDnqIIuIiIiIxKiDLCIiIiISo4AsIiIiIhKjgCwiIiIiEqOALCIiIiIS8/8BS4uZkOSgl+EAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x504 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plotting time evolution of slip-rate, stress on the fault and slip\n",
    "plt.figure(figsize=(10,7))\n",
    "plt.subplot(3,1,1)\n",
    "plt.plot(T,Vd,'r-')\n",
    "plt.title('time evolution of slip rate')\n",
    "plt.ylabel('slip-rate [m/s]')\n",
    "plt.xlabel('time [s]')\n",
    "plt.ylim([0, 8])\n",
    "\n",
    "plt.subplot(3,1,2)\n",
    "plt.plot(T, Sd, 'r-')\n",
    "plt.title('time evolution of stress on the fault')\n",
    "plt.ylabel('stress on the fault [MPa]')\n",
    "plt.xlabel('time [s]')\n",
    "plt.ylim([40, 83])\n",
    "\n",
    "plt.subplot(3,1,3)\n",
    "plt.plot(T, Slip, 'b-')\n",
    "plt.xlabel('time [s]')\n",
    "plt.ylabel('slip [m]')\n",
    "plt.title('time evolution of slip')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "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
}