{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This notebook was run with herschelhelp_internal version: \n",
      "708e28f (Tue May 8 18:05:21 2018 +0100)\n",
      "This notebook was executed on: \n",
      "2018-06-07 13:27:41.812811\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
   },
   "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
   },
   "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": {},
   "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": "markdown",
   "metadata": {},
   "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": 5,
   "metadata": {
    "collapsed": 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": 6,
   "metadata": {
    "collapsed": 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": 7,
   "metadata": {
    "collapsed": 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": {},
   "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": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "depths = Table()\n",
    "depths['hp_idx_O_13'] = list(field_moc.flattened(13))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=10</i>\n",
       "<table id=\"table140247224820792-856651\" 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(\"$('#table140247224820792-856651').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",
       "    $('#table140247224820792-856651').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": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "depths[:10].show_in_notebook()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": 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": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=10</i>\n",
       "<table id=\"table140248882108120-857202\" 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(\"$('#table140248882108120-857202').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",
       "    $('#table140248882108120-857202').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": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "depths[:10].show_in_notebook()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table masked=True length=10</i>\n",
       "<table id=\"table140248993958544-613991\" 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.2503921</td><td>1.8862768411636353</td><td>0.33097887</td><td>1.7756997346878052</td><td>0.20270073</td><td>2.926134943962097</td><td>0.2498269</td><td>2.9660435914993286</td><td>1.2030176</td><td>7.819149971008301</td><td>0.6296477</td><td>7.468801736831665</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.2503921</td><td>1.8862768411636353</td><td>0.33097887</td><td>1.7756997346878052</td><td>0.20270073</td><td>2.926134943962097</td><td>0.2498269</td><td>2.9660435914993286</td><td>1.2030176</td><td>7.819149971008301</td><td>0.6296477</td><td>7.468801736831665</td></tr>\n",
       "<tr><td>2</td><td>171728340</td><td>2683255</td><td>5.9318357</td><td>356.375198364258</td><td>12.443625</td><td>352.5657470703126</td><td>0.25039157</td><td>8.1373929977417</td><td>0.42282617</td><td>13.384291076660151</td><td>0.19149545</td><td>25.123155212402338</td><td>0.30158496</td><td>28.183526611328126</td><td>1.2030171</td><td>49.33133430480966</td><td>0.9557962</td><td>57.19062194824217</td></tr>\n",
       "<tr><td>3</td><td>171728341</td><td>2683255</td><td>5.9318357</td><td>356.375198364258</td><td>12.443625</td><td>352.5657470703126</td><td>0.25039157</td><td>8.1373929977417</td><td>0.42282617</td><td>13.384291076660151</td><td>0.19149545</td><td>25.123155212402338</td><td>0.30158496</td><td>28.183526611328126</td><td>1.2030171</td><td>49.33133430480966</td><td>0.9557962</td><td>57.19062194824217</td></tr>\n",
       "<tr><td>4</td><td>171728342</td><td>2683255</td><td>5.9318357</td><td>356.375198364258</td><td>12.443625</td><td>352.5657470703126</td><td>0.25039157</td><td>8.1373929977417</td><td>0.42282617</td><td>13.384291076660151</td><td>0.19149545</td><td>25.123155212402338</td><td>0.30158496</td><td>28.183526611328126</td><td>1.2030171</td><td>49.33133430480966</td><td>0.9557962</td><td>57.19062194824217</td></tr>\n",
       "<tr><td>5</td><td>171728343</td><td>2683255</td><td>5.9318357</td><td>356.375198364258</td><td>12.443625</td><td>352.5657470703126</td><td>0.25039157</td><td>8.1373929977417</td><td>0.42282617</td><td>13.384291076660151</td><td>0.19149545</td><td>25.123155212402338</td><td>0.30158496</td><td>28.183526611328126</td><td>1.2030171</td><td>49.33133430480966</td><td>0.9557962</td><td>57.19062194824217</td></tr>\n",
       "<tr><td>6</td><td>171728344</td><td>2683255</td><td>5.9318357</td><td>356.375198364258</td><td>12.443625</td><td>352.5657470703126</td><td>0.25039157</td><td>8.1373929977417</td><td>0.42282617</td><td>13.384291076660151</td><td>0.19149545</td><td>25.123155212402338</td><td>0.30158496</td><td>28.183526611328126</td><td>1.2030171</td><td>49.33133430480966</td><td>0.9557962</td><td>57.19062194824217</td></tr>\n",
       "<tr><td>7</td><td>171728345</td><td>2683255</td><td>5.9318357</td><td>356.375198364258</td><td>12.443625</td><td>352.5657470703126</td><td>0.25039157</td><td>8.1373929977417</td><td>0.42282617</td><td>13.384291076660151</td><td>0.19149545</td><td>25.123155212402338</td><td>0.30158496</td><td>28.183526611328126</td><td>1.2030171</td><td>49.33133430480966</td><td>0.9557962</td><td>57.19062194824217</td></tr>\n",
       "<tr><td>8</td><td>171728380</td><td>2683255</td><td>5.9318357</td><td>356.375198364258</td><td>12.443625</td><td>352.5657470703126</td><td>0.25039157</td><td>8.1373929977417</td><td>0.42282617</td><td>13.384291076660151</td><td>0.19149545</td><td>25.123155212402338</td><td>0.30158496</td><td>28.183526611328126</td><td>1.2030171</td><td>49.33133430480966</td><td>0.9557962</td><td>57.19062194824217</td></tr>\n",
       "<tr><td>9</td><td>171728371</td><td>2683255</td><td>5.9318357</td><td>356.375198364258</td><td>12.443625</td><td>352.5657470703126</td><td>0.25039157</td><td>8.1373929977417</td><td>0.42282617</td><td>13.384291076660151</td><td>0.19149545</td><td>25.123155212402338</td><td>0.30158496</td><td>28.183526611328126</td><td>1.2030171</td><td>49.33133430480966</td><td>0.9557962</td><td>57.19062194824217</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(\"$('#table140248993958544-613991').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",
       "    $('#table140248993958544-613991').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": 12,
     "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": {},
   "source": [
    "## III - Save the depth map table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "depths.write(\"{}/depths_{}_{}.fits\".format(OUT_DIR, FIELD.lower(), SUFFIX))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'90prime_g', '90prime_r', 'mosaic_z', 'ukidss_j'}"
      ]
     },
     "execution_count": 14,
     "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": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Passbands on XMM-13hr')"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAesAAAEgCAYAAACU3FvWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xl8VNXZwPHfMzPZEyCEsEOCbIIgiwiKVikqgrWolPrWXauVtlqtYFvbWlRQ6/KqrW/VgmvFDWurRVDcqohalK2yBWWRAFlYErKRbZbz/nFnYgjJZEIymTuZ5/v58MnM3Dv3PhMy88xzzrnniDEGpZRSStmXI9IBKKWUUio4TdZKKaWUzWmyVkoppWxOk7VSSillc5qslVJKKZvTZK2UUkrZnCZrZUsicqeIvNDRztVRichHInJdpONQqqPSZK1CJiK7RKRKRCpEZJ+IPCsiqZGOK9qJyBgRKRWRQfUeO0lESkQk239/l4jUiki3Bs/9r4iYevs9578/vcF+f/I/fnWQOOaLyEYR8YjInQ22fde/rUREikTkdRHp07pXrpQKlSZr1VLfN8akAmOBk4HbIxxP1DPGrAceA54USxzwDDDXGLOr3q7fAJcE7ojISCCpkUN+DVxVbz8X8ENgRzOhbAd+DSxrZNsW4FxjTBegN7ANeKKZ4zXKH49SqgU0WatjYozJA94GRgCIyDUikiMi5SKyU0RmBfYVkW4istRflRWLyEoRcfi3/UZE8vzP+0pEzqp3mkQRWezftk5ERtU75m0issO/bYuIXFRv29Ui8omI/K+IHBKRb0RkWr3tA0Rkhf+57wHd6m1LFJEX/NVjiYisFpEejf0ORGSYv/m3REQ2169m/RXuYyKyzH+ez0VkYJBf6V1AL+B64HdABfCXBvssAq6sd/8q4PlGjvUmcJqIpPvvTwU2AIVBzo8x5m/GmLeB8ka27TPG5Nd7yAsMarBbloh86n+97wZaAUQk21/VXysiu4F/B4tDKXU0TdbqmIhIP+A8YL3/of3A+UAn4BrgEREZ6982B9gLZAI9sJKREZGhwI3AycaYNOBcYFe901wA/B3oCrwEvOGvOsGqEr8DdMZKdC+ISK96z50AfIWViB8AnhYR8W97CVjr3zafelWo/3ZnoB+QAfwUqGrk9cdhJcV3ge7AL4AX/a8p4BJ/bOlYVes9DY8TYIypAa4F7vf/vq41xvga7LYK6OT/kuAE/gdorK+9GlgC/Mh//0oaT+otIiL9RaQE6/dxK9bvtb5Lsf7vuwPx/n3qOxMYhvX/rJRqAU3WqqXe8H9gfwKsAO4FMMYsM8bsMJYVWEnsO/7nuLGqxixjjNsYs9JYk9J7gQRguIjEGWN2GWPqN9WuNca8ZoxxAw8DicAp/vP93RiTb4zxGWMWYzXLjq/33FxjzJPGGC/wN//5e4hIf6zm+z8YY2qMMR9jJd0AN1aSHmSM8Rpj1hpjyhr5PZwCpAL3GWNqjTH/BpZSr5ka+Kcx5gtjjAd4ERjdzO92E+ABNhpjtjaxT6C6PgfYCuQ1sd/zwJUi0hkrSb7RzLmbZYzZ7W8G74bV/dEwxmeNMV8bY6qAVzn69d5pjDns366UagFN1qqlLjTGdDHGZBljfh744BWRaSKyyt/MXYJVdQealx/Eqizf9TeR3wZgjNkO/BK4E9gvIq+ISO9659oTuOGvMvdi9ZciIlf6B1eV+M83ot75oF6TrzGm0n8z1f/8Q8aYw/X2za13exHwDvCKiOSLyAP1qvn6egN7GlS/uUD9QVf1m50r/ecP5iGsL0B9ReRHTeyzCKuCvZog1bIx5hOslozbgaUNE6S/2b7C/+87jR6k6WMXY30B+leD/ufmXu8elFLHRJO1ajURSQD+Afwv0MNffb0FCIAxptwYM8cYcxzwfWB2oG/aGPOSMeZ0IAswWM3AAf3qncMB9AXyRSQLeBKrCT3Df75NgfM1owBIF5GUeo/1D9zwV/53GWOGAxOxmvav5Gj5QL9A33u94zRV6Qbl/31cgNXs/lPgzyLSteF+xphcrIFm5wH/bOawL2A1qR+V1I0xJxhjUv3/Vh5DyC6s5u5OLXiOLvGn1DHSZK3aQjxWc/YBwOMfzDUlsFFEzheRQf4+4zKs5m+viAwVkcn+ZF+N1RfqrXfck0Rkhr96+yVQg9Vvm4L1wX/Af/xr8A90a44/2a0B7hKReBE5HesLRCDW74rISH+fcBlWs7i3kUN9DhwGfi0icSIyyX+cV0KJoz7/F4cngV8aYw74B3m9BzzSxFOuBSY3aB1ozKNYzeUfhxhHnIgkYn0uuPyD7Zz+bTP8/18OEcnE6pZY76+ylVJhpslatZoxphy4Cauf8hBWM+2SersMBt7HGuH8H+BxY8xHWAn+PuAgVhNqd6zBZwH/whpEdQi4Apjhr3y3YDUZ/wfYB4wEPm1ByJdiDUArBu7gyMqzJ/AaVqLOwWqWPmoQlzGmFpgOTPPH/zhwZZC+5mDuBbYaY16s99gvgWkiMqXhzv6xAWuaO6gxptgY84EJfdH6J7G+MF0C/N5/+wr/tj7AcqyR4hsBH3BRI8dQSoWBhP4+VkoppVQkaGWtlFJK2Zwma6WUUsrmNFkrpZRSNqfJWimllLI5TdZKKaWUzUXd6jfdunUz2dnZkQ5DKaWiytq1aw8aYzJbeYzuLpfrKax5DbTYa1s+YJPH47nupJNO2t9wY9Ql6+zsbNasafYSU6WUUvWISG7zewXncrme6tmz57DMzMxDDodDr/ttQz6fTw4cODC8sLDwKaw5HI6g34yUUkqFakRmZmaZJuq253A4TGZmZilNzMaoyVoppVSoHJqow8f/u200L4ctWYvIMyKyX0Q2NbFdRORREdkuIhvqrX2slFJKqXrCWVk/B0wNsn0a1pzRg4HrgSfCGItSSqkOYP78+d0HDx58wqBBg06YN29ed4B9+/Y5J06cODgrK2vExIkTBx84cMDZ0uOOGTPm+LaPtu2ELVkbYz7GWiihKRcAzxvLKqCLiPQKVzxKKaWi2+rVqxOff/75zHXr1uXk5ORsXr58eZeNGzcm3HHHHb0mTZpUnpubu2nSpEnlc+fO7RnqMT0eDwDr168/lkV4Wsztdh/T8yI5GrwPRy5Gv9f/WEFkwlFKKRWqX732Zb+vC8uT2/KYQ3qmVT44c9SeprZv3LgxaezYsRVpaWk+gNNOO6188eLFXZYvX95lxYoVXwHMmjWr6MwzzxwK5M2ePbv3zp07EwoLC+MKCgrib7rppsI5c+YcXLp0adr8+fN7de/e3b1ly5bkHTt2bE5OTh5TWVm5funSpWl33XVX78zMTPeWLVuSzzvvvEMjR46sevzxx3vU1NTI66+/vuOEE06oyc/Pd11zzTVZeXl58QAPP/zw7ilTpjS6bO3s2bN7FxQUxO3evTu+a9eunjfffPOblv5uIpmspZHHGh24ICLXYzWV079//3DGpJTyeaG8ELw14K6G/PUw6CxIC7lYUSosRo8eXTVv3rw+hYWFzpSUFPPee+91HjVq1OGioiJXVlaWGyArK8tdXFxcl9tycnKS1q5dm1NeXu4cM2bM8B/84AelABs2bEhZv3795uOPP7624Xm2bt2a9Nprr+3s3r27Jysra2RCQsLBjRs35syfP7/7Qw891P2ZZ57ZM2vWrH6zZ8/ed+6551Zs27Yt/txzzx28c+fOzU3FvmHDhuTPP/98a2pq6jEN0Itkst4L9Kt3vy+Q39iOxpiFwEKAcePG6UhEpdqCzwueGoj3F0e1h+G9O+C/L4G7kQJh/PVw9p0Qn9L4sUr3QGJnSEoPZ9TKJoJVwOEyduzY6ptvvrlw8uTJQ5KTk33Dhw+vdLmCp7Fp06aVpKammtTUVM+pp55atnLlypT09HTviSeeeLixRA0wcuTIw4Hk379//5pp06aVAowaNapqxYoVaQCffvppp23btiUFnlNRUeE8dOiQIz093dfYMadOnVpyrIkaIpuslwA3isgrwASg1BijTeBKtYeSPbDwTCtBT5gF6dnw8UNQthfSekPPEeCphgFnwO5VsHc1fLHQ+verHZDSzTpO3lpYegsUbgLjtR7rNQq+9wj0Pcm676mBQ7tAHNBtcCRerepAbrnlloO33HLLQYAbb7yxT9++fWszMjI8ubm5cVlZWe7c3Ny4rl27egL7ixzZiBu4n5yc3GhSBUhISKhLqg6Hg8TERBO47fV6BcAYw5o1a3JCTcApKSlNni8UYUvWIvIyMAnoJiJ7gTuAOABjzF+Bt4DzgO1AJXBNuGJRSjWw7nmoLIIeI+DTP1uPpWfD1W9B9mmNP+fzBfD2r+HBgTBiJuz8CCoPWtuGfg9SMqC6DL5ZAU9Nht5jrSo8by24K639zvg1TP59uF+d6sDy8vJcffr08Wzbti1+2bJlXb744out33zzTcKCBQsy7r333sIFCxZkTJ06tSSw/9tvv93lnnvuKSgrK3OsWrUq7ZFHHsnbtGlTYmvjOP3008vuv//+7vPnz98H8NlnnyVNnDixqrXHbUrYkrUx5pJmthvghnCdXykVxNfLof9E+PHbVlVcXQL9TgFnkI+ECbMgYyC8+EPY9Jr12ElXw+mzIT3r2/1qymHlQ7D1LSjLgxEzoMdI65wfPwD9JsDgs8P68lTHNX369IElJSUul8tl/vSnP+3OzMz03nXXXQUXXXTRwKysrG69e/eufeONN3YE9h8zZszhs846a3B+fn78rbfeWpCdne1ui2S9cOHCPdddd13/IUOGDPd6vTJhwoTyiRMn7m7tcZsiVs6MHuPGjTM6N7hSreB1wz09YeJNcPYdx34MJHhyb8hTC49PgLgU+OlKkMbGmKpwEZG1xphxrTnGl19+uWvUqFEH2yqmcJs9e3bv1NRU77x58/ZFOpZQffnll91GjRqV3fBxnW5UqVhTvBN8HsgceuzHcMa1LFEDuOLhlJ/Dvo2wr9GJDZVSTYi6VbeUUq1Uutf62SUr+H7hMPxCq997/Ysw7b72P7+KKQ8//HCjVxiFy5///OeMJ554okf9x04++eSKRYsWtbp5XJO1UrGmwt8imNYj+H7hkJoJx02Crctg6h+1KVx1KDfffHPRzTffXBSOY2szuFKxJpCsU7pH5vwjZkLpbtj9H2vSFaVUszRZKxVryvdBfCokpEbm/MMvgIRO8Ow0eHgYVJU0/xylYpwma6VizeH9kBqhqhqsLwkzFkLX46CqGHI/i1wsSkUJTdZKxZrqUkjsEtkYhk6Da9+3bpeE7dJUpToMTdZKxZqackhIi3QUkNwVHHFQrrMMK/ux2/rWmqyVijV2SdYi1hzjlVEzx4aKIe21vnWo9NItpWJNTbk1wMsO4lOtxURU9Hnjhn7s39Km61nTfXglFz7W5GpeX331VfzUqVMHjx8/vmLdunWpw4YNq/zxj398cN68eX2Kiopczz333M7hw4fXXHbZZdm7d+9OSEpK8i1cuDB3woQJVcuWLUudM2dOf7AW8/jss8+2OhwOpk6dOqi0tNTp8Xhk7ty5+ZdffnkJQGB9a4Dbb7+9x6uvvpohIpx11lmljz/+eF7D2Hbt2hU3derUupVqtm3blpSTk7NxyJAhja7s1VKarJWKNTVlkGiTZJ2QCjUVkY5CRZE9e/YkLl68eOdJJ52Ue+KJJw578cUXM9asWbP1pZde6nLPPff06tOnT+2oUaMq33///R1LlixJu+qqqwZs3bp1y0MPPdTz0UcfzZ0yZcrh0tJSR2DVrWXLlm3v2rWrr6CgwDVhwoTjL7300hKH49tG51dffbXTsmXL0teuXbs1LS3Nt2/fPmdjcWVnZ7u3bt26BeCPf/xj5sqVK9PaKlGDJmulYosx9mkGB62so1mQCjic+vTpUzN+/PgqgCFDhlRNnjy5zOFwMHbs2Mq77767d15eXsI//vGP7QDTp08vv/76611FRUXOU045peLWW2/td/HFFxdfcsklhwYOHOirqamRX/7yl31XrVqV6nA42L9/f/zevXtd/fv3r1ti87333ut0+eWXH0xLS/MB9OjRwxssvnfffTfl+eefz1y1alWbNqNrn7VSscRdCcZno2SdArXlkY5CRZH4+PhG15p2Op14vd5GF6cSEXPvvfcWPvXUU7lVVVWOiRMnDlu/fn3iggULuhYVFbk2btyYs3Xr1i0ZGRnuqqqqI/KiMeaoNbGbkpubGzdr1qzsxYsX7+jcuXOr1q9uSJO1UrEkcJmUbZK1VtaqbZ1yyinlzz77bAbA0qVL09LT0z1du3b1bd68OWH8+PFV99xzT+HIkSMPb9q0KbG0tNTZrVs3d0JCgnnzzTfT8vPz4xseb+rUqWWLFi3qVl5e7gBoqhm8pqZGZsyYcdz8+fPzTjzxxJq2fl3aDK5ULHlmqvUzLiWycQTEJ0NtZaSjUB3I/fffn3/ppZdmDxkyZHhSUpLvueee+wbggQce6P7ZZ591cjgcZsiQIVUzZ84sLSkpcU6bNm3QiBEjhp1wwgmVAwYMOGr+25kzZ5atW7cuefTo0cPi4uLM2WefXfqXv/zlqAFm77//fsqmTZtS7r777t533313b4Dly5dvy87OdrfF69L1rJWKJXd2tn7+4GkYOTOysQAsuxU2vQa/2RXpSDq8WFzPOhrpetZKqW/FR2he8IZcCeBpswGzSnVY2gyuVCyKS4p0BBZnPHjbvHtPqbC64oor+q9evfqIb7w/+9nP9oVreUzQZK1UbLJLsnYlgM8DPh84tKFPRYdFixa1+4T2+u5QKhZlDIp0BBanf/CtVtdKBaXJWqlYMvKH1tKUyV0jHYnFlWj99GiyVioYTdZKxRJvrbXSlV24ApW1DjJTKhhN1krFEq/726ZnO3AmWD+1slYqKE3WSsUSrxucdqqs/claK2vVCuPHjx/68ccfH7UC2Jlnnjno4MGDR804Nnv27N5z587t0drz+lfaOq61xwmFjgZXKpZ4a21WWftj0cpahcGKFSu2h/P42dnZ7uXLl+8M5zkCNFkrFUvsVlkHkrWvTWZkVO3oD5/+od/2Q9vbdD3rQemDKuefNj/oetbnn3/+4G3btm0GmDt3bo+Kioq6ytnr9fLDH/4wu2/fvrWPPvpofp8+fUauWbMmp1evXp7f/OY3PRcvXtytd+/etRkZGe4xY8ZUAtx9993dn3322Uyn02mGDBlSvXTp0p2NrX2dnp5+1MIcDeMJJ03WSsUSby3E2WQtawCH/yPI6wm+n1LNcLvdcuGFFw4YPnx41f33319Yf9vKlSuTX3/99a4bN27c4na7GT169PBAsn700Ud75ubmbkxKSjKBJvOm1r6OJE3WSsUSn80GmDn8RZFPk3W0CVYBR8LPf/7zrAsvvLC4YaIG+PDDD1PPO++8ksCa1FOmTCkJbBs6dGjVRRddNGD69Okll112WQlAY2tft98raZwOMFMqltiuGdwfiyZrFQKXy2V8vm/zZnV1dV0OGzduXMXKlSs7VVZWNrr4dFNrUn/44YfbbrjhhgNr165NGTVq1HC3201ja1+39WtpKU3WSsUSu11nHWgG1z5rFYK+fft6iouLXYWFhc6qqip55513Oge2zZo16+CUKVNKzz///IFu95F/T5MnT65YtmxZl4qKCjl06JDjvffe6wJWH/eOHTviv//975c//vjje8vLy52lpaXOxta+bueXepSwJmsRmSoiX4nIdhG5rZHt/UXkQxFZLyIbROS8cMajVMyzW2Vdl6y1slbNS0hIMHPmzCkYP378sLPOOmvQoEGDjlh/+s4779w3atSoyhkzZgzwer11j59++umVF110UfGIESNOOP/88weOHz++AsDj8cill146YMiQIcNHjBgxfNasWfu6devmfeCBB7oPHjz4hKFDhw5PSkryzZw5s7SpmESkXdaZDluftYg4gceAc4C9wGoRWWKM2VJvt9uBV40xT4jIcOAtIDtcMSkV83xem1bW3uD7KeV3++2377/99tv3N7X9kUceyQ/czsvL2xi4ff/99xc21p+9du3arxo+9re//S2k/vj9+/e7Onfu3C5/vOGsrMcD240xO40xtcArwAUN9jFAYGhqZyAfpVT4+NzfDuqyg7rR4NoMrqLLxx9/nHzFFVccd+ONN+5rj/OFczR4H6D+t5O9wIQG+9wJvCsivwBSgLPDGI9SyufRZnClWuiLL75IuvLKKwfUfyw+Pt63a9euTe0VQziTdWND7xq27V8CPGeMeUhETgUWicgIY8wRw+RF5HrgeoD+/fuHJVilYoLX822CtANN1ioKjB8/vmrr1q1bmt8zfMLZDL4X6Ffvfl+Obua+FngVwBjzHyAR6NbwQMaYhcaYccaYcZmZmWEKV6kY4LNZsnZqslYqFOFM1quBwSIyQETigR8BSxrssxs4C0BEhmEl6wNhjEmp2OZz2ytZa2WtVEjClqyNMR7gRuAdIAdr1PdmEZknItP9u80BfiIiXwIvA1cbY9plGLxSMUkv3VIqKoX1K7Yx5i2sy7HqPza33u0twGnhjEEp5efzAcZmlbX/i4OOBlcqKJ3BTKlYEZglzFbJOjA3uF5nrUIzf/787oMHDz5h0KBBJ8ybN687wL59+5wTJ04cnJWVNWLixImDDxw40OLrE8eMGXN820fbdjRZKxUrAk3NtkrW2gyuQrd69erE559/PnPdunU5OTk5m5cvX95l48aNCXfccUevSZMmlefm5m6aNGlS+dy5c3uGekyPx/rbW79+/dawBd4GbPSuVUqFVaCp2ZZ91toMHm3yf/f7fjXbtrXpetYJgwdX9r73niZnD9u4cWPS2LFjKwKrZ5122mnlixcv7rJ8+fIuK1as+Apg1qxZRWeeeeZQIG/27Nm9d+7cmVBYWBhXUFAQf9NNNxXOmTPn4NKlS9Pmz5/fq3v37u4tW7Yk79ixY3NycvKYysrK9UuXLk276667emdmZrq3bNmSfN555x0aOXJk1eOPP96jpqZGXn/99R0nnHBCTX5+vuuaa67JysvLiwd4+OGHd0+ZMuVwo7+r/HzXzJkzB5SUlLhGjx5d+dFHH3Vau3ZtTq9evUL+lqqVtVKxItDUbKfKWlfdUi0wevToqs8//zytsLDQWV5e7njvvfc679mzJ76oqMiVlZXlBsjKynIXFxfX/ZHn5OQkvf/++9tWrVq19cEHH+y9a9euOIANGzakPPjgg3k7duzY3PA8W7duTXriiSf25OTkbH7ttdcyvv7668SNGzfmXHHFFQcfeuih7gCzZs3qN3v27H2bNm3Kef3113f89Kc/zW4q7ttuu633mWeeWb5ly5acGTNmHCooKGjxOrU2etcqpcLKjn3Won3W0SpYBRwuY8eOrb755psLJ0+ePCQ5Odk3fPjwSpcr+N/ztGnTSlJTU01qaqrn1FNPLVu5cmVKenq698QTTzx8/PHH1zb2nJEjRx4OJP/+/fvXTJs2rRRg1KhRVStWrEgD+PTTTztt27YtKfCciooK56FDhxzp6elHrX39xRdfpL7xxhvbAWbOnFnWqVOnFv/B2+hdq5QKq0D1aqtmcAeIQ0eDq5DdcsstB2+55ZaDADfeeGOfvn371mZkZHhyc3PjsrKy3Lm5uXFdu3ata6ppuI514H5ycvJRSTUgISGh7hJih8NBYmKiCdz2er0CYIxhzZo1Oampqc1ebtwWVyRrM7hSscJrw8oarMu3tBlchSgvL88FsG3btvhly5Z1ufbaa4vPPffckgULFmQALFiwIGPq1Kklgf3ffvvtLpWVlVJYWOhctWpV2umnn95ov3JLnX766WX3339/98D9zz77LKmpfcePH1+xaNGirgD//Oc/O5WVlbV4tLrN3rVKqbCp67O2UWUN1pcHTdYqRNOnTx9YUlLicrlc5k9/+tPuzMxM71133VVw0UUXDczKyurWu3fv2jfeeGNHYP8xY8YcPuusswbn5+fH33rrrQXZ2dnuTZs2JbY2joULF+657rrr+g8ZMmS41+uVCRMmlE+cOHF3Y/ved999+TNnzjxu+PDh6aeeempFZmamu0uXLi1qCpdomzBs3LhxZs2aNZEOQ6nosz8HHj8FZj4LI2ZEOppv/bE/jL4Upt0X6Ug6NBFZa4wZ15pjfPnll7tGjRp1sK1iCrfZs2f3Tk1N9c6bN69dlrFsSlVVlbhcLhMXF8f777+fcuONN2Y1tTDIl19+2W3UqFHZDR/XylqpWGHHPmuwFvPQS7dUB7Z9+/b4iy++eKDP5yMuLs4sWLBgV0uPoclaqVhh2z5rbQZX4fHwww83XOkxrP785z9nPPHEEz3qP3byySdXLFq0aHdOTk6rlti02btWKRU22metWs/n8/nE4XBEV/9pO7n55puLbr755qJjfb7P5xOg0VHqOhpcqVhRd511iweihpfDCV5N1lFi04EDBzr7k4pqQz6fTw4cONAZ2NTYdq2slYoVdu2zdrjA6KQo0cDj8VxXWFj4VGFh4Qi02GtrPmCTx+O5rrGNmqyVihXaZ61a6aSTTtoPTI90HLFIvxkpFSvs2mctTk3WSjVDk7VSscK2fdYunRtcqWZoslYqVti2z9qpyVqpZmiyVipWaJ+1UlFLk7VSscKO61mDJmulQqDJWqlYYcf1rEGbwZUKgSZrpWJFoBncjn3Wep21UkFpslYqVgSamm1XWWszuFLN0WStVKzQZK1U1ArpXSsiPYH+9fc3xnwWrqCUUmFg12Stk6Io1axm37Uici9wObAVCHQsGeC8MMallGprdu6z9jW60JBSyi+Ur9g/AIYYY6rDHYxSKoz00i2lolYofdbfhLifUsrObHvpliZrpZoTyru2HFgvIu8DNYEHjTGzwxaVUqrt+TxWYhSbLUXs0D5rpZoTSrJe7v+nlIpmXrf9qmrQhTyUCkGz71xjzNMi4gIG+R/abozRr8FKRRuf137LY4JOiqJUCJrtixaR7wDbgaeBZ4CvReS0UA4uIlNF5CsR2S4itzWxz8UiskVENovISy0JXinVAj63/ZbHBO2zVioEobSJPQKcZ4zZAiAiw4BFwLhgTxIRJ/AYcA6wF1gtIksCx/HvMxj4LXCaMeaQiHQ/tpehlGqWz2Ml898AAAAgAElEQVS/y7ZAr7NWKgShjPKOr59gjTE5QHwIzxuP1WS+0xhTC7wCXNBgn58AjxljDvmPvT+0sJVSLaZ91kpFrVCS9ToRWSAip/v/PQGsD+F5fYA99e7v9T9W3xBgiIh8KiKrRGRqYwcSketFZI2IrDlw4EAIp1ZKHcXOfdaarJUKKpRk/VNgB/Br4DfATmBWCM9r7PoQ0+C+CxgMTAIuAZ4SkS5HPcmYhcaYccaYcZmZmSGcWil1FO2zVipqhTIavBp4wP+vJfYC/erd7wvkN7LPKmOMG/hGRL7CSt6rW3gupVRz7NpnrclaqWY1WVmLyMv+n+tFZF3DfyEcezUwWEQGiEg88CNgSYN93gC+6z9PN6xm8Z3H8kKUUs2wbZ+1/9It07DhTSkVEOyd+yv/z5nHcmBjjEdEbgTeAZzAM8aYzSIyD1hjjFni3zZFRLZgLRLyK2NM0bGcTynVDJ/XpsnaH5PxWSPDlVJHafKda4zZ67+ZD1QbY4yIDASGAu+GcnBjzFvAWw0em1vvtgFm+/8ppcLJZ+PKGvzToWqyVqoxoQwwWwkkiUgvYAXwM6zJUZRS0cTOfdag/dZKBRFKsnYYYyqxlsr8izHm+8CJ4Q1LKdXmvB57VtZSr7JWSjUqpGQtIicDlwJL/Y9pW5VS0cZn02RdV1nrtdZKNSWUZD0buAtYZozZJCLHYTWNK6Wiie37rDVZK9WUUK6z/jfwbwAREWCfMebn4Q5MKdXGvG7ts1YqSoWy6tbzItJJRJKBzViTl+jobaWijW0v3dI+a6WaE0oz+EhjTBlwIdYlW32Bq8MZlFIqDGzbDK6VtVLNCWnVLRFxYa2Y9YZ/BS1feMNSSrU5uw8wM/qxolRTQnnnPgXsBjYBK0SkP1AR1qiUioDdRZWs+Ho/CXFOpo/qTWLckRc9fLbjIAkuJydlpUcowlby2vU6a20GV6o5oQwwewR4JHBfRPYAk8MZlFLt4d3NhXyQs5/iylrcXh+rvynmcK01Inn7/gp+d94wAJZ8mU+/9CQuffJzAL68Ywqdk2yY9Jpj1xnCtBlcqWY1maxF5BJjzMsiclMTuzwappiUCpuPvz5ARmo8ZVUeZr2wljing56dEklJcDHlhJ5ce/oArn9+DXmHqgAorXJz08tHLt/+1xU7+M3U4yMRfuv43PZcz1onRVGqWcEq60Bbny4graKeMYanP/mGu5fl1D02oFsKb/z8NDonH5nAOifHU+Ox+k+XbSioe3zqCT0RgRdX5XLzWYOPaia3PdtPN6rXWSvVlGALeTzu//mH9gtHqba1cW8pOYVlfPz1AZZuKGBcVjpjs9JxOoSfnjHwqEQNEO9y4Pb6WLahgN+9vpERfTox74IRjO7bhU+2H+TtTYV8su0gZw/vEYFX1Ap2nW5Uk7VSzWr2nesfUHYjkF1/f2PMjPCFpdSx25RXSnm1hxP7dub7f/kEAIfAnHOGcMN3B+FwSNDnxzuFWo+Ph9/7CoB7LxrJiX27ADCkRxoA+8trwvgKwsS2o8G1GVyp5oTyzl0CPA+8h16ypWysosbD/W9vZdGq3KO23feDE7l4XL+QjhPvcrD3UBW5RZXc/r1hdYkaICXB6T+Xu22Cbk+2vc5ak7VSzQnlnVtrjHk47JEo1QrVbi9XPv0563aXcNmE/gzolsLCj3cyoFsKc6YM5eTs0C+3inM66gaY9U1POmJbSrz1lqmoibImW2Ps32dtoux3qlQ7CiVZ/5+I3A68A9S1/RljNoQtKqVaoLLWwy9eWs+63SX83yVj+P6o3gD8+LQBGMDZTLN3Q/FOBx6fAaBbasIR2xwOISXeyeGaKKsCA/3Btqys9dItpZoTyjt3CHAdMI1vm8ENcEa4glKqOZW1HvaX1dCvazK/eGk9H361n/kXjqhL1ECzfdNNiXN9O7Ffw2QNkJroYtXOIrbklzG8d6djOke7CyRCWyZrXXVLqeaE8s69GMg2xkThiBrVkewuquSFz3NZl3uI9XtK8PqrX4B5F5zAFadktcl5EpzfJuvunY5O1ikJLjbnl3HeoyvZdd/32uScYefz97HbMVnrddZKNSuUd+4GII16TeBKtSevz3Df2zn87T+5YOD4XmlccUoWL36ei9truHB07zZL1AC7ig7X3U6OP/otcmz1eoQFEqGd+6w1WSvVpFCSdQawVUQ+58g+a710S4VdfkkVE+/7NwDTR/Xmd+cNo2fnRMCajWznwcPMGNsXa6n1tlHRTH90WXUUJhWvnZvB9TprpZoTyjv3nrBHoVQDn+04yLINBXy2owiAqydmc+f0E47YZ2xWOjsPHmZoz7Q2PXe5Pxn/64bTGt1++/eGcfMr/6Vbanybnjes6vqsbTjrmiZrpZoVSrL+DKg2xhgRGQgMxVrXWqk2Z4zh0Q+288j7X5OW4GJYr07ccs4QptcbOBYw/4IRXH5KFj06JbZpDIFk3atz48e9YHQfVu8q5q2NhW163rCq67O2YzO4f4yANoMr1aRQkvVK4AwR6QysANYDPwKuDGdgKvbklVRx71s5LNtQwIyxfbj3opFB599Oincyul+XJrcfqxP7duazHUWkpzRdOSfHu6isjaLk4vUna+2zVioqhZKsHcaYShH5MfAXY8x9IvLfcAemOq6KGg8v+QeHzTypL9f9bQ0b80rrtv922vFcf8ZxbdoP3RJ/veIkdhdVEldvVHhDSXFOqt0+vD7T4uu4I8LWl27ppChKNSekZC0iJwOXAtf7H7Nhx5eKBpW1Hi587FO2768A4MF3vqrblhzv5PHLxjJpaPdIhQdAp8Q4RvTpHHSf5HjrLVDl9pKaYMME2FA0JGutrJVqUijv3NnAXcAyY8wmETkOq2lcqRZ77MPtbN9fwbPXnMymvaU89N7XAHww50wGZqZGOLrQBZJ1Za0nOpJ1VDSDa2WtVFOa/ZQxxvwb+He9+zuBn4czKNUxvbE+j8c/2sGMsX347tDunDawG9UeL2cMzoyqRA2Q5L/+uqo2ShJMXWVtw2QtOsBMqeaEskTmIKzqOpsjl8icEr6wVEfz4db9zPn7l5wyIIN7LxoJWKtb/erc4yMc2bFJqausoyRZ11XWNmwF0MpaqWaF8s59DXgaeAHQd5NqsTfW5/Hr1zZwfM80nrxqXNAR3tEiKdqStfZZKxXVQnnn+owx/xf2SFSHtGZXMbNf/S/jB3TlsUvHRkf/bgiSo64Z3M7XWWtlrVRzmr425Vv/EpHrRSRTRDoF/oVycBGZKiJfich2EbktyH4zRcSIyLiQI1dR4fGPdtAtNYGnrzqZjEZWsIpW9QeYRQWvnecG14U8lGpOKGXOdf6ff6j3mAH6B3uSiDiBx4BzgL3AahFZYozZ0mC/NOAm4PNQg1bRoarWyyfbD3L5hCxSOkhFHaDN4G1IxBpkptdZK9WkUEaD9zvGY48HtvtHjyMirwAXAFsa7DcfeAC49RjPo2yqoLSKWo+PEX2iZM3nFkiOumRt4yUywYpLK2ulmhTSO1dEjgeGA3WTJRtjXmrmaX2APfXu7wUmNDjuGKCfMWapiDSZrEXkevwTsvTvH7SgVzYSWL0qLdGGTa+tlOiyknWNJ0qStZ2vswZN1ko1I5RLt24HpgDHA+8A5wKfAM0l68bmYDT1jusAHgGubi4GY8xCYCHAuHHjTDO7K5sILIiRlmjTaq4VXE7rz9vjjZI/Rzs3gwOIUweYKRVEKAPM/gf4LlBgjLkCGEVoFfleoH4Tel8gv979NGAE8JGI7AJOAZboILOOI5CsO8oI8PoC84a7fb4IRxIin40HmIE1yEwra6WaFEqyrjLGeAGPfzBYIXBcCM9bDQwWkQEiEo+1UteSwEZjTKkxppsxJtsYkw2sAqYbY9a0+FUoWwo0g3fqgM3gLkeUVdbeaOiz1spaqaaEkqzXi0gX4BlgDfAFsK65JxljPMCNWE3nOcCrxpjNIjJPRKa3ImYVQTU7v6HkH/8Iad/yaitBpHbAZnBnXbKOssrajtdZg/ZZK9WMoJ+iYq1ReKcxpgR4TETeAToZY5pN1gDGmLeAtxo8NreJfSeFFLFqd57iYgCc6ekU/P73VH35JZ1nzGh2CcuKDtwMLiLEOx3URltlbcfpRsHfDK6VtVJNCfrONcYYEVkKnOS/v71dolK2YNxuip55lgOPPgoipE2eTNX69da26mokKSno8ytqPCS4HMS7QmnAiT4up0RhZW3nZK2VtVJNCeWd+4WIjA21mlbRy/h8HP7kE9x5eXhLSjj08it49u8n7ZxzqNm5k/J3363b11ddjaOZZF1W7emQI8EDXA7B44uSytrO042C9SVCJ0VRqklNfpKKiMvf73w68BMR2QEcxrokyxhjxrZTjKodeMvL2XvDjVR+8UXdY/HHHUef3/2WtHPPxdTUcOiFF6j++mvKlryJqayE9PSgx6yoiZK1no9RnNOBO1oqaztPNwraZ61UM4J9kn4BjAUubKdYVIRUbd5MwW2/peabb+h5550kTxhPXK9eOBLr5sBBEhPJuO46yt56i7Ilb+Krqmr2uBXV7g45IUpAVCXrusrapiueabJWKqhgyVoAjDE72ikW1c5K/vFPil94gZqcHBydO9N/4QJSJk4M+pxAP3Uoybq8umNX1lafdZQ0g3vd4LTxQio6KYpSQQX7JM0UkdlNbTTGPByGeFQYGa+XQy++SNFzzyFOF+49e0g4/nh6/P73dJ7+fZydOzd7jEA/ta8yhMq6xkP/rsmtjtuu4pwO3NHSZ+11gzM+0lE0TUeDKxVUsGTtBFJpfNpQFUXc+/ZT/s5ySpe8SfWmTSSNHo3ExdH5wgvodv31SFzoTdWBZG2qQ6ysO/oAs2hpBvfW2re/GrQZXKlmBPskLTDGzGu3SFSbq9q0mUMvvUT5O+/gO3yYuKz+9P7f/6XT985r9hrpprSsGdxNWoduBnfgjppm8FqbV9aarJUKptk+axVdTG0tBjj00ksceOhhjNtN6tlnkXnTTSQOGdLq4zuSrWbtqvXr6TR1atB9K2u9HW4d6/rinRI9A8y8bptX1toMrlQwwT5Jz2q3KFSr+WprOfjEExQ99TS4rZG/KaedRu/7/ogrM7PNzhNoBi/+2/P0+O1vm47HZ/D4TIedEAWsytoTLQt52L4Z3Ame2khHoZRtNZmsjTHF7RmIOjbuwkJ85eXk/fo31OTkkDppEkmjRxHXt1+rmrub4uzaNaT9av0VZ4dO1g7RZvC24nCBab5rRalY1XHbKDu42r15HHzsMUpffx0AR1oafR55mLRzz0Uc4UuQIkL65ZdT+uabweMLJGtnx03WcU4HlbVR0s/q89i8stY+a6WC0WQdRYwxVHz4EZVr1lC6ZAnekhLSr7yCuB49STtrMvHZ2e0Sh8TFYfxN7U2p9cRAZe2MoulG7V5Zi84NrlQwmqyjgPF68ZaVsf+++yn917+Q+HgSR4yg59NPkTh0aLvHI/HxmNrg/Yt1ybqDV9aB12l7dk/WOsBMqaA0Wduc8fnY85OfcPiz/4DTScasWWTe8HMkPnIfvBIfBx4Pxudrssk9MEo6rkMn62iqrO0+GtylyVqpIDRZ25AxhsOffELRk09RvWkTvspK4gcNpM8DD5A4fHikw0PirC8Kxu1GEhqfwjImmsEdjuiaFCU+JdJRNE37rJUKSpO1zfiqq8m79VYq3v8AZ0YGadOmkjR6NF1mzmzzkd3HSuKtCs3U1kITybrG0/Era5dTR4O3GU3WSgWlydpGqr/+moLb/0D1xo10/9WvSL/ichwRbO5uSmB60mCDzALN4AkduLKOj6ZVt2zfDK591koFo8naBqo2buTQy69Qvnw5kpBA7/99kM7f+16kw2pSoL882CCzmGgGj6o+a7tX1k4wmqyVaoom6wgybjeFd99DyeLFAKSdczY9fvtb4nr3jnBkwYVWWVtJrEM3gzuirLJ22Lmy1mZwpYLRZB0hnuJi8m+7jcMfr6TzD2bQ9cqrSBza+rm724MjlMraa1VJHbmyjouq9aztPt2oJmulgtFk3c6MMZS+8S8OPPww3pISut/2GzKuvjrSYbVMCJV1LFxnHX1zg9u4GVy0z1qpYDRZtyPPgQMUzL2Dig8/JHHECHrffx8pEydGOqwWC62ytirOeJc9RrCHQ5x/iUxjjG1G6jfJ67Z3stYBZkoFpcm6HRifj0OvvMKBPz+Kqa6mx29vI/2KK8I6h3c4hdJn/W1l7WyXmCIhzmElaI/PEOeMhmStzeBKRStN1mFmvF52X3sdlatWkXzKKfSc+wcSjjsu0mG1io4Gt7j8TfweryHOzt9JjLF/M7gma6WC0mQdRr7KSr6Z+UNqd+6k24030u3nP4vaarq+llxnbfuKsxUCr83t85GEjbO1zwsYmydr/6VbxoDduxSUigBN1mF0cMFCanfupPuvf03Xq6/qEIkatLIOcAWawe0+Itzr/3+yezM4WF8snPqxpFRD+q4IA19NDYVz76D0X/+i03nTyPjxNZEOqU3VVdZBB5h1/GQd539ttr/Wui5Z27yyBv/EKPqxpFRD+q5oY56DByn4w1wqPvyQ9CuuIPOGn0c6pDZXV1mHMMAsroO0JjQm8Nrsn6z9/09RUVl7gMbnm1cqlmmybkOH//Mf8n55C97SUrr/5jdkXHN1pEMKi0Cy9jVTWcc5BYej4/Y/upzR1gxu48pa/JW1DjJTqlGarNtI9ZYt7P7J9cT17k3fJ54geeyYSIcUNqFeutWRJ0SBeqPB7T4xSjQk6/p91jHA6/NSXltOl8QukQ5FRYmwfpqKyFQR+UpEtovIbY1sny0iW0Rkg4h8ICJZ4YwnXGp27mT3rFm40tPJXvxKh07UENoAM7fXV9en21HF+yvrWo/dK+toaAaPncr6xZwXGb1oNN9Z/B0mLZ5Erbfp95FSAWH7NBURJ/AYMA0YDlwiIsMb7LYeGGeMORF4DXggXPGES82OHeReeRUY6P/cs7jS0yMdUthpZW1xOaKtsrZxsnb5+6k9NZGNo419vPdjznr1LDYc2ADAvsP7uO+L++q2F1UXMWPJDLYUbWHk30Yy5bUp+IzN/55URITz03Q8sN0Ys9MYUwu8AlxQfwdjzIfGmEr/3VVA3zDG0+Yq161j1//8CASynv8bCQMHRjqkdhHqaPCOvOIWfNtn7dY+69ZzJVo/O1CV+Y+v/8ENH9zA/qr9XPbWZRRUFHD2a2fXbd9w5QZ+N+F35Jbl8j9L/weAgsMFjHp+VKRCVjYWzk/TPsCeevf3+h9ryrXA241tEJHrRWSNiKw5cOBAG4Z47DyHDlFw+x9wdEoj++WXo35WspYQpxOczmYr64QO3gweVzeDmc0roWhoBg98kfBURzaOY/THz//Iz97/GRsPbATgP/n/Yf6q+UfsM+UfU+puzzlpDiLCJcdf0ujxCg8Xhi9YFZXCOcCssWHAjZYgInI5MA44s7HtxpiFwEKAcePGRbSM8ZaUcOjllyl5/Q08hYX0W/BX4vtGVYNAm5D4eExtM83gHTxZu+rNDW5rgWrVzutZByrrKEzWZy4+k+LqYgA+yfuEqdlTWbF3BQM6D+D5ac9z6bJL2VW264jnXD3i6rrbG6/ayENrHuLlrS9T47W6AbaXbKdnSs/2egkqCoTz03Qv0K/e/b5AfsOdRORs4PfAdGOMrTusDq/6nB3nf58Dj/4fzrQ0+j/7DCmnnhrpsCJC4uKaH2DWwZvBAwPoau1eWQf6geOSIhtHMHV91tHVDP76ttfrEvWPR/yYM/qewaqCVYzrMY6F5ywkLT6NNy9684jnLLto2VHHmTNuDmsuX1N3/1D1ofAGrqJOOCvr1cBgERkA5AE/Ai6tv4OIjAEWAFONMfvDGEur1e7dy55Zs4jr15f+Ty4kcdiwSIcUURIfH7wZ3NvxK+vApCi2v87a7R8WEhXJ2t6VtdvrptJTSeeEzhRVFTH3s7kAfPDDD+ie3D2kY6TGpza57Y/f+SO/XflbTdbqKGH7NDXGeIAbgXeAHOBVY8xmEZknItP9uz0IpAJ/F5H/isiScMXTGu59+8mbPQdE6P/00zGfqKH5yjomRoPXTYpi88raXWX9jEuObBzB2Hw0eGC09tgXxnL6K6cz8m8jmfTqpLrtzSXqB898sO62y9F0jTQtexoAT296unUBqw4nrJOiGGPeAt5q8NjcerfPPupJNlO5Zg17f3kLvspKej9wP3E9ekQ6JFuQ+LhmKmtDUnzHTtbfrrqllXWrOf3J2mu/ZP1izotHXG7V0JdXftnsMSb3m1x3Oy7I2AGn/3rzQNO6UgEd+9O0FYwxFD+/iNyrr8GZksKAxa/QacqU5p8YI7SyrnedddRU1jZO1nUDzOyVrI0xRyXqy4ddXnd7wdkLcEjzf+fx9S6bC1ZZK9UU/atphPF6KfjDXEr/+U9SJ0+m9/334UxLi3RYttJsn7XH2/Ev3YqWVbfqKmttBm+phpdQbbzKujTrF2N+wZcHvuTU3i0fYOqS0D52fcYX0hcBFRv0L6EB4/WSf9tvKf3nP8n42U/p+5f/00TdiOYq62q3j8Q4ZztG1P7iHFEyKYq7CsRh80lR7DnA7MmNT9bdnnPSnLrbyXHJx5SoAUSCL27zq3G/AqC8tvyYjq86Jq2s66lct479Dz1M1dq1ZP7yZrr99KeRDsm2HHHBK+vKWg/J8R07WbuiZVIUT7VVVTeTJCLKppV1lwRroY2PLv6IjKSMdjlneqI1ZXFJTQmdEzq3yzmV/Wll7Vf+0UfkXnY57t276Tl/nibqZliTojRdWVe5vSR18GQdFy3Tjbor7d1fDfWmG7VXsg5U1l0Tu7b6WE4J7f0QSNAFhwtafU7VcWhl7Vf87HO4evVk4Jtv4khJiXQ4thesGdznM7HRDO6vrN12X8jDXQUumyfruulG7ZWsA5prug7Fmxe9SW5ZbrP7Bar55d8s55Rep7T6vKpj0GQNVG3cROXnn5N5yy2aqEMUbIBZjcdKXh29GbwuWdt9icxoqKxFrMu3bNZnnRqXyqR+k9rkWP3S+tEvrV+z+43oNgKAzOTMNjmv6hi0GRwoevJJnF26kH7Zpc3vrIDglXVlrbUmcVIHr6ydDsHpEGq93kiHEpy7yv7JGqx+axtNN1rlqaLCXcFxndt3kZ7ACPD1+9e363mVvcV8svbV1FDxySekTZuKM7XpaQDVkYJV1lVuK3l19GQNEO90REGfdZW9L9sKcNmrst53eB9AxBbUWF24OiLnVfYU883gJX9/DVNZSdrksyIdSlSRuDh87saroOpAsu7gzeBgDTKr9di8z7q2ApLbZyRzq7gS7ZWsKyOXrMd0H4MxNv8SqNpVzFfWxc8+S/K4caScflqkQ4kqEh8PTSyRWVkbQ5W1y2H/SVGqSyGxS6SjaF58qvXFwib2V1prC4W6QEdb6p7cnZKaknY/r7KvmE7Wxhg8Bw6QeOKJbTLaM5ZYlXUTzeC1sVRZR0uyjoLrdRNSocY+E4EEKuvMpPYf6NUloQuHanTlLfWtmE7WvsOVmNpaXF3TIx1K1AmpzzoGknW8y2HvZnBjoihZp0GNvSrrtLg0kiPQ35+emE5ZTRken6fdz63sKaaTtbfE+ubqTG/9hAexRuLiwO3GNHKNcVUMNYPH2X2AWe1h8HmiI1nbrBn8YNVBuiV3i8i50xPSMRjKassicn5lP7GdrIutZeicWlm3mMRbk1g0Vl3H0mjwOKeDGo8Xn12XyawutX5GQ7JO6GSrZvCiqiIyEiMzMC8w5eiham0KV5aYTtYef7J2pWuybimJs9bkDZasO/qkKADxTuH9nP2ccMc7kQ6lcVGVrFNt1QxeXF3cbvOBNxSYxUyTtQqI6WTtLfY3g3fVZvCWqqusG5kYJdAMnhgLydq/TGbgC4rtBJJ1UrSMBi+3+tltIJKVdWAuch0RrgJiO1kfCjSDa7JuqaCVdYz1WdtaVFXWaWB8366/HUE13hrK3eVtsoDHsQhU1sXVxRE5v7Ifm3/ShJenuBiJi9P5wI9B0Mra7SXOKfZPZG3A9k39dck6CirrBP8MgjZoCg80P0eqGbz+MplKQYwna2/xIZxdu+o11sdA4puurCtrvR1+xa2Afl1tPo1nZZH1MykKxmUkdLJ+1kR+BHRRlfV7i1QzeLwznpS4FO2zVnViPFkXaxP4MQo0g9ds237Utmq3NyaawAGyM2zeKlOWZ03jGQ3JOjAl6uGDkY0DKKwsBCK78lWCM4Gc4pyInV/ZS0wna8+hQzoS/Bg5kqxVnPY/8MBR26rcXvs3D7eRrAybV9ZledCpj7UEpd2l9rB+Ht4f2TiAvPI8APqk9olYDMXVxazdtzZi51f2EtPJWivrY5d88skAxPU5+sNMm8FtpCwfOvWOdBShSfXPwV0R+WSdfzifZFdy3UCvSKq0wYA7FXmxnawPHdIJUY6RIzGR5HHjoJEZzPaXVZOZlhCBqNpfv3SbJ+vSPOjcN9JRhCY5A8QBFfsiHQn5Ffn0Tu1ti/Es/9z2z0iHoGwgZpO1r7YW3+HDuLSyPmaSnISvquqox/ccqqKv3ZNYG4l3Ofj5pIG4HJH/UD+KzwvlBdFTWTuckJJpi8q64HABvVJ6RTSGhecsBCDBFRtffFVwMZusvUXWaE+dF/zYOZKSj0rWFTUeig/X0t/uzcNtyOV04PEZ+60/XLoXjDd6Kmuw+q3L8iMdRV1lHUkDOg8AYN5/5kU0DmUPMZusPQetEaeuzMhM1N8ROJKS8FUfmaz3FFv9a/26JkUipIiI81fVHrvND56/zvrZa1Rk42iJjEFQtC2iIew7vI+y2rK6ZBkpdugvV/YRu8n6wAEAXJmRuzQj2jmSkzCVjSfrWKusATx2W31r7xpwJkCPkZGOJHTdhsChXHBXRyyEXWW7ABjYZWDEYgBIdCXW3XZ7G1+OVsWO2E3W+zVZtwSQOEsAAAxRSURBVJYkHdln/c7mQq5fZF1qYvuBV20ozmlV1u5GBttFVN5a6HUiuOIjHUnoMocABg5+HbEQNhdtBiC7U3bEYgi4buR1ALyb+26EI1GRFrvJ+sABEMGVEZkZijoCR1IypqYG47XmAv/7mr1127okx0UqrHYXGFxmq8q69jDkr4e+J0c6kpbpM876ufs/EQvhkbWPANAzpWfEYggY030MALetvC3CkahIi+lk7UxPr5uJS7VcYGIUX1U1Pp9h3e5DnDO8Bx/MOdMWl7y0l2+bwW1UWW9ZAp5qOP78SEfSMulZ0PU4+OqtiJx+32HrsrHhGcMjcv6Gzuh7Rt3tak/kugZU5IU1WYvIVBH5SkS2i8hRXw1FJEFEFvu3fy4i2eGMpz53YYE2gbeSI8Vq6vaVl7Hi6wMUH67l3BN6MjAzNcKRta9vm8FtUllXFsP7d0KPEdD/1EhH03IjL4adK6xm/Hb28NqHAXjwjAfb/dxNGZo+FICTX4yyVhLVpsKWrEXECTwGTAOGA5eISMOvq9cCh4wxg4BHgPvDFU99ngMHqPz8C2tSD9VixhgOVtTw3xTr0pY//uFJrnluNT07JfLdobH3BSgwW9sHOZGfzIMDX8Piy6HyIFz4ODiisPFswizrcrOXL4VN/2h04p22Vu2p5umNT/PWN2+R3Smb/p36h/2coXpm6jN1t1fuXRnBSFQkucJ47PHAdmPMTgAReQW4ANhSb58LgDv9t18D/iIiYsJwwWrZu+9S8spianfvxp1nzfvb+cIL2/o0HdKCFTv44ptiSqrclFTWUlBaTWWtF4zhkfT+XLDqH4y5LovJP7uMlIRw/knZ09nDenDaoAzm/mszH27dz53TTyCrPRf4qD0My+bA/hwo+K+1cMeFf42uS7bqS+4Kl/0dXvgBvPZjiLvRei0n/g+cdHVY5jmv8lTxp3V/AuCBM46e7z6SOsV3IiMxg6LqItbtX8d3+n4n0iGpCAjnJ2sfYE+9+3uBCU3tY4zxiEgpkAEcseyOiFwPXA/Qv/+xfeP1VVbiLS8nafRouvxgBimnnUbSyBHHdKxYU1BaTWFZNV2S4xjaM40zhmTSv2syQ3qkMermcRTfMZfBZ55IcgwmaoCUBBfP/3gCT3y0nRc/393+X1hcSbDnC+jcB86+E0ZfDqlR3sLRfRjc8IXVnF95EPZtgV2fWMk6DNIT01l8/mIGdB5Akst+cwR8ePGHlNSU1K1zrWKPhGvWJRH5IXCuMeY6//0rgPHGmF/U22ezf5+9/vs7/PsUNXXccePGmTVr1oQlZqVayxgTU4Pr2o3PBz5PdF2GZjMistYYo31/USqcHVp7gX717vcFGs4jWLePiLiAzkBxGGNSKqw0UYeJw6GJWsW0cCbr1cBgERkgIvHAj4AlDfZZAlzlvz0T+Hc4+quVUkqpaBa2zjV/H/SNwDuAE3jGGLNZROYBa4wxS4CngUUish2rov5RuOJRSimlolVYR8IYY94C3mrw2Nx6t6uBH4YzBqWUUiraReFFmEoppVRs0WStlFJK2Zwma6WUUsrmNFkrpZRSNhe2SVHCRUQOALlhOHRnoDQMx7WDaHttdo030nFF4vztdc5uNJi5ULW5LGNMlE9tF7uiLlmHi4gsNMZcH+k4wiHaXptd4410XJE4f3udU0TW6OxaSjVNm8G/9WakAwijaHttdo030nFF4vyRfs1KKbSyVkrZgFbWSgWnlbVSyg4WRjoApexMK2ullFLK5rSyVkoppWxOk7VSSillc5qsW0FEhonIX0XkNRH5WaTjiWX6f6GU6sg6TLIWEaeIrBeRpa04xjMisl9ENjWybaqIfCUi20XkNgBjTI4x5qfAxYCOZAVEpIs/YW4VkRwROfUYj6P/FzFORFJEZK2InB/pWJSKtA6TrIGbgZzGNohIdxFJa/DYoEZ2fQ6Y2sjzncBjwDRgOHCJiAz3b5sOfAJ80JrgO5A/A8uNMccDo2jwf6L/F7GrqS9gjX358vsN8Gr7RqmUPXWIZC0ifYHvAU81scuZwL9EJNG//0+ARxvuZIz5GChu5Pnjge3GmJ3GmFrgFeAC/3OWGGMmApe1+oVEORHpBJwBPA1gjKk1xpQ02E3/L2LXczT4AtbUly8RORvYAuxr7yCVsiNXpANoI38Cfg2kNbbRGPN3ERkAvCIifwd+DJzTguP3AfbUu78XmCAik4AZQALw1jHE3dEcx/+3d6+xdhUFFMf/yxKgtIRasfEBcpUQ5BVagcQWTBD4JJREQSARI48EEwJqQBOaSDCKRsRHwBCjII9EA5qmFkyIlACloUIp5VFaHx+IJFRUXqUUqCLX5YeZ02x2jr33tLV3n9P1+3TOnJk9s3uSu/bMntMNLwK3SDoaWAN82fYbvQr5LnZftldIGmsVb734ApDUu/iaCcygBPgWSXfb/s8uHG5Epwx9WNf7WS/YXlP/YPdl+3v1D8FPgINtvz5IN/0P6eXA8gGOM+r2AD4GXGp7laTrgCuAK5uV8l1EQ9+LL9uXAEg6D3gpQR27u1FYBj8eOF3Ss5Ql0ZMk/aJdSdIngCOB3wBXDdjHBuDAxvsDgOe3a7SjbQOwwfaq+n4xJbzfId9FNPS9+Nr6wr7V9nZvGo0YFUMf1rYX2T7A9hhwDnC/7XObdSTNA26kLK+dD8yWdPUA3awGDpH0YUl71n7u2iknMEJs/x14TtKhtehkyn3HrfJdREsuviImYejDepL2AT5r+5m6nPYF+jwTW9LtwMPAoZI2SLoQwPbbwCXAPZTdzb+2vX6XjX64XAr8UtJaYC7wndbn+S6iKRdfEZOQ/xs8InaJegF2IrA/ZZf3VbZ/LulTlE2i04CbbX976kYZ0U0J64iIiI7bXZbBIyIihlbCOiIiouMS1hERER2XsI6IiOi4hHVERETHJawjIiI6LmEdERHRcQnriIiIjktYR4wIST+W9Lik46Z6LBGxcyWsI0aApBnAHOCLwGlTPJyI2MkS1jF0JP1I0lca7++RdFPj/Q8kXbaT+xzkmduTOd4sSRc33o9JWjfJttMlPShpWq/M9hvA+ynP9L5e0p6SVkga+mfWR0TCOobT74EFAJLeRXkwxBGNzxcAK6dgXIOYBVw8Ya3+LgCW2B7vFUh6D+WJZpuBcdtvAfcBZ+/oQCNi6iWsYxitpIY1JaTXAZslvVvSXsBhwBOSlkpaI2m9pIt6jSVd05rVfkPS5ZLOlfSopCcl/bQ5c23U7Vunzoz/KOnG2t8ySdPrZ1dK+pOkeyXdLumrwHeBg+txrq2Hn9avfR+fA+5slX0d+D6wHji8li2tdSNiyCWsY+jYfh54W9KHKKH9MLAKmA8cC6ytM8sLbB9Ty75UZ58Ad/DOGedZwGO17Hjbc4FxWkEn6bAJ6hwC3GD7COBV4AxJxwJnAPOAz9SxAFwBPGN7ru2v/a/27XOvz3z+iO1nG2Vj9d/hV5RnfPdWGdYB2WwWMQJyPyuGVW92vQD4IfDB+noTZZkcSkB/ur4+kBKGL9t+QtIcSR8A3gtsBI4CjgFWSwKYDrzQ6vPkCer8xfaT9fUaYIyyRH+n7S0Akn67jXPq175tf0qQN10NfNO2JW0Na9vjkt6StK/tzdvoNyI6LmEdw6p33/ooygzyOeBy4DXgZkknAqcA822/KWk5sHej/WLgTOB9lJm2gNtsL9pGnxPV+Vfj9TglzDXAOfVr37aFxnlImkuZsZ8g6Yb62dON+nsB/xxgDBHRQVkGj2G1kvITpVdsj9t+hbJpaz5lWXw/YGMN6o8CH2+1vwM4hxLYiymbsc6UNAdA0mxJB7XaTKZO20PAQkl7S5oJnFrLNwP7DnrStjdS7m33AvsaYKHtMdtjwNHUmXVd9n/R9r8H7SciuiVhHcPqacqS8COtsk22XwJ+B+whaS3wrVY9bK+nhOVfbf/N9h8om7SW1Tb3Un4K1WwzYZ0226uBu4CngCWUe+ObbL8MrJS0rrHBbLKWUWbSJwEzbN/X6O8fwAxJs4FPAncPeOyI6CDZnuoxRIw0STNtvy5pH2AFcJHtx3fgePOAy2x/foJ6S4BFtv+8vX1FRDfknnXE/9/PJB1OuZ98244ENUDdIPeApGnN31o31V3jSxPUEaMhM+uIiIiOyz3riIiIjktYR0REdFzCOiIiouMS1hERER2XsI6IiOi4hHVERETHJawjIiI6LmEdERHRcf8FQqBcejdZXTgAAAAASUVORK5CYII=\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": {},
   "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": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ukidss_j: mean flux error: 6.260478496551514, 3sigma in AB mag (Aperture): 20.715678042766548\n",
      "90prime_g: mean flux error: 0.1966152936220169, 3sigma in AB mag (Aperture): 24.473153622729576\n",
      "90prime_r: mean flux error: 0.20376788079738617, 3sigma in AB mag (Aperture): 24.43435754125472\n",
      "mosaic_z: mean flux error: 0.880610466003418, 3sigma in AB mag (Aperture): 22.84523725640468\n",
      "ukidss_j: mean flux error: 12.898414611816406, 3sigma in AB mag (Total): 19.930856030796242\n",
      "90prime_g: mean flux error: inf, 3sigma in AB mag (Total): -inf\n",
      "90prime_r: mean flux error: 0.5266921520233154, 3sigma in AB mag (Total): 23.4033047451995\n",
      "mosaic_z: mean flux error: inf, 3sigma in AB mag (Total): -inf\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": 17,
   "metadata": {
    "collapsed": 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": 18,
   "metadata": {
    "collapsed": 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": 19,
   "metadata": {},
   "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": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Depths on XMM-13hr')"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfoAAAEgCAYAAABcsLEAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XlcVXX+P/DX+7IjiKIoighuqCgqZlhqo2OlaGpurWr9tHKbJk3NLCebNFvGtHLK0tbRsW+WWy5l6oxbOpXgyipqkAkoiiAIIsv798e92A1ZlctyeD0fjx7de87nnM/7Xnr0Pp/P+dzzFlUFERERGZOpugMgIiIi22GiJyIiMjAmeiIiIgNjoiciIjIwJnoiIiIDY6InIiIyMCZ6qrVEJF5E7qnuOIyE3ymR8TDRU6WwJIhsEckQkTQROSAik0WkUv4bE5HPReTVyjiXrYjIUBFJFhFPq233i8hZEfGwvFcROSci9lZt7EXkvIio1bbdlrZdi/Sx0bK9XylxrBCRWBEpEJH/V2Tfw5Z96ZY+/yUi9W/90xNRTcVET5VpqKq6A/AD8AaA5wF8Ur0hVR1V3QzgvwDeBgARaQDgAwBTVDXdqmkagEFW7wcDuFTMKU8AeKzwjYg0AnAHgJQyQjkKYCqAQ8Xs2w+gt6p6AGgNwB7ATV1AWV+sEFHNxURPlU5V01V1E4CHADwuIp0BQEScROQtEfnVMqr9UERcLPv6ichvIvKiiFywzBCMseybCGAMgNkikikim6266yYixywj1DUi4mw5prGIbLHMLqSKyL6SZhdEpJeIHLSc46CI9LLat1tEFojIfstsxXYRaVzKx38GwCARGQhzwt9j+S6srYJVAre8XlnMuVYDeEhE7CzvHwGwAcC1UvqHqr6vqv8BcLWYfWdU9YLVpnwAbYs0K+k7LfwbPS8iyQA+Ky0OIqoZmOjJZlT1ZwC/AbjLsulNAAEAusGcXHwAzLM6xBtAY8v2xwGsEJH2qroC5qT3D1V1U9WhVsc8CCAUQCsAXQD8P8v2mZa+vQA0BfAigBue92yZZt8KYCmARgCWANhqGT0XehTAeABNADgCmFXKZ74AYJol3iEwJ/6iNgL4k4g0sIz67wLwTTHtEgFEARhgeV/SBUGFiEgfEUkHkAFgFIB3ijQp6TsFzH8jT5hnbSbeaixEZHtM9GRriQA8RUQAPAXgWVVNVdUMAK8BeLhI+5dUNUdV98CcgB8s4/xLVTVRVVMBbIb5IgIAcgE0A+Cnqrmquk+LL+xwH4A4VV2lqnmq+n8AYgBYX0x8pqonVDUbwFdWfZTkRwAeALaranHT7FctsT4E8+ffhGJG3xYrATwmIu0BNFDV/5XRd5lU9QfL1H0LAIsAxBdpUtJ3CgAFAF62/I2ybzUWIrI9JnqyNR8AqTCPrF0BhFum09MAbLNsL3RJVa9YvU8A0LyM8ydbvc4C4GZ5vQjASQDbReS0iMwp4fjmln6sJVjiLquPkqyAOUEPtr4NUMRKmEfoZY3S1wPoD+CvME/5/4HlVkbhPy3LiOsPVPUszH+DL4vsKu3zpqhqSRclRFQDcTEN2YyI3A5zwvwBwAUA2QA6WRJMcRqKSD2rZN8SQITldYXKLFpmDGYCmCkinQDsEpGDlnvX1hJhnoa21hLmBFhhIvIEAF+Yp+3DAHwkIsGqWvS++j6YZxwU5u+nTQmfI0tEvgMwpbg2qlrWRUdZ7EvquwQsd0lUy3BET5VOROqLyBCYR4r/VtXjqloA4CMAb4tIE0s7H8uiNWuviIijiNwFc7L82rL9HMyrxMsbwxARaWu5ZXAZ5kVn+cU0/RZAgIg8avmZ20MAAgFsKf8nvt5nc5hnEp5S1RwAHwK4CGBu0baW2whDAQwr4ZaCtRcB9FXV+HLG4WhZQCcAHETEuXAhooiMEZGWYuYHYCGAohc/RGQgTPRUmTaLSAaAMzAntyUwL2Ir9DzM0+k/ishlADsBtLfanwzzz8wSYV7MNllVYyz7PgEQaJn231iOWNpZzp8J4H8Alqnq7qKNVPUizBcUM2FOyrMBDCmyMr28lgH4UlX3Wc6tMK9LmG6ZVSjad6SqRpZ1Usv98h8qEMd2mGdPesF8GyEbwJ8s+wIBHID5e9kPINYSIxEZlJQ9mCCyPcsDYP6tqi2qOxYiIiPhiJ6IiMjAmOiJiIgMjFP3REREBsYRPRERkYEx0RMRERlYrXhgTuPGjdXf37+6wyAiqlXCw8MvqKpX2S1LPL6Jvb39xwA6gwPDmqoAQEReXt6Tt9122/niGtSKRO/v74+wsLDqDoOIqFYRkaKPd64Qe3v7j729vTt6eXldMplMXNBVAxUUFEhKSkpgcnLyxwCGFdeGV2hERFSSzl5eXpeZ5Gsuk8mkXl5e6TDPuhTfxladi4iviOwSkWgRiRSRaUX2zxIRLaO2NxERVR8Tk3zNZ/kblZjPbTl1nwdgpqoeEhF3mKuW7VDVKBHxBXAvgF9t2D8REVGdZ7MRvaomqeohy+sMANH4vfTn2zA/U5xXikREVC2WLl3a6LHHHruhvPM//vEPr/fee69R0e2xsbGO7dq1u6Fuxc146KGH/MLDw50r41xlqZLFeCLiDyAYwE8iMgzAWVU9ai4sVuIxEwFMBICWLStUZpuIiOimzZ49O8XWfaxZs+aWFkpWhM0TvYi4AVgHYDrM0/lzAQwo6zhVXQFz5S306NGDI38iomr03NqjvieSM1wr85wB3u5Zi0Z3PVNWu3vuuadNUlKSY05Ojmny5MnnZs2adcHV1TV4zJgxKfv373f38PDIX7du3enmzZvnFXd8SEhI+7feeuvMn/70p6ykpCT7Hj16dDx79uxx6zZffvmlxxtvvNHsu+++O7lo0aImbm5u+fPnzz+3b98+1yeffNLfxcWloGfPnpmF7cPCwpzHjx/fKjc3VwoKCrBu3bpTfn5+ucOGDWudlJTkWFBQILNnz0586qmnLpUVU0W/t4qyaaIXEQeYk/xqVV0vIkEAWgEoHM23AHBIREJUNbnSA/huDpB8vOx2dZV3EDDojeqOgoioVKtXr45v2rRpfmZmpgQHBweOHTv2UnZ2tql79+5ZH3300W+zZs1qNmfOnOYrV668qXVfK1eubPDuu+823bFjR5yXl1e+9b4nnnjC/+233/71vvvuy5w0adL16pr//Oc/vaZOnXpuypQpqVevXpW8vDysXbvWw9vbO3f37t0nAeDixYt2t/bJK4fNEr2YM/knAKJVdQkAqOpxAE2s2sQD6HGTtb+JiKiKlGfkbStvvvlm061btzYAgOTkZIfIyEhnk8mEJ598MhUAJkyYcHHkyJFtb+bcBw4ccD969Kjrrl27Tnh6ehZY77t48aJdRkaG3X333ZdZ2M9///tfDwC48847r7z11lvNfvvtN8eHH374UlBQUE737t2z586d6ztlyhSf+++/Pz00NDSzuD6rmi1H9L0BjANwXESOWLa9qKrf2rDPP+Jo9Za9+fObiEmNqe4wqlwHzw54PuT56g6DqM7bsmWL+549e9zDwsJi3N3dC0JCQtpnZ2ffsJC8tDVf9vb2mp9vHqhnZWX9oWHLli1zfv31V6eIiAjnotPoqlrieSdPnpx61113XdmwYYPHoEGDApYtWxY/bNiwjEOHDkWtW7fOY+7cuT47d+68/NZbbyXdxMeuVLZcdf+DqoqqdlHVbpZ/vi3Sxp+jeSIiKklaWpqdh4dHvru7e8Hhw4edjx49Wg8ACgoK8NlnnzUEgM8//7xRSEhIRknn8PX1zfn555/rAcDq1asbWu9r0aLFtXXr1p0cP358q7CwsD+sgm/cuHG+m5tb/vfff+9m6cezcF9UVJRjx44dc/72t7+dHzBgQNqRI0dc4uPjHdzd3QumTp2aOn369HNHjhyp1DUNN6tWPAKXqg9HtURUnUaNGpW+YsUKr4CAgMA2bdpc7dq16xUAcHFxKYiMjHTp1KmTt7u7e/769etPl3SOOXPmnHvooYdaf/nll43uuuuuy0X3d+3aNWflypWnH3rooTabNm06ab3vk08+iS9cjNe/f//rx65atcrz66+/bmRvb69eXl65r7/+euIPP/xQ74UXXmhhMplgb2+vy5YtK3VlfWmzEJWpVtSj79Gjh/JZ90REFSMi4ara42aPP3r0aHzXrl1r5Kyrq6trcFZW1uHqjuNmBQQEBG7atOlkhw4drlXG+Y4ePdq4a9eu/sXt47PuiYiIqlCvXr3atW/fPruyknxZOHVPRES1TnGj+XHjxrU8ePCgm/W2KVOmnJs2bdrFqovsj+699942Z86ccbLetnDhwt9GjRp1wy0EW2GiJyIiQ1i1alWNq5+yY8eOU9UdA6fuiYiIDIyJnoiIyMCY6ImIiAyMiZ6IiMjAmOiJiMjQ/ve//7l069atQ0BAQGD//v3bpqamXs99L7zwgnfLli07+/v7d163bl39ip57+vTpzTdu3OheuRFXLiZ6IiIytKeeesp/4cKFv504cSJq2LBhl1555RVvAAgPD3dev369Z2xsbOS2bdtOTJ8+vWVeXrGVbouVl5eHd955J3H48OElPn63JuDP64iIqGwb/+KL81GV++z2JoFZGP6+zevRx8fHOw8aNCgTAIYMGXJ54MCBAe+++27i2rVrG4wcOTLVxcVFO3TocM3Pzy9n9+7d9Xx9fXNDQ0PbBQcHX4mIiHBt3br11a+//jre3d29wMfHJ+iRRx65sGvXrvqTJk06//3333sMGTIkffz48Zd8fHyCRowYkfrDDz+45+XlyYcffpgwZ84cn4SEBKe//vWv52bPnp0CAC+99FLTDRs2eF67dk3uu+++tLfffjuxpM/+3HPPNVu7dq1ns2bNrjVq1CgvODg4a/78+ecq8jVzRE9ERDXa6tWr4yMjI6OPHDkStXz58qbJycl2hfXoo6Kionv37p0xZ86c5iUd365du+wvvviiAQD8+9//9kxOTnYEgLNnzzr6+vpefzpd8+bNr505c8YRMF8cTJ48OeXEiRNR7u7uBYsWLfIqbOfs7FwQHh4eO3HixEtF+/L19b125MiRmJ49e2ZOmDDBf/Pmzad++umnmDfeeKM5AKxfv77+yZMnnY8dOxYdHR0ddeTIEdfvvvvOreh5AGDv3r2umzdvbnj8+PGorVu3njp27Fi9m/n+OKInIqKylWPkbSu3Wo/+008/jX/66ad9X3/99WahoaFpDg4OCpjL0BYlIgoA3t7e1wYMGHAFAMaNG3dx6dKlTQCcA4DHHnvshgRf6MEHH0wDgKCgoKwrV66YGjZsWNCwYcMCJyenggsXLtht27at/t69e+sHBgYGAkBWVpYpJibm+oyDtd27d7sNGjQozc3NTQHovffem1bOr+wPmOiJiKjGqox69MHBwVf3798fBwDHjh1z2r59ewPAXKK2cAQPAImJiY4tWrTILe581u/d3d0LSurL2dlZAcBkMsHR0fH6lYTJZEJubq6oKqZPn5703HPPlVksqLKKznHqnoiIaqzKqEd/9uxZewDIz8/Hyy+/3OyJJ544DwCjRo1KW79+vWd2drbExMQ4xsfHO/fr1+8KACQlJTnu3LmzHgB88cUXnr169bphxH0zBg0adHnVqlWN09PTTQDwyy+/OBTGV1S/fv0yv//+e4+srCxJT0837dy5s8HN9MkRPRER1ViVUY/+008/9fzkk0+aAMDgwYMvPfPMMxcBoEePHleHDx+eGhAQ0MnOzg5LlixJsLc3p8XWrVtf/fTTTxtNnTrVr1WrVjmzZs1KqYzPM3LkyMuRkZHOt99+ewcAcHV1LVi9evUvPj4+Nywk7Nu3b1ZoaGh6YGBgJx8fn5wuXbpc8fDwyK9on6xHT0RkUKxHf3NiY2MdhwwZ0i4uLi7SFueviPT0dJOHh0dBRkaG6c4772z/4YcfJvTp0yeraLvS6tFzRE9ERFRDjR071i8uLs4lJydHHn744YvFJfmyMNETEVGtY8t69O3bt79WlaP55ORku379+rUvun337t2xmzdv/uVWz89ET0REhlAT69GXh7e3d35MTEyUrc7PVfdEREQGxkRPRERkYEz0REREBsZET0REZGBM9EREZGisR09ERGRgRqhHX1BQgPz8Cj8UDwB/XkdEROXw0v6XfE9eOlmp9ejbNmybtaD3AtajL6EefWxsrOOgQYPa9erVKyM8PNztm2++ORkQEHCtuLal4YieiIhqtLpaj74wjvHjx1+Mjo6OupkkD3BET0RE5VCekbet1NV69ADQrFmza3ffffeV8nxPJWGiJyKiGqsu16MHzNXtytOuNJy6JyKiGqsu16OvLBzRExFRjVWX69FXFtajJyIyKNajvzk1qR59eZVWj95mU/ci4isiu0QkWkQiRWSaZfsCETkmIkdEZLuIlLhSkoiIiG6NLafu8wDMVNVDIuIOIFxEdgBYpKovAYCIPANgHoDJNoyDiIgMpq7Uo/f29r65p+RYsVmiV9UkAEmW1xkiEg3AR1Wta+7WA1Dz7x0QEVGNx3r0xauSxXgi4g8gGMBPlvcLATwGIB3An0s4ZiKAiQDQsmXLqgiTiIjIcGz+8zoRcQOwDsB0Vb0MAKo6V1V9AawG8HRxx6nqClXtoao9vLy8imtCREREZbBpohcRB5iT/GpVXV9Mky8AjLJlDERERHWZLVfdC4BPAESr6hKr7e2smg0DEGOrGIiIiOo6W47oewMYB6C/5ad0R0RkMIA3RCRCRI4BGABgmg1jICIiqjI1sT69LVfd/wCguIcPf2urPomIiKrTO++8U2zJ2erER+BSqdI2n8K1xFsqnFQrOTavhwZD21R3GEQ1RuKLc31z4uIqtR69U7t2Wc1fW2jzevQhISHtg4KCso4ePeqamppq/9lnn/2ycOHCZrGxsS73339/6tKlSxMB4O9//3vT1atXNwaAcePGpcybN+/85cuXTcOGDWudlJTkWFBQILNnz0586qmnLs2aNavZtm3bGuTk5Jh69OiRuXr16gSTyYRRo0b5F9an37Nnj+v06dNbZmVlmRwdHXXv3r2xDRs2vKFIzUMPPeRX+Az/c+fOOUyYMOH84sWLk27t2/0di9oQEVGNdqv16AHA0dGxICwsLHb8+PEpDzzwQNuPPvro15iYmMg1a9Y0Tk5Ottu3b5/rF1980Sg8PDw6LCwseuXKlV779+93Wb9+fX1vb+/c2NjYqLi4uMiRI0deBoDnnnvufERERHRcXFxkdna26csvv/Sw7u/q1asyZsyYNu+8886vsbGxUXv27Il1c3MrthLdmjVrEmJiYqI2bdp0skGDBnmTJk2q0AN+ysIRPZWKo1oiAoDyjLxt5Vbr0QPAiBEj0gCga9eu2W3bts328/PLBQBfX9+c06dPO+7evdtt8ODBafXr1y8AgPvuu+/Srl273IcNG5Y+d+5c3ylTpvjcf//96aGhoZkA8N1337kvWbLE++rVq6a0tDT7wMDAbJifDQMAOHbsmHOTJk1y+/btmwUAnp6epZabzcrKklGjRrV5++23fw0ICLh2019WMTiiJyKiGsu6Hn1sbGxUx44dsytajx74Y514JyenP9SJz8vLK7HAW5cuXXIOHToUFRQUlD137lyfWbNmNcvKypKZM2f6rV+//tSJEyeixo4de+Hq1at/iElVISLlfvLruHHj/IYOHXpp+PDhJZbbvVlM9EREVGNVRj368ujfv3/mt99+2yAjI8N0+fJl07ffftvwz3/+c0Z8fLyDu7t7wdSpU1OnT59+7siRI65ZWVkmAPD29s5LT083bd68uWHR83Xt2vXquXPnHPfs2eMKAJcuXTLl5uYW2/frr7/ulZmZaffaa68l38pnKAmn7omIqMaqjHr05dGnT5+sRx999GL37t07AubFeL17985et25d/RdeeKGFyWSCvb29Llu2LKFx48b5Y8aMSQkMDOzUokWLa4UxWXN2dtbVq1efeuaZZ1pevXrV5OzsXLB3794THh4eN0zhv/fee94ODg7aoUOHQACYMGFCyuzZs1Nu5fNYYz16IiKDYj36uqNa6tETERFR9ePUPRER1Tq2rEdvK+vWras/d+7cFtbbfH19c3bs2HHKlv0y0RMRkSHU9Hr0o0aNujxq1Cib1Z0vCafuiYiIDIyJnoiIyMCY6ImIiAyMiZ6IiMjAmOiJiKjWCQkJab93794bqun17du37YULF+yKbp8xY0bzefPmNb3VfuPj4x1CQ0Nb3+p5qhJX3RMRkWHs2bPnpC3P7+/vn7tt27ZbegpfVWOiJyKiMv1nZbRv6tnMSq1H7+njlnX3Yx1LrYoXGxvrOGTIkHZxcXGRADBv3rymmZmZ10fs+fn5eOCBB/xbtGhxbenSpYk+Pj5BYWFh0c2aNct7/vnnvdesWdO4efPm1xo1apQbHBycBQCvvvpqk88++8zLzs5OAwICrm7ZsuX01q1b3WbOnNkSMBfIOXDgQExxteOLxlMbMNETEVGtlJubK8OHD28VGBiY/eabb/6hIMy+fftcN2zY4Hn8+PGo3NxcdOvWLbAw0S9dutQ7ISHhuIuLixZO8y9evNh76dKlCQMGDLiSnp5ucnV1LbWsbG3CRE9ERGUqa+RdHaZOneo3fPjw1KJJHgB27drlNnjw4DR3d/cCABgwYEBa4b727dtnjxgxotWwYcPSxowZkwYAd9xxR+asWbN8H3zwwdRHHnnkUps2bQyT6LkYj4iIaix7e3stKPg951rXfe/Ro0fmvn376mdlZRVbjL6kGvW7du2K+8tf/pISHh5er2vXroG5ubl47bXXkj/++OOE7OxsU69evToePnzYubI/S3VhoiciohqrRYsWeampqfbJycl22dnZ8v3333sU7ps0adKFAQMGpA8ZMqRN0Vrv/fv3z9y6dWuDzMxMuXTpkmnHjh0NAPM9/VOnTjkOHTo0Y9myZb9lZGTYpaen20VGRjqFhIRkL1y4MDkoKOhKRESEYRI9p+6JiKjGcnJy0pkzZyaFhIR0bNGiRU7btm2vWu//+9//fu7ZZ5+1GzlyZKuNGzf+Uri9T58+WSNGjEjt3LlzJx8fn5yQkJBMAMjLy5NHH320VUZGhp2qyqRJk841btw4f+bMmc0PHDhQ32QyaUBAQPbo0aPTS4pJRGp+fXcrrEdPRGRQRq5HX1327dvnOmPGDN+DBw/GVncs1liPnoiI6Bbt3bvXddy4ca2ffvrpc9UdS0Vw6p7IQJJfew050THVHUaZnDp2gPeLL1Z3GEQl+vnnn10ee+yxVtbbHB0dC+Lj4yOqK6abxURPRERUREhISHZMTEyV1463BSZ6IgPhKJmIiuI9eiIiIgNjoiciIjIwJnoiIiIDY6InIqIabcGCBU3atWvXqW3btp3mz5/fBADOnTtn16tXr3Z+fn6de/Xq1S4lJeWGGvRlCQ4O7lD50dY8TPRERFRjHTx40HnlypVehw4dio6Ojo7ctm1bg+PHjzu9/PLLzfr165eRkJAQ0a9fv4x58+Z5l/eceXl5AIDDhw/X/N+iVgKuuiciojJ9/8E7vhfOJFRqPfrGvn5ZA6dML7Uq3vHjx126d++eWViFrnfv3hlr1qxpsG3btgZ79uyJBYBJkyZd7Nu3b3sAZ2fMmNH89OnTTsnJyQ5JSUmOzzzzTPLMmTMvbNmyxX3BggXNmjRpkhsVFeV66tSpSFdX1+CsrKzDW7ZscX/llVeae3l55UZFRbkOHjz4UlBQUPayZcua5uTkyIYNG0516tQpJzEx0X78+PF+Z8+edQSAJUuW/DpgwIArxcWdmJhoP3r06FZpaWn23bp1y9q9e3f98PDw6GbNmuVV5ndYHhzRExFRjdWtW7fsn376yT05OdkuIyPDtGPHDo8zZ844Xrx40d7Pzy8XAPz8/HJTU1OvD1yjo6Nddu7cGffjjz/GLFq0qHl8fLwDABw7dqzeokWLzp46dSqyaD8xMTEuH3zwwZno6OjItWvXNjpx4oTz8ePHo8eNG3dh8eLFTQBg0qRJvjNmzDgXERERvWHDhlOTJ0/2LynuOXPmNO/bt29GVFRU9MiRIy8lJSU5VvqXU04c0RMRUZnKGnnbSvfu3a9OmzYtuX///gGurq4FgYGBWfb2paeuQYMGpbm5uambm1venXfeeXnfvn31GjZsmN+lS5crHTp0uFbcMUFBQVcKLxxatmyZM2jQoHQA6Nq1a/aePXvcAWD//v314+LiXAqPyczMtLt06ZKpYcOGN9Su//nnn902btx4EgBGjx59uX79+vk3/SXcIpuN6EXEV0R2iUi0iESKyDTL9kUiEiMix0Rkg4g0sFUMRERU+z377LMXoqKiosPCwmI9PT3z27Vrd7VRo0Z5CQkJDgCQkJDg4OnpeX1KvGgd+sL3rq6uNyTkQk5OTtcrvJlMJjg7O2vh6/z8fAEAVUVYWFh0TExMVExMTNT58+ePFZfkC9vWFLacus8DMFNVOwK4A8BfRCQQwA4AnVW1C4ATAF6wYQxERFTLnT171h4A4uLiHLdu3drgiSeeSB04cGDa8uXLGwHA8uXLG4WGhqYVtv/uu+8aZGVlSXJyst2PP/7o3qdPn2Lvo1dUnz59Lr/55ptNCt8fOHDApaS2ISEhmatWrfIEgPXr19e/fPlyhX8VUFlsNnWvqkkAkiyvM0QkGoCPqm63avYjgNG2ioGIiGq/YcOGtUlLS7O3t7fXd95551cvL6/8V155JWnEiBFt/Pz8Gjdv3vzaxo0bTxW2Dw4OvnL33Xe3S0xMdJw1a1aSv79/bkREhPOtxrFixYozTz75ZMuAgIDA/Px86dmzZ0avXr1+La7tG2+8kTh69OjWgYGBDe+8885MLy+v3AYNGlTL9H2V1KMXEX8Ae2EeyV+22r4ZwBpV/Xcxx0wEMBEAWrZseVtCQoLN4yQiMpK6WI9+xowZzd3c3PLnz59fraVks7Ozxd7eXh0cHLBz5856Tz/9tJ8ti+SUVo/e5ovxRMQNwDoA04sk+bkwT++vLu44VV0BYAUA9OjRo+bc7CAiIirDyZMnHR988ME2BQUFcHBw0OXLl8dXVyw2TfQi4gBzkl+tquuttj8OYAiAu7UmrVg5hHpUAAAUcElEQVQgIqJabcmSJYlV2d+7777b6IMPPmhqve3222/PXLVq1a/R0dE1osxtuRK9iDgBGAXA3/oYVZ1fyjEC4BMA0aq6xGp7KIDnAfRV1aybC5uIiKj6TZs27eK0adMuVnccpSnviP4bAOkAwgHklPOY3gDGATguIkcs214EsBSAE4Adlp88/Kiqk8sdMREREZVbeRN9C1UNrciJVfUHAFLMrm8rch4iIiK6eeX9Hf0BEQmyaSRERERU6Uod0YvIcQBqaTdeRE7DPHUvANTy0BsiIiKqocoa0Q8BMBTAIABtAQywvC/cTkREZFOsR39rSh3Rq2oCAIjIKlUdZ71PRFbBvNiOiIgMLnXtCd/c5CuVWqbWwbtelufogFKL5VjXo3d2di7o27dvwIgRI9Lff/99r379+mW89tprcS+++KL3vHnzvD/44IOz5ek3Ly8P9vb2VVaPPjc3Fw4ODlXRVbHKe4++k/UbEbEDcFvlh0NERPQ763r0Dg4Of6hHP2nSpIuAuR79d9991xAwPxlv+PDhre64444APz+/zosXL24MAFu2bHHv2bNnwNChQ1u1b9++EwC4uroGF+67/fbb2w8ePLi1v79/56lTp/p88MEHnkFBQR0DAgICIyMjnQBzjfmBAwe26dy5c8fOnTt33L59e72S4p4xY0bzRx55xK93797tRo4c2crW31NpyrpH/wLMP4lzEZHL+H0V/TVYnlpHRETGV9bI21a6deuWPX/+fJ/k5GS7evXq6Y4dOzy6du16pax69OHh4dEZGRl2wcHBgaNGjUoHzPXoDx8+HFlcqdqYmBiXtWvXnm7SpEmen59fkJOT04Xjx49HL1iwoMnixYubfPrpp2cK69EPHDgwMy4uznHgwIHtTp8+fUNt+0LHjh1z/emnn2Lc3Nyq9cFwZU3dvw7gdRF5XVVZZY6IiKpUba1HDwChoaFp1Z3kgfL/jv5FERkJoA/Mq/D3qepG24VFRERk9uyzz1549tlnLwDA008/7dOiRYtrhfXo/fz8cqu6Hn15k3e9evVK7K8qlfce/fsAJgM4DiACwGQRed9mUREREVnUxnr0NUl5R/R9YS4xqwAgIv+COekTERHZVG2sR1+TlKsevYisB/Cs1c/t/AC8oaqP2Dg+AOYytWFhYVXRFRGRYbAefd1RGfXoGwGIFpGfLe9vB/A/EdkEAKo67JajJCIiokpX3kQ/z6ZREBERVYKaVI++KuMoTbkSvarusUzXt1PVnSLiAsBeVTNsGx4REVHNVRvq0Zdr1b2IPAVgLYDllk0tAPDndURERDVceX9e9xcAvQFcBgBVjQPQpNQjiIiIqNqVN9HnqOr1pwmJiD3MD84hIiKiGqy8iX6PiBQ+8/5eAF8D2Gy7sIiIiKqWUcvWljfRzwGQAvNDciYB+BbA32wVFBERUVWrqrK1Va28q+4LRGQjgI2qmmLjmIiIqIbZuHGj7/nz5yu1Hn2TJk2yhg8fXmpVvNjYWMfQ0NB2ISEhmYcOHXLr2LFj1oQJEy7Mnz/f5+LFi/aff/756cDAwJwxY8b4//rrr04uLi4FK1asSOjZs2f21q1b3WbOnNkSMD/v/sCBAzEmkwmhoaFt09PT7fLy8mTevHmJY8eOTQPMZWuzsrIOA8Df/va3pl999VUjEcHdd9+dvmzZshtq3cfHxzuEhoa2K3wfFxfnEh0dfTwgIKDYwjnVpawytQLgZQBPw1yiVkQkH8A/VXV+FcRHRER13JkzZ5zXrFlz+rbbbkvo0qVLx9WrVzcKCwuL+eKLLxosXLiwmY+Pz7WuXbtm7dy589SmTZvcH3/88VYxMTFRixcv9l66dGnCgAEDrqSnp5sKi9ps3br1pKenZ0FSUpJ9z549Ozz66KNpJtPvE9xfffVV/a1btzYMDw+PcXd3Lzh37pxdcXH5+/vnxsTERAHA66+/7rVv3z73mpbkgbJH9NNhXm1/u6r+AgAi0hrAByLyrKq+besAiYio+pU18rYlHx+fnJCQkGwACAgIyO7fv/9lk8mE7t27Z7366qvNz54967Ru3bqTADBs2LCMiRMn2l+8eNHujjvuyJw1a5bvgw8+mPrII49catOmTUFOTo5Mnz69xY8//uhmMplw/vx5x99++82+ZcuW16vf7dixo/7YsWMvuLu7FwBA06ZN80uLb/v27fVWrlzp9eOPP9bIqf+y7tE/BuCRwiQPAKp6GsBYyz4iIiKbcnR0LLaErJ2dHfLz84ut2SIi+tprryV//PHHCdnZ2aZevXp1PHz4sPPy5cs9L168aH/8+PHomJiYqEaNGuVmZ2f/IReq6g2lbkuSkJDgMGnSJP81a9ac8vDwqBFlaYsqK9E7qOoNBQ0s9+kdbBMSERFR+d1xxx0Zn332WSMA2LJli3vDhg3zPD09CyIjI51CQkKyFy5cmBwUFHQlIiLCOT093a5x48a5Tk5OunnzZvfExETHoucLDQ29vGrVqsYZGRkmAChp6j4nJ0dGjhzZesGCBWe7dOmSY9tPefPKmrov7V5DjbsPQUREdc+bb76Z+Oijj/oHBAQEuri4FHz++ee/AMA//vGPJgcOHKhvMpk0ICAge/To0elpaWl2gwYNatu5c+eOnTp1ymrVqtXVoucbPXr05UOHDrl269ato4ODg95zzz3p77333g2L8Xbu3FkvIiKi3quvvtr81VdfbQ4A27Zti/P398+1/acuv1LL1FoW3l0pbhcAZ1WtklE9y9QSEVVcXSxTW1fddJlaVS12uoKIiIhqh/KWqSUiIqqzxo0b1/LgwYNu1tumTJlyrqZXrgOY6ImIiMpUk+rLV1R5H4FLRER1T0FBQUH5fmdG1cbyNyrxp31M9EREVJKIlJQUDyb7mqugoEBSUlI8AESU1IZT90REVKy8vLwnk5OTP05OTu4MDgxrqgIAEXl5eU+W1ICJnoiIinXbbbedBzCsuuOgW8MrNCIiIgOzWaIXEV8R2SUi0SISKSLTLNsfsLwvEJGbfpADERERlc2WU/d5AGaq6iERcQcQLiI7YF4wMBLAchv2TURERLBholfVJABJltcZIhINwEdVdwAod2UgIiIiunlVco9eRPwBBAP4qQLHTBSRMBEJS0lJsVVoREREhmbzRC8ibgDWAZiuqpfLe5yqrlDVHqraw8vLy3YBEhERGZhNE72IOMCc5Fer6npb9kVEREQ3suWqewHwCYBoVV1iq36IiIioZLZcdd8bwDgAx0XkiGXbiwCcAPwTgBeArSJyRFUH2jAOIiKiOsuWq+5/AFDS0voNtuqXiIiIfscn4xERERkYEz0REZGBMdETEREZGBM9ERGRgTHRExERGRgTPRERkYEx0RMRERkYEz0REZGBMdETEREZGBM9ERGRgTHRExERGRgTPRERkYEx0RMRERkYEz0REZGBMdETEREZGBM9ERGRgTHRExERGRgTPRERkYEx0RMRERkYEz0REZGBMdETEREZGBM9ERGRgTHRExERGRgTPRERkYEx0RMRERkYEz0REZGBMdETEREZGBM9ERGRgTHRExERGRgTPRERkYHZV3cAREQ1zSubIxGVeLm6wwAABDavj5eHdqruMKgW44ieiIjIwDiiJyIqgiNoMhKO6ImIiAyMiZ6IiMjAbJboRcRXRHaJSLSIRIrINMt2TxHZISJxln83tFUMREREdZ0tR/R5AGaqakcAdwD4i4gEApgD4D+q2g7AfyzviYiIyAZsluhVNUlVD1leZwCIBuAD4H4A/7I0+xeA4baKgYiIqK6rknv0IuIPIBjATwCaqmoSYL4YANCkhGMmikiYiISlpKRURZhERESGY/NELyJuANYBmK6q5X4ChaquUNUeqtrDy8vLdgESEREZmE0TvYg4wJzkV6vqesvmcyLSzLK/GYDztoyBiIioLrPlqnsB8AmAaFVdYrVrE4DHLa8fB/CNrWIgIiKq62z5ZLzeAMYBOC4iRyzbXgTwBoCvROQJAL8CeMCGMRAREdVpNkv0qvoDAClh99226peIiIh+xyfjERERGRgTPRERkYEx0RMRERkYEz0REZGBMdETEREZGBM9ERGRgTHRExERGRgTPRERkYEx0RMRERkYEz0REZGBMdETEREZGBM9ERGRgTHRExERGRgTPRERkYEx0RMRERkYEz0REZGBMdETEREZGBM9ERGRgTHRExERGRgTPRERkYHZV3cAREQ1zb6vTuDCmczqDgMA0NjXDXc9GFDdYVAtxhE9ERGRgXFET0RUBEfQZCQc0RMRERkYEz0REZGBMdETEREZGBM9ERGRgTHRExERGRgTPRERkYEx0RMRERkYEz0REZGBiapWdwxlEpEUAAlWmzwApFdTOFWptn/O2hB/TYqxOmOp6r4bA7hQhf3VVX6q6lXdQVD1qhWJvigRWaGqE6s7Dlur7Z+zNsRfk2Kszliqum8RCVPVHlXVH1FdVlun7jdXdwBVpLZ/ztoQf02KsTpjqUnfAxFVolo5oiei2o0jeqKqU1tH9ERUu62o7gCI6gqO6ImIiAyMI3oiIiIDY6InIiIysDqT6EWko4h8KCJrRWRKdcdDJePfioio8lRLohcROxE5LCJbbuEcn4rIeRGJKGZfqIjEishJEZkDAKoaraqTATwIgKt9y0FEGliSbYyIRIvInTd5Hv6tqFQiUk9EwkVkSHXHQmQ01TWinwYgurgdItJERNyLbGtbTNPPAYQWc7wdgPcBDAIQCOAREQm07BsG4AcA/7mV4OuQdwFsU9UOALqiyN+MfysqSUkXd8Vd2Fk8D+Crqo2SqG6o8kQvIi0A3Afg4xKa9AXwjYg4W9o/BWBp0UaquhdAajHHhwA4qaqnVfUagC8B3G85ZpOq9gIw5pY/iMGJSH0AfwLwCQCo6jVVTSvSjH8rKsnnKHJxV9KFnYjcAyAKwLmqDpKoLrCvhj7fATAbgHtxO1X1axFpBeBLEfkawAQA91bg/D4Azli9/w1ATxHpB2AkACcA395E3HVNawApAD4Tka4AwgFMU9UrhQ34t6KSqOpeEfEvsvn6hR0AiEjhhZ0bgHowJ/9sEflWVQuqMFwiQ6vSRG+5/3ZeVcMt/zMvlqr+w/I/gQ8AtFHVzIp0U/wpdTeA3RU4T11nD6A7gL+q6k8i8i6AOQBesm7EvxVVQLEXdqr6NACIyP8DcIFJnqhyVfXUfW8Aw0QkHuZp2v4i8u+ijUTkLgCdAWwA8HIF+/gNgK/V+xYAEm8q2rrtNwC/qepPlvdrYU78f8C/FVVAsRd211+ofq6qN71Al4iKV6WJXlVfUNUWquoP4GEA/1XVsdZtRCQYwEcwT+mNB+ApIq9WoJuDANqJSCsRcbT0s6lSPkAdoqrJAM6ISHvLprthvo96Hf9WVEG8sCOqBjXxd/SuAB5Q1VOWKbzH8cda9AAAEfk/AP8D0F5EfhORJwBAVfMAPA3ge5hXiX+lqpFVFr2x/BXAahE5BqAbgNeK7OffiiqCF3ZE1YDPuieiSme5uOsHoDHMq+lfVtVPRGQwzAty7QB8qqoLqy9KorqBiZ6IiMjAauLUPREREVUSJnoiIiIDY6InIiIyMCZ6IiIiA2OiJyIiMjAmeiIiIgNjoiciIjIwJnoiIiIDY6InqsVE5J8ickhEbq/uWIioZmKiJ6qlRKQegCYAJgEYUs3hEFENxURPtYKIvC0i063efy8iH1u9XywiMyq5z8xKPl8DEZlq9d5fRCLKeayLiOwREbvCbap6BUAzALsBLBURRxHZKyL2lRk3EdVuTPRUWxwA0AsARMQEc7GUTlb7ewHYXw1xVUQDAFPLbFW8CQDWq2p+4QYRaQRzBcEMAPmqeg3AfwA8dKuBEpFxMNFTbbEflkQPc4KPAJAhIg1FxAlARwCHRWSjiISLSKSITCw8WETeLDKa/ruIzBSRsSLys4gcEZHl1iNmq7bFtrGMyKNF5CNLf9tFxMWy7yURiRGRHSLyfyIyC8AbANpYzrPIcnq74o4vxhgA3xTZ9jcAbwGIBBBo2bbR0paICAATPdUSqpoIIE9EWsKc8P8H4CcAdwLoAeCYZUQ7QVVvs2x7xjLqBYAv8ceR7oMAwizbeqtqNwD5KJIkRaRjGW3aAXhfVTsBSAMwSkR6ABgFIBjASEssADAHwClV7aaqz5V0fNHPbqnd3lpV4622+Vu+hzUAovH77EYEAC7MI6LreC+PapPCUX0vAEsA+Fhep8M8tQ+Yk/sIy2tfmBPpRVU9LCJNRKQ5AC8AlwAEAbgNwEERAQAXAOeL9Hl3GW1+UdUjltfhAPxhvq3wjapmA4CIbC7lMxV3fFGNYb4IsPYqgPmqqiJyPdGrar6IXBMRd1XNKKVfIqojmOipNim8Tx8E88j1DICZAC4D+FRE+gG4B8CdqpolIrsBOFsdvxbAaADeMI/wBcC/VPWFUvosq02O1et8mC8EpAKfqbjji8qG1ecQkW4wzxT0EZH3LfuOW7V3AnC1AjEQkYFx6p5qk/0w/4wsVVXzVTUV5gVud8I8le8B4JIlyXcAcEeR478E8DDMyX4tzAvXRotIEwAQEU8R8StyTHnaFPUDgKEi4iwibgDus2zPAOBe0Q+tqpdgvpdfmOzfBDBUVf1V1R9AV1hG9JZbFSmqmlvRfojImJjoqTY5DvM09o9FtqWr6gUA2wDYi8gxAAuKtIOqRsKcaM+qapKqRsG8oG275ZgdMP9czfqYMtsUpaoHAWwCcBTAepjXAqSr6kUA+0UkwmoxXnlth3kE3x9APVX9j1V/5wDUExFPAH8G8G0Fz01EBiaqWt0xEBmOiLipaqaIuALYC2Ciqh66hfMFA5ihquPKaLcewAuqGnuzfRGRsfAePZFtrBCRQJjvn//rVpI8AFgWE+4SETvr39Jbs6zO38gkT0TWOKInIiIyMN6jJyIiMjAmeiIiIgNjoiciIjIwJnoiIiIDY6InIiIyMCZ6IiIiA2OiJyIiMjAmeiIiIgP7/8yl8XJBE134AAAAAElFTkSuQmCC\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": {},
   "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": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Depths (5 $\\\\sigma$) vs coverage on XMM-13hr')"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfoAAAEYCAYAAAC0mTTAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd8VGX2P/DPSSMJGXpoSUgIpAdCMyCwgqi0FVbKWlDcVVGQ5aeAuBZ22V0FO+iyrl30C4KiUgQRFJQqKxh6KiEQpCQYWhhIKCHn98e9g0NISJ1kMvm8X695OfPMvc89907kzH2eO/eIqoKIiIhck1tNB0BERESOw0RPRETkwpjoiYiIXBgTPRERkQtjoiciInJhTPREREQujImeiIjIhTHRExERuTAmehcnIpkicmsV9fWiiEysir6K6XuriMQ4om+iqvz/gKi2YaKvZuY/OPkiYhWR0yKyWUTGiUiVfBaO+gdNRPwB3A/gXbu2dSJyXkTOmo+0SmziNQDPVTZOqnoi4mf+XY2ya7OIyC8iMtJ8nSkiF0WkWZF1d4qIikhIeZYrIY4JIpIgIhdE5ONi3v9ERLJE5IyI7BWRMZXYbSKXwURfM4aoqgVAMICXADwF4MOaDalUfwbwjarmF2mfoKp+5iOiEv0vA3CziLSqRB9OTUQ8ajqGilDVswAeAfBv8wsfALwCIEFVv7Rb9ACAe2wvRKQDAJ9iuizrckUdBTAdwJwS3n8RQIiqNgAwFMB0Eelahn6vUls/J6KSMNHXIFXNVdVlAO4C8CcRiQUAEWktIotEJEdEDojIY7Z1zDOiZ0QkWUROichHIuJtvjcPQBsAy80z7L+aq3USkd0ikisiC23Lm+s8JSJHzBGGNBG5pYRwBwFYX9F9FRFPEZlhxn/JPHtTEdllHovzALYB6F/Muk+LyJdF2v4tIrPLuQ8QkSARWWwe2xMi8qbZHmWOUJwWkSQRGVqWbZuvS/u8nhKR3QDOiYiH2WeGGW+yiAyzW76LiOww3/vC/Lyml2VbxexrsftkF9eUkv4uilLV7wCsADBbRPoCuBPAX4osNg/GqI/NnwDMLaa7si5XNIbFqroUwIkS3k9S1Qu2l+ajnd0i1/v/4JrPqbR4iGoNVeWjGh8AMgHcWkz7LwAehfHlaxuAaQC8AIQC2A9ggN36iQCCADQB8COA6SX1b77eCqC1uXwKgHHmexEADgFobb4OAdCuhLhzANxQpG2d2X7cjKPvdfb7ZQA/mXHXB7AGwGIAoXbLzAYwq5h1gwHkAWhgvnYHkAWgRzn3wR3ALgCvmzF4A+gNwBPAPgDPmse8HwCr2XeJ2zZfl+Xz2mnut4/Z9kfz83CD8SXvHIBW5voHATxuxjQcwEXb51vatorsa4n7VNrfxXU+w8bmvh8H8EBxf9cA0gBEmcfpkHn8FMaZdpmXKyWO6QA+LuG9t8zPSwFsB+BXlv0t7nPigw9XefCM3nkchfEP0A0A/FX1OVW9qKr7AbwP4G67Zd9U1UOqehLADNgNg5ZgtqoeNZdfDqCT2X4ZQD0A0SLiqaqZqppRQh+NYCQKe0/BSDYBAN6DMZLQruiKImIB8BiA0Wbc5wAsAtDE3D8bq7mdq6jqQRj/aN9hNvUDkKeqP5VzH+Jh/EP/pKqeU9XzqroJxhcGPwAvmcf8BwBfA7inlG0DZfu8Zpv7nW/uzxfm51GoqgsBpJux9QDgYS5/SVUXw0hOKMe2bErcpyJxFfd3USxVPQUgCYAvjC9pxbGdrd8GIBXAkUouVy6qOh6ABcDvzBgv2L1d2v5e9TkRuQomeucRAOAkjDOb1uZw62kROQ3jrKyF3bKH7J4fhJG8rifb7nkejAQAVd0HYCKAfwL4VUQ+E5GS+joF4x/QK1R1i6paVfWCqv4fjLP6wcWsexOA/aqabtfWuEhcMPs/XcL2F+C3JDXKfF3efQgCcFBVC4q0twZwSFUL7doOwvhMSty2qbyfF0TkfjEuPrMtHwugmRnHEVXVEtYty7bKuk9ACX8XJRGR+2CMmKyBMUJTnHkwjtGfcf3h+OsuJyL3ym8Xea68XlxFqepl8wtcIIxRMpvS9vcQiFwQE70TEJEbYPwDvAnGPzYHVLWR3cOiqvYJNMjueRsYowE29kmiVKq6QFV747eh05L+Ad8NILy07gBIMe3+ML4oAABERAAMg3GGaS8KxtB6cb4A0FdEAs11ryTbcuzDIQBtipl/PQogSK7+5UMb/HaWWeK2UbbP68pnIiLBMM7CJwBoqqqNYEzFCIxh8QDz+NjYf9Zl2VZZ96lcRKQ5jCmPhwGMBXCniNxUdDlzBOQAjC98JZ31l7qcqs7X3y7yHFSRmGGMjlwzwnQd5fp/h6i2YKKvQSLSQERuB/AZgE9UdQ+Modoz5oVBPiLiLiKx5pcBm7+ISKCINIFxRrfQ7r1jMIbTy7L9CBHpJyL1AJwHkA9jKLw43wDoY7duIxEZICLe5gVm98I4c/+2mHUTAXQRkU4i4gPj6mi1j9uMoSuA1cVtXFVzYFwT8BGMZJdSgX3YCiOZviQi9c3YewHYAmOe/K9iXDTYF8AQGJ9Lidu267O0z8tefXPfc8z4H4BxRg8A/zNjn2Ae0z/AGNKvyLauu08V8CaApaq6VlWzAPwVwPvmcS/qIQD9zCma6ynrcgCMq+HNC+jcAbjb/vbM95qLyN1i/BTQXUQGwBiF+aGM+0fkspjoa8ZyEbHCOEObCmAWgAcAY9gRxj/InWCc8RwH8AGAhnbrLwDwHYwLsfbDuDjJ5kUAfzOHdqeUEkc9GD/vOw5jWLM5jC8OxZkLYLCZqAHjYq/p+O1ivP8H4A5Vvea39KqaAONagm/MeFsCGKyql+wWGwpgnaoeLbq+nQUwLuSyP6Mu8z7YHdv2MC5+PAzgLlW9aG5/kNnPWwDuV9XUUrZd1s/LfvlkADNhJPVjADrAmPKAGcdwGAnwNID7YIx6XCjvtsq4T2UiInfAuGjxSbv+P4Bx/KYVs+0M8zO/rrIuZ+dvML7IPQ3j2OSbbYDx5elRM6ZTMO7LMFFVvypH/0QuSa6eDiRnJyKZAMao6poa2PYLAH5V1Tcc0PcWAA+pamJV912bmcflHVX9qKZjIaLaib8VpTJT1ZLO9qui7+6O6rs2EZE+MH56dhzAvQA6AlhVo0ERUa3GRE/kXCIAfA7jivAMACPNOXEiogrh0D0REZEL48V4RERELqxWDN03a9ZMQ0JCajoMIqJaZdu2bcdV1b/0JUtcv7mHh8cHMH4CyhND51QIILGgoGBM165dfy1ugVqR6ENCQpCQUJ5f4RARkYgcrMz6Hh4eH7Rs2TLK39//lJubG+d5nVBhYaHk5OREZ2dnfwDjJ7XX4Dc0IiIqSay/v/8ZJnnn5ebmpv7+/rn47cZb1y5TjfEQEVHt4sYk7/zMz6jEfO6wRC9G3e+1IpJi1sJ+vMj7U8SoR97MUTEQERHVdY6coy8A8ISqbjfLlG4TkdWqmiwiQTDKU/7iwO0TERHVeQ47o1fVLFXdbj63AkjBbyUyX4dRFINDQkRELuC1b9NarEk5dlUp6zUpxyyvfZtWXBllpzB79uym999/f5ui7a+88or/m2++2bRoe1pamldYWFhMVWz7rrvuCt62bZt3VfRVmmqZoxeREACdAWwRkaEwam6XVI7Uts4jIpIgIgk5OTnl2+CmN4ADG65uO7DBaCcioirXqU2jvMmf7wy1Jfs1Kccskz/fGdqpTaO8mo6tvP7617/mTJgw4YQjt7Fw4cKDXbt2Pe/Ibdg4PNGLiB+ARQAmwhjOn4piKl4VparvqWo3Ve3m71/On4EGdAG++PNvyf7ABuN1QJfy9UNERGVya1QL66w7O+2f/PnO0H8tT2o9+fOdobPu7LT/1qgW1kr3feut7WJiYqLat28f89prrzUDAF9f384PP/xwYHR0dNSNN94YfvTo0RKnouPj4yM2bNjgCwBZWVkeAQEBHYou89lnnzXs1KlTZFZWlsfkyZNbT5s2rQUAbNy40TciIiK6U6dOkbNmzWpuWz4hIcG7Q4cOUZGRkdHh4eHRe/bsqXfmzBm3vn37to+IiIgOCwuLef/99xuXJSZHc2iiFxFPGEl+vqouBtAOQFsAu8wqbIEAtotIy6rc7pxz+7D1tmeN5P7DDOCLP2Prbc9izrl9VbkZIiKyc2tUC+uILoE5H/2Y2WpEl8CcqkjyADB//vzMpKSklJ07dya/++67LbKzs93z8/PdunTpkpecnJzSq1cv69NPP926ov3PnTu30auvvtpy9erV6a1atSqwf++hhx4KmTVr1i87d+68qsTzf/7zH//x48cfS01NTd69e3dK27ZtLy5evLhBy5YtL6WlpSWnp6cnDR8+/ExFY6pKjrzqXgB8CCBFVWcBgKruUdXmqhqiqiEwakd3UdXsqtx2bNNYTEmbi62xtwMbXsHW2NsxJW0uYpuW+DNDIiKqpDUpxyyLth/2f6BXSNai7Yf9i87ZV9TLL7/cIiIiIrpr165R2dnZnklJSd5ubm4YM2bMSQB48MEHT2zdutWvIn1v3rzZMnPmzJarV69O9/f3v2z/3okTJ9ytVqv773//+7O27djeu/HGG8/NnDmz1dSpU1ump6d7+fn5aZcuXfI3btzY4NFHHw1YtWqVX9OmTS8X3V5NcOQZfS8AowH0E5Gd5mOwA7d3RXyreLwWcT+mHPsBb8YNwpRjP+C1iPsR3yq+OjZPRFTn2ObkZ93Zaf8/hsQctQ3jVzbZf/3115b169dbEhISUtPS0pKjoqLy8/Pzr8ldxrll8Tw8PPTyZSPn5uXlXbVgmzZtLpw7d849MTHxmgvjVLXEfseNG3fyq6++2ufj41M4aNCg8GXLllk6dux4Yfv27ckdOnTInzp1asCUKVNalXN3HcKRV91vUlVR1Y6q2sl8fFNkmRBVPV7lGz+wAfGrX8CdIYPx7pkk3BkyGPGrX7j2Aj0iIqoSO3857Ws/J2+bs9/5y+lKzUOfPn3avWHDhpctFkvhjh07vHft2lUfAAoLC/HRRx81BoCPP/64aXx8fInTBEFBQRe2bt1aHwDmz59/1bx5YGDgxUWLFu174IEH2iYkJFyV7Js1a3bZz8/v8rfffutnbqeJ7b3k5GSvqKioC3/7299+7d+//+mdO3f6ZGZmeloslsLx48efnDhx4rGdO3dWyxx8aWrFve7L7ch2bL3tWXyeNhdjO47F52mfI/62ZxF/ZDvQ9qaajo6IyOVMGRBxrGjbrVEtrJWdpx8xYkTue++95x8eHh7drl2783FxcecAwMfHpzApKcknJiampcViubx48eL9JfXx9NNPH7vrrrtCP/vss6a/+93vrpk3j4uLuzB37tz9d911V7tly5ZddTHXhx9+mDlmzJgQHx+fwn79+l1Zd968eU2++OKLph4eHurv73/pxRdfPLpp06b6zzzzTKCbmxs8PDz0rbfeum6tgeuNQlSlWlGPvlu3blqeojZbs7ZiyvopeK3Pa4hvFX/NayKiukBEtqlqt4quv2vXrsy4uLiqH3WtAr6+vp3z8vJ21HQcFRUeHh69bNmyfZGRkReror9du3Y1i4uLCynuPZe8133iicSrknp8q3i81uc1JJ5IrOHIiIioruvZs2dYREREflUl+dK45ND9g7EPXtMW3yqeZ/NERC6iuLP50aNHt/n555+vuvr+0UcfPfb444879OY313Pbbbe1O3ToUD37thkzZhweMWJEtf30ziUTPRER1T3z5s1zuvopq1evzqjpGFxy6J6IiIgMTPREREQujImeiIjIhTHRExERuTAmeiIiqrzvn2+BtJVX3+42baUF3z9f4/Xo//e///l06tQpMjw8PLpfv37tT548eSX3PfPMMy3btGkTGxISErto0aIG5e174sSJrZcuXVol9/R3FCZ6IiKqvMBueVgyLvRKsk9bacGScaEI7Fbj9egffvjhkBkzZhzeu3dv8tChQ0/961//agkA27Zt8168eHGTtLS0pFWrVu2dOHFim4KCgtK6u6KgoABvvPHG0TvuuKNKqvQ5ChM9ERFVXsQgK4a9sx9LxoVi5dOtsWRcKIa9sx8Rg2q8Hn1mZqb3oEGDzgLA7bfffubrr79uDABffvllo+HDh5/08fHRyMjIi8HBwRfWrVtXPy0tzatt27Yxw4cPDwkPD48eOHBgqNVqdQOAgICADlOmTGnVtWvXiDlz5jQeMWJEiO2e+wEBAR0mTJgQ0KlTp8jY2NioTZs2+fbu3TssKCgo9pVXXvG3xfP3v/+9RWxsbFR4eHj0pEmTrlte98knn2zVtm3bmJ49e4YNGTKk7bRp08o9QsJET0REVSNikBVx9+Rgy9utEHdPTlUkeaDy9ejDwsLyFyxY0AgAPvnkkybZ2dleAHDkyBGvoKCgK3ena9269cVDhw55AcaXg3HjxuXs3bs32WKxFL766qtXErW3t3fhtm3b0h555JFTRbcVFBR0cefOnandu3c/++CDD4YsX748Y8uWLakvvfRSawBYvHhxg3379nnv3r07JSUlJXnnzp2+K1euLLbE7oYNG3yXL1/eeM+ePckrVqzI2L17d/2KHD8meiIiqhppKy3Y9ak/uj+ahV2f+l8zZ19Bla1HP2fOnMy3337bPyYmJspqtbp5enoqYJShLUpEFABatmx5sX///ucAYPTo0Sc2b958pf/777//mgRvc+edd54GgA4dOuR16dLlXOPGjQtbt25dUK9evcLjx4+7r1q1qsGGDRsaREdHR8fExERnZGR4p6amXlMiFwDWrVvnN2jQoNN+fn7auHHjwttuu+10mQ5YEbwzHhERVZ5tTt42XB/ax1oVw/f29egtFkthfHx8RHnr0Xfu3Pn8jz/+mA4Au3fvrvfdd981AowStbYzeAA4evSoV2Bg4KXi+rN/bbFYCkvalre3twKAm5sbvLy8rnyTcHNzw6VLl0RVMXHixKwnn3yy1GJBVVV0jmf0RERUeYcTfK9K6rY5+8MJNV6P/siRIx4AcPnyZfzjH/9o9dBDD/0KACNGjDi9ePHiJvn5+ZKamuqVmZnp3bdv33MAkJWV5bVmzZr6ALBgwYImPXv2PFuZ/bAZNGjQmXnz5jXLzc11A4ADBw542uIrqm/fvme//fbbhnl5eZKbm+u2Zs2aRhXZJs/oiYio8m75+zX16BExyFrZefqqqEc/Z86cJh9++GFzABg8ePCpxx577AQAdOvW7fwdd9xxMjw8PMbd3R2zZs066OFhpMXQ0NDzc+bMaTp+/Pjgtm3bXpgyZUpOZfbDZvjw4WeSkpK8b7jhhkgA8PX1LZw/f/6BgICAay7379OnT97AgQNzo6OjYwICAi507NjxXMOGDS+Xd5suWY+eiIhYj76i0tLSvG6//faw9PT0JEf0Xx65ubluDRs2LLRarW433nhjxDvvvHOwd+/e1/xk8Xr16HlGT0RE5KTuu+++4PT0dJ8LFy7I3XfffaK4JF8aJnoiIqp1HFmPPiIi4mJ1ns1nZ2e79+3bN6Jo+7p169KWL19+oLL9M9ETEZFLcMZ69GXRsmXLy6mpqcmO6p9X3RMREbkwJnoiIiIXxkRPRETkwpjoiYiIXBgTPRERVdrs7bNbrDu07qp72687tM4ye/ts1qOvYUz0RERUaR39O+ZN3TQ11Jbs1x1aZ5m6aWpoR/+OrEdfBQoLC3H5crlvigeAiZ6IiKpA36C+1hm9Z+yfumlq6EtbX2o9ddPU0Bm9Z+zvG9SX9egrWI8+LS3NKzQ0NOa+++5rY1a68ypp2ethoicioirRN6ivdUi7ITnzU+a3GtJuSE5VJHmg7tajt8XxwAMPnEhJSUkODw+/WNJy18NET0REVWLdoXWW5RnL/e+NujdrecZy/6Jz9hVVV+vRA0CrVq0u3nLLLedKPUjXwTvjERFRpdnm5G3D9T1a9bBWxfB9Xa5HDxjV7cqy3PXwjJ6IiCptd85uX/ukbpuz352zm/Xo7ZSnHn1V4Rk9ERFV2mNdHrumHn3foL7Wys7T1+V69FWF9eiJiFwU69FXjDPVoy+r69Wjd9jQvYgEichaEUkRkSQRedxsf15EdovIThH5TkRKvFKSiIiIKseRQ/cFAJ5Q1e0iYgGwTURWA3hVVf8OACLyGIBpAMY5MA4iInIxdaUefcuWLSt2lxw7Dkv0qpoFIMt8bhWRFAABqmpfc7c+AOefOyAiIqfHevTFq5aL8UQkBEBnAFvM1zMA3A8gF8DNJazzCIBHAKBNmzbVESYREZHLcfjP60TED8AiABNV9QwAqOpUVQ0CMB/AhOLWU9X3VLWbqnbz9/cvbhEiIiIqhUMTvYh4wkjy81V1cTGLLAAwwpExEBER1WWOvOpeAHwIIEVVZ9m1h9ktNhRAqqNiICIiqusceUbfC8BoAP3Mn9LtFJHBAF4SkUQR2Q2gP4DHq3rD1vWHcD7j9FVt5zNOw7r+UFVvioiIAPz6xhstrGvXXnVve+vatZZf33ijxuvRVydnrE/vsESvqptUVVS1o6p2Mh/fqOoIVY0124eo6pGq3rZnoAUnF6RcSfbnM07j5IIUeAY61bEnInIZPnFxeUefejrUluyta9dajj71dKhPXFyN16OvTtVVn748XPJe997tGqHJqCicXJCC3O8ycXJBCpqMioJ3u0Y1HRoRkUuy3HyztfXLL+0/+tTTodkvvND66FNPh7Z++aX9lptvrvF69PHx8REPPfRQULdu3SJCQ0Nj1q9f79u/f/92wcHBsY899tiVm7b985//bBEWFhYTFhYW89xzzzUHgDNnzrj17du3fURERHRYWFjM+++/3xgApkyZ0io2NjYqLCws5p577gkuLDRqz9jXp1+/fr1v586dIyMiIqI7dOgQderUqWJz7l133RUcGRkZHRkZGd24ceO4J554olVlj5k9l0z0gJHs63dvBesPh1C/eysmeSIiB7PcfLO14R1/yDk1d16rhnf8IacqkjxQ+Xr0AODl5VWYkJCQ9sADD+T88Y9/bP/+++//kpqamrRw4cJm2dnZ7hs3bvRdsGBB023btqUkJCSkzJ071//HH3/0Wbx4cYOWLVteSktLS05PT08aPnz4GQB48sknf01MTExJT09Pys/Pd/vss88a2m/v/Pnzcu+997Z74403fklLS0tev359mp+fX7GV6BYuXHgwNTU1edmyZfsaNWpUMHbs2HLd4Kc0Lpvoz2ecxrktWbD0C8K5LVnXzNkTEVHVsq5da8ld+pV/4/tHZ+Uu/cq/6Jx9RVW2Hj0ADBs27DQAxMXF5bdv3z4/ODj4ko+PjwYFBV3Yv3+/17p16/wGDx58ukGDBoUNGzYs/P3vf39q7dq1li5duuRv3LixwaOPPhqwatUqv6ZNm14GgJUrV1o6duwYGR4eHr1582ZLYmKij/32du/e7d28efNLffr0yQOAJk2aFHp6epYYX15enowYMaLd66+//kt4ePjFSh6yq7hkorfNyTcZFYWG/UOuDOMz2RMROYZtTr71yy/tb/nss0dtw/iVTfb29ejT0tKSo6Ki8stbjx64uk58vXr1rqoTX1BQUGKBt44dO17Yvn17cocOHfKnTp0aMGXKlFZ5eXnyxBNPBC9evDhj7969yffdd9/x8+fPXxWTqkJEynzn19GjRwcPGTLklCPm910y0V86bL1qTt42Z3/psFNdH0FE5DLyd+3ytZ+Tt83Z5+/aVeP16MuiX79+Z7/55ptGVqvV7cyZM27ffPNN45tvvtmamZnpabFYCsePH39y4sSJx3bu3Ombl5fnBgAtW7YsyM3NdVu+fHnjov3FxcWdP3bsmNf69et9AeDUqVNuly5dKnbbL774ov/Zs2fdX3jhhezK7ENJXLIevaVP0DVt3u0acZ6eiMhBmk+ceE09esvNN1srO09fFfXoy6J37955o0aNOtGlS5coABg9enROr1698hctWtTgmWeeCXRzc4OHh4e+9dZbB5s1a3b53nvvzYmOjo4JDAy8aIvJnre3t86fPz/jsccea3P+/Hk3b2/vwg0bNuxt2LDhNfP0b775ZktPT0+NjIyMBoAHH3ww569//WtOZfbHHuvRExG5KNajrztqpB49ERER1TyXHLonIiLX5sh69I6yaNGiBlOnTg20bwsKCrqwevXqDEdul4meiIhcgrPXox8xYsSZESNGOKzufEk4dE9EROTCmOiJiIhcGBM9ERGRC2OiJyKiSvvpq4wWB3Yfv+oueAd2H7f89FVGnSpT64yY6ImIqNJatG2Y9/3HyaG2ZH9g93HL9x8nh7Zo29AhZWrj4+MjNmzYcM1d9/r06dP++PHj7kXbJ0+e3HratGmV/tKRmZnpOXDgwNDK9lOdeNU9ERFVWtuOzay3/Dl6//cfJ4dG9GiZk/ZTtv8tf47e37Zjs2q99/j69ev3ObL/kJCQS6tWrarUXfiqG8/oiYioSrTt2Mwa0aNlzu4fDreK6NEypyqSfFpamldYWFiM7fW0adNaTJ48+UpJ2suXL2P48OEhtrryAQEBHbKysjwA4KmnnmoZEhIS27Nnz/D09PR6tnWmT5/evF27djHh4eHRt99+eygArFixws9WEz4qKiq6pNrxReOpDXhGT0REVeLA7uOWtJ+y/Tv2C8xK+ynbPzCyidWRZ/SXLl2SO+64o210dHT+yy+/fFVBmI0bN/ouWbKkyZ49e5IvXbqETp06RXfu3DkPAGbPnt3y4MGDe3x8fNQ2zD9z5syWs2fPPti/f/9zubm5br6+vsXWjq+NeEZPRESVZpuTv+XP0ft/d2f4UdswftEL9KrS+PHjg4tL8gCwdu1av8GDB5+2WCyFTZo0Kezfv/+VOuURERH5w4YNa/vWW2818fT0VADo0aPH2SlTpgRNnz69+fHjx92vVzu+tmGiJyKiSjt2INfXfk7eNmd/7EBupcrUenh4aGHhbyfX9nXfu3Xrdnbjxo0N8vLyii1GX1KN+rVr16b/5S9/ydm2bVv9uLi46EuXLuGFF17I/uCDDw7m5+e79ezZM2rHjh3elYnbmTDRExFRpfX4Q7tjRYfp23ZsZu3xh3bXlK8tj8DAwIJIXNRtAAAe3ElEQVSTJ096ZGdnu+fn58u3337b0Pbe2LFjj/fv3z/39ttvb1e01nu/fv3OrlixotHZs2fl1KlTbqtXr24EGHP6GRkZXkOGDLG+9dZbh61Wq3tubq57UlJSvfj4+PwZM2Zkd+jQ4VxiYqLLJHrO0RMRkdOqV6+ePvHEE1nx8fFRgYGBF9q3b3/e/v1//vOfxyZNmuQ+fPjwtkuXLj1ga+/du3fesGHDTsbGxsYEBARciI+PPwsABQUFMmrUqLZWq9VdVWXs2LHHmjVrdvmJJ55ovXnz5gZubm4aHh6eP3LkyNySYhIR56/vbscl69Gf+OADeMd2QP0e3a+0nftpC84n7kHTMWMcESIRkdNx5Xr0NWXjxo2+kydPDvr555/TajoWe3WuHr13bAccmTQJ537aAsBI8kcmTYJ3bIcajoyIiGqrDRs2+I4ePTp0woQJlZqOqG4uOXRfv0d3BLz+Oo5MmoTG99yNU59+hoDXX7/qDJ+IiKgkW7du9bn//vvb2rd5eXkVZmZmJtZUTBXlkokeMJJ943vuxvG33kaz8Y8yyRMRUZnFx8fnp6amVnvteEdwyaF7wBiuP/XpZ2g2/lGc+vSzK8P4REREdYlLJnrbnHzA66/D/7HHrgzjM9kTEVFd45KJ/nzinqvm5G1z9ucT99RwZERERNXLJRN90zFjrpmTr9+jO39aR0TkIJs+m9siY9vWq253m7Ftq2XTZ3NZj76GuWSiJyKi6tUqLDJv5X9nhtqSfca2rZaV/50Z2iosstL16J9//vnmYWFhMe3bt4957rnnmgPAsWPH3Hv27BkWHBwc27Nnz7CcnJxratCXpnPnzpGVja02YKInIqJKa9c13jroL0/sX/nfmaFrP36v9cr/zgwd9Jcn9rfrGl+p6nU///yz99y5c/23b9+ekpKSkrRq1apGe/bsqfePf/yjVd++fa0HDx5M7Nu3r3XatGkty9pnQUEBAGDHjh2plYmttmCiJyKiKtGua7w15qZbcravXNYq5qZbciqb5AFgz549Pl26dDlrsVgKPT090atXL+vChQsbrVq1qtHYsWNPAMDYsWNPrFy5sjEATJ48ufUdd9zRtkePHuHBwcGxM2fObAYAX3/9taV79+7hQ4YMaRsREREDAL6+vp1t791www0RgwcPDg0JCYkdP358wNtvv92kQ4cOUeHh4dFJSUn1AODo0aMeAwYMaBcbGxsVGxsb9d1339UvKe6jR4969OzZMyw6Ojpq1KhRwa1bt+6QlZVVIz9pZ6InIqIqkbFtqyVpw/f+XQYNzUra8L1/0Tn7iujUqVP+li1bLNnZ2e5Wq9Vt9erVDQ8dOuR14sQJj+Dg4EsAEBwcfOnkyZNXkmhKSorPmjVr0n/66afUV199tXVmZqYnAOzevbv+q6++eiQjIyOp6HZSU1N93n777UMpKSlJX375ZdO9e/d679mzJ2X06NHHZ86c2RwAxo4dGzR58uRjiYmJKUuWLMkYN25cSElxP/3006379OljTU5OThk+fPiprKwsr8oei4py2RvmEBFR9bHNyduG69t06GStiuH7Ll26nH/88cez+/XrF+7r61sYHR2d5+Fx/dQ1aNCg035+furn51dw4403ntm4cWP9xo0bX+7YseO5yMjIi8Wt06FDh3O2Lw5t2rS5MGjQoFwAiIuLy1+/fr0FAH788ccG6enpPrZ1zp49637q1Cm3xo0bFxbtb+vWrX5Lly7dBwAjR44806BBg8sVPQaV5bAzehEJEpG1IpIiIkki8rjZ/qqIpIrIbhFZIiKNHBUDERFVj6z0VF/7pG6bs89KT61UPXoAmDRp0vHk5OSUhISEtCZNmlwOCws737Rp04KDBw96AsDBgwc9mzRpUmBbvmgdettrX1/faxKyTb169a5UeHNzc4O3t7fanl++fFkAQFWRkJCQkpqampyampr866+/7i4uyduWdRaOHLovAPCEqkYB6AHgLyISDWA1gFhV7QhgL4BnHBgDERFVg95333+s6Jl7u67x1t5331/pAjBHjhzxAID09HSvFStWNHrooYdODhgw4PS7777bFADefffdpgMHDjxtW37lypWN8vLyJDs72/2nn36y9O7d+1xlYwCA3r17n3n55Zeb215v3rzZp6Rl4+Pjz86bN68JACxevLjBmTNnyv2rgKrisKF7Vc0CkGU+t4pICoAAVf3ObrGfAIx0VAxERFT7DR06tN3p06c9PDw89I033vjF39//8r/+9a+sYcOGtQsODm7WunXri0uXLs2wLd+5c+dzt9xyS9jRo0e9pkyZkhUSEnIpMTHRu7JxvPfee4fGjBnTJjw8PPry5cvSvXt3a8+ePX8pbtmXXnrp6MiRI0Ojo6Mb33jjjWf9/f0vNWrUqEaG76ulHr2IhADYAONM/oxd+3IAC1X1k2LWeQTAIwDQpk2brgcPHnR4nERErqQu1qOfPHlyaz8/v8vPPfdcjZaSzc/PFw8PD/X09MSaNWvqT5gwIdiRRXKuV4/e4RfjiYgfgEUAJhZJ8lNhDO/PL249VX0PwHsA0K1bN+eZ7CAiIirFvn37vO688852hYWF8PT01HfffTezpmJxaKIXEU8YSX6+qi62a/8TgNsB3KLOdMUCERHVarNmzTpandv797//3fTtt9++6ja/N9xww9l58+b9kpKS4hRlbsuU6EWkHoARAELs11HV566zjgD4EECKqs6yax8I4CkAfVS10rdGJCIiqimPP/74iccff/xETcdxPWU9o/8KQC6AbQAulHGdXgBGA9gjIjvNtmcBzAZQD8Bq8ycPP6nquDJHTERERGVW1kQfqKoDy9Oxqm4CIMW89U15+iEiIqKKK+vv6DeLSAeHRkJERERV7rpn9CKyB4Cayz0gIvthDN0LADVvekNERHVc7reZLbzaWPJ8oppeuWlOfsoJy8VfrL4NB4TU6E/d6rrSzuhvBzAEwCAA7QH0N1/b2omIiODVxpJ38vO9ofkpJyyAkeRPfr431KuNhfXoa9h1E72qHlTVgwCm257bt1VPiERE5Ox8oppam9wZvv/k53tDTy/PaH3y872hTe4M329/hl8RrlCP/tKlS9WxmRKVdY4+xv6FiLgD6Fr14RARUW3lE9XUWr9L85yzPx5tVb9L85zKJnmg9tajnzx5cut77rknuFevXmHDhw9vW9njUBmlzdE/A+MncT4icga/XUV/EeZd64iIiABjuP7c9l/9/Xq1zjq3/Vf/eu0bWSub7Dt16pT/3HPPBWRnZ7vXr19fV69e3TAuLu5cafXot23blmK1Wt07d+4cPWLEiFzAqEe/Y8eOpOJK1aampvp8+eWX+5s3b14QHBzcoV69esf37NmT8vzzzzefOXNm8zlz5hyy1aMfMGDA2fT0dK8BAwaE7d+//5ra9ja7d+/23bJlS6qfn1+N3hjuuoleVV8E8KKIvKiqrDJHRETFss3J24br67VvZK2K4fvaWo8eAAYOHHi6ppM8UPbf0T8rIsMB9IZxFf5GVV3quLCIiKg2ufiL1dc+qdvm7C/+YvWt7Fn9pEmTjk+aNOk4AEyYMCEgMDDwoq0efXBw8KXqrkdf1uRdv379ErdXnco6R/9fAOMA7AGQCGCciPzXYVEREVGt0nBAyLGiCd0nqqm1Kn5aVxvr0TuTsp7R94FRYlYBQET+D0bSJyIicqjaWI/emZSpHr2ILAYwyfxZHUQkGMBLqnqPg+MDYJSpTUhIqI5NERG5DNajrzuqoh59UwApIrLVfH0DgP+JyDIAUNWhlY6SiIiIqlxZE/00h0ZBRERUBZypHn11xnE9ZUr0qrreHK4PU9U1IuIDwENVK30zBCIiotqqNtSjL9NV9yLyMIAvAbxrNgUC4M/riIiInFxZf173FwC9AJwBAFVNB9D8umsQERFRjStror+gqlfuJiQiHjBunENEREROrKyJfr2I2O55fxuALwAsd1xYRERUm3z//fct0tLSLPZtaWlplu+//75FSes4G1ctW1vWRP80gBwYN8kZC+AbAH9zVFBERFS7BAYG5i1ZsiTUluzT0tIsS5YsCQ0MDKx0PfrqUl1la6tbmRK9qhbCuPhuvKqOVNX3tSx32iEiojohIiLCOmzYsP1LliwJXblyZeslS5aEDhs2bH9ERESlfp2Vlpbm1bZt25i77rorOCwsLGbo0KFtly5daunSpUtkcHBw7Nq1a32PHTvmfuutt7YLDw+PjouLi9yyZYsPAKxYscIvMjIyOjIyMjoqKir61KlTbrm5uW433nhjeHR0dFR4eHj0J5980si2LVvZWgD429/+1iI8PDw6IiIievz48QHFxZaZmelp6z8yMjLa3d296969e70qs7+OUFqZWgHwDwATYJSoFRG5DOA/qvpcNcRHRES1REREhDUuLi5ny5Ytrbp3755V2SRvc+jQIe+FCxfu79q168GOHTtGzZ8/v2lCQkLqggULGs2YMaNVQEDAxbi4uLw1a9ZkLFu2zPKnP/2pbWpqavLMmTNbzp49+2D//v3P5ebmutmK2qxYsWJfkyZNCrOysjy6d+8eOWrUqNNubr+d937++ecNVqxY0Xjbtm2pFoul8NixY+7FxRUSEnIpNTU1GQBefPFF/40bN1rCw8OLrY5Xk0o7o58I42r7G1S1qao2AdAdQC8RmeTw6IiIqNZIS0uz7Nq1y7979+5Zu3bt8i86Z19RAQEBF+Lj4/Pd3d0RHh6e369fvzNubm7o0qVL3uHDh+tt3brV8tBDD50AgKFDh1pPnz7tceLECfcePXqcnTJlStD06dObHz9+3N3T0xOFhYUyceLEwPDw8Oibb745/Ndff/U6fPjwVSe9q1evbnDfffcdt1gshQDQokWLy9eL77vvvqs/d+5c/08//TSzKva3qpWW6O8HcI+qHrA1qOp+APeZ7xEREV2Zkx82bNj+QYMGHbUN41dFsvfy8iq2hKy7uzsuX75cbM0WEdEXXngh+4MPPjiYn5/v1rNnz6gdO3Z4v/vuu01OnDjhsWfPnpTU1NTkpk2bXsrPz78qF6rqNaVuS3Lw4EHPsWPHhixcuDCjYcOGTlGWtqjSEr2nql5T0EBVcwB4OiYkIiKqbQ4fPuxrPydvm7M/fPiwr6O33aNHD+tHH33UFAC+/vprS+PGjQuaNGlSmJSUVC8+Pj5/xowZ2R06dDiXmJjonZub696sWbNL9erV0+XLl1uOHj16zZz6wIEDz8ybN6+Z1Wp1A4CShu4vXLggw4cPD33++eePdOzY8YJj97LiSrsF7vXmGpxuHoKIiGrGLbfcck21uIiICGtVzdNfz8svv3x01KhRIeHh4dE+Pj6FH3/88QEAeOWVV5pv3ry5gZubm4aHh+ePHDky9/Tp0+6DBg1qHxsbGxUTE5PXtm3b80X7Gzly5Jnt27f7durUKcrT01NvvfXW3DfffPNI0eXWrFlTPzExsf706dNbT58+vTUArFq1Kj0kJOSSo/e5PK5bpta88O5ccW8B8FbVajmrZ5laIqLyq4tlauuqCpepVdVihyuIiIiodihrmVoiIqI6a/To0W1+/vlnP/u2Rx999JizV64DmOiJiIhK5Uz15currLfAJSKiuqewsLCwbL8zoxpjfkYl/rSPiZ6IiEqSmJOT05DJ3nkVFhZKTk5OQwCJJS3DoXsiIipWQUHBmOzs7A+ys7NjwRNDZ1UIILGgoGBMSQsw0RMRUbG6du36K4ChNR0HVQ6/oREREbkwhyV6EQkSkbUikiIiSSLyuNn+R/N1oYhU+EYOREREVDpHDt0XAHhCVbeLiAXANhFZDeOCgeEA3nXgtomIiAgOTPSqmgUgy3xuFZEUAAGquhpAmSsDERERUcVVyxy9iIQA6AxgSznWeUREEkQkIScnx1GhERERuTSHJ3oR8QOwCMBEVT1T1vVU9T1V7aaq3fz9/R0XIBERkQtzaKIXEU8YSX6+qi525LaIiIjoWo686l4AfAggRVVnOWo7REREVDJHXnXfC8BoAHtEZKfZ9iyAegD+A8AfwAoR2amqAxwYBxERUZ3lyKvuNwEo6dL6JY7aLhEREf2Gd8YjIiJyYUz0RERELoyJnoiIyIUx0RMREbkwJnoiIiIXxkRPRETkwpjoiYiIXBgTPRERkQtjoiciInJhTPREREQujImeiIjIhTHRExERuTAmeiIiIhfGRE9EROTCmOiJiIhcGBM9ERGRC2OiJyIicmFM9ERERC6MiZ6IiMiFMdETERG5MCZ6IiIiF8ZET0RE5MKY6ImIiFwYEz0REZELY6InIiJyYUz0RERELoyJnoiIyIUx0RMREbkwJnoiIifzzvoMbM44flXb5ozjeGd9Rg1FRLUZEz0RkZPpGNgQExbsuJLsN2ccx4QFO9AxsGENR0a1kUdNB0BERFfr2a4Z3hzVGRMW7MB93dvgky2/4M1RndGzXbOaDo1qIZ7RExE5oZ7tmuG+7m0w+4d9uK97GyZ5qjAmeiIiJ7Q54zg+2fILHuvXHp9s+eWaOXuismKiJyJyMrY5+TdHdcbk/hFXhvGZ7KkimOiJiJzM7sO5V83J2+bsdx/OreHIqDYSVXVMxyJBAOYCaAmgEMB7qvpvEWkCYCGAEACZAO5U1VPX66tbt26akJDgkDiJiFyViGxT1W41HQfVLEee0RcAeEJVowD0APAXEYkG8DSA71U1DMD35msiIiJyAIclelXNUtXt5nMrgBQAAQD+AOD/zMX+D8AdjoqBiIiorquWOXoRCQHQGcAWAC1UNQswvgwAaF7COo+ISIKIJOTk5FRHmERERC7H4YleRPwALAIwUVXPlHU9VX1PVbupajd/f3/HBUhEROTCHJroRcQTRpKfr6qLzeZjItLKfL8VgF8dGQMREVFd5rBELyIC4EMAKao6y+6tZQD+ZD7/E4CvHBUDERFRXefIe933AjAawB4R2Wm2PQvgJQCfi8hDAH4B8EcHxkBERFSnOSzRq+omAFLC27c4artERET0G94Zj4iIyIUx0RMREbkwJnoiIiIXxkRPRETkwpjoiYiIXBgTPRERkQtjoiciInJhTPREREQujImeiIjIhTHRExERuTAmeiIiIhfGRE9EROTCmOiJiIhcGBM9ERGRC2OiJyIicmFM9ERERC6MiZ6IiMiFMdETERG5MCZ6IiIns/3bgzicduqqtsNpp7D924M1FBHVZkz0REROpnlIA3z7fuKVZH847RS+fT8RzUMa1HBkVBt51HQARER0tcCIxhjwcCy+fT8RsTcFIHHDEQx4OBaBEY1rOjSqhXhGT0TkhAIjGiP2pgAkfJOJ2JsCmOSpwpjoiYic0OG0U0jccATdBocgccORa+bsicqKiZ6IyMnY5uQHPByL7kNDrwzjM9lTRTDRExE5mV8zz1w1J2+bs/8180wNR0a1ES/GIyJyMl0GBF/TFhjRmPP0VCE8oyciInJhTPREREQujImeiIjIhTHRExERuTAmeiIiIhcmqlrTMZRKRHIAVKaaQzMAx6sonOpUW+MGam/stTVuoPbGzrgdJ1hV/Ws6CKpZtSLRV5aIJKhqt5qOo7xqa9xA7Y29tsYN1N7YGTeRY3HonoiIyIUx0RMREbmwupLo36vpACqotsYN1N7Ya2vcQO2NnXETOVCdmKMnIiKqq+rKGT0REVGdxERPRETkwmpFoheRgSKSJiL7ROTpYt4PFpHvRWS3iKwTkUC7914RkSQRSRGR2WLwFZEVIpJqvveS3fJ/FpEcEdlpPsY4S9xm+zqzT1t8zc32eiKy0NzWFhEJqWjcjohdRCx2Me8UkeMi8oa5vLMc85dFJNF83GXX3tY8punmMfYy253pmJcU+3yzz0QRmSMinmZ7XxHJtTvm05ws7o9F5IBdfJ3MdjH/pvaZ/XWpaNwOjH2jXdxHRWSp2V5lx5yoXFTVqR8A3AFkAAgF4AVgF4DoIst8AeBP5vN+AOaZz3sC+NHswx3A/wD0BeAL4GZzGS8AGwEMMl//GcCbzhi3+d46AN2K2d54AO+Yz+8GsNDZYi+y/jYANznRMf89gNUwSjfXB5AAoIH53ucA7jafvwPgUSc75teLfTAAMR+f2sXeF8DXTnzMPwYwspjtDQaw0tyfHgC2OFvsRdZfBOD+qjzmfPBR3kdtOKOPB7BPVfer6kUAnwH4Q5FlogF8bz5fa/e+AvCG8T9xPQCeAI6pap6qrgUAs8/tAAJRtao87lK29wcA/2c+/xLALSLGKICzxS4iYQCaw/iCVZUqE3c0gPWqWqCq52D8oz/QPIb9YBxTwDjGd5jPneWYFxs7AKjqN2oCsBXO9XdeYtzX8QcAc81d+glAIxFp5Yyxi4gFxt/O0grGR1QlakOiDwBwyO71YbPN3i4AI8znwwBYRKSpqv4Pxv+cWebjW1VNsV9RRBoBGILf/mcGgBHmUN2XIhLkhHF/ZA79/d0usVzZnqoWAMgF0NQJYweAe2Cc/dr/5KNGj7nZPkiMaZ1mAG4GEATjGJ42j2nRPp3imF8n9ivMIfvRAFbZNd8oIrtEZKWIxDhh3DPMv4nXRaReObbnDLHblv9eVc/YtVXFMScql9qQ6Is7Qyr6m8ApAPqIyA4AfQAcAVAgIu0BRME4iwkA0E9EbrrSsYgHjOHM2aq632xeDiBEVTsCWIPfzticJe57VbUDgN+Zj9Hl2F5Nx25zN4zjblPjx1xVvwPwDYDNZmz/A1BQSp9OccyvE7u9twBsUFXbKMp2GPdBjwPwH1T8rNNRcT8DIBLADQCaAHiqHNur6dht7sHVf+dVdcyJyqU2JPrDuPqbciCAo/YLqOpRVR2uqp0BTDXbcmF8o/5JVc+q6lkYc3s97FZ9D0C6qr5h19cJVb1gvnwfQFdniltVj5j/tQJYAGP48artmV9gGgI46Uyxm7HFAfBQ1W12fTnDMYeqzlDVTqp6G4wkkA6jaEkj85gW7dNZjnlJscOM7R8A/AFMtuvrjPn5QFW/AeBpnpk6RdyqmmUOz18A8BGK+TsvaXs1HTsAmGf98QBW2PVVVcecqFxqQ6L/GUCYGFc+e8E4G1xmv4CINBMR2748A2CO+fwXGN/GPcyhyz4AUsx1psP4h3likb7s5/uG2pZ3hrjN183MdT0B3A4g0VxnGYA/mc9HAvihyNB4jcZut2rRsxynOOYi4m7+4wwR6QigI4DvzGO4FsYxBYxj/JX53CmOeUmxm6/HABgA4B5VLbTrq6Vt2kdE4mH8W3DCieJuZf5XYFwTYf93fr8YegDIVdWsCsTtsNhNf4Rx4d15u76q6pgTlY86wRWBpT1gXGm7F8YVslPNtucADDWfj4TxbXovgA8A1DPb3QG8CyNxJAOYZbYHwhiiSwGw03yMMd97EUASjDm4tQAinSju+jCuVt9txvhvAO7me94wrhDeB+Oiq1BnOuZ2/e4vekyd5Jh7m/EmA/gJQCe7PkPNY7rPPMb26zjDMb9e7AVmf7a/82lm+wS7Y/4TgJ5OFvcPAPbASPCfAPAz2wXAf81t7UExv0Cp6djN99cBGFikrcqOOR98lOfBW+ASERG5sNowdE9EREQVxERPRETkwpjoiYiIXBgTPRERkQtjoiciInJhTPRUJ5i/Yf5MRDJEJFlEvhGR8JqOi4jI0ZjoyeWZNylZAmCdqrZT1WgAzwJo4YBtuVd1n0RElcFET3XBzQAuqeo7tgZV3Qlgk4i8KkY98T1i1hQXo8b8YNuyYtRGH2HeDe1VEfnZLLYy1ny/r4isFZEFMG7iAhFZKiLbRCRJRB6x6+shEdkrRm3z90XkTbPdX0QWmX3/LCK9quXIEJHL8yh9EaJaLxbGHQWLGg6gE4A4AM0A/CwiG2CUK70LwDfmrVFvAfAogIdg3HL1BjGqqf0oIrbbnsYDiFXVA+brB1X1pIj4mP0uglG29+8AugCwwrj72y5z+X8DeF1VN4lIGwDfwigORERUKUz0VJf1BvCpql4GcExE1sOolrYSwGwzmQ+EUfUtX0T6A+goIrb73jcEEAbgIoCtdkkeAB4TkWHm8yBzuZYwapifBAAR+QKA7TqBWwFEy2/l7BuIiEWN4kVERBXGRE91QRJ+K0pjr7gypVDV8yKyDkYxmLvwWxEeAfD/VPXbqzoR6QvgXJHXtwK4UVXzzL68S9qeyc1cPr/03SEiKjvO0VNd8AOAeiLysK1BRG4AcArAXebcuz+Am2AUpwGM4fsHAPwOxjA6zP8+alblg4iEi0j9YrbXEMApM8lH4rcyvVthVPZrLEZZ2xF263wHo+iJLb5OldpjIiITz+jJ5amqmsPob4jI0wDOA8iEUaLYD8Y8uQL4q6pmm6t9B2AugGWqetFs+wBACIDt5pX8OTBKqBa1CsA4EdkNIA1GpTKo6hEReQHAFhh1z5MB5JrrPAbgv+Y6HgA2ABhXJQeAiOo0Vq8jqkYi4qeqZ80z+iUA5qjqkpqOi4hcF4fuiarXP0VkJ4w66wcALK3heIjIxfGMnoiIyIXxjJ6IiMiFMdETERG5MCZ6IiIiF8ZET0RE5MKY6ImIiFzY/welvb0oO1f4lQAAAABJRU5ErkJggg==\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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}