{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Boundary conditions tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial aims to demonstrate how users can implement various boundary conditions in Devito, building on concepts introduced in previous tutorials. Over the course of this notebook we will go over the implementation of both free surface boundary conditions and perfectly-matched layers (PMLs) in the context of the first-order acoustic wave equation. This tutorial is based on a simplified version of the method outlined in Liu and Tao's 1997 paper (https://doi.org/10.1121/1.419657).\n",
    "\n",
    "We will set up our domain with PMLs along the left, right, and bottom edges, and free surface boundaries at the top as shown below.\n",
    "\n",
    "<img src=\"figures/boundary_conditions.png\" style=\"width: 220px;\"/>\n",
    "\n",
    "Note that whilst in practice we would want the PML tapers to overlap in the corners, this requires additional subdomains. As such, they are omitted for simplicity.\n",
    "\n",
    "As always, we will begin by specifying some parameters for our `Grid`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "shape = (101, 101)\n",
    "extent = (2000., 2000.)\n",
    "nbpml = 10  # Number of PMLs on each side"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will need to use subdomains to accomodate the modified equations in the PML regions. As `Grid` objects cannot have subdomains added retroactively, we must define our subdomains beforehand."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from devito import SubDomain\n",
    "\n",
    "s_o = 6  # Space order\n",
    "\n",
    "class MainDomain(SubDomain):  # Main section with no damping\n",
    "    name = 'main'\n",
    "    def __init__(self, PMLS, S_O):\n",
    "        super().__init__()\n",
    "        self.PMLS = PMLS\n",
    "        self.S_O = S_O\n",
    "\n",
    "    def define(self, dimensions):\n",
    "        x, y = dimensions\n",
    "        return {x: ('middle', self.PMLS, self.PMLS),\n",
    "                y: ('middle', self.S_O//2, self.PMLS)}\n",
    "\n",
    "\n",
    "class Left(SubDomain):  # Left PML region\n",
    "    name = 'left'\n",
    "    def __init__(self, PMLS):\n",
    "        super().__init__()\n",
    "        self.PMLS = PMLS\n",
    "\n",
    "    def define(self, dimensions):\n",
    "        x, y = dimensions\n",
    "        return {x: ('left', self.PMLS), y: y}\n",
    "\n",
    "\n",
    "class Right(SubDomain):  # Right PML region\n",
    "    name = 'right'\n",
    "    def __init__(self, PMLS):\n",
    "        super().__init__()\n",
    "        self.PMLS = PMLS\n",
    "\n",
    "    def define(self, dimensions):\n",
    "        x, y = dimensions\n",
    "        return {x: ('right', self.PMLS), y: y}\n",
    "    \n",
    "class Base(SubDomain):  # Base PML region\n",
    "    name = 'base'\n",
    "    def __init__(self, PMLS):\n",
    "        super().__init__()\n",
    "        self.PMLS = PMLS\n",
    "\n",
    "    def define(self, dimensions):\n",
    "        x, y = dimensions\n",
    "        return {x: ('middle', self.PMLS, self.PMLS), y: ('right', self.PMLS)}\n",
    "\n",
    "class FreeSurface(SubDomain):  # Free surface region\n",
    "    name = 'freesurface'\n",
    "    def __init__(self, PMLS, S_O):\n",
    "        super().__init__()\n",
    "        self.PMLS = PMLS\n",
    "        self.S_O = S_O\n",
    "        \n",
    "    def define(self, dimensions):\n",
    "        x, y = dimensions\n",
    "        return {x: ('middle', self.PMLS, self.PMLS), y: ('left', self.S_O//2)}\n",
    "\n",
    "main_domain = MainDomain(nbpml, s_o)\n",
    "left = Left(nbpml)\n",
    "right = Right(nbpml)\n",
    "base = Base(nbpml)\n",
    "freesurface = FreeSurface(nbpml, s_o)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And create the grid, attaching our subdomains:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from devito import Grid\n",
    "\n",
    "grid = Grid(shape=shape, extent=extent,\n",
    "            subdomains=(main_domain, left, right, base, freesurface))\n",
    "x, y = grid.dimensions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then begin to specify our problem starting with some parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "density = 1.  # 1000kg/m^3\n",
    "velocity = 4.  # km/s\n",
    "gamma = 0.0002  # Absorption coefficient"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need a `TimeFunction` object for each of our wavefields. As particle velocity is a vector, we will choose a `VectorTimeFunction` object to encapsulate it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from devito import TimeFunction, VectorTimeFunction, NODE\n",
    "\n",
    "p = TimeFunction(name='p', grid=grid, time_order=1,\n",
    "                 space_order=s_o, staggered=NODE)\n",
    "v = VectorTimeFunction(name='v', grid=grid, time_order=1,\n",
    "                       space_order=s_o)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A `VectorTimeFunction` is near identical in function to a standard `TimeFunction`, albeit with a field for each grid dimension. The fields associated with each component can be accessed as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[0. 0. 0. ... 0. 0. 0.]\n",
      "  [0. 0. 0. ... 0. 0. 0.]\n",
      "  [0. 0. 0. ... 0. 0. 0.]\n",
      "  ...\n",
      "  [0. 0. 0. ... 0. 0. 0.]\n",
      "  [0. 0. 0. ... 0. 0. 0.]\n",
      "  [0. 0. 0. ... 0. 0. 0.]]\n",
      "\n",
      " [[0. 0. 0. ... 0. 0. 0.]\n",
      "  [0. 0. 0. ... 0. 0. 0.]\n",
      "  [0. 0. 0. ... 0. 0. 0.]\n",
      "  ...\n",
      "  [0. 0. 0. ... 0. 0. 0.]\n",
      "  [0. 0. 0. ... 0. 0. 0.]\n",
      "  [0. 0. 0. ... 0. 0. 0.]]]\n"
     ]
    }
   ],
   "source": [
    "print(v[0].data)  # Print the data attribute associated with the x component of v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may have also noticed the keyword `staggered` in the arguments when we created these functions. As one might expect, these are used for specifying where derivatives should be evaluated relative to the grid, as required for implementing formulations such as the first-order acoustic wave equation or P-SV elastic. Passing a function `staggered=NODE` specifies that its derivatives should be evaluated at the node. One can also pass `staggered=x` or `staggered=y` to stagger the grid by half a spacing in those respective directions. Additionally, a tuple of dimensions can be passed to stagger in multiple directions (e.g. `staggered=(x, y)`). `VectorTimeFunction` objects have their associated grids staggered by default.\n",
    "\n",
    "We will also need to define a field for integrating pressure over time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_i = TimeFunction(name='p_i', grid=grid, time_order=1,\n",
    "                   space_order=1, staggered=NODE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we prepare the source term:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from examples.seismic import TimeAxis, RickerSource\n",
    "\n",
    "t0 = 0.  # Simulation starts at t=0\n",
    "tn = 400.  # Simulation length in ms\n",
    "dt = 1e2*(1. / np.sqrt(2.)) / 60.  # Time step\n",
    "\n",
    "time_range = TimeAxis(start=t0, stop=tn, step=dt)\n",
    "\n",
    "f0 = 0.02\n",
    "src = RickerSource(name='src', grid=grid, f0=f0,\n",
    "                   npoint=1, time_range=time_range)\n",
    "\n",
    "# Position source centrally in all dimensions\n",
    "src.coordinates.data[0, :] = 1000."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "tags": [
     "nbval-skip"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAg0AAAF9CAYAAACH0lvIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAABA4UlEQVR4nO3deZxkZXn3/8/V+zrTs+PAbA0DMqCCGRFFBdGIUR9IHhKNioKJookxmoRE+ZnHGESDG/ozGgMmccMYI8EHogaRPSJoRjYZtplhFoYZmJ6t1+qu6urr+eOc01NTU119urvWM9/361Wv6j51zqn7TPVUXXXd133f5u6IiIiITKeh2g0QERGR+qCgQURERGJR0CAiIiKxKGgQERGRWBQ0iIiISCwKGkRERCSWpmo3oNYtXrzYV69eXe1miIiIVMSvfvWrve6+pNBjChqmsXr1ajZs2FDtZoiIiFSEmW2f6jF1T4iIiEgsChpEREQkFgUNIiIiEouCBhEREYlFQYOIiIjEoqBBREREYqmLoMHMjjOzvzeze81sxMzczFbHPLbBzC43s21mNmpmD5nZhWVusoiISOLURdAAnAC8GTgA/PcMj/0E8HHgy8BvAfcB3zezN5SygSIiIklXL5M73e3uywDM7N3A6+IcZGZLgcuAq9z9c+HmO8zsBOAq4MflaKyIiEgS1UWmwd0nZnnoeUALcF3e9uuAF5jZmjk1TERE5ChSF0HDHJwCjAGb87ZvDO/XVbY5IiIi9SvpQcNC4KC7e972/TmPywxlJ5xHnuknO5H/zyoiIkmW9KBhVszsUjPbYGYb+vr6qt2cmvODB57hTX//M879/J3sOpiqdnNERKRCkh40HAB6zMzytkcZhv0U4O7Xuvt6d1+/ZEnB1UGPavfvOEBzo7F93wi3PfZctZsjIiIVkvSgYSPQChyftz2qZXi0ss1Jho27Bli/aiHdrU08+dxQtZsjIiIVkvSg4WYgA7w9b/tFwCPuvrXyTapv49kJHt89wCnL57F2WReb9gxWu0kiIlIh9TJPA2b2u+GPvxHe/5aZ9QF97n5XuM848E13/0MAd99jZlcDl5vZIHA/8BbgXOD8il5AQmzpG2ZsfIJTjp3H4Og4t6p7QkTkqFE3QQPw/bzf/yG8vws4J/y5Mbzl+igwBHwQOAZ4Anizu/+wPM1Mto27+gE4dfl89g2l+d6Gp9k3NMairtYqt0xERMqtboIGd88vZoy1j7tngSvDm8zRY7sHaG1qoHdJF7v7RwHYtGdIQYOIyFEg6TUNUmK7+kdZ3tNOY4OxdlkXAJueU12DiMjRQEGDzEjf4BhLwqzCMfPaaG1qYMf+kSq3SkREKkFBg8zI3sExlswLggYzY3FXK/uG0lVulYiIVIKCBpmR3EwDwOLuVvqGxqrYIhERqRQFDRJbKp1lcGycpfMOBQ1LulrYq0yDiMhRQUGDxNY3GGQUcjMNizpb2adMg4jIUUFBg8TWNxQMsVzSnds90cK+4TQTWvFSRCTxFDRIbHsGgozC0u62yW2Lu1rJTjgHU5lqNUtERCpEQYPEFhU85mYaokmd1EUhIpJ8Choktr7BMRobjIWdLZPbFncFP2sEhYhI8ilokNj6BsdY1NlCY8Oh2bqjokiNoBARST4FDRLbnsGxw7omQN0TIiJHEwUNEluh1Sx72ptpbDD2KmgQEUk8BQ0S28DoOD3tzYdta2gwFnW2sHdQ3RMiIkmnoEFi609lmNd+5Grqi7pa2TesTIOISNIpaJBY3J3+VIb5eZkGCLooBlLjVWiViIhUkoIGiWU4nSU74QWDhnntTfRrcicRkcRT0CCxDIRBwby2AkFDWzMDowoaRESSTkGDxBJlEgpnGpongwoREUkuBQ0SS9Ggoa2Z4XSW8exEpZslIiIVpKBBYomChnlT1DQADI6qGFJEJMkUNEgsA9NkGgDVNYiIJJyCBomleKYhDBo07FJEJNEUNEgsA6kMZtDdeuTkTvPagm3KNIiIJJuCBomlP5VhXlszDTkrXEbmd0SZBgUNIiJJpqBBYhkYHS84hTSopkFE5GihoEFimWoKaVBNg4jI0UJBg8RSLGjobGmkwdBU0iIiCaegQWIpFjSYWTArpLonREQSrS6CBjNbYWbXm1m/mQ2Y2Q1mtjLmsSvN7JtmtsPMUmb2pJldaWad5W53kgyEhZBTmdemqaRFRJKucGVbDTGzDuB2YAy4GHDgSuAOM3uhuw8XObYTuBVoBv4PsAN4CfC3wFrgLeVtfXIMjGYKztEQmdfexIBmhBQRSbSaDxqA9wC9wEnuvhnAzB4GNgHvBa4ucuxZBMHBee5+S7jtDjNbCFxmZh3uPlK+pidDJjvBaGaCrgJzNESUaRARSb566J44H7gvChgA3H0rcA9wwTTHtoT3A3nbDxJc+5GTDsgRhseCDEKxoGG+ahpERBKvHoKGU4BHCmzfCKyb5thbCTISnzazdWbWZWbnAh8E/rFY14YcMhQjaOhua9LoCRGRhKuHoGEhcKDA9v3AgmIHuvso8AqC69wIDAK3AT8E/mSq48zsUjPbYGYb+vr6ZtvuxBgeywLQWSRo6GptntxPRESSqR6Chlkzszbge8BS4B3A2cBfEhRAfmWq49z9Wndf7+7rlyxZUpG21rKhsSCD0NU2ddDQ2drIcHocd69Us0REpMLqoRDyAIUzClNlIHL9IXAOcIK7bwm33W1m/cC1ZvaP7v5QyVqaUENhBqGrtXHKfTpbm3CHkXS2aEZCRETqVz1kGjYS1DXkWwc8Os2xLwAO5AQMkV+G9yfPsW1HhagQslgwED0W7SsiIslTD0HDTcCZZtYbbTCz1QTDKW+a5thngQVmdkLe9peG98+UqpFJNhTOv9DZUqymIchCDCloEBFJrHoIGr4GbANuNLMLzOx84EbgaeCaaCczW2Vm42b2sZxjv0FQ/PhjM7vYzF5tZn8JfA74FcGwTZlGFAh0F6tpaIkyDSqGFBFJqpoPGsJhkecCTwLfBr4DbAXOdfehnF0NaCTnmtx9G3Am8CDBLJI/Jpgs6lrgN919ovxXUP/idE9EwzGVaRARSa66qFhz9x3AhdPss40CkzW5+6PAm8vTsqPDUHqclqYGmhunjjGjgGIkraBBRCSpaj7TINU3NDpO9zQjIjqVaRARSTwFDTKt4bHxaYdRdrWqpkFEJOkUNMi0hsamn3uhMxw9oSGXIiLJpaBBpjU0lik6sRMcGj2h7gkRkeRS0CDTGh7LFl2sCqChwehoaVSmQUQkwRQ0yLTi1DRAUAw5rNETIiKJpaBBpjU4Nj5tpgGCYsghFUKKiCSWggaZ1nDMoKGzVd0TIiJJpqBBispOeOyVKztbmlQIKSKSYAoapKioRiFu94QyDSIiyaWgQYqKs+5EpFNBg4hIoilokKIOBQ3F52kI9lEhpIhIkilokKKiaaGjyZuK6dQ8DSIiiaagQYoaSQdBQ0fMTEMqkyU74eVuloiIVIGCBikqWuq6I0amYXLRKk3wJCKSSAoapKgo09DZEi/TADCiugYRkURS0CBFRZmG9hhBQ0e4z4gyDSIiiaSgQYo6lGmYvnuifTJoUKZBRCSJFDRIUTMphGxvDvZJZRQ0iIgkkYIGKWokPU5jg9HSOP2fSocyDSIiiaagQYoaHsvS0dKImU27b9Q9kVJNg4hIIilokKJS6exkBmE60bBMZRpERJJJQYMUNZwej1UECeqeEBFJOgUNUtRIOhtruCXkdk8oaBARSSIFDVLUyEwyDc3KNIiIJJmCBilqJpmGpsYGWhobGMmoEFJEJIkUNEhRI+lsrGWxI+0tjeqeEBFJKAUNUtTI2DjtzfG6JyAohlT3hIhIMtVN0GBmK8zsejPrN7MBM7vBzFbO4PiTzez7ZrbXzFJm9oSZfbCcbU6CkYwyDSIiEoj/FbKKzKwDuB0YAy4GHLgSuMPMXujuw9Mcvz48/k7g3UA/sBboKmOzE2FkLBtrWexIR0ujppEWEUmouggagPcAvcBJ7r4ZwMweBjYB7wWunupAM2sAvgXc5u6/k/PQHeVrbjJkshOksxOxJ3cC6Ghu0iqXIiIJVS/dE+cD90UBA4C7bwXuAS6Y5thzgJMpElhIYZOLVc0gaFD3hIhIctVL0HAK8EiB7RuBddMc+4rwvs3M7jOzjJntMbMvmVl7SVuZMKnJoEGFkCIiUj9Bw0LgQIHt+4EF0xy7PLz/HnAL8JvAZwhqG/610AFmdqmZbTCzDX19fbNrcQIMh90MMy2EVNAgIpJM9VLTMBdRYHSdu38s/PlOM2sErjKzk939sdwD3P1a4FqA9evXe+WaWltGxoIP//bmGQQNzSqEFBFJqnrJNBygcEZhqgxErn3h/U/ztt8S3p8+h3Yl2shkpmGm3RMqhBQRSaJ6CRo2EtQ15FsHPBrj2GImZtWio0DUzRB3Gulg3yZGMxNMTBy1CRoRkcSql6DhJuBMM+uNNpjZauCs8LFi/otgfofz8ra/PrzfUKI2Js5sRk9E+6qLQkQkeeolaPgasA240cwuMLPzgRuBp4Frop3MbJWZjZtZVLuAu+8D/g54n5l9ysxea2YfAT4GfDN3GKccbrJ7YoajJ4JjFTSIiCRNXRRCuvuwmZ0LfAH4NmDAbcCH3H0oZ1cDGjkyGLoCGAT+GLgM2A18FvhEmZte16JsQdsMCyEBzdUgIpJAdRE0ALj7DuDCafbZRhA45G93gsmdNMHTDKRm1T0R/ElpeWwRkeSpl+4JqYLZZBrUPSEiklwKGmRKqXSW1qYGGhuOSN5MKRppoe4JEZHkUdAgU0plsjMabgnKNIiIJJmCBpnSSDpLxwy6JkBDLkVEkkxBg0wplcnSNsNMQ3tYCJnSrJAiIomjoEGmlEpnZzRyApjMTKh7QkQkeRQ0yJRS6eyMFquCQ4WQChpERJJHQYNMaSSTnexuiKu1qYEG0+gJEZEkUtAgUxpNZ2lvntmfiJnR3tyoTIOISAIpaJApjWTGJ2d4nIn2liZSmhFSRCRxFDTIlFLpiRnNBhnpaFGmQUQkiRQ0yJRS6fEZF0KCggYRkaRS0CAFuTupzMyHXEIwgkKFkCIiyaOgQQoaG59gwpnxNNIQZRpU0yAikjQKGqSg0XAa6Nl0T7Q3N6l7QkQkgRQ0SEHR2hGzzTRo7QkRkeRR0CAFRZmC2dQ0qBBSRCSZFDRIQVEh42yGXKoQUkQkmRQ0SEFR98JsMw2pTBZ3L3WzRESkihQ0SEFRpmB28zQ0kZ1w0tmJUjdLRESqSEGDFBTVJMymEDIKNNRFISKSLLGDBgucb2afM7Ovm9mqcPvZZra8fE2UapjLkMsOLY8tIpJIsVYjMrMFwI+BlwKDQBfw98B24D3AfuBPy9RGqYJDoydms2CVggYRkSSKm2n4LLACOAtYBFjOY7cCrylxu6TKUnPKNASBhronRESSJe7XyAuAy9z9XjPL/xTZQRBQSIKkwmmg51LToKmkRUSSJW6moQt4ZorH2jg88yAJkMpkaWwwmhtn/tJOdk9oVkgRkUSJGzQ8AbxuisfOBn5dmuZIrRhJZ2lvbsRs5kFDVAip7gkRkWSJ2z3xD8CXzawf+NdwW4+ZvQv4E+DScjROqmc0k51V1wRo9ISISFLFyjS4+7XA1cDfApvDzT8FrgW+6O7fKU/zDjGzFWZ2vZn1m9mAmd1gZitncZ6PmJmb2c/K0c6kiDINs9E+mWlQTYOISJLEHk/n7h8xs68CvwksBfYBP3X3p8rVuIiZdQC3A2PAxYADVwJ3mNkL3X045nl6gb8G9pSrrUmRSmdnNYU0HBo9oUyDiEiyzGgQvrtvB/6pTG0p5j1AL3CSu28GMLOHgU3AewmyIHF8FfgOcBIzvPajTSqTndViVQAd0YyQKoQUEUmUKT84Z5r6d/cdc2/OlM4H7osChvD5tprZPQTDQacNGszsbcCLgbcCN5SroUkxl0xDQ4PR2tSgQkgRkYQp9m17G0E3QFyz+4SJ5xTgxgLbNwK/N93B4YyWXwD+yt33z2ZEwNEmlckyv7151se3hytdiohIchQLGv6AQ0FDK0EtwADw78BzwDHAm4Fu4BNlbCPAQuBAge37gQUxjv8s8CTwjThPZmaXEo4IWblyxrWWiZBKz370BAQTPCnTICKSLFMGDe7+jehnM/sicD/wO+7uOduvAP4vsK5sLZwjM3sl8E7gxbltLyYcLXItwPr162eSbUmMVGb2oycgyDRocicRkWSJO7nTW4Fr8j90w9//EXhbqRuW5wCFMwpTZSByXQP8M7DTzHrMrIcgWGoMf28taUsTYmQONQ0QZBpGlWkQEUmUuCMIuoAlUzy2FOgsTXOmtJGgriHfOuDRaY49Oby9r8BjB4A/A744l8YlUSqTpW0OQUNHS6OGXIqIJEzcoOFO4FNm9pi7/0+00czOAD4ZPl5ONwGfM7PeaF4IM1tNsOrmR6Y59tUFtn2RoHDzAxyarEpC2QknPT5BR/PsR6W2NTcyOKrJnUREkiTup8KfECyBfZ+ZPU1QCLmMYHXLreHj5fS18DluNLO/JijQ/ATwNEH3AwBmtgrYAlzh7lcAuPud+Sczs4NAU6HHJGdZ7Ja4vVdHam9upG9wrFRNEhGRGhB3GumtwPMJUvy3EcwGeRvBxEonu/u2cjUwfP5h4FyCERDfJpigaStwrrsP5exqBBmE2X/ayeSS1u0ts880qHtCRCR5ZjKNdIbgG//Xytecos+/A7hwmn22EWOZbnc/pzStSqbR9ATAnEdPaJ4GEZFk0TdyOcJIJsg0zG30RJPmaRARSZhYmQYz20rx2SHd3Y8vTZOk2qIP+7llGhpIZbK4O5qBU0QkGeJ2T9zFkUHDIuDlwBDBCpSSEFHQMNsFqyAIOLITTjo7QWtTOWcYFxGRSokVNLj7JYW2hxMl3UwwskISIqpFmFP3RFhEOZpW0CAikhRzqmlw94ME6zp8rCStkZoQjXqY69oToOWxRUSSpBSFkKPAcSU4j9SIyXka5tA9EWUpouGbIiJS/2Y9EN/MmoBTgY8TTPMsCTGamXumoU2ZBhGRxIk7emKCqUdPDABvLFmLpOqi7om51DREx2rYpYhIcsTNNFzBkUHDKLAd+C937y9pq6SqJkdPzKGAMcpSKNMgIpIccUdPfLzM7ZAakspkaWtuoKFh9vMrTBZCKtMgIpIYsQohzex2M3v+FI+daGaapyFBUunsnIogQZkGEZEkijt64hxg3hSPdQNnl6Q1UhNG0lk65rBYFSjTICKSRDMZcjlVIeTxBLNCSkKMht0Tc3FoyKWCBhGRpJjy66SZvQt4V/irA9ea2WDebu0Ewy5vK0/zpBpG0uNzzjRoyKWISPIU+zo5AWTDm+X9Ht32AV8F/rC8zZRKSmXmXtPQ2tSAmbonRESSZMqvk+7+TeCbAGZ2B/BH7v54pRom1ZNKZ+npaJnTOcyMjuZGZRpERBIk7pDLV5e7IVI7RtJZlvfMfZGp9hYFDSIiSVKspuGdwI/cfV/4c1Hu/q2StkyqZqQEQy4hDBrUPSEikhjFMg3fAM4kqFv4xjTncUBBQ0KMZrJzWnci0t6soEFEJEmKBQ1rgN05P8tRIpinoRSZhiZG1D0hIpIYxQohtxf6WZJtYsKD0RNzHHIJ0N7cwKgyDSIiiTG3GXwkcUbHw2WxS1HToNETIiKJUqwQcitTzwKZz939+NI0SaqpFMtiRzpamhhJj8z5PCIiUhuK5aDvIn7QIAkRFS6WohCyrbmR0czEnM8jIiK1oVhNwyUVbIfUiKg7oTSZhkZG0uNzPo+IiNQG1TTIYaLuiZLN06CaBhGRxIgdNJjZWjP7ppk9aWbD4f03zOyEcjZQKivKDJSye2JiQr1cIiJJECtoMLNzgIeANwH3Af8Q3v8v4NdmdnaZ2icVNjrZPTH3IZdRF0c0IkNEROpb3EzD54EHgFXu/k53/0t3fyewGngwfLyszGyFmV1vZv1mNmBmN5jZyhjHrTeza83scTMbMbMdZvYdM9OEVQWUcvRE1MWhWSFFRJIhbtCwDvi0uw/lbnT3QeDTwCmlblguM+sAbgeeD1wMvANYC9xhZp3THP77Yfu+BPwW8BHgxcAGM1tRtkbXqVLXNOSeU0RE6lvcHPROYKq1kluAZ0rTnCm9B+gFTnL3zQBm9jCwCXgvcHWRYz/t7n25G8zsHmBreN6PlaXFdaqUQy6jwGNUxZAiIokQN9PwaeBvzWx57kYzOxb4G+BTpW5YnvOB+6KAAcDdtwL3ABcUOzA/YAi3bQf6gGNL3M66V+ohl6BMg4hIUsTNNJwNzAOeMrP7gOeAZQSrYD4HnBMWS0IwO+TFJW7nKcCNBbZvBH5vpiczs5OBpcBjc2xX4kQf8G1NJaxpUKZBRCQR4gYNrwDGCVa9XBXe4NAqmK/M2bcc4+sWAgcKbN8PLJjJicysCfhHgkzDP0+xz6XApQArV05ba5koqfQ4bc0NNDTYnM/V1qKgQUQkSWIFDe6epJEGXwZeDrzR3QsFIrj7tcC1AOvXrz+qJhkIlsWe+3BLONQ9odETIiLJUJpPh/I7QOGMwlQZiILM7CqCDMLF7n5LidqWKKlMtiQjJ0BDLkVEkmZGQUM4RHEF0Jb/mLvfXqpGFbCRwsM61wGPxjmBmX0U+DDwAXf/dgnbliipdLYkRZCQM+RS3RMiIokQK2gws17gO8AZ0abw3sOfHSjNJ01hNwGfM7Ned38qbNNq4CyCeReKMrM/Ba4EPuruXy5jO+veSCmDhmjIpTINIiKJEDfT8E/ASuBDwONAulwNmsLXgD8BbjSzvyYIUj4BPA1cE+1kZquALcAV7n5FuO33gS8CNwO3m9mZOecdcPdYmYqjRSqdpa3E3RMacikikgxxg4aXAJe4+3+UszFTcfdhMzsX+ALwbYLsxm3Ah/JmqTSCjEfu/BOvD7e/Przlugs4p0zNrkupTJbFXVPN4zUzTY0NtDQ2aPSEiEhCzGRGyEpnFw7j7juAC6fZZxuHuk6ibZcAl5SrXUkzkh6no6WjZOdra27QjJAiIgkRd0bITwEfjrHOg9S5VDpbkimkIx0tTZPLbYuISH2LO0/Dt83s+cC2cEbI/GGO5ZgFUqpgpIRDLiEYQZHKTJTsfCIiUj1xR09cAlwOZAlWiMzvqjiqJkBKslKOnoCgGDKlTIOISCLErWn4W+AHwB+6+8HyNUeqKTvhpMcnSto9EWQaVNMgIpIEcWsaFgH/oIAh2Uq5wmUkyDQoaBARSYK4QcPPgJPL2RCpvqhgsdQ1DZqnQUQkGeJ2T3wQ+HczO0AwSdIR6z24u6rd6lyUEWgv0YJVEGYa1D0hIpIIcT8dHgvvv1Vkn3JOIy0VUI7uiQ5lGkREEiNu0HAFGiGReCOTmYYSF0IqaBARSYS48zR8fKrHzOwc4J2laY5U02T3RAlrGjrDyZ3cHTOb/gAREalZcQshD2NmJ5jZFWa2lWANiDeXtllSDVGmoaSjJ1oamXAYG1fJi4hIvYsdNJjZfDO71MzuAZ4APkpQEPnHwPIytU8qqFw1DYC6KEREEqBo0GBmDWb2BjP7HrAb+EdgFfCVcJcPufs17j5Q5nZKBUQzN5Zy9EQUNIxoBIWISN2b8tPBzD4PvA1YCowSzAj5TeBWYB7wJ5VooFTOSBlqGqIAZGRMU0mLiNS7Yl8p/4xgxMSPgUvcfV/0gJlpJEUClaOmoTPKNKh7QkSk7hXrnvhnYBB4I/CEmX3ZzM6oTLOkGkYzWcygtWlW9bEFtStoEBFJjCk/Hdz9PcAxwNuBDcB7gXvN7DHgw2jehsQZSWfpaG4s6dDIjrB7IpVR94SISL0r+pXS3Ufd/bvu/npgJYeWx/4IYMBVZnaRmbWVv6lSbiPpbEmLICGnEFKZBhGRuhc7D+3uu939M+5+KnAGwQiKtQRTS+8uU/ukglLpcdpbStc1ATlBw5iCBhGRejerTwh33+DuHyCYn+FC4M5SNkqqI+ieKHWmIRw9kVb3hIhIvZvTJ4S7ZwiGYv6gNM2RakplsiVddwI0T4OISJKUNhctdS2VzpZ0uCUEIzHMNCOkiEgSKGiQSSPpbEkndgIwMzqatTy2iEgSKGiQSeXongDoaG1STYOISAIoaJBJ5eiegKCuQZkGEZH6p6BBJo2kxydHO5RSu7onREQSQUGDTEplsrSVuKYBgkyDCiFFROqfggYBIJOdIJP1snRPdKqmQUQkEeomaDCzFWZ2vZn1m9mAmd1gZitjHttmZp81s91mljKze83sVeVucz1JZUq/wmVE3RMiIslQF0GDmXUAtwPPBy4G3kEwhfUdZtYZ4xT/DLwH+BjwJoJpr39iZqeVpcF1KOo+KMvoCRVCiogkQumr3srjPUAvcJK7bwYws4eBTQSrb1491YFm9iLgbcAfuPvXw213ARuBK4Dzy9v0+hB9qJcl09DSpKBBRCQB6iLTQPDBfl8UMAC4+1bgHuCCGMdmgO/lHDsO/Btwnpm1lr659Wd4LKg56CzD6ImgEFI1DSIi9a5eMg2nADcW2L4R+L0Yx25195ECx7YAJ4Q/l922vcP0DY2xrLuNpfNayzJSYbYmg4bW0v9JdLY0MpLJ4u6YWcnPP1sj6XF2HkgxkMowODrOwGiGscwEWXcm3JmYcLITzoQT/O5e7SaLiBxhUWcrF/7GcRV5rnoJGhYCBwps3w8smMOx0eOHMbNLgUsBVq6MVWsZy/W/2smX75hMltC7pJPTjuvhzN5FnHfqMcxvby7Zc81Uubsn3GE0M1GWmom4dh4Y4Y7H93DXk308vLOfPYNjVWuLiEiprHvePAUN1eTu1wLXAqxfv75kXy/f9tKVnLFmIc8OjLLrYIpHnhng7k17ueGBZ/jYTY/w1jNW8qHXnMj8jsoHD8PpMmYaWsOVLtPjVQkanuob4tM3P84tjz6HO6xY2M6rTlzCmsWdrFjYwYKOZrpam+hua6KtuZHGBqPRDDM79HMDNJhRO3kSEZFAQwUzuPUSNBygcEZhqixC/rGrpjgWDmUcym55TzvLe9oP2+buPLyzn2/du51v3budHz68my/9/um87PhFlWoWACNjQaahHEFDNMvk8FiWRV0lP/2U3J1v/nwbn/rx47Q2NfD+c07gwt84jtWLOmqqm0REpF7USyHkRoLahHzrgEdjHLsmHLaZf2wa2HzkIZVjZrxoRQ+ff/OLuPH9ZzG/vZl3/ssv+OHDuyrajqHJQsjSZwK6wkzDcAWLId2dT/7oMT7+n4/yyrWLue2ys7nsvJNYs7hTAYOIyCzVS9BwE3CmmfVGG8xsNXBW+Fgx/wk0k1MwaWZNwFuAW9y9Zjq2Tz12Pv/xRy/n9BUL+LPvPcjPN++t2HNHMzaWY+2JQ5mGygUNX7ljM//0s61c8vLVfO2d61na3Vax5xYRSap6CRq+BmwDbjSzC8zsfILRFE8D10Q7mdkqMxs3s49F29z9AYLhll80s3eb2WsIhluuAf6mgtcQy/z2Zr528XrWLO7kA999gL1DlYlphtNZmhuNlqbS/0lEXR7DFZqr4edb9vL5nz7Jb5+2nI+9aR0NDcosiIiUQl0EDe4+DJwLPAl8G/gOsBU4192HcnY1oJEjr+tdwNeBK4EfASuA17v7/WVu+qzMb2/my297MYNj4/x/N/y6Is85MjZelnoGOFQIWYlMw0h6nL/494dYs6iTT/7OCxQwiIiUUL0UQuLuO4ALp9lnGxxZ4O7uKeDPw1tdOHFZNx967Vo+c/MT3LN5L2edsLiszzeczpZlYic4NGFUJYKGa+56it39o1z/vpeVLQgSETla1UWm4Wj1B2et4diedv7uvx5jYqK8EwsNj42XZY4GyOmeKHPQsGdwlGvu3sIbX/A81q8+YvoNERGZIwUNNaytuZEPvXYtjzwzwF2b+sr6XMPpLB1l+mYeBSPlrmn41s+3MzY+wWXnnVTW5xEROVopaKhxF5x2LEu7W/nGPdvK+jwjY+OTQyNLrbWpgaYGK2umYSQ9znW/2M7r1i1jzeI4C5+KiMhMKWiocS1NDVx05iruerKPLX1D0x8wS8PpbFmGW0IwF0W5l8e+6cFdHBzJ8O5X9k6/s4iIzIqChjrw1jNW0thg3HD/zrI9x/DYeFkmdop0tTZNTiBVDjfc/wzHL+lk/arpliIREZHZUtBQB5Z0t3LWCYu58cFdeJlWWhxJj5etpgGgo7VpcgKpUnt6/wi/3Laf//3i4zTbo4hIGSloqBMXvGg5Ow+kuH/HdEttzM7wWLasmYbO1iaGxsrTPXHjg88AcMFpy8tyfhERCShoqBOvO2UZLU0N/OjhZ0t+7uyEk8pkyzqvQWdLIyNl6p645dHnOH1lD8ctyF9eRERESklBQ53obmvmZb2LuOOJPSU/dyoTrnBZpkJICNafKEdNw56BUR7e2c9rT15W8nOLiMjhFDTUkXOfv5Ste4fZune4pOeNhkJ2lGnIJQQrXZZj9MTtjwdB1GtOXlryc4uIyOEUNNSRV58UfDDe8Xhpsw3Dk8ti118h5G2P7+HYnnZOWtZd8nOLiMjhFDTUkZWLOuhd0smdT5Z2dsgoA1DOmoZyDLnMTjj3PbWPV524WKMmREQqQEFDnXn58Yv41bb9jGcnSnbOQ5mG8nVPdLQ0MpqZIFvCNTQe3TXA4Og4Z/YuKtk5RURkagoa6sxL1yxiOJ1l466Bkp0zyjSUc56GrmjRqhJ2Udz31D4AXqagQUSkIhQ01JmX9garN/5i676SnXOwIpmGIGgYKeFcDfc+tY/eJZ0snddWsnOKiMjUFDTUmaXdbfQu7uQXT+0v2TmHRoOgobutuWTnzNcZjswoVV1DdsL5n637eekaZRlERCpFQUMdOmPNQv5n2/6STSk9NJYBoLutnJM7hZmGEnVPbOkbYnBsXGtNiIhUkIKGOnT6yh4GRsdLNl/D0Og4ZkGxYrlEIzOirMZcPRBOp33ayp6SnE9ERKanoKEOvWhFDwAP7TxYkvMNjo3T1dpU1mGLURZjsETdEw8+fZD57c2sWdRZkvOJiMj0FDTUobVLu+loaeShp/tLcr6h0XG6yzhyAg6NnihdpuEgL1rRQ0OD5mcQEakUBQ11qLHBOPXY+Tz49MGSnG9wdJyuMtYzQE6mYTQz53MNj43z5HODnBZmXEREpDIUNNSp01f08OiuAcbG5z6EcSjsniinKCgpxeiJjbsGmHA4bcX8OZ9LRETiU9BQp049dj7p7ASbnhua87kGx8bpKuNwS4DWpkZamhpKUtPw6K6gW+aU5QoaREQqSUFDnVq3fB4Aj+2e+8yQQ6OZstc0AMxra2KwBDUNG3cNsKizhaXdrSVolYiIxKWgoU6tXtRJe3Mjj5YiaKhA9wSEi1aVKGhYt3yeFqkSEakwBQ11qrHBeP7zunm0BGtQDI2Ol3Vip0h3W/OcCyHT4xNs2jM4mWkREZHKUdBQx9Y9bx6P7R6Y08yQ2QlnOJ0t++gJKM3y2Jv2DJLJuuoZRESqQEFDHTv5efMYGB3nmYOpWZ8j+hCvRPdEdwlqGh7fPQjAuud1l6JJIiIyA3URNJhZg5ldbmbbzGzUzB4yswtjHDfPzD5mZj83s31mdjD8+bcr0OyyO/l5UTHk4KzPEQUNleie6CpB0LClb4imBmOVZoIUEam4uggagE8AHwe+DPwWcB/wfTN7wzTHrQT+GLgLuAh4C/Ak8AMze3/ZWlsha5d1AbB5z+yHXUaFiV2t5R1yCTCvBDUNW/qGWLWog+bGevnTFRFJjvJ/vZwjM1sKXAZc5e6fCzffYWYnAFcBPy5y+Fag191Hcrb9xMxWAB8GvlKONlfKvLZmls1rnVvQEK5wWcmaBnef9ciHLX3DHL+kq8QtExGROOrh69p5QAtwXd7264AXmNmaqQ509+G8gCGyAVheuiZWz9ql3WzeM/vuicHRytY0TDiMpGc3i2UmO8H2fcMcv1RBg4hINdRD0HAKMAZsztu+MbxfN4tzvgp4fC6NqhUnLO1i856hWY+giIKGStU0wOynkn56/wiZrNO7WPUMIiLVUA9Bw0LgoB/5qbg/5/HYzOxS4Ezg74rtY2YbzGxDX1/fjBpbaScs7WI4nWV3/+isjq/s6ImgbmK2dQ1b+oYBlGkQEamSigcNZvZaM/MYtzvL8NznAF8CvuXu35lqP3e/1t3Xu/v6JUuWlLoZJXXC0rkVQw5VMNMQTVU92xEUW/qCazx+sYIGEZFqqEYh5M+Bk2PsF9UiHAB6zMzysg1RhmE/MZjZS4CbgNuBd8dsa81bGwYNm/YM8aoTZx7gRAtIdbZUpqYBZh80PNU3xOKuVuZ3lH+kh4iIHKniQUNYmDiTeoKNQCtwPIfXNUS1DI9OdwIzewHwE+BB4EJ3n9u4vxqyqKuVBR3Ns840DKQydLc10dBQ/nUc5lrTEIycUD2DiEi11ENNw81ABnh73vaLgEfcfWuxg81sLfBT4CngTe4+++kTa9RcRlAcHEkzv70y39znUtPg7mzeM6R6BhGRKqr5eRrcfY+ZXQ1cbmaDwP0EkzSdC5yfu6+Z3QascvcTwt+XEgQMLcDfAOvy5gd4wN3Hyn8V5XX80i7+65Hds5r/oD+VqWDQMPvuif3DafpTGc3RICJSRTUfNIQ+CgwBHwSOAZ4A3uzuP8zbr5HDr2kdsCr8OX9fgDXAtpK2tArWLu3iuyMZ9g2nWdzVOqNj+1MZeipUI9DV0kSDBc85U5MjJ9Q9ISJSNXURNLh7FrgyvBXb75y83+8Eyt9ZX2W5IyhmEzQcM7+tHM06QkODMb+9mYMjswkawpETyjSIiFRNPdQ0yDSiNSg2zaIYsj81XrHuCYCejhYOzibTsGeI1qYGlve0l6FVIiISh4KGBDhmXhtdrU1smWHQ4O70p9LMq2DQEGQa0jM+7qm9w6xZ3EljBUZ5iIhIYQoaEsDMOH5J54yHXaYyWTJZp6e9pUwtO1JPR/Msaxo0ckJEpNoUNCRE75Iutu4dntEx0Yd3Jbsn5rfPPGgYzWR5ev+I6hlERKpMQUNCrFncyTMHU6RmsIJkNYKGnlkUQm7fN8KEa+SEiEi1KWhIiN7wA3Um2Yb+kSpkGjpaGBjNkJ2IvyqnRk6IiNQGBQ0J0Rsu4jSToCEaxVCpeRogyDS4z2xWyKjAs1eZBhGRqlLQkBCrF3cAwaJOcVWleyIMUGbSRbGlb4jl89voqMCiWiIiMjUFDQnR0dLE8vltPDWDTMNAGDRUcsjlZNAwg2LIp/YOa+SEiEgNUNCQIL1LumYUNPSnMphBd2vlvsHPD4d3xp2rwd3ZsmdI9QwiIjVAQUOCrFncyVN9Q7jHKzKMFquqxLLYkagrJO6wy+cGxhhOZzVyQkSkBihoSJDeJZ0Mjo6zdyjet/iDI5Vb4TISdU/EDRo0ckJEpHYoaEiQ3iUzG0FxMJWhp8JBQxSkxC2EnAwaVNMgIlJ1ChoSpHdxkMKPO4Ji39AYi2a4KuZcNTc20NXaFD9o2DNEV2sTS7sr204RETmSgoYEWd7TTktTQ+xiyL1DYyzuqty6E5GZLFq1pW+Y45d0YqaFqkREqk1BQ4I0NhirF3XwVN/0QYO7s28oXfFMA8Dirhb2DscLGp7qG5rsdhERkepS0JAwvYu7eGrv9N0T/akM4xPO4qoEDa30DY5Nu9/w2Di7+kc1ckJEpEYoaEiY3iWd7Ng3QiY7UXS/aIRFNbonlnS3sndo+qAhKujUyAkRkdqgoCFh1izuZHzC2XkgVXS/feGH9qLO6mQa9g2NTbtolUZOiIjUFgUNCRP1/083gmIy09BdnUzDhMOBaYoht+wZosFg1aKOCrVMRESKUdCQMFH//3TFkPuGq5tpAKata9jSN8zKhR20NjVWolkiIjINBQ0J09PRwoKO5mmHXe4dHMMMFnZWJ9MATFvXsKVPa06IiNQSBQ0J1Luka/ruieE0CztaaKzguhORKGgolmnITjhbtbqliEhNUdCQQL2LO2NlGqox3BIOjdgolmnYdTDF2PjE5CyXIiJSfQoaEmjNkk76BscYHJ16quZ9w2kWVWG4JUBXaxNtzQ1FMw2bw0yJJnYSEakdChoSqHdxNIJi6mxDMIV0dTINZsbirtaiq3Fuem4QgBOXKWgQEakVChoSaG34QbtpT+G6Bnfn2f5Rls2r3iJQS7qLzwr5xLNDLO1upaejOtkQERE5koKGBFq1sIOWpgaeeHag4ON9Q2OMjU+wYmH15j9Y0tXKnsHRKR/ftGeQE5d1V7BFIiIynboIGsyswcwuN7NtZjZqZg+Z2YWzOE+vmY2YmZvZCeVoay1oamxg7dIunniucKYhmi3yuAXtlWzWYY5d0M7OAyncj5wVcmLC2fTckIIGEZEaUxdBA/AJ4OPAl4HfAu4Dvm9mb5jhef4B6C9t02rTScu6p8w0PL1/BIDjFlQv07BqYQcj6WzBuoadB1KkMlnVM4iI1JiaDxrMbClwGXCVu3/O3e9w9/cCdwBXzeA8bwNOBz5dnpbWlpOO6ea5gTEOFpiquRYyDSvDqaF3hAFMriejIshjlGkQEaklNR80AOcBLcB1eduvA15gZmumO4GZLQCuJgg+Dpa6gbXopPAD94lnB494bOeBFIs6W+hoaap0syatDOspni4QNDwRBg1rNbGTiEhNqYeg4RRgDNict31jeL8uxjk+Azzu7t8uZcNqWRQ0PF4waBipapYBDnWNbN93ZNDw6539rFrUQXdbc6WbJSIiRVTvq2Z8C4GDfmTF3P6cx6dkZq8E3knQNRGLmV0KXAqwcuXK+C2tIcfMa2NxVwsP7Tx4xGM7D6RY97x5lW9UjrbmRo6Z11awe+LhnQf5jdVFX1YREamCimcazOy14eiF6W53luC5WoBrgC+4+6Nxj3P3a919vbuvX7JkyVybURVmxmkrenjw6YOHbZ+YcJ45kOK4hdXNNEBQ17Bj/+ETUPUNjrGrf5QXHTe/Sq0SEZGpVCPT8HPg5Bj7RV9BDwA9ZmZ52Yboq+h+pvYhYAHwJTPrCbdFQwa6zazb3Y/M3yfE6SsXcOtje+gfyTC/I0j1P3MwRTo7wYoqjpyIrFzYwX9v6jts28NhZuSFx/VUvkEiIlJUxYMGdx8BHp/BIRuBVuB4Dq9riGoZimUQ1gHHAM8UeOx+4CHgtBm0pa6ctqIHgAd3HuTsE4OMyUOTH8rV/ya/ZnEn1/9qJwOjGeaF9QsP7eynweDUY6vbfSIiIkeqh0LIm4EM8Pa87RcBj7j71iLHXgW8Ou8WDbm8CHh3aZtaW1543HzM4MEdBye3PbjjIK1NDTz/mOp/KEdBzf3bD0xu+9X2/Zy4rLuqIztERKSwmg8a3H0PwXDJy83sz83sHDP7KnAucHnuvmZ2m5ltzjn2cXe/M/fGoSzHL9x9Q4Uuoyq625p5/jHz+NnmQ10ADz59kFOPnU9LU/Vf+tNX9tDUYPxya9DDNDia4Zdb909mRUREpLZU/5Mjno8CVwIfBH4CnAW82d1/mLdfI/UxIqRiXn/KMWzYfoDnBkbJZCf49TP9k9/wq62jpYlTj53P/2wLgoa7n9xLJuu8dt2yKrdMREQKqYugwd2z7n6lu69y91Z3f6G7X19gv3PcffU05/qGu5u758/7kEhvfOExuMPNjzzLI8/0MzY+wekre6rdrElnrFnIQ0/3M5rJcutjz7Ggo5kXr1xQ7WaJiEgBdRE0yOydsLSbtUu7+O4vd/C5W56gu7WJlx+/uNrNmvSy3kWksxN86seP8V+P7OY1Jy+jscGq3SwRESlAQcNR4LLzTmLzniHu2byPv3jdiSzsbKl2kyadfeISzjtlGd+6dzvz2pr5q/NOqnaTRERkCur/Pwqcd8oxfP1dL+HuJ/u46MxV1W7OYRoajM+/+TQW//gx3vKSFSyd11btJomIyBTsyNmZJdf69et9w4ZED7IQERGZZGa/cvf1hR5T94SIiIjEoqBBREREYlHQICIiIrEoaBAREZFYFDSIiIhILAoaREREJBYFDSIiIhKLggYRERGJRUGDiIiIxKKgQURERGJR0CAiIiKxKGgQERGRWBQ0iIiISCxa5XIaZtYHbC/hKRcDe0t4vmrStdSepFwHJOdaknIdoGupReW4jlXuvqTQAwoaKszMNky15Gi90bXUnqRcByTnWpJyHaBrqUWVvg51T4iIiEgsChpEREQkFgUNlXdttRtQQrqW2pOU64DkXEtSrgN0LbWootehmgYRERGJRZkGERERiUVBQwWY2Qozu97M+s1swMxuMLOV1W5XMWZ2jpl5gdvBvP0WmNk/mdleMxs2s1vN7AVVajZmdpyZ/b2Z3WtmI2GbVxfYr83MPmtmu80sFe7/qgL7NZjZ5Wa2zcxGzewhM7uwxq6l0OvkZnZaLVyLmf2umf2HmW0P/62fMLO/M7PuvP1i/S3Ffe2qcR1mtrrI69FTC9cRPvd5Zna7mT1rZmNmttPM/t3M1uXtF+u9q5rvA3GuJe77WbWvpUBbbg7beeVs2liWvzF3162MN6AD2AQ8Avw2cAHwa2AL0Fnt9hVp9zmAAx8Azsy5rc/Zx4CfATuBtwKvB+4iGDN8XBXb/RzwY+An4TWsLrDfd4CDwHuA1wA3ACngtLz9PgmMAZcBrwauASaAN9TQtTjw9bzX6UygoxauBbgP+Hfg7cDZwIfCf/v7gIaZ/i3Ffe2qdB2rw9fjUwVej8ZauI7wud8KfBb43fBa3gFsBAYIxuhDzPeumbx2VbyWc5jm/awWrqXAde0O233lbNpYjr+xiv0DHK034INAFjghZ9saYBz482q3r0i7o/9kry2yzwXhPq/O2TYf2A98qUrtbsj5+d0U+KAFXhRuf1fOtibgCeCmnG1LCT5k/zbv+NuAh2vhWsLHDntTmeJcVbsWYEmBbe8M233uTP6W4r52VbyO1eHv757mXFW7jiJtOils01+Ev8d676rR94H8a5n2/ayWrgVYADxLEBTkBw1V/b+i7onyOx+4z903RxvcfStwD8GLX8/OB3a5+x3RBnfvB/6TKl2bu0/E2O18IAN8L+e4ceDfgPPMrDXcfB7QAlyXd/x1wAvMbM3cWzy1mNcSV9Wuxd37Cmz+n/D+2PA+7t9S3Neu5GJeR1xVu44i9oX34+F93Peumnsf4MhriatWruXTwCPu/t0Cj1X1/4qChvI7hSC9l28jsK7A9lrzHTPLmtk+M/vXvP7MYte20sy6KtPEGTsF2OruI3nbNxJ8sJ6Qs98YsLnAflBbr98fhf25I2H/7ivzHq+1azk7vH8svI/7txT3tauU/OuI/J2ZjYe1ADcV6G+uiesws0YzazGztQTdVc8C0QdV3PeumngfmOZaIsXez6AGrsXMXkGQwXr/FLtU9f+KgobyWwgcKLB9P0EKqlb1A58nSIufC3wCeC1wr5ktDfcpdm1Qu9c3XbsX5twf9DCvV2S/arsO+GOC1+dSYBFwu5mdk7NPzVyLmR0LXAHc6u4bcp4/zt9S3Neu7Ka4jjGCD6z3EtSNXAa8APi5mZ2cc3itXMcvCNr8JPBCgm6WPTltiPPeVSvvA8WuJc77GVT5WsysheDv53Pu/sQUu1X1/0rTbA6S5HP3B4AHcjbdZWZ3A78E/hT466o0TI7g7u/I+fW/zexGgm8iVwKvqE6rCgu/Bd1IkDZ+V5WbM2tTXYe77wbel7Prf5vZzQTf7j4KXFTJdsbwDmAe0EsQ4PzUzF7h7tuq2qrZmfJa6uj97K+AdoKi5ZqkTEP5HaBwdDpVFFiz3P1+gij+JeGmYtcWPV6Lpmv3/pz9eszMptmvprj7IPAjDr1OUAPXYmbtBP2uvcB57r4zr31x/pbivnZlM811HMHdnyaods9/Pap6HQDu/pi7/yLsO38N0AV8JHw47ntXTbwPTHMthfbPfz+DKl5L2FXyUeD/AK1m1pMzTDf6vXEGbSzL35iChvLbSNC3lG8d8GiF21IqUYq72LXtcPehyjVpRjYCa8ysI2/7OiDNoX7/jUArcHyB/aD2X7/croiqXouZNQPXA+sJhnj+Om+XuH9LcV+7sohxHcXkvx5Vu45C3P1g+LxRX3fc966aex8ocC1Fd8/5uZrX0gu0EXQ3Hsi5QZA5OUDQ1VXV/ysKGsrvJuBMM+uNNlgwQc9Z4WN1w8zWEwxl+mW46SbgWDM7O2efecD/orav7T+BZuD3og1m1gS8BbjF3cfCzTcTVB+/Pe/4iwgqm7dWoK0zFr4Gb+LQ6wRVvBYzayAYL34u8Nvufl+B3eL+LcV97Uou5nUUOm4lQTdR7utRteuYipktA55PMA8DxH/vqrn3gQLXUmif/PczqO61PEhQB5N/gyCQeDXBB311/69Uatzp0XoDOsMX+tcEw2HOBx4CngK6qt2+Iu3+DkGf+P8meJP8C4LJQ3YAi8N9GoCfA08Dv08wrO9OgrTXiiq2/XfD21cJvkX8Ufj72Tn7/BtB5P5uglTm9cAo8OK8c10Vbv9zgrHeXyWYEOlNtXAtBN9Avga8LWzfxeHfWhp4ZS1cS07br+TICY+Om+nfUtzXrkrX8XngC8CbCd7k3wdsJ5hg56RauI7wuX9AkAa/IGzne4HHw3aeGO4T671rJq9dFa9l2vezWriWKa4vf56Gqv5fqfg/wNF4A1YC/0EwQ9kg8H8pMElPLd2Ay4GHCaqOM+Ef6LXA8/L2Wwj8S/gHO0IwWdCLqtx2n+J2Z84+7cDVBMOyRgkqr88pcK5GgiKp7QSV2Q8Dv1sr10Lw7eKe8A0wQzA+/SbgjFq5FmBbkev4+Ez/luK+dtW4DuAPCOZuOBC+Hs8C/0pewFDN6wif+8PArwg+WEcIJvy5Jv99Ke57VzXfB+JcCzHfz6p9LVNc32FBw0zaWI6/Ma1yKSIiIrGopkFERERiUdAgIiIisShoEBERkVgUNIiIiEgsChpEREQkFgUNIiIiEouCBhEREYlFQYPIUcrMPMZtm5mtDn++pNptjpjZsWY2HE4FXInnMzN7wMz+qhLPJ1KrNLmTyFHKzM7M2/QDgmmCP56zbYxgcaLTgS3u3leZ1hVnZv8CLHX3N1XwOX+HYBa+4929Jlc4FSk3BQ0iAoCZbQN+5u4XVbstxYSLET0N/I67/6iCz9sI7AS+4O6fqdTzitQSdU+ISFGFuifM7BtmttPM1pvZz80sZWZPmNkbw8f/POzaGDCzG81sSd45m8zscjN73MzGzGyXmX3ezNpiNOkSgnUQfpJ3zjvN7Gdm9nozezBs0wNm9tLw+T5lZrvNbH/Y/s689nzCzLaY2aiZ7Q3P9YpoH3fPAt8nWPxH5KjUVO0GiEjdmgd8C/gcsAv4KPAfZvYV4ETg/cAy4IvAVwhWfoxcR7DY1qcJVuw7GfgEsBq4cJrnfT1wr7uPF3jsBOCzwCeBIeAzBAt43UTwfndJ+FyfBfYAUY3Ch4E/C6/hwfDa1hMsDJTrbuADZtbr7k9N006RxFHQICKz1Q28z93vBjCzXQQ1EW8C1oXfzDGzUwk+aBvdPWtmrwTeAlzs7t8Kz3Wrme0HrjOz09z9wUJPaGYGvJRg+elCFgEvjz7QzawBuBFY4+6vDff5iZm9Cvg9DgUNLwNucff/P+dc/1ng/A+E92cSLBEtclRR94SIzNZwFDCEHg/vb40ChpztTcDzwt9fD6SB68NugSYzawJuCR9/VZHn7CFY7neqgswn8zIAUZt+krff48BxYRACwXLWbzCzT5rZK8ysZYrzR8+7vEgbRRJLQYOIzNbB3F/cPR3+eCBvv2h7VK+wFGgBhoFMzm1P+PiiIs8ZnWNsiseneu5C25uAxvD3TwF/A5wP/Dewz8y+bmaL845LhfftRdookljqnhCRStsHjAKvnOLxXdMcC7CglA1y9wxBfcWnzewYgi6Wq4EOgq6USFTjsLeUzy9SLxQ0iEil3UxQeDjf3W+byYHunjazrUBvWVoWPMezwD+Z2RuAU/MeXhPeP1Gu5xepZQoaRKSi3P1OM/suQU3D1cAvgQmCkRNvAD7s7k8WOcXdwBmlbJOZ3UhQxHk/QVfG6QS1F9fk7fpSgq6U+0r5/CL1QkGDiFTDRcAHgD8gGOY4BmwjKFh8bppjvwe808xWu/u2ErXnboLRFO8n6JLYQTBc85N5+70JuMndR0r0vCJ1RTNCikhdCYdRbgK+7u5XVvB5lxPMRPm6mXariCSFggYRqTtm9naCQsU1lfrWb2ZfAF7k7udW4vlEapG6J0SkHv0rcCxBHcSj5X6ycD6HZ4Fry/1cIrVMmQYRERGJRZM7iYiISCwKGkRERCQWBQ0iIiISi4IGERERiUVBg4iIiMTy/wDndax0Nx5NvgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "src.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For our PMLs, we will need some damping parameter. In this case, we will use a quadratic taper over the absorbing regions on the left and right sides of the domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Damping parameterisation\n",
    "d_l = (1-0.1*x)**2  # Left side\n",
    "d_r = (1-0.1*(grid.shape[0]-1-x))**2  # Right side\n",
    "d_b = (1-0.1*(grid.shape[1]-1-y))**2  # Base edge"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now for our main domain equations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from devito import Eq, grad, div\n",
    "\n",
    "eq_v = Eq(v.forward, v - dt*grad(p)/density,\n",
    "          subdomain=grid.subdomains['main'])\n",
    "\n",
    "eq_p = Eq(p.forward, p - dt*velocity**2*density*div(v.forward),\n",
    "           subdomain=grid.subdomains['main'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will also need to set up `p_i` to calculate the integral of `p` over time for out PMLs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "eq_p_i = Eq(p_i.forward, p_i + dt*(p.forward+p)/2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And add the equations for our damped region:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Left side\n",
    "eq_v_damp_left = Eq(v.forward,\n",
    "                    (1-d_l)*v - dt*grad(p)/density,\n",
    "                    subdomain=grid.subdomains['left'])\n",
    "\n",
    "eq_p_damp_left = Eq(p.forward,\n",
    "                    (1-gamma*velocity**2*dt-d_l*dt)*p\n",
    "                    - d_l*gamma*velocity**2*p_i\n",
    "                    - dt*velocity**2*density*div(v.forward),\n",
    "                    subdomain=grid.subdomains['left'])\n",
    "\n",
    "# Right side\n",
    "eq_v_damp_right = Eq(v.forward,\n",
    "                     (1-d_r)*v - dt*grad(p)/density,\n",
    "                     subdomain=grid.subdomains['right'])\n",
    "\n",
    "eq_p_damp_right = Eq(p.forward,\n",
    "                     (1-gamma*velocity**2*dt-d_r*dt)*p\n",
    "                     - d_r*gamma*velocity**2*p_i\n",
    "                     - dt*velocity**2*density*div(v.forward),\n",
    "                     subdomain=grid.subdomains['right'])\n",
    "\n",
    "# Base edge\n",
    "eq_v_damp_base = Eq(v.forward,\n",
    "                    (1-d_b)*v - dt*grad(p)/density,\n",
    "                    subdomain=grid.subdomains['base'])\n",
    "\n",
    "eq_p_damp_base = Eq(p.forward,\n",
    "                    (1-gamma*velocity**2*dt-d_b*dt)*p\n",
    "                    - d_b*gamma*velocity**2*p_i\n",
    "                    - dt*velocity**2*density*div(v.forward),\n",
    "                    subdomain=grid.subdomains['base'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We next add our free surface boundary conditions. Note that this implementation is based off that found in `examples/seismic/acoustic/operators.py`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "from devito import sign, norm\n",
    "from devito.symbolics import retrieve_functions, INT\n",
    "\n",
    "def free_surface_top(eq, subdomain, update):\n",
    "    \"\"\"\n",
    "    Modify a stencil such that it is folded back on\n",
    "    itself rather than leaving the model domain. This is\n",
    "    achieved by replacing the symbolic indices for some\n",
    "    function of those indices. Depending on how this is\n",
    "    done, this can be used to implement a pressure or\n",
    "    velocity free-surface. This is the MPI-safe method\n",
    "    of implementing a free-surface boundary condition\n",
    "    in Devito.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    eq : Eq\n",
    "        The update stencil to modify\n",
    "    subdomain : FreeSurface\n",
    "        The subdomain in which the modification is\n",
    "        applied\n",
    "    update : str\n",
    "        The field being updated: 'pressure' or 'velocity'\n",
    "    \"\"\"\n",
    "    lhs, rhs = eq.evaluate.args\n",
    "    \n",
    "    # Get vertical subdimension and its parent\n",
    "    yfs = subdomain.dimensions[-1]\n",
    "    y = yfs.parent\n",
    "    \n",
    "    # Functions present in stencil\n",
    "    funcs = retrieve_functions(rhs)\n",
    "    mapper = {}\n",
    "    for f in funcs:\n",
    "        # Get the y index\n",
    "        yind = f.indices[-1]\n",
    "        if (yind - y).as_coeff_Mul()[0] < 0:\n",
    "            # If point position in stencil is negative\n",
    "            # Substitute the dimension for its subdomain companion\n",
    "            # Get the symbolic sign of this\n",
    "            s = sign(yind.subs({y: yfs, y.spacing: 1}))\n",
    "            if update == 'velocity':\n",
    "                # Antisymmetric mirror\n",
    "                # Substitute where index is negative for -ve where index is positive\n",
    "                mapper.update({f: s*f.subs({yind: INT(abs(yind))})})\n",
    "            elif update == 'pressure':\n",
    "                # Symmetric mirror\n",
    "                # Substitute where index is negative for +ve where index is positive\n",
    "                mapper.update({f: f.subs({yind: INT(abs(yind))})})\n",
    "    return Eq(lhs, rhs.subs(mapper), subdomain=subdomain)\n",
    "    \n",
    "fs_p = free_surface_top(eq_p, grid.subdomains['freesurface'], 'pressure')\n",
    "fs_v = free_surface_top(eq_v, grid.subdomains['freesurface'], 'velocity')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And our source terms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "src_term = src.inject(field=p.forward, expr=src)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct our operator and run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "tags": [
     "nbval-ignore-output"
    ]
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Kernel` ran in 0.05 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "PerformanceSummary([(PerfKey(name='section0', rank=None),\n",
       "                     PerfEntry(time=4.9e-05, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
       "                    (PerfKey(name='section1', rank=None),\n",
       "                     PerfEntry(time=0.0036749999999999908, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
       "                    (PerfKey(name='section2', rank=None),\n",
       "                     PerfEntry(time=0.013221000000000004, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
       "                    (PerfKey(name='section3', rank=None),\n",
       "                     PerfEntry(time=0.0037539999999999913, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
       "                    (PerfKey(name='section4', rank=None),\n",
       "                     PerfEntry(time=0.013060000000000006, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
       "                    (PerfKey(name='section5', rank=None),\n",
       "                     PerfEntry(time=0.003929000000000017, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
       "                    (PerfKey(name='section6', rank=None),\n",
       "                     PerfEntry(time=0.003147999999999996, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from devito import Operator\n",
    "\n",
    "op = Operator([eq_v, fs_v, eq_v_damp_left, eq_v_damp_base, eq_v_damp_right,\n",
    "               eq_p, fs_p, eq_p_damp_left, eq_p_damp_base, eq_p_damp_right,\n",
    "               eq_p_i]\n",
    "               + src_term)\n",
    "\n",
    "op(time=time_range.num-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is important to remember that the ordering of equations when an `Operator` is created dictates the order of loops within the generated c code. As such, the `v` equations all need to be placed before the `p` ones otherwise the operator will attempt to use the updated `v` fields before they have been updated.\n",
    "\n",
    "Now let's plot the wavefield."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "tags": [
     "nbval-skip"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfUAAAGDCAYAAAAyM4nNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAABiAElEQVR4nO2dfZgdRZX/PyeEEELAgOFVCC+iYhBFjSyurAZERUXUFRUVRRFQf4K6iqss6iLgiisiqy4r8Q0XXHVhfQF1QREERCJGXpQoKEKWd0lMEEggb5zfH31n+nTldk9P3zszPTPfz/Pc51Z3V1VX1+2+p+vUqXPM3RFCCCHE+GfKWDdACCGEEP1BQl0IIYSYIEioCyGEEBMECXUhhBBigiChLoQQQkwQJNSFEEKICcKkEupmtpOZXWBmfzWzB83sO2Y2Z6zbJYQQYnQxsx3N7PNmdo2ZrTIzN7NdapadYmYnmNkSM3vUzG40s9eU5D3azG42s9VmdouZvbOvF5IwaYS6mc0ALgP2AI4A3gw8CbjczDYby7YJIYQYdXYHXgesAK4aZtlTgJOALwAvBRYC55vZy2ImMzsaOBv4H+Ag4HzgLDN7V08tr8Ami/MZM3svcAbwFHe/tbNvV+CPwD+6+xlj2T4hhBCjh5lNcffHOumjgC8Bu7r7kiHKbQPcCZzm7v8c9v8U2Nrdn97ZngrcA/yvux8R8n0VOATY3t3X9veqJtFInawTFw4IdAB3vx24GnjlmLVKCCHEqDMg0BvwEmAacF6y/zxgr85gEeC5wNZd8p0LPB7Yr+H5K5lMQn1P4KYu+xcDc0e5LUIIIcYnewKrgVuT/Ys733NDPthQ7qT5+spkEupbkc2dpCwHthzltgghhBifbAU84BvOXS8Px+N3KnfSfH1l6khUOlEws2OAYwA222yzZ++xxx7DKh9/8nXryvNtPGV9vvHoo7EBeXqTTQplHnhoo8H0X/6S738sKJS22aZ4nsc9Fu6tWGhquA0e//hCmfUzNh9Mrw2zP7FpG+VNYUrymhjzxbQQY0WZGVG6Pz5L8VhMT9s4bDzwQLGCFeF5i5Vtt91g8s8rZxaK3H9/nt54465FmLXp6uJ5HnmErkyfPphcu9H0wqHYnPj4x2e5DkuWLGHZsmV9f7J3N/NVPdZxbzYiDn+oLHD3BT1W23omk1BfQfcRedkIns4NsABg3rx5fu21i4Z1wijI4/O9fn0x3w4zH8w3br45T8enbffdC2Uu/NkWg+lzzsn3x+f7uOOK53nZw/+db5wXpnlmz87Thx9eKPPgvAMG03fdle8P/xfMmtV9PxQvYWrNu20KTae6xGTmsZqKx7IX7HR/fL+O6ZhvznZr8o3vfa9YwQUXdK/g+OMHk2csen6hyBe+kKejIP/wh/P0IU+7rXiem4J2NzbuaU8bTP75cU8uFFm5Mk/Hx39m8R2jQHwuB/p6n33mlRfogUeAXs3DPwaPuvtINHAFMMvMLBmtD4y8l4d8kMmdeyvy9ZXJJNQXk89xROYCvxvpk8e37k03TQ7GBz5KxV12GUz+ZkkuxKH4HEehOn9+nn7Zvsk985Gf5emHH87T88J9v/fepeeJA5H4h7NZWBCY/ilMmzp8AV33z1mISN2XwWmUSPWp0wqb8bFcsiRP33dfnr51Vl7mgAMPLNYXH56f/CRPf+tbg8n3n7Z3ocjatflzvnhxvv/WMHt729N2K5TZ7Wl0zxguYPMdi01bnQz26zDaz2WL/wUWA5sAT6Q4rz4wR/67kA8yuXNvRb6+0uJ+6zsXAvua2eAT0XE08LzOMSGEEC3AyIRTL58R5GJgLfCmZP/hwE2dVVUA1wDLSvItJ1t51Xcm00j9S8CxwPfN7COAkzkQuJPMOYAQQohJhJkd2kk+u/P9UjNbCix19ys6edYBX3f3twO4+/1mdgZwgpk9BFwHvB44gGzpNJ18a83so2TOZu4GLu3kORI4zt3D3E3/mDRC3d1XmtkBwGfJ1gka8FPgfe7+cGXhPhC16jOmJr/lkqDPC3Nia6bnqrhbk8UTUTUYps449NCQKc7pQVEdGBsUDADveLhokBnPG6fropo92ghsMG9eNoFZd4JdiApKVe7pfRe3S+7JqTPL1e9R5R4fo8isWcVn51lxKiuq3y+9NE+ffnqhzAdPOnkwHW1lli7N09HsBmC3g3bJN+I8wbJlg8kZuwS7HWDmzPy/pa2P4iipkc9Pts/qfF8BzO+kN+p8IicCDwPvBbYDbgFe5+4/iJnc/Ytm5sAHgA8CdwDHuvtZjBAt/TlHBne/A+jqn1cIIUR7GA2h7u5DWu53y+Pu64FTO5+hyp/NKGqDJ5VQF0II0X4G5tTF8FG/CSGEEBMEjdRHicK8VeqkIk5ch3mwaaHQvHlPLxSJa0vjEvYdbr4s3/j5z4vniZOEcb4vpNO5+3vDQozNcz80pet20+nKwiyl5tfFSFJ1I5bdsMG2JC0SV32Gx7LgryFWuyhxY/Gsg/fNN3bdNU/HtWrnF6d0p4T1qUe+M4/Q+bslMwbTiX+o4iR/TEfDl8TPxdRZxSWybUQjzmbo31QIIUSrkPq9ORLqQgghWoeEejMk1MeCVN0cdXvR1eQP8tURcxJPb3OCtzkWPZCno/ot6gmh4KEuepFbtWPuQnLp9V1bvAFRVVnmThNg2kzdYmKUqFK/l039BPX7isRZdHwsoyq+7N6Py94Alk/fYTC91b5BFX99eMjiEjSAL34xT99ww2Bybnx2Y2OguMYtHjv44Dyd+G+O19DW2S+FimiGXoaEEEKICUJL39GEEEJMVowNvb2IekiojxIFtfTMouepaU94Qr4RzM/v+OMfB9OrLrqoUGZGSM8K6S2imi1R2RdczwV1YNQAphHkosV7GoFtgKjxSzWD06fnyqBozV8Vi7ZbNCghBih4kasbci3ee+FGXrUuX58RoxFDcZFKDDscjcpjROQkOnLhuSqo36OZfKJ+vy888w/G5z/kSa84xlDaI0Z3OvbYvK7pxTjMj4bntOy5Hmv05DdDQl0IIUSrkPV7cyTUhRBCtA4J9WZIqI8SjzySp1OD1933f/lgesbBPx1MP/rZzw6mf5HUF7XcIbQ5ewc9/5NTJzcxCHpIrwvtido7gI3CxFaMCR+JUwvpKWP5LbcsUcU3oEq7KiYpFVM6Uce8Zmo+eXVfWCCSWq/Hezne+/ExiqrruB+Se3J2OBgt2ZP5riUhfV1I3x/SqduY/eLGS186mFzzqtcNpm9YWCyzYxJfXUwc9FcohBCidWik3gwJdSGEEK1Cc+rNkVAXQgjROiTUmyGhPgwGltLUXWZVNs+bOnq79NI8fdBRZwym54b1Yau+9KVCmRg74oGQjtOCu0VPU8DU6G0uHHvavH0G03G5DhTnFet41Urn1Mv6YPPN8z6su6Smasq06thQbRFjR88xfmp4ioNkHj08JHdVzKnH+zo+F3EF6vbb5+k4VQ4wZ7s1+cbPQ6SkcNKH41o5is9yYRVsSO9dPA377L//YPrPn/2vwfRXTs/zpHPosa1lfV1YOoiWl44X9DcnhBCiVRhyE9sUCXUhhBCtQx7lmiGh3iOpiioS1VWbbprvT9VdMez56UFldthhCwbTZyx+X6HM3ocdNpi+47e/HUw/GPLEZTAAO5x7br4R9ORT9ssXxewWAr0AhTjMf34kX0zz5z/nWWLgizSgSzxW1/lXHXV8HXU7FPu6bpk6dYne+7Nu3ZX9Hg+G9BqmFbKVqdzvvTdPb+gNMU/H5WrRMeOcmcvzjRCABYALwnZ8yMN8W3xeAaKW/MkhvdtOO+Ub55xTKHPGDQcMpr/w3Hx/VLG/4x3F88TphGlTy//Dyqj63+sHMpRrjv6mhBBCtA4J9Wao34QQQogJgkbqDaireor5Zs7M359SK9nZs/P0vffmVrKf/ewNIb1Xoczb3vabwfRXfx6UeO9852By3Te/WShzf1C5Twuq+Fkxhvt+Bf9UMH/+YHLbcGzbEKDiD7fm1xbV7VBUoz70EF1JvdiNRaznOucZSXXzRKZJv1X9HgUr7Km5mj1O/aT3YVSzR1X8ypV5Op32ic9pnJWasejKfOPii/N0XMYCrPvVrwbTy+nODun20UfnG1/4wmDyLUfl13nuC/8vKXVhSOfTZfPmzR1MP/WpxRJbzWqfyn3D84kmSKgLIYRoFZpTb46EuhBCiNYhod4MCfVRYsq63BHF7rsXLXNf9ao8fcEFufpsxYofhFzfLpT52td2DOm3DqY/+tHc+cTJV7wzFmGbo44aTN8VYjXfE/TiM//3fwtldgzbU0KwCF6eB6F58oteNJje7mnRZrdoaRzV71FVmlodl6nfN6q5xiUG32irKr/t9HOqoUlddX+3mC86PkrV7zFWevT3Eqe+wkIPAObODutHvnZ+nj7vvMHkgwvzSCmJT6lCDPRwGnZ50pO61gXwnvNyR1Cf3+TqcCTmC3MJADxvMGV2yGA6LJDh6U9LVOfxoYvzDi25eTVSb476TQghhJggtOO1TAghhAjIo1wzJNSFEEK0DnmUa4a5+1i3oStmdijwBmAesA1wB/Ad4F/c/aFOnl2A20uq2NLdHwj1TQdOAQ4HZgE3AB9y9yu7FU6ZN2+eL7r22gZX0qEq6klY0/Vf38sDT7zpTdE92yeTCq8rOdGLB1Nbb31c4UiY/mO3M98zmH74858fTN9KkTgdulVI7xYnrqNRwMEHFysIS+JWzZ4zmK7jaQ76O8WX1hXn6OPl9Hr+lkxL9kSv8+CRJGZJgfXr83SZzUTZbwPlthkrVpSfJ3pTi97hpvzvD4uFvvGNwWRcHhrDJEWvjakjxGhdMjssVbvtw7mnyL33LpZ56KFPh62yv6anJNv/OJg6++xtBtPHHB5m9VMjg7iONI3iNAzm7bMPixYt6vug+slm/vmhs1VyEPza3edV5TGznYDPAi8iUw5cCrzP3e8YotxJwD+XHF7t7tND3iXAzl3yvdrdv1d1nia0+e/neDJB/k9kNijPBE4C9jezv3X3aPnxSYqLNQHSVdFfAV4OfBC4DXg3cImZPdfdb+h764UQQrQWM5sBXAasBo4AHDgVuNzMnu7uKyuKfxm4ONm3WWdfKosALiGTX5FbGjR7SNos1F/h7kvD9hVmthz4OjCf7McY4DZ3X0gJZvYM4I3Ake7+tc6+K4DFwMnAIWVlhRBCjD6jYMV9NLAb8BR3vxXAzH4D/BF4B3BGWUF3v4tkwYOZvZlMpn69S5FlVTKqn7RWqCcCfYABF01P6HKsikOAtYR1Ye6+zsy+BXzYzDZx99XNWtqAVE8ZAkG8ce9Zg+knLsw9Qh199McLRX7726gq/H5I/3gwtXTpbYUyT3zi6wfTb3vb5wbTX737w4PpvQ8/vFDmnssvH0zHO3hZ0Kluc36+3GeXiy4qlOfVrx5MzjjwwMH0nOiiK3Gx99jMPHBMVKk+8kiejupU2DCQzAB1l0aVHasqU2eJXZXquAn9DFCTUqUmHyDt90jdgD2Rsv6tWmVVEsOlUCaNHz5jXfC6uGhRnj4+LBv92tcKZW4L02RRFxtVhIUALKku/dv5MtRjTs+V8V96YvxvLy5pg+gtbsuQzpeQbr/964l8+ct5+mV7hGd+UXhiY0SaccAoLWk7BFg4INAB3P12M7saeCUVQr2EI4A/k43Kx4zxtqTtBZ3v3yf7P2lm68zsr2Z2oZntlRzfE7jd3Vcl+xcD04i+FYUQQow5U3r81GBP4KYu+xcDc7vsL6UzN78/8A137/Ya+wozW2Vmq81soZm9ajj1D4fWjtRTzOwJZKryS9194FV7NXA22fB0KbAH2Rz8L8xsH3cfEP5bAYnpDJC7ZN6qyzHM7BjgGIA5c+Z0yyKEEGIE6MOIc7aZBbUMC9x9QdiukgtbdtlfxeFkTe6mer+ITMt8O7AtcCzwXTN7s7unqpqeaa31e8TMZgI/I4t/sE9nPqMs705kb1oXuvvhnX0/BrZw932TvAcCPwGe7+5XVbWhZ+v3qINMLVFvDva0twb786hPTKzKz/pW/h7yvvfl+9eujfdIaj0bddTxRTRXuX/840W95ccODlb2IVjMPSFYxX2UEz1pFV6L9t8/T+9b+Flgzz3z9K675umoU03UiavW5V76oue61WFSJVUDl6mFe7VeH+vy0LtqfqSC16TXtskmeXrzzfP0jOlByZ2uFonbcX4mPlfxOYKiyj0EYbn/zjsH0+mfSnxa4n385BjbPOi+z7jpxUQ+8IHo+S0+lzeEdGoz//yQzp/Lz3wmn+t5/7FrKPCDMIUQ+yY+L3vsUSwTn58ebriRsn5/illB+jZh/hDW72a2BjjD3T+c7D8V+LC71+4YM/s98Ii7P6tG3o2AhcB27r7TUPmHS+vV72a2Kdmbzm7AS6oEOoC73wn8HHhO2L2C7m9eA5KxLIiSEEKIUWZgTn2E1e9VcqHbCL57W832IdMSdxulb4C7rwfOB3Y0s+3rnqcurRbqZrYxcAHZWvWXuftvh1E8qiAWA7t2ljBE5gJr2HB5thBCiDHEevzUYDHZvHrKXOB3w2jqEWSG2P81VMYu9F1V3lr1u5lNAb4FvAI42N1/WrPcHDLjh++5+1s6+55J5q3lre7+9c6+qcBvgVvd/RVD1duz+j2SRjBZsiRPRzXhTd1sODoES9vHDn/LYPrDQZH06U+HKBZAMSjMr8sqTrbfOJj6xCceP5j+p1eFe/744weTDyYBYaJqPioNo7uL1KBhi+gYY69g8xiti5+SOOCIAaOjNX1QM66ZWnyniz9DtJ6PqmfFUK9P1OJGtfqmm+bpLWYmgUWiyjxG/4npVJUej8VnJ+YLAYsAlocf+IGwP7ZmVvEszI7322c+M5j86pIDBtNvf3tU0qcDtbL/i2eH9KGFI//wD7nzmNNOy/dPOz93hMPixSX1UoxKE9ufrDApOJ/pgZFSvz/VzL/aYx1/O7T6/X3A6cCT3f22zr5dyJa0fdjdP1NWNtQxjSzCzs/d/ZV12tWRPb8EZrt7N6c0PdFmQ7l/B14LfAJYaWZx4vUud7/LzD5Dpm24hsxQ7inACWTP6icGMrv79Wb2beDMzuj/duBdwK7Am0bjYoQQQrSKL5EZrX3fzD5CNmo+BbiTzAAbADPbGfgTcLK7n5zUcTDZuKSr6t3M3kC2PO5HnXq3JXN89iwyj6l9p81CfSDO54mdT+TjZN55FpMJ57eSDf7+QuaU5uPunnrreRuZoD+V7KX8RuAgdy/ztyqEEGKMGOm5YXdfaWYHkLmJPZdMa/9TMjexUZ1qZK7ouzXpCDKbrB90OQbZAHIb4NNkwn8lsIhM9ozIevbWCnV336VGnq8CtbQ07v4I8P7ORwghRIsZDYOvjo/31wyRZwkl0/RDqdw7XuQOqMrTb1o7p942+jqnnk7Uli3FifPrcdkbFCeE4/KUF75wMHnbU19OJM63n39+nKeMHulS+5DoNi369MnnAo87LjfgPPPMYukpJ30s3zjnnMHksrCU6IHkjHHufVpIbxHS6Tz81LhkJ87Dx4gdT3xisVBcLhf7sCTYxZpCa8rn4aNnttQDW5vm6OsGuCnz2jaNZGlV2fKy+4JlRTo/Hu/raEMS0uvCvQIQh1AljgQ3GK3Ee2fak56Ubxwa5rSDbQjAv3wxv8tOPDG6Ab8gpKOiL7g8BIqBV/L//te+Np/3PvXUYokn/zF4ivxpMCOK/ZkGYIn3fpxHj/Prs+OiPPoWdWik5tTnmvm5PdYxr0ZAl4lIa0fqQgghJi+tXprVYtRvQgghxARB6vea9FX9nhL1uFHNFpfr/DZZov+nP+XpqN6MRPUbFLzS/dL3GUyfdFKe5eKLU/vCH4d02dLNeJ79C0ee+czcwdJHPpLv//s9Ql1f/GKxunNzxdt9wUNW9MOXKH4LKqey5XKz0ugqUW0Z+you/4leubbdtlh+661D5bNCA2Z2T0N5pJKydLftblS5yytLp1Fwovo8pqOXsqUhztKf/1wsH+/XmA4q9seSpWYxHnk4CzFIQzryiKr0HUJ6+nOCv6moVgd461sHk9/5eb5sLE4XXXXV7cmZLg/p+FzEOZXdQvqFRF70olz9Hp+xv50a/kd+kNhXpdMTA8TpodQ7XNyO925Uucf7ro+MlPp9TzNvsug7srfU70IIIUQ7kBq5GRLqQgghWkffh/+TBKnfazKi6vdIVIlGtWeqYi/zpBX3p4FjogouemR7Re5Q7xcPP71Q5PTT8/R3vxvd7kcHf1EtH+JWA0XXyjGIzN8Npvbfv+hUKRohv2yXUPcXvpCnv/nNQpn7S9T0ZdbRUHyjnVGSnlmSBpgS1e9RPRrVnqnnrjI1ffxtqlSldQOqx/sopsvU6ul2ifX6Y0H9nvhFLPzy8VhUpafTJpHY7zFczzY7JTEvjjoqT4eb5axz8hrSGZ3f/jaqtaM3xahyT+/dqOiPU0x50JVXvzpX5SfG8/ztzN/kGyGITMHiP50CifdHVKWXTQ9B+cqNEVK5R0ZK/f40Mz+/xzrmTlL1uzQcQgghxARB6nchhBCtQyPOZkj9XpNRU79HylSoUO7cIzrqiBbyUAyEEdWwUU2XOmh50YsGk3+YlVvMR014TG8Ydyc65/i/kI6K2DT6YVTT56rO5zwndxYTY8gDvPHgoDoN1vP8MDjz+PnPC2WWh8DrD4T9ZeriJBRJ4U8npuObctFdTfFYWTr9M6vz55a2LW6vq5GG4rXGY4+VpNN2lU1nRCX2NiT8XT4NU7BYDyr2z325GIgnBjq5996rw5GFIX1HcqLoGCbeb1GtHgOtwNZb5ys3osY/pndbFv4Tfvaz4inL1OxRRf6EJxTLRNV6HedIUH/qZgQYKfX7Xmb+nR7rePIkVb9rpC6EEKJ1aKTeDPWbEEIIMUHQSF0IIUTr0IizGZpTr8mYzKlHqjyG1QmkAcU59bj0rWyuHYpzdNG72t/8zWDysZfmgWPOO69Y/MtfztNXXRWXxMW50BuLhbg3pKP3rseHdOJVi3z+c+ONc89icZr2lUk8pbCSjxmLrowNzdO//GWeTrz6PRb6MPZaVcCRsjntsnnrbtvdqJqHL5vvT9/o42xs2VK+KXHONwbOAZgXpi/33TdPH3jgYPJHFxdbekGIjRLvnbVrw29A+G2A8iWUMfhQ4k2R54V07vntta/NZ/yD0zkAXjY/WFfEJWk33JCn47OTPqNxyWKcEy/zWJhuly1/bOJxcIQYqTn1p5v5RT3WsYvm1IUQQoh2oJF6MyTUhRBCtA55lGuG1O81GXP1exVlQTpSVXqZh7qoQkxiV3P33d3LR6IHteipDgpLlm6bmXuri6rWEGYdgNtvj+rVX4V09AqWBBMpqOk3Demosk9UnYVgHLm6dsstcw93++2X54gaZSiGro7pHWYGlXAaoCP2dfwNYt+mv1u6nLEb6VKmqK6NatyoBk5Vv8Fr2T0P52rpGOZ80aI8vTCuIEu2ly6NoVpiEJc0aMptIR2nXaJfwHQCIv6m8X7Llz/uumsI7kJxCWRckjbj5yFg0a/ivUZxSWj8TaK6O977acCf6Akv9nWZ90EoX542hir2KkZS/f6jHuvYSep3IYQQYuwxihYSoj4S6kIIIVqH5tSbIfV7TVqtfi8jtcatE9gjDQIT1fFRFR9Vx7FMqiqOasOognza0/L0/PmFItctmzOY/t738v3nhwgPN98c1btQtIiOKt1ocb8iKbOS7kQ/cNEfWur5LqqBZ9fYn9aR173xxrmqtW4I9rqh0deujQejtXjaH/G3/0uN/WkAlDhtsDak45gr9Xi2WUjHvto+pItTOptvnnscPOywfP/b356n/2bLPxRPE+OWR+v1dIVIpGzaosx6PeZJt8sCraTTJi1Vs5cxUur3vc38Jz3WsY3U70IIIUQ70Ei9Geo3IYQQYoKgkfpEJlXllcXvjvtTa9yoQowBJqLa8t5gtRxV9Gm+shjwl15aKPKsoN58VlDTn/zWXJP2582i5TpcfXUeKuRnP5s/mI4xXK6/PqqEAULAjYJVdpkVdjI1UQhQs74kXUWull67Nk+vSLXitUyGqs7ZW9uK0xFRXZxOR4T7o7DSYLeSPLDXXvl54gqC6Mcm+DnKtp8ZQs/EICo/DD92NNmHDaeVBoj3e8VqgMKxmI7l03mTMocxNVXsj5WMuabUckc0vjE04myKhLoQQojWIaHeDAl1IYQQrUIj9eZIqE9WogqwTC2fHitTVZY5tUm3oyX9X4IVderUJqrmo8o+qOm3TaYJ/j5ME/z9/NC2d+Y+4h/bI8ZphxtuyH2XL1qUp6+/Ps+zeHGejuGxAZYujdbzdS3E43aM8R0t1FMVedwui2iequjLLM6jY55o2Z9uRwv+3KnK9tvndcUFDFDujCem5+6RqI6XhJUK8V6J6UVxOgT44Z+754sq9nTlR5mVerRkj+m0TLzfyizZK9TqZar0ukwGlbvoDxLqQgghWodG6s2QUBdCCNE6JNSb0dp+M7P5ZuZdPg8k+bY0sy+b2TIzW2lml5rZXl3qm25mnzaze83sETO7xsyen+YTQggxtgzMqffyqXUes53M7AIz+6uZPWhm3zGzOUOXhBL55Ga2d5JvipmdYGZLzOxRM7vRzF5Ts4nDZjyM1N9DMarH4GSZmRlwEbALcByZi6wTgMvNbG93j+7EvgK8HPggmcuxdwOXmNlz3f2GkbyAcUVVrOayZXBxjjFdEhfnKeMSoTj/mS43itvRNVpZ3HgoBkqJkUVC+6ckS46eFdr6rBiM46UhEMfRJUuZgDWz8mV0d92Ve0a7776dQ7rYzNjsuHRtZZieT73DNYnnErc3C07btgyr0KpWL5bFfZm27J58I9o7ANwelgXeFZY23lDhcTBSN5hJ3C5rdLq8rE5Qm7RDasyd9zpXLsoZ6Z41sxnAZcBq4AjAgVPJ5MfT3b3M5WTkHODsZF/izpBTgOOBE4FfA4cB55vZwe49x63ZgPEg1H/v7gtLjh0CPA84wN0vBzCza8gWHf8j2QsBZvYM4I3Ake7+tc6+K4DFwMmdeoQQQkwejiZzovAUd78VwMx+QxZW8B3AGTXquLtCPmFm25AJ9NPc/fTO7svNbHfgNKDvQn28v2YeAtwzINAB3P2vZKP3Vyb51gLfDvnWAd8CXmJmm4xOc4UQQgzFKKnfDwEWDgh0AHe/HbiaovzohZeQeW86L9l/HrCXme26YZHeGA8j9W+Y2WzgAeAS4MPufkfn2J7ATV3KLAbeYmYz3f3hTr7b3X1Vl3zTyIJpL2aS0kiFODX3MjZlZomKHspV82VqdSiPLV4WhCbdjjruqjjl8VhcRlfm/StRA08L17pbUPfuVjZNAbBJeH/cOiwv2zHs33jjYpk6HsjSJVxrg/e81avz9CNhGd2tiXbxppKAP2WBgNL+jG2I6dgf8X6Aoio8ToFsvXWejvMHddl00+J2PG9ZOv2tpGYfU0ahx/cEvt9l/2LgtTXreJeZfZBs3elC4J/d/arkHKuBW5NyA/JmLkV3lj3TZqH+V+AzwBVki3ufCfwTcI2ZPdPd7we2ApZ0Kbu8870lWeiordgwJFXMt1W3BpjZMcAxAHPm1LKdEEII0Qcyk6kecJ9tZovCngXuviBsV8mF1AdyN84DfgDcA+xMZq91mZm9yN1/Fs7xgG8YDrVS9vRCa4W6u18PBDcgXGFmVwLXks2Vf2QU2rAAWABZ6NWRPp8QQgjArPcwtGvXLhvJ0Kvu/uaweZWZfZ9Mc3wqsN9InXcoWivUu+Hu15nZH4DndHatoPsb1Vbh+MD3zhX5lnc5JmpSUE1OnVY4NqVMfV2lko1q+jKVbqr6LVPNl6XTMmXq5piuG58+zVdGWb5+xNTute6yfHF/ai1e5qUw5osqdoCdwkqDMuv1tC1lgeSr+r1sSqTiPFK5T3iq5Ee3EXwl7v6Qmf0QeHtyjllmZslofcRkz3i9awc6ZzHZnEXKXOCOznz6QL5dO0sY0nxr2HC+QwghxFgydWpvn6Gpkh+/66HlUXgvBjYBntjlHPR4nq6MK6FuZvOAp5Cp4AEuBJ5gZi8IebYAXtE5NsBFwMYE4wczmwq8HvixuwdLIiGEEGPKgPp9ZIX6hcC+ZjYYG9jMdiFbJn1hWaHyJtsWwMHk8gngYrKVV29Ksh8O3NSxtu8rrVW/m9k3yKwCryOzfH8mmWOZu4HPdbJdCFwDnNexQBxwPmPAvw7U5e7Xm9m3gTPNbONOve8iC+6cdvakIwaL6LfKMdZXqopP1aapBX030jJlKtkqa+0yNXuddN18Ve2MFurr13fP05TYvxuF4C7Rsr6uo6E6UyjpsTL1e4XHmzUhbnvsgrSZ09aFRSx1pz1qOLZpg7pdgVs69GNOfWi+BBwLfN/MPkI2wj4FuJPgUMbMdgb+BJzs7id39h1PNsC8nNxQ7nhgO4JMcff7zewM4AQze4hMnr0eOIAR8o/SWqFOZnDwBjJPcTOA+4DvkC0ZWAbg7o+Z2cHA6cBZZOGorgH2d/c7k/reBnyCzIhhFnAjcJC7XzfylyKEEKJNuPtKMzsA+CxwLtlg8KfA+8LULZ39G1HUbN8CvLrzeRzZCq2rgbe7exypQ+ZJ7mHgvWRC/xbgde7+g75fFC0W6u7+SeCTNfItB47sfKryPQK8v/MRQgjRVkZnpE7H50mlH3Z3X0Im2OO+i8imdeucYz3ZYPLUZq0cHq0V6kIIISYpoyTUJyLqNTFqFObX49xhg4f3sXTpXKyvbMlS1fx2WbosP5TPnZfN7w+n7rrHBqjqwzrL09LtCk96tcqXLFlcM7W4AKVsJWE0A9h88+JpptWxzah7bS2j9BmZjLT4d2oz6jUhhBDtQiP1xoy9uacQQggh+oJehUSBKpXfWC/5qTp/2bEpUQ3cj6ViZdRR39dtQ7+XtPWSJ6WqbSXq97hUrcoRYCSq3zdoZllfV6nYS651rO/plEmvch9AI/XGqNeEEEK0Cwn1xqjXhBBCtAsJ9cao18SYMCZWvlWBQarydahU1cb48hNNhVp3lUCJ+ntdsGqPTvRSYvEqx3c82uP0RIuFhazfAy3+ndpMuyaUhBBCCNEYvQoJIYRoF1K/N0a9JkSfaWJRPZKq1hFT6VZYmMdzRi19jFuTEi3e46KFRm2uaf0uWoqEemPUa0IIIdqFhHpjNKcuhBBCTBD0KiSEEKJdaKTeGPWaqE2c22ybJ67xzkTrz7oO9mrFiqkKxFNWqEIgTLS+nrBIqDdCvSaEEKJdaKTeGPWaEEKIdiGh3hj1mmhEP1XxafnJ4kmr1ddZV38eKcQ2H/4pS73I9epBboMK20N6D2hqQPRKO+90IYQQkxeN1BujXhNCCNEuJNQbo14TrWPMVZBlMbr7zEQO3hG7sCqIS2Xc9G6VpZRZvI8TgTBZp56GREK9Meo1IYQQ7UNCvRGyyhBCCCEmCHoVEj3TZqc0E1nFPeYkI6myIC5Nqpssv1WV9ftk6YOuSP3eGPWaEEKIdiGh3hj1mhBCiHYhod6YdulKhRBCCNEYvQqJvjKRPWSN5LVN1rnUUi9yvU7KjxO0pK0EjdQb09p/XDP7mZl5yefiTp5dKvLMSuqbbmafNrN7zewRM7vGzJ4/JhcnhBCimqlTe/vUwMx2MrMLzOyvZvagmX3HzObUKDfPzBaY2c1mtsrM7jCzb5jZrl3yLimRUa8afqcMTZtfhf4fsEWy77nAGcCFyf5Pdtn3ULL9FeDlwAeB24B3A5eY2XPd/YZ+NFgIIUQfGIWRupnNAC4DVgNHAA6cClxuZk9395UVxQ8D9gQ+BywGngB8FFhkZnu7+51J/kuAk5J9t/R8EV1orVB399+l+8zsaGAN8K3k0G3uvrCsLjN7BvBG4Eh3/1pn3xVkP8bJwCH9arcQk5UmMWCGXXFKrYDs7Z0Gkrq9hNFRvx8N7AY8xd1vzU5rvwH+CLyDbABZxqfcfWncYWZXA7d36v1Ykn9ZlYzqJ+2807vQeat6LXCRuy8fZvFDgLXAtwd2uPs6speDl5jZJn1rqBBCiPHAIcDCAYEO4O63A1cDr6wqmAr0zr7/A5aSjdrHjHEj1IFXA5sDX+9y7JNmtq4zL3Khme2VHN8TuN3dVyX7FwPTgN3731whhBCNGBipj+yc+p7ATV32LwbmDr/J9lRgG+D3XQ6/ojP3vtrMFo7UfDq0WP3ehbcA9wP/G/atBs4Gfkz2hrQH8E/AL8xsH3cf6NytgBVd6lwejm+AmR0DHAMwZ86QthOiC2XqxbaqQyc1w9WZp3+cJcXXr8/TMYBLWkXfY6iPA2T9XkJ/1O+zzWxR2F7g7gvCdpVc2HI4JzKzqcAXyeTQV5LDFwG/IlPNbwscC3zXzN7s7ucN5zx1GBdC3cx2AA4E/q2jNgfA3e8F3hmyXtWxjF8MnAgc3st5OzfAAoB58+Z5L3UJIYQYBr0L9WXuPq8fTanBF4C/BV7u7oUXBXc/Lm6b2XeBhWQG3n0X6uNluHQ4WVu7qd4LdKwOfw48J+xeQfc3r4ER+nDn6IUQQowUo6N+r5IL3UbwJU2108g0uke6+4+Hyu/u64HzgR3NbPu656nLuBipky03uNHdbxxGmTiyXgy82sxmJPPqc8ms6W9FCJHTwHy9rEiqch8xxomzkjYHQJpkLCabV0+ZC2yw+qobZnYi8CHgOHc/t0Eb+q4Bbv0dZWbzyDp5yFF6J/8cYD/g2rD7ImBjMuv5gXxTgdcDP3b31X1rsBBCiN4YnZH6hcC+ZrZbflrbBXgeG/o96dJEew/ZuvYT3f0L9S9tUPbc4e731S1Xl/HwavsWMhOcb6QHzOwzZC8m15AZKDwFOAF4DPjEQD53v97Mvg2caWYbkxksvAvYFXjTSF+AEEKIYTA669S/RGa09n0z+wjZqPkU4E4yA+xOU2xn4E/Aye5+cmffYcCZwMXAZWa2b6j3wQE/K2b2BrLlcT/q1LstmeOzZwFvGImLarVQ7wjgNwAXu/v9XbIsJhPObwVmAn8h8xD0cXdPvfW8jUzQnwrMAm4EDnL360ak8aKScaOCjDrlcaLeHVFCH1T9bmvXDru6Zpbf4/A3afX93hZGQai7+0ozOwD4LHAuYMBPgfe5+8OxNcBGFDXbB3X2H9T5RK4A5nfSt5Mtc/s02Vz9SmARmey5pJ/XM0Crnwh3XwtsXXH8q8BXa9b1CPD+zkcIIcQkx93vAF4zRJ4lZAI87nsr2WByqPoXAgc0bmADWi3UhRBCTFLGoRamDajXhBBCtAuFXm2Mek2MORM5Bvtkos4quBg/vXFl4/zPftzYk4wlEuqNUa8JIYRoFxLqjdFrohBCCDFB0KuQmBBMpJVnkynIR+G36msQ9vYilXsNNFJvjHpNCCFEu5BQb4x6TQghRPuQUG+Eek1MOMa7Kn68qNtbpy0fJz+2rN/FSDI+ngIhhBCTB6nfG6NeE0II0S4k1BujXhOiZTSxfh8vFvOj9T/dZrV2m9vWGiTUG6NeE0II0S4k1BujV0YhhBBiglD7VcjMppEFdt8B2BRYBtzSCUsnhBBC9A+N1BtR2WtmthHwauAo4AXANIpxZd3M7ga+CXzJ3W8dqYaKycNkX/LTZD581ObQ4x9t25a0jRMm+/1dC6nfG1Paa2Z2KPBJYCfgEuAjwPXAUuARYCtgV+BvyAT/+83sHOAj7v7nkW22EEKICYuEemOqeu1zwL8C57j7AyV5rgW+TSbQ/wb4EHAMcEo/GymEEGISIaHemKpe283dH61bkbv/Evh7M5vee7OEGAUKquTh65JHSo06Xpe0NfEwVwzo0uCk4/CPXyp3MZKUPhHDEej9KCeEEEIAGqn3wLB6zcwM2B7YYDTu7rf1q1FCCCEmORLqjajVa2b2eODfyQziysps1K9GCSGGR1s9yIkNkfV7DTRSb0zdXvsKsD/wBeBmYM2ItUgIIYQQjagr1PcH3uvu54xgW4QQQgiN1Hugbq8tB7T2XIiW0gbr9xFjgv25x99qQv1O/URCvTF1e+3zwDvN7GJ395FskBBCiEmOhHpjavWau59hZjsAvzOzS4EVG2bxf+5764QQQkxOJNQbUcv00sxeBrwbeErn+yNdPrUwsx3N7PNmdo2ZrTIzN7NduuSbbmafNrN7zeyRTv7nd8k3xcxOMLMlZvaomd1oZq8pOffRZnazma02s1vM7J112y2EEGJiYWY7mdkFZvZXM3vQzL5jZnNqlu27jOoHdddTnAH8CngGsIm7T0k+w1nOtjvwOrLR/lUV+b4CHA18DDgYuBe4xMz2TvKdApxEZpn/UmAhcH7nRWQQMzsaOBv4H+Ag4HzgLDN71zDaLsYxjzGl8JlITEmublyybl3+mcDE32ki35M9MaB+7+Uz5ClsBnAZsAdwBPBm4EnA5Wa2WY1W9lVG9Yu6+o05wHvc/bd9OOeV7r4tgJkdBbw4zWBmzwDeCBzp7l/r7LsCWAycDBzS2bcNcDxwmruf3il+uZntDpwG/KiTbyrwCeBcdz8x5NsBOMXMvuzua/twbUIIIXpldObUjwZ2A54yEGHUzH4D/BF4B9lgtqR5/ZVR/aTuq+H1ZHHUe8bd6wwjDgHWkgWLGSi3DvgW8BIz26Sz+yVk4WDPS8qfB+xlZrt2tp8LbN0l37nA44H9hnMNQgghRpBRGKmTyZmFMWS4u98OXA28skbZfsqovlFXqL8HON7MntfvBpSwJ3C7u69K9i8m66DdQ77VQBrHfXHne27IB3DTEPmEGJdIjTt+iL/ThJg2GQlGR6jvyYYyATK5MJRM6LeM6ht19RvfA7YArjSzlcADyXF395372K6t2NDCHrL18gPHB74f6LLMrls+utSZ5hNCCDE5qJIzW/ZQduD4wHcdGdU36gr1nwKTbn26mR1DFh+eOXNqGUQKIYToA33QOM02s0Vhe4G7L+i10rZTd536W0e4HSkrgG4j/4G3muUh3ywzs+RNqFs+yN6+7q3IV6BzAywAmDdv3qR7qRHjB6luxw9VAV30O2a492URxDJ3n1dxfAXdR+Rlo/C0bD9lVN9o6+TbYmDXzpKDyFyyYDK3hnybAE/skg/gdyEf5HPrZfmEEEKMMQNCvZdPDRazoUyATC4MJRP6LaP6RqlQN7O/H25lZra9me3bW5MAuAjYGHhtqHsq8Hrgx+6+urP7YjILxDcl5Q8HbupYMgJcAywrybeczNpRCCHE5OFCYF8z221gR8cR2vM6x6rot4zqG1Xq98+b2ceALwL/7e6lagIz+zuyhftvAv6BbHF9KWZ2aCf57M73S81sKbDU3a9w9+vN7NvAmWa2MXA78C5gV0LnuPv9ZnYGcIKZPQRcR9apB9BZJ9jJt9bMPkrmbOZu4NJOniOB49xdoWQnARNBtVkWDGRCqHGjxfIEdkCjgC5D0yf1+1B8CTgW+L6ZfYTMbuwU4E4yR2UAmNnOwJ+Ak9395Kx9/ZVR/aRKqD+JbNH8yWQC/vfAjcBSMhP9LckW7s8DHgdcCbzI3X9R47znJ9tndb6vAOZ30m8jcxhzKjCrc+6D3P26pOyJwMPAe4HtgFuA17n7D2Imd/+imTnwAeCDwB3Ase5+FkIIIVrDaAh1d19pZgcAnyXzWWJkRuHvc/eHQ1YDNmJDzXZfZVS/sKGCrpnZNODVZIvo9yVzQjMd+AtwM5kw/7a73zwSDWwL8+bN80XXXjvWzZh01LWALfsDiIO/ylFRnX+QirWvI7k2vGxE3veReuyDsv6YPn0wuWZd8ZwPh7/BtSX+GTfdtLg9c2aenrIuKMwefbSqpV3bE3+f8bJWf7yP1Oftsw+LFi2yfte7997z/Mc/XjR0xgq23dZ+PYSh3IRkSOv3jmr62wTPOUIIIcRIMUrq9wmJYtsJ0TKqRm9lx8b7iK+S9N99nIfknHC2EKJVjO+nQwghxIRDI/XmSKgLIYRoFRLqzZFQF6IFSO06eRgvRnxjiYR6cyTUhRBCtA4J9WbolVEIIYSYINQaqZvZL4D/IPMst3qo/EKMC8Z4KDDmKvcm118oM61wKBqll61Tr6puWnm2ehWME6v4Kut3kSH1e3Pq3lFrgK8D95jZGWa2xwi2SQghxCRmlAK6TEhqCXV3n08WVebrwFuAxWb2MzN7fcfvrRBCCNEXJNSbU1v34+43u/v7gScAbyXzhftfwF1mdlqMdCPEaDN1avdPber8M1Qcm8Jjw/6MCVX/esP8d6zq6403zj+iSLwLhOg3w76r3H21u59L5pz+KmBr4B+BP5jZ+Wa2XZ/bKIQQYhKhkXpzhiXUzWxTMzvSzK4FfgVsQybcdyALO/e3wDf63kohhBCTCgn1ZtS1ft8LeAdZnNjNgO8DH3L3y0O2L5nZfWwYVlUIIYSojazfm1N31vFG4B7gTGCBu99bku9W4Jo+tEuIdlC1ZKrOcqo2/zN1sxmok6/DlOSap07NFX+xSNW8eqGKXrsqnDRtW5vmr7WkbWgk1JtTV6gfCnzf3ddXZXL33wP799wqIYQQQgybWkLd3b8z0g0RQgghQCP1XhgfLpjEpKXu0q+oxiwtU6VuHq1/kDb9U1W1pU47kzxTp04L6e7ZKpcZlv0edftsnHiUk8p9aCTUmzM+ngIhhBCTCgn1ZkioCyGEaBUaqTdHQl1MCAoq97pW3HVUvP1W6ZbppceCKmv+Mir6bOr07ur3vlOmz6/Q87fJ4rxNbRETDwl1IYQQrUIj9eZIqAshhGgVEurNkVAXo8dYPKV1raib6ItjmV7Ljxa9xlDv9ZoTSqdN6v5ujz6ap6dPLy8T2trPYDpN1OdSuQ+NhHpzdHcJIYQQEwSN1IUQQrQKjdSbI6EueqdXlW4/KFMF96oirlA3RzXq+PkDKrdQn1JmSV5T/d5Ee75mXd6H03r1nx/zRbU8FFXzZTS4V/qpyhdFxs8z1S4k1IUQQrQKjdSbM+pz6ma2o5l93syuMbNVZuZmtkuSZ56ZLTCzmzt57jCzb5jZrl3qW9KpI/28qkveozt1rjazW8zsnSN3pUIIIZowINTbFk/dzKaY2QkdufOomd1oZq+pUW4LM/uYmf3CzP5iZg900q/qkvekEpn2vTptHIuR+u7A64BfA1cBL+6S5zBgT+BzwGLgCcBHgUVmtre735nkvwQ4Kdl3S9wws6OBs4FPApcCLwTOMjNz9//o5YKEEEJMCk4BjgdOJJNhhwHnm9nB7v6jinJzgP8HfK1Tx2PAG4Dvmtmx7v7vXcrsB8TIqMvrNHAshPqV7r4tgJkdRXeh/il3Xxp3mNnVwO3A0cDHkvzL3H1h2QnNbCrwCeBcdz+xs/tyM9sBOMXMvuzua5tdziSiQZCPRvmazI+XHCsEeqkoX7bMqDA3/GjFsXGiKqzuwikhHebeg6e4SHrNcRr74YeH37ZpM0uWpDW5p5p4y6vLOAkcM55po/rdzLYhE+inufvpnd2Xm9nuwGlAlVC/HdjN3VeFfZeY2U7Ah4BuQv2X7j7sXhh19bu7D2lZkgr0zr7/A5aSjdqHy3OBrYHzkv3nAo8neyMSQgjRAlqqfn8JmaVpKkfOA/bqNj2cX4+vTAT6AIuAHfrXxHG0Tt3MngpsA/y+y+FXdObeV5vZwi7zFHt2vm9K9i/ufM/tX0uFEEL0SguF+p7AauDWZH8vcuT5wM0lx+40s/Vm9n9m9ikz27ROheNCj9RRn3+RbKT+leTwRcCvyNQb2wLHks1TvNndB96otup8r0jKLk+Op+c9BjgGYM6cOb1cwvihSYztJk9QXVV6SFctIUtV492p9w5bV/O7dhxO2MQ2b7xx8dhwHcelfRNV7lEVX//2yH+fLWbOLM+WLlcboFe1eD/U/P1szySmT+r32Wa2KGwvcPcFPdS3FfCAu3uyv1KOlNGRL/sChyeHbgU+DFwPONkU9T8AzwJeNFS94+Wu+wLwt8DL3b0gmN39uLhtZt8FFpIZxKVqkmHRuQEWAMybNy/9IYUQQrSXZe4+r+ygmR0I/KRGPVe4+/y+tSo793wyQ/D/dPdvxGNhMDrAT8zsLuBMMzvQ3S+tqrv1Qt3MTiMbLR/h7j8eKr+7rzez84FPmdn27n4v+Qh9S+DekH3gzaqWVaEQQoiRZ5QM5X4BPLVGvoG58BXArM6KqTjIG5YcMbPnABcClwFH1WzrN4EzgeeQrd4qpdVC3cxOJLMMPM7dz21QxUDHD8x57ElRqA/MgfyuWQvHMU1U6f1Uszfw1Fal0m2iCk/Vz92oqnf9+vJjA2y0Uf32jDbptcX+6HWhQ5Vztzp1REv8Gakqvp/eA/t9v9dBavkhGQ2h3jFcK5vP7sZiYBPgiRTn1WvLETPbi2wJ9g3AaxqsuhpSY9xaQzkzew9wKnCiu39hGOWmAq8H7nD3+zq7rwGWAW9Ksh9O9nZ1de8tFkII0Q9aav1+MbCW7nLkJne/vaqwmT2JTN1/G3Cwuz8yjHMPnPPaoTKOySujmR3aST678/1SM1sKLHX3K8zsMDJVw8XAZWa2byj+oLv/rlPPG4BXkq0PvJPMUO7dZAYFbxgo4O5rzeyjZM5m7iZTXxwAHEmmBVgzMlcqhBBiIuDu95vZGcAJZvYQcB3ZAPIA4JCY18x+Cuzs7rt3trchE+jTgH8G5ppZLHK9u6/u5L0e+E8yB2pOZhx3HHCxu182VDvHSg90frJ9Vuf7CmA+cBBgne+DkrwDeSCzeN8G+DTZvMZKsnV/B7n7JbGQu3/RzBz4APBB4A7gWHc/i4nKBibifbRej1SZSpdYrxdU6RVOXaKKuErdXa7GzdOpKrzKErzb/lRdXaZar1LrjxfN60i185FkbBL7KlrPF/u2qFCcPn1G17obzQ4FxzpT1tV8t6+7PKKulXyhQePkBhkF2uZ8psOJwMPAe4HtyATv69z9B0m+jSjK17nAzp10mhdgV2BJJ30L2Squ7clu/tuAk4F/rdPAMbmD3N2GOP5W4K016llI9pZU97xnk7mKFUII0VLa6FEOMkNssmnhU4fINz/Z/hnZQLXOOQ5r2Dyg5YZyQgghJh9tFerjAQl1IYQQrUJCvTkS6hOBJnN8ZfubeHoLwT82CPIR5kn7OVdeRdk8ejrXXcdrWmT69KHzVLVlOMdGiiZxdCJTyEM3rFk3fA99Vb9BJN4TVUviys5TteQxtqF4DxQD10yP8+3huhu5yytrnObaRZ/RnSKEEKJVaKTeHAl1IYQQrUNCvRkS6uOFumr1JnGo6y5JC+rJsiVp6ZKlqEat27QmDsPK1OxVl1amTi+oWnul8nfrQx0D1FTPdo+MPnCeGucM+6cl5yzEQy/xa1UVZ72fjuLS5YerV3dPV50n3h/Tp08J6XxJXalaHvofE76MCaia10i9Oa31KCeEEEKI4THxXvGEEEKMazRSb46EepupUt/1U81eomKHokYxWrKXGQD3I+ZLHev1qjIxPW1qUI9uYJpfo+H9DhIyUvT7nD0G/Inx0GNwltQLX9nKhxjDJZ0mKfymgWiNnzazzJq+6nJivtiGeA2bb95dLQ+Jh7omU2TDDXA/1LFxhIR6cybGHSCEEGLCIKHeHAl1IYQQrUJCvTkS6m2jrpquiUV00CGWWbI/tKJYJFoK91Pjv8kmeTp1SBJVnWXlN7BQL2tcHRV7eiyeqIk1clm9demHCrXXf8Q66vea92QxHnrRNjfee/Gy4z0wjSTQysPdvdFEa/xZs4qq8Ace6FqkUhVfdiy2Mz4f8Z4G2HTT/BmbObMkcEyVxXyd/VVMEFW8GB761YUQQrQOjdSbIaEuhBCiVUj93hwJdSGEEK1CQr05EuptoM48el0XbHF5WrLEpmzuPM4L1p3iq7E6DijOM266aZ4uzJmmS5TqzI9X0eu8ZK/z6JGxmtfs53mb/LuWTJany77KKNwT6Rx6dEtX8ltNmV4ss9WsWYPpVY/m8/oPPZTnST3NDXfZZvrsxPqip8XNN8/n16fPLC4hnfLoqqFPOgnm1yXUmyOPckIIIcQEYXy9vgkhhJjwaKTeHAn1tlE37nLJ+p+CanFpsUiZmr1Ks1e2zKhsedrmmxfLxzKlHrZS9WofA5g0WntXtr/fKsy2BVQvI7azbnDzkaRw75REiEnd0IUycYnd9K3z6YC0qqgyr/PspF1T5pFu5co8vdlmxTKbb563Z8bMBjHc+zl1NMZIqDdjfP/qQgghJhwaqTdHc+pCCCHEBEEj9TZQ9kpaDOhcOBSDV6wIavao2qtryV52ynS7LMhGqYodiqr1fsZ9r3usiQpyIgfPqGpzHa968ceu8ihX0od1HSPG+zuN2156DVUq6pL7bcrMEIQmufmjF7jCypESi/kq9XtU7VfNZhRV83kflKrlu1UyQN0VMy1EI/XmtPuXFUIIMemQUG+OhLoQQohWIaHeHAn1saDKxLwktnkakCJuF2Keh/TatcUyMXBKmVq9EHsjyTdjeok17gM11Z5NqGN+X1Wm18gzgccmuAnKlDpTGFXm3g1UurG6eL/G/WlwloJjmbLpgKroLGUNSNTvU8L2jHBtM7YuWW0S1PJQrmaP6Whhn9YRVfZlanmALbfM+2eD4DfjGAn1ZkzsfykhhBBiEqGRuhBCiFYh9XtzJNSFEEK0Cgn15oy6UDezHYEPAfOAZwCbAru6+5Ikn5dU8Ux3vyHkm9Kp7x3AdsAtwMnu/j9dzn008AFgV2AJ8Fl3/2JvV9SAirnHuJRn2X35/nROPc69rV/fva50eVqcL3/c4/J09AJXmDeH4sTespperYZL1ZKlJnPqZXXVXOIz0efOyyi77imEe6Ls94BG90ScR4/Lw+I9nd77MThLefCfmvP9VWXKluiFfDPivPvWxf5YtXnen/Ea4iO1IgRWgmIfxHQsk9q9xL6aNSu3wyk8y+NMQrZVqA9H3nQpew5wRJdD/+bu70vy7gf8K/BM4K/AfwEnuvsjGxYvMhYj9d2B1wG/Bq4CXlyR9xzg7GTfH5LtU4DjgRM7dR4GnG9mB7v7jwYydQT62cAngUuBFwJnmZm5+380vhohhBB9pa1CnZrypoKlwCHJvnvjhpk9HfgJcAlwMNkg9NPAE4DXD3WCsRDqV7r7tgBmdhTVQv1ud19YdtDMtiHr4NPc/fTO7svNbHfgNOBHnXxTgU8A57r7iSHfDsApZvZld09sxYUQQoiMuvJmCNZUybQOHwfuAl47IJfMbA3wdTP7lLtfV1V41IW6uz82dK7avASYBpyX7D8P+KqZ7erutwPPBbbuku9c4G3AfsDlfWzXsCh4hwvquKiyS1WQZZrBquVps2fn6S2iV6oyFTvUC+BRd3lZXa9vNZb41XWWVao6TihTPY9jp1x9I/ZNoT/7QFxmWeWdLRKXkcXgLD0vX6zrba9sfVpSvhA4Zrt82Vl8ltNTlqnp68ZziUydXeGVbxzQwpF6XXnTGDPbGDgIOCMZaP438CXglUClUG/75OG7zGy1ma0ys8vM7O+S43sCq4Fbk/2LO99zQz6Am4bIJ4QQYowZUL/38hkB6sqbKrYxs2Vmts7M/mBmHzKzjcLxJwLTSWSVuz8K/KnOOdr8+nYe8APgHmBn4IPAZWb2Inf/WSfPVsAD7p4a1S0Px+P3iiHyFTCzY4BjAObMmdPgEoQQQjShD0rd2Wa2KGwvcPcFPdRXV96UcQPZPPxiMsH9ajIbrycBRyV1pLJq4DxDnaO9Qt3d3xw2rzKz75O9vZxKpi4fjTYsABYAzJs3r8waf9ikqt6oZvvrX/N06qEqEmOYb7llno5W7XE/wLR1q/KNZQ10e71apddVAYb6Yl81efvuVXU8DrWW7SP8cFOnTyscamKUHp+LGA89BmdpEmymchqo7Oar8mIXA8cEL3hbRbV80h/xUYqq+KrVLtErXfwvKXqNLP7nTJva32mUlrLM3eeVHTSzA8kM0obiCnef32tj3P3MZNePzOxh4H2dufI/9noOaLFQT3H3h8zsh8Dbw+4VwKyOBXsUugNvM8tDPoAtKVoapvmEEEKMOQ6UrNXtH78Anloj38BoqK68GQ7fBN5HtsT7jxRlVcpW5Kr+UsaNUA/EzlwMbEI2DxHnOQbmHX4X8kE2J3JvRT4hhBCtYGSFuruvAm4eRpG68qZRczrffyKbt98zHjSz6cBuwPlDVTRuhLqZbUG2Zu/asPtiYC3wJrJlAAMcDtwULBGvAZZ18l2a5FsOXD0SbS6zqE7ViWWBLDbdNE+nfj7K1OyFgA6pyXwdNXuvjmB6jV9Oucq9iSV6v621RaCutXhgyrpiwJGpYUVDHVU8lDtlKcRDTx+YMiv1Kmc6Zer3Oul0u8SUfUZyzhmzc9V8dCQT1e8xuEtSXem0xYaXNnIrGvrDqIzUh0tdeTMc3kR2sb8CcPc1ZnYx8DozO8ndB37RQ8leKC4cqsIxEepmdmgn+ezO90vNbCmw1N2vMLPjgaeQLTMbMJQ7nsyDz5sG6nH3+83sDOAEM3uIzNT/9cABhAX+7r7WzD5K5mzmbjLBfgBwJHCcu0+c0EZCCDEhaNfLRl15A2BmPwV2dvfdO9s7ky2h/hbZKH8TMkO5twJnu/ufQvGTgIXAf5vZvwO7kDmfucDdfz1UO8dqpJ6qEM7qfF8BzCdzvffqzudxwINko+m3u/u1SdkTgYeB95K77Xudu/8gZnL3L3Zcz36AzJL+DuBYdz8LIYQQYmhqyRtgI4ry9SEyrfCHgG3J3lhuBt5DLv8AcPcbzOzFwKeAH5K5if1P4J/qNHBMhLq72xDHLwIuqlnXejKL+FNr5D2bDd3OCiGEaBWtVL/Xljeptby7LwdeNYzzXEnmNG3YjJs59fFCk2AgcYlKdJBVtmwNYNqjD+YbvS5Pq1qSVmfuvO48egtdRIlRoOJ3j0u64q1W5rQt3Y7LuWbOzOuaUjfYTN0lbbG+ukFk6nihi3PtULCDmRHOX/Cct8usQpHlD+T/OWl13ZqSbrdzfr2dQn08IKEuhBCihUioN0FCXQghRMvQSL0pEuoN6He87TIN9xbTg1H+fSG4OhSXq9VdnhZVeGXLf1K1ZYl3t7o0UefVXcYWKQ3iIlpN/K2mT58S0nmeKtVxWRCYGXW9w1XdbGWeDQnL8KZXqPzLotJEHXl6znQZ6uCJSqI2AVttt13e5NlbDHn6KkYyeI8YHSTUhRBCtBC9VDRBQl0IIUTLkPq9KRLqLSBq0wpBV+4KKvcq73CRmqr0MlV8jFkOwzdY38AQvkYFTdT6CrTSQoYbAAUK9/H06Xlwlo1CMMr0ty4zKo9e12ZsXaEKL3t2KvT8UeVeprGPVuQA05pYzEfKTNkryszINfFMnZmr4nsNhpQy8qp5CfWmtD2euhBCCCFqovGOEEKIFqKRehMk1MeAVCteCHIR1exVVrJlTjOiWj2xki1TzUeVe5WlcZ3T11bLVejPm6gKm6gD+72KYTJQ2zq6bnSWGHM81Lf55vl5ooV7VXUx36pHi7/tjDLPNjUDskQr9/qrM/I2lKri02e0rMIqp1LxfyL8f0zbLi6rKU6rtR+p35sioS6EEKKFyPq9CRLqQgghWoZG6k2RUG9AVBP2xSFLmWVuVaznOs5jqsqMlPl4tZPpWkXqFK9sfo1zihZSYgkfYyAk2UrT0RIeEmv4Oo5okgqjY5nUyr0WVc9yGXWnCSIlMdynzZpVyKapp4mL/vGEEEK0EI3UmyChLoQQomVI/d4UCXUhhBAtQ0K9KRLqo0RhHr1i7q507q3unHrNiee6c2plVZQuY6u5Hq0QIKPmNHzl0jnFam8XdZe0lRyLyzw33bS4HKssiEvZYwSwKiyRK8Qmr7mkLVZeWJ5W8zkqLAWsml9vEnimTr7EC92UHgM1jQ6yfm9CW39NIYQQQgwTjdSFEEK0DKnfmyKhPoKUetyq8g5XpkqvGee8oLqu29CSplRRuLa6gZtD5b3GTBfjlKroLCXpmTOL6vdHHsnT8bGIjtXSWzI6apy+Xb5cbsr0kDENoFIj8Mu0imWipd73aq/NLKFq+q5usJpAr0t0Rw4J9Sbob1IIIUTL0Ei9KW16LRNCCCFED2ik3md6jjNcI+gK1NN4V2n5y/ZXtr+OlW16krC9Zl3+DlkVn6KWdrJKZ19SqN+qxTK1ZZU6c6KVqUWVirlG3PVCwCNg883zZ2Hlyu6nSZ+PkpgnbBUt4dNCj5ao5svu/WSKbEqJ57iqPqwzLTU1+S+YNrO3lQax49qlitdIvQkS6kIIIVqGoyVtzZBQF0II0UI0Um+ChHqP9Kxuh1KL96o455G68UsaxTkZblCJ1AI4XkOFyr2MRrHaR4ky9WSV2rJOmfQ666jCq8r0s219p46TmkQtPj1Yw2+2Wfds0UIeYMWK7qecHuKkz0jV5XVU7nVv5BJVfF0/MlVlmF4jbnsDxvZ5k6FcU8Z60kQIIYQQfUIjdSGEEC1DI/WmjPpI3cx2NLPPm9k1ZrbKzNzMdknynNTZ3+3zaJJ3SUm+V3U599FmdrOZrTazW8zsnSN7tUIIIZrxWI+f/mNmU8zshI7cedTMbjSz19Qot0uFTHMzOyzkLZN/36vTxrEYqe8OvA74NXAV8OIueb4MXJzs26yz78Iu+S8BTkr23RI3zOxo4Gzgk8ClwAuBs8zM3P0/hncJPVI1od3Aw1RfPa018VZV4VWrrPjatXl6442L1fUaLGakGKs5xrLzVrWnbfYHg/S4vA1gyqOrBtObb557h4vL2x56qFhFDAJT8C4Xbt0Zs8PytjRj2fx6XBJXsZwz0iQ2S1W1BVuZYCMwpSzTuKC1I/VTgOOBE8lk2GHA+WZ2sLv/qKLcvcBzu+w/FdiPTIal7EexE5bXaeBYCPUr3X1bADM7ii5C3d3vAu6K+8zszWTt/XqXOpe5+8KyE5rZVOATwLnufmJn9+VmtgNwipl92d3XlpUXQggx2rRLqJvZNmQC/TR3P72z+3Iz2x04DSgV6u6+GijIKDObAewDXOTuK7oU+6W7D/ttbNTV7+7edAhxBPBnur/RDMVzga2B85L95wKPJ3sjEkIIIcp4CTCNDeXIecBeZrbrMOv7e2Bzug9UGzMuDOXMbCdgf+DMkjeXV5jZKmAj4HqyN6nvheN7dr5vSsot7nzPBS7vX4tHnrKwyz3TRB8Y2MDzXUmQjY02ytM1tZbVjJAXudaqsScidW/qcGzGzPz32Wyz/LdOY7NE9XuZd7nZsxNPbWXtKZuGSu/BkmVskarHrarqWozrCEitVL/vCawGbk32Rzly+zDqOwK4nw2nmge4s6MduAv4FnCSuz9SkneQ8fKrH06mVej2RnMR8CuyztwWOBb4rpm92d0H3qi26nynKo7lyfECZnYMcAzAnDlzGjdeCCHEcOlZqM82s0Vhe4G7L+ihvq2AB9zdk/2VcqQbZvYE4ADg37oMVG8FPkw2QHWyKep/AJ4FvGiouseLUH8LcL27/yY94O7HxW0z+y7Z3MUn2VBNMiw6N8ACgHnz5qU/pBBCiBGhL25il7n7vLKDZnYg8JMa9Vzh7vN7bUzCm8kGquekB8JgdICfmNldwJlmdqC7X1pVceuFupntA+wBvK9Ofndfb2bnA58ys+3d/V7yEfqWZFaIAwy8WdWyKhwxohe5BuriXlXxlRbmdQK3lMRJh6KVexmV6vdRstqVmn2MqBNlqOqeDKrwMkt4KI+1Hq3kVyR6vG3LLNvLVPGpWr7k2Vm3rqjmL6NO11TlE0PyC+CpNfINLLdYAczqrJiKg7wmcuQtwA3dBqolfBM4E3gO2eqtUsbD7XAEsBb4rwZlBzp+YM5jT4pCfW7n+3fNmiaEEKL/jPycuruvAm4eRpHFwCbAEynOqw9LjpjZc8heJv5hGOceYEiNcavdxJrZNLJ1gP/r7ktrlpkKvB64w93v6+y+BlgGvCnJfjjZ29XV/WmxEEKI/rC+x0/fuZhsgNlNjtzk7nWN5I4A1jG8gerAOa8dKuOYjNTN7NBO8tmd75ea2VJgqbtfEbIeTKba6Gryb2ZvAF5Jtj7wTjJDuXeTGRS8YSCfu681s4+SOZu5m0x9cQBwJHCcu69J6x5R+uBgpo6VbO2Q43W9YURK4r6nGsj1Jc9WdDiTXlstVXhFH459HGjRM3Wd1NSwhIeij5hYPN6fqcqeXWbl6VkhXWb9nlISDGlqRaCmuir3iU/7rN/d/X4zOwM4wcweAq4jG0AeABwS85rZT4Gd3X33ZH8cqN7f7Txmdj3wn2QO1JzMOO444GJ3v2yodo7VrXJ+sn1W5/sKYH7YfwTZSPoHJfXcDmwDfJpM+K8EFgEHuXthPbu7f9HMHPgA8EHgDuBYdz8LIYQQLaNdQr3DicDDwHuB7cgE7+vcPZVRG9Fdvr6czDdK1dr0W8hWcW1Ppk2/DTgZ+Nc6DRwToe7uVjPfK4c4vpDsLanuec8mcxUrhBBCDAt3X0/m2vXUIfLNL9n/XaBS/rn7YVXHh2LSKXWEEEK0nb4saZuUSKi3jDifXDU3XGfuvGrpS+kytqpJvpgO8+iPVcwRRs9xcR69IgZMz8vYyvpQy9ZaTq8Tx2EOe8stZxQOldl2PBL8c6VT5csfyO+drWbPzg/E+zN1XRcpWfo2bVa0Ryk+43Vu/cnjXa6V6vfWMx5/aSGEEBOa9hnKjRck1IUQQrQMCfWmSKhPAOpq1grq55LlNpU6+5gOa4SqVviUqdxjurZavOaFSuU+Aah7U5foq6dRXKU6a1Z3L25V2vN4bPrsLQbTM7YLmWJEmKoHIR4LFU9Lgr7EeOjjLgS6aAUS6kIIIVqIRupNkFAXQgjRMqR+b4qEeovpWXXcRH9XFRM6HFuzrrtlfhpCukx7X3ltPVrqSuU+wai6H2q6UJwxPb8nps7O7914v9bVnk+dmavip20XKkh1+WXtqYjNPiVUF+O5T04viXqOmzAZ7xQhhBBiQqKRuhBCiJYh9XtTJNQnGnVV7mV68UTVWUftV1c7KrW46AtNpmfCcxHV2jNndlfFJ0XK9wfHS9Ni0BcoX2FSs52RwlPY4BmtHTmmgtFfVSKh3gQJdSGEEC1DI/WmSKgLIYRoIdLsNUGGckIIIcQEQSP1iUaPy8HqLp1p5MVOiJGkwb0/bWpY6ja1/N6vMyWePjtT0kn64VQ26ZH6vSkS6kIIIVqGhHpTJNSFEEK0EAn1Jkio98gGKrc66uYq9dsYxD3up7cqqdvFeCW9d+Nz0eSxLF0C1udnvNYz1+MSNjF+kFAXQgjRMqR+b4qEuhBCiBYirV8TJNQbUKWW6tnrUlVs85YilbuYiMT7uldV9Ih6YytEmyn5z6j4L2mnml0j9aaMD6khhBBikiGh3oQ2vqIJIYQQogEaqbeBMpX7CKrie1W5SeUuJhOtVsWXqd97DNQ0tkj93hQJdSGEEC1DQr0pEupCCCFaiIR6EyTUhRBCtAxHS9qaMaoTK2Z2qJn9j5n9n5k9Yma3mNknzWzzJN+WZvZlM1tmZivN7FIz26tLfdPN7NNmdm+nvmvM7Pld8k0xsxPMbImZPWpmN5rZa5pexxQeG/xU8VjIWZt16/JPzbqbnKfsGrrXuuFHCDF6VD7j8T+j5v+HaIaZvd/MLurIHDezk4ZZfj8z+0VHXt1nZmeY2aZd8u1pZj82s4fN7C9m9jUz26rOOUbbWuJ4Mp3KPwEHAf8BvAv4iZlNATAzAy7qHD8OeA2wMXC5me2Y1PcV4GjgY8DBwL3AJWa2d5LvFOAk4AvAS4GFwPlm9rL+Xp4QQoj+sL7Hz4hwNLAN8L3hFjSzpwM/Ae4nk1cfAd4GnJPk2wH4GbApcCjwbuBA4AcDcrKK0Va/v8Ldl4btK8xsOfB1YD5wGXAI8DzgAHe/HMDMrgFuB/4ReE9n3zOANwJHuvvXOvuuABYDJ3fqwcy2IXuZOM3dT++c93Iz2x04DfjRiF2tEEKIBrTWUG5Pd3/MzKYC7xxm2Y8DdwGvdfe1AGa2Bvi6mX3K3a/r5Psg2UD2Fe7+QCffPcAVwKuA71SdZFRH6olAH+BXne8ndL4PAe4ZEOidcn8lG72/MpQ7BFgLfDvkWwd8C3iJmW3S2f0SYBpwXnLe84C9zGzXZlfTA1OnDv8zgkitLkR96k5Rjdo01hj/f4wMA0K9XSN1d2/0o5nZxmTa5/8eEOgd/htYw4ay7YcDAr1z3iuBO5J8XWnDYsUXdL5/3/neE7ipS77FwBwzmxny3e7uq7rkmwbsHvKtBm7tkg9gbsN2CyGEEHV4IjCdRLa5+6PAn+jIoc78+q5pvg6LqSGvxlSom9kTyFTll7r7os7urYAVXbIv73xvWTPfVuH7AXf3IfIJIYRoDY/1+GkVA3KmTGYNHN8SsBr5ShkzvUxnxP19YB2ZsUDrMLNjgGM6m6tto426vT2J+swGlo11IyYA6sfeUR/2h6eMTLV/vQQumt1jJdPNbFHYXuDuCwY2zOxAMsO1objC3ef32JZRY0yEekfFcBGwG/ACd78rHF5BPhqPpG86K4CdK/ItD/lmmZklo/U03wZ0boAFnTYvcvd5ZXnF0KgP+4P6sXfUh/0hEZp9w90PGol6E34BPLVGvnSKtwkDcqtMtg1MBz9AZlBQlq9UXg0w6kK9YzBwATAPeJG7/zbJshh4cZeic4E73P3hkO/VZjYjmVefS2Z4cGvItwnZnMatST6A3zW9FiGEEOOTjty4eZRO9ycy2649404zm042uD1/oE1mtiTN12EumQV8JaPtfGYK8A3gAOBV7r6wS7YLgSeY2QtCuS2AV3SODXARmdn/a0O+qcDrgR+7++rO7ovJrOTflJzncOAmd7+9p4sSQgghKnD3NWSy6HUdOTXAoWSDzijbLgRebmaPG9hhZvuRaaZjvq6M9kj938mE8CeAlWa2bzh2V0cNfyFwDXCemX2QTG1xApnxwL8OZHb3683s28CZndH/7WSObHYlCHB3v9/MzgBOMLOHgOvIBP8BdNay12TB0FnEEKgP+4P6sXfUh/1hUvWjmc0DdiEfEM81s0M76R8NaI3N7CvAEe4eZexJZI7P/tvM/r1Tz6eBC9z91yHfp8kGnRea2SeBx5HJvl8C3x2yke4+ah9gCdl8QbfPSSHfVsBXyeYPVgE/BZ7Rpb5NgTOA+4BHOxc9v0u+jci89/wfmQrkN8Cho3nt+uijjz76jO8Pmfe3Mhm2S5qvS/nnkw1aHwX+DJwJzOiSby8yI76VZAPbc4DH12mjdSoQQgghxDinDc5nWouZ7WRmF5jZX83sQTP7jpnNGet2jTVmNr8TzCD9PJDk62tgnvGMme1oZp/vXNuqTn/t0iVf34MUmdnRZnazma22LIjScN1btoZh9GO3+9PTuBCTsR9tggTWEiWMtTqjrR9gBvBHMs8+ryJzz/dbMivGzca6fWPcN/PJ1E3HAfuGz7yQx4Cfk/k6fgOZi8QryNYG75jU9w2ypRxHAy8k8238CLD3WF9rn/vsz2SxBi4hUdcNty/I7FJWk8U12B84m8zjxsuSfEd39n+ik+/Uzva7xrpPRrgfHfhacn/uS6LqnIz9SGdel8z26AXA+zr33EJgSidP35/fun2tT4+/71g3oK0f4L1kDoR3D/t2JXOW8/6xbt8Y9838zp/mgRV5XtnJs3/Y9zgyO4nPhX3P6OR7W9g3FbgFuHCsr7WPfTYlpI/qJozq9gVZlKjVwMeT8j8FfpOUvR/4epLvq50/543Hul9Goh87xxw4dYi6JmU/Alt32feWTp8d0Nnu6/Nbt6/16f0j9Xs5hwAL3X1wbbtny9+upoZTfdH3wDzjGq8XCKLfQYqeC2zdJd+5wOOB/YZzDW2gZj/WZVL2oyuw1oRGQr2cqsAyCgKT8Q0zW29mfzGz/0rsDfodmGcy0O8gRQMOLNLfYbIEM3pXZ/57lZldZmZ/lxxXP+YosNYEQUK9nKqAMd1c+E0m/gp8hkz9eQBwCnAgcI1l8euh/4F5JgP9DlJUFkRiMvTtecD/I7svjyEbUV9mZvNDHvUjCqw10RiPgXbFGOPu1wPXh11XmNmVwLXAe8h8AggxZrj7m8PmVWb2fbKR56mMM3X5SGLjILCWGB4aqZezgnKn+t3eTCc17n4d8AfgOZ1dVf03cLxOviEDGEwg6vbFCjpBimrko0udk65v3f0h4Ifk9ydM8n60YmCtl3jzwFr9vGdFj0iol7OYcqf6CgJTzoB6rar/0sA8u5rZjC75YmCeyUDdvohBitJ8kN+fA/OV6e8wmYMZRfXvpO1HKwbWepl3D6zVz+e3bl+LHpFQL+dCYF8z221gR8fJxfOo4VR/smGZT+SnkKngof+BeSYD/Q5SdA3Zkqtu+ZaTreSYFHTuvYPJ70+YpP1oCqw1odGcejlfAo4Fvm9mHyF7wz8FuJPMacKkxcy+QRZA5zoypxPPJAu6czfwuU62vgbmmQhYHvjh2Z3vl5rZUmCpu19Rty+8ZpAid19rZh8FzjKzu4FLO3mOBI7zLHLUuGOofjSz48leMC8H7iGLbnU8sB3qRxjfgbXEUIz1Qvk2f4A5wP8ADwIPAd+ji6OLyfYhe7h/Q2YFv5bsRWcBsH2Sr6+Becb7h/JAED8bbl8wjCBFwDvI7B1Wk3lJ/H9j3Rcj2Y9ko8mryUbXa4G/kAmpfdSPDgqsNaE/CugihBBCTBA0py6EEEJMECTUhRBCiAmChLoQQggxQZBQF0IIISYIEupCCCHEBEFCXQghhJggSKgLIYQQEwQJdSFajpltZmb3BE9qvda3qZnda2av60d9Qoj2IKEuRPv5AJl3tP/pR2Xu/giZq89/6bj2FEJMECTUhWgxZrYJcBxwtvfX/eM5wE7Aq/tYpxBijJFQF2IE6ajObzaza+Oo2MxebGaPmdm7h6jiVWQ+uL+d1HuOmd1lZvPM7Bdm9oiZ3WJmL+8cf7+ZLTGzB83s+2a2dSzv7iuAS4Cj+nGdQoh2IKEuxAji7iuBNwDPIIvyh5ltC/wncJG7//sQVRwE/N7dl3U5tkWnni+TjbjvB/7HzD4D7A+8G3hfJ93tPFcCLzCz6cO8LCFES1HoVSFGGM/CU34YON3MLiULA7oeeHuN4vuShansxubAO939SgAzuwe4kSxu+Fx3X9/Z/zTgODPbaGBfh+uBacCzgF8M/8qEEG1DQl2I0eFM4EXAD8gE6YtKRt8pO5CpybuxckCgd7i5831pIrxvJnvWtwfuCvuXhnMIISYAUr8LMQp0jNzOBTYBbnT3n9YsOp0s9nQ3HkjOsaaTXJHkG9ifqtkf6XxvWrMtQoiWI6EuxChgZtsB/0amSn+Gmb23ZtG/AFuOULO26nzX0RgIIcYBEupCjDBmZsDXyUbcB5Kp4j9lZk+vUfxmYLcRatqune9bRqh+IcQoI6EuxMjzfjJhfnhnKdmHgd8B3zSzoVTfVwLzzGwkntW/Ae5299tGoG4hxBggoS7ECGJmzwL+Bfiku18Bg3PfbwB2Ac4YoopvA48D/m4Emncw8K0RqFcIMUZYf51UCSH6jZn9DLjV3fvmKMbM/oZsGdtT3f0P/apXCDG2SKgL0XLM7HnApcDu7n53n+r8LrDC3Y/sR31CiHYg9bsQLcfdrwb+Adi5H/V15vFvAE7sR31CiPagkboQQggxQdBIXQghhJggSKgLIYQQEwQJdSGEEGKCIKEuhBBCTBAk1IUQQogJwv8Hxxj+oeSPEB4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "scale = np.max(p.data[1])\n",
    "plt.imshow(p.data[1].T/scale,\n",
    "           origin=\"upper\",\n",
    "           vmin=-1, vmax=1,\n",
    "           extent=[0, grid.extent[0], grid.extent[1], 0],\n",
    "           cmap='seismic')\n",
    "plt.colorbar()\n",
    "plt.xlabel(\"x (m)\")\n",
    "plt.ylabel(\"y (m)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, the wave is effectively damped at the edge of the domain by the 10 layers of PMLs, with diminished reflections back into the domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert(np.isclose(norm(v[0]), 0.1955, atol=0, rtol=1e-4))\n",
    "assert(np.isclose(norm(v[1]), 0.4596, atol=0, rtol=1e-4))\n",
    "assert(np.isclose(norm(p), 2.0043, atol=0, rtol=1e-4))"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}