{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# SA13 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-25 18:19:01.563563\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 = 'SA13'\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_sa13_20180501.fits\n"
     ]
    }
   ],
   "source": [
    "TMP_DIR = os.environ.get('TMP_DIR', \"./data_tmp\")\n",
    "OUT_DIR = os.environ.get('OUT_DIR', \"./data\")\n",
    "SUFFIX = find_last_ml_suffix()\n",
    "#SUFFIX = \"20171016\"\n",
    "\n",
    "master_catalogue_filename = \"master_catalogue_{}_{}.fits\".format(FIELD.lower(), SUFFIX)\n",
    "master_catalogue = Table.read(\"{}/{}\".format(OUT_DIR, master_catalogue_filename))\n",
    "\n",
    "print(\"Depth maps produced using: {}\".format(master_catalogue_filename))\n",
    "\n",
    "ORDER = 10\n",
    "#TODO write code to decide on appropriate order\n",
    "\n",
    "field_moc = MOC(filename=\"../../dmu2/dmu2_field_coverages/{}_MOC.fits\".format(FIELD))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Remove sources whose signal to noise ratio is less than five as these will have been selected using forced \n",
    "# photometry and so the errors will not refelct the RMS of the map \n",
    "for n,col in enumerate(master_catalogue.colnames):\n",
    "    if col.startswith(\"f_\"):\n",
    "        err_col = \"ferr{}\".format(col[1:])\n",
    "        errs = master_catalogue[err_col]\n",
    "        fluxes = master_catalogue[col]\n",
    "        mask = fluxes/errs < 5.0\n",
    "        master_catalogue[col][mask] = np.nan\n",
    "        master_catalogue[err_col][mask] = np.nan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## I - Group masterlist objects by healpix cell and calculate depths\n",
    "We add a column to the masterlist catalogue for the target order healpix cell <i>per object</i>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#Add a column to the catalogue with the order=ORDER hp_idx\n",
    "master_catalogue.add_column(Column(data=coords_to_hpidx(master_catalogue['ra'],\n",
    "                                                       master_catalogue['dec'],\n",
    "                                                       ORDER), \n",
    "                                   name=\"hp_idx_O_{}\".format(str(ORDER))\n",
    "                                  )\n",
    "                           )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Convert catalogue to pandas and group by the order=ORDER pixel\n",
    "\n",
    "group = master_catalogue.group_by([\"hp_idx_O_{}\".format(str(ORDER))])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#Downgrade the groups from order=ORDER to order=13 and then fill out the appropriate cells\n",
    "#hp.pixelfunc.ud_grade([2599293, 2599294], nside_out=hp.order2nside(13))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## II Create a table of all Order=13 healpix cells in the field and populate it\n",
    "We create a table with every order=13 healpix cell in the field MOC. We then calculate the healpix cell at lower order that the order=13 cell is in. We then fill in the depth at every order=13 cell as calculated for the lower order cell that that the order=13 cell is inside."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "depths = Table()\n",
    "depths['hp_idx_O_13'] = list(field_moc.flattened(13))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=10</i>\n",
       "<table id=\"table140650591834576-197676\" 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>177665210</td></tr>\n",
       "<tr><td>1</td><td>177656680</td></tr>\n",
       "<tr><td>2</td><td>177665538</td></tr>\n",
       "<tr><td>3</td><td>177665204</td></tr>\n",
       "<tr><td>4</td><td>177656718</td></tr>\n",
       "<tr><td>5</td><td>177665539</td></tr>\n",
       "<tr><td>6</td><td>177656666</td></tr>\n",
       "<tr><td>7</td><td>177665205</td></tr>\n",
       "<tr><td>8</td><td>177665206</td></tr>\n",
       "<tr><td>9</td><td>177665207</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(\"$('#table140650591834576-197676').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",
       "    $('#table140650591834576-197676').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=\"table140650591834632-55216\" 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>177665210</td><td>2776018</td></tr>\n",
       "<tr><td>1</td><td>177656680</td><td>2775885</td></tr>\n",
       "<tr><td>2</td><td>177665538</td><td>2776024</td></tr>\n",
       "<tr><td>3</td><td>177665204</td><td>2776018</td></tr>\n",
       "<tr><td>4</td><td>177656718</td><td>2775886</td></tr>\n",
       "<tr><td>5</td><td>177665539</td><td>2776024</td></tr>\n",
       "<tr><td>6</td><td>177656666</td><td>2775885</td></tr>\n",
       "<tr><td>7</td><td>177665205</td><td>2776018</td></tr>\n",
       "<tr><td>8</td><td>177665206</td><td>2776018</td></tr>\n",
       "<tr><td>9</td><td>177665207</td><td>2776018</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(\"$('#table140650591834632-55216').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",
       "    $('#table140650591834632-55216').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=\"table140650967946296-567138\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>idx</th><th>hp_idx_O_13</th><th>hp_idx_O_10</th><th>ferr_ap_ukidss_j_mean</th><th>f_ap_ukidss_j_p90</th><th>ferr_ukidss_j_mean</th><th>f_ukidss_j_p90</th><th>ferr_ap_90prime_g_mean</th><th>f_ap_90prime_g_p90</th><th>ferr_90prime_g_mean</th><th>f_90prime_g_p90</th><th>ferr_ap_90prime_r_mean</th><th>f_ap_90prime_r_p90</th><th>ferr_90prime_r_mean</th><th>f_90prime_r_p90</th><th>ferr_ap_mosaic_z_mean</th><th>f_ap_mosaic_z_p90</th><th>ferr_mosaic_z_mean</th><th>f_mosaic_z_p90</th></tr></thead>\n",
       "<thead><tr><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th>uJy</th><th>uJy</th><th></th><th></th><th>uJy</th><th>uJy</th><th></th><th></th><th>uJy</th><th>uJy</th></tr></thead>\n",
       "<tr><td>0</td><td>177655159</td><td>2775861</td><td>6.322884</td><td>607.8445449829104</td><td>11.099306</td><td>832.7705078125</td><td>0.21581411</td><td>6.456100940704344</td><td>0.35495663</td><td>38.10558996200552</td><td>0.30403265</td><td>12.514632987976078</td><td>0.44061643</td><td>18.10475807189941</td><td>0.76093984</td><td>27.51869926452637</td><td>0.46060133</td><td>23.69996738433838</td></tr>\n",
       "<tr><td>1</td><td>177655134</td><td>2775861</td><td>6.322884</td><td>607.8445449829104</td><td>11.099306</td><td>832.7705078125</td><td>0.21581411</td><td>6.456100940704344</td><td>0.35495663</td><td>38.10558996200552</td><td>0.30403265</td><td>12.514632987976078</td><td>0.44061643</td><td>18.10475807189941</td><td>0.76093984</td><td>27.51869926452637</td><td>0.46060133</td><td>23.69996738433838</td></tr>\n",
       "<tr><td>2</td><td>177655135</td><td>2775861</td><td>6.322884</td><td>607.8445449829104</td><td>11.099306</td><td>832.7705078125</td><td>0.21581411</td><td>6.456100940704344</td><td>0.35495663</td><td>38.10558996200552</td><td>0.30403265</td><td>12.514632987976078</td><td>0.44061643</td><td>18.10475807189941</td><td>0.76093984</td><td>27.51869926452637</td><td>0.46060133</td><td>23.69996738433838</td></tr>\n",
       "<tr><td>3</td><td>177655165</td><td>2775861</td><td>6.322884</td><td>607.8445449829104</td><td>11.099306</td><td>832.7705078125</td><td>0.21581411</td><td>6.456100940704344</td><td>0.35495663</td><td>38.10558996200552</td><td>0.30403265</td><td>12.514632987976078</td><td>0.44061643</td><td>18.10475807189941</td><td>0.76093984</td><td>27.51869926452637</td><td>0.46060133</td><td>23.69996738433838</td></tr>\n",
       "<tr><td>4</td><td>177655164</td><td>2775861</td><td>6.322884</td><td>607.8445449829104</td><td>11.099306</td><td>832.7705078125</td><td>0.21581411</td><td>6.456100940704344</td><td>0.35495663</td><td>38.10558996200552</td><td>0.30403265</td><td>12.514632987976078</td><td>0.44061643</td><td>18.10475807189941</td><td>0.76093984</td><td>27.51869926452637</td><td>0.46060133</td><td>23.69996738433838</td></tr>\n",
       "<tr><td>5</td><td>177655158</td><td>2775861</td><td>6.322884</td><td>607.8445449829104</td><td>11.099306</td><td>832.7705078125</td><td>0.21581411</td><td>6.456100940704344</td><td>0.35495663</td><td>38.10558996200552</td><td>0.30403265</td><td>12.514632987976078</td><td>0.44061643</td><td>18.10475807189941</td><td>0.76093984</td><td>27.51869926452637</td><td>0.46060133</td><td>23.69996738433838</td></tr>\n",
       "<tr><td>6</td><td>177655133</td><td>2775861</td><td>6.322884</td><td>607.8445449829104</td><td>11.099306</td><td>832.7705078125</td><td>0.21581411</td><td>6.456100940704344</td><td>0.35495663</td><td>38.10558996200552</td><td>0.30403265</td><td>12.514632987976078</td><td>0.44061643</td><td>18.10475807189941</td><td>0.76093984</td><td>27.51869926452637</td><td>0.46060133</td><td>23.69996738433838</td></tr>\n",
       "<tr><td>7</td><td>177655166</td><td>2775861</td><td>6.322884</td><td>607.8445449829104</td><td>11.099306</td><td>832.7705078125</td><td>0.21581411</td><td>6.456100940704344</td><td>0.35495663</td><td>38.10558996200552</td><td>0.30403265</td><td>12.514632987976078</td><td>0.44061643</td><td>18.10475807189941</td><td>0.76093984</td><td>27.51869926452637</td><td>0.46060133</td><td>23.69996738433838</td></tr>\n",
       "<tr><td>8</td><td>177655132</td><td>2775861</td><td>6.322884</td><td>607.8445449829104</td><td>11.099306</td><td>832.7705078125</td><td>0.21581411</td><td>6.456100940704344</td><td>0.35495663</td><td>38.10558996200552</td><td>0.30403265</td><td>12.514632987976078</td><td>0.44061643</td><td>18.10475807189941</td><td>0.76093984</td><td>27.51869926452637</td><td>0.46060133</td><td>23.69996738433838</td></tr>\n",
       "<tr><td>9</td><td>177655167</td><td>2775861</td><td>6.322884</td><td>607.8445449829104</td><td>11.099306</td><td>832.7705078125</td><td>0.21581411</td><td>6.456100940704344</td><td>0.35495663</td><td>38.10558996200552</td><td>0.30403265</td><td>12.514632987976078</td><td>0.44061643</td><td>18.10475807189941</td><td>0.76093984</td><td>27.51869926452637</td><td>0.46060133</td><td>23.69996738433838</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(\"$('#table140650967946296-567138').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",
       "    $('#table140650967946296-567138').dataTable({\n",
       "        order: [],\n",
       "        pageLength: 50,\n",
       "        lengthMenu: [[10, 25, 50, 100, 500, 1000, -1], [10, 25, 50, 100, 500, 1000, 'All']],\n",
       "        pagingType: \"full_numbers\",\n",
       "        columnDefs: [{targets: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18], type: \"optionalnum\"}]\n",
       "    });\n",
       "});\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for col in master_catalogue.colnames:\n",
    "    if col.startswith(\"f_\"):\n",
    "        errcol = \"ferr{}\".format(col[1:])\n",
    "        depths = join(depths, \n",
    "                      group[\"hp_idx_O_{}\".format(str(ORDER)), errcol].groups.aggregate(np.nanmean),\n",
    "                     join_type='left')\n",
    "        depths[errcol].name = errcol + \"_mean\"\n",
    "        depths = join(depths, \n",
    "                      group[\"hp_idx_O_{}\".format(str(ORDER)), col].groups.aggregate(lambda x: np.nanpercentile(x, 90.)),\n",
    "                     join_type='left')\n",
    "        depths[col].name = col + \"_p90\"\n",
    "\n",
    "depths[:10].show_in_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## III - Save the depth map table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "depths.write(\"{}/depths_{}_{}.fits\".format(OUT_DIR, FIELD.lower(), SUFFIX), overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## IV - Overview plots\n",
    "\n",
    "### IV.a - Filters\n",
    "First we simply plot all the filters available on this field to give an overview of coverage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'90prime_g', '90prime_r', 'mosaic_z', 'ukidss_j'}"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tot_bands = [column[2:] for column in master_catalogue.colnames \n",
    "             if (column.startswith('f_') & ~column.startswith('f_ap_'))]\n",
    "ap_bands = [column[5:] for column in master_catalogue.colnames \n",
    "            if column.startswith('f_ap_') ]\n",
    "bands = set(tot_bands) | set(ap_bands)\n",
    "bands"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Passbands on SA13')"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAesAAAEgCAYAAACU3FvWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xl8VNX5+PHPM1t2AglhDSQoYZNdREVadwR3LbVqa9XWSr9KRcFWf61FBW2rfsWlVb+gVat1wVqlCkrV1iJV2VG2gCyyhBAIIWRfZjm/P2YmhJBlApnMnczzfr14zXLv3PtMyOSZ55xzzxFjDEoppZSyLlukA1BKKaVU8zRZK6WUUhanyVoppZSyOE3WSimllMVpslZKKaUsTpO1UkopZXGarFWHIyIPiMhfO9q5lFKxS5O1ahcislNEqkSkXET2i8hLIpIc6bg6AhH5qYhsFpGywM92kYikNNjnARExIjK2wfM9ReQ9EckPbM9usP1REdkjIqUisktEfhP+d6SUakiTtWpPlxljkoHRwGnAfRGOJ+qJyNnA74DrjDEpwGDgrQb7CHADcAi4scEhfMBi4HtNnOLPwCBjTCdgHHC9iFzddu9AKRUKTdaq3Rlj9gIfAkMBRORmEckNVIY7RGRKcF8R6SoiC0XksIgcEpGlImILbLtHRPYGXrdFRM6vd5p4EZkf2LZGREbUO+a9IrI9sG2TiFxVb9tNIvJfEflfESkWkW9FZFK97f1EZEngtR8DXettixeRv4pIUSDelSLSvbGfgYgMFpH/BPbbKCKX19v2sog8E6iQy0RkuYic3MSP8zTgS2PM2sDP9pAx5i/GmLJ6+3wH6AVMA64VEVe9/4v9xphngZWNHdwYs8UYU1HvKR/Qv4lYlFJhoslatTsR6QNcDKwNPHUAuBToBNwMPCEiowPbZgB5QAbQHfg1YERkIDAVOC1QUV4E7Kx3miuAvwFpwOvAAhFxBrZtx5/AUoEHgb+KSM96rz0d2II/ET8K/DlQnRI41urAttkcXaneGDhmHyAd+DlQ1cj7dwLvAx8B3YBfAK8F3lPQdYHYugDbgIcbHidgOXCRiDwoImeJSFwj+9wYON/8wONLmzhWowJfbsrx/z8k4f8ZKKXakSZr1Z4WiMhh4L/AEvzNtxhjFhljthu/JfiT2HcCr3EDPYEsY4zbGLPU+Ce09wJxwBARcRpjdhpjttc712pjzNvGGDcwB4gHzgic72/GmHxjjM8YMx/YCtTvy91ljHneGOMF/hI4f3cR6Yu/kv2tMabGGPMZ/iQY5MafpPsbY7zGmNXGmNJGfg5nAMnAH4wxtcaYfwML8SfooHeMMSuMMR7gNWBkYz9QY8xS4Gr8XQuLgCIRmSMidgARSQS+D7we+Fm8zbFN4c0yxvwBSAmc41WgpDWvV0qdOE3Wqj1daYzpbIzJMsbcZoypAhCRSSKyLNDMfRh/1R1sXn4Mf2X5UaCJ/F4AY8w24E7gAeCAiLwpIr3qnWtP8I4xxoe/KuwVON+PReSrQBP0YfzN8V3rvbag3msrA3eTA68vbtAsvKve/VeBfwJvBgZsPVqvmq+vF7AnEFf94/RuLAagMnD+RhljPjTGXIa/FeEK4CbglsDmqwAP8EHg8WvAJBHJaOp4TZzDBJraq/BX/EqpdqTJWkVUoNn278D/At2NMZ3xJxYBMMaUGWNmGGNOAi4Dpgf7po0xrxtjxgNZgAEeqXfoPvXOYQMygXwRyQKex9+Enh4434bg+VqwD+giIkn1nusbvBOo/B80xgzBPxjrUuDHjRwnH+gT7Huvd5y9IcTQpEBLwb+AfxMYD4C/ik4GdotIAf6uASdHV/Gt4QCa6j9XSoWJJmsVaS78zdmFgCcwmGtCcKOIXCoi/QN9xqX4m7+9IjJQRM4LJPtq/BWft95xTxWRq0XEgb8CrwGW4e9zNYHzISI3cySxNcsYswtYBTwoIi4RGY//C0Qw1nNFZFigCboUf7O4t5FDLQcqgF+JiFNEzgkc581Q4qhPRK4QkWtFpIv4jQXOBpaJSG/gfPxfGkYG/o3A/6XmxnrHiMf/fwAQF3iMiNhEZEqDY98O/Ku1cSqlTowmaxVRgVHLd+C/3KgYuB54r94uOcAnQDnwJfCsMeY/+JPLH4CD+JuMu+EffBb0D+AHgWPeAFwdqHw3AY8HjrUfGAZ83oqQr8c/AO0QcD/wSr1tPfD3CZcCufj75Y+ZMMUYUwtcDkwKxP8s8GNjzOZWxBFUDPwMf797aeB8jxljXsP/vr8yxnxkjCkI/gOeBoaLSPBLShX+ny/AZo4eFHcV/gF5ZYFj/zHwTynVjsQ/VkcppZRSVqWVtVJKKWVxmqyVUkopi9NkrZRSSlmcJmullFLK4jRZK6WUUhbniHQArdW1a1eTnZ0d6TCUUiqqrF69+qAxplUz1zVyjG4Oh+MF/HMTaLHXtnzABo/Hc8upp556oOHGqEvW2dnZrFq1KtJhKKVUVBGRXS3v1TyHw/FCjx49BmdkZBTbbDa97rcN+Xw+KSwsHFJQUPAC/nkYjqLfjJRSSoVqaEZGRqkm6rZns9lMRkZGCU3MqKjJWimlVKhsmqjDJ/CzbTQvhy1Zi8iLInJARDY0sV1E5GkR2SYi6+qtX6yUUkqpesJZWb8MTGxm+yT88z7nALcCz4UxFqWUUh3A7Nmzu+Xk5JzSv3//U2bNmtUNYP/+/fZx48blZGVlDR03blxOYWGhvbXHHTVq1KC2j7bthC1ZG2M+w7/YQVOuAF4JrJO7DOgsIj3DFY9SSqnotnLlyvhXXnklY82aNbm5ubkbFy9e3Hn9+vVx999/f89zzjmnbNeuXRvOOeecspkzZ/YI9ZgejweAtWvXHs9COu0mkqPBewN76j3OCzy3LzLhKKWUCtUv3/66zzcFZYltecwBPVIqH5s8Yk9T29evX58wevTo8pSUFB/AWWedVTZ//vzOixcv7rxkyZItAFOmTCk6++yzBwJ7p0+f3mvHjh1xBQUFzn379rnuuOOOghkzZhxcuHBhyuzZs3t269bNvWnTpsTt27dvTExMHFVZWbl24cKFKQ8++GCvjIwM96ZNmxIvvvji4mHDhlU9++yz3WtqauTdd9/dfsopp9Tk5+c7br755qy9e/e6AObMmbN7woQJFY3FnZ+f75g8eXK/w4cPO0aOHFn5n//8p9Pq1atze/bs6Qn1ZxPJAWbSyHONDlwQkVtFZJWIrCosLAxzWErFNq/PS0FFAXtK97C1eCsLti2gsFI/dyryRo4cWbV8+fKUgoICe1lZme3jjz9O3bNnj6uoqMiRlZXlBsjKynIfOnSorhDNzc1N+OSTT7YuW7Zs82OPPdZr586dToB169YlPfbYY3u3b9++seF5Nm/enPDcc8/tyc3N3fj222+nf/PNN/Hr16/PveGGGw4+/vjj3QCmTJnSZ/r06fs3bNiQ++67727/+c9/nt1U3Pfee2+vs88+u2zTpk25V199dfG+fftcrX3vkays84A+9R5nAvmN7WiMmQfMAxgzZoyORFSqDXh9Xmp9tSQ4EgCodFfyxOon+Mf2f1DlqTpm/+sGXcedo+8k0XlsMeX1edlXsY8UVwqpcalhj11FXnMVcLiMHj26etq0aQXnnXfegMTERN+QIUMqHY7m09ikSZMOJycnm+TkZM+ZZ55ZunTp0qQuXbp4hw8fXjFo0KDaxl4zbNiwimDy79u3b82kSZNKAEaMGFG1ZMmSFIDPP/+809atWxOCrykvL7cXFxfbunTp4mt4vBUrViQvWLBgG8DkyZNLO3Xq5G3te49ksn4PmCoibwKnAyXGGG0CV6od7Cvfxw8W/oBKTyXXD76ezORMnl//PAUVBXRL7MbALgOp8dZwes/TWXNgDesK1/HG5jd4Y/MbLPnBEtLi0wDYcHADs76cxTfF3+A1/r8/g9MG89szfsuwjGEA1HprySvLQ0Tol9ovYu9ZdQx33XXXwbvuuusgwNSpU3tnZmbWpqene3bt2uXMyspy79q1y5mWllbXvCxydCNu8HFiYuIxSTUoLi6urii02WzEx8eb4H2v1ysAxhhWrVqVm5yc3GIBacyJ15hhS9Yi8gZwDtBVRPKA+wEngDHm/4APgIuBbUAlcHO4YlFKHe3vW/9OcU0xA7oM4KUNLwGQmZzJSxe9xJgeYxp9zWu5r/GHFX/g7PlnM6nfJJbvW86hav8Y0nP7nEtafBpltWUsL1jO9R9cz9D0oSQ4E9hwcENdpT5l+BSmjpraPm9SdUh79+519O7d27N161bXokWLOq9YsWLzt99+Gzd37tz03/3udwVz585Nnzhx4uHg/h9++GHnhx9+eF9paalt2bJlKU888cTeDRs2xJ9oHOPHjy995JFHus2ePXs/wBdffJEwbty4Y5ukgLFjx5a/+uqraQ8//HDBO++806m0tLTVo9XDlqyNMde1sN0At4fr/Eqppn2W9xmju43mL5P+wpZDWyitLWVUt1E4bE3/Sfjh4B+S3Smb2/51Gx9++yEAkwdM5pZht9A7uXfdfhXuCp5f9zyf7vmUgsoCJmZPZGDaQJbsWcLcdXMZ2W0k43uPD/t7VB3T5ZdffvLhw4cdDofDPPnkk7szMjK8Dz744L6rrrrq5KysrK69evWqXbBgwfbg/qNGjao4//zzc/Lz81133333vuzsbHdbJOt58+btueWWW/oOGDBgiNfrldNPP71s3Lhxuxvb9w9/+EP+5MmTTxoyZEiXM888szwjI8PduXPnVjWFS1uU5+1pzJgxRucGV+r4uX1uxv51LDcNvYlpo6cd9zEEaTa5H/Mar5sr/3ElCY4E/nbZ345pnlThJSKrjTGNN5uE6Ouvv945YsSIg20VU7hNnz69V3JysnfWrFn7IxlHVVWVOBwO43Q6+eSTT5KmTp2atXnz5k2N7fv11193HTFiRHbD56NuIQ+l1InZU7oHj/FwUupJx30Mp83Z+tfYndww5AYeXv4w3xR/w8C0gcd9fqWiybZt21zXXHPNyT6fD6fTaebOnbuztcfQZK1UjCmoKAA4qum6vVyYdSG/X/F7FmxbwD1j72n386vYMmfOnEavMAqXp556Kv25557rXv+50047rfzVV1/dnZub22glHSpN1krFmMIq/zXTGQkntLTxcUlPSOeMnmfw793/5len/UqbwlWHMm3atKJp06YVhePYuuqWUjHmYJW/yzE9IT0i55/UbxL5FfmsObCGGm9NRGJQKtposlYqxhysOkiiI7HRyU3aw4SsCSQ7k7lp8U1c8LcLKK0tjUgcSkUTTdZKxZiiqiK6JnSN2PkTnYn8/ju/p29KXw7XHGZ1weqIxaJUtNBkrVSMKXWX0snVKaIxnNPnHP568V8ByK9o1zFASkUlTdZKxZiK2gqSXEmRDoPOcZ1x2BwcqDwQ6VBUlBs7duzAzz777Jh+nbPPPrv/wYMHj5ktbPr06b1mzpzZveHzrbVz507nxIkTj/8ayFbQ0eBKxZhyd3nEBpfVJyKkxaVRXF0c6VBUB7VkyZJt4Tx+dna2e/HixTvCeY4gTdZKxZgKdwVJzshX1uDvv670VEY6DHU8FtzehwOb2naUYrchlVz5TJOreW3ZssV16aWX5mzdunUjwMyZM7uXl5fXVc5er5fvf//72ZmZmbVPP/10fu/evYetWrUqt2fPnp577rmnx/z587v26tWrNj093T1q1KhKgIceeqjbSy+9lGG3282AAQOqFy5cuGPRokXJM2bM6Av+L5VffPHF5sZW02oYTzhpslYqxpS7y0lxpUQ6DACSnElUuCsiHYbqANxut1x55ZX9hgwZUvXII48U1N+2dOnSxHfffTdt/fr1m9xuNyNHjhwSTNZPP/10j127dq1PSEgwwSbzxx9/vMfTTz+9a8KECRUlJSW25lboai+arJWKIcYY61XWbq2so1IzFXAk3HbbbVlXXnnloYaJGuDTTz9Nvvjiiw+npKT4ACZMmFC3KtfAgQOrrrrqqn6XX3754R/+8IeHAc4444zyu+++u88111xz6Lrrris++eSTI56sdYCZUjGkylOFz/hIdiZHOhQAEh3aDK5C53A4jM93JG9WV1fX5bAxY8aUL126tFNlZWWj0+I1NVvep59+uvX2228vXL16ddKIESOGuN1ufve73xW88MILu6qqqmzjxo0bvHbt2hNepetEabJWKobkl/svk9LKWkWjzMxMz6FDhxwFBQX2qqoq+ec//5ka3DZlypSDEyZMKLn00ktPdrvdR73uvPPOK1+0aFHn8vJyKS4utn388cedwd/HvX37dtdll11W9uyzz+aVlZXZS0pK7Bs3bowbO3Zs1cMPP1wwbNiwirZYUvNEaTO4UjHkxsU3ApDgSIhwJH6JjkSqPFWRDkNFibi4ODNjxox9Y8eOHZyZmVnTv3//6vrbH3jggf133XWX/eqrr+63YMGCb4PPjx8/vvKqq646NHTo0FN69+5dM3bs2HIAj8cj119/fb+ysjK7MUamTJmyv2vXrt4ZM2b0+uKLLzrZbDYzYMCAqsmTJ5c0FZOItMs605qslYohwak9W7MOdTg5bU5qfbWRDkNFkfvuu+/Afffd1+TF+U888UTdLDt79+5dH7z/yCOPFDTWn7169eotDZ/7y1/+ElJ//IEDBxypqaneUPY9UdoMrlQMSnREZl7whlx2F7VeTdYq+nz22WeJN9xww0lTp07d3x7ns8bXa6VUu4p3RLwLDvAna7fX3fKOSkXQihUrEn784x/3q/+cy+Xy7dy5c0N7xaDJWqkYFGePi3QIALhsLjzGg8/4sIk29ClrGjt2bNXmzZs3RTIG/XQoFYOyO2VHOgQAnHYngDaFK9UCTdZKxZCL+11M35S+dI7vHOlQgCMVvg4yU6p5mqyViiFun9syI8HB3wwOWlkr1RJN1krFELfPjdPmjHQYdVx2f7LWQWZKNU+TtVIxxOPzWCpZ1/VZazO4sphRo0YNinQM9WmyViqGuH3uugRpBdoMrqxq7dq1myMdQ33W6bxSSoWd22utPutgle/xeSIciWqt337+2z7bire16ew6/bv0r5x91uxm17OeOHFiztixY8vXrFmTPHjw4Mqf/OQnB2fNmtW7qKjI8fLLL+8YMmRIzQ9/+MPs3bt3xyUkJPjmzZu36/TTT69qbI1qm83GxIkT+5eUlNg9Ho/MnDkz/0c/+tFhgMTExFGVlZVrAe67777ub731VrqIcP7555c8++yzexvGtnPnTufEiRNzgo+3bt2akJubu37AgAFt8k3UOp9apVTYeXwey0yIAmC32QFN1ip0e/bsiZ8/f/6OU089ddfw4cMHv/baa+mrVq3a/Prrr3d++OGHe/bu3bt2xIgRlZ988sn29957L+XGG2/st3nz5k1NrVG9aNGibWlpab59+/Y5Tj/99EHXX3/9YZvtSKPzW2+91WnRokVdVq9evTklJcW3f/9+e2NxZWdnu4PXYv/+97/PWLp0aUpbJWrQZK1UTLHaADOH+P8EeU27TK+s2lBzFXA4BRbiqAIYMGBA1XnnnVdqs9kYPXp05UMPPdRr7969cX//+9+3AVx++eVlt956q6OoqMje2BrVNTU1cuedd2YuW7Ys2WazceDAAVdeXp6jb9++dd8eP/74404/+tGPDgbXwu7evXuzv6wfffRR0iuvvJKxbNmyNm1G1z5rpWKI5ZJ1oEleK2sVKpfLVbfKlc1mIz4+3gDY7Xa8Xq8Yc+wiWCJiGlujeu7cuWlFRUWO9evX527evHlTenq6u6qq6qi8aIxpci3shnbt2uWcMmVK9vz587enpqb6Wn5F6DRZKxVDPD6Ppfqs65rBjSZr1TbOOOOMspdeeikdYOHChSldunTxpKWl+Rpbo7qkpMTetWtXd1xcnHn//fdT8vPzXQ2PN3HixNJXX321a1lZmQ2gqWbwmpoaufrqq0+aPXv23uHDh9e09fsK66dWRCYCTwF24AVjzB8abO8L/AXoHNjnXmPMB+GMSalYZrVJUezi/7vn9WkzuGobjzzySP7111+fPWDAgCEJCQm+l19++VuARx99tFvDNaoPHz5snzRpUv+hQ4cOPuWUUyr79etX3fB4kydPLl2zZk3iyJEjBzudTnPBBReU/OlPfzpmgNknn3yStGHDhqSHHnqo10MPPdQLYPHixVuzs7PbZBKBRpsM2uTAInbgG+BCIA9YCVxnjNlUb595wFpjzHMiMgT4wBiT3dxxx4wZY1atWhWWmJXq6C58+0LO6HkGs8+aHelQANh4cCPXLrqWZ85/hu9mfjfS4XRoIrLaGDPmRI7x9ddf7xwxYsTBtopJHevrr7/uOmLEiOyGz4ezGXwssM0Ys8MYUwu8CVzRYB8DdArcTwXyUUqFjcfnqatmrSDYDO726QxmSjUnnO1hvYH6owXzgNMb7PMA8JGI/AJIAi4IYzxKxTyvz6vN4EqdoBtuuKHvypUrk+s/9z//8z/7p02bVhSuc4bzU9vY8LmGbe7XAS8bYx4XkTOBV0VkqDHmqFF0InIrcCtA3759wxKsUrHAqgPM9NItFU1effXV3e19znA2g+cBfeo9zuTYZu6fAm8BGGO+BOKBrg0PZIyZZ4wZY4wZk5GREaZwler4PMZTd22zFThFZzBTKhThTNYrgRwR6SciLuBa4L0G++wGzgcQkcH4k3VhGGNSKqZ5fJ66atYKdAYzpUITtmRtjPEAU4F/ArnAW8aYjSIyS0QuD+w2A/iZiHwNvAHcZMI1PF0pZb1mcNFmcKVCEdZPbeCa6Q8aPDez3v1NwFnhjEEp5eczPgzGUs3gOoOZUqHRGcyUihHBhGilyjoYi1bWKlSzZ8/ulpOTc0r//v1PmTVrVjfwzyo2bty4nKysrKHjxo3LKSwsbHVfj9XWr25Ik7VSMSKYrC3VZy3aZ61Ct3LlyvhXXnklY82aNbm5ubkbFy9e3Hn9+vVx999/f89zzjmnbNeuXRvOOeecspkzZ/YI9Zgej/93r73Wr3a7j29OAet8xVZKhVVw/m0rNYPrALPolf/r3/Sp2bq1TdezjsvJqez1u4ebXM1r/fr1CaNHjy4ProB11llnlc2fP7/z4sWLOy9ZsmQLwJQpU4rOPvvsgcDe6dOn99qxY0dcQUGBc9++fa477rijYMaMGQcXLlyYMnv27J7dunVzb9q0KXH79u0bg+tXL1y4MOXBBx/slZGR4d60aVPixRdfXDxs2LCqZ599tntNTY28++6720855ZSa/Px8x80335y1d+9eF8CcOXN2T5gwoaKxuKdPn95r3759zt27d7vS0tI877///ret/dlY51OrlAqr4MQjVqqstRlctcbIkSOrZs2a1bugoMCelJRkPv7449QRI0ZUFBUVObKystwAWVlZ7kOHDtXlttzc3ITVq1fnlpWV2UeNGjXke9/7XgnAunXrktauXbtx0KBBx6w5vXnz5oS33357R7du3TxZWVnD4uLiDq5fvz539uzZ3R5//PFuL7744p4pU6b0mT59+v6LLrqofOvWra6LLrooZ8eOHRubin3dunWJy5cv35ycnHxcg6g1WSsVI4LVq5WWyNQZzKJXcxVwuIwePbp62rRpBeedd96AxMRE35AhQyodjubT2KRJkw4nJyeb5ORkz5lnnlm6dOnSpC5duniHDx9e0ViiBhg2bFhFMPn37du3ZtKkSSUAI0aMqFqyZEkKwOeff95p69atCcHXlJeX24uLi21dunRpdGnMiRMnHj7eRA2arJWKGcHq1UoDzGxiwyY2nRtcheyuu+46eNdddx0EmDp1au/MzMza9PR0z65du5xZWVnuXbt2OdPS0ur6VRquRR18nJiY2OR603FxcY2umW2z2fB6vQL+da5XrVqVG2oCTkpKOqH1rXWAmVIxIpgQrbSQB/j70LUZXIVq7969DoCtW7e6Fi1a1PmnP/3poYsuuujw3Llz0wHmzp2bPnHixMPB/T/88MPOlZWVUlBQYF+2bFnK+PHjG+1Xbq3x48eXPvLII92Cj7/44ouE5vY/Udb5iq2UCqtgU7OVKmvw96FrM7gK1eWXX37y4cOHHQ6Hwzz55JO7MzIyvA8++OC+q6666uSsrKyuvXr1ql2wYMH24P6jRo2qOP/883Py8/Ndd999977s7Gz3hg0b4k80jnnz5u255ZZb+g4YMGCI1+uV008/vWzcuHFhmzPcWp9apVTYWPHSLdDKWrXO6tWrtzR8rkePHt4vv/zym8b2z8nJqX7jjTd21X/u0ksvLbv00kvL6j9XWVm5trFtK1as2NLY63r27OlZtGjRjlBinjNnzgkv/6zN4ErFiGBCDC6eYRUOm0P7rJVqgVbWSsUIq1bWdptdK2sVFm1R0bbGU089lf7cc891r//caaedVt4WS2pqslYqRtRNimK1PmvRPuso4vP5fGKz2XTBpUZMmzataNq0aUXH+3qfzydAo6PGtRlcqRhRV1lbbTS4zaEzmEWPDYWFhamBpKLakM/nk8LCwlRgQ2PbrfUVWykVNpYdDS7aDB4tPB7PLQUFBS8UFBQMRYu9tuYDNng8nlsa22itT61SKmysuOoWaJ91NDn11FMPAJdHOo5YpN+MlIoRVlzIA7TPWqlQaLJWKkZYdjS42Ou+SCilGqfJWqkYYcW5wcH/5cFnTmjaZKU6PE3WSsWIuj5rizWDO8ShzeBKtUCTtVIxwsoDzLQZXKnmabJWKkZYtc/aJjatrJVqgSZrpWKEVUeDO8ShfdZKtUCTtVIxIli9Wq2y1mZwpVqmyVqpGBFsBnfarLXqll5nrVTLQmoPE5EeQN/6+xtjvghXUEqpthe8dMtqc4PrdKNKtazFZC0ivwN+BGwGgp8oA1wcxriUUm0suGa0FUeDa7JWqnmhfGq/BwwwxlSHOxilVPhYts9am8GValEofdbfhrifUsrCrDoaXCtrpVoWyqe2DFgrIp8ANcEnjTHTwxaVUqrNeX1eHOJAxFpLEdvFrutZK9WCUJL14sA/pVQU8/g8lmsCBx1gplQoWkzWxpg/i4gD6B94apsxelGkUtHGYzyWG1wG/gFvPp9OiqJUc1rsixaR7wDbgD851yqrAAAgAElEQVQDLwLfiMhZoRxcRCaKyBYR2SYi9zaxzzUisklENorI660JXikVOo/PY7nLtkCXyFQqFKF8zX4CuNgYswlARAYDrwJjmnuRiNiBZ4ALgTxgpYi8FzxOYJ8c4P8BZxljikWk2/G9DaVUS7w+ryUra5vYtBlcqRaEMsrbVT/BGmNyAVcIrxuLv8l8hzGmFngTuKLBPj8DnjHGFAeOfSC0sJVSreUxHsuNBAd/M7heuqVU80JJ1mtEZK6IjA/8ew5YG8LregN76j3OCzxX3wBggIh8LiLLRGRiYwcSkVtFZJWIrCosLAzh1Eqphjw+a/ZZ6wAzpVoWSrL+ObAd+BVwD7ADmBLC6xq7PsQ0eOwAcoBzgOuAF0Sk8zEvMmaeMWaMMWZMRkZGCKdWSjVk2dHgNp0URamWhDIavBp4NPCvNfKAPvUeZwL5jeyzzBjjBr4VkS34k/fKVp5LKdUCr7Fmn7VW1kq1rMnKWkTeCNyuFZE1Df+FcOyVQI6I9BMRF3At8F6DfRYA5wbO0xV/s/iO43kjSqnmWXY0eGAGM2MaNrwppYKa+5r9y8Dt5OM5sDHGIyJTgX8CduBFY8xGEZkFrDLGvBfYNkFENuFfJOSXxpii4zmfUqp5Xp/XcstjwpFVwHzGZ8kvE0pZQZPJ2hiTF7ibD1QbY4yInAwMBD4K5eDGmA+ADxo8N7PefQNMD/xTSoWR27gtmQyDTfNe48WO9eJTygpCGWC2FEgQkZ7AEuB/8E+OopSKIla9zjr4BULnB1eqaaEka5sxphL/Upl/MsZcBgwPb1hKqbZm1dHgNvH/GdJBZko1LaRkLSKnAdcDCwPPWe8Tr5RqllVHgwdj8hmdH1yppoSSrKcDDwKLjDEbROQk/E3jSqko4vFZcwYzbQZXqmWhXGf9b+DfAOJfCHe/Mea2cAemlGpblp3BLNA0r83gSjUtlFW3XhGRTiKSCGzEP3mJjt5WKsp4jEWvsw7EpLOYKdW0UJrBhxljSoEr8V+ylQncFM6glFJtz7KVdbAZXJfJVKpJIa26JSIO/CtmLQisoKUjQZSKMl6f15KjwYMx6QAzpZoWytfsF4DdwAZgiYj0BcrDGpVSEVC7Zw/lS5dii4un0yUXY4uPP2p7xbLlSJyLxFGjIhThifH4PJacwSw46E2bwZVqWigDzJ4Angg+FpE9wHnhDEqp9lD2r39R9umneIsPY9y1VK1aja+yEoCa7dvp/iv/jLslixbhysxk9003ATBgxXLsnTpFKuzjZtk+a5s2gyvVkiaTtYhcZ4x5Q0TuaGKXp8MUk1JhU/7fz3Gkp+EtLSNv6i8QpxNH9+7YkpJIvuB80m68kbypv8Cd718gzltaSv6Mu486RtHzL9BtRvSNsbRqn3XdpChaWSvVpOY+uV0Ct7qAtIp6xhgOvfwXDjzySN1zrqwsst+ajz019ah97ampmNpaAEo/XFz3fMqFF4IIxfPn03Xq7dji4ton+DZi2UlRRCdFUaolzS3k8Wzg9rftF45Sbatqw0Zqtmyh4vP/UvrBhySMHk3CyJGI3U76LT89JlEDiMuFcbspXbyYgvvvJ37IEHrM/C3xw4dT8cWXlH30ERWff0HKeedG4B0dPysvkQnaDK5Uc1r8mh0YUDYVyK6/vzHm6vCFpdTxq960CW9ZOQlDT2Hn5MAKrzYbGdPuIH3KFMTW/EUQ4nRiamspfPqPAPR48EEShg0FIC6nPwCewsLwvYEwsfoSmdoMrlTTQmkTew94BfgYvWRLWZi3vILCOY9T/Pobx2zrOXs2nb8X2vdLcTlx783HvXs33e69py5RA9iSkgDwlUffBRFWXcijLlnrDGZKNSmUZF1rjJkT9kiUOgG+mhr23HILVV99Redrf4ArK5tDL76IKzubjGl3kHDqqSEfS5zOugFmzl69jtpmS0z0n6+iou2CbwfGGDzGmgPMdLpRpVoWyif3jyJyH/BPoCb4pDFmXdiiUqoVfJWV7J0+g6qvvqL3nMfpdPHFAKTd+GMwBrG3rpoUlws8/v5TR9euR2+z2bAlJkZdsg4mQkv2WWszuFItCiVZDwBuASZxpBncAN8NV1BKtcRXWYmnsBBnZiZ7p8+g/LPP6HH/zLpEDbTYN90UcR7p13Wkpx+z3ZacTMXKFVRv3kz8oEHHdY72FkzWVqysgzFpZa1U00L55F4DZBtjalrcU6kwqt2zh+I33qTqq6+o+vpr8B754979t/fR5brr2uQ8Nper7r4j49grF21JSdRsyuXbK69i8ObcNjlnuAWXn7TiEpnB66x1iUylmhZK6bEOSAl3IEo1xXi97H/kUXZccinFr76Kqa31J+ZABdzpssvocv31bXa+2l276+4H+6iPItJm52ovdcnagpW1DjBTqmWhfHLTgc0ispyj+6z10i0Vdu59+9h2rn92206XXEK3X/0SZ/fuAFT897/U7txJ6hVXIG2YQH0VzY/09paVttm52kswWVtxNLg2gyvVslCS9cNhj0KpBiqWLad08YdUfrkMgC433ECP3/z6qH0SRo2idudO4nJy2vTc3jJ/ss7+21uNbu/+q3vI/+UvsTfSn21VOsBMqegWSrL+Aqg2xhgRORkYiH9da6XanDGGg88+y8E//glbcjJxgwbS645fkHrJJcfs22Pmb+ly3bU4u3dr0xh8ZWUAOAIVfEOpl11K5ZrVlC3+Z5ueN5yClbWlJ0XRylqpJoWSrJcC3xWRVGAJsBa4FvhxOANTscedn8/+xx6j7MPFpF5xBT1mPdjs/Nu2hAQShg9v8zjihw2jctkyHF26NLmPLSERX1VVm587XCzdZ23TylqploTyybUZYypF5CfAn4wxfxCRr8IdmOq4vOUVHJ4/H+N2k3rVVeTddhvVGzfWbe/2y7tJ+8lP2rQfujUy//g0tbt3H3UJV0O2hARMdTXG6231ddyRYOlkrZW1Ui0KKVmLyGnA9cCtgees/9dJWZKvspKdP/gBtdu3A1D45JN12yQxkcynniT5O9+JVHgA2FNSSDjllGb3sSUmAOCrqsaenNQeYZ2Q4CIZluyz1spaqRaFkqynAw8Ci4wxG0TkJPxN40q12sG586jdvp0+8+ZSvXEjhU/5l0U/6YMPiDupX4SjC50k+JO1qaqEaEjWFq6sg9d+66pbSjWtxU+uMebfwL/rPd4B3BbOoFTHVPL++xTNm0fqFVeQ/N3vknTGGfiqa0gef1ZUJWrw91kDUdNvbeVkbQvMNKeVtVJNC2WJzP74q+tsjl4ic0L4wlIdTfmSJeTf+/9IHDuWHg8+APjn4O52152RDew41S3oocn6hAUra5/RRf2Uakoon9y3gT8DfwX0q69qtZL332ffr39D3MABZD7zDLb4+EiHdMLq+qwrKyMcSWjq5ga34HSjwT5rbQZXqmmhfHJ9xpg/hj0S1SFVrllD/j33kjhmDL2ffCIqBmOFwlbXZx0dlbXb5wasWVnrpChKtSyUucH/ISK3ikiGiHQK/gvl4CIyUUS2iMg2Ebm3mf0mi4gRkTEhR66iQtHceTjS0+nz3LM40tIiHU6bCQ4w02bwE6eXbinVslA+ubcEbn9b7zkD9G3uRSJiB54BLgTygJUi8p4xZlOD/VKAO4DloQatooOvqoqKL7+ky3XXYkvqGBV1UN0As2hpBvdZd4lMEcEmNk3WSjUjlNHgfY7z2GOBbYHR44jIm8AVwKYG+80GHgXuPs7zKItyFxRgamuJHzIk0qG0uSN91lFSWVv4Omvwx6XN4Eo1LaSv2SIyCBgC1I0MMsa83sLLegN76j3OA05vcNxRQB9jzEIRaTJZi8itBCZk6du32YJeWYivvAIAW0rHW2FVAtOgmproWObdynODg7/i18paqaaFcunWfcAEYBDwT+Ai4L9AS8m6sbkiTb3j2oAngJtaisEYMw+YBzBmzBjTwu7KInzl/gUxbMnJEY6k7QWnIjWe6BjBbOU+awCb2OpiVEodK5QBZj8AzgX2GWNuAEYQWkWeB9RvQs8E8us9TgGGAv8RkZ3AGcB7Osis4/CW+5eatHfEZO3wfwQ0WbcNu9i1slaqGaEk6ypjjBfwBAaDFQAnhfC6lUCOiPQTERf+lbreC240xpQYY7oaY7KNMdnAMuByY8yqVr8LZUkduhm8Llm7IxxJaKzeZ+2wOXRSFKWaEUqyXisinYEXgVXACmBNSy8yxniAqfibznOBt4wxG0VklohcfgIxqwjaXljOWyv3tLwjR9aF7ojN4ARX2tLKuk3Yxa7N4Eo1o9lPrvjXKHzAGHMYeEZE/gl0Msa0mKwBjDEfAB80eG5mE/ueE1LEqt0VlfsHUaUlufjV2+tYu7uY74/JbHEJS19FoBm8g122Bf7LjcTpxLijpLK2eLLWS7eUal6zn1xjjBGRhcCpgcfb2iUqZQlur495n+1gzsffIMCFQ7qzelcxANVuHwmu5ptUveXlSFwc4nK1Q7QR4HRi3NFRDVr5OmsIjAbXS7eUalIon9wVIjI61GpaRS+fz7BkayF5xVUcrqjlr8t3sb+0homn9GBbYTkfbiio27fK7W0xWfvKyjtmE3iAOBzRM8DMWLuy1gFmSjWvyU+uiDgC/c7jgZ+JyHagAv8lWcYYM7qdYlTtoLTaza2vrGLZjkN1z52ckcTM60/h4mE9qPH4ePmLnWwpKOPdtXuprPWQltR8xewrL8fWQeYCb4w/WUdHM3jd3OAWXMgD/It5aLJWqmnNfXJXAKOBK9spFhUhG/aWMOOtr9leWM7DVw3lzJPS6dU5gXjnkco53mnn52efzPtf5/Pu2r1U1bb8h9VbUY49ueONBA+Kxj7r4ApXVqMzmCnVvOaStQAYY7a3Uyyqnb21ag8vf76TTftKSU1w8vLNYxmf07XZ1yQGmr4rQ0jWsdAMHk2jwV02644dsItdl8hUqhnNJesMEZne1EZjzJwwxKPCyOsz/OWLnfz5v99itwm7D1UyuGcnHrhsCFeNyiQ1seWpKBOcrUjW5eU4+x7v1PLWJw5H1Awwc/vcOO3WnGoU/BW/XmetVNOaS9Z2IJnGpw1VUWR/aTWL1u1jwVd7WZdXwui+nXHabXxvdCa3nXsyTnsol9v7BQeVVbtDS9b2pI5bWeOMngFmbq/bsvOCg78vXZvBlWpac8l6nzFmVrtFotrc+rwSXvlyJx9uKKC8xkN2eiJPXTuSy0f0avEa6aYkuvy/MqFU1t7yjt4M7oyeZO2zdrK227QZXKnmtNhnraJLrceHwfDql7t4dPEWar0+JgzpzowJAxnY48QHewX7rFfvKuaS4T2b3ddXWdnh1rGuL5oGmLl9bstetgWBSVG0slaqSc19es9vtyjUCavxePnTv7fxf0u24/b6Fyb7Tk5XHr9mBN1S4lt4deiCzeAvfv4tMy9rep1q4/OBx9NxJ0Qh+i7dsnJl7RAHtb7aSIehlGU1mayNMYea2qasY19JFWXVHu588ys27Svl/EHdGNW3M33SEk+oubsp6S1cWx1kav1/eDt6siZKBph5fB5LJ2u7zY7Xq5W1Uk2xbruYataeQ5U89a+tvL06D4CUeAd/un4UFw/tic0Wvh4MEeGmcdm8u3Zvs/sFm4eD6z53ROJ04KusinQYIbH8aHC9zlqpZmmyjiLGGP6Ve4AVOw/xzpq9HK6s5eazsumZGs+FQ3rQr2v79A877YLb2/xlNkcqa+smiBMWRdONWr0ZXKcbVap5mqyjgNdnKKly89CiTbyzZi8uh43hvVOZ/dOxDO7Zqd3jcTls1HpCTdYduBnc6YqaAWYer/WbwXWJTKWapsna4nw+w40vruC/2w5itwm3n3sy084fgMsR+rXRbc1lt+PxGXw+02STe0w0g2tl3WbsopOiKNUcTdYWZIxhyTeF/N+S7azPK6Gi1ktOt2Se+MFIhvZOjXR4OB3+BF3r9RHfxFzTwcra1pEr6ygbDZ7gTIh0GE3ShTyUap4ma4updnu54421fLRpP12TXVwyvCej+3bhB6f1afOR3cfLFZjxrNbrO2qxj/qCyZoOXllHy2hwq1fWDnFoM7hSzdBkbSFbCsr41d/XsS7vML++eBA3jssmzmG9VZKCTfDuZvqtg83gHbqydkXRpCgWn27UJjatrJVqhiZrC/h6z2H+umwXH6zfR5zTzlPXjuLyEb0iHVaT6lfWTYmFAWY6GrztOGwOfD7ts1aqKZqsI8jt9XH/ext5ffluAC46pTszLzuF3p2t27cI1C384faYJveJjQFm0TU3uJWnG9UlMpVqnnU/vR1cUXkNM/72Nf/ZUsg1YzL5yfh+DOrR/pdhHY9gM3htMzNO+WKgstbR4G1HB5gp1TxN1u3MGMPf1+zl0cWbOVzp5r5LBnPLd06KdFitEqysa5urrDVZW4rlk7XOYKZUszRZt6MDZdX8+p31fJJ7gOGZqcy5ZiTjc7pGOqxWi3OE0GcdC83gTie43RhjLDNSvylubxRMN6qVtVJN0mTdDnw+w2vLd/G/H31DtdvLby8dws3jssM6h3c41fVZNzvALJCsO3Jl7Qx8fDwey1+iFhULeWhlrVSTNFmHmddnuOHPy/liexHjTk5n1hVD6d8tOdJhnZC6PuvmLt2KgWZwHP6Pj/F4LN2CYIyJjmZwrayVapIm6zCqrPVw2R//y/bCCu68IIc7zsuJ2mq6Pqf9yAxmTYmJZnCH/71Zvd/aa7wYjLWTdWCAWTR0KSgVCZqsw+iZT7exvbCC31w8mJ+O79chEjVoZR0k9SprK3P7/F+crN5nDf4vFg7RP0tKNaSfijCodnv59bvreWfNXi4d3pOffTe6Rnu3pG5SlFhP1oFWg2D/vFXVJWsLV9bBa8B1MQ+lGqfJuo0VltXw/95Zxye5B7hpXDZ3XpAT6ZDaXN10o6E0gzs67q/YkQFmFk/WXusn62Bl7fF5cNk77hc8pY5Xx/1LGgGfbzvI7a+vidrrp0MVcjO404nYIreUZ7hFXTO4hZO1Tfy/JzrITKnGabJuIxv2lnDjiyvo3SWBP984hlOz0iIdUtiEdulWLbYOPLgMOGo0uJVFQ591sBk8Zi7f8nmhugQSO+7fCdW2wlr2iMhEEdkiIttE5N5Gtk8XkU0isk5E/iUiWeGMJ1y2HSjn5pdXkpbkYsFtZ3XoRA1HKuuaFlbd6sgjwaFen7XFV94KJmsrD9yqawaPhfnBl/0fzEqDR/vBYzngqYl0RCoKhC1Zi4gdeAaYBAwBrhORIQ12WwuMMcYMB94GHg1XPOGy7UAZ185bhjHw+s9Op0tSx+9vc9VV1s0t5FHboQeXQb1Ltyy+pnVdn7WFK+tgP3Uw1g7jm4/g8UGQt8r/uDQfFt9zZHvFAXj2TMj/Ch5IhSeGgq4+phoRzsp6LLDNGLPDGFMLvAlcUX8HY8ynxpjKwMNlQGYY42lzq3Ye4qpnvkAE3rz1DPp3S4l0SO3CGdJo8FiorIPN4NZOMB6f/8uElfusg8m61lcb4Uja0Oq/wOvfh7J98ML5cHgPzBl8ZPv9h+Hi/4VD22He2f7nSvbArC6RiVdZWjiTdW9gT73HeYHnmvJT4MPGNojIrSKySkRWFRYWtmGIx+9QRS33/H0dnRKcvPM/46J+VrLWsNsEu01aGA0eC5V1velGLSwaBpi5bP7flRpvlDYJf/Ar+OtkyFvtf7z9U1h419H7PDn0yP0LZ4MIjP1Z48cr2RueOFXUCmcnVmMzgDTabioiPwLGAGc3tt0YMw+YBzBmzJim217bweHKWl79chdvr8ljX0k1L910Gn3SEiMZUkS47LZmZzDz1cZOso6WAWZWXs86zh4HQK03Civrx/pDRaCI2PYxnHI1fLMYMgbCTxbD8+dB0bajX3PWHUfuP1ACH90HK54HT7X/uQO5kNpcbaNiTTgr6zygT73HmUB+w51E5ALgN8DlxhhLf63+YvtBLpjzGXM++YZO8U5eu+V0zuoffatmtQWnXZpvBtcBZpYRrFaDCdGK6prBoy1Zr/3rkUR91p2QcxHs+A9kj4cb3oX4VPjF6qNf84s1xx5nwkNw3/4jjyuLwhayik7h/Kq9EsgRkX7AXuBa4Pr6O4jIKGAuMNEYcyCMsZywPYcqufmllfRNS+QvPzmNU3qlRjqkiHI57M3PDR4DlTVRMjd4daBaS3AkRDiSpgWTteWbwT214K6AhC5QXgj/uN3//PTN0KlnaMeIb+Zvx9XPwzs/g8qDJx6r6lDClqyNMR4RmQr8E7ADLxpjNorILGCVMeY94DEgGfhbYPL+3caYy8MV0/HaX1rN1DfWIgKv/vR0eqTGRzqkiHO1VFnXujt8sq4bYGbx0eBVnioA4h3W/b21fDN4/ldHBoE1pqVEPfklePtm//3muiOGfs+frJfOgTNvb32cqsMKayeWMeYD4IMGz82sd/+CcJ6/Laz49hC3vbaGyloPT/5gpCbqAJfD1uJ0o7YE61ZybeFIn7W1m8Grvf7KOt5u3d/d4OA3S44GXz4XPvxV09tnFrd8jEGXHLnf3HSqNv/15lpZq4Y67lyQJ8gYw0uff8v1zy8jJd7BgtvPYuLQEJu5YoDTbmtxutEOX1lHyWjwYDN4NFTWlmsGN+bYRH3GbUfu/+gdCGVKXUe98QIWvt5dWZd1h4dGkNdnuPfv6/jb6jwuGNydOT8YQad4/YDV12JlHQvJOkoGmEVDn3UwWVtuUpSSvKMfP1Divz33N5C3Ek4+t/XHDHVUvs8X2hcBFRP0N6EBr88w462v+NvqPH5xXn/m3XCqJupGOO22Zqcb9VVXYYu3biXXJqLk0q0qTxU2sVn6Ouvg7GqWq6yXPn7k/oWzj9yPSz6+RA3+66ubc9Hv/LfVh4/v+KpD0sq6nlU7D/HI4s2s3FnM3RMGMPW8jre8ZVtpsbKurMKWaN1Kri0cqaytnaxrvDXE2+ORlpJEBFm2GTy40Mbd2yA5o53Ome6/rSrWhT5UHa2sA/69eT/fn/slu4oq+cPVwzRRtyDO0Xyfta+6GomPlWRtsabbBqo91Zbur4Z6zeA+i/0sg5V1UhvMpxBYrKRFCYEEXbKn+f1UTNHKOuD5z76lV2oCH931XZLi9MfSEmczM5gZnw9TXd3hm8HrkrXFm8GrvdWWHgkOR0aDW66yDmqLVolfrIai7S3vF6ymN7wDJ51z4udVHYJmJWBd3mG+3FHELy8aqIk6RC67Dben8ZlfTY3/D26HbwYP9lm7LXi5UT1VnirLV9Yigsvmsl6yjusEAye1zbHS+vn/taTXaP9tSo+2Oa/qELQZHHjuP9vpkujkx2dG5XLaEeF0NF1Z+6r8k3B0+GZwux3sdkytxZpuG4iGZnDwz2JmqdHgtZVQUwpdB7TveYMjwHcva9/zKkuL+WRd7fay5JtCLhnekxQd9R0yVzPXWZtAsrYlWD9BnChxuazfZx0FzeDgT9aWqqxLA0sZpEZo5d6dSyNzXmVJMd/mO3/lHiprvVw4RJucWsPlkKYr62r/db0dfQYz8PdbWz1ZV7or6RzfOdJhtCjOHmexZB1YprJTBFa/6nMGTSxSqGJUzFfWzy/dwdh+aXw3JzZXzzpeLnvTl275KmOjGRyCydrafdZltWV0cnWKdBgtSnImUemujHQYR5Tt89926tX+5+7UEyoPtf95lWXFdLI2xnCgrIaRfTpb+hpUK2puulFTHUPN4FFQWUdLsk50JlLuLo90GEcEm8EjMdArMV3nB1dHielkXVHrpdbjIy2pY0+LGQ7NTYoSU83gLpelB5gZYyirLSPFlRLpUFqU7Ey2XmUdlwqupPY/d2I6VB0Gr7UvC1TtJ6aTdXGFv/lSk3XrOe023F6Dz3dsv1pdM3gsJGuLV9ZVnio8xhMVyTrJmUSFuyLSYRxRVgAp3SNz7sSugNEpR1WdmE7WRYFkna7JutVcDv+vTmODzOqawTv4pCgQSNY1NRhf07O5RVJpbSlA1CRrSzWDVxyEpG6ROXdwYpQKbQpXfjGdrA9V+EeedtFk3Wouu/9Xp7GmcF+Vvxk8Virr8k8/ZcupYyIdSqPKasuA6EnWlmoGrzjQfvOBNxScH7yyKDLnV5YT08m6qFwr6+NVV1k3MsjMV3eddQwk68AyoMFry60mmKyjYoCZI5EKTwXGWOSSpfLCyFXWwbnINVmrgJhO1oe0z/q4Oesq62P/sMZaM7iVRVOyTnYl4zM+qjwW+OLjroaaEkiKdGWtzeDKL+aTtctuI1nnA2+15ivranA6LZ/I2oLVWw+iqs/a4R91XemxQFN4MElqM7iyiJhO1kUVtaQlufQa6+PQ3AAzX1VVTFTVAM4+EZqKMkSHa/yjiVNdqRGOpGVJgUukgq0BEVV+wH8bqWZwRxy4UnRiFFUnppP1oUCyVq3nsvu/4Gzdf+wfVlMdO8na1dfai7/sr9hPnD2O1DjrJ+u0OP8I6OLq4ghHQr0JUSJ06RaAMx72fR258ytLielkXVRRS3qyJuvjEe+0A/DwB7nHbPNVVSMdfHnMIFdW30iH0KyCygJ6JPWIitaj9AR/029RtQWafg/v8t92zo5cDBWFsOvzyJ1fWUpMJ+tDFTVaWR+nM07y/2HN7HJsUvY3g8dGsnb2tnYz+P6K/XRPjGB12ArBZH2wygKDqg7vBlfykeudI6nWQhPFqIiJ6WRdXOHWZH2c4p12xvZLo5EJzPAcOIAjPb39g4oAV2YEVmRqhWBlHQ26xHXBJjaLJOs90LkvWKFFYs0rkY5AWUDMJusaj5fyGo9eY30CEl12qmq9xzzvzsvDmWntirOtiMtF+q23gsN6VxR4fV4KKwujprK22+ykxadRVGWBZvCS3ZFbxzrohgX+W0dsjP9QzYvZZB2cECUtKS7CkUSvRJedytqjFxrwllfgLS62/CjptiQOB3g81pnMI6CgsgCv8UZNZQ2QkZDB/sr9kQ7DX1mn9prxLoIAAA0cSURBVIlsDF1z/LcL74xsHMoSYjZZF5b5pxrNSNFkfbwSnA6q3UdfuuXemweAK0YqawBxBqpqj7VWSNpwcAMAQ9KHRDiS0GV1ymJnyc7IBlGa719Ao+uAyMaRGBtdSSo0MZusDwSSdTdN1setscranedP1s7MCFcl7SnQBG4slqzXF67HZXMxsMvASIcSsn6p/dhbvpcab03kgji41X/bbVDkYgBw1huk6amNXBzKEmI4WfsXm+jWSZP18fIn6yN91mWffELe7VMBcFp84FVbEod/pjbLJeuD6xmUPginPXpmkjsp9SQMhm9Lvo1cEPlr/bfpOZGLIWj8dP/tpn9ENg4VcbGbrEtrEIGuyZqsj1eCy06Nx4c3MCT88Dvv1m2zd+4cqbDanViwsq50V7KxaCPDuw6PdCitMixjGACr96+OXBCf3O+/TbXAF86+Z/hv37klsnGoiIvdZF1WQ1qiq25BCtV6iS7/xChVbi/G56Nq7VqSzz+fkz74ICom4WgrwT5r43ZHOJIjPtn9CTXeGs7ve36kQ2mV3sm96ZvSl//s+U9kAgjOXNZzZGTO39CAi47cd1tggRMVMWHNVCIyUUS2iMg2Ebm3ke1xIjI/sH25iGSHM5769pVU6eCyE5To8iep0io3FUuX4i0uJuWCC4g7qV+EI2tfwcraKgPMSmpKeHL1kwzoMoDR3UdHOpxWu+SkS1i+b3ndALl29fFM/+3kF9v/3E3p7m9t4OHoGdWv2l7YkrWI2IFngEnAEOA6EWk4LPWnQLExpj/wBPBIuOKp70BZNV9uL+L0fhaYnSgKGWPwFBUxcvc6pn71d/Z972r2TPk5ju7dST77u5EOr91JYLa2sk8/jXAksKNkB9M+nUZxdTEPnfUQNom+lqMfDv4hPZJ6cMe/72Dxt4vxmWMXi2lz7ir47xOw/m+Q3h/STw7/OUN108Ij97d+HLk4VESFcyaHscA2Y8wOABF5E7gC2FRvnyuABwL33wb+JCJiwnDB6uIN+3ht+W52FlWQV+xvTvreqbFzedGJKPrzn6lcuQpvSQnekhLcBQWYykrswCSELelZdL35F4yfeiO2pKRIh9vuUs49h8Qzz2D/7Ico/+wzevzmN7j6tt+c4ZXuSh5e/jDbDm9jU9Em4uxxPDT+IQanD263GNpSalwqz57/LD//5Of88rNfMvOLmQxOG8ylJ1/K5JzJ4eliqa2ETx7w37dSVQ2Q0Nm/+lfFAdj1BeRcGOmIVASEM1n3BvbUe5wHnN7UPsYYj4iUAOnAUfMNisitwK0AfY/zj2BFjZfSKjej+3bhmlP78J0BGQzPjJ1BUCfCva8A94H92FNTieuWQ9L4s3Bl9iEupz/l/Yfw6PtbuOP8nJhM1AC2pCT6vvACRc8/T/Gb89v95xDviOfrwq/pkdiDO0ffyZX9r6ybZzta9e/Sn/eufI8nVj9BcU0xW4u3srJgJZNzJofnhEnpcOsS/7XVrsTwnONE3P2Nf7nMpOj+f1XHT8I165KIfB+4yBhzS+DxDcBYY8wv6u2zMbBPXuDx9sA+Tc43OGbMGLNq1aqwxKzUiTLGxNTguvbiMz68Pm9UXYZmNSKy2hgzJtJxqOMTzg6tPKD+zBiZQH5T+4iIA0gFdLV1FbU0UYeHTWyaqFVMC2eyXgnkiEg/EXEB1wLvNdjnPeDGwP3JwL/D0V+tlFJKRbOw9VkH+qCnAv8E7MCLxpiNIjILWGWMeQ/4M/CqiGzDX1FfG654lFJKqWgV1nX9jDEfAB80eG5mvfvVwPfDGYNSSikV7aLvIkyllFIqxmiyVkoppSxOk7VSSillcZqslVJKKYsL26Qo4SIihcCuMBw6FSgJw3GtINrem1XjjXRckTh/e52zKw1mLlRtLssYkxHpINTxibpkHS4iMs8Yc2uk4wiHaHtvVo030nFF4vztdU4RWaWzaynVNG0GP+L9SAcQRtH23qwab6TjisT5I/2elVJoZa2UsgCtrJVqnlbWSikrmBfpAJSyMq2slVJKKYvTyloppZSyOE3WSimllMVpsj4BIjJYRP5PRN4Wkf+JdDyxTP8vlFIdWYdJ1iJiF5G1IrLwBI7xoogcEJENjWybKCJbRGSbiNwLYIzJNcb8HLgG0JGsgIh0DiTMzSKSKyJnHudx9P8ixolIkoisFpFLIx2LUpHWYZI1MA3IbWyDiHQTkZQGz/VvZNeXgYmNvN4OPANMAoYA14nIkMC2y4H/Av86keA7kKeAxcaYQcAIGvyf6P9F7GrqC1hjX74C7gHeat8olbKmDpGsRSQTuAR4oYldzgb+ISLxgf1/BjzdcCdjzGfAoUZePxbYZozZYYypBd4Ergi85j1jzDjghyf8RqKcyP9v715j7agKKI7/l20opa1CxUZ56AVCkFdoBQMtairgFxUTAZUEDFgSTAioATU0gWAUjfh+hBhFEBINBZpSMGmkBCgNBUotYGlRPxBJekV5llqwULkuP+w5zTA59N7TXnvnnK7fpzP77Jm9p5PcNXvPPh29HfgIcD2A7W22X25Uy7XYc91I4wbsrW6+JJ0GPAk8u7s7GdFGkye6A+PkJ8DXgRndvrR9m6RDgEWSbgMWAB/r4fgHAhtr28PAiZLmA2cAU4BlO9HvQXMo8DzwG0nHAWuBL9t+tVMh12LPZXulpKFG8fabLwBJnZuv6cA0SoBvlbTM9n93Y3cjWqXvw7p6nvWc7bXVH+yubH+v+kPwC+Aw26/00kz3Q3oFsKKH4wy6ycAHgEtsr5b0U+By4Mp6pVyLqOl682X7YgBJ5wMvJKhjTzcI0+AnA5+S9DRlSvQUSb9tVpL0YeAY4Hbgqh7bGAYOrm0fBDyzU70dbMPAsO3V1fZiSni/Sa5F1HS9+dr+wb7R9k4vGo0YFH0f1rYX2j7I9hBwNnCv7XPrdSTNAa6jTK99AZgp6eoemlkDHC7pEEl7Ve3cOS4nMEBs/xPYKOmIquhUynPH7XItoiE3XxFj0PdhPUb7AJ+x/VQ1nXYeXd6JLelm4CHgCEnDki4AsP0GcDFwF2V18622N+y23veXS4DfSVoHzAa+0/g+1yLqcvMVMQb5v8EjYreobsDmA/tTVnlfZft6SR+nLBKdBNxg+9sT18uIdkpYR0REtNyeMg0eERHRtxLWERERLZewjoiIaLmEdURERMslrCMiIlouYR0REdFyCeuIiIiWS1hHRES0XMI6YkBI+rmkRyV9cKL7EhHjK2EdMQAkTQNmAV8EPjnB3YmIcZawjr4j6ceSvlLbvkvSr2vbP5R06Ti32cs7t8dyvH0lXVTbHpK0foz7TpV0v6RJnTLbrwLvobzT+2eS9pK0UlLfv7M+IhLW0Z8eBOYBSHob5cUQR9e+nwesmoB+9WJf4KJRa3W3AFhie6RTIOmdlDeabQFGbG8D7gE+t6sdjYiJl7COfrSKKqwpIb0e2CJpP0lTgCOBxyQtlbRW0gZJF3Z2lnRNY1T7DUmXSTpX0iOSHpf0y/rItVa3a51qZPxnSddV7S2XNLX67kpJf5F0t6SbJX0V+C5wWHWc71eHn9Rt/y7OAe5olF0B/ADYABxVlS2t6kZEn0tYR9+x/QzwhqT3UkL7IWA1MBc4AVhXjSwX2D6+KvtSNfoEWMSbR5yfBf5YlZ1sezYwQiPoJB05Sp3DgWttHw28DJwp6QTgTGAOcEbVF4DLgadsz7b9tbfav3nu1TufD7X9dK1sqPp3uIXyju/OLMN6IIvNIgZAnmdFv+qMrucBPwIOrD5vpkyTQwnoT1efD6aE4Yu2H5M0S9IBwLuATcCxwPHAGkkAU4HnGm2eOkqdv9l+vPq8FhiiTNHfYXsrgKTf7+Ccuu3ftD8lyOuuBr5p25K2h7XtEUnbJM2wvWUH7UZEyyWso191nlsfSxlBbgQuA/4F3CBpPnAaMNf2vyWtAPau7b8YOAt4N2WkLeAm2wt30OZodV6vfR6hhLl6OKdu+zdtpXYekmZTRuwfknRt9d0TtfpTgNd66ENEtFCmwaNfraL8ROkl2yO2X6Is2ppLmRZ/B7CpCur3Ayc19l8EnE0J7MWUxVhnSZoFIGmmpPc19hlLnaYHgNMl7S1pOvCJqnwLMKPXk7a9ifJsuxPY1wCn2x6yPQQcRzWyrqb9n7f9n17biYh2SVhHv3qCMiX8cKNss+0XgD8AkyWtA77VqIftDZSw/Lvtf9h+krJIa3m1z92Un0LV9xm1TpPtNcCdwJ+AJZRn45ttvwiskrS+tsBsrJZTRtKnANNs31Nr71lgmqSZwEeBZT0eOyJaSLYnug8RA03SdNuvSNoHWAlcaPvRXTjeHOBS258fpd4SYKHtv+5sWxHRDnlmHfH/9ytJR1GeJ9+0K0ENUC2Qu0/SpPpvreuqVeNLE9QRgyEj64iIiJbLM+uIiIiWS1hHRES0XMI6IiKi5RLWERERLZewjoiIaLmEdURERMslrCMiIlouYR0REdFy/wNGC4z2m5f8MAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for b in bands:\n",
    "    plt.plot(Table(data = parse_single_table(FILTERS_DIR + b + '.xml').array.data)['Wavelength']\n",
    "            ,Table(data = parse_single_table(FILTERS_DIR + b + '.xml').array.data)['Transmission']\n",
    "            , label=b)\n",
    "plt.xlabel('Wavelength ($\\AA$)')\n",
    "plt.ylabel('Transmission')\n",
    "plt.xscale('log')\n",
    "plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
    "plt.title('Passbands on {}'.format(FIELD))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### IV.a - Depth overview\n",
    "Then we plot the mean depths available across the area a given band is available"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ukidss_j: mean flux error: 7.015462398529053, 3sigma in AB mag (Aperture): 20.59205611010301\n",
      "90prime_g: mean flux error: 0.20077410340309143, 3sigma in AB mag (Aperture): 24.450427625315577\n",
      "90prime_r: mean flux error: 0.31159937381744385, 3sigma in AB mag (Aperture): 23.973205422566515\n",
      "mosaic_z: mean flux error: 1.0063319206237793, 3sigma in AB mag (Aperture): 22.70034374162011\n",
      "ukidss_j: mean flux error: 14.975190162658691, 3sigma in AB mag (Total): 19.768765998223763\n",
      "90prime_g: mean flux error: 2.4989206790924072, 3sigma in AB mag (Total): 21.7128156858491\n",
      "90prime_r: mean flux error: 1.897189736366272, 3sigma in AB mag (Total): 22.011919946975105\n",
      "mosaic_z: mean flux error: 1.6178725957870483, 3sigma in AB mag (Total): 22.184836066186342\n"
     ]
    }
   ],
   "source": [
    "average_depths = []\n",
    "for b in ap_bands:\n",
    "    \n",
    "    mean_err = np.nanmean(depths['ferr_ap_{}_mean'.format(b)])\n",
    "    print(\"{}: mean flux error: {}, 3sigma in AB mag (Aperture): {}\".format(b, mean_err, flux_to_mag(3.0*mean_err*1.e-6)[0]))\n",
    "    average_depths += [('ap_' + b, flux_to_mag(1.0*mean_err*1.e-6)[0], \n",
    "                        flux_to_mag(3.0*mean_err*1.e-6)[0], \n",
    "                        flux_to_mag(5.0*mean_err*1.e-6)[0])]\n",
    "    \n",
    "for b in tot_bands:\n",
    "    \n",
    "    mean_err = np.nanmean(depths['ferr_{}_mean'.format(b)])\n",
    "    print(\"{}: mean flux error: {}, 3sigma in AB mag (Total): {}\".format(b, mean_err, flux_to_mag(3.0*mean_err*1.e-6)[0]))\n",
    "    average_depths += [(b, flux_to_mag(1.0*mean_err*1.e-6)[0], \n",
    "                        flux_to_mag(3.0*mean_err*1.e-6)[0], \n",
    "                        flux_to_mag(5.0*mean_err*1.e-6)[0])]\n",
    "    \n",
    "average_depths = np.array(average_depths,  dtype=[('band', \"<U16\"), ('1s', float), ('3s', float), ('5s', float)])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def FWHM(X,Y):\n",
    "    \n",
    "    half_max = max(Y) / 2.\n",
    "    #find when function crosses line half_max (when sign of diff flips)\n",
    "    #take the 'derivative' of signum(half_max - Y[])\n",
    "    d = half_max - Y\n",
    "    #plot(X,d) #if you are interested\n",
    "    #find the left and right most indexes\n",
    "    low_end = X[np.where(d < 0)[0][0]]\n",
    "    high_end = X[np.where(d < 0)[0][-1]]\n",
    "    return low_end, high_end, (high_end - low_end)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "data = []\n",
    "\n",
    "for b in ap_bands:\n",
    "    data += [['ap_' + b, Table(data = parse_single_table(FILTERS_DIR + b +'.xml').array.data)]]\n",
    "    \n",
    "for b in tot_bands:\n",
    "    data += [[b, Table(data = parse_single_table(FILTERS_DIR + b +'.xml').array.data)]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ap_ukidss_j (11695.0, 13280.0, 1585.0)\n",
      "ap_90prime_g (4180.0, 5470.0, 1290.0)\n",
      "ap_90prime_r (5680.0, 7150.0, 1470.0)\n",
      "ap_mosaic_z (8552.0, 10018.0, 1466.0)\n",
      "ukidss_j (11695.0, 13280.0, 1585.0)\n",
      "90prime_g (4180.0, 5470.0, 1290.0)\n",
      "90prime_r (5680.0, 7150.0, 1470.0)\n",
      "mosaic_z (8552.0, 10018.0, 1466.0)\n"
     ]
    }
   ],
   "source": [
    "wav_range = []\n",
    "for dat in data:\n",
    "    print(dat[0], FWHM(np.array(dat[1]['Wavelength']), np.array(dat[1]['Transmission'])))\n",
    "    wav_range += [[dat[0], FWHM(np.array(dat[1]['Wavelength']), np.array(dat[1]['Transmission']))]]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Depths on SA13')"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfoAAAEgCAYAAABcsLEAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XlcVeXWB/DfOswIDiiIIoITKkqKGZbaq9dKscycskGtq5XT7ZapmVbXbpaVt7TylpWNr157s5xySE3vdUpvAzgyiho4geIAgiAyrPePc7AjMsuBw+b3/Xz6dM7ez97POoc+rf08+zl7iaqCiIiIjMlU0wEQERGR7TDRExERGRgTPRERkYEx0RMRERkYEz0REZGBMdETEREZGBM91Ukikigid9d0HEREtsZETzXOknSzRSRDRNJEZI+ITBSRKvnvU0S+EpHXq+JctiQiL4rI7yKSKSInRWR5MW2+EpE8EWleZHtnEdksIudE5IaHY4jIv0QkWUQuichhEXnSlp+FiOwHEz3Zi/tV1RNAAIC3ALwA4POaDan6iMjjAMYAuFtVPQB0B/DvIm3qARgOIB3AqCKnyAXwLYAnSujiTQCBqlofwGAAr4vIrVX3CYjIXjHRk11R1XRVXQvgIQCPi0hnABARFxF5R0SOi8gZEflYRNws+/paRsAvWka0iSIyyrJvPMxJcYZlpLzOqruuInJQRNJFZLmIuFqOaSIi6y2zCxdEZFdJswsi0lNEfrOc4zcR6Wm1b7uIvCYiuy2zFT+KSJMSPvptADar6lHL95CiqouLtBkOIA3AHACPF/ne4lX1cwDRJXyv0aqaU/jW8k+bEmIhIgNhoie7pKq/AjgJ4E7LpnkAggB0BdAWgB+A2VaH+AJoYtn+OIDFItLekiyXAfiHqnqo6v1Wx4wEEA6gFYBbAPzZsn2apW9vAE0BvAhzYryOiHgB2ABgIYDGABYA2CAija2aPQpgLAAfAM4AppfwkX8G8JiIPC8i3UXEoZg2jwP4PwDfAOggIt1KOFexRGSRiGQBiAOQDOCHihxPRLUTEz3Zs9MAvEREADwF4DlVvaCqGQDeAPBwkfZ/U9UcVd0BcwIeWcb5F6rqaVW9AGAdzBcRgHkavBmAAFXNVdVdWnxRiPsAJKjqUlXNU9X/gzmJWl9MfKmqh1U1G+ap9a7FnAeq+i8AfwUwAMAOAGdFZGbhfhFpCeBPAL5W1TMwT+s/Xty5SqKqkwF4wnzxtApATulHEJERMNGTPfMDcAHmkbU7gEjLdHoagE2W7YUuquplq/dJAK5bsFaMFKvXWQA8LK/fBnAEwI8icsw64RbR3NKPtSRL3GX1cQNVXaaqdwNoCGAigDkiMsCyewyAWFXdb3m/DMCjIuJU0vlK6CNfVX8C0ALApIocS0S1ExM92SURuQ3mhPkTgHMAsgF0UtWGln8aWBatFWpkWaxWqCXMMwJAMdPupVHVDFWdpqqtYR6dTxWRu4ppehrmxYPWWgI4VZH+iuk/V1W/A3AQQGfL5scAtBaRFBFJgfk2QRMAAyvZjSN4j56oTmCiJ7siIvVFZBDM96H/paqHVLUAwKcA3hURH0s7P6vRbqFXRcRZRO4EMAjAd5btZwC0rkAMg0SkreWWwSUA+ZZ/ivoBQJCIPCoijiLyEIBgAOvL/4mv9flnEblPRDxFxCQiAwF0AvCLiNwBc1IOg3nqvyvMFwBfwzJ9L2auMK8DgIi4ioiL5bWPiDwsIh4i4mD53h4B8J+KxklEtY9jTQdAZLFORPIAFACIgXnE+rHV/hdgXnz3s2Xl+ikAHwHYbNmfAuAizKPsLAATVTXOsu9zAN9Zpvy3q+qQMmJpB+ADmG8NXASwSFW3F22kquctFyXvW2I5AmCQqp6ryAe3uATzor9/AXCA+RbAJFX9SUQ+BvC9qh6yPkBE3gewy7IosD6A3612Z1vOEQjzjMYkmL9Pk2X7FFX9vhJxElEtI8WvMSKqPUSkL8yj/xY1HQsRkb3h1D0REZGBMdETEREZGKfuiYiIDIwjeiIiIgNjoiciIjKwWvHzuiZNmmhgYGBNh0FEVKtERkaeU1XvsluWeLyPo6PjZzA/t4EDQ/tUACAqLy/vyVtvvfVscQ1qRaIPDAxERERETYdBRFSriEjRRzRXiKOj42e+vr4dvb29L5pMJi7oskMFBQWSmpoanJKS8hnMJahvYLMrNBHxF5FtIhIrItEi8myR/dNFREsp20lERDWrs7e39yUmeftlMpnU29s7HX88LvsGthzR5wGYpqp7RcQT5oIkW1Q1RkT8AdwD4LgN+ycioptjYpK3f5a/UYkDd5uN6FU1WVX3Wl5nAIjFH1W93gUwAxUsNkJEREQVUy2LK0QkEEAozAU6BgM4paoHyjhmvIhEiEhEampqNURJRER1ycKFCxs/9thjLYtu/8c//uH9wQcfNC66PT4+3rldu3adqqLvhx56KCAyMtK1Ks5VFpsvxhMRDwArAUyBeTr/JQD9yzpOVRcDWAwA3bt358ifiIiqxYwZM2w+uly+fPlNLZSsCJsmehFxgjnJL1PVVSISAqAVgAPmCqBoAWCviISpaootYyEiosp7fsUB/8MpGe5Vec4gX8+st0d0OVFWu7vvvrtNcnKyc05OjmnixIlnpk+ffs7d3T101KhRqbt37/Zs0KBB/sqVK481b948r7jjw8LC2r/zzjsn/ud//icrOTnZsXv37h1PnTp1XTXIb775psFbb73VbOPGjUfefvttHw8Pj/w5c+ac2bVrl/uTTz4Z6ObmVtCjR4/MwvYRERGuY8eObZWbmysFBQVYuXLl0YCAgNzBgwe3Tk5Odi4oKJAZM2acfuqppy6WFVNFv7eKslmit9Ty/hxArKouAABLmU0fqzaJALpXsqxn2TbOBFIOld2urvINAQa+VdNREBGVatmyZYlNmzbNz8zMlNDQ0ODRo0dfzM7ONnXr1i3r008/PTl9+vRmM2fObL5kyZJKLfBesmRJw/fff7/pli1bEry9vfOt9z3xxBOB77777vH77rsvc8KECdcqZP7zn//0njx58plJkyZduHLliuTl5WHFihUNfH19c7dv334EAM6fP+9wc5+8athyRN8LwBgAh0Rkv2Xbi6r6gw37JCIiGyjPyNtW5s2b13TDhg0NASAlJcUpOjra1WQy4cknn7wAAOPGjTs/bNiwtpU59549ezwPHDjgvm3btsNeXl4F1vvOnz/vkJGR4XDfffdlFvbzn//8pwEA3HHHHZffeeedZidPnnR++OGHL4aEhOR069Yt+6WXXvKfNGmS3wMPPJAeHh6eWVyf1c1miV5VfwIgZbQJtFX/ADhaJSKq5davX++5Y8cOz4iIiDhPT8+CsLCw9tnZ2TcsJLfcDi6Wo6Oj5uebB+pZWVnXNWzZsmXO8ePHXaKiolyLTqOraonnnThx4oU777zz8urVqxsMHDgwaNGiRYmDBw/O2Lt3b8zKlSsbvPTSS35bt2699M477yRX4mNXKT7SkIiI7FZaWppDgwYN8j09PQv27dvneuDAgXoAUFBQgC+//LIRAHz11VeNw8LCMko6h7+/f86vv/5aDwCWLVvWyHpfixYtrq5cufLI2LFjW0VERFy3Cr5Jkyb5Hh4e+Zs3b/aw9ONVuC8mJsa5Y8eOOS+//PLZ/v37p+3fv98tMTHRydPTs2Dy5MkXpkyZcmb//v1VuqahsmrFI3CJiKhuGj58ePrixYu9g4KCgtu0aXOlS5culwHAzc2tIDo62q1Tp06+np6e+atWrTpW0jlmzpx55qGHHmr9zTffNL7zzjsvFd3fpUuXnCVLlhx76KGH2qxdu/aI9b7PP/88sXAxXr9+/a4du3TpUq/vvvuusaOjo3p7e+e++eabp3/66ad6s2bNamEymeDo6KiLFi0qdWV9abMQValW1KPv3r278ln3REQVIyKRqtq9sscfOHAgsUuXLrZZLH2T3N3dQ7OysvbVdByVFRQUFLx27dojHTp0uFoV5ztw4ECTLl26BBa3j1P3RERE1ahnz57t2rdvn11VSb4snLonIqJap7jR/JgxY1r+9ttvHtbbJk2adObZZ589X32RXe+ee+5pc+LECRfrbXPnzj05fPjwG24h2AoTPRERGcLSpUvtrlDali1bjtZ0DJy6JyIiMjAmeiIiIgNjoiciIjIwJnoiIiIDY6InIiJD++9//+vWtWvXDkFBQcH9+vVre+HChWu5b9asWb4tW7bsHBgY2HnlypX1K3ruKVOmNF+zZo1n1UZctZjoiYjI0J566qnAuXPnnjx8+HDM4MGDL7766qu+ABAZGem6atUqr/j4+OhNmzYdnjJlSsu8vGIr3RYrLy8P77333ukhQ4aU+Phde8Cf1xERUdnW/MUfZ2Oq9tntPsFZGPKhzevRJyYmug4cODATAAYNGnRpwIABQe+///7pFStWNBw2bNgFNzc37dChw9WAgICc7du31/P3988NDw9vFxoaejkqKsq9devWV7777rtET0/PAj8/v5BHHnnk3LZt2+pPmDDh7ObNmxsMGjQofezYsRf9/PxChg4deuGnn37yzMvLk48//jhp5syZfklJSS5//etfz8yYMSMVAP72t781Xb16tdfVq1flvvvuS3v33XdPl/TZn3/++WYrVqzwatas2dXGjRvnhYaGZs2ZM+dMRb5mjuiJiMiuLVu2LDE6Ojp2//79MZ988knTlJQUh8J69DExMbG9evXKmDlzZvOSjm/Xrl32119/3RAA/vWvf3mlpKQ4A8CpU6ec/f39rz2drnnz5ldPnDjhDJgvDiZOnJh6+PDhGE9Pz4K3337bu7Cdq6trQWRkZPz48eMvFu3L39//6v79++N69OiROW7cuMB169Yd/eWXX+Leeuut5gCwatWq+keOHHE9ePBgbGxsbMz+/fvdN27c6FH0PACwc+dO93Xr1jU6dOhQzIYNG44ePHiwXmW+P47oiYiobOUYedvKzdaj/+KLLxKffvpp/zfffLNZeHh4mpOTkwLmMrRFiYgCgK+v79X+/ftfBoAxY8acX7hwoQ+AMwDw2GOP3ZDgC40cOTINAEJCQrIuX75satSoUUGjRo0KXFxcCs6dO+ewadOm+jt37qwfHBwcDABZWVmmuLi4azMO1rZv3+4xcODANA8PDwWg99xzT1o5v7LrMNETEZHdqop69KGhoVd2796dAAAHDx50+fHHHxsC5hK1hSN4ADh9+rRzixYtcos7n/V7T0/PgpL6cnV1VQAwmUxwdna+diVhMpmQm5srqoopU6YkP//882UWC6qqonOcuiciIrtVFfXoT5065QgA+fn5eOWVV5o98cQTZwFg+PDhaatWrfLKzs6WuLg458TERNe+ffteBoDk5GTnrVu31gOAr7/+2qtnz543jLgrY+DAgZeWLl3aJD093QQAv//+u1NhfEX17ds3c/PmzQ2ysrIkPT3dtHXr1oaV6ZMjeiIisltVUY/+iy++8Pr88899AODee++9+Mwzz5wHgO7du18ZMmTIhaCgoE4ODg5YsGBBkqOjOS22bt36yhdffNF48uTJAa1atcqZPn16alV8nmHDhl2Kjo52ve222zoAgLu7e8GyZct+9/Pzu2EhYZ8+fbLCw8PTg4ODO/n5+eXccsstlxs0aJBf0T5Zj55KNe/XeYi7EFfTYVS7Dl4d8ELYCzUdBtFNYT36yomPj3ceNGhQu4SEhGhbnL8i0tPTTQ0aNCjIyMgw3XHHHe0//vjjpN69e2cVbVdaPXqO6ImIiOzU6NGjAxISEtxycnLk4YcfPl9cki8LEz2ViqNaIrJHtqxH3759+6vVOZpPSUlx6Nu3b/ui27dv3x6/bt2632/2/Ez0RERkCPZYj748fH198+Pi4mJsdX6uuiciIjIwJnoiIiIDY6InIiIyMCZ6IiIiA2OiJyIiQ2M9eiIiIgMzQj36goIC5OdX+KF4APjzOiIiKoe/7f6b/5GLR6q0Hn3bRm2zXuv1GuvRl1CPPj4+3nngwIHtevbsmREZGenx/fffHwkKCrpaXNvScERPRER2ra7Woy+MY+zYsedjY2NjKpPkAY7oiYioHMoz8raVulqPHgCaNWt29a677rpcnu+pJEz0RERkt+pyPXrAXN2uPO1Kw6l7IiKyW3W5Hn1VsdnJRcQfwBIAvgAKACxW1fdF5DUAD1i2nQXwZ1UtdiECERHVbXW5Hn1VsVk9ehFpBqCZqu4VEU8AkQCGADipqpcsbZ4BEKyqE0s7F+vRExFVHOvRV4491aMvr9Lq0dts6l5Vk1V1r+V1BoBYAH6FSd6iHgDbXGkQERFR9SzGE5FAAKEAfrG8nwvgMQDpAP5UwjHjAYwHgJYtW1ZHmEREVEvUlXr0vr6+lXtKjhWbTd1f60DEA8AOAHNVdVWRfbMAuKrqK6Wdg1P3REQVZ+Spe7pejUzdA4CIOAFYCWBZ0SRv8TWA4baMgYiIqC6zWaIX848OPwcQq6oLrLa3s2o2GECcrWIgIiKq62x5j74XgDEADonIfsu2FwE8ISLtYf55XRKAUlfcExERUeXZLNGr6k8AintU0Q+26pOIiIiuxyfjERERVRF7rE/PZ90TERFVkffee8/unvTKRE9ERGU6/eJL/jkJCVVaj96lXbus5m/MtXk9+rCwsPYhISFZBw4ccL9w4YLjl19++fvcuXObxcfHuz3wwAMXFi5ceBoA/v73vzddtmxZEwAYM2ZM6uzZs89eunTJNHjw4NbJycnOBQUFMmPGjNNPPfXUxenTpzfbtGlTw5ycHFP37t0zly1blmQymTB8+PDAwvr0O3bscJ8yZUrLrKwsk7Ozs+7cuTO+UaNGNxSpeeihhwIKn+F/5swZp3Hjxp2dP39+8s19u3/g1D0REdm1m61HDwDOzs4FERER8WPHjk198MEH23766afH4+LiopcvX94kJSXFYdeuXe5ff/1148jIyNiIiIjYJUuWeO/evdtt1apV9X19fXPj4+NjEhISoocNG3YJAJ5//vmzUVFRsQkJCdHZ2dmmb775poF1f1euXJFRo0a1ee+9947Hx8fH7NixI97Dw6PYSnTLly9PiouLi1m7du2Rhg0b5k2YMKFCD/gpC0f0RERUpvKMvG3lZuvRA8DQoUPTAKBLly7Zbdu2zQ4ICMgFAH9//5xjx445b9++3ePee+9Nq1+/fgEA3HfffRe3bdvmOXjw4PSXXnrJf9KkSX4PPPBAenh4eCYAbNy40XPBggW+V65cMaWlpTkGBwdnw/y0VwDAwYMHXX18fHL79OmTBQBeXl6llpvNysqS4cOHt3n33XePBwUFXa30l1UMjuiJiMhuWdejj4+Pj+nYsWN2RevRA9fXiXdxcbmuTnxeXl6JT4m95ZZbcvbu3RsTEhKS/dJLL/lNnz69WVZWlkybNi1g1apVRw8fPhwzevToc1euXLkuJlWFiJT70bNjxowJuP/++y8OGTKkxHK7lcVET0REdqsq6tGXR79+/TJ/+OGHhhkZGaZLly6Zfvjhh0Z/+tOfMhITE508PT0LJk+efGHKlCln9u/f756VlWUCAF9f37z09HTTunXrGhU9X5cuXa6cOXPGeceOHe4AcPHiRVNubm6xfb/55pvemZmZDm+88UbKzXyGknDqnoiI7FZV1KMvj969e2c9+uij57t169YRMC/G69WrV/bKlSvrz5o1q4XJZIKjo6MuWrQoqUmTJvmjRo1KDQ4O7tSiRYurhTFZc3V11WXLlh195plnWl65csXk6upasHPnzsMNGjS4YQr/gw8+8HVyctIOHToEA8C4ceNSZ8yYkXozn8eazYvaVAUWtSEiqjgjF7WxZT362qjGitoQERFRzeLUPRER1Tq2rEdvKytXrqz/0ksvtbDe5u/vn7Nly5ajtuyXiZ6IiAxh6dKlx2s6htIMHz780vDhw2Oqu19O3RMRERkYEz0REZGBMdETEREZGBM9ERGRgTHRExFRrRMWFtZ+586dN1TT69OnT9tz5845FN0+derU5rNnz256s/0mJiY6hYeHt77Z81QnrronIiLD2LFjxxFbnj8wMDB306ZNN/UUvurGRE9ERGX695JY/wunMqu0Hr2Xn0fWXY91LLUqXnx8vPOgQYPaJSQkRAPA7Nmzm2ZmZl4bsefn5+PBBx8MbNGixdWFCxee9vPzC4mIiIht1qxZ3gsvvOC7fPnyJs2bN7/auHHj3NDQ0CwAeP31132+/PJLbwcHBw0KCrqyfv36Yxs2bPCYNm1aS8BcIGfPnj1xxdWOLxpPbcBET0REtVJubq4MGTKkVXBwcPa8efOuKwiza9cu99WrV3sdOnQoJjc3F127dg0uTPQLFy70TUpKOuTm5qaF0/zz58/3XbhwYVL//v0vp6enm9zd3UstK1ubMNETEVGZyhp514TJkycHDBky5ELRJA8A27Zt87j33nvTPD09CwCgf//+aYX72rdvnz106NBWgwcPThs1alQaANx+++2Z06dP9x85cuSFRx555GKbNm0Mk+i5GI+IiOyWo6OjFhT8kXOt67537949c9euXfWzsrKKLUZfUo36bdu2JfzlL39JjYyMrNelS5fg3NxcvPHGGymfffZZUnZ2tqlnz54d9+3b51rVn6WmMNETEZHdatGiRd6FCxccU1JSHLKzs2Xz5s0NCvdNmDDhXP/+/dMHDRrUpmit9379+mVu2LChYWZmply8eNG0ZcuWhoD5nv7Ro0ed77///oxFixadzMjIcEhPT3eIjo52CQsLy547d25KSEjI5aioKMMkek7dExGR3XJxcdFp06Ylh4WFdWzRokVO27Ztr1jv//vf/37mueeecxg2bFirNWvW/F64vXfv3llDhw690Llz505+fn45YWFhmQCQl5cnjz76aKuMjAwHVZUJEyacadKkSf60adOa79mzp77JZNKgoKDsESNGpJcUk4jYf313K6xHT0RkUEauR19Tdu3a5T516lT/3377Lb6mY7HGevREREQ3aefOne5jxoxp/fTTT5+p6VgqglP3RERERfz6669ujz32WCvrbc7OzgWJiYlRNRVTZTHRExERFREWFpYdFxdX7bXjbYFT90RERAbGET2RgaS88QZyYuNqOowyuXTsAN8XX6zpMIjqBI7oiYiIDIwjeiID4SiZiIriiJ6IiOzaa6+95tOuXbtObdu27TRnzhwfADhz5oxDz5492wUEBHTu2bNnu9TU1Btq0JclNDS0Q9VHa39sluhFxF9EtolIrIhEi8izlu1vi0iciBwUkdUi0tBWMRARUe3222+/uS5ZssR77969sbGxsdGbNm1qeOjQIZdXXnmlWd++fTOSkpKi+vbtmzF79mzf8p4zLy8PALBv3z77X9BSBWw5dZ8HYJqq7hURTwCRIrIFwBYAs1Q1T0TmAZgF4AUbxkFERDdp80fv+Z87kVSl9eib+AdkDZg0pdSqeIcOHXLr1q1bZmEVul69emUsX7684aZNmxru2LEjHgAmTJhwvk+fPu0BnJo6dWrzY8eOuaSkpDglJyc7P/PMMynTpk07t379es/XXnutmY+PT25MTIz70aNHo93d3UOzsrL2rV+/3vPVV19t7u3tnRsTE+N+7733XgwJCcletGhR05ycHFm9evXRTp065Zw+fdpx7NixAadOnXIGgAULFhzv37//5eLiPn36tOOIESNapaWlOXbt2jVr+/bt9SMjI2ObNWuWV5XfYXnYbESvqsmqutfyOgNALAA/Vf1RVQs/6M8AWtgqBiIiqt26du2a/csvv3impKQ4ZGRkmLZs2dLgxIkTzufPn3cMCAjIBYCAgIDcCxcuXBu4xsbGum3dujXh559/jnv77bebJyYmOgHAwYMH67399tunjh49Gl20n7i4OLePPvroRGxsbPSKFSsaHz582PXQoUOxY8aMOTd//nwfAJgwYYL/1KlTz0RFRcWuXr366MSJEwNLinvmzJnN+/TpkxETExM7bNiwi8nJyc5V/uWUU7UsxhORQAChAH4psmscgOUlHDMewHgAaNmypQ2jIyKispQ18raVbt26XXn22WdT+vXrF+Tu7l4QHByc5ehYeuoaOHBgmoeHh3p4eOTdcccdl3bt2lWvUaNG+bfccsvlDh06XC3umJCQkMuFFw4tW7bMGThwYDoAdOnSJXvHjh2eALB79+76CQkJboXHZGZmOly8eNHUqFGjG2rX//rrrx5r1qw5AgAjRoy4VL9+/fxKfwk3yeaJXkQ8AKwEMEVVL1ltfwnm6f1lxR2nqosBLAbMRW1sHScREdmn55577txzzz13DgCefvppvxYtWlxt3LhxXlJSklNAQEBuUlKSk5eX17Up8aJ16Avfu7u735CQC7m4uFzLMyaTCa6urlr4Oj8/XwBAVRERERHr4eFRZk6yp4JxNl11LyJOMCf5Zaq6ymr74wAGARil9vRtEBGR3Tl16pQjACQkJDhv2LCh4RNPPHFhwIABaZ988kljAPjkk08ah4eHpxW237hxY8OsrCxJSUlx+Pnnnz179+5d7H30iurdu/elefPm+RS+37Nnj1tJbcPCwjKXLl3qBQCrVq2qf+nSpQr/KqCq2GxEL+ZLqM8BxKrqAqvt4TAvvuujqlm26p+IiIxh8ODBbdLS0hwdHR31vffeO+7t7Z3/6quvJg8dOrRNQEBAk+bNm19ds2bN0cL2oaGhl++66652p0+fdp4+fXpyYGBgblRUlOvNxrF48eITTz75ZMugoKDg/Px86dGjR0bPnj2PF9f2rbfeOj1ixIjWwcHBje64445Mb2/v3IYNG9bI9L3N6tGLSG8AuwAcAlA4XfIigIUAXACct2z7WVUnlnYu1qMnIqq4uliPfurUqc09PDzy58yZU6OlZLOzs8XR0VGdnJywdevWek8//XSALYvklFaP3mYjelX9CYAUs+sHW/VJRERkD44cOeI8cuTINgUFBXByctJPPvkksaZi4SNwiYjIMBYsWHC6Ovt7//33G3/00UdNrbfddtttmUuXLj0eGxtrF2Vuy5XoRcQFwHAAgdbHqOoc24RFRERk/5599tnzzz777PmyW9ac8o7ovweQDiASQI7twiHeUfAyAAAUqElEQVQiIqKqVN5E30JVw20aCREREVW58v6Ofo+IhNg0EiIiIqpypY7oReQQALW0Gysix2CeuhcAqqq32D5EIiIiqqyyRvSDANwPYCCAtgD6W94XbiciIrIp1qO/OaWO6FU1CQBEZKmqjrHeJyJLAYwp9kAiIjKUCysO++emXK7SMrVOvvWyvEYElVosx7oevaura0GfPn2Chg4dmv7hhx969+3bN+ONN95IePHFF31nz57t+9FHH50qT795eXlwdHSstnr0ubm5cHJyqo6uilXee/SdrN+IiAOAW6s+HCIioj9Y16N3cnK6rh79hAkTzgPmevQbN25sBJifjDdkyJBWt99+e1BAQEDn+fPnNwGA9evXe/bo0SPo/vvvb9W+fftOAODu7h5auO+2225rf++997YODAzsPHnyZL+PPvrIKyQkpGNQUFBwdHS0C2CuMT9gwIA2nTt37ti5c+eOP/74Y72S4p46dWrzRx55JKBXr17thg0b1srW31NpyrpHPwvmx9a6icgl/PGku6uwVJYjIvuxceNGpKSk1HQYZfL19cXAgQNrOgyqgLJG3rbStWvX7Dlz5vilpKQ41KtXT7ds2dKgS5cul8uqRx8ZGRmbkZHhEBoaGjx8+PB0wFyPft++fdHFlaqNi4tzW7FixTEfH5+8gICAEBcXl3OHDh2Kfe2113zmz5/v88UXX5worEc/YMCAzISEBOcBAwa0O3bs2A217QsdPHjQ/ZdffokrT7U7Wypr6v5NAG+KyJuqOquaYiIiIgJQe+vRA0B4eHhaTSd5oPy/o39RRIYB6A3zKvxdqrrGdmERUWVwlExGVBvr0QNAvXr1SuyvOpX3Hv2HACbCXIkuCsBEEfnQZlERERFZ1MZ69PakvCP6PgA6q6WmrYj8L8xJn4iIyKZqYz16e1KuevQisgrAc1Y/twsA8JaqPmLj+ACwHj0RUWWwHn3dURX16BsDiBWRXy3vbwPwXxFZCwCqOvimoyQiIqIqV95EP9umUZDdSlt3FFdPV8ntrVrFuXk9NLy/TU2HQUQVZE/16KszjtKUK9Gr6g7LdH07Vd0qIm4AHFU1w7bhERER2S/D1KMXkacAjAfgBaANgBYAPgZwl+1CI3vAUS0RUe1W3p/X/QVALwCXAEBVEwD4lHoEERER1bjyJvocVb32NCERcYT5wTlERERkx8qb6HeISOEz7+8B8B2AdbYLi4iIqHoZtWxteRP9TACpMD8kZwKAHwC8bKugiIiIqlt1la2tbuVddV8gImsArFHVVBvHREREdmbNmjX+Z8+erdJ69D4+PllDhgwptSpefHy8c3h4eLuwsLDMvXv3enTs2DFr3Lhx5+bMmeN3/vx5x6+++upYcHBwzqhRowKPHz/u4ubmVrB48eKkHj16ZG/YsMFj2rRpLQHz8+737NkTZzKZEB4e3jY9Pd0hLy9PZs+efXr06NFpgLlsbVZW1j4AePnll5t+++23jUUEd911V/qiRYtuqHWfmJjoFB4e3q7wfUJCgltsbOyhoKCgYgvn1JSyytQKgFcAPA1ziVoRkXwA/1TVOdUQHxER1XEnTpxwXb58+bFbb7016ZZbbum4bNmyxhEREXFff/11w7lz5zbz8/O72qVLl6ytW7ceXbt2refjjz/eKi4uLmb+/Pm+CxcuTOrfv//l9PR0U2FRmw0bNhzx8vIqSE5OduzRo0eHRx99NM1k+mOC+9tvv62/YcOGRpGRkXGenp4FZ86ccSgursDAwNy4uLgYAHjzzTe9d+3a5WlvSR4oe0Q/BebV9rep6u8AICKtAXwkIs+p6ru2DpCIiGpeWSNvW/Lz88sJCwvLBoCgoKDsfv36XTKZTOjWrVvW66+/3vzUqVMuK1euPAIAgwcPzhg/frzj+fPnHW6//fbM6dOn+48cOfLCI488crFNmzYFOTk5MmXKlBY///yzh8lkwtmzZ51Pnjzp2LJly2vV77Zs2VJ/9OjR5zw9PQsAoGnTpvmlxffjjz/WW7JkiffPP/9sl1P/Zd2jfwzAI4VJHgBU9RiA0ZZ9RERENuXs7FxsCVkHBwfk5+cXW7NFRPSNN95I+eyzz5Kys7NNPXv27Lhv3z7XTz75xOv8+fOOhw4dio2Li4tp3LhxbnZ29nW5UFVvKHVbkqSkJKcJEyYELl++/GiDBg3soixtUWUleidVvaGggeU+vZNtQiIiIiq/22+/PePLL79sDADr16/3bNSoUZ6Xl1dBdHS0S1hYWPbcuXNTQkJCLkdFRbmmp6c7NGnSJNfFxUXXrVvnefr0aeei5wsPD7+0dOnSJhkZGSYAKGnqPicnR4YNG9b6tddeO3XLLbfk2PZTVl5ZU/el3Wuwu/sQRERU98ybN+/0o48+GhgUFBTs5uZW8NVXX/0OAP/4xz989uzZU99kMmlQUFD2iBEj0tPS0hwGDhzYtnPnzh07deqU1apVqytFzzdixIhLe/fude/atWtHJycnvfvuu9M/+OCDGxbjbd26tV5UVFS9119/vfnrr7/eHAA2bdqUEBgYmGv7T11+pZaptSy8K66iiQBwVdVqGdVXtkzttq8W42zSMRtEZAw+Aa3xpz+Pr+kwiMhG6mKZ2rqq0mVqVbXY6QoiIiKqHcpbprZW4miViIiqwpgxY1r+9ttvHtbbJk2adMbeK9cBBk/0REREVcGe6stXVHkfgVthIuIvIttEJFZEokXkWcv2By3vC0Sk0veOiIjI5goKCgrK9zszqjGWv1GJP+2zWaIHkAdgmqp2BHA7gL+ISDCAKADDAOy0Yd9ERHTzolJTUxsw2duvgoICSU1NbQBzbi2WzabuVTUZQLLldYaIxALwU9UtAMr9MAIiIqoZeXl5T6akpHyWkpLSGbYdGFLlFQCIysvLe7KkBtVyj15EAgGEAvilOvojIqKbd+utt54FMLim46CbY/MrNBHxALASwBRVvVSB48aLSISIRKSmsmAeERFRZdg00YuIE8xJfpmqrqrIsaq6WFW7q2p3b29v2wRIRERkcLZcdS8APgcQq6oLbNUPERERlcyW9+h7ARgD4JCI7LdsexGAC4B/AvAGsEFE9qvqABvGQUREVGfZctX9TzA/E784q23VLxEREf2BP5cgIiIyMCZ6IiIiA2OiJyIiMjAmeiIiIgNjoiciIjIwJnoiIiIDY6InIiIyMCZ6IiIiA2OiJyIiMjAmeiIiIgNjoiciIjIwJnoiIiIDY6InIiIyMCZ6IiIiA2OiJyIiMjAmeiIiIgNjoiciIjIwJnoiIiIDY6InIiIyMCZ6IiIiA2OiJyIiMjAmeiIiIgNjoiciIjIwJnoiIiIDY6InIiIyMCZ6IiIiA2OiJyIiMjAmeiIiIgNjoiciIjIwJnoiIiIDY6InIiIyMCZ6IiIiA2OiJyIiMjAmeiIiIgNztNWJRcQfwBIAvgAKACxW1fdFxAvAcgCBABIBjFTVi7aKg4iool5dF42Y05dqOgwAQHDz+njl/k41HQbVYrYc0ecBmKaqHQHcDuAvIhIMYCaAf6tqOwD/trwnIiIiG7DZiF5VkwEkW15niEgsAD8ADwDoa2n2vwC2A3jBVnEQEVUUR9BkJNVyj15EAgGEAvgFQFPLRUDhxYBPCceMF5EIEYlITU2tjjCJiIgMx+aJXkQ8AKwEMEVVy33TS1UXq2p3Ve3u7e1tuwCJiIgMzKaJXkScYE7yy1R1lWXzGRFpZtnfDMBZW8ZARERUl9ks0YuIAPgcQKyqLrDatRbA45bXjwP43lYxEBER1XU2W4wHoBeAMQAOich+y7YXAbwF4FsReQLAcQAP2jAGIiKiOs2Wq+5/AiAl7L7LVv0SERHRH/hkPCIiIgNjoiciIjIwJnoiIiIDY6InIiIyMCZ6IiIiA2OiJyIiMjAmeiIiIgNjoiciIjIwJnoiIiIDY6InIiIyMCZ6IiIiA2OiJyIiMjAmeiIiIgNjoiciIjIwJnoiIiIDY6InIiIyMCZ6IiIiA2OiJyIiMjAmeiIiIgNjoiciIjIwJnoiIiIDY6InIiIyMCZ6IiIiA2OiJyIiMjAmeiIiIgNjoiciIjIwx5oOgIjI3uz69jDOncis6TAAAE38PXDnyKCaDoNqMY7oiYiIDIwjeiKiIjiCJiPhiJ6IiMjAmOiJiIgMjImeiIjIwJjoiYiIDMxmiV5EvhCRsyISZbWti4j8V0QOicg6Ealvq/6JiIjItiP6rwCEF9n2GYCZqhoCYDWA523YPxERUZ1ns0SvqjsBXCiyuT2AnZbXWwAMt1X/REREVP336KMADLa8fhCAf0kNRWS8iESISERqamq1BEdERGQ0oqq2O7lIIID1qtrZ8r4DgIUAGgNYC+AZVW1cjvOkAkiy2tQAQHpVx2uHavvnrA3x21OMNRlLdffdBMC5auyvrgpQVe+aDoJqVrU+GU9V4wD0BwARCQJwXzmPu+4/VBFZrKrjqz5C+1LbP2dtiN+eYqzJWKq7bxGJUNXu1dUfUV1WrVP3IuJj+bcJwMsAPq7kqdZVWVD2rbZ/ztoQvz3FWJOx2NP3QERVyGZT9yLyfwD6wjxFdwbAKwA8APzF0mQVgFlqy3sHRGSXOKInqj42vUdPRFQcERmvqotrOg6iuoCJnoiIyMD4CFwiIiIDY6InIiIysDqT6EWko4h8LCIrRGRSTcdDJePfioio6tRIohcRBxHZJyLrb+IcNxTNsdoXLiLxInJERGYCgKrGqupEACMBcLVvOYhIQ0uyjRORWBG5o5Ln4d+KSiUi9UQkUkQG1XQsREZTUyP6ZwHEFrdDRHxExLPItrbFNP0KNxbNgYg4APgQwEAAwQAeEZFgy77BAH4C8O+bCb4OeR/AJlXtAKALivzN+LeikpR0cVfchZ3FCwC+rd4oieqGak/0ItIC5ififVZCkz4AvhcRV0v7p2B+bO51SiiaAwBhAI6o6jFVvQrgGwAPWI5Zq6o9AYy66Q9icJYSwv8D4HMAUNWrqppWpBn/VlSSr1Dk4q6kCzsRuRtADMzP2yCiKlatj8C1eA/ADACexe1U1e9EpBWAb0TkOwDjANxTgfP7AThh9f4kgB4i0hfAMAAuAH6oRNx1TWsAqQC+FJEuACIBPKuqlwsb8G9FJVHVnZZaF9auXdgBgIgUXth5AKgHc/LPFpEfVLWgGsMlMrRqTfSW+29nVTXS8j/zYqnqPyz/E/gIQBtVzaxIN8WfUrcD2F6B89R1jgC6Afirqv4iIu8DmAngb9aN+LeiCij2wk5VnwYAEfkzgHNM8kRVq7qn7nsBGCwiiTBP0/YTkX8VbSQidwLoDGA1zI/OrYiTuL78bQsApysVbd12EsBJVf3F8n4FzIn/OvxbUQUUe2F37YXqV6pa6QW6RFS8ak30qjpLVVuoaiCAhwH8R1VHW7cRkVAAn8I8pTcWgJeIvF6Bbn4D0E5EWomIs6WftVXyAeoQVU0BcEJE2ls23QXzfdRr+LeiCuKFHVENsMff0bsDeFBVj1qm8B7H9bXoAVwrmvNfAO1F5KSIPAEAqpoH4GkAm2FeJf6tqkZXW/TG8lcAy0TkIICuAN4osp9/K6oIXtgR1QA+656Iqlxx1StV9XMRuRfmBbkOAL5Q1bk1FyVR3cBET0REZGD2OHVPREREVYSJnoiIyMCY6ImIiAyMiZ6IiMjAmOiJiIgMjImeiIjIwJjoiYiIDIyJnoiIyMCY6IlqMRH5p4jsFZHbajoWIrJPTPREtZSI1APgA2ACgEE1HA4R2SkmeqoVRORdEZli9X6ziHxm9X6+iEyt4j4zq/h8DUVkstX7QBGJKuexbiKyQ0QcCrep6mUAzQBsB7BQRJxFZKeIOFZl3ERUuzHRU22xB0BPABARE8zFUjpZ7e8JYHcNxFURDQFMLrNV8cYBWKWq+YUbRKQxzBUEMwDkq+pVAP8G8NDNBkpExsFET7XFblgSPcwJPgpAhog0EhEXAB0B7BORNSISKSLRIjK+8GARmVdkNP13EZkmIqNF5FcR2S8in1iPmK3aFtvGMiKPFZFPLf39KCJuln1/E5E4EdkiIv8nItMBvAWgjeU8b1tO71Dc8cUYBeD7ItteBvAOgGgAwZZtayxtiYgAMNFTLaGqpwHkiUhLmBP+fwH8AuAOAN0BHLSMaMep6q2Wbc9YRr0A8A2uH+mOBBBh2dZLVbsCyEeRJCkiHcto0w7Ah6raCUAagOEi0h3AcAChAIZZYgGAmQCOqmpXVX2+pOOLfnZL7fbWqppotS3Q8j0sBxCLP2Y3ogBwYR4RXcN7eVSbFI7qewJYAMDP8jod5ql9wJzch1pe+8OcSM+r6j4R8RGR5gC8AVwEEALgVgC/iQgAuAE4W6TPu8po87uq7re8jgQQCPNthe9VNRsARGRdKZ+puOOLagLzRYC11wHMUVUVkWuJXlXzReSqiHiqakYp/RJRHcFET7VJ4X36EJhHricATANwCcAXItIXwN0A7lDVLBHZDsDV6vgVAEYA8IV5hC8A/ldVZ5XSZ1ltcqxe58N8ISAV+EzFHV9UNqw+h4h0hXmmoLeIfGjZd8iqvQuAKxWIgYgMjFP3VJvshvlnZBdUNV9VL8C8wO0OmKfyGwC4aEnyHQDcXuT4bwA8DHOyXwHzwrURIuIDACLiJSIBRY4pT5uifgJwv4i4iogHgPss2zMAeFb0Q6vqRZjv5Rcm+3kA7lfVQFUNBNAFlhG95VZFqqrmVrQfIjImJnqqTQ7BPI39c5Ft6ap6DsAmAI4ichDAa0XaQVWjYU60p1Q1WVVjYF7Q9qPlmC0w/1zN+pgy2xSlqr8BWAvgAIBVMK8FSFfV8wB2i0iU1WK88voR5hF8PwD1VPXfVv2dAVBPRLwA/AnADxU8NxEZmKhqTcdAZDgi4qGqmSLiDmAngPGquvcmzhcKYKqqjimj3SoAs1Q1vrJ9EZGx8B49kW0sFpFgmO+f/+/NJHkAsCwm3CYiDta/pbdmWZ2/hkmeiKxxRE9ERGRgvEdPRERkYEz0REREBsZET0REZGBM9ERERAbGRE9ERGRgTPREREQGxkRPRERkYEz0REREBvb/kJri0BSPz5EAAAAASUVORK5CYII=\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 SA13')"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfoAAAEYCAYAAAC0mTTAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XlcVXX+P/DXm01AroqKooAsyiqKW2jqpGm5NG1qq2VNWWmOv1KzbZxv01S2a43T1LQ3mpaVS5qpaeNao4YbIouK4gqKG14FV96/P865dkFQ4HK5cH09H4/74N6zvs+9wPt+Pp9zzltUFUREROSePFwdABERETkPEz0REZEbY6InIiJyY0z0REREboyJnoiIyI0x0RMREbkxJnoiIiI3xkRPRETkxpjoqVwikiMiN1TTtl4TkTHVsa0ytr1ORNo6Y9tERHUdE30dYibeIhGxishxEflVREaKSLV8jtWZ2EttNwjAAwA+tJu2XEROi8hJ85HlwC7eBvCSo3GSc4hIT/N3tUBEjorILyJyTalllovIMRGpV8b6o0UkRUTOiMgXpeZ9KSK5InJCRLaJyCNOPhyiOoeJvu65RVUtAMIBvA7gWQCfujakK/oTgB9VtajU9NGqGmA+Yh3Y/jwA14tICwe2UauJiJerY6gKEWkA4AcA/wTQGEAIgL8DOGO3TASAPwBQALeWsZkDAF4B8FkZ814DEKGqDcx1XxGRztV3BER1HxN9HaWqBao6D8DdAB4UkUQAEJGWIjJLRPJFZJeIPGFbx2yxPy8i6Wbr6XMR8TXnTQPQCsB8s4X9jLlaBxFJNVtjM23Lm+s8KyL7zR6GLBHpW064AwGsqOqxioi3iEw04z8nImo+NpvvxWkA6wH0K2Pd50Tku1LT/iEiUyp5DBCRMBGZbb63R0TkPXN6vNkiPS4iW0Xk1ors23x9pc/rWRFJBXBKRLzMbWab8aaLyCC75TuJyEZz3rfm5/VKRfZVxrGWeUx2cY0v7/eilBgAUNWvVPWCqhap6k+qmmq3zAMA1gD4AsCDpTegqrNVdS6AI2XM26qqti8Naj5al3dcRFclVeWjjjwA5AC4oYzpewA8DuOL23oALwDwARAFYCeA/nbrpwEIg9G6+gXAK+Vt33y9DkBLc/kMACPNebEA9gJoab6OANC6nLjzAVxTatpyc/phM47elznuN2AkgjAA9QEsBTAbQJTdMlMATC5j3XAAhQAamK89AeQC6FbJY/AEsBnAO2YMvgB6AvAGsAPAX8z3vA8Aq7ntcvdtvq7I57XJPG4/c9qd5ufhAeNL3ikALcz1dwN40oxpMICzts/3SvsqdazlHtOVfi/K2FYDGAn6PzC+8AWWscwOAKMAdAZwDkDzcrb1CoAvypj+vvk+K4ANAAJc/bfKBx+16cEWvXs4AOMf7jUAglT1JVU9q6o7AXwM4B67Zd9T1b2qehTARAD3XmHbU1T1gLn8fAAdzOkXANQDkCAi3qqao6rZ5WyjEYxEYe9ZGMkmBMBHMHoSLmmJiYgFwBMAhplxnwIwC0Bj8/hsrOZ+SlDV3TD++d9uTuoDoFBV11TyGJJhJLanVfWUqp5W1dUwvjAEAHjdfM//C6Or+t4r7Buo2Oc1xTzuIvN4vjU/j2JVnQlguxlbNwBe5vLnVHU2jGSMSuzLptxjKhVXWb8XJajqCRhfiNTcX76IzBOR5oAxfg/jC9E3qroeQDaAoWVtqzyqOgqABUb3/2zYDQsQEbvu3UUIgKMw/mG2NLtbj4vIcRitsuZ2y+61e74bRvK6nDy754UwEgBUdQeAMQBeBHBIRL4WkfK2dQzGP+KLVHWtqlpV9Yyq/gdGq/6mMta9DsBOVd1uNy2wVFwwt3+8nP3PwO9Jaqj5urLHEAZgt6qeLzW9JYC9qlpsN203jM+k3H2bKvt5QUQeEJFNdssnAmhqxrFfVbWcdSuyr4oeE1DO70VZVDVDVf+kqqFmvC0BvGvOfhDAT6p62Hw9A2V031+JGsMCqwGEwujdIiITE30dJ8bZyyEAVsP4x75LVRvZPSyqap9Aw+yet4LRG2BjnySuSFVnqKqtRaYwutjLkgpzrPZymwMgZUwPgvFFAQAgIgJgEIwWpr14GF3rZfkWQG8RCTXXvZhsK3EMewG0KuOkuAMAwqTklQ+tAOy/0r5Rsc/r4mciIuEwWsWjATRR1UYwhmIExpBAiPn+2Nh/1hXZV0WPqcpUNRPGWHyiiPgBuAtALxHJE5E8AGMBJIlIUhV34QWO0ROVwERfR4lIAxG5GcDXAL5U1S0wumpPmCdw+YmIp4gkSslLmf4sIqEi0hhGi26m3byDMLrTK7L/WBHpI8blUKcBFMHoCi/LjwB62a3bSET6i4iveYLZfTBa7ovLWDcNQCcR6WAmhtdgJL+ZdturB2N8d0lZO1fVfBjnBHwOI9llVOEY1sFIpq+LSH0z9h4A1sIYJ39GjJMGewO4BcbnUu6+7bZ5pc/LXn3z2PPN+B+C0UIGgP+ZsY8239PbYHTpV2Vflz2myhCROBF5yvyiAxEJg9HDsQbGkMYFAAkwuv47wPjCtgrGCXq2bXiZJ/t5AvC0+71pJiL3iEiAeTz9zW3/t7JxErk1V58kwEfFHzBOgiqCMR5dAOOf+58BeNot0xLAVzC6Vo/B+Id6g936zwNIh9HN/R8A/nbr3gbjxL7jAMbj0pPzXoTxpQIA2sNIHlYYwwY/wDyprYy4mwLYh99PKAsC8Ju57nEzxhsvc9wTYLQyc2G0BpuWmn8ngNlXeO+GwUiST9tNq/AxmMu3AmA7+/swjHFqAGgL46qCAvO9HXSlfVfi87qh1PITzVgPA5hs7vcRc14XGCfvnYTRkzAbwP9VZF9lxFXuMV3u96KM7YQA+AZGb8Ap8+eHME7SWwRgUhnr3GXG6GW3fS31eNH8PVph/g6dALAFwKOu/jvlg4/a9hDVSvXWUh0mIjkwksJSF+z7VQCHVPXdKy5c+W2vBTBcVdOqe9t1mfm+/FtVP3d1LETkOnXyJhxU96jqX5y47a7O2nZdIiK9AGTBaO3fB6PHYpFLgyIil2OiJ3IfsTC6yQNgXKZ2h6rmujYkInI1dt0TERG5MZ51T0RE5MbqRNd906ZNNSIiwtVhEBHVKevXrz+sqkEOrN/My8vrExiXcbJhWDsVA0g7f/78I507dz5U1gJ1ItFHREQgJSXF1WEQEdUpIrLbkfW9vLw+CQ4Ojg8KCjrm4eHBcd5aqLi4WPLz8xPy8vI+QdnVH/kNjYiIypUYFBR0gkm+9vLw8NCgoKAC/H7zrEuXcdbOxSjpuUxEMswyl0+Wmj9ejFKjTZ0VAxEROcSDSb72Mz+jcvO5M7vuzwN4SlU3mBXI1ovIElVNN2+DeSOMu7ARERGRkzitRa+quaq6wXxuhVGz2lb96h0Az6CSRVSIiIiocmpkjF5EIgB0BLBWRG6FUU6zvEpjtnUeE5EUEUnJz8+vgSiJiKiq3l6c1XxpxsES5aiXZhy0vL04q6xSyLXClClTmjzwwAOtSk9/8803g957770mpadnZWX5REdHt62Ofd99993h69ev962ObV2J0xO9iAQAmAWj7vd5GAVKXrjSeqr6kap2UdUuQUGVvDpk9bvArpUlp+1aaUwnIqJq16FVo8Jx32yKsiX7pRkHLeO+2RTVoVWjQlfHVlnPPPNM/ujRo484cx8zZ87c3blz59PO3IeNUxO9iHjDSPLTVXU2jDrRkQA2mwVWQgFsEJHgat1xSCfg2z/9nux3rTReh3Sq1t0QEZHhhvjm1sl3ddg57ptNUX+fv7XluG82RU2+q8POG+KbWx3e9g03tG7btm18mzZt2r799ttNAcDf37/jo48+GpqQkBB/7bXXxhw4cKDcc86Sk5NjV65c6Q8Aubm5XiEhIe1KL/P111837NChQ1xubq7XuHHjWr7wwgvNAWDVqlX+sbGxCR06dIibPHlyM9vyKSkpvu3atYuPi4tLiImJSdiyZUu9EydOePTu3btNbGxsQnR0dNuPP/44sCIxOZszz7oXAJ8CyFDVyQCgqltUtZmqRqhqBIzSpZ1UNa9adx55HXDnF0Zy/+9E4+edXxjTiYjIKW6Ib24d0ik0//NfcloM6RSaXx1JHgCmT5+es3Xr1oxNmzalf/jhh83z8vI8i4qKPDp16lSYnp6e0aNHD+tzzz3Xsqrbnzp1aqO33noreMmSJdtbtGhx3n7e8OHDIyZPnrxn06ZNmfbT//nPfwaNGjXqYGZmZnpqampGZGTk2dmzZzcIDg4+l5WVlb59+/atgwcPPlHVmKqTM1v0PWDU4e4jIpvMx01O3F9JkdcBXYYDK980fjLJExE51dKMg5ZZG/YFPdQjInfWhn1Bpcfsq+qNN95oHhsbm9C5c+f4vLw8761bt/p6eHjgkUceOQoADz/88JF169YFVGXbv/76q2XSpEnBS5Ys2R4UFHTBft6RI0c8rVar5x//+MeTtv3Y5l177bWnJk2a1GLChAnB27dv9wkICNBOnToVrVq1qsHjjz8esmjRooAmTZpcKL0/V3DmWferVVVUtb2qdjAfP5ZaJkJVDzslgF0rgZRPgeueMX6WHrMnIqJqYxuTn3xXh51/u6XtAVs3vqPJ/ocffrCsWLHCkpKSkpmVlZUeHx9fVFRUdEnuMjqRy+bl5aUXLhg5t7CwsMSCrVq1OnPq1CnPtLS0S06MU9Vytzty5Mij33///Q4/P7/igQMHxsybN8/Svn37Mxs2bEhv165d0YQJE0LGjx/fopKH6xTueWc825j8nV8AfSb83o3PZE9E5BSb9hz3tx+Tt43Zb9pz3KFx6OPHj3s2bNjwgsViKd64caPv5s2b6wNAcXExPv/880AA+OKLL5okJyeXO0wQFhZ2Zt26dfUBYPr06SXGzUNDQ8/OmjVrx0MPPRSZkpJSItk3bdr0QkBAwIXFixcHmPtpbJuXnp7uEx8ff+avf/3roX79+h3ftGmTX05OjrfFYikeNWrU0TFjxhzctGlTjYzBX0mduNd9pe3fUHJM3jZmv38Du/CJiJxgfP/Yg6Wn3RDf3OroOP2QIUMKPvroo6CYmJiE1q1bn05KSjoFAH5+fsVbt271a9u2bbDFYrkwe/bsneVt47nnnjt49913R3399ddN/vCHP1wybp6UlHRm6tSpO+++++7W8+bN22E/79NPP8155JFHIvz8/Ir79Olzcd1p06Y1/vbbb5t4eXlpUFDQuddee+3A6tWr6z///POhHh4e8PLy0vfff/+ytQYu1wtRnepEPfouXbooi9oQEVWOiKxX1S5VXX/z5s05SUlJzhledZC/v3/HwsLCja6Oo6piYmIS5s2btyMuLu5sdWxv8+bNTZOSkiLKmueeXfdERES1VPfu3aNjY2OLqivJX4l7dt0TEZFbK6s1P2zYsFa//fZbibPvH3/88YNPPvmkU29+czk33nhj671799aznzZx4sR9Q4YMqbFL75joiYjILUybNq3WFUpbsmRJtqtjYNc9ERGRG2OiJyIicmNM9ERERG6MiZ6IiMiNMdETEZHjfn65ObIWlrzdbdZCC35+2eX16P/3v//5dejQIS4mJiahT58+bY4ePXox9z3//PPBrVq1SoyIiEicNWtWg8pue8yYMS3nzp1bLff0dxYmeiIiclxol0LMGRl1MdlnLbRgzsgohHZxeT36Rx99NGLixIn7tm3bln7rrbce+/vf/x4MAOvXr/edPXt246ysrK2LFi3aNmbMmFbnz5+/0uYuOn/+PN59990Dt99+e7VU6XMWJnoiInJc7EArBv17J+aMjMLC51pizsgoDPr3TsQOdHk9+pycHN+BAweeBICbb775xA8//BAIAN99912jwYMHH/Xz89O4uLiz4eHhZ5YvX14/KyvLJzIysu3gwYMjYmJiEgYMGBBltVo9ACAkJKTd+PHjW3Tu3Dn2s88+CxwyZEiE7Z77ISEh7UaPHh3SoUOHuMTExPjVq1f79+zZMzosLCzxzTffDLLF83//93/NExMT42NiYhLGjh172fK6Tz/9dIvIyMi23bt3j77lllsiX3jhhUr3kDDRExFR9YgdaEXSvflY+0ELJN2bXx1JHnC8Hn10dHTRjBkzGgHAl19+2TgvL88HAPbv3+8TFhZ28e50LVu2PLt3714fwPhyMHLkyPxt27alWyyW4rfeeutiovb19S1ev3591mOPPXas9L7CwsLObtq0KbNr164nH3744Yj58+dnr127NvP1119vCQCzZ89usGPHDt/U1NSMjIyM9E2bNvkvXLiwzBK7K1eu9J8/f37gli1b0hcsWJCdmppavyrvHxM9ERFVj6yFFmz+KghdH8/F5q+CLhmzryJH69F/9tlnOR988EFQ27Zt461Wq4e3t7cCRhna0kREASA4OPhsv379TgHAsGHDjvz6668Xt//AAw9ckuBt7rrrruMA0K5du8JOnTqdCgwMLG7ZsuX5evXqFR8+fNhz0aJFDVauXNkgISEhoW3btgnZ2dm+mZmZl5TIBYDly5cHDBw48HhAQIAGBgYW33jjjccr9IaVwjvjVcJnaZ8hsUkiklskX5y2Lncd0o6k4eHEh10YGRGRi9nG5G3d9VG9rNXRfW9fj95isRQnJyfHVrYefceOHU//8ssv2wEgNTW13k8//dQIMErU2lrwAHDgwAGf0NDQc2Vtz/61xWIpLm9fvr6+CgAeHh7w8fG5+E3Cw8MD586dE1XFmDFjcp9++ukrFguqrqJzbNFXQmKTRIxfMR7rctcBMJL8+BXjkdgk0cWRERG52L4U/xJJ3TZmvy/F5fXo9+/f7wUAFy5cwN/+9rcWw4cPPwQAQ4YMOT579uzGRUVFkpmZ6ZOTk+Pbu3fvUwCQm5vrs3Tp0voAMGPGjMbdu3c/6chx2AwcOPDEtGnTmhYUFHgAwK5du7xt8ZXWu3fvk4sXL25YWFgoBQUFHkuXLm1UlX2yRV8JyS2S8XavtzF+xXjcFXsXvsn6Bm/3ertEC5+I6KrU9/8uqUeP2IFWR8fpq6Me/Weffdb4008/bQYAN91007EnnnjiCAB06dLl9O233340JiamraenJyZPnrzby8tIi1FRUac/++yzJqNGjQqPjIw8M378+HxHjsNm8ODBJ7Zu3ep7zTXXxAGAv79/8fTp03eFhIRccrp/r169CgcMGFCQkJDQNiQk5Ez79u1PNWzY8EJl98l69FXw3sb38GHqhxjRfgRGdxzt6nCIiMrEevRVk5WV5XPzzTdHb9++fasztl8ZBQUFHg0bNiy2Wq0e1157bey///3v3T179rzkksXL1aNni76S1uWuwzdZ32BE+xH4JusbJAcns0VPREROcf/994dv377d78yZM3LPPfccKSvJXwkTfSXYxuRt3fXJwcklXhMRUc1wZj362NjYszXZms/Ly/Ps3bt3bOnpy5cvz5o/f/4uR7fPRF8JaUfSSiR125h92pE0JnoiIherjfXoKyI4OPhCZmZmurO2z0RfCWVdQpfcgl33RERUe/HyOiIiIjfGRE9EROTGmOiJiIjcGBM9ERE5bMqGKc2X711e4t72y/cut0zZMIX16F2MiZ6IiBzWPqh94YTVE6JsyX753uWWCasnRLUPas969NWguLgYFy5U+qZ4AJjoiYioGvQO622d2HPizgmrJ0S9vu71lhNWT4ia2HPizt5hvVmPvor16LOysnyioqLa3n///a3MSnc+5S17OUz0RERULXqH9bbe0vqW/OkZ01vc0vqW/OpI8sDVW4/eFsdDDz10JCMjIz0mJuZsectdDhM9ERFVi+V7l1vmZ88Pui/+vtz52fODSo/ZV9XVWo8eAFq0aHG2b9++p674Jl0Gb5hDREQOs43J27rru7XoZq2O7vuruR49YFS3q8hyl8MWPREROSw1P9XfPqnbxuxT81NZj95OZerRVxenbVxEwgBMBRAMoBjAR6r6DxF5GcBt5rRDAP6kqgecFQcRETnfE52euKQefe+w3lZHx+mv5nr01cVp9ehFpAWAFqq6QUQsANYDuB3APlU9YS7zBIAEVR15uW3Vtnr0RER1AevRV01tqkdfUZerR++0rntVzVXVDeZzK4AMACG2JG+qD8A53zSIiIioZk7GE5EIAB0BrDVfTwTwAIACANeXs85jAB4DgFatWtVEmEREVEdcLfXog4ODq3aXHDtO67q/uAORAAArAExU1dml5j0PwFdV/3a5bbDrnoio8ty5655KcknXPQCIiDeAWQCml07yphkAhjgzBiIioquZ0xK9GBcdfgogQ1Un202PtlvsVgCZzoqBiIjoaufMMfoeAIYB2CIim8xpfwEwXERiYVxetxvAZc+4JyIioqpzWqJX1dUAyrpV0Y/O2icRERGVxDvjERGRww69+25z67JlJe5tb122zHLo3XddXo++JtXG+vRM9ERE5DC/pKTCA88+F2VL9tZlyywHnn0uyi8pyeX16GtSTdWnrwwmeiIicpjl+uutLd94feeBZ5+Lynv11ZYHnn0uquUbr++0XH+9y+vRJycnxw4fPjysS5cusVFRUW1XrFjh369fv9bh4eGJTzzxxMXyti+++GLz6OjottHR0W1feumlZgBw4sQJj969e7eJjY1NiI6Obvvxxx8HAsD48eNbJCYmxkdHR7e99957w4uLjdoz9vXpV6xY4d+xY8e42NjYhHbt2sUfO3aszJx79913h8fFxSXExcUlBAYGJj311FMtHH3P7DHRExFRtbBcf7214e235R+bOq1Fw9tvy6+OJA84Xo8eAHx8fIpTUlKyHnroofw777yzzccff7wnMzNz68yZM5vm5eV5rlq1yn/GjBlN1q9fn5GSkpIxderUoF9++cVv9uzZDYKDg89lZWWlb9++fevgwYNPAMDTTz99KC0tLWP79u1bi4qKPL7++uuG9vs7ffq03Hfffa3ffffdPVlZWekrVqzICggIKLMS3cyZM3dnZmamz5s3b0ejRo3OjxgxolI3+LkSJnoiIqoW1mXLLAVzvw8KfGBYbsHc74NKj9lXlaP16AFg0KBBxwEgKSmpqE2bNkXh4eHn/Pz8NCws7MzOnTt9li9fHnDTTTcdb9CgQXHDhg2L//jHPx5btmyZpVOnTkWrVq1q8Pjjj4csWrQooEmTJhcAYOHChZb27dvHxcTEJPz666+WtLQ0P/v9paam+jZr1uxcr169CgGgcePGxd7e3uXGV1hYKEOGDGn9zjvv7ImJiTnr4FtWAhM9ERE5zDYm3/KN13cG/+UvB2zd+I4me/t69FlZWenx8fFFla1HD5SsE1+vXr0SdeLPnz9f7l1i27dvf2bDhg3p7dq1K5owYULI+PHjWxQWFspTTz0VPnv27Oxt27al33///YdPnz5dIiZVhYhU+Nazw4YNC7/llluOOWN8n4meiIgcVrR5s7/9mLxtzL5o82aX16OviD59+pz88ccfG1mtVo8TJ054/Pjjj4HXX3+9NScnx9tisRSPGjXq6JgxYw5u2rTJv7Cw0AMAgoODzxcUFHjMnz8/sPT2kpKSTh88eNBnxYoV/gBw7Ngxj3PnzpW579deey3o5MmTnq+++mqeI8dQnhopakNERO6t2Zgxl9Sjt1x/vdXRcfrqqEdfET179iwcOnTokU6dOsUDwLBhw/J79OhRNGvWrAbPP/98qIeHB7y8vPT999/f3bRp0wv33XdffkJCQtvQ0NCztpjs+fr66vTp07OfeOKJVqdPn/bw9fUtXrly5baGDRteMk7/3nvvBXt7e2tcXFwCADz88MP5zzzzTL4jx2PP6UVtqgOL2hARVZ47F7VxZj36ushlRW2IiIjItdh1T0REdY4z69E7y6xZsxpMmDAh1H5aWFjYmSVLlmQ7c79M9ERE5BamTZu2x9UxXM6QIUNODBkyJL2m98uueyIiIjfGRE9EROTGmOiJiIjcGBM9ERE5bM332c13pR4ucRe8XamHLWu+z76qytTWRkz0RETksOaRDQt//iI9ypbsd6Uetvz8RXpU88iGTilTm5ycHLty5cpL7rrXq1evNocPH/YsPX3cuHEtX3jhBYe/dOTk5HgPGDAgytHt1CSedU9ERA6LbN/U2vdPCTt//iI9KrZbcH7Wmrygvn9K2BnZvmmN1mZfsWLFDmduPyIi4tyiRYscugtfTWOLnoiIqkVk+6bW2G7B+an/3dcitltwfnUk+aysLJ/o6Oi2ttcvvPBC83Hjxl0sSXvhwgUMHjw4wlZXPiQkpF1ubq4XADz77LPBERERid27d4/Zvn17Pds6r7zySrPWrVu3jYmJSbj55pujAGDBggUBtprw8fHxCeXVji8dT13AFj0REVWLXamHLVlr8oLa9wnNzVqTFxQa19jqzBb9uXPn5Pbbb49MSEgoeuONN0oUhFm1apX/nDlzGm/ZsiX93Llz6NChQ0LHjh0LAWDKlCnBu3fv3uLn56e2bv5JkyYFT5kyZXe/fv1OFRQUePj7+5dZO74uYoueiIgcZhuT7/unhJ1/uCvmgK0bv/QJetVp1KhR4WUleQBYtmxZwE033XTcYrEUN27cuLhfv37HbfNiY2OLBg0aFPn+++839vb2VgDo1q3byfHjx4e98sorzQ4fPux5udrxdQ0TPREROezgrgJ/+zF525j9wV0FDpWp9fLy0uLi3xvX9nXfu3TpcnLVqlUNCgsLyyxGX16N+mXLlm3/85//nL9+/fr6SUlJCefOncOrr76a98knn+wuKiry6N69e/zGjRt9HYm7NmGiJyIih3W7rfXB0t30ke2bWrvd1vqS8rWVERoaev7o0aNeeXl5nkVFRbJ48eKGtnkjRow43K9fv4Kbb765dela73369Dm5YMGCRidPnpRjx455LFmypBFgjOlnZ2f73HLLLdb3339/n9Vq9SwoKPDcunVrveTk5KKJEyfmtWvX7lRaWprbJHqO0RMRUa1Vr149feqpp3KTk5PjQ0NDz7Rp0+a0/fwXX3zx4NixYz0HDx4cOXfu3F226T179iwcNGjQ0cTExLYhISFnkpOTTwLA+fPnZejQoZFWq9VTVWXEiBEHmzZteuGpp56zGwbfAAAfnklEQVRq+euvvzbw8PDQmJiYojvuuKOgvJhEpPbXd7fDevRERG7KnevRu8qqVav8x40bF/bbb79luToWe6xHT0RE5KCVK1f6Dxs2LGr06NEODUfUNHbdExERlbJu3Tq/Bx54INJ+mo+PT3FOTk6aq2KqKib6SjjyySfwTWyH+t26Xpx2as1anE7bgiaPPOLCyIiIqDolJycXZWZm1njteGdg130l+Ca2w/6xY3FqzVoARpLfP3YsfBPbuTgyIiKisrFFXwn1u3VFyDvvYP/YsQi89x4c++prhLzzTokWPhERUW3CFn0l1e/WFYH33oPD73+AwHvvYZInIqJajYm+kk6tWYtjX32NpqMex7Gvvr7YjU9EdDVb/fXU5tnr15W43W32+nWW1V9PZT16F2OirwTbmHzIO+8g6IknLnbjM9kT0dWuRXRc4cJ/TYqyJfvs9essC/81KapFdJzD9ehffvnlZtHR0W3btGnT9qWXXmoGAAcPHvTs3r17dHh4eGL37t2j8/PzL6lBfyUdO3aMczS2usBpiV5EwkRkmYhkiMhWEXnSnP6WiGSKSKqIzBGRRs6KobqdTttSYkzeNmZ/Om2LiyMjInKt1p2TrQP//NTOhf+aFLXsi49aLvzXpKiBf35qZ+vOyQ5Vr/vtt998p06dGrRhw4aMjIyMrYsWLWq0ZcuWen/7299a9O7d27p79+603r17W1944YXgim7z/PnzAICNGzdmOhJbXeHMFv15AE+pajyAbgD+LCIJAJYASFTV9gC2AXjeiTFUqyaPPHLJmHz9bl15aR0REYxk3/a6vvkbFs5r0fa6vvmOJnkA2LJli1+nTp1OWiyWYm9vb/To0cM6c+bMRosWLWo0YsSIIwAwYsSIIwsXLgwEgHHjxrW8/fbbI7t16xYTHh6eOGnSpKYA8MMPP1i6du0ac8stt0TGxsa2BQB/f/+OtnnXXHNN7E033RQVERGROGrUqJAPPvigcbt27eJjYmIStm7dWg8ADhw44NW/f//WiYmJ8YmJifE//fRT/fLiPnDggFf37t2jExIS4ocOHRresmXLdrm5uS45Ad5piV5Vc1V1g/ncCiADQIiq/qSq583F1gAIdVYMRERUc7LXr7NsXflzUKeBt+ZuXflzUOkx+6ro0KFD0dq1ay15eXmeVqvVY8mSJQ337t3rc+TIEa/w8PBzABAeHn7u6NGjF5NoRkaG39KlS7evWbMm86233mqZk5PjDQCpqan133rrrf3Z2dlbS+8nMzPT74MPPtibkZGx9bvvvmuybds23y1btmQMGzbs8KRJk5oBwIgRI8LGjRt3MC0tLWPOnDnZI0eOjCgv7ueee65lr169rOnp6RmDBw8+lpub6+Poe1FVNfLtQkQiAHQEUHow+2EAM8tZ5zEAjwFAq1atnBgdERE5yjYmb+uub9Wug7U6uu87dep0+sknn8zr06dPjL+/f3FCQkKhl9flU9fAgQOPBwQEaEBAwPlrr732xKpVq+oHBgZeaN++/am4uLizZa3Trl27U7YvDq1atTozcODAAgBISkoqWrFihQUAfvnllwbbt2/3s61z8uRJz2PHjnkEBgYWl97eunXrAubOnbsDAO64444TDRo0uFDV98BRTj8ZT0QCAMwCMEZVT9hNnwCje396Weup6keq2kVVuwQFBTk7TCIickDu9kx/+6RuG7PP3Z7pUD16ABg7duzh9PT0jJSUlKzGjRtfiI6OPt2kSZPzu3fv9gaA3bt3ezdu3NjWU3xJHXrba39//0sSsk29evUuVnjz8PCAr6+v2p5fuHBBAEBVkZKSkpGZmZmemZmZfujQodSykrxt2drCqYleRLxhJPnpqjrbbvqDAG4GcJ/WpneDiIiqpOc9Dxws3XJv3TnZ2vOeBxwuALN//34vANi+fbvPggULGg0fPvxo//79j3/44YdNAODDDz9sMmDAgOO25RcuXNiosLBQ8vLyPNesWWPp2bPnKUdjAICePXueeOONN5rZXv/6669+5S2bnJx8ctq0aY0BYPbs2Q1OnDhR6asCqovTuu7F+Ar1KYAMVZ1sN30AgGcB9FJVhy+7ICIi93brrbe2Pn78uJeXl5e+++67e4KCgi78/e9/zx00aFDr8PDwpi1btjw7d+7cbNvyHTt2PNW3b9/oAwcO+IwfPz43IiLiXFpamq+jcXz00Ud7H3nkkVYxMTEJFy5ckK5du1q7d+++p6xlX3/99QN33HFHVEJCQuC11157Migo6FyjRo1c0n3vtHr0ItITwCoAWwDYujb+AmAKgHoAjpjT1qjqyMtti/XoiYgq72qsRz9u3LiWAQEBF1566SWXlpItKioSLy8v9fb2xtKlS+uPHj063JlFci5Xj95pLXpVXQ1Aypj1o7P2SUREVBvs2LHD56677mpdXFwMb29v/fDDD3NcFQuL2hARkduYPHnygZrc3z/+8Y8mH3zwQYnb/F5zzTUnp02bticjI6NWlLmtUKIXkXoAhgCIsF9HVV9yTlhERES135NPPnnkySefPHLlJV2noi367wEUAFgP4IzzwiEiIqLqVNFEH6qqA5waCREREVW7il5H/6uItHNqJERERFTtLtuiF5EtANRc7iER2Qmj614AqFmYhoiIrnIFi3Oa+7SyFPrFN7l405yijCOWs3us/g37R7j0Urer3ZVa9DcDuAXAQABtAPQzX9umExERwaeVpfDoN9uiijKOWAAjyR/9ZluUTysL69G72GUTvaruVtXdAF6xPbefVjMh1h6rV6/Grl27SkzbtWsXVq9e7aKIiIhqB7/4JtbGd8XsPPrNtqjj87NbHv1mW1Tju2J22rfwq8Id6tGfO3euJnZTroqO0be1fyEingA6V384tVtISAi+/fbbi8l+165d+PbbbxESEuLiyIiIXM8vvom1fqdm+Sd/OdCifqdm+Y4meaDu1qMfN25cy3vvvTe8R48e0YMHD4509H1wxJXG6J+HcdtaPxE5gd/vdHcWwEdOjq3WiYyMxJ133olvv/0WXbp0QUpKCu68805ERrr0MyQiqhWKMo5YTm04FBTQo2XuqQ2Hguq1aWR1NNl36NCh6KWXXgrJy8vzrF+/vi5ZsqRhUlLSqSvVo1+/fn2G1Wr17NixY8KQIUMKAKMe/caNG7eWVao2MzPT77vvvtvZrFmz8+Hh4e3q1at3eMuWLRkvv/xys0mTJjX77LPP9trq0ffv3//k9u3bffr37x+9c+fOS2rb26SmpvqvXbs2MyAgwKXF2y6b6FX1NQCvichrqvp8DcVUq0VGRqJLly5YuXIlrrvuOiZ5IiL8PiZv666v16aRtTq67+tqPXoAGDBgwHFXJ3mg4tfR/0VEBgPoCeMs/FWqOtd5YdVeu3btQkpKCq677jqkpKQgMjKSyZ6Irnpn91j97ZO6bcz+7B6rv6Ot+rFjxx4eO3bsYQAYPXp0SGho6FlbPfrw8PBzNV2PvqLJu379+uXuryZVdIz+XwBGwqhElwZgpIj8y2lR1VK2Mfk777wTffr0udiNX/oEPSKiq03D/hEHSyd0v/gm1uq4tK4u1qOvTSraou8FIFHNmrYi8h8YSf+qsn///hJj8rYx+/3797NVT0TkJHWxHn1tUqF69CIyG8BY87I6iEg4gNdV9V4nxweg9tSjt67YC+9QC3xbN7o47XT2cZzbZ4WlV5gLIyMiuhTr0V89LlePvqJd900AZIjIchFZDiAdQJCIzBORedUTZu3nHWrB0RkZOJ1t9BCdzj6OozMy4B1qcXFkREREZato1/0LTo2ijvBt3QiNh8bj6IwM1O/aAqfW5qLx0PgSLXwiInKd2lSPvibjuJwKJXpVXWF210er6lIR8QPgpaoO3wyhrvFt3Qj1u7aA9b97YekTxiRPRHQVqwv16CvUdS8ijwL4DsCH5qRQAFfl5XWns4/j1NpcWPqE4dTa3Ivd+ERERLVRRcfo/wygB4ATAKCq2wE0u+wabsg2Jt94aDwa9ou42I3PZE9ERLVVRRP9GVW9eDchEfGCceOcq8q5fdYSY/K2Mftz+666EQwiIqojKproV4iI7Z73NwL4FsB854VVO1l6XTom79u6ES+tI6Kr3s8//9w8KyurxCVIWVlZlp9//rl5eevUNu5atraiif45APkwbpIzAsCPAP7qrKCIiKhuCQ0NLZwzZ06ULdlnZWVZ5syZExUaGupwPfqaUlNla2tahRK9qhbDOPlulKreoaofa0XutENERFeF2NhY66BBg3bOmTMnauHChS3nzJkTNWjQoJ2xsbEOjW1mZWX5REZGtr377rvDo6Oj2956662Rc+fOtXTq1CkuPDw8cdmyZf4HDx70vOGGG1rHxMQkJCUlxa1du9YPABYsWBAQFxeXEBcXlxAfH59w7Ngxj4KCAo9rr702JiEhIT4mJibhyy+/vNhNaytbCwB//etfm8fExCTExsYmjBo1qsxa5Dk5Od627cfFxSV4enp23rZtm48jx+sMVypTKwD+BmA0jBK1IiIXAPxTVV+qgfiIiKiOiI2NtSYlJeWvXbu2RdeuXXMdTfI2e/fu9Z05c+bOzp07727fvn389OnTm6SkpGTOmDGj0cSJE1uEhIScTUpKKly6dGn2vHnzLA8++GBkZmZm+qRJk4KnTJmyu1+/fqcKCgo8bEVtFixYsKNx48bFubm5Xl27do0bOnTocQ+P39u933zzTYMFCxYErl+/PtNisRQfPHjQs6y4IiIizmVmZqYDwGuvvRa0atUqS0xMTJnV8VzpSi36MTDOtr9GVZuoamMAXQH0EJGxTo+OiIjqjKysLMvmzZuDunbtmrt58+ag0mP2VRUSEnImOTm5yNPTEzExMUV9+vQ54eHhgU6dOhXu27ev3rp16yzDhw8/AgC33nqr9fjx415Hjhzx7Nat28nx48eHvfLKK80OHz7s6e3tjeLiYhkzZkxoTExMwvXXXx9z6NAhn3379pVo9C5ZsqTB/ffff9hisRQDQPPmzS9cLr6ffvqp/tSpU4O++uqrnOo43up2pUT/AIB7VfVieTZV3QngfnMeERHRxTH5QYMG7Rw4cOABWzd+dSR7Hx+fMkvIenp64sKFC2XWbBERffXVV/M++eST3UVFRR7du3eP37hxo++HH37Y+MiRI15btmzJyMzMTG/SpMm5oqKiErlQVS8pdVue3bt3e48YMSJi5syZ2Q0bNqwVZWlLu1Ki91bVSwoaqGo+AG/nhOS4dd9/hz1pqSWm7UlLxbrvv3NRRERE7m3fvn3+9mPytjH7ffv2+Tt73926dbN+/vnnTQDghx9+sAQGBp5v3Lhx8datW+slJycXTZw4Ma9du3an0tLSfAsKCjybNm16rl69ejp//nzLgQMHLhlTHzBgwIlp06Y1tVqtHgBQXtf9mTNnZPDgwVEvv/zy/vbt259x7lFW3ZUS/eXGGmrdOIRNcOsY/PDu6xeT/Z60VPzw7usIbh3j4siIiNxT3759D5Yek4+NjbX27dvX6VXk3njjjQMbNmzwj4mJSZgwYULIF198sQsA3nzzzWbR0dFtY2NjE/z8/IrvuOOOgkceeeTo5s2b6ycmJsZ/+eWXjSMjI0+X3t4dd9xxYuDAgcc7dOgQHxcXl/Dyyy8Hl7XfpUuX1k9LS6v/yiuvtLSdkJeTk1PrGsGXLVNrnnh3qqxZAHxVtUYOqCplam3JPanfTdj804+4ecxzaJXY3kkREhHVPldjmdqr1eXK1F72rHtVLbO7oi5oldgeSf1uwppZX6PbkHuY5ImI6KpU0TK1dc6etFRs/ulHdBtyDzb/9CPCEtoz2RMRUZUMGzas1W+//RZgP+3xxx8/WNsr1wFumuht3fa27vqwhPYlXhMREVVGbaovX1kVvQVupYlImIgsE5EMEdkqIk+a0+80XxeLSJXHji4nL3tbiaTeKrE9bh7zHPKytzljd0RE7qq4uLi4YteZkcuYn1G5l/Y5s0V/HsBTqrpBRCwA1ovIEgBpAAbj99r21S75tjsumdYqkV33RESVlJafn58QFBRU4OHhwdue10LFxcWSn5/fEEZuLZPTEr2q5gLINZ9bRSQDQIiqLgFQ4ZsREBGRa5w/f/6RvLy8T/Ly8hLhxB5gckgxgLTz588/Ut4CNTJGLyIRADoCWFsT+yMiIsd17tz5EIBbXR0HOcbp39BEJADALABjVPVEJdZ7TERSRCQlPz/feQESERG5MacmehHxhpHkp6vq7Mqsq6ofqWoXVe0SFBTknACJiIjcnDPPuhcAnwLIUNXJztoPERERlc+ZY/Q9AAwDsEVENpnT/gKgHoB/AggCsEBENqlqfyfGQUREdNVy5ln3q2HcE78sc5y1XyIiIvodL5cgIiJyY0z0REREboyJnoiIyI0x0RMREbkxJnoiIiI3xkRPRETkxpjoiYiI3BgTPRERkRtjoiciInJjTPRERERujImeiIjIjTHRExERuTEmeiIiIjfGRE9EROTGmOiJiIjcGBM9ERGRG2OiJyIicmNM9ERERG6MiZ6IiMiNMdETERG5MSZ6IiIiN8ZET0RE5MaY6ImIiNwYEz0REZEbY6InIiJyY0z0REREboyJnoiIyI0x0RMREbkxJnoiIiI3xkRPRETkxpjoiYiI3BgTPRERkRtjoiciqmX+vSIbv2YfLjHt1+zD+PeKbBdFRHUZEz0RUS3TPrQhRs/YeDHZ/5p9GKNnbET70IYujozqIqclehEJE5FlIpIhIltF5ElzemMRWSIi282fgc6KgYioLureuineG9oRo2dsxOSfsjB6xka8N7Qjurdu6urQqA5yZov+PICnVDUeQDcAfxaRBADPAfhZVaMB/Gy+JiIiO91bN8X9XVthyn934P6urZjkqcqcluhVNVdVN5jPrQAyAIQAuA3Af8zF/gPgdmfFQERUV/2afRhfrt2DJ/q0wZdr91wyZk9UUTUyRi8iEQA6AlgLoLmq5gLGlwEAzcpZ5zERSRGRlPz8/JoIk4ioVrCNyb83tCPG9Yu92I3PZE9V4fRELyIBAGYBGKOqJyq6nqp+pKpdVLVLUFCQ8wIkIqplUvcVlBiTt43Zp+4rcHFkVBd5OXPjIuINI8lPV9XZ5uSDItJCVXNFpAWAQ86MgYiorhnZq/Ul07q3bspxeqoSZ551LwA+BZChqpPtZs0D8KD5/EEA3zsrBiIioqudM1v0PQAMA7BFRDaZ0/4C4HUA34jIcAB7ANzpxBiIiIiuak5L9Kq6GoCUM7uvs/ZLREREv+Od8YiIiNwYEz0REZEbY6InIiJyY0z0REREboyJnoiIyI0x0RMREbkxJnoiIiI3xkRPRETkxpjoiYiI3BgTPRERkRtjoiciInJjTPRERERujImeiIjIjTHRExERuTEmeiIiIjfGRE9EROTGmOiJiIjcGBM9ERGRG2OiJyIicmNM9ERERG6MiZ6IiMiNMdETERG5MSZ6IiIiN8ZET0RE5MaY6ImIapkNi3djX9axEtP2ZR3DhsW7XRQR1WVM9EREtUyziAZY/HHaxWS/L+sYFn+chmYRDVwcGdVFXq4OgIiISgqNDUT/RxOx+OM0JF4XgrSV+9H/0USExga6OjSqg9iiJyKqhUJjA5F4XQhSfsxB4nUhTPJUZUz0RES10L6sY0hbuR9dbopA2sr9l4zZE1UUEz0RUS1jG5Pv/2giut4adbEbn8meqoKJnoioljmUc6LEmLxtzP5QzgkXR0Z1EU/GIyKqZTr1D79kWmhsIMfpqUrYoiciInJjTkv0IvKZiBwSkTS7aUki8j8R2SIi80WEF4USERE5kTNb9F8AGFBq2icAnlPVdgDmAHjaifsnIiK66jkt0avqSgBHS02OBbDSfL4EwBBn7Z+IiIhqfow+DcCt5vM7AYSVt6CIPCYiKSKSkp+fXyPBERERuRtRVedtXCQCwA+qmmi+jgMwBUATAPMAPKGqTSqwnXwAta2aQ1MAh10dRBUxdtdg7K5Rl2MHHIs/XFWDqjMYqntq9PI6Vc0E0A8ARCQGwB8ruF6t+0UVkRRV7eLqOKqCsbsGY3eNuhw7UPfjJ9er0a57EWlm/vQA8FcA/67J/RMREV1tnHl53VcA/gcgVkT2ichwAPeKyDYAmQAOAPjcWfsnIiIiJ3bdq+q95cz6h7P2WcM+cnUADmDsrsHYXaMuxw7U/fjJxZx6Mh4RERG5Fm+BS0RE5MaY6ImIiNwYE30pIjJARLJEZIeIPFfG/FYiskxENopIqojcZE73FpH/mPfxzxCR52th7OEi8rMZ93IRCbWb96CIbDcfD9Zs5FWPXUQ6mPUTtprz7q7p2M04qvzem/MbiMh+EXmv5qK+uG9Hfm9aichP5u98unnvjLoS+5vm702GiEwREanh2C+pB1Jqvphx7TDj72Q3z6V/r1THqCof5gOAJ4BsAFEAfABsBpBQapmPADxuPk8AkGM+Hwrga/O5P4AcABG1LPZvATxoPu8DYJr5vDGAnebPQPN5YB2JPQZAtPm8JYBcAI1q4e9NmfHbzf8HgBkA3qtLsQNYDuBG83kAAP+6EDuA7gB+MbfhCeMKod41/N5fB6ATgLRy5t8EYCEAAdANwFpzukv/Xvmoew+26EtKBrBDVXeq6lkAXwO4rdQyCsBWda8hjMsEbdPri4gXAD8AZwGccH7IF1Uk9gQAP5vPl9nN7w9giaoeVdVjMOoQlC5I5ExVjl1Vt6nqdvP5AQCHANT0DZYcee8hIp0BNAfwUw3EWlqVYxeRBABeqroEAFT1pKoW1kzYABx73xWAL4wvCPUAeAM46PSI7WjZ9UDs3QZgqhrWAGgkIi3g+r9XqmOY6EsKAbDX7vU+c5q9FwHcLyL7APwI4P+Z078DcApGi3IPgLdV9XJ/xNWtIrFvxu+FhAYBsIhIkwqu60yOxH6RiCTD+Med7aQ4y1Pl+M2bR02C6yo5OvLexwA4LiKzzaGst0TE0+kR/67Ksavq/2Ak/lzzsVhVM5wcb2WVd3yu/nulOoaJvqSyxuhKX394L4AvVDUURtfaNPOfdTKACzC6jyMBPCUiUc4MtpSKxD4eQC8R2QigF4D9AM5XcF1nciR2YwNGS2cagIdUtdhZgZbDkfhHAfhRVffCNRyJ3QvAH8z518DoQv+T0yK9VJVjF5E2AOIBhMJIkn1E5DpnBlsF5R2fq/9eqY6p0Xvd1wH7ULKiXih+75q3GQ6zm0xV/ycivjCKTgwFsEhVzwE4JCK/AOgCY/ysJlwxdrNrezAAiEgAgCGqWmD2TvQute5yZwZbSpVjN183ALAAwF/NLs6a5sh7fy2AP4jIKBhj3D4iclJVLzmxzEkc/b3ZqKo7zXlzYYwlf1oTgcOx2B8DsEZVT5rzFsKIfSVqj/KOz9V/r1THsEVf0m8AokUkUkR8ANwDo8qevT0A+gKAiMTDGOfLN6f3Mc+UrQ/jn0ZmjUVegdhFpKnZ+wAAzwP4zHy+GEA/EQkUkUAYhYcW11DcgAOxm8vPgTGW+W0NxmyvyvGr6n2q2kpVI2C0PqfWYJIHHPu9+Q1AoIjYzonoAyC9BmK2cST2PTBa+l4i4g2jtV/buu7nAXjA/J/SDUCBqubC9X+vVNe4+mzA2vaA0R2/DcY47wRz2ksAbjWfJ8A4W3czgE0A+pnTA2Cc4bsVxj+7p2th7HcA2G4u8wmAenbrPgxgh/l4qK7EDuB+AOfMz8L26FBX4i+1jT+hhs+6r4bfmxsBpALYAuALAD51IXYYZ9p/CCO5pwOY7IL3/SsY5wecg9FKHw5gJICR5nwB8C/z2LYA6GK3rkv/XvmoWw/eApeIiMiNseueiIjIjTHRExERuTEmeiIiIjfGRE9EROTGmOiJiIjcGBM9XRVEJFhEvhaRbLPK2o8iEuPquIiInI2JntyeWX50DoDlqtpaVRMA/AVGIZnq3ldN3uudiOiKmOjpanA9gHOq+m/bBFXdBGC1WYglTUS2iFnLXkRmishNtmVF5AsRGSIinubyv5n1wUeY83uLyDIRmQHjxiYQkbkisl6MeueP2W1ruIhsE6M2+sdi1p8XkSARmWVu+zcR6VEj7wwRuT3e656uBokA1pcxfTCADgCSYNQr+E1EVsIod3o3gB/NW6v2BfA4jDuXFajqNSJSD8AvImIrLZsMIFFVd5mvH1bVoyLiZ253FoxyqP8Howa5FcB/YdxhETDq0b+jqqtFpBWMW5rGV99bQERXKyZ6upr1BPCVql4AcFBEVsCowrYQwBQzmQ8AsFJVi0SkH4D2InKHuX5DANEAzgJYZ5fkAeAJERlkPg8zlwsGsELN8sUi8i2MUq8AcAOABGOUAQDQQEQsqmqt/sMmoqsJEz1dDbbCuOd5aWWV+4SqnhaR5QD6w2jZf2W3/P9T1RIFRESkN4BTpV7fAOBaVS00t+Vb3v5MHubyRVc+HCKiiuMYPV0N/gugnog8apsgItcAOAbgbnPsPQjAdQDWmYt8DeAhGPXWbYl9MYDHzWpnEJEYs1JhaQ0BHDOTfByMSoYwt93LrDrmBWCI3To/ARhtF18Hh46YiMjEFj25PVVVsxv9XRF5DsBpADkAxsCoOrgZgAJ4RlXzzNV+AjAVwDxVPWtO+wRABIAN5pn8+QBuL2OXiwCMFJFUAFkA1phx7BeRVwGshVFXPB1AgbnOEwD+Za7jBaMu+shqeQOI6KrG6nVENUhEAlT1pNminwPgM1Wd4+q4iMh9seueqGa9KCKbAKQB2AVgrovjISI3xxY9ERGRG2OLnoiIyI0x0RMREbkxJnoiIiI3xkRPRETkxpjoiYiI3Nj/B6S+KGnfqUk1AAAAAElFTkSuQmCC\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
}