tutorials/tutorial_03/pysnack_1_03.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set up environment\n",
    "Let's load the basic modules and read the tables from disk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import astropy.units as u\n",
    "from astropy.table import Table, QTable\n",
    "from astropy.coordinates import SkyCoord\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now load the two tables we generated in the last session. `data1` contains the Gaia query in a rectangle around our region of interest, excluding negative parallaxes. `cluster1` is our selection on the previous table to select sources with parallaxes between 5.0 and 5.7 mas."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "⚠️ If you want to use the \"official\" `data1` and `cluster1` tables instead of the ones created by yourself do:\n",
    "- open a terminal\n",
    "- Navigate to the data folder  \n",
    "`cd pysnacks_1/tutorials/data`\n",
    "- Delete the files `data1.ecsv` and `cluster1.ecsv`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "for my_file in ['data1', 'cluster1']:\n",
    "    if not os.path.isfile(f'../data/{my_file}.ecsv'):\n",
    "        print(f'Using official {my_file} file')\n",
    "        os.system(f'unzip ../data/backup/{my_file}.ecsv.zip')\n",
    "        os.system(f'mv {my_file}.ecsv ../data')\n",
    "    else:\n",
    "        print(f'Doing nothing because ../data/{my_file}.ecsv already exists')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data1 = Table.read('../data/data1.ecsv')\n",
    "cluster1 = Table.read('../data/cluster1.ecsv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ADQL searches\n",
    "\n",
    "Astronomical Data Query Language (ADQL) is a language based on the Structured Query Language (SQL) but developed for astronomy to query catalogs and tables following Virtual Observatory standards. Full specifications can be found [here](https://www.ivoa.net/documents/ADQL/20180112/PR-ADQL-2.1-20180112.html).\n",
    "\n",
    "- This [ADQL Cheat Sheet](http://simbad.u-strasbg.fr/simbad/tap/help/adqlHelp.html) is very useful.\n",
    "- Many ready-to-use ADQL examples on Gaia data: https://gaia.ari.uni-heidelberg.de/tap.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astroquery.gaia import Gaia"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gaia provides a large amount of tables. If you want to see all available tables (including many cross-matches with other catalogs like SDSS, WISE, UCAC4, Hipparcos, etc) you can run `Gaia.load_tables(only_names=True)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tables = Gaia.load_tables(only_names=True)\n",
    "for table in tables:\n",
    "    print(table.name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this tutorial we are going to use table the full Gaia EDR3 catalog:  \n",
    "`gaiaedr3.gaia_source`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic ADQL query\n",
    "\n",
    "There are 3 main stages for a query:\n",
    "1. Write the query as a text string\n",
    "2. Launch the job with the function `Gaia.launch_job()`\n",
    "3. Convert the job into an astropy Table with `job.get_results()`.\n",
    "\n",
    "The query will be written as a python string, which we will send to the Gaia server as a job. The main ingredients to write a query are these statements:\n",
    "- `SELECT`: this statement is used to select data from a database.\n",
    "- `TOP`: to limit the query to a certain number of rows. Be careful not to download too many rows!\n",
    "- `FROM`: to indicate the table we want to query.\n",
    "- `WHERE`: to include a filter to select specific data. `AND` to include additional filters.\n",
    "\n",
    "Although it is not a formal requirement, these statements are usually included in CAPITAL letters, while our values are written in lowercase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query1 = \"\"\"\n",
    "SELECT\n",
    "TOP 10\n",
    "    ra, dec, pmra, pmdec, parallax\n",
    "FROM\n",
    "    gaiaedr3.gaia_source\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that dividing in lines and adding spaces is just to make it nice. 4-space identation is not relevant. The query could be written as:\n",
    "```python\n",
    "query1 = \"SELECT TOP 10 ra, dec, pmra, pmdec, parallax FROM gaiaedr3.gaia_source\"\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job1 = Gaia.launch_job(query1)   # Here is where the data are downloaded!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "table1 = job1.get_results()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "table1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's add a filter. This will be filtered in the input catalog and we only retrieve the resulting filtered table, which makes it fast:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query2 = \"\"\"\n",
    "SELECT\n",
    "TOP 10\n",
    "    ra, dec, pmra, pmdec, parallax\n",
    "FROM\n",
    "    gaiaedr3.gaia_source\n",
    "WHERE\n",
    "    parallax > 250\n",
    "\"\"\"\n",
    "job2 = Gaia.launch_job(query2)\n",
    "table2 = job2.get_results()\n",
    "table2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And adding additional filters is also easy:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query3 = \"\"\"\n",
    "SELECT\n",
    "TOP 10\n",
    "    ra, dec, pmra, pmdec, parallax, phot_g_mean_mag\n",
    "FROM\n",
    "    gaiaedr3.gaia_source\n",
    "WHERE\n",
    "    parallax > 10\n",
    "    AND dec BETWEEN 55 AND 60\n",
    "    AND pmra > 100\n",
    "\"\"\"\n",
    "job3 = Gaia.launch_job(query3)\n",
    "table3 = job3.get_results()\n",
    "table3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do we know which columns are available, their units and the description? We can get the table metadata with `Gaia.load_table()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "table_metadata = Gaia.load_table('gaiaedr3.gaia_source')\n",
    "for column in table_metadata.columns:\n",
    "    print(f\"{column.name:31}: {column.unit:16} :  {column.description}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ⛏ Exercise 3.1\n",
    "Construct and run a query named `job3_1` using the following information:\n",
    "- Select 10 rows.\n",
    "- Query table `gaiaedr3.gaia_source`.\n",
    "- Select columns: right ascension, declination and phot_g_mean_mag.\n",
    "- Use filters to select stars with declination between 50 and 80 deg and magnitude brighter than 10."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can rename columns and give a name to the table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query4 = \"\"\"\n",
    "SELECT\n",
    "TOP 10\n",
    "    ra AS right_ascension, dec AS declination, pmra, pmdec, parallax AS plx, phot_g_mean_mag\n",
    "FROM\n",
    "    gaiaedr3.gaia_source AS gaia\n",
    "WHERE\n",
    "    parallax > 10\n",
    "    AND dec > 60\n",
    "    AND pmra > 100\n",
    "\"\"\"\n",
    "job4 = Gaia.launch_job(query4)\n",
    "table4 = job4.get_results()\n",
    "table4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can create new columns as combinations of other ones including complex operations.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query5 = \"\"\"\n",
    "SELECT\n",
    "TOP 10\n",
    "    ra, dec, pmra, pmdec, pmra_error, pmdec_error,\n",
    "    sqrt(power(pmra_error, 2) + power(pmdec_error, 2)) AS pm_error\n",
    "FROM\n",
    "    gaiaedr3.gaia_source AS gaia\n",
    "WHERE\n",
    "    parallax > 10\n",
    "    AND dec > 60\n",
    "    AND pmra > 100\n",
    "    AND sqrt(power(pmra_error, 2) + power(pmdec_error, 2)) < 0.1\n",
    "\"\"\"\n",
    "job5 = Gaia.launch_job(query5)\n",
    "table5 = job5.get_results()\n",
    "table5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computing the distance in kpc is as easy as creating a new column using 1/parallax[mas]. Note: you cannot call the column \"distance\" because that is a reserved ADQL statement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query6 = \"\"\"\n",
    "SELECT\n",
    "TOP 10\n",
    "    ra, dec, pmra, pmdec, parallax, 1/parallax AS distance_kpc\n",
    "FROM\n",
    "    gaiaedr3.gaia_source AS gaia\n",
    "WHERE\n",
    "    parallax > 1\n",
    "    AND dec > 60\n",
    "    AND pmra > 100\n",
    "\"\"\"\n",
    "job6 = Gaia.launch_job(query6)\n",
    "table6 = job6.get_results()\n",
    "table6['distance_kpc'].unit = u.kpc   # It is always good to define the units of new columns as soon as they are created.\n",
    "table6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, as we did in previous session, we can directly compute the absolute magnitude and create a new column:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query7 = \"\"\"\n",
    "SELECT\n",
    "TOP 10\n",
    "    ra, dec, pmra, pmdec, parallax, phot_g_mean_mag,\n",
    "    gaia.phot_g_mean_mag + 5*log10(gaia.parallax) - 10 AS g_mag_abs\n",
    "FROM\n",
    "    gaiaedr3.gaia_source AS gaia\n",
    "WHERE\n",
    "    parallax > 10\n",
    "    AND dec > 60\n",
    "    AND pmra > 100\n",
    "\"\"\"\n",
    "job7 = Gaia.launch_job(query7)\n",
    "table7 = job7.get_results()\n",
    "table7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Select by region"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query8 = \"\"\"\n",
    "SELECT\n",
    "TOP 10\n",
    "    ra, dec, pmra, pmdec, parallax, phot_g_mean_mag\n",
    "FROM\n",
    "    gaiaedr3.gaia_source AS gaia\n",
    "WHERE\n",
    "    CONTAINS(POINT(ra, dec), CIRCLE(15, 60, 0.1)) = 1\n",
    "\"\"\"\n",
    "job8 = Gaia.launch_job(query8)\n",
    "table8 = job8.get_results()\n",
    "table8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a query programatically\n",
    "\n",
    "This is useful when we want to separate the query template from the actual variables we want to use, for example by defining the columns or the table name as independent variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_columns = 'ra, dec, phot_g_mean_mag'\n",
    "my_table = 'gaiaedr3.gaia_source'\n",
    "min_dec = 50\n",
    "max_dec = 80\n",
    "max_mag = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query9 = f\"\"\"\n",
    "SELECT\n",
    "TOP 10\n",
    "    {my_columns}\n",
    "FROM\n",
    "    {my_table}\n",
    "WHERE\n",
    "    dec BETWEEN {min_dec} AND {max_dec}\n",
    "    AND phot_g_mean_mag < {max_mag}\n",
    "\"\"\"\n",
    "print(query9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job9 = Gaia.launch_job(query9)\n",
    "table9 = job9.get_results()\n",
    "table9"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can verify that `query9` is exactly the same as `query3_1` and the `table9` is the same as `table3_1`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(query9, query3_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(table9)\n",
    "print(table3_1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ⛏ Exercise 3.2\n",
    "Find all the stars around Vega in a radius of 50 arcsec. (Don't use `TOP`).\n",
    "- Use `SkyCoord.from_name()` to get the coordinates of the star Vega and store them in `vega_coord`. You can get the R.A and Dec in degress with `vega_coord.ra.deg` and `vega_coord.dec.deg`, respectively.\n",
    "- Run an ADQL query to find stars within 50 arcsec of that position to obtain columns 'ra, dec, phot_g_mean_mag, distance_kpc'. Create a new column `distance_kpc` by transforming the parallax.\n",
    "- Is Vega included in the resulting table? Read more [I can not find data for a bright star. Why is that?](https://www.cosmos.esa.int/web/gaia/faqs#BrightStar). Also, read about the VEGAMAG system [here](https://gea.esac.esa.int/archive/documentation/GDR2/Data_processing/chap_cu5pho/sec_cu5pho_calibr/ssec_cu5pho_calibr_extern.html#SSS6):\n",
    "> The VEGAMAG system is defined so that Vega’s colours are all zero: this is equivalent to normalizing all the observed fluxes to the flux of Vega."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 🌪 Additional fun\n",
    "Repeat the query around Vega but with a radius of 100 arcsec. Select 50 rows. Make a plot of the R.A. and Dec. of the resulting stars."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Query the region of interest around our star cluster\n",
    "Now we will repeat the search we did in sessions 1 and 2: all Gaia stars in a rectangle of 3 deg x 2 deg centred at R.A. 130.025 deg and Dec. 19.98333 deg. Then, we created `data1` selecting only positive parallaxes and we added the distance and absolute magnetiude columns. We will do all steps in a single go with ADQL."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_table = 'gaiaedr3.gaia_source'\n",
    "min_ra = 130.025 - 3/2\n",
    "max_ra = 130.025 + 3/2\n",
    "min_dec = 19.98333 - 2/2\n",
    "max_dec = 19.98333 + 2/2\n",
    "min_parallax = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We defined gaia_columns as a python list in tutorial 1. Here we need a single string\n",
    "gaia_columns = ['ra','dec','parallax','parallax_error','pmra','pmra_error','pmdec','pmdec_error','astrometric_matched_transits','phot_g_mean_flux','phot_g_mean_mag','phot_bp_mean_mag','phot_rp_mean_mag','bp_rp','bp_g','g_rp','dr2_radial_velocity','dr2_radial_velocity_error','dr2_rv_template_teff']\n",
    "my_columns = ', '.join(gaia_columns)  # This creates a comma-separated string from a list\n",
    "my_columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apart from `my_columns` we will create two new columns: `distance_kpc` and `g_mag_abs`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query_final = f\"\"\"\n",
    "SELECT\n",
    "TOP 35000\n",
    "    {my_columns},\n",
    "    1/parallax AS distance_kpc,\n",
    "    phot_g_mean_mag + 5*log10(parallax) - 10 AS g_mag_abs\n",
    "FROM\n",
    "    {my_table}\n",
    "WHERE\n",
    "    ra BETWEEN {min_ra} AND {max_ra}\n",
    "    AND dec BETWEEN {min_dec} AND {max_dec}\n",
    "    AND parallax > {min_parallax}\n",
    "\"\"\"\n",
    "print(query_final)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job_final = Gaia.launch_job(query_final)\n",
    "data1_adql = job_final.get_results()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data1_adql"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Are the table obtained with ADQL exactly the same as the one obtained with astroquery and modified as a table? Let's see if the number of entries and some average values match exacty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(len(data1),               len(data1_adql))                   # The number of entries\n",
    "print(data1['ra'].mean(),       data1_adql['ra'].mean())           # Same mean R.A.?\n",
    "print(data1['parallax'].mean(), data1_adql['parallax'].mean())     # Same mean parallax?\n",
    "print(data1['distance'].mean(), data1_adql['distance_kpc'].mean()) # Same mean computed distance?\n",
    "print(data1['Mg'].mean(),       data1_adql['g_mag_abs'].mean())    # Same mean computed absolute magnitude?\n",
    "print('Yes!')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advanced plotting with multiple panels\n",
    "\n",
    "For this section we have two main objectives:\n",
    "- Describe how to manage plots with multiple panels.\n",
    "- Write functions to simplify the plot management."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiple subplots in a single figure\n",
    "\n",
    "You can use  \n",
    "`fig, ax = plt.subplots(ncols=XX, nrows=YY, figsize=(14,6))`  \n",
    "to create XX times YY panels in a single plot. Now `ax` will be an array of size (XX, YY).\n",
    "\n",
    "Here are some example of how `ax` looks like when creating multiple panels:\n",
    "\n",
    "- When you use ncols=2, nrows=1:\n",
    "  - `ax[0]` left panel\n",
    "  - `ax[1]` right panel\n",
    "\n",
    "- When you use ncols=1, nrows=2:\n",
    "  - `ax[0]` top panel\n",
    "  - `ax[1]` bottom panel\n",
    "  \n",
    "- When you use ncols=2, nrows=3:\n",
    "  - `ax[0, 0]` top    left  panel\n",
    "  - `ax[1, 0]` middle left  panel\n",
    "  - `ax[2, 0]` bottom left  panel\n",
    "  - `ax[0, 1]` top    right panel\n",
    "  - `ax[1, 1]` middle right panel\n",
    "  - `ax[2, 1]` bottom right panel\n",
    "\n",
    "In summary, the general nomenclature is [row, column]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "col1 = 'parallax'\n",
    "col2 = 'distance'\n",
    "\n",
    "fig, ax = plt.subplots(ncols=2, figsize=(14,6))# We can select ncols, nrows, or both.\n",
    "\n",
    "ax[0].hist(data1[col1], bins=np.arange(0, 8, 0.05), label=data1[col1].description)   # histogram of column 1\n",
    "ax[1].hist(data1[col2], bins=np.arange(0, 4, 0.02), label=data1[col2].description)   # histogram of column 2\n",
    "\n",
    "ax[0].legend()                                 # Activate the panel legend, which will show the \"label\" we specified\n",
    "ax[1].legend()                                 # Activate the panel legend, which will show the \"label\" we specified\n",
    "\n",
    "description1 = data1[col1].description         # Obtain column 1 description\n",
    "unit1 = data1[col1].unit                       # Obtain column 1 units\n",
    "description2 = data1[col2].description         # Obtain column 2 description\n",
    "unit2 = data1[col2].unit                       # Obtain column 2 units\n",
    "\n",
    "ax[0].set_xlabel(f\"{description1} [{unit1}]\")  # Set the axis label in the form \"Variable description [units]\"\n",
    "ax[1].set_xlabel(f\"{description2} [{unit2}]\")  # Set the axis label in the form \"Variable description [units]\"\n",
    "\n",
    "ax[0].set_ylabel('Number of stars')\n",
    "ax[1].set_ylabel('Number of stars');"
   ]
  },
  {
   "attachments": {
    "067abd68-8a10-4b42-b130-7252cc7e0013.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmAAAAJICAIAAABJ9s/cAAAAA3NCSVQICAjb4U/gAAAgAElEQVR4Xuy9CfxN1f7/fw1FGSqFCCGueZ7SlUhERKlEKkVXpUE/pXSr23DvzYcmKg3K7bpNSgOaSMauciU0UJIQcpGhJGT6P+9nff/nce45+7M/+5yzzz57eJ3H53Ee+7P22mu938+13ue9pr1WkcOHD/9OHxEQAREQAREQgf8lUFRAREAEREAEREAEkgnIQSYzUYgIiIAIiIAI/E4OUpVABERABERABCwIyEFaQFGQCIiACIiACMhBqg6IgAiIgAiIgAWB4hZhOQo64YQTqlevnqPMla0I5IzAunXrtm7dmrPsC85YJlkwG90JLYF4e/SRg8Q7Ll68OLTUpZgIFECgZcuWBdzJcbBMMscFoOxzQSDeHjXEmosSUJ4iIAIiIAK+JyAH6fsikoAiIAIiIAK5ICAHmQvqylMEREAERMD3BHw0B+l7VpETcP/+/Rs2bNi7d2/kNM+mwiVLlqxSpcoRRxyRzUyUtghYE5BRJ3Cxt0c5SOtqpFAI4B3LlCnDSo0iRYoIiCsE2Pp427ZtgK1Ro4YrCSoREUiJgIw6Hleh9qgh1pRqV7Qi03c8/vjj5R1dLHVgglSdcheRKqmUCMio43EVao9ykCnVrshFlnd0vciF1HWkSjAlAqqBCT7Shp4cpA0c3co6gWLFijVt2rRhw4bnnnvuzp07ye+HH3648MILC8q4dOnSCbemTJmyYsWKguJnO/ycc84xYmc7I6UvAkEhYIy6QYMGTZo0efjhhw8dOuRQ8iuuuOK1114j8lVXXeXcqOfOnfvRRx+ZLJ566ql//vOfDrNzEk1zkE4oKc7vqo94J20Ka/O6F/TsUUcdtWzZMu4OGDBg3Lhxd9xxR+XKlY2RFPRIQjgOskePHvXr13cY391o7777rrsJKjUR8IxAto16y5Ytl1xyyU8//XTvvfempNSzzz7rPD4OknbzaaedxiPXXHON8wedxFQP0gklxck6gbZt227cuJFs1q5dS4eSi19//bVPnz6NGze++OKL27RpE9tlCSdKy/TUU0/dvHkzLcdp06YNHz6cbujq1atjUr711ls80qxZs7POOotohM+bN484fAjctWsXIaNHj27UqBFJjRgxgn95vGvXri1atDj99NO//vprQmjP3njjjRhezZo1jc/etGlT+/btTZf3ww8/JA4rmH788UcuaCkjNp8xY8bwL1rUq1fvj3/8I+3oLl267NmzJyabLkQgIgQqVKgwfvz4xx9/nLUwBw8exE5btWqFRT/99NOGQIINxrB06NDB2DueL97eiZBg2hgavcZHHnkEq8Qk77nnngcffJBoNLv5iSCv888/f8eOHYSQ5m233da6devf//73xnidlIJ6kE4oeRQnvkFn0+vySBoPs8F4Zs2aNWjQoPg8n3jiieOOO+7zzz//8ssvqf3m1u7du6n3f/vb32699dZnnnnmzjvv7NmzJz3IhFHZdu3aLVy4kLkWmqIY4UMPPYTZ0EP9wx/+8Msvv7Cw+7333qPr+e9///voo4/evn07iQ8ePBhLq127NoFDhgyZPXs2gXjEf/3rX/hLciGLl1566eyzz8ZiERj/HZP2008/fe6553iQHwIc8xlnnIHkq1atevnllxESN//6669feumlHhJVVk4JODc65zGd5h2BeDQuGWKlKzl16tRjjjnmk08+2bdvH2ZIqxGzSrDBZB7J9p5s2vQa8aO33HILj/MzYhK5/PLLH3vsMSzxz3/+M/1X02w9cODAokWLGPUh5IMPPkjOLjlEDjKZiUK8I0DXCudHM5CuW+fOneMzxjMNHTqUELpltATNrSOPPBJ3yDXxZ86cWZCgrGWn34l7++2338wLFdjksGHD+vfv37t3b15DxDyuvPJKvCMplCtXDq9JZ/Siiy4yCWLD5uK8884rWrQo47emG0r7d+DAgbxJRnjMZxMTUWmolipVimvSp32KQyVfEwdRUbAgUV0Pp1PLyznMAxUvXpxmOO4fFAhA+KuvvornJseRI0dOmDCBOI8++igu33UZlGCuCCSMmvqhnU2rERrvv/8+jV0zEsOgK83HBBu0JJZs78mmnfwg6bMyAO/ILeZuYnaNbRKSkj1qiDUZr0K8I2DmINk+H09GDy8+Y2NXCaLwfr1Zg8ePO+3BggS94YYbrr/++i+++ILBHPNOBeOo9Cbxx3RAabqSePxaPhq5xx57LMMy5vPVV1+ZlEuUKGEujDCMr86fP/+kk0667LLL4tcCWIoae9Ze1IJUyCR8zpw5KGIGqfLy8jp16sTvEd9ckyzLHyZNmrR8+fLp06fTV6Y3nEleelYEbAh899131H/GWrERunTGvtasWUMPMsEGLRNJtvdk07Z80DLQmGRK9uiCg6RlylwOjWWzCTotVroCDFXxbQZ/kZUWa61aterUqTNjxgxL0RUYZQKMvdCVYRSUzlmMA2Mp9Hj4lx90XF1BfOgtmTnF+Ag0IXFjhEycONGEM8VILWUSglqKg8Q+//73v5thUmps2bJl6fBNnjyZyNjtZ599ZpkdjhxTZ2aR0eAlS5bE4uA4GSwiNUaE3nzzTWYxLR/PSSBDWzSiyZpvhOSCkL59+/JjgcpYJYNOORFMmaZNgG5i7C/tRDx4kEPcGP+kqUpjlIGKJ5980hj4N998g6Uk2KBDeZJNO/kXgN8TRkrMROPzzz9vupIO00+I5oKDJEW1WNOjr6diBFg7w3oZejaxEDo3GBiDq6NGjeKbSm+Ji9/6Bx54gMfjF+kwV8+4Co6KEw3NU0xCMFRLFvRZu3XrxnocRkFxljTszKz+iy++yKgjEVhWgwuxzIv1csQnL+YUzfCvida8eXNW9DD/zwQkK9SJYPm4N4H8GPHTwzgS6yPIkcHhSpUqccE3U0FcsBiqatWqRhhGm83aKG9kUy5RIGDmTbAjlshRFe+++260xi6YqsBSMMOrr76a4Z9kG3QCJ9m0eUOMVimGGb/0hpYxa4L43aDPyjSkk5Qt4xSxHB2yjFpQID1IBnNiv0R0E/kdwRqZAWLh0MqVK+k+8uztt9/ON+0INGTJYnJq/FpF/DxIv60CYKSRpZjJJeVNCEN/tDdZUIPnY3iQVicTEt5kne1cEsC6W/N5kZRXZfCFDOEwqEU7IPamJs1qBnWuu+46DNAsGqIrzKucF1xwQbzKeFbjXGmg0GnONo3cpu/c6JzHzLZGNpIkzEHGS2LmI3Nr1Nkmk176NvbowiId02Llm3YBSwEtW6xM/BjRs9ditZmdtrmVHlA95QEBRiw7duyIj6QNx+BMaLxjttHhHcmCoWDWDTF8WrFiRZqqpsFKILewwfXr1xsxWPJg4sdLhRXzIcS3JzknM5SNJzNRSOYEXBhiXbBgAfMxrJtnkQVLGJJlSuikJmx0RFsVO+RDczX5WYVElgBTC4woMB3I4jcGRSPLISXFmdoxM7JcsG6Q4Sx6kGYilu9evXqRGiEMZbNSl7USLN5hZDilLBRZBKJDwIUeZIYt1iA2V92qHzbjIW5loXQiRYDxGzqOqMwcD5uYMM3Dqym8iMn0arVq1cwqJCaHCGFCiPdAaNSyqC9SiKSsCDgnkKmDpKHKEnka+6bFynSoabGyqj6+xYqt8hYasyNqsTovG4cxbQaXbG45TNzJUmyHSSmaIZD5rL8NSd7LTliCy+EhsbenYw+y1wEfm3R0K8QEZNTxhWtvj5k6SLVYQ2xILJDh8EKdeOViEWONIAWsi2kqKRFwTkBGneAd7e0xUwepFqvzqhm4mKzmYBGH5obdLThzgrm7aSo1EXBIQEadAMreHjN1kA5LJaDRbJZTB1SjlMRmGwveJc98nJZMXUkkJeEVWQQCTcDhj0+q6xiMUQeajJfCy0F6SVt5iYAIiECACUStpSsHGeDKKtFFQASCRcBhvzBYSoVYWhfegwwxHakmAiIgAiIQWQIB60FGp/0VtaGMyFqgFBcBvxHQj0+sRALmIP1WkySPCIiACIgABELpVuUgVbdFQAR8SiCUv7lpsE51qWoaWegRSwJykJZYshWY24puk7sfTh7PFnSlG3kC0Zma8U9Rh6NxIwfptEaFo7wdamvjSh2moGgiIAIiEHQCcpD/M3SujlTQK7TkDwQBmy6dGmeBKEFvhMx5t0QO0v2Czkah2vyguK+AbYrZ0M42Q90UAY8IOK/bNl5cjWyPSss2G7d+MOUg/wezTb23LQ4f3QyBCj6iKVFCR0AGEroizaJCcpBpwnXe2EwzAz3mCQGVoyeY3clEheUOR6XimIAcpGNUiugeAf3Succylym5Uo7q0uWyCJW3LYHQOsgQW12IVbOtq7opAgEg4EqjIQB6RkPE0DrIaBSftBQBERCBgBEIUBsiig5SPTBv7EmcveHsz1wC9CNYEEBXKrAriRQkocKzTSCKDjLbTJW+iwTif1+0gN5FsNlISs4gG1TdSjNkTRZvfg3kIN2pfvppcIdj0pbH8cl6YOHyx26Vo9IRAS8JZOkXWA7Sy0IMYV6ueJQsVe6CcDvMzgN/XJCECo8CAVUw/5dygB2kw585/5eBHyQMIkz9vvih5kgGEciQQLYNOZP0A+wgMywVPe46gdx6WZvcvZmucJ1nmBK0KZ0wqSldQkZADjLrBaqfhqwjVgaRJCDLyqTY4+l50IIMaGEVzQSxnhUBERABERCBsBJQDzKsJSu9REAERCB4BHzV15SDDF4FksR+IODxCJUfVJYMYSVg75Ns7trcCgcrOchwlKO0sCMQejO2U173REAE0iUgB5kuOT0nAiLgVwJBbBIFUeYclr83uOQgc1jEyjokBBJs1YM1gSEBJzVEwN8EtIrV3+Uj6URABERABHJEQA4yR+CVrQiIgAiIgL8JyEH6u3wknQiIgAiIQI4IyEHmCLyyFQEREAER8DcBOUh/l4+kEwEREAERyBEBOcgcgVe2IiACIiAC/iYgB+nv8pF0ASTAWx+xv5yIf/DgwWbNmvXo0YPct2/f3rlz59q1a/O9Y8cOI8/IkSNr1apVp06dGTNm5ERCZSoCgSAgBxmIYpKQIpACgbFjx9arV888kJeX16lTp1WrVvHNNYErVqyYNGnS8uXLp0+fPmTIELxpCkkrqghEiYA7DlIt1ijVGenqawIbNmx45513rrrqKiPl1KlTBwwYwDXfU6ZM4YKQvn37lihRokaNGvQjFy1a5Gt9JJwI5I6AOw5SLdbclaByFoH/IXDTTTeNHj26aNH/M+3NmzdXqlSJGHxv2bKFi40bN1atWtU8U6VKFf4VQREQAUsCLjhItVgtySowZATiZxa92QcyDYBvv/12hQoVWrRoYfPs4cOH4+8WKVIkIfL48eNb5n+2bt1qk45uiUDoCbjgINViDX0tkYJBIbBgwYJp06ZVr16dQdTZs2dfeumlFStW3LRpE/Lzje/kgl7j+vXrjUa0bitXrpyg3eDBgxfnf8qXLx8UxSWnCGSDQKYOMvMWq5qr2ShXpRlNAixPxeetXbuWZThnnnnmCy+80LNnz4kTJ0KD7169enFBCHf37du3Zs0aFu+0bt06mqyktQgUSiDT0zxMi/Xdd9/du3fvzz//HGuxMuHhsMVKc5UPgjKkU6i4iiACIpASgREjRvTp02fChAnVqlWbPHkyzzZo0ICQ+vXrFy9efNy4ccWKFUspQUUWgegQyLQHqRZrdOqKNA0QgQ4dOjC6g8DHH3/8rFmz6CnyXa5cOaPCHXfcsXr16pUrV3br1i1ASklUEfCYQKY9yGRx1WJNZqIQERABERCBwBFwzUHSYuWD/qbFmgCCFiufwNGRwCIgAiIgApElkOkQa2TBSXEREAEREIFwE5CDDHf5SjsREAEREIE0CchBpglOj4mACIiACISbgBxkuMtX2omACIiACKRJQA4yTXB6TAREQAREINwE5CDDXb7STgREQAREIE0CcpBpgtNjIiACIiAC4SYgBxnu8pV2IiACIiACaRKQg0wTnB4TAREQAREINwE5yHCXr7QTAREQARFIk4AcZJrg9JgIiIAIiEC4CchBhrt8pZ0IiIAIiECaBOQg0wSnx0RABERABMJNQA4y3OUr7URABERABNIkIAeZJjg9JgIiIAIiEG4C1g7y1ltv/fnnn/fv39+pU6cTTjjhhRdeCDcFaScCPicgk/R5AUm8UBKwdpDvv/9+2bJl33777SpVqnzzzTcPPPBAKJWXUiIQFAIyyaCUlOQMEwFrB0nfESXffffdfv36lStXLkwKSxcRCCIBmWQQS00yB51AcUsFevbsWbdu3aOOOuqJJ57YunVryZIlLaMpUAREwBsCMklvOCsXEYgnYNGDPHTo0Lnnnvvxxx8vXrz4iCOOOProo6dOnSpqIiACuSIgk8wVeeUbcQIWDrJo0aI333zzcccdV6xYMeiUKlXqxBNPjDgmqS8COSQgk8whfGUdZQIWDhIcXbp0ef311w8fPhxlNNJdBPxDQCbpn7KQJNEhYD0H+fDDD+/evbt48eLMPuImixQpwlsf0YEiTUXAbwRkkn4rEckTBQLWDnLXrl1RUF46ikBQCMgkg1JSkjNMBKwdJBru2LFj1apVe/fuNdq2b98+TGpLFxEIHAGZZOCKTAIHnYC1g3z22WfHjh27YcOGpk2bLly4sG3btrNnzw66qpJfBIJLQCYZ3LKT5MElYL1IB+/4ySefnHzyyXPmzFm6dGn58uWDq6EkF4EQEJBJhqAQpULgCFg7SNbmmM0B9u3bx44BK1euDJxiElgEwkRAJhmm0pQuQSFgPcTKFqw7d+4877zzOnfuzAuRlStXDoo+klMEQklAJhnKYpVSPidg7SDffPNN5L7nnns6duz4008/devWzedqSDwRCDcBmWS4y1fa+ZOA9RDrZZddZsQ944wz2ARy4MCB/pReUolARAjIJCNS0FLTVwSsHeTy5ctjUh48ePDTTz/1ldASRgSiRkAmGbUSl75+IJDoIEeOHFmmTJnPP/+c8yD5cF2hQoVevXr5QVbJIAIRJCCTjGChS2WfEEh0kLfffjt7dgwfPpy95fhwvW3bNkzUJ+JKDBGIGgGZZNRKXPr6h0CigzSS9ejRg71YuX7hhReGDRu2bt06/0gsSUQgggRkkhEsdKmccwLWDvLaa6/lGMjPPvts9OjRbBdw+eWX51xQCSACUSbg0CTZG7J169ZNmjRp0KDB3XffDbHt27fzslbt2rX5ZrM6w5AxoVq1atWpU2fGjBlRpirdRcCegLWD5BwPTvDgnOSh+R+bjZJlkPZ8dVcEXCHg0CRLlCjBrpA0bZctWzZ9+nT2iczLy+vUqRP7KvPNNcKsWLFi0qRJrPohwpAhQ1iF54qESkQEwkfA2kGyNoc2JuOr3bt3x372799fkOYyyILIKFwEXCTg0CRp15YuXZp8sVk+ppk7YMAAQvieMmUKFzR8+/bti+XWqFGDfuSiRYtclFNJiUCYCFg7yFdeeQX7mTBhwoknnrhx40bW7BSkswyyIDIKFwEXCTg3SVq0nDHA4nPGVNu0abN58+ZKlSohCd9btmzhAouuWrWqkY0NevjXRTmVlAiEiYD1Tjr4RdbmGD2rVatmPweJQbZo0eLbb7+97rrrCjLIU089VQYZpnojXTwm4NwkixUrxvgqW0Wef/75X375ZbKcHIEeH0gbNyHO+PwPgVu3bk1+XCEiEB0C1j3IlPQ3BsnZWIzVpGGQmGLL/I+sMSXsiiwCNgSOPfbYDh06MMtYsWLFTZs2EZNvupVc0Gtcv369eRazTd5pefDgwYvzPzrGx4awbkWBgAsO0mBK2yBljVGoZ9LRGwK0Muk7kteePXs++OADjuJhq8iJEycSwrfZ8YMQFulwUM+aNWtYvMOqV29kUy4iEDgCiQ6SpW7ocNtttznURAbpEJSiiUB6BFIySbqJHDDQuHHjVq1aMQfJ25MjRoyYOXMmr3nwzTUy8AZInz596tev37Vr13HjxjEClJ5gekoEQk8gcQ4SA5s3b960adNY5xY/V9G8eXNLFsRndRzTkIcOHcLqMMi2bdtywQIfJi8nT54cb5AsVZdBWmJUoAgURCAlk8Q1csJ5fFLHH3/8rFmzEhK/I/9TUI4KFwERMAQSHeR9993Hy1LMTMQW6RCPaXxerrJEJoO0xKJAEXCLQKom6Va+SkcERCDRQV6Y//nLX/5y1113iY4IiEDOCcgkc14EEiCyBBIdpAGBd2SUdf78+fzLWjgGTiMLSIqLgB8IyCT9UAqSIWoEEhfpGP05QGDs2LFM4/Phgn+jxkX6ioCvCMgkfVUcEiYiBKx7kO+88w7vGhct+l/3yRqcZs2a6cSriFQIqelPAjJJf5aLpAo3AeseJDqbt6m4+Omnn8KNQNqJQCAIyCQDUUwSMkwErHuQjOfQa+SFKt70YCZS3ccwFbl0CSIBmWQQS00yB52AtYPs168fa3M++eQTHOSoUaPYBzLoekp+EQg0AZlkoItPwgeUgLWDRBn2/mdLqoBqJbFFIHwEZJLhK1Np5HMCBc5B+lxuiScCIiACIiACWSUgB5lVvEpcBERABEQgqAQsHCS7qjZs2DCoCkluEQgdAZlk6IpUCgWDgIWD5PXHJk2afP/998HQQFKKQNgJyCTDXsLSz6cErBfpcIAAZ+JwUFypUqWM4Ow851MNJJYIRICATDIChSwVfUfA2kHefffdvpNUAolAhAnIJCNc+FI9ZwSsHeQZZ5yxbt06Ths/66yzfv31V457zJmAylgEROB3v5NJqhaIgPcELOYgEeKZZ57hkJ2rr76a640bN5533nneS6YcRUAEYgRkkqoMIuA9AWsHOW7cuAULFpQtWxaBateuvWXLFu8lU44iIAIxAjJJVQYR8J6AtYMsUaLEkUceaaQ5cOBAkSJFvJdMOYqACMQIyCRVGUTAewLWDpIJj/vvv3/Pnj0zZ8686KKLzj33XO8lU44iIAIxAjJJVQYR8J6AtYPMy8srX758o0aNnn766XPOOeevf/2r95IpRxEQgRgBmaQqgwh4T8B6FSsvJnNOcps2bRhcrVOnjoZYvS8Y5SgC8QRkkqoPIuA9AWsHyfHl11xzzSmnnMJxV2vWrKEf2a1bN++FU44iIAKGgExSNUEEvCdg7SBvvvnmOXPm1KpVC4FWr17dvXt3OUjvy0Y5ikCMgExSlUEEvCdgPQdZoUIF4x0RqGbNmvzrvWTKUQREIEZAJqnKIALeE0jsQb7xxhsIwUasrM3p06cPs4+TJ09u1aqV95IpRxEQAQjIJFUNRCBXBBId5FtvvWVEqVix4rx587hmOeuOHTtyJZ/yFYGIE5BJRrwCSP0cEkh0kM8991wOpVHWIiACCQRkkqoSIpArAokO0sjBytXHHnts7dq1bKNjQnTcVa5KSPmKAARkkqoGIuA9AWsHye7kgwYNYgMd3r7yXiblKAIikEBAJqkqIQLeE7B2kCVLlrzxxhu9l0Y5ioAIWBKQSVpiUaAIZJWAtYMcOnTovffe26VLF7ZINtk3b948q3IocREQARsCMkkbOLolAlkiYO0gv/jii+eff3727NlmiJWXPbjOkgRKVgREoFACMslCESmCCLhOwNpBvvnmm999913sxCvXc1WCIiACKRGQSaaES5FFwBUC1mtwmjRpsnPnTlcyUCIiIAKZE5BJZs5QKYhAqgSse5CbN2+uW7cuG+jE5iD1mkeqZBVfBFwkIJN0EaaSEgGHBKwdJCt0HD6vaCIgAh4QkEl6AFlZiEACAWsHyfHlIiUCIuAfAg5Ncv369Zdffvl//vMfltcNHjyYta/bt2+/+OKL2fSjevXqr7766nHHHYdSI0eOnDBhQrFixR599NGzzz7bP2pKEhHwFQHrOcgyZcqUzf/w9hVWxEVBQmOQHTt2rFevHvubjx07lmgYZOfOnWvXrs13bBNXDJLjQTh7ecaMGQUlpXAREIGCCDg0yeLFiz/00ENfffXVwoULx40bt2LFiry8vE6dOq1atYpvrkmfwEmTJi1fvnz69OlDhgw5ePBgQZkqXAQiTsDaQe7atevn/M/evXtff/3166+/viBMMsiCyChcBFwk4NAkK1WqZF5ZxqHSbN24cePUqVMHDBiAJHxPmTKFC0L69u3L8oIaNWrQbF20aJGLciopEQgTAWsHGa8he1zZvAQpgwxTbZAugSBgb5JGBcZUly5d2qZNG1b3YKQE8r1lyxYu8JpVq1Y10apUqcK/gdBaQoqA9wSs5yDNEXRIc+jQocWLF7NRQKGS2RvkqaeeKoMslKEiiEBBBFIyyV9++eWCCy4YM2aM5eTI4cOH43NJtu7x+R/ibN26tSB5FC4CUSBg7SBjR9AxgsrcPmMy9iwyMUhZoz1b3RUBCDg3yf379+Md+/fv37t3bx7kYNdNmzbRfeS7QoUKhNBrZOmAobphw4bKlSsnEGZ1Dx8CW7ZsKfgiEGUC1g4ypSPoMjRIWWOU6590d0jAoUnSO+QcHmYfhw0bZlLu2bPnxIkTR4wYwXevXr0IJOSSSy4hwg8//MDindatWzuUQdFEIGoEEh3kfffdl4yAQZi77rorOZwQGaQlFgWKgFsEUjLJBQsWsItyo0aNmjZtigD3338/rrFPnz681FGtWrXJkycTyIJzQurXr8/4ECtdWabulqhKRwRCRiDRQZYqVSpew927d2Na27ZtK8hByiBDViGkjt8IpGSS7dq1S5hiRJ1Zs2YlKHVH/sdvmkoeEfAbgUQHefPNNxsRWVbOe40M7LAiPBaYLL0MMpmJQkTARQKpmqSLWSspEYg4AYvXPHjT/84772zcuPGBAweWLFkyatQoM7cfcVJSXwRyRUAmmSvyyjfiBBJ7kMOHD2dBOQtnOH+udOnSEacj9UUg5wRkkjkvAgkQWQKJPUj2qWJt21//+lcWf5vd5sweV5EFJMVFILcEZJK55a/co0wgsQfJzgBRxiHdRZ9qbBYAACAASURBVMBvBGSSfisRyRMdAok9yOhoLk1FQAREQAREwIaAHKQNHN0SAREQARGILgE5yOiWvTQXAREQARGwISAHaQNHt0RABERABKJLQA4yumUvzUVABERABGwIyEHawNEtERABERCB6BKQg4xu2UtzERABERABGwJykDZwdEsEREAERCC6BOQgo1v20lwEREAERMCGgBykDRzdEgEREAERiC4BOcjolr00FwEREAERsCEgB2kDR7dEQAREQASiS0AOMrplL81FQAREQARsCMhB2sDRLREQAREQgegSkIOMbtlLcxEQAREQARsCcpA2cHRLBERABEQgugTkIKNb9tJcBERABETAhoAcpA0c3RIBERABEYguATnI6Ja9NBcBERABEbAhIAdpA0e3REAEREAEoktADjK6ZS/NRUAEREAEbAjIQdrA0S0REAEREIHoEpCDjG7ZS3MREAEREAEbAnKQNnB0SwREQAREILoE5CCjW/bSXAREQAREwIaAHKQNHN0SAREQARGILgE5yOiWvTQXAREQARGwISAHaQNHt0RABERABKJLQA4yumUvzUVABERABGwIyEHawNEtERABERCB6BKQg4xu2UtzERABERABGwJykDZwdEsEREAERCC6BOQgo1v20jyUBAYOHFihQoWGDRsa7bZv3965c+fatWvzvWPHDhM4cuTIWrVq1alTZ8aMGaGEIKVEwBUCmTpIWaMrxaBERMAtAldcccX06dNjqeXl5XXq1GnVqlV8c034ihUrJk2atHz5cqINGTLk4MGDbmWtdEQgZAQydZCyxpBVCKkTdALt27cvV65cTIupU6cOGDCAf/meMmUKF4T07du3RIkSNWrUoB+5aNGioKss+UUgSwQydZCyxiwVjJIVAVcIbN68uVKlSiTF95YtW7jYuHFj1apVTeJVqlThX1cyUiIiED4CmTrIBCKyxvBVEWkUMgKHDx+O16hIkSIJCo4fP75l/mfr1q0h013qiEBKBFx2kMl5yxqTmShEBDwjULFixU2bNpEd3yze4YJe4/r1640AGzZsqFy5coIwgwcPXpz/KV++vGdyKiMR8CEBlx2krNGHZSyRokygZ8+eEydOhADfvXr14oIQFuns27dvzZo1LN5p3bp1lPlIdxGwIeCyg5Q12rDWLRHwgEC/fv3atm27cuVKeooTJkwYMWLEzJkzec2Db64RoEGDBn369Klfv37Xrl3HjRtXrFgxD6RSFiIQRALFMxQaa5w7d+6PP/6INd57771YILaHWVarVm3y5Mnx1li8eHFZY4a09bgIFErg5ZdfTogza9ashJA78j+FJqUIIhBxApk6SFljxCuQ1BcBERCBsBJweYg1rJiklwiIgAiIQNQIyEFGrcSlrwiIgAiIgCMCcpCOMCmSCIiACIhA1AjIQUatxKWvCIiACIiAIwJykI4wKZIIiIAIiEDUCMhBRq3Epa8IiIAIiIAjAnKQjjApkgiIgAiIQNQIyEFGrcSlrwiIgAiIgCMCcpCOMCmSCIiACIhA1AjIQUatxKWvCIiACIiAIwJykI4wKZIIiIAIiEDUCMhBRq3Epa8IiIAIiIAjAnKQjjApkgiIgAiIQNQIyEFGrcSlrwiIgAiIgCMCcpCOMCmSCIiACIhA1AjIQUatxKWvCIiACIiAIwJykI4wKZIIiIAIiEDUCMhBRq3Epa8IiIAIiIAjAnKQjjApkgiIgAiIQNQIyEFGrcSlrwiIgAiIgCMCcpCOMCmSCIiACIhA1AjIQUatxKWvCIiACIiAIwJykI4wKZIIiIAIiEDUCMhBRq3Epa8IiIAIiIAjAnKQjjApkgiIgAiIQNQIyEFGrcSlrwiIgAiIgCMCcpCOMCmSCIiACIhA1AjIQUatxKWvCIiACIiAIwJykI4wKZIIiIAIiEDUCMhBRq3Epa8IiIAIiIAjAnKQjjApkgiIgAiIQNQIyEFGrcSlrwiIgAiIgCMCcpCOMCmSCIiACIhA1AjIQUatxKWvCIiACIiAIwJykI4wKZIIiIAIiEDUCMhBRq3Epa8IiIAIiIAjAt45yOnTp9epU6dWrVp5eXmORFMkERCBrBGQPWYNrRIODwGPHOTBgwevu+669957b8WKFS+//DLf4UEoTUQgaARkj0ErMcmbGwIeOchFixbRd6xZs+aRRx7Zt2/fqVOn5kZd5SoCIvC738keVQtEwAkBjxzkxo0bq1atagSqUqUK/zoRTnFEQASyQUD2mA2qSjN8BIp7o9Lhw4fjMypSpEjs3/H5H/5duXJly5YtC5Lnxx9/PIFPQbezEG5yzELCdkl6n6n3OaK/95l6n6NRs+UHd9uVd/69devWFRrH3Qg29khGDk2SmMYec8I2bSABkjZAolIcQZG2Zcv/mqSNtPH26JGDpNe4fv16U6c3bNhQuXLlWP0enP8ptLrjOxcvXlxoNBcjeJ8jwnufqfc5Sk0Xa2l6SdnYIwk6NMlY1jmpQukpnpO6FwVRgwXWubQeDbG2atVq1apVa9as+e233yZNmtSzZ8+0K40eFAERyJCA7DFDgHo8IgQ86kEWL1788ccfP/vss1k+N3DgwAYNGkSEr9QUAR8SkD36sFAkkg8JFLvnnnu8Eat27do33HDD0KFD27dvn16OLVq0SO/BtJ/yPkdE9T5T73OUmmnXSbcezNwe4yXJSRVKG0WApA2QqDkx6rTrgENpiyRM12eSn54VAREQAREQgdAQ8GgOMjS8pIgIiIAIiEBECOTeQTIlWaFChYYNGxriw4cPr1u3buPGjc8///ydO3cSuH///gEDBjRq1KhevXojR45MKBiGiE866aSm+Z93333XSbEVmiMria688kpybNKkydy5cxPS3L59e+fOnRmh4nvHjh1OciROhplmruZdd90FVSh16dLlhx9+MGLDkw0c2AJwxowZmauZoGOqOaahYzLY5Ey3bdvWsWPH0qVLX3/99cmFlUZpFqqmfY7pqZksuZchCSp/9tlnbdu2xUDOPffcn3/+OUGSnG9il5K01atXRxHswuYds+yhZm0/lZNfNpZljB07lozsK2Ru2aYqbQ7ZJos6efJkIBctWtTybYgCwTLEmtvPvHnzPv30U0Q3YvBLjUfk+tb8DxcvvvjixRdfzMXu3btPPvlklsLGC3z33Xc/8MADKalQaI6sJ7riiitIc/Pmzc2bN2dhUXz6uHD8CiF8Gwmd5J5hppmr+dNPPxk5scOrr76a6+XLl+My9+7d+91337HJ0YEDBzJUM0HHVHNMQ0cELjTTX3755cMPP3zyySfZ7DC5pNIozQxzTE/NZMm9DElQGV9CwxEBJkyYcOedd8ZLQi2iLq1evXrfvn3ULuqYl3KavJxLS3x+UrZu3eq9kCZHmqr8+nFNO4M2N7hsKmTO2aYkbW7ZJovK/qZff/31GWec8cknnyQUtw3Y3PcgWbNTrly5WBuN/g1L7Pj31FNP5Y1JLthVANeIDnv27GGnurJly2bYoCs0R1B26tSJXOjaHnvssQktDvbJo0fLXb6nTJniUJgMM3WYS3y0hBxj3IBpNmpAEbb9K1GiRI0aNehHsv1Y/ONpqJlhjmnoyCOFZlqqVKl27dqVLFnSMv1sqGmfo6UYPg9MgMyeHmapHYMor7/+erzwftjEzrm0OcdeqVIlmuCIUaZMGfqR7HBkUyFzzjYlaXPLNllU8DJUZimVDdjcO0hLiQn8+9//3q1bNy4uvPBCfnFQuFq1arfccku8NzXP0uGjrcq4ivMBT8tMYzkysko1xSXTW6V9F9viwDxFtxJhuOZ7y5Ytlkk5D3SYKQlmruYdd9zBhn/0yO+77z4StN9vzBU1U8rRFR1JJCFT+7LIhpr2ObqlZqG5ZC8CEyLTpk0jfYatEqzDvlJlTySblG2k5SkaizTKWSxq9vOySSert9auXbt06dI2bdrYVEj/sHUirU/YxkS1KT4bsD51kH/729/oR/bv3x+tcO/FihWjy4y7euihhxgPjFf12muvZTxn2bJluKubb77ZhoL9rfgc8bVsNcI40k033XTaaaeZHq394+nddZ6pK2qSHT9nUMXXIjDjDPFix+//l546yU+llKMrOiJDQqbJUrkeklKObqnpuhbOE6RVN27cOJzKrl27GNSJf9CDSuVcThPTRloiLFiwYMmSJRw0hEbz589PNXFX4jMLcMEFF4wZM8Z+eMwnbB1K6we2DkW1AetHBzlx4sS3336bjo75yX7ppZe6du16xBFHMOD5hz/8IWHAs2LFirhPpl7/+Mc/JgwSOq+7CTniER955BGcLv1IFgoxNxCfFDlu2rSJEL4RyXkuCTFTzTRzNY0Al1xyiRkWs99vzC01ydR5jm7pGJ+pfQFlQ81Cc3RRTfu8snSXZXTvv/8+gyv9+vU75ZRT4nOxr1RZksc+WRtpedBse4khsyow7R8QewHs77LkAu9Is7V3797EtKmQfmDrXNqcs00Q1aYUbMD6zkGymmjUqFEM4Bx99NFGJUZWZ8+ejZNn8mzhwoVU93hVja8i5M0334wthbVhkXwrOcdff/2VvIg5c+ZMnGX9+vXjn2KfPHwbIXz36tUrOUEnIalmmrmabPVnBIOtYYgibPvHYgq65txt3bq1u2qmmmPmOiJ/cqb2xZF5aaaaoytq2iuV7btmZuHQoUN//etfr7nmmvjsfLiJnY20mDmdYOTnApef3g9IJrT5WRs0aBDTY8OGDTPp2FTInLNNSdrcsk0W1aaY7MCSUG4/rBM58cQT8UO8rfHss8/SIMWfMwvIxyy2pAYzDYmXohqNHj3aSEutMouRLr30Uqq1WXHOMKwTXQrNEYfx+9//Hi/CUh2GsBNyZBv4M888k1UtfLOg30mOxMkw08zVpH3KUmFA9ejRg9VPRmx+4FhziLK8IZO5mgk6pppjGjomg7XMlJWKxx13HDPZ1DGzrjJWf9IoTSdq2uSYnpoOq1mWoiWozGAgwyp8brvtNtwkmTKLw4oBk/s777zDLeoVtStL8tgn61xaZmdYvsCHn5ecSMv6an64zRtl/OKBLrlC+odtStLmlm2yqG+88Qbmz4wAowXMOjustNpJx6ZhoVsiIAIiIALRJeC7IdboFoU0FwEREAER8BMBOUg/lYZkEQEREAER8A0BOUjfFIUEEQEREAER8BMBOUg/lYZkEQEREAER8A0BOUjfFIUEEQEREAER8BMBOUiPSoN3wzkxgDdSOADBnFLCSym8vlJQ9pw+kXCLfV/ZJLag+NkOP+ecc4zYKWXEnu9s9PrUU08lPJWsnU2ybN/Mi0APPvigTRzdEgF3CRiD5c0o3r54+OGHeZvFYfrU+ddee43IV111lXODZfP3jz76yGSBvfzzn/90mF0sGikcc8wx2GnCgzF5nCTI2xG88eL9+6BOZMtJnP9uC66PBwSOOuootuYhI7Y4Z1Mr9gtlCw9jSA5zx0HyCmPCrgUOn808msOjxJIz4qwVm3ZAcvzkEFLgFcbkcIWIQPYIxAyWfQbYCoqjae69996UsuOtbufxcW+0GtnYkkcS9l5wnsjpp5/OHmTO4yfHJAUsnd+Z5FvRDFEP0uty5yA93vwlV7YgMC01Nu7p06cPbytzqhe7Fcf20sOJ0nrlVBP2L6Z1yQ449KXohvIGbkzot956i0eaNWt21llnEY1wzvoxp2MSaHYJYXcF8ybyiBEjiMDjbN3HRpoYA+e/EEIb88Ybb8Q4ebnb+Gx2e+FIBNPlNa/ccrQbrzATmdY0YvPhbXGjBRs4sM8fbW1ev+XElYKAsv0CurNpBUc2xuLg/AhBdw6BMoF/+ctf2KKBYyLYxky9xoJgKtwzArxXzj7mbF/Mq+WcfIcNmhr79NNPGxkS7CsmWIcOHYwt4/nibZkICWbLTwG9Rra3xOIwN04MNTWfJjXmj3WwDZ45iYE02ZyBTa/Y3MMYpiUHROX0UxrT3bt3jx2owNaAHPaE4Z999tlmOyf2WiFxrBKl1Gu0JPnfHav18YAAfSBy4YQQulPsjMw1DsOcgomTGDx4MCFffPEFAztmhyBKC49IIHUXn8EFXU8OT0gQleNVzVYmzzzzDLtVcZfW37/+9S8u8I7sRkh7EANg2ydCzL4/bAD0zTff8C/79nFYq0kZqTB+dplhJyNCsE+zsQgCc1IdIebYPAweQ2ILYBLH/NjlGS2QmYMIiHPRRRc9//zzXMQ+8TIztszmfNzit8bQ4OxPPCvykzWWjGtHd9oEtBjIlL2KYid9BvEMxXgOug4cAVNFY2Jz7N1//vMfnKIxRk5RxdNwcEKyfcXqfOzowWRbTjbb+Boeu6Zda87dpE05dOhQhCFNY+bsucM+X/FU58yZgxGZEDZbpsWM8dIWZ9yV3w0Oged3AGdJBDaY5EB4Lvj9Ya92LnC6sRN5Y79LgSuybAisIVbrdoProXStaB7SVMSu6B7Fp48/o/YTgu+hQWdusSWSGeggPlvCFiQPm8bR76Q9iAEw20c09nPHhMzex2za98EHH2AMZmNbTgrDt9EZxZOZBNmI1Vycd955bPiOzzPdUNrIHGmCfyUcsWO5IyqNWTPgyaZutGHZOpJ8TRxERcGCRMUUzSbpl112GQZJNHa/5ENPl2sEY1NT/C7b2zK6RQgOtaCkFC4CHhPgx9fU2M8//9yMsjDoSo1NsC9LqZJtOdlskx8kfWb98YjcwunGbNbsaW5vaxxLwgAMLVfmcWgQkwJHeH755Zfml4f2KGcfkTjmZgZ1GUPOcGw2Wf5whGiI1aNyNFMa69atw5MxBxmfq7G9BDk4vcQcZkItpyVYkJQ33HADYyl0PWnb0qolGuOoTH7gjxmcYQSVxOPPsaK7RluYoRvz+eqrr0zKHJtsLowwjK9iY2xdiDOLXy9gKWrsWXtRST/hRC1Su/32240k3377LfujWqZfkO4KFwFvCNBTpG4z1kr9fOyxx0yNpadltvQs9Jy4ZFtONlvnihhzS8PW6CMayfm5oGEqW3PCXA7SCSXX4jDc8eijjzKASecslijn3b/66qv8y5o36m5BmXHmuJlTjI9AMxM3Rog5YIQLphgZmaGLxnmWOEhsmPPwGLTkFgM7HDhHh48hF/7FQj777DPL7HDk/Bww/onTYhw1FgfHyVohc9oJx6cwi2n5uGUgXVvGdrjFQWYmAnMhyEbfkX8ZC2L8BxRMz+DpCWQQyTIdBYqAlwS2bt3KqhmaoThCauyTTz5pjJd5CmYuEuzLoWDJZpts3fxWsMO+mWhk5sJ0JR2mj51ia/QUGVti6JWn6tSpgyIff/wx18jPZAqJkynzLIQYw3SYeKSiyUF6XdyMKDLNFl8jhwwZQt1lcJVzvvjGMCxl4owC5uR4PH6RDvP5jL3gqE444QTzFGtnGKolC/qsHLDAehxGQXGWjIKamX/804QJE4hAi5IDLy3zYuaD+OTFoKgZ/jXRmjdvzooe1giwMohV7GZ01DKF5MCxY8fSdWbwll8Hc5cfF8Z2mBrBozMJivvnLtIiG+NIyFwQiuTEFSIC7hIwcyLYCJN5VFSziIw6zzQEVoCJcdYQQzvJ9uVEjGSzZUKBFqdZpBNLgVYvSxD4TaDn9+c//9lJyiYO8yAcqIJZcTq38awM8zIyTLsZ4yIX80oJvwOsfsAAaSvL1izx6jQPSyyeBtLQo01XsmRJPB8T77RME05p91QaVzPDmzKTmtJrHvQdWfVHJ5VWMKsH+TFCIn5QCLzllltclU6JiUB4CNCopQWc0lSisTUQ5OXl0dekCcs1ywiwWSYsw4MmA03Ug8wAnkuP4gwYWqRlR7uPAZzQeEfw0CxlAV7yRgE25GjS0sLFL3LMuvGONKJfeOEFvQppA023RIDfDbxa8kYBNmSYxcDW6A0zkHvnnXcSkwv6srHhKJtnI3JLPciIFLTUFAEREAERSI2AepCp8VJsERABERCBiBCQg4xIQUtNERABERCB1AjIQabGS7FFQAREQAQiQkAOMiIFLTVFQAREQARSIyAHmRovxRYBERABEYgIATnIiBS01BQBERABEUiNgBxkarwUWwREQAREICIE5CAjUtBSUwREQAREIDUCcpCp8VJsERABERCBiBCQg4xIQUtNERABERCB1AjIQabGS7FFQAREQAQiQkAOMiIFLTVFQAREQARSIyAHmRovxRYBERABEYgIATnIiBS01BQBERABEUiNgBxkarwUWwREQAREICIE5CAjUtBSUwREQAREIDUCcpCp8VJsERABERCBiBCQg4xIQUtNERABERCB1AgUTy16NmOfcMIJ1atXz2YOSlsE/Ehg3bp1W7dudUuynTt3XnXVVV9++WWRIkX+/ve/16lT5+KLL167di3G9eqrrx533HFkNHLkyAkTJhQrVuzRRx89++yzC8paJlkQGYWHmEC8PRY5fPiwT1Rt2bLl4sWLfSKMxBABzwi4W/MHDBhw+umn4yN/++23X3/99f777y9XrtyIESPy8vJ27NgxatSoFStW9OvXb9GiRT/88MNZZ531zTff4CktlXVXMMssFCgCfiMQX+01xOq30pE8IpA+gZ9//nn+/PmDBg0iiSOPPPLYY4+dOnUqLpN/+Z4yZQoXhPTt27dEiRI1atSoVasWnjL9/PSkCISagBxkqItXykWMwHfffVe+fPkrr7yyWbNmdCJ37969efPmSpUqgYHvLVu2cLFx48aqVasaMFWqVOHfiEGSuiLglIAcpFNSiicC/idw4MCBJUuWXHvttUuXLi1VqhTDqskyJ8yqMFWZEGf8+PGMMvFxcWY0WQyFiID/CfhokY7/YUlCDwjs379/w4YNe/fu9SCvnGRRsmRJ+m1HHHFENnInZT5t2rQh8QsvvBAHWbFixU2bNtF95LtChQqEE2H9+vUmd1BXrlw5QZLB+R8C8ZHZEDKyaYa+bgexZO3tUQ4yiGUaZpn5yS5TpgxLLpN7NiFQm97btm3b0JH5v2yoc+KJJzJ8unLlShavzpo1q37+Z+LEiSzS4btXr15k2rNnz0suuWTYsGEs0lm1alXr1q2zIYnSTCYQ7rqdrK//Qwq1RzlI/xditCSk7xhW70hB4vWPP/74rA5dPvbYY/3792cJa82aNZ977rlDhw716dOHlzqqVas2efJkZGjQoAEhOM7ixYuPGzeuoCWs0ap2nmgb7rrtCUKXMynUHuUgXSau5DInEMq+YwxLtrVr2rRpwutSdCUTCuWO/E/mJaUUUiWQ7dJPVR7Fty8RFxbpPPLII7RJGzZsyMtVNJG2b9/euXPn2rVr8817V6YAeDGZBeUM+8yYMUNFIgK5JUCfCS9Cjb3ooot4UzBVYXjpnmd5au7cuT169Ej1ccUXgewRMHWbH+QmTZo8/PDDjB+QFw2mG2+8MTlTavJLL72UHO5WCNtQ1KtXj/EMVxLs0KEDHgTL5cP8eqFpLlu27N133zXR7rnnngcffLDQR5IjZNqDZI04FHj1+KijjmLcZtKkSVx36tTJvJjMGgHzYjLhy5cvL/TF5GT5FBJlAtVHvJO2+mvzuhf0LHUV4+EupvvUU08xG1dQTBaFMg5Z0F2Fi0DaBLJdt3mlh5nmn3766d577zVrkpNFNQ6SaMm3XAl54okn3nvvvfjp9gwN6sUXX3S4cIyMsHFaBuecc04murjQg0SUPXv28E1jnBVx/nwxOZPqmAlfPetnAuw48+2337711lss++TFQbaV4a1BBKa9yTLOLl26XH755fyIEK15/uejjz6yVId37U877TRS4JsFMsSh8T5w4EAuvvjiC7qbafRTLTMKXKDsLldFxoplXtd5/PHHWYcSG+qYN2+e6YFRV3ft2kU35sMPPySEUcDkes5TdNroq9WtW5empHk76JNPPqGS0z1lbRcpHDx4cPjw4a1atWrcuPHTTz8dr+w111zDW7msCCPxeINiIzd6UMTn+/vvv+eRK664gheTOnbsyKw5EmI49DsJLBSdveX++c9/fuWVV9COb5Ki54Y6ZEGPrtCUYxEybR2fdNJJt9xyC/P/tMr5QeFj+WLyqaeearLUi8nOy0Yxs0qAJh3N265du7Zr127hwoVMRTz77LOjR49+6KGHyPfTTz/917/+Ra3Gt82cOZO14Cz4ZBLBcjdEfkHYv4a+5gcffPCnP/3p9ddfv+mmm7DGN998829/+xs/HEcffXRWdVHiIpBMAGfAEKvZHcLcZZiRZVl/+MMffvnlF6o0I3yEvP3229y1rOe8TcvIH90eHlmwYAFOkX198Td4RPZswjpY/HXMMcfgNfft20ccfv9j/UXGZqZPnz5nzhx29MVBxgzq3HPPpd3Jvk5sFMzAr9ndicm42bNnT5s2jbtkhCWSBV1A3Fu8XvhpMiWE+bsHHnjA3nL/8Y9/YK00EYiPAF9//TXC4NQZp8UfO3zPKlMHiWJ0GdesWcOmVszovPDCC8nlZP9iMs0cPjyV1aV9yVIpJLIEGPAwhkfXkF3Z6PNh9rwmyMrPmHnT8jWmyLtr119/PbbK7A7bllpCYyALg8eD4mWJT5yiRYtinzSTr776an44LJ9SoAhkm0DCby9VkQkF3Ezv3r3pq8TnblnP8YgmGvZCFxNfyAu1uC4eLFu2LN/vv//+559//tprr3GNFWACBb2/FDOojz/++I033iD+ZZddduuttxoZ8IvYTqNGjXhtl28CmUYlxwQHmTDEymszNpabwLZ79+5sr8iHvjW9uAT1EyLH/s10iJUmM0TY3QqHDHTGoMyLyWTg8MVkxrLw83xIpCApFS4CLhIwc5D4PN6IYMPSG264ARfIWChdvdgGBWxDY3JkgIgq/dlnn1FF8aCWYtx1110MEHGABmM+sRT4sShdujTz7paPKFAEsk2AEU5adWZ3CJMXY6p0zmggMqRHjypeAMt6jjsxcUiHERfcbcKaT0IwIkyJD90kepAFKRUzXZBjVwAAIABJREFUqPgIsdRMRjQrYzlyTY4FpWbC7S034dkEXexTjt3N1EEyuMrwFN1zSLGanLFjWgq8kkwG8S8ms0iHPjgE9WKyw4JRNM8I0PJlpsDU2ORMuUurGXN9/vnnmXFJjkBILAV6jSYCIUOHDmXclW0BTPva8kEFikCWCDAgxywgLb94l7Z69Wr6Z7fddhtLXXCQ7MjBkGOsxhZaz5lKoMHHgCqP8CAOjLPSnnzySTNqwvgKe/8Wqg5TmLgDotEdZIy00Pg2EewtN147m0Tsb2U6xMrqBmZxWb/ABAwTv3QHGd3Wi8n20HXXVwSYn2B2AB9Js5o2XIJsQ4YMueCCC3jFnj6iZSuY+IwUMcTKwpwzzzzTPP7//t//48Hf//73TNLwYPv27eMb8r5SX8KEiYCZPsBj8YPMGGbCCu0xY8YwD0d3kG0iunXrRrOPaKy4YUWMk3rOcAsTkPTbyIVhGIYP2RCfgVB+/+kgMQRoJhTtebJGhmU4zCASn40s7CMn3I3NQTKvSe72lovdMcnKIO3tt9+eUi7xkaNyHiSr6WzW/aeNTw+6TuCrr75iHML1ZH2VYIKONOct1/7kXObMBZPdxRdiFOp2zittGgLY2GOmQ6xpSOPBI1pc7gFkZSECIiAC4SYQTgcZ7jKTdiIgAiIgAh4QkIP0ALKyEAEREAERCB4BOcjglVnoJU54eStk+oZbu5AVluvqqPRdR5phgvYlIgeZIV497jIBNvjg1Qj7Wutylh4mh15oh44e5qms/EIg3HXbL5RTkaNQe8z0NY9UhPEirpbneEE5m3mwwwUbZIR4WyVzgnk2ESptnxIIfd32KXdbseztMWwOMhmFXGYyEz+HsCVTQbtV+VlsySYChRJQ3S4Ukd8ihGqIVb7Qb9VL8oiACIhAcAmEykEGtxgkuQiIgAiIgN8IyEH6rUQkjwiIgAiIgC8IhNZBMtyqEVdfVDEJIQIiIALBJBBaB2mKQz4ymNVSUouACIhA7gmE3EHGA1afMvfVTRKIgAiIQHAIRMhBBqdQJKkIiIAIiEDuCYTHQWo0Nfe1SRKIgAiIQIgIhMdBhqhQpIoIiIAIiEDuCchB5r4MJIEIiIAIiIAPCchB+rBQJJIIiIAIiEDuCYRhL1bNPua+HkkCERABEQgdAfUgQ1ekUkgEREAERMANApk6yJUrVzb9/z9ly5YdM2bM9u3bO3fuXLt2bb537NhhhBw5cmStWrXq1KkzY8YMN8RWGiIgAiIgAiKQXQKZOkh83rL8z6effnr00Ueff/75eXl5nTp1WrVqFd9cI/6KFSsmTZq0fPny6dOnDxky5ODBg9nVSamLgAiIgAiIQMYEMnWQMQFmzZp1yimnnHzyyVOnTh0wYADhfE+ZMoULQvr27VuiRAnO+aMfuWjRoozFVgIiIALWBKpXr96oUSOGdVq2bEkMjehYY1KoCDgg4JqDpI/Yr18/cty8eXOlSpW44HvLli1cbNy4sWrVqkYYztTmXweCKYoIiECaBObMmcOgzuLFi3leIzppQtRjIvC737njIH/77bdp06ZddNFFlkgPHz4cH16kSJH4f8ePH09Tl8/WrVstH1egCIhA2gQ0opM2Oj0oAu44yPfee6958+YVK1YEKN+bNm3igu8KFSpwQa9x/fr1hvWGDRsqV64cz33w4ME0dfmUL19e5SECIpAhARqgXbp0adGiBU1PktKIToY89XiUCbjjIF9++WUzvgrKnj17Tpw4kQu+e/XqZUIYgN23b9+aNWtYvNO6desoE5fuIpBVAgsWLFiyZAlt1nHjxs2fPz85L/sRHeJrUCcZmkKiScAFB/nrr7/OnDmzd+/ehuCIESP4l9c8+OaawAYNGvTp06d+/fpdu3bFaIsVKxZN1tJaBDwgYEZoGLxhSTkL4lId0UFCDep4UEzKIhAEXHCQvN2xbdu2Y445xih8/PHHs6KVniLf5cqVM4F33HHH6tWreWmyW7dugeAiIUUgiAR27969a9cuJOfi/fffb9iwoUZ0gliOktknBMKw1VxKKNmXbm1e95QeUWQRCAoBZhzpOCLtgQMHLrnkEsZsWrVqxfjNhAkTqlWrNnnyZG7FRnSKFy+uEZ2glKzkzAmByDnInFBWpiLgDYGaNWt+9tln8XmZEZ2E3BnR4eONSMpFBIJLwIUh1sApTydS+5sHrtQksAiIgAh4TCCKDtJjxMpOBERABEQgiATkIINYapJZBERABEQg6wTkILOOWBmIgAiIgAgEkYAcZBBLTTKLgAiIgAhknYAcZNYRKwMREAEREIEgEpCDDGKpSWYREAEREIGsEwi8g9QLG1mvI8pABERABCJJIPAOMpKlJqVFQAREQASyTkAOMuuIlYEIiIAIiEAQCchBBrHUJLMIiIAIiEDWCchBZh2xMhABERABEQgigeg6SK3uCWJ9lcwiIAIi4BmB6DpIzxArIxEQAREQgSASkIMMYqlJZhEQAREQgawTkIPMOmJlIAIiIAIiEEQCcpBBLDXJLAIiIAIikHUCcpBZR6wMREAEREAEgkhADjKIpSaZRUAEREAEsk7ABQe5c+fOCy+8sG7duvXq1fv444+3b9/euXPn2rVr871jxw6jwciRI2vVqlWnTp0ZM2ZkXSdlIAIiIAIiIAIZE3DBQQ4dOrRr165ff/31Z599ho/My8vr1KnTqlWr+OYaCVesWDFp0qTly5dPnz59yJAhBw8ezFhsdxLQq5DucFQqIiACIhBGApk6yJ9//nn+/PmDBg0CzpFHHnnsscdOnTp1wIAB/Mv3lClTuCCkb9++JUqUqFGjBv3IRYsWhZGkdBIBERABEQgVgUwd5HfffVe+fPkrr7yyWbNmV1111e7duzdv3lypUiUg8b1lyxYuNm7cWLVqVYOtSpUq/BsqhFJGBERABEQgjAQydZAHDhxYsmTJtddeu3Tp0lKlSpkx1QRQhw8fjg8pUqRI/L/jx49vmf/ZunVrSoQZINUYaUrEFFkEREAERMA5gUwdJD1CPm3atCFLlurgLCtWrLhp0yb+5btChQpcEGH9+vVGpg0bNlSuXDlevsGDBy/O/9ATdS63YoqACIiACIhAVglk6iBPPPFEhk9XrlyJlLNmzapfv37Pnj0nTpzIv3z36tWLC0JYpLNv3741a9aweKd169ZZVUmJi4AIiIAIiEDmBIpnnsRjjz3Wv3//3377rWbNms8999yhQ4f69OkzYcKEatWqTZ48mfQbNGhACL6zePHi48aNK1asWOaZKgUREAEREAERyCoBFxxk06ZNGSCNl5KuZILQd+R/sqqJEhcBERABERABFwlkOsTqoihKSgREQAREQAT8Q0AO0j9lIUlEQAREQAR8REAO0keFIVFEQAREQAT8Q0AO0j9lIUlEQAREQAR8REAO0keFIVFEwBUCbHfMzlY9evQgNR0e4ApSJRJNAnKQ0Sx3aR1mAmPHjuXYAKNhgA4PCHORSLdgEpCDDGa5SWoRKIAAm1W98847bIxs7uvwgAI4KVgECicgB1k4I8UQgQARuOmmm0aPHl206P+Ztg4PCFDZSVS/EbB2kLfeeivnWO3fv58zHU844YQXXnjBb3JLHhGIFAGHJvn222+zAXKLFi1s4NgfHsCDaZ8fYJOpbolAEAlYO8j333+/bNmyGBv7jH/zzTcPPPBAEHWTzCIQGgIOTXLBggXTpk2rXr06J7DOnj370ksvTfXwAIjp/IDQVBspkiEBawdJ35F033333X79+pUrVy7DPPS4CIhAhgQcmuTIkSOZg1y7di3HA5x55pmM/ejwgAzJ6/EoE7DeixWjqlu37lFHHfXEE09wTGPJkiWjzEi6i0DOCaRtkiNGjNDhATkvPgkQUAJFEiYkUIPjOBYuXMgycUZZOXlj9+7du3bt4lirbGvImckJm57b5+jKaclr87rb56K7IpBtAoXW/KCYZDIojFQmloxFIX4mEG+PFkOsrH+7+eabjzvuOHMuValSpTzwjn7mJdlEILcEZJK55a/cI0vAwkHCokuXLq+//npy5zKsmFzpjIYVjvTyA4GomaQfmEsGEbCeg3z44YcZWeV8Y2YfcZNFihThrQ/BEgERyBUBmWSuyCvfKBOwdpBMOkYZinQXAb8RkEn6rUQkTxQIWDtINN+xY8eqVav27t1rKLRv395XODQo6qvikDAeEPC5SXpAQFmIgMcErB3ks88+y37HvFDVtGlTVrS2bduWl449lkzZiYAIxAjIJFUZRMB7AtaLdPCOn3zyycknnzxnzpylS5eWL1/ee8mUowiIQIyATFKVQQS8J2DtIFmbYzYH2LdvHzsGrFy50kYy9rVq1KgRfU1eHyGazp+zYaVbIpAegZRMMr0s9JQIiEACAWsHyRasO3fuPO+88zp37tyrV6/KlSvbg6OjuWzZMvOav86fs2eluyKQBoFUTTKNLPSICIhAAgHrOcg333yTePfcc0/Hjh1/+umnbt26OQfH+XNz584l/oABAzp06DBq1ChC2Dq5RIkSNWrUqFWr1qJFi5jUdJ6gYoqACGRikqInAiKQHgHrHuRll11mkjvjjDPYBHLgwIE2qfOWJG8xc8IOp+QQTefP2bDSLRFIj0BKJpleFnpKBEQggYB1D3L58uWxeAcPHvz0009twHHCDmOwW7ZsYTyWCcvkmPbnz+FWjWdlV/TkZxUiAiIAgZRMUsREQARcIZDYg+S4nDJlynz++efsVM6Ha85fZRrSJjMzQ0m0888/n+HTVM+fy+3hc7xPqVcqbQpXt3JOIA2TzLnMEkAEwkEg0UHefvvt7NkxfPhw9pbjw/W2bdsw0YK0NWd9cJcLznRt2LChzp8riJXCRSANAqmaZBpZ6BEREAFLAtZDrD169MDhcY4HB64uWbJk6NChvBNp+TwzjnQcuXXgwIFLLrmka9eurVq10vlzlqwUKAJpE3BukmlnoQdFQAQSCFg7yGuvvfaz/M/o0aMHDRp0+eWXz5s3z5JdzZo1iRZ/6/jjj581a1ZC5DvyP5YpKFAERKBQAs5NstCkFEEERMAhgcQhVvMY53iwNpXXM+g78tFGyQ5pKpoIZImATDJLYJWsCNgQsO5BsjaHeUfGV+fPn88q1v3799skoVsiIALZJiCTzDZhpS8CyQSse5CvvPIK7/VPmDDhxBNP3LhxI2t2kp9UiAiIgGcEZJKeoVZGIhAjYN2DxC8OGzbMRKpWrRpzkEImAiKQQwIyyRzCV9aRJWDdg4wsDikuAiIgAiIgAoaAdQ/Sz3Sy9F6/SXZtXnc/6y7ZREAEREAEPCOQ2IPs1KkTed92222eSaCMREAEbAjIJG3g6JYIZJVAYg9y06ZNvPI4bdo0zt+I30O1efPmWZVDiYuACFgSkElaYlGgCHhAINFB3nfffRzouGHDhtgiHYTgncjZs2d7IE3Os2CgVaOsOS8FCRBPIOImqcogAjkkkOggL8z//OUvf7nrrrtyKJayFgERMARkkqoJIpArAokO0siBd2SUlV0C+JdDj9kHMlfyKV8REAEIyCRVDUTAewKJi3SMBBwgMHbs2Pr5Hy7413vJlKMIiECMgExSlUEEvCdg3YN85513li1bVrTof93ngAEDmjVrZnPilfdCK0cRiBoBmWTUSlz6+oGAdQ8SyXbu3Gnk++mnn/wgqGQQgYgTkElGvAJIfe8JWPcgGc+h19ixY0fe9GAmUt1H7wtGOYpAPAGZpOqDCHhPwLoH2a9fv4ULF/bO/3z88ce8E+m9ZMpRBEQgRsChSe7du7d169ZNmjRp0KDB3XffzePbt2/v3Llz7dq1+d6xY4dJkCZvrVq16tSpM2PGDEEWAREoiIB1D5LYlSpV6tmzZ0GPKVwERMBjAk5MkkN4eGW5dOnSHFHXrl27bt26vfHGG+zFM2LECN5v5jNq1KgVK1ZMmjRp+fLlP/zww1lnnfXNN98UK1bMY12UnQgEgoB1DzIQoktIERCBBALs6YF3JBAHyccce846O0L4njJlChcchM6YEK60Ro0a9CMXLVokjCIgApYE5CAtsShQBIJKgBPOmzZtWqFCBcZU27Rps3nzZrqeKMP3li1buOCE16pVqxr1qlSpwr9BVVVyi0CWCVg4yEOHDjVs2DDL+Sp5ERABpwRSMknGS3lHi90i6Rp++eWXyXnE77HMXXqZCXHGjx/fMv+zdevW5McVIgLRIWDhIHn9kUn+77//3jkFGq2sejUb7mhRgHNuiikCTgikYZLHHnsse2BNnz69YsWKbHdOLnzTreSCXuP69etNvvjRypUrJ8gwePDgxfmf8uXLOxFPcUQgrAQsHKSxJVbBMbfPOh3zsdef3Xbq1atn4rAQgAdXrVrFN9cExhYFYK5DhgzBm9qnprsiIAIJBHBvTkySPp95XXLPnj0ffPBB3bp1sd+JEyeSGt+9evXighAW6ezbt2/NmjXYKateRVsERMCSgPUqVrNA3PKB5EAaoWzzcccddzz88MPcZQnA3LlzuWBRAG1YVs0lLwpo27ZtcjoKEQERKIiAQ5PEj2J3tEEZle3Tpw+DOtgaFxMmTKhWrdrkyZNJH0dLCLtIFi9efNy4cVrCWhBzhYuAtYM844wz1q1bR+uSVeC//vqrfZ/vpptuGj169K5duwxNy0UBp556qrmrRQGqcyKQBgGHJtm4ceOlS5fGp3/88cfPmjUrIUeas3zSEEOPiECkCFgPsT7zzDMcsnP11VfDgkVu5513XkFQ3n77bSY2WrRoUVAEwu0XBfhtRQBHQvJno45uiYD3BJybpPeyKUcRCCsBawfJwMuCBQvKli2L2uzBYVaHWyIgGgdjVa9enTereEP50ksvTXVRgFYEWIJVoAjEE3BukuImAiLgFgFrB8lLxEceeaTJ48CBA8kLwWPZs2cVc5Br165l2v/MM8984YUXtCjArbJROiIQI+DcJAVNBETALQIFzkHef//9LISbOXPmE088ce655zrPj02ttCjAOS7FFAEnBJiDTNsknaSvOCIgAskEiiRMEJoYLIFj2dv777/P3bPPPvuqq66y6UQmJ5peCC8m8+qV/bNezg6uzetuL4zuioArBJzUfN+apD0BDFZ2ZI9Id/1GIN4erXuQvJjMYnH2qcIvsuW/B97Rb4wkjwj4ioBM0lfFIWEiQsB6DpL3Gk855ZQbb7zx+uuvZzvj9957LyI4pKYI+JNAREzSyyEifxa0pPIVAese5M033zxnzhxcI7KuXr26e/funJvjK7kljAhEioBMMlLFLWV9QsC6B8mrjcY7ImXNmjXNFo4+kVhiiEAECcgkI1joUjnnBBJ7kByvikxsRnXOOeewGJXZR7anatWqVc4FlQAiEE0C0TRJM9aqBT7RrPP+0TrRQb711ltGON73nzdvHtfs6L9jxw7/SCxJRCBSBKJmklr4Gqnq7XNlEx3kc88953OJJZ4IRIqATJLiTvCacqKRMoEcKpvoII0onIPz2GOPsT8O2+iYEPaTy6GUyloEIk5AJhnxCiD1c0LA2kGyO/mgQYPYQIe3r3IiljIVARGIJyCTVH0QAe8JWDvIkiVL8hKk99L4KkctE/BVcURcmEiZZPzbkPGjqRpZjbgVeK++tYMcOnTovffe26VLF7ZINjI1b97ce+GUowiIgCEgk1RNEAHvCVg7yC+++OL555/n+CozxMrLHlx7L5xyFAERMASibJLxHUrVBxHwkoC1g3zzzTe/++672IlXXgqkvERABJIJhNUkNWqaXNYK8Q8B6zU4TZo02blzp3+klCQiEHECoTfJNLqJaTwS8Vok9VMlYN2D3Lx5c926ddlAJzYHqdc8UiWr+CLgIgGZpIswlZQIOCRg7SBZoePweUUTARHwgIBM0kC26TXajNba3PKg7JRFcAlYO0iOLw+uSpJcBMJHIJomaeMOw1fE0siHBKwdZJkyZcwhyb/99tv+/ftLlSr1888/+1B6iSQCESEgk7Qp6AQ/qv6iDSvdSomAtYPctWtXLJUpU6YsWrQopUQVWQREwF0CITZJdRPdrSpKzUUC1qtY4zNgjyu9BOkicSUlAhkSCKtJZslTZinZDAtRjweCgHUP0hxBhwKHDh1avHixGW611Gfv3r3t27fft28f25pfeOGFLCXYvn37xRdfzEbn1atXf/XVV4877jgeHDly5IQJE4oVK/boo4+effbZlkn5M9BYl86l82fpREcq5yYZESaWbi9+cFWWG5GakFU1rR1k7Ai64sWL4+emTp1akBC8B0L/snTp0kxVtmvXrlu3blhyp06dRowYkZf/GTVq1IoVKyZNmrR8+fIffvjhrLPO+uabb/CUBSWocBEQgWQCzk0y+VmFWBLQVKUlFgXGE7B2kM6PoKNziXckRRwkH/7Fm86dO5eQAQMGdOjQAQdJSN++fXGlNWrUqFWrFjOabdu2DVYxqDUarPIKn7TOTTJ8uhuNLLuMyco6jJb8oEJEIJlAooO87777kiPh9u66667kcBNy8ODBFi1afPvtt9ddd12bNm14o7lSpUrc4nvLli1cbNy48dRTTzWRq1Spwr8FJaVwERCBBAJpmGSkGMojRqq4PVY20UHyRke8BLt372bucNu2bTYOkvHSZcuWsTXd+eef/+WXXyYrcPjw4fjAhBnN8fkfImzdujX5WYWIQMQJpGGS0SEm7xidss6JpokO8uabbzZysKx87NixDOwwOhoLtBHx2GOPZUB1+vTpFStW3LRpE91HvitUqMAj9BrXr19vnt2wYUPlypXj0xmc/yGkZcuWNunrlghEk0DaJul/XHJv/i+jiEto8ZoHy1DvvPPOxo0bszB1yZIlTCIaP2dJim6f2dZ8z549H3zwATu49uzZc+LEiUTmu1evXlwQwiIdVrquWbNm1apVrVu3tkxKgSIgApYEUjJJyxQUKAIikAaBRAc5fPhw9ihn2w7On7vnnnvMSxo26dJN7NixI96Upzp37tyjRw/Wr86cObN27dp8c82zDRo06NOnT/369bt27Tpu3LjgLmFVg9emJuhWlgikZJIM1WCP9erVw+gYAUIknCuGiT3yvWPHDiMkr12xXK5OnTozZszIktg+STZmszJen5RIsMQokjBByAnJLDfl7Y7YTCERuPZgqzmGWHnn0h5fzmu5Xoi0LyDdTYOAfc1PySRpsPJp3rw5UyQsnWMbrH/84x/lypUzr13hIM1rV/369WMxeaGvXdkL5kRTDNbGZLJnzmSanHiCJPayOdFOcUJJIL7aJ85BsjNAKHWWUiIQUAIpmSRz/2YNOYNA9CNZMR7K164KLcpk72jziDylDZyI30ocYo04DqkvAuEgwFZWS5cuLei1q6pVqxo1o/baVbLjTA4JRwWQFq4QSOxBupKoEhEBEcghgV9++eWCCy4YM2ZM2bJlk8Wwf+2K+HrzKhmaQqJJQD3IaJa7tA4tATa0wjv279+/d+/eKGleu+LC4WtXxOS1K1YD8ClfvnxoMUkxEXBAQA7SASRFEYGAEKB3OGjQIGYfhw0bZkT252tXGtgMSIWKupgaYo16DZD+YSKwYMGC559/vlGjRk2bNkWv+++/n/WrvGTFfljVqlWbPHkygbHXrlisnpPXrox3zK2PTMjd/GuWuWrNTpgsIkNd5CAzBKjHRcBHBDhRJ2GKEeFmzZqVIOId+R8fye2JKLl1yZ6oqExcJiAH6TJQJScCIuBnAnKTfi4dv8mmOUi/lYjkEQER8IiA5UCrR3krmyAQkIMMQilJRhEQAW8JqKPpLW+f5iYHmVrBYDaynNSQKbYIBIdAQdZdUHhwNJOk6RDQHGQ61PSMCIhAGgTkZtKApkdySEA9yBzCV9YiIAIiIAL+JSAHmWbZaKw1TXB6TAR8T0A9Xd8XkUcCykF6BFrZiIAIhIyA/GjICjRZHTnIZCYKEQEREIFEAnKHiUQi8L8cZDqFLFNJh5qeEYEQEdCPQIgKs0BVtIq1QDS6IQIiEGUCMRdo9miNMorI6q4eZEZFr1ZkRvj0cHgJaBVbeMs2QprJQUaosKWqCIhAJgTUIM6EXhCfzdRBrl+/vmPHjpw/xxk6Y8eOBcH27ds7d+5cu3Ztvnfs2GGgjBw5slatWnXq1JkxY0YQMdnIrJayDRzdEgEREIHgEsjUQXKk3EMPPfTVV18tXLiQs+VWrFiRl5fXqVOnVatW8c01aAicNGnS8uXLp0+fPmTIkIMHDwaXlyQXAREQARGICIFMHWSlSpWaN28OrDJlytCP3Lhx49SpUwcMGEAI31OmTOGCkL59+5YoUaJGjRr0IxctWhQRuFJTBERABEQguAQydZAxzdeuXbt06dI2bdps3rwZr0k431u2bOECr1m1alUTs0qVKvwbXF6SXAREQAREICIE3HGQv/zyywUXXDBmzJiyZcsmg0s44rxIkSLxccaPH98y/7N169bkZ4MSotn7oJSU5BQBERABhwRccJD79+/HO/bv3793797kWrFixU2bNnHBd4UKFbig18haHiPQhg0bKleuHC/c4MGDF+d/ypcv71BoRRMBERABERCBbBPI1EHSOxw0aBCzj8OGDTOy9uzZc+LEiVzz3atXLy4IYZHOvn371qxZw+Kd1q1bZ1srpS8CIiACWSKQvHA9NoCkkaQsMc9VspnupLNgwYLnn3++UaNGTZs2RYf7779/xIgRffr0mTBhQrVq1SZPnkwgb4AQUr9+fZa8stK1WLFiudJW+YqACHhMAJ8Rpp1okl1gcojHhJVd9ghk6iDbtWuXMMWIrLNmzUqQ+I78T/bU8E/KIfs58A9YSSICIiACHhPIdIjVY3GVnQiIgAh4TEB9RI+B+ye7THuQ/tEk55LIinJeBBLAzwRkIH4uHclmSUA9SEssChQBERABEYg6ATnIqNcA6S8CIiACImBJQA7SEktGgRpKygifHg4jgdAbBQqGXscwVsxCdJKDLASQbos3i7m4AAAVK0lEQVSACGRIQJ4jQ4B6PFcE5CCzQl7NyaxgVaIi4DMCBfn+5PDkEJ+pInEsCMhBWkBxK0gm4RZJpSMCIiAC3hPQax7eM1eOIiACoSVg2Sy2DAwtghApph5kdgtTY63Z5avUkwgMHDiQQwIaNmxo7mzfvr1z5861a9fme8eOHSZw5MiRnMxap06dGTNmJCWgABEQgf8jIAepqiACoSJwxRVXTJ8+PaZSXl5ep06dOCSAb64JX7FiBYcHLF++nGhDhgw5ePBgqPT3jTLqNfqmKNIXRA4yfXbOn5SpOGelmBkSaN++fbly5WKJTJ06dcCAAfzL95QpU7ggpG/fviVKlKhRowb9yEWLFmWYox4XgbASkIMMa8lKLxH4L4HNmzdXqlSJC763bNnCxcaNG6tWrWrocFYr/4pUlghohiVLYD1LVg7SM9TKSAR8QSDh+J0iRYokiDV+/PiW+Z+tW7f6QuKwCKGRpMCVpBykR0WmtqRHoJXN/xKoWLHipk2bCOObxTtc0Gtcv369ibVhw4bKlSsnMBs8ePDi/E/58uWFM3MCBfnFgsIzz1EpuEVADtItkkpHBPxIoGfPnhMnTkQyvnv16sUFISzS2bdv35o1a1i807p1az/KLZlEwAcE9B6kDwpBIoiAewT69es3d+7cH3/8kZ7ivffeO2LEiD59+kyYMKFatWqTJ08mnwYNGhBSv3794sWLjxs3rlixYu5lrpREIFQE5CBDVZxSRgRefvnlBAizZs1KCLkj/yNW3hNgWHVtXvdYvmaUNT7Ee5GUow0BDbHawHH/VmzWQdMP7sNViiIQBAKy/SCU0v/JKAcZoMKSqCIgAlEhID/qh5KWg/RDKURaBn4I9FsQ6Rog5X/3O0sTsAwULS8JZOogtfGjl6UVvrz0ExC+MpVGmRAwFiG7yIShi89m6iC18WPahSEbSBudHhSBQBNwPmpiYuq3IlfFnekqVjZ+XLt2bUx6tnlkiTn/svFjhw4dRo0albzxY9u2bXOlrR/yja/r5po1bLELP0goGURABLwkYMxfXtBL5g7zyrQHmZCNNn50yF3RREAEREAEfE7AZQeZrK02fkxmohAREAERSJWAfRfT/m6qeSm+IeCyg9TGj6pYIiACIuABAXlEDyC77CC18aMHZaYsREAEokkg2SkSkhwYTTjZ0DrTRTra+DEbpaI0RUAERCCeQLwXlEf0rG5k6iC18WPmRRWr7lxoV8bMeSoFEYgIAXnKbBe0y0Os2RY3Iumr3kekoKWmCCQTyND8zePOE3EeM1nU0Idk2oMMPSCPFVRl9Ri4shOB0BBI/vUgRINSmZSvepCZ0Mvis9Ts5OqexfyUtAi4R0BV1z2WiSlZsrUMNE/qlySRYCr/y0GmQktxRUAERCAsBOLdqo2LDYu66eghB5kONS+fUQPQS9rKSwRCQMB4uzR8XuzXptBnI/K7JAcZAnMIqgqFGmFQFZPcIuAzAmn7s4gbqRykzyry/4pjWTstA32thgPhjAGHUjUH2iuKCHhBwN6+7O96IZ//8pCD9F+ZWEmkumtFRWEiIALpEDDtUZ60/GGxDEwnm+A/IwcZmDK0qbU2twKjngQVARHwhIDznwsTM+ZNPZHOX5nIQfqrPOyliW/3Oa/l9mnqrgi4SEDV0kWY2UjKsoAsA5NzL9RTmgjxqTlMOTkvn4TIQfqkINIUo9Aqm2a6WX4soGJnmUpIkg/6b2JIisFtNSzdnkNDjkUziQSohshBul2PcpReQhXMkRTKVgREIHIELN2ex14wS9nJQQa4NhdUJ+KdZUFxcqi2D0XKIQ1lLQL+JFCQnRYUHq+Fkzix+JaRLQNTAhX7GUzpqYTIcpCZ0PPds5nXKt+pJIFEQARCQcD8OsV/x9SydGZ++DWTgwxF1ctACT/UwgzE16MiIAI+JWD/22KcYnyc2LXlhX1qNgjSfpA05SBtwIbnVnJd9LNuRlo/SyjZREAEHBJIyT8l237y4/HuM/luglQmQqHRCtJFDrIgMoEPt6kTNrdyrrafZcs5HAkgAhEkUOhvQoZe0AapzoO0gROSW/HVq9CqllWdc5t7VlVT4hSujh6MbDXwxrRjudhUtoIkiYWnVEvVg4xslS5c8fiqxrX5K/yx/BjxVTmlBx2mr2g+JFDQb5MPRZVIQSfgTWVTDzLo9SRN+ROcn0kluW1loiWExwIT6ijR4uMn3E21QsfHTxYsTbX1mAiIQLgIWP6UuaWidz3I6dOn16lTp1atWnl5eW5Jr3TcJUBVM382yabq52yS0q0cEpA95hC+sg4KAY96kAcPHrzuuutmzpxZpUqVVq1a9ezZs379+kFhFEE5nXT+kj1lLCT5VoYMTYJZ7Ud6kEWGEFx8XPboIkwlFWICHjnIRYsW0XesWbMmKPv27Tt16lQ5yBDXqiypVpDfNY6TuwkXCWLEXGCkfKFlWcgeLbEoUAQSCHjkIDdu3Fi1alWTN53If//73yoJEXCLQHLPNd6VxrymyS7+VoLLTIhgOfNq4sRmW/k3uV8byyL5llsqZ5iO6/YYTzVD2fS4CPiHgEcO8vDhw/E6FylSJPbv+PwP/65cubJly5aFoPnxxxNOOKGQOL65/aOkzVpZOGfbsuXdSJFcaQg3gbELoiVcx4sfn4JNtPi8TNaE2Eu7bt26rHGyTtjGHnnAuUnG9ErGa51xTkPtSyGnollkLmktoLgRhFXas423R48cJL3G9evXG+02bNhQuXLlmKaD8z8OFceDLl682GHknEeTtNkrArHNhK2NPZKsc5NUKWRSCvbPiq09n0zuOmfr0SpWFuasWrVqzZo1v/3226RJk1ikk4l6elYERCATArLHTOjp2egQ8KgHWbx48ccff/zss89m+dzAgQMbNGgQHcTSVAT8RkD26LcSkTz+JFDsnnvu8Uay2rVr33DDDUOHDm3fvn0mObZo0SKTxz1+VtJmD7jYZsJW9pgJPW+eVQ3PHmeHbIskTNdnTyClLAIiIAIiIAIBIuDRHGSAiEhUERABERABEYBAYBxkgHbGYr1ux44d69Wrx1Tr2LFjA1HPmBtu1qxZjx49/C/tzp07L7zwwrp160L4448/9rnAjzzyCNWgYcOG/fr127t3r8+lTUm8AJkk6x4qVKhAKaSkYK4iB+sHhFrdunXrJk2aUM/vvvu/r1T5/JPabx1DrP7/HDhwgF14Vq9evW/fvsaNGy9fvtzPMv/www+ffvopEv7888/M9PhcWkPyoYce4he8e/fufgZrZLv88sufeeYZrqkMO3bs8LPAvNFUvXr1X3/9FSEvuuii5557zs/SpiRbsExy3rx5mCS/4CnpmKvIwfoBOXTo0K5du2DFGwp4StqsueLmMN+UfuuC0YOM7Yx15JFHmp3q/NxIqVSpUvPmzZGwTJky9HLYtcTP0iIbv+PvvPPOVVdd5XM5EY82x/z58wcNGsQ1leHYY4/1ucw4kj179vCNm4x//dfnYhcqXrBMkoWB5cqVK1Qpn0QI1g8Iu76ULl0adPvzP/GbwPiEZ7wYqf7WBcNBJuyM5X+XY4pk7dq1S5cubdOmjQ8rSrxIN9100+jRo4sWDUBl+O6778qXL3/llVcyIIxH3717t5/ZnnTSSbfccku1atX4yTvmmGO6dOniZ2lTki2gJpmSjjmPHJQfEAYtmzZtyiB2586dff5zl+pvXQB+E6mm9J3jK6vPGylG1F9++eWCCy4YM2ZM2bJlc25pNgK8/fbb1GyHi55t0vHmFl2xJUuWXHvttbQ8SpUq5fOj0xgBZl9+9sdg0Axf/sILL3hDyYNcgmiSHmBxMYug/ICgcrFixZYtW0bnjHGFL7/80kUI7iaVxm9dMByk/c5Y7kJ0JTWGGvCO/fv37927tysJZi+RBQsWTJs2jakyxq5nz5596aWXZi+vzFOmJvAxrVSW6uAsM08zeyl88MEHNWrUoMt7xBFHUBM++uij7OXlccqBM0mP+WSYXYB+QGKaMt/RoUMHlm5lqHv2Hk/jty4YDjJYO2PRuGaSjNnHYcOGZa+w3Up55MiRNP0YzGELwDPPPNPnvZwTTzyRY2HY1x71Z82a5fND0xhcXbhwoVmkg7RUCbdKLefpBMskc44rJQGC9QOydetWFpajIHPttAhZXp6Ssl5GTuO3LhgOMrYzFj8xffr08flOdbRTnn/+eXpjjMvzeffdd72sBKHP67HHHqNrzmJmRnX+9Kc/+Vlferp0c1mx1ahRIxb7Od+U389KGdmCZZKs0G7bti3tKjq+EyZM8DneYP2AbNq0ibfasEfaTMxBBuJVMecVQDvpOGelmCIgAiIgAhEiEIweZIQKRKqKgAiIgAj4g4AcpD/KQVKIgAiIgAj4jIAcpM8KROKIgAiIgAj4g4AcpD/KQVKIgAiIgAj4jIAcpM8KROKIgAiIgAj4g4AcpD/KQVKIgAi4SoDtXXjJilfCOGji4Ycf5jUbkl+8ePGNN96YnA/vAb/00kvJ4W6FPProo7yixutJriTI+/h16tQxb5HxHlGhafJCVOxls3vuuefBBx8s6BEDjY2fiGB2WC0opmX4K6+8UqtWrTC96SEHaVnQfgk09ZVjejgLgvfNUxULszdH/MydOzeNWvuPf/yDXWBc2cT8ww8/5KX+oJw3lCpnxfchgaOOOgrHwFk6M2fOxD3ce++9CNmyZUt8VbK02XaQTzzxBDK8+OKLsazZNDFZDOchJIV2fF577TX7p8go3kHaRzbQ0t5V/+KLL3722WftswjWXTlIX5eXqa9sb8jJFU899ZSNrBnaW0Epu1XjTz/9dG2YUBBkhWeVAFsNjx8//vHHH2eHmlhLkQOwTA+MXe85rWnEiBG04Qjh/E6cJdWV7R34mN0BeYpOmzmFlF6g2Yf2k08+Oe200+iecsYTKbBh9/Dhw3lZnlfmn3766XiNrrnmGjbZ79mzJ4nTgWO/CLat59S2devWderUifh8f//99zxyxRVXsM8w791zuh8Scowl/U4CC+Xz1ltvsSsFupx11lmbN28mfnxGf/7zn+nboR3f3FqxYgXqkIVlWyEhrx9//JE9FjjtBwgcinL++efT0kUj0yNnYzkoAQEVChUykBEcnqGlaDkhwH7cJt8nn3wSy2HTVKyRik51/M9//sMtTij94x//yAYW7BXCptjt2rXDSPiwGQd3CTFn4M2ZM8ec9fjvf/+b6k4KfH/99dfE4XQ0Dsfg4vPPPycye2rHNOX8wuuuu878y3WvXr3ohrJrK3vZ8BSJYJPbtm0jAj9AtM0x9f+vvXvXxS2I4gDenIRoSZQK4RkkCiUKjyBKiRDR0NC4FKJSaCWUxKWQCLUXEIWEN1BIvMH5JSuZ7Oz97c1xcGxnvkLGzOy5rG+v+a/bN8uNozHC0dGR0dQ4a2KEtJg0fi5kCnweBRLvxBRuCsUyiRG8yTc3N5pgm4tPU70aL7Bb0xQeHh5c4q+gVcoBeYyh4MjICCiVi9Qtu+7m1vry8mIEoLixseFfCYQ9BRGLWxsYGHAlmxoMC1Hi9kFrYKRR6XIfzKUwMzNDJIU95+fnkuVhSTPq72r+4mhjY2PDw8NgyUe6GE3Pz8+eUpAq1Q2XpYmKjGwBeN8irUcKMEkciyMXiaaMYg6c6+trfRChq6tLUl6yOBg+Pj5+enpyM1HsNM6B6NaKtLLFXTeUf7US1f+zRXsjLy8vJyYm4J+7PSUzYceQoApKoYRMsFidronrWJO6u7sfHx/hJXdLlU5uSpRP0T1hbk10T9vJyYn8L8TJs7Ozra0tTN7T01N9KmooshgVa3EzbG9vKy8tLR0eHhoBLsJp3VZXV3H7wsLC+vr61dWVfE9xT2PdmLk+U+BrKOAQLE40OjoKSCKdgFO+2ATt5ufnmSU5OGBkNMGJ6EYupGJKXiaFGX1Ra6TrgSLwLAyeIBMPQtCOW6NK4lZNcgufnp4qTE9PLy8vR+epqSkM7m7C/v5+f1USNM1o3uJoTKxE0lTjOmXI6tY3gJfmTROVlgHAQJ0P3Zq6Wdp+6owOBPG9vT14HJWIQO9Udrw4c4xAp4zpWpRus0SN5n8zQDbT5x+3kmSDMehhLkB3mWQDG3Rk7NIGsC4pFfdiQv21SgNJjKXqzc7OOjUaNszyQ6r1cTpgYz0xsENBAXaCRlgoR8/4+LgaQzENuTj3++czadhybvoZFKDlQDt4cH9/HztiU4UTzP40QsJicZsMocDp9vaWTkbcjCZgEAXjEFjBbSnpnhqWlXj5m4lGM6t2SKPFRLgyzaj8qgOFSArvISJDKONqjN9xIk2lvVQXEzXEaKowMTcBZHHLylUi1A3V3vrsg/zW3134IAmzeI8bEhuQbe/u7qh6NLkSGyTGpjsSJDtubG1tDc7BM06LNAK8FLEWoWsdn4rKIscmNg7WhYV8PBbGgBPD8phubm6ySgF45peGYXNTpsCnUoAtkc8M4xTPd6ZC4t3Kygo9jK+B2MfWGssgRNIOwZKUAyycHdfGEoNfuCG1ehAXgEZ+kBA66Z1vSeXNhSmFjhGog4xDHSd6Y6U1s9bofHBwUH2kuLtqa10Ncu3v7yNOyrrKpMxRQm7gy7RgplqOUjVGYOOtG6fV9Rkg2/T1NbPBWxg7jUBrjJ2rWVxcZHcFY69GxNURyxnhTHE6pDg9BxAPJUNrX18fmKx7MNdnCnwSBcL6wj7JYSYohuhWnEgmczHVfHhk0MnJSRYUCpN/SZlzc3NghmYJ5+qUMNIqkCCwekQEAKFQsLfoFf5CwzLGvKrzWYwYGd5BU0Pi3d3dP6ID4zDR08fuPEhrFOjOzoTdquOQiQXm6BxBOtUOdTXUZRDO+ygKVx+ISPO2QWZV0TpC3AUfMBEhAstW3SDtrm/wT+amf06BUqAB171Xk+zGM8/uYXnYfmdnJ9aJnwnFYMlLHA9Wg3RE5Q0NDRFdWUQFDnhQhA7mVBBHNzg4yCeRdl0K0kkBOyniIHXAP4J3LImczoRrBPyDkRxPfnYW4QM5SCcRNhcyBb4nBUoHTnGRxTimhsW/sVvDCN+qKae7ard886mrp2Wy1rKdfsgsAg2E7bHufshoeZBMgUyBD6eAX0Dy1HLNVn8KybvphoGLi4uGSWmofm/Kc0knbujWoqYMkC36sr56qV53ka7sM3//41+R8SxXvb292Oyrt5HnyxTIFMgUeBcFMkC+i2z5oUyBTIFMgUyBn06BHKTz07/hvL9MgUyBTIFMgXdR4DdA+gcs8+nfIAAAAABJRU5ErkJggg=="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ✨ Exercise 3.3\n",
    "\n",
    "- Prepare a 2x2 figure with histograms for `ra`, `dec`, `parallax`, `distance`. Use `plt.subplots(ncols=2, nrows=2` and select a figsize in accordance to the number of plots.\n",
    "- You can access each axis (panel) with `ax[0,0]`, `ax[0,1]`, `ax[1,1]`, `ax[1,0]`, clockwise.\n",
    "- For labels (which will be shown with legend()) use the column description.\n",
    "- Remember to add labels and a legend for all axes.\n",
    "- Remember to include `xlabel` (automatically) and `ylabel` (manually) for all axes.\n",
    "\n",
    "In summary, obtain a figure similar to:\n",
    "\n",
    "![image.png](attachment:067abd68-8a10-4b42-b130-7252cc7e0013.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 🌪 Exercise\n",
    "Redo the previous plot show each panel shows both the Total Sample and the Cluster subset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simplify repetitive tasks with loops over axes/columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's see a basic example with the use of `zip`, which iterates over multiple lists at the same time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = [1,   2,  3]\n",
    "y = [10, 20, 30]\n",
    "\n",
    "for xi, yi in zip(x, y):\n",
    "    print(xi,yi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "col1 = 'pmdec'\n",
    "col2 = 'dec'\n",
    "cols = [col1, col2]   # We define a list of columns to plot\n",
    "\n",
    "fig, ax = plt.subplots(ncols=2, figsize=(10,5))   # We can select ncols, nrows, or both.\n",
    "\n",
    "for axi, col in zip(ax, cols):        # zip allows you to iterate over multiple variables\n",
    "    axi.hist(data1[col], bins=80)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "col1 = 'parallax'\n",
    "col2 = 'distance'\n",
    "cols = [col1, col2]   # We define a list of columns to plot\n",
    "\n",
    "fig, ax = plt.subplots(ncols=2, figsize=(14,6))   # This plot will have two panels in a single horizontal row.\n",
    "\n",
    "for axi, col in zip(ax, cols):\n",
    "    axi.hist(data1[col],    bins=np.arange(0, 7, 0.05), label='Full sample')  # This plots the blue histogram\n",
    "    axi.hist(cluster1[col], bins=np.arange(0, 7, 0.05), label='Cluster')      # This plots the orange histogram\n",
    "    axi.legend()\n",
    "    \n",
    "    description = data1[col].description\n",
    "    unit = data1[col].unit\n",
    "\n",
    "    axi.set_xlabel(f\"{description} [{unit}]\")\n",
    "    axi.set_ylabel('Number of stars')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Function to make simple plots. The `**kwargs` trick\n",
    "\n",
    "We can simplify the process even more the process by writting a function that does the repetitive syntax by itself. To be able to use different properties for each panel (for example symbol colors, sizes, labels, etc) we will use `**kwargs`, which basically translates into \"any other variable included when calling the function will be passed along\".\n",
    "\n",
    "This function will create a single `scatter` plot of two variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_variables(data, col1, col2, ax, **kwargs):\n",
    "    \n",
    "    ax.scatter(data[col1], data[col2], **kwargs)\n",
    "\n",
    "    description1 = data1[col1].description        # Obtain column 1 description\n",
    "    unit1 = data1[col1].unit                      # Obtain column 1 units\n",
    "    description2 = data1[col2].description        # Obtain column 2 description\n",
    "    unit2 = data1[col2].unit                      # Obtain column 2 units\n",
    "\n",
    "    ax.set_xlabel(f\"{description1} [{unit1}]\") # Set the axis label in the form \"Variable description [units]\"\n",
    "    ax.set_ylabel(f\"{description2} [{unit2}]\") # Set the axis label in the form \"Variable description [units]\"\n",
    "\n",
    "    ax.legend()\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "col1 = 'phot_g_mean_flux'\n",
    "col2 = 'parallax_error'\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8,6))   # We can select ncols, nrows, or both.\n",
    "\n",
    "plot_variables(data1,    col1, col2, ax, s=1, c='k',      label='Full sample', alpha=0.2)\n",
    "plot_variables(cluster1, col1, col2, ax, s=2, c='tomato', label='Cluster')\n",
    "\n",
    "ax.loglog();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now will create a figure with two panels, the left one for `parallax_error` vs `phot_g_mean_flux` and the right one for `pmra_error` vs `pmra`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "col1 = 'phot_g_mean_flux'\n",
    "col2 = 'parallax_error'\n",
    "col3 = 'pmra'\n",
    "col4 = 'pmra_error'\n",
    "\n",
    "fig, ax = plt.subplots(ncols=2, figsize=(14,6))   # We can select ncols, nrows, or both.\n",
    "\n",
    "plot_variables(data1,    col1, col2, ax[0], s=1, c='k', alpha=0.2, label='Full sample')\n",
    "plot_variables(cluster1, col1, col2, ax[0], s=2, c='tomato', label='Cluster')\n",
    "\n",
    "plot_variables(data1,    col3, col4, ax[1], s=2, c='k', alpha=0.2, label='Full sample')\n",
    "plot_variables(cluster1, col3, col4, ax[1], s=3, c='tomato', label='Cluster')\n",
    "\n",
    "# Other options specific to the panel need to be specified individually\n",
    "ax[0].loglog()\n",
    "ax[1].set_xlim(-60, 60);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Velocity distribution analysis\n",
    "We will try to find some clustering pattern in the velocity space. These are the variables:  \n",
    "\n",
    "Proper motion in right ascension is $\\mu_{\\alpha}\\cos\\delta$  \n",
    "Proper motion in declination is $\\mu_{\\delta}$  \n",
    "Spectroscopic radial velocity in the solar barycentric reference frame is $v_{r}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ✨ Exercise 3.4\n",
    "\n",
    "- Make a plot with a single panel to show the proper motion in declination vs the proper motion in right ascension for the full sample and for the cluster. The full sample table is `data1` and the cluster table is`cluster1`. Use the function `plot_variables()`.\n",
    "- Use different colors for the full sameple and the cluster.\n",
    "- Use `ax.set_aspect('equal')` so both axes represent the same scales.\n",
    "- Repeat the plot to zoom with until you see the relevant part of the plot by setting `ax.set_xlim` and `ax.set_ylim` \n",
    "- Can you approximately guess the mean proper motion of the cluster?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great, this is much clearer! The cluster is a moving group, where most of the stars have a very similar 2D velocity, even if they are positionally scattered in a wide region of the sky of several degrees. This indicates that the cluster of stars, probably very young, was formed at a similar time from a single gas cloud, so on average they have the same space velocity as the original cloud.\n",
    "\n",
    "Note that the selection based on distance alone gives a large contamination of stars that don't look like stars from the cluster. A better selection algorithm is needed to avoid those outliers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 🌪 Exercise\n",
    "\n",
    "Repeat the previous plot and include a cross on the average proper motion of the cluster, and the median proper motion of the cluster. Use values from the `cluster1` table to plot the crosses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 🌪 Exercise\n",
    "\n",
    "Make a three-panel plot to represent all the stars in `data1` and show:\n",
    "- `pmdec` vs `pmra`\n",
    "- `dr2_radial_velocity` vs `pmra`\n",
    "- `dr2_radial_velocity` vs `pmdec`\n",
    "\n",
    "- Draw horizontal and vertical lines at the center of the cluster using `axvline()` and `axhline()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 🌪 Exercise\n",
    "\n",
    "- Use a `for` loop to make a 3 panel plot showing the histograms for `pmra`, `pmdec` and `radial_velocity`\n",
    "- Include both the full sample and the cluster subset.\n",
    "- Make sure each panel has an adequate `bins` distribution, and that it is the same for the full and cluster samples. Tip: you can create a list with 3 items, each one being a `np.linespace` object.\n",
    "- Tip2: for `radial_velocity` select the `bins` that do not select value 0, for example use `np.linspace(1, 60, 50)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Velocity field\n",
    "Finally, we can plot the velocity field of the stars. Given the large number of stars it is difficult to have a clear picture of the velocity directions and magnitudes of the whole sample. But at least we can see that the great majority of the sample stars have random velocities on all directions while the cluster stars have a predominant 2D direction and sense."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(ncols=1, figsize=(12,8))\n",
    "\n",
    "ax.quiver(data1['ra'], data1['dec'], data1['pmra'], data1['pmdec'],\n",
    "          alpha=0.2, scale=1200, label='Full sample');\n",
    "\n",
    "ax.quiver(cluster1['ra'], cluster1['dec'], cluster1['pmra'], cluster1['pmdec'],\n",
    "           alpha=1, color='tomato', scale=1200, label='Cluster');\n",
    "ax.set_aspect('equal')\n",
    "\n",
    "ax.set_xlabel('Right Ascension [deg]')\n",
    "ax.set_ylabel('Declination [deg]');\n",
    "\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis of the cluster kinematics in the Galaxy\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`SkyCoord` function is much more powerful than just storing positions. In particular we can include also velocities (both in the plane of the sky and radial), or also cartesian coordinates. Let's start by simply adding the proper motions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For the full sample\n",
    "data1_coord = SkyCoord(ra  = data1['ra'],\n",
    "                       dec = data1['dec'],\n",
    "                       pm_ra_cosdec = data1['pmra'],\n",
    "                       pm_dec = data1['pmdec'],\n",
    "                       frame = 'icrs')\n",
    "\n",
    "data1_galactic = data1_coord.galactic\n",
    "\n",
    "# For the selected cluster\n",
    "cluster1_coord = SkyCoord(ra  = cluster1['ra'],\n",
    "                          dec = cluster1['dec'],\n",
    "                          pm_ra_cosdec = cluster1['pmra'],\n",
    "                          pm_dec = cluster1['pmdec'],\n",
    "                          frame = 'icrs')\n",
    "cluster1_galactic = cluster1_coord.galactic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster1_galactic[0:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are going to repeat the three basic plots (position, velocities and quiver plot) but using galactic coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,8))\n",
    "\n",
    "ax.plot(data1_galactic.l,\n",
    "        data1_galactic.b,\n",
    "        color='grey', marker='.', ms=0.5, ls='')\n",
    "\n",
    "ax.plot(cluster1_galactic.l,\n",
    "        cluster1_galactic.b,\n",
    "        color='tomato', marker='.', ls='')\n",
    "\n",
    "ax.set_aspect('equal')\n",
    "\n",
    "ax.set_xlabel('l [deg]')\n",
    "ax.set_ylabel('b [deg]');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,8))\n",
    "\n",
    "ax.plot(data1_galactic.pm_l_cosb,\n",
    "        data1_galactic.pm_b,\n",
    "        color='grey', marker='.', ms=0.5, ls='')\n",
    "\n",
    "ax.plot(cluster1_galactic.pm_l_cosb,\n",
    "        cluster1_galactic.pm_b,\n",
    "        color='tomato', marker='.', ls='')\n",
    "\n",
    "ax.set_xlim(-50, 60)\n",
    "ax.set_ylim(-60, 50)\n",
    "ax.set_aspect('equal')\n",
    "\n",
    "ax.set_xlabel('pm l [mas/yr]')\n",
    "ax.set_ylabel('pm b [mas/yr]');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(ncols=1, figsize=(15,10))\n",
    "\n",
    "ax.quiver(data1_galactic.l.value, data1_galactic.b.value,\n",
    "          data1_galactic.pm_l_cosb.value, data1_galactic.pm_b.value,\n",
    "          alpha=0.2, scale=1200, label='Full sample');\n",
    "\n",
    "ax.quiver(cluster1_galactic.l.value, cluster1_galactic.b.value,\n",
    "          cluster1_galactic.pm_l_cosb.value, cluster1_galactic.pm_b.value,\n",
    "           alpha=1, color='tomato', scale=1200, label='Cluster');\n",
    "ax.set_aspect('equal')\n",
    "\n",
    "\n",
    "ax.set_xlabel('Right Ascension [deg]')\n",
    "ax.set_ylabel('Declination [deg]');\n",
    "\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Galactocentric coordinates\n",
    "\n",
    "We will now use the full functionality of `SkyCoord` to make 3D transformations. Here we are including the radial velocity, which we now many stars are missing, so we will just focus on average values.\n",
    "\n",
    "The description of the Galactocentric reference frame in astropy can be found [here](https://docs.astropy.org/en/stable/api/astropy.coordinates.Galactocentric.html).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.coordinates import Galactocentric"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compute median values as a representative position and movement of the cluster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ra_median    = np.median(cluster1['ra'])*cluster1['ra'].unit\n",
    "dec_median   = np.median(cluster1['dec'])*cluster1['dec'].unit\n",
    "pmra_median  = np.median(cluster1['pmra'])*cluster1['pmra'].unit\n",
    "pmdec_median = np.median(cluster1['pmdec'])*cluster1['pmdec'].unit\n",
    "radial_velocity_median = np.median(cluster1['dr2_radial_velocity'].compressed())*cluster1['dr2_radial_velocity'].unit\n",
    "distance_median = np.median(cluster1['parallax'].to(u.pc, u.parallax())) # We could also have used the \"distance\" column directly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we construct the average full-coordinate object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster1_3d = SkyCoord(ra  = ra_median,\n",
    "                       dec = dec_median,\n",
    "                       pm_ra_cosdec    = pmra_median,\n",
    "                       pm_dec          = pmdec_median,\n",
    "                       radial_velocity = radial_velocity_median,\n",
    "                       distance        = distance_median,\n",
    "                       frame = 'icrs')\n",
    "cluster1_3d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we transform the 3D coordinates to the Galactocentric reference frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster1_galactocentric = cluster1_3d.transform_to(Galactocentric)\n",
    "cluster1_galactocentric"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tip: you can set your own values of the Galactic Centre distance and the Sun 3D velocity, see example in [Transforming positions and velocities to and from a Galactocentric frame](https://docs.astropy.org/en/stable/generated/examples/coordinates/plot_galactocentric-frame.html#sphx-glr-generated-examples-coordinates-plot-galactocentric-frame-py)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### ✨ Exercise 3.5 (part 1)\n",
    "- Find the Vizier catalog associated with this publication: [Open cluster kinematics with Gaia DR2](https://www.aanda.org/articles/aa/full_html/2018/11/aa34020-18/aa34020-18.html). (Tip: Search in that page the button \"Tables at CDS\", you will see the catalog code at the top of the page). \n",
    "- Download the first table of the catalog and call it `soubiran18` using the function `Vizier.get_catalogs()`.\n",
    "- Print the first few rows of the table to see what is the name of the columns for Right Ascension and Declination.\n",
    "- Construct the full coordinates of the clusters in the catalog with:\n",
    "```python\n",
    "soubiran18_coord = SkyCoord(ra=soubiran18['RA_ICRS'],\n",
    "                            dec=soubiran18['DE_ICRS'],\n",
    "                            frame='icrs')\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### ✨ Exercise 3.5 (part 2)\n",
    "- Calculate the separation between `cluster1_3d` and the coordinates in Soubiran et al. 2018. Store the separation in degrees in the variable `separations`.\n",
    "- Find which item in the catalog is the closest to the center of `cluster1`. You can find the row index with the function `np.argmin` (If you need, run `np.argmin?` to know more about it).\n",
    "```python\n",
    "idx_min = np.argmin(separations)\n",
    "```\n",
    "- Print item `soubiran18[idx_min]` and explore the table:\n",
    "  - What is the NGC name of our cluster?\n",
    "  - Compare our values to the values in the catalog (VR, X, Y, Z, vZ, etc) with the ones we found in the variable `cluster1_galactocentric`\n",
    "  \n",
    "Tip: for the comparison, print these values:\n",
    "```python\n",
    "print(cluster1_3d.ra.deg)\n",
    "print(cluster1_3d.dec.deg)\n",
    "print(cluster1_3d.distance)\n",
    "print(cluster1_3d.radial_velocity)\n",
    "```\n",
    "and also:\n",
    "```python\n",
    "print(cluster1_galactocentric.x)\n",
    "print(cluster1_galactocentric.y)\n",
    "print(cluster1_galactocentric.z)\n",
    "print(cluster1_galactocentric.v_z)\n",
    "```"
   ]
  },
  {
   "attachments": {
    "5298dde0-215c-4d62-9561-c1d07b0a28ef.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtQAAAFqCAYAAAA6Bc1oAAAABHNCSVQICAgIfAhkiAAAIABJREFUeF7tnYuSHTeORCV5P3ZiPsgxP2t5RcmU0dngs0gWWXUmdkPmC0gcsG6nylfy1y/8DwIQgAAEugh8+/bt766DGx/6/v37143lIQ0CEIDAlgT44NyyLYiCAARmEHiiAZ7BaWRMDPpImsSCAAR2JYCh3rUz6IIABLIEMMdZPEcvYsKPbh/iIfBKAhjqV7adoiGwJwFM8p592VkV5nvn7qANAu8hgKF+T6+pFAK3EsAs34r/1ckx3a9uP8VDYAkBDPUSzCSBwLMJYJaf3d83VIfpfkOXqREC8whgqOexJTIEHkUA0/yodlJMAwHMdgMstkLgpQQw1C9tPGVDwCOAafaoMAeBNAHMdpoNKxB4EwEM9Zu6Ta0Q+EEA08w1gMAaApjtNZzJAoEdCGCod+gCGiAwgQDGeQJUQkJgAAGM9gCIhIDAZgQw1Js1BDkQaCWAcW4lxn4I7EkAo71nX1AFgRoCGOoaSuyBwCYEMM+bNAIZEFhEAJO9CDRpIHCRAIb6IkCOQ2AWAczzLLLEhcDZBDDZZ/cP9c8kgKF+Zl+p6jACmOfDGoZcCGxGAJO9WUOQ8zoCGOrXtZyC7yaAeb67A+SHwDsIYLLf0Weq3IMAhnqPPqDioQQwzw9tLGVB4FACmOxDG4fs7QlgqLdvEQJPIoCBPqlbaIUABDDY3AEIjCGAoR7DkSgvJYCBfmnjKRsCDyWAwX5oYylrOgEM9XTEJHgSAQz0k7pJLRCAQIkABrtEiHUI/CKAoeYmQCBDAAOdgcMSBCDwOgIY7Ne1nIIrCWCoK0Gx7R0EMNDv6DNVQgACYwhgsMdwJMr5BDDU5/eQCi4QwEBfgMdRCEAAAkIAg82VeCsBDPVbO//SujHQL208ZUMAArcQwGDfgp2kNxDAUN8AnZRrCWCi1/ImGwQgAAGPAObao8LcUwhgqJ/SSer4TQADzWWAAAQgsD8BDPb+PUJhPQEMdT0rdm5MABO9cXOQBgEIQKBAAHNdAMTy9gQw1Nu3CIEeAQy0R4U5CEAAAs8ggMF+Rh/fVAWG+k3dPrxWTPThDUQ+BCAAgQ4CmOsOaBxZTgBDvRw5CVsIYKJbaLEXAhCAwLMJYK6f3d+Tq8NQn9y9h2rHRD+0sZQFAQhAYCABzPVAmIS6TABDfRkhAUYQwESPoEgMCEAAAu8kgLl+Z993qhpDvVM3XqYFE/2yhlMuBCAAgQUEMNcLIJPiEwEM9SckTMwkgImeSZfYEIAABCBgCWCuuQ+rCGCoV5F+cR5M9IubT+kQgAAENiGAud6kEQ+VgaF+aGPvLgsTfXcHyA8BCEAAAikCmOsUGeZ7CWCoe8lxziWAkXaxMAkBCEAAAhsSwFhv2JRDJWGoD23cTrIx0Tt1Ay0QgAAEINBDAHPdQ40zkQCGmrvQRQAT3YWNQxCAAAQgcAABzPUBTdpMIoZ6s4bsLgcjvXuH0AcBCEAAAqMIYKxHkXx+HAz183t8uUJM9GWEBIAABCAAgcMJYK4Pb+Bk+RjqyYBPDo+RPrl7aIcABCAAgRkEMNYzqJ4fE0N9fg+HVoCJHoqTYBCAAAQg8GACmOsHN7exNAx1I7CnbsdIP7Wz1AUBCEAAArMJYKxnE94/PoZ6/x5NVYiRnoqX4BCAAAQg8CICGOsXNVtKxVC/sPeY6Bc2nZIhAAEIQGApAcz1Uty3J8NQ396CdQIw0utYkwkCEIAABCAQCGCs33EPMNQv6DNG+gVNpkQIQAACENiaAMZ66/ZcFoehvoxw3wAY6X17gzIIQAACEHgnAYz1M/uOoX5YXzHRD2so5UAAAhCAwGMJYK6f01oM9UN6iZF+SCMpAwIQgAAEXkcAY31+yzHUh/cQI314A5EPAQhAAAIQ+IcAxvrcq4ChPrR3GOlDG4dsCEAAAhCAQIEAxroAaMNlDPWGTclJwkjn6LAGAQhAAAIQeA4BjPU5vcRQH9IrjPQhjUImBCAAAQhAYDABjPVgoBPCYagnQB0ZEiM9kiaxIAABCEAAAucSwFjv2zsM9aa9wUhv2hhkQQACEIAABG4mgLG+uQFOegy1A+XOKYz0nfTJDQEIQAACEDiHAMZ6n15hqDfpBUZ6k0YgAwIQgAAEIHAYAYz1/Q3DUN/fgy+Y6Q2agAQIQAACEIDAwQQw1fc2D0N9I3+M9I3wSQ0BCEAAAhB4IAGM9T1NxVDfwB0jfQN0UkIAAhCAAAReRABjvbbZGOqFvDHSC2GTCgIQgAAEIACBLxjrNZcAQ72AM0Z6AWRSQAACEIAABCCQJICxTqIZsoChHoIxHQQznWbDCgQgAAEIQAAC6whgquexxlBPYouRngSWsBCAAAQgAAEIXCKAsb6Ezz2MoXax9E9ipPvZcRICEIAABCAAgXUEMNbjWGOox7Hk75MeyJJQEIAABCAAAQjMJ4CpHsMYQz2AI2+lB0AkBAQgAAEIQAACtxHAWF9Dj6G+wA8jfQEeRyEAAQhAAAIQ2I4AxrqvJRjqPm58vaOTG8cgAAEIQAACENibAKa6vT8Y6kZmvJVuBMZ2CEAAAhCAAASOJICxrm8bhrqeFW+lG1ixFQIQgAAEIACB8wlgqut6iKGu4MRb6QpIbIEABCAAAQhA4LEEMNb51mKo83x4K13gwzIEIAABCEAAAu8ggKlO9xlDnWDDW+kEGKYhAAEIQAACEHg1AYz15/ZjqD8z4a20w4QpCEAAAhCAAAQgEAlgqj/eBQy14cFbaT4oIAABCEAAAhCAQD0BjPUvVhjqf+4MZrr+4WEnBCAAAQhAAAIQiAQw1V++YKh/3AbMNB8KEIAABCAAAQhAoJ/A2031qw01Rrr/weEkBCAAAQhAAAIQUAJvNdavNdSYaX0EGEMAAhCAAAQgAIHrBN5oql9pqDHT1x8WIkAAAhCAAAQgAIEUgbeZ6lcZaox06tozDwEIQAACEIAABMYTeIuxfo2hxkyPf0iICAEIQAACEIAABEoE3mCqX2GoMdOlq846BCAAAQhAAAIQmEfg6ab60YYaIz3vwSAyBCAAAQhAAAIQaCXwVGP9WEONmW694uyHAAQgAAEIQAAC8wk80VQ/0lBjpuc/DGSAAAQgAAEIQAACvQSeZqofZ6gx071Xm3MQgAAEIAABCEBgHYEnmepHGWrM9LqHgEwQgAAEIAABCEDgKoGnmOpHGGqM9NXrzHkIQAACEIAABCBwH4HTjfW3+9CNyYyZHsORKBCAAAQgAAEIPJPA97/+3L6w0/3c0W+oT4e//e1GIAQgAAEIQAACRxOwZvrbH//dvpZT31Qfa6gx09s/EwiEAAQgAAEIQOBGAt6baUz1nIYc+ZUPzPScy0BUCEAAAhCAAAT2JuCZ5JRiNc86Tp27e/5En3fcG+oTId99MckPAQhAAAIQgMD5BHq/vhHOnWKmbZdO+vrHUYYaM33+hwEVQAACEIAABCDQTsB7M32iSW6t/BRTfYyhxky3XkH2QwACEIAABCDwJAK9b6hPZ3CCqT7CUGOmT38U0A8BCEAAAhCAwAgCp35942rtu5vq7Q01ZvrqFeQ8BCAAAQhAAAIQOJ/AzqZ6a0ONmT7/8lMBBCAAAQhAAAIQGEVgV1O9raHGTI+6esSBAAQgAAEIQAACzyGwo6ne0lBjpp9z6akEAhCAAAQgAAEIjCawm6nezlBjpkdfOeJBAAIQgAAEIACB5xHYyVRvZagx08+77FQEAQhAAAIQgAAEZhHYxVRvY6gx07OuGnEhAAEIQAACEIDAcwnsYKq3MNSY6edeciqDAAQgAAEIQAACswncbapvN9SY6dlXjPgQgAAEIAABCEDg+QTuNNW3GmrM9PMvNxVCAAIQgAAEIACBVQTuMtW3GWrM9KqrRR4IQAACEIAABCDwHgJ3mOpbDDVm+j2XmkohAAEIQAACEIDAagKrTfW31QWSDwIQgAAEIAABCEAAAk8isPwNNW+nn3R9qAUCEIAABCAAAQjsSWDlW+qlhhozveeFQxUEIAABCEAAAhB4IoFVpnqZocZMP/GaUhMEIAABCEAAAhDYm8AKU73EUGOm975oqIMABCAAAQhAAAJPJjDbVPOHEp98e6gNAhCAAAQgAAEIQGA6gelvqHk7Pb2HJIAABCAAAQhAAAIQKBCY+ZZ6qqHGTBc6yzIEIAABCEAAAhCAwDICs0z1NEONmV52N0gEAQhAAAIQgAAEIFBJYIap5jvUlfDZBgEIQAACEIAABCAAAY/AlDfUvJ32UDMHAQhAAAIQgAAEILADgdFvqYe/ocZM73BNvnz5/tefewhBBQQgAAEIQAACENiMwGi/OtRQjxa3Gftj5EQzjak+pmUIhQAEIAABCLyawB2eZaRvHWqoX30TNileL6SON5GJDAhAAAIQgAAEIPCTwBNeBA77DvVIl8/9ukbAmuhvf/z3WjBOQwACEIAABCAAgUkEvBd/q73LiO9TDzHUmOlJt+xC2HBBV1/IC3I5CgEIQAACEIDASwns8CLwqqnGUL/08lI2BCAAAQhAAAIQ2IXA3S8CbzfUpbfTdwPa5aKgAwIQgAAEIAABCEBgXwJXTPXUP5T4hC+Z79t2lEEAAhCAAAQgAAEI7EDg0lc+cm+nd/iS+Q6A0QABCEAAAhCAAAQgcAaB3rfU3YY6Z6Yjsh2+ZH5G+1AJAQhAAAIQgAAEILADgR5TPfUrH/FvmeBvm9jheqABAhCAAAQgAAEIQGAGga431DVvp2eIJSYEIAABCEAAAhCAAARmE2h9Sz31DfXsYokPAQhAAAIQgAAEIACBuwk0v6Hm7fTdLSM/BCAAAQhAAAIQgMBsAi1vqW99Q+39TSCz4YyO/4QaRjMhHgQgAAEIQAACEHgTgSZDPfLt9BP+juon1PCmy06tEIAABCAAAQhAoPZlaIvvbTLUo1qgheh4VJ6ZcVSzjmfmJjYEIAABCEAAAhDIEcCX+HRmvQyt/g51i0v3S/g4axt96l+r94QaanrFHghAAAIQgAAEziGAP/F75f0mo8aD1nyX+pY31KHMJ/wd1U+owb9yzEIAAhCAwBsJeIbjjRxOrll7qOOTa7uqXc2zjq/Er3pDPfrt9BXBnIUABN5DIPwgGPmB9x5yVAqBdgK81WxntusJepnvTM/PltJbagx1njmrEIDATQT4gXATeNK+koD3FpPfzJ59FXpM49kVz1V/2VDzdnpug4gOAQh8JsAP989MmIHAbAL8JnY2YeKfTiBnqm/7DvXpUNEPgVYCnklsjfGW/fpmTMdv4UCdEFhJgD8XtJI2uZ5GIPuVD95Or2s3/2pmHes7MvHmp486z0UfN05BAAIQgMAcAqm31LyhnsO7KWo0W7zBbMJ2zGbtq46PKeQGoaU307C8oSmkhAAEXkdg9Wft6nwjGoqhHkHxQgy9NDq+EJqjmxBQU6jjTWQeJyM+Kzwzx7UOwRCAwEEEVn/Wrs43qhXJr3zwdY9RiMtxrCHAbJV5nboj9Jn+jumeZ6JhO4YtUSAAAQhEAqs/a1fn6+2097UP3lD30hx4jj8IMhDmxqEwfOOaoyx1PC4TkSAAgdMJeCbt9JpW6dfPVh2P1qHxdTw6X4g36n5gqGd0pyPmikvTIWvqkVGXeKpIgncRWNHb+My88dnpagqHIPBCAvGzaMVn0lPxrv6sXZlv5P1wDTVf93jqY7FPXSMv8ZWq+JC9Qu/fs5bjyN6W+oOZHtM/ouQJlO5h/jSrdxHQvun4Ll0j866qafVn7Yp8yk7HuT55Ptn9DrW3MReYNQi0EPAu7YqHRzVaHXfkVz2njr1+2lp62dKfU2/Es3RzD8/rZ+hZ/Nx5cv+eXNuqW3eFoX6P2n1DvaoQ8ryTgBosHa+goiZQxyM0zIg5QtfIGKUae3urcXU8sgZiQSBFQO+djlPnmL+PQOxR/DV+BvV+Ft1XST6z3kUdh9PeXD5q35k7YmrOnlpDjJH349Mbat5Oa5vOHIfLtfsHyN0a7QM4gpWtZ3TsnW+hV+uI3npxd+aAtmcS4B6e01fbq6h6xGf7rgRydzO3lqqn50wqVpyfEVNzrsihOX/X9/37bx+NoU5RcuZD00Y+nKPjrbzADp4jp0b1wD7QHoiR98aLf/fcKI5ax6y4mocxBHIEuIc5Onut3Wmu7iDh3U3v51HpZ1DPmVK9M2JqzhU5NKcd2699YKhzpMza6Id0dLwoddbl8h7aSnSP3+Yxt0WXPsgeD4gCIQABCCwkwM+rj1/dqP0ZNNqX6M/GWh2tV2W07pb81lB/+A717l/30Oa0FH1lr+bVcWtsPa/j1nh2v15YHffEjvpG6uzRYc/spEUZh3Gc07WrdXMeAhCAAATyBPjc7ftu8MifW/oz+qk9sb75mD+UONPUaeP1UdWLoGPdXxrreR2XzpfWZz4UJVYlbSPWZ96FXn0e89F97dXGOQhAAAIQmPMH7nbm2vMzqOeMMljpEzSXjlXbzPGHr3zs+obaAzSi6QGsjV2KGfaW9rQ0a3S8ltwte1sYtcTt2TvzLvTo2fXMKXdrV37oggAEnkVgp59jzyLrV7OS98pcXrXxax9HvKFWE6tjr8CaOTVnOtYYo/LGuKPjqd5R46hzB72qQcepmku9TZ0L81fO5uLOWot6T9M9iwdxIQCBdxPQz0Idv5vOnOpX+oaVuXK0fhvq3NvpGZdvRsxcod6amjEde2feOrcTm9aH54rB1LM73NvcHVR9Os6dnb22k5bZtRIfAhDYh4D+/NLxPkrblez8ubqS88pc2qXon4tvqNVQaKCecWtMvTB2rGutelrNWWt89s8hUPvw6P3QcU6d7m29t7nYs9aUi45n5S3FPYGdV4PeAW/PbnMnat6NIXqeR+CJP+t3+lzlc+fHHwTNPTYKSMe5s6k1jaFj71zKFIy6TKn4npYd52oYjtC9Kk9Oa6sG7W0Y18bQs1aXF8Oby9Uya223HxzKRcezOFyNG3WeojfUe6Lmq33iPARqCeQ+02tj6L67Ph80r45V58wxnzu/6GYNtV4+Hfc0SGPouDamXh4dxzip+do8o/bN0rHqIq/Kk+Pdq8EazNYYKXOq97Y1bq7OEWuqb0TM3hiqRce9cWee0+dVxzNz98ZWjTrujcs5CDyNQO+zoefu+NyPOfVzVMerepZi0pJfY7Scbd07M1fWUAehKUPRWoS3P8auKVAvS2kc8t1x2b06W3TUsIg5dK+OPS09cxpXxz0xW89oTh3n4oW94b7oGR17MeLZsJZ6FjSOjr24b5tLsduVQ83ny27aT9S8G8Mn6uHz6GNXI48Ul9J86nzq3Mg7pbmvfK6u0Ku1ezm1Jj0zcjw718+/Ni/3BxJHFhNiWaBqcvQHgpc7nLf7dBzPeI2rie/lvDLXokPZ1OTtOVMTV/esyqN57bhHg57R8ah8LXFzOVnbg0DoZ/x8uuNz4wqFqP1KjN6zd+bu0Xya3tYa+Vz6SMzyiCvqJ0rzdn0l35L2lrsxSnd8fmrieXtG1pSr38sT9o/+bF9qqFNFWRAjC/QamIM+a61Gh8emxKLnzJUaQ76Spivxa862aEjxqYmROpvTWBM3d561PQjY3t993/cgUqfiNG6n6a3rwr+7ej7DWnOcuD/V9xKv3LlVnxMpDS19KNVZG0u1hHGKQy6nrqVi1OrSfRo/ro/OE+IWv/IRk6dEqfjcuFRAaT0X21uL8a7GvVp7jQ7VqONcfXGt5owXp3auNf5Vbp6uGg0xr+6t6UOKpcbq1eade8PcjLswg5vq1PGMnE+IqZx0vFuNqk/Hu+nt0aOfWTruifmEM6mfAyU+Ledm3aeUhit9KdXtxdb6wjgXR9fsWNc0tpe/di4VS3PWxivtqzLUUZQVlxJqE9bsifunFfjju7O9msI5r/YSVG+9pr64p2avciudaemFp7/1/ChunpbcnObtYdrKNqfniWun3IUe9voc6bgn5hvOKCcd78ZA9cVx692+s64arVc+/67WltOXW7uat+a89j+e0XnVqeteLv0Z5O25MlejIcZX/Tru1aEadOzFbbmLM3XWaPX018y5X/kIxYSkNUWlxNmzdo/GTJ2vEd+6pzW37g/5VuptrS+1P9WL1H6drzkf70w4exe3nrxWt9Z99zjWs9Odq7kLlltPT+7mHu/wTtx3YFKj4e7nqTW/3d96t2t4zNqzi9YU75y+3NosXho3pTs++7o/jGs+D2xtMUbNOS/f1TnlHGvW+St5chw1bmnvSF2tuXV/7/iTofYuRCm4XhgvRtjjzYfYer6Ur2e9N/fMJvfU0XrGq7uFd8153aO9bskX6wsxe8/FGKXzO/fWY9ra+9H7VVOIX2Ic9uzMeTQj4t1HoPWe2c+Y3rt9R7W7aE3xzunLrfWw1B6O+jzydAZ9V+NbvT311p5J6Y81rNIR9Xp3xdPgzdXWvMO+H3/Bx7e/vaI9ceEy6YXScWyYPR/2pBrsnfdyh7lUjNR+1eDtK8WM+lp0ennumlPdOi7p0v069vjZh0L3l/KF9RjTi106X9svja1jq6OUc/R6rRZv30zd2ksdpzjU9iR1nvk5BFL3Z062uVG1Fh1r9rgef9W7bMelWBp79jildaVOzWXHKX2BS27N45bLY3uo/fRihbmaeJ7OGE/Pax5dt+NajRqzZ6ycbYygI7feky93xmOSYhF16Zlc/JlrQcdf379Xp/jwHeoU5DBv1+I/p/aH7LrH2+vNpZSnGpDa7817+bw5baa3x4tvz2mMuD8178UbNae9aI2bO59iE+pMreXyKx8d587GtVzeGE/36Djuy+XPraV01pxRLSGWzqX0peZTelrnc3chF0v15/aevlbT4xE1lvLk1mffkxH1tcTQ+6VjG0u5xLF3t0dx0pwttXl7VesInS0alW9qrPOhFtXu1RfmtCY7zmnNrameONZcUWdqf0pzar9q0nEqXu/87PgtupSJnlWtXi/smdK6xu8dhzx/f/3Pz+O1pvqnobYF6WVPwUjNW/G6x8bWtVzRKeC5M6m1kFdrtHt7m2XPpWKk5lNaS/PKxasjzrXw9vLmzntr3lyIm9Pcmtfbn5ur5a8a7Vhj6N4R+UMMy09Zak7VFDXovpy2ljXVkzs7S0Mu551rqV6M1lTKk1vXnujY01qzxzu3ci73uW516P1NPWtas45ra8v1IhdD8+k46tZ5HedyxLUejco7xLC5lXOuB6pRa9Cx7m+J7em25zWX7ldmubPKJO7NscnVVtMn1a/xenNrnJpx1GIZ5vKr9tTY46B7c31Kabdxo5mOe2tM9TdPWOoCpUS0zMfYqeK9WNoAHbfEivE1RpjXODqOZ3Vex7YGj6+XS+v2Ytq5VFwb24uheUaNA8/SvfE0ezWN0uT1Ic4pGzvWu6F3NnVW52vye2finGVqY6X0peZn8KyJ6fW75type7SXOh5Vl8ZtHes9CbpsDO2bjkt1qJ7S/pHrXm1efLuvVq/GTp3zWEYNqTOqUZmnxuGc6tKxxtaxatKx7rfjmCtXc+58XNP6tAZvHOfCr/afa/LZeBo7nFcGukf1xpy6r0aLt0fzW01ebrt/lAZPV8uc6swxj3tVux17TCIXzRV1luZtPXZvyPv17/99KPePb+W/FO/rl7//9/s71OF0LCAE9/65BNSeS+21YBRg6kyY92L3xkrl8eIpi3i2ptmpPGFea495UhpqYtmzns5cjJlrnq7afMpZudXGyWmwd13jhbXSWa9nGkf36Djs9+Y0ThyHvR6L1HwqTszrxcqdya15vEbGz+W+c62lf1d0lvLk1uP98PZ4ffN0pnrpxfTOr5yL9dqcXp36OWP3a71enRrT+0yJnyUar6QtxcvTnIudihO0hzdywUR45z2GMZbWHee9OF5+77xlF5m1xrX6arTkdJT64/XBq9XOpTTV3K1U7BodmjfX21Se0nwNy9yelCbvjKdF70zco/ORRS5ufC5CDGWnuV1DnQpeCuZdBE3oxS7F1RhxPDJWiBmbaJvp5bB69AJ7DdImhvNacylPioGNFf6VRPhdVE0fLEPV4uVSJqkzdp/m6KkxslOGqfyedjvXo6EU01tXfZGL/dU7p3MaR9c93ronjlN7LRPLOxWndl7j1p6btS9Vv81X2lNaD7Fq9oyosZTHW9ee2D12rUaf3k3vvO5JxfW0pva2zGu92us4tjpLdXjrLZq8nHq+JUfUnqtV4+vY/utsa6pVR6qfui/ET+21ub1zysfbUxM7xMkxCWseO82vrHJxPa3e+VQO73xOoxdb2WjMVDw9F+v05r28Oqfsr/C2sTVubT9CHcoinPV45Oa0zjgOsT8Y6lTCVAALOic0BcMWk8pRmvfgls54614crybvrF64eHHir/GMxvOa5sXXOXsu/rP9QIymWnVpnJSe0r64rvFLDOP9ir/aPFqTZah6wjj8rrHmX8F4Z8Oc1p7ad2Xe669lVsMrl987n9qf2pvjoP1Nxc7Nxz7m9qxYS9Vvc5f2lNZX1HElh9dr7bG3J5czd17XUnFmcfVqUU2p+1nS5MXW+rznP7VH58P4ag6tNcbUee+7ofqvu60+PZ/T6u3VWr069VypHxqzpMmLp3OpuxFz6bpXh6crznk1hjnVkarFu18a02rNafE0eTpqYtg9kZEXy5uria/nPO6WQ9TgsVBeujey132qM2r4Fjfqr3rAG2thdo8nwCvci1uai3Fs8SFfa3wbx+YM8zZ2rZ6wLxXTixH3KqtcL3RviKsfiGFs91kusTY7F798H9eUhafd1qr/HMc2R5xL9UlZpPaFOK1/8jalf+S89iWOPQYxr56xenJrkaXdr3lya3ZvLk8uZi27XPzaGLlavBiqOzW286k9MX5p3dOxci7oU42av6YXNXtS/Qj57fmoJ6dL11rHWmNu7NXmzYUYcV7rUX1h3YsR5ywTb19Yr+Gkf0BKe6Cxw1i1pvKoedZximkqnmqzNaZipZjH/V4/crG0dtXkjW2OVG02rsfcxg2HahmNAAAgAElEQVTruT02lpcv9jCs5eIoI+WSY6F641jP6Djs8+Y0d+4e9vRUc+o45LescuyUqT1r43r7bJ1279fw91Br0rghBvJE24BagBYV92qcklBtThhrjJirBYDG8Zoec+uajm2tnja77uWNeZS57o37YjytV99Qa1x73v6z/aC2H6ReTzWGxyLu8e5Obn/uXEpvnG99U+31ydPm1aAMwlhZ5c7ZveGsp0VZeDn1rMbVM3pfvHVlUIoZY4TYtXs1b8tYWXk5U3XqvI6VZ9SlvbV6vfyt9XgxUjxz8y26PI42trLR/V6u3B6PpY2h+eJaLma8qx6/cN47m9pb0zMvntaVYlgT3+7RO5f6jK6NG1npfp233xXVvSPHIW/4eaXGPdcfy7ZWS65nKd5h3junrKJWT5d3PsTVGLV1xH3xvOVkc5X4lfLpea2jJ2/MqbFUi+bW9dJ5y8jLWRM/7qllajXGM9mvfKgIr6jcnlQDrBA9ryA90bk9Hljdn6rDmw9ng8awFn/VeHFP+LUUQ89qzBpmGiPmjx9SGtPbH+a8tx45Ux21pWqMebSnkV1YL51NabW5vR8wmjMVJ85bHbm6cmshViqvV6fu9faobj2j68o2tT/msvXoP3t6UvFyHFVj7djWomc8bcrf2xP023kdxzyWhc5ZLZqjxEfriGPVNGLe5op1evq0hpRGy0H12l61xFM9LWdTOjVmiaXGsbV4azrnja2GkTXpv30MudWIenpKc6nnoHRuxLr3+R3jer3Uu9eioXRPvXwhfmsPY5zWcy21hL2le2Z1aG1RW05rTXyvRs1l6/L2e3VrjNi7nvPeGa0tx8HqU12e9jD3U2+ucE9ULpHut+OUKD2TEhvmUzH0TNiXi6tx4jj8mltLxYxnfgL9EUP/Z+MrP40Zxzqv5zRH2B8/aHNn7Tn9YE6Z6XDG1qi5S9osE4+Pntc9ej6+kbZ6a2uOubyelPJq3bpf60iNo9bcec2VGscYuXujbOxe+8+qR8eqIRVX99WMc/rD+ZQWq0H3hLGn0dsXNcY13ZOqQeOHfd6cPa/rqdpz82FN49gacmdra7PclYuN0RJP+bSe9fqQ4uDt1bkUJ9Wp53Sci6N7wzhXt13Tf/uW+4z28qTmrjBLxayZ15c4OtYYqlPHYb83F+MElnE9d4dtXi9erl9WQ2mfzRP2pvZ78zrnjYP21F2M+bz6lHtqHM56eVP7w7zu9/bqnlQN3tnaHPGsjX2FhWr5yddOalF2bBPHxngBa+PFfZqzVGDcr+ds3ppm1MQp1WcbpDm9+PYyeus2nub26svt8dYss/jP8QNazbV3PszZGPphmOuJjWdrV02Ro8fH3g39QVOb29ORqlXnVauu6zhXQ65OG6f0PIS9useOc2uq17ufuseOlXkce2/VcnFqNFptOQ2ax9NYyqdncvnCXo+35ijFCOuaV8daW2qsuePY/pqKHebtmo2VOpPbo2c8Lak6aue9HCkGNmZuj66VtEQNqiV1TuPrebvuvUAIccOe2nwpHT3zmlPHuZj6c0bHKS4xpubS+xS5xP26Hs7nuGn+EEd749Vnddn9qtc7WzuX06Z1a/067ollY9gavVhakzIMY51TjRojjvWsx9ibSzGyeWzsWj32fPI71FaQB8wTHC9q6sLaeW+PzePFt8I9TbkGpNZ03osbtXhrel6b4zXFq601di5vbk1zxz60sA/xve9sez3NabFrWr+nM+4Pa7o/rOmZ2txx34yYuRo9zSkNJba5/nkxPTa9/Kw27154ueJc2F/6qx5z+q1mb599dkfcGy9Hrr4c01RP7XxrvpyWmrVWRp4+rdnbU6OlZk9tLt0XY1ttusdb07maHoZcJQY2t+awHFJrGj/2MdVPG9Pbo3Maz44ty/g97NT5uDf1fW3LQevWNa3Z7q/5Z41n64j/bPek8qXiqAbvvHIKZ1I5vTxeTM3r5fD2ePXrPjvWuJ6+3HmbL56tqSecq82l8VKaw74U96izJufPN9S6Uce1UMI+FRbPxsK0QF1Pja0GL0bUrNp1XKpF9+u4dN7TqXo9FlZ/b05loHF0HLR6c6Ua9Q1kGHs1ad2puN6+MOfFDDHCvOrWcSpXbt6L4WnLxUiteXG8fDoXxikONlc8p+fDHp2ze3PnUrXofIzh3Qvda8dxf/g1pcPjZmPY9VSdlkFuTyqundfzufpKa14svdu6R8elHLF2PRfG3pxlFWPrvhwPb68315LH0+rVnborqfxWQ9ij571zcU5/VT01Zz2+qkHHXl7d48VVfTq2cT0WMYf9VWuM/8Yy/Kqawtjm0H/DqHp0rLls73Rv7Vg1Kjebs7S3lNPWb/dq3Jgzzpf67XFRLZojVWdun+YJY91fGquuMI5nwq963tvfMufFC3MeU61P85TW4/6fb6j1sB17ojR4aY+3HnPYWHaf5lCNpb0WXDirY42n49T+XC0aI1dbSX/UrDFT4xyvVC0xlldTLl44Z82T/ms71ZiL5eX2zqd4lWrTWDVjT1OuBo2pmnrixRg9Z1VPHKuu1L4r87VvqNV8h5ypH7Ieg6hR+2L36pqtK8eiJoanKeRLzdcwzeVVvXas/2xzRQZxj8YJe7253LxXSyqGx1z3euN4TnsY9qbWomZPn8bRPam4qXk9nxunYmjdMUZqf20Oj4+NaePY+6HnUme8897z7P1c0D701KocanTqmTBWLd6eMFeKX4qjNebi6WeIxtZYNfq8unJxvf12TjXqmmrSXHF/jkNOQyqePZOKHc+G9VQc76zuTZ0vGupeOJ6AFCRbZNijZ1PnUkXpfu8S6p6Wsc0bY1vNKV1xvqZhUY+3V7XW8tJzdtzDKHyIeh+aXp6cxtoaU1y9fLVzqZg9PELO1LnUfI3OK2dr4o/eE+5FyhzbXLXmW7mqXv38sOu5e6dxUnm8GHpn9fmPsb2zV/J6Z+OcagrzIX9qPhfrrrWc1prnIHc+VVPqTGo+FcebL8UI694dSc17OeJciY9qqXlu9IzNr7q959k7X5M3V6e3FnnV5lPtlmHqWfbypuJ4PYlz8ZmMv+q85snx0r2lWGFdNXvMUnF75zVnjOPdWa+XlllJQ6oej6PqSp213DzNv7mXxNlAmtye1TUVFtd1n1dkjSarK7dfdeg4dza1FpsbG5+qSc+nGIR9KV21sTVX6zinLaUvmCZPnzeX06P7dW9Om8fNm9OYtiZvfy6nF8tjlIubiuHNe1q82N7ZO+ZqzHTQFffV7A8MvHti53RdxzUsSmcC9xJ7r1+l3JpXx6XzYV3PhHFKa2q+Jk/tnp4cXg0hn8bScdSUOp/TnDqTms/F0jWNYddjDV4t2jtvTypXKmeYj2t2j+7XNV0PeeOcrcF7nvWsnrM15Na01jBWfh4jrcVq15gaL7e3tBZja/1Rd5gP+eK6/uqd11hhnJtLrWvdGkPXdXx1v+2TV7fO6Vj11IxjDL0jOs7VFvbqfh1XvaGuERwviu5NCbSXSUWFGKlzGr80HhU7p7dXq9VWimHze6xL50uccus1OlVf1NiqS/uVO+/p8uZsbVGn5gl7crlyfDR+HMd4JU2l2Mr2arxSvhPWlYn24EovU3fDmw95r+RS1rm6dK+O41mN4ekeqVl1hPHVO6o1tMb0zns67VzqTGq+FC+ue/y9s7YnqTMz+1aqU9c9jSl99mz8Z3u+9bPSy+0xDXMpTdp7Pa/9sBprYoZ4NTq9WMpateqZ0v4Uh1RfQnxPu9enFDdPk9dzPT9irNq1lzGHcozzet7yy9WQNNQRhgclV3Au2YxzuZgeHA9sCqp3PpWvFEPP9XLKaWrVoJq8celizcpZiuvp8vR7/Q77wvzVHnj5wpx9ZjydpdpsXNV4NV5KM/MfCdRwDydaejmTserVXPFO2rsZ9uhYz/WMZ97RGXp7amw94/Unxcmbt/niZ9eMu1fLN6cxp8s7F/fn1mz93r5UP3Ja9IzXI93TMp6lUzXkeparSdfiWJml4nvn9Wypb7n9WmfrWPnbXKmaNEcrk68hgP7BRBXyc88PA1L7v1qxGq/3nMbxxhpba0zVp/u82Kmz3t4w58WsiZGroeZ8Sk9p3uqNeby5UpzR654Gby7HXJmO1pjLXcqVuiepGkvxnrw+o48ac1fuqXtS6vfMembG1rq0T7o+c9yS29ub6p03H+oIn7+z2LbEzelL8U6diXWFX2s16L4wVjY2bkqTzsc4Ol8a67mUHi9O6We3xrYxlENqLc57P797NMUzOW1eXKtXNXn7r8ytzPW7FhXsiQh7UvN6PoxLF8Q7c+VcKp6dt5q8Wrw5L26IY2P11KpndOzljfqszniuVJsXr2VO8ygrHYfY3pzNWVqv0ae6whlvzs7HuKl9NXlb93j9ralfz6lmXW/V9ZT93rMxojbL13vuRuS4K4bePx1f1aV3NRVP8+o4dS7O1/S+NWYpZ0tuG8t7XnVOuem61qLjWu26T+Po2O5PralWzZFat/Nav8aIY90Xxp6uMOfNl+LW1Bv3xPjer6maU/l1XmOGdW/OzscYmtvjrPm0ptR6KkfrftVYOt+yrrF13BKrdu+UN9S1ye/cpw+ZB1v3BL3evt46QvyaeLU67L6auEF3rQatMZcrtxZzxni1OjV/77i33t589lyJSyrHnZpTmnaZr302ruhdkeOKvnC25271nLmq057X/Dou5arpS2vMUs64XpO7Nlbsn/dZ6OWxcb0zLXljbj1TinuFazgb4sdfNXdqXLM/x6tUk+bVWPa8rUHP6Tie03hhX06Tt9/GjgzjXCpWjluqjlQsra13nNPUGzN1bmUu11AHYRb0bLgpELPn44XN1WcvdW7fKq0hj6fDe/i8fVbn1dq8i1rSUVofxdHTNir2lTi76rpS091nr97jGv0rctToyO3puVs9Z3Iaate8zwE9W/r8CvtzffFy1MRUHTqeFVfzxLFX4+i+eTlSeqyuETxLecJ6iz6vPzFHrd5UjHBetZTGtr7Yt9r+pXTYelKxUvMp3jFXLaNUnNnzrXXN1mPjJw31ShG752ppYMve1rpLse3DV3oo7N6oo3SmVm9JR2m9Nk9q3+z4Xt5Sb7wzzI0jsIL/ihzjiOwfSZ9THddWkOtLb8xUbhsv7hn1uZnKGeZzNebOtaytyNGiJ+7tZa7navuk52yfvbVosm38kSy9nEFTrh57Jrevpx93ninVNZJ7T52XDPXd4msLXqWz1OxavVf2tdQ6U29JR2m9l4GtKcaY/YEyk2MvB85dIxB7OvvuXFO55+mWZ1v36nhEhaNj8ryP6EpbDP1cb3kue/qfy5dba6vq4+6Sztp7p/pClhZeV2qYebZUVy2fmRp/GurwP/2bPlJJY9N3EJ/SaOdX6Sw1u0brHXtiP1fmnp1zVc8Ds1P7vrLfp+XSnt79w2j28zKyPyufvZG6W2Od1JPW2mbv72Gnz2TQ2PJc9uZM5VA9qX21LEvPTWld87Tu1/O7jZW39r+0vqKe79+/f20y1J7oKPTqhZpRsKd3ps6nXeLZPZndi5nxLRv6PuOm3BPT+8zQD++Vyk66Wx67Vc/gyp6Qq5/Alfvce7b3XK7Kq3c9nI/PRilWaV11zqhXc6wc19Z/d93BUH+rBeMVFc/u+qGpunRcW3vtvhh/dp5aPbvt0zuk45F6V/aAvo/s3L2xVt6bUqX6fOi4dH7Gek6DstPxDD3EXEsg1/+SEj2r49L5ns9ZzaHjUk5v3YvRctfj+firnm0dW42qTcdePbPmRuUu8Yj6e+7H6Nq731AH8QGYFjta4Ih4p+gcUWttjLuY2IfshLtTy5N9zyLg/TC4477u9LzUarnrs+VZN3C/amr7n1M+IkYuvl2L93BGzt6Y9lzUGj9XSs9Nbc7afbUce/bN0FDi06Nz5Jnmr3yE5LsXNRLQU2PNuOwtrLhDLbTYexcB+5xEDXeZ6jvyWu67sLjrLrw978j+r/j8159xM3L2xlRtNXerlX+vthotpT2tWkO8O/WW6qld/2Cow6HaP5hYm4B9+xHouez7VVGn6AkPaV2l9+96Ius3PSs1N8jyWG3wn3i/apjvtOfO/rdwOOG5tfe59m6fwj/0qkVry169B7Xs9NzocTDTIWb1d6hHC7grnvew3aXljrz6g1DHd2iakTP2+e39nsFWYz6VtT4bOlYOTx/H+ldzeOr9Ou2+3NX/Vk56P3XcGm/G/qip527vWI8yqr0r+vNZxxrXjnvY5eKNWPv9HeoQ7ClvqANo79LZZnnrI4CeEiPF6BT9OZ3eQ/n2fud4XVl7A+snPytXer/i7Bvu1wqOb8yx+3Pbcref7F16amtht+LuP/YNdQStwEvjFdB3yvFkg6m16XhlH/Tercy9Ipey1fEKDbNzPLGm2cxGxVf2Oh6VhzjPI7D7XVF9Oo4d0Z8hOh7Zuauxe87HulP1e/XpXh17Z1bMffjKR3TZKxLPyKHNtGMFruMZeoh5H4Geh3S02nj/9F6OznN3vB1Y382A/PMIcL/msSXyeAItn/c1d1u9io5TFbToCDFafl55sVvOq+bamuy5GnaaZ8bY+uYPX/kIyU7/2kfORMdL09O8GY14UszAHa7/dtT7wIHPHjeeu7pHH1ABgacRKPmPK/W2fG616mj5eeXFbjl/hcGOZx9tqAPwlou3Y4NO0+Q9YKfVMEMvXGZQvRaTnlzjx2kIQMAncIep9LxOr46az0YvtkfjTS+PrKH+9Ld8nP61j9DcNzXTu8wr5/QB0/FKLbvl2uVfSe3G5S49ejd1fJcu8kIAAucTUN+h49EVxs8v/RzTvDpO6aj5eaWxwtibS+UI86o3t7e0NjJWKZe3rn7501c+wqHTv/bhFc7cPAL2UuvDNS8rkSHQToC72s6MExDYmUB4pnf6ubNCj2cklcFMHV5sb07vzcjP35GxVGftGENdS+pl+2oehhySq+dzsVmDwEgC3NWRNNfHon/rme+acQdTdRebXWtPPZ9Wb2SmvwmoZTkyVm1Ob1+VoQ4H41vqFBwv+GlzT66tpRe7PpgtNbAXAhB4PgE+q57f49oKdzFVtXpn7NvNw5Sez9J6C6ORsVryxr1qpsP8p+9Q28BRsHdxewTsdOa02mb1QOPqeKeeoQUCEHgvAf1s0vF7ybyzcn27qeM3UNmpZn0edRz6EfWO0D0y1qi7kjTUCkPHowTcEUdr0fEdmnI5o74ZOvVi6zinizUIQAACqwjoZ5OOV+kgzz4EekzVjJ+j+xC5T4k+jzqOyrz53p5orN44o6i5fyjxd/C///d3DsIoEXfEseC1KXfoSeX0LsgMvSHPjLipupiHAAQg0EOAz6oeapwJBLyf+9ynsXejlafXkx5Fo+LU5Pa+7hHOZQ11+B51K5waMbvsOaW2lRdll96gAwIQ2I/AKZ+Z+5FD0d0E7M9RT8uuL5Se/Mx5Penpw6g43r3w5roNtReMufUEnvxQradJRghAoJUAv7FvJcb+HIE7fqZ5xstq7DFzuRqvrr3hmRtV46g4NT3rMtQhMH8ndQ1e9kAAAhB4LgHPiOxmPp5L/3mVrTQ/Ss8a+Tt1qC4dv+mZsz1RDi3jnjitZ1JmOujMfuUjbMBQt7STvRCAAASeSWBn8/FM4s+sajej2GqoVnaFZ24u7R6+lwx1KAdTPbepRIcABCBwAoGdzccJ/ND4i0CPkXkrO565OZ3v+Y1dzkwHlcm/Nm9OCedF9aCfVwWKIQABCFwnwNc8rjMkwti/j/jpPEc8c/iYz7dEuer484nyTPErHzHESW+pR/2Ojt9Fly8QOyAAAQhAAAIQ2JMAPibfl1q/WHo7HbI8zlCPujze7+hG/A4m31pWIQABCEAAAhCAwHUC+JjrDGOEGkNd/ZWPmmDjpPdF0suj45aoap513BKLvRCAAAQgAAEIQMASuOJRakiqb9FxTQz2/PgPAn3/XvXyudpQnwBVL4uOW2uI56/Gac3LfghA4BeB2T9w4AwBCEDgDgLxs232Zxw+Zl13q1y3lXPCd6nDBcUEr7tEZILADAL2Bw3P8wzCxIQABO4g4JloPuP+7cROHq727XRQ/0hDfccDQk4IQGAcAX7gjGNJJAhAYD8CvDDwe7Ibl6mGOiA44S213ypmIQCBUwjs9sF6Cjd0QgACZxDY6U3sDsR2e5HSYqYDv0d9h1ovhNcc3cMYAhDYkwDf/duzL6janwA/+/bvUVDI1zw+9kl56Hj3rnYZ6lbXPgtC7kMjruX2zNJFXAhAYAyB0z5Qx1RNFAj0E+BnXz87Ts4jUOvF7nyRYjX2+Nzm71BH3Hd/7cMWrj90vcbpnnnXhsgQgAAEIACB9QT42beeORnLBHJ+rXx6zQ7V2GOou95Qh/J6ko3Coh8aOlbzrONROogDAQhAAAIQ2IWA/qzT8S460fEeAurPSuM7yIzS1P2GOhZ915tq/d2E14Swhw8UjwxzEIAABCDwVAL87HtqZ8+sK+XXUvN3VPnBVH/9T5c37jpki73LUAcNfGjcce3ICQEIQAACEIAABMoEUqb5g4H9J8zdL0B/auo006GEy4Y6BLnTVJfbyQ4IQAACEIAABCAAgZUESqY5ZbZXarS5rn6Vufs71HcVTF4IQAACEIAABCAAgb0J6Bvn1Fjn964qrW7IG+oQfte31HwtJN18ViAAAQhAAAIQgMBMAif4sKtvpwO/YYY6BNvNVO/2rxNmXlhiQwACEIAABCAAAQi0ERhhpkPGx37lQ7+7o+M23OyGAAQgAAEIQAACEICAT2CooR7l8n2pbbP6nRwdt0WbsxuTP4crUSEAAQhAAAIQgECJwEjfOtRQB+EjxZVAlNajid7ZTGOqS11kHQIQgAAEIAABCHwkcNU/jfarQ79DbUvd7fvUO11E7xLsaPp3YoYWCEAAAhCAAAQgEAhc/TNyo8100DT8DTWtLhNQ86zjcgR2QAACEIAABCAAgfcR0JeSOr6LyLQ31KEg3lLn2xouAWY6z4hVCEAAAhCAAAQgYAlceUM94+100DbVUIcEmGoeAghAAAIQgAAEIACBkQR6XkrOMtOhrumGOiTBVI+8QsSCAAQgAAEIQAACEGghMNNMBx18h7qlG+yFAAQgAAEIQAACEICAEFhiqGf/riDUtMuX0rlhEIAABCAAAQhAAAL7EFjhQ5d85SMinfXVjytfTt+n3SiBAAQgAAEIQAACEBhJYIWZDnqXGuqQcLSp9t5M8zdnjLyKxIIABCAAAQhAAALnEVhlpgOZ5YY6JJ1pqjHT5114FEMAAhCAAAQgAIGRBFaa6aB7yXeoRwLyYu38nxj39DIHAQhAAAIQgAAEIPAcAre8oQ74Rr+lfk5LqAQCEIAABCAAAQhAoJfA6rfTQedthjokx1T3XhXOQQACEIAABCAAAQgogTvMdNBwq6EOAjDVehUYQwACEIAABCAAAQi0ErjLTAedtxvqIAJT3Xpl2A8BCEAAAhCAAAQgEAncaaaDhi3+UOLdELiOEIAABCAAgTcT8P4K2jfzoPazCOzgI7d4Qx3bxpvqsy4waiEAAQhA4HwC/MfRzu/hmyvYwUwH/lsZ6iAIU/3mx4LaIQABCEBgJQHvzTT/PYeVHSDXFQK7mOlQw3aGOojCVF+5XpyFAAQgAAEI1BPgDXU9K3buQ2AnMx2obGmogzBM9T6XFiUQgAAEIPBsAsFU82b62T1+UnW7menAdltDHcRhqp90/akFAhCAAAQgAAEIXCOwo5kOFW1tqINATPW1i8dpCEAAAhCAAAQg8AQCu5rpwHZ7Qx1EYqqf8BhQAwQgAAEIQAACEOgjsLOZDhUdYaiDUEx13wXkFAQgAAEIQAACEDiZwO5mOrA9xlAHsZjqkx8HtEMAAhCAAAQgAIE2AieY6VDRUYY6CMZUt11EdkMAAhCAAAQgAIETCZxipgPb4wx1EI2pPvGxQDMEIAABCEAAAhCoI3CSmQ4VHWmog3BMdd2FZBcEIAABCEAAAhA4icBpZjqwPdZQB/GY6pMeD7RCAAIQgAAEIACBPIETzXSo6GhDHQrAVOcvJqsQgAAEIAABCEDgBAKnmunA9nhDHYrAVJ/wmKARAhCAAAQgAAEI+ARONtOhokcY6tgajLV/SZmFAAQgAAEIQAACOxI43UhHpo8y1KEoTPWOjwuaIAABCEAAAhCAwEcCTzHToarHGepQFKaaRxYCEIAABCAAAQjsS+BJZjpQfqShDoVhqvd9iFAGAQhAAAIQgMB7CTzNTIdOPtZQh+Iw1e99WKkcAhCAAAQgAIH9CDzRTAfKjzbU8RphrPd7oFAEAQhAAAIQgMB7CDzVSMcOvsJQh2Ix1e95aKkUAhCAAAQgAIF9CDzdTAfSrzHUoVhM9T4PF0ogAAEIQAACEHg+gTeY6dDFVxnqeG0x1s9/gKkQAhCAAAQgAIH7CLzFSEfCrzTUoXhM9X0PGZkhAAEIQAACEHgugbeZ6dDJ1xrqUDym+rkPM5VBAAIQgAAEILCewBvNdKD8akMdrxnGev0DR0YIQAACEIAABJ5D4K1GOnYQQ/0PCUz1cx5qKoEABCAAAQhAYB2Bt5vpQBpDbe4bpnrdw0cmCEAAAhCAAATOJ4CZ/tVDDLVzlzHWDhSmIAABCEAAAhCAwD8EMNIfrwKGOvFoYKoTYJiGAAQgAAEIQODVBDDTn9uPof7M5MMMxroAiGUIQAACEIAABF5BACOdbjOGOs3m9wqmugISWyAAAQhAAAIQeCwBzHS+tRjqPJ8PqxjrBlhshQAEIAABCEDgeAIY6boWYqjrOP3ehaluBMZ2CEAAAhCAAASOJICZrm8bhrqe1YedGOtOcByDAAQgAAEIQGBrAhjp9vZgqNuZ/T6Bqb4Aj6MQgAAEIAABCGxHADPd1xIMdR+3D6cw1gMgEgICEIAABCAAgdsIYKSvocdQX+P34TTGeiBMQkEAAhCAAAQgMJ0ARnoMYgz1GI6/o21+/qsAAAhYSURBVGCqBwMlHAQgAAEIQAACUwhgpsdhxVCPY/khEsZ6EljCQgACEIAABCBwiQBG+hI+9zCG2sUybhJjPY4lkSAAAQhAAAIQ6CeAke5nVzqJoS4RGrCOqR4AkRAQgAAEIAABCHQTwEx3o6s6iKGuwjRmE8Z6DEeiQAACEIAABCBQRwAjXcfp6i4M9VWCHecx1h3QOAIBCEAAAhCAQDUBjHQ1qiEbMdRDMPYFwVj3ceMUBCAAAQhAAAI+AYy0z2X2LIZ6NuGK+BjrCkhsgQAEIAABCEAgSQAjnUSzZAFDvQRzOQmmusyIHRCAAAQgAAEIfCaAmf7MZPUMhno18UI+jHUBEMsQgAAEIAABCPwkgJHe5yJgqPfpxQclGOtNG4MsCEAAAhCAwM0EMNI3N8BJj6F2oOw0hbHeqRtogQAEIAABCNxHACN9H/tSZgx1idAm6xjrTRqBDAhAAAIQgMBiAhjpxcA70mGoO6DdeQRjfSd9ckMAAhCAAATWEcBIr2N9NdO3qwE4v5ZAeLh4wNYyJxsEIACBtxH4/tefbyt5q3r5Wb9VO6rE8Ia6CtO+m3hjvW9vUAYBCEDgRALWTH/7478nlnCsZl6YHdu6Lxjqc3v3QTnG+iGNpAwIQAACNxLw3kxjquc3BCM9n/HsDBjq2YQXx8dYLwZOOghAAAIPI8Ab6nUNxUivYz07E4Z6NuEb42Oub4RPaghAAAIHEwimmjfTcxqIiZ7D9e6oGOq7O7AgP8Z6AWRSQAACEIAABDIEMNIZOA9YwlA/oIm1JWCsa0mxDwIQgAAEIDCGAEZ6DMfdo2Cod+/QBH0Y6wlQCQkBCEAAAhAwBDDS77oOGOp39ftTtZjrT0iYgAAEIAABCHQRwER3YXvEIQz1I9p4vQiM9XWGRIAABCAAgXcSwEi/s++2agw1d+ADAYw1FwICEIAABCBQRwAjXcfpDbsw1G/ocmeNmOtOcByDAAQgAIHHEsBEP7a1lwrDUF/C947DGOt39JkqIQABCEAgTQAjnWbDyhf+0+NcgjYCmOs2XuyGAAQgAIFzCWCiz+3dauW8oV5N/CH5MNYPaSRlQAACEIDAJwIY6U9ImCgQwFAXALFcJoC5LjNiBwQgAAEI7E0AE713f3ZXh6HevUOH6cNcH9Yw5EIAAhB4MQFM9IubP7h0DPVgoIT7RQBjzU2AAAQgAIFdCWCkd+3Mubow1Of27hjlmOtjWoVQCEAAAo8lgIl+bGu3KAxDvUUb3iMCc/2eXlMpBCAAgbsJYKLv7sB78mOo39Pr7SrFXG/XEgRBAAIQOJ4AJvr4Fh5ZAIb6yLY9TzTm+nk9pSIIQAACqwhgoleRXp/n7y9ffvzfr//9MK3dvnVUnBSBbmGpgMxD4CoBzPVVgpyHAAQg8HwCmOhn93iGAY4xrxjzFHUMdYoM81sQwFxv0QZEQAACENiCACZ6izZMFTHDSEfBGOqprSP4KQQw16d0Cp0QgAAExhHARI9juXOkmUY61p0y1DZ32BvfYOu8XVOWvKFWIoyPIIC5PqJNiIQABCDQRQAT3YXt2EOecb1STOorHZ6h1jk7zq2pvv/TCcYQOIGAfthisE/oGhohAAEI+AT0M93fxSwE9iWAod63NyhrIGA/jDHXDeDYCgEIQOAmApjom8BvmDa8Uda31Km3zBvK/ykJQ71rZ9DVTUA/pDHY3Sg5CAEIQGAYAf1sHhaYQI8goN9bDgb7JFPNd6gfcQ0popYA5rqWFPsgAAEIXCeAib7O8K0R7Bvrq8Za335HpmriU/O6z9ODoX7rTaXunwQw2FwECEAAAuMIYKDHsSTSLwIjjfVMpt9mBic2BHYnED787f/vrhd9zyXw/a8/n1sclT2WAJ+hj23tNoWFt8Hx/7cR5QjhDbUDhSkIRAK8weYurCBgzfS3P/67IiU5INBFgDfQXdg49AICGOoXNJkSxxHAYI9jSaRfBLw305hqbscuBDDQu3QCHbsTwFDv3iH0bU0Ag711e44RxxvqY1r1eKEY6Me3eNsC7/iutP6HW67AwVBfocdZCAgBDDZXopdAMNW8me6lx7leAhjoXnKcSxHQv1HD+xsxSmdbzqRihfmSYS6t52LrGoZaiTCGwEACGOyBMAkFAQhcJoCBvoyQABkCalB1nDn6c6l1/+p4uXwY6hwd1iAwgQAmewJUQkIAAp8IYJ4/IWFiMoGSIS69vU6dz51Lrel8KN2++bbr+kZcz+q6hxFD7VFhDgKLCWCyFwMnHQQeRgDz/LCGHlpOyhB75Xh7vTk9a/fo/tI4FyuutcaI5/hPjytdxhC4gYD3wxCTfUMjSAmBAwh4nxcHyEbiywjYt7zxDa+da8FRe67mTXJL3rC3NiaGupUs+yGwiID3QxOTvQg+aSCwCQHvc2ATaciAQJZAMKLWCKfe/GaD/FjsPVeKO3odQz2aKPEgMJGA98MVkz0ROKEhsJCA93wvTE8qCBxPQM23LSi3liu89hzfoc5RZA0CBxPAaB/cPKQ/mgDG+dHtfX1x3tczcl/50DfZFmDpXNir+fQrGnY9Fy/ESq1rTK/JGGqPCnMQeDABjPaDm0tpWxHAOG/VDsRAYCoBDPVUvASHwFkEMNtn9Qu19xPANN/fAxRAYAcCGOoduoAGCBxAALN9QJOQOIUApnkKVoJC4FEEMNSPaifFQOAeApjte7iTdRwBTPM4lkSCwBsJYKjf2HVqhsANBDDdN0An5U8CmGUuAgQgMJsAhno2YeJDAALVBDDd1ajY+A8BzDJXAQIQ2IEAhnqHLqABAhBoJoD5bkZ2zAFM8jGtQigEIPAPAQw1VwECEHgNAUz4+lZjjtczJyMEILCewP8D8lwgKT4dTL8AAAAASUVORK5CYII="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 🌪 Exercise\n",
    "Reproduce the following plot from the Soubiran et al. 2018 CDS table: http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/619/A155#/browse\n",
    "\n",
    "![image.png](attachment:5298dde0-215c-4d62-9561-c1d07b0a28ef.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}