{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# XMM-13hr 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 12:32:49.540731\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 = 'XMM-13hr'\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_xmm-13hr_20180501.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=\"table140260369485608-992409\" 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>171770621</td></tr>\n",
       "<tr><td>1</td><td>171770622</td></tr>\n",
       "<tr><td>2</td><td>171770623</td></tr>\n",
       "<tr><td>3</td><td>171770786</td></tr>\n",
       "<tr><td>4</td><td>171770787</td></tr>\n",
       "<tr><td>5</td><td>171770790</td></tr>\n",
       "<tr><td>6</td><td>171770791</td></tr>\n",
       "<tr><td>7</td><td>171770792</td></tr>\n",
       "<tr><td>8</td><td>171770793</td></tr>\n",
       "<tr><td>9</td><td>171770794</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(\"$('#table140260369485608-992409').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",
       "    $('#table140260369485608-992409').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=\"table140262023737528-455938\" 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>171770621</td><td>2683915</td></tr>\n",
       "<tr><td>1</td><td>171770622</td><td>2683915</td></tr>\n",
       "<tr><td>2</td><td>171770623</td><td>2683915</td></tr>\n",
       "<tr><td>3</td><td>171770786</td><td>2683918</td></tr>\n",
       "<tr><td>4</td><td>171770787</td><td>2683918</td></tr>\n",
       "<tr><td>5</td><td>171770790</td><td>2683918</td></tr>\n",
       "<tr><td>6</td><td>171770791</td><td>2683918</td></tr>\n",
       "<tr><td>7</td><td>171770792</td><td>2683918</td></tr>\n",
       "<tr><td>8</td><td>171770793</td><td>2683918</td></tr>\n",
       "<tr><td>9</td><td>171770794</td><td>2683918</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(\"$('#table140262023737528-455938').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",
       "    $('#table140262023737528-455938').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 masked=True length=10</i>\n",
       "<table id=\"table140260366492000-201898\" 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_ukidss_j_mean</th><th>f_ap_ukidss_j_p90</th><th>ferr_ukidss_j_mean</th><th>f_ukidss_j_p90</th><th>ferr_ap_90prime_g_mean</th><th>f_ap_90prime_g_p90</th><th>ferr_90prime_g_mean</th><th>f_90prime_g_p90</th><th>ferr_ap_90prime_r_mean</th><th>f_ap_90prime_r_p90</th><th>ferr_90prime_r_mean</th><th>f_90prime_r_p90</th><th>ferr_ap_mosaic_z_mean</th><th>f_ap_mosaic_z_p90</th><th>ferr_mosaic_z_mean</th><th>f_mosaic_z_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>uJy</th><th>uJy</th><th></th><th></th><th>uJy</th><th>uJy</th><th></th><th></th><th>uJy</th><th>uJy</th></tr></thead>\n",
       "<tr><td>0</td><td>171728253</td><td>2683253</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>0.2503924</td><td>2.0666188478469847</td><td>0.37544167</td><td>2.0267515182495117</td><td>0.20270073</td><td>2.926134943962097</td><td>0.2498269</td><td>2.9660435914993286</td><td>1.2030196</td><td>10.519013404846191</td><td>0.60006213</td><td>8.704996156692506</td></tr>\n",
       "<tr><td>1</td><td>171728255</td><td>2683253</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>0.2503924</td><td>2.0666188478469847</td><td>0.37544167</td><td>2.0267515182495117</td><td>0.20270073</td><td>2.926134943962097</td><td>0.2498269</td><td>2.9660435914993286</td><td>1.2030196</td><td>10.519013404846191</td><td>0.60006213</td><td>8.704996156692506</td></tr>\n",
       "<tr><td>2</td><td>171728340</td><td>2683255</td><td>5.9318357</td><td>356.375198364258</td><td>12.612087</td><td>413.83007202148434</td><td>0.25039154</td><td>18.096925163269056</td><td>0.4938725</td><td>29.27996597290041</td><td>0.19070134</td><td>26.524150848388686</td><td>0.2986061</td><td>33.34638252258304</td><td>1.203017</td><td>98.89707946777344</td><td>0.95266587</td><td>111.70701904296901</td></tr>\n",
       "<tr><td>3</td><td>171728341</td><td>2683255</td><td>5.9318357</td><td>356.375198364258</td><td>12.612087</td><td>413.83007202148434</td><td>0.25039154</td><td>18.096925163269056</td><td>0.4938725</td><td>29.27996597290041</td><td>0.19070134</td><td>26.524150848388686</td><td>0.2986061</td><td>33.34638252258304</td><td>1.203017</td><td>98.89707946777344</td><td>0.95266587</td><td>111.70701904296901</td></tr>\n",
       "<tr><td>4</td><td>171728342</td><td>2683255</td><td>5.9318357</td><td>356.375198364258</td><td>12.612087</td><td>413.83007202148434</td><td>0.25039154</td><td>18.096925163269056</td><td>0.4938725</td><td>29.27996597290041</td><td>0.19070134</td><td>26.524150848388686</td><td>0.2986061</td><td>33.34638252258304</td><td>1.203017</td><td>98.89707946777344</td><td>0.95266587</td><td>111.70701904296901</td></tr>\n",
       "<tr><td>5</td><td>171728343</td><td>2683255</td><td>5.9318357</td><td>356.375198364258</td><td>12.612087</td><td>413.83007202148434</td><td>0.25039154</td><td>18.096925163269056</td><td>0.4938725</td><td>29.27996597290041</td><td>0.19070134</td><td>26.524150848388686</td><td>0.2986061</td><td>33.34638252258304</td><td>1.203017</td><td>98.89707946777344</td><td>0.95266587</td><td>111.70701904296901</td></tr>\n",
       "<tr><td>6</td><td>171728344</td><td>2683255</td><td>5.9318357</td><td>356.375198364258</td><td>12.612087</td><td>413.83007202148434</td><td>0.25039154</td><td>18.096925163269056</td><td>0.4938725</td><td>29.27996597290041</td><td>0.19070134</td><td>26.524150848388686</td><td>0.2986061</td><td>33.34638252258304</td><td>1.203017</td><td>98.89707946777344</td><td>0.95266587</td><td>111.70701904296901</td></tr>\n",
       "<tr><td>7</td><td>171728345</td><td>2683255</td><td>5.9318357</td><td>356.375198364258</td><td>12.612087</td><td>413.83007202148434</td><td>0.25039154</td><td>18.096925163269056</td><td>0.4938725</td><td>29.27996597290041</td><td>0.19070134</td><td>26.524150848388686</td><td>0.2986061</td><td>33.34638252258304</td><td>1.203017</td><td>98.89707946777344</td><td>0.95266587</td><td>111.70701904296901</td></tr>\n",
       "<tr><td>8</td><td>171728380</td><td>2683255</td><td>5.9318357</td><td>356.375198364258</td><td>12.612087</td><td>413.83007202148434</td><td>0.25039154</td><td>18.096925163269056</td><td>0.4938725</td><td>29.27996597290041</td><td>0.19070134</td><td>26.524150848388686</td><td>0.2986061</td><td>33.34638252258304</td><td>1.203017</td><td>98.89707946777344</td><td>0.95266587</td><td>111.70701904296901</td></tr>\n",
       "<tr><td>9</td><td>171728371</td><td>2683255</td><td>5.9318357</td><td>356.375198364258</td><td>12.612087</td><td>413.83007202148434</td><td>0.25039154</td><td>18.096925163269056</td><td>0.4938725</td><td>29.27996597290041</td><td>0.19070134</td><td>26.524150848388686</td><td>0.2986061</td><td>33.34638252258304</td><td>1.203017</td><td>98.89707946777344</td><td>0.95266587</td><td>111.70701904296901</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(\"$('#table140260366492000-201898').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",
       "    $('#table140260366492000-201898').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], 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": [
       "{'90prime_g', '90prime_r', 'mosaic_z', 'ukidss_j'}"
      ]
     },
     "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 XMM-13hr')"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAesAAAEgCAYAAACU3FvWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xl4VNX5wPHvO1t2QoCwL0EBWQWEoqKtO4J1L7buu2CtS4tW/VXF3Vata13qUtfWrbYigqJQN0SRRVRWZZElhEBYErJnlvP7486EISSTSchk7mTez/PwzHLv3PtOwuSd95xzzxFjDEoppZSyL0e8A1BKKaVUZJqslVJKKZvTZK2UUkrZnCZrpZRSyuY0WSullFI2p8laKaWUsjlN1sqWROQOEflnWztXWyUin4rI5fGOQ6m2SpO1ipqIrBeRShEpE5GtIvKiiGTGO65EJyIjRaRERPqFPTdKRIpFJC/4eL2I1IhIpzqv/VZETNh+LwUfn1pnv0eDz18cIY67RWSpiPhE5I46244JbisWkR0i8o6I9Ni/d66UipYma9VUpxhjMoFDgJ8Bt8Y5noRnjFkCPAk8JxY38AIw1RizPmzXn4BzQg9EZBiQVs8hfwQuCtvPBZwFrG0klDXAjcDMeratAE40xrQHugOrgacbOV69gvEopZpAk7VqFmPMZuADYCiAiFwiIitFpFRE1onI5NC+ItJJRGYEq7KdIjJXRBzBbTeJyObg634QkePCTpMqIm8Gt30jIsPDjnmziKwNblshImeEbbtYRL4Qkb+KyC4R+UlEJoRt7ysinwVfOxvoFLYtVUT+Gawei0VkoYh0qe9nICKDgs2/xSKyPLyaDVa4T4rIzOB5vhaRAyP8SO8EugGTgD8BZcATdfZ5Fbgw7PFFwCv1HOs94AgRyQk+Hg98DxRGOD/GmJeNMR8ApfVs22qMKQh7yg/0q7NbHxGZF3y/H4VaAUQkL1jVXyYiG4GPI8WhlNqXJmvVLCLSCzgJWBJ8ahtwMtAOuAR4REQOCW67HsgHcoEuWMnIiMhBwNXAz4wxWcCJwPqw05wG/BvoALwGTAtWnWBViT8HsrES3T9FpFvYaw8FfsBKxA8A/xARCW57DVgc3HY3YVVo8H420AvoCFwJVNbz/t1YSfEjoDNwDfCv4HsKOScYWw5W1Xpv3eOEGGOqgcuA+4M/r8uMMYE6u80H2gW/JDiB3wD19bVXAdOBs4OPL6T+pN4kItJbRIqxfh43YP1cw52L9bvvDHiC+4Q7ChiE9XtWSjWBJmvVVNOCf7C/AD4D7gMwxsw0xqw1ls+wktjPg6/xYlWNfYwxXmPMXGNNSu8HUoDBIuI2xqw3xoQ31S42xrxtjPECDwOpwGHB8/3bGFNgjAkYY97EapYdE/baDcaY54wxfuDl4Pm7iEhvrOb724wx1caYz7GSbogXK0n3M8b4jTGLjTG76/k5HAZkAn8xxtQYYz4GZhDWTA381xizwBjjA/4FjGjkZ7sM8AFLjTGrGtgnVF2fAKwCNjew3yvAhSKSjZUkpzVy7kYZYzYGm8E7YXV/1I3xRWPMj8aYSuAt9n2/dxhjyoPblVJNoMlaNdXpxpj2xpg+xpirQn94RWSCiMwPNnMXY1XdoeblB7Eqy4+CTeQ3Axhj1gC/B+4AtonIGyLSPexcm0J3glVmPlZ/KSJyYXBwVXHwfEPDzgdhTb7GmIrg3czg63cZY8rD9t0Qdv9V4EPgDREpEJEHwqr5cN2BTXWq3w1A+KCr8GbniuD5I3kI6wtQTxE5u4F9XsWqYC8mQrVsjPkCqyXjVmBG3QQZbLYvC/77eb0HafjYO7G+AL1bp/+5sfe7CaVUs2iyVvtNRFKA/wB/BboEq6/3AQEwxpQaY643xhwAnAJMCfVNG2NeM8YcCfQBDFYzcEivsHM4gJ5AgYj0AZ7DakLvGDzfstD5GrEFyBGRjLDneofuBCv/O40xg4GxWE37F7KvAqBXqO897DgNVboRBX8ep2E1u18JPCYiHeruZ4zZgDXQ7CTgv40c9p9YTer7JHVjzBBjTGbw39xmhOzCau5u14TX6BJ/SjWTJmvVEjxYzdlFgC84mGtcaKOInCwi/YJ9xruxmr/9InKQiBwbTPZVWH2h/rDjjhKRM4PV2++Baqx+2wysP/xFweNfQnCgW2OCyW4RcKeIeETkSKwvEKFYjxGRYcE+4d1YzeL+eg71NVAO3CgibhE5OnicN6KJI1zwi8NzwO+NMUXBQV6zgUcaeMllwLF1Wgfq8zhWc/nnUcbhFpFUrL8LruBgO2dw25nB35dDRHKxuiWWBKtspVSMabJW+80YUwpci9VPuQurmXZ62C79gTlYI5y/Ap4yxnyKleD/AmzHakLtjDX4LORdrEFUu4ALgDODle8KrCbjr4CtwDBgXhNCPhdrANpO4Hb2rjy7Am9jJeqVWM3S+wziMsbUAKcCE4LxPwVcGKGvOZL7gFXGmH+FPfd7YIKIjKu7c3BswKLGDmqM2WmM+Z+JftH657C+MJ0D3BK8f0FwWw9gFtZI8aVAADijnmMopWJAov8cK6WUUioetLJWSimlbE6TtVJKKWVzmqyVUkopm9NkrZRSStmcJmullFLK5hJu9ZtOnTqZvLy8eIehlFIJZfHixduNMbn7eYzOLpfreax5DbTYa1kBYJnP57t81KhR2+puTLhknZeXx6JFjV5iqpRSKoyIbGh8r8hcLtfzXbt2HZSbm7vL4XDodb8tKBAISFFR0eDCwsLnseZw2It+M1JKKRWtobm5ubs1Ubc8h8NhcnNzS2hgNkZN1koppaLl0EQdO8Gfbb15OWbJWkReEJFtIrKsge0iIo+LyBoR+T5s7WOllFJKhYllZf0SMD7C9glYc0b3ByYBT8cwFqWUUipqI0eOHBjvGMLFLFkbYz7HWiihIacBrxjLfKC9iHSLVTxKKaVUtJYsWdKcRXliJp6jwXuw92L0+cHntsQnHKWUUtH649vf9fqxsDS9JY85oGtWxYMTh29qaPsPP/zgGT9+fP8xY8aUffPNN5mDBg2quPTSS7ffddddPXbs2OF66aWX1g0ePLj6vPPOy9u4cWNKWlpa4Nlnn91w6KGHVs6cOTPz+uuv7w0gInz55ZerHA4H48eP71dSUuL0+XwyderUgvPPP78YID09fWRFRcUSgFtvvbXLW2+91VFEOO6440qeeuqpfdatX79+vXv8+PH9Q49Xr16dtnLlyqUDBgyoaYmfTTyTtdTzXL0DF0RkElZTOb17945lTEolPX/AsK20ihpfgCpvgO/zizlqQC6d26XGOzSl2LRpU+qbb765btSoURsOPvjgQf/61786Llq0aNVrr73W/t577+3Wo0ePmuHDh1fMmTNn7fTp07MuuuiivqtWrVrx0EMPdX388cc3jBs3rrykpMSRnp4eAJg5c+aaDh06BLZs2eI69NBDB5577rnFDseeRue33nqr3cyZM3MWL168KisrK7B161ZnfXHl5eV5V61atQLgz3/+c+7cuXOzWipRQ3yTdT7QK+xxT6Cgvh2NMc8CzwKMHj1aRyIq1QL8AUONL0Cax/rbU1Hj4y8frOLtxflU1Pj32f+iw/tw04SBpHv2/bPhDxgKiitpl+omO90d89hV/EWqgGOpR48e1WPGjKkEGDBgQOWxxx672+FwcMghh1Tcc8893Tdv3pzyn//8Zw3AqaeeWjpp0iTXjh07nIcddljZDTfc0OvXv/71znPOOWfXgQceGKiurpbf//73PefPn5/pcDjYtm2bJz8/39W7d29f6HyzZ89ud/7552/PysoKAHTp0mXfD0eYjz76KOOVV17JnT9/fos2o8czWU8HrhaRN4BDgRJjjDaBK9UKNhdXcsrfvqC82sfFR+TRu0M6T368hoKSKrq2S2VQt3ZU+/yMPbATC9fvZMnGYl7+agMvf7WBxbceT8fMFAC+21TMLdOWsnJLKf6A9T16aI923HP6MEb0ag9Atc/Ppp0ViAgH5mbG7T2rtsHj8dQWbA6Hg9TUVAPgdDrx+/3idDr3KehExNx3332Fp59+esm7776bPXbs2EGzZs36ce7cuRk7duxwLV26dGVKSorp0aPHsMrKyr3GchljEKmvIXhfGzZscE+ePDnv3XffXZOdnR3Yz7e6l5glaxF5HTga6CQi+cDtgBvAGPN34H3gJGANUAFcEqtYlFJ7e3PBRnaW1zCwaxbPfLYOgN4d0nlz0mEcekDHel/z0ryfuOO9FYy6Zw6nDu/OvDXb2VFutfKdMLgLHTM8lFb5mLd2O6c/OY/hPbNJ8zj5blMJlV6rGLn22H5MGXdQ67xJlZQOO+yw0hdffLHjgw8+uGXGjBlZOTk5vg4dOgSWL1+eMmbMmMoxY8ZUfv311xnLli1LLSkpcXbq1MmbkpJi3nvvvayCggJP3eONHz9+97333tv9iiuu2BlqBq+vuq6urpYzzzzzgLvvvnvzwQcfXN3S7ytmydoYc04j2w3wu1idXynVsP+t2saYvA68deXhrNyym5JKL6P75OByNnyByMVH9KVvbiaXvLiA6d9ZPVbnjOnNVUcfSK8Oe8YZlVX7ePKTNcxesZWCkipOGd6NQd3a8fGqbTz+8RoO6ZPD0Qd1jvl7VMnp/vvvLzj33HPzBgwYMDgtLS3w0ksv/QTwwAMPdP7yyy/bORwOM2DAgMqJEyeWFBcXOydMmNBv6NChg4YMGVLRt2/fqrrHmzhx4u5vvvkmfcSIEYPcbrc5/vjjS5544ol9BpjNmTMnY9myZRn33HNP93vuuac7wKxZs1bn5eV5W+J9iZUzE8fo0aONzg2uVPN5/QEG3TaLSb84gBvHN+9SUq8/gEDE5F5XjS/AuEc+I83j4v1rj4y6aVG1DBFZbIwZvT/H+O6779YPHz58e0vFpPb13XffdRo+fHhe3ed1ulGlksyGHeX4AoZ+nZvff+x2OpqUqAE8LgeXHdmXlVt2s3JLabPPrVQy0mStVJIpKLZa+sKbrlvLhGHdcAj8e3FcBhIr1SIuuOCC3gMHDhwc/u+xxx6rf7BHC0m4JTKVUvtnW6k19qVzVkqrn7tTZgpH9OvER8u3MvXkwdoUrhLSq6++urG1z6mVtVJJpiiYrDtltn6yBjh1eHc2F1eycP0uqrwRL1lVSgVpslYqyWwrrSLD4yQjJT4NaycN60ZWiotfP/MVh//5f5RUtshgWaXaNE3WSiWZ7WU15MahCTwkI8XFI78ZQV7HdHZVeFnwU6T1fpRSoMlaqaSzu9JLdlp8pwQ9fnAX/nvVEQDk76qIayxKJQJN1kolmbJqH5mp8R9bmpPuxu0Utu5u8cmeVJIZM2bMQZ9//vk+lzccddRR/bZv377PwhtTpkzpPnXq1C77e97gSlsH7O9xohH/T6xSqlWVVfnolNn6l23VJSJ0yPCws1yTtYqNzz77bE0sj5+Xl+edNWvWulieI0STtVJJpqzaR2aKPVbGykhxUV7PCl8qAUz7XS+2rWjZb32dB1dw+pMR17M++eST+69evXo5wNSpU7uUlZXVVs5+v5+zzjorr2fPnjWPP/54QY8ePYYtWrRoZbdu3Xw33XRT1zfffLNT9+7dazp27OgdOXJkBcA999zT+cUXX8x1Op1mwIABVTNmzFhX39rXOTk5+yzMUTeeWNJkrVSSKa3ykmWDZnCAzBQX5dW+xndUqhFer1dOP/30voMHD668//77C8O3zZ07N/2dd97psHTp0hVer5cRI0YMDiXrxx9/vOuGDRuWpqWlmVCTeUNrX8eTPT6xSqlWYYyhrNpnm2Sd7nFSUa2VdUKKUAHHw1VXXdXn9NNP31k3UQN88sknmSeddFJxaE3qcePGFYe2HXTQQZVnnHFG31NPPbX4vPPOKwaob+3r1nsn9dMBZkolkUqvn4CxKlo7yPC4KNPKWkXJ5XKZQGBP3qyqqqrNYaNHjy6bO3duu4qKinqnxWtotrxPPvlk9e9+97uixYsXZwwfPnyw1+vlvvvuK3z++ec3VFZWOsaOHTtoyZIlqS39XppKk7VSSSR/VyWALUaDg9VnXVGjyVpFp2fPnr6dO3e6CgsLnZWVlfLhhx9mh7ZNnjx5+7hx40pOPvnkA73evSfaOfbYY8tmzpzZvqysTHbt2uWYPXt2e7D6uNeuXes55ZRTSp966qn80tJSZ0lJiTO09vW9995bOGzYsPJly5bFPVnb4xOrlGoVZ/39K8BqfraDdI+TCh1gpqKUkpJirr/++i1jxowZ1LNnz+p+/frttf70HXfcsfUPf/iD88wzz+w7bdq0n0LPH3nkkRVnnHHGzqFDhw7p0aNH9ZgxY8oAfD6fnHvuuX1LS0udxhiZPHny1k6dOvmvv/767nXXvm4oJhFplXWmNVkrlURCU3s6HfZoVPO4HNT4494dqBLIrbfeuu3WW2/d1tD2Rx55pCB0f/PmzUtD9++///7C+vqzFy9e/EPd515++eWo+uO3bdvmys7ObpVvm/b4xCqlWlWGTSprj9NBjU+TtUo8n3/+efoFF1xwwNVXX721Nc6nlbVSSSjNbZNk7dJkrexvwYIFaRdeeGHf8Oc8Hk9g/fr1y1orBk3WSiWhFBsla1/AEAgYHA5d21rZ05gxYypXrVq1Ip4xaDO4UknogE4Z8Q4BsJI1oP3WSjVCk7VSSeS0Ed3J65hOToYn3qEAkOKyKvxqbQpXKiJN1kolEa8/gMtpn499bWWtyVqpiOzzqVVKxVyNz+C2UbJOcWozuFLRsM+nVikVc75AAI/TPgO5tLJWTXX33Xd37t+//5B+/foNueuuuzoDbN261Tl27Nj+ffr0GTp27Nj+RUVFTR5BOXLkyIEtH23L0WStVBLx+gO2qqw1WaumWLhwYeorr7yS+80336xcuXLl8lmzZrVfunRpyu23397t6KOPLt2wYcOyo48+unTq1Kldoz2mz2dNd7tkyZJVMQu8BeilW0olEa/P4LJRZR364uDVZvCEc9u823qt2bWmRdez7pfTr+LuI+5ucPawpUuXph1yyCFlodWzjjjiiNI333yz/axZs9p/9tlnPwBMnjx5x1FHHXUQsHnKlCnd161bl1JYWOjesmWL59prry28/vrrt8+YMSPr7rvv7ta5c2fvihUr0teuXbs8PT19ZEVFxZIZM2Zk3Xnnnd1zc3O9K1asSD/ppJN2DRs2rPKpp57qUl1dLe+8887aIUOGVBcUFLguueSSPps3b/YAPPzwwxvHjRtXXl/cBQUFrokTJ/YtLi52jRgxouLTTz9tt3jx4pXdunWLemJ8+3zFVkrFXI3NKmtX8NpqX6BVpldWCW7EiBGVX3/9dVZhYaGztLTUMXv27OxNmzZ5duzY4erTp48XoE+fPt6dO3fWFqIrV65MmzNnzur58+evevDBB7uvX7/eDfD9999nPPjgg5vXrl27vO55Vq1alfb0009vWrly5fK33367448//pi6dOnSlRdccMH2hx56qDPA5MmTe02ZMmXrsmXLVr7zzjtrr7zyyryG4r755pu7H3XUUaUrVqxYeeaZZ+7asmVLky/H0MpaqSRi9VnbJ1k7g8naH9DKOtFEqoBj5ZBDDqm67rrrCo899tgB6enpgcGDB1e4XJHT2IQJE4ozMzNNZmam7/DDD989d+7cjJycHP/BBx9cPnDgwJr6XjNs2LDyUPLv3bt39YQJE0oAhg8fXvnZZ59lAcybN6/d6tWr00KvKSsrc+7atcuRk5Ozz3/mBQsWZE6bNm0NwMSJE3e3a9euyfOJa7JWKol4bTYaPNQk7/NrZa2i84c//GH7H/7wh+0AV199dY+ePXvWdOzY0bdhwwZ3nz59vBs2bHB36NChtnm57jrWocfp6ekNfkNMSUmp/Q/pcDhITU01oft+v18AjDEsWrRoZWZmZqP/eY3Z///f9vnUKqVizrrO2j591q7g6l/aDK6itXnzZhfA6tWrPTNnzmx/2WWX7TzxxBOLn3nmmY4AzzzzTMfx48cXh/b/4IMP2ldUVEhhYaFz/vz5WUceeWS9/cpNdeSRR+6+//77O4cef/nll2kN7TtmzJiyV199tQPAf//733a7d+9u8mj1mFbWIjIeeAxwAs8bY/5SZ3tv4GWgfXCfm40x78cyJqWSmTdgrz5rp/ZZqyY69dRTDywuLna5XC7z6KOPbszNzfXfeeedW84444wD+/Tp06l79+4106ZNWxvaf+TIkeXHHXdc/4KCAs8NN9ywJS8vz7ts2bLU/Y3j2Wef3XT55Zf3HjBgwGC/3y+HHnpo6dixYzfWt+9f/vKXgokTJx4wePDgnMMPP7wsNzfX2759+yY1hccsWYuIE3gSOAHIBxaKyHRjTPhk6LcCbxljnhaRwcD7QF6sYlIq2fn9pnZQlx24ndpnrZqmvvWnu3bt6v/qq69+rG///v37V73++usbwp87+eSTS08++eTS8OcqKiqW1LdtwYIFP9T3um7duvlmzpy5LpqYO3To4P/8889/dLvdzJkzJ2PevHlZaWlpTfqGGsvKegywxhizDkBE3gBOA8KTtQHaBe9nAwUopWLGG7DXpVuhytqrfdaqDVuzZo3n17/+9YGBQAC3222eeeaZ9U09RiyTdQ8gfLRgPnBonX3uAD4SkWuADOD4GMajVNLzB0xtP7EdhGLxazO4ioGHH364VQvAxx57rOPTTz/dJfy5n/3sZ2WvvvrqxpUrV+7XEpuxTNb1fX2v+4k8B3jJGPOQiBwOvCoiQ40xe7WJicgkYBJA7969YxKsUsnA6w/UVrN2oH3Wqi257rrrdlx33XU7YnHsWH7Fzgd6hT3uyb7N3JcBbwEYY74CUoFOdQ9kjHnWGDPaGDM6Nzc3RuEq1fb5A6a2n9gO3LWXbmmftVKRxDJZLwT6i0hfEfEAZwPT6+yzETgOQEQGYSXrohjGpFRS8/kNThs1g2tlrVR0YvapNcb4gKuBD4GVWKO+l4vIXSJyanC364ErROQ74HXgYtMSV48rpeplXbpln8pa+6yVik5Mr7MOXjP9fp3npobdXwEcEcsYlFKWQMBgDLbqs3ZpM7hSUbFPe5hSKqa8wWuZ7TQpii7koZpK17NWSrVpoaZmO1XWexby0GStGtcW1rP2er3Nep0u5KFUkghNPGKnGcxCfdY6KUriKfjTLb2qV69u0fWsU/r3r+h+371tbj3rKVOmdN+yZYt748aNng4dOvjee++9n5r6s9FkrVSSCFWvtkrWOt2oaoIRI0ZU3nXXXT0KCwudGRkZZvbs2dnDhw8vb2w968WLF68sLS11jhw5cvCvfvWrErDWs16yZMny+pbJXLVqVdrbb7+9rnPnzr4+ffoMS0lJ2b506dKVd999d+eHHnqo8wsvvLAptJ71iSeeWLZ69WrPiSee2H/dunX7rI0d8v3336d//fXXq6JZpas+mqyVShKhQVwuG/VZO0X7rBNVpAo4VhJ1PWuA8ePHFzc3UYMma6WSRigh2unSLYdDcIiuZ62il4jrWQNkZGTsV/ORfb5iK6ViKpQQ7TQpCliVvlbWKlqJuJ51S9DKWqkk4au9dMs+lTVYfejaZ62ilYjrWbcESbQJw0aPHm0WLVoU7zCUSjg/bi1l3COf88S5Izn54O7xDqfWsDs+ZOKontx+ypB4h9KmichiY8zo/TnGd999t3748OHbWyqmWJsyZUr3zMxM/1133bU13rFE67vvvus0fPjwvLrP26s9TCkVM77aS7fs9bF3Ox3aZ61UI7QZXKkkEWoGt9OlW2BNjKJ91ioW7LSe9f4eW5O1UkkilBBd2metmi8QCATE4XDot6t67O961oFAQIB6Pwz2ag9TSsWMXZvBnQ7RZvDEsayoqCg7mFRUCwoEAlJUVJQNLKtvu1bWSiWJ2mZwO1bWCTbQNVn5fL7LCwsLny8sLByKFnstLQAs8/l8l9e3UZO1UknCZ8O5wUH7rBPJqFGjtgGnxjuOZKTfjJRKErVzg9toulGwkrVfm8GVishen1qlVMx4/XYdDa4zmCnVGE3WSiUJv41Hgwe0z1qpiDRZK5UkvDZcIhO0z1qpaGiyVipJ+GsnRbHXx16vs1aqcfb61CqlYsZbu+qWvSprh15nrVSjNFkrlSRCCdFts9Hg2metVOPs9alVSsVMqKnZbpW19lkr1ThN1kolCW9tZW2vZG31WWuyViqSqGYwE5GuQO/w/Y0xX8YqKKVUywslRFtW1tpnrVREjSZrEbkPOB9YBfiDTxvgpBjGpZRqYd5gM7jd+qyd2metVKOiqax/BQwwxlTFOhilVOz4bToa3KUzmCnVqGi+Yv8U5X5KKRuz86Qo2metVGTRVNalwBIRmQNUh540xkyJWVRKqRbnDwRwOQQR+yVrn06KolRE0STrWcF/SqkE5vMb2zWBg666pVQ0Gk3Wxph/iIgL6Bd8ao0xxhfbsJRSLc0XMLYbXAbBS7d0gJlSETX6yRWRnwNrgH8ALwA/isgR0RxcRMaLyA8iskZEbm5gn1+LyAoRWS4irzUleKVU9Hz+gH0ra+2zViqiaJrBHwFOMsasABCRQcCrwOhILxIRJ/AkcAKQDywUkemh4wT36Q/8H3CEMWaXiHRu3ttQSjXGqqztmax1NLhSkUXTJuYJT7DGmJWAJ4rXjcFqMl9njKkB3gBOq7PPFcCTxphdwWNviy5spVRTaZ+1UokrmmT9jYg8IyJHBv89DSyJ4nU9gE1hj/ODz4UbAAwQkXkiMl9Extd3IBGZJCKLRGRRUVFRFKdWStXlCxjbLY8J2metVDSi+eReCawFbgRuAtYBk6N4XX1f4et+Il1Af+Bo4BzgeRFpv8+LjHnWGDPaGDM6Nzc3ilMrperyBQK4bNkMrpOiKNWYaEaDVwEPBP81RT7QK+xxT6Cgnn3mG2O8wE8i8gNW8l7YxHMppRphVdb2S9a6kIdSjWuwshaR14O3S0Tkm7r/ojj2QqC/iPQVEQ9wNjC9zj7TgGOC5+mE1Sy+rjlvRCkVmc8fsGUzuCOYrI02hSvVoEiV9R+DtxObc2BjjE9ErgY+BJzAC8aY5SJyF7DIGDM9uG2ciKzAWiTkj8aYHc05n1IqMn+6ICghAAAgAElEQVTA2LIZPFTtBwzYMDylbKHBZG2MyQ/eLQCqjDFGRA4EDgI+iubgxpj3gffrPDc17L4BpgT/KaViyOu3ZzN4aIS6LxDA6XDGORql7CmaNrG5QJqIdAM+A36LNTmKUiqBWJW1/ZrBQ18gtN9aqYZF88l1GGMqsJbKfMIYcwpwcGzDUkq1NK+NZzADdES4UhFElaxF5GfAucCM4HPaVqVUgvHbeAYzgIAma6UaFE2yngLcCcw0xiwTkQOwmsaVUgnEGzA4bTga3KWVtVKNiuY664+BjwHEWgh3qzHmqlgHppRqWT5/ALctm8GtLxDaZ61Uw6JZdesVEWknIunAcqzJS3T0tlIJxh+w69zg1q1W1ko1LJo2sWHGmN3A6ViXbPUELo5lUEqpluf1B2y5nnVtZa2LeSjVoKhW3RIRF9aKWdOCK2gFYhuWUqql2bWyrr10S2cwU6pB0axn/TywEVgGfCYivYGymEalVBzUbNpE2dy5OFJSaffLk3Ckpu61vXz+10iKh/SRI+MU4f7x+u05g5mz9jprrQGUakg0A8weAR4JPRaRTcCxsQxKqdZQ+r//UfrJJ/h3FWO8NVQuWkygogKA6rVr6XKjNeNuycyZeHr2ZOPFFwMwYMHXONu1i1fYzea38UIeoH3WSkXSYLIWkXOMMa+LyLUN7PJ4jGJSKmbKvpiHq2MH/LtLyb/6GsTtxtWlC46MDDKPP44OF11E/tXX4C2wFojz795NwfU37HWMHc89T+frE2+MpbVEpv36rB2hZK191ko1KFJlnRO81QWkVcIzxrDzpZfZdv/9tc95+vQh7603cWZn77WvMzsbU1MDwO4PZtU+n3XCCSDCrjffpNPVv8ORktI6wbcQX8DY8tKtPQt5aLJWqiGRFvJ4Knh7W+uFo1TLqly2nOoffqB83hfsfv8D0g45hLQRIxCnk46XX7ZPogYQjwfj9bJ71iwKb7+d1MGD6Tr1NlIPPpjyL7+i9KOPKJ/3JVnHHhOHd9R8Pr89J0XR6UaValyjfdbBAWVXA3nh+xtjzoxdWEo1X9WKFfhLy0gbOoT1E4MrvDoc5F53LR0nT0YaSVjidmNqaih6/G8AdL3zTtKGDQUgpX8/AHxFRbF7AzHiCwRsOd2oSydFUapR0YwGnw68AsxGL9lSNuYvK6fo4YfY9drr+2zrdvfdtP9VdN8vxePGu7kA78aNdL75ptpEDeDIyAAgUJZ4F0RYlbX9knXou5P2WSvVsGiSdY0x5uGYR6LUfghUV7Pp8sup/PZb2p/9Gzx98tj5wgt48vLIve5a0kaNivpY4nbXDjBzd+++1zZHerp1vvLylgu+FRhj8Nl2iUwrJu2zVqph0STrv4nIrcCHQHXoSWPM9zGLSqkmCFRUsHnK9VR++y09Hn6IdiedBECHiy4EYxBn0xaJE48HfD4AXJ067b3N4cCRnp5wyTrUxGzHS7e0z1qpxkWTrAcAlwMT2NMMboBfxCoopRoTqKjAV1SEu2dPNk+5nrLPP6fr7VNrEzXQaN90Q8Ttrr3v6thxn+2OzEzKFy6gatUqUgcObNY5WlsoEdpxUhSXToqiVKOiSda/BvKMMdWN7qlUDNVs2sSu19+g8ttvqfzuO/D7a7d1ue1Wcs45p0XO4/B4au+7cve9ctGRkUH1ipX8dPoZDFq1skXOGWu+RKistc9aqQZFU3p8D2TFOhClGmL8frbe/wDrfnkyu159FVNTYyXmYAXc7pRTyDn33BY7X82GjbX3Q33UexH7JbzGhBbJcNn40i0dDa5Uw6KprDsCq0Tka/bus9ZLt1TMebdsYc0x1uy27X75Szrf+EfcXboAUP7FF9SsX0/2aachLZhAA+WRR3r7S3e32LlaizfYxGzrZnAdYKZUg6JJ1vfGPAql6iif/zW7Z31AxVfzAci54AK63vKnvfZJGzmSmvXrSenfv0XP7S+1knXev9+qd3uXG2+i4I9/xFlPf7ZdhapWO166pZW1Uo2LJll/CVQZY4yIHAgchLWutVItzhjD9qeeYvvfnsCRmUnKwIPofu01ZP/yl/vs23XqbeScczbuLp1bNIZAaSkArmAFX1f2KSdT8c1iSmd92KLnjSWv36qs3TZuBtc+a6UaFk2yngv8QkSygc+AJcDZwIWxDEwlH29BAVsffJDSD2aRfdppdL3rzojzbzvS0kg7+OAWjyN12DAq5s/HlZPT4D6OtHQClZUtfu5YCSVCOzaDa2WtVOOiSdYOY0yFiFwKPGGM+YuIfBvrwFTb5S8rp/jNNzFeL9lnnEH+VVdRtXx57fbOf7yBDpde2qL90E3R82+PU7Nx416XcNXlSEvDVFVh/P4mX8cdD77aPmv7Vda1041qn7VSDYoqWYvIz4BzgUnB5+z/10nZUqCigvW/+Q01a9cCUPToo7XbJD2dno89SubPfx6v8ABwZmWRNmRIxH0c6WkABCqrcGZmtEZY+yUhLt3SylqpBkWTrKcAdwIzjTHLROQArKZxpZps+zPPUrN2Lb2efYaq5cspesxaFv2A998n5YC+cY4uepJmJWtTWQGJkKz99k3WtaPB/TopilINaTRZG2M+Bj4Oe7wOuCqWQam2qeS999jx7LNkn3Yamb/4BRmHHUagqprMI49IqEQNVp81kDD91rUDzGzYDO7QylqpRkWzRGY/rOo6L3x/Y8y42IWl2pqyzz6j4Ob/I33MGLreeQdgzcHd+Q+/j29gzVS7oEeCJOtEmG5UF/JQqmHRNIO/DfwD+Cfgb2RfpfZR8t57bPnTLaQcNICeTz6JIzU13iHtt9o+64qKOEcSnVAzuJ2vs9bKWqmGRZOsA8aYv8U8EtUmVXzzDQU33Uz66NH0ePSRhBiMFQ1HbZ91olTW9m0G39NnrclaqYZE88l9V0QmiUiuiLQL/Yvm4CIyXkR+EJE1InJzhP0miogRkdFRR64Swo5nnsXVsSO9nn4KV4cO8Q6nxYQGmCVMM7iNB5hpZa1U46KprC8P3t4W9pwBekd6kYg4gSeBE4B8YKGITDfGrKizXxZwLfB1tEGrxBCorKT8q6/IOedsHBlto6IOqR1glijN4MFEaMfKWkRwiPZZKxVJNKPBezXz2GOANcHR44jIG8BpwIo6+90NPADc0MzzKJvyFhZiampIHTw43qG0uD191olSWVvN4HbsswZrYhStrJVqWDSVNSIyEBgM1I4MMsa81sjLegCbwh7nA4fWOe5IoJcxZoaINJisRWQSwQlZeveOWNArGwmUlQPgyGp7K6xKcBpUU50Yy7x7aytreyZrp0N0ulGlIojm0q1bgXHAQOBD4ETgC6CxZF3fX4XaT6OIOIBHgIsbi8EY8yzwLMDo0aP1E50gAmXWghiOzMw4R9LyQlORGp8vzpFEJ1RZ23E9a7CStS7koVTDovnk/gY4BthijLkAGE50FXk+EN6E3hMoCHucBQwFPhWR9cBhwHQdZNZ2+MuspSadbTFZu6yPQMIkaxtfZw2hylpnMFOqIdEk60pjjB/wBQeDFQIHRPG6hUB/EekrIh6slbqmhzYaY0qMMZ2MMXnGmDxgPnCqMWZRk9+FsqU23Qxem6y9cY4kOntGg9uzsnY5RBfyUCqCaD65S0SkPfACsAhYAHzT2IuMMT7gaqym85XAW8aY5SJyl4icuh8xqzj6qeQn3ln9TlT7htaFbovN4IRW2kqYyjq06padK2tN1ko1JGJztlhrFN5hjCkGnhSRD4F2xphGkzWAMeZ94P06z01tYN+jo4pYtbqdVTsByEnJYeq8qXy//XtO73d6o0tYBsqDzeBt7LItsC43Ercb402MytobrKzdNq2stc9aqcgiJmtjjBGRGcCo4OM1rRKVsgVvwMvLy1/miSVPIAjH9D6Gb4uspcyr/FWkudIivt5fVoakpCAeT2uE2/rcbow3MSrrUH+wUytrpRJSNAPFFojIIdFW0ypxBUyAeZvnUVBWQHF1MW/98BbbKrdxfO/jWVeyjtkbZtfuW+VrPFkHSsvaZhN4kLhcCTPAzGvjGcxA+6yVakyDyVpEXMF+5yOBK0RkLVCOdUmWMcYc0koxqlZQWlPKdZ9cx8LChbXP9c3uy1/H/JVxfcZR7a/mtVWvsXrXamasm0Glr5IcciIeM1BWhqONzAVeHytZJ0YzeKiJ2Y4zmEGwGVwra6UaFKmyXgAcApzeSrGoOFmxYwW3fHEL60vWc9thtzGm6xi6ZnQl1bVndaxUVyqXDr2UWT/Nqk3WjfGXl+HMbHsjwUMSqc86NMDMzjOY6UIeSjUsUrIWAGPM2laKRbWyd1a/w2urXmPVzlW087TjqeOf4vDuh0d8TajpO5pknQzN4IkyGtzrN3hc9qyqARxaWSsVUaRknSsiUxraaIx5OAbxqBjyB/y8vup1XlnxCk5xkl+Wz0E5B3HzmJs5+YCTyU7JbvQYoWo7qmRdVoa7d3Onlrc/cbkSZoCZ1x/AY9MmcLD6rHUhD6UaFilZO4FM6p82VCWQbRXb+Gj9R8xYN4PlO5YzPHc4boebU/udyuXDLsftcEd9rCZV1mVlODPabmWNO5EGmAVsOy84aJ+1Uo2JlKy3GGPuarVIVItbvmM5b6x6g9kbZlPuLad3Vm/u//n9TOg7odFrpBvSlGTtL2vrzeDuBEvW9q6sdbpRpRrWaJ+1SixevxeD4Y1Vb/DoN4/iDXg5ttexXD3yavrn9N/v44eS9bfbvuXEvBMj7huoqGhz61iHS6QBZjU+Y+tk7dBJUZSKKFKyPq7VolD7rcZfwzPfP8MLy17AF7CqvbHdx3LvkffSKa1Ti50nlKz/ufKf3DTmpgb3M4EA+Hxtd0IUEuvSLbs3g7scQo1PK2ulGtJgsjbG7GzNQFTzFJYXUlZTxv998X+s2rmKo3oexcG5B9Mzs+d+NXc3pENqh6j2MzU1AG0+WZNAA8zsXFk7dVIUpSKKZgYzZUObyzbz9LdP8+7adwHIcmfx4FEPMq7POBwSuz/KIsK5A89lxroZEfcLNQ+H1n1ui8TtIlDReN+9HXj99m4Gd+l0o0pFpMk6gRhj+Cz/MxZvXcx7a9+jpLqE8wedT5f0LhzT+xj6tOvTKnG4HW68gcjNv3sq67abrEmo6UYDuG18nbUu5KFUZJqsE4A/4Gd3zW7+uuivTF87HY/Dw5BOQ7jl0Fs4qMNBrR6Px+nB6482WbfhZnC3J2EGmFnXWdu3z1oX8lAqMk3WNhcwAa6ccyXzt8zHKU6uGHYFvx3+W9zO+FWsbqcbn/ERMIEGm9yTohk80SprWzeDO7TPWqkINFnbkDGGeQXzeGHZCyzfvpwKXwUHZh/In3/+ZwZ1HBTv8GonUfEGvKQ4U+rdJ1RZO9pyZZ1Ao8Fr/IZ0j32TtVbWSkWmydpmqnxV3PT5TXy86WM6pHbgxLwTGZ47nDP7n9niI7uby+OwEnCNv6bRZE0br6wTZjS4z+6VtdQuNqKU2pcmaxtZvWs1t395O8u2L+P6Uddz7qBz8TjtV5mGmuAjDTILNYO36crakziTonj9ATwue3zZq4/DIbrqllIRaLK2gWXbl/HmD2/y0fqPSHGmcP8vrClB7Sq8sm5IMgwwS7jR4DavrLXPWqmGabKOI2/Ay5+//jP//vHfABzX+zhu+tlNdMvsFufIImtKZd22B5gl0tzgBpfDvsla+6yVikyTdZzsrNrJLV/cwhebv+CMfmdw/uDzGZAzIN5hRSVUWUe6fCuQBJV1Io0Gr7F5M7hLV91SKiJN1q3MGMP0tdN57JvHKK4u5o+j/8iFQy6Md1hNEj4avCHJ0AyeSMna7s3g2metVGSarFvR9srt3PnlnXya/ylDOg7h3iPv5fDuh8c7rCYLNYNH7LNOhmZwtxu8Xowxthmp35BEGA2ufdZKNUyTdSsImABv/fAWf1vyN6r91dz4sxs5b9B5MZ3DO5aiq6yDybotV9bu4MfH57P9JWregL3nBnc6HNoMrlQEmqxjzB/wM3n2ZL4u/JpDux7Knw77EwdkHxDvsPZL6HKymoCOBgcwPp+tWxCMMbafblQX8lAqMk3WMVThreDsmWfzU8lPXDX8KiYPn5yw1XS42so6wgCzpGgGd1nvze791v6AwRhsXVk7gsk6EboUlIoHTdYx9PzS5/mp5CduGH0DFwy+oE0katDKOkTCKms78wYHbtl51S2Xw0rQ/oDBZeMWAKXiRZN1DFT7q7nrq7uYvnY64/PGc9GQi+IdUouKqrJOhmQdbDUI9c/bVY3fmsbTzpW1M5SsjdE/SkrVQz8XLSx8xPd5g87jt8N/G++QWlztddbRTIriarv/xfYMMLN3svYGk7Xd+6wB7bdWqgFt9y9pHMzfMp8bPruBkuoSbhh9Q5urqEOiunSrpgbcbsTGs2btr8RpBk+cylpHhCtVP03WLWTljpX8dvZv6ZbZjSeOfYIRnUfEO6SYiXZSFEcbHlwG7DUa3M68vmCfdQIk66SZGCXgh6oSSO8Q70hUgojpp1dExovIDyKyRkRurmf7FBFZISLfi8j/RKRPLOOJlXUl67jqf1eRk5rDaye91qYTNYQNMGtkUpS2PBIcwvqsbb7yVqjP2s4Dt1zJVFnP/zvc1QEe6AsP9gdfdbwjUgkgZslaRJzAk8AEYDBwjogMrrPbEmC0MeZg4G3ggVjFEyvritdx6axLMcbw/LjnaZ/aPt4hxVxUlbW3pk0PLoOwS7dsvqb1nj5r+1bWnuBI9dAXizbjx4/goYGQv8h6vLsAZt20Z3v5NnjqcCj4Fu7IhkeGgq7rreoRy0/vGGCNMWadMaYGeAM4LXwHY8wnxpiK4MP5QM8YxtPilmxbwnnvn4eI8ML4FzigfWJPdhKtULKOfOlWMlTWoWZwe1fWidBnneJyAlDja0OJavHL8NpZULoFnj8OijfBw4P2bL+9GE76K+xcC88eZT1XsgnuyolPvMrWYvnp7QFsCnucH3yuIZcBH9S3QUQmicgiEVlUVFTUgiE2366qXdz+5e1kebJ4dcKrCT8rWVM4HU6c4mxkUpRkqKzDphu1sdpkbePrrEOVdbXPH+dImun9G+GfEyF/sfV47Scw4w977/Po0D33T7gbRGDMFfUfr2RzbOJUCSuWA8zq6yCrt0NKRM4HRgNH1bfdGPMs8CzA6NGj49qpVVJdwhur3uDdte+ytXwrTx7/JD2zEqpBoEV4nJ6IzeCBmuRJ1nYfYFYTGmDmsG+fdUooWXsTsLJ+sB+UB4uINbNhyJnw4yzIPQgunQXPHQs71uz9miOu3XP/jhL46FZY8Bz4qqzntq2E7Ei1jUo2sfyqnQ/0CnvcEyiou5OIHA/cApxqjLH1SIsFWxZw2rTTePLbJ8nyZPHcuOc4rNth8Q4rLlwOlw4wS5ABZqFqNcXtjHMkDattBk+0Pusl/9yTqI/4PfQ/EdZ9CnlHwgXvQGo2XLN479dc882+xxl3D9y6dc/jih0xC1klplhW1guB/iLSF9gMnA2cG76DiIwEngHGG2O2xTCW/ZZfms9V/7uKnpk9+fsJf2dgh4HxDimuPI7IlbVJgsqaBJkbvMprJes0GydrT6JU1r4a8JZDWg6UFcG7v7Oen7IK2nWL7hip2Q1vO/M5+O8VULF9/2NVbUrMkrUxxiciVwMfAk7gBWPMchG5C1hkjJkOPAhkAv8OTt6/0Rhzaqxiaq5tFdu48fMbEYRnTniGLhld4h1S3Lmd7kYmRfG2+WRdO8DM5qPBK0PJ2mPfZJ1i9z7rgm/3DAKrT2OJeuKL8PYl1n1HhD+7Q39lJeu5D8Phv2t6nKrNiumkKMaY94H36zw3Nez+8bE8f0tYvHUx1396PRW+Cv788z9rog5qtLL2enGkpbViRK1vT5+1vZvBK2usajURKmtbjgb/+hn44MaGt0/d1fgxBv5yz31nhC+xjuDvSCtrVYd9h4fGmTGGf638F5d/eDmZnkxeO+k1ju9j++8WrcbtcCd9M3iijAavTIBm8D2Vtc2StTH7JurDrtpz//z/QjRT6rpS9tx3tu2xHCo2dLrRevgDfu746g6mrZnG0b2O5r4j7yPLkxXvsGzF4/Q0uupWm0/WCTLALNRnneqx73fz0OA321XWJfl7P76jxLo95hbIXwgHHtP0Y0ZqBg8XCET3RUAlBf2fUIc/4OeWebcwbc00Jh08iceOeUwTdT3cDnfESVECVZU4UlNbMaI4SJBLtypr/DjE5jOYOW3aZz33oT33T7h7z/2UzOYlarCur47kxPus26ri5h1ftUlaWYdZsm0Jjy5+lG+2fcM1I69h0sGT4h2SbbmdjTSDV1TiSG/jfdbuxJhutMrrJ83tRBpLEnGU4rZpM3hooY0b1kBmbiuds6N1W7lLF/pQtez7VbuVfZ7/ORd9cBGbSjdxx+F3aKJuhMfhiTgaPFBVhaQmS7K2dzN4pddv65HgYOM+61BlndFp/48lUf4O0oIJumRT5P1UUtHKOujl5S/TNaMr006bRro7Pd7h2J7b6aamuv5kbQIBTFVVm28Gr03Wdm8G9/prJx2xqz3N4DZL1iEt0SpxzWLYsbbx/ULV9LL/wgFH7/95VZugyRpYvn05CwoXcN0h12mijlKkS7dMtTURXZtvBg/1WXsbbmGwg6oEqKxFBI/LYb8+65R2cNCEljlWh77Wv8Z0P8S6zeraMudVbYI2gwP/WPYP2qe05+yDzo53KAkj0qVbgcpKgLbfDO50gtOJqbF5M3iN39aXbYWkOB32Gg1eUwHVu6HTgNY9b2gE+Mb5rXteZWtJn6yr/dV8sfkLTsw7kUxPZrzDSRhup7vBS7dMMFk70tp2MziAeDyJ0WedCMna7bBXM/ju4FIG2XFaqGf93PicV9lS0jeD/+fH/1Dpq+SYXs28DCNJRbp0K1BlrRzU1mcwA6vf2u7JuqLGT066/a95T3E5a68Jt4XdwWUq28Vh9ateh9HAIoUqSSV9Zf3KilcY1WUUY7uPjXcoCSXSEpmBiuRoBodQsrZ3n/XuSi/ZafafNSsjxUlFtY2SdekW67Zd99Y/d7tuULGz9c+rbCupk7UxhqKKIoZ1Gmbra1DtyO1oeCEPU5VEzeAJUFnvrvLRLs3+jWgZKS7Kqm00sj7UDB6PgV7pHXV+cLWXpE7WFb4KagI15KTmxDuUhBOxsk6mZnCPx9YDzIwx7K700i7V/pV1pt2SdekWSMkGT0brnzu9I1QWg99GPw8VV0mdrHdVWavl5KRosm4qt8ONL+AjYPYdEFTbDJ4MydrmlXVFjR9fwNAuAZrBM1NclNsqWRdCVpxW2UvvBBidclTV0mQNdEjVKf2ayhNc5q++6rq2GbyNT4oCwWRdXY0J2GgUc5jdVdbvRyvrZijfDhmd43Pu0MQo5doUrizJnayrrWTdPrV9nCNJPG6H9ce/vsu3ApVWM3iyVNZln3zCD6NGxzuUeu2utJKf9lk3Q/m21psPvK7Q/OAVO+JzfmU7SZ2sd1ZZoy07pGhl3VShZF3f5VuB2uuskyBZB5cBDV1bbjehyjoRRoOHmsGNscklS2VF8ausQ3ORa7JWQUmdrGv7rHWAWZO5nQ1X1snWDG5nuysTqBk81UXAWJO4xJ23CqpLICPelbU2gytL0idrt8NNhjsOoz0TnMdhVZT1V9ZV4HbbPpG1BLu3HtT2WSdAZZ2RYjXV26IpPJQktRlc2URSJ+udVTvJSc3Ra6ybIdIAs0BlZVJU1QDuXnGaijJKO8ut30/7BEjWWcFkXVplg2Rdts26jVczuCsFPFk6MYqqldTJelf1Lh0J3kyhPuu1xfsu+WeqkidZe3r3iXcIERWWVJLictA+3f7JukOG9QVwZ7kNZoSrnRAlTpduAbhTYct38Tu/spXkTtZVu/Qa62ZKdVnJ+KFFD+2zLVBZhbTx5TFDPH16xzuEiApKqujePi0hWo9ys1IA2F5aHedIgOIN1m37vPjFUF4EG+bF7/zKVpI6WYeawVXTje5iXarUPXPfeZOtZvDkSNbuHvZuBi8sqaJru8Ro5eiUaSXrojI7JOuN4Mncc71zPNWUxzsCZQNJnayLq4u1GbyZUl2pjOoyqt4ZzHzbtuHq2DEOUbU+T884rMjUBFuKK+nWPjGSdYcMDw6BIltU1pugfW+wQ4vEN6/EOwJlA0mbrGv8NZR7y7Wy3g9prjQqffteX+zNz8fd094VZ0sRj4eOkyaBy36TjvgDhq2l1XTLToxk7XQIHTNT2G6HyrpkY/zWsQ65YJp160qM35+KraRN1qEJUTRZN199ydpfVo5/1y7bj5JuSeJygc9Gk3kEFRRX4g8YurdPnC6JzlkpbCmpincYVmWd3Su+MXTqb93O+H1841C2kLTJenuldR1lp9ROcY4kcaW50qjy7f2H1bs5HwBPklTWAOIOVtU+G1xyFOb7/BIAhnbPjnMk0evbKYN1RXHuo91dYC2g0WlAfONIT46uJBWdpE3WRRVFAOSmx2nSgzagvsram28la3fPOFclrSnYBG5slqy/3bQLj8vBoG7t4h1K1A7MzWTTrgqq4jmL2fbV1m3ngfGLAcAd1iLis8HlbCqukjdZV1rJulOaVtbNle5K3ytZl86ZQ/7vrgbAbfOBVy1JXNY1zPZL1sUM6d4OjytxPub9OmdiDKwtKotfEAVLrNuO/eMXQ8iRU6zbFe/GNw4Vd4nzKW5h2yu3Iwgd07SpqbnSXGlU+6vxB6wqqPi/79Ruc7ZPnpXMxIaVdUWNj+/zSxjZK7HGZIzoZf2/WfhTHGfumnO7dZttgy+cvQ+zbv97eXzjUHGXtMm6qLKInNSc2pm4VNOluaxmuip/FSYQoHLJEjKPO44D3n8/ISbhaCmhPmvj3Xfq1Xj5YGkh1b4AJw6J4wxczdCrQzp5HdOZs3JbfAIIzVzWbUR8zl/XgBP33Pfac2U31TpimqxFZLyI/CAia0Tk5igk3EgAAAtYSURBVHq2p4jIm8HtX4tIXizjCVdYXqhN4Psp3Z0OQGlNKeVz5+LftYus448n5YC+cY6sdYUqa7sMMCuuqOH+WasY2DWLn+Ul3jwCp43owby12/luU3Hrn3z2VOt24gutf+6GdBlm3d7bNb5xqLiKWbIWESfwJDABGAycIyKD6+x2GbDLGNMPeAS4P1bxhNteuZ2FhQsZ1WVUa5yuzTHG4Nuxg6ErK7h8lp+CM37DpslX4urShcyjfhHv8FqdBGdrK/3kkzhHAmu2lTHp1cXsLK/hr2cNx+FIvBaOS47Io3t2Gle8soj3visgEGiFS+K8lfDFI7D039CxH3Q8MPbnjNbFM/bcXz07fnGouIrlTA5jgDXGmHUAIvIGcBqwImyf04A7gvffBp4QETExuGB1zoY5vPXDW2ws3UhBmdXUddqBp7X0adqkHf/4BxULF+EvKcFfUoK3sBBTUYEAJwis7rGdTr89k8Mv/xOOjORbbjTrmKNJP/wwtt59D2Wff07XW27B07v15gyvqPFx67RlrN5axtLNJaS4HDz06+EM7ZE4l2yFa5/u4cVLfsZFLyzgmteXcOPb3zO0RzvOGNmTc8b0ik0XS00FzLnDum+nqhogrb21+lf5NtjwJfQ/Id4RqTiIZbLuAWwKe5wPHNrQPsYYn4iUAB2BvVZcF5FJwCSA3s38I1jhq6C0ppThucM5o98ZjO0+liGdhjTrWMnGu6UQ77atOLOzSencn4wjj8DTsxcp/ftRPrAXjy++hysPPispEzWAIyOD3s8/z47nnmPXG2+2+s8h1eVkycZiumWnctP4gZw1umftPNuJakCXLOZMOYr7Z61iR3kNPxSWMn/dDs4ZE6NLAjM6wqTPrGurPemxOcf+uOFHa7nMDB0Qm6wkVrMuichZwInGmMuDjy8AxhhjrgnbZ3lwn/zg47XBfRpccX306NFm0aJFMYlZqf1ljEmqwXWtJRAw+AImoS5DsxsRWWyMGR3vOFTzxPJ/fj4Q/jW4J1DQ0D4i4gKyAV1tXSUsTdSx4XCIJmqV1GL5v38h0F9E+oqIBzgbmF5nn+nARcH7E4GPY9FfrZRSSiWymPVZB/ugrwY+BJzAC8aY5SJyF7DIGDMd+Afwqoiswaqoz45VPEoppVSiium6fsaY94H36zw3Nex+FXBWLGNQSimlEp12AimllFI2p8laKaWUsjlN1koppZTNabJWSimlbC5mk6LEiogUARticOhsoCQGx7WDRHtvdo033nHF4/ytdc5O1Jm5ULW4PsaY3HgHoZon4ZJ1rIjIs8aYSfGOIxYS7b3ZNd54xxWP87fWOUVkkc6upVTDtBl8j/fiHUAMJdp7s2u88Y4rHueP93tWSqGVtVLKBrSyVioyrayVUnbwbLwDUMrOtLJWSimlbE4ra6WUUsrmNFkrpZRSNqfJej+IyCAR+buIvC0iv413PMlMfxdKqbaszSRrEXGKyBIRmbEfx3hBRLaJyLJ6to0XkR9EZI2I3AxgjFlpjLkS+DWgI1kBEWkfTJirRGSliBzezOPo7yLJiUiGiCwWkZPjHYtS8dZmkjVwHbCyvg0i0llEsuo816+eXV8CxtfzeifwJDABGAycIyKDg9tOBb4A/rc/wbchjwGzjDEDgeHU+Z3o7yJ5NfQFrL4vX0E3AW+1bpRK2VObSNYi0hP4JfB8A7scBbwrIqnB/a8AHq+7kzHmc2BnPa8fA6wxxqwzxtQAbwCnBV8z3RgzFjhvv99IghORdsAvgH8AGGNqjDHFdXbT30Xyeok6X8Aa+vIlIscDK4CtrR2kUnbkincALeRR4EYgq76Nxph/i0hf4A0R+TdwKXBCE47fA9gU9jgfOFREjgbOBFKA95sRd1tzAFAEvCgiw4HFwHXGmPLQDvq7SF7GmM9FJK/O07VfvgBEJPTlKxPIwErglSLyvjEm0IrhKmUrCZ+sg/1Z24wxi4N/sOtljHkg+IfgaeBAY0xZU05T/yHNp8CnTThOW+cCDgGuMcZ8LSKPATcDt4XvpL8LFabeL1/GmKsBRORiYLsmapXs2kIz+BHAqSKyHqtJ9FgR+WfdnUTk58BQ4B3g9iaeIx/oFfa4J1DQrGjbtnwg3xjzdfDx21jJey/6u1Bh6v3yVXvHmJeMMc0eNKpUW5HwydoY83/GmJ7GmDzgbOBjY8z54fuIyEjgOazmtUuADiJyTxNOsxDoLyJ9RcQTPM/0FnkDbYgxphDYJCIHBZ86DqvfsZb+LlQd+uVLqSgkfLKOUjpwljFmbbA57SLqWRNbRF4HvgIOEpF8EbkMwBjjA64GPsQa3fyWMWZ5q0WfWK4B/iUi3wMj+P/27iXEqyqA4/j310RpKpmZ9KAaiiAVSdEgrYU9VpWLUiqoNi4MXERoQUJR9CLpBYWLXoKrDES0IEqRSppe5iMdeyyiICp6qNlUljX8Wpw7cfszOTNpzf1Pv8/q/s89555zZxa/Oeee/1x4oOV8fhdRlz++IgYh/xs8Iv4T1R9gc4GJlF3ed9l+VtLllE2iHcBK2/cP3ygjmilhHRER0XD/l2XwiIiItpWwjoiIaLiEdURERMMlrCMiIhouYR0REdFwCeuIiIiGS1hHREQ0XMI6IiKi4RLWESOEpCckbZN0/nCPJSKOrIR1xAggaQwwCbgJuHKYhxMRR1jCOtqOpMck3VL7/IqkZ2qfH5G05Aj3OZR3bg/meuMlLa597pTUPci2oyW9Lqmjr8z2T8AplHd6Py7pGEmbJbX9O+sjImEd7elNYA6ApKMoL4aYWjs/B+gahnENxXhg8YC1+rcQWGu7t69A0omUN5r1AL22DwKbgGsPd6ARMfwS1tGOuqjCmhLS3UCPpBMkHQtMBrZLWidpq6Tdkhb1NZa0vGVWe7ekpZJukPSupB2SnqzPXGt1+61TzYw/lPR01d8GSaOrc3dK+kjSRknPSboVeBA4u7rOQ9XlO/pr34/rgfUtZXcADwO7gSlV2bqqbkS0uYR1tB3bXwK/SzqDEtpvAe8As4FZwM5qZrnQ9syq7OZq9gmwmr/OOK8B3qvKLrQ9HeilJegkTR6gzjnACttTge+B+ZJmAfOBGcDV1VgAbgc+sT3d9m1/17713qt3Pp9l+7NaWWf1c3ie8o7vvlWGbiCbzSJGgDzPinbVN7ueAzwKnFYd76csk0MJ6Kuq49MpYbjH9nZJkySdCpwE7AOmATOBLZIARgPftPR56QB1PrW9ozreCnRSlujX2z4AIOnFQ9xTf+1bTaQEed19wD22LenPsLbdK+mgpHG2ew7Rb0Q0XMI62lXfc+tplBnk58BS4AdgpaS5wGXAbNs/S3oNGFVrvwZYAJxMmWkLWGV72SH6HKjOr7XjXkqYawj31F/7Vgeo3Yek6ZQZ+0WSVlTndtXqHwv8MoQxREQDZRk82lUX5StKe2332t5L2bQ1m7Isfjywrwrqc4ELWtqvBq6jBPYaymasBZImAUiaIOnMljaDqdPqDWCepFGSxgJXVOU9wLih3rTtfZRn232BvRyYZ7vTdidwHtXMulr2/9b2b0PtJyKaJWEd7WoXZUn47Zay/ba/A14Gjpa0E7i3pR62d1PC8gvbX9n+gLJJa0PVZiPlq1D1NgPWaWV7C/AC8D6wlvJsfL/tPUCXpO7aBrPB2kCZSV8CjLG9qdbf18AYSROAi4GXhnjtiGgg2R7uMUSMaJLG2v5R0nHAZmCR7W2Hcb0ZwBLbNw5Qby2wzPbH/7SviGiGPLOO+Pc9JWkK5XnyqsMJaoBqg9yrkjrq37Wuq3aNr0tQR4wMmVlHREQ0XJ5ZR0RENFzCOiIiouES1hEREQ2XsI6IiGi4hHVERETDJawjIiIaLmEdERHRcAnriIiIhvsD6F+yIicqwdMAAAAASUVORK5CYII=\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": [
      "ukidss_j: mean flux error: 6.34442663192749, 3sigma in AB mag (Aperture): 20.701215914476926\n",
      "90prime_g: mean flux error: 0.19539351761341095, 3sigma in AB mag (Aperture): 24.479921484593298\n",
      "90prime_r: mean flux error: 0.20366299152374268, 3sigma in AB mag (Aperture): 24.434916566562514\n",
      "mosaic_z: mean flux error: 0.8766658902168274, 3sigma in AB mag (Aperture): 22.850111590395464\n",
      "ukidss_j: mean flux error: 14.469086647033691, 3sigma in AB mag (Total): 19.806094069728168\n",
      "90prime_g: mean flux error: 8.714170455932617, 3sigma in AB mag (Total): 20.35663173615155\n",
      "90prime_r: mean flux error: 0.6011230945587158, 3sigma in AB mag (Total): 23.25978832955969\n",
      "mosaic_z: mean flux error: 8.997316360473633, 3sigma in AB mag (Total): 20.321914385057077\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_ukidss_j (11695.0, 13280.0, 1585.0)\n",
      "ap_90prime_g (4180.0, 5470.0, 1290.0)\n",
      "ap_90prime_r (5680.0, 7150.0, 1470.0)\n",
      "ap_mosaic_z (8552.0, 10018.0, 1466.0)\n",
      "ukidss_j (11695.0, 13280.0, 1585.0)\n",
      "90prime_g (4180.0, 5470.0, 1290.0)\n",
      "90prime_r (5680.0, 7150.0, 1470.0)\n",
      "mosaic_z (8552.0, 10018.0, 1466.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 XMM-13hr')"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfoAAAEgCAYAAABcsLEAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XlcVXX6B/DPc9kRRFEURQQ3VBQVMyy10bFSNDW3VrVGK7dp0tTMcsYmy8oxrZyytHV07JfllktaOuOWTiW4sooaZAKKIgiCyPL8/rgXuyGrclkOn/fr1at7z/me833upVfP+X7P955HVBVERERkTKbqDoCIiIhsh4meiIjIwJjoiYiIDIyJnoiIyMCY6ImIiAyMiZ6IiMjAmOip1hKReBG5p7rjMBJ+p0TGw0RPlcKSILJFJENE0kTkgIhMFpFK+W9MRD4TkVcr41y2IiJDRSRZRDyttt0vImdFxMPyXkXknIjYW7WxF5HzIqJW23Zb2nYt0sdGy/Z+pcSxQkRiRaRARP5UZN/Dln3plj7/JSL1b/3TE1FNxURPlWmoqroD8APwBoDnAXxcvSFVHVXdDOC/AN4CABFpAOB9AFNUNd2qaRqAQVbvBwO4VMwpTwB4rPCNiDQCcAeAlDJCOQpgKoBDxezbD6C3qnoAaA3AHsBNXUBZX6wQUc3FRE+VTlXTVXUTgIcAPC4inQFARJxE5E0R+cUyqv1ARFws+/qJyK8i8qKIXLDMEIyx7JsIYAyA2SKSKSKbrbrrJiLHLCPUNSLibDmmsYhsscwupIrIvpJmF0Skl4gctJzjoIj0stq3W0ReEZH9ltmK70SkcSkf/xkAg0RkIMwJf4/lu7C2ClYJ3PJ6ZTHnWg3gIRGxs7x/BMAGANdK6R+q+p6q/gfA1WL2nVHVC1ab8gG0LdKspO+08G/0vIgkA/i0tDiIqGZgoiebUdWfAPwK4C7LpoUAAgB0gzm5+ACYZ3WIN4DGlu2PA1ghIu1VdQXMSe8fquqmqkOtjnkQQCiAVgC6APiTZftMS99eAJoCeBHADc97tkyzbwWwFEAjAEsAbLWMngs9CmA8gCYAHAHMKuUzXwAwzRLvEJgTf1EbAfxBRBpYRv13Afi6mHaJAKIADLC8L+mCoEJEpI+IpAPIADAKwNtFmpT0nQLmv5EnzLM2E281FiKyPSZ6srVEAJ4iIgCeAvCsqqaqagaA1wA8XKT931Q1R1X3wJyAHyzj/EtVNVFVUwFshvkiAgByATQD4Kequaq6T4sv7HAfgDhVXaWqear6fwBiAFhfTHyqqidUNRvAl1Z9lOQHAB4AvlPV4qbZr1pifQjmz78JxYy+LVYCeExE2gNooKr/K6PvMqnq95ap+xYAFgGIL9KkpO8UAAoAvGT5G2XfaixEZHtM9GRrPgBSYR5ZuwIIt0ynpwHYbtle6JKqXrF6nwCgeRnnT7Z6nQXAzfJ6EYCTAL4TkdMiMqeE45tb+rGWYIm7rD5KsgLmBD3Y+jZAESthHqGXNUpfD6A/gL/APOX/O5ZbGYX/tCwjrt9R1bMw/w2+KLKrtM+boqolXZQQUQ3ExTRkMyJyO8wJ83sAFwBkA+hkSTDFaSgi9aySfUsAEZbXFSqzaJkxmAlgpoh0ArBLRA5a7l1bS4R5GtpaS5gTYIWJyBMAfGGetg8D8KGIBKtq0fvq+2CecVCYv582JXyOLBHZBmBKcW1UtayLjrLYl9R3CVjukqiW4YieKp2I1BeRITCPFP+tqsdVtQDAhwDeEpEmlnY+lkVr1l4WEUcRuQvmZPmVZfs5mFeJlzeGISLS1nLL4DLMi87yi2n6DYAAEXnU8jO3hwAEAthS/k98vc/mMM8kPKWqOQA+AHARwNyibS23EYYCGFbCLQVrLwLoq6rx5YzD0bKATgA4iIhz4UJEERkjIi3FzA/AAgBFL36IyECY6KkybRaRDABnYE5uS2BexFboeZin038QkcsAdgJob7U/GeafmSXCvJhtsqrGWPZ9DCDQMu2/sRyxtLOcPxPA/wAsU9XdRRup6kWYLyhmwpyUZwMYUmRlenktA/CFqu6znFthXpcw3TKrULTvSFWNLOuklvvl31cgju9gnj3pBfNthGwAf7DsCwRwAObvZT+AWEuMRGRQUvZggsj2LA+A+beqtqjuWIiIjIQjeiIiIgNjoiciIjIwTt0TEREZGEf0REREBsZET0REZGC14oE5jRs3Vn9//+oOg4ioVgkPD7+gql5ltyzx+Cb29vYfAegMDgxrqgIAEXl5eU/edttt54trUCsSvb+/P8LCwqo7DCKiWkVEij7euULs7e0/8vb27ujl5XXJZDJxQVcNVFBQICkpKYHJyckfARhWXBteoRERUUk6e3l5XWaSr7lMJpN6eXmlwzzrUnwbW3UuIr4isktEokUkUkSmFdk/S0S0jNreRERUfUxM8jWf5W9UYj635dR9HoCZqnpIRNxhrlq2Q1WjRMQXwL0AfrFh/0RERHWezUb0qpqkqocsrzMAROO30p9vwfxMcV4pEhFRtVi6dGmjxx577Ibyzv/4xz+83n333UZFt8fGxjq2a9fuhroVN+Ohhx7yCw8Pd66Mc5WlShbjiYg/gGAAP4rIMABnVfWoubBYicdMBDARAFq2rFCZbSIiops2e/bsFFv3sWbNmltaKFkRNk/0IuIGYB2A6TBP588FMKCs41R1BcyVt9CjRw+O/ImIqtFza4/6nkjOcK3McwZ4u2ctGt31TFnt7rnnnjZJSUmOOTk5psmTJ5+bNWvWBVdX1+AxY8ak7N+/393DwyN/3bp1p5s3b55X3PEhISHt33zzzTN/+MMfspKSkux79OjR8ezZs8et23zxxRceb7zxRrNt27adXLRoURM3N7f8+fPnn9u3b5/rk08+6e/i4lLQs2fPzML2YWFhzuPHj2+Vm5srBQUFWLdu3Sk/P7/cYcOGtU5KSnIsKCiQ2bNnJz711FOXyoqpot9bRdk00YuIA8xJfrWqrheRIACtABSO5lsAOCQiIaqaXOkBbJsDJB8vu11d5R0EDHqjuqMgIirV6tWr45s2bZqfmZkpwcHBgWPHjr2UnZ1t6t69e9aHH37466xZs5rNmTOn+cqVK29q3dfKlSsbvPPOO0137NgR5+XllW+974knnvB/6623frnvvvsyJ02adL265j//+U+vqVOnnpsyZUrq1atXJS8vD2vXrvXw9vbO3b1790kAuHjxot2tffLKYbNEL+ZM/jGAaFVdAgCqehxAE6s28QB63GTtbyIiqiLlGXnbysKFC5tu3bq1AQAkJyc7REZGOptMJjz55JOpADBhwoSLI0eObHsz5z5w4ID70aNHXXft2nXC09OzwHrfxYsX7TIyMuzuu+++zMJ+/vvf/3oAwJ133nnlzTffbPbrr786Pvzww5eCgoJyunfvnj137lzfKVOm+Nx///3poaGhmcX1WdVsOaLvDWAcgOMicsSy7UVV/caGff4eR6u3bOFPCxGTGlPdYVS5Dp4d8HzI89UdBlGdt2XLFvc9e/a4h4WFxbi7uxeEhIS0z87OvmEheWlrvuzt7TU/3zxQz8rK+l3Dli1b5vzyyy9OERERzkWn0VW1xPNOnjw59a677rqyYcMGj0GDBgUsW7YsftiwYRmHDh2KWrduncfcuXN9du7cefnNN99MuomPXalsuer+e1UVVe2iqt0s/3xTpI0/R/NERFSStLQ0Ow8Pj3x3d/eCw4cPOx89erQeABQUFODTTz9tCACfffZZo5CQkIySzuHr65vz008/1QOA1atXN7Te16JFi2vr1q07OX78+FZhYWG/WwXfuHHjfDc3t/xvv/3WzdKPZ+G+qKgox44dO+b89a9/PT9gwIC0I0eOuMTHxzu4u7sXTJ06NXX69Onnjhw5UqlrGm5WrXgELlUfjmqJqDqNGjUqfcWKFV4BAQGBbdq0udq1a9crAODi4lIQGRnp0qlTJ293d/f89evXny7pHHPmzDn30EMPtf7iiy8a3XXXXZeL7u/atWvOypUrTz/00ENtNm3adNJ638cffxxfuBivf//+149dtWqV51dffdXI3t5evby8cl9//fXE77//vt4LL7zQwmQywd7eXpctW1bqyvrSZiEqU62oR9+jRw/ls+6JiCpGRMJVtcfNHn/06NH4rl271shZV1dX1+CsrKzD1R3HzQoICAjctGnTyQ4dOlyrjPMdPXq0cdeuXf2L28dn3RMREVWhXr16tWvfvn12ZSX5snDqnoiIap3iRvPjxo1refDgQTfrbVOmTDk3bdq0i1UX2e/de++9bc6cOeNkvW3BggW/jho16oZbCLbCRE9ERIawatWqGlc/ZceOHaeqOwZO3RMRERkYEz0REZGBMdETEREZGBM9ERGRgTHRExGRof3vf/9z6datW4eAgIDA/v37t01NTb2e+1544QXvli1bdvb39++8bt26+hU99/Tp05tv3LjRvXIjrlxM9EREZGhPPfWU/4IFC349ceJE1LBhwy69/PLL3gAQHh7uvH79es/Y2NjI7du3n5g+fXrLvLxiK90WKy8vD2+//Xbi8OHDS3z8bk3An9cREVHZNv7ZF+ejKvfZ7U0CszD8PZvXo4+Pj3ceNGhQJgAMGTLk8sCBAwPeeeedxLVr1zYYOXJkqouLi3bo0OGan59fzu7du+v5+vrmhoaGtgsODr4SERHh2rp166tfffVVvLu7e4GPj0/QI488cmHXrl31J02adP7bb7/1GDJkSPr48eMv+fj4BI0YMSL1+++/d8/Ly5MPPvggYc6cOT4JCQlOf/nLX87Nnj07BQD+9re/Nd2wYYPntWvX5L777kt76623Ekv67M8991yztWvXejZr1uxao0aN8oKDg7Pmz59/riJfM0f0RERUo61evTo+MjIy+siRI1HLly9vmpycbFdYjz4qKiq6d+/eGXPmzGle0vHt2rXL/vzzzxsAwL///W/P5ORkRwA4e/aso6+v7/Wn0zVv3vzamTNnHAHzxcHkyZNTTpw4EeXu7l6waNEir8J2zs7OBeHh4bETJ068VLQvX1/fa0eOHInp2bNn5oQJE/w3b9586scff4x54403mgPA+vXr6588edL52LFj0dHR0VFHjhxx3bZtm1vR8wDA3r17XTdv3tzw+PHjUVu3bj117Nixejfz/XFET0REZSvHyNtWbrUe/SeffBL/9NNP+77++uvNQkND0xwcHBQwl6EtSkQUALy9va8NGDDgCgCMGzfu4tKlS5sAOAcAjz322A0JvtCDDz6YBgBBQUFZV65cMTVs2LCgYcOGBU5OTgUXLlyw2759e/29e/fWDwwMDASArKwsU0xMzPUZB2u7d+92GzRoUJqbm5sC0HvvvTetnF/Z7zDRExFRjVUZ9eiDg4Ov7t+/Pw4Ajh075vTdd981AMwlagtH8ACQmJjo2KJFi9zizmf93t3dvaCkvpydnRUATCYTHB0dr19JmEwm5Obmiqpi+vTpSc8991yZxYIqq+gcp+6JiKjGqox69GfPnrUHgPz8fLz00kvNnnjiifMAMGrUqLT169d7ZmdnS0xMjGN8fLxzv379rgBAUlKS486dO+sBwOeff+7Zq1evG0bcN2PQoEGXV61a1Tg9Pd0EAD///LNDYXxF9evXL/Pbb7/1yMrKkvT0dNPOnTsb3EyfHNETEVGNVRn16D/55BPPjz/+uAkADB48+NIzzzxzEQB69Ohxdfjw4akBAQGd7OzssGTJkgR7e3NabN269dVPPvmk0dSpU/1atWqVM2vWrJTK+DwjR468HBkZ6Xz77bd3AABXV9eC1atX/+zj43PDQsK+fftmhYaGpgcGBnby8fHJ6dKlyxUPD4/8ivbJevRERAbFevQ3JzY21nHIkCHt4uLiIm1x/opIT083eXh4FGRkZJjuvPPO9h988EFCnz59soq2K60ePUf0RERENdTYsWP94uLiXHJycuThhx++WFySLwsTPRER1Tq2rEffvn37a1U5mk9OTrbr169f+6Lbd+/eHbt58+afb/X8TPRERGQINbEefXl4e3vnx8TERNnq/Fx1T0REZGBM9ERERAbGRE9ERGRgTPREREQGxkRPRESGxnr0REREBmaEevQFBQXIz6/wQ/EA8Od1RERUDn/b/zffk5dOVmo9+rYN22a90vsV1qMvoR59bGys46BBg9r16tUrIzw83O3rr78+GRAQcK24tqXhiJ6IiGq0ulqPvjCO8ePHX4yOjo66mSQPcERPRETlUJ6Rt63U1Xr0ANCsWbNrd99995XyfE8lYaInIqIaqy7XowfM1e3K0640nLonIqIaqy7Xo68sHNETEVGNVZfr0VcWm9WjFxFfACsBeAMoALBCVd8RkVcA3G/Zdh7An1S12BWHhViPnoio4liP/ubUpHr05VVaPXpbTt3nAZipqh0B3AHgzyISCGCRqnZR1W4AtgCYZ8MYiIiI6jSbTd2rahKAJMvrDBGJBuCjqtal+OoBsM2UAhERGVZdqUfv7e19c0/JsVIl9+hFxB9AMIAfLe8XAHgMQDqAP1ZFDEREZGysR188m6+6FxE3AOsATFfVywCgqnNV1RfAagBPl3DcRBEJE5GwlJRKWQNBRERU59g00YuIA8xJfrWqri+myecARhV3rKquUNUeqtrDy8uruCZERERUBpslejE/XeBjANGqusRqezurZsMAxNgqBiIiorrOlvfoewMYB+C4iByxbHsRwBMi0h7mn9clAJhswxiIiIjqNFuuuv8eQHHPJPzGVn0SERFVp+nTpzfv169fRlWUri0vPhmPiIiokrz99tulPgCuOjDRExFRmRJfnOubExdXqfXondq1y2r+2gKb16MPCQlpHxQUlHX06FHX1NRU+08//fTnBQsWNIuNjXW5//77U5cuXZoIAH//+9+brl69ujEAjBs3LmXevHnnL1++bBo2bFjrpKQkx4KCApk9e3biU089dWnWrFnNtm/f3iAnJ8fUo0ePzNWrVyeYTCaMGjXKv7A+/Z49e1ynT5/eMisry+To6Kh79+6Nbdiw4Q1Fah566CG/wmf4nzt3zmHChAnnFy9enHRr3+5vWNSGiIhqtFutRw8Ajo6OBWFhYbHjx49PeeCBB9p++OGHv8TExESuWbOmcXJyst2+fftcP//880bh4eHRYWFh0StXrvTav3+/y/r16+t7e3vnxsbGRsXFxUWOHDnyMgA899xz5yMiIqLj4uIis7OzTV988YWHdX9Xr16VMWPGtHn77bd/iY2NjdqzZ0+sm5tbsZXo1qxZkxATExO1adOmkw0aNMibNGlShR7wUxaO6ImIqEzlGXnbyq3WoweAESNGpAFA165ds9u2bZvt5+eXCwC+vr45p0+fdty9e7fb4MGD0+rXr18AAPfdd9+lXbt2uQ8bNix97ty5vlOmTPG5//7700NDQzMBYNu2be5Llizxvnr1qiktLc0+MDAwG+aHwAEAjh075tykSZPcvn37ZgGAp6dnqeVms7KyZNSoUW3eeuutXwICAq7d9JdVDI7oiYioxrKuRx8bGxvVsWPH7IrWowd+Xyfeycnpd3Xi8/LySizw1qVLl5xDhw5FBQUFZc+dO9dn1qxZzbKysmTmzJl+69evP3XixImosWPHXrh69ervYlJViEi5H/E+btw4v6FDh16yxSI+JnoiIqqxKqMefXn0798/85tvvmmQkZFhunz5sumbb75p+Mc//jEjPj7ewd3dvWDq1Kmp06dPP3fkyBHXrKwsEwB4e3vnpaenmzZv3tyw6Pm6du169dy5c4579uxxBYBLly6ZcnNzi+379ddf98rMzLR77bXXkm/lM5SEU/dUqrTNp3At8Up1h1HlHJvXQ4Ohbao7DKI6rzLq0ZdHnz59sh599NGL3bt37wiYF+P17t07e926dfVfeOGFFiaTCfb29rps2bKExo0b548ZMyYlMDCwU4sWLa4VxmTN2dlZV69efeqZZ55pefXqVZOzs3PB3r17T3h4eNwwhf/uu+96Ozg4aIcOHQIBYMKECSmzZ8+utGe/26wefWViPfrqw0RPVHuxHn3dUVo9eo7oqVRMdkREtRsTPRER1Tq2rEdvK+vWras/d+7cFtbbfH19c3bs2HHKlv0y0RMRkSHU9Hr0o0aNujxq1Cib1Z0vCVfdExERGRgTPRERkYEx0RMRERkYEz0REZGBMdETEVGtExIS0n7v3r03VNPr27dv2wsXLtgV3T5jxozm8+bNa3qr/cbHxzuEhoa2vtXzVCWuuiciIsPYs2fPSVue39/fP3f79u239BS+qsZET0REZfrPymjf1LOZlVqP3tPHLevuxzqWWhUvNjbWcciQIe3i4uIiAWDevHlNMzMzr4/Y8/Pz8cADD/i3aNHi2tKlSxN9fHyCwsLCops1a5b3/PPPe69Zs6Zx8+bNrzVq1Cg3ODg4CwBeffXVJp9++qmXnZ2dBgQEXN2yZcvprVu3us2cObMlYC6Qc+DAgZjiascXjac2YKInIqJaKTc3V4YPH94qMDAwe+HChb8rCLNv3z7XDRs2eB4/fjwqNzcX3bp1CyxM9EuXLvVOSEg47uLiooXT/IsXL/ZeunRpwoABA66kp6ebXF1dSy0rW5sw0RMRUZnKGnlXh6lTp/oNHz48tWiSB4Bdu3a5DR48OM3d3b0AAAYMGJBWuK99+/bZI0aMaDVs2LC0MWPGpAHAHXfckTlr1izfBx98MPWRRx651KZNG8Mkei7GIyKiGsve3l4LCn7LudZ133v06JG5b9+++llZWcUWoy+pRv2uXbvi/vznP6eEh4fX69q1a2Bubi5ee+215I8++ighOzvb1KtXr46HDx92ruzPUl2Y6ImIqMZq0aJFXmpqqn1ycrJddna2fPvttx6F+yZNmnRhwIAB6UOGDGlTtNZ7//79M7du3dogMzNTLl26ZNqxY0cDwHxP/9SpU45Dhw7NWLZs2a8ZGRl26enpdpGRkU4hISHZCxYsSA4KCroSERFhmETPqXsiIqqxnJycdObMmUkhISEdW7RokdO2bdur1vv//ve/n3v22WftRo4c2Wrjxo0/F27v06dP1ogRI1I7d+7cycfHJyckJCQTAPLy8uTRRx9tlZGRYaeqMmnSpHONGzfOnzlzZvMDBw7UN5lMGhAQkD169Oj0kmISkZpf390K69ETGUjya68hJzqmusMok1PHDvB+8cXqDsPwjFyPvrrs27fPdcaMGb4HDx6Mre5YrJVWj55T90REROWwd+9e13HjxrV++umnz1V3LBXBqXsiA+Eomahy/PTTTy6PPfZYK+ttjo6OBfHx8RHVFdPNYqInIiIqIiQkJDsmJqbKa8fbAqfuiYiIDIyJnoiIyMCY6ImIiAyMiZ6IiMjAmOiJiKhGe+WVV5q0a9euU9u2bTvNnz+/CQCcO3fOrlevXu38/Pw69+rVq11KSsoNNejLEhwc3KHyo615mOiJiKjGOnjwoPPKlSu9Dh06FB0dHR25ffv2BsePH3d66aWXmvXr1y8jISEhol+/fhnz5s3zLu858/LyAACHDx+u+U+XqgT8eR0REZXp2/ff9r1wJqFS69E39vXLGjhleqlV8Y4fP+7SvXv3zMIqdL17985Ys2ZNg+3btzfYs2dPLABMmjTpYt++fdsDODtjxozmp0+fdkpOTnZISkpyfOaZZ5Jnzpx5YcuWLe6vvPJKsyZNmuRGRUW5njp1KtLV1TU4Kyvr8JYtW9xffvnl5l5eXrlRUVGugwcPvhQUFJS9bNmypjk5ObJhw4ZTnTp1yklMTLQfP36839mzZx0BYMmSJb8MGDDgSnFxJyYm2o8ePbpVWlqafbdu3bJ2795dPzw8PLpZs2Z5lfkdlofNRvQi4isiu0QkWkQiRWSaZfsiEYkRkWMiskFEGtgqBiIiqt26deuW/eOPP7onJyfbZWRkmHbs2OFx5swZx4sXL9r7+fnlAoCfn19uamrq9YFrdHS0y86dO+N++OGHmEWLFjWPj493AIBjx47VW7Ro0dlTp05FFu0nJibG5f333z8THR0duXbt2kYnTpxwPn78ePS4ceMuLF68uAkATJo0yXfGjBnnIiIiojds2HBq8uTJ/iXFPWfOnOZ9+/bNiIqKih45cuSlpKQkx0r/csrJliP6PAAzVfWQiLgDCBeRHQB2AHhBVfNEZCGAFwA8b8M4iIjoFpU18raV7t27X502bVpy//79A1xdXQsCAwOz7O1LT12DBg1Kc3NzUzc3t7w777zz8r59++o1bNgwv0uXLlc6dOhwrbhjgoKCrhReOLRs2TJn0KBB6QDQtWvX7D179rgDwP79++vHxcW5FB6TmZlpd+nSJVPDhg1vqF3/008/uW3cuPEkAIwePfpy/fr182/6S7hFNkv0qpoEIMnyOkNEogH4qOp3Vs1+ADDaVjEQEVHt9+yzz1549tlnLwDA008/7dOiRYtrjRo1yktISHDw8/PLTUhIcPD09Lw+JV60Dn3he1dX1xsSciEnJ6frFd5MJhOcnZ218HV+fr4AgKoiLCws2s3NrcxqcDWpYFyVLMYTEX8AwQB+LLJrAoBtJRwzUUTCRCQsJSXFtgESEVGNdfbsWXsAiIuLc9y6dWuDJ554InXgwIFpy5cvbwQAy5cvbxQaGppW2H7btm0NsrKyJDk52e6HH35w79OnT7H30SuqT58+lxcuXNik8P2BAwdcSmobEhKSuWrVKk8AWL9+ff3Lly9X+FcBlcXmi/FExA3AOgDTVfWy1fa5ME/vry7uOFVdAWAFYC5Ta+s4iYioZho2bFibtLQ0e3t7e3377bd/8fLyyn/55ZeTRowY0cbPz69x8+bNr23cuPFUYfvg4OArd999d7vExETHWbNmJfn7++dGREQ432ocK1asOPPkk0+2DAgICMzPz5eePXtm9OrV65fi2r7xxhuJo0ePbh0YGNjwzjvvzPTy8spt0KBBtUzf27QevYg4ANgC4FtVXWK1/XEAkwHcrapZZZ2H9eiJiCquLtajnzFjRnM3N7f8+fPnV2sp2ezsbLG3t1cHBwfs3Lmz3tNPP+1nyyI5pdWjt9mIXsw3RT4GEF0kyYfCvPiub3mSPBH8yfUsAAAUkklEQVQRUW1z8uRJxwcffLBNQUEBHBwcdPny5fHVFYstp+57AxgH4LiIHLFsexHAUgBOAHZYFkj8oKqTbRgHERHVEUuWLEmsyv7eeeedRu+//35T622333575qpVq36Jjo6uEWVuy5XoRcQJwCgA/tbHqOr8ko5R1e8BSDG7vqlYiERERDXTtGnTLk6bNu1idcdRmvKO6L8GkA4gHECO7cIhIiKiylTeRN9CVUNtGgkRERFVuvL+jv6AiATZNBIiIiKqdKWO6EXkOAC1tBsvIqdhnroXAKqqXWwfIhEREd2sskb0QwAMBTAIQFsAAyzvC7cTERHZFOvR35pSR/SqmgAAIrJKVcdZ7xORVTD/fI6IiAwude0J39zkK5VaptbBu16W5+iAUovlWNejd3Z2Lujbt2/AiBEj0t977z2vfv36Zbz22mtxL774ove8efO833///bPl6TcvLw/29vZVVo8+NzcXDg4OVdFVscp7j76T9RsRsQNwW+WHQ0RE9BvrevQODg6/q0c/adKki4C5Hv22bdsaAuYn4w0fPrzVHXfcEeDn59d58eLFjQFgy5Yt7j179gwYOnRoq/bt23cCAFdX1+DCfbfffnv7wYMHt/b39+88depUn/fff98zKCioY0BAQGBkZKQTYK4xP3DgwDadO3fu2Llz547fffddvZLinjFjRvNHHnnEr3fv3u1GjhzZytbfU2nKukf/AswPuXERkcv47Xfx12B5Dj0RERlfWSNvW+nWrVv2/PnzfZKTk+3q1aunO3bs8OjateuVsurRh4eHR2dkZNgFBwcHjho1Kh0w16M/fPhwZHGlamNiYlzWrl17ukmTJnl+fn5BTk5OF44fPx79yiuvNFm8eHGTTz755ExhPfqBAwdmxsXFOQ4cOLDd6dOnb6htX+jYsWOuP/74Y0x5qt3ZUllT968DeF1EXlfVF6ooJiIiIgC1tx49AISGhqZVd5IHyv87+hdFZCSAPjCvwt+nqhttFxYREZFZbaxHDwD16tUrsb+qVN579O/BXG3uOIAIAJNF5D2bRUVERGRRG+vR1yTlHdH3BdBZLTVtReRfMCd9IiIim6qN9ehrknLVoxeR9QCetfq5nR+AN1T1ERvHB4D16ImIbgbr0dcdlVGPvhGAaBH5yfL+dgD/E5FNAKCqw245SiIiIqp05U3082waBRERUSWoSfXoqzKO0pQr0avqHst0fTtV3SkiLgDsVTXDtuERERHVXLWhHn25Vt2LyFMA1gJYbtnUAgB/XkdERFTDlffndX8G0BvAZQBQ1TgATUo9goiIiKpdeRN9jqpef5qQiNjD/OAcIiIiqsHKm+j3iEjhM+/vBfAVgM22C4uIiKhqGbVsbXkT/RwAKTA/JGcSgG8A/NVWQREREVW1qipbW9XKu+q+QEQ2Atioqik2jomIiGqYjRs3+p4/f75S69E3adIka/jw4aVWxYuNjXUMDQ1tFxISknno0CG3jh07Zk2YMOHC/PnzfS5evGj/2WefnQ4MDMwZM2aM/y+//OLk4uJSsGLFioSePXtmb9261W3mzJktAfPz7g8cOBBjMpkQGhraNj093S4vL0/mzZuXOHbs2DTAXLY2KyvrMAD89a9/bfrll182EhHcfffd6cuWLbuh1n18fLxDaGhou8L3cXFxLtHR0ccDAgKKLZxTXcoqUysAXgLwNMwlakVE8gH8U1XnV0F8RERUx505c8Z5zZo1p2+77baELl26dFy9enWjsLCwmM8//7zBggULmvn4+Fzr2rVr1s6dO09t2rTJ/fHHH28VExMTtXjxYu+lS5cmDBgw4Ep6erqpsKjN1q1bT3p6ehYkJSXZ9+zZs8Ojjz6aZjL9NsH95Zdf1t+6dWvD8PDwGHd394Jz587ZFReXv79/bkxMTBQAvP7661779u1zr2lJHih7RD8d5tX2t6vqzwAgIq0BvC8iz6rqW7YOkIiIql9ZI29b8vHxyQkJCckGgICAgOz+/ftfNplM6N69e9arr77a/OzZs07r1q07CQDDhg3LmDhxov3Fixft7rjjjsxZs2b5Pvjgg6mPPPLIpTZt2hTk5OTI9OnTW/zwww9uJpMJ58+fd/z111/tW7Zseb363Y4dO+qPHTv2gru7ewEANG3aNL+0+L777rt6K1eu9Prhhx9q5NR/WffoHwPwSGGSBwBVPQ1grGUfERGRTTk6OhZbQtbOzg75+fnF1mwREX3ttdeSP/roo4Ts7GxTr169Oh4+fNh5+fLlnhcvXrQ/fvx4dExMTFSjRo1ys7Ozf5cLVfWGUrclSUhIcJg0aZL/mjVrTnl4eNSIsrRFlZXoHVT1hoIGlvv0DrYJiYiIqPzuuOOOjE8//bQRAGzZssW9YcOGeZ6engWRkZFOISEh2QsWLEgOCgq6EhER4Zyenm7XuHHjXCcnJ928ebN7YmKiY9HzhYaGXl61alXjjIwMEwCUNHWfk5MjI0eObP3KK6+c7dKlS45tP+XNK2vqvrR7DTXuPgQREdU9CxcuTHz00Uf9AwICAl1cXAo+++yznwHgH//4R5MDBw7UN5lMGhAQkD169Oj0tLQ0u0GDBrXt3Llzx06dOmW1atXqatHzjR49+vKhQ4dcu3Xr1tHBwUHvueee9HffffeGxXg7d+6sFxERUe/VV19t/uqrrzYHgO3bt8f5+/vn2v5Tl1+pZWotC++uFLcLgLOqVsmonmVqiYgqri6Wqa2rbrpMraoWO11BREREtUN5y9QSERHVWePGjWt58OBBN+ttU6ZMOVfTK9cBTPRERERlqkn15SuqvI/AJSKiuqegoKCgfL8zo2pj+RuV+NM+JnoiIipJREpKigeTfc1VUFAgKSkpHgAiSmpjs6l7EfEFsBKAN8xXGitU9R0ReQDA3wF0BBCiqlxOT0RUA+Xl5T2ZnJz8UXJycmdwYFhTFQCIyMvLe7KkBra8R58HYKaqHhIRdwDhIrID5quOkQCW27BvIiK6Rbfddtt5AMOqOw66NTZL9KqaBCDJ8jpDRKIB+KjqDgDlfrwgERER3bwqmYoREX8AwQB+rMAxE0UkTETCUlJYGZeIiOhm2DzRi4gbgHUApqvq5fIep6orVLWHqvbw8vKyXYBEREQGZtNELyIOMCf51aq63pZ9ERER0Y1slujFfBP+YwDRqrrEVv0QERFRyWy56r43gHEAjovIEcu2FwE4AfgnAC8AW0XkiKoOtGEcREREdZYtV91/D3OVu+JssFW/RERE9Bs+AIGIiMjAmOiJiIgMjImeiIjIwJjoiYiIDIyJnoiIyMCY6ImIiAyMiZ6IiMjAmOiJiIgMjImeiIjIwJjoiYiIDIyJnoiIyMCY6ImIiAyMiZ6IiMjAmOiJiIgMjImeiIjIwJjoiYiIDIyJnoiIyMCY6ImIiAyMiZ6IiMjAmOiJiIgMjImeiIjIwJjoiYiIDIyJnoiIyMCY6ImIiAyMiZ6IiMjAmOiJiIgMjImeiIjIwJjoiYiIDIyJnoiIyMDsqzsAIqKa5uXNkYhKvFzdYQAAApvXx0tDO1V3GFSLcURPRERkYBzRExEVwRE0GQlH9ERERAZms0QvIr4isktEokUkUkSmWbZ7isgOEYmz/LuhrWIgIiKq62w5os8DMFNVOwK4A8CfRSQQwBwA/1HVdgD+Y3lPRERENmCzRK+qSap6yPI6A0A0AB8A9wP4l6XZvwAMt1UMREREdV2V3KMXEX8AwQB+BNBUVZMA88UAgCYlHDNRRMJEJCwlJaUqwiQiIjIcmyd6EXEDsA7AdFUt9w9TVXWFqvZQ1R5eXl62C5CIiMjAbJroRcQB5iS/WlXXWzafE5Fmlv3NAJy3ZQxERER1mS1X3QuAjwFEq+oSq12bADxuef04gK9tFQMREVFdZ8sH5vQGMA7AcRE5Ytn2IoA3AHwpIk8A+AXAAzaMgYiIqE6zWaJX1e8BSAm777ZVv0RERPQbPhmPiIjIwJjoiYiIDMzQRW12fbYC5xNOV3cYNVYTv9b4458mVncYRERkQ4ZO9ER1zbZt25CcnFzdYZTJ29sbgwYNqu4wiOoEQyd6jlaJiKiuM3SiJ6prOEomoqK4GI+IiMjAmOiJiIgMjImeiIjIwJjoiYiIDIyJnoiIyMCY6ImIiAyMiZ6IiMjAmOiJiIgMjImeiIjIwJjoiYiIDIyJnoiIyMCY6ImIiAyMiZ6IiMjAmOiJiIgMjImeiIjIwJjoiYiIDIyJnoiIyMCY6ImIiAyMiZ6IiMjA7Ks7ACKimmbflydw4UxmdYcBAGjs64a7Hgyo7jCoFuOInoiIyMA4oiciKoIjaDISjuiJiIgMjImeiIjIwJjoiYiIDIyJnoiIyMCY6ImIiAyMiZ6IiMjAmOiJiIgMjImeiIjIwERVqzuGMolICoAEq00eANKrKZyqVNs/Z22IvybFWJ2xVHXfjQFcqML+6io/VfWq7iCoetWKRF+UiKxQ1YnVHYet1fbPWRvir0kxVmcsVd23iISpao+q6o+oLqutU/ebqzuAKlLbP2dtiL8mxVidsdSk74GIKlGtHNETUe3GET1R1amtI3oiqt1WVHcARHUFR/REREQGxhE9ERGRgTHRExERGVidSfQi0lFEPhCRtSIypbrjoZLxb0VEVHmqJdGLiJ2IHBaRLbdwjk9E5LyIRBSzL1REYkXkpIjMAQBVjVbVyQAeBMDVvuUgIg0syTZGRKJF5M6bPA//VlQqEaknIuEiMqS6YyEymuoa0U8DEF3cDhFpIiLuRba1LabpZwBCizneDsB7AAYBCATwiIgEWvYNA/A9gP/cSvB1yDsAtqtqBwBdUeRvxr8VlaSki7viLuwsngfwZdVGSVQ3VHmiF5EWAO4D8FEJTfoC+FpEnC3tnwKwtGgjVd0LILWY40MAnFTV06p6DcAXAO63HLNJVXsBGHPLH8TgRKQ+gD8A+BgAVPWaqqYVaca/FZXkMxS5uCvpwk5E7gEQBeBcVQdJVBfYV0OfbwOYDcC9uJ2q+pWItALwhYh8BWACgHsrcH4fAGes3v8KoKeI9AMwEoATgG9uIu66pjWAFACfikhXAOEApqnqlcIG/FtRSVR1r4j4F9l8/cIOAESk8MLODUA9mJN/toh8o6oFVRgukaFVaaK33H87r6rhlv+ZF0tV/2H5n8D7ANqoamZFuin+lLobwO4KnKeuswfQHcBfVPVHEXkHwBwAf7NuxL8VVUCxF3aq+jQAiMifAFxgkieqXFU9dd8bwDARiYd5mra/iPy7aCMRuQtAZwAbALxUwT5+BeBr9b4FgMSbirZu+xXAr6r6o+X9WpgT/+/wb0UVUOyF3fUXqp+p6k0v0CWi4lVpolfVF1S1har6A3gYwH9Vdax1GxEJBvAhzFN64wF4isirFejmIIB2ItJKRBwt/WyqlA9Qh6hqMoAzItLesulumO+jXse/FVUQL+yIqkFN/B29K4AHVPWUZQrvcfy+Fj0AQET+D8D/ALQXkV9F5AkAUNU8AE8D+BbmVeJfqmpklUVvLH8BsFpEjgHoBuC1Ivv5t6KK4IUdUTXgs+6JqNJZLu76AWgM82r6l1T1YxEZDPOCXDsAn6jqguqLkqhuYKInIiIysJo4dU9ERESVhImeiIjIwJjoiYiIDIyJnoiIyMCY6ImIiAyMiZ6IiMjAmOiJiIgMjImeiIjIwJjoiWoxEfmniBwSkdurOxYiqpmY6IlqKRGpB6AJgEkAhlRzOERUQzHRU60gIm+JyHSr99+KyEdW7xeLyIxK7jOzks/XQESmWr33F5GIch7rIiJ7RMSucJuqXgHQDMBuAEtFxFFE9oqIfWXGTUS1GxM91RYHAPQCABExwVwspZPV/l4A9ldDXBXRAMDUMlsVbwKA9aqaX7hBRBrBXEEwA0C+ql4D8B8AD91qoERkHEz0VFvshyXRw5zgIwBkiEhDEXEC0BHAYRHZKCLhIhIpIhMLDxaRhUVG038XkZkiMlZEfhKRIyKy3HrEbNW22DaWEXm0iHxo6e87EXGx7PubiMSIyA4R+T8RmQXgDQBtLOdZZDm9XXHHF2MMgK+LbPsrgDcBRAIItGzbaGlLRASAiZ5qCVVNBJAnIi1hTvj/A/AjgDsB9ABwzDKinaCqt1m2PWMZ9QLAF/j9SPdBAGGWbb1VtRuAfBRJkiLSsYw27QC8p6qdAKQBGCUiPQCMAhAMYKQlFgCYA+CUqnZT1edKOr7oZ7fUbm+tqvFW2/wt38MaANH4bXYjAgAX5hHRdbyXR7VJ4ai+F4AlAHwsr9NhntoHzMl9hOW1L8yJ9KKqHhaRJiLSHIAXgEsAggDcBuCgiACAC4DzRfq8u4w2P6vqEcvrcAD+MN9W+FpVswFARDaX8pmKO76oxjBfBFh7FcB8VVURuZ7oVTVfRK6JiLuqZpTSLxHVEUz0VJsU3qcPgnnkegbATACXAXwiIv0A3APgTlXNEpHdAJytjl8LYDQAb5hH+ALgX6r6Qil9ltUmx+p1PswXAlKBz1Tc8UVlw+pziEg3mGcK+ojIe5Z9x63aOwG4WoEYiMjAOHVPtcl+mH9Glqqq+aqaCvMCtzthnsr3AHDJkuQ7ALijyPFfAHgY5mS/FuaFa6NFpAkAiIiniPgVOaY8bYr6HsBQEXEWETcA91m2ZwBwr+iHVtVLMN/LL0z2CwEMVVV/VfUH0BWWEb3lVkWKquZWtB8iMiYmeqpNjsM8jf1DkW3pqnoBwHYA9iJyDMArRdpBVSNhTrRnVTVJVaNgXtD2neWYHTD/XM36mDLbFKWqBwFsAnAUwHqY1wKkq+pFAPtFJMJqMV55fQfzCL4/gHqq+h+r/s4BqCcingD+COCbCp6biAxMVLW6YyAyHBFxU9VMEXEFsBfARFU9dAvnCwYwQ1XHldFuPYAXVDX2ZvsiImPhPXoi21ghIoEw3z//160keQCwLCbcJSJ21r+lt2ZZnb+RSZ6IrHFET0REZGC8R09ERGRgTPREREQGxkRPRERkYEz0REREBsZET0REZGBM9ERERAbGRE9ERGRgTPREREQG9v/PfQqiu4xWDQAAAABJRU5ErkJggg==\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 XMM-13hr')"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfoAAAEYCAYAAAC0mTTAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd8VGX2P/DPSSMJGXogkIQUSA8JzYDACoLSFlkBK4prR1l+isiuKLvuroJiAV3WtbcvCCsqoCACglJlAQNSUgmhQ4KhhYGEEnJ+f9w7OKmEJJMy+bxfr3kx89x27p0hZ+7z3LlHVBVERETknFxqOwAiIiJyHCZ6IiIiJ8ZET0RE5MSY6ImIiJwYEz0REZETY6InIiJyYkz0REREToyJnoiIyIkx0Ts5EdkvIjdV07peFpEJ1bGuUta9RURiHLFuour8f0BU3zDR1zDzD06+iFhF5LSIbBSRx0SkWt4LR/1BExFfAPcBeM+ubY2InBeRs+YjvQqbeB3AC1WNk6qfiPiYn6vRdm0WETkoIreZr/eLyEURaVVs2e0ioiISfC3zlRHHeBFJFJELIvJpKdM/E5EsETkjIrtF5OEq7DaR02Cirx23qKoFQBCA6QCeAfBR7YZ0VfcD+E5V84u1j1dVH/MRUYX1LwZwo4i0rcI66jQRcavtGCpDVc8CeBTAv8wvfADwKoBEVf3KbtZ9AO62vRCRTgC8SlllRecr7iiAqQA+LmP6ywCCVbUJgOEApopItwqst4j6+j4RlYWJvhapaq6qLgZwJ4A/ikgsAIhIOxFZICI5IrJPRJ6wLWOeET0rIikickpEPhERT3PaHADtASwxz7D/Yi7WWUR2ikiuiMy3zW8u84yIHDF7GNJFZEAZ4Q4BsLay+yoi7iIyzYz/knn2piKywzwW5wFsBTCwlGUni8hXxdr+JSKzrnEfICKBIrLQPLYnROQtsz3K7KE4LSLJIjK8Its2X1/t/XpGRHYCOCcibuY6M814U0RkhN38XUXkF3Pal+b7NbUi2yplX0vdJ7u4JpX1uShOVb8HsBTALBHpB+AOAH8qNtscGL0+Nn8EMLuU1VV0vuIxLFTVrwGcKGN6sqpesL00Hx3sZinv/0GJ9+lq8RDVG6rKRw0+AOwHcFMp7QcBPA7jy9dWAM8D8AAQCmAvgEF2yycBCATQAsBPAKaWtX7z9RYA7cz5UwE8Zk6LAHAIQDvzdTCADmXEnQPgumJta8z242Yc/crZ71cAbDLjbgxgFYCFAELt5pkFYGYpywYByAPQxHztCiALQM9r3AdXADsAvGHG4AmgDwB3AHsAPGce8/4ArOa6y9y2+boi79d2c7+9zLbbzffDBcaXvHMA2prLHwDwpBnTSAAXbe/v1bZVbF/L3KerfS7KeQ+bm/t+HMADpX2uAaQDiDKP0yHz+CmMM+0Kz3eVOKYC+LSMaW+b75cC2AbApyL7W9r7xAcfzvLgGX3dcRTGH6DrAPiq6guqelFV9wL4AMBddvO+paqHVPUkgGmw6wYtwyxVPWrOvwRAZ7P9MoBGAKJFxF1V96tqZhnraAYjUdh7Bkay8QfwPoyehA7FFxQRC4AnAIwx4z4HYAGAFub+2VjN7RShqgdg/NG+1WzqDyBPVTdd4z4kwPhD/2dVPaeq51V1A4wvDD4AppvH/EcA3wK4+yrbBir2fs0y9zvf3J8vzfejUFXnA8gwY+sJwM2c/5KqLoSRnHAN27Ipc5+KxVXa56JUqnoKQDIAbxhf0kpjO1u/GUAagCNVnO+aqOo4ABYAvzNjvGA3+Wr7W+R9InIWTPR1hz+AkzDObNqZ3a2nReQ0jLOyNnbzHrJ7fgBG8ipPtt3zPBgJAKq6B8AEAP8A8KuIfC4iZa3rFIw/oFeo6mZVtarqBVX9Pxhn9UNLWfYGAHtVNcOurXmxuGCu/3QZ25+H35LUaPP1te5DIIADqlpQrL0dgEOqWmjXdgDGe1Lmtk3X+n5BRO4T4+Iz2/yxAFqZcRxRVS1j2Ypsq6L7BJTxuSiLiNwLo8dkFYwemtLMgXGM7kf53fHlzici98hvF3kuKy+u4lT1svkFLgBGL5nN1fb3EIicEBN9HSAi18H4A7wBxh+bfarazO5hUVX7BBpo97w9jN4AG/skcVWqOk9V++C3rtOy/oDvBBB+tdUBkFLafWF8UQAAiIgAGAHjDNNeFIyu9dJ8CaCfiASYy15JttewD4cAtC9l/PUogEAp+suH9vjtLLPMbaNi79eV90REgmCchY8H0FJVm8EYihEY3eL+5vGxsX+vK7Ktiu7TNRGR1jCGPB4BMBbAHSJyQ/H5zB6QfTC+8JV11n/V+VR1rv52keeQysQMo3ekRA9TOa7p/w5RfcFEX4tEpImIDAPwOYDPVHUXjK7aM+aFQV4i4ioiseaXAZs/iUiAiLSAcUY3327aMRjd6RXZfoSI9BeRRgDOA8iH0RVemu8A9LVbtpmIDBIRT/MCs3tgnLmvKGXZJABdRaSziHjBuDpa7eM2Y+gGYGVpG1fVHBjXBHwCI9mlVmIftsBIptNFpLEZe28Am2GMk/9FjIsG+wG4Bcb7Uua27dZ5tffLXmNz33PM+B+AcUYPAP8zYx9vHtM/wOjSr8y2yt2nSngLwNequlpVswD8BcAH5nEv7iEA/c0hmvJUdD4AxtXw5gV0rgBcbZ89c1prEblLjJ8CuorIIBi9MD9WcP+InBYTfe1YIiJWGGdoUwDMBPAAYHQ7wviD3BnGGc9xAB8CaGq3/DwA38O4EGsvjIuTbF4G8Feza3fSVeJoBOPnfcdhdGu2hvHFoTSzAQw1EzVgXOw1Fb9djPf/ANyqqiV+S6+qiTCuJfjOjNcPwFBVvWQ323AAa1T1aPHl7cyDcSGX/Rl1hffB7th2hHHx42EAd6rqRXP7Q8z1vA3gPlVNu8q2K/p+2c+fAmAGjKR+DEAnGEMeMOMYCSMBngZwL4xejwvXuq0K7lOFiMitMC5a/LPd+j+EcfyeL2XbmeZ7Xq6KzmfnrzC+yE2GcWzyzTbA+PL0uBnTKRj3ZZigqt9cw/qJnJIUHQ6kuk5E9gN4WFVX1cK2XwLwq6q+6YB1bwbwkKomVfe66zPzuLyrqp/UdixEVD/xt6JUYapa1tl+day7h6PWXZ+ISF8YPz07DuAeAHEAltdqUERUrzHRE9UtEQC+gHFFeCaA28wxcSKiSmHXPRERkRPjxXhEREROrF503bdq1UqDg4NrOwwionpl69atx1XV9+pzlrl8azc3tw9h/ASUJ4Z1UyGApIKCgoe7dev2a2kz1ItEHxwcjMTEa/kVDhERiciBqizv5ub2oZ+fX5Svr+8pFxcXjvPWQYWFhZKTkxOdnZ39IYyf1JbAb2hERFSWWF9f3zNM8nWXi4uL+vr65uK3G2+VnKcG4yEiovrFhUm+7jPfozLzucMSvRh1v1eLSKpZC/vJYtMniVGPvJWjYiAiImroHDlGXwDgaVXdZpYp3SoiK1U1RUQCYZSnPOjA7RMRETV4DjujV9UsVd1mPrcCSMVvJTLfgFEUg11CRERO4PUV6W1WpR4rUsp6Veoxy+sr0ksro1wnzJo1q+V9993Xvnj7q6++6vvWW2+1LN6enp7uERYWFlMd277zzjuDtm7d6lkd67qaGhmjF5FgAF0AbBaR4TBqbpdVjtS2zKMikigiiTk5OdUb0IY3gX3rirbtW2e0ExHRNevcvlnexC+2h9qS/arUY5aJX2wP7dy+WV5tx3at/vKXv+SMHz/+hCO3MX/+/APdunU778ht2Dg80YuID4AFACbA6M6fglIqXhWnqu+randV7e7rW+mfgZbOvyvw5f2/Jft964zX/l2rdztERA3ETVFtrDPv6Lx34hfbQ/+5JLndxC+2h868o/Pem6LaWKu87ptu6hATExPVsWPHmNdff70VAHh7e3d55JFHAqKjo6Ouv/768KNHj5Y5FJ2QkBCxbt06bwDIyspy8/f371R8ns8//7xp586dI7OystwmTpzY7vnnn28DAOvXr/eOiIiI7ty5c+TMmTNb2+ZPTEz07NSpU1RkZGR0eHh49K5duxqdOXPGpV+/fh0jIiKiw8LCYj744IPmFYnJ0Rya6EXEHUaSn6uqCwF0ABACYIdZhS0AwDYR8XNkHMV9fG4Pttz8nJHcf5wGfHk/ttz8HD4+t6cmwyAicio3RbWxjuoakPPJT/vbjuoakFMdSR4A5s6duz85OTl1+/btKe+9916b7Oxs1/z8fJeuXbvmpaSkpPbu3ds6efLkdpVd/+zZs5u99tprfitXrsxo27Ztgf20hx56KHjmzJkHt2/fXqTE87///W/fcePGHUtLS0vZuXNnakhIyMWFCxc28fPzu5Senp6SkZGRPHLkyDOVjak6OfKqewHwEYBUVZ0JAKq6S1Vbq2qwqgbDqB3dVVWzHRVHaWJbxmJS+mxsiR0GrHsVW2KHYVL6bMS2LPNniEREdBWrUo9ZFmw77PtA7+CsBdsO+xYfs6+sV155pU1ERER0t27dorKzs92Tk5M9XVxc8PDDD58EgAcffPDEli1bfCqz7o0bN1pmzJjht3LlygxfX9/L9tNOnDjharVaXX//+9+ftW3HNu36668/N2PGjLZTpkzxy8jI8PDx8dGuXbvmr1+/vsnjjz/uv3z5cp+WLVteLr692uDIM/reAMYA6C8i283HUAdur8IS2ibg9Yj7MOnYj3grfggmHfsRr0fch4S2CbUdGhFRvWQbk595R+e9f78l5qitG7+qyf7bb7+1rF271pKYmJiWnp6eEhUVlZ+fn18idxnnlqVzc3PTy5eNnJuXl1dkxvbt2184d+6ca1JSUokL41S1zPU+9thjJ7/55ps9Xl5ehUOGDAlfvHixJS4u7sK2bdtSOnXqlD9lyhT/SZMmtb3G3XUIR151v0FVRVXjVLWz+fiu2DzBqnrcUTGUad86JKx8CXcED8V7Z5JxR/BQJKx8qeQFekREVCHbD572th+Tt43Zbz94ukrj0KdPn3Zt2rTpZYvFUvjLL7947tixozEAFBYW4pNPPmkOAJ9++mnLhISEMocJAgMDL2zZsqUxAMydO7fIuHlAQMDFBQsW7HnggQdCEhMTiyT7Vq1aXfbx8bm8YsUKH3M7LWzTUlJSPKKioi789a9//XXgwIGnt2/f7rV//353i8VSOG7cuJMTJkw4tn379hoZg7+aenGv+2p3ZBu23PwcvkifjbFxY/FF+hdIuPk5JBzZBoTcUNvRERHVO5MGRRwr3nZTVBtrVcfpR40alfv+++/7hoeHR3fo0OF8fHz8OQDw8vIqTE5O9oqJifGzWCyXFy5cuLesdUyePPnYnXfeGfr555+3/N3vfldi3Dw+Pv7C7Nmz9955550dFi9eXORirY8++mj/ww8/HOzl5VXYv3//K8vOmTOnxZdfftnSzc1NfX19L7388stHN2zY0PjZZ58NcHFxgZubm7799tvl1hoorxeiOtWLevTdu3fX6ixqsyVrCyatnYTX+76OhLYJJV4TETkDEdmqqt0ru/yOHTv2x8fH13yvawV4e3t3ycvL+6W246is8PDw6MWLF++JjIy8WB3r27FjR6v4+Pjg0qY1yHvdJ51IKpLUE9om4PW+ryPpRFItR0ZERM6uV69eYREREfnVleSvpkF23T8Y+2CJtoS2CTybJyKqJ0o7mx8zZkz7n3/+ucjV948//vixJ5980qE3vynPzTff3OHQoUON7NumTZt2eNSoUTX207sGmeiJiMj5zJkzp87VT1m5cmVmbcfQILvuiYiIGgomeiIiIifGRE9EROTEmOiJiIicGBM9ERFV3Q8vtkH6sqK3u01fZsEPL9Z6Pfr//e9/Xp07d44MDw+P7t+/f8eTJ09eyX3PPvusX/v27WODg4NjFyxY0ORa1z1hwoR2X3/9dbXc099RmOiJiKjqArrnYdFjoVeSffoyCxY9FoqA7rVej/6RRx4JnjZt2uHdu3enDB8+/NQ///lPPwDYunWr58KFC1ukp6cnL1++fPeECRPaFxQUXG11VxQUFODNN988euutt1ZLlT5HYaInIqKqixhixYh392LRY6FYNrkdFj0WihHv7kXEkFqvR79//37PIUOGnAWAYcOGnfn222+bA8BXX33VbOTIkSe9vLw0MjLyYlBQ0IU1a9Y0Tk9P9wgJCYkZOXJkcHh4ePTgwYNDrVarCwD4+/t3mjRpUttu3bpFfPzxx81HjRoVbLvnvr+/f6fx48f7d+7cOTI2NjZqw4YN3n369AkLDAyMffXVV31t8fztb39rExsbGxUeHh791FNPlVte989//nPbkJCQmF69eoXdcsstIc8///w195Aw0RMRUfWIGGJF/N052PxOW8TfnVMdSR6oej36sLCw/Hnz5jUDgM8++6xFdna2BwAcOXLEIzAw8Mrd6dq1a3fx0KFDHoDx5eCxxx7L2b17d4rFYil87bXXriRqT0/Pwq1bt6Y/+uijp4pvKzAw8OL27dvTevTocfbBBx8MXrJkSebmzZvTpk+f3g4AFi5c2GTPnj2eO3fuTE1NTU3Zvn2797Jly0otsbtu3TrvJUuWNN+1a1fK0qVLM3fu3Nm4MsePiZ6IiKpH+jILdvzXFz0ez8KO//qWGLOvpKrWo//444/3v/POO74xMTFRVqvVxd3dXQGjDG1xIqIA4Ofnd3HgwIHnAGDMmDEnNm7ceGX99913X4kEb3PHHXecBoBOnTrlde3a9Vzz5s0L27VrV9CoUaPC48ePuy5fvrzJunXrmkRHR0fHxMREZ2ZmeqalpZUokQsAa9as8RkyZMhpHx8fbd68eeHNN998ukIHrBjeGY+IiKrONiZv664P7Wutju57+3r0FoulMCEhIeJa69F36dLl/E8//ZQBADt37mz0/fffNwOMErW2M3gAOHr0qEdAQMCl0tZn/9pisRSWtS1PT08FABcXF3h4eFz5JuHi4oJLly6JqmLChAlZf/7zn69aLKi6is7xjJ6IiKrucKJ3kaRuG7M/nFjr9eiPHDniBgCXL1/G3//+97YPPfTQrwAwatSo0wsXLmyRn58vaWlpHvv37/fs16/fOQDIysryWLVqVWMAmDdvXotevXqdrcp+2AwZMuTMnDlzWuXm5roAwL59+9xt8RXXr1+/sytWrGial5cnubm5LqtWrWpWmW3yjJ6IiKpuwN9K1KNHxBBrVcfpq6Me/ccff9zio48+ag0AQ4cOPfXEE0+cAIDu3bufv/XWW0+Gh4fHuLq6YubMmQfc3Iy0GBoaev7jjz9uOW7cuKCQkJALkyZNyqnKftiMHDnyTHJysud1110XCQDe3t6Fc+fO3efv71/icv++ffvmDR48ODc6OjrG39//Qlxc3LmmTZtevtZtNsh69EREDQHr0VdOenq6x7Bhw8IyMjKSHbH+a5Gbm+vStGnTQqvV6nL99ddHvPvuuwf69OlT4ieL5dWj5xk9ERFRHXXvvfcGZWRkeF24cEHuuuuuE6Ul+athoicionrHkfXoIyIiLtbk2Xx2drZrv379Ioq3r1mzJn3JkiX7qrp+JnoiInIKdbEefUX4+fldTktLS3HU+nnVPRERkRNjoiciInJiTPREREROjImeiIjIiTHRExFRlc3aNqvNmkNritzbfs2hNZZZ22axHn0tY6InIqIqi/ONy5uyYUqoLdmvObTGMmXDlNA43zjWo68GhYWFuHz5mm+KB4CJnoiIqkG/wH7WaX2m7Z2yYUro9C3T203ZMCV0Wp9pe/sF9mM9+krWo09PT/cIDQ2Nuffee9uble48ypq3PEz0RERULfoF9rPe0uGWnLmpc9ve0uGWnOpI8kDDrUdvi+OBBx44kZqamhIeHn6xrPnKw0RPRETVYs2hNZYlmUt874m6J2tJ5hLf4mP2ldVQ69EDQNu2bS8OGDDg3FUPUjl4ZzwiIqoy25i8rbu+Z9ue1urovm/I9egBo7pdReYrD8/oiYioynbm7PS2T+q2MfudOTtZj97OtdSjry48oycioip7ousTJerR9wvsZ63qOH1DrkdfXRxWj15EAgHMBuAHoBDA+6r6LxF5EcAfzLZfAdyvqkfLWxfr0RMRXTvWo6+culSPvqLKq0fvyK77AgBPq2oUgJ4A/iQi0QBeU9U4Ve0M4FsAzzswBiIiogbNYV33qpoFIMt8bhWRVAD+qmpfiq8xAMd0KRARkdNqKPXo/fz8KneXHDs1MkYvIsEAugDYbL6eBuA+ALkAbixjmUcBPAoA7du3r4kwiYioHmM9+tI5/Kp7EfEBsADABFU9AwCqOkVVAwHMBTC+tOVU9X1V7a6q3X19fUubhYiIiK7CoYleRNxhJPm5qrqwlFnmARjlyBiIiIgaMoclejHuLvARgFRVnWnXHmY323AAaY6KgYiIqKFz5Bh9bwBjAOwSke1m23MAHhKRCBg/rzsA4DEHxkBERNSgOfKq+w0ASrsn4XeO2iYREdWOX998s41XfHye5cYbr9wgx7p6tSV/xw7v1hMmlLiZjrOaMGFCu379+llronRtRfEWuEREVGVe8fF5R5+ZHGpdvdoCGEn+6DOTQ73i42u9Hn1Nqqn69NeCiZ6IiKrMcuON1navTN979JnJodkvvdTu6DOTQ9u9Mn2v/Rl+ZVW1Hn1CQkLEQw89FNi9e/eI0NDQmLVr13oPHDiwQ1BQUOwTTzxxpbztP/7xjzZhYWExYWFhMS+88EJrADhz5oxLv379OkZERESHhYXFfPDBB80BYNKkSW1jY2OjwsLCYu6+++6gwkKj9ox9ffq1a9d6d+nSJTIiIiK6U6dOUadOnSo15955551BkZGR0ZGRkdHNmzePf/rpp9tW9ZjZY6InIqJqYbnxRmvTW/+Qc2r2nLZNb/1DTnUkeaDq9egBwMPDozAxMTH9gQceyLn99ts7fvDBBwfT0tKS58+f3yo7O9t1/fr13vPmzWu5devW1MTExNTZs2f7/vTTT14LFy5s4ufndyk9PT0lIyMjeeTIkWcA4M9//vOvSUlJqRkZGcn5+fkun3/+eVP77Z0/f17uueeeDm+++ebB9PT0lLVr16b7+PiUWolu/vz5B9LS0lIWL168p1mzZgVjx469phv8XA0TPRERVQvr6tWW3K+/8W1+35is3K+/8bV141dVVevRA8CIESNOA0B8fHx+x44d84OCgi55eXlpYGDghb1793qsWbPGZ+jQoaebNGlS2LRp08Lf//73p1avXm3p2rVr/vr165s8/vjj/suXL/dp2bLlZQBYtmyZJS4uLjI8PDx648aNlqSkJC/77e3cudOzdevWl/r27ZsHAC1atCh0d3cvM768vDwZNWpUhzfeeONgeHj4xSoesiIaZKK3rj2E85mni7SdzzwN69pDtRQREVH9ZhuTb/fK9L1+zz131NaNX9Vkb1+PPj09PSUqKir/WuvRA0XrxDdq1KhInfiCgoIyC7zFxcVd2LZtW0qnTp3yp0yZ4j9p0qS2eXl58vTTTwctXLgwc/fu3Sn33nvv8fPnzxeJSVUhIhW+xfuYMWOCbrnlllOOGN9vkInePcCCk/NSryT785mncXJeKtwDquXLJxFRg5O/Y4e3/Zi8bcw+f8eOWq9HXxH9+/c/+9133zWzWq0uZ86ccfnuu++a33jjjdb9+/e7WyyWwnHjxp2cMGHCse3bt3vn5eW5AICfn19Bbm6uy5IlS5oXX198fPz5Y8eOeaxdu9YbAE6dOuVy6dKlUrf98ssv+549e9b1pZdeyq7KPpSlQdaj9+zQDC1GR+HkvFQ07tEW5zZnocXoKHh2aFbboRER1Uul/YTOcuON1qqO01dHPfqK6NOnT97o0aNPdO3aNQoAxowZk9O7d+/8BQsWNHn22WcDXFxc4Obmpm+//faBVq1aXb7nnntyoqOjYwICAi7aYrLn6empc+fOzXziiSfanz9/3sXT07Nw3bp1u5s2bVpinP6tt97yc3d318jIyGgAePDBB3P+8pe/5FRlf+w5rB59dXJUPfrc7/fD+uMhWPoHounA4GpfPxFRbWI9+oajturR12nnM0/j3OYsWPoH4tzmrBJj9kRERM6gQXbd28bkbd31jTo0K/KaiIjqNkfWo3eUBQsWNJkyZUqAfVtgYOCFlStXZjpyuw0y0V86bC2S1G1j9pcOW5noiYjqqbpej37UqFFnRo0a5bC682VpkIne0jewRJtnh2ZM8kRE5HQa7Bg9ERFRQ8BET0RE5MSY6ImIqMo2fZPZZt/O40XuOrZv53HLpm8y29RWTGRgoicioiprE9I074dPU0JtyX7fzuOWHz5NCW0T0tQhZWoTEhIi1q1bV+Kue3379u14/Phx1+LtEydObPf8889X+UvH/v373QcPHhxa1fXUpAZ5MR4REVWvkLhW1gH3R+/94dOU0Iiefjnpm7J9B9wfvTckrlWN1mZfu3btHkeuPzg4+NLy5curdBe+msYzeiIiqhYhca2sET39cnb+eLhtRE+/nOpI8unp6R5hYWExttfPP/98m4kTJ14pSXv58mWMHDky2FZX3t/fv1NWVpYbADzzzDN+wcHBsb169QrPyMhoZFtm6tSprTt06BATHh4ePWzYsFAAWLp0qY+tJnxUVFR0WbXji8dTH/CMnoiIqsW+ncct6ZuyfeP6B2Slb8r2DYhsYXXkGf2lS5fk1ltvDYmOjs5/5ZVXihSEWb9+vfeiRYta7Nq1K+XSpUvo3LlzdJcuXfIAYNasWX4HDhzY5eXlpbZu/hkzZvjNmjXrwMCBA8/l5ua6eHt7l1o7vj7iGT0REVWZbUx+wP3Re393R/hRWzd+8Qv0qtO4ceOCSkvyALB69WqfoUOHnrZYLIUtWrQoHDhw4JX7nEdEROSPGDEi5O23327h7u6uANCzZ8+zkyZNCpw6dWrr48ePu5ZXO76+YaInIqIqO7Yv19t+TN42Zn9sX26VytS6ublpYeFvJ9f2dd+7d+9+dv369U3y8vJKLUZfVo361atXZ/zpT3/K2bp1a+P4+PjoS5cu4aWXXsr+8MMPD+Tn57v06tUr6pdffvGsStx1SYNI9Cc+/BDnNm0u0nZu02ac+PDDWopqlXTLAAAgAElEQVSIiMi59PxDh2PFu+lD4lpZe/6hQ4nytdciICCg4OTJk27Z2dmu+fn5smLFiqa2aWPHjj0+cODA3GHDhnUoXuu9f//+Z5cuXdrs7NmzcurUKZeVK1c2A4wx/czMTI9bbrnF+vbbbx+2Wq2uubm5rsnJyY0SEhLyp02blt2pU6dzSUlJTPT1iWdsJxx56qkryf7cps048tRT8IztVMuRERFReRo1aqRPP/10VkJCQtSAAQM6duzY8bz99H/84x/H4uPj80aOHBly+fLlK+19+vTJGzFixMnY2NiYYcOGdUhISDgLAAUFBTJ69OiQ8PDw6NjY2OixY8cea9Wq1eVXX321dVhYWExERES0l5dX4W233ZZbVkwiUvfru9tpMPXobcm9+d134dR/P4f/G2+gcc8e1RQhEVHd48z16GvL+vXrvSdOnBj4888/p9d2LPZYjx5A45490Pzuu3D87XfQ/O67mOSJiOiarFu3znvMmDGh48ePr9JwRE1rMD+vO7dpM07993O0Gvc4Tv33c3gn9GCyJyKiUm3ZssXrvvvuC7Fv8/DwKNy/f39SbcVUWQ0i0du67W3d9d4JPYq8JiIispeQkJCflpZW47XjHaFBdN2fT9pVJKk37tkD/m+8gfNJu2o5MiIiIsdqEGf0LR9+uERb457suiciIufXIM7oiYiIGiomeiIiqrINn89uk7l1S5Hb3WZu3WLZ8Pls1qOvZUz0RERUZW3DIvOW/WdGqC3ZZ27dYln2nxmhbcMiq1yP/sUXX2wdFhYW07Fjx5gXXnihNQAcO3bMtVevXmFBQUGxvXr1CsvJySlRg/5qunTpElnV2OoDJnoiIqqyDt0SrEP+9PTeZf+ZEbr60/fbLfvPjNAhf3p6b4duCVWqXvfzzz97zp4923fbtm2pqampycuXL2+2a9euRn//+9/b9uvXz3rgwIGkfv36WZ9//nm/iq6zoKAAAPDLL7+kVSW2+oKJnoiIqkWHbgnWmBsG5GxbtrhtzA0Dcqqa5AFg165dXl27dj1rsVgK3d3d0bt3b+v8+fObLV++vNnYsWNPAMDYsWNPLFu2rDkATJw4sd2tt94a0rNnz/CgoKDYGTNmtAKAb7/91tKjR4/wW265JSQiIiIGALy9vbvYpl133XURQ4cODQ0ODo4dN26c/zvvvNOiU6dOUeHh4dHJycmNAODo0aNugwYN6hAbGxsVGxsb9f333zcuK+6jR4+69erVKyw6Ojpq9OjRQe3ateuUlZVVKxfAOyzRi0igiKwWkVQRSRaRJ83210QkTUR2isgiEWnmqBiIiKjmZG7dYkle94Nv1yHDs5LX/eBbfMy+Mjp37py/efNmS3Z2tqvVanVZuXJl00OHDnmcOHHCLSgo6BIABAUFXTp58uSVJJqamuq1atWqjE2bNqW99tpr7fbv3+8OADt37mz82muvHcnMzEwuvp20tDSvd95551BqamryV1991XL37t2eu3btSh0zZszxGTNmtAaAsWPHBk6cOPFYUlJS6qJFizIfe+yx4LLinjx5cru+fftaU1JSUkeOHHkqKyvLo6rHorIc+e2iAMDTqrpNRCwAtorISgArATyrqgUi8gqAZwE848A4iIjIwWxj8rbu+vadOluro/u+a9eu55988sns/v37h3t7exdGR0fnubmVn7qGDBly2sfHR318fAquv/76M+vXr2/cvHnzy3FxceciIyMvlrZMp06dztm+OLRv3/7CkCFDcgEgPj4+f+3atRYA+Omnn5pkZGR42ZY5e/as66lTp1yaN29eWHx9W7Zs8fn666/3AMBtt912pkmTJpeLz1NTHHZGr6pZqrrNfG4FkArAX1W/V9UCc7ZNAAIcFQMREdWMrIw0b/ukbhuzz8pIq1I9egB46qmnjqekpKQmJiamt2jR4nJYWNj5li1bFhw4cMAdAA4cOODeokULW14pUYfe9trb27tEQrZp1KjRlQpvLi4u8PT0VNvzy5cvCwCoKhITE1PT0tJS0tLSUn799dedpSV527x1RY2M0YtIMIAuADYXm/QggGVlLPOoiCSKSGJOTo5jAyQioirpc9d9x4qfuXfolmDtc9d9VS4Ac+TIETcAyMjI8Fi6dGmzhx566OSgQYNOv/feey0B4L333ms5ePDg07b5ly1b1iwvL0+ys7NdN23aZOnTp8+5qsYAAH369DnzyiuvtLa93rhxo1dZ8yYkJJydM2dOCwBYuHBhkzNnzlzzrwKqi8MvDBARHwALAExQ1TN27VNgdO/PLW05VX0fwPuAUabW0XESEVHdNHz48A6nT592c3Nz0zfffPOgr6/v5X/+859ZI0aM6BAUFNSqXbt2F7/++utM2/xdunQ5N2DAgLCjR496TJo0KSs4OPhSUlKSZ1XjeP/99w89/PDD7cPDw6MvX74sPXr0sPbq1etgafNOnz796G233RYaHR3d/Prrrz/r6+t7qVmzZrXSfe/QevQi4g7gWwArVHWmXfsfATwGYICqXvU3ltVRj56IqKFpiPXoJ06c2M7Hx+fyCy+8UKulZPPz88XNzU3d3d2xatWqxuPHjw9yZJGc8urRO+yMXoxBkY8ApBZL8oNhXHzXtyJJnoiIqL7Zs2ePxx133NGhsLAQ7u7u+t577+2vrVgc2XXfG8AYALtEZLvZ9hyAWQAaAVhpXiCxSVUfc2AcRETUQMycOfNoTW7vX//6V8t33nmnyG1+r7vuurNz5sw5mJqaWifK3FYo0YtIIwCjAATbL6OqL5S1jKpuACClTPru2kIkIiKqm5588skTTz755InajqM8FT2j/wZALoCtAC44LhwiIiKqThVN9AGqOtihkRAREVG1q+jv6DeKSCeHRkJERETVrtwzehHZBUDN+R4Qkb0wuu4FgKpqnONDJCKiui53xf42Hu0teV5RLa/cNCc/9YTl4kGrd9NBwbX6U7eG7mpn9MMA3AJgCICOAAaar23tRERE8GhvyTv5xe7Q/NQTFsBI8ie/2B3q0d7CevS1rNxEr6oHVPUAgKm25/ZtNRMiERHVdV5RLa0t7gjfe/KL3aGnl2S2O/nF7tAWd4TvtT/DrwxnqEd/6dKlmthMmSo6Rh9j/0JEXAF0q/5wiIiovvKKamlt3LV1ztmfjrZt3LV1TlWTPFB/69FPnDix3d133x3Uu3fvsJEjR4ZU9ThUxdXG6J+FcZMbLxE5g99+F38R5n3oiYiIAKO7/ty2X319erfLOrftV99GHZtZq5rsO3funP/CCy/4Z2dnuzZu3FhXrlzZND4+/tzV6tFv3bo11Wq1unbp0iV61KhRuYBRj/6XX35JLq1UbVpamtdXX321t3Xr1gVBQUGdGjVqdHzXrl2pL774YusZM2a0/vjjjw/Z6tEPGjTobEZGhsegQYPC9u7dW6K2vc3OnTu9N2/enObj41Or9VrKTfSq+jKAl0XkZVV9toZiIiKiesY2Jm/rrm/UsZm1Orrv62s9egAYPHjw6dpO8kDFf0f/nIiMBNAHxlX461X1a8eFRURE9cnFg1Zv+6RuG7O/eNDqXdWz+qeeeur4U089dRwAxo8f7x8QEHDRVo8+KCjoUk3Xo69o8m7cuHGZ26tJFR2j/w+ManO7ACQBeExE/uOwqIiIqF5pOij4WPGE7hXV0lodP62rj/Xo65KKntH3BRCrZk1bEfk/GEmfiIjIoepjPfq6pEL16EVkIYCnzJ/VQUSCAExX1bsdHB8A1qMnIqoM1qNvOKqjHn1LAKkissV8fR2A/4nIYgBQ1eFVjpKIiIiqXUUT/fMOjYKIiKga1KV69DUZR3kqlOhVda3ZXR+mqqtExAuAm6pW+WYIRERE9VV9qEdfoavuReQRAF8BeM9sCgDAn9cRERHVcRX9ed2fAPQGcAYAVDUDQOtylyAiIqJaV9FEf0FVr9xNSETcYNw4h4iIiOqwiib6tSJiu+f9zQC+BLDEcWEREVF98sMPP7RJT0+32Lelp6dbfvjhhzZlLVPXOGvZ2oom+skAcmDcJGcsgO8A/NVRQRERUf0SEBCQt2jRolBbsk9PT7csWrQoNCAgoMr16GtKTZWtrWkVSvSqWgjj4rtxqnqbqn6gFbnTDhERNQgRERHWESNG7F20aFHosmXL2i1atCh0xIgReyMiIqr066z09HSPkJCQmDvvvDMoLCwsZvjw4SFff/21pWvXrpFBQUGxq1ev9j527JjrTTfd1CE8PDw6Pj4+cvPmzV4AsHTpUp/IyMjoyMjI6KioqOhTp0655Obmulx//fXh0dHRUeHh4dGfffZZM9u2bGVrAeCvf/1rm/Dw8OiIiIjocePG+ZcW2/79+91t64+MjIx2dXXttnv3bo+q7K8jXK1MrQD4O4DxMErUiohcBvBvVX2hBuIjIqJ6IiIiwhofH5+zefPmtj169MiqapK3OXTokOf8+fP3duvW7UBcXFzU3LlzWyYmJqbNmzev2bRp09r6+/tfjI+Pz1u1alXm4sWLLX/84x9D0tLSUmbMmOE3a9asAwMHDjyXm5vrYitqs3Tp0j0tWrQozMrKcuvRo0fk6NGjT7u4/Hbe+8UXXzRZunRp861bt6ZZLJbCY8eOuZYWV3Bw8KW0tLQUAHj55Zd9169fbwkPDy+1Ol5tutoZ/QQYV9tfp6otVbUFgB4AeovIUw6PjoiI6o309HTLjh07fHv06JG1Y8cO3+Jj9pXl7+9/ISEhId/V1RXh4eH5/fv3P+Pi4oKuXbvmHT58uNGWLVssDz300AkAGD58uPX06dNuJ06ccO3Zs+fZSZMmBU6dOrX18ePHXd3d3VFYWCgTJkwICA8Pj77xxhvDf/31V4/Dhw8XOelduXJlk3vvvfe4xWIpBIA2bdpcLi++77//vvHs2bN9//vf/+6vjv2tbldL9PcBuFtV99kaVHUvgHvNaURERFfG5EeMGLF3yJAhR23d+NWR7D08PEotIevq6orLly+XWrNFRPSll17K/vDDDw/k5+e79OrVK+qXX37xfO+991qcOHHCbdeuXalpaWkpLVu2vJSfn18kF6pqiVK3ZTlw4ID72LFjg+fPn5/ZtGnTOlGWtrirJXp3VS1R0EBVcwC4OyYkIiKqbw4fPuxtPyZvG7M/fPiwt6O33bNnT+snn3zSEgC+/fZbS/PmzQtatGhRmJyc3CghISF/2rRp2Z06dTqXlJTkmZub69qqVatLjRo10iVLlliOHj1aYkx98ODBZ+bMmdPKarW6AEBZXfcXLlyQkSNHhr744otH4uLiLjh2LyvvarfALW+soc6NQxARUe0YMGBAiWpxERER1uoapy/PK6+8cnT06NHB4eHh0V5eXoWffvrpPgB49dVXW2/cuLGJi4uLhoeH59922225p0+fdh0yZEjH2NjYqJiYmLyQkJDzxdd32223ndm2bZt3586do9zd3fWmm27Kfeutt44Un2/VqlWNk5KSGk+dOrXd1KlT2wHA8uXLM4KDgy85ep+vRbllas0L786VNgmAp6rWyFk9y9QSEV27hlimtqGqdJlaVS21u4KIiIjqh4qWqSUiImqwxowZ0/7nn3/2sW97/PHHj9X1ynUAEz0REdFV1aX68teqorfAJSKihqewsLCwYr8zo1pjvkdl/rSPiZ6IiMqSlJOT05TJvu4qLCyUnJycpgCSyprHYV33IhIIYDYAPxjfNN5X1X+JyO0A/gEgCkCCqvJyeiKiOqigoODh7OzsD7Ozs2PBE8O6qhBAUkFBwcNlzeDIMfoCAE+r6jYRsQDYKiIrYXzrGAngPQdum4iIqqhbt26/Ahhe23FQ1Tgs0atqFoAs87lVRFIB+KvqSgAVvr0gERERVV6NdMWISDCALgA2X8Myj4pIoogk5uTkOCo0IiIip+bwRC8iPgAWAJigqmcqupyqvq+q3VW1u6+vr+MCJCIicmIOTfQi4g4jyc9V1YWO3BYRERGV5LBEL8Yg/EcAUlV1pqO2Q0RERGVz5FX3vQGMAbBLRLabbc8BaATg3wB8ASwVke2qOsiBcRARETVYjrzqfgOMKnelWeSo7RIREdFveAMEIiIiJ8ZET0RE5MSY6ImIiJwYEz0REZETY6InIiJyYkz0REREToyJnoiIyIkx0RMRETkxJnoiIiInxkRPRETkxJjoiYiInBgTPRERkRNjoiciInJiTPREREROjImeiIjIiTHRExEROTEmeiIiIifGRE9EROTEmOiJiIicGBM9ERGRE2OiJyIicmJM9ERERE6MiZ6IiMiJMdETERE5MSZ6IiIiJ8ZET0RE5MSY6ImIiJwYEz0RUR3z7tpMbMw8XqRtY+ZxvLs2s5YiovqMiZ6IqI6JC2iK8fN+uZLsN2Yex/h5vyAuoGktR0b1kVttB0BEREX16tAKb43ugvHzfsG9Pdrjs80H8dboLujVoVVth0b1EM/oiYjqoF4dWuHeHu0x68c9uLdHeyZ5qjQmeiKiOmhj5nF8tvkgnujfEZ9tPlhizJ6oopjoiYjqGNuY/Fuju2DiwIgr3fhM9lQZDkv0IhIoIqtFJFVEkkXkSbO9hYisFJEM89/mjoqBiKg+2nk4t8iYvG3Mfufh3FqOjOojUVXHrFikLYC2qrpNRCwAtgK4FcD9AE6q6nQRmQyguao+U966unfvromJiQ6Jk4jIWYnIVlXtXttxUO1y2Bm9qmap6jbzuRVAKgB/AH8A8H/mbP8HI/kTERGRA9TIGL2IBAPoAmAzgDaqmgUYXwYAtC5jmUdFJFFEEnNycmoiTCIiIqfj8EQvIj4AFgCYoKpnKrqcqr6vqt1Vtbuvr6/jAiQiInJiDk30IuIOI8nPVdWFZvMxc/zeNo7/qyNjICIiasgcedW9APgIQKqqzrSbtBjAH83nfwTwjaNiICIiaugceQvc3gDGANglItvNtucATAfwhYg8BOAggNsdGAMREVGD5rBEr6obAEgZkwc4artERET0G6e8M96Wb77CwaSdRdoOJu3Elm++qqWIiIiIaodTJnq/DuH49s3pV5L9waSd+PbN6fDrEA4A2LBhA/bt21dkmX379mHDhg01HisREZEjOWWibx8bh2ETJuPbN6fjpy8+w7dvTsewCZPRPjYOAODv748vv/zySrLft28fvvzyS/j7+9dm2ERERNXOaevRt4+NQ/zAodi04HP0HHXXlSQPACEhIbj99tvx5Zdfonv37khMTMTtt9+OkJCQWoyYiIio+jnlGT1gdNfv+P479Bx1F3Z8/12JMfuQkBB0794d69atQ/fu3ZnkiYjIKTlloreNyQ+bMBm977j3Sje+fbLft28fEhMTccMNNyAxMbHEmD0REZEzcMpEn525u8iYvG3MPjtzN4DfxuRvv/129O/f/0o3PpM9ERE5G4eVqa1O1V2mdsOGDfD39y/SXb9v3z4cOXIEffr0qbbtEBHVJpapJcCJL8YrT2nJPCQkhOP0RETkdJyy656IiIgMTPREREROjImeiIjIiTHRExEROTEmeiIiIifGRE9EROTEmOiJiIicGBM9ERGRE2OiJyIicmJM9ERERE6MiZ6IiMiJMdETERE5MSZ6IiIiJ8ZET0RE5MSY6ImI6phtKw7gcPqpIm2H009h24oDtRQR1WdM9EREdUzr4CZY8UHSlWR/OP0UVnyQhNbBTWo5MqqP3Go7ACIiKiogojkGPRKLFR8kIfYGfyStO4JBj8QiIKJ5bYdG9RDP6ImI6qCAiOaIvcEfid/tR+wN/kzyVGlM9EREddDh9FNIWncE3YcGI2ndkRJj9kQVxURPRFTH2MbkBz0Six7DQ6904zPZU2Uw0RMR1TG/7j9TZEzeNmb/6/4ztRwZ1Ue8GI+IqI7pOiioRFtARHOO01Ol8IyeiIjIiTHRExEROTEmeiIiIifGRE9EROTEmOiJiIicmKhqbcdwVSKSA8BR1RxaATjuoHU7Un2NG6i/sTPumldfY68rcQepqm9tB0G1q14kekcSkURV7V7bcVyr+ho3UH9jZ9w1r77GXl/jJufErnsiIiInxkRPRETkxJjogfdrO4BKqq9xA/U3dsZd8+pr7PU1bnJCDX6MnoiIyJnxjJ6IiMiJMdETERE5sXqf6EVksIiki8geEZlcyvQgEflBRHaKyBoRCbCb9qqIJItIqojMEoO3iCwVkTRz2nS7+e8XkRwR2W4+Hq5LsZvta8x12mJsbbY3EpH55rY2i0hwXYlbRCx28W4XkeMi8qY5f7Ud8yrG/YqIJJmPO+3aQ8zjmWEeXw+zva4c77LinmuuM0lEPhYRd7O9n4jk2h3v5ysbtwNj/1RE9tnF2NlsF/MztcdcX9c6Fvd6u5iPisjXZnu1HnOiElS13j4AuALIBBAKwAPADgDRxeb5EsAfzef9Acwxn/cC8JO5DlcA/wPQD4A3gBvNeTwArAcwxHx9P4C36mrs5rQ1ALqXsr1xAN41n98FYH5dirvY8lsB3FCdx7yKcf8ewEoYZZ0bA0gE0MSc9gWAu8zn7wJ4vA4d7/LiHgpAzMd/7eLuB+DbOvAZLy/2TwHcVsr2hgJYZu5TTwCb61LcxZZfAOC+6j7mfPBR2qO+n9EnANijqntV9SKAzwH8odg80QB+MJ+vtpuuADxh/EduBMAdwDFVzVPV1QBgrnMbgABUv2qP/Srb+wOA/zOffwVggIjRC1CX4haRMACtYXzBqk5ViTsawFpVLVDVczD+8A82j19/GMcTMI7vrebzunC8S40bAFT1OzUB2IK69xkvM/Zy/AHAbHO3NgFoJiJt61rcImKB8bn5uhKxEV2z+p7o/QEcsnt92GyztwPAKPP5CAAWEWmpqv+D8R80y3ysUNVU+wVFpBmAW/Dbf2gAGGV2130lIoF1NPZPzC7Av9kllyvbU9UCALkAWtaxuAHgbhhnv/Y/B6mOY17puM32IWIM67QCcCOAQBjH77R5PIuvs9aPdzlxX2F22Y8BsNyu+XoR2SEiy0QkphIx10Ts08zPxBsi0ugatlfbcdvm/0FVz9i1VdcxJyqhvif60s6Qiv9ecBKAviLyC4C+AI4AKBCRjgCiYJzJ+APoLyI3XFmxiBuMLs1ZqrrXbF4CIFhV4wCswm9nbHUp9ntUtROA35mPMdewvdqM2+YuGMfdprqOeaXjVtXvAXwHYKMZ2/8AFFxlnbV+vMuJ297bANapqq0HZRuM+6PHA/g3qnbW6ajYnwUQCeA6AC0APHMN26vNuG3uRtHPeHUec6IS6nuiP4yi35YDABy1n0FVj6rqSFXtAmCK2ZYL41v1JlU9q6pnYYzt9bRb9H0AGar6pt26TqjqBfPlBwC61bXYVfWI+a8VwDwY3ZBFtmd+iWkK4GRdiduMKx6Am6putVtXdR3zqsQNVZ2mqp1V9WYYiSADRtGSZubxLL7OunC8y4obZlx/B+ALYKLdus6Y7w1U9TsA7uaZaWU4JHZVzTK75y8A+ASlfMbL2l5txg0A5ll/AoClduuqzmNOVEJ9T/Q/AwgT48pnDxhng4vtZxCRViJi289nAXxsPj8I4xu5m9l92RdAqrnMVBh/mCcUW5f9eN9w2/x1JXbzdStzWXcAwwAkmcssBvBH8/ltAH4s1j1ea3HbLVr8TKc6j3ml4xYRV/MPNEQkDkAcgO/N47caxvEEjOP7jfm81o93WXGbrx8GMAjA3apaaLcuP9twj4gkwPgbcaIScTsy9rbmvwLjmgj7z/h9YugJIFdVs+pK3KbbYVx4d95uXdV5zIlK0jpwRWBVHjCutN0N4yrZKWbbCwCGm89vg/GNejeADwE0MttdAbwHI3GkAJhptgfA6KZLBbDdfDxsTnsZQDKMcbjVACLrWOyNYVyxvtOM818AXM1pnjCuFN4D4+Kr0LoSt9169xY/ptV5zKsQt6cZbwqATQA6260z1Dyee8zja79MbR/v8uIuMNdn+4w/b7aPtzvemwD0qqXPeHmx/whgF4wE/xkAH7NdAPzH3NYulPLrk9qM25y+BsDgYm3Vesz54KP4g7fAJSIicmL1veueiIiIysFET0RE5MSY6ImIiJwYEz0REZETY6InIiJyYkz01CCYv1X+XEQyRSRFRL4TkfDajouIyNGY6MnpmTcjWQRgjap2UNVoAM8BaOOAbblW9zqJiKqCiZ4aghsBXFLVd20NqrodwAYReU2MuuG7xKwdLkYd+aG2ecWofz7KvOvZayLys1lQZaw5vZ+IrBaReTBu1AIR+VpEtopIsog8areuh0Rktxg1zD8QkbfMdl8RWWCu+2cR6V0jR4aInJ7b1WchqvdiYdwxsLiRADoDiAfQCsDPIrIORlnSOwF8Z94CdQCAxwE8BOO2qteJUTHtJxGx3d40AUCsqu4zXz+oqidFxMtc7wIYpXn/BqArACuMO7ztMOf/F4A3VHWDiLQHsAJGASAioiphoqeGrA+A/6rqZQDHRGQtjIpoywDMMpP5YBjV3fJFZCCAOBGx3du+KYAwABcBbLFL8gDwhIiMMJ8HmvP5wahVfhIARORLALbrBG4CEC2/laxvIiIWNYoTERFVGhM9NQTJ+K3wjL3SypFCVc+LyBoYRV/uxG+FdgTA/1PVFUVWItIPwLlir28CcL2q5pnr8ixreyYXc/78q+8OEVHFcYyeGoIfATQSkUdsDSJyHYBTAO40x959AdwAowANYHTfPwDgdzC60WH++7hZeQ8iEi4ijUvZXlMAp8wkH4nfSvFugVG9r7kYpWtH2S3zPYziJrb4Oldpj4mITDyjJ6enqmp2o78pIpMBnAewH0YZYh8Y4+QK4C+qmm0u9j2A2QAWq+pFs+1DAMEAtplX8ufAKJNa3HIAj4nITgDpMCqSQVWPiMhLADbDqG+eAiDXXOYJAP8xl3EDsA7AY9VyAIioQWP1OqIaJCI+qnrWPKNfBOBjVV1U23ERkfNi1z1RzfqHiGyHUUt9H4CvazkeInJyPKMnIn4EAPEAAAArSURBVCJyYjyjJyIicmJM9ERERE6MiZ6IiMiJMdETERE5MSZ6IiIiJ/b/AVD8MH2/xD6HAAAAAElFTkSuQmCC\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
}