{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 01 - Introduction to seismic modelling\n", "\n", "This notebook is the first in a series of tutorials highlighting various aspects of seismic inversion based on Devito operators. In this first example we aim to highlight the core ideas behind seismic modelling, where we create a numerical model that captures the processes involved in a seismic survey. This forward model will then form the basis for further tutorials on the implementation of inversion processes using Devito operators.\n", "\n", "## Modelling workflow\n", "\n", "The core process we are aiming to model is a seismic survey, which consists of two main components:\n", "\n", "- **Source** - A source is positioned at a single or a few physical locations where artificial pressure is injected into the domain we want to model. In the case of land survey, it is usually dynamite blowing up at a given location, or a vibroseis (a vibrating engine generating continuous sound waves). For a marine survey, the source is an air gun sending a bubble of compressed air into the water that will expand and generate a seismic wave.\n", "- **Receiver** - A set of microphones or hydrophones are used to measure the resulting wave and create a set of measurements called a *Shot Record*. These measurements are recorded at multiple locations, and usually at the surface of the domain or at the bottom of the ocean in some marine cases.\n", "\n", "In order to create a numerical model of a seismic survey, we need to solve the wave equation and implement source and receiver interpolation to inject the source and record the seismic wave at sparse point locations in the grid.\n", "\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The acoustic seismic wave equation\n", "The acoustic wave equation for the square slowness $m$, defined as $m=\\frac{1}{c^2}$, where $c$ is the speed of sound in the given physical media, and a source $q$ is given by:\n", "\n", "\\begin{cases}\n", " &m \\frac{d^2 u(x,t)}{dt^2} - \\nabla^2 u(x,t) = q \\ \\text{in } \\Omega \\\\\n", " &u(.,t=0) = 0 \\\\\n", " &\\frac{d u(x,t)}{dt}|_{t=0} = 0 \n", "\\end{cases}\n", "\n", "with the zero initial conditions to guarantee unicity of the solution.\n", "The boundary conditions are Dirichlet conditions:\n", "\\begin{equation}\n", " u(x,t)|_\\delta\\Omega = 0\n", "\\end{equation}\n", "\n", "where $\\delta\\Omega$ is the surface of the boundary of the model $\\Omega$.\n", "\n", "\n", "# Finite domains\n", "\n", "The last piece of the puzzle is the computational limitation. In the field, the seismic wave propagates in every direction to an \"infinite\" distance. However, solving the wave equation in a mathematically/discrete infinite domain is not feasible. In order to compensate, Absorbing Boundary Conditions (ABC) or Perfectly Matched Layers (PML) are required to mimic an infinite domain. These two methods allow to approximate an infinite media by damping and absorbing the waves at the limit of the domain to avoid reflections.\n", "\n", "The simplest of these methods is the absorbing damping mask. The core idea is to extend the physical domain and to add a Sponge mask in this extension that will absorb the incident waves. The acoustic wave equation with this damping mask can be rewritten as:\n", "\n", "\\begin{cases} \n", " &m \\frac{d^2 u(x,t)}{dt^2} - \\nabla^2 u(x,t) + \\eta \\frac{d u(x,t)}{dt}=q \\ \\text{in } \\Omega \\\\\n", " &u(.,0) = 0 \\\\\n", " &\\frac{d u(x,t)}{dt}|_{t=0} = 0 \n", "\\end{cases}\n", "\n", "where $\\eta$ is the damping mask equal to $0$ inside the physical domain and increasing inside the sponge layer. Multiple choice of profile can be chosen for $\\eta$ from linear to exponential." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Seismic modelling with devito\n", "\n", "We describe here a step by step setup of seismic modelling with Devito in a simple 2D case. We will create a physical model of our domain and define a single source and an according set of receivers to model for the forward model. But first, we initialize some basic utilities." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "#NBVAL_IGNORE_OUTPUT\n", "# Adding ignore due to (probably an np notebook magic) bug\n", "import numpy as np\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Define the physical problem\n", "\n", "The first step is to define the physical model:\n", "\n", "- What are the physical dimensions of interest\n", "- What is the velocity profile of this physical domain\n", "\n", "We will create a simple velocity model here by hand for demonstration purposes. This model essentially consists of two layers, each with a different velocity: $1.5km/s$ in the top layer and $2.5km/s$ in the bottom layer. We will use this simple model a lot in the following tutorials, so we will rely on a utility function to create it again later." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "output_type": "stream", "name": "stderr", "text": [ "Operator `initdamp` ran in 0.01 s\n", "Operator `pad_vp` ran in 0.01 s\n" ] }, { "output_type": "display_data", "data": { "text/plain": "
", "image/svg+xml": "\n\n\n \n \n \n \n 2021-04-15T19:25:24.169700\n image/svg+xml\n \n \n Matplotlib v3.4.1, https://matplotlib.org/\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAc8AAAGDCAYAAABN4ps8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAA0TklEQVR4nO3debgcRb3/8feHBGTzhyACsoaIYoIbGjWCF5IoAiLggoiCLLKJG6iIyBIg4AVFAblcgQDKei/IJosgyJKgKELcgABhSwjhBklI2EMg4fv7o3rIMJlzpvvMTJ8zcz6v55mnz3R3VddUTs53qrq6ShGBmZmZ5bdMfxfAzMys0zh4mpmZFeTgaWZmVpCDp5mZWUEOnmZmZgU5eJqZmRVUevCUtJ6kyyQ9K+k5SVdIWj9n2uUlnShptqQFkv4iaYt2l9nMzKxaqcFT0orALcC7gT2ArwLvBG6VtFKOLM4B9gXGA58BZgM3SPpAWwpsZmZWh8qcJEHSgcBJwMYR8XC2b0PgIeCQiDipl7TvB/4JfC0ifp3tGwpMBaZFxA5tLr6ZmRlQfrftDsAdlcAJEBHTgduBHXOkfRW4pCrtIuBiYGtJb2p9cc3MzJZWdvDcBLi3zv6pwMgcaadHxEt10i4HbNR88czMzBorO3iuBsyvs38esGoTaSvHzczM2m5ofxeg3STtB+yX3i37IVi9X8tjZjawPEPESwLYSFqqa6+o2XBDRGzTgoINaGUHz/nUb2H21KqsTbtBD2lhSQv0DSJiIjARQFo7Xo+jZmZG9ucRgAXAN5vM7YhB0kIpO3hOJd27rDUSuC9H2s9JWrHmvudI4BXg4frJzMwsDwHL9nchOkTZ9zyvBkZLGl7ZIWkYsHl2rDfXkP5dv1iVdijwJeDGiFjY8tKamQ0iIrWomnkNFmUHz7OAGcBVknaUtANwFfA4cGblJEkbSFokaXxlX0T8g/SYyimS9pH0CdJjKhsCR5X4GczMbJAr9YtCRLwoaRxwMnAB6YvOzcBBEfFC1akChrB0cN8L+DFwHPAW4F/ANhHx9zYX3cys67nbNr/SW9kRMRP4QoNzZpD+HWv3LwC+l73MzKyFKt221pjryczMALc8i/CSZGZmZgW55WlmZoC7bYtwPZmZGeBu2yIcPM3MDHDLswjXk5mZAW55FuEBQ2ZmZgW55WlmZoBbnkU4eJqZ2escFPJxPZmZGeCWZxEOnmZmBni0bREeMGRmZlaQv2SYmRngbtsiHDzNzAxwt20RriczMwPc8izC9zzNzMwKcsvTzMwAd9sW4XoyMzPA3bZFOHiamRnglmcRriczMwPc8izCA4bMzKw0knaSdLmkxyQtkDRN0vGS3lwwn0MlhaQ/tausvXHL08zMgNJangcDM4HDgFnApsDRwFhJm0XEa40ykDQcOAJ4qo3l7JWDp5mZva6EoLB9RMypej9Z0jzgPGAMcEuOPE4HLgI2pp/imIOnmZkBWcuz2aiwqPfDNYGz4q5su06j7CV9Bfgg8GXgioKlaxkHTzMzA0CCoW0Onj3YMtve39tJklYFTgYOiYh5kvp0sVZw8DQzs1ZaXdKUqvcTI2JiTydLWgeYANwUEVN6Oi9zIvAgcG7TpWySg6eZmQGp5bnskKazmRsRo/JdTysDV5Haq3s1OPc/gN2BD0ZENF3KJjl4mpkZ0KJu29zX0grANcBwYMuImNUgyZnAOcAsSW/J9g0FhmTvF0TEwjYVdykOnmZmBrRowFCe60jLApcBo4CtIuKeHMlGZK+v1zk2H/gucEqrytiIg6eZmZVG0jKkx0zGAZ+JiDtyJh1bZ98pwBDg28DDLSlgTg6eZmaWiBSK2uu/gS8CPwZelDS66tisiJglaQPgEWBCREwAiIhJSxVXegYYWu9Yuzl4mplZUs7M8Ntm28OzV7VjSLMNVcL4gJ1C1sHTzMySEoJnRAzLcc6MrDSNzhvTfIn6xsHTzMyWcFTIZcA2ic3MzAYqf8cwM7OknAFDXcHB08zMknIGDHUFV5OZmSUOnrn5nqeZmVlB/o5hZmZL+J5nLg6eZmaWuNs2N1eTmZklDp65uZrMzGwJd9vm4gFDZmZmBbnlaWZmibttc3M1mZlZ4uCZm6vJzMwSB8/cXE1mZraEBwzl4gFDZmZmBbnlaWZmibttc3M1mZlZ4uCZm6vJzMwSr+eZW+n3PCWtJ+kySc9Kek7SFZLWz5FulKSJkh6Q9JKkmZIukrRhGeU2MzOrKLXlKWlF4BZgIbAHEMBxwK2S3hcRL/aSfBdgE+BUYCqwDnAkMEXSByLi8bYW3sys27nbNreyq2lfYDiwcUQ8DCDpbuAhYH/gpF7S/iQi5lTvkHQ7MD3Ld3xbSmxmNpg4eOZSdrftDsAdlcAJEBHTgduBHXtLWBs4s32PAXNIrVAzM2tG5Z5nM69BouzguQlwb539U4GRRTOTNAJYA7i/yXKZmVml27aZ1yBRdvBcDZhfZ/88YNUiGUkaCpxBanme03zRzMzM8unk7wmnAZsB20VEvYAMgKT9gP3Su1VKKZiZWUfygKHcyq6m+dRvYfbUIq1L0gmkgLhHRNzY27kRMRGYmNKtHfmLamY2CDl45lJ2NU0l3fesNRK4L08Gkg4Hfgh8OyIuaGHZzMwGN0+SkFvZ9zyvBkZLGl7ZIWkYsHl2rFeSvkN6LvTwiDitXYU0MxuUPGAot7KD51nADOAqSTtK2gG4CngcOLNykqQNJC2SNL5q3y7AKcDvgVskja56FR6pa2Zm1lelfk+IiBcljQNOBi4gfc+5GTgoIl6oOrXSeVAd3LfJ9m+TvapNBsa0qdhmZoODBwzlVno1RcRM4AsNzplB+mes3rcnsGe7ymVmZvieZ07+jmFmZolbnrmVvqqKmZlZp/N3DDMzS9zyzM3VZGZmiYNnbq4mMzNbwgOGcnHwNDOzxC3P3DxgyMzMrCB/xzAzs8Qtz9xcTWZmtoTveebi4GlmZolbnrm5mszMLHHwzM0DhszMzArydwwzM0u8GHZuDp5mZpa42zY3V5OZmS3hqJCL73mamZkV5O8YZmaWuNs2N1eTmZklHjCUm4OnmZklbnnm5moyM7MlHBVy8YAhMzMrjaSdJF0u6TFJCyRNk3S8pDc3SDdK0kRJD0h6SdJMSRdJ2rCsslfzdwwzM0vKued5MDATOAyYBWwKHA2MlbRZRLzWQ7pdgE2AU4GpwDrAkcAUSR+IiMfbXfBqDp5mZpaUc89z+4iYU/V+sqR5wHnAGOCWHtL9pCYdkm4HpgP7AuPbUNYeOXiamVlSQvCsDYCZu7LtOkXSRcRjkub0lq5dHDzNzGyJ/nlUZctse3+RRJJGAGsUTdcKDp5mZtZKq0uaUvV+YkRM7OlkSesAE4CbImJKT+fVSTcUOAOYA5zT18L2lYOnmZklrem2nRsRo3JdTloZuApYBOxV8DqnAZsB20XE/IJpm+bgaWZmSYmTJEhaAbgGGA5sGRGzCqQ9AdgP2CMibmxTEXvl4GlmZklJwVPSssBlwChgq4i4p0Daw4EfAt+OiAvaVMSGHDzNzKw0kpYBLgLGAZ+JiDsKpP0OcBxweESc1qYi5uLgaWZmS7R/tO1/A18Efgy8KGl01bFZETFL0gbAI8CEiJgAIGkX4BTg98AtNemei4j72l7yKg6eZmaWlNNtu222PTx7VTuGNNtQZa6j6ilkt8n2b5O9qk0mTbBQGgdPMzNLypkkYViOc2ZkpanetyewZzvK1BcOnmZmtoTX88zFwdPMzLpeNhnDVsBoYG1gBWAuMI3U7Tu5l0npl+LgaWZmSRcuhi1pS+AHwNakdvUs0qxEC4D3A9uTJpWfLeks4KSIeK5Rvl1WTWZm1mddFjwl/Q4YS5qMYWfgjxExt+acZYD3kILoV4BvStotIm7oLe8uqiYzM2tKlwVP4EFg74h4sqcTsq7au7PXjyXtAKzSKOPuqiYzM2tKdNGAoYj4bh/SXJ3nvGUan2JmZtadJL21L+kcPM3MDIAQLB7a3GugkrSvpB9UvX+vpFnAU5KmSFqrSH4OnmZmlnRx8AS+TRphW3ES8AxwEOke54QimQ3sj2pmZqUJwaIhzbapcj8qWbYNgAcAJK0CbAl8NiKuk/Q0cHyRzNzyNDOzwWAZlkT2jwMBTMrePw6sUSQztzzNzAyAkFg8tNmw8EpLytIGDwHbAbcAuwB/joiXsmNrA/OKZObgaWZmr1s8pIueVXmjnwEXSNoDWJW0LFrFWNJznrk5eJqZGQCBWNylM8NHxP9Imgl8FLgrIm6rOvxvINfznRWFgmc2lLd6Qt3pETFg2+hmZpZfIBZ1UfCUdBhwZUTcDxARfwL+VHteRBxVNO+GA4YkjZJ0hqTpwBPAXcBtwH3As5Juk3SApDcXvbiZmVkbfRW4V9I0ST+V9LFWZdxjy1PSKFIf8RbAPcC1wD9YMhv9asCGpCbwCcAJkn4K/DwiXm5VAc3MrDyLu+huXkSMkPRu4LPAjsD3JT1Fmij+SuDmvvae9lZLk4GzgAMqTd6eSFo+K9ghpNbssX0pjJmZ9Z9uvOcZEQ+wpIG3FksC6ZXAQkk3ZD//Ls9SZBW9Bc939DYTfU3hXgYuAS6RtGbei5uZ2cDRjcGzWhbTzgDOyG41bkcKpKcDK0iaFBFb58mrx+CZN3DWSffvvqQzM7P+183Bs1pEPA9cDFwsaVngE6RAmkufZhiStEztq0Da9SRdJulZSc9JukLS+n0ow6GSQtJSI6fMzMzyiohXI+L3EXFA3jS57gxLWgE4ivRQ6bp10kWevCStSJrdYSGwR5buOOBWSe+LiBdzlmc4cATwVJ7zzcyssW57VKVa1sjbjxTH1gOWrzklImKDvPnlHVb1S2BX0gili+n7/Ev7AsOBjSPiYQBJd5OmTdqfNMt9HqcDFwEb44kezMxaIt3z7No/qT8Fvkd6auQumpxHMG8t7QAcHBGnNnOxLJ87KoETICKmS7qd1NfcMHhK+grwQeDLwBVNlsfMzKp08T3P3YBj+zIhQj1571UuBHp9XCWnTYB76+yfCoxslFjSqsDJwCERUWgSXzMzG9SGkib4aYm8wfNc0iz0zVoNmF9n/zzSRL2NnAg8mJUnF0n7ZauET4GXGicwMxukKo+qNPMawC4Dcj2GkkfebtsjgdMl3QjcQJ0AGBG/alWh6pH0H8DuwAcjIvKmi4iJwMSUx9q505mZDTYBXTtgiHS/8yJJE+k5jt2SN7O8wfNDpPuVawCfrHM8gDzBcz71W5g9tUirnQmcA8yS9JZs31BgSPZ+QUQszFEGMzOrq6sHDL2dNGB1R2Cfqv0BKNvm/uaQt5bOAJ4mjZZ9gL6PUppKuu9ZayRpovnejMheX69zbD7wXeCUPpbLzGzQ6/IZhn4NrA4cSHNxDMgfPN8N7BQR1zVzMdJ6aT+TNDwiHgWQNAzYHDi0QdqxdfadQvqm8G3g4TrHzczMAEYBu0fEZa3ILG/wnAas1ILrnQV8C7hK0hGkZvKxwOOkblkAJG0APAJMiIgJABExqTYzSc8AQ+sdMzOz4rq45TmTJlub1fKOtj0UOCILan2WzSA0jjRi9gLSRAfTgXER8ULVqSK1KPs0faCZmRXX5aNtjwN+KGnlVmSWt+V5BGmw0IOSHmTpwT0REVvmySgiZgJfaHDODFIAbZTXmDzXNDOzxrp5ej7SYyrrAjMk/YX6cWyPvJnlDZ6LSTdYzczMOtHHgdeA54H31Dle6FHGXMHTLTwzs8GhWx9ViYgNW5lfrnuKktZtcDxXl62ZmQ1c3XzPU9IqDY43nCK2Wt4BOTdUTUxQe8H/AK4tclEzMxt4ujl4AtdKelO9A5JGADcXySxv8HwB+J2kN6x/JunjwHWk5zfNzKzDLWJIU68BbC3gYklvGIwq6d2kdaanFsksb/DcDngrcGm2oCiSNiMFzt+RlnoxMzMbqLYGRvPGOQXeRQqc04Dti2SWd8DQXEnbALcD52QT615Pmlx31yITtZuZ2cDUzYthR8Sjkj4NTJL0b+A84FbShDzbRcSCIvnlrqWImCFpW2AysCtwDbBLRCwuckEzMxuYunxuWyLiH5I+T+ox/SZpneptsgl8CukxeEr6Wg+Hrga2BW4E9qh0H7d7STIzM2u/bgqeksbV2R2kNaE/T1oj+qNVcawlS5Kd3SDt6TWFcfA0M+tgXTjD0E0sWXKsovr95dm2pUuStfSBUjMzs5LVW42rJXoMnhHxWLsuamZmA0+3DRiKiMntyrt7asnMzJrWTfc826nH5zwl/VPS52ofKO3l/HUlnSrpkNYVz8zMytJtMwxJulrSpgXOX17S9yR9vdG5vU2ScD5p8epZkk6W9HlJ75D0/yS9SdJakjaTdJCkm4EZwMbAb/MW1MzMrI1mAHdI+quk70j6oKQ39LhKWlvSZyWdA8wG9gb+3ijj3u55npRltk+W2YEsvWSLgIXAVcAn2tm/bGZm7dVtz3lGxHck/QI4CDgaWAUISc+RYtdbgOVIsezO7LwL88xf0Os9z4h4Fvg58HNJ65OmNlobWB54mrTG550RsbAPn8vMzAaYLntUhYh4BPi2pO8DHwM+ytJx7Laig2SLzDA0E5hZJHMzM+sc3TbatlpEvEKaIa8lPaTdWUtmZlZYt3XbtlPeVVXMzMyaJmknSZdLekzSAknTJB0v6c050i4v6URJs7O0f5G0RRnlruWWp5mZva6ElufBpFuAhwGzgE1Jg3nGStosIl7rJe05pCUyfwA8Sprc/QZJH4uIf7az0LUcPM3MDChtbtvtI2JO1fvJkuaRlggbQ1pfcymS3g98BfhaRPw62zeZtIj1BGCHdha6loOnmZkB5QwYqgmcFXdl23V6SboD8CpwSVVeiyRdDBwq6U1lPvnh4GlmZq/rpwFDW2bb+3s5ZxNgekS8VLN/KulZzY2yn+uS9K6IeLCpUlbJHTwlDQd2BtYnPR9TLSJi71YVyszMBgdJ65C6XW+KiCm9nLoaML/O/nlVx3vzgKRbgTOAKyNiUeHCVskVPCV9FvgNaXTuU6SZGarVzjxkZmYdpkWPqqwuqToIToyIifVOlLQyaYa6RcBezV64ga8B+5G6fZ+S9CvgrIiY3pfM8rY8jwUmAbv20F/dEd7ObPbjmP4uhpnZgFEb1VoQPOdGxKhGJ0laAbgGGA5sGRGzGiSZD2xQZ3+lxTmvzrHXRcS5wLmS3gfsD3wDOETSTcDpwDUNRvq+Qd7nPIcDP+vkwGlmZr2rjLZt5pWHpGWBy4BRwKcj4p4cyaYCG0pasWb/SOAV4OFcnzHi7oj4JmmKvv2BNYErgJmSjpa0Zp588gbPB4C35jzXzMysLknLABcB44DPRsQdOZNeAywLfLEqr6HAl4Ab+zDSdhjwvmz7CnAv8D3gYUmfa5Q4b7ftIcApkv4aEY8WLKCZmXWAkua2/W9SAPwx8KKk0VXHZkXELEkbAI8AEyJiAkBE/EPSJaRYtCwwHTgA2BDYNc+FJS2XXXt/YHPgMeAE4JyImCtpVVJP9knAlb3l1WMtSbqtZtdbgfslPcTSfcsREVtiZmYdrYRHVbbNtodnr2rHkGYbEjCEpXtH9yIF3eNIy4n9C9gmIhquvynp58DuwKrADaTnRq+LiNcHvEbE/GwJs9r4t5TevmK8xhtH0U5rlJmZmXWuMiaGj4hhOc6ZQQqgtfsXkLpWv9eHS38V+BVwRoMRtg+QY+Rvb4thjylcNDMz61glTc/XX9bNliXrVUTMJU0V2KtcA4Yk7S6p7oAhSatJ2j1PPmZmZv1kgaSP1Dsg6UOSFhfJLO9o218D7+jh2IbZcTMz63CLGdrUawBbqhu4yhAKTvaT95P2dtGVSLNDmJlZB+vGxbCzR2MqMWyZ7H21FUiDmOYWybe30bYfAD5YtWt7Se+pc9FdgIeKXNTMzAaebgueko4CxmdvA7i9l9N/WSTv3lqeOwJHVV20dkhxxdOAJ4U3M+sC3RQ8SdPKQmp5jictpl07DeBC4D7g2iIZ9xY8TwHOzS76KPB54B91Lvrv6udkzMzMBoKImAxMBpAUwNkR8UQr8u7tUZVngWezi24IzM4zzNfMzDpTNz+qEhEtXRUk14ChiHgMQNJY4GOk1b6fAP4SEbe2skBmZtY/SpqerzTZsmPHRsT07OfeFFqXOu96nqsBlwJjSTMPzSdNcaRscdGdI6LX5WDMzGzg67J7nmOBX2Q/j6P3x1Ha8qjKqcCHgd2ASyPi1Wxi3p1JI5R+QZr6yMzMbECIiA2rfh7WyrzzBs/tgR9FxP9UFeRV4KKsVXpcKwtlZmbl67ZHVdopb/BcTM/Pck7LjpuZWQfr5gFDkvYCNoiIo+scOxqYHhEN57StyDs931WkBUfr2QX4bd4LmpnZwNXF0/MdSJqXoJ6ngIOKZJb3k14DnCzpd6SBQ/8G1iTd89wEOFDSuMrJEXFLkUKYmVn/6/Ju242AqT0cu5+e52+vK2/wvCzbrseShUyrXZ5tRRqx1LW1b2ZmHWkRsHoPx95WNLO8wXNs0YzNzKyzdHnL807g68Bv6hz7OnBXkczyTpIwuUimZmbWmbp1wBDwY+AmSX8FziZN9LMOsA9pEZStimRW6O6upNWB0cBbgWsiYp6k5YFXIuK1InmZmdnA0m0zDFWLiMmSdiLN235m1aEZwBciYlKR/PLOMCTgp8C3geVI9zU/DMwjjcT9E3BskQubmdnA0uXdtkTEVcBVkjYmNQLnRsSDfckr76MqPwK+BUwAPsobF8e+BvhM3gtKWk/SZZKelfScpCskrV8g/QhJl0qaK2mBpGmSDsyb3szMBreImBYRf+5r4IT83bb7ABMi4nhJtV9LHibnEF9JKwK3kJYy24PUgj0OuFXS+yLixQbpR2XpJ2VlehZ4J7Byzs9hZma96OaWp6T3ktap3pI0P/t84FbS5PH3FMkrb/BcB7ijh2OvACvlzGdfYDiwcUQ8DCDpbtLsRfsDJ/WUUNIywPnAzRHxuapDXtXFzKwFurnbVtKHSWt7LgCuBp4E1iJNP7udpC0i4m9588sbPJ8A3kP9QPV+YHrOfHYA7qgEToBsqZjbgR3pJXgCY4ARpCBrZmYtFnT1aNvjgXuBT0TE85Wdkt4M3JQd/1TezPLe87wUGC9p86p9IeldwPeBi3Pmswmp8LWmAiMbpP14tl1e0h2SXpX0lKRTJa2Q8/pmZjY4jQaOrw6cANn7n5DWqs4tb8vzaGAz4DbgsWzfpaQZh/4MnJAzn9VIfcy15pH6n3uzdra9BDgNOBQYRRrEtB7wuXqJJO0H7AewSs5CmpkNTt37qAqN1+ts/XqeEbFA0hjgK8DWpEFCT5MeT7koIhYVuWgfVVrJF0bE+OznSdkAphMkjYiI+2sTRcREYCLA2lKhyjEzG0y6+Z4n8FfgMEk31XTbrgT8kJ7H9dSV+ytGRCwGLshefTWf+i3Mnlqk1Sqz4f+hZv+NpJbvpqTJfc3MrI+6OHgeRnpS4zFJ1wKzSQOGPg2sSBpXk1veSRKWJ3WRvp3UtJ0N/C0iXi5yMdK9zU3q7B8J3JcjbW88w5GZWRO6eT3PiLhT0mhgPKkHdTXSLcPWP6oi6U2kmYX2Bd7EkskRAnhZ0unAYRHxSs7rXQ38TNLwiHg0u8YwYHPSPczeXE96PnRr0sQMFdtk2yk5y2BmZoNQRNwN7NSKvBq1PK8FxpGm4LsOmEkKoOuRZhX6LqnV+Omc1zuLNFPRVZKOIAXhY4HHqZprUNIGwCOkiRkmAETE05KOB46U9BxpsoRRpG8R51U//mJmZsV189y2rdZjLUn6Imkpsp0i4so6p5wt6QvAJZI+HxFXNLpYRLyYLZp9MuneqYCbgYMi4oXqy5PWBK19lGYC8DzwDeBgUvfxiXheXTOzluime56SflXg9IiIvfOe3NtXjC8Dv+khcFaudLmkS4FdgYbBM0szE/hCg3Nm8Mb5cyv7gzSRQm+TKZiZWR904WjbceR/BKVlj6psChyRI49rSfPTmplZBwvE4te6J3hGxLB25d3bDENvI93jbGQmsEZrimNmZjbw9dbyXJE0urWRV4DlW1McMzPrNwGLFnVPy7NWNiHC3sAWpPU894uIhyTtAvwzIh7Im1ejYVXrSBre4Jx1817MzMwGrgixeFF3jraVtB5pkoR1gQdIi528OTs8FvgkaanLXBrV0mV5ykTBG61mZjbwpODZtS3Pn5N6U99FWimsen6CyaR1PnPrLXjuVbhoZmZmA9NWpG7ax7I50as9QVq3Orceg2dEnNeHwpmZWacKurnluRxpnoB6VgEKLXDSnZ3bZmZWWIRY9GrXBs+7SXMM/L7OsW2BvxXJzMHTzMwy4rXFXRsWTgQukwTwP9m+kZJ2JI3A3aFIZl1bS2ZmVlAAXdptGxFXSPoGaQnLr2W7zyd15X4rIuq1SHvk4GlmZl0pm9v23Ii4DSAizpB0AfAx0uQ+TwN/rl4cOy8HTzMzS0Ld1vL8ErCHpJmkVub5EfEIcFOzGfc2PZ+ZmQ0mASxSc6+BZU3SxAczSHO1Pyjpdkn7SlqlmYwdPM3MbIlFTb4GkIh4ISJ+HRFjgWHAkcCqpPWjZ0u6WNK2kgrHQgdPMzPrehHxeET8Z0SMBEYDvyItWXYt8ISknxXJz8HTzMySoKtanj2JiDsj4lukWYVOJg0e+m6RPDxgyMzMkkrw7HKSNgJ2B3Yjdec+B/ymSB4OnmZmlgTwan8Xoj0krQrsQgqaHyF92j8AhwG/jYiXi+Tn4GlmZkkAi/u7EK0jaVngM6SAuS1pftv7gEOBCyNidl/zdvA0M7Nu9W/SpO/zgInAeRFRaA7bnnjAkJmZLVHCgCFJ60r6L0l/kfSSpJA0LGfat0r6haRHJS2QNF3SaZLeVuf0yaTJ4NeOiO+0KnCCW55mZlZR3oChjYCdSSuZ/BH4VJ5ESrO6X01a0Ho8cD8wEpgAjJL0sYiIyvkR8bkWl/t1Dp5mZpaUFzxvi4g1ASTtQ87gCbwT2AzYPyImZvsmSXoNOJ0UVKe1urD1OHiamVlSUvCMiNf6mHS5bPtczf5nsm1ptyIdPM3MrFNMBW4DjpT0MPAAqdt2PHB9RNxfVkEcPM3MLGlNy3N1SVOq3k+s6mJtSkSEpE8DFwB3VR36HfDFVlwjLwdPMzNbovngOTciRrWgJD05izQ37ddJA4ZGAMcAl0navoku4UIcPM3MLBngMwxJ2g74MvDJiLg5232bpEeBG4HtgavKKIuf8zQzs07x3mx7V83+O7PtiLIK4uBpZmZJZXq+Zl7t9WS2/UjN/o9m2yfaXoKMu23NzCwpcVUVSTtlP34o224raQ4wJyImZ+csIk2pt3d2zhXAj4HzJR1LGm37buAo4HHgynJK7+BpZmYV5S5JdmnN+19m28nAmOznIdkLgIh4TtJo4GjgEODtwGzgGuDoiHihjeV9AwdPMzNLSgyeEaG+nBMRjwN71zm9VL7naWZmVpBbnmZmtkR53bYdzcHTzMyScu95djQHTzMzSxw8c3PwNDOzZIDPMDSQeMCQmZlZQW55mplZUplhyBpy8DQzsyV8zzMXB08zM0s8YCg33/M0MzMryC1PMzNL3PLMzcHTzMwSP6qSm4OnmZklHm2bm4OnmZkt4W7bXDxgyMzMrCC3PM3MLPGAodwcPM3MLPGAodwcPM3MLPGAodwcPM3MLHG3bW6lDxiStJ6kyyQ9K+k5SVdIWj9n2vUlnSdppqQFkh6UdJykldpdbjMzs4pSW56SVgRuARYCe5C+5xwH3CrpfRHxYi9pVwJuApYFjgRmAh8GjgHeCXypvaU3MxsE3PLMpexu232B4cDGEfEwgKS7gYeA/YGTekm7OSlIbh0RN2b7bpW0GnCwpBUj4qX2Fd3MrMt5wFBuZXfb7gDcUQmcABExHbgd2LFB2uWy7XM1+58hfQ61qIxmZoNTZcBQM69BouzguQlwb539U4GRDdLeRGqh/kTSSEkrSxoHHAic0VuXr5mZWSuV3W27GjC/zv55wKq9JYyIlyV9HLicFGwrzga+1bISmpkNVh5tm1vHPKoiaXngEmAN4KukAUMfAcaT/rkP6CHdfsB+AKuUUlIzsw7l4Jlb2cFzPvVbmD21SKvtDYwBNoqIR7J9t0l6Fpgo6YyI+FdtooiYCEwEWFuKvhbczKzrecBQbmXf85xKuu9ZayRwX4O07wXmVwXOijuz7Ygmy2ZmZh4wlEvZwfNqYLSk4ZUdkoaRHkO5ukHaJ4FVJW1Us/+j2faJVhXSzMysN2UHz7OAGcBVknaUtANwFfA4cGblJEkbSFokaXxV2nOB54HrJO0haaykHwA/A/5GetzFzMz6qnLPs5nXIFFq8MweJxkHPAhcAFwETAfGRcQLVacKGFJdvoiYAYwG/kmaleg60qQLE4GtIuK19n8CM7Mu5uCZW+mjbSNiJvCFBufMoM6kBxFxH7Bze0pmZjbIecBQbh3zqIqZmbWZlyTLrfRVVczMzDqdW55mZrbEILpv2QwHTzMzSzzDUG4OnmZmlnjAUG6+52lmZlaQW55mZpZ4tG1uDp5mZpb4nmduDp5mZraEg2cuDp5mZpZ4wFBuHjBkZmZWkFueZmaWeMBQbg6eZmaWeMBQbg6eZmaWOHjm5uBpZmaJBwzl5gFDZmZmBbnlaWZmS3jAUC4OnmZmtkT0dwE6g7ttzczMCnLwNDMzK8jB08zMSiVpXUn/Jekvkl6SFJKGFUi/jqRfSXpS0kJJ0yUd38YiL8X3PM3MrGwbATsDfwP+CHwqb8IsyN4OTAe+A/wbGJblWRoHTzMzy5T2oOdtEbEmgKR9KBA8gTOAJ4CxEVEp7OQWl68hB08zM8uUM8VQRLzWl3SS3gFsDexeFTj7he95mplZp9g82y6Q9Ifsfud8SedLemuZBXHwNDOzTKXbtpkXq0uaUvXar4UFXDvb/gp4ENgW+CGwHXCDpNJimrttzcws05Ju27kRMaoFhamnEhwnRcQ3s59vkfQscDGpS/f6Nl37DRw8zcwsM+Bnhn862/6hZv+N2XZTHDzNzKxcAz54Tm1wvE8DkfrC9zzNzKxT3AE8SeqerbZNtr2rrIK45WlmZlXKWQ1b0k7Zjx/KtttKmgPMiYjJ2TmLgPMiYm+AiFgk6VDgXElnAFeQJkf4MTAJuKWUwuPgaWZmryu12/bSmve/zLaTgTHZz0Oy1+si4jxJr5FG2e4FzAMuBH4UEaWtCePgaWZmmXImSQCICPX1nIi4ALig5YUqwPc8zczMCnLL08zMMgN+tO2A4eBpZmaZ8rptO52Dp5mZZdzyzMvB08zMMm555uUBQ2ZmZgW55WlmZhl32+bl4GlmZhl32+bl4GlmZhm3PPPyPU8zM7OC3PI0M7Mq7rbNw8HTzMwy7rbNy8HTzMwyDp55OXiamVnGo23z8oAhMzOzgtzyNDOzjLtt83LwNDOzjLtt83LwNDOzjFueeTl4mplZxi3PvEofMCRpXUn/Jekvkl6SFJKG5Uy7jKQfSZoh6WVJ/5L0hTYX2czM7A36Y7TtRsDOwHzgjwXTHgscDZwGbAvcAVwq6dOtLKCZ2eBU6bZt5jU49Ee37W0RsSaApH2AT+VJJGkN4GDghIj4Wbb7VkkbAScA17WjsGZmg4e7bfMqveUZEa/1MenWwHLAhTX7LwTeK2nDpgpmZjboueWZVydNkrAJsBB4uGb/1Gw7stzimJnZYNVJo21XA56JiKjZP6/quJmZ9Zm7bfPqpODZJ5L2A/bL3i48Bu7tz/IMAKsDc/u7EP3MdeA6qHA9wMZLfpx9Axy9epP5DYr67KTgOR94iyTVtD4rLc55ddIQEROBiQCSpkTEqPYWc2BzHbgOwHVQ4XpIdVD5OSK26c+ydJJOuuc5FXgT8I6a/ZV7nfeVWxwzMxusOil4/p40lGvXmv27AfdGxPTyi2RmZoNRv3TbStop+/FD2XZbSXOAORExOTtnEXBeROwNEBFPSToJ+JGk54G/A18CxgE75Lz0xFZ9hg7mOnAdgOugwvXgOugTLT14tYSLSj1ddHJEjKk657yI2LMq3RDgR8C+wFrANGBCRFzW1gKbmZlV6ZfgaWZm1sk66Z5nXZLWk3SZpGclPSfpCknr50y7vKQTJc2WtCCbrH6Ldpe51fpaB5JGSZoo6YFskv6Zki7qxNmamvk9qMnn0Gyxgj+1o5zt1mw9SBoh6VJJc7P/E9MkHdjOMrdak38T1pd0XvZ/YYGkByUdJ2mldpe7lbwAR/t1dPCUtCJwC/BuYA/gq8A7SXPe5vllP4fUBTwe+AwwG7hB0gfaUuA2aLIOdiHN3HQqaaL9Q4EPAlMkrde2QrdYC34PKvkMB44AnmpHOdut2XqQNAr4K2lU+z7Ap4GfA0PaVeZWa6YOsuM3AVsAR5I+/9nA94FftbHY7eAFONotIjr2BRwILAY2qtq3IWmKjO81SPt+0nQae1XtG0q6j3p1f3+2kurgbXX2bQC8RrqX3O+fr911UJPPDcCZwCTgT/39uUr+XViG9LjXlf39OfqxDj6V/U34VM3+E7L0K/b35ytQD8tU/bxP9rmG5Ui3Bmka1GNq9t8M3N3fn2sgvTq65UkaZXtHRLw+322kR1ZuB3bMkfZV4JKqtIuAi4GtJb2p9cVtiz7XQUTMqbPvMWAOsE6Ly9lOzfweACDpK6RW94/aUsJyNFMPY4ARwEltK105mqmD5bLtczX7nyF9uVCLyth24QU42q7Tg+cm1J9ubyqNJ4rfBJgeES/VSbscqdujEzRTB0uRNIL07fP+JstVpqbqQNKqwMnAIRFRd6aqDtFMPXw82y4v6Q5Jr0p6StKpklZoaSnbq5k6uAl4CPiJpJGSVpY0jtSaPSMiXmxtUQckL8CRU6cHz9VIffq15gGrNpG2crwTNFMHbyBpKHAGqeV5TvNFK02zdXAi8CBwbgvL1B+aqYe1s+0lwI3AVsBPSV1+/9OqApagz3UQES+TvkQsQwoWz5O6K68FvtXaYg5YXoAjp06a29ba7zRgM2C7iKj3B6jrSPoPYHfgg3X+YAwmlS/SF0bE+OznSdmz1SdIGhERndQbUZik5UlfHtYgDTSaCXyENKBwEXBA/5XOBppOD57zqf9tsqdvn7VpN+ghLfQw0fwA1EwdvE7SCaTVZ/aIiBtbVLayNFMHZ5Ja2bMkvSXbNxQYkr1fEBELW1TOdmumHp7Otn+o2X8jacDMpnRGV34zdbA36d7vRhHxSLbvNknPAhMlnRER/2pZSQemPi3AMRh1erftVFIffa2RNJ4ofiqwYTa0vTbtKyzd5z9QNVMHAEg6HPgh8J2IuKCFZStLM3UwAvg66Y9G5bU5MDr7uZNaG83+f+hNXweglK2ZOngvML8qcFbcmW1HNFm2TuAFOHLq9OB5NTA6ez4PgOxB4M2zY725BlgW+GJV2qGk+XJv7KDWRjN1gKTvAMcBh0fEae0qZJs1Uwdj67z+RRp0MhbopKkfm6mH60kDRbau2V9ZomoKnaGZOngSWFVS7WDBj2bbJ1pVyAHMC3Dk1d/PyjTzAlYitRDvIQ1D34H0h+9RYOWq8zYg3bMYX5P+YlLrYh/gE6Q/lC+T7n/1++drdx2QJkl4jfSHc3TNa2R/f7ayfg/q5DeJznzOs9n/D0dl+/8T+CRp0owFwLn9/dnKqANgGOkxlQdJEyyMBX6Q7ZtC1bOTnfACdspep5Oe8zwge79l1TmLgHNq0p2Q/R38Hqkb+/Ts78Rn+vszDaRXvxegBb8g6wOXZ7/gzwO/peZh4Ow/RQBH1+xfgfRc25PZL8tfgTH9/ZnKqgPS6NLo4TWpvz9XWb8HdfLqyODZbD2QnmP8XhZ8XgEeAyYAy/b35yqxDkYCvwEeJ31xeBD4GbBqf3+uPtRDw//b2ftza9INIc209RipN+JuYKf+/jwD7eWJ4c3MzArq9HueZmZmpXPwNDMzK8jB08zMrCAHTzMzs4IcPM3MzApy8DQzMyvIwdMGBEmXSJonaa2a/UMk3SXpoYG0NJakYZJC0p5V+/aU9LU65+6ZnTuszDJm115G0j8lHVy17+isPG2b21rSQZLukeS/MdaV/IttA8W3SQ9s/7Jm/8HAh4B9ImJB6aXq2WzgY8DvqvbtCSwVPLNzPpalKdtuwNtZul7b7UzgbaSZesy6joOnDQgR8RTwXeBzkr4IIOldwNHAmRExuR+Lt5SIWBgRd0TEnBznzsnO7Y/5kg8Gzo+lF31vq+yLzvnZ9c26joOnDRgRcT5pYurTJK1OWipsDnBIo7RVXaNbSPqtpBckPS3pv2u7eyW9XdL5kuZKWijpbkm71ZyzlqTzJP1fds5sSddKWiM7/oZuW0mTgC2BzbP9ke2r220raVlJx0maIemVbHucpGWrzqlcY39JE7IyPCPpGknr5qiTj5JWCmm4mLWkbbI6Oy3r6q1c++uSjpf0pKTnJV0oaUVJG0m6IUvzsKR6LcyLgZGSNmt0fbNO0+nreVr32Z+0LNJfgeGkhbmfL5D+QtLcpL9kyULGK5G6VJG0EjCZtObjYaQ5THcDLpC0YkRMzPK5gDR5+A+yc9YkLR5Qu4RdxTeyaw/JPgOkuVV7ch6wM2kS9j+RFiE/PPvMX6k590fAn0ldwmsAP8+uNaaX/CGtiPI8aWL0HknaHTgbmBARx2X7qq89idT9OhL4KWmS8E2Bs0jzvh4A/FrSlIioXtrsn9n1t8nKb9Y9+ntyXb/8qn0Bx5Puf15eIM2eWZozavYfDiwG3pW9/1Z23pia824CngKGZO9fIK1v2tP1hmX57Fm1bxJ1JpSvKtuw7P17qD8p+RHZ/vfVXGNSzXkHZ/vXblAn1wO319l/dJZ+KKlV/yrpnnK9z3dLzf4rsv27Ve1blbQ6x1F1rvVH0hJ//f575ZdfrXy529YGFEn/D/gq6Q/0hyW9uWAWv6l5fzHp9sRHsvdbAE9ExKSa8y4kDXCpLPp7F/ADSQdKeq+qmmItsEXVNWvLAKn7t9p1Ne/vybbrN7jO2qRu756cDBxDWjHj7B7Oub7m/QPZ9obKjoiYT/risV6d9HOycph1FQdPG2hOJLVktiN1UR5fMP2/e3i/TrZdjfqjXp+sOg5pUfSrSS2zu4EnJI1v0aMXlWvUlqO2DBXzat5XBh4t3+A6y1edW8+XSYt+39TLOfNr3r/Sy/565VlAWvrPrKs4eNqAIWkMsC9wRERcDxwHHFBwwMmaPbx/ItvOA9ZiaWtVHScinoqIb0bEOsC7SWufHsOS+5nNqATD2nKsVXO8WU+Tvoj05BOk1uv1klZu0TVrrQbMbVPeZv3GwdMGhGxE7Fmk7tJfZLt/Qho8dLak5XJmtXPN+11IA1z+mr2fDKwrafOa875C6nq8rzbDiJgWEYeRWlvv6eXaC8nXyrqtqmzVds22k3LkkccDpAFIPZlKGnT0TtoXQDcEprUhX7N+5eBpA8UE0ujWfSLiNYCIeBXYB9iYNPAnj09LOlHSVpIOB44iPef4UHb8XOAh4ApJ+2SPaFwAbAUcGRGLJa2SzWp0UHb8E5JOJbXibuzl2vcB75H0JUmjJG1c76SIuBf4X+BoSUdlZR1PGsjzvxFxT710fXAb8A5Jb+3phIi4nxRA3wHc0Id7zD2S9BbgXSz5smDWNfyoivU7SaNIEyT8Z23giIg7Jf0COFTSb+KNj0LUsxvwfdLjE6+QWrOvP6gfES9K2pL0yMUJwJtJLaOvRkRlwM7LwN9JXcgbkFqu04BdI+KqXq79E1KgPxtYmdTKHdPDuXsCj5IePzkC+L8s/TENPl8RV5E+y2dIj8bUFRHTsjq5FbhR0tYtuv52pH+DK1uUn9mAoYjo7zKYNS2brODXwDsj4uF+Ls6AIelcYN2I+GQ/XPt6YG5EfLXsa5u1m1ueZt3tGOB+SaMiYkpZF5X0AWAcsElZ1zQrk+95mnWxiJhO6iJeo+RLr0WaQMK9ANaV3G1rZmZWkFueZmZmBTl4mpmZFeTgaWZmVpCDp5mZWUEOnmZmZgU5eJqZmRX0/wFYazX2t7L1GQAAAABJRU5ErkJggg==\n" }, "metadata": { "needs_background": "light" } } ], "source": [ "#NBVAL_IGNORE_OUTPUT\n", "from examples.seismic import Model, plot_velocity\n", "\n", "# Define a physical size\n", "shape = (101, 101) # Number of grid point (nx, nz)\n", "spacing = (10., 10.) # Grid spacing in m. The domain size is now 1km by 1km\n", "origin = (0., 0.) # What is the location of the top left corner. This is necessary to define\n", "# the absolute location of the source and receivers\n", "\n", "# Define a velocity profile. The velocity is in km/s\n", "v = np.empty(shape, dtype=np.float32)\n", "v[:, :51] = 1.5\n", "v[:, 51:] = 2.5\n", "\n", "# With the velocity and model size defined, we can create the seismic model that\n", "# encapsulates this properties. We also define the size of the absorbing layer as 10 grid points\n", "model = Model(vp=v, origin=origin, shape=shape, spacing=spacing,\n", " space_order=2, nbl=10, bcs=\"damp\")\n", "\n", "plot_velocity(model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Acquisition geometry\n", "\n", "To fully define our problem setup we also need to define the source that injects the wave to model and the set of receiver locations at which to sample the wavefield. The source time signature will be modelled using a Ricker wavelet defined as\n", "\n", "\\begin{equation}\n", " q(t) = (1-2\\pi^2 f_0^2 (t - \\frac{1}{f_0})^2 )e^{- \\pi^2 f_0^2 (t - \\frac{1}{f_0})}\n", "\\end{equation}\n", "\n", "To fully define the source signature we first need to define the time duration for our model and the timestep size, which is dictated by the CFL condition and our grid spacing. Luckily, our `Model` utility provides us with the critical timestep size, so we can fully discretize our model time axis as an array:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "from examples.seismic import TimeAxis\n", "\n", "t0 = 0. # Simulation starts a t=0\n", "tn = 1000. # Simulation last 1 second (1000 ms)\n", "dt = model.critical_dt # Time step from model grid spacing\n", "\n", "time_range = TimeAxis(start=t0, stop=tn, step=dt)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The source is positioned at a $20m$ depth and at the middle of the $x$ axis ($x_{src}=500m$), with a peak wavelet frequency of $10Hz$." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "output_type": "display_data", "data": { "text/plain": "
", "image/svg+xml": "\n\n\n \n \n \n \n 2021-04-15T19:25:24.541945\n image/svg+xml\n \n \n Matplotlib v3.4.1, https://matplotlib.org/\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAg0AAAF9CAYAAACH0lvIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAA9HUlEQVR4nO3deZxkdXnv8c/TXb13zz7DNhvDoDAIbqMSRUHUiGjAXNyiRtEIGo3R5GWiXHONIjG4BI3BBUyuKJhoNBpINIpsGlGujiKEYR2YFYZhmO6ZXqu6luf+cU5199T0cqr7VHWdPt/361Wv6j51ll/9GLqeen7P73fM3RERERGZSdN8N0BERESSQUGDiIiIRKKgQURERCJR0CAiIiKRKGgQERGRSBQ0iIiISCSZ+W5Ao1uxYoWvX79+vpshIiJSF7/+9a+fdPeVk72moGEG69evZ8uWLfPdDBERkbows51TvabhCREREYlEQYOIiIhEoqBBREREIlHQICIiIpEoaBAREZFIFDSIiIhIJAoaREREJJJEBA1mttrM/sHMfmFmw2bmZrY+4rFNZnaJme0ws6yZ3WVmF9S4ySIiIgtOIoIGYCPwOqAP+O8qj/048FHgSuAVwB3At83s3DgbKCIistAlZUXIn7r7UQBm9g7gd6McZGargA8Al7v7Z8LNt5rZRuBy4Ae1aKyIiMhClIhMg7uXZnnoy4FW4LqK7dcBp5rZ8XNqmIiISIokImiYg1OAHLCtYvvW8HlTfZvT+B4/lKVvaHS+myEiIg1ooQcNy4CD7u4V23snvH4EM7vYzLaY2Zb9+/fXtIGNpFAs8b++eDv/+3v/M99NERGRBrTQg4ZZcfer3X2zu29euXLSu4MuSDfdt4/HDmX59c6++W6KiIg0oIUeNPQBS8zMKraXMwy9yJhv/mo3AE8M5Hj8UHaeWyMiIo1moQcNW4E24ISK7eVahnvr25zGdv/eAU5Y2QXA3XsOzm9jRESk4Sz0oOGHQB54U8X2NwP3uPv2+jepMRVLzv7BHGeftIpMk3H3nkPz3SQREWkwSVmnATN7Tfjjs8PnV5jZfmC/u/8k3KcAfM3d/wjA3Z8wsyuAS8xsAPgN8HrgbOC8ur6BBvfkYI5iyVm7vIvjV3Tx4L6B+W6SiIg0mMQEDcC3K37/Yvj8E+Cs8Ofm8DHRh4FB4H3A0cADwOvc/T9r08xk2tcf1DAcvaidlT1t9GrapYiIVEhM0ODulcWMkfZx9yJwWfiQKZQLH49e1M6yrla2PtY/zy0SEZFGs9BrGiSicqbhqEVtrOhu48Bgbp5bJCIijUZBgwCwrz9Hc5OxvLuNZV2t9GcLjBZmu3q3iIgsRAoaBIDH+7Os6mkLA4dWANU1iIjIYRQ0CBAMT6xa1A7A8q4gaDgwpCEKEREZp6BBANg/kGNVTxsAy7uD5wODyjSIiMg4BQ0CQP9IniUdLQAsU6ZBREQmoaBBAOjPFlgUBg0rupRpEBGRIyloEArFEoO5Aj3twbIdizoyZJqMAyqEFBGRCRQ0CIO5AgCL2oNMg5mxrKuVXmUaRERkAgUNwkA2DBrC4QmApZ2t9A0raBARkXEKGoRDI3kAFrWPryre3Z5haLQwX00SEZEGpKBB6M8GQUNP+3imobstw2BWQYOIiIxT0CD0j5SHJw7PNAzkFDSIiMg4BQ0ylmlYNCHT0KNMg4iIVFDQIPSXaxo6KoYnlGkQEZEJFDQI/WFGobvt8OGJ4dEihaLudCkiIgEFDcJANk9PW4bmJhvbVg4ghnLF+WqWiIg0GAUNQv9I4bChCWBsdciBXH4+miQiIg1IQYPQn82PBQll3W1BEKG6BhERKVPQIPSP5I/INHSHQYRmUIiISJmCBmEgW6Cn7fBMw/jwhIIGEREJKGgQhkcLdFUGDW3KNIiIyOEUNAhDo0W62poP2zY2PKFMg4iIhBQ0CMO5Ah0tlYWQyjSIiMjhFDSknLsznD8y09DVqpoGERE5nIKGlMvmS7hDZ+vhmYamJtOdLkVE5DAKGlJuaDQICiozDVC+/4QWdxIRkYCChpQbDpeJ7mg5MmjoamvWMtIiIjImEUGDma0xs++Y2SEz6zez75rZ2ojHrjWzr5nZLjMbMbMHzewyM+uqdbuTYDhfzjRkjnitszXD8KiGJ0REJHDkJ0WDMbNO4BYgB7wVcOAy4FYzO83dh6Y5tgu4CWgB/g+wC3gO8DHgROD1tW194ytnEjpbj8w0dLQ0MzyqTIOIiAQaPmgALgI2AE91920AZnY38BDwTuCKaY59AUFw8HJ3vzHcdquZLQM+YGad7j5cu6Y3vnImobIQEqCjtZm+4dF6N0lERBpUEoYnzgPuKAcMAO6+HbgdOH+GY1vD5/6K7QcJ3ruRcuVMwmSZhs5WZRpERGRcEoKGU4B7Jtm+Fdg0w7E3EWQkPmlmm8ys28zOBt4HfHm6oY20GB6duqaho7WZEQUNIiISSkLQsAzom2R7L7B0ugPdPQucQfA+twIDwM3AfwJ/MtVxZnaxmW0xsy379++fbbsToVzT0DVlpkGFkCIiEkhC0DBrZtYOfAtYBfwhcCbwFwQFkF+Y6jh3v9rdN7v75pUrV9alrfOlnEnomDRoyGh4QkRExiShELKPyTMKU2UgJvoj4Cxgo7s/HG77qZkdAq42sy+7+12xtTSBhqYrhGxpJlcoUSo5TU2pL/8QEUm9JGQathLUNVTaBNw7w7GnAn0TAoayX4bPJ8+xbYk3PFqkvaWJ5kmCgnJx5Ehe2QYREUlG0HADcLqZbShvMLP1BNMpb5jh2MeBpWa2sWL788LnR+NqZFIN5QpjN6eqVA4aNEQhIiKQjKDhK8AO4HozO9/MzgOuB3YDV5V3MrN1ZlYws49MOPYaguLHH5jZW83sxWb2F8BngF8TTNtMtZHR4qT1DAAdYTChGRQiIgIJCBrCaZFnAw8C1wLfALYDZ7v74IRdDWhmwnty9x3A6cBvCVaR/AHBYlFXAy9z91Lt30FjGxqdOtNQvh9FealpERFJtyQUQuLuu4ALZthnB5Ms1uTu9wKvq03Lkm94tEjnJHe4BA1PiIjI4Ro+0yC1NTxanHQ1SBifhqnhCRERAQUNqTc8WqSjRYWQIiIyMwUNKZfNT10IOR40qKZBREQUNKTeyGiRjpbJ/xlo9oSIiEykoCHlRvLFsVkSlTpbtLiTiIiMU9CQctl8kfYpgoYO1TSIiMgEChpSrFRycoXSlEFDW6aJJtPwhIiIBBQ0pFi2MPUdLgHMTHe6FBGRMQoaUiybDxbEnKqmAaC9pZkRrQgpIiIoaEi1coFj+xSzJyCYdqlMg4iIgIKGVCvXKkxV0xC81kRWsydERAQFDalWDgZmGp4oD2OIiEi6KWhIsbGgYYpCSID2TDO5gjINIiKioCHVxmsapg4a2lqalGkQERFAQUOqlWsaZh6eUKZBREQUNKRalExDe0szuYIyDSIioqAh1XLldRqmrWnQ7AkREQkoaEixsUxDZup/BhqeEBGRMgUNKTYSZfaECiFFRCSkoCHFxhZ3ykwzeyLTTLZQxN3r1SwREWlQChpSLFso0pppoqnJptynvaUJdxgtKtsgIpJ2ChpSLDtanHa6JYzPrNAQhYiIKGhIsZH8zEFDW/h6TsWQIiKpp6AhxbL50rRFkDA+s0KZBhERUdCQYiP5Im3TTLeE8eEJ3X9CREQUNKRYNl+cOdOgmgYREQkpaEixkUiFkOHwhDINIiKpp6AhxbKF4rT3nYCJmQYFDSIiaZeYoMHM1pjZd8zskJn1m9l3zWxtFcefbGbfNrMnzWzEzB4ws/fVss2NLpsvjWUSplJe+EnDEyIikpnvBkRhZp3ALUAOeCvgwGXArWZ2mrsPzXD85vD424B3AIeAE4HuGja74eUKRdqmWQ0SJgxPKNMgIpJ6iQgagIuADcBT3X0bgJndDTwEvBO4YqoDzawJ+Dpws7v//oSXbq1dc5MhFyHT0JbR8ISIiASSMjxxHnBHOWAAcPftwO3A+TMcexZwMtMEFmmVzVeRaShoeEJEJO2SEjScAtwzyfatwKYZjj0jfG43szvMLG9mT5jZ582sI9ZWJkyuUJpxnQatCCkiImVJCRqWAX2TbO8Fls5w7LHh87eAG4GXAZ8iqG3458kOMLOLzWyLmW3Zv3//7Frc4Nw9CBqiTrlU0CAiknpJqWmYi3JgdJ27fyT8+TYzawYuN7OT3f2+iQe4+9XA1QCbN29ekPeEzoXDDTNlGlqbmzDT7AkREUlOpqGPyTMKU2UgJjoQPv+4YvuN4fMz59CuxMqFQcBM6zSYGe2ZZi0jLSIiiQkathLUNVTaBNwb4djppPIrdDkImCnTAMEQhTINIiKSlKDhBuB0M9tQ3mBm64EXhK9N578I1nd4ecX2c8LnLTG1MVGyETMN5X1U0yAiIkkJGr4C7ACuN7Pzzew84HpgN3BVeSczW2dmBTMr1y7g7geAvwXeZWafMLOXmtmHgI8AX5s4jTNNqss0NGvKpYiIJKMQ0t2HzOxs4LPAtYABNwPvd/fBCbsa0MyRwdClwADwbuADwF7g08DHa9z0hlVNpqEt06RMg4iIJCNoAHD3XcAFM+yzgyBwqNzuBIs7aYGnUNWZBgUNIiKpl5ThCYlZecpltJqGprHZFiIikl4KGlKqnDmIkmloyzST1ZRLEZHUU9CQUmOLO81wwyooT7lU0CAiknYKGlKqHAS0z3DDKijXNGh4QkQk7RQ0pFRVmYaMCiFFRERBQ2pVl2loGgsyREQkvRQ0pFR1NQ3KNIiIiIKG1BqfPRFhcaeWZnKFEsFyFyIiklYKGlIqVyjR0mw0Nx2xFtYR2sNshIYoRETSTUFDSuXypUj1DDBe96AhChGRdFPQkFLZQjFSPQOMrxqpaZciIummoCGlcvlSpHoGGB+eUKZBRCTdFDSkVDWZhnJwoaWkRUTSTUFDSs0u06DhCRGRNFPQkFK5QnEsGJjJeE2DMg0iImmmoCGlgkxD1KBBNQ0iIqKgIbWyheJYBmEmYzUNGp4QEUk1BQ0pVV2mIQgaciqEFBFJNQUNKZWrItMwtiKkMg0iIqmmoCGlsrPINGjKpYhIuiloSKnqMg2aPSEiIlUEDRY4z8w+Y2ZfNbN14fYzzezY2jVRaqGqTENG6zSIiAhkouxkZkuBHwDPAwaAbuAfgJ3ARUAv8Kc1aqPEzN3JFYqRF3fKNDeRaTJlGkREUi5qpuHTwBrgBcByYOL9lG8CXhJzu6SG8kWn5ERe3AmCIQplGkRE0i1SpgE4H/iAu//CzCq/nu4iCCgkIcpTJ6NmGoJ9m1QIKSKSclG/anYDj07xWjuHZx6kwZUzBtVnGhQ0iIikWdRPjQeA353itTOB/4mnOVIPs8o0tDSRK2h4QkQkzaIOT3wRuNLMDgH/HG5bYmZvA/4EuLgWjZPaKGcaot4aG6A900xOmQYRkVSL9Knh7lcDVwAfA7aFm38MXA18zt2/UZvmjTOzNWb2HTM7ZGb9ZvZdM1s7i/N8yMzczH5Wi3YmgTINIiIyG1EzDbj7h8zsS8DLgFXAAeDH7v5IrRpXZmadwC1ADngr4MBlwK1mdpq7D0U8zwbgr4AnatXWJCh/+FdV05BRTYOISNpFDhoA3H0n8I81ast0LgI2AE91920AZnY38BDwToIsSBRfAr4BPJUq3/tCUv7wrzbTMDRUqFWTREQkAab84Kw29e/uu+benCmdB9xRDhjC6203s9sJpoPOGDSY2RuBZwF/AHy3Vg1NgnKmodqaBmUaRETSbbpv2zsIhgGiiv61tXqnANdPsn0r8NqZDg5XtPws8Jfu3muW7hmi5YLGdtU0iIhIFaYLGt7OeNDQRlAL0A/8K7APOBp4HdADfLyGbQRYBvRNsr0XWBrh+E8DDwLXRLmYmV1MOCNk7dqqay0bnjINIiIyG1MGDe5+TflnM/sc8Bvg993dJ2y/FPh3YFPNWjhHZvZC4C3Asya2fTrhbJGrATZv3lxNtiURyh/+Ue9yCco0iIhI9MWd/gC4qvJDN/z9y8Ab425YhT4mzyhMlYGY6Crgn4A9ZrbEzJYQBEvN4e9tsbY0AcYyDRHvcglaEVJERKLPIOgGVk7x2iqgK57mTGkrQV1DpU3AvTMce3L4eNckr/UBfwZ8bi6NS5pZZRoyQabB3Ul7TYiISFpFDRpuAz5hZve5+6/KG83sucDfhK/X0g3AZ8xsQ3ldCDNbT3DXzQ/NcOyLJ9n2OYLCzfcyvlhVauTys8s0uMNosVTVVE0REVk4ogYNf0JwC+w7zGw3QSHkUQR3t9wevl5LXwmvcb2Z/RVBgebHgd0Eww8AmNk64GHgUne/FMDdb6s8mZkdBDKTvZYGuUKJ5iajpTl60FAOMHIFBQ0iImkVdRnp7cBJBCn+mwlWg7yZYGGlk919R60aGF5/CDibYAbEtQQLNG0Hznb3wQm7GkEGIfqnYQpl88WqsgwAbeFQhuoaRETSq5plpPME3/i/UrvmTHv9XcAFM+yzgwi36Xb3s+JpVTIF2YIqg4ZypiGvGRQiImmlb+QplM0XqyqChPGiyfLNrkREJH0iZRrMbDvTrw7p7n5CPE2SWptLpiGrTIOISGpFHZ74CUcGDcuB5wODBHeglIRQpkFERGYjUtDg7hdOtj1cKOmHBDMrJCFmk2loV02DiEjqzammwd0PEtzX4SOxtEbqIpsvjs2GiGps9oQyDSIiqRVHIWQWWB3DeaROZpVpaFFNg4hI2s06aDCzjJk9A/gowTLPkhCzqWkoL+ikmgYRkfSKOnuixNSzJ/qBV8bWIqm5UWUaRERkFqLOnriUI4OGLLAT+C93PxRrq6SmcoXS7DMNWhFSRCS1os6e+GiN2yF1NJtlpMcyDQVlGkRE0irSJ4eZ3WJmJ03x2lPMTOs0JMhsbjo1nmlQ0CAiklZRv26eBSya4rUe4MxYWiN1ERRCVpdpCO6KaZpyKSKSYtV8ckxVCHkCwaqQkgCFYolCyWd1e+u2TLMyDSIiKTZlTYOZvQ14W/irA1eb2UDFbh3A0whuky0JkAtrEqrNNJSPUaZBRCS9pvvkKAHF8GEVv5cfB4AvAX9U22ZKXMpBQ7WFkMExyjSIiKTZlJkGd/8a8DUAM7sV+GN3v79eDZPayIZTJqudcgnQpkyDiEiqRZ1y+eJaN0TqYyzTMIvhCWUaRETSbbqahrcA33f3A+HP03L3r8faMqmJ8jLQ7bMohGxvadIy0iIiKTZdpuEa4HSCuoVrZjiPAwoaEqC8DPTsMg1NyjSIiKTYdEHD8cDeCT/LAlBeBno2Uy7bW5rpHRqNu0kiIpIQ0xVC7pzsZ0m27BymXCrTICKSbrO+NbYk01wzDZo9ISKSXtMVQm5n6lUgK7m7nxBPk6SWlGkQEZHZmq6m4SdEDxokIZRpEBGR2ZqupuHCOrZD6iQ7h3Ua2luaxxaHEhGR9FFNQ8rMJdPQlmkiVyjhrgSUiEgaRQ4azOxEM/uamT1oZkPh8zVmtrGWDZR4ze2GVc24w2hRdQ0iImkUaRlpMzsL+AEwAnwf2AccBfwe8HozO8fdf1KjNkqMypmG1ubZFUJCEHjMJlMhIiLJFvWT4++AO4F17v4Wd/8Ld38LsB74bfh6TZnZGjP7jpkdMrN+M/uuma2NcNxmM7vazO43s2Ez22Vm3zCzVC5YFXzgN2FmVR/bFt7kSnUNIiLpFDVo2AR80t0HJ2509wHgk8ApcTdsIjPrBG4BTgLeCvwhcCJwq5l1zXD4G8L2fR54BfAh4FnAFjNbU7NGN6hsvjirO1zChEyDpl2KiKRSpOEJYA/QOsVrrcCj8TRnShcBG4Cnuvs2ADO7G3gIeCdwxTTHftLd90/cYGa3A9vD836kJi1uUOVMw2yUgw3dtEpEJJ2ifnp8EviYmR07caOZHQf8NfCJuBtW4TzgjnLAAODu24HbgfOnO7AyYAi37QT2A8fF3M6GF0emIatMg4hIKkXNNJwJLAIeMbM7GC+EPD38+aywWBKC1SHfGnM7TwGun2T7VuC11Z7MzE4GVgH3zbFdiaNMg4iIzFbUoOEMoEBw18t14QPG74L5wgn71mIS/zKgb5LtvcDSak5kZhngywSZhn+aYp+LgYsB1q6dsdYyUVTTICIisxUpaHD3hTTT4Erg+cAr3X2yQAR3vxq4GmDz5s0LaiWjODINWkpaRCSdkrIiZB+TZxSmykBMyswuJ8ggvN3db4ypbYmSzRdntYQ0KNMgIpJ2UYcngGCtBGAN0F75mrvfElejJrGVyad1bgLujXICM/sw8EHgve5+bYxtS5RsvsSyrqkmwkxPmQYRkXSLuiLkBuAbwHPLm8JnD392oJZLBN4AfMbMNrj7I2Gb1gMvIFh3YVpm9qfAZcCH3f3KGraz4WULxbFFmqqlTIOISLpFzTT8I7AWeD9wPzBaqwZN4SvAnwDXm9lfEQQpHwd2A1eVdzKzdcDDwKXufmm47Q3A54AfAreY2ekTztvv7pEyFQtFLl+ifZZLQLdrRUgRkVSLGjQ8B7jQ3f+tlo2ZirsPmdnZwGeBawmyGzcD769YpdIIMh4TB+3PCbefEz4m+glwVo2a3ZCC2RNzrGkoKNMgIpJG1awIWe/swmHcfRdwwQz77GB86KS87ULgwlq1K2m0uJOIiMxW1K+cnwA+GOE+D9LA3J1soTTrTEOmuYlMk2lxJxGRlIq6TsO1ZnYSsCNcEbJymmMtVoGUmOWLTrHks65pgKCuQZkGEZF0ijp74kLgEqBIcIfIyqGKBbUA0kJVnio52+GJ4NgmTbkUEUmpqDUNHwO+B/yRux+sXXOklsqzHmY7PAHQlmnWlEsRkZSK+umxHPiiAoZkK3/YzyXT0KZMg4hIakUNGn4GnFzLhkjtjWca5hA0KNMgIpJaUYcn3gf8q5n1ESySdMT9HtxdnyQNLhtDpqG9pUmzJ0REUipq0HBf+Pz1afap5TLSEoPxQsi51DQ0KdMgIpJSUYOGS9EMicSLY3iivaWZ3qF5XedLRETmSdR1Gj461WtmdhbwlniaI7U0Njwxh3UalGkQEUmvWeWpzWyjmV1qZtsJ7gHxunibJbUwEsOUy/aWZs2eEBFJqcifHma22MwuNrPbgQeADxMURL4bOLZG7ZMYxTN7QpkGEZG0mjZoMLMmMzvXzL4F7AW+DKwDvhDu8n53v8rd+2vcTolBLgwa2pRpEBGRWZiypsHM/g54I7AKyBKsCPk14CZgEfAn9WigxKdc09ChTIOIiMzCdIWQf0YwY+IHwIXufqD8gplpJkUCxTV7Ilso4u6Y2cwHiIjIgjFdnvqfgAHglcADZnalmT23Ps2SWsgWijQ3GS3Nc1unwT24Y6aIiKTLlJ8e7n4RcDTwJmAL8E7gF2Z2H/BBtG5D4mTzJdozsw8YYDxLoboGEZH0mfYTxN2z7v4v7n4OsJbx22N/CDDgcjN7s5m1176pMlfZfHFOQxMQZBoA1TWIiKRQ5K+d7r7X3T/l7k8Dnkswg+JEgqWl99aofRKjbL4096ChnGnIK9MgIpI2s8pVu/sWd38vwfoMFwC3xdkoqY1svjin6ZYwIdNQUKZBRCRtot57YlLunieYivm9eJojtZTNF+e0hDRMqGlQpkFEJHXm9rVTEiVbKM5pCWlQpkFEJM0UNKRIHDUN5eNzyjSIiKSOgoYUyeaLc1oNEjTlUkQkzRQ0pIimXIqIyFwoaEiRbL4059kTyjSIiKSXgoYUyRWUaRARkdlT0JAiwTLSmnIpIiKzk5igwczWmNl3zOyQmfWb2XfNbG3EY9vN7NNmttfMRszsF2b2olq3udGM5DXlUkREZi8RQYOZdQK3ACcBbwX+kGAJ61vNrCvCKf4JuAj4CPAqgmWvf2Rmz6hJgxtQvliiWPLYhieyGp4QEUmdOa0IWUcXARuAp7r7NgAzuxt4iODum1dMdaCZPR14I/B2d/9quO0nwFbgUuC82ja9MZSHE+aaacg0N5FpMnIqhBQRSZ1EZBoIPtjvKAcMAO6+HbgdOD/CsXngWxOOLQDfBF5uZm3xN7fxlDMDc800lM+hTIOISPokJdNwCnD9JNu3Aq+NcOx2dx+e5NhWYGP4c80dHB7FMBZ3ttTjcocZzzTMPWhoyzTN25TLfLHEoZE8h0byDGQL5IslSiWn6E6pRPDsDj4vzRMRqbvu9gzPWb+sLtdKStCwDOibZHsvsHQOx5ZfP4yZXQxcDLB2baRay0i+evsO/v7mh1jUnuGEVd2csXEFLzn5KJ6xZkls15hKeTghvkxDfYKG4dECP753Hzdu3ce9e/vZeWCIkgICEZExm45ZxA/e98K6XCspQUNdufvVwNUAmzdvju0j6iUnr6K7LcPuvmHuefQQX7h1G/9wyzaeu34ZH3/103jq0T1xXeoIY8MTmbmPSHW01j5oKBRL/OPPtnPlLdsYzBU4alEbz1q7lN877RhW9LSxuKOFnvYMLc1NNJnRZEZzk9HcBBb+LiKSBnO9PUA1khI09DF5RmGqLELlseumOBbGMw41d9rqJZy2esnY74dG8nzvN3v4/C3bOO/Kn/HZ1z+Dc089pibXjnN4oqPGNQ2DuQLvvHYLt287wEtPXsVFL9zAc9Yvo6lJgYCIyHxKStCwlaA2odIm4N4Ix/6+mXVW1DVsAkaBbZMfVnuLO1q48AXH86qnH8s7r/017/2XO+lsbeasp66K/VpxFkJ2tDQzMlqbTMNoocQ7vvYrfrWjj09dcBqve86amlxHRESql5TZEzcAp5vZhvIGM1sPvCB8bTr/AbQwoWDSzDLA64Eb3T0Xe2urtKK7ja+//bmcdHQP7/3nO9l7aCT2a4zENOUSoL21eex8cbvixw9yxyO9ChhERBpQUoKGrwA7gOvN7HwzO49gNsVu4KryTma2zswKZvaR8jZ3v5NguuXnzOwdZvYSgumWxwN/Xcf3MK2utgxffNOzyJdK/J9/j38yR7zDE001qWnY+tghrvrpw7zhOWu44NmrYz+/iIjMTSKCBncfAs4GHgSuBb4BbAfOdvfBCbsa0MyR7+ttwFeBy4DvA2uAc9z9NzVuelXWLe/ifS95Cjfdt48tO+IttRgLGuZ47wkIhydqEDT87Q/uZ0lHC5ece3Ls5xYRkblLRNAA4O673P0Cd1/k7j3u/mp331Gxzw53N3f/aMX2EXf/c3c/2t3b3f157n5bHZsf2Vufv44V3a187qaHYj1vtlCuaYhn9kTcNQ137urjZ9ue5D0v3sjijvqvYyEiIjNLTNCQFp2tGd5+xvH8bNuTbHtiILbz5sqZhtZ41mmIO9Nwzc930NOW4Q3PjW9dDBERiZeChgb02mevIdNkfPOXu2M7Z9zDE3HWNPQOjfL9u/fyms2r6W5LyoQeEZH0UdDQgFb2tPGyTUfx3TsfpVCMZz2EbL5Ek0FL89zXOmhvaSZfdPIxte2/7tlLoeS8RsWPIiINTUFDg/q9px9L79AoW3bOtHZVNNl8kfaWZiyGlRLLq4/FlW34/t172bCii03HLIrlfCIiUhsKGhrUmU9ZSWumiRu37ovlfNlCMZbpljBeFxFHXUPv0Ch3PHKAV552TCwBjYiI1I6ChgbV1ZbhhRtX8OP7Ho/lfNl8KZb7TsCETMPo3Icnbt/2JCWHs0+KfxVMERGJl4KGBvaip6xkd+8Iu3sr7+pdvZF8fJmGctAQR6bhpw/uZ1F75rB7coiISGNS0NDAnn/CcgB+/vCTcz5XLl+kLa6goTX4ZzPXoMHd+e+HnuSME1fQrJtRiYg0PAUNDWzjqm5WdLfx84cPzPlc2XwploWdYHwp6rku8LTzwDCP92d5/gkr4miWiIjUmIKGBmZm/M4Jy7njkbkHDSP5Ymz3XI9r9sSvw5khz1m/bIY9RUSkEShoaHDPXruEff25Od/5cni0SGcMq0FCsIw0zH144te7+uhpy3Diqu44miUiIjWmoKHBPX3NEgDu2n1wTucZGS3Q0RrPaosdMQ1P/GZnH89ct5Qm1TOIiCSCgoYGd/Ixi2hpNu6cY9AwPFqks4FmTwzmCjywb4BnrV0SS5tERKT2FDQ0uPaWZjYdsyiGTENxbFhhzm1qnXtNw317+3GHU49bHEubRESk9hQ0JMCpqxdzz6P9uPusjnd3hvMx1jTEUAi59dFDAJxyrIIGEZGkUNCQACcfs4jBXIE9fbMrhhwtliiWPLagoaW5iUyTzWl44t69/SzvauWoRW2xtElERGpPQUMCnBzeyOn+xwdmdXy5YDGuQkgIsg0jc1hGeutj/Ww6dpHuNyEikiAKGhLgqUf1AEEdwGyUMwJxZRogqGsYyRdmdWy+WOKhfYNsOlZ3tRQRSRIFDQnQ1ZZh3fJO7n98dkHD8Gj8QUNna/PYeau188AQo8USJx3dE1t7RESk9hQ0JMRJR/dw/945Dk/ENOUSoLM1w1BudkHDQ/sGAThxlYIGEZEkUdCQECeu6mFn7zCjherrCMYzDfHVNHS1NjM8OrvhiYeeGMQMTliplSBFRJJEQUNCbFzVTbHk7DwwVPWx5Q/3uNZpKJ9rtsMTDz0xyOqlHbG2R0REak9BQ0KUv5U/vH+w6mNHalDT0NWamX2mYd+AhiZERBJIQUNCbFjZBcDD+2eTaahBIWRb86xqGool55Enh9iom1SJiCSOgoaE6GrLcMzidh5+ovpMw3C+vE5DvLMnZrO402MHRxgtlNiwoiu2toiISH0oaEiQE1Z2s21WwxPBMEK8hZAZhnLVD0/sPDAMwNrlnbG1RURE6kNBQ4Icv6KLHU/Ofngi7imXuUKwPHU1dvYG7V+3XJkGEZGkUdCQIOuWd9KfLXBoOF/VcSOjRdoyTTQ3xbdkc7k+otpiyF0HhmltbuLoRe2xtUVEROojEUGDmTWZ2SVmtsPMsmZ2l5ldEOG4RWb2ETP7uZkdMLOD4c+vrkOzY7dmWZDS39U7XNVxw6Px3eGyrLOtHDRUV9ew88Awq5d1xBrAiIhIfSQiaAA+DnwUuBJ4BXAH8G0zO3eG49YC7wZ+ArwZeD3wIPA9M3tPzVpbI2vnFDTEV88A45mGausadvYOs15DEyIiiRTvJ0kNmNkq4APA5e7+mXDzrWa2Ebgc+ME0h28HNrj7xE/ZH5nZGuCDwBdq0eZaKWcaynUBUY3kC7EvpFQOQqrJNLg7uw4M8bzjl8XaFhERqY8kZBpeDrQC11Vsvw441cyOn+pAdx+qCBjKtgDHxtfE+uhuy7C8q5XdDTA80TWLoOHJwVGGRous08wJEZFESkLQcAqQA7ZVbN8aPm+axTlfBNw/l0bNlzXLOqsfnsgVY505AeNrPgxVUQi5a2zmhIIGEZEkSkLQsAw46O6Vc/t6J7wemZldDJwO/O10+5jZFjPbsn///qoaW2trZxE0DOYKdLfFOxLVFRZCjlSRaRhbo2GZahpERJKo7kGDmb3UzDzC47YaXPss4PPA1939G1Pt5+5Xu/tmd9+8cuXKuJsxJ2uXdfLYwSz5YvS7XQ6NFuhujzloCIcnqimE3HlgGDNYs6wj1raIiEh9zEch5M+BkyPsV/463QcsMTOryDaUMwy9RGBmzwFuAG4B3hGxrQ1n7bJOiiVn78Fs5FUVB7MFumLONHS0Vj/lclfvMMcu7qAto7tbiogkUd2DhrAwsZp6gq1AG3ACh9c1lGsZ7p3pBGZ2KvAj4LfABe5e3epIDaQcKOzqHY4eNNRieGIWhZA7DwyNTRsVEZHkSUJNww+BPPCmiu1vBu5x9+3THWxmJwI/Bh4BXuXuIzVpZZ1Uu1ZDvlgiVyjFHjS0tzRhVt2KkDsPDKsIUkQkwRp+nQZ3f8LMrgAuMbMB4DcEizSdDZw3cV8zuxlY5+4bw99XEQQMrcBfA5vMDluJ8E53z9X+XcTnqEXttDY3RQ4ayjUHcQ9PmBmdLc2RMw3DowUODI2OrTUhIiLJ0/BBQ+jDwCDwPuBo4AHgde7+nxX7NXP4e9oErAt/rtwX4HhgR6wtrbHmJmP10o7IazUMhkFDd1v8dQRdbRkGs9EyDY/2BQme1UtVBCkiklSJCBrcvQhcFj6m2++sit9vAxbcTQ7WLOuMvCrkeNDQEns7utszY+efyZ6DQdBw3BIFDSIiSZWEmgapsGZZB3v6opVmjA9PxJ9p6GlvYSBi0FDONBynTIOISGIpaEig1Us7OTicZyA78ySQwVxQc9AT8zoNAD1tmUhtAHj04AiZJmNVj26JLSKSVAoaEqhcF/DowZmzDeWag7gLISEIRKqpaThmSbtuiS0ikmAKGhKoXBewp3fmoGFseCLmW2NDcAOtgahBw8ER1TOIiCScgoYEWr00mLYYKdMQBg01GZ5ob4lcCPlo3wjHLdF0SxGRJFPQkEArultpyzSxp2/maZeDNVqnAcZnTxRLlfcSO9xoocS+gayKIEVEEk5BQwKZBWs1RJlBMZQr0JppoqU5/v/Ui8LsxUy3x97Xn8UdVmt4QkQk0RQ0JNTqpZ2RgobBXIGeGmQZgLGlqWcqhtyj6ZYiIguCgoaECjIN0YYnajE0AUFNAzBjMeSjWthJRGRBUNCQUMct7aBvOD82O2IqQzUMGrrD4YnB3PRrNZQXdjp6sdZoEBFJMgUNCRV1BsVAtnbDE+UZGf0zZhqGWdnTRntL/KtSiohI/ShoSKjyAk8zDVEcGsmzuDP++04AY8HITDUNWqNBRGRhUNCQUONBw/SZhkMjeRZ31ChoiFrT0DeiIkgRkQVAQUNCrexuC9dqmL+gIUpNQ6nkPHYoq+mWIiILgIKGhDIzjpthBsVoocTwaJElNQoaulqbMZs+0/DkYI7RQmksMyIiIsmloCHBjlsy/QJPh0aCDECtahrMbMb7T+wO21cu3BQRkeRS0JBgMy3wNBY01CjTALCks4WDw6NTvl7OhKxZpkyDiEjSKWhIsNVLO+gdGmV4imWcD40EH+a1DBqWdbbSNzx1TcPYapC6WZWISOIpaEiwcp3Ao1NkG8qZhiWdrTVrw5LOVvpmyDSs6G6lo1VrNIiIJJ2ChgQr1wlMNURRj+GJpZ0tMwQNIxynegYRkQVBQUOCrZlhgaeD4bBBrWZPACztauXg0PTDE5o5ISKyMChoSLAV3W20TrNWQznTsKimmYZWBnIFRgulI14rlZxHFTSIiCwYChoSrKnJpp12eXA4T097huYmq1kblnYF9RIHR44cotg/mGO0WNJ0SxGRBUJBQ8KtXtrBniluWtVfw9Ugy5aGa0AcnGQGRXnYRJkGEZGFQUFDwq1e2sGjU9Q01HIJ6bKl4cyM3qEjMw3lDMgaZRpERBYEBQ0Jt3ppJ08OjjIyWjzitScHcyzvbqvp9ctBw2QLPO0ZWw1SmQYRkYVAQUPCja3VcPDIbMMTAzlW9dQ4aOgKMhmTLfC0u3eYFd1ttLdojQYRkYVAQUPClYOG3RXFkKWSs78eQcMMwxPKMoiILByJCBrMrMnMLjGzHWaWNbO7zOyCWZxng5kNm5mb2cZatLXeysszV86g6BsepVDymgcN7S3NdLQ00zdp0DCsoEFEZAFJRNAAfBz4KHAl8ArgDuDbZnZulef5InAo3qbNr1U9bbQ02xFLST8xkAteX9Re+zYsamNfeL2yYsl59OCIpluKiCwgDR80mNkq4APA5e7+GXe/1d3fCdwKXF7Fed4IPBP4ZG1aOj/G12o4vKZhLGiocaYB4JjF7eytmPa5p2+YfNHZsKKr5tcXEZH6aPigAXg50ApcV7H9OuBUMzt+phOY2VLgCoLg42DcDZxvq5d2HlHTsK8/C8CqntpnGo5d0sFjFUHDI/uHADhhlYIGEZGFIglBwylADthWsX1r+Lwpwjk+Bdzv7tfG2bBGcfyKLh55YhB3H9u2f2x4ovaZhmMXd7BvIEexNH79h/cPArBhRXfNry8iIvWRhKBhGXDQJ34iBnonvD4lM3sh8Bbg3VEvaGYXm9kWM9uyf//+qho7H048qpuBXIF9/eN1BU/0Z+lpz9RluuMxS9oplpwnBrJj2x7eP8TSzpaxZaZFRCT56h40mNlLw9kLMz1ui+FarcBVwGfd/d6ox7n71e6+2d03r1y5cq7NqLkTV/UA8OC+gbFt9VijoezYxcEMiYlDFA/vH2TDSmUZREQWksw8XPPnwMkR9itX9vUBS8zMKrIN5QxDL1N7P7AU+LyZLQm3lcv5e8ysx90HJjswSZ5yVPDh/OC+AV70lCDI2XlgmOPqNHPhmCVB3cRjB7M8e12w7ZH9Q5x9UuMHXCIiEl3dgwZ3Hwbur+KQrUAbcAKH1zWUaxmmyyBsAo4GHp3ktd8AdwHPqKItDWl5dxvLu1p5aF9QR1AsOQ/vH+QFG5fX5frHLgkyDXsPBZmGff1ZnhzM8ZSjeupyfRERqY/5yDRU64dAHngT8LEJ298M3OPu26c59nLgmopt5wAfDI9/IL5mzq8Tj+rm/nB4YlfvMLlCiRPr9KG9qL2FRe0ZdhwIkkN37T4IwDPXLqnL9UVEpD4aPmhw9yfM7ArgEjMbIMgQvB44Gzhv4r5mdjOwzt03hsfeT0VWw8zWhz/+P3evnJGRWM9Ys5R/+tkjDI8Wxmob6vlN/9TVi7l7z0EA7tpzkEyTccqxi+t2fRERqb0kzJ4A+DBwGfA+4EfAC4DXuft/VuzXTAICoVr4nROWky86v97Zx0Nh0HDiqvoVIj5zzVLu2zvAyGiRu3Yf4qRjenSjKhGRBSYRQYO7F939Mndf5+5t7n6au39nkv3Ocvf1M5zrGne3hZRlANi8bimZJuMXDx/gvr0DHLekg662+sVPz1y7hGLJ+eWOXn67+yBPX72kbtcWEZH6SOW38oWoqy3DaasX85937+Xx/iyvefbqul7/GWuWAHDpf2xlMFfg1c88rq7XFxGR2ktEpkGi+dOXnMjuvmGKJeddLzqhrtde3t3GS08+iof3D/H01YvZvG5pXa8vIiK1p0zDAnLWU1fxqQtOI5svsnZ5/e8u+cU3PYuv/PcjvOjElZhZ3a8vIiK1paBhgXnt5jXzdu3WTBPvefHGebu+iIjUloYnREREJBIFDSIiIhKJggYRERGJREGDiIiIRKKgQURERCJR0CAiIiKRKGgQERGRSBQ0iIiISCQKGkRERCQSBQ0iIiISiYIGERERiURBg4iIiESioEFEREQiMXef7zY0NDPbD+yM8ZQrgCdjPF8aqQ/nTn04d+rDuVMfzl0t+nCdu6+c7AUFDXVmZlvcffN8tyPJ1Idzpz6cO/Xh3KkP567efajhCREREYlEQYOIiIhEoqCh/q6e7wYsAOrDuVMfzp36cO7Uh3NX1z5UTYOIiIhEokyDiIiIRKKgoQ7MbI2ZfcfMDplZv5l918zWzne75puZvcbM/s3MdprZiJk9YGZ/a2Y9FfstNbN/NLMnzWzIzG4ys1MnOV+7mX3azPaG5/uFmb2ofu+oMZjZD83Mzeyyiu3qx2mY2blm9lMzGwz/P91iZmdPeF39Nw0ze4GZ3WhmT5jZgJn9xszeXrFPpL4xsyYzu8TMdphZ1szuMrML6vduas/MVpvZP4R9MBz+P7t+kv1i7zMzu8jM7jezXPh3912RG+7uetTwAXQCDwH3AK8Gzgf+B3gY6Jrv9s1z39wB/CvwJuBM4P3AwXB7U7iPAT8D9gB/AJwD/IRgXvLqivN9Izz+IuAlwHeBEeAZ8/1e69infwDsBRy4bMJ29eP0/fZOIA98FngZ8HLgg8Cr1H+R+u+08D3eGv6NexlwVfjv8I+r7Rvgb4Ac8AHgxeG5SsC58/1eY+yzs4B9wA+AH4V9tX6S/WLts/A8pXD/FwOXhb//caR2z3fHLfQH8D6gCGycsO14oAD8+Xy3b577ZuUk294S/s9zdvj7+eHvL56wz2KgF/j8hG1PD/d724RtGeAB4Ib5fq916s+lwOPhh1pl0KB+nLrf1od/hN8/zT7qv+n78BPAKNBdsf0XwC+q6RtgVfjh97GKc90M3D3f7zXGPmua8PM7Jgsa4u6z8NgngK9V7Pd/CQLglpnareGJ2jsPuMPdt5U3uPt24HaCP0Sp5e77J9n8q/D5uPD5POAxd791wnGHgP/g8P47j+Cb4rcm7FcAvgm83MzaYmx6o/okcI+7/8skr6kfp/Z2gm9aX55mH/Xf9FoJ3vdIxfZDjA+DR+2bl4fnu67iXNcBp5rZ8fE2fX64eynCbnH32e8AKyfZ71pgOXDGTA1S0FB7pxAMTVTaCmyqc1uS4Mzw+b7webr+W2tm3RP22+7uw5Ps1wpsjLuhjcTMziDI0rxnil3Uj1M7A7gfeIOZPWxmBTPbZmYT+1L9N71rwufPm9mxZrbEzMrp9M+Gr0Xtm1MIvjVvm2Q/SNffzbj77JTwufLfcuS+VdBQe8uAvkm29xKkkyVkZscBlwI3ufuWcPN0/QfjfTjTfsviamejMbNWgvHLz7j7A1Pspn6c2rHAicCngcuB3wV+DFxpZu8L91H/TcPd7yEYoz8feJSgD74AvMvdvxnuFrVvlgEHPcybT7NfGsTdZ+XnynNG7tvMTDuI1EP4Te16glqPt81zc5LmL4EOgsImqV4T0ANc6O7fDbfdElayX2Jmn5+3liWEmZ0I/BvBN9Z3EQxTnA982cyy7v6N+WyfxEdBQ+31MXlGYaoIMnXMrINgbHgDcKa775nw8nT9V369/Lxumv16J3kt8cKpux8mKKRqqxgzbzOzJcAA6sfpHCDINPy4YvuNBLMkjkH9N5NPEIy9v8rd8+G2m81sOfD3ZvYvRO+bPmCJmVnFN+eF3oeTibvPyv9OlxLMsppqvylpeKL2tjI+jjTRJuDeOrel4ZhZC/AdYDPB1KD/qdhluv7b5e6DE/Y73sw6J9lvlCPH+haKDUA7QWFT34QHBFOv+oBTUT9OZ+sMr5dQ/83kVOCuCQFD2S8JCuxWEb1vtgJtwAmT7Afp+rsZd5+V/61X/luO3LcKGmrvBuB0M9tQ3hCmPV8QvpZaZtZEMAf5bODV7n7HJLvdABxnZmdOOG4R8Hsc3n//AbQAr52wXwZ4PXCju+fifwcN4bcEc60rHxAEEi8m+MOifpza98Lnl1dsPwfY4+6Po/6byePAM8L6momeB2QJvsFG7ZsfEmQt3lRxrjcTzA7aHn/zG1bcffYLgqmVk+3XSzCrb3rzPVd1oT+ALoI/2v9DMMZ3HnAX8AgVc5rT9gC+RLieAHB6xWN1uE8T8HNgN/AGgj/st4X/wNdUnO+bBN+s30FQtf0dgj9Yz5rv9zoPfVu5ToP6ceq+MuAWgmGKdxEUQn4l7MML1X+R+vA1YX/9KPw797vAleG2K6rtG4KC1Czw5wQFll8iyPi8ar7faw367TUT/hb+cfj7mbXqs/DfeCn8u3sWQfF5CXhPpDbPd6el4QGsJSgS6icYX/53Jln5K20PYEf4P8pkj49O2G8ZweIjvcAwwYIlT5/kfB3AFQTferLA/wPOmu/3OU99e1jQoH6csb8WEVT77yNI+94NvFH9V1UfvoIgkNof/p37LfBuoLnavgGagb8CdhJMJbwbeM18v8ca9NlUf/9uq2WfEayA+mC430PAu6O2WXe5FBERkUhU0yAiIiKRKGgQERGRSBQ0iIiISCQKGkRERCQSBQ0iIiISiYIGERERiURBg4iIiESioEEkpczMIzx2mNn68OcL57vNZWZ2nJkNmdnmOl3PzOxOM/vLelxPpFFpcSeRlDKz0ys2fY9gifOPTtiWI7iJzTOBh919f31aNz0z+7/AKnd/VR2v+fsEK0Ke4O5putOiyBgFDSICgJntAH7m7m+e77ZMx8yOIrgHxO+7+/freN1mYA/wWXf/VL2uK9JINDwhItOabHjCzK4xsz1mttnMfm5mI2b2gJm9Mnz9z8OhjX4zu97MVlacM2Nml5jZ/WaWM7PHzOzvzKw9QpMuJLi3wY8qznmbmf3MzM4xs9+GbbrTzJ4XXu8TZrbXzHrD9ndVtOfjZvawmWXN7MnwXGeU93H3IvBtghsHiaRSZr4bICKJtQj4OvAZ4DHgw8C/mdkXgKcA7wGOAj5HcDOo10049jqC20p/kuDukScDHwfWAxfMcN1zgF+4e2GS1zYCnwb+BhgEPkVw6+obCP7eXRhe69PAE0C5RuGDwJ+F7+G34XvbTHCTqol+CrzXzDa4+yMztFNkwVHQICKz1QO8y91/CmBmjxHURLwK2BR+M8fMnkbwQdvs7kUzeyHweuCt7v718Fw3mVkvcJ2ZPcPdfzvZBc3MgOcBn52iTcuB55c/0M2sCbgeON7dXxru8yMzexHwWsaDht8BbnT3v59wrv+Y5Px3hs+nE9zeXiRVNDwhIrM1VA4YQveHzzeVA4YJ2zPAMeHv5xDcfvo74bBAxswywI3h6y+a5ppLCG4VPFVB5oMVGYBym35Usd/9wOowCAH4FXCumf2NmZ1hZq1TnL983WOnaaPIgqWgQURm6+DEX9x9NPyxr2K/8vZyvcIqoBUYAvITHk+Ery+f5prlc+SmeH2qa0+2PQM0h79/Avhr4Dzgv4EDZvZVM1tRcdxI+NwxTRtFFiwNT4hIvR0AssALp3j9sRmOBVgaZ4PcPU9QX/FJMzuaYIjlCqCTYCilrFzj8GSc1xdJCgUNIlJvPyQoPFzs7jdXc6C7j5rZdmBDTVoWXONx4B/N7FzgaRUvHx8+P1Cr64s0MgUNIlJX7n6bmf0LQU3DFcAvgRLBzIlzgQ+6+4PTnOKnwHPjbJOZXU9QxPkbgqGMZxLUXlxVsevzCIZS7ojz+iJJoaBBRObDm4H3Am8nmOaYA3YQFCzum+HYbwFvMbP17r4jpvb8lGA2xXsIhiR2EUzX/JuK/V4F3ODuwzFdVyRRtCKkiCRKOI3yIeCr7n5ZHa97LMFKlL9b7bCKyEKhoEFEEsfM3kRQqHh8vb71m9lngae7+9n1uJ5II9LwhIgk0T8DxxHUQdxb64uF6zk8Dlxd62uJNDJlGkRERCQSLe4kIiIikShoEBERkUgUNIiIiEgkChpEREQkEgUNIiIiEsn/B7DXXIngA6krAAAAAElFTkSuQmCC\n" }, "metadata": { "needs_background": "light" } } ], "source": [ "#NBVAL_IGNORE_OUTPUT\n", "from examples.seismic import RickerSource\n", "\n", "f0 = 0.010 # Source peak frequency is 10Hz (0.010 kHz)\n", "src = RickerSource(name='src', grid=model.grid, f0=f0,\n", " npoint=1, time_range=time_range)\n", "\n", "# First, position source centrally in all dimensions, then set depth\n", "src.coordinates.data[0, :] = np.array(model.domain_size) * .5\n", "src.coordinates.data[0, -1] = 20. # Depth is 20m\n", "\n", "# We can plot the time signature to see the wavelet\n", "src.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Similarly to our source object, we can now define our receiver geometry as a symbol of type `Receiver`. It is worth noting here that both utility classes, `RickerSource` and `Receiver` are thin wrappers around the Devito's `SparseTimeFunction` type, which encapsulates sparse point data and allows us to inject and interpolate values into and out of the computational grid. As we have already seen, both types provide a `.coordinates` property to define the position within the domain of all points encapsulated by that symbol. \n", "\n", "In this example we will position receivers at the same depth as the source, every $10m$ along the x axis. The `rec.data` property will be initialized, but left empty, as we will compute the receiver readings during the simulation." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "output_type": "display_data", "data": { "text/plain": "
", "image/svg+xml": "\n\n\n \n \n \n \n 2021-04-15T19:25:24.928129\n image/svg+xml\n \n \n Matplotlib v3.4.1, https://matplotlib.org/\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAc8AAAGDCAYAAABN4ps8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAA2IElEQVR4nO3deZhcRb3/8fcnC7L+EETAsIWIYoIbGAXBC0kUARFQQUTZZZOriAsiAgYIeEFRQC5XIIACgXtBNlkEQYQERdncgABhSwjBIAkJeyBM8v39UadJp9Mzfc70Mj09n9fz9HOm65yqU30ymW9XnTpVigjMzMwsv0F9XQEzM7P+xsHTzMysIAdPMzOzghw8zczMCnLwNDMzK8jB08zMrKCWB09J60m6UtKLkl6SdLWk9XPmXV7SqZJmS1og6S+Stm52nc3MzMq1NHhKWhG4DXgfsC+wN/Ae4HZJK+Uo4gLgIGA88FlgNnCzpA83pcJmZmZVqJWTJEg6HDgN2DgiHs/SNgQeA46MiNN6yPsh4B/AVyPiV1naEGAqMC0idm5y9c3MzIDWd9vuDNxVCpwAETEduBPYJUfeN4HLy/J2AZcB20l6W+Ora2ZmtqxWB89NgAerpE8FRuXIOz0iXquSdzlgo/qrZ2ZmVlurg+fqwPwq6fOA1erIW9pvZmbWdEP6ugLNJulg4OD0buhHYI0+rY+ZWXt5gYjXBLCRtEzXXlGz4eaI2L4BFWtrrQ6e86newuyuVVmZd4Nu8sKSFuhSImIiMBFAGhZvxVEzMyP78wjAAuDrdZZ27ABpobS623Yq6d5lpVHAQznybpg97lKZdyHw+LJZKqz4KgzpylHNzJAu+Ng9zuM8ztOqPO1ev07Ms+Krb70VMLTO10DR6uB5HbCFpBGlBEnDga2yfT25nvRv88WyvEOALwG3RMQbNc++6kuw96R8v1xDutKx293sPM7TsDyDWMyOPMqxTGFHHmUQi9umbn2ep93r16l5Vn3prSSRuiPreQ0UrQ6e5wEzgGsl7SJpZ+Ba4Gng3NJBkjaQ1CVpfCktIv5OekzlDEkHSvok6TGVDYHjctdg2L9q/3KVfqmG/QsGL3Ye52lInkEs5mYu4f+4iuOZzP9xFTdzCYMGL+zzuvV5nnavXyfnsV5pafCMiFeBccCjwCTgUmA6MC4iXik7VMDgKvXbH/gVcBLwW2A9YPuI+FvuSgzt6vmXq/yXamiX8zhPw/LsMGgam/MMq7CQwcAqLGRzZrHD2PP7vG59mqfd6zcQ8mTcbZtfy+e2jYiZEbFrRPy/iFglIj4XETMqjpkREYqI4yvSF0TEdyJi7YhYPiI2j4jJhSsxtAvWnQWbVYm5m/0t7av4pXIe56k3z6bv+jsrsnCpXSvyJh/umtvndevTPO1ev4GSB3fbFjEwV1V5cwjMWhf+ttmy+/62Wdr35hDncZ6G5vn77E15jeWW2vUaQ/nHkDX6vG59mqfd6zdQ8uCWZxEDL3i+OQT+NQwm7Q1dVb4ndQ1J+/41bMkvl/M4TwPy3LR4Y+5mHV5mORYBL7Mcd7MuN91+YJ/XrU/ztHv9BkIeK2zgBc9a/4lh6V+uRYOcx3kakmcxg9iOvfgyu3IcY/kyu7Ide7F40XJ9Xrc+z9Pu9evkPGXcbZtfS1dV6WtaadVg4WG1/xOXDOlK9wf+tpnzOI/ztCJPu9evE/Ms99/Eqy8KYIQUJ+bL2a294K8RMbrOYtrewAqenmHIzKzCRCL+JYB3S/FfdZa2xwAJngOplW1mZj0oDRiy2gbePU8zM7M6ueVpZmaAW55FOHiamdlbHBTy8XUyMzPALc8iHDzNzAxY8pyn1eYBQ2ZmZgX5S4aZmQHuti3CwdPMzAB32xbh62RmZoBbnkX4nqeZmVlBbnmamRngbtsifJ3MzAxwt20RDp5mZga45VmEr5OZmQFueRbhAUNmZtYyknaTdJWkpyQtkDRN0smSVilYzlGSQtKfmlXXnrjlaWZmQMtankcAM4GjgVnApsDxwFhJW0bE4loFSBoBHAs818R69sjB08zM3tKCoLBTRMwpez9F0jzgImAMcFuOMs4GLgU2po/imIOnmZkBWcuz3qjQ1fPuisBZcm+2XadW8ZK+AmwGfBm4umDtGsbB08zMAJBgSJODZze2ybYP93SQpNWA04EjI2KepF6drBEcPM3MrJHWkHRf2fuJETGxu4MlrQNMAG6NiPu6Oy5zKvAocGHdtayTg6eZmQGp5Tl0cN3FzI2I0fnOp5WBa0nt1f1rHPsfwD7AZhERddeyTg6eZmYGNKjbNve5tAJwPTAC2CYiZtXIci5wATBL0tuztCHA4Oz9goh4o0nVXYaDp5mZAQ0aMJTnPNJQ4EpgNLBtRDyQI9vI7PW1KvvmA98GzmhUHWtx8DQzs5aRNIj0mMk44LMRcVfOrGOrpJ0BDAYOAx5vSAVzcvA0M7NEpFDUXP8DfBH4EfCqpC3K9s2KiFmSNgCeACZExASAiJi8THWlF4Ah1fY1m4OnmZklrZkZfodse0z2KncCabahUhhv2ylkHTzNzCxpQfCMiOE5jpmR1abWcWPqr1HvOHiamdkSjgq5tG2T2MzMrF35O4aZmSWtGTDUERw8zcwsac2AoY7gy2RmZomDZ26+52lmZlaQv2OYmdkSvueZi4OnmZkl7rbNzZfJzMwSB8/cfJnMzGwJd9vm4gFDZmZmBbnlaWZmibttc/NlMjOzxMEzN18mMzNLHDxz82UyM7MlPGAoFw8YMjMzK8gtTzMzS9xtm5svk5mZJQ6eufkymZlZ4vU8c2v5PU9J60m6UtKLkl6SdLWk9XPkGy1poqRHJL0maaakSyVt2Ip6m5mZlbS05SlpReA24A1gXyCAk4DbJX0wIl7tIfsewCbAmcBUYB3gh8B9kj4cEU83tfJmZp3O3ba5tfoyHQSMADaOiMcBJN0PPAYcApzWQ94fR8Sc8gRJdwLTs3LHN6XGZmYDiYNnLq3utt0ZuKsUOAEiYjpwJ7BLTxkrA2eW9hQwh9QKNTOzepTuedbzGiBaHTw3AR6skj4VGFW0MEkjgTWBh+usl5mZlbpt63kNEK0OnqsD86ukzwNWK1KQpCHAOaSW5wX1V83MzCyf/vw94SxgS2DHiKgWkAGQdDBwcHq3aksqZmbWL3nAUG6tvkzzqd7C7K5FWpWkU0gBcd+IuKWnYyNiIjAx5RsW+atqZjYAOXjm0urLNJV037PSKOChPAVIOgb4PnBYRExqYN3MzAY2T5KQW6vveV4HbCFpRClB0nBgq2xfjyR9k/Rc6DERcVazKmlmNiB5wFBurQ6e5wEzgGsl7SJpZ+Ba4Gng3NJBkjaQ1CVpfFnaHsAZwO+A2yRtUfYqPFLXzMyst1r6PSEiXpU0DjgdmET6nvMH4FsR8UrZoaXOg/Lgvn2Wvn32KjcFGNOkapuZDQweMJRbyy9TRMwEdq1xzAzSP2N52n7Afs2ql5mZ4XueOfk7hpmZJW555tbyVVXMzMz6O3/HMDOzxC3P3HyZzMwscfDMzZfJzMyW8IChXBw8zcwsccszNw8YMjMzK8jfMczMLHHLMzdfJjMzW8L3PHNx8DQzs8Qtz9x8mczMLHHwzM0DhszMzArydwwzM0u8GHZuDp5mZpa42zY3XyYzM1vCUSEX3/M0MzMryN8xzMwscbdtbr5MZmaWeMBQbg6eZmaWuOWZmy+TmZkt4aiQiwcMmZlZy0jaTdJVkp6StEDSNEknS1qlRr7RkiZKekTSa5JmSrpU0oatqns5f8cwM7OkNfc8jwBmAkcDs4BNgeOBsZK2jIjF3eTbA9gEOBOYCqwD/BC4T9KHI+LpZle8nIOnmZklrbnnuVNEzCl7P0XSPOAiYAxwWzf5flyRD0l3AtOBg4DxTahrtxw8zcwsaUHwrAyAmXuz7TpF8kXEU5Lm9JSvWRw8zcxsib55VGWbbPtwkUySRgJrFs3XCA6eZmbWSGtIuq/s/cSImNjdwZLWASYAt0bEfd0dVyXfEOAcYA5wQW8r21sOnmZmljSm23ZuRIzOdTppZeBaoAvYv+B5zgK2BHaMiPkF89bNwdPMzJIWTpIgaQXgemAEsE1EzCqQ9xTgYGDfiLilSVXskYOnmZklLQqekoYCVwKjgW0j4oECeY8Bvg8cFhGTmlTFmhw8zcysZSQNAi4FxgGfjYi7CuT9JnAScExEnNWkKubi4GlmZks0f7Tt/wBfBH4EvCppi7J9syJilqQNgCeACRExAUDSHsAZwO+A2yryvRQRDzW95mUcPM3MLGlNt+0O2faY7FXuBNJsQ6W5jsqnkN0+S98+e5WbQppgoWUcPM3MLGnNJAnDcxwzI6tNedp+wH7NqFNvOHiamdkSXs8zFwdPMzPreNlkDNsCWwDDgBWAucA0UrfvlB4mpV+Gg6eZmSUduBi2pG2A7wHbkdrVs0izEi0APgTsRJpUfrak84DTIuKlWuV22GUyM7Ne67DgKem3wFjSZAy7A3+MiLkVxwwC3k8Kol8Bvi5pr4i4uaeyO+gymZlZXToseAKPAgdExLPdHZB11d6fvX4kaWdg1VoFd9ZlMjOzukQHDRiKiG/3Is91eY4bVPsQMzOzziTpHb3J5+BpZmYAhGDRkPpe7UrSQZK+V/b+A5JmAc9Juk/S2kXKc/A0M7Okg4MncBhphG3JacALwLdI9zgnFCmsvT+qmZm1TAi6Btfbpsr9qGSrbQA8AiBpVWAb4HMRcaOk54GTixTmlqeZmQ0Eg1gS2T8BBDA5e/80sGaRwtzyNDMzAEJi0ZB6w8LChtSlCR4DdgRuA/YA/hwRr2X7hgHzihTm4GlmZm9ZNLiDnlVZ2k+BSZL2BVYjLYtWMpb0nGduDp5mZgZAIBZ16MzwEfG/kmYCmwP3RsQdZbv/DeR6vrOkUPDMhvKWT6g7PSLato1uZmb5BaKrg4KnpKOBayLiYYCI+BPwp8rjIuK4omXXHDAkabSkcyRNB54B7gXuAB4CXpR0h6RDJa1S9ORmZmZNtDfwoKRpkn4i6eONKrjblqek0aQ+4q2BB4AbgL+zZDb61YENSU3gU4BTJP0E+FlEvN6oCpqZWess6qC7eRExUtL7gM8BuwDflfQcaaL4a4A/9Lb3tKerNAU4Dzi01OTtjqTls4odSWrNntibypiZWd/pxHueEfEISxp4a7MkkF4DvCHp5uzn3+ZZiqykp+D57p5moq+o3OvA5cDlktbKe3IzM2sfnRg8y2Ux7RzgnOxW446kQHo2sIKkyRGxXZ6yug2eeQNnlXz/7k0+MzPre50cPMtFxMvAZcBlkoYCnyQF0lx6NcOQpEGVrwJ515N0paQXJb0k6WpJ6/eiDkdJCknLjJwyMzPLKyLejIjfRcShefPkujMsaQXgONJDpetWyRd5ypK0Iml2hzeAfbN8JwG3S/pgRLyasz4jgGOB5/Icb2ZmtXXaoyrlskbewaQ4th6wfMUhEREb5C0v77CqXwB7kkYoXUbv5186CBgBbBwRjwNIup80bdIhpFnu8zgbuBTYGE/0YGbWEOmeZ8f+Sf0J8B3SUyP3Uuc8gnmv0s7AERFxZj0ny8q5qxQ4ASJiuqQ7SX3NNYOnpK8AmwFfBq6usz5mZlamg+957gWc2JsJEarJe6/yDaDHx1Vy2gR4sEr6VGBUrcySVgNOB46MiEKT+JqZ2YA2hDTBT0PkDZ4Xkmahr9fqwPwq6fNIE/XWcirwaFafXCQdnK0Sfh+8VjuDmdkAVXpUpZ5XG7sSyPUYSh55u21/CJwt6RbgZqoEwIj4ZaMqVY2k/wD2ATaLiMibLyImAhNTGcNy5zMzG2gCOnbAEOl+56WSJtJ9HLstb2F5g+dHSPcr1wQ+VWV/AHmC53yqtzC7a5GWOxe4AJgl6e1Z2hBgcPZ+QUS8kaMOZmZWVUcPGHoXacDqLsCBZekBKNvm/uaQ9yqdAzxPGi37CL0fpTSVdN+z0ijSRPM9GZm9vlZl33zg28AZvayXmdmA1+EzDP0KWAM4nPriGJA/eL4P2C0ibqznZKT10n4qaUREPAkgaTiwFXBUjbxjq6SdQfqmcBjweJX9ZmZmAKOBfSLiykYUljd4TgNWasD5zgO+AVwr6VhSM/lE4GlStywAkjYAngAmRMQEgIiYXFmYpBeAIdX2mZlZcR3c8pxJna3NcnlH2x4FHJsFtV7LZhAaRxoxO4k00cF0YFxEvFJ2qEgtyl5NH2hmZsV1+Gjbk4DvS1q5EYXlbXkeSxos9KikR1l2cE9ExDZ5CoqImcCuNY6ZQQqgtcoak+ecZmZWWydPz0d6TGVdYIakv1A9ju2bt7C8wXMR6QarmZlZf/QJYDHwMvD+KvsLPcqYK3i6hWdmNjB06qMqEbFhI8vLdU9R0ro19ufqsjUzs/bVyfc8Ja1aY3/NKWLL5R2Qc3PZxASVJ/wP4IYiJzUzs/bTycETuEHS26rtkDQS+EORwvIGz1eA30paav0zSZ8AbiQ9v2lmZv1cF4PrerWxtYHLJC01GFXS+0jrTE8tUlje4Lkj8A7gimxBUSRtSQqcvyUt9WJmZtautgO2YOk5Bd5LCpzTgJ2KFJZ3wNBcSdsDdwIXZBPr3kSaXHfPIhO1m5lZe+rkxbAj4klJnwEmS/o3cBFwO2lCnh0jYkGR8nJfpYiYIWkHYAqwJ3A9sEdELCpyQjMza08dPrctEfF3SV8g9Zh+nbRO9fbZBD6FdBs8JX21m13XATsAtwD7lrqPm70kmZmZNV8nBU9J46okB2lN6C+Q1ojevCyONWRJsvNr5D27ojIOnmZm/VgHzjB0K0uWHCspf39Vtm3okmQNfaDUzMysxaqtxtUQ3QbPiHiqWSc1M7P202kDhiJiSrPK7pyrZGZmdeuke57N1O1znpL+IenzlQ+U9nD8upLOlHRk46pnZmat0mkzDEm6TtKmBY5fXtJ3JH2t1rE9TZJwMWnx6lmSTpf0BUnvlvT/JL1N0tqStpT0LUl/AGYAGwO/yVtRMzOzJpoB3CXpbknflLSZpKV6XCUNk/Q5SRcAs4EDgL/VKrine56nZYUdmBV2OMsu2SLgDeBa4JPN7F82M7Pm6rTnPCPim5J+DnwLOB5YFQhJL5Fi19uB5Uix7J7suEvyzF/Q4z3PiHgR+BnwM0nrk6Y2GgYsDzxPWuPznoh4oxefy8zM2kyHPapCRDwBHCbpu8DHgc1ZNo7dUXSQbJEZhmYCM4sUbmZm/UenjbYtFxELSTPkNaSHtDOvkpmZFdZp3bbNlHdVFTMzs7pJ2k3SVZKekrRA0jRJJ0taJUfe5SWdKml2lvcvkrZuRb0rueVpZmZvaUHL8wjSLcCjgVnApqTBPGMlbRkRi3vIewFpiczvAU+SJne/WdLHI+Ifzax0JQdPMzMDWja37U4RMafs/RRJ80hLhI0hra+5DEkfAr4CfDUifpWlTSEtYj0B2LmZla7k4GlmZkBrBgxVBM6Se7PtOj1k3Rl4E7i8rKwuSZcBR0l6Wyuf/HDwNDOzt/TRgKFtsu3DPRyzCTA9Il6rSJ9KelZzo+znqiS9NyIerauWZXIHT0kjgN2B9UnPx5SLiDigUZUyM7OBQdI6pG7XWyPivh4OXR2YXyV9Xtn+njwi6XbgHOCaiOgqXNkyuYKnpM8BvyaNzn2ONDNDucqZh8zMrJ9p0KMqa0gqD4ITI2JitQMlrUyaoa4L2L/eE9fwVeBgUrfvc5J+CZwXEdN7U1jelueJwGRgz276q/uFdzGbgzmhr6thZtY2KqNaA4Ln3IgYXesgSSsA1wMjgG0iYlaNLPOBDaqkl1qc86rse0tEXAhcKOmDwCHAfwJHSroVOBu4vsZI36Xkfc5zBPDT/hw4zcysZ6XRtvW88pA0FLgSGA18JiIeyJFtKrChpBUr0kcBC4HHc33GiPsj4uukKfoOAdYCrgZmSjpe0lp5yskbPB8B3pHzWDMzs6okDQIuBcYBn4uIu3JmvR4YCnyxrKwhwJeAW3ox0nY48MFsuxB4EPgO8Likz9fKnLfb9kjgDEl3R8STBStoZmb9QIvmtv0fUgD8EfCqpC3K9s2KiFmSNgCeACZExASAiPi7pMtJsWgoMB04FNgQ2DPPiSUtl537EGAr4CngFOCCiJgraTVST/ZpwDU9ldXtVZJ0R0XSO4CHJT3Gsn3LERHbYGZm/VoLHlXZIdsek73KnUCabUjAYJbtHd2fFHRPIi0n9k9g+4iouf6mpJ8B+wCrATeTnhu9MSLeGvAaEfOzJcwq498yevqKsZilR9FOq1WYmZn1X62YGD4ihuc4ZgYpgFamLyB1rX6nF6feG/glcE6NEbaPkGPkb0+LYY8pXDUzM+u3WjQ9X19ZN1uWrEcRMZc0VWCPcg0YkrSPpKoDhiStLmmfPOWYmZn1kQWSPlZth6SPSFpUpLC8o21/Bby7m30bZvvNzKyfW8SQul5tbJlu4DKDKTjZT95P2tNJVyLNDmFmZv1YJy6GnT0aU4phg7L35VYgDWKaW6TcnkbbfhjYrCxpJ0nvr3LSPYDHipzUzMzaT6cFT0nHAeOztwHc2cPhvyhSdk8tz12A48pOWjmkuOR5wJPCm5l1gE4KnqRpZSG1PMeTFtOunAbwDeAh4IYiBfcUPM8ALsxO+iTwBeDvVU767/LnZMzMzNpBREwBpgBICuD8iHimEWX39KjKi8CL2Uk3BGbnGeZrZmb9Uyc/qhIRDV0VJNeAoYh4CkDSWODjpNW+nwH+EhG3N7JCZmbWN1o0PV/LZMuOnRgR07Ofe1JoXeq863muDlwBjCXNPDSfNMWRssVFd4+IHpeDMTOz9tdh9zzHAj/Pfh5Hz4+jNOVRlTOBjwJ7AVdExJvZxLy7k0Yo/Zw09ZGZmVlbiIgNy34e3siy8wbPnYAfRMT/llXkTeDSrFV6UiMrZWZmrddpj6o0U97guYjun+Wclu03M7N+rJMHDEnaH9ggIo6vsu94YHpE1JzTtiTv9HzXkhYcrWYP4Dd5T2hmZu2rg6fnO5w0L0E1zwHfKlJY3k96PXC6pN+SBg79G1iLdM9zE+BwSeNKB0fEbUUqYWZmfa/Du203AqZ2s+9hup+/vaq8wfPKbLseSxYyLXdVthVpxFLHXn0zM+uXuoA1utn3zqKF5Q2eY4sWbGZm/UuHtzzvAb4G/LrKvq8B9xYpLO8kCVOKFGpmZv1Tpw4YAn4E3CrpbuB80kQ/6wAHkhZB2bZIYYXu7kpaA9gCeAdwfUTMk7Q8sDAiFhcpy8zM2kunzTBULiKmSNqNNG/7uWW7ZgC7RsTkIuXlnWFIwE+Aw4DlSPc1PwrMI43E/RNwYpETm5lZe+nwblsi4lrgWkkbkxqBcyPi0d6UlfdRlR8A3wAmAJuz9OLY1wOfzXtCSetJulLSi5JeknS1pPUL5B8p6QpJcyUtkDRN0uF585uZ2cAWEdMi4s+9DZyQv9v2QGBCRJwsqfJryePkHOIraUXgNtJSZvuSWrAnAbdL+mBEvFoj/+gs/+SsTi8C7wFWzvk5zMysB53c8pT0AdI61duQ5mefD9xOmjz+gSJl5Q2e6wB3dbNvIbBSznIOAkYAG0fE4wCS7ifNXnQIcFp3GSUNAi4G/hARny/b5VVdzMwaoJO7bSV9lLS25wLgOuBZYG3S9LM7Sto6Iv6at7y8wfMZ4P1UD1QfAqbnLGdn4K5S4ATIloq5E9iFHoInMAYYSQqyZmbWYEFHj7Y9GXgQ+GREvFxKlLQKcGu2/9N5C8t7z/MKYLykrcrSQtJ7ge8Cl+UsZxNS5StNBUbVyPuJbLu8pLskvSnpOUlnSloh5/nNzGxg2gI4uTxwAmTvf0xaqzq3vC3P44EtgTuAp7K0K0gzDv0ZOCVnOauT+pgrzSP1P/dkWLa9HDgLOAoYTRrEtB7w+WqZJB0MHAywas5KmpkNTJ37qAq11+ts/HqeEbFA0hjgK8B2pEFCz5MeT7k0IrqKnLSXSq3kSyJifPbz5GwA0ymSRkbEw5WZImIiMBFgmFTo4piZDSSdfM8TuBs4WtKtFd22KwHfp/txPVXl/ooREYuASdmrt+ZTvYXZXYu0XGk2/N9XpN9CavluSprc18zMeqmDg+fRpCc1npJ0AzCbNGDoM8CKpHE1ueWdJGF5Uhfpu0hN29nAXyPi9SInI93b3KRK+ijgoRx5e+IZjszM6tDJ63lGxD2StgDGk3pQVyfdMmz8oyqS3kaaWegg4G0smRwhgNclnQ0cHRELc57vOuCnkkZExJPZOYYDW5HuYfbkJtLzoduRJmYo2T7b3pezDmZmNgBFxP3Abo0oq1bL8wZgHGkKvhuBmaQAuh5pVqFvk1qNn8l5vvNIMxVdK+lYUhA+EXiasrkGJW0APEGamGECQEQ8L+lk4IeSXiJNljCa9C3iovLHX8zMrLhOntu20bq9SpK+SFqKbLeIuKbKIedL2hW4XNIXIuLqWieLiFezRbNPJ907FfAH4FsR8Ur56UlrglY+SjMBeBn4T+AIUvfxqXheXTOzhuike56Sflng8IiIA/Ie3NNXjC8Dv+4mcJbOdJWkK4A9gZrBM8szE9i1xjEzWHr+3FJ6kCZS6GkyBTMz64UOHG07jvyPoDTsUZVNgWNzlHEDaX5aMzPrxwKxaHHnBM+IGN6ssnuaYeidpHuctcwE1mxMdczMzNpfTy3PFUmjW2tZCCzfmOqYmVmfCejq6pyWZ6VsQoQDgK1J63keHBGPSdoD+EdEPJK3rFrDqtaRNKLGMevmPZmZmbWvCLGoqzNH20pajzRJwrrAI6TFTlbJdo8FPkVa6jKXWlfpyjx1ouCNVjMzaz8peHZsy/NnpN7U95JWCiufn2AKaZ3P3HoKnvsXrpqZmVl72pbUTftUNid6uWdI61bn1m3wjIiLelE5MzPrr4JObnkuR5onoJpVgUILnHRm57aZmRUWIbre7NjgeT9pjoHfVdm3A/DXIoU5eJqZWUYsXtSxYeFU4EpJAP+bpY2StAtpBO7ORQrr2KtkZmYFBdCh3bYRcbWk/yQtYfnVLPliUlfuNyKiWou0Ww6eZmbWkbK5bS+MiDsAIuIcSZOAj5Mm93ke+HP54th5OXiamVkS6rSW55eAfSXNJLUyL46IJ4Bb6y24p+n5zMxsIAmgS/W92stapIkPZpDman9U0p2SDpK0aj0FO3iamdkSXXW+2khEvBIRv4qIscBw4IfAaqT1o2dLukzSDpIKx0IHTzMz63gR8XRE/FdEjAK2AH5JWrLsBuAZST8tUp6Dp5mZJUFHtTy7ExH3RMQ3SLMKnU4aPPTtImV4wJCZmSWl4NnhJG0E7APsRerOfQn4dZEyHDzNzCwJ4M2+rkRzSFoN2IMUND9G+rS/B44GfhMRrxcpz8HTzMySABb1dSUaR9JQ4LOkgLkDaX7bh4CjgEsiYnZvy3bwNDOzTvVv0qTv84CJwEURUWgO2+54wJCZmS3RggFDktaV9N+S/iLpNUkhaXjOvO+Q9HNJT0paIGm6pLMkvbPK4VNIk8EPi4hvNipwglueZmZW0roBQxsBu5NWMvkj8Ok8mZRmdb+OtKD1eOBhYBQwARgt6eMREaXjI+LzDa73Wxw8zcwsaV3wvCMi1gKQdCA5gyfwHmBL4JCImJilTZa0GDibFFSnNbqy1Th4mplZ0qLgGRGLe5l1uWz7UkX6C9m2ZbciHTzNzKy/mArcAfxQ0uPAI6Ru2/HATRHxcKsq4uBpZmZJY1qea0i6r+z9xLIu1rpEREj6DDAJuLds12+BLzbiHHk5eJqZ2RL1B8+5ETG6ATXpznmkuWm/RhowNBI4AbhS0k51dAkX4uBpZmZJm88wJGlH4MvApyLiD1nyHZKeBG4BdgKubUVd/JynmZn1Fx/ItvdWpN+TbUe2qiIOnmZmlpSm56vn1VzPZtuPVaRvnm2faXoNMu62NTOzpIWrqkjaLfvxI9l2B0lzgDkRMSU7pos0pd4B2TFXAz8CLpZ0Imm07fuA44CngWtaU3sHTzMzK2ntkmRXVLz/RbadAozJfh6cvQCIiJckbQEcDxwJvAuYDVwPHB8RrzSxvktx8DQzs6SFwTMi1JtjIuJp4IAqh7eU73mamZkV5JanmZkt0bpu237NwdPMzJLW3vPs1xw8zcwscfDMzcHTzMySNp9hqJ14wJCZmVlBbnmamVlSmmHIanLwNDOzJXzPMxcHTzMzSzxgKDff8zQzMyvILU8zM0vc8szNwdPMzBI/qpKbg6eZmSUebZubg6eZmS3hbttcPGDIzMysILc8zcws8YCh3Bw8zcws8YCh3Bw8zcws8YCh3Bw8zcwscbdtbi0fMCRpPUlXSnpR0kuSrpa0fs6860u6SNJMSQskPSrpJEkrNbveZmZmJS1teUpaEbgNeAPYl/Q95yTgdkkfjIhXe8i7EnArMBT4ITAT+ChwAvAe4EvNrb2Z2QDglmcure62PQgYAWwcEY8DSLofeAw4BDith7xbkYLkdhFxS5Z2u6TVgSMkrRgRrzWv6mZmHc4DhnJrdbftzsBdpcAJEBHTgTuBXWrkXS7bvlSR/gLpc6hBdTQzG5hKA4bqeQ0QrQ6emwAPVkmfCoyqkfdWUgv1x5JGSVpZ0jjgcOCcnrp8zczMGqnV3barA/OrpM8DVuspY0S8LukTwFWkYFtyPvCNhtXQzGyg8mjb3PrNoyqSlgcuB9YE9iYNGPoYMJ70z31oN/kOBg4GWLUlNTUz66ccPHNrdfCcT/UWZnct0nIHAGOAjSLiiSztDkkvAhMlnRMR/6zMFBETgYkAw6TobcXNzDqeBwzl1up7nlNJ9z0rjQIeqpH3A8D8ssBZck+2HVln3czMzAOGcml18LwO2ELSiFKCpOGkx1Cuq5H3WWA1SRtVpG+ebZ9pVCXNzMx60urgeR4wA7hW0i6SdgauBZ4Gzi0dJGkDSV2SxpflvRB4GbhR0r6Sxkr6HvBT4K+kx13MzKy3Svc863kNEC0NntnjJOOAR4FJwKXAdGBcRLxSdqiAweX1i4gZwBbAP0izEt1ImnRhIrBtRCxu/icwM+tgDp65tXy0bUTMBHatccwMqkx6EBEPAbs3p2ZmZgOcBwzl1m8eVTEzsybzkmS5tXxVFTMzs/7OLU8zM1tiAN23rIeDp5mZJZ5hKDcHTzMzSzxgKDff8zQzMyvILU8zM0s82jY3B08zM0t8zzM3B08zM1vCwTMXB08zM0s8YCg3DxgyMzMryC1PMzNLPGAoNwdPMzNLPGAoNwdPMzNLHDxzc/A0M7PEA4Zy84AhMzOzgtzyNDOzJTxgKBcHTzMzWyL6ugL9g7ttzczMCnLwNDMzK8jB08zMWkrSupL+W9JfJL0mKSQNL5B/HUm/lPSspDckTZd0chOrvAzf8zQzs1bbCNgd+CvwR+DTeTNmQfZOYDrwTeDfwPCszJZx8DQzs0zLHvS8IyLWApB0IAWCJ3AO8AwwNiJKlZ3S4PrV5OBpZmaZ1kwxFBGLe5NP0ruB7YB9ygJnn/A9TzMz6y+2yrYLJP0+u985X9LFkt7Ryoo4eJqZWabUbVvPizUk3Vf2OriBFRyWbX8JPArsAHwf2BG4WVLLYpq7bc3MLNOQbtu5ETG6AZWpphQcJ0fE17Ofb5P0InAZqUv3piadeykOnmZmlmn7meGfz7a/r0i/JdtuioOnmZm1VtsHz6k19vdqIFJv+J6nmZn1F3cBz5K6Z8ttn23vbVVF3PI0M7MyrVkNW9Ju2Y8fybY7SJoDzImIKdkxXcBFEXEAQER0SToKuFDSOcDVpMkRfgRMBm5rSeVx8DQzs7e0tNv2ior3v8i2U4Ax2c+Ds9dbIuIiSYtJo2z3B+YBlwA/iIiWrQnj4GlmZpnWTJIAEBHq7TERMQmY1PBKFeB7nmZmZgW55WlmZpm2H23bNhw8zcws07pu2/7OwdPMzDJueebl4GlmZhm3PPPygCEzM7OC3PI0M7OMu23zcvA0M7OMu23zcvA0M7OMW555+Z6nmZlZQW55mplZGXfb5uHgaWZmGXfb5uXgaWZmGQfPvBw8zcws49G2eXnAkJmZWUFueZqZWcbdtnk5eJqZWcbdtnk5eJqZWcYtz7wcPM3MLOOWZ14tHzAkaV1J/y3pL5JekxSShufMO0jSDyTNkPS6pH9K2rXJVTYzM1tKX4y23QjYHZgP/LFg3hOB44GzgB2Au4ArJH2mkRU0MxuYSt229bwGhr7otr0jItYCkHQg8Ok8mSStCRwBnBIRP82Sb5e0EXAKcGMzKmtmNnC42zavlrc8I2JxL7NuBywHXFKRfgnwAUkb1lUxM7MBzy3PvPrTJAmbAG8Aj1ekT822o1pbHTMzG6j602jb1YEXIiIq0ueV7Tczs15zt21e/Sl49oqkg4GDs7dvnAAP9mV92sAawNy+rkQf8zXwNSjxdYCNl/w4+2Y4fo06yxsQ17M/Bc/5wNslqaL1WWpxzquSh4iYCEwEkHRfRIxubjXbm6+BrwH4GpT4OqRrUPo5Irbvy7r0J/3pnudU4G3AuyvSS/c6H2ptdczMbKDqT8Hzd6ShXHtWpO8FPBgR01tfJTMzG4j6pNtW0m7Zjx/JtjtImgPMiYgp2TFdwEURcQBARDwn6TTgB5JeBv4GfAkYB+yc89QTG/UZ+jFfA18D8DUo8XXwNegVLTt4tQUnlbo76ZSIGFN2zEURsV9ZvsHAD4CDgLWBacCEiLiyqRU2MzMr0yfB08zMrD/rT/c8q5K0nqQrJb0o6SVJV0taP2fe5SWdKmm2pAXZZPVbN7vOjdbbayBptKSJkh7JJumfKenS/jhbUz2/BxXlHJUtVvCnZtSz2eq9DpJGSrpC0tzs/8Q0SYc3s86NVuffhPUlXZT9X1gg6VFJJ0laqdn1biQvwNF8/Tp4SloRuA14H7AvsDfwHtKct3l+2S8gdQGPBz4LzAZulvThplS4Ceq8BnuQZm46kzTR/lHAZsB9ktZrWqUbrAG/B6VyRgDHAs81o57NVu91kDQauJs0qv1A4DPAz4DBzapzo9VzDbL9twJbAz8kff7zge8Cv2xitZvBC3A0W0T02xdwOLAI2KgsbUPSFBnfqZH3Q6TpNPYvSxtCuo96XV9/thZdg3dWSdsAWEy6l9znn6/Z16CinJuBc4HJwJ/6+nO1+HdhEOlxr2v6+nP04TX4dPY34dMV6adk+Vfs689X4DoMKvv5wOxzDc+Rb03SNKgnVKT/Abi/rz9XO736dcuTNMr2roh4a77bSI+s3AnskiPvm8DlZXm7gMuA7SS9rfHVbYpeX4OImFMl7SlgDrBOg+vZTPX8HgAg6SukVvcPmlLD1qjnOowBRgKnNa12rVHPNVgu275Ukf4C6cuFGlTHpgsvwNF0/T14bkL16famUnui+E2A6RHxWpW8y5G6PfqDeq7BMiSNJH37fLjOerVSXddA0mrA6cCREVF1pqp+op7r8Ilsu7ykuyS9Kek5SWdKWqGhtWyueq7BrcBjwI8ljZK0sqRxpNbsORHxamOr2pa8AEdO/T14rk7q0680D1itjryl/f1BPddgKZKGAOeQWp4X1F+1lqn3GpwKPApc2MA69YV6rsOwbHs5cAuwLfATUpff/zaqgi3Q62sQEa+TvkQMIgWLl0ndlTcA32hsNduWF+DIqT/NbWvNdxawJbBjRFT7A9RxJP0HsA+wWZU/GANJ6Yv0JRExPvt5cvZs9SmSRkZEf+qNKEzS8qQvD2uSBhrNBD5GGlDYBRzad7WzdtPfg+d8qn+b7O7bZ2XeDbrJC91MNN+G6rkGb5F0Cmn1mX0j4pYG1a1V6rkG55Ja2bMkvT1LGwIMzt4viIg3GlTPZqvnOjyfbX9fkX4LacDMpvSPrvx6rsEBpHu/G0XEE1naHZJeBCZKOici/tmwmranXi3AMRD1927bqaQ++kqjqD1R/FRgw2xoe2XehSzb59+u6rkGAEg6Bvg+8M2ImNTAurVKPddgJPA10h+N0msrYIvs5/7U2qj3/0NPejsApdXquQYfAOaXBc6Se7LtyDrr1h94AY6c+nvwvA7YIns+D4DsQeCtsn09uR4YCnyxLO8Q0ny5t/Sj1kY91wBJ3wROAo6JiLOaVckmq+cajK3y+idp0MlYoD9N/VjPdbiJNFBku4r00hJV99E/1HMNngVWk1Q5WHDzbPtMoyrZxrwAR159/axMPS9gJVIL8QHSMPSdSX/4ngRWLjtuA9I9i/EV+S8jtS4OBD5J+kP5Oun+V59/vmZfA9IkCYtJfzi3qHiN6uvP1qrfgyrlTaZ/PudZ7/+H47L0/wI+RZo0YwFwYV9/tlZcA2A46TGVR0kTLIwFvpel3UfZs5P94QXslr3OJj3neWj2fpuyY7qACyrynZL9HfwOqRv77OzvxGf7+jO106vPK9CAX5D1gauyX/CXgd9Q8TBw9p8igOMr0lcgPdf2bPbLcjcwpq8/U6uuAWl0aXTzmtzXn6tVvwdVyuqXwbPe60B6jvE7WfBZCDwFTACG9vXnauE1GAX8Gnia9MXhUeCnwGp9/bl6cR1q/t/O3l9YkW8waaatp0i9EfcDu/X152m3lyeGNzMzK6i/3/M0MzNrOQdPMzOzghw8zczMCnLwNDMzK8jB08zMrCAHTzMzs4IcPK0tSLpc0jxJa1ekD5Z0r6TH2mlpLEnDJYWk/crS9pP01SrH7pcdO7yVdczOPUjSPyQdUZZ2fFafps1tLelbkh6Q5L8x1pH8i23t4jDSA9u/qEg/AvgIcGBELGh5rbo3G/g48NuytP2AZYJndszHszytthfwLpa9rs12LvBO0kw9Zh3HwdPaQkQ8B3wb+LykLwJIei9wPHBuREzpw+otIyLeiIi7ImJOjmPnZMf2xXzJRwAXx7KLvjdV9kXn4uz8Zh3HwdPaRkRcTJqY+ixJa5CWCpsDHFkrb1nX6NaSfiPpFUnPS/qfyu5eSe+SdLGkuZLekHS/pL0qjllb0kWS/pUdM1vSDZLWzPYv1W0raTKwDbBVlh5ZWtVuW0lDJZ0kaYakhdn2JElDy44pneMQSROyOrwg6XpJ6+a4JpuTVgqpuZi1pO2za3ZW1tVbOvfXJJ0s6VlJL0u6RNKKkjaSdHOW53FJ1VqYlwGjJG1Z6/xm/U1/X8/TOs8hpGWR7gZGkBbmfrlA/ktIc5P+giULGa9E6lJF0krAFNKaj0eT5jDdC5gkacWImJiVM4k0efj3smPWIi0eULmEXcl/ZucenH0GSHOrduciYHfSJOx/Ii1Cfkz2mb9ScewPgD+TuoTXBH6WnWtMD+VDWhHlZdLE6N2StA9wPjAhIk7K0srPPZnU/ToK+AlpkvBNgfNI874eCvxK0n0RUb602T+y82+f1d+sc/T15Lp++VX5Ak4m3f+8qkCe/bI851SkHwMsAt6bvf9GdtyYiuNuBZ4DBmfvXyGtb9rd+YZn5exXljaZKhPKl9VtePb+/VSflPzYLP2DFeeYXHHcEVn6sBrX5Cbgzirpx2f5h5Ba9W+S7ilX+3y3VaRfnaXvVZa2Gml1juOqnOuPpCX++vz3yi+/Gvlyt621FUn/D9ib9Af6o5JWKVjEryveX0a6PfGx7P3WwDMRMbniuEtIA1xKi/7eC3xP0uGSPqCyplgDbF12zso6QOr+LXdjxfsHsu36Nc4zjNTt3Z3TgRNIK2ac380xN1W8fyTb3lxKiIj5pC8e61XJPyerh1lHcfC0dnMqqSWzI6mL8uSC+f/dzft1su3qVB/1+mzZfkiLol9HapndDzwjaXyDHr0onaOyHpV1KJlX8b408Gj5GudZvuzYar5MWvT71h6OmV/xfmEP6dXqs4C09J9ZR3HwtLYhaQxwEHBsRNwEnAQcWnDAyVrdvH8m284D1mZZa5ftJyKei4ivR8Q6wPtIa5+ewJL7mfUoBcPKeqxdsb9ez5O+iHTnk6TW602SVm7QOSutDsxtUtlmfcbB09pCNiL2PFJ36c+z5B+TBg+dL2m5nEXtXvF+D9IAl7uz91OAdSVtVXHcV0hdjw9VFhgR0yLiaFJr6/09nPsN8rWy7iirW7k9s+3kHGXk8QhpAFJ3ppIGHb2H5gXQDYFpTSjXrE85eFq7mEAa3XpgRCwGiIg3gQOBjUkDf/L4jKRTJW0r6RjgONJzjo9l+y8EHgOulnRg9ojGJGBb4IcRsUjSqtmsRt/K9n9S0pmkVtwtPZz7IeD9kr4kabSkjasdFBEPAv8HHC/puKyu40kDef4vIh6olq8X7gDeLekd3R0QEQ+TAui7gZt7cY+5W5LeDryXJV8WzDqGH1WxPidpNGmChP+qDBwRcY+knwNHSfp1LP0oRDV7Ad8lPT6xkNSafetB/Yh4VdI2pEcuTgFWIbWM9o6I0oCd14G/kbqQNyC1XKcBe0bEtT2c+8ekQH8+sDKplTumm2P3A54kPX5yLPCvLP8JNT5fEdeSPstnSY/GVBUR07Jrcjtwi6TtGnT+HUn/Btc0qDyztqGI6Os6mNUtm6zgV8B7IuLxPq5O25B0IbBuRHyqD859EzA3IvZu9bnNms0tT7POdgLwsKTREXFfq04q6cPAOGCTVp3TrJV8z9Osg0XEdFIX8ZotPvXapAkk3AtgHcndtmZmZgW55WlmZlaQg6eZmVlBDp5mZmYFOXiamZkV5OBpZmZWkIOnmZlZQf8f8inbFNX8lWUAAAAASUVORK5CYII=\n" }, "metadata": { "needs_background": "light" } } ], "source": [ "#NBVAL_IGNORE_OUTPUT\n", "from examples.seismic import Receiver\n", "\n", "# Create symbol for 101 receivers\n", "rec = Receiver(name='rec', grid=model.grid, npoint=101, time_range=time_range)\n", "\n", "# Prescribe even spacing for receivers along the x-axis\n", "rec.coordinates.data[:, 0] = np.linspace(0, model.domain_size[0], num=101)\n", "rec.coordinates.data[:, 1] = 20. # Depth is 20m\n", "\n", "# We can now show the source and receivers within our domain:\n", "# Red dot: Source location\n", "# Green dots: Receiver locations (every 4th point)\n", "plot_velocity(model, source=src.coordinates.data,\n", " receiver=rec.coordinates.data[::4, :])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Finite-difference discretization\n", "\n", "Devito is a finite-difference DSL that solves the discretized wave-equation on a Cartesian grid. The finite-difference approximation is derived from Taylor expansions of the continuous field after removing the error term.\n", "\n", "## Time discretization\n", "\n", "We only consider the second order time discretization for now. From the Taylor expansion, the second order discrete approximation of the second order time derivative is:\n", "\\begin{equation}\n", "\\begin{aligned}\n", " \\frac{d^2 u(x,t)}{dt^2} = \\frac{\\mathbf{u}(\\mathbf{x},\\mathbf{t+\\Delta t}) - 2 \\mathbf{u}(\\mathbf{x},\\mathbf{t}) + \\mathbf{u}(\\mathbf{x},\\mathbf{t-\\Delta t})}{\\mathbf{\\Delta t}^2} + O(\\mathbf{\\Delta t}^2).\n", "\\end{aligned}\n", "\\end{equation} \n", "\n", "where $\\mathbf{u}$ is the discrete wavefield, $\\mathbf{\\Delta t}$ is the discrete\n", "time-step (distance between two consecutive discrete time points) and $O(\\mathbf{\\Delta\n", " t}^2)$ is the discretization error term. The discretized approximation of the\n", "second order time derivative is then given by dropping the error term. This derivative is represented in Devito by `u.dt2` where u is a `TimeFunction` object.\n", "\n", "## Spatial discretization \n", "\n", "We define the discrete Laplacian as the sum of the second order spatial\n", "derivatives in the three dimensions:\n", "\\begin{equation}\n", "\\begin{aligned}\n", "\\Delta \\mathbf{u}(\\mathbf{x},\\mathbf{y},\\mathbf{z},\\mathbf{t})= \\sum_{j=1}^{j=\\frac{k}{2}} \\Bigg[\\alpha_j \\Bigg(&\n", "\\mathbf{u}(\\mathbf{x+jdx},\\mathbf{y},\\mathbf{z},\\mathbf{t})+\\mathbf{u}(\\mathbf{x-jdx},\\mathbf{y},\\mathbf{z},\\mathbf{t}) + \\\\\n", "&\\mathbf{u}(\\mathbf{x},\\mathbf{y+jdy},\\mathbf{z},\\mathbf{t})+\\mathbf{u}(\\mathbf{x},\\mathbf{y-jdy},\\mathbf{z}\\mathbf{t}) + \\\\\n", "&\\mathbf{u}(\\mathbf{x},\\mathbf{y},\\mathbf{z+jdz},\\mathbf{t})+\\mathbf{u}(\\mathbf{x},\\mathbf{y},\\mathbf{z-jdz},\\mathbf{t})\\Bigg) \\Bigg] + \\\\\n", "&3\\alpha_0 \\mathbf{u}(\\mathbf{x},\\mathbf{y},\\mathbf{z},\\mathbf{t}).\n", "\\end{aligned}\n", "\\end{equation}\n", "\n", "This derivative is represented in Devito by `u.laplace` where u is a `TimeFunction` object.\n", "\n", "## Wave equation\n", "\n", "With the space and time discretization defined, we can fully discretize the wave-equation with the combination of time and space discretizations and obtain the following second order in time and $k^{th}$ order in space discrete stencil to update one grid point at position $\\mathbf{x}, \\mathbf{y},\\mathbf{z}$ at time $\\mathbf{t}$, i.e.\n", "\\begin{equation}\n", "\\begin{aligned}\n", "\\mathbf{u}(\\mathbf{x},\\mathbf{y},\\mathbf{z},\\mathbf{t+\\Delta t}) = &2\\mathbf{u}(\\mathbf{x},\\mathbf{y},\\mathbf{z},\\mathbf{t}) - \\mathbf{u}(\\mathbf{x},\\mathbf{y}, \\mathbf{z},\\mathbf{t-\\Delta t}) +\\\\\n", "& \\frac{\\mathbf{\\Delta t}^2}{\\mathbf{m(\\mathbf{x},\\mathbf{y},\\mathbf{z})}} \\Big(\\Delta \\mathbf{u}(\\mathbf{x},\\mathbf{y},\\mathbf{z},\\mathbf{t}) + \\mathbf{q}(\\mathbf{x},\\mathbf{y},\\mathbf{z},\\mathbf{t}) \\Big). \n", "\\end{aligned}\n", "\\end{equation}" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "damp(x, y)*Derivative(u(t, x, y), t) - Derivative(u(t, x, y), (x, 2)) - Derivative(u(t, x, y), (y, 2)) + Derivative(u(t, x, y), (t, 2))/vp(x, y)**2" ], "text/latex": "$\\displaystyle \\operatorname{damp}{\\left(x,y \\right)} \\frac{\\partial}{\\partial t} u{\\left(t,x,y \\right)} - \\frac{\\partial^{2}}{\\partial x^{2}} u{\\left(t,x,y \\right)} - \\frac{\\partial^{2}}{\\partial y^{2}} u{\\left(t,x,y \\right)} + \\frac{\\frac{\\partial^{2}}{\\partial t^{2}} u{\\left(t,x,y \\right)}}{\\operatorname{vp}^{2}{\\left(x,y \\right)}}$" }, "metadata": {}, "execution_count": 6 } ], "source": [ "# In order to represent the wavefield u and the square slowness we need symbolic objects \n", "# corresponding to time-space-varying field (u, TimeFunction) and \n", "# space-varying field (m, Function)\n", "from devito import TimeFunction\n", "\n", "# Define the wavefield with the size of the model and the time dimension\n", "u = TimeFunction(name=\"u\", grid=model.grid, time_order=2, space_order=2)\n", "\n", "# We can now write the PDE\n", "pde = model.m * u.dt2 - u.laplace + model.damp * u.dt\n", "\n", "# The PDE representation is as on paper\n", "pde" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "Eq(u(t + dt, x, y), (-(-2.0*u(t, x, y)/dt**2 + u(t - dt, x, y)/dt**2)/vp(x, y)**2 + Derivative(u(t, x, y), (x, 2)) + Derivative(u(t, x, y), (y, 2)) + damp(x, y)*u(t, x, y)/dt)/(damp(x, y)/dt + 1/(dt**2*vp(x, y)**2)))" ], "text/latex": "$\\displaystyle u{\\left(t + dt,x,y \\right)} = \\frac{- \\frac{- \\frac{2.0 u{\\left(t,x,y \\right)}}{dt^{2}} + \\frac{u{\\left(t - dt,x,y \\right)}}{dt^{2}}}{\\operatorname{vp}^{2}{\\left(x,y \\right)}} + \\frac{\\partial^{2}}{\\partial x^{2}} u{\\left(t,x,y \\right)} + \\frac{\\partial^{2}}{\\partial y^{2}} u{\\left(t,x,y \\right)} + \\frac{\\operatorname{damp}{\\left(x,y \\right)} u{\\left(t,x,y \\right)}}{dt}}{\\frac{\\operatorname{damp}{\\left(x,y \\right)}}{dt} + \\frac{1}{dt^{2} \\operatorname{vp}^{2}{\\left(x,y \\right)}}}$" }, "metadata": {}, "execution_count": 7 } ], "source": [ "# This discrete PDE can be solved in a time-marching way updating u(t+dt) from the previous time step\n", "# Devito as a shortcut for u(t+dt) which is u.forward. We can then rewrite the PDE as \n", "# a time marching updating equation known as a stencil using customized SymPy functions\n", "from devito import Eq, solve\n", "\n", "stencil = Eq(u.forward, solve(pde, u.forward))\n", "stencil" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Source injection and receiver interpolation\n", "\n", "With a numerical scheme to solve the homogenous wave equation, we need to add the source to introduce seismic waves and to implement the measurement operator, and interpolation operator. This operation is linked to the discrete scheme and needs to be done at the proper time step. The semi-discretized in time wave equation with a source reads:\n", "\n", "\\begin{equation}\n", "\\begin{aligned}\n", "\\mathbf{u}(\\mathbf{x},\\mathbf{y},\\mathbf{z},\\mathbf{t+\\Delta t}) = &2\\mathbf{u}(\\mathbf{x},\\mathbf{y},\\mathbf{z},\\mathbf{t}) - \\mathbf{u}(\\mathbf{x},\\mathbf{y}, \\mathbf{z},\\mathbf{t-\\Delta t}) +\\\\\n", "& \\frac{\\mathbf{\\Delta t}^2}{\\mathbf{m(\\mathbf{x},\\mathbf{y},\\mathbf{z})}} \\Big(\\Delta \\mathbf{u}(\\mathbf{x},\\mathbf{y},\\mathbf{z},\\mathbf{t}) + \\mathbf{q}(\\mathbf{x},\\mathbf{y},\\mathbf{z},\\mathbf{t}) \\Big). \n", "\\end{aligned}\n", "\\end{equation}\n", "\n", "It shows that in order to update $\\mathbf{u}$ at time $\\mathbf{t+\\Delta t}$ we have to inject the value of the source term $\\mathbf{q}$ of time $\\mathbf{t}$. In Devito, it corresponds the update of $u$ at index $t+1$ (t = time implicitly) with the source of time $t$.\n", "On the receiver side, the problem is either as it only requires to record the data at the given time step $t$ for the receiver at time $time=t$.\n" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# Finally we define the source injection and receiver read function to generate the corresponding code\n", "src_term = src.inject(field=u.forward, expr=src * dt**2 / model.m)\n", "\n", "# Create interpolation expression for receivers\n", "rec_term = rec.interpolate(expr=u.forward)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Devito operator and solve\n", "After constructing all the necessary expressions for updating the wavefield, injecting the source term and interpolating onto the receiver points, we can now create the Devito operator that will generate the C code at runtime. When creating the operator, Devito's two optimization engines will log which performance optimizations have been performed:\n", "* **DSE:** The Devito Symbolics Engine will attempt to reduce the number of operations required by the kernel.\n", "* **DLE:** The Devito Loop Engine will perform various loop-level optimizations to improve runtime performance.\n", "\n", "**Note**: The argument `subs=model.spacing_map` causes the operator to substitute values for our current grid spacing into the expressions before code generation. This reduces the number of floating point operations executed by the kernel by pre-evaluating certain coefficients." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "#NBVAL_IGNORE_OUTPUT\n", "from devito import Operator\n", "\n", "op = Operator([stencil] + src_term + rec_term, subs=model.spacing_map)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can execute the create operator for a number of timesteps. We specify the number of timesteps to compute with the keyword `time` and the timestep size with `dt`." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "output_type": "stream", "name": "stderr", "text": [ "Operator `Kernel` ran in 0.01 s\n" ] }, { "output_type": "execute_result", "data": { "text/plain": [ "PerformanceSummary([(PerfKey(name='section0', rank=None),\n", " PerfEntry(time=0.004936999999999981, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n", " (PerfKey(name='section1', rank=None),\n", " PerfEntry(time=2.600000000000001e-05, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n", " (PerfKey(name='section2', rank=None),\n", " PerfEntry(time=0.0007560000000000057, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])" ] }, "metadata": {}, "execution_count": 10 } ], "source": [ "#NBVAL_IGNORE_OUTPUT\n", "op(time=time_range.num-1, dt=model.critical_dt)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After running our operator kernel, the data associated with the receiver symbol `rec.data` has now been populated due to the interpolation expression we inserted into the operator. This allows us the visualize the shot record:" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "output_type": "display_data", "data": { "text/plain": "
", "image/svg+xml": "\n\n\n \n \n \n \n 2021-04-15T19:24:51.857808\n image/svg+xml\n \n \n Matplotlib v3.4.1, https://matplotlib.org/\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n", "image/png": "iVBORw0KGgoAAAANSUhEUgAAATAAAAEKCAYAAACVGgk4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAABlYUlEQVR4nO29eZhk11Un+DsRGbnvS2VlVla5qlTWaoEsy5IlYbwBn1AL2zN4GmEwYywsTLf7azBuBpoZtTHTPfaYNmM+bGhh/BloMNvMYOEFN4wtC+OWrJJkLKtkSVWlWrIqs3KNzNgjMuLOHxHn1omb9773YsuMyHy/74vvLfGW+9677/d+59xzzyWlFEKECBGiExHZ7QKECBEiRL0ICSxEiBAdi5DAQoQI0bEICSxEiBAdi5DAQoQI0bEICSxEiBAdi5YSGBHdQ0QvENFpIvpVy/89RPQXlf+fIKKjrSxPiBAh9hZaRmBEFAXwSQA/CuBGAD9JRDcamz0AYF0pdQLAbwP4aKvKEyJEiL2HViqw2wGcVkqdVUrlAfw5gLcZ27wNwB9V5v8awFuIiFpYphAhQuwhdLXw2IcAXBTL8wDucG2jlNoiog0AEwBW5EZE9CCABwEgFou9ZnJyslVlDhEiBICFhYUVpdQUANxzzz1qZWXFbxcAwFNPPfUVpdQ9LS2cQCsJrGlQSj0M4GEAmJ2dVT//8z+/yyUKYYOtWxqvcwnrUHC3Jz70oQ+d5/mVlRU8+eSTgfaLRCI7qi5aSWCXABwWy3OVdbZt5omoC8AIgNUWlilEk+HVl1YphWKxCACIRqNWsvIjuBDtgXbtM91KH9iTAF5JRMeIqBvA/QAeMbZ5BMD/XJl/B4Cvqna9UyE0lFL65/q/VCqhWCyiUCigUCigWCyiVCp57uN1zBC7C/l8vH47jZYpsIpP6/0AvgIgCuAzSqnniOjDAE4qpR4B8IcA/oSITgNYQ5nkQrQhglZOrshbW1solUrIZrMAygorEomgq6tLLwc5V6jMdh/t/GFpqQ9MKfUlAF8y1j0k5rMA/qdWliFE46iHvIrFIvL5PDY2NgAAkUgE3d3dABCIxOQxQxLbfZRKpd0ughUd4cQPsXsIQl68DZuNW1tbyGazSKfTWFxcBADEYjH09fWhr68PQNknFomUPRh+BBX6yXYf+1KBhehM1FJZpf+DfV65XA6JRALr6+t44YUXEIlE0NfXh9HRURARenp69P5MSkHVWNBtQzQXIYGF6AjUQ16lUgmlUgn5fB65XA6pVArLy8tYWFjAuXPnAABjY2MoFAogIpRKJfT09KC7uxuRSCSwEpPnDUls57BvfWAhOge1EhcA3aq4tbWFra0tpNNppNNprK2t4fz58zh9+rTe59SpU8jlcohGo9rUBMr+sK6uLu3kB0I1tptwfRxCAgvRlqi1YpomY6lUQi6XqzIbL126hFOnTlU5fjOZDE6dOoXu7m4Ui0WMj48DgDYnTfKqRY3Vsn0Ib7juY0hgIdoK9VRI02QsFArY2tpCKpVCJpPBysoKLl++jLNnz1pbrUqlEl566SXdSsnqjYgQjUYBoC6T0ryekMyaj7AVMkTboB7VBVSbjMViEZlMBoVCARsbG9jY2MClS5dw5swZZDIZ57ESiQReeuklfayRkREQkQ6xiEajdZmUIVqH0AcWoi3QiOoyTcatrS0kk0lks1ksLS1hZWXFl7wYmUwGZ8+eBQDkcjkQEXp7e7cpMaB2k5LLXOs+IbwREliIXUMziAsA8vk8SqUS0uk0CoUC1tbWkEwmMT8/j8uXLwciL0Y6ncbZs2eRzWZBRBgaGkIkEkEsFqsiMqDs6Gfnckhku4OQwELsOOqtdCZ5SeVVKBSwubmJbDaLK1euYGNjw+nz8kM6ncbLL7+MaDSK0dFRRCIRrcRisRiAq0pMElqthBQSWeMICSzEjqJe1cVTdtRzlyA2GXO5HFZWVpBOpzE/P4/V1dWGHLzFYhHnzp3DxMQEotEoBgYGqoJdu7q6tH9MOviB+ggpjCGrHVwf2hEhge1BNMNkZGd9oVBANptFPp9HPB5HNpvF4uIiEokELl++3JTy5nI5XL58GbFYDENDQ4hGo1qJyf6TXV1dKJVKVal56jEpQxKrHaECC9FSNGIu8lQ66vP5PAqFAvL5PBKJBLLZLJaXl5FKpTA/P687aXuhp6cH3/nOd1AqlXDzzTdja2vLc/uLFy9iZGQEXV1dWon19vZCKYXu7m7EYrGq1kqJeszKkMSCo10JLBxWbQ+gGb4uDo3Y2tpCPp9HNptFKpXSwamrq6tYWFjA5cuXEY/Hfc/53ve+F6dOncK1116L66+/Ht/97nfxcz/3c577lEolrK+vY2FhAQsLC1hdXcX6+joSiQRSqZRWgrKsjeSiaufwgHZDM/OBEVGUiJ4hoi80Wq5QgXUwGnn5zKBUqbzy+TwymQw2NzeRyWSwtLSEZDKJS5cuIZfLeR73xIkTeNe73oWHHnqoav11112HP/iDP8DMzAw+97nPVXUzMnHlyhX09vYiFothYGAAANDX14fh4WGtxpRSiEaj2Nraqjv4Vd6LevfdD2gB0f9bAM8DGG70QCGBdSAaJS7galAqExdH1WcyGWQyGaRSKayvryOZTGJhYUE78F3o6urCfffdhwceeAD33Xefc7sPf/jDuP322/HpT38aX/ziF51mZTabxcLCgjYlBwcHtVLs6+uDUgpdXV2IxWJVBBQSWWvQLAIjojkA/wLAfwTwgUaPFxJYh6GZvi6lFAqFAkqlko6qTyaTSCaT2NjYwMrKCpLJJBYXFz3Pe9ttt+Guu+7Ce97zHnz/93+/b1nuu+8+zM3N4RWveAW++c1v4uTJk9btNjc3kUgkEI1GkUqldLoeNh3NUAubYz9sqWwOamiFnCQi+UAfVuVBeRj/F4BfATDUjHKFBNYhaIa5CECbjJzymaPq2ce0sbGBzc1NxONxbTq6zn3ixAkcPXoU73znO/HWt74VExMTgct0yy234KGHHsLnP/95fO5zn8PLL7+MM2fOWMu+vLyMdDoNALpltFQqobe3FwC2dT1qVrhFvfvuRdRQ/1aUUrfZ/iCi+wAsKaWeIqI3NqNcIYG1ORqV7qaTldULR9WnUilsbW1hY2MD6XQa6+vriMfj2NjYwNramvWYBw8exOTkJN7+9rfjlltuwY//+I/XVbaJiQm85z3vwfDwML797W/j85//PJaXl3HlypWq7TjbRSQS0QqsVCqhv78fwFUCi0aj2j8GNK7EgJDIgKb6wO4G8FYiuhdAL4BhIvqvSqmfrveAIYHtUdiCUtlJr5RCNpvVwan5fF53C4rH49r3ZWJkZASjo6N485vfjGPHjuGBBx7A7Oxsw2V9xzvegde97nXo7u7GuXPn8NWvfhXr6+vY3Nys2m59fV2PcLS1tYXBwUEdKyYHDIlGo1BKaSXWqBoDQiJrBoEppX4NwK8BQEWBfbAR8gJCAms7sP+lGSajbGnkpIP5fB7FYhGpVEoPupHNZnW4AmeWkOB00LfddhtmZ2dx//3347rrrsPMzEyjl6sxNzeH9773vXj++ecRi8Vw+fJlPPXUU4jH47qPZbFY1GXjhgfufmT2o+zu7tbBrwCq+lE2QkL7lcjaNdwkJLA2RLODUqXzO51OY2trC5ubm8jlclhbW9MElkwmqzpkd3d3Y2xsDCdOnMDs7CzuvfdeHD16FG984xubcZnbMDMzg5mZGRARzp07h97eXiwsLOCll15CPB5HPp8HUO5DyeYwB7vyqEdEpEmLW1p5AJFGg1/3M5pNYEqpRwE82uhxQgJrEzTL1wWgKraLTa5cLod8Pq9Nxng8jlwupxUOp8bhfoijo6OYnp7GoUOHcNddd+Ho0aO4//77NTkExWOPPQYiwutf//rA+7zpTW9CoVBAV1cXzp8/j6GhIVy+fBmLi4vY2NjQJjCrqr6+PkSjUfT09FR1P5LxYrLFMjQpa0PYFzKEE03yLWzzdbGJxd2BMpkMcrkcNjY2qggskUhov1hXVxeGh4cxODiIw4cP48Ybb8Q111yDe++9F9dcc03N5PX444/jDW94AwDgW9/6Fl772tcG3jcWi+EnfuIncPbsWfT39+Ps2bM4deoULl68iEQigUQioVP7FItFTWCRSEQTWaFQqAq3YL8YX0czEibuFyILTcgQVWgWcQFXTSUZIiG7BLHy4jCJXC6Hzc3NqpiqwcFBdHd3Y25uDnNzc7juuutw55134rrrrsOxY8dqLtvHPvYx/M3f/I1e/uAHP4gf+7Efwwc/+MHAx+ju7sb111+P7u5uvPjiixgfH8eLL76I+fl5zM/Pa2Lm+DVuqWQC6+np0ZH7pVJJ+8WAavIKicwfIYGF0GiGuchT9nUB0MTFyiubzSKTySCbzeocXhsbG9ja2kIul9OxVNFoFAcPHsTo6ChuvPFG3HDDDbj55ptx55131ly2kydP4hvf+AY+9alP6SHVgLIpeeHCBXR3d+Ouu+7CbbdZQ4WsOH78OI4fP47h4WHMzc3h1KlTGBgYQDwex8LCgiZtpVSVKcwdwXt7ezW5x2IxlEqlKgc/0JyQC2DvBsGGBBaiqeYiUO3r4tguzlXP6oQH3EgkEsjlcshms/oFjkQimJiYQH9/P66//nrMzs7iNa95DW699VbMzc3VVK54PI5vfvOb+MIXvoCnn366irwY586dw+c+9zl873vfw/LyMu68806Mjo4GPsddd92FI0eO4ODBgxgZGcGlS5cQi8V0tycOE+HwEB6HkgNfeTxKvncAmhr8ytiLaiwksH2MZpqLMjzC7A7EsV3pdBrZbBbpdFpH2KdSKf3Scifp3t5eXHPNNRgdHcXtt9+O48eP44477qgaOTsInnjiCTzxxBN45pln8Oijj1rJi/H4449jcXER2WwWL730El73utfh9ttvD3yuubk5TE5OYnR0FGfOnEF3dzfW19d1eupkMllFZJLAZAYLaUq6gl9ty7VgrxBZ6MQP0RBc/Rg5rYxUXuwPkgTGpmU0GsXg4CB6e3tx6NAhDA0N4dZbb8XBgwdx991348iRIzWV69SpU3jhhRfw+OOP49vf/jZOnz7tSV6Mc+fO4bHHHsOlS5dw+fJlXLp0CTfccAOuv/76QOft7e3FD/7gD+Lo0aPo6enB4uIiurq6sLm5qTNmsH+PSV3GwwHYlnufp3Idx+N1OgE1A6EC22do1gO3BaXK1DfFYlG3xHFwKsdzcVJCoPzS9/T0YHJyEoODg7jpppswMTGBu+++G3Nzc5ieng5cpvPnz+PChQv4+te/rlsHX3rpJWfXIxvOnDmD9fV1bGxsYGlpCd/73vfw+te/HkeOHAlMpEeOHMFb3vIWXLx4EdFoFKurq4jFYkgkElhdXUUul9MKNZvNagUGYFtiRHM0pEgkEga/CoQEto/QCpORVRe/hMViUb+UmUxGd8jO5XLacc/qobe3F4ODgxgeHsaxY8cwPj6OO+64AwcOHMCrX/3qwGXigWu/9rWv4dKlS3jmmWd0i6Ct65Ef1tbW8Nxzz2FjYwPz8/NYX1/H3Nwc3vSmN2FmZgaTk5O+x5iensb09DQikQiWlpbQ1dWFtbU1nD17FolEQqtRvmcMbqlkX6DshgRcdcY3M/i1k4ksJLB9gGYTF4CquC4mL25BlMTFvi/2/QDllzQWi2F4eBgHDx7E1NQUbr75Zhw4cAA/8AM/oDtD+yGTyWBhYQFf/epXsbi4iJMnT+LKlSt4+eWXsbq66psq2gvJZBKnT59GPB5HIpHAwYMHkUwmMT09jTe/+c2YmZnRWSe8cOuttyKZTGoi6+7uxsrKChYXF3WvAzalM5kMisWiHv1IRu9HIhHtI7ONFs7YTy2VTezM3XSEBNYENNNc5KnZHYgj6qU5xN2CeMrkxaNbc+frAwcO4MSJE5iensbtt9+OqampwOR1+fJl/OM//iOuXLmCxx9/HGtra3jhhRewsbGB9fX1plz31tYWFhcXkcvlsLi4iEwmg/HxcaTTaUxPT+P1r399oH6Xg4ODuOuuu7QSW1pawksvvYTl5WWdYYN9hwC2EZhUXDx4iJkwkaeNkpCsM51AZiGB7VE0m7y4tcd01HNEfaFQ0MqLiYsVGVBOLTM0NIS+vj4cPHgQMzMzmJubw6te9SocPHgQR48eDVSe1dVVPPnkk7h06RKeeOIJbe4lk0lcuXJF+9aaifX1da0mh4eHkc/nMT4+jmQyidnZWdx+++0YHx/3PMbg4KDOUnHlyhX09fXpxoIrV67okBIA+qMgCYyj9jlmDNju4OeI/mYpqU5QZGEr5B5Dq1QX+2nMVM/5fF63rjGBcUwXvwCDg4OIxWI4cOAARkdH8YpXvALHjx/H4cOHccstt2wzg2xIJpN47rnncPbsWTz11FNYXl7WxCVb+Ezwyx2Lxao6TpvqxVSU5n3M5/N6wNxisYjBwUHk83lMTU0hHo/j2LFjuOmmmzA4OOh5HceOHcMrXvEKfT+Gh4fR39+PeDyuCTibzepzcm59vjYupxlyYV6zbb4etDuJ7UsFRkT3APgEgCiATyulPmL8/wEAPwdgC8AygPcopc63skzNQLNbGIFqXxe3nHGIBKuvbDZbZUIy2XE0PQelHj16FNPT07jmmmtw7bXX4tChQ74vRz6fx4ULF/Dd735X9zl8/vnntYM9l8shlUpVXXs0GtWpa2KxmM5RzwRmntM0ibnXALem8vXwf5cvX0ZPTw+KxSJGRkaQSqVw9uxZLC4u4qabbtIE5UIkEsGtt96K6elpDA0NYWhoCFeuXEFXVxfS6TRWV1e3dXjn2DEZesFmpFRifPxmBb/y/WnWsZqJfekDI6IogE8C+GEA8wCeJKJHlFKnxGbPALhNKZUmol8A8H8C+IlWlalRtNLXZaZ65pc6m81q1cPExWllpK+rp6cHR44cwcjICK677jrMzc3h+PHjviYXACwvL+vo+RdeeAFnzpzBysoKLl68qLshlUolHVbQ3d2tO0/zWI1MXJxc0HzZ+ZplCAgTGRNHLperuj7ugH7p0iWsrq6iUCjg0qVLWFtbw+LiIo4dO4Zbb73Vs7WSiDA3N6dHNbp48SJisRji8TgikYju4M6merFYRCQS0QS2tbWlGxGkU9/VDcm2vBew7wgMwO0ATiulzgIAEf05gLcB0ASmlPqa2P5xAA1lZ2wlWm0yMnGZaoAJjM1I9r9Eo1H09/eju7sbMzMzGBwcxLXXXovJyUlcd911OHLkiKc6AcoEcfLkSSwsLODZZ5/FhQsXcOnSJczPzyOVSiEej+vy8WjZkUgEfX196OrqQm9vr1ZdMiuqVCY8ZR+KjajNoNNMJqPvBfdvzGQyUEphdXUVqVQKy8vLeoDdmZkZvOY1r0FfX5/zWicmJjA8PIyxsTH09PRgeXkZ0WhUt1xy+IkMUyEibeZKU9JGZM0Ofm23kIv9SGCHAFwUy/MA7vDY/gEAX7b9QUQPAngQKKc13kk088HZglJl6hseZIOVF6sSVmNsxnBQKpuMx48f1x2xDxw44NuPUSmFp59+GisrK3j66aexuLiIM2fO4MqVKzonPr/ErKqi0ajOu9Xf34+uri709PRoAuOwAzYdXf42OZwbkwUTGBM1m40moSQSCe3/29jYwOrqKhKJBKanp5FKpTA5OYnXvOY1zuuOxWK45ppr0N3djaWlJUSjUa3E0uk0VlZWtNqVz4RhIzBz2qzgV/msmnWsZpSj3dAWTnwi+mkAtwF4g+1/VR6W6WEAmJ2d3ZE72Wzi4qkMSuXYLlZcSqkqc4pfbFYD3d3d6O7uxsjICAYGBnD48GEMDQ3hhhtuwNjYGK6//nrfmKnvfe97WFlZwVNPPYXV1VU8//zzWF1d1Y5zJkxWWkxcXV1dmrik8orFYtqUlTno5Yss74NJ3LKVla+3t7dXE5mM22LzktNgp9NppNNpTExMIJPJYHJyUk+9uiUdPnwYk5OTiEajWF9fRyQSQSKRQCQSQSaT0TnTZFJIpa7mFeMWS74+Ju5WBL/Ke9esY9Vz7v3YCnkJwGGxPFdZVwUi+iEAvw7gDUop72GfdwitMBcBVDnpWVmYueqln4vjldhR3tvbi6GhIRw4cAAjIyM4ceIERkZGcMMNN2B0dNSTvObn57G0tIRnnnkGa2trePbZZ7GxsYELFy4gmUzqKH7OoxWLxbTCGhgYQFdXV5XpyC2O/PIyaZnqxIQZJiJJ3EVg3d3dVaSulNJmZj6fx/r6OvL5PEZHR5HL5TShTU1NOdVoX18frr/+esTjcQDAxsYGIpEINjc3EYvFdAQ/m+18XgBVzyUSiaBYLOr0PNFoVCvXZga/MnartXI/KrAnAbySiI6hTFz3A3in3ICIXg3gvwC4Rym11MKyBEIrVZdUHtIsYuJi05FjuthpzubZ4OAgRkdHMTY2hiNHjmB0dBTXX3+9jrJ3IR6P4/Lly3juueewvLysievll19GOp3WuebZx8NJ/3p7e9Hb24vu7m709/drJcYJA6W/S5qNrExcL5l5P2TLJN+Tnp4enfqGl2VWWSb8UqmkU2TzMGuFQkET2eTkJG666SbMzs5a0/b09vbi4MGDUEphY2MDRIR4PK4zXMTjcSSTyaqxKFmdMllL0ubrkgqsmcGv8h7KY+8E9h2BKaW2iOj9AL6CchjFZ5RSzxHRhwGcVEo9AuBjAAYB/FXlYVxQSr21VWXyKW/Tj2WqDSYu2QonW95ssUdDQ0Po7+/HxMQEpqenMTExgePHj2NkZMQzzXOhUMD58+dx9uxZXL58GadOnUI8Hsfp06e1v4cVD5unJnH19fUhFotVEZc0F1mB8AsN+KdpNu+NSeySuDj5YHd3NwqFAnp7e3U4CfsI+Rri8bge7KO/vx9bW1sYGxtDOp3GzMwMrrnmGmfDBve7ZCKLRqMYHR3FwsIC1tfXkU6ndfCrbEU1CYyDW1mRyf+bHfwq7+VOENm+IzAAUEp9CcCXjHUPifkfavH5fR9uq31dQHXLm8yaKqPo2URhAhgYGEAsFsPExARGRkYwPT2Nubk5TWCDg4NO8rp48SIWFhZw7tw5nD59GktLSzhz5gxSqRSuXLmCfD6vVZd00LOZ2tPTg56enqrBMsxRf1hx2VRX0HtuU6ccusA+J47DYiWWy+XQ29urzTvp6M/lctqxz/cwn8/j0qVLiMfjWF1dxczMDA4fPrytTLFYDMeOHUMymQQRYW1tDb29vejv79ekls/n9QjhTGSywYKfh3wurQx+lfez1SS2Lwlst+H1UFtFXEC1r4tfSLPFjX1c0rfCmSO6urowNjaG/v5+3Ql7enpax3m54p6uXLmCeDyOF198EZcuXcKFCxdw7tw5nXq5UCjoF5DjtqSZGIvFtAJjH5iMrpfEJV9cqcCCvkjynskQC1afPJUEFovFdDxcT0+PHrlI9krgD8Ty8jLW19ehlMKVK1eQTCaxtraGubk5pNNpjI2N4cCBA1Vl4msvlUoYHx9Hd3c3BgYGsLy8rINfOUaMo/hZwQLVfj1pTvJ943Wt9I21gsz2qxO/bbETvi6pJGRHbBmBLrsBsQk3PDyMnp4eHDx4EMPDwzh06BCmp6cxNTWF2dlZq6Oex3R84YUXsLKygtOnT+Py5ctYWlrSnaM53Q2bf0yUTFwcU9bT06NVmYztMk0lqbrqURXyvrGJJe8fcHWEbSayrq4urcTYnGRS436hTGRMMIuLi9rPF4/HsbKygmQyiampKVx77bUYGhra5h/jrkeRSET3rezu7sbm5qaOGYvH41XpvGW3KElgfC0mwbfSN9YKhAqsTbATJqPsiC2j6pnAuOIDV5vg2WScnp5Gf38/5ubmMDY2htnZWd3qaJJXMplEOp3Giy++iM3NTZ1rnlscNzc3kUgk9Lmkk56Ji8/LqkuGRpgDwtZjLpqBrAxbiIW8l1KpRKNRlEol7UPiwNlCoaCnTGSsyFjhckptjiNLJBJIp9OYnJxELpfD6OgoTpw4gf7+/qr+lb29vZidncXAwIA2rznkIpVKASj7GVOpVBWJEVEVgcl+lCbZtyrzayt8YyGB7SKaffP9fDhsFkrlJf1fwNXOz6x4xsfH0dfXp1M9Hz58GGNjY5iamsL4+HhVZeSsq+fOndPEFY/HdZbT1dVVbG5uauKMRCKasNgxzz40Ji4mNnY+m62LwHZ/Ti39AM3YKPNeAtA+JbnM52GlxkosFotha2tLE5c55fTabL5zayX7zVZWVnToRaFQwPDwMI4ePYqBgQGdrbWvr6/q/oyMjCASiejo/Uwmg0gkojvas1/MbKCQLZMm8duIrR2JLCSwXUArbrpfeIStpZGJTL6QHMk+NDSkv/b9/f04cuQIhoaGMDs7qweZ5QpYKpWQSCSwsLCAjY0NvPDCC9jc3NQKbHFxUSsMHpGHg1/5ReSYroGBAU2gkrhMJ72s/LYWxkZfDtlX0sxgYW7DRLa1taXNSZPIWIFJf5ls8WVzkgltZGQEW1tbGBkZ0cszMzMYGhrS1z86OqrvGQe9EhFSqRSICNlsVo/6xM+e/XF8LayCpao07ytfa7P9WI0ej+t5O2JPE1gz4eWoN7sDmVNJXmzG9ff3o7+/H5OTk+jv78fhw4d1dP3AwIBWZHzORCKB9fV1rK2t4eWXX9bO+mQyiYsXLyKVSumXiM/Djnh2RkviksGoZkiEqbi8zMVWtH6Z5zOVRCQS0QkH+RokoUki43ECZBCsHOhkfX0dxWIRQ0NDmsg4oeLY2JgmsuHhYU3w7E9kAkun09rJz1lxpRLjMks1K1UaE0yrgl/5PI2SWDtizxLYTjjqbR2T2XQ0v8AchsAqaGRkBIODg5iZmdHE1d/fj6mpKa2WlFLaz8X56FdXV3HmzBlsbm7i/PnzyGQyWFtbq0oFw+YgO7s5JIJ9XdJB74rlspGW3wtQ6wsin5HXvqxK2KnPZZRmpSQ1k8h4PEwZBCu7JCml0NfXh1KphKGhIZ1IcXZ2Vn9gBgcH0dPToz8sSindosvdkLgDPJOYmVtMEpi819ywYt6LdvKNha2QO4BWmowyKJWnMqpeEhlPGVxBBwcH0d/fj+HhYUxOTmJoaAiHDh1Cf38/Dhw4gN7eXh17xbm/lpeXdfqYCxcuYH19HefPn0c6ncbS0pJOR8MBn9yHkcmLyZCJS0bRmyERrpYyiWYqLr8wF5M8eZ1UZlK5mITGAaVM4kxiTHYcO7a2tqYDXPv7+1EsFjE6OqoDWMfHxzE1NaUVLX9omPzS6TSi0SgSiQRisRg2NzeRTqf1GJXs9+SkiQC2qV1JznK52cGvfN/87r9rn3bDniGwnXbUs7+LnfRScUnHLRFpQhkZGdF56qempnQ3IPaFRaNRTUirq6s6C+rS0hKuXLmC+fl5JBIJLC4uVgVxcqyRjJznYFT2cTG5meMhemVQ2MmuKiaCKDIut3zB5TIrMvaJdXd3I5fLVTn6lbragX55eVkPcjs0NKR9ZtPT08jlchgYGMDExIQOQ+H4NO6In0gk9HPg4FdO08N1RroUpCKzOfqZaFv1TIISY+gDazFaZS4C9u5AZt89mQwPqG5h5Ngu7g40NjamB9ro7+/H6OiorricJiabzero8fn5eSwuLmJtbQ1XrlzRDmOu2PzCyPxgMhDVRVzyZWkX0vKDSxHaXkRezyTADRVsxnNXIBmUymNpcmYO7qK0ubmJTCaD0dFR3dl8ZGQEROU03tL85BAVfvZm8KsZFCoHEeH/ZGOGVGJ+96IeBFVjIYE1GTvZwigDFc3Eg9JEYDLgOKqhoSF0d3djYmICQ0NDmJiYqJpnRzu3LubzeVy+fBmpVArnz5/H+vo6FhYWdPAlv1RsInHAKae54VAJzibBL62M5bIRV9CXYafJzesZm2ak9JOZ2/H96urq0n6xQqGgCYz9Zhx2wpkoWFVlMhlkMhmMjY2hUChgYGAAxWKxqr/o8PAw+vr6QESawLq7u7V/LJ/Pg4iq0vNwfQKuNgrJUZDMVuBW+cbkffT6vx3RkQTWKnOR52VQqhlVL9UWO2n54XMoAishdvhOTU1p0uK0N/zV5q/+8vIyMpmMju2an59HPB7H2toaNjY2dO547lQt09tI01EmGLQRV5DYrXZRYa5yeDn/bftwQwUTOrdQssnORJZKpbSJl8/ntfJic53T9oyMjKBUKulnK2Po+Pnyc+DQC1Z53GuA65BU9bK88lp2IvgV8FZjIYG1KbxUl1RcptkouwHxyyH7Mfb09GBqagoDAwOYnJzUMV0c7c0O5I2NDaTTaZ3GmQlseXlZx3RxwCQTE8d2sfKSna45tEASVxC11S6kFQSyrC4yk619vI0M0GWHOjv65ahE3B1JfqTy+bwOyRgZGUGxWNRKrK+vT49LwD0ZuDW4v78fkUhEt1hyzjWZQpuvQ6avltcgG104MLmVpr9JVqbZ2wiI6DCAPwYwDUABeFgp9Yl6j7evCSyIychkJWO6TF8XO4jZZORQiAMHDmBgYACjo6O6yw6bE4lEAqlUCouLi0gmkzh//jxSqRTm5+d1VlA2Y9hJz+eRGVJlyyIfn8vlRVydRFheMK/DpiKkiQmgKnxEtliyr4pNTNnxnhXw1tYWNjc3USyWh3zjKZuW3O1IquFoNKpjxjKZDIhIt4bysWUUv3kN0kyWrazy+pqtxkw0UYFtAfhlpdTTRDQE4Cki+ntVPdhPYOxLAnMFpQLV/RhlTJdUXUD1kGLc6jc+Pq6Jq6+vDxMTE+jr69O+qWKxiFQqpTMjsKnIwajZbFbHdHHLFpMXO+dlameZLUKGE3hlh9grxOWCfOnlMs+zT0ySmklg0unOpMXZcpW6mg2WM1cMDAxga2sLQ0NDGB8f1920uEsSh7YAQDabBVE5ep87hnPPCWlWypALL5OylcGvEs0iMKXUAoCFynyCiJ5HefyMkMD8YPN1SdVlmopSgUmTUZpxvb292oHLg2xMTU2hp6cHw8PDOr6oUCggmUwikUggHo/rztbz8/NIp9NYXl7Wvhb2zxBdTe/MMV1mSmdJWvUGoO5FmGam6Udi9cVk5lJiHI8nxyxgkmGC4tix4eFh3Vo5NDSku4HxBweAPg4rsUwmo+P+ZJ42abaZZCxNS1u2i3rivPzQCh8YER0F8GoAT9R7jH1DYPIBmEGpZlS9VGAyTYp0BLOy6u/v17m7JicndbwXmw4AdGT26uqq7g7E/RaXlpZ03zw+l5namVu05ChAZhQ90FkhETsJ0xRj8L0ziYx9TfyRyOVy2m/GBMPOffZ1ySj+VCqF8fFxnc+fCYz7vXK/TKCsyKLRqO6GxOTJvjfTpOQy8vVI9WWGXDTTrKyBwCaJ6KRYfliVB+WpAhENAvi/AfyiUmqz3nLtCwLz8nXJjtjSx2VWHFZE3MI4ODiov7JsOo6NjekwBvZzFItFrK2tIZVKYWlpCSsrK1qBZbNZbGxsVPnV2J8lu/4wmcn8XMDVzJ8hcQWDK/TCjIiXioyd8uwXY8KTIRe8TyKRQKFQ0JH4HPw6Pj6u3Q1EhJGRER1KwaYkx+7JzK9mazgR6VZT4KovzzSDTRO50fpQoxN/RSl1m9cGRBRDmbz+VCn1/zRStj1NYKavy+aoNwlMRtMDV7/SrH44cHF4eFgPbzY2NqZ9YdySpZTSjnhWW8vLyzrlMQ8ay19aDsGQI//IUa/NHOtcQUPiqh0uIuNMsPI/Od4lExjn/ZL1Z3NzU+cG4xxk2WwWQ0NDKBaLVcGv/EEaHR3VKo4/WF1dXdq8lCnH5QfVVa9tDn7bddeDZpmQVC7EHwJ4Xin18UaPt2cJzHy4QHU/RtNMlAGFwNUgQg4W5TxRTFpMYH19fTpjAR+bm9uXlpZ0S+PGxgbW19c1qXFcEBMk+7Vsg8aa+blC4moOXEQmpzIshZ8Fd5pnxz63UhcKBcTjce03y2azVSEX3F91YGBAk1ZPTw+UUlWxY6lUSrdQAtgW/CpjxyRxyWsCqlP2NKrGmugDuxvAuwA8S0Tfrqz796o8fkbN2HME5uWodykvGX/DFZkDQbkFkeN82EE7MDCgu45wBeJhv9bW1pDJZLCwsIBkMomVlRVtVqTTaV3ZZCoYjuXiiszEacsW4RfTFaI22Hxk0mHOP64TksB4mVPocKtlqVTSg+8WCgUMDg7qlsvx8fGq1N0DAwNVA6fwPMeMcfCraSGwSmMSky2XpsOf19Xr4G9iK+Q3ADSt4u4pArOFRwDVgy1Ic5HnbeTFRDI8PIzu7m6MjY1p03FgYEA715VS29KzLCwsIJPJ6BgvJi9WerKfpNmyyOYFkxeHRwAhcbUa5n2V910SmQyClUTGz5dHE5fxY0xg/f39mshGRka0X5PrkkzpzSEXksi447nNspDr+IPHjRHSuV+rGpPHbTfsGQLzc9SbJqPpV2ATjSvQ4OAgYrGYjqofHR1FX1+fJi929LLDlUMhWHllMhmsrq5WjWFIRPr47Ntiv1lvb++2TtemuRgS187ARmTyAyfDMCSBMVnJlOKcPJFNTe5C1tfXh62tLd2Sbfo6o9EoMpkMgKvR+9yVjM1V+fE1o/ht4RXmtFYSa0d0PIHZHJpmP0bZwiiVGHD168r+B3amssnIBMajBbEZwRU1lUrpwU8XFhaQzWaxvLyMXC6nR3U248e4dZGVF/vaXH6ukLh2B6af0fSPcR1ixzmr666urqoAWI7tK5VKVT4vVmTj4+PaVcHxhQB0qyUTmCQyPqbMPSdbCk0T0uUbs12rDc3qStRsdCyBefm6zK+SJC6zhZFVEefO4sybo6OjOqZLRsGz45Yj6hOJBFZWVnSCQe7fKCsXf1XZ7yETDMqYLibTkLjaCzYiA67WHzMYVmaCZbOSA5m58aanpwcA9Cji3E+WO39zLwsA2pkviYw/ojJVtvlhtvm8pGI01Zh5rRKhAmsi/MjL9G/JKXD1q8SSnbNsclQ9ExhnE2ByKRaLOr6HWxWTySSWl5eRzWaxvr6um77Z5yBVFxOYGdO1F/xc9VbwTrpOqb5s/iQmMmkKcmskAK3INjc3dQxfb2+vzsnP/i2uj0RXBzpmkxSoJjJOmmgGv8oyy6kM3rVlJrGZlaEPrAWQtr6U0PInlRdwlbjk4K4cv8UdcVlxDQ4OVmV1YD9WPB7Xuc85QHVtbU07bmUQIX9FzZguW+tiJxJXMyq1PEYnXLvLP8aD73rFjnH0Pasm3lap8tgH3AeWHf4ynIa7K3GjEdHVbkis9phEZYYLU1nxOkleZm8O23MNCayJ8FJdpskoIfsXcswV/wYHB7UCY9MOgHaWbm5u6s7Wm5ub2NjY0DE/nEeKz8fqir+eZoJBAFp57Vfi8jpuJ9wL06yUqkzG7ZmhF9zww2E3vE8qldJqbXh4WGd+5frIJMONSACqYgXT6bS2EoDqD7tsZTcbu2z9KW0ICayJsDnqTZORtwOqg1K5EnBLoyQw9nXxg+auIoVCAaurq7plkcMiNjc3dRAjcDWanv1bMjTCDEa1OVfbGTtZgTtJlZlEJqP5ZX9VGTPGRMMmpRyKLZvNagc9t1hyXeW6xPVdZiOR3ZAAaKKUH3KbaWnO2xz8QEhgTYPZZGyajaZslpKeu+kMDQ1phz37vjjxHFcujqZeX19HNpvVjvqNjQ3dVYSlPJsL0mSUebr4KyzNxXZ/MRm7XXE7RZVJH5P0iZmOfhkzxk5+JjMmHc5DJmPGOEZM1i1upZQEJhMmAtDHNkMspEkpP6YyMyzDZs20CzqOwICrN5Rju4BqyQxUd8BmHxTHdkkCYzOPo+mZuDY2NpDL5bCysqJNx0wmo534XAnYTOCQCNktiM3ETjQXd5u4THQakZlmJZOa2TmbW6uZaGTAM0ffM4FxVzbZSCDrr4zeZyUmBxORH375fNnfZraCS7RbfWB0HIFJc9EMVOVKwg+TW/w4xobzczGB8VeMiZDjuuRgsSsrK8jlctjc3KyK7eFzcOsiVyb2sUkC7SRzsV0rKqMTiMxlVnKdAKoJjZ36rMZYmfGHUsaOcQprmeFX1jdJYESkGxFkdyfzvZG9C1xo13rRcQQGwPoQ+IvH5pwc9GJoaEgTF/u52Lzj5mkev48d8ysrK8jn84jH4zrehh+yHDiDZT0TpgyNCM3F1qGWKPLdgs2sNINguTsSt1iyucdkw8oLwDYi43osB3FhBcUKjP2z7Hszg1/5HZI9C2xo17rRcQQmI+vl19gMSuUYLqm8ONUvf7mkEzWVSiGdTmN1dRW5XE6PQMOSXsZ1yYysrpiuTjEX27ViBkGnOPvNgFGpyJg8eKxKVvmyQziTEA+51tPTg2KxqBufOBEAt3jz8YCrBMZhG/yxlsGvAHzJq13rSUsJjIjuAfAJAFEAn1ZKfcSx3Y8D+GsAr1VKnbRtIyFvJpOFJBYeWIF9XjJHOXdsZR8Dd/nh0YGYuBKJRFWmCv7CcSwXx3XJmC6Zu7ydXyhGu1bKetDuikyWTZpqTCByGDzuU8kxXZyiiesXuz36+/t1v0gZ/MofaABVYzuwMuPzMFFKInPViX3nxCeiKIBPAvhhAPMAniSiR5Qx+giVRyb5t6ghL7Y0GfnBy7TLnLuLiYs7TPMXj79AiUQCyWRSExhnSOXQCOm3kGpLmo9AZ8Z07SXyYrQ7iQHbQxbMaHj5QZQNVGx1JJNJHUuYTqe1k56DX2VGYPbN8gdbEhj312Q/mwx+taFd60srFdjtAE4rpc4CABH9OYC3YfvoI78J4KMA/l3QA8vUJjJ3PP/YdOTQCP5ysaM0kUggl8shHo8jmUwilUrpdMDc7aNSZl0ZZJ4u+bUMnfTthU528gOwEphSSkfxs4kJQI/qziE9xWIRPT09GBoaqooPA64qKHbux2Ixnc5aBtlK14xEu9abVhLYIQAXxfI8gDvkBkR0K4DDSqkvEpGTwIjoQQAPAqh6OBzbxVkkOMMlj5co/QBKKZ0pNR6P67EX0+m0TjwnHxIfn79icvRlSaCV8jXnjrUY7VoJW4FOUmOyHrFPzJV3jIf5y+VyOlCVuyAVCgUd/CozvzIpsmtFfvhly6VMcS6xb31gXiCiCICPA3i337aqPKrJwwAwMzOjOASCQyTY78Vkxq02bNtzZHM8HtfKK5PJaMc9x93IoFMmLnPsRf46doqfC9hfxCXRCWoMqCYyWa9kFD8TjFRKnKpHhvcwgXGLpXShyL69sk8u+3L5HbD5u9q1DrWSwC4BOCyW5yrrGEMAXgXg0coDOwjgESJ6q5cjnwNFBwYG0NXVpf1dTGqytz1/lTiGa21tDfl8Xvdr5FYf4Gpl4YcplVenmotA+1a8nUQnqDEAVcRlhlxIAisWi9oUlCNayUFA2NHPUzYnpfriDzan52FTMjQhy3gSwCuJ6BjKxHU/gHfyn0qpDQCTvExEjwL4oF8rZCQS0XmTZAJCNu/YIckhEKy4eMrrWYoDV/0OcrRrGRoRdgHqfHSiGjNjxkwiY9XFpp/skiRjxwBok5JJjBsCZE466Qsz0bGtkETUC+A+AK8HMAsgA+C7AL6olHrOtZ9SaouI3g/gKyiHUXxGKfUcEX0YwEml1CP1FJgJTGawZFnMrYfZbFbn7Mrn8zo0IplMamcnR9OzQ5NJUIZGdKK5CITk5YVOUmNShZkKjAmMwyg4PhIAkslkVVZXnrJpyV3rpMXBiRiB7fWnY31gRPQbKJPXoyiHOSwB6AVwLYCPVMjtl5VS37Htr8pDJX3JWPeQY9s3Bikwj+LC/Q3ZZOQHmc1mdeZLVlybm5taXsvIZ6nc+EHKDKmVe9ARFZ7RrhWtndApJAbYR0fiD69sDeeWSeBq7BdwdfBjqcg4il+2UrJPGbCrrXatV34K7FtKqf/g+O/jRHQAwJEml8kT7ANjomH5nMlkdGsi+7g4tTO3MMoO2ExW7A9gJ73ZDaiT0K6VrB3RaSTGakzGjslcY6ZPTLZUspUBQOfi53xjAKpcJYC9HrVr3fIkMKXUF811ldbDQaXUplJqCWVVtmOQcTMy7U0qldKtiuy05xZG6euSslmmIelUJz3QvpWr3dEpfjGGDJQ2O4izEmMHv4zp4pZKAHp0pHw+rwfb5f68wFXFZqJd61ggJz4R/RmA9wEoouycHyaiTyilPtbKwtnA9jg7LVOpFHK5HBKJhM7TlUwmq2JamJRkhDKbjtIMrVzrTl9S3WjXStVp6CQis8WOSdPSFvzKZCaJjPOOcfArqztWahKyq1G7IWgr5I1KqU0i+ikAXwbwqwCeArArBMZZK7k7EKd15syWPKgGUC2PZRCfzJbZaeZiSFytQacSmXR5sCUhR0rilknOdsFdirgvJGe3kOOUmmjXOheUwGJEFAPwdgC/q5QqENGuXBF3BeIWRXbYc8sjB/jJaP0gxMXTdn1QjHYv315AJ/jHzO5IrMRkahwZDymDXzk+kklNEpgrL1i71rugBPZfAJwD8M8AHiOiVwDYbFWhvFAsFhGPx5FIJPQAnzI1L/sHbMRl5kyywRXI1w5o13LtRbQziZnkJacykNsMguWMFNy5m01KNie7u7v3JoEppX4HwO/wMhFdAPCmVhXKC8ViEevr63oUF/7xTeeAVJmTnlssXeai+XDajcTaqSz7Ce1oUtrK4lJjLiKTaazZD1YoFPaeAiOinwbwZ0qpKg+eKl/NFhFdA2BGKfWNFpaxCsViUae+4VxdlbJWtSzKRIPsFzDJS+7brg+oXcu1n9CJaozLLNNIAVczswKoGmeyUCjoRgDbB75d66GfApsA8AwRPYWy034Z5UDWEwDeAGAFZYf+joFbHjlYTzYhc0wXtzZyTJdrqChJXCaJ7TaptWuF2a9oBzVmnttGXnJZrmM1JgO0uVcKO/pZldnQka2QSqlPENHvAngzgLsBfB/KXYmeB/AupdSF1hdxW5l0UzB/YaTqsuUH90I7klhIXu2L3VJjtZCXCe5XyfnDzAFv2YyUg+6aaNc66esDU0oVAfx95dcW4OA72bNeBqUGIS6JdiKxdq0oIa5ip0msEfIy95Pd5ICr18KxXrZ8YLxdO6LjBvUAUEVWsj+jmaLX5px3oR1IrF0rSYjt2CkSa4S8uD6ZPizZZUgSmDmivTxOu9bNjiMwTjbIpOUiLqWUls5mq0qQEAozLqzVcWLtWkFCuNFKEvMiLtuyWS45lel4JORIWkRkdeCbx2w3dByBcVcgGRohTUZ+SGY0ctAK4Ke2mq3G2rVihAiGVjj3g/htvcpjKi/zJ98R2QeYxYDruM0CBRytLAiC9oWcBvCfAMwqpX6UiG4EcKdS6g/rPXG94GwUrkh64CqJmd0sbDFgrrgaL+XVLBILyWvvoFlqrF7l5Udatv6MfCy2YLhV31Yvm9UKSQFHKwsK91ji1fgsyokJZyvLLwL4xXpO2CiIaJuz3mXecdcJfnjmw5TLtvPUslwrQvLae2j0mdZT57zqNC/LjCxmWWXgq0wx5TqH3y8A9GhlSqk8AB6trC4EJbBJpdRfAigB5WyrKGem2HEwgckRguR/lfJV7SNJTKYZ8XsAQSpUvUTWroGRIepHI3WhVvLyUlmyjptxXTaTl81H2ZpvogYCmySik+L3oHEo22hlhwLcJiuC+sBSRDQBQAEAEb0OwEa9J20ULHfZQS+bgSvlA1At67lzq82pz7CZjkGc+fWalLsdLBuieWjWhywIccl5m+loc9rzu2ASl5+LxXZeH6wopW4LunGjCEpgHwDwCIBriOifAEwBeEfLSuUBlrs2m9y23iQxeRzbQ3MRlq0crfCLhdg/qKdlUS4H9XXZ6qWLvFwf9ybWbb/RympC0M7cTxPRGwBcB4AAvKCUske87RD4hsuWRgA6KyVvA2x3sLIaMyHVWT3Kq55Qi5D4Oh+1qi9Xw5Ft2UZcwNWPsR952Y4tuxNxfXcRF5+jiV2JPEcrqxVBWyGjAO4FcLSyz49UXryP13viRmDKYCYefojckmKalBKulkoz5MKPlGwEVCsphSTWuWgFeQFu4gqqumxwqS5Xa75XeeqFcoxWVu/xgpqQfwsgC+BZVBz57QS/SuH1cG1qzOYn86qozfCLhSTWeWiUvILsL4kLsCsvP/IyVVZQv5etHM2AsoxWVi+CEticUur7mnHCZsHPNyVVmctnxpAOfpsS42P6PehGTcqQxDoHtZBXUNVlwuWkl+uDkJeX2vJTXrIs7YigYRRfJqIfaWlJGoRLFvMDdEUYM7zCK/xixsxyBFnndx0h2hO1Pp966kOQcAWZVcIFGWXvZTYGQQ1hFDuKoArscQD/L5WHVCug7MhXSqnhlpXMA1INuRz1EmaUcVCT0mVK+qmx0KTcm9hpk9FGDkFNRi/FZfPxmuc317VrXQxKYB8HcCeAZ1WbXIlXMYKQC+CW35LEbP0oazlXaFLuDeyE6jKntZqMXr4uuWw7rx86MqGhwEUA320X8goC29dPtjz6+cZs20ofGR8ziBprpJUyJLHdx06Sl4uwalVdvE6evxHXRLvWwaAEdhbAo0T0ZQA5XrlbYRRBHOrA9gfGJGSuA7zVmOt4fucz19erxmpVbiGag3ZQXUFbGL18XLX6umpZv9sISmAvV37dld+ug29oUCIzSU+2VkpfVxCzks/vcoT6tZD6rfO7jhCtRyvJy4u45HIQc7HR1kWzTK7/2rXeBY3E/41WF6ReBCWyIOYeYO+OxDDDLeo5T6MO/pDEWo9GyKtWc1HO1+LrsjnmTdKqpdGgWdvtNPyGVftdpdT7iehvUenILaGUemvLSuaAn8Stx+Eu1RjgT2K8jdd5/MrTiIM/NClbg51QXTzvIrBGWxiDXoef4qpl+92EnwL7GQDvB/BbO1CWwPAih6BEFkQlAf4mpZ8pGUSt+a3zug6gfStXp6BW53YzTcZmtjAGuZZaiYvRqa2QZwBAKfX1HShLzfAjskZIzEy9Y3uAksTkeb3S9bjWheEWu4Od9nUB27sD1errsi0HuZZ6HfSd7AObIqIPuP7crVZIEy4iC2LG8XayYtgelsustJmUrhz8rVZj7VrJ2hU7RV5ezvqd8HU1o2WxXeuWH4FFAQyiHHnf9miEyJhcbGQmfWOAtxoLCi+T02+d3zHbtbK1C1ppMrrMRTnvZzL6tTDapl5lCbreC+1ap/wIbEEp9eF6D04BRh8hon8J4EMoNxL8s1Kq7txAjHqILEilZjJxEZmfg98kStt5Q5OytWiV6jLdCDy1EVhQ4pLzu0Vczdi3lfAjsLqVFwUYfYSIXgng1wDcrZRaJ6ID9Z7PBpfJ5kdkfv4zPoaXWemnxkxScpmcfuuCHn+/o1Wqy7y/LtLi+VpUl23Z61rqddD7IUjr6G7Bj8De0sCx9egjAEBEPPqIHD7pvQA+qZRaBwCl1FID57OiHke/l0Iyt6mFxGzns/nhbOXwWueF/U5ktRKXax8/8vJTXYC3yRhEdXldTytUVyuO02x4EphSaq2BY9tGH7nD2OZaAKBynv0ogA8ppf7OPBCVRzZ5EABGRkbqKkytaqwWYgliUtbbKtoMk9JW5v2AnfZ18XQnfV2tUl2tPFYzsdsjc3cBeCWAN6Kc3P8xIrpZKRWXGymlHgbwMADMzs7WfSf91Jjtv6AmJeBNZLZ+mCb8fGM2NSbL7of9RGKN+rpc603V5UVggJ28bJ2sa/V17RRxtfKYzUArCSzI6CPzAJ5Q5QFCXiaiF1EmtCdbWK6aicwkFq+WSsBOZCaJeZmvQcph7hP6xspohuqyrfcyF23zwHby8iIuc70LO01cO3HsRhC83b926NFHiKgb5dFHHjG2+RuU1ReIaBJlk/JsC8tUBZN4zP8kXD4K1w/Y3hops77K89t+XuWwoVkvbidjN8jL5vOyZUv1yxrhct67ymGubzV5Bf3tNFqmwJRj9BEi+jCAk0qpRyr//QgRnUJ5pO9/p5RabVWZPMoa2MnvpcZs8FJjQWAzK10qqh6Tspbt2xU7TVw8bUR12ea9ytdqJ70fOrUVsiEoy+gjSqmHxLxCedBcZ7T/TsFlVgYxKW2Q2yi1PXYsiF/MJEgbkTVqUtazfTuhFeRl3gsvk1FuK19y+YHyMhVtH0jbuU3s9PNq1/qx2078tkNQIvPyVXkdW4ZdmPn3/WA7pxeJBT0ub9+uldSFZjjq5X9exMVT2zxgJ69afF3tSly7fV4/hATmgJdZCdhVUC0kBqDKV1KrGpNl8cJedPA3299nu956HPU21SXnXU76oM90N8mrXetESGAe8PJvyf9sysjPN+ZSY7WSTdBt27UCthq1kl2zVZdt3lW2diSvdjm/CyGB+cBLXfmpMT8SA7BNjQUxKYMe29zH77jm8dsRQa+5nlbWoKoL2N7v1cu/5Ude7UxcjHYph4mQwAIiKJHValLy/rZWSSYpv8oTqrFqNFt1ySnDy1x0rTPL1wnExdiXrZB7EUGIzERQ35jfej/yqUWdtNsL0izUQl4uknKZjK7zBW1hNM9b63+7hdAHtgdRqwkHNI/IzGNKpVaP6bQfUY/qkqilhbHWMrUj2rVsHUlgtZpo7VSOoITj1TLmd+wgZdqLKqzWD0QQX5frHC4CC1oOV5l2E52mDIEOJTCGV/zTbpSj2Q7mIGrMRoL1+OGCnLcdUK/KNO8ZT4MSF1C7r6uWMu00aj13u9aJjiYwG3bTjKqnZdC2j1zn1fJlOvmDnN/PgRykora6MgeJ23L5IF3X56W8zP9tmXX9wiL2OnmFTvwdhPmAdpLQ6lFjfiTGyxw35gqedO0ry1OPmbDTL1uQ5+dFYrZlPwd9kD6MPG3EZNwN4mrGOUMFtovYDVVmth56wfRHuchMkhjDJLNmXutuEZorwLOWwE9zGxtxBUk0yFO/kIggZdgNNOvcIYHtMjpFlZlOfrnOD7YMsK4Gg1rCQHbCuesV2OlSj37X4FJftZhDtnvWqA+zlWjVOXfiWojoYwB+DEAe5TFpf1YZyU1N7BsCMxHEYbvT55TnNv05QUxAhkuV2Y6z20TmdV1+xGUjK3O5HtLicwRRYLttXu/UOXfovH8P4NdUORXXR1Ee8Od/8dph3xKYH1rxwJplbphmpNc+fr4dG5HVQ2D1OnlrGUvTBhdxmeXxKru8djmYrNf52t10bDZ24lqUUv9NLD4O4B1++4QEtoOopRK4tvUyJ23rzRfORWi1KrtmwW9kahtsZXURVlDi8oN5H/cSOQXBLrRCvgfAX/htFBJYC9DMym3z4wR5WV0vnGkKNbti1mIm+6GWsnmFRchlPyc8d6h3+Q+b8Wx3OwC7HtRw3ZNEdFIsP6zKg/IAAIjoHwActOz360qpz1e2+XUAWwD+1O9kIYHViaAPtNEK7+Xnsa3zepFtAa6NxDM1cm1ePq9a9ve63iDqy6uFM2jjSRAEMdFt5WgH1OgDW1FK3eZxrB/y2pmI3g3gPgBvUQFOGhJYANTSXN+K/RlB/Tqul9hsiXTt66dSGrmeevxLLpKW834KzNYgYh7TVKXShK3neuQ2QfdvlNhbhR1qhbwHwK8AeINSKh1kn5DADPg9qFqJyrX9TvgUvMpjM48Ad2CtRLNfukaVVxB/VxDyNlHPMwoaWAwEJ7Z2ILUd8vn9LoAeAH9fub7HlVLv89ph3xNYPYQVdF0rSMqvvLWoJdO343rJ/NDql8n1ogc1F/kYfmayaW7Xg1oaJWohNr/n2upnsEMf3BO17rMvCawWFRWErII83FqJslkKxs/57Ip1MkML6ilXI5DlLZVK2xSJH/EyvAJRg5izfuv97of83y/yvxZ4EVqzn1GNPrAdxb4gsKCV0bauVrKqxYwMehw/n5S53kuZ+LXG8VTGQ/m9+I2+MH4NEfK6mMxcTvYgrY1B/IBey0H/4//9lK2tTjVCaq1QZyGB7SB2irCCVvQgD9+vktv+t1VyL8VVa98/r+VWwvRvyWU58AlPa1E2Qe6Pl9pzlbVW+DWWuOL1XPAyO5tBaCGBtRhBycNrOShZNasi16oSzOHp5XZ++9biOzMJi89da+hFUH+cV0ui2ePANGdq9Qma+9jIURKjrYFD7hskPMJrvW1/27OsRaXZiMtVzqBkFhJYC1CPgz0IYQUxweqBi2RcZOXnn3L5rHj/RssZVInVYp4xXCTCU1N5mcv1whaK4jp/NBrd9r/c3yQKG8H5wVVfbcc2z+vVq8KLeG1lN88RpJztgI4ksFpUVSsIK6gk93IgM0xV5SInP0XkRyz1wI/ImnEuF3nI+WYQl9e5gig+v31N89Pmowuq0GyEY27fTELzMpF5fZjQsMkIQlp+hNWIWejls/IiKxcR2MxDcx8XqTVCLPWYlrJsrvP6wUZYNvKop5x+5zOPbxKVbb2rbLzetp+t9dQ2NbfxU2XmvuZx/AhNwkVmtvO3IzqSwLwIyEZaQQgryAPyMwFdZORHPFJN2fat1aQzj+MFP7KXCJLoz4/IXC+xzWyqp0wuuEiSp37Kz6W4JLkFIUJzXzl1ldm1bNs+KKEFIbNazr1b2BMEVkvnZnPeD0FeUFNZceWwTf22CUJc8rxBy+gC3wtbVgrbfWrULA0CfpHMl8zroxH0Wtmvxcty6vKPBSEycyqJzWsbm7rjqSSUWkjOZcbKe+VHZn7naRd0HIHZHjqvl1PXOj/4kYEkDS8SMk0/8+ciLC8VZXtZXfOudUHuhblNEN9KMxSY3/kYQUxl23V73QsmN6+yefnH5Da2/2S99SI/P5UWpE6b1+V1L4KQmcusbQd0HIEBbjNATs15F4IQlh8J8Tbmst+Ptw1SDnM+KLzMENu9C/JxsMVLucrG11dPA4rN9OFj8j2UKqPW++NXfhvB+ZGbF3mZ95cJy1z22sd8Vq5757oOv3sR5Dm1EzqSwID6vkwMLx+OST42kpL/eW3rIirbec35Wu5BPctBzCa/j0MQxSfNGa+y2Y5pI3MXWdl8YX6qrJblIATHU5uJ6UVgpjqT817b2s5rlstVXhdcdTBshWwiXA/L9oDMSmgz00wS8iIuLwLj47Ip4mXq+JGV67r8iLvWRgyXCpPz9XwkGoGLuLw+PF4fBNu9rlf5BiU11301yScIgblUm+2YvI08tznPy/wxCIJ9qcConN/nEwCiAD6tlPqI8f8RAH8EYLSyza8qpb5U4zn0zQ1iDppTUz1Fo1ErYdnIyuW8D0pSjCBmHc8HbUXzmgbZJsg5beWvB7aPjESQXgBBCCwIMbrOWSvBBSE3mzPfRWrmtFgs6mPYlJ15fNvUVn4bZN1qN7SMwIgoCuCTAH4YwDyAJ4noEaXUKbHZ/wrgL5VSv0dENwL4EoCjAY+/bd5V8f3MwUgk4iQsP7+W1/klzBfe6wsNuH0qtm1dy+bxGUG6TPkRqW1aq5kRRAUREYrFYk332pY9w0VGQZdtHylXHQiq2MyWX/NZeykvOfVTazaVJs8bBPuOwADcDuC0UuosABDRnwN4GwBJYArAcGV+BMDlIAc2v5RynYuwALu6siksuW0thGWT6X4EVYvz10VUQbvHeJXV7zrktB7/iwlJTC6l7LrPXuTl2tZFMHw+L+Ly+8ky+8UAmueV8Pqo8Y+Vl02heak0nrrUej1ZVtoBrSSwQwAuiuV5AHcY23wIwH8jon8DYACANV82ET0I4EEAGBkZQVdXl2dl81JTJlm5VJpZAb0QlKRcJoJrW9s0qNLyIpRGHLKuc7vI1AWvHgZScdUDrywUfmrPnPeK6TOntnpo29al0rzUJd/frq7yKxukAYAJzDQ3Xa2dfB4bQie+HT8J4LNKqf9MRHcC+BMiepVSqupuqfKoJg8DwKFDhxQTmB8psTOdp17+LFcl8lIjrp/rK+mqcObL76W2zDIAzW/6DqosXUpQbs/z5gtbKpU8synIY/iVxwTff7/rYniZsq4PJe/n+ni6lgG3urfVQ7OMQVSaaUpKInORnVxnwrW+HdBKArsE4LBYnqusk3gAwD0AoJT670TUC2ASwJLroESEWCzmNAdd/izALe9tMB+Yl1/Bz0/hqlguApT/yXPbyuaqcH5wXbtJHn7HdikvF/HL80oSk+uV8k9T4ypjPddeKm1PmyNhU05+5qSr/kmFKdd7+WttZXFdu1lvXPVP/myqzXbsdkQrCexJAK8komMoE9f9AN5pbHMBwFsAfJaIbgDQC2DZ66BEhO7ubkSj0W2kZRKYq4KZCPo1kw/fJDAXcXkFJ7pIykVQ9b60LpPE9V8tkC+tvAYbEdqIIaiZ7odaSdcsl21/LpcXwXn1qHApMa9GJHPZnJrnMstjXqerTrt+IYFVoJTaIqL3A/gKyiESn1FKPUdEHwZwUin1CIBfBvAHRPRLKDv036187hQRoaenRxOYSWS2L6GlbACqScP1JWLZLeW37UsmVYhtmc9rElWtZNUs2O5LM4iEj2OqqFadwzxuPfetlg+DzRTmeZPYZD2U62zLJmlJNwhPAXtDlDynWUZ5LbawDUlcPK1X1e8GWuoDU+WYri8Z6x4S86cA3F3LMdmE7Orq0gTGlUfO87aV8+ipKa9NEuIHuLW1BQD6wdrIyqv5Wk5t55frjfsT6B4EvVeN7C8RhFylCvM7h8sUYuKrVyF6Kb+gL2EQNeenMKVqsxEZT11KzLQomMDY/yuJjH+mSe73HnBd5/lSqaTrfS33ZDex2078mkFE6Orq0gRmPnjeBrA/MABVZCW/QHKdqbx4vZepKM9lnp/RSEVolLj8/vMrW5DrMU1JG2zhAy5CqOd++ZmEQU3uWsxSc53N5LSpNT/TUiow9p+Z6yXJSRLjbVxlYgXJZOvVGBS2QjYRNtPQVDe8ztW8bLP9XcQVlLSaTViMWlRIs0w0L7+R3zV5OZqDno/XmdfTjOsLSozNMktdqtI0tW1kxg5/k8hcLhRJbkxMpg/N63q9lHE7oiMJzPySmY5Hm5loEphJTjalZds3iNJq1sOu19QLorJsHwDXtq7/5fkaIRazzK5z+ZltXvvWuo3tHI08V7/7Z86baozJaGtry6nATL+ZqzXU5SPzUlohgTUJSikUCgXPh2ASmBcZ+YU82B6sS20xmlHhdwp+xOW1bBKXa951LJNMbeRqvuC2/73O107wuiemQpMfSklgbKLLUAyb/0xOpQrj47s+Ol5E247oSALL5XLb/Fsmcdn+cznVXWrK64tkyn+bigjyUvn5WepxYtezvpYKavoabev8CMx8aeuFn4np+sDUel6/7ev1T/qRPvu0vMxQk5BsJqnpd/PbV0K+H+2GjiOwUqmEdDrtjHaX63je3L9R2B5ykJfCS621ojxe27pMMT/FY9vH7wXwOn8tL0cQhedXdi8VF/TcQeFHrkH9e/XUWdOXZXtWrl4CNoQE1iRIApN+LJuZxzArhhkv5pr3g9d2JqEF9fPUc74gL4JtG5tpGJSA5Lxc5+Vn4XMA1X5M28cmCNkHUTM2og5675vxjLxIt96Pl/mB9rMaXO+CK+7MhrAVskkolUpIpVLW2CzzRXB9dczWHmB7f0mvl9NWJnPbZiMoKTV6zFrPaRJWUBPSJBgz2r1ZqsjrGlxlkdvV6/C3LbvK61evZNlMknLlBfNqbDLfC1tPFhOhAmsSSqUS8vk8tra2nATGMF8ms3WGvzimj0ESm3xBZYyNhLnORmjmy9kMk8SvHLbzmjC7//jBj7DqMSFt62v94jd67eb2vM9OPSdX2cz74WqIMi0Rtk7kMWzvh0lgnO1CIvSBNRFKlVshC4UCAFgJzCQanpdxMSZxmeQmR6iR67kieZlKpgwP8vL4we8FdZXF/N8si4t8XT4UOW+a27W87PJjYZuaIwTVUr4gy7Z1O/WcbLARt0lUwPaMEmYXN5tLRR7TFq0vVZeLqEICayLk18kWWCq34XlJVlJNyW4Yksh4W3NfSWQ2Egv6EtcT7GmiEZPV5cy2+UCCqq16y+NSZbaWSvPD4lfWIOdl2F7SWp9T0Htg8/fZ6rX5gXb1zTW3c/nD5Hkkacn67lfedkJHEpispEwusrJ7PTSg+iWQRCaJSXbZkNsopbb1PTMVSRBC86oQjfpdbPtLwpIvv83B61JbXuvqQVDFaJKYjbxc97wW5RX0ery2C/K8bURt82FJUjK7srnmAXc+NIYtvMLvmYZO/CbBVEm2L48kNWD7zWdyMr94fGypuvicvJ2EVG62SmCSWRB4kU/Q/W2+H0kC5rwJL6XVLNXlVW5bWW3mrw1BiawW1RyEBP1gOtP5eiTx2ExBs4eIaToGUVyuuC/ThOR5r7K3GzqOwIDqTqou2WwjNGD7l5BVl0lA8jimqVgqlaqGTjMrJq/neVkpgiqyIM51cx+XsrLBJAQvhdUq0nId02Ue2uC6hiDnqee6/HxvtnLJFlZZV8xfsVjcRlxbW1vbTEazjsvj2hpZbIRlW8f7265pJwmMiH4ZwG8BmFJKrXht23EERlTORiHVle0rJtfZHrbrgZiqzfbltzUQSMKyqRs/UjDPL0lXrjfL5jqGTbkEIS1XGVtBWi7US2auY7juveu4rg9OkPth+6jIeZOwTNXF85zWxtai6GX2mx9iGzlJkrL1k9xNAiOiwwB+BOVkp77oWALjl5iHe7fJaBeRuQJe5TnM//krurW1VTXEfLFY1Kl9uKOtPAYrPJuvzHV9trJJ8pXnN5VaLRXNpUD8XtSgKq9W2Movz2WWxctcDnIuiSDxbLWqQlMpmf4sk6zMjCg2kzCIKjZJybUsf+Z9cF3TDuC3AfwKgM8H2bjjCAwok5Z0xJsVRq43v3i2/3idnEqY5qisYJK42H/GhCUzBthiyFxEJtUTn99UVDxlMjNVmZcKM89Vy3yr4fr6exGZ13ovhQLUNmguz/vVERthmcRlIzKbQ16WM8h1yevzGhow6AdVXl9ATBLRSbH8sCoPyuMLInobgEtKqX8OWuc6jsCYEEzYiMhUKTaT0vzSmcTmOo+5vUlcvMwEY4ZiyBZNvi4vAnEREpdBNkq49nHdT9u8bdl2H5oJPxJq1JwMorqCmJmyPLY6ZrYQSvPQRmQuv5Ysozy3jYjkzwzCttUxr+szYSNUD6wopW5z/UlE/wDgoOWvXwfw71E2HwNjzxCYDTay4WXbl9I2L/czlYD8n4nCRmRMWBw8K9fxNcnrM9cFuSdAdatoEPIy93ctm2ilORGEbE01Wcs1uhQIT/1ecPMDaX7AJBl5KS6TuGzX5Sozk5RtXhKY3M92LbWgWc9cKeUa+/VmAMcAsPqaA/A0Ed2ulFp0Ha/jCAywVz5znmFTZDbyssXVuAhNHptoe+4meRwA2kdmkplcNiuheU0uUjLnJYkFRdBtW0lctvMEUY1BYZqN5jo/8gpKXDazUPq4ALcv1iyPF2GZDnlb/bFdk4Sp4s11rm1bAaXUswAO8DIRnQNwm9qLrZAyjMJlFpiwOe5tRGVOzdYgWWnNc5tfUq5UXDll9kzeX64DtvfX9FMXNiIzu041AztFXuY5g5Q/6DX6EZft/ObUbBX0iteyhT4A9qBQ0/Qz1RQP5uFFXLbrC/I+2K7VdS/aDR1HYMD2UV0A71Yjpaq7yHj5vswvKfuwXMOq2ZSZrKjS9yV9YEwuxWJRN0qY28ipeW0mbP/J1koXavBteKLeCt4sgvXr8lOLme5FXOaUY7fMemNT9V5lN585z0tHvGwYksQlj2EjZNP0lvfLvEaXUttpAlNKHQ2yXccSWK1yWcJmVppkZBIXE5H8usp1LjNT+sS4gkiHe1dXlyYxXjZNSrOy+plWTIS2bc2KKF/8esmskcodVGXZ4BUgXIv5aSMsnprP1zbsnqnUveqDLLvt+UpFLglMEhdQnf6Jr6uR+i8/sl77tBs6jsBsctn18Gp5kExGZsXj9ZLQSqWSJh5TnUlykpVXzssKKv1opgKz+Tp4f79rdr2otq+xLJdZVj80o2LXSmJBGz+8juny/5iq3NaiaJqMQdWWK6DURlK2AW3lvubP73rldcp6JkmM66ENzVLqzUbHERhgj30x18up33Fcigy4SmDs8zBDJKSZaPrLuIK4VBkfD0DVIA22n/xPVjg/57/rem2VmOfNkAzed6e/wrZrqDf0gWESl6myAGwjK9Mpb/4HeHegtpmIQPUgtbZnDmCb8qpHcdnugXz+5nq/+9ZO6EgCA7a3ENXqH7L9b8pp+Z9pQsp5VmeS3GwmpqzkfC7pDLaZFeZXGaiu+KZiM6/XRWReJkMtJNYMYvMiWwnXyxvkZbYRl2nuuUZltykt88NkM9ldCTRtYzyaBGeSXb3E5fds5P9e9SEksCbD9hWxKQv5UMwH7qoAtn14amuN5HVSrbGJKbeRasfWACCdqFx5pTKTA5ZKnxr3DWUy9SOyoCabzaS0maCNkFiQZ+KnurzUp4u4ZG8K2WFa/udKYWMeV5bDFu4gQ2dcxGXze7qu3XWvzDK5zGS//2wICazJqPWG2kjJfJG91IpUO1zxTYXFy7LlEaju42aaKgzTb0ZEVaqMyUxOWYnxOWX3KvMLLq/Hi9RtcBGZSWJ87iAIorqCmou2Z2i6BWwfHX4eTGCcptzV4iyPK5U5L5uthHL0bGB7A41NoZk/874EVZjmvNfUi9xsx2wndCSBuW6wS33Z7H3A/RLa5k01J18Kk7ik/0KqpWg0qpvd+Vj8EtmODVSnwZYVn48pz8nL0qfC5TZfdPM+eN0XhklkLjXmeka2c9nWuVSHi8AkaiEuSWDmetNctJVXPheTuOQzMInMj7i8rs92rbZrb2RqQ+jEbxK87PGgJGYjNS8yk8teZGYjMunsl1leJZHxSwRc9YmZX0ZpOprl5/WuiH9guzPY9qLYFGktRMZltd0/E0FNRZ76qRHzJTT9jmaAqR9xmYpLloM/IkC1qpLEZZKYrUW5HUjLtc48R6jAmgjZQghsDw+oh8RsCg3wVgwuc8XmvJfqSBIZv0QmIcqvv6nIXGWX55AqTW5jNse7UA+Rue6XH2olLobtpZStgnzv5TOQgadm0kDznstrMlWwbEiR91sSmun78iMuLwQhnkbIC/BWWSGBtQDmDTcJzVQDrhfEj8TkPuY6L9ViHluqMKmUpBrgbfmlki+lzcSU18vH48BY+fLwC8fncMUVeV23i/hd+cm84EVatqkNpokIbCcum4NeKjG5r0nEfG55j8zQByYscz1PJXH5OeRtJFELKTWTsFzlaDd0JIG5iCYIoXmZkDz1IzVT7Zng/WSF5W3lf6YPi184GYJhtoTxeeVUrucXRRIlm7A2IrO1gPmpTq/76RoOzTyOeUzb1AWTuKSz3SQuc/zQII55k7hsHwJzan4wbKEQ5r2V8y6SMYmpEbUVhLBczy0ksCYjyAtiIzQvVWYjGnOdeS7T1OCyyW1NZ7o8rukvIyLEYjGrYjDjlkxzh/09kozYTJUhGLYXL8gLZ7vH8j74Ebtt/1pJC6ju9iXVlJmS2TQhTbVmlkn+bGagrUXRK5bL/CCYU7MOB1VVQdVWM03CkMCaiFpvJlcYF6HZ1FgQYjPPYZbLtq3tnHIbsw+mzflsUxImkXFZmbDkVCo00+Qxm/jNn+2+2u5HkA+M1zo+hvlz5ZG3TU3ikvdIPltJOiYp2dbxlJ+XjbBc5GXeGz+1VYup6CIs27Oo5R1ykX47oGUERkSfAXAfgCWl1Kss/xOATwC4F0AawLuVUk83et4gD8ZGaH5k5kViLlJznVeSnTRXgKsEZjN7+OXkUcldI9bICmcSmYu4zBY0aY76+XD8rjsobC+kvBc873LImyTvMhNN4pLmIQDEYrFtpO5SqebUdU+8SMSPpGolLZeyawT7UYF9FsDvAvhjx/8/CuCVld8dAH6vMvWE+RWtBybZ8AsuIVvpvNSZbWrOS5jrpTOdTUnAbSbJsAyzIcAVBiBfXqnEAHcYQFdX17aki17XXC+JeZlIJhmxf9BlMgLYprjM+26qTElU0iy0EZjNvHapUvP6bMu1EBjP+xFWo2rLhX1HYEqpx4joqMcmbwPwx6p8Zx4nolEimlFKLbSgLFXLNnNP/gdsD9WohcDM4wfx+ZgmBhOaVFZmi5rN8c8vN5fBVCKyYYOnrD54X56Xjmqp6Jh0g9wTE+bLaq6TilI64nkqW2yZxHgf+dzk9Zpmoqk65XXaHPau8AcveJGKF3HJ+yDXmcQV5PjNxL4jsAA4BOCiWJ6vrNtGYET0IIAHK4u5D33oQ99tffGahkkAnmlx2widVFags8rbSWUFgOvE/FdQLn8Q7Og1doQTX5WHZXoYAIjopPIY9aTd0Enl7aSyAp1V3k4qK1AuL88rpe7ZzbJ4wTsPb2txCcBhsTxXWRciRIgQgbCbBPYIgJ+hMl4HYKMV/q8QIULsXbQyjOJzAN6I8ki98wD+A4AYACilfh/Al1AOoTiNchjFzwY8dKBRftsInVTeTior0Fnl7aSyAh1SXmrX1oUQIUKE8MNumpAhQoQI0RBCAgsRIkTHom0JjIjuIaIXiOg0Ef2q5f8eIvqLyv9P+ATNthQByvoBIjpFRN8hov+PiF6xG+UU5fEsr9jux4lIEdGuNf8HKSsR/cvK/X2OiP5sp8tolMWvLhwhoq8R0TOV+nDvbpSzUpbPENESEVnjKisNbL9TuZbvENGtO11GX8huF+3yAxAFcAbAcQDdAP4ZwI3GNv8KwO9X5u8H8BdtXNY3AeivzP/CbpU1aHkr2w0BeAzA4wBua9eyotwV7RkAY5XlA+18b1F2jv9CZf5GAOd2sbw/COBWAN91/H8vgC8DIACvA/DEbpXV9WtXBXY7gNNKqbNKqTyAP0e565HE2wD8UWX+rwG8hertkNcYfMuqlPqaUipdWXwc5Zi33UKQewsAvwngowCyO1k4A0HK+l4An1RKrQOAUmpph8soEaS8CsBwZX4EwOUdLF91QZR6DMCaxya6u59S6nEAo0Q0szOlC4Z2JTBXNyPrNkqpLQAbACZ2pHSOclRgK6vEAyh/1XYLvuWtmAqHlVJf3MmCWRDk3l4L4Foi+iciepyIdjNqPEh5PwTgpyuhRV8C8G92pmh1oda6vePoiK5EewVE9NMAbgPwht0uiwtEFAHwcQDv3uWiBEUXymbkG1FWto8R0c1KqfhuFsoDPwngs0qp/0xEdwL4EyJ6lVKqPRNutTnaVYEF6WaktyGiLpTl+OqOlM5RjgqsXaKI6IcA/DqAtyqlcjtUNhv8yjsE4FUAHiWicyj7Ph7ZJUd+kHs7D+ARpVRBKfUygBdRJrTdQJDyPgDgLwFAKfXfAfQieEfpnUb7d/fbbSecw3nYBeAsgGO46gy9ydjmX6Paif+XbVzWV6Ps3H1lJ9xbY/tHsXtO/CD39h4Af1SZn0TZ5Jlo4/J+GeXknQBwA8o+MNrF+nAUbif+v0C1E/9bu1VOZ/l3uwAeN/ZelL+mZwD8emXdh1FWMED5y/VXKHdF+haA421c1n8AcAXAtyu/R9r53hrb7hqBBby3hLLJewrAswDub+d7i3LL4z9VyO3bAH5kF8v6OZTTVxVQVrIPAHgfgPeJe/vJyrU8u5v1wPULuxKFCBGiY9GuPrAQIUKE8EVIYCFChOhYhAQWIkSIjkVIYCFChOhYhAQWIkSIjkVIYHsIRHSYiF4movHK8lhl+WiLzvc+IvqZyvy7iWhW/PdpIrqxSed5OxE9VJn/LBG9o87jTBHR3zWjTCHaAyGB7SEopS6iPEDwRyqrPgLgYaXUuRad7/eVUjxw8bsBzIr/fk4pdapJp/oVAJ9q9CBKqWUAC0R0d+NFCtEOCAls7+G3AbyOiH4RwA8A+C1zAyI6SkTfI6I/JaLnieiviai/8t9bKrmqnq3ki+qprP+IyGn2W5V1HyKiD1YU0W0A/pSIvk1EfUT0KHc/IqKfrBzvu0T0UVGOJBH9RyL650pH7GlLWa8FkFNKbRtvkIh+s6LIokR0joj+j8r5TxLRrUT0FSI6Q0TvE7v9DYCfqvfmhmgvhAS2x6CUKgD4dygT2S9Wlm24DsCnlFI3ANgE8K+IqBfAZwH8hFLqZpS7xvwCEU0A+B9Q7hbzfQD+d+Ocfw3gJICfUkrdopTK8H8Vs/KjAN4M4BYAryWit1f+HgDwuFLq+1HOPfZeSznvBvC0uZKIPgZgCsDPKqWKldUXlFK3APjHynW8A+UuML8hdj0J4PWOexKiwxAS2N7Ej6LcReRVHttcVEr9U2X+v6Ks1q4D8LJS6sXK+j9COendBsp5wf6QiP5HlEeRCorXAnhUKbWsymmP/rRyTADIA/hCZf4plPvlmZgBsGys+98AjCil3qequ5I8Upk+i3LyvUTFbMwR0WjlvyUIUzdEZyMksD0GIroFwA+jrDx+ySMBndmHzNmnrEI8t6OcOPI+AM1yhBcEARVhT++UQbnfq8STAF7DjRUCnOWjJOZ5mY/dWzlmiD2AkMD2ECoZaX8PZdPxAoCPweIDq+BIJR8VALwTwDcAvADgKBGdqKx/F4CvE9EgyornSwB+CcD3W46XQDkVj4lvAXgDEU0SURTlfFhfr+Gyngdwwlj3dyg3UHyRiGzn9MK1AKw54EN0HkIC21t4L8p+oL+vLH8KwA1EZEug+AKAf01EzwMYA/B7SqksygMM/xURPYuycvl9lInpC0T0HZSJ7gOW430WwO+zE59XqvJo678K4GsoZ2B4Sin1+Rqu6TEArzbThSul/grAH6Ccq6zPuqcdbwKw25lmQzQJYTaKfYhKXNgXlFJePrK2ARF9AsDfKqX+oQnHegzA21Qlh36IzkaowEJ0Av4TgP5GD0JEUwA+HpLX3kGowEKECNGxCBVYiBAhOhYhgYUIEaJjERJYiBAhOhYhgYUIEaJjERJYiBAhOhb/P97o1qs4mfnsAAAAAElFTkSuQmCC\n" }, "metadata": { "needs_background": "light" } } ], "source": [ "#NBVAL_IGNORE_OUTPUT\n", "from examples.seismic import plot_shotrecord\n", "\n", "plot_shotrecord(rec.data, model, t0, tn)" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "assert np.isclose(np.linalg.norm(rec.data), 370, rtol=1)" ] } ], "metadata": { "hide_input": false, "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.8.9-final" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autoclose": false, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 1, "hotkeys": { "equation": "Ctrl-E", "itemize": "Ctrl-I" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": false, "user_envs_cfg": false }, "widgets": { "state": {}, "version": "1.1.2" } }, "nbformat": 4, "nbformat_minor": 4 }