{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# SSDF Selection Functions\n",
    "## Depth maps and selection functions\n",
    "\n",
    "The simplest selection function available is the field MOC which specifies the area for which there is Herschel data. Each pristine catalogue also has a MOC defining the area for which that data is available.\n",
    "\n",
    "The next stage is to provide mean flux standard deviations which act as a proxy for the catalogue's 5$\\sigma$ depth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This notebook was run with herschelhelp_internal version: \n",
      "017bb1e (Mon Jun 18 14:58:59 2018 +0100)\n",
      "This notebook was executed on: \n",
      "2018-06-26 11:40:51.154813\n"
     ]
    }
   ],
   "source": [
    "from herschelhelp_internal import git_version\n",
    "print(\"This notebook was run with herschelhelp_internal version: \\n{}\".format(git_version()))\n",
    "import datetime\n",
    "print(\"This notebook was executed on: \\n{}\".format(datetime.datetime.now()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "#%config InlineBackend.figure_format = 'svg'\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rc('figure', figsize=(10, 6))\n",
    "\n",
    "import os\n",
    "import time\n",
    "\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "from astropy.table import Column, Table, join\n",
    "import numpy as np\n",
    "from pymoc import MOC\n",
    "import healpy as hp\n",
    "#import pandas as pd #Astropy has group_by function so apandas isn't required.\n",
    "import seaborn as sns\n",
    "\n",
    "import warnings\n",
    "#We ignore warnings - this is a little dangerous but a huge number of warnings are generated by empty cells later\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "from herschelhelp_internal.utils import inMoc, coords_to_hpidx, flux_to_mag\n",
    "from herschelhelp_internal.masterlist import find_last_ml_suffix, nb_ccplots\n",
    "\n",
    "from astropy.io.votable import parse_single_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "FIELD = 'SSDF'\n",
    "#FILTERS_DIR = \"/Users/rs548/GitHub/herschelhelp_python/database_builder/filters/\"\n",
    "FILTERS_DIR = \"/opt/herschelhelp_python/database_builder/filters/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Depth maps produced using: master_catalogue_ssdf_20180221.fits\n"
     ]
    }
   ],
   "source": [
    "TMP_DIR = os.environ.get('TMP_DIR', \"./data_tmp\")\n",
    "OUT_DIR = os.environ.get('OUT_DIR', \"./data\")\n",
    "SUFFIX = find_last_ml_suffix()\n",
    "#SUFFIX = \"20171016\"\n",
    "\n",
    "master_catalogue_filename = \"master_catalogue_{}_{}.fits\".format(FIELD.lower(), SUFFIX)\n",
    "master_catalogue = Table.read(\"{}/{}\".format(OUT_DIR, master_catalogue_filename))\n",
    "\n",
    "print(\"Depth maps produced using: {}\".format(master_catalogue_filename))\n",
    "\n",
    "ORDER = 10\n",
    "#TODO write code to decide on appropriate order\n",
    "\n",
    "field_moc = MOC(filename=\"../../dmu2/dmu2_field_coverages/{}_MOC.fits\".format(FIELD))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Remove sources whose signal to noise ratio is less than five as these will have been selected using forced \n",
    "# photometry and so the errors will not refelct the RMS of the map \n",
    "for n,col in enumerate(master_catalogue.colnames):\n",
    "    if col.startswith(\"f_\"):\n",
    "        err_col = \"ferr{}\".format(col[1:])\n",
    "        errs = master_catalogue[err_col]\n",
    "        fluxes = master_catalogue[col]\n",
    "        mask = fluxes/errs < 5.0\n",
    "        master_catalogue[col][mask] = np.nan\n",
    "        master_catalogue[err_col][mask] = np.nan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## I - Group masterlist objects by healpix cell and calculate depths\n",
    "We add a column to the masterlist catalogue for the target order healpix cell <i>per object</i>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#Add a column to the catalogue with the order=ORDER hp_idx\n",
    "master_catalogue.add_column(Column(data=coords_to_hpidx(master_catalogue['ra'],\n",
    "                                                       master_catalogue['dec'],\n",
    "                                                       ORDER), \n",
    "                                   name=\"hp_idx_O_{}\".format(str(ORDER))\n",
    "                                  )\n",
    "                           )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Convert catalogue to pandas and group by the order=ORDER pixel\n",
    "\n",
    "group = master_catalogue.group_by([\"hp_idx_O_{}\".format(str(ORDER))])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#Downgrade the groups from order=ORDER to order=13 and then fill out the appropriate cells\n",
    "#hp.pixelfunc.ud_grade([2599293, 2599294], nside_out=hp.order2nside(13))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## II Create a table of all Order=13 healpix cells in the field and populate it\n",
    "We create a table with every order=13 healpix cell in the field MOC. We then calculate the healpix cell at lower order that the order=13 cell is in. We then fill in the depth at every order=13 cell as calculated for the lower order cell that that the order=13 cell is inside."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "depths = Table()\n",
    "depths['hp_idx_O_13'] = list(field_moc.flattened(13))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=10</i>\n",
       "<table id=\"table140353671788584-629698\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>idx</th><th>hp_idx_O_13</th></tr></thead>\n",
       "<tr><td>0</td><td>759169024</td></tr>\n",
       "<tr><td>1</td><td>759169025</td></tr>\n",
       "<tr><td>2</td><td>759169026</td></tr>\n",
       "<tr><td>3</td><td>759169027</td></tr>\n",
       "<tr><td>4</td><td>759169028</td></tr>\n",
       "<tr><td>5</td><td>759169029</td></tr>\n",
       "<tr><td>6</td><td>759169030</td></tr>\n",
       "<tr><td>7</td><td>759169031</td></tr>\n",
       "<tr><td>8</td><td>759169032</td></tr>\n",
       "<tr><td>9</td><td>759169033</td></tr>\n",
       "</table><style>table.dataTable {clear: both; width: auto !important; margin: 0 !important;}\n",
       ".dataTables_info, .dataTables_length, .dataTables_filter, .dataTables_paginate{\n",
       "display: inline-block; margin-right: 1em; }\n",
       ".paginate_button { margin-right: 5px; }\n",
       "</style>\n",
       "<script>\n",
       "\n",
       "var astropy_sort_num = function(a, b) {\n",
       "    var a_num = parseFloat(a);\n",
       "    var b_num = parseFloat(b);\n",
       "\n",
       "    if (isNaN(a_num) && isNaN(b_num))\n",
       "        return ((a < b) ? -1 : ((a > b) ? 1 : 0));\n",
       "    else if (!isNaN(a_num) && !isNaN(b_num))\n",
       "        return ((a_num < b_num) ? -1 : ((a_num > b_num) ? 1 : 0));\n",
       "    else\n",
       "        return isNaN(a_num) ? -1 : 1;\n",
       "}\n",
       "\n",
       "require.config({paths: {\n",
       "    datatables: 'https://cdn.datatables.net/1.10.12/js/jquery.dataTables.min'\n",
       "}});\n",
       "require([\"datatables\"], function(){\n",
       "    console.log(\"$('#table140353671788584-629698').dataTable()\");\n",
       "    \n",
       "jQuery.extend( jQuery.fn.dataTableExt.oSort, {\n",
       "    \"optionalnum-asc\": astropy_sort_num,\n",
       "    \"optionalnum-desc\": function (a,b) { return -astropy_sort_num(a, b); }\n",
       "});\n",
       "\n",
       "    $('#table140353671788584-629698').dataTable({\n",
       "        order: [],\n",
       "        pageLength: 50,\n",
       "        lengthMenu: [[10, 25, 50, 100, 500, 1000, -1], [10, 25, 50, 100, 500, 1000, 'All']],\n",
       "        pagingType: \"full_numbers\",\n",
       "        columnDefs: [{targets: [0, 1], type: \"optionalnum\"}]\n",
       "    });\n",
       "});\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "depths[:10].show_in_notebook()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "depths.add_column(Column(data=hp.pixelfunc.ang2pix(2**ORDER,\n",
    "                     hp.pixelfunc.pix2ang(2**13, depths['hp_idx_O_13'], nest=True)[0],\n",
    "                     hp.pixelfunc.pix2ang(2**13, depths['hp_idx_O_13'], nest=True)[1],\n",
    "                     nest = True),\n",
    "                     name=\"hp_idx_O_{}\".format(str(ORDER))\n",
    "                        )\n",
    "                 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=10</i>\n",
       "<table id=\"table140353671788696-980795\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>idx</th><th>hp_idx_O_13</th><th>hp_idx_O_10</th></tr></thead>\n",
       "<tr><td>0</td><td>759169024</td><td>11862016</td></tr>\n",
       "<tr><td>1</td><td>759169025</td><td>11862016</td></tr>\n",
       "<tr><td>2</td><td>759169026</td><td>11862016</td></tr>\n",
       "<tr><td>3</td><td>759169027</td><td>11862016</td></tr>\n",
       "<tr><td>4</td><td>759169028</td><td>11862016</td></tr>\n",
       "<tr><td>5</td><td>759169029</td><td>11862016</td></tr>\n",
       "<tr><td>6</td><td>759169030</td><td>11862016</td></tr>\n",
       "<tr><td>7</td><td>759169031</td><td>11862016</td></tr>\n",
       "<tr><td>8</td><td>759169032</td><td>11862016</td></tr>\n",
       "<tr><td>9</td><td>759169033</td><td>11862016</td></tr>\n",
       "</table><style>table.dataTable {clear: both; width: auto !important; margin: 0 !important;}\n",
       ".dataTables_info, .dataTables_length, .dataTables_filter, .dataTables_paginate{\n",
       "display: inline-block; margin-right: 1em; }\n",
       ".paginate_button { margin-right: 5px; }\n",
       "</style>\n",
       "<script>\n",
       "\n",
       "var astropy_sort_num = function(a, b) {\n",
       "    var a_num = parseFloat(a);\n",
       "    var b_num = parseFloat(b);\n",
       "\n",
       "    if (isNaN(a_num) && isNaN(b_num))\n",
       "        return ((a < b) ? -1 : ((a > b) ? 1 : 0));\n",
       "    else if (!isNaN(a_num) && !isNaN(b_num))\n",
       "        return ((a_num < b_num) ? -1 : ((a_num > b_num) ? 1 : 0));\n",
       "    else\n",
       "        return isNaN(a_num) ? -1 : 1;\n",
       "}\n",
       "\n",
       "require.config({paths: {\n",
       "    datatables: 'https://cdn.datatables.net/1.10.12/js/jquery.dataTables.min'\n",
       "}});\n",
       "require([\"datatables\"], function(){\n",
       "    console.log(\"$('#table140353671788696-980795').dataTable()\");\n",
       "    \n",
       "jQuery.extend( jQuery.fn.dataTableExt.oSort, {\n",
       "    \"optionalnum-asc\": astropy_sort_num,\n",
       "    \"optionalnum-desc\": function (a,b) { return -astropy_sort_num(a, b); }\n",
       "});\n",
       "\n",
       "    $('#table140353671788696-980795').dataTable({\n",
       "        order: [],\n",
       "        pageLength: 50,\n",
       "        lengthMenu: [[10, 25, 50, 100, 500, 1000, -1], [10, 25, 50, 100, 500, 1000, 'All']],\n",
       "        pagingType: \"full_numbers\",\n",
       "        columnDefs: [{targets: [0, 1, 2], type: \"optionalnum\"}]\n",
       "    });\n",
       "});\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "depths[:10].show_in_notebook()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=10</i>\n",
       "<table id=\"table140353675322256-686087\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>idx</th><th>hp_idx_O_13</th><th>hp_idx_O_10</th><th>ferr_ap_vista_j_mean</th><th>f_ap_vista_j_p90</th><th>ferr_vista_j_mean</th><th>f_vista_j_p90</th><th>ferr_ap_vista_h_mean</th><th>f_ap_vista_h_p90</th><th>ferr_vista_h_mean</th><th>f_vista_h_p90</th><th>ferr_ap_vista_ks_mean</th><th>f_ap_vista_ks_p90</th><th>ferr_vista_ks_mean</th><th>f_vista_ks_p90</th><th>ferr_ap_irac_i1_mean</th><th>f_ap_irac_i1_p90</th><th>ferr_irac_i1_mean</th><th>f_irac_i1_p90</th><th>ferr_ap_irac_i2_mean</th><th>f_ap_irac_i2_p90</th><th>ferr_irac_i2_mean</th><th>f_irac_i2_p90</th><th>ferr_ap_decam_g_mean</th><th>f_ap_decam_g_p90</th><th>ferr_decam_g_mean</th><th>f_decam_g_p90</th><th>ferr_ap_decam_r_mean</th><th>f_ap_decam_r_p90</th><th>ferr_decam_r_mean</th><th>f_decam_r_p90</th><th>ferr_ap_decam_i_mean</th><th>f_ap_decam_i_p90</th><th>ferr_decam_i_mean</th><th>f_decam_i_p90</th><th>ferr_ap_decam_z_mean</th><th>f_ap_decam_z_p90</th><th>ferr_decam_z_mean</th><th>f_decam_z_p90</th><th>ferr_ap_decam_y_mean</th><th>f_ap_decam_y_p90</th><th>ferr_decam_y_mean</th><th>f_decam_y_p90</th></tr></thead>\n",
       "<thead><tr><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th>uJy</th><th>uJy</th><th>uJy</th><th>uJy</th><th>uJy</th><th>uJy</th><th>uJy</th><th>uJy</th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th></tr></thead>\n",
       "<tr><td>0</td><td>571095615</td><td>8923368</td><td>2.3876948</td><td>14.246212005615234</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>0.14885482325566032</td><td>0.788024718842322</td><td>nan</td><td>nan</td><td>0.18116324810829695</td><td>3.170422339039513</td><td>0.21496994625413532</td><td>2.6941524288751073</td><td>0.32797459582511074</td><td>6.148384335767076</td><td>0.41297766829611154</td><td>6.059855671774965</td><td>0.561371343824901</td><td>8.726608154435011</td><td>0.7914561385330093</td><td>8.359950365213846</td><td>1.7160408967689051</td><td>11.46283892678945</td><td>1.7980520934595297</td><td>9.022755789438893</td></tr>\n",
       "<tr><td>1</td><td>571095658</td><td>8923369</td><td>2.985487</td><td>937.9022079467774</td><td>7.291202</td><td>929.9276901245117</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>6.3473563</td><td>836.1466239929199</td><td>14.499721</td><td>1034.2193603515625</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>0.15823875757826117</td><td>51.766714928481214</td><td>0.22606746166863143</td><td>44.695235724724554</td><td>0.18602492681917504</td><td>151.68227198124202</td><td>0.2731085492503393</td><td>132.60886512386293</td><td>0.3118835430831495</td><td>339.4680307424261</td><td>0.4969930345142875</td><td>331.66935920958343</td><td>0.5741958423201232</td><td>592.108303665983</td><td>1.0988309505461755</td><td>567.9553119454632</td><td>1.7439328523741413</td><td>779.379656337846</td><td>3.0913253706111683</td><td>651.4730653688407</td></tr>\n",
       "<tr><td>2</td><td>571095659</td><td>8923369</td><td>2.985487</td><td>937.9022079467774</td><td>7.291202</td><td>929.9276901245117</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>6.3473563</td><td>836.1466239929199</td><td>14.499721</td><td>1034.2193603515625</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>0.15823875757826117</td><td>51.766714928481214</td><td>0.22606746166863143</td><td>44.695235724724554</td><td>0.18602492681917504</td><td>151.68227198124202</td><td>0.2731085492503393</td><td>132.60886512386293</td><td>0.3118835430831495</td><td>339.4680307424261</td><td>0.4969930345142875</td><td>331.66935920958343</td><td>0.5741958423201232</td><td>592.108303665983</td><td>1.0988309505461755</td><td>567.9553119454632</td><td>1.7439328523741413</td><td>779.379656337846</td><td>3.0913253706111683</td><td>651.4730653688407</td></tr>\n",
       "<tr><td>3</td><td>571095703</td><td>8923370</td><td>2.5948894</td><td>311.5769348144531</td><td>6.784351</td><td>754.6826354980471</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>6.1474032</td><td>404.22176513671866</td><td>15.348894</td><td>644.5146026611328</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>0.12715930290395008</td><td>14.891955977165672</td><td>0.16338834901202023</td><td>25.60270616639062</td><td>0.14531298717198027</td><td>31.14947795692911</td><td>0.21102354498619774</td><td>46.36672924628673</td><td>0.23867347272856043</td><td>68.17725072626247</td><td>0.38003366244928194</td><td>97.83919024826011</td><td>0.4712349195624097</td><td>152.5575980504968</td><td>0.8607349629731965</td><td>195.31063208816298</td><td>1.3646868871369906</td><td>277.83456226325575</td><td>2.323780044363607</td><td>316.8299991419369</td></tr>\n",
       "<tr><td>4</td><td>571095730</td><td>8923370</td><td>2.5948894</td><td>311.5769348144531</td><td>6.784351</td><td>754.6826354980471</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>6.1474032</td><td>404.22176513671866</td><td>15.348894</td><td>644.5146026611328</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>0.12715930290395008</td><td>14.891955977165672</td><td>0.16338834901202023</td><td>25.60270616639062</td><td>0.14531298717198027</td><td>31.14947795692911</td><td>0.21102354498619774</td><td>46.36672924628673</td><td>0.23867347272856043</td><td>68.17725072626247</td><td>0.38003366244928194</td><td>97.83919024826011</td><td>0.4712349195624097</td><td>152.5575980504968</td><td>0.8607349629731965</td><td>195.31063208816298</td><td>1.3646868871369906</td><td>277.83456226325575</td><td>2.323780044363607</td><td>316.8299991419369</td></tr>\n",
       "<tr><td>5</td><td>571095738</td><td>8923370</td><td>2.5948894</td><td>311.5769348144531</td><td>6.784351</td><td>754.6826354980471</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>6.1474032</td><td>404.22176513671866</td><td>15.348894</td><td>644.5146026611328</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>0.12715930290395008</td><td>14.891955977165672</td><td>0.16338834901202023</td><td>25.60270616639062</td><td>0.14531298717198027</td><td>31.14947795692911</td><td>0.21102354498619774</td><td>46.36672924628673</td><td>0.23867347272856043</td><td>68.17725072626247</td><td>0.38003366244928194</td><td>97.83919024826011</td><td>0.4712349195624097</td><td>152.5575980504968</td><td>0.8607349629731965</td><td>195.31063208816298</td><td>1.3646868871369906</td><td>277.83456226325575</td><td>2.323780044363607</td><td>316.8299991419369</td></tr>\n",
       "<tr><td>6</td><td>571095707</td><td>8923370</td><td>2.5948894</td><td>311.5769348144531</td><td>6.784351</td><td>754.6826354980471</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>6.1474032</td><td>404.22176513671866</td><td>15.348894</td><td>644.5146026611328</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>0.12715930290395008</td><td>14.891955977165672</td><td>0.16338834901202023</td><td>25.60270616639062</td><td>0.14531298717198027</td><td>31.14947795692911</td><td>0.21102354498619774</td><td>46.36672924628673</td><td>0.23867347272856043</td><td>68.17725072626247</td><td>0.38003366244928194</td><td>97.83919024826011</td><td>0.4712349195624097</td><td>152.5575980504968</td><td>0.8607349629731965</td><td>195.31063208816298</td><td>1.3646868871369906</td><td>277.83456226325575</td><td>2.323780044363607</td><td>316.8299991419369</td></tr>\n",
       "<tr><td>7</td><td>571095710</td><td>8923370</td><td>2.5948894</td><td>311.5769348144531</td><td>6.784351</td><td>754.6826354980471</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>6.1474032</td><td>404.22176513671866</td><td>15.348894</td><td>644.5146026611328</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>0.12715930290395008</td><td>14.891955977165672</td><td>0.16338834901202023</td><td>25.60270616639062</td><td>0.14531298717198027</td><td>31.14947795692911</td><td>0.21102354498619774</td><td>46.36672924628673</td><td>0.23867347272856043</td><td>68.17725072626247</td><td>0.38003366244928194</td><td>97.83919024826011</td><td>0.4712349195624097</td><td>152.5575980504968</td><td>0.8607349629731965</td><td>195.31063208816298</td><td>1.3646868871369906</td><td>277.83456226325575</td><td>2.323780044363607</td><td>316.8299991419369</td></tr>\n",
       "<tr><td>8</td><td>571095718</td><td>8923370</td><td>2.5948894</td><td>311.5769348144531</td><td>6.784351</td><td>754.6826354980471</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>6.1474032</td><td>404.22176513671866</td><td>15.348894</td><td>644.5146026611328</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>0.12715930290395008</td><td>14.891955977165672</td><td>0.16338834901202023</td><td>25.60270616639062</td><td>0.14531298717198027</td><td>31.14947795692911</td><td>0.21102354498619774</td><td>46.36672924628673</td><td>0.23867347272856043</td><td>68.17725072626247</td><td>0.38003366244928194</td><td>97.83919024826011</td><td>0.4712349195624097</td><td>152.5575980504968</td><td>0.8607349629731965</td><td>195.31063208816298</td><td>1.3646868871369906</td><td>277.83456226325575</td><td>2.323780044363607</td><td>316.8299991419369</td></tr>\n",
       "<tr><td>9</td><td>571095709</td><td>8923370</td><td>2.5948894</td><td>311.5769348144531</td><td>6.784351</td><td>754.6826354980471</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>6.1474032</td><td>404.22176513671866</td><td>15.348894</td><td>644.5146026611328</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>0.12715930290395008</td><td>14.891955977165672</td><td>0.16338834901202023</td><td>25.60270616639062</td><td>0.14531298717198027</td><td>31.14947795692911</td><td>0.21102354498619774</td><td>46.36672924628673</td><td>0.23867347272856043</td><td>68.17725072626247</td><td>0.38003366244928194</td><td>97.83919024826011</td><td>0.4712349195624097</td><td>152.5575980504968</td><td>0.8607349629731965</td><td>195.31063208816298</td><td>1.3646868871369906</td><td>277.83456226325575</td><td>2.323780044363607</td><td>316.8299991419369</td></tr>\n",
       "</table><style>table.dataTable {clear: both; width: auto !important; margin: 0 !important;}\n",
       ".dataTables_info, .dataTables_length, .dataTables_filter, .dataTables_paginate{\n",
       "display: inline-block; margin-right: 1em; }\n",
       ".paginate_button { margin-right: 5px; }\n",
       "</style>\n",
       "<script>\n",
       "\n",
       "var astropy_sort_num = function(a, b) {\n",
       "    var a_num = parseFloat(a);\n",
       "    var b_num = parseFloat(b);\n",
       "\n",
       "    if (isNaN(a_num) && isNaN(b_num))\n",
       "        return ((a < b) ? -1 : ((a > b) ? 1 : 0));\n",
       "    else if (!isNaN(a_num) && !isNaN(b_num))\n",
       "        return ((a_num < b_num) ? -1 : ((a_num > b_num) ? 1 : 0));\n",
       "    else\n",
       "        return isNaN(a_num) ? -1 : 1;\n",
       "}\n",
       "\n",
       "require.config({paths: {\n",
       "    datatables: 'https://cdn.datatables.net/1.10.12/js/jquery.dataTables.min'\n",
       "}});\n",
       "require([\"datatables\"], function(){\n",
       "    console.log(\"$('#table140353675322256-686087').dataTable()\");\n",
       "    \n",
       "jQuery.extend( jQuery.fn.dataTableExt.oSort, {\n",
       "    \"optionalnum-asc\": astropy_sort_num,\n",
       "    \"optionalnum-desc\": function (a,b) { return -astropy_sort_num(a, b); }\n",
       "});\n",
       "\n",
       "    $('#table140353675322256-686087').dataTable({\n",
       "        order: [],\n",
       "        pageLength: 50,\n",
       "        lengthMenu: [[10, 25, 50, 100, 500, 1000, -1], [10, 25, 50, 100, 500, 1000, 'All']],\n",
       "        pagingType: \"full_numbers\",\n",
       "        columnDefs: [{targets: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42], type: \"optionalnum\"}]\n",
       "    });\n",
       "});\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for col in master_catalogue.colnames:\n",
    "    if col.startswith(\"f_\"):\n",
    "        errcol = \"ferr{}\".format(col[1:])\n",
    "        depths = join(depths, \n",
    "                      group[\"hp_idx_O_{}\".format(str(ORDER)), errcol].groups.aggregate(np.nanmean),\n",
    "                     join_type='left')\n",
    "        depths[errcol].name = errcol + \"_mean\"\n",
    "        depths = join(depths, \n",
    "                      group[\"hp_idx_O_{}\".format(str(ORDER)), col].groups.aggregate(lambda x: np.nanpercentile(x, 90.)),\n",
    "                     join_type='left')\n",
    "        depths[col].name = col + \"_p90\"\n",
    "\n",
    "depths[:10].show_in_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## III - Save the depth map table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "depths.write(\"{}/depths_{}_{}.fits\".format(OUT_DIR, FIELD.lower(), SUFFIX), overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## IV - Overview plots\n",
    "\n",
    "### IV.a - Filters\n",
    "First we simply plot all the filters available on this field to give an overview of coverage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'decam_g',\n",
       " 'decam_i',\n",
       " 'decam_r',\n",
       " 'decam_y',\n",
       " 'decam_z',\n",
       " 'irac_i1',\n",
       " 'irac_i2',\n",
       " 'vista_h',\n",
       " 'vista_j',\n",
       " 'vista_ks'}"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tot_bands = [column[2:] for column in master_catalogue.colnames \n",
    "             if (column.startswith('f_') & ~column.startswith('f_ap_'))]\n",
    "ap_bands = [column[5:] for column in master_catalogue.colnames \n",
    "            if column.startswith('f_ap_') ]\n",
    "bands = set(tot_bands) | set(ap_bands)\n",
    "bands"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Passbands on SSDF')"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAEgCAYAAACHCxsiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXl8FPX5xz/P7G421+YOSYAEAiSEcMkhHoiKVasWFItXvbVWwVq1qMVSS63+LOLZKmrxLoInqCB41IujIrSBcJgQ7hBCLsidbI7dnef3x+xudpNNsjszmwTyfb9eebE7M99jQ3aeeW5iZggEAoFAIOg9pN7egEAgEAgE/R0hjAUCgUAg6GWEMBYIBAKBoJcRwlggEAgEgl5GCGOBQCAQCHoZIYwFAoFAIOhlhDAWnHIQ0aNEtPxUW0sgEJy6CGEs6BGIqJCImoiogYjKiegtIors7X2dChDRr4mogIjqnb/bdURkcZ4bTESriOgEEdUS0W4iutV5bigRsfP/xPX/spaILmo3v+f/netnYC98VIHglEUIY0FPMpOZIwFMBHA6gEd6eT8nPUR0HoC/AfgVM1sAjALwoccl7wA4CmAIgHgANwMobzdNjPP/ZTyArwF84hLYHsxk5kiPnxL9P41A0H8RwljQ4zDzMQBfABgDAER0GxHtcWp2h4joLte1RJTg1NZqiKiKiDYRkeQ8N5+IjjnH7SWin3ksE0pEHzjPbSei8R5zPkxEB53n8onoSo9ztxLRf4joGSKqJqLDRHSpx/l0ItrgHPs1gASPc6FEtJyIKp37/R8RJfn6HRDRKCJa77wuj4gu9zj3NhG95NRw64loKxEN7+TXeTqAH5k51/m7rWLmfzFzvcf5t5m5kZntzJzLzF908v9Sxsz/APAogMWu37NAIAg+4ssm6HGIKBXAZQBynYcqAMwAEAXgNgDPE9FE57kHABQDSASQBGABACaikQDuAXC6UyP8OYBCj2WuAPARgDgA7wL4lIhMznMHAUwDEA3grwCWE1GKx9gzAOyFImifAvAGEZHz3LsAtjnPPQ7gFo9xtzjnTIWihc4B0OTj85sAfAbg3wAGAPgdgBXOz+TiV869xQI4AOCJ9vM42Qrg50T0VyKaSkTmdue3AHiJiK4jorRO5mjPx859jezuQoFAoA9CGAt6kk+JqAbAfwBsgGJeBTOvY+aDrLABipCa5hxjA5ACYAgz25h5EysF1R0AzACyicjEzIXMfNBjrW3MvJKZbQCeAxAK4Ezneh8xcwkzy8z8AYD9AKZ4jD3CzK8xswPAv5zrJzmF2ekA/szMLcy8EYpQdWGDIoRHMLODmbcxc52P38OZACIBPMnMrcz8HYC1UASwi4+Z+b/MbAewAsBpvn6hzLwJwC+hmP7XAagkoueIyOC85GoAmwD8GcBhItpBRKf7mssDlwk6zuPYp04tvoaIPu1mvEAgCBAhjAU9ySxmjmHmIcx8NzM3AQARXUpEW5xm6BooWrPL/Ps0FM3w304T9sMAwMwHANwPxaRaQUTvtwsqOup6wcwyFO16oHO9m51Cqca53hiP9QCgzGOs1fky0jm+mpkbPa494vH6HQBfAXifiEqI6CkPbdyTgQCOOvflOc8gX3sAYHWu7xNm/oKZZ0IRnlcAuBXAHc5z1cz8MDOPhmJZ2AFFsFJn83nso8rjmOv/LoaZZ3UxViAQqEAIY0Gv4jSrrgLwDIAkZo4B8DkAAgBmrmfmB5h5GICZAOa5fMPM/C4znwMlOIkBLPaYOtVjDQnAYAAlRDQEwGtQTNzxzvV+cq3XDaUAYokowuOY2/Tr1Nz/yszZAM6GYnq/2cc8JQBS2/lk0wAc82MPneLU9L8F8B2c/vh2509A+T0PhLfW254robgO9mrZj0Ag8B8hjAW9TQgUc/NxAHZnsNTFrpNENIOIRjg1uToo5mkHEY0kogucwrwZim/W4THvJCL6JREZoWjQLVD8pxFQBPdx5/y3wYfg8gUzHwGQA+CvRBRCROdAeUBw7XU6EY11mojroJitHT6m2gqgEcAfiMhEROc753nfn314QkRXOP3BsaQwBcB5zs8KIlpMRGOIyEhKutNcAAeYudLHXElEdA+AvwD4YzvNXSAQBBEhjAW9ijPq914o6TjVAK4HsMbjkgwA3wBoAPAjgJeZeT0UAf4kgBNQTLoDoAR3uVgN4FrnnDcB+KVTc80H8KxzrnIAYwH8EMCWr4cS4FUFRWgt8ziXDGAlFEG8B4pfvENBEGZuBXA5gEud+38ZwM3MXBDAPlxUA/gNFL93nXO9p5l5hfN8OIBPANQAOATFinB5uzlqiKgRwG4oLoKrmflNFXsRCAQqISUWRiAQCAQCQW8hNGOBQCAQCHoZIYwFAoFAIOhlgiaMiehNIqogop86OU9E9AIRHSCiXR5FHgQCgUAg6FcEUzN+G8AlXZy/FEpwTgaAOwG8EsS9CAQCgUDQZwmaMHZWJ6rq4pIrACxzVl3aAiCmXUlCgUAgEAj6BcZeXHsQPKokQamQNAhKYQUviOhOKNozIiIiJmVlZfXIBgUCgeBUYdu2bSeYOVHjHAOMRuPrUHLzRcxRYMgAfrLb7XdMmjSpov3J3hTGvioe+cyzYuZXAbwKAJMnT+acnJxg7ksgEAhOOYjoSPdXdY3RaHw9OTl5VGJiYrUkSSIvNgBkWabjx49nl5WVvY6Ouf69+mRTDI+ShXCWK+ylvQgEAoGge8YkJibWCUEcOJIkcWJiYi06qfjXm8J4DYCbnVHVZwKoZeYOJmqBQCAQ9BkkIYjV4/zd+ZS7QTNTE9F7AM4HkEBExVBKB5oAgJn/CaUZwGVQOvJYofSxFQgEAoGg3xE0YczMv+rmPAP4bbDWFwgEAsGpzbx58wZGRkY6HnvssfLe3otWRDScQCAQCAS9TG9GUwsEAoHgJOWhlTtT95XVh+s5Z2ayxfr0VeOPdnXN/Pnzkz/44IOEgQMHtsbHx9smTJhgzcvLM8+ZMyetqqrKGBoaKr/++utHJkyY0Hz06FHj7bffPqSoqMgMAEuWLDly0UUXNV544YXDS0tLQ1paWqQ5c+aUP/jggycAIDw8fMItt9xSsXHjxqjo6GjHE088UTx//vzUkpKSkMWLFxfdcMMNtb72NGnSpJEvvvhi0dlnn90EABMnTsx65ZVXjpxxxhlN/n52oRkLBH0Q2SGj8lhDb29DIOhTbNq0KfyTTz6J2717d/7atWsP7Ny5MwIA7rjjjiEvv/xyUV5e3p6nn366eO7cuWkAMGfOnLRp06bV7927Nz8vLy9/4sSJzQCwYsWKwry8vD07duzIX7p0aVJZWZkBAJqamqTp06fX5+Xl7YmIiHA88sgjgzZt2rTvo48+OvD4448P6mxft95664nXX389AQB27dplbm1tpUAEMSA0Y4GgT7Lpw/34acMx9/u5L50PySCenQV9h+402GDw/fffR1522WU1FotFBoCLL764prm5WcrNzY28+uqrh7uua21tJQDYvHmzZeXKlYcBwGg0Ij4+3gEAixcvTlq3bl0MAJSVlZny8vJCk5OTG00mE1911VV1ADB69Ogms9ksm81mnjJlStOxY8dCOtvXrbfeWv3000+ntLS0FP/zn/9MuP76608E+tmEMBYI+iCeghgAXvntevz2nxf00m4Egr4DkXe9KFmWYbFY7AUFBfn+jF+7dq1lw4YNlpycnAKLxSJPmTJlZFNTkwQARqORJUl56JUkCWazmQHAYDDA4XD4KlQFALBYLPK0adPq3n333Zg1a9bEbdu2za+9eCIetQWCPkbpQZ9uKVSVNvbI+myXUfV+Aay7j/fIegKBv1xwwQUN69ati2loaKDq6mrp66+/jgkPD5cHDx7c+uabb8YCinD+8ccfwwBg6tSp9U8//XQiANjtdlRVVUk1NTWG6Ohoh8VikXNzc0Ndpm6tzJkz58T8+fNTx48f35iUlOQIdLwQxgJBH+Pjp7e5X589e4T79Xt/3doj6x975AdYdxxH1YoCsIPhqGuBbLX1yNoCQVecc8451iuvvLJqzJgxo2fMmDF8ypQpDQDw3nvvHXrrrbcSRo4cmZ2RkTF61apVMQDwyiuvFG3YsMGSmZmZPWbMmOzt27eHzZ49u9Zut1NmZmb2ggULBo4fP16Xp9xp06ZZIyIiHLfddlvAJmoAICXd9+RB1KYWnOq8NOc7AMDFd4xGxuQkrHt5Fwp3Kd9vvU3VbJfRvL8add8UwRBtRtw1mSh59Eef1w5+cpquawt6FiLaxsyTtcyxc+fOwvHjx6sSNqc6hYWFpvPPP3/kwYMHfzIYDJ1et3PnzoTx48cPbX9caMYCQQ+zd0spXprzHY4WdOww2tJkd7/OmJwEALh0ztig7aX8hVxU/isftmMNaM6v7FQQCwSCzlmyZEn8mWeeOWrhwoXHuhLEXSECuASCHuabt/cAANb8fQcA4O5XpruDUl7//cYO10tSp3EjmnA0tMJeYfX7epYZFKS9CAQnA6tWrYr605/+NNjzWGpqasvXX3998J577qnUMrcQxgJBL7NswWbcsmgqPF1G46YP9nltQ3UzImNDdVm3/vvOM1MS7x6P4y/v9DrGrQ5QqLhlCPovs2fPrps9e3bAkdL+IMzUAkEP4itGo6G6BcyMl+d+7z52zjUZXtecdlEaAKCyRL+I6oYfOu9Yak6LQsojZ3gd41ZZt7UFAoE34jFXIOgBcr4oRGpWHKIHhAEAzrh8GMacNwhvPLAJALDmHzu8rm+fS5k2Kg47vi7C1tWHMGR0vOb9sOz9UJD88BRIYQbIjXYY4xTN2xAZgsFPTkPj9nJUf7gP3BpwtoZAIPAToRkLBEGmoboZW1cfwsrFOTi4vQKAEqgVGmFCyohoAEBxQbX7+rtfmd5hjkEjYwAA9ZXNuuypdt0h9+uUBVNgjDFDMhvdgtgTKUQJSJFtQjMWCIKFEMYCQZApL6xzv7Y7BVrKMEUIXzlvote102/K6qAVA3CXwmxu1Cff19NEbYgyd3ktmZS1hWYsEAQPIYwFgiBjDje5X8t2xTwcN1Ap+kMS4bK7xwEALrwtG9lTB3Y6T/r4hCDusnPIqRmzTQhjQd9i3rx5AxcuXJjU2/vQA+EzFgi6odHugESEMJWNGhwe5t3mxlYAQGRsmzaaPi4hoGIeeqYYpfzpjG6vcQtjEcAlEAQNIYwFgm4Yvmk3BppN2H72aFXj7R4apcPGMIUaYAwJvDDAoMxYHN55wu1vVovDw9RtsHTaiMaN20wtfMYCTz79bSoq8nXtZ4wB2VbMeumk62dcX18vXXvttUMPHDgQmpGR0Xz06NGQJUuWFJ177rl+J/ILM7VA4AclLep9tZ6acWuzXZUgBoAwiyKAm+pbVe8FAEof3xLQ9W4tXD65SucKTj36aj/jp59+OjEmJsaxb9++/EcffbQkPz8/4OYTQjMWCAA4mLG3sRnZkWFex+vsbVrtkaYWDAnrOtjJ59z2NiHWWNuCyJjA5wCAsEhFi22qb0Vssi6NZvzDKYzZIYSxwINuNNhg0Ff7GW/evDnyvvvuqwCA008/vTkzM9P/0nZOhDAWCAC8XFSBJw6VYsmoNFyVHOc+vmBfsfv1N5V1+PXgxIDnlh1tmnFNRZNqYRwerdwLGmu0acYuzM6I7m4xCM1Y0Hfoi/2M9Wi4JMzUgn6PzIx3SpSyst9V1Xudq/bw98ab1D27yh4aZd3xJtRVNqmaxxyurN/abO/mys4pfniT+3XC7WP8GuMyU7cvFCIQ9DR9tZ/x2Wef3fD+++/HAsC2bdtC9+3bF9bdmPYIYSzo9+yqb0JRs6Jtbqv1Ljc5PqrtO6VWFMntzLuJqRZV87h8zXaVUc22496WMzL6+fV3+YwdIoBL0Lv01X7GDz300PHKykpjZmZm9hNPPJE8cuTIptjY2IByAYWZWtDv+bGmAQBwY0o8VpRWwiYzTE4BVOuhGVtVCqP2wjjZX/NwO0xmRRjbWtTl+1q3VagaRwaXZqxquECgK4sXLy5bvHhxWfvjmzZt2t/+WGpqqv3bb7892P74xo0bO1wLAFarNdf1+rnnnivp7Fx7wsPD5Y8//vhweHg45+XlmS+++OLMjIyMgPxJQhgL+j1fV9ZhRLgZp0WFY3lpJcpbbRgcqvhna+0OWAwS6h2yemEse49Tm5ZkMEqQJIJNZSWs+vUq421ENLVA0CX19fXStGnTRtpsNmJmPP/880dCQ0MD+sIIYSzo1xy0NmNzTQP+mJ6CFLMiJEtb2oTxPmszhoaZsbuhSTfNWDKoL9hhNBtgV6kZe5L04GS/r23zGQvVWNC/6aqf8U8//bRHy9xCGAv6NStKqmAg4LqUOBxvVXKJyz1yigubWnB1UhzyGprQpFIYtRfGWupLm8wG1Zqxi8hzB8GUEEB8idtnLDRjQf8mmP2MhTAW9FtqbXYsKzmBSxOikWQ2uWVOuVMo22RGnV1GfIgR4QZJg2YsQzKSuy519jmd15/uDpMOmnH0z4cGdD1JBJCIphYIgomIphb0W9Ydr0WDQ8Zv05Q68/EmIwwEVLQqqUPVNuXfOJM2YeywMwwGCbHJSuXAkFD1z8DGEEl1AJcLUlNjWyLhMxYIgojQjAX9lq21jUgMMeI0i2KylYiQaDK5zdSVTmEcbzIiTJJg1WCmlgyE6x89U/Oe9TBTq4EkEpqxQBBEhDAW9FsOWVuQER7qVdFngNnoNlNXujVjg1MzVicEZZk1BW15YgoxoNmqruiHMSkcpsSAaxEoSCR8xgJBEBFmakG/5WBTM4aHe5emTAoxocIpjA9YWwC0mambVAojZu5Qwk8tRrMBdpWaMdtl/wt9tIMMQjMW9D1OpX7GQRXGRHQJEe0logNE9LCP82lE9D0R5RLRLiK6LJj7EQhcVNvsqLI5kB7WURiXtyia5welVQCAASEmhEmS6mhqlhk6yWLFTK3SZ8w2GVApjIXPWCAILkEzUxORAcBLAC4CUAzgf0S0hpk9w8IfAfAhM79CRNkAPgcwNFh7EghcuFoipoZ6N2IZYDai0maHXWZEGZWKV/EhRpgkQqPKfr4sc1sbQo2YQtRrxtCiGUskujYJvPjzD39OPVB9QNd+xiNiR1gfn/r4SdfP+Nprrx3iqnFdXl5uuv322yueffbZ0kA+ezA14ykADjDzIWZuBfA+gCvaXcMAopyvowGUQCDoAapa24KzPBkQYgIDOGGzo1mWMTUmEgBgJIJdZWcW5o6dZtRiCJFgU1mbmu2sWhjDIDRjQe/TV/sZf/DBB0cKCgry16xZcyAmJsZ+1113VQb62YIZwDUIgOcTTjGAM9pd8yiAfxPR7wBEALjQ10REdCeAOwEgLS1N940K+h/uSOkQ769AovP9iVYbqmx2jIwIBQCYiGBTK4xlBun02GswSF4tGQPah0MGGdU9FIhoakF7utNgg0Ff7WcMAFarlWbPnj38+eefL8rMzAy4z2kwhbGvb337b/OvALzNzM8S0VkA3iGiMczeJemZ+VUArwLA5MmTxR1BoBnPtCVPEpzvT9jsqLTZEed8b5Q0asY6maldxUMCDQpjZiUaWk2OMSCiqQV9hr7YzxgAbrrppiEzZ86snjVrVn1X13VGMM3UxQBSPd4PRkcz9K8BfAgAzPwjgFAACUHck0AAQBHGBCDWZPA6nhii1KeuaLWjxuZwC2tNZmpZv2hqg1OYyoFqqU5BqlozFtHUgj5AX+1nvGjRosSGhgbD3/72tw7dpPwlmML4fwAyiCidiEIAXAdgTbtrigD8DACIaBQUYXw8iHsSCAAAla12xJoMMLQTki4z9b7GZshAm2ZMSnlMNegZwCU5hamrtKbfe3CatlVV3wIAEpqxoPfpq/2MlyxZkrx3796wrKys7KysrOynnnoqMdA5gmamZmY7Ed0D4CsABgBvMnMeET0GIIeZ1wB4AMBrRPR7KCbsW5lVqh8CQQBUeWi9nliMBkQZJbznTGuKc2rOJo0BXJKOPmMAcNhld39jv/bgDPqSm9UVDIFEygcRCHqZvtjP+NixY7u723d3BLUCFzN/DiVdyfPYQo/X+QCmBnMPAoEvPP3B7RlkDsGexmYAwM8TogEoZmq1AVyyzNAr0djg0owD1FJbCusAAI1bShF98dDAF5ZIyGKBIIiIcpiCfkmVzY5h7Qp+uJgWa8GexmaYJUKkM9fYJJF6Ky0zJN3M1G2acSAY45Wo8OhfDFO1LhGEZizo93TVz1jr3EIYC/oNzIw/7j+GmYnRqLTZcXq077iNywfE4NXi42jx8BEbNGnGuinGMBhcmnGA6U3OzyKFqfzKk8gzFghEP2OBQAdamfH2sRN4+9gJGAidmqnHW8IRQoRZSTHuYyYi2NUGcLGeAVzOaOoA1XR3JLTahhUSwOrSmwUCgR8IYSzoNzR5aJMOBuJNvgOgTBJh29nZ7nKYgA5FP3RSjV3dnxwBRlO7tFq1DwVEIrVJIAgmQhgL+g3tGz0kOHOKfZHY7pyRCAzAwdwhHao7FM04oCGd4s4zVmmmhloNXRT9EAiCimihKOg3tLbT7BI6MVP7wugUwGrSm1jWrza1K884UM2YNWrGIAifsUAQRIQwFvQb2gvShJAAhLFTiKnxGwdFMw4wmhquy9WaqSWCKAEg6Gv0lX7GGzduDL/11ltTu7+yc4SZWtBv0KIZm5wyTI3fmGUN3ZLa0eYzDkwYa9eMqWNleYFAAAA499xzreeee65VyxxCGAv6De01486iqX3hMlOrEsYM3fKMyZXapLI2tWqfsTBTC9pRsuBPqS379+vaz9ickWEd+LcnTrp+xmvXrrU8++yzSd9///0BtZ9dmKkF/Yb2taWNAQgmk/NaNTFMrGMFLpdQDziymbUHcIloakFv01f7GeuB0IwF/QZPrXaUs0+xvxg0a8YBD/OJy0wdaM5vm5la3bokalML2tGdBhsM+nI/Y60IYSzoN7gE6X1DknD5gJhurvbGROoDuPSsTe3y+QZsptaa2kQQPmNBn6Cv9jPWijBTC/oNrgCui+KjMDoyLKCxJi2asaxjbWq3MA4wgMuhQwCXMFMLepm+2s9YD4QwFvQbXAFcRhVaqqY8Y1ZvHm4PSerM1NBYDpNE1yZBH6Cv9jPWA2GmFvQbHG5hHPhYbUU/dCyH6dKMVdamFkU/BCc7fbGf8YwZM+pnzJhR393eu0JoxoJ+Q1t2jwrNWEvRD1m/RhGkMpqabYoqrTrfmUQAl0AQTIRmLOg3aClC5WoZ4VCxLjPr1kJRUp1n7Pz0KoUxSSS6Ngn6PaKfsUCgA7JTs5Og3mfsUFubWucAroA1Y6cgVW0ulyA0Y0G/J5j9jIWZWtBv0KIZu8aoC+DS30wdqM9Ye2qTMFMLBMFECGNBv8GlGRs0aMZqYpiUAK7Ax/lCUpln7Nak1X7jSUUEt0Ag8BshjAX9Bi3lmQ0aoqll1s9MTQZ1Zmolv0q9mZokkWcsEAQTIYwF/QbZWUJKTTS1QYPPGMyq1vSFWs0YMrRVARNdmwSCoCKEsaDfwBosta5aGWoaRchBMFMHHMDFrO3bLgK4BH2UwsJC0yWXXDKss/MnTpwwPPnkk4lq5t67d29IRkbGaPW78x8hjAX9BodbMw58rDuaWo16qKeZWkMAl5Y9EImuTYK+ydChQ21ffvnloc7OV1ZWGt54440BPbknNYjUJkG/wV0RUkUAl5ZoalnHClyAoh0HLBi1NqsQXZsE7fh22Z7UqmMNuvYzjhsUaf3ZzaM67QY1d+7cQUOGDGl9+OGHjwPAvHnzBlosFse7776bsH///rycnJzQ2267Ld1ms5Esy1i1atXBP/7xj4OOHj1qzsrKyj7vvPPqnnrqqZJLLrlkRG1trcFut9PChQtLbrzxxprO1nQ4HLjuuuuG5OTkRCYlJbV+9dVXByIjI3X/MgjNWNBvcPl71cgkzdHUOn7TyEDqoqm1aOeia5OgD3DjjTdWrVq1Ks71fvXq1bFnnnmmu7b0iy++mHj33XeXFxQU5O/atWtPenp667PPPlucmpraUlBQkL906dLi8PBwed26dQfy8/P3bNiwYd+CBQsGd9V4paioKPTee++tOHDgQF50dLRj2bJlscH4bEIzFvQbXF83g5oALmhtFKGfZkxS4MIYWptVOAO4lGpiQe0kJzhJ6EqDDRZTp05tqqysNBYWFppKS0uN0dHRjmHDhrW6zp911lmNzzzzTEpxcXHIddddVz127NiW9nPIskz333//4C1btkRKkoSKioqQ4uJiY1pamt3XmoMGDWo5++yzmwBgwoQJ1sLCQnMwPpvQjAX9Bn0CuHq3UQTgNFOrKfqhYQ/uhwmRayzoZWbOnFm9fPny2BUrVsTNnj27yvPcnDlzqlavXn0gLCxMvvTSSzPXrFljaT9+6dKlcZWVlcbdu3fvKSgoyI+Pj7e5+hn7IiQkxP1lMxgMbLfbg/I0KjRjQb/BJUi1pTYFvq6eLRQBdT5jzc0qXPtnBlT43AUCvbjpppuqfvOb3wytrq42btiwYW9zc7P7DzI/Pz9k1KhRLaNHj644dOiQeceOHWFTpkyxNjY2ur+BtbW1hoSEBJvZbObPPvvMUlJSEtI7n8QboRkL+g1uM7WKsdpqU+urGavxGUOzz9g5VgRxCXqZyZMnNzc2NkpJSUmtQ4YMsXmee+edd+IyMzNHZ2VlZe/fvz/0rrvuqkxOTnZMmjSpISMjY/Rdd901+I477qjauXNnxJgxY0YtX748Lj09vbm3PosnfmnGRJQMIM3zembeHKxNCQTBQHYHcPVcNLXmPsKd7EWNz1iLMHb9zlgWerGg99m3b5+7WcPIkSNb9+/fnwcAixYtKlu0aFGHXsefffbZYc/3O3bsKPBnHc+5AeCxxx4rV7/rrulWGBPR3wDcCKAAbR3kGMBlwdqUQBAM3I0iVIx1R1MHOI5dpnGdo6lVmam1SFEvM7VAINAbfzTj2QAymblPqPICgVrcjSI0+IztKlsX6laCC4rPOHDNWCcztSj8ITgFKSsrM5x//vkj2x9fv3793uTkZDVtzAPGH2F8GCp9y0R0CYB/QHHTvc7MT/q45hoAj0LRtncy8/Vq1hIIukMidK6sAAAgAElEQVRbowjnHAGOa9OMdRTGBqnno6mFy1hwCpOcnOwoKCgISp9if/FHGNcDyCWibwC4c7aYeV5Xg4jIAOAlABcBKAbwPyJaw8z5HtdkAPgjgKnMXE1Efb5kmeDkxd0oQsVYV55xoAFcbg1WR0crqfAZs6wxolsSmrFAEEz8EcZfOn8CZQqAA8x8CACI6H0AVwDwfPr4DYCXmLkaAJi5QsU6AoFfyEoXQVUBXGq7Nrlzm3XVjCnw3sJao6ldY4UsFgiCQrfCmJnfICIjgBHOQweY2WelknYMAuBZoaUYwBntrskEACL6AYop+1Fm7iD4iehOAHcCQFpamh9LCwQdkaFeHqnt2uSOptYztalXfMauaYQ0FgiCQbeGKyKaBuAAgDcAvAlgHxFN9WNuX9/89t9kI4AMAOcD+BWA14kopsMg5leZeTIzT05MVNUJSyCAzKyqSQSgFAohqNGMXalNqpb1vReJAu7apDXXmUQAl0AQVPy5RTwP4DJmnsrMZwP4BZSgrO4oBpDq8X4wgBIf16xmZhszHwawF4pwFgh0x6FROTSQCmHsNCfrrRkHrKFqzDMWPmNBX2HChAlZwZp7xYoV0QsWLEgGgC+++CIyOzt7lNFonPTWW28FpTmEJ/74jEM8g66YeQ8R+VM+7H8AMogoHcAxANcBaB8p/SkUjfhtIkqAYrbutC+lQKAFGQzSEEllIFIdTa1r0Q+JIDsCdBrLrK3enmv7QhYLepnc3NwOBTvsdjuMRu3VnW+44YZaALUAMGzYsNa33nqr8Mknn0zSPLEf+LP77US0FMA7zvc3AMjtbhAz24noHgBfQfEHv8nMeUT0GIAcZl7jPHcxEeVDyRp5iJkr1XwQQd+i5dBhmIel9/Y2vGBu8/2qQQKp0IzVt23sDJICD+DSWpvaNVb4jAUuvnrl76knjh7RtZ9xQuoQ68/n3t9lN6jw8PAJVqs1d+3atZbHH388ZcCAAbb8/PzwgwcP5l144YXDS0tLQ1paWqQ5c+aUP/jggycAYOXKlVELFy4c5HA4KC4uzv7jjz/u8zX3Cy+8EJ+TkxOxbNmyopEjR7YCgKRnxZ4u8EcYzwFwL4A/QHk+3gjgRX8mZ+bPAXze7thCj9cMYJ7zR3CKUPvZWpQ89BBSX12KyHPP7e3tuHGANVlqjRS4ldYlu3RvoWhXoRkbNdxUhM9Y0AfZtWtXRG5ubl5WVlYrAKxYsaIwKSnJ0dDQQBMmTMi+8cYbq2VZpnvuuWfo+vXrC7KyslrLy8vVlKcPOv5EUzcDeMr5IxB4wXY7IEkgj6dHa04OAODonXdhVMEe9/Efjv2AAzUHcMvoW3p8n4AiR9QGcAGKmVp1bWpdK3CpqASmtaeyuxym+ikEpxbdabA9wbhx4xpdghgAFi9enLRu3boYACgrKzPl5eWFlpeXG6dMmVLvui4pKalHKmoFSqfCmIjeY+ZfEVEufHwFmXliUHcm6HM4GhpQ/rdFCJ80CdFXzgJJEgrGjEXUjBkY9MzT7utqPvigw9iXdryEf+78JwDg8uGXIzY06PEQHZChzVwsqQjgcqUg6VqbWpLAsj/ZhZ4b0Vabuq1RhJDGgr5DeHi420S0du1ay4YNGyw5OTkFFotFnjJlysimpiaJWd+uacGiq1vEQ85/rwJwtY8fQT/AmpMDW5nSBOXEkpdQ+/HHKP3Tn1D1r2Vu/2Hd2rUAANlqRenCv3iN//zQ5xj7r7FuQQwAOeU5PbR7b2RmSBo140BrbbgfY3XWjHuvhaL6KQSCYFJTU2OIjo52WCwWOTc3N3Tnzp0RADB9+vTGrVu3WgoKCkIAoK+aqTsVxsxc7HxZAuAQMx90vh8J4EiwNyboXaz/+x/2ZI3CkRtvwqFfzEDL4cOofu89WC6+GKHjxqH2szXg1lavMbWrV6Pmww+9js3fNN/9+qOZHwEADtX0TsC8rDGAy6jCTN2mGeuc2qQigEtbapPzX6EZB4TVahVBbz3E7Nmza+12O2VmZmYvWLBg4Pjx4xsBYODAgfYXXnih8MorrxwxcuTI7CuvvHKYP/Nt2LAhPCkpadznn38e+/vf/37IiBEjRgdz//4EcG0CcC4RRQPYACWS+joANwdzY4LepWzRIvdrubER1e+9B25pQcLcOaj7/AtUvvYa6r/6yntQF6kFa2atQXp0OuJC47Cveh9KGkowMHJgsLbvEyWAS71AkqCl6Ie+qU2B5xlri6aGiKYOmC+++AJbt27F5MmTMWPGjN7ezimD1WrNBYAZM2bUz5gxo951PCwsjDdu3Ljf15hrrrmm7pprrum2EcS9995bCaASAM477zxreXn5Lp223S3+eLIkZrZCaaW4hJlnAhgX3G0JehNHfT1a8vcg8b57kfbWmwAA21HFUGIaPBgx11wNCg1FyR/atF5ubYXc2OhzvsenPo70aCXNKSokCv8+8m/8fNXPUdNc06NastZUWwNR4NHU7qIfGhZuh6IZB2qmhqZmFSTM1AHBzNi9ezcAICcnB3Z7gD5+Qb/DL2FMRKdDKdix1nmsT9rcBfpQ/+9/AwBCs7NBJhMAwHbsGCg8HFJkJEJSUxE+cYLXmGMPPIja1Wt8zpcY1lbCNMoc5X5965e34orVV8DmsOn9EXwia9SM1VXg0l8zJjXlMHWqTS3M1P5RWVkJq9WKpCSlXsTRo70eeCzw4B//+Ed8VlZWtufPTTfd1KuND/wRxvMA/BXAOmb+iYiGQTFdC05B7MePo/RPjwAAwiZOBDlNz47qahiio9s0JINyPHr2LwEA9V9/jZY9e2BKS0Pqq0u95jxr4Fnu11EhbcL4YK0ShrD7xO7gfJh2aNWMFZ9xYGOCYaYmSYW5WGPRD7eZWghjvzhyRAmrmTlzptd7Qd/gvvvuqywoKMj3/HnnnXeKenNP3d6bmPk7Zr6MmZ8g5U5czsx398DeBEGgvj4PR4uXweFoAgCUl69FS0u5+3zL/jaXi8FicQtdR0MDpNBQ97nE396N+LvuQuI993jNb6+oAJnbrlt60VJIHl0Sos3RHfZ0pK5nblQys7sVohokIndPZH8JhplaUtO1SWZNm2gzUwth7A9HjhxBREQEBg0ahOTkZCGMBd3SbQAXES0DcA8AO4AcAAlE9CQzPxfszQn04/jxb1B8bDmqqn4AIMPhaILNVo2iotcQFTUBp09eCVtFBZp2/wQASF+zGgBAJuVPhJuaQGFtQjbstNMQdtppYNk7rJebm73emyST1/tYc8f84rLGMs2fzx+0tFAEFN9MXwjgIonAAXdt0qnoR8C5Xf2TwsJCDB06FESEIUOGYNu2bbrVTxacmvhjtRvLzHUAZgH4N5TuS7cGc1MC/ZBlO44f/zd+yrsPVVWbEB93DgCgpua/KCp6DQBQV7cDzIwD556H488/DwAIGTQIAECGtvAAKTSsw/wkSTANbIuKHvqhd8EPE3nffOJC4zrMcbS+Z/xpDmYtMUxKowjV/Yw1LNwOkgIPJFP6GWtYVERT+01DQwPq6uowyPkdGjJkCOx2O0pLS3t5Z4K+jD9fzxAiMgK4AsCnzNwK8Xx80nCkaCl27Z4LWW7GxInv47TT3oLBEI76+rYofyIJ9pYGr3FSRITywuNJ3tNM7cmQd5a5X4eOGeN1LkLyriMfE9qhXTWqmqv8+zAaYUCTmVpVOcxgpTb1lpla+Iy7paBAaSqUlqbEA6WmKp1kRRCXoCv8EcavAygCEAtgAxGlAWjoeoigL1BaugqFhS8BAEYMn4/YmNOdZyS0tlYAANLT7wezA4cX3O4eN/SD992vyUMYU1hHzRgATE4NAIBXjWoASA4f4PX+WP2xDuMbbb5TovTGwawxtUlFowjnY6uWKO72qElt0m6mdglj9VP0F0pKShAWFubWjC0WCxISErB161Y4HH2yLPJJTWFhoemSSy7ptJDHiRMnDE8++WRiZ+e7Yu/evSEZGRlBLfbhwp9GEc8DeN71noiOArggmJsSaMPhaEZNzVbk7/kDAOCMM75EZESG+zx5BFTFRE8CAFj37YLZKapCx45tu9ZDuHamGQOAZLFArq/veLydzLhtzG1446c33O+HRQ9Dg61nnu2UClwa+xkHWoHLrRmrXrYDEqkM4NIjtUmYqbultLQUKSkpXvWQzzzzTKxduxbffvstLr744l7cnX5UrdyXaitr1LWFoik5whp3VWZAJoShQ4favvzyy04LFlRWVhreeOONAQ8//PBx7TsMHp3eIojoV85/7/X8AfA7KG0VBX0QZge2516PHTvbNF1PQQwAdnsdAGBk5mOQpBDloNM1HDVjRgft1gWFmjtdd/hXX2L4V192HNNOk4o2R+M/1/0HI2NHAgBSLamw2qxdfia90Fj3QqnAFWjVC3c0tc6pTWp6OWrpoChSmzqlsbERX375Jerr6+FwOFBRUYHk5GSvayZNmoT09HTk5OTAZrOhtrYWb775JsrKeiZ48VRh7ty5gzy13Hnz5g38y1/+kuTSXnNyckLHjh07KisrKzszMzN79+7d5gceeGDw0aNHzVlZWdl33XXX4NraWumss87KzM7OHpWZmZm9fPnyjr4zH+Tn54eMGjUqe8OGDeG+1tH62brSjF1hr6rUe0HvUFr6Cerqdrrfn33W+k6vNZli3MKYnX8J0ZfP9L7IQ4iQ0Tsy2hNjXBwQ1zE4y1fLwmhzNJb8bAm2lG5BfmU+cityO51XT2TWVvTDSIRW1Zqxjj5jg7ra1JoeCCSR2tQZOTk52LJlC5gZEyZMgMPhQEpKitc1RISpU6di+fLlOHjwII4cOYKioiJs3rwZv/zlL3tp59oIVIPVgxtvvLHq/vvvT3NpuatXr45dsmTJkXfffTcBAF588cXEu+++u3zu3LlVzc3NZLfb8eyzzxbPmDEjrKCgIB8AbDYb1q1bdyAuLk4uLS01nnHGGVnXX399jdRFa7WdO3ear7vuuuFvvPHG4bPPPrvplltuSW2/jlY6FcbM/LLz3z9rXkXQIzDLOFL0T1gsozFp4vtwOKwICUno9HqDIRxEioB1CWNu79PyFMYGFaqV3bePLDkiGbNGzMKhmkNodbT6vEZvtFpqDURwyIFJwWD0MyY1ZmrdKnCpn+JUpahIqRVRVlbm1nTba8YAkJ6eDrPZjP379+Pw4cMAgIqKip7b6CnA1KlTmyorK42FhYWm0tJSY3R0tGPYsGHuG8hZZ53V+Mwzz6QUFxeHXHfdddVjx45taT+HLMt0//33D96yZUukJEmoqKgIKS4uNqalpfmUqFVVVcZZs2aN+Oijjw5Onjy52d91AqXbuysRpRHRU0T0IRF97PrRurBAf6qq/gOr9TDSUu+AwRDepSAGAMkQCsmZB8zhys09ZOjQzgcYAs+RJGtTl+dNBhNa5dYeSZmRwT41dX9R+hkHNsb1sfT0GQcawMXM2mtTCzN1pxw/rrgiS0tLceLECUiShDhfViKDASkpKSgoKEBVVRXMZjPKy8vR0qL5Pt6vmDlzZvXy5ctjV6xYETd79myvVIw5c+ZUrV69+kBYWJh86aWXZq5Zs8bSfvzSpUvjKisrjbt3795TUFCQHx8fb2tqaur0G2qxWBwpKSmt69evjwxknUDx5+66BsAyAF9DPBf3Sapr/ofDh/6O6potMJniMGDAJX6N89SMa25yIOuud2BOT/e+yFMz7sKM0xmdNY9wESKFQGYZdrbDRJ2bwfVAD81YVlubWs9+xoY2weiX+dv9QKBHP2MhjD2RZRn19fWwWCyor6/H/v37ERMTA4PBd/n+5ORkFBYWAgCmTZuGb775BiUlJUhv/70TdMpNN91U9Zvf/GZodXW1ccOGDXubm5vdf9hOv27L6NGjKw4dOmTesWNH2JQpU6yNjY3um1dtba0hISHBZjab+bPPPrOUlJSEdLWeyWTiL7/88uD06dMzIiMj5Tlz5lT5Wufyyy/vGMEaAP4I41ZRbavvUlr6CfL3POh+n5x0eVtQVjeYjNFe17pziz3xFCKd3GC6glu7bgIRYlDWtzlsHap16Y3WAC4joU/kGbsEu8x+avoubVakNnVJReEhfP3aEsz8/cOIShjQ/QAowVvMjIyMDGzfvh3l5eXIyMjo9HpP8/XYsWPxzTffoKKiQgjjAJg8eXJzY2OjlJSU1DpkyBDb3r173Texd955J+6jjz6KNxqNnJiYaFu0aFFJUlKSY9KkSQ0ZGRmjL7jggtpHH3207NJLLx0xZsyYUaNHj7amp6c3d7UeAERFRclfffXVgfPPPz8zMjJSzsvLC22/jtbP5Y8wfpGIHgHwFQC3PYWZe6zPo8A3VVU/IH/Pg4iOmoCQkAQcP/E1IiNH+j3eZIp1m6kBQAr3nUfsglQI4+567rmEcaujFeEmXbMkOqA1gMtAhECzRIPTQtE5t4P965/G2oWxe81TVDN22O1490/z4LDbsXH5W5hx//zuB0ERxoDiD961axfsdjvi4+M7vd4ljKOjoxEVFYXQ0FC3mVvgP/v27XNXLRo5cmTr/v378wBg0aJFZYsWLeoQov7ZZ58d9ny/Y8eOAn/W8Zw7ISHB8dNPP+1xnfO1jhb8EcaZAO4AcCnanosZwLl6bkQQGMwO7Nv/OAAga9QiHCn8JwC4zc7+YDJ5N22QfBb18NSMVTg+u7l5u4WxHPwgLhn+/cF3hgQNZmqda1MD8DuIi/VIrzpFzNTMjH1bfkDikHTEDWwrVlNZXASH3Q5jiBn7tvwAW0szTObO8+pdNDQoOfLR0dFITExEaWkpBgzoXKseMGAALrjgAowbNw5EhMTERCGMBQD8uzddA2AoM4sogz4Cs4wdO25DY+N+JCfNQmREBoYPfxAMh9/+Yk8im4eiIbTQd4Utj/s3qQjg6lYYO83kLY7g/3nJzDCo8Hu7MKgK4FIGSDqXw1Tm9nsTzoFaFj01zNT/Xb0S/3nvXwCA037+C+z4ah1+89JbqCxWIqLPue5mrF/2Gt790wO46akXIEldmx5cwjgyMhIzZ87Evn37MG7cuE6vlyQJ557bpsckJCRg3759Wj+WQCNlZWWG888/v4NZcf369XuTk5N7pGyaP3fXXQAs8DBRC3qXsvI1qKr+AcOHPYS0tDsAAKGhKRgz+u+q5otoGogGU2GXFbYABFUztjm69i3rgdauTUY1tald2qvOjSK85u4OHXzGbsX4JI+mdgliANjx1ToAwLZ1n8JkDgVJEsZfdCnWL3sNJ44ewZZV7+Psq2/ocj6XMI6IiEBcXBwGejRN8Yf4+Hg0NjaiubkZod19/wRBIzk52eHKQ+4t/Lm7xgMoIKJ1IrWp9ykuXoH8/AdgsYzBkCF3QpIC11YHD74Fw9J/33bAAcX3aOo4F3nlGQe+Vnc+Rpdm7GmmbnG0BKWTk8waK3BR4C0UXWnJwdCMZT/VdF1M5adA0Y8Wq+LfHTJugvuYwWjE4R3bUFl8BLHJA2EMCcENTyjxqqX79+J4USFK9u3xOR8AWK1WGAwGmM3qCjC5TNqio5PAn7vrE0HfhaBLmBm7ds/BiRPfuI9lZDziVWM6EEZmLvR6T66buuTDBqm16Ec3N2+TQfFxexb+WJH/Lp7f/hw+veJTDI8ZHvianSBDewBXoFZa98NIMDRjfwWja9N6+IxPYjN1TZki8MZfeClGn/czVJUUw2QOxaZ334a1phppY8YDAJJHZCJr6nkozt+NZQ/dAwB44IO1PudsaWnRpNGmpaXBYDBg586dqK2txaZNm3DzzTcjOjq6+8GCUwp/hPFmAM3MzEQ0HMBIKH2NBT1Eff1uL0GcPvR3Hh2YdMBZd0aW7ZCkdk/4XqlN+vuMzQZlPU+f8cr89QCAzUd36iuMNbb0NSDwRhEurTQoPmN/A7h08BmfCtHU1WVK9klMykAkpg0FABzbq2i9LdZGJA5tSy9KychCwQ8b3O+ZfZcTbWlpUa0VA0BoaChGjx6NHTt2YMeOHQCAgwcPYuLEiarnFJyc+PP13AQgjIhSAGwAMBfAm0HdlcCL8vK1IDLBaIxBqHkghg27X98F7MoNlrnr+qrB0Ix9+YxrGpSgmdzSA4Gv1wVKapP68QZVecbKv/pGUyv/+mumhh6FR07SfsbMjGevnYFnr53h1oxjBrTl+iYPH+F+PXR8mwBMGZHpNY/D5jumQQ9f7wUXXICoqCh34Jcokdk1EyZMyArW3CtWrIhesGBBMgA8+uijScOHDx+dmZmZfdZZZ2Xu27fPvwIOKvFH1ZGY2UpEtwNYwsxPEtGOYG5K0IbD0YSS0pVISLgAWSP/r6PmqgduYezjhqOm6IenwFLhMzYYlH18XfoOgD/4t6YfyNDeQjHwfsZBqMAVsJlaz6IfJ5cwri5tq8VQsjcfETGxMHkIT4PRhNN+PgMOWyuSh7cJ4MQh3kU4bK0tMIZ0vBdr1YwBICYmBvPmzQMAHD161B0UJvBNbm5uhxxhu90Oo1FL4qLCDTfcUAugFgAmTZpkfeCBB/ZYLBZ58eLFib///e8Hr1u3rtNWjVrxSxgT0ekArgdwp/OYmuoPAhWUl38Gu70WqYNvRUhIx3q3ekA25QYr+9KMvXzGKipwdSMwjM4ANE+fsQNt9azXHFyDy4dfHvC6vtAawGVQEcDVVoFLw8LtaIum9nMPegRwnaTxW5XHityvD+/YhoGZozpc87PbO3aENYaEIG3MOBT9pNQ2sre2QEkq8aa5ubnLIh+BEhkZ6S4k0tf59NNPUysqKnSt1DNgwADrrFmzuozeDA8Pn2C1WnPXrl1refzxx1MGDBhgy8/PDz948GDehRdeOLy0tDSkpaVFmjNnTvmDDz54AgBWrlwZtXDhwkEOh4Pi4uLsP/74o898shdeeCE+JycnYtmyZUUzZ850l7c855xzGj744AP9/qN94I8wngfgrwDWMfNPRDQMiulaEGRk2YbDhUsQGZmNGD19xO2xKXd1lrtpA+avMPbUAru5ebsCuB798VEMiRoC4lDU0173+eX5y3UTxqwxgEvSVIGr94p+uP8PtPiMiRSBfLJpxiXHvN5HD0jye+ysPyzE9s/X4D/vL4O9k2YOWgO42hMeHo6amhrd5jvV2bVrV0Rubm5eVlZWKwCsWLGiMCkpydHQ0EATJkzIvvHGG6tlWaZ77rln6Pr16wuysrJay8vLA9Yqli5dmnjhhRfW6v8J2uhWGDPzdwC+83h/CMDdwdyUQKG6+kc0Nx/DuLF/1vVm3oGuzNQeUDcFEHzSjfrmag5R31qP/9vyf9hx3NsDEtNN56lA0BrAZSQVAVxBrE0dcJ6x1r8hifq8arxn0/fY+N6/cPNTLyIs0oKashKER8eguaEBssOOuIGD/Z7LZA5FTLKSN9yZz1gPM7Un4eHhKCnRXOa4R+hOg+0Jxo0b1+gSxACwePHipHXr1sUAQFlZmSkvLy+0vLzcOGXKlHrXdUlJSQE9U7/88stxO3fuDF+6dOne7q9Wjz8tFEcQ0ctE9DkR/dv148/kRHQJEe0logNE9HAX111FRExEkwPZ/KlOfYMS6RkTc0ZwF3JqxrLc8YbjlWdsVCOM/UttAoBmR1u99qZj18LRnIQ95ScCX7MTHFpTm6DCTB2E2tSurk3+l8PU6YGAyG/TeG/x46r30FB5Aru+/gKA4jOOTRmI4ZOnAACGn35mQPNJzr95R/s+31A6NumtGYeFhaGpqeu2o4I2wsPD3X+Ra9eutWzYsMGSk5NTsHfv3vxRo0Y1NTU1SZ1FwvvDp59+annmmWdSPv/88wNhYWFBfRL1R1FYCWAPgP8D8GePny4hIgOAl6DUtM4G8CsiyvZxnQXAvQC2+r/t/kFjwz6YzckwmaKCu5DLTN1NNDX81YwDCOAyehQt8Yyo/sesK+GwjkCtfAA2Hw8JatCc2kSkuhxmMGpT+60Zu83U2vZAEvq8mVpypt/95/1leP76K1C85yfEpgzCpXfPw7WPPulOafIXg3M+2dHxu9HaqihkemvGdrvdPbfAf2pqagzR0dEOi8Ui5+bmhu7cuTMCAKZPn964detWS0FBQQgA+Gum/uGHH8J+97vfDVm9evWBQYMGdXNz1I4/PmOZmV9UMfcUAAecZm0Q0fsArgDQvuTY4wCeAvAgBF7Y7LUICQlqzIBzoS6EsUbNuNsALmr7E6xvVeIlHE2DcdmobLy7ezR2Nf+A4rpSpMekBbx2e5SiH+rHS6TItUC6P7m0Vy0aeXvcWUY9WA7TvXAfN1M31lS7X8tObTY1eyxMoaEYPGpMwPNJzjgJ2d5RM25x+pH1FMZhzvrwTU1NCPERvS3onNmzZ9e++uqriZmZmdnDhw9vHj9+fCMADBw40P7CCy8UXnnllSNkWUZ8fLxt8+bN+7ub76GHHkq1Wq2Gq6++erhzntbvvvtO33xLD/wRxquJ6E4An8C7hWJdN+MGAfD0KRQD8LK3EtEEAKnMvJaIOhXGzvXvBJSKNf0Flm2QqAe+kG5h7MOV4ilE/NWMVQRwAUBFk5JfOSg6EkSE0wcNx66DwP+KD+oijJmVzktqMTo/l4MDkGtBiKZ2makDjqbW+jwgUZ+uTW1vbUVzQz3OvvoGhISFIyUjE/mb1iPzzHNUzyl1oRk3NytuFb0DuAClzKaowuUbq9WaCwAzZsyonzFjhjviOSwsjDdu3OhTyF5zzTV111xzTbe1p++9995KAJUAsHnz5h7t4OGPML7D+a+naZoBdHd39PXVd3+TSanl+DyAW7vbADO/CuBVAJg8eXLfvRvojMytIKkHhHGrcqPxKYw98FszVmmmdtHoqAQAZCcOAQ4C+Sf2AZju39pdoLVRhMEtjBkmP4W6HMRoav/N1PpoxiSh24er3qSxpgoAYIlPwJjpFwGAz1SmQOjKZ+zy7ertMwYUYSzoX/gTTZ2qcu5iAJ5jBwPwDBO0ABgDYL3zRpUMYA0RXc7MOSrXPKWQ5aOOEnsAACAASURBVFYYjR1zG3XHaYJjX3daL81YjXrnv5naxego5UY6LnkY2GHGoRp98uwdzJo0Y9endwQgkYLhM5Z6K5q6j5upG6oUYRwZq18+fnufcW1tLaqqqpCenu7WjMN89gFXh0szFkFcweUf//hH/CuvvOKV53b66ac3vPPOO0WdjQk2fpUsIaIsKEFY7kdAZn63m2H/A5BBROkAjgG4DkrhENf4WgDuvBUiWg/gQSGI22DZBqkHNGO2OU1w3ZipSU2FG7lrW6rBh+l7avy1AIABllCwPQbHm/SJqJYDMS/7wKUZB2KpZZl1jaQGAs8zdkd0a30gkKhPN4poqFaEcYSOwri9z/jDDz/EsWPHsGDBAqEZn8Tcd999lffdd19lb+/DE39Smx6BYiL+J5TI6L8DuKq7caxEA90D4Cso0dgfMnMeET1GRPpUcThVYAZ2vAu0eJfBk7kVRKZOBumIzakZd+eEVKEZq2kskBQV5lyOYKZo1NmqAp7DF1q7Nrl8xoHUp2ZZX60Y8PQZBxrApW1dor7pM7Y1N6O6rMRtptZVGBu9NeNjx5QiIidOnAiKZuwZwCXoX/jz9bwWisOulJlvAjAefmrUzPw5M2cy83BmfsJ5bCEzr/Fx7fn9UitmBvI+Bj6dC/y4xOuULLf0kGbsMlP7EsYq+hl7BXAFWiRDQmhIm7YcZohCs6O+ixGBzK1NHrlkaiDpTSyz7sLY5X/2WzPWyWcMCX3STP3d20vx5n13orq0BJLBiDCLfqmALs3Y5TOWnA+kx48fdwtMPaOpXfWVv//+e93mFJwc+HNvamIlssfuzAkuAzAsuNvqJ9hbgNcuAFberrwPifA6LfeQmRpuM3VHYewpV/3u2uQVwBXYVppLroHRQ2hEmCJg4+YuRviPXgFccgACSeYgCGNXO8MANWPNQWQS9ck845++/xoAcGxvPiJiY3UNlnP7jO3Kd8SludbW1ro7NkmqYik6JzY2FoDQjvsb/vwV5RJRDJS2iTkA/gtge1B31R8o2gp8dBtQ4vGrNHs/0ctyK6QeMFO7fMbdpjap6WccoDRmeyQMHsIrMiQMMlp16aMrawzgMrg040ACuGTWtZcxEHhqk9vgoTWamqgvKsYwmRWf7fHCQ7oGbwFt0dSunGXZGQNRX1+PpqYmXf3FLs4991wAQhj3N7q8u5LyiPkoM9cAeImIvgIQxcxCGGth+zvAmnuU15NvB5rrgJ9WAu0qTXFPpzZ1E50TjH7GHS53hMPooWlEmcMAqw21TTbEhGv7XTigLaDY4PYZ+z+GHaxrjjEQuJlar9Qm9NEKXObISNhaFOuJJWGArnN75hkzs7vQR11dHWRZ1tVf7MJV7MPWST1sQRvz5s0bGBkZ6XjsscfKe3svWulSGDMzE9FaAJOc74NWfaTf0FIPfPMoEDUYuOZfwKBJQFO1IozbFRboCTO13NICudEZuelL1dKqGQcsjCPhKfOjw8JBtXaU1TZrFsbMDIMWzRhqzNTQXTN25xn7uY+22tRaF+6bZmqXZgwAUQmJus7t9hnbHXA4HF6asSRJQdGMTSbFGiaE8cmJzWZz/x8Ggj931/8S0UShDetAcy2w6g7AegK441tgsLMvhqvwRTvNuCfM1I6aGpBTBvuMpvbqZ6wimrqb1KYO19sjYPDQjGPDlLzLY7X1yErRFpij3Wes/BtwAJfOuU0u4S77uxGd8oz7qpnaszqWJT44wphlGXZ72zp1dXUIDQ1FVJT+deNdPmg5wO9OT5O/Z35qY8M+XfsZR0RmWrNHLe6yG9T8+fOTP/jgg4SBAwe2xsfH2yZMmGDNy8szz5kzJ62qqsoYGhoqv/7660cmTJjQfPToUePtt98+pKioyAwAS5YsOXLRRRc1dtb3ODw8fMItt9xSsXHjxqjo6GjHE088UTx//vzUkpKSkMWLFxfdcMMNPtsovvDCC/FffPFFdEtLi2S1WqUtW7YEXL2rU2FMREZnetI5AH5DRAcBNEIJr2VmnhjoYv2awxuBZbMAMPCL59oEMQC4SkJ6NEqQZTsAOfgBXHa7263bnZna737GngRw8zZSCACDVwBXXJgS1FZcWwulwqp69GgUAQTWuSko0dSBasY6NYroq2Zqz7rRumvGbsHYphWHhoaisbERDocjKGZqgyu3uY8L495g06ZN4Z988knc7t278202G0477bTsCRMmWO+4444hr7766pGxY8e2fPfddxFz585N27Jly745c+akTZs2rX7hwoUH7XY7amtrDYDvvsfJycmOpqYmafr06fWvvPLKsYsuumj4I488MmjTpk37tm/fHnrbbbeldyaMAWD79u2Ru3btygu0RaOLrjTj/wKYCGCWmokFHhRtBZZdoZiBpz8CnP5r7/OSUxjLbU/ert7CwfYZsywDsvMm3V3RD1XC2P+b98Nj/oX5+Ue8ArgSIiIBAGV12tObtOYZG1TlGQfBZ6w6mlrjwn20n7HssMMcEQGH3Y60seN1nZucwphl2Z3eFBMTg7KysqAFcLkeAHyV4OxLdKfBBoPvv/8+8rLLLquxWCwyAFx88cU1zc3NUm5ubqSroQMAtLa2EgBs3rzZsnLlysOAkjYWHx/vAHz3PU5OTm40mUx81VVX1QHA6NGjm8xms2w2m3nKlClNx44d6/JmPG3atDq1ghjoWhgTADDzQbWTCwD89DHw8Z2Kj/jGlUDiyI7XuKpQeWnGSgu1oGvGDkcAZmoVwjgA1dhIyo3NUzMOMymfv6Jee0UiB/sumO4vLjN1IPqKHIxoandtan83oVNt6j5a9MNhtyNr6vm48NdzdZ+bpDYztUtTjY2NRVlZGQD0azN1b9He7SPLMiwWi72goKDbRhCAd99ji8UiT5kyZWRTU5MEAEajkV2/f0mSYDabGVCsFQ6Ho8svkGdvZTV09cyeSETzOvvRsmi/Ye8XwMrbgLh04K4NvgUxoAg8yeTlM5admnGwfcYss99malLR0o19tJ5rz66bd2H7TdthJMXk56kZuxpJHK/XlubhMukatGjG6Ftmav/LYeoVTU19slGE7LDDoKZUqx94mqldmmpycrL7fFJSks9x+qwphHF7LrjggoZ169bFNDQ0UHV1tfT111/HhIeHy4MHD2598803YwHl9/bjjz+GAcDUqVPrn3766UQAsNvtqKqqkjrre9zbdCWMDQAioTR08PUj6AzZAfzn78D7Nyjvb/oECO8m/1G2Af953v3W1VuYSI02GgCyo03V81X0w/O1ikpD3NJ9wQ4igkkywe4MSPJMbTI4P//xRm3CWI9UWzUVuGRZ/2jqNs04QDO15n7G6LM+Y0mV1aZ7XGZq2eGtGbsYMEDfVCpACOOuOOecc6xXXnll1ZgxY0bPmDFj+JQpUxoA4L333jv01ltvJYwcOTI7IyNj9KpVq2IA4JVXXinasGGDJTMzM3vMmDHZ27dvD5s9e3at3W6nzMzM7AULFgx09T3ubbp6nCxl5sd6bCenEtveBr75CxASCdywEogeHPAULCtP4eSjq5GesENu04y7M1Or0IzlJv+rZzmcN3qDoW1NVyOJygaNwliH8sxGNQFcQanApU4zJoNGM7VEuhRf0RtHEDVjIgKR5OUzNhgMuOuuu3DixAl3lyU9MRg6b9soABYvXly2ePHisvbHN23a1KGXcWpqqv3bb7/t4GrtrO+xq1cyADz33HMlnZ1rj2cfZLV06zMWBMj6J4H1i4Ch04BbPvM/aibjYqChwuOASxjrHP3THtkBct9fuzFTq9CM5Sb/fb12uaMp2dVisaqpGQ6ZvUzYAe3D+cShRwBXwGbqIKU2+a0ZO3Q0U9v6ljBmWQbLctA0Y0DRjtkjmlqSJKSkpCAlJSUo6wnNuH/SlTD+WY/t4lRh29uKIAaAy54JLHyVDF7RzO7SlEE2U7NDdsvg7sphkp+J7F59kR3+31AczpuPp8B1acYyO1BtbUVCpLqi/C55pEcAV6B5xpJGjbQ9vaUZO5Matc2hM25t9f/bO/Mwuco6339+59TSazpJZyOks7AGEpQlIhlccMYFryKO4Mg8A/eqg+B278N1vHNRZ5h5HJ/BdWbuqM+96OAdB1SYEUdAuSqjAqIiiwRIIBDCYkJI0ums3enuqnPOe/8451RXuqu761Sdt6qS+n2eJ09313Letzvd9a3vb3sz9morHMchmOSMbaJi3Lrcdtttcz71qU8dEeocGBgYv/vuu+sudJ5WjI0x6Zxb1y7c+wW4529h2bnw3h9CJmFI13GPOPs3DhlLvefezYYpC1NXqs4pF+Oq31yYaT6fmThMnalQwAUBQ8O1i3EcXq3HGcdzrYME31Ng5Tzj8GOjnbE4Uv0IzgYRD/yw6oxdd0qY2iYqxq3LJZdccvCSSy6pqmo7KZZf6duEn18PP/8MnPAGeM9NyYUYwlfYCs7YegGX78/sjGvgCPOUwEl5lXLG0fcvEjA0PF7znuKXtXp+mploW8n7jC0N/ajytTq1amppvWrqeOCHU9MhJtUROuMjw9Q20aEf7Ynd6qB24IX74d7Php+/+5+ho8a+Q8cNK5sj4jYj22JsgmAiZzzbbOpqr1n2ip2k4GdGZyw+e0YKifcSE5Sccc2XKMsZJ1jXRjV10oMiAgOOHJNHKPpe1AKYsZ0zbrwz1gKu9kKdcT1svgu+eVH4+Qd+VrsQw5ScMaZBBVzlzrhSAVedYlyTM3amOuMwTF2/M64rTN1K4zCl+jC18U39rpjIGLeYGAelnLE9XxGLcaOcsYap2xMV41oZGQonay04FT70q/D0pXqY7IwbVcBVNvSjYp9xLeJVW8q4zBmX9RlHBVyOE+aMa6V0VkLNV5gIUyfRo7C1qY5Fp8FJMg0rMPUXb4WLtlwB10TO2GKY2nW1gEuxjopxrdz3BSgMh6HpxWvqv564R4hhwwq4Ar+szzidsFj+rFdy/+mxclX/ghI743ITF7c29XQ4DI3U44zTa21KkjMO/PTHYUJYGV11NbUfpOKMwzB1/ZdJE9+zL5DiOAR+43LG8RvgVuzpbkVeeOGF7IUXXnjCdPfv2bPH/exnP1vzCSJnnXXW6lqfmwQV46QYA//2XvjN/4Z174dFKf0/Oc6knHFjCriMH0zMpk6rOieb4R8vjvedJGcc4E7KbcY54zmdDntScMb1/DRLE7gSfE82wtQQDeBosDMWaT2BKDlji2Fqp0k5Y3XG1bFy5crij370o+emu39oaMi98cYbax6V9uijj26u9blJ0AKupGz9KWz6d+gbgDd/Jr3rTukzbkwBF75XFqZOqZo6uqBJ+OLtVRjqEeeMezqcOnPGKTjjuLUpUZg6hTGUFXAkQTV1Sjnj1izgsi/Gkwu4NGcccs1TvxvYPDKW6giy1d0dh//htOXTngb1oQ996PgVK1YUrr322kGAj33sY0t7e3v9b3/72wu2bNmy6eGHH+543/vet6pYLEoQBNx2221bP/GJTxy/bdu2/OrVq09//etff/Dzn//8jgsvvPCkAwcOuJ7nyXXXXbfj8ssv3z/dml1dXWfNNH0rLdQZJ8Ebh59cF465/OD9kEvx93BKzrgxs6mN55W1NqXzx18SYEmWYwwCc0QlNUzkjHvyDkN1VVOHHxs9DjOwMIELop7faveRljuX1gtTx87Ytd7aNFHAZdsZa5h6ei6//PK9t912W2nQ/+233z7vvPPOK82W/vKXv7zwwx/+8K7Nmzc/+fjjjz+1atWqwpe+9KXtAwMD45s3b37yhhtu2N7V1RX88Ic/fPbJJ5986t57733mk5/85LJWeOOjzjgJP/007N4El9wInXPTvfaUaupg4naLmOKEM57t1KbaFqj+oZWccZwz7u6Q+gq4oo/16GK8taR9xjZyxo4roeOtcg+phKlbcDb1RJ+xzZyxi/H9hoWpw3nY0vLOeCYHa4vzzz9/dGhoKPPCCy9kX3755UxfX59/wgknlF4Y1q9fP/LFL37xuO3bt+cuu+yyfWecccaUcFoQBHLNNdcse+CBB3ocx2H37t257du3Z5YvX+5NfmwjUWdcLS/8En79VTjjj2DNu9K//pQJXI1pbTKeh5johTolZxxE1zEJnbE/gzPuzgvD4x5jxdpC6aU+4zrqqeMCrkRhals5Y0ngjFMLU9NyYeqGVFNPcsa2w9TlaypTueiii/bdfPPN8771rW/Nv+SSS46YFPnBD35w7+233/5sZ2dn8Na3vvWUO+64Y8oJgzfccMP8oaGhzBNPPPHU5s2bn+zv7y/G5xk3E3XG1TB+CL7/IZi3At7+92GxVdpMM5ta6io5mh3jxWcoO+lN4IrtcMJZxqEzPvJnGxdwdeZCMRkaKXD83M7EeyoVcNU19CP8mKSAKwiMnV+XBAVcaTnjpGmHRhDnjK33GZvGFXBB+Gar1aIQrcIVV1yx9wMf+MDKffv2Ze69996nx8bGSr/cTz75ZO60004bX7Nmze7nnnsuv2HDhs5zzz338MjISOmv8MCBA+6CBQuK+Xze3Hnnnb07duyoYWRi+jT93cBRwQ//DA5sg3f+H8j32FljSjV1XMBl+b8oejETnNSqqSeuIySqpvYrOOMoTN+Vj8S4xiKuNAq4MjVM4LLljJ0k1dRpOeNk/50NIR760YgJXOqMW4N169aNjYyMOIsXLy6sWLGiWH7fTTfdNP+UU05Zs3r16tO3bNnScfXVVw8tWbLEP+ecc4ZPPvnkNVdfffWyK6+8cu9jjz3WvXbt2tNuvvnm+atWrZrxnFcbNR+VUGc8Gy8/Do/fCmf/Z1ix3t46007gsu2M4+EiR87Gruua8Tt6RzAJ+4yn5IwjZ9wRO+Ma88ZpFHDFIe5kOWM71dTiJDy1ya1fQBK1UzUI/xgc+gEqxrPxzDPPlA5rOPXUUwtbtmzZBHD99dfvvP7666ecdXznnXc+X/71hg0bqmpX2rlzp9vX19eQXLI649nY/IPw4/qP2l1nSjV1gwq4YmcsTtMLuOI+43JiZ9yZnQhT10IaBVzxBK5E4zCNzT7jKh98LFdTH4NDP+I1NEzdXF544YXseeedd9pHPvKRXY1YT53xTIzuhwe/DkvPhoWn2l1r2pyx7QKuYmmdtFqbai3g8mYo4Mplw+vUHKZOoYAr3lsxgTsMApNKhHgyScLUqU7gajGBCEoHRdg9zzjuM26EKwaOimrqY4mdO3e6F1xwwZQX+QceeOCpJUuWNOTEDhXjmfjNDTC6Fy690f5akehEx/w07ghFr6yfOa0+4xoLuAIzfWuTIwEd2dp7jScOiqjp6QBkaxiHGRZPpf+GKtHZwgFINqXWpiQJ8wZQCh036NSmRomxhqkby5IlS/zNmzdbOae4WjRMPR1DW+Hez8GJvx/+s00supEIN6qAKw5TI05pBGfd1zRlBVxJnLFfYQJX9CYlMAH93Xn21OmM3TqccSzGxYRDP2w440TV1Gk5YzecwNVK4dPGtDa5pdamRoSowzVbNkwdBEHQmIqmY5DoZ1fxXZaK8XT8x1+Bm4X/9MXGrBf/kcd540YVcBXLcsY2nHGi2dRTxdgRB0ccPOOxoCdXewFX6Xo1PR0IQ4eugJeomtrSOMyE1dSpDP2Ir9FCRVyNGfrhlIZ+aJiajYODg30qyMkJgkAGBwf7gI2V7rcaphaRC4H/RTif/5+MMZ+ddP/HgCsBDxgE3m+MedHmnqpiZAie+Qmsex/0n9iYNWPRDUJxnBDGBhRwuS7gpBamLpGwV9KvEKaGsIjLD3z6e/LsPjRjF8KM106DrEiinLHNgyKCBBO4UjnPOBJj4xvbdYVVM9FnbDdn3Axn3Ipi7HnelTt37vynnTt3rkXNXFICYKPneVdWutOaGEto6b4KvAnYDjwkIncYY8rj8o8C64wxh0XkQ8DngffY2lPVPHgD+OPhqUyNwpkUpo5nUzu2xbiIZDKIuKkN/QhKleCS+DzjSn3AGSeDF3j0d+d46uWDNe4p/Jits2cwI5I8Z2xBjN2M4Febv01r6Eec+26hvPHEqU3HXs64FcPU55xzzm7gHc3ex7GIzXc25wLPGmOeM8YUgFuAi8sfYIz5uTHmcPTlA8Ayi/upjsN74cGvwalvs19BXU7UTxuHqSfOM7ZdwOUjrhu1NqU09MPEpzZJovOMKxVwQeSMjc/8KExdy4tULKCZOsU4K5IsZ2xszaZ2CPzqfrbpO+PWcWyloR8Wc8biuJjAb6gzbuEwtWIJm2Hq44HyQeLbgVfP8Pg/Bf6fxf1Ux2PfgdF98Po/b+y6pQKuqC2Ixs2mJpsNW6imccYda9eSW7my+mvWWE3tB6Y0/7mc2Bkv6s5T8AMOjXvM6UgWlozFuNL1k5BxEjpj32Djv9DNCMXxBDnjtFqboKVyxqUwtUXH6rhOaeiHVlMrtrApxpX++iv+FYvI5cA64PXT3H8VcBXA8uXL09rfVPwibPg2LFwNS8+0t04lJhVwNaq1yXgeksmAuNMO/Vj13X+r8erJ5idGXV1TiJ1xf084QnZouJBYjOPIaqZOTcqJUKi2itkYjKElnHEa7VXxNUySCjbLBL6HOA5i0bGK4zZFjFsxTK3Yw6bt2g4MlH29DNgx+UEi8kbgU8A7jDEV+1aMMV8zxqwzxqxbuHChlc0SBPDtP4JdG+G1H7ezxkxMam2KHbJ9MY5zxmLhPONkZ7JOW8DlxGKcB2ob/JFWmDpJzjgusHIy6f+ZuRkHv1pR9E0qf+mlMHWLOWObZxlDPPQjDFNrNbViC5ti/BBwsoisEpEccBlwR/kDROQs4AZCId5tcS+zs/E22PozeNPfwBmXNn790tCPI52x9YLFyBmnOfRjcHQQiA5nSBimrljAJRMFXEBNvcZxNXUlsU9C1qk+Z+x70fjENIqnJuFmhMCr0hn7KQ0eKR1b1ToiEfi+1eItiELGfuiM272aWrGHtd8sE5YDfxT4MfAU8K/GmE0i8mkRiavxvgD0AP8mIhtE5I5pLmefx74D808MZ1A36JSOI5g89KNRzrgYhalJbzb1PdvuAaAYFBNVU09bwOW4eIHHwt7QGQ/W0GucqjOu0hnGzti1MIHLyTj41YqihdamViHwPavFWwDiulpNrVjH6m+xMeYu4K5Jt11X9vkbba5fNbs2wXP3wPqP2DmruBomO+OGFnBlEPFSC1Ofs/gcvv7E1+nMdqVSwFXKGXfnEIHBQ7WEqeNrJX7qESSppi6FqW04Y1cIqgxTmyCo/xuHlmxt8j3P6lnGEJ/aFA79yGbt9TOXo2Hq9kObto2Bu/4cOubAa/578/YxuZq6YTljD8nMXE2dlKC0dye5GFdwcBkngx/4ZFyH/u5cTWLsNyVnHP4cXAs546TOOI1q6vgaLdXa5PnWnbHjuKVTm7SaWrGFHhSx6d/hxfvhbX8HXfObt49pxmHafr9kvCLiulE1dUp9xvF1HMEE1Qv8TGLsRUNQFvTka3TGKfUZO9VP4IoLrJrpjI0x4dyftGZT04Jhats5Y9ctjcPU2dSKLdrbGftFuPuvYMkZcM57m7uXKTljDxEXsZ2/9vyyauqUJ3AtmI/38stVP8+fZkBGPA4TYNGcDgZrGImZWp9xojB1VMBVbz9VBap2xtEbh1RnU7eQGPt+A5yx65TC1FpNrdiivcV44/fgwO/gDX8xkbNtFpNyxoHxELEfuCjljHFJ6+T4khh3dpQOoqjqedPljKMCLoCFNTrjNPuME+eMLbip2BnP5p5KLjaVAq6oz7iFWpsCr2g9ZyyOW6qm1jC1Yov2FePhwfBkpoWr4eQ3N3s30zhj+8Uicc4YC6c2iZPsmtP1GWckgx/9XBb25hkcHk8cwkutmtpJXk1tI0wd9y7PdliEKUb5+zTy1q3a2mRZIB13ImesYWrFFu0rxvd9PpxD/a6vN6+CupySM44LuBrkjItFJJuNiq3SHvohiUYnBgEzHhQBsKg3T9E37D9cTLSntMLUSaqp4z5jGwVccbuUP0uvcSzGTq5+wWrF1qaGDP0oq6ZWZ6zYogVUqEm89AgsfzUc94pm7yRksjMOPOuV1ACmUEByOcBJPWcsjpPooIiwgGvq7XFrE8CiOWGv8e6EoerYzTZlApcVZyxHrDEdphC1yOXScMat19rUmKEfLhijOWPFKu0pxsUx2LkRlrSIEMOUaurAeDhOA8LUhQKSy4ZHKKZUTR3EueekrU0zjcOMC7h6OwASn2s84YwTPW0KiSZwlVqbbEzgqtIZl8Q4BWfcgq1Nvmd/6EccBtehH4pN2lOMd20Mzytevr7ZO5lgSs642JgwdaGAk8shpFdNXTq0KeELSjDDOMySM46mcO0+mMwZx2au3vOM844wnjhnbKHP2K3WGUdRihSccSuGqQPfw7U8iCM+hKLROWN1xu1Fe4rx3ufCjwtOae4+ypk8gStoXM6YbDYKiadcTZ3wPOOZnHExCHPEcZh6V43OuN4wdYfjMFalMwyszqauzhkHxfScMa3Y2tSgCVyGxjpjDVO3H+0pxi8/Bk4W5q1o9k4mmOKM/YY64zSrqUth6hpyxpWccXnOuCuXoSefSeyMYzGudP0kdDjCWJXfU0Oc8SyDP0ph6mwKYepsfIRi64hE0KACrhgNUyu2aE8x/t0DMPBqyOSbvZMJJlVTB6aI4zRGjCWbs1JNLY4kDlPPNA4zZlFv8l5jzxhyKQxQ6XAcRoPZ+3uh7KAImznjWVx6HKZ20ghTR4IeC3wr4HtFHNvO2HHD+gcaK8bqjNuL9hTjwjB0zWv2Lo6kYp+xfTEOisWJamrSeZEtFYKJpFLAlXEmcsYQ9honLeAqBIZsCoMvOiOXO1ZF3jgWSivV1EmdcVoFXBkptUu1Ao0KU8cnuTUqZ6xh6vajPcW4eBiy3c3exZFMqqYOc8Z2C1NMEEAkxiJuamGxidYmN1GYero+Y1cmJnABLJ7Tkbi1qWAM+TTE2InFePbvK/DsU5KiCAAAFaxJREFUhandbJwznvkNVJoFXBC641YS48DzrDtjcR1M9HupzlixRXuKceEwZDubvYsjmWY2tU1MMSyKkmwWQVI7takUpk4Y+g6d8dTby8dhQhim3n0w2RSuQhDUXUkN0BE50mrEuBi50mwaxVOTiK9ZHJ+lgCvFnDGAk3VK12wFfN+374ydCWesOWPFFu0pxsXDkGs1Zzx1NrVjOUxdEuNcLjq1KeVqakcSzTGe7TzjmEVz8owWfYbHq597XTCGXAohxo7YGVdRUexFopVJyZWWk83HYjyLMy76kHFSOUIRwnB3azlj+7Opy8PUWk2t2KL9xNiYKEzd1eydHEklZ2x56IcpFMKlczm84n5GRrawd+8v677unc/dGV43wdCPIIgPVZgqGlknO6mAKx78UX2ouhikV8AFMFrFC6VXCEAmQsppkomcsTeLSzWFIJXirRjJOi0lxr7n4WTs/p04rouJCriylnuaS2tqmLrtaD8xLo6GH3MtJsZTZlPbH/oxIcZZDh56AoDfbftG3dd9ZNcj4ScJWpv8GWZHu+KWzjOGicEfuw5WX8RVCAy5FNxhR3SNanqNvYJPJmfnGMyqnXHBT6fHOCIU49YIUxtjGlPA5bilmo5GirGGqduL9hXjTKvljKP/iiP6jC3njMuccWkbKax56SmXAtCR7ay6d9mfwRmXj8MEWLUwTDFs2TVc9Z4KJt1q6mqdcdZCiBomQt/F2ZxxMSj1B6eB5NxSUVizMUEAxlg/tUma4IzD88Wra6FTjg3aUIwPhx9b1hlHOePAw7FdTR2LcTYLhEIlKfxKuOIyNz830alNQeyMK4nxJGe8ZE4HC3ryPLZtf9V7GvWDUiV0PUxUU1eXM86kVDg1mbiAy2uGM26RAi7fC2seXOthaqfkjDOWXXhpzbIRnEp70IZiHDnjls8ZFxHLQz/KC7hct/PIfdRBYAIciYqGqj1UIZg+TJ1xMgQmmGiZEuHMgT42bK9ejA/7AV0ptBh1xM64ijB1seBbKd6CsOc3k3UozuJSg4KfWlsTtFbO2PfCN2i2w9Sum8U0IUwNqDNuI9pQjCNn3GpiPHk2dQOGfsTO2MnlcN0w9CvUH8oNTIAr0dSiKt/Zx+cTd1ZwcW70BqG8ovrMgbk8NzjCgSrPNT4cpCTGTvWtTV4hKOV2bZDJu1U44yCVs4xjJOe2TGtT4If7sN1n7OaypTRSI8PUoM64nWhDMY6dcavljKMXzLv/CmjMbOrynHEm0wuAH4zWfV3f+LiOi3TkCUZGqnpOLMaL53RMuS8TRQjK88ZnDoQT1B6r0h0fTilM3R3lJ4cTFHDZIpt3KczS3pV2mNrpyGBGq28ps4lfjMPUdv9Osrl805yxinH70IZi3KLOuGt++DEaeh8ERet9xsERYtwDgO9VJ54z4Qc+rrhklyzBP3CgJPoz4UUvOpkKoyNLYlzmjF8x0IcIbKgyb5xWmHpelAPeW5xdkIrj9sLUAPmuDIXDs4lxugVcTmcGUwxa4rCI4njY2pbN2Z0x72Yb74w1TN1+tKEYt6gzdrOw4FRY+RoAgmAccXKzPKk+guFQeJ3ubub0vgKATLav7uv6xg9zxtELl/FnD2vGOePMNAVcwBFTuOZ0ZDlxYU/DxTjvOPS6DkOF2cV47LBHvsvei3e+K8v4bGJcTNkZd0VvFsea7469QijGmbxdMc6UOWMt4FJs0cZi3GLOGEJB9j2MCfD9kVLo2Bberp3hsr29rFhxFQDz5p1X93VLOeMopGu8KlykP0M1tTNVjCHMG2/Ytn9W9xAYw2hKOWOA/lyGoSqc8fhIkY5um2KcYWwGMTaBIRj1cDrTE5D4WkELhKqLY2GfeTY/NbWRJplsFpywX1xzxoot2lCM4zB1izljACcDQRHfHwEM2cwcq8vt+tvrw2V7eyf6i1M4RtE3YZhaopB7NWIcO+NsBcGsVMAFoRjvHSmwbe/Mee6D0WEKfZl0HGJ/dnYxDgLD+GGPfLc9J5XvyjA+QwFbMOqBAacnPQGRjhYS48gZWxfjXB7jumRdOwNcKqFh6vajDcW4RcPUEDnjIsXiQQDrzjjG6e4uibFJ4bCIUW+UrJtFYvGrIkwd54wrOeNsNBa0vIALQjEGeHTbvhmvHQtnfzYdYezPZmbNGce53A7bYeqR6fcRDIe5ejdFMS6FqUeqq2K3SXE8dMa2w9RuLotxXDIWTt+ads0oquRX8bejHBu0oRi3aAEXgJOFwMPzDwGQseyMY8R1iX8VSmcR18GO4R0s61kGmeqd8Uv7wzdJFXPGcZjaHHmd1Ut66cg6s+aN90bjG+enJca5DEOztPcM7w9dW1efvbx/Z28W3wsYn8alekOhWLm96YlVZn7oQr29yc6TtoHXqDB1JotxMxV/N22Ri6bijY8nOypUOXppQzEeBQQydt9N14SbAb/IwQMbgMY5Y5jIUaURph7zxujKdpWFqWd/d783cnELKwhHKUw9yRlnXIczju+bVYzjYqv+XDpivCAKU/szhBAPRWLV229PKPoWhNGdg4OVw/SxGGcWp/fG0+nOInkXf6j5YlyIxbjD7t+yOA4mm0vl1K9qicW4UEUngnJs0IZiHJ3Y1KDcTyKcLARFNj/9KYDSIA5bZJcupe/ii6Ov0g1Td2Y6S2HqanqNhwseOddhSYU+45wbuQR/qks4c2Aum3YcpDBDq03aYepVnXmKxrBtbPoXyqGXwrnZ85bY+z+cG4ns/l2HK97vDY0iebcUWk4DESGzoJPiYOU1G8nYcBhB6uix/6bVZLKpzDavlnwUeldn3D60oRiPtma+GEo545iurlVWlzPFYumQCJH6w9Rb92/ljG+ewb7xfXRlukqtTYXnn5/1uftHisztylYskOnNhS+2BwsHp9x35sA8Cl7AUy9PvS9m88gonY6wOJdO7vTU7vANwzMj07vDwd8dom9hJ/kUK5kn07co/D1++sGdFe8f3TREdnFX6kVH2SXdFHeMNL24aHT4EJlc3nqfcaFQwLgZXL9xRWsqxu1Hm4pxC+aLAZwMgQnFeNWqa8im0PM7E6ZYLAmmlE6Nqj1M/R8v/kfp85PmnUT+1FNL68zGjgOj9PdUflGdUYyXh0VcM4WqHz5wmFf2dqXmbE7t7sAVeOTg9O5w94sHWbjCrmOLD6F48YmhKfcVBw8THCoQzNKHXAu5gR6CkSL+vuYKxeiB/XTOsV9XsWfPnmjBxkUDNEzdflgVYxG5UESeFpFnReTaCvfnReTW6P7fiMhKm/sBYPxg653YFJPt5JAbhd7yx1ldKjh8GH94GGdOuWA4GGoPUz++5/HS52v71+J0hs5t/NlnZ3ze4KFxfrFlD+edML/i/XNy4QvuwfGpYry0r4OFvdOf4DTqBzwxfJh1femFi3syLq+a083dew5UdIcHBg8zvHechcvth0/PfOMAAHu2H3mc5PB9LwEw9x0npr5m/oTwDdDeWzanfu0kHBjcRd/CxdbX2bFjBwDm0AHra8XEznhsrPm5eaUxWBNjCXtlvgq8FTgd+GMROX3Sw/4U2GeMOQn4e+BztvYDgFeA7Q/D/PRfoOplfHw3z84b4uGTxshm59Hf/zpra5kg4MDtt4Pn0XXOOaXbXbeLkZGZhXMy+8f284PnfsC2Q9u4b/t9pdtX9q3EnRfOjx7++c+nDWk+v2eEv73rKQBed8rCio9Z0rUEgOt+dd2U+8ITnOZWdMbGGD73/Mt4Bt4wP11hfOfieTw5MsZ7N06E4IPA8Pzje7j5Lx8A4OR19oVi9frwTdutn3mQFzcOURj1KLw0zMhDYeg6f/Lc1NfMLuqCjEPhd4fYe+vTePvHMFUelZkWXrHI4IvP079swPpaW7ZsIec6jA3tblhovqOjg3w+z969exuyntJ8bM52Oxd41hjzHICI3AJcDDxZ9piLgb+OPv8u8BUREWPjN/7h/xsewjB+AFa9NvXL18ozz/wNu3b/gEIhDIX1DHu8YsGl5POLrKx36J57ePkTn8Tft4+OV76C7vXrS/ctH3g/z7/wj2zbfhMDy66Y9Vq/3vFrrrr7qim3337x7TjiQEcHnWefzehvf8u+m7/F/CsuB8IBH+/4yv0MDRfYeTB853/eCfN57UkLKq6TdSdyvd/b8j3edfK7jrj/zIG53P3kLg4cLtLXlSUwhk88s52bdgwRAJcsnsf589IV4z85rp9rn9nOj/ccZO39G3nPormccMtL7N8RFqudtG4RvfPtttwA9B/fw4lnLWTro4P84CuPcVqHwykdYfi6+7zjrA2pWPTBV7D7Kxs4/OhuDj+6m57zlzL3osa8yb39i59h1/NbKYyOctKr1s/+hCq59dZb2bFjB11dXeTz+ZIr3blzJ6uXD/DSxgf55sc/wrnvfDenv/YNqa1bCRHhuOOOY8OGDYyMjBAEAb7vc+6553LKKadYXVtpDjbF+HhgW9nX24FXT/cYY4wnIgeAfmBP+YNE5CrgKoDly5fXtpu5y+G0i2DNH8JJf1DbNSyQ7ziO/v430N19IvPnvYbeB74DA79vbb3cwAA9r3st3eefT+9b3oKUzdpdseKDjBzeSlfnyqqu1ZPtYU3/GjYNbWJt/1reuuqtvHnlm1nSvaT0mKWf/xxb3/gmpKz9xHWE5fO7OGVxL69c1scfnLaYgfkzpw5+8Z5fcP2D17O0Z+mU+85ePo+TF/Ww8+AYfV1Zhv2AX+0fpifjcN2Jx/PHx1UOf9dD1hG2vu4M/uWlIX61f5hthSInOULfok5e9bZVnPrqJbNfJCXectVabvz4Lxgf8Yizmp1r+5n3zpOsrZlb1svia85m1z/8FgCxeFTkZPxikXnHHc9r3nMFK195dmrXXbp0KdlsltHRUcbHx+nt7aVQKLBmzRre/ra38VCHy+7nt5amY9nmggsu4Gc/+xmDg4O4rovrunhV9OwrRydiK+wiIu8G3mKMuTL6+grgXGPMfy17zKboMdujr7dGj5lakRKxbt068/DDD1vZs3LsMB4E5EQaNr7QC0xDh0IoSlJE5BFjzLpm70OpjM23eNuB8oTOMmDHdI+R8PDePkCTJErd5B2nYUIMlSeHKYqiVItNMX4IOFlEVolIDrgMuGPSY+4A/kv0+aXAz6zkixVFURSlhbGWM45ywB8Ffkw43ukbxphNIvJp4GFjzB3AjcBNIvIsoSO+zNZ+FEVRFKVVsXpStjHmLuCuSbddV/b5GPBum3tQFEVRlFan/SZwKYqiKEqLoWKsKIqiKE1GxVhRFEVRmoyKsaIoiqI0GWtDP2whIoPAiyldrg9oxPR3G+ukdc16r1Pr85M+bwGTJrMpiWnU77stWmH/jdxD2mutMMZUHgCvNJ2jTozTRES+ZoyZOlz5KFgnrWvWe51an5/0eSLysE4Pqo9G/b7bohX238g9tML3qzSOdg9T33kUr5PWNeu9Tq3Pb9TPXpngaP+Zt8L+G7mHVvh+lQbR1s5YOXpQZ6woyrFMuztj5ejha83egKIoii3UGSuKoihKk1FnrCiKoihNRsVYURRFUZqMirGiKIqiNBkVY+WoRUS6ReQREXl7s/eiKIpSDyrGSssgIt8Qkd0isnHS7ReKyNMi8qyIXFt21/8E/rWxu1QURUkfraZWWgYReR0wDPyLMWZtdJsLPAO8CdgOPAT8MbCUcERmB7DHGPODpmxaURQlBTLN3oCixBhj7hORlZNuPhd41hjzHICI3AJcDPQA3cDpwKiI3GWMCRq4XUVRlNRQMVZaneOBbWVfbwdebYz5KICIvJfQGasQK4py1KJirLQ6UuG2Um7FGPPPjduKoiiKHbSAS2l1tgMDZV8vA3Y0aS+KoihWUDFWWp2HgJNFZJWI5IDLgDuavCdFUZRUUTFWWgYR+Q7wa+BUEdkuIn9qjPGAjwI/Bp4C/tUYs6mZ+1QURUkbbW1SFEVRlCajzlhRFEVRmoyKsaIoiqI0GRVjRVEURWkyKsaKoiiK0mRUjBVFURSlyagYK4qiKEqTUTFWFEVRlCajYqwoiqIoTUbFWFGOEUTkyyLyWxF5VbP3oihKMlSMFeUYQES6gUXA1cDbm7wdRVESomKsHHWIyN+LyDVlX/9YRP6p7OsvicjHUl5zOOXrzRWRD5d9vVJENlb53E4RuVdE3Pg2Y8wIcBxwD/CPIpITkftERI9JVZSjABVj5WjkV8DvAYiIAywA1pTd/3vAL5uwryTMBT4866Mq837ge8YYP75BRPqBLuAQ4BtjCsBPgffUu1FFUeyjYqwcjfySSIwJRXgjcEhE5olIHjgNeFREvi8ij4jIJhG5Kn6yiHxukiv9axH5MxG5XEQeFJENInJDufMse2zFx0TO9ikR+Xq03k9EpDO67y9FZLOI3C0i3xGRjwOfBU6MrvOF6PJupedX4E+A2yfd9hfAF4FNwOnRbd+PHqsoSoujYqwcdRhjdgCeiCwnFOVfA78B1gPrgMcjZ/h+Y8w50W3/LXKPALdwpGP8I+Dh6LbzjTFnAj6ThExETpvlMScDXzXGrAH2A5eIyDrgEuAs4F3RXgCuBbYaY840xvyP6Z4/+XuPznQ+wRjzQtltK6Ofw62Ex0zGUYKNgBZzKcpRgOaTlKOV2B3/HvB3wPHR5wcIw9gQCvAfRp8PEIrdkDHmURFZJCJLgYXAPuAM4BzgIREB6AR2T1rzD2Z5zPPGmA3R548AKwlD6LcbY0YBROTOGb6nSs+fzAJCoS7nM8CnjTFGREpibIzxRaQgIr3GmEMzrKsoSpNRMVaOVuK88RmEDnAb8GfAQeAbInIB8EZgvTHmsIjcA3SUPf+7wKXAEkKnLMA3jTGfmGHN2R4zXva5TyjWkuB7qvT8yYxS9n2IyJmEjvs1IvLV6L4nyh6fB8YS7EFRlCagYWrlaOWXhC08e40xvjFmL2FR1HrCsHUfsC8S4tXAeZOefwtwGaEgf5ew2OlSEVkEICLzRWTFpOdU85jJ3A9cJCIdItIDvC26/RDQm/SbNsbsI8wtx4L8OeAiY8xKY8xK4JVEzjgKyw8aY4pJ11EUpbGoGCtHK08QhmwfmHTbAWPMHuBHQEZEHgf+ZtLjMMZsIhTDl4wxLxtjniQsgvpJ9Jy7CVuFyp8z62MmY4x5CLgDeAz4HmFu+oAxZgj4pYhsLCvgqpafEDrh3we6jTE/LVtvF9AtIvOBNwB3Jby2oihNQIwxzd6DohzTiEiPMWZYRLqA+4CrjDG/reN6ZwEfM8ZcMcvjvgd8whjzdK1rKYrSGDRnrCj2+ZqInE6Yz/1mPUIMEBWg/VxE3PJe43KiquvvqxArytGBOmNFURRFaTKaM1YURVGUJqNirCiKoihNRsVYURRFUZqMirGiKIqiNBkVY0VRFEVpMirGiqIoitJkVIwVRVEUpcmoGCuKoihKk/n/2/VAOZK9F8EAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for b in bands:\n",
    "    plt.plot(Table(data = parse_single_table(FILTERS_DIR + b + '.xml').array.data)['Wavelength']\n",
    "            ,Table(data = parse_single_table(FILTERS_DIR + b + '.xml').array.data)['Transmission']\n",
    "            , label=b)\n",
    "plt.xlabel('Wavelength ($\\AA$)')\n",
    "plt.ylabel('Transmission')\n",
    "plt.xscale('log')\n",
    "plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
    "plt.title('Passbands on {}'.format(FIELD))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### IV.a - Depth overview\n",
    "Then we plot the mean depths available across the area a given band is available"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "vista_j: mean flux error: 2.8761963844299316, 3sigma in AB mag (Aperture): 21.560150523177505\n",
      "vista_h: mean flux error: 4.399040222167969, 3sigma in AB mag (Aperture): 21.098802030896927\n",
      "vista_ks: mean flux error: 5.63737154006958, 3sigma in AB mag (Aperture): 20.829505216595372\n",
      "irac_i1: mean flux error: 97.5140380859375, 3sigma in AB mag (Aperture): 17.734529010508318\n",
      "irac_i2: mean flux error: 124.7995376586914, 3sigma in AB mag (Aperture): 17.466664422123507\n",
      "decam_g: mean flux error: 0.11054158604272202, 3sigma in AB mag (Aperture): 25.098382634409752\n",
      "decam_r: mean flux error: 0.13319610783392877, 3sigma in AB mag (Aperture): 24.895968027296526\n",
      "decam_i: mean flux error: 0.21353520462317221, 3sigma in AB mag (Aperture): 24.383523149421386\n",
      "decam_z: mean flux error: 0.3768913241960511, 3sigma in AB mag (Aperture): 23.76665651226154\n",
      "decam_y: mean flux error: 1.1882132369649205, 3sigma in AB mag (Aperture): 22.519960897684264\n",
      "vista_j: mean flux error: 6.796464920043945, 3sigma in AB mag (Total): 20.626489164127\n",
      "vista_h: mean flux error: 11.229823112487793, 3sigma in AB mag (Total): 20.081264574478517\n",
      "vista_ks: mean flux error: 15.062243461608887, 3sigma in AB mag (Total): 19.762472705379658\n",
      "irac_i1: mean flux error: 24.063777923583984, 3sigma in AB mag (Total): 19.25378733576563\n",
      "irac_i2: mean flux error: 211.563232421875, 3sigma in AB mag (Total): 16.893596378513998\n",
      "decam_g: mean flux error: 0.155469530744103, 3sigma in AB mag (Total): 24.72808364388211\n",
      "decam_r: mean flux error: 0.2068035400430719, 3sigma in AB mag (Total): 24.41830194145927\n",
      "decam_i: mean flux error: 0.3669096121902294, 3sigma in AB mag (Total): 23.79579913964026\n",
      "decam_z: mean flux error: 0.7181806986342703, 3sigma in AB mag (Total): 23.066612540386352\n",
      "decam_y: mean flux error: 2.195031613479145, 3sigma in AB mag (Total): 21.853594914556545\n"
     ]
    }
   ],
   "source": [
    "average_depths = []\n",
    "for b in ap_bands:\n",
    "    \n",
    "    mean_err = np.nanmean(depths['ferr_ap_{}_mean'.format(b)])\n",
    "    print(\"{}: mean flux error: {}, 3sigma in AB mag (Aperture): {}\".format(b, mean_err, flux_to_mag(3.0*mean_err*1.e-6)[0]))\n",
    "    average_depths += [('ap_' + b, flux_to_mag(1.0*mean_err*1.e-6)[0], \n",
    "                        flux_to_mag(3.0*mean_err*1.e-6)[0], \n",
    "                        flux_to_mag(5.0*mean_err*1.e-6)[0])]\n",
    "    \n",
    "for b in tot_bands:\n",
    "    \n",
    "    mean_err = np.nanmean(depths['ferr_{}_mean'.format(b)])\n",
    "    print(\"{}: mean flux error: {}, 3sigma in AB mag (Total): {}\".format(b, mean_err, flux_to_mag(3.0*mean_err*1.e-6)[0]))\n",
    "    average_depths += [(b, flux_to_mag(1.0*mean_err*1.e-6)[0], \n",
    "                        flux_to_mag(3.0*mean_err*1.e-6)[0], \n",
    "                        flux_to_mag(5.0*mean_err*1.e-6)[0])]\n",
    "    \n",
    "average_depths = np.array(average_depths,  dtype=[('band', \"<U16\"), ('1s', float), ('3s', float), ('5s', float)])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def FWHM(X,Y):\n",
    "    \n",
    "    half_max = max(Y) / 2.\n",
    "    #find when function crosses line half_max (when sign of diff flips)\n",
    "    #take the 'derivative' of signum(half_max - Y[])\n",
    "    d = half_max - Y\n",
    "    #plot(X,d) #if you are interested\n",
    "    #find the left and right most indexes\n",
    "    low_end = X[np.where(d < 0)[0][0]]\n",
    "    high_end = X[np.where(d < 0)[0][-1]]\n",
    "    return low_end, high_end, (high_end - low_end)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "data = []\n",
    "\n",
    "for b in ap_bands:\n",
    "    data += [['ap_' + b, Table(data = parse_single_table(FILTERS_DIR + b +'.xml').array.data)]]\n",
    "    \n",
    "for b in tot_bands:\n",
    "    data += [[b, Table(data = parse_single_table(FILTERS_DIR + b +'.xml').array.data)]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ap_vista_j (11670.0, 13380.0, 1710.0)\n",
      "ap_vista_h (15000.0, 17900.0, 2900.0)\n",
      "ap_vista_ks (19930.0, 23010.0, 3080.0)\n",
      "ap_irac_i1 (31754.0, 39164.8, 7410.801)\n",
      "ap_irac_i2 (39980.1, 50052.3, 10072.199)\n",
      "ap_decam_g (4180.0, 5470.0, 1290.0)\n",
      "ap_decam_r (5680.0, 7150.0, 1470.0)\n",
      "ap_decam_i (7090.0, 8560.0, 1470.0)\n",
      "ap_decam_z (8490.0, 9960.0, 1470.0)\n",
      "ap_decam_y (9510.0, 10170.0, 660.0)\n",
      "vista_j (11670.0, 13380.0, 1710.0)\n",
      "vista_h (15000.0, 17900.0, 2900.0)\n",
      "vista_ks (19930.0, 23010.0, 3080.0)\n",
      "irac_i1 (31754.0, 39164.8, 7410.801)\n",
      "irac_i2 (39980.1, 50052.3, 10072.199)\n",
      "decam_g (4180.0, 5470.0, 1290.0)\n",
      "decam_r (5680.0, 7150.0, 1470.0)\n",
      "decam_i (7090.0, 8560.0, 1470.0)\n",
      "decam_z (8490.0, 9960.0, 1470.0)\n",
      "decam_y (9510.0, 10170.0, 660.0)\n"
     ]
    }
   ],
   "source": [
    "wav_range = []\n",
    "for dat in data:\n",
    "    print(dat[0], FWHM(np.array(dat[1]['Wavelength']), np.array(dat[1]['Transmission'])))\n",
    "    wav_range += [[dat[0], FWHM(np.array(dat[1]['Wavelength']), np.array(dat[1]['Transmission']))]]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Depths on SSDF')"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfIAAAFfCAYAAAC1E07gAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xtc1FX+P/DXGxAQGVEBQRHFCyMXb3ih0jUvWdmmuxq1tqnlpdTWtjT1q1lrbv66kJrlupqaXVQqS8tMu2y75mWz2sQLytVLeEdFELnf5vz+mMGIBmaEmWE+4+v5ePh4DJ/LOeczlm/Omc98XqKUAhEREWmTW2MPgIiIiOqPhZyIiEjDWMiJiIg0jIWciIhIw1jIiYiINIyFnIiISMNYyOmmJCKZIjKsscdBRNRQLOTU6ExFtVhE8kXkqojsE5FpImKT/z5F5F0R+X+2aMueRGS+iPwsIgUiclZENlXbFy0i/xKRXNN7lCgivzftGywiBtN5Ved+JCL9arSvRKSw2nFXHX2NRGR7LOTkLEYqpXQAOgB4BcBcAOsad0iOIyKPABgPYJhSyhdAXwD/qXbI5wC+ARAEoDWAJwFcq7b/vOk8HYBbAaQB2Csid9ToqqdSytf0p4V9roaIHImFnJyKUipPKbUNwBgAj4hINwAQES8RWSIip0Xkooi8KSJNTfsGm2ah80Uk2zTDH2vaNwXAWAD/Z5qFfl6tu14ikiQieSKySUS8TecEiMh208w3R0T21rY6ICL9ReQnUxs/iUj/avt2icgiEfnOtNrwLxEJqOXS+wH4Wil1wvQ+ZCml1lSNB0BHAGuVUmWmP98ppf5r5v1TSqmzSqkFAN4CEG/1m09EmsRCTk5JKfU/AGcBDDRtigegB9ALQBcAIQAWVDslGECAafsjANaISFdTMUwA8KppFjqy2jl/AjAcxiLZA8AE0/ZZpr4DYZwBzwfwm2cZi0grADsALAfgD+A1ADtExL/aYQ8BmAjjLNoTwOxaLvkHAA+LyBwR6Ssi7tX2XQFwHMBGERklIkG1tFHTJwB6i0gzK48nIg1iISdndh5AKxERAI8BmKmUylFK5QN4CcCDNY7/m1KqVCm1G8YC+ycL7S9XSp1XSuXAuHTdy7S9HEAbAB2UUuVKqb3KfCjBvQCOKaU2KKUqlFIfwLikXf2XhXeUUhlKqWIAH1Xr41eUUhsB/BXA3QB2A7gkIvNM+xSAIQAyASwFcEFE9ohIuIXrOw9AAFRfQj9gWmm4KiLLLZxPRBrAQk7OLARADowzYx8AiVVFCMBXpu1VcpVShdV+PgWgrYX2s6q9LgLga3q9GMYZ8L9E5GRVQTWjramf6k6Zxm2pj99QSiUopYbBWHinAXhBRO427TurlHpCKdUZxvsICgGsr+viTONQAKrf1NZbKdXC9OdJC+cTkQawkJNTMt1xHQLgvwCyARQDiK5WhPxMN3dVaVljCbk9jDNSwMyyeF2UUvlKqVlKqU4wzq6fNnPTGEztd6ixrT2AczfSn5n+y5VSHwNIAtDNzP4zAP5pbl8NowEcqPELDhG5GBZycioi0lxERgD4EMBGpdQRpZQBwFoAy0Sktem4kKrZajV/FxFPERkIYASAj03bLwLodANjGCEiXUxL+tcAVJr+1PQFAL2IPCQiHiIyBkAUgO3WX/H1PieIyL0iohMRNxG5B0A0gB9FpKWI/N00JjfTzW+TYPxcvWY7YnpvngfwKIyf7xORC2MhJ2fxuYjkAzgD4FkYbxybWG3/XBiXu38QkWsA/g2ga7X9WQByYZwlJwCYppRKM+1bByDKtCy/1YqxhJvaLwDwPYCVSqldNQ9SSl2B8ReGWTDekPZ/AEYopbKtuuJfuwZj0T0N41L4qwAeN92ZXgYgzDSmawCOAijFLzfnAUBbESkwjfknAN0BDFZK/aseYyEiDRHz9/AQaYeIDIZx9t6uscdCRORonJETERFpGAs5ERGRhnFpnYiISMM4IyciItIwFnIiIiIN82jsAVgjICBAhYWFNfYwiIg0JTExMVspFWj5yFrPb+3h4fEWjA8f4sSvcRgAHK2oqHi0T58+l8wdoIlCHhYWhv379zf2MIiINEVEaj5C+IZ4eHi8FRwcHBkYGJjr5ubGG6oagcFgkMuXL0dlZWW9BeAP5o7hb1hERFSbboGBgddYxBuPm5ubCgwMzEMdj2RmISciotq4sYg3PtPfQa31moWciIhIw1jIiYjopjRmzJgOiYmJ3rXtX758uX9mZmaT+rQ9Y8aMtlu3btXVf3TW08TNbkRERLa2adOmOm8G3LhxY0CvXr2Kw8LCym+07ddff/285aNsg4WciIgsmrP5cGhGVr6PLdvUB+uKFt/f84yl44YNG9b5woULnqWlpW7Tpk27OHv27GwfH5+YsWPHXv7uu+90fn5+lVu2bDnZtm3biprnHjhwwPuRRx7peOTIkVQASE9P9xw5cmSXjIyMlNjY2K5Lliw5079//6IxY8aEJSUlNRMRNXbs2Oz27duXHz161Ofhhx/u5O3tbdi/f3/qwoULg7/66qsWpaWlbn379i1ISEg45eZmfmE7Li4ubMSIEXkTJ07MbfAbZQGX1omIyKklJCRkJicnpx46dChl9erVQVlZWe7FxcVuvXv3LkpJSUkdMGBA/rx589qaO7d3794l5eXlkpKS4gkA69evbzVq1KhfFdfvv//e58KFC02OHTuWnJGRkTJ9+vQrEydOzO3WrVvR+vXrT6alpaX4+vqqOXPmXDp69GjqsWPHkouLi90+/PBDP0dcvyUuPSP/9t01uHTqZGMPA607dMKQCVMaexhERPVmzczZXuLj44N27NjRAgCysrKaJCcne7u5ueHRRx/NAYBJkyZdue+++7rUdv6oUaNyNm7c2Oqll17K+vTTT1tu2rTpV4UhIiKi9MyZM16PPPJI6MiRI/NGjx59zVw7X375pe61114LLikpcbt69apHVFRUMYA8G15qvXBGTkRETmv79u263bt36/bv35+Wnp6eEhkZWVxcXPyb2iUitbYxfvz43K1bt7ZMSkryEhF07969tPr+wMDAyqNHj6YMGTIkf+XKla0ffPDBsJptFBUVyaxZszp88sknJzIyMlLGjRuXXVJS4hQ11KVn5JwFExFp29WrV939/PwqdTqd4eDBg96HDx9uBgAGgwHvvPNOyylTpuS+++67/rGxsfm1tREdHV3q5uaGBQsWtB09enROzf0XLlzw8PLyMkyYMOGqXq8vnTRpUkcA8PX1rczLy3MHgKKiIjcACA4OrsjLy3P7/PPPW44cOdLun39bw6ULORERaVtcXFzemjVrAvV6fVTnzp1LevbsWQgATZs2NSQnJzeNjo4O1ul0lZ988kmdn6Ped999OYsWLWoXHx9/rua+zMzMJpMnTw4zGAwCAC+88MJZAHj44Yez//rXv3aYM2eOYf/+/aljx469HBUVFd2uXbuyqnHURUQc8jAdTeSR9+3bV7nSs9avfn4CZect/jdgU55tm6HFyM4O7ZOIGpeIJCql+tb3/MOHD2f27Nkz25ZjshUfH5+YoqKig409jtoMHTq0y8yZMy+OHDmy1pWCG3H48OGAnj17hpnb5xTr+0RERK7igQceCCsuLna76667ChzRH5fWGwFnxkREDWNuNj5+/Pj2P/30k2/1bY8//vjFp5566oq9xtEYfdbEQk5ERC5hw4YNp2+GPmty6ULO75ETEZGr42fkREREGma3GbmIhAJYDyAYgAHAGqXUG9X2zwawGECgUsoud0VyFkxERK7OnkvrFQBmKaUOiIgOQKKIfKOUSjEV+TsBNPpnC0RERFpmt6V1pdQFpdQB0+t8AKkAQky7lwH4PwDO/yV2IiJySfbMI4+Nje26Z88em6bF1cYhn5GLSBiAGAA/isgfAJxTSh12RN9ERETmbNq06VSfPn1Katu/cePGgNOnT9erkDuS3e9aFxFfAFsAzIBxuf1ZAHdZcd4UAFMAoH379vYcIhERWbJ1eigupdh2htk6qgij/umSeeQA8MEHH7ScPn16h/z8fPc333wzc/jw4XZ5QIxdZ+Qi0gTGIp6glPoEQGcAHQEcFpFMAO0AHBCR4JrnKqXWKKX6KqX6BgYG2nOYRETkxLSaR15RUSFHjhxJjY+PP/PCCy+YHZ8t2POudQGwDkCqUuo1AFBKHQHQutoxmQD62uuudWfFZ60TkeZYMXO2F63mkT/wwAO5ANC/f//COXPmeNbj0q1iz6X1AQDGAzgiIodM2+Yrpb6wY59UQ05ODsrKypB/tRQZ2Xts3n5wcDDuuecem7dLRAT8Oo9cp9MZYmNju9Ynj/yBBx7o9OCDD+bWlUf+6aefNl+5cmXrTZs2tfr4448zqx9TlUf+448/pnTp0qX86aefbmspj9zb21sBgIeHByorK2sfYAPZrZArpf4LoM6BK6XC7NW/M3PkzDjxyy+RlZXlsP6IiGyJeeSWufQjWgmcLRORpmk5j9xRmEdOROSimEfuOphHTkRE5KK4tE5ERJrDPPJfsJATEZFLuFnzyLm0TkREpGEs5ERERBrGQk5ERKRhLOREREQaxkJOREQ3JVfJI+dd60REdFPatGnTqbr2b9y4MaBXr17FYWFh5Y4aU32wkBMRkUV/++5vocdzj9t0htmlZZeiRQMWuWweOQBUVlbigQceCGvXrl3Za6+9dr5mP88///yler+BJizkLu5LO4emMP2MiOwtISEhMygoqLKgoEBiYmKixo0bl1uVR7527dqzs2fPbjNv3ry269ev/813uqvnkUdFRZVZyiMHgOzsbPeAgIDKVatWtV6yZMmZ22+/vQgA5syZc2nJkiUXAGDUqFEdP/zwQ7+HHnqo1hjT8vJyGTVqVMeoqKji+Pj4rL179/6mH1u8PyzkVC/FxadQWVmES5cPIfHABpu0qfONhF7/N5u0RUS2Zc3M2V60mkf+l7/8pcOoUaNy4uPjs26knxvFQu7i7DVbzshYhPyCs3Zpm4ioipbzyPv27Vuwd+/e5kVFRRd9fHyUNf3UBws51QtnzkTkCFrOI586dWr2zp07dSNGjOj89ddfH8/OzjbbT0OxkBMRkdPSeh75woULL86cOdP9vvvu6zh//vwsc/00FPPIiYhcFPPIXQfzyImIiFwUl9aJiEhzmEf+CxZyIiJyCcwjJyIiIs1hISciItIwFnIiIiINYyEnIiLSMBZyIiJyaTNmzGi7detWnb3aj4mJiah6PXDgwHCdTtdryJAhtT773dZ41zoREbm0119//by57RUVFfDwaHgZPHjwYFrV69mzZ2cVFha6rV27NrDBDVuJhZyIiCw6P//Z0NJjx2yaR+4VHl7U9qUX7ZpHDgBxcXFhI0aMyJs4cWJuSEhI9z//+c/Z3377bfOpU6deys/Pd3/nnXcCy8vLJSwsrHTz5s0/63Q6w5kzZzwmTZrU4fTp014AsGLFilN33nmn2ceyVn/K3B//+Mf87du32232bw6X1omIyKklJCRkJicnpx46dChl9erVQVlZWe5VeeQpKSmpAwYMyJ83b15ba9vz9vY2JCYmpk+ZMiV37NixuUePHk1NT09P6dq1a/Hy5csDAGDatGntBw4cmJ+enp6SnJyc0rt37xL7XWHDcEZOREQWWTNztpeG5pHX9PDDD19PLUtMTGy6YMGCkPz8fPfCwkL3QYMG5QHAvn37dJs3b/4ZADw8PODv719p26uyHRZyqhdjHnmqTdvU+UYyHpWIfsUWeeQ16XQ6Q9XrKVOmdNy8efPx2267rXj58uX+u3fvduiyuC1waZ2IiJyWpTxyALCUR16XoqIit/bt25eXlpbKhx9+2Kpq+4ABA/IXL14cCBhvisvJyXHaeskZOdULZ85E5Ai2yiOvzbx5887HxsZGhoSElEVGRhYVFBS4A8CqVatOT5gwoYNerw9wc3PDihUrTg0bNsxiBnmfPn26njx50ru4uNg9KCiox8qVKzPj4uKu1Wds1mIeORGRi2IeuetgHjkREZGL4tI6ERFpjqPzyLOystwHDx7cteb2Xbt2pQcHBzfqHe0s5ERE5BLsmQ0eHBxcmZaWlmKv9huCS+tEREQaxkJORESkYXYr5CISKiLfikiqiCSLyFOm7YtFJE1EkkTkUxFpYa8xEBERuTp7zsgrAMxSSkUCuBXAdBGJAvANgG5KqR4AMgA8Y8cxEBERuTS7FXKl1AWl1AHT63wAqQBClFL/UkpVJdT8AKCdvcZARETkqDzyffv2Ne3Vq1dEly5dovV6fdTatWtb2qvP6hxy17qIhAGIAfBjjV2TAGxyxBiIiOjm5Kg8cl9fX8OGDRt+7t69e2lmZmaTfv36RY4ePfpaQECAXb+eZvdCLiK+ALYAmKGUulZt+7MwLr8n1HLeFABTAKB9+/b2HiY5yN+OncXRguI6j+nm2xSLwrlQQ+RM/rM+NTTnXIFN88hbhfgW3fFwpMvkkffo0aO0altYWFh5q1atKi5cuOBh70Ju17vWRaQJjEU8QSn1SbXtjwAYAWCsquUZsUqpNUqpvkqpvoGBgfYcJhEROTEt5pF/++23PuXl5RIVFVVq+eiGsduMXIyZcusApCqlXqu2fTiAuQAGKaWK7NU/OSfOtIm0yZqZs71oLY/81KlTTSZOnNhp3bp1P7u7u1t7Wr3Zc0Y+AMB4AENF5JDpz+8BrACgA/CNadubdhwDERFpWPU88vT09JTIyMhiW+eRr1ix4nRGRkbK3Llzz5eWljaoLubk5Ljdc889XRYsWHDujjvusJiWZgv2vGv9v0opUUr1UEr1Mv35QinVRSkVWm3bNHuNgYiItE1LeeQlJSVy7733dnnwwQevTJo0KdfS8bbCJ7sREZHTiouLy6uoqBC9Xh81f/78tmbyyCP37Nmje/nlly/Up/2qPPKBAwfqw8PDr38OvmrVqtO7d+/W6fX6qG7dukUdOHCgqaW23n777ZY//fST7/vvvx8QERERFREREbVv3z6L5zUU88iJiFwU88hdB/PIiYiIXBRjTImISHOYR/4LFnIiInIJzCMnIiIizWEhJyIi0jAWciIiIg1jISciItIwFnIiIrpppKene4aHh0c39jhsiYWciIhIw/j1M3Io5pETadPXq14PzT5zyqZ55AGhHYrufnyG3fPI9+7d6/Poo4+GNW3a1HDLLbcUVG2vqKjA9OnT23333Xe6srIyeeyxxy7NmTMnGwCee+65oI8++shfRHDHHXfkrVy58tzSpUsDzGWXx8XFhXl7exuOHz/ufe7cOa/Vq1f//O677wYkJiY2i4mJKdyyZUtmbde2bNmygDfeeCO4devW5Z06dSrx9PRU69evv6Gv0XFGTkRETq2heeSTJ08Oe+21104fOnQorfr2119/PcDPz6/y6NGjqYcPH0597733AtPS0jw/+uij5jt27GiZmJiYlp6envL8889nAUBt2eUAkJeX5/H9999nvPLKK2fGjBkTPmfOnIvHjh1LTktLa1rb89YzMzObLFmypM2PP/6Yunfv3oxjx4551+f94YycHIozbSJtsmbmbC8NySO/cuWKe35+vvu9995bUHXszp07/QDg3//+d/O0tDSfbdu2tQSA/Px895SUFO9vvvmm+bhx47Kr4k6DgoIqgdqzywHg3nvvverm5obevXsX+fv7l8fGxhYDgF6vLz5x4oRX//79f7MUuXfv3ma33HJLflX7o0ePzs3IyLjhYs5CTkRETqt6HrlOpzPExsZ2vZE8cqVUXftk6dKlp+Pi4q5V3/7FF180N3fOlClTOm7evPn4bbfdVrx8+XL/3bt366r2eXt7KwBwd3eHp6fn9TQyNzc3VFRUmB2ArULLuLROREROq6F55AEBAZW+vr6VX3/9ta/p2OuZ43feeWfeqlWrAktLSwUAkpKSvK5du+Y2fPjwaxs2bAjIz893A4CLFy+6A7Vnl9fXwIEDC3/88Ufd5cuX3cvLy/HZZ5+1rE87nJETEZHTiouLy1uzZk2gXq+P6ty5c4mZPPJgnU5X+cknn5ysrY1169ZlVt3sNnTo0Ouz75kzZ2ZnZmZ6de/ePVIpJa1atSr/4osvTtx///3XDhw44NOrV6/IJk2aqGHDhuWtWLHiXFV2eUhISFlkZGRRQUGBe0OurWPHjuUzZ8680K9fv8jWrVuX6/X6Yj8/vxsOYGEeORGRi2IeufPLy8tz8/PzM5SXl+Puu+/uMmHChOyHH374as3j6soj54ycXM7fP09Gyvlrlg+0IKptczw/0qWeG0FETmbOnDlt9+zZ07y0tFQGDRp0bdy4cb8p4pawkBMRkeY4Oo+8oXr06BFRVlb2q/vS1q9f//OaNWvONrRtqwq5iHgBiAMQVv0cpdQLDR0Aka1xFk10c7JnHnlDJSUlpVk+qn6snZF/BiAPQCKAUnsNhoiIiG6MtYW8nVJquF1HQkRERDfM2u+R7xOR7nYdCREREd2wOmfkInIEgDIdN1FETsK4tC4AlFKqh/2HSERERLWxtLQ+wiGjICIicoD09HTPESNGhB87diy5scdiK3UWcqXUKQAQkQ1KqfHV94nIBgDjzZ5IREREN6SiogIeHjf+rXBrz/jV93lExB1AnxvujYiINClnc0ZoeVahTfPImwQ3K2p1v/6mziP38fGJmTJlysWdO3c2X7x48dm77767oLZja1PnzW4i8oyI5APoISLXRCTf9PMlGL+SRkREZFeumkcOAMXFxW7dunUrTkpKSqtPEQcsL62/DOBlEXlZKfVMfToguql8OQ/IOmL7doO7A/e8Yvt2iaxkzczZXlw1jxwwxp5OmDAhtyHvj7VL6/NF5D4Av4PxLva9SqmtDemYiIjIElfOIwcAT09PQ30+F6/O2rP/CaALgA9MP08TkTuVUtMb1DuRq+GsmcimLOWRT5kyJdfaPPK77767wFwe+YgRI/K9vLxUUlKSV1hYWPnw4cOvvfjii20fe+yxHJ1OZ7h48aJ7UFBQZc088jZt2pQ76n2oi7WFfBCAbsqUeSoi7wGww/ohERHRL1w5j9xWrMojF5FPAMys9nW0DgBeUUr92c7jA8A8ciKi+mAeueuwRR65P4BUEfmf6ed+AL4XkW0AoJT6Q4NHSURERDfM2kK+wK6jICKbi/9fPNJyGp6cGNEqAnNj59pgRES24yp55FV3tzeEVYVcKbXbtJwerpT6t4g0BeChlDJ7cwEREZGjMY+8DiLyGIApAFoB6AygHYA3Adxhr4ERUcNwFk10c7A2xnQ6gAEArgGAUuoYgNb2GhQRERFZx9pCXqqUKqv6QUQ8YHwwDBERETUia2922y0i8wE0FZE7AfwFwOd1nSAioQDWAwgGYACwRin1hoi0ArAJQBiATAB/Uko16PF0RNX9/fNkpJy/ZvlAC6LaNsfzI6MtH0hEmuFsMaaDBg3qsmXLlp8DAgIq69uGtTPyeQAuw/gQmKkAvgDwnIVzKgDMUkpFArgVwHQRiTK19R+lVDiA/5h+JiIiuuns3r37eEOKOGD9XesGEdkKYKtS6rKV51wAcMH0Ol9EUgGEAPgjgMGmw94DsAsA78ohm+Esmsi1uHKMaUhISPf9+/entmnTxuzYrVFnIRfjU+OfB/AEADFtqgTwD6XUC9Z2IiJhAGIA/AggyFTkoZS6ICJmb5oTkSkw3imP9u3bW9sVERHZwdatW0MvXbpk0zzy1q1bF40aNcpiqlpCQkJmUFBQZUFBgcTExESNGzcutyrGdO3atWdnz57dZt68eW3Xr19v9utnkydPDlu2bNnpe++9t2Dq1KntqrZXjzEtLi6Wfv36RYwcOfJaUlKSd1WMadWz1gFjjOmsWbOyAeDJJ59su3z58oBnn332EvBLjOn777/fYsyYMeE7d+5M69OnT3GPHj0i9+3b17S29DNbsLS0PgPGu9X7KaX8lVKtANwCYICIzLSmAxHxBbAFwAyllNUfXCql1iil+iql+gYGBlp7GhERuZj4+Pigrl27RvXp0yeythjT//3vf77mzjUXY1q179///nfzjz76yD8iIiIqJiYmMjc318NSjGmfPn266vX6qC1btvgnJyd7V7VlLsbU3d39eoypPd8fS0vrDwO4Uyl1/Vm7SqmTIjIOwL8ALKvrZBFpAmMRT1BKfWLafFFE2phm420AXKr/8ImIyBGsmTnbg6vHmNqCpRl5k+pFvIrpc/ImdZ1oWpZfByBVKfVatV3bADxiev0IgM+sHy4REd1MLMWYAoC1MaamY38TY1paWioAkJSU5HXt2jW34cOHX9uwYUNAfn6+GwBULa3XjDG175Vbz9KMvKye+wDjkvx4AEdE5JBp23wArwD4SEQmAzgN4AFrBkpERDcfxphaVmeMqenGtkJzuwB4K6XqnJXbCmNMiYhuHGNMXUe9Y0yVUk7x2wYRERGZZ+2T3YiIiJwGY0x/wUJOREQu4WaNMbX2Ea1ERETkhDgjJ7KlL+cBWUds325wd+CeV2zfLhFpHmfkREREGsYZOZEtcdZMRA7GGTkREd000tPTPcPDw10qIpGFnIiISMO4tE5ERBalpM4NLSzIsGmMaTNffVFUZLzFMBZXzSNPSEjwW7RoUQgAlJSUuJWXl8u5c+du+G5ZzsiJiMipJSQkZCYnJ6ceOnQoZfXq1UFZWVnuVXnkKSkpqQMGDMifN29e29rOnzx5cthrr712+tChQ7/6Lnf1PPLDhw+nvvfee4FpaWmeH330UfOqPPL09PSU559/Pgsw5pEfPXo0NT09PaVr167Fy5cvD6hqqyqP/JVXXjkzZsyY8Dlz5lw8duxYclpaWtN9+/Y1NTeusWPH5qWlpaWkpaWlREVFFT3xxBNZ9Xl/OCMnIiKLrJk520t8fHzQjh07WgBAbXnk9913Xxdz55rLI9+5c6cfYMwjT0tL89m2bVtLAMjPz3e3lEe+YMGCkPz8fPfCwkL3QYMG5VX1Yy6PHMD1PPL+/fvX+gS35557Lsjb29vwzDPPXK7P+8NCTuSi4v8Xj7Schj9MKqJVBObGzrXBiIhunKvnkX/22We6rVu3tvrhhx/q/T8rl9aJiMhpuXIeeUZGhudTTz3VYfPmzSft9gPiAAAbRElEQVR8fX1rjyK1gDNyIhfFWTS5AlfOI1+9erV/Xl6e+6hRo7oAQFBQUNnu3buP32g7deaROwvmkRMR3TjmkbuOuvLIubRORESkYVxaJyIizWEe+S9YyImIyCUwj5yIiIg0h4WciIhIw1jIiYiINIyFnIiISMNYyImI6KbhinnkvGudiDQr66WXUJpqt5uBr/OKjEDw/Pl274eoPljIiYjIohmpp0PTCktsmkce0cy76PXI9jdtHvmyZcsCjh492nTdunVnAGDp0qUBqamp3m+99dbZG3kfWciJSLM4S745JCQkZAYFBVUWFBRITExM1Lhx43Kr8sjXrl17dvbs2W3mzZvXdv369Wa/Rz558uSwZcuWnb733nsLpk6d2q5qe/U88uLiYunXr1/EyJEjryUlJXlX5ZHrdDpDVWjK2LFjc2fNmpUNAE8++WTb5cuXBzz77LOXgF/yyN9///0WY8aMCd+5c2danz59inv06BG5b9++puZiTCdPnpwTHR0dVVpaetbLy0tt3LgxYPXq1adu9P1hISciIousmTnbi6vmkTdv3twwYMCA/E2bNvl17969pLy8XOrzpDcWciIiclqunkc+ZcqU7BdffDFYr9eXjBs3rl4BNbxrnYiInJYr55EDwNChQwsvXLjg+emnn/pPnjw5pz5tcEZOREROy5XzyKuMGjUqNykpyScwMLCyPuczj5yIyEUxj1wbhgwZ0mXGjBkX//jHP5pdVQCYR05EROR0srOz3cPCwrp5e3sb6irilnBpnYiINMdV8sgzMzOPNrRtFnIiInIJzCMnIiIizWEhJyIi0jAWciIiIg1jISciItIwuxVyEXlbRC6JyNFq23qJyA8ickhE9otIrL36JyIi15SZmdlk+PDhnWrbn52d7f7KK68E1rf9mJiYiPqe2xjsOSN/F8DwGtteBfB3pVQvAAtMPxMREVktLCys/Kuvvqr1SW5XrlxxX7duXev6tn/w4EH7h9zbkN2+fqaU2iMiYTU3A2hueu0H4Ly9+iciItuZs/lwaEZWvk3zyPXBuqLF9/esM1Xt8ccfD+nQoUPZvHnzLgPA008/3Van01W+//77AceOHUvev3+/98SJEzuWl5eLwWDAli1bTjzzzDMhZ86c8YqIiIgaNGjQtVdfffX88OHDu+Tl5blXVFTIggULzo8bN+5qbX1q7alxjv6MfAaAxSJyBsASAM/UdqCITDEtv++/fPmywwZIRETOY9y4cTlbtmy5HlDy2Weftbz11lsLq37+xz/+EfiXv/zlYlpaWkpSUlJqx44dy5YuXXo2NDS0NC0tLWX16tVnfXx8DDt27DiekpKSunv37oz58+e3MxgMjXNBduDoB8I8DmCmUmqLiPwJwDoAw8wdqJRaA2ANYHzWuuOGSERENVmaOdvLgAEDiq9cueKRmZnZ5MKFCx5+fn6VnTp1Kqvaf9tttxUuWbKkzdmzZz0ffPDB3O7du5fWbMNgMMiMGTPa/fDDD75ubm64dOmS59mzZz3at29f4dirsQ9Hz8gfAfCJ6fXHAHizGxER1WnkyJG5GzdubJmQkNAqLi7uV1Gf06ZNy/nss8+ON23a1HDPPffot23bpqt5/urVq1tduXLF48iRI6lpaWkp/v7+5eYyzbXK0TPy8wAGAdgFYCiAYw7un4iINGb8+PE5jz32WFhubq7H7t2700tKSqRqX0pKimdkZGRpdHT0pZMnT3odOnSoaWxsbFFhYeH1Qp2Xl+ceEBBQ7uXlpT7//HPd+fPnPRvnSuzDboVcRD4AMBhAgIicBfA8gMcAvCEiHgBKAEyxV/9EROQa+vbtW1JYWOgWFBRU1qFDh/L09PTrhXjDhg2tPv74Y38PDw8VGBhY/vLLL58PCgqq7NOnT0F4eHj00KFD8xYuXJh1zz33dOnWrVtkdHR0UceOHUvq6k9E6trtdJhHTkTkolw5j9xesrKy3Hv37h11/vz5I409luqYR05ERGRBZmZmk1tvvTVy+vTpFxt7LDeCMaZERHTTycrKch88eHDXmtt/+OGH1ODg4MrGGFN9sZATEdFNJzg4uDItLS2lscdhC1xaJyIi0jAWciIiIg1jISciItIwFnIiIiINYyEnIiJNsWceeXp6umd4eHh0/UfneCzkRESkKfbOI9cafv2MiIgs2zo9FJdSbJpHjtZRRRj1T6fLI6+srMSDDz7YYf/+/b5BQUFlX3/99XFfX1+nfQwqZ+REROS0GiOP/PTp095PPvnkpePHjyf7+flVrl+/vqWdL7NBOCMnIiLLLMyc7aUx8shDQkJK+/fvXwwAMTExRZmZmV72u8KG44yciIicmqPzyD09Pa8vo7u7u6uKigqnjkPjjJyINCvrpZdQmppm9368IiMQPH++3fsh85hHXjfOyImIyKnVzCOvvm/Dhg2t9Hp9dERERNSxY8e8p06deiU4OPh6HvnUqVPbPfroozmHDx9u1q1bt8iNGze2spRHrjXMIyciclHMI3cdzCMnIiJyUfyMnIiIbjq15ZHv2rUrnXnkRERETs6V8shZyImIrLT3owxknylwaJ8Bob4Y+Ce9Q/skbeFn5ERERBrGGTkRkZU4MyZnxBk5ERGRhrGQExGRpjCP/NdYyImISFOYR/5r/IyciIgs+tt3fws9nnvcpnnkXVp2KVo0YJHT5ZFXSUlJ8YyLi+vy5ptvZjZr1sxQsx9zSWuNgTNyIiJyWo2RRw4Ahw8f9oqLi+uybt26nwcNGlRkrh87XvYN4YyciIgssjRztpfGyCPPycnxGDVqVJePP/74RN++fUus7aexcEZOREROzdF55DqdrrJNmzZlu3bt8r2RfhoLZ+REROTUHJ1H3qRJE/XVV1+dGDJkSLivr69h2rRpOeb6+cMf/pBvz+u2Fgs5ERE5tZp55Onp6dcL8YYNG1p9/PHH/h4eHiowMLD85ZdfPh8UFHQ9j3zo0KF5CxcuzLrnnnu6dOvWLTI6OrrImjzy5s2bG77++uvjgwcP1vv6+hqSk5O9a/Zj36u2HvPIiYhcFPPIXQfzyImIiFwUl9aJiKzE9DPXwTxyIiIiDWMeORHRTYgzY3JG/IyciIhIw1jIiYiINIyFnIiISMPsVshF5G0RuSQiR2ts/6uIpItIsoi8aq/+iYhI+2JiYiLs1XZCQoLf/PnzgwHgyy+/9I2Kior08PDo884777S0V5/2YM+b3d4FsALA+qoNIjIEwB8B9FBKlYrITZMXS0REN+7gwYNpNbdVVFTAw6Ph5Wvs2LF5APIAoFOnTmXvvPNO5iuvvBLU4IYdzG6FXCm1R0TCamx+HMArSqlS0zGX7NU/ERHZzvn5z4aWHjtm0zxyr/DworYvvVhnqpqPj09MUVHRwe3bt+sWLVrUpnXr1uUpKSk+J06cSB42bFjnCxcueJaWlrpNmzbt4uzZs7MBYPPmzc0XLFgQUllZKa1atar4/vvvM8y1vXz5cv/9+/c3W79+/emuXbuWAYCbm/Y+cXb018/0AAaKyIsASgDMVkr9ZO5AEZkCYAoAtG/f3nEjJCIip5SUlNTs4MGDyREREWUAkJCQkBkUFFRZUFAgMTExUePGjcs1GAzyxBNPhO3atSstIiKi7OLFi+6NPW57c3Qh9wDQEsCtAPoB+EhEOikzD3xXSq0BsAYwPmvdoaMkIqJfsTRzdoQePXoUVhVxAIiPjw/asWNHCwDIyspqkpyc7H3x4kWP2NjY/KrjgoKCNPWUtvpwdCE/C+ATU+H+n4gYAAQAuFzXSYmJidkicsqK9v1g+rzDAezVly3btUVb9W2jPucFAGBAg2058v8Je3Km69DSvzMdbDUQZ+Dj42Ooer19+3bd7t27dfv370/T6XSG2NjYrsXFxW5KKYhIXc24HEcX8q0AhgLYJSJ6AJ6w4h9upVSgNY2LyBql1JSGDdE69urLlu3aoq36tlGf80Rkf0OSmui3HPn/hD0503W4wr8zruDq1avufn5+lTqdznDw4EHvw4cPNwOAIUOGFM6aNatDWlqaZ9XSuqvPyu359bMPAHwPoKuInBWRyQDeBtDJ9JW0DwE8Ym5ZvQE+t2FbjdWXLdu1RVv1bcORfxdUO1f5e3Cm63CFf2c0Ly4uLq+iokL0en3U/Pnz2/bs2bMQANq2bVuxfPnyzNGjR3fp2rVr1OjRoztZ097u3bt9goKCenzxxRctZ86c2aFLly7R9r0C29FEHjndHDgjJ3IuzCN3HswjJ61Y09gDICLSGqafkdMwfVOBiMim3njjDf9Vq1b96kEv/fr1K9iwYcPpxhqTLXFpnYiIzOLSuvPg0joREZGLYiEnIiLSMBZycmoi0kxEEkVkRGOPhYjIGbGQk0PVEW873BRve1xE5lXbNRfAR44dJRGRdrCQk6O9C2B49Q0i4g7gnwDuARAF4M8iEiUiwwCkALjo6EESkXNwVB75woULgzp37hyt1+ujbrvtNn1GRoanvfq1NX79jByqlnjbWADHlVInAUBEPoQxt94XQDMYi3uxiHyhlDKAiG4ajsoj79OnT9GsWbNSdTqdIT4+PnDmzJntduzYcbLBnTgACzk5gxAA1ZOVzgK4RSn1BACIyAQA2SziRI3nP+tTQ3POFdg0j7xViG/RHQ9HOkUe+ciRI/Ortv/ud78r2LRpk78tr9WeWMjJGZiLKrr+gAOl1LuOGwoROStH5ZGvXr06cNiwYc6StmcRCzk5g7MAQqv93A7A+UYaCxGZYWnm7AiOyCNfuXJlq8OHD/usXr063bajtx8WcnIGPwEIF5GOAM4BeBDAQ407JCJyNvbOI9+6datuyZIlbfbu3ZvetGlTzTz2lHetk0OZi7dVSlUAeALA1wBSAXyklEpuzHESkXOrK4/8xx9/1KWlpXkCgLVL6999913Tv/71rx0+++yz4yEhIRX2HLutcUZODqWU+nMt278A8IWDh0NEGhUXF5e3Zs2aQL1eH9W5c+cSc3nkBoMB/v7+5fv27Ttmqb05c+aEFhUVuT/wwAOdTe2U7dy587i9r8MWGJpCRERmMTTFeTA0hYiIyEVxaZ2IiFwa88iJiOimxKV158GldSIiIhfFQk5ERKRhLOREREQaxkJORESkYSzkRBomIv8QkQMi0q+xx0LkCE8//XTbBQsWBFk+8ubBQk6kUSLSDEBrAFMBjGjk4RBRI+H3yEkTRGQZgFNKqddNP38N4IxS6lHTz0sBnFNKvWbDPguUUr42bK8FgIeUUitNP4cB2K6U6mbFuU0BfAVgqFKqEgCUUoUi0gbALgDtRcQTwL9Nx2jqWdHk/L5e9Xpo9plTNs0jDwjtUHT34zMspqrNnTs3eNOmTQFt27Yt8/f3L4+JiSlKTk72mjZtWvucnBwPb29vw1tvvXUqJiam5MyZMx6TJk3qcPr0aS8AWLFixak777yzsLbsch8fn5hHHnnk0p49e5r7+flVvvjii2fnzp0bev78ec/4+PjTY8eONRtnmp+f7zZmzJiw48ePe4eHh5ecOXPGc8WKFadvv/32Ilu+R9bgjJy0Yh+A/gAgIm4AAgBEV9vfH8B3jTCuG9ECwF/qee4kAJ9UFXEAEBF/AD4A8gFUKqXKAPwHwJiGDpTIWezdu9fn008/bXXkyJGU7du3H68KR3n00Uc7rFy58nRycnLq4sWLzz7++OPtAWDatGntBw4cmJ+enp6SnJyc0rt37xLAmF2enJyceujQoZTVq1cHZWVluQNAcXGx25AhQ/KTk5NTmzVrVvncc8+F7N27N+Pjjz8+vmjRopDaxrV48eLAFi1aVGZkZKQsXLjwfEpKSjNHvB/mcEZOWvEdgGWm19EAjgJoIyItARQBiARwUES2wpht7g3gDaXUGgAQkXgYZ/RVs+GFMBbAiwCeBOAJ4EcAf6leLE3HjjN3jGlG/SWA/8L4i8Q5AH9UShWLyN8AjAVwBkA2gEQAfQF0FpFDAL4B8E8A7iKytub5Zq5/LH4b7focgCUAHgMQBeMvO1sBvAwgweI7SnQDrJk528O3337r+/vf//6qTqczAMBdd911taSkxO3gwYO+VQEnAFBWViYAsG/fPt3mzZt/BgAPDw/4+/tXAuazy4ODgwubNGmi7r///msAEB0dXezl5WXw8vJSsbGxxefOnfOsbVz79u3zfeqppy4BQL9+/Ur0er3DZ+JVOCMnTVBKnQdQISLtYSx638NYVG+DsUAmmWakk5RSfUzbnjTNWgHgQ/x6pvonAPtN2wYopXoBqISxYF4nIpEWjgkH8E+lVDSAqwDiRKQvgDgAMQDuM40FAOYBOKGU6qWUmlPb+TWv3bRk3kkplVltW5jpfdgEY/Rr1erEUQC88Y1cSs18cYPBAJ1OV5GWlpZS9efkyZO1Rh9Xzy5PT09PiYyMLC4uLnYDAA8PD+XmZiyFbm5u8PLyUgDg7u6OysrKWoPNnempqCzkpCXfwVi8qgr599V+3mc65kkROQzgBxhn5uEAoJQ6CKC1iLQVkZ4AcgF0B9AHwE+mWfIdADrV6PMOC8f8rJQ6ZHqdCCAMwO8AfKaUKlZK5QP4vI5rMnd+TQEwFvnq/h+AF5TxX5Prhdy0mlAmIro6+iTSjKFDhxbs2LGjRUFBgeTm5rp98803LXx8fAzt2rUre/vtt1sCxsL+/fffNwWAAQMG5C9evDgQACoqKpCTk+NWW3Z5Q/Tv37/gww8/bAkAiYmJ3hkZGU0b2mZ9cWmdtKTqc/LuMM48zwCYBeAagLdFZDCAYQBuU0oVicguGJfYq2wGcD+AYBhn6ALgPaXUM3X0aemY0mqvKwE0NZ1jLXPn11SMatchIr1gnOn/TkT+adp3pNrxXgBKbmAMRE7rd7/7XdHo0aNzunXrFh0SElIaGxtbAAAffPDByccee6xDfHx8m4qKChk9enTObbfdVrxq1arTEyZM6KDX6wPc3NywYsWKU7VllzfEnDlzLv/pT38K0+v1Ud26dSvq2rVrccuWLSstn2l7DE0hzTAVsE8AnFRKDTNtSwQQAqAbgAEAHlVKjRSRCACHAAxXSu0yHRsNYC2MM9xBAFoC+AzGZfNLItIKgE4pdcp0fAGA2NqOqXnXuYjMBuALYAeA1TD+0uEB40x7LYB3ABxQSnUwHW/2fKXUQjPXfgZAuFKqxHTH/qtKqf+Y9gUBOKiUamv6KOG/SqnIBr3ZRGBoSl0qKipQVlYmPj4+Kjk52euuu+7Snzhx4qi3t7ddimpdoSmckZOWHIGxCL9fY5uvUipbRL4CME1EkgCkw7i8fp1SKtm05HxOKXUBwAUReQ7Av0x3wpcDmA7gVLVzUiwdU5NS6icR2QbgsOm4/QDylFJXROQ7ETkK401y/7yBa/8XjDNwA4BmVUXc1N9FEWlm+iVjCIAvbqBdIqqH/Px8t4EDB3YtLy8XpRSWLVt2yl5F3BLOyInsQER8lVIFIuIDYA+AKUqpAw1oLwbA00qp8RaO+wTAM0qp9Pr2RVSFM3Jgy5YtzZ999tl21beFhoaWfvPNNyccOQ7OyIkcb42IRMH4+fV7DSnigPFmPRH5VkTca349rorp7vatLOJEthMXF3ctLi4upbHHURcWciI7UErV/M63Ldp828L+MgDrbd0vETk3fv2MiIhIw1jIiYiINIyFnIiISMNYyImISDOYR/5bLOREREQNUF5e3qj98651IiKyKGdzRmh5VqFN88ibBDcranW/XpN55MuXL/f/8ssv/UpLS92Kiorcfvjhhwxbvjc3gjNyIiJyWs6aRw4ABw4c8P3ggw9+bswiDnBGTkREVrBm5mwPzppHDgADBw68FhQU1ChBKdWxkBMRkVOrK4/cmvOr55HrdDpDbGxs14bmkQOAj4+PoT7XY2tcWiciIqflrHnkzoSFnIiInFb1PPIRI0Z0rp5H/s477wR07do1Kjw8PHrLli0tAGDVqlWnd+/erTPlhEcdOHCgaVxcXF5FRYXo9fqo+fPnt7VFHrkzYfoZERGZxfQz51FX+hln5ERERBrGm92IiIhq4Sx55HVhISciIqqFFvLIubRORESkYSzkREREGsZCTkREpGEs5ERERBrGQk5ERJrhLHnke/bs8ZkwYUJoY48D4F3rREREN+z2228vuv3224saexwACzkREVlh69atoZcuXbJpHnnr1q2LRo0apck88u3bt+uWLl0a9O233x635XtSH1xaJyIip+XMeeTOgjNyIiKyyJqZsz04cx65s2AhJyIip+aseeTOgkvrRETktJhHbhkLOREROS3mkVvGPHIiIjKLeeTOg3nkRERELoo3uxEREdWCeeREREQaxjxyIiIisisWciIiIg1jISciItIwFnIiItIMZ4kxdSYs5ERERBrGu9aJiMiilNS5oYUFGTaNMW3mqy+KiozXZIzpmDFjOlQ96vXixYtNJk2adGnp0qUXbPn+WIszciIiclrOGmO6adOmU2lpaSnbtm073qJFi4qpU6deccT7YQ5n5EREZJE1M2d7cOYY06KiIomLi+u8bNmy03q9vsw+74BlLOREROTUnDXGdPz48R1GjhyZO2rUqPz6XJetcGmdiIiclrPGmL788suBBQUF7i+99FJWQ9tqKBZyIiJyWs4aY7pixYrg9PT0phEREVERERFRr776amBD26wvxpgSEZFZjDF1HowxJSIiclG82Y2IiKgWjDElIiLSMMaYEhERkV2xkBMREWkYCzkREZGGsZATERFpGAs5ERFpBvPIf4uFnIiISMP49TMiIrJoRurp0LTCEpvmkUc08y56PbK9JvPI+/Tp0/Uf//jH6f79+xcDQO/evSNWrVp16pZbbim25XtkDc7IiYjIaTlrHvmECROy33rrrQAASEpK8iorK5PGKOIAZ+RERGQFa2bO9uCseeQTJkzIXbx4cZvS0tKzb775ZsBDDz3UaM+kZyEnIiKn5ox55DqdzjBw4MBr77//fott27a1SkxMbLSnv3FpnYiInJaz5pEDwLRp07Lnzp0b2rNnz8KgoKBKW7RZHyzkRETktJw1jxwABg4cWNSsWbPKiRMnNmrUK/PIiYjILOaR1y0zM7PJ4MGDu544ceKou7u7XftiHjkREZENrVixwv/WW2+NXLBgwTl7F3FLeLMbERFRLerKI3/iiSeuNNa4qmMhJyIiqgXzyImISMsMBoOh1q9gkWOY/g4Mte1nISciotocvXz5sh+LeeMxGAxy+fJlPwBHazuGS+tERGRWRUXFo1lZWW9lZWV1Ayd+jcUA4GhFRcWjtR3Ar58RERFpGH/DIiIi0jAWciIiIg1jISciItIwFnIiIiINYyEnIiLSsP8P0ZDSsqk/hpAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for dat in data:\n",
    "    wav_deets = FWHM(np.array(dat[1]['Wavelength']), np.array(dat[1]['Transmission']))\n",
    "    depth = average_depths['5s'][average_depths['band'] == dat[0]]\n",
    "    #print(depth)\n",
    "    plt.plot([wav_deets[0],wav_deets[1]], [depth,depth], label=dat[0])\n",
    "    \n",
    "plt.xlabel('Wavelength ($\\AA$)')\n",
    "plt.ylabel('Depth')\n",
    "plt.xscale('log')\n",
    "plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
    "plt.title('Depths on {}'.format(FIELD))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### IV.c - Depth vs coverage comparison\n",
    "\n",
    "How best to do this? Colour/intensity plot over area? Percentage coverage vs mean depth?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Depths (5 $\\\\sigma$) vs coverage on SSDF')"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfIAAAFhCAYAAACLYyg9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XlclWX6P/DPBaiIHgmBWEUWZRfXcC2XzLJy0pzM0syyzMpMy75ZzlTTTDP2S9NxmkrbTNOKLLXNslxQx9TcQEBQQVQEFBTxKLjB/fvjPMcOyCbnHM7i5/168YrzrNd5OHmd+76f575EKQUiIiJyTC62DoCIiIgaj4mciIjIgTGRExEROTAmciIiIgfGRE5EROTAmMiJiIgcGBM5ERGRA2MiJyIicmBM5FQrEckVkcEWOta/RGSqJY5Vw7G3i0icNY5NRGTvmMgdiJZYy0VELyKnRWSLiEwSEYv8HS2ZuKsd1xfAOAALTJZtEJHzInJW+8ky4xSzAbxubpxkHSLST/uslorIKRH5n4jcVN86bX29n3mTbc6a/ATa4r0S2QITueMZppTSAWgPYBaAFwF8ZNuQ6jUewI9KqfJqyycrpVprP1FmHP9bAANFJMCMY9g1EXGzdQyNISJtAHwP4D8A2gIIAvA3ABfqWlftMA35zA8z+Sy1VkrlW+s9EdkbJnIHpZQqVUp9C+B+AA+LSDwAiEigiHwtIkUickhEphj30VouL4lIhoiUiMgnIuKurVsCIATAd1qL5v+03bqISKrWYvrSuL22z4sickxrLWWJyK21hDsUQHJj36uINBORN7T4L4mI0n5StGtxHsBOAENq2HeGiCyvtuzfIjL/Gt8DRKSdiHyjXduTIvKOtjxG62E4LSLpIvKnhpxbe13f3+tFEUkFcE5E3LRjZmvxZojICJPtu4nIbm3dV9rf6x8NOVcN77XG92QS1/TaPhfVRAKAUupzpVSFUqpcKbVGKZVaz7qr1PaZJ7ruKaX44yA/AHIBDK5h+REAT8LwxWwngFcANAcQDiAHwO0m+6cBaAdDC+h/AP5R2/G119sBBGrb7wMwSVsXBeAogEDtdSiAiFriLgJwU7VlG7TlxVocA+p4328C2KrF3QrArwC+ARBuss18AG/XsG97AGUA2mivXQEUAOh1je/BFUAKgLlaDO4A+gFoBuAggJe1az4IgF47dq3n1l435O+1R3vfLbVl92l/DxcYEto5AAHa/ocBPKvFdC+Ai8a/b33nqvZea31P9X0uajhWGwAnAXwKwxc6r4asa+hnvq5t+MOf6+WHLXLnkA/DP6g3AfBVSr2ulLqolMoB8AGA0SbbvqOUOqqUOgXgDQAP1HPs+UqpfG377wB00ZZXAGgBIFZEmimlcpVS2bUc4wYYEoGpF2FIJkEAFsLQExBRfUcR0QGYAuAhLe5zAL4G0FZ7f0Z67TxVKKUOA9gFYLi2aBCAMqXU1mt8D4kwJK4XlFLnlFLnlVKbYfhC0BrALO2ar4Ohu/iBes4NNOzvNV973+Xa+/lK+3tUKqW+BHBAi60XADdt+0tKqW9gSLa4hnMZ1fqeqsVV0+eiCqXUGRi+8CjtfEUi8q2I+NW1rqZjVWP8zBut1HoPTovIygbsT+Q0mMidQxCAUzC0AANN/kE7DUOryvQfxqMmvx+GITnVpdDk9zIY/oGHUuoggKkAXgNwQkS+qOMGoxIAOtMFSqltSim9UuqCUupTGFrld9aw7y0AcpRSB0yWeVWLC9rxT9dy/mX4Iwk9qL2+1vfQDsBhpdTlassDARxVSlWaLDsMw9+k1nNrrvXvBREZJyJ7TLaPB+CjxXFMKaVq2bch52roewJq+VzURCm1Tyk1XikVrMUbCGBefevqYfzMGw1XSt2g/QyvbSciZ8RE7uDEcIdvEIDNMPzDfcjkH7QblFI6pZRpgmxn8nsIDC0bo2sqTq+UWqaU6gdDklAwdIHX5Mp4aF2HAyA1LPeF4YsAAEBEBMAIGFqIpmJg6PquyVcABohIsLbvlWR6De/hKICQGm46ywfQTqo+ORAC4Fh950bD/l5X/iYi0h6GlutkAN5KqRtgGCoRGLrsg7TrY2T6t27IuRr6nhpNKZUJYBEMSbvB60xV+8wTXfeYyB2UiLQRkbsBfAHgM6XUXhi6Us9oN0i1FBFXEYkXk8d5ADwtIsEi0haGFtmXJuuOw9Dd3ZDzR4nIIBFpAeA8gHIYuqpr8iOA/ib73iAit4uIu3YD1xgYWt4/17BvGoBuItJFRFoC+BcMye1Lk+O1ANAdwC81nVwpVQTDmPwnMCSzfY14D9thSJazRKSVFntfANtgGKf+PzHclDcAwDAY/i61ntvkmPX9vUy10t57kRb/I/gj6f2mxT5Zu6b3wNDl3phz1fmeroWIRIvI89oXGYhIOxh6KLbWta6WY9X0mSe67jGRO57vREQPQwtrJoC3ATwCAEqpChj+we0C4BAMN5J9CMDTZP9lANbAcKNTDoB/mKz7F4C/aF2v0+uJowUMjwIVw9DNeiMMXwxqshjAnVoiBgw3U/0Df9zs9gwMXaNXPUuulNoBw1j+j1q8/gDuVEpdMtnsTwA2qLofOVoGYDCqtogb/B5Mrm0HGG60ygNwv1Lqonb+odpx3gUwTmtd1nXuhv69TLfPADAHhqR9HEAnGIYkoMVxL4AJMAwxjIWh1+LCtZ6rge+pofQAegLYJiLnYEjSaQCer2edqVo/80QESNUhNXJmIpIL4DGl1K82OPc/AZxQSjVk/PNaj70NwASlVJqlj+3ItOvyvlLqE1vHQkTW45CTTJDjUUrV1lq3xLF7WuvYjkRE+gPIgqEVPQZAAoCfbBoUEVkdEzmR84gCkATDHeTZAP6slCqwbUhEZG3sWiciInJgvNmNiIjIgTGRExEROTCHGCP38fFRoaGhtg6DiMih7Ny5s1gp5WvG/je6ubl9CMN8BWz42UYlgLTLly8/1r179xM1beAQiTw0NBQ7duywdRhERA5FRA6bs7+bm9uH/v7+Mb6+viUuLi68ocoGKisrpaioKLawsPBDGOZ3uAq/YRERUW3ifX19zzCJ246Li4vy9fUtRR1TFzORExFRbVyYxG1P+xvUmq+ZyImIiBwYEzkREZlt9s9Zfr/uO16lXPGv+47rZv+c1ZD68jZx//33t9+5c6d7bevnz5/vnZub26wxx546dWrgypUrdfVvaT4mciIiMluXkBvKnkvaE25M5r/uO657LmlPeJeQG8psHVttvvzyy8Pdu3c/X9v6zz77zOfIkSONSuTz5s3LHz58uL7x0TUcEzkRkZ3Zvmo5jqSlVll2JC0V21ctt1FE9Rsc46d/e1SXnOeS9oT/7bv0wOeS9oS/PapLzuAYP7OT2eDBgyPi4uJiOnToEDd79mwfAPDw8Oj6+OOPB8fGxsb07t07Mj8/v8ansHbt2uXeqVOnGOPrrKys5pGRkbEAkJiYGLVx40aPy5cvY+TIkaEdO3aMi4yMjP3b3/524yeffOKVlpbmMW7cuPDo6OjYs2fPyvTp0wPi4+NjOnbsGPfAAw+0r6ysrDXmkSNHhn7yySde5r73hmAiJyKyM/4Rkfh+3qwryfxIWiq+nzcL/hGRNo6sboNj/PQjuwUXffK/3ICR3YKLLJHEAWDp0qW56enp+/bs2ZOxYMECv8LCQtfy8nKXbt26lWVkZOzr27evfsaMGYE17dutW7fzly5dkoyMjOYAsHjx4rbDhw8vMd3mt99+8ygoKGh24MCB9P3792c8/fTTJx955JGS+Pj4ssWLF+dkZmZmtG7dWr3wwgsn0tLS9h04cCC9vLzc5Ysvvqix5HBTYyInIrIzIfEJuHvqDHw/bxb+l/QZvp83C3dPnYGQ+ARbh1anX/cd1329K8/3kb6hBV/vyvOtPmbeWG+++aZfVFRUbPfu3WMKCwubpaenu7u4uOCxxx47BQCPPvroye3bt7eubf/hw4ef+uyzz9oCwIoVK7weeuihU6bro6OjLxw9erTFww8/3G758uVtvLy8Kmo6zurVq3UJCQnRkZGRsVu2bNGlpaW1tMT7MxcTORGRHQqJT0DnIXdi69dfoPOQOx0iiRu7018dFpdv7GY3N5l///33uuTkZN2OHTsys7KyMmJiYsrLy8uvyl0iUusxHnrooZKVK1d6paamthARdOrU6YLpel9f34q0tLSMgQMH6t99990bR48eHVr9GGVlZfL888+3/+abb7L379+fMXbs2OLz58/bRQ61iyCIiOgP+uSjOPLrLqSs+RG9Ro5GypofceTXXdAnH7V1aLXac+S0h+mYuHHMfM+R0x7mHPf06dOunp6eFTqdrnL37t3uKSkprQCgsrISxjHoRYsWeScmJtbajR8XF3fBxcUFr7zySuCIESNOVV9fUFDgVlFRgfHjx5/+xz/+cWzv3r0eANC6deuK0tJSVwAoKytzAQB/f//LpaWlLt99912TjH83hENM0UpEdD0puVSES2uKcffo6QgZ3A3t2sbi0ppilAwBdGhn6/BqNP32qOPVlw2O8dObO04+cuTI0oULF/pGRkbGRkREnO/cufM5AGjZsmVlenp6y7i4OH+dTlfxzTff5NR1nHvvvffU3//+9+A333zzWPV1ubm5zSZMmBBaWVkpAPD666/nAcC4ceOKn3nmmfYvvPBC5Y4dO/aNGTOmKDY2Ni44OPiiMY66iEiTTKbjEPXIe/TooTjXOhFdL7avWg7/VuFw++0CWvUMwLltBbjcuwUKz+Ug8Z4/N/g4IrJTKdWjsXGkpKTkdu7cubix+1uTh4dH17Kyst22jqM2gwYN6jBt2rTjw4YNs8gNfykpKT6dO3cOrWkdW+RERHbGmKxLK3OhX3cUukHt4Dk4FCHoZuPIqCHuu+++0PLycpchQ4acbYrzMZETEdmZ7auWo23LQHhsE+gGtcO5bQUoQj5OledfU4vcmdXUGn/ooYdCfv/99yp3rz/55JPHn3322ZPWisMW56yOiZyIyM60bRmInxbNwx3jpyJwSCiKkH/lNdVuyZIlR66Hc1bHRE5EZGf8WrTDHeOnYk3SOzh+Ogcpa37EHeOnwq+Ffd7oRrbFRE5EZGd0/dtBh3Y4fjoHW7/+Ar1GjkaHIX1sHRbZKas9Ry4i7URkvYjsE5F0EXm22vrpIqJExMdaMRAROaojaalVnyOvNvc6kZE1W+SXATyvlNolIjoAO0XkF6VUhoi0A3AbAJuPLRAR2Rvj3OrGaVnbxSY4zDSt1PSs1iJXShUopXZpv+sB7AMQpK2eC+D/ANj/Q+xERE2sMHt/laRtnHu9MHu/jSOrw9q/+yFrddXpWLNW67D279dlPXJjZbXGR9dwTTJFq4iEAugKYJuI/AnAMaVUSlOcm4jI0STe8+erWt4h8Qn2/ehZcI8yrJgUfiWZZ63WYcWkcAT3uC7rkTclqydyEWkN4GsAU2Hobp8J4JUG7DdRRHaIyI6ioiIrR0lEZD/0yUdxPvt0lWXns0/b9VzriBqqx4j3c7BiUjhWzwjEiknhGPF+DqKGXpf1yAHg888/9+rUqVNMaGho/E8//VRrdTZzWTWRi0gzGJL4UqXUNwAiAIQBSBGRXADBAHaJiH/1fZVSC5VSPZRSPXx9fa0ZJhGRXckqyUXRkvQryfx89mkULUlHVkmubQOrT9RQPTo/UIRt7wWg8wNFlkjigOPWI798+bLs3bt335tvvnn09ddfrzE+S7DmXesC4CMA+5RSbwOAUmqvUupGpVSoUioUQB6AbkqpQmvFQUTkaLw7B2Nts70oWpKO0jWGpL622V54dw62dWh1y1qtQ8rnvuj5ZAFSPve9asy8kRy1Hvl9991XAgB9+vQ5l5eX1/xa33dDWbNF3hfAQwAGicge7edOK56PiMgphIWFof+DtyMNR6BfdxRpOIL+D96OsLAwW4dWO+OY+Ij3czB0Vv6VbnYzk7kj1yN3d3dXAODm5oaKioraAzSTNe9a36yUEqVUglKqi/bzY7VtQpVSdllZh4jIlgIqvRBTEYRdrocQUxGEgEq7KX9ds7wdHlXGxI1j5nk7WI/cyjizGxGRnTGOiW9wz0BI70hs+C0Dty5xhe9DcXCPuMHW4dXs1r9eVY8cUUP15o6TO3I98qbCeuRERHbm8IpUrM/YcqU7/dChQ0he9jMGxvZB+xENnxCG9cidB+uRExE5kKO+Z6qMiYeFhQEP3o6jx46hvY1jI/vDRE5EZGf69et31bKwsDD7vtmtibEe+R+YyImIyClcr/XIm2SKViIiIrIOJnIiIiIHxkRORETkwJjIiYiIHBgTORERmW3+rvl+G45uqDId64ajG3Tzd81nPXIrYyInIrIzmzdvxqFDh6osO3ToEDZv3myjiOqX4JtQNnPzzHBjMt9wdINu5uaZ4Qm+CaxHbmVM5EREdsbNbRe+TFp2JZkfOnQIXyYtg5vbLhtHVrsB7Qbo3+j3Rs7MzTPDZ22fFThz88zwN/q9kTOg3YDrth45AFRUVODee+8NnTJlSmBN5zH32gBM5EREdicy6ibERG/Cl0nLsG7dOnyZtAwx0ZsQGXWTrUOr04B2A/TDIoYVLd23NGBYxLAiSyRxwHHrkV+6dEmGDx8e1rFjx/Pz58/Pr+k8lrg+TORERHamrVdv3HzLq/D3y8DGjRvh75eBm295FW29ets6tDptOLpB9132d75jYsYUfJf9nW/1MfPGctR65E899VT72NjY8jfffLPwWs5zrZjIiYjsUOlpfxQURKJdSCoKCiJRetrf1iHVyTgm/ka/N3JmJM7IN3azm5vMHbkeeY8ePc5u2rSpTVlZmTT0PI3BRE5EZGeMY+Kxsb9hwICbERv7W5Uxc3uUWpTqYTombhwzTy1KvW7rkT/xxBPFQ4YMKb377rsjLl26VOt5zMW51omI7MzBg78hJnrTle50L69eAP6Ggwftt3DKlG5TrqpHPqDdAL254+SOXo/8tddeOz5t2jTXe++9N+zll18urOk85mI9ciIiO3P48ALo2iRUGRM/VfIb9GdS0b79Ew0+DuuROw/WIyciciA1Jeu2Xr3t/mY3sg0mciIicjisR/4HJnIiInIKrEdOREREDoeJnIiIyIExkRMRETkwJnIiIiIHxkRORERmOzFvnp9+/foq07Hq16/XnZg3z+b1yKdOnRq4cuVKi8z7XpOuXbtGG3+/+eabO+p0ui4DBw7sYK3zVcdETkREZmvZuXNZ/oszwo3JXL9+vS7/xRnhLTt3tnk98nnz5uUPHz78qhnmLl++bJHj7969O9P4+/Tp0wsXLFjQpHPpMpETEZHZdAMH6gPfnJWT/+KM8MJ//jMw/8UZ4YFvzsrRDRxo03rkADBy5MhQ47zsQUFBnaZPnx7QvXv3qI8//thrzpw5PvHx8TFRUVGxt99+e4Rer3cBgKNHj7rddtttEVFRUbFRUVGxv/zyS6vaju/h4dHV+Ps999yjb9OmTf2Fyi2IiZyIyM4cPrwAp0p+q7LsVMlvOHx4gY0iahjdwIF6z+H3FJUsXhLgOfyeIkskccC8euQ1cXd3r9y5c2fWxIkTS8aMGVOSlpa2LysrKyMqKqp8/vz5PgAwadKkkJtvvlmflZWVkZ6entGtW7fzlngv1sBETkRkZ3RtEpCWNuVKMj9V8hvS0qZA1ybBxpHVTb9+va505Spfr3EPFZSuXOVbfcy8scytR17duHHjSoy/79y5s2X37t2jIiMjY7/++mvv9PR0dwDYsmWL7oUXXigCADc3N3h7e1ukdrg1cGY3IiI709arN+Lj5yMtbQqCgh7EsWPLEB8/367nWjeOiRu701v17q23RPe6aT1ynU5XmZiYGHWt9cir0+l0V7q+J06cGLZ8+fKDvXv3Lp8/f753cnKy1W6Ksxa2yImI7FBbr94ICnoQubnvICjoQbtO4gBQnpLiYZq0jWPm5SkpNq9HXpeysjKXkJCQSxcuXJAvvviirXF537599W+99ZYvYLgp7tSpU3abL9kiJyKyQ6dKfsOxY8sQGjoZx44tg5dXL7tO5jdOnXpVPXLdwIF6c8fJLVWPvDYzZszIT0xMjAkKCroYExNTdvbsWVcAeO+9946MHz++fWRkpI+Liwveeeedw4MHD663Bnn37t2jcnJy3MvLy139/PwS3n333dyRI0eeaUxsDcV65EREdsY4Jm7sTq/+uqFYj9x51FWP3G67CoiIrlf6M6lVkrZxzFx/JtXGkZE9Ytc6EZGdad/+iauWtfXqbddd602tqeuRFxYWug4YMCCq+vINGzZk+fv72/SOdiZyIiJyCtasDe7v71+RmZmZYa3jm4Nd60RERA6MiZyIiMiBWS2Ri0g7EVkvIvtEJF1EntWWvyUimSKSKiIrROQGa8VARETk7KzZIr8M4HmlVAyAXgCeFpFYAL8AiFdKJQDYD+AlK8ZARETk1KyWyJVSBUqpXdrvegD7AAQppdYopYy147YCCLZWDEREjuidw8exuaTqPCqbS/R45/BVc67Yja2rsv0OpRZXmd70UGqxbuuq7OumHvmWLVtadunSJbpDhw5xkZGRsR988IGXtc5pqknGyEUkFEBXANuqrXoUwOqmiIGIyFF0aeOBiem5V5L55hI9Jqbnoksbs2Y7tSq/MM+ytYsywo3J/FBqsW7tooxwvzDP66YeeevWrSuXLFly6ODBg+lr1qw58PLLL7crLi52tchJ6mD1RC4irQF8DWCqUuqMyfKZMHS/L61lv4kiskNEdhQVFVk7TCIiu9HPS4eFcaGYmJ6LN3MKMDE9FwvjQtHPy37reYQl+OhvHR+bs3ZRRvimpP2BaxdlhN86PjYnLMHnuqlHnpCQcKFTp04XACA0NPRS27ZtLxcUFFj9MW+rJnIRaQZDEl+qlPrGZPnDAO4GMEbVMkesUmqhUqqHUqqHr6+vNcMkIrI7/bx0eDjQB3MPH8fDgT52ncSNwhJ89FG9/ItS1+UFRPXyL7JEEgccsx75+vXrPS5duiSxsbEXrvX9Xitr3rUuAD4CsE8p9bbJ8jsAvAjgT0opm3e5EBHZo80lenyaX4xp7f3waX7xVWPm9uhQarEua2uhb8Kg4IKsrYW+1cfMG8vR6pEfPny42SOPPBL+wQcf5Lq6Wr1n3aot8r4AHgIwSET2aD93AngHgA7AL9qy960YAxGRwzGOiS+MC8WL4QFXutntOZkbx8RvHR+bc/OoyHxjN7u5ydy0HnlWVlZGTExMuaXrkb/zzjtH9u/fn/Hiiy/mX7hwway8eOrUKZehQ4d2eOWVV47deuut9VZLswRr3rW+WSklSqkEpVQX7edHpVQHpVQ7k2WTrBUDEZEj2nOmrMqYuHHMfM8Z++3EPH6o1MN0TNw4Zn78UOl1U4/8/Pnzctddd3UYPXr0yUcffbSkvu0thXOtExHZmcntr35iq5+Xzq7HyXvdE3HVs3FhCT56c8fJHake+ccff+z1+++/ty4pKXFbtmyZj7bsUJ8+fcobE1tDsR45EZGTYj1y58F65ERERE6KXetERORwWI/8D0zkRETkFFiPnIiIiBwOEzkREZEDYyInIiJyYEzkREREDoyJnIjIzjhiPfLNXyz2y965vcqMNdk7t+s2f7HY5vXITWVlZTXv2LFjnK3jsCQmciIiO+OI9cgDOkaXrf7vnHBjMs/euV23+r9zwgM6RtvvvLJOgomciMjOOGI98ojuifqhTz+fs/q/c8LXL1oYuPq/c8KHPv18TkT3xs2BbsrceuSbNm3yiIqKiu3SpUv022+/faNx+eXLl/HEE08Ex8fHx0RGRsa+9dZbPsZ1f/nLX/wiIyNjo6KiYp966qkgAKitdvnIkSNDx4wZE9KzZ8/I4ODgTj/88EPr++67LzQ8PDxu5MiRoXW9t7lz5/qEhobGJyYmRo0ePbr9uHHjQq71+jCRExHZIUesRx7RPVEfd8utRbtWfxsQd8utRZZI4oD59cgnTJgQ+vbbbx/Zs2dPpunyefPm+Xh6elakpaXtS0lJ2ffpp5/6ZmZmNk9KSmrzww8/eO3cuTMzKysr49VXXy0EgNpqlwNAaWmp22+//bZ/1qxZR++///6OL7zwwvEDBw6kZ2ZmttyyZUvLmuLKzc1tNnv27IBt27bt27Rp0/4DBw64N+b6MJETEdkhR6xHnr1zuy5941rfbkP/VJC+ca1v9THzxjKnHvnJkydd9Xq961133XXWuK1x3a+//tomKSnJOzo6OrZr164xJSUlbhkZGe6//PJLm7FjxxYby536+flVALXXLgeAu+6667SLiwu6detW5u3tfSkxMbHc1dUVkZGR5dnZ2S1qim3Tpk2tevbsqffz86to0aKFGjFiRKMqpnFmNyIiO2Naj7yflw59vVrbffe6cUzc2J0e0qmL3hLd66b1yHU6XWViYmLUtdQjV0rVtU7mzJlzZOTIkWdMl//4449tatpn4sSJYcuXLz/Yu3fv8vnz53snJydf+WO4u7srAHB1dUXz5s2vVCNzcXHB5cuXawzAUkXL2CInIrIzjliPvOBApodp0jaOmRccyLRpPXIfH5+K1q1bV/z888+ttW2v1By/7bbbSt977z3fCxcuCACkpqa2OHPmjMsdd9xxZsmSJT7GMfDjx4+7ArXXLm+sm2+++dy2bdt0RUVFrpcuXcKqVau8GnMctsiJiKzk/eRsJAR7ok/ElaFUbMkuRmpeKSb1j6h1P0esR95v9Lirno2L6J6oN3ec3BL1yD/66KPcxx57LLRly5aVgwYNutL6njZtWnFubm6LTp06xSilpG3btpd+/PHH7D//+c9ndu3a5dGlS5eYZs2aqcGDB5e+8847x2qrXd5YYWFhl6ZNm1Zw0003xdx4442XIiMjyz09Pa+5AAvrkRMRWcmW7GJMXrYb7zzYFX0ifK56bW2sR27/SktLXTw9PSsvXbqE22+/vcP48eOLx40bd7r6dnXVI2eLnIjISvpE+OCdB7ti8rLdGNszBJ9tO9JkSZwcwwsvvBC4cePGNhcuXJD+/fufGTt27FVJvD5M5EREVtQnwgdje4Zg/rqDmDKoA5O4hTR1PXJzJSQkRF+8eLHKfWmLFy8+tHC0IyqNAAAgAElEQVThwjxzj92gRC4iLQCMBBBquo9S6nVzAyAicmZbsovx2bYjmDKoAz7bdgS9IryZzK3EmvXIzZWamppZ/1aN09AW+SoApQB2ArhgrWCIiJxJ9THxXhHeTTpGTteHhibyYKXUHVaNhIjIyaTmlVZJ2sYx89S8UiZyspiGJvItItJJKbXXqtEQETmRmh4x6xPhwyROFlVnIheRvQCUtt0jIpIDQ9e6AFBKqQTrh0hERES1qa9FfneTREFERA6t9Odcv+YhurKWMd5XJoAp33dSd/GI3sPz9lC7KaSelZXV/O677+544MCBdFvHYil1TtGqlDqslDoM4B/G302XNU2IRERk75qH6MpOJe0PL993UgcYkvippP3hzUN09juvrJ25fPlyo/Zr6FzrcaYvRMQVQPdGnZGc1+Z5wKGNVZcd2mhYTkROrWWMt77tqMicU0n7w09/lx14Kml/eNtRkTmmLfTGcuZ65B4eHl2nTp0amJCQEL127doaK7jVp85ELiIviYgeQIKInBERvfb6BAyPpBH9Iagb8NX4P5L5oY2G10HdbBkVETWRljHe+lbdbiw6+7/8gFbdbiyyRBIHnLceOQCUl5e7xMfHl6empmbefvvtZxtzferrWv+XUkoH4C2lVBullE778VZKvdSYE5ITC7sFuG+RIXmve8Pw3/sWGZYTkdMr33dSd27XCd/WfQMLzu064WvsZjeXs9YjBwxlT8ePH9+oOuRGDX387GURuRdAPxjuYt+klFppzonJSYXdAvSYAGz8f8At/8ckTnSdMI6JG7vTW3S4QW+J7nVnrkcOAM2bN690czNvtvSGjpH/F8AkAHsBpAGYJCL/NevM5JwObQR2fGRI4js+unrMnIic0sUjeg/TpG0cM794RM965FbW0K8B/QHEK63mqYh8CkNSJ/qDcUzc2J0edjO714muEzU9YtYyxltv7ji5M9cjt5QG1SMXkW8ATNMeO4OItAcwSyn1gJXjA8B65A5j8zzDjW2mSfvQRuDYLqDfVNvFRXbl47SPEe8dj8SAxCvLthdsR9rJNDwa/6gNI3M+rEfuPOqqR97QrnVvAPtEZIOIbACQAcBXRL4VkW8tEyY5vH5Tr255h93CJE5VxHvHY3rydGwv2A7AkMSnJ09HvHe8jSMjckwN7Vp/xapRENF1IzEgEbP7z8b05OkYFTUKSVlJmN1/dpUWOlF9nKUeeWJiYrm5x25QIldKJWvd6R2VUr+KSEsAbkopizwjSETXl8SARIyKGoUFqQvwRMITTOJkEddrPfIGda2LyOMAlgNYoC0KBsDHz4ioUbYXbEdSVhKeSHgCSVlJV7rZiejaNXSM/GkAfQGcAQCl1AEAN9a5BxFRDYxj4rP7z8bkrpOvdLM3JJm/n5yNLdlV773akl2M95OzrRUukd1raCK/oJS6aHwhIm4wTAxDRHRN0k6mVRkTN46Zp51Mq3ffhGBPTF62+0oy35JdjMnLdiMh2NOqMRPZs4Ym8mQReRlASxG5DcBXAL6rawcRaSci60Vkn4iki8iz2vK2IvKLiBzQ/utl3lsgIkfyaPyjV42JJwYkNujRsz4RPnjnwa6YvGw33l6ThcnLduOdB7uiT4RPvfsSAYYyph07doyrf8um0b9//w7FxcVmPY/e0EQ+A0ARDJPAPAHgRwB/qWefywCeV0rFAOgF4GkRidWOtVYp1RHAWu01EVGD9InwwdieIZi/7iDG9gxhErcTa9eu9cvKyqoyt3pWVpZu7dq1fraKyREkJycf9PHxqTDnGA1K5EqpShhubntKKfVnpdQHqp6ZZJRSBUqpXdrvegD7AAQBuAfAp9pmnwIY3tjgiej6syW7GJ9tO4Ipgzrgs21HrhozJ9sIDg4uW7FiRbgxmWdlZelWrFgRHhwcbHY9cmcuYxoUFNSpoKDArMnW6ytjKiLymogUA8gEkCUiRSJyTc+Vi0gogK4AtgHwU0oVAIZkj1pumhORiSKyQ0R2FBUVXcvpiMhJGcfE33mwK54bEnWlm53J3PaioqL0I0aMyFmxYkX46tWrA1esWBE+YsSInKioKLMfU3bmMqaWUF+LfCoMd6vfpJUubQugJ4C+IjKtIScQkdYAvgYwVSl1pr7tjZRSC5VSPZRSPXx9fRu6GxE5sdS80ipj4sYx89S8UhtHRoAhmXfu3Llo27ZtAZ07dy6yRBIHnLuMqSXU15wfB+A2pdSVr7tKqRwRGQtgDYC5de0sIs1gSOJLlVLfaIuPi0iAUqpARAIAnGh8+ER0PZnUP+KqZX0ifDhObieysrJ0KSkpvj179ixISUnxDQ8P15ubzJ29jKkl1Ncib2aaxI2UUkUAmtW1oxiuwkcA9iml3jZZ9S2Ah7XfHwawquHhEhGRPTKOiY8YMSJn6NCh+cZu9uo3wF0rljGtX30t8ouNXAcYuuQfArBXRPZoy14GMAtAkohMAHAEwH0NCZSIiOxXXl6eh+mYuHHMPC8vz8OcVjnLmNavzjKmIlIB4FxNqwC4K6XqbJVbCsuYEhFdO5YxdR51lTGts0WulLKLbxtERERUM7OeXSMiIrIFljH9g3Mm8s3zgKBuQNgtfyw7tBE4tgvoN9V2cRERkdWwjKkzCeoGfDXekLwBw3+/Gm9YTkRE5EScs0Uedgtw3yJD8u4xAdjxkeG1aQudiIjICThnixwwJO0eE4CN/8/wXyZxIiJyQs6byA9tNLTEb/k/w3+N3exEREROxDkTuXFM/L5FwKCZf3SzM5kTEVlFdvYcv6LitVVmcSsqXqvLzp5jV2VM7a0euSU4ZyI/tqvqmLhxzPzYLltGRUTktNp4dinLyJgebkzmRcVrdRkZ08PbeHYxu4wp1c05E3m/qVePiYfdwkfPiIisxNfnVn1s7OycjIzp4fv3/z0wI2N6eGzs7Bxfn1vNroDmrPXIly5d6hkdHR0bHR0dGxoaGh8UFNSpMdfHORP5dejjtI+xvWB7lWXbC7bj47SPbRQREV1vfH1u1Qf431t0NG9RQID/vUWWSOKA89YjHzNmTGlmZmZGZmZmRmxsbNnkyZMLG3N9mMidRLx3PKYnT7+SzLcXbMf05OmI9463cWREdL0oKl6rKyj8xrdd8PiCgsJvfKuPmTeWs9cj/8tf/uLn7u5e+dJLLxU15vo453Pk16HEgETM7j8b05OnY1TUKCRlJWF2/9lIDEi0dWhEdB0wjokbu9O92vbRW6J73dnrka9atUq3cuXKtlu3bm30zG9skTuRxIBEjIoahQWpCzAqahSTOBE1mTOlezxMk7ZxzPxM6R4Pc47rzPXI9+/f3/zZZ59tv3z58uzWrVvXXoq0HmyRO5HtBduRlJWEJxKeQFJWEhL9E5nMiahJREQ8f7z6Ml+fW/XmjpM7cz3yBQsWeJeWlroOHz68AwD4+fldTE5OPnitx6mzHrm9YD3y+hnHxI3d6dVfE9H1h/XInUdd9cjZte4k0k6mVUnaxjHztJNpNo6MiIisiV3rTuLR+EevWpYYwK51InJOrEf+ByZyIiJyCqxHTkRERA6HiZyIiMiBMZETERE5MCZyO3Hyww9xbuu2KsvObd2Gkx9+aKOIiIjIETCR2wn3+E44Nm3alWR+bus2HJs2De7xjSqGQ0TUpP6VU+C3pri0ytzqa4pLdf/KKWA9citjIrcTrXr1RNDcuTg2bRqK5s/HsWnTEDR3Llr16mnr0IiI6tW9jUfZM/uOhBuT+ZriUt0z+46Ed2/jwXrkVsZEbkda9eoJrwdGo/jd9+D1wGgmcSJyGEN8PPX/iQnJeWbfkfC/HsgLfGbfkfD/xITkDPHxZD3yWsydO9dnwoQJ7Yyv58yZ4/PYY48FX+v1YSK3I+e2bkPJ51/A56knUfL5F1eNmRMR2bMhPp76Uf5eRR/kFQeM8vcqskQSB5y3HvmECRNOrVmzxtNYtOWzzz7zmThx4jVPXsNEbieMY+JBc+fCd8qUK93sTOZE5CjWFJfqkgpLfB8P9ilIKizxrT5m3ljOWo+8TZs2lX379tV/+eWXnrt373a/dOmSNGamN87sZifOp+2tMiZuHDM/n7aXXexEZPeMY+LG7vSbvXR6S3SvO3s98okTJxa/8cYb/pGRkefHjh3bqAI1bJHbCe/HHrsqYbfq1RPejz1mo4iIiBpu55kyD9OkbRwz33mmjPXI6zBo0KBzBQUFzVesWOE9YcKEU405BlvkRERktpfCA66qRz7Ex1Nv7ji5M9cjNxo+fHhJamqqh6+vb0Vj9mc9ciIiJ8V65I5h4MCBHaZOnXr8nnvuqfVLD+uRExER2Zni4mLX0NDQeHd398q6knh92LVOREQOx1nqkefm5qaZe2wmciIicgqsR05EREQOh4mciIjIgTGRExEROTAmciIiIgdmtUQuIh+LyAkRSTNZ1kVEtorIHhHZISKJ1jo/ERE1ndk/Z/n9uu94lbnVf913XDf75yyL1yPPzc1tdscdd4TXtr64uNh11qxZvo09fteuXaMbu68tWLNFvgjAHdWW/T8Af1NKdQHwivaaiIgcXJeQG8qeS9oTbkzmv+47rnsuaU94l5AbLF6PPDQ09NJPP/1U60xuJ0+edP3oo49urG19fXbv3m21O8ytwWqJXCm1EUD1eWMVgDba754A8q11fiIiajqDY/z0b4/qkvNc0p7wv32XHvhc0p7wt0d1yRkc42fWFK1PPvlkkGnr+rnnngt89dVX/Tp27BgHADt27HDv1KlTTHR0dGxkZGTs3r17Wzz//PPBR48ebREdHR37xBNPBJeWlrr07t07MjY2NiYyMjL2s88+u6Guc3p4eHQ1J+am1tRj5FMBvCUiRwHMBvBSbRuKyESt+31HUVFRkwVIRESNMzjGTz+yW3DRJ//LDRjZLbjI3CQOAGPHjj319ddfXylQsmrVKq9evXqdM77+z3/+4/vUU08dz8zMzEhNTd0XFhZ2cc6cOXnt2rW7kJmZmbFgwYI8Dw+Pyh9++OFgRkbGvuTk5P0vv/xycGVlpbmh2Y2mTuRPApimlGoHYBqAj2rbUCm1UCnVQynVw9e30UMdRETURH7dd1z39a4830f6hhZ8vSvPt/qYeWP07du3/OTJk265ubnNfvvtt5aenp4V4eHhF43re/fufW7OnDkBM2fO9D9w4EDz1q1bX1VApLKyUqZOnRocGRkZO3DgwMgTJ040z8vLc5oJ0Zo6kT8M4Bvt968A8GY3IiInYBwTf3tUl5xXh8XlG7vZLZHMhw0bVvLZZ595LV26tO3IkSOrDNlOmjTp1KpVqw62bNmycujQoZHffvvtVedbsGBB25MnT7rt3bt3X2ZmZoa3t/elmmqaO6qm/kaSD6A/gA0ABgE40MTnJyIiK9hz5LSH6Zi4ccx8z5HTHuZ2sT/00EOnHn/88dCSkhK35OTkrPPnz4txXUZGRvOYmJgLcXFxJ3Jyclrs2bOnZWJiYtm5c+euJOrS0lJXHx+fSy1atFDfffedLj8/v7k58dgbqyVyEfkcwAAAPiKSB+BVAI8D+LeIuAE4D2Citc5PRERNZ/rtUVfVIx8c46e3xDh5jx49zp87d87Fz8/vYvv27S9lZWVdScRLlixp+9VXX3m7ubkpX1/fS//617/y/fz8Krp37362Y8eOcYMGDSp97bXXCocOHdohPj4+Ji4uriwsLOx8XecTkbpW2x3WIyciclLOXI/cWgoLC127desWm5+fv9fWsZhiPXIiIqJ65ObmNuvVq1fM008/fVXvgj1zmrv2iIiIGqqwsNB1wIABUdWXb926dZ+/v3+FLWJqLCZyIiK67vj7+1dkZmZm2DoOS2DXOhERkQNjIiciInJgTOREREQOjImciIjIgTGRExGR+db+3Q9Zq6tOj5q1Woe1f3eoeuRZWVnNjZXVHAUTORERmS+4RxlWTAq/ksyzVuuwYlI4gns4XD1yR8NETkRE5osaqseI93OwYlI4Vs8IxIpJ4Rjxfg6ihjpcPfKKigqMHj26fYcOHeL69u3b8ezZs3Y9ZysTORERWUbUUD06P1CEbe8FoPMDReYmccA29ciPHDniPmXKlBMHDx5M9/T0rFi8eLGXue/DmpjI7cjJDz/Eua3bqiw7t3UbTn74oY0iIiK6BlmrdUj53Bc9nyxAyue+V42ZN4It6pEHBQVd6NOnTzkAdO3atSw3N7eFue/DmpjI7Yh7fCccmzbtSjI/t3Ubjk2bBvf4TjaOjIioHsYx8RHv52DorPwr3ewWSOZNXY+8efPmV74MuLq6qsuXL9t11zqnaLUjrXr1RNDcuTg2bRq8HhiNks+/QNDcuWjVq6etQyMiqlveDo8qY+LGMfO8HR7mdrGzHnndmMjtTKtePeH1wGgUv/sefJ56kkmciBzDrX+9umJY1FC9JcbJm7oeuaNhPXI7Y+xOZ4uciMzFeuTOg/XIHYQxiQfNnQvfKVOudLNXvwGOiIjIiF3rZtr182HcGNoGwVF/PJ2Ql1WCE7ln0O329td0rPNpe6u0wI1j5ufT9rJVTkRkQbXVI9+wYUMW65FfZ24MbYOfP0jD7Y/HIzjKC3lZJVdeXyvvxx67almrXj2ZxImILMyZ6pEzkZspOMoLtz8ej58/SEP8LUFI23jsSlInIiKyNo6RW0BwlBfibwnCjh9zEX9LEJM4ERE1GSZyC8jLKkHaxmPocWco0jYeQ15Wia1DIiKi6wQTuZlMx8R7/in8Sjc7kzkRETUFJnIzncg9U2VM3DhmfiL3jI0jIyJqOvN3zffbcHRDlelRNxzdoJu/az7rkVsZE7mZut3e/qox8eAor2t+9IyIyJEl+CaUzdw8M9yYzDcc3aCbuXlmeIJvAuuRWxkTORERmW1AuwH6N/q9kTNz88zwWdtnBc7cPDP8jX5v5AxoN8Dh6pEbafO4xyYnJ3vUdB5z3pclMZETEZFFDGg3QD8sYljR0n1LA4ZFDCsyN4kDtqlHDgApKSktRo4c2eGjjz461L9//7KazmPue7MUPkdOREQWseHoBt132d/5jokZU/Bd9ne+vQJ66c1N5qb1yAsKCtxqqkc+e/bsgLy8vOajR48u6dSp04XqxzDWI9+6dWtrFxcXGOuRh4SEXK7pnKdOnXIbPnx4h6+++iq7R48e5xt6Hlthi5yIiMxmHBN/o98bOTMSZ+Qbu9mr3wDXGE1dj1yn01UEBARc3LBhQ+trOY+tsEVORERmSy1K9TAdEzeOmacWpXqY2ypv6nrkzZo1Uz/99FP2wIEDO7Zu3bpy0qRJp2o6z5/+9Cezhw4sgYmciIjMNqXblKvqkQ9oN8DsrnXANvXI27RpU/nzzz8fHDBgQGTr1q0r09PT3aufx9z3ZSmsR24hlqyCRkRkCaxH7jxYj7wJGKugGWd0M874dmNoGxtHRkREzoxd6xbCKmhERI6D9cipRqZV0HrcGcokTkRkp5ypHjm71i2IVdCIiKipMZFbCKugERGRLTCRWwiroBERkS1wjNxCanrELDjKi+PkRERkVVZrkYvIxyJyQkTSqi1/RkSyRCRdRP6ftc5PRERN58S8eX769eurTFuqX79ed2LePLPqkXft2jXavMhqt3TpUs+XX37ZHwBWr17dOjY2NsbNza37J5984lAtMGt2rS8CcIfpAhEZCOAeAAlKqTgAs614fiIiaiItO3cuy39xRrgxmevXr9flvzgjvGXnzmbVI9+9e3dm9WWXL9dY6+SajRkzpvSf//xnIQCEh4df/OSTT3KHDRt20iIHb0JWS+RKqY0ATlVb/CSAWUqpC9o2J6x1fiIiajq6gQP1gW/Oysl/cUZ44T//GZj/4ozwwDdn5egGDjRrilYPD4+uAPD999/revbsGTls2LCwqKioOAAYPHhwRFxcXEyHDh3iZs+e7WPcZ/ny5W1iY2NjoqKiYnv37h1Z27Hnz5/vPW7cuBAAiIqKutizZ89yFxfHu3WsqcfIIwHcLCJvADgPYLpS6veaNhSRiQAmAkBISEjTRUhERI2iGzhQ7zn8nqKSxUsCvMY9VGBuEq8uNTW11e7du9Ojo6MvAsDSpUtz/fz8Ks6ePStdu3aNHTt2bEllZaVMnjw5dMOGDZnR0dEXjx8/7mrJGOxRUydyNwBeAHoBuAlAkoiEqxomfFdKLQSwEDDMtd6kURIR0TXTr1+vK125ytdr3EMFpStX+bbq3VtvyWSekJBwzpjEAeDNN9/0++GHH24AgMLCwmbp6enux48fd0tMTNQbt/Pz83OoWdoao6kTeR6Ab7TEvV1EKgH4ACiqa6edO3cWi8hhK8fmA8ARigMwTstinJbFOC3PnFibrGKTcUzc2J3eqndvvaW61408PDwqjb9///33uuTkZN2OHTsydTpdZWJiYlR5ebmLUgoiUtdhnE5TJ/KVAAYB2CAikQCaowEfUKWUr7UDE5Ed5lQJaiqM07IYp2UxTstzlFjLU1I8TJO2ccy8PCXFw9Jd7ABw+vRpV09PzwqdTle5e/du95SUlFYAMHDgwHPPP/98+8zMzObGrnVnb5VbLZGLyOcABgDwEZE8AK8C+BjAx9ojaRcBPFxTtzoRETmWG6dOvaoeuW7gQIt2rZsaOXJk6cKFC30jIyNjIyIiznfu3PkcAAQGBl6eP39+7ogRIzpUVlbC29v70pYtWw7Ud7zk5GSPUaNGdThz5ozr2rVrb3jjjTcCDx48mG6N2C3NIeqRNwVH+dbLOC2LcVoW47Q8W8bKeuT2g/XIG2ahrQNoIMZpWYzTshin5TlSrGQDbJETEVGNnKVF/u9//9v7vffeqzLD3E033XR2yZIlR2wV07Wqq0XORE5ERDVylkTuDNi1TkRE5KScNpGLyB1acZaDIjKjhvUhIrJeRHaLSKqI3Kktv01EdorIXu2/g0z22aAdc4/2c6MN4wwVkXKTWN432ae7Fv9BEZkvFnio0ow4x5jEuEdEKkWki7bOFtezvYis1WLcICLBJuseFpED2s/DJsttcT1rjFNEuojIb1rRoVQRud9kn0UicsjkenaxVZzaugqTWL41WR4mItu06/yliDS3VZwiMrDa5/O8iAzX1ln0ekothaRM1ov2+TqoxdnNZF2TfTbJ8Thl17qIuALYD+A2GCah+R3AA0qpDJNtFgLYrZR6T0RiAfyolAoVka4Ajiul8kUkHsDPSqkgbZ8NMEwru8MO4gwF8L1SKr6G424H8CyArQB+BDBfKbXaFnFWO04nAKuUUuHa6w1o+uv5FQzX7VMxfEl7RCn1kIi0BbADQA8ACsBOAN2VUiU2up61xRkJQCmlDohIoBZnjFLqtIgs0vZZ3tjYLBWntu6sUqp1DcdNgmFiqC/E8AU0RSn1nq3iNNmmLYCDAIKVUmVWuJ63ADgLYHEt/9/eCeAZAHcC6Ang30qpnk352ayOXev243rsWk8EcFAplaOUugjgCxiqrplSANpov3sCyAcApdRupVS+tjwdgLuItLC3OGsjIgEA2iilftOe0V8MYLidxPkAgM/NjKUuDYkzFsBa7ff1JutvB/CLUuqUUqoEwC8A7rDh9awxTqXUfqXUAe33fAAnAFhrwiRzrmeNtBbjIADG5PgpbHg9q/kzgNVKKbOqddWmlkJSpu6BIckrpdRWADdon7+m/GySA3LWRB4E4KjJ6zxtmanXAIwVw2Q1P8LwTbi6kTC0Mi+YLPtE62b7qwW6scyNM0wMXdnJInKzyTHz6jlmU8dpdD+uTuRNfT1TYPi7AsAIADoR8a5jX1tdz9rivEJEEmGYHTHbZPEbWrfsXAt8ATU3TncR2SEiW43d1QC8AZxWShnrUNrN9QQwGld/Pi15PetT12ewqT6bjbZ1VbbfodTiKvXID6UW67auynaIeuSvvfaaX0RERFxkZGRs7969I/fv32/2kE9TcdZEXlNCqD6G8ACARUqpYBi6spaIyJXrISJxAN4E8ITJPmOUUp0A3Kz9VOmaa+I4CwCEKKW6AngOwDIRadPAYzZlnIYDiPQEUKaUMh0ftMX1nA6gv4jsBtAfwDEAl+vY11bXs7Y4DQcwtMaWwNBFbJx/+iUA0TAUJGoL4EUbxxmiTWTyIIB5IhLRwGM2dZzG69kJwM8m+1j6etbnWj+D1riWjeYX5lm2dlFGuDGZH0ot1q1dlBHuF+bpEPXIu3fvXrZnz559+/fvzxg+fHjJtGnTguvb3144ayLPA9DO5HUwru7qnQAgCQCUUr8BcIehOAHEcCPMCgDjlFJXWjtKqWPaf/UAlsHQpWeTOJVSF5RSJ7XlO2FolUVqxzT9ANZ0zCaL02T9Va0dW1xPpVS+Uupe7QvQTG1ZaR372uR61hEntC9sPwD4i9YFa9ynQOuWvQDgE9j2ehq7/qGUygGwAUBXGGor3CAibrUds6nj1IwCsEIpdclkH0tfz/rU9Rlsqs9mo4Ul+OhvHR+bs3ZRRvimpP2BaxdlhN86PjYnLMHHIeqRDxs2TK/T6SoBoF+/fmcLCgrYIrex3wF0FMPdsc1hSCLfVtvmCIBbAUBEYmBIPEUicgMM/0i+pJT6n3FjEXETEWOibwbgbgA13n3aRHH6ajf5QETCAXQEkKOUKgCgF5FeWlf1OACrbBWn9toFwH0wjF1CW2aT6ykiPiY9BS/BMP8/YGiJDRERLxHxAjAEhhsdbXI9a4tT234FDGOpX1XbJ0D7r8AwVmqz66ldxxbGbQD0BZChjeWuh2E8GgAehg2vp4mr7t+wwvWsz7cAxolBLwCl2uevKT+bZglL8NFH9fIvSl2XFxDVy7/I3CReXWpqaqu33nrrWHZ2djpgqEeenp6+b8+ePRkLFizwKywsdM3Pz3ebPHly6DfffJOdlZWVsXLlyuz6jlvdggULfAcPHlxa/5b2oamrnzUJpdRlEZkMw/8ArgA+Vkqli8jrAHYopb4F8DyAD0RkGgzdUeOVUkrbrwOAv4rIX2gSw9cAAA/uSURBVLVDDgFwDsDPWtJxBfArgA9sGOctAF4XkcsAKgBMUkoZb6R5EsAiAC0BrNZ+bBKndohbAORpLTOjFrDN9RwA4F8iogBsBPC0tu8pEfk7DEkBAF638fWsMU4YWo63APAWkfHasvFKqT0AloqILwxdrnsATLJhnDEAFoihVLELgFnqj7vIXwTwhYj8A8BuAB/ZME6I4QmQdgCSqx3aotdTai4k1Ux7D+/DcG/JnTDcOV8G4BFtXZN9Ns11KLVYl7W10DdhUHBB1tZC3+DotnpLJvOmqEf+7rvvtk1JSfFYsGBBlqXitjanfPyMiIjMdy2PnxnHxI3d6dVfNzYGDw+PrmVlZbu///573Zw5c/zWr19/EDB0tb/66quBGzZsOGCsR/7KK6/kl5aWuiQlJbVdtWrVofqOPX/+fO8dO3a0Wrx48REAWLlype65554L2bRpU1ZQUJBlBuIt5Hp8/IyIiJrQ8UOlHqZJ2zhmfvxQqYc1zldXPfJt27bpMjMzmwPA8ePHXRtyvP/9738tn3nmmfarVq06aG9JvD5O2bVORERNq9c9EVfVIw9L8LFo17opS9cjf+GFF9qVlZW53nfffRHacS6uW7fuoDVitzR2rRMRUY04s5v9YNc6ERGRk2LXOhEROTVnqEdeF3at03VBRPwBzINhlq4LAHIBTFVK7bdlXET2jF3r9oNd63Rd0ybLWAFgg1IqQikVC+BlAGbNAV3LuRp0hywRkaUwkdP1YCCAS9qkGwAAbQKVzSLyloikiaGm8/0AIIYa2XcatxVDXeqRIuKqbf+7GAppPKGtHyCGWuzLAOzVlq0UQz37dBGZaHKsCSKyXww1sT8QkXe05b4i8rV27N9FpG+TXBkicngcI6frQTwMNZyruxdAFwCdYZgX/ncR2QjDVLL3A/hRm/LzVhhm0JoAw7SZN2nTj/5PRNZox0oEEK+UMk5C8ag2I1dL7bhfwzCb3V8BdAOgB7AOhqpcAPBvAHOVUptFJASGWcpiLHcJiMhZMZHT9awfgM+VUhUAjotIMgxj6KsBzNeS9R0ANiqlykVkCIAEETHOE+4Jwxz3FwFsN0niADBFREZov7fTtvMHkGycXlNEvoKh0A0ADAYQK39Ucm0jIjqtoAwRUa3YtU7Xg3QA3WtYXmP9c6XUeRiqdd0OQ8v8C5Ptn1FKddF+wpRSxhb5uSsHFRkAQ2LurZTqDMN84u61nU/jom1vPHYQkzg5ks1fLPbL3rm9Sj3y7J3bdZu/WGzRe1Gee+65wFdeecXi97c4MiZyuh6sA9BCRB43LhCRmwCUALhfG/v2haEYyXZtky9gKFpxM/6oUf0zgCe1Qi8QkUgRaVXD+TwBlCilykQkGkAvbfl2GGpie4mhjOdIk33WAJhsEl8Xs94xURML6Bhdtvq/c8KNyTx753bd6v/OCQ/oGG1WPXKqHxM5OT2tCtsIALeJSLaIpAN4DYYa6KkwjFOvA/B/SqlCbbc1MCT2X5VSxmpLHwLIALBLRNIALEDNw1M/AXATkVQAfwewVYvjGIB/AtgGQ7W3DADGUolTAPTQbqLLgJmVtoiaWkT3RP3Qp5/PWf3fOeHrFy0MXP3fOeFDn34+J6J7otk9Sy+++KJ/aGhofJ8+fSIPHDjQAgDS09Nb3HzzzR3j4uJiunfvHrV79253ADh69KjbbbfdFhEVFRUbFRUV+8svv7QCaq9d7uHh0fXJJ58MiouLi+nTp0/k+vXrPRITE6OCg4M7LV261LO2mPR6vcudd94ZHhkZGXvXXXeFJyQkRG/cuNEq88rXh8+REzUhEWmtlDqrtchXwFByc4Wt4yKqSWOeI1+/aGHgrtXfBnQb+qeCgeMn5psbw6ZNmzwmTJgQunPnzsxLly6hS5cusePHjy/65ZdfPBcuXHi4U6dOF9atW9fq5ZdfDtq6dev+u+66K7xnz55nX3nllROXL19GaWmpq7e3d8Xx48dd/fz8Ks6ePStdu3aN3bRpU6a/v3+FiHT/8ssvD4waNerMbbfdFlFWVuaybt26g7t27XJ/5JFHwjIzMzNqiuuVV17xO3jwoPuyZcsO//777+69e/eOW7du3b5bbrnFKj0QdT1HzpvdiJrWayIyGIYx8zUAVto4HiKLyd65XZe+ca1vt6F/KkjfuNY3pFMXvbkt8vXr17e+8847T+t0ukoAGDJkyOnz58+77N69u7WxwAkAXLx4UQBgy5YtuuXLlx8CADc3N3h7e1cANdcu9/f3P9esWTP15z//+QwAxMXFlbdo0aKyRYsWKjExsfzYsWPNa4try5YtrZ999tkTAHDTTTedj4yMtNkQAhM5URNSSk23dQxE1mAcEzd2p4d06qK3VPe6ydMcAIDKykrodLrLtbWWq/v+++91ycnJuh07dmQaa5eXl5e7AICbm5tycTGMMru4uKBFixYKAFxdXVFRUVHrDar21JvNMXIiIjJbwYFMD9OkbRwzLziQada48aBBg87+8MMPN5z9/+3dXUxTaR7H8QfMjlKmQqEdBAHNQVtgSJhAIr7ViCRMJtHEThPcmLipiUZM9EqNGXW9MdEoMyFxiGjWqDeKYpgmLokZCRdKVuUCIkpby1oWGXapU14WDpYVat2LnSboUGeXc7p9yfdzpzSH88SLn+f/PD2/qamk8fHx5La2tnSNRhPMzc2duXLlik6I/wT7o0ePUoQQYsOGDXJdXZ1BCCECgYAYGxtLDtddrsT69eunbt68qRNCiK6uriV9fX0pSq+5UDyRAwAU2/j7P/yqj7ygfI3i0frGjRv9FotlrKSk5PPly5e/WbNmzZQQQjQ1NfXv3bt3xdmzZ7MDgUCSxWIZW7du3XRjY+OgzWZbYTQa9cnJyaKhoeFluO5yJY4cOeKrqalZaTQai0tKSvwmk2lap9O9VXrdheCwGwBgXpSmhBcIBMTMzEySRqN553A4FldXVxs9Hk/vkiVLIhKqHHYDAEBFsiwnm81m0+zsbNK7d+9EfX39y0iF+G8hyAEACKOlpWXp8ePHc+f+XV5e3pu2tjZPb2+vK1r3NRdBDgBAGFarddJqtf5Xp+OjhVPrAADEMYIcAIA4RpADABDHCHIAAOIYQQ4AUGzix4Gsadfoe33k065R7cSPA/SRRxhBDgBQ7JN8rX+suU8Khfm0a1Q71twnfZKvTfg+8tnZ2aj+foIcAKBYSlGmnFFj7B9r7pP++WdPzlhzn5RRY+xPKcpMyD7y8+fPZ3711VfSli1bVpnNZqPSNSpBkAMAVJFSlCmnln3mm/rLP7JTyz7zqRHiHR0dGrvdnvHs2TNna2vri1DhyZ49e1ZcuHBh0OFwuOrq6ob279+fL4QQtbW1+WazWXa73U6Hw+EsKyv7lxBCXL9+fcDhcLiePHnivHTpUpbX610khBDT09PJlZWVssPhcKWmpr49ceLE8o6Ojr7bt2+/OHXq1PKP3Vt3d/enTU1Nf3v8+HGf0nUqwQthAACqmHaNal93/2z4dEPO8Ovunw2LV6XLSsM8VvvIhRDCbDZPZmVlRaUoZS6CHACgWGhPPDROX7wqXVZrvB6LfeRCCKHRaIILWY/aGK0DABSbGZQ1c0M7tGc+MygnZB95LCHIAQCKpX258tWHT94pRZly2pcrf9VT/r+Y20e+devWgrl95FevXtWbTKbi1atXf97S0pIuhBCNjY2D9+/f1/7SE17c3d2dYrVaJwKBQJLRaCw+duxYjhp95LGEPnIAwLzoI48dH+sj54kcAIA4xmE3AADC+FgfebTu6UMEOQAAYdBHDgAAIoogBwAgjhHkAADEMYIcAIA4RpADABRrb2/Pcrvd7/WRu91ubXt7e0L2kT948EBjs9nyon0fQhDkAAAV5Obm+u12uxQKc7fbrbXb7VJubm5C9pFv2rTJf+3atZ+ifR9CEOQAABWYTCbZYrH02+126e7duzl2u12yWCz9JpMpIfvIW1tbtZWVlauUrk0NBDkAQBUmk0kuLS31dXZ2ZpeWlvrUCPFY7iOPFbwQBgCgCrfbre3p6TFUVFQM9/T0GCRJkpWGeSz3kccKghwAoFhoTzw0TpckSVZrvB6rfeSxgtE6AECxoaEhzdzQDu2ZDw0N0UceYQQ5AECxqqqqVx8+eZtMJrmqqoo+8gijjxwAMC/6yGMHfeQAACQoDrsBABAGfeQAAMQx+sgBAEBEEeQAAMQxghwAgDhGkAMAFPN4vsvyjbS/V2PqG2nXejzfJWSNaSwhyAEAii1N+8LvdB6WQmHuG2nXOp2HpaVpXyRkjWksIcgBAIoZ9FVycfG3/U7nYamv71SO03lYKi7+tt+gr0rIGtMdO3asKCwsLC4sLCzW6XSlhw4dyla6zoUiyAEAqjDoq+TsZV/7fhq6lp297GufGiEeqzWmt27devn8+XPnnTt3XqSnpwf27ds3qnStC8X3yAEAqvCNtGuHvT8Y8nJtw8PeHwy6jPWy0jCP5RpTv9+fZLVaC+rr6weNRuOMknUqQZADABQL7YmHxum6jPWyWuP1WK0x3bVr14pt27aNb9++XfHkQQlG6wAAxSYnnmjmhnZoz3xy4klC1pieOXPGMDU1tej06dNepddSiiAHAChWUHDo1YdP3gZ9lVxQcCgha0wbGhqWud3ulNCBt3PnzhmUXnOhqDEFAMyLGtPYQY0pAAAJisNuAACEQY0pAABxjBpTAAAQUQQ5AABxjCAHACCOEeQAAMQxghwAoNiZ/uGseyMT7/WR3xuZ0J7pH6aPPMIIcgCAYuVLNf6DrkEpFOb3Ria0B12DUvlSDX3kEUaQAwAUq9anyd8X5fcfdA1Kf/zrUM5B16D0fVF+f7U+LSH7yMvLy00PHz5MCf25rKyssLOzMyXc5yOJIAcAqKJanybXLNP5/jQ0kl2zTOdTI8RjtY/cZrONXL58WS+EEE+fPl08MzOTVFFRMa10vQvBC2EAAKq4NzKhbfaOG/bm6oebveMGs04rKw3zWO0jt9ls43V1ddlv3rwZunjxon7nzp1Reyc9QQ4AUCy0Jx4ap5t1Wlmt8Xos9pFrtdqg2WyevHHjRvqdO3cyurq6ovb2N0brAADFuib9mrmhHdoz75r0J2QfuRBC1NbWjhw9ejSvtLT0dVZW1ls1rrkQBDkAQLFvpOxXHz55V+vT5G+k7ITsIxdCCLPZ7E9NTX27e/fuqFa90kcOAJgXfeQfNzAw8LvNmzebPB5P76JFiyL6u+gjBwBARQ0NDZlr164tOnny5N8jHeK/hcNuAACE8bE+8gMHDoxG677mIsgBAAiDPnIAQDwLBoPBsF/Bwv/HL/8GwXA/J8gBAOH0+ny+NMI8eoLBYJLP50sTQvSG+wyjdQDAvAKBwB6v13vZ6/WWCB78oiUohOgNBAJ7wn2Ar58BABDH+B8WAABxjCAHACCOEeQAAMQxghwAgDhGkAMAEMf+DUA5pYvCXzhuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for dat in data:\n",
    "    wav_deets = FWHM(np.array(dat[1]['Wavelength']), np.array(dat[1]['Transmission']))\n",
    "    depth = average_depths['5s'][average_depths['band'] == dat[0]]\n",
    "    #print(depth)\n",
    "    coverage = np.sum(~np.isnan(depths['ferr_{}_mean'.format(dat[0])]))/len(depths)\n",
    "    plt.plot(coverage, depth, 'x', label=dat[0])\n",
    "    \n",
    "plt.xlabel('Coverage')\n",
    "plt.ylabel('Depth')\n",
    "#plt.xscale('log')\n",
    "plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
    "plt.title('Depths (5 $\\sigma$) vs coverage on {}'.format(FIELD))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (herschelhelp_internal)",
   "language": "python",
   "name": "helpint"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}