{
 "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",
    "class MainDomain(SubDomain):  # Main section with no damping\n",
    "    name = 'main'\n",
    "    def __init__(self, PMLS):\n",
    "        super().__init__()\n",
    "        self.PMLS = PMLS\n",
    "            \n",
    "\n",
    "    def define(self, dimensions):\n",
    "        x, y = dimensions\n",
    "        return {x: ('middle', self.PMLS, self.PMLS), y: ('middle', 0, 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",
    "    \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",
    "\n",
    "main_domain = MainDomain(nbpml)\n",
    "left = Left(nbpml)\n",
    "right = Right(nbpml)\n",
    "base = Base(nbpml)"
   ]
  },
  {
   "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))\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=6, staggered=NODE)\n",
    "v = VectorTimeFunction(name='v', grid=grid, time_order=1,\n",
    "                       space_order=6)"
   ]
  },
  {
   "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": "iVBORw0KGgoAAAANSUhEUgAAAhEAAAGBCAYAAADYEOPMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmcZGV59//P1fs6+wIDs7LJAIJmROIGogkaDZhgNCpuiWJiTDSGJ5GfCTFIFOMSf2pUMHncMIlKMBBFRdlUBA2yybDNDLMwbNOz9VrdXV19PX+cc3qKmurqU9W1nvq+X696dfepc07dZ6qm6qrrvu77NndHREREpFgttW6AiIiINCYFESIiIlISBREiIiJSEgURIiIiUhIFESIiIlISBREiIiJSEgURIiIiUhIFESIiIlISBREiIiJSEgURIiIiUpK2Wjeg3i1btszXrVtX62aIiIhUxa9+9au97r48zr4KIuawbt067rzzzlo3Q0REpCrMbGfcfdWdISIiIiVRECEiIiIlURAhIiIiJVEQISIiIiVRECEiIiIlURAhIiIiJWmIIMLMjjazz5rZ7WY2ZmZuZutiHttiZheb2Q4zGzeze83s/Mq2WEREJPkaIogAjgVeBxwAflrksR8GPgR8DnglcAfwbTP7nXI2UEREpNk0ymRTP3H3lQBm9g7gt+McZGYrgIuAy939E+Hmm83sWOBy4PpKNFZERKQZNEQmwt2nSzz0HKADuCpn+1XAKWa2fl4NExERaWINEUTMw0nABLA1Z/vm8OfG6jZHREQkOZIeRCwBDrq752zfn3W/FCkz7dz/+CCZ6dx/VhERaSZJDyJKYmYXmtmdZnbnwMBArZtTd75z9+O8+rM/4+xP3sITB1O1bo6IiNRI0oOIA8AiM7Oc7VEGYj95uPuV7r7J3TctXx5rNdSmcteuA7S3Gjv3jXHjg0/XujkiIlIjSQ8iNgOdwDE526NaiAeq25xk2PzEEJvWLqG/s41Hnh6pdXNERKRGkh5E/ABIA2/K2X4BcL+7b69+kxrbVGaah54c4qRVCzhuZR9b9gzXukkiIlIjjTJPBGb22vDX3wh/vtLMBoABd7813GcK+Kq7/zGAu+8xs08BF5vZMHAX8HrgbODcql5AQmwbGGViapqTjlrA8PgUP1Z3hohI02qYIAL4ds7fnw9/3gqcFf7eGt6yfRAYAd4LHAE8DLzO3b9bmWYm2+YnBgE4edVC9o1M8s07H2PfyARL+zpr3DIREam2hgki3D23ODLWPu6eAS4LbzJPDz45RGdbCxuW9/Hk4DgAW/aMKIgQEWlCSa+JkDJ7YnCcVYu6aW0xjlvZB8CWp1UXISLSjBRESFEGhidYHmYdjljQRWdbC7v2j9W4VSIiUgsKIqQoe4cnWL4gCCLMjGV9newbmaxxq0REpBYUREhRsjMRAMv6OxkYmahhi0REpFYUREhsqckMwxNTrFhwKIhY3tfBXmUiRESakoIIiW1gOMg4ZGcilvZ2sk+ZCBGRpqQgQmIbGAmGdC7vz+7O6GDf6CTTWtFTRKTpKIiQ2PYMBRmHFf1dM9uW9XWSmXYOptK1apaIiNSIggiJLSqgzM5ERJNMqUtDRKT5KIiQ2AaGJ2htMZb0dsxsW9YX/K4RGiIizUdBhMQ2MDzB0t4OWlsOzS4eFVlqhIaISPNRECGx7RmeeEZXBqg7Q0SkmSmIkNjyrda5qLud1hZjr4IIEZGmoyBCYhsan2JRd/sztrW0GEt7O9g7rO4MEZFmoyBCYhtMpVnQffjq8Uv7Otk3qkyEiEizURAhsbg7g6k0C3MyERB0aQylpmrQKhERqSUFERLL6GSGzLTnDSIWdLcxqMmmRESajoIIiWUoDBIWdOUJIrraGRpXECEi0mwUREgsUaYhfyaifSbIEBGR5qEgQmIpGER0tTM6mWEqM13tZomISA0piJBYoiBiwSw1EQDD4yquFBFpJgoiJJahOTIRgOoiRESajIIIiaVwJiIMIjTMU0SkqSiIkFiGUmnMoL/z8MmmFnQF25SJEBFpLgoiJJbBVJoFXe20ZK3gGVnYE2UiFESIiDQTBRESy9D4VN4pr0E1ESIizUpBhMQy25TXoJoIEZFmpSBCYikURPR2tNJiaOprEZEm0xBBhJmtNrOrzWzQzIbM7BozWxPz2DVm9lUz22VmKTN7xMwuM7PeSrc7SQoFEWYWzFqp7gwRkaaSv5O7jphZD3ATMAG8FXDgMuBmM3u2u48WOLYX+DHQDvwdsAt4HvAPwHHA6yvb+uQYCgsrZ7OgS1Nfi4g0m7oPIoB3AhuAE9x9K4CZ3QdsAd4FfKrAsS8kCBbOcfcbwm03m9kS4CIz63H3sco1PTmGxtN554iILOhuY0gzVoqINJVG6M44F7gjCiAA3H07cBtw3hzHdoQ/h3K2HyS49sPHK8ph0plpxtPT9OWZIyKiTISISPNphCDiJOD+PNs3AxvnOPbHBBmLj5nZRjPrM7OzgfcCXyzUFSKHjE4EGYZCQcRC1USIiDSdRggilgAH8mzfDywudKC7jwMvIrjOzcAwcCPwXeA95W1mco3ECCL6u9o0OkNEpMk0Qk1EycysC/gmsAJ4M0Fh5enAJcAU8KezHHchcCHAmjWxBoEk2uhEBoDeAkFEX2f7zH4iItIcGiGIOED+jMNsGYpsfwycBRzr7tvCbT8xs0HgSjP7orvfm3uQu18JXAmwadMmL7XhSTEyEWQY+rpmf7n0drYyOjmFu2OmUhMRkWbQCN0ZmwnqInJtBB6Y49hTgANZAUTkl+HPE+fZtqYwEmYY+jpbZ92nt7MNdxibVDZCRKRZNEIQcR1whpltiDaY2TqC4ZvXzXHsU8BiMzs2Z/vzw5+Pl6mNiRYVVhbqzojui/YVEZHka4Qg4kvADuBaMzvPzM4FrgUeA66IdjKztWY2ZWaXZB37FYJiyuvN7K1m9lIz+z/AJ4BfEQwTlTmMhPM/9HYUqokIshQjCiJERJpG3QcR4TDMs4FHgK8D3wC2A2e7+0jWrga0knVN7r4DOAO4h2CWy+sJJq+6Evgtd5+uwiU0vCgw6C9UE9ERZSLUnSEi0iwaobASd98FnD/HPjvIM3mUuz8AvK4yLWsOcbozouGfykSIiDSPus9ESO2NTE7R0dZCe+vsL5cowBibVBAhItIsFETInEbGp+gvkIWAQ0GEMhEiIs1DQYTMaXRiqmBXBhzqzlBNhIhI81AQIXMamcjMGUT0hqMzNMRTRKR5KIiQOY1MpAtONAWHRmeoO0NEpHkoiJA5jU5kCi6+BdDSYvR0tCoTISLSRBREyJzi1ERAUFw5qtEZIiJNQ0GEzGl4YmrOTAQExZUjKqwUEWkaCiJkTqMxg4jeTnVniIg0EwURUlBm2hmbnHt0BgTFlSqsFBFpHgoipKCoxiFud4YyESIizUNBhBQUZ92MSK+CCBGRpqIgQgo6FEQUnici2EeFlSIizURBhBQUTWMdTSZVSK/miRARaSoKIqSgsckgiOiJmYlIpTNkpr3SzRIRkTqgIEIKipb27omRiZhZhEsTTomINAUFEVJQlIno7YiXiQAYU12EiEhTUBAhBUWZiO4YQURPuM+YMhEiIk1BQYQUdCgTMXd3RvdMEKFMhIhIM1AQIQUVU1jZ3R7sk0oriBARaQYKIqSgsckpWluMjta5Xyo9ykSIiDQVBRFS0OhEhp6OVsxszn2j7oyUaiJERJqCgggpKDWZmckwzCUaBqpMhIhIc1AQIQWNTk7FKqoEdWeIiDQbBRFS0NhkJtbwTsjuzlAQISLSDBRESEFjxWQi2pWJEBFpJgoipKBiMhFtrS10tLYwllZhpYhIM1AQIQWNTWZiLQMe6e5oVXeGiEiTaJggwsxWm9nVZjZoZkNmdo2ZrSni+BPN7NtmttfMUmb2sJm9t5JtToKxiSm62+N1Z0BQXKnuDBGR5hD/06GGzKwHuAmYAN4KOHAZcLOZPdvdR+c4flN4/C3AO4BB4Digr4LNToSxtDIRIiKSX0MEEcA7gQ3ACe6+FcDM7gO2AO8CPjXbgWbWAnwNuNHdfy/rrpsr19zkGJvIxFoGPNLT0appr0VEmkSjdGecC9wRBRAA7r4duA04b45jzwJOpECgIfmlM9NMZqZjTzYF0NPeplU8RUSaRKMEEScB9+fZvhnYOMexLwp/dpnZHWaWNrM9ZvYZM+suaysTZmbxrSKCCHVniIg0j0YJIpYAB/Js3w8snuPYVeHPbwI3AL8F/BNBbcS/l6uBSZSaCSJUWCkiIodrlJqI+YgCpavc/ZLw91vMrBW43MxOdPcHsw8wswuBCwHWrIk9ACRxRsNuiWILKxVEiIg0h0bJRBwgf8ZhtgxFtn3hzx/lbL8h/Pmc3APc/Up33+Tum5YvX15UQ5NkbCIIBrrbiwgi2lVYKSLSLBoliNhMUBeRayPwQIxjC5kuqUVNYGwmE1Fsd4YKK0VEmkGjBBHXAWeY2YZog5mtA14Y3lfI9wnmlzgnZ/srwp93lqeJyRN1S8Sd9jrYt43x9DTT016pZomISJ1olCDiS8AO4FozO8/MzgWuBR4Droh2MrO1ZjZlZlHtA+6+D/go8Cdm9hEze7mZfQC4BPhq9rBReaZSRmdE+6pLQ0Qk+RqisNLdR83sbOCfga8DBtwIvM/dR7J2NaCVw4OjS4Fh4N3ARcCTwMeBD1e46Q1tpjujyNEZwbGZorpBRESk8TTMu7y77wLOn2OfHQSBRO52J5hsShNOFSHKJnQVWVgJaK4IEZEm0CjdGVIDqZK6M4K4VMuBi4gkn4IImVUpmYjs7gwREUk2BREyq9Rkhs62FlpbDushmlU0kkPdGSIiyacgQmaVSmeKGt4JykSIiDQTBREyq7HJDD1FdGWAhniKiDQTBREyq1Q6Q1eRmYjusLAypVkrRUQST0GEzCo1mSlqZAYwk7lQd4aISPIpiJBZpSYzRS2+BYcKKxVEiIgkn4IImdVYOjPTPRFXZ1sLLabRGSIizUBBhMxqfDJDd3txLxEzo7u9VZkIEZEmoCBCZjWWnpqZgbIY3R1tpDRjpYhI4imIkFmlJqeLmq0y0tOhTISISDNQECGzSk1OFV1YCQoiRESahYIIycvdSaWLH+IJwQgNFVaKiCSfggjJa2Jqmmmn6GmvIcpEqCZCRCTpFERIXuPhtNWldGd0t7epO0NEpAkoiJC8orUvSs1EaO0MEZHkUxAheUWZhFJqIlRYKSLSHBRESF5RYWQpQzxVWCki0hwUREheUXdEqZmIVDqDu5e7WSIiUkcUREheUSahtHki2shMO5OZ6XI3S0RE6kjsIMIC55rZJ8zsy2a2Ntx+ppmtqlwTpRaimoZSCiujwENdGiIiyRZrYQQzWwxcDzwfGAb6gM8CO4F3AvuBv6hQG6UG5jPEsydrOfBFPWVtloiI1JG4mYiPA6uBFwJLAcu678fAy8rcLqmxQ6MzSlmA61AQISIiyRX3E+I84CJ3v93Mcr+a7iIIMCRBUvPKRAQvK3VniIgkW9xMRB/w+Cz3dfHMzIQkQCqctno+NRGa+lpEJNniBhEPA789y31nAr8uT3OkXqTSGVpbjPbW4uPDme4MzVopIpJocbszPg98zswGgX8Pty0ys7cD7wEurETjpHbGJjN0t7diVnwQERVWqjtDRCTZYmUi3P1K4FPAPwBbw80/Aq4EPu3u36hM8w4xs9VmdrWZDZrZkJldY2ZrSjjPB8zMzexnlWhnUoynMyV1ZcAzR2eIiEhyxS69d/cPmNkXgN8CVgD7gB+5+6OValzEzHqAm4AJ4K2AA5cBN5vZs919NOZ5NgB/C+ypVFuTIspElKJ7JhOhmggRkSQravyeu+8E/rVCbSnkncAG4AR33wpgZvcBW4B3EWRJ4vgC8A3gBIq89maTmsyUNOU1HBqdoUyEiEiyzfpBWmxXgbvvmn9zZnUucEcUQISPt93MbiMYfjpnEGFmbwSeC7wBuKZSDU2KVDpT0uJbAD3RjJUqrBQRSbRC38Z3EHQbxFXaJ048JwHX5tm+GfiDuQ4OZ9z8Z+Cv3X1/KcWCzWY+mYiWFqOzrUWFlSIiCVcoiPgjDgURnQS1BEPAt4CngSOA1wH9wIcr2EaAJcCBPNv3A4tjHP9x4BHgK2VsU6Kl0hkWdreXfHx3uJKniIgk16xBhLt/JfrdzD4N3AX8nmet72xmlwL/DWysYBvnxcxeDLwFeK7HXJvazC4kHLa6Zk3RA0ASITVZ+ugMCCacUiZCRCTZ4k429QbgitwP4fDvLwJvLHfDchwgf8ZhtgxFtiuAfwN2m9kiM1tEEDy1hn935h7g7le6+yZ337R8+fL5tr0hpdKlj86AIBOhyaZERJIt7giFPmC2T9MVQG95mjOrzQR1Ebk2Ag/MceyJ4e1P8tx3APhL4NPzal0Cjc2jJgKCTMS4MhEiIokWN4i4BfiImT3o7v8bbTSz04F/DO+vpOuAT5jZhmheCjNbR7Cq6AfmOPalebZ9mqAQ9M85NHmWZEmlM3TNI4jo6WjVEE8RkYSLG0S8h2DJ7zvM7DGCwsqVBKt3bg/vr6QvhY9xrZn9LUHB54eBxwi6KwAws7XANuBSd78UwN1vyT2ZmR0E2vLdJ5CZdianpulpL30qja72VobHNdmUiEiSxZ32ejvwLIIugRsJZqu8kWCipxPdfUelGhg+/ihwNsEIi68TTBi1HTjb3UeydjWCDEPcWg/JY2YZ8I7S/xm721sZV02EiEiiFTPtdZogI/ClyjWn4OPvAs6fY58dxFiW3N3PKk+rkmlsZhnw0jMR6s4QEUk+fWOXw4xPTgPMe3SG5okQEUm2WF81zWw7hWevdHc/pjxNklobSweZiPmNzmjTPBEiIgkXN199K4cHEUuBFwAjBCtsSkJEH/7zy0S0kEpncHc0zbiISDLFCiLc/W35tocTN/2AYOSGJEQURJS6ABcEAUhm2pnMTNPZVsllVUREpFbmVRPh7gcJ1qW4pDzNkXoQ1TLMqzsjLMqM6itERCR5ylFYOQ4cXYbzSJ2IRlXMd+0M0HLgIiJJVvIYPjNrA04GPkQwLbUkxMw8EfPozoiyGNFwURERSZ64ozOmmX10xhDwqrK1SGpuPD3/TESXMhEiIokXNxNxKYcHEePATuD77j5Y1lZJTUXdGfOpiYiO1TBPEZHkijs640MVbofUkZnRGfMYVRFlMZSJEBFJrliFlWZ2k5k9a5b7jjczzRORIKl0hq72FlpaSp/fYaawUpkIEZHEijs64yxgwSz39QNnlqU1UhdSk5l5FVWCMhEiIs2gmCGesxVWHkMwa6UkxNhkhp55LL4FykSIiDSDWT8pzOztwNvDPx240syGc3brJhjmeWNlmie1MB52Z8zHoSGeCiJERJKq0CfFNJAJb5bzd3TbB3wB+OPKNlOqaWxyat6ZCA3xFBFJvlk/Kdz9q8BXAczsZuBP3f2hajVMaieVnn9NRGdbC2bqzhARSbK4QzxfWumGSP1ITWZY1NMxr3OYGT3trcpEiIgkWKGaiLcA33P3feHvBbn718raMqmZsckMqxbNf+XN7g4FESIiSVYoE/EV4AyCuoevzHEeBxREJMRYGYZ4QhhEqDtDRCSxCgUR64Ens36XJjGezsxr3YxId7uCCBGRJCtUWLkz3++SfME8EeXIRLQxpu4MEZHEmt9kAJI409MejM6Y5xBPgO72FsaViRARSaxChZXbmX2Wylzu7seUp0lSS+NT4TLg5aiJaG9l78jkvM8jIiL1qdDXzVuJH0RIQpRjGfBIT0cbY5Nj8z6PiIjUp0I1EW+rYjukTkSFkOUorOxqb2U8PT3v84iISH1STYQ8QzSvQ3kyEa2MTU7N+zwiIlKfYgcRZnacmX3VzB4xs9Hw51fM7NhKNlCqK+rOKNs8ERqdISKSWLFK8M3sLOB6IAV8D3gaWAn8LvB6M3uFu99aqUZK9USZg3J2Z0xPOy0tNu/ziYhIfYmbifgkcDew1t3f4u7/x93fAqwD7gnvrygzW21mV5vZoJkNmdk1ZrYmxnGbzOxKM3vIzMbMbJeZfcPMNIFWHuMz3RnzH+IZdYlEIz5ERCRZ4gYRG4GPuftI9kZ3HwY+BpxU7oZlM7Me4CbgWcBbgTcDxwE3m1nvHIf/Ydi+zwCvBD4APBe408xWV6zRDaqcozOiLhHNWikikkxxv27uBmZb1rEDeLw8zZnVO4ENwAnuvhXAzO4DtgDvAj5V4NiPuftA9gYzuw3YHp73koq0uEGVuyYiOufSeZ9NRETqTdxMxMeAfzCzVdkbzewo4O+Bj5S7YTnOBe6IAggAd98O3AacV+jA3AAi3LYTGACOKnM7G145h3hGgci4iitFRBIpbibiTGAB8KiZ3cGhwsozwt/PCosvIZi98q1lbudJwLV5tm8G/qDYk5nZicAK4MF5titxyj3EEw5lN0REJFniBhEvAqYIVvVcG97g0CqfL87atxKzXC4BDuTZvh9YXMyJzKwN+CJBJuLf5t+0ZIk+8LvaylgToUyEiEgixQoi3D1JIxk+B7wAeJW75wtMMLMLgQsB1qyZcwBIoqQmp+hqbynLkMyuDgURIiJJ1igzVh4gf8ZhtgxFXmZ2OUFw8EfufsNs+7n7le6+yd03LV++vOjGNrJgGfD5D++EQ90ZGp0hIpJMRX1ahEMiVwNdufe5+03lalQem8k/jHQj8ECcE5jZB4G/Af7c3b9exrYlSiqdKcvIDNAQTxGRpIs7Y+UG4BvA6dGm8KeHvztQnk+e/K4DPmFmG9z90bBN64AXEsz7UJCZ/QVwGfBBd/9cBdvZ8FKTmbIUVULWEE91Z4iIJFLcTMS/AmuA9wEPAZMVa1F+XwLeA1xrZn9LELR8GHgMuCLayczWAtuAS9390nDbHwKfBn4A3GRmZ2Sdd8jdY2UymsVYOYOIaIinMhEiIokUN4h4HvA2d/+vSjZmNu4+amZnA/8MfJ0g+3Ej8L6cWTSNICOSXevxinD7K8JbtluBsyrU7IaUmszQVebuDA3xFBFJpmJmrKx29uEZ3H0XcP4c++zgUFdLtO1twNsq1a6kSaUzLOubbXLS4rS1ttDR2qLRGSIiCRV3dMZHgL+JsU6FNLixyamyjc4A6Gpv0YyVIiIJFXeeiK+b2bOAHeGMlbnDKisxS6XUQGoyU5YpryM9HW0zy4uLiEiyxB2d8TbgYiBDsAJmbtdGJWaplBoYK+MQTwhGaKTS02U7n4iI1I+4eet/AL4D/LG7H6xge6TGyjk6A4LiypQyESIiiRS3JmIp8HkFEMmWmXYmp6bL2p0RZCJUEyEikkRxg4ifASdWsiFSe+VcwTMSZCIURIiIJFHc7oz3At8yswMEkzYdtl6Fu6vju8FFBZDlronYOzJRtvOJiEj9iBtEPBj+/FqBfSo57bVUQZQx6C7jEM/udnVniIgkVdxPi0vRCIzEq0R3Rk9Hq2asFBFJqLjzRHxotvvM7CzgLWVqj9TQ2EwmosyFlQoiREQSKW5h5TOY2bFmdqmZbSdYw+J15W2W1MJMd0YZayJ6w8mm3JXIEhFJmthBhJktNLMLzew24GHggwQFlu8GVlWofVJFUSairKMzOlqZdpiYUt2tiEjSFAwizKzFzH7HzL4JPAl8EVgL/Eu4y/vc/Qp3H6pwO6UKKlUTAahLQ0QkgWYNIszsk8DjwP8AryaYsfIVwBrgEnJWy5TGF80sWc7RGVEQMaYRGiIiiVPo0+IvCUZkXA+8zd33RXeYmTq4E2isAjURUUAyNqGpr0VEkqZQd8a/AcPAq4CHzexzZnZ6dZoltVCJmojeKBOh7gwRkcSZNYhw93cCRwBvAu4E3gXcbmYPAn+D5o1InPF0BjPobCtp0E5e3QoiREQSq+CnhbuPu/t/uHtUCxEtB/4BgpqIy83sAjPrqnxTpdLGJjP0tLdiVr5yl56wOyOVVneGiEjSxP7K6e5Puvs/ufvJwOkEIzSOI5gK+8kKtU+qaGwyU9aiSsgqrFQmQkQkcUrKW7v7ne7+5wTzQ5wP3FLORkltpCan6O4oX1cGZAUREwoiRESSZl5fO909TTD08zvlaY7UUtCdUe5MRDg6Y1LdGSIiSVPer53S0FLpTFnXzQDNEyEikmQKImRGajJT1uGdEIz0MNOMlSIiSaQgQmaMTWbKOtEUgJnR067lwEVEkkhBhMyoRHcGQE9nm2oiREQSSEGEzKhEdwYEdRHKRIiIJI+CCJkxNjk1M5qinLrVnSEikkgKImRGKp2hq8w1ERBkIlRYKSKSPA0TRJjZajO72swGzWzIzK4xszUxj+0ys4+b2ZNmljKz283sJZVucyNJZ6ZJZ7wi3Rm9qokQEUmkhggizKwHuAl4FvBW4M0EU27fbGa9MU7xb8A7gUuAVxNM0/1DMzutMi1uPKl0+VfwjKg7Q0QkmcrfAV4Z7wQ2ACe4+1YAM7sP2EKwuuinZjvQzE4F3gj8kbt/Odx2K7AZuBQ4t7JNbwxRd0NFRmeosFJEJJEaIhNB8EF/RxRAALj7duA24LwYx6aBb2YdOwX8J3COmXWWv7mNJ/qQr0gmoqNNQYSISAI1ShBxEnB/nu2bgY0xjt3u7mN5ju0Ajp1/8xrf6ERQs9BbgdEZQWGlaiJERJKmUbozlgAH8mzfDyyex7HR/VWxY+8oAyMTrOzvYsWCzoqMhCjVTBDRWf6XRG9HK2PpDO6OmZX9/KUam5xi94EUQ6k0w+NTDI2nmUhPk3Fn2p3paScz7Uw7wd/utW6yiMhhlvZ2cv5vHF2Tx26UIKKqzOxC4EKANWtiDQCJ5epf7eZzN8/0yLBheS+nHb2IMzYs5ZyTj2Bhd3vZHqtYle7OcIfx9HRFai7i2n1gjJsf2sOtjwxw3+5B9gxP1KwtIiLlsvHIBQoi5nCA/BmH2bIMuceuneVYOJSRmOHuVwJXAmzatKlsXz/f+Pw1nL5+CU8NjfPEwRT3Pz7ET7bs5Zq7H+eS6+7nDaev4X0vO56FPdUPJkYnK5iJ6AxX8pycqkkQ8ejACB/7wUPc8MDTuMPqJd285PjlrF/Wy+olPSzuaaevs43+rja62ltpbTFazTCzQ7+3QIvXM/ReAAAgAElEQVQZ9ZNHEREJtNQww9soQcRmgtqGXBuBB2Ic+3tm1pNTF7ERmAS25j+s/FYt6mbVou5nbHN37ts9yNdu38nXbt/Jd+97ks/84XP4zWOWVqtZAIxNBJmISgQR0SyYoxMZlvaV/fSzcne++vMdfOT6h+hsa+HPzjqW83/jaNYt7amrbhURkUbVKIWV1wFnmNmGaIOZrQNeGN5XyP8A7cAfZB3bBrweuMHda5rTNjNOXb2IT77uVK79sxeysLudt/zfX/Dd+56oajtGZgory58p6AszEaNVLK50d/7xew/yof95gBcft4wbLzqTi845gfXLehVAiIiUSaMEEV8CdgDXmtl5ZnYucC3wGHBFtJOZrTWzKTO7JNrm7ncTDO/8tJm9w8xeRjC8cz3w91W8hjmdfNRC/utPX8BzVi/mL795Dz/furdqjx3NKFmJtTMOZSKqF0T8y81b+defbedtL1jHl96yiRX9XVV7bBGRZtEQQYS7jwJnA48AXwe+AWwHznb3kaxdDWjl8Ot6O/Bl4DLge8Bq4BXufleFm160hd3tfOmtm1i/rJc//4+72TtSnUTJ6GSG9lajo638L4moi2S0SnNF/HzbXj75o0d4zWmruOTVG2lpUeZBRKQSGiKIAHD3Xe5+vrsvcPd+d3+Nu+/I2WeHu5u7fyhne8rd3+/uR7h7l7s/391vqWLzi7Kwu53PvfG5DE9M8f9d8+uqPObYxFRF6iHgUGFlNTIRY5NT/NW37mX90l7+8fdOUQAhIlJBDRNENJvjV/bzvpcfxw0PPM1tVejWGJ3MVGSiKTg0gVU1gogrbn2UJwfH+afXPrtiQZGIiAQURNSxP3rheo5a1M1Hv/8g09OVnehodGKqInNEQFZ3RoWDiD3D41zxk2286pQj2bSuanOIiYg0LQURdayrvZX3vfw47n98iFu3DFT0sUYnM/RU6Jt7FJxUuibiaz/fycTUNBedc0JFH0dERAIKIurceacdxYr+Tr5y246KPs7YxNTMUMxy62xroa3FKpqJGJuc4qpf7OS3N65k/bI4q8OLiMh8KYiocx1tLVxwxlpufWSAbQMjcx9QotHJTEWGd0IwF0allwO/7p4nODiW5h0v3jD3ziIiUhYKIhrAG05fQ2uLcc1duyv2GKMTUxWZaCrS19k2M6FVJVxz1+Mcs7yXTWvnWo9NRETKRUFEA1je38kLj13Gtfc8gVdoJcmxyamK1UQA9HS2zUxoVW6P7R/jlzv28/vPPVqzUYqIVJGCiAZx3qmr2H0gxV275lpvrDSjE5mKZiJ6O9sYmahMd8a19zwOwHmnrarI+UVEJD8FEQ3it09aSUdbC9+776mynzsz7aTSmYrOq9Db0cpYhbozbnjgaZ6zZhFHL+6pyPlFRCQ/BRENor+rnd/csJSbH95T9nOn0uEKnhUqrIRg/YxK1ETsGRrnvt2DvPzElWU/t4iIFKYgooGc/awVbN87yva9o2U9bzT0sqdCQzwhWMmzEqMzbnooCKpeduKKsp9bREQKUxDRQF56QvBBefND5c1GjM4sA954hZU3PrSHoxZ1c8LK/rKfW0REClMQ0UDWLO1hw/JebnmkvLNXRhmCStZEVGKIZ2bauePRfbzk+GUalSEiUgMKIhrMC45Zyq927GcqM122cx7KRFSuO6Ono5Xx9DSZMq4B8sATQwyPT3HGhqVlO6eIiMSnIKLBPH/9UkYnM2x+Yqhs54wyEZWcJ6IvWoSrjF0adzy6D4DfVBAhIlITCiIazPM3BKtT/mL7vrKdc7gqmYggiBgr41wRtz+6jw3Le1mxoKts5xQRkfgURDSYFf1dbFjWyy8e3V+2c46MB0FEf1d72c6Zqzcc+VGuuojMtPO/2/fz/PXKQoiI1IqCiAZ0+vol/O+O/WWbAntkIg1Af1clJ5sKMxFl6s7YNjDC8MSU1soQEakhBREN6DlrFjE0PlW2+SJGxqcwC4ofKyUa+RFlPebr7nD679PWLCrL+UREpHgKIhrQqauDD857dx8sy/mGJ6bo62yr6DDJKMsxXKbujHseO8jC7nbWL+0ty/lERKR4CiIa0HEr+unpaOXexwbLcr6R8Sn6KzgyAw6NzihfJuIgp65eREuL5ocQEakVBRENqLXFOPmohdzzWJkyEeNT9FWwHgKyMhHj6Xmfa3RiikeeHua01erKEBGpJQURDeo5qxfxwBNDTEzNf8jkSNidUUlRkFKO0Rmbnxhi2uG01QvnfS4RESmdgogGdfJRC5nMTLPl6ZF5n2t4Yoq+Cg7vBOhsa6WjraUsNREPPBF045y0SkGEiEgtKYhoUBtXLQDgwSfnP3PlyHi64jURAAu62hguQ03E5ieGWNrbwYr+zjK0SkRESqUgokGtW9pLd3srD5QjiKhCdwaEi3CVKYjYuGqBFt0SEakxBRENqrXFeNaR/TxQhjU0RsanKjrRVKS/q33ehZWTU9Ns2TM8k4kREZHaURDRwDYeuYAHnxya18yVmWlndDJT8dEZUJ7lwLfsGSadcdVDiIjUgYYIIsysxcwuNrMdZjZuZvea2fkxjltgZpeY2c/NbJ+ZHQx/f0012l1pJx65gKHxKR4/mCr5HNGHejW6M/rLUBPx0JPDAGw8sr8cTRIRkXloiCAC+DDwIeBzwCuBO4Bvm9nvzHHcGuDdwK3ABcDrgUeA75jZn1WstVVy4pFRceVwyeeIgohqdGf0lSGI2DYwQluLsVYzVYqI1FzlPznmycxWABcBl7v7J8LNN5vZscDlwPUFDt8ObHD3saxtPzSz1cDfAP9SiTZXy3Er+wDYumeE39q4sqRzRIWOfZ2VHeIJsKAMNRHbBkZYu7SH9tZGiX9FRJKrEd6JzwE6gKtytl8FnGJm62c70N1HcwKIyJ3AqvI1sTYWdLWzckEnW/eUPldEtIJnNWsi5lPDsW1glGOW95WxVSIiUqpGCCJOAiaArTnbN4c/N5ZwzpcAD82nUfXiuBX9bN1TenfG8Hh1ayKmHcYmS5tlM52ZZue+UY5ZoSBCRKQeNEIQsQQ46Id/fd2fdX9sZnYhcAbw0TK0reaOXdHH1j0jJX+7j4KIatVEQOlTXz+2f4x0xtmwTPUQIiL1oOpBhJm93Mw8xu2WCjz2WcBngK+5+zcK7Hehmd1pZncODAyUuxlldeyKPkYnMzw5OF7S8dUdnRHUXZRaF7FtYBRAmQgRkTpRi8LKnwMnxtgvqmU4ACwyM8vJRkQZiP3EYGbPA64DbgLeUWhfd78SuBJg06ZNpXfgV8GxKw4VV65a1F308SNVzEREU2uXOkJj20BQ+3HMMgURIiL1oOpBRFjoWEw9wmagEziGZ9ZFRLUQD8x1AjM7BfghcA9wvrvPfz3qOnFcGERs2TPCS45fXvTx0YJYvR3VqYmA0oOIRwdGWNbXycKeyo8kERGRuTVCTcQPgDTwppztFwD3u/v2Qgeb2XHAj4BHgVe7e+kzM9WhpX2dLO5pL3mExlAqTX9XGy0tlV+HYr41EcHIDNVDiIjUi7qfJ8Ld95jZp4CLzWwYuItg0qizgXOz9zWzG4G17n5s+PcKggCiA/h7YGPOok13u/tE5a+isuYzQuPg2CQLu6vzzX4+NRHuztY9I7zq2UeWu1kiIlKiug8iQh8ERoD3AkcADwOvc/fv5uzXyjOvaSOwNvw9d1+A9cCOsra0Bo5Z0cf3738Sdy96ZcvBVLqKQUTp3Rn7RycZTKU1R4SISB1piCDC3TPAZeGt0H5n5fx9C5D49aKPW9HHf4yl2Tc6ybK+zqKOHUylWVSlGoO+jjZaLHjMYs2MzFB3hohI3WiEmgiZQ/YIjWJVMxPR0mIs7G7n4FgpQUQ4MkOZCBGRuqEgIgGiNTS2lBRETFUtiABY1NPBwVIyEXtG6GxrKWkYq4iIVIaCiAQ4YkEXfZ1tbCsyiHB3BlOTLKhiEBFkIiaLPu7RvaOsX9ZLaxVGkYiISDwKIhLAzDhmeW/R3RmpdIZ0xlnU3VGhlh1uUU97iTURI5qpUkSkziiISIgNy/vYvne0qGOiD/Nqdmcs7C4+iBhPZ3hs/5jqIURE6oyCiIRYv6yXxw+mSBWxQmYtgohFJRRW7tw3xrRrZIaISL1REJEQG8IP2GKyEYNjNchE9HQwNJ4mMx1/SRKNzBARqU8KIhJiQ7goVTFBRDRKolrzRECQiXAvbtbKqGB0gzIRIiJ1RUFEQqxb1gMEi1TFVZPujDBgKaZLY9vACKsWdtFThUXCREQkPgURCdHT0caqhV08WkQmYigMIqo5xHMmiCiiuPLRvaMamSEiUocURCTIhuV9RQURg6k0ZtDfWb1v+AvD4aRx54pwd7btGVE9hIhIHVIQkSDrl/Xy6MAI7vGKFqMpr6uxDHgk6jqJO8zz6aEJRiczGpkhIlKHFEQkyIblvQyPT7F3JN63/INj1Vs3IxJ1Z8QNIjQyQ0SkfimISJANy4sboXEwlWZRlYOIKGiJW1g5E0SoJkJEpO4oiEiQDcuClH/cERr7RiZYWuTS4fPV3tpCX2db/CBizwh9nW2s6K9uO0VEZG4KIhJk1aJuOtpaYhdX7h2ZYFlf9dbNiBSzCNe2gVGOWd6LmRbeEhGpNwoiEqS1xVi3tIdHB+YOItydfSOTVc9EACzr62DvaLwg4tGBkZluGhERqS8KIhJmw7I+Ht07d3fGYCrN1LSzrCZBRCcDwxNz7jc6McUTg+MamSEiUqcURCTMhuW97No3RjozXXC/aARHLbozlvd3sndk7iAiKhDVyAwRkfqkICJh1i/rZWra2X0gVXC/feGH+NLe2mQi9o1MzLkIl0ZmiIjUNwURCRPVD8w1QmMmE9Ffm0zEtMOBOYort+0ZocVg7dKeKrVMRESKoSAiYaL6gbmKK/eN1jYTAcxZF7FtYJQ1S3robGutRrNERKRICiISZlFPB4t72ucc5rl3eAIzWNJbm0wEMGddxLYBrZkhIlLPFEQk0IblfXN3Z4xOsqSng9YqrpsRiYKIQpmIzLSzXat3iojUNQURCbRhWW+sTEQthnfCoREhhTIRTxxMMTE1PTMLp4iI1B8FEQm0fnkvA8MTDI/PPrX0vtFJltZgeCdAX2cbXe0tBTMRW8NMiiaaEhGpXwoiEmjDsmiExuzZiGDK69pkIsyMZX2dBVcb3fL0MADHr1QQISJSrxREJNBx4Qfvlj356yLcnacGx1m5oHaLWi3vLzxr5cNPjbCiv5NFPbXJloiIyNwURCTQ2iU9dLS18PBTQ3nvHxiZYGJqmtVLajf/wvK+TvYMj896/5Y9wxy/sr+KLRIRkWI1RBBhZi1mdrGZ7TCzcTO718zOL+E8G8xszMzczI6tRFvrQVtrC8et6OPhp/NnIqLZLI9e3F3NZj3DUYu72X0ghfvhs1ZOTztbnh5RECEiUucaIogAPgx8CPgc8ErgDuDbZvY7RZ7n88BgeZtWn05Y2T9rJuKx/WMAHL24dpmItUt6GJvM5K2L2H0gRSqdUT2EiEidq/sgwsxWABcBl7v7J9z9Znd/F3AzcHkR53kj8BzgY5VpaX054Yh+nh6a4GCeqaXrIROxJpzKelcY0GR7JCqqPEKZCBGRelb3QQRwDtABXJWz/SrgFDNbP9cJzGwx8CmCYORg2VtYh04IP4Affmr4sPt2H0ixtLeDno62ajdrxpqwHuOxPEHEw2EQcZwmmhIRqWuNEEScBEwAW3O2bw5/boxxjn8CHnL3r5ezYfUsCiIeyhtEjNU0CwGHulJ27js8iPj17kHWLu2hv6u92s0SEZEi1O6raHxLgIN+eAXe/qz7Z2VmLwbeQtCVEYuZXQhcCLBmzZr4La0jRyzoYllfB/fuPjzxsvtAio1HLqhBqw7pam/liAVdebsz7tt9kN9YV/BpFRGROlD1TISZvTwcHTHX7ZYyPFYHcAXwz+7+QNzj3P1Kd9/k7puWL18+32bUhJlx2upF3PPYM4OI6Wnn8QMpjl5S20wEBHURu/Y/c0KsgeEJnhgc59SjF9aoVSIiElctMhE/B06MsV/0FfUAsMjMLCcbEX1V3c/s3gcsBj5jZovCbdGQhH4z63f3w/P9CfGcNYv58YN7GBxLs7An6Bp4/GCKycw0q2s4MiOyZkkPP90y8Ixt94WZk2cfvSjfISIiUkeqHkS4+xjwUBGHbAY6gWN4Zl1EVAtRKMOwETgCeDzPfXcB9wKnFdGWhnLa6uCD+J7dBznz+CCjcu/Mh3Ttv+mvX9bL1b/azdB4mgVh/cO9uwdpMTj5qNp2t4iIyNwaobDyB0AaeFPO9guA+919e4FjLwdemnOLhnheALyjvE2tL88+eiFmcM+uQ10a9+w6SGdbC886ovYf0lGQc9fOAzPbfrVzP8ev7K/pyBEREYmn7oMId99DMDzzYjN7v5mdZWZfAM4GLs7e18xuNLOtWcc+5O63ZN84lAX5hbvfWaXLqIn+rnaedcQCfrb1UJfBPY8d5OSjFtLRVvun/jlrFtHWYvxye9AjNTye5pfb989kTUREpL7V/pMkng8ClwHvBX4IvBB4nbt/N2e/VhpjxEnVvOKkI7hz5wGeHhonnZnm148PzmQAaq2no42Tj1rI/+4IgoifPLKXdMZ5+caVNW6ZiIjE0RBBhLtn3P0yd1/r7p3u/mx3vzrPfme5+7o5zvUVdzd3z513IpFe9ewjcIcf3P8U9z8+yMTUNM9ZUx9BBMDp65dw72ODjKcz/PjBp1nc085z1yyudbNERCSGhggipHTHrujnuBV9/Mcvd/GJGx6mv7ONFxyzrNbNmvGbG5YymZnmI9c/yPfvf5KXnbiS1hardbNERCQGBRFN4KJzTmDrnhFu27qPv/rt41nS21HrJs048/jlnHPSSr52+04WdLXz1+ecUOsmiYhITKofaALnnHQEX3778/jJIwNccMbaWjfnGVpajE++7jSWXf8gr3/ealYs6Kp1k0REJCY7fDZpybZp0ya/885ED+IQERGZYWa/cvdNcfZVd4aIiIiUREGEiIiIlERBhIiIiJREQYSIiIiUREGEiIiIlERBhIiIiJREQYSIiIiUREGEiIiIlERBhIiIiJREQYSIiIiUREGEiIiIlERBhIiIiJREQYSIiIiURKt4zsHMBoCdZTzlMmBvGc9XS7qW+pOU64DkXEtSrgN0LfWoEtex1t2Xx9lRQUSVmdmdcZdYrXe6lvqTlOuA5FxLUq4DdC31qNbXoe4MERERKYmCCBERESmJgojqu7LWDSgjXUv9Scp1QHKuJSnXAbqWelTT61BNhIiIiJREmQgREREpiYKIKjCz1WZ2tZkNmtmQmV1jZmtq3a5CzOwsM/M8t4M5+y02s381s71mNmpmPzazU2rY7qPN7LNmdruZjYVtXpdnvy4z+7iZPWlmqXD/l+TZr8XMLjazHWY2bmb3mtn5dXYt+Z4nN7PT6uFazOy1ZvZfZrYz/Ld+2Mw+amb9OfvFei3Ffe5qcR1mtq7A87GoHq4jfOxzzOwmM3vKzCbMbLeZfcvMNubsF+u9q5bvA3GuJe77Wa2vJU9bfhC287JS2liV15i761bBG9ADbAHuB14DnAf8GtgG9Na6fQXafRbgwJ8DZ2TdNmXtY8DPgN3AG4BXALcSjFk+uobtfhq4HvhheA3r8uz3DeAg8E7gZcA1QAo4LWe/fwQmgIuAlwJXANPA79TRtTjw5Zzn6Qygpx6uBbgD+BbwJuBM4H3hv/0dQEuxr6W4z12NrmNd+Hx8JM/z0VoP1xE+9huAjwOvDa/lzcBmYIhgjgCI+d5V6/eBmNdyFnO8n9XDteS5rifDdl9WShur8Rqr2j9Is96A9wIZ4NisbeuBKeD9tW5fgXZH/+leXmCf88J9Xpq1bSGwH/hMjdrdkvX7O8jzwQucGm5/e9a2NuBh4LqsbSsIPnT/Ief4G4H76uFawvue8SYzy7lqdi3A8jzb3hK2++xiXktxn7saXse68O93zHGuml1HgTadELbpr8K/Y7131en7QO61zPl+Vk/XAiwGniIIEnKDiLr6v6LujMo7F7jD3bdGG9x9O3AbwYuhkZ0LPOHuN0cb3H0Q+B9qdG3uPh1jt3OBNPDNrOOmgP8EzjGzznDzOUAHcFXO8VcBp5jZ+vm3eHYxryWuml2Luw/k2fy/4c+jwp9xX0txn7uyi3kdcdXsOgrYF/6cCn/Gfe+qu/cBDr+WuOrlWj4G3O/u/5Hnvrr6v6IgovJOIkgH5toMbMyzvd58w8wyZrbPzP49pz+00LWtMbO+6jSxaCcB2919LGf7ZoIP2mOz9psAtubZD+rr+fvTsD94LOwffnHO/fV2LWeGPx8Mf8Z9LcV97qol9zoiHzWzqbCW4Lo8/dV1cR1m1mpmHWZ2HEH31lNA9MEV972rLt4H5riWSKH3M6iDazGzFxFkuP5sll3q6v+KgojKWwIcyLN9P0HKql4NAp8kSKOfDXwYeDlwu5mtCPcpdG1Qv9c3V7uXZP086GEesMB+tXYV8G6C5+dCYClwk5mdlbVP3VyLmR0FXAr82N3vzHr8OK+luM9dxc1yHRMEH2DvIqg7uQg4Bfi5mZ2YdXi9XMcvCNr8CPBsgm6ZPVltiPPeVS/vA4WuJc77GdT4Wsysg+D18wl3f3iW3erq/0pbOU4iyePudwN3Z2261cx+AvwS+Avgb2vSMDmMu78568+fmtm1BN9ULgNeVJtW5Rd+S7qWIM389ho3p2SzXYe7Pwn8SdauPzWzHxB8+/sgcEE12xnDm4EFwAaCgOdHZvYid99R01aVZtZraaD3s78GugmKoBuCMhGVd4D80etsUWLdcve7CKL854WbCl1bdH89mqvd+7P2W2RmNsd+dcXdh4Hvceh5gjq4FjPrJui33QCc4+67c9oX57UU97mrmDmu4zDu/hhBNX3u81HT6wBw9wfd/Rdh3/vLgD7gA+Hdcd+76uJ9YI5rybd/7vsZ1PBawq6VDwJ/B3Sa2aKsYcHR361FtLEqrzEFEZW3maBvKtdG4IEqt6VcopR4oWvb5e4j1WtSUTYD682sJ2f7RmCSQ3UDm4FO4Jg8+0H9P3/ZXRc1vRYzaweuBjYRDCn9dc4ucV9LcZ+7iohxHYXkPh81u4583P1g+LhRX3nc9666ex/Icy0Fd8/6vZbXsgHoIuiePJB1gyCzcoCga6yu/q8oiKi864AzzGxDtMGCCYNeGN7XMMxsE8HQqV+Gm64DjjKzM7P2WQD8LvV9bf8DtAN/EG0wszbg9cAN7j4Rbv4BQXXzm3KOv4Cgcnp7FdpatPA5eDWHnieo4bWYWQvBePWzgde4+x15dov7Wor73JVdzOvId9wagm6l7OejZtcxGzNbCTyLYB4IiP/eVXfvA3muJd8+ue9nUNtruYegjib3BkFg8VKCD/76+r9SrXGvzXoDesMn/tcEw2/OBe4FHgX6at2+Au3+BkGf+u8TvGn+FcFkJruAZeE+LcDPgceAPyQYRngLQZpsdQ3b/trw9gWCbxl/Gv59ZtY+/0kQ2b+DIPV5NTAOPDfnXJeH299PMNb8CwQTNL26Hq6F4BvKl4A3hu17a/hamwReXA/XktX2yzh8Aqaji30txX3uanQdnwT+GXgdwZv+nwA7CSb8OaEeriN87O8QpM3PC9v5LuChsJ3Hh/vEeu8q5rmr4bXM+X5WD9cyy/XlzhNRV/9Xqv4P0ow3YA3wXwQzqA0D/02eSYPq6QZcDNxHUNWcDl+wVwJH5uy3BPi/4Qt4jGDyolNr3Haf5XZL1j7dwKcIhoGNE1R2n5XnXK0ERVc7CSq/7wNeWy/XQvDt47bwDTFNMD7+OuD0erkWYEeB6/hQsa+luM9dLa4D+COCuSMOhM/HU8C/kxNA1PI6wsf+G+BXBB+0YwQTEF2R+74U972rlu8Dca6FmO9ntb6WWa7vGUFEMW2sxmtMq3iKiIhISVQTISIiIiVRECEiIiIlURAhIiIiJVEQISIiIiVRECEiIiIlURAhIiIiJVEQIdKEzMxj3HaE+34l+r1emNlnzOy7VXy815jZ03W8vL1ITWieCJEmZGZn5Gz6DsFshB/K2jbh7neb2THAAg9WQqy5sD0PAi/wQ0twV/oxjWAVyGvd/e+r8ZgijUBBhIgQZhp+5u71tlT1Yczss8AZ7v68OXcu7+O+G/gwcJS7j1fzsUXqlbozRKSg3O4MM1sXdnf8iZl91MyeMrNhM7vKzHrM7Fgz+6GZjZjZVjN7a55znmpm15nZATNLmdltZvbiGG3pJFg07N9ztp8Vtuk1ZnaFme03s4Nm9mkzazWz55nZz8xs1Mw2m9k5Occ/z8x+ZGb7wvY8amafz3n4bwGLCNZfEBEURIhI6S4GVhEs+nUJweqAXyToGvke8HsE6xV82cxmli42s+cSLCC0BHgncD7Bmh8/NrPfmOMxzyD4IP/pLPd/GhgN2/JZ4L3htq8RrDXw+wTrDVxjZsvC9vQBPwQywNuAVwKXAm3ZJ3b3vQTdKK+Yo40iTaNt7l1ERPLa5u5RluGHYSbhzcCb3f0qADO7k2D1x9cCm8N9P06weuLZ7j4Z7vdD4H6C1RhfU+AxzyBYkOi+We6/yd3fH/7+IzN7FfAeghVNfxY+1pME9R+vAr5KsGT0YuCv3T37vF/Jc/67wzaICMpEiEjpvp/z90Phzx9GG9z9ALAHWA1gZt3AmcC3gWkzazOzNsCAHwMvmeMxVwFDUfARs02jUQCR087V4c8tBCtAXmFmF5jZamY3ELZBRFAQISKlO5Dz92SB7V3h70sIliT/O4IlmbNv7wEWm1mh96UugmXMi2nTwewNWQFIV/j3IPBS4Ang88AuM7vfzM7Pc/5U1rWIND11Z4hINR0EpoF/IahTOIy7Txc4fh9BTURZufs9wPlhVmQTQb3Ht8zsVHe/P2vXJWEbRAQFESJSRe4+amY/BU4F7pojYLLbGEoAAAEpSURBVMjnIaDDzI52990VaN8UcIeZ/R1BLceJBLUakfXAw+V+XJFGpSBCRKrt/cBPCIox/w14ElgGPBdodfcPFDj2J+HP04GyBBFm9mrgQuC/ge1AL/AXwDBwe9Z+Fj5u7tBPkaalmggRqSp3vwt4HkG3wGeAG4D/HziFQ0HCbMfuAH4J/G4Zm7SFoNbh7wgKM78MTAG/lZPteAHBKI7/LONjizQ0zVgpIg3FzN5GEHQc6e5jVXzcLwAnu/uck2KJNAtlIkSk0VxFMJLi3dV6QDM7gmBSrQ9W6zFFGoGCCBFpKGHx49uBqmUhgHXAX7l7we4WkWaj7gwREREpiTIRIiIiUhIFESIiIlISBREiIiJSEgURIiIiUhIFESIiIlISBREiIiJSkv8HlpgLNwGi9NEAAAAASUVORK5CYII=\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": [
    "Add our free surface boundary conditions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def freesurface_top(p_func, v_func):\n",
    "    time = p_func.grid.stepping_dim\n",
    "    pos = int(max(p_func.space_order, v_func.space_order)/2)\n",
    "    \n",
    "    bc_p = [Eq(p[time+1, x, pos], 0.)]\n",
    "    bc_v = [Eq(v[1][time+1, x, i], -v[1][time+1, x, 2*pos-i]) for i in range(pos)]\n",
    "    \n",
    "    return bc_p + bc_v\n",
    "\n",
    "bc = freesurface_top(p, v)"
   ]
  },
  {
   "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` run in 0.12 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "PerformanceSummary([(PerfKey(name='section0', rank=None),\n",
       "                     PerfEntry(time=0.10119199999999998, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
       "                    (PerfKey(name='section1', rank=None),\n",
       "                     PerfEntry(time=0.009614000000000003, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
       "                    (PerfKey(name='section2', rank=None),\n",
       "                     PerfEntry(time=0.005519000000000003, 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, eq_v_damp_left, eq_v_damp_base, eq_v_damp_right,\n",
    "               eq_p, eq_p_damp_left, eq_p_damp_base, eq_p_damp_right,\n",
    "               eq_p_i]\n",
    "               + src_term + bc)\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": {
    "scrolled": true,
    "tags": [
     "nbval-skip"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfUAAAGDCAYAAAAyM4nNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsnXuYHVWVt9+VNJ0mhJiEhISLISAKJoNcbBhUdBBRYUS8gIAKogiIIyoqjGDUYQA/YFRkRgcl4m3AEQfUAbwAAuEiEjFycQiCAomABkxMYkhC57q+P+qcrlW7z66urnO6u073ep+nn95VtXfVPnWqzqr6rb3XElXFcRzHcZz2Z8xwd8BxHMdxnNbgRt1xHMdxRghu1B3HcRxnhOBG3XEcx3FGCG7UHcdxHGeE4EbdcRzHcUYIo8qoi8gLReRaEfmbiKwWkR+KyMzh7pfjOI4ztIjIziLyZRG5R0TWiYiKyKyCbceIyDkiskREekTkQRE5KlL3FBF5RETWi8ijInJaKz9HyKgx6iIyHrgN2BM4ETgBeDEwX0S2Gc6+OY7jOEPO7sAxwErgrgG2PR84F/gKcDiwALhGRP7RVhKRU4DLgR8AhwHXAJeJyAeb6nkOMlqCz4jIR4FLgD1U9bHaul2BPwD/rKqXDGf/HMdxnKFDRMao6pZa+WTg68Cuqrqkn3bbA08BF6nqv5j1twLTVPVlteUO4M/Az1T1RFPvm8CRwA6qurG1n2oUvamTnMQFdYMOoKqLgbuBtwxbrxzHcZwhp27QS/BGoBO4Klh/FbBX7WUR4BXAtAb1rgS2Aw4qefxcRpNRnwM81GD9ImD2EPfFcRzHaU/mAOuBx4L1i2r/Z5t60NfuhPVaymgy6lNIfCchK4DJQ9wXx3Ecpz2ZAqzSvr7rFWa7/R/anbBeS+kYjJ2OFETkVOBUgG222eble+655zD3yHEcpzosWbKE5cuXS6v3u7uIrmtyH0uTN+Ies2qeqs5rcreVZzQZ9ZU0fiOPvcFTuwDmAXR3d+vCe+8dvN45juO0Gd0HHDAo+30eaHZ4+GehR1W7W9GfgJXAJBGR4G29/ua9wtSDxO4szanXUkaT/L6I1MdhmQ08PMR9cRzHcdqTRcA44EXB+rqP/GFTD/ranbBeSxlNRv164EAR2a2+ohZo4FW1bY7jOE5FGNPk3yByI7AReHew/njgodqsKoB7gOWReitIZl61nNEkv38dOB24TkQ+DShJAIGnSIIDOI7jOBVAGJo3ThE5ulZ8ee3/4SKyDFimqnfU6mwCvqOq7wdQ1b+IyCXAOSLyHHAfcCxwCMnUaWr1NorIZ0iCzfwJuKVW5yTgw6q6YTA+06gx6qq6VkQOAb5EMk9QgFuBM1R1zbB2znEcx8kwRDLyNcHyZbX/dwAH18pja3+WucAa4KPADOBR4BhV/bGtpKpfExEFPgGcBTwJnK6qlzFIjBqjDqCqTwIN4/M6juM4owtV7XfkfqM6qroZuKD211/7yxlCNXhUGXXHcRynPRhNA75aiRt1x3Ecp1IMlU99JOJG3XEcx6kcbtTL4UbdcRzHqRT+pl4eP2+O4ziOM0LwN/USbDHPQmMomL1v06a03NOT3dbV1Vtct6mzt7xqVbyJ3V2H+Ra32SYtT5uWbTNm+V/ShTVmFp/d2aRJ2UZTp/YWV6xKP3fYnzoTJmSXzUejs8Ocq9gHCNjiz51OQPSes9eULYeYi3L1mvT6svdbo+U69hbZeWfTryVPZCs+80zj/tgdzJqVabKaiQ2bx+43yN5j5nZlyiRzntYEs3btPWd3YMi79wr/7jWB3/nlcKPuOI7jVI6WZ4kZJbhRdxzHcSqF0Dfai1MMN+olsNJTKFFFZSmrn4VSmJHmekyKXSu/LV8e352VvK0cOH1a0Be7Q1veuDEt77BDto3ZeU/P+IbNrZIXyoRWDuwMpPlecqTSMWbnLsWPTnKlXnvB5cnvkevISuyPPZZtYq9xi1XMZ84wkT5/97tsxQceSMvr1zfeQeB6en672b3lJUvS9U8/3bgvADNmNN7dlA7zOxN+GPujYW/SHFfYUON3ezn8vDmO4zjOCKE6j2WO4ziOg09pawY36oNJbMR7qKUbJv1dKr/bJqF6lje4t5dQ5rda47JlaXmrrQrt2PYnpnqG2NH4EybYWQOGHKnUcfoQk9nzhoibayp2K4b32FNPpWV7i1i1OsPatdnlZ59Ny/bes30ORr9vOyuV320/7e26eXP2MHbwur3fMvp9OJTf+ukqihv1cvivp+M4jlM53KiXw8+b4ziO44wQ/E3dcRzHqRTuUy+PG/UWk4k2F/MNhz71pUvTNsaRtueeB/SW86Jd2dkpL3iBqZQzdY6tt07L48al5TC6lPkM9uOMNZNI82YS2Zk81kc4vsvsrNAAgfyphM7IIjONLbw+ikaOs5jrepW5/ex9FN4u1ndtb5fMbZ0Xmc3emHbn9kYw9z7A+J4VveWdd07H19ifjPAj77lnWp6+1kS1W7gwLYeRIiPT6qp0X1WnJ+2FG3XHcRynUggeUa4sbtQdx3GcyuER5crhRn0A1KWpUskMIkkTAHjwwbT8q1/1FqccnIaR+sfDDss0+fOqNLqb3fWUrnXpwjPBFB8rwVlpMBZdCtjSkSaYsSr9ttumZTvdJ8Qqjc89l5Y7OlJxrTNvClvBxC/OKKKI/B5kFlrdk17HVsoOJXeLvV222y4tR6e0hdPE5sxp3Mj22fqnAB55pLe43z779JZnzUrv91BJH/OLO9OFq3+clu2HO/jgbCNzfspI7lWS6Z0s/ivpOI7jVAofKFceN+qO4zhO5XCjXg436iUoKj1lRsLnSNyZIbhXX924fOihmSY7vupV6cKLX9x432FycysP2m1Gv9+wKfhskQHzVgKMpWYPl2Mj4TsmpNIowJhNG3BGH9ER73kXlcVcx1u6xmc2rTIJUezt9vzzaTm8Xezy9Olp2d5imXt8992zO7AVrea/PDL8HuChh9LyggW9xSm2ja0T1Mvs74gj0nIQuW7DhHRkPRX1cLlRL0eFvkLHcRzHcfm9Gfy8OY7jOM4Iwd/UBxGrEq5Zkz4/Tfq7l2Xqjdl7795yz9e/3lt+wAwX7/ze9zJtdjPLk/bfP91w9NFpOZDst+z+kt6yzc9spfC8+Bl2WyhVNtpXuGzPh5U9Q8mvTGAaS9EmVZIanRzyvlD7JZqL0iZAgXzFu044qtzmKd9pp7Q8eXLjrvUwMdP++XHp8ta7z+wtT9z5L2mlX/wie1ArpZt7/PfmRgpTq9u5OPuY8lQzen7DPgdgsclrKppO3d84S1Khr9BxHMdxXH5vBjfqjuM4TuXwiHLl8Ichx3EcxxkhVPZNXUSOBt4JdAPbA08CPwT+n6o+V6szC1gc2cVkVe31nolIF3A+cDwwCXgA+KSq3hlp31Ksb/nXv85u2/eUD/WWu269tbe85Uc/6i0vIMtvTXk/s8OXWWdi4CS0U256elJ/38qVaR07bQ2yPrZYOeZfh/jMJDu9LYxIVyTaXBjVr0yEq5irtkp+xVFF7GLJ86mbi2+1Gbfy179mq8USotjIiDZqHGT96NOmpWV77a3rSY8Z5mmKJV2aOMMMTlmyJNvI+NFvMz8aJjULJmYkEPjR99orXTj77N7i7bdn29jZrdEIecPMUISJFZEXAl8CXk8iDtwCnKGqT/bT7lzgXyKb16tql6m7BNilQb23qer/luh2LlX++TqTxJB/imRsyL7AucBrReSVqmp/1S8Erg/aPxcsfwN4E3AW8ATwIeAmEXmFqj7Q+u47juM4ZRgKn7qIjAduA9YDJwIKXADMF5GXqeranOZXADcG67aprQttEcBNJPbL8miJbvdLlY36m1XVjmG9Q0RWAN8BDib5Muo8oarhy2wvIrI38C7gJFX9Vm3dHcAi4DzgyBb33XEcx2mCIfANnwLsBuyhqo8BiMhvgT8AHwAuiTVU1acJJiKIyAkkNvU7DZosz7NRraSyRj0w6HXqOvNODbblcSSwEfi+2f8mEbkaOFtExqnq+mjrHIrKwFbye/zx7Laf/CQtn3zpD3vLr1x1SG953fz5mTaPmLKV4juMnDf7K1/JHsjoyi9561t7y3820aXyUrDbaWihTB/DStnWBZGJKBdchZnlrjHReoNFiVl0GUaCfN/sOYgRnptocqS8RD5meQONE7WE09bCqZZ17DRNO1UNstJ8JsqhOX5eN62sPXOGaf/jW9KyjRoJPGQ6+nuz3nb/JWR5q+3oA6no+KlPp/dOKLFb+T12veYlrhrshC5DNPr9SGBB3aADqOpiEbkbeAs5Rj3CicCzJG/lw0a7DZT7h9r/3wXrLxSRTSLyNxG5XkT2CrbPARarauiOWgR0AkF8R8dxHGeEMwd4qMH6RcDsgeyo5pt/LfBdVW30SPxmEVknIutFZIGIvLVBnZbQNu8UIrITiVR+i6rWx42sBy4HbgaWAXuS+OB/KSIHqGrd+E8BVtKXFWZ7o2OeCpwKMHPmzEZVHMdxnEGgBW+cU0XEjjGcp6rzzHKeXZjcYH0ex5N0uZH0fgOJyrwYmA6cDvxIRE5Q1asGeJx+aQujLiITgOtIUg+8r75eVZcCp5mqd4nIjSRPWnNJTnRpahfAPIDu7m6ty1F50pOVrOwo7rwR4j82KZAvvTQtz52bDhv45PHfzLQ58P3v7y3fZ9abQFF0/N//Zdq85Iwz0gUj0+1oEz/sY8fSwoZJ2/eWY9Hh8ojVi8mheViptIzEHfZlsGTywZKuRy3BF2WTtaw0TjoruYfye8x1tM02aTm8RzPRFXsaj8afaLo2sSM4qMmNzhW3p+Urr+wtPhmMfl9hyvZV8e2mvP1nPpNpc333eb3l04ysbvtvb33IToyxpzdPcrcUrdcMLTDqy1W1u/ndFOI9wP2q+ttwg6p+2C6LyI9IJjRdCLTcqFdefheRrUmedHYD3lgboBBFVZ8CfgGY2KmspPGTV/0NfUWDbY7jOM4wUPepN/NXgDy70OgNvnFfRQ4gUYkbvaX3QVU3A9cAO4vIDkWPU5RKG3UR2Qq4lmSu+j+q6v/108SiprwI2LU2hcEyG9gAPIbjOI5TGaTJvwIsIvGrh8wGHh5AV08kGYj93wNoU0f7rzIwKiu/i8gY4LvAIcARRacDiMhM4CDATuq/AfhX4B3UnqZEpAM4Frh5oCPfy0hUEyakz08velG2npXCnnsuDT5z9tmpKHE2R9sm/OAHJ/WW3/7Yv/WWez75yd5yeFUuNAliZnz5y73lnf/XnKq3ZsdvdJqkMJ1WmjfZLtZtSkcg58nqsXgi4Yh7y9hIBIow8UyRYw5kW52RMJK9bYgNJQ++bCutP/tsWi6StCXcnR04HiZ0GbNmdeOd27KV2BcEP1E/+1naxMjsJp1Lnx/g19gFM3vlhzukAaqOOmpj0OpHppz6E/ba6w295b/7u2wLOxp+KKT0inI98AUR2U1Vn4DegGavAs7OadeLiHQCxwE/i8zYatSmbnueVNVn+qs/UKr8k/WfJEb4c8BaETnQbHtaVZ8WkS+SqA33kAyU2wM4hyRx0efqlVX1fhH5PnBp7e1/MfBBYFfg3UPxYRzHcZxiCEMSUe7rJIPWrhORT5O8NZ8PPEUyADvpi8guwOPAeap6XrCPI0jk+obSu4i8k2R63E9r+51OEvhsP5KIqS2nykb98Nr/ubU/y7+SROdZRGKc3wtMAP5KEpTmX1U1jNbzPhJDfwFJmNgHgcNU9T4cx3GcSjHYvmFVXSsih5CEib2S5FniVpIwsVY/rD9jNOrSiSRjsn7cYBskL5DbA58nMf5rSaL+HqaqgzKfvbJGXVVnFajzTeCb/dWr1X0e+Hjtz3Ecx6kwQzHgqxbj/ah+6iwh4qZX1bf003YBiQt5yKisUR9pdJJGlNprr87MNuvGnj/fxsG5NlKGo45KvRFdXamQ8fif/rm3vN/J/5hp84Tx8dmRgU8+9VRveabxtUPgb7dT3w4+uLc43vjax++ejePz7LL01oxFp9u8mULYennT02y56PSymO+8zPS0keCHb+W0vMLnwzi77bTRcMyFdWk/80zj9eHYDusvf8EL0rKNrDZ+1Z+zjay//CETo+QXv0jLP/95b/HpwJFv92ZP5yxT3vGEEzJtfnvmf/WW997bfohPmXI4XtgOoP5Eb+m449K1B1rnJdC5ycbhitw8Tlvi36DjOI5TKYYoTOyIxI264ziOUzncqJfDjfpQYfTM8WuyMt3JJ6dR25YsSdPufulLx5paXwx2mE6f6el5c295p51O6S2fdtpPMy2+etPNveXd3pK6ghYarTKcBpeR5r/61d7yztcad0BElgeYbnS/ybunqSieNiGElgUTQax0utHM3rGScBj9y6qGsWlwITaPe15ijv7Wh4yWiHJlPmdft0njn24ruYd5yq3kvtKECLHXjY0UB9nc6Pvum5Y7H7g3XQiTjt+SJl7ZYmT2J0wVK7GHszltDpUD7Ny5G9OMne/6yiszbb639/fMUtFpz+m9/L737dFbPv30tEYf14Kd15cX7nKY8Df18rhRdxzHcSqHG/Vy+HlzHMdxnBGCqLY8St2IpLu7Wxfee2//FWNYbfCZIIiQ1Rf33LO3ePFXJ/aWzz57cbDD/zDlWJTbPYLldGT8ddelkv+Rj6QR6daYiHSQzdVunQY23q4ZQNw3h+3b3paWDzssLRuZfsOsbIZoK83Hos2FUngswlyZwbxFpXwr3xfZV9UpIqdvDIOZDZC8c2b3vXZtWg7dM3/9a1o2QRIz0eHMbQTAbh1Ppgt2Rsf3Url7dRARzt5V9o41mdEzEvt+2UMy/uKLe8s/3D299446aqmpdVHQ6gkaY0e4Z/KDcNZZaabpf7vIRIeznycMl2eH/Tchv3cfcAALFy4sGJW1OHNE9PtN7mMv+M0QJnSpDG30k+M4juOMBoYootyIxI264ziOUzncN1wOl98L0rT8brXNUH63QS5sfmUTyOXhGdmgRKeZLPJ33XWr2WJHzNrUESGv6y11daXJlu++O1trv2vToBerL7ywt2xleZu3NhtWB2aYspXmJ9hoGFaiBzjooLRsMlGsJnVHhLJ8LFlM0YQuZUa8NztKPkbR0ft5FA3oEyN23poN5hOut5/V9tl+v3lJgmxikh3tWPQfB1E7v/Wt3uJyI0vb2R5hZg37Ue11bGX2Seec01v+1Vv+X6b9q1+dljdutLNXbieOzQSazn7Zd9839ZZNnhcAXtllIl3bIDlWcg8zupiETIWzIzVgsOT3vUT0h03u4yWjVH73hyHHcRzHGSG4/O44juNUDn/jLIcbdcdxHKdyuFEvhxv1ocI6E8MpJHbZOhPN1JvZU00SCeDOS9Mobtc/nfrHTzstLS9d+j2yXGfKqR++p+f23vLLX/4mLJMnp37Cny1Iywfd0puunp5Pf7q3HOaxXRIpTzV+zd2CqUQ7TpuWLhx6aG9x4qtelZZf+tLsgWbNSsvGX7jOTL4L8m1kfLUxv22zPuRWt2mWVkSBixFLqjNuXFq2084AxneZKVj2S+gy5fCLe8xMNrvW3BfmfvnL/2WTnljfuZkxiTl6xm8OcIApW9/5fUen94G5PFl54eVkiWXjtD7sbJKvadOO7y1/7Wvp+rfvaT6BiUgHZMfo2Klq9p7IC8FYQTyiXHn8vDmO4zjOCKHaj2uO4zjOqMTfOMvhU9oK0vSUNkuo9Vr57IEH0vLChWn5sSBqnJXPTD5z3vve3uIlV6VR4wA+8Ql73G+a8s2mHM6FsvOr3mDKH+gt/eAHaZ23b/qfbHOTVeJhExrMxPfKROiCbLQ6G7FrR1PePtRxu83MFXs+7FSeINd7TLJf3ZNOzLN539evzzYvM42uXYhJ6UW9SJ1rzERHe30vDiIj2imcf/xj4/V2yiewxUjr9jqysnqQAyaDldntfKeOf//3TL2b9/xIb9nmLNq48TJT62c5R7L3zutMOU269MUvZqeTffx4Mw31299Oy3aqWoi9jmPXvpXloalpbJbBmtL2MhH9af/VcnnhKJ3S5m/qjuM4TqXwiHLlcaPuOI7jVA6X38vh8ntBWiq/h9gR7zabiZXcrCwPWXnSjg62IbaC3OZb3ntSb/nss9P1n/+8yZDRJ4fz7aYcCuV1THQ4js9s+djH0vzwl3zaSLJnntlb7DHRvgB+b8pWRrXOg/Bp1Er2U0zZOiCmhpLj/vun5ZhMH5Hokx2ac22jdxkdekuX7VlWmrdemDzJvkgSlbxEKVYyt6cgHAA9pmddumCvSXt9WSk9jIz4hz+k5UWL0rK9joNR6fb7XW3KNmBgeNXZU2V/+K03YGbQZsLhh6cLX/hCb/HiG2b3ls8+294HkHVR3Ulj7PvkwcG2E3pLZ521XW/5ggvSGp1X2WMAvzCj+W2iJ+vbCN1I9tq1mWzykra0aPT7YMnve4vozf1Xy2WGy++O4ziOM/z4lLbyuFF3HMdxKocb9XK4/F6QQZXfrfZqZc+8UcOPP56W7ehg2ybUdK0cZ6T51Ycd01s+99xsky99yYqiNsPx7aZs64TsZsrp6Plp09IgN0YNBeA9B5mc0heZfNNmNPDvA03apPLI9MaegfAJdqIp22zTUyLrJ4Yad0yat2Ur0YfLVhK1unjoJogNP88bch+LrGOvr+XBGHG7bK+jp55Ky8btE4SEyUjptmzr5eRmycjn9qyFUvqUF784XTj66LRsZn5c/8hLMm3sdX3//TbY0U9MOXBxZbBXy2tMOb13PvYxm4wFTEwmptxiZoXcfntatu42yH6/9n61kvuLXpRts+uuadlee/b6GqRgM4Mlv+8rorc1uY8pLr87juM4TjXwN/Vy+HlzHMdxnBGCy+8FGVT53RKT4sP411YetRKeHYFsg3mEbSwRWR5gw1tTedHK5FbO3LjxetNifrDzIGhOLzbzepDr2cib226bBu2wI/bN4Plkb1/7j3TBjKZfZWYN2EAlkJWIzbhvYgPRwyfgzki5K1IOl22bjki50XEbsSVY3hQp25HkoRTeEynH2oTnyfYhNio9jK9ux3FnRqh/6EO9xYd3zeYj+PKX0/I3vpGWN268y9QK79VHTdkEeMG6VKyrCOx1uNVWaR+srG6vw/E3BhnA82T2OqF7Zpd0tgjWzWBdPXkunbzpDYPAYMnv+4noHU3uY6LL747jOI5TDVxGLocbdcdxHKdyuFEvR2XPm4gcLCLa4G9VUG+yiFwhIstFZK2I3CIiezXYX5eIfF5ElorI8yJyj4i8JqznOI7jjA5E5IUicq2I/E1EVovID0UknHARa9vIPqmI7BPUGyMi54jIEhHpEZEHReSowflEFfapi8jBJA7ajwC/Nps2qerCWh0B7gJmAWcBK4FzgDnAPqr6tNnfd4E31eo9AXwIOBx4harmzWUBhtCnbsmbslQk4peNOhcuxyLShdgpMjZpiplK9P2fpdN9gpwY3HOP9anfasq/MeWl8eNnsD7P/YNtB/eWXvvadIzAB9K8Mxz7Zus5B669Ni3/2OS+viv1za425zNMEmLPmo2AFvNHQ9YPHfrB+1s/EGJP63Z9KNPFxgXYmHh2Ylfg2WWKTbJz0EFp+U3GJ35KmswE4LIr0qPaMRuLF9tbPozmZrOj/4XGbBMs22vHPvO/Oi29Opv0xLj1OfZwM1HS5G3PJF2K+c0h6/eO5TwPl+29F4lYCMSnPA4Bg+VTf7mI/rLJfXT141MXkfHAg8B64NOAAheQXPIvU9W1efsXEQW+DVwebPqtqq4z9T4HnAnMJfnhO44kq88Rqk3nrelDO8jvv1PVBZFtRwKvAg5R1fkAInIPsBj4Z5IHAkRkb+BdwEmq+q3aujuARcB5tf04juM4FWEIZORTSJ729lDVxwBE5LfAH0jSUF5SYB9/yrFPiMj2JAb9IlWtP7bOF5HdgYuAlhv1ysrvBTkS+HPdoAOo6t+AG4C3BPU2YqKnqOom4GrgjSIybmi66ziO4/RHPUxsM38FOBJYUDfoAKq6GLibrP1ohjeSiGBXBeuvAvYSkV37NmmOdnhT/66ITCVRO28CzlbV+uykOUCjRMOLgPeIyARVXVOrt9hKIqZeJ8nMmkVUjTxZLZalw0pz4dSXmOxnpfhQQrRyvpWoTeKJY00SiWO/cmim+ZNT9+stf/vb6QSmK65I6zz11IPZY3K3KVvPyBORMthod/Pnb2PK6bSg4/pMnfv73tLkye/pLR9qPoItB7P92G+SkX5/beRim8zksWBKnz3XJr88fzXJRJ57LtumJy8OW40wCp2VwrdLk4kwbVpaDqVfG7Vszpy0vO++vcUnN6VZ7W80uUcgO4PrxhvT8lM/M9EQTzeR1YCsZy38TotgJXObWOh1mVr775/WO97kHDr22LQ8/Y+Be80mVznTRG20kffsfRgm/LHL9lzb+zBsE5PZh3iqWhUYgjfOOcB1DdYvAt5RcB8fFJGzgM3AAuBfVNXOrZxDIu+Hc3vrPxKzSZTlllHlN/W/AV8ETgYOAc4HDgXuqUkakET0XNmgbT0d2OSC9aY02IaInCoiC0Vk4TL7A+w4juNUnan13+/a36nB9jy7MLnB+pCrgH8isUunAtsBt9XGg9ljrNK+g9dybU8zVPaRT1XvB+43q+4QkTtJokp8hGRgw2D3YR4wD5KBcoN9PMdxHCchGQfdBKrLBzP4jKqeYBbvEpHrSJTjC4CDGrcafCpr1BuhqveJyO9Jhz6vpPET1RSzvf5/l5x6KxpsqzZWgovJdOEoWSvtWdnPyoFLg5HoVi5+9tm0bCXIBWacSJD3faY5zmdN3ufP3n5wb/nhnr0zba6+Ol2+yniiFi+2UcGsbAtg83RbF8IDkTJYN9fKlWle7GuumWzKNiP7DkH7nSPld0Xb7LBDehx72rczavfk4IqOKa+x3OwAK837h1X2nzbOqqU/35xtlJmFYMu/iqw31wOQHYluX4ByBxFHsK6j/YJtqbS+115pPvTTTktrWIkdYOJDZiy1HbH+OaOKhhEX7Qm295L94mx5p52y7V/4wrRsXWF5I9lHoczeEJHmP3+Q9KkBefaj0Rt8Lqr6nIj8BHh/cIxJIiLB2/pAdsxNAAAgAElEQVSg2Z4qy+951E/OIhKfRchs4MmaP71eb9faFIaw3gbisUwdx3Gc4aCjo7m//smzHw83WF8Ua7wXAeOAILUe9afRZo7TkLYy6iLSDexBGtj5emAnEfkHU2ci8Obatjo3kAR5foep1wEcC9ysqusHueuO4zhOtbgeOFBEeoMYiMgskmnS10faRKnZniPIJh64kWTm1buD6scDD9VG27eUKgef+S7JqMD7SEa+70sSWGYdsJ+qLheRMcAvgBeSDT7zMmBvVX3K7O9qkukFZ9X2+0GSL+CVqnpff/0ZluAzrSam1+Yljonl2C5SDvdtj29lxp2zQT8wo+k5MB3R/OSEVGq1o6shOzD/5z9Pyz09VrK3ZUimo9axkr2VkfNyxZdhrCnbZCKx9DBQLqXLhkjZSpKh/D5QwnQ1sdA0thy6MKxXLP1+9903HbFvR6gDnGA8mTsuMbL6fJNMaFEwmSXMHV/Hyt9WFoesi8pK6zvs0LhOONvE7nuEjmQfrOAz3WPH6sJtwgBCA0Oee66/4DPbkASfeZ40+Mz5wLYkwWfW1OrtAjwOnKeq59XWnUnygjkf+DPJhVxf9zo7Al5ELgLOAD5FYs+OJZkHf6Sqml+u1lDlK+oh4J3Ah0ki/DwD/JBkysByAFXdIiJHAF8ALiP5lbkHeK016DXeB3yOZBDDJJIv87AiBt1xHMcZQlrhU+8HVV0rIocAXwKuJJkefytwhnHdUls/luzT9aPA22p/LyB5+r8beL+qhm9/c0mCTn6UJFHho8Axg2HQocJGXVUvBC4sUG8FcFLtL6/e88DHa3+O4zhOVRkCow5Qi3mSG4ddVZeQGHa77gYSt26RY2wmeZm8oFwvB0ZljbrjOI4zShkioz4S8bM2migzDS6WiML6KG1gnmeDaU52ipydW2Xbh1HXHjHRu265pbdop8edGkRDO/V044e/Ik2S9GTPHr3lhQv3sE0yM/Fs2QaH6+mxvvY/ZvuZmd5lPltmNkzok7fTu6yv285JC33dYVqYRowNlq1ffmJkfeiztPVMFDqmm7Id/5CdJTptWtrGzF7MlO1wiXDbPia31cSnzaDgBUFo7QvMlDR77djxIGOD82Ej6Vk/+C7mM9gpaGG9IlPSwqh+I8h37rQPfqU5juM41cMfhErhZ81xHMepFi6/l8bPmhOX5SEuzcei0+0aJB2ykmiR6XFhPdveJpsJE89Yidb0eabp50yb2AR4+98ZufWIWWnZJDZZPSGVm3/3u+zUu1h6+j/9Kd7NmNdirVHl7SxAyM4EtEGytjIz4sLfP/tV2ZlBVoUumu/H5nmx8rldDzBmiUnIYj+4/dDhh3vAfL+3R66VcDqavSbsB4/lIg+3xSLCFZ2SFku0VNAIbSkYHmRMn2mKowg36qXxs+Y4juNUCzfqpWmriHKO4ziO48TxRyEnnyJSo5Xow6hcsdHzVkLNi2IXqxe2sctW4rUysNXIAX5lEpVEXBATTfnvg8+WWZ5syruachgVy2ZriY2cDkdRF3ljsRo9ZCMGxqIHrgxyVsTO7wNmZP8tke8mPI4lJmOHy7HrKLymLOPGNa63QxC5zsrxtpzXtwLSelEp3Rkg/qZeGj9rjuM4TvVwo14KP2uO4zhOtfA39dL4WXPKUXTUr5VUrURspdJQtrUScUw6DkdRF5Hp8yT7mFyd194Syt8xYucq7xyWkd+LbCvaJlbOG3Jv3S7TTfCaUBa39UK3Q528vNh2CkCeZF8gYIxL6c5IwI264ziOUy38Tb00ftYcx3GcauFGvTR+1hzHcZxq4Ua9NH7WnMIU8Tn2iYJVxPce+lKtz7OIbxeK+d7z/PDPP5+WbXg3689dvz7b3rax22zfwvECRT/PQCnqk8/7Pop8V3nRB2NRBmOJUcJ92H3HxjiEy7GplXmREd2P3h64US+FX9GO4ziOM0LwRyHHcRynWrj8Xho/a05hYgkmrIRZRs4ckycDxyg6HauojGvLRSXy2HHKyO9507Y2h/nVGxDmD7fEMr+E5zm2j1j7PLeJldlNed0mm889ewrs7jonmBzyodskFrmuoGuhSpJ7eE9VqW/Dihv10vhZcxzHcaqFG/XS+FlzHMdxqoUb9dL4WXOGnVByLCLz09HZsA4Ecn6eXNzK0edlJPsi61tBmSh2ZfZr5PctEyb2lmOB+/LonFBwlH5sfRmXzjBQ9Np3nKJU80p3HMdxRi/+pl4aP2uO4zhO9XCjXgo/a07TWMmwFaN3m92HbR+V4qF4opL+1oeUSa5SNYq4JkJ3hlmO5eQJd1XolJZJltMmBiFv9PuoluL9Tb00Pn/CcRzHcUYI/ijkOI7jVAt/Uy+NnzXHcRynWrhRL42fNaelVC1CVmEfZQk/evSzFZ1uVzUG6u8PPsuGTen5KLorG8Qus7t2GXvQJD6lLYIb9dJU1qcuIreLiEb+bqzVmZVTZ1Kwvy4R+byILBWR50XkHhF5zfB8OsdxHCeXjo7m/gogIi8UkWtF5G8islpEfigiMwu06xaReSLyiIisE5EnReS7IrJrg7pLIjbqrSXOSr9U+VHon4CJwbpXAJcA1wfrL2yw7rlg+RvAm4CzgCeADwE3icgrVPWBlvTYcRzHaQtEZDxwG7AeOBFQ4AJgvoi8TFXX5jQ/DpgD/AewCNgJ+AywUET2UdWngvo3AecG6x5t+kM0oLJGXVUfDteJyCnABuDqYNMTqrogti8R2Rt4F3CSqn6rtu4Oki/jPODIVvXbcfKISfaVll1LyKA2P02ZGWmtnsY23G6gGJX+3oeToZHfTwF2A/ZQ1ceSw8pvgT8AHyB5gYxxsaousytE5G5gcW2/nw3qL8+zUa2kmld6A2pPVe8AblDVFQNsfiSwEfh+fYWqbiJ5OHijiIxrWUcdx3Gc5qgb9cGV348EFtQNOoCqLgbuBt6S1zA06LV1fwSWkby1DxttY9SBtwHbAt9psO1CEdlU84tcLyJ7BdvnAItVdV2wfhHQCeze+u46juM4pRgaoz4HeKjB+kXA7IF3WV4KbA/8rsHmN9d87+tFZMFg+dOhwvJ7A94D/AX4mVm3HrgcuJnkCWlP4FPAL0XkAFWtn9wpwMoG+1xhtvdBRE4FTgWYObPfsRNOAwolZ6kyRSLNtYDBPB9NS7yxcxCcj00FkrWEpzCWqp1mB7+3ychpH/0eoTXy+1QRWWiW56nqPLOcZxcmD+RAItIBfI3EDn0j2HwD8GsSaX46cDrwIxE5QVWvGshxitAWV76I7AgcCvx7TTYHQFWXAqeZqnfVRsYvAuYCxzdz3NoFMA+gu7tbm9mX4ziOM6QsV9XuITrWV4BXAm9S1cyDgqp+2C6LyI+ABSQDvFtu1NvkdYnjSfraSHrPUBt1+Atgf7N6JY2fvOpv6AP10TuO4ziDyeDL73l2odEbfENE5CISRfckVb25v/qquhm4BthZRHYoepyitMWbOsl0gwdV9cEBtLFv1ouAt4nI+MCvPptkNP1jOG3NEKnkQ8JokmCj31XRXPdt+GW3OgHSiGRoRr8vIvGrh8wG+sy+aoSIzAU+CXxYVa8s0YeWK8CVv6JEpJvkJPf7ll6rPxM4CLjXrL4B2Ipk9Hy9XgdwLHCzqq5vWYcdx3Gc5hiagXLXAweKyG7pYWUW8Cr6xj1p0EX5CMm89rmq+pXiH63X9jypqs8UbVeUdnjMfQ/JsJnvhhtE5IskDyb3kAxQ2AM4B9gCfK5eT1XvF5HvA5eKyFYkAxY+COwKvHuwP4DjOI5TOb5OMmjtOhH5NMlb8/nAUyQDsAEQkV2Ax4HzVPW82rrjgEuBG4HbRORAs9/V9TgrIvJOkulxP63tdzpJ4LP9gHcOxoeqtFGvGeB3Ajeq6l8aVFlEYpzfC0wA/koSIehfVTWM1vM+EkN/ATAJeBA4TFXvG5zeO3kMpgQ5kqT4ShA5ieH3ViRejI31HlLY7dBsvvthxiX3AgyB/K6qa0XkEOBLwJWAALcCZ6jqGtsbYCxZZfuw2vrDan+WO4CDa+XFJNPcPk/iq18LLCSxPTe18vPUqfRdoKobgWk5278JfLPgvp4HPl77cxzHcarKECV0UdUngaP6qbOExIDbde8leZnsb/8LgENKd7AElTbqjuM4ziilTZSXquFnzXEcx6kWnnq1NH7WnGGnajnYnRwK/tBu3lysyWj83fYpbc5gMgpvKcdxHKfS+Jt6afysOY7jONXCjXpp/Kw5TsUYtUk+iuZQb3Ncci+AG/XS+FlzHMdxqoUb9dL4I6PjOI7jjBD8UchxKka7yO2hWl5EPbf506FgQpcRNnzeR78XpA2/2yrgZ81xHMepFi6/l8bPmuM4jlMt3KiXxs+a44QMc0aYMqPfR+2I+QhVlrWr3Den/XGj7jiO41QLf1MvjZ81x3Ecp1q4US9N4bMmIp0kid13BLYGlgOP1tLSOY7jOE7rcKNeityzJiJjgbcBJwP/AHSSzSurIvIn4HvA11X1scHqqDN6GO1Tfsr4w0v50EdaBLc2MQKj/fouhL+plyZ6RYnI0cAjwFXAeuDTwOuBvYGXAAcC7wKuJTH8vxORr4vI9MHutOM4juM4fcl7FPoP4N+Ab6vqqkide4HvAx8Xkb8HPgmcCpzf0l46juM4owd/Uy9N3lnbTVV7iu5IVX8FvF1EuprvluNUn8GSUYdtSlsROd7+0LZavS/qDmjzH3uX3AvgRr000bM2EIPeinaO4ziOA7hRb4IBnTUREWAHoM/buKo+0apOOY7jOKMcN+qlKHTWRGQ74D9JBsTF2oxtVaccxxkYwxFBrhWD5zO/2yNsMH4MH/3uDCZFH4W+AbwW+ArJiPgNg9Yjx3EcZ3Tj8ntpip611wIfVdVvD2JfHMdxHMeNehMUPWsrgGcHsyOO45SnXRK6+O909ruq6vc07LhRL01Rh86XgdNqA+Ucx3Ecx6kghR6FVPUSEdkReFhEbgFW9q2i/9Ly3jmO4zijD39TL02hN3UR+UfgQ8Aetf+fbvBXCBHZWUS+LCL3iMg6EVERmdWgXpeIfF5ElorI87X6r2lQb4yInCMiS0SkR0QeFJGjIsc+RUQeEZH1IvKoiJxWtN+O4zjOENLR0dxfAUTkhSJyrYj8TURWi8gPRWRmwbYtt1GtoKj8fgnwa5K47+NUdUzwN5DpbLsDx5C87d+VU+8bwCnAZ4EjgKXATSKyT1DvfOBckpH5hwMLgGtqDyK9iMgpwOXAD4DDgGuAy0TkgwPou9PGbGFM5m8kMSb4dFVi7Nj0z8l+VyP5mmyK+pv6IBp1ERkP3AbsCZwInAC8GJgvItsU6GVLbVSrKKpvzAQ+oqr/14Jj3qmq0wFE5GTgDWEFEdmbJFnMSar6rdq6O4BFwHnAkbV12wNnAhep6hdqzeeLyO7ARcBPa/U6gM8BV6rqXFNvR+B8EblCVTe24LM5juM4zTI08vspwG7AHvUMoyLyW+APwAdIXmYj3WutjWolRR8N7yfJo940qlrkNeJIYCNJsph6u03A1cAbRWRcbfUbSdLBXhW0vwrYS0R2rS2/ApjWoN6VwHbAQQP5DI7jOE7bcySwwKYMV9XFwN3AWwq0baWNahlFjfpHgDNF5FWt7kCEOcBiVV0XrF9EcoJ2N/XWA2Ee90W1/7NNPYCH+qnnjDA2bUr/WrqziuUidxm3fbDfU5XdJsPKEMjvJHYhtAmQ2IX+bEKrbVTLKKpv/C8wEbhTRNYCYSpWVdVdWtivKfQdYQ/JfPn69vr/VaqqBerRYJ9hPcdxHGe4GRr5Pc/OTG6ibX17/X8RG9Uyip61W4GwUyMeETmVJD88M2cWGhDpOI7jtIAWKE5TRWShWZ6nqvOa3WnVKTpP/b2D3I+QlUCjN//6U80KU2+SiEjwJNSoHiRPX0tz6mWoXQDzALq7u0fdQ81IoKUP+xWeN+vSbfuQl9DFv8cE1ZZ4uJaranfO9pU0fiOPvYWHbVtpo1pGVZ1vi4Bda1MOLLNJksk8ZuqNA17UoB7Aw6YepL71WD3HcRxndLCIvjYBErvQn01otY1qGVGjLiJvH+jORGQHETmwuS4BcAOwFfAOs+8O4FjgZlVdX1t9I8kIxHcH7Y8HHqqNZAS4B1geqbeCZLSj4ziOUwHqb+rN/BXgeuBAEdmtvqIWCO1VtW15tNpGtYw8PfHLIvJZ4GvA/6hqVCYQkVeTTNx/N/Axksn1UUTk6Frx5bX/h4vIMmCZqt6hqveLyPeBS0VkK2Ax8EFgV8zJUdW/iMglwDki8hxwH8lJPYTaPMFavY0i8hmSYDN/Am6p1TkJ+LCqeirZUcBIkDZjyUCqLONu3jzcPagWntClf1okv/fH14HTgetE5NMk48bOB54iCVQGgIjsAjwOnKeq5yX9a62NaiV5Rv3FJJPmzyMx8L8DHgSWkQzRn0wycb8beAFwJ/B6Vf1lgeNeEyxfVvt/B3Bwrfw+koAxFwCTasc+TFXvC9rOBdYAHwVmAI8Cx6jqj20lVf2aiCjwCeAs4EngdFW9DMdxHKcyDIVRV9W1InII8CWSmCVCMij8DFVdY6oKMJa+ynZLbVSrkL4j7YMKIp3A20gm0R9IEoSmC/gr8AiJMf++qj4yGB2sCt3d3brw3nuHuxujjqIjYGM/AHZsW5+3olgju97uIGeg3GDODY+9kTf9ph5+/iK/ol1dvcV1Pdlj9vQ0LlsmTIgvj+kxU35jO4D4d2L61i5z9dv9Tb37gANYuHBhy7N37rNPt95228L+K+aw3Xbym34Gyo1I+h3OW5Omv4+JnOM4juM4g0nF4ju1DdWdo+M4o5S8t7fYtnZ54/Mf6vYZCzGcDJFPfUTiRt1xHMepFG7Uy+NG3XEcx6kUbtTL40bdcSpAO8qurQiwZ3+4O5vfXVvQLoP4nPbEjbrjOI5TKfxNvTxu1B3HcZzK4Ua9HIWMuoj8EvgqSWS59f3Vd5y2o8gvSFinSf250pJ7bK5+5hwMk2Ae61ubkDf63UnwN/XyFL2iNgDfAf4sIpeIyJ6D2CfHcRzHcUpQyKir6sEkWWW+A7wHWCQit4vIsbW4t47jOI7TEoYoocuIpLB2VQsD+3EROQc4BjgV+G9guYh8iyQB/ROD003HGSYKSr1NS+lV/hWKnIPwdNjlsWPTcuGELs1K6Zl+VncsvUvu/ePye3kGfHWp6npVvZIkOP1dwDTgn4Hfi8g1IjKjxX10HMdxRhH+pl6eARl1EdlaRE4SkXuBXwPbkxj3HUnSzr0S+G7Le+k4juOMKtyol6Po6Pe9gA+Q5IndBrgO+KSqzjfVvi4iz9A3rarjOI7jOENAUSfWg8CfgUtJfOdLI/UeA+5pRcccZyAUcsfmpRotkru16ON/mTZVeLUocj7M+jHBSe/oaCz8Wf964ePnrS/wZYdjHKrkx/Ypbf3jPvXyFDXqRwPXqWrukBdV/R3w2qZ75TiO44xa3KiXp5BRV9UfDnZHHMdxHAfcqDdD+4VjckYVRaeKWRkz2qaM/F6UmORetV+mvP4U6XdPT1ruI7+n08i2ikSvyJsGhz1kmX62SXQ5l9z7x416efzqchzHcZwRQns82jqO4zijCn9TL4cbdad6lLibo5JTnqQ80CQuefJuu8jAZUbmxwjad3Sl8nvsMLmnpojkn7fznM9jR+oPt/zto9/7x+X38lT418dxHMcZjbhRL48/JjqO4zjOCMHf1J2ho50CsdSxUm/e0O2Yrlw1KT5PC7fbQsm7CfJOQaGZCnnXQ85o/OgxW/idlJHPXXLvH39TL0/FfnEcx3Gc0Y4b9fK4UXccx3EqhRv18rhRd5qnWVm9TEz1ovWKDr0uIp8HdayMWuV4M9luNx6hDoEsXWSUfFdXZtEq4bbJxo1pOQxKY+PFdxb9fouMkg/6VogSsnzR4EjOwKnafdQuuHPHcRzHcUYIQ27URWRnEfmyiNwjIutEREVkVlCnW0TmicgjtTpPish3RWTXBvtbUttH+PfWBnVPqe1zvYg8KiKnDd4ndRzHccpQl9+rlk9dRMaIyDk1u9MjIg+KyFEF2k0Ukc+KyC9F5K8isqpWbmSnzo3YtP8t0sfhkN93B44BfgPcBbyhQZ3jgDnAfwCLgJ2AzwALRWQfVX0qqH8TcG6w7lG7ICKnAJcDFwK3AK8DLhMRUdWvNvOBHMdxnNZRYZ/6+cCZwFwSG3YccI2IHKGqP81pNxP4J+BbtX1sAd4J/EhETlfV/2zQ5iDAZkZdUaSDw2HU71TV6QAicjKNjfrFqrrMrhCRu4HFwCnAZ4P6y1V1QeyAItIBfA64UlXn1lbPF5EdgfNF5ApV3Rhr79Qo4xMfrDszNzNIyhaTZCTP/xmbZpSZ5bUmvs2S50OuEn1P4RhTNr53EynOEs56s8trzLkqOqxhyiTjB7c+8WYjARalRN72yk1ZHCFU0aiLyPYkBv0iVf1CbfV8EdkduAjIM+qLgd1UdZ1Zd5OIvBD4JNDIqP9KVQd8FoZcflfVfkeWhAa9tu6PwDKSt/aB8gpgGnBVsP5KYDuSJyLHcRynAlRUfn8jyUjT0I5cBezVyD2cfh5dGxj0OguBHVvXxTYaKCciLwW2B37XYPOba7739SKyoIGfYk7t/0PB+kW1/7Nb2FXHcRxn5DEHWA88Fqxvxo68Bngksu0pEdksIn8UkYtFZOsiO2wL7agmn3+N5E39G8HmG4Bfk8gb04HTSfwUJ6hq/YlqSu3/yqDtimB7eNxTgVMBZs6c2cxHqDZl5PNWPgoXnGpmpfQ+imxkOlWWgT/DxqZmAWze3LherE7VGDs2vs26DWJfT/iZreQek9/z9mHl/4mTJhVrZCkS4a+/DsU7V6xNDJfpB0wLfmKmishCszxPVec1sb8pwCpV1WB9rh2JUbMvBwLHB5seA84G7geUxEX9MWA/4PX97bddrrSvAK8E3qSqGcOsqh+2yyLyI2AByYC4UCYZELULYB5Ad3d3+EU6juM4g0CLfOrLVbU7tlFEDgV+XmA/d6jqwU33Jnvsg0kGgv+Xqn7XbjMvo3V+LiJPA5eKyKGqekvevitv1EXkIpK35RNV9eb+6qvqZhG5BrhYRHZQ1aWkb+iTgaWmev3JqtCoQsdxHGfwGaKBcr8EXlqgXt0XvhKYVJsxZV/yBmRHRGR/4HrgNuDkgn39HnApsD/J7K0olTbqIjKXZGTgh1X1yhK7qJ/4us9jDlmjXveBPFyuh21GK0evF5FA87blSKUbNkUiteVI7FYaLyp558nPRfY1VKNzm+1n0Tax48RU6NAdEYuq9/zz8WMWcZX0keKLJJspE52u2fugjMTusvywURu4FvNnN2IRMA54EVm/emE7IiJ7kUzBfgA4qsSsq34V48oOlBORjwAXAHNV9SsDaNcBHAs8qarP1FbfAywH3h1UP57k6eru5nvsOI7jtIKKjn6/EdhIYzvykKouzmssIi8mkfufAI5Q1efz6gfUj3lvfxWH5TFRRI6uFV9e+3+4iCwDlqnqHSJyHInUcCNwm4gcaJqvVtWHa/t5J/AWkvmBT5EMlPsQyYCCd9YbqOpGEfkMSbCZP5HIF4cAJ5GoABsG6aM6juM4A6SK89RV9S8icglwjog8B9xH8gJ5CHCkrSsitwK7qOruteXtSQx6J/AvwGwRsU3uV9X1tbr3A/9FEkBNSQbHfRi4UVVv66+fw6X9XBMsX1b7fwdwMHAYILX/hwV163UgGfG+PfB5Er/GWpJ5f4ep6k22kap+TUQU+ARwFvAkcLqqXsZIoqiEOFjJVcJEGnb0uhGGMsk/AjU1JteWUUfzBkTbfcekZ7s+lI5jp8C2KRN8phWKbLM/iEUGkvdNztK4Xt53GAtYk72MsoJiV9f4hscsPsrelE1gnTGbgmf72P1S9N4pk1go1n4UUjWjXmMusAb4KDCDxPAeo6o/DuqNJWtfZwO71MphXYBdgSW18qMks7h2ILn4nwDOA/6tSAeH5apRVeln+3uB9xbYzwKSp6Six72cJFSs4ziOU1Gq+KYOyUBsErfwBf3UOzhYvp3kRbXIMY4r2T2gwj51x3Ecx3EGxujWdxzHcZzKUdU39XbAjXq7UsSv1+pEK9bHZ5yeeZHebBKUor7yMvk6irgfQ795kahpZSga2CzWptU0+1VbbFIcGwEu75iZpDgFZqOF9cJhGkVc3c8XHFdsr4GuIHFNV8zfXtTXPlgX8ijwtbtRL8/Ivzocx3GctsKNenncqDuO4ziVw416Odyotwt50l6zknueXhyR2a08aiX2vBzbrZSBc2bORaeUhW3C5Tp5edcLUSZhyCD+gDXOhk7fkxvrT+TzdAbtbT70mDS/Zk3D1X26U0Zhtt0Mpx/GrsO8Y44bl5a33jo9ixMmFJDlw4PG1ud9B0VPyCiQ453i+NXgOI7jVAqX38vjRt1xHMepFG7Uy+NGvcoUHVlbVH6PyXkRiR2ycqkdUbx+fVouI7HnKYsFutkHu8226ewwUnofebTgeWsl7fJLVSbioPkSsklYxjSqAsSV6AkTGpch+E4zpMcJk81YYqPxw+vQblu7Ni3b+8DK8uHo+c4JZocxKb7V18MIkeLdqJfHg884juM4zghhZDzWOY7jOCMGf1Mvjxv1qhHTBssEtsgZ7m1ldiuxh6OTM6Pcm5TZbXdsOQwKs/XWjevZffUZoR6NdlIw8kmRhBsjRNqMEvtSi37xkW0xKR7iLhUruXcSJFpZ1XgI/Xizs65p4zPb7CyIVavSct5Hi0nzscA44Wexo+SjI+aLThfJuz5jN1kbX69u1MvTvt+64ziOM2Jxo14ON+qO4zhOpfA39fL4QDnHcRzHGSH4m3oVKDI9rehUNes378r6FWO+c1suGhSr6LSzmM/Rlsd35fjHy0wFKjodK+anLOJfD4m1z6s3HBTtW5FzUHSapfkOu4JrMtadzLS10Ie+JhLC0PRzTMaPD1PMso12Vx3aA4gAACAASURBVHQ8SZEENeFpsvvLTtFr7GsHGBNz3rcyamQb4G/q5Wmvb9pxHMcZ8bhRL48bdcdxHKdSuFEvjxv1KlBEcs+bxmK0vXU9qbS46plsk5hqaQ8TRuKKJUSx5W22Sct2OlpYLzM1yXZgVUHN39IKOTHmQyhSLrrfqlG0b80mRC94+FI/3LYPsQwxOVPF7BS7CTNSd4Cd6pZ3mNi9E+Zwt8uxfYXR8iZMSPsz0STIic4tDRkh8ju4US+LD5RzHMdxnBFC+z2+OY7jOCMal9/L40Z9OCh6tVrtOtDpYjJ73mheK63bfNM2olson9vDWpl9220bd3NMz7rsDtZEtMoy0d2azSldtP1oyWPd7Ej4vO8wcm6KDpi3IuL4cEpF7LwXjcBoymPMBT4luMe2TOo/6mJMloes/F6kPWTvtzUT0nMwaVIqy4+flJMVp+jMj4pfu27Uy1Ptb9ZxHMcZdbhRL4/71B3HcRxnhOBv6lXDSI0bSOW/lcuy1f72t7S8cmVatnnOQ2JBYqzk94IXZNtYmT0TJMbqhssHLnvmUiSyTd5sgDIUkNy3jLRn4I4g8EmYJKc/Qlm8hNvEStTWJWRdRWFylkxgmaIBWmLDz3NcQmNsfnhbnpGWM26wYPR8Eck9r824cY3rTZqUvQ6nTk3PT2eHmWHSxq+6/qZeHjfqjuM4TuVwo14ON+qO4zhOpfA39fK4UXccx3EqhRv18gy5UReRnYFPAt3A3sDWwK6quiSop5Fd7KuqD5h6Y2r7+wAwA3gUOE9Vf9Dg2KcAnwB2BZYAX1LVrzX5kZqnq7GPbvlyGpYh62OLXfx9o1WlZeuWnDo1LY+3PjkIHIAlEkwU8a3m+cdjfvTQnxujhB8/5jtvs1lBAyb2uce08IPmfR2xGXI2qiFkk7NEferhzmJT9PLaFMjcYqfbjZ+RveFWr2nsb7ddKRrFbu3atPzcc/E2M2ak4yTGdxVMMuQUZiD2pkHbbwMnNtj076p6RlD3IODfgH2BvwH/DcxV1ecbtM8wHD9LuwPHAL8B7gLekFP328DlwbrfB8vnA2cCc2v7PA64RkSOUNWf1ivVDPrlwIXALcDrgMtERFT1q6U/jeM4jtNSKvymXsje5LAMODJYt9QuiMjLgJ8DNwFHkLyEfh7YCTi2vwMMh1G/U1WnA4jIyeQb9T+p6oLYRhHZnuQEX6SqX6itni8iuwMXAT+t1esAPgdcqapzTb0dgfNF5ApVDaKeO47jOMNBFY16UXvTDxvybFqNfwWeBt5Rt0sisgH4johcrKr35TUecqOuqgOcN5PLG4FO4Kpg/VXAN0VkV1VdDLwCmNag3pXA+4CDgPkt7Fc+gZy5YVNjmW6peX6zU9ggPnvISuxWVg+Xp0yy+arNQXMSYRTKsZ0npTcpv28xU7CKKv5FpeOiknts20iT4i323BSW4iMnLq95TBUPJWqbD90mZ8k0CkO1FZlul9e5mG8gJ7nMRHMzTth5Ym/ZTlULD2ndbLFpcEWv/alTbVS+9rtAq2bUKW5vSiMiWwGHAZcEL5r/A3wdeAuQa9SrPvH2gyKyXkTWichtIvLqYPscYD3wWLB+Ue3/bFMP4KF+6jmO4zjDTP1NvZm/QaCovcljexFZLiKbROT3IvJJETGBunkR0EVgq1S1B3i8yDGq/Ph2FfBj4M/ALsBZwG0i8npVvb1WZwqwSlXDQXUrzHb7f2U/9TKIyKnAqQAzZ84s8REcx3GcYWKqiCw0y/NUdV4T+ytqb2I8QOKHX0RiuN9GMsbrxcDJwT5CW1U/Tn/HqK5RV9UTzOJdInIdydPLBSRy+VD0YR4wD6C7uzs2Gn/AWLkd4hGm8qLDFRnJHsrvnZtMspVnIpJ7mZHszUZ6y2kTk9yLjkTPSMcDjZjWDyNZcm8p5svq6MpGsbMj22O5YkLlO5OP3ORDt8lZciMbFpXfi4yYz7sobeIY8yGmT05v2G23zZ4P+xFiUnx4PuzvhP39yMr82d+czo7W3guDQQs8tctVtTu2UUQOJRmQ1h93qOrBzXZGVS8NVv1URNYAZ9R85X9o9hhQYaMeoqrPichPgPeb1SuBSbUR7Nbo1p9mVph6AJPJjjQM6zmO4zjDjgKb+63VJL8EXlqgXv1tqKi9GQjfA84gmeL9B7K2KmQKqdQfpW2MusGezEXAOBI/hPVz1P0OD5t6kPhElubUcxzHcSrB4Bp1VV0HPDKAJkXtTanu1P4/TuK3n2M3ikgXsBtwTX87ahujLiITSebs3WtW3whsBN5NMg2gzvHAQ2Yk4j3A8lq9W4J6K4C7B6PPRUdUx0a2WvnMSuyQldbttvEYiX15TmSLmGwYSpBFgr/kBYIpoVEXldyL0ErJ3eX2gPCEFMjlPWZTNrhRl5HjY0p43v1i5eZMPvS8WRxF3Ejhtti+igaviejnYa74mTPSm3nq1PTc2M+ZF7DGYmX5vu6qwXNLjWCK2puB8G4Sg/5rAFXdICI3AseIyLmqWv/mjiZ5oLi+vx0Oy8+UiBxdK7689v9wEVkGLFPVO0TkTGAPkmlm9YFyZ5JE8Hl3fT+q+hcRuQQ4R0SeIxnqfyxwCGaCv6puFJHPkASb+ROJYT8EOAn4sKoGYdQcx3Gc4WNI5PcBUdTeAIjIrcAuqrp7bXkXkinUV5O85Y8jGSj3XuByVX3cND8XWAD8j4j8JzCLJPjMtar6m/76OVzvHqGEcFnt/x3AwSSh995W+3sBsJrkbfr9qnpv0HYusAb4KGnYvmNU9ce2kqp+rRZ69hMkI+mfBE5X1ctwHMdxKkYlFYRC9gYYS9a+PkeiCn8SmE7y4R4BPkJq/wBQ1QdE5A3AxcBPSMLE/hfwqSIdHBajrqrSz/YbgBsK7mszyYj4CwrUvZy+YWcdx3GcSlG9N3Uobm/C0fKqugJ46wCOcydJ0LQB417CEuRNk4r50S15vmHrYsudntazOl14xsx9KTo9LeYTD/2KsW1FI3FZYv3JaV/U9T9YuB99AOQ5wiPrrU/djiHJm9IW27ZlUrqvMeF1XOTaKzOeJM+nHutoLFQcZOaxZZLFmIEzO+6ZHWDz7LL0N8cmfol1BcKPXUX/ejWNejtQ9YhyjuM4juMUxN9DHMdxnArib+plcKPeJEXkdshXwmNT1yZ2mUH5Tz+dbbQ8IrnnyYl2yo+VEGPrg2X7WTNSOHHKyHnNxm2ujoQ4iikSjY3sFLett07lcysj582ciynZE8PrOCaFF4ygmLn2MdPwunIk/1Bab3TMvPlptp7ddzC/dfqMGb3ldTPSKKJh3vUiFI3AWPR3rzwuv5fFjbrjOI5TQfzhvAxu1B3HcZyK4W/qZXGjPgBaJTmFcuK226blTES4JUZyf+aZbKNYGKmYrB4uR8o2mhsMPKJbn9HiBRqF57VMFLnoKPXRkvS8XchJejJhQnrtPf98WiUvUFtsUPnEGQWjwxVMyGIl99jxu7qCpCnh/ddoB2V8Czm54sfPSvfXMTnN4R77uchj8CV2ZzDwXznHcRynYvibelncqDuO4zgVxI16GdyoDxF5ym8nZpS7DSRjR8bmJRCPSe558rsd1R5JoNJouRG2K31GzJYIOFPkOGHzzHELavYuLw6c3NHRRSTunG2xkfB5A9ljavW6nux3m0mcUjQhi2ljR7nHpOw+640c3xkL4hTeo7FzZf0RGzdmt9kPbmbFdM6wvxHjaS/8Tb0s/qvmOI7jOCMEf1N3HMdxKohPaSuDG/VhoLMjuFhXNc61nJEzQ5muiJxXMI77oI02L7iDZiV/p+LkxYSPXHw2JnxemvPYKPkwpsv4Gea+sHJ1wdjtY8xBuwpI8X3Iu5dj9Yq4DMI2kekAnVOzJ7H6rieX38viP42O4zhOBXGjXgY36o7jOE7F8Df1slRdg3Ecx3EcpyD+pj4AGiU4KOqbyrTNy7tsKeEfz533FUlQEamSS2Yam5mKVCpBRo6bNfZxBnPqnDOI5EVQM/dB54S0no00B1mXuPWp23IYdG31mvR6m2jvq7wBJZHxLbZvdBW7/zNTAYuOlbGJW/J86rHfj5zxAmMiiZqqg7+pl8V/8RzHcZwK4qPfy+BG3XEcx6kY/qZeFjfqg0jhvN6xaSx5kacK5Dnvc5gChy9KVHLPm25jyrH01nm4kt5GlIkuZzEXRSi/21shltAllN/tFLcJO6eJTsbEdhDuPEJnzPVFTvS9PBdZEfLcd0Wk+ADbt1LuRKdS+M+k4ziOU0H8Tb0MbtQdx3GciuHye1ncqLeYUrJUbDRsTp7z2GDYMmpe4VHlTYae27CpmJsg1rdovwrS6lG+MdkyT84caW0GTN4FGpGOM64esnJ80ZTj9rDjxqXl6ZMjI8zDncSu/ZwZKmMiMnvRmR8xOjqyyVk6J5njrolEpwyJ/Gi09LtuGjfqZXCj7jiO41QMxUe/l2O4H8Ucx3Ecx2kR/qbeJIXl9oJBWTIyXUc8cURsdy2XtWM6f0zzDyTIzGeIBA0J2XrrxrsuNZtgEInJk3myZZE24ecsIoXntWll30qR933ErqOcwClWfp8QydMSJnSxy/Yw226b7mt8GOAppufH7olwhopdjpyDgvltMvTZVSxve5MM/wh3l9/L4EbdcRzHqRg+UK4sbtQdx3GciuFGvSxD7lMXkZ1F5Msico+IrBMRFZFZQZ1za+sb/fUEdZdE6r21wbFPEZFHRGS9iDwqIqcN7qd1HMdxyrGlyb/WIyJjROScmt3pEZEHReSoAu1m5dg0FZHjTN2Y/fvfIn0cjjf13YFjgN8AdwFvaFDnCuDGYN02tXXXN6h/E3BusO5RuyAipwCXAxcCtwCvAy4TEVHVrw7sI7SYHJ/0YBGdBpfn5AvnCdXJiapl3Y/Wj77ZPISPHZvdnd1FZ4e5OUtMY2slw+VjjB03rz/D7w8tSMyPbgmTkZg2EyaMN+W0Tl5wuOXL07JtM3OGmd4WVoxd+3nXZORezvObF7nEcxMgdaVjBMZ0DO/9MkI5HzgTmEtiw44DrhGRI1T1pzntlgKvaLD+AuAgEhsWchBZuWJFkQ4Oh1G/U1WnA4jIyTQw6qr6NPC0XSciJ5D09zsN9rlcVRfEDigiHcDngCtVdW5t9XwR2RE4X0SuUNWNpT6N4ziO02KqJ7+LyPYkBv0iVf1CbfV8EdkduAiIGnVVXQ9kbJSIjAcOAG5Q1ZUNmv1KVQf8ZDbk8ruqln2FOBF4lsZPNP3xCmAacFWw/kpgO5InIsdxHKcybG7yr+W8Eeikrx25CthLRHYd4P7eDmxL4xfV0rTFQDkReSHwWuDSyJPLm0VkHTAWuJ/kScr6H+bU/j8UtFtU+z8bmN/CLvelxDSrMpHi8uS7QvvLm9IWyS9t6RP5zqiW69c3PmSoUkbdASVOSLPTsdpGxm4X8r63oklgzHU40URTWzMh/a5t1DiIz06zCvvUqdlrNzPFzfbB+pHs+vBCtveL0fnLBGa0FP4pKeLaqCzVe1MnsSPrgceC9daOLB7A/k4E/kJfV3Odp2rqwNPA1cC5qpozGTihLYw6cDyJqtDoieYG4NckJ3M6cDrwIxE5QVXrT1RTav9DiWNFsD2DiJwKnAowc+bM0p13HMdxhpypIrLQLM9T1XlN7G8KsEpVNVifa0caISI7AYcA/97gRfUx4GySF1QlcVF/DNgPeH1/+24Xo/4e4H5V/W24QVU/bJdF5EckvosL6SuTDIjaBTAPoLu7O/wiHcdxnEGj6Tf15araHdsoIocCPy+wnztU9eBmOxNwAsmL6rfDDeZltM7PReRp4FIROVRVb8nbceWNuogcAOwJnFGkvqpuFpFrgItFZAdVXUr6hj6ZZBRinfqTVaFRhYNGJNlDnnoW21ZGcctIzEUjykWi4BVtnjNgPpDfG/e5D01GkXOZvQIUld8j0eYmTUpHwocR5dauTcsxKT5sM36SGQ3fFUmaYqX4cIS8bZ/5DFmZv5XYw3QOUWTFwWFIYr//EnhpgXrrav9XApNqM6bsS14ZO/Ie4IFGL6oRvgdcCuxPMnsrSjt86ycCG4H/LtG2fuLrPo85ZI367Nr/h8t1zXEcx2k9g+9TV9V1wCMDaLIIGAe8iKxffUB2RET2J3mY+NgAjl2nX8W40gldRKSTZB7gz1R1WcE2HcCxwJOq+kxt9T3AcuDdQfXjSZ6u7m5Njx3HcZzWULnR7zeSvGA2siMPqWrRQXInkmiQA3lRrR/z3v4qDsubuogcXSu+vPb/cBFZBixT1TtM1SNIpI2GQ/5F5J3AW0jmBz5FMlDuQyQDCt5Zr6eqG0XkMyTBZv5EIl8cAvz/9u49WJKyvOP497cssLuCXFxwuQi7hCrjEsSYTSQFyrIBQe6WSEJiYkmEiLhiEFKSmBQKiopSxAQSiAooJHIT2UVLdLksCZcYAgFFF4OChIsslwWWZS+AT/54ezh9muk5feb0nOkz8/tUdc05PW/3TL+n5zzzvv32+xwDLI6IDcV991QX3WLFTap0uVceaFx1OG5Jl3t+ZG+nxDP5TfITzhQHDZdeDpikEe82iar8TTud/LkTblZuJPyWW44+B1avHnPzV3W/b/+bue7zLUvyrnf67JSUm567Q6TToeWVffTGes7qExErJZ0NnCppNXAnqQG5CDgsX1bS9cDOEbFrYX2+obqy3etIugv4OmkCtSANjlsMfC8ibhjrffbrFLii8Pt52eNyYGFu/ftJLelrS/bzALAtcBYp+K8B7gAOjIhR97NHxD9LCuDjwCnAQ8BHIuI8zMysQRp5SxukmeSeB04E5pAC71ERUYxRG9E+vh5Mmhul073p95Hu4tqO1Jv+C+DTwBeqvMG+BPWIUMVyh4/x/O2kb0lVX/d80lSxZmbWaM0L6hHxMmlq1zPGKLewZP3VQMf4FxF/1On5sbizxszMGmZSRr8PJAf1hslfT+50bbib6+h5oxKlPF9yjbC4w/zPuWuM+Vnk1hbmOyq7jr755u3LvOo9lB1QhwMtq0PfttZwXfytR8ldw549e1bZU6PkZzkslnn8iZFz5/U77jjyRFmSo04XyHPlNtkyn8Bp9O1tE5xAccA0r6U+FXhEkZmZ2YAYyu9/ZmbWZI0dKNd4DuoDoGrX3Kju5yozxRV/z997lruNrSzVdHGTsl1V7hbv4pY2d7lPUWUnC5T3Ued+3mT66LtU58xpfxtZ8Ta2vPwsdC9s/tpXfp41d+7IE/mMMGV9/MXn8l3x+YTujM6HPuVysNTKQb1bDupmZtZADurd8DV1MzOzAeGWeoNNuOu4UyKMCslZXvV7rht0w0vtvw/OnFlp885JZMp2UJG73AdMN0O/C+fUrBkj+5g9u33e9fxI+KL8jHTku+Ln5N5b8TpU2XldNqUdMG167hJC7gPT6U6Yweymd/d7txzUzcysgfzlvBsO6mZm1jBuqXfLQb0fOk3wUsf+qqiYBSI/sUzZyxQHJ5e9TGm3+HDOrGHdKjt3O30Ocs/lu+KnT5/WrghQPpg9X27DjJFJbjaZXfggVE38Uia3/bQOn9f8MfTS5N9V4qDeDQ+UMzMzGxBuIpmZWcO4+71bDupmZtZAHijXDQf1qaKO+1YqXEevmkSmKt9eZpOmixM0n9ioeG16vLsrfnam5QebNO2+s8aPY3FLvVtN/8uamdnQcVDvlgfKmZmZDQi31LvQqYu6q+7mCeYPn6hOxzNe7m63qap47k70VrFRt4D18PPbq89cnf8XuuOWejcc1M3MrGHc/d4tB3UzM2sg9/p1w0G9ZhOedalKV3zDuMvdBlH+vJ5oV3RPZ2Mrm/qu6qyRfe9mtzpNjahhZmZDxN3v3XJQNzOzBnJQ74aDehNUmZii5q74iXa5ucvdhkmju+Lzedzz/yfyk98U/n/UeQy94ZZ6txzUzcysYRzUu+UREmZmZgPCLXUzM2sgt9S7MaktdUlHSrpK0i8lrZV0n6QzJW1eKLeVpK9IelLSGknLJO3eZn8zJJ0l6bFsf7dJekebctMknSrpQUnrJN0t6T29PFZI151aSy/33cvXaf8qvp5uNpk6fsbXrRtZXnppZJnSgnSf+kSW+kk6SdLSLOaEpNPGuf3ekm7N4tWvJJ0taWabcrtJ+r6k5yU9JelCSVtXeY3J7n4/mfT166+BA4F/Ao4HfiBpGoAkAUuz5xcD7wE2Bm6UtGNhf18FjgX+DjgEeAy4TtJbCuVOB04D/hF4F3A7cIWkg2o+PjMzq8XLE1x64lhgW+Db491Q0puBHwArSfHqk8AHgIsK5bYHbgJmAkcCJwD7Ade24mQnk939fmhEPJH7fbmkp4GLgYXADcBhwF7Aooi4EUDSbcADwF8BH83W7QH8MXBMRFyYrVsO3At8OtsPkrYlfZn4XER8MXvdGyXtCnwO+G7PjtbMzLrQ2IFyu0XEryVNBz40zm0/BTwMvDciXgSQtAG4WNLnI+LOrNwppIbsoRHxTFbuUWA5cATwrU4vMqkt9UJAb/mv7HGH7PEw4NFWQM+2e5bUej88t91hwIvAZblyLwHfBA6QtGm2+gBgE+CSwuteAuwuaV53RzMB06ePf6lZWbd6cell177ZVDHRz0FPL2NtttnI0sP/GQYR0dUfTdLGpN7ny1sBPXM5sIFXx7bvtAJ69ro3Aw8VyrXVhP/U+2SPP80edwN+3KbcvcBOkjbLlXsgIl5oU24TYNdcufXA/W3KAczv8n2bmVlPtFrqjet+79ZvADMoxLaIWAf8nCwOZdfX5xXLZe6lQrzqa1CXtAOpq3xZRNyRrd4aWNWm+NPZ41YVy22de3wmImKMcmZm1hjNGyg3Aa04UxazWs9vBahCuVJ966PJWtzXAC+RBgs0jqTjgOOyX9dro43afXuy6mYDT/b7TQwA1+PEuQ7r8cbe7PbZ62Dp7AnuZIakO3K/XxARF7R+kbQfaeDaWJZHxMIJvpdJ05egnnUxLAV2AfaJiIdzT69ipDWeV/ymswrYuUO5p3PltpSkQmu9WO5VshPgguw93xERC8rK2thch/VwPU6c67AehaBZm4g4sBf7LbgVeFOFcsVLvN1oxa2y2Na6HPwM6dpDWbnSeNUy6UE9GzBwJbAA2D8iflQoci/wzjabzgceiojnc+XeLWlW4br6fNLAg/tz5TYlXdO4v1AO4CfdHouZmU1NWdxYMUkv93PS2K7d8islzSA1bq9ovSdJDxbLZeaTRsB3NNmTz0wDLgUWAUdExO1tii0BdpC0T2671wKHZs+1LCUN+39vrtx04A+B70fE+mz190ij5P+k8DrvA34cEQ9M6KDMzMw6iIgNpFh0VBanWo4kNTrzsW0JcLCkLVorJO1N6pnOl2trslvq55KC8GeANZL2zD33cNYNvwS4DbhE0imkbotTSYMHvtAqHBF3SboMOCdr/T9AmshmHrkAHhErJZ0NnCppNXAnKfAvIruXvaILxi5iY3Ad1sP1OHGuw3oMVT1KWgDMZaRBPF/SkdnP3231Gkv6KvD+iMjH2NNIE59dLuncbD9nAVdGxH/nyp1FanQukXQmsAUp9v0ncPWYbzIiJm0BHiRdL2i3nJYrtzXwNdL1gxeA64E92uxvJnA28CtgXXbQC9uU24g0e88vSV0g9wBHTuaxe/HixYuXqb2QZn8ri2Fzi+XabP8OUqN1HfA4cA4wq0253UmD+NaQGrYXAa+r8h6V7cDMzMymuCZMPtNYkt4g6UpJz0p6TtK3JO3U7/fVb5IWZskMisszhXK1JuaZyiTtKOkfsmN7IauvuW3K1Z6kSNKxklZIWq+URGm801s2xjjqsd35GcW8EMNYjxqyxFpDp9/dGU1dgFnA/5Jm9jmCND3fj0ijGF/T7/fX57pZSOpuWgzsmVsW5MoI+A/SXMdHk6ZIXE66N3jHwv4uJd3KcSzwB6S5jdcCb+n3sdZcZ4+Tcg1cR6G7brx1QRqXsp6U12Bf4HzSjBsHFcodm63/TFbujOz34/tdJz2uxwAuLJyfe1Lo6hzGeiS7rksae7QP8LHsnLsdmJaVqf3zW7WuvUzw79vvN9DUBTiRNNfgrrl180iT5ZzU7/fX57pZmP3T3K9DmcOzMvvm1m1BGifx5dy6PbJyH8itmw7cByzp97HWWGfTcj9/sF0wqloXpCxR64FPFba/HrinsO1K4OJCua9l/5w37ne99KIes+cCOGOMfQ1lPQLbtFn3Z1mdLcp+r/XzW7WuvUx8cfd7ucOA2yPilXvbI93+dgsVJtW32hPzTGlRLRFE3UmKfh/Ypk25bwCvA/YezzE0QcV6rGoo6zGcWGugOaiX65RYxklgkkslvSzpKUn/WhhvUHdinmFQd5Ki1gQWxb/DsCQzOj67/v2CpBskvb3wvOtxhBNrDQgH9XKdEsa0m8JvmDwLfInU/bkIOB3YD7hNKX891J+YZxjUnaSoLInEMNTtJcCHSeflcaQW9Q2SFubKuB5xYq1B46S7Nm4RcRdwV27Vckk3Az8EPkqaE8CsbyLiT3O//ruka0gtzzOYYt3lvaQpkFjLxsct9XKdEsu0+2Y61CLiTuBnwO9mq8aTmKdTuTETGAyQqnXxSpKiCuVos8+hq9uIWA18h5HzE4a8HjU6sdYB0X1irTrPWZsgB/Vy91I+qb6TwJRrda91qr9iYp55kma1KZdPzDMMqtZFPklRsRyMnJ+t65XFv8MwJzPKd/8ObT1qdGKtg6J9Yq06P79V69omyEG93BJgT0m7tFZkk1zsRYVJ9YeN0pzIbyR1wUP9iXmGQd1Jim4j3XLVrtzTpDs5hkJ27h3CyPkJQ1qPcmKtgeZr6uX+BfgIcI2kT5K+4Z8O/B9p0oShJelSUgKdO0mTTvw2KenOI8CXs2K1JuYZBBpJ/PA72eO7JD0BPBERy6vWRVRMUhQRL0r6W+A8SY8Ay7IyxwCLI2WOmnLGqkdJJ5O+YN4IPErKbnUyMAfXI0ztxFo23R94uAAAAzxJREFUln7fKN/kBdgJuAp4DlgNfJs2E10M20L6cN9DGgX/IumLzgXAdoVytSbmmeoL5YkgbhpvXTCOJEXAX5DGO6wnzZL44X7XRS/rkdSavIXUun4ReIoUpH7P9RjgxFoDvTihi5mZ2YDwNXUzM7MB4aBuZmY2IBzUzczMBoSDupmZ2YBwUDczMxsQDupmZmYDwkHdzMxsQDiomzWcpNdIejQ3k9pE9zdT0mOSjqpjf2bWHA7qZs33cdLsaFfVsbOIWEua6vOz2dSeZjYgHNTNGkzSpsBi4Pyod/rHi4A3AO+ucZ9m1mcO6mY9lHWdr5D0w3yrWNI7Jf1a0glj7OII0hzclxX2e5GkhyUtkHSrpLWS7pN0cPb8SZIelPScpGskbZPfPiJWAdcBH6zlQM2sERzUzXooItYARwN7kLL8Ien1wNeBpRFx7hi7OBD4aUQ82ea512b7+Qqpxb0SuErSl4B9gROAj2U/t3udm4F9JM0Y73GZWTM59apZj0VKT/kJ4IuSlpHSgL4M/HmFzfckpalsZ3PgQxFxM4CkR4G7SXnD50fEy9n63wIWS9qotS5zF7AJ8Fbg1vEfmZk1jYO62eQ4B9gfuJYUSPcvaX0XbU/qJm9nTSugZ1Zkj8sKwXsF6bO+HfBwbv0TudcwswHg7nezSZANcvsGsClwd0RcX3HTGaTc0+08U3iNDdmPqwrlWuuL3exrs8eZFd+LmTWcg7rZJJA0B/h7Ulf6HpJOrLjpU8BWPXpbW2ePVXoMzGwKcFA36zFJAi4mtbj3I3XFf17SmytsvgLYpUdvbV72eF+P9m9mk8xB3az3TiIF8/dlt5J9AvgJ8G+Sxur6vhlYIKkXn9W3AY9ExC96sG8z6wMHdbMekvRW4LPAmRGxHF659n00MBc4e4xdXAZsAby9B2/vEOCbPdivmfWJ6p2kyszqJukm4P6IqG2iGElvI93G9qaI+Fld+zWz/nJQN2s4SXsBy4BdI+KRmvZ5NbAqIo6pY39m1gzufjdruIi4BfhLYOc69pddx/8f4G/q2J+ZNYdb6mZmZgPCLXUzM7MB4aBuZmY2IBzUzczMBoSDupmZ2YBwUDczMxsQ/w9EMhNvcStaKwAAAABJRU5ErkJggg==\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",
    "p.data[:, :, :3] = 0 #  Mute out mirrored wavefield above free surface\n",
    "\n",
    "scale = np.max(p.data[1])\n",
    "# fig = plt.figure()\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."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}