docs/tutorials/data_model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparison of data and model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial, we will compare the photometric fluxes of the brown dwarf companion PZ Tel B with a synthetic spectrum from the ATMO grid."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initiating *species*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by importing *species*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import species"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And initiating the workflow with the `SpeciesInit` class. This will create the configuration file and the HDF5 database."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initiating species v0.3.1... [DONE]\n",
      "Creating species_config.ini... [DONE]\n",
      "Database: /Users/tomasstolker/applications/species/docs/tutorials/species_database.hdf5\n",
      "Data folder: /Users/tomasstolker/applications/species/docs/tutorials/data\n",
      "Working folder: /Users/tomasstolker/applications/species/docs/tutorials\n",
      "Creating species_database.hdf5... [DONE]\n",
      "Creating data folder... [DONE]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<species.core.setup.SpeciesInit at 0x144c1b1d0>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "species.SpeciesInit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adding model spectra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create now a `Database` object which is used for importing various data into the database."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "database = species.Database()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The spectra of ATMO are downloaded and added to the database with the `add_model` method of `Database`. This requires sufficient disk storage in the *data_folder* that is set in the configuration file. The full ATMO grid is downloaded but the `teff_range` parameter can be used to only import a certain T<sub>eff</sub> range into the database."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Downloading ATMO model spectra (430 MB)... [DONE]\n",
      "Unpacking ATMO model spectra (430 MB)... [DONE]\n",
      "Adding ATMO model spectra... [DONE]                          \n",
      "Grid points stored in the database:\n",
      "   - Teff = [2500. 2600. 2700. 2800. 2900. 3000.]\n",
      "   - log(g) = [2.5 3.  3.5 4.  4.5 5.  5.5]\n",
      "Number of grid points per parameter:\n",
      "   - teff: 6\n",
      "   - logg: 7\n",
      "Fix missing grid points with a linear interpolation:\n",
      "Number of stored grid points: 42\n",
      "Number of interpolated grid points: 0\n",
      "Number of missing grid points: 0\n"
     ]
    }
   ],
   "source": [
    "database.add_model('atmo', teff_range=(2500., 3000.))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adding companion data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we add the distance and magnitudes of PZ Tel B to the database with the `add_companion`` method. This will automatically download the required filter profiles and a flux-calibrated spectrum of Vega. These are used to convert the magnitudes into fluxes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Downloading Vega spectrum (270 kB)... [DONE]\n",
      "Adding Vega spectrum... [DONE]\n",
      "Adding filter: Paranal/SPHERE.ZIMPOL_R_PRIM... [DONE]\n",
      "Adding filter: Paranal/SPHERE.ZIMPOL_I_PRIM... [DONE]\n",
      "Adding filter: Paranal/SPHERE.IRDIS_D_H23_2... [DONE]\n",
      "Adding filter: Paranal/SPHERE.IRDIS_D_H23_3... [DONE]\n",
      "Adding filter: Paranal/SPHERE.IRDIS_D_K12_1... [DONE]\n",
      "Adding filter: Paranal/SPHERE.IRDIS_D_K12_2... [DONE]\n",
      "Adding filter: Paranal/NACO.J... [DONE]\n",
      "Adding filter: Paranal/NACO.H... [DONE]\n",
      "Adding filter: Paranal/NACO.Ks... [DONE]\n",
      "Adding filter: Paranal/NACO.Lp... [DONE]\n",
      "Adding filter: Paranal/NACO.NB405... [DONE]\n",
      "Adding filter: Paranal/NACO.Mp... [DONE]\n",
      "Adding filter: Gemini/NICI.ED286... [DONE]\n",
      "Adding filter: Gemini/NIRI.H2S1v2-1-G0220... [DONE]\n",
      "Adding object: PZ Tel B\n",
      "   - Distance (pc) = 47.13 +/- 0.13\n",
      "   - Paranal/SPHERE.ZIMPOL_R_PRIM:\n",
      "      - Apparent magnitude = 17.84 +/- 0.31\n",
      "      - Flux (W m-2 um-1) = 1.83e-15 +/- 5.29e-16\n",
      "   - Paranal/SPHERE.ZIMPOL_I_PRIM:\n",
      "      - Apparent magnitude = 15.16 +/- 0.12\n",
      "      - Flux (W m-2 um-1) = 1.09e-14 +/- 1.20e-15\n",
      "   - Paranal/SPHERE.IRDIS_D_H23_2:\n",
      "      - Apparent magnitude = 11.78 +/- 0.19\n",
      "      - Flux (W m-2 um-1) = 2.54e-14 +/- 4.47e-15\n",
      "   - Paranal/SPHERE.IRDIS_D_H23_3:\n",
      "      - Apparent magnitude = 11.65 +/- 0.19\n",
      "      - Flux (W m-2 um-1) = 2.43e-14 +/- 4.27e-15\n",
      "   - Paranal/SPHERE.IRDIS_D_K12_1:\n",
      "      - Apparent magnitude = 11.56 +/- 0.09\n",
      "      - Flux (W m-2 um-1) = 1.15e-14 +/- 9.58e-16\n",
      "   - Paranal/SPHERE.IRDIS_D_K12_2:\n",
      "      - Apparent magnitude = 11.29 +/- 0.10\n",
      "      - Flux (W m-2 um-1) = 1.14e-14 +/- 1.05e-15\n",
      "   - Paranal/NACO.J:\n",
      "      - Apparent magnitude = 12.47 +/- 0.20\n",
      "      - Flux (W m-2 um-1) = 3.11e-14 +/- 5.76e-15\n",
      "   - Paranal/NACO.H:\n",
      "      - Apparent magnitude = 11.93 +/- 0.14\n",
      "      - Flux (W m-2 um-1) = 1.97e-14 +/- 2.55e-15\n",
      "   - Paranal/NACO.Ks:\n",
      "      - Apparent magnitude = 11.53 +/- 0.07\n",
      "      - Flux (W m-2 um-1) = 1.12e-14 +/- 7.25e-16\n",
      "   - Paranal/NACO.Lp:\n",
      "      - Apparent magnitude = 11.04 +/- 0.22\n",
      "      - Flux (W m-2 um-1) = 2.02e-15 +/- 4.12e-16\n",
      "   - Paranal/NACO.NB405:\n",
      "      - Apparent magnitude = 10.94 +/- 0.07\n",
      "      - Flux (W m-2 um-1) = 1.67e-15 +/- 1.08e-16\n",
      "   - Paranal/NACO.Mp:\n",
      "      - Apparent magnitude = 10.93 +/- 0.03\n",
      "      - Flux (W m-2 um-1) = 9.19e-16 +/- 2.54e-17\n",
      "   - Gemini/NICI.ED286:\n",
      "      - Apparent magnitude = 11.68 +/- 0.14\n",
      "      - Flux (W m-2 um-1) = 2.78e-14 +/- 3.60e-15\n",
      "   - Gemini/NIRI.H2S1v2-1-G0220:\n",
      "      - Apparent magnitude = 11.39 +/- 0.14\n",
      "      - Flux (W m-2 um-1) = 1.06e-14 +/- 1.37e-15\n"
     ]
    }
   ],
   "source": [
    "database.add_companion('PZ Tel B')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, the `add_object` method of `Database` can be used for manually adding magnitudes and spectra of an individual object. Before coninuing, let's check the content of the database."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Database content:\n",
      "- filters: <HDF5 group \"/filters\" (2 members)>\n",
      "\t- Gemini: <HDF5 group \"/filters/Gemini\" (2 members)>\n",
      "\t\t- NICI.ED286: <HDF5 dataset \"NICI.ED286\": shape (387, 2), type \"<f4\">\n",
      "\t\t\t- det_type: energy\n",
      "\t\t- NIRI.H2S1v2-1-G0220: <HDF5 dataset \"NIRI.H2S1v2-1-G0220\": shape (129, 2), type \"<f4\">\n",
      "\t\t\t- det_type: energy\n",
      "\t- Paranal: <HDF5 group \"/filters/Paranal\" (12 members)>\n",
      "\t\t- NACO.H: <HDF5 dataset \"NACO.H\": shape (23, 2), type \"<f4\">\n",
      "\t\t\t- det_type: energy\n",
      "\t\t- NACO.J: <HDF5 dataset \"NACO.J\": shape (20, 2), type \"<f4\">\n",
      "\t\t\t- det_type: energy\n",
      "\t\t- NACO.Ks: <HDF5 dataset \"NACO.Ks\": shape (27, 2), type \"<f4\">\n",
      "\t\t\t- det_type: energy\n",
      "\t\t- NACO.Lp: <HDF5 dataset \"NACO.Lp\": shape (31, 2), type \"<f4\">\n",
      "\t\t\t- det_type: energy\n",
      "\t\t- NACO.Mp: <HDF5 dataset \"NACO.Mp\": shape (18, 2), type \"<f4\">\n",
      "\t\t\t- det_type: energy\n",
      "\t\t- NACO.NB405: <HDF5 dataset \"NACO.NB405\": shape (67, 2), type \"<f4\">\n",
      "\t\t\t- det_type: energy\n",
      "\t\t- SPHERE.IRDIS_D_H23_2: <HDF5 dataset \"SPHERE.IRDIS_D_H23_2\": shape (113, 2), type \"<f4\">\n",
      "\t\t\t- det_type: energy\n",
      "\t\t- SPHERE.IRDIS_D_H23_3: <HDF5 dataset \"SPHERE.IRDIS_D_H23_3\": shape (180, 2), type \"<f4\">\n",
      "\t\t\t- det_type: energy\n",
      "\t\t- SPHERE.IRDIS_D_K12_1: <HDF5 dataset \"SPHERE.IRDIS_D_K12_1\": shape (175, 2), type \"<f4\">\n",
      "\t\t\t- det_type: energy\n",
      "\t\t- SPHERE.IRDIS_D_K12_2: <HDF5 dataset \"SPHERE.IRDIS_D_K12_2\": shape (191, 2), type \"<f4\">\n",
      "\t\t\t- det_type: energy\n",
      "\t\t- SPHERE.ZIMPOL_I_PRIM: <HDF5 dataset \"SPHERE.ZIMPOL_I_PRIM\": shape (189, 2), type \"<f4\">\n",
      "\t\t\t- det_type: energy\n",
      "\t\t- SPHERE.ZIMPOL_R_PRIM: <HDF5 dataset \"SPHERE.ZIMPOL_R_PRIM\": shape (169, 2), type \"<f4\">\n",
      "\t\t\t- det_type: energy\n",
      "- models: <HDF5 group \"/models\" (1 members)>\n",
      "\t- atmo: <HDF5 group \"/models/atmo\" (4 members)>\n",
      "\t\t- flux: <HDF5 dataset \"flux\": shape (6, 7, 98116), type \"<f8\">\n",
      "\t\t- logg: <HDF5 dataset \"logg\": shape (7,), type \"<f8\">\n",
      "\t\t- teff: <HDF5 dataset \"teff\": shape (6,), type \"<f8\">\n",
      "\t\t- wavelength: <HDF5 dataset \"wavelength\": shape (98116,), type \"<f8\">\n",
      "- objects: <HDF5 group \"/objects\" (1 members)>\n",
      "\t- PZ Tel B: <HDF5 group \"/objects/PZ Tel B\" (3 members)>\n",
      "\t\t- Gemini: <HDF5 group \"/objects/PZ Tel B/Gemini\" (2 members)>\n",
      "\t\t\t- NICI.ED286: <HDF5 dataset \"NICI.ED286\": shape (4,), type \"<f8\">\n",
      "\t\t\t\t- n_phot: 1\n",
      "\t\t\t- NIRI.H2S1v2-1-G0220: <HDF5 dataset \"NIRI.H2S1v2-1-G0220\": shape (4,), type \"<f8\">\n",
      "\t\t\t\t- n_phot: 1\n",
      "\t\t- Paranal: <HDF5 group \"/objects/PZ Tel B/Paranal\" (12 members)>\n",
      "\t\t\t- NACO.H: <HDF5 dataset \"NACO.H\": shape (4,), type \"<f8\">\n",
      "\t\t\t\t- n_phot: 1\n",
      "\t\t\t- NACO.J: <HDF5 dataset \"NACO.J\": shape (4,), type \"<f8\">\n",
      "\t\t\t\t- n_phot: 1\n",
      "\t\t\t- NACO.Ks: <HDF5 dataset \"NACO.Ks\": shape (4,), type \"<f8\">\n",
      "\t\t\t\t- n_phot: 1\n",
      "\t\t\t- NACO.Lp: <HDF5 dataset \"NACO.Lp\": shape (4,), type \"<f8\">\n",
      "\t\t\t\t- n_phot: 1\n",
      "\t\t\t- NACO.Mp: <HDF5 dataset \"NACO.Mp\": shape (4,), type \"<f8\">\n",
      "\t\t\t\t- n_phot: 1\n",
      "\t\t\t- NACO.NB405: <HDF5 dataset \"NACO.NB405\": shape (4,), type \"<f8\">\n",
      "\t\t\t\t- n_phot: 1\n",
      "\t\t\t- SPHERE.IRDIS_D_H23_2: <HDF5 dataset \"SPHERE.IRDIS_D_H23_2\": shape (4,), type \"<f8\">\n",
      "\t\t\t\t- n_phot: 1\n",
      "\t\t\t- SPHERE.IRDIS_D_H23_3: <HDF5 dataset \"SPHERE.IRDIS_D_H23_3\": shape (4,), type \"<f8\">\n",
      "\t\t\t\t- n_phot: 1\n",
      "\t\t\t- SPHERE.IRDIS_D_K12_1: <HDF5 dataset \"SPHERE.IRDIS_D_K12_1\": shape (4,), type \"<f8\">\n",
      "\t\t\t\t- n_phot: 1\n",
      "\t\t\t- SPHERE.IRDIS_D_K12_2: <HDF5 dataset \"SPHERE.IRDIS_D_K12_2\": shape (4,), type \"<f8\">\n",
      "\t\t\t\t- n_phot: 1\n",
      "\t\t\t- SPHERE.ZIMPOL_I_PRIM: <HDF5 dataset \"SPHERE.ZIMPOL_I_PRIM\": shape (4,), type \"<f8\">\n",
      "\t\t\t\t- n_phot: 1\n",
      "\t\t\t- SPHERE.ZIMPOL_R_PRIM: <HDF5 dataset \"SPHERE.ZIMPOL_R_PRIM\": shape (4,), type \"<f8\">\n",
      "\t\t\t\t- n_phot: 1\n",
      "\t\t- distance: <HDF5 dataset \"distance\": shape (2,), type \"<f8\">\n",
      "- spectra: <HDF5 group \"/spectra\" (1 members)>\n",
      "\t- calibration: <HDF5 group \"/spectra/calibration\" (1 members)>\n",
      "\t\t- vega: <HDF5 dataset \"vega\": shape (3, 8827), type \"<f8\">\n"
     ]
    }
   ],
   "source": [
    "database.list_content()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see the various groups, subgroups, datasets, and attributes that are stored in the HDF5 database."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading model spectra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Model spectra are read from the database by first creating an instance of `ReadModel`. The model name and optionally a wavelength range are provided as arguments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "readmodel = species.ReadModel('atmo', wavel_range=(0.5, 10.))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before extracting a spectrum, let's check which parameters are required for the ATMO model spectra."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['teff', 'logg']"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "readmodel.get_parameters()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And also the parameter boundaries of the grid that is stored in the database."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'teff': (2500.0, 3000.0), 'logg': (2.5, 5.5)}"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "readmodel.get_bounds()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The parameters are provided in a dictionary for which we have to make sure that chose values are within the grid boundaries. The radius (R<sub>Jup</sub>) and distance (pc) will scale the emitted spectrum to the observer. Without these values, the spectrum fluxes are provided at the surface of the atmosphere."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_param = {'teff': 2900., 'logg': 4.5, 'radius': 2.2, 'distance': 47.13}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now use the `get_model` method of `ReadModel` to linearly interpolate the grid of spectra and store the extracted spectrum in a `ModelBox`. The spectrum is smoothed to a spectral resolution of R = 100."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "modelbox = readmodel.get_model(model_param, spec_res=100., smooth=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading companion data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The photometric data of PZ Tel B are also read from the database and stored in an `ObjectBox`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Getting object: PZ Tel B... [DONE]\n"
     ]
    }
   ],
   "source": [
    "objectbox = database.get_object(object_name='PZ Tel B')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Synthetic photometry for all filters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For comparison, we create synthetic photometry from the extracted ATMO spectrum for all filters of PZ Tel B. The synthetic fluxes are stored in a `SynphotBox`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating synthetic photometry... [DONE]\n"
     ]
    }
   ],
   "source": [
    "synphotbox = species.multi_photometry(datatype='model',\n",
    "                                      spectrum='atmo',\n",
    "                                      filters=objectbox.filters,\n",
    "                                      parameters=model_param)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating flux residuals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `get_residuals` function is now used to calculate the difference between the observed fluxes and the synthetic fluxes from the model spectrum. The residuals are stored in a `ResidualsBox`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating synthetic photometry... [DONE]\n",
      "Calculating residuals... [DONE]\n",
      "Residuals (sigma):\n",
      "   - Gemini/NICI.ED286: 1.35\n",
      "   - Gemini/NIRI.H2S1v2-1-G0220: -0.02\n",
      "   - Paranal/NACO.H: -0.54\n",
      "   - Paranal/NACO.J: -0.19\n",
      "   - Paranal/NACO.Ks: 0.85\n",
      "   - Paranal/NACO.Lp: 0.02\n",
      "   - Paranal/NACO.Mp: 5.45\n",
      "   - Paranal/NACO.NB405: 0.31\n",
      "   - Paranal/SPHERE.IRDIS_D_H23_2: 0.54\n",
      "   - Paranal/SPHERE.IRDIS_D_H23_3: 0.00\n",
      "   - Paranal/SPHERE.IRDIS_D_K12_1: 1.00\n",
      "   - Paranal/SPHERE.IRDIS_D_K12_2: 1.00\n",
      "   - Paranal/SPHERE.ZIMPOL_I_PRIM: -7.73\n",
      "   - Paranal/SPHERE.ZIMPOL_R_PRIM: -3.96\n"
     ]
    }
   ],
   "source": [
    "res_box = species.get_residuals(datatype='model',\n",
    "                                spectrum='atmo',\n",
    "                                parameters=model_param,\n",
    "                                objectbox=objectbox,\n",
    "                                inc_phot=True,\n",
    "                                inc_spec=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Opening Box objects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `open_box` function can be used to view the content of a `Box` object. For example, the `ModelBox` contains several attributes, including the wavelengths and fluxes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Opening ModelBox...\n",
      "model = atmo\n",
      "type = None\n",
      "wavelength = [ 0.4999641   0.50001364  0.50006318 ...  9.99858711  9.99957782\n",
      " 10.00056862]\n",
      "flux = [2.05385551e-15 2.05647449e-15 2.05928733e-15 ... 5.56118275e-17\n",
      " 5.56005834e-17 5.55895317e-17]\n",
      "parameters = {'teff': 2900.0, 'logg': 4.5, 'radius': 2.2, 'distance': 47.13, 'mass': 59.04988128500266, 'luminosity': 0.003114426265448587}\n",
      "quantity = flux\n"
     ]
    }
   ],
   "source": [
    "modelbox.open_box()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, an `ObjectBox` contains a dictionary with the magnitudes and a dictionary with the fluxes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Opening ObjectBox...\n",
      "name = PZ Tel B\n",
      "filters = ['Gemini/NICI.ED286', 'Gemini/NIRI.H2S1v2-1-G0220', 'Paranal/NACO.H', 'Paranal/NACO.J', 'Paranal/NACO.Ks', 'Paranal/NACO.Lp', 'Paranal/NACO.Mp', 'Paranal/NACO.NB405', 'Paranal/SPHERE.IRDIS_D_H23_2', 'Paranal/SPHERE.IRDIS_D_H23_3', 'Paranal/SPHERE.IRDIS_D_K12_1', 'Paranal/SPHERE.IRDIS_D_K12_2', 'Paranal/SPHERE.ZIMPOL_I_PRIM', 'Paranal/SPHERE.ZIMPOL_R_PRIM']\n",
      "magnitude = {'Gemini/NICI.ED286': array([11.68,  0.14]), 'Gemini/NIRI.H2S1v2-1-G0220': array([11.39,  0.14]), 'Paranal/NACO.H': array([11.93,  0.14]), 'Paranal/NACO.J': array([12.47,  0.2 ]), 'Paranal/NACO.Ks': array([11.53,  0.07]), 'Paranal/NACO.Lp': array([11.04,  0.22]), 'Paranal/NACO.Mp': array([10.93,  0.03]), 'Paranal/NACO.NB405': array([10.94,  0.07]), 'Paranal/SPHERE.IRDIS_D_H23_2': array([11.78,  0.19]), 'Paranal/SPHERE.IRDIS_D_H23_3': array([11.65,  0.19]), 'Paranal/SPHERE.IRDIS_D_K12_1': array([11.56,  0.09]), 'Paranal/SPHERE.IRDIS_D_K12_2': array([11.29,  0.1 ]), 'Paranal/SPHERE.ZIMPOL_I_PRIM': array([15.16,  0.12]), 'Paranal/SPHERE.ZIMPOL_R_PRIM': array([17.84,  0.31])}\n",
      "flux = {'Gemini/NICI.ED286': array([2.78256301e-14, 3.59792032e-15]), 'Gemini/NIRI.H2S1v2-1-G0220': array([1.05904383e-14, 1.36936892e-15]), 'Paranal/NACO.H': array([1.96875856e-14, 2.54565177e-15]), 'Paranal/NACO.J': array([3.11068456e-14, 5.76255347e-15]), 'Paranal/NACO.Ks': array([1.12416276e-14, 7.25276730e-16]), 'Paranal/NACO.Lp': array([2.02006233e-15, 4.12126882e-16]), 'Paranal/NACO.Mp': array([9.18778503e-16, 2.53900187e-17]), 'Paranal/NACO.NB405': array([1.67046281e-15, 1.07773345e-16]), 'Paranal/SPHERE.IRDIS_D_H23_2': array([2.54130005e-14, 4.46991835e-15]), 'Paranal/SPHERE.IRDIS_D_H23_3': array([2.42699566e-14, 4.26886720e-15]), 'Paranal/SPHERE.IRDIS_D_K12_1': array([1.15478089e-14, 9.58329864e-16]), 'Paranal/SPHERE.IRDIS_D_K12_2': array([1.14281212e-14, 1.05405765e-15]), 'Paranal/SPHERE.ZIMPOL_I_PRIM': array([1.08628802e-14, 1.20305574e-15]), 'Paranal/SPHERE.ZIMPOL_R_PRIM': array([1.82633135e-15, 5.28569080e-16])}\n",
      "distance = [47.13  0.13]\n",
      "spectrum = None\n"
     ]
    }
   ],
   "source": [
    "objectbox.open_box()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The attributes in a `Box` object can be extracted for further analyis or creating plots. For example, to extract the array with wavelengths from the `ModelBox`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.4999641 ,  0.50001364,  0.50006318, ...,  9.99858711,\n",
       "        9.99957782, 10.00056862])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "modelbox.wavelength"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting model spectrum and photometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we will combine the model spectrum and the photometric fluxes in a plot with `plot_spectrum`. A list with `Box` objects is provided as an argument of `boxes`. These are interpreted accordingly by the `plot_spectrum` function. Also a list with filter names can be provided as argument of `filters` to show the filter profiles. The `ResidualsBox` is provided as arguments of `residuals`. Finally, the optional argument of `plot_kwargs` contains a list with optional dictionaries to tune the visualization of the plotted data. The number of items in the list of `plot_kwargs` should be equal to the number of items in the list of `boxes`. For the `SynphotBox`, we can set the item in `plot_kwargs` to `None` such that the marker design is based on the data from `ObjectBox`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Plotting spectrum: spectrum.png... [DONE]\n"
     ]
    }
   ],
   "source": [
    "species.plot_spectrum(boxes=[modelbox, objectbox, synphotbox],\n",
    "                      filters=objectbox.filters,\n",
    "                      residuals=res_box,\n",
    "                      plot_kwargs=[{'ls': '-', 'lw': 1., 'color': 'black'},\n",
    "                                   {'Gemini/NICI.ED286': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},\n",
    "                                    'Gemini/NIRI.H2S1v2-1-G0220': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},\n",
    "                                    'Paranal/NACO.H': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},\n",
    "                                    'Paranal/NACO.J': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},\n",
    "                                    'Paranal/NACO.Ks': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},\n",
    "                                    'Paranal/NACO.Lp': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},\n",
    "                                    'Paranal/NACO.Mp': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},\n",
    "                                    'Paranal/NACO.NB405': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},\n",
    "                                    'Paranal/SPHERE.IRDIS_D_H23_2': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},\n",
    "                                    'Paranal/SPHERE.IRDIS_D_H23_3': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},\n",
    "                                    'Paranal/SPHERE.IRDIS_D_K12_1': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},\n",
    "                                    'Paranal/SPHERE.IRDIS_D_K12_2': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},\n",
    "                                    'Paranal/SPHERE.ZIMPOL_I_PRIM': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'},\n",
    "                                    'Paranal/SPHERE.ZIMPOL_R_PRIM': {'marker': 's', 'ms': 4., 'color': 'tab:blue', 'ls': 'none'}},\n",
    "                                   None],\n",
    "                      xlim=(0.5, 5.5),\n",
    "                      ylim=(-5e-15, 5.5e-14),\n",
    "                      ylim_res=(-8, 8),\n",
    "                      scale=('linear', 'linear'),\n",
    "                      title='PZ Tel B and ATMO spectrum',\n",
    "                      offset=(-0.45, -0.04),\n",
    "                      legend={'loc': 'upper right', 'frameon': False, 'fontsize': 12},\n",
    "                      figsize=(9., 4),\n",
    "                      quantity='flux',\n",
    "                      output='spectrum.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's have a look at the result. The plot is stored in the working folder. The blue squares are the photometric fluxes of PZ Tel B and the open squares are the synthetic photometry computed from the model spectrum. The residuals are shown relative to the uncertainties on the fluxes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsIAAAFmCAYAAABnd/y8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAACOeklEQVR4nOzdd3hU1fbw8e9KIXRCD4TeRUWk2hAsCIjYrwVEUa/YFTve670E+2u5dsUOFvBnR8UCKqBYQBAQ6QKhBAIkECCEkLbeP6Y4mUySmUxLWZ/nOTCn7b3OzCRZs2efvUVVMcYYY4wxpqaJiXYAxhhjjDHGRIMlwsYYY4wxpkayRNgYY4wxxtRIlggbY4wxxpgayRJhY4wxxhhTI1kibIwxxhhjaiRLhI0xIScig0XkVxFREVkqIvNE5A8R+VJEujqPGScia5z7PJe9IpIjIi18lOt5/DJn+Z7rKWXENNx5TGoZx7wpIunOGOaJyM8islFEpotIo1A8N+URkVPLi9Pr+MdFJFVExGOb5/Ppeh1+9VifKiL/dT6fKiLnllL2gx7njvTY3lpEXvUoc5GI3C8idYK9fmOMiSSxcYSNMeEgIh2ATcAgVV0gIjHAe0Av4CjgMgBVnepxTm9gEXCXqj7jo8x5qjrE+XgIMFdVxWN9iKqmlBHTEGCqqnYo45ipQBtVPd253hJYA7yhqneUc9kh4U+czuNigWU4ns8hqjrfud3zeeqA43XoqKqpzvUUVR0nIuOAF4A/VPV4r7IbAKuBZNe5zu1tgF+Ayar6mnNbAjAF6Aycpqr5wVx/qImI4nENxhjjYi3CxpiIUNUi4B2gO9AD+Ma5ACAitYBpwE/As6UUc28ZVawApockWA+quhNYjyPuyuZ04P+A34HLPbaX9TztBJ7yWH8dOM6ZfHsaD3zo4/zngcWuJBhAVQ8D1+N4jm71M3ZjjIk6S4SNMZEU7/w/X1V3qOoOj30pQEfgSi3lqypV/aW0glU1U1XXichZIvKLiMwXkR9EZFAwAYvI0cCRwJdlHHOZiMwVkW+ddV7qse9JZ3eLZ5xdEn4WkYXOllnXMa2d3UaWiMjnwLF+hjcGx4eHN4ELRaQ2lPs8HVLV5R6bfgdmAxM94qkFjAQ+9rrOxsAo4HMf5eYCc4DRvuoVhyec1/6tiHwlIseJSJKze4WKyG0i8o2zG827zlZp1/nHOV/T+c7n8HyPfTEicrfzdf9eRH4UkbEiUktE5jkPe89Zz1ARucvZnWS6iDznfM2KROR65/apznKP9+6mIiLfiUiuiEwQkU+c3UueFZG2zvKWi8j7IhJX2mtgjKlEVDXkC9A+HOXaYostVWcBOgAKnORcrwd8h6PrQ4zXsQOBAuCaAMof4vgVVmxbHyAXONq53hvYD7T0OCe1nHKnAnuBecBC4CCOFuqYMs4ZDyQ6HycAq3B8Fe9Z5iqgoXP9NeB5j/1zgVedj2OBj/yIswHwifNxY+d1X1LG69DBx75xzmWI85jezu1XAVcDJ3meCwxwrp9RSkwPAwdL2TccR1eLGOf6BTi6aLj2q+s5cT4Hs4HXnOvJQDYw3GM9C+jlXJ8ArAQaeMS51KvsDl7xpAAZOP9eAfcDSc7tU73eZ6le56a6jgHqALudr3GcM/Y/gQui/TNoiy22lL8E3SLs/JTfUkTauRZgRrDlGmOqjeecrXLzgA3AKHV0kwDA2Yo5DZitqq8GWdd44EdVXQGgqsuAvyillbIMS1R1iKoOBFrj6IP7aRnH/w5ME5EFOLp7tAb6eR0zR1X3Ox8vBlw3DbbBkWy95oy5EP9+h16Io1sEqroXmAmM9eO8ElR1Ho6kf6KICI7n6+2KlFWGvUAb4GLna/4J8KjXMW864ykE3gAuc/YtHwvsUNWvnfvTgO9xJOsA1wHvqOoB5/5FOJLj8vykqpud5/xXVdMDuJ5PnecdAtYCy1S1wBn7cqBbAGUZY6IkqK9uRGQs8ByOlgnx2GV34BljXG5W1QVl7H8EaAGcGoK62gFHeXwdDo4Wu4YVLVBV94nIC8CHItJDVdd47nd+ff8N8LCqPuncNg9HC7invR6Pc3G0HIOjdRMcrYouGX6EdjFQX0Suc643BXqISAtV3eXH+d4ewdESvRj4UlXzRMT7mL+AIhwJrS/tcLT6lqCqC0XkPOBmHP2MZwH/ArZ5HJbp8XgXjueombPcZl6va1OP49vh6PvsWd/8UmL0tLf8Q0qV5fG4wGs9n79fX2NMJRZsH6Z/4+hHtlgdN0sAICJzgyzXGFMDiMjJOG6uGquq2z22NwVHv98Ai9wC5KqqZ//RugR/P0SB839fvzN7AE2Azzy21Qqg7DTn/y2Ajc7Hzco6QUTaAumqOs5jWy0cyeNo4OkA6nf5DEfL5kQc3SlKUNU9IvIFjt/7b3jFlIDj5r1HSom5EfCzqn4rIk1wdCV4i+IfgJri6HYA0BI4jONDwRZgkzpHwvCor7ZzdYvzeM/6BqrqwjKutzR5HuWC47U1xlRTwf5x2KKqP3kmwU5nBFmuMaaaE5H6OJKhj1X1Xa/do5xLoF4BThWRTs464nB8hT0wiDhjcSSXm3Ekit424UjYTnYe3w7HEHF+UdVtOLqN/NOjvsvLOgdHV4FiIzqoah7BdY9Q4FocNytml3HojcAAEbnStcGZhL+Ao0vAC6Wcdx5wg7OuPTiGffP+YDHWWV4cjm4P7zi70bwNdBaR45z7BXgRR/cQcAzddpnr5joRGQw85FHuAaCB8+a3x8q4NnCMENJbROKc9ZxfzvHGmKosmA7GOL7iGuZj+yfR7vxsiy22RG8BBgO/4ugmtRTHmLPexzzs3L/KeaznsgEYV0b543AkUoojiezuse9MHEOwzQMWANc7tw93npMLfFdKuW8C6fx9s9w8HF0FPgGOLCOe0TgS4u+AV3EkhGtwJH+TnWWm4hhi7GznvizgZef5rXGMSrEURzeL/zjjnOejrltxdAn4CWjisf1CYJ3zOfkBx01bwz1eh1+BwV7lrHEut/qo50JnPK5zT/XY1xpHi/CvzudoETAJqFXGc9Qdx4eS75yxzwV6eux3JeJf8fdQeA099h/nPGe+83WdxN9j4ccAdzvjmYuj20Vbj3NTgD+cr+Ug5+uQ6nxd5nnFGQ+874zhE+fz5H7P4Piwket8Lx0PvOx8LX293jdF+2fRFltsKXsJakINEdkEtAJy+Lt/lOC4Q7tuhQs2xhhTo4hNemGMiYJg+whn4WiZ8SQUH6zdGGOMMcaYSifYFuHTVfVbH9uPU9Vfg4rMGGNMtSciSTim3h6MYwi3a9Q5/J0xxoRbUIkwuMcAHYVj+JrNwCx1jKtojDHGGGNMpRVsi3BX4FscQ/1kePw/VFXXhSRCY4wxxhhjwiDYRPgLHHPRT1XVIucMQFfimFryzBDF6E8c8ThmfgLHoPRFZRxujDHGGGOqvxigufPxn6qa731AsInwPPUY4Ly87eEiIsfimOLUGGOMMcYYb31Udan3xmAn1KglInU8NzgHya9dyvHGGGOMMcZUCsEOn/Yx8KuIvItjas+WwGU4ps0MORHpD8wAHlTVqR67drseLFmyhKSkJL/LTE5OJi0trfwDg2T1WD1WT+Wupzpdi9Vj9Vg9Vk+k66mM15Kenk7fvn1dq7t9HROKUSNuAa4D2uOY730K8KwGW3DJes4D/oFjdqLnPBNhEWkNpAGkpaXRunXrQMolxKFaPVaP1VMF66lO12L1WD1Wj9UT6Xoq47Vs376d5ORk12qyqm73PibYFmFU9VngWa8gm+KYAjSUflPVT0RkXigLnTRpUiiLi3o9kVLdnrfqVk+kVKfnrTpdSyRVt+etutUTKdXteatu9URCVX3Ogm4R9lmoyPeqemrIC3aUPQ/HKBVTPbZVuEW4uonUJzJTudn7wNh7wIC9D4xDTX0fhKVFWEQ+AS5X1QMiUgR4P7PiY1vEeFww4PjkkJKSEvJ67rrrLjZs2EDdunVp1qwZxx9/PD/99BMbN25ERAIqq6ioiPvvv9+zH0uFVadPl6bi7H1g7D1gwN4HxqEmvQ9SUlKYPHmy38cH3CIsIhcD76uqisgvwCXehwAzVPX4gAr2v/55hKhFeMuWLbz55pvce++9TJkyhVNPPZWjjnIMR/z888/TpEkTLrroIuLiSn5eOOuss/joo484dOgQ27dv56effuKEE07gyCOPDPiapk+fjohw6aWXBnyuMcYYY4wpKSwtwqr6fx6rt6nqZu9jROS2QMuNtH379nHVVVdxzjnn0Lt3by6++GLuvPNOvvzyS9544w02bNgAwLBhw/jf//7HMcccU+x8ESEhIYGEhAQSExPp2bNnhWNJSEjgwIEDQV2PMcEqKiri0KFD1KtXD8D9NVqg33AYUx1kZ2eTnZ3tXi8oKGDLli1s3LiRDRs2sG3bNgoLCwFo1KgRnTp1olOnTnTu3JnExES/64mJiaFZs2bExAQ+mmleXh6pqanumDZu3EhRUZE7jk6dOtGxY0cSEhICLtuYmiLYm+U2icjJwAIgHrjLWebjwQYWTps2beLqq6/mkUceYeDAgdx0002ICC+88AIDBgxg2LBhPPHEE8TGxnLRRRdx++23U1RUxDXXXMOQIUNCnhgkJCSQkZER0jKNf1SVdevWUVhYWOaHmby8PA4dOkSjRo0iGF34FRYW8r///Y+ff/6Zw4cPk5OTw0MPPcSJJ57Is88+y9y5c/n000+jHaYxYXPw4MFi60uWLOGVV17hwIEDNG/e3L09NjaWtm3b0qlTJ4YPH07btm3d3xbu3buXjRs3snHjRr777jv279/vd/0FBQVkZGSgqsTFxZGcnEx8fHyZ8aanp6OqxMfH06FDB3fSe9pppxETE+NOir/99ltSU1M5fPgwAC1btuSss85i1KhRZdZhTE0SbCL8BLAH+Bm4D0c3iTSgMzA2yLLD4pdffuG+++7jjTfeoEOHDsDfLV433ngj//znP4t9em7RogXvvPMOW7duZeLEiRw+fJjhw4eHtNN5QkICubm5ISvPlE9VWb9+PXfeeSetW7dm48aNzJ49u8Qxr7/+Oh988AGqyqFDh7jzzjs588wzq80fkSlTphATE8O7775L3bp1ycrK4rzzzuOrr75i1qxZ5Z6flZXF//73PyZPnuzzA6KqWouyqXSys7OZMWMGH3zwAQ0bNizW/a1bt2489NBDtG/f3u/yWrRoQffu3YOOKz8/n+3bt7tbmn2pU6cOLVu2LLMFuVu3biW2qSrbt2/ngw8+YMSIERx//PFcc801tGvXLui4janKgk2EW6vqWAARGQecpqrrROS3oCPzIiJ9gSeB3sBEETlbVc8PpIz09HTuu+8+Pv7441Jb9kr7Cqlt27acd955pKenBxa4HxISEtyf2E1k3HrrreTm5jJ58mSOPfZYzjzzzBLHfPvttyxdupQvv/yS2NhYdu7cyfXXX09MTAyjRo2KQtShpap8+umnzJo1i1q1agGQmJjIjTfeyKhRo+jbty8rVqwos4y5c+fy4osvcvLJJ5OYmMidd95Jy5Ytueuuu/j++++ZOXMm/fv356mnnmLz5s3s2LGD448Py+0DxvglNTWVK664gquvvpqZM2dSp06d8k+KkPj4+IAS8ECICMnJyUyYMIFbb72VH374gSuvvJJGjRrx+OOP07lz57DUa/528OBBdu/eTUZGBrt373Y/Tk9PZ/v27ezdu9d9rIj4bEQoLCxEVTn66KM5+eSTOemkkwLqimNKCjYRjgMQkUHAVlVd59y+t/RTKkZVlwBDgjifm2++maeeeqrCX283btyYLVu2UFRUVKH+XKWxRDiyPvroI+rXr8+zz/49/HVCQgJ5eXnuhBBg3rx5XHXVVcTGxgKOrxVvv/12fvjhh2qRCC9ZsoRevXoVu2aACy64ABHhjDPOYNy4cRw+fLjUD4hr167llVde4d5776Vp06Z88skn7Nu3jyeeeIITTzyRBQsW8MILL3DnnXeyfPlyDh48yC+//BKJyzOmhJ07d3LllVcybdo09zeCNZGIsHr1ar7//nvOOeccXnzxRZ588slohxVRs2bN4pFHHqFFixb07t2bvn370q9fP1q2bOl3GarqTm49E9udO3eyY8cOtm/fXqyfeb169WjevDnNmjWjefPmNG/enGOOOYZhw4bRunVrEhMT/foGraioiD///JMffviBd955h6ysLLp3707//v3p2bMn3bt3d9/rYcoXbCK8V0ReAk4CngZ3y23QE3WE2q+//kqHDh3o1atXhcto0qQJe/fu5fDhw9SuXTtksVkiHFnvvvsub7/9drFtrtfW85fgtm3bSrTOJCUlsWvXrojEGW6vv/46N910U4ntIsIFF1wA/P28lDZteXp6OhdeeCEfffQRzZo1o27dujRu3Jjnn3/efcxNN93EjBkzuOiii3jsscfKTKyNCZesrCzGjBnDCy+8UKOTYHD8bvvwww/Zvn0711xzDUuXLq1R3Zh2797Nk08+ybfffktRURHLly9nyZIlfPDBB+zevduvhi5X98j69esXS26bNWvGgAEDaN26Na1ataJBgwYhjz8mJoZevXrRq1cvbrrpJve9LkuWLOGTTz5hzZo1HDx4kPj4eLp06cIJJ5zg/p1uSgo2Yb0Kxw1yb6vq6yKSBNwEPBZ0ZCG2Z88eOnXqFFQZjRs3Zs+ePeTm5loiXIUdPny4xKflxMTEEolwZmYmTZs2LXZckyZN2LNnT0TiDJeioiL2799PampqucP9NWnShMzMzDIT4aSkJOrXr19mOa6hAdu1a8fmzZt99mE0JlxycnIYPXo0jzzySFAj/FQX8+fPZ+zYsbRq1Yp+/fqxatUqVq5c6R4+tDpTVW655RYef/xx99/x448/vkp32RIRunfvXqKfen5+Pn/99RcpKSl07do1qIbA6iyo7/dVdY+q3quqjznX01X1SlX9OjThhc6hQ4eC7gvmah07dOhQSBPh2rVrWyIcIZ6tHp988gmjR48GoGHDhj6HsPNuIWnUqBH79u0Lf6BhdOONNzJkyBAmTpxY7rGuD3+lOXjwYLlJsKfWrVuzY8cOv483JhRuv/127rzzTvr37x/tUCqFGTNmMHLkSAAGDx5Mo0aNmDlzZpSjiox3332XI444IiQTWFV28fHxHHHEETz88MM8+uij0Q6n0goqERaRESLyhojUFpEjRGSDiKSKSKX7aBWKRLhBgwbs37+f3NzckN5gYS3CkZOdne3+qmr58uV8//33gOO19Wcs59jY2DLv6K4KtmzZwrJlyxg8eHC5xzZs2LBYH7dgNW3alMzMzJCVZ0x5CgsLSU1N5dRTT412KJXC77//TteuXWnWrBkAJ5xwAlu2bOHHH3+McmTht3XrVt566y3uvffeaIcSUZ07d6Z27dqsXLky2qFUSsHe8XUb8Jyq5gIPAHOAyThGd6hUQpEIu1oHrWtEcNasWcMbb7zBDTfcwLXXXsuiRYsiNgf6nj17aNy4MQAHDhygfv36bN68mczMTP744w9WrlxJXl5eRGKJhkCf53r16lkibKq0hQsXMnDgwGiHUWmsXLmy2PORkJBAfHw8LVq0YPv2EpNuVRtFRUXceOONPPvss9Vm+MtA3HvvvTzyyCPRDqNSCraPsKjqUhFpgGNEh7aqesg5lFqlEopE2CUciXB1Gkc4KyuLn376iSVLlrBx40aOPPJI+vXrx+bNm3nvvffo3LkzQ4YM4a677iI/P59XXnmF++67j+effz6sfUdVlbVr17Jz507+/e9/8/7775Odnc3DDz/Mzp07iYmJYdu2baxbt468vDyWL1/OP/7xD2JiYujevTv9+vWjW7duHDx4kPXr15OcnEzdunXDFm847Ny5M6C7ouvXrx/Q5ADladq0KatXrw5ZecaU58svv6wWo7yESmpqKkOHDi22LT4+nhEjRvDZZ59x3XXXRSmy8HrxxRc5/fTT6dGjR7RDiYquXbsSFxfH6tWrOeKII6IdTqUSbCLcSERqAWOAb1T1kHN7pfvuOJSJcKj7CFeXFuHCwkJeeeUVPv74Y84++2yGDh1Kx44dWbVqFYsWLaJ58+Z89NFHJW5Ue+KJJ9i2bRuXX345H3/8cUjHRHSNlTt16lQOHz5MfHw8sbGxnH/++SxduhQR4eWXX+bjjz8mKyuLq666yn3eyJEjeeuttygoKGD16tUsXryYefPmsXnzZl544QU2b97MoUOHaNasmXsIshYtWtCvXz969+7t7jtbv379gPrRhtOePXtK3ABYlnr16oW0T6+1CJtI++2337j//vujHUalsXnz5hKj4TRv3py+ffty9913V8tEeO3atXzzzTc1ph90ae69914eeugh3nrrrWiHUqkEmwh/AOwEagGDRKQh8BCO2eUqlVAlwiJCTk5OSBPhuLg4CgoKQlZeJB08eJClS5eyePFiZs2axXnnncfXX3/tHnsXHEOOldc/r02bNjzyyCNceeWVfPDBB8VmeqqoDRs2cNddd3HUUUcxffp06tWrxzfffMO6devo27cvMTEx7u4udevWLfa1oIgQExPjfs8MGDCAAQMGALBu3Tqefvpp97F79+51v37bt29n8eLFPPvssxw65PhcmJWVRXZ2Ns2aNePRRx8lOTk56GurqAMHDgQ0nE/9+vVL7RqRn58f8OtU3s13xoRSWloaSUlJIR33vapLT08v8a1QUlIS2dnZiAj79++nYcOGUYou9AoKCrjlllt49dVXa/z7wDWqxNq1a0MyE2J1EVS2oaqPi8hnwEFV3SYidYAPgUr33WeoEuGGDRuyd+/ekN4sV1XHbvz++++5//77Of300+nXrx/jxo0LqjV34MCBXHDBBUyePJkHHnggqNgWLFjA/fffzwsvvEDXrl3d23Nyctwt0iLi7jNbp04dd+IKjtZtz2S+LK4+x4B7gHRf/vjjD66//no+/fTTqP1CrkgifPDgQQ4dOsSyZctYvHgx+fn59OvXj06dOpVo6b7ttts444wzSE9PZ968eUybNq3Yfn9vSjQmFL766iufs0bWZKpa4vdPq1at2L59OyNGjODrr7/moosuilJ0ofe///2P0aNH21TSTq6+wlOnTo12KJVG0M1uqrrW4/EhYL6IPI5jfOFKI1SJcJMmTcjIyKhyfUNDSVV59tlnWbhwIbNmzQrpDDaXXXYZZ511VlBl/Pzzzzz00EN89NFHJZK+nJwc6tati6oWG0qtbt26xRLh8t4vFRl8vlevXgwfPpynnnqKO+64I6BzQ2HXrl3Mnz+fNWvWuK81PT2dDRs2oKo+E//s7Gw2b97MqlWrOOaYY+jbty9xcXEsWLCA//f//h8rVqzg8ssvp1u3bsTExDBt2jS2b99OYmIiv//+e4ny6tatS05OTtiv1RiA2bNn8/LLL0c7jEqjtFlRXYnwqFGjuPvuu6tVIjx37ly+/PLLaIdRaRxxxBEUFhayfv36Yo1ENVnAibCInA98oqoqIm+Ucthwqmki7Ppqt0mTJiGIqmo5dOgQH374Ie+++y6nnHIK7777bqVrzV62bBn3338/77//vs+Wz5ycHFq0aMG+ffvcYwKrKnXq1CmWoOXk5JT6fnGN+1yR7jHXX389F154IaeffnqpLceh8tdff/HBBx/w+++/k52dTYsWLSgsLKRnz56ccsopADRr1ozOnTuX2kK9Y8cO/vvf//Lqq68W2z5s2DDWrFnD66+/zt133+1OpmfOnMmSJUvYv38/jRo1Iisrq9i3BJXt/WKqr8OHD3Pw4MFi39jUdDt27KBVq1Yltrdu3Zrly5fTqlUrMjIyyM/PrxYjK7i6MdrvneImTpzII488whtvlJbC1SwVaRG+CscwaQeAU4CpPo6pdEMghLJFeMuWLSHtI1zZFRUV8frrrzNjxgwuvfRSPvzww0pz85e3Tz/9lAceeKDUPm6uFuHt27fTunVrDh06RH5+PrVr1y42csehQ4dKbfV3tR5X5D0gIrz44otcd911fPLJJwGf7485c+bw/PPP06RJE8aOHcvNN9/sfr2ef/55Onfu7PcsSvXq1ePgwYM+92VnZ7unF23evDkFBQUkJSXRtGlTCgsL2b17N2PHjqV58+acf/759O/fn5YtW9ofJRMRCxcu5IQTToh2GJVKamqqz+mlW7Vq5b4p9qSTTuLnn3/2a5zxym7x4sU2iYoPRx55JLm5uWzYsIHOnTtHO5yoCzgRVlXP762nqepk72OkEv6lC2WL8IoVK0LaR7gyW7NmDbfeeitnnXUWs2fPDslNbOVx3aBVkWTb1x3RnlwJ7O7du2nevDnbtm3j8OHD1KlTp0QiXNpr7Go9rmhLU8uWLcnPz6/QuWUpKiri/vvvJz09nalTp/qML9A+wmUlwq5xmF127NhB69atUVXmz59P586d+fzzz1m3bh2zZ8/mww8/ZOPGjZX2Q5SpXpYvX86xxx4b7TAqlc2bN/uVCP/yyy/VIhH++eef7cNQKVx9hV977bVohxJ1wU6xnFLKrl3BlBsOeXl5Ifmqp0mTJmRlZdWYFuGnnnqKp59+mptvvjkiSTA4EsX09PQKnetKcEvjmhUwIyOD5s2bu4eu824RdrUc+xKqfq6hnEQkJyeHSy65hFatWjFlypRSk/QDBw4EdEd4WTPpec7SB5CRkUGzZs0YMWJEsXFbu3Xrxk033cTUqVMZO3ZshV9bYwLxxx9/0KtXr2iHUamkpqb6bCiIi4tz/5z37duXxYsXRzq0sPjtt9+sRbgURx99NAcPHmTjxo3RDiXqgp1iOUZEzhORe0Tkv64FuDNE8YVUKBqqGzduzP79+2tMIrxly5aID0DesmVLdu7cWeHzy3qdXS3CmZmZNG3a1J0A165d2++b5bxHmKiIxMRE9u3bF1QZnr7++mtOOukkrr322jKPC7RFuCzerfauMYrPOeccHnvsMZ/nXHzxxaSlVbrRFU01lJaWFtWhCiujrVu30rZt2zKPcY0UE6nZPsNFVTl48GBIb+aubiZMmMCUKVOiHUbUBTuG0wzgfuBooKPHUm2zxCZNmpCdnV0jEuGCggJiY2Mj3qezoolwUVFRubG6WoRdLb6uFmHvSU3CnQgnJyeHNCH8/fffOfHEE8s9bv/+/WFLhDMzM8u9ibRhw4YkJCTw119/hSQGY3xx/S6ohL30oio9PZ2kpCSf+zyHk2zfvj1btmyJZGght2HDBrp06RLtMCq1AQMGsGTJkio7j0GoBJsIdweOUdXLVPVK1wK8EoLYKqXGjRuTnZ1dI/oIR2vQ7Yomwrt376ZFixZlHuNKhF3Jr2sECM+vBqHsrhGVMRH+888/Oeqoo8o9riItwqUlE75ahP0ZTaV79+4lRqEwJpRKuymspsvLy3PPgunN836A4447jl9//TWSoYWc9Q8un4gwbNgwZs+eHe1QoirYRHhlKduXBllupdW4ceOQzyxXWS1btiwqN5tUNBHesmVLuYOmu7pBuBLhhISEYn2DXapSi7Cquq+nPHl5eX4d5w/vm+Vc3U3K07lzZ3788UeKiopCEocx3v744w+OPvroaIdRpSQmJpKVlQVYIlyTXHbZZbz99tvRDiOqgk2E7wFeE5FbRORy1wI8EoLYQipUX5HVrl3bPdxWKMXGxla6ryeWLl1apRLhzZs3+5UIe3aH8OwS4dknzp9RI4IRykQ4LS2NNm3ahKQsX0rrK1jRFuGGDRvSp08f5s6dG7IYjfG0YsUKu1HOS2mTabgkJiayd+9ewHGD69q1a0s9tioo7cZAU1zr1q3Jzc1lz5490Q4laoJNhO8CzgL+AVzpsfjuhFRNFBYWhjwR9u6jWhmsWbMmKl0jWrRowa5dgQ88UpEWYc/RIjw/LJU3akRlahFesmQJffv2DUlZvsTFxfkc7q2iiXCDBg0YNGgQn376aSjDNMbN365CNcnevXvL/Pn0vIE3JiaG+Pj4Svc3yV+ueyGsj7h/Lr30Ut57771ohxE1wSbCQ4F2qjpIVU9xLcD7IYit0gpHIuw9fFe0qSqFhYURGzLNkytRDdSWLVvKbQEoLCwkNjaWw4cPU6tWrVI/gIS7a0RSUlLIhhELdyLsuovcm/fwaXv27PFrbOUGDRrQqFEjNm3aFNI4jXHZv39/QMME1gRZWVll/ny6Ztp0OfbYY1m2bFkEIgu9RYsWMXDgwGiHUWWcffbZzJw5M9phRE2wifAaVfWVvb0ZZLkhFephYFQ15NNPVrYW4a1bt5bbulrZbN261e8uAr66RngK981y8fHxIesKE+7+kK4JTrxlZ2cXG5rI+4NTXFycz2ts0KABOTk51kfYhEVFZ32s7rynO/fm2UcYYODAgVW2n7D1Dw5M7dq16dy5MytXlnbbV/UWbCI8X0TeEJGzReRk1wK8FIrgQqWsO2UrKtRfuVS2RDha/YOD4ZoYoyyuD0Wum8Y8W+ID6SMcbCIcSv5cNziuryLv29JahMv7ufIem9mlQYMG7ok9QjmWsjEAq1evpmfPntEOo9LJysqiUaNGpe737CMMjqG1Fi5cGInQQm7JkiX06dMn2mFUKePGjWPatGnRDiMqgk2EHwVOAZ4BpnkskZ2BoRyhml45nCpjIty7d++o1R+KG9J8cSWCvlqEPZPESCTCpfW9DcT27dtp3bq1X8fm5ORU6OegtBbh8pJq72mrXRo1asT+/fs54ogjWL16dcDxGFOWFStW2IgRPuzdu7fcFmHPD6ZNmzYtlhhXFUVFRWG5ob2669+/f40dUzjYRPgnVe3ovQCValC6cCTCoe5uUdkS4WjfdR3oyBGBDmnnPY4wFH9Nw901AqBVq1bs2LEjqDICafmo6Kxy9erV85kIe/L181Bai7CrL6IlwiYcLBH2rbzp533NdtmiRYugZvmMhtWrV3PEEUdEO4wqpyaPKRxUIqyqp5ey/dxgyg21UCfCMTExIU9aK1sifOjQoVITwUhISkoK6Bewv32aXa2Y+fn5xMfHFxtHOCYmxj2pRiRahEMxcsTvv//u941yFU2ES2sR9i7b++akslqELRE24bJ27Vq6desW7TAqHX8SYc8+wuAYT7iqdY/45ZdfOP7446MdRpVUU8cUDioRFpEBIvJfEaklIski8q2IzBORyI+5VYZQJ8JxcXE++0wGozIlwrt27fJrGKxwCrRF2J8RI6B4y6WIFHvePRO3cI8jDI7xG4NNhJcvX84xxxzj17H79u0rs49gaUpLhD2fS19Dp5X2gcGVCHfr1o3169cHHI8xZSkoKAj5zczVgT+JsHdXiKo4sYbdKFdxNXVM4WC7RkwCVgD5wEPAQWA+8GKQ5YZUOFqEq3MiPHPmTM4666yoxlCRRDiQFmEXz5vlPMcHdo037EsoxhGG0LQIB/Le3rdvX5l9BEvj62Y5764QFUmEK9N73lQPGRkZNGvWLNphVErlJcK+fl6PPvpo/vjjj3CHFlKB3DdhSqqJYwoHmwjXUtVPgFo4Jta4QlUnAaGZwzVEQj2cTrhahCvLOMJffPEFI0eOjGoM4UqEXVyJnHeLsGdLb2k3g1WWrhHp6em0bNnS7+ND2SLsnYD7ml65tOepbt267ue5Mr3vTdVn/YNL5z0Bjj/i4uLcY8pXBZmZmVH/NrOqq4ljCgebCLt+qs4FFqhqlnO98owtRehbhGNjY0OeCFd0EolQy8zMpFatWhXqSxpK4UiECwoKiI2NBf5Ocj0TYc8ErSzx8fFBj/YAwSfCgfQPhoonwr5ulvPubxxIi7DnB4yuXbta9wgTMpYIl628kV583fR63HHH8f3334crpJD69ddfrX9wkGrimMLBJsILRGQV8ArwrIjUEZE7gLLvrImwUCbCqlqt+wjPnDmTc889N9phBJwI7969u9yvRH11d/DuGuFKhMsaFSRUY0h7z+QUqEDHygymRfjAgQPFtnknwr5aYvxpObcb5kwo/fHHH1Ed7aYy82ekI19/h6699lpefLFS9XYsld0oFxo1bUzhYEeNuAu4DRiqqt8D8UAGcHcIYguZUCbCrmSquibCn3/+OaNGjYp2GBXqflBeguprSDTP592f0RFCKdiEetmyZQGN9VzegPql8TWskvcoEYF0jfBkibAJJesf6tu+ffv8mnK6SZMmJW6Ya9GiBc2bN68SLYRLly71++ZhU7r+/fuzcuVKfvnll2iHEhHBtgijqt+o6iLn4/2qOg04NejIQig3N9cSYT/s3buX2NhYv35hViZFRUV+JZW+EmHPLin+jJfrEspxpCta1sGDB4tNcVyeirYIN23alMzMzBJlWSJsKhPX74FQz/pZHezevduv+wmaNGlS4mcd4LbbbuOpp54KR2ghs379epo3b24jhoSAiDB9+nTuv/9+5s2bF+1wwi4umJNFJBY4B+hK8RvkxgEvB1N2KB06dKjMu2UD4Uqqq2MiPHPmTM4555yoxlARf/75Jz16lD+ZYWktwq6uEaVNJRxOrrE7GzduHNB55d0B7ktFE+G6deuWeF4yMjKK1e8rES5tQg1w9LPPz893zzJnTLDWrl1L165dox1GpZSVleXXiDGtW7dmx44dHHnkkcW2H3HEEezZs4f09HSSkpLCFGVwnn76aW677bZoh1FtNGrUiA8//JDRo0dz6NAhRowYEe2QwibYFuHpOIZN6w109Fgq1dyGoewakZGRQZMmTaplIvzpp59Wim4RLv6OKDB//nyGDBlS7nGeLai+Ro1o0KABBw4cKHZTXbhV9Ia5JUuWBHSjHPj/x9CbrxY272GqfA3NVtqEGgDNmzdn9+7dgCMprip3pZvK66uvvmLYsGHRDqNS8rdbVPv27dm0aZPPfTfffDPPP/98qEMLiYyMDLZu3WrdIkKsXr16/N///R+vv/46H3/8cbTDCZtgE+EjgKNU9VJVvdK1AC+FILaQCWUivHXrVpKTk6tdIrxhwwYSExMrlCiFi783zP3444+cdNJJ5R7na/iguLg49wgQrkT44MGDERs1I5hEOJAb5QDy8vJISAjNyIbeibCqEhNT/NdJWV0jkpKSSE9PB6BTp06sXbs2JHGZmmvu3Ll+fSCuifwdQ7xDhw5s3rzZ574hQ4awaNGiiH9r5o+XXnqJ66+/PtphVEu1a9dmxowZfPzxx7zzzjvRDicsgk2E16uqr6acSjVZdagT4bZt24ZkZjFP0R5P9eWXX+a6666LWv2++JMIFxUVsX//fr9aO7Kzs6lXr16xPsWefQpdozgcOHCg3PE2PadjDsaAAQMq1Adr2bJlHHvssQGdE2y/Zs/z/Zm4oKxRMVq1auVOhM877zw++OCDoGIzNdv+/fuJj48P6TCZ1Ul6ejotWrQo97iyEmER4Yorrqh0ownk5uYyb948hg8fHu1Qqq34+HimTZvGvHnzeOWVV6IdTsgFmwj/ICLTROQcETnZtVDJWoSzs7NL9A2tqK1bt9KhQ4dq1SKcm5vLsmXLGDhwYFTqL01ycnKpv5RdVq1aVaI/W2lcfWS9W0ZdCZ5rilF/Bp73NZJCRQwcOJBFixZRUFAQ0Hnew5eVJ9gk2Lv/tPcfVl/l+7rJziUpKYkdO3YAcPLJJ/Pjjz+G9AZEU7N8++23DB06NNphVFqbN2+mQ4cO5R6XmJhY5vS6F110ER988EGl6sr07rvvMmbMGLtJMsxiY2N55ZVXWLFiRaW/cTJQwSbCjwInA08D0zyW8u9ciqDU1FTat28fkrK2bdtGp06dqtWEGh988AEXXXRRpftFcuaZZ5bbL2n+/PkMHjzYr/JcLce5ubk+uwg0adKEPXv2+HVTWVlJXiBEhGHDhjF7tv9folRk9qQ9e/aUuJktEO3atSv2ocRz+DRfNyHC3zcC+uLZIhwTE0OfPn1YsmRJheMzNduXX35ZrW/mCZa/fwNdfwNK+1AaHx/PmWeeyWeffRbS+CqqqKiI6dOnM2bMmGiHUiPExMTw7LPPkp6ezoMPPlhtGi+CTYR/UtWO3guVrGtEqPtGtm3btlq1CE+fPp1LL700KnWXpXXr1uTl5ZWZcP74448MGjTIr/Jcs59lZ2cXG/rL9cu/cePG7N27l127dpX7NWJpwwxVxNixY3n77bf9Pj7QGeXAMb5qq1atAg3NrUuXLvz1118+97m6C3mLjY2lqKjI5znt2rUrdlPO6NGjeffddyscn6m5VNXvFs+a6uDBg35Pr9ymTRu2bdtW6v7x48fz6quvhiq0oHz99deccsopIfv7bsonIjz66KOoKv/617+qRTIc7IQap5ey69Zgyg2lgwcPhqxbhEs4htlq3LgxGRkZQZeTkpIS0PE///wzXbt2DWg82kgaPXo006dP97lPVQMaeszVr3X//v0+x0p2jV7gTyIcqhZhcLSO5uXl+fX6qyqvv/56uXfHe78P1q5dS7du3SocY9euXVmzZg1QsqvRli1bfCbCrnh9ad26Ndu3b3evH3PMMaxYsaJSfeVa1QX6u6CqWrp0acD95WuSSZMmBZSsnHTSSXz55Zel7m/UqBFHHHEEn3/+eSjCC8qUKVPsJjk/hfL3gYjwn//8h+bNm3PrrbeW2uBRVQQ9oQaAiDQVkXauBZgRinK96qgtIlNF5FcRWSwiZ/hz3po1azjiiCNCEoOqoqphab113WEfbCIwefJkv4/dsWMH//73vyv1H8xRo0aV+gt3zZo1fo0f7OJqFS0tEQbHa5yWllZu66l3V4FgjRkzptSE39PUqVM56qijOOqoo8o8zvt9EEhfal/69evHokWLgJJDt61evZru3bv7PM81XrA3769gRYRTTjmFuXPnVjhGU1wgvwuqsi+//JIzzzwz2mFUWvfff79fk2m4XHLJJeV2Sbv//vt57733ePrpp6PWIrhs2TLatm0bVJevmiQcvw9uv/12jjzySK699toq3YgRVCIsIseKyHpgF5DqsYTjrqsUQFT1OGA08J6IlPvTvWrVKnr27Fl6oQEkga5+liJCly5d/D7P33qOOeYY/vjjj4DKrajc3FzGjRvHCy+8EHB/U4hca9PDDz/MkUceydKlS0vsmzdvXkDDJblGjfCeFU1V3deTmJjIihUrSk3sXHr06OFuIQ1Eac/bWWedxccff1zmjSobNmzg/fff59577w243iVLltCrV6+Az3Np0KABOTk55Ofn8/3337u7o6SkpPD777+XOtVz165dWb9+vc99Xbp0KTZs2ujRo3nhhRfIy8srcWwk3m+Rek9X5g+eFRGt5y03N5cnn3ySH374gRNPPDFs9YRLJN8HgfShjY+Pd3cTK029evV45513yM3N5Z///Ce5ubkRf97CPYFGdXwfhMO1117LySefzLhx4/jPf/4TkTpD/ZxJMJ/mROQb4DrgDVU9RURqAaOAo1U1JTQhgojE4Ei2z1fVH5zbvgO+VNUnRaQDsAkcf/A9Z755+OGHOeecc0ptDQtkHNc///yTL774gokTJwZ8Df7U8/3337NhwwauueaagMsPpJ5du3aRkpLCBRdcwGmnnRa2ekIhOTmZOXPmMH36dB588EH39rlz5/LMM88wdepUv8bHPHz4MOPHj2fatGm88847NG/e3N294Oqrr+brr78mLS2Nl156iU8//ZRvvvmmzPIKCgq4/PLL/WrF9b6e0p63xYsX8/jjj3Puuedy8cUXu8fl3bNnD3/88QdTpkzh8ccfL7UbQmn15OTkcMMNNzB16tSAYvX2+uuvk5iYyIwZM3j//feJiYkhOTmZQYMG8d577/k8Z/bs2WzcuNHn0Hzz5s1j+fLl3Hrr3z2pZs6cyTvvvMMtt9xSrO93JN5vkXxPWz0Vr6ewsJCPPvqI999/n8suu4yzzz67xBjWoagn3CJRz5o1azjttNPYtm1bQDdDz5gxg40bNzJy5MhSP+S6fP/997z66qv88MMPLFq0iNatW4ftxmtVpU2bNpx33nk0a9YsrElkdXofRKKeWbNmMX78eDZs2EDt2uGdUy2Qa0lPT/f8BrOjqqZ6HxNsIjxbVc8QkXmqOsRj+0xVDdlcvSLSBVgPtFTVXc5tzwHNVPVSEbmEMHTHMMYYY4wx1cKlqlqi1SbYj9Fxzv+LRKQvgIh0Ayr+Haxvri4QngO3ZgHljxBujDHGGGOMD3HlH1Km30TkDuBl4FcRSceRtD4dbGCl8G6+dn3/8qtrw6xZs4p1jShP3759IzJ+qdVj9Vg9lbue6nQtVo/VY/VYPZGupzJeS3p6OiNHjnSt/urrmKC6RhQrSORE4Hhgjap+EZJC/y67K7COkl0jmqvqJSLSGkgDSEtLo3Xr1oGUHZG7Xq0eq8fqqdz1VKdrsXqsHqvH6ol0PZXxWrZv305ycrJrNVlVt3sfE+yoERtF5EUAVf1JVZ8IdRLstAHYA3jeyt8T+M37wOTkZETE7070kyZNCkV8laaeSKluz1t1qydSqtPzVp2uJZKq2/NW3eqJlOr2vFW3eiKhsj1nKSkpnklwqYK9WW6ZqvaucAGB1fUojhbhK50txL8CR6pqejAtwtVNpD6RmcrN3gfG3gMG7H1gHGrq+yDsLcLAWhGJ994oIv8LslxfUhxFy684Roi4VFXTw1BPlVadPl2airP3gbH3gAF7HxgHex+ULtgW4fuAocDHgOfo2xNVtfRZLELMWoSNMcYYY4wnf1qEgx014g5gGXCu13b/h20wxhhjjDEmCgJOhEXkv86HM4BZqnqZj2OmBBuYMcYYY4wx4VSRPsLnAfOATF9JMICqlpxPNUICHTXCGGOMMcZUL2EbNUJE5qrqKeUc01VV1wdUcNnlTQU6eG0+S1Wznfutj7AxxhhjjHELVx9hfzLnl4FTK1B26ZWqDgllecYYY4wxpmarSCI8REQKQx6JMcYYY4wxEVSRPsLLcbT2lracBvwRqgBdRGSKiPwoIl+JSJldM4wxVUuPHj0YMmQIQ4YMISkpiZYtW7rXe/ToEZWY8vPzefrppxkyZAiDBw/m+OOP57vvvnPv37lzJxdeeCGDBg3i6KOP5s033yx2fm5uLuPGjeO4446jX79+zJ49O6D9LnPmzKF3796ICIMHDyY9PZ1NmzbRs2dPWrVqxR133FHutdx///0kJSVV6Xsn1q9fT3x8PPPmzSv1mHHjxrnfN64lOzu73LK9n+OTTz6Zrl27MnbsWA4ePBjCqyj/fVXRY8NxHZ988om7vOnTp5fYf+DAARo1akT79u1DPk5tWa9leT97kbqGcPxceb+Ge/bsCenxkRDo+9afn9vffvuNLl26MHXq1NAHrKoBLcD3fhxzRaDlllPePUB/5+P+QDbQ22N/axxdNjQtLU2NMVXL4MGD3Y+vuOIKHTNmjM99kbRp0ybt0KGDZmVlqarq7NmztV69erpt2zZVVT3llFM0JSVFVVU3b96sDRo00F9++cV9/j333KOXX365qqquXbtWGzdurOnp6X7v9zR37lwFND8/X1VVt27dqqeffrqmpqb6fT1XXHGFTpo0ye/jK5tLLrlE69atq3Pnzi31mCuuuKLC5Xs/x3v27NEWLVqE/Dkr731V0WPDdR1z587VOnXq6MCBA0vse+6557ROnTr673//u0Jll6Ws17K8nz1v4byGcPxceb+GvqSlpWlycrLfx0dSoO/b8n5uP/74Y7300ku1T58++uabbwYUS1pamrryQ6C1+sgxg51ZrrTkelqIy/t/qvqb8/FvwCzgWl/HukaNcC1VuQXEmJrikUceqdC+cGrQoAH3338/jRo1AmDo0KHUrl2bn3/+mW3btjF37lyuuuoqANq1a8fQoUN5/fXXASgqKuK1117j6quvBqBbt24ce+yxvPPOO37tL8vWrVu58soree2112jfvn3Ir7sy+u2336hfvz7NmzePWJ2NGzdm0KBBLF68OKTllvW+CubY0oTiOi655BIWL17Mb7/95t6mqsyZM4f+/ftXuNyKKO9nrzSV6RpC4csvv2T48OHRDsOnULxvPfXv35/p06fToEEDv45PSUlx54D+jBpRkUS4lYh8LyIDK3BuqGwBfP4FSEtLK5bpWyJsyjJ94ZZoh2CA448/vkL7wqlp06aMHTvWva6q5OXl0bx5c9LS0gBo1qyZe3/Lli1ZtmwZABs3biQzM7NYt46ePXu6k5Hy9pdm69atXHXVVbz++utBJcHr169n+PDhnHzyyZxwwgl89dVXxfYvXryYfv36ceKJJ3L99ddz0kkn0aNHDz777LOAjgmVSZMm+f219XXXXcegQYMYMWIEc+fODaregoIC2rRpE1QZ3sp6XwVzbFmCvY527dpxzjnn8Mwzz7i3zZ49m6FDhyIiFS63PL5ey/J+9koTqWsIxc+WP7788kvOPPPMkMUdShV535b1cxvoezclJcWdA7reL2UJ+GY5VT0i0HOCJSJ3q+pjHptaAiWGwDAmUDMWbWH0wHbRDsOE2YQJE0r9Q9m7d2+efvrpcsuYP38+7du35+STT2bHjh2A449yly5dAEhPT2fvXsdM8zt37gRwt4gAJCYmsmrVKr/2l+aUU07h3XffpV27ir9nCwoKGDVqFBMnTmTcuHH89ddf9OnTh6VLl9K5c2fy8vI477zzeOyxx7j00ktZtmwZ/fr147XXXuPss88G8OsYbxV9Db766iuOPPJIv/4YHnHEEZx66qn079+f3377jVNOOYUFCxbQu3dvf56aYrZs2YKq8p///KfUY0L9vgrlsS7+XIc/brnlFs444wyeeOIJkpKSeOutt3j55Zf58MMPyzyvos9Raa+l631Q2s9eOK7BX6H42fJHfn4+P/zwA9OmhfTLdyA072lv5b1vQ/lzWxHBTrEcKbeLyFRV3SUiHYFzAP/fNcaYGq0iv7w95ebm8q9//YupU6cSExNDcnIyp59+Ok8//TTPPfcca9as4bvvvqNt27bFzvNuaVKvcdvL2++tffv2XHvttSxYsIDExMQKXcvChQvZuHEjl13mmA+pS5cuDBw4kHfffZf//ve//PLLL+zatYuLLroIcPzx69mzZ7Ey/DnGW0Veg6KiIv7f//t/fPzxx34df88997gf9+/fn5EjR/Lyyy/z0ksv+V3naaedRk5ODitXruSxxx4r86vVUL+vQnUsBHYd/hg8eDBHHHEEU6ZMYezYsSQlJVG/fv1yz6voc1TWa+nPz14or8FfofjZ8scPP/zAMccc43dXgUAE+5725s/7NhQ/t8GoKonwE8AnIlIA1ANuUtX5UY7JmCovnF9rupSX3AXjjjvuoHv37mRmZpKRkeF+fO+994asDlVl/Pjx3HbbbfTt29e9/b333uOuu+7ixBNPpF27dvzzn/9kzZo1ALRo0QKArKwsn4/L21+ajz76iMGDB3PeeefxzTffUKtWrYCvZ9u2bTRu3Ji4uL9//Tdv3pxt27YBsGPHDhITE4mNjXXvb9KkSbEy/DkmFKZPn86wYcMqXHa7du1YuXJlQOd89913xMXFcc8993DXXXdx0UUXlfu6VERp76tgj3UJx3XcfPPN3HfffWRmZnLrrbcGVVagPF/Lsn72yhPOawjFz5Y/Klu3iKlTp7pHcxg3bhzjxo0DKva+hYr93AajSiTCqvoEjmTYGBNC4UxSI2HZsmU8+eSTgKMFzPXYWzBf991xxx0MGDCAf/zjHxw+fJidO3fSrl07mjZtyhtvvOE+7uqrr+bUUx3zCHXu3JkmTZqwdu1ad/KxatUq9x+v8vaXJjExkS+//JLjjjuOcePG8e677wb8YaZt27bs3buXgoIC9x/s3bt3u/srt2rViqysrGL7MzMzi5XhzzHeKvIa/Pjjj/z555988803gOMr8AkTJtC+fXtmzpxZ4vjHHnuMu+++272+c+fOCs80OmnSJN58801efvnlUrsVhON9Feyx5V3HL7/8wp133skJJ5zAo48+WiwpK8+YMWO45557SE1NdXdLKE9Fn6OyXsuyfvaCvYZgnp9Q/GyVZt++fSxYsICRI0fy5ZdflvstiefxGzZs4IknnmDfvn0MGzaMK664otTzKvJ6eSa/nvx934by57ZCfA0lUdUWPIZPcy1VeZggEzlnPftjtEMwXryHT1NV3bFjh9599936zjvv6E033aTbt2/Xr7/+Wrt166Zvvvlmsceh9Oijj+rNN9+sBw4c0AMHDujKlSvdv1suu+wy3b17t6qqrly5Urt16+YeLkjVMTzauHHjVFV13bp12qRJE92xY4ff+z15D4+0ePFirVevnt51111+X4trmKeCggI94ogjdNq0aaqqumHDBm3QoIH+9ddfqqp6+PBhbdOmjb7zzjuqqrp06VJNSEgo9tz6c0w4tG/fvtjwabfeeqtOnz7dvd6yZUvduXOnqqpu3LhRGzZsqPPmzfOrbF9DUD344IPaqlUrzcvLC80FOJX1vlItfl3lHVuR67jiiivKHIbOuzzP+mbOnKnLli1zrw8ePDgsw6eV9VqW97MX7DUE8vx4Dp8Wip8tV7zer+GmTZv0iiuu0I0bN2rHjh39Pj4tLU0vvvhi3b59u+bn5+ubb76pDz/8sF/XFoyy3rcV/bkdPHhwQL9jJk2apF65oc/h06KexIZiwcYRNhV0wiPf+XXcu79uDnMkRlX1rrvu0pYtW2qLFi2KJXkTJ07UCRMm6JtvvqmTJk3Szz//XFWLjzEc6vGG165d6/1LtNiH7IkTJ2rPnj31xBNP1GHDhunKlSuLnX/o0CG94oordODAgdq3b1/95ptvAtrvMnv2bD3mmGMU0JNPPtmdLH/xxRcaGxurJ5xwgq5fv17vvfdevfnmm32WMXnyZG3ZsqW2b99eX3vtNf3rr790+PDhOmjQID3++OP1yy+/LHb8okWLtE+fPnriiSfqhAkTdNCgQTp16tSAjwmVhQsX6uDBgzUhIUGPOeYYnTx5sqqqnnXWWfrkk0+6j3v88cf1hBNO0JNPPln79u2rb731VrFySnuOvJ9j12uZlZWlDRs21L59++qPP4bmQ3N57yvP6/Ln2Ipch2eiV9b7xlVe+/bt9bHHHiuxf+zYsdqoUSNt166dPvLII0E9L97Kei3L+9kL9hq8E+HSniPvnytVDfpn64cfftCTTz5ZAT3//PP1ggsu0AsuuEBHjBihV1xxhT733HN64403+n38E088oVu3btXly5frDTfcoKqORDScynvfBvpzu3jxYh08eLA2atRIu3fvruedd57fsfgzjnDUk9hQLJYIm4rqcd9Xfh1nLcfRNXHiRHcrypYtW3TVqlWqGt5EuKo544wz9MMPPwxJWZmZmcXWe/bsWeIPuj/HVDahfI6qMs9Ez56TkrwT4Uj/bJVlxIgROmvWLL+Pf/LJJ0tMvHPLLbf4fX5VF7UJNYwxJpQmTJjAvHnzmDZtGq+++irJycl888037Nixg2+++abY45ro66+/pk6dOpx//vkhKe+yyy4jIyMDgCVLlrBjxw4GDhwY8DGVSaifo+rAnpPyReNnqyxDhgzhlFNO8fv4MWPGcPfdd/PXX3+Rk5PDSy+9VGMm4vFXlbhZzlQ90xfa+LwmdFq2bMmrr75abNuwYcNYu3ate93zcU0zaNCgkE4MMHToUIYPH069evU4fPgwH374YYm72/05pjIJ9XNUFR04cIAGDRqQn59PnTp16N+/f41/Tjx5Pz8QnZ+tsnjeVOaPli1b8uSTT/L000+TkZHBWWedxYUXXhho2NWaqFbtu8YBRKQ1UGz6kEmTJtmsclE06rkFfH7zSdEOo1xH/OdrVj9Q/jSVVeV6jDGmNPfccw9t27ZlwYIFTJs2jYSEhGiHVKnY81O9pKSkMHnyZM9NyapaYjK2apcIp6WlRXbYDeNTVUkcLRE2xhhjqqft27d7TibjMxG2PsLGGGOMMaZGskTYGGOMMcbUSHaznKmR+j04h4zsPAA6TJxFs/q1WHzf0ChHZYwxxphIshZhUyO5kuDS1o0xxhhT/VW7RDg5ORkRsREjKqG8vDyqw82ZxhhjjKncUlJSPG+UK1W1S4TT0tJQVUuEKwFV5YUXXmDKlClceeWV1K9fn/bt2/PDDz9EOzRjjDHGVGMpKSmkpaWVe1y1S4RNZO3bt48RI0bwwAMPlNj30EMP8eKLL7Jo0SLatWvH7t27ef7557nyyispKiqKQrR/a1a/VpnrxhhjjKn+7GY5U2Gqyrhx42jcuDHPPPMMY8eOpUOHDgDs3buXV157hYULF9KqVSv3OWeffTYTJ05k0aJFHHfccVGKHPeNcf6OI2yMMcaY6sdahI3fMjMzmTx5Mrt27QLg8ccfZ8eOHbz55puMHTuWl19+GYClS5eyefNmPvnkk2JJsMvpp5/OggULIhq7McYYY4y3gBNhERkkIjeIyCQRuVVEzhSR+HAEZ6Jvzpw5LFq0CID777+fl156iYsvvpjZs2fz1FNP8cEHH5CQkMB1113HG2+8wYYNGzj33HNp374dffv29Vlm3759WbJkSSQvwxhjjDGmBL8TYRE5XkRWA/8DTgc6AycBE4ENIjImPCEGxkaNCJ033niD8ePHc/bZZ/Pggw/y1ltv8fvvvxMTE8PIkSN57733aNu2LQDdu3dn+PDhdOvWjZtuuonGjZuUWm7Pnj1Zu3ZtpC7DGGOMMTWMv6NGiD/DWYlIV+BfwD2qusvH/ro4EuLfVPXzwMP1nzOWVcBQVZ3n3NYaSAPHqBGtW7cOZwg1xlFHHcVLL71EUlISI0eO5JJLLuH+++/n8OHDHDhwgGbNmhU7XlU58dHv2b4v170tObEOP008tdhxmZmZdOrUiaysLEQkItdSGn/7CI96bgGf33xSBCIyxhhjTChs377dMxlOVtXt3sf4e7NcFnC1qvq81V9Vc4D/ikjLigQaoPsBm/0gBLZv307Tpk1JSEgose/gwYNs2LCB448/nri4ONatW+fel5CQ4PMcEWH7vlxSHx3p3tZh4qwSxzVp4mgt3rt3r/txZdXjP1+Rm+9423eYOIvacTGseXBElKMyxhhjTCiUmwiLSG1gn2cSLCKtVHWH97GqujPE8XnH0h/IBnaHs56aYOrUqVxzzTWccsopzJ49u8T+575cSqtxT3PeS78GXdeo50reGNds9OP849XfqFu3XqnnXTqgHaMHtgu6fl/8nWI5N7+oRGL/8ccfc/7554clLmOMMcZETpmJsIjcBowCCkRkA3C7qh4C3gVOLevcMJkMjAdKZizGb5s2beLOO+9k+fLljBgxgj/++INevXoB8NFHH3HyySfT5nAqA/bO5b2brw2obF8twL66FJxzzuNcMfQKzj9/WMUuIkj+TLE84b1lPs+97Sf49uAiXhw7IByhGWOMMSZCymsRvkBVTwIQkaHA5yJyffjDKklERgArVXVbtPuVVnV33XUXt912Gz179mTkyJHMmTOHXr16sWrVKi688EJuuOEGGjVqRM+ePQMuOzmxTrFkODmxjs/j2rdvz5YtWyp8DcGYvtB3vaOeW1CsFXpI7VQ+peSAKKfu+ZykTR0BS4SNMcaYqqy8USNiRSQOQFXnAOOAF4FuYY6rGBGJAe4BHinvWNeoEa7FRo/4W35+Pm+99RYLFy7k9ttvB+CUU05h7ty5gOMOy7Fjx/LVV1+xatWqCiXChwsKy1x3adeuXdQS4dK6W3x+80nufYsXL+aaa66hdlwMHSbOci+142O47bbbePXVVyks9H1txhhjjImOlJQUdw7oz6gR5SXCdwLuoQFUdRswErgvqCgDNxr4RlX3lHdgWloaqupeLBF2UFXOPvtsnn/+ed5//33q1HG01A4ZMoQff/yR33//nR9++IEXX3yR3NxcZs6cyZFHHhlwPf50OQBHIrx58+bALyREyppiOSMjgwsvvJApU6aw5sERpD46kqOTG5H66EjWPDCCXr160apVK7755ptIh22MMcaYMqSkpLhzwLS0tHKPL7NrhKr+5GNbHjC1whFWzCDgKBFxdShNAp4Wkc2qek6EY6mSZs+ezdatW1m2bBlxcX+/7M2bN6dTp06MGDGCu+66i/r16zNgwAC+S83lrjkZyLfBzwDn62a5gwcbsaXlcJ/7XMJ5s1xpUywXFhYyZswY/vGPf3DBBReUev6YMWP48MMPOfPMM8MSnzHGGGPCz9/h0wAQkbtV9bFwBVMaVS12x5aIpAITXOMIm/I9+eST3HPPPcWSYJfrr7+eV155heuvv969nvzZZ7xwy6CA6/H3ZrkdO3bQu/flfL4zrAONBOy+++4jLy+PRx4puxfO2WefzUMPPURhYSGxsbERis4YY4wxoRToFMvlzzwQRiIyQETm8XeL8H+jGU9VkZGRwcKFC7nwwgt97h8/fjyLFy+mXj3HUGbDhg3jhRdeqFBdZXU58NSyZUv27dvHoUOHKlRPqDSp93d87777Lv/3f//HBx984PMDg6eOHTuSlJTEr78GP7ycMcYYY6Ij0EQ4qsM1qOoiVR2iqrVVtbeq3u99jE2xXNJnn33GGWec4e4XHE6L7xtarE+tr7F5AWJiYmjTpg1bt24Ne0xlcSXCP/74I7fddhufffZZiRnzSnP22Wfz+edhnUjRGGOMMRXg7xTLgSbC5c/HHGWum+UsEf7bxx9/zHnnnRftMEqI5sgRntatW8c//vEP3nnnHY466ii/zxs1apQlwsYYY0wllJKS4tfNclWqRdgELjc3l/nz51fKm7qiPXIEOG6OGzVqFA888ABnnHFGQOf279+fzMxMNm7cGKbojDHGGBNOgSbCN/jaKCJdQxCLCYMFCxZw9NFHk5iYGO1QSojmpBpH/OdrOkycxar0bA6f9z/+t6VtwGXExMQwcuRIvvjiizBEaIwxxphwCygRVtXVIhIjIqeJyGUicrmIXA5U7M4qE3azZ88OuKUzUqLZNeJQfiGpj450L4fyy54cY89B3+MhW/cIY4wxpuoKtEUY4DMgBTgdOMW5lN8b2UTFnDlzGDrU9w1r0dalSxdWrVoV8XrLmmK5tH2lJcJDhw5l4cKFZGZmhiw+Y4wxxkRGRRLhRFUdpKrjVPVKVb0Sxwx0lUJNGjXi22+/5ffffy91/969e9mwYQMDBgyIYFT+O+6441izZg3p6ekRrfeS/m18bvecYtlf9erV48ILL+Tll18ORWjGGGOMCQF/R40Q1cAGghCRh4A3VHWDx7abVfW5gKMMERFpDaSBY9SI1q1bRyuUiCksLCQuLo4OHTqwceNGRP6+j3HVqlUcf/zxHHXUUdStW5c5c+ZEPL5Rzy3wOZGGtxtvvJG4uDieeeYZwHFd8+bNY/DgweWO5VtRTz31FFMyunO48O/3fp342GIzzLmc+Oj3pGX9PdZxcmIdfpp4arFjVqxYwbBhw9i0aRMJCQlhidkYY4wxgdm+fbtnMpysqtu9j6lIi/BvwFIRSRORjSKyCXggiDhNBfzyyy/06tWLhIQEfvzxx2L7nnnmGcaOHcvPP//MJZdcEqUI/ZOSksKnn37KlClTSE1N5fTTT+eMM87ggQfC85ZauXIlDz/8MF9e1b3YeMe+kmCAtKxDxfoSeybFLkcffTQ9e/bko48+CkvMxhhjjAmPirQIrwZuBDbiGFdYgBRVHRfy6PyPqca1CD/++OOkpaVx7LHH8p///IdZs2Zx9NFHk5eXR+vWrfn9999p2bJl1Foo/W0RBlizZg1XXHEFq1ev5u677+aSSy7huOOOY8uWLdStW9fvOqcv3MKMRaXffJeXd5i1a9fRqlUr96QZf+3KpkuL+qWesyJtH6mPjnSvd5g4i6OTG5U4LitrL+npO+nRo4ff8QJcOqBdwN0xjDHGGFM+f1qEK/Ld8wpV/d5zg4hYi3CELV26lGHDhnHFFVeQlZXFf/7zHz799FO++eYbevbsSbt2VSe56tGjBwsXLiy27fjjj+e9997jqquu8ruc0QNLTypXrFjBqFGjmDBhAhMmnOveXl7C3mHirBLbfB1fUFBA586deeCamfTu3dvvmI0xxhgTPRXpGrFFRN4UkWs8hk97KdSBmbItW7bMnXBdcsklzJ8/n8LCQqZOncqYMWOiG1wI3HDDDbzwwgsE+o2Ft9zcXJ566ilOO+00Hn74YSZMmBDQ+cmJdegwcZZ7SU70PU11XFwc48aNY9q0aUHFa4wxxpjIqUgifKnz/xOohMOn1YRRI3JyckhNTeWII44AoGXLljRv3pyvvvqKefPmVYtEeNiwYWRlZfHzzz9X6PyioiKmTp1K165dmTdvHt9//z2jR48OuJyfJp5K6qMjqRMfS+qjI0vcKOdpzJgxzJgxg4KCggrFbIwxxpjQ8HfUiIp0jXhAVad4bhCRkaUdHGk1oY/wihUr6NGjB7Vq1XJvO/HEExkzZgxXX3019euX3ue1qoiJiWHixImkpKQEPOpFbm4uY8eOJTU1lQ8++IDjjjsuTFEW161bN9q3b8+3337L8OG+b74zxhhjTPilpKQwfvz4cpPhgFuEvZNg57aSHSlN2CxdurREP9STTjqJ/fv3c+ONN0YnqDAYN24cGzdu5IcffvD7nOzsbM4880zAMb10eUnwpQPK7kvd78E5dJg4i0P5hXSYOIt+D5adlF922WW88847fsdrjDHGmOipSNcIE2XLli3j2GOPLbZt7NixLFy4kM6dO0cpqtCLj4/nv//9L5MnTy6xb+3atSxatKjYtgMHDjBixAg6duzIe++959eIGeWN2JCRnVfmurdLLrmEL774guzs7HLrNsYYY0x0BTx8WmVU3YdPy8/PJz4+3r0+cOBAnnjiCQYNGhTFqMoWyPBpZcnLy6NTp0588cUX9O7dG1XlmWee4aGHHiI+Pp67776bCRMmkJ6ezhnXT4bOJ9GufTsEKb9wP6xI21dim6/h0zz99dd6GjduQtOmTcst34ZPM8YYY8IjpMOniUgPVV0TquCMf3799VeOP/54cnJyqFOnDgUFBfz5558cc8wx0Q4tImrVqsXNN9/MAw88wFtvvcX48eNZs2YNCxcuJDY2llNPPZU5c+awaNEibrjhBiZNupSYmNB90eHv8Gme3n9/O6+88hyff/ttyOIwxhhjTOgFkjG8IiLHiedcvpVQdRs1Yvbs2QD8+eefAKxb55gQomHDhtEMq1zl9b0NxE033URaWhpNmzYlISGBBQsW0KlTJ9q3b8/ChQsZM2YMP/30E5MnTw5pEgzQrH6tMtd9GTVqFL///jvbtm0LaSzGGGOM8Y+/o0b43TVCRN4DZgEDnJsWA1+r6s6KBhkIEbkZOBfHTHbNgVdU9TnnvmrbNeLKK69k6tSpvPnmm4wbN47p06fz8ccf8+GHH0Y7tIgqKioiKyuLJk2aRKX+QLt6jB8/ns6dO3PPPfeEMSpjjDHGlMafrhGBNJ+NUdW3VfVm4BZgNTBeRL4SkYeDD7dc/wQuVdVTgfOBp0RkQDnnVHmpqakMHjyYdevWAY4RI7xvlKsJYmJiopYEV8TYsWN56623gp4QxBhjjDHh43cirKqFACLSEEcieq1z6Qu0CUt0xY1V1V3OWNYDe4EOEag3qjZt2sQZZ5zhToQ9Z5QzldeJJ57IoUOHWLp0abRDMcYYY0wp/E6EReRfIvIDkA5MADbh6KrQUlUvD0t0HlT1D49YLgCygdnhrjeaCgoK2LFjB6eddhrr1q1DVX0OnWYqn5iYGHersDHGGGMqp0C6RhwNtAL+A1ykqg+q6mKN4He/InK0iKwC/gdcoqpZ3se4bpZzLVX5prlt27bRsmVLjjzySP766y9SU1OJjY2lVatW0Q7N+OGyyy5jxowZ5OfnRzsUY4wxpkZISUlx54D+3CwXSNeIS4HuwE84+ga/LCIPicjJIhJb8ZD9p6orVLUnjpbomb76CKelpaGq7qUqJ8JbtmyhXbt21K9fn8aNG/Ppp5/Sv39/KvnAHcapa9eudO7c2T3yhzHGGGPCKyUlxZ0DpqWllXt8QGNNqWqRqv6qqpNxtAxnAp86/48YVV2KYwSLWyJZbyR899137hZEVyIM0K1bN95++2369+8fzfBMgC6//HLrHmGMMcZUUoH0EX5JRE4UkQdFZDGwBTgbeAwYHK4AnXU3dfYL9nQQqBfOeiNt//79nH766cydOxcomQgvXbqUoUOHRjNEE6CLLrqIb775hr1790Y7FGOMMcZ4CaRF+FrgPSAJeARooapDVPVRVV0eluj+1gD4j4jUBRCR5ji6R3wX5nojatOmTQBs2LABgK1bt9K2bVsArr32WsaPH8/AgQOjFp8JXJMmTTjnnHN48cUXox2KMcYYY7wEkgi/rKptVfWfqvqRqu4PW1QlpQOfA9+JyHzgW+ANoFplF5s3bwYgIyMDKN4i3KdPH15++eWQz5xm/BPMTHn33nsvzzzzDNnZ2SGMyBhjjDHBCuRmuevDGUg5deeq6n9U9XhVHayqx6hqiqoWeR9bladY3rJlC1A8EXa1CJvoGj2w4olwjx49OOWUU3jhhRdCGFHlN33hlmiHYIwxpobyd4rlate86Bo1oiomwps3b+bII48kMzMTVWXTpk107Ngx2mGZELj//vt5/PHHSU9Pj3YoETNjkSXCxhhjoiMlJSX0o0aY8Nq8eTN9+vQhIyODXbt2Ubt2bRo1ahTtsEwIdO/enSuvvJJbbrnFpl02xhhjKglLhCuRzZs307dvXzIzM9mwYQOdO3eOdkgmhCZPnszWrVuZMGGCJcPGGGNMJRBwIiwiV4UjEPN3IpyRkWGJcDVUt25dvvzySxYuXMh1111HUVGJLu7VgvUNNsYYU1WUmQiLyBwRme1c5ojIHOBREbGpskIsNzeXvXv3ctRRR5GRkcHGjRvp1KlTtMMyIda4cWPmzJnDmjVruOaaaypVy3CoEtj/zvyTDhNnsSJtHx0mzqLfg3NCUq4xxhgTauW1CM8C0oDHVXWoqg4FlqjqGeEPrWKq6qgRW7ZsoU2bNjRq1IjDhw+zatUqaxGupho0aMBXX33F4sWLeeONN6Idjluobm4rKCqe3Gdk54WkXGOMMcZf/o4aEVfWTlV9WkTigLEiMgP4jErerzgtLY3WrVtHO4yAbd68mfbt2yMiNGvWjIULF3LjjTdGOywTJnXr1uXtt9/mjDPO4KKLLqJBgwbRDqlU0xduKXX4uOkLt/iVQI96bkGx9UsHtAtqSDpjjDGmLCkpKYwfP77cZLjMRBhAVQuAN0VkKnA+MD8kEdZAmzdvJj8/ny5duvjc1759ewCaNm3Kn3/+Sbdu3SIdoomgXr16cfrpp/Pss8/y73//O9rhlGrGotIT4dEDSya0HSbOKnHc5zefFJbYjDHGmGCUmwi7qKMz40cAItJTVVeFLapq6sILLyQjI8M9lbInz1nkDh8+DEBSUlJE4zORd99993HyySdz6623Ur9+/ajE0O/BOe7uCx0mzqJZ/Vosvm9oqdvLExcjxbpHNKtfKzyBG2OMMUHyOxH28jRQafsJV1arVq0iJyfH574tW7YwaNAgAIYNG0bTpk0jGZqJkh49enDqqafy0ksvcdddd0UlBu8+vK710raXpsd/viI3/++RMBLiYlj74IgQRWmMMcaEnpR117qI7Ad+BQRwHShAb1VtHv7w/CMirXHc1Fep+wg3b96cjIwM8vPziYsr/hnkzDPP5MYbb2TkyJEUFRUhIohIlCI1kbRixQqGDh3Kxo0bqVu3LgB79uxh+vTpLFy4kGXLlnHo0CFOOeUUBg0aRJs2bWjZsiXJyckkJib63U/XZc/BPPYc/DupPZRfWOKYOvGxPrcfnVzaBC/KirT9pD46khMf/Z60rEPuPfVqxdKpue/WbusrbIwxJly2b9/u2Uc4WVW3ex9TXovwOuAiVc3y3OgcRq1Scl3wpEmTKtXIEUVFRezbt4+kpCRSU1NL9BPeuXMnLVq0ACAmplLfj2hC7Oijj+aEE07g+eef5+yzz+bZZ59lxowZjBw5klNPPZXbbruN+Ph45s6dy9dff82OHTtIT09n27ZtNGzYkNGjRzPt3ntp0qRJher31ad39QPDy+zre/jwYZYtW8Yvv/zCggUL+OGHH6h71ZsApGUdIvXRkcXKf2pEks++8cYYY0w4pKSkMHny5HKPKy8RHgoc8N7oHEatUqqsLcK7d++mUaNG9OrVi7Vr15ZICnbt2uVOhE3N8/DDD3P66afzv//9j/Hjx7N69eoSfcSPPvpobrnlFve6qrJu3TqeffZZjjzySF577TVGjhzpXXS5mtWvVazbg6tPb+34mGJdHWrHxXDTTTexaNEiVq5cSbdu3Wg96B/sPWosxx5/A2t3HypRtsvpj82h55FpxMbEAtYSbIwxJrxCMmqEqu71XBeRu1X1sRDEV+OkpaWRnJxM165dOeuss9i7dy+JiYmAI6GxRLhm69GjB1u3bkVV/f5GQETo3r07L7zwAhdffDGXXXYZF154IY888ggJCQl+1+26AW7UcwuKje6Qm19UomW3U6dOXHLJJfTp08fdjcN9Df/5ymcrMsCgQz+TtG0jjz/+uN9xGWOMMeFWZh/hEgeLfK+qp4YxngqpCn2Ev/jiC1588UWefvpphg8fzo033sgdd9wBwL59+2jbti379++PcpSmKnL1ES4oKCA1NZW8/Dw6depE7YTafp2/cXc2B/P+7g9cJ06Q7F3k1G5eIhEuvY/w39btPMDhgr9bkuvViqVtYgIrVzlakevW+TuBtpZhY4wx4RKKPsLe7O6tCkpPTycpKYlu3boxY8YMLr30Um6++WZq1aplrcEmKJ5j+aoqL774IikplzJ+/HhOPfVUjjzySFq0aOGzpVlV6XjvlyUS3snHCpNWl6zLn/GAvVuWXV57bQNTpvyLz+bPp169egFcoTHGGBMegSbC/jcf10DLly/n0KFDHHfccSX27dq1i5YtWwIwYMAAevTowXPPPccdd9xhibAJGRHhxhtvZOjQobz66qukpKSwevVqDhw4QMuWLWnevDmx3QZzoGUvCvILOHz4MDHNOpQo58P9nalXK7tYV4fkxDpBxXbVVVfx888/M3z4cD7//HN31yBjjDEmWqpdi3A0R4246KKLyMjIIDMzs8S+Xbt2uWeOExGuuOIK3n//fUuETVh069atWH/c3Nxc0tPT2b17N5mZmeTk5NC8eXM6duzISc8vL3F+OGaCi4mJ4bXXXuP2229n8ODBzJ07t8IjXRhjjDFlCdWoEd5uqFg4wRGReOBG4FwcyXgt4D5V/c772Gj1EXbdwV+nju9Ws507d9K/f3/3eufOndm4cSNgI0aY8KtduzYdOnSgQ4cOJfYlJ64LactvWWJiYnjqqae48847Ofvss5kzZ06pPzPGGGNMRYVk1AhvqrpaRJKADs5NqaqaXrEQA5IM3IpjIo99IjIUmCki3VU1LQL1lysjI4P4+PhSJ8Hw7BoBjumTd+7cCRQfQ9iYSPtpYmTvfxURHn/8ccaOHcull17Khx9+WGKCGWOMMSYS/J65QUQ6OCfS2Ax86lw2i8i3ItIpPOG5HQD+q6r7AFR1DpALnBDmev22bt06jj32WPLy8sjPzy+x3zvZbdq0KXv27HEPneaZJBtTlV06oPxRIGJiYnjzzTfJz89n2LBhpKdH4vO0McYYU1wgU5g9DTwF1FfVJFVNAuo7tz0VhtjcVDVTVd92rYuj2bUWsDuc9QZi3bp1dO/enUaNGrFv374S+72T3Tp16hATE0NOTg67du2iefNKM2O1MUHxdzi0WrVq8dlnn3HSSSfRp08f5s2bF97AjDHGGC+BJMI7VPVLVXU3d6pqvqrOAkqMyxZmg3G0TP/gvSM5ORkRcS+RumFuy5YttG/fnsTERLKysortKywsZO/evTRt2rTY9iZNmrBnzx52795tibCpkWJjY5k8eTLTpk3jkksu4fbbb2fz5s3RDssYY0wVlZKS4s4By+sfDIElwkkicpaI1HJtEJFaIjIKSCrjvJASkdrAw8A4VS3y3p+WloaqupdIJcLp6em0atWKhg0blmgRzsjIoHHjxiX6QTZt2pTMzExLhE2NN3ToUJYsWUJMTAx9+vTh5ZdfjnZIxhhjqqCUlBR3DpiWVv5tZIEkwrcA1wEHRWSPiOwBsoFrgQkVCTZQzi4RrwBPqeqSSNTpL9eEGXXr1iUnJ6fYvtJGhbAWYWP+lpyczBNPPMFvv/3G5MmT+e67EoPCGGOMMSHl963aqroVOEtEmgDtnZs3q+qesETm25PAIlX9QEQSgJaquiWC9ZfKlQjXq1ePgwcPFtu3c+dOnzfDNWnShIyMDPbs2UOzZs0iFaoxlVqnTp145513GDNmDL/99htt2rSJdkjGGGOqqUBahAFQ1T2qutS5RCwJFpF7cCTuU0WkPtAZuCpS9Zenoi3CGzZsoH79+sTHx0cqVGMqvVNPPZUJEyZw4YUXcvjw4WiHY4wxppoKOBH2RUTODkU5ZZTfDXgUuBnHUGoHgJXhrDMQqlpmi3BpiXDTpk1Zu3atdYswxoe7776btm3bcs0116Bqs7sbY4wJvUDGEe5W2gJcGMYYUdV1qio+lhTvY12jRkRyeuWDBw8SExND3bp1fbYIl9U1whJhY3wTEaZNm8aqVat49NFHox2OMcaYKiQlJcWvUSMCmc7pd6AOjimOvSlweQBlhU00pljOysqiUaNGAKW2CHfu3LnEea5E+OSTT45InMZUNXXr1mXmzJkcd9xxdOvWjQsuuCDaIRljjKkC/J1iOZCuEZOAk1Q1xnsBXgwm2Kpu37597kTYV4vwjh07SEoqOcJckyZN2Lt3r1+fWIypqZKTk5k5cybXXXcdv/76a7TDMcYYU40EOrPckaXsezf4UKouz0TYV4twamoqHTp0KHFekyZNAGjbtm3YYzSmKuvTpw9vvvkm5557LitXVprbA4wxxlRxgQyfVgi8Vsq+Gt1Ms3///mItwjt37uTER78nLeuQ44CzH2f85zv5+aijip3nSoTbt2+PMaZsZ511Fk888QRnnHEG33//Pd27d492SMYYY6o4vxJhEekPZKnq+jKOaQgMck65XKP46hqRlnWI1EdHuo/pMLHk0+LqNzxgwIDIBGpMFXfZZZeRn5/P6aefzty5c+nSpUu0QzLGGFOFiT/DEolIPPAh8DMwB9gCHARq45he+STgYmCsqu4IW7Slx9caKDaP3qRJk8I6csT0hVuYscgxl8fu3bs5mHOQDu07kLknk31Z+ziY0LREInx0cqMS5agqjgnz4NIB7Rg9sF3YYjamunjllVd46KGHmDdvHh07dox2OMYYYyqZlJQUJk+e7LkpWVW3ex/nVyIM7mT4LmAs0B3HSBEAO3AkyQ+oamYwQVeUZyIcjVEjHn/8cXbu3MkTTzzBJ598wrRp01jW45oSibDnujEmOC+88AJPPPEE8+bNs+5FxhhjSti+fbvngAQ+E+FA+gjnAw8DDzuT4mbAPlXNKfvM6s/XzXLJiXWKdYdITqwTrfCMqZZuvPFGioqKOOGEE/jggw844YQToh2SMcaYKqZCM8upar6q7rAk2MFXIvzTxFNJfXQkw7M+5+52m/lp4qlRjtKY6ufmm2/m1Vdf5dxzz+WVV16xGeiMMcYEJCRTLNd0nolw/fr1yc7Odu/btGmTz6HTjDGhceaZZ/Ljjz/y/PPPM3r0aPbv3x/tkIwxxlQRlgiHQFmJcGljCBtjQqd79+4sXLiQRo0a0bdvX+bPnx/tkIwxxlQB1S4RTk5ORkTCOmKEt9IS4aKiIrZs2WI38hgTAXXq1GHKlCk89thjXHHFFVx44YWsX1/qiI/GGGOqsZSUFL9m7q12iXBaWhqqWikS4a1bt9K0aVPq1asXsViMqenOO+88Vq9eTZ8+fTjhhBO48sor2bhxY7TDMsYYE0EpKSmkpaWVe1xQibCIJInIYBG5RkSeFJFPgymvqtq3bx8NGzYEHK1Shw8fprCwkHXr1tGtW7coR2dMzVOnTh3+9a9/sX79etq1a0f//v355z//yfz588nPz492eMYYYyoJv4dP8yQibwHtcUymcQh4BHgLWBu60KoOzxbhmJgY9+xylggbE12JiYlMnjyZW2+9lRdffJHbb7+dv/76i/79+9O/f38aNmxIUVEReXl57qWoqIhmzZrRqlUrkpKSKCoqYufOnRw4cIATTjiB/v37uyfBMcYYU7VVKBEG1qjq5SISA4wHRgOPqmpu6EKrGlS1WCIMf3ePsETYmMqhSZMm3Hfffdx3331kZGSwaNEiFi9ezL59+4iJiaFWrVo0aNCAWrVqISJkZGTwww8/sGPHDmJjY0lKSqJ27dpMmTKFgoICLr/8ci6//HK7EdYYY6q4CiXCqvqw8/8iYIqIvA3cIyLjVfXKUAYYKFfH6HBPseySk5NDfHw8tWrVcm9zJcJr1qxh6NChYY/BGOO/Zs2aceaZZ3LmmWcGfK6qsnjxYqZNm0a/fv046qijGDduHJdeeikJCQlhiNYYY0xF+Jhi2adAplhuB+SpanoZx7T2NX1duEVziuXt27fTp08f0tP/flqOPfZYXn/9dYYPH86SJUto27ZtxOIxxkTG4cOHmTVrFq+88gp//vknd9xxB+PHj7ebY40xppLwZ4rlQG6W+wg4V0TaOZfm3geEOwkWkf4i8peIjAtnPYHw7hYBjtnl1q9fT1FREW3atIlSZMaYcEpISOD888/n66+/5rPPPuOnn36iTZs2XHDBBbz66qukpqZGO0RjjDHlCCQRXqyqU4ArgS+AYeEJyTcROQ+4DdgXyXrL4ysRrl+/Pj/++CPHHnus3VRjTA3Qp08fPvzwQ9asWcM555zD/PnzGThwIN26deOmm27i008/JTMzM9phGmOM8RJIIqwAqjoZWKmq77h2iEgkOsf9pqqjgQMRqMtv+/fv95kIL1iwgGOPPTZKURljoqFly5ZcfvnlvPPOO+zYsYP333+f9u3b89JLL9GxY0d69erFjTfeyPTp09myZUu0wzXGmBqvoqNGeHcsfhy4JchYyq5QdVs4y68ozzGEXerXr8/y5cu55557ohSVMSbaYmJi6N27N7179+auu+6ioKCApUuX8sMPP/Dhhx9y2223UatWLU488UT69etHly5d6Ny5M61atSIxMZG4uIr+ejbGGOOvQFqEzxeRn0XkZ2Co67Fz/cIwxRcw1xTLriXcI0f46hrRokULAGsRNsa4xcXF0b9/f+644w4+/vhj0tPTmTt3LsOHD2fr1q289tprXHTRRXTv3p3atWvToEED2rVrx0knncTTTz/N1q1bo30JxhhT6aWkpLhzQH+mWA6kyWEVMK2UfZcFUE5YhXvUCFUlLS3NfROcr0S4ZcuWAHTt2jVscRhjqjYRoUuXLnTp0oVx48YV26eqHDhwgKysLFavXs2MGTN44IEHSE5OpmvXrrRp04ZWrVq5J/1wPW7atCkxMUFNGGqMMVVaSkqKuxHUa9QInwJJhO9T1Z997RCRGjOj3Pfff8/pp59Ofn4+cXFxPvsIn3/++bye2ZXO//7ava1Z/Vosvs/GFDbGlE9EaNiwIQ0bNqRdu3YMGzaM/Px8li9fzsaNG0lLS2PHjh2sXLmSHTt2uJfs7GxatmzpTpCbN29O8+bNady4sXtJTEykcePGNGrUiAYNGtCgQQPq1q1rN/ZGwfSFWxg9sF20wzCmRvM7ES4jCW4BNPS1rzr6888/Adi2bRsdOnRg3759tGtX/BdZx44dySlaVWxbRnZexGI0xlQ/8fHx9OvXj379+pV6TG5uLjt37mTHjh2kp6eze/duMjIyyMzM5K+//mLv3r1kZWWxd+9e9u3bx4EDBzhw4ACHDx+mTp061KlTh9q1awf02LUeFxdHbGwscXFx7qVWrVokJCRQp04d6tevT4MGDahfvz61atUiJibGvcTGxhZbd50bFxdXrRP0GYssETYm2gK6G0NEegDHAL09lvrAMmB2SCOrpNaudTR+v/PLJhZ+vo3U2L7U39eAuc8tKPfcUeUcc+mAdvZL0RhTYbVr16Z9+/a0b98+oPMKCgo4dOgQhw4dIjc3N6DH2dnZ7N69m4KCAgoLCyksLKSgoICCggLy8vI4fPgwOTk5ZGdnk52dzYEDB8jPz6eoqIjCwkKKioqKLa7z8/PzKSgoID4+3j17Z61atdyPvf93Jdy1a9embt26Ppd69eqV2OaZ0LuW+vXrk5CQUK2TcGOMg9+JsIgsALoBmUB74G4c4/quVX+npwuCiPQFnsSRfE8UkbNV9fxw1+tt06ZNAPRpdIj7Lj2F88//H6NHj+bCC08qdlyHibNKnPv5zSeV2GaMMdEWFxfn7iZRmagq+fn55OXllfv/4cOH3cl5Tk4Ohw4d4uDBg+Tk5JCTk8OOHTvcj3Nyctz7XOd4nnvw4EEKCwupX79+sZZs11KvXj1q167tTs5dS0JCArVq1aJ27drUq1fPfaz3Muqtv8jNLwIcfytqx8ew5oERUX62jamZAmkRbgW0VtUCEbkVuBVYr6prwhNacaq6BBhS3nGuTtGTJk0Ky4gRaWlp9OzZkz179gC+xxEGR59gz+4QzerXCnksxhhTnYmIO8mMtLy8PA4ePOhuyfZcDhw4QF5enns5fPiw+/HBgwfJzMwkOzubgwcP+lwKLnyG1EdHuuvqMHEW7dq1o169ejRo0ICmTZvSrFkz9+Jab9WqFW3atKF169bEx8dH/DkxpipJSUlh8uTJ5R4n/jbmikgdVT3ksd4NeBNYDtzhuS/SRKQ1kAbhHzWiadOmDB06lBNOOIFbbrmFPn368Morr5Tab89uhjDGmMib8N4yvlmZHu0wfDqUX1giEU6IE1BQxz+oqmPB8T/O9SLn/4IgMY4homJEEIlxDBkVI8RIDI5eHb67dhzKLyyxrU58bIltw45M4ulLeofoqo2JPK9RI5JVdbv3MYHcLHfIa32diAzC0UViEXB0ELFWCa7+cN26dXO3CGdkZNC8efNSz7Ek2BhjIq8yJ3C+us6tffBMv88vKChg586dbNu2rdiyefNmNm3axMaNGzl8+DCdOnXyuZz9zibyCorc5dWKi2H1A8NDcm3GVDVBTV2kqkXAoyLyeYjiqdR27NjhHqtzw4YNqCq7d++mWbNm0Q7NGGNMFVE7PqZYMlw7PrCxn+Pi4khOTiY5OZmBAwf6PCYrK8udFG/cuJGVK1fy9s8b2dWgG9q4bYkW6RNSZjKiW0OuOfUIWrZsaTcKmhojJHN4qurKUJRT2aWlpZGcnEyTJk347bffyMnJQUSoV69etEMzxhhTRbhujBv13IKw3USdmJjIscceW2yG0+kLtzBj0RZWpO0rcfy+ffv44IOveW7CRxw8eJAOHTrQtm1b2rRpU2xxjUpSt27dsMRtTKTZZPYBcCXCTZs2JTMz01qDjTHGVNilAyLbdW70QMcQnb66Zqx86nLgcuAlDhw4wKZNm4p1u/j555/ZsmWLe2nQoIE7Kfa1NG7c2FqVTZVQ7RLhcI4a4Z0IZ2RkWCJsjDGmQqJ1D0lyYp1iyXByYp1i+xs0aECvXr3o1auXz/OLiorYtWsXmzdvdi/r16/n22+/da8XFRWVmSgnJSXZdOAmrPwdNaLaJcLhHDUiNTWVLl260KRJE/bs2cPOnTtp2bJlWOoyxhhjwuGniacGdX5MTAxJSUkkJSWV2Ud58+bNpKamupPj3377zf14//79tGvXjq5duxZbunTpQtu2bW14OBO0lJQUxo8f7zlqhE/VLhEOlW+//ZZt27Yxbtw497aNGzdyxhlnuFuEN2/eXGJ6ZWOMMaamS0xMJDExkWOOOcbn/pycHFJTU1m/fj3r16/njz/+4KOPPmLDhg2kp6eTnJzsHuWiY8eOxUa9aNKkSZXsdtHvwTkl5hdYfN/QKEZkwBJhn1SViy++mD179pRIhDt16kSjRo04cOAAGzduDHgqU2OMMaamq1u3Lj179qRnz54l9uXl5RUbCm7jxo188MEHbNq0yT1iU/v27WndunWxpXnz5jRp0sS9NG3alHr16lWapNkzCfa1bqLDEmEfVq1aRYMGDSgqKnJ3fygqKiI1NZWOHTsSGxtLo0aNWLJkCddcc020wzXGGGOqjVq1arm7Sviyd+9eNm/ezPbt293L8uXLyczMJDMzkz179rBnzx5yWh1LQs/TiIuLIy4ulti4OOJi44iNjSUmJsb9f0xsDLExrm0xxMTEEhsbg8TEECMxxMQ4Jizx/L+0yUrKvK64mBLjN496bkFFn6awuHRAuxo3/0G1S4RDcbPcd999xxlnnMHGjRtZsmQJZ555Jlu3bqVJkybuIWM6d+7MvHnzeOGFF0IVujHGGGPK0bhxYxo3bkzv3r3LPTY3N9edGO/Zs4fMzEwOHDhAdvY+n1Nou7YdOHCAQ4cOcfjwYXJzc8nNzXU/zsvLo1atWtSpU4f69etTr169Yv+Xtm1FQZcS4zc/emrjYsclJCRUmhbsqs5ulgvCt99+y2WXXcby5ctZuHAhZ555JkuXLi32Q9epUyd+++03unfvHmTExhhjjAk117jJf4sBmjsXoLZz8Rr8qY5zKZ1jquuioiKKCosoKiqisKiQosIiDhQVkeV8XFRUSGFREUUHiijcV+io3suo536k0FlGUVEhqkqMV+v0363WMRXaFxPjXKT8luzq1CJcY2+WmzVrFnXr1mXMmDEVOj8/P5/58+fz+uuvU7t2bV588UWAEonw5MmTufDCC4mNLTk/uzHGGGOiyzVucmXha/zmjVOuK7ZeUFDgbpX2/L/44wNl7Cu5npOTQ05ODgUFBdStW5c6deq4//dePvmhDt+Usq9OnTrUrl2b2rVrux/72ub5uFatWpW+hVtUNdoxBE1EWgNpAI0aNWLfvn0Eel0nPvo9aVmH3OvJiXX4+Moj6dGjB5mZmZx55pmMHz+e888/P6SxG2OMMab685VnBDuUXSAKCgo4dOiQe8nJySlz3XvJzc11/1/aY+9tBQUF7uTYM0FOSEigVq1aJCQkFHvsa1t5+0vbVqtWLfbu3ctJJ7lnb0xW1e3ez0u1axHet28ftWvXDuicZ599lrSsziX67rRs2ZLExESWL1/Ozz//zIwZM0IdrjHGGGNqgEgmvb7ExcXRoEEDGjRoELE6CwsLOXz4cIkE+fDhw+Tl5RX7v7xte/fu9eucvLw895KTk1NujNWuRbh9+/bs2rWLn376iQ8++IBrr72W9u3b++gr5FBYWMjSpUtJaNW1RCJ8dHIjtmzZwv79+4mvFU/3bo7+wNWpD40xxhhjTHW0fft2zz7CNaNFePPmzQD06dOH0047jZdeeokPP/yQxnl5fHrDGcTGxpKTk0OdOnUQEd566y0WfvEACde8VaKsz28+iTVr1jB48GDeffddTj/9pBLHGGOMMcaYysXfUSOqTIuwiNQGpgA9cCTw/1LV2c597hbhJ598kqOPPpr69etz/PHH88Ybb3D11VcDMHXqVDp37sygQYNISUlh0qRJnHXWWYwaNYont7Qhr/Dv+mrHx7DmgRGRvUhjjDHGGBMS1a1FOAVH4n6ciHQDfhWRI1R1p+dBQ4YMoU+fPu710aNHs3btWnr06MG4ceNo0qQJTz75JM888wxXXnklP/30E++++y6PPFJ8UOvc/CKMMcYYY0z1VSVahEUkBtgFnK+qPzi3fQd8qapPerYIlzaO8LsLN/P4Rz/TKLERDRs0ZOXKlRQUFNCsWTOSk5NZkbavxDlHJzfyGY/1ETbGGGOMqdyqU4twJ6ApsMZj2yqgn78FjBnYnjED27vXP/poJz/++COPPzCe+Ph4n+P7fX6z9Qk2xhhjjKmuqkoi3NL5v2ezbRbQs6IFXnDBBVxwwQXu9Wb1a5GRnVds3RhjjDHGVF8+Jvyr1Lz7cZSYriQ5ORkRcS8XXHCBexy5nJwcpk+fzrZt29zHf/LJJ6xYsYLF9w0l9dGRTBtRn2kj6rP4vqEArFixgk8++cR9/LZt25g+fbpfZbrMnz+f+fPnu9etTCvTyrQyrUwr08q0Mq3M0Jd58cUXu3PA8qZXhqrTR7grsA5oqaq7nNueA5qr6iX+9BE2xhhjjDE1hz99hKtKi/AGYA/Q3WNbT+C36IRjjDHGGGOquiqRCKtqEfAqcBW4W4h7A+9GMSxjjDHGGFOFVYlE2CkFEBH5FZgBXKqq6dENyRhjjDHGVFVVZdQIVDUXGBftOIwxxhhjTPVQlVqE/eIaNSIlJSXaoRhjjDHGmChISUmpPqNGlMdGjTDGGGOMMZ6q06gRxhhjjDHGhJQlwtWMdQkxYO8DY+8B42DvAwP2PiiLdY2oZkSE6vCamuDY+8DYe8CAvQ+MQ019H9TIrhGB3iwXqU9J1e3TWHV73qpbPZFSnZ636nQtkVTdnrfqVk+kVLfnrbrVEwmV7TmraTfLdQA2AcyaNYukpCS/z+3bty9LliwJU2RWj9Vj9VSVeqrTtVg9Vo/VY/VEup7KeC3p6emMHDnStdpRVVO9j6kuifAlOCbZMMYYY4wxxtulqvqe98Zq1zXCGGOMMcYYf1SZmeXK8avrwZIlSwLqGpGcnExaWlpYgrJ6rB6rp+rUU52uxeqxeqweqyfS9VTGa0lPT6dv376u1V99HVNdukZUeNSISN1JafVYPVZP5a6nOl2L1WP1WD1WT6TrqYzXUiNHjQjUpEmTqlU9kVLdnrfqVk+kVKfnrTpdSyRVt+etutUTKdXteatu9URCVX3OanyLcHVTU8cKNMXZ+8DYe8CAvQ+MQ019H1iLcA1UnT5dmoqz94Gx94ABex8YB3sflM5ahI0xxhhjTLXjT4twlRg1QkSmAh28Np+lqtmRj8YYY4wxJjD9HpxDRnaee71Z/Vosvm9oFCMyUEUSYQBVHRLtGIwxxhhjKsIzCfa1bqLD+ggbY4wxxpgaqcokwiIyRUR+FJGvROSUaMdjjDHGGGOqtqqSCK8GXlfVQcB/gc9FpLevA5OTkxER95KSkhLBMI0xxhhjSmpWv1aZ6yY0UlJS3Dmgx41yparwqBEiMgg4GmgOZAHrgTmqml+hAgOr+/+APap6vXPdRo0wxhhjjDFuYRlHWESOF5HVwP+A04HOwEnARGCDiIypeMh+2wK0j0A9xhhjjDGmmgpo1AgR6QqMBwar6i4f++sCE0VklKp+HqIYEZG7VfUxj00tgRJZvTHGGGOMMf4KtEU4C7haVXeJSDfvnaqao6r/BRaFIjgPt4tICwAR6QicA7wd4jqMMcYYY0wNElCLsKru9lhdLSKtgKbAOlUt9DhuZ4jic3kC+ERECoB6wE2qOj/EdRhjjDHGmBokmAk1BNgE7ALqisgTqvp4aMIqTlWfwJEMG2OMMcYYExLBDJ92AOikqh2BY4CBInJ7aMIyxhhjjDEmvIJJhDe6ukCoajowBrg0JFEZY4wxxhgTZsEkwptE5J8e60VAwyDjMcYYY4wxJiKC6SN8E/CxiNwILAV6Az+FIihjjDHGGGPCrcKJsHN2juNE5HgcfYS/BGaGKjBjwqnfg3PIyM5zrzerX4vF9w2NYkTGGGOMibRAJ9ToD2Sp6nrXNlX9BfjF45iGwCBVnRWyKI0JMc8k2Nd6TRLqDwX2IcMYY0xVEWiL8DLgQxH5GZiDY6rjg0BtIAnHVMsXA2NDGKMxJoxK+1BQ0YTWPmQYY4ypKgK6WU5V84ELgVjgXRxjCGcDGTgS4yOBi1V1R4jjRERqi8hUEflVRBaLyBmhrsMY8zdLaI0xxlR3AfcRdibDDwMPi0g80AzYp6o5oQ7OSwogqnqcc3rnX0XkiDDMYmeMMcZUW9Z9yZi/BTNqhCspDnnrrzcRiQH+CZzvrHediCwFLgOeDHf9puaqCX8wmtWvVeIaK8L7uQq2PGOqq2j/XrFve4z5W1CJcAR1ApoCazy2rQL6RSccU5po/4IPtZrwB6O01yfQBNnXc5P66MjggjOmGqoJv1eMqSqCmVAjklo6/9/nsS0LaOF9YHJyMiLiXi644AJychy9NnJycpg+fTrbtm1zH//JJ5+wYsUK9/r8+fOZP3++e33FihV88skn7vVt27Yxffp0K7OUMn39gq+McTaI12JxNqxFmWV6q8qvkb9lLr5vKKmPjmTBTcfwcK99/HD7iWWWWZaqdu1WppUZ7jJ9iVScZamqz6eVaWW6XHzxxe4cMDk5mfKIqpZ7ULSJyInAAiBBVfOc2x4ATlTVU0WkNZAGkJaWRuvWraMXbA3XYWLJUfOqcqtgdbuecLLnyhj/RPtnpbp9c2dMabZv3+6ZDCc758AoJuiuESKSBHQHugE9gM6qem6w5XrZ5fw/sZTHxpgoC1VfY2Oqu2j/rFjSa8zfKpwIi8hbQHsc4wcfAh4B3gLWhia0YjYAe3Ak3K7ktyeO2exMJRLtX/AmeuyPqzH+sZ8VYyqPYFqE16jq5c4RHcYDo4FHVTU3NKH9TVWLRORV4CrgRxHpCvQGxoS6LhOc6vYL3hJ7Y4wxpvoKWR9hEakH3AO0VdUrQ1Jo8fJrA1NwdL+IA/6lqrOd+6yPsDHGGGOMcQtLH2ERaQfkqWq653ZVPQj815mUhpyzpXlcOMo2xhhjjDE1T0WGT/sIOFdE2jmX5p47fWXbxhhjjDHGVDYVSYQXq+oU4ErgC2BYaEMyxhhjjDEm/CqSCCuAqk4GVqrqO64dIpIQqsCMMcYYY4wJp2BnlvO+0+7xIMszxhhjjDEmIioyfNr5ItLb+biriPzssa8DcEuwQRljjDHGGBNuFUmEVwHTfGwXqvC4vjblpDHGGGNMzVKRRPg+Vf3Z1w4R2RBkPFHjmQT7WjfGGGOMMdVLwImwdxIsIq1wjO97JdCpImWWRUSm4uhy4eksVc0OZT3GGGOMMaZmqVDSKiJxwDnA1cCJwKfADcBrIYvMg6oOCUe5xhhjjDGm5qrIzHJPAaOB1Tj6Cl/kap0VkaLQhhc5zerXKtFH2BhjjDHGVF8VaRG+FvgdeASYrareQ6iFnIhMAY4EsoHHVHVuqOuwG+OMMcYYY2qWiowj3BqYDjwEbBWRR0Wke2jDKmY18LqqDgL+C3zuMXxbCcnJyYiIe0lJSQljaMYYY4wxprJISUlx54DJycnlHi/BNOiKSC/gKhzDpm0EuqtqYoUL9K/O/wP2qOr1HttaA2kAaWlptG7dOpwhGGOMMcaYSm779u2eyXCyqm73PiaomeVU9Q9VnYCjlfgJwOewaqURka9EJLuUZXYpp20B2gcTtzHGGGOMMSEZ6kxV84EPnEsg540o7xgRuVtVH/PY1BIokdEbY4wxxhgTiKBahCPkdhFpASAiHXEM2/Z2dEMyxhhjjDFVXUgnvwiTJ4BPRKQAqAfcpKrzoxyTMcYYY4yp4ip9i7CqPqGqJ6rqYFXtp6rWGlwGGyXDgL0PjL0HjIO9DwzY+6AsQY0aUVnYqBF/ExGqw2tqgmPvA2PvAQP2PjAONfV9EPZRI6qDcH1K6vfgHDpMnOVeWp86Niz1REukPl1aPZVbdXreqtO1RFJ1e96qWz2RUt2et+pWTyRU1eesurQItwG2AixZsoSkpCS/z01OTiYtLS3kMQ146Nti62kvXhGWeryF63qsHqunutdTna7F6rF6rB6rJ9L1VMZrSU9Pp2/fvq7Vtqq6zfuY6pIIH4tj2mdjjDHGGGO89VHVpd4ba3zXCGOMMcYYUzNVlxbheOAo5+puoCiK4RhjjDHGmOiLAZo7H//pnACumGqRCBtjjDHGGBMo6xphjDHGGGNqJEuEqxER6S8if4nIuGjHYiJPROJFZIKIzBOR+SLyi4icFu24TGSJyM0i8p2IfC8iK0Tk5mjHZKJDRLqKSL6IDIl2LCbyRGSq8++B51I/2nFVNlVhimXjBxE5D/gHsC/asZioSQZuBXqr6j4RGQrMFJHuqhr+MW1MZfFPYKiq7hKRrsBqEVmoqouiHZiJuPuBvGgHYaJHVYdEO4bKzlqEq4/fVHU0cCDagZioOQD8V1X3AajqHCAXOCGqUZlIG6uquwBUdT2wF+gQ1YhMxIlIfyAbxw3kxphSWCJcTfgaJNrULKqaqapvu9ZFRIBa2B/CGkVV/3A9FpELcCRDs6MXkYmSyc7F1GAiMkVEfhSRr0TklGjHUxlZ1whjqq/BwGbgh2gHYiJLRI4G/g+oB1ykqlnRjchEkoiMAFaq6jbH52FTQ60GvlfV35zfEMwVkZNUdVmU46pUrEXYmGpIRGoDDwPjVNXG1a5hVHWFqvYEzsXRT3xAlEMyESIiMcA9wCPRjsVEl6r+P1X9zfn4N2AWcG10o6p8LBE2pppxdol4BXhKVZdEOx4TPc7pRGcBt0Q7FhMxo4FvVHVPtAMxlc4WoH20g6hsLBE2pvp5Elikqh+ISIKItIt2QCYyRKSps1+wp4M4ukiYmmEQcJZruCwgCXhaRGZGNywTaSJyt9emlsD2aMRSmVkibEw1IiL34Oj7P9U5XmRn4KroRmUiqAHwHxGpCyAizXF0j/gumkGZyFHVa1X1RFUd4hw6Kx2YoKrnRDk0E3m3i0gLABHpCJwDvF32KTWPTbFcTYhIXxwtgb1x/OJbparnRzUoE1Ei0g1Y62PXZFVNiXA4JgqcfcP/DZyOY/zYROAT4H7rK16zOPuFPwYcB6wBPlbV+6MblYkkEbkTOA8owPGt0DOeIwsZB0uEjTHGGGNMjWRdI4wxxhhjTI1kibAxxhhjjKmRLBE2xhhjjDE1kiXCxhhjjDGmRrJE2BhjjDHG1EiWCBtjjDHGmBrJEmFjjDHGGFMjWSJsjDHGGGNqJEuEjTHGGGNMjWSJsDHGGGOMqZEsETbGGGOMMTWSJcLGmBpJRL4QkUIRWSMi4z22vy8ib3msny0iq5zL4DDF8rSIpItISjjKDyCO70QkV0SGBHheXRH5fyJSP0xxdRKRf4lIXDjKN8bUXJYIG2NqJFU9C1gGfK6qrwCISAxwKnCax3GfAR8Do1V1fphimQB8HY6ySyMi40RknlccpwHpFShuKjBPVbNDEFoJqroR2Ac8Eo7yjTE1lyXCxpiabA5wusf6scA8oLWIHOmxvRewPIJxVRkicibQRFW/CnNVU4DzvF4XY4wJiiXCxpiabA5wjIg0d66fDrwMbAKGAohIIyBLVVVEhju7D8wRkR9FZILzmI4islJE8kTkC+e2U0VkrYisFpEuIpIgIk+IyK8iMl9E3hGRxqUFJiKtReQjZz0/icgkZ4s1IvKuiGQ5t/2fiCxxxtTY4/weznMXishnzq4LWSLynoicDUwEeovIPBH5zqv63iLysYiscHYViS/jORwDuM8XketFJFVEpjrXjxeRZSKS6nGMqwvGBBH5xNk95VkRaSsi00VkubNed1cIVS3E8SFldBmxGGNMQCwRNsbUZAuAXP7uCnEi8COOBHmoc9spwFzn4wbAP1V1KHAycK6IDFLVTcAZgAB3Aqjq98BCYISq/gX8P2AAcLKqDgZ2Aq/5CsqZ8H4OrFfVQTi6awwDbnSWPQZHt47BwGWq2td56tXO8+Oc53+tqgOBi51lLFPVS5zdPR51rg9xdonw1Be4AEcLeS/gnDKewwHAX64VVX0JR1cJ1/ovwATPEzy6YPRW1fOc9VwKPABcDvQBevqodx3Qv4xYjDEmIJYIG2NqLFU9jCMZHioiCc5tucC3wMnOltDTnOvgSD4fFpGfcCTHPYDjnOelAV8B4wGcrbONVDXVmdiOB15W1TxnWdNwfNXfxEdoA3Akg895xPke8E+v4z5X1Xzn4yVAV+fj44AuOBNtVT0E/F8AT81H6lDgvOauZRybBOwPoGxPn3rEtxZHYl7gbP1dDnTzOn6/sz5jjAkJuwPXGFPTzQFuAk4CfnJu+w6oCxwPtFXVrc7ts3Dc1Dba2VViKlDPo6xXgTdF5F5gLPC2c3szoA5wp4hc49wWC2wBWgF7vGJq5/z/XRFxbatHycaLvR6Pc4EE5+Nk5/+7PfZneF94GbJKKdeXGEADKLu0egq81vN91FuENeAYY0LIEmFjTE03B3gMuAF4GEBV94jIUuAqYCuAsx9xV+BGVXUlfrW8yvoSOISjW8HZwAjn9gzn9knObgk4y2xK8eTPZYvz/7NVdb/zWMGRUPsjzfl/C/4eBcLfcwO1E0eXEW+1PR77avWuiEbO+owxJiTsk7UxpqZbjqPl9GRgqcf2OThadec41zNxJLQng/smukGeBTm/0n8DR3/gxa5uC6paBLwCXOW6AUxEuuPoluHLIuB34FqPbZcD//Pzmn7F0W/3GmdddXAk557240xgnTfxDfSzbG+/Ax19bO8tInHOBP68CpbtrRPFXyNjjAmKJcLGmBrN2br7LTDXmbC6zMHxVfw853FFOG7o+oeILABeANYD40Tkeo/zXsPRj/VVr6ruwdEPdoGIfA88C1ykqoUi8jQw3FnWXc66RgEDRGSBc7zfM3B04UBEXgZ6AxNF5DxxTAgyDhguIpOdfXvPBoaJyKL/394dq0QMBAEY/qe2Eiws9QUE38DXOLjCSm20FKtDC3tfwUYLQTi0Eks7Gx/EzlaYK3YP4yEokjTZ/+sSQibTDZPZXeAGeOL7CMMz8Flz2QHeImJev/2q7vYw63xXN8euW8qCwlXrlCL5AXgFNpe7U/wQZzWfi07cZc4B7NV4ktSL+PrDJ0kai4jYyMz3zvUpsJuZk57jBGUk5DIzX+q9c2ArM/d7jDOl7Lhx8OvDkvRHdoQlaZyu6/gFEbFG6WY/9h2kdtQnwDSGO2J5m9IxPhni/ZLaZUdYkkaojkscAh+UHSfuM3PwI4oj4gw4oiyWu8vM46FjStJ/WQhLkiSpSY5GSJIkqUkWwpIkSWqShbAkSZKaZCEsSZKkJlkIS5IkqUkWwpIkSWrSArcRQ/K0U3KoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import Image\n",
    "Image('spectrum.png') "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}