{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# SPIRE-NEP Selection Functions\n",
    "## Depth maps and selection functions\n",
    "\n",
    "The simplest selection function available is the field MOC which specifies the area for which there is Herschel data. Each pristine catalogue also has a MOC defining the area for which that data is available.\n",
    "\n",
    "The next stage is to provide mean flux standard deviations which act as a proxy for the catalogue's 5$\\sigma$ depth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This notebook was run with herschelhelp_internal version: \n",
      "017bb1e (Mon Jun 18 14:58:59 2018 +0100)\n",
      "This notebook was executed on: \n",
      "2018-06-26 11:39:26.404448\n"
     ]
    }
   ],
   "source": [
    "from herschelhelp_internal import git_version\n",
    "print(\"This notebook was run with herschelhelp_internal version: \\n{}\".format(git_version()))\n",
    "import datetime\n",
    "print(\"This notebook was executed on: \\n{}\".format(datetime.datetime.now()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "#%config InlineBackend.figure_format = 'svg'\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rc('figure', figsize=(10, 6))\n",
    "\n",
    "import os\n",
    "import time\n",
    "\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "from astropy.table import Column, Table, join\n",
    "import numpy as np\n",
    "from pymoc import MOC\n",
    "import healpy as hp\n",
    "#import pandas as pd #Astropy has group_by function so apandas isn't required.\n",
    "import seaborn as sns\n",
    "\n",
    "import warnings\n",
    "#We ignore warnings - this is a little dangerous but a huge number of warnings are generated by empty cells later\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "from herschelhelp_internal.utils import inMoc, coords_to_hpidx, flux_to_mag\n",
    "from herschelhelp_internal.masterlist import find_last_ml_suffix, nb_ccplots\n",
    "\n",
    "from astropy.io.votable import parse_single_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "FIELD = 'SPIRE-NEP'\n",
    "#FILTERS_DIR = \"/Users/rs548/GitHub/herschelhelp_python/database_builder/filters/\"\n",
    "FILTERS_DIR = \"/opt/herschelhelp_python/database_builder/filters/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Depth maps produced using: master_catalogue_spire-nep_20180220.fits\n"
     ]
    }
   ],
   "source": [
    "TMP_DIR = os.environ.get('TMP_DIR', \"./data_tmp\")\n",
    "OUT_DIR = os.environ.get('OUT_DIR', \"./data\")\n",
    "SUFFIX = find_last_ml_suffix()\n",
    "#SUFFIX = \"20171016\"\n",
    "\n",
    "master_catalogue_filename = \"master_catalogue_{}_{}.fits\".format(FIELD.lower(), SUFFIX)\n",
    "master_catalogue = Table.read(\"{}/{}\".format(OUT_DIR, master_catalogue_filename))\n",
    "\n",
    "print(\"Depth maps produced using: {}\".format(master_catalogue_filename))\n",
    "\n",
    "ORDER = 10\n",
    "#TODO write code to decide on appropriate order\n",
    "\n",
    "field_moc = MOC(filename=\"../../dmu2/dmu2_field_coverages/{}_MOC.fits\".format(FIELD))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Remove sources whose signal to noise ratio is less than five as these will have been selected using forced \n",
    "# photometry and so the errors will not refelct the RMS of the map \n",
    "for n,col in enumerate(master_catalogue.colnames):\n",
    "    if col.startswith(\"f_\"):\n",
    "        err_col = \"ferr{}\".format(col[1:])\n",
    "        errs = master_catalogue[err_col]\n",
    "        fluxes = master_catalogue[col]\n",
    "        mask = fluxes/errs < 5.0\n",
    "        master_catalogue[col][mask] = np.nan\n",
    "        master_catalogue[err_col][mask] = np.nan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## I - Group masterlist objects by healpix cell and calculate depths\n",
    "We add a column to the masterlist catalogue for the target order healpix cell <i>per object</i>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#Add a column to the catalogue with the order=ORDER hp_idx\n",
    "master_catalogue.add_column(Column(data=coords_to_hpidx(master_catalogue['ra'],\n",
    "                                                       master_catalogue['dec'],\n",
    "                                                       ORDER), \n",
    "                                   name=\"hp_idx_O_{}\".format(str(ORDER))\n",
    "                                  )\n",
    "                           )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Convert catalogue to pandas and group by the order=ORDER pixel\n",
    "\n",
    "group = master_catalogue.group_by([\"hp_idx_O_{}\".format(str(ORDER))])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#Downgrade the groups from order=ORDER to order=13 and then fill out the appropriate cells\n",
    "#hp.pixelfunc.ud_grade([2599293, 2599294], nside_out=hp.order2nside(13))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## II Create a table of all Order=13 healpix cells in the field and populate it\n",
    "We create a table with every order=13 healpix cell in the field MOC. We then calculate the healpix cell at lower order that the order=13 cell is in. We then fill in the depth at every order=13 cell as calculated for the lower order cell that that the order=13 cell is inside."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "depths = Table()\n",
    "depths['hp_idx_O_13'] = list(field_moc.flattened(13))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=10</i>\n",
       "<table id=\"table139938461869952-818608\" 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>190654398</td></tr>\n",
       "<tr><td>1</td><td>190678020</td></tr>\n",
       "<tr><td>2</td><td>190678021</td></tr>\n",
       "<tr><td>3</td><td>190678022</td></tr>\n",
       "<tr><td>4</td><td>190678023</td></tr>\n",
       "<tr><td>5</td><td>190678028</td></tr>\n",
       "<tr><td>6</td><td>190678029</td></tr>\n",
       "<tr><td>7</td><td>190678030</td></tr>\n",
       "<tr><td>8</td><td>190678031</td></tr>\n",
       "<tr><td>9</td><td>190678032</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(\"$('#table139938461869952-818608').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",
       "    $('#table139938461869952-818608').dataTable({\n",
       "        order: [],\n",
       "        pageLength: 50,\n",
       "        lengthMenu: [[10, 25, 50, 100, 500, 1000, -1], [10, 25, 50, 100, 500, 1000, 'All']],\n",
       "        pagingType: \"full_numbers\",\n",
       "        columnDefs: [{targets: [0, 1], type: \"optionalnum\"}]\n",
       "    });\n",
       "});\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "depths[:10].show_in_notebook()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "depths.add_column(Column(data=hp.pixelfunc.ang2pix(2**ORDER,\n",
    "                     hp.pixelfunc.pix2ang(2**13, depths['hp_idx_O_13'], nest=True)[0],\n",
    "                     hp.pixelfunc.pix2ang(2**13, depths['hp_idx_O_13'], nest=True)[1],\n",
    "                     nest = True),\n",
    "                     name=\"hp_idx_O_{}\".format(str(ORDER))\n",
    "                        )\n",
    "                 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=10</i>\n",
       "<table id=\"table139938461869392-379929\" 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>190654398</td><td>2978974</td></tr>\n",
       "<tr><td>1</td><td>190678020</td><td>2979344</td></tr>\n",
       "<tr><td>2</td><td>190678021</td><td>2979344</td></tr>\n",
       "<tr><td>3</td><td>190678022</td><td>2979344</td></tr>\n",
       "<tr><td>4</td><td>190678023</td><td>2979344</td></tr>\n",
       "<tr><td>5</td><td>190678028</td><td>2979344</td></tr>\n",
       "<tr><td>6</td><td>190678029</td><td>2979344</td></tr>\n",
       "<tr><td>7</td><td>190678030</td><td>2979344</td></tr>\n",
       "<tr><td>8</td><td>190678031</td><td>2979344</td></tr>\n",
       "<tr><td>9</td><td>190678032</td><td>2979344</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(\"$('#table139938461869392-379929').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",
       "    $('#table139938461869392-379929').dataTable({\n",
       "        order: [],\n",
       "        pageLength: 50,\n",
       "        lengthMenu: [[10, 25, 50, 100, 500, 1000, -1], [10, 25, 50, 100, 500, 1000, 'All']],\n",
       "        pagingType: \"full_numbers\",\n",
       "        columnDefs: [{targets: [0, 1, 2], type: \"optionalnum\"}]\n",
       "    });\n",
       "});\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "depths[:10].show_in_notebook()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=10</i>\n",
       "<table id=\"table139939698934728-669377\" 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_gpc1_g_mean</th><th>f_ap_gpc1_g_p90</th><th>ferr_gpc1_g_mean</th><th>f_gpc1_g_p90</th><th>ferr_ap_gpc1_r_mean</th><th>f_ap_gpc1_r_p90</th><th>ferr_gpc1_r_mean</th><th>f_gpc1_r_p90</th><th>ferr_ap_gpc1_i_mean</th><th>f_ap_gpc1_i_p90</th><th>ferr_gpc1_i_mean</th><th>f_gpc1_i_p90</th><th>ferr_ap_gpc1_z_mean</th><th>f_ap_gpc1_z_p90</th><th>ferr_gpc1_z_mean</th><th>f_gpc1_z_p90</th><th>ferr_ap_gpc1_y_mean</th><th>f_ap_gpc1_y_p90</th><th>ferr_gpc1_y_mean</th><th>f_gpc1_y_p90</th></tr></thead>\n",
       "<tr><td>0</td><td>190654398</td><td>2978974</td><td>0.8508792947319282</td><td>5.987972755250806</td><td>0.8725040622072272</td><td>5.809788316065019</td><td>1.2234321270806152</td><td>30.299190007362277</td><td>1.248650958468424</td><td>28.14446545789632</td><td>1.181045124443541</td><td>60.989906475001305</td><td>1.260240301825581</td><td>55.96301407894633</td><td>3.8086216866115024</td><td>84.58289644434113</td><td>2.9351938366156145</td><td>74.56950616219889</td><td>6.336414850332582</td><td>104.77101905631234</td><td>7.799786297185613</td><td>91.98493050835323</td></tr>\n",
       "<tr><td>1</td><td>190654396</td><td>2978974</td><td>0.8508792947319282</td><td>5.987972755250806</td><td>0.8725040622072272</td><td>5.809788316065019</td><td>1.2234321270806152</td><td>30.299190007362277</td><td>1.248650958468424</td><td>28.14446545789632</td><td>1.181045124443541</td><td>60.989906475001305</td><td>1.260240301825581</td><td>55.96301407894633</td><td>3.8086216866115024</td><td>84.58289644434113</td><td>2.9351938366156145</td><td>74.56950616219889</td><td>6.336414850332582</td><td>104.77101905631234</td><td>7.799786297185613</td><td>91.98493050835323</td></tr>\n",
       "<tr><td>2</td><td>190654399</td><td>2978974</td><td>0.8508792947319282</td><td>5.987972755250806</td><td>0.8725040622072272</td><td>5.809788316065019</td><td>1.2234321270806152</td><td>30.299190007362277</td><td>1.248650958468424</td><td>28.14446545789632</td><td>1.181045124443541</td><td>60.989906475001305</td><td>1.260240301825581</td><td>55.96301407894633</td><td>3.8086216866115024</td><td>84.58289644434113</td><td>2.9351938366156145</td><td>74.56950616219889</td><td>6.336414850332582</td><td>104.77101905631234</td><td>7.799786297185613</td><td>91.98493050835323</td></tr>\n",
       "<tr><td>3</td><td>190654397</td><td>2978974</td><td>0.8508792947319282</td><td>5.987972755250806</td><td>0.8725040622072272</td><td>5.809788316065019</td><td>1.2234321270806152</td><td>30.299190007362277</td><td>1.248650958468424</td><td>28.14446545789632</td><td>1.181045124443541</td><td>60.989906475001305</td><td>1.260240301825581</td><td>55.96301407894633</td><td>3.8086216866115024</td><td>84.58289644434113</td><td>2.9351938366156145</td><td>74.56950616219889</td><td>6.336414850332582</td><td>104.77101905631234</td><td>7.799786297185613</td><td>91.98493050835323</td></tr>\n",
       "<tr><td>4</td><td>190654429</td><td>2978975</td><td>3.7878850055096756</td><td>232.2456878511544</td><td>4.426276266109323</td><td>219.44388234258776</td><td>52.249254726880224</td><td>263.8486962726817</td><td>39.83238424600856</td><td>444.56029236611124</td><td>6.440277562451329</td><td>290.75032083762835</td><td>5.1880766752061005</td><td>340.97327183677675</td><td>15.119204809341662</td><td>330.94791233667166</td><td>67.50130814821412</td><td>471.8443930103042</td><td>136.02148970223428</td><td>1359.215446869019</td><td>151.29731701815248</td><td>5548.739115285172</td></tr>\n",
       "<tr><td>5</td><td>190654456</td><td>2978975</td><td>3.7878850055096756</td><td>232.2456878511544</td><td>4.426276266109323</td><td>219.44388234258776</td><td>52.249254726880224</td><td>263.8486962726817</td><td>39.83238424600856</td><td>444.56029236611124</td><td>6.440277562451329</td><td>290.75032083762835</td><td>5.1880766752061005</td><td>340.97327183677675</td><td>15.119204809341662</td><td>330.94791233667166</td><td>67.50130814821412</td><td>471.8443930103042</td><td>136.02148970223428</td><td>1359.215446869019</td><td>151.29731701815248</td><td>5548.739115285172</td></tr>\n",
       "<tr><td>6</td><td>190654424</td><td>2978975</td><td>3.7878850055096756</td><td>232.2456878511544</td><td>4.426276266109323</td><td>219.44388234258776</td><td>52.249254726880224</td><td>263.8486962726817</td><td>39.83238424600856</td><td>444.56029236611124</td><td>6.440277562451329</td><td>290.75032083762835</td><td>5.1880766752061005</td><td>340.97327183677675</td><td>15.119204809341662</td><td>330.94791233667166</td><td>67.50130814821412</td><td>471.8443930103042</td><td>136.02148970223428</td><td>1359.215446869019</td><td>151.29731701815248</td><td>5548.739115285172</td></tr>\n",
       "<tr><td>7</td><td>190654425</td><td>2978975</td><td>3.7878850055096756</td><td>232.2456878511544</td><td>4.426276266109323</td><td>219.44388234258776</td><td>52.249254726880224</td><td>263.8486962726817</td><td>39.83238424600856</td><td>444.56029236611124</td><td>6.440277562451329</td><td>290.75032083762835</td><td>5.1880766752061005</td><td>340.97327183677675</td><td>15.119204809341662</td><td>330.94791233667166</td><td>67.50130814821412</td><td>471.8443930103042</td><td>136.02148970223428</td><td>1359.215446869019</td><td>151.29731701815248</td><td>5548.739115285172</td></tr>\n",
       "<tr><td>8</td><td>190654431</td><td>2978975</td><td>3.7878850055096756</td><td>232.2456878511544</td><td>4.426276266109323</td><td>219.44388234258776</td><td>52.249254726880224</td><td>263.8486962726817</td><td>39.83238424600856</td><td>444.56029236611124</td><td>6.440277562451329</td><td>290.75032083762835</td><td>5.1880766752061005</td><td>340.97327183677675</td><td>15.119204809341662</td><td>330.94791233667166</td><td>67.50130814821412</td><td>471.8443930103042</td><td>136.02148970223428</td><td>1359.215446869019</td><td>151.29731701815248</td><td>5548.739115285172</td></tr>\n",
       "<tr><td>9</td><td>190654415</td><td>2978975</td><td>3.7878850055096756</td><td>232.2456878511544</td><td>4.426276266109323</td><td>219.44388234258776</td><td>52.249254726880224</td><td>263.8486962726817</td><td>39.83238424600856</td><td>444.56029236611124</td><td>6.440277562451329</td><td>290.75032083762835</td><td>5.1880766752061005</td><td>340.97327183677675</td><td>15.119204809341662</td><td>330.94791233667166</td><td>67.50130814821412</td><td>471.8443930103042</td><td>136.02148970223428</td><td>1359.215446869019</td><td>151.29731701815248</td><td>5548.739115285172</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(\"$('#table139939698934728-669377').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",
       "    $('#table139939698934728-669377').dataTable({\n",
       "        order: [],\n",
       "        pageLength: 50,\n",
       "        lengthMenu: [[10, 25, 50, 100, 500, 1000, -1], [10, 25, 50, 100, 500, 1000, 'All']],\n",
       "        pagingType: \"full_numbers\",\n",
       "        columnDefs: [{targets: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22], type: \"optionalnum\"}]\n",
       "    });\n",
       "});\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for col in master_catalogue.colnames:\n",
    "    if col.startswith(\"f_\"):\n",
    "        errcol = \"ferr{}\".format(col[1:])\n",
    "        depths = join(depths, \n",
    "                      group[\"hp_idx_O_{}\".format(str(ORDER)), errcol].groups.aggregate(np.nanmean),\n",
    "                     join_type='left')\n",
    "        depths[errcol].name = errcol + \"_mean\"\n",
    "        depths = join(depths, \n",
    "                      group[\"hp_idx_O_{}\".format(str(ORDER)), col].groups.aggregate(lambda x: np.nanpercentile(x, 90.)),\n",
    "                     join_type='left')\n",
    "        depths[col].name = col + \"_p90\"\n",
    "\n",
    "depths[:10].show_in_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## III - Save the depth map table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "depths.write(\"{}/depths_{}_{}.fits\".format(OUT_DIR, FIELD.lower(), SUFFIX), overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## IV - Overview plots\n",
    "\n",
    "### IV.a - Filters\n",
    "First we simply plot all the filters available on this field to give an overview of coverage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'gpc1_g', 'gpc1_i', 'gpc1_r', 'gpc1_y', 'gpc1_z'}"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tot_bands = [column[2:] for column in master_catalogue.colnames \n",
    "             if (column.startswith('f_') & ~column.startswith('f_ap_'))]\n",
    "ap_bands = [column[5:] for column in master_catalogue.colnames \n",
    "            if column.startswith('f_ap_') ]\n",
    "bands = set(tot_bands) | set(ap_bands)\n",
    "bands"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Passbands on SPIRE-NEP')"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdoAAAEgCAYAAAAT7bLQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXd4HOXVt+8zW7RqtmVb7t0Y94ZNN72ZYiAhAQOmGhwSIOElhPClvBB4Qw8koYYWCIQWSmgGEwhgQgBTDMY27rjI3bJcZJXdmXm+P2Z2tZLVNbuz0j73demStDP7zNm1V7855zlFlFJoNBqNRqNJDYbfBmg0Go1G05HRQqvRaDQaTQrRQqvRaDQaTQrRQqvRaDQaTQrRQqvRaDQaTQrRQqvRaDQaTQrRQqvJOETkBhF5qqNdS6PRZCdaaDXNQkRWi0iliJSLyGYR+auIFPhtV0dARGaKyBIR2e2+t2+ISKF77HERibrv+3YR+ZeIjHCP1bpJEBElInvcc9eLyF0iEkg6/r6IVLnH41+vNWLXatee/KTHLhGR9xu4Zvzr2iT7Yu5jO0TkvyJysKdvnkbTDtBCq2kJ05RSBcB+wP7Ab3y2p90jIkcANwNnK6UKgZHA83VOu9193/sBW4DHG1lyvHvuEcBZwMV1jl+hlCpI+prWhIlB4GdNnDO+zpq3Jx17zrWnGPgP8JKISBPraTQdCi20mhajlFoPvAmMARCRi0TkW9cjWyUiP4qfKyLdReR116PZLiIfiojhHvul63ntFpGlInJM0mUiIvKce+xLERmftOZ1IrLSPbZYRL6XdOxCEfmPiNwpImUi8p2InJh0fLCIfOA+919A96RjERF5SkRKXXs/E5Ge9b0HIjLS9RB3iMgiETk16djjInKf65nuFpFPRWRoA2/n/sDHSqn57nu7XSn1hFJqdz3vewXwdPx9bwyl1ArgI2BCU+c2wR3ANSLSpS2LKKViwBNAL6BbG23SaNoVWmg1LUZE+gMnAfPdh7YApwCdgIuAu0VkP/fYz4ESHI+mJ/ArQInIcOAKYH/XkzsBWJ10mdOAfwBdccTlnyISco+tBA4DOgO/A54Skd5Jzz0QWIojorcDjyZ5UU8DX7jHbgIuSHreBe6a/XHE4DKgsp7XHwJeA94GegBXAn93X1Ocs13bioAVwO/rruPyKXCCiPxORA4VkZwGzsMN1Z9LzfveIG54+TD32m3hc+B94Jq2LOK+rguBEqXUtjbapNG0K7TQalrCP0VkB04I8AOckCdKqTeUUiuVwwc4AnSY+5wY0BsYqJSKKaU+VE6DbQvIAUaJSEgptVoptTLpWl8opV5wPaG7gAhwkHu9fyilNiilbKXUc8By4ICk565RSj2slLJwvKjeQE8RGYDjQf5WKVWtlJqLI5hxYjgCu49SylJKfaGU2lXP+3AQUADcqpSKKqX+DbyOI65xXlJKzVNKmcDfacCzVEp9CHwfJxz/BlBad28Vx6PcgSOaBTiC1RBfisge4Fscgby/zvE/u154/OumRtaK87/AlSJS3Mg1k9c8IenYma7t64BJwOnNuJ5G06HQQqtpCacrpboopQYqpX6ilKoEEJETReQTNzS8A8fbjYdk78ARiLfdsPJ1kAhtXgXcAGwRkWdFpE/StdbFf1BK2ThecR/3eueLyFfxP+w4odTuSc/dlPTcCvfHAvf5ZUqpPUnnrkn6+UlgDvCsiGwQkduTvOhk+gDrXLuS1+lbnw1AhXv9elFKvenulXbF8eQvBC5JOuVO933vpZQ6tc4NSV32c691Fo5nn1/n+E/dteJfvwUQkTeTkpnOrWPfQpwbiesaumadNeckHXvefayHUupopdQXjdiu0XRItNBq2oQbEnwRuBPoqZTqAswGBEAptVsp9XOl1BBgGnB1fC9WKfW0UmoKMBBQwG1JS/dPuoaBkwi0QUQGAg/jhJ27uddbGL9eE2wEiiQpixYYEP/B9bh/p5QaBRyCEw4/v551NgD943vNSeusb4YNDeJ66O8C/6YZ+7CNrKOUUs8DH+N4o815zolJyUx/r+eU64FLqX0zodFomoEWWk1bCeOEgLcCppt4dHz8oIicIiL7uHuku3BCxpaIDBeRo12hrsLZC7WS1p0kIt8XkSCO51sNfILjoSn3eojIRTRTlJRSa3D2HH8nImERmYIj/nFbjxKRsW7YdhdOKNmqZ6lPgT3AtSISEpEj3XWebY4dyYjIaSIyXUSKxOEAnIzhT1q6Vj3cCswSkV5tXciNQDwH/LTNVmk0WYYWWk2bcLNjf4pTklIGnAO8mnTKMOAdoBzHw7pfKfU+jjjfCmzDCbP2wEmUivMKTvizDDgP+L7rcS4G/uCutRkYi5Nd21zOwQmpbsfx0v6WdKwX8AKOyH6Lsw+9VzMLpVQUOBU40bX/fuB8pdSSFtgRpwzHU1zuXvcp4I4GvMoWoZT6Buc1/CLp4Xulds1rS0K5N7J3KBrg6zpr/rEtdms0HQ3Rg981Go1Go0kd2qPVaDQajSaFaKHVaDQajSaFaKHVaDQajSaFaKHVaDQajSaFaKHVaDQajSaFBP02oKV0795dDRo0yG8zNBqNpl3xxRdfbFNKNdRGs7lr9AgGg4/g1K5rR83BBhaapnnJpEmTttR3QrsT2kGDBvH555/7bYZGo9G0K0RkTdNnNU4wGHykV69eI4uLi8sMw9C1oYBt27J169ZRmzZtegSnvn4v9B2JRqPRaJrLmOLi4l1aZGswDEMVFxfvpJEOdVpoNRqNRtNcDC2ye+O+Jw3qqRZajUaj0WhSiBZajUaj0WQVV155Zd9evXqNy8vLm9jUubfffnvxvffe260t12t3yVAajUaj0bSF008/fcc111yzZeTIkU1O/rr22mu3tvV6Wmg1Go1G02J+8cLX/Zdt2p3n5Zr79iqsuOMH49c1et1f/KL3Cy+80LV3797Rbt26mRMnTqx46623uowZM6Zi/vz5+eXl5YGHHnrou6OOOqpi586dxsyZMwcsWLAgD+BXv/rVhgsvvHDHMcccs6e5Nl199dV9CgoKrBtvvHFza1+XFlqNRpMWKnZFKV1fTr8RRTjjiTWaljF37ty81157reibb75ZHIvFZMKECaMmTpxYAVBRUWHMnz9/yZtvvlkwa9aswcuXL1903XXX9e7UqZO1bNmyxQBbt24N+GG3FlqNRpNyKnZFefp3n1C9x+T4maMZtn9Pv03StJGmPM9U8P777xeceOKJOwoKChSgjjvuuB3xY+ecc852gBNPPLG8vLzc2LZtW2Du3Lmdnn322VXxc4qLi6102ww6GUqj0aSBtYtLqd5jAvDl22uwLdtnizTtkcbmp9eNkogISqmMiJ5oodVoOihR06Yq5tzAz122lR8/9QWXPfkFby3c1OgfrFSwbN5mxBCOuXAk29aV8+xN81jy8ca026Fp3xx55JHlc+bM6VxRUSE7d+403nnnnS7xY88880wRwJw5cwoKCwutbt26WUceeeSuu+66q0f8HB061mg0LWL55t3sqjKZNLCo1uPbyqu59oUFfLRiGzHLJj8cZHe1SSTk3Fe/tWgTFx86mN+eMjJxt//h8q38/Pmv+ekxw5hx0EBP7Nu0aifhSJCcvCDrFm8ntzDEiIN6U73H5Kt31/LuE9+ydvF2jrt4VEZ4HZrM54gjjqiYOnXqzlGjRo3u27dv9bhx4/Z07tzZAigqKrImTpw4Ip4MBXDLLbdsvOiiiwYMGzZstGEY6le/+tWGCy64YMdll13W7+WXX+5aVVVl9OzZc9y555677a677tqQKrulvd1RTp48Welex5pspSpm8eh/vuPMyf3Z//fvALD61pNrnXPLm9/ylw9WMXV0L7oWhAkHDDrnhrj08CFEgga/fWUhz8xbxwPn7seJY3sDcPXzX/HSl+sJGMJbPzuMYT0L22zrfZf9u9bvx144kuEHOddTtuKz2av57PXvOOq8EYw6tE+br6dpHBH5Qik1uS1rfP3116vHjx+/zSubWsPOnTuNzp0727t37zYOPvjg4Q8++OCaq6++uv+dd9657vDDD6/wy66vv/66+/jx4wfVd0x7tBpNhrN5VxV54QAFOUFmPvEZH60oJSdYs+tj2wrDkMS5f/1oNadP6MMfp9dfi3/TaWN4Zt46lm0u58SxzmNfr9vB5IFFrNhaznUvfcOzsw4iFGj9zpKy976B79avRrzFEPY/eRBrF5Xy6aurGDi6G/ldclp9PU32MGPGjIHLly/Pra6ulunTp5dOmTLFN3FtLlpoNZoMw7YV1aZNMCDM+247F/51HgO65nH+wYP4aEUpAFt2VyfO37anmh6FEQD+8sEqLFvx8+OHN7i+4YZpFSpxvXVllRwzsifnHTyQnz37Fbe+uYTfnjKq1a+hbNPef/u69yuo9buIMOWHw3jlj/N57Z6v6DeyK9XlMcYd3Z/iAW33qDUdk9dee+27uo/NmzdvaVvX/eUvf9nrlVde6Zr82Gmnnbb9tttu29TWtbXQajQZxi1vfsvDH9b+W7Jy6x6uf3UR4/t1ZleVyZrSmnr79WWV9CiMUFpezdPz1nD6hL7079pwH4G626Fby6uJmjb9i3I5bUJf5n23nb9+9B1n7d+ffeuEkLfsruKx/6zmR4cPoSg/3OA1ls6r/bepU/dIvef1GtKZYy4YxZyHF1K63nlNSz7dxHEXjWLfA3o1uL5G4zW33XbbJi9EtT600Go0GYRlq4TIDuiah2nZ/OPHh7C2tII5izbxoyOGcN2L37CgZGfiOec/No9nLj2I77btoSpmc/7BjSczxROP4ukZ67Y73mc/V5yvOX44r329geteXMADMybRs1ONSN45ZynPf15CSVkF956zX4PX2F1aVev3CccOaPDcfSb1IJI/gWA4QEFRhDkPL+Sdx7+lqHc+xf21Z6tp/+jyHo0mgyiriALwu1NHM/fao/jwl0fTt0suBw/txg2njqZ351z6FeWycWeNkO2uMjnlnv9w5TPzEYGRvTs161rxXdS1rtAOcIW2KD/MTaePYdGGXZz98CfsrIwB8MmqUp7/vASA1xds5Jsksa/Lnh01oe1L/3g4Y4/s16gt/UZ0pdeQzhQU5XDy5eMIRwLMfmABVeWxZr0WjSaT0UKr0WQQO1yh7ZIXAiBg7F320lhYOGQYhINNf6xFSLi067ZXAtC3S27i+GkT+vLkzANZU1rB3f9aBsBfPlgJwIs/PoSu+WEu+dtnbNxZWe/6Fbui9B/Vle/9fCLhSMsCZ5H8EKdcMZ6KXVHeenihFltNu0cLrUaTQeyocESlS17D+5/9imoE8fQJfRjfP1GzTyjQvHpUocajLauIUhgJEgnVruU/YHBXpu/fn6c+WcOKLbv59LvtnH/wQCYNLOKvF+5PeZXJj578glg9XZ6qK00Ku0boM6xor2PNodeQzhx5zgg2rtjBG/cvwNKdpDTtGC20Gk0GUeYKbZHr0dZH/6Iaj/bm74/llcsP5evrj+fAwV156PzmlUk67emcn6ti1l4iG+eqY/clJ2hw7F1zqYhaHDLUGcs5vn8XbjljHAtKdjL7m417PU9ZCqOZot8QIw/pzdHnj2TTqp2s+KzVg1M0mr1oyTxaL9BCq9FkEPE92qJGPNrk0HFe2AnLds4N8dyPDubQfbo36zqOR+sobbVpJ7pG1aW4MIcLDx0EwIhehRw3qiYT+JSxvRncPZ/H/7t6r+fZlt1moQXY94CedOmZx4L317d5LY0mzumnn77j008//TZd19NZxxpNBlF3j7Y+GvN2m4tITdZxtWmRE2y4BewlU4ZQHbO59PAhtfaMDUO44OCB3PDaYr5at4MJSSFs21IYbWh4UWOnMObwvvznH8vZvmEPXfvkt3lNjUf88/L+bFns6Txaeoyq4PT7MmYebVlZmTFmzJjRq1atWpiTk6O2b99ujB07NvF7c19WSj1aEZkqIktFZIWIXFfP8QEi8p6IzBeRBSJyUirt0WgynbKKGEFDKMhp+B5YRLjn7Im8fuWUVl9HkMQebVWsYY8WnCzk35wyqlaZT5wzJvWjMCfIg++vrPW47UHoOM4+k5ye8N8t2OrJepr2S/I82jfeeGPlggULEnde8Xm0f/7zn9fMmjVrMEDyPNply5YtPvnkk3e35HpFRUX2wQcfvPv555/vDPDYY491Pemkk8paIrKQQo9WRALAfcBxQAnwmYi8qpRanHTab4DnlVIPiMgoYDYwKFU2aTSZzo6KGF3yQk022Z82vo29gVvg0TZGYSTEOQcN4OG5q9hZEaNzXgillNMW0iOhze+SQ/GAQlYvKGXS1EGerKnxgCY8z1TgxzzaWbNmbb3tttt6nXfeeTueeuqp7g8//PDqlq6RSo/2AGCFUmqVUioKPAucVuccBcSL/joDKZueoNG0B3ZURBvNOPaK5D3apjzapjhmRE9sBR+vctpD2m6f44BHQgswaGw3Nn23k4pdUc/W1LQ//JhHe/zxx+8pKSnJeeONNwosy5L999+/quln1SaVQtsXSL7jKXEfS+YGYIaIlOB4s1fWt5CIzBKRz0Xk861bdfhI03Epq4h6sgfbFJJU39MWjxZg4oAu5IcDfLTCGepiW87CXuzRxhkysRgUrF7g6+AYjc/4NY92+vTppRdddNGQGTNmtOo/YCqFtr7biLq3I2cDjyul+gEnAU+KyF42KaUeUkpNVkpNLi4uToGpGk1m4ISO0+HRNn+PtilCAYNhPQv5bpuTX6ISQuudR9utbwGF3SKsnK9vtLOZ5Hm0J5100tD65tFeccUVA//yl7+sBmce7Y4dOwLDhg0bPXz48FGzZ88uBLjsssv69ezZc1x8Hu3VV1/d6F7MzJkzS3ft2hWcOXPm9tbYncqs4xKgf9Lv/dg7NDwTmAqglPpYRCJAd2BLCu3SaDKW3VUmhS3spNQanKzjeHlP2zxagF6dIqzcWg4ke7TeCa2IsO8BPfnizTWsX1ZG331b1whD0/65/vrrN911110b4vNor7322s3PPfdct7POOqvsvvvuq1UH1rlzZ/ull15aXXeNBx98sOTBBx8sae4133333cKpU6eWde/evcV7vJBaj/YzYJiIDBaRMDAdeLXOOWuBYwBEZCQQAfQtqyZriVp2rVmzqSKpA2ObPVqAnp1y2LTL2bqKd3HyMnQMMGnqIPI6h3n1z1+xY3PGjyDVpIgZM2YMHDFixKhx48aNnDZtWlmq59FecMEF/a+//vq+N954Y6tziFJ266yUMkXkCmAOEAAeU0otEpEbgc+VUq8CPwceFpH/wQkrX6ga2+3WaDo4UdMm7LFA1YdITei4OtZ2j7Zn5wi7q0wqomZKPFqAUE6AEy4Zw8t/+JLZD37D2f97QJsTXTTtj3TPo33iiSfWUTvfqMWkNEallJqNk+SU/Nj/Jv28GDg0lTZoNO2JmGUTSofQkuTRmjY5bfVo3cHzm3dV01U5axn1DERoK32GdeGo80bw3pNL2PzdLnoN6ez5NTTZSSrn0eoWjBpNBhE17WZN32kz4pT3KKWImnbb92g7x4W2CjsROk6Nt7nPfj0IBA1W6cQoTTtBC61GkyHYtsK0VVqENu7RVpuOKHqxRwtxofW+vCeZcG6Qot55lG5oVhc9jcZ3tNBqNBlC1PUE0xI6dvc2q2PONdvq0XbKdWp/d1eZbCtxso9T5dECFPXMY8dmLbSa9oEeKqDRZAhxoU1L1rFb3lNtOtUKbfVoc90xe1Uxi3f+7nRZTaXQ5nXJoWK3HgivaR9oj1ajyRBiZho9Wpw0/yqPPNr4PNvKaE2ZYXWF2aY1GyOvMIxZbRGrblVZoybL0fNoNZosJe7RpmWP1h38Hvdo2+pFhwIGoYBQGasRvr17vHlHbqETqq7crXsfa1qOnker0WQpMdNJIkqfR6sSHm3cI20LkVCAypjF6HHdWb1gG8Mm9Wzzmg1eq8BpU1lZHqNT99yUXUfTML/96Lf9V5St8HQe7T5F+1TcdOhNGTOPFmDRokU555xzzmDLsuTYY4/d+dBDD/WsqKiY35LXpT1ajSZDiFqON5iW8h7w1KMFZ5+2KmZhW4oeAwuRFNTRxsnJdXyEaGXqwtOazCPd82gBrrjiiv4/+clPtixcuPDbPn36tCoxQHu0Gk2GEHU92nAKk4jiiNTeo/XCo80NB6iMWlgmBFJ8s5CT5/zpSuU+sKZxmvI8U4Ef82jnz59f8Pbbb68AuOSSS0pvuOGGfi1dQ3u0Gk2GkM49WvB2jxYgEnRCx7ZlYwRTe7MQTvJoN67cyepv9Pi8bMCPebReoIVWo8kQYnGhDbTdu2wK52+PSjSs8ELcI+EAlTEby1Sp92hdoX3vqSW8dMcXvHHfgsTAeU3HxY95tBMmTCh//PHHiwAee+yxrk2dXx9aaDWaDCGaKO9JQ+gYZ4/WtOMJWG2/Zm7IoCrqerQpTugK5ez993LHJj3Rp6Pjxzzae+65Z90999zTc+zYsSM3btwYKigoaHH4We/RajQZQtRD77IpnIYVYLpedNDwJhmqdE8UKyYEUhw6FkOIFISoKq/JTdlVWknXPvmNPEvTEUj3PNpBgwbFvvrqqyWGYfDQQw8VjR07tsUtybTQajQZQlpbMCIoFKbblzjohUcbDlBZZmFZgZR7tABDxndn8UcbE7/v2laV8mtq/GfGjBkDly9fnltdXS3Tp08vTfU82o8++ijvZz/72QClFJ06dbIef/zx1S1dQwutRpMhxD3a9LVghJjtnbjH62ht0ybQxpaOzaFuXkx1hW7JmA2kex7tbbfdtmnp0qWL27K2FlqNJkOIpdWjdcp7Eh6tBzWv8Tpay1IE0rDPXBddU6tpC3oerUaTBaR3j9Yp74mLe9ADcc8NOXW0tmljpOE11M0xjlbpvseazEQLrUaTIcTSWkfrtGC0bO882njo2DLtlE7uSVCnnGfJxxup2qPDx5rMQwutRpMhVKdzeo8bO46X93iVDGUrZ4B9IA2voS62pXj05x9iuzcsGk2moIVWo2kBVnk5W++9j1Wnnc7WP99D5aJFlPzsKpYdOoV1l1+BVd76YeQxd780bclQJO0Le1DeE2/jqCyVFo827s8WdovUelxnH2syDS20Gk0zUbZNyU8uZ9u994JhsO3++1l9xg/YPWcOVmkp5e++y4ojjmDrn+9BWS3fL4ymdR6t057OtBSGgOFJ6NgA5WQDp3KgQJx4O74Dpg3mez+vGStavqM65dfWtG/0PFqNJkPZ8Y8XqJg3j143/o4hL7/E4Ff+SfFVVzHohRfY59136HffvUTGjmXb/fez8Te/RdktC2HGLJuAIQTSIFLx9q8x2/YkEQqcG4T4Sl4Id5O4Lq0AfYYVcdzFowAoL9MeraZxWjuPNhZrXQ6ALu/RaJqBUortjz9OZNw4uvzwhwBEhg8nMnx44pxQ374UHnMMW++9j2333kt48GC6z7q02deIWnZa2i9C7fKekEeiGApIjdCmI3RcJ+14yIRiAMq3a6FNBxt+9ev+1cuXezqPNmfYsIo+N/8+o+bRnnHGGYOKiorMb775Jm/cuHEVDz/8cLM6SiWjhVajaQYVn35K9Lvv6H3rLU1OAym+4nKqly9n2733UnjM0eQMHdqsa0RNm3Cakoji5T2m5Z1HGzSMtArtgFFdWf7ZZrr2LXCuHw4QKQixu0yHjjsqyfNoY7GYTJgwYdTEiRMroGYe7Ztvvlkwa9aswcuXL1+UPI8WWjdUAGDlypWRjz76aFkw2DrJ1EKr0TSAUorKL74gPHAgZc88i9G5M52mTm3Wc3v99jes+uQTSi6/gsITpxLq2YvCY44mWFzc4HOilp220p6ER2srT0p7wPVo4+HcNISORxzcm4Fju5FbEE48Vtg1Qvl2LbTpoCnPMxX4MY8W4Pvf/35Za0UW9B6tJouxq6qwKysbPF7+/vusmXEeyw87nN1z5lB05pkYkUiD5ycT7N6dvvf8GbuqitIHHmTTDTew6tTTqF6xosHnpNOjRUgkQ3lR2gN1PNp07NFCLZEFyOsUZu2i0lrDBjQdB7/m0RYUFLSpZkwLrSarMMvK2HrvfXx35lks3f8Alk85jMqvvtrrPLuqii133AlAp2nTKJoxg+4/vqxF18o/4ACGvf8eI75dzIC/PgaBACVXXYUdjdZ7fizNHu3rCzby76VbPJncA04tbvxPWloaVtTDxhWOg/PFW6t9ub4mtfgxj9YLtNBqsobSxx9nxVFHO+U5lkWXH5yBRCJsvPFGlFIopbB27ya2YQNLJ0wkumoVvW++mb533E6v3/waI691eR8iQv7BB9Pn9/9HdMVKSv/yUL3nRU07LaU9cZsAtu6u9iwBy8k6dtZKR+i4PgaO6QZATl7Il+trUosf82i9QO/RarKC0kcfZcsdd1Jw1FH0uObniQSlyPARbLrhBirnz6d87lxKH3qYQFFR4nl5kyd5ZkPBEUfQ6aQTKX3kEbpedCGBgoJax9Pt0cbxLhmqZo/WL4/28LOHs/zzLQTS9D5q0k+659G++OKLez2/pej/jZoOjVKKrffcy5Y/3EXhCSfQ7957amUBd552CkZhISWXX0HpXx4iMno0oV69EsdDfby90e16/vmo6mp2z3l7r2PVafVoa372KhkqmO462noI5zq+w39fangvXNO+mTFjxsARI0aMGjdu3Mhp06aVpXoerRdoj1bTodn58j/Zdt99FJ5wAn1uuRkJ1N6iMfLzKb7ySjbffDM5I0cy8PG/YuTns+3Bv1D59ddIGzIN6yMyfjyhAQPY+dprdDnj+7WOpdejrRFCr5Khkuto/QodJwu8bSvfBF+TOvyYR9vWtbXQajosyrYpffhhIqNH0/fuu5AGkn66nn+eU3rTowcScvb2ul/2o5TYJCJ0njaNbfffT2zTplrec9S0yQun5yNZ26Ntn3W0TWGZNkbYl9wXTTtEz6PVaFrB7n+9Q/S77+h68UUNimycUN++CZFNNZ1OnApKUT53bq3H01lHm4x3yVBJe7QZ4ElaMT3FR5MZaI9W06FQSrHzlVewd+2m7LnnCA8a1OwmE+kiPHQowZ492fPxxxSdeWbi8Zip0teCMcml9cqjrdVVo6rIAAAgAElEQVTr2IcxeXWxTC20msxAC62mXWPt3IldVU2op1Mqt+O559l0ww2J433/ePde+7J+IyLk7jeRqq8X1Hrc8WjTY2vtrGOvkqEkUd7jp0e774E9WfbpZu3RajIG/287NZpWUvH55yw/6mhWHHEEm2+5xWn8/9STRMaMod/999Pv/vsyzpuNkzt2HLENGzC3bUs8lt5exzU/e9eC0UgIuPi4RztwtFNLa5k2lG+F16+GmB40oPGPlH6qRWSqiCwVkRUicl0D55wpIotFZJGIPJ1KezQdBxWLsf7n1xAqLqbz6aez/Ym/se6SS4muWEnR2dMpPPooCo8+2m8zGyR3/DgAKhd8k3jM8WjTL1ApqaP10aMNhJzXY5k2/PtG+PxRWPyKb/ZoMo8OM49WRALAfcCJwCjgbBEZVeecYcD/Aw5VSo0GrkqVPZqOxZ5PPsHcvJke1/6C3jf/nk6nnMKejz4iPGQInU480W/zmiQyahQYBlULFyYei1n+eLRe7QvXqqP10aONN6uwYgqUGz42G+5prck+WjuPtrWkco/2AGCFUmoVgIg8C5wGLE4651LgPqVUGYBSaksK7dF0IHa9/gZGYSH5hx2GGAZ97ridHlf/j1Oi43HtayowcnMJ9e5NdO3axGNpbcFIKpKh/K+jhSShNW0IuEMHlr4Fky70zaaOyLt/+7b/9vXlns6j7dq3oOKY80dm1DzaESNGJBzE1atXR1588cVlJ598cnlLXlcq/yL1BZLfsBLgwDrn7AsgIh8BAeAGpdRbdRcSkVnALIABAwakxFhN+8GuqmL3O+9QeMIJGGHnD6mIeN7FKdWE+vcnVlLTBS6tDStS0RnKh+k99ZEIHcdsMNySrWVv+maPxjv8mEe7ZMmSxQBPP/105z/84Q+9jj322GaLdJxUCm19n7S6M46CwDDgSKAf8KGIjFFK7aj1JKUeAh4CmDx5csNzkjQdCmWa7H7nHaoWLSYyejQFRxyOkZtL+QdzsffsodPJJ/ltYpsI9etL+fsfAE4Xo5il0ujR1uBpZyifex1DjUdrxqya0LHGc5ryPFOBX/Nov/nmm5xf//rX/d57771lOTk5LdagVAptCdA/6fd+wIZ6zvlEKRUDvhORpTjC+1kK7dK0A8zSUkp+cjmVX3+deCwydiwDn/wbu954g0D37uQfWDdA0r4I9+uPtW0bdkUFsXCO81i6GlYk19F6JO4iQkD8nd4DUFCUgwiULCljcGi382BuUeNP0rQL/JhHu2vXLuPMM88c+sADD6wZNGhQqwYdp/JT/RkwTEQGi0gYmA68WuecfwJHAYhId5xQ8io0WU31qlWsPvMsqpYupc8dtzP8yy/oddONVH3zDWsvnkn5e+/RaerUdrEX2xih/v0AiK1fT9RtrpDjw/SekIeiGBL/62jzO+fQpWcee3ZGodoV2uoWbalpMhQ/5tFOnz590Lnnnrtt6tSprf5PlLJPtVLKBK4A5gDfAs8rpRaJyI0icqp72hygVEQWA+8Bv1BKlabKJk3mY0ejrPvxj7Grqhj45N/oPG0aRl4eRT/8Id1/8hMqv/wSgG4XXeivoR4Q7ucIbXRdSUJo/ZjeE/AoGQpqRNvvXseBkOEkQ1Xvch6wY2BGfbVJ03bSPY922bJl4bfeeqvoqaee6j5ixIhRI0aMGDV37twWJ4A1yyUQkV7AgOTzlVL/bep5SqnZwOw6j/1v0s8KuNr90mjY+oe7iK1ZS/9HHiF37Nhax7pfeQU5I0eQM2QIob59fbLQO4K9ewNgbt5EzHJCYn7Mo/Wy7WNQMkNogyEDM2pBYFfNg9FyCHZt+EmadkE659Huu+++Udu2v2irzU0KrYjcDMwAlgDxjWQFtO9MFE3GEdu4ke1PP02XH/6QgimH7nVcROh03HE+WJYagu6AeXP7dh88Wu/H5EGN0Pq5RwvJHu3umgej5ZCnhba9M2PGjIHLly/Pra6ulunTp5d2lHm0ZwD7KqV0DzNNStn++BNg23T7UWpG1GUaEgphdO6MVbqdqOUIrR8erVd1tJBJHm2Aisoo2LucRKjKMr1P20HoqPNov0P3RNakGGXb7Hz1VQqPO45wv/YfFm4uwaIizLIajzactuk9STZ46H3G1/J7TF7Co7V2Q7ehjtBGW1z+qNkb27ZtMQyjw5VZtmUerW3bAjRYS9Ycod0NzBeRd4Dq+INKKb2vqvGM6mXLsMrKKDjiCL9NSSuBrl2xtpcRS7tH6315D0AQAZTvoePEHm2wGgp6wpbFTuhY01YWbt26dVRxcfHOjii2rcG2bdm6dWtnYGFD5zRHaN9yvzSalFHxhZNvkH/A/j5bkl4CnTsT27wpETpO1x5tcuzY22Qo57vfoeNAyMCKWc5fuMJezoNaaNuMaZqXbNq06ZFNmzaNQUc649jAQtM0L2nohCaFVin1qIgEgX3ch1a4pTsajWdULVpMoGtXgu2sjWJbkWAQTItYInTsxx6td6LoNKxQvoeOg6GA0xkqF8ejBR069oBJkyZtAU5t8kRNLZqTdXwY8CSwHufz2UtEzlNKfZRq4zTZQ9XixURGjWpzF5d2RyCAsiyq0x06Tq6j9TJ0nElZxzE3shn3aJMzkDWaNNKc0PHdwElKqcUAIjISR3gnp9IwTfagLIvoqlXkH3KI36akHQkEwDQTHq0f03u87AwVwKn98/uGKRgyMGMKpUAK3MZA2qPV+ERzPtXhuMgCKKW+BcKpM0mTbcQ2bkJFo4QHD/LblLQjQcejje/Rpq0FY3LWsYfiHhBBZUBQIj7BxyYI+cWA6D1ajW805xP2pYj8RUSmuF8PAPNTbZgme4iuWQ1AeOBAfw3xg0AQZVu+tmD0MhkqANgZILRBV2hNFYacQggXOB6tUvDdXOe7RpMmmvOpvgxYCVwL/BKn6X92dBTQpIXomjUAhAcO8tcQH3BCx1bay3uSCXgaOpa9ZmH6QS2hDRdAToGzR/vNC/DENPjq7z5bqMkmmpN1XAXc7n5pNJ5jbt4CgQDB4u5+m5J+AoYTOk6zR5uMl52hDMkMjzYx/F2FIJzvfEX3wNZvnRN2rm/k2RqNtzQotCLyjFLqbBGZz94D21FK7ZdSyzRZg7l1K8Fu3RAP/+C3FyQQBNMkmuahAsmRU29Dx5ni0TrT0CxcjzZc4OzRRt22uOEWD2DRaFpNYx7tL9zvP0iHIZrsxdy6lWCPHk2f2AFJJEOluY42GS+ToQyVWR6tEzrOr9mjjblCG8jx0TpNttHgJ0wpFR8htAFYpZRa6f4+HFiTasM02YO5ZQvB4mK/zfCHQBBl277u0XpZ3mPQSMPXNJIIHRuFYARq9mjj/rapZ6Ro0kdz6mg/BA4Xkc7ABzgZx9OB81NpmKb9U7K7hBU7VgDQJacL44rHYcjeQmJu3Uru+PHpNi8jkIDhhI5Nm4AhniYmNUZy6NjLazp9ofwnkQwVLHQeiO/RGiHn91jGT1bTdCCaI7SGUqpCRC4G7lVK3SoiX6XaME37pdKs5Hcf/47Zq2ajkv7sDisaxmPHP0aXSJfEY8q2scrKCHQt8sNU/3E7Q8Us29O90pbgaegYMDNAagPBuEfr7sXG92iVO1L74/vgyOt8sk6TbTTnE2aIyP7AOcDr7mOB1Jmkae/c9PFNzF41m4vHXMzTJz3Ns6c8y02H3sSqHau47bPbUEnulF1eDkoR6NylkRU7LhIIgmVRHbN82Z8Fb5OhDJUZoeP4UANb3N468T1ay23TXr3LJ8s02UhzPNqrgd8BbyilForIEJxwskazF+t2reO1Va9x8ZiLuWrSVYnHR3cbzYbyDTzw9QMc2f9IThh0AgDWLucPXqBTJ1/s9RsJOvesZsxM6/5scqTBy/IeASzPVms9NULrJj3luB6tFXV+z+vmk2WabKTJT5hS6t9KqZOUUr8Xp4HpZqXUT9Jgm6Yd8o9l/yAgAc4Zcc5ex2aNm8XIriP5v0/+jy0VWwCwdu4EINClc1rtzBgMR2hjMdM3jzbouUfrf+jYcG8eEkIbzne+Vzn/37BM+PZ1uH0ImNX1rKDReEeTn2wR+ZuIdBKRPGAR8J2I6KHvmr2otqp5ecXLHD3gaHrm99zreNAIctvht1EeLeexhY8BYMeFNus92hihdHq0SVro5Zg8ydTQccjdq42HjK0ovHUdVJTCBp1yokktzflkj1VK7QJOB94G+gEXptIoTWazvWo7Ly9/mdU7V9d6fPaq2eyo3sGZw89s8LmDOw/myP5H8t7a94Ca0LHRKUs92oDbWCHNHm2yz+llN6pMCx1bcaENuN/jE3zsGCj3luCx49NsnSbbaM4ebdgd/H4a8IBSKioimXDTqvGQDeUbeH3V6+zXYz8m92p4AuKaXWs4/83z2V61nZARYkKPCWwo34Bpm2yp2MKYbmM4sNeBjV5rZLeRvLP2HfbE9mDtdPdoszR0LAHnI2hGY4SDIV9s8DJ0LAqsTAgdB+KhY/c9Dboh5ITQmrBLt2HUpIfmCO0jwFpgIfCBiAwA9LypDkRFrIKZc2ZSUl6CIQa3HnYrJw4+ca/zbGVz/X+vx7RNHjj2Af654p8s2raIoV2Gsju6m0GdB3HzlJubnEU6tMtQAFbuWEmfXdkdOsYVBNO009vnOFV1tG7oWCnl60zavULHdT1ajSaNNGeowN04w98BEJF1wNGpNEqTXu6Zfw8l5SXcc/Q9PLbwMa6dey1KKU4aclKt815e/jJfbP6CGw+5kSl9pzCl75RWXW9w58GA4x332l0OwSASibT5dbRHEh5tLEo4P98XGwIeCmJcaG0FPpUFA0lCGx+dHfdoyzf5ZJEmm2lsqMDZSqlnROSnDZzy5xTZpEkRlWYld31+F59s/IQpfadwzeRr+GbbN/z9279z9oizObL/kRzS5xAufftSbvzkRoYVDWNY0TB2Vu/k9VWv86cv/8T+vfbn9H1Ob5Md3SJOacX2qu3Ye/Zg5Of76v34STwZyvaxjtbwVGgVNmDaNgHDv3L7Go/WDR3r3sYaH2nMo4236snSJrQdi/JoOTPfnsni0sWM6TaGp759il75vfh4w8d0y+3GVfs5Na/hQJjfT/k9584+lx+89gO65HRhR/UObGWzX4/9uGXKLW0WxU7hTgQlSFlVmSu0WTxJJTkZyq86Wi9dT3eogGX7u08biO/REt+jDftojSbbaVBolVL3u99/mz5zNKlgW+U2Zs6Zydpda/nTUX/iqP5H8dP3fsqdn98JwJUTryQvVCN2/Qr78cK0F3h26bMsL1vOsKJhHDvgWEZ2G+mJPSJCl0gXyqodoQ34FDLNBMQVWts0ifgQaz1gUFcKI94lYYlbR2v6LLQ1oWPt0Wr8p8k9Wjf56QpgUPL5Sqnvp84sjVcs2raIaz64htKqUh447gEO6n0QANcf7CQ1RQIRLhx94V7PK84r5sqJV6bMrq6Rrm7ouAIjTwut0xkq/aHW0yb28XZBN3RsWf4KrRiCYGPH/2QFtdBq/KM5WcevAn8D/kVm1KJrmsmi0kVcNOciuuR04eHjH2Z8cc2EnO653Xng2Ad8s60oUuSGjk2MggLf7PAdNxnKMs20DhWIN6zwcn/WWdhJaPbbowUwsGqE1sf9Yo2mOUIbVUrdlXJLNJ5z52d3UhAq4OmTn6Z7bne/zalF15yuLCpfhLVHCPbcu4tUtpBIhjJNcnyYReu1tEuG7NECGGLWCK0VqzmQ183pCJVMtALCWZwroEkpzflk3yMivxGR/UVkXPwr5ZZp2sRXW77i882fc9GYizJOZAEKw4WUx8qd0HEW79GStEfrR9ax58nedk3Wsd/U8mjzkz4D/Q/a++THT06PUZqspDke7b7AJcCJ1ISOFXB4qozStJ1HvnmELjldOGPYGX6bUi95oTwqYhXYe4ysFtr4Hu3o9UsIBQ5N23Xj/qbXZVXKraP13aNVCkNMLOWGjLsMACPodITq3Hfv8zd8mV77NFlFc4T2TGCQUkqPuGgnLN2+lA9KPuDyCZfXyibOJPJCeVSZldh7yGqhDfXtB8DEDYtY0QFCx3GP1nehtWIYmNjJo7M79YUda6DrUP/s0mQlzflkLwAKU22IxjseXfgoecE8zh5xtt+mNEh+MJ+QBVgWRl5m3gykg5whg8mdMIHKQE56WzC6eN4oRIEtKgOEthpDLGyV5EvY7tD34n39sUmTtTTHo+0GLBGRT4GEV6vLezKTZWXLmLN6DheMuoDOOZnbqD8vlEfYzU8xcnP9NcZvIhFyrB3pbVjhph172ObYwVbYRgZkHVsxZ49WJb2ncaEt7O2PTZqspTlC+/vWLi4iU4E/AQHgEaXUrQ2c9wPgH8D+SqnPW3s9Ddz1xV3kh/KZOXam36Y0Sn4onxxXaCU3O/scJ8jNI9es7hDJUEplSug46mYdJ4WO45nHoTwYegysfLf2c24fAteuSp+NmqyhOUL7X6BKKaVEZCgwHGcubaOISAC4DzgOKAE+E5FXlVKL65xXCPwU+LSlxmtq88ySZ/ho/Uf8fNLPM9qbhdpCa0Sy26NVubnkmtG0erRxvKyjVW4mVEbU0VpR16NtQGjrG+VXt+RHo/GI5nyyPwRyRaQ38AHwY+CxZjzvAGCFUmqVUioKPIsz07YuNwG3A1XNM1mTzLbKbTy5+EnOev0sbv70Zg7vdzjnjjrXb7OaJDeYS9iN5BlZ7tGqSC4RqzqtQuvlaLw48SYYTh2tz+U9VowAZp3QcVxoG7mxW/hSzQvRaDyiOR6toZSqEJGLgXuVUreKyFfNeF5fYF3S7yVArYngIjIR6K+Uel1ErmloIRGZBcwCGDBgQDMu3f5ZX76e3/znN2yu2MzNU25mQo8JgDMT9vVVr7OlYguLSxfzn/X/odKsZFS3UfzqwF/xw31/SNBozj+rv9QKHWe5R2tHcsk1q9OaDBU0nGt5mQyl3LaLNmD63ILRCR1b2HbSe3r6/fDeLa5H2wAvXATKhrE/SL2NmqyhWUIrIvsD5+CKHdCcfmb1fYITnz4RMXDm3F7Y1EJKqYeAhwAmT57c4W83Y1aMq967ipLdJVjK4roPr+O1018jFAjxpy//xGMLnYBCj9weHNrnUK6YeEVimHp7IS+YR07MTcjJco/WysklYsUIpzFyHJ/Y46Vja1mOF5sRe7RmtRs6TnqBY85wvqBxr3XVezDiFAhl9/9LjXc0R2ivBn4HvKGUWigiQ3DCyU1RAvRP+r0fsCHp90JgDPC+e1fdC3hVRE7N9oSoV1e+ypLtS7j7yLsJB8Jc/u7lvLryVcYVj+OJRU/wvX2+x/878P8RCUTa7RzX/FB+InScrUPf49g5TsP7HDOatmvGQ8fiYSWt7Xqxlvg/vQcrhiEmpt3Q3Usj9s1/yvn62QIoGpgS8zTZRZNCq5T6N/DvpN9XAT9pxtqfAcNEZDCwHpiO4xXH19kJJPqiicj7wDXZLrKmbfLowkcZ1W0Uxww4BoDxxeO5+dObyQ3l0incif+Z9D/kBtt3uDUnmEPE1RUjN3vraAHMcIQgkGOmL02hJnTs3ZqWmeTR+r3PmUiGauAFDjgYVr1f8/u4s2DsmfD3pE5qW5doodV4QpPBKhHZR0TuF5HZIvJ2/Kup5ymlTJzxenOAb4HnlVKLRORGETm17aZ3TN5a/Rbrdq9j1rhZiAgiwi1TbqE4r5iYFePWw26lKFLkt5ltJhKI6GQoFzPHuWkKRdPXfC2YgmSohEeL/2Py4uU9lt3A6zz8F3DK3c7PnfvD9x+CYcfC9Gdqzvnnj1NvpyYraE7o+AXgUeApnM9Qs1FKzQZm13nsfxs498iWrN0RsZXNIwseYZ8u+3BU/6MSj/fv1J9/nvZPqswqukS6+Gihd4SMEJFEMlSWC23YCR2H0yi0AXeP1ssQrx3fo5VMKO9xG1Y0JLRGALoNc37ukuS1jjgJLnzDGTKgy300HtEcobWVUvek3BINc0vmsnLnSm457BYMqR1siAQjRIIdR5BEhDwrCESzvjNULOz8u4Zi6Qsdh1yP1rS8K8Op5dH6LrRuw4rWvLxBU+DQn8HH94FlJmYGazStpTl5jq+IyCwRKRaRTvGvlFuWhTyx6Al65ffihEEn+G1KWsiznOR1cZOBspWoK7TBdHq07h6tl56nZcaFVvk/Js+qbtyjbYqiwU7Lxt0bvbVLk5U051btEvf7b5MeU0B2FLSmiUWli/h88+dcM/kaQkbIb3PSQq5pYIYD7TZz2iuiIedGIxRNZzKU85576Xkmh47992hjbh1tYyc1YmMXt2Bi57qanzWaVtKcrGP9vywNPLHoCfJD+Xx/WPbMaogLbbYTDYUBCFSnUWgDqQ0d+79HG3XG5DXn5dV3oxfft92xFgYe4qlpmuyjWZsPIjICGAUkNgmVUk+nyqhsY2P5Rt5e/TbnjDyHwnD2TCSMmAaxdHZpyFCq3b33QHVl2q4Z92i9DR1nUMMKK0pALOwWpW8m0WUAiAGlKzw1S5OdNCm0IvIb4HhgBE6pzgnAfwAttB7xt8V/A2DGyBk+W5JeckyIhbTQVrqhY6MqfR5tfI/W09BxfI82Y7KOTezG7Ag4kQRy6kk5CeZA0SDYtiwl5mmyi+b8lTsLOArYqJQ6DxhPMz1hTeNErSgvL3+Zp5c8zalDT6VPQR+/TUorkZgQDWX3/ixAVIJYCEZV+jzaUArKe+ItGC1U4mffiPc6bqyet/+BcNyNcNq99R/vPhy2LU+NfZqsojmCWamUskTEdEfabQKGpNiuDs9/1/+XX8z9BbuiuxhfPJ5fHvBLv01KO+GYIpodeV+NYipFdTCMpHGPNpDC8h6nM5Rny7YO0+0M1ZghIk4ZT0MU7+vMrI1WQDi7u5dp2kZzPNr5ItIFZzTe58A84MuUWtXBKY+W88sPf0mPvB7cd8x9PDHVSYTKNhyh1R5tzFJUBcKQRo82GEhFeY/r0WbEmDwnGUopUK19jUOPASsKC57z1jZN1tGoRytO3cUNSqkdwH0iMgfopJTSQtsGXlr+Ejuqd/DAsQ8wpvsYv83xjXBMUZXnt+vjPzHLpjoQQtK4R5vnZnsHPCytyrisYzeh3bZU6+bvDj4cug6B5W/D5Iu8tU+TVTQqtEopJSKvA5Pc33UKXhuxbIunlzzNfj32y2qRBQhFbapCWmhNywkdqzQK7YWHDGJHRYxLDvNuF8hOHpPnd+zYiuHme2FZNoHWJN2JQK9xsOJdqC6HnAJvbdRkDc353zdPRPZLuSVZwr/X/Zv15es5b9R5fpviO8GYTWXQ5xBjBhCzbaqDYezK9IWOI6EA1504glwP65gTnaEyYkxeNCG0je7TNsXIaRDdDeuzeqiYpo00KLQiEvd2p+CI7VIR+VJE5ouIDh23kicXP0nfgr61hgZkK8GoqYUWiJmKWDAnrUKbChKCZkhG1NF6IrT9Jjvfnzm77TZpspbGQsfzgP2A09NkS4dn4baFzN8yn2v3v5aAoTsiBaotKgM6Gcq0baLBMKrdC61z0ySGEMuEZCj3/1abhLZTX+d7rMIDozTZSmOhYwFQSq2s7ytN9nUoHlv4GPmhfL63z/f8NsV3VCyGYdlUBhVWq9v3dAxilk0slN7QcSqIC5oEJAP2aJM92jaIfiAEB/3E+bl8a9vt0mQljXm0xSJydUMHlVJ3pcCeDsmaXWt4Z807/GvNv7hs/GUUhHVShV3tTKqpDkG1VU2ekb11ijFLEQvlYJe3b6GNl/eIIRmxRxsIOh232uTRAgw7Dj65HzYvhAK95aNpOY0JbQAowPVsNa3jvbXvcfUHV2PaJvv12I+ZY2b6bVJGEA+TRoOu0IayV2hNyyYWzukAoWMFAsGgkQFj8mLehI4Beo51+h5//QwM1UKraTmNCe1GpdSNabOkA1JlVvF/n/wfQzsP5aZDb2J41+F7DXTPVmy3lCXu0WYzMUthhjpCMpRNIGAQNGzMTAgdx4W2raJfUOw0r9i4wAPDNNlIk3u0mtYzt2QuWyq3cPWkqxnZbaQW2STsCkdUtNA6e7RmKAcVjaKs9rtfbZkKIyCEAkYGhI5jGG73qzZ7tAA9R0Hpcijf0va1NFlHY3/5j0mbFR2UOavn0DXSlQN6H+C3KRmHqqodOs5mYpaNFXb3EyvT17TCa2zTxggKAUM87aHcKsxq70LHAKO/B7YJK95p+1qarKNBoVVKbU+nIR2NilgFc0vmctzA4wgaethRXeKCUh2CajO7hda0VUJoVWX7LSOxbOWEjgNCzHePNprk0Xog+j3HQiAHNi9q+1qarEPHMlPEByUfUGVVMXXQVL9NyUjsqnjoWKiy2q8X5wUxy8bKcYa/t+d9Wtu0ndCxYfjv0VoxjKBTq+5JqVEgCMXDYcvitq+lyTq00KaIOavnUJxbzMQeE/02JSOJ9/WNhpy5vNlMzFLYHSB0bJkKI+h4tBnRGSrgCK0noWOAnqNhsxZaTcvRQpsCyqPlfFjyIccPOl53gGqAROg4SNZ7tKZlo1yPtj2Hjm1LEQgIQUOIZULWcdDDZCiAHqOgfBNU6F01TcvQQpsC3lv3HlE7ygmDTvDblIwlngxVHdYebcxS2B0hdGzZGAGDYCAT6mijidCxJ3u04GQeg96n1bQYLbQpYM7qOfTM68n44vF+m5Kx1PJozez2aGNJHm17Dx0Hghnk0YacJETvPNrRzvfNC71ZT5M1aKH1mJU7VjK3ZC7Thk7TdbONEE+G0nu0TtaxiuQCYLfr0LHj0YYCRgbs0cYIBD0W2sJe0KkfrPmvN+tpsgZdd+IRa3et5fmlz/PW6rfIC+Vx/qjz/TYpo1GVVRAKYRsq6/doo6YNndzynjQOf/ca23IaVgQDglmVAXW0XoeORWDIEbB0NtgW6PwLTTPRLpcHrNu1jhmzZ/D0kqcpDBdyz9H3UBQp8tusjMaurMTIdcKl2d6wwrRtJOHRtl+hNaMWwZDhf/cfMK8AAB7USURBVOhYKXeP1mOPFmDQYVBZBluXerempsOjPVoPuGf+PVRZVbx82ssM7DTQb3PaBXZVJUZOBEOqs36P1rQUkuN6tNXt970wYzaBkEHQ8DkZyrYAhREMub96KLT93S5vJZ/VJEdpNE2gPdo2UrK7hLdWv8U5I87RItsCVEUlkpdLTiAn6z3aqGUnhLY9e7RWzCYYDjihYz/3aN3/T54nQwF0HQK5XaFknndrajo8WmjbyIvLX0REmD5iut+mtCvsykqMvHwigUjWC61pKULhEBIOJ8qe2iNmzCYYcpKhfJ3e4ybXGSHHo7W87FIl4ni1az/xbk1Nh0cLbRuIWTFeWv4SR/Q7gl75vfw2p13h7NHmkhPM0aFj2yZoCJKb2649WjPm7NH6PlTAigE1QuupRwsw4GAoXaEn+WiajRbaNvDuunfZXrWdM4ef6bcp7Q67ogIjNzfrPVqlFDFLEQwYGDk52O14j9aK2gTCAUJ+h47NeOg4RUI78BDn+9qPvV1X02HRQtsG/rH0H/Qt6MshfQ7x25R2h6qswMjLIyeQk9XlPXFBCgcEyY04ZU/tEKVUInTsJENlQug4DOJheU+c3hMgmAtrtNBqmocW2layZPsS5m2axw/2/YFuTNEK7IpKjDwndJzNY/JirggEAwZGJBe7ndbRWqb7OsLumLwMCB0TCGMExHuPNhiGfpNhrW5coWkeKVUIEZkqIktFZIWIXFfP8atFZLGILBCRd0WkXaTtVplV3DrvVjqFO3HW8LP8NqddYldWIjp0nKg3DRqCEYmg2mmvYzPqCm0oQNAQn5Oh3P9PgTCBgOG90AIMPBQ2fQNVO71fW9PhSJnQikgAuA84ERgFnC0idQvP5gOTlVLjgBeA21Nlj1d8vfVrznr9LL7c/CXXHXAdheFCv01ql9gVFRh5+Tp07Hp+4aCBRCLY1e3zpsOKOa8jEHKGCvjagjHJow0EjYS37SmDpoCydfaxplmk0qM9AFihlFqllIoCzwKnJZ+glHpPKRVv7voJ0C+F9rSZLRVbuPTtS6k0K7n/2PuZNnSa3ya1S5Rto6qqnGSoYCTLQ8dxj9Zo3x5tzAKc0HHIEGJ+NqyI984OhgmEjMRNgKf0mwyBMKz+0Pu1NR2OVHaG6gusS/q9BDiwkfNnAm/Wd0BEZgGzAAYMGOCVfS3m/q/uJ2bHePT4R+nfqb9vdrR34mJiuA0rstmjrdmjdct72ukeba3QccBwuiDaioAhPhiTFDoOkhqPNpQL/faH1f/xfm1NhyOVHm19n7B640kiMgOYDNxR33Gl1ENKqclKqcnFxcUemth8Vu5YycsrXmb68OlaZNuIXeEEMeJZx9m9R+uGjgOOR2u304YVZiwutEZCXH1LiEoOHYcCqRFacMLHG7+Gql2pWV/TYUil0JYAyYrUD9hQ9yQRORb4NXCqUipj/+L+8Ys/khfMY9a4WX6b0u6JDzcXHTpOlMEEA4JEclBV7fO9sNzQcSBkEAo4QuvbPm08dBwIEwhKakLH4DSuULZux6hpklQK7WfAMBEZLCJhYDrwavIJIjIR+AuOyGZsm5XPNn3G+yXvM3PsTD2VxwPiQmvk6jraROjYaN/lPTWhY6eOFvAv8zg56zhVyVDghI4loBOiNE2SMqFVSpnAFcAc4FvgeaXUIhG5UUROdU+7AygA/iEiX4nIqw0s5xtKKe76/C565vVkxsgZfpvTIUgOHUcCEWJ2DMu2fLbKH+LJUOFgvGFFJUr5PDS9FSRCx25nKMC/hKhE6DhEMGQkbPOcnALoPU4LraZJUjomTyk1G5hd57H/Tfr52FRe3wteX/U6C0sXctOhNxEJRvw2p0NQI7S5hANhwJlJm2fk+WmWL5h1PFqUQkWjiWk+7QUzKXQc8N2jjWcd5xAIGsQqzNRda8DB8PljYEadRhYaTT3olkaNsLxsOTd9chPjiscxbYgu5fGKRNaxu0cLEI3/ccwy4h5tKGBgRNyZtO0wfFwrdOx6tL7NpE3KOjZSGToGGHAQmFWwYX7qrqFp92ihbYCd1Tv52Xs/Iz+Uz91H3k3ACPhtUochORkqJ+CIS7bu01abjicYDgoSyQVol/u0Vj2hY/882tqh45QlQ4HTIQpgjS7z0TSMFtoG+P2nv2fjno3cfeTd9Mjr4bc5HQp7Tzx0nJ8Q2mwt8dlT7QhtQU4II9fx7ttj04p6k6F826ONZx3npDYZCiC/O/QYpetpNY2ihbYeFpcu5s3v3uTiMRczoccEv83pcNhJDSvioeNsnUm7p9rZP8zPCSAR571ojx5tYo82bBA04qHjDCnvSaXQglNPu/bTGk9ao6mDFtp6eGXFK+QEcrhw9IV+m9IhsStdjzYpdJytHu1uV2gLc0IYuW7ouKL9ebSxKstJhAo4vY7B72QoASPgNKxIZegYnPBxbA9s+Cq119G0W7TQ1sFWNu+sfYdD+xyqBwakCFVZiYTDSDBIJOB4cdkqtMkerZGfD4C9Z4+fJrWKaJVJOOLkMcSTofzrDBWFYA6IpMejje/T6r7HmgbQQluHJduXsKViC0cPONpvUzos9p6KhPeWE3STobI0dFxebRJxJ94Y+QVAexVai1COI7ShxB6tTx6tGXUa/oO7R5tiOwqKoXiE3qfVNIgW2jp8stEpPj+kzyE+W9JxsSsrkTynZjbbPdryapOCnBBAkkdb7qdJrSJWbRGKOGX5iYYVqfYkG8KKQsB5TwMhA2Ur7FR71wMPhXXzwEphza6m3aKFtg4f///27jxKrrpK4Pj31tbd1Us66WxNtg5JyAIRAgEmyDBsHtncA4MHOSCOmdED6ogoODo6jriMjorKOG4orgiIiAMSkUUQEUgIS0KCSci+dtJ7V3fXdueP96q7uro6SSf96lV13c85HKre+1W9W1VJ7ru/3+/93u5nmFs/l0lRf25eUA7SPT0E3ERb7pf3dPUmqXErwUC1852UYkWbyOo6roo4/4/FfVrtKxUH989VMOT8E+d5VTvrLIh3wr5XvD2OKUmWaLOkNc0rB17htCmn+R3KmJbu7u5PtJlZx+V6Y4HuviTVFU4lGCzpMdoUEbeijUac/8cSfiXaxEBF259ova5o3R6wrU97exxTkizRZtnVuYvuRDcLJyz0O5QxLd3ZSbDWGY8s94q2J5Ei6laAEokgkQiprtLrOo73Jgm7FW3m88T6fOpGTfUNjNGG3UTr9czjuuNg4nzYuNLb45iSZIk2y/qW9QAsaFjgcyRjW6qrk0CNM6O73C/v6U2kqAwPrDoWqK4uyYq2L5akIupUkdWZita3ruPEoMlQUICKFmDBpU5F29Pq/bFMSbFEm2VDywaCEmRu/Vy/QxnT0p1dBNyKtty7jnsS6cGJtqamf+WsUqGqxGNJKqqcBJsZo+3xret4YIH/YNi9N26hEq2mYOMj3h/LlBRLtFm2dmxlRu2M/irLeCPd2UnQrWgDEiAcCNOTKr1FGkZDbyJFVW5FW2Jdx4m+FOm0UhF1Em0k5KwO1e1X13FyoOs45H63mSUiPXXcqVAzFTY86P2xTEmxRJtlW8c2ZtbN9DuMMU1TKdKxGIHagcVAqkJV9CTKM9H2xPMk2hLrOo73OAk1k2jBGacthq7jkDtG69k9abMFAjD/Itj0x4E7CBmDJdp+qsqOzh3MrLVE66VMtZaZDAVQE64hliyt7tLR0pNI9Xe1gnOJT6kl2r5YJtGG+7dFIyFicb8mQw0sWBGKZCraAiX9BZdBvAu2PFmY45mSYInWtT+2n55kD7PqZvkdypiW6nQSbaBmINFGw1G64qXVXTpaenImQwVrakqu67g/0VYVS0Wb1XUccSvaQsUy+xyI1MD63xXmeKYkWKJ1be/cDmCJ1mPprk6A/lnH4FS03cnSquJGQyqtxJPpIV3HqVhpfRd9MeeuNRXVWYm2wsdEm+zrnwzVX9EWousYnDWW570JXnsI0j59flN0LNG6tnVsAyzRei3V0QFAsG4g0VaHq4klyq/ruNedlVsVGfhrGIhWl9ys4z53jDaSXdGGfew6jneDe0OQgle0AIveDt3NsPmxwh3TFDVLtK7tHduJBCJMrZ7qdyhjWqqtDYDg+PH926LhKF2J0uouHQ2Zy1+qci7v0VgMTZVONdTX7STUyuwxWj8r2r4uiDirbIUjBZx1nDH/EohOhFU/KtwxTVGzROvKXNoTEPtKvJRqdRNtfX3/tppwDd2J0uouHQ2xPicR5V5HC5TUOO1ARTvwOXwbo1V1JiNVON9jZmWoRCFjCUVgyVXwt4ehY3fhjmuKlmUV1/aO7XZpTwH0V7RZibZcu47be5yxzXFVA5VgqGECAMmDLb7EdDT6YgnClUECwYF/TqrCIX+WYEz2OotGRJxE2z9G21fgpH/qNU4ca35W2OOaomSJFudmAjs6d9j4bAGkWluRaJRAZWX/tupwNd2JblR9un+pTzKJtj4a6d8WbGgAIHXwgC8xHY2ejjjR2sigbdUVQX9uKtDn9gS4iTYQECKVwf6qu2Aa5sDx58KanzpVtilrlmhxbiYQT8eZPW6236GMeanWVkJZ1Sw4iVbRsruWtq0nDuRUtBMnAqVV0Xa3x4mOG5xoo5EQ3X3Jwp88ZS4Tqxi4fKyiOtw/jlxQiy+Htu2w+4XCH9sUFUu0wMa2jQDMq5/ncyRjX6qtbVC3MUBdpA6Ajr4OP0LyzUBFm9117FS0yRKqaGMdcaJ1g5ctrasKkUgpvYW6rCYjk2jdyVAAldVhet1LkApqwaUQCMO63xT+2KaoWKIFNrVtAmBO/RyfIxn7ki0tg2YcAzRUOcmlpbd0qrjRkG+MNlhfDyKkDh70K6wRi3UMrWgznynzGQsm7k6qi2RVtNEQfd0+JNqq8TDnfFj7G0j5dKmTKQqWaIFNrZuYVjONaDjqdyhjXnLfPkJTpwzaNqHSmQB0sLd0kstoaI8liIQCg2YdSyhEcPx4kgdK47tIxlPEe5JE64ok0eaM0YJb0frRdQyw9L3QsRNe/Lk/xzdFwRItTtexdRt7TxMJks3NhKcMvla5P9H2lEZyGS17O3qZUjf0TlGhhgaSLaXxXcQ6nHHmokm0cWflsSFjtH50HQOccBFMPwMe/wIkev2Jwfiu7BNtIpVga/tW5o63e9B6LdncDKrDVrTl1nW8p72XxrqqIduDExtIlUhFm0m01eMGnzDUVzmJty0WL2xA/V3HWWO00RC93T5MzAIQgfM+CV174aVfFP74piiUfaLd1rGNpCbtZu8FkNi7D4Dw1MEVbTQcpSpUVXZdx3vae2isrxyyPdQwkWSJjNF2tjhVWnX94EQ72a3U93UUuIrr7zoeWOKzojqMppVEr08rVR1/Lkw/HR7/IvS0+hOD8VXZJ9rMRChLtN5L7NoJQLixcci+ydHJ7O3eW+iQfJNIpdnb3stx9UMr2nBjI4k9e9Bk8U+gadnTjQjUTxn8OSbVVBAJBtjZVuD7DGdmrlcMHqMF6OkqcHWdIQKXfBViB+DBG5375ZqyUvaJdkPLBkISsmtoC6Bv4yYIhYjMGrowyKy6Wf03digHm5u7SKSU+VNqh+yLNM2CZJLE7uJfvq91T4zaiVWEsiZ0gbNQRGN9JbtaC5xou5uhYpxzFx1XXYPTa9BxwMcx0uNOgfP+Ddb+Gn7y9oEublMWyj7Rrtm/hkUNi4gEI4dvbI5J38aNVMxuQiJDv+tZdbPY0bmDtBb4ukufrN3lVF6Ljqsbsi9zIhLfsqWgMR2N1r3dTGiszrtvWn0Vuwpd0Xbth5rJgzaNm+xU2+3NBY4l1zkfg3d8F7b9GX75boiV15yEclbWiTaeirP2wFqWTF7idyhjnqrSu349FfNOyLu/qa6JnmQP+7r3FTgyfzy1sZmG6ghzJ9UM2VexYCGEQsRWF/eKQqlEmrb9McZPzX9Z3LT6Kn8q2pxEWz2ugmA44H+iBTj5Snjb/8D2Z+D758G+V/2OyBRAWSfaF/e/SDwdZ8kUS7Rei2/dSnLvXqJnnJF3/4kNJwLw0oGXChmWL5KpNE+81sy58ycTCMiQ/cGaaqpOPpnuZ57xIbojt2dzG+mk0jhnXN79MyZE2d/ZR1chby7QuhXGTR+0SQJC3cQq2vcXyRKfS66Cax9yLvf5wYXw6gN+R2Q8VtaJ9sEtDxINRVnWuMzvUMa8rsceB6D6rPzf9QkTTqAqVMXqvasLGZYvntvSQntPgvMXTB62TfWyZfSuXUuypXi7F9f+aReRyiDT5o/Pu//0Jueyrb9sKtBykn2d0LELJg7tNZk0o4Y9m9tJp4tkgf8Zp8OKJ2DyQrj7avj5FbCruHswzNHzNNGKyEUi8pqIbBKRm/PsrxCRX7n7nxWRJi/jybY/tp+VW1dy4awLbUUoj6VjMVrvuouqJUuIzMx/K8JwIMyyxmU8su0REumxOysznVa+8ceNNFRHDplo6y56MwAtd9xRqNBGpHl7J5vXNLP4vOlEKkN52yxtGk9NRYjHNuwvTFCZRNV48pBdTW+YSG9Xgt0b2woTy5Goa4RrH4QzPwCvP+F0JX/rNHjqa9BaPhMDy4FniVZEgsDtwMXAIuDdIrIop9n7gFZVnQt8HfiyV/Fk64h3cNOfbiKZTvL+xe8vxCHLVmLvXnbd9HESu3Yx6YbrD9n2ivlXcLD3IN9e8+0CRVdYHb0Jbrr3ZZ7b2sInLl5AVSQ4bNuKuXMZ99a30HLnT2j5xS/QRHGcfKSSaVb9fiv3fWU1FdUhTr5gxrBtw8EAF580lfte2MVLOzxOcPFu+Mu3IFgBM4YOT8w6qYGq2jBP/epvbH3lAOlUkUy6C1fCxV+Cf10Hf3+js+3R/4Db3gDfWAx3XwNP3wZbn7aZyiVMvFotRUSWAZ9V1Te7z28BUNUvZrVZ6bZ5RkRCwF5gkh4iqKVLl+qqVatGHM+TO5/k4S0Pk9Qkz+55lo6+Dm49+1YuOf6SEb+XGd7uWz4J6RSaStO7YT3x17eAKpM/8XEarr32sK//3DOf456/3UNTXRMLGxYSkhDLT1jOqVNO9T74UfaVlRvY0+ZcUtLc1cfqba3E4ik+dME8Pvqm/JPCsqXa2th5w4eIPf88wQkTiC5dSqCqksk33dR/Oz0v9XTGefreTf3Pu9v7aN7RSV93kuOXTOLsy+dRO2HoghvZmjv7uPSbT9Hc1ccZTROY5l43fN3ZszlpWv6x3SN2/wedyU8HNkKrO0P7/E/BOTflbb7j1RZ+/91XSPSlqKoN0zinnnBFkEkzaw95wlBwO56HjSth+1+dKj3hJthgBN7+HVi8/KjeVkRWq+rSUYzUHKH8fT6jYxqwI+v5TuDM4dqoalJE2oEGYNCgjoisAFYAzBym6/Fw9sX28cJ+p2vplEmnsOINKzhx4olH9V5meD1r1vQvtBCeMZ3aCy+kfvlyItOnH+aVjpvPuJnG6kZW7VvFy80vA3D+zPM9i9dLa3d18PoBZ6WicDDApYsbueaspiNOMMH6emb+5E66Hn+Ctl//mt5XnRmq2tfnWczZUsk0ezYPVKKBYICZCydwwplTaVp8ZIl+Um0FKz9yDt9+fBNPvLaf3e3OzN93nXZkfx4OacezkE7B5EUw5zznTjkLLhu2+YxFE7j61mVseekAm19o5sBOZ13kUKTIpqrMON35L+PARtjxHLx6vzOma0qOlxXt5cCbVfWf3OdXA2eo6g1Zbda5bXa6zze7bYZdf+5oK1pjjClnVtH6x8tTuZ1Adn/MdCB3qZv+Nm7X8TigeKdZGmOMMSPkZaJ9HpgnIrNFJAJcCeReMPYAcI37eDnw2KHGZ40xxphS49kYrTvmej2wEggCd6jqOhH5HLBKVR8Afgj8VEQ24VSyV3oVjzHGGOMHLydDoaoPAQ/lbPv3rMe9wOVexmCMMcb4qcim2xljjDFjiyVaY4wxxkOWaI0xxhgPWaI1xhhjPOTZghVeEZFmwKsVt8cB7R69d6kYq99BKX6uYo65mGLzM5aJ5KxkV8Rmqeokv4MoRyWXaL0kIt9T1RV+x+GnsfodlOLnKuaYiyk2P2MRkVW22pI5HOs6Hux3fgdQBMbqd1CKn6uYYy6m2IopFmOGsIrWGGOOklW05khYRWuMMUfve34HYIqfVbTGGGOMh6yiNcYYYzxkidYYY4zxkCXaYyQiC0Xkf0XkXhH5gN/xmNFjv60xZjSMqUQrIkERWSMi/3cM73GHiOwXkbV59l0kIq+JyCYRuRlAVder6r8AVwA2+9ADIlLvJrsNIrJeRJYd5fvYb2s8JSLVIrJaRC7zOxZTPMZUogU+DKzPt0NEJotIbc62uXma/hi4KM/rg8DtwMXAIuDdIrLI3fdW4M/Ao8cSvBnWbcDDqroAOJmc39h+W+OV4U7O8p2YuT4B3F3YKE2xGzOJVkSmA5cCPximyT8AvxWRSrf9+4Fv5jZS1SdxbkKf6wxgk6q+rqpx4C7gbe5rHlDVs4CrjvmDmEFEpA44B/ghgKrGVbUtp5n9tsYrPybn5Gy4EzMRuRB4FdhX6CBNcfP0xu8F9g3g40Btvp2qeo+IzAbuEpF7gOuAN43g/acBO7Ke7wTOFJFzgXcCFeTc5N6MiuOBZuBHInIysBr4sKp2ZxrYb2u8oqpPikhTzub+EzMAEcmcmNUA1TjJt0dEHlLVdAHDNUVqTCRadzxkv6qudv9xzEtV/8v9S/EdYI6qdo3kMPnfUp8AnhjB+5iRCQGnAjeo6rMichtwM/Dp7Eb225oCyntipqrXA4jItcABS7ImY6x0Hb8ReKuIbMXp9jtfRH6W20hE/h44CfgN8JkRHmMnMCPr+XRg91FFa0ZiJ7BTVZ91n9+Lk3gHsd/WFFDeE7P+B6o/VtWjnpBpxp4xkWhV9RZVna6qTcCVwGOq+p7sNiKyBPg+ThfPe4EJIvL5ERzmeWCeiMwWkYh7nAdG5QOYYanqXmCHiMx3N12AMw7Wz35bU2B2YmZGZEwk2iMUBS5X1c1ul8415LmvrYj8EngGmC8iO0XkfQCqmgSuB1bizHq9W1XXFSz68nYD8HMReRk4BfhCzn77bU0h2YmZGRFb69gYY4bhnpydi3OD933AZ1T1hyJyCc4EzCBwh6re6l+UpthZojXGGGM8VE5dx8YYY0zBWaI1xhhjPGSJ1hhjjPGQJVpjjDHGQ5ZojTHGGA9ZojXGGGM8ZInWGGOM8ZAlWmOMMcZDlmiNGSNE5Fsi8oKInO53LMaYAZZojRkDRKQamAz8M3CZz+EYY7JYojUlR0S+LiIfyXq+UkR+kPX8v0Xko6N8zJHc3/ZI3q9eRD6Y9bxJRNYe4WurRORPIhLMbFPVbqAR5/653xSRiIg8KSJj4p7TxpQyS7SmFP0FOAtARAI4C76fmLX/LOBpH+IaiXrgg4dtld91wH2qmspsEJEGnLsYdQIpVY0DjwL/eKyBGmOOjSVaU4qexk20OAl2LdApIuNFpAJYCKwRkftFZLWIrBORFZkXi8iXc6rJz4rIjSLyHhF5TkReFJHvZleMWW3ztnEr0vUi8n33eH8QkSp336dFZIOIPCIivxSRjwFfAua47/MV9+2D+V6fx1XAb3O2fQr4KrAOWORuu99ta4zxkSVaU3JUdTeQFJGZOAn3GeBZYBmwFHjZreiuU9XT3G0fcqs+gLsYXOldAaxyt71RVU8BUuQkKRFZeJg284DbVfVEoA14l4gsBd4FLAHe6cYCcDOwWVVPUdWbhnt97md37396vKpuzdrW5H4Pv8K5n26mul8L2MQoY3xm4zemVGWq2rOArwHT3MftOF3L4CTXd7iPZ+AksoOqukZEJovIccAkoBVYDJwGPC8iAFXA/pxjXnCYNltU9UX38WqgCadb+7eq2gMgIr87xGfK9/pcE3GScLbPA59TVRWR/kSrqikRiYtIrap2HuK4xhgPWaI1pSozTrsYp3LbAdwIdAB3iMi5wIXAMlWNicgTQGXW6+8FlgNTcSpcAe5U1VsOcczDtenLepzCScQygs+U7/W5esj6HCJyCk6lfLaI3O7ueyWrfQXQO4IYjDGjzLqOTal6GucylhZVTalqC84Eo2U4XcnjgFY3yS4A/i7n9XcBV+Ik23txJg4tF5HJACIyQURm5bzmSNrk+jPwFhGpFJEa4FJ3eydQO9IPraqtOGO5mWT7ZeAtqtqkqk3AybgVrdtV3qyqiZEexxgzeizRmlL1Ck436l9ztrWr6gHgYSAkIi8D/5nTDlVdh5PodqnqHlV9FWdC0R/c1zyCc7lM9msO2yaXqj4PPAC8BNyHMxbcrqoHgadFZG3WZKgj9QecCvZ8oFpVH8063j6gWkQmAOcBD43wvY0xo0xU1e8YjBnTRKRGVbtEJAo8CaxQ1ReO4f2WAB9V1asP0+4+4BZVfe1oj2WMOXY2RmuM974nIotwxk/vPJYkC+BO5npcRILZ19Jmc2cn329J1hj/WUVrjDHGeMjGaI0xxhgPWaI1xhhjPGSJ1hhjjPGQJVpjjDHGQ5ZojTHGGA9ZojXGGGM8ZInWGGOM8ZAlWmOMMcZD/w//yOboNmWYRQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for b in bands:\n",
    "    plt.plot(Table(data = parse_single_table(FILTERS_DIR + b + '.xml').array.data)['Wavelength']\n",
    "            ,Table(data = parse_single_table(FILTERS_DIR + b + '.xml').array.data)['Transmission']\n",
    "            , label=b)\n",
    "plt.xlabel('Wavelength ($\\AA$)')\n",
    "plt.ylabel('Transmission')\n",
    "plt.xscale('log')\n",
    "plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
    "plt.title('Passbands on {}'.format(FIELD))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### IV.a - Depth overview\n",
    "Then we plot the mean depths available across the area a given band is available"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gpc1_g: mean flux error: 8.010168924778235, 3sigma in AB mag (Aperture): 20.448092675884915\n",
      "gpc1_r: mean flux error: 8.685284950218165, 3sigma in AB mag (Aperture): 20.360236684401578\n",
      "gpc1_i: mean flux error: 8.62719361046711, 3sigma in AB mag (Aperture): 20.36752300186817\n",
      "gpc1_z: mean flux error: 33.421224501133366, 3sigma in AB mag (Aperture): 18.897140968910712\n",
      "gpc1_y: mean flux error: 39.71267202347856, 3sigma in AB mag (Aperture): 18.70987409051647\n",
      "gpc1_g: mean flux error: 7.715165791888239, 3sigma in AB mag (Total): 20.48883370546975\n",
      "gpc1_r: mean flux error: 10.3703401166686, 3sigma in AB mag (Total): 20.167714362688677\n",
      "gpc1_i: mean flux error: 6.177025432365751, 3sigma in AB mag (Total): 20.73024838959865\n",
      "gpc1_z: mean flux error: 13.343630063584328, 3sigma in AB mag (Total): 19.89401188030599\n",
      "gpc1_y: mean flux error: 48.299453845795625, 3sigma in AB mag (Total): 18.49734131339755\n"
     ]
    }
   ],
   "source": [
    "average_depths = []\n",
    "for b in ap_bands:\n",
    "    \n",
    "    mean_err = np.nanmean(depths['ferr_ap_{}_mean'.format(b)])\n",
    "    print(\"{}: mean flux error: {}, 3sigma in AB mag (Aperture): {}\".format(b, mean_err, flux_to_mag(3.0*mean_err*1.e-6)[0]))\n",
    "    average_depths += [('ap_' + b, flux_to_mag(1.0*mean_err*1.e-6)[0], \n",
    "                        flux_to_mag(3.0*mean_err*1.e-6)[0], \n",
    "                        flux_to_mag(5.0*mean_err*1.e-6)[0])]\n",
    "    \n",
    "for b in tot_bands:\n",
    "    \n",
    "    mean_err = np.nanmean(depths['ferr_{}_mean'.format(b)])\n",
    "    print(\"{}: mean flux error: {}, 3sigma in AB mag (Total): {}\".format(b, mean_err, flux_to_mag(3.0*mean_err*1.e-6)[0]))\n",
    "    average_depths += [(b, flux_to_mag(1.0*mean_err*1.e-6)[0], \n",
    "                        flux_to_mag(3.0*mean_err*1.e-6)[0], \n",
    "                        flux_to_mag(5.0*mean_err*1.e-6)[0])]\n",
    "    \n",
    "average_depths = np.array(average_depths,  dtype=[('band', \"<U16\"), ('1s', float), ('3s', float), ('5s', float)])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def FWHM(X,Y):\n",
    "    \n",
    "    half_max = max(Y) / 2.\n",
    "    #find when function crosses line half_max (when sign of diff flips)\n",
    "    #take the 'derivative' of signum(half_max - Y[])\n",
    "    d = half_max - Y\n",
    "    #plot(X,d) #if you are interested\n",
    "    #find the left and right most indexes\n",
    "    low_end = X[np.where(d < 0)[0][0]]\n",
    "    high_end = X[np.where(d < 0)[0][-1]]\n",
    "    return low_end, high_end, (high_end - low_end)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "data = []\n",
    "\n",
    "for b in ap_bands:\n",
    "    data += [['ap_' + b, Table(data = parse_single_table(FILTERS_DIR + b +'.xml').array.data)]]\n",
    "    \n",
    "for b in tot_bands:\n",
    "    data += [[b, Table(data = parse_single_table(FILTERS_DIR + b +'.xml').array.data)]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ap_gpc1_g (4260.0, 5500.0, 1240.0)\n",
      "ap_gpc1_r (5500.0, 6900.0, 1400.0)\n",
      "ap_gpc1_i (6910.0, 8190.0, 1280.0)\n",
      "ap_gpc1_z (8190.0, 9210.0, 1020.0)\n",
      "ap_gpc1_y (9200.0, 9820.0, 620.0)\n",
      "gpc1_g (4260.0, 5500.0, 1240.0)\n",
      "gpc1_r (5500.0, 6900.0, 1400.0)\n",
      "gpc1_i (6910.0, 8190.0, 1280.0)\n",
      "gpc1_z (8190.0, 9210.0, 1020.0)\n",
      "gpc1_y (9200.0, 9820.0, 620.0)\n"
     ]
    }
   ],
   "source": [
    "wav_range = []\n",
    "for dat in data:\n",
    "    print(dat[0], FWHM(np.array(dat[1]['Wavelength']), np.array(dat[1]['Transmission'])))\n",
    "    wav_range += [[dat[0], FWHM(np.array(dat[1]['Wavelength']), np.array(dat[1]['Transmission']))]]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Depths on SPIRE-NEP')"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfIAAAEgCAYAAABPZ/H0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xt4FOXZP/DvN+EYEsIpQEAwqJyRQIm0iop9BRtElFd8ERULyuGHBwoqKqW1tlpbsWgtSisHAQ+IilAERRGpgOIxFAgJR9EoCOEMIQQ5JPfvj53oEjfJbrKb3Vm+n+vai92ZZ565Z0DvfWaenZtmBhEREXGnmHAHICIiIhWnRC4iIuJiSuQiIiIupkQuIiLiYkrkIiIiLqZELiIi4mJK5BLVSOaQ7BXuOEREQkWJXKqMk1SPkzxK8jDJj0mOIhmUf4ckZ5P8czD6CiWSE0h+TTKf5E6Sr3mtW0Hye2fdfpILSCY76344PpIpJM1pl++c2/El9lN8vvO9Xs+WEZeR3OD990HyzyRnl7LP4teNXvGddJYdJLmMZLugnjwR+Qklcqlq/cwsAcC5AB4H8CCA58MbUtUhOQTArQB6mVk8gDQAy0s0u9tZ1wZAPQB/L6PLek7bGwA8RLJ3ifX9zCze63V3OSE2AzConDb1SvT5mte6J5x4zgGwF8DscvoSkUpSIpewMLMjZrYIwI0AhpDsBAAka5KcRPJbkntIPkeytrPuCmcEO8EZreaQvMVZNxLALQAecEaEi71214VkJskjJF8jWcvZphHJt5yrAwdJflja1QGSl5D8wunjC5KXeK1bQfJRkqudqw3vkWxUyqFfBGCpmW13zkOumU0r5RwdBDAfQCc/zmcGgGwAXcprW44nAPyJZLXKdGJmBQBegR+xi0jlKJFLWJnZ5wB2ArjMWTQRnpFoFwAXAGgO4A9emzQF0MhZPgTANJJtnWQ4B86I0Mz6eW0zEEA6gFYAOgMY6iy/z9l3EoAmACYA+Mkzi0k2APA2gMkAGgJ4CsDbJBt6NbsZwG0AGgOoAWBcKYf8KYBfk7yfZBrJ2NLOjfNlYACAtaW18Wr7C3iS5pfltS3HAgB5+PEcVQjJeHi+WJUbu4hUjhK5RIJdABqQJIARAO4xs4NmdhTAX/DTS70PmdkJM1sJT4IdWE7/k81slzPCXYwfR62nACQDONfMTpnZh+a7+EBfANvM7CUzO21mcwFsBuD9ZWGWmW01s+MAXkcpI2MzexnAaAC/ArASwN6S97YBTCZ5GMB6ALsB3FvGse0neRzAJwD+CWBhifULnSsOxa8RZfQFeL7IPATgDyRrlrFP7z7be60b58T+JYB4VPILgYiUr1KXz0SCpDmAg/CMjOMArPHkdAAAAXiPWg+Z2TGvz9/Ac1+3LLle7wu82v8NwB8BvOfsb5qZPe5j+2bOfrx948Rd2j7iSwvGzOYAmEOyOoD+zvu1ZrbUafIbM5tR5hH9qBE8yXcsgJsAVAdw0mt9fzN7v+RGJLPhmacAAH3M7EOv+JaQ/BbAyNL2aWanS1k3ycx+72fsIhIEGpFLWJG8CJ6E+BGA/QCOA+hoZvWcV6IzeapYfZJ1vD63hGdED/i4LF4WMztqZveZ2XnwjK7vJXmlj6a78GPS897vd4Hsz8f+T5nZPACZqMS9ZDMrNLMnAXwP4E4/t+noNVntQx9Nfg/gd/B8sRKRCKZELmFBsi7JawC8CuBlM9tgZkUApgP4O8nGTrvmJH9VYvM/kaxB8jIA1wCY5yzfA+C8AGK4huQFziX9PACFzqukJQDakLyZZDXn51YdALzl/xH/sM+hJPuSTCAZQ7IPgI4APgu0Lx8eh2eyX63KdmRmKwBsgGcegohEMCVyqWqLSR4FsAOeEd9T8EwSK/YgPPdXPyWZB+B9AG291ucCOATPKHkOgFFmttlZ9zyADs5925L3in1p7fSfD+ces5PAzmBmB+D5wnAfgAMAHgBwjZnt9+uIz5QHz6S6bwEchmeW+B1m9lEF+irpbXjOjfd98MUlfvP97wD6+z2ABj6WHy7RZ1n38EUkxOh7bo9I5CF5BTyj93PCHYuISKTQiFxERMTFlMhFRERcTJfWRUREXEwjchERERdTIhcREXGxqHqyW6NGjSwlJSXcYYiIuMaaNWv2m1lSJftoXK1atRnwPNhIA8TgKgKQdfr06eHdunXb66tBVCXylJQUZGRkhDsMERHXIFny8cMBq1at2oymTZu2T0pKOhQTE6OJV0FUVFTEffv2dcjNzZ0B4FpfbfTNSUREKqtTUlJSnpJ48MXExFhSUtIRlPEYZyVyERGprBgl8dBxzm2p+VqJXERExMWUyEVERILonXfeie/QoUP7atWqdZs1a1b9UO9PiVxERCSIzjvvvJOzZs3K6dev34Gq2F9UzVoXEZHwuv+N9S225h4Nah37Nk0TCv52Q+qO8tr16tXr/N27d9c4ceJEzKhRo/aMGzduf1xcXNdbbrll3+rVqxMSExML58+f/1WzZs1O+9p+5cqVcSNGjEiJi4sr+vnPf57/n//8J3Hbtm3ZkydPbvjmm2/WO3nyZMyOHTtqDhgw4MCTTz65GwCeffbZhpMnT25CEu3btz++cOHCr9u2bXsSAGJiyh8rFxYWYsiQIS0//fTThBYtWpwoKirC0KFDD9x2222H/D0/SuQiLvfOO+8gNzc33GGERdOmTdGnT59whyERYs6cOTlNmjQpzM/PZ9euXTsMHjz40PHjx2N+9rOfFUyfPn3nuHHjksePH9/sxRdf/NbX9sOHD2/1z3/+M6d3797H7rzzzube6zIzM+ts2LAhOz4+vqhr164drrvuuiNxcXFFkyZNSv7kk082Jycnn96zZ09soDG/+OKL9Xfs2FFjy5Yt2d999121Tp06dRo6dGhAI3klchERCRp/Rs6hMnHixCZvv/12PQDIzc2tnp2dXSsmJgbDhw8/CAC33377geuvv/4CX9vu378/9tixYzG9e/c+BgBDhgw5uGzZsnrF6y+99NK8pk2bFgJA3759D61YsSI+NjYW/fr1O5ScnHwaAJo0aVIYaMwffvhh/PXXX38oNjYWLVu2PP2LX/ziaKB9KJGLuJxGpCLAW2+9lbBy5cqEjIyMzQkJCUXdu3dve/z48Z9c2ybpc/vyCoiV3I4kzAwkK/Wzu2AULtNkNxERcb3Dhw/HJiYmFiYkJBStXbu21vr16+sAQFFREYpnjs+ePbth9+7dfY54k5KSCuvUqVO0fPnyOgDw0ksvNfBe/9FHH9Xds2dPbH5+PpcsWVKvZ8+e+enp6XmLFi1qkJubGwsAFbm0ftlll+UvXLiwfmFhIXbs2FHts88+Swi0D43IRUTE9QYMGHBk2rRpSW3atOlw/vnnf5+amnoMAGrXrl2UnZ1du2PHjk0TEhIKFyxY8FVpfUydOjVn1KhR58bFxRX16NHjaEJCwg+XytPS0vJvvPHGVjk5ObUGDBhw4PLLLy8AgPvuu2/3ZZdd1i4mJsY6depUMH/+/JyVK1fGDRw48IK8vLzY5cuX13vssceaffnll9m+9jlkyJBD77//fkKbNm06tmrV6vvU1NRj9erVC+gSfVTVI09LSzM9a11ExH8k15hZWmX6WL9+fU5qaur+YMUUTHFxcV0LCgrW+tP2yJEjMYmJiUUAMGHChKa7d++uPmvWrB2TJ09umJGRUae0SXKVVbzf3Nzc2Isuuqj96tWrN7ds2fKMmfXr169vlJqamuJre43IRUREALz++uuJTz75ZHJhYSGbN29+4pVXXsmpiv327t27dV5eXuypU6d4//337y6ZxMujRC4iIlHL12j81ltvbfnFF1/Eey+744479owZM+bAiBEjfvL77d/85jcHAFTq4S6ff/557V//+tetvJfVqFGjKDMzc/Pnn3++pTJ9K5GLiMhZ5aWXXgrJJfKydO/e/fjmzZs3hqJvzVoXERFxMSVyERERF1MiFxERcTElchERERdTIhcREQmiqKlHTrIFyQ9IbiKZTXKMs7wByWUktzl/+jxIkkOcNttIDglVnCIiIsFU0Xrkp08H9PPxH4Ty52enAdxnZv8lmQBgDcllAIYCWG5mj5McD2A8gAe9NyTZAMDDANIAmLPtIjPzuz6riIiEwcK7WmDvxqDWI0fjDgXoPyUq65G/9dZbCY8++mhy48aNT23cuDFu+/btPh/lWpaQjcjNbLeZ/dd5fxTAJgDNAVwH4AWn2QsA+vvY/FcAlpnZQSd5LwOQHqpYRUTE/ebMmZOTnZ29ad26dRunTp3aJDc3N7a4HvnGjRs39ejR4+j48eOblbb98OHDW02ZMuWbdevWbY6NjT3j+eWZmZl15s2b91VWVlb2okWLGqxatSouIyOj1qRJk5JXrly5dcuWLRunTp1aod+nZ2Zm1vnb3/72XUWSOFBFD4QhmQKgK4DPADQxs92AJ9mTbOxjk+YAvL997XSW+ep7JICRANCyZcvgBS0iIoHzY+QcKm6sRw4AnTt3PtauXbuTFdkWqILJbiTjAcwHMNbM8vzdzMcyn9VdzGyamaWZWVpSUlJFwxQRERfzrke+ZcuWje3btz/uhnrkABAXF1dUme1DmshJVocnic8xswXO4j0kk531yQD2+th0J4AWXp/PAbArlLGKiIh7ubUeeTCE7NI6PV9fngewycye8lq1CMAQAI87f77pY/OlAP7iNaP9KgC/DVWsIiLibm6tRx4MIatHTvJSAB8C2ACg+LLBBHjuk78OoCWAbwH8n5kdJJkGYJSZDXe2v91pDwCPmdms8vapeuQiIoFRPfIfhaseuT/CUo/czD6C73vdAHClj/YZAIZ7fZ4JYGZoojvTB7OnYe83pX5JkxBpfO55+OXQkeEOQ0QEQPjqkVeWypiKiEjUckM98sr0CyiRA4BGhSIiZ5Foq0euRA7gT4uzsXGXv7+Mk2Dp0KwuHu7XMdxhiIi4moqmiIiIuJhG5IBGhSIi4loakYuIiLiYErmIiEgQBVKPPCcnp3p6evp5ldmfErmIiEgQBVKPPCUl5dS7775bqQeZ6B65iIgEzUOrH2rx5aEvg1qP/IL6FxQ82uPRqKxHvmXLlhrXXHNN623btlX4Ea4akYuISFRwaz3yytKIXMLjnfFA7oZwRxEVJvIQNteoDjSo1G02V2rXoB0e7P5guMMQL/6MnEPFrfXIK0uJXEREXM+7HnlCQkJR9+7d27qlHnllKZFLePR5PNwRRA2NR0XKr0c+cuTIQ/7WI7/yyiuPlVaPvE6dOkVLliypN2PGjJy4uLiiG2644YIJEybsadq0aeGePXtiwzEqVyIXERHXUz3yKKF65CIigVE98h+pHrmIiIiLqR65iIhIhFE9chERkSgTbfXI9UAYERERF1MiFxERcTElchERERdTIhcREXExJXIREZEgCqQeeTAokYuIiARRIPXIg0E/PxMRkaDZNeF3LU5s2xbUeuQ1W7cuaPaXx6KyHvnYsWObvfvuu/UA4ODBg9Uuv/zyvDfeeCMnkPOjRC4iIlFhzpw5OU2aNCnMz89n165dOwwePPhQcT3y6dOn7xw3blzy+PHjm5X2qNXhw4e3+uc//5nTu3fvY3feeWdz73WZmZl1NmzYkB0fH1/UtWvXDtddd92RuLi4okmTJiV/8sknm5OTk0/v2bMnNtCYn3766V0Adh04cCD24osvbjtmzJi9gfahRC4iIkHjz8g5VNxaj7yoqAg33HBDq7vuumvPZZddVhDo9krkEhaHF2/HyV3Hwh1G1KjRrA7q9Ts/3GGIhI2b65Hfd999zZKTk0+OGTOmQvfUNdlNRERcr7x65ADgbz1yACitHnl+fj6XLFlSr2fPnvnp6el5ixYtapCbmxsLABW5tD537tzEFStW1J05c2aFr2RoRC5hodGjiASTW+uRP/3000327t1bvUuXLu0BID09/bBz39xvqkcuInIWUz3yH6keuYiIiIupHrmIiEiEUT1yERGRKKN65CIiIhIxQpbISc4kuZdklteyVJKfkNxAcjHJuqVsm+O0WUdSs9dERERKEcoR+WwA6SWWzQAw3swuBPBvAPeXsf0vzaxLZWdTioiIRLOQJXIzWwXgYInFbQGsct4vAzAgVPsXERE5G1T1PfIsANc67/8PQItS2hmA90iuITmyrA5JjiSZQTJj3759QQxVREQkcNFej/x2AHeRXAMgAcDJUtr1MLOfAejjtL+8tA7NbJqZpZlZWlJSUvAjFhERCUBU1yM3s80ArgIAkm0A9C2l3S7nz70k/w2gO368JC8iIhFq+YubWhz8Lj+o9cgbNI8vuPLX7aOyHnn//v1b3XDDDYcGDx58GACuvfbaVjfeeOPBW2655Yi/56dKR+QkGzt/xgD4PYDnfLSpQzKh+D08iT+rZDsRERFvc+bMycnOzt60bt26jVOnTm2Sm5sbW1yPfOPGjZt69OhxdPz48c1K23748OGtpkyZ8s26des2x8bGnvH88szMzDrz5s37KisrK3vRokUNVq1aFZeRkVFr0qRJyStXrty6ZcuWjVOnTg349+kjRozYN3v27IYAcODAgdg1a9bEDxw40O8kDoRwRE5yLoArADQiuRPAwwDiSd7lNFkAYJbTthmAGWZ2NYAmAP7tlIyrBuAVM3s3VHGKiEjw+DNyDhU31iPv27dv/tixY8/97rvvqs2ZM6d+3759D1WvXj2gPkKWyM3splJW/cNH210ArnbefwUgNVRxiYhI9HFzPfKBAwcemDFjRoP58+c3mDlzZk6g2+vJbiIi4npurUcOAKNGjdo/derUJgCQlpb2faDb61nrIiLiem6tRw4ALVq0OH3++ed/369fv8MVOXYlchERcb3atWvbqlWrtvla949//GMXgF3l9dGtW7fjW7du3Qh46pEXfxkAgEaNGp32VY989OjRB0aPHn3Gz8x69uxZsGfPnkx/Yz969GhMTk5OzWHDhpV8iJpfdGldREQEnnrk7dq169C6deuOH3/8cfxjjz22O9T7XLhwYUKbNm06jhgxYm/Dhg0DniwHaEQuIiJRzA31yPv377+hMn0rkYuIyFlF9chFREQkYmhELiKutHXroziavyncYfgtIb492rR5KNxhSBTSiFxERMTFNCIXEVfS6FbEQyNyERERx+jRo5s3bdq0c1xcXNdwx+IvJXIRERFH//79D3/22WfumXwBXVoXEZEgWvqvp1vs3/FNUOuRN2pxbsGv7hhbblW1+++/P/mNN95okJycfLJhw4anu3btWvDuu+/W69SpU8HatWvr5Ofnx06bNu3rX/7ylwVHjhyJGTZsWMvMzMw4AJgwYcKuoUOHHr7yyiuPlbefYtnZ2TVvvvnmVoWFhezVq9eRadOmNfH1u/VQ04hcRERcb9WqVXGLFy+uv2HDho1vv/329szMzDrF6woKCmLWrl27efLkyd+MHDmyFQCMHz8+uW7duoVbt27duHXr1o19+/b1WUylLHfffXeLO++8c29WVtamZs2anQrm8QRCI3IREQkaf0bOobBixYr4Pn36HI6PjzcA1rt37x8KkNx8880HAaBPnz75+fn5Mfv3749dtWpV3VdfffWHAipJSUkBPx517dq18e+9996XADB8+PADf/zjH88JwqEETCNyERFxvbLqiZdRSzzUYVUJJXIREXG9K664In/p0qWJBQUFPHLkSMz7779fr3jd3Llz6wPA0qVL4xMSEgobNmxYeMUVV+Q99dRTjYvb7Nu3L+Ba4l26dMmfPXt2fQCYOXNmg/Lah4oSuYiIuF7Pnj0L0tPTj3To0KHj1VdffX7nzp2PJSYmFgJA/fr1C7t27dru7rvvPnfq1Kk5APDXv/519+HDh2Nbt27dsW3bth2WLFmSAACjRo06p0mTJp2///77mCZNmnS+9957m5W2z2eeeWbHM8880+TCCy9sv3v37urx8fEVql5WWSzrcoTbpKWlWUZGRrjDEBFxDZJrzCytMn2sX78+JzU1dX+wYqqoI0eOxCQmJhYdPXo05uKLL2773HPPfXPvvfe2mDRp0o7LL7+8INj7O3r0aEydOnWKYmJiMG3atPqvvfZag+XLl28P9n4AYP369Y1SU1NTfK3TZDcREYkKgwcPPnfbtm21T5w4wUGDBh249NJLg568va1evTpuzJgxLc0MdevWLZw9e3ZOKPdXGiVyERGJCosXL/665LLPP/98S2X7ffDBB5u++eabZ9wDv+666w5OnDgxd8uWLSEpTRoIvxI5yZoABgBI8d7GzB4JTVgiIiKRYeLEibkTJ07MDXccpfF3RP4mgCMA1gA4EbpwREREJBD+JvJzzCw9pJGIiIhIwPz9+dnHJC8MaSQiIiISsDJH5CQ3ADCn3W0kv4Ln0joBmJl1Dn2IIiIiUpryLq1fUyVRiIiIRIDRo0c3nzdvXsO8vLzYcFQyq4gyE7mZfQMAJF8ys1u915F8CcCtPjcUEZGz0sE3trY4lXssqGVMqzetU9DghjZVUoylf//+h8eNG7e3ffv2nQLZ7tSpU6hevXqowiqTv5PdOnp/IBkLoFvwwxEREamYqq5HPmDAgJT69euf3rBhQ1znzp0Lpk+fvjN0R1e68u6R/xbABAC1SebBc28cAE4CmBbi2ERExGWqauRcknc98lOnTrFLly4dunbtWgD8WI/8nXfeiR85cmSrbdu2ZXvXIwcqVjQFALZv315r9erVW6tVC9/z1cq7tP5XAH8l+Vcz+20VxSQiIhKQcNQjB4Drr7/+UDiTOOD/pfUJJK8HcCk8s9g/NLOFoQtLRETEf+GqRx4fH19U6U4qyd/fkU8BMArABgBZAEaRnBKyqERERAIQjnrkkcLfRN4TwK/MbJaZzQJwNYArQhaViIhIAMJRjzxS+FWPnOQCAPd4/RztXACPm9lNIY4vIKpHLiISGNUjd4ey6pH7OyJvCGATyRUkVwDYCCCJ5CKSi3xtQHImyb0ks7yWpZL8hOQGkotJ1i1l23SSW0h+SXK8nzGKiMhZbPDgwee2a9euQ+fOndv369fvUKjrkUcKfye7/aECfc8G8CyAF72WzQAwzsxWkrwdwP0AHvLeyPmN+hQAvQHsBPAFyUVmFvaaryIiErnCUY+8sn0Hg1+J3Em85wJobWbvk6wNoJqZHS1jm1UkU0osbgtglfN+GYClKJHIAXQH8KWZfQUAJF8FcB08VwFERESqVKTXI/fr0jrJEQDeADDVWXQOgIr8/CwLwLXO+/8D0MJHm+YAvB8osNNZJiIiIiX4e4/8LgA9AOQBgJltA9C4zC18ux3AXSTXAEiA5wlxJfn6YV+pM/JIjiSZQTJj3759FQhJRETEvfxN5CfM7IekS7IaykiupTGzzWZ2lZl1AzAXwHYfzXbizJH6OQB2ldHnNDNLM7O0pKSkQEMSERFxNX8T+UqSxc9c7w1gHoDFge6MZGPnzxgAvwfwnI9mXwBoTbIVyRoABgHwOTNeREQkmEaPHt28adOmnePi4rqW1/aJJ55IevbZZxtWRVxl8TeRjwewD54nu/0/AEvgScSlIjkXwCcA2pLcSXIYgJtIbgWwGZ5R9iynbTOSSwDAzE4DuBueiXCbALxuZtmBHpiIiEig+vfvf/izzz7b5E/bBx54YN/dd999INQxlcffWetFJBcCWGhmft2ILuNhMf/w0XYXPE+LK/68BJ4vCyIi4iILFy5ssXfv3qDWI2/cuHFB//79y62qVtVlTO+9995m8fHxhY888sieyhxfZZVXxpQAHoZnhExnUSGAZ8zskSqIT0REpFzhKmMaCcobkY+FZ7b6RWb2NQCQPA/Av0jeY2Z/D3WAIiLiHv6MnEMhXGVMI0F598h/DeCm4iQOAM6DWgY760RERMIuXGVMI0F5iby6mf3kQfjOffLqoQlJREQkMCpjWjpfD2zxZ52IiEiVURnT0lZ6Jrb5msFHALXMLKJG5SpjKiISGJUxdYeyypiWOdnNzFx7qUFERM4ugwcPPnfbtm21T5w4wUGDBh1QGVMREREXURlTERER+YmoKGMqIiIikUmJXERExMWUyEVERFxMiVxERMTFlMhFREQcgdQjjxRK5CIiIo5A6pFHCv38TEREgmbjpgdbHMvfGtR65HXi2xR0aD8x4uqRt2vXrkPx+5ycnFrz58/f2rdv3/yKHWXFKZGLiIjrhaMe+ebNmzcCwCuvvJL45JNPNu3Vq5ffXwKCSYlcRESCxp+RcyiEqx75hg0bav7ud78754MPPthas2bN0ouXhJDukYuIiOuFox55Xl5ezMCBA8//17/+9U1KSsqpSnVWCUrkIiLieuGoRz5o0KCUW265ZX96enqV3xf3pkQuIiKuV9X1yLdu3Vrj3Xffrf/yyy83ateuXYd27dp1WLVqVVAn+fmrzHrkbqN65CIigVE9cneocD1yERERt1A9chERERdTPXIRERH5CdUjFxERkZBRIhcREXExJXIREREXUyIXERFxMSVyERERh+qRi4iIuJjqkYuIyFlt7KZvW2w+9n1QH1Xark6tgqfbt4yoeuSHDh2K6dSpU8evvvoqq2bNmnbw4MGYCy+88IfPlT3mQCiRi4iI61V1PfL69esXXXzxxUdff/31xFtvvfXwzJkzG1x99dWHwlHKVIlcRESCxp+RcyiEox75yJEj902cOLHprbfeevjll19uNH369JxgHEugdI9cRERcLxz1yK+66qpjO3furPn222/HFxYW8qKLLvq+Uh1WUMgSOcmZJPeSzPJa1oXkpyTXkcwg2b2UbQudNutILgpVjCIiEh3CUY8cAAYNGnTgtttuO2/w4MFhq/4WyhH5bADpJZY9AeBPZtYFwB+cz74cN7MuzuvaEMYoIiJRoKrrkRcbNmzYgby8vGrDhg07GPKDLEXI7pGb2SqSKSUXA6jrvE8EsCtU+xcRkbPLww8/nPvUU0/tKq5H/sADD+x57bXXGt54442HpkyZ8p1328TExKIFCxbklOzjueee2/ncc8/t9Hefy5cvT0hPTz/UqFGjgO+xB0tVT3YbC2ApyUnwXA24pJR2tUhmADgN4HEzW1hVAYqIiDtVdT3yIUOGtPjggw8S33pYAfWAAAANNUlEQVTrrW2h3E95qjqR3wHgHjObT3IggOcB9PLRrqWZ7SJ5HoD/kNxgZtt9dUhyJICRANCyZctQxS0iIhGuquuRv/DCCzsAhGWWvjeWNdOv0p17Lq2/ZWadnM9HANQzM6NnuuARM6tbRhcgOdvp443y9peWlmYZGRmVjltE5GxBco2ZpVWmj/Xr1+ekpqaGbbLX2WD9+vWNUlNTU3ytq+qfn+0C0NN5/z8AfnI5gmR9kjWd940A9ACwscoiFBGRQBUVFRVV7rdcUirn3BaVtj6UPz+bC+ATAG1J7iQ5DMAIAE+SXA/gL3AuiZNMIznD2bQ9gAynzQfw3CNXIhcRiVxZ+/btS1QyD76ioiLu27cvEUBWaW1COWv9plJWdfPRNgPAcOf9xwAuDFVcIiISXKdPnx6em5s7Izc3txP0oLFgKwKQdfr06eGlNdAjWkVEpFK6deu2F4Ce+REm+uYkIiLiYkrkIiIiLqZELiIi4mJK5CIiIi6mRC4iIuJiSuQiIiIupkQuIiLiYvoduYi4Uu5f/oITmzaHOwy/1WzfDk0nTAh3GBKFNCIXERFxMY3IRcSVNLoV8dCIXERExMWUyEVERFxMiVxERMTFlMhFRERcTIlcRETExZTIRUREXEyJXERExMWUyEVERFxMiVxERMTFlMhFRERcTIlcRETExZTIRUREXEyJXERExMVU/UxEJIQ+fH0r9u/ID+k+GrWIx2UD24R0HxK5NCIXERFxMY3IRURCSCNlCTWNyEVERFxMiVxERMTFlMhFRERcTIlcRETExZTIRUREXEyJXERExMWUyEVERFxMiVxERMTFlMhFRERcLKSJnORMkntJZnkt60LyU5LrSGaQ7F7KtkNIbnNeQ0IZp4iIiFuFekQ+G0B6iWVPAPiTmXUB8Afn8xlINgDwMICfA+gO4GGS9UMbqoiIiPuENJGb2SoAB0suBlDXeZ8IYJePTX8FYJmZHTSzQwCW4adfCERERM564SiaMhbAUpKT4PkicYmPNs0B7PD6vNNZ9hMkRwIYCQAtW7YMbqQiIiIRLhyT3e4AcI+ZtQBwD4DnfbShj2XmqzMzm2ZmaWaWlpSUFMQwRUREIl84EvkQAAuc9/PguQde0k4ALbw+nwPfl+BFRETOauFI5LsA9HTe/w+AbT7aLAVwFcn6ziS3q5xlIiIi4iWk98hJzgVwBYBGJHfCMxN9BIB/kKwG4Hs497dJpgEYZWbDzewgyUcBfOF09YiZlZw0JyIS8R7athNZ+cdDuo9O8bXxaOtzQroPiVwhTeRmdlMpq7r5aJsBYLjX55kAZoYoNBERkagQjlnrIiJnDY2UJdT0iFYREREXUyIXERFxMSVyERERF1MiFxERcTElchERERdTIhcREXExJXIREREXUyIXERFxMZr5LCrmSiT3Afgm3HFEuEQAR8IdRISI9nPh9uNzU/yRHGsjAPvLWH+umal0pItFVSKX8pGcZmYjwx1HJIj2c+H243NT/JEcK8kMM0sLdxwSOrq0fvZZHO4AIki0nwu3H5+b4ndTrBJlNCIXEYliGpFHP43IRUSi27RwByChpRG5iIiIi2lELiIi4mJK5CIiIi6mRC7lItme5HMk3yB5R7jjCadoPxfRcHxuOgY3xSqRS4ncRUjmkNxAch3JjEr0M5PkXpJZPtalk9xC8kuS4wHAzDaZ2SgAAwFExOxXkvWc//ltJrmJ5MUV7CfizgXJts7fcfErj+TYCvYVtuMjeQ/JbJJZJOeSrBWpx0ByjBNndkXPdVXFWlkk65BcQ/KaqtifVAEz08slLwA5ABqVsb4xgIQSyy7w0e5yAD8DkFVieSyA7QDOA1ADwHoAHZx11wL4GMDN4T4PTjwvABjuvK8BoF40ngsnjlx4nr7lmuMD0BzA1wBqO59fBzA0Eo8BQCcAWQDiAFQD8D6A1pEYaynxzwSw18c+0wFsAfAlgPFeyx8B8CCAa0L971evqnlpRB5degJ4s3jkQ3IEgMklG5nZKgAHfWzfHcCXZvaVmZ0E8CqA65xtFpnZJQBuCVXw/iJZF57/YT4PAGZ20swOl2gWLefiSgDbzazko4fdcHzVANQmWQ2eJLkrQo+hPYBPzazAzE4DWAngfyM0Vl9mw5O0f0AyFsAUAH0AdABwE8kOJHsB2AhgTwX3JRGoWrgDkIAYgPdIGoCpZnbG70PNbB7JVgBeJTkPwO0AegfQf3MAO7w+7wTwc5JXALgeQE0ASyoRf7CcB2AfgFkkUwGsATDGzI4VN4iiczEIwNySCyP9+MzsO5KTAHwL4DiA98zsvQg9hiwAj5Fs6MR6NYAzbl1FUKw/YWarSKaUWPzDFwcAIFn8xSEeQB14kvtxkkvMrKgi+5XIoUTuLj3MbBfJxgCWkdzsjAB+YGZPOP/R/gvA+WaWH0D/9LHMzGwFgBUVDToEqsFz+XK0mX1G8h8AxgN4yLuR288FyRrwXHb9ra/1kXx8JOvDkzhaATgMYB7JwWb2cokdhv0YzGwTyYkAlgHIh+ey92kf7cIeawB8fnEws7sBgORQAPuVxKODLq27iJntcv7cC+Df8HzrPgPJy+C55/dvAA8HuIudAFp4fT4HP70cGgl2AthpZp85n9+AJ7GfIQrORR8A/zUzn5dBI/z4egH42sz2mdkpAAsAXFKyUaQcg5k9b2Y/M7PL4bk0vi1SY/WTzy8OP7wxm21mb4Vw/1KFlMhdwplpmlD8HsBV8FwS9G7TFcB0eEZCtwFoQPLPAezmCwCtSbZyRoODACwKRvzBZGa5AHaQbOssuhKe+34/iJJzcRN8XFYHXHF83wL4Bck4koTn72iTd4NIOgbnKhdItoTnUvfcEusjJlY/RcIXUakq4Z5tp5d/L3juC693XtkAfuejTQ8AF3p9rg5ghI92cwHsBnAKnv/gh3mtuxrAVnhm2P5kH5HyAtAFnvuYmQAWAqgfTecCnslhBwAklrI+4o8PwJ8AbIbnC+dLAGpG6jEA+BCeL4PrAVzptvMNIAVes9bhuf30FTy3NopnyXesyn/DelXdS89aFxFxMZJzAVwBoBE8s9EfNrPnSV4N4Gl4fvo208weC1+UEkpK5CIiIi6me+QiIiIupkQuIiLiYkrkIiIiLqZELiIi4mJK5CIiIi6mRC4iIuJiSuQiIiIupkQuIiLiYkrkIlGO5DMk/0vyonDHIiLBp0QuEsWcAjuNAfw/ANeEORwRCQElcokaJP9OcqzX56UkZ3h9fpLkvUHeZyA1qf3prx7JO70+p5DMKmsbr7a1Sa4kGVu8zMyOAUiGp+b1ZJI1SK4iWS2YcYtI+CiRSzT5GE7Na5Ix8BSR6Oi1/hIAq8MQVyDqAbiz3Fa+3Q5ggZkVFi8g2RCeSmpHARSa2UkAywHcWNlARSQyKJFLNFkNJ5HDk8CzABwlWZ9kTQDtAawluZDkGpLZJEcWb0xyYonR8B9J3kdyMMnPSa4jOdV7xOvV1mcbZ0S9ieR0Z3/vkaztrHuI5GaSy0jOJTkOwOMAznf6+ZvTfayv7X24BcCbJZb9HsAkeErfdnCWLXTaikgUUCKXqGFmuwCcJtkSnoT+CYDPAFwMIA1ApjMivd3MujnLfuOMWgHgVZw5Uh0IT83zGwH0MLMuAApRIgmSbF9Om9YApphZRwCHAQwgmQZgAICuAK53YgGA8QC2m1kXM7u/tO1LHjvJGgDOM7Mcr2Upznl4DcAm/Hh1IguAJr6JRAndJ5NoUzwqvwTAUwCaO++PwHPpHfAk7/913reAJ1EeMLO1JBuTbAYgCcAhABcC6AbgC5IAUBvA3hL7vLKcNl+b2Trn/RoAKfBc9n/TzI4DAMnFZRyTr+1LagRPkvf2ZwCPmJmR/CGRm1khyZMkE8zsaBn7FREXUCKXaFN8n/xCeEaeOwDcByAPwEySVwDoBeBiMysguQJALa/t3wBwA4Cm8IzQCeAFM/ttGfssr80Jr/eF8CR6BnBMvrYv6Ti8joNkF3hG+peSnOKs2+DVviaA7wOIQUQilC6tS7RZDc/PrA6aWaGZHYRnAtnF8FxqTwRwyEni7QD8osT2rwIYBE8yfwOeiWE3kGwMACQbkDy3xDb+tCnpIwD9SNYiGQ+gr7P8KICEQA/azA7Bcy+9OJlPBNDPzFLMLAVAKpwRuXMrYZ+ZnQp0PyISeZTIJdpsgOcy86cllh0xs/0A3gVQjWQmgEdLtIOZZcOTSL8zs91mthGeCWPvOdssg+fnXN7blNumJDP7AsAiAOsBLIDnXvwRMzsAYDXJLK/Jbv56D54R+P8AqGNmy732twdAHZINAPwSwJIA+xaRCEUzC3cMImclkvFmlk8yDsAqACPN7L+V6K8rgHvN7NZy2i0A8Fsz21LRfYlI5NA9cpHwmUayAzz3r1+oTBIHAGey3gckY71/S+7Nmd2+UElcJHpoRC4iIuJiukcuIiLiYkrkIiIiLqZELiIi4mJK5CIiIi6mRC4iIuJiSuQiIiIupkQuIiLiYkrkIiIiLvb/AQJJOi/dEI0UAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for dat in data:\n",
    "    wav_deets = FWHM(np.array(dat[1]['Wavelength']), np.array(dat[1]['Transmission']))\n",
    "    depth = average_depths['5s'][average_depths['band'] == dat[0]]\n",
    "    #print(depth)\n",
    "    plt.plot([wav_deets[0],wav_deets[1]], [depth,depth], label=dat[0])\n",
    "    \n",
    "plt.xlabel('Wavelength ($\\AA$)')\n",
    "plt.ylabel('Depth')\n",
    "plt.xscale('log')\n",
    "plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
    "plt.title('Depths on {}'.format(FIELD))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### IV.c - Depth vs coverage comparison\n",
    "\n",
    "How best to do this? Colour/intensity plot over area? Percentage coverage vs mean depth?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Depths (5 $\\\\sigma$) vs coverage on SPIRE-NEP')"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfIAAAEYCAYAAACnTnQ0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XtcVHX6B/DPg5qKjIhCiCICXlA00CTN8ppmWmtpllm5WXlZ2ywzu1jtbr9fu221a62Z/VatVbup6ZalpVmRl9UtDVMQL2ggKgKKijiKVsLz++OcsQEHuQwzw4HP+/WalzNnvud7ni9Qz5zv98x5RFVBRERE1uTn6wCIiIio6pjIiYiILIyJnIiIyMKYyImIiCyMiZyIiMjCmMiJiIgsjImciIjIwpjIiYiILIyJnCAimSIyuJr6eklEHquOvlz0vVVEuniibyIiq2Iir4HMxHpOROwickpE/isik0WkWn5f1Zm4S/UbAuA+APOctq0XkfMicsZ8pLlxiJkAXnA3TvIMEelj/q0WiMhJEdksIteY7zn+ps+IyFERWSgiAU77XvybLNU2V0QWObd10cbxmFNGXJnmMZs4bZsgIusr2l958RP5EhN5zTVcVW0A2gJ4GcDTAP7l25DKdT+A1ap6rtT2KaoaYD5i3Oh/JYCBIhLmRh81mojU93UMVSEiTQF8BuANAM0BtAbwvwB+cmo2XFUDAFwN4BoAf7hMl4623QB0B/BMWW2cHlMu0199AFPLGUZ5/VUmfiKvYSKv4VS1QFVXArgLwDgR6QoAItJKRD4SkTwROSAijzr2Mc8enhGR3SKSb549NDLfew9ABIBV5tnFU+Zu3UQkxTyb+tDR3tznaRE5Ys4QpInIoDLCHQZgQ1XHKiINRORFM/5fRETNR7L5szgPYBuAIS72nSEi/y617XURmV3JMUBE2ojIx+bP9oTjzExEOpszDKdEZJeI3FqRY5uvy/t9PS0iKQDOikh9s890M97dIjLSqf3VIrLdfG+5+fv6S0WO5WKsLsfkFNcTZf1dlNIRAFR1iaoWqeo5Vf1SVVNKN1TVIwDWAOhaVlxObXMBrIWR0N3xdwBPiEgzN/upVPxE3sBEbhGquhVAFoC+YkyxrwKQDOPMZxCAx0TkJqdd7gVwE4B2MP4n+wezn98COIRfzz7+ZrYfDWAogCgAcTDOriEiMQCmALjGnCG4CUBmGWFeBcDV1PlLInJcjKnWAZcZ5l/MsfQF0AxAIoAVAEY6tdkDIN7FvksA3GyeGUJE6pljWlyZMZj7fQbgIIBIGD/fpSLSAMbP/EsAVwJ4BMAHZt9lHtt8XZHf190AbgHQTFUvAEg3fw6BMM5s3xeRMBG5wvyZLIJx5rvE+edTwWM52l5uTA4u/y5c2AegSETeEZFhIhJURjuISBsANwPYXlYbp7bhMD4g/lhe23IkAVgP4Ak3+6lU/EReoap81LAHjCQz2MX27wA8B6AXgEOl3nsGwEKn/Sc7vXczgPSy+jdfj3V6/TcAc83n7QEcAzAYQINy4v4FQKdS23oBsAFoCGAcADuAdi72tQE4B6CD07aHAKwv1e5FAAvKOP4mAPeZz290jLmSY+gNIA9A/VLb+wLIBeDntG0JgP+53LGdfgbl/b4eLCeuHQBuA9APwBEAUmrcf6nIsSo5pjL/LsqIsTOMDxhZAC7AWAoJderrDIBTMD4k/R+Axq7+Jp3a2gEojA90zVz8N+Loz/GYeLn/nmCcQRcACAEwwflvq7z+youfDz58+eAZubW0BnASxrp5K3M69JSInALwLIBQp7aHnZ4fBNCqnL5znZ4XAggAAFX9EcBjAP4HwDERWSoiZfWVDyMhX6SqW1TVrqo/qeo7ADbD+GBRWj8AGaq632lbUKm4YPZ/qozjL4ZxZgsA95ivKzuGNgAOqnFW7KwVgMOqWuy07SCM30mZxzZV9vcFEblPRHY4te8KINiM44iqahn7VuRYFR0TUMbfhSuqukdV71fVcDPeVgBmOTUZoarNVLWtqv5eL72WAqXa2gAMANAJxthdtWnm9HhLRO6VXy9WW1MqvlQYsy0zLnPMEv2V8X5F4ifyGiZyixDj6t/WMM6+DgM4UOp/OjZVdU6QbZyeRwDIdnpdqSL0qrpYVfvASBIK4JUymqbAXCu9XHcAxMX2EBgfBAAAIiIwpow/K9WuM4xpY1eWAxhgTseOhFMyrcQYDgOIkEsvOssG0EZKfnMgAsbZ8WWPjYr9vi7+TkSkLYC3YCwHtFDVZgBSYfzccgC0Nn8+Ds6/64ocq6JjqjJV3Qvj7NytdWRV3WD2M7OC7T/QXy9WG+aiyfMAJqLkhxUiS2Mir+FEpKmI/AbAUgDvq+pOAFsBnDYvkGosIvVEpKuZ7B0eFpFwEWkO44zsQ6f3jgKIruDxY0TkBhFpCOA8jOnvojKarwbQ32nfZiJyk4g0Mi/guhfGmfdaF/umArhaRLqJSGMAL8FIbh869dcQQA8AX7k6uKrmwVgHXQgjme2pwhi2wkiWL4tIEzP26wFsAXAWwFNiXJQ3AMBwGL+XMo/t1Gd5vy9nTcyx55nxP4BfE+K3ZuxTzJ/pbQB6VvFYlx1TZYhIJxGZbn6Qcawj3w1jOchdswDcKCLuXvDmmJ35EECZFwASWQ0Tec21SkTsMM6wngPwGoAHAEBVi2D8D7cbgAMAjgN4G8aFUQ6LYVzElGE+/uL03ksA/mBOvZZ38U9DGF9/Ow5jmvVKGB8MXHkXxkVfjc3XDczj5pn7PwJjevKSC+JUNQnG+vdqM96WAG5W1V+cmt0KY10zu/T+ThbDWA91PiOu8BicfrbtYVwUmAXgLlX92Tz+MLOf/4OxJr63nGNX9Pfl3H43gFdhJO2jMC4i3Gy+9zOA2wGMh7HEMBbGrMVPlT1WBcdUUXYY6/NbROQsjASeCmB6FfoqHWcejL+tP5Z6y/HNC8djRQW7fAHGh6XSqtofkU9JyaU2qg1EJBPABFX92gfH/iuAY6o6q9zGle97C4Dx5lonmcyfy1xVXejrWIjI+yx58wmquVS1rLP16ui7l6f6thIR6Q/ja37HYXzNMA7AFz4Nioh8homcyHpiACyDcQV5OoA7VDXHtyERka9wap2IiMjCeLEbERGRhdWqqfXg4GCNjIz0dRhERJaxbdu246oa4mYfV9avX/9tGF+T5Ali9SoGkHrhwoUJPXr0OOaqQa1K5JGRkUhKSvJ1GEREliEiB93to379+m+3bNmyc0hISL6fnx/Xa6tRcXGx5OXlxebm5r4N4+uil+AnJyIiclfXkJCQ00zi1c/Pz09DQkIKcJm7JDKRExGRu/yYxD3H/NmWma+ZyImIiCyMiZyIiMjCmMiJiMhrZq5NC/16z9ES5Y6/3nPUNnNtmqtSu5a0Zs2agNjY2M7169fvsXDhwiBPH4+JnMhDNm3ahAMHDpTYduDAAWzatMlHERH5XreIZoWPL9sR7UjmX+85ant82Y7obhHNCn0dW3WJjo7+eeHChZnDhw8/4Y3jMZETeUjr1q2xfPnyi8n8wIEDWL58OVq3ZilsqrsGdw61vza6W8bjy3ZE/++qXa0eX7Yj+rXR3TIGdw61u9334MHtunTp0rl9+/ZdZs6cGQwA/v7+3SdOnBgeGxvbuXfv3h2zs7PL/Nr1hg0b/Dt27BjbrVu3Tr/73e/CO3To0AUAZs+e3WLQoEHt+vbt2yEyMrLr9OnTwxz7zJkzp0XHjh1jY2JiYkeMGBEFADExMT/36tXrnJ9f+Sm2qKgIY8eOjWjfvn2XgQMHtu/fv3/7yp7FM5ETeUhUVBTuvPNOLF++HN988w2WL1+OO++8E1FRUb4OjcinBncOtY+6Ojxv4ebMsFFXh+dVRxIHgA8++CBz165de3bs2LF73rx5obm5ufXOnTvnd/XVVxfu3r17z/XXX2+fMWNGq7L2nzBhQtSbb755cMeOHXvr1atX4ir8lJSUJsuXL89ITU3dtXLlyuYbN270T0pKajRz5sywDRs27EtLS9s9b968Q5WN+d133w06fPjwFWlpabveeeedzO3btwdUtg8mciIPioqKQkJCAjZu3IiEhAQmcSIY0+kf/ZAV8sD1kTkf/ZAVUnrNvKpeeeWV0JiYmNgePXp0zs3NbbBr165Gfn5+mDBhwkkAePDBB09s3brVZaI8fvx4vbNnz/rdeOONZwFg3LhxJ53f79Onz+mWLVsWBQQE6C233JK/fv36gLVr1zYdPnx4flhY2AUACA0NLapszP/5z38Cbr/99vx69eohIiLiwrXXXlvpDzVM5EQedODAASQlJaFfv35ISkq6ZM2cqK5xrIm/NrpbxvPDu2Q7ptndTeafffaZbcOGDbakpKS9aWlpuzt37nzu3Llzl+Q4EXG5f3kFxErvJyJQVYiIW9+fr47CZUzkRB7iWBO/8847ccMNN1ycZmcyp7psx6FT/s5r4o418x2HTvm70++pU6fqBQYGFtlstuLt27c3Sk5ObgIAxcXFcKw5L1q0qEXPnj1dnvGGhIQUNWnSpDgxMbEJALz33nvNnd/ftGlT06NHj9Y7c+aMrF69uln//v3PDB069PTKlSub5+bm1gOAo0eP1qts3H379j3zySefBBUVFeHw4cP1t2zZUukPNLXqXutENcmRI0dKrIk71syPHDnCKXaqs564KeZo6W2DO4fa3V0nHzVqVMH8+fNDOnbsGNuuXbvz8fHxZwGgcePGxbt27WrcpUuXljabrejjjz/OKKuPefPmZU6ePLmtv79/8fXXX2+32WwXp8oTEhLO3HXXXVGZmZmNRo0adaJfv36FADB9+vScvn37dvLz89OuXbsWfvTRR5kbNmzwHz16dPvTp0/XS0xMbPbiiy+2+vHHH3e5Oua4cePyv/76a1vHjh27REVFnY+Pjz/brFmzSk3R16p65AkJCcqiKUREFSci21Q1wZ0+kpOTM+Pj449XV0zVyd/fv3thYeH2irQtKCjwCwwMLAaAZ599tmVOTk6DhQsXHp49e3aLpKSkJu+++26lL2arzHFzc3PrXXPNNZ03b968NyIi4oJzm+Tk5OD4+PhIV/vzjJyIiAjAsmXLAl999dWwoqIiad269U+LFy/O9MZxb7zxxg6nT5+u98svv8iTTz6ZUzqJl4eJnIiIai1XZ+O//e1vI77//vsSV68/9NBDR6dOnXpi4sSJ+aXbP/rooycAuHVzl61btza+7777SqypXXHFFcUpKSl7t27dmuZO30zkRERUp7z33nsemSK/nJ49e57bu3fvbk/0zavWiYiILIyJnIiIyMKYyImIiCyMiZyIiMjCmMiJiMh7Ev8cirQ1Je9elrbGhsQ/sx55FXkskYtIGxFZJyJ7RGSXiEw1tzcXka9EZL/5r8tBisg4s81+ERnnqTiJiMiLwhMKsWJy9MVknrbGhhWToxGeUOfrkV+4UKmvj1/kyTPyCwCmq2pnANcCeFhEYgHMAJCoqh0AJJqvSxCR5gCeB9ALQE8Az5eV8ImIyEJihtkxcm4GVkyOxpoZrbBicjRGzs1AzLA6WY/8s88+s/Xq1avj8OHDo2JiYrpUZdweS+SqmqOqP5jP7QD2AGgN4DYA75jN3gEwwsXuNwH4SlVPqmo+gK8ADPVUrERE5EUxw+yIvzsPW/4Zhvi786ojiQPWrEfu6Pvvf//7kfT0dJf3Yy+PV9bIRSQSQHcAWwCEqmoOYCR7AFe62KU1gMNOr7PMba76niQiSSKSlJeXV51hExGRJ6StsSF5SQh6PZSD5CUhl6yZV5EV65EDQFxc3NlOnTr9XJV9AS8kchEJAPARgMdU9XRFd3OxzWV1F1Wdr6oJqpoQEhJS1TCJiMgbHGviI+dmYNjL2Ren2d1M5latRw4A/v7+xe7s79FELiINYCTxD1T1Y3PzUREJM98PA3DMxa5ZANo4vQ4HkO3JWImIyAuykvxLrIk71syzkupkPfLq4LF7rYvx8eVfAPao6mtOb60EMA7Ay+a/n7rYfS2Avzpd4DYEwDOeipXIE7Z++m+0bNcREV3jLm47lJqC3PR96HnbHT6MjMiHBv3xknrkiBlmd3ed3Kr1yKuDx+qRi0gfAP8BsBOAY9rgWRjr5MsARAA4BOBOVT0pIgkAJqvqBHP/B832APCiqi4s75isR041yaHUFHw262X85rEZiOgad8lropqA9ch/5at65BXhk3rkqroJrte6AWCQi/ZJACY4vV4AYIFnoiPyvIiucfjNYzPw2ayXET/kZiR/uZpJnKgG81U9cnexjCmRh8zdkI648FaIH3IzvvtoKa4dNQZZjVth9YZ0TO7fztfhEdUJVqhH7k6/ABM5kcfEhQfiz/NXYuixr3DtqDFIWvM5vtj2M/446VZfh0ZUp9W2euRM5EQeEn4uG0OPfYUvrrwRRc2uwVdX/oyhx75C+LkEAMG+Do+IagkmciIPyU3fh5HTn0FRdkPM/uZHPDroeoxslYDc9H1cJyeiasPqZ0Qe0vO2O5DVuBXe33IIj97QHu9vOYSsxq341TMiqlZM5EQe8t/045iyeDvm3NMdjw+JwZx7umPK4u34b3qN/JYOEVkUEzmRh6RkFWDOPd1xXTtjPfy6dsGYc093pGQV+DgyIt+Z/cPs0PWH15e4Hev6w+tts3+YXSfrkWdmZjYYOnRotDvHYyIn8pArWmxAff+SN5Gq75+BK1ps8FFERL4XFxJX+Nym56IdyXz94fW25zY9Fx0XElcn65FHRkb+8sUXX5R5t7mKYCIn8pCuOXvxxLqp2JqzFQCwNWcrnlg3FV1z3P7aKJFlDWgzwP5inxczntv0XPTLW19u9dym56Jf7PNixoA2A+pkPfK0tLQrHMepKiZyIg/p2f43mHk0D0+sm4o52+fgiXVTMfNoHnq2/42vQyPyqQFtBtiHtxue98GeD8KGtxueVx1JHLBuPXJ3MZETeUpUP/QcsQCj809iXso8jM4/iZ4jFgBR/XwdGZFPrT+83rYqfVXIvZ3vzVmVviqk9Jp5VVm1Hrm7mMiJPGhro0ZY1tSG3+UXYFlTG7Y2auTrkIh8yrEm/mKfFzNm9JyR7ZhmdzeZW7keubuYyIk8xLEmPvP4KUyJ+x1mHj9VYs2cqC5KyUvxd14Td6yZp+SlsB55FfHObkQekrp/lbEmbk6n94zqi5mfPIjU/avQM6ynr8Mj8olHr370knrkA9oMsLu7Ts565LUE65FTjbJpFtD66pJr4gc2Akd+APo85ru4iJywHvmvWI+ciEpylayj+vFiN6IaivXIiYiIahjWIyciIqplals9cl61TkREZGFM5ERERBbGRE5ERGRhTOREREQWxkRORERec2zWrFD7unUlbsdqX7fOdmzWrDpZj7w6MJETEZHXNI6PL8x+eka0I5nb162zZT89I7pxfHydrEdeHZjIiYjIa2wDB9pbvfJyRvbTM6Jz//rXVtlPz4hu9crLGbaBA+tkPfLHHnusVadOnWI7deoUe+WVV8bdcccdkZUdNxM5kYfYNxzG+fRTJbadTz8F+4bDPoqIqGawDRxoDxxxW17+u++FBY64La86kjhgzXrks2bNyt67d+/uzZs3pzVr1uzC1KlTj1W2DyZyIg9pEG7DycV7Libz8+mncHLxHjQIr5bSy0SWZV+3zlbwyachQff9Nqfgk09DSq+ZV5VV65EXFxfjjjvuiHr44YeP9u3bt9JLDLyzG5GHNGrXDM3v6YyTi/egSa8wnN2Sg+b3dEajds18HRqRzzjWxB3T6U1697ZXx/S6cz1ym81W3LNnzxir1COfPn16q7CwsJ+nTp1apTV1npETeVCjds3QpFcY7N8cRpNeYUziVOedS072d07ajjXzc8nJdbIe+ZIlSwLXr1/fdMGCBVVec+MZOZEHnU8/hbNbcmC7oQ3ObslBw3bNmMypTrvysccuqUduGzjQ7u46uVXrkc+aNSv02LFjDbp169YZAIYOHXpq1qxZ2ZUZO+uRE3mIY03cMZ1e+jVRTcB65L9iPXIiKuGXLHuJpO1YM/8ly85ETlQDsR45EZVg69/mkm2NOLVO5FWsR05ERFTLsB45ERER1RgeS+QiskBEjolIqtO2eBH5VkR2isgqEWlaxr6ZZpsdIsKr14iIiMrgyTPyRQCGltr2NoAZqnoVgBUAnrzM/gNVtZu7V1MSERHVZh5L5Kq6EcDJUptjAGw0n38FYJSnjk9ERDXPd5+mhx5IOV7ilqwHUo7bvvs0vdaUMfU2b6+RpwK41Xx+J4BLL+s1KIAvRWSbiEy6XIciMklEkkQkKS8vrxpDJSKi6hYaFViYuGh3tCOZH0g5bktctDs6NCqw1pQxre31yB8E8LCIbANgA/BzGe2uV9WrAQwz2/crq0NVna+qCaqaEBISUv0RExFRtYmKC7YPuj82I3HR7uj/LNvXKnHR7uhB98dmRMUFV0sFtJqgVtcjV9W9qjpEVXsAWAIgvYx22ea/x2Cspff0XpRERORJUXHB9phrW+alfJMVFnNty7zqSuJWrEc+YsSIqPfff//izSVuvfXWqA8++CCwMuP2aiIXkSvNf/0A/AHAXBdtmoiIzfEcwBAYU/JERFQLHEg5bkv7Ljck7obwnLTvckNKr5lXlRXrkU+cODFv0aJFLQDgxIkT9bZt2xYwevTogsr04cmvny0B8C2AGBHJEpHxAO4WkX0A9gLIBrDQbNtKRFabu4YC2CQiyQC2AvhcVb/wVJxEROQ9jjXxQffHZvQd3THbMc1eHcncivXIb7nlljMHDx5sdOTIkfr/+te/mt9yyy35DRo0qFQfHruzm6reXcZbr7tomw3gZvN5BoB4T8VFRES+c/RAgb/zmrhjzfzogQJ/d6bYrVyPfPTo0Sfefvvt5h999FHzBQsWZFZ2f97ZjYiIvOba29odLZ2wo+KC7dfe1u6S8qaVYdV65AAwefLk4/PmzQsFgISEhPOV3Z/3WiciIsuzaj1yAGjTps2Fdu3anR8+fPipqoydiZyIiCyvcePGunHjxv2u3nv99dezYVyXdVk9evQ4t2/fvt2AUY/c8WEAAIKDgy+4qkf+yCOPnHjkkUdKfM2sf//+hUePHk2paOx2u90vMzOz4fjx40vfRK1COLVOREQEox55p06dYjt06NDlv//9b8CLL76Y4+ljfvLJJ7aOHTt2mThx4rEWLVpU+mI5gGfkRERUi1mhHvmIESN2utM3EzmRhxw8OA+2pnFoHtT74raT+d/CfjoFbdv+zoeREdVtrEdORBViaxqH1NRHcTL/WwBGEk9NfRS2pnE+joyIahOekRN5SPOg3ujadTZSUx9F69b34MiRxejadXaJM3QiInfxjJzIg5oH9Ubr1vcgM3MOWre+h0mciKodEzmRB53M/xZHjixGZOQUHDmy+OI0O1FdtWnpu6Hp27aWuB1r+rattk1L32U98ipiIifyEMeaeNeus9EuetrFaXYmc6rLwjp0Klzz5qvRjmSevm2rbc2br0aHdehUI+qRP/LII61btmwZ5+/v393XsVQUEzmRh9hPp5RYE3esmdtPV/g+EUS1TrsePe3DHp6esebNV6PXLZrfas2br0YPe3h6Rrserm+d6m0jRow4tWXLlj2+jqMyeLEbkYe4+opZ86DeXCenOq9dj572Lv0G5f2wZmXY1cNuzamuJP7kk0+G/fvf/24eFhb2c4sWLS5079698IsvvmjWtWvXwu3btzc5c+ZMvfnz5x8YOHBgYUFBgd/48eMjUlJS/AHg2Wefzb7//vtPDRo06Gx5x3HYtWtXw3vuuSeqqKhIBg8eXDB//vxQV99b9zSekRMRkVelb9tq27UxMeTqYbfm7NqYGFJ6zbwqNm7c6L9q1aqgnTt37v7888/TU1JSmjjeKyws9Nu+ffve2bNnH5w0aVIUAMyYMSOsadOmRfv27du9b9++3bfcckulP0xMmTKlze9///tjqampe1q1avWLu2OoKiZyIiLyGsea+LCHp2cMvH9StmOa3d1kvn79+oBhw4adCggI0KCgoOIbb7zxYgGSe+655yQADBs27MyZM2f8jh8/Xm/jxo1Np02bdszRJiQkpNK3R92+fXvAgw8+eBIAJkyY4Nad39zBRE5ERF6Ts3+vv/OauGPNPGf/Xn93+r1cPfHL1BJ355A1BhM5ERF5TZ8x9x0tvSberkdPe58x97lVj3zAgAFn1q5dG1hYWCgFBQV+X3/9dTPHe0uWLAkCgLVr1wbYbLaiFi1aFA0YMOD0a6+9dqWjTV5eXqVriXfr1u3MokWLggBgwYIFzctr7ylM5EREZHn9+/cvHDp0aEFsbGyXm2++uV1cXNzZwMDAIgAICgoq6t69e6cpU6a0nTdvXiYAvPTSSzmnTp2q16FDhy4xMTGxq1evtgHA5MmTw0NDQ+POnz/vFxoaGvf444+3KuuYb7zxxuE33ngj9Kqrruqck5PTICAgoErVy9wll5uOsJqEhARNSkrydRhERJYhIttUNcGdPpKTkzPj4+OPV1dMVVVQUOAXGBhYbLfb/Xr37h0zd+7cg48//nibmTNnHu7Xr1+1f0/dbrf7NWnSpNjPzw/z588P+vDDD5snJiamV/dxACA5OTk4Pj4+0tV7/PoZERHVCmPHjm27f//+xj/99JOMGTPmRJ8+fTx6k5nNmzf7T506NUJV0bRp06JFixZlevJ4ZWEiJyKiWmHVqlUHSm/bunVrmrv9Pv300y0//fTTEmvgt91228lXXnklNy0tzSOlSSujQolcRBoCGAUg0nkfVX3BM2ERERHVDK+88kruK6+8kuvrOMpS0TPyTwEUANgG4CfPhUNERESVUdFEHq6qQz0aCREREVVaRb9+9l8RucqjkRAREVGlXfaMXER2AlCz3QMikgFjal0AqKrGeT5EIiKqLQrWZoZeEWErbNy5xcWbwpzbc8L28yG7f+BNkW7dFKauKu+M/DcAhgMYBqA9gCHma8d2IiKiCrsiwlZ4ctm+6HN7TtgAI4mfXLYv+ooIG+uRV9FlE7mqHlTVgwD+4njuvM07IRIRUW3RuHMLe/PRHTNOLtsXfWpVequTy/ZFNx/dMcP5DN2XqlqP/JdffFYne9D/AAAWtElEQVT8rMIXu3VxfiEi9QD0qP5wiIiotmvcuYW9ydVX5p3ZnB0WcH2rnOpK4t6uRz5q1KjIoKCgCzt37vSPi4srfOutt7KqYxyVVd4a+TMAngXQWEROw1gbB4CfAcz3cGxERFQLndtzwnb2h2MhAde3yjn7w7GQhu2b2d1N5s71yH/55Rfp1q1bbPfu3QuBX+uRr1mzJmDSpElR+/fv3+VcjxyoWtEUAEhPT2+0efPmffXr++7+apc9sqq+BOAlEXlJVZ/xUkxERFRLOdbEHdPpDds3s1fH9LpzPXIAWpF65EuXLs1wtKlKPXIAuP322/N9mcSBik+tPysitwPoA+Mq9v+o6ieeC4uIiGqjnw/Z/Z2TtmPN/OdDdn93Ermv6pEHBAQUu92Jmyr6PfI3AUwGsBNAKoDJIvKmx6IiIqJaKfCmyKOlE3bjzi3s7n71zBf1yGuKiiby/gBuUtWFqroQwM0ABngsKiIiokrwRT3ymqJC9chF5GMA08yvnUFE2gJ4WVXv9nB8lcJ65ERElcN65NZwuXrkFT0jbwFgj4isF5H1AHYDCBGRlSKy0tUOIrJARI6JSKrTtngR+VZEdorIKhFpWsa+Q0UkTUR+FJEZFYyRiIjqsLFjx7bt1KlTbFxcXOfhw4fne7oeeU1R0Yvd/lSFvhcBmAPgXadtbwN4QlU3iMiDAJ4E8EfnnczvqL8J4EYAWQC+F5GVqurzmq9ERFRz+aIeubt9V4cKJXIz8bYF0EFVvxaRxgDqq2qZVxiq6kYRiSy1OQbARvP5VwDWolQiB9ATwI+qmgEAIrIUwG0wZgGIiIi8qqbXI6/Q1LqITATwbwDzzE3hAKry9bNUALeaz+8E0MZFm9YADju9zjK3lRXbJBFJEpGkvLy8KoRERERkXRVdI38YwPUATgOAqu4HcOVl93DtQQAPi8g2ADYYd4grzdUX+8q8Ik9V56tqgqomhISEVCEkIiIi66roGvlPqvqz48vzIlIfl0muZVHVvTAqqEFEOgK4xUWzLJQ8Uw8HkF3ZYxEREdUFFT0j3yAijnuu3whgOYBVlT2YiFxp/usH4A8A5rpo9j2ADiISJSJXABgDwOWV8UREZC2JiYmhaWlpNudtaWlptsTExFBfxeSsMmVM//a3v4XMmTOnhTfiupyKJvIZAPJg3NntdwBWw0jEZRKRJQC+BRAjIlkiMh7A3SKyD8BeGGfZC822rURkNQCo6gUAU2BcCLcHwDJV3VXZgRERUc0THh5euGLFimhHMk9LS7OtWLEiOjw8vEZ8VawyZUyfeuqpvClTppzwdEzlqehV68Ui8gmAT1S1QleUXeZmMa+7aJsN425xjterYXxYICKiWiQmJsY+cuTIjBUrVkTHx8fnJScnh4wcOTIjJibG7VKm3i5j+vjjj7cKCAgoeuGFF9y6vay7yitjKgCeh3GGLOamIgBvqOoLXoiPiIhqmZiYGHt8fHzeli1bwnr16pVTHUncV2VMa4LyptYfg3G1+jWq2kJVmwPoBeB6EZnm8eiIiKjWSUtLsyUnJ4f06tUrJzk5OaT0mnlVOJcxDQoKKq5IGdNp06Ydc7SpahnTmqC8RH4fgLtV9eLdcswbtYw13yMiIqowx5r4yJEjM4YNG5btmGZ3N5n7qoxpTVBeIm+gqpfcCN9cJ2/gmZCIiKi2ysrK8ndeE3esmWdlZfm70y/LmJbN1Q1bKvIeERHRJQYNGnS09Jp4TEyMfdCgQW5dMMYypmW9aVzY5uoKPgHQSFVr1Fk5y5gSEVUOy5haw+XKmF72qnVVtexUAxER1S1jx45tu3///sY//fSTjBkz5gTLmBIREVkIy5gSERHRJWpFGVMiIiKqmZjIiYiILIyJnIiIyMKYyImIiCyMiZyIiLwmPf3V0LzjiSVux5p3PNGWnv6q5eqR1xRM5ERE5DVNA7sV7t79RLQjmecdT7Tt3v1EdNPAbjXiO9+VqUdeU/DrZ0RE5DUhwYPssbEzM3bvfiI6rOXteTm5H4fExs7MCAkeZLl65J06dYp1PM/MzGz00Ucf7bvlllvOuDuOymIiJyIirwoJHmQPa3l73uGsRWFtwu/PqY4k7ot65Hv37t0NAIsXLw589dVXWw4ePLjCHwKqExM5ERF5Vd7xRFtO7schbcLvz8nJ/TgkqPl1dneTuXM9cgBakXrkS5cuzXC0qWo98p07dzZ87rnnwtetW7evYcOGZRcv8SCukRMRkdc41sRjY2dmdOz4x2zHNHvpC+Aqyxf1yE+fPu03evTodv/85z8PRkZG/uJWZ25gIiciIq85XbDD33lN3LFmfrpgh+XqkY8ZMyby3nvvPT506FCvr4s7YyInIiKvaddu+tHS0+ghwYPs7dpNt1Q98n379l3xxRdfBL3//vvBnTp1iu3UqVPsxo0b3fowUlWXrUduNaxHTkRUOaxHbg1VrkdORERkFaxHTkREZGGsR05ERESXYD1yIiIi8hgmciIiIgtjIiciIrIwJnIiIiILYyInIiKveSkjJ/TL4wUlbsf65fEC20sZOaxHXkVM5ERE5DU9mvoXPrLnULQjmX95vMD2yJ5D0T2a+teI73xbsR45EzkREXnNkOBA+xudIzIe2XMo+o/7s1o9sudQ9BudIzKGBAdWSz3yqKioLtddd12H4cOHR/3pT38K7dmzZ8yDDz7Ypnv37p06dOjQZd26df6AcRe4O+64I7Jjx46xHTt2jF20aFEzABg0aNDZtm3bllsAJT8/369169ZX/fTTTwIAJ0+eLPHam/g9ciIi8qohwYH20S2D8t7KOh42MTw4pzqSuLfrkQcFBRX37t3bvmzZssDf/va3pxYsWND85ptvzvdFKVOekRMRkVd9ebzAtiw3P2RieHDOstz8kNJr5lXhXI88KCiouCL1yKdNm3bM0aYq9cgnTZqUt2jRohYA8P777wdPmjTJJ/ebZyInIiKvcayJv9E5IuPPHcKzHdPs7iZzX9QjHzJkyNmsrKyGn3/+eUBRUZFcc801593qsIo8lshFZIGIHBORVKdt3UTkOxHZISJJItKzjH2LzDY7RGSlp2IkIiLv2na60N95TdyxZr7tdKHl6pEDwJgxY0488MAD0WPHjvVZ9TdPnpEvAjC01La/AfhfVe0G4E/ma1fOqWo383GrB2MkIiIveiY67GjpNfEhwYH2Z6LDLFWP3GH8+PEnTp8+XX/8+PEn3YnfHR672E1VN4pIZOnNAJqazwMBZHvq+EREVLc8//zzua+99lq2ox75U089dfTDDz9scdddd+W/+eabR5zbBgYGFn/88ceZpfuYO3du1ty5c7MqeszExETb0KFD84ODgyu9xl5dvH3V+mMA1orITBizAdeV0a6RiCQBuADgZVX9xFsBEhGRNXm7Hvm4ceParFu3LvCzzz7b78njlMfbifwhANNU9SMRGQ3gXwAGu2gXoarZIhIN4BsR2amq6a46FJFJACYBQEREhKfiJiKiGs7b9cjfeeedwwAOu9u/u+RyV/q53bkxtf6ZqnY1XxcAaKaqKsblggWq2vQyXUBEFpl9/Lu84yUkJGhSUpLbcRMR1RUisk1VE9zpIzk5OeOqq67K9/Pz8/p3qOuC4uJi2blzZ1B8fHy0q/e9/fWzbAD9zec3ALhkOkJEgkSkofk8GMD1AHZ7LUIiIqqs1Ly8vMDi4mKv39WstisuLpa8vLxAAKlltfHY1LqILAEwAECwiGQBeB7ARACvi0h9AOdhTomLSAKAyao6AUBnAPNEpBjGB42XVZWJnIiohrpw4cKE3Nzct3Nzc7uC9yepbsUAUi9cuDChrAYenVr3Nk6tExFVTnVMrZNv8ZMTERGRhTGRExERWRgTORERkYUxkRMREVkYEzkREZGFMZETERFZGBM5kYecePttnP1uS4ltZ7/bghNvv+2jiIioNmIiJ/KQRl2vwpFp0y4m87PfbcGRadPQqOtVPo6MiGoTbxdNIaozmlzbC63/8Q8cmTYNQXePQf6SpWj9j3+gybW9fB0aEdUiPCMn8qAm1/ZC0N1jcPz//omgu8cwiRNRtWMiJ/Kgs99tQf6SpQj+/UPIX7L0kjVzIiJ3MZETeYhjTbz1P/6BkEcfvTjNzmRORNWJiZzIQ86n7iyxJu5YMz+futPHkRFRbcLqZ0REdRirn1kfz8iJiIgsjImciIjIwpjIiYiILIyJnIiIyMKYyImIiCyMiZyIiMjCmMiJiIgsjImciIjIwpjIiTzkh7UHkZWWX2JbVlo+flh70EcREVFtxERO5CFXRjbF2rdSLybzrLR8rH0rFVdGNvVxZERUm7AeOZGHhMcE4aaJXbH2rVR07dcaqRuP4KaJXREeE+Tr0IioFuEZOZEHhccEoWu/1khanYmu/VoziRNRtWMiJ/KgrLR8pG48goSbI5G68cgla+ZERO5iIifyEMea+E0Tu6LXrdEXp9mZzImoOjGRE3nIsczTJdbEHWvmxzJP+zgyIqpNeLEbkYdcfVPbS7aFxwRxnZyIqhXPyImIiCyMiZyIiMjCmMiJiIgsjImciIjIwpjIiYiILIyJnIiIyMI8mshFZIGIHBORVKdt3UTkOxHZISJJItKzjH3Hich+8zHOk3ESERFZlafPyBcBGFpq298A/K+qdgPwJ/N1CSLSHMDzAHoB6AngeRHhl2+JiIhK8WgiV9WNAE6W3gzAUccxEEC2i11vAvCVqp5U1XwAX+HSDwRERER1ni/u7PYYgLUiMhPGB4nrXLRpDeCw0+ssc9slRGQSgEkAEBERUb2REhER1XC+uNjtIQDTVLUNgGkA/uWijbjYpq46U9X5qpqgqgkhISHVGCYREVHN54tEPg7Ax+bz5TDWwEvLAtDG6XU4XE/BE9VYcw4exaZ8e4ltm/LtmHPwqI8iIqLayBeJPBtAf/P5DQD2u2izFsAQEQkyL3IbYm4jsoxuTf0xaVfmxWS+Kd+OSbsy0a2pv48jI6LaxKNr5CKyBMAAAMEikgXjSvSJAF4XkfoAzsNc3xaRBACTVXWCqp4UkT8D+N7s6gVVLX3RHFGN1ifIhvldIjFpVybGtQrGO9nHMb9LJPoE2XwdGhHVIqLqcunZkhISEjQpKcnXYRCV8EpGDv5x8CimtQ3F09Fhvg6HqAQR2aaqCb6Og6qOd3Yj8qBN+Xa8k30c09qG4p3s45esmRMRuYuJnMhDHGvi87tE4unosIvT7EzmRFSdmMiJPGTH6cISa+KONfMdpwt9HBkR1Sa+uCEMUZ0wpW3oJdv6BNl4sRsRVSuekRMREVkYEzkREZGFMZETERFZGBM5ERGRhTGRExERWViturObiOQBOOjrOCopGMBxXwfhZRxz3cAxW0NbVWXpSAurVYncikQkqa7dHpFjrhs4ZiLv4NQ6ERGRhTGRExERWRgTue/N93UAPsAx1w0cM5EXcI2ciIjIwnhGTkREZGFM5ERERBbGRO5BIjJURNJE5EcRmeHi/bYikigiKSKyXkTCnd6LEJEvRWSPiOwWkUhvxl5Vbo75byKyyxzzbBER70ZfeSKyQESOiUhqGe+LOZYfzTFf7fTeOBHZbz7GeS9q91R1zCLSTUS+NX/HKSJyl3cjrzp3fs/m+01F5IiIzPFOxFSnqCofHngAqAcgHUA0gCsAJAOILdVmOYBx5vMbALzn9N56ADeazwMA+Pt6TJ4cM4DrAGw2+6gH4FsAA3w9pgqMuR+AqwGklvH+zQDWABAA1wLYYm5vDiDD/DfIfB7k6/F4eMwdAXQwn7cCkAOgma/H48kxO73/OoDFAOb4eix81L4Hz8g9pyeAH1U1Q1V/BrAUwG2l2sQCSDSfr3O8LyKxAOqr6lcAoKpnVLXQO2G7pcpjBqAAGsH4ANAQQAMARz0esZtUdSOAk5dpchuAd9XwHYBmIhIG4CYAX6nqSVXNB/AVgKGej9h9VR2zqu5T1f1mH9kAjgGwxB3F3Pg9Q0R6AAgF8KXnI6W6iIncc1oDOOz0Osvc5iwZwCjz+UgANhFpAePM5ZSIfCwi20Xk7yJSz+MRu6/KY1bVb2Ek9hzzsVZV93g4Xm8o62dSkZ+VVZU7NhHpCeNDW7oX4/Ikl2MWET8ArwJ40idRUZ3ARO45rtZ3S3/X7wkA/UVkO4D+AI4AuACgPoC+5vvXwJiqvt9jkVafKo9ZRNoD6AwgHMb/FG8QkX6eDNZLyvqZVORnZVWXHZt5pvoegAdUtdhrUXlWWWP+PYDVqnrYxftE1aK+rwOoxbIAtHF6HQ4g27mBOb14OwCISACAUapaICJZALaraob53icw1t3+5Y3A3eDOmCcB+E5Vz5jvrYEx5o3eCNyDyvqZZAEYUGr7eq9F5Vll/h2ISFMAnwP4gzkFXVuUNebeAPqKyO9hXOtyhYicUdVLLgQlqiqekXvO9wA6iEiUiFwBYAyAlc4NRCTYnHoDgGcALHDaN0hEHOuHNwDY7YWY3eXOmA/BOFOvLyINYJyt14ap9ZUA7jOvar4WQIGq5gBYC2CIiASJSBCAIea22sDlmM2/iRUw1pKX+zbEaudyzKp6r6pGqGokjNmod5nEqbrxjNxDVPWCiEyB8T/negAWqOouEXkBQJKqroRxRvaSiCiMM8+HzX2LROQJAInmV7C2AXjLF+OoDHfGDODfMD6w7IQxJfmFqq7y9hgqS0SWwBhTsDmT8jyMC/WgqnMBrIZxRfOPAAoBPGC+d1JE/gzjww8AvKCql7uYqsao6pgBjIZx9XcLEbnf3Ha/qu7wWvBV5MaYiTyOt2glIiKyME6tExERWRgTORERkYUxkRMREVkYEzkREZGFMZETERFZGBM5kRMRaSkiS0UkXYyqc6tFpKOv4yIiKgsTOZHJ/M7+CgDrVbWdqsYCeBZGwYvqPpYV7p1PRBbARE70q4EAfjFv8AEAMG9WssksXJMqIjsddbRF5EMRudnRVkQWicgoEalntv/erE39O/P9ASKyTkQWw7jxDUTkExHZZtbonuTU13gR2SdGzfa3HHWsRSRERD4y+/5eRK73yk+GiGos3tmN6FddYdxFr7TbAXQDEA8gGMD3IrIRRpnWuwCsNm8/OgjAQwDGw7hF5zUi0hDAZhFxlLDsCaCrqh4wXz9o3uWtsdnvRzDKuP4RRv1rO4BvYFSNA4y61v9Q1U0iEgHjLnqdq+9HQERWw0ROVL4+AJaoahGAoyKyAUZVujUAZpvJeiiAjap6TkSGAIgTkTvM/QMBdADwM4CtTkkcAB4VkZHm8zZmu5YANjhu2Soiy2GUtgWAwQBijVUAAEBTEbGpqr36h01EVsBETvSrXQDucLHdVYlKqOp5EVkP4CYYZ+ZLnNo/oqoliqCIyAAAZ0u9Hgygt6oWmn01Kut4Jj+z/bnyh0NEdQHXyIl+9Q2AhiIy0bFBRK4BkA/gLnPtOwRG4Y+tZpOlMApk9MWv1cvWAnjIrOIGEekoIk1cHC8QQL6ZxDvBKNsKs+/+ZmW0+gBGOe3zJYApTvF1c2vERGR5PCMnMqmqmtPcs0RkBoDzADIBPAajlnQyjMpsT6lqrrnblwDeBbBSVX82t70NIBLAD+aV8HkARrg45BcAJotICoA0AN+ZcRwRkb8C2AKjpvVuAAXmPo8CeNPcpz6MCnKTq+UHQESWxOpnRDWQiASo6hnzjHwFjJKwK3wdFxHVPJxaJ6qZ/kdEdgBIBXAAwCc+joeIaiiekRMREVkYz8iJiIgsjImciIjIwpjIiYiILIyJnIiIyMKYyImIiCzs/wENAqxQWJE6DQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for dat in data:\n",
    "    wav_deets = FWHM(np.array(dat[1]['Wavelength']), np.array(dat[1]['Transmission']))\n",
    "    depth = average_depths['5s'][average_depths['band'] == dat[0]]\n",
    "    #print(depth)\n",
    "    coverage = np.sum(~np.isnan(depths['ferr_{}_mean'.format(dat[0])]))/len(depths)\n",
    "    plt.plot(coverage, depth, 'x', label=dat[0])\n",
    "    \n",
    "plt.xlabel('Coverage')\n",
    "plt.ylabel('Depth')\n",
    "#plt.xscale('log')\n",
    "plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
    "plt.title('Depths (5 $\\sigma$) vs coverage on {}'.format(FIELD))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (herschelhelp_internal)",
   "language": "python",
   "name": "helpint"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}