{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Accessing DC2 data in PostgreSQL at NERSC part 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Owner: **Joanne Bogart [@jrbogart](https://github.com/LSSTDESC/DC2-analysis/issues/new?body=@jrbogart)**  \n",
    "Last Verified to Run: **2020-08-03**\n",
    "\n",
    "This notebook introduces some additional features of the PostgreSQL database at NERSC.   \n",
    "\n",
    "__Learning objectives__:\n",
    "\n",
    "After going through this notebook, you should be able to:\n",
    "  1. Discover which object catalogs are available\n",
    "  2. Query on native quantities in those catalogs\n",
    "  3. Make use of custom functions, in particular for area searches\n",
    "\n",
    "__Logistics__: This notebook is intended to be run through the JupyterHub NERSC interface available here: https://jupyter-dev.nersc.gov. To setup your NERSC environment, please follow the instructions available here: https://confluence.slac.stanford.edu/display/LSSTDESC/Using+Jupyter-dev+at+NERSC\n",
    "\n",
    "### Prerequisites\n",
    "Please see the first notebook in this series for instructions on how to gain access to the database.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "import psycopg2\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "%matplotlib inline \n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make the db connection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "dbname = 'desc_dc2_drp'\n",
    "dbuser = 'desc_dc2_drp_user'\n",
    "dbhost = 'nerscdb03.nersc.gov'\n",
    "dbconfig = {'dbname' : dbname, 'user' : dbuser, 'host' : dbhost}\n",
    "dbconn = psycopg2.connect(**dbconfig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finding Data\n",
    "Tables for the Run1.2i data as well as a view to make dpdd quantities more easily accessible are in the `schema` (acts like a namespace) `run12i`.  To reference, say, a table called `position` for Run1.2i use `run12i.position`. \n",
    "\n",
    "### Finding Datasets\n",
    "To find out which datasets are available and by what schema names, query the table `run_provenance`. It's in a special schema known as `public` which does not normally need to be specified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "schema_name    run_desig  sim_prog  db_ingest       remarks\n",
      "run12p_v4      1.2        phosim    obj+fs_complete Fixes problems with previous Run1.2p versions\n",
      "run21i_dr1b_v1_pg2 2.1        imsim     obj_complete    Representative run2.1i 6-month dataset; exclude most non-dpdd\n",
      "run21i_dr1b_v1 2.1        imsim     obj_complete    Representative run2.1i 6-month dataset\n",
      "run12i         1.2        imsim     obj_complete    Output from production Imsim run 1.2\n",
      "run12p_dpdd    1.2        phosim    superseded      Flaws include 1. Problems with calibrations 2. Wrong pixel scale. DROPPED 2020-07-07\n",
      "star_truth     2.2i                 star_truth_complete Source: /global/cfs/cdirs/lsst/shared/DC2-prod/Run2.2i/truth/startruth\n",
      "sne_truth      2.2i                 sn_truth_complete Source: /global/cfs/cdirs/lsst/shared/DC2-prod/Run2.2i/truth/sntruth\n",
      "minion_test    2.2i                 obs_db_complete Source is minion16 sqlite db\n",
      "run22i_dr6c_object 2.2i_dr6c  ImSim     obj_complete    Source is GCR-generated dpdd-only parquet\n",
      "run21i_v1      2.1        imsim     obj_complete    First run2.1i dataset DROPPED 2020-07-07\n",
      "run211i_v1     2.1.1      imsim     obj_complete    None\n"
     ]
    }
   ],
   "source": [
    "cols = ['schema_name', 'run_designation','simulation_program', 'db_ingest', 'remarks']\n",
    "hdrs = ['schema_name', 'run_desig', 'sim_prog','db_ingest', 'remarks']\n",
    "# Additional columns in run_provenance store software and input data versions\n",
    "prov_query = 'SELECT '  + ','.join(cols) + ' from run_provenance'\n",
    "with dbconn.cursor() as cursor:\n",
    "    cursor.execute(prov_query)\n",
    "    fmt = '{0!s:14} {1!s:10} {2!s:9} {3!s:15} {4!s}'\n",
    "    print(fmt.format(hdrs[0], hdrs[1], hdrs[2], hdrs[3], hdrs[4]))\n",
    "    for record in cursor:\n",
    "        print(fmt.format(record[0], record[1], record[2], record[3], record[4]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Normally only datasets where the `db_ingest` field contains 'complete' are of interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pick one of the supported datasets\n",
    "schema = 'run12p_v4'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Querying on Native Quantities\n",
    "Unlike DPDD quantities (all in a single view), native quantities are split across several tables. The first notebook in the PostgreSQL collection shows how to find out which tables belong to a schema and what columns a table has. Alternatively, if you know the column names you want, you can query for the table name.  The following looks for the table containing `ext_shapeHSM_HsmShapeRegauss_resolution`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "select table_name from information_schema.columns where table_schema='run12p_v4' and column_name='ext_shapehsm_hsmshaperegauss_flag'\n",
      "dpdd_ref\n"
     ]
    }
   ],
   "source": [
    "column = 'ext_shapehsm_hsmshaperegauss_flag'\n",
    "find_table_query = \"select table_name from information_schema.columns where table_schema='{}' and column_name='{}'\"\n",
    "find_table_query = find_table_query.format(schema, column)\n",
    "print(find_table_query)\n",
    "with dbconn.cursor() as cursor:\n",
    "    cursor.execute(find_table_query)\n",
    "    for record in cursor:\n",
    "        print(record[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is some necessary fussiness here:\n",
    "* Note `ext_shapeHSM_HsmShapeRegauss_flag` has been transformed to all lower-case in the query. This is required when querying information_schema, where this string is a __value__ in the database (not a column name). \n",
    "* In the query single quotes are used around literals like `run12p_v4`. Double quotes won't work."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now suppose we wanted to combine a cut on this quantity, in table dpdd_ref, with cuts on DPDD quantities like `clean`.  Then the query has to be made on a join of these two tables, where we specify that the value of dpdd_ref.object_id = dpdd.objectid. This causes the corresponding rows from each table (or view) to be treated as if they were assembled into one long row. Here is a simple query showing how this is done. A more realistic one would have more conditions in the `where` clause and might join more than two tables.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "select count(object_id) from run12i.dpdd join run12i.dpdd_ref on run12i.dpdd_ref.object_id = run12i.dpdd.objectid where (ext_shapeHSM_HSmShapeRegauss_flag = 'f') and clean\n",
      "CPU times: user 2.38 ms, sys: 1.38 ms, total: 3.76 ms\n",
      "Wall time: 21.9 s\n",
      "1526574\n"
     ]
    }
   ],
   "source": [
    "schema = 'run12i'\n",
    "join_query = 'select count(object_id) from {schema}.dpdd join {schema}.dpdd_ref '\n",
    "join_query += 'on {schema}.dpdd_ref.object_id = {schema}.dpdd.objectid'\n",
    "join_query = join_query.format(**locals())\n",
    "where = \" where (ext_shapeHSM_HSmShapeRegauss_flag = 'f') and clean\"\n",
    "join_query += where\n",
    "print(join_query)                # confirm the query looks reasonable\n",
    "with dbconn.cursor() as cursor:\n",
    "    %time cursor.execute(join_query)\n",
    "    for record in cursor:\n",
    "        print(record[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adjustments for Larger Datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above query I switched to a different, significantly smaller dataset (2.8 million objects rather than 13.7 million).  For the much larger 2.1i_dr1b_v1 (over 78 million objects) we have to pay some attention to performance.\n",
    "The dpdd view is formed from a join of 3 tables. **dpdd_ref** is one of them, so in the above query that table is joined to itself, increasing the resources needed to make the query unnecessarily.  It's more efficient to just get all quantities from a join of tables only, but the dpdd view is more than just the result of a join with some of the columns renamed.  `clean` is formed by doing logical operations on several native-quantity flags. Starting with run 2.1i_dr1b_v1 it can be expressed as ```good and not deblend_skipped``` where both `good` and `deblend_skipped` are columns in **dpdd_ref**. (In earlier runs, `good` existed only in the dpdd view as the result of logical operations on several native-quantity flags).  We also have to exclude non-primary objects, as the dpdd view does.  The flag `detect_isprimary` is in the **position** table. The query for run 2.1i_dr1b_v1 can be written as shown in the next cell. Skip it if you're in a hurry; even with these techniques it still takes about 13 minutes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "select count(position.object_id) from run21i_dr1b_v1.position left join run21i_dr1b_v1.dpdd_ref on run21i_dr1b_v1.position.object_id = run21i_dr1b_v1.dpdd_ref.object_id where detect_isprimary and good and (not deblend_skipped) and (ext_shapeHSM_HSmShapeRegauss_flag = 'f')\n",
      "CPU times: user 15 ms, sys: 25.8 ms, total: 40.8 ms\n",
      "Wall time: 17min 16s\n",
      "52412800\n"
     ]
    }
   ],
   "source": [
    "schema = 'run21i_dr1b_v1'\n",
    "join_query = 'select count(position.object_id) from {schema}.position left join {schema}.dpdd_ref '\n",
    "join_query += 'on {schema}.position.object_id = {schema}.dpdd_ref.object_id'\n",
    "join_query = join_query.format(**locals())\n",
    "where = \" where detect_isprimary and good and (not deblend_skipped) and (ext_shapeHSM_HSmShapeRegauss_flag = 'f')\"\n",
    "join_query += where\n",
    "print(join_query)                # confirm the query looks reasonable\n",
    "with dbconn.cursor() as cursor:\n",
    "    %time cursor.execute(join_query)\n",
    "    for record in cursor:\n",
    "        print(record[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In all such joins, **position** should be the first table in the join and left joins should be used, as in the example. This will in general result in better performance since it forces the join(s) to be done in the order specified. All tables are indexed on `object_id`. Only the **position** table has additional indexes (on `detect_isprimary` and on `coord`, a special column used in area searches)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## User-defined Functions (UDFs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many math functions from the c library have been wrapped and incorporated in an extension module installed in the database. They have their normal c library names with the prefix c_. Functions with a floating point argument or return value usually have two versions, such as c_log (double precision natural logarithm) and c_logf (single precision). They can be incorporated in queries as in this example using the command-line interface program psql:\n",
    "```\n",
    "desc_dc2_drp=> select c_asin(1.0);\n",
    "\n",
    "     c_asin\n",
    "-----------------\n",
    " 1.5707963267949\n",
    "```\n",
    "There are also functions specially crafted for HSC or LSST catalogs with suggestive names like `patch_contains`, `tract_from_object_id` (used in q3 in the first notebook of this series), `sky_to_pixel`,..\n",
    "```\n",
    "desc_dc2_drp=> select count(*) from run12i.dpdd where tractsearch(objectId, 5063);\n",
    " count\n",
    "--------\n",
    " 233982\n",
    "(1 row)\n",
    "```\n",
    "### Restricting by tract or patch\n",
    "\n",
    "Let's try the last query from the previous section restriced to tract 3446 (picked at random), which has about 570,000 objects. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "select count(position.object_id) from run21i_dr1b_v1.position left join run21i_dr1b_v1.dpdd_ref on run21i_dr1b_v1.position.object_id = run21i_dr1b_v1.dpdd_ref.object_id where detect_isprimary and tractsearch(position.object_id, 3446) and  good and (not deblend_skipped) and (ext_shapeHSM_HSmShapeRegauss_flag = 'f')\n",
      "CPU times: user 2.09 ms, sys: 423 µs, total: 2.51 ms\n",
      "Wall time: 7.33 s\n",
      "381378\n"
     ]
    }
   ],
   "source": [
    "schema = 'run21i_dr1b_v1'\n",
    "join_query = 'select count(position.object_id) from {schema}.position left join {schema}.dpdd_ref '\n",
    "join_query += 'on {schema}.position.object_id = {schema}.dpdd_ref.object_id'\n",
    "join_query = join_query.format(**locals())\n",
    "where = \" where detect_isprimary and tractsearch(position.object_id, 3446) and  \"\n",
    "where += \"good and (not deblend_skipped) and (ext_shapeHSM_HSmShapeRegauss_flag = 'f')\"\n",
    "join_query += where\n",
    "print(join_query)                # confirm the query looks reasonable\n",
    "with dbconn.cursor() as cursor:\n",
    "    %time cursor.execute(join_query)\n",
    "    for record in cursor:\n",
    "        print(record[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Tract 3446 has about 1/136 of the objects in the run21i_dr1b_v1 dataset. The elapsed time for the tract query is less than 1/136 of the time taken by the original query. It pays to restrict queries to a tract when possible, even if it means issuing the same query for several tracts. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Area Searches\n",
    "The  **dpdd** view has one extra column, `coord`, which is not formally a DPDD quantity. `coord` is an alternate way (other than `ra` and `dec`) to express location.  A `coord` value is a triple of doubles representing a position on a sphere in units of arcseconds. This column is indexed, which can make certain calculations faster. In particular, using the functions `conesearch` and `boxsearch` (which take a `coord` as input) rather than starting with `ra` and `dec` makes queries much faster.  There are also functions to translate between `coord` and `(ra, dec)`.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Cone search\n",
    "Find all stars satisfying quality cuts within a fixed radius of a particular coordinate.  The function `coneSearch` returns true if `coord` is within the cone centered at (ra, dec) of the specified radius, measured in arcseconds. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SELECT ra, dec, mag_i from run21i_dr1b_v1.dpdd  where (magerr_i < 0.04) and clean and (extendedness < 1.0) and coneSearch(coord, 54.5, -31.4, 240.0)\n",
      "CPU times: user 429 µs, sys: 1.98 ms, total: 2.41 ms\n",
      "Wall time: 4.01 s\n",
      "67 objects found \n"
     ]
    }
   ],
   "source": [
    "schema = 'run21i_dr1b_v1'\n",
    "band = 'i'\n",
    "mag_col = 'mag_' + band\n",
    "min_SNR = 25   \n",
    "max_err = 1/min_SNR\n",
    "pop = 'Stars'\n",
    "ra = 54.5\n",
    "decl = -31.4\n",
    "radius = 240.0\n",
    "where = ' where (magerr_{band} < {max_err}) and clean and (extendedness < 1.0) and coneSearch(coord, {ra}, {decl}, {radius})'\n",
    "qcone = ('SELECT ra, dec, mag_{band} from {schema}.dpdd ' + where).format(**locals())\n",
    "print(qcone)\n",
    "with dbconn.cursor() as cursor:\n",
    "    %time cursor.execute(qcone)\n",
    "    records = cursor.fetchall()\n",
    "    nObj = len(records)\n",
    "    print('{} objects found '.format(nObj))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgkAAAH3CAYAAAAxEswMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de5zkdX3n+9enIZYOEY0wdEJkpm0DnVVkjWkvJLrijWRzAT1uLruzG3XXTIga3WwI0dNHNPH0LtETSTzGyywb9cT25OKJ6wZviKGNJm1Mg4CojJp2ZgS0HRQlgDRqf84fv19DTfPr7qruqvrV5fV8POoxXb/6/X71/XbXVL3re/tFZiJJkrTeWN0FkCRJ/cmQIEmSKhkSJElSJUOCJEmqZEiQJEmVDAmSJKmSIUGSahAR8xFxU93lkDZjSJDaFBE/GhF/GhFfjIi7I+JoRCxGxKUR8UNN+z0uIl4dERP1lXZzEfGgiPjNiFiIiG9GxD0RcSQi5iLiqXWXT1K9jq+7ANIgiYgnAVcB3wTeBiwBu4GzgF8F3gt8pdz9ccCrgHngUI+LuqWI+GHgA8CZwOXAnwO3AxPAc4D5iPjJzPz72gopqVaGBKk9rwS+Bzw+M49pKo6I7weO60UhIuJBmfntHRw/BrwbOAP46cy8Yt3jrwJeAHx3RwUdIRFxHHBcZt5Td1mkTrG7QWrPI4Evrg8IAJl5R2Z+CyAiXg389/KhqyIiy9vzy8efEhHviogvlV0Wt0bEuyPijPXnLY97Z0T8bNmtcTfwO+Vjj42Iv46I5fI8R8rz7NmiHs8GngT8/vqAUNYlM/NPMvOTTeV4aES8ISJuioiViPiniPi9iGisK+/byzKfUv58W0T8c0T8eUT8QEX9HhER/09Zh5WIOBgRv10GmU2Vz/HWiDhcHrscEX8TEU9ft9/uiPjjiPhy2aVyKCL+W0XZnxcRH4qIW5q6Xt4QEQ9et985ZR1/NSIujIh/AlaAnygfP6H83dxY/l2+FhFXRMSTK+qwJyL+Z/k7+kZEvGV9uaS62JIgtedLwNMj4imZ+bFN9vsr4OHAfwL+K/C5cvta0/0vAqcC7wBuBvYC+4GPRcSZmXl03fl+HPgZ4M3AW4EvR8Ru4EqKro/XA7eW5zy3fO4jm5TvOeW/b9u0tqXyQ+sjwI8B/wP4FPCvKFpWfgz4+YrD3kfRzfK/A/8CeDFwD/Afms77I8ACcAfwfwNHgXOA11J0e7x4i6L9Zfn8fwx8EXgY8ESKrp6/KZ/jJOATwIOBAxS/lx8Hfht4DPBzTed7aXmePwJuK8/za+V+T6t4/pdRvI8eAO4CvhIRD6Lokno88BfAG4EHUgSIpwIfbzr+QRS/14+W5Tm7fL6jFL9bqV6Z6c2btxZvFB+M3wESuAb4Q2AfcErFvi8s9zun4rETKradQfFt9OXrtmd5O3vd9vPL7Y/fRj2uBr7Vxv4vKp/rt9Ztv7Tc/nNN295ebvvDdfu+gaL74sSmbe+nGNdx4rp9Xw+sAqdvUqaHlM/z21uU/U3AN4DT1m1/aXn8M7f4uzyv3O9JTdvOKbctAw9Zt///UT724opzRdPP8+V+/3ndPv8L+Frdr3Vv3jLT7gapHZn5txTf9v4ceATFN8l3AjeXzdLf1+J57lz7OSK+v/y2+w3gIMU30PU+lZkL67Z9q/z35yPiAe3VhBMpBim26jzgTopv7M1e2/T4em9ad/8qijEbe6HovgB+Gvj/gAdExMlrN+BDQABPZ2PfpmiZOKf8/d1PRATwS8AVwLfXPceHy92esbb/2t8lIsYi4iHlfn9bPlz1d5nLsoupyS9StFa8ef3Ombn+srurFC1Dza4Cdq/v4pDqYEiQ2pSZi5n5yxRN2/8CeAlwGPgN4BWtnCMifigi3hER3wD+maKr4ChFs/b9+u0pujnW+yhFc/Yrga9HxPsj4qXlB9tWbqcICq2aAL6UmXc3b8zMr1B0dzyi4pjD6+7fVv77sPLfMyiCwIUUdW++fbDc55SNCpTFAMGLKLpXvhoRf1+OA/jRpt12l8/3SxXP8dn1zxERT4iID1MEom+W+y2VD7f6dzkduCEzVzcqe5Ov5f0HoK7/PUm1cUyCtE3lt8IbgRsj4s+Af6Jomv69zY4rB+RdAfwwRXP9Zyj65Fcpui+qwvv9ZjKUz/9LEfFaivEKz6Ropr84Ip6ZmdduUozPAo+LiInMPLRZeZufcoPtscFj39tkf7ivnm+haE2osrTB9qJAmX8UEe+haMl4BvCbwCsiYn9mvq3pOd7D/Vs21twCEMV6FldRjKO4qPz32xStHx+kxb/LWtE2K3eTjX5HcN/vSaqNIUHqgMz8ejnC/dHNmzfY/TEUaxO8IDPf3vxARDyMolWhnee+mmKMwWsi4qzy54uAf7fJYf8T+PfAfwQubuFpDgE/GREPbG5NiIgfpBgbcKidMpf+iXK8RWZeuY3joTj4CMXgwDeWsycWgFmKQZlHKVpNHtjCc5wP7AJ+tjk4RcRUm0X6AnBmRIy12Jog9S27G6Q2RMQzo5gPv377JEXXw+eaNt9R/ru+mXrt2+Mx//8i4nnAD9GiiPiBss+92ecovt1u1VT9HuAfgYsi4twNzv+CiHhCefevge8Hfn3dbr/d9HhbspjB8RHgeeUsh/XPf+JmUwEjYlc5k6D5nLdRBJYfKO9/j6JL5qc2mH74wKa+/8q/C/Dy1mp0r7/gvtkq65/P1gENFFsSpPb8IfDQiHgvcAPFaP0zKLoZHkAxsn3NIsU35VeUg/S+DfwDZRcF8AcRsZeiufuJFGsXbNq8vs7zgN8om9u/SPH/+ZcppvrNbXZgZq5GxP9GMUDwgxHx1xRTBv8Z2FOW5V9SzvsHLqNodfiDss//WuApwL+lWK3x8jbK3ezXgb8DPhURl1GEnIdQtLQ8t/z30AbHnkGxBsW7KbpP7qCYffJTwJ807ffycvvfRMTbKaZvPhCYAn6hfJ55ii6Fu4H3R8RbKULDzwOVgyI38QcUU0zfHMXS1h+neG38RPnc/7XN80n1qXt6hTdvg3Sj+AB6K0VAuI1iOuTNFKsXnl2x/69TNKt/lyIwPL/c/kiKqW5rAxc/DDyW4sNqft05Enhnxbl/jGJmxZcoAsjXgY8Bz2mjPruA36JYR+BbFLMFjpTnffK6fR9KsZbBzeV+S8BrgMa6/d5elvn4ddvPoWJKKMXYjDeXz3sPxbTCj1G0Ujxwk7KfRDGt8gaKLoU7gE9TDIR8QEXZX0vRFbBC0aXzSYplsx/WtN8zKILcnRRdFW+nGPyYwKsr6vLCDcr2YOCS8m9/D/A1ihDyk037zAM3VRz7/PLcE3W/3r15i8xWx9dIkqRR4pgESZJUyZAgSZIqGRIkSVIlQ4IkSapkSJAkSZUMCZIkqZIhQZIkVTIkSJKkSi7LDJx88sk5MTFRdzEkSeqZq6+++tbM3L3ZPoYEYGJigsXFxbqLIUlSz0TE4a32sbtBkiRVMiRIkqRKhgRJklTJkCBJkioZEiRJUiVDgiRJqmRIkCRJlQwJkiSpkiFBkiRVMiRIkqRKhgRJklTJkCBJkioZEiRJUiVDgiRJqmRIkCRJlQwJkrpmeXmOhYUJ5ufHWFiYYHl5ru4iSWrD8XUXQNJwWl6e4+DB/ayu3gXAysphDh7cD8D4+L46iyapRbYkSOqKpaWZewPCmtXVu1hamqmpRO2zJUSjzpYESV2xsnKkre39xpYQyZYESV3SaOxpa3u/GYaWEGmnDAmSumJycpaxsV3HbBsb28Xk5GxNJWrPoLeESJ1gSJDUFePj+5iaOkCjsRcIGo29TE0dGJim+kFvCZE6wTEJkrpmfHzfwISC9SYnZ48ZkwCD1RIidYItCZJUYdBbQqROsCVBkjYwyC0hUifU0pIQEa+JiOsj4tqIuCIiTi23nxQRV0XEHRHxxnaPb3p8T3mOC7tdF0mShlVd3Q2vy8yzMvOxwOXAxeX2u4FXAlt9uG90/JpLgQ90ssCSJI2aWrobMvP2prsnAFluvxP4eET8yHaOB4iIZwNLwJ0dK7AkSSOotjEJETEL/ArwLeBpnTg+Ik4Afgd4Flu3RkiSpE10rbshIq6MiBsqbucDZOZMZp4GzAEvaff8Gxz/u8ClmXlHC+XbHxGLEbF49OjRdp9ekqShF5m59V7dLEDEXuB9mXlm07bnA9OZuWV4aD4+Ij4GnFY+9FBgFbg4MzccBAkwPT2di4uL262CJEkDJyKuzszpzfappbshIk7PzC+Ud88DbuzE8Zn5lKZ9Xg3csVVAkCRJ1eoak3BJRExRfNM/DFyw9kBEHAJOBB5QDkI8NzM/GxGXAW/JzMXNjpckSZ1Re3dDP7C7QZI0alrpbnBZZkmSVMmQIEmSKhkSJElSJUOCJEmqZEiQJEmVDAmSJKmSIUGSJFUyJEiSpEqGBEmSVMmQIEmSKhkSJElSJUOCJHXJ8vIcCwsTzM+PsbAwwfLyXN1FktpS11UgJWmoLS/PcfDgflZX7wJgZeUwBw/uB2B8fF+dRZNaZkuCJHXB0tLMvQFhzerqXSwtzdRUIql9hgRJ6oKVlSNtbZf6kSFBkrqg0djT1napHxkSJKkLJidnGRvbdcy2sbFdTE7O1lQiqX2GBEnqgvHxfUxNHaDR2AsEjcZepqYOOGhRA8XZDRpqy8tzLC3NsLJyhEZjD5OTs75Jq2fGx/f5etNAMyRoaDkFTZJ2xu4GDS2noEnSzhgSNLScgiZJO2NI0NByCpok7YwhQUPLKWiStDOGBA0tp6BJ0s44u0FDzSlokrR9tiRIkqRKhgRJklTJkCBJkioZEiRJUiVDgiTVbHl5joWFCebnx1hYmGB5ea7uIkmAsxskqVZeY0T9zJYESaqR1xhRPzMkSFKNvMaI+pkhQZJq5DVG1M8MCZJUI68xon5mSJCkGnmNEfUzZzdIUs28xoj6lS0JkiSpkiFBkiRVMiRIkqRKhgRJklTJkCBJkioZEiRJUiVDgiRJqmRIkEaMlyWW1CoXU5JGiJclltQOWxKkEeJliSW1w5AgjRAvSyypHYYEaYR4WWJJ7TAkSCPEyxJLaochQRohXpa4Nc4AkQq1zG6IiNcA5wOrwNeA52fmLRFxEvBu4PHA2zPzJe0cXz52FvBW4MTy8cdn5t1drpI0MLws8eacASLdp66WhNdl5lmZ+VjgcuDicvvdwCuBC7dzfEQcD7wTuCAzHw2cA3ynC+WXNKScASLdp5aQkJm3N909Achy+52Z+XGKsND28cC5wPWZeV2539cz83sdK7ikoecMEOk+tS2mFBGzwK8A3wKe1qHjzwAyIj4E7Ab+LDNf25kSa73l5TmWlmZYWTlCo7GHyclZm2M18BqNPaysHK7cLo2arrUkRMSVEXFDxe18gMycyczTgDmgcuzBZjY4/njgycC+8t/nRMQzNijf/ohYjIjFo0ePbqOGo22t37Z4M817+20d4KVB5wwQ6T5dCwmZ+czMPLPi9t51u74LeO4Onqr5+JuAj2bmrZl5F/B+4HEblO9AZk5n5vTu3bt38PSjyX5bDStngPQXZ5rUq67ZDadn5hfKu+cBN3bo+A8BF0XELuAe4KnApR0ostax31bDzBkg/cGZJvWra3bDJWXXw/UUgw1ftvZARBwCXg88PyJuiohHldsvi4jpzY7PzNvKY/8RuBa4JjPf16tKjRJX7pPUbbZY1q+WloTM3LB7ITMnNtj+whaPfyfFNEh10eTk7DEJH+y3ldRZtljWzxUXtS3220rqNlss61fbFEgNvkHtt3XqpjQYbLGsnyFBI8WBUNLgWPs/aaivT2Tm1nsNuenp6VxcXKy7GOqBhYWJDRbK2cvZZx/qfYEkqSYRcXVmTm+2j2MSNFIcCCVJrTMkaKQ4EGo0uSCPtD2GBI0Ul9wdPS4hLm2fIUEjxambo8cFeaTtc3aDRs6gTt3U9jgORdo+WxIkDTXHoUjbZ0iQNNQchyJtnyFB0lBzHIq0fY5JkDT0HIcibY8tCdI6zqmXpIIhQWrinPr7GJYkGRKkJs6pLxiWJIEhQTqGc+oLhiVJYEiQjuGc+oJhSRIYEqRjOKe+YFiSBIYE6RjOqS8YliSB6yRI9+Oceu6t/9LSDCsrR2g09jA5OTvyvxdp1BgSJFUyLEmyu0GSJFUyJEiSpEqGBEmSVMmQIEmSKhkSJElSJUOCJEmqZEiQJEmVDAmSJG1ilC+b7mJKkiRtYO2y6WtXRV27bDowEouN2ZIgSdIGRv2y6YYESZI2MOqXTTckSJK0gVG/bLohQZKkDYz6ZdMNCZIkbWB8fB9TUwdoNPYCQaOxl6mpAyMxaBGc3aABs7w8x9LSDCsrR2g09jA5OTsy/1kl1WOUL5tuSNDAGPWpSJLUa3Y3aGCM+lQkSeo1Q4IGxqhPRZKkXjMkaGCM+lQkSeo1Q4IGxqhPRZKkXjMkaGCM+lQkSeo1ZzdooIzyVCRJ6jVbEiRJUiVDgiRJqmRIkCRJlQwJkiSpkiFBkiRVMiRIkqRKhgRJklTJkCBJkirVEhIi4jURcX1EXBsRV0TEqeX2kyLiqoi4IyLeuI3jvy8i3hERn46Iz0XEK3pVJ0mShk1dLQmvy8yzMvOxwOXAxeX2u4FXAhdu8/hfABqZ+Rjgx4Ffi4iJThdekqRRUEtIyMzbm+6eAGS5/c7M/DhFWGj7+PLfEyLieOBBwD3A7UiSpLbVdu2GiJgFfgX4FvC0Dh3/buB84CvALuA3M/MbHSmwJEkjpmstCRFxZUTcUHE7HyAzZzLzNGAOeEm759/g+CcA3wNOBR4B/FZETG5Qvv0RsRgRi0ePHt1GDSVJGm5dCwmZ+czMPLPi9t51u74LeO4Onqr5+H8HfDAzv5OZXwP+DpjeoHwHMnM6M6d37969g6eXJGk41TW74fSmu+cBN3bo+CPA06NwAvCkds8tSZIKdY1JuCQipoBV4DBwwdoDEXEIOBF4QEQ8Gzg3Mz8bEZcBb8nMxU2O/2PgbcANQABvy8zre1QnSZKGSi0hITM37F7IzIkNtr9wq+Mz8w6KaZCSJGmHXHFRkiRVMiRIkqRKhgRJklTJkCBJkioZEiRJUiVDgiRJqmRIkCRJlQwJkiSpkiFBGgHLy3MsLEwwPz/GwsIEy8tzdRdJ0gAwJEhDbnl5joMH97OychhIVlYOc/DgfoNCDxjONOgMCdKQW1qaYXX1rmO2ra7exdLSTE0lGg2GMw0DQ4I05FZWjrS1XZ1hONMwMCRIQ67R2NPWdnWG4UzDwJAgDbnJyVnGxnYds21sbBeTk7M1lWg0GM40DAwJ0pAbH9/H1NQBGo29QNBo7GVq6gDj4/u6+ryjPmjPcKZhcHzdBVB3LC/PsbQ0w8rKERqNPUxOznb9Q0H9a3x8X09Cwdpr7vjjH8Z3v3s78B2AewftrZVlFKzV0/+HGmSRmXWXoXbT09O5uLhYdzE6Zm1UdfOgqbGxXT359qjRVPWaq9Jo7OXssw/1plCSNhURV2fm9Gb72N0whBxVrV6res1VcdDeaBn1LqdhYEgYQo6qVq+1+tpy0N7ocJ2Izqg7aBkShpCjqtVrrby2HLQ3Wlpt0az7Q7Cf9UPQMiQMIUdVq9eqXnMRD+C4406ilzMq1D9aadHshw/BftYPXcfObhhCjqpWr/ma03qNxp7yw//+29ds9iHoa6c/uo4NCUOqF1PepGa+5tRscnK2cpZVc4tmP3wI9rNWgla32d0wxOzrk1SXVhbxcvzU5vqh69iWhCG1ft76KC5mI6leW7UutdLaMEqqFsGbmjpQazeeIWFI2dcnqd85luU+G32xm5o6UOsCZIaEIWVfn6RB4FiWQr9+sXNMwpCyr0+SBke/frEzJAypfhjwIklqTb9+sTMkDKm6Lg8sSWpfv36xc0zCELOvrz95GW9J6/XrIE5DgtRDTk2VtJF+/GJnd4PUQ/2wFrsktcqQIPVQv45glqQqhgSph/p1BLMkVTEkDACvwTA8+nUEsyRVceBin3Og23Dp1xHMklQlMrPuMtRueno6FxcX6y5GpYWFiQ0uFbq31vW8JUmDLSKuzszpzfaxu6HPOdBNklQXQ0Kfc6CbOsWxLZLaZUjocw50UyesjW0puq7y3rEtBgVJmzEk9DmvwaBOcBEnSdvh7IYB0I9LdWr76rh2g2NbJG2HLQlSD9XV7O/YFknbYUiQeqiuZn/HtkjaDkOC1KJOzA6oq9nfsS2StsMxCVILOrXyZaOxZ4PFsbrf7O/YFkntsiVhBDlfvn2d6iYY9mZ/X1vScLElYcR4LYjt6VQ3wTBfu8HXljR8DAkjZrNvxL6Rb6yT3QTD2uzva0saPnY3jBjny2/PsHcTdIKvLWn4GBJGjPPlt8fZAVvztSUNH7sbBkAnV+ibnJw9pt8Y/EbcqmHtJugUX1vS8GmpJSEinhQRD266/+CIeOJ2nzQiXhMR10fEtRFxRUScWm4/KSKuiog7IuKNLZznwojIiDi5adsrIuKLEXEwIn5qu2XsF51eoc9vxOoWX1vS8InM3HqniE8Bj8ty54gYAxYz83HbetKIEzPz9vLnlwKPyswLIuIE4MeAM4EzM/Mlm5zjNOAy4EeBH8/MWyPiUcD/CzwBOBW4EjgjM7+3WXmmp6dzcXFxO1XpuoWFiQ0GzO3l7LMP9b5AkqShEBFXZ+b0Zvu0OiYhsilNZOYqO+iqWAsIpROALLffmZkfB+5u4TSXAhetHVs6H/izzFzJzC8BX6QIDAPLwWCSpLq0GhKWIuKlEfF95e1lwNJOnjgiZiPiy8A+4OI2jz0PuDkzr1v30A8DX266f1O5reoc+yNiMSIWjx492s7T95SDwSRJdWk1JFwA/ARwM8UH7xOB/ZsdEBFXRsQNFbfzATJzJjNPA+aADbsVKs67C5ihOlhExbbK/pTMPJCZ05k5vXv37lafvueceidJqktLXQaZ+TXgl9s5cWY+s8Vd3wW8D3hVi/s/EngEcF1EADwcuCYinkARYE5r2vfhwC0tnrcvDfMKfZKk/tZSSIiIM4A3A+OZeWZEnAWcl5n/53aeNCJOz8wvlHfPA25s9djM/DRwStO5DgHT5cDF/wW8KyJeTzFw8XTgk9spYz9x6p0kqQ6tdjf8d+AVwHcAMvN62mxZWOeSsuvheuBc4GVrD5Qf+q8Hnh8RN5UzFoiIyyJi01GYmfkZ4C+AzwIfBF681cwGSZJUrdUZCrsy85Nl8/6a7273STPzuZs8NrHB9he2sn9mzgJ22EuStEOttiTcGhGPpBwEGBH/BvhK10olSZJq12pLwouBA8CPRsTNwJcopi5KkqQhtWlIiIj/0nT3/cBVFK0PdwLPpRg7IElST3TyWjba2lYtCWvXa5gCHg+8l2Itgv8A/G0XyyVJ0jHWrmWzdhGxtWvZAAaFLtk0JGTm7wJExBUU12745/L+q4G/7HrpJEkqLS3NHHOVUYDV1btYWpoxJHRJqwMX9wD3NN2/B5joeGkkSdqA17LpvVYHLv4p8MmIeA/FDIfnAO/oWqkkSVqn0dizwVVxvZZNt7TUklCuPfAC4Dbgm8ALMvO/dbNgkiQ181o2vdfy5Z4z8xrgmi6WRZKkDXktm95rOSRIklQ3r2XTW60OXJQkSSPGkKCRtbw8x8LCBPPzYywsTLC8PFd3kSSpr9jdoJHkoiyStDVbEjSSNluURZJUMCRoJLkoiyRtzZCgkbTR4isuyiJJ9zEkaCS5KIskbc2QoJE0Pr6PqakDNBp7gaDR2MvU1AEHLaqrnFGjQePsBo0sF2VRLzmjRoPIlgRJ6gFn1GgQGRIkqQecUaNBZEjoA/ZTSsPPGTUaRIaEmq31UxbXSM97+ykNCtJwcUaNBpEhoWb2U0qjwRk1GkTObqiZ/ZSdsbw85zXm1TPbfb05o0aDxpBQs0ZjT9nVcP/tao1Ty9RLvt40SuxuqJn9lDtnl416ydebRokhoWb2U+6cXTbqJV9vGiV2N3TYdvoq7afcGbts1Eu+3jRKbEnoIKcz1sMuG/WSrzeNEkNCB9lXWQ+7bNRLvt40Suxu6CD7Kutjl416ydebRoUtCR3ksquSpGFiSOgg+yo77/OffxHz88czPx/Mzx/P5z//orqLJEkjw5DQQfZVdtbnP/8ibrnlzcD3yi3f45Zb3mxQkKQeicysuwy1m56ezsXFxbqLoXXm54/nvoDQ7DjOOee7vS6OJA2ViLg6M6c328eWBPWxqoCw2XZJUicZEtTHjmtzuySpkwwJ6lunnrq/re2SpM5ynQT1rTPOeBMAt9xygKKL4ThOPXX/vdslSd1lSFBfO+OMN/VdKNjO9TkkaRDZ3SC1wetzjKbl5TkWFiaYnx9jYWHCv7dGhiFBaoPX5xg9BsPhZwjcmCFBaoPX5xg9BsPhZgjcnCFBaoPX5xg9BsPhZgjcnCFBaoPX5xg9BsPhZgjcnCFBaoPX5xg9BsPhZgjcnFMgpTaNj+8zFIyQtb+1016H0+TkLAcP7j+my8EQeB9DgiRtwWA4vAyBmzMkSJJGmiFwY45JkNRxzjuXhoMtCZI6am3e+Vof79q8c8Bva9KAqaUlISJeExHXR8S1EXFFRJxabj8pIq6KiDsi4o0tnOfCiMiIOLm8/6yIuDoiPl3++/Ru10XSsZx3Lg2PurobXpeZZ2XmY4HLgYvL7XcDrwQu3OoEEXEa8CygeTLrrcDPZ+ZjgOcBf9rRUkvakvPOpeFRS0jIzNub7p4AZLn9zsz8OEVY2MqlwEVrx5bHfyozbynvfgZ4YEQ0OlNqSa1w3rk0PGobuBgRsxHxZWAf97UktHrsecDNmXndJrs9F/hUZq7soJiS2uTiQ9Lw6FpIiIgrI+KGitv5AJk5k5mnAXPAS9o47y5ghk2CRUQ8Gvh94Nc22Wd/RCxGxOLRo0dbfXpJW3BVSml4RGZuvVc3CxCxF3hfZp7ZtO35wHRm3i88RMRjgI8AayOjHg7cAjwhM78aEQ8H/gZ4QWb+XStlmJ6ezsXFxZ1VRJKkARIRV2fm9Gb71DIFMiJOz8wvlHfPA25s9djM/DRwStO5DlEEilsj4qHA+4BXtBoQJElStbrGJFxSdj1cD5wLvGztgfJD/xWk8xAAABCNSURBVPXA8yPipoh4VLn9sojYNPFQdFv8CPDKcnrltRFxyhbHSJKkCrV3N/QDuxskSaOmle4Gl2WWJEmVDAmSJKmSIUGSJFUyJEiSpEqGBEmSWjCKl0D3UtGSJG1hVC+BbkuCJElbGNVLoBsSJEnawqheAt2QIEnSFkb1EuiGBEmStjCql0A3JEiStIVRvQS6sxskSWrB+Pi+oQ8F69mSIEmSKhkSJElSJUOCJEmqZEiQJEmVDAmSJKmSIUGSJFUyJEiSpEqGBEmSVMmQIEmSKhkSNHKWl+dYWJhgfn6MhYUJlpfn6i6SJPUll2XWSFlenuPgwf33Xhd+ZeUwBw/uBxi55VYlaSu2JGikLC3N3BsQ1qyu3sXS0kxNJZKk/mVI0EhZWTnS1nZJGmWGBI2URmNPW9slaZQZEjRSJidnGRvbdcy2sbFdTE7O1lQiSepfhgSNlPHxfUxNHaDR2AsEjcZepqYOOGhRkio4u0EjZ3x8n6FAklpgS4IkSapkSJAkSZUMCZIkqZIhQZIkVTIkSJKkSoYESZJUyZAgSZIqGRIkSVIlQ4IkaSgsL8+xsDDB/PwYCwsTLC/P1V2kgeeKi5Kkgbe8PMfBg/vvvRT8ysphDh7cD+AKqztgS4K6xlQvqVeWlmbuDQhrVlfvYmlppqYSDQdbEtQVpnpJvbSycqSt7WqNLQnqClO9pF5qNPa0tV2tMSSoK0z1w8kupNb4e+q9yclZxsZ2HbNtbGwXk5OzNZVoOBgS1BWm+uGz1oW0snIYyHu7kPwAPJa/p3qMj+9jauoAjcZeIGg09jI1dcDuzR2KzKy7DLWbnp7OxcXFuosxVNaPSYAi1fufdnAtLEyUH3zHajT2cvbZh3pfoD7l70mDIiKuzszpzfaxJUFdYaofPnYhtcbfk4aJsxvUNePj+wwFQ6TR2LPBN2S7kJr5e9IwsSVBUkscGNYaf08aJoYEbZsjuEeLXUit8fdU8P1hODhwEQcubocDEyVtxPeHweDARXWNiyVJ/avub/G+PwwPBy5qWxzBLfWnflgS3feH4VFLS0JEvCYiro+IayPiiog4tdx+UkRcFRF3RMQbWzjPhRGREXHyuu17ynNc2K06jDoXS5L6Uz98i/f9YXjU1d3wusw8KzMfC1wOXFxuvxt4JbDlh3tEnAY8C6iKppcCH+hQWVXBEdxSf+rGt/h2uy98fxgetYSEzLy96e4JQJbb78zMj1OEha1cCly0duyaiHg2sAR8pjOlVRVHcEv9qdPf4rezzLTvD8OjtjEJETEL/ArwLeBpbR57HnBzZl4XEc3bTwB+h6KFwa6GLnOxJKn/TE7OVs4s2O63+M26Lzb7/+/7w3DoWktCRFwZETdU3M4HyMyZzDwNmANe0sZ5dwEz3NdF0ex3gUsz844WzrM/IhYjYvHo0aOtPr0k9bVOf4t3EOJo61pLQmY+s8Vd3wW8D3hVi/s/EngEsNaK8HDgmoh4AvBE4N9ExGuBhwKrEXF3Zt5vEGRmHgAOQLFOQovPLUl9r5Pf4l1merTV0t0QEadn5hfKu+cBN7Z6bGZ+Gjil6VyHgOnMvBV4StP2VwN3VAUESVJrOt19ocFS15iESyJiClgFDgMXrD1QfuifCDygHIR4bmZ+NiIuA96SmS6NKEk9stYisbQ0w8rKERqNPUxOzjreYES4LDMuyyxJGj0uyyxJkrbNkCBJkioZEiRJUiVDgiRJqmRIkCRJlQwJkiSpkiFBkiRVMiRIkqRKhgRJklTJkCBpZCwvz7GwMMH8/BgLCxMsL8/VXSSpr9V17QZJ6qnl5bljLlS0snKYgwf3A3gdAmkDtiRIGglLSzPHXMkQYHX1LpaWZmoqkdT/DAmSRsLKypG2tksyJEgaEY3Gnra2SzIkSBoRk5OzjI3tOmbb2NguJidnayqR1P8MCZJGwvj4PqamDtBo7AWCRmMvU1MHHLQobcLZDdq25eU5lpZmWFk5QqOxh8nJWd9w1dfGx/f5GpXaYEjQtjidTJKGn90N2hank0nS8DMkaFucTiZJw8+QoG1xOpkkDT9Dgtqytvb9ysphII55zOlkkjRcHLiolq0frAhJERSSRmOvsxskacgYEtSyqsGKawHh7LMP1VEkSVIX2d2gljlYUZJGiyFBLXOwoiSNFkOCWuba95I0WgwJaplr30vSaHHgotri2veSNDpsSZAkSZUMCZKG1triX/PzYywsTLC8PFd3kaSBYkiQWuQHzmBZW/yrWB00771SqX83qXWGBKkFfuAMHq9UKu2cIUFqgR84g8fFv6SdMySor/Rrk74fOIPHxb+knXMKpGqzvDzH0tIMKytHaDT2cNJJP8NXv/qOe7+xrzXpA7VPu2w09pRdDfffrv40OTm77oJkLv4ltcuWBNWiqo//llve0rdN+q422X2dbkVy8S9p52xJUC02uqJklX5o0l/7YGlu+fDS2J2z/jLknWpFcvEvaWcMCapFOx/8/dKk7wdO92w2MNTfuVQfuxtUi40/+OOYezbpjwYHhkr9yZCgWmzUx3/qqRfYhzyCnIkg9Se7G1QL+/jVzJkIUn8yJKg29vFrjaFR6k+GBEl9wdAo9R/HJEiSpEqGBEmSVMmQIEmSKhkSJElSJUOCJEmqZEiQJEmVDAmSJKmSIUGSJFWqJSRExGsi4vqIuDYiroiIU8vtJ0XEVRFxR0S8sYXzXBgRGREnN207KyIWIuIzEfHpiHhgN+siSdKwqqsl4XWZeVZmPha4HLi43H438Ergwq1OEBGnAc8CjjRtOx54J3BBZj4aOAf4TmeLLknSaKglJGTm7U13TwCy3H5nZn6cIixs5VLgorVjS+cC12fmdeX5vp6Z3+tMqSVJGi21jUmIiNmI+DKwj/taElo99jzg5rUw0OQMICPiQxFxTURc1KHiSrVZXp5jYWGC+fkxFhYmWF6eq7tIkkZE10JCRFwZETdU3M4HyMyZzDwNmANe0sZ5dwEzVAeL44EnUwSPJwPPiYhnbHCe/RGxGBGLR48ebbN2Um8sL89x8OB+VlYOA8nKymEOHtxvUJDUE10LCZn5zMw8s+L23nW7vgt4bhunfiTwCOC6iDgEPBy4JiJ+ELgJ+Ghm3pqZdwHvBx63QfkOZOZ0Zk7v3r273epJPbG0NMPq6l3HbFtdvYulpZmaSiRplNQ1u+H0prvnATe2emxmfjozT8nMicycoAgGj8vMrwIfAs6KiF3lIManAp/tYNGlnlpZOdLWdknqpONret5LImIKWAUOAxesPVC2DpwIPCAing2cm5mfjYjLgLdk5uJGJ83M2yLi9cA/UgxofH9mvq+L9ZC6qtHYU3Y13H+7JHVbLSEhMzfsXihbB6q2v7CV/TPznRTTIKWBNzk5y8GD+4/pchgb28Xk5GyNpZI0KlxxUepj4+P7mJo6QKOxFwgajb1MTR1gfHxf3UWTNALq6m6Q1KLx8X2GAkm1sCVBkiRVMiRIkqRKhgRpCLgqo6RucEyCNODWVmVcmwGxtioj4FgGSTtiS4I04FyVUVK3GBKkAeeqjJK6xZAgDbiNVl90VUZJO2VIkAbc5OQsY2O7jtk2CKsyOthS6n8OXJQG3NrgxKWlGVZWjtBo7GFycravBy062FIaDJGZdZehdtPT07m4uOF1oyR12MLCxAYXrtrL2Wcf6n2BpBEUEVdn5vRm+9jdIKnnHGwpDQZDgqSec7ClNBgMCZJ6blAHW0qjxpAgqee8BLY0GJzdIKkWXgJb6n+2JEiSpEqGBEmSVMmQIEmSKhkSJElSJUOCJEmqZEiQJEmVDAmSJKmSIUGSJFUyJEiSpEqGBEmSVMmQIEmSKhkSJElSJUOCJEmqZEiQJEmVDAmSJKlSZGbdZahdRBwFDtddji45Gbi17kJ0gfUaPMNat2GtFwxv3axXYW9m7t5sB0PCkIuIxcycrrscnWa9Bs+w1m1Y6wXDWzfr1Tq7GyRJUiVDgiRJqmRIGH4H6i5Al1ivwTOsdRvWesHw1s16tcgxCZIkqZItCZIkqZIhYYBExKGI+HREXBsRi+seuzAiMiJO3uT44yLiUxFxecVjWx7fLd2qV0T8RkQcjIjPRMRru1X+zXSjbhHx2Ij4xNo5I+IJ3azDBuXadr02OjYiHhYRH46IL5T//kC369FO+crHtlu310XEjRFxfUS8JyIe2u16tFq28rFt1avV47upW/Ua9PePTV6Lbb1/HN+ZqqiHnpaZx8yDjYjTgGcBR7Y49mXA54ATt3l8N3W0XhHxNOB84KzMXImIUzpc3nZ0+m/2WuB3M/MDEfEz5f1zOlfclu2kXvc7Fng58JHMvCQiXl7e/52OlbY9na7bh4FXZOZ3I+L3gVdQT906Xa9heP+oOnZY3j+q/mZtvX/YkjAcLgUuAjYcYBIRDwd+FrhsO8fXZCf1+nXgksxcAcjMr3WrkNu0k7ol94WGhwC3dKOA27ST19L5wDvKn98BPLtTheqQbdctM6/IzO+Wdz8BPLyTBduhnf7/H9j3j00M/PvHJtp6/zAkDJYEroiIqyNiP0BEnAfcnJnXbXHsH1K8qFabN7ZxfDd1vF7AGcBTIuIfIuKjEfH4jpe6Nd2o238GXhcRXwb+L4pvpb22k3rd79jSeGZ+BaD8t65vb92oW7P/CHygc8VtWcfrNQTvHxv9vYbh/WOjurX3/pGZ3gbkBpxa/nsKcB3wr4B/AB5Sbj8EnFxx3M8Bbyp/Pge4vPx5VyvHD1q9yvs3AG8AAngC8CXK2TxDULc3AM8tf/5F4MpBqddGx5b3v7luv9t6Xa9u1a3p8RngPYP0Wtzk2IF+/9jitTjQ7x9b1K2t9w9bEgZIZt5S/vs1ijeapwKPAK6LiEMUTZjXRMQPrjv0J4Hzyn3+DHh6RLwTeGSLx3dVF+oFcBPwV1n4JMW38Z4PqupS3Z4H/FX5819SvIn11A7qVXXsWvmXI+KHAMp/a2ni7VLdiIjnUYS/fVm+Q/dSF+o16O8fm/29Bv39Y7O6tff+0etk5G3bifIE4MFNP/898NPr9jnEFkmedd9K2z1+UOoFXAD8XvnzGcCX6fE3gS7W7XPAOeXPzwCuHpR6bXYs8Drg5eXPLwde28t6dbluPw18Ftjd6zp1s17tvpYHpV6D/v6xRd3aev9wdsPgGAfeExFQzEp5V2Z+cKOdI+JU4LLM/JkelW+7ulWvPwH+JCJuAO4Bnpfl/4oe6lbdfhX4o4g4Hrgb2Kjvu1t2Uq/Njr0E+IuI+E8Uo7Z/oXtV2FC36vZGoAF8uHz8E5l5QddqcX/dqlfdulWvQX//2OzYtt4/XHFRkiRVckyCJEmqZEiQJEmVDAmSJKmSIUGSJFUyJEiSpEqGBEmSVMmQIKkWUfA9SOpjrpMgqWciYoLi4kZXAWcD1wKPAR4EvDszX1Vb4STdjysuSuq1KeAFmfmiiHhYZn4jIo4DPhIRZ2Xm9XUXUFLBpj5JvXY4Mz9R/vyLEXEN8Cng0cCj6iuWpPVsSZDUa3cCRMQjgAuBx2fmbRHxduCBdRZM0rFsSZBUlxMpAsO3ImIc+Nc1l0fSOrYkSKpFZl4XEZ8CPgMsAX9Xc5EkrePsBkmSVMnuBkmSVMmQIEmSKhkSJElSJUOCJEmqZEiQJEmVDAmSJKmSIUGSJFUyJEiSpEr/P3m4wS7tF9r8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cmags = pd.DataFrame(records, columns=['ra', 'dec', mag_col])\n",
    "\n",
    "plt.figure(figsize=(8, 8))\n",
    "plt.xlabel('ra')\n",
    "plt.ylabel('dec')\n",
    "plt.suptitle(pop + ' Cone search', size='xx-large', y=0.92)\n",
    "p = plt.scatter(cmags['ra'], cmags['dec'], color='y')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how fast the query is.   Compare time, # objects found and scatter plot after increasing radius, e.g. by a factor of 10 to 2400.0. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Box search\n",
    "Find all stars, subject to quality cuts, with the specified ra and dec bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SELECT ra, dec, mag_i from run21i_dr1b_v1.dpdd  where (magerr_i < 0.04) and clean and (extendedness < 1.0) and boxSearch(coord, 54.4, 54.8,-31.6, -31.3)\n",
      "CPU times: user 2.23 ms, sys: 0 ns, total: 2.23 ms\n",
      "Wall time: 6.44 s\n",
      "573 objects found \n"
     ]
    }
   ],
   "source": [
    "ra1 = 54.4\n",
    "ra2 = 54.8\n",
    "decl1 = -31.6\n",
    "decl2 = -31.3\n",
    "\n",
    "where = ' where (magerr_{band} < {max_err}) and clean and (extendedness < 1.0) and boxSearch(coord, {ra1}, {ra2},{decl1}, {decl2})'\n",
    "qbox = ('SELECT ra, dec, mag_{band} from {schema}.dpdd ' + where).format(**locals())\n",
    "print(qbox)\n",
    "with dbconn.cursor() as cursor:\n",
    "    %time cursor.execute(qbox)\n",
    "    records = cursor.fetchall()\n",
    "    nObj = len(records)\n",
    "    print('{} objects found '.format(nObj))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgYAAAH3CAYAAADAGZeBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOy9e5xeR3nn+Xu62xK28QW3TAuBLSGitkUcYhJBohCISAhknIkNy+6EHeHIEOI4xEGBOGwYhUvCeJfNBaJZYhyHIehj5GQZcoGAMyFOaMbOKLAyV+M2EggZjFBbbhM7NiC7u2v/OOdtvf32uZ+qU0/V+X0/n/cj9Xkvp6pOXZ56biXGGBBCCCGEAMCY7wIQQgghRA8UDAghhBCyDAUDQgghhCxDwYAQQgghy1AwIIQQQsgyFAwIIYQQsgwFA0IIcYiIzIjIfb7LQUhVKBgQYgkRuVhEbhaRr4jI90TkhIgcFJF3ichThj73QyLyNhHZ5K+02YjIJhExI6/viMghEfkDEXmy7zISQtwy4bsAhMSAiPwogE8A+FcAfwbgCIDzATwLwC8B+DCAb6Uf/yEAbwUwA+Box0WtyocBfCj9/xkAfhjArwN4iYj8oDFmyVvJCCFOoWBAiB3eDGARwHOMMSvUxiLyRADjXRRCRE43xnzXwk/dZYz5wMhvfweJcPCDAD5r4R5BIiLjAMaNMY/5LgshLqApgRA7PAPAV0aFAgAwxjxijHkIAETkbQD+NH3rE0Pq+qvS958vIreIyNdSc8QDIvIhEZke/d30ex8QkZ9NTRbfA/B/pO9dKiJ/KyJz6e98Pf2dC1vU8Xj67+Mj5XiqiLw/vddJEblbRF4vIjL0mY+JyKMictHId28QkSUR+amiG4vIk0XkT0Tk3vQecyLyTyLykyOfO19E/lhEviEij4nIURH5v0Rk7cjndonI34vIsfRzXxeR/yIiZ418bkfazr8kIteJyFcBnATwY+n7Z4rI74rIPWk73y8iHxeRH8+ow4Ui8jci8m8i8qCI3DhaLkI0QI0BIXb4GoCfFJHnG2NuL/jcXwF4GoBfBPB/AphNr//P9N//AGADgH0AvglgI4CrAdwuIpcYY06M/N4PA7gMwHsA/AmAb4jI+QBuQ2LWeCeAB9LffHF6769XqM/pIrJu8H8A2wC8HsCnAXxp8CERmUzLvh7AHyMxofz79L7PAHBt+tFXA/gigA+IyI8ZYx4XkcsA/AqAdxpj/rGkPP8NwLPTe3wFwHkAfgSJWeafhsryLwDOAnBTWs8fBvCbAH4gLdeA16W/sxfAt9Pf+eX0cy/MuP9uJPPlTQC+A+BbInI6EvPRcwB8EMC7ATwBidDwEwDuGPr+6QD+EcAn0/JsT+93Aom2iRA9GGP44ouvli8AL0CykzYAPgPgjwDsBPDkjM++Jv3cjoz3zsy4No1kl/pbI9dN+to+cv2K9PpzGtRj09Dvjr4+CeC8kc//Xvrey4euCRIByAD4gaHrl6fX/jOAdUg0EF8AsLakTOek3/vNks/dAOBBABeMXH9d+v0XlbTzrvRzPzp0bUd6bQ7AOSOf/+30vV/N+C0Z+v9M+rlfH/nMRwDc77vv8sXX6IumBEIsYIz5H0h2gf8vgKcj2WF+AMA3UxX1aRV/59HB/0Xkieku+EEAX0ayMx3ls8aYAyPXHkr//TkRWVOvJsvsB/DT6etyJM6SlwL4axF5wtDnLkdiQvnLoToYAL8/KMPQ9Y8g2XH/FoBbAZwL4D8aY06WlOW7AB4DsCNtj1WkZoufB/BxAN8VkXWDF4B/SD+2bK4YtLOIjInIOenn/kf6dlY77zepOWiI/4BEK/Ge0Q+nbTDMEhKNzjCfAHD+qPmCEN9QMCDEEsaYg8aYVyBRc29Foka/F8CvAXhTld8QkaeIyD4ReRDAvyExA5xAouJ+UsZXvpZx7ZNIVNtvBjAvIreKyOuGTANVOGqMuS19/a0x5neRRFe8AInGY8AmAPdkfP/u9N+nj1x/PZLojOcA2GOMuausICZx8nsjElPIcRH5n6ld/+Khj52PpN1/Hkl7Db8GZVkOtRSR54rIPwB4FInJ5QQSMwhQvZ23IHHSrBKhcb9Z7RT67fTf8yp8n5DOoGBAiGVMwj3GmD9GYgd/CImauhARGUOy4/05AO8C8HIAL0Gya/8SssfrqgiE9P4/j8Qv4PcAnInE5n+PiFzaqFIJg533T4zesuA7o+89G8Agp8MPVr2xMWYvEp+F1yNR678ewBdF5FXpRwZt89c4pekYff0hkORqQLJb34BE4Lg8ff9nRn5rmLxIj6K6D7NY8J4UvEdI59D5kBCHGGPmU0/27x++nPPxHwBwCYBXGWPeP/yGiJyHRHtQ5953ArgTwNtF5Fnp/98I4D/W+Z0hBiGXTxy6dhSJdmSUrUPvAwBSlfnN6bW/BnCdiPytMea/Vbm5MebrSBz83i0iTwJwAMD1SPJGnADwMIAnGGNuK/mpK5DkZvhZY8xw+S7K/UY2hwFcIiJjFbUGhAQBNQaEWEBEXpTGt49e34xkkZwduvxI+u+oynqwq1wxLkVkF07tsquU5UnDoYIps0h2vW3U1pen/w7nMPhbAN8nIi8bur8AuG7o/QH/BcCFAF6JxLTyaQA3isiGopuKyBlpBMAyxphvIxEwnpT+vYjEfPKSnFDBJwzZ8jPbGYnvQx0+iFNRI6P3oxaABAs1BoTY4Y8AnCsiHwZwF4AFJNEEuwCsQeLBPuAgEq3Bm0TkXCQL9qeQ2OrvAfCHIrIRwDEkpoiX4pT9uwq7APyaiPw1kpC8CQCvQBLGt7/ib1wiIq9M/386EhPALyJR4/8/Q597BxInvD8XkUG44s8C+HcA/tgY80UAEJH/BcBVAH7XGPMv6bVXIhEy3i8iL8lw2BswjSTnw4eQ+As8gsTX4SUA3jf0ud9Kr/+TiLw//e0nALgIwP+GxDQzA+C/A/gegFtF5E+QCAo/ByDTsbGAPwTwMgDvEZFBeOIaJOGKn0USjkpIePgOi+CLrxheSBapP0EiFHwbSejiN5GkFd6e8flfAfBVJAKEAXBVev0ZSMLYBs6H/4AkGmAGwMzIbxgAH8j47WcjiYj4GhKhYx7A7QBeVqEem7A6THEhrcv7AWzM+M5TkeRdOIEkrHIWwBuQhuwh0XY8gET4mRj57tXpPV5XUKZJJNqGu5CYCx5BkhPhOgBrRj57LhK/isNpWR5Aopl4K4ZCLZFEKHwKifPhibRu56dledvQ53ak116TU7azkAhHX0USOXE/EsHjeUOfmQFwX8Z3r0p/e5Pv/ssXX8OvwcAlhBBCCKGPASGEEEJOQcGAEEIIIctQMCCEEELIMhQMCCGEELIMBQNCCCGELEPBgBBCCCHLUDAghBBCyDIUDAghhBCyDFMiA1i3bp3ZtGmT72IQQgghnXHnnXc+YIw5f/Q6BQMAmzZtwsGDB30XgxBCCOkMEbk36zpNCYQQQghZhoIBIYQQQpahYEAIIYSQZSgYEEIIIWQZCgaEEEIIWYaCASGEEEKWoWBACCGEkGUoGBBCCCFkGQoGhBBCCFmGggEhhBBClqFgQAghhJBlKBgQQgghZBkKBoQQQghZhoIBIYQQQpbxIhiIyNtF5Asi8jkR+biIbEivT4rIJ0TkERF5d93vp++9SUS+IiJfFpGXdFEfQgghJBZ8aQx+3xjzLGPMpQA+CuAt6fXvAXgzgOuafF9EngngFQC+H8DPALhBRMZdVIDkMze3HwcObMLMzBgOHNiEubn9vosUBWxXQkgXeBEMjDEPD/15JgCTXn/UGHMHEgGh9vcBXAHgL4wxJ40xXwPwFQDPtVZwUsrc3H58+ctX4+TJewEYnDx5L7785aujW8S6XqT70q6EEP9M+LqxiFwP4BcAPATghZa+/1QA/zL0sfvSa1nfvxrA1QBw4YUX1r09yeHIkT1YWvrOimtLS9/BkSN7MDW101Op2jE3tx9HjuzByZNfx9q1F2Jy8jIcP75vuZ6DRRqAszrWadfR8m7efH2wbU8I6R5nGgMRuU1E7sp4XQEAxpg9xpgLAOwHcG3d38/5vmR9NOf7Nxljthljtp1//vl1b09yOHny67Wuaydrp37s2I25i3Sb+xRpIKq2KzULhCTQ9NYcZ4KBMeZFxphLMl4fHvnoLQBe3uJWw9+/D8AFQ+89DcCxFr9NarJ2bbb2Je+6drJ26jmyZmPhp8piXrVdDx/ebV1oIUQjRQu/FgE5VOHEV1TClqE/Lwdwj6XvfwTAK0RkrYg8HcAWAJ9uU1ZSj82br8fY2Bkrro2NnYHNm6/3VKJ21Fnsmwo/RWaCAVXadW5uPxYW5jPvEarGhtgl1IVqlLKFv8qY8l1GzfiKSnhHalb4AoAXA9g9eENEjgJ4J4CrROS+NNIAIvJeEdlW9H1jzJcAfBDA3QD+O4BfNcYsdlUpktjYL7roJqxduxGAYO3ajbjoopuCtXHnL/YrrVZthJ8qZoIq7Vo06YWqsSH2CHmhGqVs4ddg0tQgnDTFi/OhMSbXdGCM2ZRz/TUVv389gDC3p5EwNbUzWEFglM2br8eXv3z1igE+NnYG1q/fhfn5W604+K1de2E6Wa++PkxZuxZNeqFqbIg9bDoG+3ZwLVv4q44pl2gQTpriLSqBkBAYTHYuJ8E84aPuYp43GY6PT0YjqJHm2FqoBpqHLqNyRilb+G2NKZdl1AxTIhNSwtTUTmzffhQ7dixh+/aj1ic/W+aXLD8EQLC4OB+0PZnYwZZjcBMVuW3fhjKfGw0mzZD9rcSYbA/rPrFt2zZz8OBB38UgkTGsbp2YOA/GAIuLDzpVvZ66571I/CBOje+xsTOC9vcg7Rjd6QPN+sTMzBiyI3MEO3YsVbqvyBqMjZ3Vajy0NWd0YQ7xbXIpQ0TuNMZsW3WdggEFA2KfrMlwGNeL9IEDm3LUmBuxfftRJ/ck+rGxUNXtW3mfH6ZrodWWkBQ6eYIBTQmKiSW0qI9k5z84hWvv5JAdn7qmT+PMhlmsroq8Sp/r2ls/5IiBLqBgoJSYQov6SJXJsOgzbRer2BJNuULTOAtFQKlrv6/a57oUWik4F0PBQCmUaMOmymSY9xkbi1XIjk9NabKwahlndZ65BgGijuYh2yl2NV0KrRSci6Fg4AlbufF9lI2UUzYZFi3SNhYrDV7ZXdJUmNKyc6z6zDVpOKoy2hcnJiYBnLbiM10LrX0UnOvAPAYeqBIH7CsGVkOMcgyM5j+oE5Vga7GKKdFUGU2T92iJNa/6zEM9vXS0L/r21u8iP0nIUDDwQJXB7StBR6gTj0aaLsy+Fyvfk3YT8hfWezE3tz+3/BoS4QDVn7kWDUdbNAitGsqgdazRlOABW7nxfZWNuMWnmjNEVTVQLDQVlV+LyaXqM6dtPB40jzXmMUD3eQzuuGNd5il4GmLMGf+uA187iVCff1neCO3lB6o9c8bfx4OGsZaXx4CmhI5JjsZ9eNV1kTUqHF+yVKuAYHLyMmf31KpO84kvNWeoGqNBW83OvjLz/arl99kXqzxz2sbjQfNYo2DQMYmX8eOrro+NnaVicE9N7cRDD/0zjh27EafSnhocP74P55zzPCcpQ+nsqAff/g1tmJraOZQOeiVVyh9KX9RgGyft0TzW6GPQMXnS4OLigx2XJJ/5+VsxmgvdVWy3ljhykhB6GFeb8rMvki7RPNYoGHRMCM5DNlRcVXMhaFan9REtznhNaVN+9kVShO38LprHGk0JHZNnwz958l4cOLBJhb2wrYqrjkpWszqtr4Suqg41TJToxZWZSetYo8agY1ZKicDw0bi2wlXaSrZtVVx5KtnDh3dbvxchtmBfJHn0zcxEwcADgzzjiXBg15ZvIza2rYorT/W6sDC/qhya1WnEPprTbbMvkjz6ZmZiHgN0n8dgwMzMGEYFgwTBjh1LjX5TQ2xs0fnrIcSTEzcwBj/Bd3iu7/uHiIZ51QV5eQyoMfDIxMR5mdfb2DQ1SLZFqtdYJWxSTt/UsVkUafS60Ka4yLanWQtki76ZmSgYeKJuoqOqg09D1MPU1E6Mj096LwfRNWlrEFp9kyccHTq0u5P0uLaFM81pfetQNk5smpk0jck8KBh4ok6iozqDT4tkOz29V0U5+oy2SVuD0GqDNhN7fh6T+U60KbaFsxi0QFXHycA3bMeOJWzffjRzni7rF2X30iI0UDDwRJ1ER3UGnxYHKi3l6DPaJm0tQmsb2gpbdYUg29oU28JZDFogG+Okar8oupcmQZ6CgSfqDNC6g69Msu0KLeXoK9om7RiExbaLSJ5wNDHRjenNtnAWgxbIxjip2i+K7qVJkKdg4Ik6AzSGwUe6R2O/CV1YbLuI5AlHW7Z0Y3qzLZzFoAWyMU6q9ouie2kS5Jn50BNVT0lLnBQfWfX90AZf1zAkKzvLJvtNO2xkRyzKdtdFn7WZbS+G0x5tjJOq/WJy8jIcO/aeVZ+bnLwM8/O3qsm8ScHAI2UDNO+M+YmJSWzZsjeowVeE7UW8TfrSmASKKpO2j/qG3MYuha3R5zVQIWtvG61pfatiQ7ip2i+SA+pWMz9/qypBngmO4C/BURldJNXwPUm7SHrTtN36loDHR31jaGNXYyaGtukzRf3i1HvZid8GSe26no/zEhxRMIBewcBFZsRhNExELoSfpu0Wa3azPHzUt29tXAe2jS5sLdJ5mt9hfD1jZj4MENfOYxq8YF043DRtN03OP13go759a+M6tGkbLfHvsWAzdDBrnh1Go98PBQPFuPb41TBJuxB+mrabRi9+l/iob9/auA5N20ZT/Hss2Nw0Fc2nWkN2KRgoxnXct4ZJ2oXw07TdYgi9qoOP+vatjevQtG2qLmJ91Co0rbPNTVP+mTgb1YbsMirBMW3tVC49fjV4wboKd2rSbjGEXtXBR3371sZ1aNo2VRaxNpE6oTA6105OXobjx/c1qrONsNRBmeqciaMFOh/CnfOhBue+MnxHJRBCiikbo1WcFrU4NnYZ0QEIspyQq9TZ1tyd1+7j45N4/vMfqPw7rshzPqTGwCFFKj4ti2/oMcjkFBTy4qPKTr+K5k+DP5FLrUW2g1/2prdKnW1ptuqciaMJCgYO0TAYST/og6q4j1TZXFRZxGypxtvgcqNUZ06tWmcbmyYN7d4EOh86RINzH+kHGkJPbdBHB7ki8hLijC6EZWdQuHL6rPO8XG6U8udUWfFX1z5UoTrbUjBwSKidgoRHDNopht2tJKm3ZL5Xd3PhIsKp7vOqs1GqKyDmzbUbNlzj9TTPUE8UpSnBIfTAJl0RqspymDytx+zsLgD9M4kk2p7sDJ5NNhe2/YnqmgaqRkE1MYtpnmtD9OOiYOCYEDtFbPTBKU9D6Glb8rUbi730l8hvD6OiHepqqaou3k19ETjX2oOCAYkam055mgUMzTumquRpPQB90TxdkK8F2uihNKtpoqWqsnjHYBYLHfoYRAadt1ZiyykvBPt3mQOadrLsxMP0bWGo66PU9dh35UNFp23/UDCIiBAWr66xtfuIxetfMwNHLWA88/2+LQx1HNdsj/0qQoYrxzo6bfuHpoSICCGhUtfYcsqjerMbBv00dH8JW1S1m9sc+3XMby7s+jGYxUKHgkFEcPFajS2nvBi8/kOBC0N9bI59DRsMOhL6hYJBRHDxWo2tRSYGr/+Q4MJQD5tjnxsMQh8DxdhK8jG8ePXROdGGU16oiUpiJOQ+7KrsNu3yeccE510n8UGNgVJcJPlgPv12cBfrn5D7sMuy2zS/5B242/YgXs3hvq4Itc48dhnujl1uQ9tjUrM6ZPK3/6NXCWmKluODmxBK2WdmxpCXcXHHjqVGvxnCEfS2CaHOeccu05SglDZ2vrzQpaoHshCilZDt36GU3UUegT6G+4ZcZwoGSmkzOPM6JOPDSeiEnPwmlLK7yCPQVigK0a8kFEEwCwoGSmkzOItyzjNxCNFKlck/5OQ3oZTdhaNtG6Eo1MRtoQiCWVAwUEqbwZnfITfSs94BIe5mtFF18g85OiSksttOr91GKApVJR+KIJgFnQ+h0/mwDSE4vcQC29oOoTjmkeY09dB34QzZFdqjEvKcDxmuGCHMHNcdGrLExUDI9lhSjabhvkXJm3wtvFXvG2qIMwWDQKg7AHx2SO1Ssk24oNmBWTtJHnlZRycnL/OS0yLkXBpV8eJjICJvF5EviMjnROTjIrIhvT4pIp8QkUdE5N0Nvr9JRL6bXv+ciNzYVZ1cEpLzTUhltUHIDkYuaOpvEbI9lrglzzdjfv5WL74Hofo81MGLj4GInG2MeTj9/+sAPNMYc42InAng2QAuAXCJMebamt/fBOCjxphL6pRHu49BSPbXkMpqA/oYnKJtW/RJ00Ta48v3IGSfh1FU+RgMFvWUM5G2sjHmUQB3iMj3Nfl+rISkrg6prDagP8cp8nZShw7trtQeodpjiR98mZ/6YPbyFq4oIteLyDcA7ATwFovff7qIfFZEPikiz7dUXK+EpK4Oqay2sB3aFSp5wt/i4ny0piTiDxvmpyamrz6YvZwJBiJym4jclfG6AgCMMXuMMRcA2A8g02RQRM73vwXgQmPMswG8AcAtInJ2TvmuFpGDInLwxIkTTarYGSF1xJDK6oq+5jUoEv5isr8SHbTNC9HUHyqkfBRN8Z7HQEQ2AvjYsF+AiFwFYFuej0HZ94femwFwnTGm0IFAu48BABw69FocO3YTgEUA49iw4WpMT9/gu1iZ9NlW3Gefg7m5/ZidfWXOu+HZX0nchOIP5XI+VXWIkohsGfrzcgD32Pi+iJwvIuPp/zcD2ALgSLvS+mdubj+OH9+HRCgAgEUcP75P7U60z6r1Pngs5zE1tRMTE5OZ78VsSiJhEoI/lK8oL18+Bu9IzQpfAPBiALsHb4jIUQDvBHCViNwnIs9Mr79XRLaVfP8FAL4gIp8H8CEA1xhjHuymSu7o82ITGiFMNi7ZsmVv701JJAxC8IfyNff7ikp4ecF7m3Kuv6bs+8aYvwTwl23Lp42+LzYh0QeP5SIYpUFCIS9xkiYh1tfcz8yHAdD3xSYkQphsXMOwQxICbYXYLnypfM39PF0xAOjpn41G7/8+eCwTEgtN/aHa2v6rzl2+5n7vUQkaCCEqoc+e/ln02ftfI+yfpE/UjWgYHh8TE+dhYeFhAI8vv180d/mISqBggDAEA7KSUEKN+gCFNNI36qRFzhofWfiYu1SFKxLSlq6dcjSaLbTAqBnShJDHVJ2IhqzxkYUmZ3IKBiRIugw16tuJkXVh1AypS+hjqo7tv+o4GB8/z0rZbEDBgGSiXZrv0ilH645YyzMKIR6c6ELrmKpKHSfjquNgaenf1MyzFAzIKkKQ5rv0/te4I9b0jBg1Q+qicUzVpWpEQ9b4yMKYx9QIRhQMyCpCkea7Sr2scUes6RkxRJPUReOYckXW+MhDi2DEBEfK8REGZkuajyWETWPSIm07LiY1InXQOKZcMjo+8qOqdAhG1Bgoxpe62IY0r0nV3YRh+/2RI3uwfv0uVTviPu24SHz41jL59s/Rbn5jHgN0k8egye7ZV6y+jbj0kPMMhBCXH0IZbRKL9on4R8vY0dCn8/IY0JTQAaMdcbB7BlDYEXypi20chKNN1V2HIvu9lsWoT4cVNR0/mtGwKPQVLeNbs/mNgkEHNO2IPg9PattpQz74KRShRvPEYhMtE7ktYhR0QiKU8e0T+hh0QNOOqN0OVUR2iI7g5Ml7ndn0bNkNab/XRWwTuaaIkj7C8V0OBYMOaNoRfTvotGFl2QFAMMgt7sIR0aazY8gCWYzENpHHJuiEBsd3ORQMOqBNR+wqVt8Fg7InwsFKJ1fbOySbu7CQBbIYiW0ij03QCQ2O73LoY9ABfXIUy6KLHZLte/TFfh8CIY+fLCfDvsXwa4TjuxiGK4LHLrumi9DFkMMjQ4Re9eUUhcUBYQo62mA/bAfDFYk3utghcRfWHaF61Xe9iBSZt0IzC2pkbm4/ZmevxLDvUvK37n4YAvQxIM7pwqZHu2F3uPSqd5WRzkcmTjoZuuWee34Ro75LgEmv68V31sUqUGNAOqELmx7thsXY2jHXXfCq3telJsJHLoSQc3mEgDEna13XQCjaNmoMCOkBNnfMdbzq69zXpSbCx+49tmgK24Swc7ZNKDksKBgQ0gNsTkh1Frw693W5ePsIERw1b01MTELkdMzOXtmbhTCP0A9Za0oo5iUKBqQ1fZT8Q8PmhFTHn6POfV0u3r5274NcHlu33oylpe9icXEefVoI87AhqJ577k/Vuq6BUHJYUDAgreir5B8atiekqamd2Lz5+tSO/nUcObIn85nXua/Lxdu3c2ooKuSusCGoXnrpbauEgHPP/SlceultrcrmklDMS3Q+JK2I7YCbWLEdzlnViarOfV0nMvLpnFq2EPYtHt+WY6ZmISCLUJJ1UTAgrQjFZtZ3bE9IVQXCuveNNbKkaCEMxVPdJn3OOxJCH6dg4Ii+7AAYkhUONiekOgJhCBOha4oWwj5q3ULZOfcVCgYO6NMOoM+Sf5+hQFiPooVwkK1vlNi1bhQY9ULnQwf0ydHIt1MX8UMoTlSayDsp1bZjKKOE2AZtocbAAX2zu1Py7x9UBdvDptatT9rKPNgG7aFg4ACqWUkfoEBoB5tCVp628vDh3dE+q1F/roWFR3rns2EbCgYOoN2dkLBo6ixsy8nYlpCVp5VcWJjH3Nz+6BbGLO1AHrFqbF1AHwMH0O5OSDg0TdKlMblXkVYyRh+nLA1JHlU0thp9E3yUiYKBI/IcjQghumjqLGzLydjmxF+klYxxx1y1TlU0tjYFPVvP1JfwScGAENJrmjoL23Aytj3xT03txPj4ZOZ7Mfo45dVpfPyJtTW2NgU9W8/UV4QbBQNCSK9pGi5oI8zQxcQ/Pb23N6GkmzdfD5E1q64vLp7E5s3X19LY2ooms/lMfUW4UTAgQaLRFkjCpGlOBhu5HFxM/H3ycZqa2omxsbMy3nm89kJsK5+EzWfq6zRGCgYkOLQ5fVFICZumC6mNBdjVxD/q4wQg2j66uPhg5vW6C7GtpF02n6mvRGIMVyTBoSm3PJOpxEHTcMG2YYZdhDbH3kdt5Y2xlS1LRvIAACAASURBVE/C5jP1lUhMjDFObxAC27ZtMwcPHnR6Dy2HKmkpRxtmZsYAZPVbwY4dS8t/zc3tx+HDu7GwMA8AGB+fxPT0Xqv1PXBgU86ktHF5p0ZIEa7HZOx9dFTwAZKF2Kf5JJR5VkTuNMZsG71OjUEHaJHYtZSjLVV2CHNz+3HPPa+GMY8tX1tcnMfs7KsA2KtvF85BoUwyMeCjrV1nkIw9RbvG9NyhZwWlj0EHaDlUSUs52lLF7nbkyJ4VQsEp6jslFeHaOSjLn2J29krMzIhzW3HffCe0+a7YwpcDW5cwb4xdKBh0gBaJXUs52lLF6auoTjbr69o5KDuzW2JGcblwxbpIFhGL4DwKT8LUjUYBnKaEDtByqFJeOcbHz+u0HDYoU9Xl1XXwns1yAO7UmGVCjAuny7m5/Zid3QVg0fm9NBGL4DyKRlU7SdBq3qXGoAO0SOzJ/U5bdX1p6d9USKk2yUt8Apxmvd1dqjGrCDG2/RmSiWkx8/3QF8kiNKncbe8iY1K1a9xhN0WrloqCQQdoSTgyNbUTExNnr7puzGPeO6JtpqZ24uKL34eJiVPpYcfHJ7F1658FNSlmCZWj2Fy4yg6lickuPYoWAb6PZpyqxNY2WrVUNCV0hBYv1YUFO8lAXGDbI1xLm7dhpRr4XgCC4VBN2wtXUT+I3S5dpHLvMlpBU54OF7Rpy9jaRouZeRQKBj1Da0fUamvTwLCA43qByvfNGI82re4wWcJk131T6y7SBm3bMra26SLBVRNoSoicUXvc5ORlKtSlo2i1tWnDta04T52+deu+6IWCPLrum5p8HWzTti1ja5s8MzPgN4U1BYOIybLHHT++D+vX7/Lu7zBK3Z1ATA5ImtDiD6MJV7vUvD6sxdfBBW3bMsa2yTrXwrcfBU0JEZMnnc/P36ouFWodEwfNDm6JwTfDJi7Mb1X6cIzhhW3bMua2GaDBj4KCQcSEZI/bvPn6VSmMRdZk7gQ0DBzSH1zYgcv6cKzCmY22jLVtBmiYt2lKUIZNFblGe1xR/UYP9Mo74EvDwCH9wYV5pa99mKaqcjTM29QYKMK2irxMOu/6wJii+iXOR4+PfOPxTC1AVXUkDx8itrC9S9UaHVSVNmMrxh2/zblGQ6QCNQaKsO39XCSd+0gUUlS/OjuoKg5IsSVCCQE6hFYnZCc6jq2V2G4PDVoVL4KBiLxdRL4gIp8TkY+LyIb0+qSIfEJEHhGRd1f4netExIjIuqFrbxKRr4jIl0XkJS7rYZsu1Ys+wgOL6ldHfVZl4DD8sVuqTo4UHhI0TP5NqTK2bD9nzf3GxVzjO4W1L1PC7xtj3gwAIvI6AG8BcA2A7wF4M4BL0lcuInIBgJ8G8PWha88E8AoA3w9gA4DbRGTaGJOd+F0ZttWLRap7HzbOovrVVZ+VqSP7asP1RRWHUEaTrCRUlXrZ2LL9nLX3mxjnGi8aA2PMw0N/nok0x6sx5lFjzB1IBIQy3gXgjRjODwtcAeAvjDEnjTFfA/AVAM+1U2r32FYvFk3WPhxciupnewelwYEnRvJ2blUmR9s7K827yJgpG1u2n7N27V+Mc403HwMRuV5EvgFgJxKNQZ3vXg7gm8aYz4+89VQA3xj6+770WhDYXhyLJmsfNs6y+tlUn4Vsw9VKkbmgyuRoc2dFO7d9qgpaZWMr77jzvOtlaN+RxzjXODMliMhtANZnvLXHGPNhY8weAHtE5E0ArgXw1oq/ewaAPQBenPV2xrXMmDcRuRrA1QBw4YV6JDub6sUi1b2vRCFdqU/7kAila4p2blVMQTZNZcxlYZc66vrysTWO7GO7xxuVresIjqIIg6L3YpprnAkGxpgXVfzoLQA+hoqCAYBnAHg6gM+LCAA8DcBnROS5SDQEFwx99mkAjuWU7yYANwHAtm3bsgPmHdFVGF3ZZB2qjbMqsddvQFf9qWjnVmVytBmGpX0XGRp1Ba3isZXn0tXM1avL8L0iAQlAofAU01zjxflQRLYYYw6nf14O4J6q3zXGfBHAk4d+6yiAbcaYB0TkIwBuEZF3InE+3ALg09YK3pK5uf04dGg3Fhfnl6+5dKSJUZIlK+nSMats51Y2Odrsj0VlYf6K+tgUtNau3ZjzbDbW/i2g23mszJ+hL1oqycsu5/SmIn8J4CIASwDuBXCNMeab6XtHAZwNYA2AfwXwYmPM3SLyXgA3GmMOjvzWUaSCQfr3HgCvBrAA4NeNMX9XVp5t27aZgwcPln2sFaMT+Chr125Ud34B0c+BA5tyJ2Hb/SmrD4+NneElzC6vLOvX78Lx4/tUlDEkbPYjTf2kLjMzY8i2Pg+s1NnrZSIMhSeIisidxphto9d9RSW83BhziTHmWcaYnxsIBel7m4wx5xljnmiMeZox5u70+mtGhYKhzz8w9Pf1xphnGGMuqiIUdEWWJDoMVaCkCV2q1DXF3ueVZX7+VtUe7Fqx7UAncvry/ycmJoMQCoDiCIMin4ZhJ9jZ2VcF7wTLlMgdUTZRhxzaUgTVum7p2jFLky01qyyzs1dmfpaCdzGj6vqJifNgTNKeA+fSKs89S1uwtPRdZ+W2TZk/Q5HW9xSP49Ch3WrGSROYErkjiibq0ENb8mBImXtiDJVqQ4wx5V0xCBfeuvVmLC19N/WFqjduteccKKNIK5b1Xh7DfmQhQsGgI7ImcCAsNVtdQp8kQkCTer8KrpMSZY8zweTkZVbvEzNtxm0M0SJF+VRG34sVmhI6oo8RAjFMElXxaTLRpN4voosIiqmpnXjooX/GsWM34pSjmMHx4/twzjnPC6KdfNNm3IZ+amRdJiYmsbCwWjswMTHpoTT2oMagQ+pk9osh3Wtf1Lo0mVSjKw3S/PytGPUer3qfGMZdW9qM276ZtrZs2QuRNSuuiazBli17PZXIDhQMFNJ2odEyucU6SYy27+HDu2kyqUBXGqSm96GAl9Bm3IZm2mrL1NROXHzx+1bU9+KL3xd8fb3kMdBGF3kMRilSPbeJKfYZQ5xVJyAu80lZPoqVCHbsWHJeplDoKudC0/t0mRNCO4wm6gd5eQwoGKB7wSB7cREAJjdr2OAzowvN6ABeXHwk0+blenILOalJHfIWjyz6uKAUUaWP2FiQmvbFouQ2FPDCgAJNPVQlOOo72cmOkgkpWXSyzoJabePLUn1mCQXJ77p1+OtLBELVdozBZGKbMjWzLVV+U3V2X3xiYoWmIHswKsED5YuLwUCDMCBroSnLpjiMi8ltWDrPSxUaWwRCntf1+PgkJiaeyJ1KCUURFDZPTGwSqdHlYT3EPjxx0x4UDDyQt7isZGBWyF9ofO5e5+b2Y3b2VQAeL/xcbLutvMVjenovJ5+W+A5v7WNIcUz47j8xQcHAA1mLyyhV7NM+d6+HDu1GmVAQ426Li4c7NMTAh5ITgqxGQ/+JBQoGHli5uAx8CorNBln43L0Wp/yUwgUzdAchLh5uoCqftIH9xx4UDDwxvLg0XSi17l6LPLi7yH7nihAFmpDKrLU/ayGkZ+kD9h97MFwRfvIYhM4dd6zLTQX64z/+QMY3EkKNFQ8xHDPEMpNstD9LCi1hwnBFYpWmqUBDdRAKMRwzxDKTbDQ/y9DDBLVkitUETQmkEU3VdqE6CIUo0IRYZpKN5mcZYpjgKQ3HSh+vkEybLqHGgDSmzqFQA7o+P8HWbiDE5DchlplkU/Qsy/q46x2xZqEli5UaDqDpgVsxQ8GAdEqXh6zYVHGGeCCUqzJT9do9ec9ycvKywj5uawwUPfPQBNBDh1YfejaKVqGmKygYkM5pomlogk27bIinxrkoc+j25FDJe5bz87cW9nEbY6DsmYckNM/N7S8JtU7QKtR0BX0MFFDVo5eev/WwreIMMX+B7TJrtCf3YVzk1XF29srMzw/6uI0xUPbMQwoTrCIQaRVquoSCgWeqxvWHHP/vi1AdHTWjzZ7ch3FRVMeyPm5jDFR55l0JzW2FwLJ+unbtRrVCTZfQlOCZqqo+zeFKWglJxRkKNuzJNn0UQh4XVduhqI5lfdzGGNDiQ2DDjJVX5omJSezYYZyaNkOCgoFnqu7AtO3UQiBEvwDXtF2U6yw0Wfey7aMQ6rio0w5FdSzr4zbGgBYB24YQmFeXsvwrTQnVUZemBM9UVfVRLd6MEP0CXGFD7V7Vnpx3L5HTrfoohDou6vhqlNWxrI+3HQNafAhsCIFd1iVkMxcFA89UPfiDB4SQtthyHKyy0OTdC8gOE2u6ww91XNRZ5GzXsYmdXoOAbUsIzKqLCwdWjY66VaEpwTNVVX1Ui5O2dKl2r/ubTXf4IY2LYbVy3tSb1Q426xhyuKnLvBwu2iRUMxfAQ5QA8BAl0g+6PMAq714TE5NYWvpu4WFAmsIPbZUl6xCkUbo4FCnUQ8wGuOgbLtpkbm4/Zmd3AVi0+ru24SFKhPScLp3Iipy8ina/mna0NsuSpVZOGEeXmo6Qd7FAs+RoZQ6Atttk0G+yhIIq402DwyJ9DAjpCV06XpXdK++emuyyNsuSv8gsYceOpYYlrE+ozppNqeIAaLtNioTAMuFPi8MiNQaE9Iiu0lE3vZemHa2Nsgx2f6MH9QzoekHO0uQAgpMn7w0qnK4KA3V+WYijbU1akRDY1Gm367wcFAyIMzSoxEhYaEmmU3TPqmVZfYrfSnxET6x0ZASyjhy2MU59j/0idT6wOmujTQfWNv1Gi2BMwYA4QZOtmISDlmQ6NsqSr1KG1+iJgSYnWQjtHzmsYewXtT2wepG2qUlr02+0CMYUDCzjW1LWghaVGAkLTeGHbcuSv8sTFal3Xe1ONYz9ojq4FjSz+s369btw5Mie0nVBi2BM50OLaHEc0YAWlRgJDw3JdAYMyjIIk5udvXL5jIKyMmp39MsrHzCGmZmxxs6pGsZ+ft3KHQDbMBpOuXXrzQBQeV3QkmWSGgOLaJCUtaBFJUZIFYo0fU1V41p2f3lMTl6GxMdglEW0MQHYGvtttK95bb916z6nQkFWPzl8eHetdaFLB+E8KBhYRIOkrAXtkyJpR0wms8R7/VUrJvTZ2Vct16mpwG/LLOKirefm9uP48X3Ii5YY0GRjY2Pst/VTqNr2XZz0ubAwn/l5zesCTQkWqao61JTZzRVaVGLEPrGZzA4d2g3g8ZGrj+PQod2YmtrZSuBvaxZx1dZlznnD1F3AbIx9Gzkkytredtt2lQa8CygYWKTKYSexTapFaLIV94GuBE5NSYhssLiYvaMbXPfpK+CqressYk3q2Xbsd6F9td22ef1kfHwSxqxOA65Ze0pTgkWqqK9C9EOISW0cK12GiPXNZObTLOaqrasv9uJlAevCR8l22+b1k+np4jTgGqHGwDJlknJok2poGo6uzTRazEJd7uKb7KC1tFMWExOTmXbgiYlJAH7NYq60FVnazdUINmy4xstzsnXUdFG/s922TdOAa4Qag44JzVs/JA1H14lVNCRyGdClwFl3B62pnbLYsmUvRNasuCayBlu27F3+25enuCttRZZ2c8OGX1nx99atN2N6+oZW9xmlqvbRhvNgWb8ratumWlINEQU24LHL6PbY5azjV7s4brUupyTt7HSugHR6+EsVuj5Stov7Vd1pd133vHJlXc/rR5qOn9Ws0dBctjrYnvvKfq/KmMhqWwBBzNE2yDt2mYIBuhUMgOKBrmESqHJ2vKZJfcDMzBiyw6/cCDFN7lfn+daZSDUInHllyO9H+oTLqmgYp6FhW3gt+72m80HXQnYeXfSxPMGAPgYdUeUht7Xn2+pIZaFMWj1qu/Yer3u/suc7+vwWFh6p7DegITw0z+wEjCPrMBut5rMyQvO70YJtc1fZ7zWdDzT4gfnuY/Qx6ICqNtY29nybdtyiAaDZo7Zr7/G69yt6vlnPLy+MLu/5+LZv5vebxaiSXYXkd6MJ2/5V+d8zOHBgEyYnL2vU7zT4gfnuYxQMOqDqQ24jqdrsSPkDY6Nqh5quD+Cpe7+i51sn4YzWnXZRvwktXKsIDTvKELEtuGf93oCTJ+/F8eP7sH79rtr9TkPWVt99jKaEDqj6kNuowm12JFuhQj7oOqlSnfsVPd+qz0nzcyjqNzElu9J+OJJWbJu7Vv7e6uextPQdzM/fWtsvoEo5Xdv/ffcxagw6oKpqSss53pqOvo2Joueb95wmJiaDeQ596Tdtxmnfk4XZNncNfi/7MKjmO+yicnYRfutba0GNQQdU3YG3kaht7/Jj2uFpoez5Zj2/LVv2BvUc+tBvmo5T3w5lIVJ1Z97lDruLZGK+nYkZrohuwhW7CD1hCFXY8PmFR51npiUMLhS0hut2HRbtEuYxKKDrPAaEkPCpuxjFtKB0QV1BqivBOiYBL08woI8BIRnEaAuOsU4+qRsJpCEMLiTqOlR3Fa7r2/7fBRQMCBlBe27/JlSpEwWHetRduEJdUHz1C62CVB+cbCsJBiLyoyJy1tDfZ4nIj7grFiH+8J1cxAVldYpRGHJN3YUrxAXFZ7/QLEj5Tibmmqoag/cAeGTo70fTayRS+rx79J1cxAVldYpRGHJNk4UrtAXFZ78IUZAapmwO1TzHVg1XFDPkpWiMWRIRhjpGSt/DqnwnF3FBWZ1iFIZc4zukrAt894tQw1+rnIuieY6tqjE4IiKvE5HT0tduAEea3lRE3i4iXxCRz4nIx0VkQ3p9UkQ+ISKPiMi7K/zOdSJiRGRd+vcmEflu+rufE5Ebm5axzxw+vLvXu0fNKsymlNVJqz1XO6FpAOrCftGMMk2Ldg1dVcHgGgA/BuCbAO4D8CMArm5x3983xjzLGHMpgI8CeEt6/XsA3gzgurIfEJELAPw0gFHR9avGmEvT1zUtythL5ub2Y2Gh3uE9sRG6CjOLsjrFKAw1QbN61wfsF80o07T41sSUUckcYIy5H8ArbN3UGPPw0J9nIg3uNcY8CuAOEfm+Cj/zLgBvBPBhW+UiKJRY+7RLCFWFWURRnfqgFi9Du3rXB+wXzSgz3Wk3V1YSDERkGomz4ZQx5hIReRaAy40x/7npjUXkegC/AOAhAC+s+d3LAXzTGPN5kVU5sp8uIp8F8DCA3zbG3N60jH2kSGLlLiFuYhSG6tBFqtsQ6Xu/aEJZinrtB9VVNSX8KYA3AXgcAIwxX0CJBkFEbhORuzJeV6S/sccYcwGA/QCurVpgETkDwB6cMj8M8y0AFxpjng3gDQBuEZGzc37nahE5KCIHT5w4UfX20ZMnsY6PT1qZHEJQ1YZQRmIf7epdEg5lpjvt5sqqkQVnGGM+PbI7Xyj6gjHmRRV/+xYAHwPw1oqffwaApwMYaAueBuAzIvJcY8xxACfT+98pIl8FMA1gVb5jY8xNAG4CkpTIFe8dPXmS7PT03ta/HYKqNoQyEjdoV++SsCjTtGjWxFTVGDwgIs9A6gsgIv8rkt15I0Rky9CflwO4p+p3jTFfNMY82RizyRizCYkz5A8ZY46LyPkiMp7eYzOALWgRPdFHsiTZ9et34ciRPZiZGcPtt6/DJz/5RMzMCGZmBHfcsa7yjlq7Jy4QRhmJG+hoR0hCVY3BryLZXV8sIt8E8DUAbUSdd4jIRQCWANyLJOoBACAiRwGcDWCNiLwUwIuNMXeLyHsB3GiMKTrt6AUAfldEFgAsArjGGPNgi3L2kmFJdnQHvbi4MmJhYWEe99zz6uXvFaFFVVt02IqWMo7CkxfdQ0c7QhIKT1cUkTeMXDodiZbhUQAwxrzTXdG6g6cr5pN3ktgoVU4W03AqWdmJeBrKOEqXR8pWKQsXTkLioOnpimelr20AfgXAkwCci2SH/0zbhST6qLpTrvI5DaraMlOBhjKOosW8wfMUCHGHJqfnQsHAGPM7xpjfAbAOiR3/OmPMbwD4YSROfyRyqjpeVfmcBk/cMlOBhjLmla3qdVdoEVAIqYOmBTcPbUJ3VR+DCwE8NvT3YwA2WS8NUUdWlMIoImsq76h9e+JW8Tz3XcZRtHjLaxFQiF1iNg+FEmWkLYdG1aiEmwF8WkTeJiJvBfApAPvcFYtoYXQHPT4+CZEzl9+fmJjExRe/T9UgK0KjqaCMNmW2uVti3vz40LZTtU0oWi5tQnfVlMjXi8jfAXh+eulVxpjPuisW0SDFj5Zh69abgxEA8gjR87xpmW3vlrRnayP10bZTtY22BTcPLVrBAZWPTjbGfAbAZxyWhaRoUH9pKIMrtJkKqtCkzLYn/VCEKg1CdSiEsnA2RduCm4c2obuyYEC6Q4MUr6EMACf5NriY9LULVTELtC4IZeGsw/CcMTFxHoDTkGbzB6BTy6VN6K7qY0A6RIMUr6EMsds/XdNHn4BQbMpaqOq/EoJnP7B6zlhYmIdI4hulIcqoqB2npnZi+/aj2LFjCdu3H/UqyFIwUIiGCV1DGTjJtyNER8u2aBBoQ6JKeG5IAnrWnGHMY5iYeKL3BTekdqRg4IC20rWGCV1DGTjJt0NjTgbXaBBoQ6NspxqCgD6Yc/OytGqYM0JoxwH0MbCMDRunBnuThjLk2T8BgwMHNtHfoALafQJso82JKwa0C+hZKcNH0SAYam/HYSgYWMaW056GCd1XGU45D90LQJAe6rkCOpWRLDQItLGh3UExa84dRotgqL0dh6EpwTJNpMJQHHu6YKUdDkiEAsn8bBM1HNs6fjQ5ccWABrNiEUVzqybzWdN29DFnUWNgmbpSIcOrVpIt/eefAFpHDce2JmUwPHY12rUw+XOuvxNRs2jSjr7mrMJjl/uCzWOXs+xdImswNnYWFhcfXNUZNB7z65OZmTEUCQKj1GkntjUpQtPx1qQ6MT8313NW02OXSU1GPcEnJiZhjMHi4jyyQlRCckjpgjzNyvj4ZGt1Jtu6O0I02ZR5jYdYpz4Qc/SNrzmLgoEDhm2c4+NPxHDWLWDlZNOX8Kqqk2qeHW56em/rwd+XtvZNSPHawxRNwqHWqS/E6lfia86iYOCYMolPu2OPDepMqkXSf9vB34e2tknTHXKX8dpdnR4ZUgw6iQdfcxYFA8eUSXwhqMHaTr51J1VX0v+grZP0qAljY6db+W1ttH1mbXbIXak/be/iiyZhmqGID3ytD4xKcEyVhCsachbkYcMrVtukasx3l/+/sDAfZGRCkfe8jWfWJh9HV/HaXZ4eeSqvxkpohiKu8bE+UGPgmBA0AkXYUKFqsu23qY8W57OynbKNZ9ZGmOtK/enq9MgsbRXNUATQMwe4hoJBB4TsGGNj8tU0qTatz6FDr8Xs7JUqnM/KFn4bz6yNMGdLGC6bhLsUOEMX8ENE2yLcJwdUmhJIITbUwpoSpDSpz9zcfhw7diNG8ysML8Zd1q1s4bfxzLJMYIDg5Ml7K51T0Vb9WcUc0vW5CJpNfrFR1xzWRWKqMtNV2zJoSq5FjQEpxNZuX4vWpEl9Dh3ajbykS4MJq8tdRNlO2cYzW7lDBobPrOiijlXMIdzFx0sdc1hXO3mX4azatBEUDEghWZPv+vW7cOTIHjUqvjrUXUzm5vanyanyGO88jK1s4be1YA6EueR38rUlLqhqDtEicBK71DGHdRVK6jKcVVs4LE0JpJRhFWoM5w3UUQkXD0wBsJj5jsuIiyqmGZtqbx9RJSGdRFeGJhVxKNR5/l31zyLT1ezsla3KoC1yixoDUgttkq1rigbmhg3XDKnaV+J6Aetyp+wjqkSTw2obtKmIQ6HO8++qfxZp4tqWQVPkFkDBgNREm2TrmryBOTExienpG6JZwIrwUcdY/Af6Jkjbos7z77J/ugpn1TaP0JRAMslTf8ak4q1Cnvpwy5a9AHRFXLjCVx1jiALomyBtk6rPX8MYbFsGDXUYhscuw+6xyzFQdIwpgGiPOM2DNmKSRZV+waO+iWbyjl2mxoCsokj9OZjM+rRQxrBzJXap6oTbda4FQmxAwcAx2nebWeUrU39yoewX2vuwD6qe06BNRUxIFSgYOER7aF9e+cbHz8uM3Y/Vj4Dko70P+6KO7wAFaXdQaHUDoxIcot0jOa98IlDlIUv8ob0P+6Lr8DJt5wZogKGg7qBg4BBXHsm2Jom8ciwsPBhFqBhpD73qs+kyvCxrAZydvRKHDr3W+r1CgkKrO2hKcIiL0D6bqt288gEGR47soVqOVO7DfVPpduk7kLUAAgbHjt2Ic855XtTtXASFVndQY+AQF7sKm1JyVvkGUC1HgGp9uK8q3a6yT+YvdCbI3bEtjae2bIExQcHAIS6yt9mUklefoLeSPqvlaNNNqNKHqdJ1S9FCp2l3XGXM2BQitWULjAkmOEJYCY5cJUyZmRlD9tHCgh07lhr/bogUJXjqq9q2CPYdt8zN7U8P6VndxnXGvUtzT9UxY3v+6psJyzZ5CY6oMQgMV1Iy1XKn4A64Huw7bpma2okNG65BcprnKeqMe9fmnqpjxrZfAI/ddgMFg8AoUu22UX9TLXeKOpMXTQ5h9p3Qntv09A3YuvXmxmZJ18Ju1TFDITIMGJUQIFkJU9pGKzBD2ynqeOLXbfMYVZ+u+o6rtgo1aVObREmuPfirjpnQUkTHOF6rQB8DhOVjkAcPa6lO1mAHTi1sExPnYWHhYQCPL3/Hhr2UvgvVcdlWfRwrrutc53mFstj2Ybzm+RhQMEAcgoEGBzANA76sDFmDHTgNIgJjHlu+IrIGY2NnYXHxwdy61G3zPi5ITXHZVhrGStdk93sBYLB27UbVGh5f9GG88nTFyHGRTKkOGtSzVcqQnSzmcYzKx8Y8homJJ+L5z38g935125wJWarjsq18jxUfrDT33IuBUADYG6uxnQnR5/FK58NI8O0ApsGTv0oZ6gzqss/WbXM6XlXHZVv5Hiu+GHjwJw6MKyVhRt2sps/jlYJBJLhIplQHDdJ1lTLUGdRln63b5n1dkJrgsq18jxXfaBirIdDn8UpTQkT4VOVpUM9WKUOWV3SWj0HVCaBOmzPyozqu2yo2tXcdNIzVEOjzeKXzIeJwnZSMzgAAIABJREFUPvSNBg/eqmUoi0ro0wRA+oeGsUp0QOfDHtOFt7AG6bpqGfJ2i5wUSR/QMFaJbqgxQNwag1B2B7GFOmkgxjaNsU6EjNJVP+dZCQFhM12rhmiBMvp6bG8d6vaJGNs0xjoRMoqGfk7BQBm2O0UIHsghCC8+yesThw69NldYiLFNY6wTIaNo6OcUDJRhu1OEEIvrS3gJ5SCdvD5x7NiNuQJkCAJhXTTVKZS+M0yIZe4jGvo5BQOPZA1U250ihFjcMuHFxYSmQV1Xlfxnn5+kpolAqH3h0CLkhtR3BoRY5r6ioZ9TMPBE3kAdHz8z8/MTE+c1uk8IyVyKhBdXE5oGdV1V6kwIAyGirkAYwsKhRcgNqe8MCLHMfUVDP6dg4Im8gbq4+Gjm59sEjwxSoe7YsYTt24+qEgqAYuHF1YSmQV1XlayJIsl1v5qBEFFXIAxh4dAi5LroO661NSH1976joZ97yWMgIm8HcAWAJQD3A7jKGHNMRCYBfAjAcwC83xhzbc733wbglwCcSC/9J2PMrel7bwLwiwAWAbzOGPP3LuvSlKrq4QGLiw+6K4wC8nILuJrQQsr+lhV3Pjl5GY4f31d4rn2d7H6hLBwaMhba7jtdHEAWUn8n/vu5L43B7xtjnmWMuRTARwG8Jb3+PQBvBnBdhd94lzHm0vQ1EAqeCeAVAL4fwM8AuEFExu0Xvz11B2TMA7hot+TK3qZBXVeHUa3P9PQNVncVGuyaoWC773ShrWlSZu0+J02IsU4u8CIYGGMeHvrzTKTbZGPMo8aYO5AICE24AsBfGGNOGmO+BuArAJ7bqrAWyOqMyYDMVgePXte8YLWlzLbtagHvQl3nehKyaSIKTVDyie2+04W2pm6Zm/qcaF54Q/Cj0YK3zIcicj2AXwDwEIAXGmNODL13FYBtJaaEqwA8DOAggN8wxnxbRN4N4F+MMR9IP/dfAfydMeZDRWVxmfmwKPPg7Owrc7+3du3GXmR3O3BgU46KcyO2bz8KQGe2u7IyhZJxchiN7dwHqoyBrmlSJu19XmM7+6bzsxJE5DYA6zPe2mOM+bAxZg+APalPwLUA3lrj598D4O1INA1vB/CHAF6N7C14puQjIlcDuBoALrzQnbq0SE2YLP797qhVdku+7W2jVLEJFz13LXXJEgT60u80kXXip29tTRMthvY+H4ofjQacmRKMMS8yxlyS8frwyEdvAfDymr89Z4xZNMYsAfhTnDIX3AfggqGPPg3AsZzfuMkYs80Ys+3888+vc/taZC38g+shq29tqQxDtG1XsQlrn4SoVtWDBi/0UZqMS+19PsS5xhdefAxEZMvQn5cDuKfm958y9OfLANyV/v8jAF4hImtF5OkAtgD4dJuytifP93Fc5YRQBZuLSojCUZUJUPskFEJ4Yp/QFlLcZFxq7/MhzjW+8HXs8jtE5CIk4Yr3Arhm8IaIHAVwNoA1IvJSAC82xtwtIu8FcKMx5iCA3xORS5GYCY4C+GUAMMZ8SUQ+COBuAAsAftUYs9hdtbLIu31yXZuavAo2VYYhHgFbJfRLo3p4GO27O+KXJuNSe58Pca7xBY9dhlvnwxgdXmZmxpDtuiHYsWOp6+J0TpaTVeLeYrB27cblyUazM1+M/ZL4R3OfJ6vp3PmQJHQlRXc5IPueLGXlzmPQDomgNOqIqHVS1L67I2Giuc+T6jAlsmO6ipfv0pGMtrrkueblogjBVt+0X9ZxOtUc004IyYemBLg1JXSBD7UwVYb57Z4Qn1kl24QCjI9PYnp6b/B5HAjpGzQldEyXC6cPRzKqDIvb14VZxbcwluV0CgCLi/NB5nEghGRDU4IDulbt5y9CY1bvSdXwSvLbXZz4kPjOO1AkCIWWx4H0D85f1aFg4ICuY8Szj+UFgMXGi8foIDp06LW5C1NfB1zeccgbNlxjfVesIe9AmRYkpDwOpF9oEKxDgoKBA7reLQ0cybKSKTVZPLIG0bFjN2YuTIcP7+7tgMty4Nu69WZMT99g/V4aduD5AmjCaB6HvjuoEj1oEKxDgoKBA3zslpIdarazW93FI9uWnO2kurAw3+sB11XGui7MRWWan6mpnVi/flfmd0XWrFj0XUbj9FVDRZqjQbAOCQoGDvC1W7IlkNgYLFm/wQm9OS7MRcNUVbXOz9+a+f2xsbNWLfouhCaqhEkTXG3WYp3TKBg4wNcZCFmLh8gaLCw8UqvjFjnVDTM2dgbGxycr/QYn9HbYNheNUlXVmic0Li4+2Or+VaFKmDTBxWYt5jmNgoEjfByKMiqQTExMwhiDxcV51Om4eYNow4ZrVgk709N7Kw04TujtsWkuqvr90eu+nQqpEiZNcLFZi3lOYx4DB/iMNx/OL3DgwCYsLMyveL9KLHmTw0bKPssJ3Q6u0lFX/V3fqZT7no6bNMd27pWY5zQKBpYZzfg2mju/S9p03DqDqMpnOaHbwdXCXPV3fZ9Q51swsYHvRFXEDjHPaTQlWEaTesm32ncYhq/ZwZX/Sp3f9WEma1JOjcRsl+4bMc9pPCsBds9KqHskscvdg7Z89dwpkb7D467jIvQ5jWcldEQd9ZJrs4NvtW9WeUIaNIS0ZXThyDt0Kwa7dB+JdU6jKcEyddRLXZgdhtW+mzdfjyNH9kQXc0uIbWzEp2eZDbKO6QbisEuTeKDGwDJ1dulderVqcookRDO2xkp+BlHBsLkxFrs0iQcKBg6oql7q0quVx+C6J3R7I0mwNVbyBXyDtWs3sp8QtdCU4JEuvVpjjrnVAL3N48HWWMmPCtq4IqoDQJRpdUm4UDDwSF7oFWB/otAUuhgjmsJUSTtsjZUqgj8FSqIRCgaOyXNiGlyfnb0SALB1683LuwcXE0XMMbcaoEYmHmyNlSo5FyhQEo3Qx8AheU5MDz30zzh+fF+mc5MrXwBtoYsu8GnjjzkLWt+wOVbK/I0oUBKNMMER7CY4GiYvmUlyQt7iqqsDh6Q6CZJIgu9kTr7vT8KECY/CJ2Sn47wERzQlOCRf6l8tFAw+b9MXYNiMcfvt63DHHeuidXDqSiWbZxoKPVUv8QNNfGETq48ITQkOyc90lqcxuNDaITGjO9jk6OUETTkMbEnbXahky+LbY82CRtzRBxPfKCHvsEeJNQycGgOH5O0GNmy4OneXYGvnmZ1c5RQaHJxsSttdRF3QUYy4wOehVF0T2w47Vh8RCgYOyVvkp6dvKFz8iyaKqqlaq3RM353XxkI7aI+sdLO2VbKxTgKEdEVswnWsYeA0JTgmT73cRO1cJ1Vr0YEtw5/xSduFdrXD36l0s2vXbrSuomTkAalCTKpy29Qd89rb0pbpVxvUGAREHWk7y4wxjMga7523rbSdl4t+4NFtewLR4ihm44Af4obYVOW2qTPmQ2jLWJ2OKRgERB1pe9BhE0fH1YyNneW987ZdaLtW7WuYBKpMli4FB1e/HYuwE5uq3DbaTp+1QYw+IjQlBERdVfbU1M7lzIqjLC4+aLVsTWjrkZ3XHuPj51kt5zC+Iw/KvKBdnqLp6rdjOvmTfijFaD19lqyEGoOAaLLD1u4c00baTup92qrrS0v/FuyOs4yyydLlLsvVb9v+XZ/aB+3jTQNVxzzb0h8UDAKiiSpbi13cBVNTOzExcfaq68Y8pk7daIuyydLlLsvVb9v8Xd926ZjHW9f4bss+m80oGARG3R22Bru4SxYWsk0irtSNvgd12WTpcpeV/xtjrdrBZpl926VjH29d4rMtXQmYvgXXqvCsBLg7K4G4p8tc81rOQygK4XJZxqzftnEPm2WemRkDzxohbXE1r2g7G4NnJRAV2N5xd6luzNuNzs7u6lSDUKQ1crnLKop0abMrt1lm2qWJDUIwm7mEUQmkM1x4n3eZa77sUCwt3vQuIyeKIl3aTG62yhxrwpkQ0J6MqA6ukpmFkiSNGgPSGa7sv13FEVcZvBrjrG2jeVdOG78fQrGdV8WVJtK3Q2VVKBhYxrdzmmZCUaPlUZZNckAo9WmK9sktxoQz2vHt9GkbVwJmKIIrTQkWiSlRiwtCUaPlMWq2SOTq7OOzY6aPRwWTYkIX+rNwZZLznSStChQMLGLrbO5hW93ExHkwJslUGPoEHIP9d3hQ53nTh1SfpoQwuZHuCF3oJyuhKcEiNqTmubn9mJ191bKtbmFhHouL84jBbheKGq0qsdWHkKZoNy+RejCPAezlMbARo3r77etSQSAfXzGvpDoxeWj3ET6/+rDNwiMvjwFNCRaxoSovEwqAsO12odBmkqOvSdh0+fxiWkxpXooHmhIs0pVqWYvdLtYIjDqhV1ltUNdDO9Z2DJWuPOxjC/Ej8UDBwDKDUKmtW28GAMzOXllrsp+YmCx8X4vdLuZJrerCkNUGs7NXZpqTgGxNT8ztGCpdedjHFuJH4oGCgQPaTPZbtuyFyJrM9zQ5t8U8qVVdGLLaIDtPf0KWpifmdgyVJgmcmmh9YgzxI3FAwcABbSb7qamduPji960wR2zd+gHs2GFUJWuJeVKrujAU11VW/JWn6WnTjjRBuKGuh33TjYDmDJKk31AwcEDbRTOEzG0hTWp1F9CqC0NxXU0lX5Om7UgThDvq+go13QgwxI9ohYKBA0JaNJsSyqTWZAGtujAkdZXM3xiElJYJd03bkSYIO+QJjXWE86YbAebBKKcPWjGNdWQeA9jLYzDA5vnymgkh1Mr1+eeHDr0Wx47diGHfgrrPukmmy5mZbIEEAHbs4Jiugq1x6rqP9ZU+zKO+65iXx4AaAwf0ZScQgsnDtS/E9PQN2Lr15lbPejiSZWnpuxUzXY7n/FredTJKG63L8C5vYeGRVQ7DGrVnbfCxq+2DVkxrHZngyBGjyT4GA0vz7jpGusjhbiuxS72zNlYf3lR8nYzSVGgc3eUlgtxpmJiYxMJC+GeajOIrYVfMDs4DtNaRGoMOoKOYP0LxhQDqTRKJhmI1edfJapr6AmWHqT6O8fEnqtaeNcXXrrYPvlpa60jBoAO0qov6QEhmnTqTREgCj1aatqHWXZ4rfNW3D31cax1pSuiAvk0k2gglh3udszYG9dHu/KmZpm1YZp4KwSm3Dr6OVO5DH9daRy9RCSLydgBXAFgCcD+Aq4wxx0RkEsCHADwHwPuNMdfmfP9tAH4JwIn00n8yxtwqIpsAzAL4cnr9X4wx15SVx3ZUwij0WiZViW1RiZEiT3IA0XnS+/acD5UQxnJeVIIvweBsY8zD6f9fB+CZxphrRORMAM8GcAmAS0oEg0eMMX8wcn0TgI8aYy6pUx7XgkEIAyuETkzq4+K5jv7m5ORlmJ+/tVd9J69dY9gEZNUN0Ler1UwIcz6g7NjlgVCQcibSIHBjzKMA7hCR7/NRLldoVRcN4DHBceLiuWb95rFj71l+vy99J888FbrZMK/PXHTRTdi+/eiy0DA7eyUOH95dKedGH6kXYaQPbz4GInI9gF8A8BCAFzb4iWtF5BcAHATwG8aYb6fXny4inwXwMIDfNsbcbqXALdFs5w69E5NsXDzXbI/8lfS57/iyx9uizFF6WGhYWJhf/kxfBMKqhC4gOotKEJHbROSujNcVAGCM2WOMuQDAfgCZJoMC3gPgGQAuBfAtAH+YXv8WgAuNMc8G8AYAt4jI2Tnlu1pEDorIwRMnTmR9pHN8pcYMvROTbFw816rf7Wvf0eplXpWiPlMmFDLS6hRawxCr4kwwMMa8yBhzScbrwyMfvQXAy2v+9pwxZtEYswTgTwE8N71+0hgzn/7/TgBfBTCd8xs3GWO2GWO2nX/++XWrZ52quQ7m5vbjjjvWYWZGMDMjuP32da0FiNA7McnGxXOt+t2+9p2QwmOzKOozVYS9pgKhxvMC2hC6gOglj4GIbBn683IA99T8/lOG/nwZgLvS6+eLyHj6/80AtgA40q603VAl18Hc3H7cc8+rV6jwFhfnMTv7qlYDKfROTLJx8VyzfnOUvvedEFKF51HUZ6oIe00EwhgTwIUuIPryMXiHiFyEJFzxXgDLIYUichTA2QDWiMhLAbzYGHO3iLwXwI3GmIMAfk9ELkXitHgUwC+nX38BgN8VkQUkuWGvMcY82FGdWlFF7XvkyB4Y81jGpx5vZdPV7hxJmuHiuWb9Zh+jEmKlrM+MetoP01QgjNXHSbNfWRk8XRHuwxWrUCXMaWZmDMOn+K1EsGPHkrPyEUJIk5NAy8if1zinuUZVuCJZTZWsd3kez4P3CKkDc1eQOrjqL20iOdiH3cCzEpRQxSa1efP1q453TTit1zbd0GjiaGXbOStGuy5xh8v+0tQXJsQ+HIqTJU0J0GFKqMrc3H4cPrx72QFxfHwS09N7KSUHQpOMaDazqJ3aYeVpnsLJ0Ee6w3VGxyY7/9CyTGrMhqgqJbI2QhIMSNg0mcxsTYBZE9NqaNclq9HoB6CxTEVoFGTyBAOaEgjpkCZJh+p8p0hVWSVrIX1VqhGKStgWGnOdaCxTESElkqNgQEiHNJnMqn6nzOZaNgFpzD+gcQEO0bbdFo25TjSWqYiQBBkKBoR0SJPJrOp3ypJkFQsf+hKwaF2ADx/eXZqMLDY0JuzRWKYiQhJkGK5ISIc0STpU9Ttlqsq8kFgfk2kVZzONiW/m5vavyDw6jEaVsE00JuzRWKY8QkokR8GAkI5pMplV+U5ZPLiWianqcdAabbJFWgGNKmGii7xxrC0fAwUD4gxtnT12qiTJ0rDDqqoJ0HiEcZFQolElTPRTVVDuEvoYECdotQ/HTCg216qaAJ822TynxzyhZHx8Ul07x4JGB1SbVDlAr2uoMSBO0Ggf7gMaNAJlVNUE+DJ9FO3g8rQy09N7nZapr2jcTdtGo8mMGgNinbm5/bmZ9WJ30AqRrndkdTQBPo4wLhNqQ9DKxEKT3XRoGgaNYYzUGBCrDCT8PEJz0IrdT8LHjkyLE2QeZTu4ELQysVB3Nx2ihqGKb1DXUDAgVinKrue7s9clxEmmLr5MPpoXV41Oj32l7rMI0YSpUVCmKYFYpchUEJrKVaNTkG002jd9E1IimjY0Vbl3qaqv+yxC7c8+TGZFUDAgVsm3l2303tnrEuokUweN9k3f9MGPoGnUUNfRRnWfBfuzHXi6Ini6ok00Hi3aFI2nodkmpudFqpPXt4Gkf+epsrWPCfbnevB0RdIJMe22+qBSjul5keoUn+aZrwXQrkVjf7YDNQagxoDkE3tUAuknRRqDAVlaAO0aA1KPPI0BoxIIKUCz9zwhTckKkRslSwugMbSO2IemBEKUozlhi+aykXxWqtyzyXLYo6q+H9CUAJoSiF40O1NpLlsfsGXm4nPsL3Q+jBzu3OJEcy4FzWWLHZthg9QCkFHoYxABXWfoc+mQR2e/lWj2AtdcNg247Mu2M/zV8aXhGI0fagwioMudm8sEJzyqeTWaE7ZoLptvXPdlX0IZx2g/oGBgGR8q/S4nCZdCCFXTq9GcS0Fz2Xzjui/7Eso4RvsBBQOL+JKmu5wkXAohVE2vRrP9V3PZfOO6L/sSyjhG+wF9DCzi62SvLmOLXZ48F8Opdi7sr5pzKWgum09c9+WpqZ146KF/xrFjNwFYBDCO9et3lT6Ltv0zhjFKyqHGwCK+pOkud24udyqhq6ZDs782MXsx+qUarvvy3Nx+HD++D4lQAACLOH58X+HzsNE/Qx+jpBrMYwB7eQz6ki6UUQnZhPT8m8SuM969Hi77cpO+Zqt/hjxGyUry8hhQMIA9wYATZ7+ZmRkDkDWeBDt2LHVdnEJsLizAOIAlLhId0qSvhdQ/QyF0IYkJjjqAzlj9JqTwvSZmr/z3FhGC6cQHrkwvTfpaSP2za5qa1UIyHdaBgoFlpqZ2Yvv2o9ixY2l550WbbD8Iyf5qc2EZhqFrp3C5cDTpayH1zy5p+pxiDt2kYOCQmCVKspqQNEa2FpYsGLqW4HLhaNLXQuqfXdL0OcUcuslwRYf4Cl8k/hg814HdcTC5aHveo+WsYh8d/U6yr1jM+KTBgQObgrO32sb1wtEkVJThpatp8pwOHXotsv014jDNUDBwSMwSJcmm63Mr2tB2Yclyth2gud5dwZj/MKj7nA4dei2OHXtP5nuxmGZoSnBEYi7Ibl5ODPESs91xlJWq6dXEWu+q0KYfBnWfU5JUKptYTDMUDBww2EllqVk5McRN37REA2dbQDLfj6HeTSMLaNMPg/rPKct8ljA7e2UUTuY0JTgga9eYMM6JIXL6qj6Otd5tTUO06YdBlec0yFlQjInCjEaNgQOK4r1jkShJNn1VH8da7z6Zhkg+KyPMygm9j1AwcEDxLolhizHTV/VxrPXum2mIZJOvBc4n5D5CU4IDsk47HKVN2GLoaThjp6/q4xjrHauJhNQjf5GXKPsINQYOGN095dFEomTSpLDgaYRhE6uJxDehjYuiTKEx9hEKBo4YTo2cF86V1dnKBgxtnuHgS4jzNemGNtlXYWpqJ9av34XkoCgAGMf69bui04x0SYibm6LFP0YzGk0JHZBlWsiSKKt4QNPmqZ9Tpp7V6kXXmS99JVgKKbFTHebm9uP48X04FaK2iOPH9+Gcc54XdL2KcGGqHP7NrIyZ2jPClmUKjc2MRo1BB1SVKKtoA3hCmm6qeC+7FOJ8aZRs3FejxqFvGjoXu/nR38zLA6B9czN6QF5MgsAoFAw6okqnqqINiNGeZQsNC0sV72WXQpwvjVLb+7ZZkFw+975p6FwIQlU9+n1vbjTMH1qgYKCIKtqAGO1ZNtBityxbMFwLcb40Sm3v23RBqvLc20z4fdPQuRCEqnzX9+ZGy/yhBQoGiqiqDeiTSqsqWlS+RQtGF0KcL41S2/s2XZDKnnvbCT90DV1dociFIJT/3XFo2dzk9aPZ2V29FA4oGCiC2oBT1J3QtKh88xaSrVs/sCzEuVRZ+upDbe/bdEEqe+55E/7hw7srPYOQx2QTociFIJQ/Jvap2dwUZavto+ZAjMk+U7pPbNu2zRw8eLCTezE5UTlZx/mOjZ1ROCEfOLApJ8nIxvSQn+4oesZN6tYHmrbL7bevw+Li/Krr4+OTeP7zH8DMzBgSh7diYnwGTceE66gEjfNeXlsNcD2P+GofEbnTGLNt1XUKBt0JBlwUqtFkQgulbTUJMNpoMjneccc6LCysFgwmJibx4z/+QOmEP8zoM9C+mJWRLxQJduxY6ro4qsmaP1birs3y7j0+Ponp6b1O+1yeYMA8Bh1SZA8NacJxTROzQFmcsRa0mDw00iQWfGHhwcLrVdKTDxh+BjHkZYgxVa8rBs90dnYXssIpXbZZXtTG4uK8tz5HH4MO4aJQjab25hCcMvvm5e6asvbM8hEYH58s/S0tzqxtCN1xsmumpnZi69Z9nbdZ0fzvq89RMOgQLgrViHlCi7luPqjSnqMC4/T03tLvxCDEh+w4WQebzrw+2qypg61LaErokKqpkftOKGaBPFamRB4HsIi1azeuqEOoddNGk/as8p1Y1PCxpeodxYXJp+s2KzN3+ehzXpwPReTtAK4AsATgfgBXGWOOicgkgA8BeA6A9xtjri34jV8DcC2ABQAfM8a8Mb3+JgC/iMRQ9DpjzN+XlYdRCcQWRU5MGp0hSTahOLP2Ha3OvHXn+bm5/Th8ePcqR1rXfU5VVIKInG2MeTj9/+sAPNMYc42InAng2QAuAXBJnmAgIi8EsAfAzxpjTorIk40x94vIMwH8OYDnAtgA4DYA08aY7OTcKV0KBiRufIc9EXtQiNePxsiLNkJl131OVVTCQChIORPpkzXGPArgDhH5vpKf+BUA7zDGnEy/d396/QoAf5Fe/5qIfAWJkHDAZvkJyaPMHhiSjbrvxK6GjwGNJp820Wda+pw350MRuV5EvgFgJ4C31Pz6NIDni8inROSTIvKc9PpTAXxj6HP3pdcI6YSyCSk0GzUhmtHozBuD46ozwUBEbhORuzJeVwCAMWaPMeYCAPuR+ArUYQLAkwD8KIDfBPBBEREAkvHZTFuJiFwtIgdF5OCJEydq3p6QbLImqgG+JyxCYkNj5EUM0WfOTAnGmBdV/OgtAD4G4K01fv4+AH9lEgeJT4vIEoB16fULhj73NADHcsp3E4CbgMTHoMa9Ccllpcd7flQCIcQOWtTvA2KIPvPiYyAiW4wxh9M/LwdwT82f+BsAPwlgRkSmAawB8ACAjwC4RUTeicT5cAuAT9spNSHV0DZREUK6I4aQZF95DN4hIhchCVe8F8A1gzdE5CiAswGsEZGXAnixMeZuEXkvgBuNMQcBvA/A+0TkLgCPAdiVag++JCIfBHA3kjDGXy2LSCCEEEJsEvrmgIcogeGKhBBC+kdeuCJTIhNCCCFkGQoGhBBCCFmGggEhhBBClqFgQAghhJBlKBgQQgghZBkKBoQQQghZhoIBIYQQQpahYEAIIYSQZSgYEEIIIWQZCgaEEEIIWYaCASGEEEKWoWBACCGEkGUoGBBCCCFkGQoGhBBCCFmGxy4DEJETAO61/LPrADxg+Td9EEs9ANZFK6yLTmKpSyz1AOzXZaMx5vzRixQMHCEiB7POuQ6NWOoBsC5aYV10EktdYqkH0F1daEoghBBCyDIUDAghhBCyDAUDd9zkuwCWiKUeAOuiFdZFJ7HUJZZ6AB3VhT4GhBBCCFmGGgNCCCGELEPBoAIiclREviginxORgyPvXSciRkTWFXx/XEQ+KyIfHbp2noj8g4gcTv99kss6DN3XRV3eJiLfTH/zcyJymcs6DN23cV3yvhvicymoS+fPpWU9zhWRD4nIPSIyKyLb0+shPpO8ugQ1VkTkoqGyfk5EHhaRX0/fC+q5lNQlqOeSvv96EfmSiNwlIn8uIk9Ir7d/LsYYvkpeAI4CWJdx/QIAf48kB8Kq94c+9wYAtwD46NC13wPwW+n/fwvA/x1wXd4G4LqQnkvBd4N7LgXf7fy5tKzHPgCvSf+/BsC5AT+TvLoEN1Z3D26OAAAEXklEQVSGPjsO4DiS2Pcgn0tBXYJ6LgCeCuBrAE5P//4ggKtsPRdqDNrxLgBvBJDrqCEiTwPwswDeO/LWFUgmD6T/vtRFAWvQpi7aKK1LAcE9l0AorIeInA3gBQD+KwAYYx4zxvxr+nZQz6SkLtqo079+CsBXjTGDZHBBPZcRRuuijSp1mQBwuohMADgDwLH0euvnQsGgGgbAx0XkThG5GgBE5HIA3zTGfL7ku3+E5AEvjVyfMsZ8CwDSf59sucx5uKgLAFwrIl8Qkfd1pVJEu7qs+m5KiM8lry5A98+laT02AzgB4M8kMVW9V0TOTN8L7ZkU1QUIb6wMeAWAPx/6O7TnMsxoXYCAnosx5psA/gDA1wF8C8BDxpiPp2+3fy5dq05CfAHYkP77ZACfR7Ib+BSAc0rUQf8ewA3p/3dgpfr9X0c+++2A6zKFRDU3BuB6AO/TXJe874b4XErq0vlzadG/tgFYAPAj6d97Abw9xGdSUpfgxkr6/hokqXinhq4F9VxK6hLUcwHwJAD/BOB8AKcB+BsAr7T1XKgxqIAx5lj67/0A/hrATwB4OoDPi8hRAE8D8BkRWT/y1ecBuDz9zF8A+EkR+UD63pyIPAUA0n/vd10PwE1djDFzxphFY8wSgD8F8Fzldcn67qDMoT2X3Lr4eC4t6nEfgPuMMZ9K//4QgB9K/x/aM8mtS4hjJeXfAfiMMWZu6Fpoz2XAqroE+FxeBOBrxpgTxpjHAfwVgB9L32v9XCgYlCAiZ4rIWYP/A3gxgP/PGPNkY8wmY8wmJBPBDxljjg9/1xjzJmPM09LPvALAPxljXpm+/REAu9L/7wLw4VDrMuiEKS8DcJfmuuR8d1DmoJ5LUV26fi4t+9dxAN8QkYvSSz8F4O70/0E9k6K6hDZWhvjfsVr1HtRzGWJVXQJ8Ll8H8KMicoaICJI+Npu+1/652FSLxPhCYi/8fPr6EoA9GZ85ilTdA2ADgFszPrMDK9XvkwD+EcDh9N/zAq7LzQC+COALaad8iua6FH03tOdSUpdOn0vb/gXgUgAH0/L+DYAnhfhMSuoS1FhJ/z4DwDxS9fbQ9RCfS15dQnwuvwPgnv+/vTtUySwI4zj8fzGoxWDZqmnBRdsGb2S71Wz3VoStdpOYFkzqegGCN6AYzK/hyOAFyDeCz5NOfOHA4ccMZyZLxPxNsv5Z78XJhwDAYCsBABiEAQAwCAMAYBAGAMAgDACAQRgAAIMwAKaohW8QfDHOMQBWpqp2klwkuUpymOQuyX6SzSTn3X06bTggyXJtI8Aq/Uxy1N3HVbXd3U9VtZbksqoOuvt+9oDwnVnGA1btsbuv35//VNVNktskv5LszRsLSKwYAKv3miRVtZvkJMnv7n6uqrMkGzMHA6wYAPNsZYmEl6r6keU6XGAyKwbAFN39v6pus9ws95Dk3+SRgPgrAQD4wFYCADAIAwBgEAYAwCAMAIBBGAAAgzAAAAZhAAAMwgAAGN4ABLpsnsgi/ioAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "bmags = pd.DataFrame(records, columns=['ra', 'dec', mag_col])\n",
    "\n",
    "plt.figure(figsize=(8, 8))\n",
    "plt.xlabel('ra')\n",
    "plt.ylabel('dec')\n",
    "plt.suptitle(pop + ' Box search', size='xx-large', y=0.92)\n",
    "p = plt.scatter(bmags['ra'], bmags['dec'], color='y')\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "desc-python",
   "language": "python",
   "name": "desc-python"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}