{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# AKARI-SEP 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-24 15:15:11.208479\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 = 'AKARI-SEP'\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_akari-sep_20180221.fits\n"
     ]
    }
   ],
   "source": [
    "TMP_DIR = os.environ.get('TMP_DIR', \"./data_tmp\")\n",
    "OUT_DIR = os.environ.get('OUT_DIR', \"./data\")\n",
    "SUFFIX = find_last_ml_suffix()\n",
    "#SUFFIX = \"20171016\"\n",
    "\n",
    "master_catalogue_filename = \"master_catalogue_{}_{}.fits\".format(FIELD.lower(), SUFFIX)\n",
    "master_catalogue = Table.read(\"{}/{}\".format(OUT_DIR, master_catalogue_filename))\n",
    "\n",
    "print(\"Depth maps produced using: {}\".format(master_catalogue_filename))\n",
    "\n",
    "ORDER = 10\n",
    "#TODO write code to decide on appropriate order\n",
    "\n",
    "field_moc = MOC(filename=\"../../dmu2/dmu2_field_coverages/{}_MOC.fits\".format(FIELD))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Remove sources whose signal to noise ratio is less than five as these will have been selected using forced \n",
    "# photometry and so the errors will not refelct the RMS of the map \n",
    "for n,col in enumerate(master_catalogue.colnames):\n",
    "    if col.startswith(\"f_\"):\n",
    "        err_col = \"ferr{}\".format(col[1:])\n",
    "        errs = master_catalogue[err_col]\n",
    "        fluxes = master_catalogue[col]\n",
    "        mask = fluxes/errs < 5.0\n",
    "        master_catalogue[col][mask] = np.nan\n",
    "        master_catalogue[err_col][mask] = np.nan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## I - Group masterlist objects by healpix cell and calculate depths\n",
    "We add a column to the masterlist catalogue for the target order healpix cell <i>per object</i>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#Add a column to the catalogue with the order=ORDER hp_idx\n",
    "master_catalogue.add_column(Column(data=coords_to_hpidx(master_catalogue['ra'],\n",
    "                                                       master_catalogue['dec'],\n",
    "                                                       ORDER), \n",
    "                                   name=\"hp_idx_O_{}\".format(str(ORDER))\n",
    "                                  )\n",
    "                           )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Convert catalogue to pandas and group by the order=ORDER pixel\n",
    "\n",
    "group = master_catalogue.group_by([\"hp_idx_O_{}\".format(str(ORDER))])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#Downgrade the groups from order=ORDER to order=13 and then fill out the appropriate cells\n",
    "#hp.pixelfunc.ud_grade([2599293, 2599294], nside_out=hp.order2nside(13))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## II Create a table of all Order=13 healpix cells in the field and populate it\n",
    "We create a table with every order=13 healpix cell in the field MOC. We then calculate the healpix cell at lower order that the order=13 cell is in. We then fill in the depth at every order=13 cell as calculated for the lower order cell that that the order=13 cell is inside."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "depths = Table()\n",
    "depths['hp_idx_O_13'] = list(field_moc.flattened(13))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=10</i>\n",
       "<table id=\"table140626227758920-283524\" 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>556531712</td></tr>\n",
       "<tr><td>1</td><td>556531713</td></tr>\n",
       "<tr><td>2</td><td>556531714</td></tr>\n",
       "<tr><td>3</td><td>556531715</td></tr>\n",
       "<tr><td>4</td><td>556531716</td></tr>\n",
       "<tr><td>5</td><td>556531717</td></tr>\n",
       "<tr><td>6</td><td>556531718</td></tr>\n",
       "<tr><td>7</td><td>556531719</td></tr>\n",
       "<tr><td>8</td><td>556531720</td></tr>\n",
       "<tr><td>9</td><td>556531721</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(\"$('#table140626227758920-283524').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",
       "    $('#table140626227758920-283524').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=\"table140624562097176-770017\" 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>556531712</td><td>8695808</td></tr>\n",
       "<tr><td>1</td><td>556531713</td><td>8695808</td></tr>\n",
       "<tr><td>2</td><td>556531714</td><td>8695808</td></tr>\n",
       "<tr><td>3</td><td>556531715</td><td>8695808</td></tr>\n",
       "<tr><td>4</td><td>556531716</td><td>8695808</td></tr>\n",
       "<tr><td>5</td><td>556531717</td><td>8695808</td></tr>\n",
       "<tr><td>6</td><td>556531718</td><td>8695808</td></tr>\n",
       "<tr><td>7</td><td>556531719</td><td>8695808</td></tr>\n",
       "<tr><td>8</td><td>556531720</td><td>8695808</td></tr>\n",
       "<tr><td>9</td><td>556531721</td><td>8695808</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(\"$('#table140624562097176-770017').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",
       "    $('#table140624562097176-770017').dataTable({\n",
       "        order: [],\n",
       "        pageLength: 50,\n",
       "        lengthMenu: [[10, 25, 50, 100, 500, 1000, -1], [10, 25, 50, 100, 500, 1000, 'All']],\n",
       "        pagingType: \"full_numbers\",\n",
       "        columnDefs: [{targets: [0, 1, 2], type: \"optionalnum\"}]\n",
       "    });\n",
       "});\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "depths[:10].show_in_notebook()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table masked=True length=10</i>\n",
       "<table id=\"table140624562142008-35025\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>idx</th><th>hp_idx_O_13</th><th>hp_idx_O_10</th><th>ferr_ap_vista_j_mean</th><th>f_ap_vista_j_p90</th><th>ferr_vista_j_mean</th><th>f_vista_j_p90</th><th>ferr_ap_vista_h_mean</th><th>f_ap_vista_h_p90</th><th>ferr_vista_h_mean</th><th>f_vista_h_p90</th><th>ferr_ap_vista_ks_mean</th><th>f_ap_vista_ks_p90</th><th>ferr_vista_ks_mean</th><th>f_vista_ks_p90</th><th>ferr_ap_irac_i1_mean</th><th>f_ap_irac_i1_p90</th><th>ferr_irac_i1_mean</th><th>f_irac_i1_p90</th><th>ferr_ap_irac_i2_mean</th><th>f_ap_irac_i2_p90</th><th>ferr_irac_i2_mean</th><th>f_irac_i2_p90</th><th>ferr_ap_decam_g_mean</th><th>f_ap_decam_g_p90</th><th>ferr_decam_g_mean</th><th>f_decam_g_p90</th><th>ferr_ap_decam_r_mean</th><th>f_ap_decam_r_p90</th><th>ferr_decam_r_mean</th><th>f_decam_r_p90</th><th>ferr_ap_decam_i_mean</th><th>f_ap_decam_i_p90</th><th>ferr_decam_i_mean</th><th>f_decam_i_p90</th><th>ferr_ap_decam_z_mean</th><th>f_ap_decam_z_p90</th><th>ferr_decam_z_mean</th><th>f_decam_z_p90</th><th>ferr_ap_decam_y_mean</th><th>f_ap_decam_y_p90</th><th>ferr_decam_y_mean</th><th>f_decam_y_p90</th></tr></thead>\n",
       "<tr><td>0</td><td>555969277</td><td>8687019</td><td>1.8040026</td><td>44.76484680175781</td><td>4.574854</td><td>91.49818572998048</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>5.3524575</td><td>78.49656219482422</td><td>13.370395</td><td>121.66616973876953</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>0.08941301241009876</td><td>2.0896830548623084</td><td>0.14499704199338645</td><td>3.4747648085262863</td><td>0.10159930664410594</td><td>5.127417225502492</td><td>0.15886723307782732</td><td>5.5199399946376575</td><td>0.16957028237319405</td><td>9.106179154843108</td><td>0.2912562074742357</td><td>13.683663715168695</td><td>0.3360431102221078</td><td>17.998391820257527</td><td>0.5985401711559399</td><td>24.46687227379226</td><td>1.1923333949792894</td><td>43.28051183680029</td><td>2.082680118606706</td><td>82.80657419981466</td></tr>\n",
       "<tr><td>1</td><td>555969275</td><td>8687019</td><td>1.8040026</td><td>44.76484680175781</td><td>4.574854</td><td>91.49818572998048</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>5.3524575</td><td>78.49656219482422</td><td>13.370395</td><td>121.66616973876953</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>0.08941301241009876</td><td>2.0896830548623084</td><td>0.14499704199338645</td><td>3.4747648085262863</td><td>0.10159930664410594</td><td>5.127417225502492</td><td>0.15886723307782732</td><td>5.5199399946376575</td><td>0.16957028237319405</td><td>9.106179154843108</td><td>0.2912562074742357</td><td>13.683663715168695</td><td>0.3360431102221078</td><td>17.998391820257527</td><td>0.5985401711559399</td><td>24.46687227379226</td><td>1.1923333949792894</td><td>43.28051183680029</td><td>2.082680118606706</td><td>82.80657419981466</td></tr>\n",
       "<tr><td>2</td><td>555969274</td><td>8687019</td><td>1.8040026</td><td>44.76484680175781</td><td>4.574854</td><td>91.49818572998048</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>5.3524575</td><td>78.49656219482422</td><td>13.370395</td><td>121.66616973876953</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>0.08941301241009876</td><td>2.0896830548623084</td><td>0.14499704199338645</td><td>3.4747648085262863</td><td>0.10159930664410594</td><td>5.127417225502492</td><td>0.15886723307782732</td><td>5.5199399946376575</td><td>0.16957028237319405</td><td>9.106179154843108</td><td>0.2912562074742357</td><td>13.683663715168695</td><td>0.3360431102221078</td><td>17.998391820257527</td><td>0.5985401711559399</td><td>24.46687227379226</td><td>1.1923333949792894</td><td>43.28051183680029</td><td>2.082680118606706</td><td>82.80657419981466</td></tr>\n",
       "<tr><td>3</td><td>555969273</td><td>8687019</td><td>1.8040026</td><td>44.76484680175781</td><td>4.574854</td><td>91.49818572998048</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>5.3524575</td><td>78.49656219482422</td><td>13.370395</td><td>121.66616973876953</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>0.08941301241009876</td><td>2.0896830548623084</td><td>0.14499704199338645</td><td>3.4747648085262863</td><td>0.10159930664410594</td><td>5.127417225502492</td><td>0.15886723307782732</td><td>5.5199399946376575</td><td>0.16957028237319405</td><td>9.106179154843108</td><td>0.2912562074742357</td><td>13.683663715168695</td><td>0.3360431102221078</td><td>17.998391820257527</td><td>0.5985401711559399</td><td>24.46687227379226</td><td>1.1923333949792894</td><td>43.28051183680029</td><td>2.082680118606706</td><td>82.80657419981466</td></tr>\n",
       "<tr><td>4</td><td>555969262</td><td>8687019</td><td>1.8040026</td><td>44.76484680175781</td><td>4.574854</td><td>91.49818572998048</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>5.3524575</td><td>78.49656219482422</td><td>13.370395</td><td>121.66616973876953</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>0.08941301241009876</td><td>2.0896830548623084</td><td>0.14499704199338645</td><td>3.4747648085262863</td><td>0.10159930664410594</td><td>5.127417225502492</td><td>0.15886723307782732</td><td>5.5199399946376575</td><td>0.16957028237319405</td><td>9.106179154843108</td><td>0.2912562074742357</td><td>13.683663715168695</td><td>0.3360431102221078</td><td>17.998391820257527</td><td>0.5985401711559399</td><td>24.46687227379226</td><td>1.1923333949792894</td><td>43.28051183680029</td><td>2.082680118606706</td><td>82.80657419981466</td></tr>\n",
       "<tr><td>5</td><td>555969263</td><td>8687019</td><td>1.8040026</td><td>44.76484680175781</td><td>4.574854</td><td>91.49818572998048</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>5.3524575</td><td>78.49656219482422</td><td>13.370395</td><td>121.66616973876953</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>0.08941301241009876</td><td>2.0896830548623084</td><td>0.14499704199338645</td><td>3.4747648085262863</td><td>0.10159930664410594</td><td>5.127417225502492</td><td>0.15886723307782732</td><td>5.5199399946376575</td><td>0.16957028237319405</td><td>9.106179154843108</td><td>0.2912562074742357</td><td>13.683663715168695</td><td>0.3360431102221078</td><td>17.998391820257527</td><td>0.5985401711559399</td><td>24.46687227379226</td><td>1.1923333949792894</td><td>43.28051183680029</td><td>2.082680118606706</td><td>82.80657419981466</td></tr>\n",
       "<tr><td>6</td><td>555969276</td><td>8687019</td><td>1.8040026</td><td>44.76484680175781</td><td>4.574854</td><td>91.49818572998048</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>5.3524575</td><td>78.49656219482422</td><td>13.370395</td><td>121.66616973876953</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>0.08941301241009876</td><td>2.0896830548623084</td><td>0.14499704199338645</td><td>3.4747648085262863</td><td>0.10159930664410594</td><td>5.127417225502492</td><td>0.15886723307782732</td><td>5.5199399946376575</td><td>0.16957028237319405</td><td>9.106179154843108</td><td>0.2912562074742357</td><td>13.683663715168695</td><td>0.3360431102221078</td><td>17.998391820257527</td><td>0.5985401711559399</td><td>24.46687227379226</td><td>1.1923333949792894</td><td>43.28051183680029</td><td>2.082680118606706</td><td>82.80657419981466</td></tr>\n",
       "<tr><td>7</td><td>555969278</td><td>8687019</td><td>1.8040026</td><td>44.76484680175781</td><td>4.574854</td><td>91.49818572998048</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>5.3524575</td><td>78.49656219482422</td><td>13.370395</td><td>121.66616973876953</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>0.08941301241009876</td><td>2.0896830548623084</td><td>0.14499704199338645</td><td>3.4747648085262863</td><td>0.10159930664410594</td><td>5.127417225502492</td><td>0.15886723307782732</td><td>5.5199399946376575</td><td>0.16957028237319405</td><td>9.106179154843108</td><td>0.2912562074742357</td><td>13.683663715168695</td><td>0.3360431102221078</td><td>17.998391820257527</td><td>0.5985401711559399</td><td>24.46687227379226</td><td>1.1923333949792894</td><td>43.28051183680029</td><td>2.082680118606706</td><td>82.80657419981466</td></tr>\n",
       "<tr><td>8</td><td>555969279</td><td>8687019</td><td>1.8040026</td><td>44.76484680175781</td><td>4.574854</td><td>91.49818572998048</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>5.3524575</td><td>78.49656219482422</td><td>13.370395</td><td>121.66616973876953</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>0.08941301241009876</td><td>2.0896830548623084</td><td>0.14499704199338645</td><td>3.4747648085262863</td><td>0.10159930664410594</td><td>5.127417225502492</td><td>0.15886723307782732</td><td>5.5199399946376575</td><td>0.16957028237319405</td><td>9.106179154843108</td><td>0.2912562074742357</td><td>13.683663715168695</td><td>0.3360431102221078</td><td>17.998391820257527</td><td>0.5985401711559399</td><td>24.46687227379226</td><td>1.1923333949792894</td><td>43.28051183680029</td><td>2.082680118606706</td><td>82.80657419981466</td></tr>\n",
       "<tr><td>9</td><td>555969272</td><td>8687019</td><td>1.8040026</td><td>44.76484680175781</td><td>4.574854</td><td>91.49818572998048</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>5.3524575</td><td>78.49656219482422</td><td>13.370395</td><td>121.66616973876953</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>0.08941301241009876</td><td>2.0896830548623084</td><td>0.14499704199338645</td><td>3.4747648085262863</td><td>0.10159930664410594</td><td>5.127417225502492</td><td>0.15886723307782732</td><td>5.5199399946376575</td><td>0.16957028237319405</td><td>9.106179154843108</td><td>0.2912562074742357</td><td>13.683663715168695</td><td>0.3360431102221078</td><td>17.998391820257527</td><td>0.5985401711559399</td><td>24.46687227379226</td><td>1.1923333949792894</td><td>43.28051183680029</td><td>2.082680118606706</td><td>82.80657419981466</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(\"$('#table140624562142008-35025').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",
       "    $('#table140624562142008-35025').dataTable({\n",
       "        order: [],\n",
       "        pageLength: 50,\n",
       "        lengthMenu: [[10, 25, 50, 100, 500, 1000, -1], [10, 25, 50, 100, 500, 1000, 'All']],\n",
       "        pagingType: \"full_numbers\",\n",
       "        columnDefs: [{targets: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42], type: \"optionalnum\"}]\n",
       "    });\n",
       "});\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for col in master_catalogue.colnames:\n",
    "    if col.startswith(\"f_\"):\n",
    "        errcol = \"ferr{}\".format(col[1:])\n",
    "        depths = join(depths, \n",
    "                      group[\"hp_idx_O_{}\".format(str(ORDER)), errcol].groups.aggregate(np.nanmean),\n",
    "                     join_type='left')\n",
    "        depths[errcol].name = errcol + \"_mean\"\n",
    "        depths = join(depths, \n",
    "                      group[\"hp_idx_O_{}\".format(str(ORDER)), col].groups.aggregate(lambda x: np.nanpercentile(x, 90.)),\n",
    "                     join_type='left')\n",
    "        depths[col].name = col + \"_p90\"\n",
    "\n",
    "depths[:10].show_in_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## III - Save the depth map table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "depths.write(\"{}/depths_{}_{}.fits\".format(OUT_DIR, FIELD.lower(), SUFFIX), overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## IV - Overview plots\n",
    "\n",
    "### IV.a - Filters\n",
    "First we simply plot all the filters available on this field to give an overview of coverage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'decam_g',\n",
       " 'decam_i',\n",
       " 'decam_r',\n",
       " 'decam_y',\n",
       " 'decam_z',\n",
       " 'irac_i1',\n",
       " 'irac_i2',\n",
       " 'vista_h',\n",
       " 'vista_j',\n",
       " 'vista_ks'}"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tot_bands = [column[2:] for column in master_catalogue.colnames \n",
    "             if (column.startswith('f_') & ~column.startswith('f_ap_'))]\n",
    "ap_bands = [column[5:] for column in master_catalogue.colnames \n",
    "            if column.startswith('f_ap_') ]\n",
    "bands = set(tot_bands) | set(ap_bands)\n",
    "bands"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Passbands on AKARI-SEP')"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAEgCAYAAACHCxsiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXl4VOXZuO/nzEz2kIQtbAEBCbsURXBDxapVq60Wq9Zdf1axtWqx/WyttVbrglZtFbW2Wv3crfqpFFxqXRAVWkE2E3ZkJywJgeyznPf3x5mZTPaZM+ckI3nv68qVyVne9z2TmfOcZxelFBqNRqPRaLoOo6sXoNFoNBpNd0cLY41Go9FouhgtjDUajUaj6WK0MNZoNBqNpovRwlij0Wg0mi5GC2ONRqPRaLoYLYw1KYeI3C4izx9sc2k0Gk1baGGsiQsR2SQidSJSLSK7RORpEcnp6nUdLIhIdvi9fbuVfZtE5OSYvy8QkX0ickIC50f+d2Ui8kzs/y789x/aWdtYEflXeM5KEVkiImeE950oImZ47Nifo8P7PxaR+vC2vSLyfyLS3/47pdEcnGhhrEmEs5RSOcDhwJHArV28noOJc4EG4NT2hJWIXAY8CnxXKTU/gfMj/7tvAROBXyewtn8C7wOFQF/geuBAzP4dSqmcZj8LY/ZfF567GMgHHkpgbo2mW6CFsSZhlFLbgXeAcQAicoWIrBKRKhHZKCLXRI4Vkd4iMjesUVWIyAIRMcL7bhaR7eHz1ojIt2OmyRCRV8L7vhSRCTFj/kpENoT3lYrIOTH7LheRT0Xkj2FN7msROT1m/1ARmR8+932gd8y+DBF5XkTKw+v9QkQKW3sPRGR0WOurFJESEflezL5nRORREZkXnuc/IjK8g7f1MuAvwArgojbmvBp4APiOUurzRM8HUEqVAe9hCeUOEZHewFDgb0opf/jnM6XUp/Gc32zuCuB1wp8bjUbTiBbGmoQRkSLgDGBpeNNu4EygB3AF8JCIHB7edxOwDeiDpVndAigRGQlcBxyplMoFvgNsipnm+8CrQE/gReBNEfGF920ApgJ5wO+B55tpg1OANViC9j7gKRGR8L4XgSXhfXdiCbEIl4XHLAJ6ATOAulau34elLf4LS1P8GfBC+Joi/Ci8tgJgPXBX83FixhsMnAi8EP65tJXDrg2v99tKqcU2zo8cOwg4PbymeCgPH/u8iJzd1sNJPIQF+3QaPzcajSaMFsaaRHhTRCqBT4H5wN0ASql5SqkNymI+lpCaGj4nAPQHhiilAkqpBcoqiB4C0oExIuJTSm1SSm2ImWuJUuo1pVQAeBDIAI4Kz/eqUmqHUspUSr0CrAMmx5y7WSn1N6VUCPjf8PyFYaF1JPBbpVSDUuoTLKEaIYAlhA9VSoWUUkuUUrHm2AhHATnAvWFN8UNgLpYAjvB/Sqn/KqWCWAKyPU30UmCFUqoUeAkYKyITmx1zCrAIWGnz/DdFpArYivXw9Lt21hMl/L+ahvWg9ACwU0Q+EZERMYcNCFsIYn+yY/Y/HP7cLAd2AjPjmVuj6U5oYaxJhLOVUvlKqSFKqZ8opeoAROR0EVkUNkNXYmnNEfPv/Via1b/CJuxfASil1gM3ArcDu0XkZREZEDPX1sgLpZSJpV0PCM93qYgsi9z4scyevWPOLYs5tzb8Mid8/j6lVE3MsZtjXj+HZcJ9WUR2iMh9Mdp4LAOAreF1xY4zsLU1ALXh+dviUiyBjVJqB9aDzmXNjpmB5XN9MkbLT+T8s8MWiBOBUTR9v6KIyF9igrBuCY+5TSl1nVJqODAEqAGejTltR/hzEfsT+x5fH942UCl1kVJqTzvvhUbTLdHCWJMUIpKO5Qf8I1ColMoH3gYEQClVpZS6SSk1DDgLmBnxDSulXlRKHYd1g1fArJihi2LmMIBBwA4RGQL8DcvE3Ss831eR+TpgJ1DQTGsbHHkR1tx/r5QaAxyDZXpvzeS7AyiK+L5jxtkexxqaICLHACOAX4sV6VyGZWb/kYh4Yw7dDXwby+LwmI3zI9c4H3gG6//VAqXUjJggrLtb2b8VK4BM+301GgfRwliTLGlY5uY9QDAcLHVqZKeInCkih4a1uQNY5umQiIwUkZPCwrweyzcbihn3CBH5QVig3IgVKbwIyMYS3HvC419BnIJBKbUZWAz8XkTSROQ4rAeEyFqnich4EfGE1xpotqYI/8HSDv9HRHwicmJ4nJfjWUczLsOKVB6DZcr+Vvh6srB8u7Hr3wGcBJwmIg8len4MfwJOEZEOg7hEpEBEfh/+Hxphv++VWP8LjUbjEFoYa5JCKVWFleryD2AfcCEwJ+aQEcC/gWpgIfCYUupjLAF+L7AXy6TbFyu4K8JbwPnhMS8BfhDWXEuxfJcLgV3AeOCzBJZ8IZbmWIHlN401t/YDXsMSxKuwzL0tCoIopfzA97CE3V4sTfVSpdTqBNaBiGQA5wGPKKXKYn6+xjKZNzc1RzTTk4BzwwI5ofPDY+wJX/dv41imHzgE6394AMsK0QBcHnPMAGmZZzw9nvdAo9FYiBWfodFoNBqNpqvQmrFGo9FoNF2MFsYajUaj0XQxrgljEfm7iOwWka/a2C8i8rCIrBeRFTFFIjQajUaj6Va4qRk/A5zWzv7TsYJ7RgBXA4+7uBaNRqPRaFIW14RxuLpRRTuHfB94Nly1aRGQL7qbi0aj0Wi6IS2KAnQiA4mpsoRVYWkgVmGGJoQL5F8NkJ2dfcSoUaM6ZYEajUZzsLBkyZK9Sqk+SY7R1+v1PomVy65jjhLDBL4KBoNXHXHEEbub7+xKYdxaxaRW86yUUn8F/gowadIktXjx4tYO02g0Gk0biMjmjo9qH6/X+2S/fv1G9+nTZ59hGDovNgFM05Q9e/aMKSsrexKrTkETuvLJZhsxJQ8JlzvsorVoNBqNpmPG9enT54AWxIljGIbq06fPftqoGNiVwngOcGk4qvooYL9SqoWJWqPRaDQpg6EFsX3C712rctc1M7WIvITVIaa3iGzDKj3oA1BK/QWrmcAZWB19arH64Go0Go1G0+1wTRgrpX7UwX4F/NSt+TUajUZzcDNz5swBOTk5oTvuuGNXV68lWXQ0nEaj0Wg0XUxXRlNrNBqN5hvKL19bXrS2rCrLyTGL++XW3n/uhK3tHXPzzTf3e+WVV3oPGDDA36tXr8DEiRNrS0pK0mfMmDG4oqLCm5GRYT755JObJ06cWL9161bvlVdeOWTLli3pALNnz958yimn1Jx88snDd+7cmdbQ0GDMmDFj1y9+8Yu9AFlZWRMvu+yy3Z988kmPvLy80F133bXt5ptvLtqxY0farFmztlx00UX7W1tTVVWVcf755x+yfv36jBEjRtRv3bo1bfbs2VuOP/742nivXWvGGk0KYoZC7N2yqauXodGkFAsWLMh64403eq5cubJ07ty565cvX54NcNVVVw157LHHtpSUlKy6//77t1177bWDAWbMmDF46tSpVWvWrCktKSkpPfzww+sBXnjhhU0lJSWrli1bVvrEE08UlpWVeQDq6uqMadOmVZWUlKzKzs4O3XrrrQMXLFiw9tVXX11/5513DmxrXffff3+f/Pz80Nq1a0tvv/32HaWlpdmJXpvWjDWaFOTDZ/7K8n/Ni/798xffwvB4unBFGk1TOtJg3eCjjz7KOeOMMypzc3NNgFNPPbWyvr7eWLp0ac4Pf/jD4ZHj/H6/AHz++ee5r7322tcAXq+XXr16hQBmzZpVOG/evHyAsrIyX0lJSUa/fv1qfD6fOvfccw8AjB07ti49Pd1MT09XkydPrtu+fXtaW+v6/PPPc2644YbdAEceeWR9cXFx3BpxBC2MNZoUJFYQAzx04fe56ZW5XbQajSZ1EGlaL8o0TXJzc4OrV68ujef8uXPn5s6fPz938eLFq3Nzc83JkyePrKurMwC8Xq8yDMtgbBgG6enpCsDj8RAKhVorVAWAFY+cHNpMrdGkGNvXrGp1e/m2zlFEVNCk4uXV1K7c0ynzaTTxctJJJ1XPmzcvv7q6Wvbt22e8//77+VlZWeagQYP8f//73wvAEs4LFy7MBDj22GOr7r///j4AwWCQiooKo7Ky0pOXlxfKzc01ly5dmhExdSfDMcccU/3yyy8XACxZsiRj7dq1mYmOoYWxRpNivHzbL6OvT7j4yujrZ266tlPm337rZ9Qu20PFC6tRIUXoQANmbaBT5tZo2uO4446rPeeccyrGjRs39swzzxw+efLkaoCXXnpp49NPP9175MiRY0aMGDH29ddfzwd4/PHHt8yfPz+3uLh4zLhx48Z8+eWXmdOnT98fDAaluLh4zC233DJgwoQJNcmu65e//OWe8vJyb3Fx8Zi77rqr38iRI+sKCgpCiYwhTqjXnYmuTa052Hng/DMB+O4N/8OoY47njfvuYOOS/wI4bqpWQZP6dfs48O8tePLS6XleMTtuX9jqsYPunero3JrORUSWKKUmJTPG8uXLN02YMGGvU2s6WAgGg/j9fsnKylIlJSXpp556avGGDRu+ysjIaCFgly9f3nvChAmHNN+uNWONppMp/eRDHjj/TDavXNZiX0Nt40P6qGOOB+D7v/iNa2vZ9fBSyv+3lMD2aupLy9sUxBqNpm2qqqqMyZMnjxo5cuSYc845Z/hDDz20uTVB3B46gEuj6WTeefRBAF77w60AzHz5n9GglNlXnN/ieMNwJ4o6VO0nuDv+oE9lKsRoM4ZFoznoef3113v85je/GRS7raioqOH999/f8NVXX7Ue7BEnWhhrNF3MX396Bdc89kyTiMyJp5/V6rFV5XvJ7dXbkXmrPmo7IKzPTyaw57HlTbYpfwjJ0LcMTfdl+vTpB6ZPnx5X1HaiaDO1RtOJtBajUV2+F6UUD17QKICnXXZ1k2MmnfUDAPZuTbolbeO8n7XdsTR9cA/63zqlyTblNx2bW6PRNEULY42mE/jPG/9g5/o1NNRYPuFjz7+Enzz1UnT/a3f9tsnxzXMphxw2EYBPX37WkfUos+lDQb9fTWbA74+m3/8cGQ3U8uSkMejeqRScV2yd408oOFSj0SSAtjlpNC5TVb7XEqIvwyk/vg6wArUyc3IZOGos21eXsCUmmGvmy/9sMcbgsYcBcGCvM7m/++dtjL7uf8tkPD3SATDSW94SjDTLZ20GtGas0biF1ow1GpcpW782+joY8AMwoHgUAOf/7p4mx5464/oWWjEQLYVZX3XAkTXFmqgjgrgtxGfdJrRmrNG4hxbGGo3LpGfnRF+HgkEAehcNAUAMg7P/xzJRn3HdTYyfdmqb4xx65FEurrJtJKwZq4AWxprUYubMmQNuu+22wq5ehxNoM7VG0wH++iCGIXjT7KUYhcLaMEBdWLPNiYmIHn7ElISKeSjTRAxnnqP7/2ZKh8dEhbEO4NJoOiQQCODz+RI+TwtjjaYD/nbjJ+QUpHPZPcfaOj/obxTGoUCAtMxMfGntm4ZbY9Do8az/YhH1YX+zXUI1jaUtPbltNqKJEjVTa5+xJpY3f1rE7lJH+xnTd0wtZz/6jetn/PDDD/d655138hoaGoza2lpj0aJFa1s7rj20mVqjiYPqfQ22zw3GaMb+ulq8NgQxQFZeHgB1B1q9H8TNzjsXJXR8tNCH+c0qnas5+EjVfsYAX375Zc5LL730tR1BDFoz1mgAME1FxY4aeg/KabK9oS4YfX1gbx09eifcjCXqJwao3ldBTs9ettaY2cMSxrWVlfQcMKiDox0kLIxVSAtjTQwdaLBukKr9jAGmTp16oLCw0HZghRbGGg2w9F+bWfTmRk6+Ygwjp/SLbv/k5TXR15tWlnPYtMSFoBlqFMaVO3fYFsY5+QUAVO8rt3V+c9KH5cV3oEdrxprUIRX7GQNkZWUl5cfRZmpNt0eZitJPrVSfzV81FXQNNY2CNDMn8aAMADPY+LBcuWsn+/fssjVOJCrbX19n63yAbb9aEH3d+8pxcZ0TMVM3LxSi0XQ2qdrP2Am0MNZ0e/ZsreLA3noAdn3d1B/bZ0hjoJTCnjCK1YwBCocOb+PI9vGlZwAQqLfnvw7sadoUQrxxfv0jPuOQDuDSdC2p2s/YCbSZWtPt2b62EoAxxw2g9LMdhEImHo8lqBpqGwVpsMGeMAqFmrqR+o8YZWscX4YV+BVoqLd1fu2S3bbOE09EM7Z1ukbjKLNmzSqbNWtWWfPtCxYsWNd8W1FRUfCDDz7Y0Hz7J5980uJYgNra2qWR1w8++OCOtvY15/rrry8HkvIfac1Y0+3ZvHIv+YVZ9B2SCwpq9zdGPzfUBkjLsPJsAw32YjPMYFPNOMNmWpLH68PweGwL46qPbcbb6GhqjcZ1tGas6dZU7qpl+9pKpnx/GDkFlhm4el8DuT2t1/t21tKjTyZ7t1YTsFkO0mymGXu89r92vvQM28I4lsJfTIr72EafsVaNNd2b9voZJzu2Fsaabk3ppzsQQxh9TH9qD1gace3+Rp/s/j11jDyqH+XbqgnaFsZNNeO6JOpL+zIybPuMI+QcPxBfIilaUZ+x1ow13Rs3+xlrYazptjTUBvhqwXaGfas32Xnp0ZSJmrCZOhQy8dcFyczx4U33JOUz9ni90Xzjw04+zfaandCM875zSELHiyEgOppao3ET7TPWdFs2LN1DoD7ExFOtpg2ZOT7EkKhmXF8diG73pXnsm6mDAQyvj54DiwBIy0i8cEgEX3oGgSRSmwDEY+Nrb4j2GWs0LqI1Y023Zef6SjJ7pFmBW1gaYFauj5qwuToijDNy0vCmGbYDuELBEB6PhysefDzpNfsy0gk0JGemtoMYojVjjcZFtDDWdFsqd9XRs19Wk4o+WXnpUc24LlYzTvfY9hkrM4R47HV8ao4vPYP6mmpb53oLs/D1samVG6J9xhqNi2gztabbUrm7lrzCpk1nsvPSoj7jyjKrFkBGjg9vWjLCWGE41PLQMlPb8xmroBl/oY9miEdrxprU42DqZ+yqMBaR00RkjYisF5FftbJ/sIh8JCJLRWSFiJzh5no0mgj1NQHqqwPk92kqjGM141ULy8Lb0sLC2F4Al2maLerp2sWXkWHbTK0CJtgUxtpnrNG4i2tmahHxAI8CpwDbgC9EZI5SKjYs/FbgH0qpx0VkDPA2cIhba9JoIlTvs7TL3F4ZTbZn5aVRVx3ADJmkZ1qm5cycNDweIdBgTxgp00QM58zUtqOpk9GMDdFdmzRN+O1nvy1av2+9o/2MDy04tPbOY+/8xvUzPv/884dEalzv2rXLd+WVV+5+4IEHdiZy7W5qxpOB9UqpjUopP/Ay8P1mxyigR/h1HrADjaYTiPqDc5s2f8jukQYK6qoCBAMmA0fmA2B4BNNmbWalzMaewEniTUsjaFczDirbwhiP1ow1XU+q9jN+5ZVXNq9evbp0zpw56/Pz84PXXHNNwqUx3QzgGgjEPuFsA6Y0O+Z24F8i8jMgGzi5tYFE5GrgaoDBgwc7vlBN96O+KhKc1bRFaWYP6+/aKj/11QF6DrAauhgeA9OmZmhpxs4893q83hZFROJeR8hEvPYeCnQ0taY5HWmwbpDK/Yxra2tl+vTpwx966KEtxcXF/vaObQ03hXFr3/rm3+YfAc8opR4QkaOB50RknFJNS9Irpf4K/BVg0qRJ+o6gSZq2NOPMXOv7Vlflp646QEZYWFuasU1hrJwL4DK8PkLBIEqphPzQSikrGtpOjjHoaGpNypCq/YwvueSSIWeddda+s88+u8rOdblppt4GFMX8PYiWZuj/B/wDQCm1EMgAeru4Jo0GgLpqPwikZzcVxllhYVx7wE9DTSDaw1gMwbSpGZqmCeKcZgwt6113SFiQ2taMdTS1JgVI1X7G99xzT5/q6mrP3Xff3aKbVLy4KYy/AEaIyFARSQMuAOY0O2YL8G0AERmNJYz3uLgmjQawzNQZWT6MZr7crLCZet/OGpSCjLCw9iTjMzZDjmnGUWEcTMxUrcJrt1V9C0C0ZqzpelK1n/Hs2bP7rVmzJnPUqFFjRo0aNea+++7rk+gYrpmplVJBEbkOeA/wAH9XSpWIyB3AYqXUHOAm4G8i8nMsE/blSin9jde4Tl11oIWJGiAt00tappdVn1uBkBlhzTgpM7WpHPUZA4SCQVquvp01hNOyzHp7/mYMAf3V1KQAqdjPePv27Ss7WndHuFqBSyn1Nla6Uuy222JelwLHurkGjaY16qv9UUHbnNye6ZRvtx6Wh06wvCZJBXAp5/KMDa+15lAwkNB5DZusTlE1i3aSd+ohNiYWLYs1GhfR5TA13ZK66gD5fVtPkRw0sifl22vweA3SMqyviJGEz9TpaGog2gEqXrzhfOq87w6zNa8IWjPWdHt0P2ONxgGUUnzy8loOPbwvddUB+g1vXTM+dFJfln+4lVCw0UcsHiGUTGqT0wFcCQrjSI6wkWnzKy86z1ij0f2MNRoHMIOKr+Zv56v52xFDyMxuXRj3GZKL4RWKJzWWvE0mgMs0TQdTm8KacYK5xlGt3mPTXG6Asnf5Go0mDrQw1nQbYvsRK1NFc4qb4/EYXHb3saTHaJHJ+YwdDODyJKcZ260EJqJTmzQaN9HCWNNtaN7oobVo6giRFKcIhkdAgWmqFulQHWH5jJ0K4LLnM46amO2uQxf90GhcRbdQ1HQbmpuZ29KMWyMiTO00S7BqU3dtAJdKUjNG0D5jjcZFtDDWdBuam5kTEcZG2NcasuE3dieAK7HUJiLLtmumNgRdAkCTauh+xhrNN5DY6Gho30zdHE+4cpUdv7HpYGqTEfYZhwKJCePkNWNpWVleo9E4hvYZa7oNLTTjNop+tEZEM7YVxGU62CjC4wkPmaCGHkrSZ6zN1Jpm7LjlN0UN69Y52s84fcSI2gF33/WN62d8xBFHjHzkkUe2HHPMMXUAhx9++KjHH39885QpU+rivXatGWu6Dc1NzEYCdZojwthORLFpOliBKyKMEzWXq+QDuHQ0taarSdV+xpdffvneJ598sjfAihUr0v1+vyQiiEFrxppuRKxW22tgYo1aIuZdWz5jBwO4IuMoM7GuTY1marvz6trUmqZ0pMG6Qar2M7788sv33X///f0bGhq2/eUvf+l94YUX7k302rQw1nQbzLDP+IjThnDopL4JnetJwkztZDlM25pxsqlNgvYZa1KCVOxnnJuba06dOvXAiy++mD9nzpyeS5YsSbhKlzZTa7oNkXKWhxzWm96DchM610gigMvJaOqIUDcT1YxDDgRwaTO1potJ1X7GADNmzNh78803F02YMKGmsLAwwYbjWhhruhERQWrYKAmZTACXUs4HcCnTpmZssxym6K5NmhQgVfsZA0ydOrU2Ozs7dMUVVyRsogZtptZ0I5JJ74kW/bAZwIVTwtiImKnt+ox1NLXmm00q9jMG2LRpk08pJeecc86B9o5rC60Za7oNUYFkI7I5uaIfIcc046iZOlFhHLDWLV67EVw6gEujaYvZs2f3Ouqoo0bfdttt2z1h61WiaM1Y022IyBI77lsjqhnbmVc5ntqUaDQ1kYcIm8JYDNFdmzTdnvb6GV933XXlyYythbGm25CMZixJ5Bm7Ek2doM84IkhtPxQYaM1Y0+1xs5+xNlNrug2R2sp25JERPsm0JYyda6EY8RmrBM3Uyac2aTO1RuMmWhhrug2NZurO1oxDzqc22Qzgsv2NF3smeo1GEx9aGGu6DUkFcBn2NWPTNB1MbYrkGdsohyn2zdRi6DxjjcZNtDDWdBuiZmobn/pkUpuUctBMbTvPGHv2+Qi6a5NG4ypaGGu6DVEzdTKase1ymA5FU9vNM1YquW+7DuDSpCibNm3ynXbaacPa2r93717Pvffe28fu+BMnThxl99xE0MJY021IKpo6ohnbEEhuNIpItBwmpkrqgUBEd23SpCaHHHJI4N13393Y1v7y8nLPU089lVgx+hiWLl262u65iaBTmzTdhmTyjCPndHVtarBM1bbKYSZjptZdmzTN+ODZVUUV26sd7Wfcc2BO7bcvHd1mN6hrr7124JAhQ/y/+tWv9gDMnDlzQG5ubujFF1/svW7dupLFixdnXHHFFUMDgYCYpsnrr7++4de//vXArVu3po8aNWrMCSeccOC+++7bcdpppx26f/9+TzAYlNtuu23HxRdfXNnWnFlZWRM7qr7lBFoz1nQbnAjgsqMZmw7mGVtr8diLpk7GVK67NmlSgIsvvrji9ddf7xn5+6233io46qijorWlH3nkkT4/+clPdq1evbp0xYoVq4YOHep/4IEHthUVFTWsXr269IknntiWlZVlzps3b31paemq+fPnr73lllsGJRwQ6QJaM9Z0GyKC1EimNrUNzRjTuUYR1lqMhIUxyn4vY2tSK4DLyWpimm827WmwbnHsscfWlZeXezdt2uTbuXOnNy8vLzRs2DB/ZP/RRx9d88c//rH/tm3b0i644IJ948ePb2g+hmmacuONNw5atGhRjmEY7N69O23btm3ewYMHBzv3apqiNWNNtyGaJ2un6Ec0tSnxc03TdFSAGR6PLZ9xMmbqqL+56xUITTfnrLPO2vf8888XvPDCCz2nT59eEbtvxowZFW+99db6zMxM8/TTTy+eM2dOi16pTzzxRM/y8nLvypUrV61evbq0V69egUg/465Ea8aabkNjalPndm2yArjsFY9vdS02fMYqyQCu6GO7Uth6mtFoHOKSSy6p+PGPf3zIvn37vPPnz19TX18f/UCWlpamjR49umHs2LG7N27cmL5s2bLMyZMn19bU1ESF7f79+z29e/cOpKenq3/+85+5O3bsSOuaK2lKlz8NaDSdRWN95sTPjfYztl2b2kHN2I6ZOmmfcfhcHcSl6WImTZpUX1NTYxQWFvqHDBkSiN333HPP9SwuLh47atSoMevWrcu45ppryvv16xc64ogjqkeMGDH2mmuuGXTVVVdVLF++PHvcuHGjn3/++Z5Dhw6tb2++znLLxKUZi0g/YHDs8Uqpz91alEbjBklpxmJPM45osI5GUxsGZqKtHBVJCePG69d6sabrWbt2bbRZw8iRI/3r1q0rAbjnnnvK7rnnnha9jv/5z39+Hfv3smXL4kpXKisr8+Tl5XWMaiIZAAAgAElEQVSKL7lDYSwidwMXA6uByOO4As5wcV0ajeMk1SjCpmbcGDTmYACXXTN1MlK0iZlaozn42bRpk+/EE08c+dOf/nRXZ8wXj2Y8HShWSrWryms0qU7UTN2JPuNIyoSjqU12AriUQ2ZqXfhDcxBSVlbmOfHEE0c2375o0aJV/fr1S/DLZo94hPHX2PQti8hpwJ8BD/CkUureVo45D7gdS9terpS60M5cGk1HNGrGndcoQikXhLGNPOOko6m1y1hzENOvX7/Q6tWrXelTHC/xCOMqYKmI/BuI5mwppWa2d5KIeIBHgVOAbcAXIjJHKVUac8wI4NfAsUqpfSJiu2SZRtMRjUU/Ej834vK17zN2ztMqhpFwP2NlJplnbGjNWKNxk3iE8bvhn0SZDKxXSm0EEJGXge8DsU8fPwYeVUrtA1BK7bYxj0YTF5GsnKRqU9sVxo6bqW2Uw0wqtSmiGtsfQqPRtE2Hwlgp9ZSIeIFDw5vWK6XiiS4bCMRWaNkGTGl2TDGAiHyGZcq+XSnVQvCLyNXA1QCDBw+OY2qNpiVWEJM9gWSIvaIfEeHtZACXYXSFzzgyjJbGGo0bdHiHEJGpwHrgKeDvwFoROTaOsVv75jf/JnuBEcCJwI+AJ0Ukv8VJSv1VKTVJKTWpTx/bnbA03RyVRElIMQTETgBXKDKAvYlbXYsdM3VyZSxFB3BpNK4Szx3iIeAMpdSxSqljgO9iBWV1xDagKObvQcCOVo55SykVUEp9DazBEs4ajeMkK5AMkcSLfriQ2mR4jMTN1EnmGWufsSZVcbOf8Zo1a9JGjBgx1v7q4ieeO0RabNCVUmoVEE/5sC+AESIyVETSgAuAOc2OeROYBiAivbHM1m32pdRoksFqcmD/fDES7+nbmNrkZACXzWjqpBpFhH9rWaxJMdzuZ9xZxBPA9aWIPAE8F/77IqDD3o5KqaCIXAe8h+UP/rtSqkRE7gAWK6XmhPedKiKlWAVFfqmUKrdzIZrUYu/evfTu3burl9EEy0ydTLOEVAngMjq9NrUk0UJSc3Dy3uN/Ktq7dbOj/Yx7Fw2p/c61N6ZUP+NQKMQFF1wwZPHixTmFhYX+9957b31OTo7jX4R47hAzgA3A/wA3Y2mu18QzuFLqbaVUsVJquFLqrvC228KCGGUxUyk1Rik1Xin1sr3L0KQSK1asYPbs2axbt66rl9KEpM3UhjR2fop3TlfyjG2YqZ2qTa3N1JoupCv6GW/ZsiXj+uuv371+/fqSvLy80LPPPlvgxrXFE01dD9wX/tFomhAKhRCRJj7RzZs3A/DCCy9w++23R7d/tv0z1leu57Kxl3X2MoHkArjA0g4TLvrhQm1qMTwos0Wb1vbXkaRVoLEcpv0hNAcX7WmwbtEV/YwHDhzYcMwxx9QBTJw4sXbTpk3pblxbm3cIEXkp/HupiHzZ/MeNxWhSm/r6et58802WLl0a1czuvPNO3njjjSbHLVmypMW5jy57lBn/nsEfF/+RffX7OmW9zbF8xsmZau36jJ1NbUrcTE2StantNsrQaJyms/sZp6WlRT/0Ho9HBYNBV3qltKcZ/zL8+1w3JtZ8M9i8eTP5+fnk5eUxf/58li1bxrJly6irq+Poo48GYOXKlUyfPh2/3897773X5Py3N77NzQtubrJt8a7FnDLklE67hgjJNksQQzAT9Jk2Vv1ytgJXwl2bHGuhaH8IjcYJDtZ+xm0KY6XUtvDLHUC9UkqJyHBgJPCvzlicpuvYtGkTzzzzDABpaWlcffXVfPHFF4wePZoDBw6wYsUKjjzyyCbnLF++vIVWHCuIXz3rVX74zx+ysXIjDHH9ElqQrKnWMAQV6voALstMnXiecXKpTeHfWjNOiH2BIPleT6f1xO0ONO9nvGbNmqgwfe6553q++uqrvbxer+rTp0/gnnvu2VFYWBjtZ3zSSSftv/3228tOP/30Q8eNGzd67NixtR31M+4s4ommXgAcLyJ5wHysSOoLgEvdXJima3n33cZCaH6/ny+++IJgMMjxxx9PSUkJn376KaWlTeuqt2eKnXP2HIbmDaVnRk/W7lvLjuodDMgZ4Nr6WyPpwhd2oqmjAVwe2/M2x16ecXLR1Oho6oS5dd02nty2l0sH9OK+kUUdn6CJm87qZxw7NsAdd9zhWjvFeB7XDaVULVYrxdlKqbOAw9xakKbrqa+vp6ysjGnTpnHppdYz1759lp+3oKCAww8/HK/X28RXHAwG8fv9rY5357F3MjRvKAA90nrwr83/4juvf4fK+kpLS+4klFJJBXAZtszU7uQZJ+4zpvWaePHOqc3UCaGU4v92Wd+ZZ3eU05Do/0vT7YhHMzZE5EjgQsL1obHyhjUHKRGNt3///ng81r+6srISn89Heno6GRkZDB48mI0bGwXp66+/HhXYzemT2Vj8pkd6j+jry9+9nA37N/DlxV/i8/jcuJQmKJWc71Yk8dQmV/oZG0bCtamVQ7WptZk6PjbUNVARCDEmO4PSmnoW76/h2IIWsUSaFKGtfsYff/zxms7qZxzPHWIm8HtgnlLqKxEZhmW61hyEVFVVMWeOVSht8ODBUdNzbW0tmZmZUWEW2T5x4kQAVq1aRVlZGQUFBVx00UVNxjx6wNHR1z3SGoXxhv0bAFi5d6VLV9OUZAO4DI9gJugzdqMcptiNpnbCTK2FcVz8p9JKff3jyCIEWFhZ0/4Jmi4l0s+4+U9nCWKIQxgrpT5USp2hlLpLrDvxLqXUTzphbRo32Lkc/vNX8Ndaf3/1OhzYGd29e3djF8uMjIyoEGloaMDna9ReTzjhBKZOncqJJ57YZPiqqiq83kaDyxOnPIERYxvOS89rsaTNBzYndUnxknwFLknYZxrRYJ3MM7a6NtmIpnaiUYT2GcfFwspq+qR5mdgji7E5mSyqrO7qJWlSnA7N1CLyLHAdEAQWA71F5F6l1INuL07jIKvfhi/+Bhs/tjrNB2qgtgI+fxgGHQlX/Zuqqip27LB6eVx77bUAUTN1IBBoIoyLioooKipqIRSCwaZ58z6jqfm5IL1l8ZqymhbxFq6QdJ6xJJ5nHE1tclozTrhrk0NFP7TrMy4WVlZzdH4OIsLR+dk8v6Mcv2mS5uDnQHNwEc8nY7xS6gBwNlZK0yDgcjcXpXGQUBBWzYXXroQNH8Lwk6ztmz+3BDHAtsUopXjggQf44IMPAMjPtzpZxppXY4VxBMMwyMtr1HavuuqqJvu90vR5r2dGT5qztapzCvkkn2ecRG1qB1NbbJXDVEk2itDR1HGzxx9ge0OAiblW2eYpeTnUmYqVVXVdvDJNKhNX1yYR8QLfB95USvnRz8ffHD57CF65CIJ1cMU7cPHr4MuGnSsajxGDUH1Tn1Z6ulXxrSNhDHDFFVdEXw8Y0DRdKdub3eTv/IwW7aqpqK9osc0NnMgztl0O00nN2GMnmtohM7X2GXfIu3v3AzAl3/rsT86zfn+xX/uNNW0Tzx3iSWALUADMF5HBgHaAfBNY9iJ88kfr9Sl3wJBjrNeGB6rDpuETbwEV4oXn/zd6Wqx2G48wjmjRzY8HKMwqbPL39qrtLc6vCXTOTSr5PGMbZmqVGuUwkzdTR4Sx/SG6C8sO1FLg9UQ1477pPkZkpfO3bXsI6IcZxzlY+hnH0yjiIeChyN8ishU4yc1FaZIkUAebPoM3Lb8vP/kP9B3VuD9WIA0+ChP4evue6KZY7TYeYQyWJt3Q0LJ5gTRLbr1i3BU89dVT0b+H5Q2jOtA5z3bJNoqwNOPEzjFdaRSReGpT8uUww7+1mbpDVlbXMT43s8mD348H9eF/1m7j7o07+N2hA7twdc5R8draokBZjaMtFH39smt7nluckN8q3n7GkbaLqUp7jSJ+FP59fewP8DOstoqaVMQMwTPfhRemN26LFcQA9ZYZje8+AN50GmhsQjJ+/Pg2tbj2hPH111/Pz372s5Y7mt2789Lz+PSCTxlZYKX0FeUWURuobft6HKQrGkVEfKxOB3B1ts9YdGpTm+z1B/nduu3sbggQMBWrq+sZl9NURl0yoBdTC3J4bkc59SGTHfV+zv5yHaXV2o+cCNdee+3AWC135syZA373u98VRrTXxYsXZ4wfP370qFGjxhQXF49ZuXJl+k033TQo0s/4mmuuGbR//37j6KOPLh4zZszo4uLiMc8//3xL31krhOtej5k/f35Wa/Mke23tacaRsFdb6r2mi1j+EmyPqQ99w/K2j83sCR4f9THC+LDD2i6u1p6pNTs7m+zs7BbbjVYkQF56HrO/PZtFOxdRWl7K0t1L216jgygzKbep1SgimKh52PkKXIYNn3GyJvqoVq014xY8u2MvT2zbg4nigv698CvF+NzMJseICD8p6suPVmxk/r4qFlZWs2h/DY9t2c3sMV1QqN0BEtVgneDiiy+uuPHGGwdHtNy33nqrYPbs2ZtffPHF3tDYz/jaa6+tqK+vl2AwyAMPPLDtzDPPzFy9enUpWJkh8+bNW9+zZ09z586d3ilTpoy68MILK9u7vy1fvjz9ggsuGP7UU099fcwxx9RddtllRc3nSZb2GkU8Fv7926Rn0XQOpgmfPgT9J1jBWv5ayGnnWSotGzxp+PHFDNH0Jh97A7fj92wr+rZfdj/OPvRsNlZuxB9qvYym06gk6zMbBgRTIIDLMAzMBFObcKwCl/0hDlYigVkrq+oYl2NpuuNyMlscd1xBLj28Bv8uP8Cn+6oAWF2TEj0KvjF0RT/jiooK79lnn33oq6++umHSpEn18c6TKB3eIURksIjcJyL/EJH/i/wkO7HGBTZ+COXr4eifWYK2PUEM4MsETxqBGGHcq1evNg+3I4yD/vafGH0eH37T3ykpM10bwOVcBdlEK3AppZKvTa3N1G2yNixQv6quY31tPV6BQzJbWi19hjA+J4t39+7n6zo/PbwGpdV11AQ7rcjTQUFn9zPOzc0N9e/f3//xxx/nJDJPosRTm3oO8CzwPvq5ODXZ/Dl8dDdsWgBZvWHM9+M7z5cNHl9UM77sssvo3bt3m4fbEWStBXXFkmakYSqToAriE3frUycbwGVV4ErsnKilwdE8Y0uwK9OMT+MOrzkpU7muwNUqIaXY2RCgf7qPnQ0BPig/wOCMdHxtvNfjcjL5LFyN6/rBhfxh406WVdXqutUJ0Nn9jH0+n3r33Xc3TJs2bUROTo45Y8aMitbm+d73vleVzHXFI4z9utpWCrP8ZXjjmsa/x/8QvHH2ys7MB096VDNOS2t5XrJm6lAH5tQ0jzVnIBRoUa3LcZIM4DIMG7WpTXdqU4Ml6D3xjBvRZnVqU7tUVZWyes2tjB83m4yM+Np7lvuDmMC3e/bg+Z3llNbUc3KvHm0ePzbGl/yDwgL+sHEnq2vqtTBOgK7oZ9yjRw/zvffeW3/iiScW5+TkmCUlJRnN50n2uuIRxo+IyK3Ae0BUzVFKrWj7FE2nsPFjSxAPmgw5fWH1XCgcE//5WT3DPmPrY9CaMI7FSYESISKM/SE/WT5HsyRaoJKre2GzNrXzFbgahXEITzxfYZW8MI5YFA7WClymGeCLxdNRys+69fcyftzDcZ23J2C5YY4ryOG1XRXUm4phrZioI0R8yYMyfPRP95Hn9UTN3Jr46Yp+xr179w599dVXqyL7WpsnGeIRxsXAVcDpND4XK+B4JxeiSRAzBO/cbL3+3iPwadh44YlTKwbItALmI5pxe6lLYE+gdHTzjgpj0/0gLmUqxGv/gcJebepwowgHfcaRh6J4/cYqainXZmqlFLv3vENuzmiysoZGt9fUrEMpP4aRwe7d7xAK1eHxtAzCas7uhgAAA9J9FGdlsKK6jlE5GW0ePyo7g18P7c8P+hUgIhRnZbBGC2MN8Qnj84BDlFJJR4tpHMI04fkfwJ7VcNj5Vh7xt39nCejR30t4OH/PUVDhjpm6I9IMa86GkPsfL6UUyVyCYZB4OUwX8owNT6PPOM5FhE9MZtKDw0y9efMTbNh4PwCDBl7Ctu3PcewxC6ipWQ/A8OG/YN26P/DF4h8wZfJcRNp/iNodDlDsm+7j/lFF/HvvAaYXtmyGEsEQ4YZDGqvSjchO5/3yA8leliZJUqGfcTzCeAWQS4yJWtPFrHzVMlF/+3dwzPXWtryBcO5T7Z7WFsEeQ6BCNWl92BpOpjZFiPUZu41lpk7CVOuxoxk7n2cc6zOOCwd8xlHF+BseTR0RxADbtj8HwJatf8djZCLiYeCAC1m37g/U1Kzl669nM2zYDe2Ot9tvfW77+LwckpnOhNzEXC3DMtPZ4w9yIBiih9c564kmMSL9jLtyDfHcXXsBq0Vknk5tSgG+eBLeuBr6fwuOvRE88TxPNWPKDJh2a/TPENZNwNuBsHVTM441UzeEGlzp5KSSbJZgSOKNIiIC040ArnjbKDa2cezeRT+CQSvYtWfPqdFtImmUl8+npmYdmZmH4PGkc+SkNwDYf2Ap1dVr2L//yzbH3BcIkW4I2TYF6aiwD3lFVedUodOkLvHcye9yfRWa9lEKXr4Q1rzduO20e7Btcz19VpM/QwiCiUGI5h+JWE3SDZ+xz2P5qWMLf7zy1Qv87dMHeeaitxiePzzhOdteiwOpTQlXoXS+NnUktSl+zTj82wmf8TfYTF1buwmAgQN+RP9+P6CmdgMeTxYbNtyH37+HggKrkUqPHodRWPg9Kvf9h//89wwAvn3ShlbHrAqFyPXY12gn52WTbgivlu1je32AR7bs4h8ThjMgI4HYD81BQTzC+HOgXimlRGQ4MBKrr7Gms9jxZVNBfMLNjR2YHCCEgQcTQgHwth0J6oaZOt1jzRfrMw4+9yp/nRtizaSFDJ/ipDDugqIfbrRQTDSAywGf8cEQTV1XtxmArKxDyMmx3IOV+63SscFgFbk5o6PH5vX4Frt2zYn+3dZnpypJ83IPr4ez+uTzSlkFr5RZ9Ss+rqjiwgFtF9/RHJzE8/VcAGSKSH9gPnAt8HdXV6Vpylf/B4bPin7OK4Jptzg6fEgZeAiB2b7f1s3Uplif8bBVViOL6hXO1qxOtja1vX7GLgRwRXzG8ZbEjKzBEc34myWMlVJ88OFwPvhwOLVhYZyZWRTd3yN3fPR1z16NCSI9ekxoMo7ZRrT/gWCI3CQi9AF+Naw/A9J9nBsO/NIlMuNn5syZA2677bbCjo90l08++STr8ssvL+r4yLaJRzM2lFK1InIlMFspda+ILEtmUk0C+Gth6fMw8jQ488/taq52iQrjUMvSlW7Vpo7Qms+4PsPSNL5c9i7n8UDCc7a3lmT8pklpxk5W4Eo0mtrRoh/fLGFcV7cp+nr//iWkpfXB42kMsjKMNAYNugQz1ECP3MYmKTkxWjKAadbh8bT87lUFzaQDrwZlpPHlMVbL3MUHatjjdz+YUeMsxx9/fO3xxx+flOM/LmEsIkcCFwJXh7fpsL/O4qvXoL4SplwL2e6YrkJKLDN1B5qxkwIlgtewPoKxPuMA1kPB4N2KORvm8L3hiadrtUbyRT8SjyY2XTRTx+szdiSA6xsav1VTsy76urx8Pnl5h7c4ZmTx7S22eTzpFBQczb59CwEImfW0loV/IBSiT5pzD8h9fD72dFDPPVV48803i3bv3u1opZ6+ffvWnn322e1Gb9588839Xnnlld4DBgzw9+rVKzBx4sTakpKS9BkzZgyuqKjwZmRkmE8++eTmiRMn1m/dutV75ZVXDtmyZUs6wOzZszefcsopNSeffPLwnTt3pjU0NBgzZszY9Ytf/GIvQFZW1sTLLrts9yeffNIjLy8vdNddd227+eabi3bs2JE2a9asLRdddNH+1tY0d+7c3AceeKDwo48+Wm/32uO5Q8wEfg/MU0p9JSLDsEzXGrcJBWD+/dBvvKM+4hbTIGHN2HkzdbwBXLcvvJ01FWvYVraWQ0sqAZi2UvF86fMJz9nOYpL2GSesGLrUKMIaOt484/DvZHzGIpZA/oZpxrW1TQovkZkRvyVxwmF/Y/iwmwAwQ633Ha4Ohsh1MCWpZ5qHisA3Qxh3BQsWLMh64403eq5cubJ07ty565cvX54NcNVVVw157LHHtpSUlKy6//77t1177bWDAWbMmDF46tSpVWvWrCktKSkpPfzww+sBXnjhhU0lJSWrli1bVvrEE08UlpWVeQDq6uqMadOmVZWUlKzKzs4O3XrrrQMXLFiw9tVXX11/5513DnTz2jrUjJVSHwIfxvy9EfiJm4vShPl6PuzfAqe/6GijgeZYmnHrwthtM3WkOUSVv4o/LPoDt8xc3GT/kFBcfb/jXEvyPmOVYG1q041+xjHlMONbRKRTRJJrMCTlVeOysrdYv+E+pkyeh8+XT23tJny+XgSDB1AqQFbWsLjH8ngyycw6BGjfZ9wjSZ9xLD19XlZUtS74U42ONFg3+Oijj3LOOOOMytzcXBPg1FNPrayvrzeWLl2a88Mf/jAa7en3+wXg888/z33ttde+BvB6vfTq1SsEMGvWrMJ58+blA5SVlflKSkoy+vXrV+Pz+dS55557AGDs2LF16enpZnp6upo8eXLd9u3bXQ1x71AYi8ihWNrxIbHHK6VOjePc04A/Y5m1n1RK3dvGcecCrwJHKqUWt3ZMt6RspfV7yLGuThNSxGWmdiOAK6IZA9SHWgauHPXqKjjbmbnMJPOM7bVQdD6ASxL0GTtipgaQxFO7OpuNXz9MQ0MZ27e/yCGH/ITauk1kZQ0lLa03e/a8S58+pyQ0niHWLU+pltqqqRTVIdNRzbjA62Wf1ozbpbl1yzRNcnNzg/EW7Zg7d27u/PnzcxcvXrw6NzfXnDx58shIC0Wv16si9znDMEhPT1cAHo+HUCjknkZEfIar14BVwB+A38b8tItYdeQexappPQb4kYi06GIgIrnA9cB/4l92N2H3Kugx0Oqu5CLtacaxuJFnHPEZQ9OI6sy/hmtt760g0MFDQgKLwUhCIBmGYCaoGboSwGXbTJ3cGsQg5c3URrjz14aND/DhRyOprPwvWVlDGTvmfg4//OVoSlO8SNhyY7YijKtDJgqSyjNuToHPQ72pqA2l+FNPF3HSSSdVz5s3L7+6ulr27dtnvP/++/lZWVnmoEGD/H//+98LwBLOCxcuzAQ49thjq+6///4+AMFgkIqKCqOystKTl5cXys3NNZcuXZoRMXV3NfEIY1Mp9YhS6nOl1H8iP3GcNxlYr5TaqJTyAy8DrTXavRO4D9Dx/M2pq4TstvsLO0XIpM3UJrdrU3ulURhXNTTW6B0y9TT2fedIhu0wKavc5shcyaY2iQAqsSCuxjxj533GiaY2JSuMkdQ3Uzc07Im+jmizBflT8HiyKMg/MuHxJKIZt/LdqApa77+TZSx7+qz5tHbcOscdd1ztOeecUzFu3LixZ5555vDJkydXA7z00ksbn3766d4jR44cM2LEiLGvv/56PsDjjz++Zf78+bnFxcVjxo0bN+bLL7/MnD59+v5gMCjFxcVjbrnllgETJkyo6dqrsognmvotEbkaeIOmLRQ7qm4+EIj1KWwDpsQeICITgSKl1FwR+UVbA4Xnvxpg8ODBcSz5ICHkh1bSKRyfJmqmbv9p3M0ALoDymt0AzD9jIKNFSJ88iYz3vmBXyWKKpg5ta4jE1pKMZuyxzjWVwkN84ygXymEmWoEraqZOVjm3YabvTEKhBoLBSoYOvRGvN4e8HhPYWfYWffueYXtMMdo2Ux8IC2NHzdQ+a6x9gSADdRWuVpk1a1bZrFmzWrQvXLBgwbrm24qKioIffPBBi/Jpn3zySYtjAWpra6PFDR588MEdbe1rzplnnll15plnVnW09vaIRxhfFf4da5pWQEdSsbWvfvSbLFZ9wIeAyztagFLqr8BfASZNmpS6dwOnCfldyStuTjBkkkYIVPualiuacYyZ2hOWLftD1oNq72GjqQXKvi6BqT9Mei4n+hlDWLjFef91s5+xijeAy4F+xta8xHyDUw+/33qYy0jvz4AB5wK0msqUCEY7ZurKsDDOd1QYW9+HikCnNArSpBDxRFPbrSqyDYg9dxAQ+6SRC4wDPg7fqPoBc0TkezqIK0ywATJ6uD5NSCnLTN1KdI7btaljzdSe8P1nbF+rKlL/YYexATiwdWPC87a6FjPJ1CaJCOME5nSjhWKiPmOnoqlT3Ezd0LALgPT0vo6N2dxMvb3ez6a6Bo4tyGV/WGDm+ZzXjHV6U2ry+uuv9/jNb34zKHZbUVFRw/vvv9968fIEiKvlj4iMwgrCinbNVkq92MFpXwAjRGQosB24AKtwSOT8/UDUISoiHwO/0II4hs4yU5sqbKZ2XjPuSBh7Ynyp3rBsOXqI1VXH16cPAa8gZXtaO9XWWpKtwAWJ+4ydbBJhrSPiM47XTB05z4HUphSOK2oIa8ZpTgrjZmbqH5ds4ssDtWw8/jAqg9Y2JzXjnuE2pvuCWjNORaZPn35g+vTprrRa7PAuISK3YpmI/4IVGf0n4NyOzlPWp/c64D2saOx/KKVKROQOEXGmpNJBglKK15Zso6ah2dNwyA+e1ur+OEso1LZmHIsbZupYeoY9Lr4+VqlZMQyqCtJJ29Nq0ZuEcSLPGEioPrUyQ45qxRBbDjPRAK7k5hVJTZ9xKFRLbe0m/A2WME5Pc04YNzdTf3nAqni4rrae/WGBmeegMM6P8RlruhfxfD3PB6YBO5VSlwATiFOjVkq9rZQqVkoNV0rdFd52m1JqTivHntgdtWKlFHNX7OQXry7nbwuamWODDZ3iM7Y049Z9xm4X/YjFF77/GJlRAwwNBdlkVjpTBMGJrk2QmGZsmqbjDzFRzThOe7lyyGeMQUqaqdesvYOFi75Nbd0mRHz4fAWOjd3cTO0Lf37W1tRTGQghOBtNnRb+3973dYv4JM1BTjxCtU4pFRKRYDgnuAyIv4yNpk0agiHO+8tClm+zNL/stGb/jlCgEwr74awAACAASURBVM3UHWvGbtSmjsUTuc/H5G36e2aTs6bSmQmS1oyt34lpxqbzmnEkmjreXFQnujZB2EydesJ4585XAaisXEJ6Wh9ng+WaFf3I83rYGwiyoz7A/nD7RMPh78WQjDQ21/upDATJ98Wl92gOAuK5SywVkXystomLgf8CX7q6qm7Aks0VXPfi0qggBsjNaC6MGzrHTG2abaY2daZmbER8m57G98FfkEP+gZAjfXSTDuAybARwmWbUrOwUCUdTR9abbDS1SCoqxtEuTNXVpaSlO9tNL1JEJCKMQ+E3YKffEsZOmqgj/PwQ6xr2a79xt6Ldu6tYd67blVKVSqlHge8C1yilLu2U1R2k/OOLrUx/fCHvl+7ioimD+d6EAQAEmmsdwc5JbQqFzDbN1LG47TP2RDS4mBtcsGcPfCEw9yfvNzYVyTVLiPqM45fGphs+44QrcDlopk5Bzdjrbcw4yMjo7+jYEc3YVEGUUlSFC62UNfipDIQcDd6KkBV+eNNVuFpn4sSJo9wa+4UXXsi75ZZb+gG88847OWPGjBnt9XqPePrpp53zfbRBuzYQpZQSkbnAEeG/bbeH0lhUNwSZ9e5qBuRl8OhFh/OtonwqawPMWb6DYPMvX8gPHncT/wOBAA2BIFnUuRLAlYhG64lqcI03uFDvPGudu3bjyU+yLKipkjIpGjY1Y6eFsdgt+pHsMiQ1zdSx/YkzMgY4Onasz9ivFMHw5e9sCOATcTStKUJm+HNWl8BDX3di6dKlq5tvCwaDeL3Jm/TDLRL3AwwbNsz/9NNPb7r33nudNbe0QTyr/6+IHK6U0qbpJDlQH+CGl5ZSXuPnjZ8cw8TB1sOWN1zZKdi8I1CowXVhXFdnBUdlU9dqapPbecaxRM3UMdqG2csSwHW7tpMxsjjh+ZuuJTmBlGoBXCrRcphJ+jZT1UytzMbI44x0l4QxJv6Y/3tZQ4AeXg/90p13I/nCn7NgCj74xFK66uaimuq1jvYzzs4prh0zela73aCysrIm1tbWLp07d27unXfe2b9v376B0tLSrA0bNpS01af4tdde63HbbbcNDIVC0rNnz+DChQvXtjb2ww8/3Gvx4sXZzz777JaRI0f6wX2LYIQ2hbGIeMPpSccBPxaRDUANVmUtpZRKrrRNN+PzDXu55Kn/opTiD2ePiwpiAJ/H+mcHYp+EQ0FLBXPZTB3RrgzMLk9tysAHNDQJ4KJ3TwDqd+5o/aQEUEnmNkUEeSIPGK4EcHkS1IwdahSRqmZqUzXWjXZeMw4/+KgggfAbmef1sMdv/V3gQoCVN/wZDabeW51yrFixInvp0qUlo0aN8oPVp7iwsDBUXV0tEydOHHPxxRfvM01TrrvuukM+/vjj1aNGjfLv2rXLeXOGA7T3SfovcDiONbDrvizZXMHFT/4HU8FNpxRz8VFDmuz3Gq1oxqFw/1SXA7gahbFyxWeciOD67eRbqX71t0iMucno3QuAhl3Jp3okn2ccKbbR1dHU9ipwOVGbOhVVY6WCeL09MM0ABQXHODq21XwOlDKjmurgjDRWVtdREXAngKtRGKfeex1LRxpsZ3DYYYfVRAQxtN6neNeuXd7JkydXRY4rLCxMyci49oSxACilki7z1Z2Zu2IHP39lGf3zMvnfK4/k0L65LY7xRIVxrGYc7snhcmpTtFwjZqs3Wre7NsWSjpdqGvv1Angzs6jKgLQ9u5MeP/lo6vA4CdwkzVDItWjquAPJHOralKpFP0wzQGHh9xg18veOjx2tnqZCUc14cKYljAH6u2Gm/oYI41QgKysr+iVoq09xsvUFOov2hHEfEZnZ1k6l1IMurOeg4t+lu7juxaUM75PNazOOoSC7df+viODzSNNo6khvX5c144hgMVAdlsP02BAq8QiMFZeuIKiC1M57LzJRdJ/X8LIvF/J370147lgaa0R3bm1qN3zGiVbgUg4JY0szTm4IN1AqGE1Bcp6IZhyKCsdxOZnMC1eFG52d6fyMWhjboq0+xdOmTau56aabhqxevTotYqZORe24vbuEB8jBaujQ2o+mDUKm4i/zN3DN80sAeO7/TWlTEEcIhBSPfxxjhIgEpRjuJv1HuwpZjXrbPdZOtGIg0LIPbHNEBJ/hQ4Wsa441U3vEQ0WOYO5NVhhH5rI/hr1ymC7Wpk7UTJ10P2NS0mesVCAaaOU0jWbqRs14SGajtWp0Tkar5yWDN/xv0sI4MdrqUzxgwIDgww8/vOmcc845dOTIkWPOOeecuIpWzZ8/P6uwsPCwt99+u+DnP//5kEMPPXSsm+tv7xO8Uyl1h5uTH6y89N8t3PvOarLTPDxz5WQG5Nt4eu4kYdxopm7dZxxr3nFLGEcJRwfHmqk9hod9uaB2lic8dyzR60zGTO2x2SjCrWjqRFObPMm2UBRHiq84jVJBDNeEsQCGJYzD76NXhPcnFbO+tsGdAK7wQ1MgBd/rVCDSV7h5D+HMzEzVVp/i884778B5553XYYOH66+/vhwoBzjhhBNqd+3atcKhZXdIhz5jTWL86d9r+dO/13H0sF68+OMpcd/8p43sw95qf+OGiAnScDfwL5FoareFsYpUHIo1U4uXihygYh8qFGoiqBMi2rnI3ulANEc5Ec3YFTN1tBxmnJa2kINm6kBqCQilTJQKIeKeO0fEg8KMaqo+gfG5WYzPdTSrJ0okgCuBOEHNQUB7d9dvd9oqDhJe+u8W/vRv68Hsju+PTUgL8xhCKPYmH+17564wbqIZO+QzjtWeEqlWFTVTN9OMK3MEQiahykq8vXrFPV4sphOacSSWJ9GuTQ4HcCUaTe2UZhxOakxuDIdR4bQm93zGVhBXc83YTXQAl/v8+c9/7vX44483KeZx5JFHVj/33HNbumpNbQpjpVRFZy7km84jH6zjoX+v5fDB+bx89dGkeRPThgyRqMAAOk0zbmqmbr82dbyCTNmN8olEkzcL4NofVkCC5eW2hXFUgCZxH021oh9xP+g4pBmLIQlZBToDM+zOcctnHB1bmVGzcZrbZWEjwjjF3uuDiRtuuKH8hhtuSM735TC6JYgDPPT+Wv78wTqOL+7DH889LGFBDK1pxmFh7HDwT3Oa5hk7VH4v9jISKZDRRgDX/uyw2a4iiefDiDxKQiDZDuBKEZ9x0mZqSb1o6kgDB3FZMzZjin54XXbgac24e6KFcZIs2ljOnz+wTNOPXjiR3Ax7NwXDkGhHGKALNOOOfcZxj2lbM24ZwOU1vByI0Yxtr8kRM7Wd2tShaC1pp2gsPpJAOUxDDsoWipHqW25qxuABFeszdrmVqOgAru5I5xTdPEh5v3QXF/5tEQBv/fRY24IYrC+g2apm3FkBXK37jG3Vo44RxglpxpEArmaacUQYh8rta8bRZRwsZmoRVJxPBSqkkteKCSvGKSaMlRn2GbsawGWgiPEZO/BetodObeqeaM3YJhU1fn7+yjIO7ZvDwz+ayKh+PTo+qR08Xa4Zd+wzdpu2AriqM0EZQrCiazVju2Zqw+P818wwjITKYSYdvAUpWQ6z0Uztrs84tuiH25qxjqbunmjN2CaPfLiOGn+QRy88PGlBDOEArth7a1dEU3dQmzpeBhYNZEv2libjx0VEM47RJL3iRYlg5uUmpxlH3s7O1oxdKIcJVnpTvGZqFTId0YwtM3XywzhJRBi7rhmrUDSry+1oaiMqjLU0bo3O6md8++23Fw4fPnxscXHxmKOPPrp47dq1rrbQ08I4QZRS/PTFL3n6s01cNGUwIwqdKUbmMWgawGW2FExu4EYAl2EYfNH3CyDRAK4QeDxNC42ENR4zP4dgEgFcjpTDDP8rujqAy1qLkVAFLic0Y5EEH646AdN032cs4mmqGXeSmTrFPAIpQ1v9jJ3goosu2n/33XeXARxxxBG1y5YtW7V27drSs88+e9/Pf/7zQY5M0gbaTJ0gn6zby7wVOxmYn8lvzhjj2LgtzNSd7jM2O8wzjhf7AVzBFjm5nvD1h/JyCCUTwOWAZhwxUyfUKMIFnzFYwjjuaGqHfMapGMDVKXnG4QAuf/j97qwArlT3Gd+4akvR6pp6RyufjMrOqP3T6MEp0c/4rLPOilb3Ou6446pfeeUVe3mVcaKFcQI0BEPc8/YqstM8vH3DVDLTnBOURvMArk6uTW200bXJDrHCKjHN2GwSvAWWzxggmJdN8GsHfMZJacZ2oqnd0YwT9hk7IYwl9czUZsRn7KKZmrCZOtJf2G1hHPm0hFItjywF6ax+xk888USfk08+eb/zV9CIFsYJ8Mf31rC6rIqHfzSRvExnv/xdFcAVCvsdPYQc8xlHSfCepVrRjL1h82MgP4tQeasPs/GNHQ3gsj1EYwBXov2MXcgVF08CPmOnzNQpWJtadYqZ2tukUYTbZmoRwZN6RogWdKTBdgad0c/4scce67l8+fKsJ554Yo2zq2+KFsZx8p+N5Tz56dec/a0BnDm+v+PjG9JW0Y/UqU0d95gx4yQawNXCTB3RjHtkY9bUYNbX8//bO/MwOc76zn9+VdUz03NodIzuw7IuB9mSDRYyNnbAAbywARxsvBiWPLsPOIJkYUk4doHdsDxk9+EIuwkkzm4M4QHyPIthHa8P1uZ0jB3ZkQ9syxKybCHbWJJ1ja65u6vq3T+qqrs10zPTR1V1yf37/KPpOt737Z5Rf+t3vlZX/TvllN3UMVjG9bqpE0ngqsMyjs1NTeYUIp1s6qC0yU2pHSaAjWTeTZ0Fkt7P+I477uj76le/uvTBBx/cm8/nE/2FaAJXDQxPuHzytqdYOa+b//auTU11cZoOe3KrweiLNjXLmJbHjI3nTXFTRwlchf5g56tGu3DFu59xnb2pW5zAFZdlHGZwNT9OjPip1BnbmIp2mEm7qQFs0WzqeplpP+MdO3b0PfPMMx0Atbqpt2/fnv/oRz963p133rlv+fLl8WSIzYBaxjXwp3fs4uDJMb7/ocvp6UzmI5s+gSulbGoh/g5cdWbfVnNTRwlchTmBGLuDJ8gtW9bIooIlNeOmbmALxaQSuCzbTt8yFtq2HWZlNnUqlrFI1pwQmef6668/fcsttyzcsGHDxrVr145X28/Y930WLFhQfOihh6putVjJpz71qZWjo6P2DTfcsDYcp3DfffftS2r9KsazsPvQaf7vEwe58bUree3q+YnNM6XOOO0tFIXYYsaVDTaadVNHlvHEnMA17TXY+COedpjBv9kpbao9Zozd/BrEkux14IqyqRPfKKLcgSvpmDEEYqxu6uqktZ/xQw891HiSSgOom3oWfrz7CAA3XXV+ovPYFi0pbSq5qS2Jb6OIBonqjCuJLOOJvkCM3RMnGxs7emspJ3D5mk2dKOVs6gTFmLJlLJRLj5LEFog5nVLJOGoZz8DpsSJ///ALXLyin3WL4mnuMR325ASuVljGdew9POOYjSZwVXNTh+9/rC9oftNKy9iyow0a6myHmYRlLHXEjOPswJUxa63Um9pKsDlS2PSjYAwdKVjFELmps/VZv5I4p/YzVuDb21/g5GiRr7/3gsTnqux7bFVaqWlZxhKfZVwZM67rPs8HZ5rSpk4L6erCbbAlZjmBq6HbgXLMuC7LOKl2mLYdeBJqWgRILqbSpow1TE5j1yYRG4OP65tU4sWg2dRJk8X9jNVNPQ3PHx/h6/c9x1XrB7hq/cLE57Mn96NN3TKOrzd1400/XMSu3vTDMz7O/PkNd+GKo7SpLMa1P7RkoR1mbJZxWPyapVpj46fgpg4t46IxqWRSQ5RNncpU9eL7vp/e7jGvMMLPrup/XBXjafjyvc/gWMIXrr0olfkiy7jkqk4pm/qsmHHMpU1Sr2nselMSjSyxsMTCNS72ggUN96eOq7RJLKm/6Udi7TBr/H15MTX9iMbIUBJXOtnU5d7UqVnGIlktbdp17NixfhXk+vF9X44dO9YP7Kp2PlE3tYi8FfgaQRnrN40xX5p0/uPATYALHAM+YIx5Mck11cKJkQL37T3K+7au4vyBnlTmtCeLcYqWcSAyVvwJXPWWNvkeUuX92mLj+R7O/Pm4x441tJSSZdzQ3WUsuz4xTrK0qZ4OXLHsZxyVdnkm6ehJzfgmpf2MI8s4tZhxNuuMXde96fDhw988fPjwRagxVy8+sMt13ZuqnUxMjEXEBm4G3gIcAB4VkbuMMZXp5U8AW4wxoyLyh8BXgPcktaZa+fZDL1Bwfd7/ulWpzTnVTZ1Ob2rP87BtO4hNx+Sm9hsVdc+vWoLjWA6uH1jG43sb60gXPRRYTVqI9YpxYBnHr1y24+B7NfYhiKvpR/S7yZD/NJV2mOFGEcUUY8aOSCazqS+99NKjwDtbvY5XIkk+2WwF9hlj9htjCsCtwLWVFxhj/tEYMxq+/Gcg0S2qauHUaIHvPvwCb9m4OPEM6krsigQuINV2mJZlBe7wuBK4GqwzntEyNh7O/Hl4g4MNxSwjAW1WlAIxridm7CViGduOg1fjtnHxW8bZqW8q7Wec5ENr1A4zxZixlV03tZIQSZpdy4HKRuIHgMtmuP6DwL0Jrqcmbnv8AKdGi3zsTetTnXeqmzqddpglMbbsaUubli1bxoIFte8eVtkOs77Splks4/mLMMUi/vAwdl99D0rRQ06zrUwt2zq7BG22eb1kLGPLyVEYH6/tYi+mOmMrezHjNHZtqtwoIr1s6my6qZXkSFKMq/3VVv3rEpH3A1uAN0xzfhuwDWDVquRcx0XP57bHD7B+US8XLe9PbJ5qlBK4TPoJXIGbenrLeNu2bY0NXm9p02yW8YKgA5o3OFi3GEclOc2KsW0LvltjFrMxGOOX6pPjpF7LWOLowBWOYdzsiITxC4jYieyMFVGKGfvp1Rk7IlmKBigpkOQ3/QFgZcXrFcChyReJyJuB/wS80xgzUW0gY8wtxpgtxpgtCxcmU2bk+4YPfPtRnjk8xEd+Z10ic8xEFDMuGacpJnCV3dStLW2azjK2rUCM7fmBdd5IRnVkGcfjpq7tPUUxXduJ32qznRx+jWIc9KZufk5poDd30vjGTTReDK0pbbK0HWbbkaQYPwqsF5HzRaQDuBG4q/ICEXk18LcEQnw0wbXMyt07D/Hgc8f57L/8Ld55cQMbETRJOTcm3ZhxyTK27NhixsfGgoznemO701nGjgRu6sgydo8fr3tNcbqpaxXjyHJNoulHYBkXa7o2yH6O4b969CCTsZhxki5qqBBjP71saucc2M9YiZfExNgEmRUfAX4M7AF+YIzZLSJfEJEoG+/PgV7g/4jIkyJy1zTDJc4//PIg5w/0cNOVa5pqDNEolkxK4GqFZRxTnfH9L90PQNEvxmYZB9nUA8FlDTT+KLmp07SM3bCG24nfcrPqcFOTQGlTVjB+GpaxU6ozTq/phyZwtRuJ/hUbY+4B7pl07HMVP785yflr5ZnDZ9i+7zg3XXl+InsV18KUBK60LWPis4wvXXwp33j6G+Rz+bruqylmLIJ7rHHLuNlEpnqyqSPL1bLj/29mO07Nbmrj+2WrtqlJs1fa5JsCVoINP6BsGReMYU4CmfHVsAR1U7cZbV+0bYzhv9y5m74uhz9849qWrcOenMCVdjZ1jKVNUZ1xI1soTu5NDUE2ted7iONgz5/foJs6bPuZops6asphJWAZ206uLss4jmzqaIxMlTalZhm7uKm6qSVrG2QpCdP2Yvz/nn6ZHc+f4JPXXMDc7gR3fpmFKW7qFLOpg9KmGOuMo3aYMdUZO5aDG5awOAMDDYlxXG5q2xa8OmPGrXRTG2OCvj9x9aYmY25q4ybafQvOtoxTc1Ojbup2o63FuOj5fPGeZ9i4dA7v3Zpet61qTLWM3cBFnfB/ft/3y6VNMcWMI8s435Pn9OnTtd84Xcw4bIcJ4Cxc2FBLzFa4qaNs6mQsYwe/lgSumLLIzxojQ2LsmyKScJe6ypixk5JlrG7q9qOtxfiHOw9x8NQYn7hmQ0kMW0VkGZdixl4R7GSf+KHSTR1fzDgSYyfnlDaiqIVpY8ZhAhc0bhn7pTrj5v7kLaf+bGo7sWxqd1bPQ8mKjSWBK6wzzlCar/GLKVnGLgXf0JFiO8wMfcxKCrStGB8fnuBL9z7D+kW9XH3BolYvp6IdZnjAdyHhxBSY3PQj5l2b6o0ZT9eBSxw8E1nGgRjXXzbVimzqqLQpmZgxMGt/alMM4/fOK7O0yTduupaxZlMrCdG2YvxXP3+Ok6NF/vLGS1qWQV3JlDpj34UEvsQnk0Sdcdy9qaN2mBC4qSkW8U6dqmtNXkwduCyrnmzqZGPGlXNMRyTGVkfz1nkmS5tSqzMOLeMUd23KUKMzJQXaVoyfPHCaLefN48Jl6ba9nI6qbuqEn/hhkmUcc8w4Lss4Km2CUIyh7rixia0DVz3Z1ElaxjWKcSH43KQjDss4e6VNxi+mkk0NJnXL2K/ePVh5hdKWYjxe9Njz8hkuXDan1UspUXJTlyzjYipuatd1cRwnjBnH85/fpyzG9TBTzLgygQvqF2M/Lss4I9nU0Ziz1RqXxTgGyziDpU2+KWIlLMZWWOtf8P0U9zNWN3W70ZZivOflMxRcny2r57d6KSXsKZZxym5qkdhixpRyhqwEYsYNinFMMWMnZ+EVa82mjuqMk9i1KbKMZ86oNoXwwSgGyzirbmrLSrYkUUIxdg0pduBSN3W70ZZi/OJgsIXy2oW9LV5JGauUwJWuZZxEzDhyUyP1bhQxvWVcDDeRL4nx0cbc1M1axk7Owi3W9tBStoyT2Siico7p8IvxWcZksbTJLyKJi7GDgdTrjH21jNuKthTjXQdPk7OFlfPra9eYJFXrjFOIGZfd1DHGjEM3db2WcdC2cfp2mABWdzdWT0/Dbupm64ztDhuvUGedcUKlTVCHmzoXg5s6F22hmB03tTHplDb54Vdlem7qTD3zKCnQlmL82Isnec2qeXQm4D5slKkJXG4qdcZlN3Xrs6kDy7iKmzpsh1l63UDjD98zWE7zX6SBZezX9L5KpU2tzKYuRNnUMbipQ0GPBD4L+H4BSbw3tYMbtvFPc6MIbfrRXrSlGI9MuMxrYevLalRP4ErZMo65zhjqc1NPZxk7VjlmDI2Jsef62DFsI+iEolZL3Ljc9CPJbOrZYsYxJ3A5UiqXygJB04/k9zP2WiDGmk3dXrSlGI8VPbrjiKHFSDmBKzyQQgcu3/fL7TATiBlbVr0JXF7VfXdtKXfgAnAWLarfMnZ97BgaXzihdejWIcZJbRQB4BXSS+CCwDrOkhgH7TBTtIxTrTNWMW4n2lOMCx5dGRPjyDtbclOnEDOOWlWW64xjdlNTb9MPv+ouVZXtMKFsGdcztuf68bipQ1Fza4gbu4UJAHKdnU3PO5lcZxcAxYnxGa/zY4wZA1g5qzRmFgg2ikg+m9oj+PzS3c84lamUjNCWYjxa8OiO6cspLqa6qdMT47KbOuYtFK0GYsazNP2AQIzN2Bj+yEgdQ5uYLONQjGvIqC5OBGLsJCHGXbWJsSl64FixbKEIgbs7U5axn45lXCSYIzXLWLOp2462E2NjTMbd1OltFHGWZTx2Eo7tgf33Nz3u3fvvBgLLuFbMDPs356zc2Qlci+ovb4rLTW1HbupaLWMRnFz8llvZMp6Y8TpT8GNJ3oqQnJUpMTamkEo2dYHgd9gdQ95BLThaZ9x2tJ0Yj4dfJPmO5JOj6sFqQQcuN4xpOo4Dh54IDj58c9PjPn7kcaDO0qbwwWDamLE5200N4B49WvOaPNdgx+mmrtEyznV01t2JrBYi13dxfBbLuODFU2McEohxNtzUxpjULOMCweedb3LXr1rRDlztR9uJ8Vj4RdKVy9Zbn2IZ+166MeMIaf6L+90b3g1AV66rZjGeyTKubIcJ0LF6NQAT+/bVvCbP87FiyaYO1ldLrbE7MZGIixrqcVP7pfrgOJAOu5QU1mqM8QCTSjb1RMpibAn41Nk0RzmnyZYipUAkxplzU1vV3NTJfsm4Z+23G1pvVcSwXmyxmds5F6jjy6Qey3jxYuyBAcaf3lnzmtyCV7Jqm6EcM67NTZ1E8hZAriO0jGcT4yQs44wkcBkTZJIn34HLphi6qbtSc1NHTYBSmU7JAO0nxuEXibqpJyVwdXQHB6X5Pwnf+Fhi1dX0YybL2LEcfOOftRtUftMmxnY+XfOaihM+uRhEqZxNXYObenwcpyMZMRbLwunonDVm7Be82MqaIFsxYz9skZp0zNiyOpgIxTif4kYRAJ7WGrcNbSfG46FlnM9aNvWUOuPkO3Cd5abuCPt0xxDf9I2PLXZdYhztT2x1T21Raoeu88qM6vzFmyk8/zze6dM1je8WPJzOGMS4jjrjYmGilGiVBLnOzhpixn4sexlHSIedmdKmsmWctBh3lmPGKVnG0SxqGbcPbSfGYxkV4ygU9cV79wQ/+G4sLuOZOCuBqyvcTrI41vS4nvGwLZtcLkehUKjtnlOBqDqLFk0554Sx88q4cX7zZgDGnt5V0/jFCS8WyzgXCnpxfOY2lBDEjJNyUwPkuvIUxmf+fcXtpra6HMzY7O89DXw/+NtK3jLuLGVTd6UUMy67qVWN24X2E+OSmzpbbz1qz5mLnrxTdFPbtg2dfcHBieHmx/U9bLGZM2cOY2NjJdGf+abgGqnSraokxhWWcdemTSDC2M6nalpTXJZxV0/wOxkbnrnzFUBhfDyxBK5gLb1MjMz8+zKFeBO4rLyDKfqZ2CzC9wOvgG0nu+GLZXWWE7hSsoxtFeO2I1uKlALlbOpsWcY522Ldol5etybcY9mdACe5L3KAiTDe2NnZCcteExzMz216XM94WGKVsrT9Gjp7mfDBYLpdm4CzunDZfX10rFnD+FO1JXEFlnEcdcYWHV12TWI8MTJMV09y23R29fYwPkvjE1OM2TLuDneLqsEzkDSeF3gFLDu5UACAbXeV3NRpWcYZ3K1SSZi2E+OsxowBHEsoeiZoSzkxJyYolwAAEKhJREFUBJ1zEp3vzJkzAHR1dcHrPxYcXH1V0+NGMWMr/OKqSYyLoWVcJYPctqaKMQSu6rGdO2eNSxvf4Bb8WCxjgK6+DsaHZne/jw2foau3L5Y5q9HZPbNlbHyDP+Zi5eNLVozG8jPgqva8YF9y2+pOdB7L6mSMPDYm/QQutYzbhrYT47KbOntinLMtXM+HwhBgoKs/0fl+9KMfAaFlHNU0x9AS0zNe3WJcclPnqohxlQQuCJK4vJMnKR44MOPQE6FwdMYkSvne3KyWse97TIyMJCvGPb2MzyDG/pgLBqze+MId0pUhMU7NTd3FCL3Msd1EGrhUQ7Op24/2E+MsW8a24PoGxsMM4a5kLeOIQIzDzyOGbRTH3DFydq4+y7hUZ1y9HSacncAFFUlcs7iqx0PhzPfFU4+a7+uYVYwnQvdxvjdJN3Uv48ND0573hwPr3Y5RjEtu6pHZ3fRJ40eWcQox41F66LXSewCJNqQo+CrG7ULbinHWYsYAOcui6PkwHriPk7aMIyzLKtcXx2AZHxo+xIreFSUxjhLFZqJ48FDww0xuanP2l2Hnhg1IV9esSVyRcHbFJEr53tysburhE4MA9MybH8uc1eie049XLDIxWj1u7A6GlmNffLkHzvwgPuuemLmkKg2imLFtJ+2m7mCEHvokvQeQnjBRbMRrfaKckg5tJ8bjBQ8R6Ixh04C4cWzB9QwcfCw4kHDM+CxiFONxd5zuXHddlrF38gQAzsKBKedKbupJlrE4Dl0XXjhrEtdYKJz5uMQ4tIz9GayWM8eDvtlzBqaWasVF/+IlAJw6/HLV85EYO4vjEyurJ4d02niDrRdj1wseQqyELWMRi9Myj3lWeu+5N/x+GnazUdOtJE/2FClhxooe+ZydWuynHhzbougbuDtMpupMLt4I0N/fz8UXXxy8iHpSx7Cn8Zg7Rt7Jl8S4llpjf2QEyeVwFi+ecq7DDtzLE97UblP5zZsZ37MHM8Mccbup+xfl8T3D0OD0Nb7Hf/MiAPOXr4hlzmrMW7ocgJMvH6x63h0cQzrtkms5DkQEZyBP8dhobGM2ilsMwjk5p/kKgNkYZIABu/YtO5ulLwzXDKll3Da0rRhnkZwlQQJXxIK1ic7neV55k4gYLONfn/o1m76ziZMTJ+l2uktjHz9+fNZ73VOnsOfOrfqQ1NcRPJScKZyZci5/8WZMocD43r3Tjj14aBgnZ9HdH48Yz1/aA8CJl6cXpCP79zF3yVI6u3timbMa85YsA2DPP91f9fzY7kFyi7tjf/DMLemheGik5ZsYFIsnsawu7IRLm0Y9n2F6GeBUovNU0hNZxjWEeJRXBu0nxgU/k/FiCNzU4oYW3hs/C/l5ic53lhhbzYvxz178WenndfPWsTi0cmuJGbsvH8ZesKDquRnFuIYkrsP7z7Bo9RzsmBo2zF/Wg1jCkf3Tt+I8vP85Fp+/Lpb5psPpCB4u9v/y0SnnisdG8YcK+KPxJx11rOzFHyninZy5L3bSFIqDdOSSi8lH7BsN3NMDproHIgl6w/+XIxlorqKkQ6JiLCJvFZG9IrJPRD5d5XyniHw/PL9DRFYnuR6AofFi5nZsiujK2ayYCC28/uWJzlUoFJiYmAhqjCPEaiqbeufxsiBetOAiOkKxOHbs2Iz3ucePM7J9O91bX1v1/JyOIHZ+ZmKqGDtLl2IvnH4HJ7fgcfw3QyxZE1/8vaPLYcmaObzw9GBV6/DU4ZcZOn6MxWuSFWOAS9/+LgCOvfj8WceHHwiEY+474/eudK4J3MInbn0m9rHrYWzsAF355MIAEU8NBR6Q1f6exOeK6AsfHM9ozLhtSEyMRcQGbgbeBmwE3isiGydd9kHgpDFmHfAXwJeTWg9AwfV58qVTrB5IznXYMEOHed/Qt/lfE5/BdA/AujcnNpXv+zz11FP4vs+qVavKJ3I9cKy+L9hT46f44f4f8tLQSzxw4IHS8dX9q+nuDhKH9u7dO61Ls/DCCxz5ylcA6L3yyqrXLOkOEpU+99DnppwTEfKbL65qGRtj2HHXfnzfsGpjdau7UTa8djGDB4e553+Wd47yfY9fP76Dv/vYHwBwwRXNN1CZjQvf8CYAvvsfPsrzTzxGYXSMwsFhRh49DEDn+vjjqblF3eBYFH4zxInv78U9NY5JuQTH9ycYHn6Gnp7kH3h+NniGAWuMeYXdqbnm5zg2cxyL/WOt9T4o6ZHkPoJbgX3GmP0AInIrcC3wq4prrgU+H/58G/DXIiImgb/4/73jN3zx3j0MjbtcvibeL+amuPfTsOsfYOQolwG/8s9j92v+hhv6liQy3bPPPssdd9zB6Ogoy5cvZ82aNeWTl/87+MWX4JFvwNY/mHWshw89zLafbpty/M5r78QSCytnsXLlSl566SUeeeQRLrvsMiCoKX7+hhvwBk/gHjkCQPfWrfRccUXVeXIVu1fd/tztXLf+urPO5zdvZvjnP8c7fRq7vx/jGx649Vl2P3gQY2DD1sUsvyBel/+rrlzGL773LC/sPM63PvUgG7YuYN+Omxl8KUjcuuDyqxLNpI5YuGo16y+7gud2PMTtX/o8m+b9NhvnXg5Az+uWJpaouOjDmzn6108y+sRRRp84Su/rlzH3HcnmOETs3PlhhoZ243nDLBy4JrZxP7jreZ4aGmW+49Dr2Ay5QcuNXcNjfHDBGN7xIXY88jbOO+/DLF3ye7HNWw0R4aLebn5w+ASDRRfPGAq+4QMrFvLmBSlWWSipkaQYLwdeqnh9ALhsumuMMa6InAYWAGdl/IjINmAbcLYlVwcr5uV564VL+N3NS3nDhoUNjZEI/cthwzUwcAGsvZo7n3C4Zn0yQgwwb9481q1bx9q1a9m4cWM5Zgxw5Z/A8Wdh/prpB6igN9fLhQsuZPfgbi5acBFvO/9tXLP6Gpb0lNd/3XXX8bWvfS3YGSpEbJuOFSuRdevIb9pM79VX07FiZrf8g+95kC8+8kWW9S6bci5/ySV0rFtL8cgR7P5+ChMeB589Sa7L4Yrr1vKq10+9p1ls22Lb197ArgcOcujZkwyf9LAsi3lLl/G669/Lxquujn3O6XjHn3yGv7npfYwPDzHiBa78/EULmPd7yVmNHSv6WPzHr+HIX/4SAImp1Wgt+KZAd/f5rFnzcRYsiM/7cHFfN3nL4mTRY9jzWNyZY8TzeOeiuXx2wwW8/OKHGBr+FSLpvNdPrl7Cl59/mWdHxumwhJxYTMRQ7aBkE0nK7SIiNwD/whhzU/j694GtxpiPVlyzO7zmQPj61+E1g9ONu2XLFvPYY48lsmbllYNX9LEcSa2Ezfd8rJR29FGURhCRx40xW1q9DqU6SX57HABWVrxeARya7hoRcYB+4ESCa1LaBDtnpVpLrkKsKEozJPkN8iiwXkTOF5EO4EbgrknX3AX8m/DndwP3JREvVhRFUZQsk1jMOIwBfwT4MWAD3zLG7BaRLwCPGWPuAv4O+HsR2UdgEd+Y1HoURVEUJaskmcCFMeYe4J5Jxz5X8fM4cEOSa1AURVGUrKOBLkVRFEVpMSrGiqIoitJiVIwVRVEUpcWoGCuKoihKi0ms6UdSiMgx4MWYhusHpt96Jz6SmCeuMZsdp9H7671vgEmd2ZS6SevvPSmysP401xD3XOcZYzLUflCp5JwT4zgRkVuMMVObK58D88Q1ZrPjNHp/vfeJyGPaPag50vp7T4osrD/NNWTh/Srp0e5u6rvP4XniGrPZcRq9P63PXilzrn/mWVh/mmvIwvtVUqKtLWPl3EEtY0VRXsm0u2WsnDvc0uoFKIqiJIVaxoqiKIrSYtQyVhRFUZQWo2KsKIqiKC1GxVhRFEVRWoyKsXLOIiI9IvK4iLy91WtRFEVpBhVjJTOIyLdE5KiI7Jp0/K0isldE9onIpytO/UfgB+muUlEUJX40m1rJDCLy28Aw8F1jzEXhMRt4FngLcAB4FHgvsIygRWYXcNwY88OWLFpRFCUGnFYvQFEijDEPiMjqSYe3AvuMMfsBRORW4FqgF+gBNgJjInKPMcZPcbmKoiixoWKsZJ3lwEsVrw8AlxljPgIgIv+WwDJWIVYU5ZxFxVjJOlLlWCm2Yoz5dnpLURRFSQZN4FKyzgFgZcXrFcChFq1FURQlEVSMlazzKLBeRM4XkQ7gRuCuFq9JURQlVlSMlcwgIt8DHgYuEJEDIvJBY4wLfAT4MbAH+IExZncr16koihI3WtqkKIqiKC1GLWNFURRFaTEqxoqiKIrSYlSMFUVRFKXFqBgriqIoSotRMVYURVGUFqNirCiKoigtRsVYURRFUVqMirGiKIqitBgVY0V5hSAifyUivxSR17Z6LYqi1IeKsaK8AhCRHmAR8CHg7S1ejqIodaJirJxziMhfiMgfV7z+sYh8s+L1fxeRj8c853DM480VkT+qeL1aRHbVeG9eRH4hInZ0zBgzAiwF7ge+LiIdIvKAiOg2qYpyDqBirJyLPARcASAiFjAAXFhx/gpgewvWVQ9zgT+a9arqfAC43RjjRQdEZAHQDQwBnjGmAPwceE+zC1UUJXlUjJVzke2EYkwgwruAIRGZJyKdwKuAJ0TkDhF5XER2i8i26GYR+fIkq/TzIvIJEXm/iDwiIk+KyN9WWp4V11a9JrRs94jIN8L5fiIi+fDcn4rIMyLyUxH5noh8EvgSsDYc58/D4e1q91fhXwN3Tjr2n4GvAruBjeGxO8JrFUXJOCrGyjmHMeYQ4IrIKgJRfhjYAVwObAF2hpbhB4wxl4bH/n1oPQLcytkW478CHguPvd4YcwngMUnIRORVs1yzHrjZGHMhcAq4XkS2ANcDrwauC9cC8Gng18aYS4wxn5ru/snvPdzTeY0x5oWKY6vDz+H7BNtMRl6CXYAmcynKOYDGk5Rzlcg6vgL4H8Dy8OfTBG5sCAT4XeHPKwnEbtAY84SILBKRZcBC4CSwCbgUeFREAPLA0UlzvmmWa543xjwZ/vw4sJrAhX6nMWYMQETunuE9Vbt/MgMEQl3JfwW+YIwxIlISY2OMJyIFEekzxgzNMK+iKC1GxVg5V4nixpsILMCXgE8AZ4BvicgbgTcDlxtjRkXkfqCr4v7bgHcDSwgsZQG+Y4z5zAxzznbNRMXPHoFYSx3vqdr9kxmj4n2IyCUEFveVInJzeO7pius7gfE61qAoSgtQN7VyrrKdoITnhDHGM8acIEiKupzAbd0PnAyF+LeA1026/1bgRgJBvo0g2endIrIIQETmi8h5k+6p5ZrJ/BPwDhHpEpFe4HfD40NAX71v2hhzkiC2HAnyl4F3GGNWG2NWAxcTWsahW/6YMaZY7zyKoqSLirFyrvI0gcv2nycdO22MOQ78CHBEZCfwZ5Ouwxizm0AMDxpjXjbG/IogCeon4T0/JSgVqrxn1msmY4x5FLgLeAq4nSA2fdoYMwhsF5FdFQlctfITAkv4d4AeY8zPK+Y7AvSIyHzgauCeOsdWFKUFiDGm1WtQlFc0ItJrjBkWkW7gAWCbMeaXTYz3auDjxpjfn+W624HPGGP2NjqXoijpoDFjRUmeW0RkI0E89zvNCDFAmID2jyJiV9YaVxJmXd+hQqwo5wZqGSuKoihKi9GYsaIoiqK0GBVjRVEURWkxKsaKoiiK0mJUjBVFURSlxagYK4qiKEqLUTFWFEVRlBajYqwoiqIoLUbFWFEURVFazP8HCsieEnksLh0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for b in bands:\n",
    "    plt.plot(Table(data = parse_single_table(FILTERS_DIR + b + '.xml').array.data)['Wavelength']\n",
    "            ,Table(data = parse_single_table(FILTERS_DIR + b + '.xml').array.data)['Transmission']\n",
    "            , label=b)\n",
    "plt.xlabel('Wavelength ($\\AA$)')\n",
    "plt.ylabel('Transmission')\n",
    "plt.xscale('log')\n",
    "plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
    "plt.title('Passbands on {}'.format(FIELD))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### IV.a - Depth overview\n",
    "Then we plot the mean depths available across the area a given band is available"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "vista_j: mean flux error: 3.786496639251709, 3sigma in AB mag (Aperture): 21.261602923956396\n",
      "vista_h: mean flux error: 5.977168083190918, 3sigma in AB mag (Aperture): 20.76595819133626\n",
      "vista_ks: mean flux error: 6.494208812713623, 3sigma in AB mag (Aperture): 20.67588124149986\n",
      "irac_i1: mean flux error: 1.3417919874191284, 3sigma in AB mag (Aperture): 22.38798387786226\n",
      "irac_i2: mean flux error: 1.0861581563949585, 3sigma in AB mag (Aperture): 22.61746419360326\n",
      "decam_g: mean flux error: 0.1085532901417139, 3sigma in AB mag (Aperture): 25.118089385646975\n",
      "decam_r: mean flux error: 0.12857820618118967, 3sigma in AB mag (Aperture): 24.934278456842115\n",
      "decam_i: mean flux error: 0.21076208454265227, 3sigma in AB mag (Aperture): 24.39771564992062\n",
      "decam_z: mean flux error: 0.39704004934858017, 3sigma in AB mag (Aperture): 23.71011107278023\n",
      "decam_y: mean flux error: 1.4089685117666761, 3sigma in AB mag (Aperture): 22.334943644655247\n",
      "vista_j: mean flux error: 8.843372344970703, 3sigma in AB mag (Total): 20.340652085418803\n",
      "vista_h: mean flux error: 14.800461769104004, 3sigma in AB mag (Total): 19.78150869960445\n",
      "vista_ks: mean flux error: 17.018739700317383, 3sigma in AB mag (Total): 19.62987837357435\n",
      "irac_i1: mean flux error: 2.187960386276245, 3sigma in AB mag (Total): 21.85709822647123\n",
      "irac_i2: mean flux error: 1.7050905227661133, 3sigma in AB mag (Total): 22.127828261926147\n",
      "decam_g: mean flux error: 0.17185977202471453, 3sigma in AB mag (Total): 24.61926128480865\n",
      "decam_r: mean flux error: 0.20815990768290052, 3sigma in AB mag (Total): 24.411204146671473\n",
      "decam_i: mean flux error: 0.3731635695088045, 3sigma in AB mag (Total): 23.777448766495205\n",
      "decam_z: mean flux error: 0.7493719872714292, 3sigma in AB mag (Total): 23.020453227111126\n",
      "decam_y: mean flux error: 2.7687564515644643, 3sigma in AB mag (Total): 21.601484974372745\n"
     ]
    }
   ],
   "source": [
    "average_depths = []\n",
    "for b in ap_bands:\n",
    "    \n",
    "    mean_err = np.nanmean(depths['ferr_ap_{}_mean'.format(b)])\n",
    "    print(\"{}: mean flux error: {}, 3sigma in AB mag (Aperture): {}\".format(b, mean_err, flux_to_mag(3.0*mean_err*1.e-6)[0]))\n",
    "    average_depths += [('ap_' + b, flux_to_mag(1.0*mean_err*1.e-6)[0], \n",
    "                        flux_to_mag(3.0*mean_err*1.e-6)[0], \n",
    "                        flux_to_mag(5.0*mean_err*1.e-6)[0])]\n",
    "    \n",
    "for b in tot_bands:\n",
    "    \n",
    "    mean_err = np.nanmean(depths['ferr_{}_mean'.format(b)])\n",
    "    print(\"{}: mean flux error: {}, 3sigma in AB mag (Total): {}\".format(b, mean_err, flux_to_mag(3.0*mean_err*1.e-6)[0]))\n",
    "    average_depths += [(b, flux_to_mag(1.0*mean_err*1.e-6)[0], \n",
    "                        flux_to_mag(3.0*mean_err*1.e-6)[0], \n",
    "                        flux_to_mag(5.0*mean_err*1.e-6)[0])]\n",
    "    \n",
    "average_depths = np.array(average_depths,  dtype=[('band', \"<U16\"), ('1s', float), ('3s', float), ('5s', float)])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def FWHM(X,Y):\n",
    "    \n",
    "    half_max = max(Y) / 2.\n",
    "    #find when function crosses line half_max (when sign of diff flips)\n",
    "    #take the 'derivative' of signum(half_max - Y[])\n",
    "    d = half_max - Y\n",
    "    #plot(X,d) #if you are interested\n",
    "    #find the left and right most indexes\n",
    "    low_end = X[np.where(d < 0)[0][0]]\n",
    "    high_end = X[np.where(d < 0)[0][-1]]\n",
    "    return low_end, high_end, (high_end - low_end)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "data = []\n",
    "\n",
    "for b in ap_bands:\n",
    "    data += [['ap_' + b, Table(data = parse_single_table(FILTERS_DIR + b +'.xml').array.data)]]\n",
    "    \n",
    "for b in tot_bands:\n",
    "    data += [[b, Table(data = parse_single_table(FILTERS_DIR + b +'.xml').array.data)]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ap_vista_j (11670.0, 13380.0, 1710.0)\n",
      "ap_vista_h (15000.0, 17900.0, 2900.0)\n",
      "ap_vista_ks (19930.0, 23010.0, 3080.0)\n",
      "ap_irac_i1 (31754.0, 39164.8, 7410.801)\n",
      "ap_irac_i2 (39980.1, 50052.3, 10072.199)\n",
      "ap_decam_g (4180.0, 5470.0, 1290.0)\n",
      "ap_decam_r (5680.0, 7150.0, 1470.0)\n",
      "ap_decam_i (7090.0, 8560.0, 1470.0)\n",
      "ap_decam_z (8490.0, 9960.0, 1470.0)\n",
      "ap_decam_y (9510.0, 10170.0, 660.0)\n",
      "vista_j (11670.0, 13380.0, 1710.0)\n",
      "vista_h (15000.0, 17900.0, 2900.0)\n",
      "vista_ks (19930.0, 23010.0, 3080.0)\n",
      "irac_i1 (31754.0, 39164.8, 7410.801)\n",
      "irac_i2 (39980.1, 50052.3, 10072.199)\n",
      "decam_g (4180.0, 5470.0, 1290.0)\n",
      "decam_r (5680.0, 7150.0, 1470.0)\n",
      "decam_i (7090.0, 8560.0, 1470.0)\n",
      "decam_z (8490.0, 9960.0, 1470.0)\n",
      "decam_y (9510.0, 10170.0, 660.0)\n"
     ]
    }
   ],
   "source": [
    "wav_range = []\n",
    "for dat in data:\n",
    "    print(dat[0], FWHM(np.array(dat[1]['Wavelength']), np.array(dat[1]['Transmission'])))\n",
    "    wav_range += [[dat[0], FWHM(np.array(dat[1]['Wavelength']), np.array(dat[1]['Transmission']))]]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Depths on AKARI-SEP')"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfIAAAFfCAYAAAC1E07gAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xlc1HX+B/DXGxAQGQ8OQRHFg+FQMzxo06y0S0tbjVrb1FIrtGM3XXM1K3PXrWTNLHM1LTs0K0tLTbNjf6W5aW3igXKqhSd4och9zef3xwzuRAMzwHxhvuPr+XjwWPgen8/nO2y++XzmO9+XKKVARERE+uTR3AMgIiKihmMhJyIi0jEWciIiIh1jISciItIxFnIiIiIdYyEnIiLSMRZy0iURyRaRm5t7HEREzY2FnBrNUlRLRKRARC6KyE4RmSIiTvn/l4i8IyL/cEZbWhORriJiEpGlNvYpEelh9fOTIpIjIj3rcX6RiBSKyEkReVlEPK32bxORh+oY23WW302+iOSJyPciMsCyb4KIVFnatv7qaNlf/TsuFJHTIvK2iPg3/JUiImdhISdnGamUMgDoAmA+gJkAVjbvkJrF/QAuALhXRHxqO0hEngEwFcANSqnUepzfRynlD+AGAGMATHJkUCLSGsBmAK8BCAAQBuBvAMqsDtullPKv8XXKav9IS999AQwA8IwjfRORtljIyamUUvlKqU0wF5kHRKQXAIiIj4i8JCLHLDO610WkpWXfjSJyQkRmi8g5y+xvrGVfIoCxAP5qmQ1+ZtXd1SKSYplhrhURX8s5QSKy2bI6kCciO2pbHRCRgSLyk6WNn0RkoNW+bSIyzzJzLRCRr0QkyM5LcD/MBa4CwMha+vwHgIcAXK+Uyqrv+QCglDoM4HsAV9sZTzWj5bwPlFJVSqkSpdRXSqkUB8+37vskgK0AetX3XCJyPhZy0oRS6r8ATgAYbNmUBHMxuRpAD5hnhHOsTgkFEGTZ/gCAFSISpZRaAWANgH9aZojWxe0PAIYB6ArgKgATLNunW/oOBhACYDaA3zyLWEQCAGwBsBhAIICXAWwRkUCrw+4DMBFAewDeAJ6s7ZpFZDCATgA+BPARzEW5pvkw/5FzvVLq5wacX31sNMyv7eHajqkhC0CViLwrIsNFpJ2D59nqOxzA7QD2NrQNInIeFnLS0ikAASIiAB4GME0plaeUKgDwAoB7axz/rFKqTCm1HeYC+wc77S9WSp1SSuUB+Az/m51WAOgAoItSqkIptUPZDhW4A8AhpdRqpVSlUuoDABn49Uz4baVUllKqBObiWtcM+AEAW5VSFwC8D2C4iLSvccytAL5QSh1r4Pl7RKQIQDqAbQB+8166LUqpSwCug/kPmjcAnBWRTSISYnXY7yyrGNVfR2o0s0FELgL4D4DtMP8OiaiZsZCTlsIA5ME8M/YDkFxdJAB8Ydle7YJSqsjq56MAOtppP9fq+2IA1TdfLYB5pvqViPwsIrNqOb+jpR9rRy3jttfHr1jeJrgH5tUDKKV2ATgG84ze2r0A7haRvzXw/L6WMYwBcA2AVrWMZ6vVDWtjLW2mK6UmKKU6wbws3hHAK1an/aCUamv11b1Gs6Ms27sopR61/HFDRM2MhZw0YbkbOgzm2ds5ACUAeloViTaWG6eqtRMR66LUGeYZPWBjWbwuSqkCpdR0pVQ3mGfXfxGRm2wcegrmm/OsdQZwsj79WYwG0BrAUhHJFZFcmK+/5vJ4FoCbATxa4w8MR8+HMvsIwC78+u0J62OGW92wtsbG/gwA74DvcxPpHgs5OZWItBaRETC/z/ueUuqAUsoE83LuouqlYhEJE5Hbapz+NxHxtrxXPALAx5btpwF0q8cYRohID8uS/iUAVZavmj4HYBSR+0TES0TGAIiF+e7u+noAwFsAesO8/H41gEEw35DX2/pAy13qNwOYISJT63u+lfkAEkUk1N7gRCRaRKaLSCfLz+EA/gjgh3pdJRG5HBZycpbPRKQAwHEAT8N849hEq/0zYV7u/kFELgH4N4Aoq/25MH/s6hTMy8tTLLNGwPwxtljLsvwGB8YSaWm/EOZZ61Kl1LaaBymlzsP8B8N0AOcB/BXACKXUOYeu2EJEwgDcBOAVpVSu1VcyzG8hPGCj7/0AbgPwnIg8Xd/zLW0cgPm96hkODLMA5qX4Hy3vsf8A4KDl2qtda+Nz5AMcfBmIqJmI7XuAiJqOiNwI8+y9U3OPhYhIbzgjJyIi0jEWciIiIh3j0joREZGOcUZORESkYyzkREREOubV3ANwRFBQkIqIiGjuYRAR6UpycvI5pVSw/SNrPb+9l5fXmzA/OIgTv+ZhAnCwsrLyoX79+p2xdYAuCnlERAR2797d3MMgItIVEan5COJ68fLyejM0NDQmODj4goeHB2+oagYmk0nOnj0bm5ub+yaAO20dw7+wiIioNr2Cg4MvsYg3Hw8PDxUcHJyPOh6nzEJORES18WARb36W30Gt9ZqFnIiISMdYyImI6Io0ZsyYLsnJyb617V+8eHFgdnZ2i4a0PXXq1I4bNmwwNHx0jtPFzW5ERETOtnbt2jpvBnzvvfeCrr766pKIiIiK+rb9yiuvnLJ/lHOwkBMRkV0z1u0Pz8ot8HNmm8ZQQ/GCu/sct3fczTff3D0nJ8e7rKzMY8qUKaeffPLJc35+fnFjx449+/333xvatGlTtX79+p87duxYWfPcPXv2+D7wwANdDxw4kA4AmZmZ3iNHjuyRlZWVFh8fH/XSSy8dHzhwYPGYMWMiUlJSWomIGjt27LnOnTtXHDx40O/+++/v5uvra9q9e3f63LlzQ7/44ou2ZWVlHv379y9cs2bNUQ8P2wvbCQkJESNGjMifOHHihUa/UHZwaZ2IiFzamjVrslNTU9P37duXtnz58pDc3FzPkpISj759+xanpaWlDxo0qGDWrFkdbZ3bt2/f0oqKCklLS/MGgFWrVgWMGjXqV8V1165dfjk5OS0OHTqUmpWVlfbYY4+dnzhx4oVevXoVr1q16ueMjIw0f39/NWPGjDMHDx5MP3ToUGpJSYnHhx9+2KYprt8et56Rf/vOCpw5+nNzDwPtu3TDkAmJzT0MIqIGc2TmrJWkpKSQLVu2tAWA3NzcFqmpqb4eHh546KGH8gBg0qRJ5++6664etZ0/atSovPfeey/ghRdeyP3000/brV279leFITo6uuz48eM+DzzwQPjIkSPzR48efclWO1u3bjW8/PLLoaWlpR4XL170io2NLQGQ78RLbRDOyImIyGVt3rzZsH37dsPu3bszMjMz02JiYkpKSkp+U7tEpNY2xo8ff2HDhg3tUlJSfEQEvXv3LrPeHxwcXHXw4MG0IUOGFCxdurT9vffeG1GzjeLiYpk+fXqXTz755EhWVlbauHHjzpWWlrpEDXXrGTlnwURE+nbx4kXPNm3aVBkMBtPevXt99+/f3woATCYT3n777XaJiYkX3nnnncD4+PiC2tro2bNnmYeHB+bMmdNx9OjReTX35+TkePn4+JgmTJhw0Wg0lk2aNKkrAPj7+1fl5+d7AkBxcbEHAISGhlbm5+d7fPbZZ+1Gjhyp+fvfjnDrQk5ERPqWkJCQv2LFimCj0RjbvXv30j59+hQBQMuWLU2pqakte/bsGWowGKo++eSTOt9Hveuuu/LmzZvXKSkp6WTNfdnZ2S0efPDBCJPJJADw97///QQA3H///ef+9Kc/dZkxY4Zp9+7d6WPHjj0bGxvbs1OnTuXV46iLiDTJw3R0kUfev39/xWetExHVj4gkK6X6N/T8/fv3Z/fp0+ecM8fkLH5+fnHFxcV7m3sctRk6dGiPadOmnR45cmStKwX1sX///qA+ffpE2NrnEuv7RERE7uKee+6JKCkp8bj11lsLm6I/Lq03g4ufHUH5KburMk7l3bEV2o7s3qR9EhFpxdZsfPz48Z1/+uknf+ttjzzyyOknnnjivFbjaI4+a2IhJyIit7B69epjV0KfNbGQNwPOjImIyFn4HjkREZGOsZATERHpGAs5ERGRjrGQExHRFUnLPPL4+Pio7777zqlpcbXhzW5ERHRF0jKPvCmxkBMRkX0bHgvHmTTnzjDbxxZj1L/cMo8cAD744IN2jz32WJeCggLP119/PXvYsGGaPCCGS+tEROTS9JpHXllZKQcOHEhPSko6/ve//93m+JyBM3IiIrLPgZmzVvSaR37PPfdcAICBAwcWzZgxw7sBl+4QzWbkIhIuIt+KSLqIpIrIEzX2PykiSkSCtBoDERHpm57zyH19fRUAeHl5oaqqqvYBNpKWM/JKANOVUntExAAgWUS+VkqliUg4gFsAaPpou2/fWYEzR+tMtmsS7bt0YzY6EVEDMI/cPs0KuVIqB0CO5fsCEUkHEAYgDcAiAH8FsFGr/omISP/0nEfeVJokj1xEIgB8B6AXgBsB3KSUekJEsgH0V0rVmXfLPHIiovpjHrn7qCuPXPOb3UTEH8B6AFNhXm5/GsCtDpyXCCARADp37qzlEImIiHRL00IuIi1gLuJrlFKfiEhvAF0B7LfcmNAJwB4RiVdK5Vqfq5RaAWAFYJ6RazlOIiLSF+aR/49mhVzMlXolgHSl1MsAoJQ6AKC91THZcGBpnYiIyJ4rNY9cywfCDAIwHsBQEdln+bpdw/6IiIiuOFretf4fAHV+bk4pFaFV/67s4mdHUH6qaW949O7YCm1Hdm/SPomISHt8spuby8vLQ3l5OQouliHr3HdObz80NBTDhw93ertEROQYFvJm0JQz4+StW5Gbm2v/QCIi0iUWcjfH2TIRkW1jxozp8te//vV0v379Sm3tX7x4ceCdd955qSExptXJatdff31x40daNxZyIiK6IjGPnIiIrhjPfv9s+OELh52aR96jXY/ieYPmuW0eOQBUVVXhnnvuiejUqVP5yy+/fKpmP88999yZBr+AFswjJyIil6bXPPKKigoZNWpU18jIyNLFixefstWPM14fzsiJiMguR2bOWtFrHvmjjz7aZdSoUXlJSUm59emnvjgjJyIil6XnPPL+/fsX7tixo3VxcbE42k9DsJATEZHLspdHDgDOyCOvqqrChAkTLv7jH/84eeDAAT/Afh65vbFPnjz53K233po/YsSI7hUVFbX201hcWiciIpel9zzyuXPnnp42bZrnXXfd1XX27Nm5tvpprCbJI28s5pETEdUf88jdR1155FxaJyIi0jEurRMRke4wj/x/WMiJiMgtMI+ciIiIdIeFnIiISMdYyImIiHSMhZyIiEjHWMiJiMitTZ06teOGDRsMWrUfFxcXXf394MGDIw0Gw9VDhgyp9dnvzsa71omIyK298sorp2xtr6yshJdX48vg3r17M6q/f/LJJ3OLioo83njjjeBGN+wgFnIiIrLr1Oynw8sOHXJqHrlPZGRxxxee1zSPHAASEhIiRowYkT9x4sQLYWFhvf/4xz+e+/bbb1tPnjz5TEFBgefbb78dXFFRIREREWXr1q37xWAwmI4fP+41adKkLseOHfMBgCVLlhy95ZZbbD6W1fopc7///e8LNm/erNns3xYurRMRkUtrTB65Lb6+vqbk5OTMxMTEC2PHjr1w8ODB9MzMzLSoqKiSxYsXBwHAlClTOg8ePLggMzMzLTU1Na1v376l2l1h43BG7ua2bt2K3NxczdoPDQ3F8OHDNWufiFyDIzNnrTQ2j7ym+++//0L198nJyS3nzJkTVlBQ4FlUVOR5ww035APAzp07DevWrfsFALy8vBAYGFjl3KtyHs7IiYjIZTkjj7wmg8Fgqv4+MTGx65IlS45lZWWlzZw581RZWZnu6iJn5G5Oq9lyVtY8FBSmAwCS96x2SpsG/xgYjc86pS0icg/28sgTExMv2Msjr0txcbFH586dK8rKyuTDDz8M6NChQwUADBo0qGDBggXBc+bMOVNZWYlLly55BAQEmOy11xx095cHERFdORISEvIrKyvFaDTGzp49u6ONPPKY7777zvDiiy/mNKT9WbNmnYqPj48ZPHiwMTIy8vL74MuWLTu2fft2g9FojO3Vq1fsnj17WjrSXr9+/aLGjx/fbdeuXa1DQkKuWr9+feuGjKs+mEdOROSmmEfuPpolj1xEwkXkWxFJF5FUEXnCsn2eiKSIyD4R+UpEHL7TkIiIiH5Ny/fIKwFMV0rtEREDgGQR+RrAAqXUswAgIn8GMAfAFA3HQUREbqap88hzc3M9b7zxxqia27dt25YZGhrarHe0a1bIlVI5AHIs3xeISDqAMKVUmtVhrQC4/to+ERG5PC2zwUNDQ6syMjLS7B/Z9JrkrnURiQAQB+BHy8/PA7gfQD6AIU0xBiIiInek+V3rIuIPYD2AqUqpSwCglHpaKRUOYA2Ax2s5L1FEdovI7rNnz2o9TCIiIl3StJCLSAuYi/gapdQnNg55H0CCrXOVUiuUUv2VUv2Dg5vs2fNERES6ouVd6wJgJYB0pdTLVtsjrQ67E0BGzXOJiIjIMVrOyAcBGA9gqOWjZvtE5HYA80XkoIikALgVwBMajoGIiK5wTZVHvnPnzpZXX311dI8ePXoajcbYN954o51WfVrT8q71/wCw9fDbz7Xqk4iIqKamyiP39/c3rV69+pfevXuXZWdntxgwYEDM6NGjLwUFBWn68TQ+a52IiOz6v1Xp4XknC52aRx4Q5l980/0xbpNHftVVV5VVb4uIiKgICAiozMnJ8dK6kPNZ60RE5NL0mEf+7bff+lVUVEhsbGyZ/aMbhzNyIiKyy5GZs1b0lkd+9OjRFhMnTuy2cuXKXzw9PR09rcE4IyciIpeltzzyvLw8j+HDh/eYM2fOyZtuusnmUryzsZATEZHLspdHDgDOzCOv3l6dRw6Yb4rLy8uzWy9LS0vljjvu6HHvvfeenzRp0gV7xzsLCzkREbksPeWRv/XWW+1++ukn//fffz8oOjo6Njo6Onbnzp0O5Zg3BvPIiYjcFPPI3Uez5JETERGR9njXOhER6Q7zyP+HhZyIiNzClZpHzqV1IiIiHWMhJyIi0jEWciIiIh3je+TUIFlZ81BQmO7UNg3+MTAan3Vqm0RE7o4zciIiumJkZmZ6R0ZG9mzucTgTZ+TUIJw5ExG5BhZyIiKy68tlr4SfO37UqXnkQeFdim97ZKrmeeQ7duzwe+ihhyJatmxpuuaaawqrt1dWVuKxxx7r9P333xvKy8vl4YcfPjNjxoxzAPDMM8+EfPTRR4Eigptuuil/6dKlJxcuXBhkK7s8ISEhwtfX13T48GHfkydP+ixfvvyXd955Jyg5OblVXFxc0fr167Nru7ZFixYFvfrqq6Ht27ev6NatW6m3t7datWpVvT5Gx6V1IiJyaY3NI3/wwQcjXn755WP79u3LsN7+yiuvBLVp06bq4MGD6fv3709/9913gzMyMrw/+uij1lu2bGmXnJyckZmZmfbcc8/lAkBt2eUAkJ+f77Vr166s+fPnHx8zZkzkjBkzTh86dCg1IyOjZW3PW8/Ozm7x0ksvdfjxxx/Td+zYkXXo0CHfhrw+nJETEZFdjsyctdKYPPLz5897FhQUeN5xxx2F1cd+8803bQDg3//+d+uMjAy/TZs2tQOAgoICz7S0NN+vv/669bhx485Vx52GhIRUAbVnlwPAHXfccdHDwwN9+/YtDgwMrIiPjy8BAKPRWHLkyBGfgQMHltQc244dO1pdc801BdXtjx49+kJWVla9izkLORERuSzrPHKDwWCKj4+Pqk8euVKqrn2ycOHCYwkJCZest3/++eetbZ2TmJjYdd26dYevvfbaksWLFwdu377dUL3P19dXAYCnpye8vb0vp5F5eHigsrLS5gCcFVrGpXUiInJZjc0jDwoKqvL396/68ssv/S3HXs4cv+WWW/KXLVsWXFZWJgCQkpLic+nSJY9hw4ZdWr16dVBBQYEHAJw+fdoTqD27vKEGDx5c9OOPPxrOnj3rWVFRgY0bN7ZrSDuckRMRkctKSEjIX7FiRbDRaIzt3r17qY088lCDwVD1ySef/FxbGytXrsyuvtlt6NChl2ff06ZNO5edne3Tu3fvGKWUBAQEVHz++edH7r777kt79uzxu/rqq2NatGihbr755vwlS5acrM4uDwsLK4+JiSkuLCz0bMy1de3atWLatGk5AwYMiGnfvn2F0WgsadOmTb0DWJhHTkTkpphH7vry8/M92rRpY6qoqMBtt93WY8KECefuv//+izWPYx45ERGRC5oxY0bH6OjoWKPR2LNz585l48aN+00Rt4dL60REpDtNnUfeWFdddVV0eXn5rybPq1at+mXFihUnGtu2Q4VcRHwAJACIsD5HKfX3xg6AiIjIGbTMI2+slJSUDPtHNYyjM/KNAPIBJAMo02owREREVD+OFvJOSqlh9WlYRMIBrAIQCsAEYIVS6lURWQBgJIByAEcATFRK1fs9ASIiInK8kO8Ukd5KqQP1aLsSwHSl1B4RMQBIFpGvAXwN4CmlVKWIJAF4CsDM+g2biKjp7fgoC+eOF9o/0ImCwv0x+A/GJu2T9KXOQi4iBwAoy3ETReRnmJfWBYBSSl1V27lKqRwAOZbvC0QkHUCYUuorq8N+AHB34y6BiIjoymVvRj7CGZ2ISASAOAA/1tg1CcBaZ/RBRKQ1zoz1LzMz03vEiBGRhw4dSm3usThLnYVcKXUUAERktVJqvPU+EVkNYLzNE399nD+A9QCmKqUuWW1/Gubl9zW1nJcIIBEAOnfubK8bIiIiXausrISXV/0/Fe7oGT2tfxARTwD97J0kIi1gLuJrlFKfWG1/AObZ/k2qlkfLKaVWAFgBmJ/s5uA4iYhIA3nrssIrcoucmkfeIrRVccDdxis6j9zPzy8uMTHx9DfffNN6wYIFJ2677bZ634RR55PdROQpESkAcJWIXBKRAsvPZ2D+SFpd5wqAlQDSlVIvW20fBvPNbXcqpYrrO2AiIrqyuGseOQCUlJR49OrVqyQlJSWjIUUcsL+0/iKAF0XkRaXUU/VsexDMS+8HRGSfZdtsAIsB+AD42hIT94NSako92yYioibkyMxZK+6aRw6YY08nTJhwoTGvj6NL67NF5C4A18F8F/sOpdSGuk5QSv0H5rvba/q8fkMkIqIrlTvnkQOAt7e3qSHvi1tz9Ox/AegB4APLz1NE5Bal1GON6p2IqBFyX3gBZemaPfnyMp+YaITOnq15P/Rb9vLIExMTLziaR37bbbcV2sojHzFiRIGPj49KSUnxiYiIqBg2bNil559/vuPDDz+cZzAYTKdPn/YMCQmpqplH3qFDh4qmeh3q4mghvwFAr+ob00TkXQD1eTgMERFRvblzHrmzOJRHLiKfAJhm9XG0LgDmK6X+qPH4ADCP3J08e+gEDhbafKvosl7+LTEvslMTjYjIfTGP3H3UlUfu6Iw8EEC6iPzX8vMAALtEZBMAKKXubPQoiYiIqN4cLeRzNB0FXTE40yYiZ3CXPPLqu9sbw6FCrpTabllOj1RK/VtEWgLwUkrZvLmAiIioqV2peeR1PhCmmog8DGAdgOWWTZ0A1PnxMyIiItKeQ4UcwGMwP+DlEgAopQ4BaK/VoIiIiMgxjhbyMqVUefUPIuIF84NhiIiIqBk5erPbdhGZDaCliNwC4FEAn2k3LCIi17PjoyycO96gx2E3WFC4P+NTncjVYkxvuOGGHuvXr/8lKCioqqFtOFrIZwF4EOaHwEyG+TGrbza0UyIiIgK2b99+uLFtOHrXuklENgDYoJQ629hOiYj0iDPj5uHOMaZhYWG9d+/end6hQwebY3dEnYXcEkX6HIDHYQ5AERGpAvCaUurvDe2UiIj0ZcOGDeFnzpxxah55+/bti0eNGmU3VW3NmjXZISEhVYWFhRIXFxc7bty4C9Uxpm+88caJJ598ssOsWbM6rlq1yubHzx588MGIRYsWHbvjjjsKJ0+efPlhFtYxpiUlJTJgwIDokSNHXkpJSfGtjjGtftY6YI4xnT59+jkA+POf/9xx8eLFQU8//fQZ4H8xpu+//37bMWPGRH7zzTcZ/fr1K7nqqqtidu7c2bK29DNnsHez21SY71YfoJQKVEoFALgGwCARmabVoIiIiKolJSWFREVFxfbr1y+mthjT//73v/62zrUVY1q979///nfrjz76KDA6Ojo2Li4u5sKFC172Ykz79esXZTQaY9evXx+YmprqW92WrRhTT0/PyzGmWr4+9pbW7wdwi1Lq8rN2lVI/i8g4AF8BWKTl4IiIyDU4MnPWgrvHmDqDvRl5C+siXs3yPnkLbYZERERkZi/GFAAcjTG1HPubGNOysjIBgJSUFJ9Lly55DBs27NLq1auDCgoKPACgemm9ZoyptlfuOHsz8vIG7iMiImo0xpjaV2eMqeXGtiJbuwD4KqWaZFbOGFMiovpjjKn7aHCMqVLKJf7aICIiItscfSAMEZHLyX3hBZSlaxYqdZlPTDRCZ8/WvB9yHGNM/4eFnIiI3MKVGmPKQk5EusVZMpHj6WdERETkgljIiYiIdIyFnIiISMdYyImI6IqRmZnpHRkZ2bO5x+FMLOREREQ6ptld6yISDmAVgFAAJgArlFKvisg9AOYCiAEQr5TiI9uIiFxcWvrM8KLCLKfGmLbyNxbHxiTZDWNx1zzyNWvWtJk3b14YAJSWlnpUVFTIyZMnD9T3ddTy42eVAKYrpfaIiAFAsoh8DeAggLsALNewb3JRzx46gYOFdT//oJd/S8yL7FTnMUR05XDXPPKxY8fmjx07Nh8Abr/99m6DBw+2Gfxij2aFXCmVAyDH8n2BiKQDCFNKfQ3UHjlHRESux5GZs1aSkpJCtmzZ0hYAassjv+uuu3rYOtdWHvk333zTBjDnkWdkZPht2rSpHQAUFBR42ssjnzNnTlhBQYFnUVGR5w033JBf3Y+tPHIAl/PIbRXyas8880yIr6+v6amnnjrbkNenSR4IIyIRAOIA/NgU/ZHr4kybiOrD3fPIN27caNiwYUPADz/80OAnv2l+s5uI+ANYD2CqUuqSveOtzksUkd0isvvs2Qb9kUJERDrnznnkWVlZ3k888USXdevWHfH39689itQOTWfkItIC5iK+RinqLGfMAAAczklEQVT1SX3OVUqtALACMMeYajA8IiJyce6cR758+fLA/Px8z1GjRvUAgJCQkPLt27cfrm87deaRN4aY1yXeBZCnlJpqY/82AE86ctc688iJiOqPeeTuo8F55I00CMB4AAdEZJ9l22wAPgBeAxAMYIuI7FNK3abhOIiIiNyWlnet/wdAbW/wf6pVv0RE5P6YR/4/jDElIiK3cKXmkfMRrURERDrGQk5ERKRjLOREREQ6xvfIye387bNUpJ1y+NlDtYrt2BrPjXSrtEMickOckRMR0RXDHfPIOSMnt8NZNBFdSVjIiYjIrqnpx8Izikqdmkce3cq3+JWYzldsHvmiRYuCDh482HLlypXHAWDhwoVB6enpvm+++eaJ+ryOXFonIiKXtmbNmuzU1NT0ffv2pS1fvjwkNzfXszqPPC0tLX3QoEEFs2bN6ljb+Q8++GDEyy+/fGzfvn2/+iy3dR75/v370999993gjIwM748++qh1dR55ZmZm2nPPPZcLmPPIDx48mJ6ZmZkWFRVVsnjx4qDqtqrzyOfPn398zJgxkTNmzDh96NCh1IyMjJY7d+5sWcu48r766qs21aEt7733XlBiYmK9H17DGTkREdnlyMxZK+6aR966dWvToEGDCtauXdumd+/epRUVFdKQJ72xkBMRkcty9zzyxMTEc88//3yo0WgsHTduXIMCari0TkRELsud88gBYOjQoUU5OTnen376aeCDDz6Y15A2OCMnIiKX5c555NVGjRp1ISUlxS84OLiqIedrlkfuTMwjJyKqP+aR68OQIUN6TJ069fTvf/97m6sKQN155FxaJyIiagbnzp3zjIiI6OXr62uqq4jbw6V1IiLSHXfJI8/Ozj7Y2LZZyImIyC0wj5yIiIh0h4WciIhIx1jIiYiIdIyFnIiISMdYyImISFeys7NbDBs2rFtt+8+dO+c5f/784Ia2HxcXF93Qc5sD71oncqats4DcA85vN7Q3MHy+89sl0qGIiIiKL774otYnuZ0/f95z5cqV7WfNmnW2Ie3v3btXszvMtcBCTkREds1Ytz88K7fAqXnkxlBD8YK7+9SZqvbII4+EdenSpby6KP/lL3/paDAYqt5///2gQ4cOpe7evdt34sSJXSsqKsRkMmH9+vVHnnrqqbDjx4/7REdHx95www2X/vnPf54aNmxYj/z8fM/KykqZM2fOqXHjxl2srU+9PTWOhZzImThrJnKqcePG5U2dOrVzdSHfuHFjuyVLlhx9//33gwDgtddeC3700UdPP/LII3mlpaVSWVmJhQsXnhgxYkTLjIyMNACoqKjAli1bDgcEBJhycnK8rrnmmuj77rvvooeHe7y7zEJORER22Zs5a2XQoEEl58+f98rOzm6Rk5Pj1aZNm6pu3bqVV++/9tpri1566aUOJ06c8L733nsv9O7du6xmGyaTSaZOndrphx9+8Pfw8MCZM2e8T5w44dW5c+fKpr0abbCQE7mppP8mISOv8W/1RQdEY2b8TCeMiKhhRo4ceeG9995rl5ub2yIhIeFXUZ9TpkzJGzx4cNGnn37aZvjw4calS5dmR0VF/aqYL1++POD8+fNeBw4cSPfx8VFhYWG9bWWa65VmhVxEwgGsAhAKwARghVLqVREJALAWQASAbAB/UEpd0GocRESkb+PHj897+OGHIy5cuOC1ffv2zNLSUqnel5aW5h0TE1PWs2fPMz///LPPvn37WsbHxxcXFRVdLtT5+fmeQUFBFT4+Puqzzz4znDp1yrt5rkQbWs7IKwFMV0rtEREDgGQR+RrABAD/p5SaLyKzAMwCwD/3iZyMs2hyF/379y8tKiryCAkJKe/SpUtFZmbm5UK8evXqgI8//jjQy8tLBQcHV7z44ounQkJCqvr161cYGRnZc+jQoflz587NHT58eI9evXrF9OzZs7hr166ldfUnInXtdjlNlkcuIhsBLLF83aiUyhGRDgC2KaWi6jqXeeRERPXnznnkWsnNzfXs27dv7KlTpzT4HGnDNXseuYhEAIgD8COAEKVUDgBY/rd9U4yBiIioLtnZ2S1+97vfxTz22GOnm3ss9aH5zW4i4g9gPYCpSqlLji5ZiEgigEQA6Ny5s3YDJCKiK05ubq7njTfe+JvV4B9++CE9NDS0qjnG1FCaFnIRaQFzEV+jlPrEsvm0iHSwWlo/Y+tcpdQKACsA89K6luMkIqIrS2hoaFX158z1TrOldTFPvVcCSFdKvWy1axOAByzfPwBgo1ZjICIicndazsgHARgP4ICI7LNsmw1gPoCPRORBAMcA3KPhGIiIiNyaZoVcKfUfALW9IX6TVv0SERFdSdzmyTZERERXIhZyIiLSFS3zyDMzM70jIyN7Nnx0TY/PWie387fPUpF26lKj24nt2BrPjdTVf89EVwSt88j1hoWciIjs2/BYOM6kOTWPHO1jizHqXy6XR15VVYV77723y+7du/1DQkLKv/zyy8P+/v4u+zFoFnJyO5xFE7mP5sgjP3bsmO97773388CBA4/efvvt3VatWtXu0UcfzbN5sAtgISciIvvszJy10hx55GFhYWUDBw4sAYC4uLji7OxsH+2usPF4sxsREbm06jzyNWvWBNjKI9+4cePhli1bmoYPH27ctGmToeb51nnkGRkZaYGBgRV15ZF7e3tfXkb39PRUlZWVLh2Hxhk5ERG5NOaR140zciIicmk188it961evTrAaDT2jI6Ojj106JDv5MmTz4eGhl7OI588eXKnhx56KG///v2tevXqFfPee+8F2Msj15smyyNvDOaRExHVH/PI3Uez55ETERGRNvgeORERXXFqyyPftm1bJvPIiYiIXBzzyImIiMglsJATERHpGAs5ERGRjrGQExER6RgLORER6QrzyH+NhZyIiHTF0TzyphxTc+LHz4iIyK5nv382/PCFw07NI+/RrkfxvEHzXC6PvFpaWpp3QkJCj9dffz27VatWppr92Epaaw6ckRMRkcsaN25c3vr16wOqf964cWO73/3ud0XVP1fnkWdkZKSlpKSkd+3atXzhwoUnwsPDyzIyMtKWL19+ws/Pz7Rly5bDaWlp6du3b8+aPXt2J5PJVGe/+/fv90lISOixcuXKX2644YZiW/1oeNn1whk5ERHZZW/mrJXmyCPPy8vzGjVqVI+PP/74SP/+/Usd7ae5sJATOdPWWUDuAee3G9obGD7f+e0S6UB1Hnlubm4LW3nkgwcPLvr000/bDB8+3Lh06dLsqKioXxVZ6zxyHx8fFRYW1ruuPHKDwVDVoUOH8m3btvlXF3Jb/dx5550F2lxx/bCQExGRS2vqPPIWLVqoL7744siQIUMi/f39TVOmTMmz1Q8LOZE74qyZyOlq5pFnZmZeLsSrV68O+PjjjwO9vLxUcHBwxYsvvngqJCTkch750KFD8+fOnZs7fPjwHr169Yrp2bNnsSN55K1btzZ9+eWXh2+88Uajv7+/KTU11bdmP9peteOYR05E5KaYR+4+mEdORETkpjRbWheRtwCMAHBGKdXLsq0PgNcB+APIBjBWKXVJqzEQXcmS/puEjLyMRrcTHRCNmfEznTAiItfBPHLHvANgCYBVVtveBPCkUmq7iEwCMAPAsxqOgYiI6DfcKY9cs0KulPpORCJqbI4C8J3l+68BfAkWciJNcBZNdGVo6vfIDwK40/L9PQDCm7h/IiIit9LUhXwSgMdEJBmAAUCtj7gTkUQR2S0iu8+ePdtkAyQiItKTJi3kSqkMpdStSql+AD4AcKSOY1copforpfoHBzcojY6IiMjtNWkhF5H2lv/1APAMzHewExER2RQXFxetVdtr1qxpM3v27FAA2Lp1q39sbGyMl5dXv7fffrudVn1qQcuPn30A4EYAQSJyAsBzAPxF5DHLIZ8AeFur/omISP/27t37m89QVlZWwsur8eVr7Nix+QDyAaBbt27lb7/9dvb8+fNDGt1wE9PyrvU/1rLrVa36JCIibZya/XR42aFDTs0j94mMLO74wvN1pqr5+fnFFRcX7928ebNh3rx5Hdq3b1+Rlpbmd+TIkdSbb765e05OjndZWZnHlClTTj/55JPnAGDdunWt58yZE1ZVVSUBAQGVu3btyrLV9uLFiwN3797datWqVceioqLKAcDDQ3/PSeOz1omISBdSUlJa7d27NzU6OrocANasWZMdEhJSVVhYKHFxcbHjxo27YDKZ5PHHH4/Ytm1bRnR0dPnp06c9m3vcWmMhJyIiu+zNnJvCVVddVVRdxAEgKSkpZMuWLW0BIDc3t0Vqaqrv6dOnveLj4wuqjwsJCdHVU9oaQheFPDk5+ZyIHHXg0DawvN/RBLTqy5ntOqOthrbRkPOCADCgwbma8r8JLbnSdejp35kuzhqIK/Dz8zNVf79582bD9u3bDbt3784wGAym+Pj4qJKSEg+lFESkrmbcji4KuVLKoc+ficgKpVSi1uPRsi9ntuuMthraRkPOE5HdjUlqot9qyv8mtORK1+EO/864g4sXL3q2adOmymAwmPbu3eu7f//+VgAwZMiQounTp3fJyMjwrl5ad/dZuf7e1a/bZ27QlzPbdUZbDW2jKX8XVDt3+T240nW4w78zupeQkJBfWVkpRqMxdvbs2R379OlTBAAdO3asXLx4cfbo0aN7REVFxY4ePbqbI+1t377dLyQk5KrPP/+83bRp07r06NGjp7ZX4Dy6yCOnKwNn5ESuhXnkroN55KQXK5p7AEREeqOL98jpyqCUYiEnIqd79dVXA5ctW/arB70MGDCgcPXq1ceaa0zOxKV1IiKyiUvrroNL60RERG6KhZyIiEjHWMjJpYlIKxFJFpERzT0WIiJXxEJOTUpE3hKRMyJysMb2YSKSKSKHRWSW1a6ZAD5q2lESEekHCzk1tXcADLPeICKeAP4FYDiAWAB/FJFYEbkZQBqA0009SCJyDU2VRz537tyQ7t279zQajbHXXnutMSsry1urfp2NHz+jJqWU+k5EImpsjgdwWCn1MwCIyIcAfg/AH0ArmIt7iYh8rpQygYiuGE2VR96vX7/i6dOnpxsMBlNSUlLwtGnTOm3ZsuXnRnfSBFjIyRWEAbBOVjoB4Bql1OMAICITAJxjESdqPv+3Kj0872ShU/PIA8L8i2+6P8Yl8shHjhxZUL39uuuuK1y7dm2gM69VSyzk5ApsRRVdfsCBUuqdphsKEbmqpsojX758efDNN9/sKml7drGQkys4ASDc6udOAE4101iIyAZ7M+em0BR55EuXLg3Yv3+/3/LlyzOdO3rtsJCTK/gJQKSIdAVwEsC9AO5r3iERkavROo98w4YNhpdeeqnDjh07Mlu2bKmbx57yrnVqUiLyAYBdAKJE5ISIPKiUqgTwOIAvAaQD+Egpldqc4yQi11ZXHvmPP/5oyMjI8AYAR5fWv//++5Z/+tOfumzcuPFwWFhYpZZjdzbOyKlJKaX+WMv2zwF83sTDISKdSkhIyF+xYkWw0WiM7d69e6mtPHKTyYTAwMCKnTt3HrLX3owZM8KLi4s977nnnu6Wdsq/+eabw1pfhzMwNIWIiGxiaIrrYGgKERGRm+LSOhERuTXmkRMR0RWJS+uug0vrREREboqFnIiISMdYyImIiHSMhZyIiEjHWMiJdExEXhORPSIyoLnHQtQU/vKXv3ScM2dOiP0jrxws5EQ6JSKtALQHMBnAiGYeDhE1E36OnHRBRBYBOKqUesXy85cAjiulHrL8vBDASaXUy07ss1Ap5e/E9toCuE8ptdTycwSAzUqpXg6c2xLAFwCGKqWqAEApVSQiHQBsA9BZRLwB/NtyjK6eFU2u78tlr4SfO37UqXnkQeFdim97ZKrdVLWZM2eGrl27Nqhjx47lgYGBFXFxccWpqak+U6ZM6ZyXl+fl6+trevPNN4/GxcWVHj9+3GvSpEldjh075gMAS5YsOXrLLbcU1ZZd7ufnF/fAAw+c+e6771q3adOm6vnnnz8xc+bM8FOnTnknJSUdGzt2rM0404KCAo8xY8ZEHD582DcyMrL0+PHj3kuWLDl2/fXXFzvzNXIEZ+SkFzsBDAQAEfEAEASgp9X+gQC+b4Zx1UdbAI828NxJAD6pLuIAICKBAPwAFACoUkqVA/g/AGMaO1AiV7Fjxw6/Tz/9NODAgQNpmzdvPlwdjvLQQw91Wbp06bHU1NT0BQsWnHjkkUc6A8CUKVM6Dx48uCAzMzMtNTU1rW/fvqWAObs8NTU1fd++fWnLly8Pyc3N9QSAkpISjyFDhhSkpqamt2rVquqZZ54J27FjR9bHH398eN68eWG1jWvBggXBbdu2rcrKykqbO3fuqbS0tFZN8XrYwhk56cX3ABZZvu8J4CCADiLSDkAxgBgAe0VkA8zZ5r4AXlVKrQAAEUmCeUZfPRueC3MBPA3gzwC8AfwI4FHrYmk5dpytYywz6q0A/gPzHxInAfxeKVUiIs8CGAvgOIBzAJIB9AfQXUT2AfgawL8AeIrIGzXPt3H9Y/HbaNdnALwE4GEAsTD/sbMBwIsA1th9RYnqwZGZsxa+/fZb/9tvv/2iwWAwAcCtt956sbS01GPv3r3+1QEnAFBeXi4AsHPnTsO6det+AQAvLy8EBgZWAbazy0NDQ4tatGih7r777ksA0LNnzxIfHx+Tj4+Pio+PLzl58qR3bePauXOn/xNPPHEGAAYMGFBqNBqbfCZejTNy0gWl1CkAlSLSGeaitwvmonotzAUyxTIjnaSU6mfZ9mfLrBUAPsSvZ6p/ALDbsm2QUupqAFUwF8zLRCTGzjGRAP6llOoJ4CKABBHpDyABQByAuyxjAYBZAI4opa5WSs2o7fya125ZMu+mlMq22hZheR3Wwhz9Wr06cRAAb3wjt1IzX9xkMsFgMFRmZGSkVX/9/PPPtUYfW2eXZ2ZmpsXExJSUlJR4AICXl5fy8DCXQg8PD/j4+CgA8PT0RFVVVa3B5q70VFQWctKT72EuXtWFfJfVzzstx/xZRPYD+AHmmXkkACil9gJoLyIdRaQPgAsAegPoB+Anyyz5JgDdavR5k51jflFK7bN8nwwgAsB1ADYqpUqUUgUAPqvjmmydX1MQzEXe2j8A/F2Z/zW5XMgtqwnlImKoo08i3Rg6dGjhli1b2hYWFsqFCxc8vv7667Z+fn6mTp06lb/11lvtAHNh37VrV0sAGDRoUMGCBQuCAaCyshJ5eXketWWXN8bAgQMLP/zww3YAkJyc7JuVldWysW02FJfWSU+q3yfvDfPM8ziA6QAuAXhLRG4EcDOAa5VSxSKyDeYl9mrrANwNIBTmGboAeFcp9VQdfdo7pszq+yoALS3nOMrW+TWVwOo6RORqmGf614nIvyz7Dlgd7wOgtB5jIHJZ1113XfHo0aPzevXq1TMsLKwsPj6+EAA++OCDnx9++OEuSUlJHSorK2X06NF51157bcmyZcuOTZgwoYvRaAzy8PDAkiVLjtaWXd4YM2bMOPuHP/whwmg0xvbq1as4KiqqpF27dlX2z3Q+hqaQblgK2CcAflZK3WzZlgwgDEAvAIMAPKSUGiki0QD2ARimlNpmObYngDdgnuHeAKAdgI0wL5ufEZEAAAal1FHL8YUA4ms7puZd5yLyJAB/AFsALIf5jw4vmGfabwB4G8AepVQXy/E2z1dKzbVx7ccBRCqlSi137P9TKfV/ln0hAPYqpTpa3kr4j1IqplEvNhEYmlKXyspKlJeXi5+fn0pNTfW59dZbjUeOHDno6+urSVGtKzSFM3LSkwMwF+H3a2zzV0qdE5EvAEwRkRQAmTAvr1+mlEq1LDmfVErlAMgRkWcAfGW5E74CwGMAjlqdk2bvmJqUUj+JyCYA+y3H7QaQr5Q6LyLfi8hBmG+S+1c9rv0rmGfgJgCtqou4pb/TItLK8kfGEACf16NdImqAgoICj8GDB0dVVFSIUgqLFi06qlURt4czciINiIi/UqpQRPwAfAcgUSm1pxHtxQH4i1JqvJ3jPgHwlFIqs6F9EVXjjBxYv35966effrqT9bbw8PCyr7/++khTjoMzcqKmt0JEYmF+//rdxhRxwHyznoh8KyKeNT8eV81yd/sGFnEi50lISLiUkJCQ1tzjqAsLOZEGlFI1P/PtjDbfsrO/HMAqZ/dLRK6NHz8jIiLSMRZyIiIiHWMhJyIi0jEWciIi0g3mkf8WCzkREVEjVFRUNGv/vGudiIjsyluXFV6RW+TUPPIWoa2KA+426jKPfPHixYFbt25tU1ZW5lFcXOzxww8/ZDnztakPzsiJiMhluWoeOQDs2bPH/4MPPvilOYs4wBk5ERE5wJGZsxZcNY8cAAYPHnwpJCSkWYJSrLGQExGRS6srj9yR863zyA0Ggyk+Pj6qsXnkAODn52dqyPU4G5fWiYjIZblqHrkrYSEnIiKXZZ1HPmLEiO7WeeRvv/12UFRUVGxkZGTP9evXtwWAZcuWHdu+fbvBkhMeu2fPnpYJCQn5lZWVYjQaY2fPnt3RGXnkroTpZ0REZBPTz1xHXelnnJETERHpGG92IyIiqoWr5JHXhYWciIioFnrII+fSOhERkY6xkBMREekYCzkREZGOsZATERHpGAs5ERHphqvkkX/33Xd+EyZMCG/ucQC8a52IiKjerr/++uLrr7++uLnHAbCQExGRAzZs2BB+5swZp+aRt2/fvnjUqFG6zCPfvHmzYeHChSHffvvtYWe+Jg3BpXUiInJZrpxH7io4IyciIrscmTlrwZXzyF0FCzkREbk0V80jdxVcWiciIpfFPHL7WMiJiMhlMY/cPuaRExGRTcwjdx3MIyciInJTvNmNiIioFswjJyIi0jHmkRMREZGmWMiJiIh0jIWciIhIx1jIiYhIN1wlxtSVsJATERHpGO9aJyIiu9LSZ4YXFWY5Nca0lb+xODYmSZcxpmPGjOlS/ajX06dPt5g0adKZhQsX5jjz9XEUZ+REROSyXDXGdO3atUczMjLSNm3adLht27aVkydPPt8Ur4ctnJETEZFdjsycteDKMabFxcWSkJDQfdGiRceMRmO5Nq+AfSzkRETk0lw1xnT8+PFdRo4ceWHUqFEFDbkuZ+HSOhERuSxXjTF98cUXgwsLCz1feOGF3Ma21Vgs5ERE5LJcNcZ0yZIloZmZmS2jo6Njo6OjY//5z38GN7bNhmKMKRER2cQYU9fBGFMiIiI3xZvdiIiIasEYUyIiIh1jjCkRERFpioWciIhIx1jIiYiIdIyFnIiISMdYyImISDeYR/5bLOREREQ6xo+fERGRXVPTj4VnFJU6NY88upVv8SsxnXWZR96vX7+o11577djAgQNLAKBv377Ry5YtO3rNNdeUOPM1cgRn5ERE5LJcNY98woQJ5958880gAEhJSfEpLy+X5ijiAGfkRETkAEdmzlpw1TzyCRMmXFiwYEGHsrKyE6+//nrQfffd12zPpGchJyIil+aKeeQGg8E0ePDgS++//37bTZs2BSQnJzfb09+4tE5ERC7LVfPIAWDKlCnnZs6cGd6nT5+ikJCQKme02RAs5ERE5LJcNY8cAAYPHlzcqlWrqokTJzZr1CvzyImIyCbmkdctOzu7xY033hh15MiRg56enpr2xTxyIiIiJ1qyZEng7373u5g5c+ac1LqI28Ob3YiIiGpRVx75448/fr65xmWNhZyIiKgWzCMnIiI9M5lMplo/gkVNw/I7MNW2n4WciIhqc/Ds2bNtWMybj8lkkrNnz7YBcLC2Y7i0TkRENlVWVj6Um5v7Zm5ubi9w4tdcTAAOVlZWPlTbAfz4GRERkY7xLywiIiIdYyEnIiLSMRZyIiIiHWMhJyIi0jEWciIiIh37f+MOrbM8aglJAAAAAElFTkSuQmCC\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 AKARI-SEP')"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfIAAAFhCAYAAACLYyg9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XlclOX6P/DPBa7ISCgECCKLsotruJZLltnJcjl1LM12tbLSsm+ePKc6dTzpLy2PxxZtM7fKNLXNrMwlj9txA1k1EBUBBUUcBRfk/v0xz9CwD8wMMwOf9+s1L2ee7b7mYeTiue9n7kuUUiAiIiLn5GLvAIiIiKj+mMiJiIicGBM5ERGRE2MiJyIicmJM5ERERE6MiZyIiMiJMZETERE5MSZyIiIiJ8ZE3oiJSKaIDLPSsd4UkWnWOFYVx94rItG2ODYRUWPHRN6AtMRaLCJ6ETkvIjtFZIqIWOXnYM3EXeG43gAmAlhssmyriFwWkYvaI82CJuYBeN3SOMm2tJ95gYi0rLC83OdORMZp2w2qw/7F2ucoV0SWioh7dcevIq6B2v+lQhE5JyL/FZGbqjm+8bGoinWnReRT07aJnAETecMbqZTSAegEYA6AlwB8bN+QavUwgB+UUsUVlk9VSrlrj3ALjv8NgCEi4mfBMRyaiDSzdwyWEJEgADcDUADurmG7hwC8C+BPSqltddh/pFLKHUB3AD0A/NXMuNoC+A7AfwC0A+AP4B8ArlR1fJPH1Cra7gngJgB/M6dtIkfBRG4nSqlCpdQ3AP4C4CERiQEAEekgImtFJE9EjonIs8Z9tKuHv4pIsnZl86mItNLWLQcQCOBb7eri/7TduotIgna18qVxe22fl0TklNZDkCYit1YT7ggA26pZVysRaS4is7X4r4mI0h7x2rm4DGA/gNur2HemiKypsOzfIrKwju8BItJRRL7Wzu1Zk6uySO1q8byIJInI3ea0rb2u7ef1kogkALgkIs20Y6Zr8SaLyGiT7XuKyEFt3Vfaz+uf5rRVxXut8j2ZxDWjus9FNSYC2A1gKYCHqmlzEoD5AIYrpXbWdX8AUErlAtgEQ0I3R5i23+dKqetKqWKl1E9KqQQz9zdt+xSAjQBi6rovkV0ppfhooAeATADDqlh+AsCTMPxhtR/AKwBaAAgBkAHDL0bj/okAOsJw9fFfAP+s7vja670AOmjbpwCYoq0LB3ASQAftdRCA0GrizgNwU4VlW7Xl+Vocg2t433Nh+CXeEUAbAL8A+BpAiMk2CwG8XcW+nQAUAWirvXYFkAOgbx3fgyuAeADvaDG0AjAQQHMAvwN4WTvnQwHotWNX27b22pyf1yHtfbfWlt2r/TxcYPgj7hIAP23/4wCe02IaA+Cq8edbW1sV3mu176m2z0UNP8PfATwFoBeAawB8KnzO1gI4DaBbPfcfpj0PAHAYwL9r+3+jrWsL4CyAz2D4g9PT3P93VbTdEUASgDfs/buCDz7q8rB7AE3pUd0vFBiS3CwAfQCcqLDurwA+Ndl/ism6OwGkV3d87fUEk9f/D8AH2vPOAM4AGAageS1xXwMQUWFZHwA6AC1huMLSo4okqm1TDKCLybInAWytsN1sAJ9U0/4OABO157cZ33Md30M/GP7waFZh+c0AcgG4mCz7HMBrNbVtcg5q+3k9WktchwDcA+AWAKcASIX3/U9z2qrje6r2c1FNjAO1z4CX9joVwPQKn7MLADaYtlnH/S9qnyEFYDOAG2r7f2OyPhKGK/0sACUwDNVU/EPhIoDzJo8nqlh3HMB70P7o4oMPZ3mwa90x+AM4B8MVYAetO/S8iJyH4arKx2TbkybPj8NwVVWTXJPnRQDcAUAp9TuAaQBeA3BGRL4QkeqOVQBDQi6jlNqjlNIrpa4opT6D4ar8zir2vQVAhlLqqMkyzwpxQTv++WraXwXgfu35A9rrur6HjgCOK6VKKizvAOCkUqrUZNlxGH4m1batqevPCyIyUUQOmWwfA8BLi+OUUkpVs685bZn7noBqPhfVeAjAT0qpfO31KlTuHp8CQzf3RyIi9dh/lDLcOzIYQAQM56QSERlvcsPaRgBQSqUopR5WSgXAcD47AFhQxfFvMHl8WMW6Tkqpp1Tle0GIHBoTuZ2J4e5afxiuvk4COFbhF45OKWWaIDuaPA8EkG3yuk7F5ZVSq5RSA2FIEgqGLvCqJEAbi6zpcAAq/gIHAG8Y/hAAAGi/5EfDcIOSqUgYur6r8hWAwSISoO1blkzr8B5OAgis4qazbAAdpfw3BwJhuDqusW2Y9/Mq+5mISCcAHwKYCqC9UuoGGIZKBIYue/8KSdD0Z21OW+a+J7OJSGsA9wEYJIY7ynMBTAfQTUS6mWx6BsCtMPQGvFeP/QEAynCD3FIYvslQiVJqpfrjhrURVaxP1fbnODc1GUzkdiIibUXkLgBfAFihlDoMw7jlBe0GqdYi4ioiMWLyVRoAT4tIgIi0g+GK7EuTdadhGDs1p/1wERkqhq8CXYah+/t6NZv/AGCQyb43iMhwEWml3cA1HoYr701V7JsIoKeIdNd+qb8JQ3L70uR4LWEYO/25qsaVUnkwjMl/CkMyS6nHe9gLQ7KcIyJttNgHANgDwzj1/4nhprzBAEbC8HOptm2TY9b28zLVRnvveVr8j+CPhLNLi32qdk7vARBXz7ZqfE91NEqLKwqGG9C6w/BH128w3MBWRimVDcN4/B0i8k5d9zexAMBtIlLrDW8iEiEiL2h/aEFEOsLQg7K7Du+RyKkxkTe8b0VED8MV1iwAbwN4BACUUtdh+IXbHcAxGG4k+wiAh8n+qwD8BMONThkA/mmy7k0Af9O6XmfUEkdLGL7+lg9DN+uNMPxhUJVlAO7UEjFguJnqn/jjZrdnYOierPRdcqXUPhjGv3/Q4vUFcKdS6prJZnfDMGaeXXF/E6tgGAs3vSI2+z2YnNvOMNxcmAXgL0qpq1r7I7TjvAfDmHhqLW2b+/My3T4Zhru6d8HwR1dXGIYkoMUxBsBjMAwxTICh1+JKXdsy8z2Z6yEYxuFPKKVyjQ8AiwCMr9jDoZQ6CUMy/7OIvFnX/bVj5MHwmfu7GfHpYbh/YI+IXIIhgScCeKHCdsZvcxgf6+pyEogcmZQfkiNHJiKZAB5XSv1ih7b/BeCMUqri2KM1jr0HwGNKqURrH9uZaeflA6XUp/aOhYgcl1NPUkENRylV3dW6NY7dx1bHdiZimAktDYar6PEAYgH8aNegiMjhMZETOY5wAKthuIM8HcCflVI59g2JiBwdu9aJiIicGG92IyIicmJM5ERERE7MKcbIvby8VFBQkL3DICJyKvv3789XSnlbsP+NzZo1+wiG+Q544WcfpQASS0pKHu/Vq9eZqjZwikQeFBSEffv22TsMIiKnIiLHLdm/WbNmH/n6+kZ6e3sXuLi48IYqOygtLZW8vLyo3Nzcj1BNCWH+hUVERNWJ8fb2vsAkbj8uLi7K29u7EDVMO8xETkRE1XFhErc/7WdQbb5mIiciInJiTORERGSxeZvSfH5JOV2u3PEvKad18zalVVVq1yH85S9/6bR///5W1a1fuHBh+8zMzOb1Ofa0adM6rF+/Xlf7lpZjIiciIot1D7yh6PnVh0KMyfyXlNO651cfCukeeEORvWOrzpdffnm8V69el6tbv2LFCq8TJ07UK5EvWLAge9SoUfr6R2c+JnIiIgezd8ManEhMKLfsRGIC9m5YY6eIajcs0kf/9n3dM55ffSjkH98mdXh+9aGQt+/rnjEs0sfiZDZs2LDQ6OjoyM6dO0fPmzfPCwDc3Nx6PPHEEwFRUVGR/fr1C8vOzq7yW1gHDhxo1bVr10jj67S0tBZhYWFRABAXFxe+fft2t5KSEowdOzaoS5cu0WFhYVH/+Mc/bvz00089ExMT3SZOnBgSERERdfHiRZkxY4ZfTExMZJcuXaLvv//+TqWlpdXGPHbs2KBPP/3U09L3bg4mciIiB+MbGobvFswpS+YnEhPw3YI58A0Ns3NkNRsW6aMf2zMg79P/ZvqN7RmQZ40kDgArV67MTEpKSjl06FDy4sWLfXJzc12Li4tdevbsWZScnJwyYMAA/cyZMztUtW/Pnj0vX7t2TZKTk1sAwLJly9qNGjWqwHSbXbt2ueXk5DQ/evRo0pEjR5Kffvrps4888khBTExM0bJlyzJSU1OT3d3d1YsvvngmMTEx5ejRo0nFxcUuX3zxRZUlixsaEzkRkYMJjInFXdNm4rsFc/Df1Svw3YI5uGvaTATGxNo7tBr9knJat/ZAlvcjA4Jy1h7I8q44Zl5fc+fO9QkPD4/q1atXZG5ubvOkpKRWLi4uePzxx88BwKOPPnp279697tXtP2rUqHMrVqxoBwDr1q3zfPDBB8+Zro+IiLhy8uTJlg899FDHNWvWtPX09Lxe1XE2btyoi42NjQgLC4vauXOnLjExsbU13p+lmMiJiBxQYEwsut1+J3av/QLdbr/TKZK4sTv91ZHR2cZudkuT+Xfffafbtm2bbt++falpaWnJkZGRxcXFxZVyl4hUe4wHH3ywYP369Z4JCQktRQRdu3a9Yrre29v7emJiYvKQIUP077333o3jxo0LqniMoqIieeGFFzp9/fXX6UeOHEmeMGFC/uXLlx0ihzpEEEREVN6JxATE//QD+o4dh/iffqg0Zu5oDp0472Y6Jm4cMz904rybJcc9f/68q4eHx3WdTld68ODBVvHx8W0AoLS0FMYx6KVLl7aPi4urths/Ojr6iouLC1555ZUOo0ePPldxfU5OTrPr16/j4YcfPv/Pf/7z1OHDh90AwN3d/XphYaErABQVFbkAgK+vb0lhYaHLt99+2yDj3+ZwiilaiYiaEuOYuLE7vWNUrMN3r88YHn664rJhkT56S8fJx44dW7hkyRLvsLCwqNDQ0MvdunW7BACtW7cuTUpKah0dHe2r0+muf/311xk1HWfMmDHn3njjjYC5c+eeqrguMzOz+WOPPRZUWloqAPD6669nAcDEiRPzn3nmmU4vvvhi6b59+1LGjx+fFxUVFR0QEHDVGEdNRKRBJtNxinrkvXv3VpxrnYiaiuT3foB7mC8Ch/UsW3bilwO4eCQXUU/dafZxRGS/Uqp3feOIj4/P7NatW35997clNze3HkVFRQftHUd1hg4d2nn69OmnR44caZUb/uLj4726desWVNU6dq0TETmYkOH90WzXFVxOPw8AuJx+Hs12XUHI8P52jozMce+99wYVFxe73H777Rcboj12rRMROZhWoTeg3QOROLcqBW36+OHSnhy0eyASrUJvsHdoDqOqq/EHH3ww8H//+1+5u9effPLJ088999xZW8VhjzYrYiInInJArUJvQJs+ftD/ehK6oR2ZxM2wfPnyE02hzYrYtU5E5IAup5/HpT050A3tiEt7csq62YkqYiInInIwl9PP49yqFLR7IBIetweVdbMzmVNVmMiJiBzMtSx9uTFx45j5tawGqcFBToZj5EREDkY3qGOlZa1Cb+A4OVWJV+RERGS5zW/4IG1j+elY0zbqsPmNJlmP3FhZrf7RmY+JnIiILBfQuwjrpoSUJfO0jTqsmxKCgN5Nsh55Q2IiJyIiy4WP0GP0BxlYNyUEG2d2wLopIRj9QQbCRzTJeuQA8Pnnn3t27do1MigoKObHH3+stjqbpZjIiYjIOsJH6NHt/jzsed8P3e7Ps0YSB5y3HnlJSYkcPnw4Ze7cuSdff/31KuOzBiZyIiIHs3fDmkrVzk4kJmDvhjV2ishMaRt1iP/cG32ezEH8596VxszryVnrkd97770FANC/f/9LWVlZLer6vs1ls0QuIh1FZIuIpIhIkog8V2H9DBFRIuJlqxiIiJyRb2gYvlswpyyZG6uh+YaG2TmyGhjHxEd/kIERc7LLutktTObOXI+8VatWCgCaNWuG69evVx+ghWx5RV4C4AWlVCSAvgCeFpEowJDkAdwGwO5T2xEROZrAmFjcNW0mvlswB/9dvcLhS5gCALL2uZUbEzeOmWftYz1yG7PZ98iVUjkAcrTnehFJAeAPIBnAOwD+D8AGW7VPROTMAmNi0e32O7F77RfoO3acYydxALj175XqkSN8hN7ScXJnrkfeUBqkHrmIBAHYDiAGwGAAtyqlnhORTAC9lVI11rtlPXIiamqM3endbr8T8T/9UK8rctYjbzxqqkdu85ndRMQdwFoA02Dobp8F4HYz9psEYBIABAYG2jJEIiKHYkzixuTdMSrWObrXyS5smshFpDkMSXylUuprEekKIBhAvHZjQgCAAyISp5TKNd1XKbUEwBLAcEVuyziJiBxJbvqRcknbOGaem36EiVzDeuR/sFkiF0Om/hhAilLqbQBQSh0GcKPJNpkwo2udiKgpibyhD5q3Ln+z942tA+F5g8PcX+WQWI/c+gYAeBDAUBE5pD3utGF7RESNQlpBJvKWJ5WVLb2cfh55y5OQVpBp38DIIdnyrvUdAGr83pxSKshW7RMROav23QKwOWkTbl0OtO3vjws7T2Fz88MY1G24vUMjB8SZ3YiIHExwcDAGPTAciTgB/a8nkYgTGPTAcAQHB9s7NHJArEdOROSA/Eo9Idf9ccD1GGKvB8G3lOPjVDVekRMRORjjmPjWVslwvzUAW1sllxszd0QLDyz02Xpya7k79Lae3KpbeGAh65HbGBM5EZGDOZ1wwjAm/sBwDB06FIMeGI7NzQ/jdILdb5CuVqx3bNGsHbNCjMl868mtulk7ZoXEeseyHrmNMZETETmYk94Xyo2JG8fMT3pfsHNk1RvccbB+9sDZGbN2zAqZs3dOh1k7ZoXMHjg7Y3DHwU22HjkAXL9+HWPGjAl69tlnO1TVjqXnBmAiJyJyOAMHDqx0Y1twcDAGDhxop4jMM7jjYP3I0JF5K1NW+o0MHZlnjSQOOG898mvXrsmoUaOCu3TpcnnhwoXZVbVjjfPDRE5ERFax9eRW3bfp33qPjxyf8236t94Vx8zry1nrkT/11FOdoqKiiufOnZtbl3bqiomciIgsZhwTnz1wdsbMuJnZxm52S5O5M9cj792798XffvutbVFRkZjbTn0wkRMRkcUS8hLcTMfEjWPmCXkJTbYe+eTJk/Nvv/32wrvuuiv02rVr1bZjKX6PnIiILPZsz2cr1SMf3HGw3tJxcmevR/7aa6+dnj59uuuYMWOCX3755dyq2rFUg9QjtxTrkRMR1R3rkTceNdUjZ9c6ERGRE2PXOhEROR3WI/8DEzkRETUKrEdOREREToeJnIiIyIkxkRMRETkxJnIiIiInxkROREQWO7NggY9+y5Zy07Hqt2zRnVmwwO71yKdNm9Zh/fr1Vpn3vSo9evSIMD6/+eabu+h0uu5DhgzpbKv2KmIiJyIii7Xu1q0o+6WZIcZkrt+yRZf90syQ1t262b0e+YIFC7JHjRpVaYa5kpISqxz/4MGDqcbnM2bMyF28ePExqxzYTEzkREQOZseOHTh2rHwuOHbsGHbs2GGniGqnGzJE32HunIzsl2aG5P7rXx2yX5oZ0mHunAzdkCF2rUcOAGPHjg0yzsvu7+/fdcaMGX69evUK/+STTzznz5/vFRMTExkeHh41fPjwUL1e7wIAJ0+ebHbbbbeFhoeHR4WHh0f9/PPPbao7vpubWw/j83vuuUfftm3b2guVWxETORGRg/H398dXX31VlsyPHTuGr776Cv7+/naOrGa6IUP0HqPuyStYttzPY9Q9edZI4oBl9cir0qpVq9L9+/enTZo0qWD8+PEFiYmJKWlpacnh4eHFCxcu9AKAKVOmBN588836tLS05KSkpOSePXtetsZ7sQVOCENE5GBcXH7C8Dui8NVXX6F3797Yt28fht8RBReXnwBMtnd41dJv2aIrXL/B23PigzmF6zd4t+nXT2+NZD537lyf77///gYAqK4e+ZgxY8wek544cWKB8fn+/ftbv/LKK/56vd710qVLroMGDSoEgJ07d+rWrFlzDACaNWuG9u3bW6V2uC3wipyIyMHo2saisPAfiInxxfbt2xET44vCwn9A1zbW3qFVyzgm3mHunAzfl1/ONnazV7wBrq6sUY+8Ip1OV9b1PWnSpOBFixadOHLkSPJLL72UfeXKFafLi04XMBFRY9fOsx88PF7FwYMp6BrbEgcPpsDD41W08+xn79CqVRwf72Y6Jm4cMy+Oj7d7PfKaFBUVuQQGBl67cuWKfPHFF+2MywcMGKB/6623vAHDTXHnzp1z2HzJrnUiIgdz7NgxbPoxGYOHtMPlyx9g8JAp2PRjMtrqohEcHGzv8Kp047RpleqR64YMsbhr3Vr1yKszc+bM7Li4uEh/f/+rkZGRRRcvXnQFgPfff//Eww8/3CksLMzLxcUFixYtOj5s2LBaa5D36tUrPCMjo1VxcbGrj49P7HvvvZc5duzYC/WJzVysR05E5GB27NgBXdvzKCz8B/z9H8CpU6vg4fEq9BduwMCBA80+DuuRNx52qUcuIh1FZIuIpIhIkog8py1/Q0QSROSQiPwkImbfaUhE1BRERbtqY+QLERoyHTExC1FY+A9ERbvaOzRyQLbsWi8B8IJS6oCI6ADsF5GfAbyllPo7AIjIswBeATDFhnEQETkV/YUExMQsLBsTb+fZDzExC6G/kODQ4+QNqaHrkefm5roOHjw4vOLyrVu3pvn6+tr1jnabJXKlVA6AHO25XkRSAPgrpZJNNmsDwPH79omIGlCnTpW/YtbOsx+TeC1sWRvc19f3empqanLtWza8BrnZTUSCAPQAsEd7PRvARACFAIY0RAxERESNkc1vpxcRdwBrAUxTSl0AAKXULKVURwArAUytZr9JIrJPRPbl5eXZOkwiIiKnZNNELiLNYUjiK5VSX1exySoAY6vaVym1RCnVWynV29vb25ZhEhEROS1b3rUuAD4GkKKUettkeReTze4GkFpxXyIiIjKPLa/IBwB4EMBQ7atmh0TkTgBzRCRRRBIA3A7gORvGQEREDWD3hnSfYwn55aZjPZaQr9u9Ib3J1CPfuXNn6+7du0d07tw5OiwsLOrDDz/0tFWbpmx51/oOAFVNfvuDrdokIiL78An2KNq8NDnk1oejMoJjvfTHEvJ1xtf2jm3BggXZVS0vKSlBs2aWp0FjPXJ3d/fS5cuXH+vateuVzMzM5jfddFPk6NGjL3h5edn062kOO3csERE5j+BYL/2tD0dlbF6aHPLb6iMdTJO6pcd2lnrksbGxV7p27XoFAIKCgq61a9euJCcnx+bfDmMiJyIiqwiO9dKH9/XNS/g1yy+8r2+eNZI44Jz1yLds2eJ27do1iYqKulLX91tXLJpCRERWcSwhX5e2O9c7dmhATtruXO+AiHZ6ayRzZ6tHfvz48eaPPPJIyMcff3zM1dX20+ryipyIiCxmOiZ+831h2cZu9oo3wNWVs9UjP3funMuIESM6v/LKK6duvfXWWqulWQMTORERWez0sUI30zFx45j56WOFTaYe+eXLl+VPf/pT53Hjxp199NFHC2rb3lrYtU5ERBbre09opXrkwbFeFnetO1M98k8++cTzf//7n3tBQUGzVatWeWnLjvXv37+4PrGZi/XIiYgaKdYjbzzsUo+ciIiIbI9d60RE5HRYj/wPTORERNQoNNV65OxaJyIicmJM5EREDuT48cU4V7Cr3LJzBbtw/PhiO0VEjo6JnIjIgejaxiIx8dmyZH6uYBcSE5+Frm2snSMjR8UxciIiB9LOsx9iYhYiMfFZ+Ps/gFOnViEmZiHaefazd2jkoHhFTkTkYNp59oO//wPIzFwEf/8HnCKJ7/himU/6/r3lpmNN379Xt+OLZXavR24qLS2tRZcuXaLtHYc1MZETETmYcwW7cOrUKgQFTcWpU6sqjZk7Ir8uEUUb350fYkzm6fv36ja+Oz/Er0tEkb1ja+yYyImIHIhxTDwmZiFCQ6aXdbM7ejIP7RWnH/H0Cxkb350fsmXpkg4b350fMuLpFzJCe9VvDnRTltYj/+2339zCw8OjunfvHvH222/faFxeUlKCyZMnB8TExESGhYVFvfXWW17GdX/72998wsLCosLDw6OeeuopfwCornb52LFjg8aPHx/Yp0+fsICAgK7ff/+9+7333hsUEhISPXbs2KCa3ts777zjFRQUFBMXFxc+bty4ThMnTgys6/lhIiciciD6CwnlxsSNY+b6Cwl2jqx2ob3i9NG33Jp3YOM3ftG33JpnjSQOWF6P/LHHHgt6++23Txw6dCjVdPmCBQu8PDw8ricmJqbEx8enfPbZZ96pqaktVq9e3fb777/33L9/f2paWlryq6++mgsA1dUuB4DCwsJmu3btOjJnzpyTf/nLX7q8+OKLp48ePZqUmpraeufOna2riiszM7P5vHnz/Pbs2ZPy22+/HTl69Gir+pwfJnIiIgfSqdPkSmPi7Tz7oVOnyXaKyHzp+/fqkrZv9u454u6cpO2bvSuOmdfX3LlzfcLDw6N69eoVWV098r1797pXte/Zs2dd9Xq965/+9KeLxm2N63755Ze2q1evbh8RERHVo0ePyIKCgmbJycmtfv7557YTJkzIN5Y79fHxuQ4Yapf36tUrPCwsLGrt2rXtk5KSyhLvn/70p/MuLi7o2bNnUfv27a/FxcUVu7q6IiwsrDg9Pb1lVbH99ttvbfr06aP38fG53rJlSzV69Oh6VUzjXetERGQx45i4sTs9sGt3vTW6103rket0utK4uLjwutQjV0rVtE7mz59/YuzYsRdMl//www9tq9pn0qRJwWvWrPm9X79+xQsXLmy/bdu2sj9UWrVqpQDA1dUVLVq0KKtG5uLigpKSkioDsFbRMl6RExGRxXKOprqZJm3jmHnO0VS71iP38vK67u7ufn3Tpk3u2rZlNcdvu+22wvfff9/7ypUrAgAJCQktL1y44HLHHXdcWL58uZdxDPz06dOuQPW1y+vr5ptvvrRnzx5dXl6e67Vr17BhwwbP+hyHV+RERGSxgeMmVqpHHtorTm/pOLk16pF//PHHmY8//nhQ69atS4cOHVp29T19+vT8zMzMll27do1USkm7du2u/fDDD+l//vOfLxyu1Ix4AAAgAElEQVQ4cMCte/fukc2bN1fDhg0rXLRo0anqapfXV3Bw8LXp06fn3HTTTZE33njjtbCwsGIPD486F2BhPXIiokaK9cgdX2FhoYuHh0fptWvXMHz48M4PP/xw/sSJE89X3I71yImIiBzQiy++2CEiIiIqLCwsOjAw8MqECRMqJfHasGudiIicTkPXI7dUbGxsxNWrV8tdPC9btuzYkiVLsiw9tlmJXERaAhgLIMh0H6XU65YGQEREZA22rEduqYSEhNTat6ofc7vWNwC4B0AJgEsmD6d0YNNxZKWV/7peVloBDmw6bqeIiIiI6sfcrvUApdQddTmwiHQEsAyAL4BSAEuUUv8WkbcAjARwFUA6gEeUUnUeE7DEjUFtsenDRAx/IgYB4Z7ISisoe01ERORMzL0i3ykiXet47BIALyilIgH0BfC0iEQB+BlAjFIqFsARAH+t43EtFhDuieFPxGDTh4nY801GuaRORETkTGq8IheRwwCUtt0jIpIB4AoAAaC0ZFwlpVQOgBztuV5EUgD4K6V+MtlsN4A/W/YW6icg3BMxt/hj3w+Z6H1nEJM4ERE5pdq61u+yRiMiEgSgB4A9FVY9CuBLa7RRV1lpBUjcfgq97wxC4vZT8A/3ZDInIqqnwk2ZPi0CdUWtI9uXTQBTnHJWd/WE3s1jeFClyWLsJS0trcVdd93V5ejRo0n2jsVaauxaV0odV0odB/BP43PTZeY0ICLuANYCmKaUumCyfBYM3e8rq9lvkojsE5F9eXl55r4fs5iOife5O6Ssm73iDXBERGSeFoG6onOrj4QUp5zVAYYkfm71kZAWgTrWIzdTSUlJvfYzd4w82vSFiLgC6FXbTiLSHIYkvlIp9bXJ8odguNofr6qZWk4ptUQp1Vsp1dvb29vMMM1zJvNCuTFx45j5mcwLtexJRERVaR3ZXt/uvrCMc6uPhJz/Nr3DudVHQtrdF5ZheoVeX425Hrmbm1uPadOmdYiNjY3YvHlzlRXcalNjIheRv4qIHkCsiFwQEb32+gwMX0mraV8B8DGAFKXU2ybL7wDwEoC7lVJ2+Uut5/BOlbrRA8I90XN4J3uEQ0TUKLSObK9v0/PGvIv/zfZr0/PGPGskcaDx1iMHgOLiYpeYmJjihISE1OHDh1+sz/mprWv9TaWUDsBbSqm2Simd9mivlKrtbvMBAB4EMFREDmmPOwEsAqAD8LO27IP6BO4ozn70ES7tLj/0f2n3Hpz96CM7RUREZB/FKWd1lw6c8XYf0CHn0oEz3sZudks11nrkgKHs6cMPP2zRuK653yN/WUTGABgIw13svyml1te0g1JqBwx3t1f0Q91CdGytYrri1PTp8H/nHbTp2weXdu8pe01E1FQYx8SN3ektO9+gt0b3emOuRw4ALVq0KG3WzLLZ0s0dI38XwBQAhwEkApgiIu9a1HIj0aZvH/i/8w5OTZ+OvIULyyV1IqK6WnT8NHYUlM97Owr0WHTcYW78rtLVE3o306RtHDO/ekLPeuQ2Zu6fAYNgmMRFAYCIfAZDUicYkrnn/eOQ/9778HrqSSZxIqq37m3dMCkpE0uigzDQU4cdBfqy146sqq+YtY5sr7d0nLwx1yO3FrPqkYvI1wCma187g4h0AjBHKXW/jeMD4Pj1yI3d6Z73j0PB51/wipyILGJM3g918MJn2fllSb2uWI+88aipHrm5V+TtAaSIyF7t9U0AdonINwCglLrb4iidlOmYeJu+feAW14fd60RkkYGeOjzUwQvvHD+N6Z186pXEqekwN5G/YtMonNjlxMPlkrZxzPxy4mEmciKqlx0FenyWnY/pnXzwWXY+Bni6M5lX0FjqkcfFxRVbemyzErlSapvWnd5FKfWLiLQG0EwpZZXvCDqz9o8/XmlZm759mMSJqF5Mx8QHeuowwNO93GuqHuuR10BEngCwBsBibVEAgBq/fkZERHV36EJRuaQ90FOHJdFBOHSBM51S1cz9+tnTMEzwcgEAlFJHAdxY4x4O6sCm45XmVM9KK8CBTcftFBER0R+mVjEmPtBTh6mdfOwUETk6cxP5FaXUVeMLEWkGw8QwTufGoLblCqQYC6jcGNTWzpERERHVnbmJfJuIvAygtYjcBuArAN/aLizbMRZI2fRhIvZ8k1FWBY0lTImIGr+0tLQWXbp0ia59y4YxaNCgzvn5+RZ9H93cRD4TQB4Mk8BMhmGa1b9Z0rA9BYR7IuYWf+z7IRMxt/gziRMRWWjz5s0+aWlp5cYE0tLSdJs3b+aYQA22bdv2u5eX13VLjmFWIldKlcJwc9tTSqk/K6U+rK78qDPISitA4vZT6H1nEBK3n2IdciIiCwUEBBStW7cuxJjM09LSdOvWrQsJCAiw+C69xlzG1N/fv2tOTo5Fk63XVsZUROQ1EckHkAogTUTyRMRpv1duHBMf/kQM+twdUtbNzmRORFR/4eHh+tGjR2esW7cuZOPGjR3WrVsXMnr06Izw8HCLv6bcmMuYWkNtV+TTYLhb/SatdGk7AH0ADBCR6bYMzFbOZF4oNyZuHDM/k3mhlj2JiKgm4eHh+m7duuXt2bPHr1u3bnnWSOJA4y5jag21Xc5PBHCbUqpsrl2lVIaITADwEwCnq9XZc3inSssCwj05Tk5EZKG0tDRdfHy8d58+fXLi4+O9Q0JC9JYm88ZextQaarsib26axI2UUnkAmtsmJCIicjbGMfHRo0dnjBgxItvYzV7xBri6YhnT2tV2RX61nuuIiKgJycrKcjMdEzeOmWdlZblZclXOMqa1q7GMqYhcB3CpqlUAWimlGuSq3FHLmJ796CO0iulabl71S7v34HLi4SrnYCciakgsY9p41FTGtMaudaWUq1KqbRUPXUMlcUfWKqYrTk2fjku79wD4o6Rpq5iudo6MiIiaCou+u9bUGUuWnpo+HZ73j0PB51+wDjkRUQNgGdM/MJFbqE3fPvC8fxzy33sfXk89ySRORGQnLGNK9XJp9x4UfP4FvJ56EgWff1HWzU5ERNQQmMgtYBwT93/nHXg/+2xZNzuTORERNRQmcgtcTjxcbkzcOGZ+OfGwnSMjIqKmgmPkFqjqK2Zt+vbhODkRETUYXpETEZHF0tPn++Tlby43i1te/mZdevp8hypj6mj1yK2BiZyIyIEsOn4aOwrKT4S2o0CPRcdP2yki87T16F6UnDwjxJjM8/I365KTZ4S09ehucRlTqpnNErmIdBSRLSKSIiJJIvKctvxe7XWpiNR7xiEiosaoe1s3TErKLEvmOwr0mJSUie5t3ewcWc28vW7VR0XNy0hOnhFy5MgbHZKTZ4RERc3L8Pa61eIKaI21HvnKlSs9IiIioiIiIqKCgoJi/P396zWbmC2vyEsAvKCUigTQF8DTIhIFIBHAGADbbdg2EZFTGuipw5LoIExKysTcjBxMSsrEkuggDPS0qPZIg/D2ulXv5zsm72TWUj8/3zF51kjiQOOtRz5+/PjC1NTU5NTU1OSoqKiiqVOn5tbn/NgskSulcpRSB7TnegApAPyVUilKqTRbtUtE5OwGeurwUAcvvHP8NB7q4OUUSRwwdKfn5H7t3THg4Zyc3K+9K46Z11djr0f+t7/9zadVq1alf/3rX/Pqc34a5K51EQkC0AMAv2BNRFSLHQV6fJadj+mdfPBZdj4GeLo7fDI3jokbu9M92/XXW6N7vbHXI9+wYYNu/fr17Xbv3l3vmd9sfrObiLgDWAtgmlLqQm3bm+w3SUT2ici+vLx6/ZFCROR0jGPiS6KD8FKIX1k3e8Ub4BzNhcJDbqZJ2zhmfqHwkEWD+425HvmRI0daPPfcc53WrFmT7u7uXn0p0lrY9IpcRJrDkMRXKqW+rsu+SqklAJYAhjKmNgiPiMjhHLpQVG5M3DhmfuhCkUNflYeGvlDptnpvr1v1lo6TN+Z65IsXL25fWFjoOmrUqM4A4OPjc3Xbtm2/1/U4NdYjt4QY+iU+A3BOKTWtivVbAcxQStVaaNxR65ETETky1iNvPGqqR27LK/IBAB4EcFhEDmnLXgbQEsB/AHgD+F5EDimlhtswDiIiokbLZolcKbUDQHUD/Ots1S4RETV+rEf+B861TkREjQLrkRMREZHTYSInIrKzD7alY2d6+XvKdqbn44Nt6XaKiJwJEzkRkZ3FBnhg6qqDZcl8Z3o+pq46iNgADztHRs6AY+RERHbWP9QLix7ogamrDmJCn0Cs2HMCix7ogf6hXrXvTE0er8iJiBxA/1AvTOgTiIW//o4JfQKdLom/mZHj81N+YbkZa37KL9S9mZHDeuQ2xkROROQAdqbnY8WeE3h2aGes2HOi0pi5o+vV1q3omZQTIcZk/lN+oe6ZlBMhvdq6sR65jTGRExHZmXFMfNEDPfD87eFl3ezOlMxv9/LQ/ycyMOOZlBMhfz+a1eGZlBMh/4kMzLjdy4P1yKvxzjvveD322GMdja/nz5/v9fjjjwfU9fwwkRMR2VlCVmG5MXHjmHlCVqGdI6ub27089Pf5euZ9mJXvd5+vZ541kjjQeOuRP/bYY+d++uknD2PRlhUrVnhNmjSpzpPX8GY3IiI7mzIotNKy/qFeTjdO/lN+oW51boH3EwFeOatzC7xv9tTprZHM586d6/P999/fAADV1SMfM2ZM56r2raoe+a+//uoBGOqRp6amun3zzTeeAKDX611rq0f+yiuv+Ov1etdLly65Dho0qOwvrarqkQMoq0fev3//SjO4tW3btnTAgAH6L7/80qNr166Xr127JvWZ6Y2JnIiILGYcEzd2p9/sqdNbo3u9sdcjnzRpUv7s2bN9w8LCLk+YMKFeYynsWiciIovtv1DkZpq0jWPm+y8UsR55DYYOHXopJyenxbp169o/9thj5+pzDF6RExGRxf4a4lepHvntXh4Wd6035nrkRqNGjSpISEhw8/b2vl6f/W1Wj9yaWI+ciKjuWI/cOQwZMqTztGnTTt9zzz3V/tFTUz1ydq0TERHZQX5+vmtQUFBMq1atSmtK4rVh1zoRETmdxlKPPDMzM9HSYzORExFRo8B65EREROR0mMiJiIicGBM5ERGRE2MiJ6L62bEAOLa9/LJj2w3LiajBMJETUf349wS+eviPZH5su+G1f097RkV2Mm9Tms8vKafL1SP/JeW0bt6mNKvXI8/MzGx+xx13hFS3Pj8/33XOnDne9T1+jx49Iuq7rz0wkRNR/QTfAty71JC8f51t+PfepYbl1OR0D7yh6PnVh0KMyfyXlNO651cfCukeeIPV65EHBQVd+/HHH6udye3s2bOuH3/88Y3Vra/NwYMHbXaHuS0wkRNR/QXfAvR+DNj+/wz/Mok3WcMiffRv39c94/nVh0L+8W1Sh+dXHwp5+77uGcMifSyaovXJJ5/0N726fv755zu8+uqrPl26dIkGgH379rXq2rVrZERERFRYWFjU4cOHW77wwgsBJ0+ebBkRERE1efLkgMLCQpd+/fqFRUVFRYaFhUWtWLHihpradHNz62FJzA2NiZys6pPET7A3Z2+5ZXtz9uKTxE/sFBHZ1LHtwL6PgVv+z/BvxTFzalKGRfrox/YMyPv0v5l+Y3sG5FmaxAFgwoQJ59auXVtWoGTDhg2effv2vWR8/Z///Mf7qaeeOp2ampqckJCQEhwcfHX+/PlZHTt2vJKampq8ePHiLDc3t9Lvv//+9+Tk5JRt27YdefnllwNKS0stDc1hMJGTVcW0j8GMbTPKkvnenL2YsW0GYtrH2DkysjrjmPi9S4Ghs/7oZmcyb7J+STmtW3sgy/uRAUE5aw9keVccM6+PAQMGFJ89e7ZZZmZm8127drX28PC4HhISctW4vl+/fpfmz5/vN2vWLN+jR4+2cHd3r1RApLS0VKZNmxYQFhYWNWTIkLAzZ860yMrKajQTojGRk1XF+cVh3qB5mLFtBhYdXIQZ22Zg3qB5iPOLs3doZG2nDpQfEzeOmZ86YM+oyE6MY+Jv39c949WR0dnGbnZrJPORI0cWrFixwnPlypXtxo4dW67U55QpU85t2LDh99atW5eOGDEi7JtvvqnU3uLFi9udPXu22eHDh1NSU1OT27dvf62qmubOymZ/kYhIRwDLAPgCKAWwRCn1bxFpB+BLAEEAMgHcp5QqsFUc1PDi/OJwX/h9WJywGJNjJzOJN1YDp1VeFnwLx8mbqEMnzruZjokbx8wPnTjvZmkX+4MPPnjuiSeeCCooKGi2bdu2tMuXL4txXXJycovIyMgr0dHRZzIyMloeOnSodVxcXNGlS5fKEnVhYaGrl5fXtZYtW6pvv/1Wl52d3cKSeByNLf8iKQHwglIqEkBfAE+LSBSAmQA2K6W6ANisvaZGZG/OXqxOW43JsZOxOm11pTFzImp8ZgwPP10xYQ+L9NHPGB5eqU55XfXu3fvypUuXXHx8fK526tTpmum65cuXtwsLC4uOiIiIOnr0aKvJkyef9fX1vd6rV6+LXbp0iZ48eXLA448/fi4+Pr5NTExM5IoVK9oFBwdfrqk9EalptcNpsHrkIrIBwCLtMVgplSMifgC2KqXCa9qX9cidh3FM3NidXvE1ETWcxlyP3FZyc3Nde/bsGZWdnX3Y3rGYsns9chEJAtADwB4APkqpHADQ/q33d/3I8SSeTSyXtI1j5olnLa7UR0RkU5mZmc379u0b+fTTT1vci9CQbH7Xnoi4A1gLYJpS6oK5XRYiMgnAJAAIDAy0XYBkVY/GPFppWZxfHK/Gicih5Obmug4ePLhSb/Du3btTfH19r9sjpvqyaSIXkeYwJPGVSqmvtcWnRcTPpGv9TFX7KqWWAFgCGLrWbRknERE1Lb6+vtdTU1OT7R2HNdisa10Ml94fA0hRSr1tsuobAA9pzx8CsMFWMRARETV2trwiHwDgQQCHReSQtuxlAHMArBaRxwCcAHCvDWMgIiJq1GyWyJVSOwBUNyB+q63aJSIiakoazcw2RETm+mBbOnaml/9W1c70fHywLd1OERHVHxM5ETU5sQEemLrqYFky35mej6mrDiI2wMPOkTmxzW/4IG1j+elR0zbqsPkNp6pHnpaW1sJYWc1ZMJETUZPTP9QLix7ogamrDuLtn9IwddVBLHqgB/qHetk7NOcV0LsI66aElCXztI06rJsSgoDeTleP3NkwkRNRk9Q/1AsT+gRi4a+/Y0KfQCZxS4WP0GP0BxlYNyUEG2d2wLopIRj9QQbCRzhdPfLr169j3LhxnTp37hw9YMCALhcvXnToOVuZyImoSdqZno8Ve07g2aGdsWLPiUpj5lQP4SP06HZ/Hva874du9+dZmsQB+9QjP3HiRKtnn332zO+//57k4eFxfdmyZZ6Wvg9bajT1WImIzGUcEzd2p/cNbc/udWtI26hD/Ofe6PNkDuI/90bIIL2lydy0HnlOTk6zquqRz5s3zy8rK6vFuHHjCrp27Xql4jGM9ch3797t7uLiAmM98sDAwJKq2vT397/Sv3//YgDo0aNHUWZmZktL3oOt8YqciJqchKzCcknbOGaekFVo58icmHFMfPQHGRgxJ7usm73iDXD10ND1yFu0aFE2m6irq6sqKSlx6K51XpETUZMzZVBopWX9Q714NW6JrH1u5cbEjWPmWfvcLL0qZz3ymjGRExGR5W79e+WKYeEjLO5aByrXI09LSytLxMuXL2/31VdftW/WrJny9va+9uabb2b7+PiU1SMfOnRo4WuvvZY7YsSIzjExMZHR0dFFtdUjdzYNVo/cEqxHTkRUd6xH3njYvR45ERER2Qa71omIqMmprh751q1b01iPnIiIyMGxHjkRERE5BCZyIiIiJ8ZETkRE5MSYyImIiJwYEzkREVls4YGFPltPbi03PerWk1t1Cw8sZD1yG2MiJyIii8V6xxbN2jErxJjMt57cqpu1Y1ZIrHcs65HbGBM5ERFZbHDHwfrZA2dnzNoxK2TO3jkdZu2YFTJ74OyMwR0HO109ciNtHveobdu2uVXVjiXvy5qYyIka0o4FwLHt5Zcd225YTuTkBnccrB8ZOjJvZcpKv5GhI/MsTeKAfeqRA0B8fHzLsWPHdv7444+PDRo0qKiqdix9b9bCRE7UkPx7Al89/EcyP7bd8Nq/pz2jIrKKrSe36r5N/9Z7fOT4nG/Tv/WuOGZeH6b1yHft2tW6qnrk8+fP95s1a5bv0aNHW7i7u1cqIGKsRx4WFhY1ZMiQMGM98uraPHfuXLNRo0Z1Xr58eYaxLrk57dgLEzlRQwq+Bbh3qSF5/zrb8O+9Sw3LiZyYcUx89sDZGTPjZmYbu9mtkcwbuh65Tqe77ufnd3Xr1q3udWnHXjhFK1FDC74F6P0YsP3/Abf8H5M4NQoJeQlupmPixjHzhLwEN0u72Bu6Hnnz5s3Vjz/+mD5kyJAu7u7upVOmTDlXVTt33323xUMH1sBETtTQjm0H9n1sSOL7PgaCb2YyJ6f3bM9nK9UjH9xxsN4a4+T2qEfetm3b0k2bNv0+ePDgMHd399KkpKRWFdux9H1ZS6OsR/5J4ieIaR+DOL+4smV7c/Yi8WwiHo151BYhEpnHOCZu7E6v+JrIiliPvPFocvXIY9rHYMa2GdibsxeAIYnP2DYDMe1j7BwZNXmnDpRP2sYx81MH7BkVETkxm3Wti8gnAO4CcEYpFaMt6wbgAwDuADIBjFdKXbB223F+cZg3aB5mbJuB+8Lvw+q01Zg3aF65K3Qiuxg4rfKy4Ft4NU7UwFiP3DxLASwCsMxk2UcAZiiltonIowBeBPB3WzQe5xeH+8Lvw+KExZgcO5lJnIiIyrAeuRmUUtsBnKuwOByAcTaMnwGMtVX7e3P2YnXaakyOnYzVaavLutmJiIgak4YeI08EcLf2/F4AHW3RiHFMfN6geZjaY2pZNzuTORERNTYNncgfBfC0iOwHoANQ7RR3IjJJRPaJyL68vLw6NZJ4NrHcmLhxzDzxbKIFoRMRETmeBv0euVIqFcDtACAiYQD+VMO2SwAsAQxfP6tLO1V9xSzOL47j5ERE1Og06BW5iNyo/esC4G8w3MFORERO7syCBT76LVvKTVuq37JFd2bBAovqkffo0SPCssiqt3LlSo+XX37ZFwA2btzoHhUVFdmsWbNen376qaet2rQFmyVyEfkcwC4A4SKSJSKPAbhfRI4ASAWQDeBTW7VPREQNp3W3bkXZL80MMSZz/ZYtuuyXZoa07tbNonrkBw8eTK24rKSkxJJDlhk/fnzhv/71r1wACAkJufrpp59mjhw58qxVDt6AbNa1rpS6v5pV/7ZVm0REZB+6IUP0HebOych+aWaIx6h78grXb/DuMHdOhm7IEIumaHVzc+tRVFR08LvvvtO98cYbfjfeeOO15ORkt/T09KRhw4aF5uTktLhy5YrLlClTTs+YMSMfANasWdP2lVde8b9+/bq0a9euZNeuXUeqOvbChQvb79u3r82yZctOhIeHXwUAFxfnmyeNc60TEZFV6IYM0XuMuievYNlyP8+JD+ZYmsQrSkhIaHPw4MGkiIiIqwCwcuXKTB8fn+sXL16UHj16RE2YMKGgtLRUpk6dGrR169bUiIiIq6dPn3a1ZgyOiImciIisQr9li65w/QZvz4kP5hSu3+Ddpl8/vTWTeWxs7CVjEgeAuXPn+nz//fc3AEBubm7zpKSkVqdPn24WFxenN27n4+PjVLO01YdTJPL9+/fni8jxBmrOCwCLBNSM58g8PE+14zmqnSXnqJM1A6mJcUzc2J3epl8/velra7Th5uZWanz+3Xff6bZt26bbt29fqk6nK42LiwsvLi52UUpBRGo6TKPjFIlcKeXdUG2JyD5LqgU1BTxH5uF5qh3PUe2c5RwVx8e7mSZt45h5cXy8m7W72AHg/Pnzrh4eHtd1Ol3pwYMHW8XHx7cBgCFDhlx64YUXOqWmprYwdq039qtyp0jkRETk2G6cNq1SPXLdkCFW7Vo3NXbs2MIlS5Z4h4WFRYWGhl7u1q3bJQDo0KFDycKFCzNHjx7dubS0FO3bt7+2c+fOo7Udb9u2bW733Xdf5wsXLrhu3rz5htmzZ3f4/fffk2wRu7U5RT3yhuQsf/3aE8+ReXieasdzVDt7niPWI3ccTa4euYWW2DsAJ8BzZB6ep9rxHNWO54hqxK71CrSpYakGPEfm4XmqHc9R7XiOLPfvf/+7/fvvv19uhrmbbrrp4vLly0/YKyZrYtc6ERFViV3rjoNd60RERI1Uk0rkInKHiKSJyO8iMrOK9YEiskVEDopIgojcabLur9p+aSIyvGEjbzj1PUciEiQixSJySHs02oI4ZpyjTiKyWTs/W0UkwGTdQyJyVHs81LCRNxwLz9F1k8/RNw0becMRkU9E5IyIVFlfWQwWaucwQUR6mqxrEp8jMk+TSeQi4grgXQAjAETBUMAlqsJmfwOwWinVA8A4AO9p+0Zpr6MB3AHgPe14jYol50iTrpTqrj2mNEjQDczMczQPwDKlVCyA1wG8qe3bDsCrAPoAiAPwqog4VZUlc1hyjjTFJp+juxskaPtYCsPvk+qMANBFe0wC8D7QdD5HZL4mk8hh+MD/rpTKUEpdBfAFgHsqbKMAtNWee8BQoQ3adl8opa4opY4B+F07XmNjyTlqKsw5R1EANmvPt5isHw7gZ6XUOaVUAYCfUfMvcmdlyTlqMpRS2wGcq2GTe2D4Y0cppXYDuEFE/NB0PkdkpqaUyP0BnDR5naUtM/UagAkikgXgBwDP1GHfxsCScwQAwVqX+zYRudmmkdqPOecoHsBY7floADoRaW/mvo2BJecIAFqJyD4R2S0io2wbqkOr7jw65Odo94Z0n2MJ+eXqkR9LyNft3pDuFPXIX3vtNZ/Q0NDosLCwqH79+oUdOXKkha3atbamlMirmny34i379wNYqpQKAHAngOUi4mLmvo2BJecoB0Cg1uX+PIBVItIWjVNWKcYAABH1SURBVI8552gGgEEichDAIACnAJSYuW9jYMk5Agyfo94AHgCwQERCbRapY6vuPDrk58gn2KNo89LkEGMyP5aQr9u8NDnEJ9jDKeqR9+rVq+jQoUMpR44cSR41alTB9OnTA2rb31E0pUSeBaCjyesAVO4WfgzAagBQSu0C0AqGggXm7NsY1PscacMOZ7Xl+wGkAwizecQNr9ZzpJTKVkqN0f6omaUtKzRn30bCknMEpVS29m8GgK0AejRAzI6ouvPokJ+j4Fgv/a0PR2VsXpoc8tvqIx02L00OufXhqIzgWC+L65EDhiIpffr0CRs5cmRweHh4NAAMGzYsNDo6OrJz587R8+bN8zLus2bNmrZRUVGR4eHhUf369av299DChQvbT5w4MRAARo4cqdfpdKUAMHDgwIs5OTm8IndA/wPQRUSCRaQFDDdqVbwj9gSAWwFARCJhSFJ52nbjRKSliATDcPPJ3gaLvOHU+xyJiLfxBkARCYHhHGU0WOQNp9ZzJCJeWi8FAPwVwCfa800AbhcRT+3mpNu1ZY1Nvc+Rdm5aGrcBMABAcoNF7li+ATBRu3u9L4BCpVQOHPhzFBzrpQ/v65uX8GuWX3hf3zxLk3hFCQkJbd56661T6enpSYChHnlSUlLKoUOHkhcvXuyTm5vrmp2d3Wzq1KlBX3/9dXpaWlry+vXr0+vazuLFi72HDRtWaM3YbanJzOymlCoRkakwfOBdAXyilEoSkdcB7FNKfQPgBQAfish0GLqqHlaGGXOSRGQ1DL9QSgA8rZRqdNV0LDlHInILgNdFpATAdQBTlFI13cjjlMw8R4MBvCkiCsB2AE9r+54TkTdgSHQA8DrPUflzBCASwGIRKYXhQmOOUqpRJnIR+RyG8+Cl3XPyKoDmAKCU+gCGe1DuhOHm2iIAj2jrHPZzdCwhX5e2O9c7dmhATtruXO+AiHZ6aybzhqhH/t5777WLj493W7x4cZq14ra1JpPIAUAp9QMM/zlMl71i8jwZhiuAqvadDWC2TQN0APU9R0qptQDW2jxAB2DGOVoDYE01+36CP67QG636niOl1E4AXW0eoANQSt1fy3qFP/7AqbjO4T5HxjFxY3d6QEQ7vbW6141sXY98/fr1unnz5vn99ttvaa1bt7b7fQfmakpd60REZCOnjxW6mSZt45j56WOFbrZor6Z65Hv27NGlpqa2AIDTp0+bNefHf//739bPPPNMpw0bNvzu7+9vnbvpGkiTuiInIiLb6HtPaKV65MGxXlbtWjdl7XrkL774YseioiLXe++9N1Q7ztVff/31d1vEbm0smkJERFVi0RTHwaIpREREjRS71omIqFFjPXKiRkBEfAEsAHATgCsAMgFMU0odsWdcRI6MXeuOg13r1KSJ4bso6wBsVUqFKqWiALwMwKI5oKtpq9FVxSMix8ZETk3BEADXtEk2AABKqUMAdojIWyKSKCKHReQvACAiX0r5WvRLRWSsiLhq2/9Pqw89WVs/WAw12lcBOKwtWy8i+0UkSUQmmRzrMRE5IoYa3B+KyCJtubeIrNWO/T8RqXI+AyKiijhGTk1BDID9VSwfA6A7gG4wzKn/PxHZDkPZzb8A+EGbYvRWAE/CMM98oVLqJm0a0f+KyE/aseIAxGhlbgHgUW0GrtbacdcCaAng7wB6AtAD+BWGKmAA8G8A7yildohIIAyzokVa7xQQUWPFRE5N2UAAn2vT7Z4WkW0wjKFvBLBQS9Z3ANiulCoWkdsBxIrIn7X9PWCYU/4qgL0mSRwAnhWR0drzjtp2vgC2GafTFJGv8EdhmWEAokxmpGorIjqllE2+g0tEjQe71qkpSALQq4rlVc7jqJS6DEPVreEwXJl/YbL9M0qp7tojWCllvCK/VHZQkcEwJOZ+SqluAA7CUFympnkjXbTtjcf2ZxInZ7Lji2U+6fv3lqtHnr5/r27HF8usei/K888/3+GVV16x+v0tzoyJnJqCXwG0FJEnjAtE5CYABQD+oo19ewO4BX9UtfsChiIVN+OPylKbADwpIs21Y4SJSJsq2vMAUKCUKhKRCAB9teV7YajB7SkizQCMNdnnJwBTTeLrbtE7Jmpgfl0iija+Oz/EmMzT9+/VbXx3fohflwiL6pFT7ZjIqdHTik+MBnCbiKSLSBKA1wCsApAAwzj1rwD+TymVq+32EwyJ/RellLHa0kcwVMA7ICKJABaj6uGpHwE0E5EEAG8A2K3FcQrAvwDsAfCLdixjqcRnAfTWbqJLBjDFSm+fqEGE9orTj3j6hYyN784P2bJ0SYeN784PGfH0CxmhveIs7ll66aWXfIOCgmL69+8fdvTo0ZYAkJSU1PLmm2/uEh0dHdmrV6/wgwcPtgKAkydPNrvttttCw8PDo8LDw6N+/vnnNkD1tcvd3Nx6PPnkk/7R0dGR/fv3D9uyZYtbXFxceEBAQNeVK1d6VBeTXq///+3df0zTdx7H8XeLDilWQOgKiugV1yLqcJjJ/FEz5MaibssYC7ss2aVmenPJzP7QZZnbkUtMZja2mOzI3Dyz+c90urlePBczCZhBjuktMN0oUASmDAdaxMG3FoVS7g/XBBl1d3xb+yPPx39A8+XzDX+8+Pz4fl/aDRs2mMxmc+7GjRtN999/f05tbW1I3iv/e3iOHLiLNBrNzLGxMfevM3K73Kr4tId7XMBkpvIc+akD++Y0njiWkb/+iZ5C219+VjuGuro63fPPP7+goaGhdWRkRJYtW5Zrs9lcVVVVSfv27bu4dOnSmzU1NYk7d+6ce/r06baNGzeaCgoK3OXl5Ve8Xq8MDAzEpaamjl6+fDnOaDSOut1uzQMPPJBbV1fXmp6ePqrRaJYfPnz4fFlZ2eAjjzyS7fF4tDU1Ne2NjY0zNm3a9IfW1tZJa3TLy8uN7e3tMw4ePHjx22+/nbFy5crFNTU1LWvXrg3JCsSdniPnsBtwd/1No9H8UW7tmZ8UkX+GeTxA0HQ0/EfvqK025K9/osdRW23IWrpMUTsjP3Xq1MwNGzb8otfrfSIixcXFv9y4cUP73XffzfQXnIiIDA8Pa0RE6uvr9Z9//vmPIiLTpk2T1NTUUZHJu8vT09OvT58+fezpp58eFBFZvHjxUHx8vC8+Pn5sxYoVQ5cuXbon0Ljq6+tnvvzyy1dERB588MEbZrM5bFsIBDlwF42Nje0I9xiAUPDvifuX07OWLlOCtbw+sV/c5/OJXq/3BpotTxSou1xEZNq0aWNa7a1dZq1WK/Hx8WMiInFxcTI6OhrwgGokrWazRw4AUK3nfKtufGj798x7zreq2jdet26d+8svv0x2u92aa9euaauqqpJ1Op0vMzNz+KOPPkoRuRXs33zzTYKIyOrVq5WKigqDiIjX65X+/n5toO5yNVatWuX+9NNPU0REGhoaZrS1tSWoveZUMSMHAKi25k9//k0fefbyFaqX1tesWeMpKSnpX7JkyeK5c+feXLFihVtE5NChQ51btmyZ/9Zbb2V4vV5NSUlJ/8qVK4f27t3bZbPZ5pvN5jStViuVlZUXA3WXq/HKK6+4ysrKFpjN5twlS5Z4LBbLUEpKyqja604Fh90AAJOiNCUwr9crw8PDGp1ON+ZwOOKLi4vNHR0dTTNmzAhJqHLYDQCAIFIURWu1Wi0jIyOasbEx2bNnz8VQhfjvIcgBAAjg6NGjs15//fXM8d+bN2/ezaqqqo6mpqaWcI1rPIIcAIAASktLB0tLS/+n0/Hhwql1AACiGEEOAEAUI8gBAIhiBDkAAFGMIAcAqDbw1QXjUMvV2/rIh1qu6ge+ukAfeYgR5AAA1e7J0nv6j7SZ/GE+1HJV33+kzXRPlj7m+8hHRkbC+vsJcgCAagmLUpXZZebO/iNtpl/+1TGn/0ibaXaZuTNhUWpM9pG/9957qevXrzetW7duodVqNau9RzUIcgBAUCQsSlUS8+91uf/9c0Zi/r2uYIR4XV2dzm63z/7hhx+ajx8/3u4vPNm8efP8999/v8vhcLRUVFR0v/jii1kiIlu3bs2yWq2K0+lsdjgczfn5+TdERD755JMLDoej5ezZs80ffvihsbe3N05EZGhoSFtYWKg4HI6WxMTE0TfeeGNuXV1d22effda+a9euuXcaW2Nj48xDhw79ePr06Ta196kGL4QBAATFUMtV/fXGK4aZq+f0XG+8YohfmKyoDfNI7SMXEbFarYNGozEsRSnjEeQAANX8e+L+5fT4hclKsJbXI7GPXEREp9P5pnI/wcbSOgBAteEuRTc+tP175sNdSkz2kUcSghwAoFrSowsuT5x5JyxKVZIeXfCbnvL/x/g+8sceeyx7fB/5xx9/nGaxWHLvu+++xUePHk0WEdm7d2/X119/rf+1Jzy3sbExobS0dMDr9WrMZnPuzp075wSjjzyS0EcOAJgUfeSR40595MzIAQCIYhx2AwAggDv1kYdrTBMR5AAABEAfOQAACCmCHACAKEaQAwAQxQhyAACiGEEOAFCturra6HQ6b+sjdzqd+urq6pjsI6+trdXZbLZ54R6HCEEOAAiCzMxMj91uN/nD3Ol06u12uykzMzMm+8jXrl3rOXDgwE/hHocIQQ4ACAKLxaKUlJR02u1204kTJ+bY7XZTSUlJp8Viick+8uPHj+sLCwsXqr23YCDIAQBBYbFYlLy8PNeZM2cy8vLyXMEI8UjuI48UvBAGABAUTqdTf+7cOUNBQUHPuXPnDCaTSVEb5pHcRx4pCHIAgGr+PXH/crrJZFKCtbweqX3kkYKldQCAat3d3brxoe3fM+/u7qaPPMQIcgCAakVFRZcnzrwtFotSVFREH3mI0UcOAJgUfeSRgz5yAABiFIfdAAAIgD5yAACiGH3kAAAgpAhyAACiGEEOAEAUI8gBAKp1dLxrdPVV31Zj6uqr1nd0vBuTNaaRhCAHAKg2K2mZp7l5h8kf5q6+an1z8w7TrKRlMVljGkkIcgCAaoa0IiU3953O5uYdpra2XXOam3eYcnPf6TSkFcVkjekzzzwzPycnJzcnJyc3JSUlb/v27Rlq73OqCHIAQFAY0oqUjPSnXD91H8jISH/KFYwQj9Qa08OHD19sbW1tPnbsWHtycrL3hRdeuKr2XqeK58gBAEHh6qvW9/R+YZiXaevp6f3CkDJ7laI2zCO5xtTj8WhKS0uz9+zZ02U2m4fV3KcaBDkAQDX/nrh/OT1l9iolWMvrkVpj+txzz81//PHHrz355JOqVx7UYGkdAKDa4MBZ3fjQ9u+ZDw6cjcka0927dxvcbnfcm2++2av2WmoR5AAA1bKzt1+eOPM2pBUp2dnbY7LGtLKyMt3pdCb4D7y9/fbbBrXXnCpqTAEAk6LGNHJQYwoAQIzisBsAAAFQYwoAQBSjxhQAAIQUQQ4AQBQjyAEAiGIEOQAAUYwgBwCotruzx3iyb+C2PvKTfQP63Z099JGHGEEOAFBt+SydZ1tLl8kf5if7BvTbWrpMy2fp6CMPMYIcAKBacVqS8vdFWZ3bWrpMfz3fPWdbS5fp74uyOovTkmKyj3z58uWW+vr6BP/X+fn5OWfOnEkI9PlQIsgBAEFRnJaklKWnuP7R3ZdRlp7iCkaIR2ofuc1m69u/f3+aiMj3338fPzw8rCkoKBhSe79TwQthAABBcbJvQH+k95phS2Zaz5HeawZril5RG+aR2kdus9muVVRUZNy8ebP7gw8+SHv22WfD9k56ghwAoJp/T9y/nG5N0SvBWl6PxD5yvV7vs1qtgwcPHkw+duzY7IaGhrC9/Y2ldQCAag2DHt340PbvmTcMemKyj1xEZOvWrX2vvvrqvLy8vOtGo3E0GNecCoIcAKDaa6aMyxNn3sVpScprpoyY7CMXEbFarZ7ExMTRTZs2hbXqlT5yAMCk6CO/swsXLkx/+OGHLR0dHU1xcXEh/V30kQMAEESVlZWpDz300KLy8vJLoQ7x38NhNwAAArhTH/lLL710NVzjGo8gBwAgAPrIAQDRzOfz+QI+goW749e/gS/QzwlyAEAgTS6XK4kwDx+fz6dxuVxJItIU6DMsrQMAJuX1ejf39vbu7+3tXSJM/MLFJyJNXq93c6AP8PgZAABRjP+wAACIYgQ5AABRjCAHACCKEeQAAEQxghwAgCj2X6foyxAjPG45AAAAAElFTkSuQmCC\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
}
