access_and_analyse.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Access and analyze seismograms\n",
    "\n",
    "_This notebook was created by Léonard Seydoux (seydoux@ipgp.fr) in 2022 for the course \"Scientific Computing for Geophysical Problems\" at the [institut de physique du globe de Paris](http://www.ipgp.fr)._ \n",
    "\n",
    "<img width=200px src=\"https://upload.wikimedia.org/wikipedia/fr/thumb/1/16/Institut_de_physique_du_globe_de_paris_logo.svg/1200px-Institut_de_physique_du_globe_de_paris_logo.svg.png\" />\n",
    "\n",
    "## Table of contents\n",
    "\n",
    "1. [Introduction](#1.-Introduction)\n",
    "    1. [Goals](#1.1.-Goals)\n",
    "    2. [Requirements](#1.2.-Requirements)\n",
    "    3. [Searching events of interest](#1.3.-Searching-events-of-interest)\n",
    "2. [Seismic waveforms](#2.-Seismic-waveforms)\n",
    "    1. [Downloading seismic waveforms](#2.1.-Downloading-seismic-waveforms)\n",
    "    2. [Inspection of the metadata](#2.2.-Inspection-of-the-metadata)\n",
    "    3. [Inspection of the waveform](#2.3.-Inspection-of-the-waveform)\n",
    "    4. [Saving and reading waveforms locally](#2.4.-Saving-and-reading-waveforms-locally)\n",
    "    5. [Instrument response correction](#2.5.-Instrument-response-correction)\n",
    "3. [Spectral analysis](#3.-Spectral-analysis)\n",
    "    1. [Fourier transform](#3.1.-Fourier-transform)\n",
    "    2. [Filtering](#3.2.-Filtering)\n",
    "    3. [Spectrogram](#3.3.-Spectrogram)\n",
    "4. [Analysis with a selected event](#4.-Analysis-with-a-selected-event)\n",
    "\n",
    "\n",
    "## 1. Introduction\n",
    "\n",
    "### 1.1. Goals\n",
    "\n",
    "This notebook shows how to download, plot, and analyze seismic waveforms in time and frequency. It goes from identifying events of interest by exploring online resources with reported events to extracting the seismic waveforms from suited data centers and analyzing their properties. \n",
    "\n",
    "### 1.2. Requirements\n",
    "\n",
    "This notebook relies on Python libraries listed below. If you are running this notebook on the virtual machine provided for the course, you should have all the required libraries installed. If you are running this notebook on your own machine, you will need to install the following libraries:\n",
    "\n",
    "- [ObsPy](https://github.com/obspy/obspy/wiki), an open-source project that provides a Python framework for processing seismological data. It provides parsers for standard file formats, clients to access data centers, and signal processing routines that allow the manipulation of seismological time series.\n",
    "- [NumPy](https://numpy.org), an open-source project aiming to enable numerical computing with Python.\n",
    "- [SciPy](https://numpy.org), fundamental algorithms for scientific computing in Python.\n",
    "- [Matplotlib](https://matplotlib.org), a comprehensive library for creating static, animated, and interactive visualizations in Python.\n",
    "\n",
    "If unavailable on your machine, you can use the [Anaconda](https://www.anaconda.com) package manager to install these libraries. Please ensure you can run the following cell without any error before proceeding further."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import obspy \n",
    "from obspy.clients.fdsn import Client\n",
    "from obspy.geodetics import locations2degrees\n",
    "from scipy.signal import spectrogram"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3. Searching datasets of interest\n",
    "\n",
    "#### Motivations\n",
    "\n",
    "Looking for seismic data is generally motivated by a scientific question or by curiosity to look at some recent events. Here we deliver the sources where this information is accessible. You may also find crucial information from other regional and national earthquake information services (volcano observatories, glacial area observatories).\n",
    "\n",
    "- The [United States Geological Survey (USGS) Advanced National Seismic System Comprehensive Catalog of Earthquakes Events and Products](https://earthquake.usgs.gov/earthquakes/search) provides one of the most comprehensive global earthquake catalogs. It delivers earthquake source parameters (e.g. hypocenters, magnitudes, phase picks, and amplitudes) and other products (e.g. moment tensor solutions, macroseismic information, tectonic summaries, maps) produced by contributing seismic networks. The webpage https://earthquake.usgs.gov/data/comcat/ provides a list of contributing seismic networks.\n",
    "- The [ISC-GEM Global Instrumental Earthquake Catalogue](http://www.isc.ac.uk/iscbulletin/search/catalogue/) (1904-2018) is the result of a special effort to adapt and substantially extend and improve currently existing bulletin data for large earthquakes (magnitude 5.5 and above, plus continental events down to magnitude 5.0) to serve the requirements of the specific user group who assess and model seismic hazard and risk.\n",
    "- The [European-Mediterranean Seismological Centre (EMSC), France](https://www.emsc-csem.org) was founded in 1975, following a recommendation from the European Seismological Commission (ESC). The ESC is a regional commission of the International Association of Seismology and Physics of the Earth's Interior (IASPEI), itself a specialized association of the International Union of Geodesy and Geophysics (IUGG).\n",
    "\n",
    "#### Sensors\n",
    "\n",
    "Seismic stations or networks are also listed in several websites. The [International Federation of Digital Seismograph Networks (FDSN)](https://www.fdsn.org/about/) goals related to station siting and instrumentation are to provide stations with good geographic distribution, recording data with 24 bits of resolution in continuous time series with at least a 20 sample per second sampling rate. You can find a [list of seismic networks with links to station maps and descriptions](https://www.fdsn.org/networks/) online to fetch seismic data on specific regions and temporal coverage. Also, the [Incorporated Research Institutions for Seismology (IRIS)](https://ds.iris.edu/gmap/) is a consortium of universities dedicated to the operation of science facilities for the acquisition, management, and distribution of seismological data. IRIS delivers seismic waveforms and metadata based on the FDSN inventory, which facilitates the [search of the temporary or permanent seismic network](https://ds.iris.edu/gmap/). Here are some famous ones:\n",
    "\n",
    "- [IU, the Global Seismograph Network (GSN) operated by IRIS and USGS](https://www.fdsn.org/networks/detail/IU/), described by a global distribution of high-gain, co-located broadband velocity instruments and accelerometers.\n",
    "- [II, the Global Seismograph Network (GNS) operated by IRIS and IDA](https://www.fdsn.org/networks/detail/II/).\n",
    "- [G, GEOSCOPE](https://www.fdsn.org/networks/detail/G/) is a global seismological network created in 1982 to respond to the challenge of instrumenting isolated sites all around the world. More information is available at the [Institut de Physique du Globe de Paris website](http://geoscope.ipgp.fr/index.php/en/stations/station-map).\n",
    "- [FR: RESIF](https://www.fdsn.org/networks/detail/FR/) and other broadband and permanent accelerometric networks in metropolitan France allow monitoring of the seismicity in metropolitan France and the development of seismic risk and hazards maps.\n",
    "- [PF: Piton de la Fournaise Volcano Observatory Network (Reunion Island)](https://www.fdsn.org/networks/detail/PF/) allows monitoring of the volcano-seismic activity of the Piton de la Fournaise volcano. \n",
    "- [GL: Guadeloupe Seismic and Volcano Observatory Network](https://www.fdsn.org/networks/detail/GL/) allows monitoring of the volcano-seismic activity of La Soufrière volcano. \n",
    "- [HV: Hawaiian Volcano Observatory Network](https://www.fdsn.org/networks/detail/HV/) allows monitoring of the volcano-seismic activity in Hawaii. \n",
    "- [TA: USArray Transportable Array](https://www.fdsn.org/networks/detail/TA/) is a transportable array that was deployed between 2004 and 2019 in the US.\n",
    "\n",
    "#### Data centers\n",
    "\n",
    "Finally, a comprehensive list of datacentersthat support the FDSN webservices can be found at https://www.fdsn.org/webservices/datacenters/. Among others, you can find here a shortlist of international and French data centers:\n",
    "- [IRIS Data Management Center](https://ds.iris.edu/ds/nodes/dmc/)\n",
    "- [EIDA: European Integrated Data Archive](http://www.orfeus-eu.org/data/eida/)\n",
    "- [RESIF: French seismological and geodetic network](http://seismology.resif.fr)\n",
    "- [IPGP: data center of the Institut de Physique du Globe de Paris](http://datacenter.ipgp.fr/data.php)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Seismic waveforms\n",
    "\n",
    "### 2.1. Downloading seismic waveforms\n",
    "\n",
    "Several tools allow downloading seismic waveforms provided by the datacenters via the [FDSN Web Services](https://www.fdsn.org/webservices/) which is [natively included in the obspy Python library](https://docs.obspy.org/packages/obspy.clients.fdsn.html). Here we show an example of how to download seismic data from the IRIS datacenter with the `Client` interface of `obspy`. Once the interface is connected to a client of interest, we can request a seismic waveform by specifying several metadata such as the network code, the station code, the location of the sensor in the station, the channel name, the start and end datetimes, and other information. Please check the [documentation of the `get_waveform` method](https://docs.obspy.org/packages/autogen/obspy.clients.fdsn.client.Client.get_waveforms.html).\n",
    "\n",
    "Seismic records are organized in a specific type of object in Python called a `trace`. Since several traces can be downloaded at once, several traces are gathered into a `stream` object. By default, the `get_waveform` method returns a `stream` event if only one trace was downloaded. You can print a succinct amount of information about the stream with the `print` function.\n",
    "\n",
    "> __Question 1.__ What kind of seismic network is the IU network (global, local)? Where is the seismic station CTAO located? From the date and time used in the above example waveform query, which earthquake was recorded by the seismic sensor CTAO? What is the approximate epicentral distance of the earthquake to the seismic station? Please use the different websites provided in the previous section together with the obspy documentation to answer that question. You can also use the obspy method [`location2degrees`](https://docs.obspy.org/packages/autogen/obspy.geodetics.base.locations2degrees.html) to answer the last question."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2 Trace(s) in Stream:\n",
      "IU.CTAO.00.BHZ | 2021-03-04T19:28:33.019538Z - 2021-03-04T21:29:59.994538Z | 40.0 Hz, 291480 samples\n",
      "IU.CTAO.10.BHZ | 2021-03-04T19:28:33.019538Z - 2021-03-04T21:29:59.994538Z | 40.0 Hz, 291480 samples\n"
     ]
    }
   ],
   "source": [
    "# Connect to IRIS\n",
    "client = Client(\"IRIS\")\n",
    "\n",
    "# Get the data\n",
    "stream = client.get_waveforms(\n",
    "    network=\"IU\",\n",
    "    station=\"CTAO\",\n",
    "    location=\"*\",\n",
    "    channel=\"BHZ\",\n",
    "    starttime=obspy.UTCDateTime(\"2021-03-04T19:28:33\"),\n",
    "    endtime=obspy.UTCDateTime(\"2021-03-04T21:30:00\"),\n",
    "    attach_response=True,\n",
    ")\n",
    "\n",
    "# Print the content of the stream\n",
    "print(stream)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here to answer question 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2. Inspection of the metadata\n",
    "\n",
    "A `stream` can be seen as a list (actually, an iterable) of `trace` objects, plus methods. More precisely, almost all of the methods available in a `trace` object is available in a `stream` object, to be applied on all traces individually. In the following examples, we will narrow the focus on the first trace of the downloaded stream, and take a look at its metadata. A trace's metadata are stored in a dictionnary named `stats`. The entries of this dictionnary may depend on the format of the seismic trace. The most common seismic formats are `SAC` and `MSEED` (visible in the `_format` entry of the `stats` dictionnary).\n",
    "\n",
    "> __Question 2.__ Take time to discover the `stats` dictionnary of a selected `trace` and understand the different entries. More information about this dictionnary can be found in the [obspy documentation](https://docs.obspy.org/packages/autogen/obspy.core.trace.Stats.html). In particular, what are the `npts`, `delta` and `sampling_rate` entries? Please define everything with respect to the digital signal $u[n] = u(n\\Delta t)$ and other quantities seen during the class. In which format was the trace downloaded? Considering the online [documentation about the data formats on IRIS](http://ds.iris.edu/ds/nodes/dmc/data/formats/), what does `number_of_records` mean in that case? What information is presented in the `response` entry? What is the signification of a sensor sensitivity? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3. Inspection of the waveform\n",
    "\n",
    "In the following example, we define a time vector in number of samples $t = 0\\ldots N$. Modify the cell below to turn this time vectors into a relative time (in seconds) from the begining of the record based on the entries of the `stats` dictionary, and update the x-axis name accordingly. \n",
    "\n",
    "> __Question 3.__ How would you proceed to use an absolute time array with date ticks? What is the type and unit of the raw `data`? Please name the axis accordingly, and explain your answer. What would you need to turn the waveform into physical units? Modify the code below to name the axes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkcAAAHFCAYAAAD40125AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABiiElEQVR4nO3dd3hTVQMG8Ddt0z2gdFE62aNsEIpgGVL2EEVkCYJ+LhDEicoSBFRwKyoyRGUoS1RWkT1lU/Yqs2UU6IY2bc73x6XpTZO0SZs0Sfv+nqcPyc3JzcnppXlzzrnnKoQQAkREREQEAHCwdgWIiIiIbAnDEREREZEMwxERERGRDMMRERERkQzDEREREZEMwxERERGRDMMRERERkQzDEREREZEMwxERERGRDMMRUTmxcOFCKBQKHDhwAAAwfPhweHp6Gizv6emJ4cOHG7XvtLQ0fPTRR2jRogW8vb3h4uKCiIgIjBgxAocOHQIAKBQKo362bt2q2e9XX30FhUKBqKioIl//1KlTGD58OMLCwuDs7Aw/Pz90794d69atM6r+J06cwCuvvILo6Gh4eHjo1KOwpUuXokmTJnB1dUVwcDDGjh2LjIyMYl/n0qVLOu/X29sbjRs3xhdffIG8vDyt8u3btzf43pOTk6FQKDB58mTNtoiIiGLb19jfKREZ5mTtChCRbbtw4QJiY2Nx69YtvPTSS5gyZQo8PT1x6dIl/P7772jevDlSUlKwZ88eredNnToVW7ZswebNm7W2169fX3N7/vz5AKTwsm/fPrRq1Urn9VeuXIlBgwahevXqmDBhAurUqYObN29iwYIF6N69O9566y188sknRb6HAwcOYPXq1WjatCk6deqEv/76y2DZ3377DUOGDMHzzz+Pzz//HGfPnsU777yDkydPYuPGjcW2FwCMHj0agwYNAgCkpKRgzZo1eP3113H16lXMnj3bqH3os2rVKmRnZ+t97LXXXsOBAwfQt2/fEu+fiB4SRFQuLFiwQAAQ+/fvF0IIMWzYMOHh4WGwvIeHhxg2bFiR+8zNzRUNGzYU3t7eIj4+Xm+ZtWvXiszMTJ3txb3+/v37BQDRo0cPAUC88MILOmXOnz8v3N3dRYsWLURGRobO4y+99JIAIJYsWVLk+8jLy9Pc/uOPPwQAsWXLFp1yubm5omrVqiI2NlZr+2+//SYAiLVr1xb5OgkJCQKA+PTTT3Uea9eunahatarWtpiYGNGgQQO9+7p9+7YAICZNmlTkawohxOzZswUAMX78+GLLElHxOKxGRAatXr0a8fHxGD9+vMHhn27dusHd3d3kfc+bNw8AMHPmTLRp0wZLly5FVlaWVpnPP/8cWVlZ+Prrr+Hh4aGzj9mzZ6NSpUr46KOPinwtBwfj/tTt3bsXSUlJeO6557S29+/fH56enli1apVR+9HHx8cHSqWyxM83ZMuWLXjnnXcQGxuLadOmmX3/RBURwxERGZQ/jGTuoZr79+9jyZIlaNmyJaKiojBixAikp6fjjz/+0CoXFxeHwMBAtG7dWu9+3N3dERsbi+PHj+PGjRulrtfx48cBAI0aNdLarlQqUbduXc3jxVGr1cjNzUVubi7u3LmD+fPnY/369Rg6dKje8vll5T+F5yfpc+XKFQwYMAAhISFYsmSJ0SGQiIrG/0lEZNCVK1cAAJGRkWbd7/Lly5GamoqRI0cCAAYMGABPT09Nb5L89Yt77fzH8+taGnfu3AEA+Pr66jzm6+urebw477zzDpRKJZRKJfz8/DBy5EgMGTIEU6ZM0Sl74sQJTVn5T1BQUJGv8eDBA/Tr1w8ZGRlYuXKl3joTUclUyHC0fft29OrVC8HBwVAoFFi9erXJ+xBCYNasWahduzZcXFwQGhqK6dOnm7+yROXQvHnz4ObmhmeeeQaAdOZc//79sWPHDpw7d86kfQkhAEhny5mLoX0Z+xpjxozB/v37sX//fmzZsgXTp0/H77//joEDB+qUrVGjhqas/GfTpk1FvsZLL72EgwcP4ocffkDTpk2NqhcRGadCnq2WmZmJxo0b47nnnsOTTz5Zon2MGTMGGzduxKxZs9CwYUOkpqYiOTnZzDUlKjknJ6cih2Zyc3OLnQMTFhYGAEhISEDdunXNUq/z589j+/btePLJJyGEQEpKCgDgqaeewoIFCzB//nzMmDFD8/oJCQlF7u/SpUsAgNDQ0FLXrUqVKgCkHqTAwECtx+7evWt070xISAhatGihud++fXsoFAqMHz8eGzZsQJcuXTSPubq6apXNV9Tfk6+//ho///wzRo0aZXCojohKrkL2HHXr1g3Tpk1Dv3799D6ek5ODt99+G9WqVYOHhwdatWqltSbKqVOnMGfOHPz555/o3bs3IiMj0aRJEzz++ONl9A6IihcYGIgHDx7g7t27Oo/duXMH2dnZOgGgsPwP8ZL0rhoyf/58CCGwfPlyVK5cWfPTo0cPAMDPP/+sCXWdO3fGzZs3sXfvXr37ysrKQlxcHKKiooodhjJGw4YNAQDx8fFa23Nzc3H69Oli12MqSv48pqNHj5a8ggB27NiBN954A23btsVnn31Wqn0RkX4VMhwV57nnnsOuXbuwdOlSHDt2DP3790fXrl013f1//fUXqlevjr///huRkZGIiIjA888/r/dDiMha8sP6smXLdB77/ffftcoY0qdPHzRs2BAzZswwOBl5w4YNOmeZGZKXl4eff/4ZNWrUwJYtW3R+3njjDSQlJWkWd3z99dfh5uaG0aNHIzMzU2d/b775Ju7du4cPPvjAqNcvTqtWrVC1alUsXLhQa/vy5cuRkZFh8AuVMY4cOQIACAgIKPE+rl+/jv79+8PPzw9//PGHRc5+I6IKOqxWlAsXLmDJkiW4du0agoODAUh/gNevX48FCxZg+vTpuHjxIi5fvow//vgDixYtQl5eHl5//XU89dRTOgveEVlLhw4d0Lt3b4wZMwaXLl1CTEwMhBDYvn07Pv/8c/Tu3Rvt27fXlL98+TJq1KiBYcOGaSZGOzo6YtWqVYiNjUV0dDRefvlldOjQAR4eHrh8+TKWL1+Ov/76C/fu3TOqTuvWrUNiYiI+/vhjrdfOFxUVhW+++Qbz5s1Dz549UaNGDfzyyy8YPHgwWrZsiXHjxmkWgZw/fz7WrVuHN998EwMGDNDso1OnTti2bRtyc3M127KysrB27VoA0PRCbdu2DcnJyfDw8EC3bt007/eTTz7B0KFD8eKLL2LgwIE4d+4c3n77bXTu3Bldu3bV7HPbtm3o1KkTJk6ciIkTJ2q9jytXrmheJzMzE3v27MGMGTMQHh5e4oCVk5ODfv364ebNm5g9ezYuXbqkGVKU8/b21lpok4hKwKqrLNkAAGLVqlWa+7///rsAIDw8PLR+nJycxNNPPy2EEOKFF14QAMSZM2c0zzt48KAAIE6fPl3Wb4FICKG7CKQQQuTk5Ijp06eLBg0aCBcXF+Hi4iIaNGggpk+fLnJycrSen7+Aob6FIVNSUsTUqVNFs2bNhKenp1AqlSIsLEwMGTJE7Nq1S2999C0C2bdvX+Hs7Cxu3bpl8H0888wzwsnJSdy4cUOz7cSJE2LYsGEiJCREKJVK4evrK7p27Sr++ecfnefHxMSIwn/a8t+bvp/w8HCdfSxevFg0atRIODs7i6CgIPHaa6+J9PR0rTJbtmzRWaRR3+u4urqK2rVri7Fjx4qkpCSduhq7CGRR70H+ExMTY7Bticg4CiEenupRQSkUCqxatUqzjsuyZcswePBgnDhxAo6OjlplPT09ERQUhEmTJmH69OlQqVSax+7fvw93d3ds3LgRnTt3Lsu3QERERGbEYbVCmjZtiry8PNy6dQvt2rXTW+bRRx9Fbm4uLly4gBo1agAAzp49CwAIDw8vs7oSERGR+VXInqOMjAycP38egBSGPvvsM3To0AG+vr4ICwvDkCFDsGvXLsyePRtNmzZFcnIyNm/ejIYNG6J79+5Qq9Vo2bIlPD098cUXX0CtVuPVV1+Ft7e30RemJCIiIttUIcPR1q1b0aFDB53tw4YNw8KFC6FSqTBt2jQsWrQI169fR5UqVRAdHY0pU6ZoTvVNTEzE6NGjsXHjRs2EztmzZ3OVWiIiIjtXIcMRERERkSFc54iIiIhIhuGIiIiISKZCna2mVquRmJgILy8vs16kkoiIiCxHCIH09HQEBwfDwcHy/ToVKhwlJiaa5eKUREREVPauXr2KkJAQi79OhQpHXl5eAKTG9fb2Ntt+VSoVNm7ciNjYWF7ryARsN9OxzUzHNisZtpvp2GamM7bN0tLSEBoaqvkct7QKFY7yh9K8vb3NHo7c3d3h7e3N/xAmYLuZjm1mOrZZybDdTMc2M52pbVZWU2I4IZuIiIhIhuGIiIiISIbhiIiIiEiG4YiIiIhIhuGIiIiISIbhiIiIiEiG4YiIiIhIhuGIiIiISIbhiIiIiEiG4YiIiIhIhuGIiIiISIbhiIiIiEiG4YisKicPUKuFtatBRESkwXBEVpOSpcLb/zli8Pz91q4KERGRBsMRWc3Ws7choMCByynWrgoREZEGwxGVqZSsHNzLzAEAKBQKK9eGiIhIl5O1K0AVR55aoMmHcQCA01O7Wrk2RERE+rHniMrMA1We5vbt9Gwr1oSIiMgwhiMqM/JRNCEADqoREZEtYjiiMqMoFIc45YiIiGwRwxFZhQDXNiIiItvEcERWw44jIiKyRQxHZBWCHUdERGSjGI6ozHAojYiI7AHDEVmFABeBJCIi28RwRFYhOK5GREQ2iuGIykzhPMR+IyIiskUMR2QV7DciIiJbxXBEZYaBiIiI7AHDEVmFEFwhm4iIbBPDEZUZTsImIiJ7wHBEVsKgREREtonhiKxCGlbjuBoREdkeuwlHM2bMQMuWLeHl5YWAgAD07dsXZ86csXa1yATsKyIiIntgN+Fo27ZtePXVV7F3717ExcUhNzcXsbGxyMzMtHbVqATYaURERLbKydoVMNb69eu17i9YsAABAQE4ePAgHnvsMSvVikzB+dhERGQP7CYcFZaamgoA8PX1NVgmOzsb2dnZmvtpaWkAAJVKBZVKZba65O/LnPssj3Jl7aNS5SIvL1d2n21nDB5rpmOblQzbzXRsM9MZ22Zl3aYKYYfnVwsh0KdPH9y7dw87duwwWG7y5MmYMmWKzvbFixfD3d3dklUkPbJygfH7pTw+vnEubtxXYMFZRwDAl9G5RT2ViIgqsKysLAwaNAipqanw9va2+OvZZTh69dVX8c8//2Dnzp0ICQkxWE5fz1FoaCiSk5PN2rgqlQpxcXHo3LkzlEql2fZb3qTeV6HF9C0AgHWj2+DMjVSM/eMEAODc1FhrVs1u8FgzHdusZNhupmObmc7YNktLS4Ofn1+ZhSO7G1YbPXo01qxZg+3btxcZjADAxcUFLi4uOtuVSqVFDlxL7be8UMp6RZVKJzg6Osnus91MwWPNdGyzkmG7mY5tZrri2qys29NuwpEQAqNHj8aqVauwdetWREZGWrtKZCLBk/mJiMgO2E04evXVV7F48WL8+eef8PLywo0bNwAAPj4+cHNzs3LtiIiIqLywm3WO5syZg9TUVLRv3x5Vq1bV/CxbtszaVSMj2d/sNiIiqojspufIDueNU5EU/J0SEZFNspueI7J/8ii0cHeC1epBRERUFLvpOaLy5de9V1DVW/dMQiIiImtjzxFZTdp9LvxIRES2h+GIykzhOUY/7ODQGhER2R6GI7IJeWqBMzfSOUmbiIisjuGIykxRsWfin8fR5Yvt+HzTuTKrDxERkT4MR2QTftt3BQDw1b8MR0REZF0MR1RmOGJGRET2gOGIiIiISIbhiEplw4kb6PrFdpy9mW7tqhAREZkFwxGVyou/HMTpG+l4+deDxZYVRU7JJiIisg0MR2QWqVzQkYiIygmGIzILo9YnYscRERHZAYYjMos8nopGRETlBMMRmYVaXXw4YnwiIiJ7wHBEZsGOIyIiKi8YjoiIiIhkGI7IrLJycpFnYIiNvUtERGQPGI7ILASAm2kP0GLaJrz4S/FrHhEREdkqhiMym7iTN5GVk4dNp27qfZyLQBIRkT1gOCKzUABwdFBYuxpERESlxnBEZpGenYsJq48XWYZzjoiIyB4wHJHZ5Bqx1hEREZGtYzgiIiIikmE4ojLDfiUiIrIHDEdEREREMgxHVGYEZ2QTEZEdYDgiIiIikmE4ojLDjiMiIrIHDEdkkxbsSsC2s7etXQ0iIqqAnKxdAaLCPo87iy//PQcAuDSzh5VrQ0REFQ17jsjm5AcjIiIia2A4ojIhhMCoxYesXQ0iIqJiMRxRmThzMx1Hr6VauxpERETFYjiiMpHH664REZGdYDgiixi37AjSH6isXQ0iIiKTMRyRRaw8fB3fbrmgua+Awoq1ISIiMh7DEVlMUup9zW0FsxEREdkJhiOymFzOMyIiIjvEcESWI8tG7DkiIiJ7wXBEREREJMNwRERERCTDcERlgmerERGRvWA4IsthHiIiIjvEcERlghOyiYjIXjAckcVkq/I0t5mNiIjIXjAckcVsOnULi/ddsXY1iIiITMJwRBb13qp4a1eBiIjIJAxHVCY454iIiOwFwxERERGRDMMRERERkQzDEZURjqsREZF9YDiiMiKKL0JERGQDGI6IiIiIZBiOiIiIiGQYjqhMCI6qERGRnWA4IiIiIpJhOKIywY4jIiKyFwxHRERERDIMR2RxV+5kcc4RERHZDYYjsrjHPt2CHedul+i5gqmKiIjKmF2Fo+3bt6NXr14IDg6GQqHA6tWrrV0lMtK0f06V6HlL9181c02IiIiKZlfhKDMzE40bN8Y333xj7apQGRm/Mt7aVSAiogrGydoVMEW3bt3QrVs3a1eDiIiIyjG7Ckemys7ORnZ2tuZ+WloaAEClUkGlUpntdfL3Zc59UgG2awEea6Zjm5UM2810bDPTGdtmZd2mCmGnM14VCgVWrVqFvn37GiwzefJkTJkyRWf74sWL4e7ubsHaVRxj9lg+X38ZnWvx1yAiItuVlZWFQYMGITU1Fd7e3hZ/vXIdjvT1HIWGhiI5OdmsjatSqRAXF4fOnTtDqVSabb/2oNaEjRZ/jXNTYy3+GvaiIh9rJcU2Kxm2m+nYZqYzts3S0tLg5+dXZuGoXA+rubi4wMXFRWe7Uqm0yIFrqf1WdGxTXTzWTMc2Kxm2m+nYZqYrrs3Kuj3t6mw1IiIiIkuzq56jjIwMnD9/XnM/ISEBR44cga+vL8LCwqxYMyIiIiov7CocHThwAB06dNDcHzduHABg2LBhWLhwoZVqRUREROWJXYWj9u3b83ISFdD2s7fxWG1/a1eDiIgqCM45Ipv37Pz/rF0FIiKqQBiOyC4cunLP2lUgIqIKguGI7MKP2y5auwpERFRBMByRXchVc64ZERGVDYYjsgs5eWr8sO0CTiSmWrsqRERUztnV2WpUcW0/exvbz94G1gGXZvawdnWIiKgcY88RERERkQzDEREREZEMwxERERGRjOnh6P59ICur4P7ly8AXXwAbN5qvVkRERERWYno46tMHWLRIup2SArRqBcyeLW2fM8e8tSMiIiIqY6aHo0OHgHbtpNvLlwOBgVLv0aJFwFdfmbl6RERERGXL9HCUlQV4eUm3N24E+vUDHByA1q2lkEQVBi8CTERE5ZHp4ahmTWD1auDqVWDDBiA2Vtp+6xbg7W3e2pFNYzYiIqLyyPRwNHEi8OabQESENN8oOlravnEj0LSpeWtHNo3ZiIiIyiPTV8h+6imgbVsgKQlo3Lhge6dO0hAbVRgcViMiovLI9J6jESMADw+pl8hB9vQGDYCPPzZj1cjWMRoREVF5ZHo4+vlnaa2jwu7fLzjFnyoENXuOiIioHDJ+WC0tTZqBKwSQng64uhY8lpcHrF0LBARYoIpkq5iNiIioPDI+HFWqBCgU0k/t2rqPKxTAlCnmqxkRERGRFRgfjrZskboKOnYEVqwAfH0LHnN2BsLDgeBgC1SRbBV7joiIqDwyPhzFxEj/JiQAYWFSTxFVaIJTsomIqBwyLhwdOwZERUlnp6WmAvHxhss2amSmqpGtUzMbERFROWRcOGrSBLhxQ5pw3aSJ1Gukb0xFoZAmZ1OFwHWOiIioPDIuHCUkAP7+BbeJwHWOiIiofDIuHIWH679NFRo7joiIqDwy/fIhAHD2LLB1q3SxWbVa+7GJE0tfK7ILHFYjIqLyyPRwNHcu8PLLgJ8fEBSkfdaaQsFwVIEwGxERUXlkejiaNg346CPgnXcsUB2yJ8xGRERUHpl+bbV794D+/S1QFbJlW87cwrL9V7S2cViNiIjKI9PDUf/+wMaNFqgK2bLnFuzHOyvice5mumabNaOREAL3MnOsWAMiIiqvTB9Wq1kTmDAB2LsXaNgQUCq1H3/tNTNVjWyFvIfodno2agV6AQDUVuw5endFPJYduIpfRj6CdrX8rVYPIiIqf0wPRz/+CHh6Atu2ST9yCgXDUTmUnVtwRqLSSdbZaMWuo2UHrgIAPo87y3BERERmZXo44iKQFc4DVcGq586OBeHIFmYcnb6RXnwhIiIiE5g+54jKndw8NbJzDV/2JSevoOdIvnKDLczHzsrJw4FLd822v8NX7mHG2lPIysk12z6JiMi+mN5zNGJE0Y/Pn1/CqpC1xH6+HckZ2dj/weNwcXLUeVy+zqf8YrPCJvqOgKn/nIKzowITetZHo5BKpdrXE9/tlm4ogPHd6pW+ckREZHdKdiq//OfWLWDzZmDlSiAlxfw1JItSqwUuJmci7UEuLtzK1FsmV5aO8h6mo7QHKtzPsY2LDB+9moL9l+6h//d7zLbP00kcriPbte/iHYxfGY/0ByprV4WoXDK952jVKt1tajXwyitA9epmqBKVpVxZV5Chs8/kPUdCCNzPyUOjyba3nIN84nhp5alto1eMSJ8BP+4FADg5KDC1b5SVa0NU/phnzpGDA/D668Dnn5tld1R25CFgz4U7+svIQlOeWuDcLdvvVbl2LwsT/zyOS8n6e8OKs/N8MgMS2byLyRnWrgJRuWS+CdkXLgC5nMRqb+RDZh+tPaW3TJ5W7xKggEJvOVsghMDS/66g7cdbsGjPZTzz8Bu2Icevp6L9p1uwLj5J57Ea763F+6viLVVVIiKyUaYPq40bp31fCCApCfjnH2DYMDNVi8pKbl7xvSN5hYbeFLabjRB38ibeXVkQaG6kPdBbbu72i1h/4gau3s3CrfRsvPzbIXSuH6hT7rd9V/DREw0tVl+i0th1/g7uZubA18PZ2lUhKldMD0eHD2vfd3AA/P2B2bOLP5ONbE6uEUNHhcORLTt707ghP329ZHEnb5q7OkQWN+HP4/h2UDNrV4OoXDE9HG3ZYoFqkLXIh9UMUReac2TLHBxsuFuLyAISbpdsXh0RGWZ6OKJyxdRhNSFsu/fIoZgxP1WeGqsPXzdpn0IIKGx5LJHKvb0X7yDAywWL9lzGtrO3tR6z5f+PRPaK4aiceKDKw7FrqWgWVglOjsbPszemJ0g+9JanFkYNxVnSn0cMhxvHIkKMEAJvLz+GVSaHI5RqntUDVR5+2XMZHer6o2aAl94y9zJzcDE5E83DK5f8hahcupSciWd+3ItAbxfcTMvWeZzZiMj8ePmQcmLU4sN4+oc9+PLfcyY9z9RhNbUQVh9aG7P0iEnl1Q/ru+XMLZODEVD6b+YLdl3CR2tP4fHPthssE/vFdjw5ZzdWHb4GUYafdpnZucjNM9/6UGR++fPo9AUjQDo+1WqBNC4ISWQ2DEflxKZT0mTi+TtNuzCwvBeohYFei8ITso0ZirMWfZc0mb8rAV//ew4jFh4o0T5LmwX/S9C/fpTc7XTpg+/1ZUfxmonhr6TSHqjQYNIGdPnCcGgj63M0Yh7dkHn70GjyRly5k1UGNSIq/0oXjh7oP02aLOtUUhpS7+v/lmjq57g86BiazFx4nSNr9xwVRV/Vpv1zCrPjzpZin6V7v8XNgyrsr6OJAIBvt5zH1yb2BJri4KV7AIALnNBr04o7fBQKYPfDBVxXHLpWBjUiKv9MD0dqNTB1KlCtGuDpCVy8KG2fMAGYN8/M1aPCDl6+i25f7kCTD/VfvsPUD3J5z5Gh4ZU8nTlHtjsMY4kRKWP2KYTA8eupyMwuWAg1NUuF49dTSzRfqe+3u/DphjOYHXcWP+24iF5f70Riyn2j6vHz7ks4ePlusWWN6ZEg6yvu+Lt+r+C4sOUvLkT2xPRwNG0asHAh8MkngLNs4bGGDYGffjJfzUivTaduAZD+YP5zTHdVZ1P/NubJgo6hidbyy4fcz8nDR//oX0nbFljizJ38obqi5ubEnbyJnl/vxBPf7dJs6/TZVvT8eif2JRQfVAo7cjVFc3vaP6cQfz0Vk9ecKPZ5W8/cxqQ1J/DknOIvwitvqbKc50SmuXav6FCcKbsAtLVPliAqL0wPR4sWAT/+CAweDDg6Fmxv1Ag4fdqMVSN95N/1J605rlvAxL+Nqjx5z5GhC88WbP95zyWcu1WxruekFsCF2xmoN3E9Hv9sm94hzT8fDoWdvVnQNskZOQCA9AeGL6sTfy0V3b7cYVQ97mbmFFvm8h3jh8jkgYjZyDatP56ESUaE4nyclE1kHqaHo+vXgZo1dber1YCK/zEtrbghmjwTP+Xk3fCGhsvk30bP3bTtYPTphjNm36daCHy6/gxUeQLnb2Vg5ML9OmWUJRyiemHRAZxKSjO6HoUlJGdi/6Xie6bWH7+Br/89Z7CH6MDle+j+5Q7svag9efzzuLP4ZrPl5j2RflfvZiH9gcrk43nxviu4bsTwKxEVzfR1jho0AHbsAMLDtbf/8QfQtKmZqkXGyO+ZkCvVnCMDXfJqO7p8iCUcv5aKE0mpmvsHLt8r8b7yF5Q8fSMNr/52yOC13/Q+V8+2DrO2AgA2jYtBzQBPvc+bteEMvtlyHgDQPLwy2tT009nf0z9Iw3DP/LgX0/pGIczXHdX9PTRLQwyNjoCPm9LoutqydfFJCKvijgbBPtauil4nElPR46ud6FQ3oERziJYfuIYxj9eyQM2IKg7Tw9GkScDQoVIPkloNrFwJnDkjDbf9/bcFqkhyh6+kFPm4qdlFPo/G0LCavDfK1J6p8mDQT/sMPnYnIxunktKNHs2MHL/WPJUq5ERiqsFwlB+MAOm9fPlME/RpUs3gEOwHq6Xh2jBfd8226/fu471V8VAAmNGvIbxcLRuU5u1MQKSfOzrW1b0YcGkcunIPL/92CABwaWYPs+7bXDYcvwEA+Pf0rRI9Pzs3z+BjD1R5cFU6GnyciCSmD6v16gUsWwasXSuN8UycCJw6Bfz1F9C5swWqSHL5p+yai76z1fLUAgcu3cUDVZ7mfr4KmI30upORjT+PXEfzaZswZN4+rHk45yjf5tOWuYhttkr/B5/KQLDVN4yWv4imvjWh5K7cLVgzp/tXO/DPsST8fSwJ09dadm7h8eupmPr3Sc26VLvPJ5s0l6ooZ28Yd2Hisrb34h08NWc3zpihfnlCIP2BCu+uOIad55I127/dch51J6zHC4sO8Kw2omKUbJ2jLl2AbduAjAwgKwvYuROIjTVz1chYJxJTte5fuG38vCB5b1F+UHph0QE89f0e9PtuNwCeHqxP82mbtFbqlmeQtfFJJV5wsiiHr6Qg6sN/sT3J8Pwm+TXgsnPz0PULw5O9Sxp0E5ItO+8sf0FMAJiz9QIG/bQPMZ9uxb3MnFKfVWfMs2+kPsC4ZUdwVHbGoKW99OtBHLh8D88t+K/0OxNSuy3dfxVD5hX0eubPX4o7eRO/7btc+tchKsdKvghkTg5w7Rpw5Yr2D1mMWk9ISb2vwtB52n9QN500vtciV8+p/JsfdueffDhRmOHINK88HLaxlBWXjBsWqfPBepy5abgnoqQ5Qx5eAOkYVKsFsnIMn5WnT3JGdrFh5+P1Bb1UTafG4fVlR5CUeh/T/j6pWQ36esp9oy+BctmIFaTfWn4UKw9fR59vd2ltT72vwo1Uyyx8m5IlncySmPrALGeDJhVTzxPXjTsJoCi/7r2MccuOICfXdtc9Iyop08PRuXNAu3aAm5s0KTsyUvqJiJD+JYvR96E79e+TOqd46/u4OXYtBQN+2KO1fg5Q6Gw1Ax8wFXEStq3bf+leoeFO039HJf2tXridqVns8vKdTDSfGocmH25E/YkbtMKMzuvJ6rj34h20mLYJ76w4ZtJrrz6SiOgZm/HTzgQ89ukWrDp8DY/O3IxXFxcfSLNz1fh+24Viy52XhROV7P9Ep9nb0HrGvzrh0NzWPZxzVBqFe70eFBqONTSkuuPcbfx55LrWYqaGfLD6OFYevo7aH6zDZ6VYgd6QGetOYcAPe0wO3UTmYHo4Gj4ccHCQJl8fPAgcOiT9HD4s/UsWs/6E7h/NwqdeA/p7BIbO+w/7Eu7imR/34PKdTCw/eE1a7Vo2rGaoh4jXJbU9g+btx0rZpSISkjPx046LRvWM5MsrxUrnV+5m4cyNdMR8uhW5aoG0h2s5zdkqhY/TN9Kw8lDBRXR/3n0JTafG4WRiGq6n3MczP+4FAPx+oHSXu3h92VEAwIYTUm9pUZORrxdaTFFfj0eeWmj1usxcVxD2kjOkUDRm6WF8K5vknm/FwWtYGy8tzJqZnWuwLhduZ+DDv07izyOmXwTZGFfuZuFisvYcrcILt/5+4BoSHpYRQmhWeB867z+MWXrE5OvtfWWBy9z8sO0i9iXcxdr40odFIlOZfrbakSNSKKpb1/y1IZPpWz334/Wn8fH607g0swdSsnLwzebzmoULH6jUiPl0KwDg4u0MhFQuOCNJpRZIzdJeqyorJxezNpp/7SAqvbeWF/S6fLe1+B4Rud3nk/HSryX/MlPcwpX5c53G/X4U/74Ro1nI8N2Vx3QC3I5zt+HlqoS3qxOq++s/484YKw9dw7jfj6JX42B8OaCJ1rUC1QLo8c1urfJjlx3Gd4Oba+4nptzXhJt883YmYELP+lpfHHZfuIPdF+6gd+NghD48o2/3+WS88YcU1E5M6YLWM/6Fp4sTdr/bUTMP7E5GNib/dVJz7TwACPByRXSNKiV+z/oUPstt+9nb+GWv7hyjDrO2In5yLNp/uhV3CvU+X7t3H5P/OoWMmwoEX02Bh6sLagd6wsnRvNcq33rmFn7dexnT+zVEgJcrhBBIyVKhskfB1RfuZOj21KVk5WDrmdvo0iAIbs6O2HL6FlyVjmZpy70X7yDQ2xWRfh4mPU8IYfTwLtk+08NR/fpAcnLx5cisklJNX9gt4t1/inz8u60X0Eb2xyRPLbDxpPa3tA6zthq1MjPZl6KWJyitm4XWbnphUcHk9GPXUgsX15ozd3F69xK/7rjfpXDy19FEXL6TiV9GtoKXixN6fbMbp286ofBA4tr4Gzh/KwOL9lxCTq4aS/dfNbhvfb1AqfdVCH14W96eHWZtRfqDXKQ/yEXag1woHRX4YJU0BFXYwLl7cerDrqa/2SIU7hF7dr7hSd6L913RCUb5fvvvKgBH/Pmj9PzejYPx5TNNcC9Lhd/0hK38NbzyXb2bhVkbz+CFdtURVU1aU2r3hWTM35mAiT0bIKyKO4YvkBZUFSviMW94S3z57zl8sUm7F2rGutMY1iZCawmCV347hN0X7mBI6zC81aUunnu4MOv3Q5oj0s8DS/67glEda8LP0wUAcP5WOlYfTsSLMdWLXIbi4u0MTa/m8DYRmNy7gcGygDRcqXR0gKODAs8t3I+tZ27jpboKmHoU30p/gLFLj+DpFqHo27Saic8mSzAuHKXJJu99/DHw9tvA9OnS9dSUhQ40b28zVk/Xd999h08//RRJSUlo0KABvvjiC7Rr186ir2kLzHGKrz7ypQHy1EJnrsHNNMvOr6Dyp9X0f7XuX7xt/Gn41d8zzzpQx66losOsrZgzuBlOF7Gq++OfbSt2X6OXHIafp7PO9ly1QOp9lc48q1uyOUlf/XsO83YmFLn/lYdLN7RYGjPWGb8sw5qjiajh74nPN+mfX5S/hte73erCTemo6S38+1gSzn/UDZk5eRg0VwqR+y/d07oMz7+nbyElK0cnGOU7fCUFIZXdCnrqHv7dWrb/KtrXDtCUe+nXg3B0UCBPLXDtXhZ+GtYSAPDcwv24evc+vtlyHgkzuiMx9QEys3NRO9BL63XirxeE94W7L2nC0cnENLg7O6KSuxKHrtyDl6sS/b+XFk4d3CoM0/pGYeuZ2wCA7087YpyJcwC/+vecpkfSVemIrlFBOHo1BRFVPODjrhvmHqjykKsW8HB2RK5aYOuZ22hYzQcbTtxAdX8PtKvlj39P3cS4348iwMsFvz3fCu+tioe/lytm9GuI2+nZ+PtYItrXCdD0kCWl3seUNScxrE2Epgcu7YEK7kpHnR7D/OHy/DB8NzMHLk4O8HAxvb/FVimEMTM5HRy0r1shhO51LPK35Rke8y+tZcuWYejQofjuu+/w6KOP4ocffsBPP/2EkydPIiwsrNjnp6WlwcfHB6mpqfA2Y4hTqVRYu3YtunfvDmXhsFiEnFw1HuTmwVv2TSYrJxcbT9zEhdsZGNe5tubgm/LXCSzYdclsdTbkgx71MM2GLyxLZCvqBnnhtI2um1RefTOoKTKzc/HOinijyp+d1g0rDl3D+JX6yz/fNhLNwyvj7eXHMDQ6XGd4es/4jnBXOqHxhxuLfJ1fRj6i1QPaIyoIQ6IjsO3sba2TAH56tgUerx+oWU7h4ycbYdvZ23h7ueETE5b+rzXqBXnjvioP/l4uyMjOReMpUn2ah1fGQT0r9ifM6K614Gykn4dmjtlfo9pi8l8nNM/b8XYH7Eu4izcfDgsDwI9Dm+PqvfuY+vdJAIDSUYGPn2yEvk2qYdeFZM17fb5tJOpV9dYMKQPA90Oa4dGafkYvFGvs56elPr8NMS4cbSv+25VGTEwpqlO0Vq1aoVmzZpgzZ45mW7169dC3b1/MmDGj2OdbqnH/PZmEkYsOoVGIN45dK+hlC6/iDrUQuHpXGhJrV8sPkX4eWLTnMr4a2BSvLTlc5H5Hto3E6RtpyMjOK7M1V+oEehV5+jcRUUXCv4kFoqp547gRy0DMfbYFOtc3bnV7Ww1HxvWByQPPlStAaKj+nqOrhsfsSysnJwcHDx7Eu+++q7U9NjYWu3fv1vuc7OxsZGcXdHGnPRweVKlUUJnxIrmfbJC6meXBCNBdU2XHuWTseLhibXHBCECx3fGWwD8CREQF+DexgDHBCAAquzka/RmbX6648ub8zDaG6QOEkZFAUhIQEKC9/e5d6TELDaslJycjLy8PgYHaaTQwMBA3bug/1XPGjBmYMmWKzvaNGzfC3d1dzzNK5vwtRwAluyo7ERFReXLywC5cN20JM8TFxRX5eFaW8cuUmIPp4UjffCNAupSIq6sZqlQ0RaHXLnyGhNz48eMxbtw4zf20tDSEhoYiNjbWrN1yLR7NRNtZOyEYkIiI7FL9ql44mWSbvUSuSgeMfDQC3269WOav/UpMdXy3zbTXfaZPN4Ofy4WpVCrExcWhc+fOxQ6rlSXjw1F+yFAogAkTAHnPS14esG8f0KSJeWsn4+fnB0dHR51eolu3bun0JuVzcXGBi4uLznalUmnSxOniBPh44ONH8tC9axd4ursiJ1eNCauPI8DbBd2iqiLCzx3uztpNfSvtAU4kpcFBocBjtfw0B9LF2xnYl3AXX2w6iw/7RKFLgyAA0qVDzHUmDxFRRdAivDIOyCYsPxLhi7e71sFTD880yze1bxRyctU4+XACslz+xPupfaMwtHV4sUukFNa1QRCq+3vgesp9/HmkYI2rD3rUg7+Xi9Y1GvW5ML07HB0UuJ+TpxWOqng4ay3D8EK7SMzdUTAVo3fjYJ0LYgNA7UBP1KvqjWv37mPgI2FaE7ELW/taO9QP9tYJR41DK6G6nwdWHb6OSu5KfDGgCXw9nDFpzQm0rekHZ2fdszuLU9znsjk/s41hfDg6/HCOjBBAfDwgf/POzkDjxsCbb5q5evKXcEbz5s0RFxeHJ554QrM9Li4Offr0sdjrGsvFEXB5uA6Hs5MDPn6qUZHlA7xdEeCt29NW3d8T1f09MfAR7bPv5AvaEZG2htV84OSowOErKQCANzrXxtL9V3E9xfT1wUwx8JEwLPnP8DUlB7cKw2/7ysc1J0d3rImvN+uuDC7Xo1FV/HOsYCHNN2Nr42JyJmr4e+LTDWfQr2k1XLmbpRVY8oVUdtO7qC0A/PNaW3i7KtHuky0AgNc61cKVO5lISM7EUdnaWYXP0lr+cht8vP60ZuX2uc+2gI+7Ei+3r6HZBgCDHglDRnau5uysfG1qVMGiEY/g2r37iHh4yvvUPg3w8foz+GJAE5y/naFZRf3i9O5aX2C7NAhEg2AfvNaplmZbRBUPfL35HH4a1gId60pf6k8lpePo1RR8P7Q57mXmwN3ZEWoBbDlzC080rQbHh3/73Zwd8cPQ5njxl4MY3iYCk3rVx9hlR7Dr/B2sG9MOVTyckZGdh7+PJuKrgU3xWG1/DGsTjifnFATB01O7aq0XlV/PMzfS8fexJBy7loJvBzdDVR83rTLrxrTDpxvOaK67+ePQ5vBxU2LEo5GIquat+XK/6pVH9f7+7JFxZ6vJPfcc8OWXFl/PSJ/8U/m///57REdH48cff8TcuXNx4sQJhIeHF/t8WzuV31SmfmMhMmRyr/qY/Jfut2RL+P3FaDz9w57iC5bCSzE1UK+ql+ZbeMKM7lAoFLh2LwttP96i9zlPNQ/BO13r4pvN5/Dznsvo16waejcO1ixMWJxfRj6CtjX9cF+Vh5xcNZp8qDtn4tSHXVFv4nqD+1gwvCWu3cvChD9PGPWaJVHD3wMXTFhrypCjk2Jx/lYG3lp+FJfvZOm93NBfo9qi1zc7NfcvzeyhuZ2apdKs2fP7/qt4W3Zdval9GqBv02poOFn/KfP5+8m/1p1Stu6O/O/ipZk98O+pmxj58wGt583fmYAqns7o00RaYDE7Nw+/7LmM5IwcjH28liYwXL2bhWHz/8PF5Ew4Oiiwfkw71Cq0FhIg9eTnf2H9+1giAr1d0TLCF1/GncGX/57DryMfQZtaATrPE0JaG6uSu+k9K/oUXm9In9w8NU4mpaF+VW+zr3BuDvZ9tprcggUWqIZxBgwYgDt37uDDDz9EUlISoqKisHbtWqOCEQE/DG0OBYDHavtjx7lkrZWLCytuDRcXJwdk82rcdstSfyTf6lIH/VuE4PytDLz4y0G8370e6sg+XB6r7Y/tZ29rPWdyr/ro3CAIAV4uuHI3C51mm7B0iEyYb8FQf/6HRUhld6x5JRqzVu7CwbtKZOYUnDCiAODv5YIpfaIwuXcDvR8wXw1sisvJmdiXcBc7z2tfGaBdLX8AgLuzE9ydgcfrBWDTKe1Ld7g5O2LXux1x5EoK5u28iOgaVZCTq9YMf3SoG4AVBy27CGR1f0+94WjH2x2Qel+FCX8e1/S4GbJh7GPwcVOieXhlbH6jPS7cztD5PT1eLxANQ3wwoEUolh24iv7NQ7Qely9m+HTLUDzdMhSFrRvTDo4OCvx9LEnv9dqUeo7bj56IwvurjuPzAY0BAB3rBuCDHvXQ8OGq3AAwoq32RdFdnBzxfLvqOvsK9XXH5jfb62kBbfKe/J6NgjW3X2lfHVXTT6NlRGW9z1MoFGYLRvn7K46TowMahVQy22tWFHa3nOUrr7yCV155xdrVsAtuSkfcl12NO3/+EoBi16Dwci360HB3dqzQ4Si2fiA2nrxp7WqUmNKxdMO0gd4u2PxGezSYtEFru6+HMwK8XBHg5YqjE2M1HyJ/vvooXJWOeKDK0wlHjo4OqFZJ6savUcJrq9Wr6oWmYZXx1cCmiKziofPYk5FqpDp5ag3ByD9X5B8yw9tEYOHuS9L79HJB78bBiL2Rjmd+3IN7WYZPJ/56YDOcSEzF8AX7kSFbab5aJTdUq+SGHo2qAgC+LLQKdOsaVaB0VECVZ1onvrGiq1fB9rO3df6/hvq6IxTSUEjDSRuQXmh1/BOTHse0Xzfg2e7tUCdIu/ekhr8nqvt7aK1+7vTwd/3xU40wtnMt+HvqzvcsTr2qUo/As9HOmnA0slCwKWzQI2Ho2SgYPm5S+FIoFHqDT1lxcSy+DNk+2+tjI5MF++g/S3Bk20gMaW145fBn9Hxzy9e9ofSH3N9L/x84fd/gKpKX29eAs5P9toGTQ+nqrhbQuVTAk81C8ITsulDyb9eNQyuhTpAXctW6gdqx0LffukG6wxj5hreJQIvwgm/lHs6OmNSrPno3lr69924cjIYhPnqf+1Gf+lr3FQbOLp3Uq6Cc8uHvuE6QFw5N6Ix/XmsLXw9nTNFzzS03Z0e0iPDFs9FST/YjEb569/9k82rwcnHCwEek/3/VKrlhx9sdtXrYzCmmjj8WPvdIkWV2vtMR/73XSXN/QItQODs5oHWAQM0A/YH1s6ebaN3v/jD8AUBVH7dS9U76ebrg0swe+O+9Tni/e70iyyoUCk0wIjIX+/3rThqLRur/w+fooMCLj9UAADQJraTzeFHXwRnaOhzfD2mOf15rq/dxI8/SLLe8XJ3wbte6Jj/PxYKBat6wFkaXDa7kVnyhIuR/YH38ZEOE+rph07jHMPvpxjqTPQur4qEbtqsUum7Z3GcNv48JPevjR9njRybF4rlHI40aXijc+2HoKQqFAq0ifRFexR1RwT5a2xsE++DgB49jWJsIg6/zWqda+HZQM4PvI6SyOw5O6IwZ/QpO2gjyccUXzzRBVLXSz6VwVRYcY90bBqGGvydqBWoHnAk9tYOij7sSAd6uWD+2HV5pXwPv9yw6kADS35Sd73TAiEcj8e2gZpqAak4B3q48GYWswvRhtcxMwMOj+HJkdq2r+2Lvxbs62x0M/JWvV9Ubob7uODoxFh56+no9nA1/kDk5OqBrVJDBxys6Z0fDbffpU43wloFrJQ1tHY6fSrnyebeQPASERqJBcCW8s/IYhACcHR3QqZ5xy/W/1aUOWlfX36thjMfrBWiuHD6gZRgGtCz+uob5Ivw8MKt/Y1TxdMbdjBycSExD50L1DvV118xb6dGwKv6Jl85++qBHPTg6KODr4YzPnm4MTxenUvVgFpWnlrzQGmoh9PZ+FBfEXJWOmiE0Q/T1Otar6o2/R7fDicRU9Phqp55nGSfu9RhcuJ2BdfE3NBdO9fN0wbDocPy85zIAw0NVdYO8UberFNCMWZE4pLI7JvaqX2w5IntjejgKDASefhoYMQJoq79XgSwjx8Acn8Lh6O/RbRF/PRVdGkgfOvqu6gwA7qW4gnLT0MpYf0L/yuQVgZ+Xc5E9D5882UjrbJx89YO9sX5sO2Rm5+HJOfove2NI/+YhmNijDjbHbUD37nWhVCrxeP1A/LD9Avo3NzxEKnd0UmyRQxD+Xi64LbuqvD4/DDW+h0qfp2QTdZ9srr/M+z3r4Ylm1dAq0hf/u1YdO87dxnBZb02/ZiH6n2iClgaGvQBpSNDBSou6NgjWPyxorFBfd4T6uqN9He2zpZqGVdaEIyIqmulfu5YsAVJTgU6dgNq1gZkzgUTdhabI/FL0TAad3b8x5NM4awd6IqqaDwY+ElbsN1xnA9+6i5sA+V73uvD1NN8ZF7auVaT2h+jQ1uFwd3bSTEAtTOmoQJ+mhocY6gZ5o3l4Zax8pY1J9fj4yUY6w1a+Hs4Y362eZl6IfJ7LnvEd8f2QZlrl5cFI31y17wY3Q3T1Kjrb3+pSBwDQOMRHs+6KJXm7KtG6ehUoFAo0Dq2EUR1rmeUMO/mwZv5p3bbo8Xq6p4GXVs9GVTG8TQTmDG5WfGGiCs70vza9egErVkiB6OWXpbAUHg707AmsXAnk5ha/DyqRwEKLRq4Z9SiebB6i9f12+cumfeDq82JM0Wd6BHq76kyirUiahVcCAIPza2LrB8HFyVEnVOnsJ0z/6b6GGDP3Yu6wFvh+SHOcmdYVVX3c0KVBEL4e2BQ9GlbF7y9Ga5Vd+r9onefXDfLCkv+11tn+aoeaODutm90v8vbX6LZ4vm0kDk3oXCYhr6RmP90Erz9e26z7dHJ0wOTeDdCtYdFDfkRUmgnZVaoAr78OHD0KfPYZsGkT8NRTQHAwMHEiUMYXiasIqvsXzPX67flWmrUr5DnFUG+GPoaKBngVfY08JwcHm/5gsQQnPeua9GocjLY1/TS9KgAwLDocbg/ncpm4vKpZ+Lgp0TVKCmeANMTXq3Ewvh3cDI8UCmthVXQvvuz2MPDp67lwdnKw+8mxtQO98EHP+vD1sO2eTx83JXo2ZoghspaSTzq5cQNYtEhaFPLKFSkYjRwp9SjNnAns3Qts1L/aKZWMfBKnvBdJfkqyodOT9fE28vRXB4V06nY+J0dFhQtH/74RgwOX7uGJptU0AcFV6Yhfn28FQBpKW37wGkbLLhXQtpYf/rt0V2e9qZKY0a9hqZ5vjD9eitYMXc19tgXWHE3EmKVH8EGP4s9cIvMzNMeQiCzP9HC0cqUUiDZsAOrXB159FRgyBKhUqaBMkyZA06ZmqyRJXo6pgQW7LgEovIBdyfbXs1Ewxv1+FH2bBGP1EcPzxpb+LxrD5v+n+YBXVsBwFF7FA+FVDJ+l+b/HauB/D5dNKNhWHZU9nNG+tr/mmlAl6U3ydnXSudaeuT0S4as1QVmhUKBPk2roWDcAXq5cQ8YabL13i6g8M31Y7bnngGrVgF27gCNHgFGjtIMRAFSvDrz/vlkqSAUqy/5Yyj9k5X9ETVn52NnJAZdm9sAXzxQEWX3Dco9E+moteufo4GBw+QBrK25l75Io6eiYq9IRQ1uHI9RXd/gqn7WbMf8MsJc71ND7OIOR9QR6u5q0dhURmY9x4ahZM+Dew6sov/iiNMeoZUvD5d3cgEmTzFA9kpNPgpYvze/h4oT1Y9th07jHSnxGT/2Hy/YbOoNHPtdE6aCArS6QrW/CcGkDiDk7yQovTH3wg854VU8wKXyWmaVM6Fkfe8Z3RIc65j87ikrP2LWr8uk7lojIdMZ9zT51Slr8sXJlYPZs4K23AHfD34bJMhwcFFg/th2yVWqdtYvqBpVuZd2Fz7XEhhM3DK4fIw9Djg4Kmz1bTd8aPq+0r4Fvt1wo8nljOtXCl3oudAkA7xVz+QJTFO5x8/VwRq0A3ctGyC/vEWTg8jDm4OigQFWf0q2WTdb32dON0bCaj94ryBOR6YwLR02aSMNpbdtK4zmffgp4GrhA5MSJ5qsd6ShtCDIkwNsVQ6MjDD4u/1B3cnSAYymvzWUp+uZC1TbiAyPA2/BFMs15RWt9p/93jQpCy32V0bp6FXy9+TwA6VIXi0Y8gq/+PYeZTzbSeQ5VHG/G1sasjWeLLGOORTGJqIBx4WjhQmmY7O+/pTGKdesAJz1PVSgYjsopeehQOiq0rt9kS/T1aHWqF4iPnohCVLAP+ny7S+/zDC2IaS5jH6+FE4lp6FRXd/jKVemIP16S1qfq27QaUrJyNKscP1bb36L1Its3qmMtrXDUvo4/tp65bcUaEZV/xoWjOnWApUul2w4OwL//AgGco1CRyHuOHB0UxV5g1Fr0dWg5KhQY3Cq8yOdZejhirJEL+tXwN9AjSwTphIm2Nf20wtG0vlFWrBFR+WT6qT1qrr1REcnDkdLRwXZ7jvQMq+mbHtWwmg/ir6eiflVvvPBYJJqEVrJ85YhKKUzPmY9DWhcd/InIdMaFozVrgG7dAKVSul2U3r3NUC2yNfLQYdM9R3qSkL4hs8dq++HXka3g7eZU7DXoiKxt+UvR+HrzeUzsVR//nrpp7eoQlXvGhaO+faUVsQMCpNuGKBRAXulWAibbJM8XSgcHzeUpbEENfw9cuJ0JQH840nfJCz9PF50z/ohsVYsIX/w84hEAwP0c/o0lsjTjxkbU6oI5Rmq14R8Go3JL+2w125qQ/clTjTW35T1cldyVOle+n9GvIbo0CMSgVpZdcZrIUqKq+WDs47WKL0hEJWb+5YSpXJKHDicbG1YTsuXC5Z1E7Wv761z5fuAjYQYvxfHj0OZ4e8UxpGSpNNs4F4ls0WsdayE3T6BRiI+1q0JULpn+9f+114CvvtLd/s03wNixpa8R2TwnRwfbCkey2/L5Q6bOJYptEIQBLUI194e3icCcMlqpmsgUDg4KvNmlDmIbBFm7KkTlkunhaMUK4FHdSzSgTRtg+XIzVIlskfxabtKEbOsOq/VoVFVzu7aeFaYBoCTTrOVBa3LvBlw9moioAjJ9WO3OHcBHT1eutzeQnGyGKpEtkocGpaMCrlackP1qhxpQymaI+7gr8d/7neBWqDerJGeh8cQ1IiIy/et/zZrA+vW629etA6pXN0OVyBapZV1HTg7WHVZTQKHVkwUAAV6uOleQD6lseq/PyEcjEeDlgpfb8wKeREQVlek9R+PGAaNGAbdvAx07Stv+/Ve6IO0XX5i3dmQzhFY4su6wWnG9Owuea4mNJ26WKOAEeLti33uduPYREVEFZno4GjECyM4GPvoImDpV2hYRAcyZAzz7rHlrRzZDvjC6g4MCzk62cyp/YR3qBKBDnZJf3obBiIioYivZqfwvvyz93L4NuLkBnrweVHmnLjSOZekLtRaF0YWIiCypZJ9wubnApk3AypUFpzElJgIZGWasGtkyJyuGI86aJiIiSzK95+jyZaBrV+DKFWl4rXNnwMsL+OQT4MED4PvvLVBNsrYWEb5wcXJAeBXdC1+WNQW0z54jIiIyJ9PD0ZgxQIsWwNGjQJUqBdufeAJ4/nkzVo1sia+HM3a83QFuzraz+CMREZElmB6Odu4Edu0CnJ21t4eHA9evm6laZIsCvF2tXQUA0qhabp66+IJEREQlYPrEEUMXmL12TRpeI7IwBRRoFVml+IJEREQlYHo46txZez0jhUKaiD1pEtC9u/lqRlSEdrX88OvIVtgzvqO1q0JEROWM6cNqn38OdOgA1K8vTcAeNAg4dw7w8wOWLLFAFYl0KRQKtK3lZ+1qEBFROWR6OAoOBo4ckYLQoUPSMNvIkcDgwdKaR0QWxjP5iYjIkkq2CKSbm7RS9ogRZq4OUfGYjYiIyJJKFo7OnAG+/ho4dUr6Gl+3rnS9tbp1zVw9Il2V3JXFFyIiIioh0ydkL18OREUBBw8CjRsDjRpJw2sNGwJ//GGBKhJp698i1NpVICKicsz0nqO33wbGjwc+/FB7+6RJwDvvAP37m6lqRLp6NKoKVyUXoiQiIssxvefoxg3g2Wd1tw8ZIj1GREREZMdMD0ft2wM7duhu37kTaNeu9DUiu/HryFZoHOJj7WoQERGZlenDar17S8NnBw8CrVtL2/buleYbTZkCrFmjXZbKrba1/NC6ehvUfH+dtatCRERkNqaHo1dekf797jvpR99jgHQWm77LjBCVhrB2BYiIqLwzPRypecFPKqDgioxERFTOmD7niEiG0YiIiMob48PRvn3AukJzSxYtAiIjgYAA4H//A7KzzVw9IiIiorJlfDiaPBk4dqzgfny8dE21xx8H3n0X+OsvYMYM89eQbBpH1YiIqLwxPhwdOQJ06lRwf+lSoFUrYO5cYNw44KuvgN9/N38NyaaV9ZwjwRnZRERkYcaHo3v3gMDAgvvbtgFduxbcb9kSuHrVjFUj0sUJ4EREZGnGh6PAQCAhQbqdkyNdTy06uuDx9HRAyQuCkmU5MhwREZGFGR+OunaV5hbt2CFdW83dXXtF7GPHgBo1LFBFogJVK7lauwpERFTOGb/O0bRpQL9+QEwM4OkJ/Pwz4Oxc8Pj8+UBsrAWqSATMfbYF/j6WiNEda1m7KkREVM4ZH478/aVeo9RUKRw5Froy+h9/SNuJLKBz/UB0rh9YfEEiIqJSMn2FbB8DFxr19S1lVYiIiIisjytkExEREckwHBERERHJMBwRERERyTAcEREREckwHBERERHJMBwRERERyTAcEREREcnYTTj66KOP0KZNG7i7u6NSpUrWrg6VoQXDW1q7CkREVIHYTTjKyclB//798fLLL1u7KlSGujQIRIe6AdauBhERVSCmr5BtJVOmTAEALFy40LoVISIionLNbnqOqGJSQGHtKhARUQVjNz1HJZGdnY3s7GzN/bS0NACASqWCSqUy2+vk78uc+ySJWqjZrjI81kzHNisZtpvp2GamM7bNyrpNFUIIUaavKDN58mTNcJkh+/fvR4sWLTT3Fy5ciLFjxyIlJaXE+1+8eDHc3d1Nri/pN2aP5TJ2I181RtZRW2z/RERk+7KysjBo0CCkpqbC29vb4q9n1XCUnJyM5OTkIstERETA1dVVc9+UcKSv5yg0NBTJyclmbVyVSoW4uDh07twZSqXSbPu1F7UmbLTYvmPrB+DbgU0stn97U9GPtZJgm5UM2810bDPTGdtmaWlp8PPzK7NwZNVhNT8/P/j5+Vls/y4uLnBxcdHZrlQqLXLgWmq/9mrJC60xcO7eUu3DQeHANtWDx5rp2GYlw3YzHdvMdMW1WVm3p93MObpy5Qru3r2LK1euIC8vD0eOHAEA1KxZE56entatHOkVXaOKtatARERkMrsJRxMnTsTPP/+sud+0aVMAwJYtW9C+fXsr1YqIiIjKG7s5lX/hwoUQQuj8MBiVbwqeyU9ERGXMbsIR2acvn2lSqudb73QBIiKqqBiOyKL6NKlm7SoQERGZhOGILM7JoeRjYxxWIyKissZwRBa3dkw7a1eBiIjIaAxHZHG1A72sXQUiIiKjMRwRERERyTAcEREREckwHBERERHJMBwRERERyTAcEREREckwHFGZqB3IiwMTEZF9YDiiMrH4hdbWrgIREZFRGI6oTPh5uli7CkREREZhOCKbxsuHEBFRWWM4IpvmwHRERERljOGIbFpwJTdrV4GIiCoYhiOyaa91qmXtKhARUQXDcEQ2zdPFydpVICKiCobhiIiIiEiG4YiIiIhIhuGIiIiISIbhiIiIiEiG4YgsollYJWtXgYiIqEQYjsjs2tXyw0/DWlq7GkRERCXCcERmN+LRSPh6OFu7GkRERCXCcEREREQkw3BE5sfLoRERkR1jOCKb5cCQRUREVsBwRGZnrkzjoGA6IiKissdwRGanKGWoebpFCABgDC86S0REVsCrepLNmdGvEUa0jUSdQC9rV4WIiCoghiMyu9IOhjk6KFA3yNssdSEiIjIVh9WIiIiIZBiOyOw4j5qIiOwZwxGZnYILHRERkR1jOCIiIiKSYTgis+OwGhER2TOGI7IpK16OtnYViIiogmM4IrMrTcdRgJer2epBRERUEgxHRERERDIMR2R+nHNERER2jOGIzI6n8hMRkT1jOCIiIiKSYTgiIiIikmE4IrPjOkdERGTPGI7IpjBYERGRtTEckdmVJt84MB0REZGVMRyRTfFydbJ2FYiIqIJjOCKzU5Sw9+fXka3g5ao0c22IiIhMw3BENqFlRGW0reVn7WoQERExHJH5cdoQERHZM4YjMruSZCOuqk1ERLaC4YiIiIhIhuGIiIiISIbhiMyOc46IiMieMRwRERERyTAckQWUoOuIvU1ERGQjGI7I7DisRkRE9ozhiGwC8xQREdkKhiMyu5IEnSAfV7PXg4iIqCQYjqjUXoyprnXfoQTjah/0qG+u6hAREZWKXYSjS5cuYeTIkYiMjISbmxtq1KiBSZMmIScnx9pVIwDvdq2LQxM6a+47OZoWjqr7ucPfy8Xc1SIiIioRJ2tXwBinT5+GWq3GDz/8gJo1a+L48eN44YUXkJmZiVmzZlm7ehWeQqGAr4czhkWHIyn1AepX9bZ2lYiIiErMLsJR165d0bVrV8396tWr48yZM5gzZw7DkQ2Z0ieqyMcfq+2P7Wdvl1FtiIiISsYuhtX0SU1Nha+vr7WrQSb4ZlBTa1eBiIioWHbRc1TYhQsX8PXXX2P27NlFlsvOzkZ2drbmflpaGgBApVJBpVKZrT75+zLnPssjN0dA6aiAKk9obReCbWcsHmumY5uVDNvNdGwz0xnbZmXdpgohhCi+mGVMnjwZU6ZMKbLM/v370aJFC839xMRExMTEICYmBj/99FOJ9r948WK4u7uXrNJUKuP2OiJPaE/YDnAVeL9pnpVqREREti4rKwuDBg1CamoqvL0tP6/VquEoOTkZycnJRZaJiIiAq6u0Bk5iYiI6dOiAVq1aYeHChXBwKHpUUF/PUWhoKJKTk83auCqVCnFxcejcuTOUSqXZ9lse1Z8cp9NzFFnFHRvHtrVSjewLjzXTsc1Khu1mOraZ6Yxts7S0NPj5+ZVZOLLqsJqfnx/8/PyMKnv9+nV06NABzZs3x4IFC4oNRgDg4uICFxfdU8SVSqVFDlxL7bciYLuZhsea6dhmJcN2Mx3bzHTFtVlZt6ddzDlKTExE+/btERYWhlmzZuH27YIznoKCgqxYMzIHXouNiIhsiV2Eo40bN+L8+fM4f/48QkJCtB6z4qggmQ3TERER2Q67OJV/+PDhEELo/SH7x54jIiKyJXYRjqj80JdnmY2IiMiWMByR1bHniIiIbAnDEVmdgn1HRERkQxiOqExxlhgREdk6hiOyOg6rERGRLWE4IiIiIpJhOCKrY8cRERHZEoYjKlN616biuBoREdkQhiOyOkYjIiKyJQxHZHXsOCIiIlvCcERlSt+p/AxHRERkSxiOiIiIiGQYjsjquEI2ERHZEoYjspp2NasAAJ5tHWblmhARERVgOKIyJT+T/4chTTG+cS76NqlqvQoREREVwnBEVqN0dECQO6DgjGwiIrIhDEdEREREMgxHVKacnXjIERGRbeMnFZWp355vhfAq7ljwXEtrV4WIiEgvJ2tXgCqWlhG+2PZWBwCASqWycm2IiIh0seeIiIiISIbhiIiIiEiG4YiIiIhIhuGIiIiISIbhiIiIiEiG4YiIiIhIhuGIiIiISIbhiIiIiEiG4YiIiIhIhuGIiIiISIbhiIiIiEiG4YiIiIhIhuGIiIiISIbhiIiIiEjGydoVKEtCCABAWlqaWferUqmQlZWFtLQ0KJVKs+67PGO7mY5tZjq2Wcmw3UzHNjOdsW2W/7md/zluaRUqHKWnpwMAQkNDrVwTIiIiMlV6ejp8fHws/joKUVYxzAao1WokJibCy8sLCoXCbPtNS0tDaGgorl69Cm9vb7Ptt7xju5mObWY6tlnJsN1MxzYznbFtJoRAeno6goOD4eBg+RlBFarnyMHBASEhIRbbv7e3N/9DlADbzXRsM9OxzUqG7WY6tpnpjGmzsugxyscJ2UREREQyDEdEREREMgxHZuDi4oJJkybBxcXF2lWxK2w307HNTMc2Kxm2m+nYZqaz1TarUBOyiYiIiIrDniMiIiIiGYYjIiIiIhmGIyIiIiIZhiMiIiIiGYYjM/juu+8QGRkJV1dXNG/eHDt27LB2lSxi8uTJUCgUWj9BQUGax4UQmDx5MoKDg+Hm5ob27dvjxIkTWvvIzs7G6NGj4efnBw8PD/Tu3RvXrl3TKnPv3j0MHToUPj4+8PHxwdChQ5GSkqJV5sqVK+jVqxc8PDzg5+eH1157DTk5ORZ778bavn07evXqheDgYCgUCqxevVrrcVtro/j4eMTExMDNzQ3VqlXDhx9+WGbXLspXXJsNHz5c57hr3bq1VpmK1mYzZsxAy5Yt4eXlhYCAAPTt2xdnzpzRKsNjTZsxbcZjTdecOXPQqFEjzSKN0dHRWLdunebxcnucCSqVpUuXCqVSKebOnStOnjwpxowZIzw8PMTly5etXTWzmzRpkmjQoIFISkrS/Ny6dUvz+MyZM4WXl5dYsWKFiI+PFwMGDBBVq1YVaWlpmjIvvfSSqFatmoiLixOHDh0SHTp0EI0bNxa5ubmaMl27dhVRUVFi9+7dYvfu3SIqKkr07NlT83hubq6IiooSHTp0EIcOHRJxcXEiODhYjBo1qmwaoghr164V77//vlixYoUAIFatWqX1uC21UWpqqggMDBTPPPOMiI+PFytWrBBeXl5i1qxZlmsgPYprs2HDhomuXbtqHXd37tzRKlPR2qxLly5iwYIF4vjx4+LIkSOiR48eIiwsTGRkZGjK8FjTZkyb8VjTtWbNGvHPP/+IM2fOiDNnzoj33ntPKJVKcfz4cSFE+T3OGI5K6ZFHHhEvvfSS1ra6deuKd99910o1spxJkyaJxo0b631MrVaLoKAgMXPmTM22Bw8eCB8fH/H9998LIYRISUkRSqVSLF26VFPm+vXrwsHBQaxfv14IIcTJkycFALF3715NmT179ggA4vTp00II6cPUwcFBXL9+XVNmyZIlwsXFRaSmpprt/ZZW4Q96W2uj7777Tvj4+IgHDx5oysyYMUMEBwcLtVptxpYwnqFw1KdPH4PPqehtJoQQt27dEgDEtm3bhBA81oxRuM2E4LFmrMqVK4uffvqpXB9nHFYrhZycHBw8eBCxsbFa22NjY7F7924r1cqyzp07h+DgYERGRuKZZ57BxYsXAQAJCQm4ceOGVlu4uLggJiZG0xYHDx6ESqXSKhMcHIyoqChNmT179sDHxwetWrXSlGndujV8fHy0ykRFRSE4OFhTpkuXLsjOzsbBgwct9+ZLydbaaM+ePYiJidFafK1Lly5ITEzEpUuXzN8ApbB161YEBASgdu3aeOGFF3Dr1i3NY2wzIDU1FQDg6+sLgMeaMQq3WT4ea4bl5eVh6dKlyMzMRHR0dLk+zhiOSiE5ORl5eXkIDAzU2h4YGIgbN25YqVaW06pVKyxatAgbNmzA3LlzcePGDbRp0wZ37tzRvN+i2uLGjRtwdnZG5cqViywTEBCg89oBAQFaZQq/TuXKleHs7GzT7W5rbaSvTP59W2rHbt264bfffsPmzZsxe/Zs7N+/Hx07dkR2djYAtpkQAuPGjUPbtm0RFRWlVRcea/rpazOAx5oh8fHx8PT0hIuLC1566SWsWrUK9evXL9fHmZNJpUkvhUKhdV8IobOtPOjWrZvmdsOGDREdHY0aNWrg559/1kxaLElbFC6jr3xJytgqW2ojfXUx9FxrGTBggOZ2VFQUWrRogfDwcPzzzz/o16+fwedVlDYbNWoUjh07hp07d+o8xmNNP0NtxmNNvzp16uDIkSNISUnBihUrMGzYMGzbtq3Ietr7ccaeo1Lw8/ODo6OjTiK9deuWTnotjzw8PNCwYUOcO3dOc9ZaUW0RFBSEnJwc3Lt3r8gyN2/e1Hmt27dva5Up/Dr37t2DSqWy6Xa3tTbSVyZ/CMGW27Fq1aoIDw/HuXPnAFTsNhs9ejTWrFmDLVu2ICQkRLOdx5phhtpMHx5rEmdnZ9SsWRMtWrTAjBkz0LhxY3z55Zfl+jhjOCoFZ2dnNG/eHHFxcVrb4+Li0KZNGyvVquxkZ2fj1KlTqFq1KiIjIxEUFKTVFjk5Odi2bZumLZo3bw6lUqlVJikpCcePH9eUiY6ORmpqKv777z9NmX379iE1NVWrzPHjx5GUlKQps3HjRri4uKB58+YWfc+lYWttFB0dje3bt2udCrtx40YEBwcjIiLC/A1gJnfu3MHVq1dRtWpVABWzzYQQGDVqFFauXInNmzcjMjJS63Eea7qKazN9eKzpJ4RAdnZ2+T7OTJq+TTryT+WfN2+eOHnypBg7dqzw8PAQly5dsnbVzO6NN94QW7duFRcvXhR79+4VPXv2FF5eXpr3OnPmTOHj4yNWrlwp4uPjxcCBA/We0hkSEiI2bdokDh06JDp27Kj3lM5GjRqJPXv2iD179oiGDRvqPaWzU6dO4tChQ2LTpk0iJCTEJk7lT09PF4cPHxaHDx8WAMRnn30mDh8+rFnawZbaKCUlRQQGBoqBAweK+Ph4sXLlSuHt7V3mpwoX1Wbp6enijTfeELt37xYJCQliy5YtIjo6WlSrVq1Ct9nLL78sfHx8xNatW7VOO8/KytKU4bGmrbg247Gm3/jx48X27dtFQkKCOHbsmHjvvfeEg4OD2LhxoxCi/B5nDEdm8O2334rw8HDh7OwsmjVrpnVqaHmSv36FUqkUwcHBol+/fuLEiROax9VqtZg0aZIICgoSLi4u4rHHHhPx8fFa+7h//74YNWqU8PX1FW5ubqJnz57iypUrWmXu3LkjBg8eLLy8vISXl5cYPHiwuHfvnlaZy5cvix49egg3Nzfh6+srRo0apXX6prVs2bJFAND5GTZsmBDC9tro2LFjol27dsLFxUUEBQWJyZMnl/lpwkW1WVZWloiNjRX+/v5CqVSKsLAwMWzYMJ32qGhtpq+9AIgFCxZoyvBY01Zcm/FY02/EiBGazzd/f3/RqVMnTTASovweZwohyni5TSIiIiIbxjlHRERERDIMR0REREQyDEdEREREMgxHRERERDIMR0REREQyDEdEREREMgxHRERERDIMR0TWsHUroFAAKSnWronxIiKAL74ouszkyUCTJmX3eraqfXtg7Fhr10KXQgGsXm3Z17DHY5uoEIYjIrnhw6U/7AoF4OQEhIUBL78MFLpoYoW0fz/wv/8V3DfXB+3ChUClSsW/nj1ZuRKYOtXatbCONm2ApCTAx0e6b+j3S2TDnKxdASKb07UrsGABkJsLnDwJjBghfQtessTaNbMuf//y/Xrm5Otr7RpYj7Mz8PBq7UT2ij1HRIW5uEh/3ENCgNhYYMAAYONG7TILFgD16gGurkDdusB33xU8dumS1KuydKn0LdrVFWjQQBpuMOTOHWDgQOk13d2Bhg11w5haDXz8MVCzplTHsDDgo48KHr9+Xapr5cpAlSpAnz5SXQxp3hyYPbvgft++Um9ZWpp0/8YN6X2cOSPdlw9z5V/h+oknpDKFr3j9yy/SNh8f4JlngPR0/XXYuhV47jkgNbWgx27yZN3XA6THfvgB6NlTaqN69YA9e4Dz56VhLA8PIDoauHBB+zX++kt6r66uQPXqwJQpUvA1ZP9+oHNnwM9Pqn9MDHDokHadnZ2BHTsKts2eLZXPv2J44WG1774DatWS6hAYCDz1lOHXN+ZYaN8eeO014O23pSAWFFTQbvnOnQMee0x6zfr1AdlV0Q3SN5TZpIn2vhUK4KefpN+9u7v0vtasKXhcPqxW1O/XlDYhKmMMR0RFuXgRWL8eUCoLts2dC7z/vhRMTp0Cpk8HJkwAfv5Z+7lvvQW88QZw+LAUknr3lj749HnwQPoA//tv4PhxaThp6FBg376CMuPHS+FowgSpR2vxYulDBQCysoAOHQBPT2D7dmDnTul2165ATo7+12zfviCwCSF92FeuLD0XALZskT5069TRfe7+/dK/CxZIgSD/PiCFk9Wrpffy99/Atm3AzJn669CmjfRh7O0t7ScpCXjzTf1lAWmo6tlngSNHpFA6aBDw4otS2xw4IJUZNaqg/IYNwJAhUpA4eVIKVwsXaofKwtLTgWHDpPbYu1f6AO/evSDg5QefoUOlD/2jR6XjYe5coGpV3f0dOCC9/ocfSkFz/XoptBhizLEASMebh4e0/ZNPpP3nByC1GujXD3B0lN7D998D77xj+DVNNWUK8PTTwLFjUtsMHgzcvatbztDv19Q2ISprJl+qlqg8GzZMCEdHITw8hHB1FUKKDUJ89llBmdBQIRYv1n7e1KlCREdLtxMSpOfMnFnwuEolREiIEB9/LN3fskUqU+iq01q6dxfijTek22lpQri4CDF3rv6y8+YJUaeOEPKrT2dnC+HmJsSGDfqfs2aNED4+QuTlCXHkiBD+/kK8/roQb70lPf6//wkxYEBB+fBwIT7/vOA+IMSqVdr7nDRJCHd3qb753npLiFatDL/PBQukehSm7/U++KDg/p490rZ58wq2LVki/d7ytWsnxPTp2vv95RchqlY1XJ/CcnOF8PIS4q+/CrZlZwvRtKkQTz8tRIMGQjz/vPZzYmKEGDNGur1ihRDe3tptYir5sZC//7Zttcu0bCnEO+9ItzdskI7jq1cLHl+3Tv/vTK5wmwshROPG0u81X+HfQ0aGEAqFtH8hdI9tfb9fc7QJkQVxzhFRYR06AHPmSL0xP/0EnD0LjB4tPXb7NnD1KjByJPDCCwXPyc0tmICaLzq64LaTE9CihdTTpE9entS7smyZNDyWnS39eHhIj586Jd3v1En/8w8elIaXvLy0tz94oDvMlO+xx6TekMOHgV27pOGjDh2AadOkx7duLdkZVxER2vWoWhW4dcv0/ejTqFHB7fxes4YNtbc9eCANDXp7S+2yf792T1FenlQmK0saFirs1i1g4kRg82bg5k2pfFYWcOVKQRlnZ+DXX6X6hIcXfVZd585SmerVpZ68rl0LhqT0Ke5Y0NcWgHY7nzolDbuGhBQ8Lj8eS0v+2h4e0u/blN+xqW1CVMYYjogK8/CQ5vUAwFdfSYFhyhRpSEetlrbPnQu0aqX9PEfH4vetUOjfPns28Pnn0odsw4ZSHcaOLRgSc3Mrer9qtTQU89tvuo8Zmtjs4yPNJ9m6Fdi9G+jYEWjXThqyOndOCoXt2xf/ngqTD0EC0nvOb7fSku87vy31bct/PbVa+t3166e7L1dX/a8xfLgUgr/4QvoAd3GRgkXh4cndu6V/796VfgqHl3xeXtKcpa1bpblrEydK827279d/Fldxx0K+otpZCN39Gjr25BwcdJ+rUumWK+3v2NQ2ISpjnHNEVJxJk4BZs4DERKlnolo1aS5SzZraP5GR2s/bu7fgdm6u1ItRt67+19ixQ5pAPWQI0Lix9I363LmCx2vVkgLSv//qf36zZlL5gADdehXu0ZJr316aW7R9u3S7UiVp8u60adK+6tUz/FylUurlKC1nZ/PsR59mzaQ5LYXbpGZNKQjos2OHNB+me3dpIr2LC5CcrF3mwgXg9delkNy6tTQPqqhw4OQEPP64NDfo2DFpovzmzYZfv6hjwRj160s9XYmJBdv27Cn+ef7+BZPKAakHLiHBtNcuzNDv15Q2ISpjDEdExWnfXvqQnD5duj95MjBjBvDll1LvSny8NDH5s8+0n/ftt8CqVcDp08Crr0prJY0Yof81ataUJtPu3i0Nibz4onS2WD5XV2lC7dtvA4sWSR/Oe/cC8+ZJjw8eLJ0t1aeP9OGakCBNhB4zBrh2rej3tn699M2/fv2Cbb/9Jg2zFSUiQgprN26Ubh2oiAggI0PaV3KyNIRlLhMnSu01eTJw4oTUtsuWAR98YPg5NWtKZ9udOiVNdh48WLvnLi9PmiAdGyudibVggTRxWn7mn9zff0s9kEeOAJcvS/VRq/VPdM9//aKOBWM8/ri0/2eflSaM79ghTRovTseO0nvfsUN6T8OGGdcjWhR9v19T24SojDEcERlj3Dipl+DqVeD556W5SAsXSsMeMTHS7cI9RzNnSmeXNW4sfdj8+acUYPSZMEHq5ejSRQonQUHSqfWFy7zxhvSBX6+edNp+/jwPd3ep9ycsTBpCqldPCmL370tzbwzJP0MoJqZg2CUmRgoAxYWj2bOlD/HQUKBp06LLFqVNG+Cll6T34+8v9SSYS5cu0gdxXBzQsqXUy/PZZ9JwmSHz50thr2lTKQS99prUi5bvo4+kXo4ff5TuBwVJx8MHH0gf9oVVqiQtCtmxo/R7+f576dT8Bg30v74xx0JxHBykYJ6dDTzyiHTMFnWGXr7x46VjomdPqeesb1+gRg3TXrswfb9fU9uEqIwphNA3OE1EJXbpkhSUDh8236U0iIiozLDniIiIiEiG4YiIiIhIhsNqRERERDLsOSIiIiKSYTgiIiIikmE4IiIiIpJhOCIiIiKSYTgiIiIikmE4IiIiIpJhOCIiIiKSYTgiIiIikmE4IiIiIpL5P7cUn55nynlTAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Select a seismic trace\n",
    "trace = stream[1]\n",
    "\n",
    "# Waveform data\n",
    "waveform = trace.data\n",
    "\n",
    "# Plot\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(waveform)\n",
    "\n",
    "# Axes\n",
    "ax.set_xlabel(\"Replace with time axis and units\", color=\"red\")\n",
    "ax.set_ylabel(\"Specify the units\", color=\"red\")\n",
    "ax.set_title(trace.id)\n",
    "ax.grid()\n",
    "\n",
    "# Show\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4. Saving and reading waveforms locally\n",
    "\n",
    "Of course, in order to avoid re-downloading the waveforms everytime, one should be able to write it down to a local or remote disk. This can be achieved at the `stream` or the `trace` level. Here again, you can decide the format you would like to use. You can look at the specificities of every format in the documentation. Then, local waveforms can be loaded directly in with the `obspy.read` method. This method allows to read almost all existing seismic data formats. Once the previous cell was ran, you can execute the next cell to read back the data.\n",
    "\n",
    "> __Question 4.__ Use the code below to save the seismogram locally. Change the format between SAC and MSEED. What main differences exist between these two formats? Load the locally-saved stream from MSEED format. What missing in the trace metadata after reading?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write\n",
    "stream.write(filename=\"seismogram.sac\")\n",
    "\n",
    "# Read the saved stream below into a new stream and explore its metadata\n",
    "# (compare with the original stream)\n",
    "stream_sac = obspy.read(\"seismogram01.sac\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.5. Instrument response correction\n",
    "\n",
    "For a physical interpretation of the waveform amplitude, we need to convert the raw waveforms into physical units such as ground displacement, velocity, or acceleration. Because we asked the program to _attach the response_ in the `get_waveform` method, we can use the `trace.remove_response()` method to remove the instrument response. Please check the documentation to understand how this method works.\n",
    "\n",
    "> __Question 5.__ Why do we need for such a complete tool to remove the instrument response, and why don't we simply infer how many velocity is 1 count, and divide the counts by this value? Add code to the cell below to show the corrected waveform. What changed compared with the raw waveform? We can use a `pre_filt` tuple that defines a bandpass filter in frequency domain to the data before deconvolution. The list or tuple defines the four corner frequencies (f1, f2, f3, f4) of a cosine taper which is one between f2 and f3 and tapers to zero for f1 < f < f2 and f3 < f < f4. Can you use this filter to remove the very-low-frequency signature? Plot the displacement $u[n]$, velocity $\\dot u[n]$, and acceleration $\\ddot u[n]$ on three separate subplots. What are the main differences between the signals?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prefilter, remove response and get physical units\n",
    "filter_frequencies = 0.01, 0.015, 15, 20\n",
    "\n",
    "# Loop over the units\n",
    "units = {\n",
    "    \"ACC\": \"$\\ddot u(t)$ in $m/s^2$\",\n",
    "    \"VEL\": \"$\\dot u(t)$ in $m/s$\",\n",
    "    \"DISP\": \"$u(t)$ in $m$\",\n",
    "}\n",
    "for index, (unit_name, unit_desc) in enumerate(units.items()):\n",
    "    # Remove response\n",
    "    trace_physical = trace.copy()\n",
    "    trace_physical.remove_response(\n",
    "        output=unit_name, pre_filt=filter_frequencies\n",
    "    )\n",
    "\n",
    "    # Perform the plot below\n",
    "    times = trace_physical.times()\n",
    "    waveform_physical = trace_physical.data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Spectral analysis\n",
    "\n",
    "### 3.1. Fourier transform\n",
    "\n",
    "The time-domain representation of a waveform gives information about the temporal evolution of the rupture, the different kinds of seismic phases recorded, etc. Another critical tool of the seismologist is to observe the waveform's properties as a function of frequency. The Fourier transform is the most robust and stable tool to calculate the waveform's spectrum. Note that many other transforms allow visualizing the spectral properties of signals against frequencies. Other methods exist because of the robustness to noise, capability to be applied on very-short duration time series, or other specificities.\n",
    "\n",
    "> __Question 6.__ Use the `numpy.fft` module to calculate the Fourier spectrum $\\hat u(f)$ of the seismic trace $u[n]$ in the cell below. Withing several subplots, show the raw, displacement, velocity, and acceleration spectrum $|\\hat u(f)|$  in a log-log scale. What is the effect of the pre-filter in every case? Give a very simple interpretation of the observed spectra: where is the energy located for each physical unit?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2. Filtering\n",
    "\n",
    "Thanks to digital processing tools, it is possible to filter out the waveform content within specific frequency bands. The cell below shows how to filter the seismic signal in a band of interest.\n",
    "\n",
    "> __Question 7.__  Show the filtered trace in both time and frequency domains. What is the effect of the filter? Using a low-pass filter below 0.1 Hz, what features of the waveform do you better highlight? And with a high-pass filter above 2 Hz? By playing with the `corner` value in the filter above, how do you better select the frequencies of interest?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABp2klEQVR4nO3deXhM5/s/8Pdkm0lIBtlDRKg91uSLRAlKbKW0ilJbtRVqDW2p1t6m2lKlEh97tZQPqh+tFKmtIbGF2KIoiQSJSJDEkv35/eGXqTGTmImZTObk/bquuS7zzHPOuZ8k53bPOc85RyaEECAiIiIis2dh6gCIiIiIyDBY2BERERFJBAs7IiIiIolgYUdEREQkESzsiIiIiCSChR0RERGRRLCwIyIiIpIIFnZEREREEsHCjoiIiEgiWNjRCzl69CjefPNNuLu7w8bGBm5ubhgwYABiYmL0Ws+cOXMgk8nKFMPBgwchk8lw8ODBMi2vq06dOqFTp07P7RcWFob169cbNRaiymz9+vWQyWRaX9OmTUNiYiJkMpnafli8TGJioqpt06ZNWLJkSbnHrwttY9AmPj4ec+bMURsXVW5Wpg6AzNeyZcswefJktGnTBl999RW8vLyQlJSE5cuX4+WXX8Z3332H8ePH67Sud999Fz169ChTHK1bt0ZMTAyaNGlSpuUNLSwsDE5OThg5cqSpQyGStHXr1qFRo0ZqbR4eHnB1dUVMTAzq1atX6vKbNm3C+fPnMXnyZCNGaVzx8fGYO3cuOnXqhDp16pg6HKoAWNhRmRw5cgSTJ09Gr169sGPHDlhZ/funNHjwYPTv3x+TJk1Cq1at0L59+xLX8+jRI9jZ2aFWrVqoVatWmWJxcHBAu3btyrSsqeXn50Mmk6n9/IhINz4+PvDz89P6mSlzwuPHj2Fra2uy7ZemOOeSdPFULJVJaGgoZDIZwsPDNYoSKysrhIWFQSaT4csvv1S1F59uPXXqFAYMGIDq1aurvlFrOxWbm5uLqVOnws3NDXZ2dujYsSNiY2NRp04dtaNh2k7Fjhw5ElWrVsU///yDXr16oWrVqvD09MTUqVORm5urtp25c+eibdu2qFGjBhwcHNC6dWusWbMGQgi9fy516tTBhQsXcOjQIdWpoeJv0cVx/vjjj5g6dSpq1qwJuVyOf/75B3fu3MG4cePQpEkTVK1aFS4uLujSpQuioqI0tpGbm4t58+ahcePGUCgUcHR0ROfOnREdHa3qI4RAWFgYWrZsCVtbW1SvXh0DBgzAtWvX9B4TkbnR5TRmp06dsGvXLly/fl3tVG6xvLw8LFiwAI0aNYJcLoezszNGjRqFO3fuqK2nTp06ePXVV/HLL7+gVatWUCgUmDt3LgAgNTUVY8aMQa1atWBjYwNvb2/MnTsXBQUFauu4desWBg4cCHt7eyiVSgwaNAipqanPHef69evx5ptvAgA6d+6sGkPxuDt16gQfHx/89ddfCAgIgJ2dHd555x0AwJYtWxAUFAR3d3fY2tqicePGmD59Oh4+fKixnWPHjqFPnz5wdHSEQqFAvXr1NI5yXrlyBUOGDIGLiwvkcjkaN26M5cuXP3cMZHg8TEB6KywsxIEDB+Dn51fiUTZPT0/4+vpi//79KCwshKWlpeqz119/HYMHD0ZwcLDWJFJs1KhR2LJlCz766CN06dIF8fHx6N+/P7KysnSKMz8/H3379sXo0aMxdepU/PXXX5g/fz6USiVmzZql6peYmIgxY8agdu3aAJ7MG5wwYQJu3ryp1k8XO3bswIABA6BUKhEWFgYAkMvlan1mzJgBf39/rFixAhYWFnBxcVH9ZzF79my4ubnhwYMH2LFjBzp16oR9+/ap5vYVFBSgZ8+eiIqKwuTJk9GlSxcUFBTg6NGjSEpKQkBAAABgzJgxWL9+PSZOnIiFCxfi7t27mDdvHgICAnDmzBm4urrqNS6iiqiwsFCjSNL16HdYWBjef/99XL16FTt27FD7rKioCK+99hqioqLw0UcfISAgANevX8fs2bPRqVMnnDx5Uu2I3KlTp3Dx4kV8+umn8Pb2RpUqVZCamoo2bdrAwsICs2bNQr169RATE4MFCxYgMTER69atA/Dk6F7Xrl1x69YthIaGokGDBti1axcGDRr03DH07t0bX3zxBT755BMsX74crVu3BgC1U9ApKSl4++238dFHH+GLL76AhcWT4zlXrlxBr169MHnyZFSpUgV///03Fi5ciOPHj2P//v2q5ffs2YM+ffqgcePGWLx4MWrXro3ExETs3btX1Sc+Ph4BAQGoXbs2Fi1aBDc3N+zZswcTJ05Eeno6Zs+erdPvhAxEEOkpNTVVABCDBw8utd+gQYMEAHH79m0hhBCzZ88WAMSsWbM0+hZ/VuzChQsCgPj444/V+v38888CgBgxYoSq7cCBAwKAOHDggKptxIgRAoD473//q7Z8r169RMOGDUuMubCwUOTn54t58+YJR0dHUVRUpPosMDBQBAYGljpmIYRo2rSp1n7FcXbs2PG56ygoKBD5+fnilVdeEf3791e1b9iwQQAQq1atKnHZmJgYAUAsWrRIrT05OVnY2tqKjz766LnbJ6rI1q1bJwBofeXn54uEhAQBQKxbt05jmYSEBFVb7969hZeXl8b6i/PM9u3b1dpPnDghAIiwsDBVm5eXl7C0tBSXLl1S6ztmzBhRtWpVcf36dbX2b775RgAQFy5cEEIIER4eLgCI//3vf2r93nvvPY0xaLN161aN/FcsMDBQABD79u0rdR1FRUUiPz9fHDp0SAAQZ86cUX1Wr149Ua9ePfH48eMSl+/evbuoVauWyMzMVGsfP368UCgU4u7du6VunwyLp2L18Ndff6FPnz7w8PCATCbDr7/+avRt3rx5E2+//TYcHR1hZ2eHli1bIjY21ujbNQTx/09lPnuK9Y033njusocOHQIADBw4UK19wIABOn8jl8lk6NOnj1pb8+bNcf36dbW2/fv3o2vXrlAqlbC0tIS1tTVmzZqFjIwMpKWl6bQtfZQ0/hUrVqB169ZQKBSwsrKCtbU19u3bh4sXL6r6/PHHH1AoFKrTKdr8/vvvkMlkePvtt1FQUKB6ubm5oUWLFka/epiovGzYsAEnTpxQexlivurvv/+OatWqoU+fPmr7UMuWLeHm5qaxDzVv3hwNGjTQWEfnzp3h4eGhto6ePXsC+DfHHThwAPb29ujbt6/a8kOGDHnhcQBA9erV0aVLF432a9euYciQIXBzc1PlvcDAQABQ5ZzLly/j6tWrGD16NBQKhdb15+TkYN++fejfvz/s7OzUxtqrVy/k5OTg6NGjBhkL6YanYvXw8OFDtGjRAqNGjdKpOHlR9+7dQ/v27dG5c2f88ccfcHFxwdWrV1GtWjWjb7s0Tk5OsLOzQ0JCQqn9EhMTYWdnhxo1aqi1u7u7P3cbGRkZAKBxytDKygqOjo46xWlnZ6eRjORyOXJyclTvjx8/jqCgIHTq1AmrVq1SzYX59ddf8fnnn+Px48c6bUsf2sa/ePFiTJ06FcHBwZg/fz6cnJxgaWmJzz77TK2wu3PnDjw8PFSnU7S5ffs2hBAlnm6tW7fuiw+CqAJo3LhxiRdPvIjbt2/j/v37sLGx0fp5enq62ntt+/Tt27fx22+/wdrautR1ZGRkaN1X3dzc9A1bK22xPXjwAB06dIBCocCCBQvQoEED2NnZITk5Ga+//roq7xVPESntwraMjAwUFBRg2bJlWLZsmdY+z/68yLhY2OmhZ8+eqm9b2uTl5eHTTz/Fxo0bcf/+ffj4+GDhwoU63ftMm4ULF8LT01M1FwNAhbic3dLSEp07d8bu3btx48YNrTv9jRs3EBsbi549e6rNrwM0j+BpU1y83b59GzVr1lS1FxQUqIo+Q9i8eTOsra3x+++/qxWBxjwaq238P/30Ezp16oTw8HC19uzsbLX3zs7OOHz4MIqKikos7pycnCCTyRAVFaUxvw/QnPNHROqcnJzg6OiI3bt3a/3c3t5e7b22fdrJyQnNmzfH559/rnUdHh4eAJ7kuuPHj2t8rsvFE7rQFtv+/ftx69YtHDx4UHWUDgDu37+v1s/Z2RnAk3xekurVq8PS0hLDhg3DBx98oLWPt7d3GSKnsuKpWAMaNWoUjhw5gs2bN+Ps2bN488030aNHD1y5cqVM69u5cyf8/Pzw5ptvwsXFBa1atcKqVasMHHXZzJgxA0IIjBs3DoWFhWqfFRYWYuzYsRBCYMaMGWVaf8eOHQE8uXLradu2bdOYLP0iim818nTx+fjxY/z4449lXqdcLtf7SJ9MJtMouM6ePatxo+eePXsiJyen1Kv9Xn31VQghcPPmTfj5+Wm8mjVrpldsRFJV0r766quvIiMjA4WFhVr3oYYNGz533a+++irOnz+PevXqaV1HcWHXuXNnZGdnY+fOnWrLb9q0SecxANAr5xQXe8/mnP/85z9q7xs0aIB69eph7dq1GncTKGZnZ4fOnTvj9OnTaN68udax6nqWhQyDR+wM5OrVq/j5559x48YN1Q47bdo07N69G+vWrcMXX3yh9zqvXbuG8PBwhISE4JNPPsHx48cxceJEyOVyDB8+3NBD0Ev79u2xZMkSTJ48GS+//DLGjx+P2rVrq25QfOzYMSxZskR1laa+mjZtirfeeguLFi2CpaUlunTpggsXLmDRokVQKpWlnorUR+/evbF48WIMGTIE77//PjIyMvDNN9+80FGtZs2aYfPmzdiyZQvq1q0LhULx3GLq1Vdfxfz58zF79mwEBgbi0qVLmDdvHry9vdUK2bfeegvr1q1DcHAwLl26hM6dO6OoqAjHjh1D48aNMXjwYLRv3x7vv/8+Ro0ahZMnT6Jjx46oUqUKUlJScPjwYTRr1gxjx44t8/iIpKJZs2b45ZdfEB4eDl9fX1hYWMDPzw+DBw/Gxo0b0atXL0yaNAlt2rSBtbU1bty4gQMHDuC1115D//79S133vHnzEBkZiYCAAEycOBENGzZETk4OEhMTERERgRUrVqBWrVoYPnw4vv32WwwfPhyff/456tevj4iICOzZs0enMfj4+AAAVq5cCXt7eygUCnh7e5daTAUEBKB69eoIDg7G7NmzYW1tjY0bN+LMmTMafZcvX44+ffqgXbt2mDJliirP79mzBxs3bgQAfPfdd3j55ZfRoUMHjB07FnXq1EF2djb++ecf/Pbbb2pX2VI5MOmlG2YMgNixY4fq/X//+18BQFSpUkXtZWVlJQYOHCiEEKortUp7ffDBB6p1WltbC39/f7XtTpgwQbRr165cxqiLmJgYMWDAAOHq6iqsrKyEi4uLeP3110V0dLRG3+IrX+/cuVPiZ0/LyckRISEhwsXFRSgUCtGuXTsRExMjlEqlmDJliqpfSVfFVqlSRaftrF27VjRs2FDI5XJRt25dERoaKtasWaNxBZ2uV8UmJiaKoKAgYW9vLwCorrorjnPr1q0ay+Tm5opp06aJmjVrCoVCIVq3bi1+/fVXMWLECI2r9h4/fixmzZol6tevL2xsbISjo6Po0qWLxs987dq1om3btqJKlSrC1tZW1KtXTwwfPlycPHnyuWMgqsiKr3A9ceKE1s91vSr27t27YsCAAaJatWpCJpOp5Yb8/HzxzTffiBYtWgiFQiGqVq0qGjVqJMaMGSOuXLmi6ufl5SV69+6tNY47d+6IiRMnCm9vb2FtbS1q1KghfH19xcyZM8WDBw9U/W7cuCHeeOMNUbVqVWFvby/eeOMNER0drdNVsUIIsWTJEuHt7S0sLS3VlgkMDBRNmzbVukx0dLTw9/cXdnZ2wtnZWbz77rvi1KlTWrcZExMjevbsKZRKpZDL5aJevXpqOViIJz/zd955R9SsWVNYW1sLZ2dnERAQIBYsWPDc+MmwZEKU4S6sBJlMhh07dqBfv34AnpwyHDp0KC5cuKAxp6xq1apwc3NDfn4+rl69Wup6q1evrppI6+XlhW7dumH16tWqz8PDw7FgwQLcvHnTsAMyE9HR0Wjfvj02btxosKvGiIiIpIKnYg2kVatWKCwsRFpaGjp06KC1j7W1tcZzDUvTvn17XLp0Sa3t8uXL8PLyeqFYzUVkZCRiYmLg6+sLW1tbnDlzBl9++SXq16+P119/3dThERERVTgs7PTw4MED/PPPP6r3CQkJiIuLQ40aNdCgQQMMHToUw4cPx6JFi9CqVSukp6dj//79aNasGXr16qX39qZMmYKAgAB88cUXGDhwII4fP46VK1di5cqVhhxWheXg4IC9e/diyZIlyM7OhpOTE3r27InQ0NAS76lERERUmfFUrB4OHjyIzp07a7SPGDEC69evR35+PhYsWIANGzbg5s2bcHR0hL+/P+bOnVvmKxF///13zJgxA1euXIG3tzdCQkLw3nvvvehQiIiISIJY2BERERFJBO9jR0RERCQRLOyIiIiIJIIXTzxHUVERbt26BXt7e50ehUVE5kUIgezs7Oc+g7e8MfcQSZuxcg8Lu+e4desWPD09TR0GERlZcnJyqQ87L2/MPUSVg6FzDwu75yh+2HNycjIcHBxMHA0RGVpWVhY8PT01Huxuasw9RNJmrNzDwu45ik+BODg4MLkSSVhFO93J3ENUORg691ScCSVERERE9EJY2BERERFJBAs7IiIiIolgYUdEREQkESzsTOzW/cfoseQvbDmRZOpQiKgSWX7gH/T9/jAe5BaYOhQiMiAWdia2YFc8/k7Nxsfbz5k6FCKqRL7ecwlnb2RiQ0yiqUMhIgNiYWdij/IKTR0CEVVi+QXC1CEQkQGxsCMiIiKSCBZ2Jib4ZZmIiIgMhIUdEVElVsEeuEFEL4iFHREREZFEsLAzsSKeiyUiE2IKIpIWFnYmFnUl3dQhEFEl9mvcTVOHQEQGxMKOiKgSS0h/aOoQiMiAWNgRERERSQQLOyKiSib9Qa6pQyAiI2FhR0RUyeQVFJk6BCIyEhZ2RERERBJhpe8CiYmJiIqKQmJiIh49egRnZ2e0atUK/v7+UCgUxoiRiIiIiHSg8xG7TZs2oV27dqhbty4+/PBD/Prrr4iKisLq1avRo0cPuLq6Yty4cbh+/brRgv3rr7/Qp08feHh4QCaT4ddff33uMocOHYKvry8UCgXq1q2LFStWGC0+IiIiIlPSqbBr3bo1Fi9ejLfffhuJiYlITU1FbGwsDh8+jPj4eGRlZeF///sfioqK4Ofnh61btxol2IcPH6JFixb4/vvvdeqfkJCAXr16oUOHDjh9+jQ++eQTTJw4Edu3bzdKfERERESmpNOp2Pnz56N3794lfi6Xy9GpUyd06tQJCxYsQEJCgsECfFrPnj3Rs2dPnfuvWLECtWvXxpIlSwAAjRs3xsmTJ/HNN9/gjTfeMEqMREQVHR82QSRdOhV2pRV1z3JycoKTk1OZAzKkmJgYBAUFqbV1794da9asQX5+PqytrTWWyc3NRW7uv7cCyMrKMnqcRETMPURkCHpfFXvq1CmcO3dO9f5///sf+vXrh08++QR5eXkGDe5FpaamwtXVVa3N1dUVBQUFSE/X/iiv0NBQKJVK1cvT07M8QiWiSo65h4gMQe/CbsyYMbh8+TIA4Nq1axg8eDDs7OywdetWfPTRRwYP8EXJZDK19+L/P/H62fZiM2bMQGZmpuqVnJxs9BiJiJh7iMgQ9L7dyeXLl9GyZUsAwNatW9GxY0ds2rQJR44cweDBg1Xz2SoCNzc3pKamqrWlpaXBysoKjo6OWpeRy+WQy+XlER4RkQpzDxEZgt5H7IQQKCp6ctfyP//8E7169QIAeHp6lnh601T8/f0RGRmp1rZ37174+flpnV9HREREZM70Luz8/PywYMEC/Pjjjzh06JDqwoqEhASN+WyG9uDBA8TFxSEuLk61zbi4OCQlJQF4cipj+PDhqv7BwcG4fv06QkJCcPHiRaxduxZr1qzBtGnTjBpnWd3J5vMbicj4iqekEJH06F3Yffvttzh16hTGjx+PmTNn4qWXXgIAbNu2DQEBAQYP8GknT55Eq1at0KpVKwBASEgIWrVqhVmzZgEAUlJSVEUeAHh7eyMiIgIHDx5Ey5YtMX/+fCxdurTC3uok83HFuviEiIiIzIvec+xatGihdlVssa+//hpWVnqvTi+dOnUq9Zvm+vXrNdoCAwNx6tQpI0ZlOPwSTURERC9C7yN2devWRUZGhkZ7Tk4OGjRoYJCgiIiIiEh/ehd2iYmJKCws1GjPzc3FjRs3DBJUZcUDdkRERPQidD53unPnTtW/9+zZA6VSqXpfWFiIffv2wdvb27DRVTI8FUtEREQvQufCrl+/fgCe3Nh3xIgRap9ZW1ujTp06WLRokUGDIyIiw+OXSCLp0rmwK753nbe3N06cOFFhngdLRERERE/ofRlrQkKCMeIgIiIiohdUpvuT7Nu3D/v27UNaWprqSF6xtWvXGiSwykjw8gkiIiJ6AXoXdnPnzsW8efPg5+cHd3d3yGQyY8RVKT1TIxMRERHpRe/CbsWKFVi/fj2GDRtmjHgqtbM37qOJh4OpwyAiIiIzpfd97PLy8oz+6LDKqqCIp2KJiIio7PQu7N59911s2rTJGLEQERER0QvQ+1RsTk4OVq5ciT///BPNmzeHtbW12ueLFy82WHBEREREpDu9C7uzZ8+iZcuWAIDz58+rfcYLKV5MS89qpg6BiIiIzJjehd2BAweMEQcBsGBhTERERC9A7zl2RERk3vhIMSLp0vuIXefOnUs95bp///4XCqgy4wE7IiIiehF6F3bF8+uK5efnIy4uDufPn8eIESMMFRcRERER6Unvwu7bb7/V2j5nzhw8ePDghQOqzNKyc9HY3dRREBERkbky2By7t99+m8+JfUHf779i6hCIiIjIjBmssIuJiYFCoTDU6iolPniCiIiIXoTep2Jff/11tfdCCKSkpODkyZP47LPPDBZYScLCwvD1118jJSUFTZs2xZIlS9ChQwetfQ8ePIjOnTtrtF+8eBGNGjUydqh647UTRFQeBPgtkkiq9C7slEql2nsLCws0bNgQ8+bNQ1BQkMEC02bLli2YPHkywsLC0L59e/znP/9Bz549ER8fj9q1a5e43KVLl+Dg4KB67+zsbNQ4y+rqHc5RJCIiorLTu7Bbt26dMeLQyeLFizF69Gi8++67AIAlS5Zgz549CA8PR2hoaInLubi4oFq1auUUZdnde5Rv6hCIiIjIjJV5jl1sbCx++uknbNy4EadPnzZkTFrl5eUhNjZW46hgUFAQoqOjS122VatWcHd3xyuvvMInZxAREZFk6X3ELi0tDYMHD8bBgwdRrVo1CCGQmZmJzp07Y/PmzUY7zZmeno7CwkK4urqqtbu6uiI1NVXrMu7u7li5ciV8fX2Rm5uLH3/8Ea+88goOHjyIjh07al0mNzcXubm5qvdZWVmGGwQRUQmYe4jIEPQ+YjdhwgRkZWXhwoULuHv3Lu7du4fz588jKysLEydONEaMap596oUQosQnYTRs2BDvvfceWrduDX9/f4SFhaF379745ptvSlx/aGgolEql6uXp6WnQ+ImItDFl7sl4kPv8TkRkFvQu7Hbv3o3w8HA0btxY1dakSRMsX74cf/zxh0GDe5qTkxMsLS01js6lpaVpHMUrTbt27XDlSsn3i5sxYwYyMzNVr+Tk5DLHTESkq/LMPc8+Kzb53mOjbYuIypfehV1RURGsra012q2trVFUVGSQoLSxsbGBr68vIiMj1dojIyMREBCg83pOnz4Nd/eSH+8gl8vh4OCg9iIiMjZT5p6CQuPlbiIqX3oXdl26dMGkSZNw69YtVdvNmzcxZcoUvPLKKwYN7lkhISFYvXo11q5di4sXL2LKlClISkpCcHAwgCffeIcPH67qv2TJEvz666+4cuUKLly4gBkzZmD79u0YP368UeN8EZHxt00dAhFVMgNWxJg6BCIyEL0vnvj+++/x2muvoU6dOvD09IRMJkNSUhKaNWuGn376yRgxqgwaNAgZGRmYN28eUlJS4OPjg4iICHh5eQEAUlJSkJSUpOqfl5eHadOm4ebNm7C1tUXTpk2xa9cu9OrVy6hxvoj3NpxE4pe9TR0GERERmSGZEM/OttBNZGQk/v77bwgh0KRJE3Tt2tXQsVUIWVlZUCqVyMzMNMqpkTrTd2m0sbAjKj/G3sfLyphxXb6djaBv/1JrY94hKl/G2sf1PmJXrFu3bujWrZvBAiEiovJx4O80U4dAREai9xy7iRMnYunSpRrt33//PSZPnmyImIiIyIj4pFgi6dK7sNu+fTvat2+v0R4QEIBt27YZJKjK7nFeIZb8eRnxt3iDUiIqP0evZSD84FWUcYYOEVUAehd2GRkZUCqVGu0ODg5IT083SFCV3ZJ9l7HkzyvotTTK1KEQkQRpv6U7MHjlUSzc/Tcizml/mg8RVXx6F3YvvfQSdu/erdH+xx9/oG7dugYJqrI7m5xp6hCISMIe5hWW+nnyvUflFAkRGZreF0+EhIRg/PjxuHPnDrp06QIA2LdvHxYtWoQlS5YYOr5KqYinQYjIiLbH3ij1c4uSDukRUYWnd2H3zjvvIDc3F59//jnmz58PAKhTpw7Cw8PVbg5MZce6joiMqbBIM8lkPspX/duihOdvE1HFV6bbnYwdOxZjx47FnTt3YGtri6pVqxo6rkqNR+yIyJgKteSYnWduqv7Nwo7IfJX5PnYA4OzsbKg46Cnaki4RkaFou+o1K6dA9W+eiiUyXzpdPNGjRw9ER0c/t192djYWLlyI5cuXv3BgldnppPumDoGIJEzbqdiv91xS/duClR2R2dLpiN2bb76JgQMHwt7eHn379oWfnx88PDygUChw7949xMfH4/Dhw4iIiMCrr76Kr7/+2thxExFRGd17aj6dNjKeiiUyWzoVdqNHj8awYcOwbds2bNmyBatWrcL9+/cBPEkATZo0Qffu3REbG4uGDRsaM14iIjKyx3kFz+9ERBWSznPsbGxsMGTIEAwZMgQAkJmZicePH8PR0RHW1tZGC5CIiMrXV7sv4f2O9UwdBhGVQZkvnlAqlVqfQEFEROatQMscPCIyD3o/eYLKVxETLBEREemIhV0Fd/L6PVOHQERERGaChV0Fl3zX+M9sLCoSyM4p/So5IiJDy8rJ13pPPSIqOxZ2FdzUrWeMvo0R646j2Zy9uHbngdG3RUTm4eClNKOu/2JKFprP2Yv3f4w16naIKpsyFXb379/H6tWrMWPGDNy9excAcOrUKdy8efM5S1JZpGbmGHX9UVfSAQBbn/NgcCKqPEauO2HU9f8QnQgAiIy/bdTtEFU2ehd2Z8+eRYMGDbBw4UJ88803qvvZ7dixAzNmzDB0fAQg8uJtRP/zpPjKLSjEicS7KCgsMnFURP/KfJSPDTGJyHiQa+pQyID+F3dTNU3j2p0H5TI1hEgfR/5Jx5/8cqBG79udhISEYOTIkfjqq69gb2+vau/Zs6fqHneV1eO8QsSnZKKVZ3WDPpLns1/PAwDeblcbG48loXhKSuKXvQ22DSqZEALD1x6HlYUMa0f+H+/Kr8X7P57EsYS72Hw8GRGTOpg6nErnyu1sONhaw9VBYdD1TtocBwBYOcxXdcr0wtzuqCJ/oceMk44OX0nHh9vO4Ms3miOwAZ/N/iwhBIauPgYAOPlpVzhVlZs4oopB7yN2J06cwJgxYzTaa9asidTUVIMEVZqwsDB4e3tDoVDA19cXUVFRpfY/dOgQfH19oVAoULduXaxYscJosY1afxxvhMdg7ZEEo6z/p6P/FnUA8GNMIs4k33/ucsl3H2Ho6qM48Jw5M5V1DnNRkUDM1YwSLyC5kvYAUVfSceDSHcSnZKnak+8+wt2HeeUVZoV2LOHJlIynfz5UPlIyH6Pbt3+h7Rf7jLaNp+fBNZ29B7HX7yLz8fMvuFr511UE/xjLMwwluJ2Vg7hScvjba44hJTMHI9YeV7UVFBbhwq1MXnQC4GFeoerfd7J5tqCY3oWdQqFAVpZm8r506RKcnY37jWLLli2YPHkyZs6cidOnT6NDhw7o2bMnkpKStPZPSEhAr1690KFDB5w+fRqffPIJJk6ciO3btxslvqPXnvzntum49ngM7bP/XcBry4+gzvRd2HU2BVk5+bhwK1P1uRACf12+gwk/n8aRfzIw6jlzZlYcumrskCukH2IS8daqo2g2Z6/GfQMf5BYg6Nu/VO97Lz2MhPSHyHiQiw5fHUDr+ZG81yCZ1EUTFNNvhMegxdy96B92RFVoPF3o3c7Kwfmbmfgi4m/svpCKA5fuaKwjt+DfYq+yFiltv9iHfsuP4Js9lzQ+23NB/UBJv+VHUFgkMHHzafReelh1NLUye5Dz76PvVv11zYSRVCwyoece9f777+POnTv473//ixo1auDs2bOwtLREv3790LFjRyxZssRIoQJt27ZF69atER4ermpr3Lgx+vXrh9DQUI3+H3/8MXbu3ImLFy+q2oKDg3HmzBnExMTotM2srCwolUpkZmbCwcGh1L51pu8CALzkUhV/hgTqtP7iZUzBXm6FtaP+D2+u+Pdn0djdAbsmvIzrdx9BaWsNB4UV9sbfxriNp9C/VU18NaA5ioTAgb/TENTETa9TzkVFAuduZqKgSKBmNVucvXEf3k5V8JJLVdXpzeI/R11Pd+YXFsFSJtMpjqSMR1DaWWPuzgsY1d4bzWo9eXJKp68PIDHj37lDU7o2wLd/XtZ5XE+L/bQr1h5JwPIDV3Hy066ws7GErbVluZ6+zS0ohI2lxXO3mZadg3M3MvFKY1e19oe5BcjOKYCbsuTTegnpD1Gjig2Utk8eJ/js3/GELi+hf6uacFfawtbGEoVFAo/zC2FlIYPC2rLUuFYcuoov//gbZ2YFIa+wCFduZyPgJadSl3kR+uzj5UnXuPZdvI3RP5wEoPv0DFPmHQAY16kewg6qf5GMnt4FNarY4Nqdh6jvWhVWFjJ4z4gAAByY1gl1HO1wKuk+nKvKUdvRTq/tpWXnICnjEewV1rCylOH8zUz0auYOKwtZmXNPbkEh5Fal/y0X97t1PwfJdx9h/99pmNm7MawtnxxTefr3MKRtbWw6VraDAvNfa4rezT3Qen4kggPr4aPuDZFbUARbm+fHZ0g5+YXP3b+FEDh0+Q4auTlo5JjrGQ9Rs5otrCy1H3PKLSjEtTsP0cjNHjKZDBkPcuG74E/V51VsLPHVgBbo2MAJttaWsLK0wMPcAshkgJ1N6dMH7j/KQ8t5kZj4Sn1M6Vofh/9Jh4+HEtWr2Og4ev0ZK/foXdhlZWWhV69euHDhArKzs+Hh4YHU1FT4+/sjIiICVapUMVhwT8vLy4OdnR22bt2K/v37q9onTZqEuLg4HDp0SGOZjh07olWrVvjuu+9UbTt27MDAgQPx6NEjnZ5xW5bCjohMr3+rmpj7WlM4KErfz829sPsz/jbe3XCyHCMjotL8Z5gvujd1e24/Y+UevWfAOjg44PDhw9i/fz9OnTqFoqIitG7dGl27djVYUNqkp6ejsLAQrq7qRxdcXV1LnNuXmpqqtX9BQQHS09Ph7u6usUxubi5yc/89V6/ttDMRVXw7Tt/E1KAGzy3sKoqy5p7H+YXP70RE5Wb5gX90KuyMpcyXNnXp0gVdunQxZCw6efYwuRCi1EPn2vpray8WGhqKuXPnvmCURGRqSltr2D7ntFBFUtbcY2PF+8wTVSRvtalt0u3rVNgtXbpU5xVOnDixzMGUxsnJCZaWlhpH59LS0jSOyhVzc3PT2t/KygqOjo5al5kxYwZCQkJU77OysuDp6alXrC72chyfqXkEs6hIQODJvLDieQimPn0777WmmPW/C2ptERM74HF+IeRWFqhmZ40tJ5KxbP8/AIBDH3aCs70c8bey0MTDQWPeQlGRgIWFDEIIXM94BDelAjaWFrCwkOFRXgG2x95A05pKbIu9gU3HkvBh94bo28IDj/IK8ZJLVVhayFBQWKR1joUQAkUCsPz/65fJZMh8nA+FtQWsLZ5sI7+wSDV/5dmiv/iq19eWH8G4Ti/hlUYusLSUofmcvWrbsbaUIb+wbJO5D3/cGb2XHkbm43wse6sValSxQava1WBlYQGZDCgsEs+dg6IPbfNHsnLyUcXGCpYW/84d0vZFZtHeS4i/lYVVw/0AQDVPMbegENczHqGBq73GMgWFRbCQybBk3xX0a+mBOo5VIJNBNRfqaXP6NMHL9Z3wkos9bt1/jIspWfDzqgGl3b9H0Ir/Xp5WPMfuyPQuyM0vxO9nUxAcWA82VhalfpErLBKqMZujsuaeGk/NAUoI7aX155NfWAQrCxnyCosgt7I0ed5p6GqPS7ez1dreblcbId0a4nhCBpq4P5n/2vHrAwCAN31r4aMejfAorwAFRQJ1napojLOoSEAmAx7lFeLuwzzUrGarykXnbmbi/M0s1HG0w8h1J1AkBDa91w4OtlZwd7CF0s5alZ+1/Q0V56Tiv7/CIoHMx/moZmut2kZBkYC1pebfqBAC2bkF+DP+Npbuu4L/jX8Zj/IK4GKvQL1PNPebshju74U+LTxU86V3TXwZMshQ17kKbCwt8Di/ELbWlga9Dded7Fw42/97exEhBLIeF8DB1goymazEgygPcwsQ+PUBLB3cCu3qOqrF9E9aNmpWs9OYG1h8gdrF1CxEXUnH+x3qAgDSsnPRLlT9avDuTV0xMsAbLT2rwdpShiNXM+BcVY5GbvaqbWnLO+kPcuG34E+0f8kRP77TFj8evQ6fmg5oXbu6ajza9q2KmHd0mmPn7e2t9v7OnTt49OgRqlWrBuDJkyjs7Ozg4uKCa9eMd2VK27Zt4evri7CwMFVbkyZN8Nprr5V48cRvv/2G+Ph4VdvYsWMRFxdn1IsnalazxZHpuh3NNESCffdlb7xc3wmBDZzV/vAe5xVi/99p6NDACT/GXMdrLT1Qq7rmpOOnY6iM98Y7e+M++n5/5Mm/5wSpnborKhKo+0zyfdO3FsZ2qocui57M6zw+8xW42Bv2/mHmxpz/hsx9jl301XQMWfXkXl7lffHEnD5N0Ku5u8bf/4Vbmari5052Ll5vXavUGC7O61HuE/0rgtL2m9TMHI2iZffkDui99DAK/3+hY277mqHdfZiH1vMjAQDjO7+Ead0bmjgi/Zh0jl1Cwr/3Zdu0aRPCwsKwZs0aNGz45Id46dIlvPfee1rvb2dIISEhGDZsGPz8/ODv74+VK1ciKSkJwcHBAJ5847158yY2bNgA4MkVsN9//z1CQkLw3nvvISYmBmvWrMHPP/9s1DjLq3rv0dQN3w5qWWJCtLWxRO/mT+YRftD5peeuz8W+ct7csXmtaoic0hGFQmjMx7KwkGFYOy/8ePQ6AODD7g1VP8vfJ7wMpa11pS/qyLQKynh0+UXsmxqIes5VS/y8qYfyueuwtbZUzQ+sjEUdAFxe0BPRV9PR1lvzDNKzV4zGz+sOOxsrnJsThNNJ99HWu0Z5hVlhWVv++3/tW21Ne/qzItF7jt1nn32Gbdu2qYo6AGjYsCG+/fZbDBgwAEOHDjVogE8bNGgQMjIyMG/ePKSkpMDHxwcRERHw8vICAKSkpKjd087b2xsRERGYMmUKli9fDg8PDyxduhRvvPGG0WIEACvL8insVgzzNej6ejXTvJiksqiv5bRjseBO9VSF3dMFsk/N5//nVVnUda6Ca3ceonYN/W5DQS+usJzvo/hnSOlFna7e9KuFDTHXDRCR+bKxskCnhi4lfm5p8eS0b4f6TqrpFnY2VmhvxNv/mJOnbzdjXU7/75oDvQu7lJQU5Odr3nG8sLAQt28b/3lt48aNw7hx47R+tn79eo22wMBAnDp1yshRPTGhy0tYtv8fzOnT1Cjr/3t+DyisLZGYbtj/QH1qOuD8zSy8oeV0CT05tb7hnTaqe7aRph9GtcHqqGt452Xv53cmg/o/7xqwtJChvsuLF1vaDPSrha8GtMCjvAJkPS79/ob6GOD7pLDzqVlxTn9XNCdndsX+v9PQs5nprrCsyGysLBD6ejM8yivkmZOn6H0fuz59+iApKQlr1qyBr68vZDIZTp48iffeew+enp7YuXOnsWI1CX3PgT/KK3jujRCfpstcl4vzesDGysJop3jzCopw92GewRI2kTkx9zl2wJMbw1pb6p4jdJ1jd+XznqqLkYwhJfMxnKrKjboNoorKWLlH771p7dq1qFmzJtq0aQOFQgG5XI62bdvC3d0dq1evNlhg5kqfok5XtjaWRp23Z2NlwaKOyIwprA2fI+JmdTN6weWutGVRR2Rgelchzs7OiIiIwOXLl/H3339DCIHGjRujQYMGxoiv0vuPgefRERHpopqd8R6lRETGU+bDSw0aNGAxVw6acYI+ERER6Ujvwu6dd94p9fO1a9eWORjS5M5TpERERKQjvQu7e/fuqb3Pz8/H+fPncf/+fZM8YkzqSntcGhEREdHT9C7sduzYodFWVFSEcePGoW7dugYJioiIjMfOxhKP8gpNHQYRGYFBLkeysLDAlClT8O233xpidUREZEQ8D0AkXQa7zvzq1asoKCgw1OqIiMhIxunwiEEiMk96n4oNCQlRey+EQEpKCnbt2oURI0YYLDAiIjKO5z0SjI+GIzJfehd2p0+fVntvYWEBZ2dnLFq06LlXzJJuWteuhlNJ900dBhFVIhveaYPha48DAPq19DBxNERUVnoXdgcOHDBGHPSU6rwxKBGVsw71/32wvIURn3RDRMal9xy7Ll264P79+xrtWVlZvN2JgXRu5AIAsJcb/vFkRETa7qL09K2VfL2ql2M0RGRIelcOBw8eRF5enkZ7Tk4OoqKiDBJUZfdWm9pwtpejlWc1U4dCRBJU0vG4qI8645+0B+hQ37lc4yEiw9G5sDt79qzq3/Hx8UhNTVW9LywsxO7du1GzZk3DRldJWVrI0L2pm6nDICKJKunG55417ODJCyeIzJrOhV3Lli0hk8kgk8m0nnK1tbXFsmXLDBocEREREelO58IuISEBQgjUrVsXx48fh7Pzv4fqbWxs4OLiAktLS6MESUREhsNLI4ikS+fCzsvLC8CTx4cREZH54iOoiaRLp8Ju586d6NmzJ6ytrbFz585S+/bt29cggRERERGRfnQq7Pr164fU1FS4uLigX79+JfaTyWQoLOSDpYmIKjIesSOSLp0Ku6dPv/JUrHE1cXcwdQhEJHElXRVLROZP7xsUm8q9e/cwbNgwKJVKKJVKDBs2TOuNkp82cuRI1ZW8xa927dqVT8BlZGvDC1CIyLhY1hFJl05H7JYuXarzCidOnFjmYEozZMgQ3LhxA7t37wYAvP/++xg2bBh+++23Upfr0aMH1q1bp3pvY8PHdREREZE06VTYffvttzqtTCaTGaWwu3jxInbv3o2jR4+ibdu2AIBVq1bB398fly5dQsOGDUtcVi6Xw82NN/slIiIi6dOpsEtISDB2HKWKiYmBUqlUFXUA0K5dOyiVSkRHR5da2B08eBAuLi6oVq0aAgMD8fnnn8PFxaXE/rm5ucjNzVW9z8rKMswgiIhKUZ65h3PsiKTrhebYCSEghDBULCUqviL3WS4uLmqPNntWz549sXHjRuzfvx+LFi3CiRMn0KVLF7Xk+azQ0FDVPD6lUglPT0+DjIGIqDTlmXua11Qabd1EZFplKuzWrFkDHx8fKBQKKBQK+Pj4YPXq1XqvZ86cORoXNzz7OnnyJADt3zCFEKV+8xw0aBB69+4NHx8f9OnTB3/88QcuX76MXbt2lbjMjBkzkJmZqXolJyfrPa4XUR6FMhFVPOWZeywseMSOSKp0fvJEsc8++wzffvstJkyYAH9/fwBPTpVOmTIFiYmJWLBggc7rGj9+PAYPHlxqnzp16uDs2bO4ffu2xmd37tyBq6urzttzd3eHl5cXrly5UmIfuVwOuVyu8zqJiAyBuYeIDEHvwi48PByrVq3CW2+9pWrr27cvmjdvjgkTJuhV2Dk5OcHJyem5/fz9/ZGZmYnjx4+jTZs2AIBjx44hMzMTAQEBOm8vIyMDycnJcHd313kZIiIiInOh96nYwsJC+Pn5abT7+vqioKDAIEE9q3HjxujRowfee+89HD16FEePHsV7772HV199Ve3CiUaNGmHHjh0AgAcPHmDatGmIiYlBYmIiDh48iD59+sDJyQn9+/c3SpxEREREpqR3Yff2228jPDxco33lypUYOnSoQYLSZuPGjWjWrBmCgoIQFBSE5s2b48cff1Trc+nSJWRmZgIALC0tce7cObz22mto0KABRowYgQYNGiAmJgb29vZGi/NFTQ0q+QpfIiIiotLofSoWeHLxxN69e1VPcTh69CiSk5MxfPhwhISEqPotXrzYMFECqFGjBn766adS+zx94YGtrS327NljsO2Xl/YvPf/UNBHRC3nmGq1twf6miYOIDE7vwu78+fNo3bo1AODq1asAAGdnZzg7O+P8+fOqfrxPEhGRefCoZmvqEIjIQPQu7A4cOGCMOAhA54bOpg6BiCohFnZE0vFCNygmw6pehc+xJSIiorLT+4hdTk4Oli1bhgMHDiAtLQ1FRUVqn586dcpgwRERERGR7vQu7N555x1ERkZiwIABaNOmDefSERGZGfHs1RNEJBl6F3a7du1CREQE2rdvb4x4KrWRAXVMHQIRERGZMb3n2NWsWbNC3wfOnDVxdzB1CERERGTG9C7sFi1ahI8//hjXr183RjxEREREVEZ6n4r18/NDTk4O6tatCzs7O1hbW6t9fvfuXYMFV9lwviIRERG9CL0Lu7feegs3b97EF198AVdXVxYjRERmRvDaCSLJ0ruwi46ORkxMDFq0aGGMeIiIiIiojPSeY9eoUSM8fvzYGLEQERER0QvQu7D78ssvMXXqVBw8eBAZGRnIyspSexERERGRaeh9KrZHjx4AgFdeeUWtXQgBmUyGwsJCw0RWCVlwuiIRERG9AL0LuwMHDhgjDgKviiWi8sFrJ4ikS+/CLjAwsMTP4uLiXiQWIiIiInoBes+xe1ZmZibCwsLQunVr+Pr6GiImIiIiIiqDMhd2+/fvx9tvvw13d3csW7YMvXr1wsmTJw0ZGxERERHpQa9TsTdu3MD69euxdu1aPHz4EAMHDkR+fj62b9+OJk2aGCtGIiIyEl60RSQtOh+x69WrF5o0aYL4+HgsW7YMt27dwrJly4wZGxERERHpQecjdnv37sXEiRMxduxY1K9f35gxERGREQk+U4xIsnQ+YhcVFYXs7Gz4+fmhbdu2+P7773Hnzh1jxqbm888/R0BAAOzs7FCtWjWdlhFCYM6cOfDw8ICtrS06deqECxcuGDdQIiIiIhPRubDz9/fHqlWrkJKSgjFjxmDz5s2oWbMmioqKEBkZiezsbGPGiby8PLz55psYO3aszst89dVXWLx4Mb7//nucOHECbm5u6Natm9FjJSIiIjIFva+KtbOzwzvvvIPDhw/j3LlzmDp1Kr788ku4uLigb9++xogRADB37lxMmTIFzZo106m/EAJLlizBzJkz8frrr8PHxwc//PADHj16hE2bNhktTiIiczK4TW1Th0BEBvRC97Fr2LAhvvrqK9y4cQM///yzoWIyiISEBKSmpiIoKEjVJpfLERgYiOjo6BKXy83N5fNviajcmSr3NPVwKJftEFH5eOEbFAOApaUl+vXrh507dxpidQaRmpoKAHB1dVVrd3V1VX2mTWhoKJRKperl6elp1DiJiIDyzT28dIJIugxS2JXVnDlzIJPJSn296E2Pn33+qhCi1GeyzpgxA5mZmapXcnLyC22fiEgXpso9MvBGdkRSovezYg1p/PjxGDx4cKl96tSpU6Z1u7m5AXhy5M7d3V3VnpaWpnEU72lyuRxyubxM2yQiKivmHiIyBJMWdk5OTnBycjLKur29veHm5obIyEi0atUKwJMraw8dOoSFCxcaZZtEROamrnMVU4dARAZk0lOx+khKSkJcXBySkpJQWFiIuLg4xMXF4cGDB6o+jRo1wo4dOwA8OQU7efJkfPHFF9ixYwfOnz+PkSNHws7ODkOGDDHVMIiIKpQGrvamDoGIDMikR+z0MWvWLPzwww+q98VH4Q4cOIBOnToBAC5duoTMzExVn48++giPHz/GuHHjcO/ePbRt2xZ79+6FvT0TGRFVXnzwBJF0mU1ht379eqxfv77UPs8+Jkcmk2HOnDmYM2eO8QIjIjJjvHSCSFrM5lQsEREREZWOhR0RERGRRLCwIyIiIpIIFnYm5urw5L5V9V2qmjgSIqos5Nb/pn5rK/43QCQlZnPxhFRted8faw4n4P2OdU0dChFVEg4Ka8zt2xQAUFXO/waIpIR7tInVcaqC+f18TB0GEVUyIwLqmDoEIjICHoMnIiIikggWdkREREQSwcKOiIiISCJY2BERERFJBC+eeI7ix5RlZWWZOBIiMobiffvZRxKaGnMPkbQZK/ewsHuO7OxsAICnp6eJIyEiY8rOzoZSqTR1GCrMPUSVg6Fzj0xUtK+pFUxRURFu3boFe3t7yGSlPy47KysLnp6eSE5OhoODQzlFWDFU1rFz3OY/biEEsrOz4eHhAQuLijM7hbnn+TjuyjVuQFpjN1bu4RG757CwsECtWrX0WsbBwcHs/+DKqrKOneM2bxXpSF0x5h7dcdyVj1TGbozcU3G+nhIRERHRC2FhR0RERCQRLOwMSC6XY/bs2ZDL5aYOpdxV1rFz3JVr3BVVZf19cNyVa9xA5R67rnjxBBEREZFE8IgdERERkUSwsCMiIiKSCBZ2RERERBLBwo6IiIhIIljYGVBYWBi8vb2hUCjg6+uLqKgoU4eks9DQUPzf//0f7O3t4eLign79+uHSpUtqfYQQmDNnDjw8PGBra4tOnTrhwoULan1yc3MxYcIEODk5oUqVKujbty9u3Lih1ufevXsYNmwYlEollEolhg0bhvv37xt7iDoJDQ2FTCbD5MmTVW1SHffNmzfx9ttvw9HREXZ2dmjZsiViY2NVn0t13FLDvGP+f4eVKe8AzD1GJ8ggNm/eLKytrcWqVatEfHy8mDRpkqhSpYq4fv26qUPTSffu3cW6devE+fPnRVxcnOjdu7eoXbu2ePDggarPl19+Kezt7cX27dvFuXPnxKBBg4S7u7vIyspS9QkODhY1a9YUkZGR4tSpU6Jz586iRYsWoqCgQNWnR48ewsfHR0RHR4vo6Gjh4+MjXn311XIdrzbHjx8XderUEc2bNxeTJk1StUtx3Hfv3hVeXl5i5MiR4tixYyIhIUH8+eef4p9//lH1keK4pYZ55wlz/jusTHlHCOae8sDCzkDatGkjgoOD1doaNWokpk+fbqKIXkxaWpoAIA4dOiSEEKKoqEi4ubmJL7/8UtUnJydHKJVKsWLFCiGEEPfv3xfW1tZi8+bNqj43b94UFhYWYvfu3UIIIeLj4wUAcfToUVWfmJgYAUD8/fff5TE0rbKzs0X9+vVFZGSkCAwMVCVYqY77448/Fi+//HKJn0t13FLDvGPef4eVLe8IwdxTHngq1gDy8vIQGxuLoKAgtfagoCBER0ebKKoXk5mZCQCoUaMGACAhIQGpqalqY5TL5QgMDFSNMTY2Fvn5+Wp9PDw84OPjo+oTExMDpVKJtm3bqvq0a9cOSqXSpD+rDz74AL1790bXrl3V2qU67p07d8LPzw9vvvkmXFxc0KpVK6xatUr1uVTHLSXMO+b/d1jZ8g7A3FMeWNgZQHp6OgoLC+Hq6qrW7urqitTUVBNFVXZCCISEhODll1+Gj48PAKjGUdoYU1NTYWNjg+rVq5fax8XFRWObLi4uJvtZbd68GadOnUJoaKjGZ1Id97Vr1xAeHo769etjz549CA4OxsSJE7FhwwYA0h23lDDvmPffYWXMOwBzT3mwMnUAUiKTydTeCyE02szB+PHjcfbsWRw+fFjjs7KM8dk+2vqb6meVnJyMSZMmYe/evVAoFCX2k9q4i4qK4Ofnhy+++AIA0KpVK1y4cAHh4eEYPny4qp/Uxi1FzDvaVeS/w8qadwDmnvLAI3YG4OTkBEtLS41vAWlpaRrfOiq6CRMmYOfOnThw4ABq1aqlandzcwOAUsfo5uaGvLw83Lt3r9Q+t2/f1tjunTt3TPKzio2NRVpaGnx9fWFlZQUrKyscOnQIS5cuhZWVlSomqY3b3d0dTZo0UWtr3LgxkpKSAEj39y0lzDvm+3dYWfMOwNxTHljYGYCNjQ18fX0RGRmp1h4ZGYmAgAATRaUfIQTGjx+PX375Bfv374e3t7fa597e3nBzc1MbY15eHg4dOqQao6+vL6ytrdX6pKSk4Pz586o+/v7+yMzMxPHjx1V9jh07hszMTJP8rF555RWcO3cOcXFxqpefnx+GDh2KuLg41K1bV5Ljbt++vcZtJS5fvgwvLy8A0v19Swnzjvn+HVbWvAMw95SL8rtOQ9qKbzuwZs0aER8fLyZPniyqVKkiEhMTTR2aTsaOHSuUSqU4ePCgSElJUb0ePXqk6vPll18KpVIpfvnlF3Hu3Dnx1ltvab0EvVatWuLPP/8Up06dEl26dNF6CXrz5s1FTEyMiImJEc2aNatQl6A/fXWaENIc9/Hjx4WVlZX4/PPPxZUrV8TGjRuFnZ2d+Omnn1R9pDhuqWHeeUIKf4eVIe8IwdxTHljYGdDy5cuFl5eXsLGxEa1bt1Zdsm8OAGh9rVu3TtWnqKhIzJ49W7i5uQm5XC46duwozp07p7aex48fi/Hjx4saNWoIW1tb8eqrr4qkpCS1PhkZGWLo0KHC3t5e2Nvbi6FDh4p79+6Vwyh182yCleq4f/vtN+Hj4yPkcrlo1KiRWLlypdrnUh231DDvSOPvsLLkHSGYe4xNJoQQpjlWSERERESGxDl2RERERBLBwo6IiIhIIljYEREREUkECzsiIiIiiWBhR0RERCQRLOyIiIiIJIKFHREREZFEsLAjIiIikggWdqRVYmIiZDIZ4uLiKsR6ymLkyJHo169fuW+32LBhw/DFF1+YbPuGsH79elSrVk2nvr///jtatWqFoqIi4wZFksW88+KYd4iFnQSNHDkSMpkMMpkMVlZWqF27NsaOHYt79+4ZfbvPJjRPT0+kpKTAx8fHaNstKYl/9913WL9+vdG2W5qzZ89i165dmDBhgkm2bwqvvvoqZDIZNm3aZOpQyASYd55g3ilfzDuaWNhJVI8ePZCSkoLExESsXr0av/32G8aNG1fucVhaWsLNzQ1WVlblvm2lUqnztz5D+/777/Hmm2/C3t7eJNs3lVGjRmHZsmWmDoNMhHmHeccUmHfUsbCTKLlcDjc3N9SqVQtBQUEYNGgQ9u7dq9Zn3bp1aNy4MRQKBRo1aoSwsLAS11dYWIjRo0fD29sbtra2aNiwIb777jvV53PmzMEPP/yA//3vf6pv7QcPHlT7VltUVIRatWphxYoVaus+deoUZDIZrl27BgDIzMzE+++/DxcXFzg4OKBLly44c+ZMibF5e3sDAFq1agWZTIZOnToB0Pwm36lTJ0yYMAGTJ09G9erV4erqipUrV+Lhw4cYNWoU7O3tUa9ePfzxxx9q64+Pj0evXr1QtWpVuLq6YtiwYUhPTy8xnqKiImzduhV9+/ZVaw8LC0P9+vWhUCjg6uqKAQMGqD4TQuCrr75C3bp1YWtrixYtWmDbtm1qy1+4cAG9e/eGg4MD7O3t0aFDB1y9elW1zXnz5qFWrVqQy+Vo2bIldu/erVq2+Pfwyy+/oHPnzrCzs0OLFi0QExOjto3169ejdu3asLOzQ//+/ZGRkaH2+ZkzZ9C5c2fY29vDwcEBvr6+OHnypOrzvn374vjx46rfJVUuzDvMO8w7FYAgyRkxYoR47bXXVO+vXr0qmjRpIlxdXVVtK1euFO7u7mL79u3i2rVrYvv27aJGjRpi/fr1QgghEhISBABx+vRpIYQQeXl5YtasWeL48ePi2rVr4qeffhJ2dnZiy5YtQgghsrOzxcCBA0WPHj1ESkqKSElJEbm5uRrrmTp1qnj55ZfV4p06darw9/cXQghRVFQk2rdvL/r06SNOnDghLl++LKZOnSocHR1FRkaG1vEeP35cABB//vmnSElJUfV79ucQGBgo7O3txfz588Xly5fF/PnzhYWFhejZs6dYuXKluHz5shg7dqxwdHQUDx8+FEIIcevWLeHk5CRmzJghLl68KE6dOiW6desmOnfuXOLP//Tp0wKASE1NVbWdOHFCWFpaik2bNonExERx6tQp8d1336k+/+STT0SjRo3E7t27xdWrV8W6deuEXC4XBw8eFEIIcePGDVGjRg3x+uuvixMnTohLly6JtWvXir///lsIIcTixYuFg4OD+Pnnn8Xff/8tPvroI2FtbS0uX76s9vts1KiR+P3338WlS5fEgAEDhJeXl8jPzxdCCHH06FEhk8lEaGiouHTpkvjuu+9EtWrVhFKpVMXZtGlT8fbbb4uLFy+Ky5cvi//+978iLi5ObfwuLi6qvyOqPJh3mHeYdyoGFnYSNGLECGFpaSmqVKkiFAqFACAAiMWLF6v6eHp6ik2bNqktN3/+fFWiezYxajNu3DjxxhtvqG336YSmbT2nTp0SMplMJCYmCiGEKCwsFDVr1hTLly8XQgixb98+4eDgIHJyctTWU69ePfGf//xHaxwlxaotwT6d3AsKCkSVKlXEsGHDVG0pKSkCgIiJiRFCCPHZZ5+JoKAgtfUmJycLAOLSpUta49mxY4ewtLQURUVFqrbt27cLBwcHkZWVpdH/wYMHQqFQiOjoaLX20aNHi7feeksIIcSMGTOEt7e3yMvL07pNDw8P8fnnn6u1/d///Z8YN26cEOLfn9Hq1atVn1+4cEEAEBcvXhRCCPHWW2+JHj16qK1j0KBBagnW3t7+ucmzVatWYs6cOaX2Ielh3tEeD/MO805546lYiercuTPi4uJw7NgxTJgwAd27d1dNqL1z5w6Sk5MxevRoVK1aVfVasGCB6hC7NitWrICfnx+cnZ1RtWpVrFq1CklJSXrF1apVKzRq1Ag///wzAODQoUNIS0vDwIEDAQCxsbF48OABHB0d1WJLSEgoNTZdNW/eXPVvS0tLODo6olmzZqo2V1dXAEBaWpoqngMHDqjF0qhRIwAoMZ7Hjx9DLpdDJpOp2rp16wYvLy/UrVsXw4YNw8aNG/Ho0SMAT0655OTkoFu3bmrb2bBhg2obcXFx6NChA6ytrTW2l5WVhVu3bqF9+/Zq7e3bt8fFixdLHL+7u7vaWC9evAh/f3+1/s++DwkJwbvvvouuXbviyy+/1PozsLW1VY2NKhfmHe2Yd5h3ylP5zyylclGlShW89NJLAIClS5eic+fOmDt3LubPn6+6LHzVqlVo27at2nKWlpZa1/ff//4XU6ZMwaJFi+Dv7w97e3t8/fXXOHbsmN6xDR06FJs2bcL06dOxadMmdO/eHU5OTgCezNlwd3fHwYMHNZYzxITkZxOUTCZTaytOisU/o6KiIvTp0wcLFy7UWFdxgnqWk5MTHj16hLy8PNjY2AAA7O3tcerUKRw8eBB79+7FrFmzMGfOHJw4cUK1rV27dqFmzZpq65LL5QCeJK3neTqhA0/mzzzbVtpYhRDP3cacOXMwZMgQ7Nq1C3/88Qdmz56NzZs3o3///qo+d+/ehbOz83PXRdLDvKMd8w7zTnliYVdJzJ49Gz179sTYsWPh4eGBmjVr4tq1axg6dKhOy0dFRSEgIEDtCrdnvzXZ2NigsLDwuesaMmQIPv30U8TGxmLbtm0IDw9Xfda6dWukpqbCysoKderU0Sm24iSmy7b11bp1a2zfvh116tTR+Qq7li1bAnjyjbj43wBgZWWFrl27omvXrpg9ezaqVauG/fv3o1u3bpDL5UhKSkJgYKDWdTZv3hw//PAD8vPzNf6TcHBwgIeHBw4fPoyOHTuq2qOjo9GmTRudx9qkSRMcPXpUre3Z9wDQoEEDNGjQAFOmTMFbb72FdevWqRJsTk4Orl69ilatWum8XZIu5p2yYd5h3nkRPBVbSXTq1AlNmzZV3bhyzpw5CA0NxXfffYfLly/j3LlzWLduHRYvXqx1+ZdeegknT57Enj17cPnyZXz22Wc4ceKEWp86derg7NmzuHTpEtLT05Gfn691Xd7e3ggICMDo0aNRUFCA1157TfVZ165d4e/vj379+mHPnj1ITExEdHQ0Pv30U7WroJ7m4uICW1tb7N69G7dv30ZmZmZZfkRaffDBB7h79y7eeust1VVXe/fuxTvvvFNiQnd2dkbr1q1x+PBhVdvvv/+OpUuXIi4uDtevX8eGDRtQVFSEhg0bwt7eHtOmTcOUKVPwww8/4OrVqzh9+jSWL1+OH374AQAwfvx4ZGVlYfDgwTh58iSuXLmCH3/8EZcuXQIAfPjhh1i4cCG2bNmCS5cuYfr06YiLi8OkSZN0HuvEiROxe/dufPXVV7h8+TK+//57tSvcHj9+jPHjx+PgwYO4fv06jhw5ghMnTqBx48aqPkePHoVcLtc4lUKVE/NO2TDvMO+8ENNO8SNj0DaZWAghNm7cKGxsbERSUpLqfcuWLYWNjY2oXr266Nixo/jll1+EEJoTg3NycsTIkSOFUqkU1apVE2PHjhXTp08XLVq0UK0/LS1NdOvWTVStWlUAEAcOHChxgvHy5csFADF8+HCNOLOyssSECROEh4eHsLa2Fp6enmLo0KGquLVZtWqV8PT0FBYWFiIwMFDrzyEwMFBMmjRJbTkvLy/x7bffqrUBEDt27FC9v3z5sujfv7+oVq2asLW1FY0aNRKTJ09Wm6T8rBUrVoh27dqp3kdFRYnAwEBRvXp1YWtrK5o3b666sk+IJ1flfffdd6Jhw4bC2tpaODs7i+7du4tDhw6p+pw5c0YEBQUJOzs7YW9vLzp06CCuXr0qhHgyGXzu3LmiZs2awtraWrRo0UL88ccfqmW1/R7u3bun+j0VW7NmjahVq5awtbUVffr0Ed98841qEnNubq4YPHiw8PT0FDY2NsLDw0OMHz9ePH78WLX8+++/L8aMGVPiz4Wki3knUOvPgXmHeae8yYTQ4QQ3EeklJycHDRs2xObNmyvNt8g7d+6gUaNGOHnypOoeX0RUfph3mHcAnoolMgqFQoENGzaUekNRqUlISEBYWBiTK5GJMO8QAPCIHREREZFE8IgdERERkUSwsCMiIiKSCBZ2RERERBLBwo6IiIhIIljYEREREUkECzsiIiIiiWBhR0RERCQRLOyIiIiIJIKFHREREZFEsLAjIiIikggWdkREREQSwcKOiIiISCJY2BERERFJBAs7IiIiIolgYUdEREQkESzsiIiIiCTCrAq7v/76C3369IGHhwdkMhl+/fXX5y5z6NAh+Pr6QqFQoG7dulixYoXxAyUiIiIyAbMq7B4+fIgWLVrg+++/16l/QkICevXqhQ4dOuD06dP45JNPMHHiRGzfvt3IkRIRERGVP5kQQpg6iLKQyWTYsWMH+vXrV2Kfjz/+GDt37sTFixdVbcHBwThz5gxiYmLKIUoiIiKi8mNl6gCMKSYmBkFBQWpt3bt3x5o1a5Cfnw9ra2uNZXJzc5Gbm6t6X1RUhLt378LR0REymczoMRNR+RJCIDs7Gx4eHrCwMN1JDOYeosrFWLlH0oVdamoqXF1d1dpcXV1RUFCA9PR0uLu7aywTGhqKuXPnlleIRFRBJCcno1atWibbPnMPUeVk6Nwj6cIOgMY33eIzzyV9A54xYwZCQkJU7zMzM1G7dm0kJyfDwcHBeIESkUlkZWXB09MT9vb2Jo2DuYeocjFW7pF0Yefm5obU1FS1trS0NFhZWcHR0VHrMnK5HHK5XKPdwcGByZVIwkx9upO5h6hyMnTuMaurYvXl7++PyMhItba9e/fCz89P6/w6IiIiInNmVoXdgwcPEBcXh7i4OABPbmcSFxeHpKQkAE9OZQwfPlzVPzg4GNevX0dISAguXryItWvXYs2aNZg2bZopwiciIiIyKrM6FXvy5El07txZ9b54PsqIESOwfv16pKSkqIo8APD29kZERASmTJmC5cuXw8PDA0uXLsUbb7xR7rETERERGZvZ3seuvGRlZUGpVCIzM5PzXIgkqKLu4xU1LiIyDGPt42Z1KpaIiIiISsbCjoiIiEgiWNgRERERSQQLOyIiIiKJYGFHREREJBEs7IiIiIgkgoUdERERkUSwsCMiIiKSCBZ2RERERBLBwo6IiIhIIljYEREREUkECzsiIiIiiWBhR0RERCQRLOyIiIiIJIKFHREREZFEsLAjIiIikggWdkREREQSwcKOiIiISCJY2BERERFJBAs7IiIiIolgYUdEREQkESzsiIiIiCSChR0RERGRRJhdYRcWFgZvb28oFAr4+voiKiqq1P4bN25EixYtYGdnB3d3d4waNQoZGRnlFC0RERFR+TGrwm7Lli2YPHkyZs6cidOnT6NDhw7o2bMnkpKStPY/fPgwhg8fjtGjR+PChQvYunUrTpw4gXfffbecIyciIiIyPrMq7BYvXozRo0fj3XffRePGjbFkyRJ4enoiPDxca/+jR4+iTp06mDhxIry9vfHyyy9jzJgxOHnyZDlHTkRERGR8ZlPY5eXlITY2FkFBQWrtQUFBiI6O1rpMQEAAbty4gYiICAghcPv2bWzbtg29e/cucTu5ubnIyspSexERGRtzDxEZgtkUdunp6SgsLISrq6tau6urK1JTU7UuExAQgI0bN2LQoEGwsbGBm5sbqlWrhmXLlpW4ndDQUCiVStXL09PToOMgItKGuYeIDMFsCrtiMplM7b0QQqOtWHx8PCZOnIhZs2YhNjYWu3fvRkJCAoKDg0tc/4wZM5CZmal6JScnGzR+IiJtmHuIyBCsTB2ArpycnGBpaalxdC4tLU3jKF6x0NBQtG/fHh9++CEAoHnz5qhSpQo6dOiABQsWwN3dXWMZuVwOuVxu+AEQEZWCuYeIDMFsjtjZ2NjA19cXkZGRau2RkZEICAjQusyjR49gYaE+REtLSwBPjvQRERERSYnZFHYAEBISgtWrV2Pt2rW4ePEipkyZgqSkJNWp1RkzZmD48OGq/n369MEvv/yC8PBwXLt2DUeOHMHEiRPRpk0beHh4mGoYREREREZhNqdiAWDQoEHIyMjAvHnzkJKSAh8fH0RERMDLywsAkJKSonZPu5EjRyI7Oxvff/89pk6dimrVqqFLly5YuHChqYZAREREZDQywXOSpcrKyoJSqURmZiYcHBxMHQ4RGVhF3ccralxEZBjG2sfN6lQsEREREZWMhR0RERGRRLCwIyIiIpIIFnZEREREEsHCjoiIiEgiWNgRERERSQQLOyIiIiKJYGFHREREJBEs7IiIiIgkgoUdERERkUSwsCMiIiKSCBZ2RERERBLBwo6IiIhIIljYEREREUkECzsiIiIiiWBhR0RERCQRLOyIiIiIJIKFHREREZFEsLAjIiIikggWdkREREQSwcKOiIiISCJY2BERERFJBAs7IiIiIokwu8IuLCwM3t7eUCgU8PX1RVRUVKn9c3NzMXPmTHh5eUEul6NevXpYu3ZtOUVLREREVH6sTB2APrZs2YLJkycjLCwM7du3x3/+8x/07NkT8fHxqF27ttZlBg4ciNu3b2PNmjV46aWXkJaWhoKCgnKOnIiIiMj4ZEIIYeogdNW2bVu0bt0a4eHhqrbGjRujX79+CA0N1ei/e/duDB48GNeuXUONGjXKtM2srCwolUpkZmbCwcGhzLETUcVUUffxihoXERmGsfZxszkVm5eXh9jYWAQFBam1BwUFITo6WusyO3fuhJ+fH7766ivUrFkTDRo0wLRp0/D48ePyCJmIiIioXJnNqdj09HQUFhbC1dVVrd3V1RWpqalal7l27RoOHz4MhUKBHTt2ID09HePGjcPdu3dLnGeXm5uL3Nxc1fusrCzDDYKIqATMPURkCGZzxK6YTCZTey+E0GgrVlRUBJlMho0bN6JNmzbo1asXFi9ejPXr15d41C40NBRKpVL18vT0NPgYiIiexdxDRIZgNoWdk5MTLC0tNY7OpaWlaRzFK+bu7o6aNWtCqVSq2ho3bgwhBG7cuKF1mRkzZiAzM1P1Sk5ONtwgiIhKwNxDRIZgNoWdjY0NfH19ERkZqdYeGRmJgIAArcu0b98et27dwoMHD1Rtly9fhoWFBWrVqqV1GblcDgcHB7UXEZGxMfcQkSGYTWEHACEhIVi9ejXWrl2LixcvYsqUKUhKSkJwcDCAJ994hw8fruo/ZMgQODo6YtSoUYiPj8dff/2FDz/8EO+88w5sbW1NNQwiIiIiozCbiycAYNCgQcjIyMC8efOQkpICHx8fREREwMvLCwCQkpKCpKQkVf+qVasiMjISEyZMgJ+fHxwdHTFw4EAsWLDAVEMgIiIiMhqzuo+dKfBeUkTSVlH38YoaFxEZRqW/jx0RERERlY6FHREREZFEsLAjIiIikggWdkREREQSwcKOiIiISCJY2BERERFJBAs7IiIiIolgYUdEREQkESzsiIiIiCSChR0RERGRRLCwIyIiIpIIFnZEREREEsHCjoiIiEgiWNgRERERSQQLOyIiIiKJYGFHREREJBEs7IiIiIgkgoUdERERkUSwsCMiIiKSCBZ2RERERBLBwo6IiIhIIljYEREREUkECzsiIiIiiTC7wi4sLAze3t5QKBTw9fVFVFSUTssdOXIEVlZWaNmypXEDJCIiIjIRsyrstmzZgsmTJ2PmzJk4ffo0OnTogJ49eyIpKanU5TIzMzF8+HC88sor5RQpERERUfkzq8Ju8eLFGD16NN599100btwYS5YsgaenJ8LDw0tdbsyYMRgyZAj8/f3LKVIiIiKi8mc2hV1eXh5iY2MRFBSk1h4UFITo6OgSl1u3bh2uXr2K2bNnGztEIiIiIpOyMnUAukpPT0dhYSFcXV3V2l1dXZGamqp1mStXrmD69OmIioqClZVuQ83NzUVubq7qfVZWVtmDJiLSEXMPERmC2RyxKyaTydTeCyE02gCgsLAQQ4YMwdy5c9GgQQOd1x8aGgqlUql6eXp6vnDMRETPw9xDRIYgE0IIUwehi7y8PNjZ2WHr1q3o37+/qn3SpEmIi4vDoUOH1Prfv38f1atXh6WlpaqtqKgIQghYWlpi79696NKli8Z2tH1r9vT0RGZmJhwcHIwwMiIypaysLCiVSpPv48w9RJWLsXKP2ZyKtbGxga+vLyIjI9UKu8jISLz22msa/R0cHHDu3Dm1trCwMOzfvx/btm2Dt7e31u3I5XLI5XLDBk9E9BzMPURkCGZT2AFASEgIhg0bBj8/P/j7+2PlypVISkpCcHAwAGDGjBm4efMmNmzYAAsLC/j4+Kgt7+LiAoVCodFOREREJAVmVdgNGjQIGRkZmDdvHlJSUuDj44OIiAh4eXkBAFJSUp57TzsiIiIiqTKbOXamUlHm3xCRcVTUfbyixkVEhmGsfdzsroolIiIiIu1Y2BERERFJBAs7IiIiIolgYUdEREQkESzsiIiIiCSChR0RERGRRLCwIyIiIpIIFnZEREREEsHCjoiIiEgiWNgRERERSQQLOyIiIiKJYGFHREREJBEs7IiIiIgkgoUdERERkUSwsCMiIiKSCBZ2RERERBLBwo6IiIhIIljYEREREUkECzsiIiIiiWBhR0RERCQRLOyIiIiIJIKFHREREZFEsLAjIiIikggWdkREREQSYXaFXVhYGLy9vaFQKODr64uoqKgS+/7yyy/o1q0bnJ2d4eDgAH9/f+zZs6ccoyUiIiIqP2ZV2G3ZsgWTJ0/GzJkzcfr0aXTo0AE9e/ZEUlKS1v5//fUXunXrhoiICMTGxqJz587o06cPTp8+Xc6RExERERmfTAghTB2Ertq2bYvWrVsjPDxc1da4cWP069cPoaGhOq2jadOmGDRoEGbNmqVT/6ysLCiVSmRmZsLBwaFMcRNRxVVR9/GKGhcRGYax9nErg63JyPLy8hAbG4vp06ertQcFBSE6OlqndRQVFSE7Oxs1atQosU9ubi5yc3NV77OyssoWMBGRHph7iMgQzOZUbHp6OgoLC+Hq6qrW7urqitTUVJ3WsWjRIjx8+BADBw4ssU9oaCiUSqXq5enp+UJxExHpgrmHiAzBbAq7YjKZTO29EEKjTZuff/4Zc+bMwZYtW+Di4lJivxkzZiAzM1P1Sk5OfuGYiYieh7mHiAzBbE7FOjk5wdLSUuPoXFpamsZRvGdt2bIFo0ePxtatW9G1a9dS+8rlcsjl8heOl4hIH8w9RGQIZnPEzsbGBr6+voiMjFRrj4yMREBAQInL/fzzzxg5ciQ2bdqE3r17GztMIiIiIpMxmyN2ABASEoJhw4bBz88P/v7+WLlyJZKSkhAcHAzgyamMmzdvYsOGDQCeFHXDhw/Hd999h3bt2qmO9tna2kKpVJpsHERERETGYFaF3aBBg5CRkYF58+YhJSUFPj4+iIiIgJeXFwAgJSVF7Z52//nPf1BQUIAPPvgAH3zwgap9xIgRWL9+fXmHT0RERGRUZnUfO1PgvaSIpK2i7uMVNS4iMgxj7eNmM8eOiIiIiErHwo6IiIhIIljYEREREUkECzsiIiIiiWBhR0RERCQRLOyIiIiIJIKFHREREZFEsLAjIiIikggWdkREREQSwcKOiIiISCJY2BERERFJBAs7IiIiIolgYUdEREQkESzsiIiIiCSChR0RERGRRLCwIyIiIpIIFnZEREREEsHCjoiIiEgiWNgRERERSQQLOyIiIiKJYGFHREREJBEs7IiIiIgkgoUdERERkUSYXWEXFhYGb29vKBQK+Pr6IioqqtT+hw4dgq+vLxQKBerWrYsVK1aUU6RERERE5cusCrstW7Zg8uTJmDlzJk6fPo0OHTqgZ8+eSEpK0to/ISEBvXr1QocOHXD69Gl88sknmDhxIrZv317OkRMREREZn0wIIUwdhK7atm2L1q1bIzw8XNXWuHFj9OvXD6GhoRr9P/74Y+zcuRMXL15UtQUHB+PMmTOIiYnRaZtZWVlQKpXIzMyEg4PDiw+CiCqUirqPV9S4iMgwjLWPWxlsTUaWl5eH2NhYTJ8+Xa09KCgI0dHRWpeJiYlBUFCQWlv37t2xZs0a5Ofnw9raWmOZ3Nxc5Obmqt5nZmYCePILICLpKd63Tf0dl7mHqHIxVu4xm8IuPT0dhYWFcHV1VWt3dXVFamqq1mVSU1O19i8oKEB6ejrc3d01lgkNDcXcuXM12j09PV8geiKq6DIyMqBUKk22feYeosrJ0LnHbAq7YjKZTO29EEKj7Xn9tbUXmzFjBkJCQlTv79+/Dy8vLyQlJZk06b+orKwseHp6Ijk52axP63AcFYsUxpGZmYnatWujRo0aJo1DirlHCn8fAMdR0UhlHMbKPWZT2Dk5OcHS0lLj6FxaWprGUblibm5uWvtbWVnB0dFR6zJyuRxyuVyjXalUmvUfUDEHBweOowLhOCoOCwvTXksm5dwjhb8PgOOoaKQyDkPnHrO5KtbGxga+vr6IjIxUa4+MjERAQIDWZfz9/TX67927F35+flrn1xERERGZM7Mp7AAgJCQEq1evxtq1a3Hx4kVMmTIFSUlJCA4OBvDkVMbw4cNV/YODg3H9+nWEhITg4sWLWLt2LdasWYNp06aZaghERERERmM2p2IBYNCgQcjIyMC8efOQkpICHx8fREREwMvLCwCQkpKidk87b29vREREYMqUKVi+fDk8PDywdOlSvPHGGzpvUy6XY/bs2VpPkZgTjqNi4Tgqjoo6hooalz6kMAaA46hoOI7SmdV97IiIiIioZGZ1KpaIiIiISsbCjoiIiEgiWNgRERERSQQLOyIiIiKJYGEHICwsDN7e3lAoFPD19UVUVFSp/Q8dOgRfX18oFArUrVsXK1asKKdIS6fPOH755Rd069YNzs7OcHBwgL+/P/bs2VOO0ZZM399HsSNHjsDKygotW7Y0boA60HcMubm5mDlzJry8vCCXy1GvXj2sXbu2nKItmb7j2LhxI1q0aAE7Ozu4u7tj1KhRyMjIKKdotfvrr7/Qp08feHh4QCaT4ddff33uMuWxjzPvMO8YA3NPxcg9Js07opLbvHmzsLa2FqtWrRLx8fFi0qRJokqVKuL69eta+1+7dk3Y2dmJSZMmifj4eLFq1SphbW0ttm3bVs6Rq9N3HJMmTRILFy4Ux48fF5cvXxYzZswQ1tbW4tSpU+UcuTp9x1Hs/v37om7duiIoKEi0aNGifIItQVnG0LdvX9G2bVsRGRkpEhISxLFjx8SRI0fKMWpN+o4jKipKWFhYiO+++05cu3ZNREVFiaZNm4p+/fqVc+TqIiIixMyZM8X27dsFALFjx45S+5fHPs68w7xjDMw9FSf3mDLvVPrCrk2bNiI4OFitrVGjRmL69Ola+3/00UeiUaNGam1jxowR7dq1M1qMutB3HNo0adJEzJ0719Ch6aWs4xg0aJD49NNPxezZs02eYPUdwx9//CGUSqXIyMgoj/B0pu84vv76a1G3bl21tqVLl4patWoZLUZ96ZJgy2MfZ975F/OO4TD3/Ksi5Z7yzjuV+lRsXl4eYmNjERQUpNYeFBSE6OhorcvExMRo9O/evTtOnjyJ/Px8o8VamrKM41lFRUXIzs426YPQyzqOdevW4erVq5g9e7axQ3yusoxh586d8PPzw1dffYWaNWuiQYMGmDZtGh4/flweIWtVlnEEBATgxo0biIiIgBACt2/fxrZt29C7d+/yCNlgjL2PM+/8i3nHcJh7zDv3GHIfN6snTxhaeno6CgsL4erqqtbu6uqK1NRUrcukpqZq7V9QUID09HS4u7sbLd6SlGUcz1q0aBEePnyIgQMHGiNEnZRlHFeuXMH06dMRFRUFKyvT/zmXZQzXrl3D4cOHoVAosGPHDqSnp2PcuHG4e/euyea6lGUcAQEB2LhxIwYNGoScnBwUFBSgb9++WLZsWXmEbDDG3seZd/7FvGM4zD3mnXsMuY9X6iN2xWQymdp7IYRG2/P6a2svb/qOo9jPP/+MOXPmYMuWLXBxcTFWeDrTdRyFhYUYMmQI5s6diwYNGpRXeDrR53dRVFQEmUyGjRs3ok2bNujVqxcWL16M9evXm/SbM6DfOOLj4zFx4kTMmjULsbGx2L17NxISElTPcjYn5bGPM+8w7xgDc4/55h5D7eMV46uGiTg5OcHS0lLjW0BaWppG5VzMzc1Na38rKys4OjoaLdbSlGUcxbZs2YLRo0dj69at6Nq1qzHDfC59x5GdnY2TJ0/i9OnTGD9+PIAniUoIASsrK+zduxddunQpl9iLleV34e7ujpo1a0KpVKraGjduDCEEbty4gfr16xs1Zm3KMo7Q0FC0b98eH374IQCgefPmqFKlCjp06IAFCxaY5KhSWRh7H2feYd4xBuYe8849htzHK/UROxsbG/j6+iIyMlKtPTIyEgEBAVqX8ff31+i/d+9e+Pn5wdra2mixlqYs4wCefGMeOXIkNm3aVCHmIug7DgcHB5w7dw5xcXGqV3BwMBo2bIi4uDi0bdu2vEJXKcvvon379rh16xYePHigart8+TIsLCxQq1Yto8ZbkrKM49GjR7CwUE8plpaWAP795mkOjL2PM+8w7xgDc4955x6D7uN6X24hMcWXVa9Zs0bEx8eLyZMniypVqojExEQhhBDTp08Xw4YNU/UvviR5ypQpIj4+XqxZs6ZC3XZA13Fs2rRJWFlZieXLl4uUlBTV6/79+6YaghBC/3E8qyJcnabvGLKzs0WtWrXEgAEDxIULF8ShQ4dE/fr1xbvvvmuqIQgh9B/HunXrhJWVlQgLCxNXr14Vhw8fFn5+fqJNmzamGoIQ4snP9/Tp0+L06dMCgFi8eLE4ffq06tYJptjHmXeYd4yBuafi5B5T5p1KX9gJIcTy5cuFl5eXsLGxEa1btxaHDh1SfTZixAgRGBio1v/gwYOiVatWwsbGRtSpU0eEh4eXc8Ta6TOOwMBAAUDjNWLEiPIP/Bn6/j6eVlESrL5juHjxoujatauwtbUVtWrVEiEhIeLRo0flHLUmfcexdOlS0aRJE2Frayvc3d3F0KFDxY0bN8o5anUHDhwo9W/dVPs48w7zjjEw91SM3GPKvCMTwkyOUxIRERFRqSr1HDsiIiIiKWFhR0RERCQRLOyIiIiIJIKFHREREZFEsLAjIiIikggWdkREREQSwcKOiIiISCJY2BERERFJBAs7IiIiIolgYUdEREQkESzsiIiIiCSChR0RERGRRPw/4e4SSvt8LbgAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Copy the data\n",
    "trace_filter = trace_physical.copy()\n",
    "\n",
    "# Filter\n",
    "trace_filter.filter(\"bandpass\", freqmin=1, freqmax=3, corners=4)\n",
    "\n",
    "# Spectrum\n",
    "spectrum_filter = np.fft.rfft(trace_filter.data)\n",
    "\n",
    "# Plot\n",
    "fig, ax = plt.subplots(2, 2, sharey=\"row\")\n",
    "\n",
    "# Plot the original trace\n",
    "ax[0, 0].plot(trace_filter.times(), trace_filter.data)\n",
    "ax[0, 0].set_ylabel(\"Amplitude (counts)\")\n",
    "ax[0, 1].plot(trace_filter.times(), trace_filter.data)\n",
    "\n",
    "# Plot the spectrum\n",
    "# Your code goes here\n",
    "\n",
    "# Labels\n",
    "ax[0, 0].set_xlabel(\"Relative time (seconds)\")\n",
    "ax[0, 1].set_xlabel(\"Relative time (seconds)\")\n",
    "ax[0, 0].set_title(\"Original trace\")\n",
    "ax[0, 1].set_title(\"Filtered trace\")\n",
    "\n",
    "# Show\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3. Spectrogram\n",
    "\n",
    "> __Question 8.__ The following cell performs a spectrogram —a time-dependant Fourier spectrum analysis $|\\tilde u(f, t)|$— of the raw trace. What are the `nperseg` and the `noverlap` parameters useful for? Describe at least four different signal patterns seen in the spectrogram, depending on their duration, frequency bandwidth, etc. Two waves exhibit a typical signature in the spectrogram, with an increasing frequency against time. What does this imply for the speed of the seismic waves? How does this spectrogram changes with respect to the input waveform units?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract waveform\n",
    "waveform = trace.data\n",
    "times = trace.times()\n",
    "\n",
    "# Calculate spectrogram\n",
    "frequencies, spectrogram_times, psd = spectrogram(\n",
    "    waveform, fs=trace.stats.sampling_rate, nperseg=4000, noverlap=2500\n",
    ")\n",
    "\n",
    "# What's next?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Analysis with a selected event.\n",
    "\n",
    "The time has come. Re-run the entire notebook with a event of your choice. Motivate your choice and anwser most of the questions with the new event."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.9"
  },
  "vscode": {
   "interpreter": {
    "hash": "a41af5d3c6566f72458df644874e80e2f0fb14f6983016ff4dfaa9ae44ef2b04"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}